animation and mesuerements

2 files changed, 56 insertions, 16 deletions

diff --git a/space_wolff.hpp b/space_wolff.hpp
index 211e647..81ff776 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -245,22 +245,18 @@ public:
 
 template <class U, int D, class S> class Model;
 
-/*
 template <class U, int D, class S>
   class measurement {
     public:
       virtual void pre_cluster(const Model<U, D, S>&, unsigned, const Euclidean<U, D>&) {};
-      virtual void plain_bond_visited(const Model<U, D, S>&, const X_t&, double) {};
-      virtual void plain_site_transformed(const system<R_t, X_t, G_t>&, const typename G_t::vertex&
-v, const X_t&) {};
+      virtual void plain_bond_visited(const Model<U, D, 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, S>&, const Spin<U, D, S>*, const Spin<U, D, S>&) {};
 
-      virtual void ghost_bond_visited(const system<R_t, X_t, G_t>&, const typename G_t::vertex& v,
-const X_t&, const X_t&, double) {}; virtual void ghost_site_transformed(const system<R_t, X_t,
-G_t>&, const R_t&) {};
+      virtual void ghost_bond_visited(const Model<U, D, S>&, const Spin<U, D, S>&, const Spin<U, D, S>&, double) {};
+      virtual void ghost_site_transformed(const Model<U, D, S>&, const Euclidean<U, D>&) {};
 
-      virtual void post_cluster(unsigned, unsigned, const system<R_t, X_t, G_t>&) {};
+      virtual void post_cluster(const Model<U, D, S>&) {};
   };
-  */
 
 template <class U, int D, class S> class Model {
 public:
@@ -275,6 +271,7 @@ public:
   std::vector<Matrix<U, D>> mats;
   std::vector<Vector<U, D>> steps;
   double E;
+  measurement<U, D, S> &A;
 
   void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
     if (level > 0) {
@@ -290,8 +287,8 @@ 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, unsigned dDepth, unsigned nDepth)
-      : L(L), s0(L), dict(L, dDepth), Z(Z), B(B), dDepth(dDepth), nDepth(nDepth) {
+        std::function<double(const Spin<U, D, S>&)> B, unsigned dDepth, unsigned nDepth, measurement<U, D, S>& A)
+      : L(L), s0(L), dict(L, dDepth), Z(Z), B(B), dDepth(dDepth), nDepth(nDepth), A(A) {
     std::array<double, D> ini_sequence;
     ini_sequence.fill(1);
     std::list<std::array<double, D>> sequences;
@@ -383,11 +380,13 @@ public:
             Spin<U, D, S> s0s_old = s0.inverse().act(ss);
             Spin<U, D, S> 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 (dist(rng) < p) {
               queue.push(&ss);
             }
-            s0 = s0_new;
           }
+          A.ghost_site_transformed(*this, s0_new);
+          s0 = s0_new;
         } else {
           Spin<U, D, S> si_new = r.act(*si);
           std::set<Spin<U, D, S>*> all_neighbors;
@@ -404,14 +403,17 @@ public:
                 Spin<U, D, S> s0s_old = s0.inverse().act(*si);
                 Spin<U, D, S> 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_new, *sj) - Z(*si, *sj)) / T);
+                p = 1.0 - exp(-(Z(*si, *sj) - Z(si_new, *sj)) / T);
+                A.plain_bond_visited(*this, si, sj, si_new, p);
               }
               if (dist(rng) < p) {
                 queue.push(sj);
               }
             }
           }
+          A.plain_site_transformed(*this, si, si_new);
           dict.remove(si);
           *si = si_new;
           dict.insert(si);
@@ -452,7 +454,13 @@ public:
 
       Euclidean<U, D> g(L, t, m);
 
-      this->step(T, ind_dist(rng), g, rng);
+      unsigned ind = ind_dist(rng);
+
+      A.pre_cluster(*this, ind, g);
+
+      this->step(T, ind, g, rng);
+
+      A.post_cluster(*this);
 
       this->update_energy();
     }
diff --git a/spheres.cpp b/spheres.cpp
index 633a26a..25c5a67 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -1,5 +1,36 @@
 
 #include "space_wolff.hpp"
+#include <GL/glut.h>
+
+class animation : public measurement<double, 2, double> {
+  public:
+    animation(double L, unsigned w, int argc, char *argv[]) {
+      glutInit(&argc, argv);
+      glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
+      glutInitWindowSize(w, w);
+      glutCreateWindow("wolff");
+      glClearColor(0.0,0.0,0.0,0.0);
+      glMatrixMode(GL_PROJECTION);
+      glLoadIdentity();
+      gluOrtho2D(-1, L + 1, -1 , L + 1);
+    }
+
+    void post_cluster(const Model<double, 2, double>& m) override {
+      glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
+      glClear(GL_COLOR_BUFFER_BIT);
+      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));
+        }
+        glEnd();
+      }
+      glFlush();
+      getchar();
+    }
+};
 
 int main(int argc, char* argv[]) {
   const unsigned D = 2;
@@ -54,7 +85,8 @@ int main(int argc, char* argv[]) {
       return H * sin(2 * M_PI * 3 * s.x(0) / L);
     };
 
-  Model<double, D, double> sphere(L, Z, B, std::floor(log2(L)), 2);
+  animation A(L, 750, argc, argv);
+  Model<double, D, double> sphere(L, Z, B, std::floor(log2(L)), 2, A);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
@@ -65,7 +97,7 @@ int main(int argc, char* argv[]) {
 
   for (unsigned i = 0; i < n; i++) {
     Vector<double, D> pos = {dist(rng), dist(rng)};
-    sphere.s.push_back({pos, 0.5});
+    sphere.s.push_back({pos, dist(rng) / L});
     sphere.dict.insert(&sphere.s.back());
   }