diff --git a/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp b/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp index 60f42f354..19128eb93 100644 --- a/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp +++ b/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp @@ -1,2224 +1,2225 @@  /*************************************************************************** File : XYFitCurve.cpp Project : LabPlot Description : A xy-curve defined by a fit model -------------------------------------------------------------------- Copyright : (C) 2014-2017 Alexander Semke (alexander.semke@web.de) Copyright : (C) 2016-2017 Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ /*! \class XYFitCurve \brief A xy-curve defined by a fit model \ingroup worksheet */ #include "XYFitCurve.h" #include "XYFitCurvePrivate.h" #include "backend/core/AbstractColumn.h" #include "backend/core/column/Column.h" #include "backend/lib/commandtemplates.h" #include "backend/lib/macros.h" #include "backend/gsl/errors.h" #include "backend/gsl/ExpressionParser.h" extern "C" { #include #include #include #include #include #include #include "backend/gsl/parser.h" #include "backend/nsl/nsl_sf_stats.h" #include "backend/nsl/nsl_stats.h" } #include #include #include XYFitCurve::XYFitCurve(const QString& name) : XYAnalysisCurve(name, new XYFitCurvePrivate(this)) { } XYFitCurve::XYFitCurve(const QString& name, XYFitCurvePrivate* dd) : XYAnalysisCurve(name, dd) { } XYFitCurve::~XYFitCurve() { //no need to delete the d-pointer here - it inherits from QGraphicsItem //and is deleted during the cleanup in QGraphicsScene } void XYFitCurve::recalculate() { Q_D(XYFitCurve); d->recalculate(); } void XYFitCurve::initStartValues(const XYCurve* curve) { Q_D(XYFitCurve); XYFitCurve::FitData& fitData = d->fitData; initStartValues(fitData, curve); } void XYFitCurve::initStartValues(XYFitCurve::FitData& fitData, const XYCurve* curve) { DEBUG("XYFitCurve::initStartValues()"); if (!curve) { DEBUG(" no curve given"); return; } const Column* tmpXDataColumn = dynamic_cast(curve->xColumn()); const Column* tmpYDataColumn = dynamic_cast(curve->yColumn()); if (!tmpXDataColumn || !tmpYDataColumn) { DEBUG(" data columns not available"); return; } DEBUG(" x data rows = " << tmpXDataColumn->rowCount()); nsl_fit_model_category modelCategory = fitData.modelCategory; int modelType = fitData.modelType; int degree = fitData.degree; DEBUG(" fit model type = " << modelType << ", degree = " << degree); QVector& paramStartValues = fitData.paramStartValues; //QVector* xVector = static_cast* >(tmpXDataColumn->data()); //double xmean = gsl_stats_mean(xVector->constData(), 1, tmpXDataColumn->rowCount()); double xmin = tmpXDataColumn->minimum(); double xmax = tmpXDataColumn->maximum(); //double ymin = tmpYDataColumn->minimum(); //double ymax = tmpYDataColumn->maximum(); double xrange = xmax-xmin; //double yrange = ymax-ymin; DEBUG(" x min/max = " << xmin << ' ' << xmax); //DEBUG(" y min/max = " << ymin << ' ' << ymax); switch (modelCategory) { case nsl_fit_model_basic: switch (modelType) { case nsl_fit_model_polynomial: // not needed (works anyway) break; //TODO: handle basic models case nsl_fit_model_power: case nsl_fit_model_exponential: case nsl_fit_model_inverse_exponential: case nsl_fit_model_fourier: break; } break; case nsl_fit_model_peak: // use equidistant mu's and (xmax-xmin)/(10*degree) as sigma for (int d = 0; d < degree; d++) { paramStartValues[3*d+1] = xmin + (d+1.)*xrange/(degree+1.); paramStartValues[3*d] = xrange/(10.*degree); } break; case nsl_fit_model_growth: switch (modelType) { case nsl_fit_model_atan: case nsl_fit_model_tanh: case nsl_fit_model_algebraic_sigmoid: case nsl_fit_model_erf: case nsl_fit_model_gudermann: case nsl_fit_model_sigmoid: // use (xmax+xmin)/2 as mu and (xmax-xmin)/10 as sigma paramStartValues[1] = (xmax+xmin)/2.; paramStartValues[0] = xrange/10.; break; case nsl_fit_model_hill: paramStartValues[0] = xrange/10.; break; case nsl_fit_model_gompertz: //TODO break; } break; case nsl_fit_model_distribution: switch (modelType) { case nsl_sf_stats_gaussian: case nsl_sf_stats_laplace: case nsl_sf_stats_rayleigh_tail: case nsl_sf_stats_lognormal: case nsl_sf_stats_logistic: case nsl_sf_stats_sech: case nsl_sf_stats_cauchy_lorentz: case nsl_sf_stats_levy: // use (xmax+xmin)/2 as mu and (xmax-xmin)/10 as sigma paramStartValues[1] = (xmin+xmax)/2.; paramStartValues[0] = xrange/10.; break; //TODO: other types default: break; } break; case nsl_fit_model_custom: // not possible break; } } /*! * sets the parameter names for given model category, model type and degree in \c fitData for given action */ void XYFitCurve::initFitData(PlotDataDialog::AnalysisAction action) { if (!action) return; Q_D(XYFitCurve); XYFitCurve::FitData& fitData = d->fitData; if (action == PlotDataDialog::FitLinear) { //Linear fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_polynomial; fitData.degree = 1; } else if (action == PlotDataDialog::FitPower) { //Power fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_power; fitData.degree = 1; } else if (action == PlotDataDialog::FitExp1) { //Exponential (degree 1) fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_exponential; fitData.degree = 1; } else if (action == PlotDataDialog::FitExp2) { //Exponential (degree 2) fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_exponential; fitData.degree = 2; } else if (action == PlotDataDialog::FitInvExp) { //Inverse exponential fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_inverse_exponential; } else if (action == PlotDataDialog::FitGauss) { //Gauss fitData.modelCategory = nsl_fit_model_peak; fitData.modelType = nsl_fit_model_gaussian; fitData.degree = 1; } else if (action == PlotDataDialog::FitCauchyLorentz) { //Cauchy-Lorentz fitData.modelCategory = nsl_fit_model_peak; fitData.modelType = nsl_fit_model_lorentz; fitData.degree = 1; } else if (action == PlotDataDialog::FitTan) { //Arc tangent fitData.modelCategory = nsl_fit_model_growth; fitData.modelType = nsl_fit_model_atan; } else if (action == PlotDataDialog::FitTanh) { //Hyperbolic tangent fitData.modelCategory = nsl_fit_model_growth; fitData.modelType = nsl_fit_model_tanh; } else if (action == PlotDataDialog::FitErrFunc) { //Error function fitData.modelCategory = nsl_fit_model_growth; fitData.modelType = nsl_fit_model_erf; } else { //Custom fitData.modelCategory = nsl_fit_model_custom; fitData.modelType = 0; } XYFitCurve::initFitData(fitData); } /*! * sets the model expression and the parameter names for given model category, model type and degree in \c fitData */ void XYFitCurve::initFitData(XYFitCurve::FitData& fitData) { nsl_fit_model_category modelCategory = fitData.modelCategory; int modelType = fitData.modelType; QString& model = fitData.model; QStringList& paramNames = fitData.paramNames; QStringList& paramNamesUtf8 = fitData.paramNamesUtf8; int degree = fitData.degree; QVector& paramStartValues = fitData.paramStartValues; QVector& paramLowerLimits = fitData.paramLowerLimits; QVector& paramUpperLimits = fitData.paramUpperLimits; QVector& paramFixed = fitData.paramFixed; + DEBUG(" category = " << modelCategory); DEBUG("XYFitCurve::initFitData() for model category = " << nsl_fit_model_category_name[modelCategory] << ", model type = " << modelType << ", degree = " << degree); if (modelCategory != nsl_fit_model_custom) paramNames.clear(); paramNamesUtf8.clear(); // 10 indices used in multi degree models QStringList indices = { UTF8_QSTRING("₁"), UTF8_QSTRING("₂"), UTF8_QSTRING("₃"), UTF8_QSTRING("₄"), UTF8_QSTRING("₅"), UTF8_QSTRING("₆"), UTF8_QSTRING("₇"), UTF8_QSTRING("₈"), UTF8_QSTRING("₉"), UTF8_QSTRING("₁₀")}; switch (modelCategory) { case nsl_fit_model_basic: model = nsl_fit_model_basic_equation[fitData.modelType]; switch (modelType) { case nsl_fit_model_polynomial: paramNames << "c0" << "c1"; paramNamesUtf8 << UTF8_QSTRING("c₀") << UTF8_QSTRING("c₁"); if (degree == 2) { model += " + c2*x^2"; paramNames << "c2"; paramNamesUtf8 << UTF8_QSTRING("c₂"); } else if (degree > 2) { for (int i = 2; i <= degree; ++i) { QString numStr = QString::number(i); model += "+c" + numStr + "*x^" + numStr; paramNames << "c" + numStr; paramNamesUtf8 << "c" + indices[i-1]; } } break; case nsl_fit_model_power: if (degree == 1) { paramNames << "a" << "b"; } else { paramNames << "a" << "b" << "c"; model = "a + b*x^c"; } break; case nsl_fit_model_exponential: if (degree == 1) { paramNames << "a" << "b"; } else { for (int i = 1; i <= degree; i++) { QString numStr = QString::number(i); if (i == 1) model = "a1*exp(b1*x)"; else model += " + a" + numStr + "*exp(b" + numStr + "*x)"; paramNames << "a" + numStr << "b" + numStr; paramNamesUtf8 << "a" + indices[i-1] << "b" + indices[i-1]; } } break; case nsl_fit_model_inverse_exponential: paramNames << "a" << "b" << "c"; break; case nsl_fit_model_fourier: paramNames << "w" << "a0" << "a1" << "b1"; paramNamesUtf8 << UTF8_QSTRING("ω") << UTF8_QSTRING("a₀") << UTF8_QSTRING("a₁") << UTF8_QSTRING("b₁"); if (degree > 1) { for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); model += "+ (a" + numStr + "*cos(" + numStr + "*w*x) + b" + numStr + "*sin(" + numStr + "*w*x))"; paramNames << "a" + numStr << "b" + numStr; paramNamesUtf8 << "a" + indices[i-1] << "b" + indices[i-1]; } } break; } break; case nsl_fit_model_peak: model = nsl_fit_model_peak_equation[fitData.modelType]; switch (modelType) { case nsl_fit_model_gaussian: switch (degree) { case 1: paramNames << "s" << "mu" << "a"; paramNamesUtf8 << UTF8_QSTRING("σ") << UTF8_QSTRING("μ") << "A"; break; case 2: model = "1./sqrt(2*pi) * (a1/s1 * exp(-((x-mu1)/s1)^2/2) + a2/s2 * exp(-((x-mu2)/s2)^2/2))"; paramNames << "s1" << "mu1" << "a1" << "s2" << "mu2" << "a2"; paramNamesUtf8 << UTF8_QSTRING("σ₁") << UTF8_QSTRING("μ₁") << UTF8_QSTRING("A₁") << UTF8_QSTRING("σ₂") << UTF8_QSTRING("μ₂") << UTF8_QSTRING("A₂"); break; case 3: model = "1./sqrt(2*pi) * (a1/s1 * exp(-((x-mu1)/s1)^2/2) + a2/s2 * exp(-((x-mu2)/s2)^2/2) + a3/s3 * exp(-((x-mu3)/s3)^2/2))"; paramNames << "s1" << "mu1" << "a1" << "s2" << "mu2" << "a2" << "s3" << "mu3" << "a3"; paramNamesUtf8 << UTF8_QSTRING("σ₁") << UTF8_QSTRING("μ₁") << UTF8_QSTRING("A₁") << UTF8_QSTRING("σ₂") << UTF8_QSTRING("μ₂") << UTF8_QSTRING("A₂") << UTF8_QSTRING("σ₃") << UTF8_QSTRING("μ₃") << UTF8_QSTRING("A₃"); break; default: model = "1./sqrt(2*pi) * ("; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += "a" + numStr + "/s" + numStr + "* exp(-((x-mu" + numStr + ")/s" + numStr + ")^2/2)"; paramNames << "s" + numStr << "mu" + numStr << "a" + numStr; paramNamesUtf8 << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1] << "A" + indices[i-1]; } model += ")"; } break; case nsl_fit_model_lorentz: switch (degree) { case 1: paramNames << "g" << "mu" << "a"; paramNamesUtf8 << UTF8_QSTRING("γ") << UTF8_QSTRING("μ") << "A"; break; case 2: model = "1./pi * (a1 * g1/(g1^2+(x-mu1)^2) + a2 * g2/(g2^2+(x-mu2)^2))"; paramNames << "g1" << "mu1" << "a1" << "g2" << "mu2" << "a2"; paramNamesUtf8 << UTF8_QSTRING("γ₁") << UTF8_QSTRING("μ₁") << UTF8_QSTRING("A₁") << UTF8_QSTRING("γ₂") << UTF8_QSTRING("μ₂") << UTF8_QSTRING("A₂"); break; case 3: model = "1./pi * (a1 * g1/(g1^2+(x-mu1)^2) + a2 * g2/(g2^2+(x-mu2)^2) + a3 * g3/(g3^2+(x-mu3)^2))"; paramNames << "g1" << "mu1" << "a1" << "g2" << "mu2" << "a2" << "g3" << "mu3" << "a3"; paramNamesUtf8 << UTF8_QSTRING("γ₁") << UTF8_QSTRING("μ₁") << UTF8_QSTRING("A₁") << UTF8_QSTRING("γ₂") << UTF8_QSTRING("μ₂") << UTF8_QSTRING("A₂") << UTF8_QSTRING("γ₃") << UTF8_QSTRING("μ₃") << UTF8_QSTRING("A₃"); break; default: model = "1./pi * ("; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += "a" + numStr + " * g" + numStr + "/(g" + numStr + "^2+(x-mu" + numStr + ")^2)"; paramNames << "g" + numStr << "mu" + numStr << "a" + numStr; paramNamesUtf8 << UTF8_QSTRING("γ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1] << "A" + indices[i-1]; } model += ")"; } break; case nsl_fit_model_sech: switch (degree) { case 1: paramNames << "s" << "mu" << "a"; paramNamesUtf8 << UTF8_QSTRING("σ") << UTF8_QSTRING("μ") << "A"; break; case 2: model = "1/pi * (a1/s1 * sech((x-mu1)/s1) + a2/s2 * sech((x-mu2)/s2))"; paramNames << "s1" << "mu1" << "a1" << "s2" << "mu2" << "a2"; paramNamesUtf8 << UTF8_QSTRING("σ₁") << UTF8_QSTRING("μ₁") << UTF8_QSTRING("A₁") << UTF8_QSTRING("σ₂") << UTF8_QSTRING("μ₂") << UTF8_QSTRING("A₂"); break; case 3: model = "1/pi * (a1/s1 * sech((x-mu1)/s1) + a2/s2 * sech((x-mu2)/s2) + a3/s3 * sech((x-mu3)/s3))"; paramNames << "s1" << "mu1" << "a1" << "s2" << "mu2" << "a2" << "s3" << "mu3" << "a3"; paramNamesUtf8 << UTF8_QSTRING("σ₁") << UTF8_QSTRING("μ₁") << UTF8_QSTRING("A₁") << UTF8_QSTRING("σ₂") << UTF8_QSTRING("μ₂") << UTF8_QSTRING("A₂") << UTF8_QSTRING("σ₃") << UTF8_QSTRING("μ₃") << UTF8_QSTRING("A₃"); break; default: model = "1/pi * ("; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += "a" + numStr + "/s" + numStr + "* sech((x-mu" + numStr + ")/s" + numStr + ")"; paramNames << "s" + numStr << "mu" + numStr << "a" + numStr; paramNamesUtf8 << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1] << "A" + indices[i-1]; } model += ")"; } break; case nsl_fit_model_logistic: switch (degree) { case 1: paramNames << "s" << "mu" << "a"; paramNamesUtf8 << UTF8_QSTRING("σ") << UTF8_QSTRING("μ") << "A"; break; case 2: model = "1/4 * (a1/s1 * sech((x-mu1)/2/s1)**2 + a2/s2 * sech((x-mu2)/2/s2)**2)"; paramNames << "s1" << "mu1" << "a1" << "s2" << "mu2" << "a2"; paramNamesUtf8 << UTF8_QSTRING("σ₁") << UTF8_QSTRING("μ₁") << UTF8_QSTRING("A₁") << UTF8_QSTRING("σ₂") << UTF8_QSTRING("μ₂") << UTF8_QSTRING("A₂"); break; case 3: model = "1/4 * (a1/s1 * sech((x-mu1)/2/s1)**2 + a2/s2 * sech((x-mu2)/2/s2)**2 + a3/s3 * sech((x-mu3)/2/s3)**2)"; paramNames << "s1" << "mu1" << "a1" << "s2" << "mu2" << "a2" << "s3" << "mu3" << "a3"; paramNamesUtf8 << UTF8_QSTRING("σ₁") << UTF8_QSTRING("μ₁") << UTF8_QSTRING("A₁") << UTF8_QSTRING("σ₂") << UTF8_QSTRING("μ₂") << UTF8_QSTRING("A₂") << UTF8_QSTRING("σ₃") << UTF8_QSTRING("μ₃") << UTF8_QSTRING("A₃"); break; default: model = "1/4 * ("; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += "a" + numStr + "/s" + numStr + "* sech((x-mu" + numStr + ")/2/s" + numStr + ")**2"; paramNames << "s" + numStr << "mu" + numStr << "a" + numStr; paramNamesUtf8 << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1] << "A" + indices[i-1]; } model += ")"; } break; } break; case nsl_fit_model_growth: model = nsl_fit_model_growth_equation[fitData.modelType]; switch (modelType) { case nsl_fit_model_atan: case nsl_fit_model_tanh: case nsl_fit_model_algebraic_sigmoid: case nsl_fit_model_erf: case nsl_fit_model_gudermann: paramNames << "s" << "mu" << "a"; paramNamesUtf8 << UTF8_QSTRING("σ") << UTF8_QSTRING("μ") << "A"; break; case nsl_fit_model_sigmoid: paramNames << "k" << "mu" << "a"; paramNamesUtf8 << "k" << UTF8_QSTRING("μ") << "A"; break; case nsl_fit_model_hill: paramNames << "s" << "n" << "a"; paramNamesUtf8 << UTF8_QSTRING("σ") << "n" << "A"; break; case nsl_fit_model_gompertz: paramNames << "a" << "b" << "c"; break; } break; case nsl_fit_model_distribution: model = nsl_sf_stats_distribution_equation[fitData.modelType]; switch (modelType) { case nsl_sf_stats_gaussian: case nsl_sf_stats_laplace: case nsl_sf_stats_rayleigh_tail: case nsl_sf_stats_lognormal: case nsl_sf_stats_logistic: case nsl_sf_stats_sech: paramNames << "s" << "mu" << "a"; paramNamesUtf8 << UTF8_QSTRING("σ") << UTF8_QSTRING("μ") << "A"; break; case nsl_sf_stats_gaussian_tail: paramNames << "s" << "mu" << "A" << "a"; paramNamesUtf8 << UTF8_QSTRING("σ") << UTF8_QSTRING("μ") << "A" << "a"; break; case nsl_sf_stats_exponential: paramNames << "l" << "mu" << "a"; paramNamesUtf8 << UTF8_QSTRING("λ") << UTF8_QSTRING("μ") << "A"; break; case nsl_sf_stats_exponential_power: paramNames << "s" << "mu" << "b" << "a"; paramNamesUtf8 << UTF8_QSTRING("σ") << UTF8_QSTRING("μ") << "b" << "A"; break; case nsl_sf_stats_cauchy_lorentz: case nsl_sf_stats_levy: paramNames << "g" << "mu" << "a"; paramNamesUtf8 << UTF8_QSTRING("γ") << UTF8_QSTRING("μ") << "A"; break; case nsl_sf_stats_rayleigh: paramNames << "s" << "a"; paramNamesUtf8 << UTF8_QSTRING("σ") << "A"; break; case nsl_sf_stats_landau: paramNames << "a"; paramNamesUtf8 << "A"; break; case nsl_sf_stats_levy_alpha_stable: // unused distributions case nsl_sf_stats_levy_skew_alpha_stable: case nsl_sf_stats_bernoulli: break; case nsl_sf_stats_gamma: paramNames << "t" << "k" << "a"; paramNamesUtf8 << UTF8_QSTRING("θ") << "k" << "A"; break; case nsl_sf_stats_flat: paramNames << "a" << "b" << "A"; break; case nsl_sf_stats_chi_squared: paramNames << "n" << "a"; paramNamesUtf8 << "n" << "A"; break; case nsl_sf_stats_fdist: paramNames << "n1" << "n2" << "a"; paramNamesUtf8 << UTF8_QSTRING("ν₁") << UTF8_QSTRING("ν₂") << "A"; break; case nsl_sf_stats_tdist: paramNames << "n" << "a"; paramNamesUtf8 << UTF8_QSTRING("ν") << "A"; break; case nsl_sf_stats_beta: case nsl_sf_stats_pareto: paramNames << "a" << "b" << "A"; break; case nsl_sf_stats_weibull: paramNames << "k" << "l" << "mu" << "a"; paramNamesUtf8 << "k" << UTF8_QSTRING("λ") << UTF8_QSTRING("μ") << "A"; break; case nsl_sf_stats_gumbel1: paramNames << "s" << "b" << "mu" << "a"; paramNamesUtf8 << UTF8_QSTRING("σ") << UTF8_QSTRING("β") << UTF8_QSTRING("μ") << "A"; break; case nsl_sf_stats_gumbel2: paramNames << "a" << "b" << "mu" << "A"; paramNamesUtf8 << "a" << "b" << UTF8_QSTRING("μ") << "A"; break; case nsl_sf_stats_poisson: paramNames << "l" << "a"; paramNamesUtf8 << UTF8_QSTRING("λ") << "A"; break; case nsl_sf_stats_binomial: case nsl_sf_stats_negative_binomial: case nsl_sf_stats_pascal: paramNames << "p" << "n" << "a"; paramNamesUtf8 << "p" << "n" << "A"; break; case nsl_sf_stats_geometric: case nsl_sf_stats_logarithmic: paramNames << "p" << "a"; paramNamesUtf8 << "p" << "A"; break; case nsl_sf_stats_hypergeometric: paramNames << "n1" << "n2" << "t" << "a"; paramNamesUtf8 << UTF8_QSTRING("n₁") << UTF8_QSTRING("n₂") << "t" << "A"; break; case nsl_sf_stats_maxwell_boltzmann: paramNames << "s" << "a"; paramNamesUtf8 << UTF8_QSTRING("σ") << "A"; break; case nsl_sf_stats_frechet: paramNames << "g" << "mu" << "s" << "a"; paramNamesUtf8 << UTF8_QSTRING("γ") << UTF8_QSTRING("μ") << UTF8_QSTRING("σ") << "A"; break; } break; case nsl_fit_model_custom: break; } if (paramNamesUtf8.isEmpty()) paramNamesUtf8 << paramNames; //resize the vector for the start values and set the elements to 1.0 //in case a custom model is used, do nothing, we take over the previous values if (modelCategory != nsl_fit_model_custom) { const int np = paramNames.size(); paramStartValues.resize(np); paramFixed.resize(np); paramLowerLimits.resize(np); paramUpperLimits.resize(np); for (int i = 0; i < np; ++i) { paramStartValues[i] = 1.0; paramFixed[i] = false; paramLowerLimits[i] = -std::numeric_limits::max(); paramUpperLimits[i] = std::numeric_limits::max(); } // set some model-dependent start values if (modelCategory == nsl_fit_model_distribution) { nsl_sf_stats_distribution type = (nsl_sf_stats_distribution)modelType; if (type == nsl_sf_stats_flat) paramStartValues[0] = -1.0; else if (type == nsl_sf_stats_frechet || type == nsl_sf_stats_levy || type == nsl_sf_stats_exponential_power) paramStartValues[1] = 0.0; else if (type == nsl_sf_stats_weibull || type == nsl_sf_stats_gumbel2) paramStartValues[2] = 0.0; else if (type == nsl_sf_stats_binomial || type == nsl_sf_stats_negative_binomial || type == nsl_sf_stats_pascal || type == nsl_sf_stats_geometric || type == nsl_sf_stats_logarithmic) paramStartValues[0] = 0.5; } } } /*! Returns an icon to be used in the project explorer. */ QIcon XYFitCurve::icon() const { return QIcon::fromTheme("labplot-xy-fit-curve"); } //############################################################################## //########################## getter methods ################################## //############################################################################## BASIC_SHARED_D_READER_IMPL(XYFitCurve, const AbstractColumn*, xErrorColumn, xErrorColumn) BASIC_SHARED_D_READER_IMPL(XYFitCurve, const AbstractColumn*, yErrorColumn, yErrorColumn) const QString& XYFitCurve::xErrorColumnPath() const { Q_D(const XYFitCurve); return d->xErrorColumnPath; } const QString& XYFitCurve::yErrorColumnPath() const { Q_D(const XYFitCurve); return d->yErrorColumnPath; } BASIC_SHARED_D_READER_IMPL(XYFitCurve, XYFitCurve::FitData, fitData, fitData) const XYFitCurve::FitResult& XYFitCurve::fitResult() const { Q_D(const XYFitCurve); return d->fitResult; } //############################################################################## //################# setter methods and undo commands ########################## //############################################################################## STD_SETTER_CMD_IMPL_S(XYFitCurve, SetXErrorColumn, const AbstractColumn*, xErrorColumn) void XYFitCurve::setXErrorColumn(const AbstractColumn* column) { Q_D(XYFitCurve); if (column != d->xErrorColumn) { exec(new XYFitCurveSetXErrorColumnCmd(d, column, i18n("%1: assign x-error"))); handleSourceDataChanged(); if (column) { connect(column, SIGNAL(dataChanged(const AbstractColumn*)), this, SLOT(handleSourceDataChanged())); //TODO disconnect on undo } } } STD_SETTER_CMD_IMPL_S(XYFitCurve, SetYErrorColumn, const AbstractColumn*, yErrorColumn) void XYFitCurve::setYErrorColumn(const AbstractColumn* column) { Q_D(XYFitCurve); if (column != d->yErrorColumn) { exec(new XYFitCurveSetYErrorColumnCmd(d, column, i18n("%1: assign y-error"))); handleSourceDataChanged(); if (column) { connect(column, SIGNAL(dataChanged(const AbstractColumn*)), this, SLOT(handleSourceDataChanged())); //TODO disconnect on undo } } } STD_SETTER_CMD_IMPL_F_S(XYFitCurve, SetFitData, XYFitCurve::FitData, fitData, recalculate); void XYFitCurve::setFitData(const XYFitCurve::FitData& fitData) { Q_D(XYFitCurve); exec(new XYFitCurveSetFitDataCmd(d, fitData, i18n("%1: set fit options and perform the fit"))); } //############################################################################## //######################### Private implementation ############################# //############################################################################## XYFitCurvePrivate::XYFitCurvePrivate(XYFitCurve* owner) : XYAnalysisCurvePrivate(owner), xErrorColumn(nullptr), yErrorColumn(nullptr), residualsColumn(nullptr), residualsVector(nullptr), q(owner) { } XYFitCurvePrivate::~XYFitCurvePrivate() { //no need to delete xColumn and yColumn, they are deleted //when the parent aspect is removed } // data structure to pass parameter to fit functions struct data { size_t n; //number of data points double* x; //pointer to the vector with x-data values double* y; //pointer to the vector with y-data values double* weight; //pointer to the vector with weight values nsl_fit_model_category modelCategory; int modelType; int degree; QString* func; // string containing the definition of the model/function QStringList* paramNames; double* paramMin; // lower parameter limits double* paramMax; // upper parameter limits bool* paramFixed; // parameter fixed? }; /*! * \param paramValues vector containing current values of the fit parameters * \param params * \param f vector with the weighted residuals weight[i]*(Yi - y[i]) */ int func_f(const gsl_vector* paramValues, void* params, gsl_vector* f) { DEBUG("func_f"); size_t n = ((struct data*)params)->n; double* x = ((struct data*)params)->x; double* y = ((struct data*)params)->y; double* weight = ((struct data*)params)->weight; nsl_fit_model_category modelCategory = ((struct data*)params)->modelCategory; unsigned int modelType = ((struct data*)params)->modelType; QByteArray funcba = ((struct data*)params)->func->toLatin1(); // a local byte array is needed! const char *func = funcba.constData(); // function to evaluate QStringList* paramNames = ((struct data*)params)->paramNames; double *min = ((struct data*)params)->paramMin; double *max = ((struct data*)params)->paramMax; // set current values of the parameters for (int i = 0; i < paramNames->size(); i++) { double v = gsl_vector_get(paramValues, (size_t)i); // bound values if limits are set QByteArray paramnameba = paramNames->at(i).toLatin1(); assign_variable(paramnameba.constData(), nsl_fit_map_bound(v, min[i], max[i])); QDEBUG("Parameter"<n; double* xVector = ((struct data*)params)->x; double* weight = ((struct data*)params)->weight; nsl_fit_model_category modelCategory = ((struct data*)params)->modelCategory; unsigned int modelType = ((struct data*)params)->modelType; unsigned int degree = ((struct data*)params)->degree; QStringList* paramNames = ((struct data*)params)->paramNames; double *min = ((struct data*)params)->paramMin; double *max = ((struct data*)params)->paramMax; bool *fixed = ((struct data*)params)->paramFixed; // calculate the Jacobian matrix: // Jacobian matrix J(i,j) = df_i / dx_j // where f_i = w_i*(Y_i - y_i), // Y_i = model and the x_j are the parameters double x; switch (modelCategory) { case nsl_fit_model_basic: switch (modelType) { case nsl_fit_model_polynomial: // Y(x) = c0 + c1*x + ... + cn*x^n for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < (unsigned int)paramNames->size(); ++j) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_polynomial_param_deriv(x, j, weight[i])); } } break; case nsl_fit_model_power: // Y(x) = a*x^b or Y(x) = a + b*x^c. if (degree == 1) { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_power1_param_deriv(j, x, a, b, weight[i])); } } } else if (degree == 2) { const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double c = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_power2_param_deriv(j, x, b, c, weight[i])); } } } break; case nsl_fit_model_exponential: { // Y(x) = a*exp(b*x) + c*exp(d*x) + ... double *p = new double[2*degree]; for (unsigned int i = 0; i < 2*degree; i++) p[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, i), min[i], max[i]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2*degree; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exponentialn_param_deriv(j, x, p, weight[i])); } } delete[] p; break; } case nsl_fit_model_inverse_exponential: { // Y(x) = a*(1-exp(b*x))+c const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_inverse_exponential_param_deriv(j, x, a, b, weight[i])); } } break; } case nsl_fit_model_fourier: { // Y(x) = a0 + (a1*cos(w*x) + b1*sin(w*x)) + ... + (an*cos(n*w*x) + bn*sin(n*w*x) //parameters: w, a0, a1, b1, ... an, bn double* a = new double[degree]; double* b = new double[degree]; double w = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); a[0] = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); b[0] = 0; for (unsigned int i = 1; i < degree; ++i) { a[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, 2*i), min[2*i], max[2*i]); b[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, 2*i+1), min[2*i+1], max[2*i+1]); } for (size_t i = 0; i < n; i++) { x = xVector[i]; double wd = 0; //first derivative with respect to the w parameter for (unsigned int j = 1; j < degree; ++j) { wd += -a[j]*j*x*sin(j*w*x) + b[j]*j*x*cos(j*w*x); } gsl_matrix_set(J, i, 0, weight[i]*wd); gsl_matrix_set(J, i, 1, weight[i]); for (unsigned int j = 1; j <= degree; ++j) { gsl_matrix_set(J, (size_t)i, (size_t)(2*j), nsl_fit_model_fourier_param_deriv(0, j, x, w, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(2*j+1), nsl_fit_model_fourier_param_deriv(1, j, x, w, weight[i])); } for (unsigned int j = 0; j <= 2*degree+1; j++) if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } delete[] a; delete[] b; break; } } break; case nsl_fit_model_peak: switch (modelType) { case nsl_fit_model_gaussian: case nsl_fit_model_lorentz: case nsl_fit_model_sech: case nsl_fit_model_logistic: { for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < degree; ++j) { const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 3*j), min[3*j], max[3*j]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3*j+1), min[3*j+1], max[3*j+1]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 3*j+2), min[3*j+2], max[3*j+2]); switch (modelType) { case nsl_fit_model_gaussian: gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_gaussian_param_deriv(0, x, s, mu, a, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+1), nsl_fit_model_gaussian_param_deriv(1, x, s, mu, a, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+2), nsl_fit_model_gaussian_param_deriv(2, x, s, mu, a, weight[i])); break; case nsl_fit_model_lorentz: // s,t,a gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_lorentz_param_deriv(0, x, s, mu, a, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+1), nsl_fit_model_lorentz_param_deriv(1, x, s, mu, a, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+2), nsl_fit_model_lorentz_param_deriv(2, x, s, mu, a, weight[i])); break; case nsl_fit_model_sech: gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_sech_param_deriv(0, x, s, mu, a, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+1), nsl_fit_model_sech_param_deriv(1, x, s, mu, a, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+2), nsl_fit_model_sech_param_deriv(2, x, s, mu, a, weight[i])); break; case nsl_fit_model_logistic: gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_logistic_param_deriv(0, x, s, mu, a, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+1), nsl_fit_model_logistic_param_deriv(1, x, s, mu, a, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+2), nsl_fit_model_logistic_param_deriv(2, x, s, mu, a, weight[i])); break; } } for (unsigned int j = 0; j < 3*degree; j++) if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } break; } } break; case nsl_fit_model_growth: switch (modelType) { case nsl_fit_model_atan: case nsl_fit_model_tanh: case nsl_fit_model_algebraic_sigmoid: case nsl_fit_model_sigmoid: // Y(x) = a/(1+exp(-k*(x-mu))) case nsl_fit_model_erf: case nsl_fit_model_hill: case nsl_fit_model_gompertz: // Y(x) = a*exp(-b*exp(-c*x)); case nsl_fit_model_gudermann: { const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { if ((nsl_fit_model_type_growth)modelType == nsl_fit_model_atan) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_atan_param_deriv(j, x, s, mu, a, weight[i])); else if ((nsl_fit_model_type_growth)modelType == nsl_fit_model_tanh) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_tanh_param_deriv(j, x, s, mu, a, weight[i])); else if ((nsl_fit_model_type_growth)modelType == nsl_fit_model_algebraic_sigmoid) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_algebraic_sigmoid_param_deriv(j, x, s, mu, a, weight[i])); else if ((nsl_fit_model_type_growth)modelType == nsl_fit_model_sigmoid) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_sigmoid_param_deriv(j, x, s, mu, a, weight[i])); else if ((nsl_fit_model_type_growth)modelType == nsl_fit_model_erf) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_erf_param_deriv(j, x, s, mu, a, weight[i])); else if ((nsl_fit_model_type_growth)modelType == nsl_fit_model_hill) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_hill_param_deriv(j, x, s, mu, a, weight[i])); else if ((nsl_fit_model_type_growth)modelType == nsl_fit_model_gompertz) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gompertz_param_deriv(j, x, s, mu, a, weight[i])); else if ((nsl_fit_model_type_growth)modelType == nsl_fit_model_gudermann) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gudermann_param_deriv(j, x, s, mu, a, weight[i])); } } } break; } } break; case nsl_fit_model_distribution: switch (modelType) { case nsl_sf_stats_gaussian: case nsl_sf_stats_exponential: case nsl_sf_stats_laplace: case nsl_sf_stats_cauchy_lorentz: case nsl_sf_stats_rayleigh_tail: case nsl_sf_stats_lognormal: case nsl_sf_stats_logistic: case nsl_sf_stats_sech: case nsl_sf_stats_levy: { const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { switch (modelType) { case nsl_sf_stats_gaussian: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gaussian_param_deriv(j, x, s, mu, a, weight[i])); break; case nsl_sf_stats_exponential: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exponential_param_deriv(j, x, s, mu, a, weight[i])); break; case nsl_sf_stats_laplace: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_laplace_param_deriv(j, x, s, mu, a, weight[i])); break; case nsl_sf_stats_cauchy_lorentz: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_lorentz_param_deriv(j, x, s, mu, a, weight[i])); break; case nsl_sf_stats_rayleigh_tail: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_rayleigh_tail_param_deriv(j, x, s, mu, a, weight[i])); break; case nsl_sf_stats_lognormal: if (x > 0) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_lognormal_param_deriv(j, x, s, mu, a, weight[i])); else gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); break; case nsl_sf_stats_logistic: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_logistic_param_deriv(j, x, s, mu, a, weight[i])); break; case nsl_sf_stats_sech: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_sech_dist_param_deriv(j, x, s, mu, a, weight[i])); break; case nsl_sf_stats_levy: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_levy_param_deriv(j, x, s, mu, a, weight[i])); break; } } } } break; } case nsl_sf_stats_gaussian_tail: { const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gaussian_tail_param_deriv(j, x, s, mu, A, a, weight[i])); } } break; } case nsl_sf_stats_exponential_power: { const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exp_pow_param_deriv(j, x, s, mu, b, a, weight[i])); } } break; } case nsl_sf_stats_rayleigh: { const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_rayleigh_param_deriv(j, x, s, a, weight[i])); } } break; } case nsl_sf_stats_gamma: { const double t = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double k = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gamma_param_deriv(j, x, t, k, a, weight[i])); } } break; } case nsl_sf_stats_flat: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_flat_param_deriv(j, x, a, b, A, weight[i])); } } break; } case nsl_sf_stats_chi_squared: { const double nu = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_chi_square_param_deriv(j, x, nu, a, weight[i])); } } break; } case nsl_sf_stats_tdist: { const double nu = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_students_t_param_deriv(j, x, nu, a, weight[i])); } } break; } case nsl_sf_stats_fdist: { const double n1 = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double n2 = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_fdist_param_deriv(j, x, n1, n2, a, weight[i])); } } break; } case nsl_sf_stats_beta: case nsl_sf_stats_pareto: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { switch (modelType) { case nsl_sf_stats_beta: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_beta_param_deriv(j, x, a, b, A, weight[i])); break; case nsl_sf_stats_pareto: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_pareto_param_deriv(j, x, a, b, A, weight[i])); break; } } } } break; } case nsl_sf_stats_weibull: { const double k = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double l = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { if (x > 0) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_weibull_param_deriv(j, x, k, l, mu, a, weight[i])); else gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } } } break; } case nsl_sf_stats_gumbel1: { const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gumbel1_param_deriv(j, x, s, b, mu, a, weight[i])); } } break; } case nsl_sf_stats_gumbel2: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gumbel2_param_deriv(j, x, a, b, mu, A, weight[i])); } } break; } case nsl_sf_stats_poisson: { const double l = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_poisson_param_deriv(j, x, l, a, weight[i])); } } break; } case nsl_sf_stats_maxwell_boltzmann: { // Y(x) = a*sqrt(2/pi) * x^2/s^3 * exp(-(x/s)^2/2) const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_maxwell_param_deriv(j, x, s, a, weight[i])); } } break; } case nsl_sf_stats_frechet: { const double g = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_frechet_param_deriv(j, x, g, mu, s, a, weight[i])); } } break; } case nsl_sf_stats_landau: { // const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); for (size_t i = 0; i < n; i++) { x = xVector[i]; if (fixed[0]) gsl_matrix_set(J, (size_t)i, 0, 0.); else gsl_matrix_set(J, (size_t)i, 0, nsl_fit_model_landau_param_deriv(0, x, weight[i])); } break; } case nsl_sf_stats_binomial: case nsl_sf_stats_negative_binomial: case nsl_sf_stats_pascal: { const double p = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double N = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { switch (modelType) { case nsl_sf_stats_binomial: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_binomial_param_deriv(j, x, p, N, a, weight[i])); break; case nsl_sf_stats_negative_binomial: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_negative_binomial_param_deriv(j, x, p, N, a, weight[i])); break; case nsl_sf_stats_pascal: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_pascal_param_deriv(j, x, p, N, a, weight[i])); break; } } } } break; } case nsl_sf_stats_geometric: case nsl_sf_stats_logarithmic: { const double p = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { switch (modelType) { case nsl_sf_stats_geometric: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_geometric_param_deriv(j, x, p, a, weight[i])); break; case nsl_sf_stats_logarithmic: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_logarithmic_param_deriv(j, x, p, a, weight[i])); break; } } } } break; } case nsl_sf_stats_hypergeometric: { const double n1 = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double n2 = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double t = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_hypergeometric_param_deriv(j, x, n1, n2, t, a, weight[i])); } } break; } // unused distributions case nsl_sf_stats_levy_alpha_stable: case nsl_sf_stats_levy_skew_alpha_stable: case nsl_sf_stats_bernoulli: break; } break; case nsl_fit_model_custom: QByteArray funcba = ((struct data*)params)->func->toLatin1(); const char* func = funcba.data(); QByteArray nameba; double value; const unsigned int np = paramNames->size(); for (size_t i = 0; i < n; i++) { x = xVector[i]; assign_variable("x", x); for (unsigned int j = 0; j < np; j++) { for (unsigned int k = 0; k < np; k++) { if (k != j) { nameba = paramNames->at(k).toLatin1(); value = nsl_fit_map_bound(gsl_vector_get(paramValues, k), min[k], max[k]); assign_variable(nameba.data(), value); } } nameba = paramNames->at(j).toLatin1(); const char *name = nameba.data(); value = nsl_fit_map_bound(gsl_vector_get(paramValues, j), min[j], max[j]); assign_variable(name, value); const double f_p = parse(func); const double eps = 1.e-9 * std::abs(f_p); // adapt step size to value value += eps; assign_variable(name, value); const double f_pdp = parse(func); // qDebug()<<"evaluate deriv"<* >(xColumn->data()); yVector = static_cast* >(yColumn->data()); residualsVector = static_cast* >(residualsColumn->data()); xColumn->setHidden(true); q->addChild(xColumn); yColumn->setHidden(true); q->addChild(yColumn); q->addChild(residualsColumn); q->setUndoAware(false); q->setXColumn(xColumn); q->setYColumn(yColumn); q->setUndoAware(true); } else { xVector->clear(); yVector->clear(); residualsVector->clear(); } // clear the previous result fitResult = XYFitCurve::FitResult(); //fit settings const unsigned int maxIters = fitData.maxIterations; //maximal number of iterations const double delta = fitData.eps; //fit tolerance const unsigned int np = fitData.paramNames.size(); //number of fit parameters if (np == 0) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("Model has no parameters."); emit (q->dataChanged()); sourceDataChangedSinceLastRecalc = false; return; } //determine the data source columns const AbstractColumn* tmpXDataColumn = 0; const AbstractColumn* tmpYDataColumn = 0; if (dataSourceType == XYAnalysisCurve::DataSourceSpreadsheet) { //spreadsheet columns as data source DEBUG(" spreadsheet columns as data source"); tmpXDataColumn = xDataColumn; tmpYDataColumn = yDataColumn; } else { //curve columns as data source DEBUG(" curve columns as data source"); tmpXDataColumn = dataSourceCurve->xColumn(); tmpYDataColumn = dataSourceCurve->yColumn(); } if (!tmpXDataColumn || !tmpYDataColumn) { emit (q->dataChanged()); sourceDataChangedSinceLastRecalc = false; return; } //check column sizes if (tmpXDataColumn->rowCount() != tmpYDataColumn->rowCount()) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("Number of x and y data points must be equal."); emit (q->dataChanged()); sourceDataChangedSinceLastRecalc = false; return; } if (yErrorColumn) { if (yErrorColumn->rowCount() < tmpXDataColumn->rowCount()) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("Not sufficient weight data points provided."); emit (q->dataChanged()); sourceDataChangedSinceLastRecalc = false; return; } } //copy all valid data point for the fit to temporary vectors QVector xdataVector; QVector ydataVector; QVector xerrorVector; QVector yerrorVector; double xmin, xmax; if (fitData.autoRange) { xmin = tmpXDataColumn->minimum(); xmax = tmpXDataColumn->maximum(); } else { xmin = fitData.xRange.first(); xmax = fitData.xRange.last(); } for (int row = 0; row < tmpXDataColumn->rowCount(); ++row) { //only copy those data where _all_ values (for x and y and errors, if given) are valid if (!std::isnan(tmpXDataColumn->valueAt(row)) && !std::isnan(tmpYDataColumn->valueAt(row)) && !tmpXDataColumn->isMasked(row) && !tmpYDataColumn->isMasked(row)) { // only when inside given range if (tmpXDataColumn->valueAt(row) >= xmin && tmpXDataColumn->valueAt(row) <= xmax) { if (dataSourceType == XYAnalysisCurve::DataSourceCurve || (!xErrorColumn && !yErrorColumn) || !fitData.useDataErrors) { // x-y xdataVector.append(tmpXDataColumn->valueAt(row)); ydataVector.append(tmpYDataColumn->valueAt(row)); } else if (!xErrorColumn) { // x-y-dy if (!std::isnan(yErrorColumn->valueAt(row))) { xdataVector.append(tmpXDataColumn->valueAt(row)); ydataVector.append(tmpYDataColumn->valueAt(row)); yerrorVector.append(yErrorColumn->valueAt(row)); } } else { // x-y-dx-dy if (!std::isnan(xErrorColumn->valueAt(row)) && !std::isnan(yErrorColumn->valueAt(row))) { xdataVector.append(tmpXDataColumn->valueAt(row)); ydataVector.append(tmpYDataColumn->valueAt(row)); xerrorVector.append(xErrorColumn->valueAt(row)); yerrorVector.append(yErrorColumn->valueAt(row)); } } } } } //number of data points to fit const size_t n = xdataVector.size(); DEBUG("number of data points: " << n); if (n == 0) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("No data points available."); emit (q->dataChanged()); sourceDataChangedSinceLastRecalc = false; return; } if (n < np) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("The number of data points (%1) must be greater than or equal to the number of parameters (%2).", n, np); emit (q->dataChanged()); sourceDataChangedSinceLastRecalc = false; return; } double* xdata = xdataVector.data(); double* ydata = ydataVector.data(); double* xerror = xerrorVector.data(); // size may be 0 double* yerror = yerrorVector.data(); // size may be 0 DEBUG("x errors: " << xerrorVector.size()); DEBUG("y errors: " << yerrorVector.size()); double* weight = new double[n]; switch(fitData.xErrorsType) { case nsl_fit_error_no: case nsl_fit_error_direct: break; case nsl_fit_error_inverse: for(int i = 0; i < xerrorVector.size(); i++) xerror[i] = 1./xerror[i]; break; } for (size_t i = 0; i < n; i++) weight[i] = 1.; switch (fitData.yWeightsType) { case nsl_fit_weight_no: case nsl_fit_weight_statistical_fit: case nsl_fit_weight_relative_fit: break; case nsl_fit_weight_instrumental: for(int i = 0; i < (int)n; i++) if (i < yerrorVector.size()) weight[i] = 1./gsl_pow_2(yerror[i]); break; case nsl_fit_weight_direct: for(int i = 0; i < (int)n; i++) if (i < yerrorVector.size()) weight[i] = yerror[i]; break; case nsl_fit_weight_inverse: for(int i = 0; i < (int)n; i++) if (i < yerrorVector.size()) weight[i] = 1./yerror[i]; break; case nsl_fit_weight_statistical: for (int i = 0; i < (int)n; i++) weight[i] = 1./ydata[i]; break; case nsl_fit_weight_relative: for (int i = 0; i < (int)n; i++) weight[i] = 1./gsl_pow_2(ydata[i]); break; } /////////////////////// GSL >= 2 has a complete new interface! But the old one is still supported. /////////////////////////// // GSL >= 2 : "the 'fdf' field of gsl_multifit_function_fdf is now deprecated and does not need to be specified for nonlinear least squares problems" for (unsigned int i = 0; i < np; i++) DEBUG("parameter " << i << " fixed: " << fitData.paramFixed.data()[i]); //function to fit gsl_multifit_function_fdf f; DEBUG("model = " << fitData.model.toStdString()); struct data params = {n, xdata, ydata, weight, fitData.modelCategory, fitData.modelType, fitData.degree, &fitData.model, &fitData.paramNames, fitData.paramLowerLimits.data(), fitData.paramUpperLimits.data(), fitData.paramFixed.data()}; f.f = &func_f; f.df = &func_df; f.fdf = &func_fdf; f.n = n; f.p = np; f.params = ¶ms; DEBUG("initialize the derivative solver (using Levenberg-Marquardt robust solver)"); const gsl_multifit_fdfsolver_type* T = gsl_multifit_fdfsolver_lmsder; gsl_multifit_fdfsolver* s = gsl_multifit_fdfsolver_alloc(T, n, np); DEBUG("set start values"); double* x_init = fitData.paramStartValues.data(); double* x_min = fitData.paramLowerLimits.data(); double* x_max = fitData.paramUpperLimits.data(); DEBUG("scale start values if limits are set"); for (unsigned int i = 0; i < np; i++) x_init[i] = nsl_fit_map_unbound(x_init[i], x_min[i], x_max[i]); DEBUG(" DONE"); gsl_vector_view x = gsl_vector_view_array(x_init, np); DEBUG("Turn off GSL error handler to avoid overflow/underflow"); gsl_set_error_handler_off(); DEBUG("initialize solver with function f and initial guess x"); gsl_multifit_fdfsolver_set(s, &f, &x.vector); DEBUG("Iterate ..."); int status; unsigned int iter = 0; fitResult.solverOutput.clear(); writeSolverState(s); do { iter++; // update weights for Y-depending weights if (fitData.yWeightsType == nsl_fit_weight_statistical_fit) { for (size_t i = 0; i < n; i++) weight[i] = 1./(gsl_vector_get(s->f, i) + ydata[i]); // 1/Y_i } else if (fitData.yWeightsType == nsl_fit_weight_relative_fit) { for (size_t i = 0; i < n; i++) weight[i] = 1./gsl_pow_2(gsl_vector_get(s->f, i) + ydata[i]); // 1/Y_i^2 } status = gsl_multifit_fdfsolver_iterate(s); writeSolverState(s); if (status) { DEBUG("iter " << iter << ", status = " << gsl_strerror(status)); break; } status = gsl_multifit_test_delta(s->dx, s->x, delta, delta); } while (status == GSL_CONTINUE && iter < maxIters); // second run for x-error fitting if (xerrorVector.size() > 0) { DEBUG("Rerun fit with x errors"); // y'(x) double *yd = new double[n]; for (size_t i = 0; i < n; i++) { size_t index = i; if (index == n-1) index = n-2; yd[i] = gsl_vector_get(s->f, index+1) + ydata[index+1] - gsl_vector_get(s->f, index) - ydata[index]; yd[i] /= (xdata[index+1] - xdata[index]); } switch (fitData.yWeightsType) { case nsl_fit_weight_no: break; case nsl_fit_weight_instrumental: for (size_t i = 0; i < n; i++) { double sigma; if (yerrorVector.size() > 0) // x- and y-error // sigma = sqrt(sigma_y^2 + (y'(x)*sigma_x)^2) sigma = sqrt(gsl_pow_2(yerror[i]) + gsl_pow_2(yd[i] * xerror[i])); else // only x-error sigma = yd[i] * xerror[i]; weight[i] = 1./gsl_pow_2(sigma); } break; // other weight types: y'(x) considered correctly? case nsl_fit_weight_direct: for (size_t i = 0; i < n; i++) { weight[i] = xerror[i]/yd[i]; if (yerrorVector.size() > 0) weight[i] += yerror[i]; } break; case nsl_fit_weight_inverse: for (size_t i = 0; i < n; i++) { weight[i] = yd[i]/xerror[i]; if (yerrorVector.size() > 0) weight[i] += 1./yerror[i]; } break; case nsl_fit_weight_statistical: case nsl_fit_weight_relative: break; case nsl_fit_weight_statistical_fit: for (size_t i = 0; i < n; i++) weight[i] = 1./(gsl_vector_get(s->f, i) + ydata[i]); // 1/Y_i break; case nsl_fit_weight_relative_fit: for (size_t i = 0; i < n; i++) weight[i] = 1./gsl_pow_2(gsl_vector_get(s->f, i) + ydata[i]); // 1/Y_i^2 break; } delete[] yd; do { iter++; status = gsl_multifit_fdfsolver_iterate(s); writeSolverState(s); if (status) break; status = gsl_multifit_test_delta(s->dx, s->x, delta, delta); } while (status == GSL_CONTINUE && iter < maxIters); } delete[] weight; // unscale start values for (unsigned int i = 0; i < np; i++) x_init[i] = nsl_fit_map_bound(x_init[i], x_min[i], x_max[i]); //get the covariance matrix //TODO: scale the Jacobian when limits are used before constructing the covar matrix? gsl_matrix* covar = gsl_matrix_alloc(np, np); #if GSL_MAJOR_VERSION >= 2 // the Jacobian is not part of the solver anymore gsl_matrix *J = gsl_matrix_alloc(s->fdf->n, s->fdf->p); gsl_multifit_fdfsolver_jac(s, J); gsl_multifit_covar(J, 0.0, covar); #else gsl_multifit_covar(s->J, 0.0, covar); #endif //write the result fitResult.available = true; fitResult.valid = true; fitResult.status = gslErrorToString(status); fitResult.iterations = iter; fitResult.dof = n - np; //gsl_blas_dnrm2() - computes the Euclidian norm (||r||_2 = \sqrt {\sum r_i^2}) of the vector with the elements weight[i]*(Yi - y[i]) //gsl_blas_dasum() - computes the absolute sum \sum |r_i| of the elements of the vector with the elements weight[i]*(Yi - y[i]) fitResult.sse = gsl_pow_2(gsl_blas_dnrm2(s->f)); if (fitResult.dof != 0) { fitResult.rms = fitResult.sse/fitResult.dof; fitResult.rsd = sqrt(fitResult.rms); } fitResult.mse = fitResult.sse/n; fitResult.rmse = sqrt(fitResult.mse); fitResult.mae = gsl_blas_dasum(s->f)/n; //needed for coefficient of determination, R-squared fitResult.sst = gsl_stats_tss(ydata, 1, n); fitResult.rsquare = nsl_stats_rsquare(fitResult.sse, fitResult.sst); fitResult.rsquareAdj = nsl_stats_rsquareAdj(fitResult.rsquare, np, fitResult.dof); fitResult.chisq_p = nsl_stats_chisq_p(fitResult.sse, fitResult.dof); fitResult.fdist_F = nsl_stats_fdist_F(fitResult.sst, fitResult.rms); fitResult.fdist_p = nsl_stats_fdist_p(fitResult.fdist_F, np, fitResult.dof); fitResult.aic = nsl_stats_aic(fitResult.sse, n, np); fitResult.bic = nsl_stats_bic(fitResult.sse, n, np); //parameter values const double c = GSL_MIN_DBL(1., sqrt(fitResult.rms)); //limit error for poor fit fitResult.paramValues.resize(np); fitResult.errorValues.resize(np); fitResult.tdist_tValues.resize(np); fitResult.tdist_pValues.resize(np); fitResult.tdist_marginValues.resize(np); for (unsigned int i = 0; i < np; i++) { // scale resulting values if they are bounded fitResult.paramValues[i] = nsl_fit_map_bound(gsl_vector_get(s->x, i), x_min[i], x_max[i]); // use results as start values if desired if (fitData.useResults) { fitData.paramStartValues.data()[i] = fitResult.paramValues[i]; DEBUG("saving parameter " << i << ": " << fitResult.paramValues[i] << ' ' << fitData.paramStartValues.data()[i]); } fitResult.errorValues[i] = c*sqrt(gsl_matrix_get(covar, i, i)); fitResult.tdist_tValues[i] = nsl_stats_tdist_t(fitResult.paramValues.at(i), fitResult.errorValues.at(i)); fitResult.tdist_pValues[i] = nsl_stats_tdist_p(fitResult.tdist_tValues.at(i), fitResult.dof); fitResult.tdist_marginValues[i] = nsl_stats_tdist_margin(0.05, fitResult.dof, fitResult.errorValues.at(i)); } // fill residuals vector. To get residuals on the correct x values, fill the rest with zeros. residualsVector->resize(tmpXDataColumn->rowCount()); if (fitData.evaluateFullRange) { // evaluate full range of residuals xVector->resize(tmpXDataColumn->rowCount()); for (int i = 0; i < tmpXDataColumn->rowCount(); i++) (*xVector)[i] = tmpXDataColumn->valueAt(i); ExpressionParser* parser = ExpressionParser::getInstance(); bool rc = parser->evaluateCartesian(fitData.model, xVector, residualsVector, fitData.paramNames, fitResult.paramValues); for (int i = 0; i < tmpXDataColumn->rowCount(); i++) (*residualsVector)[i] = tmpYDataColumn->valueAt(i) - (*residualsVector)[i]; if (!rc) residualsVector->clear(); } else { // only selected range size_t j = 0; for (int i = 0; i < tmpXDataColumn->rowCount(); i++) { if (tmpXDataColumn->valueAt(i) >= xmin && tmpXDataColumn->valueAt(i) <= xmax) residualsVector->data()[i] = - gsl_vector_get(s->f, j++); else // outside range residualsVector->data()[i] = 0; } } residualsColumn->setChanged(); //free resources gsl_multifit_fdfsolver_free(s); gsl_matrix_free(covar); //calculate the fit function (vectors) ExpressionParser* parser = ExpressionParser::getInstance(); if (fitData.evaluateFullRange) { // evaluate fit on full data range if selected xmin = tmpXDataColumn->minimum(); xmax = tmpXDataColumn->maximum(); } xVector->resize((int)fitData.evaluatedPoints); yVector->resize((int)fitData.evaluatedPoints); bool rc = parser->evaluateCartesian(fitData.model, QString::number(xmin), QString::number(xmax), (int)fitData.evaluatedPoints, xVector, yVector, fitData.paramNames, fitResult.paramValues); if (!rc) { xVector->clear(); yVector->clear(); } fitResult.elapsedTime = timer.elapsed(); //redraw the curve emit (q->dataChanged()); sourceDataChangedSinceLastRecalc = false; } /*! * writes out the current state of the solver \c s */ void XYFitCurvePrivate::writeSolverState(gsl_multifit_fdfsolver* s) { QString state; //current parameter values, semicolon separated double* min = fitData.paramLowerLimits.data(); double* max = fitData.paramUpperLimits.data(); for (int i = 0; i < fitData.paramNames.size(); ++i) { const double x = gsl_vector_get(s->x, i); // map parameter if bounded state += QString::number(nsl_fit_map_bound(x, min[i], max[i])) + '\t'; } //current value of the chi2-function state += QString::number(gsl_pow_2(gsl_blas_dnrm2(s->f))); state += ';'; fitResult.solverOutput += state; } //############################################################################## //################## Serialization/Deserialization ########################### //############################################################################## //! Save as XML void XYFitCurve::save(QXmlStreamWriter* writer) const { Q_D(const XYFitCurve); writer->writeStartElement("xyFitCurve"); //write the base class XYAnalysisCurve::save(writer); //write xy-fit-curve specific information //fit data - only save model expression and parameter names for custom model, otherwise they are set in XYFitCurve::initFitData() writer->writeStartElement("fitData"); WRITE_COLUMN(d->xErrorColumn, xErrorColumn); WRITE_COLUMN(d->yErrorColumn, yErrorColumn); writer->writeAttribute("autoRange", QString::number(d->fitData.autoRange)); writer->writeAttribute("xRangeMin", QString::number(d->fitData.xRange.first(), 'g', 15)); writer->writeAttribute("xRangeMax", QString::number(d->fitData.xRange.last(), 'g', 15)); writer->writeAttribute("modelCategory", QString::number(d->fitData.modelCategory)); writer->writeAttribute("modelType", QString::number(d->fitData.modelType)); writer->writeAttribute("xErrorsType", QString::number(d->fitData.xErrorsType)); writer->writeAttribute("weightsType", QString::number(d->fitData.yWeightsType)); writer->writeAttribute("degree", QString::number(d->fitData.degree)); if (d->fitData.modelCategory == nsl_fit_model_custom) writer->writeAttribute("model", d->fitData.model); writer->writeAttribute("maxIterations", QString::number(d->fitData.maxIterations)); writer->writeAttribute("eps", QString::number(d->fitData.eps, 'g', 15)); writer->writeAttribute("evaluatedPoints", QString::number(d->fitData.evaluatedPoints)); writer->writeAttribute("evaluateFullRange", QString::number(d->fitData.evaluateFullRange)); writer->writeAttribute("useDataErrors", QString::number(d->fitData.useDataErrors)); writer->writeAttribute("useResults", QString::number(d->fitData.useResults)); if (d->fitData.modelCategory == nsl_fit_model_custom) { writer->writeStartElement("paramNames"); foreach (const QString &name, d->fitData.paramNames) writer->writeTextElement("name", name); writer->writeEndElement(); } writer->writeStartElement("paramStartValues"); foreach (const double &value, d->fitData.paramStartValues) writer->writeTextElement("startValue", QString::number(value, 'g', 15)); writer->writeEndElement(); // use 16 digits to handle -DBL_MAX writer->writeStartElement("paramLowerLimits"); foreach (const double &limit, d->fitData.paramLowerLimits) writer->writeTextElement("lowerLimit", QString::number(limit, 'g', 16)); writer->writeEndElement(); // use 16 digits to handle DBL_MAX writer->writeStartElement("paramUpperLimits"); foreach (const double &limit, d->fitData.paramUpperLimits) writer->writeTextElement("upperLimit", QString::number(limit, 'g', 16)); writer->writeEndElement(); writer->writeStartElement("paramFixed"); foreach (const double &fixed, d->fitData.paramFixed) writer->writeTextElement("fixed", QString::number(fixed)); writer->writeEndElement(); writer->writeEndElement(); //"fitData" //fit results (generated columns and goodness of the fit) writer->writeStartElement("fitResult"); writer->writeAttribute("available", QString::number(d->fitResult.available)); writer->writeAttribute("valid", QString::number(d->fitResult.valid)); writer->writeAttribute("status", d->fitResult.status); writer->writeAttribute("iterations", QString::number(d->fitResult.iterations)); writer->writeAttribute("time", QString::number(d->fitResult.elapsedTime)); writer->writeAttribute("dof", QString::number(d->fitResult.dof)); writer->writeAttribute("sse", QString::number(d->fitResult.sse, 'g', 15)); writer->writeAttribute("sst", QString::number(d->fitResult.sst, 'g', 15)); writer->writeAttribute("rms", QString::number(d->fitResult.rms, 'g', 15)); writer->writeAttribute("rsd", QString::number(d->fitResult.rsd, 'g', 15)); writer->writeAttribute("mse", QString::number(d->fitResult.mse, 'g', 15)); writer->writeAttribute("rmse", QString::number(d->fitResult.rmse, 'g', 15)); writer->writeAttribute("mae", QString::number(d->fitResult.mae, 'g', 15)); writer->writeAttribute("rsquare", QString::number(d->fitResult.rsquare, 'g', 15)); writer->writeAttribute("rsquareAdj", QString::number(d->fitResult.rsquareAdj, 'g', 15)); writer->writeAttribute("chisq_p", QString::number(d->fitResult.chisq_p, 'g', 15)); writer->writeAttribute("fdist_F", QString::number(d->fitResult.fdist_F, 'g', 15)); writer->writeAttribute("fdist_p", QString::number(d->fitResult.fdist_p, 'g', 15)); writer->writeAttribute("aic", QString::number(d->fitResult.aic, 'g', 15)); writer->writeAttribute("bic", QString::number(d->fitResult.bic, 'g', 15)); writer->writeAttribute("solverOutput", d->fitResult.solverOutput); writer->writeStartElement("paramValues"); foreach (const double &value, d->fitResult.paramValues) writer->writeTextElement("value", QString::number(value, 'g', 15)); writer->writeEndElement(); writer->writeStartElement("errorValues"); foreach (const double &value, d->fitResult.errorValues) writer->writeTextElement("error", QString::number(value, 'g', 15)); writer->writeEndElement(); writer->writeStartElement("tdist_tValues"); foreach (const double &value, d->fitResult.tdist_tValues) writer->writeTextElement("tdist_t", QString::number(value, 'g', 15)); writer->writeEndElement(); writer->writeStartElement("tdist_pValues"); foreach (const double &value, d->fitResult.tdist_pValues) writer->writeTextElement("tdist_p", QString::number(value, 'g', 15)); writer->writeEndElement(); writer->writeStartElement("tdist_marginValues"); foreach (const double &value, d->fitResult.tdist_marginValues) writer->writeTextElement("tdist_margin", QString::number(value, 'g', 15)); writer->writeEndElement(); //save calculated columns if available if (d->xColumn && d->yColumn && d->residualsColumn) { d->xColumn->save(writer); d->yColumn->save(writer); d->residualsColumn->save(writer); } writer->writeEndElement(); //"fitResult" writer->writeEndElement(); //"xyFitCurve" } //! Load from XML bool XYFitCurve::load(XmlStreamReader* reader, bool preview) { Q_D(XYFitCurve); QString attributeWarning = i18n("Attribute '%1' missing or empty, default value is used"); QXmlStreamAttributes attribs; QString str; while (!reader->atEnd()) { reader->readNext(); if (reader->isEndElement() && reader->name() == "xyFitCurve") break; if (!reader->isStartElement()) continue; if (reader->name() == "xyAnalysisCurve") { if ( !XYAnalysisCurve::load(reader, preview) ) return false; } else if (!preview && reader->name() == "fitData") { attribs = reader->attributes(); READ_COLUMN(xErrorColumn); READ_COLUMN(yErrorColumn); READ_INT_VALUE("autoRange", fitData.autoRange, bool); READ_DOUBLE_VALUE("xRangeMin", fitData.xRange.first()); READ_DOUBLE_VALUE("xRangeMax", fitData.xRange.last()); READ_INT_VALUE("modelCategory", fitData.modelCategory, nsl_fit_model_category); READ_INT_VALUE("modelType", fitData.modelType, unsigned int); READ_INT_VALUE("xErrorsType", fitData.xErrorsType, nsl_fit_error_type); READ_INT_VALUE("weightsType", fitData.yWeightsType, nsl_fit_weight_type); READ_INT_VALUE("degree", fitData.degree, int); if (d->fitData.modelCategory == nsl_fit_model_custom) { READ_STRING_VALUE("model", fitData.model); DEBUG("read model = " << d->fitData.model.toStdString()); } READ_INT_VALUE("maxIterations", fitData.maxIterations, int); READ_DOUBLE_VALUE("eps", fitData.eps); READ_INT_VALUE("fittedPoints", fitData.evaluatedPoints, size_t); // old name READ_INT_VALUE("evaluatedPoints", fitData.evaluatedPoints, size_t); READ_INT_VALUE("evaluateFullRange", fitData.evaluateFullRange, bool); READ_INT_VALUE("useDataErrors", fitData.useDataErrors, bool); READ_INT_VALUE("useResults", fitData.useResults, bool); //set the model expression and the parameter names (can be derived from the saved values for category, type and degree) XYFitCurve::initFitData(d->fitData); } else if (!preview && reader->name() == "name") { // needed for custom model d->fitData.paramNames << reader->readElementText(); } else if (!preview && reader->name() == "startValue") { d->fitData.paramStartValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "fixed") { d->fitData.paramFixed << (bool)reader->readElementText().toInt(); } else if (!preview && reader->name() == "lowerLimit") { bool ok; double x = reader->readElementText().toDouble(&ok); if (ok) // -DBL_MAX results in conversion error d->fitData.paramLowerLimits << x; else d->fitData.paramLowerLimits << -std::numeric_limits::max(); } else if (!preview && reader->name() == "upperLimit") { bool ok; double x = reader->readElementText().toDouble(&ok); if (ok) // DBL_MAX results in conversion error d->fitData.paramUpperLimits << x; else d->fitData.paramUpperLimits << std::numeric_limits::max(); } else if (!preview && reader->name() == "value") { d->fitResult.paramValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "error") { d->fitResult.errorValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "tdist_t") { d->fitResult.tdist_tValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "tdist_p") { d->fitResult.tdist_pValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "tdist_margin") { d->fitResult.tdist_marginValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "fitResult") { attribs = reader->attributes(); READ_INT_VALUE("available", fitResult.available, int); READ_INT_VALUE("valid", fitResult.valid, int); READ_STRING_VALUE("status", fitResult.status); READ_INT_VALUE("iterations", fitResult.iterations, int); READ_INT_VALUE("time", fitResult.elapsedTime, int); READ_DOUBLE_VALUE("dof", fitResult.dof); READ_DOUBLE_VALUE("sse", fitResult.sse); READ_DOUBLE_VALUE("sst", fitResult.sst); READ_DOUBLE_VALUE("rms", fitResult.rms); READ_DOUBLE_VALUE("rsd", fitResult.rsd); READ_DOUBLE_VALUE("mse", fitResult.mse); READ_DOUBLE_VALUE("rmse", fitResult.rmse); READ_DOUBLE_VALUE("mae", fitResult.mae); READ_DOUBLE_VALUE("rsquare", fitResult.rsquare); READ_DOUBLE_VALUE("rsquareAdj", fitResult.rsquareAdj); READ_DOUBLE_VALUE("chisq_p", fitResult.chisq_p); READ_DOUBLE_VALUE("fdist_F", fitResult.fdist_F); READ_DOUBLE_VALUE("fdist_p", fitResult.fdist_p); READ_DOUBLE_VALUE("aic", fitResult.aic); READ_DOUBLE_VALUE("bic", fitResult.bic); READ_STRING_VALUE("solverOutput", fitResult.solverOutput); } else if (reader->name() == "column") { Column* column = new Column("", AbstractColumn::Numeric); if (!column->load(reader, preview)) { delete column; return false; } if (column->name() == "x") d->xColumn = column; else if (column->name() == "y") d->yColumn = column; else if (column->name() == "residuals") d->residualsColumn = column; } } if (preview) return true; // new fit model style (reset model type of old projects) if (d->fitData.modelCategory == nsl_fit_model_basic && d->fitData.modelType >= NSL_FIT_MODEL_BASIC_COUNT) { DEBUG("reset old fit model"); d->fitData.modelType = 0; d->fitData.degree = 1; // reset size of fields not touched by initFitData() d->fitData.paramStartValues.resize(2); d->fitData.paramFixed.resize(2); d->fitResult.paramValues.resize(2); d->fitResult.errorValues.resize(2); d->fitResult.tdist_tValues.resize(2); d->fitResult.tdist_pValues.resize(2); d->fitResult.tdist_marginValues.resize(2); } // not present in old projects if (d->fitResult.tdist_tValues.size() == 0) d->fitResult.tdist_tValues.resize(d->fitResult.paramValues.size()); if (d->fitResult.tdist_pValues.size() == 0) d->fitResult.tdist_pValues.resize(d->fitResult.paramValues.size()); if (d->fitResult.tdist_marginValues.size() == 0) d->fitResult.tdist_marginValues.resize(d->fitResult.paramValues.size()); // wait for data to be read before using the pointers QThreadPool::globalInstance()->waitForDone(); if (d->xColumn && d->yColumn && d->residualsColumn) { d->xColumn->setHidden(true); addChild(d->xColumn); d->yColumn->setHidden(true); addChild(d->yColumn); addChild(d->residualsColumn); d->xVector = static_cast* >(d->xColumn->data()); d->yVector = static_cast* >(d->yColumn->data()); d->residualsVector = static_cast* >(d->residualsColumn->data()); XYCurve::d_ptr->xColumn = d->xColumn; XYCurve::d_ptr->yColumn = d->yColumn; } return true; } diff --git a/tests/analysis/fit/FitTest.cpp b/tests/analysis/fit/FitTest.cpp index b2bfacf02..8d765a0e5 100644 --- a/tests/analysis/fit/FitTest.cpp +++ b/tests/analysis/fit/FitTest.cpp @@ -1,854 +1,908 @@ /*************************************************************************** File : FitTest.cpp Project : LabPlot Description : Tests for data fitting -------------------------------------------------------------------- Copyright : (C) 2017 Alexander Semke (alexander.semke@web.de) Copyright : (C) 2018 Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #include "FitTest.h" #include "backend/core/column/Column.h" #include "backend/worksheet/plots/cartesian/XYFitCurve.h" extern "C" { #include "backend/nsl/nsl_sf_stats.h" #include "backend/nsl/nsl_stats.h" } void FitTest::initTestCase() { // needed in order to have the signals triggered by SignallingUndoCommand, see LabPlot.cpp //TODO: redesign/remove this qRegisterMetaType("const AbstractAspect*"); qRegisterMetaType("const AbstractColumn*"); } //############################################################################## //################# linear regression with NIST datasets ###################### //############################################################################## void FitTest::testLinearNorris() { //NIST data for Norris dataset QVector xData = {0.2,337.4,118.2,884.6,10.1,226.5,666.3,996.3,448.6,777.0,558.2,0.4,0.6,775.5,666.9,338.0,447.5,11.6,556.0,228.1, 995.8,887.6,120.2,0.3,0.3,556.8,339.1,887.2,999.0,779.0,11.1,118.3,229.2,669.1,448.9,0.5}; QVector yData = {0.1,338.8,118.1,888.0,9.2,228.1,668.5,998.5,449.1,778.9,559.2,0.3,0.1,778.1,668.8,339.3,448.9,10.8,557.7,228.3, 998.0,888.8,119.6,0.3,0.6,557.6,339.3,888.0,998.5,778.9,10.2,117.6,228.9,668.4,449.2,0.2}; //data source columns Column xDataColumn("x", AbstractColumn::Numeric); xDataColumn.replaceValues(0, xData); Column yDataColumn("y", AbstractColumn::Numeric); yDataColumn.replaceValues(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_polynomial; fitData.degree = 1; XYFitCurve::initFitData(fitData); fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); const int np = fitData.paramNames.size(); QCOMPARE(np, 2); QCOMPARE(fitResult.paramValues.at(0), -0.262323073774029); QCOMPARE(fitResult.errorValues.at(0), 0.232818234301152); QCOMPARE(fitResult.paramValues.at(1), 1.00211681802045); QCOMPARE(fitResult.errorValues.at(1), 0.429796848199937e-3); QCOMPARE(fitResult.rsd, 0.884796396144373); QCOMPARE(fitResult.rsquare, 0.999993745883712); QCOMPARE(fitResult.sse, 26.6173985294224); QCOMPARE(fitResult.rms, 0.782864662630069); // DEBUG(std::setprecision(15) << fitResult.fdist_F); // result: 5436419.54079774 FuzzyCompare(fitResult.fdist_F, 5436385.54079785, 34.); } void FitTest::testLinearPontius() { //NIST data for Pontius dataset QVector xData = {150000,300000,450000,600000,750000,900000,1050000,1200000,1350000,1500000,1650000,1800000,1950000,2100000, 2250000,2400000,2550000,2700000,2850000,3000000,150000,300000,450000,600000,750000,900000,1050000,1200000,1350000,1500000, 1650000,1800000,1950000,2100000,2250000,2400000,2550000,2700000,2850000,3000000}; QVector yData = {.11019,.21956,.32949,.43899,.54803,.65694,.76562,.87487,.98292,1.09146,1.20001,1.30822,1.41599,1.52399, 1.63194,1.73947,1.84646,1.95392,2.06128,2.16844,.11052,.22018,.32939,.43886,.54798,.65739,.76596,.87474, .98300,1.09150,1.20004,1.30818,1.41613,1.52408,1.63159,1.73965,1.84696,1.95445,2.06177,2.16829}; //data source columns Column xDataColumn("x", AbstractColumn::Integer); xDataColumn.replaceInteger(0, xData); Column yDataColumn("y", AbstractColumn::Numeric); yDataColumn.replaceValues(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_polynomial; fitData.degree = 2; XYFitCurve::initFitData(fitData); fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); const int np = fitData.paramNames.size(); QCOMPARE(np, 3); QCOMPARE(fitResult.paramValues.at(0), 0.673565789473684e-3); QCOMPARE(fitResult.errorValues.at(0), 0.107938612033077e-3); QCOMPARE(fitResult.paramValues.at(1), 0.732059160401003e-6); QCOMPARE(fitResult.errorValues.at(1), 0.157817399981659e-9); QCOMPARE(fitResult.paramValues.at(2), -0.316081871345029e-14); QCOMPARE(fitResult.errorValues.at(2), 0.486652849992036e-16); QCOMPARE(fitResult.rsd, 0.205177424076185e-3); QCOMPARE(fitResult.rsquare, 0.999999900178537); QCOMPARE(fitResult.sse, 0.155761768796992e-5); QCOMPARE(fitResult.rms, 0.420977753505385e-7); DEBUG(std::setprecision(15) << fitResult.fdist_F); // result: 370661768.991551 //FuzzyCompare(fitResult.fdist_F, 185330865.995752, 1.); } void FitTest::testLinearNoInt1() { //NIST data for NoInt1 dataset QVector xData = {60,61,62,63,64,65,66,67,68,69,70}; QVector yData = {130,131,132,133,134,135,136,137,138,139,140}; //data source columns Column xDataColumn("x", AbstractColumn::Integer); xDataColumn.replaceInteger(0, xData); Column yDataColumn("y", AbstractColumn::Integer); yDataColumn.replaceInteger(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_custom; XYFitCurve::initFitData(fitData); fitData.model = "c * x"; fitData.paramNames << "c"; const int np = fitData.paramNames.size(); fitData.paramStartValues << 1.; fitData.paramLowerLimits << -std::numeric_limits::max(); fitData.paramUpperLimits << std::numeric_limits::max(); //fitData.eps = 1.e-15; fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); QCOMPARE(np, 1); DEBUG(std::setprecision(15) << fitResult.paramValues.at(0)); // result: 2.07438016513166 FuzzyCompare(fitResult.paramValues.at(0), 2.07438016528926, 1.e-9); DEBUG(std::setprecision(15) << fitResult.errorValues.at(0)); // result: 0.00463315628245255 // QCOMPARE(fitResult.errorValues.at(0), 0.165289256198347e-1); QCOMPARE(fitResult.rsd, 3.56753034006338); // QCOMPARE(fitResult.rsquare, 0.999365492298663); QCOMPARE(fitResult.sse, 127.272727272727); QCOMPARE(fitResult.rms, 12.7272727272727); DEBUG(std::setprecision(15) << fitResult.fdist_F); // result: 370661768.991551 // FuzzyCompare(fitResult.fdist_F, 15750.25, 1.); } void FitTest::testLinearNoInt1_2() { //NIST data for NoInt1 dataset QVector xData = {60,61,62,63,64,65,66,67,68,69,70}; QVector yData = {130,131,132,133,134,135,136,137,138,139,140}; //data source columns Column xDataColumn("x", AbstractColumn::Integer); xDataColumn.replaceInteger(0, xData); Column yDataColumn("y", AbstractColumn::Integer); yDataColumn.replaceInteger(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_polynomial; fitData.degree = 1; XYFitCurve::initFitData(fitData); fitData.paramStartValues[0] = 0; fitData.paramFixed[0] = true; fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); const int np = fitData.paramNames.size(); QCOMPARE(np, 2); QCOMPARE(fitResult.paramValues.at(0), 0.); QCOMPARE(fitResult.paramValues.at(1), 2.07438016528926); DEBUG(std::setprecision(15) << fitResult.errorValues.at(1)); // result: 0.0046331562857966 //QCOMPARE(fitResult.errorValues.at(1), 0.165289256198347e-1); // QCOMPARE(fitResult.rsd, 3.56753034006338); // QCOMPARE(fitResult.rsquare, 0.999365492298663); QCOMPARE(fitResult.sse, 127.272727272727); // QCOMPARE(fitResult.rms, 12.7272727272727); DEBUG(std::setprecision(15) << fitResult.fdist_F); // result: 7.77857142857144 // FuzzyCompare(fitResult.fdist_F, 15750.25, 1.); } void FitTest::testLinearNoInt2() { //NIST data for NoInt2 dataset QVector xData = {4,5,6}; QVector yData = {3,4,4}; //data source columns Column xDataColumn("x", AbstractColumn::Integer); xDataColumn.replaceInteger(0, xData); Column yDataColumn("y", AbstractColumn::Integer); yDataColumn.replaceInteger(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_custom; XYFitCurve::initFitData(fitData); fitData.model = "c * x"; fitData.paramNames << "c"; const int np = fitData.paramNames.size(); fitData.paramStartValues << 1.; fitData.paramLowerLimits << -std::numeric_limits::max(); fitData.paramUpperLimits << std::numeric_limits::max(); //fitData.eps = 1.e-15; fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); QCOMPARE(np, 1); DEBUG(std::setprecision(15) << fitResult.paramValues.at(0)); // result: 0.727272727152573 FuzzyCompare(fitResult.paramValues.at(0), 0.727272727272727, 1.e-9); DEBUG(std::setprecision(15) << fitResult.errorValues.at(0)); // result: 0.0420827316561797 FuzzyCompare(fitResult.errorValues.at(0), 0.420827318078432E-01, 1.e-8); QCOMPARE(fitResult.rsd, 0.369274472937998); // QCOMPARE(fitResult.rsquare, 0.993348115299335); QCOMPARE(fitResult.sse, 0.272727272727273); QCOMPARE(fitResult.rms, 0.136363636363636); DEBUG(std::setprecision(15) << fitResult.fdist_F); // result: 370661768.991551 // FuzzyCompare(fitResult.fdist_F, 298.666666666667, 1.); } void FitTest::testLinearNoInt2_2() { //NIST data for NoInt2 dataset QVector xData = {4,5,6}; QVector yData = {3,4,4}; //data source columns Column xDataColumn("x", AbstractColumn::Integer); xDataColumn.replaceInteger(0, xData); Column yDataColumn("y", AbstractColumn::Integer); yDataColumn.replaceInteger(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_polynomial; fitData.degree = 1; XYFitCurve::initFitData(fitData); fitData.paramStartValues[0] = 0; fitData.paramFixed[0] = true; fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); const int np = fitData.paramNames.size(); QCOMPARE(np, 2); QCOMPARE(fitResult.paramValues.at(0), 0.); QCOMPARE(fitResult.paramValues.at(1), 0.727272727272727); DEBUG(std::setprecision(15) << fitResult.errorValues.at(1)); // result: 0.0595139700643615 // QCOMPARE(fitResult.errorValues.at(1), 0.420827318078432e-1); // QCOMPARE(fitResult.rsd, 0.369274472937998); // result: 0.522233 // QCOMPARE(fitResult.rsquare, 0.993348115299335); // result: 0.590909 QCOMPARE(fitResult.sse, 0.272727272727273); // QCOMPARE(fitResult.rms, 0.136363636363636); // result: 0.272727 DEBUG(std::setprecision(15) << fitResult.fdist_F); // result: 2.44444444444445 // FuzzyCompare(fitResult.fdist_F, 298.666666666667, 1.); } void FitTest::testLinearFilip() { //NIST data for Filip dataset QVector xData = {-6.860120914,-4.324130045,-4.358625055,-4.358426747,-6.955852379,-6.661145254,-6.355462942,-6.118102026, -7.115148017,-6.815308569,-6.519993057,-6.204119983,-5.853871964,-6.109523091,-5.79832982,-5.482672118,-5.171791386,-4.851705903, -4.517126416,-4.143573228,-3.709075441,-3.499489089,-6.300769497,-5.953504836,-5.642065153,-5.031376979,-4.680685696,-4.329846955, -3.928486195,-8.56735134,-8.363211311,-8.107682739,-7.823908741,-7.522878745,-7.218819279,-6.920818754,-6.628932138,-6.323946875, -5.991399828,-8.781464495,-8.663140179,-8.473531488,-8.247337057,-7.971428747,-7.676129393,-7.352812702,-7.072065318,-6.774174009, -6.478861916,-6.159517513,-6.835647144,-6.53165267,-6.224098421,-5.910094889,-5.598599459,-5.290645224,-4.974284616,-4.64454848, -4.290560426,-3.885055584,-3.408378962,-3.13200249,-8.726767166,-8.66695597,-8.511026475,-8.165388579,-7.886056648,-7.588043762, -7.283412422,-6.995678626,-6.691862621,-6.392544977,-6.067374056,-6.684029655,-6.378719832,-6.065855188,-5.752272167,-5.132414673, -4.811352704,-4.098269308,-3.66174277,-3.2644011}; QVector yData = {0.8116,0.9072,0.9052,0.9039,0.8053,0.8377,0.8667,0.8809,0.7975,0.8162,0.8515,0.8766,0.8885,0.8859,0.8959,0.8913, 0.8959,0.8971,0.9021,0.909,0.9139,0.9199,0.8692,0.8872,0.89,0.891,0.8977,0.9035,0.9078,0.7675,0.7705,0.7713,0.7736,0.7775,0.7841, 0.7971,0.8329,0.8641,0.8804,0.7668,0.7633,0.7678,0.7697,0.77,0.7749,0.7796,0.7897,0.8131,0.8498,0.8741,0.8061,0.846,0.8751,0.8856, 0.8919,0.8934,0.894,0.8957,0.9047,0.9129,0.9209,0.9219,0.7739,0.7681,0.7665,0.7703,0.7702,0.7761,0.7809,0.7961,0.8253,0.8602, 0.8809,0.8301,0.8664,0.8834,0.8898,0.8964,0.8963,0.9074,0.9119,0.9228}; //data source columns Column xDataColumn("x", AbstractColumn::Numeric); xDataColumn.replaceValues(0, xData); Column yDataColumn("y", AbstractColumn::Numeric); yDataColumn.replaceValues(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_polynomial; fitData.degree = 10; XYFitCurve::initFitData(fitData); const int np = fitData.paramNames.size(); fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); QCOMPARE(np, 11); DEBUG(std::setprecision(15) << fitResult.paramValues.at(0)); // result: -1467.48962615175 FuzzyCompare(fitResult.paramValues.at(0), -1467.48961422980, 1.e-7); DEBUG(std::setprecision(15) << fitResult.errorValues.at(0)); // result: 298.084524514884 FuzzyCompare(fitResult.errorValues.at(0), 298.084530995537, 1.e-7); DEBUG(std::setprecision(15) << fitResult.paramValues.at(1)); // result: -2772.1796150428 FuzzyCompare(fitResult.paramValues.at(1), -2772.17959193342, 1.e-7); DEBUG(std::setprecision(15) << fitResult.errorValues.at(1)); // result: 559.779853249694 FuzzyCompare(fitResult.errorValues.at(1), 559.779865474950, 1.e-7); DEBUG(std::setprecision(15) << fitResult.paramValues.at(2)); // result: -2316.37110148409 FuzzyCompare(fitResult.paramValues.at(2), -2316.37108160893, 1.e-7); DEBUG(std::setprecision(15) << fitResult.errorValues.at(2)); // result: 466.477561928144 FuzzyCompare(fitResult.errorValues.at(2), 466.477572127796, 1.e-7); DEBUG(std::setprecision(15) << fitResult.paramValues.at(3)); // result: -1127.97395097195 FuzzyCompare(fitResult.paramValues.at(3), -1127.97394098372, 1.e-7); DEBUG(std::setprecision(15) << fitResult.errorValues.at(3)); // result: 227.204269523115 FuzzyCompare(fitResult.errorValues.at(3), 227.204274477751, 1.e-7); DEBUG(std::setprecision(15) << fitResult.paramValues.at(4)); // result: -354.478236951913 FuzzyCompare(fitResult.paramValues.at(4), -354.478233703349, 1.e-7); DEBUG(std::setprecision(15) << fitResult.errorValues.at(4)); // result: 71.6478645361214 FuzzyCompare(fitResult.errorValues.at(4), 71.6478660875927, 1.e-7); DEBUG(std::setprecision(15) << fitResult.paramValues.at(5)); // result: -75.1242024539908 FuzzyCompare(fitResult.paramValues.at(5), -75.1242017393757, 1.e-7); DEBUG(std::setprecision(15) << fitResult.errorValues.at(5)); // result: 15.289717547564 FuzzyCompare(fitResult.errorValues.at(5), 15.2897178747400, 1.e-7); DEBUG(std::setprecision(15) << fitResult.paramValues.at(6)); // result: -10.875318143236 FuzzyCompare(fitResult.paramValues.at(6), -10.8753180355343, 1.e-7); DEBUG(std::setprecision(15) << fitResult.errorValues.at(6)); // result: 2.23691155110776 FuzzyCompare(fitResult.errorValues.at(6), 2.23691159816033, 1.e-7); DEBUG(std::setprecision(15) << fitResult.paramValues.at(7)); // result: -1.06221499687347 FuzzyCompare(fitResult.paramValues.at(7), -1.06221498588947, 1.e-7); DEBUG(std::setprecision(15) << fitResult.errorValues.at(7)); // result: 0.221624317377432 FuzzyCompare(fitResult.errorValues.at(7), 0.221624321934227, 1.e-7); DEBUG(std::setprecision(15) << fitResult.paramValues.at(8)); // result: -0.0670191161850038 FuzzyCompare(fitResult.paramValues.at(8), -0.670191154593408E-01, 1.e-7); DEBUG(std::setprecision(15) << fitResult.errorValues.at(8)); // result: 0.0142363760310402 FuzzyCompare(fitResult.errorValues.at(8), 0.142363763154724E-01, 1.e-7); DEBUG(std::setprecision(15) << fitResult.paramValues.at(9)); // result: -0.00246781081080665 FuzzyCompare(fitResult.paramValues.at(9), -0.246781078275479E-02, 1.e-7); DEBUG(std::setprecision(15) << fitResult.errorValues.at(9)); // result: 0.000535617398555022 FuzzyCompare(fitResult.errorValues.at(9), 0.535617408889821E-03, 1.e-7); DEBUG(std::setprecision(15) << fitResult.paramValues.at(10)); // result: -4.02962529900222e-05 FuzzyCompare(fitResult.paramValues.at(10), -0.402962525080404E-04, 1.e-7); DEBUG(std::setprecision(15) << fitResult.errorValues.at(10)); // result: 8.96632820770946e-06 FuzzyCompare(fitResult.errorValues.at(10), 0.896632837373868E-05, 1.e-7); DEBUG(std::setprecision(15) << fitResult.rsd); // result: 0.00334801050105949 FuzzyCompare(fitResult.rsd, 0.334801051324544E-02, 1.e-8); DEBUG(std::setprecision(15) << fitResult.rsquare); // result: 0.996727416209443 FuzzyCompare(fitResult.rsquare, 0.996727416185620, 1.e-9); DEBUG(std::setprecision(15) << fitResult.sse); // result: 0.00079585137637953 FuzzyCompare(fitResult.sse, 0.795851382172941E-03, 1.e-7); DEBUG(std::setprecision(15) << fitResult.rms); // result: 1.12091743152047e-05 FuzzyCompare(fitResult.rms, 0.112091743968020E-04, 1.e-7); DEBUG(std::setprecision(15) << fitResult.fdist_F); // result: 21695.3956090808 // FuzzyCompare(fitResult.fdist_F, 2162.43954511489, 1.e-9); } void FitTest::testLinearWampler1() { //NIST data for Wampler1 dataset QVector xData = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20}; QVector yData = {1,6,63,364,1365,3906,9331,19608,37449,66430,111111, 177156,271453,402234,579195,813616,1118481,1508598,2000719,2613660,3368421}; //data source columns Column xDataColumn("x", AbstractColumn::Integer); xDataColumn.replaceInteger(0, xData); Column yDataColumn("y", AbstractColumn::Integer); yDataColumn.replaceInteger(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_polynomial; fitData.degree = 5; XYFitCurve::initFitData(fitData); fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); const int np = fitData.paramNames.size(); QCOMPARE(np, 6); for (int i = 0; i < np; i++) { const double paramValue = fitResult.paramValues.at(i); const double errorValue = fitResult.errorValues.at(i); QCOMPARE(paramValue, 1.0); QCOMPARE(errorValue, 0.0); } QCOMPARE(fitResult.rsd, 0.0); QCOMPARE(fitResult.rsquare, 1.0); QCOMPARE(fitResult.sse, 0.0); QCOMPARE(fitResult.rms, 0.0); } void FitTest::testLinearWampler2() { //NIST data for Wampler2 dataset QVector xData = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20}; QVector yData = {1.00000,1.11111,1.24992,1.42753,1.65984,1.96875,2.38336,2.94117,3.68928,4.68559, 6.00000,7.71561,9.92992,12.75603,16.32384,20.78125,26.29536,33.05367,41.26528,51.16209,63.00000}; //data source columns Column xDataColumn("x", AbstractColumn::Integer); xDataColumn.replaceInteger(0, xData); Column yDataColumn("y", AbstractColumn::Numeric); yDataColumn.replaceValues(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_polynomial; fitData.degree = 5; XYFitCurve::initFitData(fitData); fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results qDebug() << "STATUS " << fitResult.status; QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); const int np = fitData.paramNames.size(); QCOMPARE(np, 6); QCOMPARE(fitResult.paramValues.at(0), 1.0); QCOMPARE(fitResult.paramValues.at(1), 0.1); QCOMPARE(fitResult.paramValues.at(2), 0.01); QCOMPARE(fitResult.paramValues.at(3), 0.001); QCOMPARE(fitResult.paramValues.at(4), 0.0001); QCOMPARE(fitResult.paramValues.at(5), 0.00001); for (int i = 0; i < np; i++) { const double errorValue = fitResult.errorValues.at(i); FuzzyCompare(errorValue, 0., 1.); } DEBUG(std::setprecision(15) << fitResult.rsd); // result: 2.32458538254974e-15 FuzzyCompare(fitResult.rsd, 0., 1.); QCOMPARE(fitResult.rsquare, 1.); DEBUG(std::setprecision(15) << fitResult.sse); // result: 8.1055458011459e-29 FuzzyCompare(fitResult.sse, 0., 1.); DEBUG(std::setprecision(15) << fitResult.rms); // result: 5.40369720076393e-30 FuzzyCompare(fitResult.rms, 0., 1.); DEBUG(std::setprecision(15) << fitResult.fdist_F); // result: 1.22192608844148e+33 // FuzzyCompare(fitResult.fdist_F, INF, 1.); } void FitTest::testLinearWampler3() { //NIST data for Wampler3 dataset QVector xData = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20}; QVector yData = {760.,-2042.,2111.,-1684.,3888.,1858.,11379.,17560.,39287.,64382.,113159., 175108.,273291.,400186.,581243.,811568.,1121004.,1506550.,2002767.,2611612.,3369180.}; //data source columns Column xDataColumn("x", AbstractColumn::Integer); xDataColumn.replaceInteger(0, xData); Column yDataColumn("y", AbstractColumn::Numeric); yDataColumn.replaceValues(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_polynomial; fitData.degree = 5; XYFitCurve::initFitData(fitData); fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); const int np = fitData.paramNames.size(); QCOMPARE(np, 6); for (int i = 0; i < np; i++) { const double paramValue = fitResult.paramValues.at(i); QCOMPARE(paramValue, 1.0); } // TODO: actual errors are smaller // QCOMPARE(fitResult.errorValues.at(0), 2152.32624678170); // QCOMPARE(fitResult.errorValues.at(1), 2363.55173469681); // QCOMPARE(fitResult.errorValues.at(2), 779.343524331583); // QCOMPARE(fitResult.errorValues.at(3), 101.475507550350); // QCOMPARE(fitResult.errorValues.at(4), 5.64566512170752); // QCOMPARE(fitResult.errorValues.at(5), 0.112324854679312); QCOMPARE(fitResult.rsd, 2360.14502379268); QCOMPARE(fitResult.rsquare, 0.999995559025820); QCOMPARE(fitResult.sse, 83554268.0000000); QCOMPARE(fitResult.rms, 5570284.53333333); // QCOMPARE(fitResult.fdist_F, 675524.458240122); } void FitTest::testLinearWampler4() { //NIST data for Wampler4 dataset QVector xData = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20}; QVector yData = {75901,-204794,204863,-204436,253665,-200894,214131,-185192,221249,-138370, 315911,-27644,455253,197434,783995,608816,1370781,1303798,2205519,2408860,3444321}; //data source columns Column xDataColumn("x", AbstractColumn::Integer); xDataColumn.replaceInteger(0, xData); Column yDataColumn("y", AbstractColumn::Integer); yDataColumn.replaceInteger(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_polynomial; fitData.degree = 5; XYFitCurve::initFitData(fitData); fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); const int np = fitData.paramNames.size(); QCOMPARE(np, 6); for (int i = 0; i < np; i++) { const double paramValue = fitResult.paramValues.at(i); QCOMPARE(paramValue, 1.0); } // TODO: actual errors are smaller // QCOMPARE(fitResult.errorValues.at(0), 215232.624678170); // QCOMPARE(fitResult.errorValues.at(1), 236355.173469681); // QCOMPARE(fitResult.errorValues.at(2), 77934.3524331583); // QCOMPARE(fitResult.errorValues.at(3), 10147.5507550350); // QCOMPARE(fitResult.errorValues.at(4), 564.566512170752); // QCOMPARE(fitResult.errorValues.at(5), 11.2324854679312); QCOMPARE(fitResult.rsd, 236014.502379268); QCOMPARE(fitResult.rsquare, 0.957478440825662); QCOMPARE(fitResult.sse, 835542680000.000); QCOMPARE(fitResult.rms, 55702845333.3333); // QCOMPARE(fitResult.fdist_F, 67.5524458240122); } void FitTest::testLinearWampler5() { //NIST data for Wampler5 dataset QVector xData = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20}; QVector yData = {7590001,-20479994,20480063,-20479636,25231365,-20476094,20489331,-20460392,18417449,-20413570, 20591111,-20302844,18651453,-20077766,21059195,-19666384,26348481,-18971402,22480719,-17866340,10958421}; //data source columns Column xDataColumn("x", AbstractColumn::Integer); xDataColumn.replaceInteger(0, xData); Column yDataColumn("y", AbstractColumn::Integer); yDataColumn.replaceInteger(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = nsl_fit_model_polynomial; fitData.degree = 5; XYFitCurve::initFitData(fitData); fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); const int np = fitData.paramNames.size(); QCOMPARE(np, 6); for (int i = 0; i < np; i++) { const double paramValue = fitResult.paramValues.at(i); QCOMPARE(paramValue, 1.0); } // TODO: actual errors are smaller // QCOMPARE(fitResult.errorValues.at(0), 21523262.4678170); // QCOMPARE(fitResult.errorValues.at(1), 23635517.3469681); // QCOMPARE(fitResult.errorValues.at(2), 7793435.24331583); // QCOMPARE(fitResult.errorValues.at(3), 1014755.07550350); // QCOMPARE(fitResult.errorValues.at(4), 56456.6512170752); // QCOMPARE(fitResult.errorValues.at(5), 1123.24854679312); QCOMPARE(fitResult.rsd, 23601450.2379268); QCOMPARE(fitResult.rsquare, 0.224668921574940E-02); QCOMPARE(fitResult.sse, 0.835542680000000E+16); QCOMPARE(fitResult.rms, 557028453333333.); // QCOMPARE(fitResult.fdist_F, 0.675524458240122E-02); } //############################################################################## //############# non-linear regression with NIST datasets ##################### //############################################################################## void FitTest::testNonLinearMisra1a() { //NIST data for Misra1a dataset QVector xData = {77.6E0,114.9E0,141.1E0,190.8E0,239.9E0,289.0E0,332.8E0,378.4E0,434.8E0,477.3E0,536.8E0,593.1E0,689.1E0,760.0E0}; QVector yData = {10.07E0,14.73E0,17.94E0,23.93E0,29.61E0,35.18E0,40.02E0,44.82E0,50.76E0,55.05E0,61.01E0,66.40E0,75.47E0,81.78E0}; //data source columns Column xDataColumn("x", AbstractColumn::Numeric); xDataColumn.replaceValues(0, xData); Column yDataColumn("y", AbstractColumn::Numeric); yDataColumn.replaceValues(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_custom; XYFitCurve::initFitData(fitData); fitData.model = "b1*(1.-exp(-b2*x))"; fitData.paramNames << "b1" << "b2"; fitData.eps = 1.e-12; const int np = fitData.paramNames.size(); fitData.paramStartValues << 500. << 0.0001; fitData.paramLowerLimits << -std::numeric_limits::max() << -std::numeric_limits::max(); fitData.paramUpperLimits << std::numeric_limits::max() << std::numeric_limits::max(); fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); QCOMPARE(np, 2); DEBUG(std::setprecision(15) << fitResult.paramValues.at(0)); // result: 238.942305251573 FuzzyCompare(fitResult.paramValues.at(0), 2.3894212918E+02, 1.e-6); DEBUG(std::setprecision(15) << fitResult.errorValues.at(0)); // result: 2.70651225218243 FuzzyCompare(fitResult.errorValues.at(0), 2.7070075241E+00, 1.e-3); DEBUG(std::setprecision(15) << fitResult.paramValues.at(1)); // result: 0.000550155958419367 FuzzyCompare(fitResult.paramValues.at(1), 5.5015643181E-04, 1.e-6); DEBUG(std::setprecision(15) << fitResult.errorValues.at(1)); // result: 7.26565480949189e-06 FuzzyCompare(fitResult.errorValues.at(1), 7.2668688436E-06, 1.e-3); DEBUG(std::setprecision(15) << fitResult.rsd); // result: 0.101878763320394 FuzzyCompare(fitResult.rsd, 1.0187876330E-01, 1.e-9); DEBUG(std::setprecision(15) << fitResult.sse); // result: 0.124551388988316 FuzzyCompare(fitResult.sse, 1.2455138894E-01, 1.e-9); } void FitTest::testNonLinearMisra1b() { //NIST data for Misra1b dataset QVector xData = {77.6E0,114.9E0,141.1E0,190.8E0,239.9E0,289.0E0,332.8E0,378.4E0,434.8E0,477.3E0,536.8E0,593.1E0,689.1E0,760.0E0}; QVector yData = {10.07E0,14.73E0,17.94E0,23.93E0,29.61E0,35.18E0,40.02E0,44.82E0,50.76E0,55.05E0,61.01E0,66.40E0,75.47E0,81.78E0}; //data source columns Column xDataColumn("x", AbstractColumn::Numeric); xDataColumn.replaceValues(0, xData); Column yDataColumn("y", AbstractColumn::Numeric); yDataColumn.replaceValues(0, yData); XYFitCurve fitCurve("fit"); fitCurve.setXDataColumn(&xDataColumn); fitCurve.setYDataColumn(&yDataColumn); //prepare the fit XYFitCurve::FitData fitData = fitCurve.fitData(); fitData.modelCategory = nsl_fit_model_custom; XYFitCurve::initFitData(fitData); fitData.model = "b1*(1.-1./(1.+b2*x/2)^2)"; fitData.paramNames << "b1" << "b2"; fitData.eps = 1.e-12; const int np = fitData.paramNames.size(); fitData.paramStartValues << 500. << 0.0001; fitData.paramLowerLimits << -std::numeric_limits::max() << -std::numeric_limits::max(); fitData.paramUpperLimits << std::numeric_limits::max() << std::numeric_limits::max(); fitCurve.setFitData(fitData); //perform the fit fitCurve.recalculate(); const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); //check the results QCOMPARE(fitResult.available, true); QCOMPARE(fitResult.valid, true); QCOMPARE(np, 2); DEBUG(std::setprecision(15) << fitResult.paramValues.at(0)); // result: 337.99775062098 FuzzyCompare(fitResult.paramValues.at(0), 3.3799746163E+02, 1.e-6); DEBUG(std::setprecision(15) << fitResult.errorValues.at(0)); // result: 3.16358581006192 FuzzyCompare(fitResult.errorValues.at(0), 3.1643950207E+00, 1.e-3); DEBUG(std::setprecision(15) << fitResult.paramValues.at(1)); // result: 0.000390390523934039 FuzzyCompare(fitResult.paramValues.at(1), 3.9039091287E-04, 1.e-6); DEBUG(std::setprecision(15) << fitResult.errorValues.at(1)); // result: 4.25373670682006e-06 FuzzyCompare(fitResult.errorValues.at(1), 4.2547321834E-06, 1.e-3); DEBUG(std::setprecision(15) << fitResult.rsd); // result: 0.0793014720259488 FuzzyCompare(fitResult.rsd, 7.9301471998E-02, 1.e-9); DEBUG(std::setprecision(15) << fitResult.sse); // result: 0.0754646815857881 FuzzyCompare(fitResult.sse, 7.5464681533E-02, 1.e-9); } +void FitTest::testNonLinearMisra1c() { + //NIST data for Misra1b dataset + QVector xData = {77.6E0,114.9E0,141.1E0,190.8E0,239.9E0,289.0E0,332.8E0,378.4E0,434.8E0,477.3E0,536.8E0,593.1E0,689.1E0,760.0E0}; + QVector yData = {10.07E0,14.73E0,17.94E0,23.93E0,29.61E0,35.18E0,40.02E0,44.82E0,50.76E0,55.05E0,61.01E0,66.40E0,75.47E0,81.78E0}; + + //data source columns + Column xDataColumn("x", AbstractColumn::Numeric); + xDataColumn.replaceValues(0, xData); + + Column yDataColumn("y", AbstractColumn::Numeric); + yDataColumn.replaceValues(0, yData); + + XYFitCurve fitCurve("fit"); + fitCurve.setXDataColumn(&xDataColumn); + fitCurve.setYDataColumn(&yDataColumn); + + //prepare the fit + XYFitCurve::FitData fitData = fitCurve.fitData(); + fitData.modelCategory = nsl_fit_model_custom; + XYFitCurve::initFitData(fitData); + fitData.model = "b1*(1.-1./sqrt(1.+2.*b2*x))"; + fitData.paramNames << "b1" << "b2"; + fitData.eps = 1.e-12; + const int np = fitData.paramNames.size(); + fitData.paramStartValues << 500. << 0.0001; + fitData.paramLowerLimits << -std::numeric_limits::max() << -std::numeric_limits::max(); + fitData.paramUpperLimits << std::numeric_limits::max() << std::numeric_limits::max(); + fitCurve.setFitData(fitData); + + //perform the fit + fitCurve.recalculate(); + const XYFitCurve::FitResult& fitResult = fitCurve.fitResult(); + + //check the results + QCOMPARE(fitResult.available, true); + QCOMPARE(fitResult.valid, true); + + QCOMPARE(np, 2); + + DEBUG(std::setprecision(15) << fitResult.paramValues.at(0)); // result: 636.427904767969 + FuzzyCompare(fitResult.paramValues.at(0), 6.3642725809E+02, 1.e-5); + DEBUG(std::setprecision(15) << fitResult.errorValues.at(0)); // result: 4.66168062054875 + FuzzyCompare(fitResult.errorValues.at(0), 4.6638326572E+00, 1.e-3); + DEBUG(std::setprecision(15) << fitResult.paramValues.at(1)); // result: 0.000208136026420746 + FuzzyCompare(fitResult.paramValues.at(1), 2.0813627256E-04, 1.e-5); + DEBUG(std::setprecision(15) << fitResult.errorValues.at(1)); // result: 1.77209416674174e-06 + FuzzyCompare(fitResult.errorValues.at(1), 1.7728423155E-06, 1.e-3); + + DEBUG(std::setprecision(15) << fitResult.rsd); // result: 0.0584286153041661 + FuzzyCompare(fitResult.rsd, 5.8428615257E-02, 1.e-9); + DEBUG(std::setprecision(15) << fitResult.sse); // result: 0.0409668370363468 + FuzzyCompare(fitResult.sse, 4.0966836971E-02, 1.e-8); +} + //TODO //############################################################################## //######################### Fits with weights ################################# //############################################################################## //TODO QTEST_MAIN(FitTest) diff --git a/tests/analysis/fit/FitTest.h b/tests/analysis/fit/FitTest.h index 847fad0f9..c1e60fcd3 100644 --- a/tests/analysis/fit/FitTest.h +++ b/tests/analysis/fit/FitTest.h @@ -1,62 +1,63 @@ /*************************************************************************** File : FitTest.h Project : LabPlot Description : Tests for data fitting -------------------------------------------------------------------- Copyright : (C) 2017 Alexander Semke (alexander.semke@web.de) Copyright : (C) 2018 Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #include class FitTest : public QObject { Q_OBJECT private slots: void initTestCase(); // compare floats with given delta (could be useful for other tests too) // delta - relative error (set to 1. if expected == 0.) static inline void FuzzyCompare(double actual, double expected, double delta = 1.e-12) { qDebug() << qSetRealNumberPrecision(15) << actual - fabs(actual)*delta << "<=" << expected << "<=" << actual + fabs(actual)*delta; QVERIFY(actual - fabs(actual)*delta <= expected && actual + fabs(actual)*delta >=expected); } //linear regression (see NIST/linear data) void testLinearNorris(); void testLinearPontius(); void testLinearNoInt1(); // using custom model void testLinearNoInt1_2(); // using polynomial model with fixed parameter void testLinearNoInt2(); // using custom model void testLinearNoInt2_2(); // using polynomial model with fixed parameter void testLinearFilip(); void testLinearWampler1(); void testLinearWampler2(); void testLinearWampler3(); void testLinearWampler4(); void testLinearWampler5(); //non-linear regression void testNonLinearMisra1a(); void testNonLinearMisra1b(); + void testNonLinearMisra1c(); //fits with weights };