mostly implemented the ability to find dead bonds using topological properties

6 files changed, 528 insertions, 83 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) {