revised measurements to bin correctly, renamed some. not compatible with earlier measurements!

3 files changed, 28 insertions, 56 deletions

diff --git a/src/fracture.cpp b/src/fracture.cpp
index 6712ac3..ec00b89 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -57,7 +57,7 @@ int main(int argc, char* argv[]) {
   cholmod_common c;
   CHOL_F(start)(&c);
 
-  ma meas(Lx, Ly, 2*ceil(Lx), 2*ceil(Ly), beta, 10);
+  ma meas(Lx, Ly, 4*ceil(Lx), 4*ceil(Ly), beta, 4);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 81b13a2..86c7dfd 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -119,35 +119,31 @@ unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, u
 ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncum) :
   Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(2 * (unsigned)ceil(Lx * Ly / 2)),
   sc(2 * (unsigned)ceil(Lx * Ly / 2), 0),
-  sa(2 * (unsigned)ceil(Lx * Ly / 2), 0),
   ss(2 * (unsigned)ceil(Lx * Ly / 2), 0),
-  sC(2 * (unsigned)ceil(Lx * Ly / 2), 0),
-  sA(2 * (unsigned)ceil(Lx * Ly / 2), 0),
+  sm(2 * (unsigned)ceil(Lx * Ly / 2), 0),
+  sa(2 * (unsigned)ceil(Lx * Ly / 2), 0),
+  sl(2 * (unsigned)ceil(Lx * Ly / 2), 0),
   sd(3 * (unsigned)ceil(Lx * Ly / 2), 0),
   sD(3 * (unsigned)ceil(Lx * Ly / 2), 0),
   Ccc(Ncum),
   Css(Ncum),
   Cmm(Ncum),
   Caa(Ncum),
-  Cdd(Ncum),
   Cll(Ncum),
-  Cee(Ncum),
+  Cdd(Ncum),
   CDD(Ncum),
   CsD(Ncum)
 {
   N = 0;
   Nc = 0;
   Na = 0;
-  NC = 0;
-  NA = 0;
   for (unsigned i = 0; i < Ncum; i++) {
     Ccc[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
     Css[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
     Cmm[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
     Caa[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
-    Cdd[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
     Cll[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
-    Cee[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
+    Cdd[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
     CDD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
     CsD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
   }
@@ -176,21 +172,20 @@ ma::~ma() {
   fftw_destroy_plan(reverse_plan);
   fftw_cleanup();
 
-  update_distribution_file("sa", sa, Na, Lx, Ly, beta);
   update_distribution_file("sc", sc, Nc, Lx, Ly, beta);
   update_distribution_file("ss", ss, N, Lx, Ly, beta);
-  update_distribution_file("sA", sA, NA, Lx, Ly, beta);
-  update_distribution_file("sC", sC, NC, Lx, Ly, beta);
+  update_distribution_file("sm", sm, N, Lx, Ly, beta);
+  update_distribution_file("sa", sa, Na, Lx, Ly, beta);
+  update_distribution_file("sl", sl, N, Lx, Ly, beta);
   update_distribution_file("sd", sd, N, Lx, Ly, beta);
   update_distribution_file("sD", sD, N, Lx, Ly, beta);
 
   update_field_file("Ccc", Ccc, Nc, Lx, Ly, beta, Mx, My);
   update_field_file("Css", Css, N, Lx, Ly, beta, Mx, My);
   update_field_file("Cmm", Cmm, N, Lx, Ly, beta, Mx, My);
-  update_field_file("Cdd", Cdd, N, Lx, Ly, beta, Mx, My);
   update_field_file("Caa", Caa, Na, Lx, Ly, beta, Mx, My);
   update_field_file("Cll", Cll, N, Lx, Ly, beta, Mx, My);
-  update_field_file("Cee", Cee, N, Lx, Ly, beta, Mx, My);
+  update_field_file("Cdd", Cdd, N, Lx, Ly, beta, Mx, My);
   update_field_file("CDD", CDD, N, Lx, Ly, beta, Mx, My);
   update_field_file("CsD", CsD, N, Lx, Ly, beta, Mx, My);
 
@@ -207,14 +202,12 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
   long double c = logl(cur.conductivity / fabs(cur.currents[i])) + net.thresholds[i];
   if (c > lv && avalanches.back().size() > 0) {
     sa[avalanches.back().size() - 1]++;
-    sA[avalanches.back().size() - 1]++;
     Na++;
-    NA++;
 
     std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
     for (auto e : avalanches.back()) {
-      fftw_forward_in[edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0;
+      fftw_forward_in[edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My)] = 1.0;
     }
 
     autocorrelation(Mx, My, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
@@ -238,8 +231,8 @@ void ma::post_fracture(network &n) {
   std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
   for (auto edge : crack) {
-    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r, Lx, Ly, Mx, My)] = 0.5;
-    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My)] = 0.5;
+    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r, Lx, Ly, Mx, My)] = 1.0;
+    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My)] = 1.0;
   }
 
   ss[crack.size()]++;
@@ -261,14 +254,12 @@ void ma::post_fracture(network &n) {
   for (unsigned i = 0; i < num; i++) {
     if (i != crack_component && components[i].size() > 0) {
       for (auto it = components[i].begin(); it != components[i].end(); it++) {
-        fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] += 1.0;
+        fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] = 1.0;
       }
 
       sc[components[i].size() - 1]++;
-      sC[components[i].size() - 1]++;
       autocorrelation(Mx, My, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
       Nc++;
-      NC++;
 
       for (auto it = components[i].begin(); it != components[i].end(); it++) {
         fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] = 0.0;
@@ -279,28 +270,13 @@ void ma::post_fracture(network &n) {
   // spanning cluster
   std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
-  sC[components[crack_component].size() - 1]++;
-  NC++;
+  sm[components[crack_component].size() - 1]++;
   for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) {
-    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] += 1.0;
+    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] = 1.0;
   }
 
   autocorrelation(Mx, My, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
-  std::function<bool(unsigned)> inCrack = [&](unsigned i) -> bool {
-    return component[n.G.dual_edges[i].v[0]] == crack_component;
-  };
-
-  for (auto avalanche : avalanches) {
-    if (avalanche.end() != std::find_if(avalanche.begin(), avalanche.end(), inCrack)) {
-      for (auto edge : avalanche) {
-        fftw_forward_in[edge_r_to_ind(n.G.edges[edge].r, Lx, Ly, Mx, My)] += 1.0;
-      }
-    }
-  }
-
-  autocorrelation(Mx, My, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
-
   /// damage at end
   std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
@@ -309,7 +285,7 @@ void ma::post_fracture(network &n) {
   for (unsigned i = 0; i < n.G.edges.size(); i++) {
     if (n.fuses[i]) { 
       final_broken++;
-      fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0;
+      fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] = 1.0;
     }
   }
 
@@ -326,12 +302,11 @@ void ma::post_fracture(network &n) {
   std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
   // rewind the last avalanche
-  sA[avalanches.back().size() - 1]++;
-  NA++;
+  sl[avalanches.back().size() - 1]++;
   for (auto e : avalanches.back()) {
     boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G);
     n.break_edge(e, true);
-    fftw_forward_in[edge_r_to_ind(n.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0;
+    fftw_forward_in[edge_r_to_ind(n.G.edges[e].r, Lx, Ly, Mx, My)] = 1.0;
   }
 
   autocorrelation(Mx, My, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
@@ -344,7 +319,7 @@ void ma::post_fracture(network &n) {
   for (unsigned i = 0; i < n.G.edges.size(); i++) {
     if (n.fuses[i]) { 
       total_broken++;
-      fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0;
+      fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] = 1.0;
     }
   }
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index aca6dce..1d263c3 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -35,26 +35,23 @@ class ma : public hooks {
 //    std::ofstream stress_file;
     
     // measurement storage
-    std::vector<uint64_t> sc; // cluster size distribution
-    std::vector<uint64_t> sa; // avalanche size distribution
-    std::vector<uint64_t> ss; // avalanche size distribution
-    std::vector<uint64_t> sC; // cluster size distribution
-    std::vector<uint64_t> sA; // avalanche size distribution
-    std::vector<uint64_t> sd; // avalanche size distribution
-    std::vector<uint64_t> sD; // avalanche size distribution
+    std::vector<uint64_t> sc; // non-spanning cluster size distribution
+    std::vector<uint64_t> ss; // minimal spanning cluster size distribution
+    std::vector<uint64_t> sm; // spanning cluster size distribution
+    std::vector<uint64_t> sa; // non-final avalanche size distribution
+    std::vector<uint64_t> sl; // final avalanche size distribution
+    std::vector<uint64_t> sd; // pre-fracture damage distribution
+    std::vector<uint64_t> sD; // post-fracture damage distribution
     std::vector<std::vector<uint64_t>> Ccc; // cluster-cluster correlations
     std::vector<std::vector<uint64_t>> Css; // surface-surface correlations
     std::vector<std::vector<uint64_t>> Cmm; // surface-surface correlations
     std::vector<std::vector<uint64_t>> Caa; // avalanche-avalanche correlations
-    std::vector<std::vector<uint64_t>> Cdd; // damage-damage distribution
     std::vector<std::vector<uint64_t>> Cll; // damage-damage distribution
-    std::vector<std::vector<uint64_t>> Cee; // damage-damage distribution
+    std::vector<std::vector<uint64_t>> Cdd; // damage-damage distribution
     std::vector<std::vector<uint64_t>> CDD; // damage-damage distribution
     std::vector<std::vector<uint64_t>> CsD; // damage-damage distribution
     uint64_t Nc;
-    uint64_t NC;
     uint64_t Na;
-    uint64_t NA;
 
   public:
     long double lv;