Rolled code into single file to trim things down.

1 files changed, 488 insertions, 7 deletions

diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
index ac68776..e5a6405 100644
--- a/biroli-mezard.cpp
+++ b/biroli-mezard.cpp
@@ -1,11 +1,492 @@
 #include <iostream>
 #include <fstream>
+#include <limits>
+#include <list>
+#include <queue>
+#include <vector>
 
-#include "glass.hpp"
+#include <eigen3/Eigen/Dense>
 
-void print(const BiroliSystem<2>& s) {
-  for (const Vertex<2, BiroliState>& v : s.vertices) {
-      std::cerr << v.state.type;
+#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>;
+
+const unsigned Empty = std::numeric_limits<unsigned>::max();
+
+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 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;
+  }
+};
+
+template <int D> class Vertex {
+public:
+  Vector<D> position;
+  std::vector<std::reference_wrapper<Vertex<D>>> neighbors;
+  unsigned occupiedNeighbors;
+  unsigned maximumNeighbors;
+  bool marked;
+
+  Vertex() {
+    occupiedNeighbors = 0;
+    maximumNeighbors = Empty;
+    marked = false;
+  }
+
+  bool empty() const { return maximumNeighbors == Empty; }
+  bool frustrated() const { return occupiedNeighbors > maximumNeighbors; }
+};
+
+template <int D> class System {
+public:
+  const unsigned L;
+  std::vector<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, unsigned T) : L(L), N(T + 1, 0), vertices(iPow(L, D)), orientation(L) {
+    N[0] = size();
+
+    for (unsigned i = 0; i < iPow(L, D); i++) {
+      vertices[i].position = indexToVector(i);
+      vertices[i].neighbors.reserve(2 * D);
+    }
+
+    for (unsigned d = 0; d < D; d++) {
+      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].neighbors.push_back(vertices[j]);
+        vertices[j].neighbors.push_back(vertices[i]);
+      }
+    }
+  }
+
+  std::list<std::reference_wrapper<Vertex<D>>> overlaps(Vertex<D>& v) {
+    std::list<std::reference_wrapper<Vertex<D>>> o;
+
+    if (v.empty()) { // an empty site cannot be frustrating anyone
+      return o;
+    }
+
+    bool selfFrustrated = v.frustrated();
+    bool anyNeighborFrustrated = false;
+
+    for (Vertex<D>& vn : v.neighbors) {
+      if (!vn.empty()) {
+        bool thisNeighborFrustrated = vn.frustrated();
+
+        if (thisNeighborFrustrated) {
+          anyNeighborFrustrated = true;
+        }
+
+        if (selfFrustrated || thisNeighborFrustrated) {
+          o.push_back(vn);
+        }
+      }
+    }
+
+    if (selfFrustrated || anyNeighborFrustrated) {
+      o.push_back(v);
+    }
+
+    return o;
+  }
+
+  bool canReplaceWith(const Vertex<D>& v, unsigned t) const {
+    if (t == Empty) {
+      return true;
+    }
+
+    if (t < v.occupiedNeighbors) {
+      return false;
+    }
+    
+    if (v.empty()) {
+      for (const Vertex<D>& vn : v.neighbors) {
+        if (vn.maximumNeighbors < vn.occupiedNeighbors + 1) {
+          return false;
+        }
+      }
+    }
+
+    return true;
+  }
+
+  bool insert(Vertex<D>& v, unsigned t, bool force = false) {
+    if (force || (v.empty() && canReplaceWith(v, t))) {
+      v.maximumNeighbors = t;
+      for (Vertex<D>& vn : v.neighbors) {
+        vn.occupiedNeighbors++;
+      }
+      N[t]++;
+      N[0]--;
+      return true;
+    } else {
+      return false;
+    }
+  }
+
+  bool remove(Vertex<D>& v) {
+    if (v.empty()) {
+      return false;
+    } else {
+      N[v.maximumNeighbors]--;
+      N[0]++;
+      v.maximumNeighbors = Empty;
+      for (Vertex<D>& vn : v.neighbors) {
+        vn.occupiedNeighbors--;
+      }
+      return true;
+    }
+  }
+
+  bool randomLocalExchange(Rng& r) {
+    Vertex<D>& v = r.pick(vertices);
+
+    return exchange(v, r.pick(v.neighbors));
+  }
+
+  void swap(Vertex<D>& v1, Vertex<D>& v2) {
+    if (v1.maximumNeighbors != v2.maximumNeighbors) {
+      if (!v1.empty() && !v2.empty()) {
+        std::swap(v1.maximumNeighbors, v2.maximumNeighbors);
+      } else if (v1.empty()) {
+        unsigned t = v2.maximumNeighbors;
+        remove(v2);
+        insert(v1, t, true);
+      } else {
+        unsigned t = v1.maximumNeighbors;
+        remove(v1);
+        insert(v2, t, true);
+      }
+    }
+  }
+
+  bool exchange(Vertex<D>& v1, Vertex<D>& v2) {
+    if (canReplaceWith(v1, v2.maximumNeighbors) && canReplaceWith(v2, v1.maximumNeighbors)) {
+      swap(v1, v2);
+      return true;
+    } else {
+      return false;
+    }
+  }
+
+  int overlap(const System<D>& s) const {
+    int o = 0;
+
+    for (unsigned i = 0; i < size(); i++) {
+      unsigned t2 = 
+          s.vertices[vectorToIndex(orientation.inverse().apply(indexToVector(i)))].maximumNeighbors;
+      if (t2 == vertices[i].maximumNeighbors) {
+        o++;
+      } else {
+        o--;
+      }
+    }
+
+    return o;
+  }
+
+  bool randomExchange(Rng& r) {
+    Vertex<D>& v1 = r.pick(vertices);
+    Vertex<D>& v2 = r.pick(vertices);
+
+    return exchange(v1, v2);
+  }
+
+  bool compatible() const {
+    for (const Vertex<D>& v : vertices) {
+      if (v.frustrated()) {
+        return false;
+      }
+    }
+
+    return true;
+  }
+
+  unsigned occupancy() const {
+    unsigned m = 0;
+
+    for (unsigned i = 1; i < N.size(); i++) {
+      m += N[i];
+    }
+
+    return m;
+  }
+
+  unsigned types() const {
+    return N.size() - 1;
+  }
+
+  double density() const { return (double)occupancy() / size(); }
+
+  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 / (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 {
+
+        double pDel = density() / z;
+
+        if (pDel > r.uniform(0.0, 1.0)) {
+          remove(r.pick(vertices));
+        }
+      }
+
+      randomLocalExchange(r);
+    }
+  }
+
+  void sweepLocal(Rng& r) {
+    for (unsigned i = 0; i < iPow(L, D); i++) {
+      randomLocalExchange(r);
+    }
+  }
+
+  void sweepSwap(Rng& r) {
+    for (unsigned i = 0; i < iPow(L, D); i++) {
+      randomExchange(r);
+    }
+  }
+
+  unsigned flipCluster(const Transformation<D>& R, Vertex<D>& v0) {
+    unsigned n = 0;
+
+    Vertex<D>& v0New = vertices[vectorToIndex(R.apply(v0.position))];
+
+    if (&v0New != &v0) {
+      std::queue<std::reference_wrapper<Vertex<D>>> q;
+
+      v0.marked = true;
+      v0New.marked = true;
+
+      swap(v0, v0New);
+
+      for (Vertex<D>& vn : overlaps(v0)) {
+        if (!vn.marked) {
+          q.push(vn);
+        }
+      }
+      for (Vertex<D>& vn : overlaps(v0New)) {
+        if (!vn.marked) {
+          q.push(vn);
+        }
+      }
+
+      while (!q.empty()) {
+        Vertex<D>& v = q.front();
+        q.pop();
+
+        if (!v.marked && !overlaps(v).empty()) {
+          v.marked = true;
+          Vertex<D>& vNew = vertices[vectorToIndex(R.apply(v.position))];
+
+          if (&vNew != &v) {
+            vNew.marked = true;
+
+            swap(v, vNew);
+
+            for (Vertex<D>& vn : overlaps(v)) {
+              if (!vn.marked) {
+                q.push(vn);
+              }
+            }
+            for (Vertex<D>& vn : overlaps(vNew)) {
+              if (!vn.marked) {
+                q.push(vn);
+              }
+            }
+
+            n += 2;
+          } else {
+            n += 1;
+          }
+        }
+        if (q.empty()) {
+          for (Vertex<D>& 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>& v : vertices) {
+      v.marked = false;
+    }
+
+    return n;
+  }
+
+};
+
+void print(const System<2>& s) {
+  for (const Vertex<2>& v : s.vertices) {
+    if (!v.empty()) {
+      std::cerr << v.maximumNeighbors;
+    } else {
+      std::cerr << " ";
+    }
 
     if (v.position(0) == s.L - 1) {
       std::cerr << std::endl;
@@ -15,14 +496,14 @@ void print(const BiroliSystem<2>& s) {
 
 int main() {
   const unsigned D = 3;
-  unsigned L = 30;
+  unsigned L = 15;
   unsigned Nmin = 2e2;
   unsigned Nmax = 2e5;
   double Tmin = 0.04;
   double Tmax = 0.2;
   double δT = 0.02;
 
-  BiroliSystem<D> s(L);
+  System<D> s(L, 3);
 
   Rng r;
 
@@ -44,7 +525,7 @@ int main() {
 
   std::cerr << "Found state with appropriate density." << std::endl;
 
-  BiroliSystem<D> s0 = s;
+  System<D> s0 = s;
 
   std::vector<Matrix<D>> ms = generateTorusMatrices<D>();