added another qualification for finishing to minimze bias

2 files changed, 113 insertions, 31 deletions

diff --git a/ising.cpp b/ising.cpp
index 645c759..e4ab65c 100644
--- a/ising.cpp
+++ b/ising.cpp
@@ -1,18 +1,25 @@
 
 #include "space_wolff.hpp"
 
+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));
+  };
+}
+
 int main(int argc, char* argv[]) {
   const unsigned D = 2;
 
   unsigned L = 32;
   unsigned N = 1000;
+  unsigned mod = 0;
   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:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:L:T:H:e:m:")) != -1) {
     switch (opt) {
       case 'N': 
         N = (unsigned)atof(optarg);
@@ -29,6 +36,9 @@ int main(int argc, char* argv[]) {
       case 'e':
         ε = atof(optarg);
         break;
+      case 'm':
+        mod = atoi(optarg);
+        break;
       default:
         exit(1);
     }
@@ -62,11 +72,19 @@ int main(int argc, char* argv[]) {
       }
     };
 
-  std::function<double(spin<signed, D, signed>)> B =
+  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];
     };
 
+  std::function<double(spin<signed, D, signed>)> B;
+
+  if (mod > 0) {
+    B = B_sin(L, mod, H);
+  } else {
+    B = B_face;
+  }
+
   std::function<std::set<unsigned>(model<signed, D, signed>&, unsigned, spin<signed, D, signed>)> neighbors =
     [] (model<signed, D, signed>& m, unsigned i0, spin<signed, D, signed> s1) -> std::set<unsigned> {
       std::set<unsigned> nn;
@@ -127,7 +145,7 @@ int main(int argc, char* argv[]) {
     ising.wolff(T, N, rng);
     std::array<double, 2> τ = ising.Eq.τ();
     std::cout << τ[0] << " " << τ[1] << " " << τ[1] / τ[0] << "\n";
-    if (τ[1] / τ[0] < ε) {
+    if (τ[1] / τ[0] < ε && τ[0] * 1e4 < ising.Eq.num_added()) {
       break;
     }
   }
@@ -140,20 +158,11 @@ int main(int argc, char* argv[]) {
   }
 
   std::ofstream outfile;
-  outfile.open("snap.dat");
-
-  for (unsigned i = 0; i < L; i++) {
-    for (unsigned j = 0; j < L; j++) {
-      unsigned out = output[L * i + j] == 1 ? 1 : 0;
-      outfile << out << " ";
-    }
-    outfile << "\n";
-  }
-  outfile.close();
+  outfile.open("out.dat", std::ios::app);
 
   std::array<double, 2> act = ising.Eq.τ();
 
-  std::cout << act[0] << " " << act[1] << "\n";
+  outfile << L << " " << T << " " << mod << " " << H << " " << ising.Cq.avg() << " " << ising.Cq.serr() << " " << act[0] << " " << act[1] << "\n";
 
   return 0;
 }
diff --git a/space_wolff.hpp b/space_wolff.hpp
index b60e274..9a4b4b4 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -237,6 +237,10 @@ class quantity {
     double serr() const {
       return sqrt(this->σ());
     }
+
+    unsigned num_added() const {
+      return n - wait;
+    }
 };
 
 template <class U, int D, class state>
@@ -249,16 +253,87 @@ class model {
     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>)> B;
+    std::vector<matrix<U, D>> mats;
+    std::vector<vector<U, D>> steps;
     long double E;
     quantity Eq;
     quantity Cq;
 
+    void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
+      if (level > 0) {
+        unsigned new_level = level - 1;
+        unsigned old_length = sequences.size();
+        for (std::array<double, D>& sequence : sequences) {
+          std::array<double, D> new_sequence = sequence;
+          new_sequence[new_level] = -1;
+          sequences.push_front(new_sequence);
+        }
+        one_sequences(sequences, new_level);
+      }
+    }
+
     model(unsigned L, std::function<double(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) {
+        std::array<double, D> ini_sequence;
+        ini_sequence.fill(1);
+        std::list<std::array<double, D>> sequences;
+        sequences.push_back(ini_sequence);
+
+        one_sequences(sequences, D);
+
+        sequences.pop_front(); // don't want the identity matrix!
+
+        for (std::array<double, D> sequence : sequences) {
+          matrix<U, D> m;
+          for (unsigned i = 0; i < D; i++) {
+            for (unsigned j = 0; j < D; j++) {
+              if (i == j) {
+                m(i, j) = sequence[i];
+              } else {
+                m(i, j) = 0;
+              }
+            }
+          }
+
+          mats.push_back(m);
+
+          vector<U, D> v;
+          for (unsigned i = 0; i < D; i++) {
+            if (sequence[i] == 1) {
+              v(i) = 0;
+            } else {
+              v(i) = L / 2;
+            }
+          }
+
+          steps.push_back(v);
         }
 
+        for (unsigned i = 0; i < D; i++) {
+          for (unsigned j = 0; j < D; j++) {
+            if (i != j) {
+              matrix<U, D> m;
+              for (unsigned k = 0; k < D; k++) {
+                for (unsigned l = 0; l < D; l++) {
+                  if ((k == i && l == j) || (k == j && l == i)) {
+                    if (i < j) {
+                      m(k, l) = 1;
+                    } else {
+                      m(k, l) = -1;
+                    }
+                  } else {
+                    m(k, l) = 0;
+                  }
+                }
+              }
+              mats.push_back(m);
+            }
+          }
+        }
+      }
+
     void update_energy() {
       E = 0;
       for (unsigned i = 0; i < s.size(); i++) {
@@ -347,31 +422,29 @@ class model {
       std::uniform_real_distribution<double> t_dist(0, L);
       std::uniform_int_distribution<unsigned> r_dist(0, D - 1);
       std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1);
+      std::uniform_int_distribution<unsigned> coin(0, mats.size() + steps.size() - 1);
 
       for (unsigned i = 0; i < N; i++) {
         vector<U, D> t;
         matrix<U, D> m;
-        for (unsigned j = 0; j < D; j++) {
-          t(j) = (U)t_dist(rng);
-        }
-
-        unsigned flip_D1 = r_dist(rng);
-        unsigned flip_D2 = r_dist(rng);
-
-        for (unsigned j = 0; j < D; j++) {
-          for (unsigned k = 0; k < D; k++) {
-            if ((j == flip_D1 && k == flip_D2) || (j == flip_D2 && k == flip_D1)) {
-              if (flip_D1 <= flip_D2) {
-                m(j, k) = -1;
-              } else {
+        unsigned flip = coin(rng);
+        if (flip < mats.size()) {
+          for (unsigned j = 0; j < D; j++) {
+            t(j) = (U)t_dist(rng);
+          }
+          m = mats[flip];
+        } else {
+          for (unsigned j = 0; j < D; j++) {
+            for (unsigned k = 0; k < D; k++) {
+              if (j == k) {
                 m(j, k) = 1;
+              } else {
+                m(j, k) = 0;
               }
-            } else if ((j == k && j != flip_D1) && j != flip_D2) {
-              m(j, k) = 1;
-            } else {
-              m(j, k) = 0;
             }
           }
+
+          t = steps[flip - mats.size()];
         }
 
         euclidean<U, D> g(L, t, m);