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

7 files changed, 144 insertions, 254 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;
               }
             }
           }
diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp
index 5f180ee..2965915 100644
--- a/src/analysis_tools.cpp
+++ b/src/analysis_tools.cpp
@@ -1,125 +1,72 @@
 
 #include "analysis_tools.hpp"
 
-template <class T>
-bool is_shorter(const std::list<T> &l1, const std::list<T> &l2) {
-  return l1.size() < l2.size();
-}
-
-/*
-bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>) {
-  return true;
-}
-
-std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> find_minimal_crack(const Graph& G, const network& n) {
-  Graph Gtmp(n.G.vertices.size());
-  std::list<unsigned> removed_edges;
-
-  class add_tree_edges : public boost::default_dfs_visitor {
-    public:
-      Graph& G;
-      std::list<unsigned>& E;
-
-      add_tree_edges(Graph& G, std::list<unsigned>& E) : G(G), E(E) {}
-
-      void tree_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
-        boost::add_edge(boost::source(e, g), boost::target(e, g), g[e], G);
-      }
-
-      void back_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
-        if (!(boost::edge(boost::source(e, g), boost::target(e, g), G).second)) {
-          E.push_back(g[e].index);
-        }
-      }
-  };
-
-  add_tree_edges ate(Gtmp, removed_edges);
-  boost::depth_first_search(G, visitor(ate));
-
-  class find_cycle : public boost::default_dfs_visitor {
-    public:
-      std::list<unsigned>& E;
-      unsigned end;
-      struct done{};
-
-      find_cycle(std::list<unsigned>& E, unsigned end) : 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 e, const Graph& g) {
-        E.push_back(g[e].index);
-      }
-
-      void finish_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
-        E.erase(std::find(E.begin(), E.end(), g[e].index));
-      }
-  };
-
-  std::list<std::list<unsigned>> cycles;
+std::pair<std::array<unsigned, 2>, std::set<unsigned>> find_minimal_crack(const network& n, unsigned i0) {
+  std::set<unsigned> removed_edges = n.get_cycle_edges(n.G.dual_edges[i0].v[0]);
+
+  std::vector<std::pair<std::array<unsigned, 2>, std::set<unsigned>>> cycles;
+  cycles.reserve(removed_edges.size());
+
+  for (unsigned edge : removed_edges) {
+    cycles.push_back(n.get_cycle(removed_edges, n.G.dual_edges[edge].v[0], n.G.dual_edges[edge].v[1])); 
+    cycles.back().second.insert(edge);
+    if (n.G.dual_edges[edge].crossings[0])
+      cycles.back().first[0]++;
+    if (n.G.dual_edges[edge].crossings[1])
+      cycles.back().first[1]++;
+  }
 
-  for (auto edge : removed_edges) {
-    std::list<unsigned> cycle = {edge};
-    find_cycle vis(cycle, n.G.dual_edges[edge].v[1]);
-    std::vector<boost::default_color_type> new_color_map(boost::num_vertices(Gtmp));
-    try {
-    boost::depth_first_visit(Gtmp, n.G.dual_edges[edge].v[0], vis, boost::make_iterator_property_map(new_color_map.begin(), boost::get(boost::vertex_index, Gtmp), new_color_map[0]));
-    } catch(find_cycle::done const&) {
-      cycles.push_back(cycle);
-    }
+  if (cycles.size() == 1) {
+    return cycles.front();
   }
 
-  std::list<std::valarray<uint8_t>> bool_cycles;
-  for (auto cycle : cycles) {
-    std::valarray<uint8_t> bool_cycle(n.G.edges.size());
-    for (auto v : cycle) {
-      bool_cycle[v] = 1;
+  bool all_good = true;
+  unsigned not_good;
+  for (unsigned i = 0; i < 2; i++) {
+    if (!((cycles[i].first[0] % 2 == 0 && cycles[i].first[1] % 2 == 1) || (cycles[i].first[0] % 2 == 1 && cycles[i].first[1] % 2 == 0))) {
+      all_good = false;
+      not_good = i;
     }
-    bool_cycles.push_back(bool_cycle);
   }
 
-  // generate all possible cycles by taking xor of the edge sets of the known cycles
-  for (auto it1 = bool_cycles.begin(); it1 != std::prev(bool_cycles.end()); it1++) {
-    for (auto it2 = std::next(it1); it2 != bool_cycles.end(); it2++) {
-      std::valarray<uint8_t> new_bool_cycle = (*it1) ^ (*it2);
-      std::list<unsigned> new_cycle;
-      unsigned pos = 0;
-      for (uint8_t included : new_bool_cycle) {
-        if (included) {
-          new_cycle.push_back(pos);
-        }
-        pos++;
-      }
-      cycles.push_back(new_cycle);
+  if (all_good) {
+    if (cycles[0].second.size() > cycles[1].second.size()) {
+      return cycles[1];
+    } else {
+      return cycles[0];
     }
   }
 
-  std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> output;
+  const auto& cycle1 = cycles[0];
+  const auto& cycle2 = cycles[1];
 
-  // find the cycle representing the crack by counting boundary crossings
-  for (auto cycle : cycles) {
-    std::array<unsigned, 2> crossing_count{0,0};
+  unsigned sum_sig_0 = cycle1.first[0] + cycle2.first[0];
+  unsigned sum_sig_1 = cycle1.first[1] + cycle2.first[1];
 
-    for (auto edge : cycle) {
-      if (n.G.dual_edges[edge].crossings[0]) {
-        crossing_count[0]++;
-      }
-      if (n.G.dual_edges[edge].crossings[1]) {
-        crossing_count[1]++;
-      }
-    }
+  std::valarray<uint8_t> bool_cycle1(n.G.edges.size());
+  std::valarray<uint8_t> bool_cycle2(n.G.edges.size());
+  for (auto v : cycle1.second) {
+    bool_cycle1[v] = 1;
+  }
+  for (auto v : cycle2.second) {
+    bool_cycle2[v] = 1;
+  }
 
-    if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) {
-      output.push_back({{1, 0}, cycle});
-    } else if (crossing_count[0] % 2 == 0 && crossing_count[1] % 2 == 1) {
-      output.push_back({{0, 1}, cycle});
+  std::valarray<uint8_t> new_bool_cycle = bool_cycle1 ^ bool_cycle2;
+
+  std::set<unsigned> new_cycle;
+  unsigned pos = 0;
+  for (uint8_t included : new_bool_cycle) {
+    if (included) {
+      new_cycle.insert(pos);
     }
+    pos++;
   }
 
-  return output;
+  if (cycles[!not_good].second.size() > new_cycle.size()) {
+    return {{sum_sig_0, sum_sig_1}, new_cycle};
+  } else {
+    return cycles[!not_good];
+  }
 }
-*/
 
diff --git a/src/analysis_tools.hpp b/src/analysis_tools.hpp
index 2a4c5d2..72d1e82 100644
--- a/src/analysis_tools.hpp
+++ b/src/analysis_tools.hpp
@@ -7,11 +7,5 @@
 
 #include <network.hpp>
 
-
-template<class T>
-bool is_shorter(const std::list<T> &, const std::list<T> &);
-
-//bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>);
-
-//std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> find_minimal_crack(const Graph &, const network &);
+std::pair<std::array<unsigned, 2>, std::set<unsigned>> find_minimal_crack(const network &, unsigned i0);
 
diff --git a/src/animate.cpp b/src/animate.cpp
index 6c844d2..81bf4ac 100644
--- a/src/animate.cpp
+++ b/src/animate.cpp
@@ -1,5 +1,6 @@
 
 #include "animate.hpp"
+#include <unistd.h>
 #include <iostream>
 
 animate::animate(double Lx, double Ly, unsigned window_size, int argc, char *argv[]) {
@@ -139,65 +140,27 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
       }
  glEnd();
  glFlush();
- if (nw > 2 && !seen_guy) {seen_guy = true; getchar();}
+// if (nw > 2 && !seen_guy) {seen_guy = true; getchar();}
 }
 
 void animate::post_fracture(network &n) {
 
-  /*
-  std::list<unsigned> crack;
-//  unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]];
-  unsigned crack_component = 10000;
-
-  std::default_random_engine gen;
-  std::uniform_real_distribution<double> dis(0.0,1.0);
-
-  auto av_it = avalanches.rbegin();
-
-  while (true) {
-    for (unsigned e : *av_it) {
-      boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G);
-      n.fuses[e] = false;
-    }
-
-    auto cracks = find_minimal_crack(G, n);
-
-    if (cracks.size() == 0) {
-      break;
-    }
+  auto crack = find_minimal_crack(n, avalanches.back().back());
 
-    av_it++;
-  }
-
-  std::vector<unsigned> component(boost::num_vertices(G));
-  unsigned num = boost::connected_components(G, &component[0]);
-
-  std::vector<std::list<unsigned>> components(num);
-
-  for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
-    components[component[i]].push_back(i);
-  }
-
-
-  char key;
-  while ((key = getchar()) != 'n') {
-    glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
-    glClear(GL_COLOR_BUFFER_BIT);
-    glBegin(GL_LINES);
-    glColor3f(0.0f, 0.0f, 0.0f);
-    for (unsigned i = 0; i < n.G.edges.size(); i++) {
-      if (!n.fuses[i]) {
-        graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r;
-        graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r;
+  for (unsigned j : crack.second) {
+  glBegin(GL_LINES);
+  glColor3d(1.0, 0.0, 0.0);
+    graph::coordinate r1 = n.G.dual_vertices[n.G.dual_edges[j].v[0]].r;
+    graph::coordinate r2 = n.G.dual_vertices[n.G.dual_edges[j].v[1]].r;
 
-    if (n.G.edges[i].crossings[0]) {
+    if (n.G.dual_edges[j].crossings[0]) {
       if (r1.x < r2.x) {
         r1.x += n.G.L.x;
       } else {
         r2.x += n.G.L.x;
       }
     }
-    if (n.G.edges[i].crossings[1]) {
+    if (n.G.dual_edges[j].crossings[1]) {
       if (r1.y < r2.y) {
         r1.y += n.G.L.y;
       } else {
@@ -205,82 +168,12 @@ void animate::post_fracture(network &n) {
       }
     }
 
-        glVertex2d(r1.x, r1.y);
-        glVertex2d(r2.x, r2.y);
-      }
-    }
-    glEnd();
-
-    switch (key) {
-      case 's' :
-        for (auto edge : crack) {
-          glColor3f(1.0f, 0.0f, 0.0f);
-         glBegin(GL_POLYGON);
-          for (const graph::coordinate &r : n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon) {
-            glVertex2d(r.x, r.y);
-          }
-          glEnd();
-         glBegin(GL_POLYGON);
-          for (const graph::coordinate &r : n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon) {
-            glVertex2d(r.x, r.y);
-          }
-          glEnd();
-        }
-        glFlush();
-        break;
-
-      case 'c' :
-        for (unsigned i = 0; i < num; i++) {
-          if (i == crack_component) {
-            glColor3d(1.0, 0.0, 0.0);
-          } else {
-            glColor3d(dis(gen), dis(gen), dis(gen));
-          }
-
-          for (auto it = components[i].begin(); it != components[i].end(); it++) {
-        glBegin(GL_POLYGON);          // Each set of 3 vertices form a triangle
-            for (const graph::coordinate &r: n.G.dual_vertices[*it].polygon) {
-              glVertex2d(r.x, r.y);
-            }
-        glEnd();
-          }
-        }
-        glFlush();
-        break;
-
-      case 'C' :
-        for (unsigned i = 0; i < num; i++) {
-          if (components[i].size() > 1) {
-            if (i == crack_component) {
-              glColor3d(1.0, 0.0, 0.0);
-            } else {
-              glColor3d(dis(gen), dis(gen), dis(gen));
-            }
-
-            for (auto it = components[i].begin(); it != components[i].end(); it++) {
-        glBegin(GL_POLYGON);          // Each set of 3 vertices form a triangle
-              for (const graph::coordinate &r :n.G.dual_vertices[*it].polygon) {
-                glVertex2d(r.x, r.y);
-              }
-        glEnd();
-            }
-          }
-        }
-        glFlush();
-        break;
-
-      case 'm' :
-        for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) {
-        glBegin(GL_POLYGON);
-        glColor3d(1.0, 0.0, 0.0);
-          for (const graph::coordinate &r :n.G.dual_vertices[*it].polygon) {
-            glVertex2d(r.x, r.y);
-          }
-        glEnd();
-        }
-        glFlush();
-    }
+    glVertex2d(r1.x, r1.y);
+    glVertex2d(r2.x, r2.y);
+  glEnd();
   }
-*/
+  glFlush();
+
+  getchar();
 }
 
diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp
index 9d875e6..2b5e55c 100644
--- a/src/animate_fracture.cpp
+++ b/src/animate_fracture.cpp
@@ -57,7 +57,7 @@ int main(int argc, char* argv[]) {
   cholmod_common c;
   CHOL_F(start)(&c);
 
-  animate meas(sqrt(2*n *a), sqrt( 2*n / a), 700, argc, argv);
+  animate meas(sqrt(2*n *a), sqrt( 2*n / a), 1000, argc, argv);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
@@ -66,7 +66,7 @@ int main(int argc, char* argv[]) {
     graph G(n, a, rng);
     elastic_network elastic_network(G, &c);
     elastic_network.set_thresholds(beta, rng);
-    elastic_network.fracture(meas, true);
+    elastic_network.fracture(meas);
 
     if (quit.load())
       break;
diff --git a/src/animate_fracture_square.cpp b/src/animate_fracture_square.cpp
index 12de7da..ad43254 100644
--- a/src/animate_fracture_square.cpp
+++ b/src/animate_fracture_square.cpp
@@ -57,7 +57,7 @@ int main(int argc, char* argv[]) {
   cholmod_common c;
   CHOL_F(start)(&c);
 
-  animate meas(Lx, Ly, 512, argc, argv);
+  animate meas(Lx, Ly, 1000, argc, argv);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
@@ -68,7 +68,7 @@ int main(int argc, char* argv[]) {
   for (unsigned trial = 0; trial < N; trial++) {
     elastic_network tmp_network(perm_network);
     tmp_network.set_thresholds(beta, rng);
-    tmp_network.fracture(meas, true);
+    tmp_network.fracture(meas);
 
     if (quit.load())
       break;