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 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-54-g00ecf