1 files changed, 43 insertions, 67 deletions

diff --git a/src/perc_meas.cpp b/src/perc_meas.cpp
index 8270e1f..23b46b3 100644
--- a/src/perc_meas.cpp
+++ b/src/perc_meas.cpp
@@ -75,32 +75,19 @@ pm::pm(unsigned n, double a) :
   rank(2 * n),
   parent(2 * n),
   ds(&rank[0], &parent[0]),
-  sm(2 * n),
-  sl(2 * n),
   sn(3 * n),
-  sN(3 * n),
   ss(2 * n),
+  sm(2 * n),
+  sl(2 * n),
   sb(3 * n),
-  sB(3 * n),
-  sp(3 * n),
-  sr(3 * n),
-  sf(3 * n)
+  sd(3 * n)
 {
   model_string = "percolation_" + std::to_string(n) + "_" + std::to_string(a) + "_";
   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 : 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);
+    x.resize(n);
   }
 }
 
@@ -109,100 +96,73 @@ pm::pm(unsigned Lx, unsigned Ly) :
   rank(Lx * Ly / 2),
   parent(Lx * Ly / 2),
   ds(&rank[0], &parent[0]),
-  sm(Lx * Ly / 2),
-  sl(Lx * Ly / 2),
   sn(Lx * Ly),
-  sN(Lx * Ly),
   ss(Lx * Ly / 2),
+  sm(Lx * Ly / 2),
+  sl(Lx * Ly / 2),
   sb(Lx * Ly),
-  sB(Lx * Ly),
-  sp(Lx * Ly),
-  sr(Lx * Ly),
-  sf(Lx * Ly)
+  sd(Lx * Ly)
 {
   model_string = "percolation_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_";
   for (std::vector<uint64_t> &x : sn) {
     x.resize(Lx * Ly / 2);
   }
   for (std::vector<uint64_t> &x : sb) {
-    x.resize(Lx * Ly);
-  }
-  for (std::vector<uint64_t> &x : sN) {
     x.resize(Lx * Ly / 2);
   }
-  for (std::vector<uint64_t> &x : sB) {
-    x.resize(Lx * Ly);
-  }
-  for (std::vector<uint64_t> &x : sf) {
-    x.resize(Lx * Ly);
-  }
 }
 
 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("sm", sm, model_string);
+  update_distribution_file("sl", sl, 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);
+  update_distribution_file("sd", sd, model_string);
 }
 
 void pm::pre_fracture(const network&) {
   boost::remove_edge_if(trivial, G);
   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 j = 0; j < components.size(); j++) {
+  BOOST_FOREACH(VertexIndex current_index, components) {
     unsigned comp_size = 0;
-    BOOST_FOREACH(VertexIndex child_index, components[j]) {
+    BOOST_FOREACH(VertexIndex child_index, components[current_index]) {
       comp_size++;
     }
 
     sn[r][comp_size - 1]++;
-    for (unsigned k = p_old; k <= p; k++) {
-      sN[k][comp_size - 1]++;
-    }
   }
 
-  unsigned bb_size = 0;
+  std::vector<bool> vertex_in(net.G.vertices.size());
 
-  for (unsigned j = 0; j < net.G.edges.size(); j++) {
-    if (cur.currents[j] > 1.0 / net.G.edges.size()) {
-      bb_size++;
+  for (unsigned i = 0; i < net.G.edges.size(); i++) {
+    if (cur.currents[i] > 1.0 / pow(net.G.edges.size(), 2)) {
+      vertex_in[net.G.edges[i].v[0]] = true;
+      vertex_in[net.G.edges[i].v[1]] = true;
     }
   }
 
-  sb[r][bb_size - 1]++;
-  for (unsigned k = p_old; k <= p; k++) {
-    sB[k][bb_size - 1]++;
+  unsigned bb_size = 0;
+
+  for (unsigned i = 0; i < net.G.vertices.size(); i++) {
+    if (vertex_in[i]) bb_size++;
   }
 
-  sf[r][p]++;
+  sb[r][bb_size - 1]++;
 
-  p++;
   r++;
+  last_cur = cur;
 }
 
 void pm::post_fracture(network &n) {
@@ -232,12 +192,28 @@ void pm::post_fracture(network &n) {
     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]++;
+  }
+
+
+  std::vector<bool> vertex_in(n.G.vertices.size());
+
+  for (unsigned i = 0; i < n.G.edges.size(); i++) {
+    if (last_cur.currents[i] > 1.0 / n.G.edges.size()) {
+      vertex_in[n.G.edges[i].v[0]] = true;
+      vertex_in[n.G.edges[i].v[1]] = true;
     }
   }
 
-  sp[p - 1]++;
-  sr[r - 1]++;
+  unsigned bb_size = 0;
+
+  for (unsigned i = 0; i < n.G.vertices.size(); i++) {
+    if (vertex_in[i]) bb_size++;
+  }
+
+  for (unsigned i = r; i < sb.size(); i++) {
+    sb[i][bb_size - 1]++;
+  }
+
+  sd[r - 1]++;
 }