diff options
Diffstat (limited to 'src/silx/math/fit/fittheories.py')
| -rw-r--r-- | src/silx/math/fit/fittheories.py | 1580 |
1 files changed, 1580 insertions, 0 deletions
diff --git a/src/silx/math/fit/fittheories.py b/src/silx/math/fit/fittheories.py new file mode 100644 index 0000000..f20cbf1 --- /dev/null +++ b/src/silx/math/fit/fittheories.py @@ -0,0 +1,1580 @@ +# /*########################################################################## +# +# Copyright (c) 2004-2023 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 modules provides a set of fit functions and associated +estimation functions in a format that can be imported into a +:class:`silx.math.fit.FitManager` instance. + +These functions are well suited for fitting multiple gaussian shaped peaks +typically found in spectroscopy data. The estimation functions are designed +to detect how many peaks are present in the data, and provide an initial +estimate for their height, their center location and their full-width +at half maximum (fwhm). + +The limitation of these estimation algorithms is that only gaussians having a +similar fwhm can be detected by the peak search algorithm. +This *search fwhm* can be defined by the user, if +he knows the characteristics of his data, or can be automatically estimated +based on the fwhm of the largest peak in the data. + +The source code of this module can serve as template for defining your own +fit functions. + +The functions to be imported by :meth:`FitManager.loadtheories` are defined by +a dictionary :const:`THEORY`: with the following structure:: + + from silx.math.fit.fittheory import FitTheory + + THEORY = { + 'theory_name_1': FitTheory( + description='Description of theory 1', + function=fitfunction1, + parameters=('param name 1', 'param name 2', …), + estimate=estimation_function1, + configure=configuration_function1, + derivative=derivative_function1), + + 'theory_name_2': FitTheory(…), + } + +.. note:: + + The order of the provided dictionary is taken into account. + +Theory names can be customized (e.g. ``gauss, lorentz, splitgauss``…). + +The mandatory parameters for :class:`FitTheory` are ``function`` and +``parameters``. + +You can also define an ``INIT`` function that will be executed by +:meth:`FitManager.loadtheories`. + +See the documentation of :class:`silx.math.fit.fittheory.FitTheory` +for more information. + +Module members: +--------------- +""" +import numpy +import logging + +from silx.math.fit import functions +from silx.math.fit.peaks import peak_search, guess_fwhm +from silx.math.fit.filters import strip, savitsky_golay +from silx.math.fit.leastsq import leastsq +from silx.math.fit.fittheory import FitTheory + +_logger = logging.getLogger(__name__) + +__authors__ = ["V.A. Sole", "P. Knobel"] +__license__ = "MIT" +__date__ = "15/05/2017" + + +DEFAULT_CONFIG = { + "NoConstraintsFlag": False, + "PositiveFwhmFlag": True, + "PositiveHeightAreaFlag": True, + "SameFwhmFlag": False, + "QuotedPositionFlag": False, # peak not outside data range + "QuotedEtaFlag": False, # force 0 < eta < 1 + # Peak detection + "AutoScaling": False, + "Yscaling": 1.0, + "FwhmPoints": 8, + "AutoFwhm": True, + "Sensitivity": 2.5, + "ForcePeakPresence": True, + # Hypermet + "HypermetTails": 15, + "QuotedFwhmFlag": 0, + "MaxFwhm2InputRatio": 1.5, + "MinFwhm2InputRatio": 0.4, + # short tail parameters + "MinGaussArea4ShortTail": 50000.0, + "InitialShortTailAreaRatio": 0.050, + "MaxShortTailAreaRatio": 0.100, + "MinShortTailAreaRatio": 0.0010, + "InitialShortTailSlopeRatio": 0.70, + "MaxShortTailSlopeRatio": 2.00, + "MinShortTailSlopeRatio": 0.50, + # long tail parameters + "MinGaussArea4LongTail": 1000.0, + "InitialLongTailAreaRatio": 0.050, + "MaxLongTailAreaRatio": 0.300, + "MinLongTailAreaRatio": 0.010, + "InitialLongTailSlopeRatio": 20.0, + "MaxLongTailSlopeRatio": 50.0, + "MinLongTailSlopeRatio": 5.0, + # step tail + "MinGaussHeight4StepTail": 5000.0, + "InitialStepTailHeightRatio": 0.002, + "MaxStepTailHeightRatio": 0.0100, + "MinStepTailHeightRatio": 0.0001, + # Hypermet constraints + # position in range [estimated position +- estimated fwhm/2] + "HypermetQuotedPositionFlag": True, + "DeltaPositionFwhmUnits": 0.5, + "SameSlopeRatioFlag": 1, + "SameAreaRatioFlag": 1, + # Strip bg removal + "StripBackgroundFlag": True, + "SmoothingFlag": True, + "SmoothingWidth": 5, + "StripWidth": 2, + "StripIterations": 5000, + "StripThresholdFactor": 1.0, +} +"""This dictionary defines default configuration parameters that have effects +on fit functions and estimation functions, mainly on fit constraints. +This dictionary is accessible as attribute :attr:`FitTheories.config`, +which can be modified by configuration functions defined in +:const:`CONFIGURE`. +""" + +CFREE = 0 +CPOSITIVE = 1 +CQUOTED = 2 +CFIXED = 3 +CFACTOR = 4 +CDELTA = 5 +CSUM = 6 +CIGNORED = 7 + + +class FitTheories(object): + """Class wrapping functions from :class:`silx.math.fit.functions` + and providing estimate functions for all of these fit functions.""" + + def __init__(self, config=None): + if config is None: + self.config = DEFAULT_CONFIG + else: + self.config = config + + def ahypermet(self, x, *pars): + """ + Wrapping of :func:`silx.math.fit.functions.sum_ahypermet` without + the tail flags in the function signature. + + Depending on the value of `self.config['HypermetTails']`, one can + activate or deactivate the various terms of the hypermet function. + + `self.config['HypermetTails']` must be an integer between 0 and 15. + It is a set of 4 binary flags, one for activating each one of the + hypermet terms: *gaussian function, short tail, long tail, step*. + + For example, 15 can be expressed as ``1111`` in base 2, so a flag of + 15 means all terms are active. + """ + g_term = self.config["HypermetTails"] & 1 + st_term = (self.config["HypermetTails"] >> 1) & 1 + lt_term = (self.config["HypermetTails"] >> 2) & 1 + step_term = (self.config["HypermetTails"] >> 3) & 1 + return functions.sum_ahypermet( + x, + *pars, + gaussian_term=g_term, + st_term=st_term, + lt_term=lt_term, + step_term=step_term, + ) + + def poly(self, x, *pars): + """Order n polynomial. + The order of the polynomial is defined by the number of + coefficients (``*pars``). + + """ + p = numpy.poly1d(pars) + return p(x) + + @staticmethod + def estimate_poly(x, y, n=2): + """Estimate polynomial coefficients for a degree n polynomial.""" + pcoeffs = numpy.polyfit(x, y, n) + constraints = numpy.zeros((n + 1, 3), numpy.float64) + return pcoeffs, constraints + + def estimate_quadratic(self, x, y): + """Estimate quadratic coefficients""" + return self.estimate_poly(x, y, n=2) + + def estimate_cubic(self, x, y): + """Estimate coefficients for a degree 3 polynomial""" + return self.estimate_poly(x, y, n=3) + + def estimate_quartic(self, x, y): + """Estimate coefficients for a degree 4 polynomial""" + return self.estimate_poly(x, y, n=4) + + def estimate_quintic(self, x, y): + """Estimate coefficients for a degree 5 polynomial""" + return self.estimate_poly(x, y, n=5) + + def strip_bg(self, y): + """Return the strip background of y, using parameters from + :attr:`config` dictionary (*StripBackgroundFlag, StripWidth, + StripIterations, StripThresholdFactor*)""" + remove_strip_bg = self.config.get("StripBackgroundFlag", False) + if remove_strip_bg: + if self.config["SmoothingFlag"]: + y = savitsky_golay(y, self.config["SmoothingWidth"]) + strip_width = self.config["StripWidth"] + strip_niterations = self.config["StripIterations"] + strip_thr_factor = self.config["StripThresholdFactor"] + return strip( + y, w=strip_width, niterations=strip_niterations, factor=strip_thr_factor + ) + else: + return numpy.zeros_like(y) + + def guess_yscaling(self, y): + """Estimate scaling for y prior to peak search. + A smoothing filter is applied to y to estimate the noise level + (chi-squared) + + :param y: Data array + :return: Scaling factor + """ + # ensure y is an array + yy = numpy.array(y, copy=False) + + # smooth + convolution_kernel = numpy.ones(shape=(3,)) / 3.0 + ysmooth = numpy.convolve(y, convolution_kernel, mode="same") + + # remove zeros + idx_array = numpy.fabs(y) > 0.0 + yy = yy[idx_array] + ysmooth = ysmooth[idx_array] + + # compute scaling factor + chisq = numpy.mean((yy - ysmooth) ** 2 / numpy.fabs(yy)) + if chisq > 0: + return 1.0 / chisq + else: + return 1.0 + + def peak_search(self, y, fwhm, sensitivity): + """Search for peaks in y array, after padding the array and + multiplying its value by a scaling factor. + + :param y: 1-D data array + :param int fwhm: Typical full width at half maximum for peaks, + in number of points. This parameter is used for to discriminate between + true peaks and background fluctuations. + :param float sensitivity: Sensitivity parameter. This is a threshold factor + for peak detection. Only peaks larger than the standard deviation + of the noise multiplied by this sensitivity parameter are detected. + :return: List of peak indices + """ + # add padding + ysearch = numpy.ones((len(y) + 2 * fwhm,), numpy.float64) + ysearch[0:fwhm] = y[0] + ysearch[-1 : -fwhm - 1 : -1] = y[len(y) - 1] + ysearch[fwhm : fwhm + len(y)] = y[:] + + scaling = ( + self.guess_yscaling(y) + if self.config["AutoScaling"] + else self.config["Yscaling"] + ) + + if len(ysearch) > 1.5 * fwhm: + peaks = peak_search(scaling * ysearch, fwhm=fwhm, sensitivity=sensitivity) + return [ + peak_index - fwhm + for peak_index in peaks + if 0 <= peak_index - fwhm < len(y) + ] + else: + return [] + + def estimate_height_position_fwhm(self, x, y): + """Estimation of *Height, Position, FWHM* of peaks, for gaussian-like + curves. + + This functions finds how many parameters are needed, based on the + number of peaks detected. Then it estimates the fit parameters + with a few iterations of fitting gaussian functions. + + :param x: Array of abscissa values + :param y: Array of ordinate values (``y = f(x)``) + :return: Tuple of estimated fit parameters and fit constraints. + Parameters to be estimated for each peak are: + *Height, Position, FWHM*. + Fit constraints depend on :attr:`config`. + """ + fittedpar = [] + + bg = self.strip_bg(y) + + if self.config["AutoFwhm"]: + search_fwhm = guess_fwhm(y) + else: + search_fwhm = int(float(self.config["FwhmPoints"])) + search_sens = float(self.config["Sensitivity"]) + + if search_fwhm < 3: + _logger.warning("Setting peak fwhm to 3 (lower limit)") + search_fwhm = 3 + self.config["FwhmPoints"] = 3 + + if search_sens < 1: + _logger.warning( + "Setting peak search sensitivity to 1. " + + "(lower limit to filter out noise peaks)" + ) + search_sens = 1 + self.config["Sensitivity"] = 1 + + npoints = len(y) + + # Find indices of peaks in data array + peaks = self.peak_search(y, fwhm=search_fwhm, sensitivity=search_sens) + + if not len(peaks): + forcepeak = int(float(self.config.get("ForcePeakPresence", 0))) + if forcepeak: + delta = y - bg + # get index of global maximum + # (first one if several samples are equal to this value) + peaks = [numpy.nonzero(delta == delta.max())[0][0]] + + # Find index of largest peak in peaks array + index_largest_peak = 0 + if len(peaks) > 0: + # estimate fwhm as 5 * sampling interval + sig = 5 * abs(x[npoints - 1] - x[0]) / npoints + peakpos = x[int(peaks[0])] + if abs(peakpos) < 1.0e-16: + peakpos = 0.0 + param = numpy.array([y[int(peaks[0])] - bg[int(peaks[0])], peakpos, sig]) + height_largest_peak = param[0] + peak_index = 1 + for i in peaks[1:]: + param2 = numpy.array([y[int(i)] - bg[int(i)], x[int(i)], sig]) + param = numpy.concatenate((param, param2)) + if param2[0] > height_largest_peak: + height_largest_peak = param2[0] + index_largest_peak = peak_index + peak_index += 1 + + # Subtract background + xw = x + yw = y - bg + + cons = numpy.zeros((len(param), 3), numpy.float64) + + # peak height must be positive + cons[0 : len(param) : 3, 0] = CPOSITIVE + # force peaks to stay around their position + cons[1 : len(param) : 3, 0] = CQUOTED + + # set possible peak range to estimated peak +- guessed fwhm + if len(xw) > search_fwhm: + fwhmx = numpy.fabs(xw[int(search_fwhm)] - xw[0]) + cons[1 : len(param) : 3, 1] = param[1 : len(param) : 3] - 0.5 * fwhmx + cons[1 : len(param) : 3, 2] = param[1 : len(param) : 3] + 0.5 * fwhmx + else: + shape = [max(1, int(x)) for x in (param[1 : len(param) : 3])] + cons[1 : len(param) : 3, 1] = min(xw) * numpy.ones(shape, numpy.float64) + cons[1 : len(param) : 3, 2] = max(xw) * numpy.ones(shape, numpy.float64) + + # ensure fwhm is positive + cons[2 : len(param) : 3, 0] = CPOSITIVE + + # run a quick iterative fit (4 iterations) to improve + # estimations + fittedpar, _, _ = leastsq( + functions.sum_gauss, + xw, + yw, + param, + max_iter=4, + constraints=cons.tolist(), + full_output=True, + ) + + # set final constraints based on config parameters + cons = numpy.zeros((len(fittedpar), 3), numpy.float64) + peak_index = 0 + for i in range(len(peaks)): + # Setup height area constrains + if not self.config["NoConstraintsFlag"]: + if self.config["PositiveHeightAreaFlag"]: + cons[peak_index, 0] = CPOSITIVE + cons[peak_index, 1] = 0 + cons[peak_index, 2] = 0 + peak_index += 1 + + # Setup position constrains + if not self.config["NoConstraintsFlag"]: + if self.config["QuotedPositionFlag"]: + cons[peak_index, 0] = CQUOTED + cons[peak_index, 1] = min(x) + cons[peak_index, 2] = max(x) + peak_index += 1 + + # Setup positive FWHM constrains + if not self.config["NoConstraintsFlag"]: + if self.config["PositiveFwhmFlag"]: + cons[peak_index, 0] = CPOSITIVE + cons[peak_index, 1] = 0 + cons[peak_index, 2] = 0 + if self.config["SameFwhmFlag"]: + if i != index_largest_peak: + cons[peak_index, 0] = CFACTOR + cons[peak_index, 1] = 3 * index_largest_peak + 2 + cons[peak_index, 2] = 1.0 + peak_index += 1 + + return fittedpar, cons + + def estimate_agauss(self, x, y): + """Estimation of *Area, Position, FWHM* of peaks, for gaussian-like + curves. + + This functions uses :meth:`estimate_height_position_fwhm`, then + converts the height parameters to area under the curve with the + formula ``area = sqrt(2*pi) * height * fwhm / (2 * sqrt(2 * log(2))`` + + :param x: Array of abscissa values + :param y: Array of ordinate values (``y = f(x)``) + :return: Tuple of estimated fit parameters and fit constraints. + Parameters to be estimated for each peak are: + *Area, Position, FWHM*. + Fit constraints depend on :attr:`config`. + """ + fittedpar, cons = self.estimate_height_position_fwhm(x, y) + # get the number of found peaks + npeaks = len(fittedpar) // 3 + for i in range(npeaks): + height = fittedpar[3 * i] + fwhm = fittedpar[3 * i + 2] + # Replace height with area in fittedpar + fittedpar[3 * i] = ( + numpy.sqrt(2 * numpy.pi) + * height + * fwhm + / (2.0 * numpy.sqrt(2 * numpy.log(2))) + ) + return fittedpar, cons + + def estimate_alorentz(self, x, y): + """Estimation of *Area, Position, FWHM* of peaks, for Lorentzian + curves. + + This functions uses :meth:`estimate_height_position_fwhm`, then + converts the height parameters to area under the curve with the + formula ``area = height * fwhm * 0.5 * pi`` + + :param x: Array of abscissa values + :param y: Array of ordinate values (``y = f(x)``) + :return: Tuple of estimated fit parameters and fit constraints. + Parameters to be estimated for each peak are: + *Area, Position, FWHM*. + Fit constraints depend on :attr:`config`. + """ + fittedpar, cons = self.estimate_height_position_fwhm(x, y) + # get the number of found peaks + npeaks = len(fittedpar) // 3 + for i in range(npeaks): + height = fittedpar[3 * i] + fwhm = fittedpar[3 * i + 2] + # Replace height with area in fittedpar + fittedpar[3 * i] = height * fwhm * 0.5 * numpy.pi + return fittedpar, cons + + def estimate_splitgauss(self, x, y): + """Estimation of *Height, Position, FWHM1, FWHM2* of peaks, for + asymmetric gaussian-like curves. + + This functions uses :meth:`estimate_height_position_fwhm`, then + adds a second (identical) estimation of FWHM to the fit parameters + for each peak, and the corresponding constraint. + + :param x: Array of abscissa values + :param y: Array of ordinate values (``y = f(x)``) + :return: Tuple of estimated fit parameters and fit constraints. + Parameters to be estimated for each peak are: + *Height, Position, FWHM1, FWHM2*. + Fit constraints depend on :attr:`config`. + """ + fittedpar, cons = self.estimate_height_position_fwhm(x, y) + # get the number of found peaks + npeaks = len(fittedpar) // 3 + estimated_parameters = [] + estimated_constraints = numpy.zeros((4 * npeaks, 3), numpy.float64) + for i in range(npeaks): + for j in range(3): + estimated_parameters.append(fittedpar[3 * i + j]) + # fwhm2 estimate = fwhm1 + estimated_parameters.append(fittedpar[3 * i + 2]) + # height + estimated_constraints[4 * i, 0] = cons[3 * i, 0] + estimated_constraints[4 * i, 1] = cons[3 * i, 1] + estimated_constraints[4 * i, 2] = cons[3 * i, 2] + # position + estimated_constraints[4 * i + 1, 0] = cons[3 * i + 1, 0] + estimated_constraints[4 * i + 1, 1] = cons[3 * i + 1, 1] + estimated_constraints[4 * i + 1, 2] = cons[3 * i + 1, 2] + # fwhm1 + estimated_constraints[4 * i + 2, 0] = cons[3 * i + 2, 0] + estimated_constraints[4 * i + 2, 1] = cons[3 * i + 2, 1] + estimated_constraints[4 * i + 2, 2] = cons[3 * i + 2, 2] + # fwhm2 + estimated_constraints[4 * i + 3, 0] = cons[3 * i + 2, 0] + estimated_constraints[4 * i + 3, 1] = cons[3 * i + 2, 1] + estimated_constraints[4 * i + 3, 2] = cons[3 * i + 2, 2] + if cons[3 * i + 2, 0] == CFACTOR: + # convert indices of related parameters + # (this happens if SameFwhmFlag == True) + estimated_constraints[4 * i + 2, 1] = ( + int(cons[3 * i + 2, 1] / 3) * 4 + 2 + ) + estimated_constraints[4 * i + 3, 1] = ( + int(cons[3 * i + 2, 1] / 3) * 4 + 3 + ) + return estimated_parameters, estimated_constraints + + def estimate_pvoigt(self, x, y): + """Estimation of *Height, Position, FWHM, eta* of peaks, for + pseudo-Voigt curves. + + Pseudo-Voigt are a sum of a gaussian curve *G(x)* and a lorentzian + curve *L(x)* with the same height, center, fwhm parameters: + ``y(x) = eta * G(x) + (1-eta) * L(x)`` + + This functions uses :meth:`estimate_height_position_fwhm`, then + adds a constant estimation of *eta* (0.5) to the fit parameters + for each peak, and the corresponding constraint. + + :param x: Array of abscissa values + :param y: Array of ordinate values (``y = f(x)``) + :return: Tuple of estimated fit parameters and fit constraints. + Parameters to be estimated for each peak are: + *Height, Position, FWHM, eta*. + Constraint for the eta parameter can be set to QUOTED (0.--1.) + by setting :attr:`config`['QuotedEtaFlag'] to ``True``. + If this is not the case, the constraint code is set to FREE. + """ + fittedpar, cons = self.estimate_height_position_fwhm(x, y) + npeaks = len(fittedpar) // 3 + newpar = [] + newcons = numpy.zeros((4 * npeaks, 3), numpy.float64) + # find out related parameters proper index + if not self.config["NoConstraintsFlag"]: + if self.config["SameFwhmFlag"]: + j = 0 + # get the index of the free FWHM + for i in range(npeaks): + if cons[3 * i + 2, 0] != 4: + j = i + for i in range(npeaks): + if i != j: + cons[3 * i + 2, 1] = 4 * j + 2 + for i in range(npeaks): + newpar.append(fittedpar[3 * i]) + newpar.append(fittedpar[3 * i + 1]) + newpar.append(fittedpar[3 * i + 2]) + newpar.append(0.5) + # height + newcons[4 * i, 0] = cons[3 * i, 0] + newcons[4 * i, 1] = cons[3 * i, 1] + newcons[4 * i, 2] = cons[3 * i, 2] + # position + newcons[4 * i + 1, 0] = cons[3 * i + 1, 0] + newcons[4 * i + 1, 1] = cons[3 * i + 1, 1] + newcons[4 * i + 1, 2] = cons[3 * i + 1, 2] + # fwhm + newcons[4 * i + 2, 0] = cons[3 * i + 2, 0] + newcons[4 * i + 2, 1] = cons[3 * i + 2, 1] + newcons[4 * i + 2, 2] = cons[3 * i + 2, 2] + # Eta constrains + newcons[4 * i + 3, 0] = CFREE + newcons[4 * i + 3, 1] = 0 + newcons[4 * i + 3, 2] = 0 + if self.config["QuotedEtaFlag"]: + newcons[4 * i + 3, 0] = CQUOTED + newcons[4 * i + 3, 1] = 0.0 + newcons[4 * i + 3, 2] = 1.0 + return newpar, newcons + + def estimate_splitpvoigt(self, x, y): + """Estimation of *Height, Position, FWHM1, FWHM2, eta* of peaks, for + asymmetric pseudo-Voigt curves. + + This functions uses :meth:`estimate_height_position_fwhm`, then + adds an identical FWHM2 parameter and a constant estimation of + *eta* (0.5) to the fit parameters for each peak, and the corresponding + constraints. + + Constraint for the eta parameter can be set to QUOTED (0.--1.) + by setting :attr:`config`['QuotedEtaFlag'] to ``True``. + If this is not the case, the constraint code is set to FREE. + + :param x: Array of abscissa values + :param y: Array of ordinate values (``y = f(x)``) + :return: Tuple of estimated fit parameters and fit constraints. + Parameters to be estimated for each peak are: + *Height, Position, FWHM1, FWHM2, eta*. + """ + fittedpar, cons = self.estimate_height_position_fwhm(x, y) + npeaks = len(fittedpar) // 3 + newpar = [] + newcons = numpy.zeros((5 * npeaks, 3), numpy.float64) + # find out related parameters proper index + if not self.config["NoConstraintsFlag"]: + if self.config["SameFwhmFlag"]: + j = 0 + # get the index of the free FWHM + for i in range(npeaks): + if cons[3 * i + 2, 0] != 4: + j = i + for i in range(npeaks): + if i != j: + cons[3 * i + 2, 1] = 4 * j + 2 + for i in range(npeaks): + # height + newpar.append(fittedpar[3 * i]) + # position + newpar.append(fittedpar[3 * i + 1]) + # fwhm1 + newpar.append(fittedpar[3 * i + 2]) + # fwhm2 estimate equal to fwhm1 + newpar.append(fittedpar[3 * i + 2]) + # eta + newpar.append(0.5) + # constraint codes + # ---------------- + # height + newcons[5 * i, 0] = cons[3 * i, 0] + # position + newcons[5 * i + 1, 0] = cons[3 * i + 1, 0] + # fwhm1 + newcons[5 * i + 2, 0] = cons[3 * i + 2, 0] + # fwhm2 + newcons[5 * i + 3, 0] = cons[3 * i + 2, 0] + # cons 1 + # ------ + newcons[5 * i, 1] = cons[3 * i, 1] + newcons[5 * i + 1, 1] = cons[3 * i + 1, 1] + newcons[5 * i + 2, 1] = cons[3 * i + 2, 1] + newcons[5 * i + 3, 1] = cons[3 * i + 2, 1] + # cons 2 + # ------ + newcons[5 * i, 2] = cons[3 * i, 2] + newcons[5 * i + 1, 2] = cons[3 * i + 1, 2] + newcons[5 * i + 2, 2] = cons[3 * i + 2, 2] + newcons[5 * i + 3, 2] = cons[3 * i + 2, 2] + + if cons[3 * i + 2, 0] == CFACTOR: + # fwhm2 connstraint depends on fwhm1 + newcons[5 * i + 3, 1] = newcons[5 * i + 2, 1] + 1 + # eta constraints + newcons[5 * i + 4, 0] = CFREE + newcons[5 * i + 4, 1] = 0 + newcons[5 * i + 4, 2] = 0 + if self.config["QuotedEtaFlag"]: + newcons[5 * i + 4, 0] = CQUOTED + newcons[5 * i + 4, 1] = 0.0 + newcons[5 * i + 4, 2] = 1.0 + return newpar, newcons + + def estimate_splitpvoigt2(self, x, y): + """Estimation of *Height, Position, FWHM1, FWHM2, eta1, eta2* of peaks, for + asymmetric pseudo-Voigt curves. + + This functions uses :meth:`estimate_height_position_fwhm`, then + adds an identical FWHM2 parameter and a constant estimation of + *eta1* and *eta2* (0.5) to the fit parameters for each peak, and the corresponding + constraints. + + Constraint for the eta parameter can be set to QUOTED (0.--1.) + by setting :attr:`config`['QuotedEtaFlag'] to ``True``. + If this is not the case, the constraint code is set to FREE. + + :param x: Array of abscissa values + :param y: Array of ordinate values (``y = f(x)``) + :return: Tuple of estimated fit parameters and fit constraints. + Parameters to be estimated for each peak are: + *Height, Position, FWHM1, FWHM2, eta1, eta2*. + """ + fittedpar, cons = self.estimate_height_position_fwhm(x, y) + npeaks = len(fittedpar) // 3 + newpar = [] + newcons = numpy.zeros((5 * npeaks, 3), numpy.float64) + # find out related parameters proper index + if not self.config["NoConstraintsFlag"]: + if self.config["SameFwhmFlag"]: + j = 0 + # get the index of the free FWHM + for i in range(npeaks): + if cons[3 * i + 2, 0] != 4: + j = i + for i in range(npeaks): + if i != j: + cons[3 * i + 2, 1] = 4 * j + 2 + for i in range(npeaks): + # height + newpar.append(fittedpar[3 * i]) + # position + newpar.append(fittedpar[3 * i + 1]) + # fwhm1 + newpar.append(fittedpar[3 * i + 2]) + # fwhm2 estimate equal to fwhm1 + newpar.append(fittedpar[3 * i + 2]) + # eta1 + newpar.append(0.5) + # eta2 + newpar.append(0.5) + # constraint codes + # ---------------- + # height + newcons[6 * i, 0] = cons[3 * i, 0] + # position + newcons[6 * i + 1, 0] = cons[3 * i + 1, 0] + # fwhm1 + newcons[6 * i + 2, 0] = cons[3 * i + 2, 0] + # fwhm2 + newcons[6 * i + 3, 0] = cons[3 * i + 2, 0] + # cons 1 + # ------ + newcons[6 * i, 1] = cons[3 * i, 1] + newcons[6 * i + 1, 1] = cons[3 * i + 1, 1] + newcons[6 * i + 2, 1] = cons[3 * i + 2, 1] + newcons[6 * i + 3, 1] = cons[3 * i + 2, 1] + # cons 2 + # ------ + newcons[6 * i, 2] = cons[3 * i, 2] + newcons[6 * i + 1, 2] = cons[3 * i + 1, 2] + newcons[6 * i + 2, 2] = cons[3 * i + 2, 2] + newcons[6 * i + 3, 2] = cons[3 * i + 2, 2] + + if cons[3 * i + 2, 0] == CFACTOR: + # fwhm2 constraint depends on fwhm1 + newcons[6 * i + 3, 1] = newcons[6 * i + 2, 1] + 1 + # eta1 constraints + newcons[6 * i + 4, 0] = CFREE + newcons[6 * i + 4, 1] = 0 + newcons[6 * i + 4, 2] = 0 + if self.config["QuotedEtaFlag"]: + newcons[6 * i + 4, 0] = CQUOTED + newcons[6 * i + 4, 1] = 0.0 + newcons[6 * i + 4, 2] = 1.0 + # eta2 constraints + newcons[6 * i + 5, 0] = CFREE + newcons[6 * i + 5, 1] = 0 + newcons[6 * i + 5, 2] = 0 + if self.config["QuotedEtaFlag"]: + newcons[6 * i + 5, 0] = CQUOTED + newcons[6 * i + 5, 1] = 0.0 + newcons[6 * i + 5, 2] = 1.0 + return newpar, newcons + + def estimate_apvoigt(self, x, y): + """Estimation of *Area, Position, FWHM1, eta* of peaks, for + pseudo-Voigt curves. + + This functions uses :meth:`estimate_pvoigt`, then converts the height + parameter to area. + + :param x: Array of abscissa values + :param y: Array of ordinate values (``y = f(x)``) + :return: Tuple of estimated fit parameters and fit constraints. + Parameters to be estimated for each peak are: + *Area, Position, FWHM, eta*. + """ + fittedpar, cons = self.estimate_pvoigt(x, y) + npeaks = len(fittedpar) // 4 + # Assume 50% of the area is determined by the gaussian and 50% by + # the Lorentzian. + for i in range(npeaks): + height = fittedpar[4 * i] + fwhm = fittedpar[4 * i + 2] + fittedpar[4 * i] = 0.5 * (height * fwhm * 0.5 * numpy.pi) + 0.5 * ( + height * fwhm / (2.0 * numpy.sqrt(2 * numpy.log(2))) + ) * numpy.sqrt(2 * numpy.pi) + return fittedpar, cons + + def estimate_ahypermet(self, x, y): + """Estimation of *area, position, fwhm, st_area_r, st_slope_r, + lt_area_r, lt_slope_r, step_height_r* of peaks, for hypermet curves. + + :param x: Array of abscissa values + :param y: Array of ordinate values (``y = f(x)``) + :return: Tuple of estimated fit parameters and fit constraints. + Parameters to be estimated for each peak are: + *area, position, fwhm, st_area_r, st_slope_r, + lt_area_r, lt_slope_r, step_height_r* . + """ + yscaling = self.config.get("Yscaling", 1.0) + if yscaling == 0: + yscaling = 1.0 + fittedpar, cons = self.estimate_height_position_fwhm(x, y) + npeaks = len(fittedpar) // 3 + newpar = [] + newcons = numpy.zeros((8 * npeaks, 3), numpy.float64) + main_peak = 0 + # find out related parameters proper index + if not self.config["NoConstraintsFlag"]: + if self.config["SameFwhmFlag"]: + j = 0 + # get the index of the free FWHM + for i in range(npeaks): + if cons[3 * i + 2, 0] != 4: + j = i + for i in range(npeaks): + if i != j: + cons[3 * i + 2, 1] = 8 * j + 2 + main_peak = j + for i in range(npeaks): + if fittedpar[3 * i] > fittedpar[3 * main_peak]: + main_peak = i + + for i in range(npeaks): + height = fittedpar[3 * i] + position = fittedpar[3 * i + 1] + fwhm = fittedpar[3 * i + 2] + area = (height * fwhm / (2.0 * numpy.sqrt(2 * numpy.log(2)))) * numpy.sqrt( + 2 * numpy.pi + ) + # the gaussian parameters + newpar.append(area) + newpar.append(position) + newpar.append(fwhm) + # print "area, pos , fwhm = ",area,position,fwhm + # Avoid zero derivatives because of not calculating contribution + g_term = 1 + st_term = 1 + lt_term = 1 + step_term = 1 + if self.config["HypermetTails"] != 0: + g_term = self.config["HypermetTails"] & 1 + st_term = (self.config["HypermetTails"] >> 1) & 1 + lt_term = (self.config["HypermetTails"] >> 2) & 1 + step_term = (self.config["HypermetTails"] >> 3) & 1 + if g_term == 0: + # fix the gaussian parameters + newcons[8 * i, 0] = CFIXED + newcons[8 * i + 1, 0] = CFIXED + newcons[8 * i + 2, 0] = CFIXED + # the short tail parameters + if ((area * yscaling) < self.config["MinGaussArea4ShortTail"]) | ( + st_term == 0 + ): + newpar.append(0.0) + newpar.append(0.0) + newcons[8 * i + 3, 0] = CFIXED + newcons[8 * i + 3, 1] = 0.0 + newcons[8 * i + 3, 2] = 0.0 + newcons[8 * i + 4, 0] = CFIXED + newcons[8 * i + 4, 1] = 0.0 + newcons[8 * i + 4, 2] = 0.0 + else: + newpar.append(self.config["InitialShortTailAreaRatio"]) + newpar.append(self.config["InitialShortTailSlopeRatio"]) + newcons[8 * i + 3, 0] = CQUOTED + newcons[8 * i + 3, 1] = self.config["MinShortTailAreaRatio"] + newcons[8 * i + 3, 2] = self.config["MaxShortTailAreaRatio"] + newcons[8 * i + 4, 0] = CQUOTED + newcons[8 * i + 4, 1] = self.config["MinShortTailSlopeRatio"] + newcons[8 * i + 4, 2] = self.config["MaxShortTailSlopeRatio"] + # the long tail parameters + if ((area * yscaling) < self.config["MinGaussArea4LongTail"]) | ( + lt_term == 0 + ): + newpar.append(0.0) + newpar.append(0.0) + newcons[8 * i + 5, 0] = CFIXED + newcons[8 * i + 5, 1] = 0.0 + newcons[8 * i + 5, 2] = 0.0 + newcons[8 * i + 6, 0] = CFIXED + newcons[8 * i + 6, 1] = 0.0 + newcons[8 * i + 6, 2] = 0.0 + else: + newpar.append(self.config["InitialLongTailAreaRatio"]) + newpar.append(self.config["InitialLongTailSlopeRatio"]) + newcons[8 * i + 5, 0] = CQUOTED + newcons[8 * i + 5, 1] = self.config["MinLongTailAreaRatio"] + newcons[8 * i + 5, 2] = self.config["MaxLongTailAreaRatio"] + newcons[8 * i + 6, 0] = CQUOTED + newcons[8 * i + 6, 1] = self.config["MinLongTailSlopeRatio"] + newcons[8 * i + 6, 2] = self.config["MaxLongTailSlopeRatio"] + # the step parameters + if ((height * yscaling) < self.config["MinGaussHeight4StepTail"]) | ( + step_term == 0 + ): + newpar.append(0.0) + newcons[8 * i + 7, 0] = CFIXED + newcons[8 * i + 7, 1] = 0.0 + newcons[8 * i + 7, 2] = 0.0 + else: + newpar.append(self.config["InitialStepTailHeightRatio"]) + newcons[8 * i + 7, 0] = CQUOTED + newcons[8 * i + 7, 1] = self.config["MinStepTailHeightRatio"] + newcons[8 * i + 7, 2] = self.config["MaxStepTailHeightRatio"] + # if self.config['NoConstraintsFlag'] == 1: + # newcons=numpy.zeros((8*npeaks, 3),numpy.float64) + if npeaks > 0: + if g_term: + if self.config["PositiveHeightAreaFlag"]: + for i in range(npeaks): + newcons[8 * i, 0] = CPOSITIVE + if self.config["PositiveFwhmFlag"]: + for i in range(npeaks): + newcons[8 * i + 2, 0] = CPOSITIVE + if self.config["SameFwhmFlag"]: + for i in range(npeaks): + if i != main_peak: + newcons[8 * i + 2, 0] = CFACTOR + newcons[8 * i + 2, 1] = 8 * main_peak + 2 + newcons[8 * i + 2, 2] = 1.0 + if self.config["HypermetQuotedPositionFlag"]: + for i in range(npeaks): + delta = self.config["DeltaPositionFwhmUnits"] * fwhm + newcons[8 * i + 1, 0] = CQUOTED + newcons[8 * i + 1, 1] = newpar[8 * i + 1] - delta + newcons[8 * i + 1, 2] = newpar[8 * i + 1] + delta + if self.config["SameSlopeRatioFlag"]: + for i in range(npeaks): + if i != main_peak: + newcons[8 * i + 4, 0] = CFACTOR + newcons[8 * i + 4, 1] = 8 * main_peak + 4 + newcons[8 * i + 4, 2] = 1.0 + newcons[8 * i + 6, 0] = CFACTOR + newcons[8 * i + 6, 1] = 8 * main_peak + 6 + newcons[8 * i + 6, 2] = 1.0 + if self.config["SameAreaRatioFlag"]: + for i in range(npeaks): + if i != main_peak: + newcons[8 * i + 3, 0] = CFACTOR + newcons[8 * i + 3, 1] = 8 * main_peak + 3 + newcons[8 * i + 3, 2] = 1.0 + newcons[8 * i + 5, 0] = CFACTOR + newcons[8 * i + 5, 1] = 8 * main_peak + 5 + newcons[8 * i + 5, 2] = 1.0 + return newpar, newcons + + def estimate_stepdown(self, x, y): + """Estimation of parameters for stepdown curves. + + The functions estimates gaussian parameters for the derivative of + the data, takes the largest gaussian peak and uses its estimated + parameters to define the center of the step and its fwhm. The + estimated amplitude returned is simply ``max(y) - min(y)``. + + :param x: Array of abscissa values + :param y: Array of ordinate values (``y = f(x)``) + :return: Tuple of estimated fit parameters and fit newconstraints. + Parameters to be estimated for each stepdown are: + *height, centroid, fwhm* . + """ + crappyfilter = [-0.25, -0.75, 0.0, 0.75, 0.25] + cutoff = len(crappyfilter) // 2 + y_deriv = numpy.convolve(y, crappyfilter, mode="valid") + + # make the derivative's peak have the same amplitude as the step + if max(y_deriv) > 0: + y_deriv = y_deriv * max(y) / max(y_deriv) + + fittedpar, newcons = self.estimate_height_position_fwhm( + x[cutoff:-cutoff], y_deriv + ) + + data_amplitude = max(y) - min(y) + + # use parameters from largest gaussian found + if len(fittedpar): + npeaks = len(fittedpar) // 3 + largest_index = 0 + largest = [ + data_amplitude, + fittedpar[3 * largest_index + 1], + fittedpar[3 * largest_index + 2], + ] + for i in range(npeaks): + if fittedpar[3 * i] > largest[0]: + largest_index = i + largest = [ + data_amplitude, + fittedpar[3 * largest_index + 1], + fittedpar[3 * largest_index + 2], + ] + else: + # no peak was found + largest = [ + data_amplitude, # height + x[len(x) // 2], # center: middle of x range + self.config["FwhmPoints"] * (x[1] - x[0]), + ] # fwhm: default value + + # Setup constrains + newcons = numpy.zeros((3, 3), numpy.float64) + if not self.config["NoConstraintsFlag"]: + # Setup height constrains + if self.config["PositiveHeightAreaFlag"]: + newcons[0, 0] = CPOSITIVE + newcons[0, 1] = 0 + newcons[0, 2] = 0 + + # Setup position constrains + if self.config["QuotedPositionFlag"]: + newcons[1, 0] = CQUOTED + newcons[1, 1] = min(x) + newcons[1, 2] = max(x) + + # Setup positive FWHM constrains + if self.config["PositiveFwhmFlag"]: + newcons[2, 0] = CPOSITIVE + newcons[2, 1] = 0 + newcons[2, 2] = 0 + + return largest, newcons + + def estimate_slit(self, x, y): + """Estimation of parameters for slit curves. + + The functions estimates stepup and stepdown parameters for the largest + steps, and uses them for calculating the center (middle between stepup + and stepdown), the height (maximum amplitude in data), the fwhm + (distance between the up- and down-step centers) and the beamfwhm + (average of FWHM for up- and down-step). + + :param x: Array of abscissa values + :param y: Array of ordinate values (``y = f(x)``) + :return: Tuple of estimated fit parameters and fit constraints. + Parameters to be estimated for each slit are: + *height, position, fwhm, beamfwhm* . + """ + largestup, cons = self.estimate_stepup(x, y) + largestdown, cons = self.estimate_stepdown(x, y) + fwhm = numpy.fabs(largestdown[1] - largestup[1]) + beamfwhm = 0.5 * (largestup[2] + largestdown[1]) + beamfwhm = min(beamfwhm, fwhm / 10.0) + beamfwhm = max(beamfwhm, (max(x) - min(x)) * 3.0 / len(x)) + + y_minus_bg = y - self.strip_bg(y) + height = max(y_minus_bg) + + i1 = numpy.nonzero(y_minus_bg >= 0.5 * height)[0] + xx = numpy.take(x, i1) + position = (xx[0] + xx[-1]) / 2.0 + fwhm = xx[-1] - xx[0] + largest = [height, position, fwhm, beamfwhm] + cons = numpy.zeros((4, 3), numpy.float64) + # Setup constrains + if not self.config["NoConstraintsFlag"]: + # Setup height constrains + if self.config["PositiveHeightAreaFlag"]: + cons[0, 0] = CPOSITIVE + cons[0, 1] = 0 + cons[0, 2] = 0 + + # Setup position constrains + if self.config["QuotedPositionFlag"]: + cons[1, 0] = CQUOTED + cons[1, 1] = min(x) + cons[1, 2] = max(x) + + # Setup positive FWHM constrains + if self.config["PositiveFwhmFlag"]: + cons[2, 0] = CPOSITIVE + cons[2, 1] = 0 + cons[2, 2] = 0 + + # Setup positive FWHM constrains + if self.config["PositiveFwhmFlag"]: + cons[3, 0] = CPOSITIVE + cons[3, 1] = 0 + cons[3, 2] = 0 + return largest, cons + + def estimate_stepup(self, x, y): + """Estimation of parameters for a single step up curve. + + The functions estimates gaussian parameters for the derivative of + the data, takes the largest gaussian peak and uses its estimated + parameters to define the center of the step and its fwhm. The + estimated amplitude returned is simply ``max(y) - min(y)``. + + :param x: Array of abscissa values + :param y: Array of ordinate values (``y = f(x)``) + :return: Tuple of estimated fit parameters and fit constraints. + Parameters to be estimated for each stepup are: + *height, centroid, fwhm* . + """ + crappyfilter = [0.25, 0.75, 0.0, -0.75, -0.25] + cutoff = len(crappyfilter) // 2 + y_deriv = numpy.convolve(y, crappyfilter, mode="valid") + if max(y_deriv) > 0: + y_deriv = y_deriv * max(y) / max(y_deriv) + + fittedpar, cons = self.estimate_height_position_fwhm(x[cutoff:-cutoff], y_deriv) + + # for height, use the data amplitude after removing the background + data_amplitude = max(y) - min(y) + + # find params of the largest gaussian found + if len(fittedpar): + npeaks = len(fittedpar) // 3 + largest_index = 0 + largest = [ + data_amplitude, + fittedpar[3 * largest_index + 1], + fittedpar[3 * largest_index + 2], + ] + for i in range(npeaks): + if fittedpar[3 * i] > largest[0]: + largest_index = i + largest = [ + fittedpar[3 * largest_index], + fittedpar[3 * largest_index + 1], + fittedpar[3 * largest_index + 2], + ] + else: + # no peak was found + largest = [ + data_amplitude, # height + x[len(x) // 2], # center: middle of x range + self.config["FwhmPoints"] * (x[1] - x[0]), + ] # fwhm: default value + + newcons = numpy.zeros((3, 3), numpy.float64) + # Setup constrains + if not self.config["NoConstraintsFlag"]: + # Setup height constraints + if self.config["PositiveHeightAreaFlag"]: + newcons[0, 0] = CPOSITIVE + newcons[0, 1] = 0 + newcons[0, 2] = 0 + + # Setup position constraints + if self.config["QuotedPositionFlag"]: + newcons[1, 0] = CQUOTED + newcons[1, 1] = min(x) + newcons[1, 2] = max(x) + + # Setup positive FWHM constraints + if self.config["PositiveFwhmFlag"]: + newcons[2, 0] = CPOSITIVE + newcons[2, 1] = 0 + newcons[2, 2] = 0 + + return largest, newcons + + def estimate_periodic_gauss(self, x, y): + """Estimation of parameters for periodic gaussian curves: + *number of peaks, distance between peaks, height, position of the + first peak, fwhm* + + The functions detects all peaks, then computes the parameters the + following way: + + - *distance*: average of distances between detected peaks + - *height*: average height of detected peaks + - *fwhm*: fwhm of the highest peak (in number of samples) if + field ``'AutoFwhm'`` in :attr:`config` is ``True``, else take + the default value (field ``'FwhmPoints'`` in :attr:`config`) + + :param x: Array of abscissa values + :param y: Array of ordinate values (``y = f(x)``) + :return: Tuple of estimated fit parameters and fit constraints. + """ + yscaling = self.config.get("Yscaling", 1.0) + if yscaling == 0: + yscaling = 1.0 + + bg = self.strip_bg(y) + + if self.config["AutoFwhm"]: + search_fwhm = guess_fwhm(y) + else: + search_fwhm = int(float(self.config["FwhmPoints"])) + search_sens = float(self.config["Sensitivity"]) + + if search_fwhm < 3: + search_fwhm = 3 + + if search_sens < 1: + search_sens = 1 + + if len(y) > 1.5 * search_fwhm: + peaks = peak_search(yscaling * y, fwhm=search_fwhm, sensitivity=search_sens) + else: + peaks = [] + npeaks = len(peaks) + if not npeaks: + fittedpar = [] + cons = numpy.zeros((len(fittedpar), 3), numpy.float64) + return fittedpar, cons + + fittedpar = [0.0, 0.0, 0.0, 0.0, 0.0] + + # The number of peaks + fittedpar[0] = npeaks + + # The separation between peaks in x units + delta = 0.0 + height = 0.0 + for i in range(npeaks): + height += y[int(peaks[i])] - bg[int(peaks[i])] + if i != npeaks - 1: + delta += x[int(peaks[i + 1])] - x[int(peaks[i])] + + # delta between peaks + if npeaks > 1: + fittedpar[1] = delta / (npeaks - 1) + + # starting height + fittedpar[2] = height / npeaks + + # position of the first peak + fittedpar[3] = x[int(peaks[0])] + + # Estimate the fwhm + fittedpar[4] = search_fwhm + + # setup constraints + cons = numpy.zeros((5, 3), numpy.float64) + cons[0, 0] = CFIXED # the number of gaussians + if npeaks == 1: + cons[1, 0] = CFIXED # the delta between peaks + else: + cons[1, 0] = CFREE + j = 2 + # Setup height area constrains + if not self.config["NoConstraintsFlag"]: + if self.config["PositiveHeightAreaFlag"]: + # POSITIVE = 1 + cons[j, 0] = CPOSITIVE + cons[j, 1] = 0 + cons[j, 2] = 0 + j += 1 + + # Setup position constrains + if not self.config["NoConstraintsFlag"]: + if self.config["QuotedPositionFlag"]: + # QUOTED = 2 + cons[j, 0] = CQUOTED + cons[j, 1] = min(x) + cons[j, 2] = max(x) + j += 1 + + # Setup positive FWHM constrains + if not self.config["NoConstraintsFlag"]: + if self.config["PositiveFwhmFlag"]: + # POSITIVE=1 + cons[j, 0] = CPOSITIVE + cons[j, 1] = 0 + cons[j, 2] = 0 + j += 1 + return fittedpar, cons + + def configure(self, **kw): + """Add new / unknown keyword arguments to :attr:`config`, + update entries in :attr:`config` if the parameter name is a existing + key. + + :param kw: Dictionary of keyword arguments. + :return: Configuration dictionary :attr:`config` + """ + if not kw.keys(): + return self.config + for key in kw.keys(): + notdone = 1 + # take care of lower / upper case problems ... + for config_key in self.config.keys(): + if config_key.lower() == key.lower(): + self.config[config_key] = kw[key] + notdone = 0 + if notdone: + self.config[key] = kw[key] + return self.config + + +fitfuns = FitTheories() + +THEORY = dict( + ( + ( + "Gaussians", + FitTheory( + description="Gaussian functions", + function=functions.sum_gauss, + parameters=("Height", "Position", "FWHM"), + estimate=fitfuns.estimate_height_position_fwhm, + configure=fitfuns.configure, + ), + ), + ( + "Lorentz", + FitTheory( + description="Lorentzian functions", + function=functions.sum_lorentz, + parameters=("Height", "Position", "FWHM"), + estimate=fitfuns.estimate_height_position_fwhm, + configure=fitfuns.configure, + ), + ), + ( + "Area Gaussians", + FitTheory( + description="Gaussian functions (area)", + function=functions.sum_agauss, + parameters=("Area", "Position", "FWHM"), + estimate=fitfuns.estimate_agauss, + configure=fitfuns.configure, + ), + ), + ( + "Area Lorentz", + FitTheory( + description="Lorentzian functions (area)", + function=functions.sum_alorentz, + parameters=("Area", "Position", "FWHM"), + estimate=fitfuns.estimate_alorentz, + configure=fitfuns.configure, + ), + ), + ( + "Pseudo-Voigt Line", + FitTheory( + description="Pseudo-Voigt functions", + function=functions.sum_pvoigt, + parameters=("Height", "Position", "FWHM", "Eta"), + estimate=fitfuns.estimate_pvoigt, + configure=fitfuns.configure, + ), + ), + ( + "Area Pseudo-Voigt", + FitTheory( + description="Pseudo-Voigt functions (area)", + function=functions.sum_apvoigt, + parameters=("Area", "Position", "FWHM", "Eta"), + estimate=fitfuns.estimate_apvoigt, + configure=fitfuns.configure, + ), + ), + ( + "Split Gaussian", + FitTheory( + description="Asymmetric gaussian functions", + function=functions.sum_splitgauss, + parameters=("Height", "Position", "LowFWHM", "HighFWHM"), + estimate=fitfuns.estimate_splitgauss, + configure=fitfuns.configure, + ), + ), + ( + "Split Lorentz", + FitTheory( + description="Asymmetric lorentzian functions", + function=functions.sum_splitlorentz, + parameters=("Height", "Position", "LowFWHM", "HighFWHM"), + estimate=fitfuns.estimate_splitgauss, + configure=fitfuns.configure, + ), + ), + ( + "Split Pseudo-Voigt", + FitTheory( + description="Asymmetric pseudo-Voigt functions", + function=functions.sum_splitpvoigt, + parameters=("Height", "Position", "LowFWHM", "HighFWHM", "Eta"), + estimate=fitfuns.estimate_splitpvoigt, + configure=fitfuns.configure, + ), + ), + ( + "Split Pseudo-Voigt 2", + FitTheory( + description="Asymmetric pseudo-Voigt functions", + function=functions.sum_splitpvoigt2, + parameters=( + "Height", + "Position", + "LowFWHM", + "HighFWHM", + "LowEta", + "HighEta", + ), + estimate=fitfuns.estimate_splitpvoigt2, + configure=fitfuns.configure, + ), + ), + ( + "Step Down", + FitTheory( + description="Step down function", + function=functions.sum_stepdown, + parameters=("Height", "Position", "FWHM"), + estimate=fitfuns.estimate_stepdown, + configure=fitfuns.configure, + ), + ), + ( + "Step Up", + FitTheory( + description="Step up function", + function=functions.sum_stepup, + parameters=("Height", "Position", "FWHM"), + estimate=fitfuns.estimate_stepup, + configure=fitfuns.configure, + ), + ), + ( + "Slit", + FitTheory( + description="Slit function", + function=functions.sum_slit, + parameters=("Height", "Position", "FWHM", "BeamFWHM"), + estimate=fitfuns.estimate_slit, + configure=fitfuns.configure, + ), + ), + ( + "Atan", + FitTheory( + description="Arctan step up function", + function=functions.atan_stepup, + parameters=("Height", "Position", "Width"), + estimate=fitfuns.estimate_stepup, + configure=fitfuns.configure, + ), + ), + ( + "Hypermet", + FitTheory( + description="Hypermet functions", + function=fitfuns.ahypermet, # customized version of functions.sum_ahypermet + parameters=( + "G_Area", + "Position", + "FWHM", + "ST_Area", + "ST_Slope", + "LT_Area", + "LT_Slope", + "Step_H", + ), + estimate=fitfuns.estimate_ahypermet, + configure=fitfuns.configure, + ), + ), + # ('Periodic Gaussians', + # FitTheory(description='Periodic gaussian functions', + # function=functions.periodic_gauss, + # parameters=('N', 'Delta', 'Height', 'Position', 'FWHM'), + # estimate=fitfuns.estimate_periodic_gauss, + # configure=fitfuns.configure)) + ( + "Degree 2 Polynomial", + FitTheory( + description="Degree 2 polynomial" "\ny = a*x^2 + b*x +c", + function=fitfuns.poly, + parameters=["a", "b", "c"], + estimate=fitfuns.estimate_quadratic, + ), + ), + ( + "Degree 3 Polynomial", + FitTheory( + description="Degree 3 polynomial" "\ny = a*x^3 + b*x^2 + c*x + d", + function=fitfuns.poly, + parameters=["a", "b", "c", "d"], + estimate=fitfuns.estimate_cubic, + ), + ), + ( + "Degree 4 Polynomial", + FitTheory( + description="Degree 4 polynomial" + "\ny = a*x^4 + b*x^3 + c*x^2 + d*x + e", + function=fitfuns.poly, + parameters=["a", "b", "c", "d", "e"], + estimate=fitfuns.estimate_quartic, + ), + ), + ( + "Degree 5 Polynomial", + FitTheory( + description="Degree 5 polynomial" + "\ny = a*x^5 + b*x^4 + c*x^3 + d*x^2 + e*x + f", + function=fitfuns.poly, + parameters=["a", "b", "c", "d", "e", "f"], + estimate=fitfuns.estimate_quintic, + ), + ), + ) +) +"""Dictionary of fit theories: fit functions and their associated estimation +function, parameters list, configuration function and description. +""" + + +def test(a): + from silx.math.fit import fitmanager + + x = numpy.arange(1000).astype(numpy.float64) + p = [1500, 100.0, 50.0, 1500, 700.0, 50.0] + y_synthetic = functions.sum_gauss(x, *p) + 1 + + fit = fitmanager.FitManager(x, y_synthetic) + fit.addtheory( + "Gaussians", + functions.sum_gauss, + ["Height", "Position", "FWHM"], + a.estimate_height_position_fwhm, + ) + fit.settheory("Gaussians") + fit.setbackground("Linear") + + fit.estimate() + fit.runfit() + + y_fit = fit.gendata() + + print("Fit parameter names: %s" % str(fit.get_names())) + print("Theoretical parameters: %s" % str(numpy.append([1, 0], p))) + print("Fitted parameters: %s" % str(fit.get_fitted_parameters())) + + try: + from silx.gui import qt + from silx.gui.plot import plot1D + + app = qt.QApplication([]) + + # Offset of 1 to see the difference in log scale + plot1D(x, (y_synthetic + 1, y_fit), "Input data + 1, Fit") + + app.exec() + except ImportError: + _logger.warning("Unable to load qt binding, can't plot results.") + + +if __name__ == "__main__": + test(fitfuns) |