all threshold comparisons now done in log, so as to prevent underflows at small beta

3 files changed, 9 insertions, 8 deletions

diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index 849975f..d955e43 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -5,6 +5,7 @@
 #include <functional>
 #include <utility>
 #include <random>
+#include <limits>
 
 #include <cholmod.h>
 
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index bc7a0c6..9cb1007 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -89,10 +89,10 @@ void network::set_thresholds(double beta, std::mt19937& rng) {
   std::uniform_real_distribution<long double> d(0.0, 1.0);
 
   for (long double& threshold : thresholds) {
-    threshold = 0.0;
+    threshold = std::numeric_limits<long double>::lowest();
 
-    while (threshold == 0.0) {
-      threshold = expl(logl(d(rng)) / (long double)beta);
+    while (threshold == std::numeric_limits<long double>::lowest()) {
+      threshold = logl(d(rng)) / (long double)beta;
     }
   }
 }
@@ -202,11 +202,11 @@ void network::fracture(hooks& m, double cutoff) {
     }
 
     unsigned int max_pos = UINT_MAX;
-    long double max_val = 0;
+    long double max_val = std::numeric_limits<long double>::lowest();
 
     for (unsigned int i = 0; i < G.edges.size(); i++) {
       if (!fuses[i] && fabs(ci.currents[i]) > cutoff) {
-        long double val = (long double)fabs(ci.currents[i]) / thresholds[i];
+        long double val = logl(fabs(ci.currents[i])) - thresholds[i];
         if (val > max_val) {
           max_val = val;
           max_pos = i;
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 5dc6a64..98795f1 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -103,11 +103,11 @@ std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) {
       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);
+        double dx = fabs(n.G.dual_vertices[n.G.dual_edges[edge][0]].r.x - n.G.dual_vertices[n.G.dual_edges[edge][1]].r.x);
         if (dx > n.G.L.x / 2) {
           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);
+        double dy = fabs(n.G.dual_vertices[n.G.dual_edges[edge][0]].r.y - n.G.dual_vertices[n.G.dual_edges[edge][1]].r.y);
         if (dy > n.G.L.y / 2) {
           crossing_count[1]++;
         }
@@ -292,7 +292,7 @@ void ma::pre_fracture(const network &) {
 }
 
 void ma::bond_broken(const network& net, const current_info& cur, unsigned int i) {
-  long double c = cur.conductivity / fabs(cur.currents[i]) * net.thresholds[i];
+  long double c = logl(cur.conductivity / fabs(cur.currents[i])) + net.thresholds[i];
   if (c > lv && avalanches.back().size() > 0) {
     sa[avalanches.back().size() - 1]++;
     Na++;