Brownian Translational diffusion ∗ Resolution with bumps
Introduction
The objective of this notebook is to show how to use the Brownian Translational diffusion model to perform some fits using bumps .
The reference data were generated data corresponding to a Brownian Translational diffusion model with self-diffusion coefficient = 0.145 Å\(^2\times\)meV.
The model was convoluted with a Gaussian resolution function of FWHM = 0.1 meV, centered randomly in the range [-0.01, +0.01] meV.
Finally the data are sampled randomly from a Poisson distribution.
The data do not have a background.
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 = {'D': "meV.Angstrom^2",
'scale': "unit_of_signal.meV",
'center': "meV"}
Import libraries
[2]:
import ipywidgets
import numpy as np
import h5py
import QENSmodels
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 + 'BrownianDiff_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 + 'BrownianDiff_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]:
my_linestyle = [
('solid', 'solid'),
('dotted', 'dotted'),
('dashed', 'dashed'),
('dashdot', 'dashdot'),
('long dash with offset', (5, (10, 3))),
('loosely dashdotted', (0, (3, 10, 1, 10))),
('densely dashdotted', (0, (3, 1, 1, 1))),
('dashdotdotted', (0, (3, 5, 1, 5, 1, 5))),
('loosely dashdotdotted', (0, (3, 10, 1, 10, 1, 10))),
('densely dashdotdotted', (0, (3, 1, 1, 1, 1, 1)))]
[5]:
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}", ls=my_linestyle[i][1])
ax[1].semilogy(hw, res[:,i], label=f"q={q[i]:.1f}", ls=my_linestyle[i][1])
ax[1].legend(bbox_to_anchor=(1, 1.15), labelspacing=0.25, handlelength=4, loc='upper right', shadow=True)
ax[0].set_title('Signal')
ax[0].grid()
ax[1].set_ylabel('$S(Q, \hbar \omega)$ (arb. u.)')
ax[0].set_ylabel('$S(Q, \hbar \omega)$ (arb. u.)')
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
[6]:
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
[7]:
# Fitting model
def model_convol(x, q, scale=1, center=0, D=1, resolution=None):
model = QENSmodels.sqwBrownianTranslationalDiffusion(x, q, scale, center, D)
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,
D=0.1,
resolution=res[:, i]
)
model_q.scale.range(0, 1e5)
model_q.center.range(-0.1, 0.1)
model_q.D.range(1e-12, 1)
# Q-independent parameters
if i == 0:
D_q = model_q.D
else:
model_q.D = D_q
model_all_qs.append(model_q)
problem = bmp.FitProblem(model_all_qs)
Display initial configuration: experimental data, fitting model with initial guesses
[8]:
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])
[8]:
[9]:
problem.summarize().splitlines()
[9]:
[' D |......... 0.1 in (1e-12,1)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1000 in (0,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1000 in (0,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1000 in (0,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1000 in (0,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1000 in (0,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1000 in (0,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1000 in (0,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1000 in (0,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1000 in (0,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1000 in (0,100000)']
Choice of minimizer for bumps
[10]:
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
[10]:
Setting for running bumps
[11]:
steps_fitting = ipywidgets.IntText(
value=100,
step=100,
description='Number of steps when fitting',
style={'description_width': 'initial'})
steps_fitting
[11]:
Running the fit
Run the fit using the minimizer defined above with a number of steps also specified above
[12]:
# Preview of the settings
print('Initial chisq', problem.chisq_str())
Initial chisq 429.0449(59)
[13]:
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]
[14]:
print((f"With {w_choice_minimizer.value} optimizer, "
f"you can use {settings_selected_optimizer(w_choice_minimizer)} as arguments of `fit`"))
With Levenberg-Marquardt optimizer, you can use ['steps', 'ftol', 'xtol'] as arguments of `fit`
[15]:
result = fitters.fit(
problem,
starts=10,
keep_best=True,
method=options_dict[w_choice_minimizer.value],
steps=int(steps_fitting.value)
)
Showing the results
[16]:
problem.summarize().splitlines()
[16]:
[' D .|........ 0.144129 in (1e-12,1)',
' center ....|..... -0.000482045 in (-0.1,0.1)',
' scale |......... 116.48 in (0,100000)',
' center .....|.... 0.000249078 in (-0.1,0.1)',
' scale |......... 155.009 in (0,100000)',
' center ....|..... -0.000935206 in (-0.1,0.1)',
' scale |......... 229.015 in (0,100000)',
' center ....|..... -0.000353388 in (-0.1,0.1)',
' scale |......... 332.641 in (0,100000)',
' center .....|.... 0.000167603 in (-0.1,0.1)',
' scale |......... 484.244 in (0,100000)',
' center .....|.... 0.000586575 in (-0.1,0.1)',
' scale |......... 676.617 in (0,100000)',
' center ....|..... -0.00213954 in (-0.1,0.1)',
' scale |......... 904.543 in (0,100000)',
' center .....|.... 0.000547481 in (-0.1,0.1)',
' scale |......... 1175.9 in (0,100000)',
' center .....|.... 0.000560597 in (-0.1,0.1)',
' scale |......... 1475.37 in (0,100000)',
' center ....|..... -0.000484356 in (-0.1,0.1)',
' scale |......... 1823.93 in (0,100000)']
[17]:
# 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.9995(59)
D : 0.14413 +/- 0.00029 meV.Angstrom^2
center : -0.48e-3 +/- 0.56e-3 meV
scale : 116.5 +/- 1.1 unit_of_signal.meV
center : 0.25e-3 +/- 0.56e-3 meV
scale : 155.0 +/- 1.3 unit_of_signal.meV
center : -0.94e-3 +/- 0.67e-3 meV
scale : 229.0 +/- 1.5 unit_of_signal.meV
center : -0.35e-3 +/- 0.85e-3 meV
scale : 332.6 +/- 1.9 unit_of_signal.meV
center : 0.2e-3 +/- 1.1e-3 meV
scale : 484.2 +/- 2.3 unit_of_signal.meV
center : 0.6e-3 +/- 1.1e-3 meV
scale : 676.6 +/- 2.7 unit_of_signal.meV
center : -2.1e-3 +/- 1.4e-3 meV
scale : 904.5 +/- 3.2 unit_of_signal.meV
center : 0.5e-3 +/- 1.7e-3 meV
scale : 1175.9 +/- 3.7 unit_of_signal.meV
center : 0.6e-3 +/- 2.0e-3 meV
scale : 1475.4 +/- 4.2 unit_of_signal.meV
center : -0.5e-3 +/- 2.0e-3 meV
scale : 1823.9 +/- 4.8 unit_of_signal.meV
Display final configuration: experimental data, fitting model with output of fitting for the refined parameters
[18]:
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])
[18]:
[ ]: