summaryrefslogtreecommitdiff
path: root/src/silx/math/fit/leastsq.py
diff options
context:
space:
mode:
Diffstat (limited to 'src/silx/math/fit/leastsq.py')
-rw-r--r--src/silx/math/fit/leastsq.py901
1 files changed, 901 insertions, 0 deletions
diff --git a/src/silx/math/fit/leastsq.py b/src/silx/math/fit/leastsq.py
new file mode 100644
index 0000000..3df1a35
--- /dev/null
+++ b/src/silx/math/fit/leastsq.py
@@ -0,0 +1,901 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2004-2020 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+"""
+This module implements a Levenberg-Marquardt algorithm with constraints on the
+fitted parameters without introducing any other dependendency than numpy.
+
+If scipy dependency is not an issue, and no constraints are applied to the fitting
+parameters, there is no real gain compared to the use of scipy.optimize.curve_fit
+other than a more conservative calculation of uncertainties on fitted parameters.
+
+This module is a refactored version of PyMca Gefit.py module.
+"""
+__authors__ = ["V.A. Sole"]
+__license__ = "MIT"
+__date__ = "15/05/2017"
+__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
+
+import numpy
+from numpy.linalg import inv
+from numpy.linalg.linalg import LinAlgError
+import time
+import logging
+import copy
+
+_logger = logging.getLogger(__name__)
+
+# codes understood by the routine
+CFREE = 0
+CPOSITIVE = 1
+CQUOTED = 2
+CFIXED = 3
+CFACTOR = 4
+CDELTA = 5
+CSUM = 6
+CIGNORED = 7
+
+def leastsq(model, xdata, ydata, p0, sigma=None,
+ constraints=None, model_deriv=None, epsfcn=None,
+ deltachi=None, full_output=None,
+ check_finite=True,
+ left_derivative=False,
+ max_iter=100):
+ """
+ Use non-linear least squares Levenberg-Marquardt algorithm to fit a function, f, to
+ data with optional constraints on the fitted parameters.
+
+ Assumes ``ydata = f(xdata, *params) + eps``
+
+ :param model: callable
+ The model function, f(x, ...). It must take the independent
+ variable as the first argument and the parameters to fit as
+ separate remaining arguments.
+ The returned value is a one dimensional array of floats.
+
+ :param xdata: An M-length sequence.
+ The independent variable where the data is measured.
+
+ :param ydata: An M-length sequence
+ The dependent data --- nominally f(xdata, ...)
+
+ :param p0: N-length sequence
+ Initial guess for the parameters.
+
+ :param sigma: None or M-length sequence, optional
+ If not None, the uncertainties in the ydata array. These are used as
+ weights in the least-squares problem
+ i.e. minimising ``np.sum( ((f(xdata, *popt) - ydata) / sigma)**2 )``
+ If None, the uncertainties are assumed to be 1
+
+ :param constraints:
+ If provided, it is a 2D sequence of dimension (n_parameters, 3) where,
+ for each parameter denoted by the index i, the meaning is
+
+ - constraints[i][0]
+
+ - 0 - Free (CFREE)
+ - 1 - Positive (CPOSITIVE)
+ - 2 - Quoted (CQUOTED)
+ - 3 - Fixed (CFIXED)
+ - 4 - Factor (CFACTOR)
+ - 5 - Delta (CDELTA)
+ - 6 - Sum (CSUM)
+
+
+ - constraints[i][1]
+
+ - Ignored if constraints[i][0] is 0, 1, 3
+ - Min value of the parameter if constraints[i][0] is CQUOTED
+ - Index of fitted parameter to which it is related
+
+ - constraints[i][2]
+
+ - Ignored if constraints[i][0] is 0, 1, 3
+ - Max value of the parameter if constraints[i][0] is CQUOTED
+ - Factor to apply to related parameter with index constraints[i][1]
+ - Difference with parameter with index constraints[i][1]
+ - Sum obtained when adding parameter with index constraints[i][1]
+ :type constraints: *optional*, None or 2D sequence
+
+ :param model_deriv:
+ None (default) or function providing the derivatives of the fitting function respect to the fitted parameters.
+ It will be called as model_deriv(xdata, parameters, index) where parameters is a sequence with the current
+ values of the fitting parameters, index is the fitting parameter index for which the the derivative has
+ to be provided in the supplied array of xdata points.
+ :type model_deriv: *optional*, None or callable
+
+
+ :param epsfcn: float
+ A variable used in determining a suitable parameter variation when
+ calculating the numerical derivatives (for model_deriv=None).
+ Normally the actual step length will be sqrt(epsfcn)*x
+ Original Gefit module was using epsfcn 1.0e-5 while default value
+ is now numpy.finfo(numpy.float64).eps as in scipy
+ :type epsfcn: *optional*, float
+
+ :param deltachi: float
+ A variable used to control the minimum change in chisq to consider the
+ fitting process not worth to be continued. Default is 0.1 %.
+ :type deltachi: *optional*, float
+
+ :param full_output: bool, optional
+ non-zero to return all optional outputs. The default is None what will give a warning in case
+ of a constrained fit without having set this kweyword.
+
+ :param check_finite: bool, optional
+ If True, check that the input arrays do not contain nans of infs,
+ and raise a ValueError if they do. Setting this parameter to
+ False will ignore input arrays values containing nans.
+ Default is True.
+
+ :param left_derivative:
+ This parameter only has an influence if no derivative function
+ is provided. When True the left and right derivatives of the
+ model will be calculated for each fitted parameters thus leading to
+ the double number of function evaluations. Default is False.
+ Original Gefit module was always using left_derivative as True.
+ :type left_derivative: *optional*, bool
+
+ :param max_iter: Maximum number of iterations (default is 100)
+
+ :return: Returns a tuple of length 2 (or 3 if full_ouput is True) with the content:
+
+ ``popt``: array
+ Optimal values for the parameters so that the sum of the squared error
+ of ``f(xdata, *popt) - ydata`` is minimized
+ ``pcov``: 2d array
+ If no constraints are applied, this array contains the estimated covariance
+ of popt. The diagonal provides the variance of the parameter estimate.
+ To compute one standard deviation errors use ``perr = np.sqrt(np.diag(pcov))``.
+ If constraints are applied, this array does not contain the estimated covariance of
+ the parameters actually used during the fitting process but the uncertainties after
+ recalculating the covariance if all the parameters were free.
+ To get the actual uncertainties following error propagation of the actually fitted
+ parameters one should set full_output to True and access the uncertainties key.
+ ``infodict``: dict
+ a dictionary of optional outputs with the keys:
+
+ ``uncertainties``
+ The actual uncertainty on the optimized parameters.
+ ``nfev``
+ The number of function calls
+ ``fvec``
+ The function evaluated at the output
+ ``niter``
+ The number of iterations performed
+ ``chisq``
+ The chi square ``np.sum( ((f(xdata, *popt) - ydata) / sigma)**2 )``
+ ``reduced_chisq``
+ The chi square ``np.sum( ((f(xdata, *popt) - ydata) / sigma)**2 )`` divided
+ by the number of degrees of freedom ``(M - number_of_free_parameters)``
+ """
+ function_call_counter = 0
+ if numpy.isscalar(p0):
+ p0 = [p0]
+ parameters = numpy.array(p0, dtype=numpy.float64, copy=False)
+ if deltachi is None:
+ deltachi = 0.001
+
+ # NaNs can not be handled
+ if check_finite:
+ xdata = numpy.asarray_chkfinite(xdata)
+ ydata = numpy.asarray_chkfinite(ydata)
+ if sigma is not None:
+ sigma = numpy.asarray_chkfinite(sigma)
+ else:
+ sigma = numpy.ones((ydata.shape), dtype=numpy.float64)
+ ydata.shape = -1
+ sigma.shape = -1
+ else:
+ ydata = numpy.asarray(ydata)
+ xdata = numpy.asarray(xdata)
+ ydata.shape = -1
+ if sigma is not None:
+ sigma = numpy.asarray(sigma)
+ else:
+ sigma = numpy.ones((ydata.shape), dtype=numpy.float64)
+ sigma.shape = -1
+ # get rid of NaN in input data
+ idx = numpy.isfinite(ydata)
+ if False in idx:
+ # xdata must have a shape able to be understood by the user function
+ # in principle, one should not need to change it, however, if there are
+ # points to be excluded, one has to be able to exclude them.
+ # We can only hope that the sequence is properly arranged
+ if xdata.size == ydata.size:
+ if len(xdata.shape) != 1:
+ msg = "Need to reshape input xdata."
+ _logger.warning(msg)
+ xdata.shape = -1
+ else:
+ raise ValueError("Cannot reshape xdata to deal with NaN in ydata")
+ ydata = ydata[idx]
+ xdata = xdata[idx]
+ sigma = sigma[idx]
+ idx = numpy.isfinite(sigma)
+ if False in idx:
+ # xdata must have a shape able to be understood by the user function
+ # in principle, one should not need to change it, however, if there are
+ # points to be excluded, one has to be able to exclude them.
+ # We can only hope that the sequence is properly arranged
+ ydata = ydata[idx]
+ xdata = xdata[idx]
+ sigma = sigma[idx]
+ idx = numpy.isfinite(xdata)
+ filter_xdata = False
+ if False in idx:
+ # What to do?
+ try:
+ # Let's see if the function is able to deal with non-finite data
+ msg = "Checking if function can deal with non-finite data"
+ _logger.debug(msg)
+ evaluation = model(xdata, *parameters)
+ function_call_counter += 1
+ if evaluation.shape != ydata.shape:
+ if evaluation.size == ydata.size:
+ msg = "Supplied function does not return a proper array of floats."
+ msg += "\nFunction should be rewritten to return a 1D array of floats."
+ msg += "\nTrying to reshape output."
+ _logger.warning(msg)
+ evaluation.shape = ydata.shape
+ if False in numpy.isfinite(evaluation):
+ msg = "Supplied function unable to handle non-finite x data"
+ msg += "\nAttempting to filter out those x data values."
+ _logger.warning(msg)
+ filter_xdata = True
+ else:
+ filter_xdata = False
+ evaluation = None
+ except:
+ # function cannot handle input data
+ filter_xdata = True
+ if filter_xdata:
+ if xdata.size != ydata.size:
+ raise ValueError("xdata contains non-finite data that cannot be filtered")
+ else:
+ # we leave the xdata as they where
+ old_shape = xdata.shape
+ xdata.shape = ydata.shape
+ idx0 = numpy.isfinite(xdata)
+ xdata.shape = old_shape
+ ydata = ydata[idx0]
+ xdata = xdata[idx]
+ sigma = sigma[idx0]
+ weight = 1.0 / (sigma + numpy.equal(sigma, 0))
+ weight0 = weight * weight
+
+ nparameters = len(parameters)
+
+ if epsfcn is None:
+ epsfcn = numpy.finfo(numpy.float64).eps
+ else:
+ epsfcn = max(epsfcn, numpy.finfo(numpy.float64).eps)
+
+ # check if constraints have been passed as text
+ constrained_fit = False
+ if constraints is not None:
+ # make sure we work with a list of lists
+ input_constraints = constraints
+ tmp_constraints = [None] * len(input_constraints)
+ for i in range(nparameters):
+ tmp_constraints[i] = list(input_constraints[i])
+ constraints = tmp_constraints
+ for i in range(nparameters):
+ if hasattr(constraints[i][0], "upper"):
+ txt = constraints[i][0].upper()
+ if txt == "FREE":
+ constraints[i][0] = CFREE
+ elif txt == "POSITIVE":
+ constraints[i][0] = CPOSITIVE
+ elif txt == "QUOTED":
+ constraints[i][0] = CQUOTED
+ elif txt == "FIXED":
+ constraints[i][0] = CFIXED
+ elif txt == "FACTOR":
+ constraints[i][0] = CFACTOR
+ constraints[i][1] = int(constraints[i][1])
+ elif txt == "DELTA":
+ constraints[i][0] = CDELTA
+ constraints[i][1] = int(constraints[i][1])
+ elif txt == "SUM":
+ constraints[i][0] = CSUM
+ constraints[i][1] = int(constraints[i][1])
+ elif txt in ["IGNORED", "IGNORE"]:
+ constraints[i][0] = CIGNORED
+ else:
+ #I should raise an exception
+ raise ValueError("Unknown constraint %s" % constraints[i][0])
+ if constraints[i][0] > 0:
+ constrained_fit = True
+ if constrained_fit:
+ if full_output is None:
+ _logger.info("Recommended to set full_output to True when using constraints")
+
+ # Levenberg-Marquardt algorithm
+ fittedpar = parameters.__copy__()
+ flambda = 0.001
+ iiter = max_iter
+ #niter = 0
+ last_evaluation=None
+ x = xdata
+ y = ydata
+ chisq0 = -1
+ iteration_counter = 0
+ while (iiter > 0):
+ weight = weight0
+ """
+ I cannot evaluate the initial chisq here because I do not know
+ if some parameters are to be ignored, otherways I could do it as follows:
+ if last_evaluation is None:
+ yfit = model(x, *fittedpar)
+ last_evaluation = yfit
+ chisq0 = (weight * pow(y-yfit, 2)).sum()
+ and chisq would not need to be recalculated.
+ Passing the last_evaluation assumes that there are no parameters being
+ ignored or not between calls.
+ """
+ iteration_counter += 1
+ chisq0, alpha0, beta, internal_output = chisq_alpha_beta(
+ model, fittedpar,
+ x, y, weight, constraints=constraints,
+ model_deriv=model_deriv,
+ epsfcn=epsfcn,
+ left_derivative=left_derivative,
+ last_evaluation=last_evaluation,
+ full_output=True)
+ n_free = internal_output["n_free"]
+ free_index = internal_output["free_index"]
+ noigno = internal_output["noigno"]
+ fitparam = internal_output["fitparam"]
+ function_calls = internal_output["function_calls"]
+ function_call_counter += function_calls
+ #print("chisq0 = ", chisq0, n_free, fittedpar)
+ #raise
+ nr, nc = alpha0.shape
+ flag = 0
+ #lastdeltachi = chisq0
+ while flag == 0:
+ alpha = alpha0 * (1.0 + flambda * numpy.identity(nr))
+ deltapar = numpy.dot(beta, inv(alpha))
+ if constraints is None:
+ newpar = fitparam + deltapar [0]
+ else:
+ newpar = parameters.__copy__()
+ pwork = numpy.zeros(deltapar.shape, numpy.float64)
+ for i in range(n_free):
+ if constraints is None:
+ pwork [0] [i] = fitparam [i] + deltapar [0] [i]
+ elif constraints [free_index[i]][0] == CFREE:
+ pwork [0] [i] = fitparam [i] + deltapar [0] [i]
+ elif constraints [free_index[i]][0] == CPOSITIVE:
+ #abs method
+ pwork [0] [i] = fitparam [i] + deltapar [0] [i]
+ #square method
+ #pwork [0] [i] = (numpy.sqrt(fitparam [i]) + deltapar [0] [i]) * \
+ # (numpy.sqrt(fitparam [i]) + deltapar [0] [i])
+ elif constraints[free_index[i]][0] == CQUOTED:
+ pmax = max(constraints[free_index[i]][1],
+ constraints[free_index[i]][2])
+ pmin = min(constraints[free_index[i]][1],
+ constraints[free_index[i]][2])
+ A = 0.5 * (pmax + pmin)
+ B = 0.5 * (pmax - pmin)
+ if B != 0:
+ pwork [0] [i] = A + \
+ B * numpy.sin(numpy.arcsin((fitparam[i] - A)/B)+ \
+ deltapar [0] [i])
+ else:
+ txt = "Error processing constrained fit\n"
+ txt += "Parameter limits are %g and %g\n" % (pmin, pmax)
+ txt += "A = %g B = %g" % (A, B)
+ raise ValueError("Invalid parameter limits")
+ newpar[free_index[i]] = pwork [0] [i]
+ newpar = numpy.array(_get_parameters(newpar, constraints))
+ workpar = numpy.take(newpar, noigno)
+ yfit = model(x, *workpar)
+ if last_evaluation is None:
+ if len(yfit.shape) > 1:
+ msg = "Supplied function does not return a 1D array of floats."
+ msg += "\nFunction should be rewritten."
+ msg += "\nTrying to reshape output."
+ _logger.warning(msg)
+ yfit.shape = -1
+ function_call_counter += 1
+ chisq = (weight * pow(y-yfit, 2)).sum()
+ absdeltachi = chisq0 - chisq
+ if absdeltachi < 0:
+ flambda *= 10.0
+ if flambda > 1000:
+ flag = 1
+ iiter = 0
+ else:
+ flag = 1
+ fittedpar = newpar.__copy__()
+ lastdeltachi = 100 * (absdeltachi / (chisq + (chisq == 0)))
+ if iteration_counter < 2:
+ # ignore any limit, the fit *has* to be improved
+ pass
+ elif (lastdeltachi) < deltachi:
+ iiter = 0
+ elif absdeltachi < numpy.sqrt(epsfcn):
+ iiter = 0
+ _logger.info("Iteration finished due to too small absolute chi decrement")
+ chisq0 = chisq
+ flambda = flambda / 10.0
+ last_evaluation = yfit
+ iiter = iiter - 1
+ # this is the covariance matrix of the actually fitted parameters
+ cov0 = inv(alpha0)
+ if constraints is None:
+ cov = cov0
+ else:
+ # yet another call needed with all the parameters being free except those
+ # that are FIXED and that will be assigned a 100 % uncertainty.
+ new_constraints = copy.deepcopy(constraints)
+ flag_special = [0] * len(fittedpar)
+ for idx, constraint in enumerate(constraints):
+ if constraints[idx][0] in [CFIXED, CIGNORED]:
+ flag_special[idx] = constraints[idx][0]
+ else:
+ new_constraints[idx][0] = CFREE
+ new_constraints[idx][1] = 0
+ new_constraints[idx][2] = 0
+ chisq, alpha, beta, internal_output = chisq_alpha_beta(
+ model, fittedpar,
+ x, y, weight, constraints=new_constraints,
+ model_deriv=model_deriv,
+ epsfcn=epsfcn,
+ left_derivative=left_derivative,
+ last_evaluation=last_evaluation,
+ full_output=True)
+ # obtained chisq should be identical to chisq0
+ try:
+ cov = inv(alpha)
+ except LinAlgError:
+ _logger.critical("Error calculating covariance matrix after successful fit")
+ cov = None
+ if cov is not None:
+ for idx, value in enumerate(flag_special):
+ if value in [CFIXED, CIGNORED]:
+ cov = numpy.insert(numpy.insert(cov, idx, 0, axis=1), idx, 0, axis=0)
+ cov[idx, idx] = fittedpar[idx] * fittedpar[idx]
+
+ if not full_output:
+ return fittedpar, cov
+ else:
+ sigma0 = numpy.sqrt(abs(numpy.diag(cov0)))
+ sigmapar = _get_sigma_parameters(fittedpar, sigma0, constraints)
+ ddict = {}
+ ddict["chisq"] = chisq0
+ ddict["reduced_chisq"] = chisq0 / (len(yfit)-n_free)
+ ddict["covariance"] = cov0
+ ddict["uncertainties"] = sigmapar
+ ddict["fvec"] = last_evaluation
+ ddict["nfev"] = function_call_counter
+ ddict["niter"] = iteration_counter
+ return fittedpar, cov, ddict #, chisq/(len(yfit)-len(sigma0)), sigmapar,niter,lastdeltachi
+
+def chisq_alpha_beta(model, parameters, x, y, weight, constraints=None,
+ model_deriv=None, epsfcn=None, left_derivative=False,
+ last_evaluation=None, full_output=False):
+
+ """
+ Get chi square, the curvature matrix alpha and the matrix beta according to the input parameters.
+ If all the parameters are unconstrained, the covariance matrix is the inverse of the alpha matrix.
+
+ :param model: callable
+ The model function, f(x, ...). It must take the independent
+ variable as the first argument and the parameters to fit as
+ separate remaining arguments.
+ The returned value is a one dimensional array of floats.
+
+ :param parameters: N-length sequence
+ Values of parameters at which function and derivatives are to be calculated.
+
+ :param x: An M-length sequence.
+ The independent variable where the data is measured.
+
+ :param y: An M-length sequence
+ The dependent data --- nominally f(xdata, ...)
+
+ :param weight: M-length sequence
+ Weights to be applied in the calculation of chi square
+ As a reminder ``chisq = np.sum(weigth * (model(x, *parameters) - y)**2)``
+
+ :param constraints:
+ If provided, it is a 2D sequence of dimension (n_parameters, 3) where,
+ for each parameter denoted by the index i, the meaning is
+
+ - constraints[i][0]
+
+ - 0 - Free (CFREE)
+ - 1 - Positive (CPOSITIVE)
+ - 2 - Quoted (CQUOTED)
+ - 3 - Fixed (CFIXED)
+ - 4 - Factor (CFACTOR)
+ - 5 - Delta (CDELTA)
+ - 6 - Sum (CSUM)
+
+
+ - constraints[i][1]
+
+ - Ignored if constraints[i][0] is 0, 1, 3
+ - Min value of the parameter if constraints[i][0] is CQUOTED
+ - Index of fitted parameter to which it is related
+
+ - constraints[i][2]
+
+ - Ignored if constraints[i][0] is 0, 1, 3
+ - Max value of the parameter if constraints[i][0] is CQUOTED
+ - Factor to apply to related parameter with index constraints[i][1]
+ - Difference with parameter with index constraints[i][1]
+ - Sum obtained when adding parameter with index constraints[i][1]
+ :type constraints: *optional*, None or 2D sequence
+
+ :param model_deriv:
+ None (default) or function providing the derivatives of the fitting function respect to the fitted parameters.
+ It will be called as model_deriv(xdata, parameters, index) where parameters is a sequence with the current
+ values of the fitting parameters, index is the fitting parameter index for which the the derivative has
+ to be provided in the supplied array of xdata points.
+ :type model_deriv: *optional*, None or callable
+
+
+ :param epsfcn: float
+ A variable used in determining a suitable parameter variation when
+ calculating the numerical derivatives (for model_deriv=None).
+ Normally the actual step length will be sqrt(epsfcn)*x
+ Original Gefit module was using epsfcn 1.0e-10 while default value
+ is now numpy.finfo(numpy.float64).eps as in scipy
+ :type epsfcn: *optional*, float
+
+ :param left_derivative:
+ This parameter only has an influence if no derivative function
+ is provided. When True the left and right derivatives of the
+ model will be calculated for each fitted parameters thus leading to
+ the double number of function evaluations. Default is False.
+ Original Gefit module was always using left_derivative as True.
+ :type left_derivative: *optional*, bool
+
+ :param last_evaluation: An M-length array
+ Used for optimization purposes. If supplied, this array will be taken as the result of
+ evaluating the function, that is as the result of ``model(x, *parameters)`` thus avoiding
+ the evaluation call.
+
+ :param full_output: bool, optional
+ Additional output used for internal purposes with the keys:
+ ``function_calls``
+ The number of model function calls performed.
+ ``fitparam``
+ A sequence with the actual free parameters
+ ``free_index``
+ Sequence with the indices of the free parameters in input parameters sequence.
+ ``noigno``
+ Sequence with the indices of the original parameters considered in the calculations.
+ """
+ if epsfcn is None:
+ epsfcn = numpy.finfo(numpy.float64).eps
+ else:
+ epsfcn = max(epsfcn, numpy.finfo(numpy.float64).eps)
+ #nr0, nc = data.shape
+ n_param = len(parameters)
+ if constraints is None:
+ derivfactor = numpy.ones((n_param, ))
+ n_free = n_param
+ noigno = numpy.arange(n_param)
+ free_index = noigno * 1
+ fitparam = parameters * 1
+ else:
+ n_free = 0
+ fitparam = []
+ free_index = []
+ noigno = []
+ derivfactor = []
+ for i in range(n_param):
+ if constraints[i][0] != CIGNORED:
+ noigno.append(i)
+ if constraints[i][0] == CFREE:
+ fitparam.append(parameters [i])
+ derivfactor.append(1.0)
+ free_index.append(i)
+ n_free += 1
+ elif constraints[i][0] == CPOSITIVE:
+ fitparam.append(abs(parameters[i]))
+ derivfactor.append(1.0)
+ #fitparam.append(numpy.sqrt(abs(parameters[i])))
+ #derivfactor.append(2.0*numpy.sqrt(abs(parameters[i])))
+ free_index.append(i)
+ n_free += 1
+ elif constraints[i][0] == CQUOTED:
+ pmax = max(constraints[i][1], constraints[i][2])
+ pmin =min(constraints[i][1], constraints[i][2])
+ if ((pmax-pmin) > 0) & \
+ (parameters[i] <= pmax) & \
+ (parameters[i] >= pmin):
+ A = 0.5 * (pmax + pmin)
+ B = 0.5 * (pmax - pmin)
+ fitparam.append(parameters[i])
+ derivfactor.append(B*numpy.cos(numpy.arcsin((parameters[i] - A)/B)))
+ free_index.append(i)
+ n_free += 1
+ elif (pmax-pmin) > 0:
+ print("WARNING: Quoted parameter outside boundaries")
+ print("Initial value = %f" % parameters[i])
+ print("Limits are %f and %f" % (pmin, pmax))
+ print("Parameter will be kept at its starting value")
+ fitparam = numpy.array(fitparam, numpy.float64)
+ alpha = numpy.zeros((n_free, n_free), numpy.float64)
+ beta = numpy.zeros((1, n_free), numpy.float64)
+ #delta = (fitparam + numpy.equal(fitparam, 0.0)) * 0.00001
+ delta = (fitparam + numpy.equal(fitparam, 0.0)) * numpy.sqrt(epsfcn)
+ nr = y.size
+ ##############
+ # Prior to each call to the function one has to re-calculate the
+ # parameters
+ pwork = parameters.__copy__()
+ for i in range(n_free):
+ pwork [free_index[i]] = fitparam [i]
+ if n_free == 0:
+ raise ValueError("No free parameters to fit")
+ function_calls = 0
+ if not left_derivative:
+ if last_evaluation is not None:
+ f2 = last_evaluation
+ else:
+ f2 = model(x, *parameters)
+ f2.shape = -1
+ function_calls += 1
+ for i in range(n_free):
+ if model_deriv is None:
+ #pwork = parameters.__copy__()
+ pwork[free_index[i]] = fitparam [i] + delta [i]
+ newpar = _get_parameters(pwork.tolist(), constraints)
+ newpar = numpy.take(newpar, noigno)
+ f1 = model(x, *newpar)
+ f1.shape = -1
+ function_calls += 1
+ if left_derivative:
+ pwork[free_index[i]] = fitparam [i] - delta [i]
+ newpar = _get_parameters(pwork.tolist(), constraints)
+ newpar=numpy.take(newpar, noigno)
+ f2 = model(x, *newpar)
+ function_calls += 1
+ help0 = (f1 - f2) / (2.0 * delta[i])
+ else:
+ help0 = (f1 - f2) / (delta[i])
+ help0 = help0 * derivfactor[i]
+ pwork[free_index[i]] = fitparam [i]
+ #removed I resize outside the loop:
+ #help0 = numpy.resize(help0, (1, nr))
+ else:
+ help0 = model_deriv(x, pwork, free_index[i])
+ help0 = help0 * derivfactor[i]
+
+ if i == 0:
+ deriv = help0
+ else:
+ deriv = numpy.concatenate((deriv, help0), 0)
+
+ #line added to resize outside the loop
+ deriv = numpy.resize(deriv, (n_free, nr))
+ if last_evaluation is None:
+ if constraints is None:
+ yfit = model(x, *fitparam)
+ yfit.shape = -1
+ else:
+ newpar = _get_parameters(pwork.tolist(), constraints)
+ newpar = numpy.take(newpar, noigno)
+ yfit = model(x, *newpar)
+ yfit.shape = -1
+ function_calls += 1
+ else:
+ yfit = last_evaluation
+ deltay = y - yfit
+ help0 = weight * deltay
+ for i in range(n_free):
+ derivi = numpy.resize(deriv[i, :], (1, nr))
+ help1 = numpy.resize(numpy.sum((help0 * derivi), 1), (1, 1))
+ if i == 0:
+ beta = help1
+ else:
+ beta = numpy.concatenate((beta, help1), 1)
+ help1 = numpy.inner(deriv, weight*derivi)
+ if i == 0:
+ alpha = help1
+ else:
+ alpha = numpy.concatenate((alpha, help1), 1)
+ chisq = (help0 * deltay).sum()
+ if full_output:
+ ddict = {}
+ ddict["n_free"] = n_free
+ ddict["free_index"] = free_index
+ ddict["noigno"] = noigno
+ ddict["fitparam"] = fitparam
+ ddict["derivfactor"] = derivfactor
+ ddict["function_calls"] = function_calls
+ return chisq, alpha, beta, ddict
+ else:
+ return chisq, alpha, beta
+
+
+def _get_parameters(parameters, constraints):
+ """
+ Apply constraints to input parameters.
+
+ Parameters not depending on other parameters, they are returned as the input.
+
+ Parameters depending on other parameters, return the value after applying the
+ relation to the parameter wo which they are related.
+ """
+ # 0 = Free 1 = Positive 2 = Quoted
+ # 3 = Fixed 4 = Factor 5 = Delta
+ if constraints is None:
+ return parameters * 1
+ newparam = []
+ #first I make the free parameters
+ #because the quoted ones put troubles
+ for i in range(len(constraints)):
+ if constraints[i][0] == CFREE:
+ newparam.append(parameters[i])
+ elif constraints[i][0] == CPOSITIVE:
+ #newparam.append(parameters[i] * parameters[i])
+ newparam.append(abs(parameters[i]))
+ elif constraints[i][0] == CQUOTED:
+ newparam.append(parameters[i])
+ elif abs(constraints[i][0]) == CFIXED:
+ newparam.append(parameters[i])
+ else:
+ newparam.append(parameters[i])
+ for i in range(len(constraints)):
+ if constraints[i][0] == CFACTOR:
+ newparam[i] = constraints[i][2] * newparam[int(constraints[i][1])]
+ elif constraints[i][0] == CDELTA:
+ newparam[i] = constraints[i][2] + newparam[int(constraints[i][1])]
+ elif constraints[i][0] == CIGNORED:
+ # The whole ignored stuff should not be documented because setting
+ # a parameter to 0 is not the same as being ignored.
+ # Being ignored should imply the parameter is simply not accounted for
+ # and should be stripped out of the list of parameters by the program
+ # using this module
+ newparam[i] = 0
+ elif constraints[i][0] == CSUM:
+ newparam[i] = constraints[i][2]-newparam[int(constraints[i][1])]
+ return newparam
+
+
+def _get_sigma_parameters(parameters, sigma0, constraints):
+ """
+ Internal function propagating the uncertainty on the actually fitted parameters and related parameters to the
+ final parameters considering the applied constraints.
+
+ Parameters
+ ----------
+ parameters : 1D sequence of length equal to the number of free parameters N
+ The parameters actually used in the fitting process.
+ sigma0 : 1D sequence of length N
+ Uncertainties calculated as the square-root of the diagonal of
+ the covariance matrix
+ constraints : The set of constraints applied in the fitting process
+ """
+ # 0 = Free 1 = Positive 2 = Quoted
+ # 3 = Fixed 4 = Factor 5 = Delta
+ if constraints is None:
+ return sigma0
+ n_free = 0
+ sigma_par = numpy.zeros(parameters.shape, numpy.float64)
+ for i in range(len(constraints)):
+ if constraints[i][0] == CFREE:
+ sigma_par [i] = sigma0[n_free]
+ n_free += 1
+ elif constraints[i][0] == CPOSITIVE:
+ #sigma_par [i] = 2.0 * sigma0[n_free]
+ sigma_par [i] = sigma0[n_free]
+ n_free += 1
+ elif constraints[i][0] == CQUOTED:
+ pmax = max(constraints [i][1], constraints [i][2])
+ pmin = min(constraints [i][1], constraints [i][2])
+ # A = 0.5 * (pmax + pmin)
+ B = 0.5 * (pmax - pmin)
+ if (B > 0) & (parameters [i] < pmax) & (parameters [i] > pmin):
+ sigma_par [i] = abs(B * numpy.cos(parameters[i]) * sigma0[n_free])
+ n_free += 1
+ else:
+ sigma_par [i] = parameters[i]
+ elif abs(constraints[i][0]) == CFIXED:
+ sigma_par[i] = parameters[i]
+ for i in range(len(constraints)):
+ if constraints[i][0] == CFACTOR:
+ sigma_par [i] = constraints[i][2]*sigma_par[int(constraints[i][1])]
+ elif constraints[i][0] == CDELTA:
+ sigma_par [i] = sigma_par[int(constraints[i][1])]
+ elif constraints[i][0] == CSUM:
+ sigma_par [i] = sigma_par[int(constraints[i][1])]
+ return sigma_par
+
+
+def main(argv=None):
+ if argv is None:
+ npoints = 10000
+ elif hasattr(argv, "__len__"):
+ if len(argv) > 1:
+ npoints = int(argv[1])
+ else:
+ print("Usage:")
+ print("fit [npoints]")
+ else:
+ # expected a number
+ npoints = argv
+
+ def gauss(t0, *param0):
+ param = numpy.array(param0)
+ t = numpy.array(t0)
+ dummy = 2.3548200450309493 * (t - param[3]) / param[4]
+ return param[0] + param[1] * t + param[2] * myexp(-0.5 * dummy * dummy)
+
+
+ def myexp(x):
+ # put a (bad) filter to avoid over/underflows
+ # with no python looping
+ return numpy.exp(x * numpy.less(abs(x), 250)) -\
+ 1.0 * numpy.greater_equal(abs(x), 250)
+
+ xx = numpy.arange(npoints, dtype=numpy.float64)
+ yy = gauss(xx, *[10.5, 2, 1000.0, 20., 15])
+ sy = numpy.sqrt(abs(yy))
+ parameters = [0.0, 1.0, 900.0, 25., 10]
+ stime = time.time()
+
+ fittedpar, cov, ddict = leastsq(gauss, xx, yy, parameters,
+ sigma=sy,
+ left_derivative=False,
+ full_output=True,
+ check_finite=True)
+ etime = time.time()
+ sigmapars = numpy.sqrt(numpy.diag(cov))
+ print("Took ", etime - stime, "seconds")
+ print("Function calls = ", ddict["nfev"])
+ print("chi square = ", ddict["chisq"])
+ print("Fitted pars = ", fittedpar)
+ print("Sigma pars = ", sigmapars)
+ try:
+ from scipy.optimize import curve_fit as cfit
+ SCIPY = True
+ except ImportError:
+ SCIPY = False
+ if SCIPY:
+ counter = 0
+ stime = time.time()
+ scipy_fittedpar, scipy_cov = cfit(gauss,
+ xx,
+ yy,
+ parameters,
+ sigma=sy)
+ etime = time.time()
+ print("Scipy Took ", etime - stime, "seconds")
+ print("Counter = ", counter)
+ print("scipy = ", scipy_fittedpar)
+ print("Sigma = ", numpy.sqrt(numpy.diag(scipy_cov)))
+
+if __name__ == "__main__":
+ main()