code: # coding: utf-8
  #/*##########################################################################
  #
  # Copyright (c) 2004-2020 European Synchrotron Radiation Facility
  #
  # Permission is hereby granted, free of charge, to any person obtaining a copy
  # of this software and associated documentation files (the "Software"), to deal
  # in the Software without restriction, including without limitation the rights
  # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  # copies of the Software, and to permit persons to whom the Software is
  # furnished to do so, subject to the following conditions:
  #
  # The above copyright notice and this permission notice shall be included in
  # all copies or substantial portions of the Software.
  #
  # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  # THE SOFTWARE.
  #
  ########################################################################### */
  """This 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::
  
      Consider using an OrderedDict instead of a regular dictionary, when
      defining your own theory dictionary, if the order matters to you.
      This will likely be the case if you intend to load a selection of
      functions in a GUI such as :class:`silx.gui.fit.FitManager`.
  
  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
  from collections import OrderedDict
  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.,
      '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.,
      '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.float)
          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.
          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. / 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.float)
          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.float)
  
              # 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.float)
                  cons[1:len(param):3, 2] = max(xw) * numpy.ones(
                                                          shape,
                                                          numpy.float)
  
              # 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.float)
          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.float)
          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.float)
          # 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.float)
          # 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_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.float)
          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.float)
          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.float)
          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.float)
          # 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.float)
          # 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.float)
              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.float)
          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 = OrderedDict((
      ('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)),
      ('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.float)
      p = [1500, 100., 50.0,
           1500, 700., 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)
  
apis: numpy.sqrt numpy.zeros_like numpy.ones numpy.fabs numpy.append numpy.zeros numpy.take numpy.concatenate numpy.arange numpy.poly1d numpy.log numpy.array numpy.polyfit numpy.convolve numpy.nonzero