generalized code somewhat to accomodate simulation of infinite space, sample central field sim

3 files changed, 340 insertions, 24 deletions

diff --git a/space_wolff.hpp b/space_wolff.hpp
index 615fb93..7345726 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -9,6 +9,25 @@
 #include <random>
 #include <set>
 #include <vector>
+#include <unordered_map>
+#include <array>
+#include <unordered_map>
+
+namespace std {
+template <typename T, size_t N> struct hash<array<T, N>> {
+  typedef array<T, N> argument_type;
+  typedef size_t result_type;
+
+  result_type operator()(const argument_type& a) const {
+    hash<T> hasher;
+    result_type h = 0;
+    for (result_type i = 0; i < N; ++i) {
+      h = h * 31 + hasher(a[i]);
+    }
+    return h;
+  }
+};
+} // namespace std
 
 #include "randutils/randutils.hpp"
 
@@ -56,6 +75,54 @@ public:
 };
 
 template <class T, int D> class Euclidean {
+public:
+  Vector<T, D> t;
+  Matrix<T, D> r;
+
+  Euclidean(T L) {
+    for (unsigned i = 0; i < D; i++) {
+      t(i) = 0;
+      r(i, i) = 1;
+      for (unsigned j = 1; j < D; j++) {
+        r(i, (i + j) % D) = 0;
+      }
+    }
+  }
+
+  Euclidean(Vector<T, D> t0, Matrix<T, D> r0) {
+    t = t0;
+    r = r0;
+  }
+
+  template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
+    Spin<T, D, S> s_new;
+
+    s_new.x = t + r * s.x;
+    s_new.s = s.s;
+
+    return s_new;
+  }
+
+  Euclidean act(const Euclidean& x) const {
+    Vector<T, D> tnew = r * x.t + t;
+    Matrix<T, D> rnew = r * x.r;
+
+    Euclidean pnew(tnew, rnew);
+
+    return pnew;
+  }
+
+  Euclidean inverse() const {
+    Vector<T, D> tnew = -r.transpose() * t;
+    Matrix<T, D> rnew = r.transpose();
+
+    Euclidean pnew(tnew, rnew);
+
+    return pnew;
+  }
+};
+
+template <class T, int D> class TorusGroup {
 private:
   T L;
 
@@ -63,9 +130,9 @@ public:
   Vector<T, D> t;
   Matrix<T, D> r;
 
-  /** brief Euclidean - default constructor, constructs the identity
+  /** brief TorusGroup - default constructor, constructs the identity
    */
-  Euclidean(T L) : L(L) {
+  TorusGroup(T L) : L(L) {
     for (unsigned i = 0; i < D; i++) {
       t(i) = 0;
       r(i, i) = 1;
@@ -75,7 +142,7 @@ public:
     }
   }
 
-  Euclidean(T L, Vector<T, D> t0, Matrix<T, D> r0) : L(L) {
+  TorusGroup(T L, Vector<T, D> t0, Matrix<T, D> r0) : L(L) {
     t = t0;
     r = r0;
   }
@@ -93,7 +160,7 @@ public:
     return s_new;
   }
 
-  Euclidean act(const Euclidean& x) const {
+  TorusGroup act(const TorusGroup& x) const {
     Vector<T, D> tnew = r * x.t + t;
     Matrix<T, D> rnew = r * x.r;
 
@@ -101,16 +168,16 @@ public:
       tnew(i) = fmod(L + tnew(i), L);
     }
 
-    Euclidean pnew(this->L, tnew, rnew);
+    TorusGroup pnew(this->L, tnew, rnew);
 
     return pnew;
   }
 
-  Euclidean inverse() const {
+  TorusGroup inverse() const {
     Vector<T, D> tnew = -r.transpose() * t;
     Matrix<T, D> rnew = r.transpose();
 
-    Euclidean pnew(this->L, tnew, rnew);
+    TorusGroup pnew(this->L, tnew, rnew);
 
     return pnew;
   }
@@ -118,6 +185,67 @@ public:
 
 template <class T, int D, class S> class Octree {
 private:
+  unsigned N;
+  T L;
+  std::unordered_map<std::array<signed, D>, std::set<Spin<T, D, S>*>> data;
+
+public:
+  Octree(T L, unsigned N) {
+    L = L;
+    N = N;
+  }
+
+  std::array<signed, D> ind(Vector<T, D> x) const {
+    std::array<signed, D> ind;
+
+    for (unsigned i = 0; i < D; i++) {
+      ind[i] = (signed)std::floor(N * x(i) / L);
+    }
+
+    return ind;
+  }
+
+  void insert(Spin<T, D, S>* s) { data[ind(s->x)].insert(s); };
+
+  void remove(Spin<T, D, S>* s) {
+    data[ind(s->x)].erase(s);
+    if (data[ind(s->x)].empty()) {
+      data.erase(ind(s->x));
+    }
+  };
+
+  std::set<Spin<T, D, S>*> neighbors(const Vector<T, D>& x) const {
+    std::array<signed, D> i0 = ind(x);
+    std::set<Spin<T, D, S>*> ns;
+
+    nearest_neighbors_of(i0, D + 1, ns);
+
+    return ns;
+  };
+
+  void nearest_neighbors_of(std::array<signed, D> i0, unsigned depth, std::set<Spin<T, D, S>*>& ns) const {
+    if (depth == 0) {
+      auto it = data.find(i0);
+      if (it != data.end()) {
+        ns.insert(it->second.begin(), it->second.end());
+      }
+    } else {
+      for (signed j : {-1, 0, 1}) {
+        std::array<signed, D> i1 = i0;
+        if (N < 2) {
+          i1[depth - 1] += j;
+        } else {
+          i1[depth - 1] = (N + i1[depth - 1] + j) % N;
+        }
+        nearest_neighbors_of(i1, depth - 1, ns);
+      }
+    }
+  };
+};
+
+/*
+template <class T, int D, class S> class Octree {
+private:
   T L;
   unsigned N;
   std::vector<std::set<Spin<T, D, S>*>> data;
@@ -168,6 +296,7 @@ public:
     }
   }
 };
+*/
 
 class quantity {
 private:
@@ -360,18 +489,16 @@ public:
   R s0;
   std::vector<Spin<U, D, S>> s;
   Octree<U, D, S> dict;
-  unsigned dDepth;
-  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;
   randutils::mt19937_rng rng;
 
 
   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), rng() {}
+        std::function<double(const Spin<U, D, S>&)> B)
+      : L(L), s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {}
 
-  void step(double T, unsigned ind, Euclidean<U, D> r, measurement<U, D, R, S>& A) {
+  void step(double T, unsigned ind, R& r, measurement<U, D, R, S>& A) {
     unsigned cluster_size = 0;
 
     std::queue<Spin<U, D, S>*> queue;
@@ -387,7 +514,7 @@ public:
         visited.insert(si);
 
         if (si == NULL) {
-          Euclidean<U, D> s0_new = r.act(s0);
+          R s0_new = r.act(s0);
           for (Spin<U, D, S>& ss : s) {
             Spin<U, D, S> s0s_old = s0.inverse().act(ss);
             Spin<U, D, S> s0s_new = s0_new.inverse().act(ss);
@@ -402,8 +529,8 @@ public:
         } else {
           Spin<U, D, S> si_new = r.act(*si);
           std::set<Spin<U, D, S>*> all_neighbors;
-          std::set<Spin<U, D, S>*> current_neighbors = dict.neighbors(si->x, nDepth);
-          std::set<Spin<U, D, S>*> new_neighbors = dict.neighbors(si_new.x, nDepth);
+          std::set<Spin<U, D, S>*> current_neighbors = dict.neighbors(si->x);
+          std::set<Spin<U, D, S>*> 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());
diff --git a/spheres.cpp b/spheres.cpp
index f3dde42..97c1e18 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -3,9 +3,9 @@
 #include <GL/glut.h>
 
 const unsigned D = 2;
-typedef Model<double, D, Euclidean<double, D>, double> model;
+typedef Model<double, D, TorusGroup<double, D>, double> model;
 
-class animation : public measurement<double, 2, Euclidean<double, D>, double> {
+class animation : public measurement<double, 2, TorusGroup<double, D>, double> {
   private:
     uint64_t t1;
     uint64_t t2;
@@ -27,7 +27,7 @@ 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 {
+    void pre_cluster(const model&, unsigned, const TorusGroup<double, D>&) override {
       tmp = 0;
     }
 
@@ -65,8 +65,8 @@ class animation : public measurement<double, 2, Euclidean<double, D>, double> {
     }
 };
 
-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> {
+std::function<TorusGroup<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) -> TorusGroup<double, 2> {
     Vector<double, D> t;
     Matrix<double, D> m;
     unsigned flip = rng.uniform((unsigned)0, (unsigned)(mats.size() + vecs.size() - 1));
@@ -91,7 +91,7 @@ std::function<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(
       t = vecs[flip - mats.size()];
     }
 
-    Euclidean<double, D> g(M.L, t, m);
+    TorusGroup<double, D> g(M.L, t, m);
     return g;
   };
 
@@ -130,8 +130,8 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  double k = 100.0;
-  double a = 0.2;
+  double k = 1e2;
+  double a = 0.05;
 
   std::function<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z =
     [L, a, k] (const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
@@ -158,7 +158,7 @@ int main(int argc, char* argv[]) {
   std::vector<Vector<double, D>> vecs = torus_vecs<double, D>(L);
   auto g = eGen(mats, vecs, L, L);
   animation A(L, 750, argc, argv);
-  model sphere(L, Z, B, std::floor(log2(L)), 2);
+  model sphere(L, Z, B);
 
   randutils::mt19937_rng rng;
 
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
new file mode 100644
index 0000000..f2f998a
--- /dev/null
+++ b/spheres_infinite.cpp
@@ -0,0 +1,189 @@
+
+#include "space_wolff.hpp"
+#include <GL/glut.h>
+
+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);
+      glutCreateWindow("wolff");
+      glClearColor(0.0,0.0,0.0,0.0);
+      glMatrixMode(GL_PROJECTION);
+      glLoadIdentity();
+      gluOrtho2D(-L, L, -L , L);
+    }
+
+    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);
+      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();
+
+      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(double ε) {
+  return [ε] (const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
+    Vector<double, D> t;
+    Matrix<double, D> m;
+
+    double θ = rng.uniform((double)0.0, 2 * M_PI);
+    m(0, 0) = cos(θ);
+    m(1, 1) = cos(θ);
+    m(0, 1) = sin(θ);
+    m(1, 0) = -sin(θ);
+
+    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, ε);
+    }
+
+    Euclidean<double, D> g(t, m);
+    return g;
+  };
+
+}
+
+int main(int argc, char* argv[]) {
+  const unsigned D = 2;
+
+  double L = 32;
+  unsigned N = 1000;
+  double T = 2.0 / log(1.0 + sqrt(2.0));
+  double H = 1.0;
+  unsigned n = 25;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "n:N:L:T:H:")) != -1) {
+    switch (opt) {
+      case 'n': 
+        n = (unsigned)atof(optarg);
+        break;
+      case 'N': 
+        N = (unsigned)atof(optarg);
+        break;
+      case 'L':
+        L = atof(optarg);
+        break;
+      case 'T':
+        T = atof(optarg);
+        break;
+      case 'H':
+        H = atof(optarg);
+        break;
+      default:
+        exit(1);
+    }
+  }
+
+  double k = 1e2;
+  double a = 0.05;
+
+  std::function<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z =
+    [L, a, k] (const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
+      Vector<double, D> d = s1.x - s2.x;
+
+      double σ = s1.s + s2.s;
+      double δ = σ - sqrt(d.transpose() * d);
+
+      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;
+      }
+    };
+
+  std::function<double(Spin<double, D, double>)> B =
+    [L, H] (Spin<double, D, double> s) -> double {
+      return H * s.x.norm();
+    };
+
+  auto g = eGen(1);
+  animation A(L, 750, argc, argv);
+  model sphere(1, Z, B);
+
+  randutils::mt19937_rng rng;
+
+  sphere.s.reserve(n);
+
+  unsigned nx = floor(sqrt(n));
+  for (unsigned i = 0; i < n; i++) {
+    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(ε);
+    sphere.wolff(T, gn, A, N);
+    outfile << ε << " " << A.var() / sphere.s.size() << std::endl;
+    std::cout << ε << " " << A.var() / sphere.s.size() << std::endl;
+  }
+
+  outfile.close();
+
+  std::ofstream snapfile;
+  snapfile.open("sphere_snap.dat");
+
+  for (Spin<double, D, double> s : sphere.s) {
+    Spin<double, D, double> rs = sphere.s0.inverse().act(s);
+    snapfile << rs.x.transpose() << "\n";
+  }
+
+  snapfile.close();
+
+  return 0;
+}