added measurements of spanning cluster sizes and conducting backbone sizes

2 files changed, 48 insertions, 28 deletions

diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index 8f7bd12..dd11342 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -54,6 +54,8 @@ class network {
     network(const network&);
 
     void set_thresholds(double, std::mt19937&);
+    virtual void break_edge(unsigned e, bool unbreak = false) {fuses[e] = !unbreak;};
+    virtual current_info get_current_info() {current_info empty; return empty;};
 };
 
 class fuse_network : public network {
@@ -63,10 +65,11 @@ class fuse_network : public network {
     fuse_network(const graph&, cholmod_common*);
 
     void fracture(hooks&, double cutoff = 1e-13);
-    current_info get_current_info();
 };
 
 class elastic_network : public network {
+  private:
+    double weight;
   public:
     problem hook_x;
     problem hook_y;
@@ -75,6 +78,7 @@ class elastic_network : public network {
     elastic_network(const elastic_network&);
 
     void fracture(hooks&, double weight = 0.5, double cutoff = 1e-10);
-    current_info get_current_info();
+    void break_edge(unsigned, bool unbreak = false) override;
+    current_info get_current_info() override;
 };
 
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index e885fcb..ed51c8e 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -290,36 +290,16 @@ elastic_network::elastic_network(const elastic_network& n) : network(n), hook_x(
 }
 
 void elastic_network::fracture(hooks& m, double weight, double cutoff) {
+  this->weight = weight;
   m.pre_fracture(*this);
 
   while (true) {
-    current_info cx = hook_x.solve(fuses);
-    current_info cy = hook_y.solve(fuses);
+    current_info ctot = this->get_current_info();
 
-    bool done_x = cx.conductivity < 1.0 / G.edges.size();
-    bool done_y = cy.conductivity < 1.0 / G.edges.size();
-
-    if (done_x && done_y) {
+    if (ctot.conductivity == 0) {
       break;
     }
 
-    current_info ctot;
-    ctot.currents.resize(G.edges.size());
-
-    if (done_x) {
-      for (unsigned i = 0; i < G.edges.size(); i++) {
-        ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity;
-      }
-    } else if (done_y) {
-      for (unsigned i = 0; i < G.edges.size(); i++) {
-        ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity;
-      }
-    } else {
-      for (unsigned i = 0; i < G.edges.size(); i++) {
-        ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity, 2) + pow(weight * cy.currents[i] / cy.conductivity, 2));
-      }
-    }
-
     unsigned max_pos = UINT_MAX;
     long double max_val = std::numeric_limits<long double>::lowest();
 
@@ -337,9 +317,7 @@ void elastic_network::fracture(hooks& m, double weight, double cutoff) {
       throw nofuseex;
     }
 
-    fuses[max_pos] = true;
-    hook_x.break_edge(max_pos);
-    hook_y.break_edge(max_pos);
+    this->break_edge(max_pos);
 
     m.bond_broken(*this, ctot, max_pos);
   }
@@ -347,3 +325,41 @@ void elastic_network::fracture(hooks& m, double weight, double cutoff) {
   m.post_fracture(*this);
 }
 
+current_info elastic_network::get_current_info() {
+  current_info cx = hook_x.solve(fuses);
+  current_info cy = hook_y.solve(fuses);
+
+  bool done_x = cx.conductivity < 1.0 / G.edges.size();
+  bool done_y = cy.conductivity < 1.0 / G.edges.size();
+
+  current_info ctot;
+  ctot.currents.resize(G.edges.size());
+
+  if (done_x && done_y) {
+    ctot.conductivity = 0;
+  } else if (done_x) {
+    ctot.conductivity = cy.conductivity;
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity;
+    }
+  } else if (done_y) {
+    ctot.conductivity = cx.conductivity;
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity;
+    }
+  } else {
+    ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); 
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity, 2) + pow(weight * cy.currents[i] / cy.conductivity, 2));
+    }
+  }
+
+  return ctot;
+}
+
+void elastic_network::break_edge(unsigned e, bool unbreak) {
+  fuses[e] = !unbreak;
+  hook_x.break_edge(e, unbreak);
+  hook_y.break_edge(e, unbreak);
+}
+