Delta and Lorentzian ∗ Resolution with bumps
Introduction
The objective of this notebook is to show how to use the model delta_lorentz corresponding to a delta representing a fraction \(p\) of fixed atoms and a Lorentzian corresponding to a Brownian Translational diffusion model for the remaining \((1-p)\) atoms
The reference data were generated using the above function with the following parameters \(A_0 = 0.3\) and \(D = 0.145\) Å\(^2\times\)meV.
The model is convoluted with a Gaussian resolution function of Full Width Half Maximum (FWHM) equal to 0.1 meV, centered randomly in the range [-0.01, +0.01] meV.
Finally the data are sampled randomly from a Poisson distribution.
There is no background.
The data are fitted with a general model of a delta + a Lorentzian, so the fitted parameters are not \(p\) and \(D\), but \(p\) and a \(Q\)-dependent HWHM.
Physical units
For information about unit conversion, please refer to the jupyter notebook called Convert_units.ipynb in the tools folder.
The dictionary of units defined in the cell below specify the units of the refined parameters adapted to the convention used in the experimental datafile.
[1]:
# Units of parameters for selected QENS model and experimental data
dict_physical_units = {'omega': "meV",
'q': "1/Angstrom",
'scale': "unit_of_signal.meV",
'center': "meV",
'hwhm': "meV"}
Import required libraries
[2]:
import ipywidgets
import h5py
import QENSmodels
import numpy as np
from scipy.integrate import simps
import bumps.names as bmp
from bumps import fitters
from bumps.formatnum import format_uncertainty_pm
import matplotlib.pyplot as plt
%matplotlib widget
Setting of fitting
Load reference data
[3]:
path_to_data = './data/'
# Read the sample
with h5py.File(path_to_data + 'DeltaBrownianDiff_Sample.hdf', 'r') as f:
hw = f['entry1']['data1']['X'][:]
q = f['entry1']['data1']['Y'][:]
unit_w = f['entry1']['data1']['X'].attrs['long_name']
unit_q = f['entry1']['data1']['Y'].attrs['long_name']
sqw = np.transpose(f['entry1']['data1']['DATA'][:])
err = np.transpose(f['entry1']['data1']['errors'][:])
# Read resolution
with h5py.File(path_to_data + 'DeltaBrownianDiff_Resol.hdf', 'r') as f:
res = np.transpose(f['entry1']['data1']['DATA'][:])
# Force resolution function to have unit area
for i in range(len(q)):
area = simps(res[:,i], hw)
res[:,i] /= area
[4]:
fig, ax = plt.subplots(nrows=2, sharex=True)
for i in range(len(q)):
ax[0].semilogy(hw, sqw[:,i], label=f"q={q[i]:.1f}")
ax[1].semilogy(hw, res[:,i], label=f"q={q[i]:.1f}")
ax[1].legend(bbox_to_anchor=(1.1,1.3), loc='upper right', shadow=True)
ax[0].set_title('Signal')
ax[0].grid()
ax[1].set_title('Resolution')
ax[1].set_xlabel(f"$\hbar \omega ({dict_physical_units['center']})$")
ax[1].grid()
Display units of input data
Just for information in order to determine if a conversion of units is required before using the QENSmodels
[5]:
print(f"The names and units of `w` (`x`axis) and `q` are: {unit_w[0].decode()} and {unit_q[0].decode()}, respectively.")
The names and units of `w` (`x`axis) and `q` are: X and Y, respectively.
Create fitting model
[6]:
# Fitting model
def model_convol(x, q, scale=1, center=0, A0=0, hwhm=1, resolution=None):
model = QENSmodels.sqwDeltaLorentz(x, q, scale, center, A0, hwhm)
return np.convolve(model, resolution/resolution.sum(), mode='same')
# Fit
model_all_qs = []
for i in range(len(q)):
# Bumps fitting model
model_q = bmp.Curve(
model_convol,
hw,
sqw[:,i],
err[:,i],
name=f'q_{q[i]:.2f}',
q=q[i],
scale=1000,
center=0.0,
A0=0.5,
hwhm=0.01,
resolution=res[:, i]
)
model_q.scale.range(0, 1e5)
model_q.center.range(-0.1, 0.1)
model_q.A0.range(0, 1)
model_q.hwhm.range(0, 2)
# Q-independent parameters
if i == 0:
A0_q = model_q.A0
else:
model_q.A0 = A0_q
model_all_qs.append(model_q)
problem = bmp.FitProblem(model_all_qs)
Display initial configuration: experimental data, fitting model with initial guesses
[7]:
slider = ipywidgets.IntSlider(value=0, min=0, max=len(q)-1, continuous_update=False)
output = ipywidgets.Output()
def fig_q(model, ax, q_index=0):
"""
Plot of experimental data, fitting model and residual for a selected q value
Parameters
----------
model: list of bumps.curve.Curves for all q
ax: matplotlib.axes to be updated when changing the ipywidgets
q_index: int
index of q to be plotted
"""
model = model[q_index]
ax[0].errorbar(model.x,
model.y,
yerr=model.dy,
label='experimental data',
color='C0')
ax[0].plot(model.x,
model.theory(),
label='theory (model)',
color='C1')
ax[0].set_title(f'Model {model.name} - $\chi^2$={problem.chisq_str()}')
ax[0].legend()
ax[1].plot(model.x, model.residuals(), marker='o', linewidth=0, markersize=3, color='C0')
with output:
fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True)
ax[0].grid(); ax[1].grid()
ax[1].set_ylabel('Residual')
ax[1].set_xlabel(f"$\hbar \omega ({dict_physical_units['center']})$")
fig_q(model_all_qs, ax, 0)
def update_profile(change):
"""
Update plots for a new q-value
"""
with output:
ax[0].clear(); ax[1].lines.clear()
ax[0].grid()
fig_q(model_all_qs, ax, change['new'])
slider.observe(update_profile, names="value")
slider_label = ipywidgets.Label("q value to display")
slider_comp = ipywidgets.HBox([slider_label, slider])
ipywidgets.VBox([slider_comp, output])
[7]:
Choice of minimizer for bumps
[8]:
options_dict = {}
for item in fitters.__dict__.keys():
if item.endswith('Fit') and fitters.__dict__[item].id in fitters.FIT_AVAILABLE_IDS:
options_dict[fitters.__dict__[item].name] = fitters.__dict__[item].id
w_choice_minimizer = ipywidgets.Dropdown(
options=list(options_dict.keys()),
value='Levenberg-Marquardt',
description='Minimizer:',
layout=ipywidgets.Layout(height='40px'))
w_choice_minimizer
[8]:
Setting for running bumps
[9]:
steps_fitting = ipywidgets.IntText(
value=100,
step=100,
description='Number of steps when fitting',
style={'description_width': 'initial'})
steps_fitting
[9]:
[10]:
# Preview of the settings
print('Initial chisq', problem.chisq_str())
Initial chisq 1474.2217(86)
Running the fit
Run the fit using the minimizer defined above with a number of steps also specified above.
[11]:
def settings_selected_optimizer(chosen_minimizer):
"""
List the settings available for the selected optimizer
This list can be used as arguments for the `fit` function
"""
assert type(chosen_minimizer) == ipywidgets.widgets.widget_selection.Dropdown
for item in fitters.__dict__.keys():
if item.endswith('Fit') and \
fitters.__dict__[item].id == options_dict[chosen_minimizer.value]:
return [elt[0] for elt in fitters.__dict__[item].settings]
print((f"With {w_choice_minimizer.value} optimizer, "
f"you can use {settings_selected_optimizer(w_choice_minimizer)} as arguments of `fit`"))
result = fitters.fit(
problem,
starts=10,
keep_best=True,
method=options_dict[w_choice_minimizer.value],
steps=int(steps_fitting.value))
With Levenberg-Marquardt optimizer, you can use ['steps', 'ftol', 'xtol'] as arguments of `fit`
Showing the results
[12]:
problem.summarize().splitlines()
[12]:
[' A0 ...|...... 0.302293 in (0,1)',
' center .....|.... 0.00201417 in (-0.1,0.1)',
' hwhm |......... 0.00530068 in (0,2)',
' scale |......... 112.518 in (0,100000)',
' center ....|..... -0.00054969 in (-0.1,0.1)',
' hwhm |......... 0.0208184 in (0,2)',
' scale |......... 135.404 in (0,100000)',
' center ....|..... -0.00102857 in (-0.1,0.1)',
' hwhm |......... 0.0485962 in (0,2)',
' scale |......... 166.099 in (0,100000)',
' center ....|..... -0.0022506 in (-0.1,0.1)',
' hwhm |......... 0.0853624 in (0,2)',
' scale |......... 200.407 in (0,100000)',
' center ....|..... -0.00110115 in (-0.1,0.1)',
' hwhm |......... 0.137085 in (0,2)',
' scale |......... 233.516 in (0,100000)',
' center ....|..... -0.00120486 in (-0.1,0.1)',
' hwhm |......... 0.191934 in (0,2)',
' scale |......... 256.507 in (0,100000)',
' center ....|..... -0.00499987 in (-0.1,0.1)',
' hwhm .|........ 0.267982 in (0,2)',
' scale |......... 271.513 in (0,100000)',
' center .....|.... 0.00367865 in (-0.1,0.1)',
' hwhm .|........ 0.36417 in (0,2)',
' scale |......... 292.197 in (0,100000)',
' center ....|..... -0.000534443 in (-0.1,0.1)',
' hwhm ..|....... 0.4487 in (0,2)',
' scale |......... 297.888 in (0,100000)',
' center ....|..... -0.00119294 in (-0.1,0.1)',
' hwhm ..|....... 0.552777 in (0,2)',
' scale |......... 308.972 in (0,100000)']
[13]:
# Other method to display the results of the fit (chi**2 and parameters' values)
print("final chisq", problem.chisq_str())
for k, v, dv in zip(problem.labels(), result.x, result.dx):
if k in dict_physical_units.keys():
print(k, ":", format_uncertainty_pm(v, dv), dict_physical_units[k])
else:
print(k, ":", format_uncertainty_pm(v, dv))
final chisq 0.9764(86)
A0 : 0.3023 +/- 0.0017
center : 0.00201 +/- 0.00072 meV
hwhm : 0.00530 +/- 0.00039 meV
scale : 112.5 +/- 1.1 unit_of_signal.meV
center : -0.55e-3 +/- 0.79e-3 meV
hwhm : 0.02082 +/- 0.00058 meV
scale : 135.4 +/- 1.2 unit_of_signal.meV
center : -0.00103 +/- 0.00099 meV
hwhm : 0.04860 +/- 0.00093 meV
scale : 166.1 +/- 1.3 unit_of_signal.meV
center : -2.3e-3 +/- 1.3e-3 meV
hwhm : 0.0854 +/- 0.0014 meV
scale : 200.4 +/- 1.4 unit_of_signal.meV
center : -1.1e-3 +/- 1.7e-3 meV
hwhm : 0.1371 +/- 0.0020 meV
scale : 233.5 +/- 1.6 unit_of_signal.meV
center : -1.2e-3 +/- 2.1e-3 meV
hwhm : 0.1919 +/- 0.0027 meV
scale : 256.5 +/- 1.6 unit_of_signal.meV
center : -5.0e-3 +/- 2.9e-3 meV
hwhm : 0.2680 +/- 0.0037 meV
scale : 271.5 +/- 1.7 unit_of_signal.meV
center : 3.7e-3 +/- 3.6e-3 meV
hwhm : 0.3642 +/- 0.0049 meV
scale : 292.2 +/- 1.8 unit_of_signal.meV
center : -0.5e-3 +/- 5.2e-3 meV
hwhm : 0.4487 +/- 0.0062 meV
scale : 297.9 +/- 1.9 unit_of_signal.meV
center : -1.2e-3 +/- 5.3e-3 meV
hwhm : 0.5528 +/- 0.0077 meV
scale : 309.0 +/- 2.0 unit_of_signal.meV
Display final configuration: experimental data, fitting model with output of fitting for the refined parameters
[14]:
slider1 = ipywidgets.IntSlider(value=0, min=0, max=len(q)-1, continuous_update=False)
output1 = ipywidgets.Output()
with output1:
fig1, ax1 = plt.subplots(nrows=2, ncols=1, sharex=True)
ax1[0].grid(); ax1[1].grid()
ax1[1].set_ylabel('Residual')
ax1[1].set_xlabel(f"$\hbar \omega ({dict_physical_units['center']})$")
fig_q(model_all_qs, ax1, 0)
def update_profile1(change):
"""
Update plots for a new q-value
"""
with output1:
ax1[0].clear(); ax1[1].lines.clear()
ax1[0].grid()
fig_q(model_all_qs, ax1, change['new'])
slider1.observe(update_profile1, names="value")
slider1_label = ipywidgets.Label("q value to display")
slider1_comp = ipywidgets.HBox([slider1_label, slider1])
ipywidgets.VBox([slider1_comp, output1])
[14]:
[ ]: