Testing new strategy for judging tuning convergence.

1 files changed, 24 insertions, 5 deletions

diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
index a57c670..90fa661 100644
--- a/hadamard_pt.hpp
+++ b/hadamard_pt.hpp
@@ -85,11 +85,13 @@ public:
       }
 
       std::vector<double> f(Ms.size());
+      std::vector<double> δf(Ms.size());
 
       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]);
+        δf[i] = f[i] * sqrt(1.0/nu[i] + 1.0/(nu[i]+nd[i]));
         if (f[i] >= f_last) {
           f_keep.push_back(i);
           f_last = f[i];
@@ -101,9 +103,21 @@ public:
           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]].β);
+          δf[j] = sqrt(pow(δf[f_keep[i]], 2) + pow(sqrt(pow(δf[f_keep[i + 1]], 2) + pow(δf[f_keep[i]], 2)) /
+                                    (Ms[f_keep[i + 1]].β - Ms[f_keep[i]].β) *
+                                    (Ms[j].β - Ms[f_keep[i]].β), 2));
         }
       }
 
+      // measure the difference between f and a straight line
+      double Δf² = 0;
+      double δΔf⁴ = 0;
+
+      for (unsigned j = 1; j < f.size(); j++) {
+        Δf² += pow(f[j] - (j + 1.0) / (Ms.size() + 1.0), 2);
+        δΔf⁴ += pow(2 * δf[j] * f[j], 2);
+      }
+
       std::vector<double> η(Ms.size() - 1);
 
       for (unsigned i = 0; i < Ms.size() - 1; i++) {
@@ -132,18 +146,23 @@ public:
       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];
+        Ms[Ms.size() - i - 2].β = 1.0 / T1[i];
       }
 
-      double relErr = err / T1.size() * Ms.size() / (1 / Ms.front().β - 1 / Ms.back().β);
+      double Δf = sqrt(Δf²);
+      double δΔf = sqrt(δΔf⁴) / 2 / Δf;
+      double relErr = Δf / Ms.size();
+      // double relErr = err / T1.size() * Ms.size() / (1 / Ms.front().β - 1 / Ms.back().β);
 
       if (relErr < ε2) {
         return f;
       } else {
-        std::cout << "Difference from linear in transit flow: " << relErr << ", continuing tuning.\n";
+        std::cout << "RMS difference from ideal transit flow is " << relErr << " ± " << δΔf / Ms.size() << ", continuing tuning.\n";
+        if (5 * δΔf > Δf || δΔf != δΔf) {
+          n *= 2;
+          std::cout << "Error in RMS difference too close to difference, increasing tuning to " << n << ".\n";
+        }
       }
-
-      n += n0;
     }
   }