From 1586d6b2247c510c11d9a4847bde67cb0947243d Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 14 Aug 2019 11:45:12 -0400
Subject: fixed bugs in Ising, p for bonds pre-computed

---
 space_wolff.hpp | 42 +++++++++++++++++++++++++++++-------------
 1 file changed, 29 insertions(+), 13 deletions(-)

(limited to 'space_wolff.hpp')

diff --git a/space_wolff.hpp b/space_wolff.hpp
index 9a4b4b4..c87697e 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -7,6 +7,7 @@
 #include <functional>
 #include <random>
 #include <queue>
+#include <cmath>
 #include <eigen3/Eigen/Dense>
 #include "randutils/randutils.hpp"
 
@@ -36,6 +37,20 @@ using vector = Eigen::Matrix<T, D, 1>;
 template <class T, int D>
 using matrix = Eigen::Matrix<T, D, D>;
 
+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));
+    if (v(i) > L / 2) {
+      v(i) = L - v(i);
+    }
+  }
+
+  return v;
+}
+
 template <class T, int D, class state>
 class spin {
   public:
@@ -46,12 +61,12 @@ class spin {
 template <class T, int D>
 class euclidean {
   private:
-    unsigned L;
+    T L;
 
   public:
     vector<T, D> t;
     matrix<T, D> r;
-    euclidean(unsigned L) : L(L) {
+    euclidean(T L) : L(L) {
       for (unsigned i = 0; i < D; i++) {
         t(i) = 0;
         r(i, i) = 1;
@@ -61,7 +76,7 @@ class euclidean {
       }
     }
 
-    euclidean(unsigned L, vector<T, D> t0, matrix<T, D> r0) : L(L) {
+    euclidean(T L, vector<T, D> t0, matrix<T, D> r0) : L(L) {
       t = t0;
       r = r0;
     }
@@ -84,6 +99,10 @@ class euclidean {
       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);
+      }
+
       euclidean pnew(this->L, tnew, rnew);
 
       return pnew;
@@ -99,7 +118,6 @@ class euclidean {
     }
 };
 
-// to-do: clean up nnn... code
 template <int D>
 class dictionary {
   private:
@@ -246,12 +264,12 @@ class quantity {
 template <class U, int D, class state>
 class model {
   public:
-    unsigned L;
+    U L;
     euclidean<U, D> s0;
     std::vector<spin<U, D, state>> s;
     dictionary<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>, spin<U, D, state>, spin<U, D, state>)> Z;
     std::function<double(spin<U, D, state>)> B;
     std::vector<matrix<U, D>> mats;
     std::vector<vector<U, D>> steps;
@@ -272,7 +290,7 @@ class model {
       }
     }
 
-    model(unsigned L, std::function<double(spin<U, D, state>, spin<U, D, state>)> Z,
+    model(U L, std::function<double(spin<U, D, state>, 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(L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000), Cq(1000, 1000) {
@@ -338,7 +356,7 @@ class model {
       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]));
@@ -374,7 +392,7 @@ class model {
 
           for (unsigned j : neighbors(*this, i, si_new)) {
             if (j != i) {
-              double dE;
+              double p;
               bool neighbor_is_ghost = j == s.size();
 
               if (we_are_ghost || neighbor_is_ghost) {
@@ -391,13 +409,11 @@ class model {
                   s0s_new = s0_new.inverse().act(s[j]);
                 }
 
-                dE = B(s0s_old) - B(s0s_new);
+                p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T);
               } else {
-                dE = Z(s[i], s[j]) - Z(si_new, s[j]);
+                p = Z(s[i], s[j], si_new);
               }
 
-              double p = 1.0 - exp(-dE / T);
-
               if (dist(rng) < p) {
                 queue.push(j);
               }
-- 
cgit v1.2.3-54-g00ecf