1 files changed, 42 insertions, 11 deletions

diff --git a/hadamard.cpp b/hadamard.cpp
index 4630ca0..7276039 100644
--- a/hadamard.cpp
+++ b/hadamard.cpp
@@ -2,6 +2,7 @@
 #include "hadamard_pt.hpp"
 #include "matrices.hpp"
 
+#include <execution>
 #include <chrono>
 #include <fstream>
 #include <iostream>
@@ -63,9 +64,12 @@ int main(int argc, char* argv[]) {
   unsigned M = 10;
   unsigned m = 1e4;
 
+  bool loadBetasFromFile = false;
+  std::string betaFilename;
+
   int opt;
 
-  while ((opt = getopt(argc, argv, "k:b:c:n:t:N:M:e:m:f:")) != -1) {
+  while ((opt = getopt(argc, argv, "k:b:c:n:t:N:M:e:m:f:F:")) != -1) {
     switch (opt) {
     case 'k':
       k = atoi(optarg);
@@ -98,6 +102,10 @@ int main(int argc, char* argv[]) {
     case 'M':
       M = (unsigned)atof(optarg);
       break;
+    case 'F':
+      loadBetasFromFile = true;
+      betaFilename = optarg;
+      break;
     default:
       exit(1);
     }
@@ -121,17 +129,40 @@ int main(int argc, char* argv[]) {
   std::vector<double> f;
 
   std::cout << "Beginning simulation of " << n << ".\n";
-  std::cout << "Beginning " << n_tuning << " tuning tempering updates of " << M
-            << " sweeps each.\n";
-  f = p.tune(n_tuning, M, ε, ε2);
-  std::cout << "Finished tuning, beginning " << m << " measurement tempering updates of " << M
-            << " sweeps each.\n";
-  std::cout << "βs: ";
-  for (const MCMC& M : p.Ms) {
-    std::cout << M.β << " ";
+  if (loadBetasFromFile) {
+    std::cout << "Loading βs from file.\n";
+    std::ifstream betaFile(betaFilename);
+    for (unsigned i = 0; i < N; i++) {
+      double βtmp;
+      betaFile >> βtmp;
+      p.Ms[i].β = βtmp;
+    }
+    std::cout << "βs: ";
+    for (const MCMC& M : p.Ms) {
+      std::cout << M.β << " ";
+    }
+    std::cout << std::endl;
+    std::cout << "Beginning " << n_tuning << " tuning tempering updates.\n";
+    Rng rng;
+    for (unsigned i = 0; i < n_tuning; i++) {
+      std::for_each(std::execution::par_unseq, p.Ms.begin(), p.Ms.end(),
+                    [m, ε](MCMC& M) { M.tune(m, ε); });
+      for (unsigned j = 0; j < p.Ms.size(); j++) {
+        unsigned k = rng.uniform((unsigned)0, (unsigned)(p.Ms.size() - 2));
+        p.step(k, k + 1, true);
+      }
+    }
+  } else {
+    std::cout << "Beginning " << n_tuning << " tuning tempering updates of " << M
+              << " sweeps each.\n";
+    f = p.tune(n_tuning, M, ε, ε2);
+    std::cout << "βs: ";
+    for (const MCMC& M : p.Ms) {
+      std::cout << M.β << " ";
+    }
+    std::cout << std::endl;
   }
-  std::cout << std::endl;
-  std::cout << "Running " << m << " PT swaps of " << M << " sweeps each.";
+  std::cout << "Running " << m << " PT swaps of " << M << " sweeps each.\n";
   p.run(m, M); 
   std::cout << "Finished " << n << ".\n";