half implemented new tuning for pt

3 files changed, 132 insertions, 41 deletions

diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp
index c9fa2da..145454b 100644
--- a/hadamard_mcmc.hpp
+++ b/hadamard_mcmc.hpp
@@ -153,17 +153,13 @@ private:
   double θ0;
 
 public:
-  const double β;
+  double β;
   double E;
   Orthogonal M;
-  unsigned tag;
-  color c;
 
-  MCMC(unsigned n, double β0, Measurement& A, unsigned tag = 0) : A(A), M(n), β(β0) {
+  MCMC(unsigned n, double β0, Measurement& A) : A(A), M(n), β(β0) {
     θ0 = M_PI;
     E = M.energy();
-    tag = tag;
-    c = none;
   }
 
   bool step(Givens& g, bool dry = false) {
diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
index fba97b5..33dff71 100644
--- a/hadamard_pt.hpp
+++ b/hadamard_pt.hpp
@@ -2,12 +2,12 @@
 #pragma once
 #include "hadamard_mcmc.hpp"
 #include <list>
+#include <fstream>
+#include <iostream>
 
 void swap(MCMC& s1, MCMC& s2) {
   std::swap(s1.M, s2.M);
   std::swap(s1.E, s2.E);
-  std::swap(s1.tag, s2.tag);
-  std::swap(s1.c, s2.c);
 }
 
 class ParallelMeasurement {
@@ -22,6 +22,7 @@ typedef struct range {
   unsigned N;
 } range;
 
+
 class PT {
 private:
   randutils::mt19937_rng rng;
@@ -31,41 +32,13 @@ public:
   ParallelMeasurement& B;
   std::vector<Measurement*>& As;
 
-  PT(double β₀, double β₁, unsigned n, ParallelMeasurement& B, std::vector<Measurement*>& As)
+  PT(double β₀, double β₁, unsigned N, unsigned n, ParallelMeasurement& B, std::vector<Measurement*>& As)
       : B(B), As(As) {
-    Ms.reserve(n);
-    for (unsigned i = 1; i <= n; i++) {
-      double β = β₀ + i * (β₁ - β₀) / n;
-      Ms.push_back(MCMC(n, β, *As[i - 1], i - 1));
-    }
-    Ms[0].c = down;
-    Ms[n - 1].c = up;
-  }
-
-  void tune(unsigned n, unsigned m, double ε) {
-    std::vector<unsigned> nu(Ms.size());
-    std::vector<unsigned> nd(Ms.size());
-
-    for (unsigned i = 0; i < n; i++) {
-      Ms[0].c = down;
-      Ms[Ms.size() - 1].c = up;
-
-#pragma omp parallel for
-      for (unsigned j = 0; j < Ms.size(); j++) {
-        Ms[j].tune(m, ε);
-      }
-      this->sweep(true);
-
-      for (unsigned j = 0; j < Ms.size(); j++) {
-        if (Ms[j].c == up) {
-          nu[j]++;
-        } else if (Ms[j].c == down) {
-          nd[j]++;
-        }
-      }
+    Ms.reserve(N);
+    for (unsigned i = 0; i < N; i++) {
+      double β = β₀ + i * (β₁ - β₀) / (N - 1);
+      Ms.push_back(MCMC(n, β, *As[i]));
     }
-
-    for (
   }
 
   bool step(unsigned i, unsigned j, bool dry = false) {
@@ -90,6 +63,107 @@ 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);
+
+    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, ε);
+      }
+
+      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];
+      }
+    }
+
+    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;
+
+    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(); i++) {
+      T0[i] = 1 / Ms[Ms.size() - i - 1].β;
+    }
+    T1.front() = 1 / Ms.back().β;
+    T1.back() = 1 / Ms.front().β;
+
+    for (unsigned i = 0; i < η.size(); i++) {
+      C += η[i] * (T0[i + 1] - T0[i]);
+    }
+
+    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++;
+      }
+      x = xnew;
+    }
+
+    for (unsigned i = 0; i < T1.size(); i++) {
+      Ms[Ms.size() - i - 1].β = 1 / T1[i];
+    }
+
+    for (double ff : f) {
+      outfile << ff << " ";
+    }
+    outfile << "\n";
+    outfile.close();
+
+  }
+
   void run(unsigned n, unsigned m, bool dry = false) {
     for (unsigned i = 0; i < n; i++) {
 #pragma omp parallel for
diff --git a/hadamard_test.cpp b/hadamard_test.cpp
new file mode 100644
index 0000000..c1dd8cd
--- /dev/null
+++ b/hadamard_test.cpp
@@ -0,0 +1,21 @@
+#include "hadamard_pt.hpp"
+
+int main() {
+  double β0 = 1.0;
+  double β1 = 10.0;
+  unsigned N = 8;
+  unsigned n = 4;
+
+  ParallelMeasurement m1;
+  std::vector<Measurement> m2(N);
+  std::vector<Measurement*> m3;
+  for (Measurement& m : m2) {
+    m3.push_back(&m);
+  }
+
+  PT a(β0, β1, N, n, m1, m3);
+
+  while (true) {
+    a.tune(1000, 10, 0.1);
+  }
+}