RectangleModule implements 'linear' incorrectly?

[arunpersaud]

First Time Issue Code

Yes, I read the instructions and I am sure this is a GitHub Issue.

Description

The RectangleModule linear case seems have a bug in its implementation.

The docs say that for the RectangleModule the fitting function returns A/2 at u1 and u2.

The docs give the following function:

and the code here defines the function slightly different:

        arg1[arg1 < 0] = 0.0
        arg1[arg1 > 1] = 1.0
        arg2[arg2 > 0] = 0.0
        arg2[arg2 < -1] = -1.0
        out = arg1 + arg2

It probably should be:

    arg1 = (x-u1)/s1
    arg2 = -(x-u2)/s2
    # should go from mu-sigma to mu+sigma so -1 to +1
    arg1[arg1 < -1] = -1.0
    arg1[arg1 > 1] = 1.0
    arg2[arg2 > 1] = 1.0
    arg2[arg2 < -1] = -1.0
    out = arg1 + arg2
    out = out/2  # needs to be scaled by x2, since we now go from -1 to 1
    return out

A Minimal, Complete, and Verifiable example

Here is a short code snippet that plots all 3 cases

from matplotlib import pyplot as plt
import numpy as np
from scipy.special import erf

u1 = 0
u2 = 10
s1 = 2
s2 = 2
A = 1

x = np.linspace(-10, 20, 101)

a1 = (x-u1)/s1
a2 = -(x-u2)/s2

def expected(x):
    arg1 = (x-u1)/s1
    arg2 = -(x-u2)/s2
    # should go from -sigma to +sigma so -1 to +1
    arg1[arg1 < -1] = -1
    arg1[arg1 > 1] = 1.0
    arg2[arg2 > 1] = 1
    arg2[arg2 < -1] = -1.0
    out = arg1 + arg2
    out = out/2  # needs to be scaled by x2, since we now go from -1 to 1
    return A*out

def docs(x):
    """https://lmfit.github.io/lmfit-py/builtin_models.html#rectanglemodel"""
    return A*(np.minimum(1, np.maximum(0, a1)) + np.minimum(-1, np.maximum(0, a2)))
        
def code(x):
    """lineshapes.py: lines 478-482"""
    arg1 = (x-u1)/s1
    arg2 = -(x-u2)/s2
    arg1[arg1 < 0] = 0.0
    arg1[arg1 > 1] = 1.0
    arg2[arg2 > 0] = 0.0
    arg2[arg2 < -1] = -1.0
    out = arg1 + arg2
    return A*out

fig, ax = plt.subplots()

ax.plot(x, expected(x), label="expected")
ax.plot(x, code(x), label="current code")
ax.plot(x, docs(x), label="docs")

ax.scatter([u1, u2],[A/2, A/2], label="points according to docs")
plt.legend()
plt.show()

Which results in the following plot:

Notes

Happy to supply a PR for this if the above is correct. If tests are requried should lineshape.py be tested? if so what should be tested? f(u1,2) = A/2 and integral == 1?

Version information

Code taken from the current (2022-10-23) web page and commit 86a2bdc

[newville]

@arunpersaud Thanks! Yep, both code and docs are just plain wrong. I have no explanation other than stupidity ;).

I might suggest that code and doc would be best served by changing the block at

lmfit-py/lmfit/lineshapes.py

Line 477 in 86a2bdc

elif form == 'linear':

to be closer to the (corrected, of course) documentation:

    elif form == 'linear':
        out = (np.minimum(1, np.maximum(-1, arg1)) +
               np.minimum(1, np.maximum(-1, arg2)))/ 2.0

If you're willing to make a Pull Request (including doc change), that would be great. If not, that's fine too and we can make this change.

Thanks again!

[arunpersaud]

Sure, I'll give it a try...