added a stage for burnout steps to address superheating

3 files changed, 27 insertions, 14 deletions

diff --git a/hadamard.cpp b/hadamard.cpp
index 74a8f3f..1d99a2e 100644
--- a/hadamard.cpp
+++ b/hadamard.cpp
@@ -62,10 +62,11 @@ int main(int argc, char* argv[]) {
 
   unsigned M = 10;
   unsigned N = 1e4;
+  unsigned m = 1e4;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "k:b:c:n:t:N:M:e:")) != -1) {
+  while ((opt = getopt(argc, argv, "k:b:c:n:t:N:M:e:m:")) != -1) {
     switch (opt) {
     case 'k':
       k = atoi(optarg);
@@ -91,6 +92,9 @@ int main(int argc, char* argv[]) {
     case 'M':
       M = (unsigned)atof(optarg);
       break;
+    case 'm':
+      m = (unsigned)atof(optarg);
+      break;
     default:
       exit(1);
     }
@@ -133,7 +137,10 @@ int main(int argc, char* argv[]) {
 
   std::cout << "Beginning " << n_tuning << " tuning steps for " << n << ".\n";
   p.tune(n_tuning, ε);
-  std::cout << "Finished tuning, beginning " << N << " tempering updates of " << M
+  std::cout << "Finished tuning, beginning " << N << " dry tempering updates of " << M
+            << " sweeps each.\n";
+  p.run(m, M, true);
+  std::cout << "Finished tuning, beginning " << N << " measurement tempering updates of " << M
             << " sweeps each.\n";
   p.run(N, M);
   std::cout << "Finished " << n << ".\n";
diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp
index 5ef062f..9f0d8c0 100644
--- a/hadamard_mcmc.hpp
+++ b/hadamard_mcmc.hpp
@@ -156,7 +156,7 @@ public:
     E = M.energy();
   }
 
-  bool step(Givens& g) {
+  bool step(Givens& g, bool dry = false) {
     double ΔE = g.tryRotation();
 
     bool accepted = ΔE < 0 || exp(-β * ΔE) > rng.uniform((double)0.0, 1.0);
@@ -166,11 +166,13 @@ public:
       g.acceptRotation();
     }
 
-    A.after_step(accepted, g, θ0, E, ΔE, M);
+    if (!dry) {
+      A.after_step(accepted, g, θ0, E, ΔE, M);
+    }
     return accepted;
   }
 
-  double sweep() {
+  double sweep(bool dry = false) {
     unsigned total = 0;
     unsigned accepted = 0;
 
@@ -179,9 +181,9 @@ public:
         Givens g1(M, false, i, j, θ0, rng);
         Givens g2(M, true, i, j, θ0, rng);
 
-        if (this->step(g1))
+        if (this->step(g1, dry))
           accepted++;
-        if (this->step(g2))
+        if (this->step(g2, dry))
           accepted++;
 
         total += 2;
@@ -193,7 +195,7 @@ public:
 
   void tune(unsigned N, double ε) {
     for (unsigned i = 0; i < N; i++) {
-      double ratio_accepted = this->sweep();
+      double ratio_accepted = this->sweep(true);
       if (ratio_accepted > 0.5) {
         θ0 *= 1 + ε;
       } else {
@@ -202,10 +204,12 @@ public:
     }
   }
 
-  void run(unsigned N) {
+  void run(unsigned N, bool dry = false) {
     for (unsigned i = 0; i < N; i++) {
-      this->sweep();
-      A.after_sweep(θ0, E, M);
+      this->sweep(dry);
+      if (!dry) {
+        A.after_sweep(θ0, E, M);
+      }
     }
   }
 };
diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
index dad635f..bddd8b7 100644
--- a/hadamard_pt.hpp
+++ b/hadamard_pt.hpp
@@ -69,14 +69,16 @@ public:
     }
   }
 
-  void run(unsigned n, unsigned m) {
+  void run(unsigned n, unsigned m, bool dry = false) {
     for (unsigned i = 0; i < n; i++) {
 #pragma omp parallel for
       for (unsigned j = 0; j < Ms.size(); j++) {
-        Ms[j].run(m);
+        Ms[j].run(m, dry);
       }
       this->sweep();
-      B.after_sweep(this->Ms);
+      if (!dry) {
+        B.after_sweep(this->Ms);
+      }
     }
   }
 };