#include "fracture.h"

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

	// defining variables to be (potentially) set by command line flags
	unsigned int num, width, filename_len;
	double crack_len, crack_width, beta, inf, cutoff;
	boundary_type periodic;
	bool include_breaking, supplied_bound, save_clusters, voronoi,
			side_bounds, voltage_bound, repeat_voronoi, network_out, save_avg_stress, save_avg_ztress,
			save_app_stress, save_corr, save_cond, use_crit, use_first, save_damage, save_homo_damage;
	char *bound_filename;

	filename_len = 100;

	num = 100;
	width = 16;
	crack_len = 16;
	beta = .3;
	inf = 1e-8;
	cutoff = 1e-8;
	periodic = FREE_BOUND;
	include_breaking = false;
	supplied_bound = false;
	save_clusters = false;
	voronoi = false;
	side_bounds = false;
	voltage_bound = false;
	repeat_voronoi = false;
	network_out = false;
	save_avg_stress = false;
	save_avg_ztress = false;
	save_app_stress = false;
	save_damage = false;
	save_homo_damage = false;
	save_corr = false;
	save_cond = false;
	use_crit = true;
	use_first = false;

	int periodic_int;

	while ((opt = getopt(argc, argv, "n:w:b:l:f:ocp:vVsrNCaSFtzdH")) != -1) {
		switch (opt) {
		case 'n':
			num = atoi(optarg);
			break;
		case 'w':
			width = atoi(optarg);
			break;
		case 'b':
			beta = atof(optarg);
			break;
		case 'l':
			crack_len = atof(optarg);
			break;
		case 'p':
			periodic_int = atoi(optarg);
			switch (periodic_int) {
				case 0:
					periodic = FREE_BOUND;
					break;
				case 1:
					periodic = CYLINDER_BOUND;
					break;
				case 2:
					periodic = TORUS_BOUND;
					break;
			}
			break;
		case 'C':
			save_avg_stress = true;
			use_crit = true;
			break;
		case 'd':
			save_damage = true;
			break;
		case 'H':
			save_homo_damage = true;
			break;
		case 'z':
			save_avg_ztress = true;
			use_crit = true;
			break;
		case 'v':
			voronoi = true;
			break;
		case 'F':
			use_first = true;
			break;
		case 'r':
			repeat_voronoi = true;
			break;
		case 'V':
			voltage_bound = true;
			break;
		case 's':
			side_bounds = true;
			break;
		case 'c':
			save_clusters = true;
			break;
		case 'o':
			include_breaking = true;
			break;
		case 'N':
			network_out = true;
			break;
		case 'a':
			save_app_stress = true;
			break;
		case 'S':
			save_corr = true;
			break;
		case 't':
			save_cond = true;
			inf = 1;
			break;
		case 'f':
			supplied_bound = true;
			bound_filename = optarg;
			break;
		default: /* '?' */
			exit(EXIT_FAILURE);
		}
	}

	FILE *break_out;
	if (include_breaking) {
		char *break_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(break_filename, filename_len, "breaks_%u_%g_%g_%u.txt", width,
						 crack_len, beta, num);
		break_out = fopen(break_filename, "w");
		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 (voltage_bound) {
		(&c)->supernodal = CHOLMOD_SUPERNODAL;
	} else {
		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
	}

	graph_t *network;
	finst *perm_instance;
	unsigned int c_dist_size;
	unsigned int a_dist_size;
	unsigned int voronoi_num = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2);

	if (voronoi) {
		crack_width = 1. / (2 * width);
		c_dist_size = 2 * voronoi_num;
		a_dist_size = 2 * voronoi_num;
		if (repeat_voronoi) {
			while ((network = ini_voronoi_network(width, periodic,
																						genfunc_uniform, &c)) == NULL)
				;
			perm_instance = create_instance(network, inf, voltage_bound, false, &c);
			gen_crack(perm_instance, crack_len, &c);
			finish_instance(perm_instance, &c);
		}
	} else {
		network = ini_square_network(width, periodic, side_bounds, &c);
		crack_width = 1;
		perm_instance = create_instance(network, inf, voltage_bound, false, &c);
		gen_crack(perm_instance, crack_len, &c);
		finish_instance(perm_instance, &c);
		c_dist_size = network->dnv;
		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_clusters) {
		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, "cluster_%d_%g_%g.txt", width, crack_len,
						 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 *app_stress;
	double *app_stress_crit;
	if (save_app_stress) {
		app_stress = (double *)malloc(num * sizeof(double));
		app_stress_crit = (double *)malloc(num * sizeof(double));
	}

	double *avg_corr;
	unsigned int **dists = NULL;
	if (save_corr) {
		avg_corr = (double *)calloc(2 * voronoi_num, sizeof(double));
		if (!voronoi)
			dists = get_dists(network);
	}

	double *conductivity;
	if (save_cond) {
		conductivity = (double *)malloc(num * sizeof(double));
	}

	double *avg_stress;
	unsigned int *num_stress;
	if (save_avg_stress) {
		avg_stress = (double *)calloc(pow(width, 2), sizeof(double));
		num_stress = (unsigned int *)calloc(pow(width, 2), sizeof(unsigned int));
	}
	double *avg_ztress;
	if (save_avg_ztress) {
		avg_ztress = (double *)calloc(pow(width, 2), sizeof(double));
	}

	double *damage;
	if (save_damage) {
		damage = (double *)calloc(pow(width, 2), sizeof(double));
	}

	double *homo_damage;
	if (save_homo_damage) {
		homo_damage = (double *)calloc(num, sizeof(double));
	}

	double *square_bound;
	unsigned int square_num_verts = 2 * ((width + 1) / 2) * (width / 2 + 1);
	if (supplied_bound) {
		FILE *bound_file = fopen(bound_filename, "r");
		square_bound = (double *)calloc(square_num_verts, sizeof(double));
		for (unsigned int i = 0; i < square_num_verts; i++) {
			double tmp;
			fscanf(bound_file, "%lg ", &tmp);
			square_bound[i] = tmp;
		}
		fclose(bound_file);
		if (!voronoi) {
			double total = 0;
			for (unsigned int i = 0; i < square_num_verts; i++) {
				((double *)perm_instance->boundary_cond->x)[i] = square_bound[i];
				total += square_bound[i];
			}
			((double *)perm_instance->boundary_cond->x)[square_num_verts] = 0;
			((double *)perm_instance->boundary_cond->x)[square_num_verts + 1] = 0;
		}
	}

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

		data_t *breaking_data = NULL;
		while (breaking_data == NULL) {
			if (voronoi && !repeat_voronoi) {
				while ((network = ini_voronoi_network(width, periodic,
																							genfunc_uniform, &c)) == NULL)
					;
				perm_instance = create_instance(network, inf, voltage_bound, false, &c);
				gen_crack(perm_instance, crack_len, &c);
				finish_instance(perm_instance, &c);
				if (supplied_bound) {
					voronoi_bound_ini(perm_instance, square_bound, width);
				}
			}
			double *fuse_thres = gen_fuse_thres(
					network->ne, network->ex, beta, beta_scaling_flat);
			finst *instance = copy_instance(perm_instance, &c);
			breaking_data = fracture_network(instance, fuse_thres, &c, cutoff);
			free_instance(instance, &c);
			free(fuse_thres);
		}

		unsigned int min_pos = 0;
		double min_val = DBL_MAX;

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


			if (val < min_val) {
				min_pos = j;
				min_val = val;
			}
		}

		if (save_app_stress) {
			app_stress[DUMB] =
					1 / fabs(breaking_data->extern_field[breaking_data->num_broken - 1]);
			app_stress_crit[DUMB] = 1 / fabs(breaking_data->extern_field[min_pos]);
		}

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

		int stop_at = breaking_data->num_broken;
		if (use_crit)
			stop_at = min_pos;
		if (use_first)
			stop_at = 0;

		unsigned int av_size = 0;
		double cur_val = DBL_MAX;
		for (unsigned int i = 0; i < min_pos; i++) {
			double val = fabs(breaking_data->extern_field[i]);
			if (save_clusters) {
				if (val > cur_val) {
					av_size++;
				}

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

		for (unsigned int i = 0; i < stop_at; i++) {
			break_edge(tmp_instance, breaking_data->break_list[i], &c);
		}

		if (save_cond) {
			double *tmp_voltage = get_voltage(tmp_instance, &c);
			conductivity[DUMB] = fabs(tmp_voltage[tmp_instance->graph->nv + 1] - tmp_voltage[tmp_instance->graph->nv]);
			free(tmp_voltage);
		}

		if (save_homo_damage) {
			homo_damage[DUMB] = ((double)min_pos) / tmp_instance->graph->ne;
		}

		if (save_avg_stress || save_avg_ztress) {
			double *tmp_stress = get_current(tmp_instance, &c);
			if (voronoi) {
				double *tmp_stress_2 =
						bin_values(tmp_instance->graph, width, tmp_stress);
				free(tmp_stress);
				tmp_stress = tmp_stress_2;
			}
			for (unsigned int i = 0; i < pow(width, 2); i++) {
				double sigma = tmp_stress[i];
				if (sigma != 0 && sigma == sigma && save_avg_stress) {
					avg_stress[i] += sigma;
					num_stress[i]++;
				}
				if (sigma == sigma && save_avg_ztress) {
					avg_ztress[i] += sigma;
				}
			}
			free(tmp_stress);
		}

		if (save_damage) {
			double *tmp_damage = (double *)calloc(tmp_instance->graph->ne, sizeof(double));
			for (unsigned int i = 0; i < stop_at; i++) {
				tmp_damage[breaking_data->break_list[i]] += 1;
			}
			if (voronoi) {
				double *tmp_damage_2 = bin_values(tmp_instance->graph, width, tmp_damage);
				free(tmp_damage);
				tmp_damage = tmp_damage_2;
			}
			for (unsigned int i = 0; i < pow(width, 2); i++) {
				damage[i] += tmp_damage[i];
			}
			free(tmp_damage);
		}

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

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

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

		free_instance(tmp_instance, &c);
		if (voronoi && !repeat_voronoi) {
			free_net(network, &c);
			free_instance(perm_instance, &c);
		}
		if (include_breaking) {
			for (unsigned int i = 0; i < breaking_data->num_broken; i++) {
				fprintf(break_out, "%u %f ", breaking_data->break_list[i],
								breaking_data->extern_field[i]);
			}
			fprintf(break_out, "\n");
		}
		free_break_data(breaking_data);
	}

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

	if (save_clusters) {
		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(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_%d_%g_%g.txt", width, crack_len,
						 beta);
		FILE *corr_file = fopen(corr_filename, "w");
		for (unsigned int i = 0; i < 2 * voronoi_num; i++) {
			fprintf(corr_file, "%g ", avg_corr[i]);
		}
		fclose(corr_file);
		free(corr_filename);
	}

	if (save_cond) {
		char *cond_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(cond_filename, filename_len, "conductivity_%d_%g_%g.txt", width, crack_len, beta);
		FILE *cond_file = fopen(cond_filename, "a");
		for (unsigned int i = 0; i < num; i++) {
			fprintf(cond_file, "%g ", conductivity[i]);
		}
		fclose(cond_file);
		free(cond_filename);
		free(conductivity);
	}

	if (save_homo_damage) {
		char *hdam_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(hdam_filename, filename_len, "damagep_%d_%g.txt", width, beta);
		FILE *hdam_file = fopen(hdam_filename, "a");
		for (unsigned int i = 0; i < num; i++) {
			fprintf(hdam_file, "%g ", homo_damage[i]);
		}
		fclose(hdam_file);
		free(hdam_filename);
		free(homo_damage);
	}

	if (!voronoi || repeat_voronoi) {
		free_instance(perm_instance, &c);
		//free_net(network, &c);
	}

	if (include_breaking) {
		fclose(break_out);
	}

	if (save_avg_stress) {
		char *stress_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(stress_filename, filename_len, "current_%d_%g_%g_%u.txt", width,
						 crack_len, beta, num);
		FILE *stress_file = fopen(stress_filename, "w");
		for (unsigned int i = 0; i < pow(width, 2); i++) {
			if (num_stress[i] != 0)
				avg_stress[i] /= num_stress[i];
			fprintf(stress_file, "%g ", avg_stress[i]);
		}
		fclose(stress_file);
		free(stress_filename);
		free(avg_stress);
	}
	if (save_avg_ztress) {
		char *ztress_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(ztress_filename, filename_len, "zurrent_%d_%g_%g_%u.txt", width,
						 crack_len, beta, num);
		FILE *stress_file = fopen(ztress_filename, "w");
		for (unsigned int i = 0; i < pow(width, 2); i++) {
			fprintf(stress_file, "%g ", avg_ztress[i] / num);
		}
		fclose(stress_file);
		free(ztress_filename);
		free(avg_ztress);
	}

	if (save_app_stress) {
		char *a_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(a_filename, filename_len, "astress_%d_%g_%g.txt", width, crack_len,
						 beta);
		FILE *a_file = fopen(a_filename, "a");
		for (int i = 0; i < num; i++) {
			fprintf(a_file, "%g %g\n", app_stress[i], app_stress_crit[i]);
		}
		fclose(a_file);
		free(a_filename);
		free(app_stress);
		free(app_stress_crit);
	}

	if (save_damage) {
		char *damage_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(damage_filename, filename_len, "damage_%d_%g_%g_%u.txt", width,
						 crack_len, beta, num);
		FILE *damage_file = fopen(damage_filename, "w");
		for (unsigned int i = 0; i < pow(width, 2); i++) {
			fprintf(damage_file, "%g ", damage[i] / num);
		}
		fclose(damage_file);
		free(damage_filename);
		free(damage);

	}

	CHOL_F(finish)(&c);

	return 0;
}