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.

44 files changed, 1297 insertions, 1419 deletions

diff --git a/doc/compile.rst b/doc/compile.rst
index e4960a3..d335dff 100644
--- a/doc/compile.rst
+++ b/doc/compile.rst
@@ -23,7 +23,7 @@ Building With Bond Dependence
 
    When this flag is defined Wolff will ask for the indices of the spins on each side a bond, allowing the implementation of random bonds or anisotropic interactions.
 
-   When defined, the bond coupling must be a function of the form :cpp:any:`double Z(v_t, const X_t&, v_t, const X_t&)`, where the first argument is the index of the first spin and the third argument is the index of the second spin. A function of this type is passed to :cpp:class:`system` in place of the original bond coupling.
+   When defined, the bond coupling must be a function of the form :cpp:any:`double Z(const G_t::halfedge&, const X_t&, const X_t&)`, where the first argument is the edge being considered. A function of this type is passed to :cpp:class:`system` in place of the original bond coupling.
 
 
 Building With Site Dependence
@@ -33,5 +33,6 @@ Building With Site Dependence
 
    When this flag is defined Wolff will ask for the indices of the spin when measuring the external field, allowing the implementation of random fields or to emulate boundaries.
 
-   When defined, the field coupling must be a function of the form :cpp:any:`double B(v_t, const X_t&)`, where the first argument is the index of the spin. A function of this type is passed to :cpp:class:`system` in place of the original field coupling.
+   When defined, the field coupling must be a function of the form :cpp:any:`double B(const G_t::vertex&, const X_t&)`, where the first argument is the vertex the spin is on. A function of this type is passed to :cpp:class:`system` in place of the original field coupling.
 
+   An example of a system of this type can be found in :file:`examples/ising_random_field.cpp`, which uses a non-trivial vertex property to communicate vertex dependence to the field coupling function.
diff --git a/doc/finite_states.rst b/doc/finite_states.rst
index c36881e..66e34d2 100644
--- a/doc/finite_states.rst
+++ b/doc/finite_states.rst
@@ -5,7 +5,7 @@
 Finite States
 *************
 
-One of the slower steps in running the algorithm is taking the exponent involved in computing the bond activation probabilities for each prospective bond. When the spins in your system have a finite number of possible states, the algorithm can be sped up considerably by precomputing the bond activation probabilities for every possible pair of spins. Once the appropriate things have been defined for your model, the compile definition :c:macro:`WOLFF_USE_FINITE_STATES` can be set to automate this process. The provided model headers :file:`wolff/models/ising.hpp` and :file:`wolff/models/potts.hpp` demonstrate the expected usage.
+One of the slower steps in running the algorithm is taking the exponent involved in computing the bond activation probabilities for each prospective bond. When the spins in your system have a finite number of possible states, the algorithm can be sped up considerably by precomputing the bond activation probabilities for every possible pair of spins. Once the appropriate things have been defined for your model, the compile definition :c:macro:`WOLFF_USE_FINITE_STATES` can be set to automate this process. The provided model headers :file:`wolff_models/ising.hpp` and :file:`wolff_models/potts.hpp` demonstrate the expected usage.
 
 Required Definitions
 ====================
diff --git a/doc/graph.rst b/doc/graph.rst
index 8a96e50..dfe684c 100644
--- a/doc/graph.rst
+++ b/doc/graph.rst
@@ -1,50 +1,83 @@
 
 .. default-domain:: cpp
+.. _graphs:
 
 ******
 Graphs
 ******
 
-This class is defined in the header :file:`wolff/graph.hpp`.
+This class is defined in the header :file:`lib/wolff.hpp`.
 
-.. class:: graph
+.. class:: \template <class vertex_prop = std::tuple\<>, class edge_prop = std::tuple\<>> graph
 
-   Lattices are described by objects of class :class:`graph`, a minimal implementation of graphs.
+   Lattices are described by objects of class :class:`graph`, a minimal implementation of graphs. Can be called with `graph<>` if no properties need to be associated with vertices or edges. Otherwise, those properties can be supplied as classes.
 
-   .. member:: D_t graph::D
+   .. member:: unsigned graph::D
 
       The dimension of the graph. This property is unused by the core library.
 
-   .. member:: L_t graph::L
+   .. member:: unsigned graph::L
 
       The linear size of the graph. This property is unused by the core library.
 
-   .. member:: v_t graph::ne
+   .. member:: unsigned graph::ne
 
       The number of edges in the graph. This property is unused by the core library.
 
-   .. member:: v_t graph::nv
+   .. member:: unsigned graph::nv
 
       The number of vertices in the graph.
 
-   .. member:: std::vector<std::vector<v_t>> graph::adjacency
+   .. class:: graph::vertex
 
-      The adjacency list for the graph. The :math:`i\text{th}` element of :member:`adjacency` is a standard library vector containing the indices of all vertices adjacent to vertex :math:`i`.
+      This class describes the vertices on the graph.
 
-   .. member:: std::vector<std::vector<double>> graph::coordinate
+      .. member:: unsigned ind
 
-      The coordinates of the graph vertices. The :math:`i\text{th}` element of :var:`coordinate` is a standard library vector of length :var:`D` containing the spatial coordinates of vertex :math:`i`. This property is unused by the core library.
+         The index of the vertex, which is also its position in the list :var:`vertices`.
+
+      .. member:: std::list<halfedge> edges
+
+         The list of edges emanating from this vertex.
+
+      .. member:: vertex_prop prop
+
+         Template-defined class which stores optional properties of the vertex.
+
+   .. class:: graph::halfedge
+
+      This class describes the halfedges on the graph.
+
+      .. member:: vertex& self
+
+         A reference to the vertex this halfedge is emanating from.
+
+      .. member:: vertex& neighbor
+
+         A reference to the vertex this halfedge is going to.
+
+      .. member:: edge_prop prop
+
+         Template-defined class which stores optional properties of the edge.
+
+      .. function:: halfedge(vertex &self, vertex &neighbor)
+
+         Constructor which sets self and neighbor from supplied vertices.
+
+   .. member:: std::vector<vertex> graph::vertices
+
+      A list of all vertices in the graph.
 
    .. function:: graph::graph()
 
-      The default constructor. Initializes an empty graph, i.e., :var:`D`, :var:`L`, :var:`ne`, and :var:`nv` are all zero and :var:`adjacency` and :var:`coordinate` are uninitialized.
+      The default constructor. Initializes an empty graph, i.e., :var:`D`, :var:`L`, :var:`ne`, and :var:`nv` are all zero and :var:`vertices` is uninitialized.
 
-   .. function:: graph::graph(D_t D, L_t L)
+   .. function:: graph::graph(unsigned D, unsigned L)
 
       Initializes a graph of a :var:`D`-dimensional hypercubic lattice with :var:`L` vertices per side. This is the only nontrivial graph constructor supplied by the core library. The library will work with arbitrary graphs, and if a different lattice is needed consider calling the default constructor and populating the member objects youself before handing the graph to the :class:`system` constructor.
 
-      :param D_t D: The dimension of space.
-      :param L_t L: The number of vertices per edge of the hypercube.
+      :param unsigned D: The dimension of space.
+      :param unsigned L: The number of vertices per edge of the hypercube.
 
    .. function:: void graph::add_ghost()
 
diff --git a/doc/index.rst b/doc/index.rst
index 7b156d3..c379dd7 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -11,7 +11,6 @@ Wolff Library
    system
    measurement
    graph
-   types
    finite_states
    compile
 
diff --git a/doc/measurement.rst b/doc/measurement.rst
index 1b4381c..1ca41cc 100644
--- a/doc/measurement.rst
+++ b/doc/measurement.rst
@@ -5,64 +5,63 @@
 Measurements
 ************
 
-One of the most complicated parts of using this library is setting up measurements to take during a run, but once understood they provide a powerful way to quickly measure arbitrary properties. The class is defined in the header :file:`wolff/measurement.hpp`.
+One of the most complicated parts of using this library is setting up measurements to take during a run, but once understood they provide a powerful way to quickly measure arbitrary properties. The class is defined in the header :file:`lib/wolff.hpp`.
 
-.. class:: template<R_t, X_t> measurement
+.. class:: template<R_t, X_t, G_t = graph<>> measurement
 
    This class defines several virtual member functions. To use the library, you must supply your own measurement class that inherits this one and defines those functions, which may be trivial.
 
    Each member function defines a hook that is run at a specific time during execution. These hooks can be used to modify member objects of your inheritor measurement class, and thereby extract information from the simulation.
 
-   .. function:: void measurement::pre_cluster(N_t n, N_t N, const system<R_t, X_t>& S, v_t i0, const R_t& r)
+   .. function:: void measurement::pre_cluster(unsigned n, unsigned N, const system<R_t, X_t, G_t>& S, const G_t::vertex& v, const R_t& r)
 
       This hook is run prior to cluster formation.
 
-      :param N_t n: The number of clusters already flipped.
-      :param N_t N: The total number of clusters to flip.
-      :param const system<R_t, X_t>& S: The current system state.
-      :param v_t i0: The index of the seed vertex for the cluster.
+      :param unsigned n: The number of clusters already flipped.
+      :param unsigned N: The total number of clusters to flip.
+      :param const system<R_t, X_t, G_t>& S: The current system state.
+      :param const G_t\:\:vertex& v: The seed vertex for the cluster.
       :param const R_t& r: The transformation by which the cluster will be flipped.
 
-   .. function:: void measurement::plain_bond_visited(const system<R_t, X_t>& S, v_t i1, const X_t& s1_new, v_t i2, double dE)
+   .. function:: void measurement::plain_bond_visited(const system<R_t, X_t, G_t>& S, const G_t::halfedge& e, const X_t& s1_new, double dE)
 
          This hook is run when an ordinary edge is visited during cluster formation, whether it is activated or not.
 
-         :param const system<R_t, X_t>& S: The current system state.
-         :param v_t i1:                    The index of the vertex soon to be flipped.
-         :param const X_t& s1_new:         The state vertex :code:`i1` will be flipped to.
-         :param v_t i2:                    The other vertex sharing the edge.
-         :param double dE:                 The change in energy that will result when the spin at site :code:`i1` is flipped.
+         :param const system<R_t, X_t, G_t>& S: The current system state.
+         :param const G_t\:\:halfedge& e:   The edge being considered, with e.self referencing the vertex soon to be flipped.
+         :param const X_t& s1_new:         The state the vertex will be flipped to.
+         :param double dE:                 The change in energy that will result when the spin is flipped.
 
-   .. cpp:function:: void measurement::plain_site_transformed(const system<R_t, X_t>& S, v_t i, const X_t& si_new)
+   .. cpp:function:: void measurement::plain_site_transformed(const system<R_t, X_t, G_t>& S, const G_t::vertex& v, const X_t& si_new)
 
          This hook is run when an ordinary site is about to be flipped.
 
-         :param const system<R_t, X_t>& S: The current system state.
-         :param v_t i:                     The index of the vertex soon to be flipped.
-         :param const X_t& si_new:          The state vertex :code:`i` will be flipped to.
+         :param const system<R_t, X_t, G_t>& S: The current system state.
+         :param const G_t\:\:vertex& v:                     The vertex soon to be flipped.
+         :param const X_t& si_new:          The state that vertex will be flipped to.
 
-   .. cpp:function:: void measurement::ghost_bond_visited(const system<R_t, X_t>& S, v_t i, const X_t& s0si_old, const X_t& s0si_new, double dE)
+   .. cpp:function:: void measurement::ghost_bond_visited(const system<R_t, X_t, G_t>& S, const G_t::vertex& v, const X_t& s0si_old, const X_t& s0si_new, double dE)
 
          This hook is run when an edge containing the ghost site is visited during cluster formation, whether activated or not.
 
-         :param const system<R_t, X_t>& S: The current system state.
-         :param v_t i:                     The index of the non-ghost site in this edge.
+         :param const system<R_t, X_t, G_t>& S: The current system state.
+         :param const G_t\:\:vertex& v:                     The non-ghost site in this edge.
          :param const X_t& s0si_old:       The state that the spin on the non-ghost site has before the transformation is applied, rotated by the inverse action of the ghost site.
          :param const X_t& s0si_new:       The state that the spin on the non-ghost site will have after the transformation is applied, rotated by the inverse action of the ghost site.
          :param double dE:                 The change in energy that will result when one site is transformed.
 
-   .. cpp:function:: void measurement::ghost_site_transformed(const system<R_t, X_t>& S, const R_t& R_new)
+   .. cpp:function:: void measurement::ghost_site_transformed(const system<R_t, X_t, G_t>& S, const R_t& R_new)
 
          This hook is run when the ghost site is about to be flipped.
 
-         :param const system<R_t, X_t>& S: The current system state.
+         :param const system<R_t, X_t, G_t>& S: The current system state.
          :param const R_t& R_new:          The state the ghost site will be flipped to.
 
-   .. cpp:function:: void measurement::post_cluster(N_t n, N_t N, const system<R_t, X_t>& S) 
+   .. cpp:function:: void measurement::post_cluster(unsigned n, unsigned N, const system<R_t, X_t, G_t>& S) 
 
          This hook is run after a cluster has been flipped.
 
-         :param N_t n:                     The number of clusters already flipped.
-         :param N_t N:                     The total number of clusters to flip.
-         :param const system<R_t, X_t>& S: The current system state.
+         :param unsigned n:                     The number of clusters already flipped.
+         :param unsigned N:                     The total number of clusters to flip.
+         :param const system<R_t, X_t, G_t>& S: The current system state.
 
diff --git a/doc/models.rst b/doc/models.rst
index e376e2e..8ba7e51 100644
--- a/doc/models.rst
+++ b/doc/models.rst
@@ -3,7 +3,7 @@
 Defining a Model
 ****************
 
-To define your own model, you need a class for spin states, a class for spin transformations, a bond and site coupling, and a generator of transformations. Many examples of models can be found in the headers contained in :file:`wolff/models/`.
+To define your own model, you need a class for spin states, a class for spin transformations, a bond and site coupling, and a generator of transformations. Many examples of models can be found in the headers contained in :file:`lib/wolff_models/`.
  
 Spin States
 ===========
diff --git a/doc/system.rst b/doc/system.rst
index 928d85b..9cd4a90 100644
--- a/doc/system.rst
+++ b/doc/system.rst
@@ -5,25 +5,25 @@
 Constructing a System & Running Wolff
 *************************************
 
-This class and associated member functions are defined in the header files :file:`wolff.hpp`, :file:`wolff/system.hpp`, and :file:`wolff/cluster.hpp`.
+This class and associated member functions are defined in the header file :file:`lib/wolff.hpp`.
 
-.. class:: template\<R_t, X_t> system
+.. class:: template\<R_t, X_t, G_t = graph\<>> system
 
    The core of the library lies in the :class:`system` class and its member functions. Here, the state of your model is stored and cluster flip Monte Carlo can be carried out in various ways.
 
    Note that the member objects and functions described here will change when compiled with certain compiler flags active, as described in :ref:`compile`.
 
-   .. member:: v_t system::nv
+   .. member:: unsigned system::nv
 
       Stores the number of ordinary sites in the model.
 
-   .. member:: v_t system::ne
+   .. member:: unsigned system::ne
 
       Stores the number of ordinary bonds in the model.
 
-   .. member:: graph system::G
+   .. member:: G_t system::G
 
-      Stores the graph describing the lattice, including the ghost site.
+      Stores the graph describing the lattice, including the ghost site. See :ref:`graphs` for details on this class. The template parameter has a default option corresponding to a graph type with no special vertex or edge properties.
 
    .. member:: double system::T
 
@@ -45,33 +45,33 @@ This class and associated member functions are defined in the header files :file
 
       The function that returns the coupling to the field.
 
-   .. function:: system::system(graph G, double T, std::function <double(const X_t&, const X_t&)> Z, std::function <double(const X_t&)> B)
+   .. function:: system::system(G_t G, double T, std::function <double(const X_t&, const X_t&)> Z, std::function <double(const X_t&)> B)
 
          The constructor for systems.
          
-         :param graph G: A lattice graph *without* the ghost spin added.
+         :param G_t G: A lattice graph *without* the ghost spin added.
          :param double T: The temperature.
          :param std\:\:function<double(const X_t&, const X_t&)> Z: The bond coupling.
          :param std\:\:function<double(const X_t&)> B: The field coupling.
          
          The states of the spins and ghost site are initialized using the default constructors for :type:`X_t` and :type:`R_t`, respectively. :any:`nv` and :any:`ne` are taken directly from :any:`G`, after which the ghost site is added to :any:`G`.
 
-   .. function:: system::flip_cluster(v_t i0, const R_t& r, std::mt19937& rng, measurement<R_t, X_t>& A)
+   .. function:: system::flip_cluster(const G_t::vertex& v0, const R_t& r, std::mt19937& rng, measurement<R_t, X_t, G_t>& A)
 
          Performs one Wolff cluster flip to the system. 
          
-         :param v_t i0: The index of the seed site.
+         :param const G_t\:\:vertex& v0: The vertex of the seed site.
          :param const R_t& r: The transformation by which the cluster is flipped.
          :param std\:\:mt19937& rng: A random number generator.
-         :param measurement<R_t, X_t>& A: Object whose class inherits :class:`measurement` and provides relevant measurement hooks.
+         :param measurement<R_t, X_t, G_t>& A: Object whose class inherits :class:`measurement` and provides relevant measurement hooks.
 
-   .. function:: system::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)
+   .. function:: system::run_wolff(unsigned N, std::function <R_t(std::mt19937&, const system<R_t, X_t, G_t>&, const G_t::vertex&)> r_gen, measurement<R_t, X_t, G_t>& A, std::mt19937& rng)
 
          Performs :any:`N` Wolff cluster flips to the system.
          
-         :param N_t N: Number of clusters to flip.
-         :param std\:\:function <R_t(std\:\:mt19937&, const system<R_t, X_t>&, v_t>)> r_gen: Generator of transformations for the cluster flips.
-         :param measurement<R_t, X_t>& A: Object whose class inherits :class:`measurement` and provides relevant measurement hooks.
+         :param unsigned N: Number of clusters to flip.
+         :param std\:\:function <R_t(std\:\:mt19937&, const system<R_t, X_t, G_t>&, const G_t\:\:vertex&>)> r_gen: Generator of transformations for the cluster flips.
+         :param measurement<R_t, X_t, G_t>& A: Object whose class inherits :class:`measurement` and provides relevant measurement hooks.
          :param std\:\:mt19937& rng: A random number generator.
 
 
diff --git a/doc/types.rst b/doc/types.rst
deleted file mode 100644
index b729ac3..0000000
--- a/doc/types.rst
+++ /dev/null
@@ -1,36 +0,0 @@
-
-************
-Custom Types
-************
-
-The Wolff library uses several custom types, defined to promote memory and cache effiency, promote shorter lines, and associate natural scope to certain parameters. They are all aliases for `standard C99 fixed width integer types`_. These are defined in the header file :file:`wolff/types.h`. Here is a list of the types and what they are used to hold:
-
-.. cpp:type:: uint_fast32_t v_t
-
-   Holds indicies for lattice vertices and edges.
-
-.. cpp:type:: uint_fast8_t q_t
-
-   Holds indicies for enumerating states, as in :math:`q`-state Potts.
-
-.. cpp:type:: uint_fast8_t D_t
-
-   Holds the dimension of space.
-
-.. cpp:type:: uint_fast16_t L_t
-
-   Holds the linear extent of the lattice.
-
-.. cpp:type:: uint_fast64_t N_t
-
-   Holds a count of Monte Carlo steps.
-
-Other Definitions
-=================
-
-All types are unsigned integers of a minimum size appropriate for most realistic purposes, with the :c:type:`fast` directives meaning that your compiler may choose a larger type for efficiency's sake.
-
-For each type :c:type:`x_t`, the macro :c:macro:`MAX_x` supplies the maximum value the type may hold, the macro :c:macro:`PRIx` supplies the format constants for use with the :c:func:`fprintf` family of functions, and the macro :c:macro:`SCNx` supplies the format constants for use with the :c:func:`fscanf` family of functions.
-
-.. _standard C99 fixed width integer types: https://en.cppreference.com/w/c/types/integer
-
diff --git a/examples/On.cpp b/examples/On.cpp
index 6c3f0fd..6885d2e 100644
--- a/examples/On.cpp
+++ b/examples/On.cpp
@@ -3,23 +3,21 @@
 #include <iostream>
 #include <chrono>
 
-#include "simple_measurement.hpp"
-
-#include <wolff/models/vector.hpp>
-#include <wolff/models/orthogonal.hpp>
+#include <wolff_models/vector.hpp>
+#include <wolff_models/orthogonal.hpp>
 
-#include <wolff.hpp>
+#include "simple_measurement.hpp"
 
 int main(int argc, char *argv[]) {
 
   // set defaults
-  N_t N = (N_t)1e4;
-  D_t D = 2;
-  L_t L = 128;
+  unsigned N = (unsigned)1e4;
+  unsigned D = 2;
+  unsigned L = 128;
   double T = 0.8;
   vector_t<WOLFF_N, double> H;
   H.fill(0.0);
-  q_t Hi = 0;
+  unsigned Hi = 0;
 
   int opt;
 
@@ -27,7 +25,7 @@ int main(int argc, char *argv[]) {
   while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) {
     switch (opt) {
     case 'N': // number of steps
-      N = (N_t)atof(optarg);
+      N = (unsigned)atof(optarg);
       break;
     case 'D': // dimension
       D = atoi(optarg);
@@ -58,12 +56,12 @@ int main(int argc, char *argv[]) {
   };
 
   // initialize the lattice
-  graph G(D, L);
+  graph<> G(D, L);
 
   // initialize the system
   system<orthogonal_t<WOLFF_N, double>, vector_t<WOLFF_N, double>> S(G, T, Z, B);
 
-  std::function <orthogonal_t<WOLFF_N, double>(std::mt19937&, const system<orthogonal_t<WOLFF_N, double>, vector_t<WOLFF_N, double>>&, v_t)> gen_R = generate_rotation_uniform<WOLFF_N>;
+  std::function <orthogonal_t<WOLFF_N, double>(std::mt19937&, const system<orthogonal_t<WOLFF_N, double>, vector_t<WOLFF_N, double>, graph<>>&, const graph<>::vertex)> gen_R = generate_rotation_uniform<WOLFF_N, graph<>>;
 
   // initailze the measurement object
   simple_measurement A(S);
diff --git a/examples/clock.cpp b/examples/clock.cpp
index 0cba8f6..3403f23 100644
--- a/examples/clock.cpp
+++ b/examples/clock.cpp
@@ -5,26 +5,24 @@
 
 #define WOLFF_USE_FINITE_STATES
 #define WOLFF_FINITE_STATES_N WOLFF_POTTSQ
-
-#include <wolff/models/potts.hpp>
-#include <wolff/models/dihedral.hpp>
+#include <wolff_models/potts.hpp>
+#include <wolff_models/dihedral.hpp>
+#include <wolff.hpp>
 
 #include "simple_measurement.hpp"
 
-#include <wolff.hpp>
-
 using namespace wolff;
 
 int main(int argc, char *argv[]) {
 
   // set defaults
-  N_t N = (N_t)1e4;
-  D_t D = 2;
-  L_t L = 128;
+  unsigned N = (unsigned)1e4;
+  unsigned D = 2;
+  unsigned L = 128;
   double T = 2.26918531421;
   vector_t<2, double> H;
   H.fill(0.0);
-  q_t Hi = 0;
+  unsigned Hi = 0;
 
   int opt;
 
@@ -32,7 +30,7 @@ int main(int argc, char *argv[]) {
   while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) {
     switch (opt) {
     case 'N': // number of steps
-      N = (N_t)atof(optarg);
+      N = (unsigned)atof(optarg);
       break;
     case 'D': // dimension
       D = atoi(optarg);
@@ -63,22 +61,22 @@ int main(int argc, char *argv[]) {
   };
 
   // initialize the lattice
-  graph G(D, L);
+  graph<> G(D, L);
 
   // initialize the system
-  system<dihedral_t<WOLFF_POTTSQ>, potts_t<WOLFF_POTTSQ>> S(G, T, Z, B);
+  system<dihedral_t<WOLFF_POTTSQ>, potts_t<WOLFF_POTTSQ>, graph<>> S(G, T, Z, B);
 
   // initialize the random number generator
   auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
   std::mt19937 rng{seed};
 
   // define function that generates self-inverse rotations
-  std::function <dihedral_t<WOLFF_POTTSQ>(std::mt19937&, const system<dihedral_t<WOLFF_POTTSQ>, potts_t<WOLFF_POTTSQ>>&, v_t)> gen_r = [] (std::mt19937& r, const system<dihedral_t<WOLFF_POTTSQ>, potts_t<WOLFF_POTTSQ>>& S, v_t i0) -> dihedral_t<WOLFF_POTTSQ> {
+  std::function <dihedral_t<WOLFF_POTTSQ>(std::mt19937&, const system<dihedral_t<WOLFF_POTTSQ>, potts_t<WOLFF_POTTSQ>, graph<>>&, const graph<>::vertex&)> gen_r = [] (std::mt19937& r, const system<dihedral_t<WOLFF_POTTSQ>, potts_t<WOLFF_POTTSQ>, graph<>>& S, const graph<>::vertex& v) -> dihedral_t<WOLFF_POTTSQ> {
     dihedral_t<WOLFF_POTTSQ> rot;
     rot.is_reflection = true;
-    std::uniform_int_distribution<q_t> dist(0, WOLFF_POTTSQ - 2);
-    q_t x = dist(r);
-    rot.x = (2 * S.s[i0].x + x + 1) % WOLFF_POTTSQ;
+    std::uniform_int_distribution<unsigned> dist(0, WOLFF_POTTSQ - 2);
+    unsigned x = dist(r);
+    rot.x = (2 * S.s[v.ind].x + x + 1) % WOLFF_POTTSQ;
 
     return rot;
   };
diff --git a/examples/continuous_gaussian.cpp b/examples/continuous_gaussian.cpp
index ef08247..eda22cb 100644
--- a/examples/continuous_gaussian.cpp
+++ b/examples/continuous_gaussian.cpp
@@ -3,19 +3,18 @@
 #include <iostream>
 #include <chrono>
 
-#include "simple_measurement.hpp"
-
-#include <wolff/models/height.hpp>
-#include <wolff/models/dihedral_inf.hpp>
-
+#include <wolff_models/height.hpp>
+#include <wolff_models/dihedral_inf.hpp>
 #include <wolff.hpp>
 
+#include "simple_measurement.hpp"
+
 int main(int argc, char *argv[]) {
 
   // set defaults
-  N_t N = (N_t)1e4;
-  D_t D = 2;
-  L_t L = 128;
+  unsigned N = (unsigned)1e4;
+  unsigned D = 2;
+  unsigned L = 128;
   double T = 0.8;
   double H = 0.0;
 
@@ -25,7 +24,7 @@ int main(int argc, char *argv[]) {
   while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) {
     switch (opt) {
     case 'N': // number of steps
-      N = (N_t)atof(optarg);
+      N = (unsigned)atof(optarg);
       break;
     case 'D': // dimension
       D = atoi(optarg);
@@ -55,25 +54,25 @@ int main(int argc, char *argv[]) {
   };
 
   // initialize the lattice
-  graph G(D, L);
+  graph<> G(D, L);
 
   // initialize the system
-  system<dihedral_inf_t<double>, height_t<double>> S(G, T, Z, B);
+  system<dihedral_inf_t<double>, height_t<double>, graph<>> S(G, T, Z, B);
 
   bool odd_run = false;
 
-  std::function <dihedral_inf_t<double>(std::mt19937&, const system<dihedral_inf_t<double>, height_t<double>>&, v_t)> gen_R_IH = [&](std::mt19937& r, const system<dihedral_inf_t<double>, height_t<double>>& S, v_t i0) -> dihedral_inf_t<double> {
+  std::function <dihedral_inf_t<double>(std::mt19937&, const system<dihedral_inf_t<double>, height_t<double>, graph<>>&, const graph<>::vertex&)> gen_R_IH = [&](std::mt19937& r, const system<dihedral_inf_t<double>, height_t<double>, graph<>>& S, const graph<>::vertex& v) -> dihedral_inf_t<double> {
     dihedral_inf_t<double> rot;
     rot.is_reflection = true;
 
     if (odd_run) {
-      std::uniform_int_distribution<v_t> dist(0, S.nv - 2);
-      v_t j = i0;
+      std::uniform_int_distribution<unsigned> dist(0, S.nv - 2);
+      unsigned j = v.ind;
 
       //while (S.s[j].x == S.s[i0].x) {
-        v_t tmp = dist(r);
+        unsigned tmp = dist(r);
 
-        if (tmp < i0) {
+        if (tmp < v.ind) {
           j = tmp;
         } else {
           j = tmp + 1;
@@ -83,7 +82,7 @@ int main(int argc, char *argv[]) {
       rot.x = 2 * S.s[j].x;
     } else {
       std::normal_distribution<double> dist(0.0,1.0);
-      rot.x = 2 * S.s[i0].x + dist(r);
+      rot.x = 2 * S.s[v.ind].x + dist(r);
     }
 
     odd_run = !odd_run;
diff --git a/examples/discrete_gaussian.cpp b/examples/discrete_gaussian.cpp
index 75ea0f9..a6d6ceb 100644
--- a/examples/discrete_gaussian.cpp
+++ b/examples/discrete_gaussian.cpp
@@ -3,19 +3,18 @@
 #include <iostream>
 #include <chrono>
 
-#include "simple_measurement.hpp"
-
-#include <wolff/models/height.hpp>
-#include <wolff/models/dihedral_inf.hpp>
-
+#include <wolff_models/height.hpp>
+#include <wolff_models/dihedral_inf.hpp>
 #include <wolff.hpp>
 
+#include "simple_measurement.hpp"
+
 int main(int argc, char *argv[]) {
 
   // set defaults
-  N_t N = (N_t)1e4;
-  D_t D = 2;
-  L_t L = 128;
+  unsigned N = (unsigned)1e4;
+  unsigned D = 2;
+  unsigned L = 128;
   double T = 0.8;
   double H = 0.0;
 
@@ -25,7 +24,7 @@ int main(int argc, char *argv[]) {
   while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) {
     switch (opt) {
     case 'N': // number of steps
-      N = (N_t)atof(optarg);
+      N = (unsigned)atof(optarg);
       break;
     case 'D': // dimension
       D = atoi(optarg);
@@ -45,35 +44,35 @@ int main(int argc, char *argv[]) {
   }
 
   // define the spin-spin coupling
-  std::function <double(const height_t<int64_t>&, const height_t<int64_t>&)> Z = [] (const height_t<int64_t>& s1, const height_t<int64_t>& s2) -> double {
+  std::function <double(const height_t<long long>&, const height_t<long long>&)> Z = [] (const height_t<long long>& s1, const height_t<long long>& s2) -> double {
     return - pow(s1.x - s2.x, 2);
   };
 
   // define the spin-field coupling
-  std::function <double(const height_t<int64_t>&)> B = [&] (const height_t<int64_t>& s) -> double {
+  std::function <double(const height_t<long long>&)> B = [&] (const height_t<long long>& s) -> double {
     return - H * pow(s.x, 2);
   };
 
   // initialize the lattice
-  graph G(D, L);
+  graph<> G(D, L);
 
   // initialize the system
-  system<dihedral_inf_t<int64_t>, height_t<int64_t>> S(G, T, Z, B);
+  system<dihedral_inf_t<long long>, height_t<long long>, graph<>> S(G, T, Z, B);
 
   bool odd_run = false;
 
-  std::function <dihedral_inf_t<int64_t>(std::mt19937&, const system<dihedral_inf_t<int64_t>, height_t<int64_t>>&, v_t)> gen_R_IH = [&](std::mt19937& r, const system<dihedral_inf_t<int64_t>, height_t<int64_t>>& S, v_t i0) -> dihedral_inf_t<int64_t> {
-    dihedral_inf_t<int64_t> rot;
+  std::function <dihedral_inf_t<long long>(std::mt19937&, const system<dihedral_inf_t<long long>, height_t<long long>, graph<>>&, const graph<>::vertex&)> gen_R_IH = [&](std::mt19937& r, const system<dihedral_inf_t<long long>, height_t<long long>, graph<>>& S, const graph<>::vertex &v) -> dihedral_inf_t<long long> {
+    dihedral_inf_t<long long> rot;
     rot.is_reflection = true;
 
     if (odd_run) {
-      std::uniform_int_distribution<v_t> dist(0, S.nv - 2);
-      v_t j = i0;
+      std::uniform_int_distribution<unsigned> dist(0, S.nv - 2);
+      unsigned j = v.ind;
 
       //while (S.s[j].x == S.s[i0].x) {
-        v_t tmp = dist(r);
+        unsigned tmp = dist(r);
 
-        if (tmp < i0) {
+        if (tmp < v.ind) {
           j = tmp;
         } else {
           j = tmp + 1;
@@ -85,9 +84,9 @@ int main(int argc, char *argv[]) {
       std::uniform_int_distribution<int> dist(0, 1);
       int j = dist(r);
       if (j) {
-        rot.x = 2 * S.s[i0].x + 1;
+        rot.x = 2 * S.s[v.ind].x + 1;
       } else {
-        rot.x = 2 * S.s[i0].x - 1;
+        rot.x = 2 * S.s[v.ind].x - 1;
       }
     }
 
diff --git a/examples/ising.cpp b/examples/ising.cpp
index c6e75fb..ecb296b 100644
--- a/examples/ising.cpp
+++ b/examples/ising.cpp
@@ -5,9 +5,7 @@
 
 #define WOLFF_USE_FINITE_STATES
 
-#include <wolff/models/ising.hpp>
-
-#include <wolff.hpp>
+#include <wolff_models/ising.hpp>
 
 #include "simple_measurement.hpp"
 
@@ -16,9 +14,9 @@ using namespace wolff;
 int main(int argc, char *argv[]) {
 
   // set defaults
-  N_t N = (N_t)1e4;
-  D_t D = 2;
-  L_t L = 128;
+  unsigned N = (unsigned)1e4;
+  unsigned D = 2;
+  unsigned L = 128;
   double T = 2.26918531421;
   double H = 0.0;
 
@@ -28,7 +26,7 @@ int main(int argc, char *argv[]) {
   while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) {
     switch (opt) {
     case 'N': // number of steps
-      N = (N_t)atof(optarg);
+      N = (unsigned)atof(optarg);
       break;
     case 'D': // dimension
       D = atoi(optarg);
@@ -58,17 +56,17 @@ int main(int argc, char *argv[]) {
   };
 
   // initialize the lattice
-  graph G(D, L);
+  graph<> G(D, L);
 
   // initialize the system
-  system<ising_t, ising_t> S(G, T, Z, B);
+  system<ising_t, ising_t, graph<>> S(G, T, Z, B);
 
   // initialize the random number generator
   auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
   std::mt19937 rng{seed};
 
   // define function that generates self-inverse rotations
-  std::function <ising_t(std::mt19937&, const system<ising_t, ising_t>&, v_t)> gen_r = gen_ising;
+  std::function <ising_t(std::mt19937&, const system<ising_t, ising_t, graph<>>&, const graph<>::vertex&)> gen_r = gen_ising<graph<>>;
 
   // initailze the measurement object
   simple_measurement A(S);
diff --git a/examples/ising_animation.cpp b/examples/ising_animation.cpp
index 7bafcf8..bcaa589 100644
--- a/examples/ising_animation.cpp
+++ b/examples/ising_animation.cpp
@@ -6,19 +6,18 @@
 #include <GL/glut.h>
 
 #define WOLFF_USE_FINITE_STATES
-
-#include <wolff/models/ising.hpp>
-
-#include <wolff.hpp>
+#include <wolff_models/ising.hpp>
 
 using namespace wolff;
 
-class draw_ising : public measurement<ising_t, ising_t> {
+typedef system<ising_t, ising_t, graph<>> sys;
+
+class draw_ising : public measurement<ising_t, ising_t, graph<>> {
   private:
     unsigned int frame_skip;
-    v_t C;
+    unsigned C;
   public:
-    draw_ising(const system<ising_t, ising_t>& S, unsigned int window_size, unsigned int frame_skip, int argc, char *argv[]) : frame_skip(frame_skip){
+    draw_ising(const sys& S, unsigned int window_size, unsigned int frame_skip, int argc, char *argv[]) : frame_skip(frame_skip){
       glutInit(&argc, argv);
       glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
       glutInitWindowSize(window_size, window_size);
@@ -29,9 +28,9 @@ class draw_ising : public measurement<ising_t, ising_t> {
       gluOrtho2D(0.0, S.G.L, 0.0, S.G.L);
     }
 
-    void pre_cluster(N_t, N_t, const system<ising_t, ising_t>& S, v_t, const ising_t&) {
+    void pre_cluster(unsigned, unsigned, const sys& S, const graph<>::vertex&, const ising_t&) override {
       glClear(GL_COLOR_BUFFER_BIT);
-      for (v_t i = 0; i < pow(S.G.L, 2); i++) {
+      for (unsigned i = 0; i < pow(S.G.L, 2); i++) {
         if (S.s[i].x == S.s0.x) {
           glColor3f(0.0, 0.0, 0.0);
         } else {
@@ -43,9 +42,9 @@ class draw_ising : public measurement<ising_t, ising_t> {
       C = 0;
     }
 
-    void plain_site_transformed(const system<ising_t, ising_t>& S, v_t i, const ising_t&) {
+    void plain_site_transformed(const sys& S, const graph<>::vertex& v, const ising_t&) override {
       glColor3f(1.0, 0.0, 0.0);
-      glRecti(i / S.G.L, i % S.G.L, (i / S.G.L) + 1, (i % S.G.L) + 1);
+      glRecti(v.ind / S.G.L, v.ind % S.G.L, (v.ind / S.G.L) + 1, (v.ind % S.G.L) + 1);
       C++;
       if (C % frame_skip == 0) {
         glFlush();
@@ -56,9 +55,9 @@ class draw_ising : public measurement<ising_t, ising_t> {
 int main(int argc, char *argv[]) {
 
   // set defaults
-  N_t N = (N_t)1e4;
-  D_t D = 2;
-  L_t L = 128;
+  unsigned N = (unsigned)1e4;
+  unsigned D = 2;
+  unsigned L = 128;
   double T = 2.26918531421;
   double H = 0.0;
   unsigned int window_size = 512;
@@ -69,29 +68,29 @@ int main(int argc, char *argv[]) {
   // take command line arguments
   while ((opt = getopt(argc, argv, "N:D:L:T:H:w:f:")) != -1) {
     switch (opt) {
-    case 'N': // number of steps
-      N = (N_t)atof(optarg);
-      break;
-    case 'D': // dimension
-      D = atoi(optarg);
-      break;
-    case 'L': // linear size
-      L = atoi(optarg);
-      break;
-    case 'T': // temperature
-      T = atof(optarg);
-      break;
-    case 'H': // external field
-      H = atof(optarg);
-      break;
-    case 'w':
-      window_size = atoi(optarg);
-      break;
-    case 'f':
-      frame_skip = atoi(optarg);
-      break;
-    default:
-      exit(EXIT_FAILURE);
+      case 'N': // number of steps
+        N = (unsigned)atof(optarg);
+        break;
+      case 'D': // dimension
+        D = atoi(optarg);
+        break;
+      case 'L': // linear size
+        L = atoi(optarg);
+        break;
+      case 'T': // temperature
+        T = atof(optarg);
+        break;
+      case 'H': // external field
+        H = atof(optarg);
+        break;
+      case 'w':
+        window_size = atoi(optarg);
+        break;
+      case 'f':
+        frame_skip = atoi(optarg);
+        break;
+      default:
+        exit(EXIT_FAILURE);
     }
   }
 
@@ -106,15 +105,10 @@ int main(int argc, char *argv[]) {
   };
 
   // initialize the lattice
-  graph G(D, L);
+  graph<> G(D, L);
 
   // initialize the system
-  system<ising_t, ising_t> S(G, T, Z, B);
-
-  // define function that generates self-inverse rotations
-  std::function <ising_t(std::mt19937&, const system<ising_t, ising_t>&, v_t)> gen_R = [] (std::mt19937&, const system<ising_t, ising_t>&, v_t) -> ising_t {
-    return ising_t(true);
-  };
+  sys S(G, T, Z, B);
 
   // initailze the measurement object
   draw_ising A(S, window_size, frame_skip, argc, argv);
@@ -124,7 +118,7 @@ int main(int argc, char *argv[]) {
   std::mt19937 rng{seed};
 
   // run wolff N times
-  S.run_wolff(N, gen_R, A, rng);
+  S.run_wolff(N, gen_ising<graph<>>, A, rng);
 
   // exit
   return 0;
diff --git a/examples/ising_no_field.cpp b/examples/ising_no_field.cpp
index 28d7c70..0a5b722 100644
--- a/examples/ising_no_field.cpp
+++ b/examples/ising_no_field.cpp
@@ -5,22 +5,18 @@
 
 #define WOLFF_NO_FIELD
 #define WOLFF_USE_FINITE_STATES
-
-#include <wolff/models/ising.hpp>
-
-#include <wolff.hpp>
+#include <wolff_models/ising.hpp>
 
 #include "simple_measurement.hpp"
 
-
 using namespace wolff;
 
 int main(int argc, char *argv[]) {
 
   // set defaults
-  N_t N = (N_t)1e4;
-  D_t D = 2;
-  L_t L = 128;
+  unsigned N = (unsigned)1e4;
+  unsigned D = 2;
+  unsigned L = 128;
   double T = 2.26918531421;
 
   int opt;
@@ -29,7 +25,7 @@ int main(int argc, char *argv[]) {
   while ((opt = getopt(argc, argv, "N:D:L:T:")) != -1) {
     switch (opt) {
     case 'N': // number of steps
-      N = (N_t)atof(optarg);
+      N = (unsigned)atof(optarg);
       break;
     case 'D': // dimension
       D = atoi(optarg);
@@ -51,17 +47,17 @@ int main(int argc, char *argv[]) {
   };
 
   // initialize the lattice
-  graph G(D, L);
+  graph<> G(D, L);
 
   // initialize the system
-  system<ising_t, ising_t> S(G, T, Z);
+  system<ising_t, ising_t, graph<>> S(G, T, Z);
 
   // initialize the random number generator
   auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
   std::mt19937 rng{seed};
 
   // define function that generates self-inverse rotations
-  std::function <ising_t(std::mt19937&, const system<ising_t, ising_t>&, v_t)> gen_r = gen_ising;
+  std::function <ising_t(std::mt19937&, const system<ising_t, ising_t, graph<>>&, const graph<>::vertex&)> gen_r = gen_ising<graph<>>;
 
   // initailze the measurement object
   simple_measurement A(S);
diff --git a/examples/ising_random_field.cpp b/examples/ising_random_field.cpp
index b37b0ce..9284797 100644
--- a/examples/ising_random_field.cpp
+++ b/examples/ising_random_field.cpp
@@ -4,21 +4,18 @@
 #include <chrono>
 
 #define WOLFF_SITE_DEPENDENCE
+#include <wolff_models/ising.hpp>
 
 #include "simple_measurement.hpp"
 
-#include <wolff/models/ising.hpp>
-
-#include <wolff.hpp>
-
 using namespace wolff;
 
 int main(int argc, char *argv[]) {
 
   // set defaults
-  N_t N = (N_t)1e4;
-  D_t D = 2;
-  L_t L = 128;
+  unsigned N = (unsigned)1e4;
+  unsigned D = 2;
+  unsigned L = 128;
   double T = 2.26918531421;
   double H = 0.0;
 
@@ -28,7 +25,7 @@ int main(int argc, char *argv[]) {
   while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) {
     switch (opt) {
     case 'N': // number of steps
-      N = (N_t)atof(optarg);
+      N = (unsigned)atof(optarg);
       break;
     case 'D': // dimension
       D = atoi(optarg);
@@ -52,29 +49,29 @@ int main(int argc, char *argv[]) {
     return (double)(s1 * s2);
   };
 
-  // initialize the lattice
-  graph G(D, L);
+  // initialize the lattice. vertex_prop is set to double and will contain the
+  // value of the random field at that site
+  graph<double> G(D, L);
 
   // initialize the random number generator
   auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
   std::mt19937 rng{seed};
 
   // define the spin-field coupling
-  std::vector<double> H_vals(G.nv);
   std::normal_distribution<double> distribution(0.0, H);
-  for (v_t i = 0; i < G.nv; i++) {
-    H_vals[i] = distribution(rng);
+  for (auto &v : G.vertices) {
+    v.prop = distribution(rng); 
   }
 
-  std::function <double(v_t, const ising_t&)> B = [&] (v_t i, const ising_t& s) -> double {
-    return H_vals[i] * s;
+  std::function <double(const graph<double>::vertex&, const ising_t&)> B = [] (const graph<double>::vertex& v, const ising_t& s) -> double {
+    return v.prop * s;
   };
 
   // initialize the system
-  system<ising_t, ising_t> S(G, T, Z, B);
+  system<ising_t, ising_t, graph<double>> S(G, T, Z, B);
 
   // define function that generates self-inverse rotations
-  std::function <ising_t(std::mt19937&, const system<ising_t, ising_t>&, v_t)> gen_r = gen_ising;
+  std::function <ising_t(std::mt19937&, const system<ising_t, ising_t, graph<double>>&, const graph<double>::vertex&)> gen_r = gen_ising<graph<double>>;
 
   // initailze the measurement object
   simple_measurement A(S);
diff --git a/examples/ising_standalone.cpp b/examples/ising_standalone.cpp
index 40572fd..6863ba5 100644
--- a/examples/ising_standalone.cpp
+++ b/examples/ising_standalone.cpp
@@ -21,27 +21,30 @@ class ising_t {
     }
 };
 
+typedef graph<> G_t;
+typedef system<ising_t, ising_t> sys;
+
 class measure_clusters : public measurement<ising_t, ising_t> {
   private:
-    v_t C;
+    unsigned C;
 
   public:
     double Ctotal;
 
     measure_clusters() { Ctotal = 0; }
 
-    void pre_cluster(N_t, N_t, const system<ising_t, ising_t>&, v_t, const ising_t&) { C = 0; }
+    void pre_cluster(unsigned, unsigned, const sys&, const G_t::vertex&, const ising_t&) override { C = 0; }
 
-    void plain_site_transformed(const system<ising_t, ising_t>&, v_t, const ising_t&) { C++; }
+    void plain_site_transformed(const sys&, const G_t::vertex&, const ising_t&) override { C++; }
 
-    void post_cluster(N_t, N_t, const system<ising_t, ising_t>&) { Ctotal += C; }
+    void post_cluster(unsigned, unsigned, const sys&) override { Ctotal += C; }
 };
 
 int main(int argc, char *argv[]) {
   // set defaults
-  N_t N = (N_t)1e3;
-  D_t D = 2;
-  L_t L = 128;
+  unsigned N = (unsigned)1e3;
+  unsigned D = 2;
+  unsigned L = 128;
   double T = 2.26918531421;
   double H = 0.01;
 
@@ -58,14 +61,14 @@ int main(int argc, char *argv[]) {
     };
 
   // initialize the lattice
-  graph G(D, L);
+  G_t G(D, L);
 
   // initialize the system
-  system<ising_t, ising_t> S(G, T, Z, B);
+  sys S(G, T, Z, B);
 
   // define function that generates self-inverse rotations
-  std::function <ising_t(std::mt19937&, const system<ising_t, ising_t>&, v_t)> gen_R =
-    [] (std::mt19937&, const system<ising_t, ising_t>&, v_t) -> ising_t {
+  std::function <ising_t(std::mt19937&, const sys&, const G_t::vertex&)> gen_R =
+    [] (std::mt19937&, const sys&, const G_t::vertex&) -> ising_t {
       return ising_t(-1);
     };
 
diff --git a/examples/potts.cpp b/examples/potts.cpp
index 24be2b3..c15de8d 100644
--- a/examples/potts.cpp
+++ b/examples/potts.cpp
@@ -5,26 +5,24 @@
 
 #define WOLFF_USE_FINITE_STATES
 #define WOLFF_FINITE_STATES_N WOLFF_POTTSQ
-
-#include <wolff/models/potts.hpp>
-#include <wolff/models/symmetric.hpp>
+#include <wolff_models/potts.hpp>
+#include <wolff_models/symmetric.hpp>
+#include <wolff.hpp>
 
 #include "simple_measurement.hpp"
 
-#include <wolff.hpp>
-
 using namespace wolff;
 
 int main(int argc, char *argv[]) {
 
   // set defaults
-  N_t N = (N_t)1e4;
-  D_t D = 2;
-  L_t L = 128;
+  unsigned N = (unsigned)1e4;
+  unsigned D = 2;
+  unsigned L = 128;
   double T = 2.26918531421;
   vector_t<WOLFF_POTTSQ, double> H;
   H.fill(0.0);
-  q_t Hi = 0;
+  unsigned Hi = 0;
 
   int opt;
 
@@ -32,7 +30,7 @@ int main(int argc, char *argv[]) {
   while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) {
     switch (opt) {
     case 'N': // number of steps
-      N = (N_t)atof(optarg);
+      N = (unsigned)atof(optarg);
       break;
     case 'D': // dimension
       D = atoi(optarg);
@@ -67,30 +65,30 @@ int main(int argc, char *argv[]) {
   };
 
   // initialize the lattice
-  graph G(D, L);
+  graph<> G(D, L);
 
   // initialize the system
-  system<symmetric_t<WOLFF_POTTSQ>, potts_t<WOLFF_POTTSQ>> S(G, T, Z, B);
+  system<symmetric_t<WOLFF_POTTSQ>, potts_t<WOLFF_POTTSQ>, graph<>> S(G, T, Z, B);
 
   // initialize the random number generator
   auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
   std::mt19937 rng{seed};
 
   // define function that generates self-inverse rotations
-  std::function <symmetric_t<WOLFF_POTTSQ>(std::mt19937&, const system<symmetric_t<WOLFF_POTTSQ>, potts_t<WOLFF_POTTSQ>>&, v_t)> gen_r = [] (std::mt19937& r, const system<symmetric_t<WOLFF_POTTSQ>, potts_t<WOLFF_POTTSQ>>& S, v_t i0) -> symmetric_t<WOLFF_POTTSQ> {
+  std::function <symmetric_t<WOLFF_POTTSQ>(std::mt19937&, const system<symmetric_t<WOLFF_POTTSQ>, potts_t<WOLFF_POTTSQ>, graph<>>&, const graph<>::vertex&)> gen_r = [] (std::mt19937& r, const system<symmetric_t<WOLFF_POTTSQ>, potts_t<WOLFF_POTTSQ>, graph<>>& S, const graph<>::vertex& v) -> symmetric_t<WOLFF_POTTSQ> {
     symmetric_t<WOLFF_POTTSQ> rot;
 
-    std::uniform_int_distribution<q_t> dist(0, WOLFF_POTTSQ - 2);
-    q_t j = dist(r);
-    q_t swap_v;
-    if (j < S.s[i0].x) {
+    std::uniform_int_distribution<unsigned> dist(0, WOLFF_POTTSQ - 2);
+    unsigned j = dist(r);
+    unsigned swap_v;
+    if (j < S.s[v.ind].x) {
       swap_v = j;
     } else {
       swap_v = j + 1;
     }
 
-    rot[S.s[i0].x] = swap_v;
-    rot[swap_v] = S.s[i0].x;
+    rot[S.s[v.ind].x] = swap_v;
+    rot[swap_v] = S.s[v.ind].x;
 
     return rot;
   };
diff --git a/examples/simple_measurement.hpp b/examples/simple_measurement.hpp
index 217f4f0..140da3b 100644
--- a/examples/simple_measurement.hpp
+++ b/examples/simple_measurement.hpp
@@ -1,31 +1,31 @@
 
-#include <wolff/measurement.hpp>
+#include <wolff.hpp>
 
 using namespace wolff;
 
-template <class R_t, class X_t>
-class simple_measurement : public measurement<R_t, X_t> {
+template <class R_t, class X_t, class G_t>
+class simple_measurement : public measurement<R_t, X_t, G_t> {
   private:
-    N_t n;
+    unsigned n;
 
     double E;
     typename X_t::M_t M;
-    v_t C;
+    unsigned C;
 
     double totalE;
     typename X_t::F_t totalM;
     double totalC;
 
   public:
-    simple_measurement(const system<R_t, X_t>& S) {
+    simple_measurement(const system<R_t, X_t, G_t>& S) {
       n = 0;
       M = S.nv * S.s[0];
       E = 0;
 
 #ifdef WOLFF_BOND_DEPENDENCE
-      for (v_t i = 0; i < S.nv; i++) {
-        for (v_t j : S.G.adjacency[i]) {
-          E -= 0.5 * S.Z(i, S.s[i], j, S.s[j]);
+      for (unsigned i = 0; i < S.nv; i++) {
+        for (const typename G_t::halfedge& e : S.G.vertices[i].edges) {
+          E -= 0.5 * S.Z(e, S.s[i], S.s[j]);
         }
       }
 #else
@@ -34,8 +34,8 @@ class simple_measurement : public measurement<R_t, X_t> {
 
 #ifndef WOLFF_NO_FIELD
 #ifdef WOLFF_SITE_DEPENDENCE
-      for (v_t i = 0; i < S.nv; i++) {
-        E -= S.B(i, S.s[i]);
+      for (unsigned i = 0; i < S.nv; i++) {
+        E -= S.B(S.G.vertices[i], S.s[i]);
       }
 #else
       E -= S.nv * S.B(S.s[0]);
@@ -47,28 +47,30 @@ class simple_measurement : public measurement<R_t, X_t> {
       totalC = 0;
     }
 
-    void pre_cluster(N_t, N_t, const system<R_t, X_t>&, v_t, const R_t&) {
+    void pre_cluster(unsigned, unsigned, const system<R_t, X_t, G_t>&, const typename G_t::vertex&, const R_t&) override {
       C = 0;
     }
 
-    void plain_bond_visited(const system<R_t, X_t>&, v_t, const X_t&, v_t, double dE) {
+    void plain_bond_visited(const system<R_t, X_t, G_t>&, const typename G_t::halfedge&, const X_t&, double dE) override {
       E += dE;
     }
 
-    void ghost_bond_visited(const system<R_t, X_t>&, v_t, const X_t& s_old, const X_t& s_new, double dE) {
+#ifndef WOLFF_NO_FIELD
+    void ghost_bond_visited(const system<R_t, X_t, G_t>&, const typename G_t::vertex&, const X_t& s_old, const X_t& s_new, double dE) override {
       E += dE;
       M += s_new - s_old;
     }
+#endif
 
-    void plain_site_transformed(const system<R_t, X_t>& S, v_t i, const X_t& si_new) {
+    void plain_site_transformed(const system<R_t, X_t, G_t>& S, const typename G_t::vertex& v, const X_t& si_new) override {
       C++;
 
 #ifdef WOLFF_NO_FIELD
-      M += si_new - S.s[i];
+      M += si_new - S.s[v.ind];
 #endif
     }
 
-    void post_cluster(N_t, N_t, const system<R_t, X_t>&) {
+    void post_cluster(unsigned, unsigned, const system<R_t, X_t, G_t>&) override {
       totalE += E;
       totalM += M;
       totalC += C;
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
+