Steps with neg sigma

lmfit/lineshapes.py

@@ -1,7 +1,7 @@
 """Basic model lineshapes and distribution functions."""
 
 from numpy import (arctan, copysign, cos, exp, isclose, isnan, log, log1p,
-                   maximum, minimum, pi, real, sin, sqrt, where)
+                   maximum, minimum, pi, real, sign, sin, sqrt, where)
 from scipy.special import betaln as betalnfcn
 from scipy.special import erf, erfc
 from scipy.special import gamma as gamfcn
@@ -415,8 +415,8 @@ def thermal_distribution(x, amplitude=1.0, center=0.0, kt=1.0, form='bose'):
 def step(x, amplitude=1.0, center=0.0, sigma=1.0, form='linear'):
     """Return a step function.
 
-    Starts at 0.0, ends at `amplitude`, with half-max at `center`, and
-    rising with `form`:
+    Starts at 0.0, ends at `sign(sigma)*amplitude`, has a half-max at
+    `center`, rsing or falling with `form`:
 
     - `'linear'` (default) = amplitude * min(1, max(0, arg + 0.5))
     - `'atan'`, `'arctan'` = amplitude * (0.5 + atan(arg)/pi)
@@ -425,8 +425,10 @@ def step(x, amplitude=1.0, center=0.0, sigma=1.0, form='linear'):
 
     where ``arg = (x - center)/sigma``.
 
+    Note that ``sigma > 0`` gives a rising step, while ``sigma < 0`` gives
+    a falling step.
     """
-    out = (x - center)/max(tiny, sigma)
+    out = sign(sigma)*(x - center)/max(tiny*tiny, abs(sigma))
 
     if form == 'erf':
         out = 0.5*(1 + erf(out))
@@ -455,8 +457,11 @@ def rectangle(x, amplitude=1.0, center1=0.0, sigma1=1.0,
     - `'erf'`              = amplitude*(erf(arg1) + erf(arg2))/2.
     - `'logisitic'`        = amplitude*[1 - 1/(1 + exp(arg1)) - 1/(1+exp(arg2))]
 
-    where ``arg1 = (x - center1)/sigma1`` and
-    ``arg2 = -(x - center2)/sigma2``.
+    where ``arg1 = (x - center1)/sigma1`` and ``arg2 = -(x - center2)/sigma2``.
+
+    Note that, unlike `step`,  ``sigma1 > 0`` and ``sigma2 > 0``, so that a
+    rectangle can support a step up followed by a step down.
+    Use a constant offset and adjust amplitude if that is what you need.
 
     See Also
     --------

lmfit/models.py

@@ -1525,8 +1525,9 @@ class StepModel(Model):
       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:`\tt{sign}(\sigma)A` rising or falling to :math:`A/2` at :math:`\mu`,
+    with :math:`\sigma` setting the characteristic width and the sign up the step.
+    The functional forms are defined as:
 
     .. math::
         :nowrap:
@@ -1540,6 +1541,8 @@ class StepModel(Model):
 
     where :math:`\alpha = (x - \mu)/{\sigma}`.
 
+    Note that :math:`\sigma > 0` gives a rising step, while :math:`\sigma < 0` gives
+    a falling step.
     """
 
     valid_forms = ('linear', 'atan', 'arctan', 'erf', 'logistic')
@@ -1556,7 +1559,11 @@ def guess(self, data, x, **kwargs):
         xmin, xmax = min(x), max(x)
         pars = self.make_params(amplitude=(ymax-ymin),
                                 center=(xmax+xmin)/2.0)
-        pars[f'{self.prefix}sigma'].set(value=(xmax-xmin)/7.0, min=0.0)
+        n = len(data)
+        sigma = 0.1*(xmax - xmin)
+        if data[:n//5].mean() > data[-n//5:].mean():
+            sigma = -sigma
+        pars[f'{self.prefix}sigma'].set(value=sigma)
         return update_param_vals(pars, self.prefix, **kwargs)
 
     __init__.__doc__ = COMMON_INIT_DOC
@@ -1600,6 +1607,8 @@ class RectangleModel(Model):
     where :math:`\alpha_1 = (x - \mu_1)/{\sigma_1}` and
     :math:`\alpha_2 = -(x - \mu_2)/{\sigma_2}`.
 
+    Note that, unlike a StepModel, :math:`\sigma_1 > 0` is enforced, giving a
+    rising initial step, and  :math:`\sigma_2 > 0` gives a falling final step.
     """
 
     valid_forms = ('linear', 'atan', 'arctan', 'erf', 'logistic')
@@ -1624,9 +1633,10 @@ def guess(self, data, x, **kwargs):
         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[f'{self.prefix}sigma1'].set(value=(xmax-xmin)/7.0, min=0.0)
-        pars[f'{self.prefix}sigma2'].set(value=(xmax-xmin)/7.0, min=0.0)
+                                center2=3*(xmax+xmin)/4.0,
+                                sigma1={'value': (xmax-xmin)/10.0, 'min': 0},
+                                sigma2={'value': (xmax-xmin)/10.0, 'min': 0})
+
         return update_param_vals(pars, self.prefix, **kwargs)
 
     __init__.__doc__ = COMMON_INIT_DOC

tests/test_stepmodel.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-from lmfit.models import ConstantModel, StepModel
+from lmfit.models import ConstantModel, RectangleModel, StepModel
 
 
 def get_data():
@@ -54,3 +54,52 @@ def test_stepmodel_erf():
     assert out.params['amplitude'].value > 50
     assert out.params['sigma'].value > 0.2
     assert out.params['sigma'].value < 1.5
+
+
+def test_stepmodel_stepdown():
+    x = np.linspace(0, 50, 201)
+    y = np.ones_like(x)
+    y[129:] = 0.0
+    y[109:129] = 1.0 - np.arange(20)/20
+    y = y + 5e-2*np.random.randn(len(x))
+    stepmod = StepModel(form='linear')
+    pars = stepmod.guess(y, x)
+
+    out = stepmod.fit(y, pars, x=x)
+
+    assert out.nfev > 10
+    assert out.nvarys == 3
+    assert out.chisqr > 0.2
+    assert out.chisqr < 5.0
+    assert out.params['center'].value > 28
+    assert out.params['center'].value < 32
+    assert out.params['amplitude'].value > 0.5
+    assert out.params['sigma'].value < -2.0
+    assert out.params['sigma'].value > -8.0
+
+
+def test_rectangle():
+    x = np.linspace(0, 50, 201)
+    y = np.ones_like(x) * 2.5
+    y[:33] = 0.0
+    y[162:] = 0.0
+    y[33:50] = 2.5*np.arange(50-33)/(50-33)
+    y[155:162] = 2.5 * (1 - np.arange(162-155)/(162-155))
+    y = y + 5e-2*np.random.randn(len(x))
+    stepmod = RectangleModel(form='linear')
+    pars = stepmod.guess(y, x)
+
+    out = stepmod.fit(y, pars, x=x)
+
+    assert out.nfev > 10
+    assert out.nvarys == 5
+    assert out.chisqr > 0.2
+    assert out.chisqr < 5.0
+    assert out.params['center1'].value > 8
+    assert out.params['center1'].value < 14
+    assert out.params['amplitude'].value > 2.0
+    assert out.params['amplitude'].value < 10.0
+    assert out.params['sigma1'].value > 1.0
+    assert out.params['sigma1'].value < 5.0
+    assert out.params['sigma2'].value > 0.3
+    assert out.params['sigma2'].value < 2.5