many changes, simplification of some functions, removal of unneeded ones

8 files changed, 149 insertions, 508 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index c325cb4..66d9e2c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,8 +1,8 @@
 
-cmake_minimum_required(VERSION 3.7)
+cmake_minimum_required(VERSION 3.0)
 project(wolff)
 
-include_directories(src ~/.local/include)
+include_directories(lib ~/.local/include)
 link_directories(~/.local/lib)
 
 file(GLOB SOURCES lib/*.c)
diff --git a/lib/queue.c b/lib/queue.c
index 9a01741..aebbb2a 100644
--- a/lib/queue.c
+++ b/lib/queue.c
@@ -1,23 +1,23 @@
 
-#include "wolff.h"
+#include "queue.h"
 
 void stack_push(ll_t **q, uint32_t x) {
-	ll_t *nq = malloc(sizeof(ll_t));
-	nq->x = x;
-	nq->next = *q;
+  ll_t *nq = malloc(sizeof(ll_t));
+  nq->x = x;
+  nq->next = *q;
 
-	*q = nq;
+  *q = nq;
 }
 
 uint32_t stack_pop(ll_t **q) {
-	ll_t *old_q = *q;
+  ll_t *old_q = *q;
 
-	*q = old_q->next;
-	uint32_t x = old_q->x;
+  *q = old_q->next;
+  uint32_t x = old_q->x;
 
-	free(old_q);
+  free(old_q);
 
-	return x;
+  return x;
 }
 
 bool stack_contains(const ll_t *q, uint32_t x) {
@@ -29,4 +29,3 @@ bool stack_contains(const ll_t *q, uint32_t x) {
     return stack_contains(q->next, x);
   }
 }
-
diff --git a/lib/queue.h b/lib/queue.h
new file mode 100644
index 0000000..c7acec1
--- /dev/null
+++ b/lib/queue.h
@@ -0,0 +1,18 @@
+
+#pragma once
+
+#include <inttypes.h>
+#include <stdbool.h>
+#include <stdlib.h>
+#include <string.h>
+
+typedef struct ll_tag {
+  uint32_t x;
+  struct ll_tag *next;
+} ll_t;
+
+void stack_push(ll_t **q, uint32_t x);
+
+uint32_t stack_pop(ll_t **q);
+
+bool stack_contains(const ll_t *q, uint32_t x);
diff --git a/lib/wolff.h b/lib/wolff.h
index cec9ee3..d3bc412 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -1,19 +1,23 @@
 
-#include <string.h>
-#include <math.h>
-#include <getopt.h>
+#pragma once
+
+#include <assert.h>
+#include <fftw3.h>
 #include <float.h>
-#include <sys/types.h>
-#include <inttypes.h>
+#include <getopt.h>
 #include <gsl/gsl_randist.h>
 #include <gsl/gsl_rng.h>
+#include <inttypes.h>
+#include <math.h>
 #include <stdbool.h>
-#include <assert.h>
-#include <fftw3.h>
+#include <string.h>
+#include <sys/types.h>
 
 #include <jst/graph.h>
 #include <jst/rand.h>
 
+#include "queue.h"
+
 typedef struct {
   graph_t *g;
   bool *spins;
@@ -21,32 +25,18 @@ typedef struct {
   double H;
 } ising_state_t;
 
-typedef struct ll_tag {
-	uint32_t x;
-	struct ll_tag *next;
-} ll_t;
-
 typedef struct {
   uint32_t nv;
-  double dH;
-  int32_t dM;
+  int32_t dJb;
+  int32_t dHb;
 } cluster_t;
 
-double get_hamiltonian(graph_t *g, double *coupling, bool *x);
-
-void stack_push(ll_t **q, uint32_t x);
-
-uint32_t stack_pop(ll_t **q);
-
-bool stack_contains(const ll_t *q, uint32_t x);
-
-cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, gsl_rng *r);
+cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
+                        gsl_rng *r);
 
 graph_t *graph_add_ext(const graph_t *g);
 
-double hh(double th);
-
 double *get_bond_probs(double T, double H, ising_state_t *s);
 
-uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps);
-
+uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
+                    double *ps);
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index 016e19b..f45e1b5 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -1,4 +1,5 @@
 
+#include "queue.h"
 #include "wolff.h"
 
 graph_t *graph_add_ext(const graph_t *g) {
@@ -9,7 +10,7 @@ graph_t *graph_add_ext(const graph_t *g) {
 
   h->ev = (uint32_t *)malloc(2 * h->ne * sizeof(uint32_t));
   h->vei = (uint32_t *)malloc((h->nv + 1) * sizeof(uint32_t));
-  h->ve = (uint32_t *) malloc(2 * h->ne * sizeof(uint32_t));
+  h->ve = (uint32_t *)malloc(2 * h->ne * sizeof(uint32_t));
   h->vx = (double *)malloc(2 * h->nv * sizeof(double));
   h->bq = (bool *)malloc(h->nv * sizeof(bool));
 
@@ -45,42 +46,21 @@ graph_t *graph_add_ext(const graph_t *g) {
   return h;
 }
 
-double get_hamiltonian(graph_t *g, double *coupling, bool *x) {
-  double hamiltonian = 0;
-
-  for (uint32_t i = 0; i < g->ne; i++) {
-    uint32_t v1, v2;
-
-    v1 = g->ev[2 * i];
-    v2 = g->ev[2 * i + 1];
-
-    if (x[v1] == x[v2]) {
-      hamiltonian -= coupling[i];
-    } else {
-      hamiltonian += coupling[i];
-    }
-  }
-
-  return hamiltonian;
-}
-
-cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, gsl_rng *r) {
+cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
+                        gsl_rng *r) {
   uint32_t v0;
   int32_t n_h_bonds, n_bonds;
   bool x0;
   cluster_t *c;
-  
+
   v0 = gsl_rng_uniform_int(r, g->nv); // pick a random vertex
-  x0 = x[v0]; // record its orientation
+  x0 = x[v0];                         // record its orientation
 
-  ll_t *stack = NULL; // create a new stack
+  ll_t *stack = NULL;     // create a new stack
   stack_push(&stack, v0); // push the initial vertex to the stack
 
   // initiate the data structure for returning flip information
   c = (cluster_t *)calloc(1, sizeof(cluster_t));
-  c->nv = 0;
-
-  n_h_bonds = 0; n_bonds = 0;
 
   while (stack != NULL) {
     uint32_t v;
@@ -90,7 +70,7 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, g
     nn = g->vei[v + 1] - g->vei[v];
 
     if (x[v] == x0) { // if the vertex hasn't already been flipped
-      x[v] = !x[v]; // flip the vertex
+      x[v] = !x[v];   // flip the vertex
 
       for (uint32_t i = 0; i < nn; i++) {
         uint32_t e, v1, v2, vn;
@@ -102,9 +82,11 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, g
 
         vn = v == v1 ? v2 : v1; // distinguish neighboring site from site itself
 
-        bond_counter = (v1 == g->nv - 1 || v2 == g->nv - 1) ? &n_h_bonds : &n_bonds;
+        bond_counter =
+            (v1 == g->nv - 1 || v2 == g->nv - 1) ? &(c->dHb) : &(c->dJb);
 
-        if (x[vn] == x0) { // if the neighboring site matches the flipping cluster...
+        if (x[vn] ==
+            x0) { // if the neighboring site matches the flipping cluster...
           (*bond_counter)++;
 
           if (gsl_rng_uniform(r) < ps[e]) { // and with probability ps[e]...
@@ -115,22 +97,15 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, g
         }
       }
 
-      if (v != g->nv - 1) {
+      if (v != g->nv - 1) { // count the number of non-external sites that flip
         c->nv++;
       }
     }
   }
 
-  c->dH = n_bonds + H * n_h_bonds;
-  c->dM = n_h_bonds;
-
   return c;
 }
 
-double hh(double th) {
-  return (th - pow(th, 3) / 1.16951) * (1 - 0.222389 * pow(th, 2) - 0.043547 * pow(th, 4) - 0.014809 * pow(th, 6) - 0.007168 * pow(th, 8));
-}
-
 double *get_bond_probs(double T, double H, ising_state_t *s) {
   double p = 1 - exp(-2 / T);
   double q = 1 - exp(-2 * fabs(H) / T);
@@ -151,7 +126,8 @@ double *get_bond_probs(double T, double H, ising_state_t *s) {
   return ps;
 }
 
-uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps) {
+uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
+                    double *ps) {
   if (r == NULL) {
     r = gsl_rng_alloc(gsl_rng_mt19937);
     gsl_rng_set(r, jst_rand_seed());
@@ -161,10 +137,10 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps
     ps = get_bond_probs(T, H, s);
   }
 
-  cluster_t *c = flip_cluster(s->g, ps, H, s->spins, r);
+  cluster_t *c = flip_cluster(s->g, ps, s->spins, r);
 
-  s->M += -2 * c->dM;
-  s->H += 2 * c->dH;
+  s->M += -2 * c->dHb;
+  s->H += 2 * (c->dJb + H * c->dHb);
 
   uint32_t n_flips = c->nv;
 
@@ -172,4 +148,3 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps
 
   return n_flips;
 }
-
diff --git a/src/pt_wolff.c b/src/pt_wolff.c
deleted file mode 100644
index 128b601..0000000
--- a/src/pt_wolff.c
+++ /dev/null
@@ -1,317 +0,0 @@
-
-#include "wolff.h"
-
-#define TC 2. / log(1. + sqrt(2.))
-
-int main(int argc, char *argv[]) {
-	int opt;
-	bool use_scho, use_rand_ini;
-	lattice_t lat;
-	uint16_t L;
-	uint32_t n_temps;
-	uint64_t N;
-	double Tin, Hin, eps, dT, dH;
-
-	L = 128;
-	N = 1000;
-	lat = SQUARE_LATTICE;
-	Tin = 2.3;
-	Hin = 0;
-	eps = 1e30;
-	use_scho = false;
-	use_rand_ini = false;
-	dT = 0;
-	dH = 1;
-	n_temps = 1;
-
-	while ((opt = getopt(argc, argv, "N:L:q:T:H:e:St:h:n:r")) != -1) {
-		switch (opt) {
-			case 'S':
-				use_scho = true;
-				break;
-			case 'r':
-				use_rand_ini = true;
-				break;
-			case 'N':
-				N = (uint64_t)atof(optarg);
-				break;
-			case 'L':
-				L = atoi(optarg);
-				break;
-			case 'T':
-				Tin = atof(optarg);
-				break;
-			case 'H':
-				Hin = atof(optarg);
-				break;
-			case 't':
-				dT = atof(optarg);
-				break;
-			case 'h':
-				dH = atof(optarg);
-				break;
-			case 'n':
-				n_temps = atoi(optarg);
-				break;
-			case 'e':
-				eps= atof(optarg);
-				break;
-			default:
-				exit(EXIT_FAILURE);
-		}
-	}
-
-	gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
-	gsl_rng_set(r, jst_rand_seed());
-
-	graph_t *g = graph_create(lat, TORUS_BOUND, L, false);
-
-	double *Ts = (double *)malloc(n_temps * sizeof(double));
-	double *Hs = (double *)malloc(n_temps * sizeof(double));
-	double *Tins = (double *)malloc(n_temps * sizeof(double));
-	double *Hins = (double *)malloc(n_temps * sizeof(double));
-
-	for (uint32_t i = 0; i < n_temps; i++) {
-		double T, H;
-
-		// if Schofield coordinates are used, T and H are treated as R and theta.
-		if (use_scho) {
-			double R = Tin; double th = Hin;
-			double t = R * (1 - pow(th, 2));
-			T = TC / (1 - t);
-			double h0 = 0.940647;
-			double h = h0 * pow(R, 15. / 8.) * hh(th);
-			H = h * T;
-		} else {
-			T = Tin;
-			H = Hin;
-		}
-
-		Tins[i] = Tin;
-		Hins[i] = Hin;
-
-		Ts[i] = T;
-		Hs[i] = H;
-
-		Tin += dT;
-		if (use_scho) {
-			Hin = 1.08144 - (1.08144 - Hin) / dH;
-		} else {
-			Hin = Hin / dH;
-		}
-	}
-
-	bool **xs = (bool **)malloc(n_temps * sizeof(bool *));
-	for (uint32_t i = 0; i < n_temps; i++) {
-		xs[i] = (bool *)calloc(g->nv, sizeof(bool));
-	}
-
-	if (use_rand_ini) {
-		for (uint32_t j = 0; j < n_temps; j++) {
-			for (uint32_t i = 0; i < g->nv; i++) {
-				xs[j][i] = gsl_rng_uniform_int(r, 2);
-			}
-		}
-	}
-
-	int32_t *Ms = (int32_t *)calloc(n_temps, sizeof(int32_t));
-	for (uint32_t j = 0; j < n_temps; j++) {
-		if (use_rand_ini) {
-			for (uint32_t i = 0; i < g->nv; i++) {
-				if (xs[j][i]) {
-					Ms[j]++;
-				} else {
-					Ms[j]--;
-				}
-			}
-		} else {
-			Ms[j] = -g->nv;
-		}
-	}
-
-	double *Es = (double *)malloc(n_temps * sizeof(double));
-	for (uint32_t i = 0; i < n_temps; i++) {
-		Es[i] = get_hamiltonian(g, Hs[i], xs[i]);
-	}
-
-	double *ps = (double *)malloc(n_temps * sizeof(double));
-	for (uint32_t i = 0; i < n_temps; i++) {
-		ps[i] = 1 - exp(-2 / Ts[i]);
-	}
-
-	double *E1s = (double *)calloc(n_temps, sizeof(double));
-	double *E2s = (double *)calloc(n_temps, sizeof(double));
-	double *E4s = (double *)calloc(n_temps, sizeof(double));
-
-	double *M1s = (double *)calloc(n_temps, sizeof(double));
-	double *M2s = (double *)calloc(n_temps, sizeof(double));
-	double *M4s = (double *)calloc(n_temps, sizeof(double));
-
-	double *tE1s = (double *)calloc(n_temps, sizeof(double));
-	double *tE2s = (double *)calloc(n_temps, sizeof(double));
-	double *tE4s = (double *)calloc(n_temps, sizeof(double));
-
-	double *tM1s = (double *)calloc(n_temps, sizeof(double));
-	double *tM2s = (double *)calloc(n_temps, sizeof(double));
-	double *tM4s = (double *)calloc(n_temps, sizeof(double));
-
-	uint64_t n_runs = 0;
-
-	double diff = 1e30;
-
-	printf("\n");
-	while (diff > eps) {
-		printf("\033[F\033[JWOLFF: sweep %llu, %g\n", n_runs, diff);
-
-		for (uint32_t j = 0; j < n_temps; j++) {
-			tE1s[j] = 0;
-			tE2s[j] = 0;
-			tE4s[j] = 0;
-			tM1s[j] = 0;
-			tM2s[j] = 0;
-			tM4s[j] = 0;
-		}
-
-		for (uint64_t i = 0; i < N; i++) {
-#pragma omp parallel for
-			for (uint32_t j = 0; j < n_temps; j++) {
-				cluster_t *c = get_cluster(g, xs[j], ps[j], r);
-				double dHex = 0;
-				double s;
-
-				if (xs[j][c->vi->x]) {
-					s = 1;
-				} else {
-					s = -1;
-				}
-
-				dHex = s * Hs[j] * c->nv;
-
-				if (gsl_rng_uniform(r) < exp(-2 * dHex / Ts[j])) {
-					while (c->vi != NULL) {
-						uint32_t v = queue_del(&(c->vi));
-						xs[j][v] = !xs[j][v];
-					}
-					Es[j] += 2 * c->nb + dHex;
-					Ms[j] -= 2 * s * c->nv;
-				} else {
-					while (c->vi != NULL) {
-						queue_del(&(c->vi));
-					}
-				}
-
-				tE1s[j] += Es[j];
-				tM1s[j] += abs(Ms[j]);
-				tE2s[j] += pow(Es[j], 2);
-				tM2s[j] += pow(Ms[j], 2);
-				tE4s[j] += pow(Es[j], 4);
-				tM4s[j] += pow(Ms[j], 4);
-
-				free(c);
-			}
-
-			for (uint32_t j = 0; j < n_temps - 1; j++) {
-				if (gsl_rng_uniform(r) < 0.5) {
-					if (gsl_rng_uniform(r) < exp((Es[j + 1] - Es[j]) * (1 / Ts[j + 1] - 1 / Ts[j]))) {
-						bool *tmp_x = xs[j];
-						double tmp_E = Es[j];
-						int32_t tmp_M = Ms[j];
-
-						xs[j] = xs[j + 1];
-						xs[j + 1] = tmp_x;
-
-						Es[j] = Es[j + 1];
-						Es[j + 1] = tmp_E;
-
-						Ms[j] = Ms[j + 1];
-						Ms[j + 1] = tmp_M;
-					}
-				}
-			}
-		}
-
-		for (uint32_t j = 0; j < n_temps; j++) {
-			tE1s[j] /= N;
-			tM1s[j] /= N;
-			tE2s[j] /= N;
-			tM2s[j] /= N;
-			tE4s[j] /= N;
-			tM4s[j] /= N;
-
-			if (n_runs > 0) {
-				E1s[j] = E1s[j] * ((n_runs - 1.) / n_runs) + tE1s[j] * 1. / n_runs;
-				M1s[j] = M1s[j] * ((n_runs - 1.) / n_runs) + tM1s[j] * 1. / n_runs;
-				E2s[j] = E2s[j] * ((n_runs - 1.) / n_runs) + tE2s[j] * 1. / n_runs;
-				M2s[j] = M2s[j] * ((n_runs - 1.) / n_runs) + tM2s[j] * 1. / n_runs;
-				E4s[j] = E4s[j] * ((n_runs - 1.) / n_runs) + tE4s[j] * 1. / n_runs;
-				M4s[j] = M4s[j] * ((n_runs - 1.) / n_runs) + tM4s[j] * 1. / n_runs;
-			}
-		}
-
-		if (n_runs > 0) {
-			diff = 0;
-			for (uint32_t j = 0; j < n_temps; j++) {
-				double dd = sqrt(M2s[j] - pow(M1s[j], 2)) / sqrt(N * n_runs);
-				double t_diff = dd / M1s[j];
-				if (t_diff > diff) diff = t_diff;
-			}
-		} else {
-			diff = 1e30;
-		}
-
-		n_runs++;
-	}
-
-	FILE *outfile = fopen("out.dat", "a");
-
-	for (uint32_t j = 0; j < n_temps; j++) {
-		double C = (E2s[j] - pow(E1s[j], 2)) / pow(Ts[j], 2);
-		double X = (M2s[j] - pow(M1s[j], 2)) / Ts[j];
-
-		double dE1 = sqrt(E2s[j] - pow(E1s[j], 2)) / sqrt(N * n_runs);
-		double dM1 = sqrt(M2s[j] - pow(M1s[j], 2)) / sqrt(N * n_runs);
-		double dE2 = sqrt(E4s[j] - pow(E2s[j], 2)) / sqrt(N * n_runs);
-		double dM2 = sqrt(M4s[j] - pow(M2s[j], 2)) / sqrt(N * n_runs);
-
-		double dC = sqrt(pow(dE2, 2) + pow(2 * E1s[j] * dE1, 2)) / pow(Ts[j], 2);
-		double dX = sqrt(pow(dM2, 2) + pow(2 * M1s[j] * dM1, 2)) / Ts[j];
-
-		fprintf(outfile, "%u %u %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", L, use_scho, Tins[j], Hins[j], E1s[j] / g->nv, dE1 / g->nv, M1s[j] / g->nv, dM1 / g->nv, C / g->nv, dC / g->nv, X / g->nv, dX / g->nv);
-	}
-
-	fclose(outfile);
-
-	for (uint32_t i = 0; i < n_temps; i++) {
-		free(xs[i]);
-	}
-
-	free(xs);
-
-	gsl_rng_free(r);
-	graph_free(g);
-
-	free(Ts);
-	free(Hs);
-	free(Es);
-	free(Ms);
-
-	free(E1s);
-	free(E2s);
-	free(E4s);
-
-	free(M1s);
-	free(M2s);
-	free(M4s);
-
-	free(tE1s);
-	free(tE2s);
-	free(tE4s);
-
-	free(tM1s);
-	free(tM2s);
-	free(tM4s);
-
-
-	return 0;
-}
-
diff --git a/src/wolff.c b/src/wolff.c
index 0c369b7..7675889 100644
--- a/src/wolff.c
+++ b/src/wolff.c
@@ -1,5 +1,5 @@
 
-#include "wolff.h"
+#include <wolff.h>
 
 int main(int argc, char *argv[]) {
   int opt;
@@ -22,55 +22,56 @@ int main(int argc, char *argv[]) {
 
   while ((opt = getopt(argc, argv, "N:L:T:H:m:e:oq:D")) != -1) {
     switch (opt) {
-      case 'N':
-        N = (uint64_t)atof(optarg);
+    case 'N':
+      N = (uint64_t)atof(optarg);
+      break;
+    case 'L':
+      L = atoi(optarg);
+      break;
+    case 'T':
+      T = atof(optarg);
+      break;
+    case 'H':
+      H = atof(optarg);
+      break;
+    case 'm':
+      min_runs = atoi(optarg);
+      break;
+    case 'e':
+      eps = atof(optarg);
+      break;
+    case 'o':
+      output_state = true;
+      break;
+    case 'D':
+      use_dual = true;
+      break;
+    case 'q':
+      lattice_i = atoi(optarg);
+      switch (lattice_i) {
+      case 0:
+        lat = SQUARE_LATTICE;
         break;
-      case 'L':
-        L = atoi(optarg);
+      case 1:
+        lat = DIAGONAL_LATTICE;
         break;
-      case 'T':
-        T = atof(optarg);
+      case 2:
+        lat = TRIANGULAR_LATTICE;
         break;
-      case 'H':
-        H = atof(optarg);
+      case 3:
+        lat = VORONOI_HYPERUNIFORM_LATTICE;
         break;
-      case 'm':
-        min_runs = atoi(optarg);
+      case 4:
+        lat = VORONOI_LATTICE;
         break;
-      case 'e':
-        eps= atof(optarg);
-        break;
-      case 'o':
-        output_state = true;
-        break;
-      case 'D':
-        use_dual = true;
-        break;
-			case 'q':
-				lattice_i = atoi(optarg);
-				switch (lattice_i) {
-					case 0:
-            lat = SQUARE_LATTICE;
-						break;
-					case 1:
-						lat = DIAGONAL_LATTICE;
-						break;
-					case 2:
-						lat = TRIANGULAR_LATTICE;
-						break;
-					case 3:
-						lat = VORONOI_HYPERUNIFORM_LATTICE;
-						break;
-          case 4:
-						lat = VORONOI_LATTICE;
-            break;
-					default:
-						printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 (DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n");
-						exit(EXIT_FAILURE);
-				}
-				break;
       default:
+        printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 "
+               "(DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n");
         exit(EXIT_FAILURE);
+      }
+      break;
+    default:
+      exit(EXIT_FAILURE);
     }
   }
 
@@ -94,12 +95,26 @@ int main(int argc, char *argv[]) {
   double E1, E2, dE1, M1, M2, dM1, C, dC, X, dX, Mmu2, Mmu4, Emu2, Emu4;
   double clust_per_sweep = 0;
 
-  E1 = 0; E2 = 0; M1 = 0; M2 = 0; C = 0; dC = 0; X = 0; dX = 0;
-  dE1 = 0; dM1 = 0; Mmu2 = 0; Mmu4 = 0; Emu2 = 0; Emu4 = 0;
+  E1 = 0;
+  E2 = 0;
+  M1 = 0;
+  M2 = 0;
+  C = 0;
+  dC = 0;
+  X = 0;
+  dX = 0;
+  dE1 = 0;
+  dM1 = 0;
+  Mmu2 = 0;
+  Mmu4 = 0;
+  Emu2 = 0;
+  Emu4 = 0;
 
   printf("\n");
   while (diff > eps || diff == 0. || n_runs < min_runs) {
-    printf("\033[F\033[JWOLFF: sweep %lu, dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n", n_runs, fabs(dE1 / E1), dM1 / M1, dC / C, dX / X, clust_per_sweep);
+    printf("\033[F\033[JWOLFF: sweep %" PRIu64
+           ", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n",
+           n_runs, fabs(dE1 / E1), dM1 / M1, dC / C, dX / X, clust_per_sweep);
 
     uint32_t n_flips = 0;
     uint32_t n_clust = 0;
@@ -114,19 +129,27 @@ int main(int argc, char *argv[]) {
     E2 = E2 * (n_runs / (n_runs + 1.)) + pow(s->H, 2) * 1. / (n_runs + 1.);
     M2 = M2 * (n_runs / (n_runs + 1.)) + pow(s->M, 2) * 1. / (n_runs + 1.);
 
-    Mmu2 = Mmu2 * (n_runs / (n_runs + 1.)) + pow(s->M - M1, 2) * 1. / (n_runs + 1.);
-    Mmu4 = Mmu4 * (n_runs / (n_runs + 1.)) + pow(s->M - M1, 4) * 1. / (n_runs + 1.);
-    Emu2 = Emu2 * (n_runs / (n_runs + 1.)) + pow(s->H - E1, 2) * 1. / (n_runs + 1.);
-    Emu4 = Emu4 * (n_runs / (n_runs + 1.)) + pow(s->H - E1, 4) * 1. / (n_runs + 1.);
+    Mmu2 = Mmu2 * (n_runs / (n_runs + 1.)) +
+           pow(s->M - M1, 2) * 1. / (n_runs + 1.);
+    Mmu4 = Mmu4 * (n_runs / (n_runs + 1.)) +
+           pow(s->M - M1, 4) * 1. / (n_runs + 1.);
+    Emu2 = Emu2 * (n_runs / (n_runs + 1.)) +
+           pow(s->H - E1, 2) * 1. / (n_runs + 1.);
+    Emu4 = Emu4 * (n_runs / (n_runs + 1.)) +
+           pow(s->H - E1, 4) * 1. / (n_runs + 1.);
 
     if (n_runs > 1) {
       double Msigma2 = n_runs / (n_runs - 1) * (M2 - pow(M1, 2));
       X = Msigma2 / T;
-      dX = sqrt((Mmu4 - (n_runs - 3.) / (n_runs - 1.) * pow(Mmu2, 2)) / n_runs) / T;
+      dX =
+          sqrt((Mmu4 - (n_runs - 3.) / (n_runs - 1.) * pow(Mmu2, 2)) / n_runs) /
+          T;
 
       double Esigma2 = n_runs / (n_runs - 1) * (E2 - pow(E1, 2));
       C = Esigma2 / T;
-      dC = sqrt((Emu4 - (n_runs - 3.) / (n_runs - 1.) * pow(Emu2, 2)) / n_runs) / T;
+      dC =
+          sqrt((Emu4 - (n_runs - 3.) / (n_runs - 1.) * pow(Emu2, 2)) / n_runs) /
+          T;
 
       dE1 = sqrt(Esigma2 / n_runs);
       dM1 = sqrt(Msigma2 / n_runs);
@@ -134,14 +157,20 @@ int main(int argc, char *argv[]) {
       diff = fabs(dX / X);
     }
 
-    clust_per_sweep = clust_per_sweep * (n_runs / (n_runs + 1.)) + (n_clust * 1. / N) * 1. / (n_runs + 1.);
+    clust_per_sweep = clust_per_sweep * (n_runs / (n_runs + 1.)) +
+                      (n_clust * 1. / N) * 1. / (n_runs + 1.);
 
     n_runs++;
   }
-  printf("\033[F\033[JWOLFF: sweep %lu, dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n", n_runs, fabs(dE1 / E1), dM1 / M1, dC / C, dX / X, clust_per_sweep);
+  printf("\033[F\033[JWOLFF: sweep %" PRIu64
+         ", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n",
+         n_runs, fabs(dE1 / E1), dM1 / M1, dC / C, dX / X, clust_per_sweep);
 
   FILE *outfile = fopen("out.dat", "a");
-  fprintf(outfile, "%u %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", L, T, H, E1 / h->nv, dE1 / h->nv, M1 / h->nv, dM1 / h->nv, C / h->nv, dC / h->nv, X / h->nv, dX / h->nv);
+  fprintf(outfile,
+          "%u %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", L,
+          T, H, E1 / h->nv, dE1 / h->nv, M1 / h->nv, dM1 / h->nv, C / h->nv,
+          dC / h->nv, X / h->nv, dX / h->nv);
   fclose(outfile);
 
   free(bond_probs);
@@ -164,4 +193,3 @@ int main(int argc, char *argv[]) {
 
   return 0;
 }
-
diff --git a/src/wolff.h b/src/wolff.h
deleted file mode 100644
index cec9ee3..0000000
--- a/src/wolff.h
+++ /dev/null
@@ -1,52 +0,0 @@
-
-#include <string.h>
-#include <math.h>
-#include <getopt.h>
-#include <float.h>
-#include <sys/types.h>
-#include <inttypes.h>
-#include <gsl/gsl_randist.h>
-#include <gsl/gsl_rng.h>
-#include <stdbool.h>
-#include <assert.h>
-#include <fftw3.h>
-
-#include <jst/graph.h>
-#include <jst/rand.h>
-
-typedef struct {
-  graph_t *g;
-  bool *spins;
-  int32_t M;
-  double H;
-} ising_state_t;
-
-typedef struct ll_tag {
-	uint32_t x;
-	struct ll_tag *next;
-} ll_t;
-
-typedef struct {
-  uint32_t nv;
-  double dH;
-  int32_t dM;
-} cluster_t;
-
-double get_hamiltonian(graph_t *g, double *coupling, bool *x);
-
-void stack_push(ll_t **q, uint32_t x);
-
-uint32_t stack_pop(ll_t **q);
-
-bool stack_contains(const ll_t *q, uint32_t x);
-
-cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, gsl_rng *r);
-
-graph_t *graph_add_ext(const graph_t *g);
-
-double hh(double th);
-
-double *get_bond_probs(double T, double H, ising_state_t *s);
-
-uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps);
-