diff --git a/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp b/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp index 8f03d79c7..835d5a03f 100644 --- a/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp +++ b/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp @@ -1,2271 +1,2271 @@  /*************************************************************************** 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; if (modelCategory != nsl_fit_model_custom) { DEBUG("XYFitCurve::initFitData() for model category = " << nsl_fit_model_category_name[modelCategory] << ", model type = " << modelType << ", degree = " << degree); paramNames.clear(); } else { DEBUG("XYFitCurve::initFitData() for model category = nsl_fit_model_custom, model type = " << modelType << ", degree = " << degree); } 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]; 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: // yerror are sigmas 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: // yerror are weights for(int i = 0; i < (int)n; i++) if (i < yerrorVector.size()) weight[i] = yerror[i]; break; case nsl_fit_weight_inverse: // yerror are inverse weights 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" unsigned int nf = 0; // number of fixed parameter for (unsigned int i = 0; i < np; i++) { const bool fixed = fitData.paramFixed.data()[i]; if (fixed) nf++; DEBUG("parameter " << i << " fixed: " << fixed); } //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("Turning 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 (TODO: check correct residuals) 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"); unsigned int iter2 = 0; double chi = 0, chiOld = 0; double *fun = new double[n]; do { iter2++; chiOld = chi; //printf("iter2 = %d\n", iter2); // calculate df[i] for (size_t i = 0; i < n; i++) fun[i] = gsl_vector_get(s->f, i) * 1./sqrt(weight[i]) + ydata[i]; // calculate weight[i] for (size_t i = 0; i < n; i++) { // calculate df[i] - int index = i-1; + size_t index = i-1; if (i == 0) index = i; if (i == n-1) index = i-2; double df = (fun[index+1] - fun[index])/(xdata[index+1] - xdata[index]); //printf("df = %g\n", df); double sigmasq = 1.; switch (fitData.xWeightsType) { // x-error type: f'(x)^2*s_x^2 = f'(x)/w_x case nsl_fit_weight_no: break; case nsl_fit_weight_direct: // xerror = w_x sigmasq = df*df/xerror[i]; break; case nsl_fit_weight_instrumental: // xerror = s_x sigmasq = df*df*xerror[i]*xerror[i]; break; case nsl_fit_weight_inverse: // xerror = 1/w_x = s_x^2 sigmasq = df*df*xerror[i]; break; case nsl_fit_weight_statistical: // s_x^2 = 1/w_x = x sigmasq = xdata[i]; break; case nsl_fit_weight_relative: // s_x^2 = 1/w_x = x^2 sigmasq = xdata[i]*xdata[i]; break; case nsl_fit_weight_statistical_fit: // unused case nsl_fit_weight_relative_fit: break; } if (yerrorVector.size() > 0) { switch (fitData.yWeightsType) { // y-error types: s_y^2 = 1/w_y case nsl_fit_weight_no: break; case nsl_fit_weight_direct: // yerror = w_y sigmasq += 1./yerror[i]; break; case nsl_fit_weight_instrumental: // yerror = s_y sigmasq += yerror[i]*yerror[i]; break; case nsl_fit_weight_inverse: // yerror = 1/w_y sigmasq += yerror[i]; break; case nsl_fit_weight_statistical: // unused case nsl_fit_weight_relative: break; case nsl_fit_weight_statistical_fit: // s_y^2 = 1/w_y = Y_i sigmasq += fun[i]; break; case nsl_fit_weight_relative_fit: // s_y^2 = 1/w_y = Y_i^2 sigmasq += fun[i]*fun[i]; break; } } //printf ("sigma[%d] = %g\n", i, sqrt(sigmasq)); weight[i] = 1./sigmasq; } // update weights gsl_multifit_fdfsolver_set(s, &f, &x.vector); do { // fit iter++; writeSolverState(s); status = gsl_multifit_fdfsolver_iterate (s); //printf ("status = %s\n", gsl_strerror (status)); 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); chi = gsl_blas_dnrm2(s->f); //printf("chi = %.12g (dchi = %g)\n", chi, fabs(chi-chiOld)); } while (iter2 < maxIters && fabs(chi-chiOld) > fitData.eps); delete[] fun; } delete[] weight; // unscale start parameter 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 - nf); // samples - (parameter - fixed parameter) //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; // SST needed for coefficient of determination, R-squared fitResult.sst = gsl_stats_tss(ydata, 1, n); // for a linear model without intercept R-squared is calculated different // see https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-does-summary_0028_0029-report-strange-results-for-the-R_005e2-estimate-when-I-fit-a-linear-model-with-no-intercept_003f if (fitData.modelCategory == nsl_fit_model_basic && fitData.modelType == nsl_fit_model_polynomial && fitData.degree == 1 && x_init[0] == 0) { DEBUG("Using alternative R^2 for linear model without intercept"); fitResult.sst = gsl_stats_tss_m(ydata, 1, n, 0); } if (fitResult.sst < fitResult.sse) { DEBUG("Using alternative R^2 since R^2 would be negative (probably custom model without intercept)"); fitResult.sst = gsl_stats_tss_m(ydata, 1, n, 0); } fitResult.rsquare = nsl_stats_rsquare(fitResult.sse, fitResult.sst); fitResult.rsquareAdj = nsl_stats_rsquareAdj(fitResult.rsquare, np, fitResult.dof, 1); fitResult.chisq_p = nsl_stats_chisq_p(fitResult.sse, fitResult.dof); fitResult.fdist_F = nsl_stats_fdist_F(fitResult.sst, fitResult.rms, np, 1); fitResult.fdist_p = nsl_stats_fdist_p(fitResult.fdist_F, np, fitResult.dof); fitResult.logLik = nsl_stats_logLik(fitResult.sse, n); fitResult.aic = nsl_stats_aic(fitResult.sse, n, np, 1); fitResult.bic = nsl_stats_bic(fitResult.sse, n, np, 1); //parameter values // GSL: const double c = GSL_MAX_DBL(1., sqrt(fitResult.rms)); // increase error for poor fit // NIST: const double c = sqrt(fitResult.rms); // increase error for poor fit, decrease for good fit const double c = sqrt(fitResult.rms); 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("xWeightsType", QString::number(d->fitData.xWeightsType)); 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("xWeightsType", fitData.xWeightsType, nsl_fit_weight_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; }