added support for running percolation using the fuse network rules

6 files changed, 419 insertions, 0 deletions

diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index dd11342..15cb8ea 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -82,3 +82,18 @@ class elastic_network : public network {
     current_info get_current_info() override;
 };
 
+class percolation_network : public network {
+  private:
+    double weight;
+  public:
+    problem hook_x;
+    problem hook_y;
+
+    percolation_network(const graph&, cholmod_common*);
+    percolation_network(const percolation_network&);
+
+    void fracture(hooks&, double weight = 0.5, double cutoff = 1e-10);
+    void break_edge(unsigned, bool unbreak = false) override;
+    current_info get_current_info() override;
+};
+
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index ed51c8e..1d9b3f9 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -363,3 +363,83 @@ void elastic_network::break_edge(unsigned e, bool unbreak) {
   hook_y.break_edge(e, unbreak);
 }
 
+percolation_network::percolation_network(const graph& G, cholmod_common *c) : network(G), hook_x(G, 1, c), hook_y(G, 0, c) {
+}
+
+percolation_network::percolation_network(const percolation_network& n) : network(n), hook_x(n.hook_x), hook_y(n.hook_y) {
+}
+
+void percolation_network::fracture(hooks& m, double weight, double cutoff) {
+  this->weight = weight;
+  m.pre_fracture(*this);
+
+  while (true) {
+    current_info ctot = this->get_current_info();
+
+    if (ctot.conductivity == 0) {
+      break;
+    }
+
+    unsigned max_pos = UINT_MAX;
+    long double max_val = std::numeric_limits<long double>::lowest();
+
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      if (!fuses[i] && ctot.currents[i] > 1.0 / G.edges.size()) {
+        long double val = - thresholds[i];
+        if (val > max_val) {
+          max_val = val;
+          max_pos = i;
+        }
+      }
+    }
+
+    if (max_pos == UINT_MAX)  {
+      throw nofuseex;
+    }
+
+    this->break_edge(max_pos);
+
+    m.bond_broken(*this, ctot, max_pos);
+  }
+
+  m.post_fracture(*this);
+}
+
+current_info percolation_network::get_current_info() {
+  current_info cx = hook_x.solve(fuses);
+  current_info cy = hook_y.solve(fuses);
+
+  bool done_x = cx.conductivity < 1.0 / G.edges.size();
+  bool done_y = cy.conductivity < 1.0 / G.edges.size();
+
+  current_info ctot;
+  ctot.currents.resize(G.edges.size());
+
+  if (done_x && done_y) {
+    ctot.conductivity = 0;
+  } else if (done_x) {
+    ctot.conductivity = cy.conductivity;
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity;
+    }
+  } else if (done_y) {
+    ctot.conductivity = cx.conductivity;
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity;
+    }
+  } else {
+    ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); 
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity, 2) + pow(weight * cy.currents[i] / cy.conductivity, 2));
+    }
+  }
+
+  return ctot;
+}
+
+void percolation_network::break_edge(unsigned e, bool unbreak) {
+  fuses[e] = !unbreak;
+  hook_x.break_edge(e, unbreak);
+  hook_y.break_edge(e, unbreak);
+}
+
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index eecbf72..9a7ce20 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -5,6 +5,9 @@ target_link_libraries(analysis_tools frac)
 add_library(measurements measurements.cpp)
 target_link_libraries(measurements frac analysis_tools)
 
+add_library(perc_meas perc_meas.cpp)
+target_link_libraries(perc_meas frac analysis_tools)
+
 add_library(animate animate.cpp)
 target_link_libraries(animate frac analysis_tools GL GLU glut)
 
@@ -29,3 +32,5 @@ target_link_libraries(sample_fracture frac sample cholmod)
 add_executable(sample_fracture_square sample_fracture_square.cpp)
 target_link_libraries(sample_fracture_square frac sample cholmod)
 
+add_executable(percolation percolation.cpp)
+target_link_libraries(percolation frac perc_meas fftw3 cholmod profiler GL GLU glut)
diff --git a/src/perc_meas.cpp b/src/perc_meas.cpp
new file mode 100644
index 0000000..4220a93
--- /dev/null
+++ b/src/perc_meas.cpp
@@ -0,0 +1,167 @@
+
+#include "perc_meas.hpp"
+#include <iostream>
+#include <cstdio>
+
+class badcycleException: public std::exception
+{
+  virtual const char* what() const throw()
+  {
+    return "Could not find a valid cycle on the broken system.";
+  }
+} badcycleex;
+
+void update_distribution_file(std::string id, const std::vector<uint64_t>& data, std::string model_string) {
+  std::string filename = model_string + id + ".dat";
+  std::ifstream file(filename);
+
+  std::vector<uint64_t> data_old(data.size(), 0);
+
+  if (file.is_open()) {
+    for (unsigned i = 0; i < data.size(); i++) {
+      uint64_t num;
+      file >> num;
+      data_old[i] = num;
+    }
+
+    file.close();
+  }
+
+  std::ofstream file_out(filename);
+
+  for (unsigned i = 0; i < data.size(); i++) {
+    file_out <<std::fixed<< data_old[i] + data[i] << " ";
+  }
+
+  file_out.close();
+}
+
+void update_distribution_file(std::string id, const std::vector<std::vector<uint64_t>>& data, std::string model_string) {
+  std::string filename = model_string + id + ".dat";
+  std::ifstream file(filename);
+
+  std::vector<std::vector<uint64_t>> data_old(data.size());
+
+  for (unsigned i = 0; i < data.size(); i++) {
+    data_old[i].resize(data[0].size());
+  }
+
+  if (file.is_open()) {
+    for (unsigned i = 0; i < data.size(); i++) {
+      for (unsigned j = 0; j < data[0].size(); j++) {
+        uint64_t num;
+        file >> num;
+        data_old[i][j] = num;
+      }
+    }
+
+    file.close();
+  }
+
+  std::ofstream file_out(filename);
+
+  for (unsigned i = 0; i < data.size(); i++) {
+    for (unsigned j = 0; j < data[0].size(); j++) {
+      file_out << std::fixed << data_old[i][j] + data[i][j] << " ";
+    }
+    file_out << "\n";
+  }
+
+  file_out.close();
+}
+
+pm::pm(unsigned n, double a) :
+  G(2 * n),
+  rank(2 * n),
+  parent(2 * n),
+  ds(&rank[0], &parent[0]),
+  sm(2 * n),
+  sl(2 * n),
+  sn(3 * n),
+  ss(2 * n),
+  sb(3 * n),
+  sp(3 * n)
+{
+  model_string = "percolation_" + std::to_string(n) + "_" + std::to_string(a) + "_";
+  for (std::vector<uint64_t> &x : sn) {
+    x.resize(2 * n);
+  }
+  for (std::vector<uint64_t> &x : sb) {
+    x.resize(3 * n);
+  }
+}
+
+pm::~pm() {
+  update_distribution_file("sm", sm, model_string);
+  update_distribution_file("sl", sl, model_string);
+  update_distribution_file("sn", sn, model_string);
+  update_distribution_file("ss", ss, model_string);
+  update_distribution_file("sb", sb, model_string);
+  update_distribution_file("sp", sp, model_string);
+}
+
+void pm::pre_fracture(const network&) {
+  boost::remove_edge_if(trivial, G);
+  initialize_incremental_components(G, ds);
+  incremental_components(G, ds);
+  p = 0;
+}
+
+void pm::bond_broken(const network& net, const current_info& cur, unsigned i) {
+  boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G);
+  ds.union_set(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1]);
+
+  boost::component_index<VertexIndex> components(parent.begin(), parent.end());
+  std::vector<unsigned> counts(components.size());
+
+  for (unsigned i = 0; i < components.size(); i++) {
+    unsigned comp_size = 0;
+    BOOST_FOREACH(VertexIndex child_index, components[i]) {
+      comp_size++;
+    }
+
+    sn[p][comp_size - 1]++;
+  }
+
+  unsigned bb_size = 0;
+
+  for (unsigned i = 0; i <net.G.edges.size(); i++) {
+    if (cur.currents[i] > 1.0 / net.G.edges.size()) {
+      bb_size++;
+    }
+  }
+
+  sb[p][bb_size - 1]++;
+  p++;
+}
+
+void pm::post_fracture(network &n) {
+  auto post_cracks = find_minimal_crack(G, n);
+  std::vector<unsigned> component(boost::num_vertices(G));
+  unsigned num = boost::connected_components(G, &component[0]);
+  if (post_cracks.size() > 2 || post_cracks.size() == 0) {
+    throw badcycleex;
+  }
+  for (auto c : post_cracks) {
+    sl[c.second.size() - 1]++;
+  }
+  unsigned crack_component = component[n.G.dual_edges[post_cracks.front().second.front()].v[0]];
+
+  std::vector<std::list<graph::coordinate>> components(num);
+
+  for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
+    components[component[i]].push_back(n.G.dual_vertices[i].r);
+  }
+
+  for (unsigned i = 0; i < num; i++) {
+    if (i != crack_component) {
+      sm[components[i].size() - 1]++;
+    } else {
+      ss[components[i].size() - 1]++;
+    }
+    sn[p][components[i].size() - 1]++;
+  }
+
+  sp[p - 1]++;
+}
+
diff --git a/src/perc_meas.hpp b/src/perc_meas.hpp
new file mode 100644
index 0000000..b85e27e
--- /dev/null
+++ b/src/perc_meas.hpp
@@ -0,0 +1,52 @@
+
+#include <cstring>
+#include <fstream>
+#include <string>
+#include <cinttypes>
+#include <sstream>
+#include <functional>
+#include <iostream>
+#include <array>
+
+#include <fftw3.h>
+
+#include <hooks.hpp>
+#include <boost/graph/graph_utility.hpp>
+#include <boost/graph/incremental_components.hpp>
+#include <boost/pending/disjoint_sets.hpp>
+#include <boost/foreach.hpp>
+
+#include "analysis_tools.hpp"
+typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;
+typedef boost::graph_traits<Graph>::vertices_size_type VertexIndex;
+
+
+class pm : public hooks {
+  private:
+    unsigned p;
+    Graph G;
+    std::vector<VertexIndex> rank;
+    std::vector<Vertex> parent;
+    boost::disjoint_sets<VertexIndex*, Vertex*> ds;
+
+    // measurement storage
+    std::vector<uint64_t> sc; // non-spanning cluster size distribution
+    std::vector<std::vector<uint64_t>> sn; // non-spanning cluster size distribution
+    std::vector<uint64_t> ss; // minimal spanning cluster size distribution
+    std::vector<uint64_t> sm; // spanning cluster size distribution
+    std::vector<uint64_t> sl; // final avalanche size distribution
+    std::vector<std::vector<uint64_t>> sb; // final avalanche size distribution
+    std::vector<uint64_t> sp;
+
+  public:
+    std::string model_string;
+
+    pm(unsigned Lx, unsigned Ly);
+    pm(unsigned n, double a);
+    ~pm();
+
+    void pre_fracture(const network &) override;
+    void bond_broken(const network& net, const current_info& cur, unsigned i) override;
+    void post_fracture(network &n) override;
+};
+
diff --git a/src/percolation.cpp b/src/percolation.cpp
new file mode 100644
index 0000000..03da117
--- /dev/null
+++ b/src/percolation.cpp
@@ -0,0 +1,100 @@
+
+#include <random>
+#include <iostream>
+
+#include <cholmod.h>
+
+#include "randutils/randutils.hpp"
+
+#include <graph.hpp>
+#include <network.hpp>
+#include <hooks.hpp>
+#include "perc_meas.hpp"
+
+#include <csignal>
+#include <cstring>
+#include <atomic>
+
+std::atomic<bool> quit(false);    // signal flag
+
+void got_signal(int) {
+    quit.store(true);
+}
+
+int main(int argc, char* argv[]) {
+  struct sigaction sa;
+  memset( &sa, 0, sizeof(sa) );
+  sa.sa_handler = got_signal;
+  sigfillset(&sa.sa_mask);
+  sigaction(SIGINT, &sa, NULL);
+
+  int opt;
+
+  unsigned N = 1;
+  double Lx = 16;
+  double Ly = 16;
+  double  beta = 0.0001;
+
+  unsigned n = 128;
+  double a = 1.0;
+  bool use_aN = false;
+
+  while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:")) != -1) {
+    switch (opt) {
+      case 'N':
+        N = (unsigned)atof(optarg);
+        break;
+      case 'X':
+        Lx = atof(optarg);
+        break;
+      case 'Y':
+        Ly = atof(optarg);
+        break;
+      case 'b':
+        beta = atof(optarg);
+        break;
+      case 'n':
+        n = (unsigned)atof(optarg);
+        use_aN = true;
+        break;
+      case 'a':
+        a = atof(optarg);
+        use_aN = true;
+        break;
+      default:
+        exit(1);
+    }
+  }
+
+  cholmod_common c;
+  CHOL_F(start)(&c);
+
+  randutils::auto_seed_128 seeds;
+  std::mt19937 rng{seeds};
+
+  if (use_aN) {
+    pm meas(n, a);
+
+    for (unsigned trial = 0; trial < N; trial++) {
+      while (true) {
+        try {
+          graph G(n, a, rng);
+          percolation_network fuse_network(G, &c);
+          fuse_network.set_thresholds(beta, rng);
+          fuse_network.fracture(meas);
+          break;
+        } catch (std::exception &e) {
+          std::cout << e.what() << '\n';
+        }
+      }
+
+      if (quit.load())
+        break;
+    }
+  }
+
+  CHOL_F(finish)(&c);
+
+  return 0;
+}
+