diff --git a/notebooks/.ipynb_checkpoints/water_fitting-checkpoint.ipynb b/notebooks/.ipynb_checkpoints/water_fitting-checkpoint.ipynb
new file mode 100644
index 0000000..5a1c5b2
--- /dev/null
+++ b/notebooks/.ipynb_checkpoints/water_fitting-checkpoint.ipynb
@@ -0,0 +1,668 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Using LMFIT for Sequential and Global Fit of a Water Sample\n",
+ "\n",
+ "QENS data obtained from a sample of water at room temperature.\n",
+ "\n",
+ "Steps shown in this tutorial:\n",
+ "\n",
+ "- Create a simple model: one elastic line plus one lorentzian plus a linear background.\n",
+ "- Interactively find an initial guess of the model parameters for the spectrum with lowest Q.\n",
+ "- Automatic fit the spectrum with lowest Q and visualize results.\n",
+ "- Automatic sequential fit of the remaining spectra and visualize results.\n",
+ "- Fit the Q-dependence of the Lorentzian's HWHM to a Teixeira model.\n",
+ "- Simultaneous fit of all spectra with two global parameters: water diffusion coefficient and jump time.\n",
+ "- Visualize results from the simultaneous fit.\n",
+ "\n",
+ "### Useful links\n",
+ "- [lmfit documentation](https://lmfit.github.io/lmfit-py/index.html) Curve fitting\n",
+ "- [qef documentation](http://qef.readthedocs.io/en/latest/) Utilities for QENS fitting\n",
+ "- [matplotlib](https://matplotlib.org) Plotting with python\n",
+ "- [Post your questions](https://gitter.im/basisdoc/Lobby)\n",
+ "\n",
+ "
Table of Contents
\n",
+ "- Load data and visualize data \n",
+ "- Define the Fitting Range \n",
+ "- Define the model \n",
+ "- Obtain an initial guess \n",
+ "- Carry out the fit and look at results \n",
+ "- Sequential Fit \n",
+ "- Visualize sequential fits \n",
+ "- Q-dependence of some parameters \n",
+ "- Initial Guess for Teixeira Water model \n",
+ "- Model for Simultaneous Fit of All Spectra with Teixeira Water Model \n",
+ "- Visualize the Simultaneous Fit "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Imports for fitting"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from __future__ import (absolute_import, division, print_function)\n",
+ "\n",
+ "import os\n",
+ "from os.path import join as pjn\n",
+ "import sys\n",
+ "import functools\n",
+ "import lmfit\n",
+ "from lmfit.models import LinearModel, LorentzianModel, ConstantModel, LinearModel\n",
+ "\n",
+ "import qef\n",
+ "from qef.io.loaders import load_nexus\n",
+ "from qef.models.deltadirac import DeltaDiracModel\n",
+ "from qef.models.tabulatedmodel import TabulatedModel\n",
+ "from qef.models.resolution import TabulatedResolutionModel\n",
+ "from qef.models.teixeira import TeixeiraWaterModel\n",
+ "from qef.operators.convolve import Convolve\n",
+ "\n",
+ "# Finding out test data\n",
+ "qef_module_file = os.path.abspath(qef.__file__)\n",
+ "qef_dir = os.path.join(os.path.dirname(qef_module_file), os.pardir)\n",
+ "data_dir = os.path.join(qef_dir, 'tests', 'data', 'io')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from IPython.core.display import HTML\n",
+ "from IPython.core.display import display\n",
+ "from ipywidgets import widgets\n",
+ "\n",
+ "import numpy as np\n",
+ "\n",
+ "import matplotlib.pyplot as plt\n",
+ "%matplotlib notebook"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Imports for plotting and widgets"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "(Top)
\n",
+ "\n",
+ "Starting from the first spectrum, we iteratively fit spectra of higher q's.\n",
+ "\n",
+ "*Special instructions:* If the previous guess was not for the first spectrum, the code will not run but raise and error. Go back to the previous cell and carry out a guess for the first spectrum, then come back here.\n",
+ "\n",
+ "We do not assume any particular Q-dependence for the width of the Lorentzian function."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_spectra = len(fr['y'])\n",
+ "fits = [None,] * n_spectra # store fits for all the tried spectra\n",
+ "fits[0] = fit # store previous fit\n",
+ "for i in range(1, n_spectra):\n",
+ " y_exp = fr['y'][i]\n",
+ " e_exp = fr['e'][i]\n",
+ " fit = model.fit(y_exp, x=e_vals, params=params, weights = 1.0 / e_exp)\n",
+ " params = fit.params # update params with results from previous spectrum's fit\n",
+ " fits[i] = fit # store fit results\n",
+ "\n",
+ "# Show Chi-square versus Q\n",
+ "chi2s = [fit.redchi for fit in fits]\n",
+ "f, ax = plt.subplots()\n",
+ "ax.plot(qs, [fit.redchi for fit in fits])\n",
+ "ax.set_xlabel('Q')\n",
+ "ax.set_ylabel('Chi-squared')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "(Top)
\n",
+ "\n",
+ "The sample is liquid water, thus we expect $EISF \\ll 1$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "names = ('l_fwhm', 'e_amplitude') # fitting parameters we want to plot\n",
+ "ylabels = ('FWHM', 'EISF') # labels on the Y-axis of the plots\n",
+ "xlabels = ('Q^2', 'Q') # labels on the X-axis of the plots\n",
+ "\n",
+ "q_vals = np.asarray(qs)\n",
+ "xs = (q_vals * q_vals, q_vals) # we want to plot FWHM versus Q^2 and EISF versus Q\n",
+ "\n",
+ "f, axs = plt.subplots(1, len(names)) # as many plots as fitting parameters of interest\n",
+ "for i in range(len(names)):\n",
+ " name = names[i] # name of the fit parameter\n",
+ " y = [fit.params[name].value for fit in fits]\n",
+ " ax = axs[i] # select appropriate plotting area\n",
+ " ax.plot(xs[i], y, marker='o', linestyle='dashed')\n",
+ " ax.set_xlabel(xlabels[i])\n",
+ " ax.set_ylabel(ylabels[i])\n",
+ "\n",
+ "plt.tight_layout()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "(Top)
Initial Guess for Teixeira Water model
\n",
+ "\n",
+ "We use the previous $FWHM$ to fit $HWHM(Q^2)$ to Teixeira's water model to obtain initial diffusion $D$ and relaxation time coefficients $\\tau$\n",
+ "\n",
+ "
\n",
+ "\n",
+ "If $Q$ in Angstroms, $HHWM$ in $meV$, and $\\hbar$ in $meV \\cdot ps$, then units of $D$ are $A^2/ps$ and units of $\\tau$ are $ps$.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Collect FHWM from the fits, and associated error in estimation of these optimal values.\n",
+ "hwhms = 0.5 * np.asarray([fit.params['l_fwhm'].value for fit in fits]) # HWHM values\n",
+ "\n",
+ "\n",
+ "# Create the model\n",
+ "from qef.constants import hbar # units of meV x ps or ueV x ns\n",
+ "from lmfit.model import Model\n",
+ "\n",
+ "def teixeira(q2s, difcoef, tau):\n",
+ " r\"\"\"Calculate HWHM for a given Q, diffusion coefficient, and relaxation time\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " q2s : float\n",
+ " Q^2 values\n",
+ " difcoef : float\n",
+ " Diffusion coefficient parameter\n",
+ " tau : float\n",
+ " Relaxation time parameter\n",
+ "\n",
+ " Returns\n",
+ " -------\n",
+ " numpy.ndarray\n",
+ " HWHM values\n",
+ " \"\"\"\n",
+ " dq2 = difcoef * q2s\n",
+ " return hbar * dq2 / (1 + dq2 * tau)\n",
+ "\n",
+ "teixeira_model = Model(teixeira) # create LMFIT Model instance\n",
+ "teixeira_model.set_param_hint('difcoef', min=0) # diffusion coefficient must be positive\n",
+ "teixeira_model.set_param_hint('tau', min=0) # relaxation coefficient must be positive\n",
+ "\n",
+ "\n",
+ "# Carry out the fit\n",
+ "\n",
+ "teixeira_params = teixeira_model.make_params(difcoef=1.0, tau=1.0) # initial guess\n",
+ "teixeira_fit = teixeira_model.fit(hwhms, q2s=np.square(q_vals), params=teixeira_params)\n",
+ "\n",
+ "# Visualize fit results\n",
+ "o_p = teixeira_fit.params # optimal parameters\n",
+ "fmt = 'Chi-square = {}\\nD = {} A^2/ps\\ntau = {} ps'\n",
+ "print(fmt.format(teixeira_fit.redchi, o_p['difcoef'].value, o_p['tau'].value))\n",
+ "teixeira_fit.plot(xlabel='Q^2', ylabel='HWHM')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "(Top)
Model for Simultaneous Fit of All Spectra with Teixeira Water Model
\n",
+ "\n",
+ "We impose a Q-dependende for the FWHM of the Lorentzian, given by the Teixeira water model. Parameters $D$ and $\\tau$ are the only parameters that are same for all spectra.\n",
+ "\n",
+ "
\n",
+ "\n",
+ "We use $D$ and $\\tau$ of the previous fit as initial guesses. We use the sequential fits we did before to initialize all other parameters."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Create a model for each spectrum\n",
+ "\n",
+ "#initialize models and parameter sets for each spectrum\n",
+ "\n",
+ "\n",
+ "# create one model for each spectrum, but collect all parameters under\n",
+ "# a single instance of the Parameters class.\n",
+ "l_model = list()\n",
+ "g_params = lmfit.Parameters()\n",
+ "for i in range(n_spectra):\n",
+ " m, ps = generate_model_and_params(spectrum_index=i) # model and parameters for one of the spectra\n",
+ " l_model.append(m)\n",
+ " [g_params.add(p) for p in ps.values()]\n",
+ "\n",
+ "# Initialize parameter set with the optimized parameters from the sequential fit\n",
+ "for i in range(n_spectra):\n",
+ " optimized_params = fits[i].params # these are I_c, e_amplitude,...\n",
+ " for name in optimized_params:\n",
+ " prefix, base = name.split('_') # for instance, 'e_amplitude' splitted into 'e', and 'amplitude'\n",
+ " i_name = prefix + '_{}_'.format(i) + base # i_name is 'e_3_amplitude' for i=3\n",
+ " g_params[i_name].set(value=optimized_params[name].value)\n",
+ "\n",
+ "# Introduce global parameters diff and tau. Use previous optimized values as initial guess\n",
+ "g_params.add('difcoef', value=o_p['difcoef'].value, min=0)\n",
+ "g_params.add('tau', value=o_p['tau'].value, min=0)\n",
+ "\n",
+ "# Tie each lorentzian l_i_sigma to the teixeira expression\n",
+ "for i in range(n_spectra):\n",
+ " q2 = q_vals[i] * q_vals[i]\n",
+ " teixeira_expression = '{hbar}*difcoef*{q2}/(1+difcoef*{q2}*tau)'\n",
+ " g_params['l_{}_sigma'.format(i)].set(expr=teixeira_expression.format(hbar=hbar, q2=q2))\n",
+ "\n",
+ "print('Number of varying parameters =',len([p for p in g_params.values() if p.vary]),'!')\n",
+ "g_params.pretty_print()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "(Top)
Carry out the Simultaneous Fit
"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def residuals(params):\n",
+ " r\"\"\"Difference between model and experiment, weighted by the experimental error\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " params : lmfit.Parameters\n",
+ " Parameters for the global model\n",
+ "\n",
+ " Returns\n",
+ " -------\n",
+ " numpy.ndarray\n",
+ " 1D array of residuals for the global model\n",
+ " \"\"\"\n",
+ " l_residuals = list()\n",
+ " for i in range(n_spectra):\n",
+ " x = fr['x'] # fitting range of energies\n",
+ " y = fr['y'][i] # associated experimental intensities\n",
+ " e = fr['e'][i] # associated experimental errors\n",
+ " model_evaluation = l_model[i].eval(x=x, params=params)\n",
+ " l_residuals.append((model_evaluation - y) / e)\n",
+ " return np.concatenate(l_residuals)\n",
+ "\n",
+ "# Minimizer object using the parameter set for all models and the\n",
+ "# function to calculate all the residuals.\n",
+ "minimizer = lmfit.Minimizer(residuals, g_params)\n",
+ "g_fit = minimizer.minimize()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print('Chi-square = {:.2f}\\n'.format(g_fit.redchi))\n",
+ "fmt = 'D = {:.3f} A^2/ps, tau = {:.3f} ps'\n",
+ "#print('Before: ', fmt.format(g_fit.init_values['difcoef'], g_fit.init_values['tau']))\n",
+ "print('Before: D = 0.156 A^2/ps, tau = 1.112 ps')\n",
+ "print('After: ', fmt.format(g_fit.params['difcoef'].value, g_fit.params['tau'].value))\n",
+ "print('Teixeira: D = 0.19 A^2/ps, tau = 1.25 ps (J. Teixeira et al., Phys. Rev. A, 31(3), 1913-947 (1985))')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "(Top)
Visualize the Simultaneous Fit
"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "g_fig = plt.figure()\n",
+ "g_gs = plt.GridSpec(nrows=2, ncols=1, height_ratios=[4, 1])\n",
+ "ax_fit = g_fig.add_subplot(g_gs[0]) # host the fit\n",
+ "ax_res = g_fig.add_subplot(g_gs[1], sharex=ax_fit) # host the residuals\n",
+ "\n",
+ "def g_fit_changed(bunch):\n",
+ " i = bunch['new']\n",
+ " ax_fit.clear()\n",
+ " ax_fit.errorbar(fr['x'], fr['y'][i], yerr=fr['e'][i], color='black',\n",
+ " marker='o', markersize=1, markerfacecolor='none', label='data',\n",
+ " linestyle='none')\n",
+ " model_evaluation = l_model[i].eval(x=fr['x'], params=g_params)\n",
+ " ax_fit.plot(fr['x'], model_evaluation, color='red', linewidth=4, label='best fit')\n",
+ " ax_fit.legend()\n",
+ " ax_res.clear()\n",
+ " ax_res.set_xlabel('Energy (meV)')\n",
+ " ax_res.plot(fr['x'], model_evaluation - fr['y'][i], color='black', label='residuals')\n",
+ " ax_res.legend()\n",
+ "\n",
+ "# Widget for the spectrum index\n",
+ "gv_label = widgets.Label('spectrum index', layout=widgets.Layout(width='10%'))\n",
+ "gv_int_text = widgets.BoundedIntText(value=0, min=0, max=(n_spectra - 1),\n",
+ " layout=widgets.Layout(width='20%'))\n",
+ "gv_int_text.observe(g_fit_changed, 'value')\n",
+ "gv_hbox_l = [widgets.HBox([gv_label, gv_int_text])]\n",
+ "\n",
+ "gv_vertical_layout = widgets.VBox(gv_hbox_l)\n",
+ "display(gv_vertical_layout)"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 2",
+ "language": "python",
+ "name": "python2"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 2
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython2",
+ "version": "2.7.12"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/notebooks/.ipynb_checkpoints/widget_builder_on_the_fly-checkpoint.ipynb b/notebooks/.ipynb_checkpoints/widget_builder_on_the_fly-checkpoint.ipynb
new file mode 100644
index 0000000..4c2c3d1
--- /dev/null
+++ b/notebooks/.ipynb_checkpoints/widget_builder_on_the_fly-checkpoint.ipynb
@@ -0,0 +1,972 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from IPython.core.display import HTML\n",
+ "from IPython.core.display import display\n",
+ "from ipywidgets import widgets\n",
+ "\n",
+ "import numpy as np\n",
+ "\n",
+ "import matplotlib.pyplot as plt\n",
+ "%matplotlib notebook"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "3d9d64fbf07845d4891d88ad42c95d91",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "A Jupyter Widget"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "user_input_ui = widgets.HBox([widgets.Label(\"Define Function:\",\n",
+ " layout=widgets.Layout(width='10%')),\n",
+ " widgets.Text(value='',\n",
+ " placeholder=\"Define function here\",\n",
+ " layout=widgets.Layout(width='60%')),\n",
+ " widgets.Label(\"Ex: f(a,b,c)=a+b*c\",\n",
+ " layout=widgets.Layout(width='20%'))])\n",
+ "display(user_input_ui)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "..."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "At this point, you retrieve the name of the coefficient defined in the function, let's pretend those are \n",
+ "\n",
+ " * argument_1\n",
+ " * argument_2\n",
+ " * argument_3\n",
+ " \n",
+ " ```\n",
+ " function(x, argument_1, argument_2, argument_3) = 3*argument_1*x + argument_2 * argument_3\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "list_arguments = ['argument_1', 'argument_2', 'argument_3']"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def function(x, argument_1, argument_2, argument_3):\n",
+ " return 3*argument_1*x*x + argument_2*x + argument_3"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def value_changed(value):\n",
+ " # retrieve value of all widgets\n",
+ " list_value = []\n",
+ " for _top_widget in vertical_layout.children:\n",
+ " _text_float = _top_widget.children[1]\n",
+ " list_value.append(_text_float.value)\n",
+ "\n",
+ " x = np.linspace(0,100)\n",
+ " plt.clf()\n",
+ " plt.plot(x, function(x, list_value[0], list_value[1], list_value[2]))\n",
+ " plt.title(\"argument_1: {}, argument_2: {}, argument_3: {}\".format(list_value[0], list_value[1], list_value[2]))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "9adebb8a047744ca8a4fc9aa396ce817",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "A Jupyter Widget"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/javascript": [
+ "/* Put everything inside the global mpl namespace */\n",
+ "window.mpl = {};\n",
+ "\n",
+ "\n",
+ "mpl.get_websocket_type = function() {\n",
+ " if (typeof(WebSocket) !== 'undefined') {\n",
+ " return WebSocket;\n",
+ " } else if (typeof(MozWebSocket) !== 'undefined') {\n",
+ " return MozWebSocket;\n",
+ " } else {\n",
+ " alert('Your browser does not have WebSocket support.' +\n",
+ " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n",
+ " 'Firefox 4 and 5 are also supported but you ' +\n",
+ " 'have to enable WebSockets in about:config.');\n",
+ " };\n",
+ "}\n",
+ "\n",
+ "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n",
+ " this.id = figure_id;\n",
+ "\n",
+ " this.ws = websocket;\n",
+ "\n",
+ " this.supports_binary = (this.ws.binaryType != undefined);\n",
+ "\n",
+ " if (!this.supports_binary) {\n",
+ " var warnings = document.getElementById(\"mpl-warnings\");\n",
+ " if (warnings) {\n",
+ " warnings.style.display = 'block';\n",
+ " warnings.textContent = (\n",
+ " \"This browser does not support binary websocket messages. \" +\n",
+ " \"Performance may be slow.\");\n",
+ " }\n",
+ " }\n",
+ "\n",
+ " this.imageObj = new Image();\n",
+ "\n",
+ " this.context = undefined;\n",
+ " this.message = undefined;\n",
+ " this.canvas = undefined;\n",
+ " this.rubberband_canvas = undefined;\n",
+ " this.rubberband_context = undefined;\n",
+ " this.format_dropdown = undefined;\n",
+ "\n",
+ " this.image_mode = 'full';\n",
+ "\n",
+ " this.root = $('');\n",
+ " this._root_extra_style(this.root)\n",
+ " this.root.attr('style', 'display: inline-block');\n",
+ "\n",
+ " $(parent_element).append(this.root);\n",
+ "\n",
+ " this._init_header(this);\n",
+ " this._init_canvas(this);\n",
+ " this._init_toolbar(this);\n",
+ "\n",
+ " var fig = this;\n",
+ "\n",
+ " this.waiting = false;\n",
+ "\n",
+ " this.ws.onopen = function () {\n",
+ " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n",
+ " fig.send_message(\"send_image_mode\", {});\n",
+ " if (mpl.ratio != 1) {\n",
+ " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n",
+ " }\n",
+ " fig.send_message(\"refresh\", {});\n",
+ " }\n",
+ "\n",
+ " this.imageObj.onload = function() {\n",
+ " if (fig.image_mode == 'full') {\n",
+ " // Full images could contain transparency (where diff images\n",
+ " // almost always do), so we need to clear the canvas so that\n",
+ " // there is no ghosting.\n",
+ " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n",
+ " }\n",
+ " fig.context.drawImage(fig.imageObj, 0, 0);\n",
+ " };\n",
+ "\n",
+ " this.imageObj.onunload = function() {\n",
+ " fig.ws.close();\n",
+ " }\n",
+ "\n",
+ " this.ws.onmessage = this._make_on_message_function(this);\n",
+ "\n",
+ " this.ondownload = ondownload;\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype._init_header = function() {\n",
+ " var titlebar = $(\n",
+ " '');\n",
+ " var titletext = $(\n",
+ " '');\n",
+ " titlebar.append(titletext)\n",
+ " this.root.append(titlebar);\n",
+ " this.header = titletext[0];\n",
+ "}\n",
+ "\n",
+ "\n",
+ "\n",
+ "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n",
+ "\n",
+ "}\n",
+ "\n",
+ "\n",
+ "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n",
+ "\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype._init_canvas = function() {\n",
+ " var fig = this;\n",
+ "\n",
+ " var canvas_div = $('');\n",
+ "\n",
+ " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n",
+ "\n",
+ " function canvas_keyboard_event(event) {\n",
+ " return fig.key_event(event, event['data']);\n",
+ " }\n",
+ "\n",
+ " canvas_div.keydown('key_press', canvas_keyboard_event);\n",
+ " canvas_div.keyup('key_release', canvas_keyboard_event);\n",
+ " this.canvas_div = canvas_div\n",
+ " this._canvas_extra_style(canvas_div)\n",
+ " this.root.append(canvas_div);\n",
+ "\n",
+ " var canvas = $('');\n",
+ " canvas.addClass('mpl-canvas');\n",
+ " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n",
+ "\n",
+ " this.canvas = canvas[0];\n",
+ " this.context = canvas[0].getContext(\"2d\");\n",
+ "\n",
+ " var backingStore = this.context.backingStorePixelRatio ||\n",
+ "\tthis.context.webkitBackingStorePixelRatio ||\n",
+ "\tthis.context.mozBackingStorePixelRatio ||\n",
+ "\tthis.context.msBackingStorePixelRatio ||\n",
+ "\tthis.context.oBackingStorePixelRatio ||\n",
+ "\tthis.context.backingStorePixelRatio || 1;\n",
+ "\n",
+ " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n",
+ "\n",
+ " var rubberband = $('');\n",
+ " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n",
+ "\n",
+ " var pass_mouse_events = true;\n",
+ "\n",
+ " canvas_div.resizable({\n",
+ " start: function(event, ui) {\n",
+ " pass_mouse_events = false;\n",
+ " },\n",
+ " resize: function(event, ui) {\n",
+ " fig.request_resize(ui.size.width, ui.size.height);\n",
+ " },\n",
+ " stop: function(event, ui) {\n",
+ " pass_mouse_events = true;\n",
+ " fig.request_resize(ui.size.width, ui.size.height);\n",
+ " },\n",
+ " });\n",
+ "\n",
+ " function mouse_event_fn(event) {\n",
+ " if (pass_mouse_events)\n",
+ " return fig.mouse_event(event, event['data']);\n",
+ " }\n",
+ "\n",
+ " rubberband.mousedown('button_press', mouse_event_fn);\n",
+ " rubberband.mouseup('button_release', mouse_event_fn);\n",
+ " // Throttle sequential mouse events to 1 every 20ms.\n",
+ " rubberband.mousemove('motion_notify', mouse_event_fn);\n",
+ "\n",
+ " rubberband.mouseenter('figure_enter', mouse_event_fn);\n",
+ " rubberband.mouseleave('figure_leave', mouse_event_fn);\n",
+ "\n",
+ " canvas_div.on(\"wheel\", function (event) {\n",
+ " event = event.originalEvent;\n",
+ " event['data'] = 'scroll'\n",
+ " if (event.deltaY < 0) {\n",
+ " event.step = 1;\n",
+ " } else {\n",
+ " event.step = -1;\n",
+ " }\n",
+ " mouse_event_fn(event);\n",
+ " });\n",
+ "\n",
+ " canvas_div.append(canvas);\n",
+ " canvas_div.append(rubberband);\n",
+ "\n",
+ " this.rubberband = rubberband;\n",
+ " this.rubberband_canvas = rubberband[0];\n",
+ " this.rubberband_context = rubberband[0].getContext(\"2d\");\n",
+ " this.rubberband_context.strokeStyle = \"#000000\";\n",
+ "\n",
+ " this._resize_canvas = function(width, height) {\n",
+ " // Keep the size of the canvas, canvas container, and rubber band\n",
+ " // canvas in synch.\n",
+ " canvas_div.css('width', width)\n",
+ " canvas_div.css('height', height)\n",
+ "\n",
+ " canvas.attr('width', width * mpl.ratio);\n",
+ " canvas.attr('height', height * mpl.ratio);\n",
+ " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n",
+ "\n",
+ " rubberband.attr('width', width);\n",
+ " rubberband.attr('height', height);\n",
+ " }\n",
+ "\n",
+ " // Set the figure to an initial 600x600px, this will subsequently be updated\n",
+ " // upon first draw.\n",
+ " this._resize_canvas(600, 600);\n",
+ "\n",
+ " // Disable right mouse context menu.\n",
+ " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n",
+ " return false;\n",
+ " });\n",
+ "\n",
+ " function set_focus () {\n",
+ " canvas.focus();\n",
+ " canvas_div.focus();\n",
+ " }\n",
+ "\n",
+ " window.setTimeout(set_focus, 100);\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype._init_toolbar = function() {\n",
+ " var fig = this;\n",
+ "\n",
+ " var nav_element = $('')\n",
+ " nav_element.attr('style', 'width: 100%');\n",
+ " this.root.append(nav_element);\n",
+ "\n",
+ " // Define a callback function for later on.\n",
+ " function toolbar_event(event) {\n",
+ " return fig.toolbar_button_onclick(event['data']);\n",
+ " }\n",
+ " function toolbar_mouse_event(event) {\n",
+ " return fig.toolbar_button_onmouseover(event['data']);\n",
+ " }\n",
+ "\n",
+ " for(var toolbar_ind in mpl.toolbar_items) {\n",
+ " var name = mpl.toolbar_items[toolbar_ind][0];\n",
+ " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
+ " var image = mpl.toolbar_items[toolbar_ind][2];\n",
+ " var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
+ "\n",
+ " if (!name) {\n",
+ " // put a spacer in here.\n",
+ " continue;\n",
+ " }\n",
+ " var button = $('');\n",
+ " button.addClass('ui-button ui-widget ui-state-default ui-corner-all ' +\n",
+ " 'ui-button-icon-only');\n",
+ " button.attr('role', 'button');\n",
+ " button.attr('aria-disabled', 'false');\n",
+ " button.click(method_name, toolbar_event);\n",
+ " button.mouseover(tooltip, toolbar_mouse_event);\n",
+ "\n",
+ " var icon_img = $('');\n",
+ " icon_img.addClass('ui-button-icon-primary ui-icon');\n",
+ " icon_img.addClass(image);\n",
+ " icon_img.addClass('ui-corner-all');\n",
+ "\n",
+ " var tooltip_span = $('');\n",
+ " tooltip_span.addClass('ui-button-text');\n",
+ " tooltip_span.html(tooltip);\n",
+ "\n",
+ " button.append(icon_img);\n",
+ " button.append(tooltip_span);\n",
+ "\n",
+ " nav_element.append(button);\n",
+ " }\n",
+ "\n",
+ " var fmt_picker_span = $('');\n",
+ "\n",
+ " var fmt_picker = $('');\n",
+ " fmt_picker.addClass('mpl-toolbar-option ui-widget ui-widget-content');\n",
+ " fmt_picker_span.append(fmt_picker);\n",
+ " nav_element.append(fmt_picker_span);\n",
+ " this.format_dropdown = fmt_picker[0];\n",
+ "\n",
+ " for (var ind in mpl.extensions) {\n",
+ " var fmt = mpl.extensions[ind];\n",
+ " var option = $(\n",
+ " '', {selected: fmt === mpl.default_extension}).html(fmt);\n",
+ " fmt_picker.append(option)\n",
+ " }\n",
+ "\n",
+ " // Add hover states to the ui-buttons\n",
+ " $( \".ui-button\" ).hover(\n",
+ " function() { $(this).addClass(\"ui-state-hover\");},\n",
+ " function() { $(this).removeClass(\"ui-state-hover\");}\n",
+ " );\n",
+ "\n",
+ " var status_bar = $('');\n",
+ " nav_element.append(status_bar);\n",
+ " this.message = status_bar[0];\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.request_resize = function(x_pixels, y_pixels) {\n",
+ " // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n",
+ " // which will in turn request a refresh of the image.\n",
+ " this.send_message('resize', {'width': x_pixels, 'height': y_pixels});\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.send_message = function(type, properties) {\n",
+ " properties['type'] = type;\n",
+ " properties['figure_id'] = this.id;\n",
+ " this.ws.send(JSON.stringify(properties));\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.send_draw_message = function() {\n",
+ " if (!this.waiting) {\n",
+ " this.waiting = true;\n",
+ " this.ws.send(JSON.stringify({type: \"draw\", figure_id: this.id}));\n",
+ " }\n",
+ "}\n",
+ "\n",
+ "\n",
+ "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
+ " var format_dropdown = fig.format_dropdown;\n",
+ " var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n",
+ " fig.ondownload(fig, format);\n",
+ "}\n",
+ "\n",
+ "\n",
+ "mpl.figure.prototype.handle_resize = function(fig, msg) {\n",
+ " var size = msg['size'];\n",
+ " if (size[0] != fig.canvas.width || size[1] != fig.canvas.height) {\n",
+ " fig._resize_canvas(size[0], size[1]);\n",
+ " fig.send_message(\"refresh\", {});\n",
+ " };\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.handle_rubberband = function(fig, msg) {\n",
+ " var x0 = msg['x0'] / mpl.ratio;\n",
+ " var y0 = (fig.canvas.height - msg['y0']) / mpl.ratio;\n",
+ " var x1 = msg['x1'] / mpl.ratio;\n",
+ " var y1 = (fig.canvas.height - msg['y1']) / mpl.ratio;\n",
+ " x0 = Math.floor(x0) + 0.5;\n",
+ " y0 = Math.floor(y0) + 0.5;\n",
+ " x1 = Math.floor(x1) + 0.5;\n",
+ " y1 = Math.floor(y1) + 0.5;\n",
+ " var min_x = Math.min(x0, x1);\n",
+ " var min_y = Math.min(y0, y1);\n",
+ " var width = Math.abs(x1 - x0);\n",
+ " var height = Math.abs(y1 - y0);\n",
+ "\n",
+ " fig.rubberband_context.clearRect(\n",
+ " 0, 0, fig.canvas.width, fig.canvas.height);\n",
+ "\n",
+ " fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.handle_figure_label = function(fig, msg) {\n",
+ " // Updates the figure title.\n",
+ " fig.header.textContent = msg['label'];\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.handle_cursor = function(fig, msg) {\n",
+ " var cursor = msg['cursor'];\n",
+ " switch(cursor)\n",
+ " {\n",
+ " case 0:\n",
+ " cursor = 'pointer';\n",
+ " break;\n",
+ " case 1:\n",
+ " cursor = 'default';\n",
+ " break;\n",
+ " case 2:\n",
+ " cursor = 'crosshair';\n",
+ " break;\n",
+ " case 3:\n",
+ " cursor = 'move';\n",
+ " break;\n",
+ " }\n",
+ " fig.rubberband_canvas.style.cursor = cursor;\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.handle_message = function(fig, msg) {\n",
+ " fig.message.textContent = msg['message'];\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.handle_draw = function(fig, msg) {\n",
+ " // Request the server to send over a new figure.\n",
+ " fig.send_draw_message();\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.handle_image_mode = function(fig, msg) {\n",
+ " fig.image_mode = msg['mode'];\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.updated_canvas_event = function() {\n",
+ " // Called whenever the canvas gets updated.\n",
+ " this.send_message(\"ack\", {});\n",
+ "}\n",
+ "\n",
+ "// A function to construct a web socket function for onmessage handling.\n",
+ "// Called in the figure constructor.\n",
+ "mpl.figure.prototype._make_on_message_function = function(fig) {\n",
+ " return function socket_on_message(evt) {\n",
+ " if (evt.data instanceof Blob) {\n",
+ " /* FIXME: We get \"Resource interpreted as Image but\n",
+ " * transferred with MIME type text/plain:\" errors on\n",
+ " * Chrome. But how to set the MIME type? It doesn't seem\n",
+ " * to be part of the websocket stream */\n",
+ " evt.data.type = \"image/png\";\n",
+ "\n",
+ " /* Free the memory for the previous frames */\n",
+ " if (fig.imageObj.src) {\n",
+ " (window.URL || window.webkitURL).revokeObjectURL(\n",
+ " fig.imageObj.src);\n",
+ " }\n",
+ "\n",
+ " fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n",
+ " evt.data);\n",
+ " fig.updated_canvas_event();\n",
+ " fig.waiting = false;\n",
+ " return;\n",
+ " }\n",
+ " else if (typeof evt.data === 'string' && evt.data.slice(0, 21) == \"data:image/png;base64\") {\n",
+ " fig.imageObj.src = evt.data;\n",
+ " fig.updated_canvas_event();\n",
+ " fig.waiting = false;\n",
+ " return;\n",
+ " }\n",
+ "\n",
+ " var msg = JSON.parse(evt.data);\n",
+ " var msg_type = msg['type'];\n",
+ "\n",
+ " // Call the \"handle_{type}\" callback, which takes\n",
+ " // the figure and JSON message as its only arguments.\n",
+ " try {\n",
+ " var callback = fig[\"handle_\" + msg_type];\n",
+ " } catch (e) {\n",
+ " console.log(\"No handler for the '\" + msg_type + \"' message type: \", msg);\n",
+ " return;\n",
+ " }\n",
+ "\n",
+ " if (callback) {\n",
+ " try {\n",
+ " // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n",
+ " callback(fig, msg);\n",
+ " } catch (e) {\n",
+ " console.log(\"Exception inside the 'handler_\" + msg_type + \"' callback:\", e, e.stack, msg);\n",
+ " }\n",
+ " }\n",
+ " };\n",
+ "}\n",
+ "\n",
+ "// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n",
+ "mpl.findpos = function(e) {\n",
+ " //this section is from http://www.quirksmode.org/js/events_properties.html\n",
+ " var targ;\n",
+ " if (!e)\n",
+ " e = window.event;\n",
+ " if (e.target)\n",
+ " targ = e.target;\n",
+ " else if (e.srcElement)\n",
+ " targ = e.srcElement;\n",
+ " if (targ.nodeType == 3) // defeat Safari bug\n",
+ " targ = targ.parentNode;\n",
+ "\n",
+ " // jQuery normalizes the pageX and pageY\n",
+ " // pageX,Y are the mouse positions relative to the document\n",
+ " // offset() returns the position of the element relative to the document\n",
+ " var x = e.pageX - $(targ).offset().left;\n",
+ " var y = e.pageY - $(targ).offset().top;\n",
+ "\n",
+ " return {\"x\": x, \"y\": y};\n",
+ "};\n",
+ "\n",
+ "/*\n",
+ " * return a copy of an object with only non-object keys\n",
+ " * we need this to avoid circular references\n",
+ " * http://stackoverflow.com/a/24161582/3208463\n",
+ " */\n",
+ "function simpleKeys (original) {\n",
+ " return Object.keys(original).reduce(function (obj, key) {\n",
+ " if (typeof original[key] !== 'object')\n",
+ " obj[key] = original[key]\n",
+ " return obj;\n",
+ " }, {});\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.mouse_event = function(event, name) {\n",
+ " var canvas_pos = mpl.findpos(event)\n",
+ "\n",
+ " if (name === 'button_press')\n",
+ " {\n",
+ " this.canvas.focus();\n",
+ " this.canvas_div.focus();\n",
+ " }\n",
+ "\n",
+ " var x = canvas_pos.x * mpl.ratio;\n",
+ " var y = canvas_pos.y * mpl.ratio;\n",
+ "\n",
+ " this.send_message(name, {x: x, y: y, button: event.button,\n",
+ " step: event.step,\n",
+ " guiEvent: simpleKeys(event)});\n",
+ "\n",
+ " /* This prevents the web browser from automatically changing to\n",
+ " * the text insertion cursor when the button is pressed. We want\n",
+ " * to control all of the cursor setting manually through the\n",
+ " * 'cursor' event from matplotlib */\n",
+ " event.preventDefault();\n",
+ " return false;\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
+ " // Handle any extra behaviour associated with a key event\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.key_event = function(event, name) {\n",
+ "\n",
+ " // Prevent repeat events\n",
+ " if (name == 'key_press')\n",
+ " {\n",
+ " if (event.which === this._key)\n",
+ " return;\n",
+ " else\n",
+ " this._key = event.which;\n",
+ " }\n",
+ " if (name == 'key_release')\n",
+ " this._key = null;\n",
+ "\n",
+ " var value = '';\n",
+ " if (event.ctrlKey && event.which != 17)\n",
+ " value += \"ctrl+\";\n",
+ " if (event.altKey && event.which != 18)\n",
+ " value += \"alt+\";\n",
+ " if (event.shiftKey && event.which != 16)\n",
+ " value += \"shift+\";\n",
+ "\n",
+ " value += 'k';\n",
+ " value += event.which.toString();\n",
+ "\n",
+ " this._key_event_extra(event, name);\n",
+ "\n",
+ " this.send_message(name, {key: value,\n",
+ " guiEvent: simpleKeys(event)});\n",
+ " return false;\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.toolbar_button_onclick = function(name) {\n",
+ " if (name == 'download') {\n",
+ " this.handle_save(this, null);\n",
+ " } else {\n",
+ " this.send_message(\"toolbar_button\", {name: name});\n",
+ " }\n",
+ "};\n",
+ "\n",
+ "mpl.figure.prototype.toolbar_button_onmouseover = function(tooltip) {\n",
+ " this.message.textContent = tooltip;\n",
+ "};\n",
+ "mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Pan axes with left mouse, zoom with right\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n",
+ "\n",
+ "mpl.extensions = [\"eps\", \"jpeg\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\"];\n",
+ "\n",
+ "mpl.default_extension = \"png\";var comm_websocket_adapter = function(comm) {\n",
+ " // Create a \"websocket\"-like object which calls the given IPython comm\n",
+ " // object with the appropriate methods. Currently this is a non binary\n",
+ " // socket, so there is still some room for performance tuning.\n",
+ " var ws = {};\n",
+ "\n",
+ " ws.close = function() {\n",
+ " comm.close()\n",
+ " };\n",
+ " ws.send = function(m) {\n",
+ " //console.log('sending', m);\n",
+ " comm.send(m);\n",
+ " };\n",
+ " // Register the callback with on_msg.\n",
+ " comm.on_msg(function(msg) {\n",
+ " //console.log('receiving', msg['content']['data'], msg);\n",
+ " // Pass the mpl event to the overriden (by mpl) onmessage function.\n",
+ " ws.onmessage(msg['content']['data'])\n",
+ " });\n",
+ " return ws;\n",
+ "}\n",
+ "\n",
+ "mpl.mpl_figure_comm = function(comm, msg) {\n",
+ " // This is the function which gets called when the mpl process\n",
+ " // starts-up an IPython Comm through the \"matplotlib\" channel.\n",
+ "\n",
+ " var id = msg.content.data.id;\n",
+ " // Get hold of the div created by the display call when the Comm\n",
+ " // socket was opened in Python.\n",
+ " var element = $(\"#\" + id);\n",
+ " var ws_proxy = comm_websocket_adapter(comm)\n",
+ "\n",
+ " function ondownload(figure, format) {\n",
+ " window.open(figure.imageObj.src);\n",
+ " }\n",
+ "\n",
+ " var fig = new mpl.figure(id, ws_proxy,\n",
+ " ondownload,\n",
+ " element.get(0));\n",
+ "\n",
+ " // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n",
+ " // web socket which is closed, not our websocket->open comm proxy.\n",
+ " ws_proxy.onopen();\n",
+ "\n",
+ " fig.parent_element = element.get(0);\n",
+ " fig.cell_info = mpl.find_output_cell(\"\");\n",
+ " if (!fig.cell_info) {\n",
+ " console.error(\"Failed to find cell for figure\", id, fig);\n",
+ " return;\n",
+ " }\n",
+ "\n",
+ " var output_index = fig.cell_info[2]\n",
+ " var cell = fig.cell_info[0];\n",
+ "\n",
+ "};\n",
+ "\n",
+ "mpl.figure.prototype.handle_close = function(fig, msg) {\n",
+ " var width = fig.canvas.width/mpl.ratio\n",
+ " fig.root.unbind('remove')\n",
+ "\n",
+ " // Update the output cell to use the data from the current canvas.\n",
+ " fig.push_to_output();\n",
+ " var dataURL = fig.canvas.toDataURL();\n",
+ " // Re-enable the keyboard manager in IPython - without this line, in FF,\n",
+ " // the notebook keyboard shortcuts fail.\n",
+ " IPython.keyboard_manager.enable()\n",
+ " $(fig.parent_element).html('');\n",
+ " fig.close_ws(fig, msg);\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.close_ws = function(fig, msg){\n",
+ " fig.send_message('closing', msg);\n",
+ " // fig.ws.close()\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.push_to_output = function(remove_interactive) {\n",
+ " // Turn the data on the canvas into data in the output cell.\n",
+ " var width = this.canvas.width/mpl.ratio\n",
+ " var dataURL = this.canvas.toDataURL();\n",
+ " this.cell_info[1]['text/html'] = '';\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.updated_canvas_event = function() {\n",
+ " // Tell IPython that the notebook contents must change.\n",
+ " IPython.notebook.set_dirty(true);\n",
+ " this.send_message(\"ack\", {});\n",
+ " var fig = this;\n",
+ " // Wait a second, then push the new image to the DOM so\n",
+ " // that it is saved nicely (might be nice to debounce this).\n",
+ " setTimeout(function () { fig.push_to_output() }, 1000);\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype._init_toolbar = function() {\n",
+ " var fig = this;\n",
+ "\n",
+ " var nav_element = $('')\n",
+ " nav_element.attr('style', 'width: 100%');\n",
+ " this.root.append(nav_element);\n",
+ "\n",
+ " // Define a callback function for later on.\n",
+ " function toolbar_event(event) {\n",
+ " return fig.toolbar_button_onclick(event['data']);\n",
+ " }\n",
+ " function toolbar_mouse_event(event) {\n",
+ " return fig.toolbar_button_onmouseover(event['data']);\n",
+ " }\n",
+ "\n",
+ " for(var toolbar_ind in mpl.toolbar_items){\n",
+ " var name = mpl.toolbar_items[toolbar_ind][0];\n",
+ " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
+ " var image = mpl.toolbar_items[toolbar_ind][2];\n",
+ " var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
+ "\n",
+ " if (!name) { continue; };\n",
+ "\n",
+ " var button = $('');\n",
+ " button.click(method_name, toolbar_event);\n",
+ " button.mouseover(tooltip, toolbar_mouse_event);\n",
+ " nav_element.append(button);\n",
+ " }\n",
+ "\n",
+ " // Add the status bar.\n",
+ " var status_bar = $('');\n",
+ " nav_element.append(status_bar);\n",
+ " this.message = status_bar[0];\n",
+ "\n",
+ " // Add the close button to the window.\n",
+ " var buttongrp = $('');\n",
+ " var button = $('');\n",
+ " button.click(function (evt) { fig.handle_close(fig, {}); } );\n",
+ " button.mouseover('Stop Interaction', toolbar_mouse_event);\n",
+ " buttongrp.append(button);\n",
+ " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n",
+ " titlebar.prepend(buttongrp);\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype._root_extra_style = function(el){\n",
+ " var fig = this\n",
+ " el.on(\"remove\", function(){\n",
+ "\tfig.close_ws(fig, {});\n",
+ " });\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype._canvas_extra_style = function(el){\n",
+ " // this is important to make the div 'focusable\n",
+ " el.attr('tabindex', 0)\n",
+ " // reach out to IPython and tell the keyboard manager to turn it's self\n",
+ " // off when our div gets focus\n",
+ "\n",
+ " // location in version 3\n",
+ " if (IPython.notebook.keyboard_manager) {\n",
+ " IPython.notebook.keyboard_manager.register_events(el);\n",
+ " }\n",
+ " else {\n",
+ " // location in version 2\n",
+ " IPython.keyboard_manager.register_events(el);\n",
+ " }\n",
+ "\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
+ " var manager = IPython.notebook.keyboard_manager;\n",
+ " if (!manager)\n",
+ " manager = IPython.keyboard_manager;\n",
+ "\n",
+ " // Check for shift+enter\n",
+ " if (event.shiftKey && event.which == 13) {\n",
+ " this.canvas_div.blur();\n",
+ " event.shiftKey = false;\n",
+ " // Send a \"J\" for go to next cell\n",
+ " event.which = 74;\n",
+ " event.keyCode = 74;\n",
+ " manager.command_mode();\n",
+ " manager.handle_keydown(event);\n",
+ " }\n",
+ "}\n",
+ "\n",
+ "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
+ " fig.ondownload(fig, null);\n",
+ "}\n",
+ "\n",
+ "\n",
+ "mpl.find_output_cell = function(html_output) {\n",
+ " // Return the cell and output element which can be found *uniquely* in the notebook.\n",
+ " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n",
+ " // IPython event is triggered only after the cells have been serialised, which for\n",
+ " // our purposes (turning an active figure into a static one), is too late.\n",
+ " var cells = IPython.notebook.get_cells();\n",
+ " var ncells = cells.length;\n",
+ " for (var i=0; i= 3 moved mimebundle to data attribute of output\n",
+ " data = data.data;\n",
+ " }\n",
+ " if (data['text/html'] == html_output) {\n",
+ " return [cell, data, j];\n",
+ " }\n",
+ " }\n",
+ " }\n",
+ " }\n",
+ "}\n",
+ "\n",
+ "// Register the function which deals with the matplotlib target/channel.\n",
+ "// The kernel may be null if the page has been refreshed.\n",
+ "if (IPython.notebook.kernel != null) {\n",
+ " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n",
+ "}\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "\n",
+ "list_ui = []\n",
+ "for _argument in list_arguments:\n",
+ " \n",
+ " _arg_ui = widgets.HBox([widgets.Label(_argument,\n",
+ " \n",
+ " layout=widgets.Layout(width='10%')),\n",
+ " widgets.FloatText(value=0,\n",
+ " layout=widgets.Layout(width='20%'))])\n",
+ " _arg_ui.children[1].observe(value_changed, 'value')\n",
+ " \n",
+ " list_ui.append(_arg_ui)\n",
+ " vertical_layout = widgets.VBox(list_ui)\n",
+ " \n",
+ "display(vertical_layout) \n",
+ "plt.clf()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 2
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython2",
+ "version": "2.7.12"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/notebooks/water_fitting.ipynb b/notebooks/water_fitting.ipynb
new file mode 100644
index 0000000..5a1c5b2
--- /dev/null
+++ b/notebooks/water_fitting.ipynb
@@ -0,0 +1,668 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Using LMFIT for Sequential and Global Fit of a Water Sample\n",
+ "\n",
+ "QENS data obtained from a sample of water at room temperature.\n",
+ "\n",
+ "Steps shown in this tutorial:\n",
+ "\n",
+ "- Create a simple model: one elastic line plus one lorentzian plus a linear background.\n",
+ "- Interactively find an initial guess of the model parameters for the spectrum with lowest Q.\n",
+ "- Automatic fit the spectrum with lowest Q and visualize results.\n",
+ "- Automatic sequential fit of the remaining spectra and visualize results.\n",
+ "- Fit the Q-dependence of the Lorentzian's HWHM to a Teixeira model.\n",
+ "- Simultaneous fit of all spectra with two global parameters: water diffusion coefficient and jump time.\n",
+ "- Visualize results from the simultaneous fit.\n",
+ "\n",
+ "### Useful links\n",
+ "- [lmfit documentation](https://lmfit.github.io/lmfit-py/index.html) Curve fitting\n",
+ "- [qef documentation](http://qef.readthedocs.io/en/latest/) Utilities for QENS fitting\n",
+ "- [matplotlib](https://matplotlib.org) Plotting with python\n",
+ "- [Post your questions](https://gitter.im/basisdoc/Lobby)\n",
+ "\n",
+ "
Table of Contents
\n",
+ "- Load data and visualize data \n",
+ "- Define the Fitting Range \n",
+ "- Define the model \n",
+ "- Obtain an initial guess \n",
+ "- Carry out the fit and look at results \n",
+ "- Sequential Fit \n",
+ "- Visualize sequential fits \n",
+ "- Q-dependence of some parameters \n",
+ "- Initial Guess for Teixeira Water model \n",
+ "- Model for Simultaneous Fit of All Spectra with Teixeira Water Model \n",
+ "- Visualize the Simultaneous Fit "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Imports for fitting"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from __future__ import (absolute_import, division, print_function)\n",
+ "\n",
+ "import os\n",
+ "from os.path import join as pjn\n",
+ "import sys\n",
+ "import functools\n",
+ "import lmfit\n",
+ "from lmfit.models import LinearModel, LorentzianModel, ConstantModel, LinearModel\n",
+ "\n",
+ "import qef\n",
+ "from qef.io.loaders import load_nexus\n",
+ "from qef.models.deltadirac import DeltaDiracModel\n",
+ "from qef.models.tabulatedmodel import TabulatedModel\n",
+ "from qef.models.resolution import TabulatedResolutionModel\n",
+ "from qef.models.teixeira import TeixeiraWaterModel\n",
+ "from qef.operators.convolve import Convolve\n",
+ "\n",
+ "# Finding out test data\n",
+ "qef_module_file = os.path.abspath(qef.__file__)\n",
+ "qef_dir = os.path.join(os.path.dirname(qef_module_file), os.pardir)\n",
+ "data_dir = os.path.join(qef_dir, 'tests', 'data', 'io')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from IPython.core.display import HTML\n",
+ "from IPython.core.display import display\n",
+ "from ipywidgets import widgets\n",
+ "\n",
+ "import numpy as np\n",
+ "\n",
+ "import matplotlib.pyplot as plt\n",
+ "%matplotlib notebook"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Imports for plotting and widgets"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "(Top)
\n",
+ "\n",
+ "Starting from the first spectrum, we iteratively fit spectra of higher q's.\n",
+ "\n",
+ "*Special instructions:* If the previous guess was not for the first spectrum, the code will not run but raise and error. Go back to the previous cell and carry out a guess for the first spectrum, then come back here.\n",
+ "\n",
+ "We do not assume any particular Q-dependence for the width of the Lorentzian function."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_spectra = len(fr['y'])\n",
+ "fits = [None,] * n_spectra # store fits for all the tried spectra\n",
+ "fits[0] = fit # store previous fit\n",
+ "for i in range(1, n_spectra):\n",
+ " y_exp = fr['y'][i]\n",
+ " e_exp = fr['e'][i]\n",
+ " fit = model.fit(y_exp, x=e_vals, params=params, weights = 1.0 / e_exp)\n",
+ " params = fit.params # update params with results from previous spectrum's fit\n",
+ " fits[i] = fit # store fit results\n",
+ "\n",
+ "# Show Chi-square versus Q\n",
+ "chi2s = [fit.redchi for fit in fits]\n",
+ "f, ax = plt.subplots()\n",
+ "ax.plot(qs, [fit.redchi for fit in fits])\n",
+ "ax.set_xlabel('Q')\n",
+ "ax.set_ylabel('Chi-squared')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "(Top)
\n",
+ "\n",
+ "The sample is liquid water, thus we expect $EISF \\ll 1$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "names = ('l_fwhm', 'e_amplitude') # fitting parameters we want to plot\n",
+ "ylabels = ('FWHM', 'EISF') # labels on the Y-axis of the plots\n",
+ "xlabels = ('Q^2', 'Q') # labels on the X-axis of the plots\n",
+ "\n",
+ "q_vals = np.asarray(qs)\n",
+ "xs = (q_vals * q_vals, q_vals) # we want to plot FWHM versus Q^2 and EISF versus Q\n",
+ "\n",
+ "f, axs = plt.subplots(1, len(names)) # as many plots as fitting parameters of interest\n",
+ "for i in range(len(names)):\n",
+ " name = names[i] # name of the fit parameter\n",
+ " y = [fit.params[name].value for fit in fits]\n",
+ " ax = axs[i] # select appropriate plotting area\n",
+ " ax.plot(xs[i], y, marker='o', linestyle='dashed')\n",
+ " ax.set_xlabel(xlabels[i])\n",
+ " ax.set_ylabel(ylabels[i])\n",
+ "\n",
+ "plt.tight_layout()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "(Top)
Initial Guess for Teixeira Water model
\n",
+ "\n",
+ "We use the previous $FWHM$ to fit $HWHM(Q^2)$ to Teixeira's water model to obtain initial diffusion $D$ and relaxation time coefficients $\\tau$\n",
+ "\n",
+ "
\n",
+ "\n",
+ "If $Q$ in Angstroms, $HHWM$ in $meV$, and $\\hbar$ in $meV \\cdot ps$, then units of $D$ are $A^2/ps$ and units of $\\tau$ are $ps$.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Collect FHWM from the fits, and associated error in estimation of these optimal values.\n",
+ "hwhms = 0.5 * np.asarray([fit.params['l_fwhm'].value for fit in fits]) # HWHM values\n",
+ "\n",
+ "\n",
+ "# Create the model\n",
+ "from qef.constants import hbar # units of meV x ps or ueV x ns\n",
+ "from lmfit.model import Model\n",
+ "\n",
+ "def teixeira(q2s, difcoef, tau):\n",
+ " r\"\"\"Calculate HWHM for a given Q, diffusion coefficient, and relaxation time\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " q2s : float\n",
+ " Q^2 values\n",
+ " difcoef : float\n",
+ " Diffusion coefficient parameter\n",
+ " tau : float\n",
+ " Relaxation time parameter\n",
+ "\n",
+ " Returns\n",
+ " -------\n",
+ " numpy.ndarray\n",
+ " HWHM values\n",
+ " \"\"\"\n",
+ " dq2 = difcoef * q2s\n",
+ " return hbar * dq2 / (1 + dq2 * tau)\n",
+ "\n",
+ "teixeira_model = Model(teixeira) # create LMFIT Model instance\n",
+ "teixeira_model.set_param_hint('difcoef', min=0) # diffusion coefficient must be positive\n",
+ "teixeira_model.set_param_hint('tau', min=0) # relaxation coefficient must be positive\n",
+ "\n",
+ "\n",
+ "# Carry out the fit\n",
+ "\n",
+ "teixeira_params = teixeira_model.make_params(difcoef=1.0, tau=1.0) # initial guess\n",
+ "teixeira_fit = teixeira_model.fit(hwhms, q2s=np.square(q_vals), params=teixeira_params)\n",
+ "\n",
+ "# Visualize fit results\n",
+ "o_p = teixeira_fit.params # optimal parameters\n",
+ "fmt = 'Chi-square = {}\\nD = {} A^2/ps\\ntau = {} ps'\n",
+ "print(fmt.format(teixeira_fit.redchi, o_p['difcoef'].value, o_p['tau'].value))\n",
+ "teixeira_fit.plot(xlabel='Q^2', ylabel='HWHM')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "(Top)
Model for Simultaneous Fit of All Spectra with Teixeira Water Model
\n",
+ "\n",
+ "We impose a Q-dependende for the FWHM of the Lorentzian, given by the Teixeira water model. Parameters $D$ and $\\tau$ are the only parameters that are same for all spectra.\n",
+ "\n",
+ "
\n",
+ "\n",
+ "We use $D$ and $\\tau$ of the previous fit as initial guesses. We use the sequential fits we did before to initialize all other parameters."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Create a model for each spectrum\n",
+ "\n",
+ "#initialize models and parameter sets for each spectrum\n",
+ "\n",
+ "\n",
+ "# create one model for each spectrum, but collect all parameters under\n",
+ "# a single instance of the Parameters class.\n",
+ "l_model = list()\n",
+ "g_params = lmfit.Parameters()\n",
+ "for i in range(n_spectra):\n",
+ " m, ps = generate_model_and_params(spectrum_index=i) # model and parameters for one of the spectra\n",
+ " l_model.append(m)\n",
+ " [g_params.add(p) for p in ps.values()]\n",
+ "\n",
+ "# Initialize parameter set with the optimized parameters from the sequential fit\n",
+ "for i in range(n_spectra):\n",
+ " optimized_params = fits[i].params # these are I_c, e_amplitude,...\n",
+ " for name in optimized_params:\n",
+ " prefix, base = name.split('_') # for instance, 'e_amplitude' splitted into 'e', and 'amplitude'\n",
+ " i_name = prefix + '_{}_'.format(i) + base # i_name is 'e_3_amplitude' for i=3\n",
+ " g_params[i_name].set(value=optimized_params[name].value)\n",
+ "\n",
+ "# Introduce global parameters diff and tau. Use previous optimized values as initial guess\n",
+ "g_params.add('difcoef', value=o_p['difcoef'].value, min=0)\n",
+ "g_params.add('tau', value=o_p['tau'].value, min=0)\n",
+ "\n",
+ "# Tie each lorentzian l_i_sigma to the teixeira expression\n",
+ "for i in range(n_spectra):\n",
+ " q2 = q_vals[i] * q_vals[i]\n",
+ " teixeira_expression = '{hbar}*difcoef*{q2}/(1+difcoef*{q2}*tau)'\n",
+ " g_params['l_{}_sigma'.format(i)].set(expr=teixeira_expression.format(hbar=hbar, q2=q2))\n",
+ "\n",
+ "print('Number of varying parameters =',len([p for p in g_params.values() if p.vary]),'!')\n",
+ "g_params.pretty_print()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "(Top)
Carry out the Simultaneous Fit
"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def residuals(params):\n",
+ " r\"\"\"Difference between model and experiment, weighted by the experimental error\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " params : lmfit.Parameters\n",
+ " Parameters for the global model\n",
+ "\n",
+ " Returns\n",
+ " -------\n",
+ " numpy.ndarray\n",
+ " 1D array of residuals for the global model\n",
+ " \"\"\"\n",
+ " l_residuals = list()\n",
+ " for i in range(n_spectra):\n",
+ " x = fr['x'] # fitting range of energies\n",
+ " y = fr['y'][i] # associated experimental intensities\n",
+ " e = fr['e'][i] # associated experimental errors\n",
+ " model_evaluation = l_model[i].eval(x=x, params=params)\n",
+ " l_residuals.append((model_evaluation - y) / e)\n",
+ " return np.concatenate(l_residuals)\n",
+ "\n",
+ "# Minimizer object using the parameter set for all models and the\n",
+ "# function to calculate all the residuals.\n",
+ "minimizer = lmfit.Minimizer(residuals, g_params)\n",
+ "g_fit = minimizer.minimize()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print('Chi-square = {:.2f}\\n'.format(g_fit.redchi))\n",
+ "fmt = 'D = {:.3f} A^2/ps, tau = {:.3f} ps'\n",
+ "#print('Before: ', fmt.format(g_fit.init_values['difcoef'], g_fit.init_values['tau']))\n",
+ "print('Before: D = 0.156 A^2/ps, tau = 1.112 ps')\n",
+ "print('After: ', fmt.format(g_fit.params['difcoef'].value, g_fit.params['tau'].value))\n",
+ "print('Teixeira: D = 0.19 A^2/ps, tau = 1.25 ps (J. Teixeira et al., Phys. Rev. A, 31(3), 1913-947 (1985))')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "(Top)
Visualize the Simultaneous Fit
"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "g_fig = plt.figure()\n",
+ "g_gs = plt.GridSpec(nrows=2, ncols=1, height_ratios=[4, 1])\n",
+ "ax_fit = g_fig.add_subplot(g_gs[0]) # host the fit\n",
+ "ax_res = g_fig.add_subplot(g_gs[1], sharex=ax_fit) # host the residuals\n",
+ "\n",
+ "def g_fit_changed(bunch):\n",
+ " i = bunch['new']\n",
+ " ax_fit.clear()\n",
+ " ax_fit.errorbar(fr['x'], fr['y'][i], yerr=fr['e'][i], color='black',\n",
+ " marker='o', markersize=1, markerfacecolor='none', label='data',\n",
+ " linestyle='none')\n",
+ " model_evaluation = l_model[i].eval(x=fr['x'], params=g_params)\n",
+ " ax_fit.plot(fr['x'], model_evaluation, color='red', linewidth=4, label='best fit')\n",
+ " ax_fit.legend()\n",
+ " ax_res.clear()\n",
+ " ax_res.set_xlabel('Energy (meV)')\n",
+ " ax_res.plot(fr['x'], model_evaluation - fr['y'][i], color='black', label='residuals')\n",
+ " ax_res.legend()\n",
+ "\n",
+ "# Widget for the spectrum index\n",
+ "gv_label = widgets.Label('spectrum index', layout=widgets.Layout(width='10%'))\n",
+ "gv_int_text = widgets.BoundedIntText(value=0, min=0, max=(n_spectra - 1),\n",
+ " layout=widgets.Layout(width='20%'))\n",
+ "gv_int_text.observe(g_fit_changed, 'value')\n",
+ "gv_hbox_l = [widgets.HBox([gv_label, gv_int_text])]\n",
+ "\n",
+ "gv_vertical_layout = widgets.VBox(gv_hbox_l)\n",
+ "display(gv_vertical_layout)"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 2",
+ "language": "python",
+ "name": "python2"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 2
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython2",
+ "version": "2.7.12"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/notebooks/widget_builder_on_the_fly.ipynb b/notebooks/widget_builder_on_the_fly.ipynb
new file mode 100644
index 0000000..d28bbb4
--- /dev/null
+++ b/notebooks/widget_builder_on_the_fly.ipynb
@@ -0,0 +1,189 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from IPython.core.display import HTML\n",
+ "from IPython.core.display import display\n",
+ "from ipywidgets import widgets\n",
+ "\n",
+ "import numpy as np\n",
+ "\n",
+ "import matplotlib.pyplot as plt\n",
+ "%matplotlib notebook"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "726a1be1595849c7bc5bfc91ea5449a6",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "A Jupyter Widget"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "user_input_ui = widgets.HBox([widgets.Label(\"Define Function:\",\n",
+ " layout=widgets.Layout(width='10%')),\n",
+ " widgets.Text(value='',\n",
+ " placeholder=\"Define function here\",\n",
+ " layout=widgets.Layout(width='60%')),\n",
+ " widgets.Label(\"Ex: f(a,b,c)=a+b*c\",\n",
+ " layout=widgets.Layout(width='20%'))])\n",
+ "display(user_input_ui)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "..."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "At this point, you retrieve the name of the coefficient defined in the function, let's pretend those are \n",
+ "\n",
+ " * argument_1\n",
+ " * argument_2\n",
+ " * argument_3\n",
+ " \n",
+ " ```\n",
+ " function(x, argument_1, argument_2, argument_3) = 3*argument_1*x + argument_2 * argument_3\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "list_arguments = ['argument_1', 'argument_2', 'argument_3']"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def function(x, argument_1, argument_2, argument_3):\n",
+ " return 3*argument_1*x*x + argument_2*x + argument_3"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def value_changed(value):\n",
+ " # retrieve value of all widgets\n",
+ " list_value = []\n",
+ " for _top_widget in vertical_layout.children:\n",
+ " _text_float = _top_widget.children[1]\n",
+ " list_value.append(_text_float.value)\n",
+ "\n",
+ " x = np.linspace(0,100)\n",
+ " plt.clf()\n",
+ " plt.plot(x, function(x, list_value[0], list_value[1], list_value[2]))\n",
+ " plt.title(\"argument_1: {}, argument_2: {}, argument_3: {}\".format(list_value[0], list_value[1], list_value[2]))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "2d8bd4928ce143c1b41d1ace62a16e7c",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "A Jupyter Widget"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "\n",
+ "list_ui = []\n",
+ "for _argument in list_arguments:\n",
+ " \n",
+ " _arg_ui = widgets.HBox([widgets.Label(_argument,\n",
+ " \n",
+ " layout=widgets.Layout(width='10%')),\n",
+ " widgets.FloatText(value=0,\n",
+ " layout=widgets.Layout(width='20%'))])\n",
+ " _arg_ui.children[1].observe(value_changed, 'value')\n",
+ " \n",
+ " list_ui.append(_arg_ui)\n",
+ " vertical_layout = widgets.VBox(list_ui)\n",
+ " \n",
+ "display(vertical_layout) \n",
+ "plt.clf()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 2
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython2",
+ "version": "2.7.12"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/qef/constants.py b/qef/constants.py
index 8ee34c9..8c950de 100644
--- a/qef/constants.py
+++ b/qef/constants.py
@@ -4,4 +4,4 @@
from scipy import constants
planck_constant = constants.Planck / constants.e * 1E15 # meV*psec
-hbar = planck_constant / (2 * np.pi)
+hbar = planck_constant / (2 * np.pi) # meV*psec
diff --git a/qef/io/loaders.py b/qef/io/loaders.py
index 1b71de5..d94f684 100644
--- a/qef/io/loaders.py
+++ b/qef/io/loaders.py
@@ -72,7 +72,7 @@ def load_nexus_processed(file_name):
w = data['workspace']
x = w['axis1'].value # energy or time values
y = w['values'].value # intensities
- e = w['errors'].value # indeterminacies in the intensities
+ e = w['errors'].value # undeterminacies in the intensities
# Transform to point data
if len(x) == 1 + len(y[0]):
x = histogram_to_point_data(x)
diff --git a/qef/models/deltadirac.py b/qef/models/deltadirac.py
index e36d97c..8fdfde9 100644
--- a/qef/models/deltadirac.py
+++ b/qef/models/deltadirac.py
@@ -32,7 +32,9 @@ class DeltaDiracModel(Model):
to the center parameter.
At value-closest-to-center, the model evaluates to the amplitude
- parameter divided by the x-spacing.
+ parameter divided by the x-spacing. This last division is
+ necessary to preserve normalization with integrating the function
+ over the X-axis
Fitting parameters:
- integrated intensity ``amplitude`` :math:`A`
diff --git a/qef/models/tabulatedmodel.py b/qef/models/tabulatedmodel.py
index d2a7ad9..dd9e883 100644
--- a/qef/models/tabulatedmodel.py
+++ b/qef/models/tabulatedmodel.py
@@ -1,5 +1,6 @@
from __future__ import (absolute_import, division, print_function)
+import numpy as np
from scipy.interpolate import interp1d
from lmfit import Model, models
@@ -21,12 +22,19 @@ class TabulatedModel(Model):
"""
def __init__(self, xs, ys, *args, **kwargs):
- self._interp = interp1d(xs, ys, fill_value='extrapolate', kind='cubic')
+ x_1d = xs.reshape(xs.size)
+ y_1d = ys.reshape(ys.size)
+ y_at_xmin = y_1d[np.argmin(x_1d)]
+ y_at_xmax = y_1d[np.argmax(x_1d)]
+ self._interp = interp1d(x_1d, y_1d, fill_value=(y_at_xmin, y_at_xmax),
+ bounds_error=False, kind='linear')
def interpolator(x, amplitude, center):
return amplitude * self._interp(x - center)
super(TabulatedModel, self).__init__(interpolator, *args, **kwargs)
+ self.set_param_hint('amplitude', value=1.0)
+ self.set_param_hint('center', value=0.0)
def guess(self, data, x, **kwargs):
@@ -51,12 +59,12 @@ def guess(self, data, x, **kwargs):
"""
params = self.make_params()
- def pset(param, value, min):
- params["%s%s" % (self.prefix, param)].set(value=value, min=min)
+ def pset(param, value):
+ params["%s%s" % (self.prefix, param)].set(value=value)
x_at_max = x[models.index_of(data, max(data))]
ysim = self.eval(x=x_at_max, amplitude=1, center=x_at_max)
amplitude = max(data) / ysim
- pset("amplitude", amplitude, min=0.0)
- pset("center", x_at_max, min=-1000)
+ pset("amplitude", amplitude)
+ pset("center", x_at_max)
return models.update_param_vals(params, self.prefix, **kwargs)
diff --git a/qef/operators/convolve.py b/qef/operators/convolve.py
index 13c24ca..8ca52c8 100644
--- a/qef/operators/convolve.py
+++ b/qef/operators/convolve.py
@@ -33,6 +33,9 @@ class Convolve(CompositeModel):
It is assumed that the resolution FWHM is energy independent.
Non-symmetric energy ranges are allowed (when the range of negative values
is different than that of positive values).
+
+ The convolution requires multiplication by the X-spacing to preserve
+ normalization
"""
def __init__(self, resolution, model, **kws):
@@ -50,5 +53,6 @@ def eval(self, params=None, **kwargs):
e = np.concatenate((neg_e, e, pos_e))
kwargs.update({independent_var: e})
model_data = self.model.eval(params=params, **kwargs)
+ # Multiply by the X-spacing to preserve normalization
de = (e[-1] - e[0])/(len(e) - 1) # energy spacing
return de * convolve(model_data, res_data)
diff --git a/requirements.txt b/requirements.txt
index 0a0589f..14dd2a1 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,6 @@
+h5py
+lmfit
+matplotlib
numpy
scipy
-lmfit
six
-h5py
diff --git a/setup.py b/setup.py
index 5b85c52..cebdd6f 100644
--- a/setup.py
+++ b/setup.py
@@ -12,7 +12,7 @@
history = history_file.read()
requirements = [
- 'numpy', 'scipy', 'lmfit', 'six', 'h5py'
+ 'h5py', 'lmfit', 'numpy', 'matplotlib', 'scipy', 'six'
# TODO: put package requirements here
]
diff --git a/tests/conftest.py b/tests/conftest.py
index a5ae56d..6a55393 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -24,4 +24,8 @@ def ltz():
@pytest.fixture(scope='session')
def io_fix():
return dict(irs_res_f=pjn(data_dir, 'io', 'irs26173_graphite_res.nxs'),
- irs_red_f=pjn(data_dir, 'io', 'irs26176_graphite002_red.nxs'))
+ irs_red_f=pjn(data_dir, 'io', 'irs26176_graphite002_red.nxs'),
+ q_values=(0.525312757876, 0.7291668809127, 0.9233951329944,
+ 1.105593679447, 1.273206832528, 1.42416584459,
+ 1.556455009584, 1.668282739099, 1.758225254224,
+ 1.825094271503))
diff --git a/tests/integration/__init__.py b/tests/integration/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/integration/test_water.py b/tests/integration/test_water.py
new file mode 100644
index 0000000..090ea0c
--- /dev/null
+++ b/tests/integration/test_water.py
@@ -0,0 +1,172 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+from __future__ import (absolute_import, division, print_function)
+
+import pytest
+import os
+import numpy as np
+from numpy.testing import assert_almost_equal
+
+from lmfit.models import LinearModel, LorentzianModel, ConstantModel
+import lmfit
+from lmfit.model import Model
+
+from qef.constants import hbar # units of meV x ps or ueV x ns
+from qef.io.loaders import load_nexus
+from qef.models.deltadirac import DeltaDiracModel
+from qef.models.resolution import TabulatedResolutionModel
+from qef.operators.convolve import Convolve
+
+
+def test_water(io_fix):
+ # Load data
+ res = load_nexus(io_fix['irs_res_f'])
+ dat = load_nexus(io_fix['irs_red_f'])
+ q_vals = io_fix['q_values']
+
+ # Define the fitting range
+ e_min = -0.4
+ e_max = 0.4
+ # Find indexes of dat['x'] with values in (e_min, e_max)
+ mask = np.intersect1d(np.where(dat['x'] > e_min),
+ np.where(dat['x'] < e_max))
+ # Drop data outside the fitting range
+ fr = dict() # fitting range. Use in place of 'dat'
+ fr['x'] = dat['x'][mask]
+ fr['y'] = np.asarray([y[mask] for y in dat['y']])
+ fr['e'] = np.asarray([e[mask] for e in dat['e']])
+
+ # Create the model
+ def generate_model_and_params(spectrum_index=None):
+ r"""Produce an LMFIT model and related set of fitting parameters"""
+
+ sp = '' if spectrum_index is None else '{}_'.format(
+ spectrum_index) # prefix if spectrum_index passed
+
+ # Model components
+ intensity = ConstantModel(prefix='I_' + sp) # I_amplitude
+ elastic = DeltaDiracModel(prefix='e_' + sp) # e_amplitude, e_center
+ # l_amplitude, l_center, l_sigma (also l_fwhm, l_height)
+ inelastic = LorentzianModel(prefix='l_' + sp)
+ # r_amplitude, r_center (both fixed)
+ resolution = TabulatedResolutionModel(res['x'], res['y'],
+ prefix='r_' + sp)
+ background = LinearModel(prefix='b_' + sp) # b_slope, b_intercept
+
+ # Putting it all together
+ model = intensity * Convolve(resolution,
+ elastic + inelastic) + background
+ parameters = model.make_params() # model params are a separate entity
+
+ # Ties and constraints
+ parameters['e_' + sp + 'amplitude'].set(min=0.0, max=1.0)
+ parameters['l_' + sp + 'center'].set(
+ expr='e_' + sp + 'center') # centers tied
+ parameters['l_' + sp + 'amplitude'].set(
+ expr='1 - e_' + sp + 'amplitude')
+
+ # Some initial sensible values
+ init_vals = {'I_' + sp + 'c': 1.0, 'e_' + sp + 'amplitude': 0.5,
+ 'l_' + sp + 'sigma': 0.01,
+ 'b_' + sp + 'slope': 0, 'b_' + sp + 'intercept': 0}
+ for p, v in init_vals.items():
+ parameters[p].set(value=v)
+
+ return model, parameters
+ # Call the function
+ model, params = generate_model_and_params()
+
+ # Initial guess for first spectrum. Only set free parameters
+ for name, value in dict(I_c=4.0, e_center=0, e_amplitude=0.1,
+ l_sigma=0.03, b_slope=0, b_intercept=0).items():
+ params[name].set(value=value)
+ # Carry out the fit
+ fit = model.fit(fr['y'][0], x=fr['x'], params=params,
+ weights=1.0 / fr['e'][0])
+ assert_almost_equal(fit.redchi, 1.72, decimal=2)
+
+ # Carry out sequential fit
+ n_spectra = len(fr['y'])
+ fits = [None, ] * n_spectra # store fits for all the tried spectra
+ fits[0] = fit # store previous fit
+ for i in range(1, n_spectra):
+ y_exp = fr['y'][i]
+ e_exp = fr['e'][i]
+ fit = model.fit(y_exp, x=fr['x'], params=params, weights=1.0 / e_exp)
+ fits[i] = fit # store fit results
+ assert_almost_equal([f.redchi for f in fits],
+ [1.72, 1.15, 0.81, 0.73, 0.73, 0.75, 0.81, 0.86, 0.75,
+ 0.91],
+ decimal=2)
+
+ # Fit HWHM(Q^2) with Teixeira model
+ hwhms = 0.5 * np.asarray([fit.params['l_fwhm'].value for fit in fits])
+
+ def teixeira(q2s, difcoef, tau):
+ dq2 = difcoef * q2s
+ return hbar * dq2 / (1 + dq2 * tau)
+ teixeira_model = Model(teixeira) # create LMFIT Model instance
+ teixeira_model.set_param_hint('difcoef', min=0)
+ teixeira_model.set_param_hint('tau', min=0)
+ # Carry out the fit from an initial guess
+ teixeira_params = teixeira_model.make_params(difcoef=1.0, tau=1.0)
+ teixeira_fit = teixeira_model.fit(hwhms, q2s=np.square(q_vals),
+ params=teixeira_params)
+ assert_almost_equal([teixeira_fit.best_values['difcoef'],
+ teixeira_fit.best_values['tau']],
+ [0.16, 1.11], decimal=2)
+
+ # Model for Simultaneous Fit of All Spectra with Teixeira Water Model
+ #
+ # create one model for each spectrum, but collect all parameters under
+ # a single instance of the Parameters class.
+ l_model = list()
+ g_params = lmfit.Parameters()
+ for i in range(n_spectra):
+ # model and parameters for one of the spectra
+ m, ps = generate_model_and_params(spectrum_index=i)
+ l_model.append(m)
+ [g_params.add(p) for p in ps.values()]
+
+ # Initialize parameter set with optimized parameters from sequential fit
+ for i in range(n_spectra):
+ optimized_params = fits[i].params # these are I_c, e_amplitude,...
+ for name in optimized_params:
+ # for instance, 'e_amplitude' splitted into 'e', and 'amplitude'
+ prefix, base = name.split('_')
+ # i_name is 'e_3_amplitude' for i=3
+ i_name = prefix + '_{}_'.format(i) + base
+ g_params[i_name].set(value=optimized_params[name].value)
+
+ # Introduce global parameters difcoef and tau.
+ # Use previous optimized values as initial guess
+ o_p = teixeira_fit.params
+ g_params.add('difcoef', value=o_p['difcoef'].value, min=0)
+ g_params.add('tau', value=o_p['tau'].value, min=0)
+
+ # Tie each lorentzian l_i_sigma to the teixeira expression
+ for i in range(n_spectra):
+ q2 = q_vals[i] * q_vals[i]
+ fmt = '{hbar}*difcoef*{q2}/(1+difcoef*{q2}*tau)'
+ teixeira_expression = fmt.format(hbar=hbar, q2=q2)
+ g_params['l_{}_sigma'.format(i)].set(expr=teixeira_expression)
+
+ # Carry out the Simultaneous Fit
+ def residuals(params):
+ l_residuals = list()
+ for i in range(n_spectra):
+ x = fr['x'] # fitting range of energies
+ y = fr['y'][i] # associated experimental intensities
+ e = fr['e'][i] # associated experimental errors
+ model_evaluation = l_model[i].eval(x=x, params=params)
+ l_residuals.append((model_evaluation - y) / e)
+ return np.concatenate(l_residuals)
+ # Minimizer object using the parameter set for all models and the
+ # function to calculate all the residuals.
+ minimizer = lmfit.Minimizer(residuals, g_params)
+ g_fit = minimizer.minimize()
+ assert_almost_equal(g_fit.redchi, 0.93, decimal=2)
+
+
+if __name__ == '__main__':
+ pytest.main([os.path.abspath(__file__)])
diff --git a/tests/models/test_deltadirac.py b/tests/models/test_deltadirac.py
index b20498b..65710d0 100644
--- a/tests/models/test_deltadirac.py
+++ b/tests/models/test_deltadirac.py
@@ -25,25 +25,39 @@ def test_guess():
def test_convolution():
+ r"""Convolution of function with delta dirac should return the function"""
+
+ # Reference Lorentzian parameter values
amplitude = 42.0
sigma = 0.042
center = 0.0003
+
c1 = LorentzianModel(prefix='c1_')
p = c1.make_params(amplitude=amplitude, center=center, sigma=sigma)
c2 = DeltaDiracModel(prefix='c2_')
p.update(c2.make_params(amplitude=1.0, center=0.0))
+
e = 0.0004 * np.arange(-250, 1500) # energies in meV
+ # convolve Lorentzian with delta Dirac
y1 = Convolve(c1, c2).eval(params=p, x=e) # should be the lorentzian
+ # reverse order, convolve delta Dirac with Lorentzian
y2 = Convolve(c2, c1).eval(params=p, x=e) # should be the lorentzian
+
+ # We will fit a Lorentzian model against datasets y1 and y2
m = LorentzianModel()
all_params = 'amplitude sigma center'.split()
for y in (y1, y2):
params = m.guess(y, x=e)
- # Set initial parameters far from optimal solution
+ # Set initial model Lorentzian parameters far from optimal solution
params['amplitude'].set(value=amplitude * 10)
params['sigma'].set(value=sigma * 4)
params['center'].set(value=center * 7)
+
+ # fit Lorentzian model against dataset y
r = m.fit(y, params, x=e)
+
+ # Compare the reference Lorentzian parameters against
+ # parameters of the fitted model
assert_allclose([amplitude, sigma, center],
[r.params[p].value for p in all_params],
rtol=0.01, atol=0.00001)
diff --git a/tests/models/test_tabulatedmodel.py b/tests/models/test_tabulatedmodel.py
index d4e6b7b..5528e6c 100644
--- a/tests/models/test_tabulatedmodel.py
+++ b/tests/models/test_tabulatedmodel.py
@@ -17,7 +17,7 @@ def test_tabulatedmodel():
sigma=0.042)
model = TabulatedModel(x_sim, y_sim)
- params = model.guess(x_exp, y_exp)
+ params = model.guess(y_exp, x_exp)
fit = model.fit(y_exp, params, x=x_exp)
assert abs(fit.best_values['amplitude'] - intensity) < 0.0001
diff --git a/tests/operators/test_convolve.py b/tests/operators/test_convolve.py
index 049d091..a103dbf 100644
--- a/tests/operators/test_convolve.py
+++ b/tests/operators/test_convolve.py
@@ -22,6 +22,9 @@
def test_simplecases(ComponentModel, de, sigma):
r"""Convolution of two Lorentzians is one Lorentzian, and convolution
of two gaussians is one gaussian"""
+
+ # Convolve two Lorentzians or convolve two Gaussians. Store the result
+ # in dataset 'y'
s1 = 0.011 # narrow component
s2 = 0.163 # broad component
c1 = ComponentModel(prefix='c1_')
@@ -30,12 +33,21 @@ def test_simplecases(ComponentModel, de, sigma):
p.update(c2.make_params(amplitude=1.0, center=0.0, sigma=s2))
e = de * np.arange(-250, 1500) # energies in meV
y = Convolve(c1, c2).eval(params=p, x=e)
+
+ # Fit a Lorentzian to dataset 'y' or a Gaussian to dataset 'y'
m = ComponentModel()
params = m.guess(y, x=e)
params['amplitude'].set(value=42.0) # astray from correct value
params['sigma'].set(value=4.2*sigma(s1, s2)) # astray from correct value
r = m.fit(y, params, x=e)
+
+ # The assertions will be carried out using not the dataset 'y' but
+ # the fit we did to the dataset 'y'.
+ # 1. Assert the convolution preserves the normalization
assert_almost_equal(r.params['amplitude'], 1.0, decimal=3)
+
+ # 1. Assert the convolution results in the appropriate combination for
+ # their components' sigma values
assert_almost_equal(r.params['sigma'], sigma(s1, s2), decimal=3)