1 files changed, 26 insertions, 20 deletions

diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index 3505d6d..026486e 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -41,13 +41,13 @@ public:
   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) {
+    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));
+        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();
     }
@@ -78,14 +78,14 @@ Gen<double, D, Euclidean<double, D>, double> eGen(double ε) {
     m(0, 1) = -2 * cos(θ) * sin(θ);
     m(1, 0) = -2 * cos(θ) * sin(θ);
 
-    unsigned f_ind = rng.uniform((unsigned)0, (unsigned)M.s.size());
-    t =  M.s[f_ind].x;
+    Spin<double, D, double>* s = rng.pick(M.s);
+    t = s->x;
     for (unsigned j = 0; j < D; j++) {
       t(j) += rng.variate<double, std::normal_distribution>(0.0, ε);
     }
 
     Euclidean<double, D> g(t - m * t, m);
-    return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, M.s[f_ind]);
+    return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, s);
   };
 }
 
@@ -93,12 +93,15 @@ Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
   return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
     Matrix<double, D> m;
 
-    unsigned f_ind1 = rng.uniform((unsigned)0, (unsigned)M.s.size());
-    unsigned f_ind2 = rng.uniform((unsigned)0, (unsigned)M.s.size() - 1);
-    if (f_ind2 >= f_ind1) f_ind2++;
+    Spin<double, D, double>* s1 = rng.pick(M.s);
+    Spin<double, D, double>* s2 = rng.pick(M.s);
 
-    Vector<double, D> t1 = M.s[f_ind1].x;
-    Vector<double, D> t2 = M.s[f_ind2].x;
+    while (s1 == s2) {
+      s2 = rng.pick(M.s);
+    }
+
+    Vector<double, D> t1 = s1->x;
+    Vector<double, D> t2 = s2->x;
     Vector<double, D> t12 = t1 - t2;
     Vector<double, D> t = (t1 + t2) / 2;
 
@@ -110,7 +113,7 @@ Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
     m(1, 0) = -2 * cos(θ) * sin(θ);
 
     Euclidean<double, D> g(t - m * t, m);
-    return new PairFlip<double, D, Euclidean<double, D>, double>(M, g, M.s[f_ind1], M.s[f_ind2]);
+    return new PairFlip<double, D, Euclidean<double, D>, double>(M, g, s1, s2);
   };
 }
 
@@ -179,17 +182,19 @@ int main(int argc, char* argv[]) {
   auto g1 = eGen(0.5);
   auto g2 = mGen(0.2);
   animation A(L, 750, argc, argv);
-  model sphere(1, Z, B);
+  model sphere(1.0, Z, B);
 
   Rng rng;
 
-  sphere.s.reserve(n);
+  sphere.s.resize(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, rng.uniform<double>(0.25, 0.45)});
-    sphere.dict.insert(&sphere.s.back());
+  for (unsigned i = 0; i < sphere.s.size(); i++) {
+    Spin<double, 2, double>* ss = new Spin<double, 2, double>();
+    ss->x = {(i / nx) * L / nx, (i % nx) * L / nx};
+    ss->s = rng.uniform<double>(0.25, 0.45);
+    sphere.s[i] = ss;
+    sphere.dict.insert(ss);
   }
 
   sphere.wolff(T, {g1, g2}, A, N);
@@ -212,9 +217,10 @@ int main(int argc, char* argv[]) {
   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";
+    delete s;
   }
 
   snapfile.close();