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

10 files changed, 219 insertions, 30 deletions

diff --git a/src/anal_process.c b/src/anal_process.c
index c66aea0..0ebec84 100644
--- a/src/anal_process.c
+++ b/src/anal_process.c
@@ -37,6 +37,8 @@ int main(int argc, char *argv[]) {
 	double *errors_c2 = (double *)malloc(nc * sizeof(double));
 	double *vals_c3 = (double *)malloc(nc * sizeof(double));
 	double *errors_c3 = (double *)malloc(nc * sizeof(double));
+	double *vals_c4 = (double *)malloc(nc * sizeof(double));
+	double *errors_c4 = (double *)malloc(nc * sizeof(double));
 
 	char *out_filename = (char *)malloc(100 * sizeof(char));
 	uint_t out_filename_len = 0;
@@ -94,6 +96,11 @@ int main(int argc, char *argv[]) {
 			mom(dist_len, 3, dist, mom3);
 			vals_c3[i] = mom3[0];
 			errors_c3[i] = mom3[1];
+
+			double mom4[2];
+			mom(dist_len, 4, dist, mom4);
+			vals_c4[i] = mom4[0];
+			errors_c4[i] = mom4[1];
 		}
 		free(dist);
 	}
@@ -107,7 +114,7 @@ int main(int argc, char *argv[]) {
 	FILE *out_file = fopen(out_filename, "w");
 
 	for (uint_t i = 0; i < nc; i++) {
-		fprintf(out_file, "%u %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i]);
+		fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
 	}
 
 	fclose(out_file);
diff --git a/src/big_anal_process.c b/src/big_anal_process.c
index 8345785..0f7106f 100644
--- a/src/big_anal_process.c
+++ b/src/big_anal_process.c
@@ -53,6 +53,8 @@ int main(int argc, char *argv[]) {
 	double *errors_c2 = (double *)malloc(nc * sizeof(double));
 	double *vals_c3 = (double *)malloc(nc * sizeof(double));
 	double *errors_c3 = (double *)malloc(nc * sizeof(double));
+	double *vals_c4 = (double *)malloc(nc * sizeof(double));
+	double *errors_c4 = (double *)malloc(nc * sizeof(double));
 
 	char *out_filename = (char *)malloc(100 * sizeof(char));
 	uint_t out_filename_len = 0;
@@ -114,6 +116,11 @@ int main(int argc, char *argv[]) {
 			get_mom(dist_len, 3, dist, mom1, mom3);
 			vals_c3[i] = mom3[0];
 			errors_c3[i] = mom3[1];
+
+			double mom4[2];
+			get_mom(dist_len, 4, dist, mom1, mom4);
+			vals_c4[i] = mom4[0];
+			errors_c4[i] = mom4[1];
 		}
 		free(dist);
 	}
@@ -127,7 +134,7 @@ int main(int argc, char *argv[]) {
 	FILE *out_file = fopen(out_filename, "w");
 
 	for (uint_t i = 0; i < nc; i++) {
-		fprintf(out_file, "%u %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i]);
+		fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
 	}
 
 	fclose(out_file);
@@ -141,6 +148,8 @@ int main(int argc, char *argv[]) {
 	free(errors_c2);
 	free(vals_c3);
 	free(errors_c3);
+	free(vals_c4);
+	free(errors_c4);
 
 	return 0;
 }
diff --git a/src/break_edge.c b/src/break_edge.c
index 01570d8..538b88d 100644
--- a/src/break_edge.c
+++ b/src/break_edge.c
@@ -1,13 +1,21 @@
 
 #include "fracture.h"
 
-bool break_edge(net_t *instance, uint_t edge, cholmod_common *c) {
+bool break_edge(net_t *instance, uint_t edge, cholmod_common *c, bool refactor) {
 	instance->fuses[edge] = true;
 
 	if (instance->factor != NULL) {
-		uint_t w1 = instance->graph->ev_break[2 * edge];
-		uint_t w2 = instance->graph->ev_break[2 * edge + 1];
-		factor_update(instance->factor, w1, w2, c);
+		if (refactor) {
+			cholmod_sparse *laplacian = gen_laplacian(instance, c, true);
+			CHOL_F(free_factor)(&instance->factor, c);
+			instance->factor = CHOL_F(analyze)(laplacian, c);
+			CHOL_F(factorize)(laplacian, instance->factor, c);
+			CHOL_F(free_sparse)(&laplacian, c);
+		} else {
+			uint_t w1 = instance->graph->ev_break[2 * edge];
+			uint_t w2 = instance->graph->ev_break[2 * edge + 1];
+			factor_update(instance->factor, w1, w2, c);
+		}
 	}
 
 	uint_t v1, v2, s1, s2, dv1, dv2, ds1, ds2;
diff --git a/src/cen_anal_process.c b/src/cen_anal_process.c
new file mode 100644
index 0000000..98f68ac
--- /dev/null
+++ b/src/cen_anal_process.c
@@ -0,0 +1,154 @@
+
+#include "fracture.h"
+#include <gsl/gsl_sf_erf.h>
+#include <gsl/gsl_sf_laguerre.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+
+void get_mean(uint_t len, uint32_t *data, double result[2]) {
+	uint_t total = 0;
+	double mean = 0;
+	double mean_err = 0;
+
+	for (uint_t i = 0; i < len; i++) {
+		uint32_t datai = data[i];
+
+		total += datai;
+		mean += i * datai;
+		mean_err += gsl_pow_2(i) * datai; 
+	}
+
+	double meanf = mean / total;
+	double meanf_err = meanf * sqrt(mean_err / gsl_pow_2(mean) + 1 / total);
+
+	result[0] = meanf;
+	result[1] = meanf_err;
+}
+
+void get_mom(uint_t len, uint_t n, uint32_t *data, double mean[2], double result[2]) {
+	uint_t total = 0;
+	double mom = 0;
+	double mom_err = 0;
+
+	for (uint_t i = 0; i < len; i++) {
+		uint32_t datai = data[i];
+		double in = pow(fabs(((double)i) - mean[0]), n);
+
+		total += datai;
+		mom += in * datai;
+		mom_err += gsl_pow_2(in) * datai; 
+	}
+
+	double momf = mom / total;
+	double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total);
+
+	result[0] = momf;
+	result[1] = momf_err;
+}
+
+int main(int argc, char *argv[]) {
+	uint_t nc = argc - 1;
+	uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
+	double *betas_c = (double *)malloc(nc * sizeof(double));
+	double *vals_c1 = (double *)malloc(nc * sizeof(double));
+	double *errors_c1 = (double *)malloc(nc * sizeof(double));
+	double *vals_c2 = (double *)malloc(nc * sizeof(double));
+	double *errors_c2 = (double *)malloc(nc * sizeof(double));
+	double *vals_c3 = (double *)malloc(nc * sizeof(double));
+	double *errors_c3 = (double *)malloc(nc * sizeof(double));
+	double *vals_c4 = (double *)malloc(nc * sizeof(double));
+	double *errors_c4 = (double *)malloc(nc * sizeof(double));
+
+	char *out_filename = (char *)malloc(100 * sizeof(char));
+	uint_t out_filename_len = 0;
+
+	for (uint_t i = 0; i < nc; i++) {
+		char *fn = argv[1 + i];
+		char *buff = (char *)malloc(20 * sizeof(char));
+		uint_t pos = 0; uint_t c = 0;
+		while (c < 5) {
+			if (fn[pos] == '_') {
+				c++;
+			}
+			if (i == 0) {
+				out_filename[pos] = fn[pos];
+				out_filename_len++;
+			}
+			pos++;
+		}
+		c = 0;
+		while (fn[pos] != '_') {
+			buff[c] = fn[pos];
+			pos++;
+			c++;
+		}
+		buff[c] = '\0';
+		Ls_c[i] = atoi(buff);
+		c = 0;
+		pos++;
+		while (fn[pos] != '_') {
+			buff[c] = fn[pos];
+			pos++;
+			c++;
+		}
+		buff[c] = '\0';
+		betas_c[i] = atof(buff);
+		free(buff);
+
+		uint_t dist_len = 4 * pow(Ls_c[i], 2);
+		uint32_t *dist = malloc(dist_len * sizeof(uint32_t));
+		FILE *dist_file = fopen(fn, "rb");
+		fread(dist, sizeof(uint32_t), dist_len, dist_file);
+		fclose(dist_file);
+		{
+			double mom1[2];
+			get_mean(dist_len, dist, mom1);
+			vals_c1[i] = mom1[0];
+			errors_c1[i] = mom1[1];
+
+			double mom2[2];
+			get_mom(dist_len, 2, dist, mom1, mom2);
+			vals_c2[i] = mom2[0];
+			errors_c2[i] = mom2[1];
+
+			double mom3[2];
+			get_mom(dist_len, 3, dist, mom1, mom3);
+			vals_c3[i] = mom3[0];
+			errors_c3[i] = mom3[1];
+
+			double mom4[2];
+			get_mom(dist_len, 4, dist, mom1, mom4);
+			vals_c4[i] = mom4[0];
+			errors_c4[i] = mom4[1];
+		}
+		free(dist);
+	}
+
+	out_filename[out_filename_len-1] = '.';
+	out_filename[out_filename_len] = 't';
+	out_filename[out_filename_len+1] = 'x';
+	out_filename[out_filename_len+2] = 't';
+	out_filename[out_filename_len+3] = '\0';
+
+	FILE *out_file = fopen(out_filename, "w");
+
+	for (uint_t i = 0; i < nc; i++) {
+		fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
+	}
+
+	fclose(out_file);
+
+
+	free(Ls_c);
+	free(betas_c);
+	free(vals_c1);
+	free(errors_c1);
+	free(vals_c2);
+	free(errors_c2);
+	free(vals_c3);
+	free(errors_c3);
+
+	return 0;
+}
+
diff --git a/src/corr_test.c b/src/corr_test.c
index b6f8a18..4dc7efc 100644
--- a/src/corr_test.c
+++ b/src/corr_test.c
@@ -9,7 +9,7 @@ int main() {
 
 	graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 128, false, &c);
 	net_t *instance = net_create(network, 1e-14, 3, 0, true, &c);
-	net_fracture(instance, &c, 1e-10);
+	net_fracture(instance, &c, 1e-10, 40);
 
 	double *corr = get_corr(instance, NULL, &c);
 
diff --git a/src/fracture.c b/src/fracture.c
index 11e6efd..7ee82fb 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -8,7 +8,7 @@ int main(int argc, char *argv[]) {
 	// defining variables to be (potentially) set by command line flags
 	uint8_t filename_len;
 	uint32_t N;
-	uint_t L;
+	uint_t L, refactor_every;
 	double beta, inf, cutoff, crack_len;
 	bool crack_growth_crit, save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
 			 save_crit_stress, save_stress_field, save_voltage_field, save_toughness, save_energy, save_conductivity,
@@ -25,6 +25,7 @@ int main(int argc, char *argv[]) {
 	// set default values
 
 	N = 100;
+	refactor_every = UINT_MAX;
 	L = 16;
 	crack_len = 0.;
 	beta = .3;
@@ -45,7 +46,7 @@ int main(int argc, char *argv[]) {
 	save_damage_field = false;
 	save_conductivity = false;
 	save_toughness = false;
-	save_energy = false;
+	save_energy = true;
 	save_threshold = false;
 
 
@@ -58,11 +59,14 @@ int main(int argc, char *argv[]) {
 
 	// get commandline options
 
-	while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:gTE")) != -1) {
+	while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:gTER:")) != -1) {
 		switch (opt) {
 			case 'n':
 				N = atoi(optarg);
 				break;
+			case 'R':
+				refactor_every = atoi(optarg);
+				break;
 			case 'L':
 				L = atoi(optarg);
 				break;
@@ -303,7 +307,7 @@ int main(int argc, char *argv[]) {
 		graph_t *g = graph_create(lattice, boundary, L, use_dual, &c);
 		net_t *net = net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c);
 		net_t *tmp_net = net_copy(net, &c);
-		data_t *data = net_fracture(tmp_net, &c, cutoff);
+		data_t *data = net_fracture(tmp_net, &c, cutoff, refactor_every);
 		net_free(tmp_net, &c);
 
 		uint_t max_pos = 0;
@@ -357,7 +361,8 @@ int main(int argc, char *argv[]) {
 				break;
 			}
 
-			break_edge(net, next_broken, &c);
+			bool refactor = ((j + 1) % refactor_every) == 0;
+			break_edge(net, next_broken, &c, refactor);
 
 			double val = data->extern_field[j];
 			if (save_cluster_dist) {
diff --git a/src/fracture.h b/src/fracture.h
index 0b078ff..f122869 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -72,7 +72,7 @@ typedef struct {
 typedef struct {
 	const graph_t *graph;
 	bool *fuses;
-	double *thres;
+	long double *thres;
 	double inf;
 	cholmod_dense *boundary_cond;
 	cholmod_factor *factor;
@@ -117,7 +117,7 @@ double dual_vert_to_coord(uint_t width, bool periodic, uint_t vert,
 void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common *c);
 
 void net_notch(net_t *net, double notch_len, cholmod_common *c);
-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);
 double *net_voltages(const net_t *net, cholmod_common *c);
 double *net_currents(const net_t *net, const double *voltages, cholmod_common *c);
 double net_conductivity(const net_t *net, const double *voltages);
@@ -148,7 +148,7 @@ graph_t *ini_voro_graph(uint_t L, bound_t boundary, bool use_dual,
 													double *(*genfunc)(uint_t, bound_t, gsl_rng *, uint_t *),
 													cholmod_common *c);
 
-bool break_edge(net_t *instance, uint_t edge, cholmod_common *c);
+bool break_edge(net_t *instance, uint_t edge, cholmod_common *c, bool refactor);
 
 void finish_instance(net_t *instance, cholmod_common *c);
 
@@ -183,4 +183,4 @@ bool set_connected(const cholmod_sparse *laplacian, uint_t *marks, int vertex, i
 
 unsigned long int rand_seed();
 
-double rand_dist_pow(const gsl_rng *r, double beta);
+long double rand_dist_pow(const gsl_rng *r, double beta);
diff --git a/src/net.c b/src/net.c
index f3aeda9..444e76c 100644
--- a/src/net.c
+++ b/src/net.c
@@ -1,8 +1,8 @@
 
 #include "fracture.h"
 
-double *get_thres(uint_t ne, double beta) {
-	double *thres = (double *)malloc(ne * sizeof(double));
+long double *get_thres(uint_t ne, double beta) {
+	long double *thres = (long double *)malloc(ne * sizeof(long double));
 	assert(thres != NULL);
 
 	gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN);
@@ -39,7 +39,7 @@ void net_notch(net_t *net, double notch_len, cholmod_common *c) {
 		correct_length = v1x < notch_len && v2x < notch_len;
 
 		if (crosses_center && not_wrapping && correct_length) {
-			break_edge(net, i, c);
+			break_edge(net, i, c, false);
 		}
 	}
 }
@@ -96,8 +96,8 @@ net_t *net_copy(const net_t *net, cholmod_common *c) {
 	assert(net_copy->fuses != NULL);
 	memcpy(net_copy->fuses, net->fuses, fuses_size);
 
-	size_t thres_size = (net->graph)->ne * sizeof(double);
-	net_copy->thres = (double *)malloc(thres_size);
+	size_t thres_size = (net->graph)->ne * sizeof(long double);
+	net_copy->thres = (long double *)malloc(thres_size);
 	assert(net_copy->thres != NULL);
 	memcpy(net_copy->thres, net->thres, thres_size);
 
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;
diff --git a/src/rand.c b/src/rand.c
index 1a6a3d4..75722ac 100644
--- a/src/rand.c
+++ b/src/rand.c
@@ -9,12 +9,15 @@ unsigned long int rand_seed() {
 	return seed;
 }
 
-double rand_dist_pow(const gsl_rng *r, double beta) {
-	double x = 0;
+long double rand_dist_pow(const gsl_rng *r, double beta) {
+	long double x = 0;
 
 	// underflow means that for very small beta x is sometimes identically zero,
 	// which causes problems
-	while ((x = exp(log(gsl_rng_uniform_pos(r)) / beta)) == 0.0);
+	while (x == 0.0) {
+		long double y = logl(gsl_rng_uniform_pos(r)) / beta;
+		x = expl(y);
+	}
 
 	return x;
 }