summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--hadamard_pt.hpp29
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;
}
}