From d5a7a3f2e5808e7a3327242dc9b368afed9383bf Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 12 Feb 2020 14:59:13 -0500 Subject: Changed the RNG and implemented Transformations as the method for parameterizing movement. --- .gitmodules | 3 + pcg-cpp | 1 + space_wolff.hpp | 241 +++++++++++++++++++++++++++++++++++---------------- spheres_infinite.cpp | 30 ++++--- 4 files changed, 190 insertions(+), 85 deletions(-) create mode 160000 pcg-cpp diff --git a/.gitmodules b/.gitmodules index 0a97926..15de333 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,6 @@ [submodule "randutils"] path = randutils url = https://gist.github.com/imneme/540829265469e673d045/ +[submodule "pcg-cpp"] + path = pcg-cpp + url = https://github.com/imneme/pcg-cpp diff --git a/pcg-cpp b/pcg-cpp new file mode 160000 index 0000000..5b5cac8 --- /dev/null +++ b/pcg-cpp @@ -0,0 +1 @@ +Subproject commit 5b5cac8d61339e810c5dbb4692d868a1d7ca1b2d diff --git a/space_wolff.hpp b/space_wolff.hpp index 4c935db..8ea7dfb 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -7,6 +7,7 @@ #include #include +#include "pcg-cpp/include/pcg_random.hpp" #include "randutils/randutils.hpp" #include "euclidean.hpp" @@ -14,6 +15,8 @@ #include "quantity.hpp" #include "spin.hpp" +using Rng = randutils::random_generator; + template class Model; template class measurement { @@ -31,42 +34,171 @@ public: virtual void post_cluster(const Model&){}; }; -template using Gen = std::function&, randutils::mt19937_rng&)>; - -template class Model; template class Transformation { public: - virtual void ready() {}; - virtual double p(const Transformation&) const { return 0; } - virtual void apply() {}; - virtual std::set to_consider() const {}; + virtual std::set*> current() const {return {NULL};} + virtual std::set*> toConsider() const {return {};} + virtual double ΔE(const Spin*) const {return 0;} + virtual Transformation* createNew(Spin*) const {return new Transformation();} + virtual void apply(); }; +template using Gen = std::function*(Model&, Rng&)>; + template class SpinFlip : public Transformation { private: Model& M; - Spin& s; - const R& r; - Spin s_new; + Spin& sOld; + const R r; + Spin sNew; + public: + SpinFlip(Model& M, const R& r, Spin& s) : M(M), r(r), sOld(s) { + sNew = r.act(s); + } + + 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); + + neighbors.insert(current_neighbors.begin(), current_neighbors.end()); + neighbors.insert(new_neighbors.begin(), new_neighbors.end()); + neighbors.insert(NULL); + + return neighbors; + } + + 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); + } else { + return M.Z(sOld, *s) - M.Z(sNew, *s); + } + } + + Transformation* createNew(Spin* s) const override; + + void apply() override { + M.dict.remove(&sOld); + sOld = sNew; + M.dict.insert(&sOld); + } +}; + +template class PairFlip : public Transformation { + 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), s1Old(s1), s2Old(s2) { + s1New = r.act(s1); + s2New = r.act(s2); + } + + 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); + + 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); + + return neighbors; + } + + 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); + } else { + return M.Z(s1Old, *s) + M.Z(s2Old, *s) - M.Z(s1New, *s) - M.Z(s2New, *s); + } + } + + Transformation* createNew(Spin* s) const override; + + void apply() override { + M.dict.remove(&s1Old); + M.dict.remove(&s2Old); + s1Old = s1New; + s2Old = s2New; + M.dict.insert(&s1Old); + M.dict.insert(&s2Old); + } +}; + +template class FieldFlip : public Transformation { + private: + Model& M; + const R r; + R s0New; public: - SpinFlip(Model& M, Spin& s, const R& r) : M(M), s(s), r(r) {} + FieldFlip(Model& M, const R& r) : M(M), r(r), s0New(r.act(M.s0)) {} + + std::set*> toConsider() const override { + std::set*> neighbors; - void ready() override { - s_new = r.act(s); - }; + for (Spin& s : M.s) { + neighbors.insert(&s); + } - double p(const SpinFlip& sf, double T) const override { - return 1.0 - exp(-(M.Z(s, sf.s) - M.Z(s_new, sf.s)) / T); + return neighbors; } - double p(const SpinFlip& sf, double T) const override { - return 1.0 - exp(-(M.Z(s, sf.s) - M.Z(s_new, sf.s)) / T); + 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(Spin* s) const override { + return new SpinFlip(M, r, *s); + } + void apply() override { + M.s0 = s0New; + } }; +template Transformation* SpinFlip::createNew(Spin* s) const { + if (s == NULL) { + return new FieldFlip(M, r); + } else { + return new SpinFlip(M, r, *s); + } +} + +template Transformation* PairFlip::createNew(Spin* s) const { + if (s == NULL) { + return new FieldFlip(M, r); + } else { + return new SpinFlip(M, r, *s); + } +} + template class Model { public: U L; @@ -75,70 +207,36 @@ public: Octree dict; std::function&, const Spin&)> Z; std::function&)> B; - randutils::mt19937_rng rng; + Rng rng; 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)) {} - void step(double T, unsigned ind, R& r, measurement& A) { - std::queue*> queue; - queue.push(&(s[ind])); + void step(double T, Transformation* t0, measurement& A) { + std::queue*> queue; + queue.push(t0); std::set*> visited; while (!queue.empty()) { - Spin* si = queue.front(); + Transformation* t = queue.front(); queue.pop(); - if (!visited.contains(si)) { - visited.insert(si); - - if (si == NULL) { - R s0_new = r.act(s0); - for (Spin& ss : s) { - Spin s0s_old = s0.inverse().act(ss); - Spin s0s_new = s0_new.inverse().act(ss); - double p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T); - A.ghost_bond_visited(*this, s0s_old, s0s_new, p); - if (rng.uniform(0.0, 1.0) < p) { - queue.push(&ss); - } - } - A.ghost_site_transformed(*this, s0_new); - s0 = s0_new; - } else { - Spin si_new = r.act(*si); - std::set*> all_neighbors; - std::set*> current_neighbors = dict.neighbors(si->x); - std::set*> new_neighbors = dict.neighbors(si_new.x); - - all_neighbors.insert(current_neighbors.begin(), current_neighbors.end()); - all_neighbors.insert(new_neighbors.begin(), new_neighbors.end()); - all_neighbors.insert(NULL); - for (Spin* sj : all_neighbors) { - if (sj != si) { - double p; - if (sj == NULL) { - Spin s0s_old = s0.inverse().act(*si); - Spin s0s_new = s0.inverse().act(si_new); - p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T); - A.ghost_bond_visited(*this, s0s_old, s0s_new, p); - } else { - p = 1.0 - exp(-(Z(*si, *sj) - Z(si_new, *sj)) / T); - A.plain_bond_visited(*this, si, sj, si_new, p); - } - if (rng.uniform(0.0, 1.0) < p) { - queue.push(sj); - } - } + std::set*> c = t->current(); + if (!std::any_of(c.begin(), c.end(), [&visited](Spin* s) { return visited.contains(s); })) { + visited.insert(c.begin(), c.end()); + + for (Spin* s : t->toConsider()) { + double ΔE = t->ΔE(s); + if (ΔE > 0 && 1.0 - exp(-ΔE / T) > rng.uniform(0, 1)) { + queue.push(t->createNew(s)); } - A.plain_site_transformed(*this, si, si_new); - dict.remove(si); - *si = si_new; - dict.insert(si); } + + t->apply(); } + delete t; } } @@ -148,12 +246,9 @@ public: measurement& A, unsigned N) { for (unsigned i = 0; i < N; i++) { Gen& g = rng.pick(gs); - R r = g(*this, rng); - unsigned ind = rng.uniform((unsigned)0, (unsigned)(s.size() - 1)); - - A.pre_cluster(*this, ind, r); + Transformation *t = g(*this, rng); - this->step(T, ind, r, A); + this->step(T, t, A); A.post_cluster(*this); } diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp index a9cd7d0..3505d6d 100644 --- a/spheres_infinite.cpp +++ b/spheres_infinite.cpp @@ -68,7 +68,7 @@ public: }; Gen, double> eGen(double ε) { - return [ε](const model& M, randutils::mt19937_rng& rng) -> Euclidean { + return [ε](model& M, Rng& rng) -> Transformation, double>* { Vector t; Matrix m; @@ -85,12 +85,12 @@ Gen, double> eGen(double ε) { } Euclidean g(t - m * t, m); - return g; + return new SpinFlip, double>(M, g, M.s[f_ind]); }; } Gen, double> mGen(double ε) { - return [ε](const model& M, randutils::mt19937_rng& rng) -> Euclidean { + return [ε](model& M, Rng& rng) -> Transformation, double>* { Matrix m; unsigned f_ind1 = rng.uniform((unsigned)0, (unsigned)M.s.size()); @@ -110,7 +110,7 @@ Gen, double> mGen(double ε) { m(1, 0) = -2 * cos(θ) * sin(θ); Euclidean g(t - m * t, m); - return g; + return new PairFlip, double>(M, g, M.s[f_ind1], M.s[f_ind2]); }; } @@ -123,9 +123,12 @@ int main(int argc, char* argv[]) { double H = 1.0; unsigned n = 25; + double k = 1e2; + double a = 0.1; + int opt; - while ((opt = getopt(argc, argv, "n:N:L:T:H:")) != -1) { + while ((opt = getopt(argc, argv, "n:N:L:T:H:a:k:")) != -1) { switch (opt) { case 'n': n = (unsigned)atof(optarg); @@ -142,14 +145,17 @@ int main(int argc, char* argv[]) { case 'H': H = atof(optarg); break; + case 'a': + a = atof(optarg); + break; + case 'k': + k = atof(optarg); + break; default: exit(1); } } - double k = 1000; - double a = 0.05; - std::function&, const Spin&)> Z = [L, a, k](const Spin& s1, const Spin& s2) -> double { Vector d = s1.x - s2.x; @@ -170,19 +176,19 @@ int main(int argc, char* argv[]) { return H * s.x.norm(); }; - auto g1 = eGen(0.25); - auto g2 = mGen(0.1); + auto g1 = eGen(0.5); + auto g2 = mGen(0.2); animation A(L, 750, argc, argv); model sphere(1, Z, B); - randutils::mt19937_rng rng; + Rng rng; sphere.s.reserve(n); unsigned nx = floor(sqrt(n)); for (unsigned i = 0; i < n; i++) { Vector pos = {(i / nx) * L / nx, (i % nx) * L / nx}; - sphere.s.push_back({pos, rng.uniform(0.45, 0.45)}); + sphere.s.push_back({pos, rng.uniform(0.25, 0.45)}); sphere.dict.insert(&sphere.s.back()); } -- cgit v1.2.3-70-g09d2