{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "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", "# Water Teixeira ∗ Resolution with bumps \n", "\n", "## Introduction\n", "\n", "
\n", " \n", "The objective of this notebook is to show how to use the model of \n", "the QENSlibrary corresponding to a \n", "combination of a jump translational diffusion with an\n", "isotropic rotational diffusion to perform some \n", "fits using bumps .\n", "
\n", "\n", "The **reference data** were generated using the above model with the following parameters: \n", "- D = 0.145 Å$^2\\times$meV \n", "- Residence time = 1 meV$^{-1}$ \n", "- Radius = 1.10 Å \n", "- $D_{rot}$ = 0.125 meV \n", " \n", "The model was convoluted with a Gaussian resolution function \n", "of FWHM = 0.1 meV, centered randomly in the range \\[-0.01, +0.01\\] meV.\n", "\n", "Finally the data were sampled randomly from a Poisson distribution.\n", "\n", "The data do not have a background.\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", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Units of parameters for selected QENS model and experimental data\n", "dict_physical_units = {'omega': \"meV\", \n", " 'q': \"1/Angstrom\", \n", " 'D': \"meV.Angstrom^2\", \n", " 'resTime': \"1/meV\",\n", " 'radius': \"Angstrom\",\n", " 'DR': \"meV\",\n", " 'scale': \"unit_of_signal.meV\",\n", " 'center': \"meV\"}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Import required libraries" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "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": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Setting of fitting\n", "\n", "### Load reference data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "path_to_data = './data/'\n", "\n", "# Read the sample\n", "with h5py.File(path_to_data + 'JumpDiffIsoRot_Sample.hdf', 'r') as f:\n", " hw = f['entry1']['data1']['X'][:]\n", " q = f['entry1']['data1']['Y'][:]\n", " unit_w = f['entry1']['data1']['X'].attrs['long_name']\n", " unit_q = f['entry1']['data1']['Y'].attrs['long_name']\n", " sqw = np.transpose(f['entry1']['data1']['DATA'][:])\n", " err = np.transpose(f['entry1']['data1']['errors'][:])\n", "\n", "# Read resolution\n", "with h5py.File(path_to_data + 'JumpDiffIsoRot_Resol.hdf', 'r') as f:\n", " res = np.transpose(f['entry1']['data1']['DATA'][:])\n", "\n", "# Force resolution function to have unit area\n", "for i in range(len(q)):\n", " area = simps(res[:,i], hw)\n", " res[:,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)):\n", " ax[0].semilogy(hw, sqw[:,i], label=f\"q={q[i]:.1f}\")\n", " ax[1].semilogy(hw, res[:,i], label=f\"q={q[i]:.1f}\")\n", "\n", "ax[1].legend(bbox_to_anchor=(1.1,1.3), loc='upper right', shadow=True)\n", "ax[0].set_title('Signal')\n", "ax[0].grid()\n", "\n", "ax[1].set_title('Resolution')\n", "ax[1].set_xlabel(f\"$\\hbar \\omega ({dict_physical_units['center']})$\")\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\"The names and units of `w` (`x`axis) and `q` are: {unit_w[0].decode()} and {unit_q[0].decode()}, respectively.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create fitting model" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# Fitting model \n", "def model_convol(x, q, scale=1, center=0, D=1, resTime=1, radius=1, DR=1, resolution=None):\n", " model = QENSmodels.sqwWaterTeixeira(x, q, scale, center, D, resTime, radius, DR)\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)):\n", " # Bumps fitting model\n", " model_q = bmp.Curve(\n", " model_convol, \n", " hw, \n", " sqw[:,i], \n", " err[:,i], \n", " name=f'q_{q[i]:.2f}',\n", " q=q[i],\n", " scale=1000, \n", " center=0.0, \n", " D=0.1, \n", " resTime=0.5, \n", " radius=1.,\n", " DR=0.1, \n", " resolution=res[:, i]\n", " )\n", " model_q.scale.range(0, 1e5)\n", " model_q.center.range(-0.1, 0.1)\n", " model_q.D.range(0, 1)\n", " model_q.resTime.range(0, 5)\n", " model_q.radius.range(0, 3)\n", " model_q.DR.range(1e-12, 2)\n", " \n", " # Q-independent parameters\n", " if i == 0:\n", " D_q = model_q.D\n", " T_q = model_q.resTime\n", " R_q = model_q.radius\n", " DR_q = model_q.DR\n", " else:\n", " model_q.D = D_q\n", " model_q.resTime = T_q\n", " model_q.radius = R_q\n", " model_q.DR = DR_q\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)-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 ({dict_physical_units['center']})$\")\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": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Choice of minimizer for bumps" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "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": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Setting for running bumps" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "steps_fitting = ipywidgets.IntText(\n", " value=100,\n", " description='Number of steps when fitting',\n", " style={'description_width': 'initial'})\n", "\n", "steps_fitting" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "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": { "slideshow": { "slide_type": "subslide" } }, "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": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Showing the results" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "problem.summarize().splitlines()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "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", " if k in dict_physical_units.keys():\n", " print(k, \":\", format_uncertainty_pm(v, dv), dict_physical_units[k])\n", " else:\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)-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 ({dict_physical_units['center']})$\")\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" }, "livereveal": { "scroll": true }, "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": "656px", "left": "784px", "top": "111.133px", "width": "159.283px" }, "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 }