{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " \n", "Disclaimer: \n", " \n", "The main objective of the Jupyter notebooks is to show how to use the models of the QENS library by\n", " \n", "- building a fitting model: composition of models, convolution with a resolution function \n", "- setting and running the fit \n", "- extracting and displaying information about the results \n", "\n", "These steps have a minimizer-dependent syntax. That is one of the reasons why different minimizers have been used in the notebooks provided as examples. \n", "But, the initial guessed parameters might not be optimal, resulting in a poor fit of the reference data.\n", "\n", "
\n", "\n", "# Water Teixeira ∗ Resolution with lmfit\n", "\n", "## Introduction\n", "\n", "
\n", " \n", "This example shows how to use the water_teixeira model and fit the data using lmfit.\n", "
\n", "The data are two sets of water data measured at IN5 (ILL) at 5 Å.\n", "\n", "**Reference:** J. Qvist, H. Schober and B. Halle, *J. Chem. Phys.* **134**, 144508 (2011)\n", "\n", "### Physical units\n", "\n", "For information about unit conversion, please refer to the jupyter notebook called `Convert_units.ipynb` in the `tools` folder.\n", "\n", "## Import libraries" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import h5py\n", "from scipy.integrate import simps\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# the following line is to remove the warning about too many figures open simultaneously\n", "plt. rcParams.update({'figure.max_open_warning': 0})\n", "\n", "import lmfit\n", "from scipy.interpolate import interp1d\n", "\n", "%matplotlib widget" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setting of fitting\n", "\n", "### Import reference data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "path_to_data = './data/'\n", "\n", "with h5py.File(path_to_data + 'H2O_293K_5A.hdf', 'r') as f:\n", " data_in = f['entry1']\n", " w = data_in['data1']\n", " x = w['X'][()] # energy or time values\n", " unit_w = w['X'].attrs['long_name']\n", " unit_q = w['Y'].attrs['long_name']\n", " y = w['DATA'][()] # intensities\n", " e = w['errors'][()] # errors for the intensities\n", " # Obtain the momentum transfer values\n", " q = w['Y'][()]\n", " data_5A = dict(q=q, x=x, y=y, e=e)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# number of spectra (i.e. number of different q-values)\n", "nb_q_values = len(data_5A['q'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Display units of input data\n", "Just for information in order to determine if a conversion of units is required before using the QENSmodels" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print((f\"The names and units of `w` (`x`axis) and `q` are: \" \n", " f\"{unit_w[0].decode()} and {unit_q[0].decode()}, respectively.\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Import resolution data and normalize (unit area)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "path_to_data = './data/'\n", "\n", "with h5py.File(path_to_data + 'V_273K_5A.hdf', 'r') as f:\n", " data = f['entry1']\n", " w = data['data1'] \n", " res_5A_x = w['X'][()]\n", " res_5A = np.transpose(w['DATA'][()])\n", "\n", "# Force resolution function to have unit area \n", "for i in range(len(data_5A['q'])):\n", " area = simps(res_5A[:,i], data_5A['x'])\n", " res_5A[:,i] /= area" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Mask data according to energy range and filter negative error" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Filter according to energy-range\n", "mask = np.intersect1d(np.where(data_5A['x']>-1.), np.where(data_5A['x']<1.))\n", "\n", "f_5A_mask = dict()\n", "f_5A_mask['x'] = np.asarray([data_5A['x'][mask] for i in range(nb_q_values)])\n", "f_5A_mask['y'] = np.asarray([y[mask] for y in data_5A['y']])\n", "f_5A_mask['e'] = np.asarray([e[mask] for e in data_5A['e']])\n", "\n", "# Select resolution according to energy range\n", "res_5A_x = res_5A_x[mask]\n", "res_5A = res_5A[mask, :]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Filter according to negative error values \n", "# resolution\n", "selected_indices = np.where(f_5A_mask['e'][i] > 0.0)\n", "resol_5A_x = np.asarray([res_5A_x[selected_indices] for i in range(nb_q_values)])\n", "resol_5A = np.asarray([res_5A[selected_indices, i][0] for i in range(nb_q_values)]) \n", " \n", "# data\n", "f_5A = dict()\n", "f_5A['x'] = np.asarray([x[selected_indices] for i, x in enumerate(f_5A_mask['x'])])\n", "f_5A['y'] = np.asarray([y[selected_indices] for i, y in enumerate(f_5A_mask['y'])])\n", "f_5A['e'] = np.asarray([e[selected_indices] for i, e in enumerate(f_5A_mask['e'])])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot experimental data\n", "fig0= plt.figure()\n", "[plt.semilogy(f_5A['x'][i], f_5A['y'][i]) for i in range(nb_q_values)]\n", "plt.xlabel(r'Energy transfer (meV)')\n", "plt.title('Reference data - 5 Angstrom')\n", "plt.grid();" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot experimental data\n", "fig1 = plt.figure()\n", "[plt.semilogy(resol_5A_x[i], resol_5A[i]) for i in range(nb_q_values)]\n", "plt.xlabel(r'Energy transfer (meV)') \n", "plt.title('Resolution function - 5 Angstrom')\n", "plt.grid(); " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create function for instrument resolution data (cubic interpolation between tabulated data points)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_interp = [interp1d(resol_5A_x[i], \n", " resol_5A[i]/np.sum(resol_5A[i]), \n", " kind='cubic', \n", " bounds_error=False, \n", " fill_value='extrapolate') for i in range(nb_q_values)]\n", "\n", "def irf_gate(w, spectrum_nb=0):\n", " \"\"\" Function defined from the interpolation of instrument resolution data \n", " Used to define fitting model and plot \"\"\" \n", " return f_interp[spectrum_nb](w)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# check interpolation for first spectrum of resolution function: plot tabulated data and interpolated data\n", "indx = 0\n", "\n", "fig2 = plt.figure()\n", "plt.plot(resol_5A_x[indx], \n", " resol_5A[indx]/np.sum(resol_5A[indx]), \n", " '.', \n", " label=f\"tabulated data. q={data_5A['q'][indx]:.2}\")\n", "plt.plot(f_5A['x'][indx], \n", " irf_gate(f_5A['x'][indx], indx), \n", " '--', \n", " label=f\"extrapolated data. q={data_5A['q'][indx]:.2}\")\n", "plt.legend(bbox_to_anchor=(1.0, .95))\n", "plt.xlabel('w')\n", "plt.title(f'Instrument resolution: tabulated data and interpolated data for spectrum {indx}')\n", "plt.grid();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create fitting model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import QENSmodels" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create convolution function \n", "# code from https://lmfit.github.io/lmfit-py/model.html\n", "\n", "def convolve(arr, kernel):\n", " # simple convolution of two arrays\n", " npts = min(len(arr), len(kernel))\n", " pad = np.ones(npts)\n", " tmp = np.concatenate((pad*arr[0], arr, pad*arr[-1]))\n", " \n", " out = np.convolve(tmp, kernel, mode='valid')\n", " noff = int((len(out) - npts)/2)\n", " return out[noff:noff+npts]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model = lmfit.CompositeModel(lmfit.Model(irf_gate), lmfit.Model(QENSmodels.sqwWaterTeixeira), convolve) \n", "\n", "print(f'Names of parameters: {model.param_names}\\nIndependent variable(s): {model.independent_vars}')\n", "\n", "# Define boundaries for parameters to be refined\n", "model.set_param_hint('scale', min=0, max=100)\n", "model.set_param_hint('center', min=-0.1, max=0.1)\n", "model.set_param_hint('D', min=0.05, max=0.25)\n", "model.set_param_hint('resTime', min=0, max=1)\n", "model.set_param_hint('radius', min=0.9, max=1.1)\n", "model.set_param_hint('DR', min=0, max=1)\n", "\n", "# Fix some of the parameters\n", "model.set_param_hint('q', vary=False)\n", "model.set_param_hint('spectrum_nb', vary=False)\n", "\n", "params = model.make_params()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot of the fitting models without and convoluted with the resolution function\n", "# The values of the parameters are specified below. Therefore they could be \n", "#different from those used in the fitting.\n", "\n", "fig3, ax3 = plt.subplots(1, 2)\n", "# First subplot\n", "for i in range(nb_q_values):\n", " xx = f_5A['x'][i]\n", " ax3[0].plot(xx, \n", " QENSmodels.sqwWaterTeixeira(xx, \n", " data_5A['q'][i], \n", " scale=1, \n", " center=0, \n", " D=1, \n", " resTime=1, \n", " radius=1, \n", " DR=1), \n", " label=f\"q={data_5A['q'][i]:.2}\") \n", "\n", "ax3[0].grid(True)\n", "ax3[0].set(xlabel='Omega', ylabel='S(Q,w)', xlim=(-1, 1), title='No resolution') \n", "ax3[0].tick_params()\n", "\n", "plt.tight_layout(rect=[0, 0, 1, 0.8])\n", "\n", "ax3[0].legend(bbox_to_anchor=(0., 1.1, 2., 0.102), \n", " loc='lower right', \n", " ncol=5, \n", " mode=\"expand\", \n", " borderaxespad=0., \n", " fontsize=8)\n", "\n", "# Second subplot\n", "for i in range(nb_q_values):\n", " params_plot = model.make_params(nb_spectrum=i, \n", " q=data_5A['q'][i],\n", " scale=10.,\n", " center=0.,\n", " D=0.13,\n", " resTime=0.1,\n", " radius=1.,\n", " DR=0.3) \n", " xx = f_5A['x'][i]\n", " ax3[1].plot(xx, model.eval(params_plot, w=xx))\n", "\n", "ax3[1].grid(True)\n", "ax3[1].set(xlabel='w', \n", " ylabel='R $\\otimes$ S(Q,w)', \n", " xlim=(-1,1), \n", " title='Convoluted with resolution')\n", "ax3[1].tick_params();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running the fit using `lmfit`" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ini_values = {'scale': 10., \n", " 'center': 0., \n", " 'D': 0.13, \n", " 'resTime': 0.1, \n", " 'radius': 1., \n", " 'DR': 0.3}\n", "\n", "result_fit = [None,] * nb_q_values # store fits for all spectra\n", "\n", "for i in range(nb_q_values):\n", " params = model.make_params(nb_spectrum=i, \n", " q=data_5A['q'][i],\n", " scale=ini_values['scale'], \n", " center=ini_values['center'], \n", " D=ini_values['D'],\n", " resTime=ini_values['resTime'],\n", " radius=ini_values['radius'],\n", " DR=ini_values['DR'])\n", " \n", " # Q-independent parameters\n", " if i==0:\n", " D_value = params['D'].value\n", " resTime_value = params['resTime'].value\n", " radius_value = params['radius'].value\n", " DR_value = params['DR'].value\n", " else:\n", " params['D'].set(value=D_value)\n", " params['resTime'].set(value=resTime_value)\n", " params['radius'].set(value=radius_value)\n", " params['DR'].set(value=DR_value)\n", "\n", " result_fit[i] = model.fit(f_5A['y'][i], \n", " params,\n", " w=f_5A['x'][i])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Showing the results\n", "\n", "using methods implemented in `lmfit`" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# display result\n", "for i in range(nb_q_values):\n", " print(f'Result of fit {i}:\\n', result_fit[i].fit_report())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot results using lmfit's features" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for i in range(nb_q_values):\n", " result_fit[i].plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Other option to plot experimental data, initial fitting model and fitted model for each spectrum using `matplotlib.pyplot`" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for indx in range(nb_q_values):\n", " fig4 = plt.figure()\n", " plt.plot(f_5A['x'][indx], f_5A['y'][indx], 'bo',label='exp')\n", " plt.plot(f_5A['x'][indx], result_fit[indx].init_fit, 'k--',label='ini')\n", " plt.plot(f_5A['x'][indx], result_fit[indx].best_fit, 'r-', label='fin')\n", " plt.title(f\"q={data_5A['q'][indx]:.2}\")\n", " plt.legend()\n", " plt.grid()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "hide_input": false, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.9" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": true, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "165px" }, "toc_section_display": true, "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }