Remaved self inverse scattering in favor of overlap, worked on camparisons for Ciamarra.HEADmaster

3 files changed, 69 insertions, 99 deletions

diff --git a/ciamarra.cpp b/ciamarra.cpp
index d3acc38..1d4538f 100644
--- a/ciamarra.cpp
+++ b/ciamarra.cpp
@@ -84,38 +84,6 @@ class CiamarraSystem : public HardSystem<D, CiamarraState<D>> {
       }
     }
   }
-
-  /* For the Ciamarra system, position within sites must be specially acconted
-   * for in the scattering function. We therefore expand the lattice and add
-   * the state when setting positions, so as later to be able to take particle
-   * sublattice positions into account.
-   * */
-  void setInitialPosition() override {
-    CiamarraSystem::orientation = Transformation<D>(CiamarraSystem::L);
-    for (Vertex<D, CiamarraState<D>>& v : CiamarraSystem::vertices) {
-      v.initialPosition = 4 * v.position + v.state;
-    }
-  }
-
-  double selfIntermediateScattering(const std::vector<Matrix<D>>& ms) const override {
-    double F = 0;
-
-    std::array<Vector<D>, D> ks;
-    for (unsigned i = 0; i < D; i++) {
-      ks[i].setZero();
-      ks[i](i) = 1;
-    }
-    for (const Vertex<D, CiamarraState<D>>& v : CiamarraSystem::vertices) {
-      if (!v.empty()) {
-        Vector<D> Δx = ((4 * CiamarraSystem::orientation.inverse().apply(v.position)) + CiamarraSystem::orientation.inverse().apply(v.state)) - v.initialPosition;
-        for (const Vector<D>& k : ks) {
-          F += cos((M_PI / 4) * k.dot(Δx));
-        }
-      }
-    }
-
-    return F / (D * CiamarraSystem::N);
-  }
 };
 
 void print(const CiamarraSystem<2>& s) {
@@ -161,27 +129,38 @@ int main() {
     while (s.density() < ρ) {
       s.sweepGrandCanonical(z, r);
     }
-    std::cout << s.density() << " ";
-    s.setInitialPosition();
+
+    std::cout << s.N << std::endl;
+
+    CiamarraSystem<D> s0 = s;
+    double last_overlap = 1;
     auto start = std::chrono::high_resolution_clock::now();
-    for (unsigned i = 0; i < 1e4; i++) {
-      /*
+
+    while (last_overlap > 0.1) {
       for (unsigned j = 0; j < s.vertices.size(); j++) {
-        s.tryRandomNonlocalMove(r);
+        s.tryRandomMove(r);
       }
-      */
+
+      auto stop = std::chrono::high_resolution_clock::now();
+      auto duration = duration_cast<std::chrono::microseconds>(stop - start);
+      last_overlap = ((double)s.overlap(s0) / s.N - s.density() / 6) / (1 - s.density() / 6);
+      std::cout << duration.count() << " " << last_overlap << " ";
+    }
+    std::cout << std::endl;
+
+    CiamarraSystem<D> s1 = s;
+    last_overlap = 1;
+    auto start1 = std::chrono::high_resolution_clock::now();
+
+    while (last_overlap > 0.1) {
       Matrix<D> m = r.pick(ms);
       Vertex<D, CiamarraState<D>>& v = r.pick(s.vertices);
-      unsigned n = s.wolff(Transformation<D>(L, m, v.position - m * v.position), r);
-      if (i % 1 == 0) {
-        auto stop = std::chrono::high_resolution_clock::now();
-        auto duration = duration_cast<std::chrono::microseconds>(stop - start);
-        std::cout << duration.count() << " " << s.selfIntermediateScattering(ms) << " " << n << " ";
-      }
-    }
-    if (!s.compatible()) {
-      std::cerr << "Bad configuration" << std::endl;
-      return 1;
+      s.wolff(Transformation<D>(L, m, v.position - m * v.position), r.pick(s.vertices));
+
+      auto stop = std::chrono::high_resolution_clock::now();
+      auto duration = duration_cast<std::chrono::microseconds>(stop - start1);
+      last_overlap = ((double)s.overlap(s1) / s.N - s.density() / 6) / (1 - s.density() / 6);
+      std::cout << duration.count() << " " << last_overlap << " ";
     }
     std::cout << std::endl;
   }
diff --git a/distinguishable.cpp b/distinguishable.cpp
index 0aae0c1..de71c47 100644
--- a/distinguishable.cpp
+++ b/distinguishable.cpp
@@ -53,7 +53,7 @@ int main() {
   unsigned N = 1560;
   double Tmin = 0;
   double Tmax = 0.4;
-  double δT = 0.02;
+  double δT = 0.05;
 
   Rng r;
 
@@ -86,23 +86,29 @@ int main() {
     /*
      */
     std::cout << T << " ";
+    for (unsigned i = 0; i < 1e3; i++) {
+      for (unsigned j = 0; j < s.vertices.size(); j++) {
+        s.tryRandomSwap(T, r);
+      }
+    }
     s.setInitialPosition();
+    DistinguishableSystem<D> s0 = s;
     auto start = std::chrono::high_resolution_clock::now();
-    Quantity<double> energy(1e3);
+//    Quantity<double> energy(1e3);
     for (unsigned i = 0; i < 1e5; i++) {
+      /*
       for (unsigned j = 0; j < s.vertices.size(); j++) {
         s.tryRandomSwap(T, r);
       }
-      /*
+      */
       Matrix<D> m = r.pick(ms);
       Vertex<D, DistinguishableState>& v = r.pick(s.vertices);
       s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r);
-      */
       if (i % 10 == 0) {
-        energy.add(s.energy());
+//        energy.add(s.energy());
         auto stop = std::chrono::high_resolution_clock::now();
         auto duration = duration_cast<std::chrono::microseconds>(stop - start);
-        std::cout /*<< duration.count() << " "*/ << s.selfIntermediateScattering(ms) << " ";
+        std::cout << duration.count() << " " << s.overlap(s0) << " ";
       }
     }
     /*
@@ -114,12 +120,14 @@ int main() {
     }
     */
     std::cout << std::endl;
+    /*
     std::vector<double> rho = energy.ρ();
     for (const double& x : rho) {
       std::cout << x << " ";
     }
     std::cout << std::endl;
     std::cerr << T << " " << s.energy() / N << std::endl;
+    */
     //    s.sweepLocal(r);
     //    s.sweepSwap(r);
     //    s.swendsenWang(Transformation<D>(L, ms, r), r);
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;
-  }
 };