more fixes, and updating the interface

5 files changed, 54 insertions, 76 deletions

diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index dcc5e0e..ed23d45 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -90,6 +90,9 @@ void network::fracture(hooks& m, bool one_axis) {
 }
 
 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 get_cycle_edges : public boost::default_dfs_visitor {
     public:
       std::set<unsigned> tE;
@@ -191,7 +194,7 @@ void network::update_backbone(const std::vector<double>& c) {
           other_cycles.push_back(tcycle);
         }
 
-        if (std::any_of(other_cycles.begin(), other_cycles.end(), [](std::array<unsigned, 2> c){return c[0] % 2 == 0 && c[1] % 2 == 0;})) {
+        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);
@@ -304,6 +307,8 @@ void network::update_backbone(const std::vector<double>& c) {
           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)) {
@@ -320,33 +325,24 @@ void network::update_backbone(const std::vector<double>& c) {
           if (any_intervening2 && !any_intervening1) {
             for (unsigned k = il + 1; k <= ig; k++) {
               if (!backbone[G.vertices[v].ne[k]]) {
-                any_intervening_e_1 = true;
-                break;
+                ie1++;
               }
             }
           }
           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()]]) {
-                any_intervening_e_2 = true;
-                break;
+                ie2++;
               }
             }
           }
-          if ((any_intervening1 && any_intervening2) || (any_intervening1 && any_intervening_e_2) || (any_intervening2 && any_intervening_e_1)) {
+          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];
-            break;
-          }
-        }
-        }
-      }
-    }
 
-    if (found_pair) {
       std::list<unsigned> cedges = {};
       get_cycle_edges vis(cedges);
       std::vector<boost::default_color_type> cm1(G.dual_vertices.size());
@@ -390,7 +386,7 @@ void network::update_backbone(const std::vector<double>& c) {
           other_cycles.push_back(tcycle);
         }
 
-        if (std::any_of(other_cycles.begin(), other_cycles.end(), [](std::array<unsigned, 2> c){return c[0] % 2 == 0 && c[1] % 2 == 0;})) {
+        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 = {};
@@ -456,6 +452,10 @@ void network::update_backbone(const std::vector<double>& c) {
       }
   }
   }
+        }
+      }
+    }
+  }
 }
 
 void network::break_edge(unsigned e, bool unbreak) {
diff --git a/src/animate.cpp b/src/animate.cpp
index 532bd70..b1031dc 100644
--- a/src/animate.cpp
+++ b/src/animate.cpp
@@ -89,12 +89,14 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
     */
 
    bool weird = false;
+   unsigned nw = 0;
 
     for (unsigned j = 0; j < n.G.edges.size(); j++) {
       bool draw = false;
       if (cur.currents[j] < 1e-9 && !n.backbone[j]) {
     glColor3f(1.0f, 0.0f, 0.0f); // Blue
     weird = true;
+    nw++;
     draw = true;
       } else if (n.backbone[j] && !n.fuses[j] && j != i) {
     glColor3f(0.8f, 0.8f, 0.8f); // Blue
@@ -102,6 +104,9 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
       } else if (!n.fuses[j]) {
     glColor3f(0.0f, 0.0f, 0.0f); // Blue
     draw = true;
+      } else if (i == j) {
+    glColor3f(0.0f, 1.0f, 0.0f); // Blue
+    draw = true;
       }
 
       if (draw) {
@@ -130,7 +135,7 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
       }
  glEnd();
  glFlush();
- if (weird) {std::cout << "\n"; getchar();}
+ if (nw > 2) {std::cout << "\n"; getchar();}
 }
 
 void animate::post_fracture(network &n) {
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 483a3d2..5ba60e2 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -38,9 +38,10 @@ int main(int argc, char* argv[]) {
   unsigned n = 128;
   double a = 1.0;
   bool use_aN = false;
+  bool one_side = false;
   double w = 0.5;
 
-  while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:w:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:w:o")) != -1) {
     switch (opt) {
       case 'N':
         N = (unsigned)atof(optarg);
@@ -65,6 +66,9 @@ int main(int argc, char* argv[]) {
       case 'w':
         w = atof(optarg);
         break;
+      case 'o':
+        one_side = true;
+        break;
       default:
         exit(1);
     }
@@ -77,10 +81,7 @@ int main(int argc, char* argv[]) {
   std::mt19937 rng{seeds};
 
   if (use_aN) {
-    unsigned Mx = (unsigned)4*sqrt(2*n*a);
-    unsigned My = (unsigned)4*sqrt(2*n/a);
-
-    ma meas(n, a, beta);
+    ma meas(n, a, beta, w, one_side);
 
     for (unsigned trial = 0; trial < N; trial++) {
       while (true) {
@@ -88,7 +89,7 @@ int main(int argc, char* argv[]) {
           graph G(n, a, rng);
           elastic_network fuse_network(G, &c, w);
           fuse_network.set_thresholds(beta, rng);
-          fuse_network.fracture(meas, false);
+          fuse_network.fracture(meas, one_side);
           break;
         } catch (std::exception &e) {
           std::cout << e.what() << '\n';
@@ -99,11 +100,7 @@ int main(int argc, char* argv[]) {
         break;
     }
   } else {
-    // fourier transforms of powers of two are faster
-    unsigned Mx = pow(2, ceil(log2(4*Lx)));
-    unsigned My = pow(2, ceil(log2(4*Ly)));
-
-    ma meas(Lx, Ly, beta);
+    ma meas(Lx, Ly, beta, w, one_side);
 
     for (unsigned trial = 0; trial < N; trial++) {
       while (true) {
@@ -111,7 +108,7 @@ int main(int argc, char* argv[]) {
           graph G(Lx, Ly);
           elastic_network fuse_network(G, &c, w);
           fuse_network.set_thresholds(beta, rng);
-          fuse_network.fracture(meas, false);
+          fuse_network.fracture(meas, one_side);
           break;
         } catch (std::exception &e) {
           std::cout << e.what() << '\n';
diff --git a/src/measurements.cpp b/src/measurements.cpp
index c58fd84..96b4ccb 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -142,45 +142,55 @@ unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, u
   return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly);
 }
 
-ma::ma(unsigned n, double a, double beta) :
+ma::ma(unsigned n, double a, double beta, double weight, bool one) :
   G(2 * n),
   sn(2 * n),
   ss(2 * n),
   sm(2 * n),
   sl(2 * n),
-  sb(n),
+  sb(n + 1),
   sd(3 * n),
-  sc(2 * n),
   sa(3 * n),
   sA(3 * n),
   si(10000),
   sI(10000)
 {
   if (beta != 0.0) {
-    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_";
+    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_" + std::to_string(weight) + "_";
   } else {
-    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_INF_";
+    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_INF_" + std::to_string(weight) + "_";
+  }
+
+  if (one) {
+    model_string = model_string + "o_";
+  } else {
+    model_string = model_string + "t_";
   }
 }
 
-ma::ma(unsigned Lx, unsigned Ly, double beta) :
+ma::ma(unsigned Lx, unsigned Ly, double beta, double weight, bool one) :
   G(Lx * Ly / 2),
   sn(Lx * Ly / 2),
   ss(Lx * Ly / 2),
   sm(Lx * Ly / 2),
   sl(Lx * Ly / 2),
-  sb(Lx * Ly / 2),
+  sb(Lx * Ly / 2 + 1),
   sd(Lx * Ly),
-  sc(Lx * Ly / 2),
   sa(Lx * Ly),
   sA(Lx * Ly),
   si(10000),
   sI(10000)
 {
   if (beta != 0.0) {
-    model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_";
+    model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + std::to_string(weight) + "_";
   } else {
-    model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_INF_";
+    model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_INF_" + std::to_string(weight) + "_";
+  }
+
+  if (one) {
+    model_string = model_string + "o_";
+  } else {
+    model_string = model_string + "t_";
   }
 }
 
@@ -191,7 +201,6 @@ ma::~ma() {
   update_distribution_file("sl", sl, model_string);
   update_distribution_file("sb", sb, model_string);
   update_distribution_file("sd", sd, model_string);
-  update_distribution_file("sc", sc, model_string);
   update_distribution_file("sa", sa, model_string);
   update_distribution_file("sA", sA, model_string);
   update_distribution_file("si", si, model_string);
@@ -260,43 +269,10 @@ void ma::post_fracture(network &n) {
     sn[components[i].size() - 1]++;
   }
 
-  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.break_edge(e, true);
-    }
-
-    auto cracks = find_minimal_crack(G, n);
-
-    if (cracks.size() == 0) {
-      break;
-    }
-
-    av_it++;
-  }
-
-  num = boost::connected_components(G, &component[0]);
-
-  std::vector<std::list<graph::coordinate>> new_components(num);
-
-  for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
-    new_components[component[i]].push_back(n.G.dual_vertices[i].r);
-  }
-
-  for (unsigned i = 0; i < num; i++) {
-    sc[new_components[i].size() - 1]++;
-  }
-
-  /*
-  current_info ct = n.get_current_info();
-
-
   std::vector<bool> vertex_in(n.G.vertices.size());
 
   for (unsigned i = 0; i < n.G.edges.size(); i++) {
-    if (ct.currents[i] > CURRENT_CUTOFF) {
+    if (!n.backbone[i]) {
       vertex_in[n.G.edges[i].v[0]] = true;
       vertex_in[n.G.edges[i].v[1]] = true;
     }
@@ -308,9 +284,9 @@ void ma::post_fracture(network &n) {
     if (vertex_in[i]) bb_size++;
   }
 
-  sb[bb_size - 1]++;
-  */
+  sb[bb_size]++;
 
+  auto av_it = avalanches.rbegin();
   av_it++;
 
   while (av_it != avalanches.rend()) {
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 274a550..c211974 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -44,8 +44,8 @@ class ma : public hooks {
     std::list<unsigned> last_avalanche;
     std::string model_string;
 
-    ma(unsigned Lx, unsigned Ly, double beta);
-    ma(unsigned n, double a, double beta);
+    ma(unsigned Lx, unsigned Ly, double beta, double w, bool o);
+    ma(unsigned n, double a, double beta, double w, bool o);
     ~ma();
 
     void pre_fracture(const network &) override;