restarted all data generation to correct for errors, made certain data (stress, energy, thresholds) longs

5 files changed, 182 insertions, 92 deletions

diff --git a/src/data.c b/src/data.c
index b54236d..2047c44 100644
--- a/src/data.c
+++ b/src/data.c
@@ -10,7 +10,7 @@ data_t *data_create(uint_t ne) {
 	data->break_list = (uint_t *)malloc(ne * sizeof(uint_t));
 	assert(data->break_list != NULL);
 
-	data->extern_field = (double *)malloc(ne * sizeof(double));
+	data->extern_field = (long double *)malloc(ne * sizeof(long double));
 	assert(data->extern_field != NULL);
 
 	data->conductivity = (double *)malloc(ne * sizeof(double));
@@ -26,7 +26,7 @@ void data_free(data_t *data) {
 	free(data);
 }
 
-void data_update(data_t *data, uint_t last_broke, double strength, double conductivity) {
+void data_update(data_t *data, uint_t last_broke, long double strength, double conductivity) {
 	data->break_list[data->num_broken] = last_broke;
 	data->extern_field[data->num_broken] = strength;
 	data->conductivity[data->num_broken] = conductivity;
diff --git a/src/fracture.c b/src/fracture.c
index 7ee82fb..59bd424 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -10,8 +10,8 @@ int main(int argc, char *argv[]) {
 	uint32_t N;
 	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,
+	bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
+			 save_crit_stress, save_stress_field, save_voltage_field, save_energy, save_conductivity,
 			 save_damage, save_damage_field, save_threshold;
 	bound_t boundary;
 	lattice_t lattice;
@@ -33,7 +33,6 @@ int main(int argc, char *argv[]) {
 	cutoff = 1e-9;
 	boundary = FREE_BOUND;
 	lattice = VORONOI_LATTICE;
-	crack_growth_crit = false;
 	save_data = false;
 	save_cluster_dist = false;
 	use_voltage_boundaries = false;
@@ -45,7 +44,6 @@ int main(int argc, char *argv[]) {
 	save_damage = false;
 	save_damage_field = false;
 	save_conductivity = false;
-	save_toughness = false;
 	save_energy = true;
 	save_threshold = false;
 
@@ -59,7 +57,7 @@ int main(int argc, char *argv[]) {
 
 	// get commandline options
 
-	while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:gTER:")) != -1) {
+	while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDSvel:TER:")) != -1) {
 		switch (opt) {
 			case 'n':
 				N = atoi(optarg);
@@ -119,9 +117,6 @@ int main(int argc, char *argv[]) {
 						exit(EXIT_FAILURE);
 				}
 				break;
-			case 'g':
-				crack_growth_crit = true;
-				break;
 			case 'd':
 				save_damage = true;
 				break;
@@ -156,9 +151,6 @@ int main(int argc, char *argv[]) {
 			case 'r':
 				save_conductivity = true;
 				break;
-			case 't':
-				save_toughness = true;
-				break;
 			case 'E':
 				save_energy = true;
 				break;
@@ -222,9 +214,9 @@ int main(int argc, char *argv[]) {
 		}
 	}
 
-	double *crit_stress;
+	long double *crit_stress;
 	if (save_crit_stress) {
-		crit_stress = (double *)malloc(N * sizeof(double));
+		crit_stress = (long double *)malloc(N * sizeof(long double));
 	}
 
 	double *stress_field;
@@ -268,19 +260,14 @@ int main(int argc, char *argv[]) {
 		}
 	}
 
-	double *toughness;
-	if (save_toughness) {
-		toughness = (double *)malloc(N * sizeof(double));
-	}
-
-	double *energy;
+	long double *energy;
 	if (save_energy) {
-		energy = (double *)malloc(N * sizeof(double));
+		energy = (long double *)malloc(N * sizeof(long double));
 	}
 
-	double *thresholds;
+	long double *thresholds;
 	if (save_threshold) {
-		thresholds = (double *)malloc(N * sizeof(double));
+		thresholds = (long double *)malloc(N * sizeof(long double));
 	}
 
 
@@ -311,7 +298,7 @@ int main(int argc, char *argv[]) {
 		net_free(tmp_net, &c);
 
 		uint_t max_pos = 0;
-		double max_val = 0;
+		long double max_val = 0;
 
 		double cond0;
 		{
@@ -321,7 +308,7 @@ int main(int argc, char *argv[]) {
 		}
 
 		for (uint_t j = 0; j < data->num_broken; j++) {
-			double val = data->extern_field[j];
+			long double val = data->extern_field[j];
 
 			if (val > max_val) {
 				max_pos = j;
@@ -330,41 +317,15 @@ int main(int argc, char *argv[]) {
 		}
 
 		uint_t av_size = 0;
-		double cur_val = 0;
-		uint_t notch_mark = 0;
-
-		if (crack_growth_crit) {
-			uint_t n_edge = UINT_MAX;
-			for (uint_t j = 0; j < g->ne; j++) {
-				// right now the only fuses broken are from a notch
-				if (net->fuses[j]) {
-					n_edge = j;
-					break;
-				}
-			}
-
-			if (n_edge < UINT_MAX) {
-				uint_t n_vert = g->dev[2 * n_edge];
-				notch_mark = net->dual_marks[n_vert];
-			}
-		}
+		long double cur_val = 0;
 
 		for (uint_t j = 0; j < max_pos; j++) {
 			uint_t next_broken = data->break_list[j];
 
-			bool attached_to_notch = net->dual_marks[g->dev[2*next_broken]] == notch_mark;
-			bool grew_forward = g->ex[2 * next_broken] > crack_len;
-
-			if (crack_growth_crit && attached_to_notch && grew_forward) {
-				printf("whoops\n\n");
-				max_pos = j;
-				break;
-			}
-
 			bool refactor = ((j + 1) % refactor_every) == 0;
 			break_edge(net, next_broken, &c, refactor);
 
-			double val = data->extern_field[j];
+			long double val = data->extern_field[j];
 			if (save_cluster_dist) {
 				if (val < cur_val) {
 					av_size++;
@@ -422,12 +383,12 @@ int main(int argc, char *argv[]) {
 		}
 
 		if (save_energy) {
-			double tmp_energy = 0;
+			long double tmp_energy = 0;
 			if (max_pos > 0) {
-				double sigma1 = data->extern_field[0];
+				long double sigma1 = data->extern_field[0];
 				double cond1 = cond0;
 				for (uint_t j = 0; j < max_pos - 1; j++) {
-					double sigma2 = data->extern_field[j+1];
+					long double sigma2 = data->extern_field[j+1];
 					double cond2 = data->conductivity[j];
 					if (sigma2 > sigma1) {
 						tmp_energy += 0.5 * gsl_pow_2(sigma1) * (1 - cond2 / cond1) / cond1;
@@ -439,23 +400,6 @@ int main(int argc, char *argv[]) {
 			energy[i] = tmp_energy;
 		}
 
-		if (save_toughness) {
-			double tmp_toughness = 0;
-			if (max_pos > 0) {
-				double sigma1 = data->extern_field[0];
-				double epsilon1 = sigma1 / data->conductivity[0];
-				for (uint_t j = 0; j < max_pos - 1; j++) {
-					double sigma2 = data->extern_field[j+1];
-					double epsilon2 = sigma2 / data->conductivity[j+1];
-					if (epsilon2 > epsilon1) {
-						tmp_toughness += (sigma1 + sigma2) * (epsilon2 - epsilon1) / 2;
-						sigma1 = sigma2; epsilon1 = epsilon2;
-					}
-				}
-			}
-			toughness[i] = tmp_toughness;
-		}
-
 		if (save_threshold) {
 			thresholds[i] = net->thres[data->break_list[max_pos]];
 		}
@@ -499,7 +443,7 @@ int main(int argc, char *argv[]) {
 
 		if (save_data) {
 			for (uint_t j = 0; j < data->num_broken; j++) {
-				fprintf(data_out, "%u %g %g ", data->break_list[j],
+				fprintf(data_out, "%u %Lg %g ", data->break_list[j],
 								data->extern_field[j], data->conductivity[j]);
 			}
 			fprintf(data_out, "\n");
@@ -566,21 +510,11 @@ int main(int argc, char *argv[]) {
 		free(conductivity);
 	}
 
-	if (save_toughness) {
-		char *tough_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(tough_filename, filename_len, "tuff_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-		FILE *tough_file = fopen(tough_filename, "ab");
-		fwrite(toughness, sizeof(double), N, tough_file);
-		fclose(tough_file);
-		free(tough_filename);
-		free(toughness);
-	}
-
 	if (save_energy) {
 		char *tough_filename = (char *)malloc(filename_len * sizeof(char));
 		snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
 		FILE *tough_file = fopen(tough_filename, "ab");
-		fwrite(energy, sizeof(double), N, tough_file);
+		fwrite(energy, sizeof(long double), N, tough_file);
 		fclose(tough_file);
 		free(tough_filename);
 		free(energy);
@@ -590,7 +524,7 @@ int main(int argc, char *argv[]) {
 		char *thres_filename = (char *)malloc(filename_len * sizeof(char));
 		snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
 		FILE *thres_file = fopen(thres_filename, "ab");
-		fwrite(thresholds, sizeof(double), N, thres_file);
+		fwrite(thresholds, sizeof(long double), N, thres_file);
 		fclose(thres_file);
 		free(thres_filename);
 		free(thresholds);
@@ -612,7 +546,7 @@ int main(int argc, char *argv[]) {
 		char *str_filename = (char *)malloc(filename_len * sizeof(char));
 		snprintf(str_filename, filename_len, "strs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
 		FILE *str_file = fopen(str_filename, "ab");
-		fwrite(crit_stress, sizeof(double), N, str_file);
+		fwrite(crit_stress, sizeof(long double), N, str_file);
 		fclose(str_file);
 		free(str_filename);
 		free(crit_stress);
diff --git a/src/fracture.h b/src/fracture.h
index f122869..556e1e5 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -89,7 +89,7 @@ typedef struct {
 	uint_t num_broken;
 	uint_t *break_list;
 	double *conductivity;
-	double *extern_field;
+	long double *extern_field;
 } data_t;
 
 intptr_t *run_voronoi(uint_t num_coords, double *coords, bool periodic, double xmin, double xmax, double ymin, double ymax);
@@ -173,7 +173,7 @@ cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, cholmod_co
 
 data_t *data_create(uint_t num_edges);
 void data_free(data_t *data);
-void data_update(data_t *data, uint_t last_broke, double strength, double conductivity);
+void data_update(data_t *data, uint_t last_broke, long double strength, double conductivity);
 
 graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c);
 
diff --git a/src/long_anal_process.c b/src/long_anal_process.c
new file mode 100644
index 0000000..ba29152
--- /dev/null
+++ b/src/long_anal_process.c
@@ -0,0 +1,156 @@
+
+#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>
+#include <sys/stat.h>
+
+void get_mean(uint_t len, long double *data, long double result[2]) {
+	long double total = 0;
+
+	for (uint_t i = 0; i < len; i++) {
+		total += data[i];
+	}
+
+	long double mean = total / len;
+	long double total_sq_diff = 0;
+
+	for (uint_t i = 0; i < len; i++) {
+		total_sq_diff += pow(mean - data[i], 2);
+	}
+
+	long double mean_err = sqrt(total_sq_diff) / len;
+
+	result[0] = mean;
+	result[1] = mean_err;
+}
+
+void get_mom(uint_t len, uint_t n, long double *data, long double mean[2], long double result[2]) {
+	long double total_n_diff = 0;
+	long double total_n2_diff = 0;
+
+	for (uint_t i = 0; i < len; i++) {
+		total_n_diff += pow(fabsl(mean[0] - data[i]), n);
+		total_n2_diff += pow(fabsl(mean[0] - data[i]), 2 * n);
+	}
+
+	long double mom = total_n_diff / len;
+	long double mom_err = sqrt(total_n2_diff) / len;
+
+	result[0] = mom;
+	result[1] = mom_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));
+	long double *vals_c1 = (long double *)malloc(nc * sizeof(long double));
+	long double *errors_c1 = (long double *)malloc(nc * sizeof(long double));
+	long double *vals_c2 = (long double *)malloc(nc * sizeof(long double));
+	long double *errors_c2 = (long double *)malloc(nc * sizeof(long double));
+	long double *vals_c3 = (long double *)malloc(nc * sizeof(long double));
+	long double *errors_c3 = (long double *)malloc(nc * sizeof(long double));
+	long double *vals_c4 = (long double *)malloc(nc * sizeof(long double));
+	long double *errors_c4 = (long double *)malloc(nc * sizeof(long 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);
+
+		struct stat info;
+		stat(fn, &info);
+		uint_t num_bytes = info.st_size;
+		uint_t dist_len = (num_bytes * sizeof(char)) / sizeof(long double);
+
+		long double *dist = malloc(dist_len * sizeof(long double));
+		FILE *dist_file = fopen(fn, "rb");
+		fread(dist, sizeof(long double), dist_len, dist_file);
+		fclose(dist_file);
+		{
+			long double mom1[2];
+			get_mean(dist_len, dist, mom1);
+			vals_c1[i] = mom1[0];
+			errors_c1[i] = mom1[1];
+
+			long double mom2[2];
+			get_mom(dist_len, 2, dist, mom1, mom2);
+			vals_c2[i] = mom2[0];
+			errors_c2[i] = mom2[1];
+
+			long double mom3[2];
+			get_mom(dist_len, 3, dist, mom1, mom3);
+			vals_c3[i] = mom3[0];
+			errors_c3[i] = mom3[1];
+
+			long 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 %Lg %Lg %Lg %Lg %Lg %Lg %Lg %Lg\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);
+	free(vals_c4);
+	free(errors_c4);
+
+	return 0;
+}
+
diff --git a/src/net_fracture.c b/src/net_fracture.c
index 3bfcd83..b5ab524 100644
--- a/src/net_fracture.c
+++ b/src/net_fracture.c
@@ -46,7 +46,7 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff, uint_t refact
 
 		uint_t last_broke = get_next_broken(net, currents, cutoff);
 
-		double sim_current;
+		long double sim_current;
 
 		if (net->voltage_bound) {
 			sim_current = conductivity;
@@ -54,7 +54,7 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff, uint_t refact
 			sim_current = 1;
 		}
 
-		data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] / currents[last_broke]), conductivity);
+		data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] / ((long double)currents[last_broke])), conductivity);
 
 		free(voltages);
 		free(currents);