measuring the clusters before the last avalanche again

3 files changed, 82 insertions, 66 deletions

diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 3e460a9..98c6eb0 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -65,8 +65,8 @@ void network::fracture(hooks& m) {
       throw nofuseex;
     }
 
-    this->break_edge(max_pos);
     m.bond_broken(*this, c, max_pos);
+    this->break_edge(max_pos);
   }
 
   m.post_fracture(*this);
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 8b1c9d6..9150696 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -158,11 +158,11 @@ unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, u
 }
 
 ma::ma(unsigned n, double a, double beta, double weight)
-    : sn(2 * n), ss(2 * n), sm(2 * n), sl(2 * n), sb(n + 1), sd(3 * n), sa(3 * n), sA(3 * n),
-      si(10000), sI(10000), cn(pow((unsigned)sqrt(n), 2)), cs(pow((unsigned)sqrt(n), 2)),
-      cm(pow((unsigned)sqrt(n), 2)), cl(pow((unsigned)sqrt(n), 2)), cb(pow((unsigned)sqrt(n), 2)),
-      ca(pow((unsigned)sqrt(n), 2)), cA(pow((unsigned)sqrt(n), 2)), cp(pow((unsigned)sqrt(n), 2)),
-      cq(pow((unsigned)sqrt(n), 2)) {
+    : sc(2 * n), sn(2 * n), ss(2 * n), sm(2 * n), sl(2 * n), sb(n + 1), sd(3 * n), sa(3 * n),
+      sA(3 * n), si(10000), sI(10000), cc(pow((unsigned)sqrt(n), 2)), cn(pow((unsigned)sqrt(n), 2)),
+      cs(pow((unsigned)sqrt(n), 2)), cm(pow((unsigned)sqrt(n), 2)), cl(pow((unsigned)sqrt(n), 2)),
+      cb(pow((unsigned)sqrt(n), 2)), ca(pow((unsigned)sqrt(n), 2)), cA(pow((unsigned)sqrt(n), 2)),
+      cp(pow((unsigned)sqrt(n), 2)), cq(pow((unsigned)sqrt(n), 2)), last_clusters(2 * n) {
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" +
                    std::to_string(beta) + "_" + std::to_string(weight) + "_";
@@ -175,9 +175,13 @@ ma::ma(unsigned n, double a, double beta, double weight)
 }
 
 ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
-    : sn(Lx * Ly / 2), ss(Lx * Ly / 2), sm(Lx * Ly / 2), sl(Lx * Ly / 2), sb(Lx * Ly / 2 + 1),
-      sd(Lx * Ly), sa(Lx * Ly), sA(Lx * Ly), si(10000), sI(10000), cn((Lx / 2 + 1) * (Ly / 2 + 1)), cs((Lx / 2 + 1) * (Ly / 2 + 1)),
-      cm((Lx / 2 + 1) * (Ly / 2 + 1)), cl((Lx / 2 + 1) * (Ly / 2 + 1)), cb((Lx / 2 + 1) * (Ly / 2 + 1)), ca((Lx / 2 + 1) * (Ly / 2 + 1)), cA((Lx / 2 + 1) * (Ly / 2 + 1)), cp((Lx / 2 + 1) * (Ly / 2 + 1)), cq((Lx / 2 + 1) * (Ly / 2 + 1)) {
+    : sc(Lx * Ly / 2), sn(Lx * Ly / 2), ss(Lx * Ly / 2), sm(Lx * Ly / 2), sl(Lx * Ly / 2),
+      sb(Lx * Ly / 2 + 1), sd(Lx * Ly), sa(Lx * Ly), sA(Lx * Ly), si(10000), sI(10000),
+      cc((Lx / 2 + 1) * (Ly / 2 + 1)), cn((Lx / 2 + 1) * (Ly / 2 + 1)),
+      cs((Lx / 2 + 1) * (Ly / 2 + 1)), cm((Lx / 2 + 1) * (Ly / 2 + 1)),
+      cl((Lx / 2 + 1) * (Ly / 2 + 1)), cb((Lx / 2 + 1) * (Ly / 2 + 1)),
+      ca((Lx / 2 + 1) * (Ly / 2 + 1)), cA((Lx / 2 + 1) * (Ly / 2 + 1)),
+      cp((Lx / 2 + 1) * (Ly / 2 + 1)), cq((Lx / 2 + 1) * (Ly / 2 + 1)), last_clusters(Lx * Ly / 2) {
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" +
                    std::to_string(beta) + "_" + std::to_string(weight) + "_";
@@ -190,6 +194,7 @@ ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
 }
 
 ma::~ma() {
+  update_distribution_file("sc", sc, model_string);
   update_distribution_file("sn", sn, model_string);
   update_distribution_file("ss", ss, model_string);
   update_distribution_file("sm", sm, model_string);
@@ -200,6 +205,7 @@ ma::~ma() {
   update_distribution_file("sA", sA, model_string);
   update_distribution_file("si", si, model_string);
   update_distribution_file("sI", sI, model_string);
+  update_field_file("cc", cc, model_string);
   update_field_file("cl", cl, model_string);
   update_field_file("cm", cm, model_string);
   update_field_file("cs", cs, model_string);
@@ -220,6 +226,7 @@ void ma::pre_fracture(const network&) {
 void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
   long double c = net.thresholds[i] - logl(cur.currents[i]);
   if (c > lv) {
+    last_clusters = net.C;
     lv = c;
     avalanches.push_back({i});
   } else {
@@ -247,83 +254,88 @@ void ma::post_fracture(network& n) {
   for (unsigned e : crack.second) {
     cl_cs.push_back(n.G.dual_edges[e].r);
   }
-  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cl, cl_cs,
-                   crack.first);
+  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cl, cl_cs, crack.first);
 
-  std::vector<std::list<graph::coordinate>> components(n.G.dual_vertices.size());
+  std::vector<std::list<graph::coordinate>> crit_components(n.G.dual_vertices.size());
 
   for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
-    components[n.C.findroot(i)].push_back(n.G.dual_vertices[i].r);
+    crit_components[last_clusters.findroot(i)].push_back(n.G.dual_vertices[i].r);
   }
 
-  unsigned crack_component = n.C.findroot(n.G.dual_edges[avalanches.back().back()].v[0]);
-
   for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
-    if (components[i].size() > 0) {
-      if (i != crack_component) {
-        sm[components[i].size() - 1]++;
-        autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cm,
-                         components[i], crack.first);
-        autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cp,
-                         components[i], {0, 1});
-      } else {
-        ss[components[i].size() - 1]++;
-        autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cs,
-                         components[i], crack.first);
-      }
-      sn[components[i].size() - 1]++;
-      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cn,
-                       components[i], crack.first);
+    if (crit_components[i].size() > 0) {
+      sc[crit_components[i].size() - 1]++;
+      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cc, crit_components[i], crack.first);
     }
   }
 
-  std::vector<bool> vertex_in(n.G.vertices.size());
+    std::vector<std::list<graph::coordinate>> components(n.G.dual_vertices.size());
 
-  for (unsigned i = 0; i < n.G.edges.size(); i++) {
-    if (!n.backbone[i]) {
-      vertex_in[n.G.edges[i].v[0]] = true;
-      vertex_in[n.G.edges[i].v[1]] = true;
+    for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
+      components[n.C.findroot(i)].push_back(n.G.dual_vertices[i].r);
     }
-  }
 
-  unsigned bb_size = 0;
+    unsigned crack_component = n.C.findroot(n.G.dual_edges[avalanches.back().back()].v[0]);
+
+    for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
+      if (components[i].size() > 0) {
+        if (i != crack_component) {
+          sm[components[i].size() - 1]++;
+          autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cm, components[i], crack.first);
+          autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cp, components[i], {0, 1});
+        } else {
+          ss[components[i].size() - 1]++;
+          autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cs, components[i], crack.first);
+        }
+        sn[components[i].size() - 1]++;
+        autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cn, components[i], crack.first);
+      }
+    }
 
-  std::list<graph::coordinate> cb_co;
+    std::vector<bool> vertex_in(n.G.vertices.size());
 
-  for (unsigned i = 0; i < n.G.vertices.size(); i++) {
-    if (vertex_in[i]) {
-      bb_size++;
-      cb_co.push_back(n.G.vertices[i].r);
+    for (unsigned i = 0; i < n.G.edges.size(); i++) {
+      if (!n.backbone[i]) {
+        vertex_in[n.G.edges[i].v[0]] = true;
+        vertex_in[n.G.edges[i].v[1]] = true;
+      }
     }
-  }
 
-  sb[bb_size]++;
-  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cb, cb_co,
-                   crack.first);
+    unsigned bb_size = 0;
 
-  auto av_it = avalanches.rbegin();
-  av_it++;
+    std::list<graph::coordinate> cb_co;
 
-  while (av_it != avalanches.rend()) {
-    sa[(*av_it).size() - 1]++;
-    std::list<graph::coordinate> ca_co;
-    for (unsigned e : (*av_it)) {
-      ca_co.push_back(n.G.edges[e].r);
+    for (unsigned i = 0; i < n.G.vertices.size(); i++) {
+      if (vertex_in[i]) {
+        bb_size++;
+        cb_co.push_back(n.G.vertices[i].r);
+      }
     }
-    autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, ca, ca_co,
-                     crack.first);
-    autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cq, ca_co, {0, 1});
+
+    sb[bb_size]++;
+    autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cb, cb_co, crack.first);
+
+    auto av_it = avalanches.rbegin();
     av_it++;
-  }
 
-  sA[avalanches.back().size() - 1]++;
-  std::list<graph::coordinate> cA_co;
-  for (unsigned e : avalanches.back()) {
-    cA_co.push_back(n.G.edges[e].r);
-  }
-  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cA, cA_co,
-                   crack.first);
+    while (av_it != avalanches.rend()) {
+      sa[(*av_it).size() - 1]++;
+      std::list<graph::coordinate> ca_co;
+      for (unsigned e : (*av_it)) {
+        ca_co.push_back(n.G.edges[e].r);
+      }
+      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, ca, ca_co, crack.first);
+      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cq, ca_co, {0, 1});
+      av_it++;
+    }
 
-  sd[num - 1]++;
-}
+    sA[avalanches.back().size() - 1]++;
+    std::list<graph::coordinate> cA_co;
+    for (unsigned e : avalanches.back()) {
+      cA_co.push_back(n.G.edges[e].r);
+    }
+    autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cA, cA_co, crack.first);
+
+    sd[num - 1]++;
+  }
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 135f312..3b70eae 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -22,6 +22,7 @@ class ma : public hooks {
     unsigned num;
     
     // measurement storage
+    std::vector<uint64_t> sc;
     std::vector<uint64_t> sn; // non-spanning cluster size distribution
     std::vector<uint64_t> ss; // minimal spanning cluster size distribution
     std::vector<uint64_t> sm; // spanning cluster size distribution
@@ -34,6 +35,7 @@ class ma : public hooks {
     std::vector<uint64_t> si;
     std::vector<uint64_t> sI;
 
+    std::vector<uint64_t> cc;
     std::vector<uint64_t> cn;
     std::vector<uint64_t> cs;
     std::vector<uint64_t> cm;
@@ -44,6 +46,8 @@ class ma : public hooks {
     std::vector<uint64_t> cp;
     std::vector<uint64_t> cq;
 
+    ClusterTree last_clusters;
+
   public:
     long double lv;