1 files changed, 9 insertions, 12 deletions

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;
     }
   }