#include "fracture.h"

bool set_connected(const cholmod_sparse *laplacian, uint_t *marks,
									 int vertex, int label, int stop_at, int exclude) {
	unsigned int num_verts = (int_t)laplacian->nrow;

	int_t *ia = (int_t *)laplacian->p;
	int_t *ja = (int_t *)laplacian->i;
	double *a = (double *)laplacian->x;

	bool stopped = false;

	int *queue = (int *)malloc((num_verts) * sizeof(int));
	assert(queue != NULL);

	int queue_pos = 0;

	int *component_list = (int *)malloc((num_verts) * sizeof(int));
	assert(component_list != NULL);

	int list_pos = 0;
	queue[0] = vertex;
	component_list[0] = vertex;

	unsigned int old_label = marks[vertex];
	marks[vertex] = label;

	while (queue_pos >= 0) {
		int next_vertex = queue[queue_pos];
		queue_pos--;
		for (int j = 0; j < ia[next_vertex + 1] - ia[next_vertex]; j++) {
			if (ja[ia[next_vertex] + j] == stop_at && a[ia[next_vertex] + j] != 0) {
				// if we run into the other vertex, the components haven't changed
				stopped = true; queue_pos = -1;
				break;
			}
			if (marks[ja[ia[next_vertex] + j]] != label &&
					a[ia[next_vertex] + j] != 0 &&
					ja[ia[next_vertex] + j] < num_verts - exclude) {
				marks[ja[ia[next_vertex] + j]] = label;
				queue_pos++;
				queue[queue_pos] = ja[ia[next_vertex] + j];
				list_pos++;
				component_list[list_pos] = ja[ia[next_vertex] + j];
			}
		}
	}

	if (stopped) {
		for (unsigned int i = 0; i <= list_pos; i++) {
			marks[component_list[i]] = old_label;
		}
	}

	free(queue);
	free(component_list);

	if (stopped)
		return false;
	else
		return true;
}

int update_components(const cholmod_sparse *laplacian, unsigned int *marks,
											int old_num_components, int v1, int v2, int exclude) {
	assert(laplacian != NULL);
	assert(marks != NULL);

	if (set_connected(laplacian, marks, v1, old_num_components + 1, v2,
										exclude)) {
		return old_num_components + 1;
	} else
		return old_num_components;
}

unsigned int *find_components(const cholmod_sparse *laplacian, unsigned int skip) {
	size_t num_verts = laplacian->nrow;
	unsigned int *marks = (unsigned int *)calloc(num_verts, sizeof(unsigned int));

	int num_components = 0;

	for (int i = 0; i < num_verts; i++) {
		if (marks[i] == 0) {
			num_components++;
			set_connected(laplacian, marks, i, num_components, -1, skip);
		}
	}

	return marks;
}

uint_t find_cycles(uint_t num_edges, const bool *fuses, const uint_t *ev, const unsigned
		int *vei, const uint_t *ve, int **cycles) {
	uint_t num_cycles = 0;
	uint_t *elist = (uint_t *)malloc(num_edges * sizeof(uint_t));
	uint_t *side = (uint_t *)calloc(num_edges, sizeof(uint_t));
	uint_t *num = (uint_t *)malloc(num_edges * sizeof(uint_t));
	for (uint_t i = 0; i < num_edges; i++) {
		if (fuses[i]) {
		bool in_cycle = false;
		for (uint_t j = 0; j < num_cycles; j++) {
			uint_t k = 0;
			int e;
			while ((e = cycles[j][k]) >= 0) {
				if (e == i) {
					in_cycle = true;
					break;
				}
				k++;
			}
			if (in_cycle) {
				break;
			}
		}
		if (!in_cycle) {
			uint_t pos = 0;
			uint_t cur_e = i;
			bool *included = (bool *)calloc(num_edges, sizeof(bool));
			included[cur_e] = true;
			elist[0] = cur_e;
			side[0] = 1;
			num[0] = 0;
			while (true) {
				uint_t cv = ev[2 * cur_e + side[pos]];
				uint_t new_e = cur_e;
				uint_t j;
				for (j = num[pos]; j < vei[cv+1] - vei[cv]; j++) {
					new_e = ve[vei[cv] + j];
					if (new_e != cur_e && fuses[new_e]) break;
				}
				if (new_e != cur_e && fuses[new_e]) {
					if (new_e == elist[0] && cv == ev[2 * elist[0] + !side[0]]) {
						pos++;
						cycles[2*num_cycles] = (int *)malloc((pos + 1) * sizeof(int));
						cycles[2*num_cycles+1] = (int *)malloc((pos + 1) * sizeof(int));
						for (uint_t i = 0; i < pos; i++) {
							cycles[2*num_cycles][i] = elist[i];
							cycles[2*num_cycles+1][i] = side[i];
						}
						cycles[2*num_cycles][pos] = -1;
						cycles[2*num_cycles+1][pos] = -1;
						num_cycles++;
						pos--;
						num[pos]++;
						included[pos] = false;
						cur_e = elist[pos];
					} else if (included[new_e]) {
						if (pos == 0) break;
						included[cur_e] = false;
						pos--;
						num[pos]++;
						cur_e = elist[pos];
					} else {
						num[pos] = j;
						pos++;
						elist[pos] = new_e;
						included[new_e] = true;
						uint_t nv1 = ev[2 * new_e];
						uint_t nv2 = ev[2 * new_e + 1];
						uint_t v = nv1 == cv ? nv2 : nv1;
						side[pos] = v == nv1 ? 0 : 1;
						num[pos] = 0;
						cur_e = new_e;
					}
				} else {
					if (pos == 0) break;
					included[cur_e] = false;
					pos--;
					num[pos]++;
					cur_e = elist[pos];
				}
			}
			free(included);
		}
		}
	}

	free(elist);
	free(side);
	free(num);

	return num_cycles;
}