diff options
| -rw-r--r-- | measurement.hpp | 20 | ||||
| -rw-r--r-- | space_wolff.hpp | 47 | ||||
| -rw-r--r-- | spheres_infinite.cpp | 84 | ||||
| -rw-r--r-- | transformation.hpp | 219 | 
4 files changed, 227 insertions, 143 deletions
| diff --git a/measurement.hpp b/measurement.hpp new file mode 100644 index 0000000..02fde9c --- /dev/null +++ b/measurement.hpp @@ -0,0 +1,20 @@ +#pragma once + +#include "spin.hpp" +#include "transformation.hpp" + +template <class U, int D, class R, class S> class Model; + +template <class U, int D, class R, class S> class measurement { +public: +  virtual void pre_cluster(const Model<U, D, R, S>&, unsigned, const Transformation<U, D, R, S>*){}; +  virtual void plain_bond_visited(const Model<U, D, R, 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, R, S>&, const Transformation<U, D, R, S>&){}; + +  virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&, +                                  const Spin<U, D, S>&, double){}; +  virtual void ghost_site_transformed(const Model<U, D, R, S>&, const R&){}; + +  virtual void post_cluster(const Model<U, D, R, S>&){}; +}; diff --git a/space_wolff.hpp b/space_wolff.hpp index 9d529a3..3314ef9 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -8,32 +8,15 @@  #include <vector>  #include "euclidean.hpp" +#include "measurement.hpp"  #include "octree.hpp"  #include "quantity.hpp" -#include "spin.hpp"  #include "random.hpp" +#include "spin.hpp"  #include "transformation.hpp" -template <class U, int D, class R, class S> class Model; - -template <class U, int D, class R, class S> class measurement { -public: -  virtual void pre_cluster(const Model<U, D, R, S>&, unsigned, const R&){}; -  virtual void plain_bond_visited(const Model<U, D, R, 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, R, S>&, const Spin<U, D, S>*, -                                      const Spin<U, D, S>&){}; - -  virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&, -                                  const Spin<U, D, S>&, double){}; -  virtual void ghost_site_transformed(const Model<U, D, R, S>&, const R&){}; - -  virtual void post_cluster(const Model<U, D, R, S>&){}; -}; -  template <class U, int D, class R, class S> class Model {  public: -  U L;    R s0;    std::vector<Spin<U, D, S>*> s;    Octree<U, D, S> dict; @@ -43,42 +26,42 @@ 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) -      : L(L), s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {} +      : s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {}    void step(double T, Transformation<U, D, R, S>* t0, measurement<U, D, R, S>& A) {      std::queue<Transformation<U, D, R, S>*> queue;      queue.push(t0); -    std::set<Spin<U, D, S>*> visited; +    std::set<Spin<U, D, S>*> cluster = t0->current();      while (!queue.empty()) {        Transformation<U, D, R, S>* t = queue.front();        queue.pop(); -      std::set<Spin<U, D, S>*> c = t->current(); -      if (!std::any_of(c.begin(), c.end(), [&visited](Spin<U, D, S>* s) { return visited.contains(s); })) { -        visited.insert(c.begin(), c.end()); - -        for (Spin<U, D, S>* s : t->toConsider()) { +      for (Spin<U, D, S>* s : t->toConsider()) { +        if (!cluster.contains(s)) {            double ΔE = t->ΔE(s);            if (ΔE > 0 && 1.0 - exp(-ΔE / T) > rng.uniform<double>(0, 1)) { -            queue.push(t->createNew(s)); +            cluster.insert(s); +            queue.push(t->createNew(T, cluster, s, rng));            }          } - -        t->apply();        } + +      A.plain_site_transformed(*this, *t); +      t->apply();        delete t;      }    }    void resize(double T, double P) {} -  void wolff(double T, std::vector<Gen<U, D, R, S>> gs, -             measurement<U, D, R, S>& A, unsigned N) { +  void wolff(double T, std::vector<Gen<U, D, R, S>> gs, measurement<U, D, R, S>& A, unsigned N) {      for (unsigned i = 0; i < N; i++) {        Gen<U, D, R, S>& g = rng.pick(gs); -      Transformation<U, D, R, S> *t = g(*this, rng); +      Transformation<U, D, R, S>* t = g(*this, rng); + +      A.pre_cluster(*this, i, t);        this->step(T, t, A); diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp index 5287a84..e41ea17 100644 --- a/spheres_infinite.cpp +++ b/spheres_infinite.cpp @@ -14,9 +14,11 @@ private:    uint64_t t2;    unsigned n;    unsigned tmp; +  unsigned wait; +  Euclidean<double, D> s0_tmp;  public: -  animation(double L, unsigned w, int argc, char* argv[]) { +  animation(double L, unsigned w, int argc, char* argv[]) : s0_tmp(0), wait(1000) {      t1 = 0;      t2 = 0;      n = 0; @@ -24,37 +26,65 @@ public:      glutInit(&argc, argv);      glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);      glutInitWindowSize(w, w); -    glutCreateWindow("wolff"); +    glutCreateWindow("wolffWindow");      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 { +  void plain_site_transformed(const model& m, const Transformation<double, D, Euclidean<double, D>, double>& t) override { +    if (n > wait) { +    if (t.current().size() > 1) { +      glColor3f(0.0f, 0.0f, 1.0f); +    } else { +      glColor3f(0.0f, 1.0f, 0.0f); +    } +    for (const Spin<double, 2, double>* s : t.current()) { +      if (s != NULL) { +      glBegin(GL_POLYGON); +      unsigned n_points = 50; +      for (unsigned i = 0; i < n_points; i++) { +        glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points), +                   s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points)); +      } +      glEnd(); +      } +    } +    }      tmp++;    } -  void post_cluster(const model& m) override { +  void pre_cluster(const model& m, unsigned, const Transformation<double, D, Euclidean<double, D>, double>* t) override { +    s0_tmp = m.s0.inverse(); +    tmp = 0;      glClearColor(1.0f, 1.0f, 1.0f, 1.0f);      glClear(GL_COLOR_BUFFER_BIT); +    /* +    glBegin(GL_LINE); +    Vector<double, 2> r_center = s0_tmp.act({t->r.t, 0}); +    glEnd(); +    */      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(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points), +                   s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points));        }        glEnd();      }      glFlush(); +  } +  void post_cluster(const model& m) override { +    glFlush();      t1 += tmp;      t2 += tmp * tmp; +    if (n > wait) { +      sleep(2); +    }      n++;    } @@ -105,18 +135,47 @@ Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {      Vector<double, D> t12 = t1 - t2;      Vector<double, D> t = (t1 + t2) / 2; -    double θ = atan2(t12[1], t12[0]) + rng.variate<double, std::normal_distribution>(0.0, ε); +    double θ = atan2(t12(1), t12(0)) + rng.variate<double, std::normal_distribution>(0.0, ε) / t12.norm();      m(0, 0) = -cos(2 * θ);      m(1, 1) = cos(2 * θ);      m(0, 1) = -2 * cos(θ) * sin(θ);      m(1, 0) = -2 * cos(θ) * sin(θ); -    Euclidean<double, D> g(t - m * t, m); +    Vector<double, D> t3 = t - m * t; + +    if (t3(0) != t3(0)) { +      std::cout << t3 << "\n" << t << "\n" << m << "\n" << t12 << "\n" ; +      getchar(); +    } + +    Euclidean<double, D> g(t3, m);      return new PairFlip<double, D, Euclidean<double, D>, double>(M, g, s1, s2);    };  } +Gen<double, D, Euclidean<double, D>, double> rGen(double ε) { +  return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* { +    Vector<double, D> t; +    Matrix<double, D> m; + +    double θ = rng.uniform((double)0.0, 2 * M_PI); +    m(0, 0) = -cos(2 * θ); +    m(1, 1) = cos(2 * θ); +    m(0, 1) = -2 * cos(θ) * sin(θ); +    m(1, 0) = -2 * cos(θ) * sin(θ); + +    Spin<double, D, double>* s = rng.pick(M.s); +    t = M.s0.t; +    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, s); +  }; +} +  int main(int argc, char* argv[]) {    const unsigned D = 2; @@ -180,7 +239,8 @@ int main(int argc, char* argv[]) {    };    auto g1 = eGen(1); -  auto g2 = mGen(0.5); +  auto g2 = mGen(0.1); +  auto g3 = rGen(1);    animation A(L, 750, argc, argv);    model sphere(1.0, Z, B); diff --git a/transformation.hpp b/transformation.hpp index 33ac50f..b37b79e 100644 --- a/transformation.hpp +++ b/transformation.hpp @@ -9,25 +9,63 @@ template <class U, int D, class R, class S> class Model;  template <class U, int D, class R, class S> class Transformation {  public: +  Model<U, D, R, S>& M; +  const R r; + +  Transformation(Model<U, D, R, S>& M, const R& r) : r(r), M(M) {}    virtual std::set<Spin<U, D, S>*> current() const { return {NULL}; }    virtual std::set<Spin<U, D, S>*> toConsider() const { return {}; }    virtual double ΔE(const Spin<U, D, S>*) const { return 0; } -  virtual Transformation* createNew(Spin<U, D, S>*) const { return new Transformation(); } +  virtual Transformation* createNew(double T, std::set<Spin<U, D, S>*>&, Spin<U, D, S>*, +                                    Rng& rng) const { +    return new Transformation(M, r); +  }    virtual void apply(){};  };  template <class U, int D, class R, class S>  using Gen = std::function<Transformation<U, D, R, S>*(Model<U, D, R, S>&, Rng&)>; +template <class U, int D, class R, class S> class FieldFlip : public Transformation<U, D, R, S> { +private: +  R s0New; + +public: +  FieldFlip(Model<U, D, R, S>& M, const R& r) +      : Transformation<U, D, R, S>(M, r), s0New(r.act(M.s0)) {} + +  std::set<Spin<U, D, S>*> toConsider() const override { +    std::set<Spin<U, D, S>*> neighbors; + +    for (Spin<U, D, S>* s : this->M.s) { +      neighbors.insert(s); +    } + +    return neighbors; +  } + +  double ΔE(const Spin<U, D, S>* s) const override { +    Spin<U, D, S> s0s_old = this->M.s0.inverse().act(*s); +    Spin<U, D, S> s0s_new = s0New.inverse().act(*s); +    return this->M.B(s0s_new) - this->M.B(s0s_old); +  } + +  Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& cluster, +                                             Spin<U, D, S>* s, Rng& rng) const override { +    return pairGenerator<D, R, S>(this, T, cluster, s, rng); +  } + +  void apply() override { this->M.s0 = s0New; } +}; +  template <class U, int D, class R, class S> class SpinFlip : public Transformation<U, D, R, S> {  private: -  Model<U, D, R, S>& M;    Spin<U, D, S>* sOld; -  const R r;    Spin<U, D, S> sNew;  public: -  SpinFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s) : M(M), r(r) { +  SpinFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s) +      : Transformation<U, D, R, S>(M, r) {      sOld = s;      sNew = r.act(*s);    } @@ -35,47 +73,54 @@ public:    std::set<Spin<U, D, S>*> current() const override { return {sOld}; }    std::set<Spin<U, D, S>*> toConsider() const override { -    std::set<Spin<U, D, S>*> neighbors; -    std::set<Spin<U, D, S>*> current_neighbors = M.dict.neighbors(sOld->x); -    std::set<Spin<U, D, S>*> new_neighbors = M.dict.neighbors(sNew.x); +    std::set<Spin<U, D, S>*> neighbors = {NULL}; +    std::set<Spin<U, D, S>*> oldNeighbors = this->M.dict.neighbors(sOld->x); +    std::set<Spin<U, D, S>*> newNeighbors = this->M.dict.neighbors(sNew.x); + +    neighbors.insert(oldNeighbors.begin(), oldNeighbors.end()); +    neighbors.insert(newNeighbors.begin(), newNeighbors.end()); -    neighbors.insert(current_neighbors.begin(), current_neighbors.end()); -    neighbors.insert(new_neighbors.begin(), new_neighbors.end()); -    neighbors.insert(NULL); +    neighbors.erase(sOld);      return neighbors;    }    double ΔE(const Spin<U, D, S>* s) const override {      if (s == NULL) { -      Spin<U, D, S> s0s_old = M.s0.inverse().act(*sOld); -      Spin<U, D, S> s0s_new = M.s0.inverse().act(sNew); -      return M.B(s0s_new) - M.B(s0s_old); +      Spin<U, D, S> s0s_old = this->M.s0.inverse().act(*sOld); +      Spin<U, D, S> s0s_new = this->M.s0.inverse().act(sNew); +      return this->M.B(s0s_new) - this->M.B(s0s_old);      } else { -      return M.Z(*sOld, *s) - M.Z(sNew, *s); +      return this->M.Z(*sOld, *s) - this->M.Z(sNew, *s);      }    } -  Transformation<U, D, R, S>* createNew(Spin<U, D, S>* s) const override; +  Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& cluster, +                                             Spin<U, D, S>* s, Rng& rng) const override { +    if (s == NULL) { +      return new FieldFlip<U, D, R, S>(this->M, this->r); +    } else { +      return pairGenerator<D, R, S>(this, T, cluster, s, rng); +    } +  }    void apply() override { -    M.dict.remove(sOld); +    this->M.dict.remove(sOld);      *sOld = sNew; -    M.dict.insert(sOld); +    this->M.dict.insert(sOld);    }  };  template <class U, int D, class R, class S> class PairFlip : public Transformation<U, D, R, S> {  private: -  Model<U, D, R, S>& M;    Spin<U, D, S>* s1Old;    Spin<U, D, S>* s2Old; -  const R r;    Spin<U, D, S> s1New;    Spin<U, D, S> s2New;  public: -  PairFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s1, Spin<U, D, S>* s2) : M(M), r(r) { +  PairFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s1, Spin<U, D, S>* s2) +      : Transformation<U, D, R, S>(M, r) {      s1Old = s1;      s2Old = s2;      s1New = r.act(*s1); @@ -86,115 +131,91 @@ public:    std::set<Spin<U, D, S>*> toConsider() const override {      std::set<Spin<U, D, S>*> neighbors = {NULL}; -    std::set<Spin<U, D, S>*> current_neighbors_1 = M.dict.neighbors(s1Old->x); -    std::set<Spin<U, D, S>*> current_neighbors_2 = M.dict.neighbors(s2Old->x); -    std::set<Spin<U, D, S>*> new_neighbors_1 = M.dict.neighbors(s1New.x); -    std::set<Spin<U, D, S>*> new_neighbors_2 = M.dict.neighbors(s2New.x); + +    std::set<Spin<U, D, S>*> current_neighbors_1 = this->M.dict.neighbors(s1Old->x); +    std::set<Spin<U, D, S>*> current_neighbors_2 = this->M.dict.neighbors(s2Old->x); +    std::set<Spin<U, D, S>*> new_neighbors_1 = this->M.dict.neighbors(s1New.x); +    std::set<Spin<U, D, S>*> new_neighbors_2 = this->M.dict.neighbors(s2New.x);      neighbors.insert(current_neighbors_1.begin(), current_neighbors_1.end());      neighbors.insert(current_neighbors_2.begin(), current_neighbors_2.end());      neighbors.insert(new_neighbors_1.begin(), new_neighbors_1.end());      neighbors.insert(new_neighbors_2.begin(), new_neighbors_2.end()); +    neighbors.erase(s1Old); +    neighbors.erase(s2Old); +      return neighbors;    }    double ΔE(const Spin<U, D, S>* s) const override {      if (s == NULL) { -      Spin<U, D, S> s0s1_old = M.s0.inverse().act(*s1Old); -      Spin<U, D, S> s0s1_new = M.s0.inverse().act(s1New); -      Spin<U, D, S> s0s2_old = M.s0.inverse().act(*s2Old); -      Spin<U, D, S> s0s2_new = M.s0.inverse().act(s2New); -      return M.B(s0s1_new) + M.B(s0s2_new) - M.B(s0s1_old) - M.B(s0s2_old); +      Spin<U, D, S> s0s1_old = this->M.s0.inverse().act(*s1Old); +      Spin<U, D, S> s0s1_new = this->M.s0.inverse().act(s1New); +      Spin<U, D, S> s0s2_old = this->M.s0.inverse().act(*s2Old); +      Spin<U, D, S> s0s2_new = this->M.s0.inverse().act(s2New); +      return this->M.B(s0s1_new) + this->M.B(s0s2_new) - this->M.B(s0s1_old) - this->M.B(s0s2_old);      } else { -      return M.Z(*s1Old, *s) + M.Z(*s2Old, *s) - M.Z(s1New, *s) - M.Z(s2New, *s); +      return this->M.Z(*s1Old, *s) + this->M.Z(*s2Old, *s) - this->M.Z(s1New, *s) - this->M.Z(s2New, *s);      }    } -  Transformation<U, D, R, S>* createNew(Spin<double, D, S>* s) const; -  Transformation<U, D, R, S>* createNew(Spin<signed, D, S>* s) const; +  Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& cluster, +                                             Spin<U, D, S>* s, Rng& rng) const override { +    if (s == NULL) { +      return new FieldFlip<U, D, R, S>(this->M, this->r); +    } else { +      return pairGenerator<D, R, S>(this, T, cluster, s, rng); +    } +  }    void apply() override { -    M.dict.remove(s1Old); -    M.dict.remove(s2Old); +    this->M.dict.remove(s1Old); +    this->M.dict.remove(s2Old);      *s1Old = s1New;      *s2Old = s2New; -    M.dict.insert(s1Old); -    M.dict.insert(s2Old); +    this->M.dict.insert(s1Old); +    this->M.dict.insert(s2Old);    }  }; -template <class U, int D, class R, class S> class FieldFlip : public Transformation<U, D, R, S> { -private: -  Model<U, D, R, S>& M; -  const R r; -  R s0New; - -public: -  FieldFlip(Model<U, D, R, S>& M, const R& r) : M(M), r(r), s0New(r.act(M.s0)) {} - -  std::set<Spin<U, D, S>*> toConsider() const override { -    std::set<Spin<U, D, S>*> neighbors; -    for (Spin<U, D, S>* s : M.s) { -      neighbors.insert(s); +template <int D, class R, class S> +Transformation<double, D, R, S>* +pairGenerator(const Transformation<double, D, R, S>* tOld, double T, +                        std::set<Spin<double, D, S>*>& cluster, Spin<double, D, S>* s, Rng& rng) { +  SpinFlip<double, D, R, S>* t = new SpinFlip<double, D, R, S>(tOld->M, tOld->r, s); +  Spin<double, D, S>* sMax = NULL; +  double ΔEMax = 0.0; +  for (Spin<double, D, S>* ss : t->toConsider()) { +    if (ss != NULL && ss != s && !cluster.contains(ss)) { +      double ΔE = t->ΔE(ss); +      if (ΔE > ΔEMax) { +        sMax = ss; +        ΔEMax = ΔE; +      }      } - -    return neighbors;    } - -  double ΔE(const Spin<U, D, S>* s) const override { -    Spin<U, D, S> s0s_old = M.s0.inverse().act(*s); -    Spin<U, D, S> s0s_new = s0New.inverse().act(*s); -    return M.B(s0s_new) - M.B(s0s_old); -  } - -  Transformation<double, D, R, S>* createNew(Spin<double, D, S>* s) const { -    return new SpinFlip<U, D, R, S>(M, r, s); -  } - -  Transformation<signed, D, R, S>* createNew(Spin<signed, D, S>* s) const { -    Vector<signed, 2> v = r.act(*s).x; -    std::set<Spin<U, D, S>*> on_site = M.dict.at(v); -    if (on_site.empty()) { -      return new SpinFlip<U, D, R, S>(M, r, s); -    } else { -      return new PairFlip<U, D, R, S>(M, r, s, *on_site.begin()); -    } -  } - -  void apply() override { M.s0 = s0New; } -}; - -template <class U, int D, class R, class S> -Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(Spin<U, D, S>* s) const { -  if (s == NULL) { -    return new FieldFlip<U, D, R, S>(M, r); +  if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1) || cluster.contains(sMax)) { +    return t;    } else { -    return new SpinFlip(M, r, s); +    delete t; +    cluster.insert(sMax); +    return new PairFlip<double, D, R, S>(tOld->M, tOld->r, s, sMax);    }  } -template <class U, int D, class R, class S> -Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(Spin<double, D, S>* s) const { -  if (s == NULL) { -    return new FieldFlip<U, D, R, S>(M, r); +template <int D, class R, class S> +Transformation<signed, D, R, S>* +pairGenerator(const Transformation<signed, D, R, S>* tOld, double T, +                        std::set<Spin<signed, D, S>*>& cluster, Spin<signed, D, S>* s, Rng& rng) { +  Vector<signed, 2> v = tOld->r.act(*s).x; +  std::set<Spin<signed, D, S>*> on_site = tOld->M.dict.at(v); +  if (on_site.empty()) { +    return new SpinFlip<signed, D, R, S>(tOld->M, tOld->r, s);    } else { -    return new SpinFlip<U, D, R, S>(M, r, s); +    cluster.insert(*on_site.begin()); +    return new PairFlip<signed, D, R, S>(tOld->M, tOld->r, s, *on_site.begin());    }  } -template <class U, int D, class R, class S> -Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(Spin<signed, D, S>* s) const { -  if (s == NULL) { -    return new FieldFlip<U, D, R, S>(M, r); -  } else { -    Vector<signed, 2> v = r.act(*s).x; -    std::set<Spin<U, D, S>*> on_site = M.dict.at(v); -    if (on_site.empty()) { -      return new SpinFlip<U, D, R, S>(M, r, s); -    } else { -      return new PairFlip<U, D, R, S>(M, r, s, *on_site.begin()); -    } -  } -} | 
