labplot  2.8.2
About: LabPlot is an application for plotting and analysis of 2D and 3D functions and data. It is a complete rewrite of LabPlot1 and lacks in the first release a lot of features available in the predecessor. On the other hand, the GUI and the usability is more superior.
  Fossies Dox: labplot-2.8.2.tar.gz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

XYFitCurve.cpp
Go to the documentation of this file.
1 
2 /***************************************************************************
3  File : XYFitCurve.cpp
4  Project : LabPlot
5  Description : A xy-curve defined by a fit model
6  --------------------------------------------------------------------
7  Copyright : (C) 2014-2017 Alexander Semke (alexander.semke@web.de)
8  Copyright : (C) 2016-2020 Stefan Gerlach (stefan.gerlach@uni.kn)
9 
10  ***************************************************************************/
11 
12 /***************************************************************************
13  * *
14  * This program is free software; you can redistribute it and/or modify *
15  * it under the terms of the GNU General Public License as published by *
16  * the Free Software Foundation; either version 2 of the License, or *
17  * (at your option) any later version. *
18  * *
19  * This program is distributed in the hope that it will be useful, *
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
22  * GNU General Public License for more details. *
23  * *
24  * You should have received a copy of the GNU General Public License *
25  * along with this program; if not, write to the Free Software *
26  * Foundation, Inc., 51 Franklin Street, Fifth Floor, *
27  * Boston, MA 02110-1301 USA *
28  * *
29  ***************************************************************************/
30 
31 /*!
32  \class XYFitCurve
33  \brief A xy-curve defined by a fit model
34 
35  \ingroup worksheet
36 */
37 
38 #include "XYFitCurve.h"
39 #include "XYFitCurvePrivate.h"
43 #include "backend/lib/macros.h"
44 #include "backend/gsl/errors.h"
46 
47 extern "C" {
48 #include <gsl/gsl_blas.h>
49 #include <gsl/gsl_vector.h>
50 #include <gsl/gsl_matrix.h>
51 #include <gsl/gsl_version.h>
52 #include <gsl/gsl_cdf.h>
53 #include <gsl/gsl_statistics_double.h>
54 #include "backend/gsl/parser.h"
56 #include "backend/nsl/nsl_stats.h"
57 }
58 
59 #include <QDateTime>
60 #include <QElapsedTimer>
61 #include <QIcon>
62 #include <QThreadPool>
63 
64 XYFitCurve::XYFitCurve(const QString& name)
66 }
67 
68 XYFitCurve::XYFitCurve(const QString& name, XYFitCurvePrivate* dd)
69  : XYAnalysisCurve(name, dd, AspectType::XYFitCurve) {
70 }
71 
72 //no need to delete the d-pointer here - it inherits from QGraphicsItem
73 //and is deleted during the cleanup in QGraphicsScene
74 XYFitCurve::~XYFitCurve() = default;
75 
77  Q_D(XYFitCurve);
78  d->recalculate();
79 }
80 
81 void XYFitCurve::evaluate(bool preview) {
82  Q_D(XYFitCurve);
83  d->evaluate(preview);
84 }
85 
87  Q_D(XYFitCurve);
88  XYFitCurve::FitData& fitData = d->fitData;
89  initStartValues(fitData, curve);
90 }
91 
93  DEBUG("XYFitCurve::initStartValues()");
94  if (!curve) {
95  DEBUG(" no curve given");
96  return;
97  }
98 
99  const Column* tmpXDataColumn = dynamic_cast<const Column*>(curve->xColumn());
100  const Column* tmpYDataColumn = dynamic_cast<const Column*>(curve->yColumn());
101 
102  if (!tmpXDataColumn || !tmpYDataColumn) {
103  DEBUG(Q_FUNC_INFO << ", data columns not available");
104  return;
105  }
106 
107  DEBUG(Q_FUNC_INFO << ", x data rows = " << tmpXDataColumn->rowCount());
108 
109  nsl_fit_model_category modelCategory = fitData.modelCategory;
110  int modelType = fitData.modelType;
111  int degree = fitData.degree;
112  DEBUG(Q_FUNC_INFO << ", fit model type = " << modelType << ", degree = " << degree);
113 
114  QVector<double>& paramStartValues = fitData.paramStartValues;
115  //QVector<double>* xVector = static_cast<QVector<double>* >(tmpXDataColumn->data());
116 
117  //double xmean = gsl_stats_mean(xVector->constData(), 1, tmpXDataColumn->rowCount());
118  double xmin = tmpXDataColumn->minimum();
119  double xmax = tmpXDataColumn->maximum();
120  //double ymin = tmpYDataColumn->minimum();
121  double ymax = tmpYDataColumn->maximum();
122  double xrange = xmax - xmin;
123  //double yrange = ymax-ymin;
124  DEBUG(Q_FUNC_INFO << ", x min/max = " << xmin << ' ' << xmax);
125  //DEBUG(Q_FUNC_INFO <<", y min/max = " << ymin << ' ' << ymax);
126 
127  // guess start values for parameter
128  switch (modelCategory) {
129  case nsl_fit_model_basic:
130  switch (modelType) {
132  // not needed (works anyway)
133  break;
134  //TODO: handle basic models
135  case nsl_fit_model_power:
139  break;
140  }
141  break;
142  case nsl_fit_model_peak:
143  // use equidistant mu's and (xmax-xmin)/(10*degree) as sigma(, gamma)
144  switch (modelType) {
147  case nsl_fit_model_sech:
149  for (int d = 0; d < degree; d++) {
150  paramStartValues[3*d+2] = xmin + (d+1.)*xrange/(degree+1.); // mu
151  paramStartValues[3*d+1] = xrange/(10.*degree); // sigma
152  paramStartValues[3*d] = paramStartValues[3*d+1] * ymax; // A = sigma * ymax
153  }
154  break;
155  case nsl_fit_model_voigt:
156  for (int d = 0; d < degree; d++) {
157  paramStartValues[4*d+1] = xmin + (d+1.)*xrange/(degree+1.); // mu
158  paramStartValues[4*d+2] = xrange/(10.*degree); // sigma
159  paramStartValues[4*d+3] = xrange/(10.*degree); // gamma
160  }
161  break;
163  for (int d = 0; d < degree; d++) {
164  paramStartValues[4*d+1] = 0.5; // eta
165  paramStartValues[4*d+2] = xrange/(10.*degree); // sigma
166  paramStartValues[4*d+3] = xmin + (d+1.)*xrange/(degree+1.); // mu
167  }
168  break;
169  }
170  break;
172  switch (modelType) {
173  case nsl_fit_model_atan:
174  case nsl_fit_model_tanh:
176  case nsl_fit_model_erf:
179  // use (xmax+xmin)/2 as mu and (xmax-xmin)/10 as sigma
180  paramStartValues[1] = (xmax+xmin)/2.;
181  paramStartValues[2] = xrange/10.;
182  break;
183  case nsl_fit_model_hill:
184  paramStartValues[2] = xrange/10.;
185  break;
187  //TODO
188  break;
189  }
190  break;
192  switch (modelType) {
198  case nsl_sf_stats_sech:
200  case nsl_sf_stats_levy:
201  // mu
202  paramStartValues[2] = (xmin+xmax)/2.;
203  // sigma
204  paramStartValues[1] = xrange/10.;
205  // A = sigma * y_max
206  paramStartValues[0] = paramStartValues[1] * ymax;
207  break;
208  //TODO: other types
209  default:
210  break;
211  }
212  break;
214  // not possible
215  break;
216  }
217 }
218 
219 /*!
220  * sets the parameter names for given model category, model type and degree in \c fitData for given action
221  */
223  //TODO: exclude others too?
225  return;
226 
227  Q_D(XYFitCurve);
228  XYFitCurve::FitData& fitData = d->fitData;
230  //Linear
232  fitData.modelType = (int)nsl_fit_model_polynomial;
233  fitData.degree = 1;
235  //Power
237  fitData.modelType = (int)nsl_fit_model_power;
238  fitData.degree = 1;
240  //Exponential (degree 1)
242  fitData.modelType = (int)nsl_fit_model_exponential;
243  fitData.degree = 1;
245  //Exponential (degree 2)
247  fitData.modelType = (int)nsl_fit_model_exponential;
248  fitData.degree = 2;
250  //Inverse exponential
254  //Gauss
256  fitData.modelType = (int)nsl_fit_model_gaussian;
257  fitData.degree = 1;
259  //Cauchy-Lorentz
261  fitData.modelType = (int)nsl_fit_model_lorentz;
262  fitData.degree = 1;
264  //Arc tangent
266  fitData.modelType = (int)nsl_fit_model_atan;
268  //Hyperbolic tangent
270  fitData.modelType = (int)nsl_fit_model_tanh;
272  //Error function
274  fitData.modelType = (int)nsl_fit_model_erf;
275  } else {
276  //Custom
278  fitData.modelType = 0;
279  }
280 
281  XYFitCurve::initFitData(fitData);
282 }
283 
284 /*!
285  * sets the model expression and the parameter names for given model category, model type and degree in \c fitData
286  */
288  nsl_fit_model_category modelCategory = fitData.modelCategory;
289  int modelType = fitData.modelType;
290  QString& model = fitData.model;
291  QStringList& paramNames = fitData.paramNames;
292  QStringList& paramNamesUtf8 = fitData.paramNamesUtf8;
293  int degree = fitData.degree;
294  QVector<double>& paramStartValues = fitData.paramStartValues;
295  QVector<double>& paramLowerLimits = fitData.paramLowerLimits;
296  QVector<double>& paramUpperLimits = fitData.paramUpperLimits;
297  QVector<bool>& paramFixed = fitData.paramFixed;
298 
299  if (modelCategory != nsl_fit_model_custom) {
300  DEBUG(Q_FUNC_INFO << ", XYFitCurve::initFitData() for model category = " << nsl_fit_model_category_name[modelCategory] << ", model type = " << modelType
301  << ", degree = " << degree);
302  paramNames.clear();
303  } else {
304  DEBUG(Q_FUNC_INFO << ", XYFitCurve::initFitData() for model category = nsl_fit_model_custom, model type = " << modelType << ", degree = " << degree);
305  }
306  paramNamesUtf8.clear();
307 
308  // 10 indices used in multi degree models
309  QStringList indices = { UTF8_QSTRING("₁"), UTF8_QSTRING("₂"), UTF8_QSTRING("₃"),
310  UTF8_QSTRING("₄"), UTF8_QSTRING("₅"), UTF8_QSTRING("₆"), UTF8_QSTRING("₇"),
311  UTF8_QSTRING("₈"), UTF8_QSTRING("₉"), UTF8_QSTRING("₁₀")};
312 
313  switch (modelCategory) {
314  case nsl_fit_model_basic:
315  model = nsl_fit_model_basic_equation[fitData.modelType];
316  switch (modelType) {
318  paramNames << "c0" << "c1";
319  paramNamesUtf8 << UTF8_QSTRING("c₀") << UTF8_QSTRING("c₁");
320  if (degree == 2) {
321  model += " + c2*x^2";
322  paramNames << "c2";
323  paramNamesUtf8 << UTF8_QSTRING("c₂");
324  } else if (degree > 2) {
325  for (int i = 2; i <= degree; ++i) {
326  QString numStr = QString::number(i);
327  model += "+c" + numStr + "*x^" + numStr;
328  paramNames << 'c' + numStr;
329  paramNamesUtf8 << 'c' + indices[i-1];
330  }
331  }
332  break;
333  case nsl_fit_model_power:
334  if (degree == 1) {
335  paramNames << "a" << "b";
336  } else {
337  paramNames << "a" << "b" << "c";
338  model = "a + b*x^c";
339  }
340  break;
342  if (degree == 1) {
343  paramNames << "a" << "b";
344  } else {
345  for (int i = 1; i <= degree; i++) {
346  QString numStr = QString::number(i);
347  if (i == 1)
348  model = "a1*exp(b1*x)";
349  else
350  model += " + a" + numStr + "*exp(b" + numStr + "*x)";
351  paramNames << 'a' + numStr << 'b' + numStr;
352  paramNamesUtf8 << 'a' + indices[i-1] << 'b' + indices[i-1];
353  }
354  }
355  break;
357  paramNames << "a" << "b" << "c";
358  break;
360  paramNames << "w" << "a0" << "a1" << "b1";
361  paramNamesUtf8 << UTF8_QSTRING("ω") << UTF8_QSTRING("a₀")
362  << UTF8_QSTRING("a₁") << UTF8_QSTRING("b₁");
363  if (degree > 1) {
364  for (int i = 1; i <= degree; ++i) {
365  QString numStr = QString::number(i);
366  model += "+ (a" + numStr + "*cos(" + numStr + "*w*x) + b" + numStr + "*sin(" + numStr + "*w*x))";
367  paramNames << 'a' + numStr << 'b' + numStr;
368  paramNamesUtf8 << 'a' + indices[i-1] << 'b' + indices[i-1];
369  }
370  }
371  break;
372  }
373  break;
374  case nsl_fit_model_peak:
375  model = nsl_fit_model_peak_equation[fitData.modelType];
376  switch (modelType) {
378  switch (degree) {
379  case 1:
380  paramNames << "a" << "s" << "mu";
381  paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ");
382  break;
383  default:
384  model = "1./sqrt(2*pi) * (";
385  for (int i = 1; i <= degree; ++i) {
386  QString numStr = QString::number(i);
387  if (i > 1)
388  model += " + ";
389  model += 'a' + numStr + "/s" + numStr + "* exp(-((x-mu" + numStr + ")/s" + numStr + ")^2/2)";
390  paramNames << 'a' + numStr << 's' + numStr << "mu" + numStr;
391  paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1];
392  }
393  model += ')';
394  }
395  break;
397  switch (degree) {
398  case 1:
399  paramNames << "a" << "g" << "mu";
400  paramNamesUtf8 << "A" << UTF8_QSTRING("γ") << UTF8_QSTRING("μ");
401  break;
402  default:
403  model = "1./pi * (";
404  for (int i = 1; i <= degree; ++i) {
405  QString numStr = QString::number(i);
406  if (i > 1)
407  model += " + ";
408  model += 'a' + numStr + " * g" + numStr + "/(g" + numStr + "^2+(x-mu" + numStr + ")^2)";
409  paramNames << 'a' + numStr << 'g' + numStr << "mu" + numStr;
410  paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("γ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1];
411  }
412  model += ')';
413  }
414  break;
415  case nsl_fit_model_sech:
416  switch (degree) {
417  case 1:
418  paramNames << "a" << "s" << "mu";
419  paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ");
420  break;
421  default:
422  model = "1/pi * (";
423  for (int i = 1; i <= degree; ++i) {
424  QString numStr = QString::number(i);
425  if (i > 1)
426  model += " + ";
427  model += 'a' + numStr + "/s" + numStr + "* sech((x-mu" + numStr + ")/s" + numStr + ')';
428  paramNames << 'a' + numStr << 's' + numStr << "mu" + numStr;
429  paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1];
430  }
431  model += ')';
432  }
433  break;
435  switch (degree) {
436  case 1:
437  paramNames << "a" << "s" << "mu";
438  paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ");
439  break;
440  default:
441  model = "1/4 * (";
442  for (int i = 1; i <= degree; ++i) {
443  QString numStr = QString::number(i);
444  if (i > 1)
445  model += " + ";
446  model += 'a' + numStr + "/s" + numStr + "* sech((x-mu" + numStr + ")/2/s" + numStr + ")**2";
447  paramNames << 'a' + numStr << 's' + numStr << "mu" + numStr;
448  paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1];
449  }
450  model += ')';
451  }
452  break;
453  case nsl_fit_model_voigt:
454  switch (degree) {
455  case 1:
456  paramNames << "a" << "mu" << "s" << "g";
457  paramNamesUtf8 << "A" << UTF8_QSTRING("μ") << UTF8_QSTRING("σ") << UTF8_QSTRING("γ");
458  break;
459  default:
460  model.clear();
461  for (int i = 1; i <= degree; ++i) {
462  QString numStr = QString::number(i);
463  if (i > 1)
464  model += " + ";
465  model += 'a' + numStr + "*voigt(x-mu" + numStr + ",s" + numStr + ",g" + numStr + ')';
466  paramNames << 'a' + numStr << "mu" + numStr << 's' + numStr << 'g' + numStr;
467  paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1] << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("γ") + indices[i-1];
468  }
469  }
470  break;
472  switch (degree) {
473  case 1:
474  paramNames << "a" << "et" << "w" << "mu"; // eta function exists!
475  paramNamesUtf8 << "A" << UTF8_QSTRING("η") << "w" << UTF8_QSTRING("μ");
476  break;
477  default:
478  model.clear();
479  for (int i = 1; i <= degree; ++i) {
480  QString numStr = QString::number(i);
481  if (i > 1)
482  model += " + ";
483  model += 'a' + numStr + "*pseudovoigt1(x-mu" + numStr + ",eta" + numStr + ",w" + numStr + ')';
484  paramNames << 'a' + numStr << "eta" + numStr << 'w' + numStr << "mu" + numStr;
485  paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("η") + indices[i-1] << 'w' + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1];
486  }
487  }
488  break;
489  }
490  break;
492  model = nsl_fit_model_growth_equation[fitData.modelType];
493  switch (modelType) {
494  case nsl_fit_model_atan:
495  case nsl_fit_model_tanh:
497  case nsl_fit_model_erf:
499  paramNames << "a" << "mu" << "s";
500  paramNamesUtf8 << "A" << UTF8_QSTRING("μ") << UTF8_QSTRING("σ");
501  break;
503  paramNames << "a" << "mu" << "k";
504  paramNamesUtf8 << "A" << UTF8_QSTRING("μ") << "k";
505  break;
506  case nsl_fit_model_hill:
507  paramNames << "a" << "n" << "a";
508  paramNamesUtf8 << "A" << "n" << UTF8_QSTRING("σ");
509  break;
511  paramNames << "a" << "b" << "c";
512  break;
513  }
514  break;
517  switch (modelType) {
523  case nsl_sf_stats_sech:
524  paramNames << "a" << "s" << "mu";
525  paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ");
526  break;
528  paramNames << "A" << "s" << "a" << "mu";
529  paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << "a" << UTF8_QSTRING("μ");
530  break;
532  paramNames << "a" << "l" << "mu";
533  paramNamesUtf8 << "A" << UTF8_QSTRING("λ") << UTF8_QSTRING("μ");
534  break;
536  paramNames << "a" << "s" << "b" << "mu";
537  paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << "b" << UTF8_QSTRING("μ");
538  break;
540  case nsl_sf_stats_levy:
541  paramNames << "a" << "g" << "mu";
542  paramNamesUtf8 << "A" << UTF8_QSTRING("γ") << UTF8_QSTRING("μ");
543  break;
545  paramNames << "a" << "s";
546  paramNamesUtf8 << "A" << UTF8_QSTRING("σ");
547  break;
548  case nsl_sf_stats_landau:
549  paramNames << "a";
550  paramNamesUtf8 << "A";
551  break;
552  case nsl_sf_stats_levy_alpha_stable: // unused distributions
555  break;
556  case nsl_sf_stats_gamma:
557  paramNames << "a" << "k" << "t";
558  paramNamesUtf8 << "A"<< "k" << UTF8_QSTRING("θ");
559  break;
560  case nsl_sf_stats_flat:
561  paramNames << "A" << "b" << "a";
562  break;
564  paramNames << "a" << "n";
565  paramNamesUtf8 << "A" << "n";
566  break;
567  case nsl_sf_stats_fdist:
568  paramNames << "a" << "n1" << "n2";
569  paramNamesUtf8 << "A" << UTF8_QSTRING("ν₁") << UTF8_QSTRING("ν₂");
570  break;
571  case nsl_sf_stats_tdist:
572  paramNames << "a" << "n";
573  paramNamesUtf8 << "A" << UTF8_QSTRING("ν");
574  break;
575  case nsl_sf_stats_beta:
576  case nsl_sf_stats_pareto:
577  paramNames << "A" << "a" << "b";
578  break;
580  paramNames << "a" << "k" << "l" << "mu";
581  paramNamesUtf8 << "A" << "k" << UTF8_QSTRING("λ") << UTF8_QSTRING("μ");
582  break;
584  paramNames << "a" << "s" << "mu" << "b";
585  paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ") << UTF8_QSTRING("β");
586  break;
588  paramNames << "A" << "a" << "b" << "mu";
589  paramNamesUtf8 << "A" << "a" << "b" << UTF8_QSTRING("μ");
590  break;
592  paramNames << "a" << "l";
593  paramNamesUtf8 << "A" << UTF8_QSTRING("λ");
594  break;
597  case nsl_sf_stats_pascal:
598  paramNames << "a" << "p" << "n";
599  paramNamesUtf8 << "A" << "p" << "n";
600  break;
603  paramNames << "a" << "p";
604  paramNamesUtf8 << "A" << "p";
605  break;
607  paramNames << "a" << "n1" << "n2" << "t";
608  paramNamesUtf8 << "A" << UTF8_QSTRING("n₁") << UTF8_QSTRING("n₂") << "t";
609  break;
611  paramNames << "a" << "s";
612  paramNamesUtf8 << "A" << UTF8_QSTRING("σ");
613  break;
615  paramNames << "a" << "g" << "s" << "mu";
616  paramNamesUtf8 << "A" << UTF8_QSTRING("γ") << UTF8_QSTRING("σ") << UTF8_QSTRING("μ");
617  break;
618  }
619  break;
621  break;
622  }
623  DEBUG(Q_FUNC_INFO << ", model: " << STDSTRING(model));
624  DEBUG(Q_FUNC_INFO << ", # params: " << paramNames.size());
625 
626  if (paramNamesUtf8.isEmpty())
627  paramNamesUtf8 << paramNames;
628 
629  //resize the vector for the start values and set the elements to 1.0
630  //in case a custom model is used, do nothing, we take over the previous values
631  if (modelCategory != nsl_fit_model_custom) {
632  const int np = paramNames.size();
633  paramStartValues.resize(np);
634  paramFixed.resize(np);
635  paramLowerLimits.resize(np);
636  paramUpperLimits.resize(np);
637 
638  for (int i = 0; i < np; ++i) {
639  paramStartValues[i] = 1.0;
640  paramFixed[i] = false;
641  paramLowerLimits[i] = -std::numeric_limits<double>::max();
642  paramUpperLimits[i] = std::numeric_limits<double>::max();
643  }
644 
645  // set some model-dependent start values
646  // TODO: see initStartValues()
647  if (modelCategory == nsl_fit_model_distribution) {
648  if (modelType == (int)nsl_sf_stats_flat)
649  paramStartValues[2] = -1.0;
650  else if (modelType == (int)nsl_sf_stats_levy)
651  paramStartValues[2] = 0.0;
652  else if (modelType == (int)nsl_sf_stats_exponential_power || modelType == (int)nsl_sf_stats_weibull || modelType == (int)nsl_sf_stats_gumbel2
653  || modelType == (int)nsl_sf_stats_frechet)
654  paramStartValues[3] = 0.0;
655  else if (modelType == (int)nsl_sf_stats_binomial || modelType == (int)nsl_sf_stats_negative_binomial || modelType == (int)nsl_sf_stats_pascal
656  || modelType == (int)nsl_sf_stats_geometric || modelType == (int)nsl_sf_stats_logarithmic)
657  paramStartValues[1] = 0.5;
658  }
659  }
660 }
661 
662 /*!
663  Returns an icon to be used in the project explorer.
664 */
665 QIcon XYFitCurve::icon() const {
666  return QIcon::fromTheme("labplot-xy-fit-curve");
667 }
668 
669 //##############################################################################
670 //########################## getter methods ##################################
671 //##############################################################################
672 BASIC_SHARED_D_READER_IMPL(XYFitCurve, const AbstractColumn*, xErrorColumn, xErrorColumn)
673 BASIC_SHARED_D_READER_IMPL(XYFitCurve, const AbstractColumn*, yErrorColumn, yErrorColumn)
674 const QString& XYFitCurve::xErrorColumnPath() const { Q_D(const XYFitCurve); return d->xErrorColumnPath; }
675 const QString& XYFitCurve::yErrorColumnPath() const { Q_D(const XYFitCurve); return d->yErrorColumnPath; }
676 
678 
679 const XYFitCurve::FitResult& XYFitCurve::fitResult() const {
680  Q_D(const XYFitCurve);
681  return d->fitResult;
682 }
683 
684 //##############################################################################
685 //################# setter methods and undo commands ##########################
686 //##############################################################################
687 STD_SETTER_CMD_IMPL_S(XYFitCurve, SetXErrorColumn, const AbstractColumn*, xErrorColumn)
688 void XYFitCurve::setXErrorColumn(const AbstractColumn* column) {
689  Q_D(XYFitCurve);
690  if (column != d->xErrorColumn) {
691  exec(new XYFitCurveSetXErrorColumnCmd(d, column, ki18n("%1: assign x-error")));
693  if (column) {
694  connect(column, &AbstractColumn::dataChanged, this, [=](){ handleSourceDataChanged(); });
695  //TODO disconnect on undo
696  }
697  }
698 }
699 
700 STD_SETTER_CMD_IMPL_S(XYFitCurve, SetYErrorColumn, const AbstractColumn*, yErrorColumn)
701 void XYFitCurve::setYErrorColumn(const AbstractColumn* column) {
702  Q_D(XYFitCurve);
703  if (column != d->yErrorColumn) {
704  exec(new XYFitCurveSetYErrorColumnCmd(d, column, ki18n("%1: assign y-error")));
706  if (column) {
707  connect(column, &AbstractColumn::dataChanged, this, [=](){ handleSourceDataChanged(); });
708  //TODO disconnect on undo
709  }
710  }
711 }
712 
713 // do not recalculate (allow preview)
714 //STD_SETTER_CMD_IMPL_F_S(XYFitCurve, SetFitData, XYFitCurve::FitData, fitData, recalculate)
716 void XYFitCurve::setFitData(const XYFitCurve::FitData& fitData) {
717  Q_D(XYFitCurve);
718  exec(new XYFitCurveSetFitDataCmd(d, fitData, ki18n("%1: set fit options and perform the fit")));
719 }
720 
721 //##############################################################################
722 //######################### Private implementation #############################
723 //##############################################################################
725 
726 //no need to delete xColumn and yColumn, they are deleted
727 //when the parent aspect is removed
729 
730 // data structure to pass parameter to fit functions
731 struct data {
732  size_t n; //number of data points
733  double* x; //pointer to the vector with x-data values
734  double* y; //pointer to the vector with y-data values
735  double* weight; //pointer to the vector with weight values
738  int degree;
739  QString* func; // string containing the definition of the model/function
740  QStringList* paramNames;
741  double* paramMin; // lower parameter limits
742  double* paramMax; // upper parameter limits
743  bool* paramFixed; // parameter fixed?
744 };
745 
746 /*!
747  * \param paramValues vector containing current values of the fit parameters
748  * \param params
749  * \param f vector with the weighted residuals weight[i]*(Yi - y[i])
750  */
751 int func_f(const gsl_vector* paramValues, void* params, gsl_vector* f) {
752  //DEBUG("func_f");
753  size_t n = ((struct data*)params)->n;
754  double* x = ((struct data*)params)->x;
755  double* y = ((struct data*)params)->y;
756  double* weight = ((struct data*)params)->weight;
757  nsl_fit_model_category modelCategory = ((struct data*)params)->modelCategory;
758  unsigned int modelType = ((struct data*)params)->modelType;
759  QStringList* paramNames = ((struct data*)params)->paramNames;
760  double *min = ((struct data*)params)->paramMin;
761  double *max = ((struct data*)params)->paramMax;
762 
763  // set current values of the parameters
764  for (int i = 0; i < paramNames->size(); i++) {
765  double v = gsl_vector_get(paramValues, (size_t)i);
766  // bound values if limits are set
767  assign_symbol(qPrintable(paramNames->at(i)), nsl_fit_map_bound(v, min[i], max[i]));
768  QDEBUG("Parameter"<<i<<" (' "<<paramNames->at(i)<<"')"<<'['<<min[i]<<','<<max[i]
769  <<"] free/bound:"<<QString::number(v, 'g', 15)<<' '<<QString::number(nsl_fit_map_bound(v, min[i], max[i]), 'g', 15));
770  }
771 
773  QString func{*(((struct data*)params)->func)};
774  for (size_t i = 0; i < n; i++) {
775  if (std::isnan(x[i]) || std::isnan(y[i]))
776  continue;
777 
778  // checks for allowed values of x for different models
779  // TODO: more to check
781  if (x[i] < 0)
782  x[i] = 0;
783  }
784 
785  assign_symbol("x", x[i]);
786  //DEBUG("evaluate function \"" << func << "\" @ x = " << x[i] << ":");
787  double Yi = parse(qPrintable(func), qPrintable(numberLocale.name()));
788  //DEBUG(" f(x["<< i <<"]) = " << Yi);
789 
790  if (parse_errors() > 0)
791  return GSL_EINVAL;
792 
793  gsl_vector_set(f, i, sqrt(weight[i]) * (Yi - y[i]));
794  }
795 
796  return GSL_SUCCESS;
797 }
798 
799 /*!
800  * calculates the matrix elements of Jacobian matrix
801  * \param paramValues current parameter values
802  * \param params
803  * \param J Jacobian matrix
804  * */
805 int func_df(const gsl_vector* paramValues, void* params, gsl_matrix* J) {
806  //DEBUG("func_df");
807  const size_t n = ((struct data*)params)->n;
808  double* xVector = ((struct data*)params)->x;
809  double* weight = ((struct data*)params)->weight;
810  nsl_fit_model_category modelCategory = ((struct data*)params)->modelCategory;
811  unsigned int modelType = ((struct data*)params)->modelType;
812  unsigned int degree = ((struct data*)params)->degree;
813  QStringList* paramNames = ((struct data*)params)->paramNames;
814  double *min = ((struct data*)params)->paramMin;
815  double *max = ((struct data*)params)->paramMax;
816  bool *fixed = ((struct data*)params)->paramFixed;
817 
818  // calculate the Jacobian matrix:
819  // Jacobian matrix J(i,j) = df_i / dx_j
820  // where f_i = w_i*(Y_i - y_i),
821  // Y_i = model and the x_j are the parameters
822  double x;
823 
824  switch (modelCategory) {
825  case nsl_fit_model_basic:
826  switch (modelType) {
827  case nsl_fit_model_polynomial: // Y(x) = c0 + c1*x + ... + cn*x^n
828  for (size_t i = 0; i < n; i++) {
829  x = xVector[i];
830  for (unsigned int j = 0; j < (unsigned int)paramNames->size(); ++j) {
831  if (fixed[j])
832  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
833  else
834  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_polynomial_param_deriv(x, j, weight[i]));
835  }
836  }
837  break;
838  case nsl_fit_model_power: // Y(x) = a*x^b or Y(x) = a + b*x^c.
839  if (degree == 1) {
840  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
841  const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
842  for (size_t i = 0; i < n; i++) {
843  x = xVector[i];
844 
845  for (int j = 0; j < 2; j++) {
846  if (fixed[j])
847  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
848  else
849  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_power1_param_deriv(j, x, a, b, weight[i]));
850  }
851  }
852  } else if (degree == 2) {
853  const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
854  const double c = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
855  for (size_t i = 0; i < n; i++) {
856  x = xVector[i];
857 
858  for (int j = 0; j < 3; j++) {
859  if (fixed[j])
860  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
861  else
862  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_power2_param_deriv(j, x, b, c, weight[i]));
863  }
864  }
865  }
866  break;
867  case nsl_fit_model_exponential: { // Y(x) = a*exp(b*x) + c*exp(d*x) + ...
868  double *p = new double[2*degree];
869  for (unsigned int i = 0; i < 2*degree; i++)
870  p[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, i), min[i], max[i]);
871  for (size_t i = 0; i < n; i++) {
872  x = xVector[i];
873 
874  for (unsigned int j = 0; j < 2*degree; j++) {
875  if (fixed[j])
876  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
877  else
878  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exponentialn_param_deriv(j, x, p, weight[i]));
879  }
880  }
881  delete[] p;
882 
883  break;
884  }
885  case nsl_fit_model_inverse_exponential: { // Y(x) = a*(1-exp(b*x))+c
886  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
887  const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
888  for (size_t i = 0; i < n; i++) {
889  x = xVector[i];
890 
891  for (unsigned int j = 0; j < 3; j++) {
892  if (fixed[j])
893  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
894  else
895  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_inverse_exponential_param_deriv(j, x, a, b, weight[i]));
896  }
897  }
898  break;
899  }
900  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)
901  //parameters: w, a0, a1, b1, ... an, bn
902  double* a = new double[degree];
903  double* b = new double[degree];
904  double w = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
905  a[0] = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
906  b[0] = 0;
907  for (unsigned int i = 1; i < degree; ++i) {
908  a[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, 2*i), min[2*i], max[2*i]);
909  b[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, 2*i+1), min[2*i+1], max[2*i+1]);
910  }
911  for (size_t i = 0; i < n; i++) {
912  x = xVector[i];
913  double wd = 0; //first derivative with respect to the w parameter
914  for (unsigned int j = 1; j < degree; ++j) {
915  wd += -a[j]*j*x*sin(j*w*x) + b[j]*j*x*cos(j*w*x);
916  }
917 
918  gsl_matrix_set(J, i, 0, weight[i]*wd);
919  gsl_matrix_set(J, i, 1, weight[i]);
920  for (unsigned int j = 1; j <= degree; ++j) {
921  gsl_matrix_set(J, (size_t)i, (size_t)(2*j), nsl_fit_model_fourier_param_deriv(0, j, x, w, weight[i]));
922  gsl_matrix_set(J, (size_t)i, (size_t)(2*j+1), nsl_fit_model_fourier_param_deriv(1, j, x, w, weight[i]));
923  }
924 
925  for (unsigned int j = 0; j <= 2*degree+1; j++)
926  if (fixed[j])
927  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
928  }
929 
930  delete[] a;
931  delete[] b;
932 
933  break;
934  }
935  }
936  break;
937  case nsl_fit_model_peak:
938  switch (modelType) {
941  case nsl_fit_model_sech:
943  for (size_t i = 0; i < n; i++) {
944  x = xVector[i];
945 
946  for (unsigned int j = 0; j < degree; ++j) {
947  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 3*j), min[3*j], max[3*j]);
948  const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 3*j+1), min[3*j+1], max[3*j+1]);
949  const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3*j+2), min[3*j+2], max[3*j+2]);
950 
951  switch (modelType) {
953  gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_gaussian_param_deriv(0, x, a, s, mu, weight[i]));
954  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]));
955  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]));
956  break;
957  case nsl_fit_model_lorentz: // a,s,t
958  gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_lorentz_param_deriv(0, x, a, s, mu, weight[i]));
959  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]));
960  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]));
961  break;
962  case nsl_fit_model_sech:
963  gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_sech_param_deriv(0, x, a, s, mu, weight[i]));
964  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]));
965  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]));
966  break;
968  gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_logistic_param_deriv(0, x, a, s, mu, weight[i]));
969  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]));
970  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]));
971  break;
972  }
973  }
974 
975  for (unsigned int j = 0; j < 3*degree; j++)
976  if (fixed[j])
977  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
978  }
979  break;
980  case nsl_fit_model_voigt:
981  for (size_t i = 0; i < n; i++) {
982  x = xVector[i];
983 
984  for (unsigned int j = 0; j < degree; ++j) {
985  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j), min[4*j], max[4*j]);
986  const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+1), min[4*j+1], max[4*j+1]);
987  const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+2), min[4*j+2], max[4*j+2]);
988  const double g = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+3), min[4*j+3], max[4*j+3]);
989 
990  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]));
991  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]));
992  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]));
993  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]));
994  }
995  for (unsigned int j = 0; j < 4*degree; j++)
996  if (fixed[j])
997  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
998  }
999  break;
1001  for (size_t i = 0; i < n; i++) {
1002  x = xVector[i];
1003 
1004  for (unsigned int j = 0; j < degree; ++j) {
1005  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j), min[4*j], max[4*j]);
1006  const double eta = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+1), min[4*j+1], max[4*j+1]);
1007  const double w = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+2), min[4*j+2], max[4*j+2]);
1008  const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+3), min[4*j+3], max[4*j+3]);
1009 
1010  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]));
1011  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]));
1012  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]));
1013  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]));
1014  }
1015  for (unsigned int j = 0; j < 4*degree; j++)
1016  if (fixed[j])
1017  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1018  }
1019  break;
1020  }
1021  break;
1022  case nsl_fit_model_growth: {
1023  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1024  const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1025  const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1026 
1027  for (size_t i = 0; i < n; i++) {
1028  x = xVector[i];
1029 
1030  for (unsigned int j = 0; j < 3; j++) {
1031  if (fixed[j]) {
1032  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1033  } else {
1034  switch (modelType) {
1035  case nsl_fit_model_atan:
1036  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_atan_param_deriv(j, x, a, mu, s, weight[i]));
1037  break;
1038  case nsl_fit_model_tanh:
1039  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_tanh_param_deriv(j, x, a, mu, s, weight[i]));
1040  break;
1042  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_algebraic_sigmoid_param_deriv(j, x, a, mu, s, weight[i]));
1043  break;
1044  case nsl_fit_model_sigmoid:
1045  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_sigmoid_param_deriv(j, x, a, mu, s, weight[i]));
1046  break;
1047  case nsl_fit_model_erf:
1048  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_erf_param_deriv(j, x, a, mu, s, weight[i]));
1049  break;
1050  case nsl_fit_model_hill:
1051  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_hill_param_deriv(j, x, a, mu, s, weight[i]));
1052  break;
1054  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gompertz_param_deriv(j, x, a, mu, s, weight[i]));
1055  break;
1057  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gudermann_param_deriv(j, x, a, mu, s, weight[i]));
1058  break;
1059  }
1060  }
1061  }
1062  }
1063  }
1064  break;
1066  switch (modelType) {
1067  case nsl_sf_stats_gaussian:
1069  case nsl_sf_stats_laplace:
1073  case nsl_sf_stats_logistic:
1074  case nsl_sf_stats_sech:
1075  case nsl_sf_stats_levy: {
1076  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1077  const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1078  const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1079  for (size_t i = 0; i < n; i++) {
1080  x = xVector[i];
1081 
1082  for (unsigned int j = 0; j < 3; j++) {
1083  if (fixed[j])
1084  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1085  else {
1086  switch (modelType) {
1087  case nsl_sf_stats_gaussian:
1088  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gaussian_param_deriv(j, x, a, s, mu, weight[i]));
1089  break;
1091  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exponential_param_deriv(j, x, a, s, mu, weight[i]));
1092  break;
1093  case nsl_sf_stats_laplace:
1094  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_laplace_param_deriv(j, x, a, s, mu, weight[i]));
1095  break;
1097  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_lorentz_param_deriv(j, x, a, s, mu, weight[i]));
1098  break;
1100  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_rayleigh_tail_param_deriv(j, x, a, s, mu, weight[i]));
1101  break;
1103  if (x > 0)
1104  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_lognormal_param_deriv(j, x, a, s, mu, weight[i]));
1105  else
1106  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1107  break;
1108  case nsl_sf_stats_logistic:
1109  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_logistic_param_deriv(j, x, a, s, mu, weight[i]));
1110  break;
1111  case nsl_sf_stats_sech:
1112  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_sech_dist_param_deriv(j, x, a, s, mu, weight[i]));
1113  break;
1114  case nsl_sf_stats_levy:
1115  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_levy_param_deriv(j, x, a, s, mu, weight[i]));
1116  break;
1117  }
1118  }
1119  }
1120  }
1121  break;
1122  }
1124  const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1125  const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1126  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1127  const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1128  for (size_t i = 0; i < n; i++) {
1129  x = xVector[i];
1130 
1131  for (unsigned int j = 0; j < 4; j++) {
1132  if (fixed[j])
1133  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1134  else
1135  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]));
1136  }
1137  }
1138  break;
1139  }
1141  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1142  const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1143  const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1144  const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1145  for (size_t i = 0; i < n; i++) {
1146  x = xVector[i];
1147 
1148  for (unsigned int j = 0; j < 4; j++) {
1149  if (fixed[j])
1150  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1151  else
1152  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]));
1153  }
1154  }
1155  break;
1156  }
1157  case nsl_sf_stats_rayleigh: {
1158  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1159  const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1160  for (size_t i = 0; i < n; i++) {
1161  x = xVector[i];
1162 
1163  for (unsigned int j = 0; j < 2; j++) {
1164  if (fixed[j])
1165  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1166  else
1167  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_rayleigh_param_deriv(j, x, a, s, weight[i]));
1168  }
1169  }
1170  break;
1171  }
1172  case nsl_sf_stats_gamma: {
1173  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1174  const double k = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1175  const double t = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1176  for (size_t i = 0; i < n; i++) {
1177  x = xVector[i];
1178 
1179  for (unsigned int j = 0; j < 3; j++) {
1180  if (fixed[j])
1181  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1182  else
1183  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gamma_param_deriv(j, x, a, k, t, weight[i]));
1184  }
1185  }
1186  break;
1187  }
1188  case nsl_sf_stats_flat: {
1189  const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1190  const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1191  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1192  for (size_t i = 0; i < n; i++) {
1193  x = xVector[i];
1194 
1195  for (unsigned int j = 0; j < 3; j++) {
1196  if (fixed[j])
1197  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1198  else
1199  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_flat_param_deriv(j, x, A, b, a, weight[i]));
1200  }
1201  }
1202  break;
1203  }
1204  case nsl_sf_stats_chi_squared: {
1205  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1206  const double nu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1207  for (size_t i = 0; i < n; i++) {
1208  x = xVector[i];
1209 
1210  for (unsigned int j = 0; j < 2; j++) {
1211  if (fixed[j])
1212  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1213  else
1214  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_chi_square_param_deriv(j, x, a, nu, weight[i]));
1215  }
1216  }
1217  break;
1218  }
1219  case nsl_sf_stats_tdist: {
1220  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1221  const double nu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1222  for (size_t i = 0; i < n; i++) {
1223  x = xVector[i];
1224 
1225  for (unsigned int j = 0; j < 2; j++) {
1226  if (fixed[j])
1227  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1228  else
1229  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_students_t_param_deriv(j, x, a, nu, weight[i]));
1230  }
1231  }
1232  break;
1233  }
1234  case nsl_sf_stats_fdist: {
1235  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1236  const double n1 = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1237  const double n2 = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1238  for (size_t i = 0; i < n; i++) {
1239  x = xVector[i];
1240 
1241  for (unsigned int j = 0; j < 3; j++) {
1242  if (fixed[j])
1243  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1244  else
1245  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_fdist_param_deriv(j, x, a, n1, n2, weight[i]));
1246  }
1247  }
1248  break;
1249  }
1250  case nsl_sf_stats_beta:
1251  case nsl_sf_stats_pareto: {
1252  const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1253  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1254  const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1255  for (size_t i = 0; i < n; i++) {
1256  x = xVector[i];
1257 
1258  for (unsigned int j = 0; j < 3; j++) {
1259  if (fixed[j])
1260  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1261  else {
1262  switch (modelType) {
1263  case nsl_sf_stats_beta:
1264  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_beta_param_deriv(j, x, A, a, b, weight[i]));
1265  break;
1266  case nsl_sf_stats_pareto:
1267  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_pareto_param_deriv(j, x, A, a, b, weight[i]));
1268  break;
1269  }
1270  }
1271  }
1272  }
1273  break;
1274  }
1275  case nsl_sf_stats_weibull: {
1276  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1277  const double k = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1278  const double l = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1279  const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1280  for (size_t i = 0; i < n; i++) {
1281  x = xVector[i];
1282 
1283  for (unsigned int j = 0; j < 4; j++) {
1284  if (fixed[j])
1285  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1286  else {
1287  if (x > 0)
1288  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_weibull_param_deriv(j, x, a, k, l, mu, weight[i]));
1289  else
1290  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1291  }
1292  }
1293  }
1294  break;
1295  }
1296  case nsl_sf_stats_gumbel1: {
1297  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1298  const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1299  const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1300  const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1301  for (size_t i = 0; i < n; i++) {
1302  x = xVector[i];
1303 
1304  for (unsigned int j = 0; j < 4; j++) {
1305  if (fixed[j])
1306  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1307  else
1308  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gumbel1_param_deriv(j, x, a, s, mu, b, weight[i]));
1309  }
1310  }
1311  break;
1312  }
1313  case nsl_sf_stats_gumbel2: {
1314  const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1315  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1316  const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1317  const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1318  for (size_t i = 0; i < n; i++) {
1319  x = xVector[i];
1320 
1321  for (unsigned int j = 0; j < 4; j++) {
1322  if (fixed[j])
1323  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1324  else
1325  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gumbel2_param_deriv(j, x, A, a, b, mu, weight[i]));
1326  }
1327  }
1328  break;
1329  }
1330  case nsl_sf_stats_poisson: {
1331  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1332  const double l = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1333  for (size_t i = 0; i < n; i++) {
1334  x = xVector[i];
1335 
1336  for (unsigned int j = 0; j < 2; j++) {
1337  if (fixed[j])
1338  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1339  else
1340  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_poisson_param_deriv(j, x, a, l, weight[i]));
1341  }
1342  }
1343  break;
1344  }
1345  case nsl_sf_stats_maxwell_boltzmann: { // Y(x) = a*sqrt(2/pi) * x^2/s^3 * exp(-(x/s)^2/2)
1346  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1347  const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1348  for (size_t i = 0; i < n; i++) {
1349  x = xVector[i];
1350 
1351  for (unsigned int j = 0; j < 2; j++) {
1352  if (fixed[j])
1353  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1354  else
1355  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_maxwell_param_deriv(j, x, a, s, weight[i]));
1356  }
1357  }
1358  break;
1359  }
1360  case nsl_sf_stats_frechet: {
1361  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1362  const double g = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1363  const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1364  const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1365  for (size_t i = 0; i < n; i++) {
1366  x = xVector[i];
1367 
1368  for (unsigned int j = 0; j < 4; j++) {
1369  if (fixed[j])
1370  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1371  else
1372  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_frechet_param_deriv(j, x, a, g, s, mu, weight[i]));
1373  }
1374  }
1375  break;
1376  }
1377  case nsl_sf_stats_landau: {
1378  // const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1379  for (size_t i = 0; i < n; i++) {
1380  x = xVector[i];
1381  if (fixed[0])
1382  gsl_matrix_set(J, (size_t)i, 0, 0.);
1383  else
1384  gsl_matrix_set(J, (size_t)i, 0, nsl_fit_model_landau_param_deriv(0, x, weight[i]));
1385  }
1386  break;
1387  }
1388  case nsl_sf_stats_binomial:
1390  case nsl_sf_stats_pascal: {
1391  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1392  const double p = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1393  const double N = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1394  for (size_t i = 0; i < n; i++) {
1395  x = xVector[i];
1396 
1397  for (unsigned int j = 0; j < 3; j++) {
1398  if (fixed[j])
1399  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1400  else {
1401  switch (modelType) {
1402  case nsl_sf_stats_binomial:
1403  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_binomial_param_deriv(j, x, a, p, N, weight[i]));
1404  break;
1406  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_negative_binomial_param_deriv(j, x, a, p, N, weight[i]));
1407  break;
1408  case nsl_sf_stats_pascal:
1409  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_pascal_param_deriv(j, x, a, p, N, weight[i]));
1410  break;
1411  }
1412  }
1413  }
1414  }
1415  break;
1416  }
1418  case nsl_sf_stats_logarithmic: {
1419  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1420  const double p = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1421  for (size_t i = 0; i < n; i++) {
1422  x = xVector[i];
1423 
1424  for (unsigned int j = 0; j < 2; j++) {
1425  if (fixed[j])
1426  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1427  else {
1428  switch (modelType) {
1430  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_geometric_param_deriv(j, x, a, p, weight[i]));
1431  break;
1433  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_logarithmic_param_deriv(j, x, a, p, weight[i]));
1434  break;
1435  }
1436  }
1437  }
1438  }
1439  break;
1440  }
1442  const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]);
1443  const double n1 = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]);
1444  const double n2 = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]);
1445  const double t = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]);
1446  for (size_t i = 0; i < n; i++) {
1447  x = xVector[i];
1448 
1449  for (unsigned int j = 0; j < 4; j++) {
1450  if (fixed[j])
1451  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1452  else
1453  gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_hypergeometric_param_deriv(j, x, a, n1, n2, t, weight[i]));
1454  }
1455  }
1456  break;
1457  }
1458  // unused distributions
1462  break;
1463  }
1464  break;
1465  case nsl_fit_model_custom:
1466  double value;
1467  const unsigned int np = paramNames->size();
1468  QString func{*(((struct data*)params)->func)};
1469 
1471  for (size_t i = 0; i < n; i++) {
1472  x = xVector[i];
1473  assign_symbol("x", x);
1474 
1475  for (unsigned int j = 0; j < np; j++) {
1476  for (unsigned int k = 0; k < np; k++) {
1477  if (k != j) {
1478  value = nsl_fit_map_bound(gsl_vector_get(paramValues, k), min[k], max[k]);
1479  assign_symbol(qPrintable(paramNames->at(k)), value);
1480  }
1481  }
1482 
1483  value = nsl_fit_map_bound(gsl_vector_get(paramValues, j), min[j], max[j]);
1484  assign_symbol(qPrintable(paramNames->at(j)), value);
1485  const double f_p = parse(qPrintable(func), qPrintable(numberLocale.name()));
1486 
1487  double eps = 1.e-9;
1488  if (std::abs(f_p) > 0)
1489  eps *= std::abs(f_p); // scale step size with function value
1490  value += eps;
1491  assign_symbol(qPrintable(paramNames->at(j)), value);
1492  const double f_pdp = parse(qPrintable(func), qPrintable(numberLocale.name()));
1493 
1494 // DEBUG("evaluate deriv"<<func<<": f(x["<<i<<"]) ="<<QString::number(f_p, 'g', 15));
1495 // DEBUG("evaluate deriv"<<func<<": f(x["<<i<<"]+dx) ="<<QString::number(f_pdp, 'g', 15));
1496 // DEBUG(" deriv = " << STDSTRING(QString::number(sqrt(weight[i])*(f_pdp-f_p)/eps, 'g', 15));
1497 
1498  if (fixed[j])
1499  gsl_matrix_set(J, (size_t)i, (size_t)j, 0.);
1500  else // calculate finite difference
1501  gsl_matrix_set(J, (size_t)i, (size_t)j, sqrt(weight[i])*(f_pdp - f_p)/eps);
1502  }
1503  }
1504  }
1505 
1506  return GSL_SUCCESS;
1507 }
1508 
1509 int func_fdf(const gsl_vector* x, void* params, gsl_vector* f, gsl_matrix* J) {
1510  //DEBUG("func_fdf");
1511  func_f(x, params, f);
1512  func_df(x, params, J);
1513 
1514  return GSL_SUCCESS;
1515 }
1516 
1517 /* prepare the fit result columns */
1519  DEBUG("XYFitCurvePrivate::prepareResultColumns()")
1520  //create fit result columns if not available yet, clear them otherwise
1521  if (!xColumn) { // all columns are treated together
1522  DEBUG(" Creating columns")
1526  xVector = static_cast<QVector<double>* >(xColumn->data());
1527  yVector = static_cast<QVector<double>* >(yColumn->data());
1528  residualsVector = static_cast<QVector<double>* >(residualsColumn->data());
1529 
1530  xColumn->setHidden(true);
1531  q->addChild(xColumn);
1532 
1533  yColumn->setHidden(true);
1534  q->addChild(yColumn);
1535 
1537 
1538  q->setUndoAware(false);
1539  q->setXColumn(xColumn);
1540  q->setYColumn(yColumn);
1541  q->setUndoAware(true);
1542  } else {
1543  DEBUG(Q_FUNC_INFO << ", Clear columns")
1544  xVector->clear();
1545  yVector->clear();
1546  //TODO: residualsVector->clear();
1547  }
1548 }
1549 
1551  QElapsedTimer timer;
1552  timer.start();
1553 
1554  // prepare source data columns
1555  const AbstractColumn* tmpXDataColumn = nullptr;
1556  const AbstractColumn* tmpYDataColumn = nullptr;
1558  DEBUG(Q_FUNC_INFO << ", spreadsheet columns as data source");
1559  tmpXDataColumn = xDataColumn;
1560  tmpYDataColumn = yDataColumn;
1561  } else {
1562  DEBUG(Q_FUNC_INFO << ", curve columns as data source");
1563  tmpXDataColumn = dataSourceCurve->xColumn();
1564  tmpYDataColumn = dataSourceCurve->yColumn();
1565  }
1566 
1567  // clear the previous result
1569 
1570  if (!tmpXDataColumn || !tmpYDataColumn) {
1571  DEBUG(Q_FUNC_INFO << ", ERROR: Preparing source data columns failed!");
1572  emit q->dataChanged();
1574  return;
1575  }
1576 
1578 
1579  //fit settings
1580  const unsigned int maxIters = fitData.maxIterations; //maximal number of iterations
1581  const double delta = fitData.eps; //fit tolerance
1582  const unsigned int np = fitData.paramNames.size(); //number of fit parameters
1583  if (np == 0) {
1584  fitResult.available = true;
1585  fitResult.valid = false;
1586  fitResult.status = i18n("Model has no parameters.");
1587  emit q->dataChanged();
1589  return;
1590  }
1591 
1592  if (yErrorColumn) {
1593  if (yErrorColumn->rowCount() < tmpXDataColumn->rowCount()) {
1594  fitResult.available = true;
1595  fitResult.valid = false;
1596  fitResult.status = i18n("Not sufficient weight data points provided.");
1597  emit q->dataChanged();
1599  return;
1600  }
1601  }
1602 
1603  //copy all valid data point for the fit to temporary vectors
1604  QVector<double> xdataVector;
1605  QVector<double> ydataVector;
1606  QVector<double> xerrorVector;
1607  QVector<double> yerrorVector;
1608  Range<double> xRange{tmpXDataColumn->minimum(), tmpXDataColumn->maximum()};
1609  if (fitData.autoRange) { // auto x range of data to fit
1610  fitData.fitRange = xRange;
1611  } else { // custom x range of data to fit
1612  if (!fitData.fitRange.isZero()) // avoid problems with user specified zero range
1613  xRange.setRange(fitData.fitRange.min(), fitData.fitRange.max());
1614  }
1615  DEBUG(Q_FUNC_INFO << ", fit range = " << xRange.min() << " .. " << xRange.max());
1616  DEBUG(Q_FUNC_INFO << ", fitData range = " << fitData.fitRange.min() << " .. " << fitData.fitRange.max());
1617 
1618  //logic from XYAnalysisCurve::copyData(), extended by the handling of error columns.
1619  //TODO: decide how to deal with non-numerical error columns
1620  int rowCount = qMin(tmpXDataColumn->rowCount(), tmpYDataColumn->rowCount());
1621  for (int row = 0; row < rowCount; ++row) {
1622  // omit invalid data
1623  if (!tmpXDataColumn->isValid(row) || tmpXDataColumn->isMasked(row) ||
1624  !tmpYDataColumn->isValid(row) || tmpYDataColumn->isMasked(row))
1625  continue;
1626 
1627  double x = NAN;
1628  switch (tmpXDataColumn->columnMode()) {
1630  x = tmpXDataColumn->valueAt(row);
1631  break;
1633  x = tmpXDataColumn->integerAt(row);
1634  break;
1636  x = tmpXDataColumn->bigIntAt(row);
1637  break;
1638  case AbstractColumn::ColumnMode::Text: // not valid
1639  break;
1643  x = tmpXDataColumn->dateTimeAt(row).toMSecsSinceEpoch();
1644  }
1645 
1646  double y = NAN;
1647  switch (tmpYDataColumn->columnMode()) {
1649  y = tmpYDataColumn->valueAt(row);
1650  break;
1652  y = tmpYDataColumn->integerAt(row);
1653  break;
1655  y = tmpYDataColumn->bigIntAt(row);
1656  break;
1657  case AbstractColumn::ColumnMode::Text: // not valid
1658  break;
1662  y = tmpYDataColumn->dateTimeAt(row).toMSecsSinceEpoch();
1663  }
1664 
1665  if (x >= xRange.min() && x <= xRange.max()) { // only when inside given range
1666  if ((!xErrorColumn && !yErrorColumn) || !fitData.useDataErrors) { // x-y
1667  xdataVector.append(x);
1668  ydataVector.append(y);
1669  } else if (!xErrorColumn && yErrorColumn) { // x-y-dy
1670  if (!std::isnan(yErrorColumn->valueAt(row))) {
1671  xdataVector.append(x);
1672  ydataVector.append(y);
1673  yerrorVector.append(yErrorColumn->valueAt(row));
1674  }
1675  } else if (xErrorColumn && yErrorColumn) { // x-y-dx-dy
1676  if (!std::isnan(xErrorColumn->valueAt(row)) && !std::isnan(yErrorColumn->valueAt(row))) {
1677  xdataVector.append(x);
1678  ydataVector.append(y);
1679  xerrorVector.append(xErrorColumn->valueAt(row));
1680  yerrorVector.append(yErrorColumn->valueAt(row));
1681  }
1682  }
1683  }
1684  }
1685 
1686  //number of data points to fit
1687  const size_t n = xdataVector.size();
1688  DEBUG(Q_FUNC_INFO << ", number of data points: " << n);
1689  if (n == 0) {
1690  fitResult.available = true;
1691  fitResult.valid = false;
1692  fitResult.status = i18n("No data points available.");
1693  emit q->dataChanged();
1695  return;
1696  }
1697 
1698  if (n < np) {
1699  fitResult.available = true;
1700  fitResult.valid = false;
1701  fitResult.status = i18n("The number of data points (%1) must be greater than or equal to the number of parameters (%2).", n, np);
1702  emit q->dataChanged();
1704  return;
1705  }
1706 
1707  if (fitData.model.simplified().isEmpty()) {
1708  fitResult.available = true;
1709  fitResult.valid = false;
1710  fitResult.status = i18n("Fit model not specified.");
1711  emit q->dataChanged();
1713  return;
1714  }
1715 
1716  double* xdata = xdataVector.data();
1717  double* ydata = ydataVector.data();
1718  double* xerror = xerrorVector.data(); // size may be 0
1719  double* yerror = yerrorVector.data(); // size may be 0
1720  DEBUG(Q_FUNC_INFO << ", x error vector size: " << xerrorVector.size());
1721  DEBUG(Q_FUNC_INFO << ", y error vector size: " << yerrorVector.size());
1722  double* weight = new double[n];
1723 
1724  for (size_t i = 0; i < n; i++)
1725  weight[i] = 1.;
1726 
1727  const double minError = 1.e-199; // minimum error for weighting
1728 
1729  switch (fitData.yWeightsType) {
1730  case nsl_fit_weight_no:
1733  break;
1734  case nsl_fit_weight_instrumental: // yerror are sigmas
1735  for (int i = 0; i < (int)n; i++)
1736  if (i < yerrorVector.size())
1737  weight[i] = 1./gsl_pow_2(qMax(yerror[i], qMax(sqrt(minError), fabs(ydata[i]) * 1.e-15)));
1738  break;
1739  case nsl_fit_weight_direct: // yerror are weights
1740  for (int i = 0; i < (int)n; i++)
1741  if (i < yerrorVector.size())
1742  weight[i] = yerror[i];
1743  break;
1744  case nsl_fit_weight_inverse: // yerror are inverse weights
1745  for (int i = 0; i < (int)n; i++)
1746  if (i < yerrorVector.size())
1747  weight[i] = 1./qMax(yerror[i], qMax(minError, fabs(ydata[i]) * 1.e-15));
1748  break;
1750  for (int i = 0; i < (int)n; i++)
1751  weight[i] = 1./qMax(ydata[i], minError);
1752  break;
1754  for (int i = 0; i < (int)n; i++)
1755  weight[i] = 1./qMax(gsl_pow_2(ydata[i]), minError);
1756  break;
1757  }
1758 
1759  /////////////////////// GSL >= 2 has a complete new interface! But the old one is still supported. ///////////////////////////
1760  // 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"
1761  unsigned int nf = 0; // number of fixed parameter
1762  for (unsigned int i = 0; i < np; i++) {
1763  const bool fixed = fitData.paramFixed.data()[i];
1764  if (fixed)
1765  nf++;
1766  DEBUG(" parameter " << i << " fixed: " << fixed);
1767  }
1768 
1769  //function to fit
1770  gsl_multifit_function_fdf f;
1771  DEBUG(Q_FUNC_INFO << ", model = " << STDSTRING(fitData.model));
1773  f.f = &func_f;
1774  f.df = &func_df;
1775  f.fdf = &func_fdf;
1776  f.n = n;
1777  f.p = np;
1778  f.params = &params;
1779 
1780  DEBUG(Q_FUNC_INFO << ", initialize the derivative solver (using Levenberg-Marquardt robust solver)");
1781  const gsl_multifit_fdfsolver_type* T = gsl_multifit_fdfsolver_lmsder;
1782  gsl_multifit_fdfsolver* s = gsl_multifit_fdfsolver_alloc(T, n, np);
1783 
1784  DEBUG(Q_FUNC_INFO << ", set start values");
1785  double* x_init = fitData.paramStartValues.data();
1786  double* x_min = fitData.paramLowerLimits.data();
1787  double* x_max = fitData.paramUpperLimits.data();
1788  DEBUG(Q_FUNC_INFO << ", scale start values if limits are set");
1789  for (unsigned int i = 0; i < np; i++)
1790  x_init[i] = nsl_fit_map_unbound(x_init[i], x_min[i], x_max[i]);
1791  DEBUG(Q_FUNC_INFO << ", DONE");
1792  gsl_vector_view x = gsl_vector_view_array(x_init, np);
1793  DEBUG(Q_FUNC_INFO << ", Turning off GSL error handler to avoid overflow/underflow");
1794  gsl_set_error_handler_off();
1795  DEBUG(Q_FUNC_INFO << ", Initialize solver with function f and initial guess x");
1796  gsl_multifit_fdfsolver_set(s, &f, &x.vector);
1797 
1798  DEBUG(Q_FUNC_INFO << ", Iterate ...");
1799  int status = GSL_SUCCESS;
1800  unsigned int iter = 0;
1801  fitResult.solverOutput.clear();
1802  writeSolverState(s);
1803  do {
1804  iter++;
1805  DEBUG(Q_FUNC_INFO << ", iter " << iter);
1806 
1807  // update weights for Y-depending weights (using function values from residuals)
1809  for (size_t i = 0; i < n; i++)
1810  weight[i] = 1./(gsl_vector_get(s->f, i)/sqrt(weight[i]) + ydata[i]); // 1/Y_i
1812  for (size_t i = 0; i < n; i++)
1813  weight[i] = 1./gsl_pow_2(gsl_vector_get(s->f, i)/sqrt(weight[i]) + ydata[i]); // 1/Y_i^2
1814  }
1815 
1816  if (nf == np) { // all fixed parameter
1817  DEBUG(Q_FUNC_INFO << ", all parameter fixed. Stop iteration.")
1818  break;
1819  }
1820  DEBUG(Q_FUNC_INFO << ", run fdfsolver_iterate");
1821  status = gsl_multifit_fdfsolver_iterate(s);
1822  DEBUG(Q_FUNC_INFO << ", fdfsolver_iterate DONE");
1823  double chi = gsl_blas_dnrm2(s->f);
1824  writeSolverState(s, chi);
1825  if (status) {
1826  DEBUG(Q_FUNC_INFO << ", iter " << iter << ", status = " << gsl_strerror(status));
1827  if (status == GSL_ETOLX) // change in the position vector falls below machine precision: no progress
1828  status = GSL_SUCCESS;
1829  break;
1830  }
1831  if (qFuzzyIsNull(chi)) {
1832  DEBUG(Q_FUNC_INFO << ", chi is zero! Finishing.")
1833  status = GSL_SUCCESS;
1834  } else {
1835  status = gsl_multifit_test_delta(s->dx, s->x, delta, delta);
1836  }
1837  DEBUG(Q_FUNC_INFO << ", iter " << iter << ", test status = " << gsl_strerror(status));
1838  } while (status == GSL_CONTINUE && iter < maxIters);
1839 
1840  // second run for x-error fitting
1841  if (xerrorVector.size() > 0) {
1842  DEBUG(Q_FUNC_INFO << ", Rerun fit with x errors");
1843 
1844  unsigned int iter2 = 0;
1845  double chi = 0, chiOld = 0;
1846  double *fun = new double[n];
1847  do {
1848  iter2++;
1849  chiOld = chi;
1850  //printf("iter2 = %d\n", iter2);
1851 
1852  // calculate function from residuals
1853  for (size_t i = 0; i < n; i++)
1854  fun[i] = gsl_vector_get(s->f, i) * 1./sqrt(weight[i]) + ydata[i];
1855 
1856  // calculate weight[i]
1857  for (size_t i = 0; i < n; i++) {
1858  // calculate df[i]
1859  size_t index = i-1;
1860  if (i == 0)
1861  index = i;
1862  if (i == n-1)
1863  index = i-2;
1864  double df = (fun[index+1] - fun[index])/(xdata[index+1] - xdata[index]);
1865  //printf("df = %g\n", df);
1866 
1867  double sigmasq = 1.;
1868  switch (fitData.xWeightsType) { // x-error type: f'(x)^2*s_x^2 = f'(x)/w_x
1869  case nsl_fit_weight_no:
1870  break;
1871  case nsl_fit_weight_direct: // xerror = w_x
1872  sigmasq = df*df/qMax(xerror[i], minError);
1873  break;
1874  case nsl_fit_weight_instrumental: // xerror = s_x
1875  sigmasq = df*df*xerror[i]*xerror[i];
1876  break;
1877  case nsl_fit_weight_inverse: // xerror = 1/w_x = s_x^2
1878  sigmasq = df*df*xerror[i];
1879  break;
1880  case nsl_fit_weight_statistical: // s_x^2 = 1/w_x = x
1881  sigmasq = xdata[i];
1882  break;
1883  case nsl_fit_weight_relative: // s_x^2 = 1/w_x = x^2
1884  sigmasq = xdata[i]*xdata[i];
1885  break;
1886  case nsl_fit_weight_statistical_fit: // unused
1888  break;
1889  }
1890 
1891  if (yerrorVector.size() > 0) {
1892  switch (fitData.yWeightsType) { // y-error types: s_y^2 = 1/w_y
1893  case nsl_fit_weight_no:
1894  break;
1895  case nsl_fit_weight_direct: // yerror = w_y
1896  sigmasq += 1./qMax(yerror[i], minError);
1897  break;
1898  case nsl_fit_weight_instrumental: // yerror = s_y
1899  sigmasq += yerror[i]*yerror[i];
1900  break;
1901  case nsl_fit_weight_inverse: // yerror = 1/w_y
1902  sigmasq += yerror[i];
1903  break;
1904  case nsl_fit_weight_statistical: // unused
1906  break;
1907  case nsl_fit_weight_statistical_fit: // s_y^2 = 1/w_y = Y_i
1908  sigmasq += fun[i];
1909  break;
1910  case nsl_fit_weight_relative_fit: // s_y^2 = 1/w_y = Y_i^2
1911  sigmasq += fun[i]*fun[i];
1912  break;
1913  }
1914  }
1915 
1916  //printf ("sigma[%d] = %g\n", i, sqrt(sigmasq));
1917  weight[i] = 1./qMax(sigmasq, minError);
1918  }
1919 
1920  // update weights
1921  gsl_multifit_fdfsolver_set(s, &f, &x.vector);
1922 
1923  do { // fit
1924  iter++;
1925  writeSolverState(s);
1926  status = gsl_multifit_fdfsolver_iterate(s);
1927  //printf ("status = %s\n", gsl_strerror (status));
1928  if (nf == np) // stop if all parameters fix
1929  break;
1930 
1931  if (status) {
1932  DEBUG(" iter " << iter << ", status = " << gsl_strerror(status));
1933  if (status == GSL_ETOLX) // change in the position vector falls below machine precision: no progress
1934  status = GSL_SUCCESS;
1935  break;
1936  }
1937  status = gsl_multifit_test_delta(s->dx, s->x, delta, delta);
1938  } while (status == GSL_CONTINUE && iter < maxIters);
1939 
1940  chi = gsl_blas_dnrm2(s->f);
1941  } while (iter2 < maxIters && fabs(chi - chiOld) > fitData.eps);
1942 
1943  delete[] fun;
1944  }
1945 
1946  delete[] weight;
1947 
1948  // unscale start parameter
1949  for (unsigned int i = 0; i < np; i++)
1950  x_init[i] = nsl_fit_map_bound(x_init[i], x_min[i], x_max[i]);
1951 
1952  //get the covariance matrix
1953  //TODO: scale the Jacobian when limits are used before constructing the covar matrix?
1954  gsl_matrix* covar = gsl_matrix_alloc(np, np);
1955 #if GSL_MAJOR_VERSION >= 2
1956  // the Jacobian is not part of the solver anymore
1957  gsl_matrix *J = gsl_matrix_alloc(s->fdf->n, s->fdf->p);
1958  gsl_multifit_fdfsolver_jac(s, J);
1959  gsl_multifit_covar(J, 0.0, covar);
1960 #else
1961  gsl_multifit_covar(s->J, 0.0, covar);
1962 #endif
1963 
1964  //write the result
1965  fitResult.available = true;
1966  fitResult.valid = true;
1967  fitResult.status = gslErrorToString(status);
1968  fitResult.iterations = iter;
1969  fitResult.dof = n - (np - nf); // samples - (parameter - fixed parameter)
1970 
1971  //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])
1972  //gsl_blas_dasum() - computes the absolute sum \sum |r_i| of the elements of the vector with the elements weight[i]*(Yi - y[i])
1973  fitResult.sse = gsl_pow_2(gsl_blas_dnrm2(s->f));
1974 
1975  if (fitResult.dof != 0) {
1977  fitResult.rsd = sqrt(fitResult.rms);
1978  }
1980  fitResult.rmse = sqrt(fitResult.mse);
1981  fitResult.mae = gsl_blas_dasum(s->f)/n;
1982  // SST needed for coefficient of determination, R-squared and F test
1983  fitResult.sst = gsl_stats_tss(ydata, 1, n);
1984  // for a linear model without intercept R-squared is calculated differently
1985  // 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
1987  DEBUG(" Using alternative R^2 for linear model without intercept");
1988  fitResult.sst = gsl_stats_tss_m(ydata, 1, n, 0);
1989  }
1990  if (fitResult.sst < fitResult.sse) {
1991  DEBUG(" Using alternative R^2 since R^2 would be negative (probably custom model without intercept)");
1992  fitResult.sst = gsl_stats_tss_m(ydata, 1, n, 0);
1993  }
1994 
2003 
2004  //parameter values
2005  fitResult.paramValues.resize(np);
2006  fitResult.errorValues.resize(np);
2007  fitResult.tdist_tValues.resize(np);
2008  fitResult.tdist_pValues.resize(np);
2009  fitResult.tdist_marginValues.resize(np);
2010  // GSL: cerr = GSL_MAX_DBL(1., sqrt(fitResult.rms)); // increase error for poor fit
2011  // NIST: cerr = sqrt(fitResult.rms); // increase error for poor fit, decrease for good fit
2012  const double cerr = sqrt(fitResult.rms);
2013  // CI = 100* (1 - alpha)
2014  const double alpha = 1.0 - fitData.confidenceInterval/100.;
2015  for (unsigned int i = 0; i < np; i++) {
2016  // scale resulting values if they are bounded
2017  fitResult.paramValues[i] = nsl_fit_map_bound(gsl_vector_get(s->x, i), x_min[i], x_max[i]);
2018  // use results as start values if desired
2019  if (fitData.useResults) {
2021  DEBUG(" saving parameter " << i << ": " << fitResult.paramValues[i] << ' ' << fitData.paramStartValues.data()[i]);
2022  }
2023  fitResult.errorValues[i] = cerr * sqrt(gsl_matrix_get(covar, i, i));
2027  for (unsigned int j = 0; j <= i; j++)
2028  fitResult.correlationMatrix << gsl_matrix_get(covar, i, j)/sqrt(gsl_matrix_get(covar, i, i))/sqrt(gsl_matrix_get(covar, j, j));
2029  }
2030 
2031  // fill residuals vector. To get residuals on the correct x values, fill the rest with zeros.
2032  residualsVector->resize(tmpXDataColumn->rowCount());
2033  DEBUG(" Residual vector size: " << residualsVector->size())
2034  if (fitData.autoRange) { // evaluate full range of residuals
2035  xVector->resize(tmpXDataColumn->rowCount());
2036  auto mode = tmpXDataColumn->columnMode();
2037  for (int i = 0; i < tmpXDataColumn->rowCount(); i++)
2039  (*xVector)[i] = tmpXDataColumn->valueAt(i);
2041  (*xVector)[i] = tmpXDataColumn->integerAt(i);
2043  (*xVector)[i] = tmpXDataColumn->bigIntAt(i);
2045  (*xVector)[i] = tmpXDataColumn->dateTimeAt(i).toMSecsSinceEpoch();
2046 
2048  bool rc = parser->evaluateCartesian(fitData.model, xVector, residualsVector,
2050  if (rc) {
2051  for (int i = 0; i < tmpXDataColumn->rowCount(); i++)
2052  (*residualsVector)[i] = tmpYDataColumn->valueAt(i) - (*residualsVector)[i];
2053  } else {
2054  DEBUG(" ERROR: Failed parsing residuals")
2055  residualsVector->clear();
2056  }
2057  } else { // only selected range
2058  size_t j = 0;
2059  for (int i = 0; i < tmpXDataColumn->rowCount(); i++) {
2060  if (tmpXDataColumn->valueAt(i) >= xRange.min() && tmpXDataColumn->valueAt(i) <= xRange.max())
2061  residualsVector->data()[i] = - gsl_vector_get(s->f, j++);
2062  else // outside range
2063  residualsVector->data()[i] = 0;
2064  }
2065  }
2067 
2068  //free resources
2069  gsl_multifit_fdfsolver_free(s);
2070  gsl_matrix_free(covar);
2071 
2072  //calculate the fit function (vectors)
2073  evaluate();
2074  fitResult.elapsedTime = timer.elapsed();
2075 
2077 }
2078 
2079 /* evaluate fit function (preview == true: use start values, default: false) */
2080 void XYFitCurvePrivate::evaluate(bool preview) {
2081  DEBUG(Q_FUNC_INFO << ", preview = " << preview);
2082 
2083  // prepare source data columns
2084  const AbstractColumn* tmpXDataColumn = nullptr;
2086  DEBUG(Q_FUNC_INFO << ", spreadsheet columns as data source");
2087  tmpXDataColumn = xDataColumn;
2088  } else {
2089  DEBUG(Q_FUNC_INFO << ", curve columns as data source");
2090  if (dataSourceCurve)
2091  tmpXDataColumn = dataSourceCurve->xColumn();
2092  }
2093 
2094  if (!tmpXDataColumn) {
2095  DEBUG(" ERROR: Preparing source data column failed!");
2097  emit q->dataChanged();
2098  return;
2099  }
2100 
2101  //only needed for preview (else we have all columns)
2102  // should not harm even if not in preview now that residuals are not cleared
2103  if (preview)
2105 
2106  if (!xVector || !yVector) {
2107  DEBUG(Q_FUNC_INFO << ", xVector or yVector not defined!");
2109  emit q->dataChanged();
2110  return;
2111  }
2112 
2113  if (fitData.model.simplified().isEmpty()) {
2114  DEBUG(Q_FUNC_INFO << ", no fit-model specified.");
2116  emit q->dataChanged();
2117  return;
2118  }
2119 
2121  Range<double> xRange{tmpXDataColumn->minimum(), tmpXDataColumn->maximum()}; // full data range
2122  if (!fitData.autoEvalRange) { // use given range for evaluation
2123  if (!fitData.evalRange.isZero()) // avoid zero range
2124  xRange = fitData.evalRange;
2125  }
2126  DEBUG(Q_FUNC_INFO << ", eval range = " << STDSTRING(xRange.toString()));
2127  xVector->resize((int)fitData.evaluatedPoints);
2128  yVector->resize((int)fitData.evaluatedPoints);
2129  DEBUG(Q_FUNC_INFO << ", vector size = " << xVector->size());
2130 
2131  QVector<double> paramValues = fitResult.paramValues;
2132  if (preview) // results not available yet
2133  paramValues = fitData.paramStartValues;
2134 
2135  bool rc = parser->evaluateCartesian(fitData.model, QString::number(xRange.min()), QString::number(xRange.max()), (int)fitData.evaluatedPoints,
2136  xVector, yVector, fitData.paramNames, paramValues);
2137  if (!rc) {
2138  DEBUG(Q_FUNC_INFO << ", ERROR: Parsing fit function failed")
2139  xVector->clear();
2140  yVector->clear();
2141  residualsVector->clear();
2142  }
2143 
2145  emit q->dataChanged();
2146 }
2147 
2148 /*!
2149  * writes out the current state of the solver \c s
2150  */
2151 void XYFitCurvePrivate::writeSolverState(gsl_multifit_fdfsolver* s, double chi) {
2152  QString state;
2153 
2154  //current parameter values, semicolon separated
2155  double* min = fitData.paramLowerLimits.data();
2156  double* max = fitData.paramUpperLimits.data();
2157  for (int i = 0; i < fitData.paramNames.size(); ++i) {
2158  const double x = gsl_vector_get(s->x, i);
2159  // map parameter if bounded
2160  state += QString::number(nsl_fit_map_bound(x, min[i], max[i])) + '\t';
2161  }
2162 
2163  //current value of chi
2164  if (isnan(chi))
2165  chi = gsl_blas_dnrm2(s->f);
2166  state += QString::number(chi*chi);
2167  state += ';';
2168  DEBUG(Q_FUNC_INFO << ", chi^2 = " << chi*chi);
2169 
2170  fitResult.solverOutput += state;
2171 }
2172 
2173 
2174 //##############################################################################
2175 //################## Serialization/Deserialization ###########################
2176 //##############################################################################
2177 //! Save as XML
2178 void XYFitCurve::save(QXmlStreamWriter* writer) const {
2179  Q_D(const XYFitCurve);
2180 
2181  writer->writeStartElement("xyFitCurve");
2182 
2183  //write the base class
2184  XYAnalysisCurve::save(writer);
2185 
2186  //write xy-fit-curve specific information
2187 
2188  //fit data - only save model expression and parameter names for custom model, otherwise they are set in XYFitCurve::initFitData()
2189  writer->writeStartElement("fitData");
2190  WRITE_COLUMN(d->xErrorColumn, xErrorColumn);
2191  WRITE_COLUMN(d->yErrorColumn, yErrorColumn);
2192  writer->writeAttribute("autoRange", QString::number(d->fitData.autoRange));
2193  writer->writeAttribute("fitRangeMin", QString::number(d->fitData.fitRange.min(), 'g', 15));
2194  writer->writeAttribute("fitRangeMax", QString::number(d->fitData.fitRange.max(), 'g', 15));
2195  writer->writeAttribute("modelCategory", QString::number(d->fitData.modelCategory));
2196  writer->writeAttribute("modelType", QString::number(d->fitData.modelType));
2197  writer->writeAttribute("xWeightsType", QString::number(d->fitData.xWeightsType));
2198  writer->writeAttribute("weightsType", QString::number(d->fitData.yWeightsType));
2199  writer->writeAttribute("degree", QString::number(d->fitData.degree));
2200  if (d->fitData.modelCategory == nsl_fit_model_custom)
2201  writer->writeAttribute("model", d->fitData.model);
2202  writer->writeAttribute("maxIterations", QString::number(d->fitData.maxIterations));
2203  writer->writeAttribute("eps", QString::number(d->fitData.eps, 'g', 15));
2204  writer->writeAttribute("evaluatedPoints", QString::number(d->fitData.evaluatedPoints));
2205  writer->writeAttribute("autoEvalRange", QString::number(d->fitData.autoEvalRange));
2206  writer->writeAttribute("evalRangeMin", QString::number(d->fitData.evalRange.min(), 'g', 15));
2207  writer->writeAttribute("evalRangeMax", QString::number(d->fitData.evalRange.max(), 'g', 15));
2208  writer->writeAttribute("useDataErrors", QString::number(d->fitData.useDataErrors));
2209  writer->writeAttribute("useResults", QString::number(d->fitData.useResults));
2210  writer->writeAttribute("previewEnabled", QString::number(d->fitData.previewEnabled));
2211  writer->writeAttribute("confidenceInterval", QString::number(d->fitData.confidenceInterval));
2212 
2213  if (d->fitData.modelCategory == nsl_fit_model_custom) {
2214  writer->writeStartElement("paramNames");
2215  foreach (const QString &name, d->fitData.paramNames)
2216  writer->writeTextElement("name", name);
2217  writer->writeEndElement();
2218  }
2219 
2220  writer->writeStartElement("paramStartValues");
2221  foreach (const double &value, d->fitData.paramStartValues)
2222  writer->writeTextElement("startValue", QString::number(value, 'g', 15));
2223  writer->writeEndElement();
2224 
2225  // use 16 digits to handle -DBL_MAX
2226  writer->writeStartElement("paramLowerLimits");
2227  foreach (const double &limit, d->fitData.paramLowerLimits)
2228  writer->writeTextElement("lowerLimit", QString::number(limit, 'g', 16));
2229  writer->writeEndElement();
2230 
2231  // use 16 digits to handle DBL_MAX
2232  writer->writeStartElement("paramUpperLimits");
2233  foreach (const double &limit, d->fitData.paramUpperLimits)
2234  writer->writeTextElement("upperLimit", QString::number(limit, 'g', 16));
2235  writer->writeEndElement();
2236 
2237  writer->writeStartElement("paramFixed");
2238  foreach (const double &fixed, d->fitData.paramFixed)
2239  writer->writeTextElement("fixed", QString::number(fixed));
2240  writer->writeEndElement();
2241 
2242  writer->writeEndElement(); //"fitData"
2243 
2244  //fit results (generated columns and goodness of the fit)
2245  writer->writeStartElement("fitResult");
2246  writer->writeAttribute("available", QString::number(d->fitResult.available));
2247  writer->writeAttribute("valid", QString::number(d->fitResult.valid));
2248  writer->writeAttribute("status", d->fitResult.status);
2249  writer->writeAttribute("iterations", QString::number(d->fitResult.iterations));
2250  writer->writeAttribute("time", QString::number(d->fitResult.elapsedTime));
2251  writer->writeAttribute("dof", QString::number(d->fitResult.dof));
2252  writer->writeAttribute("sse", QString::number(d->fitResult.sse, 'g', 15));
2253  writer->writeAttribute("sst", QString::number(d->fitResult.sst, 'g', 15));
2254  writer->writeAttribute("rms", QString::number(d->fitResult.rms, 'g', 15));
2255  writer->writeAttribute("rsd", QString::number(d->fitResult.rsd, 'g', 15));
2256  writer->writeAttribute("mse", QString::number(d->fitResult.mse, 'g', 15));
2257  writer->writeAttribute("rmse", QString::number(d->fitResult.rmse, 'g', 15));
2258  writer->writeAttribute("mae", QString::number(d->fitResult.mae, 'g', 15));
2259  writer->writeAttribute("rsquare", QString::number(d->fitResult.rsquare, 'g', 15));
2260  writer->writeAttribute("rsquareAdj", QString::number(d->fitResult.rsquareAdj, 'g', 15));
2261  writer->writeAttribute("chisq_p", QString::number(d->fitResult.chisq_p, 'g', 15));
2262  writer->writeAttribute("fdist_F", QString::number(d->fitResult.fdist_F, 'g', 15));
2263  writer->writeAttribute("fdist_p", QString::number(d->fitResult.fdist_p, 'g', 15));
2264  writer->writeAttribute("aic", QString::number(d->fitResult.aic, 'g', 15));
2265  writer->writeAttribute("bic", QString::number(d->fitResult.bic, 'g', 15));
2266  writer->writeAttribute("solverOutput", d->fitResult.solverOutput);
2267 
2268  writer->writeStartElement("paramValues");
2269  foreach (const double &value, d->fitResult.paramValues)
2270  writer->writeTextElement("value", QString::number(value, 'g', 15));
2271  writer->writeEndElement();
2272 
2273  writer->writeStartElement("errorValues");
2274  foreach (const double &value, d->fitResult.errorValues)
2275  writer->writeTextElement("error", QString::number(value, 'g', 15));
2276  writer->writeEndElement();
2277 
2278  writer->writeStartElement("tdist_tValues");
2279  foreach (const double &value, d->fitResult.tdist_tValues)
2280  writer->writeTextElement("tdist_t", QString::number(value, 'g', 15));
2281  writer->writeEndElement();
2282 
2283  writer->writeStartElement("tdist_pValues");
2284  foreach (const double &value, d->fitResult.tdist_pValues)
2285  writer->writeTextElement("tdist_p", QString::number(value, 'g', 15));
2286  writer->writeEndElement();
2287 
2288  writer->writeStartElement("tdist_marginValues");
2289  foreach (const double &value, d->fitResult.tdist_marginValues)
2290  writer->writeTextElement("tdist_margin", QString::number(value, 'g', 15));
2291  writer->writeEndElement();
2292 
2293  writer->writeStartElement("correlationMatrix");
2294  foreach (const double &value, d->fitResult.correlationMatrix)
2295  writer->writeTextElement("correlation", QString::number(value, 'g', 15));
2296  writer->writeEndElement();
2297 
2298  //save calculated columns if available
2299  if (d->xColumn && d->yColumn && d->residualsColumn) {
2300  d->xColumn->save(writer);
2301  d->yColumn->save(writer);
2302  d->residualsColumn->save(writer);
2303  }
2304 
2305  writer->writeEndElement(); //"fitResult"
2306  writer->writeEndElement(); //"xyFitCurve"
2307 }
2308 
2309 //! Load from XML
2310 bool XYFitCurve::load(XmlStreamReader* reader, bool preview) {
2311  Q_D(XYFitCurve);
2312 
2313  KLocalizedString attributeWarning = ki18n("Attribute '%1' missing or empty, default value is used");
2314  QXmlStreamAttributes attribs;
2315  QString str;
2316 
2317  while (!reader->atEnd()) {
2318  reader->readNext();
2319  if (reader->isEndElement() && reader->name() == "xyFitCurve")
2320  break;
2321 
2322  if (!reader->isStartElement())
2323  continue;
2324 
2325  if (reader->name() == "xyAnalysisCurve") {
2326  if ( !XYAnalysisCurve::load(reader, preview) )
2327  return false;
2328  } else if (!preview && reader->name() == "fitData") {
2329  attribs = reader->attributes();
2330 
2331  READ_COLUMN(xErrorColumn);
2332  READ_COLUMN(yErrorColumn);
2333 
2334  READ_INT_VALUE("autoRange", fitData.autoRange, bool);
2335  READ_DOUBLE_VALUE("xRangeMin", fitData.fitRange.min()); // old name
2336  READ_DOUBLE_VALUE("xRangeMax", fitData.fitRange.max()); // old name
2337  READ_DOUBLE_VALUE("fitRangeMin", fitData.fitRange.min());
2338  READ_DOUBLE_VALUE("fitRangeMax", fitData.fitRange.max());
2339  READ_INT_VALUE("modelCategory", fitData.modelCategory, nsl_fit_model_category);
2340  READ_INT_VALUE("modelType", fitData.modelType, int);
2341  READ_INT_VALUE("xWeightsType", fitData.xWeightsType, nsl_fit_weight_type);
2342  READ_INT_VALUE("weightsType", fitData.yWeightsType, nsl_fit_weight_type);
2343  READ_INT_VALUE("degree", fitData.degree, int);
2344  if (d->fitData.modelCategory == nsl_fit_model_custom) {
2345  READ_STRING_VALUE("model", fitData.model);
2346  DEBUG("read model = " << STDSTRING(d->fitData.model));
2347  }
2348  READ_INT_VALUE("maxIterations", fitData.maxIterations, int);
2349  READ_DOUBLE_VALUE("eps", fitData.eps);
2350  READ_INT_VALUE("fittedPoints", fitData.evaluatedPoints, size_t); // old name
2351  READ_INT_VALUE("evaluatedPoints", fitData.evaluatedPoints, size_t);
2352  READ_INT_VALUE("evaluateFullRange", fitData.autoEvalRange, bool); // old name
2353  READ_INT_VALUE("autoEvalRange", fitData.autoEvalRange, bool);
2354  READ_DOUBLE_VALUE("evalRangeMin", fitData.evalRange.min());
2355  READ_DOUBLE_VALUE("evalRangeMax", fitData.evalRange.max());
2356  READ_INT_VALUE("useDataErrors", fitData.useDataErrors, bool);
2357  READ_INT_VALUE("useResults", fitData.useResults, bool);
2358  READ_INT_VALUE("previewEnabled", fitData.previewEnabled, bool);
2359  READ_DOUBLE_VALUE("confidenceInterval", fitData.confidenceInterval);
2360 
2361  //set the model expression and the parameter names (can be derived from the saved values for category, type and degree)
2362  XYFitCurve::initFitData(d->fitData);
2363  // remove default names and start values (will be read from project later)
2364  d->fitData.paramStartValues.clear();
2365 
2366  } else if (!preview && reader->name() == "name") { // needed for custom model
2367  d->fitData.paramNames << reader->readElementText();
2368  } else if (!preview && reader->name() == "startValue") {
2369  d->fitData.paramStartValues << reader->readElementText().toDouble();
2370  } else if (!preview && reader->name() == "fixed") {
2371  d->fitData.paramFixed << (bool)reader->readElementText().toInt();
2372  } else if (!preview && reader->name() == "lowerLimit") {
2373  bool ok;
2374  double x = reader->readElementText().toDouble(&ok);
2375  if (ok) // -DBL_MAX results in conversion error
2376  d->fitData.paramLowerLimits << x;
2377  else
2378  d->fitData.paramLowerLimits << -std::numeric_limits<double>::max();
2379  } else if (!preview && reader->name() == "upperLimit") {
2380  bool ok;
2381  double x = reader->readElementText().toDouble(&ok);
2382  if (ok) // DBL_MAX results in conversion error
2383  d->fitData.paramUpperLimits << x;
2384  else
2385  d->fitData.paramUpperLimits << std::numeric_limits<double>::max();
2386  } else if (!preview && reader->name() == "value") {
2387  d->fitResult.paramValues << reader->readElementText().toDouble();
2388  } else if (!preview && reader->name() == "error") {
2389  d->fitResult.errorValues << reader->readElementText().toDouble();
2390  } else if (!preview && reader->name() == "tdist_t") {
2391  d->fitResult.tdist_tValues << reader->readElementText().toDouble();
2392  } else if (!preview && reader->name() == "tdist_p") {
2393  d->fitResult.tdist_pValues << reader->readElementText().toDouble();
2394  } else if (!preview && reader->name() == "tdist_margin") {
2395  d->fitResult.tdist_marginValues << reader->readElementText().toDouble();
2396  } else if (!preview && reader->name() == "correlation") {
2397  d->fitResult.correlationMatrix << reader->readElementText().toDouble();
2398  } else if (!preview && reader->name() == "fitResult") {
2399  attribs = reader->attributes();
2400 
2401  READ_INT_VALUE("available", fitResult.available, int);
2402  READ_INT_VALUE("valid", fitResult.valid, int);
2403  READ_STRING_VALUE("status", fitResult.status);
2404  READ_INT_VALUE("iterations", fitResult.iterations, int);
2405  READ_INT_VALUE("time", fitResult.elapsedTime, int);
2412  READ_DOUBLE_VALUE("rmse", fitResult.rmse);
2414  READ_DOUBLE_VALUE("rsquare", fitResult.rsquare);
2415  READ_DOUBLE_VALUE("rsquareAdj", fitResult.rsquareAdj);
2416  READ_DOUBLE_VALUE("chisq_p", fitResult.chisq_p);
2417  READ_DOUBLE_VALUE("fdist_F", fitResult.fdist_F);
2418  READ_DOUBLE_VALUE("fdist_p", fitResult.fdist_p);
2421  READ_STRING_VALUE("solverOutput", fitResult.solverOutput);
2422  } else if (reader->name() == "column") {
2423  Column* column = new Column(QString(), AbstractColumn::ColumnMode::Numeric);
2424  if (!column->load(reader, preview)) {
2425  delete column;
2426  return false;
2427  }
2428  DEBUG("############################ reading column " << STDSTRING(column->name()))
2429  if (column->name() == "x")
2430  d->xColumn = column;
2431  else if (column->name() == "y")
2432  d->yColumn = column;
2433  else if (column->name() == "residuals")
2434  d->residualsColumn = column;
2435  }
2436  }
2437 
2438  ////////////////////////////// fix old projects /////////////////////////
2439 
2440  // reset model type of old projects due to new model style
2441  if (d->fitData.modelCategory == nsl_fit_model_basic && d->fitData.modelType >= NSL_FIT_MODEL_BASIC_COUNT) {
2442  DEBUG(Q_FUNC_INFO << ", RESET old fit model");
2443  d->fitData.modelType = 0;
2444  d->fitData.degree = 1;
2445  d->fitData.paramNames.clear();
2446  d->fitData.paramNamesUtf8.clear();
2447  // reset size of fields not touched by initFitData()
2448  d->fitData.paramStartValues.resize(2);
2449  d->fitData.paramFixed.resize(2);
2450  d->fitResult.paramValues.resize(2);
2451  d->fitResult.errorValues.resize(2);
2452  d->fitResult.tdist_tValues.resize(2);
2453  d->fitResult.tdist_pValues.resize(2);
2454  d->fitResult.tdist_marginValues.resize(2);
2455  }
2456 
2457  // older projects also save the param names for non-custom models: remove them
2458  while (d->fitData.paramNames.size() > d->fitData.paramStartValues.size())
2459  d->fitData.paramNames.removeLast();
2460 
2461  // not present in old projects
2462  if (d->fitData.paramNamesUtf8.isEmpty())
2463  d->fitData.paramNamesUtf8 << d->fitData.paramNames;
2464 
2465  // not present in old projects
2466  int np = d->fitResult.paramValues.size();
2467  if (d->fitResult.tdist_tValues.size() == 0)
2468  d->fitResult.tdist_tValues.resize(np);
2469  if (d->fitResult.tdist_pValues.size() == 0)
2470  d->fitResult.tdist_pValues.resize(np);
2471  if (d->fitResult.tdist_marginValues.size() == 0)
2472  d->fitResult.tdist_marginValues.resize(np);
2473  if (d->fitResult.correlationMatrix.size() == 0)
2474  d->fitResult.correlationMatrix.resize(np*(np+1)/2);
2475 
2476  // Loading done. Check some parameter
2477  DEBUG(Q_FUNC_INFO << ", model type = " << d->fitData.modelType);
2478  DEBUG(Q_FUNC_INFO << ", # params = " << d->fitData.paramNames.size());
2479  DEBUG(Q_FUNC_INFO << ", # start values = " << d->fitData.paramStartValues.size());
2480  //for (const auto& value : d->fitData.paramStartValues)
2481  // DEBUG("XYFitCurve::load() # start value = " << value);
2482 
2483  if (preview)
2484  return true;
2485 
2486  // wait for data to be read before using the pointers
2487  QThreadPool::globalInstance()->waitForDone();
2488 
2489  if (d->xColumn && d->yColumn && d->residualsColumn) {
2490  d->xColumn->setHidden(true);
2491  addChild(d->xColumn);
2492 
2493  d->yColumn->setHidden(true);
2494  addChild(d->yColumn);
2495 
2496  addChild(d->residualsColumn);
2497 
2498  d->xVector = static_cast<QVector<double>* >(d->xColumn->data());
2499  d->yVector = static_cast<QVector<double>* >(d->yColumn->data());
2500  d->residualsVector = static_cast<QVector<double>* >(d->residualsColumn->data());
2501 
2502  XYCurve::d_ptr->xColumn = d->xColumn;
2503  XYCurve::d_ptr->yColumn = d->yColumn;
2504 
2506  }
2507 
2508  return true;
2509 }
AspectType
int func_fdf(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J)
int func_f(const gsl_vector *paramValues, void *params, gsl_vector *f)
Definition: XYFitCurve.cpp:751
int func_df(const gsl_vector *paramValues, void *params, gsl_matrix *J)
Definition: XYFitCurve.cpp:805
void setUndoAware(bool)
void addChild(AbstractAspect *)
Add the given Aspect to my list of children.
QString name() const
void exec(QUndoCommand *)
Execute the given command, pushing it on the undoStack() if available.
AbstractAspectPrivate * d
void setHidden(bool)
Set "hidden" property, i.e. whether to exclude this aspect from being shown in the explorer.
Interface definition for data with column logic.
virtual double valueAt(int row) const
Return the double value in row 'row'.
bool isMasked(int row) const
Return whether a certain row is masked.
virtual int rowCount() const =0
Return the data vector size.
virtual double maximum(int count=0) const
virtual double minimum(int count=0) const
virtual QDateTime dateTimeAt(int row) const
Return the QDateTime in row 'row'.
void dataChanged(const AbstractColumn *source)
Data of the column has changed.
virtual qint64 bigIntAt(int row) const
Return the bigint value in row 'row'.
bool isValid(int row) const
Convenience method for mode-independent testing of validity.
virtual int integerAt(int row) const
Return the integer value in row 'row'.
virtual ColumnMode columnMode() const =0
Return the column mode.
Aspect that manages a column.
Definition: Column.h:42
double minimum(int count=0) const override
Definition: Column.cpp:1507
int rowCount() const override
Return the data vector size.
Definition: Column.cpp:1395
void setChanged()
Definition: Column.cpp:998
double maximum(int count=0) const override
Definition: Column.cpp:1653
void * data() const
Definition: Column.cpp:852
bool load(XmlStreamReader *, bool preview) override
Load the column from XML.
Definition: Column.cpp:1152
static ExpressionParser * getInstance()
Auxiliary class for a data range.
Definition: Range.h:47
XYAnalysisCurve::DataSourceType dataSourceType
const AbstractColumn * yDataColumn
QVector< double > * xVector
const AbstractColumn * xDataColumn
QVector< double > * yVector
Base class for all analysis curves.
void save(QXmlStreamWriter *) const override
Save as XML.
bool load(XmlStreamReader *, bool preview) override
Load from XML.
void handleSourceDataChanged()
bool sourceDataChangedSinceLastRecalc
const AbstractColumn * xColumn
void recalcLogicalPoints()
Definition: XYCurve.cpp:1090
const AbstractColumn * yColumn
A 2D-curve, provides an interface for editing many properties of the curve.
Definition: XYCurve.h:46
void dataChanged()
XYCurvePrivate *const d_ptr
Definition: XYCurve.h:186
void recalcLogicalPoints()
Definition: XYCurve.cpp:765
void writeSolverState(gsl_multifit_fdfsolver *, double chi=NAN)
void prepareResultColumns()
~XYFitCurvePrivate() override
const AbstractColumn * xErrorColumn
XYFitCurve::FitResult fitResult
XYFitCurve::FitData fitData
void evaluate(bool preview=false)
const AbstractColumn * yErrorColumn
XYFitCurvePrivate(XYFitCurve *)
Definition: XYFitCurve.cpp:724
QVector< double > * residualsVector
XYFitCurve *const q
A xy-curve defined by a fit model.
Definition: XYFitCurve.h:43
void recalculate() override
Definition: XYFitCurve.cpp:76
~XYFitCurve() override
void initFitData(PlotDataDialog::AnalysisAction)
Definition: XYFitCurve.cpp:222
const FitResult & fitResult() const
Definition: XYFitCurve.cpp:679
void save(QXmlStreamWriter *) const override
Save as XML.
void initStartValues(const XYCurve *)
Definition: XYFitCurve.cpp:86
QIcon icon() const override
Definition: XYFitCurve.cpp:665
const QString & yErrorColumnPath() const
Definition: XYFitCurve.cpp:675
XYFitCurve(const QString &name)
Definition: XYFitCurve.cpp:64
void evaluate(bool preview)
Definition: XYFitCurve.cpp:81
bool load(XmlStreamReader *, bool preview) override
Load from XML.
XML stream parser that supports errors as well as warnings. This class also adds line and column numb...
#define READ_DOUBLE_VALUE(name, var)
Definition: macros.h:475
#define STDSTRING(qstr)
Definition: macros.h:67
#define DEBUG(x)
Definition: macros.h:50
#define SET_NUMBER_LOCALE
Definition: macros.h:75
#define STD_SETTER_CMD_IMPL_S(class_name, cmd_name, value_type, field_name)
Definition: macros.h:207
#define UTF8_QSTRING(str)
Definition: macros.h:66
#define READ_COLUMN(columnName)
Definition: macros.h:462
#define READ_INT_VALUE(name, var, type)
Definition: macros.h:468
#define QDEBUG(x)
Definition: macros.h:47
#define BASIC_SHARED_D_READER_IMPL(classname, type, method, var)
Definition: macros.h:127
#define READ_STRING_VALUE(name, var)
Definition: macros.h:482
#define WRITE_COLUMN(column, columnName)
Definition: macros.h:451
#define i18n(m)
Definition: nsl_common.h:38
const char * nsl_fit_model_growth_equation[]
Definition: nsl_fit.c:52
double nsl_fit_model_gaussian_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight)
Definition: nsl_fit.c:166
double nsl_fit_model_sech_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight)
Definition: nsl_fit.c:190
double nsl_fit_model_gamma_param_deriv(unsigned int param, double x, double A, double k, double t, double weight)
Definition: nsl_fit.c:431
double nsl_fit_model_weibull_param_deriv(unsigned int param, double x, double A, double k, double l, double mu, double weight)
Definition: nsl_fit.c:554
double nsl_fit_model_binomial_param_deriv(unsigned int param, double k, double A, double p, double n, double weight)
Definition: nsl_fit.c:610
double nsl_fit_model_rayleigh_param_deriv(unsigned int param, double x, double A, double s, double weight)
Definition: nsl_fit.c:456
double nsl_fit_model_flat_param_deriv(unsigned int param, double x, double A, double b, double a, double weight)
Definition: nsl_fit.c:443
double nsl_fit_model_levy_param_deriv(unsigned int param, double x, double A, double g, double mu, double weight)
Definition: nsl_fit.c:478
double nsl_fit_model_logarithmic_param_deriv(unsigned int param, double k, double A, double p, double weight)
Definition: nsl_fit.c:667
double nsl_fit_model_logistic_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight)
Definition: nsl_fit.c:202
double nsl_fit_model_gumbel2_param_deriv(unsigned int param, double x, double A, double a, double b, double mu, double weight)
Definition: nsl_fit.c:596
double nsl_fit_map_unbound(double x, double min, double max)
Definition: nsl_fit.c:89
double nsl_fit_model_fdist_param_deriv(unsigned int param, double x, double A, double n1, double n2, double weight)
Definition: nsl_fit.c:515
double nsl_fit_model_poisson_param_deriv(unsigned int param, double x, double A, double l, double weight)
Definition: nsl_fit.c:409
double nsl_fit_model_hill_param_deriv(unsigned int param, double x, double A, double n, double s, double weight)
Definition: nsl_fit.c:309
double nsl_fit_model_erf_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight)
Definition: nsl_fit.c:298
const char * nsl_fit_model_basic_equation[]
Definition: nsl_fit.c:41
double nsl_fit_model_laplace_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight)
Definition: nsl_fit.c:373
double nsl_fit_model_beta_param_deriv(unsigned int param, double x, double A, double a, double b, double weight)
Definition: nsl_fit.c:528
double nsl_fit_model_gudermann_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight)
Definition: nsl_fit.c:330
double nsl_fit_model_pascal_param_deriv(unsigned int param, double k, double A, double p, double n, double weight)
Definition: nsl_fit.c:640
double nsl_fit_model_negative_binomial_param_deriv(unsigned int param, double k, double A, double p, double n, double weight)
Definition: nsl_fit.c:626
double nsl_fit_model_gompertz_param_deriv(unsigned int param, double x, double a, double b, double c, double weight)
Definition: nsl_fit.c:320
double nsl_fit_model_landau_param_deriv(unsigned int param, double x, double weight)
Definition: nsl_fit.c:490
double nsl_fit_model_fourier_param_deriv(unsigned int param, unsigned int degree, double x, double w, double weight)
Definition: nsl_fit.c:156
double nsl_fit_model_gumbel1_param_deriv(unsigned int param, double x, double A, double s, double mu, double b, double weight)
Definition: nsl_fit.c:582
double nsl_fit_model_geometric_param_deriv(unsigned int param, double k, double A, double p, double weight)
Definition: nsl_fit.c:643
double nsl_fit_model_voigt_param_deriv(unsigned int param, double x, double a, double mu, double s, double g, double weight)
Definition: nsl_fit.c:215
double nsl_fit_model_tanh_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight)
Definition: nsl_fit.c:265
double nsl_fit_model_sech_dist_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight)
Definition: nsl_fit.c:675
double nsl_fit_model_chi_square_param_deriv(unsigned int param, double x, double A, double n, double weight)
Definition: nsl_fit.c:496
double nsl_fit_model_atan_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight)
Definition: nsl_fit.c:254
const char * nsl_fit_model_category_name[]
Definition: nsl_fit.c:37
double nsl_fit_model_exponentialn_param_deriv(unsigned int param, double x, double *p, double weight)
Definition: nsl_fit.c:140
double nsl_fit_model_lorentz_param_deriv(unsigned int param, double x, double A, double s, double t, double weight)
Definition: nsl_fit.c:178
double nsl_fit_model_maxwell_param_deriv(unsigned int param, double x, double a, double s, double weight)
Definition: nsl_fit.c:399
double nsl_fit_model_pareto_param_deriv(unsigned int param, double x, double A, double a, double b, double weight)
Definition: nsl_fit.c:540
double nsl_fit_model_rayleigh_tail_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight)
Definition: nsl_fit.c:466
double nsl_fit_model_sigmoid_param_deriv(unsigned int param, double x, double A, double mu, double k, double weight)
Definition: nsl_fit.c:287
const char * nsl_fit_model_peak_equation[]
Definition: nsl_fit.c:46
double nsl_fit_model_hypergeometric_param_deriv(unsigned int param, double k, double A, double n1, double n2, double t, double weight)
Definition: nsl_fit.c:651
double nsl_fit_model_exp_pow_param_deriv(unsigned int param, double x, double a, double s, double b, double mu, double weight)
Definition: nsl_fit.c:385
double nsl_fit_map_bound(double x, double min, double max)
Definition: nsl_fit.c:62
double nsl_fit_model_exponential_param_deriv(unsigned int param, double x, double A, double l, double mu, double weight)
Definition: nsl_fit.c:359
double nsl_fit_model_power1_param_deriv(unsigned int param, double x, double a, double b, double weight)
Definition: nsl_fit.c:122
double nsl_fit_model_students_t_param_deriv(unsigned int param, double x, double A, double n, double weight)
Definition: nsl_fit.c:506
double nsl_fit_model_frechet_param_deriv(unsigned int param, double x, double A, double g, double s, double mu, double weight)
Definition: nsl_fit.c:568
double nsl_fit_model_lognormal_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight)
Definition: nsl_fit.c:419
double nsl_fit_model_algebraic_sigmoid_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight)
Definition: nsl_fit.c:276
double nsl_fit_model_power2_param_deriv(unsigned int param, double x, double b, double c, double weight)
Definition: nsl_fit.c:130
double nsl_fit_model_gaussian_tail_param_deriv(unsigned int param, double x, double A, double s, double a, double mu, double weight)
Definition: nsl_fit.c:343
double nsl_fit_model_polynomial_param_deriv(double x, int j, double weight)
Definition: nsl_fit.c:119
double nsl_fit_model_inverse_exponential_param_deriv(unsigned int param, double x, double a, double b, double weight)
Definition: nsl_fit.c:146
nsl_fit_weight_type
Definition: nsl_fit.h:56
@ nsl_fit_weight_no
Definition: nsl_fit.h:56
@ nsl_fit_weight_relative_fit
Definition: nsl_fit.h:63
@ nsl_fit_weight_inverse
Definition: nsl_fit.h:59
@ nsl_fit_weight_statistical
Definition: nsl_fit.h:60
@ nsl_fit_weight_relative
Definition: nsl_fit.h:62
@ nsl_fit_weight_direct
Definition: nsl_fit.h:58
@ nsl_fit_weight_instrumental
Definition: nsl_fit.h:57
@ nsl_fit_weight_statistical_fit
Definition: nsl_fit.h:61
@ nsl_fit_model_power
Definition: nsl_fit.h:36
@ nsl_fit_model_fourier
Definition: nsl_fit.h:37
@ nsl_fit_model_polynomial
Definition: nsl_fit.h:36
@ nsl_fit_model_inverse_exponential
Definition: nsl_fit.h:36
@ nsl_fit_model_exponential
Definition: nsl_fit.h:36
@ nsl_fit_model_pseudovoigt1
Definition: nsl_fit.h:40
@ nsl_fit_model_lorentz
Definition: nsl_fit.h:40
@ nsl_fit_model_sech
Definition: nsl_fit.h:40
@ nsl_fit_model_voigt
Definition: nsl_fit.h:40
@ nsl_fit_model_logistic
Definition: nsl_fit.h:40
@ nsl_fit_model_gaussian
Definition: nsl_fit.h:40
#define NSL_FIT_MODEL_BASIC_COUNT
Definition: nsl_fit.h:35
@ nsl_fit_model_sigmoid
Definition: nsl_fit.h:43
@ nsl_fit_model_atan
Definition: nsl_fit.h:43
@ nsl_fit_model_algebraic_sigmoid
Definition: nsl_fit.h:43
@ nsl_fit_model_gompertz
Definition: nsl_fit.h:44
@ nsl_fit_model_gudermann
Definition: nsl_fit.h:44
@ nsl_fit_model_hill
Definition: nsl_fit.h:44
@ nsl_fit_model_erf
Definition: nsl_fit.h:43
@ nsl_fit_model_tanh
Definition: nsl_fit.h:43
nsl_fit_model_category
Definition: nsl_fit.h:33
@ nsl_fit_model_growth
Definition: nsl_fit.h:33
@ nsl_fit_model_distribution
Definition: nsl_fit.h:33
@ nsl_fit_model_basic
Definition: nsl_fit.h:33
@ nsl_fit_model_custom
Definition: nsl_fit.h:33
@ nsl_fit_model_peak
Definition: nsl_fit.h:33
const char * nsl_sf_stats_distribution_equation[]
Definition: nsl_sf_stats.c:43
@ nsl_sf_stats_beta
Definition: nsl_sf_stats.h:38
@ nsl_sf_stats_flat
Definition: nsl_sf_stats.h:37
@ nsl_sf_stats_maxwell_boltzmann
Definition: nsl_sf_stats.h:41
@ nsl_sf_stats_sech
Definition: nsl_sf_stats.h:41
@ nsl_sf_stats_poisson
Definition: nsl_sf_stats.h:39
@ nsl_sf_stats_levy_skew_alpha_stable
Definition: nsl_sf_stats.h:37
@ nsl_sf_stats_lognormal
Definition: nsl_sf_stats.h:37
@ nsl_sf_stats_fdist
Definition: nsl_sf_stats.h:38
@ nsl_sf_stats_binomial
Definition: nsl_sf_stats.h:39
@ nsl_sf_stats_negative_binomial
Definition: nsl_sf_stats.h:39
@ nsl_sf_stats_gumbel2
Definition: nsl_sf_stats.h:39
@ nsl_sf_stats_levy_alpha_stable
Definition: nsl_sf_stats.h:36
@ nsl_sf_stats_logarithmic
Definition: nsl_sf_stats.h:40
@ nsl_sf_stats_gaussian
Definition: nsl_sf_stats.h:35
@ nsl_sf_stats_geometric
Definition: nsl_sf_stats.h:40
@ nsl_sf_stats_exponential_power
Definition: nsl_sf_stats.h:35
@ nsl_sf_stats_laplace
Definition: nsl_sf_stats.h:35
@ nsl_sf_stats_gaussian_tail
Definition: nsl_sf_stats.h:35
@ nsl_sf_stats_landau
Definition: nsl_sf_stats.h:36
@ nsl_sf_stats_gamma
Definition: nsl_sf_stats.h:37
@ nsl_sf_stats_bernoulli
Definition: nsl_sf_stats.h:39
@ nsl_sf_stats_pareto
Definition: nsl_sf_stats.h:38
@ nsl_sf_stats_chi_squared
Definition: nsl_sf_stats.h:37
@ nsl_sf_stats_exponential
Definition: nsl_sf_stats.h:35
@ nsl_sf_stats_logistic
Definition: nsl_sf_stats.h:38
@ nsl_sf_stats_pascal
Definition: nsl_sf_stats.h:40
@ nsl_sf_stats_levy
Definition: nsl_sf_stats.h:41
@ nsl_sf_stats_tdist
Definition: nsl_sf_stats.h:38
@ nsl_sf_stats_weibull
Definition: nsl_sf_stats.h:38
@ nsl_sf_stats_gumbel1
Definition: nsl_sf_stats.h:39
@ nsl_sf_stats_rayleigh
Definition: nsl_sf_stats.h:36
@ nsl_sf_stats_frechet
Definition: nsl_sf_stats.h:41
@ nsl_sf_stats_hypergeometric
Definition: nsl_sf_stats.h:40
@ nsl_sf_stats_rayleigh_tail
Definition: nsl_sf_stats.h:36
@ nsl_sf_stats_cauchy_lorentz
Definition: nsl_sf_stats.h:36
double nsl_stats_bic(double sse, size_t n, size_t np, int version)
Definition: nsl_stats.c:264
double nsl_stats_tdist_t(double parameter, double error)
Definition: nsl_stats.c:191
double nsl_stats_fdist_F(double rsquare, size_t np, size_t dof)
Definition: nsl_stats.c:217
double nsl_stats_rsquare(double sse, double sst)
Definition: nsl_stats.c:176
double nsl_stats_logLik(double sse, size_t n)
Definition: nsl_stats.c:231
double nsl_stats_tdist_margin(double alpha, double dof, double error)
Definition: nsl_stats.c:204
double nsl_stats_aic(double sse, size_t n, size_t np, int version)
Definition: nsl_stats.c:237
double nsl_stats_chisq_p(double t, double dof)
Definition: nsl_stats.c:209
double nsl_stats_tdist_p(double t, double dof)
Definition: nsl_stats.c:198
double nsl_stats_fdist_p(double F, size_t np, double dof)
Definition: nsl_stats.c:223
double nsl_stats_rsquareAdj(double rsquare, size_t np, size_t dof, int version)
Definition: nsl_stats.c:179
double parse(const char *string, const char *locale)
symbol * assign_symbol(const char *symbol_name, double value)
int parse_errors(void)
nsl_fit_model_category modelCategory
Definition: XYFitCurve.h:50
QVector< double > paramLowerLimits
Definition: XYFitCurve.h:59
size_t evaluatedPoints
Definition: XYFitCurve.h:65
QStringList paramNamesUtf8
Definition: XYFitCurve.h:57
QVector< bool > paramFixed
Definition: XYFitCurve.h:61
nsl_fit_weight_type yWeightsType
Definition: XYFitCurve.h:53
nsl_fit_weight_type xWeightsType
Definition: XYFitCurve.h:52
double confidenceInterval
Definition: XYFitCurve.h:69
Range< double > fitRange
Definition: XYFitCurve.h:73
QStringList paramNames
Definition: XYFitCurve.h:56
Range< double > evalRange
Definition: XYFitCurve.h:74
QVector< double > paramUpperLimits
Definition: XYFitCurve.h:60
QVector< double > paramStartValues
Definition: XYFitCurve.h:58
QVector< double > tdist_pValues
Definition: XYFitCurve.h:106
QVector< double > errorValues
Definition: XYFitCurve.h:104
QVector< double > tdist_tValues
Definition: XYFitCurve.h:105
QVector< double > tdist_marginValues
Definition: XYFitCurve.h:107
QVector< double > paramValues
Definition: XYFitCurve.h:103
QVector< double > correlationMatrix
Definition: XYFitCurve.h:108
double * y
Definition: XYFitCurve.cpp:734
nsl_fit_model_category modelCategory
Definition: XYFitCurve.cpp:736
int degree
Definition: XYFitCurve.cpp:738
bool * paramFixed
Definition: XYFitCurve.cpp:743
size_t n
Definition: XYFitCurve.cpp:732
double * weight
Definition: XYFitCurve.cpp:735
double * paramMax
Definition: XYFitCurve.cpp:742
double * paramMin
Definition: XYFitCurve.cpp:741
QString * func
Definition: XYFitCurve.cpp:739
double * x
Definition: XYFitCurve.cpp:733
int modelType
Definition: XYFitCurve.cpp:737
QStringList * paramNames
Definition: XYFitCurve.cpp:740