diff options
Diffstat (limited to 'transformation.hpp')
-rw-r--r-- | transformation.hpp | 252 |
1 files changed, 115 insertions, 137 deletions
diff --git a/transformation.hpp b/transformation.hpp index 717a033..b37b79e 100644 --- a/transformation.hpp +++ b/transformation.hpp @@ -9,26 +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(const R& r) : 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(double T, std::set<Spin<U, D, S>*>&, Spin<U, D, S>*, Rng& rng) const { return new Transformation(r); } + 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; Spin<U, D, S> sNew; public: - SpinFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s) : M(M), Transformation<U, D, R, S>(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); } @@ -36,13 +73,12 @@ 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(current_neighbors.begin(), current_neighbors.end()); - neighbors.insert(new_neighbors.begin(), new_neighbors.end()); - neighbors.insert(NULL); + neighbors.insert(oldNeighbors.begin(), oldNeighbors.end()); + neighbors.insert(newNeighbors.begin(), newNeighbors.end()); neighbors.erase(sOld); @@ -51,33 +87,40 @@ public: 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(double, std::set<Spin<U, D, S>*>&, Spin<U, D, S>* s, Rng&) 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; 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), Transformation<U, D, R, S>(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); @@ -87,17 +130,17 @@ public: std::set<Spin<U, D, S>*> current() const override { return {s1Old, s2Old}; } std::set<Spin<U, D, S>*> toConsider() const override { - std::set<Spin<U, D, S>*> neighbors; - 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>*> neighbors = {NULL}; + + 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.insert(NULL); neighbors.erase(s1Old); neighbors.erase(s2Old); @@ -107,137 +150,72 @@ public: 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(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const; - Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<signed, D, S>* s, Rng&) 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; - R s0New; - -public: - FieldFlip(Model<U, D, R, S>& M, const R& r) : M(M), Transformation<U, D, R, S>(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); - } - - 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(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const { - return new SpinFlip<U, D, R, S>(M, this->r, s); - } - - Transformation<signed, D, R, S>* createNew(double, std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const { - Vector<signed, 2> v = this->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, this->r, s); - } else { - return new PairFlip<U, D, R, S>(M, this->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(double T, std::set<Spin<U, D, S>*>& v, Spin<U, D, S>* s, Rng& rng) const { - if (s == NULL) { - return new FieldFlip<U, D, R, S>(M, this->r); - } else { - SpinFlip<U, D, R, S>* t = new SpinFlip<U, D, R, S>(M, this->r, s); - Spin<U, D, S>* sMax = NULL; - double ΔEMax = 0.0; - for (Spin<U, D, S>* ss : t->toConsider()) { - if (ss != NULL && ss != s && !v.contains(ss)) { - double ΔE = t->ΔE(ss); - if (ΔE > ΔEMax) { - sMax = ss; - ΔEMax = ΔE; - } +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; } } - if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1) || v.contains(sMax)) { - return t; - } else { - delete t; - v.insert(s); - v.insert(sMax); - return new PairFlip<U, D, R, S>(M, this->r, s, sMax); - } } -} - -template <class U, int D, class R, class S> -Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng& rng) const { - if (s == NULL) { - return new FieldFlip<U, D, R, S>(M, this->r); + if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1) || cluster.contains(sMax)) { + return t; } else { - SpinFlip<U, D, R, S>* t = new SpinFlip<U, D, R, S>(M, this->r, s); - Spin<U, D, S>* sMax = NULL; - double ΔEMax = 0.0; - for (Spin<U, D, S>* ss : t->toConsider()) { - if (ss != NULL && ss != s && !v.contains(ss)) { - double ΔE = t->ΔE(ss); - if (ΔE > ΔEMax) { - sMax = ss; - ΔEMax = ΔE; - } - } - } - if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1) || v.contains(sMax)) { - return t; - } else { - delete t; - v.insert(s); - v.insert(sMax); - return new PairFlip<U, D, R, S>(M, this->r, s, sMax); - } + 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(double, std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const { - if (s == NULL) { - return new FieldFlip<U, D, R, S>(M, this->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 { - Vector<signed, 2> v = this->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, this->r, s); - } else { - return new PairFlip<U, D, R, S>(M, this->r, s, *on_site.begin()); - } + cluster.insert(*on_site.begin()); + return new PairFlip<signed, D, R, S>(tOld->M, tOld->r, s, *on_site.begin()); } } + |