Issue with eval_uncertainty when param vary=False

[brycehenson]

I think i have discovered a bug in eval_uncertainty when a parameter vary=Flase is set.
It seems that eval_uncertainty cant handle a non varied parameter. Is this expected behaviour ?

def linear_dep(x, var1):
    return x * var1


def constant_offset(
    x,
    offset,
):
    return offset + x * 0.0


def test_fit_lm_combining_models():
    predictor_var = np.linspace(0, 1, 100)
    response_var = linear_dep(
        x=predictor_var,
        var1=100,
    )

    # create a model of the linear combined with offset
    sat = Model(linear_dep)
    sat.set_param_hint("var1", min=0.0)
    offset = Model(constant_offset, prefix="")
    model = sat + offset
    model.set_param_hint("var1", value=max(response_var) * 0.8)
    model.set_param_hint("offset", value=0.0)
    params = model.make_params()

    # setting this vary=false causes the bug
    params["offset"].set(value=0.0, vary=False)
    fit = model.fit(data=response_var, x=predictor_var, params=params)
    x_samp = np.linspace(np.min(predictor_var), np.max(predictor_var), 200)

    y_fit_val = fit.eval(x=x_samp)

    y_ci = fit.eval_uncertainty(x=x_samp)

    

pytest output

packages/measurements/tests/test_fitlm.py::test_fit_lm_combining_models failed: def test_fit_lm_combining_models():
    
        predictor_var = np.linspace(0, 1, 100)
        response_var = linear_dep(
            x=predictor_var,
            var1=100,
        )
    
        # create a model of the linear combined with offset
        sat = Model(linear_dep)
        sat.set_param_hint("var1", min=0.0)
        offset = Model(constant_offset, prefix="")
        model = sat + offset
    
        model.set_param_hint("var1", value=max(response_var) * 0.8)
        model.set_param_hint("offset", value=0.0)
    
        params = model.make_params()
    
        # setting this vary=false causes the bug
        params["offset"].set(value=0.0, vary=False)
    
        fit = model.fit(data=response_var, x=predictor_var, params=params)
        x_samp = np.linspace(np.min(predictor_var), np.max(predictor_var), 200)
    
        y_fit_val = fit.eval(x=x_samp)
    
>       y_ci = fit.eval_uncertainty(x=x_samp)

packages/measurements/tests/test_fitlm.py:75: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

self = <lmfit.model.ModelResult object at 0x7faef812ea70>
params = Parameters([('var1', <Parameter 'var1', value=100.0 +/- 0, bounds=[0.0:inf]>), ('offset', <Parameter 'offset', value=0.0 (fixed), bounds=[-inf:inf]>)])
sigma = 1
kwargs = {'x': array([0.        , 0.00502513, 0.01005025, 0.01507538, 0.0201005 ,
       0.02512563, 0.03015075, 0.03517588, 0....7, 0.95979899, 0.96482412, 0.96984925, 0.97487437,
       0.9798995 , 0.98492462, 0.98994975, 0.99497487, 1.        ])}
userkws = {'x': array([0.        , 0.00502513, 0.01005025, 0.01507538, 0.0201005 ,
       0.02512563, 0.03015075, 0.03517588, 0....7, 0.95979899, 0.96482412, 0.96984925, 0.97487437,
       0.9798995 , 0.98492462, 0.98994975, 0.99497487, 1.        ])}
nvarys = 1
feval = array([  0.        ,   0.50251256,   1.00502513,   1.50753769,
         2.01005025,   2.51256281,   3.01507538,   3.51... 96.48241206,  96.98492462,  97.48743719,  97.98994975,
        98.49246231,  98.99497487,  99.49748744, 100.        ])
ndata = 200, covar = array([[0.]])
fjac = {'0': array([[nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
        nan, nan, nan, nan, nan, nan, n... nan, nan,
        nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
        nan, nan, nan, nan, nan]])}
df2 = {'0': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0....., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])}

    def eval_uncertainty(self, params=None, sigma=1, **kwargs):
        """Evaluate the uncertainty of the *model function*.
    
        This can be used to give confidence bands for the model from the
        uncertainties in the best-fit parameters.
    
        Parameters
        ----------
        params : Parameters, optional
            Parameters, defaults to ModelResult.params.
        sigma : float, optional
            Confidence level, i.e. how many sigma (default is 1).
        **kwargs : optional
            Values of options, independent variables, etcetera.
    
        Returns
        -------
        numpy.ndarray
            Uncertainty at each value of the model.
    
        Notes
        -----
        1. This is based on the excellent and clear example from
           https://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html#confidence-and-prediction-intervals,
           which references the original work of:
           J. Wolberg, Data Analysis Using the Method of Least Squares, 2006, Springer
        2. The value of sigma is number of `sigma` values, and is converted
           to a probability. Values of 1, 2, or 3 give probabilities of
           0.6827, 0.9545, and 0.9973, respectively. If the sigma value is
           < 1, it is interpreted as the probability itself. That is,
           ``sigma=1`` and ``sigma=0.6827`` will give the same results,
           within precision errors.
        3. Also sets attributes of `dely` for the uncertainty of the model
           (which will be the same as the array returned by this method) and
           `dely_comps`, a dictionary of `dely` for each component.
    
    
        Examples
        --------
    
        >>> out = model.fit(data, params, x=x)
        >>> dely = out.eval_uncertainty(x=x)
        >>> plt.plot(x, data)
        >>> plt.plot(x, out.best_fit)
        >>> plt.fill_between(x, out.best_fit-dely,
        ...                  out.best_fit+dely, color='#888888')
    
        """
        userkws = self.userkws.copy()
        userkws.update(kwargs)
        if params is None:
            params = self.params
    
        nvarys = self.nvarys
        # ensure fjac and df2 are correct size if independent var updated by kwargs
        feval = self.model.eval(params, **userkws)
        ndata = len(feval.view('float64'))        # allows feval to be complex
        covar = self.covar
        if any(p.stderr is None for p in params.values()):
            return np.zeros(ndata)
    
        # '0' would be an invalid prefix, here signifying 'Full'
        fjac = {'0': np.zeros((nvarys, ndata), dtype='float64')}
        df2 = {'0': np.zeros(ndata, dtype='float64')}
    
        for comp in self.components:
            label = comp.prefix if len(comp.prefix) > 1 else comp._name
            fjac[label] = np.zeros((nvarys, ndata), dtype='float64')
            df2[label] = np.zeros(ndata, dtype='float64')
    
        # find derivative by hand!
        pars = params.copy()
        for i in range(nvarys):
            pname = self.var_names[i]
            val0 = pars[pname].value
            dval = pars[pname].stderr/3.0
            pars[pname].value = val0 + dval
            res1 = {'0': self.model.eval(pars, **userkws)}
            res1.update(self.model.eval_components(params=pars, **userkws))
    
            pars[pname].value = val0 - dval
            res2 = {'0': self.model.eval(pars, **userkws)}
            res2.update(self.model.eval_components(params=pars, **userkws))
    
            pars[pname].value = val0
            for key in fjac:
>               fjac[key][i] = (res1[key].view('float64')
                                - res2[key].view('float64')) / (2*dval)
E               AttributeError: 'float' object has no attribute 'view'

/usr/local/lib/python3.10/site-packages/lmfit/model.py:1619: AttributeError

my system information

Python 3.10.12 (main, Jul  4 2023, 05:58:45) [GCC 12.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import sys, lmfit, numpy, scipy, asteval, uncertainties
it.__version__}, scipy: {scipy.__version__}, numpy: {numpy.__version__},"
          f"asteval: {asteval.__version__}, uncertainties: {uncertainties.__version__}")
>>> print(f"Python: {sys.version}\n\nlmfit: {lmfit.__version__}, scipy: {scipy.__version__}, numpy: {numpy.__version__},"
...       f"asteval: {asteval.__version__}, uncertainties: {uncertainties.__version__}")
Python: 3.10.12 (main, Jul  4 2023, 05:58:45) [GCC 12.2.0]
lmfit: 1.2.2, scipy: 1.10.1, numpy: 1.24.4,asteval: 0.9.31, uncertainties: 3.1.7

My next steps

• i will simplify this MWE/test
• explore the multi parameter case
• create an issue if others aggree that this is not expected behaviour
• create a branch with a test of this behaviour
• try and fix

[brycehenson]

I have checked that in lmfit = "1.2.1" the above example runs without issue.

[brycehenson]

also check that same problem with scipy 1.11.1

Python 3.10.12 (main, Jul  4 2023, 05:58:45) [GCC 12.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import sys, lmfit, numpy, scipy, asteval, uncertainties
it.__version__}, scipy: {scipy.__version__}, numpy: {numpy.__version__},"
          f"asteval: {asteval.__version__}, uncertainties: {uncertainties.__version__}")
>>> print(f"Python: {sys.version}\n\nlmfit: {lmfit.__version__}, scipy: {scipy.__version__}, numpy: {numpy.__version__},"
...       f"asteval: {asteval.__version__}, uncertainties: {uncertainties.__version__}")
Python: 3.10.12 (main, Jul  4 2023, 05:58:45) [GCC 12.2.0]

lmfit: 1.2.2, scipy: 1.11.1, numpy: 1.24.4,asteval: 0.9.31, uncertainties: 3.1.7
>>> 

[newville]

@brycehenson Hm, maybe. You did not show the full traceback. Where is the error at?

Your example does not run for me. There are no import statements at all. What is npt supposed to be?

Definitely "No" on the async def wrapper. Yuck. FWIW, I despise the whole black nonsense and overly-typed code. Honestly, it makes me want to walk away from supporting users.

A model function that returns some array-like thing that is not numpy ndarray is definitely anti-recommended. A convenience function like eval_uncertainty() may very well complain if you fail to use ndarrays uniformly.

[brycehenson]

@newville thanks for the feedback. I have simplified the test and added the missing imports.

You are indeed correct that it was the mixing of arrays and floats from the two model functions.

Changing the above code to

def constant_offset(x,offset,
) :
    return offset+x*0.0

Fixes the issue. It seems that this is quite the corner case then requiring vary = false and the models to have mismatched types.