Lots of work.

3 files changed, 601 insertions, 728 deletions

diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
index 7e200b9..e2fb71a 100644
--- a/biroli-mezard.cpp
+++ b/biroli-mezard.cpp
@@ -5,6 +5,7 @@
 #include <queue>
 #include <vector>
 #include <chrono>
+#include <unordered_set>
 
 #include <eigen3/Eigen/Dense>
 
@@ -152,6 +153,7 @@ public:
   unsigned occupiedNeighbors;
   unsigned maximumNeighbors;
   bool marked;
+  bool visited;
 
   Vertex() {
     occupiedNeighbors = 0;
@@ -181,6 +183,7 @@ template <int D> class System {
 public:
   const unsigned L;
   std::vector<unsigned> N;
+  std::unordered_set<Vertex<D>*> occupants;
   std::vector<Vertex<D>> vertices;
   Transformation<D> orientation;
 
@@ -322,6 +325,9 @@ public:
       }
       N[t]++;
       N[0]--;
+      if (t > 0) {
+        occupants.insert(&v);
+      }
       return true;
     } else {
       return false;
@@ -332,6 +338,9 @@ public:
     if (v.empty()) {
       return false;
     } else {
+      if (v.maximumNeighbors > 0) {
+        occupants.erase(&v);
+      }
       N[v.maximumNeighbors]--;
       N[0]++;
       v.maximumNeighbors = Empty;
@@ -418,30 +427,32 @@ public:
     return exchange(v1, v2);
   }
 
-  void sweepGrandCanonical(double z, Rng& r) {
-    for (unsigned i = 0; i < size(); i++) {
-      if (0.5 < r.uniform(0.0, 1.0)) {
-        double pIns = size() * z / (occupancy() + 1);
-
-        if (pIns > r.uniform(0.0, 1.0)) {
-          while (true) {
-            Vertex<D>& v = r.pick(vertices);
-            if (v.empty()) {
-              insert(v, r.pick({1, 2, 2, 2, 2, 2, 3, 3, 3, 3}));
-              break;
-            }
-          }
-        }
-      } else {
+  const std::array<double, 3> ratios = {0.1, 0.5, 0.4};
 
-        double pDel = density() / z;
+  void stepGrandCanonical(double z, Rng& r) {
+    if (1.0 / 11.0 < r.uniform(0.0, 1.0)) {
+      unsigned t = r.pick({1, 2, 2, 2, 2, 2, 3, 3, 3, 3});
+      double pIns = size() * z / (occupancy() + 1);
 
-        if (pDel > r.uniform(0.0, 1.0)) {
-          remove(r.pick(vertices));
-        }
+      if (pIns > r.uniform(0.0, 1.0)) {
+        Vertex<D>& v = r.pick(vertices);
+        insert(v, t);
       }
+    } else {
+      double pDel = density() / z;
 
-      randomLocalExchange(r);
+      if (pDel > r.uniform(0.0, 1.0)) {
+        Vertex<D>* v = r.pick(occupants);
+        remove(*v);
+      }
+    }
+  }
+
+  void sweepGrandCanonical(double z, Rng& r) {
+    for (unsigned i = 0; i < size(); i++) {
+      stepGrandCanonical(z, r);
+
+      randomExchange(r);
     }
   }
 
@@ -457,6 +468,71 @@ public:
     }
   }
 
+  Vertex<D>& apply(const Transformation<D>& R, const Vertex<D>& v) {
+    return vertices[vectorToIndex(R.apply(v.position))];
+  }
+
+  bool eventChain(const Transformation<D>& R, Vertex<D>& v0) {
+    unsigned n = 0;
+
+    if (!v0.empty()) {
+      Vertex<D>& v0New = apply(R, v0);;
+      unsigned t0 = v0.maximumNeighbors;
+
+      std::queue<std::pair<std::reference_wrapper<Vertex<D>>, unsigned>> q;
+
+      q.push({v0New, t0});
+      remove(v0);
+
+      while (!q.empty()) {
+        auto [vr, t] = q.front();
+        Vertex<D>& v = vr;
+        q.pop();
+
+        if (!v.empty()) {
+          q.push({apply(R, v), v.maximumNeighbors});
+          remove(v);
+        }
+
+        insert(v, t, true);
+        v.marked = true;
+
+        for (Vertex<D>& vn : overlaps({v})) {
+          if (!vn.marked) {
+            q.push({apply(R, vn), vn.maximumNeighbors});
+            remove(vn);
+            vn.marked = true;
+          } else {
+            /* If a neighbor has already been moved but is still
+             * frustrated, it may be due to one of its neighbors!
+             */
+            if (&vn != &v) {
+              for (Vertex<D>& vnn : overlaps(vn)) {
+                if (!vnn.marked) {
+                  q.push({apply(R, vnn), vnn.maximumNeighbors});
+                  remove(vnn);
+                  vnn.marked = true;
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+
+    for (Vertex<D>& v : vertices) {
+      v.marked = false;
+    }
+
+    return true;
+  }
+
+  bool randomEventChain(Rng& r) {
+    Transformation<D> R(L);
+    R.v[r.uniform((unsigned)0, (unsigned)D-1)] = r.pick({-1, 1});
+    return eventChain(R, r.pick(vertices));
+  }
+
   unsigned flipCluster(const Transformation<D>& R, Vertex<D>& v0) {
     unsigned n = 0;
 
@@ -546,30 +622,35 @@ int main() {
   unsigned L = 15;
   unsigned Nmin = 2e2;
   unsigned Nmax = 2e5;
-  double Tmin = 0.04;
-  double Tmax = 0.2;
-  double δT = 0.02;
+  double μmin = -3;
+  double μmax = 10;
+  double dμ = 0.05;
 
   System<D> s(L, 3);
 
   Rng r;
 
-  double z = exp(1 / 0.15);
+  for (double μ = μmin; μ <= μmax; μ += dμ) {
+    double z = exp(μ);
+    for (unsigned i = 0; i < 1e4; i++) {
+      for (unsigned j = 0; j < s.size(); j++) {
+        s.stepGrandCanonical(z, r);
+//        s.randomEventChain(r);
+        s.randomExchange(r);
+      }
 
-  if (!s.compatible()) {
-    std::cerr << "Storted incompatible!" << std::endl;
-    return 1;
-  }
+      if (!s.compatible()) {
+        std::cerr << "Not compatible!" << std::endl;
+        return 1;
+      }
 
-  while (s.density() < 0.56) {
-    s.sweepGrandCanonical(z, r);
-  }
+    }
 
-  if (!s.compatible()) {
-    std::cerr << "Not compatible!" << std::endl;
-    return 1;
+    std::cout << μ << " " << s.density() << std::endl;
   }
 
+
+  /*
   std::cerr << "Found state with appropriate density." << std::endl;
 
   System<D> s0 = s;
@@ -586,22 +667,19 @@ int main() {
   auto start = std::chrono::high_resolution_clock::now();
   while (nC / s.size() < 1e5) {
     if (n < 20 * log(i + 1)) {
+      auto stop = std::chrono::high_resolution_clock::now();
+      auto duration = duration_cast<std::chrono::microseconds>(stop - start);
       n++;
-      std::cout << nC / s.size() << " " << s.overlap(s0) << std::endl;
+      std::cout << duration.count() << " " << s.overlap(s0) << std::endl;
     }
-    //unsigned nn = s.flipCluster(Transformation<D>(L, ms, r), r.pick(s.vertices));
-    //nC += nn;
+//    unsigned nn = s.flipCluster(Transformation<D>(L, ms, r), r.pick(s.vertices));
+//    nC += nn;
 //    clusterDist[nn]++;
     //s.sweepLocal(r);
     nC += s.size();
     s.sweepSwap(r);
     i++;
   }
-  auto stop = std::chrono::high_resolution_clock::now();
-
-  auto duration = duration_cast<std::chrono::microseconds>(stop - start);
-
-  std::cerr << duration.count() << std::endl;
 
   if (!s.compatible()) {
     std::cerr << "Not compatible!" << std::endl;
@@ -618,6 +696,7 @@ int main() {
     file << i << " ";
   }
   file.close();
+  */
 
   return 0;
 }
diff --git a/ciamarra.cpp b/ciamarra.cpp
index 6b932bf..f2f44c1 100644
--- a/ciamarra.cpp
+++ b/ciamarra.cpp
@@ -1,10 +1,477 @@
+#include <eigen3/Eigen/Dense>
+#include <eigen3/Eigen/src/Core/CwiseNullaryOp.h>
 #include <iostream>
+#include <list>
+#include <vector>
+#include <queue>
 
-#include "glass.hpp"
+#include "pcg-cpp/include/pcg_random.hpp"
+#include "randutils/randutils.hpp"
 
+using Rng = randutils::random_generator<pcg32>;
 
-void print(const CiamarraSystem<2>& s) {
-  for (const Vertex<2, CiamarraState<2>>& v : s.vertices) {
+template <int D> using Vector = Eigen::Matrix<int, D, 1>;
+template <int D> using Matrix = Eigen::Matrix<int, D, D>;
+
+int iPow(int x, unsigned p) {
+  if (p == 0)
+    return 1;
+  if (p == 1)
+    return x;
+
+  int tmp = iPow(x, p / 2);
+  if (p % 2 == 0)
+    return tmp * tmp;
+  else
+    return x * tmp * tmp;
+}
+
+unsigned mod(signed a, unsigned b) {
+  return ((a < 0) ? (a + (1 - a / (signed)b) * b) : a) % b;
+}
+
+template <int D, typename Derived> Vector<D> mod(const Eigen::MatrixBase<Derived>& v, unsigned b) {
+  Vector<D> u;
+  for (unsigned i = 0; i < D; i++) {
+    u(i) = mod(v(i), b);
+  }
+  return u;
+}
+
+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 <unsigned D> std::vector<Matrix<D>> generateTorusMatrices() {
+  std::vector<Matrix<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_back(); // don't want the identity matrix!
+
+  for (std::array<double, D> sequence : sequences) {
+    Matrix<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<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 if (k == l && (k != i && k != j)) {
+              m(k, l) = 1;
+            } else {
+              m(k, l) = 0;
+            }
+          }
+        }
+        mats.push_back(m);
+      }
+    }
+  }
+
+  return mats;
+}
+
+template <unsigned D> class State : public Vector<D> {
+public:
+  State() : Vector<D>(Vector<D>::Zero()) {}
+
+  State(const Vector<D>& v) : Vector<D>(v) {}
+
+  State(unsigned a, signed b) : Vector<D>(Vector<D>::Zero()) { this->operator()(a) = b; }
+
+  State(Rng& r) : State(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})) {}
+
+  bool isEmpty() const { return this->squaredNorm() == 0; }
+
+  void remove() { this->setZero(); }
+
+  State<D> flip() const {
+    State<D> s;
+    for (unsigned i = 0; i < D; i++) {
+      s(i) = -this->operator()(i);
+    }
+    return s;
+  }
+};
+
+template <unsigned D> class Transformation {
+  public:
+    unsigned L;
+    Matrix<D> m;
+    Vector<D> v;
+
+    Transformation(unsigned L) : L(L) {
+      m.setIdentity();
+      v.setZero();
+    }
+
+    Transformation(unsigned L, const Matrix<D>& m, const Vector<D>& v) : L(L), m(m), v(v) {}
+
+    Transformation(unsigned L, const std::vector<Matrix<D>>& ms, Rng& r) : m(r.pick(ms)), L(L) {
+      for (unsigned i = 0; i < D; i++) {
+        v[i] = r.uniform((unsigned)0, L - 1);
+      }
+
+      v = v - m * v;
+    }
+
+    Vector<D> apply(const Vector<D>& x) const {
+      return mod<D>(v + m * x, L);
+    }
+
+    Transformation<D> apply(const Transformation<D>& t) const {
+      Transformation<D> tNew(L);
+
+      tNew.m = m * t.m;
+      tNew.v = apply(t.v);
+
+      return tNew;
+    }
+
+    State<D> apply(const State<D>& x) const {
+      return State<D>(m * x);
+    }
+
+    Transformation<D> inverse() const {
+      return Transformation<D>(L, m.transpose(), -m.transpose() * v);
+    }
+};
+
+template <unsigned D> class HalfEdge;
+
+template <unsigned D> class Vertex {
+public:
+  Vector<D> position;
+  State<D> state;
+  std::vector<HalfEdge<D>> adjacentEdges;
+  bool marked;
+
+  bool isEmpty() const { return state.isEmpty(); }
+};
+
+template <unsigned D> class HalfEdge {
+public:
+  Vertex<D>& neighbor;
+  Vector<D> Δx;
+
+  HalfEdge(Vertex<D>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
+};
+
+template <unsigned D> class System {
+public:
+  const unsigned L;
+  unsigned N;
+  std::vector<Vertex<D>> vertices;
+  Transformation<D> orientation;
+
+  unsigned vectorToIndex(const Vector<D>& x) const {
+    unsigned i = 0;
+    for (unsigned d = 0; d < D; d++) {
+      i += x[d] * iPow(L, d);
+    }
+    return i;
+  }
+
+  Vector<D> indexToVector(unsigned i) const {
+    Vector<D> x;
+    for (unsigned d = 0; d < D; d++) {
+      x[d] = (i / iPow(L, d)) % L;
+    }
+    return x;
+  }
+
+  System(unsigned L) : L(L), N(0), vertices(iPow(L, D)), orientation(L) {
+    for (unsigned i = 0; i < iPow(L, D); i++) {
+      vertices[i].position = indexToVector(i);
+      vertices[i].adjacentEdges.reserve(2 * D);
+      vertices[i].marked = false;
+    }
+
+    for (unsigned d = 0; d < D; d++) {
+      Vector<D> Δx = Vector<D>::Zero();
+      Δx[d] = 1;
+      for (signed i = 0; i < iPow(L, D); i++) {
+        unsigned j = iPow(L, d + 1) * (i / iPow(L, d + 1)) + mod(i + iPow(L, d), pow(L, d + 1));
+        vertices[i].adjacentEdges.push_back(HalfEdge<D>(vertices[j], Δx));
+        vertices[j].adjacentEdges.push_back(HalfEdge<D>(vertices[i], -Δx));
+      }
+    }
+  }
+
+  std::list<std::reference_wrapper<Vertex<D>>> overlaps(Vertex<D>& v, const State<D>& s,
+                                                        bool excludeSelf = false) {
+    std::list<std::reference_wrapper<Vertex<D>>> o;
+
+    if (s.isEmpty()) {
+      return o;
+    }
+
+    if (!v.isEmpty() && !excludeSelf) {
+      o.push_back(v);
+    }
+
+    for (const HalfEdge<D>& e : v.adjacentEdges) {
+      if (!e.neighbor.isEmpty()) {
+        if (s.dot(e.Δx) == 1 || e.neighbor.state.dot(e.Δx) == -1) {
+          o.push_back(e.neighbor);
+        }
+      }
+    }
+
+    return o;
+  }
+
+  bool insert(Vertex<D>& v, const State<D>& s) {
+    if (overlaps(v, s).empty()) {
+      v.state = s;
+      N++;
+      return true;
+    } else {
+      return false;
+    }
+  }
+
+  bool tryDeletion(Vertex<D>& v) {
+    if (v.isEmpty()) {
+      return false;
+    } else {
+      v.state.remove();
+      N--;
+      return true;
+    }
+  }
+
+  bool tryRandomMove(Rng& r) {
+    Vertex<D>& v = r.pick(vertices);
+    State<D> oldState = v.state;
+
+    if (!tryDeletion(v)) {
+      return false;
+    }
+
+    if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) {
+      for (HalfEdge<D>& e : v.adjacentEdges) {
+        if (1 == e.Δx.dot(oldState)) {
+          if (insert(e.neighbor, oldState.flip())) {
+            return true;
+          }
+          break;
+        }
+      }
+    } else {
+      State<D> newState(r);
+      while (newState.dot(oldState) == 1) {
+        newState = State<D>(r);
+      }
+      if (insert(v, newState)) {
+        return true;
+      }
+    }
+    v.state = oldState;
+    N++;
+    return false;
+  }
+
+  bool trySwap(Vertex<D>& v1, Vertex<D>& v2) {
+    if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) {
+      std::swap(v1.state, v2.state);
+      return true;
+    } else {
+      return false;
+    }
+  }
+
+  bool tryRandomSwap(Rng& r) {
+    Vertex<D>& v1 = r.pick(vertices);
+    Vertex<D>& v2 = r.pick(vertices);
+
+    return trySwap(v1, v2);
+  }
+
+  void setGroundState() {
+    N = 0;
+
+    for (Vertex<D>& v : vertices) {
+      unsigned a = 0;
+      for (unsigned d = 0; d < D; d++) {
+        a += (d + 1) * v.position(d);
+      }
+      a %= 2 * D + 1;
+
+      v.state.setZero() = Vector<D>::Zero();
+
+      if (0 < a && a <= D) {
+        v.state(a - 1) = -1;
+        N++;
+      } else if (D < a) {
+        v.state(2 * D - a) = 1;
+        N++;
+      }
+    }
+  }
+
+  bool compatible() {
+    for (Vertex<D>& v : vertices) {
+      if (overlaps(v, v.state, true).size() > 0) {
+        return false;
+      }
+    }
+
+    return true;
+  }
+
+  double density() const { return N / pow(L, D); }
+
+  unsigned maxOccupation() const {
+    //      return (2 * D * iPow(L, D)) / (2 * D + 1);
+    return iPow(L, D);
+  }
+
+  void sweepGrandCanonical(double z, Rng& r) {
+    for (unsigned i = 0; i < iPow(L, D); i++) {
+      if (0.5 < r.uniform(0.0, 1.0)) {
+        double pIns = maxOccupation() * z / (N + 1);
+
+        if (pIns > r.uniform(0.0, 1.0)) {
+          while (true) {
+            Vertex<D>& v = r.pick(vertices);
+            if (v.isEmpty()) {
+              insert(v, State<D>(r));
+              break;
+            }
+          }
+        }
+      } else {
+
+        double pDel = N / (z * maxOccupation());
+
+        if (pDel > r.uniform(0.0, 1.0)) {
+          tryDeletion(r.pick(vertices));
+        }
+      }
+
+      tryRandomMove(r);
+    }
+  }
+
+  void sweepLocal(Rng& r) {
+    for (unsigned i = 0; i < iPow(L, D); i++) {
+      tryRandomMove(r);
+    }
+  }
+
+  void sweepSwap(Rng& r) {
+    for (unsigned i = 0; i < iPow(L, D); i++) {
+      tryRandomSwap(r);
+    }
+  }
+
+  unsigned flipCluster(const Transformation<D>& R, Vertex<D>& v0, Rng& r, bool dry = false) {
+    std::queue<std::reference_wrapper<Vertex<D>>> q;
+    q.push(v0);
+
+    unsigned n = 0;
+
+    while (!q.empty()) {
+      Vertex<D>& v = q.front();
+      q.pop();
+
+      if (!v.marked) {
+        Vector<D> xNew = R.apply(v.position);
+        Vertex<D>& vNew = vertices[vectorToIndex(xNew)];
+
+        v.marked = true;
+        vNew.marked = true;
+
+        State<D> s = R.apply(v.state);
+        State<D> sNew = R.apply(vNew.state);
+
+        std::list<std::reference_wrapper<Vertex<D>>> overlaps1 = overlaps(vNew, s, true);
+        std::list<std::reference_wrapper<Vertex<D>>> overlaps2 = overlaps(v, sNew, true);
+        overlaps1.splice(overlaps1.begin(), overlaps2);
+
+        for (Vertex<D>& vn : overlaps1) {
+          if (!vn.marked) {
+            q.push(vn);
+          }
+        }
+
+        if (!dry) {
+          v.state = sNew;
+          vNew.state = s;
+        }
+
+        n += 1;
+      }
+    }
+
+    return n;
+  }
+
+  void swendsenWang(const Transformation<D>& R, Rng& r) {
+    for (Vertex<D>& v : vertices) {
+      if (!v.marked) {
+        bool dry = 0.5 < r.uniform(0.0, 1.0);
+        unsigned n = flipCluster(R, v, r, dry);
+        if (n > pow(L, D) / 4 && !dry) {
+          orientation = R.apply(orientation);
+        }
+      }
+    }
+
+    for (Vertex<D>& v : vertices) {
+      v.marked = false;
+    }
+  }
+
+  int overlap(const System<D>& s) const {
+    int o = 0;
+
+    for (unsigned i = 0; i < vertices.size(); i++) {
+      State<D> s2 = orientation.apply(s.vertices[vectorToIndex(orientation.inverse().apply(indexToVector(i)))].state);
+      o += vertices[i].state.dot(s2);
+    }
+
+    return o;
+  }
+};
+
+void print(const System<2>& s) {
+  for (const Vertex<2>& v : s.vertices) {
     if (v.state(0) == 1 && v.state(1) == 0) {
       std::cerr << "▶";
     } else if (v.state(0) == -1 && v.state(1) == 0) {
@@ -42,7 +509,7 @@ int main() {
   double Tmax = 0.2;
   double δT = 0.02;
 
-  CiamarraSystem<D> s(L);
+  System<D> s(L);
 
   Rng r;
 
@@ -58,11 +525,8 @@ int main() {
   } 
 
   std::cerr << "Found state with appropriate density." << std::endl;
-    if (!s.compatible()) {
-      std::cerr <<"whoops!" << std::endl;
-    }
 
-  CiamarraSystem<D> s0 = s;
+  System<D> s0 = s;
 
   std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
 
@@ -72,9 +536,16 @@ int main() {
       n++;
       std::cout << i << " " << s.overlap(s0) << std::endl;
     }
+    Matrix<D> m = r.pick(ms);
+    Vertex<D>& v = r.pick(s.vertices);
+    unsigned nC = s.flipCluster(Transformation<D>(L, m, v.position - m * v.position), v, r);
+    std::cerr << nC << std::endl;
+    for (Vertex<D>& v : s.vertices) {
+      v.marked = false;
+    }
 //    s.sweepLocal(r);
 //    s.sweepSwap(r);
-    s.swendsenWang(Transformation<D>(L, ms, r), r);
+//    s.swendsenWang(Transformation<D>(L, ms, r), r);
   }
 
   return 0;
diff --git a/glass.hpp b/glass.hpp
deleted file mode 100644
index b35964b..0000000
--- a/glass.hpp
+++ /dev/null
@@ -1,677 +0,0 @@
-#include <limits>
-#include <list>
-#include <queue>
-#include <vector>
-
-#include <eigen3/Eigen/Dense>
-
-#include "pcg-cpp/include/pcg_random.hpp"
-#include "randutils/randutils.hpp"
-
-using Rng = randutils::random_generator<pcg32>;
-
-template <int D> using Vector = Eigen::Matrix<int, D, 1>;
-template <int D> using Matrix = Eigen::Matrix<int, D, D>;
-
-int iPow(int x, unsigned p) {
-  if (p == 0)
-    return 1;
-  if (p == 1)
-    return x;
-
-  int tmp = iPow(x, p / 2);
-  if (p % 2 == 0)
-    return tmp * tmp;
-  else
-    return x * tmp * tmp;
-}
-
-unsigned mod(signed a, unsigned b) { return ((a < 0) ? (a + (1 - a / (signed)b) * b) : a) % b; }
-
-template <int D, typename Derived> Vector<D> mod(const Eigen::MatrixBase<Derived>& v, unsigned b) {
-  Vector<D> u;
-  for (unsigned i = 0; i < D; i++) {
-    u(i) = mod(v(i), b);
-  }
-  return u;
-}
-
-template <int D> void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
-  if (level > 0) {
-    unsigned new_level = level - 1;
-    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 <int D> std::vector<Matrix<D>> generateTorusMatrices() {
-  std::vector<Matrix<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_back(); // don't want the identity matrix!
-
-  for (std::array<double, D> sequence : sequences) {
-    Matrix<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<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 if (k == l && (k != i && k != j)) {
-              m(k, l) = 1;
-            } else {
-              m(k, l) = 0;
-            }
-          }
-        }
-        mats.push_back(m);
-      }
-    }
-  }
-
-  return mats;
-}
-
-template <int D> class CiamarraState : public Vector<D> {
-public:
-  CiamarraState() : Vector<D>(Vector<D>::Zero()) {}
-
-  CiamarraState(const Vector<D>& v) : Vector<D>(v) {}
-
-  CiamarraState(unsigned a, signed b) : Vector<D>(Vector<D>::Zero()) { this->operator()(a) = b; }
-
-  CiamarraState(Rng& r) : CiamarraState(r.uniform(0, D - 1), r.pick({-1, 1})) {}
-
-  bool isEmpty() const { return this->squaredNorm() == 0; }
-
-  void remove() { this->setZero(); }
-
-  CiamarraState<D> flip() const {
-    CiamarraState<D> s;
-    for (unsigned i = 0; i < D; i++) {
-      s(i) = -this->operator()(i);
-    }
-    return s;
-  }
-};
-
-class BiroliState {
-public:
-  unsigned type;
-  unsigned occupiedNeighbors;
-
-  BiroliState() {
-    type = std::numeric_limits<unsigned>::max();
-    occupiedNeighbors = 0;
-  }
-
-  BiroliState(unsigned t, unsigned o) {
-    type = t;
-    occupiedNeighbors = o;
-  }
-
-  BiroliState(Rng& r) {
-    type = r.pick({1,2,2,2,2,2,3,3,3,3});
-  }
-
-  bool isEmpty() const { return type == std::numeric_limits<unsigned>::max(); }
-
-  void remove() { type = std::numeric_limits<unsigned>::max(); };
-};
-
-template <int D> class Transformation {
-public:
-  unsigned L;
-  Matrix<D> m;
-  Vector<D> v;
-
-  Transformation(unsigned L) : L(L) {
-    m.setIdentity();
-    v.setZero();
-  }
-
-  Transformation(unsigned L, const Matrix<D>& m, const Vector<D>& v) : L(L), m(m), v(v) {}
-
-  Transformation(unsigned L, const std::vector<Matrix<D>>& ms, Rng& r) : m(r.pick(ms)), L(L) {
-    for (unsigned i = 0; i < D; i++) {
-      v[i] = r.uniform((unsigned)0, L - 1);
-    }
-
-    v = v - m * v;
-  }
-
-  Transformation<D> inverse() const {
-    return Transformation<D>(L, m.transpose(), -m.transpose() * v);
-  }
-
-  Vector<D> apply(const Vector<D>& x) const { return mod<D>(v + m * x, L); }
-
-  Transformation<D> apply(const Transformation<D>& t) const {
-    Transformation<D> tNew(L);
-
-    tNew.m = m * t.m;
-    tNew.v = apply(t.v);
-
-    return tNew;
-  }
-
-  CiamarraState<D> apply(const CiamarraState<D>& s) const { return CiamarraState<D>(m * s); }
-
-  BiroliState apply(const BiroliState& s) const { return s; }
-};
-
-template <int D, typename State> class HalfEdge;
-
-template <int D, typename State> class Vertex {
-public:
-  Vector<D> position;
-  std::vector<HalfEdge<D, State>> adjacentEdges;
-  State state;
-  bool marked;
-
-  bool isEmpty() const { return state.isEmpty(); }
-};
-
-template <int D, class State> class HalfEdge {
-public:
-  Vertex<D, State>& neighbor;
-  Vector<D> Δx;
-
-  HalfEdge(Vertex<D, State>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
-};
-
-template <int D, class State> class System {
-public:
-  const unsigned L;
-  unsigned N;
-  std::vector<Vertex<D, State>> vertices;
-  Transformation<D> orientation;
-
-  unsigned vectorToIndex(const Vector<D>& x) const {
-    unsigned i = 0;
-    for (unsigned d = 0; d < D; d++) {
-      i += x[d] * iPow(L, d);
-    }
-    return i;
-  }
-
-  Vector<D> indexToVector(unsigned i) const {
-    Vector<D> x;
-    for (unsigned d = 0; d < D; d++) {
-      x[d] = (i / iPow(L, d)) % L;
-    }
-    return x;
-  }
-
-  System(unsigned L) : L(L), N(0), vertices(iPow(L, D)), orientation(L) {
-    for (unsigned i = 0; i < iPow(L, D); i++) {
-      vertices[i].position = indexToVector(i);
-      vertices[i].adjacentEdges.reserve(2 * D);
-      vertices[i].marked = false;
-    }
-
-    for (unsigned d = 0; d < D; d++) {
-      Vector<D> Δx = Vector<D>::Zero();
-      Δx[d] = 1;
-      for (signed i = 0; i < iPow(L, D); i++) {
-        unsigned j = iPow(L, d + 1) * (i / iPow(L, d + 1)) + mod(i + iPow(L, d), pow(L, d + 1));
-        vertices[i].adjacentEdges.push_back(HalfEdge<D, State>(vertices[j], Δx));
-        vertices[j].adjacentEdges.push_back(HalfEdge<D, State>(vertices[i], -Δx));
-      }
-    }
-  }
-
-  virtual std::list<std::reference_wrapper<Vertex<D, State>>> overlaps(Vertex<D, State>&) { return {}; };
-
-  virtual bool insert(Vertex<D, State>& v, const State& s, bool force = false) { return false; };
-
-  virtual bool remove(Vertex<D, State>& v) { return false;};
-
-  virtual bool randomMove(Rng& r) { return false;};
-
-  virtual void swap(Vertex<D, State>& v1, Vertex<D, State>& v2) {};
-
-  virtual bool exchange(Vertex<D, State>& v1, Vertex<D, State>& v2) { return false;};
-
-  bool randomExchange(Rng& r) {
-    Vertex<D, State>& v1 = r.pick(vertices);
-    Vertex<D, State>& v2 = r.pick(vertices);
-
-    return exchange(v1, v2);
-  }
-
-  bool compatible() {
-    for (Vertex<D, State>& v : vertices) {
-      if (overlaps(v).size() > 0) {
-        return false;
-      }
-    }
-
-    return true;
-  }
-
-  double density() const { return N / pow(L, D); }
-
-  unsigned size() const { return vertices.size(); }
-
-  void sweepGrandCanonical(double z, Rng& r) {
-    for (unsigned i = 0; i < iPow(L, D); i++) {
-      if (0.5 < r.uniform(0.0, 1.0)) {
-        double pIns = size() * z / (N + 1);
-
-        if (pIns > r.uniform(0.0, 1.0)) {
-          while (true) {
-            Vertex<D, State>& v = r.pick(vertices);
-            if (v.isEmpty()) {
-              insert(v, State(r));
-              break;
-            }
-          }
-        }
-      } else {
-
-        double pDel = N / (z * size());
-
-        if (pDel > r.uniform(0.0, 1.0)) {
-          remove(r.pick(vertices));
-        }
-      }
-
-      randomMove(r);
-    }
-  }
-
-  void sweepLocal(Rng& r) {
-    for (unsigned i = 0; i < iPow(L, D); i++) {
-      randomMove(r);
-    }
-  }
-
-  void sweepSwap(Rng& r) {
-    for (unsigned i = 0; i < iPow(L, D); i++) {
-      randomExchange(r);
-    }
-  }
-
-  unsigned flipCluster(const Transformation<D>& R, Vertex<D, State>& v0, bool dry = false) {
-    unsigned n = 0;
-
-    Vertex<D, State>& v0New = vertices[vectorToIndex(R.apply(v0.position))];
-
-    if (&v0New != &v0) {
-      std::queue<std::reference_wrapper<Vertex<D, State>>> q;
-
-      v0.marked = true;
-      v0New.marked = true;
-
-      swap(v0, v0New);
-
-      for (Vertex<D, State>& vn : overlaps(v0)) {
-        if (!vn.marked) {
-          q.push(vn);
-        }
-      }
-      for (Vertex<D, State>& vn : overlaps(v0New)) {
-        if (!vn.marked) {
-          q.push(vn);
-        }
-      }
-
-      while (!q.empty()) {
-        Vertex<D, State>& v = q.front();
-        q.pop();
-
-        if (!v.marked && !overlaps(v).empty()) {
-          v.marked = true;
-          Vertex<D, State>& vNew = vertices[vectorToIndex(R.apply(v.position))];
-
-          if (&vNew != &v) {
-            vNew.marked = true;
-
-            swap(v, vNew);
-
-            for (Vertex<D, State>& vn : overlaps(v)) {
-              if (!vn.marked) {
-                q.push(vn);
-              }
-            }
-            for (Vertex<D, State>& vn : overlaps(vNew)) {
-              if (!vn.marked) {
-                q.push(vn);
-              }
-            }
-
-            n += 2;
-          } else {
-            n += 1;
-          }
-        }
-        if (q.empty()) {
-          for (Vertex<D, State>& vv : vertices) {
-            if (!vv.marked && !overlaps(vv).empty()) {
-//              std::cerr << "Found bad state at end" << std::endl;
-              q.push(vv);
-            }
-          }
-        }
-      }
-
-
-    }
-
-    if (n > size() / 4) {
-      orientation = R.apply(orientation);
-    }
-
-    for (Vertex<D, State>& v : vertices) {
-      v.marked = false;
-    }
-
-    return n;
-  }
-
-  void swendsenWang(const Transformation<D>& R, Rng& r) {
-    for (Vertex<D, State>& v : vertices) {
-      if (!v.marked) {
-        bool dry = 0.5 < r.uniform(0.0, 1.0);
-        unsigned n = flipCluster(R, v, dry);
-        if (n > size() / 4 && !dry) {
-          orientation = R.apply(orientation);
-        }
-      }
-    }
-
-    for (Vertex<D, State>& v : vertices) {
-      v.marked = false;
-    }
-  }
-
-  virtual int overlap(const System<D, State>& s) const { return 0; };
-};
-
-template <int D> class CiamarraSystem : public System<D, CiamarraState<D>> {
-  public:
-
-    using System<D, CiamarraState<D>>::System;
-  std::list<std::reference_wrapper<Vertex<D, CiamarraState<D>>>>
-  overlaps(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s,
-           bool excludeSelf = false) override {
-    std::list<std::reference_wrapper<Vertex<D, CiamarraState<D>>>> o;
-
-    if (s.isEmpty()) {
-      return o;
-    }
-
-    if (!v.isEmpty() && !excludeSelf) {
-      o.push_back(v);
-    }
-
-    for (const HalfEdge<D, CiamarraState<D>>& e : v.adjacentEdges) {
-      if (!e.neighbor.isEmpty()) {
-        if (s.dot(e.Δx) == 1 || e.neighbor.state.dot(e.Δx) == -1) {
-          o.push_back(e.neighbor);
-        }
-      }
-    }
-
-    return o;
-  }
-
-  bool insert(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s, bool force = false) override {
-    if (force || overlaps(v, s).empty()) {
-      v.state = s;
-      this->N++;
-      return true;
-    } else {
-      return false;
-    }
-  }
-
-  bool remove(Vertex<D, CiamarraState<D>>& v) override {
-    if (v.isEmpty()) {
-      return false;
-    } else {
-      v.state.remove();
-      this->N--;
-      return true;
-    }
-  }
-
-  bool randomMove(Rng& r) override {
-    Vertex<D, CiamarraState<D>>& v = r.pick(this->vertices);
-    CiamarraState<D> oldState = v.state;
-
-    if (!remove(v)) {
-      return false;
-    }
-
-    if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) {
-      for (HalfEdge<D, CiamarraState<D>>& e : v.adjacentEdges) {
-        if (1 == e.Δx.dot(oldState)) {
-          if (insert(e.neighbor, oldState.flip())) {
-            return true;
-          }
-          break;
-        }
-      }
-    } else {
-      CiamarraState<D> newState(r);
-      while (newState.dot(oldState) == 1) {
-        newState = CiamarraState<D>(r);
-      }
-      if (insert(v, newState)) {
-        return true;
-      }
-    }
-    v.state = oldState;
-    this->N++;
-    return false;
-  }
-
-  void swap(Vertex<D, CiamarraState<D>>& v1, Vertex<D, CiamarraState<D>>& v2) override {
-    std::swap(v1.state, v2.state);
-  }
-
-  bool exchange(Vertex<D, CiamarraState<D>>& v1, Vertex<D, CiamarraState<D>>& v2) override {
-    if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) {
-      swap(v1, v2);
-      return true;
-    } else {
-      return false;
-    }
-  }
-
-  void setGroundState() {
-    this->N = 0;
-
-    for (Vertex<D, CiamarraState<D>>& v : this->vertices) {
-      unsigned a = 0;
-      for (unsigned d = 0; d < D; d++) {
-        a += (d + 1) * v.position(d);
-      }
-      a %= 2 * D + 1;
-
-      v.state.setZero() = Vector<D>::Zero();
-
-      if (0 < a && a <= D) {
-        v.state(a - 1) = -1;
-        this->N++;
-      } else if (D < a) {
-        v.state(2 * D - a) = 1;
-        this->N++;
-      }
-    }
-  }
-
-  int overlap(const System<D, CiamarraState<D>>& s) const override {
-    int o = 0;
-
-    for (unsigned i = 0; i < this->vertices.size(); i++) {
-      CiamarraState<D> s2 = this->orientation.apply(
-          s.vertices[this->vectorToIndex(this->orientation.inverse().apply(this->indexToVector(i)))].state);
-      o += this->vertices[i].state.dot(s2);
-    }
-
-    return o;
-  }
-};
-
-template <int D> class BiroliSystem : public System<D, BiroliState> {
-public:
-  using System<D, BiroliState>::System;
-
-  bool canReplaceWith(const Vertex<D, BiroliState>& v, const BiroliState& s) {
-    if (s.isEmpty()) {
-      return true;
-    }
-    if (s.type < v.state.occupiedNeighbors) {
-      return false;
-    } else {
-      int Δn = 0;
-      if (v.isEmpty()) {
-        Δn += 1;
-      }
-      for (const HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
-        if (e.neighbor.state.type < e.neighbor.state.occupiedNeighbors + Δn) {
-          return false;
-        }
-      }
-    }
-    return true;
-  }
-
-  std::list<std::reference_wrapper<Vertex<D, BiroliState>>>
-  overlaps(Vertex<D, BiroliState>& v) override {
-    std::list<std::reference_wrapper<Vertex<D, BiroliState>>> o;
-
-    if (v.isEmpty()) { // an empty site cannot be frustrating anyone
-      return o;
-    }
-
-    bool selfFrustrated = v.state.occupiedNeighbors > v.state.type;
-    bool anyNeighborFrustrated = false;
-
-    for (const HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
-      if (!e.neighbor.isEmpty()) {
-        bool thisNeighborFrustrated = e.neighbor.state.type < e.neighbor.state.occupiedNeighbors;
-
-        if (thisNeighborFrustrated) {
-          anyNeighborFrustrated = true;
-        }
-
-        if (selfFrustrated || thisNeighborFrustrated) {
-          o.push_back(e.neighbor);
-        }
-      }
-    }
-
-    if (selfFrustrated || anyNeighborFrustrated) {
-      o.push_back(v);
-    }
-
-    return o;
-  }
-
-  bool insert(Vertex<D, BiroliState>& v, const BiroliState& s, bool force = false) override {
-    if (force || (v.isEmpty() && canReplaceWith(v, s))) {
-      v.state.type = s.type;
-      for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
-        e.neighbor.state.occupiedNeighbors++;
-      }
-      this->N++;
-      return true;
-    } else {
-      return false;
-    }
-  }
-
-  bool remove(Vertex<D, BiroliState>& v) override {
-    if (v.isEmpty()) {
-      return false;
-    } else {
-      v.state.remove();
-      for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
-        e.neighbor.state.occupiedNeighbors--;
-      }
-      this->N--;
-      return true;
-    }
-  }
-
-  bool randomMove(Rng& r) override {
-    Vertex<D, BiroliState>& v = r.pick(this->vertices);
-
-    return exchange(v, r.pick(v.adjacentEdges).neighbor);
-  }
-
-  void swap(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) override {
-    if (v1.state.type != v2.state.type) {
-      if (!v1.isEmpty() && !v2.isEmpty()) {
-        std::swap(v1.state.type, v2.state.type);
-      } else if (v1.isEmpty()) {
-        unsigned t = v2.state.type;
-        remove(v2);
-        insert(v1, BiroliState(t, 0), true);
-      } else {
-        unsigned t = v1.state.type;
-        remove(v1);
-        insert(v2, BiroliState(t, 0), true);
-      }
-    }
-  }
-
-  bool exchange(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) override {
-    if (canReplaceWith(v1, v2.state) && canReplaceWith(v2, v1.state)) {
-      swap(v1, v2);
-      return true;
-    } else {
-      return false;
-    }
-  }
-
-  int overlap(const System<D, BiroliState>& s) const override {
-    int o = 0;
-
-    for (unsigned i = 0; i < this->vertices.size(); i++) {
-      BiroliState s2 = this->orientation.apply(
-          s.vertices[this->vectorToIndex(this->orientation.inverse().apply(this->indexToVector(i)))].state);
-      if (s2.type > 0 && s2.type == this->vertices[i].state.type) {
-        o++;
-      } else {
-        o--;
-      }
-    }
-
-    return o;
-  }
-};