From 7b4ab37cc6c728b399e9279e186c9d6c596ef523 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 26 Jul 2019 15:20:27 -0400
Subject: split into different files, added command line interface, and
 implemented hard spheres

---
 ising.cpp       | 141 ++++++++++++++++++
 space_wolff.cpp | 446 --------------------------------------------------------
 space_wolff.hpp | 333 ++++++++++++++++++++++++++++++++++++++++++
 spheres.cpp     | 104 +++++++++++++
 4 files changed, 578 insertions(+), 446 deletions(-)
 create mode 100644 ising.cpp
 delete mode 100644 space_wolff.cpp
 create mode 100644 space_wolff.hpp
 create mode 100644 spheres.cpp

diff --git a/ising.cpp b/ising.cpp
new file mode 100644
index 0000000..6939749
--- /dev/null
+++ b/ising.cpp
@@ -0,0 +1,141 @@
+
+#include "space_wolff.hpp"
+
+int main(int argc, char* argv[]) {
+  const unsigned D = 2;
+
+  unsigned L = 32;
+  unsigned N = 1000;
+  double T = 2.0 / log(1.0 + sqrt(2.0));
+  double H = 1.0;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "N:L:T:H:")) != -1) {
+    switch (opt) {
+      case 'N': 
+        N = (unsigned)atof(optarg);
+        break;
+      case 'L':
+        L = atoi(optarg);
+        break;
+      case 'T':
+        T = atof(optarg);
+        break;
+      case 'H':
+        H = atof(optarg);
+        break;
+      default:
+        exit(1);
+    }
+  }
+
+  std::function<double(spin<signed, D, signed>, spin<signed, D, signed>)> Z =
+    [] (spin<signed, D, signed> s1, spin<signed, D, signed> s2) -> double {
+      bool one_one = false;
+      bool many_ones = false;
+      bool any_two = false;
+      
+      for (unsigned i = 0; i < D; i++) {
+        unsigned diff = abs(s1.x(i) - s2.x(i));
+        if (diff == 1 && !one_one) {
+          one_one = true;
+        } else if (diff == 1 && one_one) {
+          many_ones = true;
+          break;
+        } else if (diff > 1) {
+          any_two = true;
+          break;
+        }
+      }
+
+      if (!one_one && !any_two) {
+        return -std::numeric_limits<double>::infinity();
+      } else if (one_one && !many_ones && !any_two) {
+        return s1.s * s2.s;
+      } else {
+        return 0;
+      }
+    };
+
+  std::function<double(spin<signed, D, signed>)> B =
+    [L, H] (spin<signed, D, signed> s) -> double {
+      return H * s.s * smiley[s.x(1) * 16 / L][s.x(0) * 16 / L];
+    };
+
+  std::function<std::set<unsigned>(model<signed, D, signed>&, unsigned, spin<signed, D, signed>)> neighbors =
+    [] (model<signed, D, signed>& m, unsigned i0, spin<signed, D, signed> s1) -> std::set<unsigned> {
+      std::set<unsigned> nn;
+      if (i0 < m.s.size()) {
+        std::set<unsigned> os1 = m.dict.on_site(s1.x);
+        std::set<unsigned> nn0 = m.dict.nearest_neighbors(m.s[i0].x);
+        std::set<unsigned> nn1 = m.dict.nearest_neighbors(s1.x);
+        nn.insert(nn0.begin(), nn0.end());
+        nn.insert(nn1.begin(), nn1.end());
+        nn.insert(os1.begin(), os1.end());
+        nn.insert(m.s.size());
+      } else {
+        for (unsigned i = 0; i < m.s.size(); i++) {
+          nn.insert(i);
+        }
+      }
+      return nn;
+    };
+
+  model<signed, D, signed> ising(L, Z, B, neighbors);
+
+  randutils::auto_seed_128 seeds;
+  std::mt19937 rng{seeds};
+
+  std::uniform_int_distribution<unsigned> coin(0, 1);
+
+  unsigned n = 0;
+  unsigned up = 0;
+  unsigned down = 0;
+  for (unsigned i = 0; i < L; i++) {
+    for (unsigned j = 0; j < L; j++) {
+      if ((coin(rng) && up < pow(L, 2) / 2) || down >= pow(L, 2) / 2) {
+        ising.s.push_back({{i, j}, 1});
+        up++;
+      } else {
+        ising.s.push_back({{i, j}, -1});
+        down++;
+      }
+      ising.dict.record<signed>({i, j}, n);
+      n++;
+    }
+  }
+  /*
+  for (unsigned i = 0; i < L; i++) {
+    for (unsigned j = 0; j < L; j++) {
+      if (i < L / 2) {
+        ising.s.push_back({{i, j}, 1});
+      } else {
+        ising.s.push_back({{i, j}, -1});
+      }
+      ising.dict.record<signed>({i, j}, n);
+      n++;
+    }
+  }
+  */
+
+  ising.wolff(T, N, rng);
+
+  std::vector<signed> output(pow(L, D));
+
+  for (spin<signed, D, signed> s : ising.s) {
+    spin<signed, D, signed> rs = ising.s0.inverse().act(s);
+    output[L * rs.x(1) + rs.x(0)] = s.s;
+  }
+
+  for (unsigned i = 0; i < L; i++) {
+    for (unsigned j = 0; j < L; j++) {
+      unsigned out = output[L * i + j] == 1 ? 1 : 0;
+      std::cout << out;
+    }
+    std::cout << "\n";
+  }
+
+  return 0;
+}
+
diff --git a/space_wolff.cpp b/space_wolff.cpp
deleted file mode 100644
index b88944e..0000000
--- a/space_wolff.cpp
+++ /dev/null
@@ -1,446 +0,0 @@
-
-#include <vector>
-#include <list>
-#include <set>
-#include <iostream>
-#include <functional>
-#include <random>
-#include <queue>
-#include <eigen3/Eigen/Dense>
-#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>;
-
-template <class T, int D, class state>
-class spin {
-  public:
-    vector<T, D> x;
-    state s;
-};
-
-template <class T, int D>
-class euclidean {
-  private:
-    unsigned L;
-
-  public:
-    vector<T, D> t;
-    matrix<T, D> r;
-    euclidean(unsigned 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(unsigned L, vector<T, D> t0, matrix<T, D> r0) : L(L) {
-      t = t0;
-      r = r0;
-    }
-
-    template <class state>
-    spin<T, D, state> act(spin<T, D, state> s) {
-      spin<T, D, state> 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;
-    }
-
-    euclidean act(const euclidean& x) {
-      vector<T, D> tnew = r * x.t + t;
-      matrix<T, D> rnew = r * x.r;
-
-      euclidean pnew(this->L, tnew, rnew);
-
-      return pnew;
-    }
-
-    euclidean inverse() {
-      vector<T, D> tnew = - r.transpose() * t;
-      matrix<T, D> rnew = r.transpose();
-
-      euclidean pnew(this->L, tnew, rnew);
-
-      return pnew;
-    }
-};
-
-template <int D>
-class dictionary {
-  private:
-    unsigned L;
-    std::vector<std::set<unsigned>> d;
-
-  public:
-    dictionary(unsigned Li) : L(Li), d(pow(Li, D)) {};
-
-    template <class T>
-    unsigned dictionary_index(vector<T, D> x) {
-      unsigned pos_ind = 0;
-
-      for (unsigned i = 0; i < D; i++) {
-        pos_ind += pow(L, i) * (unsigned)x(i);
-      };
-
-      return pos_ind;
-    }
-
-    template <class T>
-    void record(vector<T, D> x, unsigned ind) {
-      d[dictionary_index<T>(x)].insert(ind);
-    };
-
-    template <class T>
-    void remove(vector<T, D> x, unsigned ind) {
-      d[dictionary_index<T>(x)].erase(ind);
-    };
-
-    template <class T>
-    std::set<unsigned> on_site(vector<T, D> x) {
-      return d[dictionary_index<T>(x)];
-    };
-
-    template <class T>
-    std::set<unsigned> nearest_neighbors(vector<T, D> x) {
-      unsigned ind = dictionary_index<T>(x);
-      std::set<unsigned> ns;
-
-      for (unsigned i = 0; i < D; 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));
-          for (unsigned nii : d[ni]) {
-            ns.insert(nii);
-          }
-        }
-      }
-
-      return ns;
-    };
-
-    template <class T>
-    std::set<unsigned> next_nearest_neighbors(vector<T, D> x) {
-      unsigned ind = dictionary_index<T>(x);
-      std::set<unsigned> ns;
-
-      for (unsigned i = 0; i < D; 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));
-          for (unsigned k = 0; k < D; k++) {
-            if (k != i) {
-              for (signed l : {-1, 1}) {
-                unsigned nni = pow(L, k + 1) * (ni / ((unsigned)pow(L, k + 1))) + fmod(pow(L, k + 1) + ni + l * pow(L, k), pow(L, k + 1));
-                for (unsigned nnii : d[nni]) {
-                  ns.insert(nnii);
-                }
-              }
-            }
-          }
-        }
-      }
-
-      return ns;
-    };
-
-    template <class T>
-    std::set<unsigned> next_next_nearest_neighbors(vector<T, D> x) {
-      unsigned ind = dictionary_index<T>(x);
-      std::set<unsigned> ns;
-
-      for (unsigned i = 0; i < D; 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));
-          for (unsigned k = 0; k < D; k++) {
-            if (k != i) {
-              for (signed l : {-1, 1}) {
-                unsigned nni = pow(L, k + 1) * (ni / ((unsigned)pow(L, k + 1))) + fmod(pow(L, k + 1) + ni + l * pow(L, k), pow(L, k + 1));
-                for (unsigned m = 0; m < D; m++) {
-                  if (m != i && m != k) {
-                    for (signed n : {-1, 1}) {
-                      unsigned nnni = pow(L, m + 1) * (nni / ((unsigned)pow(L, m + 1))) + fmod(pow(L, m + 1) + nni + n * pow(L, m), pow(L, m + 1));
-                      for (unsigned nnnii : d[nnni]) {
-                        ns.insert(nnnii);
-                      }
-                    }
-                  }
-                }
-              }
-            }
-          }
-        }
-      }
-
-      return ns;
-    };
-};
-
-template <class U, int D, class state>
-class model {
-  public:
-    unsigned 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>)> Z;
-    std::function<double(spin<U, D, state>)> B;
-    double E;
-
-    model(unsigned L, 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) {
-        }
-
-    void step(double T, unsigned ind, euclidean<U, D> r, std::mt19937& rng) {
-      std::uniform_real_distribution<double> dist(0.0, 1.0);
-
-      std::queue<unsigned> queue;
-      queue.push(ind);
-
-      std::vector<bool> visited(s.size() + 1, false);
-
-      while (!queue.empty()) {
-        unsigned i = queue.front();
-        queue.pop();
-
-        if (!visited[i]) {
-          visited[i] = true;
-
-          bool we_are_ghost = i == s.size();
-
-          spin<U, D, state> si_new;
-          euclidean<U, D> s0_new(L);
-
-          if (we_are_ghost) {
-            s0_new = r.act(s0);
-          } else {
-            si_new = r.act(s[i]);
-          }
-
-          for (unsigned j : neighbors(*this, i, si_new)) {
-            if (j != i) {
-              double dE;
-              bool neighbor_is_ghost = j == s.size();
-
-              if (we_are_ghost || neighbor_is_ghost) {
-                spin<U, D, state> s0s_old, s0s_new;
-                unsigned non_ghost;
-
-                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]);
-                }
-
-                dE = B(s0s_old) - B(s0s_new);
-              } else {
-                dE = Z(s[i], s[j]) - Z(si_new, s[j]);
-              }
-
-              double p = 1.0 - exp(-dE / T);
-
-              if (dist(rng) < p) {
-                queue.push(j);
-              }
-            }
-          }
-
-          if (we_are_ghost) {
-            s0 = s0_new;
-          } else {
-            dict.remove(s[i].x, i);
-            s[i] = si_new;
-            dict.record(s[i].x, i);
-          }
-        }
-      }
-    }
-
-    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);
-
-      for (unsigned i = 0; i < N; i++) {
-        vector<U, D> t;
-        matrix<U, D> m;
-        for (unsigned j = 0; j < D; j++) {
-          t(j) = (U)t_dist(rng);
-        }
-
-        unsigned flip_D1 = r_dist(rng);
-        unsigned flip_D2 = r_dist(rng);
-
-        for (unsigned j = 0; j < D; j++) {
-          for (unsigned k = 0; k < D; k++) {
-            if ((j == flip_D1 && k == flip_D2) || (j == flip_D2 && k == flip_D1)) {
-              if (flip_D1 <= flip_D2) {
-                m(j, k) = -1;
-              } else {
-                m(j, k) = 1;
-              }
-            } else if ((j == k && j != flip_D1) && j != flip_D2) {
-              m(j, k) = 1;
-            } else {
-              m(j, k) = 0;
-            }
-          }
-        }
-
-        euclidean<U, D> g(L, t, m);
-
-        this->step(T, ind_dist(rng), g, rng);
-      }
-    }
-};
-
-int main(int argc, char* argv[]) {
-  unsigned L = 32;
-  const unsigned D = 2;
-
-  std::function<double(spin<signed, D, signed>, spin<signed, D, signed>)> Z =
-    [] (spin<signed, D, signed> s1, spin<signed, D, signed> s2) -> double {
-      bool one_one = false;
-      bool many_ones = false;
-      bool any_two = false;
-      
-      for (unsigned i = 0; i < D; i++) {
-        unsigned diff = abs(s1.x(i) - s2.x(i));
-        if (diff == 1 && !one_one) {
-          one_one = true;
-        } else if (diff == 1 && one_one) {
-          many_ones = true;
-          break;
-        } else if (diff > 1) {
-          any_two = true;
-          break;
-        }
-      }
-
-      if (!one_one && !any_two) {
-        return -std::numeric_limits<double>::infinity();
-      } else if (one_one && !many_ones && !any_two) {
-        return s1.s * s2.s;
-      } else {
-        return 0;
-      }
-    };
-
-  std::function<double(spin<signed, D, signed>)> B =
-    [] (spin<signed, D, signed> s) -> double {
-      return 100 * s.s * smiley[s.x(0) / 2][s.x(1) / 2];
-    };
-
-  std::function<std::set<unsigned>(model<signed, D, signed>&, unsigned, spin<signed, D, signed>)> neighbors =
-    [] (model<signed, D, signed>& m, unsigned i0, spin<signed, D, signed> s1) -> std::set<unsigned> {
-      std::set<unsigned> nn;
-      if (i0 < m.s.size()) {
-        std::set<unsigned> os1 = m.dict.on_site(s1.x);
-        std::set<unsigned> nn0 = m.dict.nearest_neighbors(m.s[i0].x);
-        std::set<unsigned> nn1 = m.dict.nearest_neighbors(s1.x);
-        nn.insert(nn0.begin(), nn0.end());
-        nn.insert(nn1.begin(), nn1.end());
-        nn.insert(os1.begin(), os1.end());
-        nn.insert(m.s.size());
-      } else {
-        for (unsigned i = 0; i < m.s.size(); i++) {
-          nn.insert(i);
-        }
-      }
-      return nn;
-    };
-
-  model<signed, D, signed> ising(L, Z, B, neighbors);
-
-  randutils::auto_seed_128 seeds;
-  std::mt19937 rng{seeds};
-
-  std::uniform_int_distribution<unsigned> coin(0, 1);
-
-  unsigned n = 0;
-  unsigned up = 0;
-  unsigned down = 0;
-  for (unsigned i = 0; i < L; i++) {
-    for (unsigned j = 0; j < L; j++) {
-      if ((coin(rng) && up < pow(L, 2) / 2) || down >= pow(L, 2) / 2) {
-        ising.s.push_back({{i, j}, 1});
-        up++;
-      } else {
-        ising.s.push_back({{i, j}, -1});
-        down++;
-      }
-      ising.dict.record<signed>({i, j}, n);
-      n++;
-    }
-  }
-  /*
-  for (unsigned i = 0; i < L; i++) {
-    for (unsigned j = 0; j < L; j++) {
-      if (i < L / 2) {
-        ising.s.push_back({{i, j}, 1});
-      } else {
-        ising.s.push_back({{i, j}, -1});
-      }
-      ising.dict.record<signed>({i, j}, n);
-      n++;
-    }
-  }
-  */
-
-  ising.wolff(2.0 / log(1.0 + sqrt(2.0)), 5000, rng);
-
-  std::vector<signed> output(pow(L, D));
-
-  for (spin<signed, D, signed> s : ising.s) {
-    spin<signed, D, signed> rs = ising.s0.inverse().act(s);
-    output[L * rs.x(1) + rs.x(0)] = s.s;
-  }
-
-  for (unsigned i = 0; i < L; i++) {
-    for (unsigned j = 0; j < L; j++) {
-      unsigned out = output[L * i + j] == 1 ? 1 : 0;
-      std::cout << out;
-    }
-    std::cout << "\n";
-  }
-
-  return 0;
-}
-
diff --git a/space_wolff.hpp b/space_wolff.hpp
new file mode 100644
index 0000000..aaef673
--- /dev/null
+++ b/space_wolff.hpp
@@ -0,0 +1,333 @@
+
+#include <vector>
+#include <list>
+#include <set>
+#include <iostream>
+#include <functional>
+#include <random>
+#include <queue>
+#include <eigen3/Eigen/Dense>
+#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>;
+
+template <class T, int D, class state>
+class spin {
+  public:
+    vector<T, D> x;
+    state s;
+};
+
+template <class T, int D>
+class euclidean {
+  private:
+    unsigned L;
+
+  public:
+    vector<T, D> t;
+    matrix<T, D> r;
+    euclidean(unsigned 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(unsigned L, vector<T, D> t0, matrix<T, D> r0) : L(L) {
+      t = t0;
+      r = r0;
+    }
+
+    template <class state>
+    spin<T, D, state> act(spin<T, D, state> s) {
+      spin<T, D, state> 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;
+    }
+
+    euclidean act(const euclidean& x) {
+      vector<T, D> tnew = r * x.t + t;
+      matrix<T, D> rnew = r * x.r;
+
+      euclidean pnew(this->L, tnew, rnew);
+
+      return pnew;
+    }
+
+    euclidean inverse() {
+      vector<T, D> tnew = - r.transpose() * t;
+      matrix<T, D> rnew = r.transpose();
+
+      euclidean pnew(this->L, tnew, rnew);
+
+      return pnew;
+    }
+};
+
+template <int D>
+class dictionary {
+  private:
+    unsigned L;
+    std::vector<std::set<unsigned>> d;
+
+  public:
+    dictionary(unsigned Li) : L(Li), d(pow(Li, D)) {};
+
+    template <class T>
+    unsigned dictionary_index(vector<T, D> x) {
+      unsigned pos_ind = 0;
+
+      for (unsigned i = 0; i < D; i++) {
+        pos_ind += pow(L, i) * (unsigned)x(i);
+      };
+
+      return pos_ind;
+    }
+
+    template <class T>
+    void record(vector<T, D> x, unsigned ind) {
+      d[dictionary_index<T>(x)].insert(ind);
+    };
+
+    template <class T>
+    void remove(vector<T, D> x, unsigned ind) {
+      d[dictionary_index<T>(x)].erase(ind);
+    };
+
+    template <class T>
+    std::set<unsigned> on_site(vector<T, D> x) {
+      return d[dictionary_index<T>(x)];
+    };
+
+    template <class T>
+    std::set<unsigned> nearest_neighbors(vector<T, D> x) {
+      unsigned ind = dictionary_index<T>(x);
+      std::set<unsigned> ns;
+
+      for (unsigned i = 0; i < D; 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));
+          for (unsigned nii : d[ni]) {
+            ns.insert(nii);
+          }
+        }
+      }
+
+      return ns;
+    };
+
+    template <class T>
+    std::set<unsigned> next_nearest_neighbors(vector<T, D> x) {
+      unsigned ind = dictionary_index<T>(x);
+      std::set<unsigned> ns;
+
+      for (unsigned i = 0; i < D; 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));
+          for (unsigned k = 0; k < D; k++) {
+            if (k != i) {
+              for (signed l : {-1, 1}) {
+                unsigned nni = pow(L, k + 1) * (ni / ((unsigned)pow(L, k + 1))) + fmod(pow(L, k + 1) + ni + l * pow(L, k), pow(L, k + 1));
+                for (unsigned nnii : d[nni]) {
+                  ns.insert(nnii);
+                }
+              }
+            }
+          }
+        }
+      }
+
+      return ns;
+    };
+
+    template <class T>
+    std::set<unsigned> next_next_nearest_neighbors(vector<T, D> x) {
+      unsigned ind = dictionary_index<T>(x);
+      std::set<unsigned> ns;
+
+      for (unsigned i = 0; i < D; 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));
+          for (unsigned k = 0; k < D; k++) {
+            if (k != i) {
+              for (signed l : {-1, 1}) {
+                unsigned nni = pow(L, k + 1) * (ni / ((unsigned)pow(L, k + 1))) + fmod(pow(L, k + 1) + ni + l * pow(L, k), pow(L, k + 1));
+                for (unsigned m = 0; m < D; m++) {
+                  if (m != i && m != k) {
+                    for (signed n : {-1, 1}) {
+                      unsigned nnni = pow(L, m + 1) * (nni / ((unsigned)pow(L, m + 1))) + fmod(pow(L, m + 1) + nni + n * pow(L, m), pow(L, m + 1));
+                      for (unsigned nnnii : d[nnni]) {
+                        ns.insert(nnnii);
+                      }
+                    }
+                  }
+                }
+              }
+            }
+          }
+        }
+      }
+
+      return ns;
+    };
+};
+
+template <class U, int D, class state>
+class model {
+  public:
+    unsigned 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>)> Z;
+    std::function<double(spin<U, D, state>)> B;
+    double E;
+
+    model(unsigned L, 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) {
+        }
+
+    void step(double T, unsigned ind, euclidean<U, D> r, std::mt19937& rng) {
+      std::uniform_real_distribution<double> dist(0.0, 1.0);
+
+      std::queue<unsigned> queue;
+      queue.push(ind);
+
+      std::vector<bool> visited(s.size() + 1, false);
+
+      while (!queue.empty()) {
+        unsigned i = queue.front();
+        queue.pop();
+
+        if (!visited[i]) {
+          visited[i] = true;
+
+          bool we_are_ghost = i == s.size();
+
+          spin<U, D, state> si_new;
+          euclidean<U, D> s0_new(L);
+
+          if (we_are_ghost) {
+            s0_new = r.act(s0);
+          } else {
+            si_new = r.act(s[i]);
+          }
+
+          for (unsigned j : neighbors(*this, i, si_new)) {
+            if (j != i) {
+              double dE;
+              bool neighbor_is_ghost = j == s.size();
+
+              if (we_are_ghost || neighbor_is_ghost) {
+                spin<U, D, state> s0s_old, s0s_new;
+                unsigned non_ghost;
+
+                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]);
+                }
+
+                dE = B(s0s_old) - B(s0s_new);
+              } else {
+                dE = Z(s[i], s[j]) - Z(si_new, s[j]);
+              }
+
+              double p = 1.0 - exp(-dE / T);
+
+              if (dist(rng) < p) {
+                queue.push(j);
+              }
+            }
+          }
+
+          if (we_are_ghost) {
+            s0 = s0_new;
+          } else {
+            dict.remove(s[i].x, i);
+            s[i] = si_new;
+            dict.record(s[i].x, i);
+          }
+        }
+      }
+    }
+
+    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);
+
+      for (unsigned i = 0; i < N; i++) {
+        vector<U, D> t;
+        matrix<U, D> m;
+        for (unsigned j = 0; j < D; j++) {
+          t(j) = (U)t_dist(rng);
+        }
+
+        unsigned flip_D1 = r_dist(rng);
+        unsigned flip_D2 = r_dist(rng);
+
+        for (unsigned j = 0; j < D; j++) {
+          for (unsigned k = 0; k < D; k++) {
+            if ((j == flip_D1 && k == flip_D2) || (j == flip_D2 && k == flip_D1)) {
+              if (flip_D1 <= flip_D2) {
+                m(j, k) = -1;
+              } else {
+                m(j, k) = 1;
+              }
+            } else if ((j == k && j != flip_D1) && j != flip_D2) {
+              m(j, k) = 1;
+            } else {
+              m(j, k) = 0;
+            }
+          }
+        }
+
+        euclidean<U, D> g(L, t, m);
+
+        this->step(T, ind_dist(rng), g, rng);
+      }
+    }
+};
+
diff --git a/spheres.cpp b/spheres.cpp
new file mode 100644
index 0000000..106ce4e
--- /dev/null
+++ b/spheres.cpp
@@ -0,0 +1,104 @@
+
+#include "space_wolff.hpp"
+
+int main(int argc, char* argv[]) {
+  const unsigned D = 2;
+
+  double L = 32;
+  unsigned N = 1000;
+  double T = 2.0 / log(1.0 + sqrt(2.0));
+  double H = 1.0;
+  unsigned n = 25;
+
+  int opt;
+
+  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);
+    }
+  }
+
+  std::function<double(spin<double, D, double>, spin<double, D, double>)> Z =
+    [L] (spin<double, D, double> s1, spin<double, D, double> s2) -> double {
+      vector<double, D> diff = s1.x - s2.x;
+      for (unsigned i = 0; i < D; i++) {
+        if (fabs(diff(i)) > L / 2) {
+          diff(i) = L - fabs(diff(i));
+        } else {
+          diff(i) = fabs(diff(i));
+        }
+      }
+      if (diff.transpose() * diff < pow(s1.s + s2.s, 2)) {
+        return -std::numeric_limits<double>::infinity();
+      } else {
+        return 0;
+      }
+    };
+
+  std::function<double(spin<double, D, double>)> B =
+    [L, H] (spin<double, D, double> s) -> double {
+      return H * sin(2 * M_PI * 3 * s.x(0) / L);
+    };
+
+  std::function<std::set<unsigned>(model<double, D, double>&, unsigned, spin<double, D, double>)> neighbors =
+    [] (model<double, D, double>& m, unsigned i0, spin<double, D, double> s1) -> std::set<unsigned> {
+      std::set<unsigned> nn;
+      if (i0 < m.s.size()) {
+        std::set<unsigned> os1 = m.dict.on_site(s1.x);
+        std::set<unsigned> nn0 = m.dict.nearest_neighbors(m.s[i0].x);
+        std::set<unsigned> nn1 = m.dict.nearest_neighbors(s1.x);
+        std::set<unsigned> nnn0 = m.dict.next_nearest_neighbors(m.s[i0].x);
+        std::set<unsigned> nnn1 = m.dict.next_nearest_neighbors(s1.x);
+        nn.insert(nn0.begin(), nn0.end());
+        nn.insert(nn1.begin(), nn1.end());
+        nn.insert(nnn0.begin(), nnn0.end());
+        nn.insert(nnn1.begin(), nnn1.end());
+        nn.insert(os1.begin(), os1.end());
+        nn.insert(m.s.size());
+      } else {
+        for (unsigned i = 0; i < m.s.size(); i++) {
+          nn.insert(i);
+        }
+      }
+      return nn;
+    };
+
+  model<double, D, double> sphere(L, Z, B, neighbors);
+
+  randutils::auto_seed_128 seeds;
+  std::mt19937 rng{seeds};
+
+  std::uniform_real_distribution<double> dist(0.0, L);
+
+  for (unsigned i = 0; i < n; i++) {
+    vector<double, D> pos = {dist(rng), dist(rng)};
+    sphere.s.push_back({pos, 0.5});
+    sphere.dict.record<double>(pos, i);
+  }
+
+  sphere.wolff(T, N, rng);
+
+  for (spin<double, D, double> s : sphere.s) {
+    spin<double, D, double> rs = sphere.s0.inverse().act(s);
+    std::cout << s.x.transpose() << "\n";
+  }
+
+  return 0;
+}
+
-- 
cgit v1.2.3-54-g00ecf