refactor

11 files changed, 669 insertions, 673 deletions

diff --git a/euclidean.hpp b/euclidean.hpp
new file mode 100644
index 0000000..b1a5e80
--- /dev/null
+++ b/euclidean.hpp
@@ -0,0 +1,116 @@
+
+#pragma once
+
+#include "matrix.hpp"
+#include "spin.hpp"
+#include "vector.hpp"
+
+template <class T, int D> class Euclidean {
+public:
+  Vector<T, D> t;
+  Matrix<T, D> r;
+
+  Euclidean(T L) {
+    for (unsigned i = 0; i < D; i++) {
+      t(i) = 0;
+      r(i, i) = 1;
+      for (unsigned j = 1; j < D; j++) {
+        r(i, (i + j) % D) = 0;
+      }
+    }
+  }
+
+  Euclidean(Vector<T, D> t0, Matrix<T, D> r0) {
+    t = t0;
+    r = r0;
+  }
+
+  template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
+    Spin<T, D, S> s_new;
+
+    s_new.x = t + r * s.x;
+    s_new.s = s.s;
+
+    return s_new;
+  }
+
+  Euclidean act(const Euclidean& x) const {
+    Vector<T, D> tnew = r * x.t + t;
+    Matrix<T, D> rnew = r * x.r;
+
+    Euclidean pnew(tnew, rnew);
+
+    return pnew;
+  }
+
+  Euclidean inverse() const {
+    Vector<T, D> tnew = -r.transpose() * t;
+    Matrix<T, D> rnew = r.transpose();
+
+    Euclidean pnew(tnew, rnew);
+
+    return pnew;
+  }
+};
+
+template <class T, int D> class TorusGroup {
+private:
+  T L;
+
+public:
+  Vector<T, D> t;
+  Matrix<T, D> r;
+
+  /** brief TorusGroup - default constructor, constructs the identity
+   */
+  TorusGroup(T L) : L(L) {
+    for (unsigned i = 0; i < D; i++) {
+      t(i) = 0;
+      r(i, i) = 1;
+      for (unsigned j = 1; j < D; j++) {
+        r(i, (i + j) % D) = 0;
+      }
+    }
+  }
+
+  TorusGroup(T L, Vector<T, D> t0, Matrix<T, D> r0) : L(L) {
+    t = t0;
+    r = r0;
+  }
+
+  template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
+    Spin<T, D, S> s_new;
+
+    s_new.x = t + r * s.x;
+    s_new.s = s.s;
+
+    for (unsigned i = 0; i < D; i++) {
+      s_new.x(i) = fmod(L + s_new.x(i), L);
+    }
+
+    return s_new;
+  }
+
+  TorusGroup act(const TorusGroup& x) const {
+    Vector<T, D> tnew = r * x.t + t;
+    Matrix<T, D> rnew = r * x.r;
+
+    for (unsigned i = 0; i < D; i++) {
+      tnew(i) = fmod(L + tnew(i), L);
+    }
+
+    TorusGroup pnew(this->L, tnew, rnew);
+
+    return pnew;
+  }
+
+  TorusGroup inverse() const {
+    Vector<T, D> tnew = -r.transpose() * t;
+    Matrix<T, D> rnew = r.transpose();
+
+    TorusGroup pnew(this->L, tnew, rnew);
+
+    return pnew;
+  }
+};
+
diff --git a/matrix.hpp b/matrix.hpp
new file mode 100644
index 0000000..96b63fb
--- /dev/null
+++ b/matrix.hpp
@@ -0,0 +1,7 @@
+
+#pragma once
+
+#include <eigen3/Eigen/Dense>
+
+template <class T, int D> using Matrix = Eigen::Matrix<T, D, D>;
+
diff --git a/octree.hpp b/octree.hpp
new file mode 100644
index 0000000..1be6e0a
--- /dev/null
+++ b/octree.hpp
@@ -0,0 +1,87 @@
+
+#pragma once
+
+#include <array>
+#include <set>
+#include <unordered_map>
+
+#include "spin.hpp"
+#include "vector.hpp"
+
+namespace std {
+template <typename T, size_t N> struct hash<array<T, N>> {
+  typedef array<T, N> argument_type;
+  typedef size_t result_type;
+
+  result_type operator()(const argument_type& a) const {
+    hash<T> hasher;
+    result_type h = 0;
+    for (result_type i = 0; i < N; ++i) {
+      h = h * 31 + hasher(a[i]);
+    }
+    return h;
+  }
+};
+} // namespace std
+
+template <class T, int D, class S> class Octree {
+private:
+  unsigned N;
+  T L;
+  std::unordered_map<std::array<signed, D>, std::set<Spin<T, D, S>*>> data;
+
+public:
+  Octree(T L, unsigned N) {
+    L = L;
+    N = N;
+  }
+
+  std::array<signed, D> ind(Vector<T, D> x) const {
+    std::array<signed, D> ind;
+
+    for (unsigned i = 0; i < D; i++) {
+      ind[i] = (signed)std::floor(N * x(i) / L);
+    }
+
+    return ind;
+  }
+
+  void insert(Spin<T, D, S>* s) { data[ind(s->x)].insert(s); };
+
+  void remove(Spin<T, D, S>* s) {
+    data[ind(s->x)].erase(s);
+    if (data[ind(s->x)].empty()) {
+      data.erase(ind(s->x));
+    }
+  };
+
+  std::set<Spin<T, D, S>*> neighbors(const Vector<T, D>& x) const {
+    std::array<signed, D> i0 = ind(x);
+    std::set<Spin<T, D, S>*> ns;
+
+    nearest_neighbors_of(i0, D + 1, ns);
+
+    return ns;
+  };
+
+  void nearest_neighbors_of(std::array<signed, D> i0, unsigned depth,
+                            std::set<Spin<T, D, S>*>& ns) const {
+    if (depth == 0) {
+      auto it = data.find(i0);
+      if (it != data.end()) {
+        ns.insert(it->second.begin(), it->second.end());
+      }
+    } else {
+      for (signed j : {-1, 0, 1}) {
+        std::array<signed, D> i1 = i0;
+        if (N < 2) {
+          i1[depth - 1] += j;
+        } else {
+          i1[depth - 1] = (N + i1[depth - 1] + j) % N;
+        }
+        nearest_neighbors_of(i1, depth - 1, ns);
+      }
+    }
+  };
+};
+
diff --git a/quantity.hpp b/quantity.hpp
new file mode 100644
index 0000000..55d5a28
--- /dev/null
+++ b/quantity.hpp
@@ -0,0 +1,83 @@
+
+#pragma once
+
+#include <vector>
+#include <list>
+#include <cmath>
+#include <array>
+
+class Quantity {
+private:
+  double total;
+  double total2;
+  std::vector<double> C;
+  unsigned wait;
+
+public:
+  unsigned n;
+  std::list<double> hist;
+
+  Quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) {
+    n = 0;
+    total = 0;
+    total2 = 0;
+  }
+
+  void add(double x) {
+    if (n > wait) {
+      hist.push_front(x);
+      if (hist.size() > C.size()) {
+        hist.pop_back();
+        unsigned t = 0;
+        for (double a : hist) {
+          C[t] += a * x;
+          t++;
+        }
+        total += x;
+        total2 += pow(x, 2);
+      }
+    }
+    n++;
+  }
+
+  double avg() const { return total / (n - wait); }
+
+  double avg2() const { return total2 / (n - wait); }
+
+  std::vector<double> ρ() const {
+    double C0 = C.front() / (n - wait);
+    double avg2 = pow(total / (n - wait), 2);
+
+    std::vector<double> ρtmp;
+
+    for (double Ct : C) {
+      ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
+    }
+
+    return ρtmp;
+  }
+
+  std::array<double, 2> τ() const {
+    double τtmp = 0.5;
+    unsigned M = 1;
+    double c = 8.0;
+
+    std::vector<double> ρ_tmp = this->ρ();
+
+    while (c * τtmp > M && M < C.size()) {
+      τtmp += ρ_tmp[M];
+      M++;
+    }
+
+    return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)};
+  }
+
+  double σ() const {
+    return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2));
+  }
+
+  double serr() const { return sqrt(this->σ()); }
+
+  unsigned num_added() const { return n - wait; }
+};
+
diff --git a/smiley.hpp b/smiley.hpp
new file mode 100644
index 0000000..bd58fc1
--- /dev/null
+++ b/smiley.hpp
@@ -0,0 +1,23 @@
+
+#pragma once
+
+#include <array>
+
+const std::array<std::array<unsigned, 16>, 16> smiley = {
+    {{{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}},
+     {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}},
+     {{0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0}},
+     {{0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}},
+     {{1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1}},
+     {{1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1}},
+     {{1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1}},
+     {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}},
+     {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}},
+     {{1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1}},
+     {{1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1}},
+     {{1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1}},
+     {{0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0}},
+     {{0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0}},
+     {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}},
+     {{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}}}};
+
diff --git a/space_wolff.hpp b/space_wolff.hpp
index 7345726..196b3b5 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -1,488 +1,36 @@
 
-#include <cmath>
-#include <eigen3/Eigen/Dense>
-#include <fstream>
+#pragma once
+
+#include <array>
 #include <functional>
-#include <iostream>
 #include <list>
 #include <queue>
-#include <random>
-#include <set>
 #include <vector>
-#include <unordered_map>
-#include <array>
-#include <unordered_map>
-
-namespace std {
-template <typename T, size_t N> struct hash<array<T, N>> {
-  typedef array<T, N> argument_type;
-  typedef size_t result_type;
-
-  result_type operator()(const argument_type& a) const {
-    hash<T> hasher;
-    result_type h = 0;
-    for (result_type i = 0; i < N; ++i) {
-      h = h * 31 + hasher(a[i]);
-    }
-    return h;
-  }
-};
-} // namespace std
 
 #include "randutils/randutils.hpp"
 
-const std::array<std::array<unsigned, 16>, 16> smiley = {
-    {{{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}},
-     {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}},
-     {{0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0}},
-     {{0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}},
-     {{1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1}},
-     {{1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1}},
-     {{1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1}},
-     {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}},
-     {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}},
-     {{1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1}},
-     {{1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1}},
-     {{1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1}},
-     {{0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0}},
-     {{0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0}},
-     {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}},
-     {{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}}}};
-
-template <class T, int D> using Vector = Eigen::Matrix<T, D, 1>;
-
-template <class T, int D> using Matrix = Eigen::Matrix<T, D, D>;
-
-/** brief diff - periodic subtraction of vectors
- */
-template <class T, int D> Vector<T, D> diff(T L, Vector<T, D> v1, Vector<T, D> v2) {
-  Vector<T, D> v;
-
-  for (unsigned i = 0; i < D; i++) {
-    v(i) = std::abs(v1(i) - v2(i));
-    if (v(i) > L / 2) {
-      v(i) = L - v(i);
-    }
-  }
-
-  return v;
-}
-
-template <class T, int D, class S> class Spin {
-public:
-  Vector<T, D> x;
-  S s;
-};
-
-template <class T, int D> class Euclidean {
-public:
-  Vector<T, D> t;
-  Matrix<T, D> r;
-
-  Euclidean(T L) {
-    for (unsigned i = 0; i < D; i++) {
-      t(i) = 0;
-      r(i, i) = 1;
-      for (unsigned j = 1; j < D; j++) {
-        r(i, (i + j) % D) = 0;
-      }
-    }
-  }
-
-  Euclidean(Vector<T, D> t0, Matrix<T, D> r0) {
-    t = t0;
-    r = r0;
-  }
-
-  template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
-    Spin<T, D, S> s_new;
-
-    s_new.x = t + r * s.x;
-    s_new.s = s.s;
-
-    return s_new;
-  }
-
-  Euclidean act(const Euclidean& x) const {
-    Vector<T, D> tnew = r * x.t + t;
-    Matrix<T, D> rnew = r * x.r;
-
-    Euclidean pnew(tnew, rnew);
-
-    return pnew;
-  }
-
-  Euclidean inverse() const {
-    Vector<T, D> tnew = -r.transpose() * t;
-    Matrix<T, D> rnew = r.transpose();
-
-    Euclidean pnew(tnew, rnew);
-
-    return pnew;
-  }
-};
-
-template <class T, int D> class TorusGroup {
-private:
-  T L;
-
-public:
-  Vector<T, D> t;
-  Matrix<T, D> r;
-
-  /** brief TorusGroup - default constructor, constructs the identity
-   */
-  TorusGroup(T L) : L(L) {
-    for (unsigned i = 0; i < D; i++) {
-      t(i) = 0;
-      r(i, i) = 1;
-      for (unsigned j = 1; j < D; j++) {
-        r(i, (i + j) % D) = 0;
-      }
-    }
-  }
-
-  TorusGroup(T L, Vector<T, D> t0, Matrix<T, D> r0) : L(L) {
-    t = t0;
-    r = r0;
-  }
-
-  template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
-    Spin<T, D, S> s_new;
-
-    s_new.x = t + r * s.x;
-    s_new.s = s.s;
-
-    for (unsigned i = 0; i < D; i++) {
-      s_new.x(i) = fmod(L + s_new.x(i), L);
-    }
-
-    return s_new;
-  }
-
-  TorusGroup act(const TorusGroup& x) const {
-    Vector<T, D> tnew = r * x.t + t;
-    Matrix<T, D> rnew = r * x.r;
-
-    for (unsigned i = 0; i < D; i++) {
-      tnew(i) = fmod(L + tnew(i), L);
-    }
-
-    TorusGroup pnew(this->L, tnew, rnew);
-
-    return pnew;
-  }
-
-  TorusGroup inverse() const {
-    Vector<T, D> tnew = -r.transpose() * t;
-    Matrix<T, D> rnew = r.transpose();
-
-    TorusGroup pnew(this->L, tnew, rnew);
-
-    return pnew;
-  }
-};
-
-template <class T, int D, class S> class Octree {
-private:
-  unsigned N;
-  T L;
-  std::unordered_map<std::array<signed, D>, std::set<Spin<T, D, S>*>> data;
-
-public:
-  Octree(T L, unsigned N) {
-    L = L;
-    N = N;
-  }
-
-  std::array<signed, D> ind(Vector<T, D> x) const {
-    std::array<signed, D> ind;
-
-    for (unsigned i = 0; i < D; i++) {
-      ind[i] = (signed)std::floor(N * x(i) / L);
-    }
-
-    return ind;
-  }
-
-  void insert(Spin<T, D, S>* s) { data[ind(s->x)].insert(s); };
-
-  void remove(Spin<T, D, S>* s) {
-    data[ind(s->x)].erase(s);
-    if (data[ind(s->x)].empty()) {
-      data.erase(ind(s->x));
-    }
-  };
-
-  std::set<Spin<T, D, S>*> neighbors(const Vector<T, D>& x) const {
-    std::array<signed, D> i0 = ind(x);
-    std::set<Spin<T, D, S>*> ns;
-
-    nearest_neighbors_of(i0, D + 1, ns);
-
-    return ns;
-  };
-
-  void nearest_neighbors_of(std::array<signed, D> i0, unsigned depth, std::set<Spin<T, D, S>*>& ns) const {
-    if (depth == 0) {
-      auto it = data.find(i0);
-      if (it != data.end()) {
-        ns.insert(it->second.begin(), it->second.end());
-      }
-    } else {
-      for (signed j : {-1, 0, 1}) {
-        std::array<signed, D> i1 = i0;
-        if (N < 2) {
-          i1[depth - 1] += j;
-        } else {
-          i1[depth - 1] = (N + i1[depth - 1] + j) % N;
-        }
-        nearest_neighbors_of(i1, depth - 1, ns);
-      }
-    }
-  };
-};
-
-/*
-template <class T, int D, class S> class Octree {
-private:
-  T L;
-  unsigned N;
-  std::vector<std::set<Spin<T, D, S>*>> data;
-
-public:
-  Octree(T L, unsigned depth) : L(L), N(pow(2, depth)), data(pow(N, D)){};
-
-  unsigned ind(Vector<T, D> x) const {
-    unsigned pos_ind = 0;
-
-    for (unsigned i = 0; i < D; i++) {
-      pos_ind += pow(N, i) * (unsigned)std::floor(N * x(i) / L);
-    }
-
-    assert(pos_ind < data.size());
-
-    return pos_ind;
-  }
-
-  void insert(Spin<T, D, S>* s) { data[ind(s->x)].insert(s); };
-
-  void remove(Spin<T, D, S>* s) { data[ind(s->x)].erase(s); };
-
-  std::set<Spin<T, D, S>*> neighbors(const Vector<T, D>& x, unsigned depth) const {
-    std::set<Spin<T, D, S>*> ns;
-    std::set<unsigned> seen;
-    nearest_neighbors_of(ind(x), depth, ns, seen);
-    return ns;
-  };
-
-  void nearest_neighbors_of(unsigned ind, unsigned depth, std::set<Spin<T, D, S>*>& ns,
-                            std::set<unsigned>& seen) const {
-    ns.insert(data[ind].begin(), data[ind].end());
-    seen.insert(ind);
-
-    if (depth > 0) {
-      for (unsigned i = 0; i < D; i++) {
-        for (signed j : {-1, 1}) {
-          unsigned nind = pow(N, i + 1) * (ind / ((unsigned)pow(N, i + 1))) +
-                          fmod(pow(N, i + 1) + ind + j * pow(N, i), pow(N, i + 1));
-
-          if (!seen.contains(nind)) {
-            seen.insert(nind);
-            nearest_neighbors_of(nind, depth - 1, ns, seen);
-          }
-        }
-      }
-    }
-  }
-};
-*/
-
-class quantity {
-private:
-  double total;
-  double total2;
-  std::vector<double> C;
-  unsigned wait;
-
-public:
-  unsigned n;
-  std::list<double> hist;
-
-  quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) {
-    n = 0;
-    total = 0;
-    total2 = 0;
-  }
-
-  void add(double x) {
-    if (n > wait) {
-      hist.push_front(x);
-      if (hist.size() > C.size()) {
-        hist.pop_back();
-        unsigned t = 0;
-        for (double a : hist) {
-          C[t] += a * x;
-          t++;
-        }
-        total += x;
-        total2 += pow(x, 2);
-      }
-    }
-    n++;
-  }
-
-  double avg() const { return total / (n - wait); }
-
-  double avg2() const { return total2 / (n - wait); }
-
-  std::vector<double> ρ() const {
-    double C0 = C.front() / (n - wait);
-    double avg2 = pow(total / (n - wait), 2);
-
-    std::vector<double> ρtmp;
-
-    for (double Ct : C) {
-      ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
-    }
-
-    return ρtmp;
-  }
-
-  std::array<double, 2> τ() const {
-    double τtmp = 0.5;
-    unsigned M = 1;
-    double c = 8.0;
-
-    std::vector<double> ρ_tmp = this->ρ();
-
-    while (c * τtmp > M && M < C.size()) {
-      τtmp += ρ_tmp[M];
-      M++;
-    }
-
-    return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)};
-  }
-
-  double σ() const {
-    return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2));
-  }
-
-  double serr() const { return sqrt(this->σ()); }
-
-  unsigned num_added() const { return n - wait; }
-};
+#include "euclidean.hpp"
+#include "octree.hpp"
+#include "quantity.hpp"
+#include "spin.hpp"
 
 template <class U, int D, class R, class S> class Model;
 
 template <class U, int D, class R, class S> class measurement {
 public:
   virtual void pre_cluster(const Model<U, D, R, S>&, unsigned, const R&){};
-  virtual void plain_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>*, const Spin<U, D, S>*,
-                                  const Spin<U, D, S>&, double){};
+  virtual void plain_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>*,
+                                  const Spin<U, D, S>*, const Spin<U, D, S>&, double){};
   virtual void plain_site_transformed(const Model<U, D, R, S>&, const Spin<U, D, S>*,
                                       const Spin<U, D, S>&){};
 
-  virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&, const Spin<U, D, S>&,
-                                  double){};
+  virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&,
+                                  const Spin<U, D, S>&, double){};
   virtual void ghost_site_transformed(const Model<U, D, R, S>&, const R&){};
 
   virtual void post_cluster(const Model<U, D, R, S>&){};
 };
 
-template<int D>
-void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
-  if (level > 0) {
-    unsigned new_level = level - 1;
-    unsigned old_length = sequences.size();
-    for (std::array<double, D>& sequence : sequences) {
-      std::array<double, D> new_sequence = sequence;
-      new_sequence[new_level] = -1;
-      sequences.push_front(new_sequence);
-    }
-    one_sequences<D>(sequences, new_level);
-  }
-}
-
-template <class U, int D>
-std::vector<Vector<U, D>> torus_vecs(U L) {
-  std::vector<Vector<U, D>> vecs;
-  std::array<double, D> ini_sequence;
-  ini_sequence.fill(1);
-  std::list<std::array<double, D>> sequences;
-  sequences.push_back(ini_sequence);
-  one_sequences<D>(sequences, D);
-  sequences.pop_front(); // don't want the identity matrix!
- 
-  for (std::array<double, D> sequence : sequences) {
-    Vector<U, D> v;
-    for (unsigned i = 0; i < D; i++) {
-      if (sequence[i] == 1) {
-        v(i) = 0;
-       } else {
-        v(i) = L / 2;
-       }
-     }
-
-     vecs.push_back(v);
-   }
-
-  return vecs;
-}
-
-template <class U, int D>
-std::vector<Matrix<U, D>> torus_mats() {
-  std::vector<Matrix<U, D>> mats;
-
-  std::array<double, D> ini_sequence;
-  ini_sequence.fill(1);
-  std::list<std::array<double, D>> sequences;
-  sequences.push_back(ini_sequence);
-
-  one_sequences<D>(sequences, D);
-
-  sequences.pop_front(); // don't want the identity matrix!
-
-  for (std::array<double, D> sequence : sequences) {
-    Matrix<U, D> m;
-    for (unsigned i = 0; i < D; i++) {
-      for (unsigned j = 0; j < D; j++) {
-        if (i == j) {
-          m(i, j) = sequence[i];
-        } else {
-          m(i, j) = 0;
-        }
-      }
-    }
-
-    mats.push_back(m);
-  }
-
-  for (unsigned i = 0; i < D; i++) {
-    for (unsigned j = 0; j < D; j++) {
-      if (i != j) {
-        Matrix<U, D> m;
-        for (unsigned k = 0; k < D; k++) {
-          for (unsigned l = 0; l < D; l++) {
-            if ((k == i && l == j) || (k == j && l == i)) {
-              if (i < j) {
-                m(k, l) = 1;
-              } else {
-                m(k, l) = -1;
-              }
-            } else {
-              m(k, l) = 0;
-            }
-          }
-        }
-        mats.push_back(m);
-      }
-    }
-  }
-
-  return mats;
-}
-
 template <class U, int D, class R, class S> class Model {
 public:
   U L;
@@ -493,14 +41,11 @@ public:
   std::function<double(const Spin<U, D, S>&)> B;
   randutils::mt19937_rng rng;
 
-
   Model(U L, std::function<double(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z,
         std::function<double(const Spin<U, D, S>&)> B)
       : L(L), s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {}
 
   void step(double T, unsigned ind, R& r, measurement<U, D, R, S>& A) {
-    unsigned cluster_size = 0;
-
     std::queue<Spin<U, D, S>*> queue;
     queue.push(&(s[ind]));
 
@@ -556,17 +101,15 @@ public:
           dict.remove(si);
           *si = si_new;
           dict.insert(si);
-          cluster_size++;
         }
       }
     }
   }
 
-  void resize(double T, double P) {
-  }
+  void resize(double T, double P) {}
 
   void wolff(double T, std::function<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)> g,
-        measurement<U, D, R, S>& A, unsigned N) {
+             measurement<U, D, R, S>& A, unsigned N) {
     for (unsigned i = 0; i < N; i++) {
       R r = g(*this, rng);
       unsigned ind = rng.uniform((unsigned)0, (unsigned)(s.size() - 1));
@@ -579,3 +122,4 @@ public:
     }
   }
 };
+
diff --git a/spheres.cpp b/spheres.cpp
index 97c1e18..b7a9d94 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -1,100 +1,105 @@
 
-#include "space_wolff.hpp"
 #include <GL/glut.h>
+#include <fstream>
+#include <iostream>
+
+#include "space_wolff.hpp"
+#include "torus_symmetries.hpp"
 
 const unsigned D = 2;
 typedef Model<double, D, TorusGroup<double, D>, double> model;
 
 class animation : public measurement<double, 2, TorusGroup<double, D>, double> {
-  private:
-    uint64_t t1;
-    uint64_t t2;
-    unsigned n;
-    unsigned tmp;
-  public:
-    animation(double L, unsigned w, int argc, char *argv[]) {
-      t1 = 0;
-      t2 = 0;
-      n = 0;
-
-      glutInit(&argc, argv);
-      glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
-      glutInitWindowSize(w, w);
-      glutCreateWindow("wolff");
-      glClearColor(0.0,0.0,0.0,0.0);
-      glMatrixMode(GL_PROJECTION);
-      glLoadIdentity();
-      gluOrtho2D(-1, L + 1, -1 , L + 1);
-    }
+private:
+  uint64_t t1;
+  uint64_t t2;
+  unsigned n;
+  unsigned tmp;
+
+public:
+  animation(double L, unsigned w, int argc, char* argv[]) {
+    t1 = 0;
+    t2 = 0;
+    n = 0;
+
+    glutInit(&argc, argv);
+    glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
+    glutInitWindowSize(w, w);
+    glutCreateWindow("wolff");
+    glClearColor(0.0, 0.0, 0.0, 0.0);
+    glMatrixMode(GL_PROJECTION);
+    glLoadIdentity();
+    gluOrtho2D(-1, L + 1, -1, L + 1);
+  }
 
-    void pre_cluster(const model&, unsigned, const TorusGroup<double, D>&) override {
-      tmp = 0;
-    }
+  void pre_cluster(const model&, unsigned, const TorusGroup<double, D>&) override { tmp = 0; }
 
-    void plain_site_transformed(const model&, const Spin<double, D, double>*, const Spin<double, D, double>&) override {
-      tmp++;
-    }
+  void plain_site_transformed(const model&, const Spin<double, D, double>*,
+                              const Spin<double, D, double>&) override {
+    tmp++;
+  }
 
-    void post_cluster(const model& m) override {
-      glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
-      glClear(GL_COLOR_BUFFER_BIT);
-      for (const Spin<double, 2, double>& s : m.s) {
-        glBegin(GL_POLYGON);
-        unsigned n_points = 50;
-        glColor3f(0.0f, 0.0f, 0.0f);
-        for (unsigned i = 0; i < n_points; i++) {
-          glVertex2d(m.s0.inverse().act(s).x(0) + s.s * cos(2 * i * M_PI / n_points), m.s0.inverse().act(s).x(1) + s.s * sin(2 * i * M_PI / n_points));
-        }
-        glEnd();
+  void post_cluster(const model& m) override {
+    glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
+    glClear(GL_COLOR_BUFFER_BIT);
+    for (const Spin<double, 2, double>& s : m.s) {
+      glBegin(GL_POLYGON);
+      unsigned n_points = 50;
+      glColor3f(0.0f, 0.0f, 0.0f);
+      for (unsigned i = 0; i < n_points; i++) {
+        glVertex2d(m.s0.inverse().act(s).x(0) + s.s * cos(2 * i * M_PI / n_points),
+                   m.s0.inverse().act(s).x(1) + s.s * sin(2 * i * M_PI / n_points));
       }
-      glFlush();
-
-      t1 += tmp;
-      t2 += tmp * tmp;
-      n++;
+      glEnd();
     }
+    glFlush();
 
-    void clear() {
-      t1 = 0;
-      t2 = 0;
-      n = 0;
-    }
+    t1 += tmp;
+    t2 += tmp * tmp;
+    n++;
+  }
 
-    double var() {
-      return (t2 - t1 * t1 / (double)n) / (double)n;
-    }
+  void clear() {
+    t1 = 0;
+    t2 = 0;
+    n = 0;
+  }
+
+  double var() { return (t2 - t1 * t1 / (double)n) / (double)n; }
 };
 
-std::function<TorusGroup<double, D>(const model&, randutils::mt19937_rng&)> eGen(const std::vector<Matrix<double, 2>>& mats, const std::vector<Vector<double, 2>>& vecs, double ε, double L) {
-  return [&mats, &vecs, L, ε] (const model& M, randutils::mt19937_rng& rng) -> TorusGroup<double, 2> {
-    Vector<double, D> t;
-    Matrix<double, D> m;
-    unsigned flip = rng.uniform((unsigned)0, (unsigned)(mats.size() + vecs.size() - 1));
-    if (flip < mats.size()) {
-      unsigned f_ind = rng.uniform((unsigned)0, (unsigned)M.s.size());
-      t = M.s[f_ind].x;
-      for (unsigned j = 0; j < D; j++) {
-        t(j) = fmod(10 * L + t(j) + rng.variate<double, std::normal_distribution>(0.0, ε), L);
-      }
-      m = mats[flip];
-    } else {
-      for (unsigned j = 0; j < D; j++) {
-        for (unsigned k = 0; k < D; k++) {
-          if (j == k) {
-            m(j, k) = 1;
-          } else {
-            m(j, k) = 0;
+std::function<TorusGroup<double, D>(const model&, randutils::mt19937_rng&)>
+eGen(const std::vector<Matrix<double, 2>>& mats, const std::vector<Vector<double, 2>>& vecs,
+     double ε, double L) {
+  return
+      [&mats, &vecs, L, ε](const model& M, randutils::mt19937_rng& rng) -> TorusGroup<double, 2> {
+        Vector<double, D> t;
+        Matrix<double, D> m;
+        unsigned flip = rng.uniform((unsigned)0, (unsigned)(mats.size() + vecs.size() - 1));
+        if (flip < mats.size()) {
+          unsigned f_ind = rng.uniform((unsigned)0, (unsigned)M.s.size());
+          t = M.s[f_ind].x;
+          for (unsigned j = 0; j < D; j++) {
+            t(j) = fmod(10 * L + t(j) + rng.variate<double, std::normal_distribution>(0.0, ε), L);
+          }
+          m = mats[flip];
+        } else {
+          for (unsigned j = 0; j < D; j++) {
+            for (unsigned k = 0; k < D; k++) {
+              if (j == k) {
+                m(j, k) = 1;
+              } else {
+                m(j, k) = 0;
+              }
+            }
           }
-        }
-      }
-
-      t = vecs[flip - mats.size()];
-    }
 
-    TorusGroup<double, D> g(M.L, t, m);
-    return g;
-  };
+          t = vecs[flip - mats.size()];
+        }
 
+        TorusGroup<double, D> g(M.L, t, m);
+        return g;
+      };
 }
 
 int main(int argc, char* argv[]) {
@@ -110,23 +115,23 @@ int main(int argc, char* argv[]) {
 
   while ((opt = getopt(argc, argv, "n:N:L:T:H:")) != -1) {
     switch (opt) {
-      case 'n': 
-        n = (unsigned)atof(optarg);
-        break;
-      case 'N': 
-        N = (unsigned)atof(optarg);
-        break;
-      case 'L':
-        L = atof(optarg);
-        break;
-      case 'T':
-        T = atof(optarg);
-        break;
-      case 'H':
-        H = atof(optarg);
-        break;
-      default:
-        exit(1);
+    case 'n':
+      n = (unsigned)atof(optarg);
+      break;
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    case 'L':
+      L = atof(optarg);
+      break;
+    case 'T':
+      T = atof(optarg);
+      break;
+    case 'H':
+      H = atof(optarg);
+      break;
+    default:
+      exit(1);
     }
   }
 
@@ -134,25 +139,24 @@ int main(int argc, char* argv[]) {
   double a = 0.05;
 
   std::function<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z =
-    [L, a, k] (const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
-      Vector<double, D> d = diff(L, s1.x, s2.x);
-
-      double σ = s1.s + s2.s;
-      double δ = σ - sqrt(d.transpose() * d);
-
-      if (δ > - a * σ) {
-        return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2));
-      } else if (δ > - 2 * a * σ) {
-        return 0.5 * k * pow(δ + 2 * a * σ, 2);
-      } else {
-        return 0;
-      }
-    };
+      [L, a, k](const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
+    Vector<double, D> d = diff(L, s1.x, s2.x);
+
+    double σ = s1.s + s2.s;
+    double δ = σ - sqrt(d.transpose() * d);
 
-  std::function<double(Spin<double, D, double>)> B =
-    [L, H] (Spin<double, D, double> s) -> double {
-      return H * s.x(1);
-    };
+    if (δ > -a * σ) {
+      return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2));
+    } else if (δ > -2 * a * σ) {
+      return 0.5 * k * pow(δ + 2 * a * σ, 2);
+    } else {
+      return 0;
+    }
+  };
+
+  std::function<double(Spin<double, D, double>)> B = [L, H](Spin<double, D, double> s) -> double {
+    return H * s.x(1);
+  };
 
   std::vector<Matrix<double, D>> mats = torus_mats<double, D>();
   std::vector<Vector<double, D>> vecs = torus_vecs<double, D>(L);
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index f2f998a..dba7ab3 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -1,4 +1,7 @@
 
+#include <fstream>
+#include <iostream>
+
 #include "space_wolff.hpp"
 #include <GL/glut.h>
 
@@ -6,67 +9,66 @@ const unsigned D = 2;
 typedef Model<double, D, Euclidean<double, D>, double> model;
 
 class animation : public measurement<double, 2, Euclidean<double, D>, double> {
-  private:
-    uint64_t t1;
-    uint64_t t2;
-    unsigned n;
-    unsigned tmp;
-  public:
-    animation(double L, unsigned w, int argc, char *argv[]) {
-      t1 = 0;
-      t2 = 0;
-      n = 0;
-
-      glutInit(&argc, argv);
-      glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
-      glutInitWindowSize(w, w);
-      glutCreateWindow("wolff");
-      glClearColor(0.0,0.0,0.0,0.0);
-      glMatrixMode(GL_PROJECTION);
-      glLoadIdentity();
-      gluOrtho2D(-L, L, -L , L);
-    }
+private:
+  uint64_t t1;
+  uint64_t t2;
+  unsigned n;
+  unsigned tmp;
+
+public:
+  animation(double L, unsigned w, int argc, char* argv[]) {
+    t1 = 0;
+    t2 = 0;
+    n = 0;
+
+    glutInit(&argc, argv);
+    glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
+    glutInitWindowSize(w, w);
+    glutCreateWindow("wolff");
+    glClearColor(0.0, 0.0, 0.0, 0.0);
+    glMatrixMode(GL_PROJECTION);
+    glLoadIdentity();
+    gluOrtho2D(-L, L, -L, L);
+  }
 
-    void pre_cluster(const model&, unsigned, const Euclidean<double, D>&) override {
-      tmp = 0;
-    }
+  void pre_cluster(const model&, unsigned, const Euclidean<double, D>&) override { tmp = 0; }
 
-    void plain_site_transformed(const model&, const Spin<double, D, double>*, const Spin<double, D, double>&) override {
-      tmp++;
-    }
+  void plain_site_transformed(const model&, const Spin<double, D, double>*,
+                              const Spin<double, D, double>&) override {
+    tmp++;
+  }
 
-    void post_cluster(const model& m) override {
-      glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
-      glClear(GL_COLOR_BUFFER_BIT);
-      for (const Spin<double, 2, double>& s : m.s) {
-        glBegin(GL_POLYGON);
-        unsigned n_points = 50;
-        glColor3f(0.0f, 0.0f, 0.0f);
-        for (unsigned i = 0; i < n_points; i++) {
-          glVertex2d(m.s0.inverse().act(s).x(0) + s.s * cos(2 * i * M_PI / n_points), m.s0.inverse().act(s).x(1) + s.s * sin(2 * i * M_PI / n_points));
-        }
-        glEnd();
+  void post_cluster(const model& m) override {
+    glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
+    glClear(GL_COLOR_BUFFER_BIT);
+    for (const Spin<double, 2, double>& s : m.s) {
+      glBegin(GL_POLYGON);
+      unsigned n_points = 50;
+      glColor3f(0.0f, 0.0f, 0.0f);
+      for (unsigned i = 0; i < n_points; i++) {
+        glVertex2d(m.s0.inverse().act(s).x(0) + s.s * cos(2 * i * M_PI / n_points),
+                   m.s0.inverse().act(s).x(1) + s.s * sin(2 * i * M_PI / n_points));
       }
-      glFlush();
-
-      t1 += tmp;
-      t2 += tmp * tmp;
-      n++;
+      glEnd();
     }
+    glFlush();
 
-    void clear() {
-      t1 = 0;
-      t2 = 0;
-      n = 0;
-    }
+    t1 += tmp;
+    t2 += tmp * tmp;
+    n++;
+  }
 
-    double var() {
-      return (t2 - t1 * t1 / (double)n) / (double)n;
-    }
+  void clear() {
+    t1 = 0;
+    t2 = 0;
+    n = 0;
+  }
+
+  double var() { return (t2 - t1 * t1 / (double)n) / (double)n; }
 };
 
 std::function<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(double ε) {
-  return [ε] (const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
+  return [ε](const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
     Vector<double, D> t;
     Matrix<double, D> m;
 
@@ -85,7 +87,6 @@ std::function<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(
     Euclidean<double, D> g(t, m);
     return g;
   };
-
 }
 
 int main(int argc, char* argv[]) {
@@ -101,51 +102,50 @@ int main(int argc, char* argv[]) {
 
   while ((opt = getopt(argc, argv, "n:N:L:T:H:")) != -1) {
     switch (opt) {
-      case 'n': 
-        n = (unsigned)atof(optarg);
-        break;
-      case 'N': 
-        N = (unsigned)atof(optarg);
-        break;
-      case 'L':
-        L = atof(optarg);
-        break;
-      case 'T':
-        T = atof(optarg);
-        break;
-      case 'H':
-        H = atof(optarg);
-        break;
-      default:
-        exit(1);
+    case 'n':
+      n = (unsigned)atof(optarg);
+      break;
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    case 'L':
+      L = atof(optarg);
+      break;
+    case 'T':
+      T = atof(optarg);
+      break;
+    case 'H':
+      H = atof(optarg);
+      break;
+    default:
+      exit(1);
     }
   }
 
-  double k = 1e2;
-  double a = 0.05;
+  double k = 1e8;
+  double a = 0.0;
 
   std::function<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z =
-    [L, a, k] (const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
-      Vector<double, D> d = s1.x - s2.x;
-
-      double σ = s1.s + s2.s;
-      double δ = σ - sqrt(d.transpose() * d);
-
-      if (δ > - a * σ) {
-        return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2));
-      } else if (δ > - 2 * a * σ) {
-        return 0.5 * k * pow(δ + 2 * a * σ, 2);
-      } else {
-        return 0;
-      }
-    };
+      [L, a, k](const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
+    Vector<double, D> d = s1.x - s2.x;
+
+    double σ = s1.s + s2.s;
+    double δ = σ - sqrt(d.transpose() * d);
+
+    if (δ > -a * σ) {
+      return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2));
+    } else if (δ > -2 * a * σ) {
+      return 0.5 * k * pow(δ + 2 * a * σ, 2);
+    } else {
+      return 0;
+    }
+  };
 
-  std::function<double(Spin<double, D, double>)> B =
-    [L, H] (Spin<double, D, double> s) -> double {
-      return H * s.x.norm();
-    };
+  std::function<double(Spin<double, D, double>)> B = [L, H](Spin<double, D, double> s) -> double {
+    return H * s.x.norm();
+  };
 
-  auto g = eGen(1);
+  auto g = eGen(0.25);
   animation A(L, 750, argc, argv);
   model sphere(1, Z, B);
 
@@ -187,3 +187,4 @@ int main(int argc, char* argv[]) {
 
   return 0;
 }
+
diff --git a/spin.hpp b/spin.hpp
new file mode 100644
index 0000000..d4d39cc
--- /dev/null
+++ b/spin.hpp
@@ -0,0 +1,11 @@
+
+#pragma once
+
+#include "vector.hpp"
+
+template <class T, int D, class S> class Spin {
+public:
+  Vector<T, D> x;
+  S s;
+};
+
diff --git a/torus_symmetries.hpp b/torus_symmetries.hpp
new file mode 100644
index 0000000..17772de
--- /dev/null
+++ b/torus_symmetries.hpp
@@ -0,0 +1,100 @@
+
+#pragma once
+
+#include <array>
+#include <list>
+#include <vector>
+
+#include "matrix.hpp"
+#include "vector.hpp"
+
+template <int D> void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
+  if (level > 0) {
+    unsigned new_level = level - 1;
+    unsigned old_length = sequences.size();
+    for (std::array<double, D>& sequence : sequences) {
+      std::array<double, D> new_sequence = sequence;
+      new_sequence[new_level] = -1;
+      sequences.push_front(new_sequence);
+    }
+    one_sequences<D>(sequences, new_level);
+  }
+}
+
+template <class U, int D> std::vector<Vector<U, D>> torus_vecs(U L) {
+  std::vector<Vector<U, D>> vecs;
+  std::array<double, D> ini_sequence;
+  ini_sequence.fill(1);
+  std::list<std::array<double, D>> sequences;
+  sequences.push_back(ini_sequence);
+  one_sequences<D>(sequences, D);
+  sequences.pop_front(); // don't want the identity matrix!
+
+  for (std::array<double, D> sequence : sequences) {
+    Vector<U, D> v;
+    for (unsigned i = 0; i < D; i++) {
+      if (sequence[i] == 1) {
+        v(i) = 0;
+      } else {
+        v(i) = L / 2;
+      }
+    }
+
+    vecs.push_back(v);
+  }
+
+  return vecs;
+}
+
+template <class U, int D> std::vector<Matrix<U, D>> torus_mats() {
+  std::vector<Matrix<U, D>> mats;
+
+  std::array<double, D> ini_sequence;
+  ini_sequence.fill(1);
+  std::list<std::array<double, D>> sequences;
+  sequences.push_back(ini_sequence);
+
+  one_sequences<D>(sequences, D);
+
+  sequences.pop_front(); // don't want the identity matrix!
+
+  for (std::array<double, D> sequence : sequences) {
+    Matrix<U, D> m;
+    for (unsigned i = 0; i < D; i++) {
+      for (unsigned j = 0; j < D; j++) {
+        if (i == j) {
+          m(i, j) = sequence[i];
+        } else {
+          m(i, j) = 0;
+        }
+      }
+    }
+
+    mats.push_back(m);
+  }
+
+  for (unsigned i = 0; i < D; i++) {
+    for (unsigned j = 0; j < D; j++) {
+      if (i != j) {
+        Matrix<U, D> m;
+        for (unsigned k = 0; k < D; k++) {
+          for (unsigned l = 0; l < D; l++) {
+            if ((k == i && l == j) || (k == j && l == i)) {
+              if (i < j) {
+                m(k, l) = 1;
+              } else {
+                m(k, l) = -1;
+              }
+            } else {
+              m(k, l) = 0;
+            }
+          }
+        }
+        mats.push_back(m);
+      }
+    }
+  }
+
+  return mats;
+}
+
diff --git a/vector.hpp b/vector.hpp
new file mode 100644
index 0000000..2e87acd
--- /dev/null
+++ b/vector.hpp
@@ -0,0 +1,20 @@
+
+#pragma once
+
+#include <eigen3/Eigen/Dense>
+
+template <class T, int D> using Vector = Eigen::Matrix<T, D, 1>;
+
+template <class T, int D> Vector<T, D> diff(T L, Vector<T, D> v1, Vector<T, D> v2) {
+  Vector<T, D> v;
+
+  for (unsigned i = 0; i < D; i++) {
+    v(i) = std::abs(v1(i) - v2(i));
+    if (v(i) > L / 2) {
+      v(i) = L - v(i);
+    }
+  }
+
+  return v;
+}
+