{ "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's one of the reasons why different minimizers have been used in the notebooks provided as examples. \n", "But, note that the initial guessed parameters might not be optimal, resulting in a poor fit of the reference data.\n", "\n", "
\n", "\n", "\n", "# Lorentzian + Isotropic Rotational diffusion ∗ Resolution with bumps\n", "\n", "## Introduction\n", "\n", "
\n", " \n", "The objective of this notebook is to show how to use a combination of models from the QENS library i.e. Lorentzian and IsotropicRotationalDiffusion models.\n", "
\n", "\n", "The data are a set of water data measured at IN5 (ILL).\n", "\n", "**Reference:** J. Qvist, H. Schober and B. Halle, *J. Chem. Phys.* **134**, 144508 (2011)\n", "\n", "### Physical units\n", "For information about unit conversion, please refer to the jupyter notebook called `Convert_units.ipynb` in the `tools` folder.\n", "\n", "\n", "## Import libraries" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import ipywidgets\n", "import h5py\n", "import QENSmodels\n", "import numpy as np\n", "from scipy.integrate import simps\n", "import bumps.names as bmp\n", "from bumps import fitters\n", "from bumps.formatnum import format_uncertainty_pm\n", "import matplotlib.pyplot as plt\n", "\n", "%matplotlib widget" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setting of fitting\n", "\n", "### Load reference data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "path_to_data = './data/'\n", "\n", "# Data\n", "# Wavelength 5 Angstrom\n", "with h5py.File(path_to_data + 'H2O_293K_5A.hdf', 'r') as f:\n", " hw_5A = f['entry1']['data1']['X'][:]\n", " q_5A = f['entry1']['data1']['Y'][:]\n", " unit_w5A = f['entry1']['data1']['X'].attrs['long_name']\n", " unit_q5A = f['entry1']['data1']['Y'].attrs['long_name']\n", " sqw_5A = np.transpose(f['entry1']['data1']['DATA'][:])\n", " err_5A = np.transpose(f['entry1']['data1']['errors'][:])\n", "\n", "# Resolution\n", "# Wavelength 5 Angstrom\n", "with h5py.File(path_to_data + 'V_273K_5A.hdf', 'r') as f:\n", " res_5A = np.transpose(f['entry1']['data1']['DATA'][:])\n", "\n", "# Force resolution function to have unit area\n", "# Wavelength 5 Angstrom\n", "for i in range(len(q_5A)):\n", " area = simps(res_5A[:, i], hw_5A)\n", " res_5A[:, i] /= area" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(nrows=2, sharex=True)\n", "\n", "for i in range(len(q_5A)):\n", " ax[0].semilogy(hw_5A, sqw_5A[:,i], label=f\"q={q_5A[i]:.1f}\")\n", " ax[1].semilogy(hw_5A, res_5A[:,i], label=f\"q={q_5A[i]:.1f}\")\n", "\n", "ax[0].set_title(r'Signal 5 $\\AA$')\n", "ax[0].grid()\n", "\n", "ax[1].set_title(r'Resolution 5 $\\AA$')\n", "ax[1].set_xlabel(f\"$\\hbar \\omega$\")\n", "ax[1].grid()" ] }, { "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\"At 5 Angstroms, the names and units of `w` (`x`axis) and `q` are: {unit_w5A[0].decode()} and {unit_q5A[0].decode()}, respectively.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create fitting model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Fit range -1 to +1 meV\n", "idx_5A = np.where(np.logical_and(hw_5A > -1.0, hw_5A < 1.0))\n", "\n", "# Fitting model\n", "def model_convol(x, q, scale=1, center=0, hwhm=1, radius=1, DR=1, resolution=None):\n", " model = QENSmodels.lorentzian(\n", " x, \n", " scale, \n", " center, \n", " hwhm\n", " ) + QENSmodels.sqwIsotropicRotationalDiffusion(\n", " x,\n", " q, \n", " scale, \n", " center, \n", " radius, \n", " DR\n", " )\n", " return np.convolve(model, resolution/resolution.sum(), mode='same')\n", "\n", "# Fit\n", "model_all_qs = []\n", "\n", "for i in range(len(q_5A)):\n", "\n", " x = hw_5A[idx_5A]\n", " data = sqw_5A[idx_5A, i]\n", " error = err_5A[idx_5A, i]\n", " resol = res_5A[idx_5A, i]\n", "\n", " # Select only valid data (error = -1 for Q, w points not accessible)\n", " valid = np.where(error > 0.0)\n", " x = x[valid[1]]\n", " data = data[valid]\n", " error = error[valid]\n", " resol = resol[valid]\n", "\n", " # model\n", " model_q = bmp.Curve(\n", " model_convol, \n", " x, data, error,\n", " name=f'q5A_{q_5A[i]:.2f}',\n", " q=q_5A[i], \n", " scale=15, \n", " center=0.0, \n", " hwhm=0.1, \n", " radius=1.1, \n", " DR=1., \n", " resolution=resol\n", " )\n", "\n", " # Fitted parameters\n", " model_q.scale.range(0, 1e2)\n", " model_q.center.range(-0.1, 0.1)\n", " model_q.hwhm.range(0., 1)\n", " model_q.radius.range(0.9, 1.1)\n", " model_q.DR.range(0.01, 5)\n", "\n", " # Q-independent parameters\n", " if i == 0:\n", " hwhm_q = model_q.hwhm\n", " R_q = model_q.radius\n", " DR_q = model_q.DR\n", " else:\n", " model_q.hwhm = hwhm_q\n", " model_q.radius = R_q\n", " model_q.DR = DR_q\n", "\n", " model_all_qs.append(model_q)\n", "\n", "problem = bmp.FitProblem(model_all_qs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Display initial configuration: experimental data, fitting model with initial guesses" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "slider = ipywidgets.IntSlider(value=0, min=0, max=len(q_5A)-1, continuous_update=False)\n", "output = ipywidgets.Output()\n", "\n", "def fig_q(model, ax, q_index=0):\n", " \"\"\"\n", " Plot of experimental data, fitting model and residual for a selected q value\n", " \n", " Parameters\n", " ----------\n", " model: list of bumps.curve.Curves for all q\n", " \n", " ax: matplotlib.axes to be updated when changing the ipywidgets\n", " \n", " q_index: int\n", " index of q to be plotted\n", " \n", " \"\"\"\n", " model = model[q_index]\n", " ax[0].errorbar(model.x,\n", " model.y, \n", " yerr=model.dy,\n", " label='experimental data',\n", " color='C0')\n", " ax[0].plot(model.x,\n", " model.theory(), \n", " label='theory (model)',\n", " color='C1')\n", " ax[0].set_title(f'Model {model.name} - $\\chi^2$={problem.chisq_str()}')\n", " ax[0].legend()\n", " ax[1].plot(model.x, model.residuals(), marker='o', linewidth=0, markersize=3, color='C0')\n", " \n", "\n", "with output:\n", " fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True)\n", " ax[0].grid(); ax[1].grid()\n", " ax[1].set_ylabel('Residual')\n", " ax[1].set_xlabel(f\"$\\hbar \\omega$\")\n", " fig_q(model_all_qs, ax, 0)\n", " \n", " \n", "def update_profile(change):\n", " \"\"\"\n", " Update plots for a new q-value\n", " \"\"\"\n", " with output:\n", " ax[0].clear(); ax[1].lines.clear()\n", " ax[0].grid()\n", " fig_q(model_all_qs, ax,change['new'])\n", " \n", "slider.observe(update_profile, names=\"value\")\n", "\n", "slider_label = ipywidgets.Label(\"q value to display\")\n", "slider_comp = ipywidgets.HBox([slider_label, slider])\n", "ipywidgets.VBox([slider_comp, output])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "problem.summarize().splitlines()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Choice of minimizer for bumps" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "options_dict = {} \n", "\n", "for item in fitters.__dict__.keys():\n", " if item.endswith('Fit') and fitters.__dict__[item].id in fitters.FIT_AVAILABLE_IDS:\n", " options_dict[fitters.__dict__[item].name] = fitters.__dict__[item].id\n", "\n", "w_choice_minimizer = ipywidgets.Dropdown(\n", " options=list(options_dict.keys()),\n", " value='Levenberg-Marquardt',\n", " description='Minimizer:',\n", " layout=ipywidgets.Layout(height='40px'))\n", "\n", "w_choice_minimizer" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Number of steps for running fit using bumps " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "steps_fitting = ipywidgets.IntText(\n", " value=100,\n", " step=100,\n", " description='Number of steps when fitting',\n", " style={'description_width': 'initial'})\n", "\n", "steps_fitting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running the fit\n", "\n", "Run the fit using the *minimizer* defined above with a number of *steps* also specified above." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Preview of the settings\n", "print('Initial chisq', problem.chisq_str())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def settings_selected_optimizer(chosen_minimizer):\n", " \"\"\" \n", " List the settings available for the selected optimizer\n", " \n", " This list can be used as arguments for the `fit` function\n", " \"\"\"\n", " \n", " assert type(chosen_minimizer) == ipywidgets.widgets.widget_selection.Dropdown\n", " \n", " for item in fitters.__dict__.keys():\n", " if item.endswith('Fit') and \\\n", " fitters.__dict__[item].id == options_dict[chosen_minimizer.value]:\n", " return [elt[0] for elt in fitters.__dict__[item].settings]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print((f\"With {w_choice_minimizer.value} optimizer, \"\n", " f\"you can use {settings_selected_optimizer(w_choice_minimizer)} as arguments of `fit`\"))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "result = fitters.fit(\n", " problem,\n", " starts=10,\n", " keep_best=True,\n", " method=options_dict[w_choice_minimizer.value], \n", " steps=int(steps_fitting.value)\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Showing the results" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "problem.summarize().splitlines()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Other method to display the results of the fit (chi**2 and parameters' values)\n", "print(\"final chisq\", problem.chisq_str())\n", "for k, v, dv in zip(problem.labels(), result.x, result.dx):\n", " print(k, \":\", format_uncertainty_pm(v, dv))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Display final configuration: experimental data, fitting model with output of fitting for the refined parameters" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "slider1 = ipywidgets.IntSlider(value=0, min=0, max=len(q_5A)-1, continuous_update=False)\n", "output1 = ipywidgets.Output()\n", "\n", "with output1:\n", " fig1, ax1 = plt.subplots(nrows=2, ncols=1, sharex=True)\n", " ax1[0].grid(); ax1[1].grid()\n", " ax1[1].set_ylabel('Residual')\n", " ax1[1].set_xlabel(f\"$\\hbar \\omega$\")\n", " fig_q(model_all_qs, ax1, 0)\n", " \n", "def update_profile1(change):\n", " \"\"\"\n", " Update plots for a new q-value\n", " \"\"\"\n", " with output1:\n", " ax1[0].clear(); ax1[1].lines.clear()\n", " ax1[0].grid()\n", " fig_q(model_all_qs, ax1, change['new'])\n", " \n", "slider1.observe(update_profile1, names=\"value\")\n", "\n", "slider1_label = ipywidgets.Label(\"q value to display\")\n", "slider1_comp = ipywidgets.HBox([slider1_label, slider1])\n", "ipywidgets.VBox([slider1_comp, output1])" ] }, { "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": {}, "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 }