From 698e4a89acc51f4da282d54cfac72938fe954a49 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
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<double(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z;
std::function<double(const Spin<U, D, S>&)> B;
- std::function<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)> g;
randutils::mt19937_rng rng;
- measurement<U, D, R, S>& A;
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, std::function<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)> g, unsigned dDepth, unsigned nDepth,
- measurement<U, D, R, S>& A)
- : L(L), s0(L), dict(L, dDepth), Z(Z), B(B), g(g), dDepth(dDepth), nDepth(nDepth), rng(), A(A) {}
+ std::function<double(const Spin<U, D, S>&)> 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<U, D> r) {
+ void step(double T, unsigned ind, Euclidean<U, D> r, measurement<U, D, R, S>& A) {
unsigned cluster_size = 0;
std::queue<Spin<U, D, S>*> queue;
@@ -438,14 +435,18 @@ public:
}
}
- void wolff(double T, unsigned N) {
+ void resize(double T, double P) {
+ }
+
+ void wolff(double T, std::function<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)> g,
+ measurement<U, D, R, S>& 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, D, Euclidean<double, D>, double> model;
class animation : public measurement<double, 2, Euclidean<double, D>, 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, 2, Euclidean<double, D>, double> {
gluOrtho2D(-1, L + 1, -1 , L + 1);
}
+ 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 {
+ 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, 2, Euclidean<double, D>, 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<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(const std::vector<Matrix<double, 2>>& mats, const std::vector<Vector<double, 2>>& vecs) {
- return [&mats, &vecs] (const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
+std::function<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(const std::vector<Matrix<double, 2>>& mats, const std::vector<Vector<double, 2>>& vecs, double ε, double L) {
+ return [&mats, &vecs, L, ε] (const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
Vector<double, D> t;
Matrix<double, D> m;
unsigned flip = rng.uniform((unsigned)0, (unsigned)(mats.size() + vecs.size() - 1));
@@ -44,7 +74,7 @@ std::function<Euclidean<double, D>(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<double, std::normal_distribution>(0.0, 0.1);
+ t(j) = fmod(10 * L + t(j) + rng.variate<double, std::normal_distribution>(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<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z =
- [L] (const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
+ [L, a, k] (const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
Vector<double, D> 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<Matrix<double, D>> mats = torus_mats<double, D>();
std::vector<Vector<double, D>> vecs = torus_vecs<double, D>(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<double, D> pos = {(i / nx) * L / nx + rng.uniform(0.0, L / (4 * nx)), (i % nx) * L / nx + rng.uniform(0.0, L / (4 * nx))};
+ Vector<double, D> 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-70-g09d2