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 = ¶ms; 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 }