From 9437ed889390e04a1e693a5dfa5290cd90ae150e Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 14 Feb 2020 17:04:17 -0500 Subject: Dynamic pairing method appears to be working with spheres, needs to be cleaned up and ising code needs to be updated with new transformation interfaces --- spheres_infinite.cpp | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) (limited to 'spheres_infinite.cpp') diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp index 5287a84..b39f164 100644 --- a/spheres_infinite.cpp +++ b/spheres_infinite.cpp @@ -117,6 +117,28 @@ Gen, double> mGen(double ε) { }; } +Gen, double> rGen(double ε) { + return [ε](model& M, Rng& rng) -> Transformation, double>* { + Vector t; + Matrix m; + + double θ = rng.uniform((double)0.0, 2 * M_PI); + m(0, 0) = -cos(2 * θ); + m(1, 1) = cos(2 * θ); + m(0, 1) = -2 * cos(θ) * sin(θ); + m(1, 0) = -2 * cos(θ) * sin(θ); + + Spin* s = rng.pick(M.s); + t = M.s0.t; + for (unsigned j = 0; j < D; j++) { + t(j) += rng.variate(0.0, ε); + } + + Euclidean g(t - m * t, m); + return new SpinFlip, double>(M, g, s); + }; +} + int main(int argc, char* argv[]) { const unsigned D = 2; @@ -179,8 +201,9 @@ int main(int argc, char* argv[]) { return H * s.x.norm(); }; - auto g1 = eGen(1); + auto g1 = eGen(2); auto g2 = mGen(0.5); + auto g3 = rGen(5); animation A(L, 750, argc, argv); model sphere(1.0, Z, B); @@ -197,7 +220,7 @@ int main(int argc, char* argv[]) { sphere.dict.insert(ss); } - sphere.wolff(T, {g1, g2}, A, N); + sphere.wolff(T, {g1, g2, g3}, A, N); std::ofstream outfile; outfile.open("test.dat"); -- cgit v1.2.3-70-g09d2 From 2ebf2f181edac37bfb932dbb353101e37e97223a Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sun, 16 Feb 2020 21:02:25 -0500 Subject: More work on new transformations, and animation of cluster formation --- space_wolff.hpp | 8 +++--- spheres_infinite.cpp | 69 ++++++++++++++++++++++++++++++++++++++-------------- transformation.hpp | 58 ++++++++++++++++++++++++++++--------------- 3 files changed, 94 insertions(+), 41 deletions(-) (limited to 'spheres_infinite.cpp') diff --git a/space_wolff.hpp b/space_wolff.hpp index 63bf118..dc06011 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -18,11 +18,10 @@ template class Model; template class measurement { public: - virtual void pre_cluster(const Model&, unsigned, const R&){}; + virtual void pre_cluster(const Model&, unsigned, const Transformation, double>*){}; virtual void plain_bond_visited(const Model&, const Spin*, const Spin*, const Spin&, double){}; - virtual void plain_site_transformed(const Model&, const Spin*, - const Spin&){}; + virtual void plain_site_transformed(const Model&, const Transformation&){}; virtual void ghost_bond_visited(const Model&, const Spin&, const Spin&, double){}; @@ -70,6 +69,7 @@ public: } } + A.plain_site_transformed(*this, *t); t->apply(); } delete t; @@ -84,6 +84,8 @@ public: Gen& g = rng.pick(gs); Transformation *t = g(*this, rng); + A.pre_cluster(*this, i, t); + this->step(T, t, A); A.post_cluster(*this); diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp index b39f164..77a3582 100644 --- a/spheres_infinite.cpp +++ b/spheres_infinite.cpp @@ -14,9 +14,10 @@ private: uint64_t t2; unsigned n; unsigned tmp; + Euclidean s0_tmp; public: - animation(double L, unsigned w, int argc, char* argv[]) { + animation(double L, unsigned w, int argc, char* argv[]) : s0_tmp(0) { t1 = 0; t2 = 0; n = 0; @@ -31,30 +32,58 @@ public: gluOrtho2D(-L, L, -L, L); } - void pre_cluster(const model&, unsigned, const Euclidean&) override { tmp = 0; } - - void plain_site_transformed(const model&, const Spin*, - const Spin&) override { + void plain_site_transformed(const model& m, const Transformation, double>& t) override { + if (n > 5000) { + if (t.current().size() > 1) { + glColor3f(0.0f, 0.0f, 1.0f); + } else { + glColor3f(0.0f, 1.0f, 0.0f); + } + for (const Spin* s : t.current()) { + if (s != NULL) { + glBegin(GL_POLYGON); + unsigned n_points = 50; + for (unsigned i = 0; i < n_points; i++) { + glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points), + s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points)); + } + glEnd(); + } + } + } tmp++; } - void post_cluster(const model& m) override { + void pre_cluster(const model& m, unsigned, const Transformation, double>* t) override { + s0_tmp = m.s0.inverse(); + tmp = 0; glClearColor(1.0f, 1.0f, 1.0f, 1.0f); glClear(GL_COLOR_BUFFER_BIT); + /* + glBegin(GL_LINE); + Vector r_center = s0_tmp.act({t->r.t, 0}); + glEnd(); + */ for (const Spin* s : m.s) { glBegin(GL_POLYGON); unsigned n_points = 50; glColor3f(0.0f, 0.0f, 0.0f); for (unsigned i = 0; i < n_points; i++) { - glVertex2d(m.s0.inverse().act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points), - m.s0.inverse().act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points)); + glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points), + s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points)); } glEnd(); } glFlush(); + } + void post_cluster(const model& m) override { + glFlush(); t1 += tmp; t2 += tmp * tmp; + if (n > 5000) { + sleep(2); + } n++; } @@ -95,17 +124,21 @@ Gen, double> mGen(double ε) { Spin* s1 = rng.pick(M.s); Spin* s2 = rng.pick(M.s); - - while (s1 == s2) { - s2 = rng.pick(M.s); - } - Vector t1 = s1->x; Vector t2 = s2->x; Vector t12 = t1 - t2; + + while (s1 == s2 || 1 / t12.norm() < rng.uniform(0, 1)) { + s1 = rng.pick(M.s); + s2 = rng.pick(M.s); + t1 = s1->x; + t2 = s2->x; + t12 = t1 - t2; + } + Vector t = (t1 + t2) / 2; - double θ = atan2(t12[1], t12[0]) + rng.variate(0.0, ε); + double θ = atan2(t12[1], t12[0]) + rng.variate(0.0, ε) / t12.norm(); m(0, 0) = -cos(2 * θ); m(1, 1) = cos(2 * θ); @@ -201,9 +234,9 @@ int main(int argc, char* argv[]) { return H * s.x.norm(); }; - auto g1 = eGen(2); - auto g2 = mGen(0.5); - auto g3 = rGen(5); + auto g1 = eGen(0.5); + auto g2 = mGen(0.05); + auto g3 = rGen(1); animation A(L, 750, argc, argv); model sphere(1.0, Z, B); @@ -220,7 +253,7 @@ int main(int argc, char* argv[]) { sphere.dict.insert(ss); } - sphere.wolff(T, {g1, g2, g3}, A, N); + sphere.wolff(T, {g1, g2}, A, N); std::ofstream outfile; outfile.open("test.dat"); diff --git a/transformation.hpp b/transformation.hpp index ad3b49f..717a033 100644 --- a/transformation.hpp +++ b/transformation.hpp @@ -9,10 +9,12 @@ template class Model; template class Transformation { public: + const R r; + Transformation(const R& r) : r(r) {} virtual std::set*> current() const { return {NULL}; } virtual std::set*> toConsider() const { return {}; } virtual double ΔE(const Spin*) const { return 0; } - virtual Transformation* createNew(double T, std::set*>&, Spin*, Rng& rng) const { return new Transformation(); } + virtual Transformation* createNew(double T, std::set*>&, Spin*, Rng& rng) const { return new Transformation(r); } virtual void apply(){}; }; @@ -23,11 +25,10 @@ template class SpinFlip : public Transformati private: Model& M; Spin* sOld; - const R r; Spin sNew; public: - SpinFlip(Model& M, const R& r, Spin* s) : M(M), r(r) { + SpinFlip(Model& M, const R& r, Spin* s) : M(M), Transformation(r) { sOld = s; sNew = r.act(*s); } @@ -72,12 +73,11 @@ private: Model& M; Spin* s1Old; Spin* s2Old; - const R r; Spin s1New; Spin s2New; public: - PairFlip(Model& M, const R& r, Spin* s1, Spin* s2) : M(M), r(r) { + PairFlip(Model& M, const R& r, Spin* s1, Spin* s2) : M(M), Transformation(r) { s1Old = s1; s2Old = s2; s1New = r.act(*s1); @@ -133,11 +133,10 @@ public: template class FieldFlip : public Transformation { private: Model& M; - const R r; R s0New; public: - FieldFlip(Model& M, const R& r) : M(M), r(r), s0New(r.act(M.s0)) {} + FieldFlip(Model& M, const R& r) : M(M), Transformation(r), s0New(r.act(M.s0)) {} std::set*> toConsider() const override { std::set*> neighbors; @@ -156,16 +155,16 @@ public: } Transformation* createNew(double T, std::set*>& v, Spin* s, Rng&) const { - return new SpinFlip(M, r, s); + return new SpinFlip(M, this->r, s); } Transformation* createNew(double, std::set*>&, Spin* s, Rng&) const { - Vector v = r.act(*s).x; + Vector v = this->r.act(*s).x; std::set*> on_site = M.dict.at(v); if (on_site.empty()) { - return new SpinFlip(M, r, s); + return new SpinFlip(M, this->r, s); } else { - return new PairFlip(M, r, s, *on_site.begin()); + return new PairFlip(M, this->r, s, *on_site.begin()); } } @@ -175,18 +174,37 @@ public: template Transformation* SpinFlip::createNew(double T, std::set*>& v, Spin* s, Rng& rng) const { if (s == NULL) { - return new FieldFlip(M, r); + return new FieldFlip(M, this->r); } else { - return new SpinFlip(M, r, s); + SpinFlip* t = new SpinFlip(M, this->r, s); + Spin* sMax = NULL; + double ΔEMax = 0.0; + for (Spin* 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(0, 1) || v.contains(sMax)) { + return t; + } else { + delete t; + v.insert(s); + v.insert(sMax); + return new PairFlip(M, this->r, s, sMax); + } } } template Transformation* PairFlip::createNew(double T, std::set*>& v, Spin* s, Rng& rng) const { if (s == NULL) { - return new FieldFlip(M, r); + return new FieldFlip(M, this->r); } else { - SpinFlip* t = new SpinFlip(M, r, s); + SpinFlip* t = new SpinFlip(M, this->r, s); Spin* sMax = NULL; double ΔEMax = 0.0; for (Spin* ss : t->toConsider()) { @@ -204,7 +222,7 @@ Transformation* PairFlip::createNew(double T, std::set(M, r, s, sMax); + return new PairFlip(M, this->r, s, sMax); } } } @@ -212,14 +230,14 @@ Transformation* PairFlip::createNew(double T, std::set Transformation* PairFlip::createNew(double, std::set*>&, Spin* s, Rng&) const { if (s == NULL) { - return new FieldFlip(M, r); + return new FieldFlip(M, this->r); } else { - Vector v = r.act(*s).x; + Vector v = this->r.act(*s).x; std::set*> on_site = M.dict.at(v); if (on_site.empty()) { - return new SpinFlip(M, r, s); + return new SpinFlip(M, this->r, s); } else { - return new PairFlip(M, r, s, *on_site.begin()); + return new PairFlip(M, this->r, s, *on_site.begin()); } } } -- cgit v1.2.3-70-g09d2 From 3a3f2dd3b2c47d9d89ed29b7039e39626f2dcf72 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 18 Feb 2020 18:59:05 -0500 Subject: Cleaned up Transformation system and sped things up generally --- measurement.hpp | 20 ++++ space_wolff.hpp | 49 +++------- spheres_infinite.cpp | 38 ++++---- transformation.hpp | 252 +++++++++++++++++++++++---------------------------- 4 files changed, 169 insertions(+), 190 deletions(-) create mode 100644 measurement.hpp (limited to 'spheres_infinite.cpp') diff --git a/measurement.hpp b/measurement.hpp new file mode 100644 index 0000000..02fde9c --- /dev/null +++ b/measurement.hpp @@ -0,0 +1,20 @@ +#pragma once + +#include "spin.hpp" +#include "transformation.hpp" + +template class Model; + +template class measurement { +public: + virtual void pre_cluster(const Model&, unsigned, const Transformation*){}; + virtual void plain_bond_visited(const Model&, const Spin*, + const Spin*, const Spin&, double){}; + virtual void plain_site_transformed(const Model&, const Transformation&){}; + + virtual void ghost_bond_visited(const Model&, const Spin&, + const Spin&, double){}; + virtual void ghost_site_transformed(const Model&, const R&){}; + + virtual void post_cluster(const Model&){}; +}; diff --git a/space_wolff.hpp b/space_wolff.hpp index dc06011..3314ef9 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -8,31 +8,15 @@ #include #include "euclidean.hpp" +#include "measurement.hpp" #include "octree.hpp" #include "quantity.hpp" -#include "spin.hpp" #include "random.hpp" +#include "spin.hpp" #include "transformation.hpp" -template class Model; - -template class measurement { -public: - virtual void pre_cluster(const Model&, unsigned, const Transformation, double>*){}; - virtual void plain_bond_visited(const Model&, const Spin*, - const Spin*, const Spin&, double){}; - virtual void plain_site_transformed(const Model&, const Transformation&){}; - - virtual void ghost_bond_visited(const Model&, const Spin&, - const Spin&, double){}; - virtual void ghost_site_transformed(const Model&, const R&){}; - - virtual void post_cluster(const Model&){}; -}; - template class Model { public: - U L; R s0; std::vector*> s; Octree dict; @@ -42,47 +26,40 @@ public: Model(U L, std::function&, const Spin&)> Z, std::function&)> B) - : L(L), s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {} + : s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {} void step(double T, Transformation* t0, measurement& A) { std::queue*> queue; queue.push(t0); - std::set*> visited; - std::set*> fused; + std::set*> cluster = t0->current(); while (!queue.empty()) { Transformation* t = queue.front(); queue.pop(); - std::set*> c = t->current(); - if (!std::any_of(c.begin(), c.end(), [&visited](Spin* s) { return visited.contains(s); }) && !(c.size() == 1 && std::any_of(c.begin(), c.end(), [&fused](Spin* s) { return fused.contains(s); }))) { - visited.insert(c.begin(), c.end()); - fused.insert(c.begin(), c.end()); - - for (Spin* s : t->toConsider()) { - if (!fused.contains(s)) { + for (Spin* s : t->toConsider()) { + if (!cluster.contains(s)) { double ΔE = t->ΔE(s); if (ΔE > 0 && 1.0 - exp(-ΔE / T) > rng.uniform(0, 1)) { - queue.push(t->createNew(T, fused, s, rng)); - } + cluster.insert(s); + queue.push(t->createNew(T, cluster, s, rng)); } } - - A.plain_site_transformed(*this, *t); - t->apply(); } + + A.plain_site_transformed(*this, *t); + t->apply(); delete t; } } void resize(double T, double P) {} - void wolff(double T, std::vector> gs, - measurement& A, unsigned N) { + void wolff(double T, std::vector> gs, measurement& A, unsigned N) { for (unsigned i = 0; i < N; i++) { Gen& g = rng.pick(gs); - Transformation *t = g(*this, rng); + Transformation* t = g(*this, rng); A.pre_cluster(*this, i, t); diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp index 77a3582..e41ea17 100644 --- a/spheres_infinite.cpp +++ b/spheres_infinite.cpp @@ -14,10 +14,11 @@ private: uint64_t t2; unsigned n; unsigned tmp; + unsigned wait; Euclidean s0_tmp; public: - animation(double L, unsigned w, int argc, char* argv[]) : s0_tmp(0) { + animation(double L, unsigned w, int argc, char* argv[]) : s0_tmp(0), wait(1000) { t1 = 0; t2 = 0; n = 0; @@ -25,7 +26,7 @@ public: glutInit(&argc, argv); glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); glutInitWindowSize(w, w); - glutCreateWindow("wolff"); + glutCreateWindow("wolffWindow"); glClearColor(0.0, 0.0, 0.0, 0.0); glMatrixMode(GL_PROJECTION); glLoadIdentity(); @@ -33,7 +34,7 @@ public: } void plain_site_transformed(const model& m, const Transformation, double>& t) override { - if (n > 5000) { + if (n > wait) { if (t.current().size() > 1) { glColor3f(0.0f, 0.0f, 1.0f); } else { @@ -81,7 +82,7 @@ public: glFlush(); t1 += tmp; t2 += tmp * tmp; - if (n > 5000) { + if (n > wait) { sleep(2); } n++; @@ -124,28 +125,31 @@ Gen, double> mGen(double ε) { Spin* s1 = rng.pick(M.s); Spin* s2 = rng.pick(M.s); - Vector t1 = s1->x; - Vector t2 = s2->x; - Vector t12 = t1 - t2; - while (s1 == s2 || 1 / t12.norm() < rng.uniform(0, 1)) { - s1 = rng.pick(M.s); - s2 = rng.pick(M.s); - t1 = s1->x; - t2 = s2->x; - t12 = t1 - t2; + while (s1 == s2) { + s2 = rng.pick(M.s); } + Vector t1 = s1->x; + Vector t2 = s2->x; + Vector t12 = t1 - t2; Vector t = (t1 + t2) / 2; - double θ = atan2(t12[1], t12[0]) + rng.variate(0.0, ε) / t12.norm(); + double θ = atan2(t12(1), t12(0)) + rng.variate(0.0, ε) / t12.norm(); m(0, 0) = -cos(2 * θ); m(1, 1) = cos(2 * θ); m(0, 1) = -2 * cos(θ) * sin(θ); m(1, 0) = -2 * cos(θ) * sin(θ); - Euclidean g(t - m * t, m); + Vector t3 = t - m * t; + + if (t3(0) != t3(0)) { + std::cout << t3 << "\n" << t << "\n" << m << "\n" << t12 << "\n" ; + getchar(); + } + + Euclidean g(t3, m); return new PairFlip, double>(M, g, s1, s2); }; } @@ -234,8 +238,8 @@ int main(int argc, char* argv[]) { return H * s.x.norm(); }; - auto g1 = eGen(0.5); - auto g2 = mGen(0.05); + auto g1 = eGen(1); + auto g2 = mGen(0.1); auto g3 = rGen(1); animation A(L, 750, argc, argv); model sphere(1.0, Z, B); diff --git a/transformation.hpp b/transformation.hpp index 717a033..b37b79e 100644 --- a/transformation.hpp +++ b/transformation.hpp @@ -9,26 +9,63 @@ template class Model; template class Transformation { public: + Model& M; const R r; - Transformation(const R& r) : r(r) {} + + Transformation(Model& M, const R& r) : r(r), M(M) {} virtual std::set*> current() const { return {NULL}; } virtual std::set*> toConsider() const { return {}; } virtual double ΔE(const Spin*) const { return 0; } - virtual Transformation* createNew(double T, std::set*>&, Spin*, Rng& rng) const { return new Transformation(r); } + virtual Transformation* createNew(double T, std::set*>&, Spin*, + Rng& rng) const { + return new Transformation(M, r); + } virtual void apply(){}; }; template using Gen = std::function*(Model&, Rng&)>; +template class FieldFlip : public Transformation { +private: + R s0New; + +public: + FieldFlip(Model& M, const R& r) + : Transformation(M, r), s0New(r.act(M.s0)) {} + + std::set*> toConsider() const override { + std::set*> neighbors; + + for (Spin* s : this->M.s) { + neighbors.insert(s); + } + + return neighbors; + } + + double ΔE(const Spin* s) const override { + Spin s0s_old = this->M.s0.inverse().act(*s); + Spin s0s_new = s0New.inverse().act(*s); + return this->M.B(s0s_new) - this->M.B(s0s_old); + } + + Transformation* createNew(double T, std::set*>& cluster, + Spin* s, Rng& rng) const override { + return pairGenerator(this, T, cluster, s, rng); + } + + void apply() override { this->M.s0 = s0New; } +}; + template class SpinFlip : public Transformation { private: - Model& M; Spin* sOld; Spin sNew; public: - SpinFlip(Model& M, const R& r, Spin* s) : M(M), Transformation(r) { + SpinFlip(Model& M, const R& r, Spin* s) + : Transformation(M, r) { sOld = s; sNew = r.act(*s); } @@ -36,13 +73,12 @@ public: std::set*> current() const override { return {sOld}; } std::set*> toConsider() const override { - std::set*> neighbors; - std::set*> current_neighbors = M.dict.neighbors(sOld->x); - std::set*> new_neighbors = M.dict.neighbors(sNew.x); + std::set*> neighbors = {NULL}; + std::set*> oldNeighbors = this->M.dict.neighbors(sOld->x); + std::set*> 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* s) const override { if (s == NULL) { - Spin s0s_old = M.s0.inverse().act(*sOld); - Spin s0s_new = M.s0.inverse().act(sNew); - return M.B(s0s_new) - M.B(s0s_old); + Spin s0s_old = this->M.s0.inverse().act(*sOld); + Spin 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* createNew(double, std::set*>&, Spin* s, Rng&) const override; + Transformation* createNew(double T, std::set*>& cluster, + Spin* s, Rng& rng) const override { + if (s == NULL) { + return new FieldFlip(this->M, this->r); + } else { + return pairGenerator(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 PairFlip : public Transformation { private: - Model& M; Spin* s1Old; Spin* s2Old; Spin s1New; Spin s2New; public: - PairFlip(Model& M, const R& r, Spin* s1, Spin* s2) : M(M), Transformation(r) { + PairFlip(Model& M, const R& r, Spin* s1, Spin* s2) + : Transformation(M, r) { s1Old = s1; s2Old = s2; s1New = r.act(*s1); @@ -87,17 +130,17 @@ public: std::set*> current() const override { return {s1Old, s2Old}; } std::set*> toConsider() const override { - std::set*> neighbors; - std::set*> current_neighbors_1 = M.dict.neighbors(s1Old->x); - std::set*> current_neighbors_2 = M.dict.neighbors(s2Old->x); - std::set*> new_neighbors_1 = M.dict.neighbors(s1New.x); - std::set*> new_neighbors_2 = M.dict.neighbors(s2New.x); + std::set*> neighbors = {NULL}; + + std::set*> current_neighbors_1 = this->M.dict.neighbors(s1Old->x); + std::set*> current_neighbors_2 = this->M.dict.neighbors(s2Old->x); + std::set*> new_neighbors_1 = this->M.dict.neighbors(s1New.x); + std::set*> 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* s) const override { if (s == NULL) { - Spin s0s1_old = M.s0.inverse().act(*s1Old); - Spin s0s1_new = M.s0.inverse().act(s1New); - Spin s0s2_old = M.s0.inverse().act(*s2Old); - Spin 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 s0s1_old = this->M.s0.inverse().act(*s1Old); + Spin s0s1_new = this->M.s0.inverse().act(s1New); + Spin s0s2_old = this->M.s0.inverse().act(*s2Old); + Spin 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* createNew(double T, std::set*>& v, Spin* s, Rng&) const; - Transformation* createNew(double T, std::set*>& v, Spin* s, Rng&) const; + Transformation* createNew(double T, std::set*>& cluster, + Spin* s, Rng& rng) const override { + if (s == NULL) { + return new FieldFlip(this->M, this->r); + } else { + return pairGenerator(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 FieldFlip : public Transformation { -private: - Model& M; - R s0New; - -public: - FieldFlip(Model& M, const R& r) : M(M), Transformation(r), s0New(r.act(M.s0)) {} - - std::set*> toConsider() const override { - std::set*> neighbors; - - for (Spin* s : M.s) { - neighbors.insert(s); - } - - return neighbors; - } - - double ΔE(const Spin* s) const override { - Spin s0s_old = M.s0.inverse().act(*s); - Spin s0s_new = s0New.inverse().act(*s); - return M.B(s0s_new) - M.B(s0s_old); - } - - Transformation* createNew(double T, std::set*>& v, Spin* s, Rng&) const { - return new SpinFlip(M, this->r, s); - } - - Transformation* createNew(double, std::set*>&, Spin* s, Rng&) const { - Vector v = this->r.act(*s).x; - std::set*> on_site = M.dict.at(v); - if (on_site.empty()) { - return new SpinFlip(M, this->r, s); - } else { - return new PairFlip(M, this->r, s, *on_site.begin()); - } - } - - void apply() override { M.s0 = s0New; } -}; -template -Transformation* SpinFlip::createNew(double T, std::set*>& v, Spin* s, Rng& rng) const { - if (s == NULL) { - return new FieldFlip(M, this->r); - } else { - SpinFlip* t = new SpinFlip(M, this->r, s); - Spin* sMax = NULL; - double ΔEMax = 0.0; - for (Spin* ss : t->toConsider()) { - if (ss != NULL && ss != s && !v.contains(ss)) { - double ΔE = t->ΔE(ss); - if (ΔE > ΔEMax) { - sMax = ss; - ΔEMax = ΔE; - } +template +Transformation* +pairGenerator(const Transformation* tOld, double T, + std::set*>& cluster, Spin* s, Rng& rng) { + SpinFlip* t = new SpinFlip(tOld->M, tOld->r, s); + Spin* sMax = NULL; + double ΔEMax = 0.0; + for (Spin* 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(0, 1) || v.contains(sMax)) { - return t; - } else { - delete t; - v.insert(s); - v.insert(sMax); - return new PairFlip(M, this->r, s, sMax); - } } -} - -template -Transformation* PairFlip::createNew(double T, std::set*>& v, Spin* s, Rng& rng) const { - if (s == NULL) { - return new FieldFlip(M, this->r); + if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform(0, 1) || cluster.contains(sMax)) { + return t; } else { - SpinFlip* t = new SpinFlip(M, this->r, s); - Spin* sMax = NULL; - double ΔEMax = 0.0; - for (Spin* 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(0, 1) || v.contains(sMax)) { - return t; - } else { - delete t; - v.insert(s); - v.insert(sMax); - return new PairFlip(M, this->r, s, sMax); - } + delete t; + cluster.insert(sMax); + return new PairFlip(tOld->M, tOld->r, s, sMax); } } -template -Transformation* PairFlip::createNew(double, std::set*>&, Spin* s, Rng&) const { - if (s == NULL) { - return new FieldFlip(M, this->r); +template +Transformation* +pairGenerator(const Transformation* tOld, double T, + std::set*>& cluster, Spin* s, Rng& rng) { + Vector v = tOld->r.act(*s).x; + std::set*> on_site = tOld->M.dict.at(v); + if (on_site.empty()) { + return new SpinFlip(tOld->M, tOld->r, s); } else { - Vector v = this->r.act(*s).x; - std::set*> on_site = M.dict.at(v); - if (on_site.empty()) { - return new SpinFlip(M, this->r, s); - } else { - return new PairFlip(M, this->r, s, *on_site.begin()); - } + cluster.insert(*on_site.begin()); + return new PairFlip(tOld->M, tOld->r, s, *on_site.begin()); } } + -- cgit v1.2.3-70-g09d2