From 8cb225efefa2d690ddfe37ed9221a1d4c45c383a Mon Sep 17 00:00:00 2001 From: Jose Borreguero Date: Mon, 16 Apr 2018 18:22:20 -0400 Subject: [PATCH 1/5] Refs #49 copied notebook from water fitting --- notebooks/strexp_fitting.ipynb | 714 +++++++++++++++++++++++++++++++++ 1 file changed, 714 insertions(+) create mode 100644 notebooks/strexp_fitting.ipynb diff --git a/notebooks/strexp_fitting.ipynb b/notebooks/strexp_fitting.ipynb new file mode 100644 index 0000000..e6f8b28 --- /dev/null +++ b/notebooks/strexp_fitting.ipynb @@ -0,0 +1,714 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Using LMFIT for Sequential and Global Fit of a Polymer 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 Fourier transform of the stretched exponential, plus a linear background, plus a tabulated model representing the empty can.\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 relaxation time to a power law.\n", + "- Simultaneous fit of all spectra with one global parameter: stretching exponent.\n", + "- Visualize results from the simultaneous fit.\n", + "\n", + "### Useful links\n", + "- [qef documentation](http://qef.readthedocs.io/en/latest/) (pip install qef) Utilities for QENS fitting\n", + "- [lmfit documentation](https://lmfit.github.io/lmfit-py/index.html) Curve fitting\n", + "- [matplotlib](https://matplotlib.org) Plotting with python\n", + "- [Post your questions](https://gitter.im/basisdoc/Lobby)\n", + "\n", + "

Table of Contents

\n", + "- Donwload Data \n", + "- Load Data and Visualize \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" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Imports for plotting and widgets" + ] + }, + { + "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": [ + "(Top)

Donwload Data

\n", + "\n", + "It's assumed git is installed in your system. Otherwise,\n", + "[follow instructions](http://qef.readthedocs.io/en/latest/installation.html#testing-tutorials-data)\n", + "to download and unpack your data to /tmp/qef_data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "qef_data_dir=\"/tmp/qef_data\"\n", + "git clone https://github.com/jmborr/qef_data ${qef_data_dir}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "qef_data_dir = '/tmp/qef_data'\n", + "data_dir = os.path.join(qef_data_dir, 'data', 'notebooks', 'strexp')\n", + "print(data_dir)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "(Top)

Load data and visualize data

" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "res = load_nexus(pjn(data_dir,'irs26173_graphite_res.nxs'))\n", + "emin, emax = np.min(res['x']), np.max(res['x'])\n", + "dat = load_nexus(pjn(data_dir,'irs26176_graphite002_red.nxs')) # data has 10 histograms\n", + "emin, emax = min(emin, np.min(dat['x'])), max(emax, np.max(dat['x']))\n", + "print('resolution range is ({:.4f}, {:.4f})'.format(res['x'][0], res['x'][-1]))\n", + "print('data range is ({:.4f}, {:.4f})'.format(dat['x'][0], dat['x'][-1]))\n", + "qs = (0.525312757876, 0.7291668809127, 0.9233951329944, 1.105593679447, 1.273206832528, 1.42416584459, 1.556455009584, 1.668282739099, 1.758225254224, 1.825094271503)\n", + "# Plot\n", + "f, (ax1, ax2) = plt.subplots(1, 2)\n", + "ax1.semilogy(res['x'], res['y'][0]); ax1.set_title('resolution')\n", + "[ax2.semilogy(dat['x'], data) for data in dat['y']]; ax2.set_title('sample')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "(Top)

Define the Fitting Range

" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "e_min = -0.4 ; e_max = 0.4\n", + "\n", + "# Find indexes of dat['x'] with values in (e_min, e_max)\n", + "mask = np.intersect1d(np.where(dat['x']>e_min), np.where(dat['x']Top)

Define the model

\n", + "\n", + "
$S(E) = I \\cdot R(Q,E) \\otimes \\big( eisf \\cdot \\delta(E-E_0) + (1-eisf) \\cdot L(E-E_0) \\big) + LB$
\n", + "\n", + "We will use the following component models:\n", + "- ConstantModel represents one number that can be fitted\n", + "- DeltaDiracModel\n", + "- LorentzianModel has parameter $sigma \\equiv HWHM$\n", + "- TabulatedResolutionModel to store the table of numbers representing the resolution\n", + "- LinearModel for the linear background" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create the Model. We put everything under a function which we'll later reuse\n", + "\n", + "def generate_model_and_params(spectrum_index=None):\n", + " r\"\"\"Produce an LMFIT model and related set of fitting parameters\"\"\"\n", + "\n", + " sp = '' if spectrum_index is None else '{}_'.format(spectrum_index) # prefix if spectrum_index passed\n", + "\n", + " # Model components\n", + " intensity = ConstantModel(prefix='I_'+sp) # I_amplitude\n", + " elastic = DeltaDiracModel(prefix='e_'+sp) # e_amplitude, e_center\n", + " inelastic = LorentzianModel(prefix='l_'+sp) # l_amplitude, l_center, l_sigma (also l_fwhm, l_height)\n", + " resolution = TabulatedResolutionModel(res['x'], res['y'], prefix='r_'+sp) # (fixed r_amplitude, r_center)\n", + " background = LinearModel(prefix='b_'+sp) # b_slope, b_intercept\n", + "\n", + " # Putting it all together\n", + " model = intensity * Convolve(resolution, elastic + inelastic) + background\n", + " parameters = model.make_params() # model parameters are a separate entity.\n", + "\n", + " # Ties and constraints\n", + " parameters['e_'+sp+'amplitude'].set(min=0.0, max=1.0)\n", + " parameters['l_'+sp+'center'].set(expr='e_'+sp+'center') # centers tied\n", + " parameters['l_'+sp+'amplitude'].set(expr='1 - e_'+sp+'amplitude')\n", + "\n", + " # Some initial sensible values\n", + " init_vals = {'I_'+sp+'c': 1.0, 'e_'+sp+'amplitude': 0.5, 'l_'+sp+'sigma': 0.01,\n", + " 'b_'+sp+'slope': 0, 'b_'+sp+'intercept': 0}\n", + " for p, v in init_vals.items():\n", + " parameters[p].set(value=v)\n", + "\n", + " return model, parameters\n", + "\n", + "# Call the function\n", + "model, params = generate_model_and_params()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "(Top)

Obtain an initial guess

\n", + "\n", + "A widget that compares the evaluation of the model with one of the experimental spectra. You can tweak only the free (unconstrained) parameters.\n", + "\n", + "When run, you will see two empty panels, one for comparison between experiment and model, and the second panel for residuals. Start changing the values of the parameters for the panels to be populated.\n", + "\n", + "A good initial guess for the first spectrum is :\n", + "\n", + "> spectum index = 0 \n", + "> I_c = 4 \n", + "> e_center = 0 \n", + "> e_amplitude = 0.1 \n", + "> l_sigma = 0.03 \n", + "> b_intercept = 0 \n", + "> b_slope = 0 \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Variables from previous cells and setting of other variables\n", + "# params\n", + "# e_vals\n", + "# model\n", + "free_params = [p for p in params.values() if not p.expr and p.vary] # only adjust free parameters\n", + "e_vals = fr['x']\n", + "y_exp = fr['y'][0] # associated experimental values for the first histogram\n", + "e_exp = fr['e'][0]\n", + "\n", + "f = plt.figure()\n", + "gs = plt.GridSpec(nrows=2, ncols=1, height_ratios=[4, 1])\n", + "ax1 = f.add_subplot(gs[0]) # host the fit\n", + "ax2 = f.add_subplot(gs[1], sharex=ax1) # host the residuals\n", + "\n", + "#f, (ax1, ax2) = plt.subplots(2, 1)\n", + "def plot_new_spectrum(an_axis):\n", + " global y_exp\n", + " an_axis.clear()\n", + " an_axis.plot(e_vals, y_exp, color='black', marker='o',\n", + " markersize=1.0, linewidth=0, label='experiment')\n", + " an_axis.legend()\n", + "\n", + "def plot_guess(an_axis, model_evaluation):\n", + " plot_new_spectrum(an_axis)\n", + " an_axis.plot(e_vals, model_evaluation, color='blue', label='model')\n", + " an_axis.legend()\n", + "\n", + "def plot_difference(an_axis, model_evaluation):\n", + " global y_exp\n", + " an_axis.clear()\n", + " an_axis.plot(e_vals, y_exp - model_evaluation, color='black',\n", + " markersize=1.0, label='exp - model')\n", + " an_axis.legend()\n", + "\n", + "def i_histogram_changed(bunch):\n", + " global y_exp\n", + " y_exp = fr['y'][bunch['new']]\n", + " global e_exp\n", + " e_exp = fr['e'][bunch['new']]\n", + " plot_new_spectrum(ax1)\n", + " ax2.clear()\n", + "\n", + "# Widget for the spectrum index\n", + "w_label = widgets.Label('spectrum index', layout=widgets.Layout(width='10%'))\n", + "w_int_text = widgets.BoundedIntText(value=0, min=0, max=len(dat['y']),\n", + " layout=widgets.Layout(width='20%'))\n", + "w_int_text.observe(i_histogram_changed, 'value')\n", + "p_hbox_l = [widgets.HBox([w_label, w_int_text])]\n", + "\n", + "def update_model(name, value, parameters):\n", + " parameters[name].set(value=value)\n", + " return model.eval(x=e_vals, params=parameters)\n", + "\n", + "widget_to_parameter = dict()\n", + "def parameter_changed(bunch):\n", + " w_float_text = bunch['owner']\n", + " p_name = widget_to_parameter[w_float_text]\n", + " value = bunch['new']\n", + " w_float_text.step = 0.1 * value # update the step as 10% of current value\n", + " model_evaluation = update_model(p_name, value, params)\n", + " plot_guess(ax1, model_evaluation)\n", + " plot_difference(ax2, model_evaluation)\n", + " \n", + "def p_hbox(p):\n", + " \"\"\"Generate an HBox widget for a given fitting parameter\n", + "\n", + " Parameters\n", + " ----------\n", + " p : lmfit parameter\n", + " \"\"\"\n", + " w_label = widgets.Label(p.name, layout=widgets.Layout(width='10%'))\n", + " w_float_text = widgets.FloatText(value=p.value, layout=widgets.Layout(width='20%'))\n", + " w_float_text.step = 0.1 * p.value\n", + " w_float_text.observe(parameter_changed, 'value')\n", + " widget_to_parameter[w_float_text] = p.name\n", + " return widgets.HBox([w_label, w_float_text])\n", + "\n", + "p_hbox_l.extend([p_hbox(p) for p in free_params])\n", + "\n", + "vertical_layout = widgets.VBox(p_hbox_l)\n", + "display(vertical_layout)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "(Top)

Carry out the fit and look at results

" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fit = model.fit(y_exp, x=e_vals, params=params, weights = 1.0 / e_exp)\n", + "print('Chi-square =', fit.redchi)\n", + "fit.plot(data_kws=dict(color='black', marker='o', markersize=1, markerfacecolor='none'),\n", + " fit_kws=dict(color='red', linewidth=4))\n", + "print('\\n'.join('{} = {}'.format(p.name, p.value) for p in fit.params.values()))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "(Top)

Sequential Fit

\n", + "\n", + "**Special instructions:** If the fit you carried out in the previous cell was not for the first spectrum, the sequential fit will not run but raise an error. Go back to the cell for the Initial Guess and carry out a guess and a subsequent fit for the first spectrum, then come back here.\n", + "\n", + "Starting from the first spectrum, we iteratively fit spectra of higher q's.\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)

Visualize sequential fits

\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sv_fig = plt.figure()\n", + "def fit_changed(bunch):\n", + " sv_fig.clear()\n", + " fits[bunch['new']].plot(fig=sv_fig,\n", + " data_kws=dict(color='black', marker='o', markersize=1, markerfacecolor='none'),\n", + " fit_kws=dict(color='red', linewidth=4))\n", + "\n", + "# Widget for the spectrum index\n", + "sv_label = widgets.Label('spectrum index', layout=widgets.Layout(width='10%'))\n", + "sv_int_text = widgets.BoundedIntText(value=0, min=0, max=(n_spectra - 1),\n", + " layout=widgets.Layout(width='20%'))\n", + "sv_int_text.observe(fit_changed, 'value')\n", + "sv_hbox_l = [widgets.HBox([sv_label, sv_int_text])]\n", + "\n", + "sv_vertical_layout = widgets.VBox(sv_hbox_l)\n", + "display(sv_vertical_layout)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "(Top)

Q-dependence of some parameters

\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", + "
$HWHM = \\frac{\\hbar D\\cdot Q^2}{1 + D \\cdot Q^2 \\cdot \\tau}$
\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", + "$S(Q, E) = I \\cdot R(Q,E) \\otimes \\big( eisf \\cdot \\delta(Q, E-E_0) + (1-eisf) \\cdot L(Q, E-E_0, FWHM = \\frac{2 \\hbar D\\cdot Q^2}{1 + D \\cdot Q^2 \\cdot \\tau}) \\big) + LB$\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": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From eedab22627c5f0942667e15558f918987af75b50 Mon Sep 17 00:00:00 2001 From: Jose Borreguero Date: Mon, 16 Apr 2018 18:35:41 -0400 Subject: [PATCH 2/5] Refs #1 Download data working --- notebooks/strexp_fitting.ipynb | 53 +++++++++++++++++++++++++++------- 1 file changed, 43 insertions(+), 10 deletions(-) diff --git a/notebooks/strexp_fitting.ipynb b/notebooks/strexp_fitting.ipynb index e6f8b28..4226e4c 100644 --- a/notebooks/strexp_fitting.ipynb +++ b/notebooks/strexp_fitting.ipynb @@ -48,9 +48,18 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/sw/miniconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n", + " from ._conv import register_converters as _register_converters\n" + ] + } + ], "source": [ "from __future__ import (absolute_import, division, print_function)\n", "\n", @@ -79,7 +88,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -106,24 +115,47 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "First, rewinding head to replay your work on top of it...\n", + "Fast-forwarded master to ad9365c5454e952489223a29cacd15dc003afe2a.\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "From https://github.com/jmborr/qef_data\n", + " ee29426..ad9365c master -> origin/master\n" + ] + } + ], "source": [ "%%bash\n", "qef_data_dir=\"/tmp/qef_data\"\n", - "git clone https://github.com/jmborr/qef_data ${qef_data_dir}" + "if [ -d ${qef_data_dir} ]; then\n", + " cd ${qef_data_dir}\n", + " git pull --rebase\n", + "else\n", + " git clone https://github.com/jmborr/qef_data ${qef_data_dir}\n", + "fi" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "qef_data_dir = '/tmp/qef_data'\n", "data_dir = os.path.join(qef_data_dir, 'data', 'notebooks', 'strexp')\n", - "print(data_dir)" + "if not os.path.isdir(data_dir):\n", + " raise RuntimeError('data dir', data_dir, 'does not exists')" ] }, { @@ -139,6 +171,7 @@ "metadata": {}, "outputs": [], "source": [ + "files = dict(dat='data.grp', bkg='background.grp', 'res'=resolution.grp)\n", "res = load_nexus(pjn(data_dir,'irs26173_graphite_res.nxs'))\n", "emin, emax = np.min(res['x']), np.max(res['x'])\n", "dat = load_nexus(pjn(data_dir,'irs26176_graphite002_red.nxs')) # data has 10 histograms\n", From 9c63dab96c122c0472afabeb5956d154cb8c2937 Mon Sep 17 00:00:00 2001 From: Jose Borreguero Date: Mon, 16 Apr 2018 18:54:00 -0400 Subject: [PATCH 3/5] Refs #49 Untested loading of DAVE data --- notebooks/strexp_fitting.ipynb | 45 +++++++++++++++++++++++++++++++--- 1 file changed, 41 insertions(+), 4 deletions(-) diff --git a/notebooks/strexp_fitting.ipynb b/notebooks/strexp_fitting.ipynb index 4226e4c..459c2cf 100644 --- a/notebooks/strexp_fitting.ipynb +++ b/notebooks/strexp_fitting.ipynb @@ -148,14 +148,23 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 24, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/tmp/qef_data/data/notebooks/strexp\n" + ] + } + ], "source": [ "qef_data_dir = '/tmp/qef_data'\n", "data_dir = os.path.join(qef_data_dir, 'data', 'notebooks', 'strexp')\n", "if not os.path.isdir(data_dir):\n", - " raise RuntimeError('data dir', data_dir, 'does not exists')" + " raise IOError('data dir ' + data_dir + 'does not exists')\n", + "print(data_dir)" ] }, { @@ -165,13 +174,41 @@ "(Top)

Load data and visualize data

" ] }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [], + "source": [ + "filenames = dict(dat='data.grp', bkg='background.grp', res='resolution.grp')\n", + "# Check files exists\n", + "for filename in filenames.values():\n", + " full_filename = pjn(data_dir, filename)\n", + " if not os.path.exists(full_filename):\n", + " raise IOError('File ' + filename + ' not found in ' + data_dir)\n", + " \n", + "res = load_dave(pjn(data_dir, 'resolution.grp'))\n", + "dat = load_dave(pjn(data_dir, 'data.grp'))\n", + "bkg = load_dave(pjn(data_dir, 'resolution.grp'))\n", + "\n", + "f, (ax1, ax2, ax3) = plt.subplots(1, 3)\n", + "[ax1.semilogy(res['x'], spectrum) for spectrum in res['y']]; ax1.set_title('resolution')\n", + "[ax2.semilogy(dat['x'], spectrum) for spectrum in dat['y']]; ax2.set_title('sample')\n", + "[ax3.semilogy(bkg['x'], spectrum) for spectrum in bkg['y']]; ax3.set_title('empty can')" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "files = dict(dat='data.grp', bkg='background.grp', 'res'=resolution.grp)\n", + "filenames = dict(dat='data.grp', bkg='background.grp', 'res'=resolution.grp)\n", + "# Check files exists\n", + "for filename in filenames.values():\n", + " if not os.path.isdir(os.path.join(data_dir, filename)):\n", + " raise IOError('File', filename, 'not found in', data_dir)\n", + " \n", "res = load_nexus(pjn(data_dir,'irs26173_graphite_res.nxs'))\n", "emin, emax = np.min(res['x']), np.max(res['x'])\n", "dat = load_nexus(pjn(data_dir,'irs26176_graphite002_red.nxs')) # data has 10 histograms\n", From 8a1c70e4ac4a6aaa9390b153375c10ce9775c6ea Mon Sep 17 00:00:00 2001 From: Jose Borreguero Date: Wed, 18 Apr 2018 19:53:35 -0400 Subject: [PATCH 4/5] Refs #41 load dave instead of nexus --- notebooks/strexp_fitting.ipynb | 991 +++++++++++++++++++++++++++++++-- 1 file changed, 960 insertions(+), 31 deletions(-) diff --git a/notebooks/strexp_fitting.ipynb b/notebooks/strexp_fitting.ipynb index 459c2cf..4e601dd 100644 --- a/notebooks/strexp_fitting.ipynb +++ b/notebooks/strexp_fitting.ipynb @@ -48,18 +48,9 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/sw/miniconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n", - " from ._conv import register_converters as _register_converters\n" - ] - } - ], + "outputs": [], "source": [ "from __future__ import (absolute_import, division, print_function)\n", "\n", @@ -71,7 +62,7 @@ "from lmfit.models import LinearModel, LorentzianModel, ConstantModel, LinearModel\n", "\n", "import qef\n", - "from qef.io.loaders import load_nexus\n", + "from qef.io.loaders import load_dave\n", "from qef.models.deltadirac import DeltaDiracModel\n", "from qef.models.tabulatedmodel import TabulatedModel\n", "from qef.models.resolution import TabulatedResolutionModel\n", @@ -88,7 +79,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -115,23 +106,14 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "First, rewinding head to replay your work on top of it...\n", - "Fast-forwarded master to ad9365c5454e952489223a29cacd15dc003afe2a.\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "From https://github.com/jmborr/qef_data\n", - " ee29426..ad9365c master -> origin/master\n" + "Current branch master is up to date.\n" ] } ], @@ -148,7 +130,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 5, "metadata": {}, "outputs": [ { @@ -176,9 +158,812 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 27, "metadata": {}, - "outputs": [], + "outputs": [ + { + "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.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" + }, + { + "data": { + "text/plain": [ + "Text(0.5,1,'sample')" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "e_min = -0.1 ; e_max = 0.1\n", + "\n", + "# Find indexes of dat['x'] with values in (e_min, e_max)\n", + "mask = np.intersect1d(np.where(dat['x']>e_min), np.where(dat['x']Top)

Define the model

\n", + "\n", + "
$S(E) = I \\cdot R(Q,E) \\otimes \\big( eisf \\cdot \\delta(E-E_0) + (1-eisf) \\cdot f(E-E_0, \\tau, \\beta) \\big) + ec(E, Q) + LB$
\n", + "\n", + "We will use the following component models:\n", + "- $I$, ConstantModel represents one number that can be fitted\n", + "- $R$, TabulatedResolutionModel to store the table of numbers representing the resolution\n", + "- $\\delta$, DeltaDiracModel\n", + "- $f$, StretchedExponentialFTModel has parameters $\\tau$ and $\\beta$\n", + "- $ec$, TabulatedModel to model the empty can data. One model for every Q value\n", + "- $LB$, LinearModel for the linear background" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "# Create the Model. We put everything under a function which we'll later reuse\n", + "\n", + "def generate_model_and_params(si):\n", + " r\"\"\"Produce an LMFIT model and related set of fitting parameters\n", + " @param si : spectrum index\n", + " \"\"\"\n", + "\n", + " sp = '' if si is None else '{}_'.format(si) # prefix if spectrum_index passed\n", + "\n", + " # Model components\n", + " intensity = ConstantModel(prefix='I_'+sp) # I_amplitude\n", + " elastic = DeltaDiracModel(prefix='e_'+sp) # e_amplitude, e_center\n", + " inelastic = StretchedExponentialFTModel(prefix='s_'+sp) # s_amplitude, s_center, s_tau, s_beta\n", + " resolution = TabulatedResolutionModel(res['x'], res['y'][si], prefix='r_'+sp) # (fixed r_amplitude, r_center)\n", + " emptycan = TabulatedModel(bkg['x'], bkg['y'][si], prefix='c_'+sp) # c_amplitude, c_center\n", + " background = LinearModel(prefix='b_'+sp) # b_slope, b_intercept\n", + "\n", + " # Putting it all together\n", + " model = intensity * Convolve(resolution, elastic + inelastic) + emptycan + background\n", + " parameters = model.make_params() # model parameters are a separate entity.\n", + "\n", + " # Ties and constraints\n", + " parameters['e_'+sp+'amplitude'].set(min=0.0, max=1.0)\n", + " parameters['s_'+sp+'center'].set(expr='e_'+sp+'center') # centers tied\n", + " parameters['s_'+sp+'amplitude'].set(expr='1 - e_'+sp+'amplitude')\n", + "\n", + " # Some initial sensible values\n", + " init_vals = {'I_'+sp+'c': 1.0,\n", + " 'e_'+sp+'amplitude': 0.5,\n", + " 's_'+sp+'tau': 1.0,\n", + " 's_'+sp+'beta': 1.0,\n", + " 'b_'+sp+'slope': 0,\n", + " 'b_'+sp+'intercept': 0}\n", + " for p, v in init_vals.items():\n", + " parameters[p].set(value=v)\n", + "\n", + " return model, parameters\n", + "\n", + "# Generate list of models and parameter sets for every spectrum\n", + "n_spectra = len(dat['y'])\n", + "models, paramsets = zip(*[generate_model_and_params(i) for i in range(n_spectra)])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "(Top)

Obtain an initial guess

\n", + "\n", + "A widget that compares the evaluation of the model with one of the experimental spectra. You can tweak only the free (unconstrained) parameters.\n", + "\n", + "When run, you will see two empty panels, one for comparison between experiment and model, and the second panel for residuals. Start changing the values of the parameters for the panels to be populated.\n", + "\n", + "A good initial guess for the first spectrum is :\n", + "\n", + "> spectum index = 0 \n", + "> I_c = 4 \n", + "> e_center = 0 \n", + "> e_amplitude = 0.1 \n", + "> l_sigma = 0.03 \n", + "> b_intercept = 0 \n", + "> b_slope = 0 \n" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "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.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" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "I_0_c = 0.445958372580825\n", + "r_0_amplitude = 1.0\n", + "r_0_center = 0.0\n", + "e_0_center = -1.723629445607605e-05\n", + "e_0_amplitude = 0.9582591026048619\n", + "s_0_beta = 0.3200157780798534\n", + "s_0_tau = 2163429.6396260983\n", + "s_0_center = -1.723629445607605e-05\n", + "s_0_amplitude = 0.04174089739513809\n", + "c_0_amplitude = 0.802991659199381\n", + "c_0_center = -0.00023979923706757554\n", + "b_0_intercept = -0.0010346103048792749\n", + "b_0_slope = -0.0032523714863594622\n" + ] + } + ], "source": [ "fit = model.fit(y_exp, x=e_vals, params=params, weights = 1.0 / e_exp)\n", "print('Chi-square =', fit.redchi)\n",