2 files changed, 72 insertions, 7 deletions

diff --git a/space_wolff.hpp b/space_wolff.hpp
index 196b3b5..4c935db 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -31,6 +31,42 @@ public:
   virtual void post_cluster(const Model<U, D, R, S>&){};
 };
 
+template <class U, int D, class R, class S> using Gen = std::function<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)>;
+
+template <class U, int D, class R, class S> class Model;
+
+template <class U, int D, class R, class S> class Transformation {
+  public:
+    virtual void ready() {};
+    virtual double p(const Transformation&) const { return 0; }
+    virtual void apply() {};
+    virtual std::set<Transformation*> to_consider() const {};
+};
+
+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>& s;
+    const R& r;
+    Spin<U, D, S> s_new;
+  public:
+    SpinFlip(Model<U, D, R, S>& M, Spin<U, D, S>& s, const R& r) : M(M), s(s), r(r) {}
+
+    void ready() override {
+      s_new = r.act(s);
+    };
+
+    double p(const SpinFlip& sf, double T) const override {
+      return 1.0 - exp(-(M.Z(s, sf.s) - M.Z(s_new, sf.s)) / T);
+    }
+
+    double p(const SpinFlip& sf, double T) const override {
+      return 1.0 - exp(-(M.Z(s, sf.s) - M.Z(s_new, sf.s)) / T);
+    }
+
+
+};
+
 template <class U, int D, class R, class S> class Model {
 public:
   U L;
@@ -108,9 +144,10 @@ public:
 
   void resize(double T, double P) {}
 
-  void wolff(double T, std::function<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)> g,
+  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);
       R r = g(*this, rng);
       unsigned ind = rng.uniform((unsigned)0, (unsigned)(s.size() - 1));
 
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index 9bc86d7..a9cd7d0 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -67,7 +67,7 @@ public:
   double var() { return (t2 - t1 * t1 / (double)n) / (double)n; }
 };
 
-std::function<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(double ε) {
+Gen<double, D, Euclidean<double, D>, double> eGen(double ε) {
   return [ε](const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
     Vector<double, D> t;
     Matrix<double, D> m;
@@ -89,6 +89,31 @@ std::function<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(
   };
 }
 
+Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
+  return [ε](const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
+    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++;
+
+    Vector<double, D> t1 = M.s[f_ind1].x;
+    Vector<double, D> t2 = M.s[f_ind2].x;
+    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, ε);
+
+    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);
+    return g;
+  };
+}
+
 int main(int argc, char* argv[]) {
   const unsigned D = 2;
 
@@ -122,8 +147,8 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  double k = 1e8;
-  double a = 0.0;
+  double k = 1000;
+  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 {
@@ -145,7 +170,8 @@ int main(int argc, char* argv[]) {
     return H * s.x.norm();
   };
 
-  auto g = eGen(0.05);
+  auto g1 = eGen(0.25);
+  auto g2 = mGen(0.1);
   animation A(L, 750, argc, argv);
   model sphere(1, Z, B);
 
@@ -156,14 +182,15 @@ int main(int argc, char* argv[]) {
   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.s.push_back({pos, rng.uniform<double>(0.45, 0.45)});
     sphere.dict.insert(&sphere.s.back());
   }
 
-  sphere.wolff(T, g, A, N);
+  sphere.wolff(T, {g1, g2}, 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);
@@ -172,6 +199,7 @@ int main(int argc, char* argv[]) {
     outfile << ε << " " << A.var() / sphere.s.size() << std::endl;
     std::cout << ε << " " << A.var() / sphere.s.size() << std::endl;
   }
+  */
 
   outfile.close();