Large refactoring around changes in the graph class.

  - Graphs now use lists of references instead of vectors of indicies.
  - Vertices and edges have associated classes that can be given arbitrary properties via template specification.
  - All essential library headers have been combined into one, wolff.hpp.

25 files changed, 1025 insertions, 1125 deletions

diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
index 1df9807..445e68f 100644
--- a/lib/CMakeLists.txt
+++ b/lib/CMakeLists.txt
@@ -1,14 +1,10 @@
 
-project(libwolff LANGUAGES C CXX)
+project(libwolff LANGUAGES CXX)
 
-add_library(wolff SHARED src/graph.cpp)
+add_library(wolff INTERFACE)
+target_include_directories(wolff INTERFACE .)
+target_compile_features(wolff INTERFACE cxx_std_17)
 
-target_include_directories(wolff PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
-    $<INSTALL_INTERFACE:include>)
-
-install(TARGETS wolff EXPORT wolffConfig
-    ARCHIVE  DESTINATION ${CMAKE_INSTALL_LIBDIR}
-    LIBRARY  DESTINATION ${CMAKE_INSTALL_LIBDIR}
-    RUNTIME  DESTINATION ${CMAKE_INSTALL_BINDIR})
-install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
+install(FILES "wolff.hpp" DESTINATION include)
+install(DIRECTORY "wolff_models" DESTINATION include)
 
diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp
deleted file mode 100644
index 7e909c8..0000000
--- a/lib/include/wolff.hpp
+++ /dev/null
@@ -1,32 +0,0 @@
-
-#ifndef WOLFF_H
-#define WOLFF_H
-
-#include "wolff/cluster.hpp"
-#include "wolff/measurement.hpp"
-
-namespace wolff{
-
-template <class R_t, class X_t>
-void system<R_t, X_t>::run_wolff(N_t N,
-           std::function <R_t(std::mt19937&, const system<R_t, X_t>&, v_t)> r_gen,
-           measurement<R_t, X_t>& A, std::mt19937& rng) {
-
-  std::uniform_int_distribution<v_t> dist(0, nv - 1);
-
-  for (N_t n = 0; n < N; n++) {
-    v_t i0 = dist(rng);
-    R_t r = r_gen(rng, *this, i0);
-
-    A.pre_cluster(n, N, *this, i0, r);
-
-    this->flip_cluster(i0, r, rng, A);
-
-    A.post_cluster(n, N, *this);
-  }
-}
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp
deleted file mode 100644
index 3912446..0000000
--- a/lib/include/wolff/cluster.hpp
+++ /dev/null
@@ -1,120 +0,0 @@
-
-#ifndef WOLFF_CLUSTER_H
-#define WOLFF_CLUSTER_H
-
-#include <random>
-#include <cmath>
-#include <vector>
-#include <queue>
-
-#include "system.hpp"
-
-namespace wolff {
-
-template <class R_t, class X_t>
-void system<R_t, X_t>::flip_cluster(v_t i0, const R_t& r,
-                  std::mt19937& rng, measurement<R_t, X_t>& A) {
-  std::uniform_real_distribution<double> dist(0.0, 1.0);
-
-  std::queue<v_t> queue;
-  queue.push(i0);
-
-  std::vector<bool> marks(G.nv, false);
-
-  while (!queue.empty()) {
-    v_t i = queue.front();
-    queue.pop();
-
-    if (!marks[i]) { // don't reprocess anyone we've already visited!
-      marks[i] = true;
-
-      X_t si_new;
-#ifndef WOLFF_NO_FIELD
-      R_t s0_new;
-
-      bool we_are_ghost = (i == nv);
-
-      if (we_are_ghost) {
-        s0_new = r.act(s0);
-      } else
-#endif
-      {
-        si_new = r.act(s[i]);
-      }
-
-      for (const v_t &j : G.adjacency[i]) {
-        double dE, p;
-
-#ifndef WOLFF_NO_FIELD
-        bool neighbor_is_ghost = (j == nv);
-
-        if (we_are_ghost || neighbor_is_ghost) {
-          X_t s0s_old, s0s_new;
-          v_t non_ghost;
-
-          if (neighbor_is_ghost) {
-            non_ghost = i;
-            s0s_old = s0.act_inverse(s[i]);
-            s0s_new = s0.act_inverse(si_new);
-          } else {
-            non_ghost = j;
-            s0s_old = s0.act_inverse(s[j]);
-            s0s_new = s0_new.act_inverse(s[j]);
-          }
-
-#ifdef WOLFF_SITE_DEPENDENCE
-          dE = B(non_ghost, s0s_old) - B(non_ghost, s0s_new);
-#else
-          dE = B(s0s_old) - B(s0s_new);
-#endif
-
-#ifdef WOLFF_USE_FINITE_STATES
-          p = Bp[s0s_old.enumerate()][s0s_new.enumerate()];
-#endif
-
-          // run measurement hooks for encountering a ghost bond
-          A.ghost_bond_visited(*this, non_ghost, s0s_old, s0s_new, dE);
-        } else // this is a perfectly normal bond!
-#endif
-        {
-#ifdef WOLFF_BOND_DEPENDENCE
-          dE = Z(i, s[i], j, s[j]) - Z(i, si_new, j, s[j]);
-#else
-          dE = Z(s[i], s[j]) - Z(si_new, s[j]);
-#endif
-
-#ifdef WOLFF_USE_FINITE_STATES
-          p = Zp[s[i].enumerate()][si_new.enumerate()][s[j].enumerate()];
-#endif
-
-          // run measurement hooks for encountering a plain bond
-          A.plain_bond_visited(*this, i, si_new, j, dE);
-        }
-
-#ifndef WOLFF_USE_FINITE_STATES
-        p = 1.0 - exp(-dE / T);
-#endif
-
-        if (dist(rng) < p) {
-          queue.push(j); // push the neighboring vertex to the queue 
-        }
-      }
-
-#ifndef WOLFF_NO_FIELD
-      if (we_are_ghost) {
-        A.ghost_site_transformed(*this, s0_new);
-        s0 = s0_new;
-      } else
-#endif
-      {
-        A.plain_site_transformed(*this, i, si_new);
-        s[i] = si_new;
-      }
-    }
-  }
-}
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp
deleted file mode 100644
index 65b3941..0000000
--- a/lib/include/wolff/graph.hpp
+++ /dev/null
@@ -1,30 +0,0 @@
-
-#ifndef WOLFF_GRAPH_H
-#define WOLFF_GRAPH_H
-
-#include <cmath>
-#include <vector>
-
-namespace wolff {
-
-#include "types.h"
-
-class graph {
-  public:
-    D_t D;
-    L_t L;
-    v_t ne;
-    v_t nv;
-    std::vector<std::vector<v_t>> adjacency;
-    std::vector<std::vector<double>> coordinate;
-
-    graph();
-    graph(D_t D, L_t L);
-
-    void add_ghost();
-};
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/measurement.hpp b/lib/include/wolff/measurement.hpp
deleted file mode 100644
index 6e352a5..0000000
--- a/lib/include/wolff/measurement.hpp
+++ /dev/null
@@ -1,28 +0,0 @@
-
-#ifndef WOLFF_MEASUREMENTS_H
-#define WOLFF_MEASUREMENTS_H
-
-#include "system.hpp"
-
-namespace wolff {
-
-template <class R_t, class X_t>
-class measurement {
-  public:
-    virtual void pre_cluster(N_t, N_t, const system<R_t, X_t>&, v_t, const R_t&) {};
-
-    virtual void plain_bond_visited(const system<R_t, X_t>&, v_t, const X_t&, v_t, double) {};
-    virtual void plain_site_transformed(const system<R_t, X_t>&, v_t, const X_t&) {};
-
-#ifndef WOLFF_NO_FIELD
-    virtual void ghost_bond_visited(const system<R_t, X_t>&, v_t, const X_t&, const X_t&, double) {};
-    virtual void ghost_site_transformed(const system<R_t, X_t>&, const R_t&) {};
-#endif
-
-    virtual void post_cluster(N_t, N_t, const system<R_t, X_t>&) {};
-};
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/models/dihedral.hpp b/lib/include/wolff/models/dihedral.hpp
deleted file mode 100644
index cf4b067..0000000
--- a/lib/include/wolff/models/dihedral.hpp
+++ /dev/null
@@ -1,56 +0,0 @@
-
-#ifndef WOLFF_MODELS_DIHEDRAL_H
-#define WOLFF_MODELS_DIHEDRAL_H
-
-#include "potts.hpp"
-
-namespace wolff {
-
-#include "../types.h"
-
-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);
-      }
-    }
-};
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/models/dihedral_inf.hpp b/lib/include/wolff/models/dihedral_inf.hpp
deleted file mode 100644
index 33d5cec..0000000
--- a/lib/include/wolff/models/dihedral_inf.hpp
+++ /dev/null
@@ -1,57 +0,0 @@
-
-#ifndef WOLFF_MODELS_DIHEDRAL_INF_H
-#define WOLFF_MODELS_DIHEDRAL_INF_H
-
-#include <cmath>
-#include "height.hpp"
-
-namespace wolff {
-
-#include "../types.h"
-
-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);
-    }
-  }
-};
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/models/height.hpp b/lib/include/wolff/models/height.hpp
deleted file mode 100644
index 9d7e87c..0000000
--- a/lib/include/wolff/models/height.hpp
+++ /dev/null
@@ -1,57 +0,0 @@
-
-#pragma once
-
-#include <cmath>
-#include <wolff/types.h>
-
-namespace wolff {
-
-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 typename height_t<T>::M_t operator*(v_t a, height_t<T> h) {
-  return h * a;
-}
-
-template <class T>
-inline typename height_t<T>::F_t operator*(double a, height_t<T> h) {
-  return h * a;
-}
-
-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
deleted file mode 100644
index d2f6f91..0000000
--- a/lib/include/wolff/models/ising.hpp
+++ /dev/null
@@ -1,93 +0,0 @@
-
-#ifndef WOLFF_MODELS_ISING_H
-#define WOLFF_MODELS_ISING_H
-
-#define WOLFF_FINITE_STATES_N 2
-
-#include "../system.hpp"
-
-namespace wolff {
-
-#include "../types.h"
-
-class ising_t {
-  public:
-    bool x;
-
-    ising_t() : x(false) {}
-
-    ising_t(bool x) : x(x) {} 
-    ising_t(q_t 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;
-      }
-    }
-
-    q_t enumerate() const {
-      return (q_t)x;
-    }
-};
-
-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;
-}
-
-ising_t gen_ising(std::mt19937&, const system<ising_t, ising_t>&, v_t) {
-  return ising_t(true);
-};
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/models/orthogonal.hpp b/lib/include/wolff/models/orthogonal.hpp
deleted file mode 100644
index 514c88a..0000000
--- a/lib/include/wolff/models/orthogonal.hpp
+++ /dev/null
@@ -1,208 +0,0 @@
-
-#ifndef WOLFF_MODELS_ORTHOGONAL_H
-#define WOLFF_MODELS_ORTHOGONAL_H
-
-#include <random>
-#include <cmath>
-
-#include "../system.hpp"
-#include "vector.hpp"
-
-namespace wolff {
-
-#include "../types.h"
-
-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 system<orthogonal_t<q, double>, vector_t<q, double>>&, v_t) {
-  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 system<orthogonal_t<q, double>, vector_t<q, double>>& S, v_t i0, 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] = S.s[i0][0] * cosr - S.s[i0][1] * sinr;
-    v[1] = S.s[i0][1] * cosr + S.s[i0][0] * sinr;
-
-    for (q_t i = 2; i < q; i++) {
-      v[i] = S.s[i0][i];
-    }
-  } else {
-    v = S.s[i0];
-  }
-
-  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;
-}
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/models/potts.hpp b/lib/include/wolff/models/potts.hpp
deleted file mode 100644
index a0cb368..0000000
--- a/lib/include/wolff/models/potts.hpp
+++ /dev/null
@@ -1,68 +0,0 @@
-
-#ifndef WOLFF_MODELS_POTTS_H
-#define WOLFF_MODELS_POTTS_H
-
-#include <cmath>
-
-#include "vector.hpp"
-
-namespace wolff {
-
-#include "../types.h"
-
-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;
-    }
-
-    q_t enumerate() const {
-      return x;
-    }
-};
-
-template<q_t q>
-inline typename potts_t<q>::M_t operator*(v_t a, const potts_t<q>& s) {
-  return s * a;
-}
-
-template<q_t q>
-inline typename potts_t<q>::F_t operator*(double a, const potts_t<q>& s) {
-  return s * a;
-}
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/models/symmetric.hpp b/lib/include/wolff/models/symmetric.hpp
deleted file mode 100644
index d15a7ba..0000000
--- a/lib/include/wolff/models/symmetric.hpp
+++ /dev/null
@@ -1,59 +0,0 @@
-
-#ifndef WOLFF_MODELS_SYMMETRIC_H
-#define WOLFF_MODELS_SYMMETRIC_H
-
-#include <array>
-
-#include "potts.hpp"
-
-namespace wolff {
-
-#include "../types.h"
-
-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;
-    }
-};
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/models/vector.hpp b/lib/include/wolff/models/vector.hpp
deleted file mode 100644
index e4c4a1c..0000000
--- a/lib/include/wolff/models/vector.hpp
+++ /dev/null
@@ -1,115 +0,0 @@
-
-#ifndef WOLFF_MODELS_VECTOR_H
-#define WOLFF_MODELS_VECTOR_H
-
-#include <cmath>
-#include <array>
-#include <iostream>
-
-namespace wolff {
-
-#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;
-}
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/system.hpp b/lib/include/wolff/system.hpp
deleted file mode 100644
index fa9a1b0..0000000
--- a/lib/include/wolff/system.hpp
+++ /dev/null
@@ -1,105 +0,0 @@
-
-#ifndef WOLFF_SYSTEM_H
-#define WOLFF_SYSTEM_H
-
-#include <functional>
-#include <vector>
-#include <random>
-
-#include "graph.hpp"
-
-namespace wolff {
-
-#include "types.h"
-
-template <class R_t, class X_t>
-class measurement;
-
-template <class R_t, class X_t>
-class system {
-  public:
-    v_t nv; // number of vertices
-    v_t ne; // number of edges
-    graph G; // the graph defining the lattice with ghost
-    double T; // the temperature
-    std::vector<X_t> s; // the state of the ordinary spins
-
-#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
-    R_t s0; // the current state of the ghost site
-#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
-
-#ifdef WOLFF_USE_FINITE_STATES
-    std::array<std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Zp;
-#ifndef WOLFF_NO_FIELD
-    std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Bp;
-#endif
-#endif
-
-    system(graph g, 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
-        ) : G(g), 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_ghost();
-#endif
-#ifdef WOLFF_USE_FINITE_STATES
-      this->finite_states_init();
-#endif
-    }
-
-    void flip_cluster(v_t, const R_t&, std::mt19937&, measurement<R_t, X_t>&);
-    void run_wolff(N_t, std::function <R_t(std::mt19937&, const system<R_t, X_t>&, v_t)> r_gen, measurement<R_t, X_t>& A, std::mt19937& rng);
-
-#ifdef WOLFF_USE_FINITE_STATES
-    void finite_states_init() {
-#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++) {
-          Bp[i][j] = 1.0 - exp(-(B(X_t(i)) - B(X_t(j))) / T);
-        }
-      }
-#endif
-      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++) {
-            Zp[i][j][k] = 1.0 - exp(-(Z(X_t(i), X_t(k)) - Z(X_t(j), X_t(k))) / T);
-          }
-        }
-      }
-    }
-#endif
-
-};
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/types.h b/lib/include/wolff/types.h
deleted file mode 100644
index 5bbaa6d..0000000
--- a/lib/include/wolff/types.h
+++ /dev/null
@@ -1,32 +0,0 @@
-
-#ifndef WOLFF_TYPES_H
-#define WOLFF_TYPES_H
-
-#include <inttypes.h>
-
-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_D UINT_FAST8_MAX
-#define MAX_L UINT_FAST16_MAX
-#define MAX_N UINT_FAST64_MAX
-
-#define PRIv PRIuFAST32
-#define PRIq PRIuFAST8
-#define PRID PRIuFAST8
-#define PRIL PRIuFAST16
-#define PRIN PRIuFAST64
-
-#define SCNv SCNuFAST32
-#define SCNq SCNuFAST8
-#define SCND SCNuFAST8
-#define SCNL SCNuFAST16
-#define SCNN SCNuFAST64
-
-#endif
-
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
deleted file mode 100644
index 4f91be4..0000000
--- a/lib/src/graph.cpp
+++ /dev/null
@@ -1,55 +0,0 @@
-
-#include <wolff/graph.hpp>
-
-namespace wolff {
-
-graph::graph() {
-  D = 0;
-  L = 0;
-  nv = 0;
-  ne = 0;
-}
-
-graph::graph(D_t D, L_t L) : D(D), L(L) {
-  nv = pow(L, D);
-  ne = D * nv;
-
-  adjacency.resize(nv);
-  coordinate.resize(nv);
-
-  for (std::vector<v_t> adj_i : adjacency) {
-    adj_i.reserve(2 * D);
-  }
-
-  for (v_t i = 0; i < nv; i++) {
-    coordinate[i].resize(D);
-    for (D_t j = 0; j < D; j++) {
-      coordinate[i][j] = (i / (v_t)pow(L, D - j - 1)) % L;
-
-      adjacency[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1)));
-      adjacency[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)));
-    }
-  }
-}
-
-void graph::add_ghost() {
-  for (std::vector<v_t>& adj_i : adjacency) {
-    adj_i.push_back(nv);
-  }
-
-  adjacency.resize(nv + 1);
-  coordinate.resize(nv + 1);
-  adjacency[nv].reserve(nv);
-
-  for (v_t i = 0; i < nv; i++) {
-    adjacency[nv].push_back(i);
-  }
-
-  coordinate[nv].resize(coordinate[0].size());
-
-  ne += nv;
-  nv++;
-}
-
-}
-
diff --git a/lib/wolff.hpp b/lib/wolff.hpp
new file mode 100644
index 0000000..d21d788
--- /dev/null
+++ b/lib/wolff.hpp
@@ -0,0 +1,320 @@
+
+#ifndef WOLFF_H
+#define WOLFF_H
+
+#include <functional>
+#include <vector>
+#include <random>
+#include <cmath>
+#include <iterator>
+#include <list>
+#include <tuple>
+#include <queue>
+
+namespace wolff{
+
+  template <class vertex_prop = std::tuple<>, class edge_prop = std::tuple<>>
+  class graph {
+    /* the graph class describes a lattice on which wolff is run. for most
+     * purposes, using the default square lattice constructor is sufficient,
+     * but arbitrary lattices can be constructed
+     */
+    public:
+      unsigned D;  // dimension of space
+      unsigned L;  // linear size
+      unsigned ne; // number of edges
+      unsigned nv; // number of vertices
+
+      struct _vertex;
+
+      typedef struct _halfedge {
+        struct _vertex &self;     // reference to the vertex this edge comes from
+        struct _vertex &neighbor; // reference to the vertex this edge goes to
+        edge_prop prop;           // optional properties
+
+        _halfedge(struct _vertex &v1, struct _vertex &v2) : self(v1), neighbor(v2) {};
+      } halfedge;
+
+      typedef struct _vertex {
+        unsigned ind;              // index of the vertex
+        std::list<halfedge> edges; // list of edges incident on this vertex
+        vertex_prop prop;          // optional properties
+      } vertex;
+
+      std::vector<vertex> vertices; /* vector of vertices, length nv, with 
+                                     * vertices[i].ind = i for all 0 <= i < nv
+                                     */
+      graph() {
+        // default constructor for empty graph. use it to build your own!
+        D = 0;
+        L = 0;
+        nv = 0;
+        ne = 0;
+      };
+
+      graph(unsigned D, unsigned L) : D(D), L(L) {
+        // default constructor for square lattice graph
+        nv = pow(L, D);
+        ne = D * nv;
+
+        vertices.resize(nv);
+
+        for (unsigned i = 0; i < nv; i++) {
+          vertices[i].ind = i;
+          for (unsigned j = 0; j < D; j++) {
+            unsigned n1 = pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) +
+                                              fmod(i + pow(L, j), pow(L, j + 1));
+            unsigned n2 = pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) +
+                              fmod(pow(L, j + 1) + i - pow(L, j), pow(L, j + 1));
+
+            halfedge f(vertices[i], vertices[n1]);
+            halfedge b(vertices[i], vertices[n2]);
+
+            vertices[i].edges.push_back(f);
+            vertices[i].edges.push_back(b);
+          }
+        }
+      };
+
+      void add_ghost() {
+        // adds a ghost site to any graph
+        vertices.resize(nv + 1);
+
+        for (auto it = vertices.begin(); it != std::prev(vertices.end()); it++) {
+          halfedge e1(*it, vertices[nv]);
+          it->edges.push_back(e1);
+
+          halfedge e2(vertices[nv], *it);
+          vertices[nv].edges.push_back(e2);
+        }
+
+        vertices[nv].ind = nv;
+
+        ne += nv;
+        nv++;
+      };
+  };
+
+  template <class R_t, class X_t, class G_t = graph<>>
+  class measurement;
+
+  template <class R_t, class X_t, class G_t = graph<>>
+  class system {
+    public:
+      unsigned nv; // number of vertices
+      unsigned ne; // number of edges
+      G_t G; // the graph defining the lattice with ghost
+      double T; // the temperature
+      std::vector<X_t> s; // the state of the ordinary spins
+
+#ifdef WOLFF_BOND_DEPENDENCE
+      std::function <double(const typename G_t::halfedge&, const X_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
+      R_t s0; // the current state of the ghost site
+#ifdef WOLFF_SITE_DEPENDENCE
+      std::function <double(const typename G_t::vertex&, const X_t&)> B; // coupling with the external field
+#else
+      std::function <double(const X_t&)> B; // coupling with the external field
+#endif
+#endif
+
+#ifdef WOLFF_USE_FINITE_STATES
+      std::array<std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Zp;
+#ifndef WOLFF_NO_FIELD
+      std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Bp;
+#endif
+#endif
+
+      system(G_t g, double T,
+#ifdef WOLFF_BOND_DEPENDENCE
+          std::function <double(const typename G_t::halfedge&, const X_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(const typename G_t::vertex&, const X_t&)> B
+#else
+          , std::function <double(const X_t&)> B
+#endif
+#endif
+          ) : G(g), 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_ghost();
+#endif
+#ifdef WOLFF_USE_FINITE_STATES
+        this->finite_states_init();
+#endif
+      }
+
+      void flip_cluster(typename G_t::vertex& v0, const R_t& r, std::mt19937& rng,
+                          measurement<R_t, X_t, G_t>& A) {
+        std::uniform_real_distribution<double> dist(0.0, 1.0);
+
+        std::queue<unsigned> queue;
+        queue.push(v0.ind);
+
+        std::vector<bool> visited(G.nv, false);
+
+        while (!queue.empty()) {
+          unsigned i = queue.front();
+          queue.pop();
+
+          if (!visited[i]) { // don't reprocess anyone we've already visited!
+            visited[i] = true;
+
+            X_t si_new;
+#ifndef WOLFF_NO_FIELD
+            R_t s0_new;
+
+            bool we_are_ghost = (i == nv);
+
+            if (we_are_ghost) {
+              s0_new = r.act(s0);
+            } else
+#endif
+            {
+              si_new = r.act(s[i]);
+            }
+
+            for (const typename G_t::halfedge &e : G.vertices[i].edges) {
+              double dE, p;
+              unsigned j = e.neighbor.ind;
+
+#ifndef WOLFF_NO_FIELD
+              bool neighbor_is_ghost = (j == nv);
+
+              if (we_are_ghost || neighbor_is_ghost) {
+                X_t s0s_old, s0s_new;
+                unsigned non_ghost;
+
+                if (neighbor_is_ghost) {
+                  non_ghost = i;
+                  s0s_old = s0.act_inverse(s[i]);
+                  s0s_new = s0.act_inverse(si_new);
+                } else {
+                  non_ghost = j;
+                  s0s_old = s0.act_inverse(s[j]);
+                  s0s_new = s0_new.act_inverse(s[j]);
+                }
+
+#ifdef WOLFF_SITE_DEPENDENCE
+                dE = B(G.vertices[non_ghost], s0s_old) - B(G.vertices[non_ghost], s0s_new);
+#else
+                dE = B(s0s_old) - B(s0s_new);
+#endif
+
+#ifdef WOLFF_USE_FINITE_STATES
+                p = Bp[s0s_old.enumerate()][s0s_new.enumerate()];
+#endif
+
+                // run measurement hooks for encountering a ghost bond
+                A.ghost_bond_visited(*this, G.vertices[non_ghost], s0s_old, s0s_new, dE);
+              } else // this is a perfectly normal bond!
+#endif
+              {
+#ifdef WOLFF_BOND_DEPENDENCE
+                dE = Z(e, s[i], s[j]) - Z(e, si_new, s[j]);
+#else
+                dE = Z(s[i], s[j]) - Z(si_new, s[j]);
+#endif
+
+#ifdef WOLFF_USE_FINITE_STATES
+                p = Zp[s[i].enumerate()][si_new.enumerate()][s[j].enumerate()];
+#endif
+
+                // run measurement hooks for encountering a plain bond
+                A.plain_bond_visited(*this, e, si_new, dE);
+              }
+
+#ifndef WOLFF_USE_FINITE_STATES
+              p = 1.0 - exp(-dE / T);
+#endif
+
+              if (dist(rng) < p) {
+                queue.push(j); // push the neighboring vertex to the queue 
+              }
+            }
+
+#ifndef WOLFF_NO_FIELD
+            if (we_are_ghost) {
+              A.ghost_site_transformed(*this, s0_new);
+              s0 = s0_new;
+            } else
+#endif
+            {
+              A.plain_site_transformed(*this, G.vertices[i], si_new);
+              s[i] = si_new;
+            }
+          }
+        }
+      }
+
+      void run_wolff(unsigned N,
+          std::function <R_t(std::mt19937&, const system<R_t, X_t, G_t>&, const typename G_t::vertex&)> r_gen, measurement<R_t, X_t, G_t>& A, std::mt19937& rng) {
+        std::uniform_int_distribution<unsigned> dist(0, nv - 1);
+
+        for (unsigned n = 0; n < N; n++) {
+          unsigned i0 = dist(rng);
+          R_t r = r_gen(rng, *this, G.vertices[i0]);
+
+          A.pre_cluster(n, N, *this, G.vertices[i0], r);
+
+          this->flip_cluster(G.vertices[i0], r, rng, A);
+
+          A.post_cluster(n, N, *this);
+        }
+      }
+
+#ifdef WOLFF_USE_FINITE_STATES
+      void finite_states_init() {
+#ifndef WOLFF_NO_FIELD
+        for (unsigned i = 0; i < WOLFF_FINITE_STATES_N; i++) {
+          for (unsigned j = 0; j < WOLFF_FINITE_STATES_N; j++) {
+            Bp[i][j] = 1.0 - exp(-(B(X_t(i)) - B(X_t(j))) / T);
+          }
+        }
+#endif
+        for (unsigned i = 0; i < WOLFF_FINITE_STATES_N; i++) {
+          for (unsigned j = 0; j < WOLFF_FINITE_STATES_N; j++) {
+            for (unsigned k = 0; k < WOLFF_FINITE_STATES_N; k++) {
+              Zp[i][j][k] = 1.0 - exp(-(Z(X_t(i), X_t(k)) - Z(X_t(j), X_t(k))) / T);
+            }
+          }
+        }
+      }
+#endif
+
+  };
+
+  template <class R_t, class X_t, class G_t>
+  class measurement {
+    public:
+      virtual void pre_cluster(unsigned, unsigned, const system<R_t, X_t, G_t>&, const typename G_t::vertex& v, const R_t&) {};
+
+      virtual void plain_bond_visited(const system<R_t, X_t, G_t>&, const typename G_t::halfedge& e, const X_t&, double) {};
+      virtual void plain_site_transformed(const system<R_t, X_t, G_t>&, const typename G_t::vertex& v, const X_t&) {};
+
+#ifndef WOLFF_NO_FIELD
+      virtual void ghost_bond_visited(const system<R_t, X_t, G_t>&, const typename G_t::vertex& v, const X_t&, const X_t&, double) {};
+      virtual void ghost_site_transformed(const system<R_t, X_t, G_t>&, const R_t&) {};
+#endif
+
+      virtual void post_cluster(unsigned, unsigned, const system<R_t, X_t, G_t>&) {};
+  };
+
+}
+
+#endif
+
diff --git a/lib/wolff_models/dihedral.hpp b/lib/wolff_models/dihedral.hpp
new file mode 100644
index 0000000..aef3be6
--- /dev/null
+++ b/lib/wolff_models/dihedral.hpp
@@ -0,0 +1,54 @@
+
+#ifndef WOLFF_MODELS_DIHEDRAL_H
+#define WOLFF_MODELS_DIHEDRAL_H
+
+#include "potts.hpp"
+
+namespace wolff {
+
+  template <unsigned q>
+  class dihedral_t {
+    public:
+      bool is_reflection;
+      unsigned x;
+
+      dihedral_t() : is_reflection(false), x(0) {}
+      dihedral_t(bool x, unsigned 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);
+        }
+      }
+  };
+
+}
+
+#endif
+
diff --git a/lib/wolff_models/dihedral_inf.hpp b/lib/wolff_models/dihedral_inf.hpp
new file mode 100644
index 0000000..0b09c59
--- /dev/null
+++ b/lib/wolff_models/dihedral_inf.hpp
@@ -0,0 +1,55 @@
+
+#ifndef WOLFF_MODELS_DIHEDRAL_INF_H
+#define WOLFF_MODELS_DIHEDRAL_INF_H
+
+#include <cmath>
+#include "height.hpp"
+
+namespace wolff {
+
+  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);
+      }
+    }
+  };
+
+}
+
+#endif
+
diff --git a/lib/wolff_models/height.hpp b/lib/wolff_models/height.hpp
new file mode 100644
index 0000000..ae3ad82
--- /dev/null
+++ b/lib/wolff_models/height.hpp
@@ -0,0 +1,56 @@
+
+#pragma once
+
+#include <cmath>
+
+namespace wolff {
+
+  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*(unsigned 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 typename height_t<T>::M_t operator*(unsigned a, height_t<T> h) {
+    return h * a;
+  }
+
+  template <class T>
+  inline typename height_t<T>::F_t operator*(double a, height_t<T> h) {
+    return h * a;
+  }
+
+  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/wolff_models/ising.hpp b/lib/wolff_models/ising.hpp
new file mode 100644
index 0000000..345eeff
--- /dev/null
+++ b/lib/wolff_models/ising.hpp
@@ -0,0 +1,92 @@
+
+#ifndef WOLFF_MODELS_ISING_H
+#define WOLFF_MODELS_ISING_H
+
+#define WOLFF_FINITE_STATES_N 2
+
+#include "wolff.hpp"
+
+namespace wolff {
+
+  class ising_t {
+    public:
+      bool x;
+
+      ising_t() : x(false) {}
+
+      ising_t(bool x) : x(x) {} 
+      ising_t(unsigned 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*(unsigned 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;
+        }
+      }
+
+      unsigned enumerate() const {
+        return (unsigned)x;
+      }
+  };
+
+  inline ising_t::M_t operator*(unsigned a, const ising_t& s) {
+    return s * a;
+  }
+
+  inline ising_t::F_t operator*(double a, const ising_t& s) {
+    return s * a;
+  }
+
+  template <class G_t>
+  ising_t gen_ising(std::mt19937&, const system<ising_t, ising_t, G_t>&, const typename G_t::vertex&) {
+    return ising_t(true);
+  };
+
+}
+
+#endif
+
diff --git a/lib/wolff_models/orthogonal.hpp b/lib/wolff_models/orthogonal.hpp
new file mode 100644
index 0000000..ea9af14
--- /dev/null
+++ b/lib/wolff_models/orthogonal.hpp
@@ -0,0 +1,206 @@
+
+#ifndef WOLFF_MODELS_ORTHOGONAL_H
+#define WOLFF_MODELS_ORTHOGONAL_H
+
+#include <random>
+#include <cmath>
+
+#include <wolff.hpp>
+#include "vector.hpp"
+
+namespace wolff {
+
+  template <unsigned q, class T>
+  class orthogonal_t : public std::array<std::array<T, q>, q> {
+    public :
+    bool is_reflection;
+
+    orthogonal_t() : is_reflection(false) {
+      for (unsigned 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 (unsigned i = 0; i < q; i++) {
+          prod += v[i] * (*this)[0][i];
+        }
+        for (unsigned i = 0; i < q; i++) {
+          v_rot[i] = v[i] - 2 * prod * (*this)[0][i];
+        }
+      } else {
+        for (unsigned i = 0; i < q; i++) {
+          for (unsigned 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 (unsigned i = 0; i < q; i++) {
+          double akOki = 0;
+
+          for (unsigned k = 0; k < q; k++) {
+            akOki += (*this)[0][k] * m[k][i];
+          }
+
+          for (unsigned j = 0; j < q; j++) {
+            m_rot[j][i] = m[j][i] - 2 * akOki * (*this)[0][j];
+          }
+        }
+      } else {
+        for (unsigned i = 0; i < q; i++) {
+          m_rot[i].fill(0);
+          for (unsigned j = 0; j < q; j++) {
+            for (unsigned 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 (unsigned i = 0; i < q; i++) {
+          for (unsigned 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 (unsigned i = 0; i < q; i++) {
+          m_rot[i].fill(0);
+          for (unsigned j = 0; j < q; j++) {
+            for (unsigned k = 0; k < q; k++) {
+              m_rot[i][j] += (*this)[j][i] * m[j][k];
+            }
+          }
+        }
+
+        return m_rot;
+      }
+    }
+
+  };
+
+  template <unsigned q, class G_t>
+  orthogonal_t <q, double> generate_rotation_uniform (std::mt19937& r, const system<orthogonal_t<q, double>, vector_t<q, double>, G_t>&, const typename G_t::vertex&) {
+    std::normal_distribution<double> dist(0.0,1.0);
+    orthogonal_t <q, double> ptr;
+    ptr.is_reflection = true;
+
+    double v2 = 0;
+
+    for (unsigned i = 0; i < q; i++) {
+      ptr[0][i] = dist(r);
+      v2 += ptr[0][i] * ptr[0][i];
+    }
+
+    double mag_v = sqrt(v2);
+
+    for (unsigned i = 0; i < q; i++) {
+      ptr[0][i] /= mag_v;
+    }
+
+    return ptr;
+  }
+
+  template <unsigned q, class G_t>
+  orthogonal_t <q, double> generate_rotation_perturbation (std::mt19937& r, const system<orthogonal_t<q, double>, vector_t<q, double>, G_t>& S, const typename G_t::vertex& 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] = S.s[v0.ind][0] * cosr - S.s[v0.ind][1] * sinr;
+      v[1] = S.s[v0.ind][1] * cosr + S.s[v0.ind][0] * sinr;
+
+      for (unsigned i = 2; i < q; i++) {
+        v[i] = S.s[v0.ind][i];
+      }
+    } else {
+      v = S.s[v0.ind];
+    }
+
+    double m_dot_v = 0;
+
+    for (unsigned 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 (unsigned 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 (unsigned i = 0; i < q; i++) {
+      m[0][i] /= mag_v; // normalize
+    }
+
+    v2 = 0;
+
+    double factor = epsilon * dist(r);
+
+    for (unsigned 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 (unsigned i = 0; i < q; i++) {
+      m[0][i] /= mag_v; // normalize
+    }
+
+    return m;
+  }
+
+}
+
+#endif
+
diff --git a/lib/wolff_models/potts.hpp b/lib/wolff_models/potts.hpp
new file mode 100644
index 0000000..dd56134
--- /dev/null
+++ b/lib/wolff_models/potts.hpp
@@ -0,0 +1,66 @@
+
+#ifndef WOLFF_MODELS_POTTS_H
+#define WOLFF_MODELS_POTTS_H
+
+#include <cmath>
+
+#include "vector.hpp"
+
+namespace wolff {
+
+  template <unsigned q>
+  class potts_t {
+    public:
+      unsigned x;
+
+      potts_t() : x(0) {}
+      potts_t(unsigned x) : x(x) {}
+
+      typedef vector_t<q, int> M_t;
+      typedef vector_t<q, double> F_t;
+
+      inline vector_t<q, int> operator*(unsigned 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;
+      }
+
+      unsigned enumerate() const {
+        return x;
+      }
+  };
+
+  template<unsigned q>
+  inline typename potts_t<q>::M_t operator*(unsigned a, const potts_t<q>& s) {
+    return s * a;
+  }
+
+  template<unsigned q>
+  inline typename potts_t<q>::F_t operator*(double a, const potts_t<q>& s) {
+    return s * a;
+  }
+
+}
+
+#endif
+
diff --git a/lib/wolff_models/symmetric.hpp b/lib/wolff_models/symmetric.hpp
new file mode 100644
index 0000000..98e53e3
--- /dev/null
+++ b/lib/wolff_models/symmetric.hpp
@@ -0,0 +1,57 @@
+
+#ifndef WOLFF_MODELS_SYMMETRIC_H
+#define WOLFF_MODELS_SYMMETRIC_H
+
+#include <array>
+
+#include "potts.hpp"
+
+namespace wolff {
+
+  template <unsigned q>
+  class symmetric_t : public std::array<unsigned, q> {
+    public:
+
+      symmetric_t() {
+        for (unsigned 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 (unsigned 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 (unsigned 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 (unsigned i = 0; i < q; i++) {
+          r_rot[(*this)[i]] = r[i];
+        }
+
+        return r_rot;
+      }
+  };
+
+}
+
+#endif
+
diff --git a/lib/wolff_models/vector.hpp b/lib/wolff_models/vector.hpp
new file mode 100644
index 0000000..4106e62
--- /dev/null
+++ b/lib/wolff_models/vector.hpp
@@ -0,0 +1,113 @@
+
+#ifndef WOLFF_MODELS_VECTOR_H
+#define WOLFF_MODELS_VECTOR_H
+
+#include <cmath>
+#include <array>
+#include <iostream>
+
+namespace wolff {
+
+  template <unsigned 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 (unsigned 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 (unsigned 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 (unsigned i = 0; i < q; i++) {
+          (*this)[i] -= (T)v[i];
+        }
+        return *this;
+      }
+
+      inline vector_t<q, T> operator*(unsigned x) const {
+        vector_t<q, T> result;
+        for (unsigned 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 (unsigned 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 (unsigned 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 (unsigned i = 0; i < q; i++) {
+          result[i] = (*this)[i] / a;
+        }
+
+        return result;
+      }
+  };
+
+  template<unsigned q, class T>
+  inline vector_t<q, T> operator*(unsigned a, const vector_t<q, T>&v) {
+    return v * a;
+  }
+
+  template<unsigned q, class T>
+  inline vector_t<q, double> operator*(double a, const vector_t<q, T>&v) {
+    return v * a;
+  }
+
+  template<unsigned 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;
+  }
+
+}
+
+#endif
+