#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;
}

graph_t *ini_square_network(uint_t width, bound_t boundary, bool side_bounds,
												 cholmod_common *c) {
	graph_t *network = (graph_t *)calloc(1, sizeof(graph_t));

	network->boundary = boundary;
	bool periodic = (boundary == CYLINDER_BOUND) || (boundary == TORUS_BOUND) ? true : false;

	network->ne = pow(width, 2);
	if (boundary == CYLINDER_BOUND) {
		assert(width % 2 == 0);
		assert(!side_bounds);
		network->nv = (width / 2) * (width + 1);
		network->nv_break = (width / 2) * (width + 1);
		network->dnv = (width / 2 + 1) * (width / 2) + pow(width / 2, 2);
		network->num_bounds = 2;
		network->bound_inds = (uint_t *)malloc((network->num_bounds + 1) *
																								 sizeof(uint_t));
		network->bound_inds[0] = 0;
		network->bound_inds[1] = width / 2;
		network->bound_inds[2] = width;
		network->bound_verts = (uint_t *)calloc(width, sizeof(uint_t));
		network->break_dim = network->nv + network->num_bounds;
	} else if (boundary == FREE_BOUND || boundary == EMBEDDED_BOUND) {
		network->nv = 2 * ((width + 1) / 2) * (width / 2 + 1);
		network->nv_break = 2 * ((width + 1) / 2) * (width / 2 + 1);
		network->dnv = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2);
		if (side_bounds)
			network->num_bounds = 4;
		else
			network->num_bounds = 2;
		network->bound_inds = (uint_t *)malloc((network->num_bounds + 1) *
																								 sizeof(uint_t));
		for (uint_t i = 0; i < network->num_bounds + 1; i++) {
			network->bound_inds[i] = i * ((width + 1) / 2);
		}
		network->bound_verts = (uint_t *)malloc(
				network->num_bounds * ((width + 1) / 2) * sizeof(uint_t));
		if (side_bounds) {
			for (uint_t i = 0; i < (width + 1) / 2; i++) {
				network->bound_verts[2 * ((width + 1) / 2) + i] =
						(width + 1) / 2 + i * (width + 1);
				if (width % 2) {
					network->bound_verts[3 * ((width + 1) / 2) + i] =
							(width + 1) / 2 + i * (width + 1) - 1;
				} else {
					network->bound_verts[3 * ((width + 1) / 2) + i] =
							(width + 1) / 2 + i * (width + 1) + width / 2;
				}
			}
		}
		network->break_dim = network->nv + 2;
	} else if (boundary == TORUS_BOUND) {
		network->nv = (width / 2) * (width + 1) - (width / 2);
		network->nv_break = (width / 2) * (width + 1);
		network->dnv = (width / 2 + 1) * (width / 2) + pow(width / 2, 2) - (width / 2);
		network->num_bounds = 1;
		network->bound_inds = (uint_t *)malloc((network->num_bounds + 1) *
																								 sizeof(uint_t));
		network->bound_inds[0] = 0;
		network->bound_inds[1] = width / 2;
		network->bound_verts = (uint_t *)calloc(width / 2, sizeof(uint_t));
		network->break_dim = network->nv_break;
	}
	if (boundary != TORUS_BOUND) {
		for (uint_t i = 0; i < (width + 1) / 2; i++) {
			network->bound_verts[i] = i;
			network->bound_verts[(width + 1) / 2 + i] = network->nv - (width + 1) / 2 + i;
		} 
	} else {
		for (uint_t i = 0; i < width / 2; i++) {
			network->bound_verts[i] = i;
		}
	}
	network->ev_break =
			(uint_t *)calloc(2 * network->ne, sizeof(uint_t));
	network->ev =
			(uint_t *)calloc(2 * network->ne, sizeof(uint_t));
	for (uint_t i = 0; i < network->ne; i++) {
		network->ev_break[2 * i] = edge_to_verts(width, periodic, i, 1);
		network->ev_break[2 * i + 1] = edge_to_verts(width, periodic, i, 0);
		network->ev[2 * i] = network->ev_break[2 * i] % network->nv;
		network->ev[2 * i + 1] = network->ev_break[2 * i + 1] % network->nv;
	}
	network->vei =
			(uint_t *)calloc(network->nv + 1, sizeof(uint_t));
	network->vei[0] = 0;
	uint_t pos1 = 0;
	for (uint_t i = 0; i < network->nv; i++) {
		bool in_bound = false;
		for (uint_t j = 0; j < network->num_bounds; j++) {
			for (uint_t k = 0;
					 k < network->bound_inds[j + 1] - network->bound_inds[j]; k++) {
				if (i == network->bound_verts[network->bound_inds[j] + k]) {
					in_bound = true;
					break;
				}
			}
		}
		if (in_bound)
			pos1 += 2;
		else
			pos1 += 4;

		network->vei[i + 1] = pos1;
	}
	network->ve = (uint_t *)calloc(
			network->vei[network->nv], sizeof(uint_t));
	uint_t *vert_counts =
			(uint_t *)calloc(network->nv, sizeof(uint_t));
	for (uint_t i = 0; i < network->ne; i++) {
		uint_t v0 = network->ev[2 * i];
		uint_t v1 = network->ev[2 * i + 1];
		network->ve[network->vei[v0] + vert_counts[v0]] =
				i;
		network->ve[network->vei[v1] + vert_counts[v1]] =
				i;
		vert_counts[v0]++;
		vert_counts[v1]++;
	}
	free(vert_counts);

	network->vx =
			(double *)malloc(2 * network->nv * sizeof(double));
	for (uint_t i = 0; i < network->nv; i++) {
		if (!periodic) {
		network->vx[2 * i + 1] = ((double)((2 * i + 1) / (width + 1)))/width;
		network->vx[2 * i] = ((double)((2 * i + 1) % (width + 1)))/width;
		} else {
		network->vx[2 * i + 1] = ((double)((2 * i + 1) / (width)))/width;
		network->vx[2 * i] = ((double)((((2 * i + 1) / (width)+1) % 2) + ((2 * i) % (width))))/width;
		}
	}

	network->L = width;
	network->ex = get_edge_coords(
			network->ne, network->vx, network->ev);

	network->dev =
			(uint_t *)malloc(2 * network->ne * sizeof(uint_t));
	for (uint_t i = 0; i < network->ne; i++) {
		network->dev[2 * i] =
				dual_edge_to_verts(width, periodic, i, 0) % network->dnv;
		network->dev[2 * i + 1] =
				dual_edge_to_verts(width, periodic, i, 1) % network->dnv;
	}
	network->dvx =
			(double *)malloc(2 * network->dnv * sizeof(double));
	for (uint_t i = 0; i < network->dnv; i++) {
		network->dvx[2 * i + 1] =
				dual_vert_to_coord(width, periodic, i, 0) / width;
		network->dvx[2 * i] =
				dual_vert_to_coord(width, periodic, i, 1) / width;
	}

	network->voltcurmat = gen_voltcurmat(network->ne,
																			 network->break_dim,
																			 network->ev_break, c);

	network->dvei =
			get_verts_to_edges_ind(network->dnv, network->ne,
														 network->dev);
	network->dve = get_verts_to_edges(
			network->dnv, network->ne, network->dev,
			network->dvei);
	network->spanning_edges = get_spanning_edges(network->ne, network->ev, network->vx, 0.51, &(network->num_spanning_edges));

	return network;
}

uint_t *get_voro_dual_edges(uint_t num_edges,
																	uint_t num_verts, uint_t *edges,
																	uint_t *triangles) {
	uint_t *dual_edges =
			(uint_t *)malloc(2 * num_edges * sizeof(uint_t));
	uint_t place = 0;
	for (uint_t i = 0; i < num_edges; i++) {
		uint_t v1, v2;
		v1 = edges[2 * i];
		v2 = edges[2 * i + 1];
		if (v1 < num_verts && v2 < num_verts) {
			bool found_match = false;
			for (uint_t j = 0; j < 3; j++) {
				for (uint_t k = 0; k < 3; k++) {
					uint_t t11, t12, t21, t22;
					t11 = triangles[3 * v1 + j];
					t12 = triangles[3 * v1 + ((j + 1) % 3)];
					t21 = triangles[3 * v2 + k];
					t22 = triangles[3 * v2 + ((k + 1) % 3)];
					if ((t11 == t21 && t12 == t22) || (t11 == t22 && t12 == t21)) {
						dual_edges[2 * place] = t11 < t12 ? t11 : t12;
						dual_edges[2 * place + 1] = t11 < t12 ? t12 : t11;
						place++;
						found_match = true;
						break;
					}
				}
				if (found_match)
					break;
			}
		}
	}

	return dual_edges;
}

graph_t *ini_voro_graph(uint_t L, bound_t boundary, bool use_dual,
													double *(*genfunc)(uint_t, bound_t, gsl_rng *, uint_t *),
													cholmod_common *c) {
	graph_t *g = (graph_t *)calloc(1, sizeof(graph_t));

	// generate the dual lattice
	double *lattice;
	uint_t num;
	{
		gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
		FILE *rf = fopen("/dev/urandom", "r");
		unsigned long int seed;
		fread(&seed, sizeof(unsigned long int), 1, rf);
		fclose(rf);
		gsl_rng_set(r, seed);
		lattice = genfunc(L, boundary, r, &num);
		gsl_rng_free(r);
	}

	// retrieve a periodic voronoi tesselation of the lattice
	bool run_periodic;
	if (boundary == EMBEDDED_BOUND) run_periodic = false;
	else run_periodic = true;
	intptr_t *vout = run_voronoi(num, lattice, run_periodic, 0, 1, 0, 1);

	uint_t tmp_num_verts = ((uint_t *)vout[0])[0];
	uint_t tmp_num_edges = ((uint_t *)vout[0])[1];
	double *tmp_vert_coords = (double *)vout[2];
	uint_t *tmp_edges = (uint_t *)vout[3];
	uint_t *tmp_tris = (uint_t *)vout[5];

	free((void *)vout[0]);
	free((void *)vout[1]);
	free((void *)vout[4]);
	free(vout);

	// get dual edges of the fully periodic graph
	uint_t *tmp_dual_edges =
			get_voro_dual_edges(tmp_num_edges, tmp_num_verts, tmp_edges, tmp_tris);

	// when use_dual is specificed, the edge and vertex sets are swapped with the
	// dual edge and dual vertex sets.  once formally relabelled, everything
	// works the same way
	if (use_dual) {
		uint_t *tmp_tmp_dual_edges = tmp_edges;
		double *tmp_lattice = tmp_vert_coords;
		uint_t tmp_num = tmp_num_verts;

		tmp_edges = tmp_dual_edges;
		tmp_dual_edges = tmp_tmp_dual_edges;

		tmp_vert_coords = lattice;
		lattice = tmp_lattice;

		tmp_num_verts = num;
		num = tmp_num;
	}

	// prune the edges of the lattice and assign boundary vertices based on the
	// desired boundary conditions
	uint_t num_bounds;
	uint_t num_verts;
	double *vert_coords;
	uint_t *bound_inds;
	uint_t *bound_verts;
	uint_t num_edges;
	uint_t *edges;
	uint_t *dual_edges;
	switch (boundary) {
		case FREE_BOUND: {
			num_bounds = 4;
			bound_inds = (uint_t *)malloc((1 + num_bounds) * sizeof(uint_t));
			bound_inds[0] = 0;
			vert_coords = tmp_vert_coords;
			num_verts = tmp_num_verts;
			num_edges = 0;
			edges = (uint_t *)malloc(2 * tmp_num_edges * sizeof(uint_t));
			dual_edges = (uint_t *)malloc(2 * tmp_num_edges * sizeof(uint_t));
			uint_t num_t, num_b, num_l, num_r;
			bool *bound_top, *bound_b, *bound_l, *bound_r;
			num_t = 0; num_b = 0; num_l = 0; num_r = 0;
			bound_top = (bool *)calloc(num_verts, sizeof(bool));
			bound_b = (bool *)calloc(num_verts, sizeof(bool));
			bound_l = (bool *)calloc(num_verts, sizeof(bool));
			bound_r = (bool *)calloc(num_verts, sizeof(bool));
			for (uint_t i = 0; i < tmp_num_edges; i++) {
				uint_t v1, v2;
				double v1x, v1y, v2x, v2y, dx, dy;
				v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1];
				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(dy) > 0.5) {
					if (dy > 0) {
						if (!bound_top[v1] && !bound_l[v1] && !bound_r[v1]) {
							bound_top[v1] = true; num_t++;
						}
						if (!bound_b[v2] && !bound_l[v2] && !bound_r[v2]) {
							bound_b[v2] = true; num_b++;
						}
					} else {
						if (!bound_top[v2] && !bound_l[v2] && !bound_r[v2]) {
							bound_top[v2] = true; num_t++;
						}
						if (!bound_b[v1] && !bound_l[v1] && !bound_r[v1]) {
							bound_b[v1] = true; num_b++;
						}
					}
				} else if (fabs(dx) > 0.5) {
					if (dx > 0) {
						if (!bound_r[v1] && !bound_top[v1] && !bound_b[v1]) {
							bound_r[v1] = true; num_r++;
						}
						if (!bound_l[v2] && !bound_top[v2] && !bound_b[v2]) {
							bound_l[v2] = true; num_l++;
						}
					} else {
						if (!bound_r[v2] && !bound_top[v2] && !bound_b[v2]) {
							bound_r[v2] = true; num_r++;
						}
						if (!bound_l[v1] && !bound_top[v1] && !bound_b[v1]) {
							bound_l[v1] = true; num_l++;
						}
					}
				} else {
					edges[2 * num_edges] = v1 < v2 ? v1 : v2;
					edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1;
					uint_t d1 = tmp_dual_edges[2 * i];
					uint_t d2 = tmp_dual_edges[2 * i + 1];
					dual_edges[2 * num_edges] = d1 < d2 ? d1 : d2;
					dual_edges[2 * num_edges + 1] = d1 < d2 ? d2 : d1;
					num_edges++;
				}
			}
			bound_verts = malloc((num_t + num_b + num_l + num_r) * sizeof(uint_t));
			bound_inds[1] = num_t; bound_inds[2] = num_t + num_b;
			bound_inds[3] = num_l + num_t + num_b;
			bound_inds[4] = num_t + num_b + num_r + num_l;
			uint_t pos_t, pos_b, pos_l, pos_r;
			pos_t = 0; pos_b = 0; pos_l = 0; pos_r = 0;
			for (uint_t i = 0; i < num_verts; i++) {
				if (bound_top[i]) {
					bound_verts[pos_t] = i; pos_t++;
				} else if (bound_b[i]) {
					bound_verts[num_t + pos_b] = i; pos_b++;
				} else if (bound_l[i]) {
					bound_verts[num_t + num_b + pos_l] = i; pos_l++;
				} else if (bound_r[i]) {
					bound_verts[num_t + num_b + num_l + pos_r] = i; pos_r++;
				}
			}
			free(bound_l); free(bound_r); free(bound_top); free(bound_b);
			free(tmp_edges);
			free(tmp_dual_edges);
			num_bounds = 2;
			g->ev_break = edges;
			g->ev = edges;
			g->nv_break = num_verts;
			g->nv = num_verts;
			g->break_dim = num_verts + num_bounds;
			break;
										 }
		case CYLINDER_BOUND: {
			num_bounds = 2;
			bound_inds = (uint_t *)malloc((1 + num_bounds) * sizeof(uint_t));
			bound_inds[0] = 0;
			vert_coords = tmp_vert_coords;
			num_verts = tmp_num_verts;
			num_edges = 0;
			edges = (uint_t *)malloc(2 * tmp_num_edges * sizeof(uint_t));
			dual_edges = (uint_t *)malloc(2 * tmp_num_edges * sizeof(uint_t));
			uint_t num_t, num_b;
			bool *bound_top, *bound_b;
			num_t = 0; num_b = 0;
			bound_top = (bool *)calloc(num_verts, sizeof(bool));
			bound_b = (bool *)calloc(num_verts, sizeof(bool));
			for (uint_t i = 0; i < tmp_num_edges; i++) {
				uint_t v1, v2;
				double v1y, v2y, dy;
				v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1];
				v1y = vert_coords[2 * v1 + 1]; v2y = vert_coords[2 * v2 + 1];
				dy = v1y - v2y;
				if (fabs(dy) > 0.5) {
					if (dy > 0) {
						if (!bound_top[v1]) {
							bound_top[v1] = true; num_t++;
						}
						if (!bound_b[v2]) {
							bound_b[v2] = true; num_b++;
						}
					} else {
						if (!bound_top[v2]) {
							bound_top[v2] = true; num_t++;
						}
						if (!bound_b[v1]) {
							bound_b[v1] = true; num_b++;
						}
					}
				} else {
					edges[2 * num_edges] = v1 < v2 ? v1 : v2;
					edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1;
					uint_t d1 = tmp_dual_edges[2 * i];
					uint_t d2 = tmp_dual_edges[2 * i + 1];
					dual_edges[2 * num_edges] = d1 < d2 ? d1 : d2;
					dual_edges[2 * num_edges + 1] = d1 < d2 ? d2 : d1;
					num_edges++;
				}
			}
			bound_verts = malloc((num_t + num_b) * sizeof(uint_t));
			bound_inds[1] = num_t; bound_inds[2] = num_t + num_b;
			uint_t pos_t, pos_b;
			pos_t = 0; pos_b = 0;
			for (uint_t i = 0; i < num_verts; i++) {
				if (bound_top[i]) {
					bound_verts[pos_t] = i; pos_t++;
				} else if (bound_b[i]) {
					bound_verts[num_t + pos_b] = i; pos_b++;
				}
			}
			free(bound_top); free(bound_b);
			free(tmp_edges);
			free(tmp_dual_edges);
			g->ev_break = edges;
			g->ev = edges;
			g->nv_break = num_verts;
			g->nv = num_verts;
			g->break_dim = num_verts + num_bounds;
			break;
												 }
		case TORUS_BOUND: {
			num_bounds = 1;
			bound_inds = (uint_t *)malloc((1 + num_bounds) * sizeof(uint_t));
			bound_inds[0] = 0;
			num_edges = tmp_num_edges;
			edges = (uint_t *)malloc(2* num_edges*sizeof(uint_t));
			for (uint_t i = 0; i < num_edges; i++) {
				edges[2*i] = tmp_edges[2*i];
				edges[2*i+1] = tmp_edges[2*i+1];
			}
			dual_edges = tmp_dual_edges;
			bool *bound_top = (bool *)calloc(tmp_num_verts, sizeof(bool));
			int *edge_change = (int *)calloc(num_edges, sizeof(int));
			uint_t num_t = 0;
			for (uint_t i = 0; i < num_edges; i++) {
				uint_t v1, v2;
				double v1y, v2y, dy;
				v1 = edges[2 * i]; v2 = edges[2 * i + 1];
				v1y = tmp_vert_coords[2 * v1 + 1];
				v2y = tmp_vert_coords[2 * v2 + 1];
				dy = v1y - v2y;
				if (fabs(dy) > 0.5) {
					if (dy > 0) {
						edge_change[i] = 1;
						if (!bound_top[v1]) {
							bound_top[v1] = true; num_t++;
						}
					} else {
						edge_change[i] = 2;
						if (!bound_top[v2]) {
							bound_top[v2] = true; num_t++;
						}
					}
				}
			}
			num_verts = tmp_num_verts + num_t;
			vert_coords = (double *)malloc(2 * num_verts * sizeof(double));
			bound_verts = malloc(num_t * sizeof(uint_t));
			bound_inds[1] = num_t;
			uint_t pos_t = 0;
			for (uint_t i = 0; i < tmp_num_verts; i++) {
				vert_coords[2*i] = tmp_vert_coords[2*i];
				vert_coords[2*i+1] = tmp_vert_coords[2*i+1];
				if (bound_top[i]) {
					bound_verts[pos_t] = i;
					vert_coords[2*(tmp_num_verts + pos_t)] = tmp_vert_coords[2*i];
					vert_coords[2*(tmp_num_verts + pos_t)+1] = tmp_vert_coords[2*i+1];
					pos_t++;
				}
			}
			for (uint_t i = 0; i < num_edges; i++) {
				if (edge_change[i]) {
					for (uint_t j = 0; j < num_t; j++) {
						if (edges[2*i+(edge_change[i]-1)] == bound_verts[j]) {
							edges[2*i+(edge_change[i]-1)] = tmp_num_verts + j;
							break;
						}
					}
				}
			}
			free(tmp_vert_coords);
			free(bound_top);
			free(edge_change);
			g->nv_break = num_verts;
			g->nv = tmp_num_verts;
			g->ev_break = edges;
			g->ev = tmp_edges;
			g->break_dim = num_verts;
			break;
											}
		case EMBEDDED_BOUND: {
			num_bounds = 4;
			bound_inds = (uint_t *)malloc(5 * sizeof(uint_t));
			bound_verts = (uint_t *)malloc(2 * L * sizeof(uint_t));
			for (uint_t i = 0; i < 5; i++) bound_inds[i] = i * L / 2;
			for (uint_t i = 0; i < 2 * L; i++) bound_verts[i] = i;
			uint_t num_away = 0;
			for (uint_t i = 0; i < tmp_num_edges; i++) {
				if (tmp_dual_edges[2*i] > num || tmp_dual_edges[2*i+1] > num) num_away++;
			}
			num_edges = (int)tmp_num_edges - (int)num_away;
			num_verts = tmp_num_verts;
			edges = tmp_edges;
			dual_edges = tmp_dual_edges;
			vert_coords = tmp_vert_coords;
			g->nv_break = num_verts;
			g->nv = num_verts;
			g->ev_break = edges;
			g->ev = edges;
			g->break_dim = num_verts + 2;
		}
	}

	g->boundary = boundary;
	g->num_bounds = num_bounds;
	g->bound_inds = bound_inds;
	g->bound_verts = bound_verts;
	g->ne = num_edges;
	g->dev = dual_edges;
	g->vx = vert_coords;

	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->L = L;

	g->ex = get_edge_coords(
			g->ne, g->vx, g->ev);

	free(tmp_tris);

	g->dvx = lattice;
	g->dnv = num;

	g->voltcurmat = gen_voltcurmat(g->ne,
																			 g->break_dim,
																			 g->ev_break, c);

	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->spanning_edges = get_spanning_edges(g->ne, g->ev_break, g->vx, 0.5, &(g->num_spanning_edges));

	return g;
}

graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c) {
	bool side_bounds;
	switch (lattice) {
		case (VORONOI_LATTICE):
			return ini_voro_graph(L, bound, dual, genfunc_uniform, c);
		case (SQUARE_LATTICE):
			if (bound == EMBEDDED_BOUND) side_bounds = true;
			else side_bounds = false;
			return ini_square_network(L, bound, side_bounds, c);
		default:
			printf("lattice type unsupported\n");
			exit(EXIT_FAILURE);
	}
}