a lot of cleaning and simplification

8 files changed, 69 insertions, 89 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6746aa9..a804ac5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,8 +1,8 @@
 
-cmake_minimum_required(VERSION 3.9)
+cmake_minimum_required(VERSION 3.8)
 
-set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall")
-set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto")
+set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall -stdlib=libc++")
+set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -stdlib=libc++")
 
 set (CMAKE_CXX_STANDARD 17)
 set (CMAKE_C_STANDARD 11)
diff --git a/lib/include/current_info.hpp b/lib/include/current_info.hpp
new file mode 100644
index 0000000..90dd3c7
--- /dev/null
+++ b/lib/include/current_info.hpp
@@ -0,0 +1,10 @@
+
+#pragma once
+
+#include <vector>
+
+typedef struct current_info {
+  double conductivity;
+  std::vector<double> currents;
+} current_info;
+
diff --git a/lib/include/hooks.hpp b/lib/include/hooks.hpp
index 350f318..e617ac3 100644
--- a/lib/include/hooks.hpp
+++ b/lib/include/hooks.hpp
@@ -1,12 +1,14 @@
 
 #pragma once
 
+#include "current_info.hpp"
+
 class network;
 
 class hooks {
   public:
     virtual void pre_fracture(const network&) {};
-    virtual void bond_broken(const network&, const std::pair<double, std::vector<double>>&, unsigned int) {};
+    virtual void bond_broken(const network&, const current_info&, unsigned int) {};
     virtual void post_fracture(network&) {}; // post fracture hook can be destructive
 };
 
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index 5cbba0d..e6d67b9 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -6,10 +6,11 @@
 #include <utility>
 #include <random>
 
-#include "cholmod.h"
+#include <cholmod.h>
 
 #include "graph.hpp"
 #include "hooks.hpp"
+#include "current_info.hpp"
 
 #ifdef FRACTURE_LONGINT
 
@@ -40,9 +41,8 @@ class network {
     ~network();
 
     void set_thresholds(double, std::mt19937&);
-    void break_edge(unsigned int);
-    void add_edge(unsigned int);
-    std::pair<double, std::vector<double>> currents();
+    void break_edge(unsigned int, bool unbreak = false);
+    current_info get_current_info();
     void fracture(hooks&, double cutoff = 1e-14);
 };
 
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index a27ca39..9f6cb0d 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -56,7 +56,6 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
 
   t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c);
 
-
   t->nnz = 2 * G.edges.size();
 
   for (unsigned int i = 0; i < G.edges.size(); i++) {
@@ -89,68 +88,17 @@ network::~network() {
 void network::set_thresholds(double beta, std::mt19937& rng) {
   std::uniform_real_distribution<long double> d(0.0, 1.0);
 
-  for (unsigned int i = 0; i < G.edges.size(); i++) {
-    long double x = 0.0;
+  for (long double& threshold : thresholds) {
+    threshold = 0.0;
 
-    while (x == 0.0) {
-      x = expl(logl(d(rng)) / (long double)beta);
+    while (threshold == 0.0) {
+      threshold = expl(logl(d(rng)) / (long double)beta);
     }
-
-    thresholds[i] = x;
-  }
-}
-
-void network::add_edge(unsigned int e) {
-  fuses[e] = false;
-  unsigned int v0 = G.edges[e][0];
-  unsigned int v1 = G.edges[e][1];
-
-  unsigned int n = factor->n;
-
-  cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
-
-  unsigned int s1, s2;
-  s1 = v0 < v1 ? v0 : v1;
-  s2 = v0 < v1 ? v1 : v0;
-
-  CHOL_INT *pp = (CHOL_INT *)update_mat->p;
-  CHOL_INT *ii = (CHOL_INT *)update_mat->i;
-  double *xx = (double *)update_mat->x;
-
-  for (unsigned int i = 0; i <= s1; i++) {
-    pp[i] = 0;
-  }
-
-  for (unsigned int i = s1 + 1; i <= n; i++) {
-    pp[i] = 2;
-  }
-
-  ii[0] = s1;
-  ii[1] = s2;
-  xx[0] = -1;
-  xx[1] = 1;
-
-  cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
-      update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c);
-
-  CHOL_F(updown)(false, perm_update_mat, factor, c);
-
-  CHOL_F(free_sparse)(&perm_update_mat, c);
-  CHOL_F(free_sparse)(&update_mat, c);
-
-  double v0y = G.vertices[v0].r.y;
-  double v1y = G.vertices[v1].r.y;
-
-  if (fabs(v0y - v1y) > 0.5) {
-    bool ind = v0y < v1y ? 0 : 1;
-
-    ((double *)b->x)[G.edges[e][ind]] += 1.0;
-    ((double *)b->x)[G.edges[e][!ind]] -= 1.0;
   }
 }
 
-void network::break_edge(unsigned int e) {
-  fuses[e] = true;
+void network::break_edge(unsigned int e, bool unbreak) {
+  fuses[e] = !unbreak;
   unsigned int v0 = G.edges[e][0];
   unsigned int v1 = G.edges[e][1];
 
@@ -182,7 +130,7 @@ void network::break_edge(unsigned int e) {
   cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
       update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c);
 
-  CHOL_F(updown)(false, perm_update_mat, factor, c);
+  CHOL_F(updown)(unbreak, perm_update_mat, factor, c);
 
   CHOL_F(free_sparse)(&perm_update_mat, c);
   CHOL_F(free_sparse)(&update_mat, c);
@@ -191,14 +139,14 @@ void network::break_edge(unsigned int e) {
   double v1y = G.vertices[v1].r.y;
 
   if (fabs(v0y - v1y) > 0.5) {
-    bool ind = v0y < v1y ? 0 : 1;
+    bool ind = v0y < v1y ? unbreak : !unbreak;
 
     ((double *)b->x)[G.edges[e][ind]] -= 1.0;
     ((double *)b->x)[G.edges[e][!ind]] += 1.0;
   }
 }
 
-std::pair<double, std::vector<double>> network::currents() {
+current_info network::get_current_info() {
   cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
 
   if (((double *)x->x)[0] != ((double *)x->x)[0]) {
@@ -240,16 +188,16 @@ std::pair<double, std::vector<double>> network::currents() {
   CHOL_F(free_dense)(&x, c);
   CHOL_F(free_dense)(&y, c);
 
-  return std::make_pair(total_current, currents);
+  return {total_current, currents};
 }
 
 void network::fracture(hooks& m, double cutoff) {
   m.pre_fracture(*this);
 
   while (true) {
-    auto currents = this->currents();
+    current_info ci = this->get_current_info();
 
-    if (currents.first < cutoff * G.vertices.size()) {
+    if (ci.conductivity < cutoff * G.vertices.size()) {
       break;
     }
 
@@ -257,8 +205,8 @@ void network::fracture(hooks& m, double cutoff) {
     long double max_val = 0;
 
     for (unsigned int i = 0; i < G.edges.size(); i++) {
-      if (!fuses[i] && fabs(currents.second[i]) > cutoff) {
-        long double val = (long double)fabs(currents.second[i]) / thresholds[i];
+      if (!fuses[i] && fabs(ci.currents[i]) > cutoff) {
+        long double val = (long double)fabs(ci.currents[i]) / thresholds[i];
         if (val > max_val) {
           max_val = val;
           max_pos = i;
@@ -272,7 +220,7 @@ void network::fracture(hooks& m, double cutoff) {
 
     this->break_edge(max_pos);
 
-    m.bond_broken(*this, currents, max_pos);
+    m.bond_broken(*this, ci, max_pos);
   }
 
   m.post_fracture(*this);
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 77af253..0aad56e 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -11,7 +11,23 @@
 #include <hooks.hpp>
 #include "measurements.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 int N = 1;
@@ -36,12 +52,10 @@ int main(int argc, char* argv[]) {
 
   cholmod_common c;
   CHOL_F(start)(&c);
-  c.supernodal = CHOLMOD_SUPERNODAL;
-
 
   graph G(L);
   network base_network(G, &c);
-  ma meas(N, L, beta);
+  ma meas(L, beta);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
@@ -50,6 +64,9 @@ int main(int argc, char* argv[]) {
     network tmp_network(base_network);
     tmp_network.set_thresholds(beta, rng);
     tmp_network.fracture(meas);
+
+    if (quit.load())
+      break;
   }
 
   CHOL_F(finish)(&c);
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 0bd72c4..ed96855 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -38,8 +38,8 @@ std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) {
 
   class find_cycle : public boost::default_dfs_visitor {
     public:
-      unsigned int end;
       std::list<unsigned int>& E;
+      unsigned int end;
       struct done{};
 
       find_cycle(std::list<unsigned int>& E, unsigned int end) : E(E), end(end) {}
@@ -226,7 +226,7 @@ void autocorrelation(unsigned int L, std::vector<T>& out_data, fftw_plan forward
   }
 }
 
-ma::ma(unsigned int N, unsigned int L, double beta) : L(L), G(2 * pow(L / 2, 2)), N(N), beta(beta),
+ma::ma(unsigned int L, double beta) : L(L), beta(beta), G(2 * pow(L / 2, 2)),
   sc(pow(L, 2), 0),
   sa(pow(L, 2), 0),
   sd(pow(L, 2), 0),
@@ -239,6 +239,7 @@ ma::ma(unsigned int N, unsigned int L, double beta) : L(L), G(2 * pow(L / 2, 2))
   Cll(pow(L / 2 + 1, 2), 0),
   Cee(pow(L / 2 + 1, 2), 0)
 {
+  N = 0;
   Nc = 0;
   Na = 0;
 
@@ -284,8 +285,9 @@ void ma::pre_fracture(const network &) {
   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) {
+void ma::bond_broken(const network& net, const current_info& cur, unsigned int i) {
+  double c = cur.conductivity / fabs(cur.currents[i]) * net.thresholds[i];
+  if (c > lv) {
     sa[avalanches.back().size()]++;
     Na++;
 
@@ -297,7 +299,7 @@ void ma::bond_broken(const network& net, const std::pair<double, std::vector<dou
 
     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];
+    lv = c;
     avalanches.push_back({i});
   } else {
     avalanches.back().push_back(i);
@@ -373,7 +375,7 @@ void ma::post_fracture(network &n) {
   // rewind the last avalanche
   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);
+    n.break_edge(e, true);
     fftw_forward_in[e] = 1;
   }
 
@@ -417,5 +419,6 @@ void ma::post_fracture(network &n) {
 
   sd[total_broken]++;
 
+  N++;
 }
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 621c425..ab31ccd 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -41,6 +41,9 @@ class ma : public hooks {
     double *fftw_reverse_out;
     fftw_plan forward_plan;
     fftw_plan reverse_plan;
+    unsigned int N;
+    unsigned int L;
+    double beta;
     Graph G;
     
     // measurement storage
@@ -59,19 +62,16 @@ class ma : public hooks {
     uint64_t Na;
 
   public:
-    unsigned int N;
-    unsigned int L;
     double lv;
-    double beta;
 
 
     std::list<std::list<unsigned int>> avalanches;
 
-    ma(unsigned int N, unsigned int L, double beta);
+    ma(unsigned int L, double beta);
     ~ma();
 
     void pre_fracture(const network &) override;
-    void bond_broken(const network& net, const std::pair<double, std::vector<double>>& cur, unsigned int i) override;
+    void bond_broken(const network& net, const current_info& cur, unsigned int i) override;
     void post_fracture(network &n) override;
 };