Model capabilities on multiple datasets

[Andrea Munroe]

Hello all,

I have a code on my local machine that gives me x,y values from a set of parameters and the energy. For a single energy, I wrote a black box with Model that can fit a single dataset. It worked great!

I'm currently trying to build up to an energy-dependent model, which means I have to fit several datasets at once with different energies but the same set of parameters. Based on several of the posted examples, I found that I could get it to work with the minimize() function, but I lose the functionality of the ModelResult class. I would really like to recover the eval() and eval_uncertainties() functions so I can plot uncertainty bands. Is there a clever way I'm not seeing to move back to Model? 

Here's a streamlined example of my minimize() code, and the single-dataset version with model() below.

Many thanks!

Andrea

import lmfit
import numpy as np
from scipy.interpolate import CubicSpline

import data_wrapper # user defined code, black box
import run # user defined code, black box

def modelfit(x, energy,parameters,parameter_labels):
    model_x, model_y = run.run(energy,parameters,parameter_labels) # returns a list of xs, ys
    spl = CubicSpline(model_x,model_y)
    return spl(x) # match y values to input x values

def get_pars(datasets, initial_guesses, parameter_labels, rangemin, rangemax):
    fit_params = lmfit.Parameters()
    for parameter_index, parameter in enumerate(parameter_labels):
        myval = np.array(initial_guesses)[parameter_index]
        fit_params.add(f'{parameter}', value=myval, min=rangemin[parameter_index], max=rangemax[parameter_index])
    return fit_params
   
def objective(parameters, datasets, parameter_labels):
    resid = []
    for dataset_index, dataset in enumerate(datasets):
        input_parameters = []
        parameter_indices = list(parameters.keys())
        for parameter_index, parameter in enumerate(parameters):
            input_parameters.append(parameters[parameter_indices[parameter_index]].value)
        model_y = modelfit(dataset.x, dataset.energy, input_parameters, parameter_labels)
        resid.extend((dataset.y = model.y)/dataset.errs) # append and flatten residuals simultaneously
    return resid

parameter_labels = ["V","R","A"] # shared by all datasets
rangemin = [0,1,2]
rangemax = [3,4,5]
initial_guesses = run.lookup_literature_guesses(E, parameter_labels) # list of floats
datasets = data_wrapper.get_datasets() # array of objects

parameters = get_pars(datasets, initial_guesses, parameter_labels, rangemin, rangemax)
output = lmfit.minimize(objective, parameters, args=(datasets, parameter_labels))
lmfit.report_fit(output.params)

###############################################################################

# single-dataset version with Model

import matplotlib.pyplot as plt

def modelfit2(x, V, R, A, energy):
    pars = [V,R,A]
    parlabs = ["V","R","A"]
    model_x, model_y = run.run( energy, pars, parlabs)
    spl = CubicSpline(model_x, model_y)

def evaluate_model(y,x,errs,parameters,rangemin, rangemax, energy):
    model = lmfit.Model(modelfit2)
    pars = mod.make_params(
        V={'value':parameters[0],'min':rangemin[0],'max':rangemax[0]},
        R={'value':parameters[1],'min':rangemin[1],'max':rangemax[1]},
        A={'value':parameters[2],'min':rangemin[2],'max':rangemax[2]},
        energy={'value':energy, 'vary':False}
    )
    return model.fit(y, pars, x=x, method='leastsq', weights=1/errs)
ds = datasets[0]
result = evaluate_model(ds.y, ds.x, ds.errs, initial_guesses, ds.energy)
xfine = ds.x
yfine = result.eval(t=ds.x)
efine = result.eval_uncertainty(t=xfine,sigma=1)
fitrange = plt.fill_between(xfine,yfine-efine, yfine+efine)

[Matt Newville]

Hi Andrea, 

Yes, there is not (currently) an easy way to use the Model class for simultaneous refinement with multiple datasets.   In my own codes, I've sort of worked around this, but I agree that having to use "Minimize" is not ideal.

I think it would not be too hard to add support for this.  It would require some new code - but mostly bookkeeping to keep the connection between each dataset and its Model.  I think this would be a nice addition.   

I believe the needed code would really be to create a "Dataset" that has data, independent data, a model instance, and a "label".  That label would be used as part of the prefix for Parameters.   So, one might make a dictionary of Datasets with the label as key, and then manage a set of Parameters for all the Datasets  - but that would allow shared or constrained Parameter values across multiple datasets, and still have individual Models per dataset.

I would be reluctant to promise how quickly this could be done, but if you or anyone else is interested, I think it is worth working on.  Maybe we should move this to Github Discussions. 

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/lmfit-py/c5842575-1ab5-446e-9516-d89142b4481an%40googlegroups.com.

--

--Matt Newville <newville at cars.uchicago.edu> 630-327-7411

[Laurence Lurio]

Matt,

This would definitely be a feature I would use. 

Larry 

To view this discussion visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESbpMhQO6jANzW2mg%3DG0eCs8NULwwWmyHQ-A-DEMhkPYdyQ%40mail.gmail.com.

[Matthew Newville]

Hi Larry,

 

Thanks – yes, me too.  I think we should do this.  I started a Github Discussion at https://github.com/orgs/lmfit/discussions/1012  which might be a good place for folks to make comments or suggestions.

 

--Matt

 

 

To view this discussion visit https://groups.google.com/d/msgid/lmfit-py/CAK4HuLF0QZNi%2BUeTdiWL0wP3jK5obaxg78QQkq094nB5xsOZkQ%40mail.gmail.com.