diff options
-rw-r--r-- | measurement.hpp | 20 | ||||
-rw-r--r-- | space_wolff.hpp | 47 | ||||
-rw-r--r-- | spheres_infinite.cpp | 84 | ||||
-rw-r--r-- | transformation.hpp | 219 |
4 files changed, 227 insertions, 143 deletions
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 U, int D, class R, class S> class Model; + +template <class U, int D, class R, class S> class measurement { +public: + virtual void pre_cluster(const Model<U, D, R, S>&, unsigned, const Transformation<U, D, R, S>*){}; + virtual void plain_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>*, + const Spin<U, D, S>*, const Spin<U, D, S>&, double){}; + virtual void plain_site_transformed(const Model<U, D, R, S>&, const Transformation<U, D, R, S>&){}; + + virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&, + const Spin<U, D, S>&, double){}; + virtual void ghost_site_transformed(const Model<U, D, R, S>&, const R&){}; + + virtual void post_cluster(const Model<U, D, R, S>&){}; +}; diff --git a/space_wolff.hpp b/space_wolff.hpp index 9d529a3..3314ef9 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -8,32 +8,15 @@ #include <vector> #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 U, int D, class R, class S> class Model; - -template <class U, int D, class R, class S> class measurement { -public: - virtual void pre_cluster(const Model<U, D, R, S>&, unsigned, const R&){}; - virtual void plain_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>*, - const Spin<U, D, S>*, const Spin<U, D, S>&, double){}; - virtual void plain_site_transformed(const Model<U, D, R, S>&, const Spin<U, D, S>*, - const Spin<U, D, S>&){}; - - virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&, - const Spin<U, D, S>&, double){}; - virtual void ghost_site_transformed(const Model<U, D, R, S>&, const R&){}; - - virtual void post_cluster(const Model<U, D, R, S>&){}; -}; - template <class U, int D, class R, class S> class Model { public: - U L; R s0; std::vector<Spin<U, D, S>*> s; Octree<U, D, S> dict; @@ -43,42 +26,42 @@ public: Model(U L, std::function<double(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z, std::function<double(const Spin<U, D, S>&)> 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<U, D, R, S>* t0, measurement<U, D, R, S>& A) { std::queue<Transformation<U, D, R, S>*> queue; queue.push(t0); - std::set<Spin<U, D, S>*> visited; + std::set<Spin<U, D, S>*> cluster = t0->current(); while (!queue.empty()) { Transformation<U, D, R, S>* t = queue.front(); queue.pop(); - std::set<Spin<U, D, S>*> c = t->current(); - if (!std::any_of(c.begin(), c.end(), [&visited](Spin<U, D, S>* s) { return visited.contains(s); })) { - visited.insert(c.begin(), c.end()); - - for (Spin<U, D, S>* s : t->toConsider()) { + for (Spin<U, D, S>* s : t->toConsider()) { + if (!cluster.contains(s)) { double ΔE = t->ΔE(s); if (ΔE > 0 && 1.0 - exp(-ΔE / T) > rng.uniform<double>(0, 1)) { - queue.push(t->createNew(s)); + cluster.insert(s); + queue.push(t->createNew(T, cluster, s, rng)); } } - - t->apply(); } + + A.plain_site_transformed(*this, *t); + t->apply(); delete t; } } void resize(double T, double P) {} - void wolff(double T, std::vector<Gen<U, D, R, S>> gs, - measurement<U, D, R, S>& A, unsigned N) { + void wolff(double T, std::vector<Gen<U, D, R, S>> gs, measurement<U, D, R, S>& A, unsigned N) { for (unsigned i = 0; i < N; i++) { Gen<U, D, R, S>& g = rng.pick(gs); - Transformation<U, D, R, S> *t = g(*this, rng); + Transformation<U, D, R, S>* t = g(*this, rng); + + A.pre_cluster(*this, i, t); this->step(T, t, A); diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp index 5287a84..e41ea17 100644 --- a/spheres_infinite.cpp +++ b/spheres_infinite.cpp @@ -14,9 +14,11 @@ private: uint64_t t2; unsigned n; unsigned tmp; + unsigned wait; + Euclidean<double, D> s0_tmp; public: - animation(double L, unsigned w, int argc, char* argv[]) { + animation(double L, unsigned w, int argc, char* argv[]) : s0_tmp(0), wait(1000) { t1 = 0; t2 = 0; n = 0; @@ -24,37 +26,65 @@ 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(); gluOrtho2D(-L, L, -L, L); } - void pre_cluster(const model&, unsigned, const Euclidean<double, D>&) override { tmp = 0; } - - void plain_site_transformed(const model&, const Spin<double, D, double>*, - const Spin<double, D, double>&) override { + void plain_site_transformed(const model& m, const Transformation<double, D, Euclidean<double, D>, double>& t) override { + if (n > wait) { + if (t.current().size() > 1) { + glColor3f(0.0f, 0.0f, 1.0f); + } else { + glColor3f(0.0f, 1.0f, 0.0f); + } + for (const Spin<double, 2, double>* 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, D, Euclidean<double, D>, 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<double, 2> r_center = s0_tmp.act({t->r.t, 0}); + glEnd(); + */ for (const Spin<double, 2, double>* 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 > wait) { + sleep(2); + } n++; } @@ -105,18 +135,47 @@ Gen<double, D, Euclidean<double, D>, double> mGen(double ε) { Vector<double, D> t12 = t1 - t2; Vector<double, D> t = (t1 + t2) / 2; - double θ = atan2(t12[1], t12[0]) + rng.variate<double, std::normal_distribution>(0.0, ε); + double θ = atan2(t12(1), t12(0)) + rng.variate<double, std::normal_distribution>(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<double, D> g(t - m * t, m); + Vector<double, D> t3 = t - m * t; + + if (t3(0) != t3(0)) { + std::cout << t3 << "\n" << t << "\n" << m << "\n" << t12 << "\n" ; + getchar(); + } + + Euclidean<double, D> g(t3, m); return new PairFlip<double, D, Euclidean<double, D>, double>(M, g, s1, s2); }; } +Gen<double, D, Euclidean<double, D>, double> rGen(double ε) { + return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* { + Vector<double, D> t; + Matrix<double, D> 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<double, D, double>* s = rng.pick(M.s); + t = M.s0.t; + for (unsigned j = 0; j < D; j++) { + t(j) += rng.variate<double, std::normal_distribution>(0.0, ε); + } + + Euclidean<double, D> g(t - m * t, m); + return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, s); + }; +} + int main(int argc, char* argv[]) { const unsigned D = 2; @@ -180,7 +239,8 @@ int main(int argc, char* argv[]) { }; auto g1 = eGen(1); - auto g2 = mGen(0.5); + 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 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()); - } - } -} |