working pt with tuning

3 files changed, 44 insertions, 141 deletions

diff --git a/hadamard.cpp b/hadamard.cpp
index 1afb270..8eba184 100644
--- a/hadamard.cpp
+++ b/hadamard.cpp
@@ -3,6 +3,7 @@
 
 #include <fstream>
 #include <iostream>
+#include <chrono>
 
 class MeasureEnergy : public Measurement {
 public:
@@ -33,20 +34,17 @@ public:
 
 class MeasureTransitionRates : public ParallelMeasurement {
 public:
-  std::vector<std::vector<unsigned>> nAccepted;
+  std::vector<unsigned> nAccepted;
   unsigned total_steps;
 
-  MeasureTransitionRates(unsigned n) : nAccepted(n - 1) {
+  MeasureTransitionRates(unsigned n) : nAccepted(n - 1, 0) {
     total_steps = 0;
-    for (unsigned i = 0; i < n - 1; i++) {
-      nAccepted[i].resize(i + 1);
-    }
   }
 
-  void after_step(bool accepted, unsigned i, unsigned j, double, double, const MCMC&,
+  void after_step(bool accepted, unsigned i, double, double, const MCMC&,
                   const MCMC&) override {
     if (accepted)
-      nAccepted[j - 1][i]++;
+      nAccepted[i]++;
   }
 
   void after_sweep(const std::vector<MCMC>&) override { total_steps++; }
@@ -54,14 +52,13 @@ public:
 
 int main(int argc, char* argv[]) {
   unsigned n_tuning = 1e2;
-  std::list<double> β0s;
-  std::list<double> β1s;
-  std::list<unsigned> Ns;
+  double β₀ = 0.1;
+  double β₁ = 10;
+  unsigned N = 16;
   unsigned k = 2;
   double ε = 0.01;
 
   unsigned M = 10;
-  unsigned N = 1e4;
   unsigned m = 1e4;
 
   int opt;
@@ -72,16 +69,16 @@ int main(int argc, char* argv[]) {
       k = atoi(optarg);
       break;
     case 'b':
-      β0s.push_back(atof(optarg));
+      β₀ = atof(optarg);
       break;
     case 'c':
-      β1s.push_back(atof(optarg));
+      β₁ = atof(optarg);
       break;
     case 'e':
       ε = atof(optarg);
       break;
     case 'n':
-      Ns.push_back((unsigned)atof(optarg));
+      m = (unsigned)atof(optarg);
       break;
     case 't':
       n_tuning = (unsigned)atof(optarg);
@@ -99,33 +96,13 @@ int main(int argc, char* argv[]) {
 
   unsigned n = pow(2, k);
 
-  std::list<range> rs;
-  unsigned num = 0;
-
-  if ((β0s.size() != β1s.size()) || (β0s.size() != Ns.size())) {
-    std::cout << "You need the same number of ranges!\n";
-    exit(0);
-  } else {
-    auto it0 = β0s.begin();
-    auto it1 = β1s.begin();
-    auto itN = Ns.begin();
-
-    while (it0 != β0s.end()) {
-      rs.push_back({*it0, *it1, *itN});
-      num += *itN;
-      it0++;
-      it1++;
-      itN++;
-    }
-  }
-
-  std::vector<Measurement*> As(num);
+  std::vector<Measurement*> As(N);
   for (Measurement*& A : As) {
     A = new MeasureEnergy();
   }
-  MeasureTransitionRates B(num);
+  MeasureTransitionRates B(N);
 
-  PT p(rs, n, B, As);
+  PT p(β₀, β₁, N, n, B, As);
 
   for (MCMC& sim : p.Ms) {
     sim.M = walsh(k);
@@ -135,115 +112,42 @@ int main(int argc, char* argv[]) {
   std::cout << "Beginning simulation of " << n << ".\n";
   std::cout << "Beginning " << n_tuning << " tuning tempering updates of " << M
             << " sweeps each.\n";
-  p.tune(n_tuning, M, ε);
-  std::cout << "Finished tuning, beginning " << N << " measurement tempering updates of " << M
+  std::vector<double> f = p.tune(n_tuning, M, ε, 0.05);
+  std::cout << "Finished tuning, beginning " << m << " measurement tempering updates of " << M
             << " sweeps each.\n";
-  p.run(N, M);
+  p.run(m, M);
   std::cout << "Finished " << n << ".\n";
 
-  std::string rs_string = "";
-
-  for (range r : rs) {
-    rs_string +=
-        "_" + std::to_string(r.β0) + "_" + std::to_string(r.β1) + "_" + std::to_string(r.N);
-  }
+  auto tag = std::chrono::high_resolution_clock::now();
 
-  std::string filename = "probs_" + std::to_string(n) + rs_string + ".dat";
-  std::ifstream file(filename);
+  std::string filename = "hmm_" + std::to_string(n) + "_" + std::to_string(β₀) + "_" + std::to_string(β₁) + "_" + std::to_string(N) + "_" + std::to_string(tag.time_since_epoch().count()) + ".dat";
 
-  unsigned N_old = 0;
-  std::vector<std::vector<unsigned long>> data_old(B.nAccepted.size());
+  std::ofstream file(filename);
 
-  for (unsigned i = 0; i < B.nAccepted.size(); i++) {
-    data_old[i].resize(B.nAccepted[i].size());
+  for (const MCMC& M : p.Ms) {
+    file << M.β << " ";
   }
-
-  if (file.is_open()) {
-    file >> N_old;
-
-    for (unsigned i = 0; i < B.nAccepted.size(); i++) {
-      for (unsigned j = 0; j < B.nAccepted[i].size(); j++) {
-        double num;
-        file >> num;
-        data_old[i][j] = num;
-      }
-    }
-
-    file.close();
-  }
-
-  std::ofstream file_out(filename);
-
-  file_out << N_old + B.total_steps << "\n";
+  file << std::endl;
 
   for (unsigned i = 0; i < B.nAccepted.size(); i++) {
-    for (unsigned j = 0; j < B.nAccepted[i].size(); j++) {
-      file_out << std::fixed << data_old[i][j] + B.nAccepted[i][j] << " ";
-    }
-    file_out << "\n";
+    file << std::fixed << B.nAccepted[i] / (double)B.total_steps << " ";
   }
-
-  file_out.close();
-
-  std::string efilename = "energies_" + std::to_string(n) + rs_string + ".dat";
-  std::ifstream efile(efilename);
-
-  unsigned Ne_old = 0;
-  std::vector<double> edata_old(p.As.size());
-
-  if (efile.is_open()) {
-    efile >> Ne_old;
-
-    for (unsigned i = 0; i < p.As.size(); i++) {
-      double num;
-      efile >> num;
-      edata_old[i] = num;
-    }
-
-    efile.close();
-  }
-
-  std::ofstream efile_out(efilename);
-
-  efile_out << Ne_old + ((MeasureEnergy*)As[0])->N << "\n";
-
-  for (unsigned i = 0; i < As.size(); i++) {
-    efile_out << std::fixed << edata_old[i] + ((MeasureEnergy*)As[i])->totalE << " ";
-  }
-
-  efile_out.close();
-
-  std::string ρfilename = "rhos_" + std::to_string(n) + rs_string + ".dat";
-  std::ifstream ρfile(ρfilename);
-
-  std::vector<std::vector<unsigned>> ρdata_old(As.size());
+  file << std::endl;
 
   for (unsigned i = 0; i < As.size(); i++) {
-    ρdata_old[i].resize(((MeasureEnergy*)As[0])->ρ_dist.size());
-  }
-
-  if (ρfile.is_open()) {
-    for (unsigned i = 0; i < As.size(); i++) {
-      for (unsigned j = 0; j < ((MeasureEnergy*)As[0])->ρ_dist.size(); j++) {
-        unsigned num;
-        ρfile >> num;
-        ρdata_old[i][j] = num;
-      }
-    }
-
-    ρfile.close();
+    file << std::fixed << ((MeasureEnergy*)As[i])->totalE / ((MeasureEnergy*)As[i])->N << " ";
   }
 
-  std::ofstream ρfile_out(ρfilename);
+  file << std::endl;
 
   for (unsigned i = 0; i < As.size(); i++) {
-    for (unsigned j = 0; j < ((MeasureEnergy*)As[0])->ρ_dist.size(); j++) {
-      ρfile_out << std::fixed << ρdata_old[i][j] + ((MeasureEnergy*)As[i])->ρ_dist[j] << " ";
+    for (unsigned j = 0; j < ((MeasureEnergy*)As[i])->ρ_dist.size(); j++) {
+      file << std::fixed << ((MeasureEnergy*)As[i])->ρ_dist[j] << " ";
     }
-    ρfile_out << "\n";
+    file << std::endl;
   }
 
-  ρfile_out.close();
+  file.close();
 
   return 0;
 }
diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
index f1cdf17..5500f57 100644
--- a/hadamard_pt.hpp
+++ b/hadamard_pt.hpp
@@ -2,8 +2,6 @@
 #pragma once
 #include "hadamard_mcmc.hpp"
 #include <list>
-#include <fstream>
-#include <iostream>
 
 void swap(MCMC& s1, MCMC& s2) {
   std::swap(s1.M, s2.M);
@@ -12,7 +10,7 @@ void swap(MCMC& s1, MCMC& s2) {
 
 class ParallelMeasurement {
 public:
-  virtual void after_step(bool, unsigned, unsigned, double, double, const MCMC&, const MCMC&){};
+  virtual void after_step(bool, unsigned, double, double, const MCMC&, const MCMC&){};
   virtual void after_sweep(const std::vector<MCMC>&){};
 };
 
@@ -54,7 +52,7 @@ public:
       swap(Ms[i], Ms[j]);
 
     if (!dry) {
-      B.after_step(accepted, i, j, Δβ, ΔE, Ms[i], Ms[j]);
+      B.after_step(accepted, i, Δβ, ΔE, Ms[i], Ms[j]);
     }
     return accepted;
   }
@@ -66,7 +64,7 @@ public:
     }
   }
 
-  void tune(unsigned n0, unsigned m, double ε, double ε2) {
+  std::vector<double> tune(unsigned n0, unsigned m, double ε, double ε2) {
 
     unsigned n = n0;
 
@@ -112,18 +110,12 @@ public:
       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];
         }
       }
 
-      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]].β);
@@ -134,7 +126,7 @@ public:
       double C = 0;
 
       for (unsigned i = 0; i < Ms.size() - 1; i++) {
-        η[η.size() - i - 1] = sqrt((f[i + 1] - f[i])) / (T(i) - T(i + 1));
+        η[η.size() - i - 1] = sqrt((f[i + 1] - f[i])) / (1 / Ms[i].β - 1 / Ms[i+1].β);
       }
 
       std::vector<double> T1(Ms.size() - 2);
@@ -156,9 +148,14 @@ public:
       }
 
       for (unsigned i = 0; i < T1.size(); i++) {
+        err += fabs(T1[i] - 1/ Ms[Ms.size() - i - 2].β) / T1[i];
         Ms[Ms.size() - i - 2].β = 1 / T1[i];
       }
 
+      if (err / T1.size() < ε2) {
+        return f;
+      }
+
       n *= 2;
     }
   }
diff --git a/hadamard_test.cpp b/hadamard_test.cpp
index e505d79..a27b117 100644
--- a/hadamard_test.cpp
+++ b/hadamard_test.cpp
@@ -1,9 +1,10 @@
 #include "hadamard_pt.hpp"
+#include <iostream>
 
 int main() {
   double β0 = 1.0;
   double β1 = 10.0;
-  unsigned N = 4;
+  unsigned N = 16;
   unsigned n = 4;
 
   ParallelMeasurement m1;
@@ -18,6 +19,7 @@ int main() {
   a.tune(1000, 10, 0.1, 0.1);
 
   for (MCMC& M : a.Ms) {
-    std::cout << M.β;
+    std::cout << M.β << " ";
   }
+  std::cout << "\n";
 }