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

---
 ising.cpp       | 110 ++++++++++++++++++++++++++++++++++++++++----------------
 space_wolff.hpp |  42 +++++++++++++++-------
 spheres.cpp     |  47 ++++++++++++------------
 3 files changed, 132 insertions(+), 67 deletions(-)

diff --git a/ising.cpp b/ising.cpp
index ca2e2e7..fffce19 100644
--- a/ising.cpp
+++ b/ising.cpp
@@ -3,7 +3,7 @@
 
 std::function<double(spin<signed, 2, signed>)> B_sin(unsigned L, unsigned n, double H) {
   return [n, H, L] (spin<signed, 2, signed> s) -> double {
-    return H * cos(2 * M_PI * n * s.x[0] / ((double)L));
+    return H * s.s * cos(2 * M_PI * n * s.x(0) / ((double)L));
   };
 }
 
@@ -14,13 +14,14 @@ int main(int argc, char* argv[]) {
   unsigned N = 1000;
   unsigned mod = 0;
   unsigned multi = 1e4;
+  double mag = 0.5;
   double T = 2.0 / log(1.0 + sqrt(2.0));
   double H = 1.0;
   double ε = 0.1;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M")) != -1) {
+  while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M:n:")) != -1) {
     switch (opt) {
       case 'N': 
         N = (unsigned)atof(optarg);
@@ -43,42 +44,76 @@ int main(int argc, char* argv[]) {
       case 'M':
         multi = atoi(optarg);
         break;
+      case 'n':
+        mag = atof(optarg);
+        break;
       default:
         exit(1);
     }
   }
 
-  std::function<double(spin<signed, D, signed>, spin<signed, D, signed>)> Z =
-    [] (spin<signed, D, signed> s1, spin<signed, D, signed> s2) -> double {
-      bool one_one = false;
-      bool many_ones = false;
-      bool any_two = false;
+  double pZ = 1.0 - exp(-1.0 / T);
+
+  std::function<double(spin<signed, D, signed>, spin<signed, D, signed>, spin<signed, D, signed>)> Z =
+    [L, pZ] (spin<signed, D, signed> s1, spin<signed, D, signed> s2, spin<signed, D, signed> s1_new) -> double {
+      bool old_one_one = false;
+      bool old_many_ones = false;
+      bool old_any_two = false;
+      bool new_one_one = false;
+      bool new_many_ones = false;
+      bool new_any_two = false;
+
+      vector<signed, D> old_diff = diff<signed, D>(L, s1.x, s2.x);
+      vector<signed, D> new_diff = diff<signed, D>(L, s1_new.x, s2.x);
       
       for (unsigned i = 0; i < D; i++) {
-        unsigned diff = abs(s1.x(i) - s2.x(i));
-        if (diff == 1 && !one_one) {
-          one_one = true;
-        } else if (diff == 1 && one_one) {
-          many_ones = true;
-          break;
-        } else if (diff > 1) {
-          any_two = true;
-          break;
+        if (old_diff(i) == 1 && !old_one_one) {
+          old_one_one = true;
+        } else if (old_diff(i) == 1 && old_one_one) {
+          old_many_ones = true;
+        } else if (old_diff(i) > 1) {
+          old_any_two = true;
+        }
+        if (new_diff(i) == 1 && !new_one_one) {
+          new_one_one = true;
+        } else if (new_diff(i) == 1 && new_one_one) {
+          new_many_ones = true;
+        } else if (new_diff(i) > 1) {
+          new_any_two = true;
         }
       }
 
-      if (!one_one && !any_two) {
-        return -std::numeric_limits<double>::infinity();
-      } else if (one_one && !many_ones && !any_two) {
-        return s1.s * s2.s;
+      bool were_on_someone = !old_one_one && !old_any_two;
+      bool are_on_someone = !new_one_one && !new_any_two;
+      bool were_nearest_neighbors = old_one_one && !old_many_ones && !old_any_two;
+      bool are_nearest_neighbors = new_one_one && !new_many_ones && !new_any_two;
+
+      if (were_on_someone) {
+        return 0.0;
+      } else if (are_on_someone) {
+        return 1.0;
+      } else if (were_nearest_neighbors && are_nearest_neighbors) {
+        return 0.0;
+      } else if (were_nearest_neighbors) {
+        if (s1.s * s2.s == 1) {
+          return pZ;
+        } else {
+          return 0.0;
+        }
+      } else if (are_nearest_neighbors) {
+        if (s1_new.s * s2.s == -1) {
+          return pZ;
+        } else {
+          return 0.0;
+        }
       } else {
-        return 0;
+        return 0.0;
       }
     };
 
   std::function<double(spin<signed, D, signed>)> B_face =
     [L, H] (spin<signed, D, signed> s) -> double {
-      return H * s.s * smiley[s.x(1) * 16 / L][s.x(0) * 16 / L];
+      return H * s.s * smiley[s.x(0) * 16 / L][s.x(1) * 16 / L];
     };
 
   std::function<double(spin<signed, D, signed>)> B;
@@ -118,7 +153,7 @@ int main(int argc, char* argv[]) {
   unsigned down = 0;
   for (unsigned i = 0; i < L; i++) {
     for (unsigned j = 0; j < L; j++) {
-      if ((coin(rng) && up < pow(L, 2) / 2) || down >= pow(L, 2) / 2) {
+      if ((coin(rng) && up < pow(L, 2) * mag) || down >= pow(L, 2) * mag) {
         ising.s.push_back({{i, j}, 1});
         up++;
       } else {
@@ -154,13 +189,6 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  std::vector<signed> output(pow(L, D));
-
-  for (spin<signed, D, signed> s : ising.s) {
-    spin<signed, D, signed> rs = ising.s0.inverse().act(s);
-    output[L * rs.x(1) + rs.x(0)] = s.s;
-  }
-
   std::ofstream outfile;
   outfile.open("out.dat", std::ios::app);
 
@@ -173,6 +201,28 @@ int main(int argc, char* argv[]) {
   }
   outfile << "\n";
 
+  std::ofstream snapfile;
+  snapfile.open("snap.dat");
+
+  std::vector<std::vector<signed>> snap(L);
+  for (std::vector<signed>& line : snap) {
+    line.resize(L);
+  }
+
+  for (spin<signed, D, signed> s : ising.s) {
+    spin<signed, D, signed> snew = ising.s0.inverse().act(s);
+    snap[snew.x(0)][snew.x(1)] = snew.s;
+  }
+
+  for (std::vector<signed> row : snap) {
+    for (signed s : row) {
+      snapfile << s << " ";
+    }
+    snapfile << "\n";
+  }
+
+  snapfile.close();
+
   return 0;
 }
 
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);
               }
diff --git a/spheres.cpp b/spheres.cpp
index 106ce4e..b0f4fb4 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -34,20 +34,20 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  std::function<double(spin<double, D, double>, spin<double, D, double>)> Z =
-    [L] (spin<double, D, double> s1, spin<double, D, double> s2) -> double {
-      vector<double, D> diff = s1.x - s2.x;
-      for (unsigned i = 0; i < D; i++) {
-        if (fabs(diff(i)) > L / 2) {
-          diff(i) = L - fabs(diff(i));
-        } else {
-          diff(i) = fabs(diff(i));
-        }
-      }
-      if (diff.transpose() * diff < pow(s1.s + s2.s, 2)) {
-        return -std::numeric_limits<double>::infinity();
+  std::function<double(spin<double, D, double>, spin<double, D, double>, spin<double, D, double>)> Z =
+    [L] (spin<double, D, double> s1, spin<double, D, double> s2, spin<double, D, double> s1_new) -> double {
+      vector<double, D> diff_old = diff(L, s1.x, s2.x);
+      vector<double, D> diff_new = diff(L, s1_new.x, s2.x);
+
+      double rad_sum = pow(s1.s + s2.s, 2);
+
+      bool old_overlap = diff_old.transpose() * diff_old < rad_sum;
+      bool new_overlap = diff_new.transpose() * diff_new < rad_sum;
+
+      if (new_overlap) {
+        return 1.0;
       } else {
-        return 0;
+        return 0.0;
       }
     };
 
@@ -60,16 +60,10 @@ int main(int argc, char* argv[]) {
     [] (model<double, D, double>& m, unsigned i0, spin<double, D, double> s1) -> std::set<unsigned> {
       std::set<unsigned> nn;
       if (i0 < m.s.size()) {
-        std::set<unsigned> os1 = m.dict.on_site(s1.x);
-        std::set<unsigned> nn0 = m.dict.nearest_neighbors(m.s[i0].x);
-        std::set<unsigned> nn1 = m.dict.nearest_neighbors(s1.x);
-        std::set<unsigned> nnn0 = m.dict.next_nearest_neighbors(m.s[i0].x);
-        std::set<unsigned> nnn1 = m.dict.next_nearest_neighbors(s1.x);
-        nn.insert(nn0.begin(), nn0.end());
-        nn.insert(nn1.begin(), nn1.end());
-        nn.insert(nnn0.begin(), nnn0.end());
-        nn.insert(nnn1.begin(), nnn1.end());
-        nn.insert(os1.begin(), os1.end());
+        std::set<unsigned> n_old = m.dict.neighbors(m.s[i0].x, 1);
+        std::set<unsigned> n_new = m.dict.neighbors(s1.x, 1);;
+        nn.insert(n_old.begin(), n_old.end());
+        nn.insert(n_new.begin(), n_new.end());
         nn.insert(m.s.size());
       } else {
         for (unsigned i = 0; i < m.s.size(); i++) {
@@ -94,11 +88,16 @@ int main(int argc, char* argv[]) {
 
   sphere.wolff(T, N, rng);
 
+  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);
-    std::cout << s.x.transpose() << "\n";
+    snapfile << rs.x.transpose() << "\n";
   }
 
+  snapfile.close();
+
   return 0;
 }
 
-- 
cgit v1.2.3-54-g00ecf