Support for fitting multiple datasets with Models

[newville]

There was a question and suggestion (https://groups.google.com/g/lmfit-py/c/9fxjGs5TpPo/m/J7K_7QVrDAAJ ) to be be able to use Model for fitting multiple datasets. I think we should try to support this, and propose that we make a few new classes:

class DataSet could have components of

• data: the numeric data
• model: an instance of a Model (including composites, etc)
• label: a string tag for the data that will be used for parameter names, so must match 'variable name' rules.
• appropriate independent variables as dictated by the model (say, x or theta)

This DataSet could have methods fit, make_params, eval, guess, and other methods associated with Model and ModelResult. Maybe even plot?

That is, we have been careful to have Model be generic, but if we make a DataSet that has a Model, I think that's OK.

class MultipleDataSets could then just be a dictionary of DataSets, keyed by label and also with methods ds fit, make_params, eval.

Here, parameters would be prefixed with the DataSet tag (and maybe model component tags). As a near concrete example, let's say there is a Line+Gaussian model for some temperature-dependent data, and one expects the center and sigma to move with temperature, but maybe in some expected (say, linear with temperature) way.

So, if you have Datasets with tags ('temp050', 'temp100', 'temp150', 'temp200', 'temp250', 'temp300'), you would end up with Parameters names like temp100_center , temp150_sigma, temp200_amplitude, and so forth.

By default, each of the Parameters would be independent, but the user could then impose constraints, perhaps something like:

params.add('center_0').add(value=500) # peak center extrapolated to t=0 
params.add('center_alpha').add(value=0.50)  # linear temperature coefficient for center

params['t050_center'].set(expr='center_0 +  50*center_alpha')
params['t100_center'].set(expr='center_0 + 100*center_alpha')
params['t150_center'].set(expr='center_0 + 150*center_alpha')
params['t200_center'].set(expr='center_0 + 200*center_alpha')
params['t250_center'].set(expr='center_0 + 250*center_alpha')
params['t300_center'].set(expr='center_0 + 300*center_alpha')

The result of the fit could then have a ModelResult per DataSet, as well as an overall ModelResult.

Does this seem reasonable? Comments? Suggestions?

[erik-mansson]

I actually have a working implementation of multi-signal fitting, used by me and a few colleagues since several years. I have hesitated to contribute it to lmfit because it relies on some __mro__-hacking to make my PrefixedWrapper class be able to make any Model-instance appear to have additional prefixes to its parameters, without any modification of code in the lmfit package. (That constraint of course disappears if we develop something new in the lmfit repository, but can still be a useful ambition to ensure no old functionality changes.)

Initially, my implementation was mixed with some custom models in an even larger Python file, which would make the scope of a patch huge, but I had started worked on creating a slightly cleaner branch and considering how to present it to the community. The recent mailing list question and you opening this discussion thread spurred me to spend an evening to prepare something I can share here. As it is an add-on rather than a change within the package, I simply attach two files here instead of making a pull request (zipped since .py-upload wasn't allowed). global_fit.zip contains global_fit.py which has the implementation but I would recommend to start reading demo_global_fit.py that has several examples of how it would be used.

In your draft above, I think your MultipleDataSets is similar to my MultiModel, but I did't introduce any DataSet or change to how data is passed for normal fitting of "single" models. MultiModel has a dict of underlying models (also called submodels) keyed by arbitrary string names and MultiModel.fit() takes a dict with matching keys for their respective data. If if the underlying models would use different names for their parameters and independent variables, nothing more needs to be done. To support more typical use-cases where the user wants to use e.g. a GaussianModel repeated 5 times, as on https://lmfit.github.io/lmfit-py/examples/example_fit_multi_datasets.html, I use a PrefixedWrapper that allows parameters to be renamed and any underlying model instance used differently for each dataset. To allow the submodels to use the same name for their independent variable, but different values (and array lengths), the fit()-method can take a dict of the different versions (with any specific keyword arguments that the submodels need). There is also an alternative constructor MultiModel.for_prefixes() that makes the typical use cases even easier, and which can automatically set some typical parameter constraints (expressions to share a value between submodels). See usage in last cell of the attached demo_global_fit.py.

In practice, I and my colleagues typically have data where each submodel fits 1D-signals of the same length, using the same independent variable array. Specifically, we have time-resolved mass spectra from pump–probe scans, and if you select the the ion yield for a handful of ion masses, this gives a handful of equally long signals as function of a common time coordinate. Instead of using the general MultiModel I then prefer to use a StackedModel class which requires that the individual datasets are 1D-signals of the same array length, which can be stacked to a 2D-array with a common coordinate for the common independent variable. I can then pass my data as a 2D-array and the prefixes can be automatically generated from a list (or xarray.DataArray coordinate) of the selected masses. Because the fit-parameter for calibration of the time axis (temporal overlap between pulses) is normally expected to be common for all signals, and in some cases also the sigma-paramer for an ExponentialGaussianModel may be shared by some signals, it makes better use of the statistics and gives more meaningful fits when being able to run a global fit rather than fitting each signal on its own. Most of the examples in demo_global_fit.py use StackedModel.

Parts of my hackish implementation of PrefixedWrapper and StackedModelResult would break if multiple threads would call for instance eval() in parallel on the same instance, but I don't think that's a huge problem (and it could probably be worked around if extra arguments and functionality would be included in relevant methods of some lmfit-clases). There may be some dead code and certainly some unimportant if-cases in my implementation, but I think it's worth a look as inspiration and for discussing approaches in this thread. I have some unit tests that exercise the branch I and my colleagues use, but I haven't worked on copying them to this "tidied up" branch. They would basically look similar to the sections in demo_global_fit.py, but with more assertions about the fit results and parameter constraints.

[newville]

@erik-mansson Thanks. I gave your global fit a brief read but not (yet?) a more careful study.

I think there is a lot of overlap with PrefixWrapper and what I sketched out.

I think we can make a new class or two specifically for this.

The idea of a "DataSet" that has a Model is sort of independent of multi-dataset fitting, but would make it a lot easier. And a "DataSet" might be useful even as a single thing. In my downstream codes for analysis and fitting, there is a loose definition of a "Data Group" - basically an empty container that expects naming conventions. But that allows a) multiple models to be perserved for a data set (compare Gaussian to Voigt, say) as a "fit history" for that data set.

That sort of organization isn't too hard to expect that downstream apps would just naturally do that. But, if we had a DataSet class with this partially builtin, that might be useful (say, to be saved and re-loaded).

One thing I would like to not have to assume is that Data Sets are the same length or shape, which I believe np.stack would imply.

But, again, thanks! I think that is very helpful, and I'll look more closely at your code.

[erik-mansson]

After having looked back at my code and at your ideas, I find three orthogonal topics:

1. Something that does the logistics of delegating to multiple underlying models for their respective piece of the multi-dataset, but within some wrapper that concatenates and appears like a regular 1D-array so that Model.fit() can run. The user must have provided models with parameter names that don't conflict unless they are meant to be shared.
2. Optional PrefixedWrapper so that the user can reuse a model instance defined once (also ExpressionModel which doesn't support any prefix argument), with different prefixes, which is useful in the context of global fitting of different datasets, but could have other uses too (e.g. to build a CompositeModel using standard models without parameter name conflicts).
3. Optional DataSet (and for point 1. also MultipleDataSets) approach to how arguments are prepared, how a fit is launched and how the results are kept in memory, as aid for organizing a workflow where for instance different ways of fitting are to be compared. (If I understood you right.)

My outlooks:

1. I will have a go at demonstrating a more basic function-based approach to 1. Although my previoysly uploaded implementation used custom classes, it seems to me now that a cleaner approach would just define a few new functions: multi_fit(data: dict[str,array_like], models: dict[str,Model], params: Parameters, weights:dict[str,array_like]=None, separate_kwargs:dict[str,dict]={}, **kwargs) would take string-keyed dicts for some arguments. Internally it would define an ad-hoc python function for the eval()-logistics, wrap it in a single-use Model-instance and run its fit() on suitably concatenated data. I would also add at least the utility functions multi_make_params() and multi_guess() to help the user build and merge the necessary parameter sets, but there won't be any management of prefixes or automatic assistance in setting expr on the parameters to apply grouping constraints. (It will be possible to add other utility functions for constraint management, without complicating the already mentioned functions, but they don't strictly need to be inside the lmfit library, so it could be left up to users how they want to manage that, made easier when PrefixedWrapper is available.) The somewhat cumbersome way of getting different independent variable arrays to different underlying datasets & models even when the variable name is the same (e.g. x, but different values for it) is via the separate_kwargs dict of dicts. Here, I see that a DataSet container (point 3 below) would help bundle data with its independent variables. But my personal preference would for that be to start supporting xarray.DataArray instead in the data-argument, which lets you have named coordinates (and other metadata) on top of a numpy.ndarray, and comes with nice automatic plotting that makes use of that information. I think a MultiModelResult class should be made for holding the results (although using dicts also for results would be possible).
2. While perhaps using unconventional methods, I'd argue my PrefixedWrapper solves topic 2 and could be moved out of the global_fit.py (not sure where). But that was not the main aim of this discussion thread.
3. Personally I don't see the benefit of 3 (yet) so I leave that to you or others, in this thread or elsewhere. (I'd say that ModelResult already bundles up a reference to the model, init_params and fitted params in a way that seems sufficient to me, getting a new ModelResult every time I run a different fit so there's no risk of changing DataSet.model but forgetting to run a DataSet.fit() and thus having DataSet.result refer to something outdated.)

[newville]

@erik-mansson Thanks. I am pretty sure we want some sort of Dataset container. You're right that this will have much in common with ModelResult, but it would also be intended for a user to create and then use, whereas ModelResult is generated as output (and not expected to be changed). Except for back-compatibility and that we'll continue to support the current Model/ModelResult approach, much of the content of ModelResult could be put into Dataset when a fit is finished. But, I would probably prefer Dataset to have a result attribute that is a ModelResult (as opposed to being a ModelResult). One reason for that is that a Dataset could have a list of those as a fit_history.

With Dataset, much of the logistics and naming conventions for curve-fitting become codified - a Dataset (like a Dataclass, which it might be able to use) essentially enforces naming conventions. This is important for multiple dataset fitting, as our code will have to find/set values for individual data sets. But I also think that having Dataset will make curve-fitting simpler for many folks.

A MultiDataset class could just a dictionary of Datasets, and its own "global" Parameters, where prefix-mangling would be done. This class can impose additional prefixes based on the dictionary keys for the datasets (dsetA_peak1_center, etc). After a fit, it could have a ModelResult, and so could each of its Datasets.

An example use case would be to creata Dataset A, fit that, then create Dataset B, and fit that. Then, a MultiDataset with Datasets A and B could be made, and a combined Parameters based on the best-fit Parameters for the datasets. Constraints could then be put on some of the parameters, and then the combined fit could be done. So having that ability to run a fit on a single Dataset seem important to me.

I am wary of using xarray. Independent variables are not required to be ndarrays or even array-like. DataFrames (from however many libraries now support that concept), HDF5 datasets, etc are definitely allowed as independent data. So are dicts and user-defined classes. We do not need to make any assumptions or assertions about what independent variables are.

[newville]

I started a PR (#1015), at this point only adding Dataset, which (as pointed out above) is really renaming ModelResult to Datatset and adding a label attribute.