From edc46c9a6e23da74bc80f41b79c928ce6edd06dd Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 20 Dec 2019 00:34:49 -0500
Subject: maybe fixed nonconstant problem

---
 hadamard_pt.hpp | 38 +++++++++++++++++---------------------
 1 file changed, 17 insertions(+), 21 deletions(-)

diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
index 74c3b73..3eee153 100644
--- a/hadamard_pt.hpp
+++ b/hadamard_pt.hpp
@@ -65,7 +65,6 @@ public:
   }
 
   void tune(unsigned n0, unsigned m, double ε, double ε2) {
-
     unsigned n = n0;
 
     while (true) {
@@ -73,9 +72,10 @@ public:
       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++) {
-        colors.front() = down;
-        colors.back() = up;
 #pragma omp parallel for
         for (unsigned j = 0; j < Ms.size(); j++) {
           Ms[j].tune(m, ε);
@@ -84,58 +84,53 @@ public:
         for (unsigned j = 0; j < Ms.size() - 1; j++) {
           if (this->step(j, j + 1, true)) {
             std::swap(colors[j], colors[j + 1]);
+          }
+        }
 
-            if (colors[j] == up) {
-              nu[j]++;
-            } else if (colors[j] == down) {
-              nd[j]++;
-            }
-            if (colors[j + 1] == up) {
-              nu[j + 1]++;
-            } else if (colors[j + 1] == down) {
-              nd[j + 1]++;
-            }
+        colors.front() = down;
+        colors.back() = up;
 
+        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());
 
-      double err = 0;
-      double δ = 0.5;
       std::vector<unsigned> f_keep;
       double f_last = 0;
       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)));
-        if (δerr < δ && f[i] >= f_last) {
+        if (f[i] >= f_last) {
           f_keep.push_back(i);
           f_last = f[i];
         }
       }
 
-      for (unsigned i = 0; i < f_keep.size() - 1; i++) {
+      for (signed 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]].β);
         }
       }
 
       std::vector<double> η(Ms.size() - 1);
-      double C = 0;
 
       for (unsigned i = 0; i < Ms.size() - 1; i++) {
         η[η.size() - i - 1] = sqrt((f[i + 1] - f[i])) / (1 / Ms[i].β - 1 / Ms[i+1].β);
       }
 
-      std::vector<double> T1(Ms.size() - 2);
-
+      double C = 0;
       for (unsigned i = 0; i < η.size(); i++) {
         C += η[i] * (T(i + 1) - T(i));
       }
 
       double x= 0;
       unsigned j = 1;
+      std::vector<double> T1(Ms.size() - 2);
 
       for (unsigned i = 0; i < η.size(); i++) {
         double xnew = x + η[i] * (T(i + 1) - T(i)) / C;
@@ -146,6 +141,7 @@ public:
         x = xnew;
       }
 
+      double err = 0;
       for (unsigned i = 0; i < T1.size(); i++) {
         err += fabs(T1[i] - 1/ Ms[Ms.size() - i - 2].β);
         Ms[Ms.size() - i - 2].β = 1 / T1[i];
-- 
cgit v1.2.3-54-g00ecf