more refactoring

17 files changed, 266 insertions, 309 deletions

diff --git a/src/break_edge.c b/src/break_edge.c
index daed1c5..54eaf34 100644
--- a/src/break_edge.c
+++ b/src/break_edge.c
@@ -7,7 +7,7 @@ bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) {
 	unsigned int v1 = instance->graph->ev_break[2 * edge];
 	unsigned int v2 = instance->graph->ev_break[2 * edge + 1];
 
-	if (instance->factor != NULL) update_factor(instance->factor, v1, v2, c);
+	if (instance->factor != NULL) factor_update(instance->factor, v1, v2, c);
 
 	if (instance->graph->boundary != TORUS_BOUND) {
 		unsigned int w1 = instance->graph->ev[2 * edge];
@@ -93,7 +93,7 @@ bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) {
 				lap_x[lap_p[tw2] + i] = 1;
 		}
 
-		set_connected(instance->adjacency, instance->dual_marks, dw1, instance->dual_marks[v1], -1, 0);
+		set_connected(instance->adjacency, instance->dual_marks, dw1, instance->dual_marks[dw2], -1, 0);
 	}
 
 
diff --git a/src/corr_test.c b/src/corr_test.c
index 4805fb0..3d56c07 100644
--- a/src/corr_test.c
+++ b/src/corr_test.c
@@ -7,8 +7,8 @@ int main() {
 
 	unsigned int width = 64;
 
-	graph_t *network = ini_square_network(width, true, false, &c);
-	net_t *instance = net_create(network, 1e-14, 0.001, 0, true, &c);
+	graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 64, false, &c);
+	net_t *instance = net_create(network, 1e-14, 0.5, 0, true, &c);
 	net_fracture(instance, &c, 1e-10);
 
 	double *corr = get_corr(instance, NULL, &c);
diff --git a/src/factor_update.c b/src/factor_update.c
new file mode 100644
index 0000000..f27971b
--- /dev/null
+++ b/src/factor_update.c
@@ -0,0 +1,38 @@
+
+#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);
+}
diff --git a/src/fracture.c b/src/fracture.c
index f36cfdb..f38568f 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -320,7 +320,7 @@ int main(int argc, char *argv[]) {
 		}
 
 		if (save_stress_field || save_voltage_field) {
-			double *tmp_voltages = get_voltage(net, &c);
+			double *tmp_voltages = net_voltages(net, &c);
 			if (save_voltage_field) {
 				for (uint_t j = 0; j < g->nv; j++) {
 					voltage_field[3 * voltage_pos] = g->vx[2 * j];
@@ -330,7 +330,7 @@ int main(int argc, char *argv[]) {
 				}
 			}
 			if (save_stress_field) {
-				double *tmp_currents = get_current_v(net, tmp_voltages, &c);
+				double *tmp_currents = net_currents(net, tmp_voltages, &c);
 				for (uint_t j = 0; j < g->ne; j++) {
 					stress_field[3 * stress_pos] = g->ex[2 * j];
 					stress_field[3 * stress_pos + 1] = g->ex[2 * j + 1];
diff --git a/src/fracture.h b/src/fracture.h
index 2401d25..5e1b3e4 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -29,6 +29,8 @@
 #define CINT_MAX INT_MAX
 #define CHOL_F(x) cholmod_##x
 
+#define GSL_RAND_GEN gsl_rng_mt19937
+
 typedef enum lattice_t {
 	VORONOI_LATTICE,
 	SQUARE_LATTICE
@@ -111,16 +113,13 @@ int dual_edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
 double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert,
 													bool index);
 
-bool update_factor(cholmod_factor *factor, unsigned int v1, unsigned int v2,
-									 cholmod_common *c);
-
-data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff);
-
-double *get_current(const net_t *instance, cholmod_common *c);
-double *get_current_v(const net_t *instance, double *voltages, cholmod_common *c);
-double *get_voltage(const net_t *instance, cholmod_common *c);
+void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, 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);
 
@@ -174,10 +173,12 @@ data_t *data_create(uint_t num_edges);
 void data_free(data_t *data);
 void data_update(data_t *data, uint_t last_broke, double strength, double conductivity);
 
-double get_conductivity(net_t *inst, double *current, cholmod_common *c);
-
 graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c);
 
 uint_t find_cycles(uint_t num_edges, const bool *fuses, const uint_t *ev, const uint_t *vei, const uint_t *ve, int **cycles);
 
 bool set_connected(const cholmod_sparse *laplacian, uint_t *marks, int vertex, int label, int stop_at, int exclude);
+
+unsigned long int rand_seed();
+
+double rand_dist_pow(const gsl_rng *r, double beta);
diff --git a/src/get_current.c b/src/get_current.c
deleted file mode 100644
index a83c399..0000000
--- a/src/get_current.c
+++ /dev/null
@@ -1,98 +0,0 @@
-
-#include "fracture.h"
-
-double *get_voltage(const net_t *instance, cholmod_common *c) {
-	cholmod_dense *b = instance->boundary_cond;
-	cholmod_factor *factor = instance->factor;
-
-	cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
-
-	if (((double *)x->x)[0] != ((double *)x->x)[0]) {
-		for (uint_t i = 0; i < instance->graph->ne; i++) { 
-			printf("%d ", instance->fuses[i]);
-		}
-		printf("\n");
-		printf("GET_VOLTAGE: value is NaN\n");
-		exit(EXIT_FAILURE);
-	}
-
-	double *field = (double *)x->x;
-	x->x = NULL;
-
-	CHOL_F(free_dense)(&x, c);
-
-
-	return field;
-}
-
-double *get_current(const net_t *instance, cholmod_common *c) {
-	unsigned int num_edges = instance->graph->ne;
-	unsigned int num_gverts = instance->graph->break_dim;
-	cholmod_sparse *voltcurmat = instance->graph->voltcurmat;
-
-	double *voltages = get_voltage(instance, c);
-	if (voltages == NULL) {
-		return NULL;
-	}
-	cholmod_dense *x = CHOL_F(allocate_dense)(
-			num_gverts, 1, num_gverts, CHOLMOD_REAL, c);
-	double *tmp_x = x->x;
-	x->x = voltages;
-
-	cholmod_dense *y =
-			CHOL_F(allocate_dense)(num_edges, 1, num_edges, CHOLMOD_REAL, c);
-
-	double alpha[2] = {1, 0};
-	double beta[2] = {0, 0};
-	CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
-
-	double *field = (double *)malloc(num_edges * sizeof(double));
-
-	for (int i = 0; i < num_edges; i++) {
-		if (instance->fuses[i])
-			field[i] = 0;
-		else
-			field[i] = ((double *)y->x)[i];
-	}
-
-	x->x = tmp_x;
-	free(voltages);
-	CHOL_F(free_dense)(&x, c);
-	CHOL_F(free_dense)(&y, c);
-
-	return field;
-}
-
-
-double *get_current_v(const net_t *instance, double *voltages, cholmod_common *c) {
-	unsigned int num_edges = instance->graph->ne;
-	unsigned int num_gverts = instance->graph->break_dim;
-	cholmod_sparse *voltcurmat = instance->graph->voltcurmat;
-
-	cholmod_dense *x = CHOL_F(allocate_dense)(
-			num_gverts, 1, num_gverts, CHOLMOD_REAL, c);
-	double *tmp_x = x->x;
-	x->x = voltages;
-
-	cholmod_dense *y =
-			CHOL_F(allocate_dense)(num_edges, 1, num_edges, CHOLMOD_REAL, c);
-
-	double alpha[2] = {1, 0};
-	double beta[2] = {0, 0};
-	CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
-
-	double *field = (double *)malloc(num_edges * sizeof(double));
-
-	for (int i = 0; i < num_edges; i++) {
-		if (instance->fuses[i])
-			field[i] = 0;
-		else
-			field[i] = ((double *)y->x)[i];
-	}
-
-	x->x = tmp_x;
-	CHOL_F(free_dense)(&x, c);
-	CHOL_F(free_dense)(&y, c);
-
-	return field;
-}
diff --git a/src/net.c b/src/net.c
new file mode 100644
index 0000000..9f6965b
--- /dev/null
+++ b/src/net.c
@@ -0,0 +1,130 @@
+
+#include "fracture.h"
+
+double *get_thres(uint_t ne, double beta) {
+	double *thres = (double *)malloc(ne * sizeof(double));
+	assert(thres != NULL);
+
+	gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN);
+	gsl_rng_set(r, 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->fuses = (bool *)calloc(g->ne, sizeof(bool));
+	assert(net->fuses != NULL);
+	net->thres = get_thres(g->ne, beta);
+	net->inf = inf;
+
+	net->voltage_bound = vb;
+	net->boundary_cond = bound_set(g, vb, notch_len, c);
+
+	if (g->boundary != TORUS_BOUND) net->adjacency = gen_adjacency(net, false, false, 0, c);
+	else net->adjacency = gen_adjacency(net, true, false, 0, c);
+
+	net->marks = (uint_t *)malloc((net->graph->break_dim) * sizeof(uint_t));
+	assert(net->marks != NULL);
+
+	net->dual_marks = (uint_t *)malloc((net->graph->dnv) * sizeof(uint_t));
+	assert(net->dual_marks != NULL);
+
+	for (uint_t i = 0; i < (net->graph->break_dim); i++) {
+		net->marks[i] = 1;
+	}
+	for (uint_t i = 0; i < (net->graph->dnv); i++) {
+		net->dual_marks[i] = i+1;
+	}
+	net->num_components = 1;
+
+	net_notch(net, notch_len, c);
+
+	{
+		cholmod_sparse *laplacian = gen_laplacian(net, c, true);
+		net->factor = CHOL_F(analyze)(laplacian, c);
+		CHOL_F(factorize)(laplacian, net->factor, c);
+		CHOL_F(free_sparse)(&laplacian, 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(double);
+	net_copy->thres = (double *)malloc(thres_size);
+	assert(net_copy->thres != NULL);
+	memcpy(net_copy->thres, net->thres, thres_size);
+
+	size_t marks_size = (net->graph->break_dim) * sizeof(uint_t);
+	net_copy->marks = (uint_t *)malloc(marks_size);
+	assert(net_copy->marks != NULL);
+	memcpy(net_copy->marks, net->marks, marks_size);
+
+	size_t dual_marks_size = (net->graph->dnv) * sizeof(uint_t);
+	net_copy->dual_marks = (uint_t *)malloc(dual_marks_size);
+	assert(net_copy->dual_marks != NULL);
+	memcpy(net_copy->dual_marks, net->dual_marks, dual_marks_size);
+
+	net_copy->adjacency = CHOL_F(copy_sparse)(net->adjacency, c);
+	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_sparse)(&(net->adjacency), c);
+	CHOL_F(free_factor)(&(net->factor), c);
+	free(net->marks);
+	free(net->dual_marks);
+	free(net);
+}
diff --git a/src/get_conductivity.c b/src/net_conductivity.c
index 23b7056..e9325bb 100644
--- a/src/get_conductivity.c
+++ b/src/net_conductivity.c
@@ -1,7 +1,7 @@
 
 #include "fracture.h"
 
-double get_conductivity(net_t *net, double *voltages, cholmod_common *c) {
+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
diff --git a/src/net_copy.c b/src/net_copy.c
deleted file mode 100644
index dcb4080..0000000
--- a/src/net_copy.c
+++ /dev/null
@@ -1,35 +0,0 @@
-
-#include "fracture.h"
-
-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(double);
-	net_copy->thres = (double *)malloc(thres_size);
-	assert(net_copy->thres != NULL);
-	memcpy(net_copy->thres, net->thres, thres_size);
-
-	size_t marks_size = (net->graph->break_dim) * sizeof(uint_t);
-	net_copy->marks = (uint_t *)malloc(marks_size);
-	assert(net_copy->marks != NULL);
-	memcpy(net_copy->marks, net->marks, marks_size);
-
-	size_t dual_marks_size = (net->graph->dnv) * sizeof(uint_t);
-	net_copy->dual_marks = (uint_t *)malloc(dual_marks_size);
-	assert(net_copy->dual_marks != NULL);
-	memcpy(net_copy->dual_marks, net->dual_marks, dual_marks_size);
-
-	net_copy->adjacency = CHOL_F(copy_sparse)(net->adjacency, c);
-	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;
-}
-
diff --git a/src/net_create.c b/src/net_create.c
deleted file mode 100644
index fedcb3a..0000000
--- a/src/net_create.c
+++ /dev/null
@@ -1,69 +0,0 @@
-
-#include "fracture.h"
-
-double *get_thres(uint_t ne, double beta) {
-	assert(beta > 0);
-
-	double *thres = (double *)malloc(ne * sizeof(double));
-	assert(thres != NULL);
-
-	gsl_set_error_handler_off();
-
-	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);
-	}
-
-	for (uint_t i = 0; i < ne; i++) {
-		while ((thres[i] = exp(log(gsl_ran_flat(r, 0, 1)) / beta)) == 0.0);
-	}
-
-	gsl_rng_free(r);
-
-	return thres;
-}
-
-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->fuses = (bool *)calloc(g->ne, sizeof(bool));
-	assert(net->fuses != NULL);
-	net->thres = get_thres(g->ne, beta);
-	net->inf = inf;
-
-	net->voltage_bound = vb;
-	net->boundary_cond = bound_set(g, vb, notch_len, c);
-
-	if (g->boundary != TORUS_BOUND) net->adjacency = gen_adjacency(net, false, false, 0, c);
-	else net->adjacency = gen_adjacency(net, true, false, 0, c);
-
-	net->marks = (uint_t *)malloc((net->graph->break_dim) * sizeof(uint_t));
-	net->dual_marks = (uint_t *)malloc((net->graph->dnv) * sizeof(uint_t));
-	assert(net->marks != NULL);
-
-	for (uint_t i = 0; i < (net->graph->break_dim); i++) {
-		net->marks[i] = 1;
-	}
-	for (uint_t i = 0; i < (net->graph->dnv); i++) {
-		net->dual_marks[i] = i+1;
-	}
-	net->num_components = 1;
-
-	net_notch(net, notch_len, c);
-
-	{
-		cholmod_sparse *laplacian = gen_laplacian(net, c, true);
-		net->factor = CHOL_F(analyze)(laplacian, c);
-		CHOL_F(factorize)(laplacian, net->factor, c);
-		CHOL_F(free_sparse)(&laplacian, c);
-	}
-
-	return net;
-}
-
diff --git a/src/net_currents.c b/src/net_currents.c
new file mode 100644
index 0000000..431818f
--- /dev/null
+++ b/src/net_currents.c
@@ -0,0 +1,35 @@
+
+#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->break_dim;
+	cholmod_sparse *voltcurmat = net->graph->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];
+		}
+	}
+
+	x->x = tmp_x;
+	CHOL_F(free_dense)(&x, c);
+	CHOL_F(free_dense)(&y, c);
+
+	return currents;
+}
diff --git a/src/net_fracture.c b/src/net_fracture.c
index 48e81f9..54f37dd 100644
--- a/src/net_fracture.c
+++ b/src/net_fracture.c
@@ -31,10 +31,10 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
 	data_t *data = data_create(net->graph->ne);
 
 	while (true) {
-		double *voltages = get_voltage(net, c);
-		double *currents = get_current_v(net, voltages, c);
+		double *voltages = net_voltages(net, c);
+		double *currents = net_currents(net, voltages, c);
 
-		double conductivity = get_conductivity(net, voltages, c);
+		double conductivity = net_conductivity(net, voltages);
 
 		if (conductivity < 1e-12 && net->graph->boundary == TORUS_BOUND) {
 			free(voltages);
diff --git a/src/net_free.c b/src/net_free.c
deleted file mode 100644
index 8b5af50..0000000
--- a/src/net_free.c
+++ /dev/null
@@ -1,14 +0,0 @@
-
-#include "fracture.h"
-
-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_sparse)(&(net->adjacency), c);
-	CHOL_F(free_factor)(&(net->factor), c);
-	free(net->marks);
-	free(net->dual_marks);
-	free(net);
-}
-
diff --git a/src/net_notch.c b/src/net_notch.c
deleted file mode 100644
index ccbb387..0000000
--- a/src/net_notch.c
+++ /dev/null
@@ -1,29 +0,0 @@
-
-#include "fracture.h"
-
-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);
-		}
-	}
-}
diff --git a/src/net_voltages.c b/src/net_voltages.c
new file mode 100644
index 0000000..dedf5b2
--- /dev/null
+++ b/src/net_voltages.c
@@ -0,0 +1,22 @@
+
+#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 *voltages = (double *)x->x;
+	x->x = NULL;
+
+	CHOL_F(free_dense)(&x, c);
+
+	return voltages;
+}
+
diff --git a/src/rand.c b/src/rand.c
new file mode 100644
index 0000000..1a6a3d4
--- /dev/null
+++ b/src/rand.c
@@ -0,0 +1,20 @@
+
+#include "fracture.h"
+
+unsigned long int rand_seed() {
+	FILE *f = fopen("/dev/urandom", "r");
+	unsigned long int seed;
+	fread(&seed, sizeof(unsigned long int), 1, f);
+	fclose(f);
+	return seed;
+}
+
+double rand_dist_pow(const gsl_rng *r, double beta) {
+	double x = 0;
+
+	// underflow means that for very small beta x is sometimes identically zero,
+	// which causes problems
+	while ((x = exp(log(gsl_rng_uniform_pos(r)) / beta)) == 0.0);
+
+	return x;
+}
diff --git a/src/update_factor.c b/src/update_factor.c
deleted file mode 100644
index debfe7b..0000000
--- a/src/update_factor.c
+++ /dev/null
@@ -1,44 +0,0 @@
-
-#include "fracture.h"
-
-bool update_factor(cholmod_factor *factor, unsigned int v1, unsigned int v2,
-									 cholmod_common *c) {
-	int n = factor->n;
-	assert(v1 < n);
-	assert(v2 < n);
-
-	cholmod_sparse *update_mat =
-			CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
-
-	unsigned int v3, v4;
-	v3 = v1 < v2 ? v1 : v2;
-	v4 = v1 > v2 ? v1 : v2;
-
-	for (int i = 0; i < n; i++) {
-		if (i <= v3)
-			((int_t *)update_mat->p)[i] = 0;
-		else
-			((int_t *)update_mat->p)[i] = 2;
-	}
-	((int_t *)update_mat->p)[n] = 2;
-	((int_t *)update_mat->i)[0] = v3;
-	((int_t *)update_mat->i)[1] = v4;
-	((double *)update_mat->x)[0] = 1;
-	((double *)update_mat->x)[1] = -1;
-
-	// assert(CHOL_F(check_sparse)(update_mat, c));
-
-	cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
-			update_mat, factor->Perm, factor->n, NULL, -1, true, true, c);
-
-	// assert(CHOL_F(check_sparse)(perm_update_mat, c));
-
-	CHOL_F(updown)(false, perm_update_mat, factor, c);
-
-	// assert(CHOL_F(check_factor)(factor, c));
-
-	CHOL_F(free_sparse)(&perm_update_mat, c);
-	CHOL_F(free_sparse)(&update_mat, c);
-
-	return true;
-}