removed bin counting measurement, added measurement of higher cumulants of the correlation functions

3 files changed, 63 insertions, 43 deletions

diff --git a/src/fracture.cpp b/src/fracture.cpp
index 3505166..e260e14 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, ceil(Lx), ceil(Ly), beta);
+  ma meas(Lx, Ly, 2*ceil(Lx), 2*ceil(Ly), beta, 10);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 72f074a..ffff06d 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -36,12 +36,17 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned int
   std::ifstream file(filename);
 
   uint64_t N_old = 0;
-  std::vector<uint64_t> data_old(data.size(), 0);
+  std::vector<std::vector<uint64_t>> data_old(data.size());
+  for (unsigned j = 0; j < data.size(); j++) {
+    data_old[j].resize(data[j].size(), 0);
+  }
 
   if (file.is_open()) {
     file >> N_old;
-    for (unsigned int i = 0; i < data.size(); i++) {
-      file >> data_old[i];
+    for (unsigned j = 0; j < data.size(); j++) {
+      for (unsigned int i = 0; i < data[j].size(); i++) {
+        file >> data_old[j][i];
+      }
     }
     file.close();
   }
@@ -49,8 +54,11 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned int
   std::ofstream file_out(filename);
 
   file_out <<std::fixed<< N_old + N <<  "\n";
-  for (unsigned int i = 0; i < data.size(); i++) {
-    file_out << data_old[i] + data[i] << " ";
+  for (unsigned j = 0; j < data.size(); j++) {
+    for (unsigned int i = 0; i < data[j].size(); i++) {
+      file_out << data_old[j][i] + data[j][i] << " ";
+    }
+    file_out << "\n";
   }
 
   file_out.close();
@@ -89,7 +97,7 @@ void correlation(unsigned int Mx, unsigned int My, std::vector<T>& data, const s
 }
 
 template <class T>
-void autocorrelation(unsigned int Mx, unsigned int My, std::vector<T>& out_data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
+void autocorrelation(unsigned int Mx, unsigned int My, std::vector<std::vector<T>>& out_data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
   fftw_execute(forward_plan);
 
   for (unsigned int i = 0; i < My * (Mx / 2 + 1); i++) {
@@ -99,8 +107,10 @@ void autocorrelation(unsigned int Mx, unsigned int My, std::vector<T>& out_data,
 
   fftw_execute(reverse_plan);
 
-  for (unsigned int i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) {
-    out_data[i] += (T)(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My));
+  for (unsigned j = 0; j < out_data.size(); j++) {
+    for (unsigned int i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) {
+      out_data[j][i] += (T)pow(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My), j + 1);
+    }
   }
 }
 
@@ -108,27 +118,35 @@ unsigned int edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned i
   return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly);
 }
 
-ma::ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta) :
+ma::ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta, unsigned Ncum) :
   Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(2 * (unsigned int)ceil(Lx * Ly / 2)),
   sc(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
   sa(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
   sC(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
   sA(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
   sd(3 * (unsigned int)ceil(Lx * Ly / 2), 0),
-  sb(log2(Mx < My ? Mx : My) + 1, 0),
-  Ccc((Mx / 2 + 1) * (My / 2 + 1), 0),
-  Css((Mx / 2 + 1) * (My / 2 + 1), 0),
-  Cmm((Mx / 2 + 1) * (My / 2 + 1), 0),
-  Caa((Mx / 2 + 1) * (My / 2 + 1), 0),
-  Cdd((Mx / 2 + 1) * (My / 2 + 1), 0),
-  Cll((Mx / 2 + 1) * (My / 2 + 1), 0),
-  Cee((Mx / 2 + 1) * (My / 2 + 1), 0)
+  Ccc(Ncum),
+  Css(Ncum),
+  Cmm(Ncum),
+  Caa(Ncum),
+  Cdd(Ncum),
+  Cll(Ncum),
+  Cee(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);
+  }
 
   // FFTW setup for correlation functions
   fftw_set_timelimit(1);
@@ -140,6 +158,8 @@ ma::ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta) :
 
   forward_plan = fftw_plan_dft_r2c_2d(My, Mx, fftw_forward_in, fftw_forward_out, 0);
   reverse_plan = fftw_plan_dft_c2r_2d(My, Mx, fftw_reverse_in, fftw_reverse_out, 0);
+  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_current.dat";
+  //stress_file.open(filename, std::ifstream::app);
 }
 
 ma::~ma() {
@@ -157,7 +177,6 @@ ma::~ma() {
   update_distribution_file("sA", sA, NA, Lx, Ly, beta);
   update_distribution_file("sC", sC, NC, Lx, Ly, beta);
   update_distribution_file("sd", sd, N, Lx, Ly, beta);
-  update_distribution_file("sb", sb, 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);
@@ -166,6 +185,8 @@ ma::~ma() {
   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);
+
+  //stress_file.close();
 }
 
 void ma::pre_fracture(const network&) {
@@ -243,22 +264,6 @@ void ma::post_fracture(network &n) {
     }
   }
 
-  unsigned int max_factor = log2(Mx < My ? Mx : My) + 1;
-  std::vector<std::valarray<unsigned int>> bins(max_factor);
-  for (unsigned int i = 0; i < max_factor; i++) {
-    bins[i].resize(Mx * My / (unsigned int)pow(2, 2 * i));
-  }
-
-  for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) {
-    for (unsigned int i = 0; i < max_factor; i++) {
-      bins[i][edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx / pow(2, i), My / pow(2, i))] = 1;
-    }
-  }
-
-  for (unsigned int i =0 ; i < max_factor; i++) {
-    sb[i] += bins[i].sum();
-  }
-
   // spanning cluster
   std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
@@ -311,6 +316,20 @@ void ma::post_fracture(network &n) {
 
   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];
+
+  for (unsigned int i = 0; i < n.G.edges.size(); i++) {
+    double x;
+    if (!n.fuses[i]) {
+      stress_file << std::scientific << fabs(cur.currents[i]) << " ";
+    }
+  } 
+  */
+
   sd[total_broken]++;
 
   N++;
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 3d4ff3e..aa41fd1 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -32,6 +32,7 @@ class ma : public hooks {
     unsigned int My;
     double beta;
     Graph G;
+//    std::ofstream stress_file;
     
     // measurement storage
     std::vector<uint64_t> sc; // cluster size distribution
@@ -40,13 +41,13 @@ class ma : public hooks {
     std::vector<uint64_t> sA; // avalanche size distribution
     std::vector<uint64_t> sd; // avalanche size distribution
     std::vector<uint64_t> sb; // bin size distribution
-    std::vector<uint64_t> Ccc; // cluster-cluster correlations
-    std::vector<uint64_t> Css; // surface-surface correlations
-    std::vector<uint64_t> Cmm; // surface-surface correlations
-    std::vector<uint64_t> Caa; // avalanche-avalanche correlations
-    std::vector<uint64_t> Cdd; // damage-damage distribution
-    std::vector<uint64_t> Cll; // damage-damage distribution
-    std::vector<uint64_t> Cee; // 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>> Cdd; // damage-damage distribution
+    std::vector<std::vector<uint64_t>> Cll; // damage-damage distribution
+    std::vector<std::vector<uint64_t>> Cee; // damage-damage distribution
     uint64_t Nc;
     uint64_t NC;
     uint64_t Na;
@@ -58,7 +59,7 @@ class ma : public hooks {
 
     std::list<std::list<unsigned int>> avalanches;
 
-    ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta);
+    ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta, unsigned Ncum);
     ~ma();
 
     void pre_fracture(const network &) override;