Split up functions, tired to generize in preparation for implementing Biroli-Mezard.

2 files changed, 279 insertions, 251 deletions

diff --git a/ciamarra.cpp b/ciamarra.cpp
new file mode 100644
index 0000000..8449cbb
--- /dev/null
+++ b/ciamarra.cpp
@@ -0,0 +1,79 @@
+#include <iostream>
+
+#include "glass.hpp"
+
+
+void print(const CiamarraSystem<2>& s) {
+  for (const Vertex<2, CiamarraState<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";
+    }
+
+    if (v.position(0) == s.L - 1) {
+      std::cerr << std::endl;
+    }
+  }
+}
+
+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);
+  }
+  return x;
+}
+
+int main() {
+  const unsigned D = 3;
+  unsigned L = 15;
+  unsigned Nmin = 2e2;
+  unsigned Nmax = 2e5;
+  double Tmin = 0.04;
+  double Tmax = 0.2;
+  double δT = 0.02;
+
+  CiamarraSystem<D> s(L);
+
+  Rng r;
+
+  double z = exp(1 / 0.08);
+
+  while (s.density() < 0.818) {
+    s.sweepGrandCanonical(z, r);
+  }
+
+  if (!s.compatible()) {
+    std::cerr << "Net compatible!" << std::endl;
+    return 1;
+  } 
+
+  std::cerr << "Found state with appropriate density." << std::endl;
+
+  CiamarraSystem<D> s0 = s;
+
+  std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
+
+  unsigned n = 1;
+  for (unsigned i = 0; i < 1e5; i++) {
+    if (n < 20 * log(i + 1)) {
+      n++;
+      std::cout << i << " " << s.overlap(s0) << std::endl;
+    }
+    s.sweepLocal(r);
+//    s.sweepSwap(r);
+//    s.swendsenWang(Transformation<D>(L, ms, r), r);
+  }
+
+  return 0;
+}
+
diff --git a/glass.cpp b/glass.hpp
index f2f44c1..da06542 100644
--- a/glass.cpp
+++ b/glass.hpp
@@ -1,9 +1,8 @@
-#include <eigen3/Eigen/Dense>
-#include <eigen3/Eigen/src/Core/CwiseNullaryOp.h>
-#include <iostream>
 #include <list>
-#include <vector>
 #include <queue>
+#include <vector>
+
+#include <eigen3/Eigen/Dense>
 
 #include "pcg-cpp/include/pcg_random.hpp"
 #include "randutils/randutils.hpp"
@@ -26,9 +25,7 @@ int iPow(int x, unsigned p) {
     return x * tmp * tmp;
 }
 
-unsigned mod(signed a, unsigned b) {
-  return ((a < 0) ? (a + (1 - a / (signed)b) * b) : a) % b;
-}
+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;
@@ -51,7 +48,7 @@ template <int D> void one_sequences(std::list<std::array<double, D>>& sequences,
   }
 }
 
-template <unsigned D> std::vector<Matrix<D>> generateTorusMatrices() {
+template <int D> std::vector<Matrix<D>> generateTorusMatrices() {
   std::vector<Matrix<D>> mats;
 
   std::array<double, D> ini_sequence;
@@ -105,22 +102,22 @@ template <unsigned D> std::vector<Matrix<D>> generateTorusMatrices() {
   return mats;
 }
 
-template <unsigned D> class State : public Vector<D> {
+template <int D> class CiamarraState : public Vector<D> {
 public:
-  State() : Vector<D>(Vector<D>::Zero()) {}
+  CiamarraState() : Vector<D>(Vector<D>::Zero()) {}
 
-  State(const Vector<D>& v) : Vector<D>(v) {}
+  CiamarraState(const Vector<D>& v) : Vector<D>(v) {}
 
-  State(unsigned a, signed b) : Vector<D>(Vector<D>::Zero()) { this->operator()(a) = b; }
+  CiamarraState(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})) {}
+  CiamarraState(Rng& r) : CiamarraState(r.uniform(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;
+  CiamarraState<D> flip() const {
+    CiamarraState<D> s;
     for (unsigned i = 0; i < D; i++) {
       s(i) = -this->operator()(i);
     }
@@ -128,74 +125,92 @@ public:
   }
 };
 
-template <unsigned D> class Transformation {
-  public:
-    unsigned L;
-    Matrix<D> m;
-    Vector<D> v;
+class BiroliState {
+public:
+  unsigned type;
+  unsigned occupiedNeighbors;
 
-    Transformation(unsigned L) : L(L) {
-      m.setIdentity();
-      v.setZero();
-    }
+  BiroliState() {
+    type = 0;
+    occupiedNeighbors = 0;
+  }
 
-    Transformation(unsigned L, const Matrix<D>& m, const Vector<D>& v) : L(L), m(m), v(v) {}
+  BiroliState(unsigned t, unsigned o) {
+    type = t;
+    occupiedNeighbors = o;
+  }
 
-    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);
-      }
+  bool isEmpty() const { return type == 0; }
 
-      v = v - m * v;
-    }
+  void remove() { type = 0; };
+};
 
-    Vector<D> apply(const Vector<D>& x) const {
-      return mod<D>(v + m * x, L);
-    }
+template <int D> class Transformation {
+public:
+  unsigned L;
+  Matrix<D> m;
+  Vector<D> v;
 
-    Transformation<D> apply(const Transformation<D>& t) const {
-      Transformation<D> tNew(L);
+  Transformation(unsigned L) : L(L) {
+    m.setIdentity();
+    v.setZero();
+  }
 
-      tNew.m = m * t.m;
-      tNew.v = apply(t.v);
+  Transformation(unsigned L, const Matrix<D>& m, const Vector<D>& v) : L(L), m(m), v(v) {}
 
-      return tNew;
+  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);
     }
 
-    State<D> apply(const State<D>& x) const {
-      return State<D>(m * x);
-    }
+    v = v - m * v;
+  }
 
-    Transformation<D> inverse() const {
-      return Transformation<D>(L, m.transpose(), -m.transpose() * 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 <unsigned D> class HalfEdge;
+template <int D, typename State> class HalfEdge;
 
-template <unsigned D> class Vertex {
+template <int D, typename State> class Vertex {
 public:
   Vector<D> position;
-  State<D> state;
-  std::vector<HalfEdge<D>> adjacentEdges;
+  std::vector<HalfEdge<D, State>> adjacentEdges;
+  State state;
   bool marked;
 
   bool isEmpty() const { return state.isEmpty(); }
 };
 
-template <unsigned D> class HalfEdge {
+template <int D, class State> class HalfEdge {
 public:
-  Vertex<D>& neighbor;
+  Vertex<D, State>& neighbor;
   Vector<D> Δx;
 
-  HalfEdge(Vertex<D>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
+  HalfEdge(Vertex<D, State>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
 };
 
-template <unsigned D> class System {
+template <int D, class State> class System {
 public:
   const unsigned L;
   unsigned N;
-  std::vector<Vertex<D>> vertices;
+  std::vector<Vertex<D, State>> vertices;
   Transformation<D> orientation;
 
   unsigned vectorToIndex(const Vector<D>& x) const {
@@ -226,126 +241,32 @@ public:
       Δ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));
+        vertices[i].adjacentEdges.push_back(HalfEdge<D, State>(vertices[j], Δx));
+        vertices[j].adjacentEdges.push_back(HalfEdge<D, State>(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;
+  virtual std::list<std::reference_wrapper<Vertex<D, State>>> overlaps(Vertex<D, State>&,
+                                                                       const State&, bool = false) { return {}; };
 
-    if (s.isEmpty()) {
-      return o;
-    }
+  virtual bool insert(Vertex<D, State>& v, const State& s) { return false; };
 
-    if (!v.isEmpty() && !excludeSelf) {
-      o.push_back(v);
-    }
+  virtual bool remove(Vertex<D, State>& v) { return false;};
 
-    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);
-        }
-      }
-    }
+  virtual bool tryRandomMove(Rng& r) { return false;};
 
-    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;
-    }
-  }
+  virtual bool swap(Vertex<D, State>& v1, Vertex<D, State>& v2) { 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;
+    Vertex<D, State>& v1 = r.pick(vertices);
+    Vertex<D, State>& v2 = r.pick(vertices);
 
-      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++;
-      }
-    }
+    return swap(v1, v2);
   }
 
   bool compatible() {
-    for (Vertex<D>& v : vertices) {
+    for (Vertex<D, State>& v : vertices) {
       if (overlaps(v, v.state, true).size() > 0) {
         return false;
       }
@@ -356,31 +277,28 @@ public:
 
   double density() const { return N / pow(L, D); }
 
-  unsigned maxOccupation() const {
-    //      return (2 * D * iPow(L, D)) / (2 * D + 1);
-    return iPow(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 = maxOccupation() * z / (N + 1);
+        double pIns = size() * z / (N + 1);
 
         if (pIns > r.uniform(0.0, 1.0)) {
           while (true) {
-            Vertex<D>& v = r.pick(vertices);
+            Vertex<D, State>& v = r.pick(vertices);
             if (v.isEmpty()) {
-              insert(v, State<D>(r));
+              insert(v, State(r));
               break;
             }
           }
         }
       } else {
 
-        double pDel = N / (z * maxOccupation());
+        double pDel = N / (z * size());
 
         if (pDel > r.uniform(0.0, 1.0)) {
-          tryDeletion(r.pick(vertices));
+          remove(r.pick(vertices));
         }
       }
 
@@ -400,39 +318,38 @@ public:
     }
   }
 
-  unsigned flipCluster(const Transformation<D>& R, Vertex<D>& v0, Rng& r, bool dry = false) {
-    std::queue<std::reference_wrapper<Vertex<D>>> q;
+  unsigned flipCluster(const Transformation<D>& R, Vertex<D, State>& v0, Rng& r, bool dry = false) {
+    std::queue<std::reference_wrapper<Vertex<D, State>>> q;
     q.push(v0);
 
     unsigned n = 0;
 
     while (!q.empty()) {
-      Vertex<D>& v = q.front();
+      Vertex<D, State>& v = q.front();
       q.pop();
 
       if (!v.marked) {
         Vector<D> xNew = R.apply(v.position);
-        Vertex<D>& vNew = vertices[vectorToIndex(xNew)];
+        Vertex<D, State>& vNew = vertices[vectorToIndex(xNew)];
 
         v.marked = true;
         vNew.marked = true;
 
-        State<D> s = R.apply(v.state);
-        State<D> sNew = R.apply(vNew.state);
+        CiamarraState<D> s = R.apply(v.state);
+        CiamarraState<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);
+        std::list<std::reference_wrapper<Vertex<D, State>>> overlaps1 = overlaps(vNew, s, true);
+        std::list<std::reference_wrapper<Vertex<D, State>>> overlaps2 = overlaps(v, sNew, true);
         overlaps1.splice(overlaps1.begin(), overlaps2);
 
-        for (Vertex<D>& vn : overlaps1) {
+        for (Vertex<D, State>& vn : overlaps1) {
           if (!vn.marked) {
             q.push(vn);
           }
         }
 
         if (!dry) {
-          v.state = sNew;
-          vNew.state = s;
+          swap(v, vNew);
         }
 
         n += 1;
@@ -443,7 +360,7 @@ public:
   }
 
   void swendsenWang(const Transformation<D>& R, Rng& r) {
-    for (Vertex<D>& v : vertices) {
+    for (Vertex<D, State>& v : vertices) {
       if (!v.marked) {
         bool dry = 0.5 < r.uniform(0.0, 1.0);
         unsigned n = flipCluster(R, v, r, dry);
@@ -453,101 +370,133 @@ public:
       }
     }
 
-    for (Vertex<D>& v : vertices) {
+    for (Vertex<D, State>& v : vertices) {
       v.marked = false;
     }
   }
 
-  int overlap(const System<D>& s) const {
-    int o = 0;
+  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;
+    }
 
-    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);
+    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;
   }
-};
 
-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 << " ";
+  bool insert(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s) override {
+    if (overlaps(v, s).empty()) {
+      v.state = s;
+      this->N++;
+      return true;
     } else {
-      std::cerr << "X";
-    }
-
-    if (v.position(0) == s.L - 1) {
-      std::cerr << std::endl;
+      return false;
     }
   }
-}
 
-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);
+  bool remove(Vertex<D, CiamarraState<D>>& v) override {
+    if (v.isEmpty()) {
+      return false;
+    } else {
+      v.state.remove();
+      this->N--;
+      return true;
+    }
   }
-  return x;
-}
-
-int main() {
-  const unsigned D = 3;
-  unsigned L = 15;
-  unsigned Nmin = 2e2;
-  unsigned Nmax = 2e5;
-  double Tmin = 0.04;
-  double Tmax = 0.2;
-  double δT = 0.02;
-
-  System<D> s(L);
 
-  Rng r;
+  bool tryRandomMove(Rng& r) override {
+    Vertex<D, CiamarraState<D>>& v = r.pick(this->vertices);
+    CiamarraState<D> oldState = v.state;
 
-  double z = exp(1 / 0.08);
+    if (!remove(v)) {
+      return false;
+    }
 
-  while (s.density() < 0.818) {
-    s.sweepGrandCanonical(z, r);
+    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;
   }
 
-  if (!s.compatible()) {
-    std::cerr << "Net compatible!" << std::endl;
-    return 1;
-  } 
+  bool swap(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) {
+      std::swap(v1.state, v2.state);
+      return true;
+    } else {
+      return false;
+    }
+  }
 
-  std::cerr << "Found state with appropriate density." << std::endl;
+  void setGroundState() {
+    this->N = 0;
 
-  System<D> s0 = s;
+    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;
 
-  std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
+      v.state.setZero() = Vector<D>::Zero();
 
-  unsigned n = 1;
-  for (unsigned i = 0; i < 1e5; i++) {
-    if (n < 20 * log(i + 1)) {
-      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;
+      if (0 < a && a <= D) {
+        v.state(a - 1) = -1;
+        this->N++;
+      } else if (D < a) {
+        v.state(2 * D - a) = 1;
+        this->N++;
+      }
     }
-//    s.sweepLocal(r);
-//    s.sweepSwap(r);
-//    s.swendsenWang(Transformation<D>(L, ms, r), r);
   }
 
-  return 0;
-}
+  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;
+  }
+};