Merge branch 'master' of git:research/wolff/code/space_wolff

2 files changed, 142 insertions, 161 deletions

diff --git a/space_wolff.hpp b/space_wolff.hpp
index f95c68a..211e647 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -29,13 +29,14 @@ const std::array<std::array<unsigned, 16>, 16> smiley = {
      {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}},
      {{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}}}};
 
-template <class U, unsigned D> using Vector = Eigen::Matrix<U, D, 1>;
+template <class T, int D> using Vector = Eigen::Matrix<T, D, 1>;
 
-template <class U, unsigned D> using Matrix = Eigen::Matrix<U, D, D>;
+template <class T, int D> using Matrix = Eigen::Matrix<T, D, D>;
 
-template <class U, unsigned D>
-Vector<U, D> diff(U L, Vector<U, D> v1, Vector<U, D> v2) {
-  Vector<U, D> v;
+/** brief diff - periodic subtraction of vectors
+ */
+template <class T, int D> Vector<T, D> diff(T L, Vector<T, D> v1, Vector<T, D> v2) {
+  Vector<T, D> v;
 
   for (unsigned i = 0; i < D; i++) {
     v(i) = std::abs(v1(i) - v2(i));
@@ -47,20 +48,23 @@ Vector<U, D> diff(U L, Vector<U, D> v1, Vector<U, D> v2) {
   return v;
 }
 
-template <class U, unsigned D, class state> class Spin {
+template <class T, int D, class S> class Spin {
 public:
-  Vector<U, D> x;
-  state s;
+  Vector<T, D> x;
+  S s;
 };
 
-template <class U, unsigned D> class Euclidean {
+template <class T, int D> class Euclidean {
 private:
-  U L;
+  T L;
 
 public:
-  Vector<U, D> t;
-  Matrix<U, D> r;
-  Euclidean(U L) : L(L) {
+  Vector<T, D> t;
+  Matrix<T, D> r;
+
+  /** brief Euclidean - default constructor, constructs the identity
+   */
+  Euclidean(T L) : L(L) {
     for (unsigned i = 0; i < D; i++) {
       t(i) = 0;
       r(i, i) = 1;
@@ -70,14 +74,13 @@ public:
     }
   }
 
-  Euclidean(U L, Vector<U, D> t0, Matrix<U, D> r0) : L(L) {
+  Euclidean(T L, Vector<T, D> t0, Matrix<T, D> r0) : L(L) {
     t = t0;
     r = r0;
   }
 
-  template <class state>
-  Spin<U, D, state> act(const Spin<U, D, state> &s) const {
-    Spin<U, D, state> s_new;
+  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;
@@ -89,9 +92,9 @@ public:
     return s_new;
   }
 
-  Euclidean act(const Euclidean &x) const {
-    Vector<U, D> tnew = r * x.t + t;
-    Matrix<U, D> rnew = r * x.r;
+  Euclidean act(const Euclidean& x) const {
+    Vector<T, D> tnew = r * x.t + t;
+    Matrix<T, D> rnew = r * x.r;
 
     for (unsigned i = 0; i < D; i++) {
       tnew(i) = fmod(L + tnew(i), L);
@@ -103,8 +106,8 @@ public:
   }
 
   Euclidean inverse() const {
-    Vector<U, D> tnew = -r.transpose() * t;
-    Matrix<U, D> rnew = r.transpose();
+    Vector<T, D> tnew = -r.transpose() * t;
+    Matrix<T, D> rnew = r.transpose();
 
     Euclidean pnew(this->L, tnew, rnew);
 
@@ -112,65 +115,60 @@ public:
   }
 };
 
-template <class T, unsigned D> class Dictionary {
+template <class T, int D, class S> class Octree {
 private:
-  unsigned N;
   T L;
-  std::vector<std::set<unsigned>> d;
+  unsigned N;
+  std::vector<std::set<Spin<T, D, S>*>> data;
 
 public:
-  Dictionary(unsigned Ni, double Li) : N(Ni), L(Li), d(pow(Ni, D)){};
+  Octree(T L, unsigned depth) : L(L), N(pow(2, depth)), data(pow(N, D)){};
 
-  unsigned dictionary_index(Vector<T, D> x) const {
+  unsigned ind(Vector<T, D> x) const {
     unsigned pos_ind = 0;
 
     for (unsigned i = 0; i < D; i++) {
-      pos_ind += pow(N, i) * (unsigned)std::floor(x(i) * N / L);
-    };
+      pos_ind += pow(N, i) * (unsigned)std::floor(N * x(i) / L);
+    }
+
+    assert(pos_ind < data.size());
 
     return pos_ind;
   }
 
-  void record(Vector<T, D> x, unsigned ind) {
-    d[this->dictionary_index(x)].insert(ind);
-  };
+  void insert(Spin<T, D, S>* s) { data[ind(s->x)].insert(s); };
 
-  void remove(Vector<T, D> x, unsigned ind) {
-    d[this->dictionary_index(x)].erase(ind);
-  };
+  void remove(Spin<T, D, S>* s) { data[ind(s->x)].erase(s); };
 
-  std::set<unsigned> neighbors(Vector<T, D> x, unsigned depth) const {
-    return nearest_neighbors_of(this->dictionary_index(x), depth, {});
+  std::set<Spin<T, D, S>*> neighbors(const Vector<T, D>& x, unsigned depth) const {
+    std::set<Spin<T, D, S>*> ns;
+    std::set<unsigned> seen;
+    nearest_neighbors_of(ind(x), depth, ns, seen);
+    return ns;
   };
 
-  std::set<unsigned> nearest_neighbors_of(unsigned ind, unsigned depth,
-                                          std::list<unsigned> ignores) const {
-    std::set<unsigned> ns = d[ind];
+  void nearest_neighbors_of(unsigned ind, unsigned depth, std::set<Spin<T, D, S>*>& ns,
+                            std::set<unsigned>& seen) const {
+    ns.insert(data[ind].begin(), data[ind].end());
+    seen.insert(ind);
 
     if (depth > 0) {
       for (unsigned i = 0; i < D; i++) {
-        if (std::none_of(ignores.begin(), ignores.end(), [i](unsigned k) { return i == k; })) {
-          std::list<unsigned> ignores_new = ignores;
-          ignores_new.push_back(i);
-          for (signed j : {-1, 1}) {
-            unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) +
-                          fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1));
+        for (signed j : {-1, 1}) {
+          unsigned nind = pow(N, i + 1) * (ind / ((unsigned)pow(N, i + 1))) +
+                          fmod(pow(N, i + 1) + ind + j * pow(N, i), pow(N, i + 1));
 
-            std::set<unsigned> nns = nearest_neighbors_of(ni, depth - 1, ignores_new);
-
-            for (unsigned guy : nns) {
-              ns.insert(guy);
-            }
+          if (!seen.contains(nind)) {
+            seen.insert(nind);
+            nearest_neighbors_of(nind, depth - 1, ns, seen);
           }
         }
       }
     }
-
-    return ns;
-  };
+  }
 };
 
-class Quantity {
+class quantity {
 private:
   double total;
   double total2;
@@ -181,7 +179,7 @@ public:
   unsigned n;
   std::list<double> hist;
 
-  Quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) {
+  quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) {
     n = 0;
     total = 0;
     total2 = 0;
@@ -245,22 +243,38 @@ public:
   unsigned num_added() const { return n - wait; }
 };
 
-template <class U, unsigned D, class state> class Model {
+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 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 post_cluster(unsigned, unsigned, const system<R_t, X_t, G_t>&) {};
+  };
+  */
+
+template <class U, int D, class S> class Model {
 public:
   U L;
   Euclidean<U, D> s0;
-  std::vector<Spin<U, D, state>> s;
-  Dictionary<U, D> dict;
-  std::function<std::set<unsigned>(Model<U, D, state> &, unsigned,
-                                   Spin<U, D, state>)>
-      neighbors;
-  std::function<double(Spin<U, D, state>, Spin<U, D, state>)> Z;
-  std::function<double(Spin<U, D, state>)> B;
+  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;
   std::vector<Matrix<U, D>> mats;
   std::vector<Vector<U, D>> steps;
-  long double E;
-  Quantity Eq;
-  Quantity Cq;
+  double E;
 
   void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
     if (level > 0) {
@@ -275,14 +289,9 @@ public:
     }
   }
 
-  Model(U L, unsigned N,
-        std::function<double(Spin<U, D, state>, Spin<U, D, state>)> Z,
-        std::function<double(Spin<U, D, state>)> B,
-        std::function<std::set<unsigned>(Model<U, D, state> &, unsigned,
-                                         Spin<U, D, state>)>
-            ns)
-      : L(L), s0(L), dict(N, L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000),
-        Cq(1000, 1000) {
+  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::array<double, D> ini_sequence;
     ini_sequence.fill(1);
     std::list<std::array<double, D>> sequences;
@@ -345,82 +354,71 @@ public:
     E = 0;
     for (unsigned i = 0; i < s.size(); i++) {
       for (unsigned j = 0; j < i; j++) {
-         E -= Z(s[i], s[j]);
+        E -= Z(s[i], s[j]);
       }
 
       E -= B(s0.inverse().act(s[i]));
     }
   }
 
-  void step(double T, unsigned ind, Euclidean<U, D> r, std::mt19937 &rng) {
+  void step(double T, unsigned ind, Euclidean<U, D> r, std::mt19937& rng) {
     unsigned cluster_size = 0;
     std::uniform_real_distribution<double> dist(0.0, 1.0);
 
-    std::queue<unsigned> queue;
-    queue.push(ind);
+    std::queue<Spin<U, D, S>*> queue;
+    queue.push(&(s[ind]));
 
-    std::vector<bool> visited(s.size() + 1, false);
+    std::set<Spin<U, D, S>*> visited;
 
     while (!queue.empty()) {
-      unsigned i = queue.front();
+      Spin<U, D, S>* si = queue.front();
       queue.pop();
 
-      if (!visited[i]) {
-        visited[i] = true;
-
-        bool we_are_ghost = i == s.size();
-
-        Spin<U, D, state> si_new;
-        Euclidean<U, D> s0_new(L);
+      if (!visited.contains(si)) {
+        visited.insert(si);
 
-        if (we_are_ghost) {
-          s0_new = r.act(s0);
+        if (si == NULL) {
+          Euclidean<U, D> 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);
+            double p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T);
+            if (dist(rng) < p) {
+              queue.push(&ss);
+            }
+            s0 = s0_new;
+          }
         } else {
-          si_new = r.act(s[i]);
-        }
-
-        for (unsigned j : neighbors(*this, i, si_new)) {
-          if (j != i) {
-            double p;
-            bool neighbor_is_ghost = j == s.size();
-
-            if (we_are_ghost || neighbor_is_ghost) {
-              Spin<U, D, state> s0s_old, s0s_new;
-              unsigned non_ghost;
-
-              if (neighbor_is_ghost) {
-                non_ghost = i;
-                s0s_old = s0.inverse().act(s[i]);
-                s0s_new = s0.inverse().act(si_new);
+          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);
+
+          all_neighbors.insert(current_neighbors.begin(), current_neighbors.end());
+          all_neighbors.insert(new_neighbors.begin(), new_neighbors.end());
+          all_neighbors.insert(NULL);
+          for (Spin<U, D, S>* sj : all_neighbors) {
+            if (sj != si) {
+              double p;
+              if (sj == NULL) {
+                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);
               } else {
-                non_ghost = j;
-                s0s_old = s0.inverse().act(s[j]);
-                s0s_new = s0_new.inverse().act(s[j]);
+                p = 1.0 - exp(-(Z(si_new, *sj) - Z(*si, *sj)) / T);
+              }
+              if (dist(rng) < p) {
+                queue.push(sj);
               }
-
-              p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T);
-            } else {
-              p = 1.0 - exp(-(Z(si_new, s[j]) - Z(s[i], s[j])) / T);
-            }
-
-            if (dist(rng) < p) {
-              queue.push(j);
             }
           }
-        }
-
-        if (we_are_ghost) {
-          s0 = s0_new;
-        } else {
-          dict.remove(s[i].x, i);
-          s[i] = si_new;
-          dict.record(s[i].x, i);
+          dict.remove(si);
+          *si = si_new;
+          dict.insert(si);
           cluster_size++;
         }
       }
     }
-
-    Cq.add(cluster_size);
   }
 
   void wolff(double T, unsigned N, std::mt19937& rng) {
@@ -457,7 +455,6 @@ public:
       this->step(T, ind_dist(rng), g, rng);
 
       this->update_energy();
-      Eq.add(E);
     }
   }
 };
diff --git a/spheres.cpp b/spheres.cpp
index b0f4fb4..633a26a 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -34,65 +34,49 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  std::function<double(spin<double, D, double>, spin<double, D, double>, spin<double, D, double>)> Z =
-    [L] (spin<double, D, double> s1, spin<double, D, double> s2, spin<double, D, double> s1_new) -> double {
-      vector<double, D> diff_old = diff(L, s1.x, s2.x);
-      vector<double, D> diff_new = diff(L, s1_new.x, s2.x);
+  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 {
+      Vector<double, D> d = diff(L, s1.x, s2.x);
 
       double rad_sum = pow(s1.s + s2.s, 2);
 
-      bool old_overlap = diff_old.transpose() * diff_old < rad_sum;
-      bool new_overlap = diff_new.transpose() * diff_new < rad_sum;
+      bool overlap = d.transpose() * d < rad_sum;
 
-      if (new_overlap) {
-        return 1.0;
+      if (overlap) {
+        return -1e8;
       } else {
-        return 0.0;
+        return 0;
       }
     };
 
-  std::function<double(spin<double, D, double>)> B =
-    [L, H] (spin<double, D, double> s) -> double {
+  std::function<double(Spin<double, D, double>)> B =
+    [L, H] (Spin<double, D, double> s) -> double {
       return H * sin(2 * M_PI * 3 * s.x(0) / L);
     };
 
-  std::function<std::set<unsigned>(model<double, D, double>&, unsigned, spin<double, D, double>)> neighbors =
-    [] (model<double, D, double>& m, unsigned i0, spin<double, D, double> s1) -> std::set<unsigned> {
-      std::set<unsigned> nn;
-      if (i0 < m.s.size()) {
-        std::set<unsigned> n_old = m.dict.neighbors(m.s[i0].x, 1);
-        std::set<unsigned> n_new = m.dict.neighbors(s1.x, 1);;
-        nn.insert(n_old.begin(), n_old.end());
-        nn.insert(n_new.begin(), n_new.end());
-        nn.insert(m.s.size());
-      } else {
-        for (unsigned i = 0; i < m.s.size(); i++) {
-          nn.insert(i);
-        }
-      }
-      return nn;
-    };
-
-  model<double, D, double> sphere(L, Z, B, neighbors);
+  Model<double, D, double> sphere(L, Z, B, std::floor(log2(L)), 2);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
 
   std::uniform_real_distribution<double> dist(0.0, L);
 
+  sphere.s.reserve(n);
+
   for (unsigned i = 0; i < n; i++) {
-    vector<double, D> pos = {dist(rng), dist(rng)};
+    Vector<double, D> pos = {dist(rng), dist(rng)};
     sphere.s.push_back({pos, 0.5});
-    sphere.dict.record<double>(pos, i);
+    sphere.dict.insert(&sphere.s.back());
   }
 
+
   sphere.wolff(T, N, rng);
 
   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);
+  for (Spin<double, D, double> s : sphere.s) {
+    Spin<double, D, double> rs = sphere.s0.inverse().act(s);
     snapfile << rs.x.transpose() << "\n";
   }