#include "fracture.h"

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

	// defining variables to be (potentially) set by command line flags
	int num = 100;
	int width = 16;
	double beta = .3;
	bool save_clusters = false;
	bool voltage_bound = false;
	bool output_break_data = false;

	while ((opt = getopt(argc, argv, "n:w:b:cVo")) != -1) {
		switch (opt) {
		case 'n':
			num = atoi(optarg);
			break;
		case 'w':
			width = atoi(optarg);
			break;
		case 'b':
			beta = atof(optarg);
			break;
		case 'c':
			save_clusters = true;
			break;
		case 'V':
			voltage_bound = true;
			break;
		case 'o':
			output_break_data = true;
			break;
		default: /* '?' */
			exit(EXIT_FAILURE);
		}
	}

	FILE *break_out;
	if (output_break_data) {
		char *break_filename = (char *)malloc(100 * sizeof(char));
		snprintf(break_filename, 100, "breaks_%d_%.3f_%d.txt", width, beta, num);
		break_out = fopen(break_filename, "w");
		free(break_filename);
	}

	bool periodic = true;
	double inf = 1;
	double cutoff = 1e-10;

	// 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 = ini_square_network(width, periodic, false, &c);
	graph_t *network_p = ini_square_network(width / 2, periodic, false, &c);
	net_t *perm_instance = create_instance(network, inf, voltage_bound, true, &c);
	unsigned int c_dist_size = network->dnv;
	unsigned int c_p_dist_size = network_p->dnv;

	// define arrays for saving cluster and avalanche distributions
	unsigned int *cluster_size_dist =
			(unsigned int *)calloc(c_dist_size, sizeof(unsigned int));
	unsigned int *cluster_p_size_dist =
			(unsigned int *)calloc(c_p_dist_size, sizeof(unsigned int));

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

		data_t *breaking_data = NULL;
		net_t *instance = NULL;
		while (breaking_data == NULL) {
			double *fuse_thres = gen_fuse_thres(
					network->ne, network->edge_coords, beta, beta_scaling_flat);
			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 (val > 10 * min_val)
				break;
		}

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

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

		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);

		net_t *instance_p = coursegrain_square(tmp_instance, network_p, &c);

		unsigned int *tmp_cluster_p_dist = get_cluster_dist(instance_p, &c);
		for (unsigned int i = 0; i < network_p->dnv; i++) {
			cluster_p_size_dist[i] += tmp_cluster_p_dist[i];
		}
		free(tmp_cluster_p_dist);

		free_instance(tmp_instance, &c);
		free_instance(instance_p, &c);
		if (output_break_data) {
			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(breaking_data->break_list);
		free(breaking_data->extern_field);
		free(breaking_data);
	}

	printf("\033[F\033[JCURRENT_SCALING: COMPLETE");

	if (save_clusters) {
		FILE *cluster_out = get_file("cluster", width, 0, beta, 1, 1, num, false);
		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 < c_p_dist_size; i++) {
			fprintf(cluster_out, "%u ", cluster_p_size_dist[i]);
		}
		fclose(cluster_out);
	}

	if (output_break_data) {
		fclose(break_out);
	}

	free(cluster_size_dist);
	free(cluster_p_size_dist);
	free_instance(perm_instance, &c);
	free_net(network, &c);
	free_net(network_p, &c);

	CHOL_F(finish)(&c);

	return 0;
}