changes

4 files changed, 188 insertions, 172 deletions

diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 9addfdd..a27ca39 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -1,7 +1,7 @@
 
 #include "network.hpp"
 
-network::network(const graph& G, cholmod_common *c) : G(G), c(c), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) {
+network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) {
   b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c);
   for (unsigned int i = 0; i < G.edges.size(); i++) {
     double v0y = G.vertices[G.edges[i][0]].r.y;
@@ -36,8 +36,8 @@ network::network(const graph& G, cholmod_common *c) : G(G), c(c), fuses(G.edges.
     unsigned int s0 = v0 < v1 ? v0 : v1;
     unsigned int s1 = v0 < v1 ? v1 : v0;
 
-    ((CHOL_INT *)t->i)[terms] = v1;
-    ((CHOL_INT *)t->j)[terms] = v0;
+    ((CHOL_INT *)t->i)[terms] = s1;
+    ((CHOL_INT *)t->j)[terms] = s0;
     ((double *)t->x)[terms] = -1.0;
 
     ((double *)t->x)[v0]++;
@@ -74,7 +74,7 @@ network::network(const graph& G, cholmod_common *c) : G(G), c(c), fuses(G.edges.
   CHOL_F(free_triplet)(&t, c);
 }
 
-network::network(const network &other) : G(other.G), thresholds(other.thresholds), fuses(other.fuses), c(other.c) {
+network::network(const network &other) : c(other.c), G(other.G), fuses(other.fuses), thresholds(other.thresholds) {
   b = CHOL_F(copy_dense)(other.b, c);
   factor = CHOL_F(copy_factor)(other.factor, c);
   voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c);
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 57a817b..d93783c 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -3,5 +3,5 @@ add_library(measurements measurements.cpp)
 target_link_libraries(measurements frac)
 
 add_executable(fracture fracture.cpp)
-target_link_libraries(fracture frac measurements fftw3 cholmod)
+target_link_libraries(fracture frac measurements fftw3 cholmod profiler)
 
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 578c70a..0bd72c4 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -11,6 +11,119 @@ bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>)
   return true;
 }
 
+std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) {
+  Graph Gtmp(n.G.vertices.size());
+  std::list<unsigned int> removed_edges;
+
+  class add_tree_edges : public boost::default_dfs_visitor {
+    public:
+      Graph& G;
+      std::list<unsigned int>& E;
+
+      add_tree_edges(Graph& G, std::list<unsigned int>& 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:
+      unsigned int end;
+      std::list<unsigned int>& E;
+      struct done{};
+
+      find_cycle(std::list<unsigned int>& E, unsigned int 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 int>> cycles;
+
+  for (auto edge : removed_edges) {
+    std::list<unsigned int> cycle = {edge};
+    find_cycle vis(cycle, n.G.dual_edges[edge][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][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) {
+    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_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 int> new_cycle;
+        unsigned int pos = 0;
+        for (uint8_t included : new_bool_cycle) {
+          if (included) {
+            new_cycle.push_back(pos);
+          }
+          pos++;
+        }
+        cycles.push_back(new_cycle);
+      }
+    }
+
+    // find the cycle representing the crack by counting boundary crossings
+    for (auto cycle : cycles) {
+      std::array<unsigned int, 2> crossing_count{0,0};
+
+      for (auto edge : cycle) {
+        double dx = fabs(n.G.vertices[n.G.dual_edges[edge][0]].r.x - n.G.vertices[n.G.dual_edges[edge][1]].r.x);
+        if (dx > 0.5) {
+          crossing_count[0]++;
+        }
+        double dy = fabs(n.G.vertices[n.G.dual_edges[edge][0]].r.y - n.G.vertices[n.G.dual_edges[edge][1]].r.y);
+        if (dy > 0.5) {
+          crossing_count[1]++;
+        }
+      }
+
+      if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) {
+        return cycle;
+      }
+    }
+  } else {
+    return cycles.front();
+  }
+
+  exit(5);
+}
+
 void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned int N, unsigned int L, double beta) {
   std::string filename = "fracture_" + std::to_string(L) + "_" + std::to_string(beta) + "_" + id + ".dat";
   std::ifstream file(filename);
@@ -97,12 +210,8 @@ void correlation(unsigned int L, std::vector<T>& data, const std::vector<fftw_co
   }
 }
 
-template <class T, class U>
-void autocorrelation(unsigned int L, std::vector<T>& out_data, const std::vector<U>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
-  for (unsigned int i = 0; i < pow(L, 2); i++) {
-    fftw_forward_in[i] = (double)data[i];
-  }
-
+template <class T>
+void autocorrelation(unsigned int L, std::vector<T>& out_data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
   fftw_execute(forward_plan);
 
   for (unsigned int i = 0; i < L * (L / 2 + 1); i++) {
@@ -124,8 +233,11 @@ ma::ma(unsigned int N, unsigned int L, double beta) : L(L), G(2 * pow(L / 2, 2))
   sb(log2(L) + 1, 0),
   Ccc(pow(L / 2 + 1, 2), 0),
   Css(pow(L / 2 + 1, 2), 0),
+  Cmm(pow(L / 2 + 1, 2), 0),
+  Caa(pow(L / 2 + 1, 2), 0),
   Cdd(pow(L / 2 + 1, 2), 0),
-  Caa(pow(L / 2 + 1, 2), 0)
+  Cll(pow(L / 2 + 1, 2), 0),
+  Cee(pow(L / 2 + 1, 2), 0)
 {
   Nc = 0;
   Na = 0;
@@ -159,34 +271,36 @@ ma::~ma() {
 
   update_field_file("Ccc", Ccc, Nc, L, beta);
   update_field_file("Css", Css, N, L, beta);
+  update_field_file("Cmm", Cmm, N, L, beta);
   update_field_file("Cdd", Cdd, N, L, beta);
   update_field_file("Caa", Caa, Na, L, beta);
-
+  update_field_file("Cll", Cll, N, L, beta);
+  update_field_file("Cee", Cee, N, L, beta);
 }
 
 void ma::pre_fracture(const network &) {
   lv = 0;
-  avalanche.clear();
+  avalanches = {{}};
   boost::remove_edge_if(trivial, G);
 }
 
 void ma::bond_broken(const network& net, const std::pair<double, std::vector<double>>& cur, unsigned int i) {
   if (cur.first / fabs(cur.second[i]) * net.thresholds[i] > lv) {
-    sa[avalanche.size()]++;
+    sa[avalanches.back().size()]++;
     Na++;
 
-    std::vector<bool>avalanche_damage(pow(L, 2), 0);
+    memset(fftw_forward_in, 0.0, net.G.edges.size());
 
-    for (auto it = avalanche.begin(); it != avalanche.end(); it++) {
-      avalanche_damage[*it] = 1;
+    for (auto e : avalanches.back()) {
+      fftw_forward_in[e] = 1.0;
     }
 
-    autocorrelation(L, Caa, avalanche_damage, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+    autocorrelation(L, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
     lv = cur.first / fabs(cur.second[i]) * net.thresholds[i];
-    avalanche.clear();
+    avalanches.push_back({i});
   } else {
-    avalanche.push_back(i);
+    avalanches.back().push_back(i);
   }
 
   boost::add_edge(net.G.dual_edges[i][0], net.G.dual_edges[i][1], {i}, G);
@@ -196,135 +310,9 @@ void ma::post_fracture(network &n) {
   std::vector<unsigned int> component(boost::num_vertices(G));
   unsigned int num = boost::connected_components(G, &component[0]);
 
-  std::vector<unsigned int> comp_sizes(num, 0);
-
-  for (unsigned int c : component) {
-    comp_sizes[c]++;
-  }
+  std::list<unsigned int> crack = find_minimal_crack(G, n);
 
-  unsigned int max_i = 0;
-  unsigned int max_size = 0;
-
-  for (unsigned int i = 0; i < num; i++) {
-    if (comp_sizes[i] > max_size) {
-      max_i = i;
-      max_size = comp_sizes[i];
-    }
-  }
-
-  unsigned int vertex_in_largest = 0;
-
-  while (component[vertex_in_largest] != max_i) {
-    vertex_in_largest++;
-  }
-
-  Graph Gtmp(2 * pow(L / 2, 2));
-  std::list<unsigned int> removed_edges;
-  class add_tree_edges : public boost::default_dfs_visitor {
-    public:
-      Graph& G;
-      std::list<unsigned int>& E;
-
-      add_tree_edges(Graph& G, std::list<unsigned int>& 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:
-      unsigned int end;
-      std::list<unsigned int>& E;
-      struct done{};
-
-      find_cycle(std::list<unsigned int>& E, unsigned int 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 int>> cycles;
-
-  for (auto edge : removed_edges) {
-    std::list<unsigned int> cycle = {edge};
-    find_cycle vis(cycle, n.G.dual_edges[edge][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][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);
-    }
-  }
-
-  std::list<unsigned int> crack;
-
-  if (cycles.size() > 1) {
-    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_cycles.push_back(bool_cycle);
-    }
-
-    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 int> new_cycle;
-        unsigned int pos = 0;
-        for (uint8_t included : new_bool_cycle) {
-          if (included) {
-            new_cycle.push_back(pos);
-          }
-          pos++;
-        }
-        cycles.push_back(new_cycle);
-      }
-    }
-
-    for (auto cycle : cycles) {
-      std::array<unsigned int, 2> crossing_count{0,0};
-
-      for (auto edge : cycle) {
-        double dx = fabs(n.G.vertices[n.G.edges[edge][0]].r.x - n.G.vertices[n.G.edges[edge][1]].r.x);
-        if (dx > 0.5) {
-          crossing_count[0]++;
-        }
-        double dy = fabs(n.G.vertices[n.G.edges[edge][0]].r.y - n.G.vertices[n.G.edges[edge][1]].r.y);
-        if (dy > 0.5) {
-          crossing_count[1]++;
-        }
-      }
-
-      if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) {
-        crack = cycle;
-        break;
-      }
-    }
-  } else {
-    crack = cycles.front();
-  }
+  unsigned int crack_component = component[n.G.dual_edges[crack.front()][0]];
 
   for (unsigned int be = 0; be < log2(L); be++) {
     unsigned int bin = pow(2, be);
@@ -333,7 +321,7 @@ void ma::post_fracture(network &n) {
       bool in_bin = false;
       for (unsigned int j = 0; j < pow(bin, 2); j++) {
         unsigned int edge = L * (bin * ((i * bin) / L) + j / bin) + (i * bin) % L + j % bin;
-        if (!n.fuses[edge] && max_i == component[n.G.dual_edges[edge][0]]) {
+        if (!n.fuses[edge] && crack_component == component[n.G.dual_edges[edge][0]]) {
           in_bin = true;
           break;
         }
@@ -347,62 +335,86 @@ void ma::post_fracture(network &n) {
 
   sb[log2(L)]++;
 
+  for (unsigned int i = 0; i < n.G.edges.size(); i++) {
+    if (n.fuses[i] && component[n.G.dual_edges[i][0]] == crack_component)  {
+      fftw_forward_in[i] = 1.0;
+    } else {
+      fftw_forward_in[i] = 0.0;
+    }
+  }
+
+  autocorrelation(L, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+
   // crack surface correlations
-  std::vector<double> crack_damage(pow(L, 2), 0.0);
+  memset(fftw_forward_in, 0.0, n.G.edges.size());
 
   for (auto edge : crack) {
-    crack_damage[edge] = 1.0;
+    fftw_forward_in[edge] = 1.0;
+  }
+
+  autocorrelation(L, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+
+  std::function<bool(unsigned int)> inCrack = [&](unsigned int i) -> bool {
+    return (bool)fftw_forward_in[i];
+  };
+
+  for (auto avalanche : avalanches) {
+    if (avalanche.end() != std::find_if(avalanche.begin(), avalanche.end(), inCrack)) {
+      for (auto edge : avalanche) {
+        fftw_forward_in[edge] = 1.0;
+      }
+    }
   }
 
-  std::vector<fftw_complex> t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out);
-  correlation(L, Css, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(L, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
-  std::vector<bool> last_crack(pow(L, 2), 0);
+  memset(fftw_forward_in, 0.0, n.G.edges.size());
 
   // rewind the last avalanche
-  for (auto e = avalanche.rbegin(); e != avalanche.rend(); e++) {
-    boost::remove_edge(n.G.dual_edges[*e][0], n.G.dual_edges[*e][1], G);
-    n.add_edge(*e);
-    last_crack[*e] = 1;
+  for (auto e : avalanches.back()) {
+    boost::remove_edge(n.G.dual_edges[e][0], n.G.dual_edges[e][1], G);
+    n.add_edge(e);
+    fftw_forward_in[e] = 1;
   }
 
+  autocorrelation(L, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
   // cluster size distribution and cluster-cluster correlation
   num = boost::connected_components(G, &component[0]);
 
   for (unsigned int i = 0; i < num; i++) {
-    double size = 0;
-    std::fill(crack_damage.begin(), crack_damage.end(), 0.0);
+    unsigned int size = 0;
 
-    for (unsigned int j = 0; j < pow(L, 2); j++) {
-      if (component[n.G.edges[j][0]] == i && n.fuses[j]) {
+    for (unsigned int j = 0; j < n.G.edges.size(); j++) {
+      if (component[n.G.dual_edges[j][0]] == i && n.fuses[j]) {
         size++;
-        crack_damage[j] = 1.0;
+        fftw_forward_in[j] = 1.0;
+      } else{
+        fftw_forward_in[j] = 0.0;
       }
     }
 
     if (size > 0) {
       sc[size - 1]++;
-      t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out);
-      correlation(L, Ccc, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+      autocorrelation(L, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
       Nc++;
     }
   }
 
-  // damage-damage correlations
-  std::vector<fftw_complex> t_damage = data_transform(L, n.fuses, forward_plan, fftw_forward_in, fftw_forward_out);
-  correlation(L, Cdd, t_damage, t_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out);
-
-
   // damage size distribution
   unsigned int total_broken = 0;
 
   for (unsigned int i = 0; i < n.G.edges.size(); i++) {
     if (n.fuses[i]) { 
       total_broken++;
+      fftw_forward_in[i] = 1.0;
+    } else {
+      fftw_forward_in[i] = 0.0;
     }
   }
 
+  autocorrelation(L, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+
   sd[total_broken]++;
 
 }
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 4c05d93..621c425 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -2,6 +2,7 @@
 #include <vector>
 #include <algorithm>
 #include <cmath>
+#include <cstring>
 #include <list>
 #include <fstream>
 #include <string>
@@ -49,8 +50,11 @@ class ma : public hooks {
     std::vector<uint64_t> sb; // bin size distribution
     std::vector<uint64_t> Ccc; // cluster-cluster correlations
     std::vector<uint64_t> Css; // surface-surface correlations
+    std::vector<uint64_t> Cmm; // surface-surface correlations
     std::vector<uint64_t> Caa; // avalanche-avalanche correlations
     std::vector<uint64_t> Cdd; // damage-damage distribution
+    std::vector<uint64_t> Cll; // damage-damage distribution
+    std::vector<uint64_t> Cee; // damage-damage distribution
     uint64_t Nc;
     uint64_t Na;
 
@@ -61,7 +65,7 @@ class ma : public hooks {
     double beta;
 
 
-    std::list<unsigned int> avalanche;
+    std::list<std::list<unsigned int>> avalanches;
 
     ma(unsigned int N, unsigned int L, double beta);
     ~ma();