mostly implemented the ability to find dead bonds using topological properties

13 files changed, 622 insertions, 123 deletions

diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index 45cc478..0452339 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -21,6 +21,8 @@ class graph {
 
     typedef struct vertex {
       coordinate r;
+      std::vector<unsigned> nd;
+      std::vector<unsigned> ne;
       std::vector<coordinate> polygon;
     } vertex;
 
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index 544f828..8211cfc 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -7,6 +7,7 @@
 #include <random>
 #include <limits>
 #include <valarray>
+#include <set>
 
 #include "problem.hpp"
 #include "graph.hpp"
@@ -16,60 +17,79 @@
 
 #define CURRENT_CUTOFF 1e-10
 
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/connected_components.hpp>
+#include <boost/graph/depth_first_search.hpp>
+#include <boost/range/combine.hpp>
+#include <boost/foreach.hpp>
+#include <boost/graph/graph_utility.hpp>
+#include <boost/graph/incremental_components.hpp>
+#include <boost/pending/disjoint_sets.hpp>
+
+struct EdgeProperties {
+  unsigned index;
+};
+
+typedef boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, boost::no_property, EdgeProperties> Graph;
+typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;
+typedef boost::graph_traits<Graph>::vertices_size_type VertexIndex;
+
 class network {
   public:
     const graph& G;
+
+    Graph bG;
+    Graph dG;
+    std::vector<VertexIndex> rank;
+    std::vector<Vertex> parent;
+    boost::disjoint_sets<VertexIndex*, Vertex*> ds;
+    std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs;
+
     std::vector<bool> fuses;
+    std::vector<bool> backbone;
     std::vector<long double> thresholds;
 
-    network(const graph&);
+    problem px;
+    problem py;
+
+    network(const graph&, cholmod_common*);
     network(const network&);
 
     void set_thresholds(double, std::mt19937&);
     void fracture(hooks&, bool one_axis = false);
+    void update_backbone(const std::vector<double>& c);
 
-    virtual void break_edge(unsigned e, bool unbreak = false) {fuses[e] = !unbreak;};
+    void break_edge(unsigned, bool unbreak = false);
     virtual current_info get_current_info() {current_info empty; return empty;};
 };
 
 class fuse_network : public network {
   private:
     double weight;
-    problem ohm_x;
-    problem ohm_y;
 
   public:
     fuse_network(const graph&, cholmod_common*, double weight = 1.0);
     fuse_network(const fuse_network&);
 
-    void break_edge(unsigned, bool unbreak = false) override;
     current_info get_current_info() override;
 };
 
 class elastic_network : public network {
   private:
     double weight;
-    problem hook_x;
-    problem hook_y;
 
   public:
     elastic_network(const graph&, cholmod_common*, double weight = 0.5);
     elastic_network(const elastic_network&);
 
-    void break_edge(unsigned, bool unbreak = false) override;
     current_info get_current_info() override;
 };
 
 class percolation_network : public network {
-  private:
-    problem px;
-    problem py;
-
   public:
     percolation_network(const graph&, cholmod_common*);
     percolation_network(const percolation_network&);
 
-    void break_edge(unsigned, bool unbreak = false) override;
     current_info get_current_info() override;
 };
 
diff --git a/lib/include/problem.hpp b/lib/include/problem.hpp
index c2c9e09..740751d 100644
--- a/lib/include/problem.hpp
+++ b/lib/include/problem.hpp
@@ -28,11 +28,13 @@ class problem {
     cholmod_common* c;
 
   public:
+    current_info sol;
+
     problem(const graph&, unsigned, cholmod_common*);
     problem(const graph&, unsigned, cholmod_sparse*, cholmod_common*);
     problem(const problem&);
     ~problem();
-    current_info solve(std::vector<bool>& fuses);
+    void solve(std::vector<bool>& fuses);
     void break_edge(unsigned, bool unbreak = false);
 };
 
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 974d3b3..bfdb952 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -33,6 +33,8 @@ graph::graph(unsigned Nx, unsigned Ny) {
   for (unsigned i = 0; i < nv; i++) {
     vertices[i].r.x = (double)((1 + i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
     vertices[i].r.y = (double)(i / (Nx / 2));
+    signed f = (1 + i / (Nx / 2)) % 2 == 1 ? 1 : -1;
+    vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)), (nv + i - Nx / 2) % nv};
     vertices[i].polygon = {
       {vertices[i].r.x - 1.0, vertices[i].r.y},
       {vertices[i].r.x, vertices[i].r.y - 1.0},
@@ -42,6 +44,7 @@ graph::graph(unsigned Nx, unsigned Ny) {
 
     dual_vertices[i].r.x = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
     dual_vertices[i].r.y = (double)(i / (Nx / 2));
+    dual_vertices[i].nd = {i, (i + Nx / 2) % nv, Nx * (i / (Nx / 2)) + (i % (Nx / 2)), (nv + i - Nx / 2) % nv};
     dual_vertices[i].polygon = {
       {dual_vertices[i].r.x - 1.0, vertices[i].r.y},
       {dual_vertices[i].r.x, vertices[i].r.y - 1.0},
@@ -67,6 +70,29 @@ graph::graph(unsigned Nx, unsigned Ny) {
     }
   }
 
+  for (vertex& v : vertices) {
+    v.ne.resize(v.nd.size());
+  }
+
+  for (unsigned i = 0; i < edges.size(); i++) {
+    for (unsigned vi : edges[i].v) {
+      auto it1 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[0]);
+      auto it2 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[1]);
+
+      unsigned d1 = std::distance(vertices[vi].nd.begin(), it1);
+      unsigned d2 = std::distance(vertices[vi].nd.begin(), it2);
+
+      unsigned dv1 = d1 < d2 ? d1 : d2;
+      unsigned dv2 = d1 < d2 ? d2 : d1;
+
+      if (dv2 - dv1 == 1) {
+        vertices[vi].ne[dv1] = i;
+      } else {
+        vertices[vi].ne[dv2] = i;
+      }
+    }
+  }
+
 }
 
 class eulerException: public std::exception
@@ -239,12 +265,15 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
 
         if (it1 == known_vertices.end()) {
           vi1 = dual_vertices.size();
-          dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{{site->p.x, site->p.y}}});
+          dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{i1},{},{{site->p.x, site->p.y}}});
           known_vertices[t1] = vi1;
         } else {
           vi1 = it1->second;
+          dual_vertices[vi1].nd.push_back(i1);
           dual_vertices[vi1].polygon.push_back({site->p.x, site->p.y});
         }
+
+        vertices[i1].nd.push_back(vi1);
         
         if (i1 < i2 || (i1 == i2 && !self_bonded)) {
           if (i1 == i2) {
@@ -359,6 +388,29 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
     }
   }
 
+  for (vertex& v : vertices) {
+    v.ne.resize(v.nd.size());
+  }
+
+  for (unsigned i = 0; i < edges.size(); i++) {
+    for (unsigned vi : edges[i].v) {
+      auto it1 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[0]);
+      auto it2 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[1]);
+
+      unsigned d1 = std::distance(vertices[vi].nd.begin(), it1);
+      unsigned d2 = std::distance(vertices[vi].nd.begin(), it2);
+
+      unsigned dv1 = d1 < d2 ? d1 : d2;
+      unsigned dv2 = d1 < d2 ? d2 : d1;
+
+      if (dv2 - dv1 == 1) {
+        vertices[vi].ne[dv1] = i;
+      } else {
+        vertices[vi].ne[dv2] = i;
+      }
+    }
+  }
+
   jcv_diagram_free(&diagram);
 }
 
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 3424b94..dcc5e0e 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -9,9 +9,20 @@ class nofuseException: public std::exception
   }
 } nofuseex;
 
-network::network(const graph& G) : G(G), fuses(G.edges.size()), thresholds(G.edges.size()) {}
+network::network(const graph& G, cholmod_common* c) :
+  G(G), bG(G.vertices.size()), dG(G.dual_vertices.size()),
+  rank(G.dual_vertices.size()), parent(G.dual_vertices.size()), ds(&rank[0], &parent[0]),
+  fuses(G.edges.size()), backbone(G.edges.size()), thresholds(G.edges.size()),
+  px(G, 0, c), py(G, 1, c) {
+  for (unsigned i = 0; i < G.edges.size(); i++) {
+    boost::add_edge(G.edges[i].v[0], G.edges[i].v[1], {i}, bG);
+  }
+  initialize_incremental_components(dG, ds);
+}
 
-network::network(const network& n) : G(n.G), fuses(n.fuses), thresholds(n.thresholds) {}
+network::network(const network& n) : G(n.G), bG(n.bG), dG(n.dG), rank(n.rank), parent(n.parent), ds(&rank[0], &parent[0]), fuses(n.fuses), backbone(n.backbone), thresholds(n.thresholds), px(n.px), py(n.py) {
+  initialize_incremental_components(dG, ds);
+}
 
 void network::set_thresholds(double beta, std::mt19937& rng) {
   if (beta == 0.0) {
@@ -36,10 +47,15 @@ void network::fracture(hooks& m, bool one_axis) {
   m.pre_fracture(*this);
 
   while (true) {
-    current_info c = this->get_current_info();
-
     double min_cond = 1.0 / G.edges.size();
 
+    current_info c = this->get_current_info();
+    if (px.sol.conductivity[0] > min_cond) {
+      this->update_backbone(px.sol.currents);
+    } else {
+      this->update_backbone(py.sol.currents);
+    }
+
     if (c.conductivity[0] < min_cond && c.conductivity[1] < min_cond) {
       break;
     }
@@ -52,7 +68,7 @@ void network::fracture(hooks& m, bool one_axis) {
     long double max_val = std::numeric_limits<long double>::lowest();
 
     for (unsigned i = 0; i < G.edges.size(); i++) {
-      if (!fuses[i] && c.currents[i] > CURRENT_CUTOFF) {
+      if (!backbone[i]) {
         long double val = logl(c.currents[i]) - thresholds[i];
 
         if (val > max_val) {
@@ -73,43 +89,413 @@ void network::fracture(hooks& m, bool one_axis) {
   m.post_fracture(*this);
 }
 
+void network::update_backbone(const std::vector<double>& c) {
+  class get_cycle_edges : public boost::default_dfs_visitor {
+    public:
+      std::set<unsigned> tE;
+      std::list<unsigned>& bE;
+
+      get_cycle_edges(std::list<unsigned>& bE) : bE(bE) {};
+
+      void tree_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
+        tE.insert(g[e].index);
+      }
+
+      void back_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
+        if (tE.find(g[e].index) == tE.end()) {
+          bE.push_back(g[e].index);
+        }
+      }
+  };
+
+  class find_cycle : public boost::default_dfs_visitor {
+    public:
+      const graph& G;
+      std::array<unsigned, 2>& E;
+      unsigned end;
+      struct done{};
+
+      find_cycle(const graph& G, std::array<unsigned, 2>& E, unsigned end) : G(G), E(E), end(end) {}
+
+      void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) {
+        if (v == end) {
+          throw done{};
+        }
+      }
+
+      void examine_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
+        unsigned e = g[ed].index;
+        if (G.dual_edges[e].crossings[0]) E[0]++;
+        if (G.dual_edges[e].crossings[1]) E[1]++;
+      }
+
+      void finish_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
+        unsigned e = g[ed].index;
+        if (G.dual_edges[e].crossings[0]) E[0]--;
+        if (G.dual_edges[e].crossings[1]) E[1]--;
+      }
+  };
+
+  class get_cluster_edges : public boost::default_dfs_visitor {
+    public:
+      std::set<unsigned>& E;
+
+      get_cluster_edges(std::set<unsigned>& E) : E(E) {}
+
+      void examine_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
+        unsigned e = g[ed].index;
+        E.insert(e);
+      }
+  };
+
+  for (unsigned i = 0; i < G.edges.size(); i++) {
+    if ((!backbone[i]) && boost::same_component(G.dual_edges[i].v[0], G.dual_edges[i].v[1], ds)) {
+      auto it_search = seen_pairs.find({G.dual_edges[i].v[0], G.dual_edges[i].v[1]});
+      if (it_search == seen_pairs.end()) {
+      std::list<unsigned> cedges = {};
+      get_cycle_edges vis(cedges);
+      std::vector<boost::default_color_type> cm1(G.dual_vertices.size());
+      boost::depth_first_visit(dG, G.dual_edges[i].v[0], vis, boost::make_iterator_property_map(cm1.begin(), boost::get(boost::vertex_index, dG), cm1[0]));
+
+      for (unsigned e : cedges) {
+        boost::remove_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], dG);
+      }
+
+      std::array<unsigned, 2> cycle = {0, 0};
+      if (G.dual_edges[i].crossings[0]) cycle[0]++;
+      if (G.dual_edges[i].crossings[1]) cycle[1]++;
+      find_cycle vis2(G, cycle, G.dual_edges[i].v[1]);
+      std::vector<boost::default_color_type> cm2(G.dual_vertices.size());
+      try {
+        boost::depth_first_visit(dG, G.dual_edges[i].v[0], vis2, boost::make_iterator_property_map(cm2.begin(), boost::get(boost::vertex_index, dG), cm2[0]));
+      } catch(find_cycle::done const&) {
+        std::list<std::array<unsigned, 2>> other_cycles = {cycle};
+        for (unsigned e : cedges) {
+          std::array<unsigned, 2> tcycle = {0, 0};
+          if (G.dual_edges[i].crossings[0]) tcycle[0]++;
+          if (G.dual_edges[i].crossings[1]) tcycle[1]++;
+          if (G.dual_edges[e].crossings[0]) tcycle[0]++;
+          if (G.dual_edges[e].crossings[1]) tcycle[1]++;
+          find_cycle vis3(G, tcycle, G.dual_edges[e].v[0]);
+          std::vector<boost::default_color_type> cm3(G.dual_vertices.size());
+          try {
+            boost::depth_first_visit(dG, G.dual_edges[i].v[0], vis3, boost::make_iterator_property_map(cm3.begin(), boost::get(boost::vertex_index, dG), cm3[0]));
+          } catch(find_cycle::done const&) {
+          }
+          find_cycle vis4(G, tcycle, G.dual_edges[e].v[1]);
+          std::vector<boost::default_color_type> cm4(G.dual_vertices.size());
+          try {
+            boost::depth_first_visit(dG, G.dual_edges[i].v[1], vis4, boost::make_iterator_property_map(cm4.begin(), boost::get(boost::vertex_index, dG), cm4[0]));
+          } catch(find_cycle::done const&) {
+          }
+          other_cycles.push_back(tcycle);
+        }
+
+        if (std::any_of(other_cycles.begin(), other_cycles.end(), [](std::array<unsigned, 2> c){return c[0] % 2 == 0 && c[1] % 2 == 0;})) {
+          seen_pairs[{G.edges[i].v[0], G.edges[i].v[1]}] = true;
+          backbone[i] = true;
+          boost::remove_edge(G.edges[i].v[0], G.edges[i].v[1], bG);
+
+          for (unsigned j = 0; j < 2; j++) {
+            std::set<unsigned> component_edges = {};
+            get_cluster_edges visc(component_edges);
+            std::vector<boost::default_color_type> cm2(G.vertices.size());
+            boost::depth_first_visit(bG, G.edges[i].v[j], visc, boost::make_iterator_property_map(cm2.begin(), boost::get(boost::vertex_index, bG), cm2[0]));
+
+            std::array<double, 2> total_cur = {0.0, 0.0};
+            for (unsigned e : component_edges) {
+              if (G.edges[e].crossings[0]) {
+                if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) {
+                  total_cur[0] += c[e];
+                } else {
+                  total_cur[0] -= c[e];
+                }
+              }
+              if (G.edges[e].crossings[1]) {
+                if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) {
+                  total_cur[1] += c[e];
+                } else {
+                  total_cur[1] -= c[e];
+                }
+              }
+            }
+
+
+            if (fabs(total_cur[0]) < 1.0 / G.edges.size() && fabs(total_cur[1]) < 1.0 / G.edges.size()) {
+              for (unsigned e : component_edges) {
+                backbone[e] = true;
+                boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
+              }
+            }
+          }
+        }
+      }
+      for (unsigned e : cedges) {
+        boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG);
+      }
+    }
+    }
+  }
+
+  class get_tie_edges : public boost::default_dfs_visitor {
+    public:
+      std::set<unsigned>& LE;
+      std::set<unsigned>& RE;
+      unsigned start_edge;
+      bool initialized;
+      bool done_left;
+      unsigned disc;
+      unsigned dr;
+
+      get_tie_edges(std::set<unsigned>& LE, std::set<unsigned>& RE) : LE(LE), RE(RE) {
+        initialized = false;
+        done_left = false;
+        disc = 0;
+        dr = 0;
+      }
 
-fuse_network::fuse_network(const graph& G, cholmod_common* c, double weight) :
-  network(G), ohm_x(G, 0, c), ohm_y(G, 1, c) {}
+      void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor vd, const Graph& g) {
+        if (!initialized && disc > 0) {
+          start_edge = vd;
+          initialized = true;
+        }
+        disc++;
+      }
 
-fuse_network::fuse_network(const fuse_network& n) :
-  network(n), ohm_x(n.ohm_x), ohm_y(n.ohm_y) {
+      void finish_vertex(boost::graph_traits<Graph>::vertex_descriptor vd, const Graph& g) {
+        if (vd == start_edge) {
+          done_left = true;
+        }
+      }
+
+      void examine_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
+        unsigned e = g[ed].index;
+        if (!done_left) {
+          LE.insert(e);
+        } else if (LE.find(e) == LE.end()) {
+          RE.insert(e);
+        }
+      }
+
+
+  };
+
+  std::set<unsigned> bb_verts;
+
+  for (unsigned i = 0; i < G.edges.size(); i++) {
+    if (!backbone[i]) {
+      bb_verts.insert(G.edges[i].v[0]);
+      bb_verts.insert(G.edges[i].v[1]);
+    }
+  }
+
+  for (unsigned v : bb_verts) {
+    bool found_pair = false;
+    unsigned d1, d2, p1, p2;
+    for (unsigned i = 0; i < G.vertices[v].nd.size(); i++) {
+      for (unsigned j = 2; j < G.vertices[v].nd.size() - 1; j++) {
+        unsigned l = (i+j)%G.vertices[v].nd.size();
+        if (boost::same_component(G.vertices[v].nd[i], G.vertices[v].nd[l], ds)) {
+          unsigned il = i < l ? i : l;
+          unsigned ig = i < l ? l : i;
+          auto it_search = seen_pairs.find({G.vertices[v].nd[il], G.vertices[v].nd[ig]});
+          if (it_search == seen_pairs.end()) {
+          bool any_intervening1 = false;
+          bool any_intervening_e_1 = false;
+          bool any_intervening2 = false;
+          bool any_intervening_e_2 = false;
+
+          for (unsigned k = il + 1; k < ig; k++) {
+            if (!boost::same_component(G.vertices[v].nd[i], G.vertices[v].nd[k], ds)) {
+              any_intervening1 = true;
+              break;
+            }
+          }
+          for (unsigned k = (ig + 1); k%G.vertices[v].nd.size() != il; k++) {
+            if (!boost::same_component(G.vertices[v].nd[i], G.vertices[v].nd[k%G.vertices[v].nd.size()], ds)) {
+              any_intervening2 = true;
+              break;
+            }
+          }
+          if (any_intervening2 && !any_intervening1) {
+            for (unsigned k = il + 1; k <= ig; k++) {
+              if (!backbone[G.vertices[v].ne[k]]) {
+                any_intervening_e_1 = true;
+                break;
+              }
+            }
+          }
+          if (any_intervening1 && !any_intervening2) {
+            for (unsigned k = ig + 1; k%G.vertices[v].nd.size() != il + 1; k++) {
+              if (!backbone[G.vertices[v].ne[k%G.vertices[v].nd.size()]]) {
+                any_intervening_e_2 = true;
+                break;
+              }
+            }
+          }
+          if ((any_intervening1 && any_intervening2) || (any_intervening1 && any_intervening_e_2) || (any_intervening2 && any_intervening_e_1)) {
+            found_pair = true;
+            p1 = il;
+            p2 = ig;
+            d1 = G.vertices[v].nd[il];
+            d2 = G.vertices[v].nd[ig];
+            break;
+          }
+        }
+        }
+      }
+    }
+
+    if (found_pair) {
+      std::list<unsigned> cedges = {};
+      get_cycle_edges vis(cedges);
+      std::vector<boost::default_color_type> cm1(G.dual_vertices.size());
+      boost::depth_first_visit(dG, d1, vis, boost::make_iterator_property_map(cm1.begin(), boost::get(boost::vertex_index, dG), cm1[0]));
+
+      for (unsigned e : cedges) {
+        boost::remove_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], dG);
+      }
+
+      std::array<unsigned, 2> cycle = {0, 0};
+      for (unsigned k = p1; k < p2; k++) {
+        if (G.dual_edges[G.vertices[v].ne[k]].crossings[0]) cycle[0]++;
+        if (G.dual_edges[G.vertices[v].ne[k]].crossings[1]) cycle[1]++;
+      }
+      find_cycle vis2(G, cycle, d2);
+      std::vector<boost::default_color_type> cm2(G.dual_vertices.size());
+      try {
+        boost::depth_first_visit(dG, d1, vis2, boost::make_iterator_property_map(cm2.begin(), boost::get(boost::vertex_index, dG), cm2[0]));
+      } catch(find_cycle::done const&) {
+        std::list<std::array<unsigned, 2>> other_cycles = {cycle};
+        for (unsigned e : cedges) {
+          std::array<unsigned, 2> tcycle = {0, 0};
+          for (unsigned k = p1; k < p2; k++) {
+            if (G.dual_edges[G.vertices[v].ne[k]].crossings[0]) tcycle[0]++;
+            if (G.dual_edges[G.vertices[v].ne[k]].crossings[1]) tcycle[1]++;
+          }
+          if (G.dual_edges[e].crossings[0]) tcycle[0]++;
+          if (G.dual_edges[e].crossings[1]) tcycle[1]++;
+          find_cycle vis3(G, tcycle, G.dual_edges[e].v[0]);
+          std::vector<boost::default_color_type> cm3(G.dual_vertices.size());
+          try {
+            boost::depth_first_visit(dG, d1, vis3, boost::make_iterator_property_map(cm3.begin(), boost::get(boost::vertex_index, dG), cm3[0]));
+          } catch(find_cycle::done const&) {
+          }
+          find_cycle vis4(G, tcycle, G.dual_edges[e].v[1]);
+          std::vector<boost::default_color_type> cm4(G.dual_vertices.size());
+          try {
+            boost::depth_first_visit(dG, d2, vis4, boost::make_iterator_property_map(cm4.begin(), boost::get(boost::vertex_index, dG), cm4[0]));
+          } catch(find_cycle::done const&) {
+          }
+          other_cycles.push_back(tcycle);
+        }
+
+        if (std::any_of(other_cycles.begin(), other_cycles.end(), [](std::array<unsigned, 2> c){return c[0] % 2 == 0 && c[1] % 2 == 0;})) {
+          seen_pairs[{d1, d2}] = true;
+
+      std::set<unsigned> left_edges = {};
+      std::set<unsigned> right_edges = {};
+      get_tie_edges visc(left_edges, right_edges);
+      std::vector<boost::default_color_type> cm2(G.vertices.size());
+      boost::depth_first_visit(bG, v, visc, boost::make_iterator_property_map(cm2.begin(), boost::get(boost::vertex_index, bG), cm2[0]));
+
+            std::array<double, 2> total_cur_left = {0.0, 0.0};
+            for (unsigned e : left_edges) {
+              if (G.edges[e].crossings[0]) {
+                if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) {
+                  total_cur_left[0] += c[e];
+                } else {
+                  total_cur_left[0] -= c[e];
+                }
+              }
+              if (G.edges[e].crossings[1]) {
+                if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) {
+                  total_cur_left[1] += c[e];
+                } else {
+                  total_cur_left[1] -= c[e];
+                }
+              }
+            }
+
+            std::array<double, 2> total_cur_right = {0.0, 0.0};
+            for (unsigned e : right_edges) {
+              if (G.edges[e].crossings[0]) {
+                if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) {
+                  total_cur_right[0] += c[e];
+                } else {
+                  total_cur_right[0] -= c[e];
+                }
+              }
+              if (G.edges[e].crossings[1]) {
+                if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) {
+                  total_cur_right[1] += c[e];
+                } else {
+                  total_cur_right[1] -= c[e];
+                }
+              }
+            }
+
+
+            if (fabs(total_cur_left[0]) < 1.0 / G.edges.size() && fabs(total_cur_left[1]) < 1.0 / G.edges.size()) {
+              for (unsigned e : left_edges) {
+                backbone[e] = true;
+                boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
+              }
+            }
+            if (fabs(total_cur_right[0]) < 1.0 / G.edges.size() && fabs(total_cur_right[1]) < 1.0 / G.edges.size()) {
+              for (unsigned e : right_edges) {
+                backbone[e] = true;
+                boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
+              }
+            }
+
+      }
+    }
+      for (unsigned e : cedges) {
+        boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG);
+      }
+  }
+  }
 }
 
-void fuse_network::break_edge(unsigned e, bool unbreak) {
+void network::break_edge(unsigned e, bool unbreak) {
   fuses[e] = !unbreak;
-  ohm_x.break_edge(e, unbreak);
-  ohm_y.break_edge(e, unbreak);
+  backbone[e] = !unbreak;
+  boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
+  boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG);
+  ds.union_set(G.dual_edges[e].v[0], G.dual_edges[e].v[1]);
+  px.break_edge(e, unbreak);
+  py.break_edge(e, unbreak);
 }
 
+
+fuse_network::fuse_network(const graph& G, cholmod_common* c, double weight) : network(G, c), weight(weight) {}
+
+fuse_network::fuse_network(const fuse_network& n) : network(n), weight(n.weight) {}
+
 current_info fuse_network::get_current_info() {
-  current_info cx = ohm_x.solve(fuses);
-  current_info cy = ohm_y.solve(fuses);
+  px.solve(fuses);
+  py.solve(fuses);
 
-  bool done_x = cx.conductivity[0] < 1.0 / G.edges.size();
-  bool done_y = cy.conductivity[1] < 1.0 / G.edges.size();
+  bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size();
+  bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size();
 
   current_info ctot;
   ctot.currents.resize(G.edges.size());
-  ctot.conductivity = {cx.conductivity[0], cy.conductivity[1]}; 
+  ctot.conductivity = {px.sol.conductivity[0], py.sol.conductivity[1]}; 
 
   if (done_x && !done_y) {
     for (unsigned i = 0; i < G.edges.size(); i++) {
-      ctot.currents[i] = fabs(weight * cy.currents[i] / cy.conductivity[1]);
+      ctot.currents[i] = fabs(weight * py.sol.currents[i] / py.sol.conductivity[1]);
     }
   } else if (done_y && !done_x) {
     for (unsigned i = 0; i < G.edges.size(); i++) {
-      ctot.currents[i] = fabs((1 - weight) * cx.currents[i] / cx.conductivity[0]);
+      ctot.currents[i] = fabs((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0]);
     }
   } else if (!done_x && !done_y) {
     for (unsigned i = 0; i < G.edges.size(); i++) {
-      ctot.currents[i] = fabs((1 - weight) * cx.currents[i] / cx.conductivity[0] +
-                                    weight * cy.currents[i] / cy.conductivity[1]);
+      ctot.currents[i] = fabs((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0] +
+                                    weight * py.sol.currents[i] / py.sol.conductivity[1]);
     }
   }
 
@@ -118,42 +504,35 @@ current_info fuse_network::get_current_info() {
 
 
 elastic_network::elastic_network(const graph& G, cholmod_common* c, double weight) :
-  network(G), weight(weight),
-  hook_x(G, 0, c), hook_y(G, 1, c) {}
+  network(G, c), weight(weight) {}
 
 elastic_network::elastic_network(const elastic_network& n) :
-  network(n), weight(n.weight), hook_x(n.hook_x), hook_y(n.hook_y) {}
-
-void elastic_network::break_edge(unsigned e, bool unbreak) {
-  fuses[e] = !unbreak;
-  hook_x.break_edge(e, unbreak);
-  hook_y.break_edge(e, unbreak);
-}
+  network(n), weight(n.weight) {}
 
 current_info elastic_network::get_current_info() {
-  current_info cx = hook_x.solve(fuses);
-  current_info cy = hook_y.solve(fuses);
+  px.solve(fuses);
+  py.solve(fuses);
 
-  bool done_x = cx.conductivity[0] < 1.0 / G.edges.size();
-  bool done_y = cy.conductivity[1] < 1.0 / G.edges.size();
+  bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size();
+  bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size();
 
   current_info ctot;
   ctot.currents.resize(G.edges.size());
-  ctot.conductivity[0] = cx.conductivity[0]; 
-  ctot.conductivity[1] = cy.conductivity[1]; 
+  ctot.conductivity[0] = px.sol.conductivity[0]; 
+  ctot.conductivity[1] = py.sol.conductivity[1]; 
 
   if (done_x && !done_y) {
     for (unsigned i = 0; i < G.edges.size(); i++) {
-      ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity[1];
+      ctot.currents[i] = weight * fabs(py.sol.currents[i]) / py.sol.conductivity[1];
     }
   } else if (done_y && !done_x) {
     for (unsigned i = 0; i < G.edges.size(); i++) {
-      ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity[0];
+      ctot.currents[i] = (1 - weight) * fabs(px.sol.currents[i]) / px.sol.conductivity[0];
     }
   } else if (!done_x && !done_y) {
     for (unsigned i = 0; i < G.edges.size(); i++) {
-      ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity[0], 2) +
-                              pow(weight * cy.currents[i] / cy.conductivity[1], 2));
+      ctot.currents[i] = sqrt(pow((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0], 2) +
+                              pow(weight * py.sol.currents[i] / py.sol.conductivity[1], 2));
     }
   }
 
@@ -161,23 +540,21 @@ current_info elastic_network::get_current_info() {
 }
 
 
-percolation_network::percolation_network(const graph& G, cholmod_common* c) :
-  network(G), px(G, 0, c), py(G, 1, c) {}
+percolation_network::percolation_network(const graph& G, cholmod_common* c) : network(G, c) {}
 
-percolation_network::percolation_network(const percolation_network& n) :
-  network(n), px(n.px), py(n.py) {}
+percolation_network::percolation_network(const percolation_network& n) : network(n) {}
 
 current_info percolation_network::get_current_info() {
   current_info ctot;
   ctot.currents.resize(G.edges.size(), 0.0);
 
-  current_info cx = px.solve(fuses);
-  current_info cy = py.solve(fuses);
+  px.solve(fuses);
+  py.solve(fuses);
 
-  ctot.conductivity = {cx.conductivity[0], cy.conductivity[1]};
+  ctot.conductivity = {px.sol.conductivity[0], py.sol.conductivity[1]};
 
   for (unsigned i = 0; i < G.edges.size(); i++) {
-    if (fabs(cx.currents[i]) > CURRENT_CUTOFF || fabs(cy.currents[i]) > CURRENT_CUTOFF) {
+    if (fabs(px.sol.currents[i]) > CURRENT_CUTOFF || fabs(py.sol.currents[i]) > CURRENT_CUTOFF) {
       ctot.currents[i] = 1.0;
     }
   }
@@ -185,9 +562,3 @@ current_info percolation_network::get_current_info() {
   return ctot;
 }
 
-void percolation_network::break_edge(unsigned e, bool unbreak) {
-  fuses[e] = !unbreak;
-  px.break_edge(e, unbreak);
-  py.break_edge(e, unbreak);
-}
-
diff --git a/lib/src/problem.cpp b/lib/src/problem.cpp
index f66a35c..0645244 100644
--- a/lib/src/problem.cpp
+++ b/lib/src/problem.cpp
@@ -75,6 +75,8 @@ problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_c
   factor = CHOL_F(analyze)(laplacian, c);
   CHOL_F(factorize)(laplacian, factor, c);
   CHOL_F(free_sparse)(&laplacian, c);
+
+  sol.currents.resize(G.edges.size());
 }
 
 problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) {
@@ -97,7 +99,7 @@ problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G,
   CHOL_F(free_triplet)(&t, c);
 }
 
-problem::problem(const problem& other) : G(other.G), axis(other.axis), c(other.c) {
+problem::problem(const problem& other) : G(other.G), axis(other.axis), c(other.c), sol(other.sol) {
   b = CHOL_F(copy_dense)(other.b, c);
   factor = CHOL_F(copy_factor)(other.factor, c);
   voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c);
@@ -109,7 +111,7 @@ problem::~problem() {
   CHOL_F(free_sparse)(&voltcurmat, c);
 }
 
-current_info problem::solve(std::vector<bool>& fuses) {
+void problem::solve(std::vector<bool>& fuses) {
   cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
 
   if (((double *)x->x)[0] != ((double *)x->x)[0]) {
@@ -122,15 +124,13 @@ current_info problem::solve(std::vector<bool>& fuses) {
   double beta[2] = {0, 0};
   CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
 
-  std::vector<double> currents(G.edges.size());
-
-  std::array<double, 2> total_current = {0, 0};
+  sol.conductivity = {0, 0};
 
   for (int i = 0; i < G.edges.size(); i++) {
     if (fuses[i]) {
-      currents[i] = 0;
+      sol.currents[i] = 0;
     } else {
-      currents[i] = ((double *)y->x)[i];
+      sol.currents[i] = ((double *)y->x)[i];
 
       graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
       graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
@@ -144,22 +144,20 @@ current_info problem::solve(std::vector<bool>& fuses) {
         }
 
         if (comp) {
-          currents[i] += 1.0;
-          total_current[axis] += currents[i];
+          sol.currents[i] += 1.0;
+          sol.conductivity[axis] += sol.currents[i];
         } else {
-          currents[i] -= 1.0;
-          total_current[axis] -= currents[i];
+          sol.currents[i] -= 1.0;
+          sol.conductivity[axis] -= sol.currents[i];
         }
       }
     }
   }
 
-  total_current[!axis] = 0.0;
+  sol.conductivity[!axis] = 0.0;
 
   CHOL_F(free_dense)(&x, c);
   CHOL_F(free_dense)(&y, c);
-
-  return {total_current, currents};
 }
 
 void problem::break_edge(unsigned e, bool unbreak) {
diff --git a/src/analysis_tools.hpp b/src/analysis_tools.hpp
index 659cacf..c3bffa2 100644
--- a/src/analysis_tools.hpp
+++ b/src/analysis_tools.hpp
@@ -1,10 +1,4 @@
 
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/connected_components.hpp>
-#include <boost/graph/depth_first_search.hpp>
-#include <boost/range/combine.hpp>
-#include <boost/foreach.hpp>
-
 #include <vector>
 #include <algorithm>
 #include <cmath>
@@ -13,14 +7,6 @@
 
 #include <network.hpp>
 
-struct EdgeProperties {
-  unsigned index;
-};
-
-typedef boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, boost::no_property, EdgeProperties> Graph;
-typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;
-typedef boost::graph_traits<Graph>::vertices_size_type VertexIndex;
-
 
 template<class T>
 bool is_shorter(const std::list<T> &, const std::list<T> &);
diff --git a/src/animate.cpp b/src/animate.cpp
index 9ba335b..532bd70 100644
--- a/src/animate.cpp
+++ b/src/animate.cpp
@@ -19,9 +19,10 @@ void animate::pre_fracture(const network &n) {
   boost::remove_edge_if(trivial, G);
 
   glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
+  glLineWidth(5);
   glClear(GL_COLOR_BUFFER_BIT);
   glBegin(GL_LINES);
-  glColor3f(0.0f, 0.0f, 0.0f);
+  glColor3f(0.9f, 0.9f, 0.9f);
   for (unsigned i = 0; i < n.G.edges.size(); i++) {
     graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r;
     graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r;
@@ -59,19 +60,62 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
 
   boost::add_edge(n.G.dual_edges[i].v[0], n.G.dual_edges[i].v[1], {i}, G);
 
+  glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
+  glClear(GL_COLOR_BUFFER_BIT);
    glBegin(GL_LINES);          // Each set of 3 vertices form a triangle
-    glColor3f(1.0f, 1.0f, 1.0f); // Blue
-    graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r;
-    graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r;
 
-    if (n.G.edges[i].crossings[0]) {
+  /*
+    glColor3f(0.0f, 0.0f, 0.0f); // Blue
+    graph::coordinate r1 = n.G.dual_vertices[n.G.dual_edges[i].v[0]].r;
+    graph::coordinate r2 = n.G.dual_vertices[n.G.dual_edges[i].v[1]].r;
+
+    if (n.G.dual_edges[i].crossings[0]) {
       if (r1.x < r2.x) {
         r1.x += n.G.L.x;
       } else {
         r2.x += n.G.L.x;
       }
     }
-    if (n.G.edges[i].crossings[1]) {
+    if (n.G.dual_edges[i].crossings[1]) {
+      if (r1.y < r2.y) {
+        r1.y += n.G.L.y;
+      } else {
+        r2.y += n.G.L.y;
+      }
+    }
+
+    glVertex2d(r1.x, r1.y);
+    glVertex2d(r2.x, r2.y);
+    */
+
+   bool weird = false;
+
+    for (unsigned j = 0; j < n.G.edges.size(); j++) {
+      bool draw = false;
+      if (cur.currents[j] < 1e-9 && !n.backbone[j]) {
+    glColor3f(1.0f, 0.0f, 0.0f); // Blue
+    weird = true;
+    draw = true;
+      } else if (n.backbone[j] && !n.fuses[j] && j != i) {
+    glColor3f(0.8f, 0.8f, 0.8f); // Blue
+    draw = true;
+      } else if (!n.fuses[j]) {
+    glColor3f(0.0f, 0.0f, 0.0f); // Blue
+    draw = true;
+      }
+
+      if (draw) {
+    graph::coordinate r1 = n.G.vertices[n.G.edges[j].v[0]].r;
+    graph::coordinate r2 = n.G.vertices[n.G.edges[j].v[1]].r;
+
+    if (n.G.edges[j].crossings[0]) {
+      if (r1.x < r2.x) {
+        r1.x += n.G.L.x;
+      } else {
+        r2.x += n.G.L.x;
+      }
+    }
+    if (n.G.edges[j].crossings[1]) {
       if (r1.y < r2.y) {
         r1.y += n.G.L.y;
       } else {
@@ -81,12 +125,17 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
 
     glVertex2d(r1.x, r1.y);
     glVertex2d(r2.x, r2.y);
+      }
+
+      }
  glEnd();
  glFlush();
+ if (weird) {std::cout << "\n"; getchar();}
 }
 
 void animate::post_fracture(network &n) {
 
+  /*
   std::list<unsigned> crack;
 //  unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]];
   unsigned crack_component = 10000;
@@ -223,5 +272,6 @@ void animate::post_fracture(network &n) {
         glFlush();
     }
   }
+*/
 }
 
diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp
index c3b4a69..9d875e6 100644
--- a/src/animate_fracture.cpp
+++ b/src/animate_fracture.cpp
@@ -31,20 +31,20 @@ int main(int argc, char* argv[]) {
   int opt;
 
   unsigned N = 1;
-  double Lx = 16.0;
-  double Ly = 16.0;
+  unsigned n = 16;
+  double a = 16.0;
   double  beta = 0.5;
 
-  while ((opt = getopt(argc, argv, "X:Y:N:b:")) != -1) {
+  while ((opt = getopt(argc, argv, "n:a:N:b:")) != -1) {
     switch (opt) {
       case 'N':
         N = (unsigned)atof(optarg);
         break;
-      case 'X':
-        Lx = atof(optarg);
+      case 'n':
+        n = atoi(optarg);
         break;
-      case 'Y':
-        Ly = atof(optarg);
+      case 'a':
+        a = atof(optarg);
         break;
       case 'b':
         beta = atof(optarg);
@@ -57,16 +57,16 @@ int main(int argc, char* argv[]) {
   cholmod_common c;
   CHOL_F(start)(&c);
 
-  animate meas(Lx, Ly, 700, argc, argv);
+  animate meas(sqrt(2*n *a), sqrt( 2*n / a), 700, argc, argv);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
 
   for (unsigned trial = 0; trial < N; trial++) {
-    graph G(Lx, Ly, rng);
+    graph G(n, a, rng);
     elastic_network elastic_network(G, &c);
     elastic_network.set_thresholds(beta, rng);
-    elastic_network.fracture(meas);
+    elastic_network.fracture(meas, true);
 
     if (quit.load())
       break;
diff --git a/src/animate_fracture_square.cpp b/src/animate_fracture_square.cpp
index 84cc8a5..bc4c3a2 100644
--- a/src/animate_fracture_square.cpp
+++ b/src/animate_fracture_square.cpp
@@ -68,7 +68,7 @@ int main(int argc, char* argv[]) {
   for (unsigned trial = 0; trial < N; trial++) {
     elastic_network tmp_network(perm_network);
     tmp_network.set_thresholds(beta, rng);
-    tmp_network.fracture(meas);
+    tmp_network.fracture(meas, false);
 
     if (quit.load())
       break;
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 7df3f67..483a3d2 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -38,8 +38,9 @@ int main(int argc, char* argv[]) {
   unsigned n = 128;
   double a = 1.0;
   bool use_aN = false;
+  double w = 0.5;
 
-  while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:w:")) != -1) {
     switch (opt) {
       case 'N':
         N = (unsigned)atof(optarg);
@@ -61,6 +62,9 @@ int main(int argc, char* argv[]) {
         a = atof(optarg);
         use_aN = true;
         break;
+      case 'w':
+        w = atof(optarg);
+        break;
       default:
         exit(1);
     }
@@ -82,9 +86,9 @@ int main(int argc, char* argv[]) {
       while (true) {
         try {
           graph G(n, a, rng);
-          elastic_network fuse_network(G, &c);
+          elastic_network fuse_network(G, &c, w);
           fuse_network.set_thresholds(beta, rng);
-          fuse_network.fracture(meas);
+          fuse_network.fracture(meas, false);
           break;
         } catch (std::exception &e) {
           std::cout << e.what() << '\n';
@@ -105,9 +109,9 @@ int main(int argc, char* argv[]) {
       while (true) {
         try {
           graph G(Lx, Ly);
-          elastic_network fuse_network(G, &c);
+          elastic_network fuse_network(G, &c, w);
           fuse_network.set_thresholds(beta, rng);
-          fuse_network.fracture(meas);
+          fuse_network.fracture(meas, false);
           break;
         } catch (std::exception &e) {
           std::cout << e.what() << '\n';
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 0732d24..c58fd84 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -153,7 +153,8 @@ ma::ma(unsigned n, double a, double beta) :
   sc(2 * n),
   sa(3 * n),
   sA(3 * n),
-  si(10000)
+  si(10000),
+  sI(10000)
 {
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_";
@@ -173,7 +174,8 @@ ma::ma(unsigned Lx, unsigned Ly, double beta) :
   sc(Lx * Ly / 2),
   sa(Lx * Ly),
   sA(Lx * Ly),
-  si(10000)
+  si(10000),
+  sI(10000)
 {
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_";
@@ -193,6 +195,7 @@ ma::~ma() {
   update_distribution_file("sa", sa, model_string);
   update_distribution_file("sA", sA, model_string);
   update_distribution_file("si", si, model_string);
+  update_distribution_file("sI", sI, model_string);
 }
 
 void ma::pre_fracture(const network&) {
@@ -212,8 +215,10 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
   }
 
   for (unsigned j = 0; j < cur.currents.size(); j++) {
-    if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0) {
+    if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0 && (!net.backbone[j] || j == i)) {
       si[(unsigned)(10000 * (logl(cur.currents[j]) + 100) / 100)]++;
+    } else if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0 && net.backbone[j] && j != i) {
+      sI[(unsigned)(10000 * (logl(cur.currents[j]) + 100) / 100)]++;
     }
   }
 
@@ -226,7 +231,13 @@ void ma::post_fracture(network &n) {
   std::vector<unsigned> component(boost::num_vertices(G));
   unsigned num = boost::connected_components(G, &component[0]);
   if (post_cracks.size() > 2 || post_cracks.size() == 0) {
-    std::cout << post_cracks.size() << "\n";
+    for (auto c : post_cracks) {
+      for (unsigned e : c.second) {
+        std::cout << e << " ";
+      }
+      std::cout << "\n";
+    }
+      getchar();
     throw badcycleex;
   }
   for (auto c : post_cracks) {
@@ -278,6 +289,7 @@ void ma::post_fracture(network &n) {
     sc[new_components[i].size() - 1]++;
   }
 
+  /*
   current_info ct = n.get_current_info();
 
 
@@ -297,6 +309,7 @@ void ma::post_fracture(network &n) {
   }
 
   sb[bb_size - 1]++;
+  */
 
   av_it++;
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 5b76e26..274a550 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -35,6 +35,7 @@ class ma : public hooks {
     std::vector<uint64_t> sA; // non-final avalanche size distribution
 
     std::vector<uint64_t> si;
+    std::vector<uint64_t> sI;
 
   public:
     long double lv;