#include "fracture.h"

cholmod_sparse *gen_adjacency(const net_t *instance, bool dual, bool breakv,
															uint_t pad, cholmod_common *c) {
	uint_t ne = instance->graph->ne;
	if (!breakv && instance->graph->boundary != TORUS_BOUND && !dual) {
		ne += instance->graph->bound_inds[2];
	}

	uint_t nnz = 2 * ne;

	uint_t nv;
	if (dual)
		nv = instance->graph->dnv;
	else {
		if (breakv) nv = instance->graph->nv_break;
		else nv = instance->graph->nv;
		if (instance->graph->boundary != TORUS_BOUND && !breakv) {
			nv += 2;
		}
	}

	cholmod_triplet *t =
			CHOL_F(allocate_triplet)(nv + pad, nv + pad, nnz, 0, CHOLMOD_REAL, c);
	int_t *ri = (int_t *)t->i;
	int_t *ci = (int_t *)t->j;
	double *ai = (double *)t->x;

	t->nnz = nnz;

	uint_t *etv;
	if (dual)
		etv = instance->graph->dev;
	else {
		if (breakv) etv = instance->graph->ev_break;
		else etv = instance->graph->ev;
	}

	uint_t count = 0;

	for (uint_t i = 0; i < instance->graph->ne; i++) {
		if ((instance->fuses[i] && dual) || (!instance->fuses[i] && !dual)) {
			uint_t v1 = etv[2 * i];
			uint_t v2 = etv[2 * i + 1];

			ri[2 * count] = v1;
			ci[2 * count] = v2;
			ai[2 * count] = 1;

			ri[2 * count + 1] = v2;
			ci[2 * count + 1] = v1;
			ai[2 * count + 1] = 1;

			count++;
		} else {
			uint_t v1 = etv[2 * i];
			uint_t v2 = etv[2 * i + 1];

			ri[2 * count] = v1;
			ci[2 * count] = v2;
			ai[2 * count] = 0;

			ri[2 * count + 1] = v2;
			ci[2 * count + 1] = v1;
			ai[2 * count + 1] = 0;

			count++;
		}
	}

	if (!breakv && instance->graph->boundary != TORUS_BOUND && !dual) {
		for (uint_t i = 0; i < 2; i++) {
			for (uint_t j = 0; j < instance->graph->bound_inds[i+1] - instance->graph->bound_inds[i]; j++) {
				ri[2*count] = instance->graph->nv + i;
				ci[2*count] = instance->graph->bound_verts[instance->graph->bound_inds[i] + j];
				ai[2*count] = 1;
				ri[2*count+1] = instance->graph->bound_verts[instance->graph->bound_inds[i] + j];
				ci[2*count+1] = instance->graph->nv + i;
				ai[2*count+1] = 1;
				count++;
			}
		}
	}


	cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c);
	CHOL_F(free_triplet)(&t, c);

	return s;
}

cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c,
															bool symmetric) {
	const graph_t *network = instance->graph;
	uint_t num_verts = network->nv_break;
	uint_t num_bounds = network->num_bounds;
	double inf = instance->inf;
	bool voltage_bound = instance->voltage_bound;
	uint_t *bound_inds = network->bound_inds;
	uint_t *bound_verts = network->bound_verts;

	uint_t num_gedges;
	if (voltage_bound) {
		num_gedges = network->bound_inds[num_bounds];
	} else {
		num_gedges = network->bound_inds[num_bounds] + (num_bounds - 2) * 2;
	}
	if (network->boundary == TORUS_BOUND) num_gedges = bound_inds[1];
	uint_t num_gverts = network->break_dim;

	uint_t nnz = num_gverts + 2 * num_gedges;

	cholmod_triplet *temp_m =
			CHOL_F(allocate_triplet)(num_gverts, num_gverts, nnz, 0, CHOLMOD_REAL, c);

	int_t *rowind = (int_t *)temp_m->i;
	int_t *colind = (int_t *)temp_m->j;
	double *acoo = (double *)temp_m->x;

	temp_m->nnz = nnz;

	for (uint_t i = 0; i < num_gverts; i++) {
		rowind[i] = i;
		colind[i] = i;
		acoo[i] = 0;
	}

	cholmod_sparse *adjacency = gen_adjacency(instance, false, true, (int)num_gverts - (int)num_verts, c);
	int_t *ia = (int_t *)adjacency->p;
	double *aa = (double *)adjacency->x;

	for (uint_t i = 0; i < num_verts; i++) {
		for (uint_t j = 0; j < ia[i + 1] - ia[i]; j++) {
			acoo[i] += aa[ia[i] + j];
		}
	}

	if (network->boundary == TORUS_BOUND) {
		for (uint_t i = 0; i < network->bound_inds[1]; i++) {
			uint_t vv = network->bound_verts[i];
			rowind[num_gverts + 2*i] = vv;
			colind[num_gverts + 2*i] = network->nv + i;
			acoo[num_gverts + 2*i] = -1;
			rowind[num_gverts + 2*i+1] = network->nv + i;
			colind[num_gverts + 2*i+1] = vv;
			acoo[num_gverts + 2*i+1] = -1;
			acoo[vv]++;
			acoo[network->nv + i]++;
		}
		acoo[bound_verts[0]]++;
	}

	if (network->boundary != TORUS_BOUND) {
		if (voltage_bound) {
			for (uint_t i = 0; i < num_bounds; i++) {
				for (uint_t j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) {
					acoo[bound_verts[bound_inds[i] + j]]++;

					rowind[nnz - 1 - 2 * (bound_inds[i] + j)] = bound_verts[bound_inds[i] + j];
					colind[nnz - 1 - 2 * (bound_inds[i] + j)] = num_gverts - 1;
					acoo[nnz - 1 - 2 * (bound_inds[i] + j)] = -1;
					rowind[nnz - 1 - 2 * (bound_inds[i] + j) - 1] = num_gverts - 1;
					colind[nnz - 1 - 2 * (bound_inds[i] + j) - 1] = bound_verts[bound_inds[i] + j];
					acoo[nnz - 1 - 2 * (bound_inds[i] + j) - 1] = -1;
					acoo[num_gverts - 1]++;
					acoo[bound_verts[bound_inds[i] + j]]++;
				}
			}
		} else {
			for (uint_t i = 0; i < num_bounds; i++) {
				if (network->boundary != EMBEDDED_BOUND) acoo[0]++;

				rowind[num_verts + i] = num_verts + i;
				colind[num_verts + i] = num_verts + i;
				if (i < 2)
					acoo[num_verts + i] = (network->bound_inds[i + 1] -
																 network->bound_inds[i] + network->num_bounds - 2) *
																inf;
				else
					acoo[num_verts + i] =
							(network->bound_inds[i + 1] - network->bound_inds[i] + 2) * inf;
			}

			uint_t start = num_gverts;

			for (uint_t i = 0; i < 2; i++) {
				for (uint_t j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) {
					rowind[start + bound_inds[i] + j] = bound_verts[bound_inds[i] + j];
					colind[start + bound_inds[i] + j] = num_verts + i;
					acoo[start + bound_inds[i] + j] = -inf;
					acoo[bound_verts[bound_inds[i] + j]] += inf;
					colind[start + bound_inds[num_bounds] + bound_inds[i] + j] =
							bound_verts[bound_inds[i] + j];
					rowind[start + bound_inds[num_bounds] + bound_inds[i] + j] =
							num_verts + i;
					acoo[start + bound_inds[num_bounds] + bound_inds[i] + j] = -inf;
				}
			}

			for (uint_t i = 0; i < num_bounds - 2; i++) {
				rowind[nnz - 2 * i - 1] = num_verts;
				colind[nnz - 2 * i - 1] = num_verts + 2 + i;
				acoo[nnz - 2 * i - 1] = -inf;
				rowind[nnz - 2 * i - 2] = num_verts + 1;
				colind[nnz - 2 * i - 2] = num_verts + 2 + i;
				acoo[nnz - 2 * i - 2] = -inf;
				colind[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = num_verts;
				rowind[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = num_verts + 2 + i;
				acoo[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = -inf;
				colind[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = num_verts + 1;
				rowind[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = num_verts + 2 + i;
				acoo[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = -inf;
			}
		}
	}


	for (uint_t i = 0; i < num_gverts; i++) {
		if (instance->marks[i] != instance->marks[network->nv])
			acoo[i]++;
	}

	assert(CHOL_F(check_triplet)(temp_m, c));

	cholmod_sparse *t_out = CHOL_F(triplet_to_sparse)(temp_m, nnz, c);
	assert(CHOL_F(check_sparse)(t_out, c));

	double alpha[2] = {1, 0};
	double beta[2] = {-1, 0};
	cholmod_sparse *t_laplacian =
			CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c);

	cholmod_sparse *laplacian;
	if (symmetric)
		laplacian = CHOL_F(copy)(t_laplacian, 1, 1, c);
	else
		laplacian = CHOL_F(copy)(t_laplacian, 0, 1, c);

	CHOL_F(free_sparse)(&t_out, c);
	CHOL_F(free_sparse)(&adjacency, c);
	CHOL_F(free_sparse)(&t_laplacian, c);
	CHOL_F(free_triplet)(&temp_m, c);

	return laplacian;
}