started wolff code in new repository

8 files changed, 793 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..a007fea
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+build/*
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..c325cb4
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,21 @@
+
+cmake_minimum_required(VERSION 3.7)
+project(wolff)
+
+include_directories(src ~/.local/include)
+link_directories(~/.local/lib)
+
+file(GLOB SOURCES lib/*.c)
+add_executable(wolff src/wolff.c ${SOURCES})
+
+find_package(OpenMP)
+if (OPENMP_FOUND)
+    set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+    set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+endif()
+
+target_link_libraries(wolff gsl m jst cblas fftw3 tcmalloc profiler)
+#target_link_libraries(pt_wolff gsl m jst cblas tcmalloc profiler)
+
+install(TARGETS wolff DESTINATION bin)
+
diff --git a/lib/queue.c b/lib/queue.c
new file mode 100644
index 0000000..9a01741
--- /dev/null
+++ b/lib/queue.c
@@ -0,0 +1,32 @@
+
+#include "wolff.h"
+
+void stack_push(ll_t **q, uint32_t x) {
+	ll_t *nq = malloc(sizeof(ll_t));
+	nq->x = x;
+	nq->next = *q;
+
+	*q = nq;
+}
+
+uint32_t stack_pop(ll_t **q) {
+	ll_t *old_q = *q;
+
+	*q = old_q->next;
+	uint32_t x = old_q->x;
+
+	free(old_q);
+
+	return x;
+}
+
+bool stack_contains(const ll_t *q, uint32_t x) {
+  if (q == NULL) {
+    return false;
+  } else if (q->x == x) {
+    return true;
+  } else {
+    return stack_contains(q->next, x);
+  }
+}
+
diff --git a/lib/wolff.h b/lib/wolff.h
new file mode 100644
index 0000000..6beb03e
--- /dev/null
+++ b/lib/wolff.h
@@ -0,0 +1,51 @@
+
+#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 <omp.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 {
+  int32_t nv;
+  double dH;
+} 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);
+
+int32_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
new file mode 100644
index 0000000..72cf9f7
--- /dev/null
+++ b/lib/wolff_tools.c
@@ -0,0 +1,179 @@
+
+#include "wolff.h"
+
+graph_t *graph_add_ext(const graph_t *g) {
+  graph_t *h = (graph_t *)calloc(1, sizeof(graph_t));
+
+  h->nv = g->nv + 1;
+  h->ne = g->ne + g->nv;
+
+  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->vx = (double *)malloc(2 * h->nv * sizeof(double));
+  h->bq = (bool *)malloc(h->nv * sizeof(bool));
+
+  memcpy(h->ev, g->ev, 2 * g->ne * sizeof(uint32_t));
+  memcpy(h->vx, g->vx, 2 * g->nv * sizeof(double));
+  memcpy(h->bq, g->bq, g->nv * sizeof(bool));
+  h->vx[2 * g->nv] = -1;
+  h->vx[2 * g->nv + 1] = -0.5;
+  h->bq[g->nv] = false;
+
+  for (uint32_t i = 0; i < g->nv; i++) {
+    h->ev[2 * g->ne + 2 * i] = i;
+    h->ev[2 * g->ne + 2 * i + 1] = g->nv;
+  }
+
+  for (uint32_t i = 0; i < g->nv; i++) {
+    h->vei[i] = g->vei[i] + i;
+
+    for (uint32_t j = 0; j < g->vei[i + 1] - g->vei[i]; j++) {
+      h->ve[h->vei[i] + j] = g->ve[g->vei[i] + j];
+    }
+
+    h->ve[h->vei[i] + g->vei[i + 1] - g->vei[i]] = g->ne + i;
+  }
+
+  h->vei[g->nv] = g->vei[g->nv] + g->nv;
+  h->vei[g->nv + 1] = h->vei[g->nv] + g->nv;
+
+  for (uint32_t i = 0; i < g->nv; i++) {
+    h->ve[h->vei[g->nv] + i] = g->ne + i;
+  }
+
+  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) {
+  uint32_t v0;
+  bool x0;
+  cluster_t *c;
+  
+  v0 = gsl_rng_uniform_int(r, g->nv); // pick a random vertex
+  x0 = x[v0]; // record its orientation
+
+  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;
+  c->dH = 0;
+
+  while (stack != NULL) {
+    uint32_t v;
+    uint16_t nn;
+
+    v = stack_pop(&stack);
+    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
+
+      for (uint16_t i = 0; i < nn; i++) {
+        uint32_t e, v1, v2, vn;
+
+        e = g->ve[g->vei[v] + i]; // select the ith bond connected to site
+        v1 = g->ev[2 * e];
+        v2 = g->ev[2 * e + 1];
+
+        vn = v == v1 ? v2 : v1; // distinguish neighboring site from site itself
+
+        if (x[vn] == x0) { // if the neighboring site matches the flipping cluster...
+          if (v1 == g->nv - 1 || v2 == g->nv - 1) {
+            c->dH += H;
+          } else {
+            c->dH += 1;
+          }
+
+          if (gsl_rng_uniform(r) < ps[e]) { // and with probability ps[e]...
+            stack_push(&stack, vn); // push the neighboring vertex to the stack
+          }
+        } else {
+          if (v1 == g->nv - 1 || v2 == g->nv - 1) {
+            c->dH -= H;
+          } else {
+            c->dH -= 1;
+          }
+        }
+      }
+
+      if (v != g->nv - 1) {
+        c->nv++;
+      }
+    }
+  }
+
+  if (x0) {
+    c->nv = -c->nv;
+  }
+
+  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 * H / T);
+
+  double *ps = (double *)malloc(s->g->ne * sizeof(double));
+
+  for (uint32_t i = 0; i < s->g->ne; i++) {
+    uint32_t v1, v2;
+    v1 = s->g->ev[2 * i];
+    v2 = s->g->ev[2 * i + 1];
+    if (v1 == s->g->nv - 1 || v2 == s->g->nv - 1) {
+      ps[i] = q;
+    } else {
+      ps[i] = p;
+    }
+  }
+
+  return ps;
+}
+
+int32_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());
+  }
+
+  if (ps == NULL) {
+    ps = get_bond_probs(T, H, s);
+  }
+
+  cluster_t *c = flip_cluster(s->g, ps, H, s->spins, r);
+
+  s->M += 2 * c->nv;
+  s->H += c->dH;
+
+  int32_t n_flipped = c->nv;
+
+  free(c);
+
+  return n_flipped;
+}
+
diff --git a/src/pt_wolff.c b/src/pt_wolff.c
new file mode 100644
index 0000000..128b601
--- /dev/null
+++ b/src/pt_wolff.c
@@ -0,0 +1,317 @@
+
+#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
new file mode 100644
index 0000000..80f8278
--- /dev/null
+++ b/src/wolff.c
@@ -0,0 +1,140 @@
+
+#include "wolff.h"
+
+int main(int argc, char *argv[]) {
+  int opt;
+  bool output_state;
+  lattice_t lat;
+  uint16_t L;
+  uint32_t min_runs;
+  uint64_t N;
+  double T, H, eps;
+
+  L = 128;
+  N = 1000;
+  lat = SQUARE_LATTICE;
+  T = 2.3;
+  H = 0;
+  eps = 1e30;
+  output_state = false;
+  min_runs = 10;
+
+  while ((opt = getopt(argc, argv, "N:L:T:H:m:e:o")) != -1) {
+    switch (opt) {
+      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;
+      default:
+        exit(EXIT_FAILURE);
+    }
+  }
+
+  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+  gsl_rng_set(r, jst_rand_seed());
+
+  graph_t *h = graph_create(lat, TORUS_BOUND, L, false);
+
+  ising_state_t *s = (ising_state_t *)calloc(1, sizeof(ising_state_t));
+
+  s->g = graph_add_ext(h);
+  s->spins = (bool *)calloc(h->nv + 1, sizeof(bool));
+  s->M = -h->nv;
+  s->H = -(1.0 * h->ne) - H * h->nv;
+
+  double *bond_probs = get_bond_probs(T, H, s);
+
+  double diff = 1e31;
+  uint64_t n_runs = 0;
+
+  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;
+
+  printf("\n");
+  while (diff > eps || diff == 0. || diff != diff || n_runs < min_runs) {
+    printf("\033[F\033[JWOLFF: sweep %llu, 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;
+
+    while (n_flips / h->nv < N) {
+      n_flips += abs(wolff_step(T, H, s, r, bond_probs));
+      n_clust++;
+    }
+
+    E1 = E1 * (n_runs / (n_runs + 1.)) + s->H * 1. / (n_runs + 1.);
+    M1 = M1 * (n_runs / (n_runs + 1.)) + abs(s->M) * 1. / (n_runs + 1.);
+    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(abs(s->M) - M1, 2) * 1. / (n_runs + 1.);
+    Mmu4 = Mmu4 * (n_runs / (n_runs + 1.)) + pow(abs(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;
+
+      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;
+
+      dE1 = sqrt(Esigma2 / n_runs);
+      dM1 = sqrt(Msigma2 / n_runs);
+
+      diff = fabs(dX / X);
+    }
+
+    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 %llu, 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);
+  fclose(outfile);
+
+  free(bond_probs);
+
+  if (output_state) {
+    FILE *state_file = fopen("state.dat", "a");
+    for (uint32_t i = 0; i < h->nv; i++) {
+      fprintf(state_file, "%d ", s->spins[i]);
+    }
+    fprintf(state_file, "\n");
+    fclose(state_file);
+  }
+
+  gsl_rng_free(r);
+  graph_free(s->g);
+  free(s->spins);
+  free(s);
+  free(bond_probs);
+  graph_free(h);
+
+  return 0;
+}
+
diff --git a/src/wolff.h b/src/wolff.h
new file mode 100644
index 0000000..9619a4d
--- /dev/null
+++ b/src/wolff.h
@@ -0,0 +1,52 @@
+
+#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 <omp.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 {
+  int32_t nv;
+  double dH;
+} 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);
+
+int32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps);
+