Working Biroli-Mezard clusters.

2 files changed, 289 insertions, 36 deletions

diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
new file mode 100644
index 0000000..ac68776
--- /dev/null
+++ b/biroli-mezard.cpp
@@ -0,0 +1,84 @@
+#include <iostream>
+#include <fstream>
+
+#include "glass.hpp"
+
+void print(const BiroliSystem<2>& s) {
+  for (const Vertex<2, BiroliState>& v : s.vertices) {
+      std::cerr << v.state.type;
+
+    if (v.position(0) == s.L - 1) {
+      std::cerr << std::endl;
+    }
+  }
+}
+
+int main() {
+  const unsigned D = 3;
+  unsigned L = 30;
+  unsigned Nmin = 2e2;
+  unsigned Nmax = 2e5;
+  double Tmin = 0.04;
+  double Tmax = 0.2;
+  double δT = 0.02;
+
+  BiroliSystem<D> s(L);
+
+  Rng r;
+
+  double z = exp(1 / 0.2);
+
+  if (!s.compatible()) {
+    std::cerr << "Storted incompatible!" << std::endl;
+    return 1;
+  } 
+
+  while (s.density() < 0.57) {
+    s.sweepGrandCanonical(z, r);
+  }
+
+  if (!s.compatible()) {
+    std::cerr << "Not compatible!" << std::endl;
+    return 1;
+  } 
+
+  std::cerr << "Found state with appropriate density." << std::endl;
+
+  BiroliSystem<D> s0 = s;
+
+  std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
+
+  std::vector<unsigned> clusterDist(s.size() + 1);
+
+  unsigned n = 1;
+  unsigned i = 0;
+  double nC = 0;
+  while (nC / s.size() < 1e5) {
+    if (n < 20 * log(i + 1)) {
+      n++;
+      std::cout << nC / s.size() << " " << (double)s.overlap(s0) / s.size() << std::endl;
+    }
+    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);
+//    s.swendsenWang(Transformation<D>(L, ms, r), r);
+    i++;
+  }
+
+  if (!s.compatible()) {
+    std::cerr << "Not compatible!" << std::endl;
+    return 1;
+  } 
+
+  std::ofstream file("dist.dat");
+  for (unsigned i : clusterDist) {
+    file << i << " ";
+  }
+  file.close();
+
+  return 0;
+}
+
diff --git a/glass.hpp b/glass.hpp
index e67bb51..b35964b 100644
--- a/glass.hpp
+++ b/glass.hpp
@@ -1,3 +1,4 @@
+#include <limits>
 #include <list>
 #include <queue>
 #include <vector>
@@ -130,7 +131,7 @@ public:
   unsigned occupiedNeighbors;
 
   BiroliState() {
-    type = 0;
+    type = std::numeric_limits<unsigned>::max();
     occupiedNeighbors = 0;
   }
 
@@ -139,9 +140,13 @@ public:
     occupiedNeighbors = o;
   }
 
-  bool isEmpty() const { return type == 0; }
+  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 = 0; };
+  void remove() { type = std::numeric_limits<unsigned>::max(); };
 };
 
 template <int D> class Transformation {
@@ -246,10 +251,9 @@ public:
     }
   }
 
-  virtual std::list<std::reference_wrapper<Vertex<D, State>>> overlaps(Vertex<D, State>&,
-                                                                       const State&, bool = false) { return {}; };
+  virtual std::list<std::reference_wrapper<Vertex<D, State>>> overlaps(Vertex<D, State>&) { return {}; };
 
-  virtual bool insert(Vertex<D, State>& v, const State& s) { return false; };
+  virtual bool insert(Vertex<D, State>& v, const State& s, bool force = false) { return false; };
 
   virtual bool remove(Vertex<D, State>& v) { return false;};
 
@@ -268,7 +272,7 @@ public:
 
   bool compatible() {
     for (Vertex<D, State>& v : vertices) {
-      if (overlaps(v, v.state, true).size() > 0) {
+      if (overlaps(v).size() > 0) {
         return false;
       }
     }
@@ -319,43 +323,78 @@ public:
     }
   }
 
-  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 flipCluster(const Transformation<D>& R, Vertex<D, State>& v0, bool dry = false) {
     unsigned n = 0;
 
-    while (!q.empty()) {
-      Vertex<D, State>& v = q.front();
-      q.pop();
+    Vertex<D, State>& v0New = vertices[vectorToIndex(R.apply(v0.position))];
 
-      if (!v.marked) {
-        Vertex<D, State>& vNew = vertices[vectorToIndex(R.apply(v.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);
+        }
+      }
 
-        v.marked = true;
-        vNew.marked = true;
+      while (!q.empty()) {
+        Vertex<D, State>& v = q.front();
+        q.pop();
 
-        State s = R.apply(v.state);
-        State sNew = R.apply(vNew.state);
+        if (!v.marked && !overlaps(v).empty()) {
+          v.marked = true;
+          Vertex<D, State>& vNew = vertices[vectorToIndex(R.apply(v.position))];
 
-        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);
+          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);
+              }
+            }
 
-        for (Vertex<D, State>& vn : overlaps1) {
-          if (!vn.marked) {
-            q.push(vn);
+            n += 2;
+          } else {
+            n += 1;
           }
         }
-
-        if (!dry) {
-          swap(v, vNew);
-          v.state = sNew;
-          vNew.state = s;
+        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);
+            }
+          }
         }
-
-        n += 1;
       }
+
+
+    }
+
+    if (n > size() / 4) {
+      orientation = R.apply(orientation);
+    }
+
+    for (Vertex<D, State>& v : vertices) {
+      v.marked = false;
     }
 
     return n;
@@ -365,7 +404,7 @@ public:
     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);
+        unsigned n = flipCluster(R, v, dry);
         if (n > size() / 4 && !dry) {
           orientation = R.apply(orientation);
         }
@@ -408,8 +447,8 @@ template <int D> class CiamarraSystem : public System<D, CiamarraState<D>> {
     return o;
   }
 
-  bool insert(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s) override {
-    if (overlaps(v, s).empty()) {
+  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;
@@ -506,3 +545,133 @@ template <int D> class CiamarraSystem : public System<D, CiamarraState<D>> {
     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;
+  }
+};