From 06d36f59fb115a0d9ef41d07685e7f73f460880c Mon Sep 17 00:00:00 2001 From: Matthew Newville Date: Sun, 6 Oct 2024 21:37:22 -0500 Subject: [PATCH 1/4] for step lineshape, allow sigma<0 to mean step down from 1 to 0 --- lmfit/lineshapes.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/lmfit/lineshapes.py b/lmfit/lineshapes.py index a9eb6012..28caf7ee 100644 --- a/lmfit/lineshapes.py +++ b/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 -------- From e55c877eae76e7970affac424e443263f0c2df97 Mon Sep 17 00:00:00 2001 From: Matthew Newville Date: Sun, 6 Oct 2024 21:41:20 -0500 Subject: [PATCH 2/4] allow StepModel to have a negative values for sigma --- lmfit/models.py | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/lmfit/models.py b/lmfit/models.py index c4bd7073..ecf20723 100644 --- a/lmfit/models.py +++ b/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 \gt 0` gives a rising step, while :math:`\sigma \lt 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 \gt 0` is enforced, giving a + rising initial step, and :math:`\sigma_2 \gt 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 From b257fca0f2584c2143cf96f42470fea05e65eb06 Mon Sep 17 00:00:00 2001 From: Matthew Newville Date: Sun, 6 Oct 2024 21:43:04 -0500 Subject: [PATCH 3/4] add test of step model with sigma<0, add another test of rectangle model --- tests/test_stepmodel.py | 51 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 50 insertions(+), 1 deletion(-) diff --git a/tests/test_stepmodel.py b/tests/test_stepmodel.py index 47f6dd8d..034217ca 100644 --- a/tests/test_stepmodel.py +++ b/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 From 5debf3dd5e712cb4041435fcf185e4c8c2bd43de Mon Sep 17 00:00:00 2001 From: Matthew Newville Date: Sun, 6 Oct 2024 22:33:54 -0500 Subject: [PATCH 4/4] use >< in place of \gt, \lt --- lmfit/models.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lmfit/models.py b/lmfit/models.py index ecf20723..369fc8b7 100644 --- a/lmfit/models.py +++ b/lmfit/models.py @@ -1541,7 +1541,7 @@ class StepModel(Model): where :math:`\alpha = (x - \mu)/{\sigma}`. - Note that :math:`\sigma \gt 0` gives a rising step, while :math:`\sigma \lt 0` gives + Note that :math:`\sigma > 0` gives a rising step, while :math:`\sigma < 0` gives a falling step. """ @@ -1607,8 +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 \gt 0` is enforced, giving a - rising initial step, and :math:`\sigma_2 \gt 0` gives a falling final step. + 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')