tuning now progresses until convergence

2 files changed, 82 insertions, 81 deletions

diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
index 33dff71..f1cdf17 100644
--- a/hadamard_pt.hpp
+++ b/hadamard_pt.hpp
@@ -22,7 +22,6 @@ typedef struct range {
   unsigned N;
 } range;
 
-
 class PT {
 private:
   randutils::mt19937_rng rng;
@@ -41,6 +40,10 @@ public:
     }
   }
 
+  double T(unsigned i) {
+    return 1 / Ms[Ms.size() - i - 1].β;
+  }
+
   bool step(unsigned i, unsigned j, bool dry = false) {
     double Δβ = Ms[i].β - Ms[j].β;
     double ΔE = Ms[i].E - Ms[j].E;
@@ -63,105 +66,101 @@ public:
     }
   }
 
-  void tune(unsigned n, unsigned m, double ε) {
-    std::vector<color> colors(Ms.size(), none);
-    std::vector<unsigned> nu(Ms.size(), 0);
-    std::vector<unsigned> nd(Ms.size(), 0);
+  void tune(unsigned n0, unsigned m, double ε, double ε2) {
 
-    colors.front() = down;
-    colors.back() = up;
+    unsigned n = n0;
 
-    for (unsigned i = 0; i < n; i++) {
-#pragma omp parallel for
-      for (unsigned j = 0; j < Ms.size(); j++) {
-        Ms[j].tune(m, ε);
-      }
+    while (true) {
+      std::vector<color> colors(Ms.size(), none);
+      std::vector<unsigned> nu(Ms.size(), 0);
+      std::vector<unsigned> nd(Ms.size(), 0);
 
-      for (unsigned j = 0; j < Ms.size() - 1; j++) {
-        if (this->step(j, j + 1, true)) {
-          std::swap(colors[j], colors[j + 1]);
-          colors.front() = down;
-          colors.back() = up;
+      colors.front() = down;
+      colors.back() = up;
+      for (unsigned i = 0; i < n; i++) {
+#pragma omp parallel for
+        for (unsigned j = 0; j < Ms.size(); j++) {
+          Ms[j].tune(m, ε);
+        }
 
-          if (colors[j] == up) {
-            nu[j]++;
-          } else if (colors[j] == down) {
-            nd[j]++;
-          }
-          if (colors[j + 1] == up) {
-            nu[j + 1]++;
-          } else if (colors[j + 1] == down) {
-            nd[j + 1]++;
+        for (unsigned j = 0; j < Ms.size() - 1; j++) {
+          if (this->step(j, j + 1, true)) {
+            std::swap(colors[j], colors[j + 1]);
+            colors.front() = down;
+            colors.back() = up;
+
+            if (colors[j] == up) {
+              nu[j]++;
+            } else if (colors[j] == down) {
+              nd[j]++;
+            }
+            if (colors[j + 1] == up) {
+              nu[j + 1]++;
+            } else if (colors[j + 1] == down) {
+              nd[j + 1]++;
+            }
           }
         }
       }
-    }
-
-    std::ofstream outfile("test.dat", std::ios::app);
 
-    std::vector<double> f(Ms.size());
-
-    double δ = 0.5;
-    std::vector<unsigned> f_keep;
-    double f_last = 0;
-    for (unsigned i = 0; i < Ms.size(); i++) {
-      f[i] = nu[i] / (double)(nu[i] + nd[i]);
-      if (fabs(f[i] - (1 - (Ms.size() - i - 1) / (Ms.size() - 1.0))) < δ && f[i] >= f_last) {
-        f_keep.push_back(i);
-        f_last = f[i];
+      std::vector<double> f(Ms.size());
+
+      double err = 0;
+      double δ = 0.5;
+      std::vector<unsigned> f_keep;
+      double f_last = 0;
+      for (unsigned i = 0; i < Ms.size(); i++) {
+        f[i] = nu[i] / (double)(nu[i] + nd[i]);
+        double δerr = fabs(f[i] - (1 - (Ms.size() - i - 1) / (Ms.size() - 1.0)));
+        err += δerr;
+        if (δerr < δ && f[i] >= f_last) {
+          f_keep.push_back(i);
+          f_last = f[i];
+        }
       }
-    }
 
-    for (unsigned i = 0; i < f_keep.size() - 1; i++) {
-      for (unsigned j = f_keep[i]; j < f_keep[i + 1]; j++) {
-        f[j] = f[f_keep[i]] + (f[f_keep[i + 1]] - f[f_keep[i]]) / (Ms[f_keep[i+1]].β - Ms[f_keep[i]].β) * (Ms[j].β - Ms[f_keep[i]].β);
+      std::cout << err << "\n";
+      if (err < ε2) {
+        break;
       }
-    }
 
+      for (unsigned i = 0; i < f_keep.size() - 1; i++) {
+        for (unsigned j = f_keep[i]; j < f_keep[i + 1]; j++) {
+          f[j] = f[f_keep[i]] + (f[f_keep[i + 1]] - f[f_keep[i]]) / (Ms[f_keep[i+1]].β - Ms[f_keep[i]].β) * (Ms[j].β - Ms[f_keep[i]].β);
+        }
+      }
 
+      std::vector<double> η(Ms.size() - 1);
+      double C = 0;
 
-    std::vector<double> η(Ms.size() - 1);
-    double C = 0;
-
-    for (unsigned i = 0; i < Ms.size() - 1; i++) {
-      η[η.size() - i - 1] = sqrt((f[i + 1] - f[i]) / pow((1 / Ms[i].β - 1 / Ms[i + 1].β), 2));
-    }
-
-    std::vector<double> T0(Ms.size());
-    std::vector<double> T1(Ms.size());
+      for (unsigned i = 0; i < Ms.size() - 1; i++) {
+        η[η.size() - i - 1] = sqrt((f[i + 1] - f[i])) / (T(i) - T(i + 1));
+      }
 
-    for (unsigned i = 0; i < Ms.size(); i++) {
-      T0[i] = 1 / Ms[Ms.size() - i - 1].β;
-    }
-    T1.front() = 1 / Ms.back().β;
-    T1.back() = 1 / Ms.front().β;
+      std::vector<double> T1(Ms.size() - 2);
 
-    for (unsigned i = 0; i < η.size(); i++) {
-      C += η[i] * (T0[i + 1] - T0[i]);
-    }
+      for (unsigned i = 0; i < η.size(); i++) {
+        C += η[i] * (T(i + 1) - T(i));
+      }
 
-    double x= 0;
-    unsigned j = 1;
+      double x= 0;
+      unsigned j = 1;
 
-    for (unsigned i = 0; i < η.size(); i++) {
-      double xnew = x + η[i] * (T0[i + 1] - T0[i]) / C;
-      while (j < xnew * η.size()) {
-        T1[j] = T0[i] + (j / (double)η.size() - x) / η[i] * C;
-        j++;
+      for (unsigned i = 0; i < η.size(); i++) {
+        double xnew = x + η[i] * (T(i + 1) - T(i)) / C;
+        while (j < xnew * η.size() && j < η.size()) {
+          T1[j - 1] = T(i) + (j / (double)η.size() - x) / η[i] * C;
+          j++;
+        }
+        x = xnew;
       }
-      x = xnew;
-    }
 
-    for (unsigned i = 0; i < T1.size(); i++) {
-      Ms[Ms.size() - i - 1].β = 1 / T1[i];
-    }
+      for (unsigned i = 0; i < T1.size(); i++) {
+        Ms[Ms.size() - i - 2].β = 1 / T1[i];
+      }
 
-    for (double ff : f) {
-      outfile << ff << " ";
+      n *= 2;
     }
-    outfile << "\n";
-    outfile.close();
-
   }
 
   void run(unsigned n, unsigned m, bool dry = false) {
diff --git a/hadamard_test.cpp b/hadamard_test.cpp
index c1dd8cd..e505d79 100644
--- a/hadamard_test.cpp
+++ b/hadamard_test.cpp
@@ -3,7 +3,7 @@
 int main() {
   double β0 = 1.0;
   double β1 = 10.0;
-  unsigned N = 8;
+  unsigned N = 4;
   unsigned n = 4;
 
   ParallelMeasurement m1;
@@ -15,7 +15,9 @@ int main() {
 
   PT a(β0, β1, N, n, m1, m3);
 
-  while (true) {
-    a.tune(1000, 10, 0.1);
+  a.tune(1000, 10, 0.1, 0.1);
+
+  for (MCMC& M : a.Ms) {
+    std::cout << M.β;
   }
 }