big library overhual, fixed all examples, added documentation with sphinx

27 files changed, 764 insertions, 159 deletions

diff --git a/.gitignore b/.gitignore
index bf5bfbe..e7bdd38 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,2 @@
 build/*
-wolfram_link/*.o
-wolfram_link/*.c
-wolfram_link/convexminorant
+doc/_build/*
diff --git a/doc/Makefile b/doc/Makefile
new file mode 100644
index 0000000..298ea9e
--- /dev/null
+++ b/doc/Makefile
@@ -0,0 +1,19 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS    =
+SPHINXBUILD   = sphinx-build
+SOURCEDIR     = .
+BUILDDIR      = _build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
\ No newline at end of file
diff --git a/doc/compile.rst b/doc/compile.rst
new file mode 100644
index 0000000..2a3b4a2
--- /dev/null
+++ b/doc/compile.rst
@@ -0,0 +1,51 @@
+
+.. _compile:
+
+***************
+Compile Options
+***************
+
+There are several features that are not active by default, but can be used by defining a compiler flag. These will modify the constructor of :cpp:class:`system`. These can be used in concert when applicable. Note that at the moment, there is no :ref:`finite_states` support for systems with bond or site dependence.
+
+Building Without A Field
+========================
+
+.. c:macro:: WOLFF_NO_FIELD
+
+When this flag is defined Wolff will not use a field. When a field isn't necessary, this will improve performance: no ghost site will be initialize and no time will be wasted checking the energy change with respect to a uncoupled ghost site. 
+
+When defined, the constructor for :cpp:class:`system` becomes
+
+.. cpp:namespace:: cflags
+
+.. cpp:function:: system(graph, double, std::function <double(const X_t&, const X_t&)>)
+
+The resulting :cpp:class:`system` object does not have member objects :cpp:member:`B` or :cpp:member:`s0`, and :cpp:member:`system::G` does not have a ghost site initialized.
+
+Building With Bond Dependence
+=============================
+
+.. c:macro:: WOLFF_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:function:: 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.
+
+
+Building With Site Dependence
+=============================
+
+.. c:macro:: WOLFF_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:function:: 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.
+
diff --git a/doc/conf.py b/doc/conf.py
new file mode 100644
index 0000000..112d4b7
--- /dev/null
+++ b/doc/conf.py
@@ -0,0 +1,177 @@
+# -*- coding: utf-8 -*-
+#
+# Configuration file for the Sphinx documentation builder.
+#
+# This file does only contain a selection of the most common options. For a
+# full list see the documentation:
+# http://www.sphinx-doc.org/en/master/config
+
+# -- Path setup --------------------------------------------------------------
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#
+# import os
+# import sys
+# sys.path.insert(0, os.path.abspath('.'))
+
+
+# -- Project information -----------------------------------------------------
+
+project = 'wolff'
+copyright = '2018, Jaron Kent-Dobias'
+author = 'Jaron Kent-Dobias'
+
+# The short X.Y version
+version = ''
+# The full version, including alpha/beta/rc tags
+release = ''
+
+
+# -- General configuration ---------------------------------------------------
+
+# If your documentation needs a minimal Sphinx version, state it here.
+#
+# needs_sphinx = '1.0'
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
+# ones.
+extensions = [
+    'sphinx.ext.mathjax',
+]
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# The suffix(es) of source filenames.
+# You can specify multiple suffix as a list of string:
+#
+# source_suffix = ['.rst', '.md']
+source_suffix = '.rst'
+
+# The master toctree document.
+master_doc = 'index'
+
+# The language for content autogenerated by Sphinx. Refer to documentation
+# for a list of supported languages.
+#
+# This is also used if you do content translation via gettext catalogs.
+# Usually you set "language" from the command line for these cases.
+language = None
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This pattern also affects html_static_path and html_extra_path.
+exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = None
+
+
+# -- Options for HTML output -------------------------------------------------
+
+# The theme to use for HTML and HTML Help pages.  See the documentation for
+# a list of builtin themes.
+#
+html_theme = 'alabaster'
+
+# Theme options are theme-specific and customize the look and feel of a theme
+# further.  For a list of options available for each theme, see the
+# documentation.
+#
+# html_theme_options = {}
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['_static']
+
+# Custom sidebar templates, must be a dictionary that maps document names
+# to template names.
+#
+# The default sidebars (for documents that don't match any pattern) are
+# defined by theme itself.  Builtin themes are using these templates by
+# default: ``['localtoc.html', 'relations.html', 'sourcelink.html',
+# 'searchbox.html']``.
+#
+# html_sidebars = {}
+
+
+# -- Options for HTMLHelp output ---------------------------------------------
+
+# Output file base name for HTML help builder.
+htmlhelp_basename = 'wolffdoc'
+
+
+# -- Options for LaTeX output ------------------------------------------------
+
+latex_elements = {
+    # The paper size ('letterpaper' or 'a4paper').
+    #
+    # 'papersize': 'letterpaper',
+
+    # The font size ('10pt', '11pt' or '12pt').
+    #
+    # 'pointsize': '10pt',
+
+    # Additional stuff for the LaTeX preamble.
+    #
+    # 'preamble': '',
+
+    # Latex figure (float) alignment
+    #
+    # 'figure_align': 'htbp',
+}
+
+# Grouping the document tree into LaTeX files. List of tuples
+# (source start file, target name, title,
+#  author, documentclass [howto, manual, or own class]).
+latex_documents = [
+    (master_doc, 'wolff.tex', 'wolff Documentation',
+     'Jaron Kent-Dobias', 'manual'),
+]
+
+
+# -- Options for manual page output ------------------------------------------
+
+# One entry per manual page. List of tuples
+# (source start file, name, description, authors, manual section).
+man_pages = [
+    (master_doc, 'wolff', 'wolff Documentation',
+     [author], 1)
+]
+
+
+# -- Options for Texinfo output ----------------------------------------------
+
+# Grouping the document tree into Texinfo files. List of tuples
+# (source start file, target name, title, author,
+#  dir menu entry, description, category)
+texinfo_documents = [
+    (master_doc, 'wolff', 'wolff Documentation',
+     author, 'wolff', 'One line description of project.',
+     'Miscellaneous'),
+]
+
+
+# -- Options for Epub output -------------------------------------------------
+
+# Bibliographic Dublin Core info.
+epub_title = project
+
+# The unique identifier of the text. This can be a ISBN number
+# or the project homepage.
+#
+# epub_identifier = ''
+
+# A unique identification for the text.
+#
+# epub_uid = ''
+
+# A list of files that should not be packed into the epub file.
+epub_exclude_files = ['search.html']
+
+
+# -- Extension configuration -------------------------------------------------
\ No newline at end of file
diff --git a/doc/finite_states.rst b/doc/finite_states.rst
new file mode 100644
index 0000000..5d1fce8
--- /dev/null
+++ b/doc/finite_states.rst
@@ -0,0 +1,24 @@
+
+.. _finite_states:
+
+*************
+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 combination pair spins. Once the appropriate things have been defined for your model, the header :file:`wolff/finite_states.hpp` can be invoked 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
+====================
+
+.. c:macro:: WOLFF_FINITE_STATES_N
+
+This macro must be defined and given a value equal to the number of states your model can take.
+
+.. cpp:var:: const X_t finite_states_possible[WOLFF_FINITE_STATES_N]
+
+You must supply a constant array which lists each of the possible states of your individual spins. 
+
+.. cpp:function:: q_t finite_states_enum(const X_t&)
+
+You must define a function which takes each state and returns its index in :cpp:var:`finite_states_possible`.
+
diff --git a/doc/graph.rst b/doc/graph.rst
new file mode 100644
index 0000000..e44a14f
--- /dev/null
+++ b/doc/graph.rst
@@ -0,0 +1,47 @@
+
+******
+Graphs
+******
+
+This class is defined in the header :file:`wolff/graph.hpp`.
+
+.. cpp:class:: graph
+
+   Lattices are described by objects of class :cpp:class:`graph`, a minimal implementation of graphs.
+
+.. cpp:member:: D_t graph::D
+
+   The dimension of the graph. This property is unused by the core library.
+
+.. cpp:member:: L_t graph::L
+
+   The linear size of the graph. This property is unused by the core library.
+
+.. cpp:member:: v_t graph::ne
+
+   The number of edges in the graph. This property is unused by the core library.
+
+.. cpp:member:: v_t graph::nv
+
+   The number of vertices in the graph.
+
+.. cpp:member:: std::vector<std::vector<v_t>> graph::adjacency
+
+   The adjacency list for the graph. The :math:`i\text{th}` element of :cpp:member:`adjacency` is a standard library vector containing the indices of all vertices adjacent to vertex :math:`i`.
+
+.. cpp:member:: std::vector<std::vector<double>> graph::coordinate
+
+   The coordinates of the graph vertices. The :math:`i\text{th}` element of :cpp:var:`coordinate` is a standard library vector of length :cpp:var:`D` containing the spatial coordinates of vertex :math:`i`. This property is unused by the core library.
+
+.. cpp:function:: graph::graph()
+
+   The default constructor. Initializes an empty graph, i.e., :cpp:var:`D`, :cpp:var:`L`, :cpp:var:`ne`, and :cpp:var:`nv` are all zero and :cpp:var:`adjacency` and :cpp:var:`coordinate` are uninitialized.
+
+.. cpp:function:: graph::graph(D_t D, L_t L)
+
+   Calling the constructor with :cpp:var:`D` and :cpp:var:`L` will initialize the graph of a :cpp:var:`D`-dimensional hypercubic lattice with :cpp:var:`L` vertices per side. This is the only nontrivial graph constructor supplied by the core library.
+
+.. cpp:function:: void graph::add_ghost()
+
+   Calling this function on a graph object will add a ghost site to the graph. Explicitly, a new vertex is added that is adjacent to every other vertex in the graph. This vertex will have the last index, which is equal to number of vertices in the original graph.
+
diff --git a/doc/index.rst b/doc/index.rst
new file mode 100644
index 0000000..4e58e46
--- /dev/null
+++ b/doc/index.rst
@@ -0,0 +1,15 @@
+
+Wolff Libraries
++++++++++++++++
+
+.. toctree ::
+   :maxdepth: 2
+   
+   models
+   system
+   measurement
+   graph
+   types
+   finite_states
+   compile
+
diff --git a/doc/measurement.rst b/doc/measurement.rst
new file mode 100644
index 0000000..e71010d
--- /dev/null
+++ b/doc/measurement.rst
@@ -0,0 +1,84 @@
+
+************
+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`.
+
+.. cpp:class:: template\<class R_t, class X_t> 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.
+
+.. cpp:function:: template<void measurement::pre_cluster(N_t n, N_t N, const system<R_t, X_t>& S, v_t i0, const R_t& r)
+
+   This hook is run prior to cluster formation. The arguments supplied are:
+
+   =========================== ===
+   :cpp:var:`n`                the number of clusters already flipped
+   :cpp:var:`N`                the total number of clusters to flip
+   :cpp:var:`S`                the current system state
+   :cpp:var:`i0`               the index of the seed vertex for the cluster
+   :cpp:var:`r`                the transformation by which the cluster will be flipped
+   =========================== ===
+
+.. cpp: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)
+
+   This hook is run when an ordinary edge is visited during cluster formation, whether it is activated or not. The arguments supplied are:
+
+   ================================= ===
+   :code:`const system<R_t, X_t>& S` the current system state
+   :code:`v_t i1`                    the index of the vertex soon to be flipped
+   :code:`const X_t& s1_new`         the state vertex :cpp:var:`i1` will be flipped to
+   :code:`v_t i2`                    the other vertex sharing the edge
+   :code:`double dE`                 the change in energy that will result when the spin at site :cpp:var:`i1` is flipped
+   ================================= ===
+
+.. cpp:function:: void measurement::plain_site_transformed(const system<R_t, X_t>& S, v_t i, const X_t& si_new)
+
+   This hook is run when an ordinary site is about to be flipped. The arguments supplied are:
+
+   ================================= ===
+   :code:`const system<R_t, X_t>& S` the current system state
+   :code:`v_t i`                     the index of the vertex soon to be flipped
+   :code:`const X_t& si_new`         the state vertex :cpp:var:`i` 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)
+
+   This hook is run when an edge containing the ghost site is visited during cluster formation, whether activated or not. The arguments supplied are:
+
+   ================================= ===
+   :code:`const system<R_t, X_t>& S` the current system state
+   :code:`v_t i`                     the index of the non-ghost site in this edge
+   :code:`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
+   :code:`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
+   :code:`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)
+
+   This hook is run when the ghost site is about to be flipped. The arguments supplied are:
+
+   ================================= ===
+   :code:`const system<R_t, X_t>& S` the current system state
+   :code:`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) 
+
+   This hook is run after a cluster has been flipped. The arguments supplied are:
+
+   ================================= ===
+   :code:`N_t n`                     the number of clusters already flipped
+   :code:`N_t N`                     the total number of clusters to flip
+   :code:`const system<R_t, X_t>& S` the current system state
+   ================================= ===
+
diff --git a/doc/models.rst b/doc/models.rst
new file mode 100644
index 0000000..aa80ee6
--- /dev/null
+++ b/doc/models.rst
@@ -0,0 +1,61 @@
+
+****************
+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.
+
+Spin States
+===========
+
+.. cpp:class:: X_t
+
+   The spin type, which we have been denoting with :cpp:class:`X_t`, only needs to have a default constructor. If your spins can take only finitely many values, consider following the instructions in :ref:`finite_states` to significantly speed the algorithm.
+
+Transformations
+===============
+
+.. cpp:class:: R_t
+
+   The transformation type must have a default constructor, and the following member functions:
+
+   .. cpp:function:: X_t act(const X_t&)
+
+      The action of the transformation on a spin, returning the transformed spin.
+
+   .. cpp:function:: R_t act(const R_t&)
+
+      The action of the transformation on another transformation, returning the transformed transformation.
+
+   .. cpp:function:: X_t act_inverse(const X_t&)
+
+      The action of the inverse transformation on a spin, returning the transformed spin.
+
+   .. cpp:function:: R_t act_inverse(const R_t&)
+
+      The action of the inverse transformation on another transformation, returning the transformed transformation.
+
+Couplings
+=========
+
+When a :cpp:class:`system` object is initialized it needs to be given a bond and field coupling, to resemble the Hamiltonian
+
+.. math::
+         \mathcal H = -\sum_{\langle ij\rangle}Z(s_i,s_j)-\sum_iB(s_i)
+
+Note that building with certain compile flags will change the form that these coupling functions must take, as outlined in :ref:`compile`.
+
+Bond Coupling
+-------------
+
+.. cpp:function:: double Z(const X_t&, const X_t&)
+
+   A function that takes two spins and outputs the coupling between them. For traditional spin models, this is typically something like a dot product.
+
+Field Coupling
+--------------
+
+.. cpp:function:: double B(const X_t&)
+
+   A function that takes one spin and outputs the coupling between it and an external field.
+
diff --git a/doc/system.rst b/doc/system.rst
new file mode 100644
index 0000000..1751e2b
--- /dev/null
+++ b/doc/system.rst
@@ -0,0 +1,61 @@
+
+*************************************
+Constructing a System & Running Wolff
+*************************************
+
+The core of the library lies in the :cpp: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`.
+
+.. cpp:class:: template\<class R_t, class X_t> system
+
+Member Objects
+==============
+
+.. cpp:member:: v_t system::nv
+
+   Stores the number of ordinary sites in the model.
+
+.. cpp:member:: v_t system::ne
+
+   Stores the number of ordinary bonds in the model.
+
+.. cpp:member:: graph system::G
+
+   Stores the graph describing the lattice, including the ghost site.
+
+.. cpp:member:: double system::T
+
+   Stores the temperature of the model.
+
+.. cpp:member:: std::vector<X_t> system::s
+
+   The :math:`i\text{th}` component stores the spin state of site :math:`i`.
+
+.. cpp:member:: R_t system::s0
+
+   Stores the state of the ghost site.
+
+.. cpp:member:: std::function <double(const X_t&, const X_t&)> system::Z
+
+   The function that returns the coupling between neighboring sites.
+
+.. cpp:member:: std::function <double(const X_t&)> system::B
+
+   The function that returns the coupling to the field.
+
+Member Functions
+================
+
+.. cpp:function:: system::system(graph G, double T, std::function <double(const X_t&, const X_t&)> Z, std::function <double(const X_t&)> B)
+
+   The constructor for systems. The arguments given are a graph :cpp:any:`G` *without* the ghost spin added, the temperature :cpp:any:`T`, and the coupling functions :cpp:any:`Z` and :cpp:any:`B`. The states of the spins and ghost site are initialized using the default constructors for :cpp:type:`X_t` and :cpp:type:`R_t`, respectively. :cpp:any:`nv` and :cpp:any:`ne` are taken directly from :cpp:any:`G`, after which the ghost site is added to :cpp:any:`G`.
+
+.. cpp:function:: system::flip_cluster(v_t i0, const R_t& r, std::mt19937& rng, measurement<R_t, X_t>& A)
+
+   Performs one Wolff cluster flip to the system. The cluster is seeded at vertex :cpp:any:`i0` and the spins added are transformed by :cpp:any:`r`. A random number generator :cpp:any:`rng` provides required random numbers during, and the relevant measurement hooks defined in the inherited class of :cpp:any:`A` are run during the cluster formation.
+
+.. cpp: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)
+
+   Performs :cpp:any:`N` Wolff cluster flips to the system. One must provide a function :cpp:func:`r_gen` that takes a random number generator, the system state, and the index of the seed site and returns a transformation for each flip. Also required are an object from a class which inherits :cpp:class:`measurement` to run measurements, and a random number generator :cpp:any:`rng`.
+
diff --git a/doc/types.rst b/doc/types.rst
new file mode 100644
index 0000000..b729ac3
--- /dev/null
+++ b/doc/types.rst
@@ -0,0 +1,36 @@
+
+************
+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/CMakeLists.txt b/examples/CMakeLists.txt
index 73922c6..a4a4a5c 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,11 +1,11 @@
 
 add_executable(ising ising.cpp)
-add_executable(animate_ising animate_ising.cpp)
+add_executable(ising_animation ising_animation.cpp)
 add_executable(xy On.cpp)
 
 add_compile_definitions(xy WOLFF_N=2)
 
 target_link_libraries(ising wolff)
-target_link_libraries(animate_ising wolff GL GLU glut)
+target_link_libraries(ising_animation wolff GL GLU glut)
 target_link_libraries(xy wolff)
 
diff --git a/examples/On.cpp b/examples/On.cpp
index e045f52..6c3f0fd 100644
--- a/examples/On.cpp
+++ b/examples/On.cpp
@@ -7,6 +7,7 @@
 
 #include <wolff/models/vector.hpp>
 #include <wolff/models/orthogonal.hpp>
+
 #include <wolff.hpp>
 
 int main(int argc, char *argv[]) {
@@ -56,10 +57,13 @@ int main(int argc, char *argv[]) {
     return H * s;
   };
 
+  // initialize the lattice
+  graph G(D, L);
+
   // initialize the system
-  wolff_system<orthogonal_t<WOLFF_N, double>, vector_t<WOLFF_N, double>> S(D, L, T, Z, B);
+  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 vector_t<WOLFF_N, double>&)> 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>>&, v_t)> gen_R = generate_rotation_uniform<WOLFF_N>;
 
   // initailze the measurement object
   simple_measurement A(S);
@@ -69,7 +73,7 @@ int main(int argc, char *argv[]) {
   std::mt19937 rng{seed};
 
   // run wolff N times
-  wolff<orthogonal_t<WOLFF_N, double>, vector_t<WOLFF_N, double>>(N, S, gen_R, A, rng);
+  S.run_wolff(N, gen_R, A, rng);
 
   // print the result of our measurements
   std::cout << "Wolff complete!\nThe average energy per site was " << A.avgE() / S.nv
diff --git a/examples/ising.cpp b/examples/ising.cpp
index ffcb7e4..fe21279 100644
--- a/examples/ising.cpp
+++ b/examples/ising.cpp
@@ -7,8 +7,11 @@
 
 #include <wolff/models/ising.hpp>
 #include <wolff/finite_states.hpp>
+
 #include <wolff.hpp>
 
+using namespace wolff;
+
 int main(int argc, char *argv[]) {
 
   // set defaults
@@ -53,23 +56,24 @@ int main(int argc, char *argv[]) {
     return H * s;
   };
 
+  // initialize the lattice
+  graph G(D, L);
+
   // initialize the system
-  wolff_system<ising_t, ising_t> S(D, L, T, Z, B);
+  system<ising_t, ising_t> 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 ising_t&)> gen_R = [] (std::mt19937&, const ising_t& s) -> ising_t {
-    return ising_t(true);
-  };
+  std::function <ising_t(std::mt19937&, const system<ising_t, ising_t>&, v_t)> gen_r = gen_ising;
 
   // initailze the measurement object
   simple_measurement A(S);
 
   // run wolff N times
-  wolff<ising_t, ising_t>(N, S, gen_R, A, rng);
+  S.run_wolff(N, gen_r, A, rng);
 
   // print the result of our measurements
   std::cout << "Wolff complete!\nThe average energy per site was " << A.avgE() / S.nv
diff --git a/examples/animate_ising.cpp b/examples/ising_animation.cpp
index e33736e..2104e91 100644
--- a/examples/animate_ising.cpp
+++ b/examples/ising_animation.cpp
@@ -7,14 +7,17 @@
 
 #include <wolff/models/ising.hpp>
 #include <wolff/finite_states.hpp>
+
 #include <wolff.hpp>
 
-class draw_ising : public wolff_measurement<ising_t, ising_t> {
+using namespace wolff;
+
+class draw_ising : public measurement<ising_t, ising_t> {
   private:
     unsigned int frame_skip;
     v_t C;
   public:
-    draw_ising(const wolff_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 system<ising_t, ising_t>& 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);
@@ -22,39 +25,39 @@ class draw_ising : public wolff_measurement<ising_t, ising_t> {
       glClearColor(0.0,0.0,0.0,0.0);
       glMatrixMode(GL_PROJECTION);
       glLoadIdentity();
-      gluOrtho2D(0.0, S.L, 0.0, S.L);
+      gluOrtho2D(0.0, S.G.L, 0.0, S.G.L);
     }
 
-    void pre_cluster(N_t, N_t, const wolff_system<ising_t, ising_t>& S, v_t, const ising_t&) {
+    void pre_cluster(N_t, N_t, const system<ising_t, ising_t>& S, v_t, const ising_t&) {
       glClear(GL_COLOR_BUFFER_BIT);
-      for (v_t i = 0; i < pow(S.L, 2); i++) {
+      for (v_t 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 {
           glColor3f(1.0, 1.0, 1.0);
         }
-        glRecti(i / S.L, i % S.L, (i / S.L) + 1, (i % S.L) + 1);
+        glRecti(i / S.G.L, i % S.G.L, (i / S.G.L) + 1, (i % S.G.L) + 1);
       }
       glFlush();
       C = 0;
     }
 
-    void plain_bond_visited(const wolff_system<ising_t, ising_t>&, v_t, const ising_t&, v_t, double dE) {}
+    void plain_bond_visited(const system<ising_t, ising_t>&, v_t, const ising_t&, v_t, double dE) {}
 
-    void ghost_bond_visited(const wolff_system<ising_t, ising_t>&, v_t, const ising_t& s_old, const ising_t& s_new, double dE) {}
+    void ghost_bond_visited(const system<ising_t, ising_t>&, v_t, const ising_t& s_old, const ising_t& s_new, double dE) {}
 
-    void plain_site_transformed(const wolff_system<ising_t, ising_t>& S, v_t i, const ising_t&) {
+    void plain_site_transformed(const system<ising_t, ising_t>& S, v_t i, const ising_t&) {
       glColor3f(1.0, 0.0, 0.0);
-      glRecti(i / S.L, i % S.L, (i / S.L) + 1, (i % S.L) + 1);
+      glRecti(i / S.G.L, i % S.G.L, (i / S.G.L) + 1, (i % S.G.L) + 1);
       C++;
       if (C % frame_skip == 0) {
         glFlush();
       }
     }
 
-    void ghost_site_transformed(const wolff_system<ising_t, ising_t>&, const ising_t&) {}
+    void ghost_site_transformed(const system<ising_t, ising_t>&, const ising_t&) {}
 
-    void post_cluster(N_t, N_t, const wolff_system<ising_t, ising_t>&) {}
+    void post_cluster(N_t, N_t, const system<ising_t, ising_t>&) {}
 };
 
 int main(int argc, char *argv[]) {
@@ -109,11 +112,14 @@ int main(int argc, char *argv[]) {
     return H * s;
   };
 
+  // initialize the lattice
+  graph G(D, L);
+
   // initialize the system
-  wolff_system<ising_t, ising_t> S(D, L, T, Z, B);
+  system<ising_t, ising_t> S(G, T, Z, B);
 
   // define function that generates self-inverse rotations
-  std::function <ising_t(std::mt19937&, const ising_t&)> gen_R = [] (std::mt19937&, const ising_t& s) -> ising_t {
+  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);
   };
 
@@ -125,7 +131,7 @@ int main(int argc, char *argv[]) {
   std::mt19937 rng{seed};
 
   // run wolff N times
-  wolff<ising_t, ising_t>(N, S, gen_R, A, rng);
+  S.run_wolff(N, gen_R, A, rng);
 
   // exit
   return 0;
diff --git a/examples/simple_measurement.hpp b/examples/simple_measurement.hpp
index 2287c58..518631c 100644
--- a/examples/simple_measurement.hpp
+++ b/examples/simple_measurement.hpp
@@ -1,8 +1,10 @@
 
 #include <wolff/measurement.hpp>
 
+using namespace wolff;
+
 template <class R_t, class X_t>
-class simple_measurement : public wolff_measurement<R_t, X_t> {
+class simple_measurement : public measurement<R_t, X_t> {
   private:
     N_t n;
 
@@ -15,7 +17,7 @@ class simple_measurement : public wolff_measurement<R_t, X_t> {
     double totalC;
 
   public:
-    simple_measurement(const wolff_system<R_t, X_t>& S) {
+    simple_measurement(const system<R_t, X_t>& S) {
       n = 0;
       M = S.nv * S.s[0];
       E = - (S.ne * S.Z(S.s[0], S.s[0]) + S.nv * S.B(S.s[0]));
@@ -25,26 +27,26 @@ class simple_measurement : public wolff_measurement<R_t, X_t> {
       totalC = 0;
     }
 
-    void pre_cluster(N_t, N_t, const wolff_system<R_t, X_t>&, v_t, const R_t&) {
+    void pre_cluster(N_t, N_t, const system<R_t, X_t>&, v_t, const R_t&) {
       C = 0;
     }
 
-    void plain_bond_visited(const wolff_system<R_t, X_t>&, v_t, const X_t&, v_t, double dE) {
+    void plain_bond_visited(const system<R_t, X_t>&, v_t, const X_t&, v_t, double dE) {
       E += dE;
     }
 
-    void ghost_bond_visited(const wolff_system<R_t, X_t>&, v_t, const X_t& s_old, const X_t& s_new, double 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) {
       E += dE;
       M += s_new - s_old;
     }
 
-    void plain_site_transformed(const wolff_system<R_t, X_t>&, v_t,  const X_t&) {
+    void plain_site_transformed(const system<R_t, X_t>&, v_t,  const X_t&) {
       C++;
     }
 
-    void ghost_site_transformed(const wolff_system<R_t, X_t>&, const R_t&) {}
+    void ghost_site_transformed(const system<R_t, X_t>&, const R_t&) {}
 
-    void post_cluster(N_t, N_t, const wolff_system<R_t, X_t>&) {
+    void post_cluster(N_t, N_t, const system<R_t, X_t>&) {
       totalE += E;
       totalM += M;
       totalC += C;
diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp
index b78cd4b..7e909c8 100644
--- a/lib/include/wolff.hpp
+++ b/lib/include/wolff.hpp
@@ -3,25 +3,30 @@
 #define WOLFF_H
 
 #include "wolff/cluster.hpp"
+#include "wolff/measurement.hpp"
+
+namespace wolff{
 
 template <class R_t, class X_t>
-void wolff(N_t N, wolff_system<R_t, X_t>& S,
-           std::function <R_t(std::mt19937&, X_t)> r_gen,
-           wolff_measurement<R_t, X_t>& A, std::mt19937& rng) {
+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, S.nv - 1);
+  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, S.s[i0]);
+    R_t r = r_gen(rng, *this, i0);
 
-    A.pre_cluster(n, N, S, i0, r);
+    A.pre_cluster(n, N, *this, i0, r);
 
-    wolff_cluster_flip<R_t, X_t>(S, i0, r, rng, A);
+    this->flip_cluster(i0, r, rng, A);
 
-    A.post_cluster(n, N, S);
+    A.post_cluster(n, N, *this);
   }
 }
 
+}
+
 #endif
 
diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp
index 055cdf3..b66f367 100644
--- a/lib/include/wolff/cluster.hpp
+++ b/lib/include/wolff/cluster.hpp
@@ -7,20 +7,19 @@
 #include <vector>
 #include <queue>
 
-#include "types.h"
 #include "system.hpp"
-#include "graph.hpp"
-#include "measurement.hpp"
+
+namespace wolff {
 
 template <class R_t, class X_t>
-void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r,
-                  std::mt19937& rng, wolff_measurement<R_t, X_t>& A) {
+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(S.G.nv, false);
+  std::vector<bool> marks(G.nv, false);
 
   while (!queue.empty()) {
     v_t i = queue.front();
@@ -33,21 +32,21 @@ void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r,
 #ifndef WOLFF_NO_FIELD
       R_t s0_new;
 
-      bool we_are_ghost = (i == S.nv);
+      bool we_are_ghost = (i == nv);
 
       if (we_are_ghost) {
-        s0_new = r.act(S.s0);
+        s0_new = r.act(s0);
       } else
 #endif
       {
-        si_new = r.act(S.s[i]);
+        si_new = r.act(s[i]);
       }
 
-      for (const v_t &j : S.G.adj[i]) {
+      for (const v_t &j : G.adjacency[i]) {
         double dE, p;
 
 #ifndef WOLFF_NO_FIELD
-        bool neighbor_is_ghost = (j == S.nv);
+        bool neighbor_is_ghost = (j == nv);
 
         if (we_are_ghost || neighbor_is_ghost) {
           X_t s0s_old, s0s_new;
@@ -55,18 +54,18 @@ void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r,
 
           if (neighbor_is_ghost) {
             non_ghost = i;
-            s0s_old = S.s0.act_inverse(S.s[i]);
-            s0s_new = S.s0.act_inverse(si_new);
+            s0s_old = s0.act_inverse(s[i]);
+            s0s_new = s0.act_inverse(si_new);
           } else {
             non_ghost = j;
-            s0s_old = S.s0.act_inverse(S.s[j]);
-            s0s_new = s0_new.act_inverse(S.s[j]);
+            s0s_old = s0.act_inverse(s[j]);
+            s0s_new = s0_new.act_inverse(s[j]);
           }
 
 #ifdef WOLFF_SITE_DEPENDENCE
-          dE = S.B(non_ghost, s0s_old) - S.B(non_ghost, s0s_new);
+          dE = B(non_ghost, s0s_old) - B(non_ghost, s0s_new);
 #else
-          dE = S.B(s0s_old) - S.B(s0s_new);
+          dE = B(s0s_old) - B(s0s_new);
 #endif
 
 #ifdef WOLFF_FINITE_STATES
@@ -75,28 +74,28 @@ void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r,
 #endif
 
           // run measurement hooks for encountering a ghost bond
-          A.ghost_bond_visited(S, non_ghost, s0s_old, s0s_new, dE);
+          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 = S.Z(i, S.s[i], j, S.s[j]) - S.Z(i, si_new, j, S.s[j]);
+          dE = Z(i, s[i], j, s[j]) - Z(i, si_new, j, s[j]);
 #else
-          dE = S.Z(S.s[i], S.s[j]) - S.Z(si_new, S.s[j]);
+          dE = Z(s[i], s[j]) - Z(si_new, s[j]);
 #endif
 
 #ifdef WOLFF_FINITE_STATES
-          p = finite_states_Zp[finite_states_enum(S.s[i])]
+          p = finite_states_Zp[finite_states_enum(s[i])]
                               [finite_states_enum(si_new)]
-                              [finite_states_enum(S.s[j])];
+                              [finite_states_enum(s[j])];
 #endif
 
           // run measurement hooks for encountering a plain bond
-          A.plain_bond_visited(S, i, si_new, j, dE);
+          A.plain_bond_visited(*this, i, si_new, j, dE);
         }
 
 #ifndef FINITE_STATES
-        p = 1.0 - exp(-dE / S.T);
+        p = 1.0 - exp(-dE / T);
 #endif
 
         if (dist(rng) < p) {
@@ -106,17 +105,19 @@ void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r,
 
 #ifndef WOLFF_NO_FIELD
       if (we_are_ghost) {
-        A.ghost_site_transformed(S, s0_new);
-        S.s0 = s0_new;
+        A.ghost_site_transformed(*this, s0_new);
+        s0 = s0_new;
       } else
 #endif
       {
-        A.plain_site_transformed(S, i, si_new);
-        S.s[i] = si_new;
+        A.plain_site_transformed(*this, i, si_new);
+        s[i] = si_new;
       }
     }
   }
 }
 
+}
+
 #endif
 
diff --git a/lib/include/wolff/finite_states.hpp b/lib/include/wolff/finite_states.hpp
index 4cca69e..7e54eff 100644
--- a/lib/include/wolff/finite_states.hpp
+++ b/lib/include/wolff/finite_states.hpp
@@ -8,13 +8,15 @@
 
 #include "system.hpp"
 
+namespace wolff {
+
 std::array<std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Zp;
 #ifndef WOLFF_NO_FIELD
 std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Bp;
 #endif
 
 template <class R_t, class X_t>
-void finite_states_init(const wolff_system<R_t, X_t>& S) {
+void finite_states_init(const system<R_t, X_t>& S) {
 #ifndef WOLFF_NO_FIELD
   for (q_t i = 0; i < WOLFF_FINITE_STATES_N; i++) {
     for (q_t j = 0; j < WOLFF_FINITE_STATES_N; j++) {
@@ -31,5 +33,7 @@ void finite_states_init(const wolff_system<R_t, X_t>& S) {
   }
 }
 
+}
+
 #endif
 
diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp
index fcd73fa..65b3941 100644
--- a/lib/include/wolff/graph.hpp
+++ b/lib/include/wolff/graph.hpp
@@ -5,23 +5,26 @@
 #include <cmath>
 #include <vector>
 
-#include "types.h"
+namespace wolff {
 
-typedef enum lattice_t {
-  SQUARE_LATTICE,
-  DIAGONAL_LATTICE
-} lattice_t;
+#include "types.h"
 
-class graph_t {
+class graph {
   public:
+    D_t D;
+    L_t L;
     v_t ne;
     v_t nv;
-    std::vector<std::vector<v_t>> adj;
+    std::vector<std::vector<v_t>> adjacency;
     std::vector<std::vector<double>> coordinate;
 
-    graph_t(D_t D, L_t L, lattice_t lat = SQUARE_LATTICE);
-    void add_ext();
+    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
index a6a5f79..604eaf5 100644
--- a/lib/include/wolff/measurement.hpp
+++ b/lib/include/wolff/measurement.hpp
@@ -4,21 +4,25 @@
 
 #include "system.hpp"
 
+namespace wolff {
+
 template <class R_t, class X_t>
-class wolff_measurement {
+class measurement {
   public:
-    virtual void pre_cluster(N_t, N_t, const wolff_system<R_t, X_t>&, v_t, const R_t&) = 0;
+    virtual void pre_cluster(N_t, N_t, const system<R_t, X_t>&, v_t, const R_t&) = 0;
 
-    virtual void plain_bond_visited(const wolff_system<R_t, X_t>&, v_t, const X_t&, v_t, double) = 0;
-    virtual void plain_site_transformed(const wolff_system<R_t, X_t>&, v_t, const X_t&) = 0;
+    virtual void plain_bond_visited(const system<R_t, X_t>&, v_t, const X_t&, v_t, double) = 0;
+    virtual void plain_site_transformed(const system<R_t, X_t>&, v_t, const X_t&) = 0;
 
 #ifndef WOLFF_NO_FIELD
-    virtual void ghost_bond_visited(const wolff_system<R_t, X_t>&, v_t, const X_t&, const X_t&, double) = 0;
-    virtual void ghost_site_transformed(const wolff_system<R_t, X_t>&, const R_t&) = 0;
+    virtual void ghost_bond_visited(const system<R_t, X_t>&, v_t, const X_t&, const X_t&, double) = 0;
+    virtual void ghost_site_transformed(const system<R_t, X_t>&, const R_t&) = 0;
 #endif
 
-    virtual void post_cluster(N_t, N_t, const wolff_system<R_t, X_t>&) = 0;
+    virtual void post_cluster(N_t, N_t, const system<R_t, X_t>&) = 0;
 };
 
+}
+
 #endif
 
diff --git a/lib/include/wolff/models/ising.hpp b/lib/include/wolff/models/ising.hpp
index c6a3a47..824c063 100644
--- a/lib/include/wolff/models/ising.hpp
+++ b/lib/include/wolff/models/ising.hpp
@@ -3,6 +3,9 @@
 #define WOLFF_MODELS_ISING
 
 #include "../types.h"
+#include "../system.hpp"
+
+namespace wolff {
 
 class ising_t {
   public:
@@ -73,9 +76,15 @@ 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);
+};
+
 #define WOLFF_FINITE_STATES_N 2
 const ising_t finite_states_possible[2] = {ising_t(0), ising_t(1)};
-q_t finite_states_enum(ising_t state) { return (q_t)state.x; }
+q_t finite_states_enum(const ising_t& state) { return (q_t)state.x; }
+
+}
 
 #endif
 
diff --git a/lib/include/wolff/models/orthogonal.hpp b/lib/include/wolff/models/orthogonal.hpp
index 9dd5ddd..58203b7 100644
--- a/lib/include/wolff/models/orthogonal.hpp
+++ b/lib/include/wolff/models/orthogonal.hpp
@@ -4,7 +4,8 @@
 #include <random>
 #include <cmath>
 
-#include <wolff/types.h>
+#include "../types.h"
+#include "../system.hpp"
 #include "vector.hpp"
 
 template <q_t q, class T>
@@ -113,7 +114,7 @@ class orthogonal_t : public std::array<std::array<T, q>, q> {
 };
 
 template <q_t q>
-orthogonal_t <q, double> generate_rotation_uniform (std::mt19937& r, const vector_t <q, double>& v) {
+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;
@@ -135,7 +136,7 @@ orthogonal_t <q, double> generate_rotation_uniform (std::mt19937& r, const vecto
 }
 
 template <q_t q>
-orthogonal_t <q, double> generate_rotation_perturbation (std::mt19937& r, const vector_t <q, double>& v0, double epsilon, unsigned int n) {
+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;
@@ -149,14 +150,14 @@ orthogonal_t <q, double> generate_rotation_perturbation (std::mt19937& r, const
     double cosr = cos(2 * M_PI * rotation / (double)n / 2.0);
     double sinr = sin(2 * M_PI * rotation / (double)n / 2.0);
 
-    v[0] = v0[0] * cosr - v0[1] * sinr;
-    v[1] = v0[1] * cosr + v0[0] * sinr;
+    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] = v0[i];
+      v[i] = S.s[i0][i];
     }
   } else {
-    v = v0;
+    v = S.s[i0];
   }
 
   double m_dot_v = 0;
diff --git a/lib/include/wolff/models/potts.hpp b/lib/include/wolff/models/potts.hpp
index 903f25f..10b440e 100644
--- a/lib/include/wolff/models/potts.hpp
+++ b/lib/include/wolff/models/potts.hpp
@@ -46,5 +46,5 @@ class potts_t {
 
 const potts_t<POTTSQ> states[256] = {{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, {11}, {12}, {13}, {14}, {15}, {16}, {17}, {18}, {19}, {20}, {21}, {22}, {23}, {24}, {25}, {26}, {27}, {28}, {29}, {30}, {31}, {32}, {33}, {34}, {35}, {36}, {37}, {38}, {39}, {40}, {41}, {42}, {43}, {44}, {45}, {46}, {47}, {48}, {49}, {50}, {51}, {52}, {53}, {54}, {55}, {56}, {57}, {58}, {59}, {60}, {61}, {62}, {63}, {64}, {65}, {66}, {67}, {68}, {69}, {70}, {71}, {72}, {73}, {74}, {75}, {76}, {77}, {78}, {79}, {80}, {81}, {82}, {83}, {84}, {85}, {86}, {87}, {88}, {89}, {90}, {91}, {92}, {93}, {94}, {95}, {96}, {97}, {98}, {99}, {100}, {101}, {102}, {103}, {104}, {105}, {106}, {107}, {108}, {109}, {110}, {111}, {112}, {113}, {114}, {115}, {116}, {117}, {118}, {119}, {120}, {121}, {122}, {123}, {124}, {125}, {126}, {127}, {128}, {129}, {130}, {131}, {132}, {133}, {134}, {135}, {136}, {137}, {138}, {139}, {140}, {141}, {142}, {143}, {144}, {145}, {146}, {147}, {148}, {149}, {150}, {151}, {152}, {153}, {154}, {155}, {156}, {157}, {158}, {159}, {160}, {161}, {162}, {163}, {164}, {165}, {166}, {167}, {168}, {169}, {170}, {171}, {172}, {173}, {174}, {175}, {176}, {177}, {178}, {179}, {180}, {181}, {182}, {183}, {184}, {185}, {186}, {187}, {188}, {189}, {190}, {191}, {192}, {193}, {194}, {195}, {196}, {197}, {198}, {199}, {200}, {201}, {202}, {203}, {204}, {205}, {206}, {207}, {208}, {209}, {210}, {211}, {212}, {213}, {214}, {215}, {216}, {217}, {218}, {219}, {220}, {221}, {222}, {223}, {224}, {225}, {226}, {227}, {228}, {229}, {230}, {231}, {232}, {233}, {234}, {235}, {236}, {237}, {238}, {239}, {240}, {241}, {242}, {243}, {244}, {245}, {246}, {247}, {248}, {249}, {250}, {251}, {252}, {253}, {254}, {255}};
 template <q_t q>
-q_t finite_states_enum(potts_t<q> state) { return (q_t)state.x; }
+q_t finite_states_enum(const potts_t<q>& state) { return (q_t)state.x; }
 
diff --git a/lib/include/wolff/system.hpp b/lib/include/wolff/system.hpp
index ac82f44..32ad38e 100644
--- a/lib/include/wolff/system.hpp
+++ b/lib/include/wolff/system.hpp
@@ -4,23 +4,25 @@
 
 #include <functional>
 #include <vector>
+#include <random>
 
-#include "types.h"
 #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 wolff_system {
+class system {
   public:
-    D_t D; // dimension
-    L_t L; // linear size
     v_t nv; // number of vertices
     v_t ne; // number of edges
-    graph_t G; // the graph defining the lattice with ghost
+    graph G; // the graph defining the lattice with ghost
     double T; // the temperature
     std::vector<X_t> s; // the state of the ordinary spins
-#ifndef WOLFF_NO_FIELD
-    R_t s0; // the current state of the ghost site
-#endif
 
 #ifdef WOLFF_BOND_DEPENDENCE
     std::function <double(v_t, const X_t&, v_t, const X_t&)> Z; // coupling between sites
@@ -29,6 +31,7 @@ class wolff_system {
 #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
@@ -36,7 +39,7 @@ class wolff_system {
 #endif
 #endif
 
-    wolff_system(D_t D, L_t L, double T,
+    system(graph g, double T,
 #ifdef WOLFF_BOND_DEPENDENCE
         std::function <double(v_t, const X_t&, v_t, const X_t&)> Z
 #else
@@ -49,7 +52,7 @@ class wolff_system {
         , std::function <double(const X_t&)> B
 #endif
 #endif
-        , lattice_t lat = SQUARE_LATTICE) : D(D), L(L), G(D, L, lat), T(T), Z(Z)
+        ) : G(g), T(T), Z(Z)
 #ifndef WOLFF_NO_FIELD
              , s0(), B(B)
 #endif
@@ -58,13 +61,18 @@ class wolff_system {
       ne = G.ne;
       s.resize(nv);
 #ifndef WOLFF_NO_FIELD
-      G.add_ext();
+      G.add_ghost();
 #endif
 #ifdef WOLFF_FINITE_STATES
       finite_states_init(*this);
 #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);
 };
 
+}
+
 #endif
 
diff --git a/lib/include/wolff/types.h b/lib/include/wolff/types.h
index 68bd09e..1bc8c4d 100644
--- a/lib/include/wolff/types.h
+++ b/lib/include/wolff/types.h
@@ -8,7 +8,7 @@ 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_q UINT_FAST8_MAX
 #define MAX_D UINT_FAST8_MAX
 #define MAX_L UINT_FAST16_MAX
 #define MAX_N UINT_FAST64_MAX
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index b42aa78..4f91be4 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -1,74 +1,55 @@
 
 #include <wolff/graph.hpp>
 
-graph_t::graph_t(D_t D, L_t L, lattice_t lat) {
-  switch (lat) {
-    case SQUARE_LATTICE: {
-        nv = pow(L, D);
-        ne = D * nv;
+namespace wolff {
 
-        adj.resize(nv);
-        coordinate.resize(nv);
-
-        for (std::vector<v_t> adj_i : adj) {
-          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;
-
-            adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1)));
-            adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1)));
-          }
-        }
-        break;
-      }
-    case DIAGONAL_LATTICE: {
-        nv = D * pow(L, D);
-        ne = 2 * nv;
+graph::graph() {
+  D = 0;
+  L = 0;
+  nv = 0;
+  ne = 0;
+}
 
-        adj.resize(nv);
-        coordinate.resize(nv);
+graph::graph(D_t D, L_t L) : D(D), L(L) {
+  nv = pow(L, D);
+  ne = D * nv;
 
-        for (std::vector<v_t> adj_i : adj) {
-          adj_i.reserve(4 * (D - 1));
-        }
+  adjacency.resize(nv);
+  coordinate.resize(nv);
 
-        for (D_t i = 0; i < D; i++) {
-          v_t sb = i * pow(L, D);
+  for (std::vector<v_t> adj_i : adjacency) {
+    adj_i.reserve(2 * D);
+  }
 
-          for (v_t j = 0; j < pow(L, D); j++) {
-            v_t vc = sb + j;
+  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;
 
-            adj[vc].push_back(((i + 1) % D) * pow(L, D) + j);
-            adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * (j / (v_t)pow(L, D - 1)) + (j + 1 - 2 * (i % 2)) % L);
-            adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j - i) % L);
-            adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j + 1 - i) % L);
-          }
-        }
-        break;
-      }
+      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_t::add_ext() {
-  for (std::vector<v_t>& adj_i : adj) {
+void graph::add_ghost() {
+  for (std::vector<v_t>& adj_i : adjacency) {
     adj_i.push_back(nv);
   }
 
-  adj.resize(nv + 1);
+  adjacency.resize(nv + 1);
   coordinate.resize(nv + 1);
-  adj[nv].reserve(nv);
+  adjacency[nv].reserve(nv);
 
   for (v_t i = 0; i < nv; i++) {
-    adj[nv].push_back(i);
+    adjacency[nv].push_back(i);
   }
 
   coordinate[nv].resize(coordinate[0].size());
 
   ne += nv;
-  nv += 1;
+  nv++;
+}
+
 }