File indexing completed on 2025-02-23 03:35:12

0001 /*
0002     File                 : XYFitCurve.cpp
0003     Project              : LabPlot
0004     Description          : A xy-curve defined by a fit model
0005     --------------------------------------------------------------------
0006     SPDX-FileCopyrightText: 2014-2021 Alexander Semke <alexander.semke@web.de>
0007     SPDX-FileCopyrightText: 2016-2022 Stefan Gerlach <stefan.gerlach@uni.kn>
0008 
0009     SPDX-License-Identifier: GPL-2.0-or-later
0010 */
0011 
0012 /*!
0013   \class XYFitCurve
0014   \brief A xy-curve defined by a fit model
0015 
0016   \ingroup worksheet
0017 */
0018 
0019 #include "XYFitCurve.h"
0020 #include "XYFitCurvePrivate.h"
0021 #include "backend/core/AbstractColumn.h"
0022 #include "backend/core/column/Column.h"
0023 #include "backend/gsl/ExpressionParser.h"
0024 #include "backend/gsl/errors.h"
0025 #include "backend/lib/XmlStreamReader.h"
0026 #include "backend/lib/commandtemplates.h"
0027 #include "backend/lib/macros.h"
0028 #include "backend/worksheet/plots/cartesian/Histogram.h"
0029 
0030 #include "backend/gsl/parser.h"
0031 #include "backend/nsl/nsl_sf_stats.h"
0032 #include "backend/nsl/nsl_stats.h"
0033 #include <gsl/gsl_blas.h>
0034 #include <gsl/gsl_cdf.h>
0035 #include <gsl/gsl_matrix.h>
0036 #include <gsl/gsl_multifit.h>
0037 #include <gsl/gsl_sf_erf.h>
0038 #include <gsl/gsl_sf_gamma.h>
0039 #include <gsl/gsl_statistics.h>
0040 #include <gsl/gsl_vector.h>
0041 #include <gsl/gsl_version.h>
0042 
0043 #include <QDateTime>
0044 #include <QElapsedTimer>
0045 #include <QIcon>
0046 #include <QThreadPool>
0047 
0048 XYFitCurve::XYFitCurve(const QString& name)
0049     : XYAnalysisCurve(name, new XYFitCurvePrivate(this), AspectType::XYFitCurve) {
0050 }
0051 
0052 XYFitCurve::XYFitCurve(const QString& name, XYFitCurvePrivate* dd)
0053     : XYAnalysisCurve(name, dd, AspectType::XYFitCurve) {
0054 }
0055 
0056 // no need to delete the d-pointer here - it inherits from QGraphicsItem
0057 // and is deleted during the cleanup in QGraphicsScene
0058 XYFitCurve::~XYFitCurve() = default;
0059 
0060 void XYFitCurve::recalculate() {
0061     Q_D(XYFitCurve);
0062     d->recalculate();
0063 }
0064 
0065 void XYFitCurve::evaluate(bool preview) {
0066     Q_D(XYFitCurve);
0067     if (d->evaluate(preview)) {
0068         // redraw the curve
0069         recalcLogicalPoints();
0070         Q_EMIT dataChanged();
0071     }
0072 }
0073 
0074 const XYAnalysisCurve::Result& XYFitCurve::result() const {
0075     Q_D(const XYFitCurve);
0076     return d->fitResult;
0077 }
0078 void XYFitCurve::initStartValues(const XYCurve* curve) {
0079     Q_D(XYFitCurve);
0080     XYFitCurve::FitData& fitData = d->fitData;
0081     initStartValues(fitData, curve);
0082 }
0083 
0084 void XYFitCurve::initStartValues(XYFitCurve::FitData& fitData, const XYCurve* curve) {
0085     DEBUG(Q_FUNC_INFO);
0086     // TODO: curve used for anything?
0087     if (!curve) {
0088         DEBUG(Q_FUNC_INFO << ", WARNING: no curve given");
0089         return;
0090     }
0091 
0092     Q_D(XYFitCurve);
0093     const Column* xColumn = dynamic_cast<const Column*>(d->xDataColumn);
0094     const Column* yColumn = dynamic_cast<const Column*>(d->yDataColumn);
0095     const Column* yErrorColumn = dynamic_cast<const Column*>(d->yErrorColumn);
0096 
0097     if (!xColumn || !yColumn) {
0098         DEBUG(Q_FUNC_INFO << ", data columns not available");
0099         return;
0100     }
0101 
0102     DEBUG(Q_FUNC_INFO << ", x data rows = " << xColumn->rowCount());
0103     DEBUG(Q_FUNC_INFO << ", y data rows = " << yColumn->rowCount());
0104 
0105     nsl_fit_model_category modelCategory = fitData.modelCategory;
0106     int modelType = fitData.modelType;
0107     const int degree = fitData.degree;
0108     DEBUG(Q_FUNC_INFO << ", fit model type = " << modelType << ", degree = " << degree);
0109 
0110     QVector<double>& paramStartValues = fitData.paramStartValues;
0111     // QVector<double>* xVector = static_cast<QVector<double>* >(tmpXDataColumn->data());
0112 
0113     // double xmean = gsl_stats_mean(xVector->constData(), 1, tmpXDataColumn->rowCount());
0114     double xmin = xColumn->minimum();
0115     double xmax = xColumn->maximum();
0116     // double ymin = tmpYDataColumn->minimum();
0117     double ymax = yColumn->maximum();
0118     double xrange = xmax - xmin;
0119     // double yrange = ymax-ymin;
0120     DEBUG(Q_FUNC_INFO << ", x min/max = " << xmin << ' ' << xmax);
0121     // DEBUG(Q_FUNC_INFO <<", y min/max = " << ymin << ' ' << ymax);
0122 
0123     // guess start values of parameter
0124     switch (modelCategory) {
0125     case nsl_fit_model_basic:
0126         switch (modelType) {
0127         case nsl_fit_model_polynomial: { // do a multiparameter linear regression
0128             // copy all valid data point for the fit to temporary vectors
0129             QVector<double> xdataVector;
0130             QVector<double> ydataVector;
0131             QVector<double> yerrorVector;
0132 
0133             Range<double> xRange{xmin, xmax};
0134             if (fitData.autoRange) { // auto x range of data to fit
0135                 fitData.fitRange = xRange;
0136             } else { // custom x range of data to fit
0137                 if (!fitData.fitRange.isZero()) // avoid problems with user specified zero range
0138                     xRange.setRange(fitData.fitRange.start(), fitData.fitRange.end());
0139             }
0140 
0141             const int rowCount = std::min(xColumn->rowCount(), yColumn->rowCount());
0142             for (int row = 0; row < rowCount; ++row) {
0143                 // omit invalid data
0144                 if (!xColumn->isValid(row) || xColumn->isMasked(row) || !yColumn->isValid(row) || yColumn->isMasked(row))
0145                     continue;
0146 
0147                 double x = NAN;
0148                 switch (xColumn->columnMode()) {
0149                 case AbstractColumn::ColumnMode::Double:
0150                 case AbstractColumn::ColumnMode::Integer:
0151                 case AbstractColumn::ColumnMode::BigInt:
0152                     x = xColumn->valueAt(row);
0153                     break;
0154                 case AbstractColumn::ColumnMode::Text: // not valid
0155                     break;
0156                 case AbstractColumn::ColumnMode::DateTime:
0157                 case AbstractColumn::ColumnMode::Day:
0158                 case AbstractColumn::ColumnMode::Month:
0159                     x = xColumn->dateTimeAt(row).toMSecsSinceEpoch();
0160                 }
0161 
0162                 double y = NAN;
0163                 switch (yColumn->columnMode()) {
0164                 case AbstractColumn::ColumnMode::Double:
0165                 case AbstractColumn::ColumnMode::Integer:
0166                 case AbstractColumn::ColumnMode::BigInt:
0167                     y = yColumn->valueAt(row);
0168                     break;
0169                 case AbstractColumn::ColumnMode::Text: // not valid
0170                     break;
0171                 case AbstractColumn::ColumnMode::DateTime:
0172                 case AbstractColumn::ColumnMode::Day:
0173                 case AbstractColumn::ColumnMode::Month:
0174                     y = yColumn->dateTimeAt(row).toMSecsSinceEpoch();
0175                 }
0176 
0177                 if (x >= xRange.start() && x <= xRange.end()) { // only when inside given range
0178                     if (!yErrorColumn || !fitData.useDataErrors) { // x-y
0179                         xdataVector.append(x);
0180                         ydataVector.append(y);
0181                     } else if (yErrorColumn) { // x-y-dy
0182                         if (!std::isnan(yErrorColumn->valueAt(row))) {
0183                             xdataVector.append(x);
0184                             ydataVector.append(y);
0185                             yerrorVector.append(yErrorColumn->valueAt(row));
0186                         }
0187                     }
0188                 }
0189             }
0190 
0191             const int n = xdataVector.size();
0192             double* xdata = xdataVector.data();
0193             double* ydata = ydataVector.data();
0194             double* yerror = yerrorVector.data(); // size may be 0
0195 
0196             const double np = degree + 1;
0197             gsl_matrix* X = gsl_matrix_alloc(n, np); // X matrix
0198             gsl_vector* y = gsl_vector_alloc(n); // y values
0199             gsl_vector* w = gsl_vector_alloc(n); // weights
0200 
0201             gsl_vector* c = gsl_vector_alloc(np); // best fit parameter
0202             gsl_matrix* cov = gsl_matrix_alloc(np, np);
0203 
0204             const double minError = 1.e-199; // minimum error for weighting
0205             for (int i = 0; i < n; i++) {
0206                 double xi = xdata[i];
0207                 double yi = ydata[i];
0208 
0209                 for (int j = 0; j < np; j++)
0210                     gsl_matrix_set(X, i, j, gsl_pow_int(xi, j));
0211                 gsl_vector_set(y, i, yi);
0212 
0213                 if (yErrorColumn && i < yerrorVector.size()) {
0214                     switch (fitData.yWeightsType) {
0215                     case nsl_fit_weight_no:
0216                     case nsl_fit_weight_statistical_fit:
0217                     case nsl_fit_weight_relative_fit:
0218                         gsl_vector_set(w, i, 1.);
0219                         break;
0220                     case nsl_fit_weight_instrumental: // yerror are sigmas
0221                         gsl_vector_set(w, i, 1. / gsl_pow_2(std::max(yerror[i], std::max(sqrt(minError), std::abs(yi) * 1.e-15))));
0222                         break;
0223                     case nsl_fit_weight_direct: // yerror are weights
0224                         gsl_vector_set(w, i, yerror[i]);
0225                         break;
0226                     case nsl_fit_weight_inverse: // yerror are inverse weights
0227                         gsl_vector_set(w, i, 1. / std::max(yerror[i], std::max(minError, std::abs(yi) * 1.e-15)));
0228                         break;
0229                     case nsl_fit_weight_statistical:
0230                         gsl_vector_set(w, i, 1. / std::max(yi, minError));
0231                         break;
0232                     case nsl_fit_weight_relative:
0233                         gsl_vector_set(w, i, 1. / std::max(gsl_pow_2(yi), minError));
0234                         break;
0235                     }
0236                 } else
0237                     gsl_vector_set(w, i, 1.);
0238             }
0239 
0240             auto* work = gsl_multifit_linear_alloc(n, np);
0241             double chisq;
0242             int status = gsl_multifit_wlinear(X, w, y, c, cov, &chisq, work);
0243             gsl_multifit_linear_free(work);
0244             gsl_vector_free(w);
0245 
0246             for (int i = 0; i < np; i++) {
0247                 if (!std::isnan(gsl_vector_get(c, i)))
0248                     paramStartValues[i] = gsl_vector_get(c, i);
0249             }
0250 
0251             // results
0252             d->fitResult = XYFitCurve::FitResult(); // clear result
0253 
0254             d->fitResult.available = true;
0255             d->fitResult.valid = true;
0256             d->fitResult.status = gslErrorToString(status);
0257 
0258             d->fitResult.sse = chisq;
0259             d->fitResult.dof = n - np;
0260             // SST needed for coefficient of determination, R-squared and F test
0261             d->fitResult.sst = gsl_stats_tss(y->data, 1, n);
0262             // for a linear model without intercept R-squared is calculated differently
0263             // also using alternative R^2 when R^2 would be negative
0264             if (degree == 1 || d->fitResult.sst < d->fitResult.sse)
0265                 d->fitResult.sst = gsl_stats_tss_m(y->data, 1, n, 0);
0266 
0267             d->fitResult.calculateResult(n, np);
0268 
0269             d->fitResult.paramValues.resize(np);
0270             d->fitResult.errorValues.resize(np);
0271             d->fitResult.tdist_tValues.resize(np);
0272             d->fitResult.tdist_pValues.resize(np);
0273             d->fitResult.marginValues.resize(np);
0274 
0275             const double cerr = sqrt(d->fitResult.rms);
0276             // CI = 100 * (1 - alpha)
0277             const double alpha = 1.0 - fitData.confidenceInterval / 100.;
0278             for (unsigned int i = 0; i < np; i++) {
0279                 for (unsigned int j = 0; j <= i; j++)
0280                     d->fitResult.correlationMatrix << gsl_matrix_get(cov, i, j) / sqrt(gsl_matrix_get(cov, i, i)) / sqrt(gsl_matrix_get(cov, j, j));
0281                 d->fitResult.paramValues[i] = gsl_vector_get(c, i);
0282                 d->fitResult.errorValues[i] = cerr * sqrt(gsl_matrix_get(cov, i, i));
0283                 d->fitResult.tdist_tValues[i] = nsl_stats_tdist_t(d->fitResult.paramValues.at(i), d->fitResult.errorValues.at(i));
0284                 d->fitResult.tdist_pValues[i] = nsl_stats_tdist_p(d->fitResult.tdist_tValues.at(i), d->fitResult.dof);
0285                 d->fitResult.marginValues[i] = nsl_stats_tdist_margin(alpha, d->fitResult.dof, d->fitResult.errorValues.at(i));
0286             }
0287 
0288             // residuals
0289             gsl_vector* r = gsl_vector_alloc(n);
0290             status = gsl_multifit_linear_residuals(X, y, c, r);
0291             if (!status)
0292                 d->fitResult.mae = gsl_blas_dasum(r) / n;
0293             // TODO: show residuals?
0294 
0295             gsl_matrix_free(X);
0296             gsl_vector_free(y);
0297             gsl_vector_free(c);
0298             gsl_matrix_free(cov);
0299             break;
0300         }
0301         // TODO: use regression for all basic models?
0302         case nsl_fit_model_power: // a x^b, a + b x^c
0303         case nsl_fit_model_exponential: // a e^(bx), a1 e^(b1 x) + a2 e^(b2 x), ...
0304         case nsl_fit_model_inverse_exponential: // a (1-e^(bx)) + c
0305         case nsl_fit_model_fourier: // a0 + a1*sin(x) + b1 * cos(x) + ...
0306             break;
0307         }
0308         break;
0309     case nsl_fit_model_peak:
0310         // use equidistant mu's and (xmax-xmin)/(10*degree) as sigma(, gamma)
0311         switch (modelType) {
0312         case nsl_fit_model_gaussian:
0313         case nsl_fit_model_lorentz:
0314         case nsl_fit_model_sech:
0315         case nsl_fit_model_logistic:
0316             for (int d = 0; d < degree; d++) {
0317                 paramStartValues[3 * d + 2] = xmin + (d + 1.) * xrange / (degree + 1.); // mu
0318                 paramStartValues[3 * d + 1] = xrange / (10. * degree); // sigma
0319                 paramStartValues[3 * d] = paramStartValues[3 * d + 1] * ymax; // A = sigma * ymax
0320             }
0321             break;
0322         case nsl_fit_model_voigt:
0323             for (int d = 0; d < degree; d++) {
0324                 paramStartValues[4 * d + 1] = xmin + (d + 1.) * xrange / (degree + 1.); // mu
0325                 paramStartValues[4 * d + 2] = xrange / (10. * degree); // sigma
0326                 paramStartValues[4 * d + 3] = xrange / (10. * degree); // gamma
0327             }
0328             break;
0329         case nsl_fit_model_pseudovoigt1:
0330             for (int d = 0; d < degree; d++) {
0331                 paramStartValues[4 * d + 1] = 0.5; // eta
0332                 paramStartValues[4 * d + 2] = xrange / (10. * degree); // sigma
0333                 paramStartValues[4 * d + 3] = xmin + (d + 1.) * xrange / (degree + 1.); // mu
0334             }
0335             break;
0336         }
0337         break;
0338     case nsl_fit_model_growth:
0339         switch (modelType) {
0340         case nsl_fit_model_atan:
0341         case nsl_fit_model_tanh:
0342         case nsl_fit_model_algebraic_sigmoid:
0343         case nsl_fit_model_erf:
0344         case nsl_fit_model_gudermann:
0345         case nsl_fit_model_sigmoid:
0346             // use (xmax+xmin)/2 as mu and (xmax-xmin)/10 as sigma
0347             paramStartValues[1] = (xmax + xmin) / 2.;
0348             paramStartValues[2] = xrange / 10.;
0349             break;
0350         case nsl_fit_model_hill:
0351             paramStartValues[2] = xrange / 10.;
0352             break;
0353         case nsl_fit_model_gompertz:
0354             // TODO
0355             break;
0356         }
0357         break;
0358     case nsl_fit_model_distribution:
0359         switch (modelType) {
0360         case nsl_sf_stats_gaussian:
0361         case nsl_sf_stats_laplace:
0362         case nsl_sf_stats_rayleigh_tail:
0363         case nsl_sf_stats_lognormal:
0364         case nsl_sf_stats_logistic:
0365         case nsl_sf_stats_sech:
0366         case nsl_sf_stats_cauchy_lorentz:
0367         case nsl_sf_stats_levy:
0368             // mu
0369             paramStartValues[2] = (xmin + xmax) / 2.;
0370             // sigma
0371             paramStartValues[1] = xrange / 10.;
0372             // A = sigma * y_max
0373             paramStartValues[0] = paramStartValues[1] * ymax;
0374             break;
0375         // TODO: other types
0376         default:
0377             break;
0378         }
0379         break;
0380     case nsl_fit_model_custom:
0381         // not possible
0382         break;
0383     }
0384 }
0385 
0386 /*!
0387  * sets the parameter names for given model category, model type and degree in \c fitData for given action
0388  */
0389 void XYFitCurve::initFitData(XYAnalysisCurve::AnalysisAction action) {
0390     // TODO: exclude others too?
0391     if (action == XYAnalysisCurve::AnalysisAction::DataReduction)
0392         return;
0393 
0394     Q_D(XYFitCurve);
0395     XYFitCurve::FitData& fitData = d->fitData;
0396     if (action == XYAnalysisCurve::AnalysisAction::FitLinear) {
0397         // Linear
0398         fitData.modelCategory = nsl_fit_model_basic;
0399         fitData.modelType = (int)nsl_fit_model_polynomial;
0400         fitData.degree = 1;
0401     } else if (action == XYAnalysisCurve::AnalysisAction::FitPower) {
0402         // Power
0403         fitData.modelCategory = nsl_fit_model_basic;
0404         fitData.modelType = (int)nsl_fit_model_power;
0405         fitData.degree = 1;
0406     } else if (action == XYAnalysisCurve::AnalysisAction::FitExp1) {
0407         // Exponential (degree 1)
0408         fitData.modelCategory = nsl_fit_model_basic;
0409         fitData.modelType = (int)nsl_fit_model_exponential;
0410         fitData.degree = 1;
0411     } else if (action == XYAnalysisCurve::AnalysisAction::FitExp2) {
0412         // Exponential (degree 2)
0413         fitData.modelCategory = nsl_fit_model_basic;
0414         fitData.modelType = (int)nsl_fit_model_exponential;
0415         fitData.degree = 2;
0416     } else if (action == XYAnalysisCurve::AnalysisAction::FitInvExp) {
0417         // Inverse exponential
0418         fitData.modelCategory = nsl_fit_model_basic;
0419         fitData.modelType = (int)nsl_fit_model_inverse_exponential;
0420     } else if (action == XYAnalysisCurve::AnalysisAction::FitGauss) {
0421         // Gauss
0422         fitData.modelCategory = nsl_fit_model_peak;
0423         fitData.modelType = (int)nsl_fit_model_gaussian;
0424         fitData.degree = 1;
0425     } else if (action == XYAnalysisCurve::AnalysisAction::FitCauchyLorentz) {
0426         // Cauchy-Lorentz
0427         fitData.modelCategory = nsl_fit_model_peak;
0428         fitData.modelType = (int)nsl_fit_model_lorentz;
0429         fitData.degree = 1;
0430     } else if (action == XYAnalysisCurve::AnalysisAction::FitTan) {
0431         // Arc tangent
0432         fitData.modelCategory = nsl_fit_model_growth;
0433         fitData.modelType = (int)nsl_fit_model_atan;
0434     } else if (action == XYAnalysisCurve::AnalysisAction::FitTanh) {
0435         // Hyperbolic tangent
0436         fitData.modelCategory = nsl_fit_model_growth;
0437         fitData.modelType = (int)nsl_fit_model_tanh;
0438     } else if (action == XYAnalysisCurve::AnalysisAction::FitErrFunc) {
0439         // Error function
0440         fitData.modelCategory = nsl_fit_model_growth;
0441         fitData.modelType = (int)nsl_fit_model_erf;
0442     } else {
0443         // Custom
0444         fitData.modelCategory = nsl_fit_model_custom;
0445         fitData.modelType = 0;
0446     }
0447 
0448     XYFitCurve::initFitData(fitData);
0449 }
0450 
0451 /*!
0452  * sets the model expression and the parameter names for given model category, model type and degree in \c fitData
0453  */
0454 void XYFitCurve::initFitData(XYFitCurve::FitData& fitData) {
0455     nsl_fit_model_category modelCategory = fitData.modelCategory;
0456     int modelType = fitData.modelType;
0457     QString& model = fitData.model;
0458     QStringList& paramNames = fitData.paramNames;
0459     QStringList& paramNamesUtf8 = fitData.paramNamesUtf8;
0460     int degree = fitData.degree;
0461     QVector<double>& paramStartValues = fitData.paramStartValues;
0462     QVector<double>& paramLowerLimits = fitData.paramLowerLimits;
0463     QVector<double>& paramUpperLimits = fitData.paramUpperLimits;
0464     QVector<bool>& paramFixed = fitData.paramFixed;
0465 
0466     if (modelCategory != nsl_fit_model_custom) {
0467         DEBUG(Q_FUNC_INFO << ", XYFitCurve::initFitData() for model category = " << nsl_fit_model_category_name[modelCategory] << ", model type = " << modelType
0468                           << ", degree = " << degree);
0469         paramNames.clear();
0470     } else {
0471         DEBUG(Q_FUNC_INFO << ", XYFitCurve::initFitData() for model category = nsl_fit_model_custom, model type = " << modelType << ", degree = " << degree);
0472     }
0473     paramNamesUtf8.clear();
0474 
0475     // 10 indices used in multi degree models
0476     QStringList indices = {UTF8_QSTRING("₁"),
0477                            UTF8_QSTRING("₂"),
0478                            UTF8_QSTRING("₃"),
0479                            UTF8_QSTRING("₄"),
0480                            UTF8_QSTRING("₅"),
0481                            UTF8_QSTRING("₆"),
0482                            UTF8_QSTRING("₇"),
0483                            UTF8_QSTRING("₈"),
0484                            UTF8_QSTRING("₉"),
0485                            UTF8_QSTRING("₁₀")};
0486 
0487     switch (modelCategory) {
0488     case nsl_fit_model_basic:
0489         model = QLatin1String(nsl_fit_model_basic_equation[fitData.modelType]);
0490         switch (modelType) {
0491         case nsl_fit_model_polynomial:
0492             paramNames << QStringLiteral("c0") << QStringLiteral("c1");
0493             paramNamesUtf8 << UTF8_QSTRING("c₀") << UTF8_QSTRING("c₁");
0494             if (degree == 2) {
0495                 model += QStringLiteral(" + c2*x^2");
0496                 paramNames << QStringLiteral("c2");
0497                 paramNamesUtf8 << UTF8_QSTRING("c₂");
0498             } else if (degree > 2) {
0499                 for (int i = 2; i <= degree; ++i) {
0500                     QString numStr = QString::number(i);
0501                     model += QStringLiteral("+c") + numStr + QStringLiteral("*x^") + numStr;
0502                     paramNames << QStringLiteral("c") + numStr;
0503                     paramNamesUtf8 << QStringLiteral("c") + indices[i - 1];
0504                 }
0505             }
0506             break;
0507         case nsl_fit_model_power:
0508             if (degree == 1) {
0509                 paramNames << QStringLiteral("a") << QStringLiteral("b");
0510             } else {
0511                 paramNames << QStringLiteral("a") << QStringLiteral("b") << QStringLiteral("c");
0512                 model = QStringLiteral("a + b*x^c");
0513             }
0514             break;
0515         case nsl_fit_model_exponential:
0516             if (degree == 1) {
0517                 paramNames << QStringLiteral("a") << QStringLiteral("b");
0518             } else {
0519                 for (int i = 1; i <= degree; i++) {
0520                     QString numStr = QString::number(i);
0521                     if (i == 1)
0522                         model = QStringLiteral("a1*exp(b1*x)");
0523                     else
0524                         model += QStringLiteral(" + a") + numStr + QStringLiteral("*exp(b") + numStr + QStringLiteral("*x)");
0525                     paramNames << QStringLiteral("a") + numStr << QStringLiteral("b") + numStr;
0526                     paramNamesUtf8 << QStringLiteral("a") + indices[i - 1] << QStringLiteral("b") + indices[i - 1];
0527                 }
0528             }
0529             break;
0530         case nsl_fit_model_inverse_exponential:
0531             paramNames << QStringLiteral("a") << QStringLiteral("b") << QStringLiteral("c");
0532             break;
0533         case nsl_fit_model_fourier:
0534             paramNames << QStringLiteral("w") << QStringLiteral("a0") << QStringLiteral("a1") << QStringLiteral("b1");
0535             paramNamesUtf8 << UTF8_QSTRING("ω") << UTF8_QSTRING("a₀") << UTF8_QSTRING("a₁") << UTF8_QSTRING("b₁");
0536             if (degree > 1) {
0537                 for (int i = 1; i <= degree; ++i) {
0538                     QString numStr = QString::number(i);
0539                     model += QStringLiteral("+ (a") + numStr + QStringLiteral("*cos(") + numStr + QStringLiteral("*w*x) + b") + numStr + QStringLiteral("*sin(")
0540                         + numStr + QStringLiteral("*w*x))");
0541                     paramNames << QStringLiteral("a") + numStr << QStringLiteral("b") + numStr;
0542                     paramNamesUtf8 << QStringLiteral("a") + indices[i - 1] << QStringLiteral("b") + indices[i - 1];
0543                 }
0544             }
0545             break;
0546         }
0547         break;
0548     case nsl_fit_model_peak:
0549         model = QLatin1String(nsl_fit_model_peak_equation[fitData.modelType]);
0550         switch (modelType) {
0551         case nsl_fit_model_gaussian:
0552             switch (degree) {
0553             case 1:
0554                 paramNames << QStringLiteral("a") << QStringLiteral("s") << QStringLiteral("mu");
0555                 paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("σ") << UTF8_QSTRING("μ");
0556                 break;
0557             default:
0558                 model = QStringLiteral("1/sqrt(2*pi) * (");
0559                 for (int i = 1; i <= degree; ++i) {
0560                     QString numStr = QString::number(i);
0561                     if (i > 1)
0562                         model += QStringLiteral(" + ");
0563                     model += QStringLiteral("a") + numStr + QStringLiteral("/s") + numStr + QStringLiteral("* exp(-((x-mu") + numStr + QStringLiteral(")/s")
0564                         + numStr + QStringLiteral(")^2/2)");
0565                     paramNames << QStringLiteral("a") + numStr << QStringLiteral("s") + numStr << QStringLiteral("mu") + numStr;
0566                     paramNamesUtf8 << QStringLiteral("A") + indices[i - 1] << UTF8_QSTRING("σ") + indices[i - 1] << UTF8_QSTRING("μ") + indices[i - 1];
0567                 }
0568                 model += QStringLiteral(")");
0569             }
0570             break;
0571         case nsl_fit_model_lorentz:
0572             switch (degree) {
0573             case 1:
0574                 paramNames << QStringLiteral("a") << QStringLiteral("g") << QStringLiteral("mu");
0575                 paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("γ") << UTF8_QSTRING("μ");
0576                 break;
0577             default:
0578                 model = QStringLiteral("1/pi * (");
0579                 for (int i = 1; i <= degree; ++i) {
0580                     QString numStr = QString::number(i);
0581                     if (i > 1)
0582                         model += QStringLiteral(" + ");
0583                     model += QStringLiteral("a") + numStr + QStringLiteral(" * g") + numStr + QStringLiteral("/(g") + numStr + QStringLiteral("^2+(x-mu")
0584                         + numStr + QStringLiteral(")^2)");
0585                     paramNames << QStringLiteral("a") + numStr << QStringLiteral("g") + numStr << QStringLiteral("mu") + numStr;
0586                     paramNamesUtf8 << QStringLiteral("A") + indices[i - 1] << UTF8_QSTRING("γ") + indices[i - 1] << UTF8_QSTRING("μ") + indices[i - 1];
0587                 }
0588                 model += QStringLiteral(")");
0589             }
0590             break;
0591         case nsl_fit_model_sech:
0592             switch (degree) {
0593             case 1:
0594                 paramNames << QStringLiteral("a") << QStringLiteral("s") << QStringLiteral("mu");
0595                 paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("σ") << UTF8_QSTRING("μ");
0596                 break;
0597             default:
0598                 model = QStringLiteral("1/pi * (");
0599                 for (int i = 1; i <= degree; ++i) {
0600                     QString numStr = QString::number(i);
0601                     if (i > 1)
0602                         model += QStringLiteral(" + ");
0603                     model += QStringLiteral("a") + numStr + QStringLiteral("/s") + numStr + QStringLiteral("* sech((x-mu") + numStr + QStringLiteral(")/s")
0604                         + numStr + QStringLiteral(")");
0605                     paramNames << QStringLiteral("a") + numStr << QStringLiteral("s") + numStr << QStringLiteral("mu") + numStr;
0606                     paramNamesUtf8 << QStringLiteral("A") + indices[i - 1] << UTF8_QSTRING("σ") + indices[i - 1] << UTF8_QSTRING("μ") + indices[i - 1];
0607                 }
0608                 model += QStringLiteral(")");
0609             }
0610             break;
0611         case nsl_fit_model_logistic:
0612             switch (degree) {
0613             case 1:
0614                 paramNames << QStringLiteral("a") << QStringLiteral("s") << QStringLiteral("mu");
0615                 paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("σ") << UTF8_QSTRING("μ");
0616                 break;
0617             default:
0618                 model = QStringLiteral("1/4 * (");
0619                 for (int i = 1; i <= degree; ++i) {
0620                     QString numStr = QString::number(i);
0621                     if (i > 1)
0622                         model += QStringLiteral(" + ");
0623                     model += QStringLiteral("a") + numStr + QStringLiteral("/s") + numStr + QStringLiteral("* sech((x-mu") + numStr + QStringLiteral(")/2/s")
0624                         + numStr + QStringLiteral(")**2");
0625                     paramNames << QStringLiteral("a") + numStr << QStringLiteral("s") + numStr << QStringLiteral("mu") + numStr;
0626                     paramNamesUtf8 << QStringLiteral("A") + indices[i - 1] << UTF8_QSTRING("σ") + indices[i - 1] << UTF8_QSTRING("μ") + indices[i - 1];
0627                 }
0628                 model += QStringLiteral(")");
0629             }
0630             break;
0631         case nsl_fit_model_voigt:
0632             switch (degree) {
0633             case 1:
0634                 paramNames << QStringLiteral("a") << QStringLiteral("mu") << QStringLiteral("s") << QStringLiteral("g");
0635                 paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("μ") << UTF8_QSTRING("σ") << UTF8_QSTRING("γ");
0636                 break;
0637             default:
0638                 model.clear();
0639                 for (int i = 1; i <= degree; ++i) {
0640                     QString numStr = QString::number(i);
0641                     if (i > 1)
0642                         model += QStringLiteral(" + ");
0643                     model += QStringLiteral("a") + numStr + QStringLiteral("*voigt(x-mu") + numStr + QStringLiteral(",s") + numStr + QStringLiteral(",g")
0644                         + numStr + QStringLiteral(")");
0645                     paramNames << QStringLiteral("a") + numStr << QStringLiteral("mu") + numStr << QStringLiteral("s") + numStr << QStringLiteral("g") + numStr;
0646                     paramNamesUtf8 << QStringLiteral("A") + indices[i - 1] << UTF8_QSTRING("μ") + indices[i - 1] << UTF8_QSTRING("σ") + indices[i - 1]
0647                                    << UTF8_QSTRING("γ") + indices[i - 1];
0648                 }
0649             }
0650             break;
0651         case nsl_fit_model_pseudovoigt1:
0652             switch (degree) {
0653             case 1:
0654                 paramNames << QStringLiteral("a") << QStringLiteral("et") << QStringLiteral("w") << QStringLiteral("mu"); // eta already exists as function!
0655                 paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("η") << QStringLiteral("w") << UTF8_QSTRING("μ");
0656                 break;
0657             default:
0658                 model.clear();
0659                 for (int i = 1; i <= degree; ++i) {
0660                     QString numStr = QString::number(i);
0661                     if (i > 1)
0662                         model += QStringLiteral(" + ");
0663                     model += QStringLiteral("a") + numStr + QStringLiteral("*pseudovoigt1(x-mu") + numStr + QStringLiteral(",eta") + numStr
0664                         + QStringLiteral(",w") + numStr + QStringLiteral(")");
0665                     paramNames << QStringLiteral("a") + numStr << QStringLiteral("eta") + numStr << QStringLiteral("w") + numStr
0666                                << QStringLiteral("mu") + numStr;
0667                     paramNamesUtf8 << QStringLiteral("A") + indices[i - 1] << UTF8_QSTRING("η") + indices[i - 1] << QStringLiteral("w") + indices[i - 1]
0668                                    << UTF8_QSTRING("μ") + indices[i - 1];
0669                 }
0670             }
0671             break;
0672         }
0673         break;
0674     case nsl_fit_model_growth:
0675         model = QLatin1String(nsl_fit_model_growth_equation[fitData.modelType]);
0676         switch (modelType) {
0677         case nsl_fit_model_atan:
0678         case nsl_fit_model_tanh:
0679         case nsl_fit_model_algebraic_sigmoid:
0680         case nsl_fit_model_erf:
0681         case nsl_fit_model_gudermann:
0682             paramNames << QStringLiteral("a") << QStringLiteral("mu") << QStringLiteral("s");
0683             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("μ") << UTF8_QSTRING("σ");
0684             break;
0685         case nsl_fit_model_sigmoid:
0686             paramNames << QStringLiteral("a") << QStringLiteral("mu") << QStringLiteral("k");
0687             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("μ") << QStringLiteral("k");
0688             break;
0689         case nsl_fit_model_hill:
0690             paramNames << QStringLiteral("a") << QStringLiteral("n") << QStringLiteral("a");
0691             paramNamesUtf8 << QStringLiteral("A") << QStringLiteral("n") << UTF8_QSTRING("σ");
0692             break;
0693         case nsl_fit_model_gompertz:
0694             paramNames << QStringLiteral("a") << QStringLiteral("b") << QStringLiteral("c");
0695             break;
0696         }
0697         break;
0698     case nsl_fit_model_distribution:
0699         model = QLatin1String(nsl_sf_stats_distribution_equation[fitData.modelType]);
0700         switch (modelType) {
0701         case nsl_sf_stats_gaussian:
0702         case nsl_sf_stats_laplace:
0703         case nsl_sf_stats_rayleigh_tail:
0704         case nsl_sf_stats_lognormal:
0705         case nsl_sf_stats_logistic:
0706         case nsl_sf_stats_sech:
0707             paramNames << QStringLiteral("a") << QStringLiteral("s") << QStringLiteral("mu");
0708             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("σ") << UTF8_QSTRING("μ");
0709             break;
0710         case nsl_sf_stats_gaussian_tail:
0711             paramNames << QStringLiteral("A") << QStringLiteral("s") << QStringLiteral("a") << QStringLiteral("mu");
0712             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("σ") << QStringLiteral("a") << UTF8_QSTRING("μ");
0713             break;
0714         case nsl_sf_stats_exponential:
0715             paramNames << QStringLiteral("a") << QStringLiteral("l") << QStringLiteral("mu");
0716             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("λ") << UTF8_QSTRING("μ");
0717             break;
0718         case nsl_sf_stats_exponential_power:
0719             paramNames << QStringLiteral("a") << QStringLiteral("s") << QStringLiteral("b") << QStringLiteral("mu");
0720             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("σ") << QStringLiteral("b") << UTF8_QSTRING("μ");
0721             break;
0722         case nsl_sf_stats_cauchy_lorentz:
0723         case nsl_sf_stats_levy:
0724             paramNames << QStringLiteral("a") << QStringLiteral("g") << QStringLiteral("mu");
0725             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("γ") << UTF8_QSTRING("μ");
0726             break;
0727         case nsl_sf_stats_rayleigh:
0728             paramNames << QStringLiteral("a") << QStringLiteral("s");
0729             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("σ");
0730             break;
0731         case nsl_sf_stats_landau:
0732             paramNames << QStringLiteral("a");
0733             paramNamesUtf8 << QStringLiteral("A");
0734             break;
0735         case nsl_sf_stats_levy_alpha_stable: // unused distributions
0736         case nsl_sf_stats_levy_skew_alpha_stable:
0737         case nsl_sf_stats_bernoulli:
0738             break;
0739         case nsl_sf_stats_gamma:
0740             paramNames << QStringLiteral("a") << QStringLiteral("k") << QStringLiteral("t");
0741             paramNamesUtf8 << QStringLiteral("A") << QStringLiteral("k") << UTF8_QSTRING("θ");
0742             break;
0743         case nsl_sf_stats_flat:
0744             paramNames << QStringLiteral("A") << QStringLiteral("b") << QStringLiteral("a");
0745             break;
0746         case nsl_sf_stats_chi_squared:
0747             paramNames << QStringLiteral("a") << QStringLiteral("n");
0748             paramNamesUtf8 << QStringLiteral("A") << QStringLiteral("n");
0749             break;
0750         case nsl_sf_stats_fdist:
0751             paramNames << QStringLiteral("a") << QStringLiteral("n1") << QStringLiteral("n2");
0752             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("ν₁") << UTF8_QSTRING("ν₂");
0753             break;
0754         case nsl_sf_stats_tdist:
0755             paramNames << QStringLiteral("a") << QStringLiteral("n");
0756             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("ν");
0757             break;
0758         case nsl_sf_stats_beta:
0759         case nsl_sf_stats_pareto:
0760             paramNames << QStringLiteral("A") << QStringLiteral("a") << QStringLiteral("b");
0761             break;
0762         case nsl_sf_stats_weibull:
0763             paramNames << QStringLiteral("a") << QStringLiteral("k") << QStringLiteral("l") << QStringLiteral("mu");
0764             paramNamesUtf8 << QStringLiteral("A") << QStringLiteral("k") << UTF8_QSTRING("λ") << UTF8_QSTRING("μ");
0765             break;
0766         case nsl_sf_stats_gumbel1:
0767             paramNames << QStringLiteral("a") << QStringLiteral("s") << QStringLiteral("mu") << QStringLiteral("b");
0768             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("σ") << UTF8_QSTRING("μ") << UTF8_QSTRING("β");
0769             break;
0770         case nsl_sf_stats_gumbel2:
0771             paramNames << QStringLiteral("A") << QStringLiteral("a") << QStringLiteral("b") << QStringLiteral("mu");
0772             paramNamesUtf8 << QStringLiteral("A") << QStringLiteral("a") << QStringLiteral("b") << UTF8_QSTRING("μ");
0773             break;
0774         case nsl_sf_stats_poisson:
0775             paramNames << QStringLiteral("a") << QStringLiteral("l");
0776             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("λ");
0777             break;
0778         case nsl_sf_stats_binomial:
0779         case nsl_sf_stats_negative_binomial:
0780         case nsl_sf_stats_pascal:
0781             paramNames << QStringLiteral("a") << QStringLiteral("p") << QStringLiteral("n");
0782             paramNamesUtf8 << QStringLiteral("A") << QStringLiteral("p") << QStringLiteral("n");
0783             break;
0784         case nsl_sf_stats_geometric:
0785         case nsl_sf_stats_logarithmic:
0786             paramNames << QStringLiteral("a") << QStringLiteral("p");
0787             paramNamesUtf8 << QStringLiteral("A") << QStringLiteral("p");
0788             break;
0789         case nsl_sf_stats_hypergeometric:
0790             paramNames << QStringLiteral("a") << QStringLiteral("n1") << QStringLiteral("n2") << QStringLiteral("t");
0791             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("n₁") << UTF8_QSTRING("n₂") << QStringLiteral("t");
0792             break;
0793         case nsl_sf_stats_maxwell_boltzmann:
0794             paramNames << QStringLiteral("a") << QStringLiteral("s");
0795             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("σ");
0796             break;
0797         case nsl_sf_stats_frechet:
0798             paramNames << QStringLiteral("a") << QStringLiteral("g") << QStringLiteral("s") << QStringLiteral("mu");
0799             paramNamesUtf8 << QStringLiteral("A") << UTF8_QSTRING("γ") << UTF8_QSTRING("σ") << UTF8_QSTRING("μ");
0800             break;
0801         }
0802         break;
0803     case nsl_fit_model_custom:
0804         break;
0805     }
0806     DEBUG(Q_FUNC_INFO << ", model: " << STDSTRING(model));
0807     DEBUG(Q_FUNC_INFO << ", # params: " << paramNames.size());
0808 
0809     if (paramNamesUtf8.isEmpty())
0810         paramNamesUtf8 << paramNames;
0811 
0812     // resize the vector for the start values and set the elements to 1.0
0813     // in case a custom model is used, do nothing, we take over the previous values
0814     if (modelCategory != nsl_fit_model_custom) {
0815         const int np = paramNames.size();
0816         paramStartValues.resize(np);
0817         paramFixed.resize(np);
0818         paramLowerLimits.resize(np);
0819         paramUpperLimits.resize(np);
0820 
0821         for (int i = 0; i < np; ++i) {
0822             paramStartValues[i] = 1.0;
0823             paramFixed[i] = false;
0824             paramLowerLimits[i] = -std::numeric_limits<double>::max();
0825             paramUpperLimits[i] = std::numeric_limits<double>::max();
0826         }
0827 
0828         // set some model-dependent start values
0829         // TODO: see initStartValues()
0830         if (modelCategory == nsl_fit_model_distribution) {
0831             if (modelType == (int)nsl_sf_stats_flat)
0832                 paramStartValues[2] = -1.0;
0833             else if (modelType == (int)nsl_sf_stats_levy)
0834                 paramStartValues[2] = 0.0;
0835             else if (modelType == (int)nsl_sf_stats_exponential_power || modelType == (int)nsl_sf_stats_weibull || modelType == (int)nsl_sf_stats_gumbel2
0836                      || modelType == (int)nsl_sf_stats_frechet)
0837                 paramStartValues[3] = 0.0;
0838             else if (modelType == (int)nsl_sf_stats_binomial || modelType == (int)nsl_sf_stats_negative_binomial || modelType == (int)nsl_sf_stats_pascal
0839                      || modelType == (int)nsl_sf_stats_geometric || modelType == (int)nsl_sf_stats_logarithmic)
0840                 paramStartValues[1] = 0.5;
0841         }
0842     }
0843 }
0844 
0845 void XYFitCurve::clearFitResult() {
0846     Q_D(XYFitCurve);
0847     d->fitResult = XYFitCurve::FitResult();
0848 }
0849 
0850 void XYFitCurve::FitResult::calculateResult(size_t n, unsigned int np) {
0851     if (dof != 0) {
0852         rms = sse / dof;
0853         rsd = std::sqrt(rms);
0854     }
0855 
0856     mse = sse / n;
0857     rmse = std::sqrt(mse);
0858 
0859     rsquare = nsl_stats_rsquare(sse, sst);
0860     rsquareAdj = nsl_stats_rsquareAdj(rsquare, np, dof, 1);
0861     chisq_p = nsl_stats_chisq_p(sse, dof);
0862     fdist_F = nsl_stats_fdist_F(rsquare, np, dof);
0863     fdist_p = nsl_stats_fdist_p(fdist_F, np, dof);
0864     logLik = nsl_stats_logLik(sse, n);
0865     aic = nsl_stats_aic(sse, n, np, 1);
0866     bic = nsl_stats_bic(sse, n, np, 1);
0867 }
0868 
0869 /*!
0870     Returns an icon to be used in the project explorer.
0871 */
0872 QIcon XYFitCurve::icon() const {
0873     return QIcon::fromTheme(QStringLiteral("labplot-xy-fit-curve"));
0874 }
0875 
0876 // ##############################################################################
0877 // ##########################  getter methods  ##################################
0878 // ##############################################################################
0879 const AbstractColumn* XYFitCurve::residualsColumn() const {
0880     Q_D(const XYFitCurve);
0881     return d->residualsColumn;
0882 }
0883 
0884 BASIC_SHARED_D_READER_IMPL(XYFitCurve, const AbstractColumn*, xErrorColumn, xErrorColumn)
0885 BASIC_SHARED_D_READER_IMPL(XYFitCurve, const AbstractColumn*, yErrorColumn, yErrorColumn)
0886 const QString& XYFitCurve::xErrorColumnPath() const {
0887     Q_D(const XYFitCurve);
0888     return d->xErrorColumnPath;
0889 }
0890 const QString& XYFitCurve::yErrorColumnPath() const {
0891     Q_D(const XYFitCurve);
0892     return d->yErrorColumnPath;
0893 }
0894 
0895 BASIC_SHARED_D_READER_IMPL(XYFitCurve, const Histogram*, dataSourceHistogram, dataSourceHistogram)
0896 const QString& XYFitCurve::dataSourceHistogramPath() const {
0897     Q_D(const XYFitCurve);
0898     return d->dataSourceHistogramPath;
0899 }
0900 
0901 BASIC_SHARED_D_READER_IMPL(XYFitCurve, XYFitCurve::FitData, fitData, fitData)
0902 
0903 const XYFitCurve::FitResult& XYFitCurve::fitResult() const {
0904     Q_D(const XYFitCurve);
0905     return d->fitResult;
0906 }
0907 
0908 // ##############################################################################
0909 // #################  setter methods and undo commands ##########################
0910 // ##############################################################################
0911 STD_SETTER_CMD_IMPL_F_S(XYFitCurve, SetDataSourceHistogram, const Histogram*, dataSourceHistogram, retransform)
0912 void XYFitCurve::setDataSourceHistogram(const Histogram* histogram) {
0913     Q_D(XYFitCurve);
0914     if (histogram != d->dataSourceHistogram) {
0915         exec(new XYFitCurveSetDataSourceHistogramCmd(d, histogram, ki18n("%1: data source histogram changed")));
0916         handleSourceDataChanged();
0917 
0918         connect(histogram, &Histogram::dataChanged, this, &XYFitCurve::handleSourceDataChanged);
0919         // TODO: add disconnect in the undo-function
0920     }
0921 }
0922 
0923 STD_SETTER_CMD_IMPL_S(XYFitCurve, SetXErrorColumn, const AbstractColumn*, xErrorColumn)
0924 void XYFitCurve::setXErrorColumn(const AbstractColumn* column) {
0925     Q_D(XYFitCurve);
0926     if (column != d->xErrorColumn) {
0927         exec(new XYFitCurveSetXErrorColumnCmd(d, column, ki18n("%1: assign x-error")));
0928         handleSourceDataChanged();
0929         if (column) {
0930             connect(column, &AbstractColumn::dataChanged, this, [=]() {
0931                 handleSourceDataChanged();
0932             });
0933             // TODO: disconnect on undo
0934         }
0935     }
0936 }
0937 
0938 STD_SETTER_CMD_IMPL_S(XYFitCurve, SetYErrorColumn, const AbstractColumn*, yErrorColumn)
0939 void XYFitCurve::setYErrorColumn(const AbstractColumn* column) {
0940     Q_D(XYFitCurve);
0941     if (column != d->yErrorColumn) {
0942         exec(new XYFitCurveSetYErrorColumnCmd(d, column, ki18n("%1: assign y-error")));
0943         handleSourceDataChanged();
0944         if (column) {
0945             connect(column, &AbstractColumn::dataChanged, this, [=]() {
0946                 handleSourceDataChanged();
0947             });
0948             // TODO: disconnect on undo
0949         }
0950     }
0951 }
0952 
0953 // do not recalculate (allow preview)
0954 // STD_SETTER_CMD_IMPL_F_S(XYFitCurve, SetFitData, XYFitCurve::FitData, fitData, recalculate)
0955 STD_SETTER_CMD_IMPL_S(XYFitCurve, SetFitData, XYFitCurve::FitData, fitData)
0956 void XYFitCurve::setFitData(const XYFitCurve::FitData& fitData) {
0957     Q_D(XYFitCurve);
0958     exec(new XYFitCurveSetFitDataCmd(d, fitData, ki18n("%1: set fit options and perform the fit")));
0959 }
0960 
0961 // ##############################################################################
0962 // ######################### Private implementation #############################
0963 // ##############################################################################
0964 XYFitCurvePrivate::XYFitCurvePrivate(XYFitCurve* owner)
0965     : XYAnalysisCurvePrivate(owner)
0966     , q(owner) {
0967 }
0968 
0969 // no need to delete xColumn and yColumn, they are deleted
0970 // when the parent aspect is removed
0971 XYFitCurvePrivate::~XYFitCurvePrivate() = default;
0972 
0973 // data structure to pass parameter to fit functions
0974 struct data {
0975     size_t n; // number of data points
0976     double* x; // pointer to the vector with x-data values
0977     double* y; // pointer to the vector with y-data values
0978     double* weight; // pointer to the vector with weight values
0979     nsl_fit_model_category modelCategory;
0980     int modelType;
0981     int degree;
0982     QString* func; // string containing the definition of the model/function
0983     QStringList* paramNames;
0984     double* paramMin; // lower parameter limits
0985     double* paramMax; // upper parameter limits
0986     bool* paramFixed; // parameter fixed?
0987 };
0988 
0989 /*!
0990  * \param paramValues vector containing current values of the fit parameters
0991  * \param params
0992  * \param f vector with the weighted residuals weight[i]*(Yi - y[i])
0993  */
0994 int func_f(const gsl_vector* paramValues, void* params, gsl_vector* f) {
0995     // DEBUG(Q_FUNC_INFO);
0996     size_t n = ((struct data*)params)->n;
0997     double* x = ((struct data*)params)->x;
0998     double* y = ((struct data*)params)->y;
0999     double* weight = ((struct data*)params)->weight;
1000     nsl_fit_model_category modelCategory = ((struct data*)params)->modelCategory;
1001     unsigned int modelType = ((struct data*)params)->modelType;
1002     QStringList* paramNames = ((struct data*)params)->paramNames;
1003     double* min = ((struct data*)params)->paramMin;
1004     double* max = ((struct data*)params)->paramMax;
1005 
1006     // set current values of the parameters
1007     for (int i = 0; i < paramNames->size(); i++) {
1008         double v = gsl_vector_get(paramValues, (size_t)i);
1009         // bound values if limits are set
1010         assign_symbol(qPrintable(paramNames->at(i)), nsl_fit_map_bound(v, min[i], max[i]));
1011         QDEBUG("Parameter" << i << " (' " << paramNames->at(i) << "')" << '[' << min[i] << ',' << max[i] << "] free/bound:" << QString::number(v, 'g', 15)
1012                            << ' ' << QString::number(nsl_fit_map_bound(v, min[i], max[i]), 'g', 15));
1013     }
1014 
1015     QString func{*(((struct data*)params)->func)};
1016     for (size_t i = 0; i < n; i++) {
1017         if (std::isnan(x[i]) || std::isnan(y[i]))
1018             continue;
1019 
1020         // checks for allowed values of x for different models
1021         // TODO: more to check
1022         if (modelCategory == nsl_fit_model_distribution && modelType == nsl_sf_stats_lognormal) {
1023             if (x[i] < 0)
1024                 x[i] = 0;
1025         }
1026 
1027         assign_symbol("x", x[i]);
1028         // DEBUG("evaluate function \"" << STDSTRING(func) << "\" @ x = " << x[i] << ":");
1029         double Yi = parse(qPrintable(func), qPrintable(QLocale().name()));
1030         if (parse_errors() > 0) // fallback to default locale
1031             Yi = parse(qPrintable(func), "en_US");
1032         // DEBUG("  f(x["<< i <<"]) = " << Yi);
1033 
1034         if (parse_errors() > 0)
1035             return GSL_EINVAL;
1036 
1037         // DEBUG("  weight["<< i <<"]) = " << weight[i]);
1038         gsl_vector_set(f, i, sqrt(weight[i]) * (Yi - y[i]));
1039     }
1040 
1041     return GSL_SUCCESS;
1042 }
1043 
1044 /*!
1045  * calculates the matrix elements of Jacobian matrix
1046  * \param paramValues current parameter values
1047  * \param params
1048  * \param J Jacobian matrix
1049  * */
1050 int func_df(const gsl_vector* paramValues, void* params, gsl_matrix* J) {
1051     // DEBUG(Q_FUNC_INFO);
1052     const size_t n = ((struct data*)params)->n;
1053     double* xVector = ((struct data*)params)->x;
1054     double* weight = ((struct data*)params)->weight;
1055     auto modelCategory = ((struct data*)params)->modelCategory;
1056     unsigned int modelType = ((struct data*)params)->modelType;
1057     unsigned int degree = ((struct data*)params)->degree;
1058     auto* paramNames = ((struct data*)params)->paramNames;
1059     double* min = ((struct data*)params)->paramMin;
1060     double* max = ((struct data*)params)->paramMax;
1061     bool* fixed = ((struct data*)params)->paramFixed;
1062 
1063     // calculate the Jacobian matrix:
1064     // Jacobian matrix J(i,j) = df_i / dx_j
1065     // where f_i = w_i*(Y_i - y_i),
1066     // Y_i = model and the x_j are the parameters
1067     double x;
1068 
1069     switch (modelCategory) {
1070     case nsl_fit_model_basic:
1071         switch (modelType) {
1072         case nsl_fit_model_polynomial: // Y(x) = c0 + c1*x + ... + cn*x^n
1073             for (size_t i = 0; i < n; i++) {
1074                 x = xVector[i];
1075                 for (unsigned int j = 0; j < (unsigned int)paramNames->size(); ++j) {
1076                     if (fixed[j])
1077                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1078                     else
1079                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_polynomial_param_deriv(x, j, weight[i]));
1080                 }
1081             }
1082             break;
1083         case nsl_fit_model_power: // Y(x) = a*x^b or Y(x) = a + b*x^c.
1084             if (degree == 1) {
1085                 const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1086                 const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1087                 for (size_t i = 0; i < n; i++) {
1088                     x = xVector[i];
1089 
1090                     for (int j = 0; j < 2; j++) {
1091                         if (fixed[j])
1092                             gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1093                         else
1094                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_power1_param_deriv(j, x, a, b, weight[i]));
1095                     }
1096                 }
1097             } else if (degree == 2) {
1098                 const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1099                 const double c = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1100                 for (size_t i = 0; i < n; i++) {
1101                     x = xVector[i];
1102 
1103                     for (int j = 0; j < 3; j++) {
1104                         if (fixed[j])
1105                             gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1106                         else
1107                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_power2_param_deriv(j, x, b, c, weight[i]));
1108                     }
1109                 }
1110             }
1111             break;
1112         case nsl_fit_model_exponential: { // Y(x) = a*exp(b*x) + c*exp(d*x) + ...
1113             double* p = new double[2 * degree];
1114             for (unsigned int i = 0; i < 2 * degree; i++)
1115                 p[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, i), min[i], max[i]);
1116             for (size_t i = 0; i < n; i++) {
1117                 x = xVector[i];
1118 
1119                 for (unsigned int j = 0; j < 2 * degree; j++) {
1120                     if (fixed[j])
1121                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1122                     else
1123                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exponentialn_param_deriv(j, x, p, weight[i]));
1124                 }
1125             }
1126             delete[] p;
1127 
1128             break;
1129         }
1130         case nsl_fit_model_inverse_exponential: { // Y(x) = a*(1-exp(b*x))+c
1131             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1132             const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1133             for (size_t i = 0; i < n; i++) {
1134                 x = xVector[i];
1135 
1136                 for (unsigned int j = 0; j < 3; j++) {
1137                     if (fixed[j])
1138                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1139                     else
1140                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_inverse_exponential_param_deriv(j, x, a, b, weight[i]));
1141                 }
1142             }
1143             break;
1144         }
1145         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)
1146             // parameters: w, a0, a1, b1, ... an, bn
1147             double* a = new double[degree];
1148             double* b = new double[degree];
1149             double w = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1150             a[0] = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1151             b[0] = 0;
1152             for (unsigned int i = 1; i < degree; ++i) {
1153                 a[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, 2 * i), min[2 * i], max[2 * i]);
1154                 b[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, 2 * i + 1), min[2 * i + 1], max[2 * i + 1]);
1155             }
1156             for (size_t i = 0; i < n; i++) {
1157                 x = xVector[i];
1158                 double wd = 0; // first derivative with respect to the w parameter
1159                 for (unsigned int j = 1; j < degree; ++j) {
1160                     wd += -a[j] * j * x * sin(j * w * x) + b[j] * j * x * cos(j * w * x);
1161                 }
1162 
1163                 gsl_matrix_set(J, i, 0, weight[i] * wd);
1164                 gsl_matrix_set(J, i, 1, weight[i]);
1165                 for (unsigned int j = 1; j <= degree; ++j) {
1166                     gsl_matrix_set(J, (size_t)i, (size_t)(2 * j), nsl_fit_model_fourier_param_deriv(0, j, x, w, weight[i]));
1167                     gsl_matrix_set(J, (size_t)i, (size_t)(2 * j + 1), nsl_fit_model_fourier_param_deriv(1, j, x, w, weight[i]));
1168                 }
1169 
1170                 for (unsigned int j = 0; j <= 2 * degree + 1; j++)
1171                     if (fixed[j])
1172                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1173             }
1174 
1175             delete[] a;
1176             delete[] b;
1177 
1178             break;
1179         }
1180         }
1181         break;
1182     case nsl_fit_model_peak:
1183         switch (modelType) {
1184         case nsl_fit_model_gaussian:
1185         case nsl_fit_model_lorentz:
1186         case nsl_fit_model_sech:
1187         case nsl_fit_model_logistic:
1188             for (size_t i = 0; i < n; i++) {
1189                 x = xVector[i];
1190 
1191                 for (unsigned int j = 0; j < degree; ++j) {
1192                     const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 3 * j), min[3 * j], max[3 * j]);
1193                     const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 3 * j + 1), min[3 * j + 1], max[3 * j + 1]);
1194                     const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3 * j + 2), min[3 * j + 2], max[3 * j + 2]);
1195 
1196                     switch (modelType) {
1197                     case nsl_fit_model_gaussian:
1198                         gsl_matrix_set(J, (size_t)i, (size_t)(3 * j), nsl_fit_model_gaussian_param_deriv(0, x, a, s, mu, weight[i]));
1199                         gsl_matrix_set(J, (size_t)i, (size_t)(3 * j + 1), nsl_fit_model_gaussian_param_deriv(1, x, a, s, mu, weight[i]));
1200                         gsl_matrix_set(J, (size_t)i, (size_t)(3 * j + 2), nsl_fit_model_gaussian_param_deriv(2, x, a, s, mu, weight[i]));
1201                         break;
1202                     case nsl_fit_model_lorentz: // a,s,t
1203                         gsl_matrix_set(J, (size_t)i, (size_t)(3 * j), nsl_fit_model_lorentz_param_deriv(0, x, a, s, mu, weight[i]));
1204                         gsl_matrix_set(J, (size_t)i, (size_t)(3 * j + 1), nsl_fit_model_lorentz_param_deriv(1, x, a, s, mu, weight[i]));
1205                         gsl_matrix_set(J, (size_t)i, (size_t)(3 * j + 2), nsl_fit_model_lorentz_param_deriv(2, x, a, s, mu, weight[i]));
1206                         break;
1207                     case nsl_fit_model_sech:
1208                         gsl_matrix_set(J, (size_t)i, (size_t)(3 * j), nsl_fit_model_sech_param_deriv(0, x, a, s, mu, weight[i]));
1209                         gsl_matrix_set(J, (size_t)i, (size_t)(3 * j + 1), nsl_fit_model_sech_param_deriv(1, x, a, s, mu, weight[i]));
1210                         gsl_matrix_set(J, (size_t)i, (size_t)(3 * j + 2), nsl_fit_model_sech_param_deriv(2, x, a, s, mu, weight[i]));
1211                         break;
1212                     case nsl_fit_model_logistic:
1213                         gsl_matrix_set(J, (size_t)i, (size_t)(3 * j), nsl_fit_model_logistic_param_deriv(0, x, a, s, mu, weight[i]));
1214                         gsl_matrix_set(J, (size_t)i, (size_t)(3 * j + 1), nsl_fit_model_logistic_param_deriv(1, x, a, s, mu, weight[i]));
1215                         gsl_matrix_set(J, (size_t)i, (size_t)(3 * j + 2), nsl_fit_model_logistic_param_deriv(2, x, a, s, mu, weight[i]));
1216                         break;
1217                     }
1218                 }
1219 
1220                 for (unsigned int j = 0; j < 3 * degree; j++)
1221                     if (fixed[j])
1222                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1223             }
1224             break;
1225         case nsl_fit_model_voigt:
1226             for (size_t i = 0; i < n; i++) {
1227                 x = xVector[i];
1228 
1229                 for (unsigned int j = 0; j < degree; ++j) {
1230                     const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 4 * j), min[4 * j], max[4 * j]);
1231                     const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 4 * j + 1), min[4 * j + 1], max[4 * j + 1]);
1232                     const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 4 * j + 2), min[4 * j + 2], max[4 * j + 2]);
1233                     const double g = nsl_fit_map_bound(gsl_vector_get(paramValues, 4 * j + 3), min[4 * j + 3], max[4 * j + 3]);
1234 
1235                     gsl_matrix_set(J, (size_t)i, (size_t)(4 * j), nsl_fit_model_voigt_param_deriv(0, x, a, mu, s, g, weight[i]));
1236                     gsl_matrix_set(J, (size_t)i, (size_t)(4 * j + 1), nsl_fit_model_voigt_param_deriv(1, x, a, mu, s, g, weight[i]));
1237                     gsl_matrix_set(J, (size_t)i, (size_t)(4 * j + 2), nsl_fit_model_voigt_param_deriv(2, x, a, mu, s, g, weight[i]));
1238                     gsl_matrix_set(J, (size_t)i, (size_t)(4 * j + 3), nsl_fit_model_voigt_param_deriv(3, x, a, mu, s, g, weight[i]));
1239                 }
1240                 for (unsigned int j = 0; j < 4 * degree; j++)
1241                     if (fixed[j])
1242                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1243             }
1244             break;
1245         case nsl_fit_model_pseudovoigt1:
1246             for (size_t i = 0; i < n; i++) {
1247                 x = xVector[i];
1248 
1249                 for (unsigned int j = 0; j < degree; ++j) {
1250                     const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 4 * j), min[4 * j], max[4 * j]);
1251                     const double eta = nsl_fit_map_bound(gsl_vector_get(paramValues, 4 * j + 1), min[4 * j + 1], max[4 * j + 1]);
1252                     const double w = nsl_fit_map_bound(gsl_vector_get(paramValues, 4 * j + 2), min[4 * j + 2], max[4 * j + 2]);
1253                     const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 4 * j + 3), min[4 * j + 3], max[4 * j + 3]);
1254 
1255                     gsl_matrix_set(J, (size_t)i, (size_t)(4 * j), nsl_fit_model_pseudovoigt1_param_deriv(0, x, a, eta, w, mu, weight[i]));
1256                     gsl_matrix_set(J, (size_t)i, (size_t)(4 * j + 1), nsl_fit_model_pseudovoigt1_param_deriv(1, x, a, eta, w, mu, weight[i]));
1257                     gsl_matrix_set(J, (size_t)i, (size_t)(4 * j + 2), nsl_fit_model_pseudovoigt1_param_deriv(2, x, a, eta, w, mu, weight[i]));
1258                     gsl_matrix_set(J, (size_t)i, (size_t)(4 * j + 3), nsl_fit_model_pseudovoigt1_param_deriv(3, x, a, eta, w, mu, weight[i]));
1259                 }
1260                 for (unsigned int j = 0; j < 4 * degree; j++)
1261                     if (fixed[j])
1262                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1263             }
1264             break;
1265         }
1266         break;
1267     case nsl_fit_model_growth: {
1268         const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1269         const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1270         const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1271 
1272         for (size_t i = 0; i < n; i++) {
1273             x = xVector[i];
1274 
1275             for (unsigned int j = 0; j < 3; j++) {
1276                 if (fixed[j]) {
1277                     gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1278                 } else {
1279                     switch (modelType) {
1280                     case nsl_fit_model_atan:
1281                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_atan_param_deriv(j, x, a, mu, s, weight[i]));
1282                         break;
1283                     case nsl_fit_model_tanh:
1284                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_tanh_param_deriv(j, x, a, mu, s, weight[i]));
1285                         break;
1286                     case nsl_fit_model_algebraic_sigmoid:
1287                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_algebraic_sigmoid_param_deriv(j, x, a, mu, s, weight[i]));
1288                         break;
1289                     case nsl_fit_model_sigmoid:
1290                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_sigmoid_param_deriv(j, x, a, mu, s, weight[i]));
1291                         break;
1292                     case nsl_fit_model_erf:
1293                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_erf_param_deriv(j, x, a, mu, s, weight[i]));
1294                         break;
1295                     case nsl_fit_model_hill:
1296                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_hill_param_deriv(j, x, a, mu, s, weight[i]));
1297                         break;
1298                     case nsl_fit_model_gompertz:
1299                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gompertz_param_deriv(j, x, a, mu, s, weight[i]));
1300                         break;
1301                     case nsl_fit_model_gudermann:
1302                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gudermann_param_deriv(j, x, a, mu, s, weight[i]));
1303                         break;
1304                     }
1305                 }
1306             }
1307         }
1308     } break;
1309     case nsl_fit_model_distribution:
1310         switch (modelType) {
1311         case nsl_sf_stats_gaussian:
1312         case nsl_sf_stats_exponential:
1313         case nsl_sf_stats_laplace:
1314         case nsl_sf_stats_cauchy_lorentz:
1315         case nsl_sf_stats_rayleigh_tail:
1316         case nsl_sf_stats_lognormal:
1317         case nsl_sf_stats_logistic:
1318         case nsl_sf_stats_sech:
1319         case nsl_sf_stats_levy: {
1320             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1321             const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1322             const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1323             for (size_t i = 0; i < n; i++) {
1324                 x = xVector[i];
1325 
1326                 for (unsigned int j = 0; j < 3; j++) {
1327                     if (fixed[j])
1328                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1329                     else {
1330                         switch (modelType) {
1331                         case nsl_sf_stats_gaussian:
1332                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gaussian_param_deriv(j, x, a, s, mu, weight[i]));
1333                             break;
1334                         case nsl_sf_stats_exponential:
1335                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exponential_param_deriv(j, x, a, s, mu, weight[i]));
1336                             break;
1337                         case nsl_sf_stats_laplace:
1338                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_laplace_param_deriv(j, x, a, s, mu, weight[i]));
1339                             break;
1340                         case nsl_sf_stats_cauchy_lorentz:
1341                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_lorentz_param_deriv(j, x, a, s, mu, weight[i]));
1342                             break;
1343                         case nsl_sf_stats_rayleigh_tail:
1344                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_rayleigh_tail_param_deriv(j, x, a, s, mu, weight[i]));
1345                             break;
1346                         case nsl_sf_stats_lognormal:
1347                             if (x > 0)
1348                                 gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_lognormal_param_deriv(j, x, a, s, mu, weight[i]));
1349                             else
1350                                 gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1351                             break;
1352                         case nsl_sf_stats_logistic:
1353                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_logistic_param_deriv(j, x, a, s, mu, weight[i]));
1354                             break;
1355                         case nsl_sf_stats_sech:
1356                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_sech_dist_param_deriv(j, x, a, s, mu, weight[i]));
1357                             break;
1358                         case nsl_sf_stats_levy:
1359                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_levy_param_deriv(j, x, a, s, mu, weight[i]));
1360                             break;
1361                         }
1362                     }
1363                 }
1364             }
1365             break;
1366         }
1367         case nsl_sf_stats_gaussian_tail: {
1368             const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1369             const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1370             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1371             const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1372             for (size_t i = 0; i < n; i++) {
1373                 x = xVector[i];
1374 
1375                 for (unsigned int j = 0; j < 4; j++) {
1376                     if (fixed[j])
1377                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1378                     else
1379                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gaussian_tail_param_deriv(j, x, A, s, a, mu, weight[i]));
1380                 }
1381             }
1382             break;
1383         }
1384         case nsl_sf_stats_exponential_power: {
1385             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1386             const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1387             const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1388             const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1389             for (size_t i = 0; i < n; i++) {
1390                 x = xVector[i];
1391 
1392                 for (unsigned int j = 0; j < 4; j++) {
1393                     if (fixed[j])
1394                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1395                     else
1396                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exp_pow_param_deriv(j, x, a, s, b, mu, weight[i]));
1397                 }
1398             }
1399             break;
1400         }
1401         case nsl_sf_stats_rayleigh: {
1402             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1403             const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1404             for (size_t i = 0; i < n; i++) {
1405                 x = xVector[i];
1406 
1407                 for (unsigned int j = 0; j < 2; j++) {
1408                     if (fixed[j])
1409                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1410                     else
1411                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_rayleigh_param_deriv(j, x, a, s, weight[i]));
1412                 }
1413             }
1414             break;
1415         }
1416         case nsl_sf_stats_gamma: {
1417             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1418             const double k = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1419             const double t = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1420             for (size_t i = 0; i < n; i++) {
1421                 x = xVector[i];
1422 
1423                 for (unsigned int j = 0; j < 3; j++) {
1424                     if (fixed[j])
1425                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1426                     else
1427                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gamma_param_deriv(j, x, a, k, t, weight[i]));
1428                 }
1429             }
1430             break;
1431         }
1432         case nsl_sf_stats_flat: {
1433             const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1434             const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1435             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1436             for (size_t i = 0; i < n; i++) {
1437                 x = xVector[i];
1438 
1439                 for (unsigned int j = 0; j < 3; j++) {
1440                     if (fixed[j])
1441                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1442                     else
1443                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_flat_param_deriv(j, x, A, b, a, weight[i]));
1444                 }
1445             }
1446             break;
1447         }
1448         case nsl_sf_stats_chi_squared: {
1449             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1450             const double nu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1451             for (size_t i = 0; i < n; i++) {
1452                 x = xVector[i];
1453 
1454                 for (unsigned int j = 0; j < 2; j++) {
1455                     if (fixed[j])
1456                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1457                     else
1458                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_chi_square_param_deriv(j, x, a, nu, weight[i]));
1459                 }
1460             }
1461             break;
1462         }
1463         case nsl_sf_stats_tdist: {
1464             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1465             const double nu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1466             for (size_t i = 0; i < n; i++) {
1467                 x = xVector[i];
1468 
1469                 for (unsigned int j = 0; j < 2; j++) {
1470                     if (fixed[j])
1471                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1472                     else
1473                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_students_t_param_deriv(j, x, a, nu, weight[i]));
1474                 }
1475             }
1476             break;
1477         }
1478         case nsl_sf_stats_fdist: {
1479             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1480             const double n1 = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1481             const double n2 = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1482             for (size_t i = 0; i < n; i++) {
1483                 x = xVector[i];
1484 
1485                 for (unsigned int j = 0; j < 3; j++) {
1486                     if (fixed[j])
1487                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1488                     else
1489                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_fdist_param_deriv(j, x, a, n1, n2, weight[i]));
1490                 }
1491             }
1492             break;
1493         }
1494         case nsl_sf_stats_beta:
1495         case nsl_sf_stats_pareto: {
1496             const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1497             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1498             const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1499             for (size_t i = 0; i < n; i++) {
1500                 x = xVector[i];
1501 
1502                 for (unsigned int j = 0; j < 3; j++) {
1503                     if (fixed[j])
1504                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1505                     else {
1506                         switch (modelType) {
1507                         case nsl_sf_stats_beta:
1508                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_beta_param_deriv(j, x, A, a, b, weight[i]));
1509                             break;
1510                         case nsl_sf_stats_pareto:
1511                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_pareto_param_deriv(j, x, A, a, b, weight[i]));
1512                             break;
1513                         }
1514                     }
1515                 }
1516             }
1517             break;
1518         }
1519         case nsl_sf_stats_weibull: {
1520             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1521             const double k = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1522             const double l = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1523             const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1524             for (size_t i = 0; i < n; i++) {
1525                 x = xVector[i];
1526 
1527                 for (unsigned int j = 0; j < 4; j++) {
1528                     if (fixed[j])
1529                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1530                     else {
1531                         if (x > 0)
1532                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_weibull_param_deriv(j, x, a, k, l, mu, weight[i]));
1533                         else
1534                             gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1535                     }
1536                 }
1537             }
1538             break;
1539         }
1540         case nsl_sf_stats_gumbel1: {
1541             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1542             const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1543             const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1544             const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1545             for (size_t i = 0; i < n; i++) {
1546                 x = xVector[i];
1547 
1548                 for (unsigned int j = 0; j < 4; j++) {
1549                     if (fixed[j])
1550                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1551                     else
1552                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gumbel1_param_deriv(j, x, a, s, mu, b, weight[i]));
1553                 }
1554             }
1555             break;
1556         }
1557         case nsl_sf_stats_gumbel2: {
1558             const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1559             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1560             const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1561             const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1562             for (size_t i = 0; i < n; i++) {
1563                 x = xVector[i];
1564 
1565                 for (unsigned int j = 0; j < 4; j++) {
1566                     if (fixed[j])
1567                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1568                     else
1569                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gumbel2_param_deriv(j, x, A, a, b, mu, weight[i]));
1570                 }
1571             }
1572             break;
1573         }
1574         case nsl_sf_stats_poisson: {
1575             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1576             const double l = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1577             for (size_t i = 0; i < n; i++) {
1578                 x = xVector[i];
1579 
1580                 for (unsigned int j = 0; j < 2; j++) {
1581                     if (fixed[j])
1582                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1583                     else
1584                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_poisson_param_deriv(j, x, a, l, weight[i]));
1585                 }
1586             }
1587             break;
1588         }
1589         case nsl_sf_stats_maxwell_boltzmann: { // Y(x) = a*sqrt(2/pi) * x^2/s^3 * exp(-(x/s)^2/2)
1590             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1591             const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1592             for (size_t i = 0; i < n; i++) {
1593                 x = xVector[i];
1594 
1595                 for (unsigned int j = 0; j < 2; j++) {
1596                     if (fixed[j])
1597                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1598                     else
1599                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_maxwell_param_deriv(j, x, a, s, weight[i]));
1600                 }
1601             }
1602             break;
1603         }
1604         case nsl_sf_stats_frechet: {
1605             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1606             const double g = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1607             const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1608             const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1609             for (size_t i = 0; i < n; i++) {
1610                 x = xVector[i];
1611 
1612                 for (unsigned int j = 0; j < 4; j++) {
1613                     if (fixed[j])
1614                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1615                     else
1616                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_frechet_param_deriv(j, x, a, g, s, mu, weight[i]));
1617                 }
1618             }
1619             break;
1620         }
1621         case nsl_sf_stats_landau: {
1622             // const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1623             for (size_t i = 0; i < n; i++) {
1624                 x = xVector[i];
1625                 if (fixed[0])
1626                     gsl_matrix_set(J, (size_t)i, 0, 0.);
1627                 else
1628                     gsl_matrix_set(J, (size_t)i, 0, nsl_fit_model_landau_param_deriv(0, x, weight[i]));
1629             }
1630             break;
1631         }
1632         case nsl_sf_stats_binomial:
1633         case nsl_sf_stats_negative_binomial:
1634         case nsl_sf_stats_pascal: {
1635             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1636             const double p = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1637             const double N = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1638             for (size_t i = 0; i < n; i++) {
1639                 x = xVector[i];
1640 
1641                 for (unsigned int j = 0; j < 3; j++) {
1642                     if (fixed[j])
1643                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1644                     else {
1645                         switch (modelType) {
1646                         case nsl_sf_stats_binomial:
1647                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_binomial_param_deriv(j, x, a, p, N, weight[i]));
1648                             break;
1649                         case nsl_sf_stats_negative_binomial:
1650                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_negative_binomial_param_deriv(j, x, a, p, N, weight[i]));
1651                             break;
1652                         case nsl_sf_stats_pascal:
1653                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_pascal_param_deriv(j, x, a, p, N, weight[i]));
1654                             break;
1655                         }
1656                     }
1657                 }
1658             }
1659             break;
1660         }
1661         case nsl_sf_stats_geometric:
1662         case nsl_sf_stats_logarithmic: {
1663             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1664             const double p = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1665             for (size_t i = 0; i < n; i++) {
1666                 x = xVector[i];
1667 
1668                 for (unsigned int j = 0; j < 2; j++) {
1669                     if (fixed[j])
1670                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1671                     else {
1672                         switch (modelType) {
1673                         case nsl_sf_stats_geometric:
1674                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_geometric_param_deriv(j, x, a, p, weight[i]));
1675                             break;
1676                         case nsl_sf_stats_logarithmic:
1677                             gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_logarithmic_param_deriv(j, x, a, p, weight[i]));
1678                             break;
1679                         }
1680                     }
1681                 }
1682             }
1683             break;
1684         }
1685         case nsl_sf_stats_hypergeometric: {
1686             const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1687             const double n1 = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1688             const double n2 = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1689             const double t = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1690             for (size_t i = 0; i < n; i++) {
1691                 x = xVector[i];
1692 
1693                 for (unsigned int j = 0; j < 4; j++) {
1694                     if (fixed[j])
1695                         gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1696                     else
1697                         gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_hypergeometric_param_deriv(j, x, a, n1, n2, t, weight[i]));
1698                 }
1699             }
1700             break;
1701         }
1702         // unused distributions
1703         case nsl_sf_stats_levy_alpha_stable:
1704         case nsl_sf_stats_levy_skew_alpha_stable:
1705         case nsl_sf_stats_bernoulli:
1706             break;
1707         }
1708         break;
1709     case nsl_fit_model_custom:
1710         double value;
1711         const auto np = paramNames->size();
1712         QString func{*(((struct data*)params)->func)};
1713 
1714         const auto numberLocale = QLocale();
1715         for (size_t i = 0; i < n; i++) {
1716             x = xVector[i];
1717             assign_symbol("x", x);
1718 
1719             for (auto j = 0; j < np; j++) {
1720                 for (auto k = 0; k < np; k++) {
1721                     if (k != j) {
1722                         value = nsl_fit_map_bound(gsl_vector_get(paramValues, k), min[k], max[k]);
1723                         assign_symbol(qPrintable(paramNames->at(k)), value);
1724                     }
1725                 }
1726 
1727                 value = nsl_fit_map_bound(gsl_vector_get(paramValues, j), min[j], max[j]);
1728                 assign_symbol(qPrintable(paramNames->at(j)), value);
1729                 double f_p = parse(qPrintable(func), qPrintable(numberLocale.name()));
1730                 if (parse_errors() > 0) // fallback to default locale
1731                     f_p = parse(qPrintable(func), "en_US");
1732 
1733                 double eps = 1.e-9;
1734                 if (std::abs(f_p) > 0)
1735                     eps *= std::abs(f_p); // scale step size with function value
1736                 value += eps;
1737                 assign_symbol(qPrintable(paramNames->at(j)), value);
1738                 double f_pdp = parse(qPrintable(func), qPrintable(numberLocale.name()));
1739                 if (parse_errors() > 0) // fallback to default locale
1740                     f_pdp = parse(qPrintable(func), "en_US");
1741 
1742                 //              DEBUG("evaluate deriv"<<func<<": f(x["<<i<<"]) ="<<QString::number(f_p, 'g', 15));
1743                 //              DEBUG("evaluate deriv"<<func<<": f(x["<<i<<"]+dx) ="<<QString::number(f_pdp, 'g', 15));
1744                 //              DEBUG(" deriv = " << STDSTRING(QString::number(sqrt(weight[i])*(f_pdp-f_p)/eps, 'g', 15));
1745 
1746                 if (fixed[j])
1747                     gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1748                 else // calculate finite difference
1749                     gsl_matrix_set(J, (size_t)i, (size_t)j, sqrt(weight[i]) * (f_pdp - f_p) / eps);
1750             }
1751         }
1752     }
1753 
1754     return GSL_SUCCESS;
1755 }
1756 
1757 int func_fdf(const gsl_vector* x, void* params, gsl_vector* f, gsl_matrix* J) {
1758     // DEBUG(Q_FUNC_INFO);
1759     func_f(x, params, f);
1760     func_df(x, params, J);
1761 
1762     return GSL_SUCCESS;
1763 }
1764 
1765 /* prepare the fit result columns */
1766 void XYFitCurvePrivate::prepareResultColumns() {
1767     // DEBUG(Q_FUNC_INFO)
1768     // create fit result columns if not available yet, clear them otherwise
1769 
1770     // Done also in XYAnalysisCurve, but this function will be also called directly() from evaluate()
1771     // and not over recalculate(). So this is also needed here!
1772     if (!xColumn) { // all columns are treated together
1773         DEBUG(" Creating columns")
1774         xColumn = new Column(QStringLiteral("x"), AbstractColumn::ColumnMode::Double);
1775         yColumn = new Column(QStringLiteral("y"), AbstractColumn::ColumnMode::Double);
1776 
1777         xVector = static_cast<QVector<double>*>(xColumn->data());
1778         yVector = static_cast<QVector<double>*>(yColumn->data());
1779 
1780         xColumn->setHidden(true);
1781         q->addChild(xColumn);
1782 
1783         yColumn->setHidden(true);
1784         q->addChild(yColumn);
1785 
1786         q->setUndoAware(false);
1787         q->setXColumn(xColumn);
1788         q->setYColumn(yColumn);
1789         q->setUndoAware(true);
1790     } else {
1791         DEBUG(Q_FUNC_INFO << ", Clear columns")
1792         xColumn->invalidateProperties();
1793         yColumn->invalidateProperties();
1794         xVector->clear();
1795         yVector->clear();
1796         // TODO: residualsVector->clear();
1797     }
1798 
1799     if (!residualsColumn) {
1800         residualsColumn = new Column(QStringLiteral("residuals"), AbstractColumn::ColumnMode::Double);
1801         residualsVector = static_cast<QVector<double>*>(residualsColumn->data());
1802         residualsColumn->setFixed(true); // visible in the project explorer but cannot be modified (renamed, deleted, etc.)
1803         q->addChild(residualsColumn);
1804     }
1805 }
1806 
1807 void XYFitCurvePrivate::resetResults() {
1808     fitResult = XYFitCurve::FitResult();
1809 }
1810 
1811 void XYFitCurvePrivate::prepareTmpDataColumn(const AbstractColumn** tmpXDataColumn, const AbstractColumn** tmpYDataColumn) {
1812     // prepare source data columns
1813     DEBUG(Q_FUNC_INFO << ", data source: " << ENUM_TO_STRING(XYAnalysisCurve, DataSourceType, dataSourceType))
1814     switch (dataSourceType) {
1815     case XYAnalysisCurve::DataSourceType::Spreadsheet:
1816         *tmpXDataColumn = xDataColumn;
1817         *tmpYDataColumn = yDataColumn;
1818         break;
1819     case XYAnalysisCurve::DataSourceType::Curve:
1820         *tmpXDataColumn = dataSourceCurve->xColumn();
1821         *tmpYDataColumn = dataSourceCurve->yColumn();
1822         break;
1823     case XYAnalysisCurve::DataSourceType::Histogram:
1824         switch (fitData.algorithm) {
1825         case nsl_fit_algorithm_lm:
1826             *tmpXDataColumn = dataSourceHistogram->bins(); // bins
1827             switch (dataSourceHistogram->normalization()) { // TODO: not exactly
1828             case Histogram::Normalization::Count:
1829             case Histogram::Normalization::CountDensity:
1830                 *tmpYDataColumn = dataSourceHistogram->binValues(); // values
1831                 break;
1832             case Histogram::Normalization::Probability:
1833             case Histogram::Normalization::ProbabilityDensity:
1834                 *tmpYDataColumn = dataSourceHistogram->binPDValues(); // normalized values
1835             }
1836             break;
1837         case nsl_fit_algorithm_ml:
1838             *tmpXDataColumn = dataSourceHistogram->dataColumn(); // data
1839             *tmpYDataColumn = dataSourceHistogram->binPDValues(); // normalized values
1840         }
1841         // debug
1842         /*for (int i = 0; i < dataSourceHistogram->bins()->rowCount(); i++)
1843             DEBUG("BINS @ " << i << ": " << dataSourceHistogram->bins()->valueAt(i))
1844         for (int i = 0; i < dataSourceHistogram->binValues()->rowCount(); i++)
1845             DEBUG("BINValues @ " << i << ": " << dataSourceHistogram->binValues()->valueAt(i))
1846         for (int i = 0; i < dataSourceHistogram->binPDValues()->rowCount(); i++)
1847             DEBUG("BINPDValues @ " << i << ": " << dataSourceHistogram->binPDValues()->valueAt(i))
1848         */
1849     }
1850 }
1851 
1852 bool XYFitCurvePrivate::recalculateSpecific(const AbstractColumn* tmpXDataColumn, const AbstractColumn* tmpYDataColumn) {
1853     QElapsedTimer timer;
1854     timer.start();
1855 
1856     // determine range of data
1857     Range<double> xRange{tmpXDataColumn->minimum(), tmpXDataColumn->maximum()};
1858     if (fitData.autoRange) { // auto x range of data to fit
1859         fitData.fitRange = xRange;
1860     } else { // custom x range of data to fit
1861         if (!fitData.fitRange.isZero()) // avoid problems with user specified zero range
1862             xRange.setRange(fitData.fitRange.start(), fitData.fitRange.end());
1863     }
1864     DEBUG(Q_FUNC_INFO << ", fit data range = " << xRange.start() << " .. " << xRange.end());
1865 
1866     prepareResultColumns();
1867     const size_t rowCount = tmpXDataColumn->rowCount();
1868 
1869     // fill residuals vector. To get residuals on the correct x values, fill the rest with zeros.
1870     residualsVector->resize(rowCount);
1871     // DEBUG("  Residual vector size: " << residualsVector->size())
1872 
1873     DEBUG("#######################################\nALGORITHM: " << nsl_fit_algorithm_name[fitData.algorithm])
1874     switch (fitData.algorithm) {
1875     case nsl_fit_algorithm_lm:
1876         runLevenbergMarquardt(tmpXDataColumn, tmpYDataColumn, xRange);
1877         break;
1878     case nsl_fit_algorithm_ml: {
1879         double width = xRange.size() / tmpYDataColumn->rowCount();
1880         double norm = 1.;
1881         if (dataSourceHistogram) {
1882             switch (dataSourceHistogram->normalization()) {
1883             case Histogram::Normalization::Count:
1884                 norm = rowCount * width;
1885                 break;
1886             case Histogram::Normalization::Probability:
1887                 norm = width;
1888                 break;
1889             case Histogram::Normalization::CountDensity:
1890                 norm = rowCount;
1891                 break;
1892             case Histogram::Normalization::ProbabilityDensity:
1893                 break;
1894             }
1895         } else { // spreadsheet or curve
1896             norm = ((Column*)tmpYDataColumn)->statistics().arithmeticMean * xRange.size(); // integral
1897         }
1898         runMaximumLikelihood(tmpXDataColumn, norm);
1899     }
1900     }
1901 
1902     const bool update = evaluate(); // calculate the fit function (vectors)
1903 
1904     // ML uses dataSourceHistogram->bins() as x for residuals
1905     if (dataSourceType == XYAnalysisCurve::DataSourceType::Histogram && dataSourceHistogram && fitData.algorithm == nsl_fit_algorithm_ml)
1906         tmpXDataColumn = dataSourceHistogram->bins();
1907 
1908     if (fitData.autoRange || fitData.algorithm == nsl_fit_algorithm_ml) { // evaluate residuals
1909         QVector<double> v;
1910         v.resize(rowCount);
1911         for (size_t i = 0; i < rowCount; i++)
1912             if (tmpXDataColumn->isNumeric())
1913                 v[i] = tmpXDataColumn->valueAt(i);
1914             else if (tmpXDataColumn->columnMode() == AbstractColumn::ColumnMode::DateTime)
1915                 v[i] = tmpXDataColumn->dateTimeAt(i).toMSecsSinceEpoch();
1916 
1917         auto* parser = ExpressionParser::getInstance();
1918         // fill residualsVector with model values
1919         // QDEBUG("xVector: " << v)
1920         bool rc = parser->evaluateCartesian(fitData.model, &v, residualsVector, fitData.paramNames, fitResult.paramValues);
1921         // QDEBUG("residualsVector: " << *residualsVector)
1922         if (rc) {
1923             switch (fitData.algorithm) {
1924             case nsl_fit_algorithm_lm:
1925                 for (size_t i = 0; i < rowCount; i++)
1926                     (*residualsVector)[i] = tmpYDataColumn->valueAt(i) - (*residualsVector).at(i);
1927                 break;
1928             case nsl_fit_algorithm_ml:
1929                 for (size_t i = 0; i < rowCount; i++) {
1930                     // DEBUG("y data / column @" << i << ":" << tmpXDataColumn->valueAt(i))
1931                     if (xRange.contains(tmpXDataColumn->valueAt(i)))
1932                         (*residualsVector)[i] = tmpYDataColumn->valueAt(i) - (*residualsVector).at(i);
1933                     else
1934                         (*residualsVector)[i] = 0.;
1935                 }
1936             }
1937         } else {
1938             WARN("  ERROR: Failed parsing residuals")
1939             residualsVector->clear();
1940         }
1941     } // else: see LM method
1942 
1943     residualsColumn->setChanged();
1944 
1945     fitResult.elapsedTime = timer.elapsed();
1946 
1947     return update;
1948 }
1949 
1950 void XYFitCurvePrivate::runMaximumLikelihood(const AbstractColumn* tmpXDataColumn, const double norm) {
1951     const size_t n = tmpXDataColumn->rowCount();
1952 
1953     fitResult.available = true;
1954     fitResult.valid = true;
1955     fitResult.status = i18n("Success"); // can it fail in any way?
1956 
1957     const unsigned int np = fitData.paramNames.size(); // number of fit parameters
1958     fitResult.dof = n - np;
1959     fitResult.paramValues.resize(np);
1960     fitResult.errorValues.resize(np);
1961     fitResult.tdist_tValues.resize(np);
1962     fitResult.tdist_pValues.resize(np);
1963     fitResult.marginValues.resize(np);
1964     fitResult.correlationMatrix.resize(np * (np + 1) / 2);
1965 
1966     DEBUG("DISTRIBUTION: " << fitData.modelType)
1967     fitResult.paramValues[0] = norm; // A - normalization
1968     // TODO: parameter values (error, etc.)
1969     // TODO: currently all values are used (data range not changeable)
1970     const double alpha = 1.0 - fitData.confidenceInterval / 100.;
1971     const auto statistics = ((Column*)tmpXDataColumn)->statistics();
1972     const double mean = statistics.arithmeticMean;
1973     const double var = statistics.variance;
1974     const double median = statistics.median;
1975     const double madmed = statistics.meanDeviationAroundMedian;
1976     const double iqr = statistics.iqr;
1977     switch (fitData.modelType) { // only these are supported
1978     case nsl_sf_stats_gaussian: {
1979         const double sigma = std::sqrt(var);
1980         const double mu = mean;
1981         fitResult.paramValues[1] = sigma;
1982         fitResult.paramValues[2] = mu;
1983         DEBUG("mu = " << mu << ", sigma = " << sigma)
1984 
1985         fitResult.errorValues[2] = sigma / std::sqrt(n);
1986         double margin = nsl_stats_tdist_margin(alpha, fitResult.dof, fitResult.errorValues.at(2));
1987         // DEBUG("z = " << nsl_stats_tdist_z(alpha, fitResult.dof))
1988         fitResult.marginValues[2] = margin;
1989 
1990         fitResult.errorValues[1] = sigma * sigma / std::sqrt(2 * n);
1991         margin = nsl_stats_tdist_margin(alpha, fitResult.dof, fitResult.errorValues.at(1));
1992         // WARN("sigma CONFIDENCE INTERVAL: " << fitResult.paramValues[1] - margin << " .. " << fitResult.paramValues[1] + margin)
1993         fitResult.marginValues[1] = margin;
1994 
1995         // normalization for spreadsheet or curve
1996         if (dataSourceType != XYAnalysisCurve::DataSourceType::Histogram)
1997             fitResult.paramValues[0] /=
1998                 (gsl_sf_erf((tmpXDataColumn->maximum() - mu) / sigma) - gsl_sf_erf((tmpXDataColumn->minimum() - mu) / sigma)) / (2. * std::sqrt(2.));
1999         break;
2000     }
2001     case nsl_sf_stats_exponential: {
2002         const double mu = tmpXDataColumn->minimum();
2003         const double lambda = 1. / (mean - mu); // 1/(<x>-\mu)
2004         fitResult.paramValues[1] = lambda * (1 - 1. / (n - 1)); // unbiased
2005         fitResult.paramValues[2] = mu;
2006 
2007         fitResult.errorValues[1] = lambda / std::sqrt(n);
2008         // exact method
2009         double margin = lambda * (1. - gsl_cdf_chisq_Pinv(alpha / 2., 2 * n) / (2. * n));
2010         fitResult.marginValues[1] = margin;
2011         margin = lambda * (gsl_cdf_chisq_Pinv(1. - alpha / 2., 2 * n) / (2. * n) - 1.);
2012         fitResult.margin2Values.resize(2);
2013         fitResult.margin2Values[1] = margin;
2014 
2015         DEBUG("error = " << fitResult.errorValues.at(1))
2016         DEBUG("CI: " << gsl_cdf_chisq_Pinv(alpha / 2., 2 * n) / (2. * n) * lambda << " .. " << gsl_cdf_chisq_Pinv(1. - alpha / 2., 2 * n) / (2. * n) * lambda)
2017         DEBUG("normal approx.: " << lambda * (1. - 1.96 / std::sqrt(n)) << " .. " << lambda * (1. + 1.96 / std::sqrt(n)))
2018         DEBUG("1/l = " << 1. / fitResult.paramValues.at(1))
2019         DEBUG("1/l CI: " << 2. * n / gsl_cdf_chisq_Pinv(1. - alpha / 2., 2 * n) / lambda << " .. " << 2. * n / gsl_cdf_chisq_Pinv(alpha / 2., 2 * n) / lambda)
2020         DEBUG("1/l normal approx.: " << 1. / lambda / (1. + 1.96 / std::sqrt(n)) << " .. " << 1. / lambda / (1. - 1.96 / std::sqrt(n)))
2021 
2022         // normalization for spreadsheet or curve
2023         if (dataSourceType != XYAnalysisCurve::DataSourceType::Histogram)
2024             fitResult.paramValues[0] /= std::exp(-lambda * (tmpXDataColumn->minimum() - mu)) - std::exp(-lambda * (tmpXDataColumn->maximum() - mu));
2025         break;
2026     }
2027     case nsl_sf_stats_laplace: {
2028         double sigma = madmed;
2029         if (n > 2) // bias correction
2030             sigma *= n / (n - 2.);
2031         const double mu = median;
2032         fitResult.paramValues[1] = sigma;
2033         fitResult.paramValues[2] = mu;
2034 
2035         // TODO: error + CI
2036 
2037         // normalization for spreadsheet or curve
2038         if (dataSourceType != XYAnalysisCurve::DataSourceType::Histogram) {
2039             double Fmin, Fmax;
2040             const double xmin = tmpXDataColumn->minimum(), xmax = tmpXDataColumn->maximum();
2041             if (xmin < mu)
2042                 Fmin = .5 * std::exp((xmin - mu) / sigma);
2043             else
2044                 Fmin = 1. - .5 * std::exp(-(xmin - mu) / sigma);
2045             if (xmax < mu)
2046                 Fmax = .5 * std::exp((xmax - mu) / sigma);
2047             else
2048                 Fmax = 1. - .5 * std::exp(-(xmax - mu) / sigma);
2049 
2050             fitResult.paramValues[0] /= Fmax - Fmin;
2051         }
2052         break;
2053     }
2054     case nsl_sf_stats_cauchy_lorentz: {
2055         const double gamma = iqr / 2.;
2056         const double mu = median; // better truncated mean of middle 24%
2057         fitResult.paramValues[1] = gamma;
2058         fitResult.paramValues[2] = mu;
2059 
2060         // TODO: error + CI
2061 
2062         // normalization for spreadsheet or curve
2063         if (dataSourceType != XYAnalysisCurve::DataSourceType::Histogram)
2064             fitResult.paramValues[0] /= 1. / M_PI * (atan((tmpXDataColumn->maximum() - mu) / gamma) - atan((tmpXDataColumn->minimum() - mu) / gamma));
2065         break;
2066     }
2067     case nsl_sf_stats_lognormal: {
2068         // calculate mu and sigma
2069         double mu = 0.;
2070         for (size_t i = 0; i < n; i++)
2071             mu += std::log(tmpXDataColumn->valueAt(i));
2072         mu /= n;
2073         double var = 0.;
2074         for (size_t i = 0; i < n; i++)
2075             var += gsl_pow_2(std::log(tmpXDataColumn->valueAt(i)) - mu);
2076         var /= (n - 1);
2077         const double sigma = std::sqrt(var);
2078         fitResult.paramValues[1] = sigma;
2079         fitResult.paramValues[2] = mu;
2080 
2081         // TODO: error + CI
2082 
2083         // normalization for spreadsheet or curve
2084         if (dataSourceType != XYAnalysisCurve::DataSourceType::Histogram)
2085             fitResult.paramValues[0] /= gsl_sf_erf((std::log(tmpXDataColumn->maximum()) - mu) / sigma)
2086                 - gsl_sf_erf((std::log(tmpXDataColumn->minimum()) - mu) / sigma) / (2. * std::sqrt(2.));
2087         break;
2088     }
2089     case nsl_sf_stats_poisson: {
2090         const double lambda = mean;
2091         fitResult.paramValues[1] = lambda;
2092         fitResult.errorValues[1] = std::sqrt(lambda / n);
2093 
2094         // double margin = 1.96 * fitResult.errorValues[1]; // normal approx.
2095         // DEBUG("low / high = " << nsl_stats_chisq_low(alpha, n * lambda) << ", " << nsl_stats_chisq_high(alpha, n * lambda)) // exact formula
2096         fitResult.marginValues[1] = lambda - nsl_stats_chisq_low(alpha, n * lambda) / n;
2097         fitResult.margin2Values.resize(2);
2098         fitResult.margin2Values[1] = nsl_stats_chisq_high(alpha, n * lambda) / n - lambda;
2099 
2100         // normalization for spreadsheet or curve
2101         if (dataSourceType != XYAnalysisCurve::DataSourceType::Histogram)
2102             fitResult.paramValues[0] /=
2103                 gsl_sf_gamma_inc_Q(floor(tmpXDataColumn->maximum() + 1), lambda) - gsl_sf_gamma_inc_Q(floor(tmpXDataColumn->minimum() + 1), lambda);
2104         break;
2105     }
2106     case nsl_sf_stats_binomial: {
2107         const double p = mean / n;
2108         fitResult.paramValues[1] = p;
2109         fitResult.paramValues[2] = n;
2110         fitResult.errorValues[1] = std::sqrt(p * (1 - p) / n);
2111 
2112         // Clopper-Pearson exact method
2113         const double k = p * n;
2114         fitResult.marginValues[1] = (p - gsl_cdf_beta_Pinv(alpha / 2., k, n - k + 1)) / std::sqrt(n);
2115         fitResult.margin2Values.resize(2);
2116         fitResult.margin2Values[1] = (gsl_cdf_beta_Pinv(1. - alpha / 2., k + 1, n - k) - p) / std::sqrt(n);
2117 
2118         // normalization for spreadsheet or curve
2119         const double kmin = tmpXDataColumn->minimum(), kmax = tmpXDataColumn->maximum();
2120         if (dataSourceType != XYAnalysisCurve::DataSourceType::Histogram)
2121             fitResult.paramValues[0] /= gsl_sf_beta_inc(n - kmax, kmax + 1, 1 - p) - gsl_sf_beta_inc(n - kmin, kmin + 1, 1 - p);
2122     }
2123     }
2124 
2125     fitResult.calculateResult(n, np);
2126 
2127     if (fitData.useResults) // set start values
2128         for (unsigned int i = 0; i < np; i++)
2129             fitData.paramStartValues.data()[i] = fitResult.paramValues.at(i);
2130 }
2131 
2132 void XYFitCurvePrivate::runLevenbergMarquardt(const AbstractColumn* tmpXDataColumn, const AbstractColumn* tmpYDataColumn, const Range<double> xRange) {
2133     // fit settings
2134     const unsigned int maxIters = fitData.maxIterations; // maximal number of iterations
2135     const double delta = fitData.eps; // fit tolerance
2136     const auto np = fitData.paramNames.size(); // number of fit parameters
2137     if (np == 0) {
2138         DEBUG(Q_FUNC_INFO << ", WARNING: no parameter found.")
2139         fitResult.available = true;
2140         fitResult.valid = false;
2141         fitResult.status = i18n("Model has no parameters.");
2142         return;
2143     }
2144 
2145     if (yErrorColumn && yErrorColumn->rowCount() < tmpXDataColumn->rowCount()) {
2146         fitResult.available = true;
2147         fitResult.valid = false;
2148         fitResult.status = i18n("Not sufficient weight data points provided.");
2149         return;
2150     }
2151 
2152     // copy all valid data point for the fit to temporary vectors
2153     QVector<double> xdataVector;
2154     QVector<double> ydataVector;
2155     QVector<double> xerrorVector;
2156     QVector<double> yerrorVector;
2157 
2158     // logic from XYAnalysisCurve::copyData(), extended by the handling of error columns.
2159     // TODO: decide how to deal with non-numerical error columns
2160     int rowCount = std::min(tmpXDataColumn->rowCount(), tmpYDataColumn->rowCount());
2161     for (int row = 0; row < rowCount; ++row) {
2162         // omit invalid data
2163         if (!tmpXDataColumn->isValid(row) || tmpXDataColumn->isMasked(row) || !tmpYDataColumn->isValid(row) || tmpYDataColumn->isMasked(row))
2164             continue;
2165 
2166         double x = NAN;
2167         switch (tmpXDataColumn->columnMode()) {
2168         case AbstractColumn::ColumnMode::Double:
2169         case AbstractColumn::ColumnMode::Integer:
2170         case AbstractColumn::ColumnMode::BigInt:
2171             x = tmpXDataColumn->valueAt(row);
2172             break;
2173         case AbstractColumn::ColumnMode::Text: // not valid
2174             break;
2175         case AbstractColumn::ColumnMode::DateTime:
2176         case AbstractColumn::ColumnMode::Day:
2177         case AbstractColumn::ColumnMode::Month:
2178             x = tmpXDataColumn->dateTimeAt(row).toMSecsSinceEpoch();
2179         }
2180 
2181         double y = NAN;
2182         switch (tmpYDataColumn->columnMode()) {
2183         case AbstractColumn::ColumnMode::Double:
2184         case AbstractColumn::ColumnMode::Integer:
2185         case AbstractColumn::ColumnMode::BigInt:
2186             y = tmpYDataColumn->valueAt(row);
2187             break;
2188         case AbstractColumn::ColumnMode::Text: // not valid
2189             break;
2190         case AbstractColumn::ColumnMode::DateTime:
2191         case AbstractColumn::ColumnMode::Day:
2192         case AbstractColumn::ColumnMode::Month:
2193             y = tmpYDataColumn->dateTimeAt(row).toMSecsSinceEpoch();
2194         }
2195 
2196         if (x >= xRange.start() && x <= xRange.end()) { // only when inside given range
2197             if ((!xErrorColumn && !yErrorColumn) || !fitData.useDataErrors) { // x-y
2198                 xdataVector.append(x);
2199                 ydataVector.append(y);
2200             } else if (!xErrorColumn && yErrorColumn) { // x-y-dy
2201                 if (!std::isnan(yErrorColumn->valueAt(row))) {
2202                     xdataVector.append(x);
2203                     ydataVector.append(y);
2204                     yerrorVector.append(yErrorColumn->valueAt(row));
2205                 }
2206             } else if (xErrorColumn && yErrorColumn) { // x-y-dx-dy
2207                 if (!std::isnan(xErrorColumn->valueAt(row)) && !std::isnan(yErrorColumn->valueAt(row))) {
2208                     xdataVector.append(x);
2209                     ydataVector.append(y);
2210                     xerrorVector.append(xErrorColumn->valueAt(row));
2211                     yerrorVector.append(yErrorColumn->valueAt(row));
2212                 }
2213             }
2214         }
2215     }
2216 
2217     // QDEBUG(Q_FUNC_INFO << ", data: " << ydataVector)
2218 
2219     // number of data points to fit
2220     const auto n = xdataVector.size();
2221     DEBUG(Q_FUNC_INFO << ", number of data points: " << n);
2222     if (n == 0) {
2223         fitResult.available = true;
2224         fitResult.valid = false;
2225         fitResult.status = i18n("No data points available.");
2226         return;
2227     }
2228 
2229     if (n < np) {
2230         fitResult.available = true;
2231         fitResult.valid = false;
2232         fitResult.status = i18n("The number of data points (%1) must be greater than or equal to the number of parameters (%2).", n, np);
2233         return;
2234     }
2235 
2236     if (fitData.model.simplified().isEmpty()) {
2237         fitResult.available = true;
2238         fitResult.valid = false;
2239         fitResult.status = i18n("Fit model not specified.");
2240         return;
2241     }
2242 
2243     double* xdata = xdataVector.data();
2244     double* ydata = ydataVector.data();
2245     double* xerror = xerrorVector.data(); // size may be 0
2246     double* yerror = yerrorVector.data(); // size may be 0
2247     DEBUG(Q_FUNC_INFO << ", x error vector size: " << xerrorVector.size());
2248     DEBUG(Q_FUNC_INFO << ", y error vector size: " << yerrorVector.size());
2249     double* weight = new double[n];
2250 
2251     for (auto i = 0; i < n; i++)
2252         weight[i] = 1.;
2253 
2254     const double minError = 1.e-199; // minimum error for weighting
2255 
2256     switch (fitData.yWeightsType) {
2257     case nsl_fit_weight_no:
2258     case nsl_fit_weight_statistical_fit:
2259     case nsl_fit_weight_relative_fit:
2260         break;
2261     case nsl_fit_weight_instrumental: // yerror are sigmas
2262         for (int i = 0; i < (int)n; i++)
2263             if (i < yerrorVector.size())
2264                 weight[i] = 1. / gsl_pow_2(std::max(yerror[i], std::max(sqrt(minError), std::abs(ydata[i]) * 1.e-15)));
2265         break;
2266     case nsl_fit_weight_direct: // yerror are weights
2267         for (int i = 0; i < (int)n; i++)
2268             if (i < yerrorVector.size())
2269                 weight[i] = yerror[i];
2270         break;
2271     case nsl_fit_weight_inverse: // yerror are inverse weights
2272         for (int i = 0; i < (int)n; i++)
2273             if (i < yerrorVector.size())
2274                 weight[i] = 1. / std::max(yerror[i], std::max(minError, std::abs(ydata[i]) * 1.e-15));
2275         break;
2276     case nsl_fit_weight_statistical:
2277         for (int i = 0; i < (int)n; i++)
2278             weight[i] = 1. / std::max(ydata[i], minError);
2279         break;
2280     case nsl_fit_weight_relative:
2281         for (int i = 0; i < (int)n; i++)
2282             weight[i] = 1. / std::max(gsl_pow_2(ydata[i]), minError);
2283         break;
2284     }
2285 
2286     /////////////////////// GSL >= 2 has a complete new interface! But the old one is still supported. ///////////////////////////
2287     // 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"
2288     int nf = 0; // number of fixed parameter
2289     // size of paramFixed may be smaller than np
2290     if (fitData.paramFixed.size() < np)
2291         fitData.paramFixed.resize(np);
2292     for (auto i = 0; i < np; i++) {
2293         const bool fixed = fitData.paramFixed.at(i);
2294         if (fixed)
2295             nf++;
2296         DEBUG(" parameter " << i << " fixed: " << fixed);
2297     }
2298 
2299     // function to fit
2300     gsl_multifit_function_fdf f;
2301     DEBUG(Q_FUNC_INFO << ", model = " << STDSTRING(fitData.model));
2302     struct data params = {static_cast<size_t>(n),
2303                           xdata,
2304                           ydata,
2305                           weight,
2306                           fitData.modelCategory,
2307                           fitData.modelType,
2308                           fitData.degree,
2309                           &fitData.model,
2310                           &fitData.paramNames,
2311                           fitData.paramLowerLimits.data(),
2312                           fitData.paramUpperLimits.data(),
2313                           fitData.paramFixed.data()};
2314     f.f = &func_f;
2315     f.df = &func_df;
2316     f.fdf = &func_fdf;
2317     f.n = n;
2318     f.p = np;
2319     f.params = &params;
2320 
2321     DEBUG(Q_FUNC_INFO << ", initialize the derivative solver (using Levenberg-Marquardt robust solver)");
2322     const auto* T = gsl_multifit_fdfsolver_lmsder;
2323     auto* s = gsl_multifit_fdfsolver_alloc(T, n, np);
2324 
2325     DEBUG(Q_FUNC_INFO << ", set start values");
2326     double* x_init = fitData.paramStartValues.data();
2327     double* x_min = fitData.paramLowerLimits.data();
2328     double* x_max = fitData.paramUpperLimits.data();
2329     DEBUG(Q_FUNC_INFO << ", scale start values if limits are set");
2330     for (auto i = 0; i < np; i++)
2331         x_init[i] = nsl_fit_map_unbound(x_init[i], x_min[i], x_max[i]);
2332     DEBUG(Q_FUNC_INFO << ", DONE");
2333     auto x = gsl_vector_view_array(x_init, np);
2334     DEBUG(Q_FUNC_INFO << ", Turning off GSL error handler to avoid overflow/underflow");
2335     gsl_set_error_handler_off();
2336     DEBUG(Q_FUNC_INFO << ", Initialize solver with function f and initial guess x");
2337     gsl_multifit_fdfsolver_set(s, &f, &x.vector);
2338 
2339     DEBUG(Q_FUNC_INFO << ", Iterate ...");
2340     int status = GSL_SUCCESS;
2341     unsigned int iter = 0;
2342     fitResult.solverOutput.clear();
2343     writeSolverState(s);
2344     do {
2345         iter++;
2346         DEBUG(Q_FUNC_INFO << ", iter " << iter);
2347 
2348         // update weights for Y-depending weights (using function values from residuals)
2349         if (fitData.yWeightsType == nsl_fit_weight_statistical_fit) {
2350             for (auto i = 0; i < n; i++)
2351                 weight[i] = 1. / (gsl_vector_get(s->f, i) / sqrt(weight[i]) + ydata[i]); // 1/Y_i
2352         } else if (fitData.yWeightsType == nsl_fit_weight_relative_fit) {
2353             for (auto i = 0; i < n; i++)
2354                 weight[i] = 1. / gsl_pow_2(gsl_vector_get(s->f, i) / sqrt(weight[i]) + ydata[i]); // 1/Y_i^2
2355         }
2356 
2357         if (nf == np) { // all fixed parameter
2358             DEBUG(Q_FUNC_INFO << ", all parameter fixed. Stop iteration.")
2359             break;
2360         }
2361         DEBUG(Q_FUNC_INFO << ", run fdfsolver_iterate");
2362         status = gsl_multifit_fdfsolver_iterate(s);
2363         DEBUG(Q_FUNC_INFO << ", fdfsolver_iterate DONE");
2364         double chi = gsl_blas_dnrm2(s->f);
2365         writeSolverState(s, chi);
2366         if (status) {
2367             DEBUG(Q_FUNC_INFO << ", iter " << iter << ", status = " << gsl_strerror(status));
2368             if (status == GSL_ETOLX) // change in the position vector falls below machine precision: no progress
2369                 status = GSL_SUCCESS;
2370             break;
2371         }
2372         if (qFuzzyIsNull(chi)) {
2373             DEBUG(Q_FUNC_INFO << ", chi is zero! Finishing.")
2374             status = GSL_SUCCESS;
2375         } else {
2376             status = gsl_multifit_test_delta(s->dx, s->x, delta, delta);
2377         }
2378         DEBUG(Q_FUNC_INFO << ", iter " << iter << ", test status = " << gsl_strerror(status));
2379     } while (status == GSL_CONTINUE && iter < maxIters);
2380 
2381     // second run for x-error fitting
2382     if (xerrorVector.size() > 0) {
2383         DEBUG(Q_FUNC_INFO << ", Rerun fit with x errors");
2384 
2385         unsigned int iter2 = 0;
2386         double chi = 0, chiOld = 0;
2387         double* fun = new double[n];
2388         do {
2389             iter2++;
2390             chiOld = chi;
2391             // printf("iter2 = %d\n", iter2);
2392 
2393             // calculate function from residuals
2394             for (auto i = 0; i < n; i++)
2395                 fun[i] = gsl_vector_get(s->f, i) * 1. / sqrt(weight[i]) + ydata[i];
2396 
2397             // calculate weight[i]
2398             for (auto i = 0; i < n; i++) {
2399                 // calculate df[i]
2400                 size_t index = i - 1;
2401                 if (i == 0)
2402                     index = i;
2403                 if (i == n - 1)
2404                     index = i - 2;
2405                 double df = (fun[index + 1] - fun[index]) / (xdata[index + 1] - xdata[index]);
2406                 // printf("df = %g\n", df);
2407 
2408                 double sigmasq = 1.;
2409                 switch (fitData.xWeightsType) { // x-error type: f'(x)^2*s_x^2 = f'(x)/w_x
2410                 case nsl_fit_weight_no:
2411                     break;
2412                 case nsl_fit_weight_direct: // xerror = w_x
2413                     sigmasq = df * df / std::max(xerror[i], minError);
2414                     break;
2415                 case nsl_fit_weight_instrumental: // xerror = s_x
2416                     sigmasq = df * df * xerror[i] * xerror[i];
2417                     break;
2418                 case nsl_fit_weight_inverse: // xerror = 1/w_x = s_x^2
2419                     sigmasq = df * df * xerror[i];
2420                     break;
2421                 case nsl_fit_weight_statistical: // s_x^2 = 1/w_x = x
2422                     sigmasq = xdata[i];
2423                     break;
2424                 case nsl_fit_weight_relative: // s_x^2 = 1/w_x = x^2
2425                     sigmasq = xdata[i] * xdata[i];
2426                     break;
2427                 case nsl_fit_weight_statistical_fit: // unused
2428                 case nsl_fit_weight_relative_fit:
2429                     break;
2430                 }
2431 
2432                 if (yerrorVector.size() > 0) {
2433                     switch (fitData.yWeightsType) { // y-error types: s_y^2 = 1/w_y
2434                     case nsl_fit_weight_no:
2435                         break;
2436                     case nsl_fit_weight_direct: // yerror = w_y
2437                         sigmasq += 1. / std::max(yerror[i], minError);
2438                         break;
2439                     case nsl_fit_weight_instrumental: // yerror = s_y
2440                         sigmasq += yerror[i] * yerror[i];
2441                         break;
2442                     case nsl_fit_weight_inverse: // yerror = 1/w_y
2443                         sigmasq += yerror[i];
2444                         break;
2445                     case nsl_fit_weight_statistical: // unused
2446                     case nsl_fit_weight_relative:
2447                         break;
2448                     case nsl_fit_weight_statistical_fit: // s_y^2 = 1/w_y = Y_i
2449                         sigmasq += fun[i];
2450                         break;
2451                     case nsl_fit_weight_relative_fit: // s_y^2 = 1/w_y = Y_i^2
2452                         sigmasq += fun[i] * fun[i];
2453                         break;
2454                     }
2455                 }
2456 
2457                 // printf ("sigma[%d] = %g\n", i, sqrt(sigmasq));
2458                 weight[i] = 1. / std::max(sigmasq, minError);
2459             }
2460 
2461             // update weights
2462             gsl_multifit_fdfsolver_set(s, &f, &x.vector);
2463 
2464             do { // fit
2465                 iter++;
2466                 writeSolverState(s);
2467                 status = gsl_multifit_fdfsolver_iterate(s);
2468                 // printf ("status = %s\n", gsl_strerror (status));
2469                 if (nf == np) // stop if all parameters fix
2470                     break;
2471 
2472                 if (status) {
2473                     DEBUG("     iter " << iter << ", status = " << gsl_strerror(status));
2474                     if (status == GSL_ETOLX) // change in the position vector falls below machine precision: no progress
2475                         status = GSL_SUCCESS;
2476                     break;
2477                 }
2478                 status = gsl_multifit_test_delta(s->dx, s->x, delta, delta);
2479             } while (status == GSL_CONTINUE && iter < maxIters);
2480 
2481             chi = gsl_blas_dnrm2(s->f);
2482         } while (iter2 < maxIters && fabs(chi - chiOld) > fitData.eps);
2483 
2484         delete[] fun;
2485     }
2486 
2487     delete[] weight;
2488 
2489     // unscale start parameter
2490     for (auto i = 0; i < np; i++)
2491         x_init[i] = nsl_fit_map_bound(x_init[i], x_min[i], x_max[i]);
2492 
2493     // get the covariance matrix
2494     // TODO: scale the Jacobian when limits are used before constructing the covar matrix?
2495     auto* covar = gsl_matrix_alloc(np, np);
2496 #if GSL_MAJOR_VERSION >= 2
2497     // the Jacobian is not part of the solver anymore
2498     auto* J = gsl_matrix_alloc(s->fdf->n, s->fdf->p);
2499     gsl_multifit_fdfsolver_jac(s, J);
2500     gsl_multifit_covar(J, 0.0, covar);
2501 #else
2502     gsl_multifit_covar(s->J, 0.0, covar);
2503 #endif
2504 
2505     // write the result
2506     fitResult.available = true;
2507     fitResult.valid = true;
2508     fitResult.status = gslErrorToString(status);
2509     fitResult.iterations = iter;
2510     fitResult.dof = n - (np - nf); // samples - (parameter - fixed parameter)
2511 
2512     // 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])
2513     // gsl_blas_dasum() - computes the absolute sum \sum |r_i| of the elements of the vector with the elements weight[i]*(Yi - y[i])
2514     fitResult.sse = gsl_pow_2(gsl_blas_dnrm2(s->f));
2515     fitResult.mae = gsl_blas_dasum(s->f) / n;
2516 
2517     // SST needed for coefficient of determination, R-squared and F test
2518     fitResult.sst = gsl_stats_tss(ydata, 1, n);
2519     // for a linear model without intercept R-squared is calculated differently
2520     // see
2521     // 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
2522     if (fitData.modelCategory == nsl_fit_model_basic && fitData.modelType == nsl_fit_model_polynomial && fitData.degree == 1 && x_init[0] == 0) {
2523         DEBUG(" Using alternative R^2 for linear model without intercept");
2524         fitResult.sst = gsl_stats_tss_m(ydata, 1, n, 0);
2525     }
2526     if (fitResult.sst < fitResult.sse) {
2527         DEBUG(" Using alternative R^2 since R^2 would be negative (probably custom model without intercept)");
2528         fitResult.sst = gsl_stats_tss_m(ydata, 1, n, 0);
2529     }
2530     fitResult.calculateResult(n, np);
2531 
2532     // parameter values
2533     fitResult.paramValues.resize(np);
2534     fitResult.errorValues.resize(np);
2535     fitResult.tdist_tValues.resize(np);
2536     fitResult.tdist_pValues.resize(np);
2537     fitResult.marginValues.resize(np);
2538     // GSL: cerr = GSL_MAX_DBL(1., sqrt(fitResult.rms)); // increase error for poor fit
2539     // NIST: cerr = sqrt(fitResult.rms); // increase error for poor fit, decrease for good fit
2540     const double cerr = sqrt(fitResult.rms);
2541     // CI = 100 * (1 - alpha)
2542     const double alpha = 1.0 - fitData.confidenceInterval / 100.;
2543     for (auto i = 0; i < np; i++) {
2544         // scale resulting values if they are bounded
2545         fitResult.paramValues[i] = nsl_fit_map_bound(gsl_vector_get(s->x, i), x_min[i], x_max[i]);
2546         // use results as start values if desired
2547         if (fitData.useResults) {
2548             fitData.paramStartValues.data()[i] = fitResult.paramValues.at(i);
2549             DEBUG(" saving parameter " << i << ": " << fitResult.paramValues[i] << ' ' << fitData.paramStartValues.data()[i]);
2550         }
2551         fitResult.errorValues[i] = cerr * sqrt(gsl_matrix_get(covar, i, i));
2552         fitResult.tdist_tValues[i] = nsl_stats_tdist_t(fitResult.paramValues.at(i), fitResult.errorValues.at(i));
2553         fitResult.tdist_pValues[i] = nsl_stats_tdist_p(fitResult.tdist_tValues.at(i), fitResult.dof);
2554         fitResult.marginValues[i] = nsl_stats_tdist_margin(alpha, fitResult.dof, fitResult.errorValues.at(i));
2555         for (auto j = 0; j <= i; j++)
2556             fitResult.correlationMatrix << gsl_matrix_get(covar, i, j) / sqrt(gsl_matrix_get(covar, i, i)) / sqrt(gsl_matrix_get(covar, j, j));
2557     }
2558 
2559     // residuals for selected range
2560     if (!fitData.autoRange) {
2561         size_t j = 0;
2562         for (int i = 0; i < tmpXDataColumn->rowCount(); i++) {
2563             if (xRange.contains(tmpXDataColumn->valueAt(i)))
2564                 residualsVector->data()[i] = -gsl_vector_get(s->f, j++);
2565             else // outside range
2566                 residualsVector->data()[i] = 0;
2567         }
2568     }
2569 
2570     gsl_multifit_fdfsolver_free(s);
2571     gsl_matrix_free(covar);
2572 }
2573 
2574 /* evaluate fit function (preview == true: use start values, default: false) */
2575 bool XYFitCurvePrivate::evaluate(bool preview) {
2576     DEBUG(Q_FUNC_INFO << ", preview = " << preview);
2577 
2578     // prepare source data columns
2579     DEBUG(Q_FUNC_INFO << ", data source: " << ENUM_TO_STRING(XYAnalysisCurve, DataSourceType, dataSourceType))
2580     const AbstractColumn* tmpXDataColumn = nullptr;
2581     switch (dataSourceType) {
2582     case XYAnalysisCurve::DataSourceType::Spreadsheet:
2583         tmpXDataColumn = xDataColumn;
2584         break;
2585     case XYAnalysisCurve::DataSourceType::Curve:
2586         if (dataSourceCurve)
2587             tmpXDataColumn = dataSourceCurve->xColumn();
2588         break;
2589     case XYAnalysisCurve::DataSourceType::Histogram:
2590         if (dataSourceHistogram)
2591             tmpXDataColumn = dataSourceHistogram->bins();
2592     }
2593 
2594     if (!tmpXDataColumn) {
2595         DEBUG(Q_FUNC_INFO << ", ERROR: Preparing source data column failed!");
2596         return true;
2597     }
2598 
2599     // only needed for preview (else we have all columns)
2600     //  should not harm even if not in preview now that residuals are not cleared
2601     if (preview)
2602         prepareResultColumns();
2603 
2604     if (!xVector || !yVector) {
2605         DEBUG(Q_FUNC_INFO << ", xVector or yVector not defined!");
2606         return true;
2607     }
2608 
2609     if (fitData.model.simplified().isEmpty()) {
2610         DEBUG(Q_FUNC_INFO << ", no fit-model specified.");
2611         return true;
2612     }
2613 
2614     auto* parser = ExpressionParser::getInstance();
2615     Range<double> xRange{tmpXDataColumn->minimum(), tmpXDataColumn->maximum()}; // full data range
2616     if (!fitData.autoEvalRange) { // use custom range for evaluation
2617         if (!fitData.evalRange.isZero()) // avoid zero range
2618             xRange = fitData.evalRange;
2619     }
2620     DEBUG(Q_FUNC_INFO << ", eval range = " << xRange.toStdString());
2621 
2622     xVector->resize((int)fitData.evaluatedPoints);
2623     yVector->resize((int)fitData.evaluatedPoints);
2624     DEBUG(Q_FUNC_INFO << ", vector size = " << xVector->size());
2625 
2626     auto paramValues = fitResult.paramValues;
2627     if (preview) // results not available yet
2628         paramValues = fitData.paramStartValues;
2629 
2630     bool rc = parser->evaluateCartesian(fitData.model, xRange, (int)fitData.evaluatedPoints, xVector, yVector, fitData.paramNames, paramValues);
2631     if (!rc) {
2632         DEBUG(Q_FUNC_INFO << ", ERROR: Parsing fit function failed")
2633         xVector->clear();
2634         yVector->clear();
2635         residualsVector->clear();
2636     }
2637 
2638     return true;
2639 }
2640 
2641 /*!
2642  * writes out the current state of the solver \c s
2643  */
2644 void XYFitCurvePrivate::writeSolverState(gsl_multifit_fdfsolver* s, double chi) {
2645     QString state;
2646 
2647     // current parameter values, semicolon separated
2648     double* min = fitData.paramLowerLimits.data();
2649     double* max = fitData.paramUpperLimits.data();
2650     for (int i = 0; i < fitData.paramNames.size(); ++i) {
2651         const double x = gsl_vector_get(s->x, i);
2652         // map parameter if bounded
2653         state += QString::number(nsl_fit_map_bound(x, min[i], max[i])) + QStringLiteral("\t");
2654     }
2655 
2656     // current value of chi
2657     if (std::isnan(chi))
2658         chi = gsl_blas_dnrm2(s->f);
2659     state += QString::number(chi * chi);
2660     state += QStringLiteral(";");
2661     DEBUG(Q_FUNC_INFO << ", chi^2 = " << chi * chi);
2662 
2663     fitResult.solverOutput += state;
2664 }
2665 
2666 // ##############################################################################
2667 // ##################  Serialization/Deserialization  ###########################
2668 // ##############################################################################
2669 //! Save as XML
2670 void XYFitCurve::save(QXmlStreamWriter* writer) const {
2671     Q_D(const XYFitCurve);
2672 
2673     writer->writeStartElement(QStringLiteral("xyFitCurve"));
2674 
2675     // write the base class
2676     XYAnalysisCurve::save(writer);
2677 
2678     // write xy-fit-curve specific information
2679 
2680     // fit data - only save model expression and parameter names for custom model, otherwise they are set in XYFitCurve::initFitData()
2681     writer->writeStartElement(QStringLiteral("fitData"));
2682     WRITE_COLUMN(d->xErrorColumn, xErrorColumn);
2683     WRITE_COLUMN(d->yErrorColumn, yErrorColumn);
2684     WRITE_PATH(d->dataSourceHistogram, dataSourceHistogram);
2685     writer->writeAttribute(QStringLiteral("autoRange"), QString::number(d->fitData.autoRange));
2686     writer->writeAttribute(QStringLiteral("fitRangeMin"), QString::number(d->fitData.fitRange.start(), 'g', 15));
2687     writer->writeAttribute(QStringLiteral("fitRangeMax"), QString::number(d->fitData.fitRange.end(), 'g', 15));
2688     writer->writeAttribute(QStringLiteral("modelCategory"), QString::number(d->fitData.modelCategory));
2689     writer->writeAttribute(QStringLiteral("modelType"), QString::number(d->fitData.modelType));
2690     writer->writeAttribute(QStringLiteral("xWeightsType"), QString::number(d->fitData.xWeightsType));
2691     writer->writeAttribute(QStringLiteral("weightsType"), QString::number(d->fitData.yWeightsType));
2692     writer->writeAttribute(QStringLiteral("degree"), QString::number(d->fitData.degree));
2693     if (d->fitData.modelCategory == nsl_fit_model_custom)
2694         writer->writeAttribute(QStringLiteral("model"), d->fitData.model);
2695     writer->writeAttribute(QStringLiteral("algorithm"), QString::number(d->fitData.algorithm));
2696     writer->writeAttribute(QStringLiteral("maxIterations"), QString::number(d->fitData.maxIterations));
2697     writer->writeAttribute(QStringLiteral("eps"), QString::number(d->fitData.eps, 'g', 15));
2698     writer->writeAttribute(QStringLiteral("evaluatedPoints"), QString::number(d->fitData.evaluatedPoints));
2699     writer->writeAttribute(QStringLiteral("autoEvalRange"), QString::number(d->fitData.autoEvalRange));
2700     writer->writeAttribute(QStringLiteral("evalRangeMin"), QString::number(d->fitData.evalRange.start(), 'g', 15));
2701     writer->writeAttribute(QStringLiteral("evalRangeMax"), QString::number(d->fitData.evalRange.end(), 'g', 15));
2702     writer->writeAttribute(QStringLiteral("useDataErrors"), QString::number(d->fitData.useDataErrors));
2703     writer->writeAttribute(QStringLiteral("useResults"), QString::number(d->fitData.useResults));
2704     writer->writeAttribute(QStringLiteral("previewEnabled"), QString::number(d->fitData.previewEnabled));
2705     writer->writeAttribute(QStringLiteral("confidenceInterval"), QString::number(d->fitData.confidenceInterval));
2706 
2707     if (d->fitData.modelCategory == nsl_fit_model_custom) {
2708         writer->writeStartElement(QStringLiteral("paramNames"));
2709         for (const QString& name : d->fitData.paramNames)
2710             writer->writeTextElement(QStringLiteral("name"), name);
2711         writer->writeEndElement();
2712     }
2713 
2714     writer->writeStartElement(QStringLiteral("paramStartValues"));
2715     for (const double& value : d->fitData.paramStartValues)
2716         writer->writeTextElement(QStringLiteral("startValue"), QString::number(value, 'g', 15));
2717     writer->writeEndElement();
2718 
2719     // use 16 digits to handle -DBL_MAX
2720     writer->writeStartElement(QStringLiteral("paramLowerLimits"));
2721     for (const double& limit : d->fitData.paramLowerLimits)
2722         writer->writeTextElement(QStringLiteral("lowerLimit"), QString::number(limit, 'g', 16));
2723     writer->writeEndElement();
2724 
2725     // use 16 digits to handle DBL_MAX
2726     writer->writeStartElement(QStringLiteral("paramUpperLimits"));
2727     for (const double& limit : d->fitData.paramUpperLimits)
2728         writer->writeTextElement(QStringLiteral("upperLimit"), QString::number(limit, 'g', 16));
2729     writer->writeEndElement();
2730 
2731     writer->writeStartElement(QStringLiteral("paramFixed"));
2732     for (const bool& fixed : d->fitData.paramFixed)
2733         writer->writeTextElement(QStringLiteral("fixed"), QString::number(fixed));
2734     writer->writeEndElement();
2735 
2736     writer->writeEndElement(); //"fitData"
2737 
2738     // fit results (generated columns and goodness of the fit)
2739     writer->writeStartElement(QStringLiteral("fitResult"));
2740     writer->writeAttribute(QStringLiteral("available"), QString::number(d->fitResult.available));
2741     writer->writeAttribute(QStringLiteral("valid"), QString::number(d->fitResult.valid));
2742     writer->writeAttribute(QStringLiteral("status"), d->fitResult.status);
2743     writer->writeAttribute(QStringLiteral("iterations"), QString::number(d->fitResult.iterations));
2744     writer->writeAttribute(QStringLiteral("time"), QString::number(d->fitResult.elapsedTime));
2745     writer->writeAttribute(QStringLiteral("dof"), QString::number(d->fitResult.dof));
2746     writer->writeAttribute(QStringLiteral("sse"), QString::number(d->fitResult.sse, 'g', 15));
2747     writer->writeAttribute(QStringLiteral("sst"), QString::number(d->fitResult.sst, 'g', 15));
2748     writer->writeAttribute(QStringLiteral("rms"), QString::number(d->fitResult.rms, 'g', 15));
2749     writer->writeAttribute(QStringLiteral("rsd"), QString::number(d->fitResult.rsd, 'g', 15));
2750     writer->writeAttribute(QStringLiteral("mse"), QString::number(d->fitResult.mse, 'g', 15));
2751     writer->writeAttribute(QStringLiteral("rmse"), QString::number(d->fitResult.rmse, 'g', 15));
2752     writer->writeAttribute(QStringLiteral("mae"), QString::number(d->fitResult.mae, 'g', 15));
2753     writer->writeAttribute(QStringLiteral("rsquare"), QString::number(d->fitResult.rsquare, 'g', 15));
2754     writer->writeAttribute(QStringLiteral("rsquareAdj"), QString::number(d->fitResult.rsquareAdj, 'g', 15));
2755     writer->writeAttribute(QStringLiteral("chisq_p"), QString::number(d->fitResult.chisq_p, 'g', 15));
2756     writer->writeAttribute(QStringLiteral("fdist_F"), QString::number(d->fitResult.fdist_F, 'g', 15));
2757     writer->writeAttribute(QStringLiteral("fdist_p"), QString::number(d->fitResult.fdist_p, 'g', 15));
2758     writer->writeAttribute(QStringLiteral("aic"), QString::number(d->fitResult.aic, 'g', 15));
2759     writer->writeAttribute(QStringLiteral("bic"), QString::number(d->fitResult.bic, 'g', 15));
2760     writer->writeAttribute(QStringLiteral("solverOutput"), d->fitResult.solverOutput);
2761 
2762     writer->writeStartElement(QStringLiteral("paramValues"));
2763     for (const double& value : d->fitResult.paramValues)
2764         writer->writeTextElement(QStringLiteral("value"), QString::number(value, 'g', 15));
2765     writer->writeEndElement();
2766 
2767     writer->writeStartElement(QStringLiteral("errorValues"));
2768     for (const double& value : d->fitResult.errorValues)
2769         writer->writeTextElement(QStringLiteral("error"), QString::number(value, 'g', 15));
2770     writer->writeEndElement();
2771 
2772     writer->writeStartElement(QStringLiteral("tdist_tValues"));
2773     for (const double& value : d->fitResult.tdist_tValues)
2774         writer->writeTextElement(QStringLiteral("tdist_t"), QString::number(value, 'g', 15));
2775     writer->writeEndElement();
2776 
2777     writer->writeStartElement(QStringLiteral("tdist_pValues"));
2778     for (const double& value : d->fitResult.tdist_pValues)
2779         writer->writeTextElement(QStringLiteral("tdist_p"), QString::number(value, 'g', 15));
2780     writer->writeEndElement();
2781 
2782     writer->writeStartElement(QStringLiteral("tdist_marginValues"));
2783     for (const double& value : d->fitResult.marginValues)
2784         writer->writeTextElement(QStringLiteral("tdist_margin"), QString::number(value, 'g', 15));
2785     writer->writeEndElement();
2786     writer->writeStartElement(QStringLiteral("tdist_margin2Values"));
2787     for (const double& value : d->fitResult.margin2Values)
2788         writer->writeTextElement(QStringLiteral("tdist_margin2"), QString::number(value, 'g', 15));
2789     writer->writeEndElement();
2790 
2791     writer->writeStartElement(QStringLiteral("correlationMatrix"));
2792     for (const double& value : d->fitResult.correlationMatrix)
2793         writer->writeTextElement(QStringLiteral("correlation"), QString::number(value, 'g', 15));
2794     writer->writeEndElement();
2795 
2796     // save calculated columns if available
2797     if (saveCalculations() && d->xColumn && d->yColumn && d->residualsColumn) {
2798         d->xColumn->save(writer);
2799         d->yColumn->save(writer);
2800         d->residualsColumn->save(writer);
2801     }
2802 
2803     writer->writeEndElement(); //"fitResult"
2804     writer->writeEndElement(); //"xyFitCurve"
2805 }
2806 
2807 //! Load from XML
2808 bool XYFitCurve::load(XmlStreamReader* reader, bool preview) {
2809     Q_D(XYFitCurve);
2810 
2811     QXmlStreamAttributes attribs;
2812     QString str, model;
2813 
2814     while (!reader->atEnd()) {
2815         reader->readNext();
2816         if (reader->isEndElement() && reader->name() == QLatin1String("xyFitCurve"))
2817             break;
2818 
2819         if (!reader->isStartElement())
2820             continue;
2821 
2822         if (reader->name() == QLatin1String("xyAnalysisCurve")) {
2823             if (!XYAnalysisCurve::load(reader, preview))
2824                 return false;
2825         } else if (!preview && reader->name() == QLatin1String("fitData")) {
2826             attribs = reader->attributes();
2827 
2828             READ_COLUMN(xErrorColumn);
2829             READ_COLUMN(yErrorColumn);
2830             READ_PATH(dataSourceHistogram);
2831 
2832             READ_INT_VALUE("autoRange", fitData.autoRange, bool);
2833             READ_DOUBLE_VALUE("xRangeMin", fitData.fitRange.start()); // old name
2834             READ_DOUBLE_VALUE("xRangeMax", fitData.fitRange.end()); // old name
2835             READ_DOUBLE_VALUE("fitRangeMin", fitData.fitRange.start());
2836             READ_DOUBLE_VALUE("fitRangeMax", fitData.fitRange.end());
2837             READ_INT_VALUE("modelCategory", fitData.modelCategory, nsl_fit_model_category);
2838             READ_INT_VALUE("modelType", fitData.modelType, int);
2839             READ_INT_VALUE("xWeightsType", fitData.xWeightsType, nsl_fit_weight_type);
2840             READ_INT_VALUE("weightsType", fitData.yWeightsType, nsl_fit_weight_type);
2841             READ_INT_VALUE("degree", fitData.degree, int);
2842             // older projects have custom models with category == type == 0! So read the model
2843             READ_STRING_VALUE("model", fitData.model);
2844             model = d->fitData.model;
2845             DEBUG("got model = " << STDSTRING(model));
2846             READ_INT_VALUE("algorithm", fitData.algorithm, nsl_fit_algorithm);
2847 
2848             READ_INT_VALUE("maxIterations", fitData.maxIterations, int);
2849             READ_DOUBLE_VALUE("eps", fitData.eps);
2850             READ_INT_VALUE("fittedPoints", fitData.evaluatedPoints, size_t); // old name
2851             READ_INT_VALUE("evaluatedPoints", fitData.evaluatedPoints, size_t);
2852             READ_INT_VALUE("evaluateFullRange", fitData.autoEvalRange, bool); // old name
2853             READ_INT_VALUE("autoEvalRange", fitData.autoEvalRange, bool);
2854             READ_DOUBLE_VALUE("evalRangeMin", fitData.evalRange.start());
2855             READ_DOUBLE_VALUE("evalRangeMax", fitData.evalRange.end());
2856             READ_INT_VALUE("useDataErrors", fitData.useDataErrors, bool);
2857             READ_INT_VALUE("useResults", fitData.useResults, bool);
2858             READ_INT_VALUE("previewEnabled", fitData.previewEnabled, bool);
2859             READ_DOUBLE_VALUE("confidenceInterval", fitData.confidenceInterval);
2860 
2861             // set the model expression and the parameter names (can be derived from the saved values for category, type and degree)
2862             XYFitCurve::initFitData(d->fitData);
2863             // remove default names and start values (will be read from project later)
2864             d->fitData.paramStartValues.clear();
2865 
2866         } else if (!preview && reader->name() == QLatin1String("name")) { // needed for custom model
2867             d->fitData.paramNames << reader->readElementText();
2868         } else if (!preview && reader->name() == QLatin1String("startValue")) {
2869             d->fitData.paramStartValues << reader->readElementText().toDouble();
2870         } else if (!preview && reader->name() == QLatin1String("fixed")) {
2871             d->fitData.paramFixed << (bool)reader->readElementText().toInt();
2872         } else if (!preview && reader->name() == QLatin1String("lowerLimit")) {
2873             bool ok;
2874             double x = reader->readElementText().toDouble(&ok);
2875             if (ok) // -DBL_MAX results in conversion error
2876                 d->fitData.paramLowerLimits << x;
2877             else
2878                 d->fitData.paramLowerLimits << -std::numeric_limits<double>::max();
2879         } else if (!preview && reader->name() == QLatin1String("upperLimit")) {
2880             bool ok;
2881             double x = reader->readElementText().toDouble(&ok);
2882             if (ok) // DBL_MAX results in conversion error
2883                 d->fitData.paramUpperLimits << x;
2884             else
2885                 d->fitData.paramUpperLimits << std::numeric_limits<double>::max();
2886         } else if (!preview && reader->name() == QLatin1String("value")) {
2887             d->fitResult.paramValues << reader->readElementText().toDouble();
2888         } else if (!preview && reader->name() == QLatin1String("error")) {
2889             d->fitResult.errorValues << reader->readElementText().toDouble();
2890         } else if (!preview && reader->name() == QLatin1String("tdist_t")) {
2891             d->fitResult.tdist_tValues << reader->readElementText().toDouble();
2892         } else if (!preview && reader->name() == QLatin1String("tdist_p")) {
2893             d->fitResult.tdist_pValues << reader->readElementText().toDouble();
2894         } else if (!preview && reader->name() == QLatin1String("tdist_margin")) {
2895             d->fitResult.marginValues << reader->readElementText().toDouble();
2896         } else if (!preview && reader->name() == QLatin1String("tdist_margin2")) {
2897             d->fitResult.margin2Values << reader->readElementText().toDouble();
2898         } else if (!preview && reader->name() == QLatin1String("correlation")) {
2899             d->fitResult.correlationMatrix << reader->readElementText().toDouble();
2900         } else if (!preview && reader->name() == QLatin1String("fitResult")) {
2901             attribs = reader->attributes();
2902 
2903             READ_INT_VALUE("available", fitResult.available, int);
2904             READ_INT_VALUE("valid", fitResult.valid, int);
2905             READ_STRING_VALUE("status", fitResult.status);
2906             READ_INT_VALUE("iterations", fitResult.iterations, int);
2907             READ_INT_VALUE("time", fitResult.elapsedTime, int);
2908             READ_DOUBLE_VALUE("dof", fitResult.dof);
2909             READ_DOUBLE_VALUE("sse", fitResult.sse);
2910             READ_DOUBLE_VALUE("sst", fitResult.sst);
2911             READ_DOUBLE_VALUE("rms", fitResult.rms);
2912             READ_DOUBLE_VALUE("rsd", fitResult.rsd);
2913             READ_DOUBLE_VALUE("mse", fitResult.mse);
2914             READ_DOUBLE_VALUE("rmse", fitResult.rmse);
2915             READ_DOUBLE_VALUE("mae", fitResult.mae);
2916             READ_DOUBLE_VALUE("rsquare", fitResult.rsquare);
2917             READ_DOUBLE_VALUE("rsquareAdj", fitResult.rsquareAdj);
2918             READ_DOUBLE_VALUE("chisq_p", fitResult.chisq_p);
2919             READ_DOUBLE_VALUE("fdist_F", fitResult.fdist_F);
2920             READ_DOUBLE_VALUE("fdist_p", fitResult.fdist_p);
2921             READ_DOUBLE_VALUE("aic", fitResult.aic);
2922             READ_DOUBLE_VALUE("bic", fitResult.bic);
2923             READ_STRING_VALUE("solverOutput", fitResult.solverOutput);
2924         } else if (reader->name() == QLatin1String("column")) {
2925             Column* column = new Column(QString(), AbstractColumn::ColumnMode::Double);
2926             if (!column->load(reader, preview)) {
2927                 delete column;
2928                 return false;
2929             }
2930             DEBUG("############################   reading column " << STDSTRING(column->name()))
2931             if (column->name() == QLatin1String("x"))
2932                 d->xColumn = column;
2933             else if (column->name() == QLatin1String("y"))
2934                 d->yColumn = column;
2935             else if (column->name() == QLatin1String("residuals"))
2936                 d->residualsColumn = column;
2937         } else { // unknown element
2938             reader->raiseUnknownElementWarning();
2939             if (!reader->skipToEndElement())
2940                 return false;
2941         }
2942     }
2943 
2944     ////////////////////////////// fix old projects /////////////////////////
2945 
2946     // reset model type of old projects due to new model style
2947     if (d->fitData.modelCategory == nsl_fit_model_basic && d->fitData.modelType >= NSL_FIT_MODEL_BASIC_COUNT) {
2948         DEBUG(Q_FUNC_INFO << ", RESET old fit model");
2949         d->fitData.modelType = 0;
2950         d->fitData.degree = 1;
2951         d->fitData.paramNames.clear();
2952         d->fitData.paramNamesUtf8.clear();
2953         // reset size of fields not touched by initFitData()
2954         d->fitData.paramStartValues.resize(2);
2955         d->fitData.paramFixed.resize(2);
2956         d->fitResult.paramValues.resize(2);
2957         d->fitResult.errorValues.resize(2);
2958         d->fitResult.tdist_tValues.resize(2);
2959         d->fitResult.tdist_pValues.resize(2);
2960         d->fitResult.marginValues.resize(2);
2961         d->fitResult.margin2Values.resize(2);
2962     }
2963 
2964     // older projects also save the param names for non-custom models: remove them
2965     while (d->fitData.paramNames.size() > d->fitData.paramStartValues.size())
2966         d->fitData.paramNames.removeLast();
2967 
2968     // Utf8 names missing in old projects
2969     if (d->fitData.paramNamesUtf8.isEmpty())
2970         d->fitData.paramNamesUtf8 << d->fitData.paramNames;
2971 
2972     // if we have more paramNames than the saved model type, we have a custom model
2973     if (d->fitData.paramNamesUtf8.size() < d->fitData.paramNames.size()) {
2974         d->fitData.modelCategory = nsl_fit_model_custom;
2975         d->fitData.model = model;
2976         d->fitData.paramNamesUtf8 = d->fitData.paramNames;
2977     }
2978 
2979     // not present in old projects
2980     int np = d->fitResult.paramValues.size();
2981     if (d->fitResult.tdist_tValues.size() == 0)
2982         d->fitResult.tdist_tValues.resize(np);
2983     if (d->fitResult.tdist_pValues.size() == 0)
2984         d->fitResult.tdist_pValues.resize(np);
2985     if (d->fitResult.marginValues.size() == 0)
2986         d->fitResult.marginValues.resize(np);
2987     if (d->fitResult.margin2Values.size() == 0)
2988         d->fitResult.margin2Values.resize(np);
2989     if (d->fitResult.correlationMatrix.size() == 0)
2990         d->fitResult.correlationMatrix.resize(np * (np + 1) / 2);
2991 
2992     // Loading done. Check some parameter
2993     DEBUG(Q_FUNC_INFO << ", model category = " << d->fitData.modelCategory);
2994     DEBUG(Q_FUNC_INFO << ", model type = " << d->fitData.modelType);
2995     DEBUG(Q_FUNC_INFO << ", # params = " << d->fitData.paramNames.size());
2996     DEBUG(Q_FUNC_INFO << ", # start values = " << d->fitData.paramStartValues.size());
2997     // for (const auto& value : d->fitData.paramStartValues)
2998     //  DEBUG("XYFitCurve::load() # start value = " << value);
2999 
3000     if (preview)
3001         return true;
3002 
3003     // wait for data to be read before using the pointers
3004     QThreadPool::globalInstance()->waitForDone();
3005 
3006     if (d->xColumn && d->yColumn && d->residualsColumn) {
3007         d->xColumn->setHidden(true);
3008         addChild(d->xColumn);
3009 
3010         d->yColumn->setHidden(true);
3011         addChild(d->yColumn);
3012 
3013         addChild(d->residualsColumn);
3014 
3015         d->xVector = static_cast<QVector<double>*>(d->xColumn->data());
3016         d->yVector = static_cast<QVector<double>*>(d->yColumn->data());
3017         d->residualsVector = static_cast<QVector<double>*>(d->residualsColumn->data());
3018 
3019         static_cast<XYCurvePrivate*>(d_ptr)->xColumn = d->xColumn;
3020         static_cast<XYCurvePrivate*>(d_ptr)->yColumn = d->yColumn;
3021 
3022         recalcLogicalPoints();
3023     }
3024 
3025     return true;
3026 }