does potts now, no external libraries

14 files changed, 431 insertions, 511 deletions

diff --git a/lib/convex.c b/lib/convex.c
index 096504d..df506f9 100644
--- a/lib/convex.c
+++ b/lib/convex.c
@@ -5,7 +5,7 @@ double slope(point_t *P, point_t *Q) {
   return (Q->y - P->y) / ((double)(Q->x) - (double)(P->x));
 }
 
-double *get_convex_minorant(uint64_t n, double *Gammas) {
+double *get_convex_minorant(count_t n, double *Gammas) {
   if (n < 2) {
     return Gammas;
   }
@@ -17,7 +17,7 @@ double *get_convex_minorant(uint64_t n, double *Gammas) {
 
   list_t *pos = L;
 
-  for (uint64_t i = 1; i < n; i++) {
+  for (count_t i = 1; i < n; i++) {
     pos->next = (list_t *)calloc(1, sizeof(list_t));
     pos->next->p = (point_t *)calloc(1, sizeof(point_t));
     pos->next->p->x = i;
@@ -69,7 +69,7 @@ double *get_convex_minorant(uint64_t n, double *Gammas) {
   double *g = (double *)calloc(n + 1, sizeof(double));
   double rho = 0;
 
-  for (uint64_t i = 0; i < n + 1; i++) {
+  for (count_t i = 0; i < n + 1; i++) {
     if (i > pos->next->p->x) {
       pos = pos->next;
     }
@@ -87,7 +87,7 @@ double *get_convex_minorant(uint64_t n, double *Gammas) {
     }
   }
 
-  for (uint64_t i = 0; i < n + 1; i++) {
+  for (count_t i = 0; i < n + 1; i++) {
     g[i] += rho / 2;
   }
 
diff --git a/lib/convex.h b/lib/convex.h
index 29f4ae0..5a405d4 100644
--- a/lib/convex.h
+++ b/lib/convex.h
@@ -6,8 +6,10 @@
 #include <stdlib.h>
 #include <string.h>
 
+#include "types.h"
+
 typedef struct {
-  uint64_t x;
+  count_t x;
   double y;
 } point_t;
 
@@ -17,5 +19,5 @@ typedef struct list_tag {
   point_t *p;
 } list_t;
 
-double *get_convex_minorant(uint64_t n, double *Gammas);
+double *get_convex_minorant(count_t n, double *Gammas);
 
diff --git a/lib/graph.c b/lib/graph.c
new file mode 100644
index 0000000..a1edf13
--- /dev/null
+++ b/lib/graph.c
@@ -0,0 +1,67 @@
+
+#include "graph.h"
+
+graph_t *graph_create_square(D_t D, L_t L) {
+  v_t nv = pow(L, D);
+  v_t ne = D * nv;
+
+  v_t *v_i = (v_t *)malloc((nv + 1) * sizeof(v_t));
+
+  for (v_t i = 0; i < nv + 1; i++) {
+    v_i[i] = 2 * D * i;
+  }
+
+  v_t *v_adj = (v_t *)malloc(2 * D * nv * sizeof(v_t));
+
+  for (v_t i = 0; i < nv; i++) {
+    for (D_t j = 0; j < D; j++) {
+      v_adj[v_i[i] + 2 * j] = pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1));
+      v_adj[v_i[i] + 2 * j + 1] = pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1));
+    }
+  }
+
+  graph_t *g = (graph_t *)malloc(sizeof(graph_t));
+
+  g->ne = ne;
+  g->nv = nv;
+  g->v_i = v_i;
+  g->v_adj = v_adj;
+
+  return g;
+}
+
+graph_t *graph_add_ext(const graph_t *G) {
+  graph_t *tG = (graph_t *)calloc(1, sizeof(graph_t));
+
+  tG->nv = G->nv + 1;
+  tG->ne = G->ne + G->nv;
+
+  tG->v_i = (v_t *)malloc((tG->nv + 1) * sizeof(v_t));
+  tG->v_adj = (v_t *)malloc(2 * tG->ne * sizeof(v_t));
+
+  for (v_t i = 0; i < G->nv + 1; i++) {
+    tG->v_i[i] = G->v_i[i] + i;
+  }
+
+  tG->v_i[tG->nv] = 2 * tG->ne;
+
+  for (v_t i = 0; i < G->nv; i++) {
+    v_t nn = G->v_i[i + 1] - G->v_i[i];
+
+    for (v_t j = 0; j < nn; j++) {
+      tG->v_adj[tG->v_i[i] + j] = G->v_adj[G->v_i[i] + j];
+    }
+
+    tG->v_adj[tG->v_i[i] + nn] = G->nv;
+    tG->v_adj[tG->v_i[G->nv] + i] = i;
+  }
+
+  return tG;
+}
+
+void graph_free(graph_t *g) {
+  free(g->v_i);
+  free(g->v_adj);
+  free(g);
+}
+
diff --git a/lib/graph.h b/lib/graph.h
new file mode 100644
index 0000000..eda7630
--- /dev/null
+++ b/lib/graph.h
@@ -0,0 +1,18 @@
+
+#include <inttypes.h>
+#include <math.h>
+#include <stdlib.h>
+
+#include "types.h"
+
+typedef struct {
+  v_t ne;
+  v_t nv;
+  v_t *v_i;
+  v_t *v_adj;
+} graph_t;
+
+graph_t *graph_create_square(D_t D, L_t L);
+
+graph_t *graph_add_ext(const graph_t *G);
+
diff --git a/lib/rand.c b/lib/rand.c
new file mode 100644
index 0000000..76f537d
--- /dev/null
+++ b/lib/rand.c
@@ -0,0 +1,20 @@
+
+#include <assert.h>
+#include <stdio.h>
+
+unsigned long int rand_seed() {
+  /*
+    returns a random unsigned long integer read from the standard unix
+    pseudorandom device /dev/urandom
+  */
+
+  FILE *f = fopen("/dev/urandom", "r");
+  assert(f != NULL);
+
+  unsigned long int seed;
+  fread(&seed, sizeof(unsigned long int), 1, f);
+
+  fclose(f);
+
+  return seed;
+}
diff --git a/lib/rand.h b/lib/rand.h
new file mode 100644
index 0000000..2354f6a
--- /dev/null
+++ b/lib/rand.h
@@ -0,0 +1,7 @@
+
+#pragma once
+
+#include <assert.h>
+#include <stdio.h>
+
+unsigned long int rand_seed();
diff --git a/lib/queue.c b/lib/stack.c
index 0823518..9fde1e6 100644
--- a/lib/queue.c
+++ b/lib/stack.c
@@ -1,7 +1,7 @@
 
-#include "queue.h"
+#include "stack.h"
 
-void stack_push(ll_t **q, uint32_t x) {
+void stack_push(ll_t **q, v_t x) {
   ll_t *nq = malloc(sizeof(ll_t));
   nq->x = x;
   nq->next = *q;
@@ -17,11 +17,11 @@ void stack_push_d(dll_t **q, double x) {
   *q = nq;
 }
 
-uint32_t stack_pop(ll_t **q) {
+v_t stack_pop(ll_t **q) {
   ll_t *old_q = *q;
 
   *q = old_q->next;
-  uint32_t x = old_q->x;
+  v_t x = old_q->x;
 
   free(old_q);
 
@@ -39,12 +39,3 @@ double stack_pop_d(dll_t **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/queue.h b/lib/stack.h
index 23834ba..a354ab5 100644
--- a/lib/queue.h
+++ b/lib/stack.h
@@ -6,8 +6,10 @@
 #include <stdlib.h>
 #include <string.h>
 
+#include "types.h"
+
 typedef struct ll_tag {
-  uint32_t x;
+  v_t x;
   struct ll_tag *next;
 } ll_t;
 
@@ -16,10 +18,9 @@ typedef struct dll_tag {
   struct dll_tag *next;
 } dll_t;
 
-void stack_push(ll_t **q, uint32_t x);
+void stack_push(ll_t **q, v_t x);
 void stack_push_d(dll_t **q, double x);
 
-uint32_t stack_pop(ll_t **q);
+v_t stack_pop(ll_t **q);
 double stack_pop_d(dll_t **q);
 
-bool stack_contains(const ll_t *q, uint32_t x);
diff --git a/lib/tree.c b/lib/tree.c
new file mode 100644
index 0000000..081f1d1
--- /dev/null
+++ b/lib/tree.c
@@ -0,0 +1,81 @@
+
+#include "tree.h"
+
+node_t *tree_createNode(v_t value, v_t level, node_t *left, node_t *right) {
+  node_t *node = (node_t *)malloc(sizeof(node_t));
+
+  node->value = value;
+  node->level = level;
+  node->left = left;
+  node->right = right;
+
+  return node;
+}
+
+void tree_freeNode(node_t *node) {
+  if (node != NULL) {
+    tree_freeNode(node->left);
+    tree_freeNode(node->right);
+    free(node);
+  }
+}
+
+void tree_skew(node_t **T) {
+  if (*T != NULL) {
+    if ((*T)->left != NULL) {
+      if ((*T)->left->level == (*T)->level) {
+        node_t *L = (*T)->left;
+        (*T)->left = L->right;
+        L->right = (*T);
+
+        *T = L;
+      }
+    }
+  }
+}
+
+void tree_split(node_t **T) {
+  if (*T != NULL) {
+    if ((*T)->right != NULL) {
+      if ((*T)->right->right != NULL) {
+        if ((*T)->level == (*T)->right->right->level) {
+          node_t *R = (*T)->right;
+          (*T)->right = R->left;
+          R->left = *T;
+          R->level = R->level + 1;
+          *T = R;
+        }
+      }
+    }
+  }
+}
+
+void tree_insert(node_t **T, v_t x) {
+  if (*T == NULL) {
+    *T = tree_createNode(x, 1, NULL, NULL);
+  } else if (x < (*T)->value) {
+    node_t *L = (*T)->left;
+    tree_insert(&L, x);
+    (*T)->left = L;
+  } else if (x > (*T)->value) {
+    node_t *R = (*T)->right;
+    tree_insert(&R, x);
+    (*T)->right = R;
+  }
+
+  tree_skew(T);
+  tree_split(T);
+}
+
+bool tree_contains(node_t *T, v_t x) {
+  if (T == NULL) {
+    return false;
+  } else if (x < T->value) {
+    return tree_contains(T->left, x);
+  } else if (x > T->value) {
+    return tree_contains(T->right, x);
+  } else {
+    return true;
+  }
+}
+
diff --git a/lib/tree.h b/lib/tree.h
new file mode 100644
index 0000000..dc22c2d
--- /dev/null
+++ b/lib/tree.h
@@ -0,0 +1,22 @@
+
+#pragma once
+
+#include <inttypes.h>
+#include <stdlib.h>
+#include <stdbool.h>
+
+#include "types.h"
+
+typedef struct node_t {
+  v_t value;
+  v_t level;
+  struct node_t *left;
+  struct node_t *right;
+} node_t;
+
+void tree_insert(node_t **T, v_t x);
+
+bool tree_contains(node_t *T, v_t x);
+
+void tree_freeNode(node_t *T);
+
diff --git a/lib/types.h b/lib/types.h
new file mode 100644
index 0000000..982db47
--- /dev/null
+++ b/lib/types.h
@@ -0,0 +1,16 @@
+
+#include <inttypes.h>
+
+typedef uint_fast32_t v_t;
+typedef uint_fast8_t q_t;
+typedef uint_fast8_t D_t;
+typedef uint_fast16_t L_t;
+typedef uint_fast64_t count_t;
+
+#define MAX_Q 256
+
+#define PRID PRIu8
+#define PRIL PRIu16
+#define PRIcount PRIu64
+#define PRIq PRIu8
+
diff --git a/lib/wolff.h b/lib/wolff.h
index dfe7b1c..dae907a 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -13,34 +13,26 @@
 #include <string.h>
 #include <sys/types.h>
 
-#include <jst/graph.h>
-#include <jst/rand.h>
-
-#include "queue.h"
+#include "types.h"
+#include "rand.h"
+#include "stack.h"
 #include "convex.h"
-
-typedef enum {
-  WOLFF,
-  WOLFF_GHOST,
-  METROPOLIS
-} sim_t;
+#include "graph.h"
+#include "tree.h"
 
 typedef struct {
   graph_t *g;
-  bool *spins;
-  int32_t M;
-  double H;
+  q_t *spins;
+  double T;
+  double *H;
+  double T_prob;
+  double *H_probs;
+  double E;
+  v_t *M;
+  q_t q;
 } ising_state_t;
 
 typedef struct {
-  uint32_t nv;
-  int32_t dJb;
-  int32_t dHb;
-  bool hit_ghost;
-  ll_t *spins;
-} cluster_t;
-
-typedef struct {
   uint64_t n;
   double x;
   double dx;
@@ -60,21 +52,13 @@ typedef struct {
   double O2;
 } autocorr_t;
 
-int8_t sign(double x);
-
-cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_on_ghost,
-                        gsl_rng *r);
+v_t flip_cluster(ising_state_t *s, v_t v0, q_t s1, gsl_rng *r);
 
 graph_t *graph_add_ext(const graph_t *g);
 
-uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r,
-                    double *ps);
-
 void update_meas(meas_t *m, double x);
 
 void update_autocorr(autocorr_t *OO, double O);
 
-double add_to_avg(double mx, double x, uint64_t n);
-
 double rho(autocorr_t *o, uint64_t i);
 
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index a417d4e..416a498 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -1,257 +1,94 @@
 
-#include "queue.h"
 #include "wolff.h"
 
-int8_t spin_to_sign(bool spin) {
-  /* takes a spin (represented by a bool) and converts it to an integer (plus
-   * or minus one). our convention takes true to be negative spin and false to
-   * be positive spin
-   */
-  if (spin) {
-    return -1;
-  } else {
-    return 1;
-  }
-}
-
-int8_t sign(double x) {
-  // the sign function. zero returns positive sign.
-  return x >= 0 ? 1 : -1;
-}
-
-graph_t *graph_add_ext(const graph_t *G) {
-  /* takes a graph object G and returns tG, the same graph with an extra vertex
-   * that is connected to every other vertex.
-   */
-  graph_t *tG = (graph_t *)calloc(1, sizeof(graph_t));
-
-  tG->nv = G->nv + 1;
-  tG->ne = G->ne + G->nv;
-
-  tG->ev = (uint32_t *)malloc(2 * tG->ne * sizeof(uint32_t));
-  tG->vei = (uint32_t *)malloc((tG->nv + 1) * sizeof(uint32_t));
-  tG->ve = (uint32_t *)malloc(2 * tG->ne * sizeof(uint32_t));
-  tG->vx = (double *)malloc(2 * tG->nv * sizeof(double));
-  tG->bq = (bool *)malloc(tG->nv * sizeof(bool));
-
-  memcpy(tG->ev, G->ev, 2 * G->ne * sizeof(uint32_t));
-  memcpy(tG->vx, G->vx, 2 * G->nv * sizeof(double));
-  memcpy(tG->bq, G->bq, G->nv * sizeof(bool));
-
-  tG->vx[2 * G->nv] = -1;
-  tG->vx[2 * G->nv + 1] = -0.5;
-  tG->bq[G->nv] = false;
-
-  for (uint32_t i = 0; i < G->nv; i++) {
-    tG->ev[2 * G->ne + 2 * i] = i;
-    tG->ev[2 * G->ne + 2 * i + 1] = G->nv;
-  }
-
-  for (uint32_t i = 0; i < G->nv; i++) {
-    tG->vei[i] = G->vei[i] + i;
-
-    for (uint32_t j = 0; j < G->vei[i + 1] - G->vei[i]; j++) {
-      tG->ve[tG->vei[i] + j] = G->ve[G->vei[i] + j];
-    }
-
-    tG->ve[tG->vei[i] + G->vei[i + 1] - G->vei[i]] = G->ne + i;
-  }
-
-  tG->vei[G->nv] = G->vei[G->nv] + G->nv;
-  tG->vei[G->nv + 1] = tG->vei[G->nv] + G->nv;
-
-  for (uint32_t i = 0; i < G->nv; i++) {
-    tG->ve[tG->vei[G->nv] + i] = G->ne + i;
-  }
-
-  return tG;
-}
-
-uint32_t get_neighbor(const graph_t *g, uint32_t v, uint32_t i) {
-  // returns the index of the ith neighbor of vertex v on graph g.
-  assert(i < g->vei[v + 1] - g->vei[v]); // don't request a neighbor over the total number of neighbors!
-  
-  uint32_t e, v1, v2;
-
-  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];
-
-  return v == v1 ? v2 : v1; // distinguish neighboring site from site itself
+q_t delta(q_t s0, q_t s1) {
+  return s0 == s1 ? 1 : 0;
 }
 
-cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *s, bool stop_on_ghost,
-                        gsl_rng *r) {
-  /* flips a wolff cluster on the graph g with initial state s. s is always
-   * changed by the routine. the probability of adding a normal bond to the
-   * cluster is passed by ps[0], and the probability of adding a bond with the
-   * external field spin to the cluster is passed by ps[1]. if stop_on_ghost is
-   * true, adding the external field spin to the cluster stops execution of the
-   * routine. flip_cluster returns an object of type cluster_t, which encodes
-   * information about the flipped cluster.
-   */
-  uint32_t v0;
-  int32_t n_h_bonds, n_bonds;
-  bool s0;
-  cluster_t *c;
-
-  v0 = gsl_rng_uniform_int(r, g->nv); // pick a random vertex
-  s0 = s[v0];                         // record its orientation
+v_t flip_cluster(ising_state_t *s, v_t v0, q_t step, gsl_rng *r) {
+  q_t s0 = s->spins[v0];
+  v_t nv = 0;
 
   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));
+  node_t *T = NULL;
 
   while (stack != NULL) {
-    uint32_t v;
-    uint32_t nn;
-
-    v = stack_pop(&stack);
-    nn = g->vei[v + 1] - g->vei[v];
+    v_t v = stack_pop(&stack);
 
-    if (s[v] == s0 && !(c->hit_ghost && stop_on_ghost)) { // if the vertex hasn't already been flipped
-      s[v] = !s[v];   // flip the vertex
-      if (stop_on_ghost) {
-        stack_push(&(c->spins), v);
-      }
-
-      for (uint32_t i = 0; i < nn; i++) {
-        bool is_ext;
-        uint32_t vn;
-        int32_t *bond_counter;
-        double prob;
+    if (!tree_contains(T, v)) { // if the vertex hasn't already been flipped
+      q_t s_old = s->spins[v];
+      q_t s_new = (s->spins[v] + step) % s->q;
 
-        vn = get_neighbor(g, v, i);
+      s->spins[v] = s_new;   // flip the vertex
+      tree_insert(&T, v);
 
-        is_ext = (v == g->nv - 1 || vn == g->nv - 1); // our edge contained the "ghost" spin if either of its vertices was the last in the graph
+      v_t nn = s->g->v_i[v + 1] - s->g->v_i[v];
 
-        bond_counter = is_ext ? &(c->dHb) : &(c->dJb);
-        prob = is_ext ? ps[1] : ps[0];
+      for (v_t i = 0; i < nn; i++) {
+        v_t vn = s->g->v_adj[s->g->v_i[v] + i];
+        q_t sn = s->spins[vn];
+        double prob;
 
-        if (s[vn] ==
-            s0) { // if the neighboring site matches the flipping cluster...
-          (*bond_counter)++;
+        bool is_ext = (v == s->g->nv - 1 || vn == s->g->nv - 1);
 
-          if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
-            if (is_ext && stop_on_ghost) {
-              c->hit_ghost = true;
-            }
-            stack_push(&stack, vn); // push the neighboring vertex to the stack
+        if (is_ext) {
+            q_t M_ind_0;
+            q_t M_ind_1;
+          if (vn == s->g->nv - 1) {
+            M_ind_0 = (s_old + s->q - s->spins[s->g->nv - 1]) % s->q;
+            M_ind_1 = (s_new + s->q - s->spins[s->g->nv - 1]) % s->q;
+          } else {
+            M_ind_0 = (sn + s->q - s_old) % s->q;
+            M_ind_1 = (sn + s->q - s_new) % s->q;
           }
+          prob = s->H_probs[M_ind_1 * s->q + M_ind_0];
+          s->M[M_ind_0]--;
+          s->M[M_ind_1]++;
+          s->E += - s->H[M_ind_1] + s->H[M_ind_0];
         } else {
-          (*bond_counter)--;
+          prob = s->T_prob * delta(s0, sn);
+          s->E += - ((double)delta(s->spins[v], sn)) + ((double)delta(s0, sn));
         }
-      }
-
-      if (v != g->nv - 1) { // count the number of non-external sites that flip
-        c->nv++;
-      }
-    }
-  }
-
-  return c;
-}
-
-uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r,
-                    double *ps) {
-  uint32_t n_flips;
-  bool no_r, no_ps;
-  no_r = false;
-  no_ps = false;
-
-  if (r == NULL) {
-    no_r = true;
-    r = gsl_rng_alloc(gsl_rng_mt19937);
-    gsl_rng_set(r, jst_rand_seed());
-  }
-
-  if (ps == NULL) { // computing exponentials is relatively expensive, so calling functions are given the option to supply values calculated once in the entire runtime
-    no_ps = true;
-    ps = (double *)malloc(2 * sizeof(double));
-    ps[0] = 1 - exp(-2 / T);
-    ps[1] = 1 - exp(-2 * fabs(H) / T);
-  }
-
-  switch (sim) {
-    case METROPOLIS: {
-      uint32_t v0 = gsl_rng_uniform_int(r, s->g->nv); // pick a random vertex
-      uint32_t nn = s->g->vei[v0 + 1] - s->g->vei[v0];
-
-      double dE = 0;
-      for (uint32_t i = 0; i < nn; i++) {
-        bool is_ext;
-        uint32_t vn;
-        int32_t *bond_counter;
-        double prob;
 
-        vn = get_neighbor(s->g, v0, i);
-
-        is_ext = (v0 == s->g->nv - 1 || vn == s->g->nv - 1); // our edge contained the "ghost" spin if either of its vertices was the last in the graph
-
-        if (is_ext) {
-          dE += 2 * H * spin_to_sign(s->spins[vn]) * spin_to_sign(s->spins[v0]);
-        } else {
-          dE += 2 * spin_to_sign(s->spins[vn]) * spin_to_sign(s->spins[v0]);
+        if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
+          stack_push(&stack, vn); // push the neighboring vertex to the stack
         }
       }
 
-      double p = exp(-dE / T);
-      if (gsl_rng_uniform(r) < p) {
-        s->M += - sign(H) * 2 * spin_to_sign(s->spins[v0]) * spin_to_sign(s->spins[s->g->nv - 1]);
-        s->H += dE;
-        s->spins[v0] = !s->spins[v0];
+      if (v != s->g->nv - 1) { // count the number of non-external sites that flip
+        nv++;
       }
+    }
+  }
 
-      n_flips = 1;
-      } break;
-
-    case WOLFF: {
-      cluster_t *c = flip_cluster(s->g, ps, s->spins, false, r);
-      s->M += - sign(H) * 2 * c->dHb;
-      s->H += 2 * (c->dJb + sign (H) * H * c->dHb);
-      n_flips = c->nv;
+  tree_freeNode(T);
 
-      free(c);
-      } break;
-    case WOLFF_GHOST: {
-      cluster_t *c = flip_cluster(s->g, ps, s->spins, true, r);
+  return nv;
+}
 
-      if (c->hit_ghost) {
-        while (c->spins != NULL) {
-          uint32_t v = stack_pop(&(c->spins));
-          s->spins[v] = !s->spins[v];
-          // if we hit the external spin, undo the cluster flip
-        }
-      } else {
-        while (c->spins != NULL) {
-          stack_pop(&(c->spins));
-          // we have to clear the memory on the stack anyway...
-        }
-        s->M += - sign(H) * 2 * c->dHb;
-        s->H += 2 * (c->dJb + sign (H) * H * c->dHb);
-      }
+double add_to_avg(double mx, double x, count_t n) {
+  return mx * (n / (n + 1.0)) + x * 1.0 / (n + 1.0);
+}
 
-      n_flips = c->nv;
+void update_meas(meas_t *m, double x) {
+  count_t n = m->n;
 
-      free(c);
-      } break;
-  }
+  m->x = add_to_avg(m->x, x, n);
+  m->x2 = add_to_avg(m->x2, pow(x, 2), n);
 
-  if (no_ps) {
-    free(ps);
-  }
+  m->m2 = add_to_avg(m->m2, pow(x - m->x, 2), n);
+  m->m4 = add_to_avg(m->m4, pow(x - m->x, 4), n);
 
-  if (no_r) {
-    gsl_rng_free(r);
+  if (n > 1) {
+    double s2 = n / (n - 1.) * (m->x2 - pow(m->x, 2));
+    m->dx = sqrt(s2 / n);
+    m->c = s2;
+    m->dc = sqrt((m->m4 - (n - 3.)/(n - 1.) * pow(m->m2, 2)) / n);
   }
 
-  return n_flips;
-}
-
-double add_to_avg(double mx, double x, uint64_t n) {
-  return mx * (n / (n + 1.)) + x * 1. / (n + 1.);
+  (m->n)++;
 }
 
 void update_autocorr(autocorr_t *OO, double O) {
@@ -260,7 +97,7 @@ void update_autocorr(autocorr_t *OO, double O) {
 
   dll_t *Otmp = OO->Op;
   dll_t *Osave;
-  uint64_t t = 0;
+  count_t t = 0;
 
   while (Otmp != NULL) {
     OO->OO[t] = add_to_avg(OO->OO[t], O * (Otmp->x), OO->n - t - 1);
@@ -286,25 +123,6 @@ void update_autocorr(autocorr_t *OO, double O) {
   OO->n++;
 }
 
-double rho(autocorr_t *o, uint64_t i) {
+double rho(autocorr_t *o, count_t i) {
   return (o->OO[i] - pow(o->O, 2)) / (o->O2 - pow(o->O, 2));
 }
-
-void update_meas(meas_t *m, double x) {
-  uint64_t n = m->n;
-
-  m->x = add_to_avg(m->x, x, n);
-  m->x2 = add_to_avg(m->x2, pow(x, 2), n);
-
-  m->m2 = add_to_avg(m->m2, pow(x - m->x, 2), n);
-  m->m4 = add_to_avg(m->m4, pow(x - m->x, 4), n);
-
-  if (n > 1) {
-    double s2 = n / (n - 1.) * (m->x2 - pow(m->x, 2));
-    m->dx = sqrt(s2 / n);
-    m->c = s2;
-    m->dc = sqrt((m->m4 - (n - 3.)/(n - 1.) * pow(m->m2, 2)) / n);
-  }
-
-  (m->n)++;
-}
diff --git a/src/wolff.c b/src/wolff.c
index f3f45cf..2a10b7a 100644
--- a/src/wolff.c
+++ b/src/wolff.c
@@ -2,113 +2,53 @@
 #include <wolff.h>
 
 int main(int argc, char *argv[]) {
+
+  L_t L = 128;
+  count_t N = (count_t)1e7;
+  count_t min_runs = 10;
+  count_t n = 3;
+  q_t q = 2;
+  D_t D = 2;
+  double T = 2.26918531421;
+  double *H = (double *)calloc(MAX_Q, sizeof(double));
+  double eps = 0;
+  bool pretend_ising = false;
+
   int opt;
-  sim_t sim;
-  bool output_state, use_dual, M_stop, record_autocorrelation;
-  lattice_t lat;
-  uint16_t L;
-  uint32_t min_runs, lattice_i, sim_i;
-  uint64_t N, n, W, ac_skip;
-  double T, H, eps;
-
-  sim = WOLFF;
-  L = 128;
-  N = 100000000000000;
-  W = 10;
-  ac_skip = 1;
-  n = 3;
-  lat = SQUARE_LATTICE;
-  use_dual = false;
-  M_stop = false;
-  record_autocorrelation = false;
-  T = 2.3;
-  H = 0;
-  eps = 0;
-  output_state = false;
-  min_runs = 10;
-
-  while ((opt = getopt(argc, argv, "n:N:L:T:H:m:S:e:oq:DMas:W:")) != -1) {
+  q_t H_ind = 0;
+
+  while ((opt = getopt(argc, argv, "N:n:D:L:q:T:H:m:e:I")) != -1) {
     switch (opt) {
     case 'N':
-      N = (uint64_t)atof(optarg);
-      break;
-    case 'W':
-      W = (uint64_t)atof(optarg);
-      break;
-    case 'S':
-      ac_skip = (uint64_t)atof(optarg);
+      N = (count_t)atof(optarg);
       break;
     case 'n':
-      n = (uint64_t)atof(optarg);
+      n = (count_t)atof(optarg);
+      break;
+    case 'D':
+      D = atoi(optarg);
       break;
     case 'L':
       L = atoi(optarg);
       break;
+    case 'q':
+      q = atoi(optarg);
+      break;
     case 'T':
       T = atof(optarg);
       break;
     case 'H':
-      H = atof(optarg);
+      H[H_ind] = atof(optarg);
+      H_ind++;
       break;
     case 'm':
       min_runs = atoi(optarg);
       break;
-    case 'M':
-      M_stop = true;
-      break;
     case 'e':
       eps = atof(optarg);
       break;
-    case 'o':
-      output_state = true;
-      break;
-    case 'D':
-      use_dual = true;
-      break;
-    case 'a':
-      record_autocorrelation = true;
-      break;
-    case 's':
-      sim_i = atoi(optarg);
-      switch (sim_i) {
-      case 0:
-        sim = WOLFF;
-        break;
-      case 1:
-        sim = WOLFF_GHOST;
-        break;
-      case 2:
-        sim = METROPOLIS;
-        break;
-      default:
-        printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 "
-               "(DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n");
-        exit(EXIT_FAILURE);
-      }
-      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);
-      }
+    case 'I':
+      pretend_ising = true;
       break;
     default:
       exit(EXIT_FAILURE);
@@ -116,180 +56,133 @@ int main(int argc, char *argv[]) {
   }
 
   gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
-  gsl_rng_set(r, jst_rand_seed());
+  gsl_rng_set(r, rand_seed());
 
-  graph_t *h = graph_create(lat, TORUS_BOUND, L, use_dual);
+  if (pretend_ising) {
+    q = 2;
+    T /= 2;
+    H[0] /= 2;
+    H[1] = -H[0];
+  }
 
   ising_state_t *s = (ising_state_t *)calloc(1, sizeof(ising_state_t));
 
+  graph_t *h = graph_create_square(D, L);
   s->g = graph_add_ext(h);
-  s->spins = (bool *)calloc(h->nv + 1, sizeof(bool));
-  s->M = sign(H) * h->nv;
-  s->H = -(1.0 * h->ne) - sign (H) * H * h->nv;
 
-  double *bond_probs = (double *)malloc(2 * sizeof(double));
-  bond_probs[0] = 1 - exp(-2 / T);
-  bond_probs[1] = 1 - exp(-2 * fabs(H) / T);
+  s->q = q;
+
+  s->spins = (q_t *)calloc(h->nv + 1, sizeof(q_t));
+
+  s->T = T;
+  s->H = H;
+
+  s->T_prob = 1.0 - exp(- 1.0 / T);
+  s->H_probs = (double *)calloc(pow(q, 2), sizeof(double));
+  for (q_t i = 0; i < q; i++) {
+    for (q_t j = 0; j < q; j++) {
+      s->H_probs[q * i + j] = 1.0 - exp((s->H[i] - s->H[j]) / T);
+    }
+  }
+
+  s->M = (v_t *)calloc(q, sizeof(v_t));
+  s->M[0] = h->nv;
+  s->E = - ((double)h->nv) * s->H[0];
 
   double diff = 1e31;
-  uint64_t n_runs = 0;
-  uint64_t n_steps = 0;
-  double clust_per_sweep = 0;
+  count_t n_runs = 0;
 
-  meas_t *M, *aM, *eM, *mM, *E, *eE, *mE, *clust;
+  meas_t *E, *clust, **M;
+
+  M = (meas_t **)malloc(q * sizeof(meas_t *));
+  for (q_t i = 0; i < q; i++) {
+    M[i] = (meas_t *)calloc(1, sizeof(meas_t));
+  }
 
-  M = calloc(1, sizeof(meas_t));
-  aM = calloc(1, sizeof(meas_t));
-  eM = calloc(1, sizeof(meas_t));
-  mM = calloc(1, sizeof(meas_t));
   E = calloc(1, sizeof(meas_t));
-  eE = calloc(1, sizeof(meas_t));
-  mE = calloc(1, sizeof(meas_t));
   clust = calloc(1, sizeof(meas_t));
 
-  autocorr_t *autocorr;
-  if (record_autocorrelation) {
-    autocorr = (autocorr_t *)calloc(1, sizeof(autocorr_t));
-    autocorr->W = 2 * W + 1;
-    autocorr->OO = (double *)calloc(2 * W + 1, sizeof(double));
-  }
-
   printf("\n");
   while (((diff > eps || diff != diff) && n_runs < N) || n_runs < min_runs) {
     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(E->dx / E->x), M->dx / M->x, E->dc / E->c, M->dc / M->c, clust_per_sweep);
+           n_runs, fabs(E->dx / E->x), M[0]->dx / M[0]->x, E->dc / E->c, M[0]->dc / M[0]->c, h->nv / clust->x);
 
-    uint32_t n_flips = 0;
-    uint32_t n_clust = 0;
+    count_t n_flips = 0;
 
     while (n_flips / h->nv < n) {
-      uint32_t tmp_flips = wolff_step(T, H, s, sim, r, bond_probs);
-      n_flips += tmp_flips;
-      n_clust++;
-      if (n_runs > 0){
-        n_steps++;
-      }
+      v_t v0 = gsl_rng_uniform_int(r, h->nv);
+      q_t step = 1 + gsl_rng_uniform_int(r, q - 1);
 
-      if (record_autocorrelation && n_runs > 0) {
-        if (n_steps % ac_skip == 0) {
-          update_autocorr(autocorr, s->H);
-        }
-        update_meas(clust, tmp_flips);
-      }
-    }
+      v_t tmp_flips = flip_cluster(s, v0, step, r);
+      n_flips += tmp_flips;
 
-    double HH = 1;
-    if (H < 0) {
-      HH = -1;
+      update_meas(clust, tmp_flips);
     }
 
-    update_meas(M, (double)(s->M));
-    update_meas(aM, HH * fabs((double)(s->M)));
-    update_meas(E, s->H);
-
-    if (s->M > 0) {
-      update_meas(eM, HH * fabs((double)(s->M)));
-      update_meas(eE, s->H);
-    } else {
-      update_meas(mM, - HH * fabs((double)(s->M)));
-      update_meas(mE, s->H);
+    for (q_t i = 0; i < q; i++) {
+      update_meas(M[i], s->M[i]);
     }
 
-    if (M_stop) {
-      diff = fabs(eM->dx / eM->x);
-    } else {
-      diff = fabs(eM->dc / eM->c);
-    }
+    update_meas(E, s->E);
 
-    clust_per_sweep = add_to_avg(clust_per_sweep, n_clust * 1. / n, n_runs);
+    diff = fabs(M[0]->dc / M[0]->c);
 
     n_runs++;
   }
 
   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(E->dx / E->x), M->dx / M->x, E->dc / E->c, M->dc / M->c, clust_per_sweep);
-
-  FILE *outfile = fopen("out.dat", "a");
-
-  double tau = 0;
-  double dtau = 0;
-
-  if (record_autocorrelation) {
-    double *Gammas = (double *)malloc((W + 1) * sizeof(double));
-
-    Gammas[0] = 1 + rho(autocorr, 0);
-    for (uint64_t i = 0; i < W; i++) {
-      Gammas[1 + i] = rho(autocorr, 2 * i + 1) + rho(autocorr, 2 * i + 2);
-    } 
-
-    uint64_t n;
-    for (n = 0; n < W + 1; n++) {
-      if (Gammas[n] <= 0) {
-        break;
+         n_runs, fabs(E->dx / E->x), M[0]->dx / M[0]->x, E->dc / E->c, M[0]->dc / M[0]->c, h->nv / clust->x);
+
+  FILE *outfile = fopen("out.m", "a");
+
+  if (pretend_ising) {
+    fprintf(outfile,
+            "%u %.15f %" PRIu64 "  %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n",
+            L, 2 * T, n_runs,
+            (2 * E->x - h->ne) / h->nv, 2 * E->dx / h->nv, (M[0]->x - M[1]->x) / h->nv, M[0]->dx / h->nv, 4 * E->c / h->nv, 4 * E->dc / h->nv, 4 * M[0]->c / h->nv, 4 * M[0]->dc / h->nv); 
+  } else {
+    fprintf(outfile, "{\"D\"->%" PRID ",\"L\"->%" PRIL ",\"q\"->%" PRIq ",\"T\"->%.15f,\"H\"->{", D, L, q, T);
+    for (q_t i = 0; i < q; i++) {
+      fprintf(outfile, "%.15f", H[i]);
+      if (i != q-1) {
+        fprintf(outfile, ",");
       }
     }
-
-    if (n == W + 1) {
-      printf("WARNING: correlation function never hit the noise floor.\n");
-    }
-
-    if (n < 2) {
-      printf("WARNING: correlation function only has one nonnegative term.\n");
-    }
-
-    double *conv_Gamma = get_convex_minorant(n, Gammas);
-
-    double ttau = - 0.5;
-
-    for (uint64_t i = 0; i < n + 1; i++) {
-      ttau += conv_Gamma[i];
+    fprintf(outfile, "},\"E\"->{%.15f,%.15f},\"C\"->{%.15f,%.15f},\"M\"->{", E->x / h->nv, E->dx / h->nv, E->c / h->nv, E->dc / h->nv);
+    for (q_t i = 0; i < q; i++) {
+      fprintf(outfile, "{%.15f,%.15f}", M[i]->x / h->nv, M[i]->dx / h->nv);
+      if (i != q-1) {
+        fprintf(outfile, ",");
+      }
     }
-    
-    free(Gammas);
-    free(autocorr->OO);
-    while (autocorr->Op != NULL) {
-      stack_pop_d(&(autocorr->Op));
+    fprintf(outfile, "},\"X\"->{");
+    for (q_t i = 0; i < q; i++) {
+      fprintf(outfile, "{%.15f,%.15f}", M[i]->c / h->nv, M[i]->dc / h->nv);
+      if (i != q-1) {
+        fprintf(outfile, ",");
+      }
     }
-    free(autocorr);
-    
-    tau = ttau * ac_skip * clust->x / h->nv;
-    dtau = 0;
+    fprintf(outfile, "},\"n\"->{%.15f,%.15f}}\n", clust->c / h->nv, clust->dc / h->nv);
   }
 
-  fprintf(outfile,
-          "%u %.15f %.15f %" PRIu64 " %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n",
-          L, T, H, n_runs, n_steps,
-          E->x / h->nv, E->dx / h->nv, M->x / h->nv, M->dx / h->nv, E->c / h->nv, E->dc / h->nv, M->c / h->nv, M->dc / h->nv,
-          eE->n, eE->x / h->nv, eE->dx / h->nv, eM->x / h->nv, eM->dx / h->nv, eE->c / h->nv, eE->dc / h->nv, eM->c / h->nv, eM->dc / h->nv,
-          mE->n, mE->x / h->nv, mE->dx / h->nv, mM->x / h->nv, mM->dx / h->nv, mE->c / h->nv, mE->dc / h->nv, mM->c / h->nv, mM->dc / h->nv,
-          aM->x / h->nv, aM->dx / h->nv, aM->c / h->nv, aM->dc / h->nv, tau, dtau); 
   fclose(outfile);
 
-  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);
+  free(E);
+  free(clust);
+  for (q_t i = 0; i < q; i++) {
+    free(M[i]);
   }
-
-  gsl_rng_free(r);
-  graph_free(s->g);
+  free(M);
+  free(s->H_probs);
+  free(s->M);
   free(s->spins);
+  free(s->g);
   free(s);
-  free(M);
-  free(aM);
-  free(eM);
-  free(mM);
-  free(E);
-  free(eE);
-  free(mE);
-  free(clust);
-  free(bond_probs);
-  graph_free(h);
+  free(H);
+  free(h);
 
   return 0;
 }
+