#include "fracture.h"

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

	// defining variables to be (potentially) set by command line flags
	uint8_t filename_len;
	uint32_t N;
	uint_t L;
	double beta, inf, cutoff, crack_len;
	bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
			 save_crit_stress, save_energy, save_conductivity,
			 save_damage, save_threshold, save_current_load;
	bound_t boundary;
	lattice_t lattice;


	// assume filenames will be less than 100 characters

	filename_len = 100;


	// set default values

	N = 100;
	L = 16;
	crack_len = 0.;
	beta = .3;
	inf = 1e10;
	cutoff = 1e-9;
	boundary = FREE_BOUND;
	lattice = VORONOI_LATTICE;
	save_data = false;
	save_cluster_dist = false;
	use_voltage_boundaries = false;
	use_dual = false;
	save_network = false;
	save_crit_stress = false;
	save_damage = false;
	save_conductivity = false;
	save_energy = false;
	save_threshold = false;
	save_current_load = false;


	uint8_t bound_i;
	char boundc2 = 'f';
	uint8_t lattice_i;
	char lattice_c = 'v';
	char dual_c = 'o';


	// get commandline options

	while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDl:TE")) != -1) {
		switch (opt) {
			case 'n':
				N = atoi(optarg);
				break;
			case 'L':
				L = atoi(optarg);
				break;
			case 'b':
				beta = atof(optarg);
				break;
			case 'l':
				crack_len = atof(optarg);
				break;
			case 'B':
				bound_i = atoi(optarg);
				switch (bound_i) {
					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;
					case 3:
						boundary = EMBEDDED_BOUND;
						boundc2 = 'e';
						use_dual = true;
						use_voltage_boundaries = true;
						break;
					default:
						printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), or 2 (TORUS_BOUND).\n");
						exit(EXIT_FAILURE);
				}
				break;
			case 'q':
				lattice_i = atoi(optarg);
				switch (lattice_i) {
					case 0:
						lattice = VORONOI_LATTICE;
						lattice_c = 'v';
						break;
					case 1:
						lattice = SQUARE_LATTICE;
						lattice_c = 's';
						break;
					case 2:
						lattice = VORONOI_HYPERUNIFORM_LATTICE;
						lattice_c = 'h';
						break;
					default:
						printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 (SQUARE_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n");
						exit(EXIT_FAILURE);
				}
				break;
			case 'd':
				save_damage = true;
				break;
			case 'V':
				use_voltage_boundaries = true;
				break;
			case 'D':
				use_dual = true;
				dual_c = 'd';
				break;
			case 'c':
				save_cluster_dist = true;
				break;
			case 'o':
				save_data = true;
				break;
			case 'N':
				save_network = true;
				break;
			case 's':
				save_crit_stress = true;
				break;
			case 'r':
				save_conductivity = true;
				break;
			case 'E':
				save_energy = true;
				break;
			case 'T':
				save_threshold = true;
				break;
			case 'C':
				save_current_load = true;
				break;
			default: /* '?' */
				exit(EXIT_FAILURE);
		}
	}


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

	FILE *data_out;
	if (save_data) {
		char *data_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(data_filename, filename_len, "data_%c_%c_%c_%c_%u_%g_%g.txt", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
		data_out = fopen(data_filename, "a");
		free(data_filename);
	}

	uint_t max_verts, max_edges;

	// these are very liberal estimates
	max_verts = 4 * pow(L, 2);
	max_edges = 4 * pow(L, 2);

	if (max_verts > CINT_MAX) {
		exit(EXIT_FAILURE);
	}

	// define arrays for saving cluster and avalanche distributions
	uint32_t *cluster_size_dist;
	uint32_t *avalanche_size_dist;
	char *c_filename;
	char *a_filename;
	if (save_cluster_dist) {
		cluster_size_dist =
				(uint32_t *)calloc(max_verts, sizeof(uint32_t));
		avalanche_size_dist =
				(uint32_t *)calloc(max_edges, sizeof(uint32_t));

		c_filename = (char *)malloc(filename_len * sizeof(char));
		a_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(c_filename, filename_len, "cstr_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
		snprintf(a_filename, filename_len, "avln_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);

		FILE *cluster_out = fopen(c_filename, "rb");
		FILE *avalanche_out = fopen(a_filename, "rb");

		if (cluster_out != NULL) {
			fread(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out);
			fclose(cluster_out);
		}
		if (avalanche_out != NULL) {
			fread(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out);
			fclose(avalanche_out);
		}
	}

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

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

	// define arrays for saving damage distributions
	uint32_t *damage;
	char *d_filename;
	if (save_damage) {
		damage =
				(uint32_t *)calloc(max_edges, sizeof(uint32_t));

		d_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(d_filename, filename_len, "damg_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);

		FILE *damage_out = fopen(d_filename, "rb");

		if (damage_out != NULL) {
			fread(damage, sizeof(uint32_t), max_edges, damage_out);
			fclose(damage_out);
		}
	}

	long double *energy;
	if (save_energy) {
		energy = (long double *)malloc(N * sizeof(long double));
	}

	long double *thresholds;
	if (save_threshold) {
		thresholds = (long double *)malloc(N * sizeof(long double));
	}

	long double *loads;
	if (save_current_load) {
		loads = (long double *)malloc(N * sizeof(long double));
	}


	// 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;
	}


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

		graph_t *g = graph_create(lattice, boundary, L, use_dual);
		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);
		net_free(tmp_net, &c);

		uint_t max_pos = 0;
		long double max_val = 0;

		double cond0;
		{
			double *tmp_voltages = net_voltages(net, &c);
			cond0 = net_conductivity(net, tmp_voltages);
			free(tmp_voltages);
		}

		for (uint_t j = 0; j < data->num_broken; j++) {
			long double val = data->extern_field[j];

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

		uint_t av_size = 0;
		long double cur_val = 0;

		for (uint_t j = 0; j < max_pos; j++) {
			uint_t next_broken = data->break_list[j];

			break_edge(net, next_broken, &c);

			long double val = 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_crit_stress) crit_stress[i] = data->extern_field[max_pos];

		if (save_conductivity) {
			if (max_pos > 0) {
				conductivity[i] = data->conductivity[max_pos - 1];
			} else {
				conductivity[i] = cond0;
			}
		}

		if (save_damage) {
			uint_t would_break = 0;
			double *tmp_voltage = net_voltages(net, &c);
			double *tmp_current = net_currents(net, tmp_voltage, &c);
			free(tmp_voltage);
			for (uint_t j = 0; j < g->ne; j++) {
				bool broken = net->fuses[j];
				bool under_thres = net->thres[j] < net->thres[data->break_list[max_pos]];
				bool zero_field = fabs(tmp_current[j]) < cutoff;
				if (!broken && under_thres && zero_field) {
					break_edge(net, j, &c);
				}
			}
			damage[net->num_broken]++;
			free(tmp_current);
		}

		if (save_energy) {
			long double tmp_energy = 0;
			if (max_pos > 0) {
				long double sigma1 = data->extern_field[0];
				double cond1 = cond0;
				for (uint_t j = 0; j < max_pos - 1; j++) {
					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;
						sigma1 = sigma2;
						cond1 = cond2;
					}
				}
			}
			energy[i] = tmp_energy;
		}

		if (save_threshold) {
			thresholds[i] = net->thres[data->break_list[max_pos]];
		}

		if (save_current_load) {
			loads[i] = data->extern_field[max_pos] / net->thres[data->break_list[max_pos]];
		}

		if (save_data) {
			for (uint_t j = 0; j < data->num_broken; j++) {
				fprintf(data_out, "%u %Lg %g ", data->break_list[j],
								data->extern_field[j], data->conductivity[j]);
			}
			fprintf(data_out, "\n");
		}

		data_free(data);
		if (save_network) {
			FILE *net_out = fopen("network.txt", "w");
			for (uint_t j = 0; j < g->nv; j++) {
				fprintf(net_out, "%f %f ", g->vx[2 * j],
						g->vx[2 * j + 1]);
			}
			fprintf(net_out, "\n");
			for (uint_t j = 0; j < g->ne; j++) {
				fprintf(net_out, "%u %u ", g->ev[2 * j], g->ev[2 * j + 1]);
			}
			fprintf(net_out, "\n");
			for (uint_t j = 0; j < g->dnv; j++) {
				fprintf(net_out, "%f %f ", g->dvx[2 * j],
						g->dvx[2 * j + 1]);
			}
			fprintf(net_out, "\n");
			for (uint_t j = 0; j < g->ne; j++) {
				fprintf(net_out, "%u %u ", g->dev[2 * j], g->dev[2 * j + 1]);
			}
			fprintf(net_out, "\n");
			for (uint_t j = 0; j < g->ne; j++) {
				fprintf(net_out, "%d ", net->fuses[j]);
			}
			fclose(net_out);
		}

		if (save_cluster_dist) {
			uint_t *tmp_cluster_dist = get_cluster_dist(net);
			for (uint_t j = 0; j < g->dnv; j++) {
				cluster_size_dist[j] += tmp_cluster_dist[j];
			}
			free(tmp_cluster_dist);
		}


		net_free(net, &c);
		graph_free(g);
	}

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

	if (save_cluster_dist) {
		FILE *cluster_out = fopen(c_filename, "wb");
		FILE *avalanche_out = fopen(a_filename, "wb");

		fwrite(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out);
		fwrite(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out);

		fclose(cluster_out);
		fclose(avalanche_out);

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

	if (save_conductivity) {
		char *cond_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(cond_filename, filename_len, "cond_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
		FILE *cond_file = fopen(cond_filename, "ab");
		fwrite(conductivity, sizeof(double), N, cond_file);
		fclose(cond_file);
		free(cond_filename);
		free(conductivity);
	}

	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(long double), N, tough_file);
		fclose(tough_file);
		free(tough_filename);
		free(energy);
	}

	if (save_threshold) {
		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(long double), N, thres_file);
		fclose(thres_file);
		free(thres_filename);
		free(thresholds);
	}

	if (save_current_load) {
		char *load_filename = (char *)malloc(filename_len * sizeof(char));
		snprintf(load_filename, filename_len, "load_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
		FILE *load_file = fopen(load_filename, "ab");
		fwrite(loads, sizeof(long double), N, load_file);
		fclose(load_file);
		free(load_filename);
		free(loads);
	}

	if (save_damage) {
		FILE *hdam_file = fopen(d_filename, "wb");
		fwrite(damage, sizeof(uint32_t), max_edges, hdam_file);
		fclose(hdam_file);
		free(d_filename);
		free(damage);
	}

	if (save_data) {
		fclose(data_out);
	}

	if (save_crit_stress) {
		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(long double), N, str_file);
		fclose(str_file);
		free(str_filename);
		free(crit_stress);
	}

	CHOL_F(finish)(&c);

	return 0;
}