Commit 670f5e95 authored by Stefan Gerlach's avatar Stefan Gerlach

[fit] Check if chi is zero instead of chi^2 to avoid precision problems

parent e37cefcd
......@@ -1824,16 +1824,16 @@ void XYFitCurvePrivate::recalculate() {
DEBUG(Q_FUNC_INFO << ", run fdfsolver_iterate");
status = gsl_multifit_fdfsolver_iterate(s);
DEBUG(Q_FUNC_INFO << ", fdfsolver_iterate DONE");
double chi2 = gsl_pow_2(gsl_blas_dnrm2(s->f));
writeSolverState(s, chi2);
double chi = gsl_blas_dnrm2(s->f);
writeSolverState(s, chi);
if (status) {
DEBUG(Q_FUNC_INFO << ", iter " << iter << ", status = " << gsl_strerror(status));
if (status == GSL_ETOLX) // change in the position vector falls below machine precision: no progress
status = GSL_SUCCESS;
break;
}
if (qFuzzyIsNull(chi2)) {
DEBUG(Q_FUNC_INFO << ", chi^2 is zero! Finishing.")
if (qFuzzyIsNull(chi)) {
DEBUG(Q_FUNC_INFO << ", chi is zero! Finishing.")
status = GSL_SUCCESS;
} else {
status = gsl_multifit_test_delta(s->dx, s->x, delta, delta);
......@@ -1846,11 +1846,11 @@ void XYFitCurvePrivate::recalculate() {
DEBUG(Q_FUNC_INFO << ", Rerun fit with x errors");
unsigned int iter2 = 0;
double chi2 = 0, chi2Old = 0;
double chi = 0, chiOld = 0;
double *fun = new double[n];
do {
iter2++;
chi2Old = chi2;
chiOld = chi;
//printf("iter2 = %d\n", iter2);
// calculate function from residuals
......@@ -1941,8 +1941,8 @@ void XYFitCurvePrivate::recalculate() {
status = gsl_multifit_test_delta(s->dx, s->x, delta, delta);
} while (status == GSL_CONTINUE && iter < maxIters);
chi2 = gsl_blas_dnrm2(s->f);
} while (iter2 < maxIters && fabs(chi2 - chi2Old) > fitData.eps);
chi = gsl_blas_dnrm2(s->f);
} while (iter2 < maxIters && fabs(chi - chiOld) > fitData.eps);
delete[] fun;
}
......@@ -2152,7 +2152,7 @@ void XYFitCurvePrivate::evaluate(bool preview) {
/*!
* writes out the current state of the solver \c s
*/
void XYFitCurvePrivate::writeSolverState(gsl_multifit_fdfsolver* s, double chi2) {
void XYFitCurvePrivate::writeSolverState(gsl_multifit_fdfsolver* s, double chi) {
QString state;
//current parameter values, semicolon separated
......@@ -2164,12 +2164,12 @@ void XYFitCurvePrivate::writeSolverState(gsl_multifit_fdfsolver* s, double chi2)
state += QString::number(nsl_fit_map_bound(x, min[i], max[i])) + '\t';
}
//current value of chi^2
if (isnan(chi2))
chi2 = gsl_pow_2(gsl_blas_dnrm2(s->f));
state += QString::number(chi2);
//current value of chi
if (isnan(chi))
chi = gsl_blas_dnrm2(s->f);
state += QString::number(chi*chi);
state += ';';
DEBUG(Q_FUNC_INFO << ", chi^2 = " << chi2);
DEBUG(Q_FUNC_INFO << ", chi^2 = " << chi*chi);
fitResult.solverOutput += state;
}
......
......@@ -63,7 +63,7 @@ public:
private:
void prepareResultColumns();
void writeSolverState(gsl_multifit_fdfsolver*, double chi2 = NAN);
void writeSolverState(gsl_multifit_fdfsolver*, double chi = NAN);
};
#endif
......@@ -539,6 +539,9 @@ void FitTest::testLinearWampler2() {
const int np = fitData.paramNames.size();
QCOMPARE(np, 6);
for (int i = 0; i < np; i++) {
DEBUG(std::setprecision(15) << fitResult.paramValues.at(i));
}
QCOMPARE(fitResult.paramValues.at(0), 1.0);
QCOMPARE(fitResult.paramValues.at(1), 0.1);
QCOMPARE(fitResult.paramValues.at(2), 0.01);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment