https://github.com/lmfit/lmfit-py
Raw File
Tip revision: 7fec373139313d7cc6c2b56cc3ea98e54decb82c authored by Matt Newville on 26 November 2018, 15:25:27 UTC
Merge pull request #519 from reneeotten/PR512
Tip revision: 7fec373
models.py
"""Module containing built-in fitting models."""
import time

from asteval import Interpreter, get_ast_names
import numpy as np

from . import lineshapes
from .lineshapes import (breit_wigner, damped_oscillator, dho, donaich,
                         expgaussian, exponential, gaussian, linear, lognormal,
                         lorentzian, moffat, parabolic, pearson7, powerlaw,
                         pvoigt, rectangle, skewed_gaussian, skewed_voigt,
                         step, students_t, voigt)
from .model import Model


class DimensionalError(Exception):
    """Raise exception when number of independent variables is not one."""

    pass


def _validate_1d(independent_vars):
    if len(independent_vars) != 1:
        raise DimensionalError(
            "This model requires exactly one independent variable.")


def index_of(arr, val):
    """Return index of array nearest to a value."""
    if val < min(arr):
        return 0
    return np.abs(arr-val).argmin()


def fwhm_expr(model):
    """Return constraint expression for fwhm."""
    fmt = "{factor:.7f}*{prefix:s}sigma"
    return fmt.format(factor=model.fwhm_factor, prefix=model.prefix)


def height_expr(model):
    """Return constraint expression for maximum peak height."""
    fmt = "{factor:.7f}*{prefix:s}amplitude/max(1.e-15, {prefix:s}sigma)"
    return fmt.format(factor=model.height_factor, prefix=model.prefix)


def guess_from_peak(model, y, x, negative, ampscale=1.0, sigscale=1.0):
    """Estimate amp, cen, sigma for a peak, create params."""
    if x is None:
        return 1.0, 0.0, 1.0
    maxy, miny = max(y), min(y)
    maxx, minx = max(x), min(x)
    imaxy = index_of(y, maxy)
    cen = x[imaxy]
    amp = (maxy - miny)*3.0
    sig = (maxx-minx)/6.0

    halfmax_vals = np.where(y > (maxy+miny)/2.0)[0]
    if negative:
        imaxy = index_of(y, miny)
        amp = -(maxy - miny)*3.0
        halfmax_vals = np.where(y < (maxy+miny)/2.0)[0]
    if len(halfmax_vals) > 2:
        sig = (x[halfmax_vals[-1]] - x[halfmax_vals[0]])/2.0
        cen = x[halfmax_vals].mean()
    amp = amp*sig*ampscale
    sig = sig*sigscale

    pars = model.make_params(amplitude=amp, center=cen, sigma=sig)
    pars['%ssigma' % model.prefix].set(min=0.0)
    return pars


def update_param_vals(pars, prefix, **kwargs):
    """Update parameter values with keyword arguments."""
    for key, val in kwargs.items():
        pname = "%s%s" % (prefix, key)
        if pname in pars:
            pars[pname].value = val
    return pars


COMMON_INIT_DOC = """
    Parameters
    ----------
    independent_vars : ['x']
        Arguments to func that are independent variables.
    prefix : str, optional
        String to prepend to parameter names, needed to add two Models that
        have parameter names in common.
    nan_policy : str, optional
        How to handle NaN and missing values in data. Must be one of:
        'raise' (default), 'propagate', or 'omit'. See Notes below.
    missing : str, optional
        Synonym for 'nan_policy' for backward compatibility.
    **kwargs : optional
        Keyword arguments to pass to :class:`Model`.

    Notes
    -----
    1. nan_policy sets what to do when a NaN or missing value is seen in the
    data. Should be one of:

        - 'raise' : Raise a ValueError (default)
        - 'propagate' : do nothing
        - 'omit' : (was 'drop') drop missing data

    2. The `missing` argument is deprecated in lmfit 0.9.8 and will be
    removed in a later version. Use `nan_policy` instead, as it is
    consistent with the Minimizer class.

    """

COMMON_GUESS_DOC = """Guess starting values for the parameters of a model.

    Parameters
    ----------
    data : array_like
        Array of data to use to guess parameter values.
    **kws : optional
        Additional keyword arguments, passed to model function.

    Returns
    -------
    params : Parameters

    """

COMMON_DOC = COMMON_INIT_DOC


class ConstantModel(Model):
    """Constant model, with a single Parameter: ``c``.

    Note that this is 'constant' in the sense of having no dependence on
    the independent variable ``x``, not in the sense of being non-varying.
    To be clear, ``c`` will be a Parameter that will be varied
    in the fit (by default, of course).

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})

        def constant(x, c=0.0):
            return c
        super(ConstantModel, self).__init__(constant, **kwargs)

    def guess(self, data, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = self.make_params()

        pars['%sc' % self.prefix].set(value=data.mean())
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class ComplexConstantModel(Model):
    """Complex constant model, with wo Parameters: ``re``, and ``im``.

    Note that ``re`` and ``im`` are 'constant' in the sense of having no
    dependence on the independent variable ``x``, not in the sense of
    being non-varying. To be clear, ``re`` and ``im`` will be Parameters
    that will be varied in the fit (by default, of course).

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 name=None, **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})

        def constant(x, re=0., im=0.):
            return re + 1j*im
        super(ComplexConstantModel, self).__init__(constant, **kwargs)

    def guess(self, data, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = self.make_params()
        pars['%sre' % self.prefix].set(value=data.real.mean())
        pars['%sim' % self.prefix].set(value=data.imag.mean())
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class LinearModel(Model):
    """Linear model, with two Parameters ``intercept`` and ``slope``.

    Defined as:

    .. math::

        f(x; m, b) = m x + b

    with  ``slope`` for :math:`m` and  ``intercept`` for :math:`b`.

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(LinearModel, self).__init__(linear, **kwargs)

    def guess(self, data, x=None, **kwargs):
        """Estimate initial model parameter values from data."""
        sval, oval = 0., 0.
        if x is not None:
            sval, oval = np.polyfit(x, data, 1)
        pars = self.make_params(intercept=oval, slope=sval)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class QuadraticModel(Model):
    """A quadratic model, with three Parameters ``a``, ``b``, and ``c``.

    Defined as:

    .. math::

        f(x; a, b, c) = a x^2 + b x + c

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(QuadraticModel, self).__init__(parabolic, **kwargs)

    def guess(self, data, x=None, **kwargs):
        """Estimate initial model parameter values from data."""
        a, b, c = 0., 0., 0.
        if x is not None:
            a, b, c = np.polyfit(x, data, 2)
        pars = self.make_params(a=a, b=b, c=c)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


ParabolicModel = QuadraticModel


class PolynomialModel(Model):
    r"""A polynomial model with up to 7 Parameters, specfied by ``degree``.

    .. math::

        f(x; c_0, c_1, \ldots, c_7) = \sum_{i=0, 7} c_i  x^i

    with parameters ``c0``, ``c1``, ..., ``c7``.  The supplied ``degree``
    will specify how many of these are actual variable parameters.  This
    uses :numpydoc:`polyval` for its calculation of the polynomial.

    """

    MAX_DEGREE = 7
    DEGREE_ERR = "degree must be an integer less than %d."

    def __init__(self, degree, independent_vars=['x'], prefix='',
                 nan_policy='raise', **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        if not isinstance(degree, int) or degree > self.MAX_DEGREE:
            raise TypeError(self.DEGREE_ERR % self.MAX_DEGREE)

        self.poly_degree = degree
        pnames = ['c%i' % (i) for i in range(degree + 1)]
        kwargs['param_names'] = pnames

        def polynomial(x, c0=0, c1=0, c2=0, c3=0, c4=0, c5=0, c6=0, c7=0):
            return np.polyval([c7, c6, c5, c4, c3, c2, c1, c0], x)

        super(PolynomialModel, self).__init__(polynomial, **kwargs)

    def guess(self, data, x=None, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = self.make_params()
        if x is not None:
            out = np.polyfit(x, data, self.poly_degree)
            for i, coef in enumerate(out[::-1]):
                pars['%sc%i' % (self.prefix, i)].set(value=coef)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class GaussianModel(Model):
    r"""A model based on a Gaussian or normal distribution lineshape (see
    https://en.wikipedia.org/wiki/Normal_distribution), with three Parameters:
    ``amplitude``, ``center``, and ``sigma``.
    In addition, parameters ``fwhm`` and ``height`` are included as constraints
    to report full width at half maximum and maximum peak height, respectively.

    .. math::

        f(x; A, \mu, \sigma) = \frac{A}{\sigma\sqrt{2\pi}} e^{[{-{(x-\mu)^2}/{{2\sigma}^2}}]}

    where the parameter ``amplitude`` corresponds to :math:`A`, ``center`` to
    :math:`\mu`, and ``sigma`` to :math:`\sigma`.  The full width at
    half maximum is :math:`2\sigma\sqrt{2\ln{2}}`, approximately
    :math:`2.3548\sigma`.

    """

    fwhm_factor = 2*np.sqrt(2*np.log(2))
    height_factor = 1./np.sqrt(2*np.pi)

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(GaussianModel, self).__init__(gaussian, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('sigma', min=0)
        self.set_param_hint('fwhm', expr=fwhm_expr(self))
        self.set_param_hint('height', expr=height_expr(self))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class LorentzianModel(Model):
    r"""A model based on a Lorentzian or Cauchy-Lorentz distribution function
    (see https://en.wikipedia.org/wiki/Cauchy_distribution), with three Parameters:
    ``amplitude``, ``center``, and ``sigma``.
    In addition, parameters ``fwhm`` and ``height`` are included as constraints
    to report full width at half maximum and maximum peak height, respectively.

    .. math::

        f(x; A, \mu, \sigma) = \frac{A}{\pi} \big[\frac{\sigma}{(x - \mu)^2 + \sigma^2}\big]

    where the parameter ``amplitude`` corresponds to :math:`A`, ``center`` to
    :math:`\mu`, and ``sigma`` to :math:`\sigma`.  The full width at
    half maximum is :math:`2\sigma`.

    """

    fwhm_factor = 2.0
    height_factor = 1./np.pi

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(LorentzianModel, self).__init__(lorentzian, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('sigma', min=0)
        self.set_param_hint('fwhm', expr=fwhm_expr(self))
        self.set_param_hint('height', expr=height_expr(self))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative, ampscale=1.25)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class VoigtModel(Model):
    r"""A model based on a Voigt distribution function (see
    https://en.wikipedia.org/wiki/Voigt_profile), with four Parameters:
    ``amplitude``, ``center``, ``sigma``, and ``gamma``.  By default,
    ``gamma`` is constrained to have a value equal to ``sigma``, though it
    can be varied independently.  In addition, parameters ``fwhm`` and
    ``height`` are included as constraints to report full width at half
    maximum and maximum peak height, respectively.  The definition for the
    Voigt function used here is

    .. math::

        f(x; A, \mu, \sigma, \gamma) = \frac{A \textrm{Re}[w(z)]}{\sigma\sqrt{2 \pi}}

    where

    .. math::
        :nowrap:

        \begin{eqnarray*}
            z &=& \frac{x-\mu +i\gamma}{\sigma\sqrt{2}} \\
            w(z) &=& e^{-z^2}{\operatorname{erfc}}(-iz)
        \end{eqnarray*}

    and :func:`erfc` is the complementary error function.  As above,
    ``amplitude`` corresponds to :math:`A`, ``center`` to
    :math:`\mu`, and ``sigma`` to :math:`\sigma`. The parameter ``gamma``
    corresponds  to :math:`\gamma`.
    If ``gamma`` is kept at the default value (constrained to ``sigma``),
    the full width at half maximum is approximately :math:`3.6013\sigma`.

    """
    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(VoigtModel, self).__init__(voigt, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('sigma', min=0)
        self.set_param_hint('gamma', expr='%ssigma' % self.prefix)

        fexpr = ("1.0692*{pre:s}gamma+" +
                 "sqrt(0.8664*{pre:s}gamma**2+5.545083*{pre:s}sigma**2)")
        hexpr = ("({pre:s}amplitude/({pre:s}sigma*sqrt(2*pi)))*"
                 "wofz((1j*{pre:s}gamma)/({pre:s}sigma*sqrt(2))).real")

        self.set_param_hint('fwhm', expr=fexpr.format(pre=self.prefix))
        self.set_param_hint('height', expr=hexpr.format(pre=self.prefix))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative,
                               ampscale=1.5, sigscale=0.65)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class PseudoVoigtModel(Model):
    r"""A model based on a pseudo-Voigt distribution function
    (see https://en.wikipedia.org/wiki/Voigt_profile#Pseudo-Voigt_Approximation),
    which is a weighted sum of a Gaussian and Lorentzian distribution function
    that share values for ``amplitude`` (:math:`A`), ``center`` (:math:`\mu`)
    and full width at half maximum ``fwhm`` (and so have  constrained values of
    ``sigma`` (:math:`\sigma`) and ``height`` (maximum peak height).
    A parameter ``fraction`` (:math:`\alpha`) controls the relative weight of
    the Gaussian and Lorentzian components, giving the full definition of

    .. math::

        f(x; A, \mu, \sigma, \alpha) = \frac{(1-\alpha)A}{\sigma_g\sqrt{2\pi}}
        e^{[{-{(x-\mu)^2}/{{2\sigma_g}^2}}]}
        + \frac{\alpha A}{\pi} \big[\frac{\sigma}{(x - \mu)^2 + \sigma^2}\big]

    where :math:`\sigma_g = {\sigma}/{\sqrt{2\ln{2}}}` so that the full width
    at half maximum of each component and of the sum is :math:`2\sigma`. The
    :meth:`guess` function always sets the starting value for ``fraction`` at 0.5.

    """

    fwhm_factor = 2.0

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(PseudoVoigtModel, self).__init__(pvoigt, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('sigma', min=0)
        self.set_param_hint('fraction', value=0.5, min=0.0, max=1.0)
        self.set_param_hint('fwhm', expr=fwhm_expr(self))
        fmt = ("(((1-{prefix:s}fraction)*{prefix:s}amplitude)/"
               "({prefix:s}sigma*sqrt(pi/log(2)))+"
               "({prefix:s}fraction*{prefix:s}amplitude)/"
               "(pi*{prefix:s}sigma))")
        self.set_param_hint('height', expr=fmt.format(prefix=self.prefix))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative, ampscale=1.25)
        pars['%sfraction' % self.prefix].set(value=0.5, min=0.0, max=1.0)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class MoffatModel(Model):
    r"""A model based on the Moffat distribution function
    (see https://en.wikipedia.org/wiki/Moffat_distribution), with four Parameters:
    ``amplitude`` (:math:`A`), ``center`` (:math:`\mu`), a width parameter
    ``sigma`` (:math:`\sigma`) and an exponent ``beta`` (:math:`\beta`).
    In addition, parameters ``fwhm`` and ``height`` are included as constraints
    to report full width at half maximum and maximum peak height, respectively.

    .. math::

        f(x; A, \mu, \sigma, \beta) = A \big[(\frac{x-\mu}{\sigma})^2+1\big]^{-\beta}

    the full width have maximum is :math:`2\sigma\sqrt{2^{1/\beta}-1}`.
    The :meth:`guess` function always sets the starting value for ``beta`` to 1.

    Note that for (:math:`\beta=1`) the Moffat has a Lorentzian shape.

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(MoffatModel, self).__init__(moffat, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('sigma', min=0)
        self.set_param_hint('beta')
        self.set_param_hint('fwhm', expr="2*%ssigma*sqrt(2**(1.0/%sbeta)-1)" % (self.prefix, self.prefix))
        self.set_param_hint('height', expr="%samplitude" % self.prefix)

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative, ampscale=0.5, sigscale=1.)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class Pearson7Model(Model):
    r"""A model based on a Pearson VII distribution (see
    https://en.wikipedia.org/wiki/Pearson_distribution#The_Pearson_type_VII_distribution),
    with four parameters: ``amplitude`` (:math:`A`), ``center``
    (:math:`\mu`), ``sigma`` (:math:`\sigma`), and ``exponent`` (:math:`m`).
    In addition, parameters ``fwhm`` and ``height`` are included as constraints
    to report full width at half maximum and maximum peak height, respectively.

    .. math::

        f(x; A, \mu, \sigma, m) = \frac{A}{\sigma{\beta(m-\frac{1}{2}, \frac{1}{2})}} \bigl[1 + \frac{(x-\mu)^2}{\sigma^2}  \bigr]^{-m}

    where :math:`\beta` is the beta function (see :scipydoc:`special.beta`)
    The :meth:`guess` function always gives a starting value for ``exponent``
    of 1.5.  In addition, parameters ``fwhm`` and ``height`` are included as
    constraints to report full width at half maximum and maximum peak height,
    respectively.

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(Pearson7Model, self).__init__(pearson7, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('expon', value=1.5, max=100)
        fmt = ("sqrt(2**(1/{prefix:s}expon)-1)*2*{prefix:s}sigma")
        self.set_param_hint('fwhm', expr=fmt.format(prefix=self.prefix))
        fmt = ("{prefix:s}amplitude * gamfcn({prefix:s}expon)/"
               "(gamfcn(0.5)*gamfcn({prefix:s}expon-0.5)*{prefix:s}sigma)")
        self.set_param_hint('height', expr=fmt.format(prefix=self.prefix))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative)
        pars['%sexpon' % self.prefix].set(value=1.5)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class StudentsTModel(Model):
    r"""A model based on a Student's t distribution function (see
    https://en.wikipedia.org/wiki/Student%27s_t-distribution), with three Parameters:
    ``amplitude`` (:math:`A`), ``center`` (:math:`\mu`) and ``sigma`` (:math:`\sigma`).
    In addition, parameters ``fwhm`` and ``height`` are included as constraints
    to report full width at half maximum and maximum peak height, respectively.

    .. math::

        f(x; A, \mu, \sigma) = \frac{A \Gamma(\frac{\sigma+1}{2})} {\sqrt{\sigma\pi}\,\Gamma(\frac{\sigma}{2})} \Bigl[1+\frac{(x-\mu)^2}{\sigma}\Bigr]^{-\frac{\sigma+1}{2}}


    where :math:`\Gamma(x)` is the gamma function.

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(StudentsTModel, self).__init__(students_t, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('sigma', min=0.0, max=100)
        fmt = ("{prefix:s}amplitude*gamfcn(({prefix:s}sigma+1)/2)/"
               "(sqrt({prefix:s}sigma*pi)*gamfcn({prefix:s}sigma/2))")
        self.set_param_hint('height', expr=fmt.format(prefix=self.prefix))
        fmt = ("2*sqrt(2**(2/({prefix:s}sigma+1))*"
               "{prefix:s}sigma-{prefix:s}sigma)")
        self.set_param_hint('fwhm', expr=fmt.format(prefix=self.prefix))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class BreitWignerModel(Model):
    r"""A model based on a Breit-Wigner-Fano function (see
    https://en.wikipedia.org/wiki/Fano_resonance), with four Parameters:
    ``amplitude`` (:math:`A`), ``center`` (:math:`\mu`),
    ``sigma`` (:math:`\sigma`), and ``q`` (:math:`q`).
    In addition, parameters ``fwhm`` and ``height`` are included as constraints
    to report full width at half maximum and maximum peak height, respectively.

    .. math::

        f(x; A, \mu, \sigma, q) = \frac{A (q\sigma/2 + x - \mu)^2}{(\sigma/2)^2 + (x - \mu)^2}

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(BreitWignerModel, self).__init__(breit_wigner, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('sigma', min=0.0)
        fmt = ("{prefix:s}amplitude*{prefix:s}q**2")
        self.set_param_hint('height', expr=fmt.format(prefix=self.prefix))
        fmt = ("2*(sqrt({prefix:s}q**2*{prefix:s}sigma**2*({prefix:s}q**2+2))/"
               "max(1.e-15, 2*({prefix:s}q**2)-2))")
        self.set_param_hint('fwhm', expr=fmt.format(prefix=self.prefix))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative)
        pars['%sq' % self.prefix].set(value=1.0)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class LognormalModel(Model):
    r"""A model based on the Log-normal distribution function
    (see https://en.wikipedia.org/wiki/Lognormal), with three Parameters
    ``amplitude`` (:math:`A`), ``center`` (:math:`\mu`) and ``sigma``
    (:math:`\sigma`).
    In addition, parameters ``fwhm`` and ``height`` are included as constraints
    to report full width at half maximum and maximum peak height, respectively.

    .. math::

        f(x; A, \mu, \sigma) = \frac{A}{\sigma\sqrt{2\pi}}\frac{e^{-(\ln(x) - \mu)^2/ 2\sigma^2}}{x}

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(LognormalModel, self).__init__(lognormal, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('center', min=1.e-19)
        self.set_param_hint('sigma', min=0)

        fmt = ("{prefix:s}amplitude/({prefix:s}sigma*sqrt(2*pi))"
               "*exp({prefix:s}sigma**2/2-{prefix:s}center)")
        self.set_param_hint('height', expr=fmt.format(prefix=self.prefix))
        fmt = ("exp({prefix:s}center-{prefix:s}sigma**2+{prefix:s}sigma*sqrt("
               "2*log(2)))-"
               "exp({prefix:s}center-{prefix:s}sigma**2-{prefix:s}sigma*sqrt("
               "2*log(2)))")
        self.set_param_hint('fwhm', expr=fmt.format(prefix=self.prefix))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = self.make_params(amplitude=1.0, center=0.0, sigma=0.25)
        pars['%ssigma' % self.prefix].set(min=0.0)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class DampedOscillatorModel(Model):
    r"""A model based on the Damped Harmonic Oscillator Amplitude
    (see https://en.wikipedia.org/wiki/Harmonic_oscillator#Amplitude_part), with
    three Parameters:  ``amplitude`` (:math:`A`), ``center`` (:math:`\mu`) and
    ``sigma`` (:math:`\sigma`).
    In addition, parameters ``fwhm`` and ``height`` are included as constraints
    to report full width at half maximum and maximum peak height, respectively.

    .. math::

        f(x; A, \mu, \sigma) = \frac{A}{\sqrt{ [1 - (x/\mu)^2]^2 + (2\sigma x/\mu)^2}}

    """

    height_factor = 0.5

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(DampedOscillatorModel, self).__init__(damped_oscillator, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('sigma', min=0)
        self.set_param_hint('height', expr=height_expr(self))
        fmt = ("sqrt(abs({prefix:s}center**2*(1-2*{prefix:s}sigma**2)+"
               "(2*sqrt({prefix:s}center**4*{prefix:s}sigma**2*"
               "({prefix:s}sigma**2+3)))))-"
               "sqrt(abs({prefix:s}center**2*(1-2*{prefix:s}sigma**2)-"
               "(2*sqrt({prefix:s}center**4*{prefix:s}sigma**2*"
               "({prefix:s}sigma**2+3)))))")
        self.set_param_hint('fwhm', expr=fmt.format(prefix=self.prefix))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative,
                               ampscale=0.1, sigscale=0.1)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class DampedHarmonicOscillatorModel(Model):
    r"""A model based on a variation of the Damped Harmonic Oscillator (see
    https://en.wikipedia.org/wiki/Harmonic_oscillator), following the
    definition given in DAVE/PAN (see https://www.ncnr.nist.gov/dave/) with
    four Parameters: ``amplitude`` (:math:`A`), ``center`` (:math:`\mu`),
    ``sigma`` (:math:`\sigma`), and ``gamma`` (:math:`\gamma`).
    In addition, parameters ``fwhm`` and ``height`` are included as constraints
    to report full width at half maximum and maximum peak height, respectively.

    .. math::

        f(x; A, \mu, \sigma, \gamma) = \frac{A\sigma}{\pi [1 - \exp(-x/\gamma)]}
                \Big[ \frac{1}{(x-\mu)^2 + \sigma^2} - \frac{1}{(x+\mu)^2 + \sigma^2} \Big]

    where :math:`\gamma=kT` k is the Boltzmann constant in :math:`evK^-1`
    and T is the temperature in :math:`K`.

    """

    fwhm_factor = 2.0

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(DampedHarmonicOscillatorModel, self).__init__(dho, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('sigma', min=0)
        self.set_param_hint('gamma', min=1.e-19)
        fmt = ("({prefix:s}amplitude*{prefix:s}sigma)/"
               "(pi*(1-exp(-{prefix:s}center/{prefix:s}gamma)))*("
               "1/{prefix:s}sigma**2-1/(4*{prefix:s}center**2+"
               "{prefix:s}sigma**2))")
        self.set_param_hint('height', expr=fmt.format(prefix=self.prefix))
        self.set_param_hint('fwhm', expr=fwhm_expr(self))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative,
                               ampscale=0.1, sigscale=0.1)
        pars['%sgamma' % self.prefix].set(value=1.0, min=0.0)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class ExponentialGaussianModel(Model):
    r"""A model of an Exponentially modified Gaussian distribution
    (see https://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution) with
    four Parameters ``amplitude`` (:math:`A`), ``center`` (:math:`\mu`),
    ``sigma`` (:math:`\sigma`), and  ``gamma`` (:math:`\gamma`).
    In addition, parameters ``fwhm`` and ``height`` are included as constraints
    to report full width at half maximum and maximum peak height, respectively.

    .. math::

        f(x; A, \mu, \sigma, \gamma) = \frac{A\gamma}{2}
        \exp\bigl[\gamma({\mu - x  + \gamma\sigma^2/2})\bigr]
        {\operatorname{erfc}}\Bigl(\frac{\mu + \gamma\sigma^2 - x}{\sqrt{2}\sigma}\Bigr)


    where :func:`erfc` is the complementary error function.

    """

    fwhm_factor = 2*np.sqrt(2*np.log(2))

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(ExponentialGaussianModel, self).__init__(expgaussian, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('sigma', min=0)
        self.set_param_hint('gamma', min=0, max=20)
        fmt = ("{prefix:s}amplitude*{prefix:s}gamma/2*"
               "exp({prefix:s}gamma**2*{prefix:s}sigma**2/2)*"
               "erfc({prefix:s}gamma*{prefix:s}sigma/sqrt(2))")
        self.set_param_hint('height', expr=fmt.format(prefix=self.prefix))
        self.set_param_hint('fwhm', expr=fwhm_expr(self))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class SkewedGaussianModel(Model):
    r"""A variation of the Exponential Gaussian, this uses a skewed normal distribution
    (see https://en.wikipedia.org/wiki/Skew_normal_distribution), with Parameters
    ``amplitude`` (:math:`A`), ``center`` (:math:`\mu`),  ``sigma`` (:math:`\sigma`),
    and ``gamma`` (:math:`\gamma`).
    In addition, parameters ``fwhm`` and ``height`` are included as constraints
    to report full width at half maximum and maximum peak height, respectively.

    .. math::

       f(x; A, \mu, \sigma, \gamma) = \frac{A}{\sigma\sqrt{2\pi}}
       e^{[{-{(x-\mu)^2}/{{2\sigma}^2}}]} \Bigl\{ 1 +
       {\operatorname{erf}}\bigl[
       \frac{{\gamma}(x-\mu)}{\sigma\sqrt{2}}
       \bigr] \Bigr\}

    where :func:`erf` is the error function.

    """

    fwhm_factor = 2*np.sqrt(2*np.log(2))
    height_factor = 1./np.sqrt(2*np.pi)

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(SkewedGaussianModel, self).__init__(skewed_gaussian, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('sigma', min=0)
        self.set_param_hint('height', expr=height_expr(self))
        self.set_param_hint('fwhm', expr=fwhm_expr(self))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class SkewedVoigtModel(Model):
    r"""A variation of the Skewed Gaussian, this applies the same skewing
      to a Voigt distribution (see
      https://en.wikipedia.org/wiki/Voigt_distribution).  It has Parameters
      ``amplitude`` (:math:`A`), ``center`` (:math:`\mu`), ``sigma``
      (:math:`\sigma`), and ``gamma`` (:math:`\gamma`), as usual for a
      Voigt distribution, and add a Parameter ``skew``.  In addition,
      parameters ``fwhm`` and ``height`` are included as constraints to
      report full width at half maximum and maximum peak height, of the
      Voigt distribution, respectively.

    .. math::


       f(x; A, \mu, \sigma, \gamma, \rm{skew}) = {\rm{Voigt}}(x; A, \mu, \sigma, \gamma)
       \Bigl\{ 1 +  {\operatorname{erf}}\bigl[
       \frac{{\rm{skew}}(x-\mu)}{\sigma\sqrt{2}}
       \bigr] \Bigr\}

    where :func:`erf` is the error function.

    """

    fwhm_factor = 3.60131
    height_factor = 1./np.sqrt(2*np.pi)

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(SkewedVoigtModel, self).__init__(skewed_voigt, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('sigma', min=0)
        self.set_param_hint('gamma', expr='%ssigma' % self.prefix)
        self.set_param_hint('height', expr=height_expr(self))
        self.set_param_hint('fwhm', expr=fwhm_expr(self))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class DonaichModel(Model):
    r"""A model of an Doniach Sunjic asymmetric lineshape
    (see https://www.casaxps.com/help_manual/line_shapes.htm), used in
    photo-emission, with four Parameters ``amplitude`` (:math:`A`),
    ``center`` (:math:`\mu`), ``sigma`` (:math:`\sigma`), and ``gamma``
    (:math:`\gamma`).
    In addition, parameter ``height`` is included as a constraint.

    .. math::

        f(x; A, \mu, \sigma, \gamma) = \frac{A}{\sigma^{1-\gamma}}
        \frac{\cos\bigl[\pi\gamma/2 + (1-\gamma)
        \arctan{((x - \mu)}/\sigma)\bigr]} {\bigr[1 + (x-\mu)/\sigma\bigl]^{(1-\gamma)/2}}

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(DonaichModel, self).__init__(donaich, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        fmt = ("{prefix:s}amplitude/({prefix:s}sigma**(1-{prefix:s}gamma))"
               "*cos(pi*{prefix:s}gamma/2)")
        self.set_param_hint('height', expr=fmt.format(prefix=self.prefix))

    def guess(self, data, x=None, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative, ampscale=0.5)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class PowerLawModel(Model):
    r"""A model based on a Power Law (see https://en.wikipedia.org/wiki/Power_law),
    with two Parameters: ``amplitude`` (:math:`A`), and ``exponent`` (:math:`k`), in:

    .. math::

        f(x; A, k) = A x^k

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(PowerLawModel, self).__init__(powerlaw, **kwargs)

    def guess(self, data, x=None, **kwargs):
        """Estimate initial model parameter values from data."""
        try:
            expon, amp = np.polyfit(np.log(x+1.e-14), np.log(data+1.e-14), 1)
        except TypeError:
            expon, amp = 1, np.log(abs(max(data)+1.e-9))

        pars = self.make_params(amplitude=np.exp(amp), exponent=expon)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class ExponentialModel(Model):
    r"""A model based on an exponential decay function
    (see https://en.wikipedia.org/wiki/Exponential_decay) with two Parameters:
    ``amplitude`` (:math:`A`), and ``decay`` (:math:`\tau`), in:

    .. math::

        f(x; A, \tau) = A e^{-x/\tau}

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super(ExponentialModel, self).__init__(exponential, **kwargs)

    def guess(self, data, x=None, **kwargs):
        """Estimate initial model parameter values from data."""
        try:
            sval, oval = np.polyfit(x, np.log(abs(data)+1.e-15), 1)
        except TypeError:
            sval, oval = 1., np.log(abs(max(data)+1.e-9))
        pars = self.make_params(amplitude=np.exp(oval), decay=-1.0/sval)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class StepModel(Model):
    r"""A model based on a Step function, with three Parameters:
    ``amplitude`` (:math:`A`), ``center`` (:math:`\mu`) and ``sigma`` (:math:`\sigma`).

    There are four choices for ``form``:

    - ``linear`` (the default)
    - ``atan`` or ``arctan`` for an arc-tangent function
    - ``erf`` for an error function
    - ``logistic`` for a logistic function (see https://en.wikipedia.org/wiki/Logistic_function)

    The step function starts with a value 0, and ends with a value of
    :math:`A` rising to :math:`A/2` at :math:`\mu`, with :math:`\sigma`
    setting the characteristic width. The functional forms are defined as:

    .. math::
        :nowrap:

        \begin{eqnarray*}
        & f(x; A, \mu, \sigma, {\mathrm{form={}'linear{}'}})  & = A \min{[1, \max{(0,  \alpha)}]} \\
        & f(x; A, \mu, \sigma, {\mathrm{form={}'arctan{}'}})  & = A [1/2 + \arctan{(\alpha)}/{\pi}] \\
        & f(x; A, \mu, \sigma, {\mathrm{form={}'erf{}'}})     & = A [1 + {\operatorname{erf}}(\alpha)]/2 \\
        & f(x; A, \mu, \sigma, {\mathrm{form={}'logistic{}'}})& = A [1 - \frac{1}{1 +  e^{\alpha}} ]
        \end{eqnarray*}

    where :math:`\alpha  = (x - \mu)/{\sigma}`.

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 form='linear', **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'form': form, 'independent_vars': independent_vars})
        super(StepModel, self).__init__(step, **kwargs)

    def guess(self, data, x=None, **kwargs):
        """Estimate initial model parameter values from data."""
        if x is None:
            return
        ymin, ymax = min(data), max(data)
        xmin, xmax = min(x), max(x)
        pars = self.make_params(amplitude=(ymax-ymin),
                                center=(xmax+xmin)/2.0)
        pars['%ssigma' % self.prefix].set(value=(xmax-xmin)/7.0, min=0.0)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class RectangleModel(Model):
    r"""A model based on a Step-up and Step-down function, with five
    Parameters: ``amplitude`` (:math:`A`), ``center1`` (:math:`\mu_1`),
    ``center2`` (:math:`\mu_2`), `sigma1`` (:math:`\sigma_1`) and
    ``sigma2`` (:math:`\sigma_2`).

    There are four choices for ``form``, which is used for both the Step up
    and the Step down:

    - ``linear`` (the default)
    - ``atan`` or ``arctan`` for an arc-tangent function
    - ``erf`` for an error function
    - ``logistic`` for a logistic function (see https://en.wikipedia.org/wiki/Logistic_function)

    The function starts with a value 0, transitions to a value of
    :math:`A`, taking the value :math:`A/2` at :math:`\mu_1`, with :math:`\sigma_1`
    setting the characteristic width. The function then transitions again to
    the value :math:`A/2` at :math:`\mu_2`, with :math:`\sigma_2` setting the
    characteristic width. The functional forms are defined as:

    .. math::
        :nowrap:

        \begin{eqnarray*}
        &f(x; A, \mu, \sigma, {\mathrm{form={}'linear{}'}})   &= A \{ \min{[1, \max{(0, \alpha_1)}]} + \min{[-1, \max{(0,  \alpha_2)}]} \} \\
        &f(x; A, \mu, \sigma, {\mathrm{form={}'arctan{}'}})   &= A [\arctan{(\alpha_1)} + \arctan{(\alpha_2)}]/{\pi} \\
        &f(x; A, \mu, \sigma, {\mathrm{form={}'erf{}'}})      &= A [{\operatorname{erf}}(\alpha_1) + {\operatorname{erf}}(\alpha_2)]/2 \\
        &f(x; A, \mu, \sigma, {\mathrm{form={}'logistic{}'}}) &= A [1 - \frac{1}{1 + e^{\alpha_1}} - \frac{1}{1 +  e^{\alpha_2}} ]
        \end{eqnarray*}


    where :math:`\alpha_1  = (x - \mu_1)/{\sigma_1}` and
    :math:`\alpha_2  = -(x - \mu_2)/{\sigma_2}`.

    """

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 form='linear', **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'form': form, 'independent_vars': independent_vars})
        super(RectangleModel, self).__init__(rectangle, **kwargs)

        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        self.set_param_hint('center1')
        self.set_param_hint('center2')
        self.set_param_hint('midpoint',
                            expr='(%scenter1+%scenter2)/2.0' % (self.prefix,
                                                                self.prefix))

    def guess(self, data, x=None, **kwargs):
        """Estimate initial model parameter values from data."""
        if x is None:
            return
        ymin, ymax = min(data), max(data)
        xmin, xmax = min(x), max(x)
        pars = self.make_params(amplitude=(ymax-ymin),
                                center1=(xmax+xmin)/4.0,
                                center2=3*(xmax+xmin)/4.0)
        pars['%ssigma1' % self.prefix].set(value=(xmax-xmin)/7.0, min=0.0)
        pars['%ssigma2' % self.prefix].set(value=(xmax-xmin)/7.0, min=0.0)
        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC


class ExpressionModel(Model):

    idvar_missing = "No independent variable found in\n %s"
    idvar_notfound = "Cannot find independent variables '%s' in\n %s"
    no_prefix = "ExpressionModel does not support `prefix` argument"

    def __init__(self, expr, independent_vars=None, init_script=None,
                 nan_policy='raise', **kws):
        """Generate a model from user-supplied expression.

        Parameters
        ----------
        expr : str
            Mathematical expression for model.
        independent_vars : list of str or None, optional
            Variable names to use as independent variables.
        init_script : str or None, optional
            Initial script to run in asteval interpreter.
        nan_policy : str, optional
            How to handle NaN and missing values in data. Must be one of:
            'raise' (default), 'propagate', or 'omit'. See Notes below.
        missing : str, optional
            Synonym for 'nan_policy' for backward compatibility.
        **kws : optional
            Keyword arguments to pass to :class:`Model`.

        Notes
        -----
        1. each instance of ExpressionModel will create and using its own
           version of an asteval interpreter.

        2. prefix is **not supported** for ExpressionModel.

        3. nan_policy sets what to do when a NaN or missing value is seen in
        the data. Should be one of:

            - 'raise' : Raise a ValueError (default)
            - 'propagate' : do nothing
            - 'omit' : (was 'drop') drop missing data

        4. The `missing` argument is deprecated in lmfit 0.9.8 and will be
        removed in a later version. Use `nan_policy` instead, as it is
        consistent with the Minimizer class.

        """
        # create ast evaluator, load custom functions
        self.asteval = Interpreter(max_time=3600.0)
        for name in lineshapes.functions:
            self.asteval.symtable[name] = getattr(lineshapes, name, None)
        if init_script is not None:
            self.asteval.eval(init_script)

        # save expr as text, parse to ast, save for later use
        self.expr = expr.strip()
        self.astcode = self.asteval.parse(self.expr)

        # find all symbol names found in expression
        sym_names = get_ast_names(self.astcode)

        if independent_vars is None and 'x' in sym_names:
            independent_vars = ['x']
        if independent_vars is None:
            raise ValueError(self.idvar_missing % (self.expr))

        # determine which named symbols are parameter names,
        # try to find all independent variables
        idvar_found = [False]*len(independent_vars)
        param_names = []
        for name in sym_names:
            if name in independent_vars:
                idvar_found[independent_vars.index(name)] = True
            elif name not in param_names and name not in self.asteval.symtable:
                param_names.append(name)

        # make sure we have all independent parameters
        if not all(idvar_found):
            lost = []
            for ix, found in enumerate(idvar_found):
                if not found:
                    lost.append(independent_vars[ix])
            lost = ', '.join(lost)
            raise ValueError(self.idvar_notfound % (lost, self.expr))

        kws['independent_vars'] = independent_vars
        if 'prefix' in kws:
            raise Warning(self.no_prefix)

        def _eval(**kwargs):
            for name, val in kwargs.items():
                self.asteval.symtable[name] = val
            self.asteval.start_time = time.time()
            return self.asteval.run(self.astcode)

        kws["nan_policy"] = nan_policy

        super(ExpressionModel, self).__init__(_eval, **kws)

        # set param names here, and other things normally
        # set in _parse_params(), which will be short-circuited.
        self.independent_vars = independent_vars
        self._func_allargs = independent_vars + param_names
        self._param_names = param_names
        self._func_haskeywords = True
        self.def_vals = {}

    def __repr__(self):
        """TODO: docstring in magic method."""
        return "<lmfit.ExpressionModel('%s')>" % (self.expr)

    def _parse_params(self):
        """Over-write ExpressionModel._parse_params with `pass`.

        This prevents normal parsing of function for parameter names.

        """
        pass
back to top