From 698e4a89acc51f4da282d54cfac72938fe954a49 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 10 Jan 2020 16:36:58 -0500 Subject: progress toward measuring the variation in cluster size for the purpose of optimization --- space_wolff.hpp | 17 ++++++------- spheres.cpp | 74 ++++++++++++++++++++++++++++++++++++++++++++++----------- 2 files changed, 69 insertions(+), 22 deletions(-) diff --git a/space_wolff.hpp b/space_wolff.hpp index b542fbb..615fb93 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -364,17 +364,14 @@ public: unsigned nDepth; std::function&, const Spin&)> Z; std::function&)> B; - std::function&, randutils::mt19937_rng&)> g; randutils::mt19937_rng rng; - measurement& A; Model(U L, std::function&, const Spin&)> Z, - std::function&)> B, std::function&, randutils::mt19937_rng&)> g, unsigned dDepth, unsigned nDepth, - measurement& A) - : L(L), s0(L), dict(L, dDepth), Z(Z), B(B), g(g), dDepth(dDepth), nDepth(nDepth), rng(), A(A) {} + std::function&)> B, unsigned dDepth, unsigned nDepth) + : L(L), s0(L), dict(L, dDepth), Z(Z), B(B), dDepth(dDepth), nDepth(nDepth), rng() {} - void step(double T, unsigned ind, Euclidean r) { + void step(double T, unsigned ind, Euclidean r, measurement& A) { unsigned cluster_size = 0; std::queue*> queue; @@ -438,14 +435,18 @@ public: } } - void wolff(double T, unsigned N) { + void resize(double T, double P) { + } + + void wolff(double T, std::function&, randutils::mt19937_rng&)> g, + measurement& A, unsigned N) { for (unsigned i = 0; i < N; i++) { R r = g(*this, rng); unsigned ind = rng.uniform((unsigned)0, (unsigned)(s.size() - 1)); A.pre_cluster(*this, ind, r); - this->step(T, ind, r); + this->step(T, ind, r, A); A.post_cluster(*this); } diff --git a/spheres.cpp b/spheres.cpp index 4ed329e..f3dde42 100644 --- a/spheres.cpp +++ b/spheres.cpp @@ -6,8 +6,17 @@ const unsigned D = 2; typedef Model, double> model; class animation : public measurement, double> { + private: + uint64_t t1; + uint64_t t2; + unsigned n; + unsigned tmp; public: animation(double L, unsigned w, int argc, char *argv[]) { + t1 = 0; + t2 = 0; + n = 0; + glutInit(&argc, argv); glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); glutInitWindowSize(w, w); @@ -18,6 +27,14 @@ class animation : public measurement, double> { gluOrtho2D(-1, L + 1, -1 , L + 1); } + void pre_cluster(const model&, unsigned, const Euclidean&) override { + tmp = 0; + } + + void plain_site_transformed(const model&, const Spin*, const Spin&) override { + tmp++; + } + void post_cluster(const model& m) override { glClearColor(1.0f, 1.0f, 1.0f, 1.0f ); glClear(GL_COLOR_BUFFER_BIT); @@ -31,12 +48,25 @@ class animation : public measurement, double> { glEnd(); } glFlush(); - getchar(); + + t1 += tmp; + t2 += tmp * tmp; + n++; + } + + void clear() { + t1 = 0; + t2 = 0; + n = 0; + } + + double var() { + return (t2 - t1 * t1 / (double)n) / (double)n; } }; -std::function(const model&, randutils::mt19937_rng&)> eGen(const std::vector>& mats, const std::vector>& vecs) { - return [&mats, &vecs] (const model& M, randutils::mt19937_rng& rng) -> Euclidean { +std::function(const model&, randutils::mt19937_rng&)> eGen(const std::vector>& mats, const std::vector>& vecs, double ε, double L) { + return [&mats, &vecs, L, ε] (const model& M, randutils::mt19937_rng& rng) -> Euclidean { Vector t; Matrix m; unsigned flip = rng.uniform((unsigned)0, (unsigned)(mats.size() + vecs.size() - 1)); @@ -44,7 +74,7 @@ std::function(const model&, randutils::mt19937_rng&)> eGen( unsigned f_ind = rng.uniform((unsigned)0, (unsigned)M.s.size()); t = M.s[f_ind].x; for (unsigned j = 0; j < D; j++) { - t(j) += rng.variate(0.0, 0.1); + t(j) = fmod(10 * L + t(j) + rng.variate(0.0, ε), L); } m = mats[flip]; } else { @@ -100,16 +130,20 @@ int main(int argc, char* argv[]) { } } + double k = 100.0; + double a = 0.2; + std::function&, const Spin&)> Z = - [L] (const Spin& s1, const Spin& s2) -> double { + [L, a, k] (const Spin& s1, const Spin& s2) -> double { Vector d = diff(L, s1.x, s2.x); - double rad_sum = pow(s1.s + s2.s, 2); - - bool overlap = d.transpose() * d < rad_sum; + double σ = s1.s + s2.s; + double δ = σ - sqrt(d.transpose() * d); - if (overlap) { - return -1e8; + if (δ > - a * σ) { + return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2)); + } else if (δ > - 2 * a * σ) { + return 0.5 * k * pow(δ + 2 * a * σ, 2); } else { return 0; } @@ -122,9 +156,9 @@ int main(int argc, char* argv[]) { std::vector> mats = torus_mats(); std::vector> vecs = torus_vecs(L); - auto g = eGen(mats, vecs); + auto g = eGen(mats, vecs, L, L); animation A(L, 750, argc, argv); - model sphere(L, Z, B, g, std::floor(log2(L)), 2, A); + model sphere(L, Z, B, std::floor(log2(L)), 2); randutils::mt19937_rng rng; @@ -132,13 +166,25 @@ int main(int argc, char* argv[]) { unsigned nx = floor(sqrt(n)); for (unsigned i = 0; i < n; i++) { - Vector pos = {(i / nx) * L / nx + rng.uniform(0.0, L / (4 * nx)), (i % nx) * L / nx + rng.uniform(0.0, L / (4 * nx))}; + Vector pos = {(i / nx) * L / nx, (i % nx) * L / nx}; sphere.s.push_back({pos, 0.5}); sphere.dict.insert(&sphere.s.back()); } + sphere.wolff(T, g, A, N); + std::ofstream outfile; + outfile.open("test.dat"); + + for (signed i = -10; i <= 10; i++) { + A.clear(); + double ε = pow(2, -4 + i / 2.0); + auto gn = eGen(mats, vecs, ε, L); + sphere.wolff(T, gn, A, N); + outfile << ε << " " << A.var() / sphere.s.size() << std::endl; + std::cout << ε << " " << A.var() / sphere.s.size() << std::endl; + } - sphere.wolff(T, N); + outfile.close(); std::ofstream snapfile; snapfile.open("sphere_snap.dat"); -- cgit v1.2.3-54-g00ecf