successfully implemented the backbone trimming with homespun graph algorithms, measurements have mostly been disabled and need to be migrated

3 files changed, 245 insertions, 305 deletions

diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index bf9015c..8748ea9 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -18,34 +18,46 @@
 
 #define CURRENT_CUTOFF 1e-10
 
+class fuse_network;
+class elastic_network;
+class percolation_network;
+
 class network {
-  private:
-    problem px;
-    problem py;
-    std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs;
+  friend class fuse_network;
+  friend class elastic_network;
+  friend class percolation_network;
 
-    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;
+  private:
+  problem px;
+  problem py;
+  std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs;
+
+  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;
+  bool find_cycle_helper(std::array<unsigned, 2>&, const std::set<unsigned>&, unsigned, unsigned, unsigned) const;
+  std::array<unsigned, 2> find_cycle(const std::set<unsigned>&, unsigned, unsigned) const;
+  void get_cluster_edges_helper(std::set<unsigned>&, unsigned) const;
+  std::set<unsigned> get_cluster_edges(unsigned) const;
+  void get_tie_flaps_helper(std::set<unsigned>&, unsigned, unsigned) const;
+  std::list<std::set<unsigned>> get_tie_flaps(unsigned) const;
 
   public:
-    const graph& G;
+  const graph& G;
 
-    ClusterTree C;
+  ClusterTree C;
 
-    std::vector<bool> fuses;
-    std::vector<bool> backbone;
-    std::vector<long double> thresholds;
+  std::vector<bool> fuses;
+  std::vector<bool> backbone;
+  std::vector<long double> thresholds;
 
-    network(const graph&, cholmod_common*);
-    network(const network&);
+  network(const graph&, cholmod_common*);
 
-    void set_thresholds(double, std::mt19937&);
-    void fracture(hooks&, bool one_axis = true);
+  void set_thresholds(double, std::mt19937&);
+  void fracture(hooks&, bool one_axis = true);
 
-    virtual current_info get_current_info() const {current_info empty; return empty;};
+  virtual current_info get_current_info() {current_info empty; return empty;};
 };
 
 class fuse_network : public network {
@@ -56,7 +68,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() const override;
+    current_info get_current_info() override;
 };
 
 class elastic_network : public network {
@@ -67,7 +79,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() const override;
+    current_info get_current_info() override;
 };
 
 class percolation_network : public network {
@@ -75,6 +87,6 @@ class percolation_network : public network {
     percolation_network(const graph&, cholmod_common*);
     percolation_network(const percolation_network&);
 
-    current_info get_current_info() const override;
+    current_info get_current_info() override;
 };
 
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index bfdb952..a5063de 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -44,7 +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].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (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},
@@ -93,6 +93,16 @@ graph::graph(unsigned Nx, unsigned Ny) {
     }
   }
 
+  for (vertex& v : dual_vertices) {
+    v.ne.reserve(v.nd.size());
+  }
+
+  for (unsigned i = 0; i < dual_edges.size(); i++) {
+    for (unsigned vi : dual_edges[i].v) {
+      dual_vertices[vi].ne.push_back(i);
+    }
+  }
+
 }
 
 class eulerException: public std::exception
@@ -411,6 +421,16 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
     }
   }
 
+  for (vertex& v : dual_vertices) {
+    v.ne.reserve(v.nd.size());
+  }
+
+  for (unsigned i = 0; i < dual_edges.size(); i++) {
+    for (unsigned vi : dual_edges[i].v) {
+      dual_vertices[vi].ne.push_back(i);
+    }
+  }
+
   jcv_diagram_free(&diagram);
 }
 
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 2a187c5..0cf50c7 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -77,6 +77,7 @@ void network::fracture(hooks& m, bool one_axis) {
 void network::get_cycle_edges_helper(std::set<unsigned>& cycle_edges,
                                      std::set<unsigned>& seen_vertices, unsigned v_prev,
                                      unsigned v_cur) const {
+  seen_vertices.insert(v_cur);
   for (unsigned ei : G.dual_vertices[v_cur].ne) {
     if (fuses[ei]) {
       const std::array<unsigned, 2>& e = G.dual_edges[ei].v;
@@ -84,9 +85,8 @@ void network::get_cycle_edges_helper(std::set<unsigned>& cycle_edges,
 
       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);
+        if (it != seen_vertices.end()) {
+          cycle_edges.insert(ei);
         } else {
           this->get_cycle_edges_helper(cycle_edges, seen_vertices, v_cur, vn);
         }
@@ -107,16 +107,21 @@ std::set<unsigned> network::get_cycle_edges(unsigned v0) const {
 bool network::find_cycle_helper(std::array<unsigned, 2>& sig, const std::set<unsigned>& cycle_edges, unsigned vPrev, unsigned vCur, unsigned vEnd) const {
   for (unsigned ei : G.dual_vertices[vCur].ne) {
     if (fuses[ei]) {
-      const std::array<unsigned, 2>& e = G.dual_edges[ei].v;
-      unsigned vn = e[0] == vCur ? e[1] : e[0];
-      if (vn != vPrev) {
-        if (vn == vEnd) {
-          return true;
-        } else {
-          if (find_cycle_helper(sig, cycle_edges, vCur, vn, vEnd)) {
+      auto it = cycle_edges.find(ei);
+      if (it == cycle_edges.end()) {
+        const std::array<unsigned, 2>& e = G.dual_edges[ei].v;
+        unsigned vn = e[0] == vCur ? e[1] : e[0];
+        if (vn != vPrev) {
+          if (vn == vEnd) {
             if (G.dual_edges[ei].crossings[0]) sig[0]++;
             if (G.dual_edges[ei].crossings[1]) sig[1]++;
             return true;
+          } else {
+            if (this->find_cycle_helper(sig, cycle_edges, vCur, vn, vEnd)) {
+              if (G.dual_edges[ei].crossings[0]) sig[0]++;
+              if (G.dual_edges[ei].crossings[1]) sig[1]++;
+              return true;
+            }
           }
         }
       }
@@ -127,202 +132,169 @@ bool network::find_cycle_helper(std::array<unsigned, 2>& sig, const std::set<uns
 }
 
 std::array<unsigned, 2> network::find_cycle(const std::set<unsigned>& cycle_edges, unsigned v0, unsigned v1) const {
-  std::array<unsigned, 2> sig;
-  this->find_cycle_helper(sig, cycle_edges, v0, v0, 
+  std::array<unsigned, 2> sig = {0, 0};
+  this->find_cycle_helper(sig, cycle_edges, v0, v0, v1);
+  return sig;
 }
 
-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 {};
-
-    find_cycle(const graph& G, std::array<unsigned, 2>& E, unsigned end) : G(G), E(E), end(end) {}
+void network::get_cluster_edges_helper(std::set<unsigned>& seen_edges, unsigned v) const {
+  for (unsigned ei : G.vertices[v].ne) {
+    if (!backbone[ei]) {
+      auto it = seen_edges.find(ei);
+      if (it == seen_edges.end()) {
+        const std::array<unsigned, 2>& e = G.edges[ei].v;
+        unsigned vn = e[0] == v ? e[1] : e[0];
 
-    void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) {
-      if (v == end) {
-        throw done{};
+        seen_edges.insert(ei);
+        this->get_cluster_edges_helper(seen_edges, vn);
       }
     }
+  }
+}
 
-    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]++;
-    }
+std::set<unsigned> network::get_cluster_edges(unsigned v0) const {
+  std::set<unsigned> cluster_edges;
+  this->get_cluster_edges_helper(cluster_edges, v0);
+  return cluster_edges;
+}
 
-    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 network::get_tie_flaps_helper(std::set<unsigned>& added_edges, unsigned v0, unsigned vCur) const {
+  for (unsigned ei : G.vertices[vCur].ne) {
+    if (!backbone[ei]) {
+      auto it = added_edges.find(ei);
+      if (it == added_edges.end()) {
+        const std::array<unsigned, 2>& e = G.edges[ei].v;
+        unsigned vn = e[0] == vCur ? e[1] : e[0];
+
+        if (vn != v0) {
+          added_edges.insert(ei);
+          this->get_tie_flaps_helper(added_edges, v0, vn);
+        }
+      }
     }
-  };
+  }
+}
 
-  class get_cluster_edges : public boost::default_dfs_visitor {
-  public:
-    std::set<unsigned>& E;
+std::list<std::set<unsigned>> network::get_tie_flaps(unsigned v0) const {
+  std::list<std::set<unsigned>> tie_flaps;
+
+  for (unsigned ei : G.vertices[v0].ne) {
+    if (!backbone[ei]) {
+      bool seen_edge = false;
+      for (const std::set<unsigned>& flap : tie_flaps) {
+        auto it = flap.find(ei);
+        if (it != flap.end()) {
+          seen_edge = true;
+          break;
+        }
+      }
 
-    get_cluster_edges(std::set<unsigned>& E) : E(E) {}
+      if (!seen_edge) {
+        tie_flaps.push_back({ei});
 
-    void examine_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
-      unsigned e = g[ed].index;
-      E.insert(e);
+        const std::array<unsigned, 2>& e = G.edges[ei].v;
+        unsigned vn = e[0] == v0 ? e[1] : e[0];
+
+        this->get_tie_flaps_helper(tie_flaps.back(), v0, vn);
+      }
     }
-  };
+  }
+
+  return tie_flaps;
+}
+
+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();
 
+  // First, we will check for lollipops!
   for (unsigned i = 0; i < G.edges.size(); i++) {
     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 = this->get_cycle_edges(G.dual_edges[i].v[0], G.dual_edges[i].v[1]);
-
-        std::array<unsigned, 2> cycle = {0, 0};
+        // This is a candidate lollipop stem. First, we will identify any
+        // cycles in the dual cluster that impinges on the stem and mark them
+        // by an edge that uniquely severs each.
+        std::set<unsigned> cedges = this->get_cycle_edges(G.dual_edges[i].v[0]);
+
+        // Now, we create a crossing signature for each cycle. First, we do it
+        // for the new cycle that would form by removing the stem.
+        std::array<unsigned, 2> base_sig = this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[i].v[1]);
         if (G.dual_edges[i].crossings[0])
-          cycle[0]++;
+          base_sig[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);
-          }
+          base_sig[1]++;
 
-          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];
-                  }
+        // Then, we do it for each of the existing cycles we found in the
+        // previous step.
+        std::list<std::array<unsigned, 2>> all_sigs = {base_sig};
+        for (unsigned e : cedges) {
+          std::array<unsigned, 2> new_sig_1 = this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[e].v[0]);
+          std::array<unsigned, 2> new_sig_2 = this->find_cycle(cedges, G.dual_edges[i].v[1], G.dual_edges[e].v[1]);
+          std::array<unsigned, 2> new_sig = {new_sig_1[0] + new_sig_2[0], new_sig_1[1] + new_sig_2[1]};
+
+          if (G.dual_edges[i].crossings[0])
+            new_sig[0]++;
+          if (G.dual_edges[i].crossings[1])
+            new_sig[1]++;
+          if (G.dual_edges[e].crossings[0])
+            new_sig[0]++;
+          if (G.dual_edges[e].crossings[1])
+            new_sig[1]++;
+
+          all_sigs.push_back(new_sig);
+        }
+
+        // Now, having found all cycles involving the candidate stem, we check
+        // if any of them spans the torus and therefore can carry current.
+        if (std::any_of(all_sigs.begin(), all_sigs.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));
+              })) {
+          // If it does, we remove it from the backbone!
+          seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
+          backbone[i] = true;
+
+          // We're not done yet: we've severed the stem but the pop remains! We
+          // check each side of the lollipop and sum up all of the currents
+          // that span an edge. Any component without a sufficiently large net
+          // current must be disconnected from the current-carrying cluster and
+          // can also be removed from the backbone.
+          for (unsigned j = 0; j < 2; j++) {
+            std::set<unsigned> cluster_edges = this->get_cluster_edges(G.edges[i].v[j]);
+            std::array<double, 2> total_current = {0.0, 0.0};
+
+            for (unsigned e : cluster_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_current[0] += c[e];
+                } else {
+                  total_current[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_current[1] += c[e];
+                } else {
+                  total_current[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_current[0]) < 1.0 / G.edges.size() &&
+                fabs(total_current[1]) < 1.0 / G.edges.size()) {
+              for (unsigned e : cluster_edges) {
+                backbone[e] = true;
               }
             }
           }
         }
-        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;
-      }
-
-      void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor vd, const Graph& g) {
-        if (!initialized && disc > 0) {
-          start_edge = vd;
-          initialized = true;
-        }
-        disc++;
-      }
-
-      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);
-        }
-      }
-  };
-
+  // Now, we will check for bow ties!
   std::set<unsigned> bb_verts;
 
   for (unsigned i = 0; i < G.edges.size(); i++) {
@@ -385,139 +357,75 @@ void network::update_backbone(const std::vector<double>& c) {
               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::set<unsigned> cedges = this->get_cycle_edges(d1);
 
-              std::array<unsigned, 2> cycle = {0, 0};
+              std::array<unsigned, 2> base_sig = this->find_cycle(cedges, d1, d2);
               for (unsigned k = p1; k < p2; k++) {
                 if (G.dual_edges[G.vertices[v].ne[k]].crossings[0])
-                  cycle[0]++;
+                  base_sig[0]++;
                 if (G.dual_edges[G.vertices[v].ne[k]].crossings[1])
-                  cycle[1]++;
+                  base_sig[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);
+
+              // Then, we do it for each of the existing cycles we found in the
+              // previous step.
+              std::list<std::array<unsigned, 2>> all_sigs = {base_sig};
+              for (unsigned e : cedges) {
+                std::array<unsigned, 2> new_sig_1 = this->find_cycle(cedges, d1, G.dual_edges[e].v[0]);
+                std::array<unsigned, 2> new_sig_2 = this->find_cycle(cedges, d2, G.dual_edges[e].v[1]);
+                std::array<unsigned, 2> new_sig = {new_sig_1[0] + new_sig_2[0], new_sig_1[1] + new_sig_2[1]};
+
+                for (unsigned k = p1; k < p2; k++) {
+                  if (G.dual_edges[G.vertices[v].ne[k]].crossings[0])
+                    new_sig[0]++;
+                  if (G.dual_edges[G.vertices[v].ne[k]].crossings[1])
+                    new_sig[1]++;
                 }
+                if (G.dual_edges[e].crossings[0])
+                  new_sig[0]++;
+                if (G.dual_edges[e].crossings[1])
+                  new_sig[1]++;
 
-                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];
-                      }
-                    }
-                  }
+                all_sigs.push_back(new_sig);
+              }
+
+              if (std::any_of(all_sigs.begin(), all_sigs.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::array<double, 2> total_cur_right = {0.0, 0.0};
-                  for (unsigned e : right_edges) {
+                std::list<std::set<unsigned>> flaps = this->get_tie_flaps(v);
+
+                for (const std::set<unsigned>& flap : flaps) {
+                  std::array<double, 2> total_current = {0.0, 0.0};
+                  for (unsigned e : flap) {
                     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];
+                        total_current[0] += c[e];
                       } else {
-                        total_cur_right[0] -= c[e];
+                        total_current[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];
+                        total_current[1] += c[e];
                       } else {
-                        total_cur_right[1] -= c[e];
+                        total_current[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) {
+                  if (fabs(total_current[0]) < 1.0 / G.edges.size() &&
+                      fabs(total_current[1]) < 1.0 / G.edges.size()) {
+                    for (unsigned e : flap) {
                       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);
-              }
             }
           }
         }