summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2018-12-06 23:25:46 -0500
committerJaron Kent-Dobias <jaron@kent-dobias.com>2018-12-06 23:25:46 -0500
commitf1825b6ca420f63e17ca69e3c9412b80adcdbb1c (patch)
tree5927c43e7f74899d1f333834eefee9e8193efd61 /lib
parent869df8ab1856fa36eaca09c15582c2fb1335a64c (diff)
downloadfuse_networks-f1825b6ca420f63e17ca69e3c9412b80adcdbb1c.tar.gz
fuse_networks-f1825b6ca420f63e17ca69e3c9412b80adcdbb1c.tar.bz2
fuse_networks-f1825b6ca420f63e17ca69e3c9412b80adcdbb1c.zip
all threshold comparisons now done in log, so as to prevent underflows at small beta
Diffstat (limited to 'lib')
-rw-r--r--lib/include/network.hpp1
-rw-r--r--lib/src/network.cpp10
2 files changed, 6 insertions, 5 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;