Implemented new behavior when w=0 or w=1.

5 files changed, 157 insertions, 75 deletions

diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index d562530..37f1912 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -31,9 +31,10 @@ class network {
 private:
   problem px;
   problem py;
-  std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs;
+  std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs_x;
+  std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs_y;
 
-  void update_backbone(const std::vector<double>& c);
+  void update_backbone(const std::vector<double>& c, const std::array<double, 2>& cond);
   void break_edge(unsigned, bool unbreak = false);
   void get_cycle_edges_helper(std::set<unsigned>&, std::set<unsigned>&, unsigned, unsigned) const;
   bool find_cycle_helper(std::array<unsigned, 2>&, const std::set<unsigned>&, unsigned, unsigned,
@@ -47,6 +48,7 @@ private:
   std::list<std::set<unsigned>> get_tie_flaps(unsigned) const;
 
 public:
+  bool two_sides;
   const graph& G;
 
   ClusterTree C;
@@ -55,7 +57,7 @@ public:
   std::vector<bool> backbone;
   std::vector<long double> thresholds;
 
-  network(const graph&, cholmod_common*);
+  network(const graph&, bool two_sides, cholmod_common*);
 
   void set_thresholds(double, std::mt19937&);
   void fracture(hooks&);
@@ -76,7 +78,7 @@ private:
   double weight;
 
 public:
-  fuse_network(const graph&, cholmod_common*, double weight = 1.0);
+  fuse_network(const graph&, cholmod_common*);
   fuse_network(const fuse_network&);
 
   current_info get_current_info() override;
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 561796b..8944c15 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -1,14 +1,16 @@
 
 #include "network.hpp"
 #include <sstream>
+#include <iostream>
 
 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)
-    : 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()) {}
+network::network(const graph& G, bool two_sides, cholmod_common* c)
+    : px(G, 0, c), py(G, 1, c), two_sides(two_sides), 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) {
@@ -38,14 +40,16 @@ void network::fracture(hooks& m) {
 
     current_info c = this->get_current_info();
 
-    if (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond) {
+    if ((c.conductivity[0] < min_cond || c.conductivity[1] < min_cond) && two_sides) {
+      break;
+    } else if (c.conductivity[0] < min_cond && !two_sides) {
       break;
     }
 
     if (px.sol.conductivity[0] > min_cond) {
-      this->update_backbone(px.sol.currents);
+      this->update_backbone(px.sol.currents, c.conductivity);
     } else {
-      this->update_backbone(py.sol.currents);
+      this->update_backbone(py.sol.currents, c.conductivity);
     }
 
     unsigned max_pos = UINT_MAX;
@@ -74,8 +78,8 @@ void network::fracture(hooks& m) {
 }
 
 void network::get_cycle_edges_helper(std::set<unsigned>& cycle_edges,
-                                     std::set<unsigned>& seen_vertices, unsigned v_prev,
-                                     unsigned v_cur) const {
+    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]) {
@@ -103,7 +107,7 @@ 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 {
+    unsigned vPrev, unsigned vCur, unsigned vEnd) const {
   for (unsigned ei : G.dual_vertices[vCur].ne) {
     if (fuses[ei]) {
       if (!cycle_edges.contains(ei)) {
@@ -134,14 +138,15 @@ 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 {
+    unsigned v1) const {
   std::array<unsigned, 2> sig = {0, 0};
   this->find_cycle_helper(sig, cycle_edges, v0, v0, v1);
   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 {
+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]) {
       if (!cycle_edges.contains(ei)) {
@@ -173,8 +178,8 @@ bool network::get_cycle_helper(std::array<unsigned, 2>& sig, std::set<unsigned>&
   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::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);
@@ -202,7 +207,7 @@ std::set<unsigned> network::get_cluster_edges(unsigned v0) const {
 }
 
 void network::get_tie_flaps_helper(std::set<unsigned>& added_edges, unsigned v0,
-                                   unsigned vCur) const {
+    unsigned vCur) const {
   for (unsigned ei : G.vertices[vCur].ne) {
     if (!backbone[ei]) {
       if (!added_edges.contains(ei)) {
@@ -245,14 +250,25 @@ std::list<std::set<unsigned>> network::get_tie_flaps(unsigned v0) const {
   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();
+void network::update_backbone(const std::vector<double>& c, const std::array<double, 2>& cc) {
+  bool done_x = cc[0] < 1.0 / G.edges.size();
+  bool done_y = cc[1] < 1.0 / G.edges.size();
+
+  auto cycle_check_current = [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));
+  };
+  auto cycle_check_wrap_x = [](std::array<unsigned, 2> c) {
+    return (c[0] % 2 == 0 && c[1] % 2 == 1);
+  };
+  auto cycle_check_wrap_y = [](std::array<unsigned, 2> c) {
+    return (c[0] % 2 == 1 && c[1] % 2 == 0);
+  };
 
   // 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])) {
-      if (!seen_pairs.contains({G.dual_edges[i].v[0], G.dual_edges[i].v[1]})) {
+      if ((!seen_pairs_x.contains({G.dual_edges[i].v[0], G.dual_edges[i].v[1]}) && !done_x) || (!seen_pairs_y.contains({G.dual_edges[i].v[0], G.dual_edges[i].v[1]}) && !done_y)) {
         // 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.
@@ -261,7 +277,7 @@ void network::update_backbone(const std::vector<double>& c) {
         // 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]);
+          this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[i].v[1]);
         if (G.dual_edges[i].crossings[0])
           base_sig[0]++;
         if (G.dual_edges[i].crossings[1])
@@ -272,11 +288,11 @@ void network::update_backbone(const std::vector<double>& c) {
         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]);
+            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]);
+            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]};
+            new_sig_1[1] + new_sig_2[1]};
 
           if (G.dual_edges[i].crossings[0])
             new_sig[0]++;
@@ -292,13 +308,10 @@ void network::update_backbone(const std::vector<double>& c) {
 
         // 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 (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_current)) {
           // If none does, we remove it from the backbone!
-          seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
+          seen_pairs_x[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
+          seen_pairs_y[{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
@@ -334,14 +347,13 @@ 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));
-                        })) {
+        } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_x)) {
           // 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;
+          seen_pairs_x[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
+        } else  if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_y)) {
+          seen_pairs_y[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
         }
       }
     }
@@ -371,7 +383,7 @@ void network::update_backbone(const std::vector<double>& c) {
         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;
-          if (!seen_pairs.contains({G.vertices[v].nd[il], G.vertices[v].nd[ig]})) {
+          if ((!seen_pairs_x.contains({G.vertices[v].nd[il], G.vertices[v].nd[ig]}) && !done_x) || (!seen_pairs_y.contains({G.vertices[v].nd[il], G.vertices[v].nd[ig]}) && !done_y)) {
             bool any_intervening1 = false;
             bool any_intervening2 = false;
             unsigned ie1 = 0;
@@ -449,14 +461,10 @@ void network::update_backbone(const std::vector<double>& c) {
                 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);
-                              })) {
+              if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_current)) {
                 // one of our pairs qualifies for trimming!
-                seen_pairs[{d1, d2}] = true;
+                seen_pairs_x[{d1, d2}] = true;
+                seen_pairs_y[{d1, d2}] = true;
 
                 // We separate the flaps of the bowtie (there may be more than
                 // two!) into distinct sets.
@@ -490,11 +498,10 @@ 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;
+              } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_x)) {
+                seen_pairs_x[{d1, d2}] = true;
+              } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_y)) {
+                seen_pairs_y[{d1, d2}] = true;
               }
             }
           }
@@ -520,7 +527,7 @@ std::string network::write() {
   output += "\"fuses\"->{";
   for (unsigned i = 0; i < G.edges.size(); i++) {
     if (!fuses[i]) {
-      output += std::to_string(i) +",";
+      output += std::to_string(i) + ",";
     }
   }
   output.pop_back();
@@ -536,7 +543,8 @@ std::string network::write() {
     output += std::to_string(t) + ",";
   }
   output.pop_back();
-  output += "},\"conductivity\"->{" + std::to_string(c.conductivity[0]) + "," + std::to_string(c.conductivity[1]);
+  output += "},\"conductivity\"->{" + std::to_string(c.conductivity[0]) + "," +
+            std::to_string(c.conductivity[1]);
   output += "},\"currents\"->{";
   for (const double& t : c.currents) {
     output += std::to_string(t) + ",";
@@ -547,35 +555,23 @@ std::string network::write() {
   return output;
 };
 
+fuse_network::fuse_network(const graph& G, cholmod_common* c) : network(G, false, c), weight(1.0) {}
 
-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) {}
+fuse_network::fuse_network(const fuse_network& n) : network(n), weight(1.0) {}
 
 current_info fuse_network::get_current_info() {
   px.solve(fuses);
   py.solve(fuses);
 
   bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size();
-  bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size();
 
   current_info ctot;
   ctot.currents.resize(G.edges.size());
   ctot.conductivity = {px.sol.conductivity[0], py.sol.conductivity[1]};
 
-  if (done_x && !done_y) {
-    for (unsigned i = 0; i < G.edges.size(); i++) {
-      ctot.currents[i] = fabs(weight * py.sol.currents[i] / py.sol.conductivity[1]);
-    }
-  } else if (done_y && !done_x) {
-    for (unsigned i = 0; i < G.edges.size(); i++) {
-      ctot.currents[i] = fabs((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0]);
-    }
-  } else if (!done_x && !done_y) {
+  if (!done_x) {
     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]);
+      ctot.currents[i] = fabs(px.sol.currents[i] / px.sol.conductivity[0]);
     }
   }
 
@@ -583,7 +579,7 @@ current_info fuse_network::get_current_info() {
 }
 
 elastic_network::elastic_network(const graph& G, cholmod_common* c, double weight)
-    : network(G, c), weight(weight) {}
+    : network(G, true, c), weight(weight) {}
 
 elastic_network::elastic_network(const elastic_network& n) : network(n), weight(n.weight) {}
 
@@ -617,7 +613,7 @@ 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 graph& G, cholmod_common* c) : network(G, true, c) {}
 
 percolation_network::percolation_network(const percolation_network& n) : network(n) {}
 
diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp
index 2965915..8ad8235 100644
--- a/src/analysis_tools.cpp
+++ b/src/analysis_tools.cpp
@@ -30,10 +30,22 @@ std::pair<std::array<unsigned, 2>, std::set<unsigned>> find_minimal_crack(const
   }
 
   if (all_good) {
-    if (cycles[0].second.size() > cycles[1].second.size()) {
-      return cycles[1];
+    if (n.two_sides) {
+      if (cycles[0].second.size() > cycles[1].second.size()) {
+        return cycles[1];
+      } else {
+        return cycles[0];
+      }
     } else {
-      return cycles[0];
+      if (cycles[0].first[0] % 2 == 0) {
+        return cycles[0];
+      } else {
+        return cycles[1];
+      }
+    }
+  } else if (!n.two_sides) { 
+    if (cycles[!not_good].first[0] % 2 == 0) {
+      return cycles[!not_good];
     }
   }
 
@@ -63,7 +75,7 @@ std::pair<std::array<unsigned, 2>, std::set<unsigned>> find_minimal_crack(const
     pos++;
   }
 
-  if (cycles[!not_good].second.size() > new_cycle.size()) {
+  if (cycles[!not_good].second.size() > new_cycle.size() || !n.two_sides) {
     return {{sum_sig_0, sum_sig_1}, new_cycle};
   } else {
     return cycles[!not_good];
diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp
index e8de8c1..fd1ffe2 100644
--- a/src/animate_fracture.cpp
+++ b/src/animate_fracture.cpp
@@ -69,6 +69,33 @@ int main(int argc, char* argv[]) {
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
 
+  if (w == 0.0 || w == 1.0) {
+  if (use_aN) {
+    animate meas(sqrt(2*n *a), sqrt( 2*n / a), width, argc, argv);
+
+    for (unsigned trial = 0; trial < N; trial++) {
+      graph G(n, a, rng);
+      fuse_network net(G, &c);
+      net.set_thresholds(beta, rng);
+      net.fracture(meas);
+    }
+  } else {
+    animate meas(Lx, Ly, width, argc, argv);
+    const graph G(Lx, Ly);
+    const fuse_network plain_net(G, &c);
+
+    for (unsigned trial = 0; trial < N; trial++) {
+      fuse_network net = plain_net;
+      net.set_thresholds(beta, rng);
+      try {
+      net.fracture(meas);
+      } catch (std::exception &e) {
+        std::cout << e.what() << std::endl;
+        getchar();
+      }
+    }
+  }
+  } else {
   if (use_aN) {
     animate meas(sqrt(2*n *a), sqrt( 2*n / a), width, argc, argv);
 
@@ -89,6 +116,7 @@ int main(int argc, char* argv[]) {
       net.fracture(meas);
     }
   }
+  }
 
   CHOL_F(finish)(&c);
 
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 96cdac3..22101c8 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -76,6 +76,49 @@ int main(int argc, char* argv[]) {
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
 
+  if (w == 0 || w == 1) {
+  if (use_aN) {
+    ma meas(n, a, beta, w);
+
+    for (unsigned trial = 0; trial < N; trial++) {
+      while (true) {
+        try {
+          graph G(n, a, rng);
+          fuse_network n(G, &c);
+          n.set_thresholds(beta, rng);
+          n.fracture(meas);
+          break;
+        } catch (std::exception &e) {
+          std::cout << e.what() << '\n';
+        }
+      }
+
+      if (quit.load())
+        break;
+    }
+  } else {
+    ma meas(Lx, Ly, beta, w);
+
+    const graph G(Lx, Ly);
+    const fuse_network n(G, &c);
+
+    for (unsigned trial = 0; trial < N; trial++) {
+      while (true) {
+        try {
+          fuse_network net = n;
+          net.set_thresholds(beta, rng);
+          net.fracture(meas);
+          break;
+        } catch (std::exception &e) {
+          std::cout << e.what() << '\n';
+        }
+      }
+
+      if (quit.load())
+        break;
+    }
+  }
+  } else {
   if (use_aN) {
     ma meas(n, a, beta, w);
 
@@ -117,6 +160,7 @@ int main(int argc, char* argv[]) {
         break;
     }
   }
+  }
 
   CHOL_F(finish)(&c);