Lots of work.

1 files changed, 121 insertions, 42 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;
 }