Chudley-Elliott diffusion with bumps
Introduction
The objective of this notebook is to show how to use the Chudley Elliott diffusion model to perform some fits using bumps .
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': "1/ps",
'q': "1/Angstrom",
'D': "ps.Angstrom^2",
'L': "Angstrom",
'scale': "unit_of_signal/ps",
'center': "1/ps"}
Import required libraries
[2]:
import numpy as np
import ipywidgets
import QENSmodels
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
Create reference data
[3]:
nb_points = 500
xx = np.linspace(-4, 4, nb_points)
q = np.linspace(0.2, 2, 10)
added_noise = np.random.normal(0, 1, nb_points)
chudley_elliott_noisy = QENSmodels.sqwChudleyElliottDiffusion(
xx,
q,
scale=1.,
center=0.,
D=4,
L=10
) * (1. + 0.01 * added_noise) + 0.01 * added_noise
# arbitrary values for errors
err = [0.01 for i in range(nb_points)]
[4]:
fig = plt.figure()
for i in range(len(q)):
plt.plot(xx, chudley_elliott_noisy[i,:], label=f"q={q[i]:.1f}")
plt.legend()
plt.title('Signal')
plt.xlabel(f"$\hbar \omega $")
plt.grid()
Create fitting model
[5]:
model_all_qs = []
for i in range(len(q)):
# Bumps fitting model
model_q = bmp.Curve(
QENSmodels.sqwChudleyElliottDiffusion,
xx,
chudley_elliott_noisy[i],
err,
name=f'q_{q[i]:.2f}',
q=q[i],
scale=1,
center=0.,
D=2.,
L=8.
)
model_q.scale.range(0.1, 1e5)
model_q.center.range(-0.1, 0.1)
model_q.D.range(0.1, 5)
model_q.L.range(0.1, 15)
# Q-independent parameters
if i == 0:
D_q = model_q.D
L_q = model_q.L
else:
model_q.D = D_q
model_q.L = L_q
model_all_qs.append(model_q)
problem = bmp.FitProblem(model_all_qs)
Display initial configuration: experimental data, fitting model with initial guesses
[6]:
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])
[6]:
[7]:
problem.summarize().splitlines()
[7]:
[' D ...|...... 2 in (0.1,5)',
' L .....|.... 8 in (0.1,15)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)']
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]:
Running the fit
Run the fit using the minimizer defined above with a number of steps also specified above
[10]:
# Preview of the settings
print('Initial chisq', problem.chisq_str())
Initial chisq 72.4499(49)
[11]:
problem.summarize().splitlines()
[11]:
[' D ...|...... 2 in (0.1,5)',
' L .....|.... 8 in (0.1,15)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 1 in (0.1,100000)']
[12]:
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]
[13]:
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`
[14]:
result = fitters.fit(
problem,
starts=10,
keep_best=True,
method=options_dict[w_choice_minimizer.value],
steps=int(steps_fitting.value))
Showing the results
[15]:
problem.summarize().splitlines()
[15]:
[' D .......|.. 4.00284 in (0.1,5)',
' L ......|... 9.99903 in (0.1,15)',
' center ....|..... -0.0010826 in (-0.1,0.1)',
' scale |......... 1.00293 in (0.1,100000)',
' center .....|.... 0.000180073 in (-0.1,0.1)',
' scale |......... 1.00367 in (0.1,100000)',
' center ....|..... -0.000296364 in (-0.1,0.1)',
' scale |......... 1.0036 in (0.1,100000)',
' center ....|..... -0.000763784 in (-0.1,0.1)',
' scale |......... 1.00337 in (0.1,100000)',
' center ....|..... -0.000272319 in (-0.1,0.1)',
' scale |......... 1.00353 in (0.1,100000)',
' center ....|..... -0.000302171 in (-0.1,0.1)',
' scale |......... 1.0036 in (0.1,100000)',
' center ....|..... -0.000633794 in (-0.1,0.1)',
' scale |......... 1.00344 in (0.1,100000)',
' center ....|..... -0.000384421 in (-0.1,0.1)',
' scale |......... 1.00348 in (0.1,100000)',
' center ....|..... -0.000311514 in (-0.1,0.1)',
' scale |......... 1.00358 in (0.1,100000)',
' center ....|..... -0.000566816 in (-0.1,0.1)',
' scale |......... 1.00348 in (0.1,100000)']
[16]:
# 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 1.2992(49)
D : 4.0028 +/- 0.0087 ps.Angstrom^2
L : 9.999 +/- 0.011 Angstrom
center : -0.00108 +/- 0.00021 1/ps
scale : 1.0029 +/- 0.0015 unit_of_signal/ps
center : 0.18e-3 +/- 0.68e-3 1/ps
scale : 1.0037 +/- 0.0018 unit_of_signal/ps
center : -0.30e-3 +/- 0.56e-3 1/ps
scale : 1.0036 +/- 0.0018 unit_of_signal/ps
center : -0.76e-3 +/- 0.43e-3 1/ps
scale : 1.0034 +/- 0.0015 unit_of_signal/ps
center : -0.27e-3 +/- 0.57e-3 1/ps
scale : 1.0035 +/- 0.0017 unit_of_signal/ps
center : -0.30e-3 +/- 0.56e-3 1/ps
scale : 1.0036 +/- 0.0018 unit_of_signal/ps
center : -0.63e-3 +/- 0.47e-3 1/ps
scale : 1.0034 +/- 0.0016 unit_of_signal/ps
center : -0.38e-3 +/- 0.54e-3 1/ps
scale : 1.0035 +/- 0.0017 unit_of_signal/ps
center : -0.31e-3 +/- 0.56e-3 1/ps
scale : 1.0036 +/- 0.0017 unit_of_signal/ps
center : -0.57e-3 +/- 0.49e-3 1/ps
scale : 1.0035 +/- 0.0016 unit_of_signal/ps
Display final configuration: experimental data, fitting model with output of fitting for the refined parameters
[17]:
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])
[17]:
[ ]: