various changes, including adding central moments and changing the fuse thresholds to long doubles

1 files changed, 9 insertions, 6 deletions

diff --git a/src/net_fracture.c b/src/net_fracture.c
index dcf67c9..3bfcd83 100644
--- a/src/net_fracture.c
+++ b/src/net_fracture.c
@@ -3,14 +3,14 @@
 
 uint_t get_next_broken(net_t *net, double *currents, double cutoff) {
 	uint_t max_pos = UINT_MAX;
-	double max_val = 0;
+	long double max_val = 0;
 
 	for (uint_t i = 0; i < net->graph->ne; i++) {
-		double current = fabs(currents[i]);
+		long double current = fabs(currents[i]);
 		bool broken = net->fuses[i];
 
 		if (!broken && current > cutoff) {
-			double val = current / net->thres[i];
+			long double val = current / net->thres[i];
 
 			if (val > max_val) {
 				max_val = val;
@@ -27,10 +27,12 @@ uint_t get_next_broken(net_t *net, double *currents, double cutoff) {
 	return max_pos;
 }
 
-data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
+data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff, uint_t refactor_every) {
 	data_t *data = data_create(net->graph->ne);
 
+	uint_t n = 0;
 	while (true) {
+		n++;
 		double *voltages = net_voltages(net, c);
 		double *currents = net_currents(net, voltages, c);
 
@@ -52,12 +54,13 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
 			sim_current = 1;
 		}
 
-		data_update(data, last_broke, fabs(sim_current * (net->thres)[last_broke] / currents[last_broke]), conductivity);
+		data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] / currents[last_broke]), conductivity);
 
 		free(voltages);
 		free(currents);
 
-		break_edge(net, last_broke, c);
+		bool refactor = (n % refactor_every) == 0;
+		break_edge(net, last_broke, c, refactor);
 
 		if (net->num_components > 1 && net->graph->boundary == TORUS_BOUND) {
 			break;