#include "fracture.h"

uint_t get_next_broken(net_t *net, double *currents, double cutoff) {
	uint_t max_pos = UINT_MAX;
	double max_val = 0;

	for (uint_t i = 0; i < net->graph->ne; i++) {
		double current = fabs(currents[i]);
		bool broken = net->fuses[i];

		if (!broken && current > cutoff) {
			double val = current / net->thres[i];

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

	if (max_pos == UINT_MAX) {
		printf("GET_NEXT_BROKEN: currents is zero or NaN, no max_val found\n");
		exit(EXIT_FAILURE);
	}

	return max_pos;
}

data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
	data_t *data = data_create(net->graph->ne);

	while (true) {
		double *voltages = net_voltages(net, c);
		double *currents = net_currents(net, voltages, c);

		double conductivity = net_conductivity(net, voltages);

		if (conductivity < 1e-12 && net->graph->boundary == TORUS_BOUND) {
			free(voltages);
			free(currents);
			break;
		}

		uint_t last_broke = get_next_broken(net, currents, cutoff);

		double sim_current;

		if (net->voltage_bound) {
			sim_current = conductivity;
		} else {
			sim_current = 1;
		}

		data_update(data, last_broke, fabs(sim_current * (net->thres)[last_broke] / currents[last_broke]), conductivity);

		free(voltages);
		free(currents);

		break_edge(net, last_broke, c);

		if (net->num_components > 1 && net->graph->boundary == TORUS_BOUND) {
			break;
		}

		if (net->marks[net->graph->nv] != net->marks[net->graph->nv + 1] && net->graph->boundary != TORUS_BOUND) {
			break;
		}
	}

	return data;
}