removed a lot of research code to simplify library and examples for publication

18 files changed, 845 insertions, 223 deletions

diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp
index b730c8d..b78cd4b 100644
--- a/lib/include/wolff.hpp
+++ b/lib/include/wolff.hpp
@@ -1,30 +1,27 @@
 
+#ifndef WOLFF_H
+#define WOLFF_H
+
 #include "wolff/cluster.hpp"
-#include "wolff/state.hpp"
 
 template <class R_t, class X_t>
-void wolff(count_t N, state_t <R_t, X_t>& s, std::function <R_t(std::mt19937&, X_t)> gen_R, wolff_measurement<R_t, X_t>& m, std::mt19937& r) {
-
-#ifdef FINITE_STATES
-#ifdef NOFIELD
-  initialize_probs(s.J, s.T);
-#else
-  initialize_probs(s.J, s.H, s.T);
-#endif
-#endif
+void wolff(N_t N, wolff_system<R_t, X_t>& S,
+           std::function <R_t(std::mt19937&, X_t)> r_gen,
+           wolff_measurement<R_t, X_t>& A, std::mt19937& rng) {
 
-  std::uniform_int_distribution<v_t> dist(0, s.nv);
+  std::uniform_int_distribution<v_t> dist(0, S.nv - 1);
 
-  for (count_t steps = 0; steps < N; steps++) {
-    v_t v0 = dist(r);
-    R_t step = gen_R(r, s.spins[v0]);
+  for (N_t n = 0; n < N; n++) {
+    v_t i0 = dist(rng);
+    R_t r = r_gen(rng, S.s[i0]);
 
-    m.pre_cluster(s, steps, N, v0, step);
+    A.pre_cluster(n, N, S, i0, r);
 
-    flip_cluster<R_t, X_t>(s, v0, step, r, m);
+    wolff_cluster_flip<R_t, X_t>(S, i0, r, rng, A);
 
-    m.post_cluster(s, steps, N);
+    A.post_cluster(n, N, S);
   }
-
 }
 
+#endif
+
diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp
index 805e2c3..055cdf3 100644
--- a/lib/include/wolff/cluster.hpp
+++ b/lib/include/wolff/cluster.hpp
@@ -1,120 +1,122 @@
 
-#pragma once
+#ifndef WOLFF_CLUSTER_H
+#define WOLFF_CLUSTER_H
 
 #include <random>
 #include <cmath>
 #include <vector>
-#include <stack>
+#include <queue>
 
 #include "types.h"
-#include "state.hpp"
+#include "system.hpp"
 #include "graph.hpp"
-#include "meas.h"
+#include "measurement.hpp"
 
 template <class R_t, class X_t>
-void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, std::mt19937& rand, wolff_measurement<R_t, X_t>& m) {
+void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r,
+                  std::mt19937& rng, wolff_measurement<R_t, X_t>& A) {
   std::uniform_real_distribution<double> dist(0.0, 1.0);
 
-  std::stack<v_t> stack;
-  stack.push(v0);
+  std::queue<v_t> queue;
+  queue.push(i0);
 
-  std::vector<bool> marks(s.g.nv, false);
+  std::vector<bool> marks(S.G.nv, false);
 
-  while (!stack.empty()) {
-    v_t v = stack.top();
-    stack.pop();
+  while (!queue.empty()) {
+    v_t i = queue.front();
+    queue.pop();
 
-    if (!marks[v]) { // don't reprocess anyone we've already visited!
-      marks[v] = true;
+    if (!marks[i]) { // don't reprocess anyone we've already visited!
+      marks[i] = true;
 
       X_t si_new;
-#ifndef NOFIELD
-      R_t R_new;
+#ifndef WOLFF_NO_FIELD
+      R_t s0_new;
 
-      bool v_is_ghost = (v == s.nv); // ghost site has the last index
+      bool we_are_ghost = (i == S.nv);
 
-      if (v_is_ghost) {
-        R_new = r.act(s.R); // if we are, then we're moving the transformation
+      if (we_are_ghost) {
+        s0_new = r.act(S.s0);
       } else
 #endif
       {
-        si_new = r.act(s.spins[v]); // otherwise, we're moving the spin at our site
+        si_new = r.act(S.s[i]);
       }
 
-      for (const v_t &vn : s.g.v_adj[v]) {
-        double dE, prob;
+      for (const v_t &j : S.G.adj[i]) {
+        double dE, p;
 
-#ifndef NOFIELD
-        bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost
+#ifndef WOLFF_NO_FIELD
+        bool neighbor_is_ghost = (j == S.nv);
 
-        if (v_is_ghost || vn_is_ghost) { // this is a ghost-involved bond
-          X_t rs_old, rs_new;
+        if (we_are_ghost || neighbor_is_ghost) {
+          X_t s0s_old, s0s_new;
           v_t non_ghost;
 
-          if (vn_is_ghost) {
-            // if our neighbor is the ghost, the current site is a normal
-            // spin - rotate it back!
-            rs_old = s.R.act_inverse(s.spins[v]);
-            rs_new = s.R.act_inverse(si_new);
-            non_ghost = v;
+          if (neighbor_is_ghost) {
+            non_ghost = i;
+            s0s_old = S.s0.act_inverse(S.s[i]);
+            s0s_new = S.s0.act_inverse(si_new);
           } else {
-            /* if we're the ghost, we need to rotate our neighbor back in
-               both the old and new ways */
-            rs_old = s.R.act_inverse(s.spins[vn]);
-            rs_new = R_new.act_inverse(s.spins[vn]);
-            non_ghost = vn;
+            non_ghost = j;
+            s0s_old = S.s0.act_inverse(S.s[j]);
+            s0s_new = s0_new.act_inverse(S.s[j]);
           }
 
-#ifdef SITE_DEPENDENCE
-          dE = s.H(non_ghost, rs_old) - s.H(non_ghost, rs_new);
+#ifdef WOLFF_SITE_DEPENDENCE
+          dE = S.B(non_ghost, s0s_old) - S.B(non_ghost, s0s_new);
 #else
-          dE = s.H(rs_old) - s.H(rs_new);
+          dE = S.B(s0s_old) - S.B(s0s_new);
 #endif
 
-#ifdef FINITE_STATES
-          prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)];
+#ifdef WOLFF_FINITE_STATES
+          p = finite_states_Bp[finite_states_enum(s0s_old)]
+                              [finite_states_enum(s0s_new)];
 #endif
 
           // run measurement hooks for encountering a ghost bond
-          m.ghost_bond_added(non_ghost, rs_old, rs_new, dE);
+          A.ghost_bond_visited(S, non_ghost, s0s_old, s0s_new, dE);
         } else // this is a perfectly normal bond!
 #endif
         {
-#ifdef BOND_DEPENDENCE
-          dE = s.J(v, s.spins[v], vn, s.spins[vn]) - s.J(v, si_new, vn, s.spins[vn]);
+#ifdef WOLFF_BOND_DEPENDENCE
+          dE = S.Z(i, S.s[i], j, S.s[j]) - S.Z(i, si_new, j, S.s[j]);
 #else
-          dE = s.J(s.spins[v], s.spins[vn]) - s.J(si_new, s.spins[vn]);
+          dE = S.Z(S.s[i], S.s[j]) - S.Z(si_new, S.s[j]);
 #endif
 
-
-#ifdef FINITE_STATES
-          prob = J_probs[state_to_ind(s.spins[v])][state_to_ind(si_new)][state_to_ind(s.spins[vn])];
+#ifdef WOLFF_FINITE_STATES
+          p = finite_states_Zp[finite_states_enum(S.s[i])]
+                              [finite_states_enum(si_new)]
+                              [finite_states_enum(S.s[j])];
 #endif
 
           // run measurement hooks for encountering a plain bond
-          m.plain_bond_added(v, s.spins[v], si_new, vn, s.spins[vn], dE);
+          A.plain_bond_visited(S, i, si_new, j, dE);
         }
 
 #ifndef FINITE_STATES
-        prob = 1.0 - exp(-dE / s.T);
+        p = 1.0 - exp(-dE / S.T);
 #endif
 
-        if (dist(rand) < prob) {
-          stack.push(vn); // push the neighboring vertex to the stack
+        if (dist(rng) < p) {
+          queue.push(j); // push the neighboring vertex to the queue 
         }
       }
 
-#ifndef NOFIELD
-      if (v_is_ghost) {
-        m.ghost_site_transformed(s.R, R_new);
-        s.R = R_new;
+#ifndef WOLFF_NO_FIELD
+      if (we_are_ghost) {
+        A.ghost_site_transformed(S, s0_new);
+        S.s0 = s0_new;
       } else
 #endif
       {
-        m.plain_site_transformed(v, s.spins[v], si_new);
-        s.spins[v] = si_new;
+        A.plain_site_transformed(S, i, si_new);
+        S.s[i] = si_new;
       }
     }
   }
 }
 
+#endif
+
diff --git a/lib/include/wolff/finite_states.hpp b/lib/include/wolff/finite_states.hpp
index 426edad..4cca69e 100644
--- a/lib/include/wolff/finite_states.hpp
+++ b/lib/include/wolff/finite_states.hpp
@@ -1,40 +1,35 @@
 
-#pragma once
+
+#ifndef WOLFF_FINITE_STATES
+#define WOLFF_FINITE_STATES
 
 #include <cmath>
-#include <functional>
 #include <array>
 
-#define FINITE_STATES
-
-// must have N_STATES, states[N_STATES], and state_to_ind defined before
-// invoking header
+#include "system.hpp"
 
-std::array<std::array<std::array<double, N_STATES>, N_STATES>, N_STATES> J_probs;
-#ifndef NOFIELD
-std::array<std::array<double, N_STATES>, N_STATES> H_probs;
+std::array<std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Zp;
+#ifndef WOLFF_NO_FIELD
+std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Bp;
 #endif
 
-template <class X_t>
-#ifndef NOFIELD
-void initialize_probs(std::function <double(X_t, X_t)> J, std::function <double(X_t)> H, double T) {
-  for (q_t i = 0; i < N_STATES; i++) {
-    for (q_t j = 0; j < N_STATES; j++) {
-      H_probs[i][j] = 1.0 - exp(-(H(states[i]) - H(states[j])) / T);
+template <class R_t, class X_t>
+void finite_states_init(const wolff_system<R_t, X_t>& S) {
+#ifndef WOLFF_NO_FIELD
+  for (q_t i = 0; i < WOLFF_FINITE_STATES_N; i++) {
+    for (q_t j = 0; j < WOLFF_FINITE_STATES_N; j++) {
+      finite_states_Bp[i][j] = 1.0 - exp(-(S.B(finite_states_possible[i]) - S.B(finite_states_possible[j])) / S.T);
     }
   }
-#else
-void initialize_probs(std::function <double(X_t, X_t)> J, double T) {
 #endif
-  for (q_t i = 0; i < N_STATES; i++) {
-    for (q_t j = 0; j < N_STATES; j++) {
-      for (q_t k = 0; k < N_STATES; k++) {
-        J_probs[i][j][k] = 1.0 - exp(-(J(states[i], states[k]) - J(states[j], states[k])) / T);
+  for (q_t i = 0; i < WOLFF_FINITE_STATES_N; i++) {
+    for (q_t j = 0; j < WOLFF_FINITE_STATES_N; j++) {
+      for (q_t k = 0; k < WOLFF_FINITE_STATES_N; k++) {
+        finite_states_Zp[i][j][k] = 1.0 - exp(-(S.Z(finite_states_possible[i], finite_states_possible[k]) - S.Z(finite_states_possible[j], finite_states_possible[k])) / S.T);
       }
     }
   }
 }
 
-
-
+#endif
 
diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp
index 0aeb6af..fcd73fa 100644
--- a/lib/include/wolff/graph.hpp
+++ b/lib/include/wolff/graph.hpp
@@ -1,5 +1,6 @@
 
-#pragma once
+#ifndef WOLFF_GRAPH_H
+#define WOLFF_GRAPH_H
 
 #include <cmath>
 #include <vector>
@@ -15,10 +16,12 @@ class graph_t {
   public:
     v_t ne;
     v_t nv;
-    std::vector<std::vector<v_t>> v_adj;
+    std::vector<std::vector<v_t>> adj;
     std::vector<std::vector<double>> coordinate;
 
     graph_t(D_t D, L_t L, lattice_t lat = SQUARE_LATTICE);
     void add_ext();
 };
 
+#endif
+
diff --git a/lib/include/wolff/meas.h b/lib/include/wolff/meas.h
deleted file mode 100644
index 4895eac..0000000
--- a/lib/include/wolff/meas.h
+++ /dev/null
@@ -1,19 +0,0 @@
-
-#pragma once
-
-#include "types.h"
-
-template <class R_t, class X_t>
-class wolff_measurement {
-  public:
-    virtual void pre_cluster(const state_t<R_t, X_t>&, count_t, count_t, v_t, const R_t&) = 0;
-
-    virtual void plain_bond_added(v_t, const X_t&, const X_t&, v_t, const X_t&, double) = 0;
-    virtual void ghost_bond_added(v_t, const X_t&, const X_t&, double) = 0;
-
-    virtual void plain_site_transformed(v_t, const X_t&, const X_t&) = 0;
-    virtual void ghost_site_transformed(const R_t&, const R_t&) = 0;
-
-    virtual void post_cluster(const state_t<R_t, X_t>&, count_t, count_t) = 0;
-};
-
diff --git a/lib/include/wolff/measurement.hpp b/lib/include/wolff/measurement.hpp
new file mode 100644
index 0000000..a6a5f79
--- /dev/null
+++ b/lib/include/wolff/measurement.hpp
@@ -0,0 +1,24 @@
+
+#ifndef WOLFF_MEASUREMENTS
+#define WOLFF_MEASUREMENTS
+
+#include "system.hpp"
+
+template <class R_t, class X_t>
+class wolff_measurement {
+  public:
+    virtual void pre_cluster(N_t, N_t, const wolff_system<R_t, X_t>&, v_t, const R_t&) = 0;
+
+    virtual void plain_bond_visited(const wolff_system<R_t, X_t>&, v_t, const X_t&, v_t, double) = 0;
+    virtual void plain_site_transformed(const wolff_system<R_t, X_t>&, v_t, const X_t&) = 0;
+
+#ifndef WOLFF_NO_FIELD
+    virtual void ghost_bond_visited(const wolff_system<R_t, X_t>&, v_t, const X_t&, const X_t&, double) = 0;
+    virtual void ghost_site_transformed(const wolff_system<R_t, X_t>&, const R_t&) = 0;
+#endif
+
+    virtual void post_cluster(N_t, N_t, const wolff_system<R_t, X_t>&) = 0;
+};
+
+#endif
+
diff --git a/lib/include/wolff/models/dihedral.hpp b/lib/include/wolff/models/dihedral.hpp
new file mode 100644
index 0000000..cbc5687
--- /dev/null
+++ b/lib/include/wolff/models/dihedral.hpp
@@ -0,0 +1,48 @@
+
+#pragma once
+
+#include <wolff/types.h>
+#include "potts.hpp"
+
+template <q_t q>
+class dihedral_t {
+  public:
+    bool is_reflection;
+    q_t x;
+
+    dihedral_t() : is_reflection(false), x(0) {}
+    dihedral_t(bool x, q_t y) : is_reflection(x), x(y) {}
+
+    potts_t<q> act(const potts_t<q>& s) const {
+      if (this->is_reflection) {
+        return potts_t<q>(((q + this->x) - s.x) % q);
+      } else {
+        return potts_t<q>((this->x + s.x) % q);
+      }
+    }
+
+    dihedral_t<q> act(dihedral_t<q> r) const {
+      if (this->is_reflection) {
+        return dihedral_t<q>(!(r.is_reflection), ((q + this->x) - r.x) % q); 
+      } else {
+        return dihedral_t<q>(r.is_reflection, (this->x + r.x) % q);
+      }
+    }
+
+    potts_t<q> act_inverse(potts_t<q> s) const {
+      if (this->is_reflection) {
+        return this->act(s);
+      } else {
+        return potts_t<q>(((s.x + q) - this->x) % q);
+      }
+    }
+
+    dihedral_t<q> act_inverse(dihedral_t<q> r) const {
+      if (this->is_reflection) {
+        return this->act(r);
+      } else {
+        return dihedral_t<q>(r.is_reflection, ((r.x + q) - this->x) % q);
+      }
+    }
+};
+
diff --git a/lib/include/wolff/models/dihedral_inf.hpp b/lib/include/wolff/models/dihedral_inf.hpp
new file mode 100644
index 0000000..bded387
--- /dev/null
+++ b/lib/include/wolff/models/dihedral_inf.hpp
@@ -0,0 +1,47 @@
+
+#include "../types.h"
+#include <cmath>
+#include "height.hpp"
+
+template <class T>
+class dihedral_inf_t {
+  public:
+  bool is_reflection;
+  T x;
+
+  dihedral_inf_t() : is_reflection(false), x(0) {}
+  dihedral_inf_t(bool x, T y) : is_reflection(x), x(y) {}
+
+  height_t<T> act(const height_t<T>& h) const {
+    if (this->is_reflection) {
+      return height_t(this->x - h.x);
+    } else {
+      return height_t(this->x + h.x);
+    }
+  }
+
+  dihedral_inf_t<T> act(const dihedral_inf_t<T>& r) const {
+    if (this->is_reflection) {
+      return dihedral_inf_t<T>(!r.is_reflection, this->x - r.x);
+    } else {
+      return dihedral_inf_t<T>(r.is_reflection, this->x + r.x);
+    }
+  }
+
+  height_t<T> act_inverse(const height_t<T>& h) const {
+    if (this->is_reflection) {
+      return this->act(h);
+    } else {
+      return height_t(h.x - this->x);
+    }
+  }
+
+  dihedral_inf_t<T> act_inverse(const dihedral_inf_t<T>& r) const {
+    if (this->is_reflection) {
+      return this->act(r);
+    } else {
+      return dihedral_inf_t<T>(r.is_reflection, r.x - this->x);
+    }
+  }
+};
+
diff --git a/lib/include/wolff/models/height.hpp b/lib/include/wolff/models/height.hpp
new file mode 100644
index 0000000..bd0ceb6
--- /dev/null
+++ b/lib/include/wolff/models/height.hpp
@@ -0,0 +1,43 @@
+
+#pragma once
+
+#include <cmath>
+#include <wolff/types.h>
+
+template <class T>
+struct height_t {
+  T x;
+
+  height_t() : x(0) {}
+  height_t(T x) : x(x) {}
+
+  typedef T M_t;
+  typedef double F_t;
+
+  inline T operator*(v_t a) const {
+    return x * a;
+  }
+
+  inline double operator*(double a) const {
+    return x * a;
+  }
+
+  inline T operator-(const height_t& h) const {
+    return x - h.x;
+  }
+};
+
+template <class T>
+inline T& operator+=(T& M, const height_t<T> &h) {
+  M += h.x;
+
+  return M;
+}
+
+template <class T>
+inline T& operator-=(T& M, const height_t<T> &h) {
+  M -= h.x;
+
+  return M;
+}
+
diff --git a/lib/include/wolff/models/ising.hpp b/lib/include/wolff/models/ising.hpp
new file mode 100644
index 0000000..c6a3a47
--- /dev/null
+++ b/lib/include/wolff/models/ising.hpp
@@ -0,0 +1,81 @@
+
+#ifndef WOLFF_MODELS_ISING
+#define WOLFF_MODELS_ISING
+
+#include "../types.h"
+
+class ising_t {
+  public:
+    bool x;
+
+    ising_t() : x(false) {}
+
+    ising_t(bool x) : x(x) {} 
+    ising_t(int x) : x((bool)x) {}
+
+    ising_t act(const ising_t& s) const {
+      if (x) {
+        return ising_t(!s.x);
+      } else {
+        return ising_t(s.x);
+      }
+    }
+
+    ising_t act_inverse(const ising_t& s) const {
+      return this->act(s);
+    }
+
+    typedef int M_t;
+    typedef double F_t;
+
+    inline M_t operator*(v_t a) const {
+      if (x) {
+        return -(M_t)a;
+      } else {
+        return (M_t)a;
+      }
+    }
+
+    inline F_t operator*(double a) const {
+      if (x) {
+        return -a;
+      } else {
+        return a;
+      }
+    }
+
+    inline M_t operator-(const ising_t &s) const {
+      if (x == s.x) {
+        return 0;
+      } else {
+        if (x) {
+          return -2;
+        } else {
+          return 2;
+        }
+      }
+    }
+
+    inline int operator*(const ising_t& s) const {
+      if (x == s.x) {
+        return 1;
+      } else {
+        return -1;
+      }
+    }
+};
+
+inline ising_t::M_t operator*(v_t a, const ising_t& s) {
+  return s * a;
+}
+
+inline ising_t::F_t operator*(double a, const ising_t& s) {
+  return s * a;
+}
+
+#define WOLFF_FINITE_STATES_N 2
+const ising_t finite_states_possible[2] = {ising_t(0), ising_t(1)};
+q_t finite_states_enum(ising_t state) { return (q_t)state.x; }
+
+#endif
+
diff --git a/lib/include/wolff/models/orthogonal.hpp b/lib/include/wolff/models/orthogonal.hpp
new file mode 100644
index 0000000..9dd5ddd
--- /dev/null
+++ b/lib/include/wolff/models/orthogonal.hpp
@@ -0,0 +1,199 @@
+
+#pragma once
+
+#include <random>
+#include <cmath>
+
+#include <wolff/types.h>
+#include "vector.hpp"
+
+template <q_t q, class T>
+class orthogonal_t : public std::array<std::array<T, q>, q> {
+  public :
+  bool is_reflection;
+
+  orthogonal_t() : is_reflection(false) {
+    for (q_t i = 0; i < q; i++) {
+      (*this)[i].fill(0);
+      (*this)[i][i] = (T)1;
+    }
+  }
+
+  vector_t<q, T> act(const vector_t <q, T>& v) const {
+    vector_t <q, T> v_rot;
+    v_rot.fill(0);
+
+    if (is_reflection) {
+      double prod = 0;
+      for (q_t i = 0; i < q; i++) {
+        prod += v[i] * (*this)[0][i];
+      }
+      for (q_t i = 0; i < q; i++) {
+        v_rot[i] = v[i] - 2 * prod * (*this)[0][i];
+      }
+    } else {
+      for (q_t i = 0; i < q; i++) {
+        for (q_t j = 0; j < q; j++) {
+          v_rot[i] += (*this)[i][j] * v[j];
+        }
+      }
+    }
+
+    return v_rot;
+  }
+
+  orthogonal_t<q, T> act(const orthogonal_t <q, T>& m) const {
+    orthogonal_t <q, T> m_rot;
+
+    m_rot.is_reflection = false;
+
+    if (is_reflection) {
+      for (q_t i = 0; i < q; i++) {
+        double akOki = 0;
+
+        for (q_t k = 0; k < q; k++) {
+          akOki += (*this)[0][k] * m[k][i];
+        }
+
+        for (q_t j = 0; j < q; j++) {
+          m_rot[j][i] = m[j][i] - 2 * akOki * (*this)[0][j];
+        }
+      }
+    } else {
+      for (q_t i = 0; i < q; i++) {
+        m_rot[i].fill(0);
+        for (q_t j = 0; j < q; j++) {
+          for (q_t k = 0; k < q; k++) {
+            m_rot[i][j] += (*this)[i][j] * m[j][k];
+          }
+        }
+      }
+    }
+
+    return m_rot;
+  }
+
+  vector_t <q, T> act_inverse(const vector_t <q, T>& v) const {
+    if (is_reflection) {
+      return this->act(v); // reflections are their own inverse
+    } else {
+      vector_t <q, T> v_rot;
+      v_rot.fill(0);
+
+      for (q_t i = 0; i < q; i++) {
+        for (q_t j = 0; j < q; j++) {
+          v_rot[i] += (*this)[j][i] * v[j];
+        }
+      }
+
+      return v_rot;
+    }
+  }
+
+  vector_t <q, T> act_inverse(const orthogonal_t <q, T>& m) const {
+    if (is_reflection) {
+      return this->act(m); // reflections are their own inverse
+    } else {
+      orthogonal_t <q, T> m_rot;
+      m_rot.is_reflection = false;
+
+      for (q_t i = 0; i < q; i++) {
+        m_rot[i].fill(0);
+        for (q_t j = 0; j < q; j++) {
+          for (q_t k = 0; k < q; k++) {
+            m_rot[i][j] += (*this)[j][i] * m[j][k];
+          }
+        }
+      }
+
+      return m_rot;
+    }
+  }
+
+};
+
+template <q_t q>
+orthogonal_t <q, double> generate_rotation_uniform (std::mt19937& r, const vector_t <q, double>& v) {
+  std::normal_distribution<double> dist(0.0,1.0);
+  orthogonal_t <q, double> ptr;
+  ptr.is_reflection = true;
+
+  double v2 = 0;
+
+  for (q_t i = 0; i < q; i++) {
+    ptr[0][i] = dist(r);
+    v2 += ptr[0][i] * ptr[0][i];
+  }
+
+  double mag_v = sqrt(v2);
+
+  for (q_t i = 0; i < q; i++) {
+    ptr[0][i] /= mag_v;
+  }
+
+  return ptr;
+}
+
+template <q_t q>
+orthogonal_t <q, double> generate_rotation_perturbation (std::mt19937& r, const vector_t <q, double>& v0, double epsilon, unsigned int n) {
+  std::normal_distribution<double> dist(0.0,1.0);
+  orthogonal_t <q, double> m;
+  m.is_reflection = true;
+
+  vector_t <q, double> v;
+
+  if (n > 1) {
+    std::uniform_int_distribution<unsigned int> udist(0, n);
+    unsigned int rotation = udist(r);
+
+    double cosr = cos(2 * M_PI * rotation / (double)n / 2.0);
+    double sinr = sin(2 * M_PI * rotation / (double)n / 2.0);
+
+    v[0] = v0[0] * cosr - v0[1] * sinr;
+    v[1] = v0[1] * cosr + v0[0] * sinr;
+
+    for (q_t i = 2; i < q; i++) {
+      v[i] = v0[i];
+    }
+  } else {
+    v = v0;
+  }
+
+  double m_dot_v = 0;
+
+  for (q_t i = 0; i < q; i++) {
+    m[0][i] = dist(r); // create a random vector
+    m_dot_v += m[0][i] * v[i];
+  }
+
+  double v2 = 0;
+
+  for (q_t i = 0; i < q; i++) {
+    m[0][i] = m[0][i] - m_dot_v * v[i]; // find the component orthogonal to v
+    v2 += pow(m[0][i], 2);
+  }
+
+  double mag_v = sqrt(v2);
+
+  for (q_t i = 0; i < q; i++) {
+    m[0][i] /= mag_v; // normalize
+  }
+
+  v2 = 0;
+
+  double factor = epsilon * dist(r);
+
+  for (q_t i = 0; i < q; i++) {
+    m[0][i] += factor * v[i]; // perturb orthogonal vector in original direction
+    v2 += pow(m[0][i], 2);
+  }
+
+  mag_v = sqrt(v2);
+
+  for (q_t i = 0; i < q; i++) {
+    m[0][i] /= mag_v; // normalize
+  }
+
+  return m;
+}
+
diff --git a/lib/include/wolff/models/potts.hpp b/lib/include/wolff/models/potts.hpp
new file mode 100644
index 0000000..903f25f
--- /dev/null
+++ b/lib/include/wolff/models/potts.hpp
@@ -0,0 +1,50 @@
+
+#pragma once
+
+#include <cmath>
+
+#include "../types.h"
+#include "vector.hpp"
+
+template <q_t q>
+class potts_t {
+  public:
+    q_t x;
+
+    potts_t() : x(0) {}
+    potts_t(q_t x) : x(x) {}
+
+    typedef vector_t<q, int> M_t;
+    typedef vector_t<q, double> F_t;
+
+    inline vector_t<q, int> operator*(v_t a) const {
+      vector_t<q, int> result;
+      result.fill(0);
+      result[x] = (int)a;
+
+      return result;
+    }
+
+    inline vector_t<q, double> operator*(double a) const {
+      vector_t<q, double> result;
+      result.fill(0.0);
+      result[x] = a;
+
+      return result;
+    }
+
+    inline vector_t<q, int> operator-(const potts_t<q> &s) const {
+      vector_t<q, int> result;
+      result.fill(0);
+
+      result[x]++;
+      result[s.x]--;
+
+      return result;
+    }
+};
+
+const potts_t<POTTSQ> states[256] = {{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, {11}, {12}, {13}, {14}, {15}, {16}, {17}, {18}, {19}, {20}, {21}, {22}, {23}, {24}, {25}, {26}, {27}, {28}, {29}, {30}, {31}, {32}, {33}, {34}, {35}, {36}, {37}, {38}, {39}, {40}, {41}, {42}, {43}, {44}, {45}, {46}, {47}, {48}, {49}, {50}, {51}, {52}, {53}, {54}, {55}, {56}, {57}, {58}, {59}, {60}, {61}, {62}, {63}, {64}, {65}, {66}, {67}, {68}, {69}, {70}, {71}, {72}, {73}, {74}, {75}, {76}, {77}, {78}, {79}, {80}, {81}, {82}, {83}, {84}, {85}, {86}, {87}, {88}, {89}, {90}, {91}, {92}, {93}, {94}, {95}, {96}, {97}, {98}, {99}, {100}, {101}, {102}, {103}, {104}, {105}, {106}, {107}, {108}, {109}, {110}, {111}, {112}, {113}, {114}, {115}, {116}, {117}, {118}, {119}, {120}, {121}, {122}, {123}, {124}, {125}, {126}, {127}, {128}, {129}, {130}, {131}, {132}, {133}, {134}, {135}, {136}, {137}, {138}, {139}, {140}, {141}, {142}, {143}, {144}, {145}, {146}, {147}, {148}, {149}, {150}, {151}, {152}, {153}, {154}, {155}, {156}, {157}, {158}, {159}, {160}, {161}, {162}, {163}, {164}, {165}, {166}, {167}, {168}, {169}, {170}, {171}, {172}, {173}, {174}, {175}, {176}, {177}, {178}, {179}, {180}, {181}, {182}, {183}, {184}, {185}, {186}, {187}, {188}, {189}, {190}, {191}, {192}, {193}, {194}, {195}, {196}, {197}, {198}, {199}, {200}, {201}, {202}, {203}, {204}, {205}, {206}, {207}, {208}, {209}, {210}, {211}, {212}, {213}, {214}, {215}, {216}, {217}, {218}, {219}, {220}, {221}, {222}, {223}, {224}, {225}, {226}, {227}, {228}, {229}, {230}, {231}, {232}, {233}, {234}, {235}, {236}, {237}, {238}, {239}, {240}, {241}, {242}, {243}, {244}, {245}, {246}, {247}, {248}, {249}, {250}, {251}, {252}, {253}, {254}, {255}};
+template <q_t q>
+q_t finite_states_enum(potts_t<q> state) { return (q_t)state.x; }
+
diff --git a/lib/include/wolff/models/symmetric.hpp b/lib/include/wolff/models/symmetric.hpp
new file mode 100644
index 0000000..7bcc9a6
--- /dev/null
+++ b/lib/include/wolff/models/symmetric.hpp
@@ -0,0 +1,51 @@
+
+#pragma once
+
+#include <array>
+
+#include "../types.h"
+#include "potts.hpp"
+
+template <q_t q>
+class symmetric_t : public std::array<q_t, q> {
+  public:
+
+    symmetric_t() {
+      for (q_t i = 0; i < q; i++) {
+        (*this)[i] = i;
+      }
+    }
+
+    potts_t<q> act(const potts_t<q> &s) const {
+      return potts_t<q>((*this)[s.x]);
+    }
+
+    symmetric_t<q> act(const symmetric_t<q>& r) const {
+      symmetric_t<q> r_rot;
+      for (q_t i = 0; i < q; i++) {
+        r_rot[i] = (*this)[r[i]];
+      }
+
+      return r_rot;
+    }
+
+    potts_t<q> act_inverse(const potts_t<q>& s) const {
+      for (q_t i = 0; i < q; i++) {
+        if ((*this)[i] == s.x) {
+          return potts_t<q>(i);
+        }
+      }
+
+      exit(EXIT_FAILURE);
+    }
+
+    symmetric_t<q> act_inverse(const symmetric_t<q>& r) const {
+      symmetric_t<q> r_rot;
+      for (q_t i = 0; i < q; i++) {
+        r_rot[(*this)[i]] = r[i];
+      }
+
+      return r_rot;
+    }
+};
+
diff --git a/lib/include/wolff/models/vector.hpp b/lib/include/wolff/models/vector.hpp
new file mode 100644
index 0000000..1d635c8
--- /dev/null
+++ b/lib/include/wolff/models/vector.hpp
@@ -0,0 +1,108 @@
+
+#pragma once
+
+#include <cmath>
+#include <array>
+#include <iostream>
+
+#include <wolff/types.h>
+
+template <q_t q, class T>
+class vector_t : public std::array<T, q> {
+  public: 
+
+    vector_t() {
+      this->fill((T)0);
+      (*this)[0] = (T)1;
+    }
+
+    vector_t(const T *x) {
+      for (q_t i = 0; i < q; i++) {
+        (*this)[i] = x[i];
+      }
+    }
+
+    typedef vector_t <q, T> M_t;
+    typedef vector_t <q, double> F_t;
+
+    template <class U>
+    inline vector_t<q, T>& operator+=(const vector_t<q, U> &v) {
+      for (q_t i = 0; i < q; i++) {
+        (*this)[i] += (T)v[i];
+      }
+      return *this;
+    }
+
+    template <class U>
+    inline vector_t<q, T>& operator-=(const vector_t<q, U> &v) {
+      for (q_t i = 0; i < q; i++) {
+        (*this)[i] -= (T)v[i];
+      }
+      return *this;
+    }
+
+    inline vector_t<q, T> operator*(v_t x) const {
+      vector_t<q, T> result;
+      for (q_t i = 0; i < q; i++) {
+        result[i] = x * (*this)[i];
+      }
+
+      return result;
+    }
+
+    inline vector_t<q, double> operator*(double x) const {
+      vector_t<q, double> result;
+      for (q_t i = 0; i < q; i++) {
+        result[i] = x * (*this)[i];
+      }
+
+      return result;
+    }
+
+    inline vector_t<q, T> operator-(const vector_t<q, T>& v) const {
+      vector_t<q, T> diff = *this;
+      diff -= v;
+      return diff;
+    }
+
+    inline T operator*(const vector_t<q, T>& v) const {
+      double prod = 0;
+
+      for (q_t i = 0; i < q; i++) {
+        prod += v[i] * (*this)[i];
+      }
+
+      return prod;
+    }
+
+    template <class U>
+    inline vector_t<q, T> operator/(U a) const {
+      vector_t<q, T> result;
+      for (q_t i = 0; i < q; i++) {
+        result[i] = (*this)[i] / a;
+      }
+
+      return result;
+    }
+};
+
+template<q_t q, class T>
+inline vector_t<q, T> operator*(v_t a, const vector_t<q, T>&v) {
+  return v * a;
+}
+
+template<q_t q, class T>
+inline vector_t<q, double> operator*(double a, const vector_t<q, T>&v) {
+  return v * a;
+}
+
+template<q_t q, class T>
+std::ostream& operator<<(std::ostream& os, const vector_t<q, T>&v) {
+  os << "( ";
+  for (T vi : v) {
+    os << vi << " ";
+  }
+  os << ")";
+  return os;
+}
+
diff --git a/lib/include/wolff/state.hpp b/lib/include/wolff/state.hpp
deleted file mode 100644
index 4909881..0000000
--- a/lib/include/wolff/state.hpp
+++ /dev/null
@@ -1,68 +0,0 @@
-
-#pragma once
-
-#include <functional>
-#include <vector>
-
-#include "types.h"
-#include "graph.hpp"
-
-template <class R_t, class X_t>
-class state_t {
-  public:
-    D_t D; // the dimension of the system
-    L_t L; // the linear size of the lattice
-    v_t nv; // the number of vertices in the original lattice
-    v_t ne; // the number of edges in the original lattice
-    graph_t g; // the graph defining the lattice with ghost
-    double T; // the temperature
-    std::vector<X_t> spins; // the state of the ordinary spins
-#ifndef NOFIELD
-    R_t R; // the current state of the ghost site
-#endif
-
-#ifdef BOND_DEPENDENCE
-    std::function <double(v_t, const X_t&, v_t, const X_t&)> J; // coupling between sites
-#else
-    std::function <double(const X_t&, const X_t&)> J; // coupling between sites
-#endif
-
-#ifndef NOFIELD
-#ifdef SITE_DEPENDENCE
-    std::function <double(v_t, const X_t&)> H; // coupling with the external field
-#else
-    std::function <double(const X_t&)> H; // coupling with the external field
-#endif
-#endif
-
-    state_t(D_t D, L_t L, double T,
-#ifdef BOND_DEPENDENCE
-        std::function <double(v_t, const X_t&, v_t, const X_t&)> J
-#else
-        std::function <double(const X_t&, const X_t&)> J
-#endif
-#ifndef NOFIELD
-#ifdef SITE_DEPENDENCE
-        , std::function <double(v_t, const X_t&)> H
-#else
-        , std::function <double(const X_t&)> H
-#endif
-#endif
-        , lattice_t lat = SQUARE_LATTICE) : D(D), L(L), g(D, L, lat), T(T),
-#ifndef NOFIELD
-              R(),
-#endif
-              J(J)
-#ifndef NOFIELD
-               , H(H)
-#endif
-    {
-      nv = g.nv;
-      ne = g.ne;
-      spins.resize(nv);
-#ifndef NOFIELD
-      g.add_ext();
-#endif
-    }
-};
-
diff --git a/lib/include/wolff/system.hpp b/lib/include/wolff/system.hpp
new file mode 100644
index 0000000..ac82f44
--- /dev/null
+++ b/lib/include/wolff/system.hpp
@@ -0,0 +1,70 @@
+
+#ifndef WOLFF_STATE_H
+#define WOLFF_STATE_H
+
+#include <functional>
+#include <vector>
+
+#include "types.h"
+#include "graph.hpp"
+
+template <class R_t, class X_t>
+class wolff_system {
+  public:
+    D_t D; // dimension
+    L_t L; // linear size
+    v_t nv; // number of vertices
+    v_t ne; // number of edges
+    graph_t G; // the graph defining the lattice with ghost
+    double T; // the temperature
+    std::vector<X_t> s; // the state of the ordinary spins
+#ifndef WOLFF_NO_FIELD
+    R_t s0; // the current state of the ghost site
+#endif
+
+#ifdef WOLFF_BOND_DEPENDENCE
+    std::function <double(v_t, const X_t&, v_t, const X_t&)> Z; // coupling between sites
+#else
+    std::function <double(const X_t&, const X_t&)> Z; // coupling between sites
+#endif
+
+#ifndef WOLFF_NO_FIELD
+#ifdef WOLFF_SITE_DEPENDENCE
+    std::function <double(v_t, const X_t&)> B; // coupling with the external field
+#else
+    std::function <double(const X_t&)> B; // coupling with the external field
+#endif
+#endif
+
+    wolff_system(D_t D, L_t L, double T,
+#ifdef WOLFF_BOND_DEPENDENCE
+        std::function <double(v_t, const X_t&, v_t, const X_t&)> Z
+#else
+        std::function <double(const X_t&, const X_t&)> Z
+#endif
+#ifndef WOLFF_NO_FIELD
+#ifdef WOLFF_SITE_DEPENDENCE
+        , std::function <double(v_t, const X_t&)> B
+#else
+        , std::function <double(const X_t&)> B
+#endif
+#endif
+        , lattice_t lat = SQUARE_LATTICE) : D(D), L(L), G(D, L, lat), T(T), Z(Z)
+#ifndef WOLFF_NO_FIELD
+             , s0(), B(B)
+#endif
+    {
+      nv = G.nv;
+      ne = G.ne;
+      s.resize(nv);
+#ifndef WOLFF_NO_FIELD
+      G.add_ext();
+#endif
+#ifdef WOLFF_FINITE_STATES
+      finite_states_init(*this);
+#endif
+    }
+};
+
+#endif
+
diff --git a/lib/include/wolff/types.h b/lib/include/wolff/types.h
index ec9efaf..68bd09e 100644
--- a/lib/include/wolff/types.h
+++ b/lib/include/wolff/types.h
@@ -1,36 +1,27 @@
 
 #include <inttypes.h>
 
-typedef uint_fast32_t v_t;
-typedef uint_fast8_t q_t;
-typedef uint_fast16_t R_t;
-typedef uint_fast8_t D_t;
-typedef uint_fast16_t L_t;
-typedef uint_fast64_t count_t;
-typedef int_fast64_t h_t;
+typedef uint_fast32_t v_t;     // vertex and edge indices
+typedef uint_fast8_t q_t;      // enumerated states
+typedef uint_fast8_t D_t;      // dimension
+typedef uint_fast16_t L_t;     // linear size
+typedef uint_fast64_t N_t;     // cycle iterator
 
 #define MAX_v UINT_FAST32_MAX
 #define MAX_Q UINT_FAST8_MAX
-#define MAX_R UINT_FAST16_MAX
 #define MAX_D UINT_FAST8_MAX
 #define MAX_L UINT_FAST16_MAX
-#define MAX_COUNT UINT_FAST64_MAX
-#define MAX_H INT_FAST64_MAX
-#define MIN_H INT_FAST64_MIN
+#define MAX_N UINT_FAST64_MAX
 
 #define PRIv PRIuFAST32
 #define PRIq PRIuFAST8
-#define PRIR PRIuFAST16
 #define PRID PRIuFAST8
 #define PRIL PRIuFAST16
-#define PRIcount PRIuFAST64
-#define PRIh PRIdFAST64
+#define PRIN PRIuFAST64
 
 #define SCNv SCNuFAST32
 #define SCNq SCNuFAST8
-#define SCNR SCNuFAST16
 #define SCND SCNuFAST8
 #define SCNL SCNuFAST16
-#define SCNcount SCNuFAST64
-#define SCNh SCNdFAST64
+#define SCNN SCNuFAST64
 
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index c6f0ba6..b42aa78 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -7,11 +7,11 @@ graph_t::graph_t(D_t D, L_t L, lattice_t lat) {
         nv = pow(L, D);
         ne = D * nv;
 
-        v_adj.resize(nv);
+        adj.resize(nv);
         coordinate.resize(nv);
 
-        for (std::vector<v_t> v_adj_i : v_adj) {
-          v_adj_i.reserve(2 * D);
+        for (std::vector<v_t> adj_i : adj) {
+          adj_i.reserve(2 * D);
         }
 
         for (v_t i = 0; i < nv; i++) {
@@ -19,21 +19,21 @@ graph_t::graph_t(D_t D, L_t L, lattice_t lat) {
           for (D_t j = 0; j < D; j++) {
             coordinate[i][j] = (i / (v_t)pow(L, D - j - 1)) % L;
 
-            v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1)));
-            v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1)));
+            adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1)));
+            adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1)));
           }
         }
         break;
       }
     case DIAGONAL_LATTICE: {
         nv = D * pow(L, D);
-        ne = D * nv;
+        ne = 2 * nv;
 
-        v_adj.resize(nv);
+        adj.resize(nv);
         coordinate.resize(nv);
 
-        for (std::vector<v_t> v_adj_i : v_adj) {
-          v_adj_i.reserve(4 * (D - 1));
+        for (std::vector<v_t> adj_i : adj) {
+          adj_i.reserve(4 * (D - 1));
         }
 
         for (D_t i = 0; i < D; i++) {
@@ -42,10 +42,10 @@ graph_t::graph_t(D_t D, L_t L, lattice_t lat) {
           for (v_t j = 0; j < pow(L, D); j++) {
             v_t vc = sb + j;
 
-            v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + j);
-            v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * (j / (v_t)pow(L, D - 1)) + (j + 1 - 2 * (i % 2)) % L);
-            v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j - i) % L);
-            v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j + 1 - i) % L);
+            adj[vc].push_back(((i + 1) % D) * pow(L, D) + j);
+            adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * (j / (v_t)pow(L, D - 1)) + (j + 1 - 2 * (i % 2)) % L);
+            adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j - i) % L);
+            adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j + 1 - i) % L);
           }
         }
         break;
@@ -54,16 +54,16 @@ graph_t::graph_t(D_t D, L_t L, lattice_t lat) {
 }
 
 void graph_t::add_ext() {
-  for (std::vector<v_t>& v_adj_i : v_adj) {
-    v_adj_i.push_back(nv);
+  for (std::vector<v_t>& adj_i : adj) {
+    adj_i.push_back(nv);
   }
 
-  v_adj.resize(nv + 1);
+  adj.resize(nv + 1);
   coordinate.resize(nv + 1);
-  v_adj[nv].reserve(nv);
+  adj[nv].reserve(nv);
 
   for (v_t i = 0; i < nv; i++) {
-    v_adj[nv].push_back(i);
+    adj[nv].push_back(i);
   }
 
   coordinate[nv].resize(coordinate[0].size());