now collect correlatino functions and size distributions for truncated and untruncated measurements

2 files changed, 191 insertions, 24 deletions

diff --git a/src/measurements.cpp b/src/measurements.cpp
index 86c7dfd..324313c 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -125,6 +125,13 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncu
   sl(2 * (unsigned)ceil(Lx * Ly / 2), 0),
   sd(3 * (unsigned)ceil(Lx * Ly / 2), 0),
   sD(3 * (unsigned)ceil(Lx * Ly / 2), 0),
+  Sc(2 * (unsigned)ceil(Lx * Ly / 2), 0),
+  Ss(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),
@@ -132,7 +139,15 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncu
   Cll(Ncum),
   Cdd(Ncum),
   CDD(Ncum),
-  CsD(Ncum)
+  CsD(Ncum),
+  ccc(Ncum),
+  css(Ncum),
+  cmm(Ncum),
+  caa(Ncum),
+  cll(Ncum),
+  cdd(Ncum),
+  cDD(Ncum),
+  csD(Ncum)
 {
   N = 0;
   Nc = 0;
@@ -146,6 +161,14 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncu
     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);
+    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);
+    cll[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);
   }
 
   // FFTW setup for correlation functions
@@ -180,6 +203,14 @@ ma::~ma() {
   update_distribution_file("sd", sd, N, Lx, Ly, beta);
   update_distribution_file("sD", sD, N, Lx, Ly, beta);
 
+  update_distribution_file("Sc", Sc, Nc, Lx, Ly, beta);
+  update_distribution_file("Ss", Ss, N, 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);
@@ -189,6 +220,15 @@ ma::~ma() {
   update_field_file("CDD", CDD, N, Lx, Ly, beta, Mx, My);
   update_field_file("CsD", CsD, N, Lx, Ly, beta, Mx, My);
 
+  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("caa", caa, Na, Lx, Ly, beta, Mx, My);
+  update_field_file("cll", cll, 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);
+
   //stress_file.close();
 }
 
@@ -206,12 +246,28 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
 
     std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
+    unsigned avalanches_trunc = 0;
+
     for (auto e : avalanches.back()) {
-      fftw_forward_in[edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My)] = 1.0;
+      unsigned ind = edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My);
+      if (fftw_forward_in[ind] == 0.0) {
+        fftw_forward_in[ind] = 1.0;
+        avalanches_trunc++;
+      }
     }
+    Sa[avalanches_trunc]++;
 
     autocorrelation(Mx, My, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
+    std::fill_n(fftw_forward_in, Mx * My, 0.0);
+
+    for (auto e : avalanches.back()) {
+      unsigned ind = edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My);
+      fftw_forward_in[ind] += 1.0;
+    }
+
+    autocorrelation(Mx, My, caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+
     lv = c;
     avalanches.push_back({i});
   } else {
@@ -230,16 +286,39 @@ void ma::post_fracture(network &n) {
   // crack surface correlations
   std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
+  unsigned surface_size_trunc = 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)] = 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;
+    unsigned ind1 = edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r, Lx, Ly, Mx, My);
+    unsigned ind2 = edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My);
+    if (fftw_forward_in[ind1] == 0.0) {
+      fftw_forward_in[ind1] = 1.0;
+      surface_size_trunc++;
+    }
+    if (fftw_forward_in[ind2] == 0.0) {
+      fftw_forward_in[ind2] = 1.0;
+      surface_size_trunc++;
+    }
   }
 
   ss[crack.size()]++;
 
+  fftw_complex *Tss = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out);
+
+  correlation(Mx, My, Css, Tss, Tss, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+
+  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;
+  }
+
+  Ss[surface_size_trunc]++;
+
   fftw_complex *tss = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out);
 
-  correlation(Mx, My, Css, tss, tss, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  correlation(Mx, My, css, tss, tss, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
   unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]];
 
@@ -251,50 +330,105 @@ void ma::post_fracture(network &n) {
 
   // non-spanning clusters
   std::fill_n(fftw_forward_in, Mx * My, 0.0);
+
   for (unsigned i = 0; i < num; i++) {
     if (i != crack_component && components[i].size() > 0) {
+      unsigned cluster_size_trunc = 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;
+        unsigned ind = edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My);
+        if (fftw_forward_in[ind] == 0.0) {
+          cluster_size_trunc++;
+          fftw_forward_in[ind] = 1.0;
+        }
       }
 
+      Sc[cluster_size_trunc - 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++;
 
       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;
       }
+
+      for (auto it = components[i].begin(); it != components[i].end(); it++) {
+        unsigned ind = edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My);
+        fftw_forward_in[ind] += 1.0;
+      }
+
+      autocorrelation(Mx, My, ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+
+      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;
+      }
     }
   }
 
   // spanning cluster
-  std::fill_n(fftw_forward_in, Mx * My, 0.0);
+  // std::fill_n(fftw_forward_in, Mx * My, 0.0); we already reset in the last loop
 
   sm[components[crack_component].size() - 1]++;
+  unsigned spanning_cluster_trunc = 0;
   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;
+    unsigned ind = edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My);
+    if (fftw_forward_in[ind] == 0.0) {
+      fftw_forward_in[ind] = 1.0;
+      spanning_cluster_trunc++;
+    }
   }
 
   autocorrelation(Mx, My, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
+
+  Sm[spanning_cluster_trunc - 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;
+  }
+
+  autocorrelation(Mx, My, cmm, 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);
 
   unsigned final_broken = 0;
+  unsigned final_trunc = 0;
 
   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;
+      unsigned ind = edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My);
+      if (fftw_forward_in[ind] == 0.0) {
+        final_trunc++;
+        fftw_forward_in[ind] = 1.0;
+      }
     }
   }
 
-  fftw_complex *tDD = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out);
+  fftw_complex *TDD = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out);
 
+  SD[final_trunc]++;
   sD[final_broken]++;
 
-  correlation(Mx, My, CDD, tDD, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);
-  correlation(Mx, My, CsD, tss, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  correlation(Mx, My, CDD, TDD, TDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  correlation(Mx, My, CsD, Tss, TDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+
+  free(Tss);
+  free(TDD);
+
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
+
+  for (unsigned i = 0; i < n.G.edges.size(); i++) {
+    if (n.fuses[i]) { 
+      fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0;
+    }
+  }
+
+  fftw_complex *tDD = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out);
+
+  correlation(Mx, My, cDD, tDD, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  correlation(Mx, My, csD, tss, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
   free(tss);
   free(tDD);
@@ -303,41 +437,59 @@ void ma::post_fracture(network &n) {
 
   // rewind the last avalanche
   sl[avalanches.back().size() - 1]++;
+  unsigned last_trunc = 0;
   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;
+    unsigned ind = edge_r_to_ind(n.G.edges[e].r, Lx, Ly, Mx, My);
+    if (fftw_forward_in[ind] == 0.0) {
+      fftw_forward_in[ind] = 1.0;
+      last_trunc++;
+    }
   }
+  Sl[last_trunc - 1]++;
 
   autocorrelation(Mx, My, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
+
+  // rewind the last avalanche
+  for (auto e : avalanches.back()) {
+    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);
+
   // damage size distribution
   unsigned total_broken = 0;
 
   std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
+  unsigned damage_trunc = 0;
+
   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;
+      unsigned ind = edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My);
+      if (fftw_forward_in[ind] == 0.0) {
+        fftw_forward_in[ind] = 1.0;
+        damage_trunc++;
+      }
     }
   }
+  Sd[damage_trunc]++;
 
   autocorrelation(Mx, My, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
-  /*
-  current_info cur = n.get_current_info();
-
-  unsigned ii = avalanches.back().front();
-  long double c = logl(cur.conductivity / fabs(cur.currents[ii])) + n.thresholds[ii];
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
   for (unsigned i = 0; i < n.G.edges.size(); i++) {
-    double x;
-    if (!n.fuses[i]) {
-      stress_file << std::scientific << fabs(cur.currents[i]) << " ";
+    if (n.fuses[i]) { 
+      fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0;
     }
-  } 
-  */
+  }
+
+  autocorrelation(Mx, My, cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
   sd[total_broken]++;
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 1d263c3..c735163 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -42,6 +42,13 @@ class ma : public hooks {
     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<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
@@ -50,6 +57,14 @@ class ma : public hooks {
     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
+    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>> cll; // 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 Na;