partial implementation of dynamic pt

3 files changed, 47 insertions, 22 deletions

diff --git a/hadamard.cpp b/hadamard.cpp
index 1d99a2e..1afb270 100644
--- a/hadamard.cpp
+++ b/hadamard.cpp
@@ -92,9 +92,6 @@ int main(int argc, char* argv[]) {
     case 'M':
       M = (unsigned)atof(optarg);
       break;
-    case 'm':
-      m = (unsigned)atof(optarg);
-      break;
     default:
       exit(1);
     }
@@ -135,11 +132,10 @@ int main(int argc, char* argv[]) {
     sim.E = sim.M.energy();
   }
 
-  std::cout << "Beginning " << n_tuning << " tuning steps for " << n << ".\n";
-  p.tune(n_tuning, ε);
-  std::cout << "Finished tuning, beginning " << N << " dry tempering updates of " << M
+  std::cout << "Beginning simulation of " << n << ".\n";
+  std::cout << "Beginning " << n_tuning << " tuning tempering updates of " << M
             << " sweeps each.\n";
-  p.run(m, M, true);
+  p.tune(n_tuning, M, ε);
   std::cout << "Finished tuning, beginning " << N << " measurement tempering updates of " << M
             << " sweeps each.\n";
   p.run(N, M);
diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp
index 9f0d8c0..c9fa2da 100644
--- a/hadamard_mcmc.hpp
+++ b/hadamard_mcmc.hpp
@@ -140,6 +140,12 @@ public:
   virtual void after_sweep(double, double, const Orthogonal&){};
 };
 
+typedef enum {
+  none,
+  up,
+  down
+} color;
+
 class MCMC {
 private:
   randutils::mt19937_rng rng;
@@ -150,10 +156,14 @@ public:
   const double β;
   double E;
   Orthogonal M;
+  unsigned tag;
+  color c;
 
-  MCMC(unsigned n, double β0, Measurement& A) : A(A), M(n), β(β0) {
+  MCMC(unsigned n, double β0, Measurement& A, unsigned tag = 0) : 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 a3d4657..fba97b5 100644
--- a/hadamard_pt.hpp
+++ b/hadamard_pt.hpp
@@ -6,6 +6,8 @@
 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 {
@@ -29,23 +31,41 @@ public:
   ParallelMeasurement& B;
   std::vector<Measurement*>& As;
 
-  PT(std::list<range> ranges, unsigned n, ParallelMeasurement& B, std::vector<Measurement*>& As)
+  PT(double β₀, double β₁, unsigned n, ParallelMeasurement& B, std::vector<Measurement*>& As)
       : B(B), As(As) {
-    unsigned count = 0;
-    for (range r : ranges) {
-      for (unsigned i = 1; i <= r.N; i++) {
-        double β = r.β0 + i * (r.β1 - r.β0) / r.N;
-        Ms.push_back(MCMC(n, β, *As[count]));
-        count++;
-      }
+    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, double ε) {
+  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 i = 0; i < Ms.size(); i++) {
-      Ms[i].tune(N, ε);
+      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]++;
+        }
+      }
     }
+
+    for (
   }
 
   bool step(unsigned i, unsigned j, bool dry = false) {
@@ -64,10 +84,9 @@ public:
   }
 
   void sweep(bool dry = false) {
+
     for (unsigned i = 0; i < Ms.size() - 1; i++) {
-      for (unsigned j = i + 1; j < Ms.size(); j++) {
-        this->step(i, j, dry);
-      }
+      this->step(i, i + 1, dry);
     }
   }