Merge branch 'pair_generation'

1 files changed, 120 insertions, 99 deletions

diff --git a/transformation.hpp b/transformation.hpp
index 33ac50f..b37b79e 100644
--- a/transformation.hpp
+++ b/transformation.hpp
@@ -9,25 +9,63 @@ template <class U, int D, class R, class S> class Model;
 
 template <class U, int D, class R, class S> class Transformation {
 public:
+  Model<U, D, R, S>& M;
+  const R r;
+
+  Transformation(Model<U, D, R, S>& M, const R& r) : r(r), M(M) {}
   virtual std::set<Spin<U, D, S>*> current() const { return {NULL}; }
   virtual std::set<Spin<U, D, S>*> toConsider() const { return {}; }
   virtual double ΔE(const Spin<U, D, S>*) const { return 0; }
-  virtual Transformation* createNew(Spin<U, D, S>*) const { return new Transformation(); }
+  virtual Transformation* createNew(double T, std::set<Spin<U, D, S>*>&, Spin<U, D, S>*,
+                                    Rng& rng) const {
+    return new Transformation(M, r);
+  }
   virtual void apply(){};
 };
 
 template <class U, int D, class R, class S>
 using Gen = std::function<Transformation<U, D, R, S>*(Model<U, D, R, S>&, Rng&)>;
 
+template <class U, int D, class R, class S> class FieldFlip : public Transformation<U, D, R, S> {
+private:
+  R s0New;
+
+public:
+  FieldFlip(Model<U, D, R, S>& M, const R& r)
+      : Transformation<U, D, R, S>(M, r), s0New(r.act(M.s0)) {}
+
+  std::set<Spin<U, D, S>*> toConsider() const override {
+    std::set<Spin<U, D, S>*> neighbors;
+
+    for (Spin<U, D, S>* s : this->M.s) {
+      neighbors.insert(s);
+    }
+
+    return neighbors;
+  }
+
+  double ΔE(const Spin<U, D, S>* s) const override {
+    Spin<U, D, S> s0s_old = this->M.s0.inverse().act(*s);
+    Spin<U, D, S> s0s_new = s0New.inverse().act(*s);
+    return this->M.B(s0s_new) - this->M.B(s0s_old);
+  }
+
+  Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& cluster,
+                                             Spin<U, D, S>* s, Rng& rng) const override {
+    return pairGenerator<D, R, S>(this, T, cluster, s, rng);
+  }
+
+  void apply() override { this->M.s0 = s0New; }
+};
+
 template <class U, int D, class R, class S> class SpinFlip : public Transformation<U, D, R, S> {
 private:
-  Model<U, D, R, S>& M;
   Spin<U, D, S>* sOld;
-  const R r;
   Spin<U, D, S> sNew;
 
 public:
-  SpinFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s) : M(M), r(r) {
+  SpinFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s)
+      : Transformation<U, D, R, S>(M, r) {
     sOld = s;
     sNew = r.act(*s);
   }
@@ -35,47 +73,54 @@ public:
   std::set<Spin<U, D, S>*> current() const override { return {sOld}; }
 
   std::set<Spin<U, D, S>*> toConsider() const override {
-    std::set<Spin<U, D, S>*> neighbors;
-    std::set<Spin<U, D, S>*> current_neighbors = M.dict.neighbors(sOld->x);
-    std::set<Spin<U, D, S>*> new_neighbors = M.dict.neighbors(sNew.x);
+    std::set<Spin<U, D, S>*> neighbors = {NULL};
+    std::set<Spin<U, D, S>*> oldNeighbors = this->M.dict.neighbors(sOld->x);
+    std::set<Spin<U, D, S>*> newNeighbors = this->M.dict.neighbors(sNew.x);
+
+    neighbors.insert(oldNeighbors.begin(), oldNeighbors.end());
+    neighbors.insert(newNeighbors.begin(), newNeighbors.end());
 
-    neighbors.insert(current_neighbors.begin(), current_neighbors.end());
-    neighbors.insert(new_neighbors.begin(), new_neighbors.end());
-    neighbors.insert(NULL);
+    neighbors.erase(sOld);
 
     return neighbors;
   }
 
   double ΔE(const Spin<U, D, S>* s) const override {
     if (s == NULL) {
-      Spin<U, D, S> s0s_old = M.s0.inverse().act(*sOld);
-      Spin<U, D, S> s0s_new = M.s0.inverse().act(sNew);
-      return M.B(s0s_new) - M.B(s0s_old);
+      Spin<U, D, S> s0s_old = this->M.s0.inverse().act(*sOld);
+      Spin<U, D, S> s0s_new = this->M.s0.inverse().act(sNew);
+      return this->M.B(s0s_new) - this->M.B(s0s_old);
     } else {
-      return M.Z(*sOld, *s) - M.Z(sNew, *s);
+      return this->M.Z(*sOld, *s) - this->M.Z(sNew, *s);
     }
   }
 
-  Transformation<U, D, R, S>* createNew(Spin<U, D, S>* s) const override;
+  Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& cluster,
+                                             Spin<U, D, S>* s, Rng& rng) const override {
+    if (s == NULL) {
+      return new FieldFlip<U, D, R, S>(this->M, this->r);
+    } else {
+      return pairGenerator<D, R, S>(this, T, cluster, s, rng);
+    }
+  }
 
   void apply() override {
-    M.dict.remove(sOld);
+    this->M.dict.remove(sOld);
     *sOld = sNew;
-    M.dict.insert(sOld);
+    this->M.dict.insert(sOld);
   }
 };
 
 template <class U, int D, class R, class S> class PairFlip : public Transformation<U, D, R, S> {
 private:
-  Model<U, D, R, S>& M;
   Spin<U, D, S>* s1Old;
   Spin<U, D, S>* s2Old;
-  const R r;
   Spin<U, D, S> s1New;
   Spin<U, D, S> s2New;
 
 public:
-  PairFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s1, Spin<U, D, S>* s2) : M(M), r(r) {
+  PairFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s1, Spin<U, D, S>* s2)
+      : Transformation<U, D, R, S>(M, r) {
     s1Old = s1;
     s2Old = s2;
     s1New = r.act(*s1);
@@ -86,115 +131,91 @@ public:
 
   std::set<Spin<U, D, S>*> toConsider() const override {
     std::set<Spin<U, D, S>*> neighbors = {NULL};
-    std::set<Spin<U, D, S>*> current_neighbors_1 = M.dict.neighbors(s1Old->x);
-    std::set<Spin<U, D, S>*> current_neighbors_2 = M.dict.neighbors(s2Old->x);
-    std::set<Spin<U, D, S>*> new_neighbors_1 = M.dict.neighbors(s1New.x);
-    std::set<Spin<U, D, S>*> new_neighbors_2 = M.dict.neighbors(s2New.x);
+
+    std::set<Spin<U, D, S>*> current_neighbors_1 = this->M.dict.neighbors(s1Old->x);
+    std::set<Spin<U, D, S>*> current_neighbors_2 = this->M.dict.neighbors(s2Old->x);
+    std::set<Spin<U, D, S>*> new_neighbors_1 = this->M.dict.neighbors(s1New.x);
+    std::set<Spin<U, D, S>*> new_neighbors_2 = this->M.dict.neighbors(s2New.x);
 
     neighbors.insert(current_neighbors_1.begin(), current_neighbors_1.end());
     neighbors.insert(current_neighbors_2.begin(), current_neighbors_2.end());
     neighbors.insert(new_neighbors_1.begin(), new_neighbors_1.end());
     neighbors.insert(new_neighbors_2.begin(), new_neighbors_2.end());
 
+    neighbors.erase(s1Old);
+    neighbors.erase(s2Old);
+
     return neighbors;
   }
 
   double ΔE(const Spin<U, D, S>* s) const override {
     if (s == NULL) {
-      Spin<U, D, S> s0s1_old = M.s0.inverse().act(*s1Old);
-      Spin<U, D, S> s0s1_new = M.s0.inverse().act(s1New);
-      Spin<U, D, S> s0s2_old = M.s0.inverse().act(*s2Old);
-      Spin<U, D, S> s0s2_new = M.s0.inverse().act(s2New);
-      return M.B(s0s1_new) + M.B(s0s2_new) - M.B(s0s1_old) - M.B(s0s2_old);
+      Spin<U, D, S> s0s1_old = this->M.s0.inverse().act(*s1Old);
+      Spin<U, D, S> s0s1_new = this->M.s0.inverse().act(s1New);
+      Spin<U, D, S> s0s2_old = this->M.s0.inverse().act(*s2Old);
+      Spin<U, D, S> s0s2_new = this->M.s0.inverse().act(s2New);
+      return this->M.B(s0s1_new) + this->M.B(s0s2_new) - this->M.B(s0s1_old) - this->M.B(s0s2_old);
     } else {
-      return M.Z(*s1Old, *s) + M.Z(*s2Old, *s) - M.Z(s1New, *s) - M.Z(s2New, *s);
+      return this->M.Z(*s1Old, *s) + this->M.Z(*s2Old, *s) - this->M.Z(s1New, *s) - this->M.Z(s2New, *s);
     }
   }
 
-  Transformation<U, D, R, S>* createNew(Spin<double, D, S>* s) const;
-  Transformation<U, D, R, S>* createNew(Spin<signed, D, S>* s) const;
+  Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& cluster,
+                                             Spin<U, D, S>* s, Rng& rng) const override {
+    if (s == NULL) {
+      return new FieldFlip<U, D, R, S>(this->M, this->r);
+    } else {
+      return pairGenerator<D, R, S>(this, T, cluster, s, rng);
+    }
+  }
 
   void apply() override {
-    M.dict.remove(s1Old);
-    M.dict.remove(s2Old);
+    this->M.dict.remove(s1Old);
+    this->M.dict.remove(s2Old);
     *s1Old = s1New;
     *s2Old = s2New;
-    M.dict.insert(s1Old);
-    M.dict.insert(s2Old);
+    this->M.dict.insert(s1Old);
+    this->M.dict.insert(s2Old);
   }
 };
 
-template <class U, int D, class R, class S> class FieldFlip : public Transformation<U, D, R, S> {
-private:
-  Model<U, D, R, S>& M;
-  const R r;
-  R s0New;
-
-public:
-  FieldFlip(Model<U, D, R, S>& M, const R& r) : M(M), r(r), s0New(r.act(M.s0)) {}
-
-  std::set<Spin<U, D, S>*> toConsider() const override {
-    std::set<Spin<U, D, S>*> neighbors;
 
-    for (Spin<U, D, S>* s : M.s) {
-      neighbors.insert(s);
+template <int D, class R, class S>
+Transformation<double, D, R, S>*
+pairGenerator(const Transformation<double, D, R, S>* tOld, double T,
+                        std::set<Spin<double, D, S>*>& cluster, Spin<double, D, S>* s, Rng& rng) {
+  SpinFlip<double, D, R, S>* t = new SpinFlip<double, D, R, S>(tOld->M, tOld->r, s);
+  Spin<double, D, S>* sMax = NULL;
+  double ΔEMax = 0.0;
+  for (Spin<double, D, S>* ss : t->toConsider()) {
+    if (ss != NULL && ss != s && !cluster.contains(ss)) {
+      double ΔE = t->ΔE(ss);
+      if (ΔE > ΔEMax) {
+        sMax = ss;
+        ΔEMax = ΔE;
+      }
     }
-
-    return neighbors;
   }
-
-  double ΔE(const Spin<U, D, S>* s) const override {
-    Spin<U, D, S> s0s_old = M.s0.inverse().act(*s);
-    Spin<U, D, S> s0s_new = s0New.inverse().act(*s);
-    return M.B(s0s_new) - M.B(s0s_old);
-  }
-
-  Transformation<double, D, R, S>* createNew(Spin<double, D, S>* s) const {
-    return new SpinFlip<U, D, R, S>(M, r, s);
-  }
-
-  Transformation<signed, D, R, S>* createNew(Spin<signed, D, S>* s) const {
-    Vector<signed, 2> v = r.act(*s).x;
-    std::set<Spin<U, D, S>*> on_site = M.dict.at(v);
-    if (on_site.empty()) {
-      return new SpinFlip<U, D, R, S>(M, r, s);
-    } else {
-      return new PairFlip<U, D, R, S>(M, r, s, *on_site.begin());
-    }
-  }
-
-  void apply() override { M.s0 = s0New; }
-};
-
-template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(Spin<U, D, S>* s) const {
-  if (s == NULL) {
-    return new FieldFlip<U, D, R, S>(M, r);
+  if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1) || cluster.contains(sMax)) {
+    return t;
   } else {
-    return new SpinFlip(M, r, s);
+    delete t;
+    cluster.insert(sMax);
+    return new PairFlip<double, D, R, S>(tOld->M, tOld->r, s, sMax);
   }
 }
 
-template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(Spin<double, D, S>* s) const {
-  if (s == NULL) {
-    return new FieldFlip<U, D, R, S>(M, r);
+template <int D, class R, class S>
+Transformation<signed, D, R, S>*
+pairGenerator(const Transformation<signed, D, R, S>* tOld, double T,
+                        std::set<Spin<signed, D, S>*>& cluster, Spin<signed, D, S>* s, Rng& rng) {
+  Vector<signed, 2> v = tOld->r.act(*s).x;
+  std::set<Spin<signed, D, S>*> on_site = tOld->M.dict.at(v);
+  if (on_site.empty()) {
+    return new SpinFlip<signed, D, R, S>(tOld->M, tOld->r, s);
   } else {
-    return new SpinFlip<U, D, R, S>(M, r, s);
+    cluster.insert(*on_site.begin());
+    return new PairFlip<signed, D, R, S>(tOld->M, tOld->r, s, *on_site.begin());
   }
 }
 
-template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(Spin<signed, D, S>* s) const {
-  if (s == NULL) {
-    return new FieldFlip<U, D, R, S>(M, r);
-  } else {
-    Vector<signed, 2> v = r.act(*s).x;
-    std::set<Spin<U, D, S>*> on_site = M.dict.at(v);
-    if (on_site.empty()) {
-      return new SpinFlip<U, D, R, S>(M, r, s);
-    } else {
-      return new PairFlip<U, D, R, S>(M, r, s, *on_site.begin());
-    }
-  }
-}