diff --git a/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp b/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp index e04e694f6..2ea019a89 100644 --- a/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp +++ b/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp @@ -1,2392 +1,2388 @@  /*************************************************************************** File : XYFitCurve.cpp Project : LabPlot Description : A xy-curve defined by a fit model -------------------------------------------------------------------- Copyright : (C) 2014-2017 Alexander Semke (alexander.semke@web.de) Copyright : (C) 2016-2018 Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ /*! \class XYFitCurve \brief A xy-curve defined by a fit model \ingroup worksheet */ #include "XYFitCurve.h" #include "XYFitCurvePrivate.h" #include "backend/core/AbstractColumn.h" #include "backend/core/column/Column.h" #include "backend/lib/commandtemplates.h" #include "backend/lib/macros.h" #include "backend/gsl/errors.h" #include "backend/gsl/ExpressionParser.h" extern "C" { #include #include #include #include #include #include #include "backend/gsl/parser.h" #include "backend/nsl/nsl_sf_stats.h" #include "backend/nsl/nsl_stats.h" } #include #include #include XYFitCurve::XYFitCurve(const QString& name) : XYAnalysisCurve(name, new XYFitCurvePrivate(this)) { } XYFitCurve::XYFitCurve(const QString& name, XYFitCurvePrivate* dd) : XYAnalysisCurve(name, dd) { } //no need to delete the d-pointer here - it inherits from QGraphicsItem //and is deleted during the cleanup in QGraphicsScene XYFitCurve::~XYFitCurve() = default; void XYFitCurve::recalculate() { Q_D(XYFitCurve); d->recalculate(); } void XYFitCurve::evaluate(bool preview) { Q_D(XYFitCurve); d->evaluate(preview); } void XYFitCurve::initStartValues(const XYCurve* curve) { Q_D(XYFitCurve); XYFitCurve::FitData& fitData = d->fitData; initStartValues(fitData, curve); } void XYFitCurve::initStartValues(XYFitCurve::FitData& fitData, const XYCurve* curve) { DEBUG("XYFitCurve::initStartValues()"); if (!curve) { DEBUG(" no curve given"); return; } const Column* tmpXDataColumn = dynamic_cast(curve->xColumn()); const Column* tmpYDataColumn = dynamic_cast(curve->yColumn()); if (!tmpXDataColumn || !tmpYDataColumn) { DEBUG(" data columns not available"); return; } DEBUG(" x data rows = " << tmpXDataColumn->rowCount()); nsl_fit_model_category modelCategory = fitData.modelCategory; int modelType = fitData.modelType; int degree = fitData.degree; DEBUG(" fit model type = " << modelType << ", degree = " << degree); QVector& paramStartValues = fitData.paramStartValues; //QVector* xVector = static_cast* >(tmpXDataColumn->data()); //double xmean = gsl_stats_mean(xVector->constData(), 1, tmpXDataColumn->rowCount()); double xmin = tmpXDataColumn->minimum(); double xmax = tmpXDataColumn->maximum(); //double ymin = tmpYDataColumn->minimum(); //double ymax = tmpYDataColumn->maximum(); double xrange = xmax-xmin; //double yrange = ymax-ymin; DEBUG(" x min/max = " << xmin << ' ' << xmax); //DEBUG(" y min/max = " << ymin << ' ' << ymax); switch (modelCategory) { case nsl_fit_model_basic: switch (modelType) { case nsl_fit_model_polynomial: // not needed (works anyway) break; //TODO: handle basic models case nsl_fit_model_power: case nsl_fit_model_exponential: case nsl_fit_model_inverse_exponential: case nsl_fit_model_fourier: break; } break; case nsl_fit_model_peak: // use equidistant mu's and (xmax-xmin)/(10*degree) as sigma(, gamma) switch (modelType) { case nsl_fit_model_gaussian: case nsl_fit_model_lorentz: case nsl_fit_model_sech: case nsl_fit_model_logistic: for (int d = 0; d < degree; d++) { paramStartValues[3*d+2] = xmin + (d+1.)*xrange/(degree+1.); // mu paramStartValues[3*d+1] = xrange/(10.*degree); // sigma } break; case nsl_fit_model_voigt: for (int d = 0; d < degree; d++) { paramStartValues[4*d+1] = xmin + (d+1.)*xrange/(degree+1.); // mu paramStartValues[4*d+2] = xrange/(10.*degree); // sigma paramStartValues[4*d+3] = xrange/(10.*degree); // gamma } break; case nsl_fit_model_pseudovoigt1: for (int d = 0; d < degree; d++) { paramStartValues[4*d+1] = 0.5; // eta paramStartValues[4*d+2] = xrange/(10.*degree); // sigma paramStartValues[4*d+3] = xmin + (d+1.)*xrange/(degree+1.); // mu } break; } break; case nsl_fit_model_growth: switch (modelType) { case nsl_fit_model_atan: case nsl_fit_model_tanh: case nsl_fit_model_algebraic_sigmoid: case nsl_fit_model_erf: case nsl_fit_model_gudermann: case nsl_fit_model_sigmoid: // use (xmax+xmin)/2 as mu and (xmax-xmin)/10 as sigma paramStartValues[1] = (xmax+xmin)/2.; paramStartValues[2] = xrange/10.; break; case nsl_fit_model_hill: paramStartValues[2] = xrange/10.; break; case nsl_fit_model_gompertz: //TODO break; } break; case nsl_fit_model_distribution: switch (modelType) { case nsl_sf_stats_gaussian: case nsl_sf_stats_laplace: case nsl_sf_stats_rayleigh_tail: case nsl_sf_stats_lognormal: case nsl_sf_stats_logistic: case nsl_sf_stats_sech: case nsl_sf_stats_cauchy_lorentz: case nsl_sf_stats_levy: // use (xmax+xmin)/2 as mu and (xmax-xmin)/10 as sigma paramStartValues[2] = (xmin+xmax)/2.; paramStartValues[1] = xrange/10.; break; //TODO: other types default: break; } break; case nsl_fit_model_custom: // not possible break; } } /*! * sets the parameter names for given model category, model type and degree in \c fitData for given action */ void XYFitCurve::initFitData(PlotDataDialog::AnalysisAction action) { if (!action) return; Q_D(XYFitCurve); XYFitCurve::FitData& fitData = d->fitData; if (action == PlotDataDialog::FitLinear) { //Linear fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = (int)nsl_fit_model_polynomial; fitData.degree = 1; } else if (action == PlotDataDialog::FitPower) { //Power fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = (int)nsl_fit_model_power; fitData.degree = 1; } else if (action == PlotDataDialog::FitExp1) { //Exponential (degree 1) fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = (int)nsl_fit_model_exponential; fitData.degree = 1; } else if (action == PlotDataDialog::FitExp2) { //Exponential (degree 2) fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = (int)nsl_fit_model_exponential; fitData.degree = 2; } else if (action == PlotDataDialog::FitInvExp) { //Inverse exponential fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = (int)nsl_fit_model_inverse_exponential; } else if (action == PlotDataDialog::FitGauss) { //Gauss fitData.modelCategory = nsl_fit_model_peak; fitData.modelType = (int)nsl_fit_model_gaussian; fitData.degree = 1; } else if (action == PlotDataDialog::FitCauchyLorentz) { //Cauchy-Lorentz fitData.modelCategory = nsl_fit_model_peak; fitData.modelType = (int)nsl_fit_model_lorentz; fitData.degree = 1; } else if (action == PlotDataDialog::FitTan) { //Arc tangent fitData.modelCategory = nsl_fit_model_growth; fitData.modelType = (int)nsl_fit_model_atan; } else if (action == PlotDataDialog::FitTanh) { //Hyperbolic tangent fitData.modelCategory = nsl_fit_model_growth; fitData.modelType = (int)nsl_fit_model_tanh; } else if (action == PlotDataDialog::FitErrFunc) { //Error function fitData.modelCategory = nsl_fit_model_growth; fitData.modelType = (int)nsl_fit_model_erf; } else { //Custom fitData.modelCategory = nsl_fit_model_custom; fitData.modelType = 0; } XYFitCurve::initFitData(fitData); } /*! * sets the model expression and the parameter names for given model category, model type and degree in \c fitData */ void XYFitCurve::initFitData(XYFitCurve::FitData& fitData) { nsl_fit_model_category modelCategory = fitData.modelCategory; int modelType = fitData.modelType; QString& model = fitData.model; QStringList& paramNames = fitData.paramNames; QStringList& paramNamesUtf8 = fitData.paramNamesUtf8; int degree = fitData.degree; QVector& paramStartValues = fitData.paramStartValues; QVector& paramLowerLimits = fitData.paramLowerLimits; QVector& paramUpperLimits = fitData.paramUpperLimits; QVector& paramFixed = fitData.paramFixed; if (modelCategory != nsl_fit_model_custom) { DEBUG("XYFitCurve::initFitData() for model category = " << nsl_fit_model_category_name[modelCategory] << ", model type = " << modelType << ", degree = " << degree); paramNames.clear(); } else { DEBUG("XYFitCurve::initFitData() for model category = nsl_fit_model_custom, model type = " << modelType << ", degree = " << degree); } paramNamesUtf8.clear(); // 10 indices used in multi degree models QStringList indices = { UTF8_QSTRING("₁"), UTF8_QSTRING("₂"), UTF8_QSTRING("₃"), UTF8_QSTRING("₄"), UTF8_QSTRING("₅"), UTF8_QSTRING("₆"), UTF8_QSTRING("₇"), UTF8_QSTRING("₈"), UTF8_QSTRING("₉"), UTF8_QSTRING("₁₀")}; switch (modelCategory) { case nsl_fit_model_basic: model = nsl_fit_model_basic_equation[fitData.modelType]; switch (modelType) { case nsl_fit_model_polynomial: paramNames << "c0" << "c1"; paramNamesUtf8 << UTF8_QSTRING("c₀") << UTF8_QSTRING("c₁"); if (degree == 2) { model += " + c2*x^2"; paramNames << "c2"; paramNamesUtf8 << UTF8_QSTRING("c₂"); } else if (degree > 2) { for (int i = 2; i <= degree; ++i) { QString numStr = QString::number(i); model += "+c" + numStr + "*x^" + numStr; paramNames << 'c' + numStr; paramNamesUtf8 << 'c' + indices[i-1]; } } break; case nsl_fit_model_power: if (degree == 1) { paramNames << "a" << "b"; } else { paramNames << "a" << "b" << "c"; model = "a + b*x^c"; } break; case nsl_fit_model_exponential: if (degree == 1) { paramNames << "a" << "b"; } else { for (int i = 1; i <= degree; i++) { QString numStr = QString::number(i); if (i == 1) model = "a1*exp(b1*x)"; else model += " + a" + numStr + "*exp(b" + numStr + "*x)"; paramNames << 'a' + numStr << 'b' + numStr; paramNamesUtf8 << 'a' + indices[i-1] << 'b' + indices[i-1]; } } break; case nsl_fit_model_inverse_exponential: paramNames << "a" << "b" << "c"; break; case nsl_fit_model_fourier: paramNames << "w" << "a0" << "a1" << "b1"; paramNamesUtf8 << UTF8_QSTRING("ω") << UTF8_QSTRING("a₀") << UTF8_QSTRING("a₁") << UTF8_QSTRING("b₁"); if (degree > 1) { for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); model += "+ (a" + numStr + "*cos(" + numStr + "*w*x) + b" + numStr + "*sin(" + numStr + "*w*x))"; paramNames << 'a' + numStr << 'b' + numStr; paramNamesUtf8 << 'a' + indices[i-1] << 'b' + indices[i-1]; } } break; } break; case nsl_fit_model_peak: model = nsl_fit_model_peak_equation[fitData.modelType]; switch (modelType) { case nsl_fit_model_gaussian: switch (degree) { case 1: paramNames << "a" << "s" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ"); break; default: model = "1./sqrt(2*pi) * ("; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += 'a' + numStr + "/s" + numStr + "* exp(-((x-mu" + numStr + ")/s" + numStr + ")^2/2)"; paramNames << 'a' + numStr << 's' + numStr << "mu" + numStr; paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1]; } model += ')'; } break; case nsl_fit_model_lorentz: switch (degree) { case 1: paramNames << "a" << "g" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("γ") << UTF8_QSTRING("μ"); break; default: model = "1./pi * ("; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += 'a' + numStr + " * g" + numStr + "/(g" + numStr + "^2+(x-mu" + numStr + ")^2)"; paramNames << 'a' + numStr << 'g' + numStr << "mu" + numStr; paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("γ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1]; } model += ')'; } break; case nsl_fit_model_sech: switch (degree) { case 1: paramNames << "a" << "s" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ"); break; default: model = "1/pi * ("; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += 'a' + numStr + "/s" + numStr + "* sech((x-mu" + numStr + ")/s" + numStr + ')'; paramNames << 'a' + numStr << 's' + numStr << "mu" + numStr; paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1]; } model += ')'; } break; case nsl_fit_model_logistic: switch (degree) { case 1: paramNames << "a" << "s" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ"); break; default: model = "1/4 * ("; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += 'a' + numStr + "/s" + numStr + "* sech((x-mu" + numStr + ")/2/s" + numStr + ")**2"; paramNames << 'a' + numStr << 's' + numStr << "mu" + numStr; paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1]; } model += ')'; } break; case nsl_fit_model_voigt: switch(degree) { case 1: paramNames << "a" << "mu" << "s" << "g"; paramNamesUtf8 << "A" << UTF8_QSTRING("μ") << UTF8_QSTRING("σ") << UTF8_QSTRING("γ"); break; default: model = ""; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += 'a' + numStr + "*voigt(x-mu" + numStr + ",s" + numStr + ",g" + numStr + ')'; paramNames << 'a' + numStr << "mu" + numStr << 's' + numStr << 'g' + numStr; paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1] << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("γ") + indices[i-1]; } } break; case nsl_fit_model_pseudovoigt1: switch(degree) { case 1: paramNames << "a" << "et" << "w" << "mu"; // eta function exists! paramNamesUtf8 << "A" << UTF8_QSTRING("η") << "w" << UTF8_QSTRING("μ"); break; default: model=""; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += 'a' + numStr + "*pseudovoigt1(x-mu" + numStr + ",eta" + numStr + ",w" + numStr + ')'; paramNames << 'a' + numStr << "eta" + numStr << 'w' + numStr << "mu" + numStr; paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("η") + indices[i-1] << 'w' + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1]; } } break; } break; case nsl_fit_model_growth: model = nsl_fit_model_growth_equation[fitData.modelType]; switch (modelType) { case nsl_fit_model_atan: case nsl_fit_model_tanh: case nsl_fit_model_algebraic_sigmoid: case nsl_fit_model_erf: case nsl_fit_model_gudermann: paramNames << "a" << "mu" << "s"; paramNamesUtf8 << "A" << UTF8_QSTRING("μ") << UTF8_QSTRING("σ"); break; case nsl_fit_model_sigmoid: paramNames << "a" << "mu" << "k"; paramNamesUtf8 << "A" << UTF8_QSTRING("μ") << "k"; break; case nsl_fit_model_hill: paramNames << "a" << "n" << "a"; paramNamesUtf8 << "A" << "n" << UTF8_QSTRING("σ"); break; case nsl_fit_model_gompertz: paramNames << "a" << "b" << "c"; break; } break; case nsl_fit_model_distribution: model = nsl_sf_stats_distribution_equation[fitData.modelType]; switch (modelType) { case nsl_sf_stats_gaussian: case nsl_sf_stats_laplace: case nsl_sf_stats_rayleigh_tail: case nsl_sf_stats_lognormal: case nsl_sf_stats_logistic: case nsl_sf_stats_sech: paramNames << "a" << "s" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ"); break; case nsl_sf_stats_gaussian_tail: paramNames << "A" << "s" << "a" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << "a" << UTF8_QSTRING("μ"); break; case nsl_sf_stats_exponential: paramNames << "a" << "l" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("λ") << UTF8_QSTRING("μ"); break; case nsl_sf_stats_exponential_power: paramNames << "a" << "s" << "b" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << "b" << UTF8_QSTRING("μ"); break; case nsl_sf_stats_cauchy_lorentz: case nsl_sf_stats_levy: paramNames << "a" << "g" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("γ") << UTF8_QSTRING("μ"); break; case nsl_sf_stats_rayleigh: paramNames << "a" << "s"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ"); break; case nsl_sf_stats_landau: paramNames << "a"; paramNamesUtf8 << "A"; break; case nsl_sf_stats_levy_alpha_stable: // unused distributions case nsl_sf_stats_levy_skew_alpha_stable: case nsl_sf_stats_bernoulli: break; case nsl_sf_stats_gamma: paramNames << "a" << "k" << "t"; paramNamesUtf8 << "A"<< "k" << UTF8_QSTRING("θ"); break; case nsl_sf_stats_flat: paramNames << "A" << "b" << "a"; break; case nsl_sf_stats_chi_squared: paramNames << "a" << "n"; paramNamesUtf8 << "A" << "n"; break; case nsl_sf_stats_fdist: paramNames << "a" << "n1" << "n2"; paramNamesUtf8 << "A" << UTF8_QSTRING("ν₁") << UTF8_QSTRING("ν₂"); break; case nsl_sf_stats_tdist: paramNames << "a" << "n"; paramNamesUtf8 << "A" << UTF8_QSTRING("ν"); break; case nsl_sf_stats_beta: case nsl_sf_stats_pareto: paramNames << "A" << "a" << "b"; break; case nsl_sf_stats_weibull: paramNames << "a" << "k" << "l" << "mu"; paramNamesUtf8 << "A" << "k" << UTF8_QSTRING("λ") << UTF8_QSTRING("μ"); break; case nsl_sf_stats_gumbel1: paramNames << "a" << "s" << "mu" << "b"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ") << UTF8_QSTRING("β"); break; case nsl_sf_stats_gumbel2: paramNames << "A" << "a" << "b" << "mu"; paramNamesUtf8 << "A" << "a" << "b" << UTF8_QSTRING("μ"); break; case nsl_sf_stats_poisson: paramNames << "a" << "l"; paramNamesUtf8 << "A" << UTF8_QSTRING("λ"); break; case nsl_sf_stats_binomial: case nsl_sf_stats_negative_binomial: case nsl_sf_stats_pascal: paramNames << "a" << "p" << "n"; paramNamesUtf8 << "A" << "p" << "n"; break; case nsl_sf_stats_geometric: case nsl_sf_stats_logarithmic: paramNames << "a" << "p"; paramNamesUtf8 << "A" << "p"; break; case nsl_sf_stats_hypergeometric: paramNames << "a" << "n1" << "n2" << "t"; paramNamesUtf8 << "A" << UTF8_QSTRING("n₁") << UTF8_QSTRING("n₂") << "t"; break; case nsl_sf_stats_maxwell_boltzmann: paramNames << "a" << "s"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ"); break; case nsl_sf_stats_frechet: paramNames << "a" << "g" << "s" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("γ") << UTF8_QSTRING("σ") << UTF8_QSTRING("μ"); break; } break; case nsl_fit_model_custom: break; } DEBUG("model: " << model.toStdString()); if (paramNamesUtf8.isEmpty()) paramNamesUtf8 << paramNames; //resize the vector for the start values and set the elements to 1.0 //in case a custom model is used, do nothing, we take over the previous values if (modelCategory != nsl_fit_model_custom) { const int np = paramNames.size(); paramStartValues.resize(np); paramFixed.resize(np); paramLowerLimits.resize(np); paramUpperLimits.resize(np); for (int i = 0; i < np; ++i) { paramStartValues[i] = 1.0; paramFixed[i] = false; paramLowerLimits[i] = -std::numeric_limits::max(); paramUpperLimits[i] = std::numeric_limits::max(); } // set some model-dependent start values // TODO: see initStartValues() if (modelCategory == nsl_fit_model_distribution) { if (modelType == (int)nsl_sf_stats_flat) paramStartValues[2] = -1.0; else if (modelType == (int)nsl_sf_stats_levy) paramStartValues[2] = 0.0; else if (modelType == (int)nsl_sf_stats_exponential_power || modelType == (int)nsl_sf_stats_weibull || modelType == (int)nsl_sf_stats_gumbel2 || modelType == (int)nsl_sf_stats_frechet) paramStartValues[3] = 0.0; else if (modelType == (int)nsl_sf_stats_binomial || modelType == (int)nsl_sf_stats_negative_binomial || modelType == (int)nsl_sf_stats_pascal || modelType == (int)nsl_sf_stats_geometric || modelType == (int)nsl_sf_stats_logarithmic) paramStartValues[1] = 0.5; } } } /*! Returns an icon to be used in the project explorer. */ QIcon XYFitCurve::icon() const { return QIcon::fromTheme("labplot-xy-fit-curve"); } //############################################################################## //########################## getter methods ################################## //############################################################################## BASIC_SHARED_D_READER_IMPL(XYFitCurve, const AbstractColumn*, xErrorColumn, xErrorColumn) BASIC_SHARED_D_READER_IMPL(XYFitCurve, const AbstractColumn*, yErrorColumn, yErrorColumn) const QString& XYFitCurve::xErrorColumnPath() const { Q_D(const XYFitCurve); return d->xErrorColumnPath; } const QString& XYFitCurve::yErrorColumnPath() const { Q_D(const XYFitCurve); return d->yErrorColumnPath; } BASIC_SHARED_D_READER_IMPL(XYFitCurve, XYFitCurve::FitData, fitData, fitData) const XYFitCurve::FitResult& XYFitCurve::fitResult() const { Q_D(const XYFitCurve); return d->fitResult; } //############################################################################## //################# setter methods and undo commands ########################## //############################################################################## STD_SETTER_CMD_IMPL_S(XYFitCurve, SetXErrorColumn, const AbstractColumn*, xErrorColumn) void XYFitCurve::setXErrorColumn(const AbstractColumn* column) { Q_D(XYFitCurve); if (column != d->xErrorColumn) { exec(new XYFitCurveSetXErrorColumnCmd(d, column, ki18n("%1: assign x-error"))); handleSourceDataChanged(); if (column) { connect(column, &AbstractColumn::dataChanged, this, [=](){ handleSourceDataChanged(); }); //TODO disconnect on undo } } } STD_SETTER_CMD_IMPL_S(XYFitCurve, SetYErrorColumn, const AbstractColumn*, yErrorColumn) void XYFitCurve::setYErrorColumn(const AbstractColumn* column) { Q_D(XYFitCurve); if (column != d->yErrorColumn) { exec(new XYFitCurveSetYErrorColumnCmd(d, column, ki18n("%1: assign y-error"))); handleSourceDataChanged(); if (column) { connect(column, &AbstractColumn::dataChanged, this, [=](){ handleSourceDataChanged(); }); //TODO disconnect on undo } } } // do not recalculate (allow preview) //STD_SETTER_CMD_IMPL_F_S(XYFitCurve, SetFitData, XYFitCurve::FitData, fitData, recalculate) STD_SETTER_CMD_IMPL_S(XYFitCurve, SetFitData, XYFitCurve::FitData, fitData) void XYFitCurve::setFitData(const XYFitCurve::FitData& fitData) { Q_D(XYFitCurve); exec(new XYFitCurveSetFitDataCmd(d, fitData, ki18n("%1: set fit options and perform the fit"))); } //############################################################################## //######################### Private implementation ############################# //############################################################################## -XYFitCurvePrivate::XYFitCurvePrivate(XYFitCurve* owner) : XYAnalysisCurvePrivate(owner), - xErrorColumn(nullptr), yErrorColumn(nullptr), - residualsColumn(nullptr), - residualsVector(nullptr), - q(owner) { -} +XYFitCurvePrivate::XYFitCurvePrivate(XYFitCurve* owner) : XYAnalysisCurvePrivate(owner), q(owner) {} //no need to delete xColumn and yColumn, they are deleted //when the parent aspect is removed XYFitCurvePrivate::~XYFitCurvePrivate() = default; // data structure to pass parameter to fit functions struct data { size_t n; //number of data points double* x; //pointer to the vector with x-data values double* y; //pointer to the vector with y-data values double* weight; //pointer to the vector with weight values nsl_fit_model_category modelCategory; int modelType; int degree; QString* func; // string containing the definition of the model/function QStringList* paramNames; double* paramMin; // lower parameter limits double* paramMax; // upper parameter limits bool* paramFixed; // parameter fixed? }; /*! * \param paramValues vector containing current values of the fit parameters * \param params * \param f vector with the weighted residuals weight[i]*(Yi - y[i]) */ int func_f(const gsl_vector* paramValues, void* params, gsl_vector* f) { //DEBUG("func_f"); size_t n = ((struct data*)params)->n; double* x = ((struct data*)params)->x; double* y = ((struct data*)params)->y; double* weight = ((struct data*)params)->weight; nsl_fit_model_category modelCategory = ((struct data*)params)->modelCategory; unsigned int modelType = ((struct data*)params)->modelType; QByteArray funcba = ((struct data*)params)->func->toLatin1(); // a local byte array is needed! const char *func = funcba.constData(); // function to evaluate QStringList* paramNames = ((struct data*)params)->paramNames; double *min = ((struct data*)params)->paramMin; double *max = ((struct data*)params)->paramMax; // set current values of the parameters for (int i = 0; i < paramNames->size(); i++) { double v = gsl_vector_get(paramValues, (size_t)i); // bound values if limits are set QByteArray paramnameba = paramNames->at(i).toLatin1(); assign_variable(paramnameba.constData(), nsl_fit_map_bound(v, min[i], max[i])); QDEBUG("Parameter"<n; double* xVector = ((struct data*)params)->x; double* weight = ((struct data*)params)->weight; nsl_fit_model_category modelCategory = ((struct data*)params)->modelCategory; unsigned int modelType = ((struct data*)params)->modelType; unsigned int degree = ((struct data*)params)->degree; QStringList* paramNames = ((struct data*)params)->paramNames; double *min = ((struct data*)params)->paramMin; double *max = ((struct data*)params)->paramMax; bool *fixed = ((struct data*)params)->paramFixed; // calculate the Jacobian matrix: // Jacobian matrix J(i,j) = df_i / dx_j // where f_i = w_i*(Y_i - y_i), // Y_i = model and the x_j are the parameters double x; switch (modelCategory) { case nsl_fit_model_basic: switch (modelType) { case nsl_fit_model_polynomial: // Y(x) = c0 + c1*x + ... + cn*x^n for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < (unsigned int)paramNames->size(); ++j) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_polynomial_param_deriv(x, j, weight[i])); } } break; case nsl_fit_model_power: // Y(x) = a*x^b or Y(x) = a + b*x^c. if (degree == 1) { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_power1_param_deriv(j, x, a, b, weight[i])); } } } else if (degree == 2) { const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double c = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_power2_param_deriv(j, x, b, c, weight[i])); } } } break; case nsl_fit_model_exponential: { // Y(x) = a*exp(b*x) + c*exp(d*x) + ... double *p = new double[2*degree]; for (unsigned int i = 0; i < 2*degree; i++) p[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, i), min[i], max[i]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2*degree; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exponentialn_param_deriv(j, x, p, weight[i])); } } delete[] p; break; } case nsl_fit_model_inverse_exponential: { // Y(x) = a*(1-exp(b*x))+c const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_inverse_exponential_param_deriv(j, x, a, b, weight[i])); } } break; } case nsl_fit_model_fourier: { // Y(x) = a0 + (a1*cos(w*x) + b1*sin(w*x)) + ... + (an*cos(n*w*x) + bn*sin(n*w*x) //parameters: w, a0, a1, b1, ... an, bn double* a = new double[degree]; double* b = new double[degree]; double w = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); a[0] = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); b[0] = 0; for (unsigned int i = 1; i < degree; ++i) { a[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, 2*i), min[2*i], max[2*i]); b[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, 2*i+1), min[2*i+1], max[2*i+1]); } for (size_t i = 0; i < n; i++) { x = xVector[i]; double wd = 0; //first derivative with respect to the w parameter for (unsigned int j = 1; j < degree; ++j) { wd += -a[j]*j*x*sin(j*w*x) + b[j]*j*x*cos(j*w*x); } gsl_matrix_set(J, i, 0, weight[i]*wd); gsl_matrix_set(J, i, 1, weight[i]); for (unsigned int j = 1; j <= degree; ++j) { gsl_matrix_set(J, (size_t)i, (size_t)(2*j), nsl_fit_model_fourier_param_deriv(0, j, x, w, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(2*j+1), nsl_fit_model_fourier_param_deriv(1, j, x, w, weight[i])); } for (unsigned int j = 0; j <= 2*degree+1; j++) if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } delete[] a; delete[] b; break; } } break; case nsl_fit_model_peak: switch (modelType) { case nsl_fit_model_gaussian: case nsl_fit_model_lorentz: case nsl_fit_model_sech: case nsl_fit_model_logistic: for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < degree; ++j) { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 3*j), min[3*j], max[3*j]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 3*j+1), min[3*j+1], max[3*j+1]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3*j+2), min[3*j+2], max[3*j+2]); switch (modelType) { case nsl_fit_model_gaussian: gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_gaussian_param_deriv(0, x, a, s, mu, weight[i])); 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])); 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])); break; case nsl_fit_model_lorentz: // a,s,t gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_lorentz_param_deriv(0, x, a, s, mu, weight[i])); 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])); 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])); break; case nsl_fit_model_sech: gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_sech_param_deriv(0, x, a, s, mu, weight[i])); 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])); 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])); break; case nsl_fit_model_logistic: gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_logistic_param_deriv(0, x, a, s, mu, weight[i])); 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])); 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])); break; } } for (unsigned int j = 0; j < 3*degree; j++) if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } break; case nsl_fit_model_voigt: for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < degree; ++j) { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j), min[4*j], max[4*j]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+1), min[4*j+1], max[4*j+1]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+2), min[4*j+2], max[4*j+2]); const double g = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+3), min[4*j+3], max[4*j+3]); 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])); 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])); 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])); 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])); } for (unsigned int j = 0; j < 4*degree; j++) if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } break; case nsl_fit_model_pseudovoigt1: for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < degree; ++j) { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j), min[4*j], max[4*j]); const double eta = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+1), min[4*j+1], max[4*j+1]); const double w = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+2), min[4*j+2], max[4*j+2]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+3), min[4*j+3], max[4*j+3]); gsl_matrix_set(J, (size_t)i, (size_t)(4*j), nsl_fit_model_voigt_param_deriv(0, x, a, eta, w, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(4*j+1), nsl_fit_model_voigt_param_deriv(1, x, a, eta, w, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(4*j+2), nsl_fit_model_voigt_param_deriv(2, x, a, eta, w, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(4*j+3), nsl_fit_model_voigt_param_deriv(3, x, a, eta, w, mu, weight[i])); } for (unsigned int j = 0; j < 4*degree; j++) if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } break; } break; case nsl_fit_model_growth: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) { gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } else { switch (modelType) { case nsl_fit_model_atan: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_atan_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_tanh: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_tanh_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_algebraic_sigmoid: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_algebraic_sigmoid_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_sigmoid: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_sigmoid_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_erf: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_erf_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_hill: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_hill_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_gompertz: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gompertz_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_gudermann: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gudermann_param_deriv(j, x, a, mu, s, weight[i])); break; } } } } } break; case nsl_fit_model_distribution: switch (modelType) { case nsl_sf_stats_gaussian: case nsl_sf_stats_exponential: case nsl_sf_stats_laplace: case nsl_sf_stats_cauchy_lorentz: case nsl_sf_stats_rayleigh_tail: case nsl_sf_stats_lognormal: case nsl_sf_stats_logistic: case nsl_sf_stats_sech: case nsl_sf_stats_levy: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { switch (modelType) { case nsl_sf_stats_gaussian: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gaussian_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_exponential: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exponential_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_laplace: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_laplace_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_cauchy_lorentz: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_lorentz_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_rayleigh_tail: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_rayleigh_tail_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_lognormal: if (x > 0) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_lognormal_param_deriv(j, x, a, s, mu, weight[i])); else gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); break; case nsl_sf_stats_logistic: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_logistic_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_sech: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_sech_dist_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_levy: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_levy_param_deriv(j, x, a, s, mu, weight[i])); break; } } } } break; } case nsl_sf_stats_gaussian_tail: { const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gaussian_tail_param_deriv(j, x, A, s, a, mu, weight[i])); } } break; } case nsl_sf_stats_exponential_power: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exp_pow_param_deriv(j, x, a, s, b, mu, weight[i])); } } break; } case nsl_sf_stats_rayleigh: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_rayleigh_param_deriv(j, x, a, s, weight[i])); } } break; } case nsl_sf_stats_gamma: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double k = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double t = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gamma_param_deriv(j, x, a, k, t, weight[i])); } } break; } case nsl_sf_stats_flat: { const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_flat_param_deriv(j, x, A, b, a, weight[i])); } } break; } case nsl_sf_stats_chi_squared: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double nu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_chi_square_param_deriv(j, x, a, nu, weight[i])); } } break; } case nsl_sf_stats_tdist: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double nu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_students_t_param_deriv(j, x, a, nu, weight[i])); } } break; } case nsl_sf_stats_fdist: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double n1 = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double n2 = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_fdist_param_deriv(j, x, a, n1, n2, weight[i])); } } break; } case nsl_sf_stats_beta: case nsl_sf_stats_pareto: { const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { switch (modelType) { case nsl_sf_stats_beta: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_beta_param_deriv(j, x, A, a, b, weight[i])); break; case nsl_sf_stats_pareto: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_pareto_param_deriv(j, x, A, a, b, weight[i])); break; } } } } break; } case nsl_sf_stats_weibull: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double k = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double l = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { if (x > 0) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_weibull_param_deriv(j, x, a, k, l, mu, weight[i])); else gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } } } break; } case nsl_sf_stats_gumbel1: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gumbel1_param_deriv(j, x, a, s, mu, b, weight[i])); } } break; } case nsl_sf_stats_gumbel2: { const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gumbel2_param_deriv(j, x, A, a, b, mu, weight[i])); } } break; } case nsl_sf_stats_poisson: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double l = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_poisson_param_deriv(j, x, a, l, weight[i])); } } break; } case nsl_sf_stats_maxwell_boltzmann: { // Y(x) = a*sqrt(2/pi) * x^2/s^3 * exp(-(x/s)^2/2) const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_maxwell_param_deriv(j, x, a, s, weight[i])); } } break; } case nsl_sf_stats_frechet: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double g = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_frechet_param_deriv(j, x, a, g, s, mu, weight[i])); } } break; } case nsl_sf_stats_landau: { // const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); for (size_t i = 0; i < n; i++) { x = xVector[i]; if (fixed[0]) gsl_matrix_set(J, (size_t)i, 0, 0.); else gsl_matrix_set(J, (size_t)i, 0, nsl_fit_model_landau_param_deriv(0, x, weight[i])); } break; } case nsl_sf_stats_binomial: case nsl_sf_stats_negative_binomial: case nsl_sf_stats_pascal: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double p = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double N = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { switch (modelType) { case nsl_sf_stats_binomial: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_binomial_param_deriv(j, x, a, p, N, weight[i])); break; case nsl_sf_stats_negative_binomial: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_negative_binomial_param_deriv(j, x, a, p, N, weight[i])); break; case nsl_sf_stats_pascal: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_pascal_param_deriv(j, x, a, p, N, weight[i])); break; } } } } break; } case nsl_sf_stats_geometric: case nsl_sf_stats_logarithmic: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double p = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { switch (modelType) { case nsl_sf_stats_geometric: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_geometric_param_deriv(j, x, a, p, weight[i])); break; case nsl_sf_stats_logarithmic: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_logarithmic_param_deriv(j, x, a, p, weight[i])); break; } } } } break; } case nsl_sf_stats_hypergeometric: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double n1 = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double n2 = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double t = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_hypergeometric_param_deriv(j, x, a, n1, n2, t, weight[i])); } } break; } // unused distributions case nsl_sf_stats_levy_alpha_stable: case nsl_sf_stats_levy_skew_alpha_stable: case nsl_sf_stats_bernoulli: break; } break; case nsl_fit_model_custom: QByteArray funcba = ((struct data*)params)->func->toLatin1(); const char* func = funcba.data(); QByteArray nameba; double value; const unsigned int np = paramNames->size(); for (size_t i = 0; i < n; i++) { x = xVector[i]; assign_variable("x", x); for (unsigned int j = 0; j < np; j++) { for (unsigned int k = 0; k < np; k++) { if (k != j) { nameba = paramNames->at(k).toLatin1(); value = nsl_fit_map_bound(gsl_vector_get(paramValues, k), min[k], max[k]); assign_variable(nameba.data(), value); } } nameba = paramNames->at(j).toLatin1(); const char *name = nameba.data(); value = nsl_fit_map_bound(gsl_vector_get(paramValues, j), min[j], max[j]); assign_variable(name, value); const double f_p = parse(func); double eps = 1.e-9; if (std::abs(f_p) > 0) eps *= std::abs(f_p); // scale step size with function value value += eps; assign_variable(name, value); const double f_pdp = parse(func); // DEBUG("evaluate deriv"<* >(xColumn->data()); yVector = static_cast* >(yColumn->data()); residualsVector = static_cast* >(residualsColumn->data()); xColumn->setHidden(true); q->addChild(xColumn); yColumn->setHidden(true); q->addChild(yColumn); q->addChild(residualsColumn); q->setUndoAware(false); q->setXColumn(xColumn); q->setYColumn(yColumn); q->setUndoAware(true); } else { xVector->clear(); yVector->clear(); residualsVector->clear(); } } void XYFitCurvePrivate::recalculate() { DEBUG("XYFitCurvePrivate::recalculate()"); + QElapsedTimer timer; timer.start(); // prepare source data columns const AbstractColumn* tmpXDataColumn = nullptr; const AbstractColumn* tmpYDataColumn = nullptr; if (dataSourceType == XYAnalysisCurve::DataSourceSpreadsheet) { DEBUG(" spreadsheet columns as data source"); tmpXDataColumn = xDataColumn; tmpYDataColumn = yDataColumn; } else { DEBUG(" curve columns as data source"); tmpXDataColumn = dataSourceCurve->xColumn(); tmpYDataColumn = dataSourceCurve->yColumn(); } if (!tmpXDataColumn || !tmpYDataColumn) { DEBUG("ERROR: Preparing source data columns failed!"); return; } prepareResultColumns(); // clear the previous result fitResult = XYFitCurve::FitResult(); //fit settings const unsigned int maxIters = fitData.maxIterations; //maximal number of iterations const double delta = fitData.eps; //fit tolerance const unsigned int np = fitData.paramNames.size(); //number of fit parameters if (np == 0) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("Model has no parameters."); emit q->dataChanged(); sourceDataChangedSinceLastRecalc = false; return; } //determine the data source columns if (!tmpXDataColumn || !tmpYDataColumn) { emit q->dataChanged(); sourceDataChangedSinceLastRecalc = false; return; } //check column sizes if (tmpXDataColumn->rowCount() != tmpYDataColumn->rowCount()) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("Number of x and y data points must be equal."); emit q->dataChanged(); sourceDataChangedSinceLastRecalc = false; return; } if (yErrorColumn) { if (yErrorColumn->rowCount() < tmpXDataColumn->rowCount()) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("Not sufficient weight data points provided."); emit q->dataChanged(); sourceDataChangedSinceLastRecalc = false; return; } } //copy all valid data point for the fit to temporary vectors QVector xdataVector; QVector ydataVector; QVector xerrorVector; QVector yerrorVector; double xmin, xmax; if (fitData.autoRange) { xmin = tmpXDataColumn->minimum(); xmax = tmpXDataColumn->maximum(); } else { xmin = fitData.fitRange.first(); xmax = fitData.fitRange.last(); } DEBUG("fit range = " << xmin << " .. " << xmax); for (int row = 0; row < tmpXDataColumn->rowCount(); ++row) { //only copy those data where _all_ values (for x and y and errors, if given) are valid if (!std::isnan(tmpXDataColumn->valueAt(row)) && !std::isnan(tmpYDataColumn->valueAt(row)) && !tmpXDataColumn->isMasked(row) && !tmpYDataColumn->isMasked(row)) { // only when inside given range if (tmpXDataColumn->valueAt(row) >= xmin && tmpXDataColumn->valueAt(row) <= xmax) { if (dataSourceType == XYAnalysisCurve::DataSourceCurve || (!xErrorColumn && !yErrorColumn) || !fitData.useDataErrors) { // x-y xdataVector.append(tmpXDataColumn->valueAt(row)); ydataVector.append(tmpYDataColumn->valueAt(row)); } else if (!xErrorColumn) { // x-y-dy if (!std::isnan(yErrorColumn->valueAt(row))) { xdataVector.append(tmpXDataColumn->valueAt(row)); ydataVector.append(tmpYDataColumn->valueAt(row)); yerrorVector.append(yErrorColumn->valueAt(row)); } } else { // x-y-dx-dy if (!std::isnan(xErrorColumn->valueAt(row)) && !std::isnan(yErrorColumn->valueAt(row))) { xdataVector.append(tmpXDataColumn->valueAt(row)); ydataVector.append(tmpYDataColumn->valueAt(row)); xerrorVector.append(xErrorColumn->valueAt(row)); yerrorVector.append(yErrorColumn->valueAt(row)); } } } } } //number of data points to fit const size_t n = xdataVector.size(); DEBUG("number of data points: " << n); if (n == 0) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("No data points available."); emit q->dataChanged(); sourceDataChangedSinceLastRecalc = false; return; } if (n < np) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("The number of data points (%1) must be greater than or equal to the number of parameters (%2).", n, np); emit q->dataChanged(); sourceDataChangedSinceLastRecalc = false; return; } double* xdata = xdataVector.data(); double* ydata = ydataVector.data(); double* xerror = xerrorVector.data(); // size may be 0 double* yerror = yerrorVector.data(); // size may be 0 DEBUG("x errors: " << xerrorVector.size()); DEBUG("y errors: " << yerrorVector.size()); double* weight = new double[n]; for (size_t i = 0; i < n; i++) weight[i] = 1.; switch (fitData.yWeightsType) { case nsl_fit_weight_no: case nsl_fit_weight_statistical_fit: case nsl_fit_weight_relative_fit: break; case nsl_fit_weight_instrumental: // yerror are sigmas for (int i = 0; i < (int)n; i++) if (i < yerrorVector.size()) weight[i] = 1./gsl_pow_2(yerror[i]); break; case nsl_fit_weight_direct: // yerror are weights for (int i = 0; i < (int)n; i++) if (i < yerrorVector.size()) weight[i] = yerror[i]; break; case nsl_fit_weight_inverse: // yerror are inverse weights for (int i = 0; i < (int)n; i++) if (i < yerrorVector.size()) weight[i] = 1./yerror[i]; break; case nsl_fit_weight_statistical: for (int i = 0; i < (int)n; i++) weight[i] = 1./ydata[i]; break; case nsl_fit_weight_relative: for (int i = 0; i < (int)n; i++) weight[i] = 1./gsl_pow_2(ydata[i]); break; } /////////////////////// GSL >= 2 has a complete new interface! But the old one is still supported. /////////////////////////// // GSL >= 2 : "the 'fdf' field of gsl_multifit_function_fdf is now deprecated and does not need to be specified for nonlinear least squares problems" unsigned int nf = 0; // number of fixed parameter for (unsigned int i = 0; i < np; i++) { const bool fixed = fitData.paramFixed.data()[i]; if (fixed) nf++; DEBUG("parameter " << i << " fixed: " << fixed); } //function to fit gsl_multifit_function_fdf f; DEBUG("model = " << fitData.model.toStdString()); struct data params = {n, xdata, ydata, weight, fitData.modelCategory, fitData.modelType, fitData.degree, &fitData.model, &fitData.paramNames, fitData.paramLowerLimits.data(), fitData.paramUpperLimits.data(), fitData.paramFixed.data()}; f.f = &func_f; f.df = &func_df; f.fdf = &func_fdf; f.n = n; f.p = np; f.params = ¶ms; DEBUG("initialize the derivative solver (using Levenberg-Marquardt robust solver)"); const gsl_multifit_fdfsolver_type* T = gsl_multifit_fdfsolver_lmsder; gsl_multifit_fdfsolver* s = gsl_multifit_fdfsolver_alloc(T, n, np); DEBUG("set start values"); double* x_init = fitData.paramStartValues.data(); double* x_min = fitData.paramLowerLimits.data(); double* x_max = fitData.paramUpperLimits.data(); DEBUG("scale start values if limits are set"); for (unsigned int i = 0; i < np; i++) x_init[i] = nsl_fit_map_unbound(x_init[i], x_min[i], x_max[i]); DEBUG(" DONE"); gsl_vector_view x = gsl_vector_view_array(x_init, np); DEBUG("Turning off GSL error handler to avoid overflow/underflow"); gsl_set_error_handler_off(); DEBUG("Initialize solver with function f and initial guess x"); gsl_multifit_fdfsolver_set(s, &f, &x.vector); DEBUG("Iterate ..."); int status; unsigned int iter = 0; fitResult.solverOutput.clear(); writeSolverState(s); do { iter++; DEBUG(" iter " << iter); // update weights for Y-depending weights (using function values from residuals) if (fitData.yWeightsType == nsl_fit_weight_statistical_fit) { for (size_t i = 0; i < n; i++) weight[i] = 1./(gsl_vector_get(s->f, i)/sqrt(weight[i]) + ydata[i]); // 1/Y_i } else if (fitData.yWeightsType == nsl_fit_weight_relative_fit) { for (size_t i = 0; i < n; i++) weight[i] = 1./gsl_pow_2(gsl_vector_get(s->f, i)/sqrt(weight[i]) + ydata[i]); // 1/Y_i^2 } DEBUG(" run fdfsolver_iterate"); status = gsl_multifit_fdfsolver_iterate(s); DEBUG(" fdfsolver_iterate DONE"); writeSolverState(s); if (status) { DEBUG("iter " << iter << ", status = " << gsl_strerror(status)); break; } status = gsl_multifit_test_delta(s->dx, s->x, delta, delta); DEBUG(" iter " << iter << ", test status = " << status); } while (status == GSL_CONTINUE && iter < maxIters); // second run for x-error fitting if (xerrorVector.size() > 0) { DEBUG("Rerun fit with x errors"); unsigned int iter2 = 0; double chi = 0, chiOld = 0; double *fun = new double[n]; do { iter2++; chiOld = chi; //printf("iter2 = %d\n", iter2); // calculate function from residuals for (size_t i = 0; i < n; i++) fun[i] = gsl_vector_get(s->f, i) * 1./sqrt(weight[i]) + ydata[i]; // calculate weight[i] for (size_t i = 0; i < n; i++) { // calculate df[i] size_t index = i-1; if (i == 0) index = i; if (i == n-1) index = i-2; double df = (fun[index+1] - fun[index])/(xdata[index+1] - xdata[index]); //printf("df = %g\n", df); double sigmasq = 1.; switch (fitData.xWeightsType) { // x-error type: f'(x)^2*s_x^2 = f'(x)/w_x case nsl_fit_weight_no: break; case nsl_fit_weight_direct: // xerror = w_x sigmasq = df*df/xerror[i]; break; case nsl_fit_weight_instrumental: // xerror = s_x sigmasq = df*df*xerror[i]*xerror[i]; break; case nsl_fit_weight_inverse: // xerror = 1/w_x = s_x^2 sigmasq = df*df*xerror[i]; break; case nsl_fit_weight_statistical: // s_x^2 = 1/w_x = x sigmasq = xdata[i]; break; case nsl_fit_weight_relative: // s_x^2 = 1/w_x = x^2 sigmasq = xdata[i]*xdata[i]; break; case nsl_fit_weight_statistical_fit: // unused case nsl_fit_weight_relative_fit: break; } if (yerrorVector.size() > 0) { switch (fitData.yWeightsType) { // y-error types: s_y^2 = 1/w_y case nsl_fit_weight_no: break; case nsl_fit_weight_direct: // yerror = w_y sigmasq += 1./yerror[i]; break; case nsl_fit_weight_instrumental: // yerror = s_y sigmasq += yerror[i]*yerror[i]; break; case nsl_fit_weight_inverse: // yerror = 1/w_y sigmasq += yerror[i]; break; case nsl_fit_weight_statistical: // unused case nsl_fit_weight_relative: break; case nsl_fit_weight_statistical_fit: // s_y^2 = 1/w_y = Y_i sigmasq += fun[i]; break; case nsl_fit_weight_relative_fit: // s_y^2 = 1/w_y = Y_i^2 sigmasq += fun[i]*fun[i]; break; } } //printf ("sigma[%d] = %g\n", i, sqrt(sigmasq)); weight[i] = 1./sigmasq; } // update weights gsl_multifit_fdfsolver_set(s, &f, &x.vector); do { // fit iter++; writeSolverState(s); status = gsl_multifit_fdfsolver_iterate (s); //printf ("status = %s\n", gsl_strerror (status)); if (status) { DEBUG("iter " << iter << ", status = " << gsl_strerror(status)); break; } status = gsl_multifit_test_delta(s->dx, s->x, delta, delta); } while (status == GSL_CONTINUE && iter < maxIters); chi = gsl_blas_dnrm2(s->f); //printf("chi = %.12g (dchi = %g)\n", chi, fabs(chi-chiOld)); } while (iter2 < maxIters && fabs(chi-chiOld) > fitData.eps); delete[] fun; } delete[] weight; // unscale start parameter for (unsigned int i = 0; i < np; i++) x_init[i] = nsl_fit_map_bound(x_init[i], x_min[i], x_max[i]); //get the covariance matrix //TODO: scale the Jacobian when limits are used before constructing the covar matrix? gsl_matrix* covar = gsl_matrix_alloc(np, np); #if GSL_MAJOR_VERSION >= 2 // the Jacobian is not part of the solver anymore gsl_matrix *J = gsl_matrix_alloc(s->fdf->n, s->fdf->p); gsl_multifit_fdfsolver_jac(s, J); gsl_multifit_covar(J, 0.0, covar); #else gsl_multifit_covar(s->J, 0.0, covar); #endif //write the result fitResult.available = true; fitResult.valid = true; fitResult.status = gslErrorToString(status); fitResult.iterations = iter; fitResult.dof = n - (np - nf); // samples - (parameter - fixed parameter) //gsl_blas_dnrm2() - computes the Euclidian norm (||r||_2 = \sqrt {\sum r_i^2}) of the vector with the elements weight[i]*(Yi - y[i]) //gsl_blas_dasum() - computes the absolute sum \sum |r_i| of the elements of the vector with the elements weight[i]*(Yi - y[i]) fitResult.sse = gsl_pow_2(gsl_blas_dnrm2(s->f)); if (fitResult.dof != 0) { fitResult.rms = fitResult.sse/fitResult.dof; fitResult.rsd = sqrt(fitResult.rms); } fitResult.mse = fitResult.sse/n; fitResult.rmse = sqrt(fitResult.mse); fitResult.mae = gsl_blas_dasum(s->f)/n; // SST needed for coefficient of determination, R-squared fitResult.sst = gsl_stats_tss(ydata, 1, n); // for a linear model without intercept R-squared is calculated differently // see https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-does-summary_0028_0029-report-strange-results-for-the-R_005e2-estimate-when-I-fit-a-linear-model-with-no-intercept_003f if (fitData.modelCategory == nsl_fit_model_basic && fitData.modelType == nsl_fit_model_polynomial && fitData.degree == 1 && x_init[0] == 0) { DEBUG("Using alternative R^2 for linear model without intercept"); fitResult.sst = gsl_stats_tss_m(ydata, 1, n, 0); } if (fitResult.sst < fitResult.sse) { DEBUG("Using alternative R^2 since R^2 would be negative (probably custom model without intercept)"); fitResult.sst = gsl_stats_tss_m(ydata, 1, n, 0); } fitResult.rsquare = nsl_stats_rsquare(fitResult.sse, fitResult.sst); fitResult.rsquareAdj = nsl_stats_rsquareAdj(fitResult.rsquare, np, fitResult.dof, 1); fitResult.chisq_p = nsl_stats_chisq_p(fitResult.sse, fitResult.dof); fitResult.fdist_F = nsl_stats_fdist_F(fitResult.sst, fitResult.rms, np, 1); fitResult.fdist_p = nsl_stats_fdist_p(fitResult.fdist_F, np, fitResult.dof); fitResult.logLik = nsl_stats_logLik(fitResult.sse, n); fitResult.aic = nsl_stats_aic(fitResult.sse, n, np, 1); fitResult.bic = nsl_stats_bic(fitResult.sse, n, np, 1); //parameter values // GSL: const double c = GSL_MAX_DBL(1., sqrt(fitResult.rms)); // increase error for poor fit // NIST: const double c = sqrt(fitResult.rms); // increase error for poor fit, decrease for good fit const double c = sqrt(fitResult.rms); fitResult.paramValues.resize(np); fitResult.errorValues.resize(np); fitResult.tdist_tValues.resize(np); fitResult.tdist_pValues.resize(np); fitResult.tdist_marginValues.resize(np); for (unsigned int i = 0; i < np; i++) { // scale resulting values if they are bounded fitResult.paramValues[i] = nsl_fit_map_bound(gsl_vector_get(s->x, i), x_min[i], x_max[i]); // use results as start values if desired if (fitData.useResults) { fitData.paramStartValues.data()[i] = fitResult.paramValues[i]; DEBUG("saving parameter " << i << ": " << fitResult.paramValues[i] << ' ' << fitData.paramStartValues.data()[i]); } fitResult.errorValues[i] = c*sqrt(gsl_matrix_get(covar, i, i)); fitResult.tdist_tValues[i] = nsl_stats_tdist_t(fitResult.paramValues.at(i), fitResult.errorValues.at(i)); fitResult.tdist_pValues[i] = nsl_stats_tdist_p(fitResult.tdist_tValues.at(i), fitResult.dof); fitResult.tdist_marginValues[i] = nsl_stats_tdist_margin(0.05, fitResult.dof, fitResult.errorValues.at(i)); } // fill residuals vector. To get residuals on the correct x values, fill the rest with zeros. residualsVector->resize(tmpXDataColumn->rowCount()); if (fitData.autoRange) { // evaluate full range of residuals xVector->resize(tmpXDataColumn->rowCount()); for (int i = 0; i < tmpXDataColumn->rowCount(); i++) (*xVector)[i] = tmpXDataColumn->valueAt(i); ExpressionParser* parser = ExpressionParser::getInstance(); bool rc = parser->evaluateCartesian(fitData.model, xVector, residualsVector, fitData.paramNames, fitResult.paramValues); for (int i = 0; i < tmpXDataColumn->rowCount(); i++) (*residualsVector)[i] = tmpYDataColumn->valueAt(i) - (*residualsVector)[i]; if (!rc) residualsVector->clear(); } else { // only selected range size_t j = 0; for (int i = 0; i < tmpXDataColumn->rowCount(); i++) { if (tmpXDataColumn->valueAt(i) >= xmin && tmpXDataColumn->valueAt(i) <= xmax) residualsVector->data()[i] = - gsl_vector_get(s->f, j++); else // outside range residualsVector->data()[i] = 0; } } residualsColumn->setChanged(); //free resources gsl_multifit_fdfsolver_free(s); gsl_matrix_free(covar); //calculate the fit function (vectors) evaluate(); fitResult.elapsedTime = timer.elapsed(); sourceDataChangedSinceLastRecalc = false; } /* evaluate fit function */ void XYFitCurvePrivate::evaluate(bool preview) { DEBUG("XYFitCurvePrivate::evaluate() preview = " << preview); // prepare source data columns const AbstractColumn* tmpXDataColumn = nullptr; if (dataSourceType == XYAnalysisCurve::DataSourceSpreadsheet) { DEBUG(" spreadsheet columns as data source"); tmpXDataColumn = xDataColumn; } else { DEBUG(" curve columns as data source"); tmpXDataColumn = dataSourceCurve->xColumn(); } if (!tmpXDataColumn) { DEBUG("ERROR: Preparing source data column failed!"); return; } prepareResultColumns(); if (!xVector || !yVector) { DEBUG(" xVector or yVector not defined!"); return; } ExpressionParser* parser = ExpressionParser::getInstance(); double xmin, xmax; if (fitData.autoEvalRange) { // evaluate fit on full data range xmin = tmpXDataColumn->minimum(); xmax = tmpXDataColumn->maximum(); } else { // use given range for evaluation xmin = fitData.evalRange.first(); xmax = fitData.evalRange.last(); } DEBUG(" eval range = " << xmin << " .. " << xmax); xVector->resize((int)fitData.evaluatedPoints); yVector->resize((int)fitData.evaluatedPoints); DEBUG(" vector size = " << xVector->size()); QVector paramValues = fitResult.paramValues; if (preview) // results not available yet paramValues = fitData.paramStartValues; bool rc = parser->evaluateCartesian(fitData.model, QString::number(xmin), QString::number(xmax), (int)fitData.evaluatedPoints, xVector, yVector, fitData.paramNames, paramValues); if (!rc) { xVector->clear(); yVector->clear(); residualsVector->clear(); } emit q->dataChanged(); } /*! * writes out the current state of the solver \c s */ void XYFitCurvePrivate::writeSolverState(gsl_multifit_fdfsolver* s) { QString state; //current parameter values, semicolon separated double* min = fitData.paramLowerLimits.data(); double* max = fitData.paramUpperLimits.data(); for (int i = 0; i < fitData.paramNames.size(); ++i) { const double x = gsl_vector_get(s->x, i); // map parameter if bounded state += QString::number(nsl_fit_map_bound(x, min[i], max[i])) + '\t'; } //current value of the chi2-function state += QString::number(gsl_pow_2(gsl_blas_dnrm2(s->f))); state += ';'; DEBUG(" chi = " << gsl_pow_2(gsl_blas_dnrm2(s->f))); fitResult.solverOutput += state; } //############################################################################## //################## Serialization/Deserialization ########################### //############################################################################## //! Save as XML void XYFitCurve::save(QXmlStreamWriter* writer) const { Q_D(const XYFitCurve); writer->writeStartElement("xyFitCurve"); //write the base class XYAnalysisCurve::save(writer); //write xy-fit-curve specific information //fit data - only save model expression and parameter names for custom model, otherwise they are set in XYFitCurve::initFitData() writer->writeStartElement("fitData"); WRITE_COLUMN(d->xErrorColumn, xErrorColumn); WRITE_COLUMN(d->yErrorColumn, yErrorColumn); writer->writeAttribute("autoRange", QString::number(d->fitData.autoRange)); writer->writeAttribute("fitRangeMin", QString::number(d->fitData.fitRange.first(), 'g', 15)); writer->writeAttribute("fitRangeMax", QString::number(d->fitData.fitRange.last(), 'g', 15)); writer->writeAttribute("modelCategory", QString::number(d->fitData.modelCategory)); writer->writeAttribute("modelType", QString::number(d->fitData.modelType)); writer->writeAttribute("xWeightsType", QString::number(d->fitData.xWeightsType)); writer->writeAttribute("weightsType", QString::number(d->fitData.yWeightsType)); writer->writeAttribute("degree", QString::number(d->fitData.degree)); if (d->fitData.modelCategory == nsl_fit_model_custom) writer->writeAttribute("model", d->fitData.model); writer->writeAttribute("maxIterations", QString::number(d->fitData.maxIterations)); writer->writeAttribute("eps", QString::number(d->fitData.eps, 'g', 15)); writer->writeAttribute("evaluatedPoints", QString::number(d->fitData.evaluatedPoints)); writer->writeAttribute("autoEvalRange", QString::number(d->fitData.autoEvalRange)); writer->writeAttribute("useDataErrors", QString::number(d->fitData.useDataErrors)); writer->writeAttribute("useResults", QString::number(d->fitData.useResults)); writer->writeAttribute("previewEnabled", QString::number(d->fitData.previewEnabled)); if (d->fitData.modelCategory == nsl_fit_model_custom) { writer->writeStartElement("paramNames"); foreach (const QString &name, d->fitData.paramNames) writer->writeTextElement("name", name); writer->writeEndElement(); } writer->writeStartElement("paramStartValues"); foreach (const double &value, d->fitData.paramStartValues) writer->writeTextElement("startValue", QString::number(value, 'g', 15)); writer->writeEndElement(); // use 16 digits to handle -DBL_MAX writer->writeStartElement("paramLowerLimits"); foreach (const double &limit, d->fitData.paramLowerLimits) writer->writeTextElement("lowerLimit", QString::number(limit, 'g', 16)); writer->writeEndElement(); // use 16 digits to handle DBL_MAX writer->writeStartElement("paramUpperLimits"); foreach (const double &limit, d->fitData.paramUpperLimits) writer->writeTextElement("upperLimit", QString::number(limit, 'g', 16)); writer->writeEndElement(); writer->writeStartElement("paramFixed"); foreach (const double &fixed, d->fitData.paramFixed) writer->writeTextElement("fixed", QString::number(fixed)); writer->writeEndElement(); writer->writeEndElement(); //"fitData" //fit results (generated columns and goodness of the fit) writer->writeStartElement("fitResult"); writer->writeAttribute("available", QString::number(d->fitResult.available)); writer->writeAttribute("valid", QString::number(d->fitResult.valid)); writer->writeAttribute("status", d->fitResult.status); writer->writeAttribute("iterations", QString::number(d->fitResult.iterations)); writer->writeAttribute("time", QString::number(d->fitResult.elapsedTime)); writer->writeAttribute("dof", QString::number(d->fitResult.dof)); writer->writeAttribute("sse", QString::number(d->fitResult.sse, 'g', 15)); writer->writeAttribute("sst", QString::number(d->fitResult.sst, 'g', 15)); writer->writeAttribute("rms", QString::number(d->fitResult.rms, 'g', 15)); writer->writeAttribute("rsd", QString::number(d->fitResult.rsd, 'g', 15)); writer->writeAttribute("mse", QString::number(d->fitResult.mse, 'g', 15)); writer->writeAttribute("rmse", QString::number(d->fitResult.rmse, 'g', 15)); writer->writeAttribute("mae", QString::number(d->fitResult.mae, 'g', 15)); writer->writeAttribute("rsquare", QString::number(d->fitResult.rsquare, 'g', 15)); writer->writeAttribute("rsquareAdj", QString::number(d->fitResult.rsquareAdj, 'g', 15)); writer->writeAttribute("chisq_p", QString::number(d->fitResult.chisq_p, 'g', 15)); writer->writeAttribute("fdist_F", QString::number(d->fitResult.fdist_F, 'g', 15)); writer->writeAttribute("fdist_p", QString::number(d->fitResult.fdist_p, 'g', 15)); writer->writeAttribute("aic", QString::number(d->fitResult.aic, 'g', 15)); writer->writeAttribute("bic", QString::number(d->fitResult.bic, 'g', 15)); writer->writeAttribute("solverOutput", d->fitResult.solverOutput); writer->writeStartElement("paramValues"); foreach (const double &value, d->fitResult.paramValues) writer->writeTextElement("value", QString::number(value, 'g', 15)); writer->writeEndElement(); writer->writeStartElement("errorValues"); foreach (const double &value, d->fitResult.errorValues) writer->writeTextElement("error", QString::number(value, 'g', 15)); writer->writeEndElement(); writer->writeStartElement("tdist_tValues"); foreach (const double &value, d->fitResult.tdist_tValues) writer->writeTextElement("tdist_t", QString::number(value, 'g', 15)); writer->writeEndElement(); writer->writeStartElement("tdist_pValues"); foreach (const double &value, d->fitResult.tdist_pValues) writer->writeTextElement("tdist_p", QString::number(value, 'g', 15)); writer->writeEndElement(); writer->writeStartElement("tdist_marginValues"); foreach (const double &value, d->fitResult.tdist_marginValues) writer->writeTextElement("tdist_margin", QString::number(value, 'g', 15)); writer->writeEndElement(); //save calculated columns if available if (d->xColumn && d->yColumn && d->residualsColumn) { d->xColumn->save(writer); d->yColumn->save(writer); d->residualsColumn->save(writer); } writer->writeEndElement(); //"fitResult" writer->writeEndElement(); //"xyFitCurve" } //! Load from XML bool XYFitCurve::load(XmlStreamReader* reader, bool preview) { DEBUG("XYFitCurve::load()"); Q_D(XYFitCurve); KLocalizedString attributeWarning = ki18n("Attribute '%1' missing or empty, default value is used"); QXmlStreamAttributes attribs; QString str; while (!reader->atEnd()) { reader->readNext(); if (reader->isEndElement() && reader->name() == "xyFitCurve") break; if (!reader->isStartElement()) continue; if (reader->name() == "xyAnalysisCurve") { if ( !XYAnalysisCurve::load(reader, preview) ) return false; } else if (!preview && reader->name() == "fitData") { attribs = reader->attributes(); READ_COLUMN(xErrorColumn); READ_COLUMN(yErrorColumn); READ_INT_VALUE("autoRange", fitData.autoRange, bool); READ_DOUBLE_VALUE("xRangeMin", fitData.fitRange.first()); // old name READ_DOUBLE_VALUE("xRangeMax", fitData.fitRange.last()); // old name READ_DOUBLE_VALUE("fitRangeMin", fitData.fitRange.first()); READ_DOUBLE_VALUE("fitRangeMax", fitData.fitRange.last()); READ_INT_VALUE("modelCategory", fitData.modelCategory, nsl_fit_model_category); READ_INT_VALUE("modelType", fitData.modelType, unsigned int); READ_INT_VALUE("xWeightsType", fitData.xWeightsType, nsl_fit_weight_type); READ_INT_VALUE("weightsType", fitData.yWeightsType, nsl_fit_weight_type); READ_INT_VALUE("degree", fitData.degree, int); if (d->fitData.modelCategory == nsl_fit_model_custom) { READ_STRING_VALUE("model", fitData.model); DEBUG("read model = " << d->fitData.model.toStdString()); } READ_INT_VALUE("maxIterations", fitData.maxIterations, int); READ_DOUBLE_VALUE("eps", fitData.eps); READ_INT_VALUE("fittedPoints", fitData.evaluatedPoints, size_t); // old name READ_INT_VALUE("evaluatedPoints", fitData.evaluatedPoints, size_t); READ_INT_VALUE("evaluateFullRange", fitData.autoEvalRange, bool); // old name READ_INT_VALUE("autoEvalRange", fitData.autoEvalRange, bool); READ_INT_VALUE("useDataErrors", fitData.useDataErrors, bool); READ_INT_VALUE("useResults", fitData.useResults, bool); READ_INT_VALUE("previewEnabled", fitData.previewEnabled, bool); //set the model expression and the parameter names (can be derived from the saved values for category, type and degree) XYFitCurve::initFitData(d->fitData); // remove default start values d->fitData.paramStartValues.clear(); } else if (!preview && reader->name() == "name") { // needed for custom model d->fitData.paramNames << reader->readElementText(); } else if (!preview && reader->name() == "startValue") { d->fitData.paramStartValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "fixed") { d->fitData.paramFixed << (bool)reader->readElementText().toInt(); } else if (!preview && reader->name() == "lowerLimit") { bool ok; double x = reader->readElementText().toDouble(&ok); if (ok) // -DBL_MAX results in conversion error d->fitData.paramLowerLimits << x; else d->fitData.paramLowerLimits << -std::numeric_limits::max(); } else if (!preview && reader->name() == "upperLimit") { bool ok; double x = reader->readElementText().toDouble(&ok); if (ok) // DBL_MAX results in conversion error d->fitData.paramUpperLimits << x; else d->fitData.paramUpperLimits << std::numeric_limits::max(); } else if (!preview && reader->name() == "value") { d->fitResult.paramValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "error") { d->fitResult.errorValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "tdist_t") { d->fitResult.tdist_tValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "tdist_p") { d->fitResult.tdist_pValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "tdist_margin") { d->fitResult.tdist_marginValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "fitResult") { attribs = reader->attributes(); READ_INT_VALUE("available", fitResult.available, int); READ_INT_VALUE("valid", fitResult.valid, int); READ_STRING_VALUE("status", fitResult.status); READ_INT_VALUE("iterations", fitResult.iterations, int); READ_INT_VALUE("time", fitResult.elapsedTime, int); READ_DOUBLE_VALUE("dof", fitResult.dof); READ_DOUBLE_VALUE("sse", fitResult.sse); READ_DOUBLE_VALUE("sst", fitResult.sst); READ_DOUBLE_VALUE("rms", fitResult.rms); READ_DOUBLE_VALUE("rsd", fitResult.rsd); READ_DOUBLE_VALUE("mse", fitResult.mse); READ_DOUBLE_VALUE("rmse", fitResult.rmse); READ_DOUBLE_VALUE("mae", fitResult.mae); READ_DOUBLE_VALUE("rsquare", fitResult.rsquare); READ_DOUBLE_VALUE("rsquareAdj", fitResult.rsquareAdj); READ_DOUBLE_VALUE("chisq_p", fitResult.chisq_p); READ_DOUBLE_VALUE("fdist_F", fitResult.fdist_F); READ_DOUBLE_VALUE("fdist_p", fitResult.fdist_p); READ_DOUBLE_VALUE("aic", fitResult.aic); READ_DOUBLE_VALUE("bic", fitResult.bic); READ_STRING_VALUE("solverOutput", fitResult.solverOutput); } else if (reader->name() == "column") { Column* column = new Column("", AbstractColumn::Numeric); if (!column->load(reader, preview)) { delete column; return false; } if (column->name() == "x") d->xColumn = column; else if (column->name() == "y") d->yColumn = column; else if (column->name() == "residuals") d->residualsColumn = column; } } if (d->fitData.paramNamesUtf8.isEmpty()) d->fitData.paramNamesUtf8 << d->fitData.paramNames; DEBUG("# params = " << d->fitData.paramNames.size()); if (preview) return true; // new fit model style (reset model type of old projects) if (d->fitData.modelCategory == nsl_fit_model_basic && d->fitData.modelType >= NSL_FIT_MODEL_BASIC_COUNT) { DEBUG("RESET old fit model"); d->fitData.modelType = 0; d->fitData.degree = 1; d->fitData.paramNames.clear(); d->fitData.paramNamesUtf8.clear(); // reset size of fields not touched by initFitData() d->fitData.paramStartValues.resize(2); d->fitData.paramFixed.resize(2); d->fitResult.paramValues.resize(2); d->fitResult.errorValues.resize(2); d->fitResult.tdist_tValues.resize(2); d->fitResult.tdist_pValues.resize(2); d->fitResult.tdist_marginValues.resize(2); } // not present in old projects int np = d->fitResult.paramValues.size(); if (d->fitResult.tdist_tValues.size() == 0) d->fitResult.tdist_tValues.resize(np); if (d->fitResult.tdist_pValues.size() == 0) d->fitResult.tdist_pValues.resize(np); if (d->fitResult.tdist_marginValues.size() == 0) d->fitResult.tdist_marginValues.resize(np); DEBUG("# start values = " << d->fitData.paramStartValues.size()); // wait for data to be read before using the pointers QThreadPool::globalInstance()->waitForDone(); if (d->xColumn && d->yColumn && d->residualsColumn) { d->xColumn->setHidden(true); addChild(d->xColumn); d->yColumn->setHidden(true); addChild(d->yColumn); addChild(d->residualsColumn); d->xVector = static_cast* >(d->xColumn->data()); d->yVector = static_cast* >(d->yColumn->data()); d->residualsVector = static_cast* >(d->residualsColumn->data()); XYCurve::d_ptr->xColumn = d->xColumn; XYCurve::d_ptr->yColumn = d->yColumn; } return true; } diff --git a/src/backend/worksheet/plots/cartesian/XYFitCurve.h b/src/backend/worksheet/plots/cartesian/XYFitCurve.h index 46211d93b..1655e4680 100644 --- a/src/backend/worksheet/plots/cartesian/XYFitCurve.h +++ b/src/backend/worksheet/plots/cartesian/XYFitCurve.h @@ -1,160 +1,145 @@ /*************************************************************************** File : XYFitCurve.h Project : LabPlot Description : A xy-curve defined by a fit model -------------------------------------------------------------------- Copyright : (C) 2014-2017 Alexander Semke (alexander.semke@web.de) Copyright : (C) 2016-2018 Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #ifndef XYFITCURVE_H #define XYFITCURVE_H #include "backend/worksheet/plots/cartesian/XYAnalysisCurve.h" #include "kdefrontend/spreadsheet/PlotDataDialog.h" //for PlotDataDialog::AnalysisAction. TODO: find a better place for this enum. extern "C" { #include "backend/nsl/nsl_fit.h" } class XYFitCurvePrivate; class XYFitCurve : public XYAnalysisCurve { Q_OBJECT public: struct FitData { - FitData() : modelCategory(nsl_fit_model_basic), modelType(0), - xWeightsType(nsl_fit_weight_no), - yWeightsType(nsl_fit_weight_no), - degree(1), - maxIterations(500), - eps(1.e-4), - evaluatedPoints(1000), - useDataErrors(true), - useResults(true), - previewEnabled(false), - autoRange(true), - autoEvalRange(true), - fitRange(2), - evalRange(2) {}; - - nsl_fit_model_category modelCategory; - int modelType; - nsl_fit_weight_type xWeightsType; - nsl_fit_weight_type yWeightsType; - int degree; + FitData() {}; + + nsl_fit_model_category modelCategory{nsl_fit_model_basic}; + int modelType{0}; + nsl_fit_weight_type xWeightsType{nsl_fit_weight_no}; + nsl_fit_weight_type yWeightsType{nsl_fit_weight_no}; + int degree{1}; QString model; QStringList paramNames; QStringList paramNamesUtf8; // Utf8 version of paramNames QVector paramStartValues; QVector paramLowerLimits; QVector paramUpperLimits; QVector paramFixed; - int maxIterations; - double eps; - size_t evaluatedPoints; - bool useDataErrors; // use given data errors when fitting (default) - bool useResults; // use results as new start values (default) - bool previewEnabled; // preview fit function with given start parameters - - bool autoRange; // use all data points? (default) - bool autoEvalRange; // evaluate fit function on full data range (default) - QVector fitRange; // x fit range - QVector evalRange; // x evaluation range + int maxIterations{500}; + double eps{1.e-4}; + size_t evaluatedPoints{1000}; + bool useDataErrors{true}; // use given data errors when fitting (default) + bool useResults{true}; // use results as new start values (default) + bool previewEnabled{false}; // preview fit function with given start parameters + + bool autoRange{true}; // use all data points? (default) + bool autoEvalRange{true}; // evaluate fit function on full data range (default) + QVector fitRange{0., 0.}; // x fit range + QVector evalRange{0., 0.}; // x evaluation range }; struct FitResult { - FitResult() : available(false), valid(false), iterations(0), elapsedTime(0), - dof(0), sse(0), sst(0), rms(0), rsd(0), mse(0), rmse(0), mae(0), rsquare(0), rsquareAdj(0), - chisq_p(0), fdist_F(0), fdist_p(0), logLik(0), aic(0), bic(0) {}; + FitResult() {}; - bool available; - bool valid; + bool available{false}; + bool valid{false}; QString status; - int iterations; - qint64 elapsedTime; - double dof; //degrees of freedom + int iterations{0}; + qint64 elapsedTime{0}; + double dof{0}; //degrees of freedom // residuals: r_i = y_i - Y_i - double sse; // sum of squared errors (SSE) / residual sum of squares (RSS) / sum of sq. residuals (SSR) / S = chi^2 = \sum_i^n r_i^2 - double sst; // total sum of squares (SST) = \sum_i^n (y_i - )^2 - double rms; // residual mean square / reduced chi^2 = SSE/dof - double rsd; // residual standard deviation = sqrt(SSE/dof) - double mse; // mean squared error = SSE/n - double rmse; // root-mean squared error = \sqrt(mse) - double mae; // mean absolute error = \sum_i^n |r_i| - double rsquare; - double rsquareAdj; - double chisq_p; // chi^2 distribution p-value - double fdist_F; // F distribution F-value - double fdist_p; // F distribution p-value - double logLik; // log likelihood - double aic; // Akaike information criterion - double bic; // Schwarz Bayesian information criterion + double sse{0}; // sum of squared errors (SSE) / residual sum of squares (RSS) / sum of sq. residuals (SSR) / S = chi^2 = \sum_i^n r_i^2 + double sst{0}; // total sum of squares (SST) = \sum_i^n (y_i - )^2 + double rms{0}; // residual mean square / reduced chi^2 = SSE/dof + double rsd{0}; // residual standard deviation = sqrt(SSE/dof) + double mse{0}; // mean squared error = SSE/n + double rmse{0}; // root-mean squared error = \sqrt(mse) + double mae{0}; // mean absolute error = \sum_i^n |r_i| + double rsquare{0}; + double rsquareAdj{0}; + double chisq_p{0}; // chi^2 distribution p-value + double fdist_F{0}; // F distribution F-value + double fdist_p{0}; // F distribution p-value + double logLik{0}; // log likelihood + double aic{0}; // Akaike information criterion + double bic{0}; // Schwarz Bayesian information criterion // see also http://www.originlab.com/doc/Origin-Help/NLFit-Algorithm QVector paramValues; QVector errorValues; QVector tdist_tValues; QVector tdist_pValues; QVector tdist_marginValues; QString solverOutput; }; explicit XYFitCurve(const QString& name); ~XYFitCurve() override; void recalculate() override; void evaluate(bool preview); void initFitData(PlotDataDialog::AnalysisAction); static void initFitData(XYFitCurve::FitData&); void initStartValues(const XYCurve*); static void initStartValues(XYFitCurve::FitData&, const XYCurve*); QIcon icon() const override; void save(QXmlStreamWriter*) const override; bool load(XmlStreamReader*, bool preview) override; POINTER_D_ACCESSOR_DECL(const AbstractColumn, xErrorColumn, XErrorColumn) POINTER_D_ACCESSOR_DECL(const AbstractColumn, yErrorColumn, YErrorColumn) const QString& xErrorColumnPath() const; const QString& yErrorColumnPath() const; CLASS_D_ACCESSOR_DECL(FitData, fitData, FitData) const FitResult& fitResult() const; typedef XYFitCurvePrivate Private; protected: XYFitCurve(const QString& name, XYFitCurvePrivate* dd); private: Q_DECLARE_PRIVATE(XYFitCurve) signals: void xErrorColumnChanged(const AbstractColumn*); void yErrorColumnChanged(const AbstractColumn*); void fitDataChanged(const XYFitCurve::FitData&); }; #endif diff --git a/src/backend/worksheet/plots/cartesian/XYFitCurvePrivate.h b/src/backend/worksheet/plots/cartesian/XYFitCurvePrivate.h index f0d2ae599..973ebf545 100644 --- a/src/backend/worksheet/plots/cartesian/XYFitCurvePrivate.h +++ b/src/backend/worksheet/plots/cartesian/XYFitCurvePrivate.h @@ -1,69 +1,69 @@ /*************************************************************************** File : XYFitCurvePrivate.h Project : LabPlot Description : Private members of XYFitCurve -------------------------------------------------------------------- Copyright : (C) 2014-2017 Alexander Semke (alexander.semke@web.de) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #ifndef XYFITCURVEPRIVATE_H #define XYFITCURVEPRIVATE_H #include "backend/worksheet/plots/cartesian/XYAnalysisCurvePrivate.h" #include "backend/worksheet/plots/cartesian/XYFitCurve.h" class XYFitCurve; class Column; extern "C" { #include } class XYFitCurvePrivate : public XYAnalysisCurvePrivate { public: explicit XYFitCurvePrivate(XYFitCurve*); ~XYFitCurvePrivate() override; void recalculate(); void evaluate(bool preview = false); - const AbstractColumn* xErrorColumn; //* residualsVector; + Column* residualsColumn{nullptr}; + QVector* residualsVector{nullptr}; XYFitCurve* const q; private: void prepareResultColumns(); void writeSolverState(gsl_multifit_fdfsolver*); }; #endif diff --git a/src/doc/coding_style.dox b/src/doc/coding_style.dox index 1a3a073d5..e2c6306f3 100644 --- a/src/doc/coding_style.dox +++ b/src/doc/coding_style.dox @@ -1,134 +1,136 @@ /**\page coding_style Coding Style The following rules are not used everywhere (yet), but are intended as guidelines for new code and eventually old code should be adapted as well. They apply to C++ and C code. The standards are C++11 and C99. \section files Files - Files use Unix-style line endings ('\\n'). - C++ source files use “.cpp” as extension, C source code use "*.c" and header files use “.h”. - The code is documented using Doxygen comments which are placed in the source files, not the header files. - Every file should be named exactly like the class inside and there should be only one class per file, with the exception of really short classes. Very short classes can be bundled in one file which then is named using all lower case letters. \section identifier Identifier names - Class names start with a capital letter and use CamelCase, acronyms in class names are use like normal words. Example: MySuperHtmlToPdfConverter - Function/method names start with a lower case letter and use CamelCase Example: doSomethingImportant() - Variable/object names start with a lower case letter and use CamelCase, underscores are used for special prefixes only. - Only private class member variables are prefixed with “m_” to distinguish them easily. d-pointer and UI-widgets are called d and ui, respectively, i.e. without prefix. - Property access methods use Qt style: property() and setProperty(), except for boolean properties (isVisible(), hasChanged()). Accessor functions (getter/setter) can be done using macros. - Avoid abbreviations, except for local counters and temporaries whose purpose is obvious. \section indent Indentation, spacing and line breaks - Tabs are used for indentation because they allow everyone to choose the indentation depth for him/herself. - Try to keep lines shorter than 100 characters, inserting line breaks as necessary and indent the following lines to improved readability. - included headers should be in order: own header, local header, Qt/KDE header, system header, extern header - Opening braces (‘{‘) are placed behind the statement and are preceded by a space. This also goes for function implementations, class, struct and namespace declarations, which are exceptions in other coding styles. Example: @code void MyClass::doSomething() { if (condition) { ... } ... } @endcode - Opening brackets (‘(‘) are preceded by a space in for/switch/if/while statements, but not for function calls. Example: @code if (condition) { doSomething(myData); ... } @endcode - For pointers or references, use a single space after ‘*’ or ‘&’ (i.e. specifier is bound to the data type not the name). Example: @code void doSomething(int* dataPointer, const QString& name); ... = static_cast(...) @endcode “public” and namespace enclosures are not indented. Example: @code class MyClass: public QObject { public: void doSomething(); @endcode “case” of switch is not indented. “default” should be present only if data type is not an enum. Example: @code switch (condition) { case 1: handleCaseOne(); break; case 2: { int i=0; ... break; } ... default: ... } @endcode - Each comma in a function call or semicolon in a for statement is followed by a space character; no space before the first and after the last argument. Example: @code for (int i = 0; i < 10; i++) { ... doSomething(arg1, arg2, arg3); } @endcode "else" (and "catch" if it is ever used) is put after the closing brace like this: "} else {" - Use as many brackets in conditions/math terms as you see fit for optimum readability. All operators ('=', '==', '<', '+', '-', '<<', etc.) and castings should always be surrounded by spaces. Examples: @code foo/2 + bar/4 + baz/3 for (int i = 0; i < bar+1; i++) var = (foo - 1) + (bar - 2) + (baz - 3) char *s = (char*) malloc(LENGTH * sizeof(char)); @endcode - enum and structs should be defined first in a class - parameter names in a method definition should only be used to explains the usage of the parameter - In SIGNAL() and SLOT() macros, use as little whitespace as possible. This gives a little speed up since Qt does not have to normalize the signal/slot name. \section constructs Usage of specific constructs * Use C++ casting (static_cast, const_cast, dynamic_cast) in C++ and qobject_cast in Qt classes since they include checks see https://en.wikibooks.org/wiki/C%2B%2B_Programming/Programming_Languages/C%2B%2B/Code/Statements/Variables/Type_Casting * In C++ use Qt container instead of STL container https://marcmutz.wordpress.com/effective-qt/containers/ * In C++ use range-based loops instead of foreach/Q_FOREACH https://www.kdab.com/goodbye-q_foreach/ * For integer data types int is preferred for small numbers and size_t for big, unsigned values. Use double as floating point type. * The 'auto' keyword should be used in range-based loops and for variables initialized by casting or with the 'new' operator but only for non-basic types. Do not omit '*' and '&' to keep readability. * use smart pointers unique_ptr when possible and shared_ptr otherwise. * Avoid const pass-by-value parameters in function declarations. Still make the parameter const in the same function's definition if it won't be modified. * Use the 'override' specifier when overriding virtual functions from the base class http://en.cppreference.com/w/cpp/language/override * Use braces to enclose a single statement only for readability * In C++ nullptr should be used instead of bug-prone NULL and 0. +* Use brace initializing for default values (default initialization should be avoided, like bool, int, double: false, 0, 0.0). + Attention: v{2} initializes a vector to one element of value 2, use v{0., 0.} for two values initalized to 0. * #include <...> vs. #include "...": Include headers from external libraries using angle brackets (as in #include ) and headers from LabPlot/SciDAVis using double quotes (as in #include "core/AbstractAspect.h"). Rationale: Headers of external libraries are never in the same directory as the including file, so it makes sense to use the angle bracket form (which searches only in directories specified using -I). If you work with a build system that does not include the current source directory, or disable CMAKE_INCLUDE_CURRENT_DIR, then all angle-bracket-includes referencing LabPlot/SciDAVis headers will break. Excluding the current directory from -I dirs may be desirable one day when using a new library or a new version of a library containing a header file with the same name as one of our headers. * Use DEBUG() macro for debugging code when possible and QDEBUG() only for special Qt data types. DEBUG() works on all supported systems. @code QString string; DEBUG(" string : " << string.toStdString()); @endcode * Use Qt functions for user messages: qDebug(), qWarning(), qCritical(), qFatal(). Check conditions with Q_ASSERT(cond) or Q_ASSERT_X(cond, where, what) and pointers with Q_CHECK_PTR(ptr). * Import C header (from GSL etc.) with extern statement. We use "std::" prefix (C++11) and try to avoid C header like cmath, cfloat etc. by using corresponding C++ constructs (fabs() -> std::abs(), DBL_MAX -> std::numeric_limits::max(), round() -> qRound()). Example: @code extern "C" { #include } if (std::isnan(x)) { ... } @endcode \section links Links Apart from that, the following links are recommended as guidelines for the coding style: http://techbase.kde.org/index.php?title=Policies/Library_Code_Policy http://doc.trolltech.com/qq/qq13-apis.html http://techbase.kde.org/Policies/Kdelibs_Coding_Style */