New way of doing clusters without missing sites: check next nearest neighbors...

1 files changed, 35 insertions, 42 deletions

diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
index 13f3b3f..7e200b9 100644
--- a/biroli-mezard.cpp
+++ b/biroli-mezard.cpp
@@ -4,6 +4,7 @@
 #include <list>
 #include <queue>
 #include <vector>
+#include <chrono>
 
 #include <eigen3/Eigen/Dense>
 
@@ -151,7 +152,6 @@ public:
   unsigned occupiedNeighbors;
   unsigned maximumNeighbors;
   bool marked;
-  bool visited;
 
   Vertex() {
     occupiedNeighbors = 0;
@@ -470,57 +470,45 @@ public:
 
       swap(v0, v0New);
 
-      for (Vertex<D>& vn : v0.neighbors) {
-        if ((!vn.marked && !vn.visited) && !vn.empty()) {
+      for (Vertex<D>& vn : overlaps({v0, v0New})) {
+        if (!vn.marked) {
           q.push(vn);
-          vn.visited = true;
-        }
-      }
-      for (Vertex<D>& vn : v0New.neighbors) {
-        if ((!vn.marked && !vn.visited) && !vn.empty()) {
-          q.push(vn);
-          vn.visited = true;
         }
       }
 
-      unsigned last_unsatisfied = 0;
-
-      while (!q.empty() && last_unsatisfied < q.size()) {
+      while (!q.empty()) {
         Vertex<D>& v = q.front();
         q.pop();
 
-        if (!v.marked) {
-          if (!overlaps(v).empty()) {
-            last_unsatisfied = 0;
-            v.marked = true;
-            Vertex<D>& vNew = vertices[vectorToIndex(R.apply(v.position))];
+        if (!v.marked && !overlaps(v).empty()) {
+          v.marked = true;
+          Vertex<D>& vNew = vertices[vectorToIndex(R.apply(v.position))];
 
-            if (&vNew != &v) {
-              vNew.marked = true;
+          if (&vNew != &v) {
+            vNew.marked = true;
 
-              swap(v, vNew);
+            swap(v, vNew);
 
-              for (Vertex<D>& vn : v.neighbors) {
-                if ((!vn.marked && !vn.visited) && !vn.empty()) {
-                  q.push(vn);
-                  vn.visited = true;
-                }
-              }
-              for (Vertex<D>& vn : vNew.neighbors) {
-                if ((!vn.marked && !vn.visited) && !vn.empty()) {
-                  q.push(vn);
-                  vn.visited = true;
+            for (Vertex<D>& vn : overlaps({v, vNew})) {
+              if (!vn.marked) {
+                q.push(vn);
+              } else {
+                /* If a neighbor has already been moved but is still
+                 * frustrated, it may be due to one of its neighbors!
+                 */
+                if (&vn != &v && &vn != &vNew) {
+                  for (Vertex<D>& vnn : overlaps(vn)) {
+                    if (!vnn.marked) {
+                      q.push(vnn);
+                    }
+                  }
                 }
               }
-
-              n += 2;
-            } else {
             }
+            n += 2;
           } else {
-            q.push(v);
-            last_unsatisfied++;
+            n += 1;
           }
-          n += 1;
         }
       }
     } else {
@@ -533,7 +521,6 @@ public:
 
     for (Vertex<D>& v : vertices) {
       v.marked = false;
-      v.visited = false;
     }
 
     return n;
@@ -596,19 +583,25 @@ int main() {
   unsigned n = 1;
   unsigned i = 0;
   double nC = 0;
+  auto start = std::chrono::high_resolution_clock::now();
   while (nC / s.size() < 1e5) {
     if (n < 20 * log(i + 1)) {
       n++;
       std::cout << nC / s.size() << " " << 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);
+    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;