Merge pull request #44 from jmborr/40_simultaneous_fit

Simultaneous fitting of water data

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": [
+       "<matplotlib.figure.Figure at 0x7fc003586390>"
+      ]
+     },
+     "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
+}

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

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)

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`

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)

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)

requirements.txt

@@ -1,5 +1,6 @@
+h5py
+lmfit
+matplotlib
 numpy
 scipy
-lmfit
 six
-h5py

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
 ]
 

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))