Plotting of ConstantModel Fails

Description

If I try and fit some data to a ConstantModel and then plot the results, I get a TypeError

A Minimal, Complete, and Verifiable example

import numpy as np
import lmfit

model = lmfit.models.ConstantModel()
result = model.fit(np.random.random((128,)), x=np.arange(128))

result.plot()

If you change out ConstantModel for a GaussianModel, for example, this works fine.

Error message:

Traceback (most recent call last):
  File "C:\Users\nitrocalcite\Miniconda3\envs\base\lib\site-packages\IPython\core\interactiveshell.py", line 3331, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-39-20f99b654d42>", line 7, in <module>
    result.plot()
  File "C:\Users\nitrocalcite\Miniconda3\envs\base\lib\site-packages\lmfit\model.py", line 49, in wrapper
    return function(*args, **kws)
  File "C:\Users\nitrocalcite\Miniconda3\envs\base\lib\site-packages\lmfit\model.py", line 2081, in plot
    show_init=show_init, parse_complex=parse_complex)
  File "C:\Users\nitrocalcite\Miniconda3\envs\base\lib\site-packages\lmfit\model.py", line 49, in wrapper
    return function(*args, **kws)
  File "C:\Users\nitrocalcite\Miniconda3\envs\base\lib\site-packages\lmfit\model.py", line 1863, in plot_fit
    **{independent_var: x_array_dense})),
  File "C:\Users\nitrocalcite\Miniconda3\envs\base\lib\site-packages\lmfit\model.py", line 93, in reducer
    if any(np.iscomplex(array)):
TypeError: 'numpy.bool_' object is not iterable

Version information

Python: 3.7.7 (default, Mar 23 2020, 23:19:08) [MSC v.1916 64 bit (AMD64)] lmfit: 1.0.1, scipy: 1.4.1, numpy: 1.18.1, asteval: 0.9.19, uncertainties: 3.1.4

Link(s)

I think this discussion is related. ConstantModel and ComplexConstantModel both return single floats. This is different than every other model, which returns something with the same shape as the independent variable. Thus, I suspect this is what's happening in the stack trace:

array is just a float, not an array
np.iscomplex returns a bool, instead of an array of bools
any then gets a bool, rather than an iterable

This issue gets patched over when ConstantModel is used in a CompositeModel, because numpy will coerce 5 + np.array([0,1,2]) to an array of shape (3,). But it becomes problematic when ConstantModel is used alone.

As a proposed fix, I'll suggest changing the underlying model function for ConstantModel, which look like this right now:

        def constant(x, c=0.0):
            return c
        super().__init__(constant, **kwargs)

to something like:

        def constant(x, c=0.0):
            return c * np.ones((len(x), ))
        super().__init__(constant, **kwargs)

thus making the return types consistent. I'm happy to submit such a PR & tests if this solution is desired. Thanks for your time!

@nitrocalcite

I think this discussion is related. ConstantModel and ComplexConstantModel both return single floats.

Yep.

This is different than every other model, which returns something with the same shape as the independent variable.

That is not completely correct. All of the built-in models have this feature and are limited to 1 independent variable that is assumed to be 1-D and the first argument in the call signature. But that is not assumed to always be true for user-defined models.

As a proposed fix, I'll suggest changing the underlying model function for ConstantModel, to something like:
   def constant(x, c=0.0):
       return c * np.ones((len(x), ))
   super().__init__(constant, **kwargs)
thus making the return types consistent.

-1 from me. I think you should re-read the discussion you linked to at https://groups.google.com/g/lmfit-py/c/QQYhastApsE/m/ayjKgGTeGgAJ?pli=1

The Independent variables of a Model may be things other than 1-D arrays of floats. They may be complex, they may be pandas dataframes or Series, they may be dictionaries or handles to database connections. There may be multiple independent variables. The return value of the model function does not need to be the same shape as the independent variable even when that independent variable has a meaningful shape.

The model function returns a scalar or a 1-d array of floats - that is what the minimization routines want the objective/model functions to return. They don't consider the independent variable at all.

The plotting routines are a constant source of problems. They may work for 80% or more of the cases, and they are still perilous. From nothing more than a maintenance POV, they were a mistake. Relying on matplotlib (which, to be clear, I use every day) is a mistake. The "parse_complex" argument is incomplete and should do something sensible for a morel that returns a scalar. All of the plotting functions more or less assume that the independent variable is a 1-D array and that the Model will calculate a result with a one-to-one correspondence to that independent variable. That may very well describe the vast majority of cases, but it does not describe them all.

So, we could go with: a) Shrug. Well, sometimes plotting might fail. lmfit is not a plotting library.
b) Try to fix this case, while not breaking too much else along the way and declare victory that "we fixed it". Well, until it breaks again.

I don't really have a strong preference.

@newville OK, I see what you're talking about in the linked discussion better now. I thought your concerns there were specifically about forcing all models into numpy arrays by type coercion in the Model.eval method, not about one particular model doing so:

Yeah, that is a problem. I think that Model.eval() would have to be changed to return np.asarray(func(**self.make_funcargs(params, kwargs))

Correct me if I'm wrong here, but it seems to me that you want to preserve that fact that ConstantModel doesn't assume array-structured inputs, though the downside is that ConstantModel behaves differently with array-structured inputs as compared to other models: say, GaussianModel or PolynomialModel. If this is indeed your concern I'll propose an alternate solution:

        def constant(x, c=0.0):
            return c * np.ones((len(x), )) if isinstance(x, np.ndarray) else c
        super().__init__(constant, **kwargs)

This would retain ConstantModel's compatibility with non-array types, while also bringing its array behavior in-line with other models' array behavior. Let me know if you consider this a more satisfactory solution.

@nitrocalcite

Correct me if I'm wrong here, but it seems to me that you want to preserve that fact that ConstantModel doesn't assume array-structured inputs, though the downside is that ConstantModel behaves differently with array-structured inputs as compared to other models: say, GaussianModel or PolynomialModel.

Well, it acts differently because it is different. It models a constant. I am pretty sure we never say that other models will behave consistently. We could have multi-dimensional models or models that transform the data in some complex way (say Fourier transforming) for example and return values from that calculation. We definitely want users to be able to do such things. And if they want to add a Constant background, they should be able to do that. I really see no reason to change ConstantModel.

The problem you see is not with the way the ConstantModel behaves, but with PLOTTING. If you want to fix that plotting, go ahead. The traceback you posted points out exactly what the problem is.

Fixing plotting problems is a very low priority for me. I would be more easily persuaded to remove the plotting completely than to try to make it work for every situation. Again, it will sort of work for simple 1-D cases. And it will fail other times.

@newville

We could have multi-dimensional models or models that transform the data in some complex way (say Fourier)... We definitely want users to be able to do such things.

I completely agree that you don't want to restrict the user in principle. However, 1) my proposal doesn't affect the Model base class and 2) while this may be a documentation gap, I don't see anything in the models guide, built-in models list, or examples that suggest Model is intended for arbitrary data pre-processing in the manner you suggest here.

I am pretty sure we never say that other models will behave consistently.

Again while this might be a documentation issue, as far as I can tell every piece of documentation in the models guide, built-in models list, and examples index show models behaving consistently with array inputs. Thus while not explicitly stated, I think it's fair to say that users may reasonably have this expectation.

The problem you see is...with PLOTTING.

My direct problem in usage is not plotting. It is further processing. Given an array that I feed into a built-in lmfit model, I do not know how to tell what data structure to expect back for further processing, nor how to interpret it. If this is intended to be done in some other way, please let me know. (For example, perhaps there's a nominal subclass of Model denoting 1D->1D array evals, as (to my knowledge) all the built-in models save Constant currently support)

(When I first had this issue, I looked into how lmfit might handle different data structures return internally by checking plot. Given that all the other models I've tried worked in the same manner, I assumed's ConstantModel's different behavior with arrays is not intended.)

If this isn't resonating with you, I would propose either:

An entirely separate ConstantArrayModel with the behavior Christoph & I are expecting
A note added to the documentation of ConstantModel warning that its return type when fed an array is different than many (all?) of the other builtin models

Maybe broadcasting the result against the data would work for both cases?

btw the plotting features are very useful, even in the current state, please do not remove them. They save a lot of typing and I often miss them when not working with Model.

On Tue, Dec 1, 2020 at 12:22 AM Jonathan Okasinski notifications@github.com wrote:

@newville https://github.com/newville

We could have multi-dimensional models or models that transform the data in some complex way (say Fourier)... We definitely want users to be able to do such things.

I completely agree that you don't want to restrict the user in principle. However, 1) my proposal doesn't affect the Model base class and 2) while this may be a documentation gap, I don't see anything in the models guide https://lmfit.github.io/lmfit-py/model.html, built-in models list https://lmfit.github.io/lmfit-py/builtin_models.html, or examples https://lmfit.github.io/lmfit-py/examples/index.html that suggest Model is intended for arbitrary data pre-processing in the manner you suggest here.

Huh, what arbitrary data pre-processing?. The Model class wraps a model function that takes some inputs (some of which will become variable parameters in the fit) and returns an output. That function might do complex calculations with the inputs.

I am pretty sure we never say that other models will behave consistently.

Again while this might be a documentation issue, as far as I can tell every piece of documentation in the models guide, built-in models list, and examples index show models behaving consistently with array inputs. Thus while not explicitly stated, I think it's fair to say that users may reasonably have this expectation.

I don't know what you think reasonable users will expect. Like, what should

 MyComplexModel.eval(my_test_params, x=7.3)  #assuming MyComplexModel

has independent_vars = ['x']

return? Do you want that to be an array?

The problem you see is...with PLOTTING.

My direct problem in usage is not plotting.

Um, no, this is definitely about plotting. Plotting is in the title of the Issue, and the traceback is because the evaluated model is not correctly coerced for plotting.

It is further processing. Given an array that I feed into a built-in lmfit

model, I do not know how to tell what data structure to expect back for further processing, nor how to interpret it. If this is intended to be done in some other way, please let me know. (For example, perhaps there's a nominal subclass of Model denoting 1D->1D array evals, as (to my knowledge) all the built-in models save Constant currently support)

Do we say somewhere in the document that the return type will be related to that of the independent variable? The fitting process needs a scalar float or a 1-D array of dtype 'float64' -- really, nothing else will do. Lmfit intercepts the data and the return value from Model.eval() and tries to coerce those into a 1-D float64 array and then try to handle NaNs (separate, except to emphasize that we do such interception).

The model function that the user provides may return a (float) scalar or a 1-D float64 array or something else - 2-D arrays, ints, complex values, etc. - maybe even Pandas series or other (maybe dask?) arrays. Lmfit does not reject those out of hand, but they might not all be obvious how to coerce into what the fitting process needs. We try to do such coercion, we might not always succeed. For sure, using numpy arrays are recommended.

But we do not ever say (I think) that the independent variable(s) must be 1-D or that the return value of the model will be the same shape or datatype as the independent variables(s). For sure, many simple curve-fitting problems have a single 1-d independent variable (say, x or t) and the model gives a 1-d array of the same length, implying "point by point" evaluation. With such simple models, we can readily infer how to "plot the data and fit". For sure, that's a useful feature -- we want Model.eval() to use the independent variables as they are passed in and do the calculation with those.

But, that is a feature of the individual models, not a feature of Model itself, and not a feature of all models.

(When I first had this issue, I looked into how lmfit might handle

different data structures return internally by checking plot. Given that all the other models I've tried worked in the same manner, I assumed's ConstantModel's different behavior with arrays is not intended.)

Well, the behavior of ConstantModel is intended: it returns a constant. It's kind of not rocket science ;).

The plot methods are convenient for simple cases, and such "simple 1-D cases" really do make up a lot of the usage here. But the plotting functions are limited and full of potential problems - this is not a plotting library. We chose one common plot library and use it to make the simple things simple. We cannot promise to cover every use-case, but if you want to fix a problem with plotting, that's fine.

If this isn't resonating with you, I would propose either:

An entirely separate ConstantArrayModel with the behavior Christoph & I are expecting

A note added to the documentation of ConstantModel warning that its return type when fed an array is different than many (all?) of the other builtin models

Shrug. The problem in this issue is with plotting - or with coercing data from a model so that it can be plotted. I'm confused why you keep proposing to change things other than where the problem actually is.

My deeper question is:

I do not know how to tell what data structure to expect back for further processing, nor how to interpret it.

I guess you've answered it:

The fitting process needs a scalar float or a 1-D array of dtype 'float64' -- really, nothing else will do.

To wrap things up:

Do we say somewhere in the document that the return type will be related to that of the independent variable?

No. However, as I have been saying, the documentation strongly implies that the return type will be related to that of the independent variable. This implication stems from the fact that it exclusively (to my knowledge) provides use-cases and examples in which they are the same. That's why I think that it should be noted in the documentation that this behavior is not guaranteed or canon.

Shrug. The problem in this issue is with plotting - or with coercing data from a model so that it can be plotted.

Exactly. The problem is interpreting data coming out of a Model.

To clarify this point:

Like, what should


has independent_vars = ['x']```
return?   Do you want that to be an array?

No. You've given a scalar input, getting a scalar back is fine. But this conversation has been about array inputs.

On Tue, Dec 1, 2020 at 10:46 AM Jonathan Okasinski notifications@github.com wrote:

My deeper question is:

I do not know how to tell what data structure to expect back for further processing, nor how to interpret it.

I guess you've answered it:

The fitting process needs a scalar float or a 1-D array of dtype 'float64' -- really, nothing else will do.

To wrap things up:

Do we say somewhere in the document that the return type will be related to that of the independent variable?

No. However, as I have been saying, the documentation strongly implies that the return type will be related to that of the independent variable. This implication stems from the fact that it exclusively (to my knowledge) provides use-cases and examples in which they are the same. That's why I think that it should be noted in the documentation that this behavior is not guaranteed or canon.

We have some built-in models, and we document those. I really don't think we imply anything. But, if you want to improve the documentation so that people do not incorrectly infer features that do not exist, that would be fine.

Shrug. The problem in this issue is with plotting - or with coercing data

from a model so that it can be plotted.

Exactly. The problem is interpreting data coming out of a Model.

Ok. So maybe "interpret them correctly"? Again, the problem here is really with plotting.

To clarify this point:

Like, what should MyComplexModel.eval(my_test_params, x=7.3) #assuming MyComplexModel has independent_vars = ['x'] return? Do you want that to be an array?

No. You've given a scalar input, getting a scalar back is fine.

Yup. And if you give a Pandas series or a 2-D array of complex, you might get those back as well. But that depends on what the model actually does. The return value of a model is dependent on what that model does. Whether and how it aligns with the independent variables is up to the model.

Your confusion sort of suggests to me that we should keep ConstantModel as it is. It also suggests to me that this is among the best arguments for having 2-D models. That is, for a 2-D gaussian (easily extended to lorentzian and voigt), the signature would be

def gaussian2d(y, x, amplitude, xcenter, ycenter, xsigma, ysigma, theta)

where x and y would typically be 1-d arrays of independent variables, that do not need to be the same size, with the rest being model parameters. The result would be a 2-D array, with shape (Ny, Nx). That would then give more models that did not map 1-d independent variables to a 1-D array of the same shape.

But this conversation has been about array inputs.

Nope, this issue is about plotting. It's right there in the subject line. ;)

Alright, I think I can see where you're going with this. I think the best examples would be of models using categorial data, like dictionaries, or something else very not-array, to emphasize that Models are completely generic.

In my opinion, rather than trying to emphasize this in the documentation, it might be more natural to explicit denote the return types in the Models. This would again emphasize the diversity of Models, but also allows users to know what to expect from working with any given model. This could be:

A nominal subclass of Model like Array1DModel, where all Array1DModels must map (n,) -> (n,) numpy arrays
An extra, optional field on Model like function_type that subclasses could use to describe their input & return types
Type hints provided by each model subclass (though, this might imply the models will coerce inputs)

Downstream processing code, and/or plot, could then read this information to figure out whether or not it knows how to interpret the Model outputs, without having to guess and deal with completely arbitrary data types coming from Model.

Let me know what you think about this.

@nitrocalcite

Alright, I think I can see where you're going with this.

Okay, just to be clear, this is still a discussion about an Issue that you raised about plotting of a ConstantModel. Like, I'm impressed with your general enthusiasm, but I may be more inclined to call this Issue more or less or done.

I think the best examples would be of models using categorial data, like dictionaries, or something else very not-array, to emphasize that Models are completely generic.

Well, maybe you have a natural example of such a model that would be generic enough to include in the models namespace - I don't have anything in mind that is both general and has an "obvious implementation" using a dictionary. Like, most simple curve-fitting problem really do use a 1-d array of "x" and the data is a 1-d array "y". Last I looked, scipy.optimize.curve_fit() only worked with 1 independent variable that was array-like. Apparently, lots of people use that routine - go figure ;).

In my opinion, rather than trying to emphasize this in the documentation, it might be more natural to explicit denote the return types in the Models.

Models are classes. What is 'return type'?

This would again emphasize the diversity of Models, but also allows users to know what to expect from working with any given model. This could be:

A nominal subclass of Model like Array1DModel, where all Array1DModels must map (n,) -> (n,) numpy arrays An extra, optional field on Model like function_type that subclasses could use to describe their input & return types Type hints provided by each model subclass (though, this might imply the models will coerce inputs)

Seems like a lot of work. I'm not sure I'm seeing where Array1DModel would be used. If type hints or extra fileds are optional, you'd have to handle them not being given, right?

Downstream processing code, and/or plot, could then read this information to figure out whether or not it knows how to interpret the Model outputs, without having to guess and deal with completely arbitrary data types coming from Model.

Let me know what you think about this.

I'd say "maybe". What strikes me most is that you keep proposing complex changes to the code for an Issue that I think is just not in need fixing at all. Maybe I'm just not getting what you're talking about.

lmfit / lmfit-py