added recording of normalizations for correlation functions

2 files changed, 57 insertions, 25 deletions

diff --git a/src/measurements.cpp b/src/measurements.cpp
index 24ffe77..2432198 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -36,13 +36,15 @@ void update_distribution_file(std::string id, const std::vector<uint64_t>& data,
 }
 
 template <class T>
-void update_field_file(std::string id, const std::vector<T>& data, std::string model_string) {
+void update_field_file(std::string id, uint64_t N, const std::vector<T>& data, std::string model_string) {
   std::string filename = model_string + id + ".dat";
   std::ifstream file(filename);
 
   std::vector<T> data_old(data.size(), 0);
+  uint64_t old_N = 0;
 
   if (file.is_open()) {
+    file >> old_N;
     for (unsigned j = 0; j < data.size(); j++) {
       file >> data_old[j];
     }
@@ -51,6 +53,7 @@ void update_field_file(std::string id, const std::vector<T>& data, std::string m
 
   std::ofstream file_out(filename);
 
+  file_out << old_N + N << "\n";
   for (unsigned j = 0; j < data.size(); j++) {
     file_out << data_old[j] + data[j] << " ";
   }
@@ -162,7 +165,7 @@ ma::ma(unsigned n, double a, double beta, double weight)
       sA(3 * n), si(10000), sI(10000), cc(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cn(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
       cs(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cm(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cl(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
       cb(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), ca(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cA(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
-      cp(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cq(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), last_clusters(2 * n) {
+      cq(pow(1+2*(unsigned)ceil(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) + "_";
@@ -172,6 +175,15 @@ ma::ma(unsigned n, double a, double beta, double weight)
   }
   Mx = 4 * ceil(sqrt(n)) * a;
   My = 4 * ceil(sqrt(n)) / a;
+  Nc = 0;
+  Nl = 0;
+  Nm = 0;
+  Ns = 0;
+  Nn = 0;
+  Nb = 0;
+  Na = 0;
+  NA = 0;
+  Nq = 0;
 }
 
 ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
@@ -181,7 +193,7 @@ ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
       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) {
+      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) + "_";
@@ -191,6 +203,15 @@ ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
   }
   Mx = Lx;
   My = Ly;
+  Nc = 0;
+  Nl = 0;
+  Nm = 0;
+  Ns = 0;
+  Nn = 0;
+  Nb = 0;
+  Na = 0;
+  NA = 0;
+  Nq = 0;
 }
 
 ma::~ma() {
@@ -205,16 +226,15 @@ 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);
-  update_field_file("cn", cn, model_string);
-  update_field_file("cb", cb, model_string);
-  update_field_file("ca", ca, model_string);
-  update_field_file("cA", cA, model_string);
-  update_field_file("cp", cp, model_string);
-  update_field_file("cq", cq, model_string);
+  update_field_file("cc", Nc, cc, model_string);
+  update_field_file("cl", Nl, cl, model_string);
+  update_field_file("cm", Nm, cm, model_string);
+  update_field_file("cs", Ns, cs, model_string);
+  update_field_file("cn", Nn, cn, model_string);
+  update_field_file("cb", Nb, cb, model_string);
+  update_field_file("ca", Na, ca, model_string);
+  update_field_file("cA", NA, cA, model_string);
+  update_field_file("cq", Nq, cq, model_string);
 }
 
 void ma::pre_fracture(const network&) {
@@ -254,6 +274,7 @@ void ma::post_fracture(network& n) {
   for (unsigned e : crack.second) {
     cl_cs.push_back(n.G.dual_edges[e].r);
   }
+  Nl += crack.second.size();
   autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cl, cl_cs, crack.first);
 
   std::vector<std::list<graph::coordinate>> crit_components(n.G.dual_vertices.size());
@@ -268,6 +289,7 @@ void ma::post_fracture(network& n) {
     }
     autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cc, crit_components[i], crack.first);
   }
+  Nc += n.G.dual_vertices.size();
 
     std::vector<std::list<graph::coordinate>> components(n.G.dual_vertices.size());
 
@@ -278,22 +300,20 @@ void ma::post_fracture(network& n) {
     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 (i != crack_component) {
-        if (components[i].size() > 0) {
+      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);
+        } else {
+          ss[components[i].size() - 1]++;
+          autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cs, components[i], crack.first);
+          Ns += components[i].size();
         }
-        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);
-      }
-
-      if (components[i].size() > 0) {
         sn[components[i].size() - 1]++;
+        autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cn, components[i], crack.first);
       }
-      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cn, components[i], crack.first);
     }
+    Nn += n.G.dual_vertices.size();
 
     std::vector<bool> vertex_in(n.G.vertices.size());
 
@@ -317,12 +337,15 @@ void ma::post_fracture(network& n) {
 
     sb[bb_size]++;
     autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cb, cb_co, crack.first);
+    Nb += bb_size;
 
     auto av_it = avalanches.rbegin();
     av_it++;
 
     while (av_it != avalanches.rend()) {
       sa[(*av_it).size() - 1]++;
+      Na += (*av_it).size();
+      Nq += (*av_it).size();
       std::list<graph::coordinate> ca_co;
       for (unsigned e : (*av_it)) {
         ca_co.push_back(n.G.edges[e].r);
@@ -338,6 +361,7 @@ void ma::post_fracture(network& n) {
       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);
+    NA += avalanches.back().size();
 
     sd[num - 1]++;
   }
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 3b70eae..ebf79e5 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -20,6 +20,15 @@ class ma : public hooks {
     unsigned Mx;
     unsigned My;
     unsigned num;
+    uint64_t Nc;
+    uint64_t Nl;
+    uint64_t Nm;
+    uint64_t Ns;
+    uint64_t Nn;
+    uint64_t Nb;
+    uint64_t Na;
+    uint64_t NA;
+    uint64_t Nq;
     
     // measurement storage
     std::vector<uint64_t> sc;
@@ -43,7 +52,6 @@ class ma : public hooks {
     std::vector<uint64_t> cb;
     std::vector<uint64_t> ca;
     std::vector<uint64_t> cA;
-    std::vector<uint64_t> cp;
     std::vector<uint64_t> cq;
 
     ClusterTree last_clusters;