Inconsistencies with non-parameter model function arguments

[schtandard]

Description

I have a model function that contains a keyword argument that is not part of the fit parameters (but not a compulsory independent variable either). My expectation was that I can just pass it as a keyword argument to the model's fit() method when I need to set it. Here's an example with s being the argument in question.

import lmfit
import numpy as np
from matplotlib import pyplot as plt

# Create a model.
def foo(x, a=3, b=5, s=1):
    return a + s * b * x**2
mdl = lmfit.Model(foo, param_names=['a', 'b'])
mdl.set_param_hint('a', value=3)
mdl.set_param_hint('b', value=5)

# Synthesize some noisy data.
x = np.linspace(-5, 5, 101)
s = 1 + 0.7 * np.sin(x)
y = foo(x, 2.3, 1.7, s) + np.random.normal(0, 1, x.shape)

# Fit the data and report the result
fitres = mdl.fit(y, x=x, s=s)
print(fitres.fit_report())
fitres.plot()
plt.show()

# Manually check the result.
plt.plot(x, y, '.', label='data')
plt.plot(x, fitres.eval(x=x), label='fit result eval()')
plt.plot(x, foo(x, a=fitres.params['a'], b=fitres.params['b'], s=s), label='manual check')
plt.legend()
plt.show()

Passing s to fit() like this results a warning.

C:\Users\wilde\miniconda3\lib\site-packages\lmfit\model.py:1053: UserWarning: The keyword argument s does not match any arguments of the model function. It will be ignored.

However, the fit still seems to perform as expected, i.e. the argument is passed on to the fitting function, the resulting parameters are correct. Here's the fit report (a = 2.3 and b = 1.7 were used during data synthesis).

[[Model]]
    Model(foo)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 7
    # data points      = 101
    # variables        = 2
    chi-square         = 109.127855
    reduced chi-square = 1.10230157
    Akaike info crit   = 11.8173660
    Bayesian info crit = 17.0476071
    R-squared          = 0.99636556
[[Variables]]
    a:  2.13258073 +/- 0.13627980 (6.39%) (init = 3)
    b:  1.69611547 +/- 0.01029550 (0.61%) (init = 5)
[[Correlations]] (unreported correlations are < 0.100)
    C(a, b) = -0.6421

Adding to the confusion, the fit result's plot() function draws the correct residuals and data but an incorrect best fit, i.e. it really does ignore s here.

The eval() method on the other hand seems to respect s.

Version information

Python: 3.9.18 | packaged by conda-forge | (main, Aug 30 2023, 03:40:31) [MSC v.1929 64 bit (AMD64)]

lmfit: 1.2.2, scipy: 1.11.3, numpy: 1.26.0,asteval: 0.9.31, uncertainties: 3.1.7

[reneeotten]

The plotting feature is a very rudimentary convenience function. You presumably need to pass the independent variable to the plot function as well. I can take a closer look at this, but it doesn’t necessarily appear as a real GitHub Issue. If you read the instructions you would have realized that, and at least provided all the information requested in the template.

[schtandard]

I did read the instructions and I believe that I did provide all the requested information (though I did take the liberty to remove some section headings as they would have forced a more confusing structure). But no matter.

Though the inconsistent behavior of the plot() function is incorrect and confusing and I would have deemed it to warrant a Github issue on its own (it is an error in the lmfit code, after all; feel free to educate me on why it does not warrant an issue), the more serious part of the issue, in my opinion, is the warning, as it is exclusively states wrong information:

• s does match an argument of the model function.
• It is not ignored.

[newville]

@schtandard

I did read the instructions and I believe that I did provide all the requested information

You did not provide all the requested information.

(though I did take the liberty to remove some section headings as they would have forced a more confusing structure). But no matter.

It may be no matter to you. It is a matter to us.

Again, you probably really want to specify that s is an independent variable as with

mdl = lmfit.Model(foo, independent_vars=['x', 's'])

[schtandard]

I'm unsure how to interpret the lack of response here. Are you dealing with this in the background, do you not want to engage with this issue or is time just in too short supply at the moment? Would it help you if I opened another issue, following the template to the letter?

[reneeotten]

okay, so that seems like a very specific, unconventional way of doing things. I doubt we will go out of our way to support that.

[newville]

@schtandard I think what you described for your code is exactly that x and s should be independent variables. That was approximately the first thing we said, and you sort of rejected that idea immediately. You have repeatedly said that all independent variables must always be specified but that is simply untrue. You sort of keep saying there is a bug or Issue. I don't really understand what the alleged bug is.

If you had defined your function with a signature of s=None and then put in your code if s is None: s = 1.0 (I probably would have used scale=None in the signature and then if scale is not None: result *= scale in the code), your life probably would have been a lot easier: your signature would make s an independent variable by default, and allow it to be an array sometimes (as your example), and have sensible default behavior.

But, if the goal is to make your code complicated, that is certainly possible. ;)

[schtandard]

You have repeatedly said that all independent variables must always be specified but that is simply untrue.

Ok, then I am missing something. If I set s to be an independent variable and then not specify it for fitting

import lmfit
import numpy as np

# Create a model.
def foo(x, a=3, b=5, s=1):
    return a + s * b * x**2
mdl = lmfit.Model(foo, independent_vars=['x', 's'])

# Synthesize some noisy data.
x = np.linspace(-5, 5, 101)
y = foo(x, 2.3, 1.7) + np.random.normal(0, 1, x.shape)

# Fit the data.
fitres = mdl.fit(y, x=x)

I get an error

Traceback (most recent call last):

  File ~\miniconda3\Lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File c:\users\wilde\desktop\temptest\lmfit\nonparamargs3.py:15
    fitres = mdl.fit(y, x=x)

  File ~\miniconda3\Lib\site-packages\lmfit\model.py:1086 in fit
    if not np.isscalar(kwargs[var]):

KeyError: 's'

How do I use the model for fitting without specifying s?

If you had defined your function with a signature of s=None and then put in your code if s is None: s = 1.0 (I probably would have used scale=None in the signature and then if scale is not None: result *= scale in the code), your life probably would have been a lot easier: your signature would make s an independent variable by default, and allow it to be an array sometimes (as your example), and have sensible default behavior.

But, if the goal is to make your code complicated, that is certainly possible. ;)

No, unless I am really missing quite a lot. Using s=None in the function signature keeps s from being a fitting parameter but it does not make it an independent variable:

import lmfit

# Create a model.
def foo(x, a=3, b=5, s=None):
    if s is None:
        s = 1
    return a + s * b * x**2
mdl = lmfit.Model(foo)

print(mdl.independent_vars)

yields just ['x']. In other words, it leads to (as far as I can tell) exactly the same model as

# Create a model.
def foo(x, a=3, b=5, s=1):
    return a + s * b * x**2
mdl = lmfit.Model(foo, param_names=['a', 'b'])
mdl.set_param_hint('a', value=3)
mdl.set_param_hint('b', value=5)

I, too, am a fan of clarity, which is why I prefer

def foo(x, a=3, b=5, s=1):

to

def foo(x, a=3, b=5, s=None):
    if s is None:
        s = 1

(or any equivalent phrasing).

Of course, all of this has nothing to do with the (alleged) bug(s), rather there are two differently clumsy ways of achieving what I want which might as well be the simpler and clearer

# Create a model.
def foo(x, a=3, b=5, s=1):
    return a + s * b * x**2
mdl = lmfit.Model(foo, param_names=['a', 'b'])

I wanted to keep that part of the discussion in the corresponding feature request (#919), but this seems to have failed.

You sort of keep saying there is a bug or Issue. I don't really understand what the alleged bug is.

Alright, I will try to clear this up once more. As I see it, there are two things constituting a bug here, though they are so closely related that, without looking into the code, it is hard to say if they are really separate or if there is just one underlying bug.

1. Executing the this code

   import lmfit
   import numpy as np
   from matplotlib import pyplot as plt
   
   # Create a model.
   def foo(x, a=3, b=5, s=1):
       return a + s * b * x**2
   mdl = lmfit.Model(foo, param_names=['a', 'b'])
   mdl.set_param_hint('a', value=3)
   mdl.set_param_hint('b', value=5)
   
   # Synthesize some noisy data.
   x = np.linspace(-5, 5, 101)
   s = 1 + 0.7 * np.sin(x)
   y = foo(x, 2.3, 1.7, s) + np.random.normal(0, 1, x.shape)
   
   # Fit the data and report the result
   fitres = mdl.fit(y, x=x, s=s)

   leads to the warning

   C:\Users\wilde\miniconda3\lib\site-packages\lmfit\model.py:1053: UserWarning: The keyword argument s does not match any arguments of the model function. It will be ignored.
   

   • The first part The keyword argument s does not match any arguments of the model function. is obviously wrong.
   • Inspection of the fit result

     print(fitres.fit_report())

     and manually plotting the fit result (either from the parameters or using fitres.eval()) reveals that the second part It will be ignored. is also untrue.

   To me, this incorrect warning is the first clear bug (possibly pointing to more serious, yet undetected problems with handling s).

2. Now, instead of plotting the fit result manually, one might employ the builtin method fitres.plot(), as did I when encountering the warning above. This draws a plot of the data with a curve labeled "best fit". This curve does not in fact show the best fit (as manual plotting reveals). Furthermore, a second plot is drawn showing the residual for each data point. The two plots are obviously inconsistent, i.e. they show contradictory data.

   To me, this is the second clear bug (possibly an effect of the same underlying problem as the first, possibly separate).

Now, to the second point @reneeotten said "The plotting feature is a very rudimentary convenience function." which is certainly true but should not affect if it is a bug or not. To the first point I have not seen a response other than "just don't do that then". Again, I won't tell you how to run your open source project and if this is something not worth fixing to you, so be it. But saying "bugs in parts of the code that are not important to us are not bugs" seems questionable. Maybe this is not what you are trying to say but this is the gist I am getting.

[newville]

@schtandard Well, you can certainly specify param_names if you want... but that doesn't allow you to specify default values. That seems like a problem for you.

If you have your signature for s be non-numerical to turn it into an independent variable, then your code will work much more simply. Maybe something like this:

import numpy as np
from lmfit import Model

def func(x, a=3, b=5, scale=None):
    output = a +x*b
    if scale is not None:  output *= scale
    return output

x = np.linspace(-10, 10, 51)
y = func(x, a=3.5, b=1.0) + np.random.normal(size=len(x), scale=0.5)

# variation 1 -- just do not mention `scale`:
mod = Model(func)
print(mod.eval(x=2))
print(mod.eval(x=2, scale=10))
print(mod.eval(x=2, a=0, scale=10))

params = mod.make_params(a=2.5, b=2.0)
for par in params.values():
    print(par)

out = mod.fit(y, params, x=x)
print(out.fit_report())

print("## now with scale")
out = mod.fit(y, params, x=x, scale=2.0)
print(out.fit_report())


# variation 1 -- state that `scale` is an independent parameter
mod = Model(func, independet_vars=['x', 'scale'])
print(mod.eval(x=2))
print(mod.eval(x=2, scale=10))
print(mod.eval(x=2, a=0, scale=10))

params = mod.make_params(a=2.5, b=2.0)
for par in params.values():
    print(par)

out = mod.fit(y, params, x=x, scale=None)
print(out.fit_report())

print("## now with scale")
out = mod.fit(y, params, x=x, scale=2.0)
print(out.fit_report())

Is variation 1 not acceptable?
Is the complaint that with variation 2, you cannot do out = mod.fit(y, params, x=x)?

It sort of seems like you are aiming to make your life hard and trying to find an issue.

[schtandard]

If you have your signature for s be non-numerical to turn it into an independent variable,

Again, this is not what happens. Unless you call every argument that is not a fitting parameter an "independent variable". That would be misleading, though.

Is variation 1 not acceptable?

Variant 1 (the first one, I assume) is certainly acceptable, though I feel like leaving the function signature transparent and being more verbose when deriving the model is preferable to obscuring the function signature by hiding the default value in the function body (in this example it may be natural to assume that a scale=None is equivalent to scale=1, but this not always the case; explicit is better than implicit, also in the function signature) just for the brevity at model creation.

What I would prefer, of course, is the best of both worlds: A clear function signature and a brief model derivation. All the information is there, after all. That's what #919 is about.

As I said, that is completely orthogonal to this bug report. In fact, your code produces the same erroneous warning as mine. So, even with that approach, I would have come here to file the same bug report, though the code example would have looked slightly different.

Is the complaint that with variation 2, you cannot do out = mod.fit(y, params, x=x)?

If by "complaint" you mean the reason why I dismissed your suggestion of using s as an independent variable, then yes. This is what I expressed by saying "being an independent variable it has to be specified every time", though I wouldn't call that a complaint.

But let us cease this discussion (i.e. the one pertaining to #919), I feel that we are turning in circles. Since you said that my desired change would not be rejected (reliant on being well executed, of course) I will try and prepare a PR.

It sort of seems like you are aiming to make your life hard and trying to find an issue.

I assure you that my life is already harder than I would prefer. As to finding an issue, I believe that I have already done so and outlined this in detail, though that part of my comments has not been addressed.