#include "fracture.h"

uint_t *get_spanning_edges(uint_t num_edges, uint_t *edges_to_verts, double *vert_coords, double cut, uint_t *n) {
	uint_t *spanning_edges = (uint_t *)malloc(num_edges * sizeof(uint_t));
	(*n) = 0;
	for (uint_t i = 0; i < num_edges; i++) {
		uint_t v1, v2;
		v1 = edges_to_verts[2 * i];
		v2 = edges_to_verts[2 * i + 1];
		double v1y, v2y;
		v1y = vert_coords[2 * v1 + 1];
		v2y = vert_coords[2 * v2 + 1];
		if ((fabs(v1y - v2y) < 0.5) && ((v1y <= cut && v2y > cut) || (v1y >= cut && v2y < cut))) {
			spanning_edges[*n] = i;
			(*n)++;
		}
	}
	return spanning_edges;
}

double *get_edge_coords(uint_t num_edges, double *vert_coords,
												uint_t *edges_to_verts) {
	double *output = (double *)malloc(2 * num_edges * sizeof(double));

	for (uint_t i = 0; i < num_edges; i++) {
		uint_t v1, v2;
		double v1x, v1y, v2x, v2y, dx, dy;
		v1 = edges_to_verts[2 * i];
		v2 = edges_to_verts[2 * i + 1];
		output[2 * i] = 0;
		output[2 * i + 1] = 0;
		v1x = vert_coords[2 * v1];
		v1y = vert_coords[2 * v1 + 1];
		v2x = vert_coords[2 * v2];
		v2y = vert_coords[2 * v2 + 1];
		dx = v1x - v2x;
		dy = v1y - v2y;
		if (fabs(dx) > 0.5) {
			if (dx > 0) {
				v2x = v2x + 1;
			} else {
				v1x = v1x + 1;
			}
		}
		if (fabs(dy) > 0.5) {
			if (dy > 0) {
				v2y = v2y + 1;
			} else {
				v1y = v1y + 1;
			}
		}
		output[2 * i] = (v1x + v2x) / 2;
		output[2 * i + 1] = (v1y + v2y) / 2;
	}

	return output;
}

uint_t *get_verts_to_edges_ind(uint_t num_verts,
																		 uint_t num_edges,
																		 const uint_t *edges_to_verts) {
	uint_t *output =
			(uint_t *)calloc(num_verts + 1, sizeof(uint_t));
	assert(output != NULL);

	for (uint_t i = 0; i < 2 * num_edges; i++) {
		if (edges_to_verts[i] < num_verts) {
			output[edges_to_verts[i] + 1]++;
		}
	}

	for (uint_t i = 0; i < num_verts; i++) {
		output[i + 1] += output[i];
	}

	return output;
}

uint_t *get_verts_to_edges(uint_t num_verts, uint_t num_edges,
																 const uint_t *edges_to_verts,
																 const uint_t *verts_to_edges_ind) {
	uint_t *output = (uint_t *)calloc(verts_to_edges_ind[num_verts],
																								sizeof(uint_t));
	uint_t *counts =
			(uint_t *)calloc(num_verts, sizeof(uint_t));
	for (int i = 0; i < 2 * num_edges; i++) {
		if (edges_to_verts[i] < num_verts) {
			output[verts_to_edges_ind[edges_to_verts[i]] +
						 counts[edges_to_verts[i]]] = i / 2;
			counts[edges_to_verts[i]]++;
		}
	}

	free(counts);

	return output;
}

uint_t get_cut_edges(uint_t ne, const uint_t *ev, const double *vx, bool both, uint_t *ce) {
	uint_t nce = 0;

	for (uint_t i = 0; i < ne; i++) {
		uint_t v1 = ev[2 * i];
		uint_t v2 = ev[2 * i + 1];

		double v1y = vx[2 * v1 + 1];
		double v2y = vx[2 * v2 + 1];

		if (fabs(v1y - v2y) > 0.5) {
			ce[nce] = i;
			nce++;
		} else if (both) {
			double v1x = vx[2 * v1];
			double v2x = vx[2 * v2];

			if (fabs(v1x - v2x) > 0.5) {
				ce[nce] = i;
				nce++;
			}
		}
	}

	return nce;
}

graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c) {
	graph_t *g = (graph_t *)calloc(1, sizeof(graph_t));
	frame_t *f = frame_create(lattice, L, dual);

	g->L = L;
	g->boundary = bound;
	g->ne = f->ne;

	if (bound == TORUS_BOUND) {
		g->nv = f->nv;
		g->dnv = f->dnv;
		g->nb = 1;

		g->ev = f->ev;
		f->ev = NULL;
		g->dev = f->dev;
		f->dev = NULL;
		g->vx = f->vx;
		f->vx = NULL;
		g->dvx = f->dvx;
		f->dvx = NULL;

		// the boundary for the torus consists of a list of edges required to cut
		// the torus into a cylinder
		g->bi = (uint_t *)calloc(2, sizeof(uint_t));
		g->b = (uint_t *)malloc(g->ne * sizeof(uint_t));
		g->bi[1] = get_cut_edges(g->ne, g->ev, g->vx, false, g->b);
		g->bq = (bool *)calloc(g->ne, sizeof(bool));
		for (uint_t i = 0; i < g->bi[1]; i++) {
			g->bq[g->b[i]] = true;
		}
	} else {
		uint_t *ce = (uint_t *)malloc(g->ne * sizeof(uint_t));
		uint_t nce = 0;

		if (bound == CYLINDER_BOUND) {
			g->nb = 2;
			nce = get_cut_edges(f->ne, f->ev, f->vx, false, ce);
		} else {
			g->nb = 4;
			nce = get_cut_edges(f->ne, f->ev, f->vx, true, ce);
		}

		g->nv = f->nv;
		g->dnv = f->dnv;
		g->vx = (double *)malloc(2 * (f->nv + nce) * sizeof(double));
		g->dvx = (double *)malloc(2 * (f->dnv + nce) * sizeof(double));
		g->ev = f->ev;
		g->dev = f->dev;
		f->ev = NULL;
		f->dev = NULL;
		memcpy(g->vx, f->vx, 2 * f->nv * sizeof(double));
		memcpy(g->dvx, f->dvx, 2 * f->dnv * sizeof(double));

		uint_t nbv = 0;
		uint_t *bv = (uint_t *)calloc(f->nv, sizeof(uint_t));
		uint_t *dbv = (uint_t *)calloc(f->dnv, sizeof(uint_t));
		uint_t nside = 0;
		bool *side = (bool *)calloc(f->nv, sizeof(bool));

		for (uint_t i = 0; i < nce; i++) {
			uint_t cv1 = g->ev[2 * ce[i]];
			uint_t cv2 = g->ev[2 * ce[i] + 1];
			uint_t dv1 = g->dev[2 * ce[i]];
			uint_t dv2 = g->dev[2 * ce[i] + 1];

			uint_t cin = 1;

			if (bound == FREE_BOUND && (f->vx[2 * cv2] - f->vx[2 * cv1]) > 0.5) {
				cin = 0; 
			}

			uint_t vin = f->vx[2 * cv1 + cin] < f->vx[2 * cv2 + cin] ? 0 : 1;
			uint_t dvin = f->dvx[2 * dv1 + cin] < f->dvx[2 * dv2 + cin] ? 0 : 1;

			if (bv[g->ev[2 * ce[i] + vin]] == 0) {
				nbv++;
				if (cin == 0) {
					nside++;
					side[g->ev[2 * ce[i] + vin]] = true;
				}

				bv[g->ev[2 * ce[i] + vin]] = g->nv;

				g->vx[2 * g->nv + cin] = 1 + f->vx[2 * g->ev[2 * ce[i] + vin] + cin];
				g->vx[2 * g->nv + !cin] = f->vx[2 * g->ev[2 * ce[i] + vin] + !cin];
				g->ev[2 * ce[i] + vin] = g->nv;

				g->nv++;
			} else {
				g->ev[2 * ce[i] + vin] = bv[g->ev[2 * ce[i] + vin]];
			}
			if (dbv[g->dev[2 * ce[i] + dvin]] == 0) {
				dbv[g->dev[2 * ce[i] + dvin]] = g->dnv;

				if (f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin] < 0.5) {
					g->dvx[2 * g->dnv + cin] = 1 + f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin];
				} else {
					g->dvx[2 * g->dnv + cin] = f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin];
				}
				g->dvx[2 * g->dnv + !cin] = f->dvx[2 * g->dev[2 * ce[i] + dvin] + !cin];
				g->dev[2 * ce[i] + dvin] = g->dnv;

				g->dnv++;
			} else {
				g->dev[2 * ce[i] + dvin] = dbv[g->dev[2 * ce[i] + dvin]];
			}
		}

		g->bi = (uint_t *)calloc(1 + g->nb, sizeof(uint_t));
		g->b = (uint_t *)malloc(2 * nbv * sizeof(uint_t));

		g->bi[1] = ((int_t)nbv - (int_t)nside);
		g->bi[g->nb] = 2 * nbv;

		if (bound == FREE_BOUND) {
			g->bi[2] = 2 * ((int_t)nbv - (int_t)nside);
			g->bi[3] = 2 * (int_t)nbv - (int_t)nside;
		}

		uint_t n_ins0 = 0;
		uint_t n_ins1 = 0;

		g->nbi = (uint_t *)malloc(((int_t)g->nv - (int_t)g->bi[g->nb]) * sizeof(uint_t));
		g->bni = (uint_t *)malloc((g->nv + 1) * sizeof(uint_t));
		g->bq = (bool *)calloc(g->nv, sizeof(bool));
		uint_t n_tmp = 0;

		for (uint_t i = 0; i < f->nv; i++) {
			g->bni[i] = n_tmp;
			if (bv[i] != 0) {
				g->bq[i] = true;
				g->bq[bv[i]] = true;
				if (side[i]) {
					g->b[g->bi[2] + n_ins1] = i;
					g->b[g->bi[3] + n_ins1] = bv[i];
					n_ins1++;
				} else {
					g->b[g->bi[0] + n_ins0] = i;
					g->b[g->bi[1] + n_ins0] = bv[i];
					n_ins0++;
				}
			} else {
				g->nbi[n_tmp] = i;
				n_tmp++;
			}
		}

		for (uint_t i = 0; i < g->nv - f->nv + 1; i++) {
			g->bni[f->nv + i] = n_tmp;
		}

		free(bv);
		free(dbv);
		free(side);
	}

	g->vei = get_verts_to_edges_ind(g->nv, g->ne, g->ev);
	g->ve = get_verts_to_edges(g->nv, g->ne, g->ev, g->vei);
	g->dvei = get_verts_to_edges_ind(g->dnv, g->ne, g->dev);
	g->dve = get_verts_to_edges(g->dnv, g->ne, g->dev, g->dvei);

	g->voltcurmat = gen_voltcurmat(g->ne, g->nv, g->ev, c);
	uint_t num_spanning_edges;
	g->spanning_edges = get_spanning_edges(g->ne, g->ev, g->vx, 0.5, &num_spanning_edges);
	g->num_spanning_edges = num_spanning_edges;

	frame_free(f);

	return g;
}