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.

11 files changed, 179 insertions, 197 deletions

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;