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 +}