made backbone cleaning more efficient by restricting to one-side-broken only, updated the crack finder to this new paradigm

2 files changed, 71 insertions, 15 deletions

diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index 29bf55a..5b474bb 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -35,9 +35,10 @@ private:
   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;
+  bool get_cycle_helper(std::array<unsigned, 2>&, std::set<unsigned>&, 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;
@@ -56,7 +57,10 @@ public:
   network(const graph&, cholmod_common*);
 
   void set_thresholds(double, std::mt19937&);
-  void fracture(hooks&, bool one_axis = true);
+  void fracture(hooks&);
+
+  std::set<unsigned> get_cycle_edges(unsigned) const;
+  std::pair<std::array<unsigned, 2>, std::set<unsigned>> get_cycle(const std::set<unsigned>&, unsigned, unsigned) const;
 
   virtual current_info get_current_info() {
     current_info empty;
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index eef7e9a..3e460a9 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -6,8 +6,8 @@ class nofuseException : public std::exception {
 } nofuseex;
 
 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) {}
+    : px(G, 0, c), py(G, 1, c), G(G), C(G.dual_vertices.size()), fuses(G.edges.size()),
+      backbone(G.edges.size()), thresholds(G.edges.size()) {}
 
 void network::set_thresholds(double beta, std::mt19937& rng) {
   if (beta == 0.0) {
@@ -29,25 +29,22 @@ void network::set_thresholds(double beta, std::mt19937& rng) {
   }
 }
 
-void network::fracture(hooks& m, bool one_axis) {
+void network::fracture(hooks& m) {
   m.pre_fracture(*this);
 
   while (true) {
-    double min_cond = 1.0 / G.edges.size();
+    const 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) {
+    if (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond) {
       break;
     }
 
-    if (one_axis && (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond)) {
-      break;
+    if (px.sol.conductivity[0] > min_cond) {
+      this->update_backbone(px.sol.currents);
+    } else {
+      this->update_backbone(py.sol.currents);
     }
 
     unsigned max_pos = UINT_MAX;
@@ -144,6 +141,48 @@ std::array<unsigned, 2> network::find_cycle(const std::set<unsigned>& cycle_edge
   return sig;
 }
 
+bool network::get_cycle_helper(std::array<unsigned, 2>& sig, std::set<unsigned>& edges, const std::set<unsigned>& cycle_edges,
+                                unsigned vPrev, unsigned vCur, unsigned vEnd) const {
+  for (unsigned ei : G.dual_vertices[vCur].ne) {
+    if (fuses[ei]) {
+      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) {
+            edges.insert(ei);
+            if (G.dual_edges[ei].crossings[0])
+              sig[0]++;
+            if (G.dual_edges[ei].crossings[1])
+              sig[1]++;
+            return true;
+          } else {
+            if (this->get_cycle_helper(sig, edges, cycle_edges, vCur, vn, vEnd)) {
+              edges.insert(ei);
+              if (G.dual_edges[ei].crossings[0])
+                sig[0]++;
+              if (G.dual_edges[ei].crossings[1])
+                sig[1]++;
+              return true;
+            }
+          }
+        }
+      }
+    }
+  }
+
+  return false;
+}
+
+std::pair<std::array<unsigned, 2>, std::set<unsigned>> network::get_cycle(const std::set<unsigned>& cycle_edges, unsigned v0,
+                                            unsigned v1) const {
+  std::set<unsigned> edges;
+  std::array<unsigned, 2> sig = {0, 0};
+  this->get_cycle_helper(sig, edges, cycle_edges, v0, v0, v1);
+  return {sig, edges};
+}
+
 void network::get_cluster_edges_helper(std::set<unsigned>& seen_edges, unsigned v) const {
   for (unsigned ei : G.vertices[v].ne) {
     if (!backbone[ei]) {
@@ -264,7 +303,7 @@ void network::update_backbone(const std::vector<double>& 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!
+          // If none 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;
 
@@ -301,6 +340,14 @@ void network::update_backbone(const std::vector<double>& c) {
               }
             }
           }
+        } else if (std::any_of(all_sigs.begin(), all_sigs.end(),
+                        [](std::array<unsigned, 2> c) {
+                          return ((c[0] % 2 == 0 && c[1] % 2 == 1) || (c[0] % 2 == 1 && c[1] % 2 == 0));
+                        })) {
+          // If the bond separates a wrapping path, breaking it would end the
+          // fracture and therefore it will never be removed from the backbone
+          // in this manner. We can skip it in the future.
+          seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
         }
       }
     }
@@ -450,6 +497,11 @@ void network::update_backbone(const std::vector<double>& c) {
                     }
                   }
                 }
+              } else if (std::any_of(all_sigs.begin(), all_sigs.end(),
+                              [](std::array<unsigned, 2> c) {
+                                return ((c[0] % 2 == 0 && c[1] % 2 == 1) || (c[0] % 2 == 1 && c[1] % 2 == 0));
+                              })) {
+                seen_pairs[{d1, d2}] = true;
               }
             }
           }