1 files changed, 59 insertions, 25 deletions

diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
index c1ed606..6f390e3 100644
--- a/biroli-mezard.cpp
+++ b/biroli-mezard.cpp
@@ -145,6 +145,7 @@ public:
 
 template <int D> class Vertex {
 public:
+  Vector<D> initialPosition;
   Vector<D> position;
   std::vector<std::reference_wrapper<Vertex<D>>> neighbors;
   unsigned occupiedNeighbors;
@@ -170,6 +171,9 @@ public:
 
     return o == occupiedNeighbors;
   }
+  void setInitial() {
+    initialPosition = position;
+  }
 };
 
 template <int D> class System {
@@ -293,6 +297,12 @@ public:
     return true;
   }
 
+  void setInitial() {
+    for (Vertex<D>& v : vertices) {
+      v.setInitial();
+    }
+  }
+
   bool insert(Vertex<D>& v, unsigned t, bool force = false) {
     if (force || (v.empty() && canReplaceWith(v, t))) {
       v.maximumNeighbors = t;
@@ -340,6 +350,7 @@ public:
         remove(v1);
         insert(v2, t, true);
       }
+      std::swap(v1.initialPosition, v2.initialPosition);
     }
   }
 
@@ -352,20 +363,41 @@ public:
     }
   }
 
-  int overlap(const System<D>& s) const {
+  std::complex<double> selfIntermediateScattering(unsigned k, unsigned t = 2) const {
+    std::complex<double> F = 0;
+
+    using namespace std::complex_literals;
+
+    for (const Vertex<D>& v : vertices) {
+      if (v.maximumNeighbors == t) {
+        for (unsigned d = 0; d < D; d++) {
+          F += exp(1i * 2.0 * M_PI * (double)k * (double)(v.position(d) - v.initialPosition(d)) / (double)L);
+        }
+      }
+    }
+
+    F /= D * N[t];
+
+    return F;
+  }
+
+  double overlap(const System<D>& s, unsigned t = 2) const {
     int o = 0;
 
-    for (unsigned i = 0; i < size(); i++) {
-      unsigned t2 =
-          s.vertices[vectorToIndex(orientation.inverse().apply(indexToVector(i)))].maximumNeighbors;
-      if (t2 == vertices[i].maximumNeighbors) {
+    for (const Vertex<D>& v1 : vertices) {
+      const Vertex<D>& v2 = s.vertices[vectorToIndex(orientation.inverse().apply(v1.position))];
+      if (v1.maximumNeighbors == v2.maximumNeighbors) {
         o++;
-      } else {
-        o--;
       }
     }
 
-    return o;
+    unsigned nn = 0;
+
+    for (unsigned n : N) {
+      nn += iPow(n, 2);
+    }
+
+    return ((double)o / (double)size() - (double)nn / pow((double)size(), 2)) / (1 - (double)nn / pow((double)size(), 2));
   }
 
   bool randomExchange(Rng& r) {
@@ -403,13 +435,13 @@ public:
   }
 
   void sweepLocal(Rng& r) {
-    for (unsigned i = 0; i < size(); i++) {
+    for (unsigned i = 0; i < size() / 2; i++) {
       randomLocalExchange(r);
     }
   }
 
   void sweepSwap(Rng& r) {
-    for (unsigned i = 0; i < size(); i++) {
+    for (unsigned i = 0; i < size() / 2; i++) {
       randomExchange(r);
     }
   }
@@ -420,7 +452,7 @@ public:
     Vertex<D>& v0New = vertices[vectorToIndex(R.apply(v0.position))];
 
     if (&v0New != &v0) {
-      std::queue<std::reference_wrapper<Vertex<D>>> q;
+      std::list<std::reference_wrapper<Vertex<D>>> q;
 
       v0.marked = true;
       v0New.marked = true;
@@ -429,18 +461,18 @@ public:
 
       for (Vertex<D>& vn : overlaps(v0)) {
         if (!vn.marked) {
-          q.push(vn);
+          q.push_front(vn);
         }
       }
       for (Vertex<D>& vn : overlaps(v0New)) {
         if (!vn.marked) {
-          q.push(vn);
+          q.push_front(vn);
         }
       }
 
       while (!q.empty()) {
         Vertex<D>& v = q.front();
-        q.pop();
+        q.pop_front();
 
         if (!v.marked && !overlaps(v).empty()) {
           v.marked = true;
@@ -453,13 +485,13 @@ public:
 
             for (Vertex<D>& vn : overlaps(v)) {
               if (!vn.marked) {
-                q.push(vn);
+                q.push_front(vn);
               }
             }
 
             for (Vertex<D>& vn : overlaps(vNew)) {
               if (!vn.marked) {
-                q.push(vn);
+                q.push_front(vn);
               }
             }
 
@@ -469,11 +501,13 @@ public:
           }
         }
 
-        // For some reason, the process above misses frustrated states. This is an inelegant stopgap.
+        // 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()) {
-              q.push(vv);
+              std::cerr << "Found bad vertex." << std::endl;
+              q.push_front(vv);
             }
           }
         }
@@ -528,7 +562,7 @@ int main() {
     return 1;
   }
 
-  while (s.density() < 0.58) {
+  while (s.density() < 0.56) {
     s.sweepGrandCanonical(z, r);
   }
 
@@ -545,22 +579,22 @@ int main() {
 
   std::vector<unsigned> clusterDist(s.size() + 1);
 
+  s.setInitial();
+
   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;
+      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;
     clusterDist[nn]++;
-    /*
-    s.sweepLocal(r);
-    nC += s.size();
-    s.sweepSwap(r);
-    */
+    //s.sweepLocal(r);
+    //nC += s.size();
+    //s.sweepSwap(r);
     i++;
   }