1 files changed, 64 insertions, 39 deletions

diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
index 35360d0..c1ed606 100644
--- a/biroli-mezard.cpp
+++ b/biroli-mezard.cpp
@@ -159,6 +159,17 @@ public:
 
   bool empty() const { return maximumNeighbors == Empty; }
   bool frustrated() const { return occupiedNeighbors > maximumNeighbors; }
+  bool valid() const {
+    unsigned o = 0;
+
+    for (const Vertex<D>& v : neighbors) {
+      if (!v.empty()) {
+        o++;
+      }
+    }
+
+    return o == occupiedNeighbors;
+  }
 };
 
 template <int D> class System {
@@ -168,6 +179,36 @@ public:
   std::vector<Vertex<D>> vertices;
   Transformation<D> orientation;
 
+  unsigned size() const { return vertices.size(); }
+
+  unsigned types() const { return N.size() - 1; }
+
+  unsigned occupancy() const {
+    return size() - N[0];
+  }
+
+  double density() const { return (double)occupancy() / size(); }
+
+  bool compatible() const {
+    for (const Vertex<D>& v : vertices) {
+      if (v.frustrated()) {
+        return false;
+      }
+    }
+
+    return true;
+  }
+
+  bool valid() const {
+    for (const Vertex<D>& v : vertices) {
+      if (!v.valid()) {
+        return false;
+      }
+    }
+
+    return true;
+  }
+
   unsigned vectorToIndex(const Vector<D>& x) const {
     unsigned i = 0;
     for (unsigned d = 0; d < D; d++) {
@@ -187,13 +228,13 @@ public:
   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++) {
+    for (unsigned i = 0; i < size(); 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++) {
+      for (signed i = 0; i < size(); 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]);
@@ -334,34 +375,8 @@ public:
     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++) {
+    for (unsigned i = 0; i < size(); i++) {
       if (0.5 < r.uniform(0.0, 1.0)) {
         double pIns = size() * z / (occupancy() + 1);
 
@@ -388,13 +403,13 @@ public:
   }
 
   void sweepLocal(Rng& r) {
-    for (unsigned i = 0; i < iPow(L, D); i++) {
+    for (unsigned i = 0; i < size(); i++) {
       randomLocalExchange(r);
     }
   }
 
   void sweepSwap(Rng& r) {
-    for (unsigned i = 0; i < iPow(L, D); i++) {
+    for (unsigned i = 0; i < size(); i++) {
       randomExchange(r);
     }
   }
@@ -441,6 +456,7 @@ public:
                 q.push(vn);
               }
             }
+
             for (Vertex<D>& vn : overlaps(vNew)) {
               if (!vn.marked) {
                 q.push(vn);
@@ -452,18 +468,21 @@ public:
             n += 1;
           }
         }
+
+        // For some reason, the process above misses frustrated states. This is an inelegant stopgap.
         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);
             }
           }
         }
       }
+    } else {
+      n = 1;
     }
 
-    if (n > size() / 4) {
+    if (n > occupancy() / 2) {
       orientation = R.apply(orientation);
     }
 
@@ -502,14 +521,14 @@ int main() {
 
   Rng r;
 
-  double z = exp(1 / 0.2);
+  double z = exp(1 / 0.15);
 
   if (!s.compatible()) {
     std::cerr << "Storted incompatible!" << std::endl;
     return 1;
   }
 
-  while (s.density() < 0.57) {
+  while (s.density() < 0.58) {
     s.sweepGrandCanonical(z, r);
   }
 
@@ -537,10 +556,11 @@ int main() {
     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);
+    /*
+    s.sweepLocal(r);
+    nC += s.size();
+    s.sweepSwap(r);
+    */
     i++;
   }
 
@@ -549,6 +569,11 @@ int main() {
     return 1;
   }
 
+  if (!s.valid()) {
+    std::cerr << "Not valid!" << std::endl;
+    return 1;
+  }
+
   std::ofstream file("dist.dat");
   for (unsigned i : clusterDist) {
     file << i << " ";