From 1e49049787ca9c2c138f28a611dfe7496fcbb0a0 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Sun, 21 Mar 2021 11:20:27 +0100
Subject: Got Swendsen-Wang working.

---
 glass.cpp | 570 ++++++++++++++++++++++++++++++++++++++++----------------------
 1 file changed, 368 insertions(+), 202 deletions(-)

diff --git a/glass.cpp b/glass.cpp
index e8739d0..ba983d2 100644
--- a/glass.cpp
+++ b/glass.cpp
@@ -1,294 +1,447 @@
+#include <eigen3/Eigen/Dense>
 #include <iostream>
-#include <limits>
-#include <vector>
 #include <list>
-#include <set>
-#include <eigen3/Eigen/Dense>
-#include <unordered_map>
+#include <vector>
+#include <queue>
 
 #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>;
+template <int D> using Vector = Eigen::Matrix<int, D, 1>;
+template <int D> using Matrix = Eigen::Matrix<int, D, D>;
 
-// Integer power function
 int iPow(int x, unsigned p) {
-  if (p == 0) return 1;
-  if (p == 1) return x;
+  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;
+  if (p % 2 == 0)
+    return tmp * tmp;
+  else
+    return x * tmp * tmp;
 }
 
 unsigned mod(signed a, unsigned b) {
-  return ((a %= b) < 0) ? a + b : a;
+  return ((a < 0) ? (a + (1 - a / (signed)b) * b) : a) % b;
 }
 
-template<unsigned D> class State {
-  public:
-    Vector<D> σ;
+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;
+}
 
-    State(unsigned a, signed b) : σ(Vector<D>::Zero()) {
-      σ[a] = b;
+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);
+  }
+}
 
-    State(Rng& r) : State(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})) {}
+template <unsigned D> std::vector<Matrix<D>> generateTorusMatrices() {
+  std::vector<Matrix<D>> mats;
 
-    State<D> apply(const Matrix<D>& m) const {
-      return {m * σ};
+  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;
+        }
+      }
     }
-};
 
-template <unsigned D> class Vertex;
+    mats.push_back(m);
+  }
 
-template<unsigned D> class Particle {
-  public:
-    Vertex<D>& position;
-    State<D> state;
-    bool marked;
+  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);
+      }
+    }
+  }
 
-    Particle(Vertex<D>& p, const State<D>& s) : position(p), state(s) {}
-};
+  return mats;
+}
 
-template<unsigned D> class HalfEdge {
-  public:
-    Vertex<D>& neighbor;
-    Vector<D> Δx;
+template <unsigned D> class State : public Vector<D> {
+public:
+  State() : Vector<D>(Vector<D>::Zero()) {}
 
-    HalfEdge(Vertex<D>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
-};
+  State(const Vector<D>& v) : Vector<D>(v) {}
 
-template<unsigned D> class Vertex {
-  public:
-    Vector<D> position;
-    std::vector<HalfEdge<D>> adjacentEdges;
-    typename std::unordered_map<unsigned, Particle<D>>::iterator occupant;
+  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(); }
 };
 
-template<unsigned D> class System {
+template <unsigned D> class Transformation {
   public:
-    const unsigned L;
-    std::vector<Vertex<D>> vertices;
-    std::unordered_map<unsigned, Particle<D>> particles;
+    unsigned L;
+    Matrix<D> m;
+    Vector<D> v;
 
-    unsigned vectorToIndex(const Vector<D>& x) const {
-      unsigned i = 0;
-      for (unsigned d = 0; d < D; d++) {
-        i += x[d] * iPow(L, d);
+    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);
       }
-      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;
+    Vector<D> apply(const Vector<D>& x) const {
+      return mod<D>(v + m * (x - v), L);
     }
 
-    bool isEmpty(const Vertex<D>& v) const {
-      return v.occupant == particles.end();
+    State<D> apply(const State<D>& x) const {
+      return State<D>(m * x);
     }
+};
+
+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;
 
-    unsigned N() const {
-      return particles.size();
+  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;
+  }
 
-    System(unsigned L) : L(L), vertices(iPow(L, D)) {
-      for (unsigned i = 0; i < iPow(L, D); i++) {
-        vertices[i].position = indexToVector(i);
-        vertices[i].adjacentEdges.reserve(2 * D);
-        vertices[i].occupant = particles.end();
-      }
+  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;
+  }
 
-      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));
-        }
+  System(unsigned L) : L(L), N(0), vertices(iPow(L, D)) {
+    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<typename std::unordered_map<unsigned, Particle<D>>::iterator> overlaps(const Particle<D>& p, bool excludeSelf = false) const {
-      std::list<typename std::unordered_map<unsigned, Particle<D>>::iterator> o;
+  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 (!excludeSelf && !isEmpty(p.position)) {
-        o.push_back(p.position.occupant);
-      }
+    if (s.isEmpty()) {
+      return o;
+    }
+
+    if (!v.isEmpty() && !excludeSelf) {
+      o.push_back(v);
+    }
 
-      for (const HalfEdge<D>& e : p.position.adjacentEdges) {
-        if (!isEmpty(e.neighbor)) {
-          if (p.state.σ.dot(e.Δx) == 1 || e.neighbor.occupant->second.state.σ.dot(e.Δx) == -1) {
-            o.push_back(e.neighbor.occupant);
-          } 
+    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;
+    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 insert(const Particle<D>& p, Rng& r) {
-      if (overlaps(p).empty()) {
-        typename std::unordered_map<unsigned, Particle<D>>::iterator it;
-        std::tie(it, std::ignore) = particles.insert({r.uniform((unsigned)0, (unsigned)pow(2, 28)), p});
-        p.position.occupant = it;
+  bool tryRandomMove(Rng& r) {
+    Vertex<D>& v = r.pick(vertices);
+    State<D> oldState = v.state;
+    if (!tryDeletion(v)) {
+      return false;
+    }
+
+    if (0.2 > r.uniform(0.0, 1.0)) {
+      if (insert(v, State<D>(r))) {
+        return true;
+      }
+    } else {
+      if (insert(r.pick(v.adjacentEdges).neighbor, State<D>(r))) {
         return true;
-      } else {
-        return false;
       }
     }
+    v.state = oldState;
+    N++;
+    return false;
+  }
 
-    void remove(typename std::unordered_map<unsigned, Particle<D>>::iterator ip) {
-      ip->second.position.occupant = particles.end();
-      particles.erase(ip);
+  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 randomMove(Rng& r) {
-      Particle<D>& p = *(r.pick(particles));
-      Particle<D> pN = p;
+  bool tryRandomSwap(Rng& r) {
+    Vertex<D>& v1 = r.pick(vertices);
+    Vertex<D>& v2 = r.pick(vertices);
 
-      if (0.2 < r.uniform(0.0, 1.0)) {
-        pN.position = r.pick(p.position.adjacentEdges).neighbor;
+    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);
       }
-      pN.state = randomState<D>(r);
+      a %= 2 * D + 1;
 
-      if (overlaps(pN, true).empty()) {
-        remove(&p);
-        insert(pN);
-        return true;
-      } else {
-        return false;
+      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 swap(Particle<D>& p1, Particle<D>& p2) {
-      std::swap(p1.state, p2.state);
-      if (overlaps(p1, true).empty() && overlaps(p2, true).empty()) {
-        return true;
-      } else {
-        std::swap(p1.state, p2.state);
+  bool compatible() {
+    for (Vertex<D>& v : vertices) {
+      if (overlaps(v, v.state, true).size() > 0) {
         return false;
       }
     }
 
-    bool randomSwap(Rng& r) {
-      return swap(r.pick(particles).second, r.pick(particles).second);
-    }
+    return true;
+  }
 
-    void setGroundState(Rng& r) {
-      for (typename std::unordered_map<unsigned, Particle<D>>::iterator ip = particles.begin(); ip != particles.end(); ip++) {
-        remove(ip);
-      }
+  double density() const { return N / pow(L, D); }
 
-      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;
-
-        if (a > 0) {
-          if (a <= D) {
-            State<D> s(a - 1, -1);
-            insert(Particle<D>(v, s), r);
-          } else if (D < a) {
-            State<D> s(2 * D - a, 1);
-            insert(Particle<D>(v, s), r);
+  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());
 
-    bool compatible() {
-      for (Particle<D> p : particles) {
-        if (!overlaps(p, true).empty()) {
-          return false;
+        if (pDel > r.uniform(0.0, 1.0)) {
+          tryDeletion(r.pick(vertices));
         }
       }
 
-      return true;
+      tryRandomMove(r);
+      tryRandomSwap(r);
     }
+  }
 
-    double density() const {
-      return N() / pow(L, D);
+  void sweepLocal(double z, Rng& r) {
+    for (unsigned i = 0; i < iPow(L, D); i++) {
+      tryRandomMove(r);
     }
+  }
 
-    unsigned maxOccupation() const {
-//      return (2 * D * iPow(L, D)) / (2 * D + 1);
-      return iPow(L, D);
+  void sweepSwap(double z, Rng& r) {
+    for (unsigned i = 0; i < iPow(L, D); i++) {
+      tryRandomSwap(r);
     }
+  }
 
-    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 (isEmpty(v)) {
-                insert(Particle<D>(v, State<D>(r)), r);
-                break;
-              }
-            }
-          }
-        } else {
+  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)];
 
-          double pDel = N() / (z * maxOccupation());
+        v.marked = true;
+        vNew.marked = true;
 
-          if (pDel > r.uniform(0.0, 1.0)) {
-            remove(r.choose(particles));
+        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);
           }
         }
 
-//        randomMove(r);
-        randomSwap(r);
-      }
-    }
+        if (!dry) {
+          v.state = sNew;
+          vNew.state = s;
+        }
 
-    void sweepLocal(double z, Rng& r) {
-      for (unsigned i = 0; i < iPow(L, D); i++) {
-//        randomMove(r);
+        n += 1;
       }
     }
 
-    void sweepSwap(double z, Rng& r) {
-      for (unsigned i = 0; i < iPow(L, D); i++) {
-        randomSwap(r);
+    return n;
+  }
+
+  void swendsenWang(const Transformation<D>& R, Rng& r) {
+    for (Vertex<D>& v : vertices) {
+      if (!v.marked) {
+        unsigned n = flipCluster(R, v, r, 0.5 < r.uniform(0.0, 1.0));
+        std::cerr << n << " ";
       }
     }
 
-    void flipCluster(const Matrix<D>& R, Vertex<D>& v0, Rng& r, bool dry = false) {
-      std::list<std::reference_wrapper<Vertex<D>>> q = {v0};
+    std::cerr << std::endl;
 
-      while (!q.empty()) {
-        Vertex<D>& v = q.front();
-        q.pop_front();
+    for (Vertex<D>& v : vertices) {
+      v.marked = false;
+    }
+  }
+};
 
-        Vector<D> xNew = R * v.position;
-        State<D> sNew = v.state.apply(R);
+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) {
+      std::cerr << "◀";
+    } else if (v.state(0) == 0 && v.state(1) == -1) {
+      std::cerr << "▲";
+    } else if (v.state(0) == 0 && v.state(1) == 1) {
+      std::cerr << "▼";
+    } else if (v.state(0) == 0 && v.state(1) == 0) {
+      std::cerr << " ";
+    } else {
+      std::cerr << "X";
+    }
 
-        for (Vertex<D>& vn : overlaps(vertices[vectorToIndex(xNew)], sNew)) {
-          if (!vn.marked) {
-            vn.marked = true;
-            q.push_back(vn);
-          }
-        }
-      }
+    if (v.position(0) == s.L - 1) {
+      std::cerr << std::endl;
     }
-};
+  }
+}
 
-template<unsigned D> Vector<D> randomVector(unsigned L, Rng& r) {
+template <unsigned D> Vector<D> randomVector(unsigned L, Rng& r) {
   Vector<D> x;
   for (unsigned i = 0; i < D; i++) {
     x[i] = r.uniform((unsigned)0, L - 1);
@@ -309,12 +462,32 @@ int main() {
 
   Rng r;
 
+  double z = exp(1 / 0.15);
+
+  s.setGroundState();
+  std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
+
+  for (unsigned n = 0; n < 1e3; n++) {
+    s.sweepGrandCanonical(z, r);
+  }
+
+  if (s.compatible()) {
+    std::cerr << "Still compatible!" << std::endl;
+  }
+
+
+  for (unsigned n = 0; n < 10; n++) {
+    s.swendsenWang(Transformation<D>(L, ms, r), r);
+  }
+
+  if (s.compatible()) {
+    std::cerr << "Still compatible!" << std::endl;
+  }
+
   /*
+
   for (unsigned N = Nmin; N <= Nmax; N *= 10) {
     s.setGroundState();
-    if (!s.compatible()) {
-      return 1;
-    }
     for (double T = Tmin; T < Tmax; T *= 1 + δT) {
       double z = exp(1 / T);
 
@@ -337,13 +510,6 @@ int main() {
   }
   */
 
-  s.setGroundState(r);
-
-  for (unsigned n = 0; n < 10; n++) {
-    s.sweepGrandCanonical(exp(1/0.15), r);
-  }
-
   return 0;
 }
 
-
-- 
cgit v1.2.3-54-g00ecf