started migration away from boost algorithms is the backbone trimming code

4 files changed, 411 insertions, 325 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index f5b6d1a..f7477fa 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -3,8 +3,8 @@ cmake_minimum_required(VERSION 3.8)
 
 project(fracture)
 
-set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall")
-set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -march=native")
+set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall -std=c++2a")
+set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -march=native -std=c++2a")
 
 set (CMAKE_CXX_STANDARD 17)
 set (CMAKE_C_STANDARD 11)
diff --git a/lib/include/clusters.hpp b/lib/include/clusters.hpp
new file mode 100644
index 0000000..4cc9073
--- /dev/null
+++ b/lib/include/clusters.hpp
@@ -0,0 +1,52 @@
+
+#include <vector>
+#include "graph.hpp"
+
+class ClusterTree {
+private:
+  std::vector<signed> p;
+  std::vector<unsigned> o;
+
+public:
+  std::vector<unsigned> c;
+
+  ClusterTree(unsigned n) : p(n, -1), o(n, 1), c(n, 0) { c[0] = n; }
+
+  unsigned findroot(unsigned i) {
+    if (p[i] < 0)
+      return i;
+    else
+      return p[i] = this->findroot(p[i]);
+  }
+
+  void join(unsigned i, unsigned j) {
+    c[o[i] - 1]--;
+    c[o[j] - 1]--;
+
+    if (o[i] < o[j]) {
+      p[i] = j;
+      o[j] += o[i];
+      c[o[j] - 1]++;
+    } else {
+      p[j] = i;
+      o[i] += o[j];
+      c[o[i] - 1]++;
+    }
+  }
+
+  void add_bond(const graph::edge& b) {
+    unsigned i = this->findroot(b.v[0]);
+    unsigned j = this->findroot(b.v[1]);
+
+    if (i != j) {
+      this->join(i, j);
+    }
+  }
+
+  bool same_component(unsigned v0, unsigned v1) {
+    unsigned i = this->findroot(v0);
+    unsigned j = this->findroot(v1);
+
+    return i == j;
+  }
+};
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index 8211cfc..bf9015c 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -14,53 +14,38 @@
 #include "hooks.hpp"
 #include "current_info.hpp"
 #include "array_hash.hpp"
+#include "clusters.hpp"
 
 #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;
-};
+class network {
+  private:
+    problem px;
+    problem py;
+    std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs;
 
-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;
+    void update_backbone(const std::vector<double>& c);
+    void break_edge(unsigned, bool unbreak = false);
+    void get_cycle_edges_helper(std::set<unsigned>&, std::set<unsigned>&, unsigned, unsigned) const;
+    std::set<unsigned> get_cycle_edges(unsigned) const;
+    std::array<unsigned, 2> find_cycle(const std::set<unsigned>&, unsigned, unsigned) const;
 
-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;
+    ClusterTree C;
 
     std::vector<bool> fuses;
     std::vector<bool> backbone;
     std::vector<long double> thresholds;
 
-    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);
+    void fracture(hooks&, bool one_axis = true);
 
-    void break_edge(unsigned, bool unbreak = false);
-    virtual current_info get_current_info() {current_info empty; return empty;};
+    virtual current_info get_current_info() const {current_info empty; return empty;};
 };
 
 class fuse_network : public network {
@@ -71,7 +56,7 @@ class fuse_network : public network {
     fuse_network(const graph&, cholmod_common*, double weight = 1.0);
     fuse_network(const fuse_network&);
 
-    current_info get_current_info() override;
+    current_info get_current_info() const override;
 };
 
 class elastic_network : public network {
@@ -82,7 +67,7 @@ class elastic_network : public network {
     elastic_network(const graph&, cholmod_common*, double weight = 0.5);
     elastic_network(const elastic_network&);
 
-    current_info get_current_info() override;
+    current_info get_current_info() const override;
 };
 
 class percolation_network : public network {
@@ -90,6 +75,6 @@ class percolation_network : public network {
     percolation_network(const graph&, cholmod_common*);
     percolation_network(const percolation_network&);
 
-    current_info get_current_info() override;
+    current_info get_current_info() const override;
 };
 
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index ed23d45..6b3873b 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -1,28 +1,13 @@
 
 #include "network.hpp"
 
-class nofuseException: public std::exception
-{
-  virtual const char* what() const throw()
-  {
-    return "No valid fuse was available to break.";
-  }
+class nofuseException : public std::exception {
+  virtual const char* what() const throw() { return "No valid fuse was available to break."; }
 } nofuseex;
 
-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), 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);
-}
+network::network(const graph& G, cholmod_common* c)
+    : G(G), C(G.dual_vertices.size()), fuses(G.edges.size()), backbone(G.edges.size()),
+      thresholds(G.edges.size()), px(G, 0, c), py(G, 1, c) {}
 
 void network::set_thresholds(double beta, std::mt19937& rng) {
   if (beta == 0.0) {
@@ -78,7 +63,7 @@ void network::fracture(hooks& m, bool one_axis) {
       }
     }
 
-    if (max_pos == UINT_MAX)  {
+    if (max_pos == UINT_MAX) {
       throw nofuseex;
     }
 
@@ -89,155 +74,199 @@ void network::fracture(hooks& m, bool one_axis) {
   m.post_fracture(*this);
 }
 
-void network::update_backbone(const std::vector<double>& c) {
+void network::get_cycle_edges_helper(std::set<unsigned>& cycle_edges,
+                                     std::set<unsigned>& seen_vertices, unsigned v_prev,
+                                     unsigned v_cur) const {
+  for (unsigned ei : G.dual_vertices[v_cur].ne) {
+    const std::array<unsigned, 2>& e = G.dual_edges[ei].v;
+    unsigned vn = e[0] == v_cur ? e[1] : e[0];
+
+    if (vn != v_prev) {
+      auto it = seen_vertices.find(vn);
+      if (it == seen_vertices.end()) {
+        //    if (seen_vertices.contains()) { // need to wait for better c++20 support...
+        cycle_edges.insert(vn);
+      } else {
+        this->get_cycle_edges_helper(cycle_edges, seen_vertices, v_cur, vn);
+      }
+    }
+  }
+}
 
-        bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size();
-        bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size();
-  class get_cycle_edges : public boost::default_dfs_visitor {
-    public:
-      std::set<unsigned> tE;
-      std::list<unsigned>& bE;
+std::set<unsigned> network::get_cycle_edges(unsigned v0) const {
+  std::set<unsigned> seen_vertices;
+  std::set<unsigned> cycle_edges;
 
-      get_cycle_edges(std::list<unsigned>& bE) : bE(bE) {};
+  this->get_cycle_edges_helper(cycle_edges, seen_vertices, v0, v0);
 
-      void tree_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
-        tE.insert(g[e].index);
-      }
+  return cycle_edges;
+}
 
-      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);
-        }
-      }
-  };
+std::array<unsigned, 2> network::find_cycle(const std::set<unsigned>& cycle_edges, unsigned v0, unsigned v1) const {
+  for (unsigned ei : G.dual_vertices[v0].ne) {
+    const std::array<unsigned, 2>& e = G.dual_edges[ei].v;
+    unsigned vn = e[0] == v0 ? e[1] : e[0];
+    if (vn == v1) {
+      return {ei};
+    } else {
+      std::list<unsigned> edges = this->get_cycle_edges(vn, v1);
+      edges.splice(edges.end(), {ei});
+      return edges;
+    }
+  }
+
+  return {};
+}
+
+void network::update_backbone(const std::vector<double>& c) {
+  bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size();
+  bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size();
 
   class find_cycle : public boost::default_dfs_visitor {
-    public:
-      const graph& G;
-      std::array<unsigned, 2>& E;
-      unsigned end;
-      struct done{};
+  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) {}
+    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 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 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]--;
-      }
+    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;
+  public:
+    std::set<unsigned>& E;
 
-      get_cluster_edges(std::set<unsigned>& E) : E(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);
-      }
+    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)) {
+    if ((!backbone[i]) && C.same_component(G.dual_edges[i].v[0], G.dual_edges[i].v[1])) {
       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&) {
+        std::list<unsigned> cedges = this->get_cycle_edges(G.dual_edges[i].v[0], G.dual_edges[i].v[1]);
+
+        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);
           }
-          other_cycles.push_back(tcycle);
-        }
 
-        if (std::any_of(other_cycles.begin(), other_cycles.end(), [done_x, done_y](std::array<unsigned, 2> c){return (c[0] % 2 == 0 && c[1] % 2 == 0) || ((c[0] % 2 == 0 && done_x) || (c[1] % 2 == 0 && done_y));})) {
-          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 (std::any_of(other_cycles.begin(), other_cycles.end(),
+                [done_x, done_y](std::array<unsigned, 2> c) {
+                return (c[0] % 2 == 0 && c[1] % 2 == 0) ||
+                ((c[0] % 2 == 0 && done_x) || (c[1] % 2 == 0 && done_y));
+                })) {
+            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 (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);
+              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);
+        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 {
@@ -279,8 +308,6 @@ void network::update_backbone(const std::vector<double>& c) {
           RE.insert(e);
         }
       }
-
-
   };
 
   std::set<unsigned> bb_verts;
@@ -297,161 +324,189 @@ void network::update_backbone(const std::vector<double>& c) {
     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 l = (i + j) % G.vertices[v].nd.size();
+        if (C.same_component(G.vertices[v].nd[i], G.vertices[v].nd[l])) {
           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;
-          unsigned ie1 = 0;
-          unsigned ie2 = 0;
-
-          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]]) {
-                ie1++;
+            bool any_intervening1 = false;
+            bool any_intervening_e_1 = false;
+            bool any_intervening2 = false;
+            bool any_intervening_e_2 = false;
+            unsigned ie1 = 0;
+            unsigned ie2 = 0;
+
+            for (unsigned k = il + 1; k < ig; k++) {
+              if (!C.same_component(G.vertices[v].nd[i], G.vertices[v].nd[k])) {
+                any_intervening1 = 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()]]) {
-                ie2++;
+            for (unsigned k = (ig + 1); k % G.vertices[v].nd.size() != il; k++) {
+              if (!C.same_component(G.vertices[v].nd[i],
+                    G.vertices[v].nd[k % G.vertices[v].nd.size()])) {
+                any_intervening2 = true;
+                break;
               }
             }
-          }
-          if ((any_intervening1 && any_intervening2) || (any_intervening1 && ie2 > 1) || (any_intervening2 && ie1 > 1)) {
-            found_pair = true;
-            p1 = il;
-            p2 = ig;
-            d1 = G.vertices[v].nd[il];
-            d2 = G.vertices[v].nd[ig];
-
-      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(), [done_x, done_y](std::array<unsigned, 2> c){return ((c[0] % 2 == 0 && c[1] % 2 == 0) || (c[0] % 2 == 0 && done_x)) || (c[1] % 2 == 0 && done_y);})) {
-          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 (any_intervening2 && !any_intervening1) {
+              for (unsigned k = il + 1; k <= ig; k++) {
+                if (!backbone[G.vertices[v].ne[k]]) {
+                  ie1++;
                 }
               }
-              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];
+            }
+            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()]]) {
+                  ie2++;
                 }
               }
             }
+            if ((any_intervening1 && any_intervening2) || (any_intervening1 && ie2 > 1) ||
+                (any_intervening2 && ie1 > 1)) {
+              found_pair = true;
+              p1 = il;
+              p2 = ig;
+              d1 = G.vertices[v].nd[il];
+              d2 = G.vertices[v].nd[ig];
+
+              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<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];
-                }
+              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]++;
               }
-              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];
+              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 (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 (std::any_of(other_cycles.begin(), other_cycles.end(),
+                      [done_x, done_y](std::array<unsigned, 2> c) {
+                      return ((c[0] % 2 == 0 && c[1] % 2 == 0) ||
+                          (c[0] % 2 == 0 && done_x)) ||
+                      (c[1] % 2 == 0 && done_y);
+                      })) {
+                  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);
+                    }
+                  }
+                }
               }
-            }
-            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);
               }
             }
-
-      }
-    }
-      for (unsigned e : cedges) {
-        boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG);
-      }
-  }
-  }
+          }
         }
       }
     }
@@ -461,15 +516,13 @@ void network::update_backbone(const std::vector<double>& c) {
 void network::break_edge(unsigned e, bool unbreak) {
   fuses[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]);
+  C.add_bond(G.dual_edges[e]);
   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 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) {}
 
@@ -482,7 +535,7 @@ current_info fuse_network::get_current_info() {
 
   current_info ctot;
   ctot.currents.resize(G.edges.size());
-  ctot.conductivity = {px.sol.conductivity[0], py.sol.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++) {
@@ -495,19 +548,17 @@ current_info fuse_network::get_current_info() {
   } else if (!done_x && !done_y) {
     for (unsigned i = 0; i < G.edges.size(); i++) {
       ctot.currents[i] = fabs((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0] +
-                                    weight * py.sol.currents[i] / py.sol.conductivity[1]);
+                              weight * py.sol.currents[i] / py.sol.conductivity[1]);
     }
   }
 
   return ctot;
 }
 
+elastic_network::elastic_network(const graph& G, cholmod_common* c, double weight)
+    : network(G, c), weight(weight) {}
 
-elastic_network::elastic_network(const graph& G, cholmod_common* c, double weight) :
-  network(G, c), weight(weight) {}
-
-elastic_network::elastic_network(const elastic_network& n) :
-  network(n), weight(n.weight) {}
+elastic_network::elastic_network(const elastic_network& n) : network(n), weight(n.weight) {}
 
 current_info elastic_network::get_current_info() {
   px.solve(fuses);
@@ -518,8 +569,8 @@ current_info elastic_network::get_current_info() {
 
   current_info ctot;
   ctot.currents.resize(G.edges.size());
-  ctot.conductivity[0] = px.sol.conductivity[0]; 
-  ctot.conductivity[1] = py.sol.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++) {
@@ -539,7 +590,6 @@ current_info elastic_network::get_current_info() {
   return ctot;
 }
 
-
 percolation_network::percolation_network(const graph& G, cholmod_common* c) : network(G, c) {}
 
 percolation_network::percolation_network(const percolation_network& n) : network(n) {}
@@ -558,7 +608,6 @@ current_info percolation_network::get_current_info() {
       ctot.currents[i] = 1.0;
     }
   }
-  
+
   return ctot;
 }
-