reverted change to pZ

2 files changed, 337 insertions, 383 deletions

diff --git a/ising.cpp b/ising.cpp
index fffce19..1987a8f 100644
--- a/ising.cpp
+++ b/ising.cpp
@@ -54,17 +54,13 @@ int main(int argc, char* argv[]) {
 
   double pZ = 1.0 - exp(-1.0 / T);
 
-  std::function<double(spin<signed, D, signed>, spin<signed, D, signed>, spin<signed, D, signed>)> Z =
-    [L, pZ] (spin<signed, D, signed> s1, spin<signed, D, signed> s2, spin<signed, D, signed> s1_new) -> double {
+  std::function<double(spin<signed, D, signed>, spin<signed, D, signed>)> Z =
+    [L, pZ] (spin<signed, D, signed> s1, spin<signed, D, signed> s2) -> double {
       bool old_one_one = false;
       bool old_many_ones = false;
       bool old_any_two = false;
-      bool new_one_one = false;
-      bool new_many_ones = false;
-      bool new_any_two = false;
 
       vector<signed, D> old_diff = diff<signed, D>(L, s1.x, s2.x);
-      vector<signed, D> new_diff = diff<signed, D>(L, s1_new.x, s2.x);
       
       for (unsigned i = 0; i < D; i++) {
         if (old_diff(i) == 1 && !old_one_one) {
@@ -74,38 +70,15 @@ int main(int argc, char* argv[]) {
         } else if (old_diff(i) > 1) {
           old_any_two = true;
         }
-        if (new_diff(i) == 1 && !new_one_one) {
-          new_one_one = true;
-        } else if (new_diff(i) == 1 && new_one_one) {
-          new_many_ones = true;
-        } else if (new_diff(i) > 1) {
-          new_any_two = true;
-        }
       }
 
       bool were_on_someone = !old_one_one && !old_any_two;
-      bool are_on_someone = !new_one_one && !new_any_two;
       bool were_nearest_neighbors = old_one_one && !old_many_ones && !old_any_two;
-      bool are_nearest_neighbors = new_one_one && !new_many_ones && !new_any_two;
 
       if (were_on_someone) {
-        return 0.0;
-      } else if (are_on_someone) {
-        return 1.0;
-      } else if (were_nearest_neighbors && are_nearest_neighbors) {
-        return 0.0;
+        return -std::numeric_limits<double>::infinity();;
       } else if (were_nearest_neighbors) {
-        if (s1.s * s2.s == 1) {
-          return pZ;
-        } else {
-          return 0.0;
-        }
-      } else if (are_nearest_neighbors) {
-        if (s1_new.s * s2.s == -1) {
-          return pZ;
-        } else {
-          return 0.0;
-        }
+        return s1.s * s2.s;
       } else {
         return 0.0;
       }
@@ -141,7 +114,7 @@ int main(int argc, char* argv[]) {
       return nn;
     };
 
-  model<signed, D, signed> ising(L, Z, B, neighbors);
+  model<signed, D, signed> ising(L, L, Z, B, neighbors);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
@@ -160,7 +133,7 @@ int main(int argc, char* argv[]) {
         ising.s.push_back({{i, j}, -1});
         down++;
       }
-      ising.dict.record<signed>({i, j}, n);
+      ising.dict.record({i, j}, n);
       n++;
     }
   }
diff --git a/space_wolff.hpp b/space_wolff.hpp
index c87697e..6ac9788 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -1,44 +1,39 @@
 
-#include <vector>
-#include <list>
-#include <set>
-#include <iostream>
+#include "randutils/randutils.hpp"
+#include <cmath>
+#include <eigen3/Eigen/Dense>
 #include <fstream>
 #include <functional>
-#include <random>
+#include <iostream>
+#include <list>
 #include <queue>
-#include <cmath>
-#include <eigen3/Eigen/Dense>
-#include "randutils/randutils.hpp"
+#include <random>
+#include <set>
+#include <vector>
 
-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>;
-
-template <class T, int D>
-vector<T, D> diff(T L, vector<T, D> v1, vector<T, D> v2) {
+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>;
+
+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++) {
@@ -51,425 +46,411 @@ vector<T, D> diff(T L, vector<T, D> v1, vector<T, D> v2) {
   return v;
 }
 
-template <class T, int D, class state>
-class spin {
-  public:
-    vector<T, D> x;
-    state s;
+template <class T, int D, class state> class spin {
+public:
+  vector<T, D> x;
+  state s;
 };
 
-template <class T, int D>
-class euclidean {
-  private:
-    T L;
-
-  public:
-    vector<T, D> t;
-    matrix<T, D> r;
-    euclidean(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;
-        }
+template <class T, int D> class euclidean {
+private:
+  T L;
+
+public:
+  vector<T, D> t;
+  matrix<T, D> r;
+  euclidean(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;
       }
     }
+  }
 
-    euclidean(T L, vector<T, D> t0, matrix<T, D> r0) : L(L) {
-      t = t0;
-      r = r0;
-    }
+  euclidean(T L, vector<T, D> t0, matrix<T, D> r0) : L(L) {
+    t = t0;
+    r = r0;
+  }
 
-    template <class state>
-    spin<T, D, state> act(const spin<T, D, state>& s) const {
-      spin<T, D, state> s_new;
+  template <class state> spin<T, D, state> act(const spin<T, D, state>& s) const {
+    spin<T, D, state> s_new;
 
-      s_new.x = t + r * s.x;
-      s_new.s = s.s;
+    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;
+    for (unsigned i = 0; i < D; i++) {
+      s_new.x(i) = fmod(L + s_new.x(i), L);
     }
 
-    euclidean act(const euclidean& 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);
-      }
+    return s_new;
+  }
 
-      euclidean pnew(this->L, tnew, rnew);
+  euclidean act(const euclidean& x) const {
+    vector<T, D> tnew = r * x.t + t;
+    matrix<T, D> rnew = r * x.r;
 
-      return pnew;
+    for (unsigned i = 0; i < D; i++) {
+      tnew(i) = fmod(L + tnew(i), L);
     }
 
-    euclidean inverse() const {
-      vector<T, D> tnew = - r.transpose() * t;
-      matrix<T, D> rnew = r.transpose();
+    euclidean pnew(this->L, tnew, rnew);
+
+    return pnew;
+  }
 
-      euclidean pnew(this->L, tnew, rnew);
+  euclidean inverse() const {
+    vector<T, D> tnew = -r.transpose() * t;
+    matrix<T, D> rnew = r.transpose();
 
-      return pnew;
-    }
+    euclidean pnew(this->L, tnew, rnew);
+
+    return pnew;
+  }
 };
 
-template <int D>
-class dictionary {
-  private:
-    unsigned L;
-    std::vector<std::set<unsigned>> d;
+template <class T, int D> class dictionary {
+private:
+  unsigned N;
+  T L;
+  std::vector<std::set<unsigned>> d;
 
-  public:
-    dictionary(unsigned Li) : L(Li), d(pow(Li, D)) {};
+public:
+  dictionary(unsigned Ni, double Li) : N(Ni), L(Li), d(pow(Ni, D)) {};
 
-    template <class T>
-    unsigned dictionary_index(vector<T, D> x) const {
-      unsigned pos_ind = 0;
+  unsigned dictionary_index(vector<T, D> x) const {
+    unsigned pos_ind = 0;
 
-      for (unsigned i = 0; i < D; i++) {
-        pos_ind += pow(L, i) * (unsigned)x(i);
-      };
+    for (unsigned i = 0; i < D; i++) {
+      pos_ind += pow(N, i) * (unsigned)std::floor(x(i) * N / L);
+    };
 
-      return pos_ind;
-    }
+    return pos_ind;
+  }
 
-    template <class T>
-    void record(vector<T, D> x, unsigned ind) {
-      d[dictionary_index<T>(x)].insert(ind);
-    };
+  void record(vector<T, D> x, unsigned ind) {
+    d[this->dictionary_index(x)].insert(ind);
+  };
 
-    template <class T>
-    void remove(vector<T, D> x, unsigned ind) {
-      d[dictionary_index<T>(x)].erase(ind);
-    };
+  void remove(vector<T, D> x, unsigned ind) {
+    d[this->dictionary_index(x)].erase(ind);
+  };
 
-    template <class T>
-    std::set<unsigned> neighbors(vector<T, D> x, unsigned depth) const {
-      return nearest_neighbors_of(dictionary_index<T>(x), depth, {});
-    };
+  std::set<unsigned> neighbors(vector<T, D> x, unsigned depth) const {
+    return nearest_neighbors_of(this->dictionary_index(x), depth, {});
+  };
 
-    std::set<unsigned> nearest_neighbors_of(unsigned ind, unsigned depth, std::list<unsigned> ignores) const {
-      std::set<unsigned> ns = d[ind];
+  std::set<unsigned> nearest_neighbors_of(unsigned ind, unsigned depth,
+                                          std::list<unsigned> ignores) const {
+    std::set<unsigned> ns = d[ind];
 
-      if (depth > 0) {
-        for (unsigned i = 0; i < D; i++) {
-          if (std::none_of(ignores.begin(), ignores.end(), [i](unsigned k){return i==k;})) {
-            std::list<unsigned> ignores_new = ignores;
-            ignores_new.push_back(i);
-            for (signed j : {-1, 1}) {
-              unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) + fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1));
+    if (depth > 0) {
+      for (unsigned i = 0; i < D; i++) {
+        if (std::none_of(ignores.begin(), ignores.end(), [i](unsigned k) { return i == k; })) {
+          std::list<unsigned> ignores_new = ignores;
+          ignores_new.push_back(i);
+          for (signed j : {-1, 1}) {
+            unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) +
+                          fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1));
 
-              std::set<unsigned> nns = nearest_neighbors_of(ni, depth - 1, ignores_new);
+            std::set<unsigned> nns = nearest_neighbors_of(ni, depth - 1, ignores_new);
 
-              for (unsigned guy : nns) {
-                ns.insert(guy);
-              }
+            for (unsigned guy : nns) {
+              ns.insert(guy);
             }
           }
         }
       }
+    }
 
-      return ns;
-    };
+    return ns;
+  };
 };
 
 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;
-    }
+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);
+  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);
     }
+    n++;
+  }
 
-    double avg2() const {
-      return total2 / (n - wait);
-    }
+  double avg() const { return total / (n - wait); }
 
-    std::vector<double> ρ() const {
-      double C0 = C.front() / (n - wait);
-      double avg2 = pow(total / (n - wait), 2);
+  double avg2() const { return total2 / (n - wait); }
 
-      std::vector<double> ρtmp;
+  std::vector<double> ρ() const {
+    double C0 = C.front() / (n - wait);
+    double avg2 = pow(total / (n - wait), 2);
 
-      for (double Ct : C) {
-        ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
-      }
+    std::vector<double> ρtmp;
 
-      return ρtmp;
+    for (double Ct : C) {
+      ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
     }
 
-    std::array<double, 2> τ() const {
-      double τtmp = 0.5;
-      unsigned M = 1;
-      double c = 8.0;
+    return ρtmp;
+  }
 
-      std::vector<double> ρ_tmp = this->ρ();
+  std::array<double, 2> τ() const {
+    double τtmp = 0.5;
+    unsigned M = 1;
+    double c = 8.0;
 
-      while (c * τtmp > M && M < C.size()) {
-        τtmp += ρ_tmp[M];
-        M++;
-      }
+    std::vector<double> ρ_tmp = this->ρ();
 
-      return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)};
+    while (c * τtmp > M && M < C.size()) {
+      τtmp += ρ_tmp[M];
+      M++;
     }
 
-    double σ() const {
-      return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2));
-    }
+    return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)};
+  }
 
-    double serr() const {
-      return sqrt(this->σ());
-    }
+  double σ() const {
+    return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2));
+  }
 
-    unsigned num_added() const {
-      return n - wait;
-    }
+  double serr() const { return sqrt(this->σ()); }
+
+  unsigned num_added() const { return n - wait; }
 };
 
-template <class U, int D, class state>
-class model {
-  public:
-    U L;
-    euclidean<U, D> s0;
-    std::vector<spin<U, D, state>> s;
-    dictionary<D> dict;
-    std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> neighbors;
-    std::function<double(spin<U, D, state>, spin<U, D, state>, spin<U, D, state>)> Z;
-    std::function<double(spin<U, D, state>)> B;
-    std::vector<matrix<U, D>> mats;
-    std::vector<vector<U, D>> steps;
-    long double E;
-    quantity Eq;
-    quantity Cq;
-
-    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(sequences, new_level);
+template <class U, int D, class state> class model {
+public:
+  U L;
+  euclidean<U, D> s0;
+  std::vector<spin<U, D, state>> s;
+  dictionary<U, D> dict;
+  std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> neighbors;
+  std::function<double(spin<U, D, state>, spin<U, D, state>)> Z;
+  std::function<double(spin<U, D, state>)> B;
+  std::vector<matrix<U, D>> mats;
+  std::vector<vector<U, D>> steps;
+  long double E;
+  quantity Eq;
+  quantity Cq;
+
+  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(sequences, new_level);
     }
+  }
 
-    model(U L, std::function<double(spin<U, D, state>, spin<U, D, state>, spin<U, D, state>)> Z,
+  model(U L, unsigned N, std::function<double(spin<U, D, state>, spin<U, D, state>)> Z,
         std::function<double(spin<U, D, state>)> B,
-        std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> ns) : 
-      L(L), s0(L), dict(L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000), Cq(1000, 1000) {
-        std::array<double, D> ini_sequence;
-        ini_sequence.fill(1);
-        std::list<std::array<double, D>> sequences;
-        sequences.push_back(ini_sequence);
+        std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> ns)
+      : L(L), s0(L), dict(N, L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000), Cq(1000, 1000) {
+    std::array<double, D> ini_sequence;
+    ini_sequence.fill(1);
+    std::list<std::array<double, D>> sequences;
+    sequences.push_back(ini_sequence);
 
-        one_sequences(sequences, D);
+    one_sequences(sequences, D);
 
-        sequences.pop_front(); // don't want the identity matrix!
+    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;
-              }
-            }
+    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);
-
-          vector<U, D> v;
-          for (unsigned i = 0; i < D; i++) {
-            if (sequence[i] == 1) {
-              v(i) = 0;
-            } else {
-              v(i) = L / 2;
-            }
-          }
+      mats.push_back(m);
 
-          steps.push_back(v);
+      vector<U, D> v;
+      for (unsigned i = 0; i < D; i++) {
+        if (sequence[i] == 1) {
+          v(i) = 0;
+        } else {
+          v(i) = L / 2;
         }
+      }
 
-        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;
-                  }
+      steps.push_back(v);
+    }
+
+    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);
             }
           }
+          mats.push_back(m);
         }
       }
+    }
+  }
 
-    void update_energy() {
-      E = 0;
-      for (unsigned i = 0; i < s.size(); i++) {
-        for (unsigned j = 0; j < i; j++) {
-        //  E -= Z(s[i], s[j]);
-        }
-
-        E -= B(s0.inverse().act(s[i]));
+  void update_energy() {
+    E = 0;
+    for (unsigned i = 0; i < s.size(); i++) {
+      for (unsigned j = 0; j < i; j++) {
+         E -= Z(s[i], s[j]);
       }
-    }
 
-    void step(double T, unsigned ind, euclidean<U, D> r, std::mt19937& rng) {
-      unsigned cluster_size = 0;
-      std::uniform_real_distribution<double> dist(0.0, 1.0);
+      E -= B(s0.inverse().act(s[i]));
+    }
+  }
 
-      std::queue<unsigned> queue;
-      queue.push(ind);
+  void step(double T, unsigned ind, euclidean<U, D> r, std::mt19937& rng) {
+    unsigned cluster_size = 0;
+    std::uniform_real_distribution<double> dist(0.0, 1.0);
 
-      std::vector<bool> visited(s.size() + 1, false);
+    std::queue<unsigned> queue;
+    queue.push(ind);
 
-      while (!queue.empty()) {
-        unsigned i = queue.front();
-        queue.pop();
+    std::vector<bool> visited(s.size() + 1, false);
 
-        if (!visited[i]) {
-          visited[i] = true;
+    while (!queue.empty()) {
+      unsigned i = queue.front();
+      queue.pop();
 
-          bool we_are_ghost = i == s.size();
+      if (!visited[i]) {
+        visited[i] = true;
 
-          spin<U, D, state> si_new;
-          euclidean<U, D> s0_new(L);
+        bool we_are_ghost = i == s.size();
 
-          if (we_are_ghost) {
-            s0_new = r.act(s0);
-          } else {
-            si_new = r.act(s[i]);
-          }
+        spin<U, D, state> si_new;
+        euclidean<U, D> s0_new(L);
 
-          for (unsigned j : neighbors(*this, i, si_new)) {
-            if (j != i) {
-              double p;
-              bool neighbor_is_ghost = j == s.size();
+        if (we_are_ghost) {
+          s0_new = r.act(s0);
+        } else {
+          si_new = r.act(s[i]);
+        }
 
-              if (we_are_ghost || neighbor_is_ghost) {
-                spin<U, D, state> s0s_old, s0s_new;
-                unsigned non_ghost;
+        for (unsigned j : neighbors(*this, i, si_new)) {
+          if (j != i) {
+            double p;
+            bool neighbor_is_ghost = j == s.size();
 
-                if (neighbor_is_ghost) {
-                  non_ghost = i;
-                  s0s_old = s0.inverse().act(s[i]);
-                  s0s_new = s0.inverse().act(si_new);
-                } else {
-                  non_ghost = j;
-                  s0s_old = s0.inverse().act(s[j]);
-                  s0s_new = s0_new.inverse().act(s[j]);
-                }
+            if (we_are_ghost || neighbor_is_ghost) {
+              spin<U, D, state> s0s_old, s0s_new;
+              unsigned non_ghost;
 
-                p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T);
+              if (neighbor_is_ghost) {
+                non_ghost = i;
+                s0s_old = s0.inverse().act(s[i]);
+                s0s_new = s0.inverse().act(si_new);
               } else {
-                p = Z(s[i], s[j], si_new);
+                non_ghost = j;
+                s0s_old = s0.inverse().act(s[j]);
+                s0s_new = s0_new.inverse().act(s[j]);
               }
 
-              if (dist(rng) < p) {
-                queue.push(j);
-              }
+              p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T);
+            } else {
+              p = 1.0 - exp(-(Z(si_new, s[j]) - Z(s[i], s[j])) / T);
             }
-          }
 
-          if (we_are_ghost) {
-            s0 = s0_new;
-          } else {
-            dict.remove(s[i].x, i);
-            s[i] = si_new;
-            dict.record(s[i].x, i);
-            cluster_size++;
+            if (dist(rng) < p) {
+              queue.push(j);
+            }
           }
         }
-      }
 
-      Cq.add(cluster_size);
+        if (we_are_ghost) {
+          s0 = s0_new;
+        } else {
+          dict.remove(s[i].x, i);
+          s[i] = si_new;
+          dict.record(s[i].x, i);
+          cluster_size++;
+        }
+      }
     }
 
-    void wolff(double T, unsigned N, std::mt19937& rng) {
-      std::uniform_real_distribution<double> t_dist(0, L);
-      std::uniform_int_distribution<unsigned> r_dist(0, D - 1);
-      std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1);
-      std::uniform_int_distribution<unsigned> coin(0, mats.size() + steps.size() - 1);
-
-      for (unsigned i = 0; i < N; i++) {
-        vector<U, D> t;
-        matrix<U, D> m;
-        unsigned flip = coin(rng);
-        if (flip < mats.size()) {
-          for (unsigned j = 0; j < D; j++) {
-            t(j) = (U)t_dist(rng);
-          }
-          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;
-              }
+    Cq.add(cluster_size);
+  }
+
+  void wolff(double T, unsigned N, std::mt19937& rng) {
+    std::uniform_real_distribution<double> t_dist(0, L);
+    std::uniform_int_distribution<unsigned> r_dist(0, D - 1);
+    std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1);
+    std::uniform_int_distribution<unsigned> coin(0, mats.size() + steps.size() - 1);
+
+    for (unsigned i = 0; i < N; i++) {
+      vector<U, D> t;
+      matrix<U, D> m;
+      unsigned flip = coin(rng);
+      if (flip < mats.size()) {
+        for (unsigned j = 0; j < D; j++) {
+          t(j) = (U)t_dist(rng);
+        }
+        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 = steps[flip - mats.size()];
         }
 
-        euclidean<U, D> g(L, t, m);
+        t = steps[flip - mats.size()];
+      }
 
-        this->step(T, ind_dist(rng), g, rng);
+      euclidean<U, D> g(L, t, m);
 
-        this->update_energy();
-        Eq.add(E);
-      }
+      this->step(T, ind_dist(rng), g, rng);
+
+      this->update_energy();
+      Eq.add(E);
     }
+  }
 };