diff --git a/src/backend/gsl/ExpressionParser.cpp b/src/backend/gsl/ExpressionParser.cpp index ca144ee7a..d723e14c1 100644 --- a/src/backend/gsl/ExpressionParser.cpp +++ b/src/backend/gsl/ExpressionParser.cpp @@ -1,1584 +1,1586 @@ /*************************************************************************** File : ExpressionParser.cpp Project : LabPlot -------------------------------------------------------------------- Copyright : (C) 2014 Alexander Semke (alexander.semke@web.de) Copyright : (C) 2014-2018 Stefan Gerlach (stefan.gerlach@uni.kn) Description : C++ wrapper for the bison generated parser. ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #include "backend/lib/macros.h" #include "backend/gsl/ExpressionParser.h" #include #include #include extern "C" { #include #include #include #include #include #include "backend/gsl/parser.h" } ExpressionParser* ExpressionParser::instance = NULL; ExpressionParser::ExpressionParser() { init_table(); initFunctions(); initConstants(); } void ExpressionParser::initFunctions() { //functions (sync with functions.h!) for (int i = 0; _functions[i].name != 0; i++) m_functions << _functions[i].name; m_functionsGroups << i18n("Standard Mathematical functions"); //http://www.gnu.org/software/gsl/manual/html_node/Special-Functions.html m_functionsGroups << i18n("Airy Functions and Derivatives"); m_functionsGroups << i18n("Bessel Functions"); m_functionsGroups << i18n("Clausen Functions"); m_functionsGroups << i18n("Coulomb Functions"); // m_functionsGroups << i18n("Coupling Coefficients"); m_functionsGroups << i18n("Dawson Function"); m_functionsGroups << i18n("Debye Functions"); m_functionsGroups << i18n("Dilogarithm"); // m_functionsGroups << i18n("Elementary Operations"); m_functionsGroups << i18n("Elliptic Integrals"); // m_functionsGroups << i18n("Elliptic Functions (Jacobi)"); #ifndef _MSC_VER m_functionsGroups << i18n("Error Functions and Related Functions"); #else m_functionsGroups << i18n("Error Functions"); #endif m_functionsGroups << i18n("Exponential Functions"); m_functionsGroups << i18n("Exponential Integrals"); m_functionsGroups << i18n("Fermi-Dirac Function"); m_functionsGroups << i18n("Gamma and Beta Functions"); m_functionsGroups << i18n("Gegenbauer Functions"); #if (GSL_MAJOR_VERSION > 2) || (GSL_MAJOR_VERSION == 2) && (GSL_MINOR_VERSION >= 4) m_functionsGroups << i18n("Hermite Polynomials and Functions"); #endif m_functionsGroups << i18n("Hypergeometric Functions"); m_functionsGroups << i18n("Laguerre Functions"); m_functionsGroups << i18n("Lambert W Functions"); m_functionsGroups << i18n("Legendre Functions and Spherical Harmonics"); m_functionsGroups << i18n("Logarithm and Related Functions"); // m_functionsGroups << i18n("Mathieu Functions"); m_functionsGroups << i18n("Power Function"); m_functionsGroups << i18n("Psi (Digamma) Function"); m_functionsGroups << i18n("Synchrotron Functions"); m_functionsGroups << i18n("Transport Functions"); m_functionsGroups << i18n("Trigonometric Functions"); m_functionsGroups << i18n("Zeta Functions"); // GSL random distribution functions m_functionsGroups << i18n("Gaussian Distribution"); m_functionsGroups << i18n("Exponential Distribution"); m_functionsGroups << i18n("Laplace Distribution"); m_functionsGroups << i18n("Exponential Power Distribution"); m_functionsGroups << i18n("Cauchy Distribution"); m_functionsGroups << i18n("Rayleigh Distribution"); m_functionsGroups << i18n("Landau Distribution"); m_functionsGroups << i18n("Gamma Distribution"); m_functionsGroups << i18n("Flat (Uniform) Distribution"); m_functionsGroups << i18n("Lognormal Distribution"); m_functionsGroups << i18n("Chi-squared Distribution"); m_functionsGroups << i18n("F-distribution"); m_functionsGroups << i18n("t-distribution"); m_functionsGroups << i18n("Beta Distribution"); m_functionsGroups << i18n("Logistic Distribution"); m_functionsGroups << i18n("Pareto Distribution"); m_functionsGroups << i18n("Weibull Distribution"); m_functionsGroups << i18n("Gumbel Distribution"); m_functionsGroups << i18n("Poisson Distribution"); m_functionsGroups << i18n("Bernoulli Distribution"); m_functionsGroups << i18n("Binomial Distribution"); m_functionsGroups << i18n("Pascal Distribution"); m_functionsGroups << i18n("Geometric Distribution"); m_functionsGroups << i18n("Hypergeometric Distribution"); m_functionsGroups << i18n("Logarithmic Distribution"); int index = 0; // Standard mathematical functions m_functionsNames << i18n("pseudo-random integer [0,RAND_MAX]"); m_functionsNames << i18n("nonlinear additive feedback rng [0,RAND_MAX]"); m_functionsNames << i18n("nonlinear additive feedback rng [0,1]"); m_functionsNames << i18n("Smallest integral value not less"); m_functionsNames << i18n("Absolute value"); m_functionsNames << i18n("Base 10 logarithm"); m_functionsNames << i18n("Power function [x^y]"); m_functionsNames << i18n("Nonnegative square root"); m_functionsNames << i18n("Sign function"); m_functionsNames << i18n("Heavyside theta function"); m_functionsNames << i18n("Harmonic number function"); #ifndef HAVE_WINDOWS m_functionsNames << i18n("Cube root"); m_functionsNames << i18n("Extract the exponent"); m_functionsNames << i18n("Round to an integer value"); m_functionsNames << i18n("Round to the nearest integer"); m_functionsNames << i18n("Round to the nearest integer"); #endif m_functionsNames << QString("log(1+x)"); m_functionsNames << QString("x * 2^e"); m_functionsNames << QString("x^n"); m_functionsNames << QString("x^2"); m_functionsNames << QString("x^3"); m_functionsNames << QString("x^4"); m_functionsNames << QString("x^5"); m_functionsNames << QString("x^6"); m_functionsNames << QString("x^7"); m_functionsNames << QString("x^8"); m_functionsNames << QString("x^9"); #ifndef HAVE_WINDOWS for (int i = 0; i < 27; i++) #else for (int i = 0; i < 22; i++) #endif m_functionsGroupIndex << index; // Airy Functions and Derivatives m_functionsNames << i18n("Airy function of the first kind"); m_functionsNames << i18n("Airy function of the second kind"); m_functionsNames << i18n("Scaled Airy function of the first kind"); m_functionsNames << i18n("Scaled Airy function of the second kind"); m_functionsNames << i18n("Airy function derivative of the first kind"); m_functionsNames << i18n("Airy function derivative of the second kind"); m_functionsNames << i18n("Scaled Airy function derivative of the first kind"); m_functionsNames << i18n("Scaled Airy function derivative of the second kind"); m_functionsNames << i18n("n-th zero of the Airy function of the first kind"); m_functionsNames << i18n("n-th zero of the Airy function of the second kind"); m_functionsNames << i18n("n-th zero of the Airy function derivative of the first kind"); m_functionsNames << i18n("n-th zero of the Airy function derivative of the second kind"); index++; for (int i = 0; i < 12; i++) m_functionsGroupIndex << index; // Bessel Functions m_functionsNames << i18n("Regular cylindrical Bessel function of zeroth order"); m_functionsNames << i18n("Regular cylindrical Bessel function of first order"); m_functionsNames << i18n("Regular cylindrical Bessel function of order n"); m_functionsNames << i18n("Irregular cylindrical Bessel function of zeroth order"); m_functionsNames << i18n("Irregular cylindrical Bessel function of first order"); m_functionsNames << i18n("Irregular cylindrical Bessel function of order n"); m_functionsNames << i18n("Regular modified cylindrical Bessel function of zeroth order"); m_functionsNames << i18n("Regular modified cylindrical Bessel function of first order"); m_functionsNames << i18n("Regular modified cylindrical Bessel function of order n"); m_functionsNames << i18n("Scaled regular modified cylindrical Bessel function of zeroth order exp(-|x|) I0(x)"); m_functionsNames << i18n("Scaled regular modified cylindrical Bessel function of first order exp(-|x|) I1(x)"); m_functionsNames << i18n("Scaled regular modified cylindrical Bessel function of order n exp(-|x|) In(x)"); m_functionsNames << i18n("Irregular modified cylindrical Bessel function of zeroth order"); m_functionsNames << i18n("Irregular modified cylindrical Bessel function of first order"); m_functionsNames << i18n("Irregular modified cylindrical Bessel function of order n"); m_functionsNames << i18n("Scaled irregular modified cylindrical Bessel function of zeroth order exp(x) K0(x)"); m_functionsNames << i18n("Scaled irregular modified cylindrical Bessel function of first order exp(x) K1(x)"); m_functionsNames << i18n("Scaled irregular modified cylindrical Bessel function of order n exp(x) Kn(x)"); m_functionsNames << i18n("Regular spherical Bessel function of zeroth order"); m_functionsNames << i18n("Regular spherical Bessel function of first order"); m_functionsNames << i18n("Regular spherical Bessel function of second order"); m_functionsNames << i18n("Regular spherical Bessel function of order l"); m_functionsNames << i18n("Irregular spherical Bessel function of zeroth order"); m_functionsNames << i18n("Irregular spherical Bessel function of first order"); m_functionsNames << i18n("Irregular spherical Bessel function of second order"); m_functionsNames << i18n("Irregular spherical Bessel function of order l"); m_functionsNames << i18n("Scaled regular modified spherical Bessel function of zeroth order, exp(-|x|) i0(x)"); m_functionsNames << i18n("Scaled regular modified spherical Bessel function of first order, exp(-|x|) i1(x)"); m_functionsNames << i18n("Scaled regular modified spherical Bessel function of second order, exp(-|x|) i2(x)"); m_functionsNames << i18n("Scaled regular modified spherical Bessel function of order l, exp(-|x|) il(x)"); m_functionsNames << i18n("Scaled irregular modified spherical Bessel function of zeroth order, exp(x) k0(x)"); m_functionsNames << i18n("Scaled irregular modified spherical Bessel function of first order, exp(-|x|) k1(x)"); m_functionsNames << i18n("Scaled irregular modified spherical Bessel function of second order, exp(-|x|) k2(x)"); m_functionsNames << i18n("Scaled irregular modified spherical Bessel function of order l, exp(-|x|) kl(x)"); m_functionsNames << i18n("Regular cylindrical Bessel function of fractional order"); m_functionsNames << i18n("Irregular cylindrical Bessel function of fractional order"); m_functionsNames << i18n("Regular modified Bessel function of fractional order"); m_functionsNames << i18n("Scaled regular modified Bessel function of fractional order"); m_functionsNames << i18n("Irregular modified Bessel function of fractional order"); m_functionsNames << i18n("Logarithm of irregular modified Bessel function of fractional order"); m_functionsNames << i18n("Scaled irregular modified Bessel function of fractional order"); m_functionsNames << i18n("n-th positive zero of the Bessel function J0"); m_functionsNames << i18n("n-th positive zero of the Bessel function J1"); m_functionsNames << i18n("n-th positive zero of the Bessel function Jnu"); index++; for (int i = 0; i < 44; i++) m_functionsGroupIndex << index; // Clausen Functions m_functionsNames << i18n("Clausen function"); index++; m_functionsGroupIndex << index; // Coulomb Functions m_functionsNames << i18n("Lowest-order normalized hydrogenic bound state radial wavefunction"); m_functionsNames << i18n("n-th normalized hydrogenic bound state radial wavefunction"); index++; for (int i = 0; i < 2; i++) m_functionsGroupIndex << index; // Dawson Function m_functionsNames << i18n("Dawson integral"); index++; m_functionsGroupIndex << index; // Debye Functions m_functionsNames << i18n("First-order Debye function"); m_functionsNames << i18n("Second-order Debye function"); m_functionsNames << i18n("Third-order Debye function"); m_functionsNames << i18n("Fourth-order Debye function"); m_functionsNames << i18n("Fifth-order Debye function"); m_functionsNames << i18n("Sixth-order Debye function"); index++; for (int i = 0; i < 6; i++) m_functionsGroupIndex << index; // Dilogarithm m_functionsNames << i18n("Dilogarithm for a real argument"); index++; m_functionsGroupIndex << index; // Elliptic Integrals m_functionsNames << i18n("Legendre form of complete elliptic integral K"); m_functionsNames << i18n("Legendre form of complete elliptic integral E"); m_functionsNames << i18n("Legendre form of complete elliptic integral Pi"); m_functionsNames << i18n("Legendre form of incomplete elliptic integral F"); m_functionsNames << i18n("Legendre form of incomplete elliptic integral E"); m_functionsNames << i18n("Legendre form of incomplete elliptic integral P"); m_functionsNames << i18n("Legendre form of incomplete elliptic integral D"); m_functionsNames << i18n("Carlson form of incomplete elliptic integral RC"); m_functionsNames << i18n("Carlson form of incomplete elliptic integral RD"); m_functionsNames << i18n("Carlson form of incomplete elliptic integral RF"); m_functionsNames << i18n("Carlson form of incomplete elliptic integral RJ"); index++; for (int i = 0; i < 11; i++) m_functionsGroupIndex << index; // Error Functions m_functionsNames << i18n("Error function"); m_functionsNames << i18n("Complementary error function"); m_functionsNames << i18n("Logarithm of complementary error function"); m_functionsNames << i18n("Gaussian probability density function Z"); m_functionsNames << i18n("Upper tail of the Gaussian probability function Q"); m_functionsNames << i18n("Hazard function for the normal distribution Z/Q"); int count = 6; #ifndef _MSC_VER m_functionsNames << i18n("Underflow-compensating function exp(x^2) erfc(x) for real x"); m_functionsNames << i18n("Imaginary error function erfi(x) = -i erf(ix) for real x"); m_functionsNames << i18n("Imaginary part of Faddeeva's scaled complex error function w(x) = exp(-x^2) erfc(-ix) for real x"); m_functionsNames << i18n("Dawson's integral D(z) = sqrt(pi)/2 * exp(-z^2) * erfi(z)"); m_functionsNames << i18n("Voigt profile"); count += 5; #endif + m_functionsNames << i18n("Pseudo-Voigt profile (same width)"); + count += 1; index++; for (int i = 0; i < count; i++) m_functionsGroupIndex << index; // Exponential Functions m_functionsNames << i18n("Exponential function"); m_functionsNames << i18n("exponentiate x and multiply by y"); m_functionsNames << QString("exp(x) - 1"); m_functionsNames << QString("(exp(x)-1)/x"); m_functionsNames << QString("2(exp(x)-1-x)/x^2"); m_functionsNames << i18n("n-relative exponential"); index++; for (int i = 0; i < 6; i++) m_functionsGroupIndex << index; // Exponential Integrals m_functionsNames << i18n("Exponential integral"); m_functionsNames << i18n("Second order exponential integral"); m_functionsNames << i18n("Exponential integral of order n"); m_functionsNames << i18n("Exponential integral Ei"); m_functionsNames << i18n("Hyperbolic integral Shi"); m_functionsNames << i18n("Hyperbolic integral Chi"); m_functionsNames << i18n("Third-order exponential integral"); m_functionsNames << i18n("Sine integral"); m_functionsNames << i18n("Cosine integral"); m_functionsNames << i18n("Arctangent integral"); index++; for (int i = 0; i < 10; i++) m_functionsGroupIndex << index; // Fermi-Dirac Function m_functionsNames << i18n("Complete Fermi-Dirac integral with index -1"); m_functionsNames << i18n("Complete Fermi-Dirac integral with index 0"); m_functionsNames << i18n("Complete Fermi-Dirac integral with index 1"); m_functionsNames << i18n("Complete Fermi-Dirac integral with index 2"); m_functionsNames << i18n("Complete Fermi-Dirac integral with integer index j"); m_functionsNames << i18n("Complete Fermi-Dirac integral with index -1/2"); m_functionsNames << i18n("Complete Fermi-Dirac integral with index 1/2"); m_functionsNames << i18n("Complete Fermi-Dirac integral with index 3/2"); m_functionsNames << i18n("Incomplete Fermi-Dirac integral with index zero"); index++; for (int i = 0; i < 9; i++) m_functionsGroupIndex << index; // Gamma and Beta Functions m_functionsNames << i18n("Gamma function"); m_functionsNames << i18n("Gamma function"); m_functionsNames << i18n("Logarithm of the gamma function"); m_functionsNames << i18n("Logarithm of the gamma function"); m_functionsNames << i18n("Regulated gamma function"); m_functionsNames << i18n("Reciprocal of the gamma function"); m_functionsNames << i18n("Factorial n!"); m_functionsNames << i18n("Double factorial n!!"); m_functionsNames << i18n("Logarithm of the factorial"); m_functionsNames << i18n("Logarithm of the double factorial"); m_functionsNames << i18n("Combinatorial factor"); m_functionsNames << i18n("Logarithm of the combinatorial factor"); m_functionsNames << i18n("Taylor coefficient"); m_functionsNames << i18n("Pochhammer symbol"); m_functionsNames << i18n("Logarithm of the Pochhammer symbol"); m_functionsNames << i18n("Relative Pochhammer symbol"); m_functionsNames << i18n("Unnormalized incomplete gamma function"); m_functionsNames << i18n("Normalized incomplete gamma function"); m_functionsNames << i18n("Complementary normalized incomplete gamma function"); m_functionsNames << i18n("Beta function"); m_functionsNames << i18n("Logarithm of the beta function"); m_functionsNames << i18n("Normalized incomplete beta function"); index++; for (int i = 0; i < 22; i++) m_functionsGroupIndex << index; // Gegenbauer Functions m_functionsNames << i18n("Gegenbauer polynomial C_1"); m_functionsNames << i18n("Gegenbauer polynomial C_2"); m_functionsNames << i18n("Gegenbauer polynomial C_3"); m_functionsNames << i18n("Gegenbauer polynomial C_n"); index++; for (int i = 0; i < 4; i++) m_functionsGroupIndex << index; #if (GSL_MAJOR_VERSION > 2) || (GSL_MAJOR_VERSION == 2) && (GSL_MINOR_VERSION >= 4) // Hermite Polynomials and Functions m_functionsNames << i18n("Hermite polynomials physicists version"); m_functionsNames << i18n("Hermite polynomials probabilists version"); m_functionsNames << i18n("Hermite functions"); m_functionsNames << i18n("Derivatives of Hermite polynomials physicists version"); m_functionsNames << i18n("Derivatives of Hermite polynomials probabilists version"); m_functionsNames << i18n("Derivatives of Hermite functions"); index++; for (int i = 0; i < 6; i++) m_functionsGroupIndex << index; #endif // Hypergeometric Functions m_functionsNames << i18n("Hypergeometric function 0F1"); m_functionsNames << i18n("Confluent hypergeometric function 1F1 for integer parameters"); m_functionsNames << i18n("Confluent hypergeometric function 1F1 for general parameters"); m_functionsNames << i18n("Confluent hypergeometric function U for integer parameters"); m_functionsNames << i18n("Confluent hypergeometric function U"); m_functionsNames << i18n("Gauss hypergeometric function 2F1"); m_functionsNames << i18n("Gauss hypergeometric function 2F1 with complex parameters"); m_functionsNames << i18n("Renormalized Gauss hypergeometric function 2F1"); m_functionsNames << i18n("Renormalized Gauss hypergeometric function 2F1 with complex parameters"); m_functionsNames << i18n("Hypergeometric function 2F0"); index++; for (int i = 0; i < 10; i++) m_functionsGroupIndex << index; // Laguerre Functions m_functionsNames << i18n("generalized Laguerre polynomials L_1"); m_functionsNames << i18n("generalized Laguerre polynomials L_2"); m_functionsNames << i18n("generalized Laguerre polynomials L_3"); index++; for (int i = 0; i < 3; i++) m_functionsGroupIndex << index; // Lambert W Functions m_functionsNames << i18n("Principal branch of the Lambert W function"); m_functionsNames << i18n("Secondary real-valued branch of the Lambert W function"); index++; for (int i = 0; i < 2; i++) m_functionsGroupIndex << index; // Legendre Functions and Spherical Harmonics m_functionsNames << i18n("Legendre polynomial P_1"); m_functionsNames << i18n("Legendre polynomial P_2"); m_functionsNames << i18n("Legendre polynomial P_3"); m_functionsNames << i18n("Legendre polynomial P_l"); m_functionsNames << i18n("Legendre function Q_0"); m_functionsNames << i18n("Legendre function Q_1"); m_functionsNames << i18n("Legendre function Q_l"); m_functionsNames << i18n("Associated Legendre polynomial"); m_functionsNames << i18n("Normalized associated Legendre polynomial"); m_functionsNames << i18n("Irregular spherical conical function P^1/2"); m_functionsNames << i18n("Regular spherical conical function P^(-1/2)"); m_functionsNames << i18n("Conical function P^0"); m_functionsNames << i18n("Conical function P^1"); m_functionsNames << i18n("Regular spherical conical function P^(-1/2-l)"); m_functionsNames << i18n("Regular cylindrical conical function P^(-m)"); m_functionsNames << i18n("Zeroth radial eigenfunction of the Laplacian on the 3-dimensional hyperbolic space"); m_functionsNames << i18n("First radial eigenfunction of the Laplacian on the 3-dimensional hyperbolic space"); m_functionsNames << i18n("l-th radial eigenfunction of the Laplacian on the 3-dimensional hyperbolic space"); index++; for (int i = 0; i < 18; i++) m_functionsGroupIndex << index; // Logarithm and Related Functions m_functionsNames << i18n("Logarithm"); m_functionsNames << i18n("Logarithm of the magnitude"); m_functionsNames << QString("log(1+x)"); m_functionsNames << QString("log(1+x) - x"); index++; for (int i = 0; i < 4; i++) m_functionsGroupIndex << index; // Power Function m_functionsNames << i18n("x^n for integer n with an error estimate"); index++; m_functionsGroupIndex << index; // Psi (Digamma) Function m_functionsNames << i18n("Digamma function for positive integer n"); m_functionsNames << i18n("Digamma function"); m_functionsNames << i18n("Real part of the digamma function on the line 1+i y"); m_functionsNames << i18n("Trigamma function psi' for positive integer n"); m_functionsNames << i18n("Trigamma function psi'"); m_functionsNames << i18n("Polygamma function psi^(n)"); index++; for (int i = 0; i < 6; i++) m_functionsGroupIndex << index; // Synchrotron Functions m_functionsNames << i18n("First synchrotron function"); m_functionsNames << i18n("Second synchrotron function"); index++; for (int i = 0; i < 2; i++) m_functionsGroupIndex << index; // Transport Functions m_functionsNames << i18n("Transport function"); m_functionsNames << i18n("Transport function"); m_functionsNames << i18n("Transport function"); m_functionsNames << i18n("Transport function"); index++; for (int i = 0; i < 4; i++) m_functionsGroupIndex << index; // Trigonometric Functions m_functionsNames << i18n("Sine"); m_functionsNames << i18n("Cosine"); m_functionsNames << i18n("Tangent"); m_functionsNames << i18n("Inverse sine"); m_functionsNames << i18n("Inverse cosine"); m_functionsNames << i18n("Inverse tangent"); m_functionsNames << i18n("Inverse tangent using sign"); m_functionsNames << i18n("Hyperbolic sine"); m_functionsNames << i18n("Hyperbolic cosine"); m_functionsNames << i18n("Hyperbolic tangent"); m_functionsNames << i18n("Inverse hyperbolic cosine"); m_functionsNames << i18n("Inverse hyperbolic sine"); m_functionsNames << i18n("Inverse hyperbolic tangent"); m_functionsNames << i18n("Secant"); m_functionsNames << i18n("Cosecant"); m_functionsNames << i18n("Cotangent"); m_functionsNames << i18n("Inverse secant"); m_functionsNames << i18n("Inverse cosecant"); m_functionsNames << i18n("Inverse cotangent"); m_functionsNames << i18n("Hyperbolic secant"); m_functionsNames << i18n("Hyperbolic cosecant"); m_functionsNames << i18n("Hyperbolic cotangent"); m_functionsNames << i18n("Inverse hyperbolic secant"); m_functionsNames << i18n("Inverse hyperbolic cosecant"); m_functionsNames << i18n("Inverse hyperbolic cotangent"); m_functionsNames << i18n("Sinc function sin(x)/x"); m_functionsNames << QString("log(sinh(x))"); m_functionsNames << QString("log(cosh(x))"); m_functionsNames << i18n("Hypotenuse function"); m_functionsNames << i18n("Three component hypotenuse function"); m_functionsNames << i18n("restrict to [-pi,pi]"); m_functionsNames << i18n("restrict to [0,2 pi]"); index++; for (int i = 0; i < 32; i++) m_functionsGroupIndex << index; // Zeta Functions m_functionsNames << i18n("Riemann zeta function for integer n"); m_functionsNames << i18n("Riemann zeta function"); m_functionsNames << i18n("zeta(n)-1 for integer n"); m_functionsNames << i18n("zeta(x)-1"); m_functionsNames << i18n("Hurwitz zeta function"); m_functionsNames << i18n("Eta function for integer n"); m_functionsNames << i18n("Eta function"); index++; for (int i = 0; i < 7; i++) m_functionsGroupIndex << index; // GSL Random Number Distributions: see http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Distributions.html // Gaussian Distribution m_functionsNames << i18n("Probability density for a Gaussian distribution"); m_functionsNames << i18n("Probability density for a unit Gaussian distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); m_functionsNames << i18n("Cumulative unit distribution function P"); m_functionsNames << i18n("Cumulative unit distribution function Q"); m_functionsNames << i18n("Inverse cumulative unit distribution function P"); m_functionsNames << i18n("Inverse cumulative unit distribution function Q"); m_functionsNames << i18n("Probability density for Gaussian tail distribution"); m_functionsNames << i18n("Probability density for unit Gaussian tail distribution"); m_functionsNames << i18n("Probability density for a bivariate Gaussian distribution"); index++; for (int i = 0; i < 13; i++) m_functionsGroupIndex << index; // Exponential Distribution m_functionsNames << i18n("Probability density for an exponential distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 5; i++) m_functionsGroupIndex << index; // Laplace Distribution m_functionsNames << i18n("Probability density for a Laplace distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 5; i++) m_functionsGroupIndex << index; // Exponential Power Distribution m_functionsNames << i18n("Probability density for an exponential power distribution"); m_functionsNames << i18n("cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); index++; for (int i = 0; i < 3; i++) m_functionsGroupIndex << index; // Cauchy Distribution m_functionsNames << i18n("Probability density for a Cauchy distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 5; i++) m_functionsGroupIndex << index; // Rayleigh Distribution m_functionsNames << i18n("Probability density for a Rayleigh distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); m_functionsNames << i18n("Probability density for a Rayleigh tail distribution"); index++; for (int i = 0; i < 6; i++) m_functionsGroupIndex << index; // Landau Distribution m_functionsNames << i18n("Probability density for a Landau distribution"); index++; m_functionsGroupIndex << index; // Gamma Distribution m_functionsNames << i18n("Probability density for a gamma distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 5; i++) m_functionsGroupIndex << index; // Flat (Uniform) Distribution m_functionsNames << i18n("Probability density for a uniform distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 5; i++) m_functionsGroupIndex << index; // Lognormal Distribution m_functionsNames << i18n("Probability density for a lognormal distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 5; i++) m_functionsGroupIndex << index; // Chi-squared Distribution m_functionsNames << i18n("Probability density for a chi squared distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 5; i++) m_functionsGroupIndex << index; // F-distribution m_functionsNames << i18n("Probability density for a F-distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 5; i++) m_functionsGroupIndex << index; // t-distribution m_functionsNames << i18n("Probability density for a t-distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 5; i++) m_functionsGroupIndex << index; // Beta Distribution m_functionsNames << i18n("Probability density for a beta distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 5; i++) m_functionsGroupIndex << index; // Logistic Distribution m_functionsNames << i18n("Probability density for a logistic distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 5; i++) m_functionsGroupIndex << index; // Pareto Distribution m_functionsNames << i18n("Probability density for a Pareto distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 5; i++) m_functionsGroupIndex << index; // Weibull Distribution m_functionsNames << i18n("Probability density for a Weibull distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 5; i++) m_functionsGroupIndex << index; // Gumbel Distribution m_functionsNames << i18n("Probability density for a Type-1 Gumbel distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); m_functionsNames << i18n("Probability density for a Type-2 Gumbel distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Inverse cumulative distribution function P"); m_functionsNames << i18n("Inverse cumulative distribution function Q"); index++; for (int i = 0; i < 10; i++) m_functionsGroupIndex << index; // Poisson Distribution m_functionsNames << i18n("Probability density for a Poisson distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); index++; for (int i = 0; i < 3; i++) m_functionsGroupIndex << index; // Bernoulli Distribution m_functionsNames << i18n("Probability density for a Bernoulli distribution"); index++; m_functionsGroupIndex << index; // Binomial Distribution m_functionsNames << i18n("Probability density for a binomial distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); m_functionsNames << i18n("Probability density for a negative binomial distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); index++; for (int i = 0; i < 6; i++) m_functionsGroupIndex << index; // Pascal Distribution m_functionsNames << i18n("Probability density for a Pascal distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); index++; for (int i = 0; i < 3; i++) m_functionsGroupIndex << index; // Geometric Distribution m_functionsNames << i18n("Probability density for a geometric distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); index++; for (int i = 0; i < 3; i++) m_functionsGroupIndex << index; // Hypergeometric Distribution m_functionsNames << i18n("Probability density for a hypergeometric distribution"); m_functionsNames << i18n("Cumulative distribution function P"); m_functionsNames << i18n("Cumulative distribution function Q"); index++; for (int i = 0; i < 3; i++) m_functionsGroupIndex << index; // Logarithmic Distribution m_functionsNames << i18n("Probability density for a logarithmic distribution"); index++; m_functionsGroupIndex << index; } //TODO: decide whether we want to have i18n here in the backend part of the code void ExpressionParser::initConstants() { for (int i = 0; _constants[i].name != 0; i++) m_constants << _constants[i].name; //groups m_constantsGroups << i18n("Mathematical constants"); m_constantsGroups << i18n("Fundamental constants"); m_constantsGroups << i18n("Astronomy and Astrophysics"); m_constantsGroups << i18n("Atomic and Nuclear Physics"); m_constantsGroups << i18n("Measurement of Time"); m_constantsGroups << i18n("Imperial Units"); m_constantsGroups << i18n("Speed and Nautical Units"); m_constantsGroups << i18n("Printers Units"); m_constantsGroups << i18n("Volume, Area and Length"); m_constantsGroups << i18n("Mass and Weight"); m_constantsGroups << i18n("Thermal Energy and Power"); m_constantsGroups << i18n("Pressure"); m_constantsGroups << i18n("Viscosity"); m_constantsGroups << i18n("Light and Illumination"); m_constantsGroups << i18n("Radioactivity"); m_constantsGroups << i18n("Force and Energy"); //Mathematical constants m_constantsNames << i18n("Base of exponentials"); m_constantsValues << QString::number(M_E,'g',15); m_constantsUnits << ""; m_constantsNames << i18n("Pi"); m_constantsValues << QString::number(M_PI,'g',15); m_constantsUnits << ""; m_constantsNames << i18n("Euler's constant"); m_constantsValues << QString::number(M_EULER,'g',15); m_constantsUnits << ""; for (int i = 0; i < 3; i++) m_constantsGroupIndex << 0; //Fundamental constants m_constantsNames << i18n("Speed of light"); m_constantsValues << QString::number(GSL_CONST_MKSA_SPEED_OF_LIGHT,'g',15); m_constantsUnits << "m / s"; m_constantsNames << i18n("Vacuum permeability"); m_constantsValues << QString::number(GSL_CONST_MKSA_VACUUM_PERMEABILITY,'g',15); m_constantsUnits << "kg m / A^2 s^2"; m_constantsNames << i18n("Vacuum permittivity"); m_constantsValues << QString::number(GSL_CONST_MKSA_VACUUM_PERMITTIVITY,'g',15); m_constantsUnits << "A^2 s^4 / kg m^3"; m_constantsNames << i18n("Planck constant"); m_constantsValues << QString::number(GSL_CONST_MKSA_PLANCKS_CONSTANT_H,'g',15); m_constantsUnits << "kg m^2 / s"; m_constantsNames << i18n("Reduced Planck constant"); m_constantsValues << QString::number(GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR,'g',15); m_constantsUnits << "kg m^2 / s"; m_constantsNames << i18n("Avogadro constant"); m_constantsValues << QString::number(GSL_CONST_NUM_AVOGADRO,'g',15); m_constantsUnits << "1 / mol"; m_constantsNames << i18n("Faraday"); m_constantsValues << QString::number(GSL_CONST_MKSA_FARADAY,'g',15); m_constantsUnits << "A s / mol"; m_constantsNames << i18n("Boltzmann constant"); m_constantsValues << QString::number(GSL_CONST_MKSA_BOLTZMANN,'g',15); m_constantsUnits << "kg m^2 / K s^2"; m_constantsNames << i18n("Molar gas"); m_constantsValues << QString::number(GSL_CONST_MKSA_MOLAR_GAS,'g',15); m_constantsUnits << "kg m^2 / K mol s^2"; m_constantsNames << i18n("Standard gas volume"); m_constantsValues << QString::number(GSL_CONST_MKSA_STANDARD_GAS_VOLUME,'g',15); m_constantsUnits << "m^3 / mol"; m_constantsNames << i18n("Stefan-Boltzmann constant"); m_constantsValues << QString::number(GSL_CONST_MKSA_STEFAN_BOLTZMANN_CONSTANT,'g',15); m_constantsUnits << "kg / K^4 s^3"; m_constantsNames << i18n("Gauss"); m_constantsValues << QString::number(GSL_CONST_MKSA_GAUSS,'g',15); m_constantsUnits << "kg / A s^2"; for (int i = 0; i < 12; i++) m_constantsGroupIndex << 1; // Astronomy and Astrophysics m_constantsNames << i18n("Astronomical unit"); m_constantsValues << QString::number(GSL_CONST_MKSA_ASTRONOMICAL_UNIT,'g',15); m_constantsUnits << "m"; m_constantsNames << i18n("Gravitational constant"); m_constantsValues << QString::number(GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT,'g',15); m_constantsUnits << "m^3 / kg s^2"; m_constantsNames << i18n("Light year"); m_constantsValues << QString::number(GSL_CONST_MKSA_LIGHT_YEAR,'g',15); m_constantsUnits << "m"; m_constantsNames << i18n("Parsec"); m_constantsValues << QString::number(GSL_CONST_MKSA_PARSEC,'g',15); m_constantsUnits << "m"; m_constantsNames << i18n("Gravitational acceleration"); m_constantsValues << QString::number(GSL_CONST_MKSA_GRAV_ACCEL,'g',15); m_constantsUnits << "m / s^2"; m_constantsNames << i18n("Solar mass"); m_constantsValues << QString::number(GSL_CONST_MKSA_SOLAR_MASS,'g',15); m_constantsUnits << "kg"; for (int i = 0; i < 6; i++) m_constantsGroupIndex << 2; // Atomic and Nuclear Physics; m_constantsNames << i18n("Charge of the electron"); m_constantsValues << QString::number(GSL_CONST_MKSA_ELECTRON_CHARGE,'g',15); m_constantsUnits << "A s"; m_constantsNames << i18n("Energy of 1 electron volt"); m_constantsValues << QString::number(GSL_CONST_MKSA_ELECTRON_VOLT,'g',15); m_constantsUnits << "kg m^2 / s^2"; m_constantsNames << i18n("Unified atomic mass"); m_constantsValues << QString::number(GSL_CONST_MKSA_UNIFIED_ATOMIC_MASS,'g',15); m_constantsUnits << "kg"; m_constantsNames << i18n("Mass of the electron"); m_constantsValues << QString::number(GSL_CONST_MKSA_MASS_ELECTRON,'g',15); m_constantsUnits << "kg"; m_constantsNames << i18n("Mass of the muon"); m_constantsValues << QString::number(GSL_CONST_MKSA_MASS_MUON,'g',15); m_constantsUnits << "kg"; m_constantsNames << i18n("Mass of the proton"); m_constantsValues << QString::number(GSL_CONST_MKSA_MASS_PROTON,'g',15); m_constantsUnits << "kg"; m_constantsNames << i18n("Mass of the neutron"); m_constantsValues << QString::number(GSL_CONST_MKSA_MASS_NEUTRON,'g',15); m_constantsUnits << "kg"; m_constantsNames << i18n("Electromagnetic fine structure constant"); m_constantsValues << QString::number(GSL_CONST_NUM_FINE_STRUCTURE,'g',15); m_constantsUnits << ""; m_constantsNames << i18n("Rydberg constant"); m_constantsValues << QString::number(GSL_CONST_MKSA_RYDBERG,'g',15); m_constantsUnits << "kg m^2 / s^2"; m_constantsNames << i18n("Bohr radius"); m_constantsValues << QString::number(GSL_CONST_MKSA_BOHR_RADIUS,'g',15); m_constantsUnits << "m"; m_constantsNames << i18n("Length of 1 angstrom"); m_constantsValues << QString::number(GSL_CONST_MKSA_ANGSTROM,'g',15); m_constantsUnits << "m"; m_constantsNames << i18n("Area of 1 barn"); m_constantsValues << QString::number(GSL_CONST_MKSA_BARN,'g',15); m_constantsUnits << "m^2"; m_constantsNames << i18n("Bohr Magneton"); m_constantsValues << QString::number(GSL_CONST_MKSA_BOHR_MAGNETON,'g',15); m_constantsUnits << "A m^2"; m_constantsNames << i18n("Nuclear Magneton"); m_constantsValues << QString::number(GSL_CONST_MKSA_NUCLEAR_MAGNETON,'g',15); m_constantsUnits << "A m^2"; m_constantsNames << i18n("Magnetic moment of the electron [absolute value]"); m_constantsValues << QString::number(GSL_CONST_MKSA_ELECTRON_MAGNETIC_MOMENT,'g',15); m_constantsUnits << "A m^2"; m_constantsNames << i18n("Magnetic moment of the proton"); m_constantsValues << QString::number(GSL_CONST_MKSA_PROTON_MAGNETIC_MOMENT,'g',15); m_constantsUnits << "A m^2"; m_constantsNames << i18n("Thomson cross section"); m_constantsValues << QString::number(GSL_CONST_MKSA_THOMSON_CROSS_SECTION,'g',15); m_constantsUnits << "m^2"; m_constantsNames << i18n("Electric dipole moment of 1 Debye"); m_constantsValues << QString::number(GSL_CONST_MKSA_DEBYE,'g',15); m_constantsUnits << "A s^2 / m^2"; for (int i = 0; i < 18; i++) m_constantsGroupIndex << 3; // Measurement of Time m_constantsNames << i18n("Number of seconds in 1 minute"); m_constantsValues << QString::number(GSL_CONST_MKSA_MINUTE,'g',15); m_constantsUnits << "s"; m_constantsNames << i18n("Number of seconds in 1 hour"); m_constantsValues << QString::number(GSL_CONST_MKSA_HOUR,'g',15); m_constantsUnits << "s"; m_constantsNames << i18n("Number of seconds in 1 day"); m_constantsValues << QString::number(GSL_CONST_MKSA_DAY,'g',15); m_constantsUnits << "s"; m_constantsNames << i18n("Number of seconds in 1 week"); m_constantsValues << QString::number(GSL_CONST_MKSA_WEEK,'g',15); m_constantsUnits << "s"; for (int i = 0; i < 4; i++) m_constantsGroupIndex << 4; // Imperial Units m_constantsNames << i18n("Length of 1 inch"); m_constantsValues << QString::number(GSL_CONST_MKSA_INCH,'g',15); m_constantsUnits << "m"; m_constantsNames << i18n("Length of 1 foot"); m_constantsValues << QString::number(GSL_CONST_MKSA_FOOT,'g',15); m_constantsUnits << "m"; m_constantsNames << i18n("Length of 1 yard"); m_constantsValues << QString::number(GSL_CONST_MKSA_YARD,'g',15); m_constantsUnits << "m"; m_constantsNames << i18n("Length of 1 mile"); m_constantsValues << QString::number(GSL_CONST_MKSA_MILE,'g',15); m_constantsUnits << "m"; m_constantsNames << i18n("Length of 1/1000th of an inch"); m_constantsValues << QString::number(GSL_CONST_MKSA_MIL,'g',15); m_constantsUnits << "m"; for (int i = 0; i < 5; i++) m_constantsGroupIndex << 5; // Speed and Nautical Units m_constantsNames << i18n("Speed of 1 kilometer per hour"); m_constantsValues << QString::number(GSL_CONST_MKSA_KILOMETERS_PER_HOUR,'g',15); m_constantsUnits << "m / s"; m_constantsNames << i18n("Speed of 1 mile per hour"); m_constantsValues << QString::number(GSL_CONST_MKSA_MILES_PER_HOUR,'g',15); m_constantsUnits << "m / s"; m_constantsNames << i18n("Length of 1 nautical mile"); m_constantsValues << QString::number(GSL_CONST_MKSA_NAUTICAL_MILE,'g',15); m_constantsUnits << "m"; m_constantsNames << i18n("Length of 1 fathom"); m_constantsValues << QString::number(GSL_CONST_MKSA_FATHOM,'g',15); m_constantsUnits << "m"; m_constantsNames << i18n("Speed of 1 knot"); m_constantsValues << QString::number(GSL_CONST_MKSA_KNOT,'g',15); m_constantsUnits << "m / s"; for (int i = 0; i < 5; i++) m_constantsGroupIndex << 6; // Printers Units m_constantsNames << i18n("length of 1 printer's point [1/72 inch]"); m_constantsValues << QString::number(GSL_CONST_MKSA_POINT,'g',15); m_constantsUnits << "m"; m_constantsNames << i18n("length of 1 TeX point [1/72.27 inch]"); m_constantsValues << QString::number(GSL_CONST_MKSA_TEXPOINT,'g',15); m_constantsUnits << "m"; for (int i = 0; i < 2; i++) m_constantsGroupIndex << 7; // Volume, Area and Length m_constantsNames << i18n("Length of 1 micron"); m_constantsValues << QString::number(GSL_CONST_MKSA_MICRON,'g',15); m_constantsUnits << "m"; m_constantsNames << i18n("Area of 1 hectare"); m_constantsValues << QString::number(GSL_CONST_MKSA_HECTARE,'g',15); m_constantsUnits << "m^2"; m_constantsNames << i18n("Area of 1 acre"); m_constantsValues << QString::number(GSL_CONST_MKSA_ACRE,'g',15); m_constantsUnits << "m^2"; m_constantsNames << i18n("Volume of 1 liter"); m_constantsValues << QString::number(GSL_CONST_MKSA_LITER,'g',15); m_constantsUnits << "m^3"; m_constantsNames << i18n("Volume of 1 US gallon"); m_constantsValues << QString::number(GSL_CONST_MKSA_US_GALLON,'g',15); m_constantsUnits << "m^3"; m_constantsNames << i18n("Volume of 1 Canadian gallon"); m_constantsValues << QString::number(GSL_CONST_MKSA_CANADIAN_GALLON,'g',15); m_constantsUnits << "m^3"; m_constantsNames << i18n("Volume of 1 UK gallon"); m_constantsValues << QString::number(GSL_CONST_MKSA_UK_GALLON,'g',15); m_constantsUnits << "m^3"; m_constantsNames << i18n("Volume of 1 quart"); m_constantsValues << QString::number(GSL_CONST_MKSA_QUART,'g',15); m_constantsUnits << "m^3"; m_constantsNames << i18n("Volume of 1 pint"); m_constantsValues << QString::number(GSL_CONST_MKSA_PINT,'g',15); m_constantsUnits << "m^3"; for (int i = 0; i < 9; i++) m_constantsGroupIndex << 8; // Mass and Weight m_constantsNames << i18n("Mass of 1 pound"); m_constantsValues << QString::number(GSL_CONST_MKSA_POUND_MASS,'g',15); m_constantsUnits << "kg"; m_constantsNames << i18n("Mass of 1 ounce"); m_constantsValues << QString::number(GSL_CONST_MKSA_OUNCE_MASS,'g',15); m_constantsUnits << "kg"; m_constantsNames << i18n("Mass of 1 ton"); m_constantsValues << QString::number(GSL_CONST_MKSA_TON,'g',15); m_constantsUnits << "kg"; m_constantsNames << i18n("Mass of 1 metric ton [1000 kg]"); m_constantsValues << QString::number(GSL_CONST_MKSA_METRIC_TON,'g',15); m_constantsUnits << "kg"; m_constantsNames << i18n("Mass of 1 UK ton"); m_constantsValues << QString::number(GSL_CONST_MKSA_UK_TON,'g',15); m_constantsUnits << "kg"; m_constantsNames << i18n("Mass of 1 troy ounce"); m_constantsValues << QString::number(GSL_CONST_MKSA_TROY_OUNCE,'g',15); m_constantsUnits << "kg"; m_constantsNames << i18n("Mass of 1 carat"); m_constantsValues << QString::number(GSL_CONST_MKSA_CARAT,'g',15); m_constantsUnits << "kg"; m_constantsNames << i18n("Force of 1 gram weight"); m_constantsValues << QString::number(GSL_CONST_MKSA_GRAM_FORCE,'g',15); m_constantsUnits << "kg m / s^2"; m_constantsNames << i18n("Force of 1 pound weight"); m_constantsValues << QString::number(GSL_CONST_MKSA_POUND_FORCE,'g',15); m_constantsUnits << "kg m / s^2"; m_constantsNames << i18n("Force of 1 kilopound weight"); m_constantsValues << QString::number(GSL_CONST_MKSA_KILOPOUND_FORCE,'g',15); m_constantsUnits << "kg m / s^2"; m_constantsNames << i18n("Force of 1 poundal"); m_constantsValues << QString::number(GSL_CONST_MKSA_POUNDAL,'g',15); m_constantsUnits << "kg m / s^2"; for (int i = 0; i < 11; i++) m_constantsGroupIndex << 9; // Thermal Energy and Power m_constantsNames << i18n("Energy of 1 calorie"); m_constantsValues << QString::number(GSL_CONST_MKSA_CALORIE,'g',15); m_constantsUnits << "kg m^2 / s^2"; m_constantsNames << i18n("Energy of 1 British Thermal Unit"); m_constantsValues << QString::number(GSL_CONST_MKSA_BTU,'g',15); m_constantsUnits << "kg m^2 / s^2"; m_constantsNames << i18n("Energy of 1 Therm"); m_constantsValues << QString::number(GSL_CONST_MKSA_THERM,'g',15); m_constantsUnits << "kg m^2 / s^2"; m_constantsNames << i18n("Power of 1 horsepower"); m_constantsValues << QString::number(GSL_CONST_MKSA_HORSEPOWER,'g',15); m_constantsUnits << "kg m^2 / s^3"; for (int i = 0; i < 4; i++) m_constantsGroupIndex << 10; // Pressure m_constantsNames << i18n("Pressure of 1 bar"); m_constantsValues << QString::number(GSL_CONST_MKSA_BAR,'g',15); m_constantsUnits << "kg / m s^2"; m_constantsNames << i18n("Pressure of 1 standard atmosphere"); m_constantsValues << QString::number(GSL_CONST_MKSA_STD_ATMOSPHERE,'g',15); m_constantsUnits << "kg / m s^2"; m_constantsNames << i18n("Pressure of 1 torr"); m_constantsValues << QString::number(GSL_CONST_MKSA_TORR,'g',15); m_constantsUnits << "kg / m s^2"; m_constantsNames << i18n("Pressure of 1 meter of mercury"); m_constantsValues << QString::number(GSL_CONST_MKSA_METER_OF_MERCURY,'g',15); m_constantsUnits << "kg / m s^2"; m_constantsNames << i18n("Pressure of 1 inch of mercury"); m_constantsValues << QString::number(GSL_CONST_MKSA_INCH_OF_MERCURY,'g',15); m_constantsUnits << "kg / m s^2"; m_constantsNames << i18n("Pressure of 1 inch of water"); m_constantsValues << QString::number(GSL_CONST_MKSA_INCH_OF_WATER,'g',15); m_constantsUnits << "kg / m s^2"; m_constantsNames << i18n("Pressure of 1 pound per square inch"); m_constantsValues << QString::number(GSL_CONST_MKSA_PSI,'g',15); m_constantsUnits << "kg / m s^2"; for (int i = 0; i < 7; i++) m_constantsGroupIndex << 11; // Viscosity m_constantsNames << i18n("Dynamic viscosity of 1 poise"); m_constantsValues << QString::number(GSL_CONST_MKSA_POISE,'g',15); m_constantsUnits << "kg / m s"; m_constantsNames << i18n("Kinematic viscosity of 1 stokes"); m_constantsValues << QString::number(GSL_CONST_MKSA_STOKES,'g',15); m_constantsUnits << "m^2 / s"; for (int i = 0; i < 2; i++) m_constantsGroupIndex << 12; // Light and Illumination m_constantsNames << i18n("Luminance of 1 stilb"); m_constantsValues << QString::number(GSL_CONST_MKSA_STILB,'g',15); m_constantsUnits << "cd / m^2"; m_constantsNames << i18n("Luminous flux of 1 lumen"); m_constantsValues << QString::number(GSL_CONST_MKSA_LUMEN,'g',15); m_constantsUnits << "cd sr"; m_constantsNames << i18n("Illuminance of 1 lux"); m_constantsValues << QString::number(GSL_CONST_MKSA_LUX,'g',15); m_constantsUnits << "cd sr / m^2"; m_constantsNames << i18n("Illuminance of 1 phot"); m_constantsValues << QString::number(GSL_CONST_MKSA_PHOT,'g',15); m_constantsUnits << "cd sr / m^2"; m_constantsNames << i18n("Illuminance of 1 footcandle"); m_constantsValues << QString::number(GSL_CONST_MKSA_FOOTCANDLE,'g',15); m_constantsUnits << "cd sr / m^2"; m_constantsNames << i18n("Luminance of 1 lambert"); m_constantsValues << QString::number(GSL_CONST_MKSA_LAMBERT,'g',15); m_constantsUnits << "cd sr / m^2"; m_constantsNames << i18n("Luminance of 1 footlambert"); m_constantsValues << QString::number(GSL_CONST_MKSA_FOOTLAMBERT,'g',15); m_constantsUnits << "cd sr / m^2"; for (int i = 0; i < 7; i++) m_constantsGroupIndex << 13; // Radioactivity m_constantsNames << i18n("Activity of 1 curie"); m_constantsValues << QString::number(GSL_CONST_MKSA_CURIE,'g',15); m_constantsUnits << "1 / s"; m_constantsNames << i18n("Exposure of 1 roentgen"); m_constantsValues << QString::number(GSL_CONST_MKSA_ROENTGEN,'g',15); m_constantsUnits << "A s / kg"; m_constantsNames << i18n("Absorbed dose of 1 rad"); m_constantsValues << QString::number(GSL_CONST_MKSA_RAD,'g',15); m_constantsUnits << "m^2 / s^2"; for (int i = 0; i < 3; i++) m_constantsGroupIndex << 14; // Force and Energy m_constantsNames << i18n("SI unit of force"); m_constantsValues << QString::number(GSL_CONST_MKSA_NEWTON,'g',15); m_constantsUnits << "kg m / s^2"; m_constantsNames << i18n("Force of 1 Dyne"); m_constantsValues << QString::number(GSL_CONST_MKSA_DYNE,'g',15); m_constantsUnits << "kg m / s^2"; m_constantsNames << i18n("SI unit of energy"); m_constantsValues << QString::number(GSL_CONST_MKSA_JOULE,'g',15); m_constantsUnits << "kg m^2 / s^2"; m_constantsNames << i18n("Energy 1 erg"); m_constantsValues << QString::number(GSL_CONST_MKSA_ERG,'g',15); m_constantsUnits << "kg m^2 / s^2"; for (int i = 0; i < 4; i++) m_constantsGroupIndex << 15; } ExpressionParser::~ExpressionParser() { delete_table(); } ExpressionParser* ExpressionParser::getInstance() { if (!instance) instance = new ExpressionParser(); return instance; } const QStringList& ExpressionParser::functions() { return m_functions; } const QStringList& ExpressionParser::functionsGroups() { return m_functionsGroups; } const QStringList& ExpressionParser::functionsNames() { return m_functionsNames; } const QVector& ExpressionParser::functionsGroupIndices() { return m_functionsGroupIndex; } const QStringList& ExpressionParser::constants() { return m_constants; } const QStringList& ExpressionParser::constantsGroups() { return m_constantsGroups; } const QStringList& ExpressionParser::constantsNames() { return m_constantsNames; } const QStringList& ExpressionParser::constantsValues() { return m_constantsValues; } const QStringList& ExpressionParser::constantsUnits() { return m_constantsUnits; } const QVector& ExpressionParser::constantsGroupIndices() { return m_constantsGroupIndex; } bool ExpressionParser::isValid(const QString& expr, const QStringList& vars) { for (int i = 0; i < vars.size(); ++i) { QByteArray varba = vars.at(i).toLatin1(); assign_variable(varba.constData(), 0); } QByteArray funcba = expr.toLatin1(); const char* data = funcba.constData(); gsl_set_error_handler_off(); parse(data); return !(parse_errors() > 0); } QStringList ExpressionParser::getParameter(const QString& expr, const QStringList& vars) { DEBUG("ExpressionParser::getParameter()"); QDEBUG("variables:" << vars); QStringList parameters; QStringList strings = expr.split(QRegExp("\\W+"), QString::SkipEmptyParts); QDEBUG("found strings:" << strings); for (int i = 0; i < strings.size(); ++i) { // check if token is not a known constant/function/variable or number if (constants().indexOf(strings[i]) == -1 && functions().indexOf(strings[i]) == -1 && vars.indexOf(strings[i]) == -1 && QRegExp("[0-9]*").exactMatch(strings[i]) == 0) parameters << strings[i]; else QDEBUG(strings[i] << ':' << constants().indexOf(strings[i]) << ' ' << functions().indexOf(strings[i]) << ' ' << vars.indexOf(strings[i]) << ' ' << QRegExp("[0-9]*").exactMatch(strings[i])); } parameters.removeDuplicates(); QDEBUG("parameters found:" << parameters); return parameters; } bool ExpressionParser::evaluateCartesian(const QString& expr, const QString& min, const QString& max, int count, QVector* xVector, QVector* yVector, const QStringList& paramNames, const QVector& paramValues) { QByteArray xminba = min.toLatin1(); const double xMin = parse(xminba.constData()); QByteArray xmaxba = max.toLatin1(); const double xMax = parse(xmaxba.constData()); const double step = (xMax - xMin)/(double)(count - 1); QByteArray funcba = expr.toLatin1(); const char* func = funcba.constData(); double x, y; gsl_set_error_handler_off(); for (int i = 0; i < paramNames.size(); ++i) { QByteArray paramba = paramNames.at(i).toLatin1(); assign_variable(paramba.constData(), paramValues.at(i)); } for (int i = 0; i < count; i++) { x = xMin + step * i; assign_variable("x", x); y = parse(func); if (parse_errors() > 0) return false; (*xVector)[i] = x; if (std::isfinite(y)) (*yVector)[i] = y; else (*yVector)[i] = NAN; } return true; } bool ExpressionParser::evaluateCartesian(const QString& expr, const QString& min, const QString& max, int count, QVector* xVector, QVector* yVector) { QByteArray xminba = min.toLatin1(); const double xMin = parse(xminba.constData()); QByteArray xmaxba = max.toLatin1(); const double xMax = parse(xmaxba.constData()); const double step = (xMax - xMin)/(double)(count - 1); QByteArray funcba = expr.toLatin1(); const char* func = funcba.constData(); double x, y; gsl_set_error_handler_off(); for (int i = 0; i < count; i++) { x = xMin + step*i; assign_variable("x", x); y = parse(func); if (parse_errors() > 0) return false; (*xVector)[i] = x; if (std::isfinite(y)) (*yVector)[i] = y; else (*yVector)[i] = NAN; } return true; } bool ExpressionParser::evaluateCartesian(const QString& expr, QVector* xVector, QVector* yVector) { QByteArray funcba = expr.toLatin1(); const char* func = funcba.constData(); double x, y; gsl_set_error_handler_off(); for (int i = 0; i < xVector->count(); i++) { x = xVector->at(i); assign_variable("x", x); y = parse(func); if (parse_errors() > 0) return false; if (std::isfinite(y)) (*yVector)[i] = y; else (*yVector)[i] = NAN; } return true; } bool ExpressionParser::evaluateCartesian(const QString& expr, QVector* xVector, QVector* yVector, const QStringList& paramNames, const QVector& paramValues) { QByteArray funcba = expr.toLatin1(); const char* func = funcba.constData(); double x, y; gsl_set_error_handler_off(); for (int i = 0; i < paramNames.size(); ++i) { QByteArray paramba = paramNames.at(i).toLatin1(); assign_variable(paramba.constData(), paramValues.at(i)); } for (int i = 0; i < xVector->count(); i++) { x = xVector->at(i); assign_variable("x", x); y = parse(func); if (parse_errors() > 0) return false; if (std::isfinite(y)) (*yVector)[i] = y; else (*yVector)[i] = NAN; } return true; } /*! evaluates multivariate function y=f(x_1, x_2, ...). Variable names (x_1, x_2, ...) are stored in \c vars. Data is stored in \c dataVectors. */ bool ExpressionParser::evaluateCartesian(const QString& expr, const QStringList& vars, const QVector*>& xVectors, QVector* yVector) { Q_ASSERT(vars.size() == xVectors.size()); QByteArray funcba = expr.toLatin1(); const char* func = funcba.constData(); double y, varValue; QString varName; gsl_set_error_handler_off(); bool stop = false; for (int i = 0; i < yVector->size(); i++) { //stop iterating over i if one of the x-vectors has no elements anymore. for (int n = 0; n < xVectors.size(); ++n) { if (i == xVectors.at(n)->size()) { stop = true; break; } } if (stop) break; for (int n = 0; n < vars.size(); ++n) { varName = vars.at(n); varValue = xVectors.at(n)->at(i); QByteArray varba = varName.toLatin1(); assign_variable(varba.constData(), varValue); } y = parse(func); if (parse_errors() > 0) return false; if (std::isfinite(y)) (*yVector)[i] = y; else (*yVector)[i] = NAN; } return true; } bool ExpressionParser::evaluatePolar(const QString& expr, const QString& min, const QString& max, int count, QVector* xVector, QVector* yVector) { QByteArray minba = min.toLatin1(); const double minValue = parse(minba.constData()); QByteArray maxba = max.toLatin1(); const double maxValue = parse(maxba.constData()); const double step = (maxValue - minValue)/(double)(count - 1); QByteArray funcba = expr.toLatin1(); const char* func = funcba.constData(); double r, phi; gsl_set_error_handler_off(); for (int i = 0; i < count; i++) { phi = minValue + step * i; assign_variable("phi", phi); r = parse(func); if (parse_errors() > 0) return false; if (std::isfinite(r)) { (*xVector)[i] = r*cos(phi); (*yVector)[i] = r*sin(phi); } else { (*xVector)[i] = NAN; (*yVector)[i] = NAN; } } return true; } bool ExpressionParser::evaluateParametric(const QString& expr1, const QString& expr2, const QString& min, const QString& max, int count, QVector* xVector, QVector* yVector) { QByteArray minba = min.toLatin1(); const double minValue = parse(minba.constData()); QByteArray maxba = max.toLatin1(); const double maxValue = parse(maxba.constData()); const double step = (maxValue - minValue)/(double)(count - 1); QByteArray xfuncba = expr1.toLatin1(); const char* xFunc = xfuncba.constData(); QByteArray yfuncba = expr2.toLatin1(); const char* yFunc = yfuncba.constData(); double x, y, t; gsl_set_error_handler_off(); for (int i = 0; i < count; i++) { t = minValue + step*i; assign_variable("t", t); x = parse(xFunc); if (parse_errors() > 0) return false; if (std::isfinite(x)) (*xVector)[i] = x; else (*xVector)[i] = NAN; y = parse(yFunc); if (parse_errors() > 0) return false; if (std::isfinite(y)) (*yVector)[i] = y; else (*yVector)[i] = NAN; } return true; } diff --git a/src/backend/gsl/functions.h b/src/backend/gsl/functions.h index eafa330b2..83e59834a 100644 --- a/src/backend/gsl/functions.h +++ b/src/backend/gsl/functions.h @@ -1,491 +1,492 @@ /*************************************************************************** File : functions.h Project : LabPlot Description : definition of functions -------------------------------------------------------------------- Copyright : (C) 2014 by Alexander Semke (alexander.semke@web.de) Copyright : (C) 2014-2018 by Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #ifndef FUNCTIONS_H #define FUNCTIONS_H #include #include #include #include #include #include "backend/nsl/nsl_sf_basic.h" #ifdef _MSC_VER /* avoid intrinsics */ #pragma function(ceil, floor) #endif /* sync with ExpressionParser.cpp */ struct func _functions[] = { /* standard functions */ {"rand", (func_t)nsl_sf_rand}, {"random", (func_t)nsl_sf_random}, {"drand", (func_t)nsl_sf_drand}, /* math.h */ {"ceil", (func_t)ceil}, {"fabs", (func_t)fabs}, {"log10", (func_t)log10}, {"pow", (func_t)pow}, {"sqrt", (func_t)sqrt}, {"sgn", (func_t)nsl_sf_sgn}, {"theta", (func_t)nsl_sf_theta}, {"harmonic", (func_t)nsl_sf_harmonic}, #ifndef _WIN32 {"cbrt", (func_t)cbrt}, {"logb", (func_t)logb}, {"rint", (func_t)rint}, {"round", (func_t)round}, {"trunc", (func_t)trunc}, #endif /* GSL mathematical functions: see http://www.gnu.org/software/gsl/manual/gsl-ref.html#Mathematical-Functions */ {"log1p", (func_t)gsl_log1p}, {"ldexp", (func_t)nsl_sf_ldexp}, {"powint", (func_t)nsl_sf_powint}, {"pow2", (func_t)gsl_pow_2}, {"pow3", (func_t)gsl_pow_3}, {"pow4", (func_t)gsl_pow_4}, {"pow5", (func_t)gsl_pow_5}, {"pow6", (func_t)gsl_pow_6}, {"pow7", (func_t)gsl_pow_7}, {"pow8", (func_t)gsl_pow_8}, {"pow9", (func_t)gsl_pow_9}, /* GSL special functions: see http://www.gnu.org/software/gsl/manual/html_node/Special-Functions.html */ /* Airy Functions and Derivatives */ {"Ai", (func_t)nsl_sf_airy_Ai}, {"Bi", (func_t)nsl_sf_airy_Bi}, {"Ais", (func_t)nsl_sf_airy_Ais}, {"Bis", (func_t)nsl_sf_airy_Bis}, {"Aid", (func_t)nsl_sf_airy_Aid}, {"Bid", (func_t)nsl_sf_airy_Bid}, {"Aids", (func_t)nsl_sf_airy_Aids}, {"Bids", (func_t)nsl_sf_airy_Bids}, {"Ai0", (func_t)nsl_sf_airy_0_Ai}, {"Bi0", (func_t)nsl_sf_airy_0_Bi}, {"Aid0", (func_t)nsl_sf_airy_0_Aid}, {"Bid0", (func_t)nsl_sf_airy_0_Bid}, /* Bessel Functions */ {"J0", (func_t)gsl_sf_bessel_J0}, {"J1", (func_t)gsl_sf_bessel_J1}, {"Jn", (func_t)nsl_sf_bessel_Jn}, {"Y0", (func_t)gsl_sf_bessel_Y0}, {"Y1", (func_t)gsl_sf_bessel_Y1}, {"Yn", (func_t)nsl_sf_bessel_Yn}, {"I0", (func_t)gsl_sf_bessel_I0}, {"I1", (func_t)gsl_sf_bessel_I1}, {"In", (func_t)nsl_sf_bessel_In}, {"I0s", (func_t)gsl_sf_bessel_I0_scaled}, {"I1s", (func_t)gsl_sf_bessel_I1_scaled}, {"Ins", (func_t)nsl_sf_bessel_Ins}, {"K0", (func_t)gsl_sf_bessel_K0}, {"K1", (func_t)gsl_sf_bessel_K1}, {"Kn", (func_t)nsl_sf_bessel_Kn}, {"K0s", (func_t)gsl_sf_bessel_K0_scaled}, {"K1s", (func_t)gsl_sf_bessel_K1_scaled}, {"Kns", (func_t)nsl_sf_bessel_Kns}, {"j0", (func_t)gsl_sf_bessel_j0}, {"j1", (func_t)gsl_sf_bessel_j1}, {"j2", (func_t)gsl_sf_bessel_j2}, {"jl", (func_t)nsl_sf_bessel_jl}, {"y0", (func_t)gsl_sf_bessel_y0}, {"y1", (func_t)gsl_sf_bessel_y1}, {"y2", (func_t)gsl_sf_bessel_y2}, {"yl", (func_t)nsl_sf_bessel_yl}, {"i0s", (func_t)gsl_sf_bessel_i0_scaled}, {"i1s", (func_t)gsl_sf_bessel_i1_scaled}, {"i2s", (func_t)gsl_sf_bessel_i2_scaled}, {"ils", (func_t)nsl_sf_bessel_ils}, {"k0s", (func_t)gsl_sf_bessel_k0_scaled}, {"k1s", (func_t)gsl_sf_bessel_k1_scaled}, {"k2s", (func_t)gsl_sf_bessel_k2_scaled}, {"kls", (func_t)nsl_sf_bessel_kls}, {"Jnu", (func_t)gsl_sf_bessel_Jnu}, {"Ynu", (func_t)gsl_sf_bessel_Ynu}, {"Inu", (func_t)gsl_sf_bessel_Inu}, {"Inus", (func_t)gsl_sf_bessel_Inu_scaled}, {"Knu", (func_t)gsl_sf_bessel_Knu}, {"lnKnu", (func_t)gsl_sf_bessel_lnKnu}, {"Knus", (func_t)gsl_sf_bessel_Knu_scaled}, {"J0_0", (func_t)nsl_sf_bessel_0_J0}, {"J1_0", (func_t)nsl_sf_bessel_0_J1}, {"Jnu_0", (func_t)nsl_sf_bessel_0_Jnu}, /* Clausen functions */ {"clausen", (func_t)gsl_sf_clausen}, /* Coulomb functions */ {"hydrogenicR_1", (func_t)gsl_sf_hydrogenicR_1}, {"hydrogenicR", (func_t)nsl_sf_hydrogenicR}, {"dawson", (func_t)gsl_sf_dawson}, /* Debye functions */ {"D1", (func_t)gsl_sf_debye_1}, {"D2", (func_t)gsl_sf_debye_2}, {"D3", (func_t)gsl_sf_debye_3}, {"D4", (func_t)gsl_sf_debye_4}, {"D5", (func_t)gsl_sf_debye_5}, {"D6", (func_t)gsl_sf_debye_6}, {"Li2", (func_t)gsl_sf_dilog}, /* Elliptic integrals */ {"Kc", (func_t)nsl_sf_ellint_Kc}, {"Ec", (func_t)nsl_sf_ellint_Ec}, {"Pc", (func_t)nsl_sf_ellint_Pc}, {"F", (func_t)nsl_sf_ellint_F}, {"E", (func_t)nsl_sf_ellint_E}, {"P", (func_t)nsl_sf_ellint_P}, {"D", (func_t)nsl_sf_ellint_D}, {"RC", (func_t)nsl_sf_ellint_RC}, {"RD", (func_t)nsl_sf_ellint_RD}, {"RF", (func_t)nsl_sf_ellint_RF}, {"RJ", (func_t)nsl_sf_ellint_RJ}, /* Error functions */ {"erf", (func_t)gsl_sf_erf}, {"erfc", (func_t)gsl_sf_erfc}, {"log_erfc", (func_t)gsl_sf_log_erfc}, {"erf_Z", (func_t)gsl_sf_erf_Z}, {"erf_Q", (func_t)gsl_sf_erf_Q}, {"hazard", (func_t)gsl_sf_hazard}, #ifndef _MSC_VER {"erfcx", (func_t)nsl_sf_erfcx}, {"erfi", (func_t)nsl_sf_erfi}, {"im_w_of_x", (func_t)nsl_sf_im_w_of_x}, {"dawson", (func_t)nsl_sf_dawson}, {"voigt", (func_t)nsl_sf_voigt}, #endif + {"pseudovoigt1", (func_t)nsl_sf_pseudovoigt1}, /* Exponential Functions */ {"exp", (func_t)gsl_sf_exp}, {"exp_mult", (func_t)gsl_sf_exp_mult}, {"expm1", (func_t)gsl_expm1}, {"exprel", (func_t)gsl_sf_exprel}, {"exprel2", (func_t)gsl_sf_exprel_2}, {"expreln", (func_t)nsl_sf_exprel_n}, /* Exponential Integrals */ {"E1", (func_t)gsl_sf_expint_E1}, {"E2", (func_t)gsl_sf_expint_E2}, {"En", (func_t)gsl_sf_expint_En}, {"Ei", (func_t)gsl_sf_expint_Ei}, {"Shi", (func_t)gsl_sf_Shi}, {"Chi", (func_t)gsl_sf_Chi}, {"Ei3", (func_t)gsl_sf_expint_3}, {"Si", (func_t)gsl_sf_Si}, {"Ci", (func_t)gsl_sf_Ci}, {"Atanint", (func_t)gsl_sf_atanint}, /* Fermi-Dirac Function */ {"Fm1", (func_t)gsl_sf_fermi_dirac_m1}, {"F0", (func_t)gsl_sf_fermi_dirac_0}, {"F1", (func_t)gsl_sf_fermi_dirac_1}, {"F2", (func_t)gsl_sf_fermi_dirac_2}, {"Fj", (func_t)nsl_sf_fermi_dirac_int}, {"Fmhalf", (func_t)gsl_sf_fermi_dirac_mhalf}, {"Fhalf", (func_t)gsl_sf_fermi_dirac_half}, {"F3half", (func_t)gsl_sf_fermi_dirac_3half}, {"Finc0", (func_t)gsl_sf_fermi_dirac_inc_0}, /* Gamma and Beta Functions */ {"gamma", (func_t)gsl_sf_gamma}, {"tgamma", (func_t)gsl_sf_gamma}, {"lgamma", (func_t)gsl_sf_lngamma}, {"lngamma", (func_t)gsl_sf_lngamma}, {"gammastar", (func_t)gsl_sf_gammastar}, {"gammainv", (func_t)gsl_sf_gammainv}, {"fact", (func_t)nsl_sf_fact}, {"doublefact", (func_t)nsl_sf_doublefact}, {"lnfact", (func_t)nsl_sf_lnfact}, {"lndoublefact", (func_t)nsl_sf_lndoublefact}, {"choose", (func_t)nsl_sf_choose}, {"lnchoose", (func_t)nsl_sf_lnchoose}, {"taylor", (func_t)nsl_sf_taylorcoeff}, {"poch", (func_t)gsl_sf_poch}, {"lnpoch", (func_t)gsl_sf_lnpoch}, {"pochrel", (func_t)gsl_sf_pochrel}, {"gammainc", (func_t)gsl_sf_gamma_inc}, {"gammaincQ", (func_t)gsl_sf_gamma_inc_Q}, {"gammaincP", (func_t)gsl_sf_gamma_inc_P}, {"beta", (func_t)gsl_sf_beta}, {"lnbeta", (func_t)gsl_sf_lnbeta}, {"betainc", (func_t)gsl_sf_beta_inc}, /* Gegenbauer Functions */ {"C1", (func_t)gsl_sf_gegenpoly_1}, {"C2", (func_t)gsl_sf_gegenpoly_2}, {"C3", (func_t)gsl_sf_gegenpoly_3}, {"Cn", (func_t)nsl_sf_gegenpoly_n}, #if (GSL_MAJOR_VERSION > 2) || (GSL_MAJOR_VERSION == 2) && (GSL_MINOR_VERSION >= 4) /* Hermite polynomials and functions */ {"Hen", (func_t)nsl_sf_hermite_prob}, {"Hn", (func_t)nsl_sf_hermite_phys}, {"Hfn", (func_t)nsl_sf_hermite_func}, {"Hend", (func_t)nsl_sf_hermite_prob_der}, {"Hnd", (func_t)nsl_sf_hermite_phys_der}, {"Hfnd", (func_t)nsl_sf_hermite_func_der}, #endif /* Hypergeometric Functions */ {"hyperg_0F1", (func_t)gsl_sf_hyperg_0F1}, {"hyperg_1F1i", (func_t)nsl_sf_hyperg_1F1i}, {"hyperg_1F1", (func_t)gsl_sf_hyperg_1F1}, {"hyperg_Ui", (func_t)nsl_sf_hyperg_Ui}, {"hyperg_U", (func_t)gsl_sf_hyperg_U}, {"hyperg_2F1", (func_t)gsl_sf_hyperg_2F1}, {"hyperg_2F1c", (func_t)gsl_sf_hyperg_2F1_conj}, {"hyperg_2F1r", (func_t)gsl_sf_hyperg_2F1_renorm}, {"hyperg_2F1cr", (func_t)gsl_sf_hyperg_2F1_conj_renorm}, {"hyperg_2F0", (func_t)gsl_sf_hyperg_2F0}, /* Laguerre Functions */ {"L1", (func_t)gsl_sf_laguerre_1}, {"L2", (func_t)gsl_sf_laguerre_2}, {"L3", (func_t)gsl_sf_laguerre_3}, /* Lambert W Functions */ {"W0", (func_t)gsl_sf_lambert_W0}, {"Wm1", (func_t)gsl_sf_lambert_Wm1}, /* Legendre Functions and Spherical Harmonics */ {"P1", (func_t)gsl_sf_legendre_P1}, {"P2", (func_t)gsl_sf_legendre_P2}, {"P3", (func_t)gsl_sf_legendre_P3}, {"Pl", (func_t)nsl_sf_legendre_Pl}, {"Q0", (func_t)gsl_sf_legendre_Q0}, {"Q1", (func_t)gsl_sf_legendre_Q1}, {"Ql", (func_t)nsl_sf_legendre_Ql}, {"Plm", (func_t)nsl_sf_legendre_Plm}, {"Pslm", (func_t)nsl_sf_legendre_sphPlm}, {"Phalf", (func_t)gsl_sf_conicalP_half}, {"Pmhalf", (func_t)gsl_sf_conicalP_mhalf}, {"Pc0", (func_t)gsl_sf_conicalP_0}, {"Pc1", (func_t)gsl_sf_conicalP_1}, {"Psr", (func_t)nsl_sf_conicalP_sphreg}, {"Pcr", (func_t)nsl_sf_conicalP_cylreg}, {"H3d0", (func_t)gsl_sf_legendre_H3d_0}, {"H3d1", (func_t)gsl_sf_legendre_H3d_1}, {"H3d", (func_t)nsl_sf_legendre_H3d}, /* Logarithm and Related Functions */ {"log", (func_t)gsl_sf_log}, {"logabs", (func_t)gsl_sf_log_abs}, {"logp", (func_t)gsl_sf_log_1plusx}, {"logpm", (func_t)gsl_sf_log_1plusx_mx}, /* Power Function */ {"gsl_powint", (func_t)nsl_sf_powint}, /* Psi (Digamma) Function */ {"psiint", (func_t)nsl_sf_psiint}, {"psi", (func_t)gsl_sf_psi}, {"psi1piy", (func_t)gsl_sf_psi_1piy}, {"psi1int", (func_t)nsl_sf_psi1int}, {"psi1", (func_t)gsl_sf_psi_1}, {"psin", (func_t)nsl_sf_psin}, /* Synchrotron Functions */ {"synchrotron1", (func_t)gsl_sf_synchrotron_1}, {"synchrotron2", (func_t)gsl_sf_synchrotron_2}, /* Transport Functions */ {"J2", (func_t)gsl_sf_transport_2}, {"J3", (func_t)gsl_sf_transport_3}, {"J4", (func_t)gsl_sf_transport_4}, {"J5", (func_t)gsl_sf_transport_5}, /* Trigonometric Functions */ {"sin", (func_t)gsl_sf_sin}, {"cos", (func_t)gsl_sf_cos}, {"tan", (func_t)tan}, {"asin", (func_t)asin}, {"acos", (func_t)acos}, {"atan", (func_t)atan}, {"atan2", (func_t)atan2}, {"sinh", (func_t)sinh}, {"cosh", (func_t)cosh}, {"tanh", (func_t)tanh}, {"acosh", (func_t)gsl_acosh}, {"asinh", (func_t)gsl_asinh}, {"atanh", (func_t)gsl_atanh}, {"sec", (func_t)nsl_sf_sec}, {"csc", (func_t)nsl_sf_csc}, {"cot", (func_t)nsl_sf_cot}, {"asec", (func_t)nsl_sf_asec}, {"acsc", (func_t)nsl_sf_acsc}, {"acot", (func_t)nsl_sf_acot}, {"sech", (func_t)nsl_sf_sech}, {"csch", (func_t)nsl_sf_csch}, {"coth", (func_t)nsl_sf_coth}, {"asech", (func_t)nsl_sf_asech}, {"acsch", (func_t)nsl_sf_acsch}, {"acoth", (func_t)nsl_sf_acoth}, {"sinc", (func_t)gsl_sf_sinc}, {"logsinh", (func_t)gsl_sf_lnsinh}, {"logcosh", (func_t)gsl_sf_lncosh}, {"hypot", (func_t)gsl_sf_hypot}, {"hypot3", (func_t)gsl_hypot3}, {"anglesymm", (func_t)gsl_sf_angle_restrict_symm}, {"anglepos", (func_t)gsl_sf_angle_restrict_pos}, /* Zeta Functions */ {"zetaint", (func_t)nsl_sf_zetaint}, {"zeta", (func_t)gsl_sf_zeta}, {"zetam1int", (func_t)nsl_sf_zetam1int}, {"zetam1", (func_t)gsl_sf_zetam1}, {"hzeta", (func_t)gsl_sf_hzeta}, {"etaint", (func_t)nsl_sf_etaint}, {"eta", (func_t)gsl_sf_eta}, /* GSL Random Number Distributions: see http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Distributions.html */ /* Gaussian Distribution */ {"gaussian", (func_t)gsl_ran_gaussian_pdf}, {"ugaussian", (func_t)gsl_ran_ugaussian_pdf}, {"gaussianP", (func_t)gsl_cdf_gaussian_P}, {"gaussianQ", (func_t)gsl_cdf_gaussian_Q}, {"gaussianPinv", (func_t)gsl_cdf_gaussian_Pinv}, {"gaussianQinv", (func_t)gsl_cdf_gaussian_Qinv}, {"ugaussianP", (func_t)gsl_cdf_ugaussian_P}, {"ugaussianQ", (func_t)gsl_cdf_ugaussian_Q}, {"ugaussianPinv", (func_t)gsl_cdf_ugaussian_Pinv}, {"ugaussianQinv", (func_t)gsl_cdf_ugaussian_Qinv}, {"gaussiantail", (func_t)gsl_ran_gaussian_tail_pdf}, {"ugaussiantail", (func_t)gsl_ran_ugaussian_tail_pdf}, {"gaussianbi", (func_t)gsl_ran_bivariate_gaussian_pdf}, /* Exponential Distribution */ {"exponential", (func_t)gsl_ran_exponential_pdf}, {"exponentialP", (func_t)gsl_cdf_exponential_P}, {"exponentialQ", (func_t)gsl_cdf_exponential_Q}, {"exponentialPinv", (func_t)gsl_cdf_exponential_Pinv}, {"exponentialQinv", (func_t)gsl_cdf_exponential_Qinv}, /* Laplace Distribution */ {"laplace", (func_t)gsl_ran_laplace_pdf}, {"laplaceP", (func_t)gsl_cdf_laplace_P}, {"laplaceQ", (func_t)gsl_cdf_laplace_Q}, {"laplacePinv", (func_t)gsl_cdf_laplace_Pinv}, {"laplaceQinv", (func_t)gsl_cdf_laplace_Qinv}, /* Exponential Power Distribution */ {"exppow", (func_t)gsl_ran_exppow_pdf}, {"exppowP", (func_t)gsl_cdf_exppow_P}, {"exppowQ", (func_t)gsl_cdf_exppow_Q}, /* Cauchy Distribution */ {"cauchy", (func_t)gsl_ran_cauchy_pdf}, {"cauchyP", (func_t)gsl_cdf_cauchy_P}, {"cauchyQ", (func_t)gsl_cdf_cauchy_Q}, {"cauchyPinv", (func_t)gsl_cdf_cauchy_Pinv}, {"cauchyQinv", (func_t)gsl_cdf_cauchy_Qinv}, /* Rayleigh Distribution */ {"rayleigh", (func_t)gsl_ran_rayleigh_pdf}, {"rayleighP", (func_t)gsl_cdf_rayleigh_P}, {"rayleighQ", (func_t)gsl_cdf_rayleigh_Q}, {"rayleighPinv", (func_t)gsl_cdf_rayleigh_Pinv}, {"rayleighQinv", (func_t)gsl_cdf_rayleigh_Qinv}, {"rayleigh_tail", (func_t)gsl_ran_rayleigh_tail_pdf}, /* Landau Distribution */ {"landau", (func_t)gsl_ran_landau_pdf}, /* Gamma Distribution */ {"gammapdf", (func_t)gsl_ran_gamma_pdf}, {"gammaP", (func_t)gsl_cdf_gamma_P}, {"gammaQ", (func_t)gsl_cdf_gamma_Q}, {"gammaPinv", (func_t)gsl_cdf_gamma_Pinv}, {"gammaQinv", (func_t)gsl_cdf_gamma_Qinv}, /* Flat (Uniform) Distribution */ {"flat", (func_t)gsl_ran_flat_pdf}, {"flatP", (func_t)gsl_cdf_flat_P}, {"flatQ", (func_t)gsl_cdf_flat_Q}, {"flatPinv", (func_t)gsl_cdf_flat_Pinv}, {"flatQinv", (func_t)gsl_cdf_flat_Qinv}, /* Lognormal Distribution */ {"lognormal", (func_t)gsl_ran_lognormal_pdf}, {"lognormalP", (func_t)gsl_cdf_lognormal_P}, {"lognormalQ", (func_t)gsl_cdf_lognormal_Q}, {"lognormalPinv", (func_t)gsl_cdf_lognormal_Pinv}, {"lognormalQinv", (func_t)gsl_cdf_lognormal_Qinv}, /* Chi-squared Distribution */ {"chisq", (func_t)gsl_ran_chisq_pdf}, {"chisqP", (func_t)gsl_cdf_chisq_P}, {"chisqQ", (func_t)gsl_cdf_chisq_Q}, {"chisqPinv", (func_t)gsl_cdf_chisq_Pinv}, {"chisqQinv", (func_t)gsl_cdf_chisq_Qinv}, /* F-distribution */ {"fdist", (func_t)gsl_ran_fdist_pdf}, {"fdistP", (func_t)gsl_cdf_fdist_P}, {"fdistQ", (func_t)gsl_cdf_fdist_Q}, {"fdistPinv", (func_t)gsl_cdf_fdist_Pinv}, {"fdistQinv", (func_t)gsl_cdf_fdist_Qinv}, /* t-distribution */ {"tdist", (func_t)gsl_ran_tdist_pdf}, {"tdistP", (func_t)gsl_cdf_tdist_P}, {"tdistQ", (func_t)gsl_cdf_tdist_Q}, {"tdistPinv", (func_t)gsl_cdf_tdist_Pinv}, {"tdistQinv", (func_t)gsl_cdf_tdist_Qinv}, /* Beta Distribution */ {"betapdf", (func_t)gsl_ran_beta_pdf}, {"betaP", (func_t)gsl_cdf_beta_P}, {"betaQ", (func_t)gsl_cdf_beta_Q}, {"betaPinv", (func_t)gsl_cdf_beta_Pinv}, {"betaQinv", (func_t)gsl_cdf_beta_Qinv}, /* Logistic Distribution */ {"logistic", (func_t)gsl_ran_logistic_pdf}, {"logisticP", (func_t)gsl_cdf_logistic_P}, {"logisticQ", (func_t)gsl_cdf_logistic_Q}, {"logisticPinv", (func_t)gsl_cdf_logistic_Pinv}, {"logisticQinv", (func_t)gsl_cdf_logistic_Qinv}, /* Pareto Distribution */ {"pareto", (func_t)gsl_ran_pareto_pdf}, {"paretoP", (func_t)gsl_cdf_pareto_P}, {"paretoQ", (func_t)gsl_cdf_pareto_Q}, {"paretoPinv", (func_t)gsl_cdf_pareto_Pinv}, {"paretoQinv", (func_t)gsl_cdf_pareto_Qinv}, /* Weibull Distribution */ {"weibull", (func_t)gsl_ran_weibull_pdf}, {"weibullP", (func_t)gsl_cdf_weibull_P}, {"weibullQ", (func_t)gsl_cdf_weibull_Q}, {"weibullPinv", (func_t)gsl_cdf_weibull_Pinv}, {"weibullQinv", (func_t)gsl_cdf_weibull_Qinv}, /* Gumbel Distribution */ {"gumbel1", (func_t)gsl_ran_gumbel1_pdf}, {"gumbel1P", (func_t)gsl_cdf_gumbel1_P}, {"gumbel1Q", (func_t)gsl_cdf_gumbel1_Q}, {"gumbel1Pinv", (func_t)gsl_cdf_gumbel1_Pinv}, {"gumbel1Qinv", (func_t)gsl_cdf_gumbel1_Qinv}, {"gumbel2", (func_t)gsl_ran_gumbel2_pdf}, {"gumbel2P", (func_t)gsl_cdf_gumbel2_P}, {"gumbel2Q", (func_t)gsl_cdf_gumbel2_Q}, {"gumbel2Pinv", (func_t)gsl_cdf_gumbel2_Pinv}, {"gumbel2Qinv", (func_t)gsl_cdf_gumbel2_Qinv}, /* Poisson Distribution */ {"poisson", (func_t)nsl_sf_poisson}, {"poissonP", (func_t)gsl_cdf_poisson_P}, {"poissonQ", (func_t)gsl_cdf_poisson_Q}, /* Bernoulli Distribution */ {"bernoulli", (func_t)nsl_sf_bernoulli}, /* Binomial Distribution */ {"binomial", (func_t)nsl_sf_binomial}, {"binomialP", (func_t)gsl_cdf_binomial_P}, {"binomialQ", (func_t)gsl_cdf_binomial_Q}, {"negative_binomial", (func_t)nsl_sf_negative_binomial}, {"negative_binomialP", (func_t)gsl_cdf_negative_binomial_P}, {"negative_binomialQ", (func_t)gsl_cdf_negative_binomial_Q}, /* Pascal Distribution */ {"pascal", (func_t)nsl_sf_pascal}, {"pascalP", (func_t)gsl_cdf_pascal_P}, {"pascalQ", (func_t)gsl_cdf_pascal_Q}, /* Geometric Distribution */ {"geometric", (func_t)nsl_sf_geometric}, {"geometricP", (func_t)gsl_cdf_geometric_P}, {"geometricQ", (func_t)gsl_cdf_geometric_Q}, /* Hypergeometric Distribution */ {"hypergeometric", (func_t)nsl_sf_hypergeometric}, {"hypergeometricP", (func_t)gsl_cdf_hypergeometric_P}, {"hypergeometricQ", (func_t)gsl_cdf_hypergeometric_Q}, /* Logarithmic Distribution */ {"logarithmic", (func_t)nsl_sf_logarithmic}, {0, (func_t)0} }; #endif /*FUNCTIONS_H*/ diff --git a/src/backend/nsl/nsl_fit.c b/src/backend/nsl/nsl_fit.c index 5b4402976..6f4600e12 100644 --- a/src/backend/nsl/nsl_fit.c +++ b/src/backend/nsl/nsl_fit.c @@ -1,686 +1,686 @@ /*************************************************************************** File : nsl_fit.c Project : LabPlot Description : NSL (non)linear fit functions -------------------------------------------------------------------- Copyright : (C) 2016-2018 by Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #include "nsl_fit.h" #include "nsl_common.h" #include "nsl_sf_basic.h" #include #include #include #include #include const char* nsl_fit_model_category_name[] = {i18n("Basic functions"), i18n("Peak functions"), i18n("Growth (sigmoidal)"), i18n("Statistics (distributions)"), i18n("Custom")}; const char* nsl_fit_model_basic_name[] = {i18n("Polynomial"), i18n("Power"), i18n("Exponential"), i18n("Inverse exponential"), i18n("Fourier")}; const char* nsl_fit_model_basic_equation[] = {"c0 + c1*x", "a*x^b", "a*exp(b*x)", "a*(1-exp(b*x)) + c", "a0 + (a1*cos(w*x) + b1*sin(w*x))"}; const char* nsl_fit_model_basic_pic_name[] = {"polynom", "power", "exponential", "inv_exponential", "fourier"}; const char* nsl_fit_model_peak_name[] = {i18n("Gaussian (normal)"), i18n("Cauchy-Lorentz"), i18n("Hyperbolic secant (sech)"), i18n("Logistic (sech squared)"), i18n("Voigt profile"), i18n("Pseudo-Voigt (same width)")}; const char* nsl_fit_model_peak_equation[] = {"a/sqrt(2*pi)/s * exp(-((x-mu)/s)^2/2)", "a/pi * g/(g^2+(x-mu)^2)", "a/pi/s * sech((x-mu)/s)", - "a/4/s * sech((x-mu)/2/s)**2", "a*voigt(x - mu, s, g)", "a*pseudovoigt1(x - mu, eta, w)"}; + "a/4/s * sech((x-mu)/2/s)**2", "a*voigt(x - mu, s, g)", "a*pseudovoigt1(x - mu, et, w)"}; // eta is already used as function const char* nsl_fit_model_peak_pic_name[] = {"gaussian", "cauchy_lorentz", "sech", "logistic", "voigt", "pseudovoigt1"}; const char* nsl_fit_model_growth_name[] = {i18n("Inverse tangent"), i18n("Hyperbolic tangent"), i18n("Algebraic sigmoid"), i18n("Logistic function"), i18n("Error function (erf)"), i18n("Hill"), i18n("Gompertz"), i18n("Gudermann (gd)")}; const char* nsl_fit_model_growth_equation[] = {"a * atan((x-mu)/s)", "a * tanh((x-mu)/s)", "a * (x-mu)/s/sqrt(1+((x-mu)/s)^2)", "a/(1+exp(-k*(x-mu)))", "a/2 * erf((x-mu)/s/sqrt(2))", "a * x^n/(s^n + x^n)", "a*exp(-b*exp(-c*x))", "a * asin(tanh((x-mu)/s))"}; const char* nsl_fit_model_growth_pic_name[] = {"atan", "tanh", "alg_sigmoid", "logistic_function", "erf", "hill", "gompertz", "gd"}; const char* nsl_fit_weight_type_name[] = {"No", "Instrumental (1/col^2)", "Direct (col)", "Inverse (1/col)", "Statistical (1/data)", "Statistical (Fit)", "Relative (1/data^2)", "Relative (Fit)"}; /* see http://seal.web.cern.ch/seal/documents/minuit/mnusersguide.pdf and https://lmfit.github.io/lmfit-py/bounds.html */ double nsl_fit_map_bound(double x, double min, double max) { if (max <= min) { printf("given bounds must fulfill max > min (min = %g, max = %g)! Giving up.\n", min, max); return DBL_MAX; } /* not bounded */ if (min == -DBL_MAX && max == DBL_MAX) return x; /* open bounds */ if (min == -DBL_MAX) return max + 1. - sqrt(x*x + 1.); if (max == DBL_MAX) return min - 1. + sqrt(x*x + 1.); return min + sin(x + 1.) * (max - min)/2.; /* alternative transformation for closed bounds return min + (max - min)/(1. + exp(-x)); */ } /* see http://seal.web.cern.ch/seal/documents/minuit/mnusersguide.pdf and https://lmfit.github.io/lmfit-py/bounds.html */ double nsl_fit_map_unbound(double x, double min, double max) { if (max <= min) { printf("given bounds must fulfill max > min (min = %g, max = %g)! Giving up.\n", min, max); return DBL_MAX; } if (x < min || x > max) { printf("given value must be within bounds! Giving up.\n"); return -DBL_MAX; } /* not bounded */ if (min == -DBL_MAX && max == DBL_MAX) return x; /* open bounds */ if (min == -DBL_MAX) return sqrt(gsl_pow_2(max - x + 1.) - 1.); if (max == DBL_MAX) return sqrt(gsl_pow_2(x - min + 1.) - 1.); return asin(2. * (x - min)/(max - min) - 1.); /* alternative transformation for closed bounds return -log((max - x)/(x - min)); */ } /********************** parameter derivatives ******************/ /* basic */ double nsl_fit_model_polynomial_param_deriv(double x, int j, double weight) { return sqrt(weight)*pow(x, j); } double nsl_fit_model_power1_param_deriv(unsigned int param, double x, double a, double b, double weight) { if (param == 0) return sqrt(weight)*pow(x, b); if (param == 1) return sqrt(weight)*a*pow(x, b)*log(x); return 0; } double nsl_fit_model_power2_param_deriv(unsigned int param, double x, double b, double c, double weight) { if (param == 0) return sqrt(weight); if (param == 1) return sqrt(weight)*pow(x, c); if (param == 2) return sqrt(weight)*b*pow(x, c)*log(x); return 0; } double nsl_fit_model_exponentialn_param_deriv(unsigned int param, double x, double *p, double weight) { if (param % 2 == 0) return sqrt(weight)*exp(p[param+1]*x); else return sqrt(weight)*p[param-1]*x*exp(p[param]*x); } double nsl_fit_model_inverse_exponential_param_deriv(unsigned int param, double x, double a, double b, double weight) { if (param == 0) return sqrt(weight)*(1. - exp(b*x)); if (param == 1) return -sqrt(weight)*a*x*exp(b*x); if (param == 2) return sqrt(weight); return 0; } double nsl_fit_model_fourier_param_deriv(unsigned int param, unsigned int degree, double x, double w, double weight) { if (param == 0) return sqrt(weight)*cos(degree*w*x); if (param == 1) return sqrt(weight)*sin(degree*w*x); return 0; } /* peak */ double nsl_fit_model_gaussian_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { double s2 = s*s, norm = sqrt(weight)/sqrt(2.*M_PI)/s, efactor = exp(-(x-mu)*(x-mu)/(2.*s2)); if (param == 0) return norm * efactor; if (param == 1) return A * norm/(s*s2) * ((x-mu)*(x-mu) - s2) * efactor; if (param == 2) return A * norm/s2 * (x-mu) * efactor; return 0; } double nsl_fit_model_lorentz_param_deriv(unsigned int param, double x, double A, double s, double t, double weight) { double norm = sqrt(weight)/M_PI, denom = s*s+(x-t)*(x-t); if (param == 0) return norm * s/denom; if (param == 1) return A * norm * ((x-t)*(x-t) - s*s)/(denom*denom); if (param == 2) return A * norm * 2.*s*(x-t)/(denom*denom); return 0; } double nsl_fit_model_sech_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { double y = (x-mu)/s, norm = sqrt(weight)/M_PI/s; if (param == 0) return norm/cosh(y); if (param == 1) return A/s * norm * (y*tanh(y)-1.)/cosh(y); if (param == 2) return A/s * norm * tanh(y)/cosh(y); return 0; } double nsl_fit_model_logistic_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { double y = (x-mu)/2./s, norm = sqrt(weight)/4./s; if (param == 0) return norm/cosh(y)/cosh(y); if (param == 1) return A/s * norm * (2.*y*tanh(y)-1.)/cosh(y); if (param == 2) return A/s * norm * tanh(y)/cosh(y)/cosh(y); return 0; } double nsl_fit_model_voigt_param_deriv(unsigned int param, double x, double a, double mu, double s, double g, double weight) { #if !defined(_MSC_VER) if (s <= 0 || g < 0) return 0; double y = x - mu, norm = a * sqrt(weight/2./M_PI)/(s*s*s); double v = nsl_sf_voigt(y, s, g), im_w = nsl_sf_im_w_of_z(y); if (param == 0) return sqrt(weight) * v; if (param == 1) return a*sqrt(weight)*y/(s*s) * v - norm * g * im_w; if (param == 2) // return a*sqrt(weight)*g/M_PI/(s*s*s) + norm*sqrt(2.*M_PI)*v * (y*y+mu*mu-s*s) - norm/s*im_w*2.*g*y; return a/(s*s*s)*sqrt(weight)*(g/M_PI + v*(y*y -g*g -s*s) + im_w*2.*g*y/s); if (param == 3) return -a*sqrt(weight)/M_PI/(s*s) + norm*sqrt(2.*M_PI)*s*v*g + im_w; #endif return 0; } double nsl_fit_model_pseudovoigt1_param_deriv(unsigned int param, double x, double a, double eta, double w, double mu, double weight) { double y = x - mu, norm = sqrt(weight); - double sigma = w/sqrt(2.*log(2.)); + double sigma = w/sqrt(2.*M_LN2); if (param == 0) return norm * nsl_sf_pseudovoigt1(y, eta, w); if (param == 1) - return a * norm * (gsl_ran_cauchy_pdf(x, w) - gsl_ran_gaussian_pdf(x, sigma)); + return a * norm * (gsl_ran_cauchy_pdf(y, w) - gsl_ran_gaussian_pdf(y, sigma)); if (param == 2) - return a/w * norm * (eta*(1.-2.*w*w)*gsl_ran_cauchy_pdf(x, w) + (eta-1.)*gsl_ran_gaussian_pdf(x, sigma)*(1.-2.*M_LN2*y*y/w/w)); + return a/w * norm * (eta*(1.-2.*w*w)*gsl_ran_cauchy_pdf(y, w) + (eta-1.)*gsl_ran_gaussian_pdf(y, sigma)*(1.-2.*M_LN2*y*y/w/w)); if (param == 3) - return 2.*a*y/w/w * norm * (eta*M_PI*w*gsl_pow_2(gsl_ran_cauchy_pdf(x, w)) + (1.-eta)*M_LN2*gsl_ran_gaussian_pdf(x, sigma)); + return 2.*a*y/w/w * norm * (eta*M_PI*w*gsl_pow_2(gsl_ran_cauchy_pdf(y, w)) + (1.-eta)*M_LN2*gsl_ran_gaussian_pdf(y, sigma)); return 0; } /* growth */ double nsl_fit_model_atan_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight) { double norm = sqrt(weight), y = (x-mu)/s; if (param == 0) return norm * atan(y); if (param == 1) return -A/s * norm * 1./(1+y*y); if (param == 2) return -A/s * norm * y/(1.+y*y); return 0; } double nsl_fit_model_tanh_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight) { double norm = sqrt(weight), y = (x-mu)/s; if (param == 0) return norm * tanh(y); if (param == 1) return -A/s * norm * 1./cosh(y)/cosh(y); if (param == 2) return -A/s * norm * y/cosh(y)/cosh(y); return 0; } double nsl_fit_model_algebraic_sigmoid_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight) { double norm = sqrt(weight), y = (x-mu)/s, y2 = y*y; if (param == 0) return norm * y/sqrt(1.+y2); if (param == 1) return -A/s * norm * 1./pow(1.+y2, 1.5); if (param == 2) return -A/s * norm * y/pow(1.+y2, 1.5); return 0; } double nsl_fit_model_sigmoid_param_deriv(unsigned int param, double x, double A, double mu, double k, double weight) { double norm = sqrt(weight), y = k*(x-mu); if (param == 0) return norm/(1. + exp(-y)); if (param == 1) return -A*k * norm * exp(-y)/gsl_pow_2(1. + exp(-y)); if (param == 2) return A/k * norm * y*exp(-y)/gsl_pow_2(1. + exp(-y)); return 0; } double nsl_fit_model_erf_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight) { double norm = sqrt(weight), y = (x-mu)/(sqrt(2.)*s); if (param == 0) return norm/2. * gsl_sf_erf(y); if (param == 1) return -A/sqrt(2.*M_PI)/s * norm * exp(-y*y); if (param == 2) return -A/sqrt(M_PI)/s * norm * y*exp(-y*y); return 0; } double nsl_fit_model_hill_param_deriv(unsigned int param, double x, double A, double n, double s, double weight) { double norm = sqrt(weight), y = x/s; if (param == 0) return norm * pow(y, n)/(1.+pow(y, n)); if (param == 1) return A * norm * log(y)*pow(y, n)/gsl_pow_2(1.+pow(y, n)); if (param == 2) return -A*n/s * norm * pow(y, n)/gsl_pow_2(1.+pow(y, n)); return 0; } double nsl_fit_model_gompertz_param_deriv(unsigned int param, double x, double a, double b, double c, double weight) { if (param == 0) return sqrt(weight)*exp(-b*exp(-c*x)); if (param == 1) return -sqrt(weight)*a*exp(-c*x-b*exp(-c*x)); if (param == 2) return sqrt(weight)*a*b*x*exp(-c*x-b*exp(-c*x)); return 0; } double nsl_fit_model_gudermann_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight) { double norm = sqrt(weight), y = (x-mu)/s; if (param == 0) return -asin(tanh(y)); if (param == 1) return -A/s * norm * 1./cosh(y); if (param == 2) return -A/s * norm * y/cosh(y); return 0; } /* distributions */ double nsl_fit_model_gaussian_tail_param_deriv(unsigned int param, double x, double A, double s, double a, double mu, double weight) { if (x < a) return 0; double s2 = s*s, N = erfc(a/s/M_SQRT2)/2., norm = sqrt(weight)/sqrt(2.*M_PI)/s/N, efactor = exp(-(x-mu)*(x-mu)/(2.*s2)); if (param == 0) return norm * efactor; if (param == 1) return A * norm/(s*s2) * ((x-mu)*(x-mu) - s2) * efactor; if (param == 2) return A/norm/norm * efactor * exp(-a*a/(2.*s2)); if (param == 3) return A * norm/s2 * (x-mu) * efactor; return 0; } double nsl_fit_model_exponential_param_deriv(unsigned int param, double x, double A, double l, double mu, double weight) { if (x < mu) return 0; double y = l*(x-mu), efactor = exp(-y); if (param == 0) return sqrt(weight) * l * efactor; if (param == 1) return sqrt(weight) * A * (1. - y) * efactor; if (param == 2) return sqrt(weight) * A * gsl_pow_2(l) * efactor; return 0; } double nsl_fit_model_laplace_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { double norm = sqrt(weight)/(2.*s), y = fabs((x-mu)/s), efactor = exp(-y); if (param == 0) return norm * efactor; if (param == 1) return A/s*norm * (y-1.) * efactor; if (param == 2) return A/(s*s)*norm * (x-mu)/y * efactor; return 0; } double nsl_fit_model_exp_pow_param_deriv(unsigned int param, double x, double a, double s, double b, double mu, double weight) { double norm = sqrt(weight)/2./s/gsl_sf_gamma(1.+1./b), y = (x-mu)/s, efactor = exp(-pow(fabs(y), b)); if (param == 0) return norm * efactor; if (param == 1) return norm * a/s * efactor * (b * y * pow(fabs(1./y), 1.-b) * GSL_SIGN(y) - 1.); if (param == 2) return norm * a/b * gsl_sf_gamma(1.+1./b)/gsl_sf_gamma(1./b) * efactor * (gsl_sf_psi(1.+1./b) - gsl_pow_2(b) * pow(fabs(y), b) * log(fabs(y))); if (param == 3) return norm * a*b/s * efactor * pow(fabs(y), b-1.) * GSL_SIGN(y); return 0; } double nsl_fit_model_maxwell_param_deriv(unsigned int param, double x, double a, double s, double weight) { double s2 = s*s, s3 = s*s2, norm = sqrt(weight)*sqrt(2./M_PI)/s3, x2 = x*x, efactor = exp(-x2/2./s2); if (param == 0) return norm * x2 * efactor; if (param == 1) return a * norm * x2*(x2-3.*s2)/s3 * efactor; return 0; } double nsl_fit_model_poisson_param_deriv(unsigned int param, double x, double A, double l, double weight) { double norm = sqrt(weight)*pow(l, x)/gsl_sf_gamma(x+1.); if (param == 0) return norm * exp(-l); if (param == 1) return A/l * norm *(x-l)*exp(-l); return 0; } double nsl_fit_model_lognormal_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { double norm = sqrt(weight)/sqrt(2.*M_PI)/(x*s), y = log(x)-mu, efactor = exp(-(y/s)*(y/s)/2.); if (param == 0) return norm * efactor; if (param == 1) return A * norm * (y*y - s*s) * efactor; if (param == 2) return A * norm * y/(s*s) * efactor; return 0; } double nsl_fit_model_gamma_param_deriv(unsigned int param, double x, double A, double k, double t, double weight) { double factor = sqrt(weight)*pow(x, k-1.)/pow(t, k)/gsl_sf_gamma(k), efactor = exp(-x/t); if (param == 0) return factor * efactor; if (param == 1) return A * factor * (log(x/t) - gsl_sf_psi(k)) * efactor; if (param == 2) return A * factor/t * (x/t-k) * efactor; return 0; } double nsl_fit_model_flat_param_deriv(unsigned int param, double x, double A, double b, double a, double weight) { if (x < a || x > b) return 0; if (param == 0) return sqrt(weight)/(b-a); if (param == 1) return - sqrt(weight) * A/gsl_pow_2(a-b); if (param == 2) return sqrt(weight) * A/gsl_pow_2(a-b); return 0; } double nsl_fit_model_rayleigh_param_deriv(unsigned int param, double x, double A, double s, double weight) { double y=x/s, norm = sqrt(weight)*y/s, efactor = exp(-y*y/2.); if (param == 0) return norm * efactor; if (param == 1) return A*y/(s*s) * (y*y-2.)*efactor; return 0; } double nsl_fit_model_rayleigh_tail_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { double norm = sqrt(weight)*x/(s*s), y = (mu*mu - x*x)/2./(s*s); if (param == 0) return norm * exp(y); if (param == 1) return -2. * A * norm/s * (1. + y) * exp(y); if (param == 2) return A * mu * norm/(s*s) * exp(y); return 0; } double nsl_fit_model_levy_param_deriv(unsigned int param, double x, double A, double g, double mu, double weight) { double y=x-mu, norm = sqrt(weight)*sqrt(g/(2.*M_PI))/pow(y, 1.5), efactor = exp(-g/2./y); if (param == 0) return norm * efactor; if (param == 1) return A/2.*norm/g/y * (y - g) * efactor; if (param == 2) return A/2.*norm/y/y * (3.*y - g) * efactor; return 0; } double nsl_fit_model_landau_param_deriv(unsigned int param, double x, double weight) { if (param == 0) return sqrt(weight) * gsl_ran_landau_pdf(x); return 0; } double nsl_fit_model_chi_square_param_deriv(unsigned int param, double x, double A, double n, double weight) { double y=n/2., norm = sqrt(weight)*pow(x, y-1.)/pow(2., y)/gsl_sf_gamma(y), efactor = exp(-x/2.); if (param == 0) return norm * efactor; if (param == 1) return A/2. * norm * (log(x/2.) - gsl_sf_psi(y)) * efactor; return 0; } double nsl_fit_model_students_t_param_deriv(unsigned int param, double x, double A, double n, double weight) { if (param == 0) return sqrt(weight) * gsl_ran_tdist_pdf(x, n); if (param == 1) return sqrt(weight) * A * gsl_sf_gamma((n+1.)/2.)/2./pow(n, 1.5)/sqrt(M_PI)/gsl_sf_gamma(n/2.) * pow(1.+x*x/n, - (n+3.)/2.) * (x*x - 1. - (n+x*x)*log(1.+x*x/n) + (n+x*x)*(gsl_sf_psi((n+1.)/2.) - gsl_sf_psi(n/2.)) ) ; return 0; } double nsl_fit_model_fdist_param_deriv(unsigned int param, double x, double A, double n1, double n2, double weight) { double norm = sqrt(weight) * gsl_sf_gamma((n1+n2)/2.)/gsl_sf_gamma(n1/2.)/gsl_sf_gamma(n2/2.) * pow(n1, n1/2.) * pow(n2, n2/2.) * pow(x, n1/2.-1.); double y = n2+n1*x; if (param == 0) return sqrt(weight) * gsl_ran_fdist_pdf(x, n1, n2); if (param == 1) return A/2. * norm * pow(y, -(n1+n2+2.)/2.) * (n2*(1.-x) + y*(log(n1) + log(x) - log(y) + gsl_sf_psi((n1+n2)/2.) - gsl_sf_psi(n1/2.))); if (param == 2) return A/2. * norm * pow(y, -(n1+n2+2.)/2.) * (n1*(x-1.) + y*(log(n2) - log(y) + gsl_sf_psi((n1+n2)/2.) - gsl_sf_psi(n2/2.))); return 0; } double nsl_fit_model_beta_param_deriv(unsigned int param, double x, double A, double a, double b, double weight) { double norm = sqrt(weight) * A * gsl_sf_gamma(a+b)/gsl_sf_gamma(a)/gsl_sf_gamma(b) * pow(x, a-1.) * pow(1.-x, b-1.); if (param == 0) return sqrt(weight) * gsl_ran_beta_pdf(x, a, b); if (param == 1) return norm * (log(x) - gsl_sf_psi(a) + gsl_sf_psi(a+b)); if (param == 2) return norm * (log(1.-x) - gsl_sf_psi(b) + gsl_sf_psi(a+b)); return 0; } double nsl_fit_model_pareto_param_deriv(unsigned int param, double x, double A, double a, double b, double weight) { if (x < b) return 0; double norm = sqrt(weight) * A; if (param == 0) return sqrt(weight) * gsl_ran_pareto_pdf(x, a, b); if (param == 1) return norm * pow(b/x, a) * (1. + a * log(b/x))/x; if (param == 2) return norm * a*a * pow(b/x, a-1.)/x/x; return 0; } double nsl_fit_model_weibull_param_deriv(unsigned int param, double x, double A, double k, double l, double mu, double weight) { double y = (x-mu)/l, z = pow(y, k), efactor = exp(-z); if (param == 0) return sqrt(weight) * k/l * z/y * efactor; if (param == 1) return sqrt(weight) * A/l * z/y*(k*log(y)*(1.-z) + 1.) * efactor; if (param == 2) return sqrt(weight) * A*k*k/l/l * z/y*(z-1.) * efactor; if (param == 3) return sqrt(weight) * A*k/l/l * z/y/y*(k*z + 1. - k) * efactor; return 0; } double nsl_fit_model_frechet_param_deriv(unsigned int param, double x, double A, double g, double s, double mu, double weight) { double y = (x-mu)/s, efactor = exp(-pow(y, -g)); if (param == 0) return g * sqrt(weight)/s * pow(y, -g-1.) * efactor; if (param == 1) return sqrt(weight) * A/s * pow(y, -2.*g-1.) * (g*log(y)*(1.-pow(y, g))+pow(y, g)) * efactor; if (param == 2) return A * sqrt(weight) * gsl_pow_2(g/s)*pow(y, -2.*g-1.) * (pow(y, g)-1.) * efactor; if (param == 3) return A * sqrt(weight) * g/(s*s)*pow(y, -g-2.) * (g+1.-g*pow(y, -g)) * efactor; return 0; } double nsl_fit_model_gumbel1_param_deriv(unsigned int param, double x, double A, double s, double mu, double b, double weight) { double norm = sqrt(weight)/s, y = (x-mu)/s, efactor = exp(-y - b*exp(-y)); if (param == 0) return norm * efactor; if (param == 1) return A/s * norm * (y - 1. - b*exp(-y)) * efactor; if (param == 2) return A/s * norm * (1. - b*exp(-y)) * efactor; if (param == 3) return -A * norm * exp(-y) * efactor; return 0; } double nsl_fit_model_gumbel2_param_deriv(unsigned int param, double x, double A, double a, double b, double mu, double weight) { double y = x - mu, norm = A * sqrt(weight) * exp(-b * pow(y, -a)); if (param == 0) return sqrt(weight) * gsl_ran_gumbel2_pdf(y, a, b); if (param == 1) return norm * b * pow(y, -1. -2.*a) * (pow(y, a) -a*(pow(y, a)-b)*log(y)); if (param == 2) return norm * a * pow(y, -1. -2.*a) * (pow(y, a) - b); if (param == 3) return norm * a * b * pow(y, -2.*(a + 1.)) * ((1. + a)*pow(y, a) - a*b); return 0; } double nsl_fit_model_binomial_param_deriv(unsigned int param, double k, double A, double p, double n, double weight) { if (k < 0 || k > n || n < 0 || p < 0 || p > 1.) return 0; k = round(k); n = round(n); double norm = sqrt(weight) * gsl_sf_fact((unsigned int)n)/gsl_sf_fact((unsigned int)(n-k))/gsl_sf_fact((unsigned int)(k)); if (param == 0) return sqrt(weight) * gsl_ran_binomial_pdf((unsigned int)k, p, (unsigned int)n); if (param == 1) return A * norm * pow(p, k-1.) * pow(1.-p, n-k-1.) * (k-n*p); if (param == 2) return A * norm * pow(p, k) * pow(1.-p, n-k) * (log(1.-p) + gsl_sf_psi(n+1.) - gsl_sf_psi(n-k+1.)); return 0; } double nsl_fit_model_negative_binomial_param_deriv(unsigned int param, double k, double A, double p, double n, double weight) { if (k < 0 || k > n || n < 0 || p < 0 || p > 1.) return 0; double norm = A * sqrt(weight) * gsl_sf_gamma(n+k)/gsl_sf_gamma(k+1.)/gsl_sf_gamma(n); if (param == 0) return sqrt(weight) * gsl_ran_negative_binomial_pdf((unsigned int)k, p, n); if (param == 1) return - norm * pow(p, n-1.) * pow(1.-p, k-1.) * (n*(p-1.) + k*p); if (param == 2) return norm * pow(p, n) * pow(1.-p, k) * (log(p) - gsl_sf_psi(n) + gsl_sf_psi(n+k)); return 0; } double nsl_fit_model_pascal_param_deriv(unsigned int param, double k, double A, double p, double n, double weight) { return nsl_fit_model_negative_binomial_param_deriv(param, k, A, p, round(n), weight); } double nsl_fit_model_geometric_param_deriv(unsigned int param, double k, double A, double p, double weight) { if (param == 0) return sqrt(weight) * gsl_ran_geometric_pdf((unsigned int)k, p); if (param == 1) return A * sqrt(weight) * pow(1.-p, k-2.) * (1.-k*p); return 0; } double nsl_fit_model_hypergeometric_param_deriv(unsigned int param, double k, double A, double n1, double n2, double t, double weight) { if (t > n1 + n2) return 0; double norm = sqrt(weight) * gsl_ran_hypergeometric_pdf((unsigned int)k, (unsigned int)n1, (unsigned int)n2, (unsigned int)t); if (param == 0) return norm; if (param == 1) return A * norm * (gsl_sf_psi(n1+1.) - gsl_sf_psi(n1-k+1.) - gsl_sf_psi(n1+n2+1.) + gsl_sf_psi(n1+n2-t+1.)); if (param == 2) return A * norm * (gsl_sf_psi(n2+1.) - gsl_sf_psi(n2+k-t+1.) - gsl_sf_psi(n1+n2+1.) + gsl_sf_psi(n1+n2-t+1.)); if (param == 3) return A * norm * (gsl_sf_psi(n2+k-t+1.) - gsl_sf_psi(n1+n2-t+1.) - gsl_sf_psi(t-k+1.) + gsl_sf_psi(t+1.)); return 0; } double nsl_fit_model_logarithmic_param_deriv(unsigned int param, double k, double A, double p, double weight) { if (param == 0) return sqrt(weight) * gsl_ran_logarithmic_pdf((unsigned int)k, p); if (param == 1) return A * sqrt(weight) * pow(1.-p, k-2.) * (1.-k*p); return 0; } double nsl_fit_model_sech_dist_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { double norm = sqrt(weight)/2./s, y = M_PI/2.*(x-mu)/s; if (param == 0) return norm * 1./cosh(y); if (param == 1) return -A/s * norm * (y*tanh(y)+1.)/cosh(y); if (param == 2) return A*M_PI/2./s * norm * tanh(y)/cosh(y); return 0; } diff --git a/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp b/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp index 2b53648e7..e1ed5c067 100644 --- a/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp +++ b/src/backend/worksheet/plots/cartesian/XYFitCurve.cpp @@ -1,2396 +1,2396 @@  /*************************************************************************** File : XYFitCurve.cpp Project : LabPlot Description : A xy-curve defined by a fit model -------------------------------------------------------------------- Copyright : (C) 2014-2017 Alexander Semke (alexander.semke@web.de) Copyright : (C) 2016-2018 Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ /*! \class XYFitCurve \brief A xy-curve defined by a fit model \ingroup worksheet */ #include "XYFitCurve.h" #include "XYFitCurvePrivate.h" #include "backend/core/AbstractColumn.h" #include "backend/core/column/Column.h" #include "backend/lib/commandtemplates.h" #include "backend/lib/macros.h" #include "backend/gsl/errors.h" #include "backend/gsl/ExpressionParser.h" extern "C" { #include #include #include #include #include #include #include "backend/gsl/parser.h" #include "backend/nsl/nsl_sf_stats.h" #include "backend/nsl/nsl_stats.h" } #include #include #include XYFitCurve::XYFitCurve(const QString& name) : XYAnalysisCurve(name, new XYFitCurvePrivate(this)) { } XYFitCurve::XYFitCurve(const QString& name, XYFitCurvePrivate* dd) : XYAnalysisCurve(name, dd) { } XYFitCurve::~XYFitCurve() { //no need to delete the d-pointer here - it inherits from QGraphicsItem //and is deleted during the cleanup in QGraphicsScene } void XYFitCurve::recalculate() { Q_D(XYFitCurve); d->recalculate(); } void XYFitCurve::evaluate(bool preview) { Q_D(XYFitCurve); d->evaluate(preview); } void XYFitCurve::initStartValues(const XYCurve* curve) { Q_D(XYFitCurve); XYFitCurve::FitData& fitData = d->fitData; initStartValues(fitData, curve); } void XYFitCurve::initStartValues(XYFitCurve::FitData& fitData, const XYCurve* curve) { DEBUG("XYFitCurve::initStartValues()"); if (!curve) { DEBUG(" no curve given"); return; } const Column* tmpXDataColumn = dynamic_cast(curve->xColumn()); const Column* tmpYDataColumn = dynamic_cast(curve->yColumn()); if (!tmpXDataColumn || !tmpYDataColumn) { DEBUG(" data columns not available"); return; } DEBUG(" x data rows = " << tmpXDataColumn->rowCount()); nsl_fit_model_category modelCategory = fitData.modelCategory; int modelType = fitData.modelType; int degree = fitData.degree; DEBUG(" fit model type = " << modelType << ", degree = " << degree); QVector& paramStartValues = fitData.paramStartValues; //QVector* xVector = static_cast* >(tmpXDataColumn->data()); //double xmean = gsl_stats_mean(xVector->constData(), 1, tmpXDataColumn->rowCount()); double xmin = tmpXDataColumn->minimum(); double xmax = tmpXDataColumn->maximum(); //double ymin = tmpYDataColumn->minimum(); //double ymax = tmpYDataColumn->maximum(); double xrange = xmax-xmin; //double yrange = ymax-ymin; DEBUG(" x min/max = " << xmin << ' ' << xmax); //DEBUG(" y min/max = " << ymin << ' ' << ymax); switch (modelCategory) { case nsl_fit_model_basic: switch (modelType) { case nsl_fit_model_polynomial: // not needed (works anyway) break; //TODO: handle basic models case nsl_fit_model_power: case nsl_fit_model_exponential: case nsl_fit_model_inverse_exponential: case nsl_fit_model_fourier: break; } break; case nsl_fit_model_peak: // use equidistant mu's and (xmax-xmin)/(10*degree) as sigma(, gamma) switch (modelType) { case nsl_fit_model_gaussian: case nsl_fit_model_lorentz: case nsl_fit_model_sech: case nsl_fit_model_logistic: for (int d = 0; d < degree; d++) { paramStartValues[3*d+2] = xmin + (d+1.)*xrange/(degree+1.); // mu paramStartValues[3*d+1] = xrange/(10.*degree); // sigma } break; case nsl_fit_model_voigt: for (int d = 0; d < degree; d++) { paramStartValues[4*d+1] = xmin + (d+1.)*xrange/(degree+1.); // mu paramStartValues[4*d+2] = xrange/(10.*degree); // sigma paramStartValues[4*d+3] = xrange/(10.*degree); // gamma } break; case nsl_fit_model_pseudovoigt1: for (int d = 0; d < degree; d++) { paramStartValues[4*d+1] = 0.5; // eta paramStartValues[4*d+2] = xrange/(10.*degree); // sigma paramStartValues[4*d+3] = xmin + (d+1.)*xrange/(degree+1.); // mu } break; } break; case nsl_fit_model_growth: switch (modelType) { case nsl_fit_model_atan: case nsl_fit_model_tanh: case nsl_fit_model_algebraic_sigmoid: case nsl_fit_model_erf: case nsl_fit_model_gudermann: case nsl_fit_model_sigmoid: // use (xmax+xmin)/2 as mu and (xmax-xmin)/10 as sigma paramStartValues[1] = (xmax+xmin)/2.; paramStartValues[2] = xrange/10.; break; case nsl_fit_model_hill: paramStartValues[2] = xrange/10.; break; case nsl_fit_model_gompertz: //TODO break; } break; case nsl_fit_model_distribution: switch (modelType) { case nsl_sf_stats_gaussian: case nsl_sf_stats_laplace: case nsl_sf_stats_rayleigh_tail: case nsl_sf_stats_lognormal: case nsl_sf_stats_logistic: case nsl_sf_stats_sech: case nsl_sf_stats_cauchy_lorentz: case nsl_sf_stats_levy: // use (xmax+xmin)/2 as mu and (xmax-xmin)/10 as sigma paramStartValues[2] = (xmin+xmax)/2.; paramStartValues[1] = xrange/10.; break; //TODO: other types default: break; } break; case nsl_fit_model_custom: // not possible break; } } /*! * sets the parameter names for given model category, model type and degree in \c fitData for given action */ void XYFitCurve::initFitData(PlotDataDialog::AnalysisAction action) { if (!action) return; Q_D(XYFitCurve); XYFitCurve::FitData& fitData = d->fitData; if (action == PlotDataDialog::FitLinear) { //Linear fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = (int)nsl_fit_model_polynomial; fitData.degree = 1; } else if (action == PlotDataDialog::FitPower) { //Power fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = (int)nsl_fit_model_power; fitData.degree = 1; } else if (action == PlotDataDialog::FitExp1) { //Exponential (degree 1) fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = (int)nsl_fit_model_exponential; fitData.degree = 1; } else if (action == PlotDataDialog::FitExp2) { //Exponential (degree 2) fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = (int)nsl_fit_model_exponential; fitData.degree = 2; } else if (action == PlotDataDialog::FitInvExp) { //Inverse exponential fitData.modelCategory = nsl_fit_model_basic; fitData.modelType = (int)nsl_fit_model_inverse_exponential; } else if (action == PlotDataDialog::FitGauss) { //Gauss fitData.modelCategory = nsl_fit_model_peak; fitData.modelType = (int)nsl_fit_model_gaussian; fitData.degree = 1; } else if (action == PlotDataDialog::FitCauchyLorentz) { //Cauchy-Lorentz fitData.modelCategory = nsl_fit_model_peak; fitData.modelType = (int)nsl_fit_model_lorentz; fitData.degree = 1; } else if (action == PlotDataDialog::FitTan) { //Arc tangent fitData.modelCategory = nsl_fit_model_growth; fitData.modelType = (int)nsl_fit_model_atan; } else if (action == PlotDataDialog::FitTanh) { //Hyperbolic tangent fitData.modelCategory = nsl_fit_model_growth; fitData.modelType = (int)nsl_fit_model_tanh; } else if (action == PlotDataDialog::FitErrFunc) { //Error function fitData.modelCategory = nsl_fit_model_growth; fitData.modelType = (int)nsl_fit_model_erf; } else { //Custom fitData.modelCategory = nsl_fit_model_custom; fitData.modelType = 0; } XYFitCurve::initFitData(fitData); } /*! * sets the model expression and the parameter names for given model category, model type and degree in \c fitData */ void XYFitCurve::initFitData(XYFitCurve::FitData& fitData) { nsl_fit_model_category modelCategory = fitData.modelCategory; int modelType = fitData.modelType; QString& model = fitData.model; QStringList& paramNames = fitData.paramNames; QStringList& paramNamesUtf8 = fitData.paramNamesUtf8; int degree = fitData.degree; QVector& paramStartValues = fitData.paramStartValues; QVector& paramLowerLimits = fitData.paramLowerLimits; QVector& paramUpperLimits = fitData.paramUpperLimits; QVector& paramFixed = fitData.paramFixed; if (modelCategory != nsl_fit_model_custom) { DEBUG("XYFitCurve::initFitData() for model category = " << nsl_fit_model_category_name[modelCategory] << ", model type = " << modelType << ", degree = " << degree); paramNames.clear(); } else { DEBUG("XYFitCurve::initFitData() for model category = nsl_fit_model_custom, model type = " << modelType << ", degree = " << degree); } paramNamesUtf8.clear(); // 10 indices used in multi degree models QStringList indices = { UTF8_QSTRING("₁"), UTF8_QSTRING("₂"), UTF8_QSTRING("₃"), UTF8_QSTRING("₄"), UTF8_QSTRING("₅"), UTF8_QSTRING("₆"), UTF8_QSTRING("₇"), UTF8_QSTRING("₈"), UTF8_QSTRING("₉"), UTF8_QSTRING("₁₀")}; switch (modelCategory) { case nsl_fit_model_basic: model = nsl_fit_model_basic_equation[fitData.modelType]; switch (modelType) { case nsl_fit_model_polynomial: paramNames << "c0" << "c1"; paramNamesUtf8 << UTF8_QSTRING("c₀") << UTF8_QSTRING("c₁"); if (degree == 2) { model += " + c2*x^2"; paramNames << "c2"; paramNamesUtf8 << UTF8_QSTRING("c₂"); } else if (degree > 2) { for (int i = 2; i <= degree; ++i) { QString numStr = QString::number(i); model += "+c" + numStr + "*x^" + numStr; paramNames << 'c' + numStr; paramNamesUtf8 << 'c' + indices[i-1]; } } break; case nsl_fit_model_power: if (degree == 1) { paramNames << "a" << "b"; } else { paramNames << "a" << "b" << "c"; model = "a + b*x^c"; } break; case nsl_fit_model_exponential: if (degree == 1) { paramNames << "a" << "b"; } else { for (int i = 1; i <= degree; i++) { QString numStr = QString::number(i); if (i == 1) model = "a1*exp(b1*x)"; else model += " + a" + numStr + "*exp(b" + numStr + "*x)"; paramNames << 'a' + numStr << 'b' + numStr; paramNamesUtf8 << 'a' + indices[i-1] << 'b' + indices[i-1]; } } break; case nsl_fit_model_inverse_exponential: paramNames << "a" << "b" << "c"; break; case nsl_fit_model_fourier: paramNames << "w" << "a0" << "a1" << "b1"; paramNamesUtf8 << UTF8_QSTRING("ω") << UTF8_QSTRING("a₀") << UTF8_QSTRING("a₁") << UTF8_QSTRING("b₁"); if (degree > 1) { for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); model += "+ (a" + numStr + "*cos(" + numStr + "*w*x) + b" + numStr + "*sin(" + numStr + "*w*x))"; paramNames << 'a' + numStr << 'b' + numStr; paramNamesUtf8 << 'a' + indices[i-1] << 'b' + indices[i-1]; } } break; } break; case nsl_fit_model_peak: model = nsl_fit_model_peak_equation[fitData.modelType]; switch (modelType) { case nsl_fit_model_gaussian: switch (degree) { case 1: paramNames << "a" << "s" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ"); break; default: model = "1./sqrt(2*pi) * ("; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += 'a' + numStr + "/s" + numStr + "* exp(-((x-mu" + numStr + ")/s" + numStr + ")^2/2)"; paramNames << 'a' + numStr << 's' + numStr << "mu" + numStr; paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1]; } model += ')'; } break; case nsl_fit_model_lorentz: switch (degree) { case 1: paramNames << "a" << "g" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("γ") << UTF8_QSTRING("μ"); break; default: model = "1./pi * ("; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += 'a' + numStr + " * g" + numStr + "/(g" + numStr + "^2+(x-mu" + numStr + ")^2)"; paramNames << 'a' + numStr << 'g' + numStr << "mu" + numStr; paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("γ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1]; } model += ')'; } break; case nsl_fit_model_sech: switch (degree) { case 1: paramNames << "a" << "s" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ"); break; default: model = "1/pi * ("; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += 'a' + numStr + "/s" + numStr + "* sech((x-mu" + numStr + ")/s" + numStr + ')'; paramNames << 'a' + numStr << "s" + numStr << "mu" + numStr; paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1]; } model += ')'; } break; case nsl_fit_model_logistic: switch (degree) { case 1: paramNames << "a" << "s" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ"); break; default: model = "1/4 * ("; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += 'a' + numStr + "/s" + numStr + "* sech((x-mu" + numStr + ")/2/s" + numStr + ")**2"; paramNames << 'a' + numStr << 's' + numStr << "mu" + numStr; paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1]; } model += ')'; } break; case nsl_fit_model_voigt: switch(degree) { case 1: paramNames << "a" << "mu" << "s" << "g"; paramNamesUtf8 << "A" << UTF8_QSTRING("μ") << UTF8_QSTRING("σ") << UTF8_QSTRING("γ"); break; default: model = ""; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += 'a' + numStr + "*voigt(x-mu" + numStr + ",s" + numStr + ",g" + numStr + ')'; paramNames << "a" + numStr << "mu" + numStr << "s" + numStr << "g" + numStr; paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1] << UTF8_QSTRING("σ") + indices[i-1] << UTF8_QSTRING("γ") + indices[i-1]; } } break; case nsl_fit_model_pseudovoigt1: switch(degree) { case 1: - paramNames << "a" << "eta" << "w" << "mu"; + paramNames << "a" << "et" << "w" << "mu"; // eta function exists! paramNamesUtf8 << "A" << UTF8_QSTRING("η") << "w" << UTF8_QSTRING("μ"); break; default: model=""; for (int i = 1; i <= degree; ++i) { QString numStr = QString::number(i); if (i > 1) model += " + "; model += 'a' + numStr + "*pseudovoigt1(x-mu" + numStr + ",eta" + numStr + ",w" + numStr + ')'; paramNames << "a" + numStr << "eta" + numStr << "w" + numStr << "mu" + numStr; paramNamesUtf8 << 'A' + indices[i-1] << UTF8_QSTRING("η") + indices[i-1] << 'w' + indices[i-1] << UTF8_QSTRING("μ") + indices[i-1]; } } break; } break; case nsl_fit_model_growth: model = nsl_fit_model_growth_equation[fitData.modelType]; switch (modelType) { case nsl_fit_model_atan: case nsl_fit_model_tanh: case nsl_fit_model_algebraic_sigmoid: case nsl_fit_model_erf: case nsl_fit_model_gudermann: paramNames << "a" << "mu" << "s"; paramNamesUtf8 << "A" << UTF8_QSTRING("μ") << UTF8_QSTRING("σ"); break; case nsl_fit_model_sigmoid: paramNames << "a" << "mu" << "k"; paramNamesUtf8 << "A" << UTF8_QSTRING("μ") << "k"; break; case nsl_fit_model_hill: paramNames << "a" << "n" << "a"; paramNamesUtf8 << "A" << "n" << UTF8_QSTRING("σ"); break; case nsl_fit_model_gompertz: paramNames << "a" << "b" << "c"; break; } break; case nsl_fit_model_distribution: model = nsl_sf_stats_distribution_equation[fitData.modelType]; switch (modelType) { case nsl_sf_stats_gaussian: case nsl_sf_stats_laplace: case nsl_sf_stats_rayleigh_tail: case nsl_sf_stats_lognormal: case nsl_sf_stats_logistic: case nsl_sf_stats_sech: paramNames << "a" << "s" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ"); break; case nsl_sf_stats_gaussian_tail: paramNames << "A" << "s" << "a" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << "a" << UTF8_QSTRING("μ"); break; case nsl_sf_stats_exponential: paramNames << "a" << "l" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("λ") << UTF8_QSTRING("μ"); break; case nsl_sf_stats_exponential_power: paramNames << "a" << "s" << "b" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << "b" << UTF8_QSTRING("μ"); break; case nsl_sf_stats_cauchy_lorentz: case nsl_sf_stats_levy: paramNames << "a" << "g" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("γ") << UTF8_QSTRING("μ"); break; case nsl_sf_stats_rayleigh: paramNames << "a" << "s"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ"); break; case nsl_sf_stats_landau: paramNames << "a"; paramNamesUtf8 << "A"; break; case nsl_sf_stats_levy_alpha_stable: // unused distributions case nsl_sf_stats_levy_skew_alpha_stable: case nsl_sf_stats_bernoulli: break; case nsl_sf_stats_gamma: paramNames << "a" << "k" << "t"; paramNamesUtf8 << "A"<< "k" << UTF8_QSTRING("θ"); break; case nsl_sf_stats_flat: paramNames << "A" << "b" << "a"; break; case nsl_sf_stats_chi_squared: paramNames << "a" << "n"; paramNamesUtf8 << "A" << "n"; break; case nsl_sf_stats_fdist: paramNames << "a" << "n1" << "n2"; paramNamesUtf8 << "A" << UTF8_QSTRING("ν₁") << UTF8_QSTRING("ν₂"); break; case nsl_sf_stats_tdist: paramNames << "a" << "n"; paramNamesUtf8 << "A" << UTF8_QSTRING("ν"); break; case nsl_sf_stats_beta: case nsl_sf_stats_pareto: paramNames << "A" << "a" << "b"; break; case nsl_sf_stats_weibull: paramNames << "a" << "k" << "l" << "mu"; paramNamesUtf8 << "A" << "k" << UTF8_QSTRING("λ") << UTF8_QSTRING("μ"); break; case nsl_sf_stats_gumbel1: paramNames << "a" << "s" << "mu" << "b"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ") << UTF8_QSTRING("μ") << UTF8_QSTRING("β"); break; case nsl_sf_stats_gumbel2: paramNames << "A" << "a" << "b" << "mu"; paramNamesUtf8 << "A" << "a" << "b" << UTF8_QSTRING("μ"); break; case nsl_sf_stats_poisson: paramNames << "a" << "l"; paramNamesUtf8 << "A" << UTF8_QSTRING("λ"); break; case nsl_sf_stats_binomial: case nsl_sf_stats_negative_binomial: case nsl_sf_stats_pascal: paramNames << "a" << "p" << "n"; paramNamesUtf8 << "A" << "p" << "n"; break; case nsl_sf_stats_geometric: case nsl_sf_stats_logarithmic: paramNames << "a" << "p"; paramNamesUtf8 << "A" << "p"; break; case nsl_sf_stats_hypergeometric: paramNames << "a" << "n1" << "n2" << "t"; paramNamesUtf8 << "A" << UTF8_QSTRING("n₁") << UTF8_QSTRING("n₂") << "t"; break; case nsl_sf_stats_maxwell_boltzmann: paramNames << "a" << "s"; paramNamesUtf8 << "A" << UTF8_QSTRING("σ"); break; case nsl_sf_stats_frechet: paramNames << "a" << "g" << "s" << "mu"; paramNamesUtf8 << "A" << UTF8_QSTRING("γ") << UTF8_QSTRING("σ") << UTF8_QSTRING("μ"); break; } break; case nsl_fit_model_custom: break; } DEBUG("model: " << model.toStdString()); if (paramNamesUtf8.isEmpty()) paramNamesUtf8 << paramNames; //resize the vector for the start values and set the elements to 1.0 //in case a custom model is used, do nothing, we take over the previous values if (modelCategory != nsl_fit_model_custom) { const int np = paramNames.size(); paramStartValues.resize(np); paramFixed.resize(np); paramLowerLimits.resize(np); paramUpperLimits.resize(np); for (int i = 0; i < np; ++i) { paramStartValues[i] = 1.0; paramFixed[i] = false; paramLowerLimits[i] = -std::numeric_limits::max(); paramUpperLimits[i] = std::numeric_limits::max(); } // set some model-dependent start values // TODO: see initStartValues() if (modelCategory == nsl_fit_model_distribution) { if (modelType == (int)nsl_sf_stats_flat) paramStartValues[2] = -1.0; else if (modelType == (int)nsl_sf_stats_levy) paramStartValues[2] = 0.0; else if (modelType == (int)nsl_sf_stats_exponential_power || modelType == (int)nsl_sf_stats_weibull || modelType == (int)nsl_sf_stats_gumbel2 || modelType == (int)nsl_sf_stats_frechet) paramStartValues[3] = 0.0; else if (modelType == (int)nsl_sf_stats_binomial || modelType == (int)nsl_sf_stats_negative_binomial || modelType == (int)nsl_sf_stats_pascal || modelType == (int)nsl_sf_stats_geometric || modelType == (int)nsl_sf_stats_logarithmic) paramStartValues[1] = 0.5; } } } /*! Returns an icon to be used in the project explorer. */ QIcon XYFitCurve::icon() const { return QIcon::fromTheme("labplot-xy-fit-curve"); } //############################################################################## //########################## getter methods ################################## //############################################################################## BASIC_SHARED_D_READER_IMPL(XYFitCurve, const AbstractColumn*, xErrorColumn, xErrorColumn) BASIC_SHARED_D_READER_IMPL(XYFitCurve, const AbstractColumn*, yErrorColumn, yErrorColumn) const QString& XYFitCurve::xErrorColumnPath() const { Q_D(const XYFitCurve); return d->xErrorColumnPath; } const QString& XYFitCurve::yErrorColumnPath() const { Q_D(const XYFitCurve); return d->yErrorColumnPath; } BASIC_SHARED_D_READER_IMPL(XYFitCurve, XYFitCurve::FitData, fitData, fitData) const XYFitCurve::FitResult& XYFitCurve::fitResult() const { Q_D(const XYFitCurve); return d->fitResult; } //############################################################################## //################# setter methods and undo commands ########################## //############################################################################## STD_SETTER_CMD_IMPL_S(XYFitCurve, SetXErrorColumn, const AbstractColumn*, xErrorColumn) void XYFitCurve::setXErrorColumn(const AbstractColumn* column) { Q_D(XYFitCurve); if (column != d->xErrorColumn) { exec(new XYFitCurveSetXErrorColumnCmd(d, column, i18n("%1: assign x-error"))); handleSourceDataChanged(); if (column) { connect(column, SIGNAL(dataChanged(const AbstractColumn*)), this, SLOT(handleSourceDataChanged())); //TODO disconnect on undo } } } STD_SETTER_CMD_IMPL_S(XYFitCurve, SetYErrorColumn, const AbstractColumn*, yErrorColumn) void XYFitCurve::setYErrorColumn(const AbstractColumn* column) { Q_D(XYFitCurve); if (column != d->yErrorColumn) { exec(new XYFitCurveSetYErrorColumnCmd(d, column, i18n("%1: assign y-error"))); handleSourceDataChanged(); if (column) { connect(column, SIGNAL(dataChanged(const AbstractColumn*)), this, SLOT(handleSourceDataChanged())); //TODO disconnect on undo } } } // TODO: do not recalculate STD_SETTER_CMD_IMPL_F_S(XYFitCurve, SetFitData, XYFitCurve::FitData, fitData, recalculate); void XYFitCurve::setFitData(const XYFitCurve::FitData& fitData) { Q_D(XYFitCurve); exec(new XYFitCurveSetFitDataCmd(d, fitData, i18n("%1: set fit options and perform the fit"))); } //############################################################################## //######################### Private implementation ############################# //############################################################################## XYFitCurvePrivate::XYFitCurvePrivate(XYFitCurve* owner) : XYAnalysisCurvePrivate(owner), xErrorColumn(nullptr), yErrorColumn(nullptr), residualsColumn(nullptr), residualsVector(nullptr), q(owner) { } XYFitCurvePrivate::~XYFitCurvePrivate() { //no need to delete xColumn and yColumn, they are deleted //when the parent aspect is removed } // data structure to pass parameter to fit functions struct data { size_t n; //number of data points double* x; //pointer to the vector with x-data values double* y; //pointer to the vector with y-data values double* weight; //pointer to the vector with weight values nsl_fit_model_category modelCategory; int modelType; int degree; QString* func; // string containing the definition of the model/function QStringList* paramNames; double* paramMin; // lower parameter limits double* paramMax; // upper parameter limits bool* paramFixed; // parameter fixed? }; /*! * \param paramValues vector containing current values of the fit parameters * \param params * \param f vector with the weighted residuals weight[i]*(Yi - y[i]) */ int func_f(const gsl_vector* paramValues, void* params, gsl_vector* f) { //DEBUG("func_f"); size_t n = ((struct data*)params)->n; double* x = ((struct data*)params)->x; double* y = ((struct data*)params)->y; double* weight = ((struct data*)params)->weight; nsl_fit_model_category modelCategory = ((struct data*)params)->modelCategory; unsigned int modelType = ((struct data*)params)->modelType; QByteArray funcba = ((struct data*)params)->func->toLatin1(); // a local byte array is needed! const char *func = funcba.constData(); // function to evaluate QStringList* paramNames = ((struct data*)params)->paramNames; double *min = ((struct data*)params)->paramMin; double *max = ((struct data*)params)->paramMax; // set current values of the parameters for (int i = 0; i < paramNames->size(); i++) { double v = gsl_vector_get(paramValues, (size_t)i); // bound values if limits are set QByteArray paramnameba = paramNames->at(i).toLatin1(); assign_variable(paramnameba.constData(), nsl_fit_map_bound(v, min[i], max[i])); QDEBUG("Parameter"<n; double* xVector = ((struct data*)params)->x; double* weight = ((struct data*)params)->weight; nsl_fit_model_category modelCategory = ((struct data*)params)->modelCategory; unsigned int modelType = ((struct data*)params)->modelType; unsigned int degree = ((struct data*)params)->degree; QStringList* paramNames = ((struct data*)params)->paramNames; double *min = ((struct data*)params)->paramMin; double *max = ((struct data*)params)->paramMax; bool *fixed = ((struct data*)params)->paramFixed; // calculate the Jacobian matrix: // Jacobian matrix J(i,j) = df_i / dx_j // where f_i = w_i*(Y_i - y_i), // Y_i = model and the x_j are the parameters double x; switch (modelCategory) { case nsl_fit_model_basic: switch (modelType) { case nsl_fit_model_polynomial: // Y(x) = c0 + c1*x + ... + cn*x^n for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < (unsigned int)paramNames->size(); ++j) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_polynomial_param_deriv(x, j, weight[i])); } } break; case nsl_fit_model_power: // Y(x) = a*x^b or Y(x) = a + b*x^c. if (degree == 1) { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_power1_param_deriv(j, x, a, b, weight[i])); } } } else if (degree == 2) { const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double c = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_power2_param_deriv(j, x, b, c, weight[i])); } } } break; case nsl_fit_model_exponential: { // Y(x) = a*exp(b*x) + c*exp(d*x) + ... double *p = new double[2*degree]; for (unsigned int i = 0; i < 2*degree; i++) p[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, i), min[i], max[i]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2*degree; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exponentialn_param_deriv(j, x, p, weight[i])); } } delete[] p; break; } case nsl_fit_model_inverse_exponential: { // Y(x) = a*(1-exp(b*x))+c const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_inverse_exponential_param_deriv(j, x, a, b, weight[i])); } } break; } case nsl_fit_model_fourier: { // Y(x) = a0 + (a1*cos(w*x) + b1*sin(w*x)) + ... + (an*cos(n*w*x) + bn*sin(n*w*x) //parameters: w, a0, a1, b1, ... an, bn double* a = new double[degree]; double* b = new double[degree]; double w = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); a[0] = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); b[0] = 0; for (unsigned int i = 1; i < degree; ++i) { a[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, 2*i), min[2*i], max[2*i]); b[i] = nsl_fit_map_bound(gsl_vector_get(paramValues, 2*i+1), min[2*i+1], max[2*i+1]); } for (size_t i = 0; i < n; i++) { x = xVector[i]; double wd = 0; //first derivative with respect to the w parameter for (unsigned int j = 1; j < degree; ++j) { wd += -a[j]*j*x*sin(j*w*x) + b[j]*j*x*cos(j*w*x); } gsl_matrix_set(J, i, 0, weight[i]*wd); gsl_matrix_set(J, i, 1, weight[i]); for (unsigned int j = 1; j <= degree; ++j) { gsl_matrix_set(J, (size_t)i, (size_t)(2*j), nsl_fit_model_fourier_param_deriv(0, j, x, w, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(2*j+1), nsl_fit_model_fourier_param_deriv(1, j, x, w, weight[i])); } for (unsigned int j = 0; j <= 2*degree+1; j++) if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } delete[] a; delete[] b; break; } } break; case nsl_fit_model_peak: switch (modelType) { case nsl_fit_model_gaussian: case nsl_fit_model_lorentz: case nsl_fit_model_sech: case nsl_fit_model_logistic: for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < degree; ++j) { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 3*j), min[3*j], max[3*j]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 3*j+1), min[3*j+1], max[3*j+1]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3*j+2), min[3*j+2], max[3*j+2]); switch (modelType) { case nsl_fit_model_gaussian: gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_gaussian_param_deriv(0, x, a, s, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+1), nsl_fit_model_gaussian_param_deriv(1, x, a, s, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+2), nsl_fit_model_gaussian_param_deriv(2, x, a, s, mu, weight[i])); break; case nsl_fit_model_lorentz: // a,s,t gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_lorentz_param_deriv(0, x, a, s, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+1), nsl_fit_model_lorentz_param_deriv(1, x, a, s, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+2), nsl_fit_model_lorentz_param_deriv(2, x, a, s, mu, weight[i])); break; case nsl_fit_model_sech: gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_sech_param_deriv(0, x, a, s, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+1), nsl_fit_model_sech_param_deriv(1, x, a, s, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+2), nsl_fit_model_sech_param_deriv(2, x, a, s, mu, weight[i])); break; case nsl_fit_model_logistic: gsl_matrix_set(J, (size_t)i, (size_t)(3*j), nsl_fit_model_logistic_param_deriv(0, x, a, s, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+1), nsl_fit_model_logistic_param_deriv(1, x, a, s, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(3*j+2), nsl_fit_model_logistic_param_deriv(2, x, a, s, mu, weight[i])); break; } } for (unsigned int j = 0; j < 3*degree; j++) if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } break; case nsl_fit_model_voigt: for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < degree; ++j) { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j), min[4*j], max[4*j]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+1), min[4*j+1], max[4*j+1]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+2), min[4*j+2], max[4*j+2]); const double g = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+3), min[4*j+3], max[4*j+3]); gsl_matrix_set(J, (size_t)i, (size_t)(4*j), nsl_fit_model_voigt_param_deriv(0, x, a, mu, s, g, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(4*j+1), nsl_fit_model_voigt_param_deriv(1, x, a, mu, s, g, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(4*j+2), nsl_fit_model_voigt_param_deriv(2, x, a, mu, s, g, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(4*j+3), nsl_fit_model_voigt_param_deriv(3, x, a, mu, s, g, weight[i])); } for (unsigned int j = 0; j < 4*degree; j++) if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } break; case nsl_fit_model_pseudovoigt1: for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < degree; ++j) { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j), min[4*j], max[4*j]); const double eta = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+1), min[4*j+1], max[4*j+1]); const double w = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+2), min[4*j+2], max[4*j+2]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 4*j+3), min[4*j+3], max[4*j+3]); gsl_matrix_set(J, (size_t)i, (size_t)(4*j), nsl_fit_model_voigt_param_deriv(0, x, a, eta, w, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(4*j+1), nsl_fit_model_voigt_param_deriv(1, x, a, eta, w, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(4*j+2), nsl_fit_model_voigt_param_deriv(2, x, a, eta, w, mu, weight[i])); gsl_matrix_set(J, (size_t)i, (size_t)(4*j+3), nsl_fit_model_voigt_param_deriv(3, x, a, eta, w, mu, weight[i])); } for (unsigned int j = 0; j < 4*degree; j++) if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } break; } break; case nsl_fit_model_growth: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) { gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } else { switch (modelType) { case nsl_fit_model_atan: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_atan_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_tanh: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_tanh_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_algebraic_sigmoid: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_algebraic_sigmoid_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_sigmoid: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_sigmoid_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_erf: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_erf_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_hill: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_hill_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_gompertz: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gompertz_param_deriv(j, x, a, mu, s, weight[i])); break; case nsl_fit_model_gudermann: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gudermann_param_deriv(j, x, a, mu, s, weight[i])); break; } } } } } break; case nsl_fit_model_distribution: switch (modelType) { case nsl_sf_stats_gaussian: case nsl_sf_stats_exponential: case nsl_sf_stats_laplace: case nsl_sf_stats_cauchy_lorentz: case nsl_sf_stats_rayleigh_tail: case nsl_sf_stats_lognormal: case nsl_sf_stats_logistic: case nsl_sf_stats_sech: case nsl_sf_stats_levy: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { switch (modelType) { case nsl_sf_stats_gaussian: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gaussian_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_exponential: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exponential_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_laplace: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_laplace_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_cauchy_lorentz: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_lorentz_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_rayleigh_tail: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_rayleigh_tail_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_lognormal: if (x > 0) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_lognormal_param_deriv(j, x, a, s, mu, weight[i])); else gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); break; case nsl_sf_stats_logistic: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_logistic_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_sech: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_sech_dist_param_deriv(j, x, a, s, mu, weight[i])); break; case nsl_sf_stats_levy: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_levy_param_deriv(j, x, a, s, mu, weight[i])); break; } } } } break; } case nsl_sf_stats_gaussian_tail: { const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gaussian_tail_param_deriv(j, x, A, s, a, mu, weight[i])); } } break; } case nsl_sf_stats_exponential_power: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_exp_pow_param_deriv(j, x, a, s, b, mu, weight[i])); } } break; } case nsl_sf_stats_rayleigh: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_rayleigh_param_deriv(j, x, a, s, weight[i])); } } break; } case nsl_sf_stats_gamma: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double k = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double t = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gamma_param_deriv(j, x, a, k, t, weight[i])); } } break; } case nsl_sf_stats_flat: { const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_flat_param_deriv(j, x, A, b, a, weight[i])); } } break; } case nsl_sf_stats_chi_squared: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double nu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_chi_square_param_deriv(j, x, a, nu, weight[i])); } } break; } case nsl_sf_stats_tdist: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double nu = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_students_t_param_deriv(j, x, a, nu, weight[i])); } } break; } case nsl_sf_stats_fdist: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double n1 = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double n2 = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_fdist_param_deriv(j, x, a, n1, n2, weight[i])); } } break; } case nsl_sf_stats_beta: case nsl_sf_stats_pareto: { const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { switch (modelType) { case nsl_sf_stats_beta: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_beta_param_deriv(j, x, A, a, b, weight[i])); break; case nsl_sf_stats_pareto: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_pareto_param_deriv(j, x, A, a, b, weight[i])); break; } } } } break; } case nsl_sf_stats_weibull: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double k = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double l = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { if (x > 0) gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_weibull_param_deriv(j, x, a, k, l, mu, weight[i])); else gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); } } } break; } case nsl_sf_stats_gumbel1: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gumbel1_param_deriv(j, x, a, s, mu, b, weight[i])); } } break; } case nsl_sf_stats_gumbel2: { const double A = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double b = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_gumbel2_param_deriv(j, x, A, a, b, mu, weight[i])); } } break; } case nsl_sf_stats_poisson: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double l = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_poisson_param_deriv(j, x, a, l, weight[i])); } } break; } case nsl_sf_stats_maxwell_boltzmann: { // Y(x) = a*sqrt(2/pi) * x^2/s^3 * exp(-(x/s)^2/2) const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_maxwell_param_deriv(j, x, a, s, weight[i])); } } break; } case nsl_sf_stats_frechet: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double g = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double s = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double mu = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_frechet_param_deriv(j, x, a, g, s, mu, weight[i])); } } break; } case nsl_sf_stats_landau: { // const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); for (size_t i = 0; i < n; i++) { x = xVector[i]; if (fixed[0]) gsl_matrix_set(J, (size_t)i, 0, 0.); else gsl_matrix_set(J, (size_t)i, 0, nsl_fit_model_landau_param_deriv(0, x, weight[i])); } break; } case nsl_sf_stats_binomial: case nsl_sf_stats_negative_binomial: case nsl_sf_stats_pascal: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double p = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double N = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 3; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { switch (modelType) { case nsl_sf_stats_binomial: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_binomial_param_deriv(j, x, a, p, N, weight[i])); break; case nsl_sf_stats_negative_binomial: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_negative_binomial_param_deriv(j, x, a, p, N, weight[i])); break; case nsl_sf_stats_pascal: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_pascal_param_deriv(j, x, a, p, N, weight[i])); break; } } } } break; } case nsl_sf_stats_geometric: case nsl_sf_stats_logarithmic: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double p = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 2; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else { switch (modelType) { case nsl_sf_stats_geometric: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_geometric_param_deriv(j, x, a, p, weight[i])); break; case nsl_sf_stats_logarithmic: gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_logarithmic_param_deriv(j, x, a, p, weight[i])); break; } } } } break; } case nsl_sf_stats_hypergeometric: { const double a = nsl_fit_map_bound(gsl_vector_get(paramValues, 0), min[0], max[0]); const double n1 = nsl_fit_map_bound(gsl_vector_get(paramValues, 1), min[1], max[1]); const double n2 = nsl_fit_map_bound(gsl_vector_get(paramValues, 2), min[2], max[2]); const double t = nsl_fit_map_bound(gsl_vector_get(paramValues, 3), min[3], max[3]); for (size_t i = 0; i < n; i++) { x = xVector[i]; for (unsigned int j = 0; j < 4; j++) { if (fixed[j]) gsl_matrix_set(J, (size_t)i, (size_t)j, 0.); else gsl_matrix_set(J, (size_t)i, (size_t)j, nsl_fit_model_hypergeometric_param_deriv(j, x, a, n1, n2, t, weight[i])); } } break; } // unused distributions case nsl_sf_stats_levy_alpha_stable: case nsl_sf_stats_levy_skew_alpha_stable: case nsl_sf_stats_bernoulli: break; } break; case nsl_fit_model_custom: QByteArray funcba = ((struct data*)params)->func->toLatin1(); const char* func = funcba.data(); QByteArray nameba; double value; const unsigned int np = paramNames->size(); for (size_t i = 0; i < n; i++) { x = xVector[i]; assign_variable("x", x); for (unsigned int j = 0; j < np; j++) { for (unsigned int k = 0; k < np; k++) { if (k != j) { nameba = paramNames->at(k).toLatin1(); value = nsl_fit_map_bound(gsl_vector_get(paramValues, k), min[k], max[k]); assign_variable(nameba.data(), value); } } nameba = paramNames->at(j).toLatin1(); const char *name = nameba.data(); value = nsl_fit_map_bound(gsl_vector_get(paramValues, j), min[j], max[j]); assign_variable(name, value); const double f_p = parse(func); double eps = 1.e-9; if (std::abs(f_p) > 0) eps *= std::abs(f_p); // scale step size with function value value += eps; assign_variable(name, value); const double f_pdp = parse(func); // DEBUG("evaluate deriv"<* >(xColumn->data()); yVector = static_cast* >(yColumn->data()); residualsVector = static_cast* >(residualsColumn->data()); xColumn->setHidden(true); q->addChild(xColumn); yColumn->setHidden(true); q->addChild(yColumn); q->addChild(residualsColumn); q->setUndoAware(false); q->setXColumn(xColumn); q->setYColumn(yColumn); q->setUndoAware(true); } else { xVector->clear(); yVector->clear(); residualsVector->clear(); } } void XYFitCurvePrivate::recalculate() { DEBUG("XYFitCurvePrivate::recalculate()"); QElapsedTimer timer; timer.start(); // prepare source data columns const AbstractColumn* tmpXDataColumn = 0; const AbstractColumn* tmpYDataColumn = 0; if (dataSourceType == XYAnalysisCurve::DataSourceSpreadsheet) { DEBUG(" spreadsheet columns as data source"); tmpXDataColumn = xDataColumn; tmpYDataColumn = yDataColumn; } else { DEBUG(" curve columns as data source"); tmpXDataColumn = dataSourceCurve->xColumn(); tmpYDataColumn = dataSourceCurve->yColumn(); } if (!tmpXDataColumn || !tmpYDataColumn) { DEBUG("ERROR: Preparing source data columns failed!"); return; } prepareResultColumns(); // clear the previous result fitResult = XYFitCurve::FitResult(); //fit settings const unsigned int maxIters = fitData.maxIterations; //maximal number of iterations const double delta = fitData.eps; //fit tolerance const unsigned int np = fitData.paramNames.size(); //number of fit parameters if (np == 0) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("Model has no parameters."); emit q->dataChanged(); sourceDataChangedSinceLastRecalc = false; return; } //determine the data source columns if (!tmpXDataColumn || !tmpYDataColumn) { emit q->dataChanged(); sourceDataChangedSinceLastRecalc = false; return; } //check column sizes if (tmpXDataColumn->rowCount() != tmpYDataColumn->rowCount()) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("Number of x and y data points must be equal."); emit q->dataChanged(); sourceDataChangedSinceLastRecalc = false; return; } if (yErrorColumn) { if (yErrorColumn->rowCount() < tmpXDataColumn->rowCount()) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("Not sufficient weight data points provided."); emit q->dataChanged(); sourceDataChangedSinceLastRecalc = false; return; } } //copy all valid data point for the fit to temporary vectors QVector xdataVector; QVector ydataVector; QVector xerrorVector; QVector yerrorVector; double xmin, xmax; if (fitData.autoRange) { xmin = tmpXDataColumn->minimum(); xmax = tmpXDataColumn->maximum(); } else { xmin = fitData.fitRange.first(); xmax = fitData.fitRange.last(); } DEBUG("fit range = " << xmin << " .. " << xmax); for (int row = 0; row < tmpXDataColumn->rowCount(); ++row) { //only copy those data where _all_ values (for x and y and errors, if given) are valid if (!std::isnan(tmpXDataColumn->valueAt(row)) && !std::isnan(tmpYDataColumn->valueAt(row)) && !tmpXDataColumn->isMasked(row) && !tmpYDataColumn->isMasked(row)) { // only when inside given range if (tmpXDataColumn->valueAt(row) >= xmin && tmpXDataColumn->valueAt(row) <= xmax) { if (dataSourceType == XYAnalysisCurve::DataSourceCurve || (!xErrorColumn && !yErrorColumn) || !fitData.useDataErrors) { // x-y xdataVector.append(tmpXDataColumn->valueAt(row)); ydataVector.append(tmpYDataColumn->valueAt(row)); } else if (!xErrorColumn) { // x-y-dy if (!std::isnan(yErrorColumn->valueAt(row))) { xdataVector.append(tmpXDataColumn->valueAt(row)); ydataVector.append(tmpYDataColumn->valueAt(row)); yerrorVector.append(yErrorColumn->valueAt(row)); } } else { // x-y-dx-dy if (!std::isnan(xErrorColumn->valueAt(row)) && !std::isnan(yErrorColumn->valueAt(row))) { xdataVector.append(tmpXDataColumn->valueAt(row)); ydataVector.append(tmpYDataColumn->valueAt(row)); xerrorVector.append(xErrorColumn->valueAt(row)); yerrorVector.append(yErrorColumn->valueAt(row)); } } } } } //number of data points to fit const size_t n = xdataVector.size(); DEBUG("number of data points: " << n); if (n == 0) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("No data points available."); emit q->dataChanged(); sourceDataChangedSinceLastRecalc = false; return; } if (n < np) { fitResult.available = true; fitResult.valid = false; fitResult.status = i18n("The number of data points (%1) must be greater than or equal to the number of parameters (%2).", n, np); emit q->dataChanged(); sourceDataChangedSinceLastRecalc = false; return; } double* xdata = xdataVector.data(); double* ydata = ydataVector.data(); double* xerror = xerrorVector.data(); // size may be 0 double* yerror = yerrorVector.data(); // size may be 0 DEBUG("x errors: " << xerrorVector.size()); DEBUG("y errors: " << yerrorVector.size()); double* weight = new double[n]; for (size_t i = 0; i < n; i++) weight[i] = 1.; switch (fitData.yWeightsType) { case nsl_fit_weight_no: case nsl_fit_weight_statistical_fit: case nsl_fit_weight_relative_fit: break; case nsl_fit_weight_instrumental: // yerror are sigmas for(int i = 0; i < (int)n; i++) if (i < yerrorVector.size()) weight[i] = 1./gsl_pow_2(yerror[i]); break; case nsl_fit_weight_direct: // yerror are weights for(int i = 0; i < (int)n; i++) if (i < yerrorVector.size()) weight[i] = yerror[i]; break; case nsl_fit_weight_inverse: // yerror are inverse weights for(int i = 0; i < (int)n; i++) if (i < yerrorVector.size()) weight[i] = 1./yerror[i]; break; case nsl_fit_weight_statistical: for (int i = 0; i < (int)n; i++) weight[i] = 1./ydata[i]; break; case nsl_fit_weight_relative: for (int i = 0; i < (int)n; i++) weight[i] = 1./gsl_pow_2(ydata[i]); break; } /////////////////////// GSL >= 2 has a complete new interface! But the old one is still supported. /////////////////////////// // GSL >= 2 : "the 'fdf' field of gsl_multifit_function_fdf is now deprecated and does not need to be specified for nonlinear least squares problems" unsigned int nf = 0; // number of fixed parameter for (unsigned int i = 0; i < np; i++) { const bool fixed = fitData.paramFixed.data()[i]; if (fixed) nf++; DEBUG("parameter " << i << " fixed: " << fixed); } //function to fit gsl_multifit_function_fdf f; DEBUG("model = " << fitData.model.toStdString()); struct data params = {n, xdata, ydata, weight, fitData.modelCategory, fitData.modelType, fitData.degree, &fitData.model, &fitData.paramNames, fitData.paramLowerLimits.data(), fitData.paramUpperLimits.data(), fitData.paramFixed.data()}; f.f = &func_f; f.df = &func_df; f.fdf = &func_fdf; f.n = n; f.p = np; f.params = ¶ms; DEBUG("initialize the derivative solver (using Levenberg-Marquardt robust solver)"); const gsl_multifit_fdfsolver_type* T = gsl_multifit_fdfsolver_lmsder; gsl_multifit_fdfsolver* s = gsl_multifit_fdfsolver_alloc(T, n, np); DEBUG("set start values"); double* x_init = fitData.paramStartValues.data(); double* x_min = fitData.paramLowerLimits.data(); double* x_max = fitData.paramUpperLimits.data(); DEBUG("scale start values if limits are set"); for (unsigned int i = 0; i < np; i++) x_init[i] = nsl_fit_map_unbound(x_init[i], x_min[i], x_max[i]); DEBUG(" DONE"); gsl_vector_view x = gsl_vector_view_array(x_init, np); DEBUG("Turning off GSL error handler to avoid overflow/underflow"); gsl_set_error_handler_off(); DEBUG("Initialize solver with function f and initial guess x"); gsl_multifit_fdfsolver_set(s, &f, &x.vector); DEBUG("Iterate ..."); int status; unsigned int iter = 0; fitResult.solverOutput.clear(); writeSolverState(s); do { iter++; // update weights for Y-depending weights (using function values from residuals) if (fitData.yWeightsType == nsl_fit_weight_statistical_fit) { for (size_t i = 0; i < n; i++) weight[i] = 1./(gsl_vector_get(s->f, i)/sqrt(weight[i]) + ydata[i]); // 1/Y_i } else if (fitData.yWeightsType == nsl_fit_weight_relative_fit) { for (size_t i = 0; i < n; i++) weight[i] = 1./gsl_pow_2(gsl_vector_get(s->f, i)/sqrt(weight[i]) + ydata[i]); // 1/Y_i^2 } status = gsl_multifit_fdfsolver_iterate(s); writeSolverState(s); if (status) { DEBUG("iter " << iter << ", status = " << gsl_strerror(status)); break; } status = gsl_multifit_test_delta(s->dx, s->x, delta, delta); } while (status == GSL_CONTINUE && iter < maxIters); // second run for x-error fitting if (xerrorVector.size() > 0) { DEBUG("Rerun fit with x errors"); unsigned int iter2 = 0; double chi = 0, chiOld = 0; double *fun = new double[n]; do { iter2++; chiOld = chi; //printf("iter2 = %d\n", iter2); // calculate function from residuals for (size_t i = 0; i < n; i++) fun[i] = gsl_vector_get(s->f, i) * 1./sqrt(weight[i]) + ydata[i]; // calculate weight[i] for (size_t i = 0; i < n; i++) { // calculate df[i] size_t index = i-1; if (i == 0) index = i; if (i == n-1) index = i-2; double df = (fun[index+1] - fun[index])/(xdata[index+1] - xdata[index]); //printf("df = %g\n", df); double sigmasq = 1.; switch (fitData.xWeightsType) { // x-error type: f'(x)^2*s_x^2 = f'(x)/w_x case nsl_fit_weight_no: break; case nsl_fit_weight_direct: // xerror = w_x sigmasq = df*df/xerror[i]; break; case nsl_fit_weight_instrumental: // xerror = s_x sigmasq = df*df*xerror[i]*xerror[i]; break; case nsl_fit_weight_inverse: // xerror = 1/w_x = s_x^2 sigmasq = df*df*xerror[i]; break; case nsl_fit_weight_statistical: // s_x^2 = 1/w_x = x sigmasq = xdata[i]; break; case nsl_fit_weight_relative: // s_x^2 = 1/w_x = x^2 sigmasq = xdata[i]*xdata[i]; break; case nsl_fit_weight_statistical_fit: // unused case nsl_fit_weight_relative_fit: break; } if (yerrorVector.size() > 0) { switch (fitData.yWeightsType) { // y-error types: s_y^2 = 1/w_y case nsl_fit_weight_no: break; case nsl_fit_weight_direct: // yerror = w_y sigmasq += 1./yerror[i]; break; case nsl_fit_weight_instrumental: // yerror = s_y sigmasq += yerror[i]*yerror[i]; break; case nsl_fit_weight_inverse: // yerror = 1/w_y sigmasq += yerror[i]; break; case nsl_fit_weight_statistical: // unused case nsl_fit_weight_relative: break; case nsl_fit_weight_statistical_fit: // s_y^2 = 1/w_y = Y_i sigmasq += fun[i]; break; case nsl_fit_weight_relative_fit: // s_y^2 = 1/w_y = Y_i^2 sigmasq += fun[i]*fun[i]; break; } } //printf ("sigma[%d] = %g\n", i, sqrt(sigmasq)); weight[i] = 1./sigmasq; } // update weights gsl_multifit_fdfsolver_set(s, &f, &x.vector); do { // fit iter++; writeSolverState(s); status = gsl_multifit_fdfsolver_iterate (s); //printf ("status = %s\n", gsl_strerror (status)); if (status) { DEBUG("iter " << iter << ", status = " << gsl_strerror(status)); break; } status = gsl_multifit_test_delta(s->dx, s->x, delta, delta); } while (status == GSL_CONTINUE && iter < maxIters); chi = gsl_blas_dnrm2(s->f); //printf("chi = %.12g (dchi = %g)\n", chi, fabs(chi-chiOld)); } while (iter2 < maxIters && fabs(chi-chiOld) > fitData.eps); delete[] fun; } delete[] weight; // unscale start parameter for (unsigned int i = 0; i < np; i++) x_init[i] = nsl_fit_map_bound(x_init[i], x_min[i], x_max[i]); //get the covariance matrix //TODO: scale the Jacobian when limits are used before constructing the covar matrix? gsl_matrix* covar = gsl_matrix_alloc(np, np); #if GSL_MAJOR_VERSION >= 2 // the Jacobian is not part of the solver anymore gsl_matrix *J = gsl_matrix_alloc(s->fdf->n, s->fdf->p); gsl_multifit_fdfsolver_jac(s, J); gsl_multifit_covar(J, 0.0, covar); #else gsl_multifit_covar(s->J, 0.0, covar); #endif //write the result fitResult.available = true; fitResult.valid = true; fitResult.status = gslErrorToString(status); fitResult.iterations = iter; fitResult.dof = n - (np - nf); // samples - (parameter - fixed parameter) //gsl_blas_dnrm2() - computes the Euclidian norm (||r||_2 = \sqrt {\sum r_i^2}) of the vector with the elements weight[i]*(Yi - y[i]) //gsl_blas_dasum() - computes the absolute sum \sum |r_i| of the elements of the vector with the elements weight[i]*(Yi - y[i]) fitResult.sse = gsl_pow_2(gsl_blas_dnrm2(s->f)); if (fitResult.dof != 0) { fitResult.rms = fitResult.sse/fitResult.dof; fitResult.rsd = sqrt(fitResult.rms); } fitResult.mse = fitResult.sse/n; fitResult.rmse = sqrt(fitResult.mse); fitResult.mae = gsl_blas_dasum(s->f)/n; // SST needed for coefficient of determination, R-squared fitResult.sst = gsl_stats_tss(ydata, 1, n); // for a linear model without intercept R-squared is calculated differently // see https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-does-summary_0028_0029-report-strange-results-for-the-R_005e2-estimate-when-I-fit-a-linear-model-with-no-intercept_003f if (fitData.modelCategory == nsl_fit_model_basic && fitData.modelType == nsl_fit_model_polynomial && fitData.degree == 1 && x_init[0] == 0) { DEBUG("Using alternative R^2 for linear model without intercept"); fitResult.sst = gsl_stats_tss_m(ydata, 1, n, 0); } if (fitResult.sst < fitResult.sse) { DEBUG("Using alternative R^2 since R^2 would be negative (probably custom model without intercept)"); fitResult.sst = gsl_stats_tss_m(ydata, 1, n, 0); } fitResult.rsquare = nsl_stats_rsquare(fitResult.sse, fitResult.sst); fitResult.rsquareAdj = nsl_stats_rsquareAdj(fitResult.rsquare, np, fitResult.dof, 1); fitResult.chisq_p = nsl_stats_chisq_p(fitResult.sse, fitResult.dof); fitResult.fdist_F = nsl_stats_fdist_F(fitResult.sst, fitResult.rms, np, 1); fitResult.fdist_p = nsl_stats_fdist_p(fitResult.fdist_F, np, fitResult.dof); fitResult.logLik = nsl_stats_logLik(fitResult.sse, n); fitResult.aic = nsl_stats_aic(fitResult.sse, n, np, 1); fitResult.bic = nsl_stats_bic(fitResult.sse, n, np, 1); //parameter values // GSL: const double c = GSL_MAX_DBL(1., sqrt(fitResult.rms)); // increase error for poor fit // NIST: const double c = sqrt(fitResult.rms); // increase error for poor fit, decrease for good fit const double c = sqrt(fitResult.rms); fitResult.paramValues.resize(np); fitResult.errorValues.resize(np); fitResult.tdist_tValues.resize(np); fitResult.tdist_pValues.resize(np); fitResult.tdist_marginValues.resize(np); for (unsigned int i = 0; i < np; i++) { // scale resulting values if they are bounded fitResult.paramValues[i] = nsl_fit_map_bound(gsl_vector_get(s->x, i), x_min[i], x_max[i]); // use results as start values if desired if (fitData.useResults) { fitData.paramStartValues.data()[i] = fitResult.paramValues[i]; DEBUG("saving parameter " << i << ": " << fitResult.paramValues[i] << ' ' << fitData.paramStartValues.data()[i]); } fitResult.errorValues[i] = c*sqrt(gsl_matrix_get(covar, i, i)); fitResult.tdist_tValues[i] = nsl_stats_tdist_t(fitResult.paramValues.at(i), fitResult.errorValues.at(i)); fitResult.tdist_pValues[i] = nsl_stats_tdist_p(fitResult.tdist_tValues.at(i), fitResult.dof); fitResult.tdist_marginValues[i] = nsl_stats_tdist_margin(0.05, fitResult.dof, fitResult.errorValues.at(i)); } // fill residuals vector. To get residuals on the correct x values, fill the rest with zeros. residualsVector->resize(tmpXDataColumn->rowCount()); if (fitData.autoRange) { // evaluate full range of residuals xVector->resize(tmpXDataColumn->rowCount()); for (int i = 0; i < tmpXDataColumn->rowCount(); i++) (*xVector)[i] = tmpXDataColumn->valueAt(i); ExpressionParser* parser = ExpressionParser::getInstance(); bool rc = parser->evaluateCartesian(fitData.model, xVector, residualsVector, fitData.paramNames, fitResult.paramValues); for (int i = 0; i < tmpXDataColumn->rowCount(); i++) (*residualsVector)[i] = tmpYDataColumn->valueAt(i) - (*residualsVector)[i]; if (!rc) residualsVector->clear(); } else { // only selected range size_t j = 0; for (int i = 0; i < tmpXDataColumn->rowCount(); i++) { if (tmpXDataColumn->valueAt(i) >= xmin && tmpXDataColumn->valueAt(i) <= xmax) residualsVector->data()[i] = - gsl_vector_get(s->f, j++); else // outside range residualsVector->data()[i] = 0; } } residualsColumn->setChanged(); //free resources gsl_multifit_fdfsolver_free(s); gsl_matrix_free(covar); //calculate the fit function (vectors) evaluate(); fitResult.elapsedTime = timer.elapsed(); //redraw the curve emit q->dataChanged(); sourceDataChangedSinceLastRecalc = false; } /* evaluate fit function */ void XYFitCurvePrivate::evaluate(bool preview) { DEBUG("XYFitCurvePrivate::evaluate() preview = " << preview); // prepare source data columns const AbstractColumn* tmpXDataColumn = 0; if (dataSourceType == XYAnalysisCurve::DataSourceSpreadsheet) { DEBUG(" spreadsheet columns as data source"); tmpXDataColumn = xDataColumn; } else { DEBUG(" curve columns as data source"); tmpXDataColumn = dataSourceCurve->xColumn(); } if (!tmpXDataColumn) { DEBUG("ERROR: Preparing source data column failed!"); return; } prepareResultColumns(); if (!xVector || !yVector) { DEBUG(" xVector or yVector not defined!"); return; } ExpressionParser* parser = ExpressionParser::getInstance(); double xmin, xmax; if (fitData.autoEvalRange) { // evaluate fit on full data range xmin = tmpXDataColumn->minimum(); xmax = tmpXDataColumn->maximum(); } else { // use given range for evaluation xmin = fitData.evalRange.first(); xmax = fitData.evalRange.last(); } DEBUG(" eval range = " << xmin << " .. " << xmax); xVector->resize((int)fitData.evaluatedPoints); yVector->resize((int)fitData.evaluatedPoints); DEBUG(" vector size = " << xVector->size()); QVector paramValues = fitResult.paramValues; if (preview) // results not available yet paramValues = fitData.paramStartValues; // Debug /* if (paramValues.size() == 0) DEBUG(" ERROR: No parameter defined!"); for (auto value: paramValues) DEBUG(" param value = " << value); */ bool rc = parser->evaluateCartesian(fitData.model, QString::number(xmin), QString::number(xmax), (int)fitData.evaluatedPoints, xVector, yVector, fitData.paramNames, paramValues); if (!rc) { xVector->clear(); yVector->clear(); residualsVector->clear(); } // Debug /* DEBUG(" x | y"); for (int i = 0; i < qMin(10, xVector->size()); i++) DEBUG(" " << (*xVector)[i] << " | " << (*yVector)[i]); */ // TODO: do we weed to do something to make preview work? // this should be already done by dataChanged() q->retransform(); // PREVIEW redraw the curve (this breaks context menu fit!) if (preview) { // emit q->dataChanged(); // sourceDataChangedSinceLastRecalc = false; } } /*! * writes out the current state of the solver \c s */ void XYFitCurvePrivate::writeSolverState(gsl_multifit_fdfsolver* s) { QString state; //current parameter values, semicolon separated double* min = fitData.paramLowerLimits.data(); double* max = fitData.paramUpperLimits.data(); for (int i = 0; i < fitData.paramNames.size(); ++i) { const double x = gsl_vector_get(s->x, i); // map parameter if bounded state += QString::number(nsl_fit_map_bound(x, min[i], max[i])) + '\t'; } //current value of the chi2-function state += QString::number(gsl_pow_2(gsl_blas_dnrm2(s->f))); state += ';'; fitResult.solverOutput += state; } //############################################################################## //################## Serialization/Deserialization ########################### //############################################################################## //! Save as XML void XYFitCurve::save(QXmlStreamWriter* writer) const { Q_D(const XYFitCurve); writer->writeStartElement("xyFitCurve"); //write the base class XYAnalysisCurve::save(writer); //write xy-fit-curve specific information //fit data - only save model expression and parameter names for custom model, otherwise they are set in XYFitCurve::initFitData() writer->writeStartElement("fitData"); WRITE_COLUMN(d->xErrorColumn, xErrorColumn); WRITE_COLUMN(d->yErrorColumn, yErrorColumn); writer->writeAttribute("autoRange", QString::number(d->fitData.autoRange)); writer->writeAttribute("fitRangeMin", QString::number(d->fitData.fitRange.first(), 'g', 15)); writer->writeAttribute("fitRangeMax", QString::number(d->fitData.fitRange.last(), 'g', 15)); writer->writeAttribute("modelCategory", QString::number(d->fitData.modelCategory)); writer->writeAttribute("modelType", QString::number(d->fitData.modelType)); writer->writeAttribute("xWeightsType", QString::number(d->fitData.xWeightsType)); writer->writeAttribute("weightsType", QString::number(d->fitData.yWeightsType)); writer->writeAttribute("degree", QString::number(d->fitData.degree)); if (d->fitData.modelCategory == nsl_fit_model_custom) writer->writeAttribute("model", d->fitData.model); writer->writeAttribute("maxIterations", QString::number(d->fitData.maxIterations)); writer->writeAttribute("eps", QString::number(d->fitData.eps, 'g', 15)); writer->writeAttribute("evaluatedPoints", QString::number(d->fitData.evaluatedPoints)); writer->writeAttribute("autoEvalRange", QString::number(d->fitData.autoEvalRange)); writer->writeAttribute("useDataErrors", QString::number(d->fitData.useDataErrors)); writer->writeAttribute("useResults", QString::number(d->fitData.useResults)); if (d->fitData.modelCategory == nsl_fit_model_custom) { writer->writeStartElement("paramNames"); foreach (const QString &name, d->fitData.paramNames) writer->writeTextElement("name", name); writer->writeEndElement(); } writer->writeStartElement("paramStartValues"); foreach (const double &value, d->fitData.paramStartValues) writer->writeTextElement("startValue", QString::number(value, 'g', 15)); writer->writeEndElement(); // use 16 digits to handle -DBL_MAX writer->writeStartElement("paramLowerLimits"); foreach (const double &limit, d->fitData.paramLowerLimits) writer->writeTextElement("lowerLimit", QString::number(limit, 'g', 16)); writer->writeEndElement(); // use 16 digits to handle DBL_MAX writer->writeStartElement("paramUpperLimits"); foreach (const double &limit, d->fitData.paramUpperLimits) writer->writeTextElement("upperLimit", QString::number(limit, 'g', 16)); writer->writeEndElement(); writer->writeStartElement("paramFixed"); foreach (const double &fixed, d->fitData.paramFixed) writer->writeTextElement("fixed", QString::number(fixed)); writer->writeEndElement(); writer->writeEndElement(); //"fitData" //fit results (generated columns and goodness of the fit) writer->writeStartElement("fitResult"); writer->writeAttribute("available", QString::number(d->fitResult.available)); writer->writeAttribute("valid", QString::number(d->fitResult.valid)); writer->writeAttribute("status", d->fitResult.status); writer->writeAttribute("iterations", QString::number(d->fitResult.iterations)); writer->writeAttribute("time", QString::number(d->fitResult.elapsedTime)); writer->writeAttribute("dof", QString::number(d->fitResult.dof)); writer->writeAttribute("sse", QString::number(d->fitResult.sse, 'g', 15)); writer->writeAttribute("sst", QString::number(d->fitResult.sst, 'g', 15)); writer->writeAttribute("rms", QString::number(d->fitResult.rms, 'g', 15)); writer->writeAttribute("rsd", QString::number(d->fitResult.rsd, 'g', 15)); writer->writeAttribute("mse", QString::number(d->fitResult.mse, 'g', 15)); writer->writeAttribute("rmse", QString::number(d->fitResult.rmse, 'g', 15)); writer->writeAttribute("mae", QString::number(d->fitResult.mae, 'g', 15)); writer->writeAttribute("rsquare", QString::number(d->fitResult.rsquare, 'g', 15)); writer->writeAttribute("rsquareAdj", QString::number(d->fitResult.rsquareAdj, 'g', 15)); writer->writeAttribute("chisq_p", QString::number(d->fitResult.chisq_p, 'g', 15)); writer->writeAttribute("fdist_F", QString::number(d->fitResult.fdist_F, 'g', 15)); writer->writeAttribute("fdist_p", QString::number(d->fitResult.fdist_p, 'g', 15)); writer->writeAttribute("aic", QString::number(d->fitResult.aic, 'g', 15)); writer->writeAttribute("bic", QString::number(d->fitResult.bic, 'g', 15)); writer->writeAttribute("solverOutput", d->fitResult.solverOutput); writer->writeStartElement("paramValues"); foreach (const double &value, d->fitResult.paramValues) writer->writeTextElement("value", QString::number(value, 'g', 15)); writer->writeEndElement(); writer->writeStartElement("errorValues"); foreach (const double &value, d->fitResult.errorValues) writer->writeTextElement("error", QString::number(value, 'g', 15)); writer->writeEndElement(); writer->writeStartElement("tdist_tValues"); foreach (const double &value, d->fitResult.tdist_tValues) writer->writeTextElement("tdist_t", QString::number(value, 'g', 15)); writer->writeEndElement(); writer->writeStartElement("tdist_pValues"); foreach (const double &value, d->fitResult.tdist_pValues) writer->writeTextElement("tdist_p", QString::number(value, 'g', 15)); writer->writeEndElement(); writer->writeStartElement("tdist_marginValues"); foreach (const double &value, d->fitResult.tdist_marginValues) writer->writeTextElement("tdist_margin", QString::number(value, 'g', 15)); writer->writeEndElement(); //save calculated columns if available if (d->xColumn && d->yColumn && d->residualsColumn) { d->xColumn->save(writer); d->yColumn->save(writer); d->residualsColumn->save(writer); } writer->writeEndElement(); //"fitResult" writer->writeEndElement(); //"xyFitCurve" } //! Load from XML bool XYFitCurve::load(XmlStreamReader* reader, bool preview) { Q_D(XYFitCurve); QString attributeWarning = i18n("Attribute '%1' missing or empty, default value is used"); QXmlStreamAttributes attribs; QString str; while (!reader->atEnd()) { reader->readNext(); if (reader->isEndElement() && reader->name() == "xyFitCurve") break; if (!reader->isStartElement()) continue; if (reader->name() == "xyAnalysisCurve") { if ( !XYAnalysisCurve::load(reader, preview) ) return false; } else if (!preview && reader->name() == "fitData") { attribs = reader->attributes(); READ_COLUMN(xErrorColumn); READ_COLUMN(yErrorColumn); READ_INT_VALUE("autoRange", fitData.autoRange, bool); READ_DOUBLE_VALUE("xRangeMin", fitData.fitRange.first()); // old name READ_DOUBLE_VALUE("xRangeMax", fitData.fitRange.last()); // old name READ_DOUBLE_VALUE("fitRangeMin", fitData.fitRange.first()); READ_DOUBLE_VALUE("fitRangeMax", fitData.fitRange.last()); READ_INT_VALUE("modelCategory", fitData.modelCategory, nsl_fit_model_category); READ_INT_VALUE("modelType", fitData.modelType, unsigned int); READ_INT_VALUE("xWeightsType", fitData.xWeightsType, nsl_fit_weight_type); READ_INT_VALUE("weightsType", fitData.yWeightsType, nsl_fit_weight_type); READ_INT_VALUE("degree", fitData.degree, int); if (d->fitData.modelCategory == nsl_fit_model_custom) { READ_STRING_VALUE("model", fitData.model); DEBUG("read model = " << d->fitData.model.toStdString()); } READ_INT_VALUE("maxIterations", fitData.maxIterations, int); READ_DOUBLE_VALUE("eps", fitData.eps); READ_INT_VALUE("fittedPoints", fitData.evaluatedPoints, size_t); // old name READ_INT_VALUE("evaluatedPoints", fitData.evaluatedPoints, size_t); READ_INT_VALUE("evaluateFullRange", fitData.autoEvalRange, bool); // old name READ_INT_VALUE("autoEvalRange", fitData.autoEvalRange, bool); READ_INT_VALUE("useDataErrors", fitData.useDataErrors, bool); READ_INT_VALUE("useResults", fitData.useResults, bool); //set the model expression and the parameter names (can be derived from the saved values for category, type and degree) XYFitCurve::initFitData(d->fitData); } else if (!preview && reader->name() == "name") { // needed for custom model d->fitData.paramNames << reader->readElementText(); } else if (!preview && reader->name() == "startValue") { d->fitData.paramStartValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "fixed") { d->fitData.paramFixed << (bool)reader->readElementText().toInt(); } else if (!preview && reader->name() == "lowerLimit") { bool ok; double x = reader->readElementText().toDouble(&ok); if (ok) // -DBL_MAX results in conversion error d->fitData.paramLowerLimits << x; else d->fitData.paramLowerLimits << -std::numeric_limits::max(); } else if (!preview && reader->name() == "upperLimit") { bool ok; double x = reader->readElementText().toDouble(&ok); if (ok) // DBL_MAX results in conversion error d->fitData.paramUpperLimits << x; else d->fitData.paramUpperLimits << std::numeric_limits::max(); } else if (!preview && reader->name() == "value") { d->fitResult.paramValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "error") { d->fitResult.errorValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "tdist_t") { d->fitResult.tdist_tValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "tdist_p") { d->fitResult.tdist_pValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "tdist_margin") { d->fitResult.tdist_marginValues << reader->readElementText().toDouble(); } else if (!preview && reader->name() == "fitResult") { attribs = reader->attributes(); READ_INT_VALUE("available", fitResult.available, int); READ_INT_VALUE("valid", fitResult.valid, int); READ_STRING_VALUE("status", fitResult.status); READ_INT_VALUE("iterations", fitResult.iterations, int); READ_INT_VALUE("time", fitResult.elapsedTime, int); READ_DOUBLE_VALUE("dof", fitResult.dof); READ_DOUBLE_VALUE("sse", fitResult.sse); READ_DOUBLE_VALUE("sst", fitResult.sst); READ_DOUBLE_VALUE("rms", fitResult.rms); READ_DOUBLE_VALUE("rsd", fitResult.rsd); READ_DOUBLE_VALUE("mse", fitResult.mse); READ_DOUBLE_VALUE("rmse", fitResult.rmse); READ_DOUBLE_VALUE("mae", fitResult.mae); READ_DOUBLE_VALUE("rsquare", fitResult.rsquare); READ_DOUBLE_VALUE("rsquareAdj", fitResult.rsquareAdj); READ_DOUBLE_VALUE("chisq_p", fitResult.chisq_p); READ_DOUBLE_VALUE("fdist_F", fitResult.fdist_F); READ_DOUBLE_VALUE("fdist_p", fitResult.fdist_p); READ_DOUBLE_VALUE("aic", fitResult.aic); READ_DOUBLE_VALUE("bic", fitResult.bic); READ_STRING_VALUE("solverOutput", fitResult.solverOutput); } else if (reader->name() == "column") { Column* column = new Column("", AbstractColumn::Numeric); if (!column->load(reader, preview)) { delete column; return false; } if (column->name() == "x") d->xColumn = column; else if (column->name() == "y") d->yColumn = column; else if (column->name() == "residuals") d->residualsColumn = column; } } if (preview) return true; // new fit model style (reset model type of old projects) if (d->fitData.modelCategory == nsl_fit_model_basic && d->fitData.modelType >= NSL_FIT_MODEL_BASIC_COUNT) { DEBUG("reset old fit model"); d->fitData.modelType = 0; d->fitData.degree = 1; // reset size of fields not touched by initFitData() d->fitData.paramStartValues.resize(2); d->fitData.paramFixed.resize(2); d->fitResult.paramValues.resize(2); d->fitResult.errorValues.resize(2); d->fitResult.tdist_tValues.resize(2); d->fitResult.tdist_pValues.resize(2); d->fitResult.tdist_marginValues.resize(2); } // not present in old projects if (d->fitResult.tdist_tValues.size() == 0) d->fitResult.tdist_tValues.resize(d->fitResult.paramValues.size()); if (d->fitResult.tdist_pValues.size() == 0) d->fitResult.tdist_pValues.resize(d->fitResult.paramValues.size()); if (d->fitResult.tdist_marginValues.size() == 0) d->fitResult.tdist_marginValues.resize(d->fitResult.paramValues.size()); // wait for data to be read before using the pointers QThreadPool::globalInstance()->waitForDone(); if (d->xColumn && d->yColumn && d->residualsColumn) { d->xColumn->setHidden(true); addChild(d->xColumn); d->yColumn->setHidden(true); addChild(d->yColumn); addChild(d->residualsColumn); d->xVector = static_cast* >(d->xColumn->data()); d->yVector = static_cast* >(d->yColumn->data()); d->residualsVector = static_cast* >(d->residualsColumn->data()); XYCurve::d_ptr->xColumn = d->xColumn; XYCurve::d_ptr->yColumn = d->yColumn; } return true; }