1 files changed, 28 insertions, 45 deletions

diff --git a/glass.hpp b/glass.hpp
index 8ad0b52..1b1beb7 100644
--- a/glass.hpp
+++ b/glass.hpp
@@ -153,6 +153,10 @@ public:
     CiamarraState<D> s2(m * s);
     return s2;
   }
+
+  template <class A> A apply(const A& s) const {
+    return s;
+  };
 };
 
 template <typename T> concept State = requires(T v, const T& v2) {
@@ -166,7 +170,6 @@ template <unsigned D, State S> class HalfEdge;
 template <unsigned D, State S> class Vertex {
 public:
   Vector<D> position;
-  Vector<D> initialPosition;
   S state;
   std::vector<HalfEdge<D, S>> adjacentEdges;
   bool marked;
@@ -208,7 +211,6 @@ public:
   System(unsigned L) : L(L), N(0), vertices(iPow(L, D)), orientation(L) {
     for (unsigned i = 0; i < iPow(L, D); i++) {
       vertices[i].position = indexToVector(i);
-      vertices[i].initialPosition = vertices[i].position;
       vertices[i].adjacentEdges.reserve(2 * D);
       vertices[i].marked = false;
     }
@@ -232,31 +234,22 @@ public:
     return iPow(L, D);
   }
 
-  virtual void setInitialPosition() {
-    orientation = Transformation<D>(L);
-    for (Vertex<D, S>& v : vertices) {
-      v.initialPosition = v.position;
-    }
-  }
+  int overlap(const System<D, S>& s) const {
+    unsigned o = 0;
+    Transformation<D> oRel = orientation.apply(s.orientation.inverse());
 
-  virtual double selfIntermediateScattering(const std::vector<Matrix<D>>& ms) const {
-    double F = 0;
+    for (unsigned i = 0; i < vertices.size(); i++) {
+      if (!vertices[i].empty()) {
+        S s1 = vertices[i].state; 
+        S s2 = oRel.apply(s.vertices[vectorToIndex(oRel.inverse().apply(indexToVector(i)))].state);
 
-    std::array<Vector<D>, D> ks;
-    for (unsigned i = 0; i < D; i++) {
-      ks[i].setZero();
-      ks[i](i) = 1;
-    }
-    for (const Vertex<D, S>& v : vertices) {
-      if (!v.empty()) {
-        Vector<D> Δx = v.position - orientation.apply(v.initialPosition);
-        for (const Vector<D>& k : ks) {
-          F += cos(M_PI * k.dot(Δx));
+        if (s1 == s2) {
+          o += 1;
         }
       }
     }
 
-    return F / (D * System::N);
+    return o;
   }
 };
 
@@ -294,7 +287,6 @@ public:
 
     if (exp(-ΔE / T) > r.uniform(0.0, 1.0)) {
       /* Accept the swap. */
-      std::swap(v1.initialPosition, v2.initialPosition);
       return true;
     } else {
       /* Revert the swap. */
@@ -381,7 +373,6 @@ public:
 
         if (!dry) {
           std::swap(v.state, vNew.state);
-          std::swap(v.initialPosition, vNew.initialPosition);
         }
 
         if (!v.empty()) {
@@ -453,7 +444,6 @@ public:
   bool tryRandomMove(Rng& r) {
     Vertex<D, S>& v = r.pick(HardSystem::vertices);
     S oldState = v.state;
-    Vector<D> oldPos = v.initialPosition;
 
     if (!tryDeletion(v)) {
       return false;
@@ -463,7 +453,6 @@ public:
       for (HalfEdge<D, S>& e : v.adjacentEdges) {
         if (1 == e.Δx.dot(oldState)) {
           if (insert(e.neighbor, oldState.flip())) {
-            e.neighbor.initialPosition = oldPos;
             return true;
           }
           break;
@@ -475,11 +464,9 @@ public:
         newState = S(r);
       }
       if (insert(v, newState)) {
-        v.initialPosition = oldPos;
         return true;
       }
     }
-    v.initialPosition = oldPos;
     v.state = oldState;
     HardSystem::N++;
     return false;
@@ -488,7 +475,6 @@ public:
   bool tryRandomNonlocalMove(Rng& r) {
     Vertex<D, S>& v = r.pick(HardSystem::vertices);
     S oldState = v.state;
-    Vector<D> oldPos = v.initialPosition;
 
     if (!tryDeletion(v)) {
       return false;
@@ -498,11 +484,9 @@ public:
     S newState(r);
 
     if (insert(vNew, newState)) {
-      vNew.initialPosition = oldPos;
       return true;
     }
 
-    v.initialPosition = oldPos;
     v.state = oldState;
     HardSystem::N++;
     return false;
@@ -514,7 +498,6 @@ public:
     }
     if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) {
       std::swap(v1.state, v2.state);
-      std::swap(v1.initialPosition, v2.initialPosition);
       return true;
     } else {
       return false;
@@ -578,7 +561,7 @@ public:
     }
   }
 
-  unsigned flipCluster(const Transformation<D>& R, Vertex<D, S>& v0, Rng& r, bool dry = false) {
+  unsigned flipCluster(const Transformation<D>& R, Vertex<D, S>& v0, bool dry = false) {
     std::queue<std::reference_wrapper<Vertex<D, S>>> q;
     q.push(v0);
 
@@ -612,7 +595,6 @@ public:
           if (!dry) {
             v.state = sNew;
             vNew.state = s;
-            std::swap(v.initialPosition, vNew.initialPosition);
           }
 
           if (!v.empty()) {
@@ -629,8 +611,19 @@ public:
     return n;
   }
 
+  unsigned wolff(const Transformation<D>& R, Vertex<D, S>& v0) {
+    unsigned n = flipCluster(R, v0);
+    if (n > HardSystem::N / 2) {
+      HardSystem::orientation = R.apply(HardSystem::orientation);
+    }
+    for (Vertex<D, S>& v : this->vertices) {
+      v.marked = false;
+    }
+    return n;
+  }
+
   unsigned wolff(const Transformation<D>& R, Rng& r) {
-    unsigned n = flipCluster(R, r.pick(this->vertices), r);
+    unsigned n = flipCluster(R, r.pick(this->vertices));
     if (n > HardSystem::N / 2) {
       HardSystem::orientation = R.apply(HardSystem::orientation);
     }
@@ -644,7 +637,7 @@ public:
     for (Vertex<D, S>& v : HardSystem::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 > pow(HardSystem::L, D) / 4 && !dry) {
           HardSystem::orientation = R.apply(HardSystem::orientation);
         }
@@ -656,14 +649,4 @@ public:
     }
   }
 
-  int overlap(const System<D, S>& s) const {
-    int o = 0;
-
-    for (unsigned i = 0; i < HardSystem::vertices.size(); i++) {
-      S s2 = HardSystem::orientation.apply(s.vertices[HardSystem::vectorToIndex(HardSystem::orientation.inverse().apply(HardSystem::indexToVector(i)))].state);
-      o += HardSystem::vertices[i].state.dot(s2);
-    }
-
-    return o;
-  }
 };