From b849c84315a262ac47010cdee306875b650aa918 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 8 Jan 2020 17:22:42 -0500
Subject: many updates

---
 space_wolff.hpp | 20 +++++++++-----------
 spheres.cpp     | 23 +++++++++--------------
 2 files changed, 18 insertions(+), 25 deletions(-)

diff --git a/space_wolff.hpp b/space_wolff.hpp
index 1335f07..b542fbb 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -271,7 +271,7 @@ void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level)
       new_sequence[new_level] = -1;
       sequences.push_front(new_sequence);
     }
-    one_sequences(sequences, new_level);
+    one_sequences<D>(sequences, new_level);
   }
 }
 
@@ -282,6 +282,7 @@ std::vector<Vector<U, D>> torus_vecs(U L) {
   ini_sequence.fill(1);
   std::list<std::array<double, D>> sequences;
   sequences.push_back(ini_sequence);
+  one_sequences<D>(sequences, D);
   sequences.pop_front(); // don't want the identity matrix!
  
   for (std::array<double, D> sequence : sequences) {
@@ -309,7 +310,7 @@ std::vector<Matrix<U, D>> torus_mats() {
   std::list<std::array<double, D>> sequences;
   sequences.push_back(ini_sequence);
 
-  one_sequences(sequences, D);
+  one_sequences<D>(sequences, D);
 
   sequences.pop_front(); // don't want the identity matrix!
 
@@ -375,7 +376,6 @@ public:
 
   void step(double T, unsigned ind, Euclidean<U, D> r) {
     unsigned cluster_size = 0;
-    std::uniform_real_distribution<double> dist(0.0, 1.0);
 
     std::queue<Spin<U, D, S>*> queue;
     queue.push(&(s[ind]));
@@ -396,7 +396,7 @@ public:
             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) {
+            if (rng.uniform(0.0, 1.0) < p) {
               queue.push(&ss);
             }
           }
@@ -423,7 +423,7 @@ public:
                 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) {
+              if (rng.uniform(0.0, 1.0) < p) {
                 queue.push(sj);
               }
             }
@@ -439,15 +439,13 @@ public:
   }
 
   void wolff(double T, unsigned N) {
-    std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1);
-
     for (unsigned i = 0; i < N; i++) {
-      R g = g(*this, rng);
-      unsigned ind = ind_dist(rng);
+      R r = g(*this, rng);
+      unsigned ind = rng.uniform((unsigned)0, (unsigned)(s.size() - 1));
 
-      A.pre_cluster(*this, ind, g);
+      A.pre_cluster(*this, ind, r);
 
-      this->step(T, ind, g, rng);
+      this->step(T, ind, r);
 
       A.post_cluster(*this);
     }
diff --git a/spheres.cpp b/spheres.cpp
index 9e67b87..4ed329e 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -37,17 +37,14 @@ 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) {
   return [&mats, &vecs] (const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
-    std::uniform_real_distribution<double> t_dist(0, M.L);
-    std::uniform_int_distribution<unsigned> r_dist(0, D - 1);
-    std::uniform_int_distribution<unsigned> ind_dist(0, M.s.size() - 1);
-    std::uniform_int_distribution<unsigned> coin(0, mats.size() + vecs.size() - 1);
-
     Vector<double, D> t;
     Matrix<double, D> m;
-    unsigned flip = coin(rng);
+    unsigned flip = rng.uniform((unsigned)0, (unsigned)(mats.size() + vecs.size() - 1));
     if (flip < mats.size()) {
+      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) = (double)t_dist(rng);
+        t(j) += rng.variate<double, std::normal_distribution>(0.0, 0.1);
       }
       m = mats[flip];
     } else {
@@ -120,7 +117,7 @@ int main(int argc, char* argv[]) {
 
   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);
+      return H * s.x(1);
     };
 
   std::vector<Matrix<double, D>> mats = torus_mats<double, D>();
@@ -129,16 +126,14 @@ int main(int argc, char* argv[]) {
   animation A(L, 750, argc, argv);
   model sphere(L, Z, B, g, std::floor(log2(L)), 2, A);
 
-  randutils::auto_seed_128 seeds;
-  std::mt19937 rng{seeds};
-
-  std::uniform_real_distribution<double> dist(0.0, L);
+  randutils::mt19937_rng rng;
 
   sphere.s.reserve(n);
 
+  unsigned nx = floor(sqrt(n));
   for (unsigned i = 0; i < n; i++) {
-    Vector<double, D> pos = {dist(rng), dist(rng)};
-    sphere.s.push_back({pos, dist(rng) / L});
+    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))};
+    sphere.s.push_back({pos, 0.5});
     sphere.dict.insert(&sphere.s.back());
   }
 
-- 
cgit v1.2.3-54-g00ecf