fitting a CompositeModel to multiple peaks

[francesco femi marafatto]

Hello All,

First of all, I really like this library, which I am using to fit synchrotron X-ray Emission spectra. 

I am trying to use the library to fit multiple species, where each species is adequately represented by a double pseudo-voigt peak (kalpha1 and kalpha2 fluorescence). I managed to fit a single species by a CompositeModel generated according to the documentation, but now I would like to use the fit parameters and the composite model to fit multiple species in a linear combination, but I am having trouble making lmfit work. namely, it is complaining when I try to read a composite model in a function and then add two composite models together.

here below the functions I defined:

def pseudo_voigt_doublet(x, data, comp=1, cal_shift=17.5, **kwargs):
pv1 = PseudoVoigtModel(prefix=f'comp{comp}_ka1_')
pars = pv1.guess(data, x=x)
pv2 = PseudoVoigtModel(prefix=f'comp{comp}_ka2_')
pars.update(pv2.make_params())
# change the center to be shifted to the right by (self.e_diff*pix_cal/2)
pars.add(f'comp{comp}_shift', value=cal_shift, min=12, max=25)
pars[f'comp{comp}_ka2_center'].set(
expr=f'comp{comp}_ka1_center - comp{comp}_shift') 
pars.add(f'comp{comp}_amp_ratio', value=2, min=1.7, max=2.3)
pars[f'comp{comp}_ka2_amplitude'].set(expr=f'comp{comp}_ka1_amplitude / comp{comp}_amp_ratio')
pars[f'comp{comp}_ka2_sigma'].set(value=pars[f'comp{comp}_ka1_sigma'], vary=False)
if 'sigma' in kwargs:
pars[f'comp{comp}_ka1_sigma'].set(value=kwargs.get('sigma'), vary=False)
mod = pv1 + pv2
return mod, pars

def two_components_1(model, fit_params, width):
mod1 = model
mod1.prefix('comp1_')
pars = mod1.make_params(**fit_params.valuesdict())
mod2 = model
mod2.prefix('comp2_')
pars.update(mod2.make_params(**fit_params.valuesdict()))
# set all params to not vary
for param in pars:
param.set(vary=False)
pars.add('width', value=width, min=width / 2, max=width * 2) # arbitrary values for width bounds
pars['comp1_ka1_center'].set(expr='comp1_ka1_center + comp1_ka1_center/width', vary=True)
pars['comp2_ka1_center'].set(expr='comp2_ka1_center - comp2_ka1_center/width', vary=True)
sum_mod = mod1 + mod2
return sum_mod, pars

My main issue is how to provide the variable "model" to the second function so that it can be read in as a model and then combined in the following steps. 

There is also probably an issue in assigning the parameters, in that I have not tested yet the **fit_params.valuesdict() as an input to the make_params method.

thanks and apologies if the issue is trivial, but I feel I am in a "dog biting it's tail" kind of situation...

Francesco

[francesco femi marafatto]

I forgot to add: the reason why I fit in two steps is that the first composite model should be very similar to the second one, with the only difference being the center position and the relative weight (the changes in the other parameters between different species are minor, I estimate within the uncertainty of the data).  

That, and because there may arise the possibility of there being more than 2 species, is why I did not try to include the 4 pseudo-voigt peaks in one function but tried to split the fitting procedures into two steps (if it works, the second function could potentially be generalized to n species)

[francesco femi marafatto]

And I forgot the actual code and traceback:

mod, pars = pseudo_voigt_doublet(x_arr, data)
out = mod.fit(data], pars, x=x_arr)
mod2, pars2 = two_components_1(mod, out.params, width)
out2 = mod2.fit(data, pars2, x=x_arr)

Traceback (most recent call last):

  File "<input>", line 1, in <module>

  File "C:\Program Files\JetBrains\PyCharm 2021.1.1\plugins\python\helpers\pydev\_pydev_bundle\pydev_umd.py", line 197, in runfile

    pydev_imports.execfile(filename, global_vars, local_vars)  # execute the script

  File "C:\Program Files\JetBrains\PyCharm 2021.1.1\plugins\python\helpers\pydev\_pydev_imps\_pydev_execfile.py", line 18, in execfile

    exec(compile(contents+"\n", file, 'exec'), glob, loc)

  File "C:/Users/marafatto_f/PycharmProjects/XES_1/main.py", line 27, in <module>

    mod2, pars2 = two_components_1(mod, out.params, width)

  File "C:\Users\marafatto_f\PycharmProjects\XES_1\functions.py", line 78, in two_components_1

    mod1.prefix('comp1_')

TypeError: 'str' object is not callable

[mpm...@gmail.com]

Dear Francesco,

What you're doing is creating an error at the line 

mod1.prefix('comp1_')

because

mod1.prefix
is a simple string label and not a function to set the string label

you want

mod1.prefix = 'comp1_'

Having said this, I think the approach is a bit strange. Even if the model parameters are related to one another I would still create the whole model and then assign the parameters.

You can make what I think is the model you want as follows:

import lmfit

from lmfit.models import PseudoVoigtModel

def mkPV(peak_ind, line):

    return PseudoVoigtModel(prefix=f'comp{peak_ind}_{line}_')

n = 3 # number of emission lines (pairs of ka1, ka2)

model = sum((mkPV(peak_ind, line)

             for peak_ind in range(1, n)

             for line in ['ka_1', 'ka_2']),

            mkPV(0, 'ka_1') + mkPV(0, 'ka_2'))

Best,

Mark

[francesco femi marafatto]

Dear Mark,

Thanks a lot for the answer! indeed, that was a dumb mistake higher in the line of bugs I had, and now it came back to the original issue:

Traceback (most recent call last):

  File "<input>", line 1, in <module>

  File "C:\Program Files\JetBrains\PyCharm 2021.1.1\plugins\python\helpers\pydev\_pydev_bundle\pydev_umd.py", line 197, in runfile

    pydev_imports.execfile(filename, global_vars, local_vars)  # execute the script

  File "C:\Program Files\JetBrains\PyCharm 2021.1.1\plugins\python\helpers\pydev\_pydev_imps\_pydev_execfile.py", line 18, in execfile

    exec(compile(contents+"\n", file, 'exec'), glob, loc)

  File "C:/Users/marafatto_f/PycharmProjects/XES_1/main.py", line 27, in <module>

    mod2, pars2 = two_components_1(mod, out.params, width)

  File "C:\Users\marafatto_f\PycharmProjects\XES_1\functions.py", line 81, in two_components_1

    mod2 = Model(model)

  File "C:\Users\marafatto_f\PycharmProjects\XES_1\venv\lib\site-packages\lmfit\model.py", line 276, in __init__

    self._parse_params()

  File "C:\Users\marafatto_f\PycharmProjects\XES_1\venv\lib\site-packages\lmfit\model.py", line 478, in _parse_params

    sig = inspect.signature(self.func)

  File "C:\Users\marafatto_f\AppData\Local\Programs\Python\Python39\lib\inspect.py", line 3130, in signature

    return Signature.from_callable(obj, follow_wrapped=follow_wrapped)

  File "C:\Users\marafatto_f\AppData\Local\Programs\Python\Python39\lib\inspect.py", line 2879, in from_callable

    return _signature_from_callable(obj, sigcls=cls,

  File "C:\Users\marafatto_f\AppData\Local\Programs\Python\Python39\lib\inspect.py", line 2254, in _signature_from_callable

    raise TypeError('{!r} is not a callable object'.format(obj))

TypeError: <lmfit.Model: (Model(pvoigt, prefix='comp1_ka1_') + Model(pvoigt, prefix='comp1_ka2_'))> is not a callable object

My approach was such because I need a first fit of the pseudo_voigt_doublet to constrain the parameters, and then use these constrained parameters in a linear combination for as many species are present (where then the only varying fit parameters would be the weight of each species - constituted of ka1 and ka2 - and the position of the center). I thought this approach would constrain a bit the model, which otherwise I fear would have been too unconstrained.

In summary my wished approach:

 

1. make an initial fit with 2 pseudo voigts on a dataset where I know there is only one species (not very elegant, but for now works)

2. use this initial fit to set the relationships between parameters for ka1 and ka2 --> save the fit parameters and the model

3. run a linear combination of the output model and parameters for as many species are present, allowing only the "center" position and the weight of each parameter to vary

with your solution, would it be possible to provide or fix the relationships between the ka1 and ka2 parameters for each species? I apologize if it is a trivial question.

Francesco

[francesco femi marafatto]

Ok please ignore the Traceback, I figured out the issue (clumsy code). My open question is still on the second part, i.e. the approach to the problem and how to implement it with lmfit.

thanks!

Francesco

[mpm...@gmail.com]

Dear Francesco,

You can add almost any constraint. 

I think you basically need to read this and the documentation, in general, to work out what to do. 

https://lmfit.github.io/lmfit-py/constraints.html?highlight=expression

As an example, this fixes the splitting for all ka_1 ka_2 pairs to peak_splitting_val=3. (It sounds like you want to first run a 2 PV fit and use this to determine the value for peak_splitting_val. 

import lmfit

from lmfit.models import PseudoVoigtModel

def mkPV(peak_ind, line):

    return PseudoVoigtModel(prefix=f'comp{peak_ind}_{line}_')

n = 3 # number of emission lines (pairs of ka1, ka2)

model = sum((mkPV(peak_ind, line)

             for peak_ind in range(1, n)

             for line in ['ka_1', 'ka_2']),

            mkPV(0, 'ka_1') + mkPV(0, 'ka_2'))

params = model.make_params()

peak_splitting_val = 3

params.add('peak_splitting', value=peak_splitting_val, vary=False)

for peak_ind in range(n):

    params[f'comp{peak_ind}_ka_2_center'].set(expr=f'comp{peak_ind}_ka_1_center + peak_splitting')

[francesco femi marafatto]

Dear Mark,

Thank you, I will try implementing your logic in my code, it does seem more simple this way.

best,

Francesco