Skip to main content
`Ctrl`+`K`

# Source code for silx.math.fit.fittheories

    # /*##########################################################################
    #
    # 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.asarray(y)

            # 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)