Water Teixeira ∗ Resolution with bumps (2 wavelengths)
Introduction
The objective of this notebook is to show how to use the water_teixeira model and fit the data using bumps.
The data are two sets of water data measured at IN5 (ILL) using two different wavelengths, 5 and 8 Å.
Reference: J. Qvist, H. Schober and B. Halle, J. Chem. Phys. 134, 144508 (2011)
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",
'scale': "unit_of_signal/ps",
'center': "1/ps",
'D': "Angstrom^2/ps",
'radius': "Angstrom",
'resTime': "ps"}
Import libraries
[2]:
import ipywidgets
# the following two lines are to remove the warning about too many figures open simultaneously
from matplotlib import rcParams
rcParams.update({'figure.max_open_warning': 0})
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/'
# Data
# Wavelength 5 Angstrom
with h5py.File(path_to_data + 'H2O_293K_5A.hdf', 'r') as f:
hw_5A = f['entry1']['data1']['X'][:]
q_5A = f['entry1']['data1']['Y'][:]
unit_w5A = f['entry1']['data1']['X'].attrs['long_name']
unit_q5A = f['entry1']['data1']['Y'].attrs['long_name']
sqw_5A = np.transpose(f['entry1']['data1']['DATA'][:])
err_5A = np.transpose(f['entry1']['data1']['errors'][:])
# Wavelength 8 Angstrom
with h5py.File(path_to_data + 'H2O_293K_8A.hdf', 'r') as f:
hw_8A = f['entry1']['data1']['X'][:]
q_8A = f['entry1']['data1']['Y'][:]
unit_w8A = f['entry1']['data1']['X'].attrs['long_name']
unit_q8A = f['entry1']['data1']['Y'].attrs['long_name']
sqw_8A = np.transpose(f['entry1']['data1']['DATA'][:])
err_8A = np.transpose(f['entry1']['data1']['errors'][:])
# Resolution
# Wavelength 5 Angstrom
with h5py.File(path_to_data + 'V_273K_5A.hdf', 'r') as f:
res_5A = np.transpose(f['entry1']['data1']['DATA'][:])
# Wavelength 8 Angstrom
with h5py.File(path_to_data + 'V_273K_8A.hdf', 'r') as f:
res_8A = np.transpose(f['entry1']['data1']['DATA'][:])
# Force resolution function to have unit area
# 5 Angstrom
for i in range(len(q_5A)):
area = simps(res_5A[:,i], hw_5A)
res_5A[:,i] /= area
# 8 Angstrom
for i in range(len(q_8A)):
area = simps(res_8A[:,i], hw_8A)
res_8A[:,i] /= area
[4]:
fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(8, 8), sharex=True)
for i in range(len(q_5A)):
ax[0, 0].semilogy(hw_5A, sqw_5A[:,i], label=f"q={q_5A[i]:.1f}")
ax[1, 0].semilogy(hw_5A, res_5A[:,i], label=f"q={q_5A[i]:.1f}")
ax[0, 0].set_title(r'Signal 5 $\AA$')
ax[0, 0].grid()
ax[1, 0].set_title(r'Resolution 5 $\AA$')
ax[1, 0].set_xlabel(f"$\hbar \omega$")
ax[1, 0].grid()
for i in range(len(q_8A)):
ax[0, 1].semilogy(hw_8A, sqw_8A[:,i], label=f"q={q_8A[i]:.1f}")
ax[1, 1].semilogy(hw_8A, res_8A[:,i], label=f"q={q_8A[i]:.1f}")
ax[1, 1].legend(bbox_to_anchor=(1.2,1.5), loc='upper right', shadow=True)
ax[0, 1].set_title(r'Signal 8 $\AA$')
ax[0, 1].grid()
ax[1, 1].set_title(r'Resolution 8 $\AA$')
ax[1, 1].set_xlabel(f"$\hbar \omega$")
ax[1, 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"At 5 Angstroms, the names and units of `w` ( `x`axis) and `q` are: {unit_w5A[0].decode()} and {unit_q5A[0].decode()}, respectively.")
print(f"At 8 Angstroms, the names and units of `w` ( `x`axis) and `q` are: {unit_w8A[0].decode()} and {unit_q8A[0].decode()}, respectively.")
At 5 Angstroms, the names and units of `w` ( `x`axis) and `q` are: Energy Transfer (meV) and Wavevector Transfer (A!U-1!N), respectively.
At 8 Angstroms, the names and units of `w` ( `x`axis) and `q` are: Energy Transfer (meV) and Wavevector Transfer (A!U-1!N), respectively.
Create fitting model
[6]:
# Fit range -1 to +1 meV
idx_5A = np.where(np.logical_and(hw_5A > -1.0, hw_5A < 1.0))
idx_8A = np.where(np.logical_and(hw_8A > -1.0, hw_8A < 1.0))
def model_convol(x, q, scale=1, center=0, D=1, resTime=1, radius=1, DR=1, resolution=None):
model = QENSmodels.sqwWaterTeixeira(x, q, scale, center, D, resTime, radius, DR)
return np.convolve(model, resolution/resolution.sum(), mode='same')
model_all_qs = []
# First dataset: wavelength=5 Angstrom
for i in range(len(q_5A)):
x = hw_5A[idx_5A]
data = sqw_5A[idx_5A, i]
error = err_5A[idx_5A, i]
resol = res_5A[idx_5A, i]
# Select only valid data (error = -1 for Q, w points not accessible)
valid = np.where(error > 0.0)
x = x[valid[1]]
if len(valid[1]) != len(x):
print(i, "truncate to make vectors symmetric with respect to max")
data = data[valid]
error = error[valid]
resol = resol[valid]
# Teixeira model
model_q = bmp.Curve(
model_convol,
x, data, error,
name=f'q5A_{q_5A[i]:.2f}',
q=q_5A[i],
scale=10,
center=0.0,
D=0.13,
resTime=0.1,
radius=1.0,
DR=0.3,
resolution=resol
)
# Fitted parameters
model_q.scale.range(1.e-12, 1e2)
model_q.center.range(-0.1, 0.1)
model_q.D.range(0.05, 0.25)
model_q.resTime.range(0, 1)
model_q.radius.range(0.9, 1.1)
model_q.DR.range(1e-12, 1)
# Q-independent parameters
if i == 0:
D_q = model_q.D
T_q = model_q.resTime
R_q = model_q.radius
DR_q = model_q.DR
else:
model_q.D = D_q
model_q.resTime = T_q
model_q.radius = R_q
model_q.DR = DR_q
model_all_qs.append(model_q)
# Second dataset: wavelength=8 Angstrom
for i in range(len(q_8A)):
x = hw_8A[idx_8A]
data = sqw_8A[idx_8A, i]
error = err_8A[idx_8A, i]
resol = res_8A[idx_8A, i]
# Select only valid data (error = -1 for Q, w points not accessible)
valid = np.where(error > 0.0)
if len(valid[1]) != len(x):
print(i, "truncate to make vectors symmetric with respect to max")
x = x[valid[1]]
data = data[valid]
error = error[valid]
resol = resol[valid]
model_q = bmp.Curve(
model_convol,
x, data, error,
name=f'q8A_{q_8A[i]:.2f}',
q=q_8A[i],
scale=10,
center=0.0,
D=0.13,
resTime=0.1,
radius=1.0,
DR=0.3,
resolution=resol
)
# Fitted parameters
model_q.scale.range(1e-12, 1e2)
model_q.center.range(-0.1, 0.1)
model_q.D.range(0.05, 0.25)
model_q.resTime.range(0, 1)
model_q.radius.range(0.9, 1.1)
model_q.DR.range(1e-12, 1)
# Q-independent parameters set with 8A data
model_q.D = D_q
model_q.resTime = T_q
model_q.radius = R_q
model_q.DR = DR_q
model_all_qs.append(model_q)
problem = bmp.FitProblem(model_all_qs)
0 truncate to make vectors symmetric with respect to max
9 truncate to make vectors symmetric with respect to max
10 truncate to make vectors symmetric with respect to max
Display initial configuration: experimental data, fitting model with initial guesses
[7]:
slider = ipywidgets.IntSlider(value=0, min=0, max=len(q_5A)+len(q_8A)-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$")
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]:
[8]:
problem.summarize().splitlines()
[8]:
[' D ...|...... 0.13 in (0.05,0.25)',
' DR ..|....... 0.3 in (1e-12,1)',
' center ....|..... 0 in (-0.1,0.1)',
' radius ....|..... 1 in (0.9,1.1)',
' resTime |......... 0.1 in (0,1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)',
' center ....|..... 0 in (-0.1,0.1)',
' scale |......... 10 in (1e-12,100)']
Choice of minimizer for bumps
[9]:
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
[9]:
Number of steps for running fit using bumps
[10]:
steps_fitting = ipywidgets.IntText(
value=100,
step=100,
description='Number of steps when fitting',
style={'description_width': 'initial'})
steps_fitting
[10]:
Running the fit
Run the fit using the minimizer defined above with a number of steps also specified above.
[11]:
# Preview of the settings
print('Initial chisq', problem.chisq_str())
Initial chisq 5235.127(13)
[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 ...|...... 0.129307 in (0.05,0.25)',
' DR .|........ 0.186859 in (1e-12,1)',
' center .....|.... 0.00367957 in (-0.1,0.1)',
' radius |......... 0.9 in (0.9,1.1)',
' resTime ....|..... 0.469848 in (0,1)',
' scale .|........ 10.0037 in (1e-12,100)',
' center .....|.... 7.0008e-05 in (-0.1,0.1)',
' scale .|........ 10.9983 in (1e-12,100)',
' center .....|.... 0.000404934 in (-0.1,0.1)',
' scale .|........ 11.0133 in (1e-12,100)',
' center .....|.... 0.000185203 in (-0.1,0.1)',
' scale .|........ 10.9899 in (1e-12,100)',
' center ....|..... -0.0014919 in (-0.1,0.1)',
' scale .|........ 10.8764 in (1e-12,100)',
' center .....|.... 6.32419e-05 in (-0.1,0.1)',
' scale .|........ 10.7586 in (1e-12,100)',
' center ....|..... -0.00132593 in (-0.1,0.1)',
' scale .|........ 10.9923 in (1e-12,100)',
' center ....|..... -0.000775034 in (-0.1,0.1)',
' scale .|........ 10.8514 in (1e-12,100)',
' center ....|..... -0.000329188 in (-0.1,0.1)',
' scale .|........ 11.0382 in (1e-12,100)',
' center ....|..... -0.000867221 in (-0.1,0.1)',
' scale .|........ 10.6375 in (1e-12,100)',
' center .....|.... 0.00126753 in (-0.1,0.1)',
' scale .|........ 10.8632 in (1e-12,100)',
' center .....|.... 0.00257032 in (-0.1,0.1)',
' scale .|........ 10.7416 in (1e-12,100)',
' center .....|.... 0.00654292 in (-0.1,0.1)',
' scale .|........ 10.7529 in (1e-12,100)',
' center .....|.... 0.000673311 in (-0.1,0.1)',
' scale .|........ 10.8618 in (1e-12,100)',
' center ....|..... -0.00171027 in (-0.1,0.1)',
' scale .|........ 10.7119 in (1e-12,100)',
' center .....|.... 0.00473395 in (-0.1,0.1)',
' scale .|........ 10.6957 in (1e-12,100)',
' center .....|.... 0.011566 in (-0.1,0.1)',
' scale .|........ 10.8829 in (1e-12,100)',
' center .....|.... 0.000460568 in (-0.1,0.1)',
' scale |......... 9.24631 in (1e-12,100)',
' center ....|..... -0.00104432 in (-0.1,0.1)',
' scale |......... 9.27992 in (1e-12,100)',
' center ....|..... -0.00128986 in (-0.1,0.1)',
' scale |......... 9.84988 in (1e-12,100)',
' center ....|..... -0.00313492 in (-0.1,0.1)',
' scale .|........ 10.0664 in (1e-12,100)',
' center ....|..... -0.00471765 in (-0.1,0.1)',
' scale .|........ 10.7211 in (1e-12,100)',
' center ....|..... -0.00662987 in (-0.1,0.1)',
' scale .|........ 11.3506 in (1e-12,100)',
' center ....|..... -0.00990988 in (-0.1,0.1)',
' scale .|........ 11.6836 in (1e-12,100)',
' center ....|..... -0.0130454 in (-0.1,0.1)',
' scale .|........ 11.7779 in (1e-12,100)',
' center ....|..... -0.0167299 in (-0.1,0.1)',
' scale .|........ 11.9711 in (1e-12,100)',
' center ...|...... -0.0229142 in (-0.1,0.1)',
' scale .|........ 11.8458 in (1e-12,100)',
' center ..|....... -0.0514621 in (-0.1,0.1)',
' scale .|........ 11.9089 in (1e-12,100)']
[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):
print(k, ":", format_uncertainty_pm(v, dv))
final chisq 4273.718(13)
D : 0.12931 +/- 0.00017
DR : 0.1869 +/- 0.0017
center : 0.003680 +/- 0.000085
radius : 0.9000 +/- 0.0031
resTime : 0.4698 +/- 0.0072
scale : 10.004 +/- 0.012
center : 0.07e-3 +/- 0.11e-3
scale : 10.998 +/- 0.014
center : 0.40e-3 +/- 0.13e-3
scale : 11.013 +/- 0.014
center : 0.19e-3 +/- 0.19e-3
scale : 10.990 +/- 0.014
center : -0.00149 +/- 0.00021
scale : 10.876 +/- 0.014
center : 0.06e-3 +/- 0.15e-3
scale : 10.759 +/- 0.014
center : -0.00133 +/- 0.00029
scale : 10.992 +/- 0.014
center : -0.78e-3 +/- 0.33e-3
scale : 10.851 +/- 0.014
center : -0.33e-3 +/- 0.33e-3
scale : 11.038 +/- 0.014
center : -0.87e-3 +/- 0.44e-3
scale : 10.638 +/- 0.014
center : 0.00127 +/- 0.00055
scale : 10.863 +/- 0.015
center : 0.00257 +/- 0.00067
scale : 10.742 +/- 0.015
center : 0.00654 +/- 0.00070
scale : 10.753 +/- 0.015
center : 0.67e-3 +/- 0.84e-3
scale : 10.862 +/- 0.015
center : -0.00171 +/- 0.00083
scale : 10.712 +/- 0.016
center : 0.00473 +/- 0.00096
scale : 10.696 +/- 0.016
center : 0.01157 +/- 0.00099
scale : 10.883 +/- 0.016
center : 461e-6 +/- 25e-6
scale : 9.2463 +/- 0.0093
center : -0.001044 +/- 0.000038
scale : 9.280 +/- 0.011
center : -0.001290 +/- 0.000076
scale : 9.850 +/- 0.015
center : -0.00313 +/- 0.00014
scale : 10.066 +/- 0.020
center : -0.00472 +/- 0.00021
scale : 10.721 +/- 0.024
center : -0.00663 +/- 0.00028
scale : 11.351 +/- 0.026
center : -0.00991 +/- 0.00036
scale : 11.684 +/- 0.027
center : -0.01305 +/- 0.00045
scale : 11.778 +/- 0.027
center : -0.01673 +/- 0.00052
scale : 11.971 +/- 0.026
center : -0.02291 +/- 0.00062
scale : 11.846 +/- 0.026
center : -0.05146 +/- 0.00077
scale : 11.909 +/- 0.026
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_5A)+len(q_8A)-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$")
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]:
[ ]: