Fixed Ciamarra model.

2 files changed, 27 insertions, 18 deletions

diff --git a/ciamarra.cpp b/ciamarra.cpp
index 8449cbb..6b932bf 100644
--- a/ciamarra.cpp
+++ b/ciamarra.cpp
@@ -58,6 +58,9 @@ int main() {
   } 
 
   std::cerr << "Found state with appropriate density." << std::endl;
+    if (!s.compatible()) {
+      std::cerr <<"whoops!" << std::endl;
+    }
 
   CiamarraSystem<D> s0 = s;
 
@@ -69,9 +72,9 @@ int main() {
       n++;
       std::cout << i << " " << s.overlap(s0) << std::endl;
     }
-    s.sweepLocal(r);
+//    s.sweepLocal(r);
 //    s.sweepSwap(r);
-//    s.swendsenWang(Transformation<D>(L, ms, r), r);
+    s.swendsenWang(Transformation<D>(L, ms, r), r);
   }
 
   return 0;
diff --git a/glass.hpp b/glass.hpp
index da06542..e67bb51 100644
--- a/glass.hpp
+++ b/glass.hpp
@@ -38,7 +38,6 @@ template <int D, typename Derived> Vector<D> mod(const Eigen::MatrixBase<Derived
 template <int D> void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
   if (level > 0) {
     unsigned new_level = level - 1;
-    unsigned old_length = sequences.size();
     for (std::array<double, D>& sequence : sequences) {
       std::array<double, D> new_sequence = sequence;
       new_sequence[new_level] = -1;
@@ -254,15 +253,17 @@ public:
 
   virtual bool remove(Vertex<D, State>& v) { return false;};
 
-  virtual bool tryRandomMove(Rng& r) { return false;};
+  virtual bool randomMove(Rng& r) { return false;};
 
-  virtual bool swap(Vertex<D, State>& v1, Vertex<D, State>& v2) { return false;};
+  virtual void swap(Vertex<D, State>& v1, Vertex<D, State>& v2) {};
 
-  bool tryRandomSwap(Rng& r) {
+  virtual bool exchange(Vertex<D, State>& v1, Vertex<D, State>& v2) { return false;};
+
+  bool randomExchange(Rng& r) {
     Vertex<D, State>& v1 = r.pick(vertices);
     Vertex<D, State>& v2 = r.pick(vertices);
 
-    return swap(v1, v2);
+    return exchange(v1, v2);
   }
 
   bool compatible() {
@@ -302,19 +303,19 @@ public:
         }
       }
 
-      tryRandomMove(r);
+      randomMove(r);
     }
   }
 
   void sweepLocal(Rng& r) {
     for (unsigned i = 0; i < iPow(L, D); i++) {
-      tryRandomMove(r);
+      randomMove(r);
     }
   }
 
   void sweepSwap(Rng& r) {
     for (unsigned i = 0; i < iPow(L, D); i++) {
-      tryRandomSwap(r);
+      randomExchange(r);
     }
   }
 
@@ -329,14 +330,13 @@ public:
       q.pop();
 
       if (!v.marked) {
-        Vector<D> xNew = R.apply(v.position);
-        Vertex<D, State>& vNew = vertices[vectorToIndex(xNew)];
+        Vertex<D, State>& vNew = vertices[vectorToIndex(R.apply(v.position))];
 
         v.marked = true;
         vNew.marked = true;
 
-        CiamarraState<D> s = R.apply(v.state);
-        CiamarraState<D> sNew = R.apply(vNew.state);
+        State s = R.apply(v.state);
+        State sNew = R.apply(vNew.state);
 
         std::list<std::reference_wrapper<Vertex<D, State>>> overlaps1 = overlaps(vNew, s, true);
         std::list<std::reference_wrapper<Vertex<D, State>>> overlaps2 = overlaps(v, sNew, true);
@@ -350,6 +350,8 @@ public:
 
         if (!dry) {
           swap(v, vNew);
+          v.state = sNew;
+          vNew.state = s;
         }
 
         n += 1;
@@ -364,7 +366,7 @@ public:
       if (!v.marked) {
         bool dry = 0.5 < r.uniform(0.0, 1.0);
         unsigned n = flipCluster(R, v, r, dry);
-        if (n > pow(L, D) / 4 && !dry) {
+        if (n > size() / 4 && !dry) {
           orientation = R.apply(orientation);
         }
       }
@@ -426,7 +428,7 @@ template <int D> class CiamarraSystem : public System<D, CiamarraState<D>> {
     }
   }
 
-  bool tryRandomMove(Rng& r) override {
+  bool randomMove(Rng& r) override {
     Vertex<D, CiamarraState<D>>& v = r.pick(this->vertices);
     CiamarraState<D> oldState = v.state;
 
@@ -457,9 +459,13 @@ template <int D> class CiamarraSystem : public System<D, CiamarraState<D>> {
     return false;
   }
 
-  bool swap(Vertex<D, CiamarraState<D>>& v1, Vertex<D, CiamarraState<D>>& v2) override {
+  void swap(Vertex<D, CiamarraState<D>>& v1, Vertex<D, CiamarraState<D>>& v2) override {
+    std::swap(v1.state, v2.state);
+  }
+
+  bool exchange(Vertex<D, CiamarraState<D>>& v1, Vertex<D, CiamarraState<D>>& v2) override {
     if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) {
-      std::swap(v1.state, v2.state);
+      swap(v1, v2);
       return true;
     } else {
       return false;