From 2b0d533782501d7a35d7c69f62f75447d879e22e Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 4 Jun 2019 11:13:05 -0400
Subject: added measurements for the bonafide percolation process

---
 src/perc_meas.cpp | 56 +++++++++++++++++++++++++++++++++++++++++++++++--------
 src/perc_meas.hpp |  6 ++++++
 2 files changed, 54 insertions(+), 8 deletions(-)

diff --git a/src/perc_meas.cpp b/src/perc_meas.cpp
index ec7746e..9deda41 100644
--- a/src/perc_meas.cpp
+++ b/src/perc_meas.cpp
@@ -78,9 +78,13 @@ pm::pm(unsigned n, double a) :
   sm(2 * n),
   sl(2 * n),
   sn(3 * n),
+  sN(3 * n),
   ss(2 * n),
   sb(3 * n),
-  sp(3 * n)
+  sB(3 * n),
+  sp(3 * n),
+  sr(3 * n),
+  sf(3 * n)
 {
   model_string = "percolation_" + std::to_string(n) + "_" + std::to_string(a) + "_";
   for (std::vector<uint64_t> &x : sn) {
@@ -89,15 +93,28 @@ pm::pm(unsigned n, double a) :
   for (std::vector<uint64_t> &x : sb) {
     x.resize(3 * n);
   }
+  for (std::vector<uint64_t> &x : sN) {
+    x.resize(2 * n);
+  }
+  for (std::vector<uint64_t> &x : sB) {
+    x.resize(3 * n);
+  }
+  for (std::vector<uint64_t> &x : sf) {
+    x.resize(3 * n);
+  }
 }
 
 pm::~pm() {
   update_distribution_file("sm", sm, model_string);
   update_distribution_file("sl", sl, model_string);
   update_distribution_file("sn", sn, model_string);
+  update_distribution_file("sN", sN, model_string);
   update_distribution_file("ss", ss, model_string);
   update_distribution_file("sb", sb, model_string);
+  update_distribution_file("sB", sB, model_string);
   update_distribution_file("sp", sp, model_string);
+  update_distribution_file("sr", sr, model_string);
+  update_distribution_file("sf", sf, model_string);
 }
 
 void pm::pre_fracture(const network&) {
@@ -105,34 +122,53 @@ void pm::pre_fracture(const network&) {
   initialize_incremental_components(G, ds);
   incremental_components(G, ds);
   p = 0;
+  r = 0;
+  last_thres = std::numeric_limits<long double>::lowest();
 }
 
 void pm::bond_broken(const network& net, const current_info& cur, unsigned i) {
+  unsigned p_old = p;
+  for (unsigned j = 0; j < net.G.edges.size(); j++) {
+    if (!net.fuses[j] && net.thresholds[j] < net.thresholds[i] && net.thresholds[j] > last_thres) {
+      p++;
+    }
+  }
+  last_thres = net.thresholds[i];
   boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G);
   ds.union_set(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1]);
 
   boost::component_index<VertexIndex> components(parent.begin(), parent.end());
   std::vector<unsigned> counts(components.size());
 
-  for (unsigned i = 0; i < components.size(); i++) {
+  for (unsigned j = 0; j < components.size(); j++) {
     unsigned comp_size = 0;
-    BOOST_FOREACH(VertexIndex child_index, components[i]) {
+    BOOST_FOREACH(VertexIndex child_index, components[j]) {
       comp_size++;
     }
 
-    sn[p][comp_size - 1]++;
+    sn[r][comp_size - 1]++;
+    for (unsigned k = p_old; k <= p; k++) {
+      sN[k][comp_size - 1]++;
+    }
   }
 
   unsigned bb_size = 0;
 
-  for (unsigned i = 0; i <net.G.edges.size(); i++) {
-    if (cur.currents[i] > 1.0 / net.G.edges.size()) {
+  for (unsigned j = 0; j < net.G.edges.size(); j++) {
+    if (cur.currents[j] > 1.0 / net.G.edges.size()) {
       bb_size++;
     }
   }
 
-  sb[p][bb_size - 1]++;
+  sb[r][bb_size - 1]++;
+  for (unsigned k = p_old; k <= p; k++) {
+    sB[k][bb_size - 1]++;
+  }
+
+  sf[r][p]++;
+
   p++;
+  r++;
 }
 
 void pm::post_fracture(network &n) {
@@ -159,11 +195,15 @@ void pm::post_fracture(network &n) {
     } else {
       ss[components[i].size() - 1]++;
     }
-    for (unsigned j = p; j < sn.size(); j++) {
+    for (unsigned j = r; j < sn.size(); j++) {
       sn[j][components[i].size() - 1]++;
     }
+    for (unsigned j = p; j < sn.size(); j++) {
+      sN[j][components[i].size() - 1]++;
+    }
   }
 
   sp[p - 1]++;
+  sr[r - 1]++;
 }
 
diff --git a/src/perc_meas.hpp b/src/perc_meas.hpp
index b85e27e..e310189 100644
--- a/src/perc_meas.hpp
+++ b/src/perc_meas.hpp
@@ -24,6 +24,8 @@ typedef boost::graph_traits<Graph>::vertices_size_type VertexIndex;
 class pm : public hooks {
   private:
     unsigned p;
+    unsigned r;
+    long double last_thres;
     Graph G;
     std::vector<VertexIndex> rank;
     std::vector<Vertex> parent;
@@ -32,11 +34,15 @@ class pm : public hooks {
     // measurement storage
     std::vector<uint64_t> sc; // non-spanning cluster size distribution
     std::vector<std::vector<uint64_t>> sn; // non-spanning cluster size distribution
+    std::vector<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
     std::vector<uint64_t> sl; // final avalanche size distribution
     std::vector<std::vector<uint64_t>> sb; // final avalanche size distribution
+    std::vector<std::vector<uint64_t>> sB; // final avalanche size distribution
     std::vector<uint64_t> sp;
+    std::vector<uint64_t> sr;
+    std::vector<std::vector<uint64_t>> sf; // final avalanche size distribution
 
   public:
     std::string model_string;
-- 
cgit v1.2.3-54-g00ecf