From 2cc44d8981e42bfe663e8e7dee8dbb7914a9482f Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 20 Dec 2019 17:37:01 -0500
Subject: fixed it

---
 hadamard.cpp      |  7 ++++++-
 hadamard_mcmc.hpp |  2 +-
 hadamard_pt.hpp   | 35 +++++++++++++++++++----------------
 3 files changed, 26 insertions(+), 18 deletions(-)

diff --git a/hadamard.cpp b/hadamard.cpp
index 5154677..2f45dcc 100644
--- a/hadamard.cpp
+++ b/hadamard.cpp
@@ -116,7 +116,7 @@ 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, ε, ε2);
+  std::vector<double> f = p.tune(n_tuning, M, ε, ε2);
   std::cout << "Finished tuning, beginning " << m << " measurement tempering updates of " << M
             << " sweeps each.\n";
   p.run(m, M);
@@ -133,6 +133,11 @@ int main(int argc, char* argv[]) {
   }
   file << std::endl;
 
+  for (double ff : f) {
+    file << ff << " ";
+  }
+  file << std::endl;
+
   for (unsigned i = 0; i < B.nAccepted.size(); i++) {
     file << std::fixed << B.nAccepted[i] / (double)B.total_steps << " ";
   }
diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp
index 145454b..eab9916 100644
--- a/hadamard_mcmc.hpp
+++ b/hadamard_mcmc.hpp
@@ -148,11 +148,11 @@ typedef enum {
 
 class MCMC {
 private:
-  randutils::mt19937_rng rng;
   Measurement& A;
   double θ0;
 
 public:
+  randutils::mt19937_rng rng;
   double β;
   double E;
   Orthogonal M;
diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
index 4cde211..de6590d 100644
--- a/hadamard_pt.hpp
+++ b/hadamard_pt.hpp
@@ -2,7 +2,6 @@
 #pragma once
 #include "hadamard_mcmc.hpp"
 #include <list>
-#include <iostream>
 
 void swap(MCMC& s1, MCMC& s2) {
   std::swap(s1.M, s2.M);
@@ -65,15 +64,16 @@ public:
     }
   }
 
-  void tune(unsigned n0, unsigned m, double ε, double ε2) {
+  std::vector<double> tune(unsigned n0, unsigned m, double ε, double ε2) {
     unsigned n = n0;
 
     while (true) {
       std::vector<color> colors(Ms.size(), none);
+      colors.front() = down;
+      colors.back() = up;
+
       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
@@ -81,23 +81,25 @@ public:
           Ms[j].tune(m, ε);
         }
 
-        for (unsigned j = 0; j < Ms.size() - 1; j++) {
+        for (unsigned k = 0; k < m * Ms.size() - 1; k++) {
+          unsigned j = Ms[0].rng.uniform((unsigned)0, (unsigned)(Ms.size() - 2));
+
           if (this->step(j, j + 1, true)) {
             std::swap(colors[j], colors[j + 1]);
             colors.front() = down;
             colors.back() = up;
           }
-          if (i > n / 2) {
-            for (unsigned j = 0; j < Ms.size(); j++) {
-              if (colors[j] == up) {
-                nu[j]++;
-              } else if (colors[j] == down) {
-                nd[j]++;
-              }
-            }
+        }
+
+        for (unsigned j = 0; j < Ms.size(); j++) {
+          if (colors[j] == up) {
+            nu[j]++;
+          } else if (colors[j] == down) {
+            nd[j]++;
           }
         }
 
+
       }
 
       std::vector<double> f(Ms.size());
@@ -129,10 +131,11 @@ public:
         C += η[i] * (T(i + 1) - T(i));
       }
 
-      double x= 0;
-      unsigned j = 1;
       std::vector<double> T1(Ms.size() - 2);
 
+      double x = 0;
+      unsigned j = 1;
+
       for (unsigned i = 0; i < η.size(); i++) {
         double xnew = x + η[i] * (T(i + 1) - T(i)) / C;
         while (j < xnew * η.size() && j < η.size()) {
@@ -149,7 +152,7 @@ public:
       }
 
       if (err / T1.size() * Ms.size() / (1/Ms.front().β - 1/Ms.back().β) < ε2) {
-        break;
+        return f;
       }
 
       n *= 2;
-- 
cgit v1.2.3-54-g00ecf