Fitting linecuts using the lmfit package#

Import functions#

from nxs_analysis_tools.datareduction import load_data, Scissors
from nxs_analysis_tools.fitting import LinecutModel

from lmfit.models import GaussianModel, LinearModel

Load data#

data = load_data('example_data/sample_name/15/example_hkli.nxs')

data:NXdata
  @axes = ['H', 'K', 'L']
  @signal = 'counts'
  H = float64(100)
  K = float64(150)
  L = float64(200)
  counts = float64(100x150x200)

Perform linecuts#

s = Scissors(data=data)

linecut = s.cut_data(center=(0,0,0), window=(0.1,0.5,0.1))

Linecut axis: K
Integrated axes: ['H', 'L']

The `LinecutModel` class#

The workhorse of the fitting algorithm here is the LinecutModel class which wraps both the Model and ModelResult classes of the lmfit package, among other information.

lm = LinecutModel(data=linecut)

Create `lmfit` model#

Use the .set_model_components() method to set the model to be used for fitting the linecut. The model_components parameter must be a Model or list of Model objects.

lm.set_model_components([GaussianModel(prefix='peak'), LinearModel(prefix='background')])

Set model constraints#

Set constraints on the model using .set_param_hint().

lm.set_param_hint('peakcenter', min=-0.1, max=0.1)

After the model and the hints have been specified, use the .make_params() method to initialize the parameters.

Initialize parameters#

lm.make_params()

Parameters
name	value	initial value	min	max	vary	expression
peakamplitude	1.00000000	None	-inf	inf	True
peakcenter	0.00000000	None	-0.10000000	0.10000000	True
peaksigma	1.00000000	None	0.00000000	inf	True
backgroundslope	1.00000000	None	-inf	inf	True
backgroundintercept	0.00000000	None	-inf	inf	True
peakfwhm	2.35482000	None	-inf	inf	False	2.3548200*peaksigma
peakheight	0.39894230	None	-inf	inf	False	0.3989423*peakamplitude/max(1e-15, peaksigma)
peakcorrlength	2.66822318	None	-inf	inf	False	(2 * 3.141592653589793) / peakfwhm

Perform initial guess#

Use the .guess() method to perform an initial guess.

lm.guess()

The initial values of the parameters that result can be viewed using the .print_initial_params() method.

lm.print_initial_params()

peaksigma
	min: 0
peakfwhm
	expr: 2.3548200*peaksigma
peakheight
	expr: 0.3989423*peakamplitude/max(1e-15, peaksigma)
peakcenter
	min: -0.1
	max: 0.1

These can also be accessed through the individual LinecutModel objects.

lm.params

Parameters
name	value	initial value	min	max	vary	expression
peakamplitude	219.885366	219.88536567293792	-inf	inf	True
peakcenter	0.00000000	0.0	-inf	inf	True
peaksigma	0.06040268	0.06040268456375841	0.00000000	inf	True
backgroundslope	2.3853e-14	2.3853065274578794e-14	-inf	inf	True
backgroundintercept	180.019533	180.0195329458017	-inf	inf	True
peakfwhm	0.14223745	None	-inf	inf	False	2.3548200*peaksigma
peakheight	1452.27938	None	-inf	inf	False	0.3989423*peakamplitude/max(1e-15, peaksigma)
peakcorrlength	44.1739171	None	-inf	inf	False	(2 * 3.141592653589793) / peakfwhm

Visualize the initial guesses#

lm.plot_initial_guess()

../_images/123ae7acb74e43730dabb061a68ecd4dc7bc4253dc3287ce2e952f9e74fc7a7c.png

Perform the fit#

lm.fit()

Visualize the fit#

lm.plot_fit()

../_images/a310ff86305a0a6e44e3931c50c5805a0e52e584b437f41fef7ace596f106d07.png

[[Model]]
    (Model(gaussian, prefix='peak') + Model(linear, prefix='background'))
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 25
    # data points      = 76
    # variables        = 5
    chi-square         = 1.0034e-16
    reduced chi-square = 1.4132e-18
    Akaike info crit   = -3118.82112
    Bayesian info crit = -3107.16745
    R-squared          = 1.00000000
##  Warning: uncertainties could not be estimated:
    backgroundslope:      at initial value
[[Variables]]
    peakamplitude:        183.644087 (init = 219.8854)
    peakcenter:          -6.9385e-15 (init = 0)
    peaksigma:            0.06000000 (init = 0.06040268)
    backgroundslope:     -2.8140e-11 (init = 2.385307e-14)
    backgroundintercept: -1.3796e-09 (init = 180.0195)
    peakfwhm:             0.14128920 == '2.3548200*peaksigma'
    peakheight:           1221.05658 == '0.3989423*peakamplitude/max(1e-15, peaksigma)'
    peakcorrlength:       44.4703863 == '(2 * 3.141592653589793) / peakfwhm'

<Axes: xlabel='x', ylabel='y'>

The number of x values used to evaluate the fit curve can be set using the numpoints parameter, and any keyword arguments here will be passed to the .plot() method of the underlying ModelResult object.

lm.plot_fit(numpoints=1000)

../_images/21b50b08e88197d191c1a7bb022e5e2ba6b558c6919e75d5293d0d161e41dee3.png

[[Model]]
    (Model(gaussian, prefix='peak') + Model(linear, prefix='background'))
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 25
    # data points      = 76
    # variables        = 5
    chi-square         = 1.0034e-16
    reduced chi-square = 1.4132e-18
    Akaike info crit   = -3118.82112
    Bayesian info crit = -3107.16745
    R-squared          = 1.00000000
##  Warning: uncertainties could not be estimated:
    backgroundslope:      at initial value
[[Variables]]
    peakamplitude:        183.644087 (init = 219.8854)
    peakcenter:          -6.9385e-15 (init = 0)
    peaksigma:            0.06000000 (init = 0.06040268)
    backgroundslope:     -2.8140e-11 (init = 2.385307e-14)
    backgroundintercept: -1.3796e-09 (init = 180.0195)
    peakfwhm:             0.14128920 == '2.3548200*peaksigma'
    peakheight:           1221.05658 == '0.3989423*peakamplitude/max(1e-15, peaksigma)'
    peakcorrlength:       44.4703863 == '(2 * 3.141592653589793) / peakfwhm'

<Axes: xlabel='x', ylabel='y'>

Access the fit values#

The x and y values of the fit are stored in the .x (identical to the raw x values) and the .y_fit attributes.

lm.x

array([-0.503356, -0.489933, -0.47651 , ...,  0.47651 ,  0.489933,
        0.503356])

lm.y_fit

array([-1.364786e-09, -1.361743e-09, -1.341593e-09, ..., -1.368411e-09,
       -1.389316e-09, -1.393114e-09])

The residuals are stored within the ModelResult class of the lmfit package, which is stored in the .modelresult attribute of the LinecutModel class.

lm.modelresult.residual

array([-1.367040e-09, -1.366048e-09, -1.366214e-09, ..., -1.393032e-09,
       -1.393621e-09, -1.395368e-09])

Viewing the fit report#

lm.print_fit_report()

[[Model]]
    (Model(gaussian, prefix='peak') + Model(linear, prefix='background'))
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 25
    # data points      = 76
    # variables        = 5
    chi-square         = 1.0034e-16
    reduced chi-square = 1.4132e-18
    Akaike info crit   = -3118.82112
    Bayesian info crit = -3107.16745
    R-squared          = 1.00000000
##  Warning: uncertainties could not be estimated:
    backgroundslope:      at initial value
[[Variables]]
    peakamplitude:        183.644087 (init = 219.8854)
    peakcenter:          -6.9385e-15 (init = 0)
    peaksigma:            0.06000000 (init = 0.06040268)
    backgroundslope:     -2.8140e-11 (init = 2.385307e-14)
    backgroundintercept: -1.3796e-09 (init = 180.0195)
    peakfwhm:             0.14128920 == '2.3548200*peaksigma'
    peakheight:           1221.05658 == '0.3989423*peakamplitude/max(1e-15, peaksigma)'
    peakcorrlength:       44.4703863 == '(2 * 3.141592653589793) / peakfwhm'