many changes to measurements, right now halfway through implementing a better, self-consistent record of the crack surface

2 files changed, 270 insertions, 62 deletions

diff --git a/src/measurements.cpp b/src/measurements.cpp
index 7d5d539..1159fd6 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -1,11 +1,17 @@
 
 #include "measurements.hpp"
+#include <iostream>
+
+template <class T>
+bool is_shorter(std::list<T> l1, std::list<T> l2) {
+  return l1.size() < l2.size();
+}
 
 bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>) {
   return true;
 }
 
-void update_distribution_data(std::string id, const std::vector<unsigned int>& data, unsigned int N, unsigned int L, double beta) {
+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);
 
@@ -33,58 +39,43 @@ void update_distribution_data(std::string id, const std::vector<unsigned int>& d
   file_out.close();
 }
 
-void update_field_data(std::string id, double tot, const std::vector<double>& data, unsigned int L, double beta) {
+template <class T>
+void update_field_file(std::string id, const std::vector<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);
 
   uint64_t N_old = 0;
-  uint64_t tot_old = 0;
-  uint64_t tot_2_old = 0;
   std::vector<uint64_t> data_old(data.size(), 0);
-  std::vector<uint64_t> data_old_2(data.size(), 0);
 
   if (file.is_open()) {
     file >> N_old;
-    file >> tot_old;
-    file >> tot_2_old;
     for (unsigned int i = 0; i < data.size(); i++) {
-      uint64_t num;
-      file >> num;
-      data_old[i] = num;
+      file >> data_old[i];
     }
-    for (unsigned int i = 0; i < data.size(); i++) {
-      uint64_t num;
-      file >> num;
-      data_old_2[i] = num;
-    }
-
     file.close();
   }
 
   std::ofstream file_out(filename);
 
-  file_out <<std::fixed<< N_old + 1 << " " << (uint64_t)(tot_old + tot) << " " << (uint64_t)(tot_2_old + pow(tot, 2))<< "\n";
-  for (unsigned int i = 0; i < data.size(); i++) {
-    file_out << (uint64_t)(data_old[i] + data[i]) << " ";
-  }
-  file_out <<std::fixed<< "\n";
+  file_out <<std::fixed<< N_old + N <<  "\n";
   for (unsigned int i = 0; i < data.size(); i++) {
-    file_out << (uint64_t)(data_old_2[i] + pow(data[i], 2)) << " ";
+    file_out << data_old[i] + data[i] << " ";
   }
 
   file_out.close();
 }
 
-std::vector<fftw_complex> data_transform(unsigned int L, const std::vector<double>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) {
+template <class T>
+std::vector<fftw_complex> data_transform(unsigned int L, const std::vector<T>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) {
   for (unsigned int i = 0; i < pow(L, 2); i++) {
-    fftw_forward_in[i] = data[i];
+    fftw_forward_in[i] = (double)data[i];
   }
 
   fftw_execute(forward_plan);
 
-  std::vector<fftw_complex> output(pow(L, 2));
+  std::vector<fftw_complex> output(L * (L / 2 + 1));
 
-  for (unsigned int i = 0; i < pow(L, 2); i++) {
+  for (unsigned int i = 0; i < L * (L / 2 + 1); i++) {
     output[i][0] = fftw_forward_out[i][0];
     output[i][1] = fftw_forward_out[i][1];
   }
@@ -92,26 +83,52 @@ std::vector<fftw_complex> data_transform(unsigned int L, const std::vector<doubl
   return output;
 }
 
-std::vector<double> correlation(unsigned int L, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
-  for (unsigned int i = 0; i < pow(L, 2); i++) {
+template <class T>
+void correlation(unsigned int L, std::vector<T>& data, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
+  for (unsigned int i = 0; i < L * (L / 2 + 1); i++) {
     fftw_reverse_in[i][0] = tx1[i][0] * tx2[i][0] + tx1[i][1] * tx2[i][1];
     fftw_reverse_in[i][1] = tx1[i][0] * tx2[i][1] - tx1[i][1] * tx2[i][0];
   }
 
   fftw_execute(reverse_plan);
 
-  std::vector<double> output(pow(L / 2 + 1, 2));
-
   for (unsigned int i = 0; i < pow(L / 2 + 1, 2); i++) {
-    output[i] = fftw_reverse_out[L * (i / (L / 2 + 1)) + i % (L / 2 + 1)] / pow(L, 2);
+    data[i] += (T)(fftw_reverse_out[L * (i / (L / 2 + 1)) + i % (L / 2 + 1)] / pow(L, 2));
   }
+}
 
-  return output;
+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];
+  }
+
+  fftw_execute(forward_plan);
+
+  for (unsigned int i = 0; i < L * (L / 2 + 1); i++) {
+    fftw_reverse_in[i][0] = pow(fftw_forward_out[i][0], 2) + pow(fftw_forward_out[i][1], 2);
+    fftw_reverse_in[i][1] = 0.0;
+  }
+
+  fftw_execute(reverse_plan);
+
+  for (unsigned int i = 0; i < pow(L / 2 + 1, 2); i++) {
+    out_data[i] += (T)(fftw_reverse_out[L * (i / (L / 2 + 1)) + i % (L / 2 + 1)] / pow(L, 2));
+  }
 }
 
-ma::ma(unsigned int N, unsigned int L, double beta) : L(L), G(2 * pow(L / 2, 2)), bin_counts(log2(L) + 1, 0), N(N), beta(beta) {
-  ad.resize(pow(L, 2), 0);
-  cd.resize(pow(L, 2), 0);
+ma::ma(unsigned int N, unsigned int L, double beta) : L(L), G(2 * pow(L / 2, 2)), N(N), beta(beta),
+  sc(pow(L, 2), 0),
+  sa(pow(L, 2), 0),
+  sd(pow(L, 2), 0),
+  sb(log2(L) + 1, 0),
+  Ccc(pow(L / 2 + 1, 2), 0),
+  Css(pow(L / 2 + 1, 2), 0),
+  Cdd(pow(L / 2 + 1, 2), 0),
+  Caa(pow(L / 2 + 1, 2), 0)
+{
+  Nc = 0;
+  Na = 0;
 
   // FFTW setup for correlation functions
   fftw_set_timelimit(1);
@@ -135,31 +152,44 @@ ma::~ma() {
   fftw_destroy_plan(reverse_plan);
   fftw_cleanup();
 
-  update_distribution_data("na", ad, N, L, beta);
-  update_distribution_data("nc", cd, N, L, beta);
-  update_distribution_data("bc", bin_counts, N, L, beta);
+  update_distribution_file("sa", sa, N, L, beta);
+  update_distribution_file("sc", sc, N, L, beta);
+  update_distribution_file("sd", sd, N, L, beta);
+  update_distribution_file("sb", sb, N, L, beta);
+
+  update_field_file("Ccc", Ccc, Nc, L, beta);
+  update_field_file("Css", Css, N, L, beta);
+  update_field_file("Cdd", Cdd, N, L, beta);
+  update_field_file("Caa", Caa, Na, L, beta);
 
 }
 
 void ma::pre_fracture(const network &) {
   lv = 0;
-  as = 0;
-  avalanches.clear();
+  avalanche.clear();
   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) {
-    ad[as]++;
-    as = 0;
+    sa[avalanche.size()]++;
+    Na++;
+
+    std::vector<bool>avalanche_damage(pow(L, 2), 0);
+
+    for (auto it = avalanche.begin(); it != avalanche.end(); it++) {
+      avalanche_damage[*it] = 1;
+    }
+
+    autocorrelation(L, Caa, avalanche_damage, 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];
-    avalanches.push_back({});
+    avalanche.clear();
   } else {
-    as++;
-    avalanches.back().push_back(i);
+    avalanche.push_back(i);
   }
 
-  boost::add_edge(net.G.dual_edges[i][0], net.G.dual_edges[i][1], G);
+  boost::add_edge(net.G.dual_edges[i][0], net.G.dual_edges[i][1], {i}, G);
 }
 
 void ma::post_fracture(network &n) {
@@ -182,6 +212,144 @@ void ma::post_fracture(network &n) {
     }
   }
 
+  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>& V;
+      std::list<unsigned int>& E;
+      struct done{};
+
+      find_cycle(std::list<unsigned int>& V, std::list<unsigned int>& E, unsigned int end) : V(V), E(E), end(end) {}
+
+      void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) {
+        V.push_back(v);
+        if (v == end) {
+          throw done{};
+        }
+      }
+
+      void finish_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) {
+        V.remove(v);
+      }
+
+      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));
+      }
+  };
+
+  if (removed_edges.size() == 0) {
+    std::pair<double, std::vector<double>> currents = n.currents();
+    std::cout << currents.first << "\n";
+    for (unsigned int i = 0; i < L; i++) {
+      for (unsigned int j = 0; j < L; j++) {
+        std::cout << n.fuses[i * L + j] << " ";
+      }
+      std::cout << "\n";
+    }
+  }
+
+  std::list<std::list<unsigned int>> cycles;
+
+  for (auto edge : removed_edges) {
+    std::list<unsigned int> cycle = {edge};
+    std::list<unsigned int> vcycle = {};
+    find_cycle vis(vcycle, 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();
+  }
+
+  for (auto edge : crack) {
+    std::cout << edge << " ";
+  }
+  std::cout << "\n";
+//  std::cout << max_size << " " << boost::num_edges(Gtmp) << " " << removed_edges.size() << "\n";
+
   for (unsigned int be = 0; be < log2(L); be++) {
     unsigned int bin = pow(2, be);
 
@@ -196,13 +364,14 @@ void ma::post_fracture(network &n) {
       }
 
       if (in_bin) {
-        bin_counts[be]++;
+        sb[be]++;
       }
     }
   }
 
-  bin_counts[log2(L)]++;
+  sb[log2(L)]++;
 
+  // crack surface correlations
   std::vector<double> crack_damage(pow(L, 2), 0.0);
 
   double damage_tot = 0;
@@ -214,14 +383,19 @@ void ma::post_fracture(network &n) {
   }
 
   std::vector<fftw_complex> t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out);
-  std::vector<double> Ccc = correlation(L, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  correlation(L, Css, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
-  update_field_data("Ccc", damage_tot, Ccc, L, beta);
+  std::vector<bool> last_crack(pow(L, 2), 0);
 
-  for (auto e = avalanches.back().rbegin(); e != avalanches.back().rend(); e++) {
+  // 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;
   }
 
+
+  // cluster size distribution and cluster-cluster correlation
   num = boost::connected_components(G, &component[0]);
 
   for (unsigned int i = 0; i < num; i++) {
@@ -236,13 +410,28 @@ void ma::post_fracture(network &n) {
     }
 
     if (size > 0) {
-      cd[size - 1]++;
+      sc[size - 1]++;
       t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out);
-      Ccc = correlation(L, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+      correlation(L, Ccc, t_crack_damage, t_crack_damage, 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);
+
 
-      update_field_data("Cclcl", size, Ccc, L, beta);
+  // 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++;
     }
   }
 
+  sd[total_broken]++;
+
 }
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 817481b..4c05d93 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -7,19 +7,30 @@
 #include <string>
 #include <cinttypes>
 #include <sstream>
+#include <functional>
+#include <iostream>
+#include <valarray>
+#include <array>
 
 #include <boost/graph/adjacency_list.hpp>
 #include <boost/graph/connected_components.hpp>
+#include <boost/graph/depth_first_search.hpp>
+#include <boost/range/combine.hpp>
+#include <boost/foreach.hpp>
 
 #include <fftw3.h>
 
 #include <network.hpp>
 #include <hooks.hpp>
 
+struct EdgeProperties {
+  unsigned int index;
+};
+
+typedef boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, boost::no_property, EdgeProperties> Graph;
+
 class ma : public hooks {
   // need:
-  //  - measurements for correlation functions (nice looking?)
-  //  - interface for reading and writing from files (nice looking?)
   //  - interface for turning on and off specific measurements
   //
   private:
@@ -29,20 +40,28 @@ class ma : public hooks {
     double *fftw_reverse_out;
     fftw_plan forward_plan;
     fftw_plan reverse_plan;
-    boost::adjacency_list <boost::listS, boost::vecS, boost::undirectedS> G;
+    Graph G;
+    
+    // measurement storage
+    std::vector<uint64_t> sc; // cluster size distribution
+    std::vector<uint64_t> sa; // avalanche size distribution
+    std::vector<uint64_t> sd; // avalanche size distribution
+    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> Caa; // avalanche-avalanche correlations
+    std::vector<uint64_t> Cdd; // damage-damage distribution
+    uint64_t Nc;
+    uint64_t Na;
 
   public:
     unsigned int N;
     unsigned int L;
-    unsigned int as;
     double lv;
     double beta;
 
 
-    std::vector<unsigned int> ad;
-    std::vector<unsigned int> cd;
-    std::vector<unsigned int> bin_counts;
-    std::list<std::list<unsigned int>> avalanches;;
+    std::list<unsigned int> avalanche;
 
     ma(unsigned int N, unsigned int L, double beta);
     ~ma();