#include "fracture.h"

int main(int argc, char *argv[]) {
	int opt;

	// defining variables to be (potentially) set by command line flags
	unsigned int N, L, filename_len;
	double beta, inf, cutoff;
	bound_t boundary;
	bool include_breaking, save_cluster_dist, use_voltage_boundaries, save_network,
			 save_crit_stress, save_toughness, save_corr, save_conductivity,
			 save_damage;

	filename_len = 100;

	N = 100;
	L = 16;
	beta = .3;
	inf = 1e-8;
	cutoff = 1e-10;
	boundary = FREE_BOUND;
	include_breaking = false;
	save_cluster_dist = false;
	use_voltage_boundaries = false;
	save_network = false;
	save_crit_stress = false;
	save_damage = false;
	save_corr = false;
	save_conductivity = false;
	save_toughness = false;

	int boundary_int;
	char boundc2 = 'f';

	while ((opt = getopt(argc, argv, "n:L:b:B:dVcoNsCrt")) != -1) {
		switch (opt) {
			case 'n':
				N = atoi(optarg);
				break;
			case 'L':
				L = atoi(optarg);
				break;
			case 'b':
				beta = atof(optarg);
				break;
			case 'B':
				boundary_int = atoi(optarg);
				switch (boundary_int) {
					case 0:
						boundary = FREE_BOUND;
						boundc2 = 'f';
						break;
					case 1:
						boundary = CYLINDER_BOUND;
						boundc2 = 'c';
						break;
					case 2:
						boundary = TORUS_BOUND;
						use_voltage_boundaries = true;
						boundc2 = 't';
						break;
					default:
						printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), or 2 (TORUS_BOUND).\n");
				}
				break;
			case 'd':
				save_damage = true;
				break;
			case 'V':
				use_voltage_boundaries = true;
				break;
			case 'c':
				save_cluster_dist = true;
				break;
			case 'o':
				include_breaking = true;
				break;
			case 'N':
				save_network = true;
				break;
			case 's':
				save_crit_stress = true;
				break;
			case 'C':
				save_corr = true;
				break;
			case 'r':
				save_conductivity = true;
				inf = 1;
				break;
			case 't':
				save_toughness = true;
				break;
			default: /* '?' */
				exit(EXIT_FAILURE);
		}
	}

	char boundc;
	if (use_voltage_boundaries) boundc = 'v';
	else boundc = 'c';

	FILE *break_out;
	if (include_breaking) {
		char *break_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(break_filename, filename_len, "breaks_s_%c_%c_%u_%g.txt", boundc, boundc2, L, beta);
		break_out = fopen(break_filename, "a");
		free(break_filename);
	}

	// start cholmod
	cholmod_common c;
	CHOL_F(start)(&c);

	/* if we use voltage boundary conditions, the laplacian matrix is positive
	 * definite and we can use a supernodal LL decomposition.  otherwise we need
	 * to use the simplicial LDL decomposition
	 */
	if (use_voltage_boundaries) {
		//(&c)->supernodal = CHOLMOD_SUPERNODAL;
		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
	} else {
		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
	}


	graph_t *network = ini_square_network(L, boundary, false, &c);
	net_t *perm_instance = create_instance(network, inf, use_voltage_boundaries, true, &c);
	unsigned int c_dist_size = network->dnv;
	unsigned int a_dist_size = network->nv;

	// define arrays for saving cluster and avalanche distributions
	unsigned int *cluster_size_dist;
	unsigned int *avalanche_size_dist;
	char *c_filename;
	if (save_cluster_dist) {
		cluster_size_dist =
				(unsigned int *)calloc(c_dist_size, sizeof(unsigned int));
		avalanche_size_dist =
				(unsigned int *)calloc(a_dist_size, sizeof(unsigned int));

		c_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(c_filename, filename_len, "cstr_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);

		FILE *cluster_out = fopen(c_filename, "r");

		if (cluster_out != NULL) {
			for (unsigned int i = 0; i < c_dist_size; i++) {
				unsigned int tmp;
				fscanf(cluster_out, "%u ", &tmp);
				cluster_size_dist[i] = tmp;
			}
			fscanf(cluster_out, "\n");
			for (unsigned int i = 0; i < a_dist_size; i++) {
				unsigned int tmp;
				fscanf(cluster_out, "%u ", &tmp);
				avalanche_size_dist[i] = tmp;
			}
			fclose(cluster_out);
		}
	}

	double *crit_stress;
	if (save_crit_stress) {
		crit_stress = (double *)malloc(N * sizeof(double));
	}

	double *avg_corr;
	unsigned int **dists = NULL;
	if (save_corr) {
		avg_corr = (double *)calloc(c_dist_size, sizeof(double));
	}

	double *conductivity;
	if (save_conductivity) {
		conductivity = (double *)malloc(N * sizeof(double));
	}

	double *damage;
	if (save_damage) {
		damage = (double *)malloc(N * sizeof(double));
	}

	double *toughness;
	if (save_toughness) {
		toughness = (double *)malloc(N * sizeof(double));
	}


	printf("\n");
	for (unsigned int i = 0; i < N; i++) {
		printf("\033[F\033[JFRACTURE: %0*d / %d\n", (int)log10(N) + 1, i + 1, N);

		double *fuse_thres = gen_fuse_thres(network->ne, network->edge_coords, beta, beta_scaling_flat);
		net_t *instance = copy_instance(perm_instance, &c);
		data_t *breaking_data = fracture_network(instance, fuse_thres, &c, cutoff);
		free_instance(instance, &c);
		free(fuse_thres);

		unsigned int max_pos = 0;
		double max_val = 0;

		for (unsigned int j = 0; j < breaking_data->num_broken; j++) {
			double val = breaking_data->extern_field[j];

			if (val > max_val) {
				max_pos = j;
				max_val = val;
			}
		}

		if (save_crit_stress) {
			crit_stress[i] =
					breaking_data->extern_field[max_pos];
		}

		net_t *tmp_instance = copy_instance(perm_instance, &c);

		unsigned int av_size = 0;
		double cur_val = DBL_MAX;
		for (unsigned int j = 0; j < max_pos; j++) {
			break_edge(tmp_instance, breaking_data->break_list[j], &c);

			double val = fabs(breaking_data->extern_field[j]);
			if (save_cluster_dist) {
				if (val < cur_val) {
					av_size++;
				}

				if (val > cur_val) {
					avalanche_size_dist[av_size]++;
					av_size = 0;
					cur_val = val;
				}
			}
		}

		if (save_conductivity) {
			if (!use_voltage_boundaries) {
				double *tmp_voltage = get_voltage(tmp_instance, &c);
				conductivity[i] = 1/fabs(tmp_voltage[tmp_instance->graph->nv + 1] - tmp_voltage[tmp_instance->graph->nv]);
				free(tmp_voltage);
			} else {
				conductivity[i] = breaking_data->conductivity[max_pos];
			}
		}

		if (save_toughness) {
			double tmp_toughness = 0;
			if (max_pos > 0) {
				double sigma1 = breaking_data->extern_field[0];
				double epsilon1 = sigma1 / breaking_data->conductivity[0];
				for (unsigned int j = 0; j < max_pos - 1; j++) {
					double sigma2 = breaking_data->extern_field[j+1];
					double epsilon2 = sigma2 / breaking_data->conductivity[j+1];
					if (epsilon2 > epsilon1) {
						tmp_toughness += (sigma1 + sigma2) * (epsilon2 - epsilon1) / 2;
						sigma1 = sigma2; epsilon1 = epsilon2;
					}
				}
			}
			toughness[i] = tmp_toughness;
		}

		if (save_damage) {
			damage[i] = ((double)max_pos) / tmp_instance->graph->ne;
		}

		if (save_cluster_dist) {
			unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c);
			for (unsigned int j = 0; j < tmp_instance->graph->dnv; j++) {
				cluster_size_dist[j] += tmp_cluster_dist[j];
			}
			free(tmp_cluster_dist);
		}

		if (save_corr) {
			double *tmp_corr = get_corr(tmp_instance, dists, &c);
			for (unsigned int j = 0; j < tmp_instance->graph->dnv; j++) {
				avg_corr[i] += tmp_corr[j] / N;
			}
			free(tmp_corr);
		}

		if (save_network) {
			FILE *net_out = fopen("network.txt", "w");
			for (unsigned int j = 0; j < network->nv; j++) {
				fprintf(net_out, "%f %f ", network->vert_coords[2 * j],
								tmp_instance->graph->vert_coords[2 * j + 1]);
			}
			fprintf(net_out, "\n");
			for (unsigned int j = 0; j < tmp_instance->graph->ne; j++) {
				fprintf(net_out, "%u %u ", tmp_instance->graph->ev[2 * j],
								tmp_instance->graph->ev[2 * j + 1]);
			}
			fprintf(net_out, "\n");
			for (unsigned int j = 0; j < tmp_instance->graph->dnv; j++) {
				fprintf(net_out, "%f %f ", tmp_instance->graph->dual_vert_coords[2 * j],
								tmp_instance->graph->dual_vert_coords[2 * j + 1]);
			}
			fprintf(net_out, "\n");
			for (unsigned int j = 0; j < tmp_instance->graph->ne; j++) {
				fprintf(net_out, "%u %u ", tmp_instance->graph->dev[2 * j],
								tmp_instance->graph->dev[2 * j + 1]);
			}
			fprintf(net_out, "\n");
			for (unsigned int j = 0; j < tmp_instance->graph->ne; j++) {
				fprintf(net_out, "%d ", tmp_instance->fuses[j]);
			}
			fclose(net_out);
		}

		free_instance(tmp_instance, &c);

		if (include_breaking) {
			for (unsigned int j = 0; j < breaking_data->num_broken; j++) {
				fprintf(break_out, "%u %g %g ", breaking_data->break_list[j],
								breaking_data->extern_field[j], breaking_data->conductivity[j]);
			}
			fprintf(break_out, "\n");
		}

		free_break_data(breaking_data);
	}

	printf("\033[F\033[JFRACTURE: COMPLETE\n");

	free_instance(perm_instance, &c);
	free_net(network, &c);

	if (save_cluster_dist) {
		FILE *cluster_out = fopen(c_filename, "w");

		for (int i = 0; i < c_dist_size; i++) {
			fprintf(cluster_out, "%u ", cluster_size_dist[i]);
		}
		fprintf(cluster_out, "\n");
		for (int i = 0; i < a_dist_size; i++) {
			fprintf(cluster_out, "%u ", avalanche_size_dist[i]);
		}
		fclose(cluster_out);

		free(c_filename);
		free(cluster_size_dist);
		free(avalanche_size_dist);
	}

	if (save_corr) {
		char *corr_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(corr_filename, filename_len, "corr_s_%c_%c_%d_%g.txt", boundc, boundc2, L,
						 beta);
		FILE *corr_file = fopen(corr_filename, "w");
		for (unsigned int i = 0; i < c_dist_size; i++) {
			fprintf(corr_file, "%g ", avg_corr[i]);
		}
		fclose(corr_file);
		free(corr_filename);
		free(avg_corr);
	}

	if (save_conductivity) {
		char *cond_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(cond_filename, filename_len, "cond_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
		FILE *cond_file = fopen(cond_filename, "a");
		for (unsigned int i = 0; i < N; i++) {
			fprintf(cond_file, "%g ", conductivity[i]);
		}
		fclose(cond_file);
		free(cond_filename);
		free(conductivity);
	}

	if (save_toughness) {
		char *tough_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(tough_filename, filename_len, "tuff_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
		FILE *tough_file = fopen(tough_filename, "a");
		for (unsigned int i = 0; i < N; i++) {
			fprintf(tough_file, "%g ", toughness[i]);
		}
		fclose(tough_file);
		free(tough_filename);
		free(toughness);
	}

	if (save_damage) {
		char *hdam_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(hdam_filename, filename_len, "damg_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
		FILE *hdam_file = fopen(hdam_filename, "a");
		for (unsigned int i = 0; i < N; i++) {
			fprintf(hdam_file, "%g ", damage[i]);
		}
		fclose(hdam_file);
		free(hdam_filename);
		free(damage);
	}

	if (include_breaking) {
		fclose(break_out);
	}

	if (save_crit_stress) {
		char *a_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(a_filename, filename_len, "strs_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
		FILE *a_file = fopen(a_filename, "a");
		for (int i = 0; i < N; i++) {
			fprintf(a_file, "%g ", crit_stress[i]);
		}
		fclose(a_file);
		free(a_filename);
		free(crit_stress);
	}

	CHOL_F(finish)(&c);

	return 0;
}