changed code to rely on jst

16 files changed, 1036 insertions, 0 deletions

diff --git a/lib/bound_set.c b/lib/bound_set.c
new file mode 100644
index 0000000..65f1e6f
--- /dev/null
+++ b/lib/bound_set.c
@@ -0,0 +1,102 @@
+
+#include "fracture.h"
+
+double th_p(double x, double y, double th) {
+	if (x >= 0 && y >= 0) return th;
+	else if (x <  0 && y >= 0) return M_PI - th;
+	else if (x <  0 && y <  0) return th - M_PI;
+	else return -th;
+
+}
+
+double u_y(double x, double y) {
+	double r = sqrt(pow(x, 2) + pow(y, 2));
+	double th = th_p(x, y, atan(fabs(y / x)));
+
+	return sqrt(r) * sin(th / 2);
+}
+
+void bound_set_embedded(double *bound, const graph_t *g, double notch_len) {
+	uint_t L = g->L;
+
+	for (uint_t i = 0; i < L / 2; i++) {
+		double x1, y1, x2, y2, x3, y3, x4, y4;
+		x1 = (2. * i + 1.) / L - notch_len; y1 = 0.5 - 1.;
+		x2 = (2. * i + 1.) / L - notch_len; y2 = 0.5 - 0.;
+		y3 = (2. * i + 1.) / L - 0.5; x3 = 0.5 - 1.;
+		y4 = (2. * i + 1.) / L - 0.5; x4 = 0.5 - 0.;
+
+		bound[g->b[g->bi[0] + i]] = u_y(x1, y1);
+		bound[g->b[g->bi[1] + i]] = u_y(x2, y2);
+		bound[g->b[g->bi[2] + i]] = u_y(x3, y3);
+		bound[g->b[g->bi[3] + i]] = u_y(x4, y4);
+	}
+}
+
+bool is_in(uint_t len, uint_t *list, uint_t element) {
+	for (uint_t i = 0; i < len; i++) {
+		if (list[i] == element) { 
+			return true;
+		}
+	}
+	return false;
+}
+
+cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, cholmod_common *c) {
+
+	uint_t dim = g->nv;
+
+	if (vb && g->boundary != TORUS_BOUND) {
+		dim -= g->bi[g->nb];
+	} else if (!vb) {
+		dim += 2;
+	}
+
+	cholmod_dense *boundary = CHOL_F(zeros)(dim, 1, CHOLMOD_REAL, c);
+	double *bound = (double *)boundary->x;
+
+	switch (g->boundary) {
+		case TORUS_BOUND:
+			for (uint_t i = 0; i < g->bi[1]; i++) {
+				uint_t be = g->b[i];
+				uint_t v1 = g->ev[2 * be];
+				uint_t v2 = g->ev[2 * be + 1];
+				double v1y = g->vx[2 * v1 + 1];
+				double v2y = g->vx[2 * v2 + 1];
+
+				uint_t ind = v1y < v2y ? 0 : 1;
+
+				bound[g->ev[2 * be + ind]] += 1;
+				bound[g->ev[2 * be + !ind]] -= 1;
+			}
+			break;
+			/*
+		case EMBEDDED_BOUND:
+			bound_set_embedded(bound, g, notch_len);
+			break;
+			*/
+		default:
+			if (vb) {
+				for (uint_t i = 0; i < dim; i++) {
+					uint_t v = g->nbi[i];
+					for (uint_t j = 0; j < g->vei[v+1] - g->vei[v]; j++) {
+						uint_t e = g->ve[g->vei[v] + j];
+						uint_t v0 = g->ev[2 * e];
+						uint_t v1 = g->ev[2 * e + 1];
+
+						if (g->bq[v0] || g->bq[v1]) {
+							uint_t vv = v0 == v ? v1 : v0;
+							if (is_in(g->bi[1], g->b, vv)) {
+								bound[i]++;
+							}
+						}
+					}
+				}
+			} else {
+				bound[g->nv] = 1;
+				bound[g->nv + 1] = -1;
+			}
+	}
+
+	return boundary;
+}
diff --git a/lib/break_edge.c b/lib/break_edge.c
new file mode 100644
index 0000000..2f112c2
--- /dev/null
+++ b/lib/break_edge.c
@@ -0,0 +1,51 @@
+
+#include "fracture.h"
+
+
+bool break_edge(net_t *net, uint_t edge, cholmod_common *c) {
+	net->fuses[edge] = true;
+	net->num_broken++;
+
+	double *x = (double *)net->boundary_cond->x;
+	const graph_t *g = net->graph;
+
+	uint_t v1 = net->graph->ev[2 * edge];
+	uint_t v2 = net->graph->ev[2 * edge + 1];
+
+	uint_t s1 = v1 < v2 ? v1 : v2;
+	uint_t s2 = v1 < v2 ? v2 : v1;
+
+	if (net->graph->boundary == TORUS_BOUND) {
+		if (net->graph->bq[edge]) {
+			double v1y = g->vx[2 * v1 + 1];
+			double v2y = g->vx[2 * v2 + 1];
+			uint_t ind = v1y < v2y ? 0 : 1;
+			x[g->ev[2 * edge + ind]] -= 1;
+			x[g->ev[2 * edge + !ind]] += 1;
+		}
+		if (net->factor != NULL) {
+			factor_update(net->factor, s1, s2, c);
+		}
+	} else if (net->voltage_bound) {
+		if (g->bq[v1] || g->bq[v2]) {
+			uint_t vv = g->bq[v1] ? v2 : v1;
+			uint_t uu = v1 == vv ? v2 : v1;
+			if (is_in(g->bi[1], g->b, uu)) {
+				x[g->bni[vv]] -= 1;
+			}
+			if (net->factor != NULL) {
+				factor_update2(net->factor, g->bni[vv], c);
+			}
+		} else {
+			if (net->factor != NULL) {
+				factor_update(net->factor, g->bni[s1], g->bni[s2], c);
+			}
+		}
+	} else {
+		if (net->factor != NULL) {
+			factor_update(net->factor, s1, s2, c);
+		}
+	}
+
+	return true;
+}
diff --git a/lib/correlations.c b/lib/correlations.c
new file mode 100644
index 0000000..63532ec
--- /dev/null
+++ b/lib/correlations.c
@@ -0,0 +1,47 @@
+
+#include "fracture.h"
+
+double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c) {
+	uint_t nv = instance->graph->dnv;
+	uint_t ne = instance->graph->ne;
+	uint_t *ev = instance->graph->dev;
+	bool nulldists = false;
+	if (dists == NULL) {
+		dists = graph_dists(instance->graph, false);
+		nulldists = true;
+	}
+	double *corr = calloc(nv, sizeof(double));
+	uint_t *numat = calloc(nv, sizeof(uint_t));
+
+	for (uint_t i = 0; i < ne; i++) {
+		uint_t v1 = ev[2 * i];
+		uint_t v2 = ev[2 * i + 1];
+		for (uint_t j = 0; j < ne; j++) {
+			uint_t v3 = ev[2 * j];
+			uint_t v4 = ev[2 * j + 1];
+			uint_t dist1 = dists[v1][v3];
+			uint_t dist2 = dists[v1][v4];
+			uint_t dist3 = dists[v2][v3];
+			uint_t dist4 = dists[v2][v4];
+			uint_t dist = (dist1 + dist2 + dist3 + dist4) / 4;
+			corr[dist] += instance->fuses[i] && instance->fuses[j];
+			numat[dist]++;
+		}
+	}
+
+	for (uint_t i = 0; i < nv; i++) {
+		if (numat[i] > 0) {
+			corr[i] /= numat[i];
+		}
+	}
+
+	if (nulldists) {
+		for (int i = 0; i < nv; i++) {
+			free(dists[i]);
+		}
+		free(dists);
+	}
+
+	return corr;
+}
+
diff --git a/lib/data.c b/lib/data.c
new file mode 100644
index 0000000..2047c44
--- /dev/null
+++ b/lib/data.c
@@ -0,0 +1,35 @@
+
+#include "fracture.h"
+
+data_t *data_create(uint_t ne) {
+	data_t *data = malloc(1 * sizeof(data_t));
+	assert(data != NULL);
+
+	data->num_broken = 0;
+
+	data->break_list = (uint_t *)malloc(ne * sizeof(uint_t));
+	assert(data->break_list != NULL);
+
+	data->extern_field = (long double *)malloc(ne * sizeof(long double));
+	assert(data->extern_field != NULL);
+
+	data->conductivity = (double *)malloc(ne * sizeof(double));
+	assert(data->conductivity != NULL);
+
+	return data;
+}
+
+void data_free(data_t *data) {
+	free(data->break_list);
+	free(data->extern_field);
+	free(data->conductivity);
+	free(data);
+}
+
+void data_update(data_t *data, uint_t last_broke, long double strength, double conductivity) {
+	data->break_list[data->num_broken] = last_broke;
+	data->extern_field[data->num_broken] = strength;
+	data->conductivity[data->num_broken] = conductivity;
+	data->num_broken++;
+}
+
diff --git a/lib/factor_update.c b/lib/factor_update.c
new file mode 100644
index 0000000..a51705a
--- /dev/null
+++ b/lib/factor_update.c
@@ -0,0 +1,68 @@
+
+#include "fracture.h"
+
+void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common *c) {
+	uint_t n = factor->n;
+
+	cholmod_sparse *update_mat =
+			CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
+
+	uint_t s1, s2;
+	s1 = v1 < v2 ? v1 : v2;
+	s2 = v1 > v2 ? v1 : v2;
+
+	int_t *pp = (int_t *)update_mat->p;
+	int_t *ii = (int_t *)update_mat->i;
+	double *xx = (double *)update_mat->x;
+
+	for (uint_t i = 0; i <= s1; i++) {
+		pp[i] = 0;
+	}
+
+	for (uint_t i = s1 + 1; i <= n; i++) {
+		pp[i] = 2;
+	}
+
+	ii[0] = s1;
+	ii[1] = s2;
+	xx[0] = 1;
+	xx[1] = -1;
+
+	cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
+			update_mat, factor->Perm, factor->n, NULL, -1, true, true, c);
+
+	CHOL_F(updown)(false, perm_update_mat, factor, c);
+
+	CHOL_F(free_sparse)(&perm_update_mat, c);
+	CHOL_F(free_sparse)(&update_mat, c);
+}
+
+void factor_update2(cholmod_factor *factor, uint_t v, cholmod_common *c) {
+	uint_t n = factor->n;
+
+	cholmod_sparse *update_mat =
+			CHOL_F(allocate_sparse)(n, n, 1, true, true, 0, CHOLMOD_REAL, c);
+
+	int_t *pp = (int_t *)update_mat->p;
+	int_t *ii = (int_t *)update_mat->i;
+	double *xx = (double *)update_mat->x;
+
+	for (uint_t i = 0; i <= v; i++) {
+		pp[i] = 0;
+	}
+
+	for (uint_t i = v + 1; i <= n; i++) {
+		pp[i] = 1;
+	}
+
+	ii[0] = v;
+	xx[0] = 1;
+
+	cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
+			update_mat, factor->Perm, factor->n, NULL, -1, true, true, c);
+
+	CHOL_F(updown)(false, perm_update_mat, factor, c);
+
+	CHOL_F(free_sparse)(&perm_update_mat, c);
+	CHOL_F(free_sparse)(&update_mat, c);
+}
diff --git a/lib/fracture.h b/lib/fracture.h
new file mode 100644
index 0000000..b1114fb
--- /dev/null
+++ b/lib/fracture.h
@@ -0,0 +1,123 @@
+
+#pragma once
+
+#include <assert.h>
+#include <cholmod.h>
+#include <float.h>
+#include <getopt.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_sf_erf.h>
+#include <gsl/gsl_sf_exp.h>
+#include <gsl/gsl_sf_log.h>
+#include <inttypes.h>
+#include <math.h>
+#include <omp.h>
+#include <stdbool.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <unistd.h>
+
+#include <jst/graph.h>
+#include <jst/rand.h>
+
+
+// these defs allow me to switch to long int cholmod in a sitch
+#define int_t int
+#define uint_t unsigned int
+#define CINT_MAX INT_MAX
+#define CHOL_F(x) cholmod_##x
+
+#define GSL_RAND_GEN gsl_rng_mt19937
+
+typedef struct {
+	const graph_t *graph;
+	bool *fuses;
+	long double *thres;
+	double inf;
+	cholmod_dense *boundary_cond;
+	cholmod_factor *factor;
+	bool voltage_bound;
+	uint_t num_broken;
+	uint_t dim;
+	uint_t nep;
+	uint_t *evp;
+	cholmod_sparse *voltcurmat;
+} net_t;
+
+typedef struct {
+	uint_t num_broken;
+	uint_t *break_list;
+	double *conductivity;
+	long double *extern_field;
+} data_t;
+
+intptr_t *run_voronoi(uint_t num_coords, double *coords, bool periodic, double xmin, double xmax, double ymin, double ymax);
+
+cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, bool symmetric, cholmod_common *c);
+
+cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c);
+
+int edge_to_verts(uint_t width, bool periodic, uint_t edge,
+									bool index);
+
+int dual_edge_to_verts(uint_t width, bool periodic, uint_t edge,
+											 bool index);
+
+double dual_vert_to_coord(uint_t width, bool periodic, uint_t vert,
+													bool index);
+
+void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common *c);
+void factor_update2(cholmod_factor *factor, uint_t v1, cholmod_common *c);
+
+void net_notch(net_t *net, double notch_len, cholmod_common *c);
+data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff);
+double *net_voltages(const net_t *net, cholmod_common *c);
+double *net_currents(const net_t *net, const double *voltages, cholmod_common *c);
+double net_conductivity(const net_t *net, const double *voltages);
+
+void update_boundary(net_t *instance, const double *avg_field);
+
+FILE *get_file(const char *prefix, uint_t width, uint_t crack,
+							 double beta, uint_t iter, uint_t num_iter,
+							 uint_t num, bool read);
+
+double update_beta(double beta, uint_t width, const double *stress,
+									 const double *damage, double bound_total);
+
+cholmod_sparse *gen_voltcurmat(uint_t num_edges, uint_t num_verts,
+															 uint_t *edges_to_verts, cholmod_common *c);
+
+net_t *net_copy(const net_t *net, cholmod_common *c);
+
+void net_free(net_t *instance, cholmod_common *c);
+
+net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c);
+
+bool break_edge(net_t *instance, uint_t edge, cholmod_common *c);
+
+components_t *get_clusters(net_t *instance);
+
+uint_t *get_cluster_dist(net_t *instance);
+
+void randfunc_flat(gsl_rng *r, double *x, double *y);
+void randfunc_gaus(gsl_rng *r, double *x, double *y);
+
+double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c);
+
+double *bin_values(graph_t *network, uint_t width, double *values);
+
+cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, cholmod_common *c);
+
+data_t *data_create(uint_t num_edges);
+void data_free(data_t *data);
+void data_update(data_t *data, uint_t last_broke, long double strength, double conductivity);
+
+long double rand_dist_pow(const gsl_rng *r, double beta);
+
+bool is_in(uint_t len, uint_t *list, uint_t element);
diff --git a/lib/gen_laplacian.c b/lib/gen_laplacian.c
new file mode 100644
index 0000000..1cc93a4
--- /dev/null
+++ b/lib/gen_laplacian.c
@@ -0,0 +1,137 @@
+
+#include "fracture.h"
+
+cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, bool symmetric, cholmod_common *c) {
+	const graph_t *g = net->graph;
+	uint_t nv;
+	uint_t ne;
+	uint_t nre;
+	uint_t *ev;
+
+	if (use_gp) {
+		nv = net->dim;
+		ne = net->nep;
+		nre = (int_t)net->nep - (int_t)net->num_broken;
+		ev = net->evp;
+	} else if (dual) {
+		nv = g->dnv;
+		ne = g->ne;
+		nre = net->num_broken;
+		ev = g->dev;
+	} else {
+		nv = g->nv;
+		ne = g->ne;
+		nre = (int_t)g->ne - (int_t)net->num_broken;
+		ev = g->ev;
+	}
+
+	uint_t nnz = nre;
+
+	cholmod_triplet *t = CHOL_F(allocate_triplet)(nv, nv, nnz, 1, 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 a = 0;
+
+	for (uint_t i = 0; i < ne; i++) {
+		if ((net->fuses[i] && dual) || (!net->fuses[i] && !dual)) {
+			uint_t v1 = ev[2 * i];
+			uint_t v2 = ev[2 * i + 1];
+
+			uint_t s1 = v1 < v2 ? v1 : v2;
+			uint_t s2 = v1 < v2 ? v2 : v1;
+
+			ri[a] = s2;
+			ci[a] = s1;
+			ai[a] = 1;
+
+			a++;
+		}
+	}
+
+	cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c);
+	CHOL_F(free_triplet)(&t, c);
+
+	if (!symmetric) {
+		cholmod_sparse *tmp_s = CHOL_F(copy)(s, 0, 1, c);
+		CHOL_F(free_sparse)(&s, c);
+		s = tmp_s;
+	}
+
+	return s;
+}
+
+cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c) {
+	const graph_t *g = net->graph;
+	uint_t nv = net->dim;
+	uint_t ne = net->nep;
+	uint_t *ev = net->evp;
+
+	uint_t nnz = nv;
+
+	cholmod_triplet *temp_m = CHOL_F(allocate_triplet)(nv, nv, nnz, 1, 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 < nv; i++) {
+		rowind[i] = i;
+		colind[i] = i;
+		acoo[i] = 0;
+	}
+
+	cholmod_sparse *adjacency = gen_adjacency(net, false, true, true, c);
+
+	for (uint_t i = 0; i < ne; i++) {
+		if (!net->fuses[i]){
+			uint_t v1 = ev[2 * i];
+			uint_t v2 = ev[2 * i + 1];
+
+			acoo[v1]++;
+			acoo[v2]++;
+		}
+	}
+
+	if (net->voltage_bound && g->boundary != TORUS_BOUND) {
+		for (uint_t i = 0; i < net->dim; i++) {
+			uint_t v = g->nbi[i];
+			for (uint_t j = 0; j < g->vei[v+1] - g->vei[v]; j++) {
+				uint_t e = g->ve[g->vei[v] + j];
+				uint_t v0 = g->ev[2 * e];
+				uint_t v1 = g->ev[2 * e + 1];
+
+				if (g->bq[v0] || g->bq[v1]) {
+					acoo[i]++;
+				}
+			}
+		}
+	} else {
+		acoo[0]++;
+	}
+
+	for (uint_t i = 0; i < nv; i++) {
+		if (acoo[i] == 0) 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 *laplacian = CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c);
+
+	CHOL_F(free_sparse)(&t_out, c);
+	CHOL_F(free_sparse)(&adjacency, c);
+	CHOL_F(free_triplet)(&temp_m, c);
+
+	return laplacian;
+}
diff --git a/lib/gen_voltcurmat.c b/lib/gen_voltcurmat.c
new file mode 100644
index 0000000..e870140
--- /dev/null
+++ b/lib/gen_voltcurmat.c
@@ -0,0 +1,36 @@
+
+#include "fracture.h"
+
+cholmod_sparse *gen_voltcurmat(unsigned int num_edges, unsigned int num_verts,
+															 unsigned int *edges_to_verts,
+															 cholmod_common *c) {
+
+	cholmod_triplet *t_m = CHOL_F(allocate_triplet)(
+			num_edges, num_verts, num_edges * 2, 0, CHOLMOD_REAL, c);
+	assert(t_m != NULL);
+
+	int_t *rowind = (int_t *)t_m->i;
+	int_t *colind = (int_t *)t_m->j;
+	double *acoo = (double *)t_m->x;
+
+	for (unsigned int i = 0; i < num_edges; i++) {
+		unsigned int v1 = edges_to_verts[2 * i];
+		unsigned int v2 = edges_to_verts[2 * i + 1];
+		rowind[2 * i] = i;
+		rowind[2 * i + 1] = i;
+		colind[2 * i] = v1;
+		colind[2 * i + 1] = v2;
+		acoo[2 * i] = 1;
+		acoo[2 * i + 1] = -1;
+	}
+
+	t_m->nnz = num_edges * 2;
+
+	assert(CHOL_F(check_triplet)(t_m, c));
+
+	cholmod_sparse *m = CHOL_F(triplet_to_sparse)(t_m, num_edges * 2, c);
+
+	CHOL_F(free_triplet)(&t_m, c);
+
+	return m;
+}
diff --git a/lib/geometry.c b/lib/geometry.c
new file mode 100644
index 0000000..ec788f1
--- /dev/null
+++ b/lib/geometry.c
@@ -0,0 +1,55 @@
+
+#include "fracture.h"
+
+int edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
+									bool index) {
+	assert(edge < pow(width, 2));
+
+	int x = edge / width + 1;
+	int y = edge % width + 1;
+
+	if (periodic) {
+		return (((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) - 1) %
+								width +
+						(x - index) * width) /
+					 2;
+	} else {
+		return ((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) +
+						(x - index) * (width + 1) - 1) /
+					 2;
+	}
+}
+
+int dual_edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
+											 bool index) {
+	assert(edge < pow(width, 2));
+
+	int x = edge / width + 1;
+	int y = edge % width + 1;
+
+	if (periodic) {
+		return (((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) - 1) %
+								width +
+						(x - index) * width) /
+					 2;
+	} else {
+		return ((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) +
+						(x - index) * (width + 1) - 1) /
+					 2;
+	}
+}
+
+double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert,
+													bool index) {
+	if (periodic) {
+		if (index)
+			return (2 * vert) % width + (2 * vert / width) % 2;
+		else
+			return 2 * vert / width;
+	} else {
+		if (index)
+			return (2 * vert) % (width + 1);
+		else
+			return (2 * vert) / (width + 1);
+	}
+}
diff --git a/lib/get_dual_clusters.c b/lib/get_dual_clusters.c
new file mode 100644
index 0000000..eaf7562
--- /dev/null
+++ b/lib/get_dual_clusters.c
@@ -0,0 +1,30 @@
+
+#include "fracture.h"
+
+components_t *get_clusters(net_t *instance) {
+	components_t *c = graph_components_get(instance->graph, instance->fuses, true);
+	return c;
+}
+
+unsigned int *get_cluster_dist(net_t *instance) {
+	components_t *c = get_clusters(instance);
+	unsigned int *cluster_dist = (unsigned int *)calloc(
+			instance->graph->dnv, sizeof(unsigned int));
+
+	for (uint32_t i = 1; i <= c->n; i++) {
+		unsigned int num_in_cluster = 0;
+		for (unsigned int j = 0; j < instance->graph->dnv; j++) {
+			if (c->labels[j] == i)
+				num_in_cluster++;
+		}
+
+		if (num_in_cluster == 0)
+			break;
+
+		cluster_dist[num_in_cluster - 1]++;
+	}
+
+	graph_components_free(c);
+
+	return cluster_dist;
+}
diff --git a/lib/net.c b/lib/net.c
new file mode 100644
index 0000000..eb44721
--- /dev/null
+++ b/lib/net.c
@@ -0,0 +1,144 @@
+
+#include "fracture.h"
+
+long double *get_thres(uint_t ne, double beta) {
+	long double *thres = (long double *)malloc(ne * sizeof(long double));
+	assert(thres != NULL);
+
+	gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN);
+	gsl_rng_set(r, jst_rand_seed());
+
+	for (uint_t i = 0; i < ne; i++) {
+		thres[i] = rand_dist_pow(r, beta);
+	}
+
+	gsl_rng_free(r);
+
+	return thres;
+}
+
+void net_notch(net_t *net, double notch_len, cholmod_common *c) {
+	for (uint_t i = 0; i < net->graph->ne; i++) {
+		uint_t v1, v2;
+		double v1x, v1y, v2x, v2y, dy;
+		bool crosses_center, not_wrapping, correct_length;
+
+		v1 = net->graph->ev[2 * i];
+		v2 = net->graph->ev[2 * i + 1];
+
+		v1x = net->graph->vx[2 * v1];
+		v1y = net->graph->vx[2 * v1 + 1];
+		v2x = net->graph->vx[2 * v2];
+		v2y = net->graph->vx[2 * v2 + 1];
+
+		dy = v1y - v2y;
+
+		crosses_center = (v1y >= 0.5 && v2y <= 0.5) || (v1y <= 0.5 && v2y >= 0.5);
+		not_wrapping = fabs(dy) < 0.5;
+		//correct_length = v1x + dx / dy * (v1y - 0.5) <= notch_len;
+		correct_length = v1x < notch_len && v2x < notch_len;
+
+		if (crosses_center && not_wrapping && correct_length) {
+			break_edge(net, i, c);
+		}
+	}
+}
+
+net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c) {
+	net_t *net = (net_t *)calloc(1, sizeof(net_t));
+	assert(net != NULL);
+
+	net->graph = g;
+	net->num_broken = 0;
+
+	net->fuses = (bool *)calloc(g->ne, sizeof(bool));
+	assert(net->fuses != NULL);
+
+	net->thres = get_thres(g->ne, beta);
+	net->inf = inf;
+
+	net->dim = g->nv;
+
+	if (g->boundary == TORUS_BOUND) {
+		net->nep = g->ne;
+		net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
+		memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
+	} else {
+		if (vb) {
+			net->dim -= g->bi[g->nb];
+			net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
+			net->nep = 0;
+			for (uint_t i = 0; i < g->ne; i++) {
+				if (!(g->bq[g->ev[2*i]] || g->bq[g->ev[2*i+1]])) {
+					net->evp[2*net->nep] = g->bni[g->ev[2*i]];
+					net->evp[2*net->nep + 1] = g->bni[g->ev[2*i + 1]];
+					net->nep++;
+				}
+			}
+		} else {
+			net->dim += 2;
+			net->evp = (uint_t *)malloc(2 * (g->ne + g->bi[2]) * sizeof(uint_t));
+			memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
+			net->nep = g->ne + g->bi[2];
+
+			for (uint_t i = 0; i < 2; i++) {
+				for (uint_t j = 0; j < g->bi[i+1] - g->bi[i]; j++){ 
+					net->evp[2 * (g->ne + g->bi[i] + j)] = g->b[g->bi[i] + j];
+					net->evp[2 * (g->ne + g->bi[i] + j) + 1] = g->nv + i;
+				}
+			}
+		}
+	}
+
+	net->voltage_bound = vb;
+	net->boundary_cond = bound_set(g, vb, notch_len, c);
+
+	net_notch(net, notch_len, c);
+
+	{
+		cholmod_sparse *laplacian = gen_laplacian(net, c);
+		net->factor = CHOL_F(analyze)(laplacian, c);
+		CHOL_F(factorize)(laplacian, net->factor, c);
+		CHOL_F(free_sparse)(&laplacian, c);
+	}
+
+	net->voltcurmat = gen_voltcurmat(g->ne, g->nv, g->ev, c);
+
+	return net;
+}
+
+net_t *net_copy(const net_t *net, cholmod_common *c) {
+	net_t *net_copy = (net_t *)calloc(1, sizeof(net_t));
+	assert(net_copy != NULL);
+	memcpy(net_copy, net, sizeof(net_t));
+
+	size_t fuses_size = (net->graph)->ne * sizeof(bool);
+	net_copy->fuses = (bool *)malloc(fuses_size);
+	assert(net_copy->fuses != NULL);
+	memcpy(net_copy->fuses, net->fuses, fuses_size);
+
+	size_t thres_size = (net->graph)->ne * sizeof(long double);
+	net_copy->thres = (long double *)malloc(thres_size);
+	assert(net_copy->thres != NULL);
+	memcpy(net_copy->thres, net->thres, thres_size);
+
+	size_t evp_size = 2 * net->nep * sizeof(uint_t);
+	net_copy->evp = (uint_t *)malloc(thres_size);
+	assert(net_copy->evp != NULL);
+	memcpy(net_copy->evp, net->evp, evp_size);
+
+	net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c);
+	net_copy->factor = CHOL_F(copy_factor)(net->factor, c);
+
+	return net_copy;
+}
+
+void net_free(net_t *net, cholmod_common *c) {
+	free(net->fuses);
+	free(net->thres);
+	CHOL_F(free_dense)(&(net->boundary_cond), c);
+	CHOL_F(free_factor)(&(net->factor), c);
+	CHOL_F(free_sparse)(&(net->voltcurmat), c);
+	free(net->evp);
+	free(net);
+}
diff --git a/lib/net_conductivity.c b/lib/net_conductivity.c
new file mode 100644
index 0000000..e9325bb
--- /dev/null
+++ b/lib/net_conductivity.c
@@ -0,0 +1,35 @@
+
+#include "fracture.h"
+
+double net_conductivity(const net_t *net, const double *voltages) {
+	if (net->voltage_bound) {
+		// the voltage drop across the network is fixed to one with voltage
+		// boundary conditions, so the conductivity is the total current flowing
+		double tot_cur = 0;
+		for (uint_t i = 0; i < net->graph->num_spanning_edges; i++) {
+			uint_t e = net->graph->spanning_edges[i];
+
+			if (!net->fuses[e]) {
+				uint_t v1, v2, s1, s2;
+				double v1y, v2y;
+
+				v1 = net->graph->ev[2 * e];
+				v2 = net->graph->ev[2 * e + 1];
+
+				v1y = net->graph->vx[2 * v1 + 1];
+				v2y = net->graph->vx[2 * v2 + 1];
+
+				s1 = v1y < v2y ? v1 : v2;
+				s2 = v1y < v2y ? v2 : v1;
+
+				tot_cur += voltages[s1] - voltages[s2];
+			}
+		}
+
+		return fabs(tot_cur);
+	} else {
+		// the current across the network is fixed to one with current boundary
+		// conditions, so the conductivity is the inverse of the total voltage drop
+		return 1 / fabs(voltages[net->graph->nv] - voltages[net->graph->nv + 1]);
+	}
+}
diff --git a/lib/net_currents.c b/lib/net_currents.c
new file mode 100644
index 0000000..a078336
--- /dev/null
+++ b/lib/net_currents.c
@@ -0,0 +1,51 @@
+
+#include "fracture.h"
+
+double *net_currents(const net_t *net, const double *voltages, cholmod_common *c) {
+	uint_t ne = net->graph->ne;
+	uint_t dim = net->graph->nv;
+	cholmod_sparse *voltcurmat = net->voltcurmat;
+
+	cholmod_dense *x = CHOL_F(allocate_dense)(dim, 1, dim, CHOLMOD_REAL, c);
+
+	double *tmp_x = x->x;
+	x->x = (void *)voltages;
+
+	cholmod_dense *y = CHOL_F(allocate_dense)(ne, 1, ne, CHOLMOD_REAL, c);
+
+	double alpha[2] = {1, 0};
+	double beta[2] = {0, 0};
+	CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
+
+	double *currents = (double *)malloc(ne * sizeof(double));
+
+	for (int i = 0; i < ne; i++) {
+		if (net->fuses[i]) {
+			currents[i] = 0;
+		} else {
+			currents[i] = ((double *)y->x)[i];
+		}
+	}
+
+	if (net->graph->boundary == TORUS_BOUND) {
+		for (uint_t i = 0; i < net->graph->bi[1]; i++) {
+			uint_t e = net->graph->b[i];
+			uint_t v1 = net->graph->ev[2 * e];
+			uint_t v2 = net->graph->ev[2 * e + 1];
+			double v1y = net->graph->vx[2 * v1 + 1];
+			double v2y = net->graph->vx[2 * v2 + 1];
+
+			if (v1y > v2y) {
+				currents[e] += 1;
+			} else {
+				currents[e] -= 1;
+			}
+		}
+	}
+
+	x->x = tmp_x;
+	CHOL_F(free_dense)(&x, c);
+	CHOL_F(free_dense)(&y, c);
+
+	return currents;
+}
diff --git a/lib/net_fracture.c b/lib/net_fracture.c
new file mode 100644
index 0000000..e7f18fc
--- /dev/null
+++ b/lib/net_fracture.c
@@ -0,0 +1,67 @@
+
+#include "fracture.h"
+
+uint_t get_next_broken(net_t *net, double *currents, double cutoff) {
+	uint_t max_pos = UINT_MAX;
+	long double max_val = 0;
+
+	for (uint_t i = 0; i < net->graph->ne; i++) {
+		long double current = fabs(currents[i]);
+		bool broken = net->fuses[i];
+
+		if (!broken && current > cutoff) {
+			long 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);
+
+	uint_t n = 0;
+	while (true) {
+		n++;
+		double *voltages = net_voltages(net, c);
+		double *currents = net_currents(net, voltages, c);
+
+		double conductivity = net_conductivity(net, voltages);
+
+		if (conductivity < cutoff) {
+			free(voltages);
+			free(currents);
+			break;
+		}
+
+		uint_t last_broke = get_next_broken(net, currents, cutoff);
+
+		long double sim_current;
+
+		if (net->voltage_bound) {
+			sim_current = conductivity;
+		} else {
+			sim_current = 1;
+		}
+
+		data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] / ((long double)currents[last_broke])), conductivity);
+
+		free(voltages);
+		free(currents);
+
+		break_edge(net, last_broke, c);
+	}
+
+	return data;
+}
+
diff --git a/lib/net_voltages.c b/lib/net_voltages.c
new file mode 100644
index 0000000..7b07201
--- /dev/null
+++ b/lib/net_voltages.c
@@ -0,0 +1,40 @@
+
+#include "fracture.h"
+
+double *net_voltages(const net_t *net, cholmod_common *c) {
+	cholmod_dense *b = net->boundary_cond;
+	cholmod_factor *factor = net->factor;
+
+	cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
+
+	if (((double *)x->x)[0] != ((double *)x->x)[0]) {
+		printf("GET_VOLTAGE: value is NaN\n");
+		exit(EXIT_FAILURE);
+	}
+
+	double *t_voltages = (double *)x->x;
+	x->x = NULL;
+	CHOL_F(free_dense)(&x, c);
+
+	const graph_t *g = net->graph;
+
+	if (g->boundary == TORUS_BOUND) {
+		return t_voltages;
+	} else if (net->voltage_bound) {
+		double *voltages = (double *)malloc(g->nv * sizeof(double));
+		for (uint_t i = 0; i < g->nv - g->bi[g->nb]; i++) {
+			voltages[g->nbi[i]] = t_voltages[i];
+		}
+		for (uint_t i = 0; i < 2; i++) {
+			for (uint_t j = 0; j < g->bi[i + 1] - g->bi[i]; j++) {
+				voltages[g->b[g->bi[i] + j]] = 1 - i;
+			}
+		}
+
+		free(t_voltages);
+		return voltages;
+	} else {
+		return t_voltages;
+	}
+}
+
diff --git a/lib/rand.c b/lib/rand.c
new file mode 100644
index 0000000..16aeba4
--- /dev/null
+++ b/lib/rand.c
@@ -0,0 +1,15 @@
+
+#include "fracture.h"
+
+long double rand_dist_pow(const gsl_rng *r, double beta) {
+	long double x = 0;
+
+	// underflow means that for very small beta x is sometimes identically zero,
+	// which causes problems
+	while (x == 0.0) {
+		long double y = logl(gsl_rng_uniform_pos(r)) / beta;
+		x = expl(y);
+	}
+
+	return x;
+}