progress toward measuring the variation in cluster size for the purpose of optimization

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");