Repository: jmborr/nscsim

Issue:
Refactor the incoherent scattering to emulate the coherent scattering calculations
- `incoherent2.py` will host the new incoherent function
- `incoherent2.py` tries to follow the structure of `coherent.py`

Discussion:

Patch:
diff --git a/.travis.yml b/.travis.yml
index ab323bd..dfa757e 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,7 +1,7 @@
 language: python
-
+xdist: xenial
 python:
-  - "3.4"
+  - "3.5"
   - "3.6"
 
 env: PYTHONPATH=$PYTHONPATH:$TRAVIS_BUILD_DIR/tests
@@ -23,10 +23,8 @@ install:
   - pip install -r requirements_dev.txt
 
 # Set display for tests involving matplotlib
-before_script: # configure a headless display to test plot generation
-- "export DISPLAY=:99.0"
-- "sh -e /etc/init.d/xvfb start"
-- sleep 3 # give xvfb some time to start
+services:
+- xvfb
 
 script:
   - pytest --cov=nscsim tests/unit
@@ -37,4 +35,3 @@ script:
 notifications:
   - on_success: never
   - on_failure: never
-
diff --git a/notebooks/quasielastic/README.txt b/notebooks/quasielastic/README.txt
new file mode 100644
index 0000000..fff6a2f
--- /dev/null
+++ b/notebooks/quasielastic/README.txt
@@ -0,0 +1,6 @@
+Water box of 30 Angstroms with Lithium and Choride atoms dissolved in it.
+- 154 Lithium (LIT) atoms
+- 154 Chloride (CLA) atoms
+- 1128 Water molecules
+
+Frames are separated by 1ps, for a total of 4ns (4000 frames)
diff --git a/notebooks/quasielastic/incoherent.ipynb b/notebooks/quasielastic/incoherent.ipynb
new file mode 100644
index 0000000..193f83c
--- /dev/null
+++ b/notebooks/quasielastic/incoherent.ipynb
@@ -0,0 +1,1235 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Requirements\n",
+    "\n",
+    "#### Python packages\n",
+    "\n",
+    "The following packages are not required by <i>nscsim</i> but are needed to run this notebook\n",
+    "\n",
+    "- [cython](https://pypi.org/project/Cython/)\n",
+    "- [MDAnalysis](https://pypi.org/project/MDAnalysis/)\n",
+    "- [periodictable](https://pypi.org/project/periodictable/)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Requirements"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Does everything look good? Mark the review as Approve. A member of @nscsim/gatekeepers will take care of it.\n",
+    "\n",
+    "%matplotlib notebook\n",
+    "import matplotlib.pyplot as plt\n",
+    "plt.ion()\n",
+    "import numpy as np\n",
+    "import MDAnalysis as mda\n",
+    "import periodictable\n",
+    "import nscsim\n",
+    "from nscsim.quasielastic import incoherent as incoh\n",
+    "from nscsim.qvec import moduli_linscale\n",
+    "from nscsim import plotting"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Download water simulation\n",
+    "`wget` and `tar` commands are necessary.  \n",
+    "You may skip this step if running the notebook more than once"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%capture\n",
+    "%%bash\n",
+    "wget https://www.dropbox.com/s/7imfnz6ma1kh41w/water_4000.tar.gz\n",
+    "tar zxf water_4000.tar.gz"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Create system and extract coordinates\n",
+    "\n",
+    "`u` is a `MDAnalysis.Universe` object. It contains 154 Lithium atoms, 154 Chloride atoms, and 1128 TIP3 water molecules."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "u = mda.Universe('topology.psf', 'trajectory.dcd')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We investigate the trajectory."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "number of conformations = 4000\n",
+      "Size of the box, in Angstroms = [33.6035 33.6035 33.6035]\n",
+      "It is a cubic box, as evidenced by the angles between box sides:  [90. 90. 90.]\n"
+     ]
+    }
+   ],
+   "source": [
+    "print('number of conformations =', u.trajectory.n_frames)\n",
+    "print('Size of the box, in Angstroms =', u.dimensions[0:3])\n",
+    "print('It is a cubic box, as evidenced by the angles between box sides: ',u.dimensions[3:])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Select only the hydrogens and oxygens. There should be three atoms per water molecule. A look at file <i>topology.psf</i> shows that water hydrogens have type `HT` and water oxygens have type `OT`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "number of atoms selected = 3384\n",
+      "number of atoms per water molecule =  3.0\n"
+     ]
+    }
+   ],
+   "source": [
+    "sel = u.select_atoms('type HT or type OT')\n",
+    "print('number of atoms selected =', len(sel))\n",
+    "print('number of atoms per water molecule = ', len(sel)/len(sel.residues))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Obtain the bound incoherent neutron scattering cross sections for the atoms of our selection `sel`. Notice we have to \"translate\" from NAMD atom types to the name of the atomic element. We use dictionary `type_to_element` for this. Then we iterate over our selected atoms"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Calculate the incoherent dynamic structure factor `S(q, t)`\n",
+    "\n",
+    "The liquid allows to compute scattering only as a function of Q-vector modulus so we carry out a spherical average for all orientations of vector Q (function `intermediate_spherical`)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "q_mod = moduli_linscale(0.1, q_mod_max=2.1, q_mod_delta=0.2)\n",
+    "# Extract coordinates to numpy array with shape=(#atoms, #frames, 3)\n",
+    "tr = np.asarray([sel.positions for _ in u.trajectory[0:1000:1]])  # shape=(#frames, #atoms, #3)\n",
+    "tr = tr.transpose((1,0,2))  # shape=(#atoms, #frames, 3)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[ 0.   80.26 80.26  0.   80.26 80.26  0.   80.26 80.26]\n"
+     ]
+    }
+   ],
+   "source": [
+    "type_to_element = dict(HT='H', OT='O')\n",
+    "xi = np.asarray([getattr(periodictable, type_to_element[at.type]).neutron.incoherent for at in sel])\n",
+    "print(xi[0:9])\n",
+    "bi = np.sqrt(xi)  # incoherent.intermediate_spherical module requires bi instead of xi"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "results = incoh.self_intermediate_spherical(q_mod, tr, bi, nt=500, n_cores=8)\n",
+    "sf, times = results.sf, results.t"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import nscsim.quasielastic.incoherent2 as inc2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 47,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<module 'nscsim.quasielastic.incoherent2' from '/home/39p/projects/nscsim/nscsim/quasielastic/incoherent2.py'>"
+      ]
+     },
+     "execution_count": 47,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "import importlib\n",
+    "importlib.reload(inc2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 46,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "Calculating coherent scattering for a set of q-moduli\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "\n",
+      "  0%|          | 0/10 [00:00<?, ?it/s]\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\u001b[A\n",
+      "Calculating incoherent intensities for one set of q vectors\n",
+      "\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "> /home/39p/projects/nscsim/nscsim/quasielastic/incoherent2.py(70)intermediate_vector_set()\n",
+      "-> result = np.sum(result)\n",
+      "(Pdb) obj = next(result)\n",
+      "*** ValueError: object too deep for desired array\n",
+      "(Pdb) serial_worker(tr[0], s_inc[0])\n",
+      "*** ValueError: object too deep for desired array\n",
+      "(Pdb) tr[0]\n",
+      "array([[ 15.827579 , -15.575326 ,   9.810642 ],\n",
+      "       [ 15.530334 , -15.732183 ,   9.759977 ],\n",
+      "       [ 15.9462805, -15.673354 ,  10.000448 ],\n",
+      "       ...,\n",
+      "       [ 13.059978 , -21.137667 ,  19.665363 ],\n",
+      "       [ 13.388484 , -21.850504 ,  19.279696 ],\n",
+      "       [ 14.395192 , -20.926067 ,  19.130993 ]], dtype=float32)\n",
+      "(Pdb) s_inc[0]\n",
+      "2.4093567606313515\n",
+      "(Pdb) serial_worker\n",
+      "<function intermediate_vector_set.<locals>.serial_worker at 0x7f5e484146a8>\n",
+      "(Pdb) ai = np.exp(-1j*np.tensordot(q, tr[0],axes=(1,1)))  \n",
+      "(Pdb) ai\n",
+      "array([[-0.410982  -0.91164346j, -0.40732624-0.91328272j,\n",
+      "        -0.43242458-0.90167011j, ..., -0.9922127 -0.12455505j,\n",
+      "        -0.99363402-0.11265626j, -0.98944188-0.14493019j],\n",
+      "       [-0.50359873-0.86393768j, -0.50812092-0.86128574j,\n",
+      "        -0.52079445-0.85368211j, ..., -0.99196758-0.12649239j,\n",
+      "        -0.99678871-0.08007669j, -0.98917699-0.14672724j],\n",
+      "       [ 0.09221902+0.99573875j,  0.12205932+0.99252281j,\n",
+      "         0.07428757+0.99723686j, ..., -0.11670681+0.99316641j,\n",
+      "        -0.1060775 +0.99435787j, -0.19722034+0.98035919j],\n",
+      "       ...,\n",
+      "       [ 0.60661421-0.79499635j,  0.63276009-0.77434789j,\n",
+      "         0.60065824-0.79950589j, ...,  0.89814757-0.43969414j,\n",
+      "         0.89563129-0.44479725j,  0.83205466-0.55469365j],\n",
+      "       [-0.54920121-0.83569015j, -0.52884577-0.84871795j,\n",
+      "        -0.5670394 -0.82369067j, ..., -0.81694438-0.57671647j,\n",
+      "        -0.83044394-0.55710221j, -0.86111353-0.50841272j],\n",
+      "       [ 0.55060688+0.83476468j,  0.54067208+0.84123344j,\n",
+      "         0.5367955 +0.84371239j, ..., -0.54499735+0.83843777j,\n",
+      "        -0.53415565+0.84538615j, -0.46135365+0.88721632j]])\n",
+      "(Pdb) s_inc[0] * np.correlate(ai, ai, 'full')\n",
+      "*** ValueError: object too deep for desired array\n",
+      "(Pdb) ai.shape\n",
+      "(200, 1000)\n",
+      "(Pdb) tr.shape\n",
+      "(3384, 1000, 3)\n",
+      "(Pdb) q.shape\n"
+     ]
+    },
+    {
+     "ename": "KeyboardInterrupt",
+     "evalue": "",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mBdbQuit\u001b[0m                                   Traceback (most recent call last)",
+      "\u001b[0;32m~/projects/nscsim/nscsim/quasielastic/incoherent2.py\u001b[0m in \u001b[0;36mintermediate_vector_set\u001b[0;34m(tr, q, s_inc, n_cores)\u001b[0m\n\u001b[1;32m     69\u001b[0m         \u001b[0mbreakpoint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 70\u001b[0;31m         \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     71\u001b[0m     \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m~/projects/nscsim/nscsim/quasielastic/incoherent2.py\u001b[0m in \u001b[0;36mintermediate_vector_set\u001b[0;34m(tr, q, s_inc, n_cores)\u001b[0m\n\u001b[1;32m     69\u001b[0m         \u001b[0mbreakpoint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 70\u001b[0;31m         \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     71\u001b[0m     \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m~/anaconda3/envs/nscsim/lib/python3.7/bdb.py\u001b[0m in \u001b[0;36mtrace_dispatch\u001b[0;34m(self, frame, event, arg)\u001b[0m\n\u001b[1;32m     87\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mevent\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'line'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 88\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdispatch_line\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mframe\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     89\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mevent\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'call'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m~/anaconda3/envs/nscsim/lib/python3.7/bdb.py\u001b[0m in \u001b[0;36mdispatch_line\u001b[0;34m(self, frame)\u001b[0m\n\u001b[1;32m    111\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstop_here\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mframe\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbreak_here\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mframe\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 112\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0muser_line\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mframe\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    113\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mquitting\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mBdbQuit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;31mBdbQuit\u001b[0m: ",
+      "\nDuring handling of the above exception, another exception occurred:\n",
+      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
+      "\u001b[0;32m<ipython-input-46-97aec4de631f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0msf2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minc2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mintermediate_spherical\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq_mod\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_cores\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
+      "\u001b[0;32m~/projects/nscsim/nscsim/quasielastic/incoherent2.py\u001b[0m in \u001b[0;36mintermediate_spherical\u001b[0;34m(tr, q_mod, s_inc, n_cores)\u001b[0m\n\u001b[1;32m    136\u001b[0m     sf = [\n\u001b[1;32m    137\u001b[0m         \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mintermediate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq_set\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ms_inc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_cores\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mn_cores\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 138\u001b[0;31m         \u001b[0;32mfor\u001b[0m \u001b[0mq_set\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mtqdm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mq_sets\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    139\u001b[0m     ]\n\u001b[1;32m    140\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m~/projects/nscsim/nscsim/quasielastic/incoherent2.py\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m    136\u001b[0m     sf = [\n\u001b[1;32m    137\u001b[0m         \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mintermediate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq_set\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ms_inc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_cores\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mn_cores\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 138\u001b[0;31m         \u001b[0;32mfor\u001b[0m \u001b[0mq_set\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mtqdm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mq_sets\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    139\u001b[0m     ]\n\u001b[1;32m    140\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m~/projects/nscsim/nscsim/quasielastic/incoherent2.py\u001b[0m in \u001b[0;36mintermediate\u001b[0;34m(tr, q, s_inc, n_cores, averaging)\u001b[0m\n\u001b[1;32m    107\u001b[0m         shape = (#frames)\n\u001b[1;32m    108\u001b[0m     \"\"\"\n\u001b[0;32m--> 109\u001b[0;31m     \u001b[0msf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mintermediate_vector_set\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ms_inc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_cores\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mn_cores\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    110\u001b[0m     \u001b[0mav_args\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0maveraging\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0maveraging\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    111\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0maveraging\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mav_args\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0maveraging\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m~/projects/nscsim/nscsim/quasielastic/incoherent2.py\u001b[0m in \u001b[0;36mintermediate_vector_set\u001b[0;34m(tr, q, s_inc, n_cores)\u001b[0m\n\u001b[1;32m     71\u001b[0m     \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     72\u001b[0m         \u001b[0mpool\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 73\u001b[0;31m         \u001b[0mpool\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     74\u001b[0m         \u001b[0mpool\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrestart\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     75\u001b[0m         \u001b[0mpool\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclear\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m~/anaconda3/envs/nscsim/lib/python3.7/site-packages/pathos/multiprocessing.py\u001b[0m in \u001b[0;36mjoin\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    206\u001b[0m         \u001b[0m_pool\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m__STATE\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_id\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    207\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0m_pool\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__nodes\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0m_pool\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__nodes\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 208\u001b[0;31m             \u001b[0m_pool\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    209\u001b[0m         \u001b[0;32mreturn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    210\u001b[0m     \u001b[0;31m# interface\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m~/anaconda3/envs/nscsim/lib/python3.7/site-packages/multiprocess/pool.py\u001b[0m in \u001b[0;36mjoin\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    554\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_state\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32min\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mCLOSE\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mTERMINATE\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    555\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"In unknown state\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 556\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_worker_handler\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    557\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_task_handler\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    558\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_result_handler\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m~/anaconda3/envs/nscsim/lib/python3.7/threading.py\u001b[0m in \u001b[0;36mjoin\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m   1030\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1031\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mtimeout\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1032\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_wait_for_tstate_lock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1033\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1034\u001b[0m             \u001b[0;31m# the behavior of a negative timeout isn't documented, but\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m~/anaconda3/envs/nscsim/lib/python3.7/threading.py\u001b[0m in \u001b[0;36m_wait_for_tstate_lock\u001b[0;34m(self, block, timeout)\u001b[0m\n\u001b[1;32m   1046\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mlock\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# already determined that the C code is done\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1047\u001b[0m             \u001b[0;32massert\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_is_stopped\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1048\u001b[0;31m         \u001b[0;32melif\u001b[0m \u001b[0mlock\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0macquire\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mblock\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1049\u001b[0m             \u001b[0mlock\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrelease\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1050\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_stop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
+     ]
+    }
+   ],
+   "source": [
+    "sf2 = inc2.intermediate_spherical(tr, q_mod, bi, n_cores=8)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 50,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 51,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.save('tr_ref.npy', tr)\n",
+    "np.save('q_mod_ref.npy', q_mod)\n",
+    "np.save('bi_ref.npy', bi)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Simple plotting of the curve."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "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 = $('<div/>');\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",
+       "        '<div class=\"ui-dialog-titlebar ui-widget-header ui-corner-all ' +\n",
+       "        'ui-helper-clearfix\"/>');\n",
+       "    var titletext = $(\n",
+       "        '<div class=\"ui-dialog-title\" style=\"width: 100%; ' +\n",
+       "        'text-align: center; padding: 3px;\"/>');\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 = $('<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 = $('<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 = $('<canvas/>');\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 = $('<div/>');\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 = $('<button/>');\n",
+       "        button.addClass('ui-button ui-widget ui-state-default ui-corner-all ' +\n",
+       "                        'ui-button-icon-only');\n",
+       "        button.attr('role', 'button');\n",
+       "        button.attr('aria-disabled', 'false');\n",
+       "        button.click(method_name, toolbar_event);\n",
+       "        button.mouseover(tooltip, toolbar_mouse_event);\n",
+       "\n",
+       "        var icon_img = $('<span/>');\n",
+       "        icon_img.addClass('ui-button-icon-primary ui-icon');\n",
+       "        icon_img.addClass(image);\n",
+       "        icon_img.addClass('ui-corner-all');\n",
+       "\n",
+       "        var tooltip_span = $('<span/>');\n",
+       "        tooltip_span.addClass('ui-button-text');\n",
+       "        tooltip_span.html(tooltip);\n",
+       "\n",
+       "        button.append(icon_img);\n",
+       "        button.append(tooltip_span);\n",
+       "\n",
+       "        nav_element.append(button);\n",
+       "    }\n",
+       "\n",
+       "    var fmt_picker_span = $('<span/>');\n",
+       "\n",
+       "    var fmt_picker = $('<select/>');\n",
+       "    fmt_picker.addClass('mpl-toolbar-option ui-widget ui-widget-content');\n",
+       "    fmt_picker_span.append(fmt_picker);\n",
+       "    nav_element.append(fmt_picker_span);\n",
+       "    this.format_dropdown = fmt_picker[0];\n",
+       "\n",
+       "    for (var ind in mpl.extensions) {\n",
+       "        var fmt = mpl.extensions[ind];\n",
+       "        var option = $(\n",
+       "            '<option/>', {selected: fmt === mpl.default_extension}).html(fmt);\n",
+       "        fmt_picker.append(option);\n",
+       "    }\n",
+       "\n",
+       "    // Add hover states to the ui-buttons\n",
+       "    $( \".ui-button\" ).hover(\n",
+       "        function() { $(this).addClass(\"ui-state-hover\");},\n",
+       "        function() { $(this).removeClass(\"ui-state-hover\");}\n",
+       "    );\n",
+       "\n",
+       "    var status_bar = $('<span class=\"mpl-message\"/>');\n",
+       "    nav_element.append(status_bar);\n",
+       "    this.message = status_bar[0];\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.request_resize = function(x_pixels, y_pixels) {\n",
+       "    // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n",
+       "    // which will in turn request a refresh of the image.\n",
+       "    this.send_message('resize', {'width': x_pixels, 'height': y_pixels});\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.send_message = function(type, properties) {\n",
+       "    properties['type'] = type;\n",
+       "    properties['figure_id'] = this.id;\n",
+       "    this.ws.send(JSON.stringify(properties));\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.send_draw_message = function() {\n",
+       "    if (!this.waiting) {\n",
+       "        this.waiting = true;\n",
+       "        this.ws.send(JSON.stringify({type: \"draw\", figure_id: this.id}));\n",
+       "    }\n",
+       "}\n",
+       "\n",
+       "\n",
+       "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
+       "    var format_dropdown = fig.format_dropdown;\n",
+       "    var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n",
+       "    fig.ondownload(fig, format);\n",
+       "}\n",
+       "\n",
+       "\n",
+       "mpl.figure.prototype.handle_resize = function(fig, msg) {\n",
+       "    var size = msg['size'];\n",
+       "    if (size[0] != fig.canvas.width || size[1] != fig.canvas.height) {\n",
+       "        fig._resize_canvas(size[0], size[1]);\n",
+       "        fig.send_message(\"refresh\", {});\n",
+       "    };\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.handle_rubberband = function(fig, msg) {\n",
+       "    var x0 = msg['x0'] / mpl.ratio;\n",
+       "    var y0 = (fig.canvas.height - msg['y0']) / mpl.ratio;\n",
+       "    var x1 = msg['x1'] / mpl.ratio;\n",
+       "    var y1 = (fig.canvas.height - msg['y1']) / mpl.ratio;\n",
+       "    x0 = Math.floor(x0) + 0.5;\n",
+       "    y0 = Math.floor(y0) + 0.5;\n",
+       "    x1 = Math.floor(x1) + 0.5;\n",
+       "    y1 = Math.floor(y1) + 0.5;\n",
+       "    var min_x = Math.min(x0, x1);\n",
+       "    var min_y = Math.min(y0, y1);\n",
+       "    var width = Math.abs(x1 - x0);\n",
+       "    var height = Math.abs(y1 - y0);\n",
+       "\n",
+       "    fig.rubberband_context.clearRect(\n",
+       "        0, 0, fig.canvas.width / mpl.ratio, fig.canvas.height / mpl.ratio);\n",
+       "\n",
+       "    fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.handle_figure_label = function(fig, msg) {\n",
+       "    // Updates the figure title.\n",
+       "    fig.header.textContent = msg['label'];\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.handle_cursor = function(fig, msg) {\n",
+       "    var cursor = msg['cursor'];\n",
+       "    switch(cursor)\n",
+       "    {\n",
+       "    case 0:\n",
+       "        cursor = 'pointer';\n",
+       "        break;\n",
+       "    case 1:\n",
+       "        cursor = 'default';\n",
+       "        break;\n",
+       "    case 2:\n",
+       "        cursor = 'crosshair';\n",
+       "        break;\n",
+       "    case 3:\n",
+       "        cursor = 'move';\n",
+       "        break;\n",
+       "    }\n",
+       "    fig.rubberband_canvas.style.cursor = cursor;\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.handle_message = function(fig, msg) {\n",
+       "    fig.message.textContent = msg['message'];\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.handle_draw = function(fig, msg) {\n",
+       "    // Request the server to send over a new figure.\n",
+       "    fig.send_draw_message();\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.handle_image_mode = function(fig, msg) {\n",
+       "    fig.image_mode = msg['mode'];\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.updated_canvas_event = function() {\n",
+       "    // Called whenever the canvas gets updated.\n",
+       "    this.send_message(\"ack\", {});\n",
+       "}\n",
+       "\n",
+       "// A function to construct a web socket function for onmessage handling.\n",
+       "// Called in the figure constructor.\n",
+       "mpl.figure.prototype._make_on_message_function = function(fig) {\n",
+       "    return function socket_on_message(evt) {\n",
+       "        if (evt.data instanceof Blob) {\n",
+       "            /* FIXME: We get \"Resource interpreted as Image but\n",
+       "             * transferred with MIME type text/plain:\" errors on\n",
+       "             * Chrome.  But how to set the MIME type?  It doesn't seem\n",
+       "             * to be part of the websocket stream */\n",
+       "            evt.data.type = \"image/png\";\n",
+       "\n",
+       "            /* Free the memory for the previous frames */\n",
+       "            if (fig.imageObj.src) {\n",
+       "                (window.URL || window.webkitURL).revokeObjectURL(\n",
+       "                    fig.imageObj.src);\n",
+       "            }\n",
+       "\n",
+       "            fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n",
+       "                evt.data);\n",
+       "            fig.updated_canvas_event();\n",
+       "            fig.waiting = false;\n",
+       "            return;\n",
+       "        }\n",
+       "        else if (typeof evt.data === 'string' && evt.data.slice(0, 21) == \"data:image/png;base64\") {\n",
+       "            fig.imageObj.src = evt.data;\n",
+       "            fig.updated_canvas_event();\n",
+       "            fig.waiting = false;\n",
+       "            return;\n",
+       "        }\n",
+       "\n",
+       "        var msg = JSON.parse(evt.data);\n",
+       "        var msg_type = msg['type'];\n",
+       "\n",
+       "        // Call the  \"handle_{type}\" callback, which takes\n",
+       "        // the figure and JSON message as its only arguments.\n",
+       "        try {\n",
+       "            var callback = fig[\"handle_\" + msg_type];\n",
+       "        } catch (e) {\n",
+       "            console.log(\"No handler for the '\" + msg_type + \"' message type: \", msg);\n",
+       "            return;\n",
+       "        }\n",
+       "\n",
+       "        if (callback) {\n",
+       "            try {\n",
+       "                // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n",
+       "                callback(fig, msg);\n",
+       "            } catch (e) {\n",
+       "                console.log(\"Exception inside the 'handler_\" + msg_type + \"' callback:\", e, e.stack, msg);\n",
+       "            }\n",
+       "        }\n",
+       "    };\n",
+       "}\n",
+       "\n",
+       "// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n",
+       "mpl.findpos = function(e) {\n",
+       "    //this section is from http://www.quirksmode.org/js/events_properties.html\n",
+       "    var targ;\n",
+       "    if (!e)\n",
+       "        e = window.event;\n",
+       "    if (e.target)\n",
+       "        targ = e.target;\n",
+       "    else if (e.srcElement)\n",
+       "        targ = e.srcElement;\n",
+       "    if (targ.nodeType == 3) // defeat Safari bug\n",
+       "        targ = targ.parentNode;\n",
+       "\n",
+       "    // jQuery normalizes the pageX and pageY\n",
+       "    // pageX,Y are the mouse positions relative to the document\n",
+       "    // offset() returns the position of the element relative to the document\n",
+       "    var x = e.pageX - $(targ).offset().left;\n",
+       "    var y = e.pageY - $(targ).offset().top;\n",
+       "\n",
+       "    return {\"x\": x, \"y\": y};\n",
+       "};\n",
+       "\n",
+       "/*\n",
+       " * return a copy of an object with only non-object keys\n",
+       " * we need this to avoid circular references\n",
+       " * http://stackoverflow.com/a/24161582/3208463\n",
+       " */\n",
+       "function simpleKeys (original) {\n",
+       "  return Object.keys(original).reduce(function (obj, key) {\n",
+       "    if (typeof original[key] !== 'object')\n",
+       "        obj[key] = original[key]\n",
+       "    return obj;\n",
+       "  }, {});\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.mouse_event = function(event, name) {\n",
+       "    var canvas_pos = mpl.findpos(event)\n",
+       "\n",
+       "    if (name === 'button_press')\n",
+       "    {\n",
+       "        this.canvas.focus();\n",
+       "        this.canvas_div.focus();\n",
+       "    }\n",
+       "\n",
+       "    var x = canvas_pos.x * mpl.ratio;\n",
+       "    var y = canvas_pos.y * mpl.ratio;\n",
+       "\n",
+       "    this.send_message(name, {x: x, y: y, button: event.button,\n",
+       "                             step: event.step,\n",
+       "                             guiEvent: simpleKeys(event)});\n",
+       "\n",
+       "    /* This prevents the web browser from automatically changing to\n",
+       "     * the text insertion cursor when the button is pressed.  We want\n",
+       "     * to control all of the cursor setting manually through the\n",
+       "     * 'cursor' event from matplotlib */\n",
+       "    event.preventDefault();\n",
+       "    return false;\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
+       "    // Handle any extra behaviour associated with a key event\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.key_event = function(event, name) {\n",
+       "\n",
+       "    // Prevent repeat events\n",
+       "    if (name == 'key_press')\n",
+       "    {\n",
+       "        if (event.which === this._key)\n",
+       "            return;\n",
+       "        else\n",
+       "            this._key = event.which;\n",
+       "    }\n",
+       "    if (name == 'key_release')\n",
+       "        this._key = null;\n",
+       "\n",
+       "    var value = '';\n",
+       "    if (event.ctrlKey && event.which != 17)\n",
+       "        value += \"ctrl+\";\n",
+       "    if (event.altKey && event.which != 18)\n",
+       "        value += \"alt+\";\n",
+       "    if (event.shiftKey && event.which != 16)\n",
+       "        value += \"shift+\";\n",
+       "\n",
+       "    value += 'k';\n",
+       "    value += event.which.toString();\n",
+       "\n",
+       "    this._key_event_extra(event, name);\n",
+       "\n",
+       "    this.send_message(name, {key: value,\n",
+       "                             guiEvent: simpleKeys(event)});\n",
+       "    return false;\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.toolbar_button_onclick = function(name) {\n",
+       "    if (name == 'download') {\n",
+       "        this.handle_save(this, null);\n",
+       "    } else {\n",
+       "        this.send_message(\"toolbar_button\", {name: name});\n",
+       "    }\n",
+       "};\n",
+       "\n",
+       "mpl.figure.prototype.toolbar_button_onmouseover = function(tooltip) {\n",
+       "    this.message.textContent = tooltip;\n",
+       "};\n",
+       "mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Pan axes with left mouse, zoom with right\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n",
+       "\n",
+       "mpl.extensions = [\"eps\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\"];\n",
+       "\n",
+       "mpl.default_extension = \"png\";var comm_websocket_adapter = function(comm) {\n",
+       "    // Create a \"websocket\"-like object which calls the given IPython comm\n",
+       "    // object with the appropriate methods. Currently this is a non binary\n",
+       "    // socket, so there is still some room for performance tuning.\n",
+       "    var ws = {};\n",
+       "\n",
+       "    ws.close = function() {\n",
+       "        comm.close()\n",
+       "    };\n",
+       "    ws.send = function(m) {\n",
+       "        //console.log('sending', m);\n",
+       "        comm.send(m);\n",
+       "    };\n",
+       "    // Register the callback with on_msg.\n",
+       "    comm.on_msg(function(msg) {\n",
+       "        //console.log('receiving', msg['content']['data'], msg);\n",
+       "        // Pass the mpl event to the overridden (by mpl) onmessage function.\n",
+       "        ws.onmessage(msg['content']['data'])\n",
+       "    });\n",
+       "    return ws;\n",
+       "}\n",
+       "\n",
+       "mpl.mpl_figure_comm = function(comm, msg) {\n",
+       "    // This is the function which gets called when the mpl process\n",
+       "    // starts-up an IPython Comm through the \"matplotlib\" channel.\n",
+       "\n",
+       "    var id = msg.content.data.id;\n",
+       "    // Get hold of the div created by the display call when the Comm\n",
+       "    // socket was opened in Python.\n",
+       "    var element = $(\"#\" + id);\n",
+       "    var ws_proxy = comm_websocket_adapter(comm)\n",
+       "\n",
+       "    function ondownload(figure, format) {\n",
+       "        window.open(figure.imageObj.src);\n",
+       "    }\n",
+       "\n",
+       "    var fig = new mpl.figure(id, ws_proxy,\n",
+       "                           ondownload,\n",
+       "                           element.get(0));\n",
+       "\n",
+       "    // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n",
+       "    // web socket which is closed, not our websocket->open comm proxy.\n",
+       "    ws_proxy.onopen();\n",
+       "\n",
+       "    fig.parent_element = element.get(0);\n",
+       "    fig.cell_info = mpl.find_output_cell(\"<div id='\" + id + \"'></div>\");\n",
+       "    if (!fig.cell_info) {\n",
+       "        console.error(\"Failed to find cell for figure\", id, fig);\n",
+       "        return;\n",
+       "    }\n",
+       "\n",
+       "    var output_index = fig.cell_info[2]\n",
+       "    var cell = fig.cell_info[0];\n",
+       "\n",
+       "};\n",
+       "\n",
+       "mpl.figure.prototype.handle_close = function(fig, msg) {\n",
+       "    var width = fig.canvas.width/mpl.ratio\n",
+       "    fig.root.unbind('remove')\n",
+       "\n",
+       "    // Update the output cell to use the data from the current canvas.\n",
+       "    fig.push_to_output();\n",
+       "    var dataURL = fig.canvas.toDataURL();\n",
+       "    // Re-enable the keyboard manager in IPython - without this line, in FF,\n",
+       "    // the notebook keyboard shortcuts fail.\n",
+       "    IPython.keyboard_manager.enable()\n",
+       "    $(fig.parent_element).html('<img src=\"' + dataURL + '\" width=\"' + width + '\">');\n",
+       "    fig.close_ws(fig, msg);\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.close_ws = function(fig, msg){\n",
+       "    fig.send_message('closing', msg);\n",
+       "    // fig.ws.close()\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.push_to_output = function(remove_interactive) {\n",
+       "    // Turn the data on the canvas into data in the output cell.\n",
+       "    var width = this.canvas.width/mpl.ratio\n",
+       "    var dataURL = this.canvas.toDataURL();\n",
+       "    this.cell_info[1]['text/html'] = '<img src=\"' + dataURL + '\" width=\"' + width + '\">';\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype.updated_canvas_event = function() {\n",
+       "    // Tell IPython that the notebook contents must change.\n",
+       "    IPython.notebook.set_dirty(true);\n",
+       "    this.send_message(\"ack\", {});\n",
+       "    var fig = this;\n",
+       "    // Wait a second, then push the new image to the DOM so\n",
+       "    // that it is saved nicely (might be nice to debounce this).\n",
+       "    setTimeout(function () { fig.push_to_output() }, 1000);\n",
+       "}\n",
+       "\n",
+       "mpl.figure.prototype._init_toolbar = function() {\n",
+       "    var fig = this;\n",
+       "\n",
+       "    var nav_element = $('<div/>');\n",
+       "    nav_element.attr('style', 'width: 100%');\n",
+       "    this.root.append(nav_element);\n",
+       "\n",
+       "    // Define a callback function for later on.\n",
+       "    function toolbar_event(event) {\n",
+       "        return fig.toolbar_button_onclick(event['data']);\n",
+       "    }\n",
+       "    function toolbar_mouse_event(event) {\n",
+       "        return fig.toolbar_button_onmouseover(event['data']);\n",
+       "    }\n",
+       "\n",
+       "    for(var toolbar_ind in mpl.toolbar_items){\n",
+       "        var name = mpl.toolbar_items[toolbar_ind][0];\n",
+       "        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
+       "        var image = mpl.toolbar_items[toolbar_ind][2];\n",
+       "        var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
+       "\n",
+       "        if (!name) { continue; };\n",
+       "\n",
+       "        var button = $('<button class=\"btn btn-default\" href=\"#\" title=\"' + name + '\"><i class=\"fa ' + image + ' fa-lg\"></i></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 = $('<span class=\"mpl-message\" style=\"text-align:right; float: right;\"/>');\n",
+       "    nav_element.append(status_bar);\n",
+       "    this.message = status_bar[0];\n",
+       "\n",
+       "    // Add the close button to the window.\n",
+       "    var buttongrp = $('<div class=\"btn-group inline pull-right\"></div>');\n",
+       "    var button = $('<button class=\"btn btn-mini btn-primary\" href=\"#\" title=\"Stop Interaction\"><i class=\"fa fa-power-off icon-remove icon-large\"></i></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<ncells; i++) {\n",
+       "        var cell = cells[i];\n",
+       "        if (cell.cell_type === 'code'){\n",
+       "            for (var j=0; j<cell.output_area.outputs.length; j++) {\n",
+       "                var data = cell.output_area.outputs[j];\n",
+       "                if (data.data) {\n",
+       "                    // IPython >= 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": [
+       "<IPython.core.display.Javascript object>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/html": [
+       "<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAyAAAAJYCAYAAACadoJwAAAgAElEQVR4XuzdB3RURdsH8D/ppBAg9N47CSEURT5REBERlaY0AZX+AjYUBRVFUUEFpCggKkgHwYYCIqIiRSD0LgQpUkInEELqd56JuwZI2U1md+/e/c8571HDvVN+M/rOk7kzky8tLS0NTBSgAAUoQAEKUIACFKAABZwgkI8BiBOUWQQFKEABClCAAhSgAAUooAQYgHAgUIACFKAABShAAQpQgAJOE2AA4jRqFkQBClCAAhSgAAUoQAEKMADhGKAABShAAQpQgAIUoAAFnCbAAMRp1CyIAhSgAAUoQAEKUIACFGAAwjFAAQpQgAIUoAAFKEABCjhNgAGI06hZEAUoQAEKUIACFKAABSjAAIRjgAIUoAAFKEABClCAAhRwmgADEKdRsyAKUIACFKAABShAAQpQgAEIxwAFKEABClCAAhSgAAUo4DQBBiBOo2ZBFKAABShAAQpQgAIUoAADEI4BClCAAhSgAAUoQAEKUMBpAgxAnEbNgihAAQpQgAIUoAAFKEABBiAcAxSgAAUoQAEKUIACFKCA0wQYgDiNmgVRgAIUoAAFKEABClCAAgxAOAYoQAEKUIACFKAABShAAacJMABxGjULogAFKEABClCAAhSgAAUYgJh0DJw7dw4rV65EhQoVkD9/fpO2ks2iAAUoQAEKUIACFHC0wPXr1/H333+jVatWKFKkSJ6LYwCSZ0JjZjB37lx0797dmJVjrShAAQpQgAIUoAAF3E5gzpw56NatW57rzQAkz4TGzGDdunVo2rQpZKDUrFnTmJVkrShAAQpQgAIUoAAFDC+wb98+9YvtP/74A3fddVee68sAJM+Exsxg69atiIqKQnR0NOrXr2/MSrJWFKAABShAAQpQgAKGF9A9r2QAYvguz10FdQ+U3NWCb1GAAhSgAAUoQAEKuLuA7nklAxB3HxFZ1F/3QDEpE5tFAQpQgAIUoAAFKJCDgO55JQMQkw453QPFpExsFgUoQAEKUIACFKAAAxCOAR0CDEB0KDIPClCAAhSgAAUoQAHd80qugJh0TOkeKCZlYrMoQAEKUIACFHCyQFJSEi5evIirV68iNTXVyaWzOItAvnz54O/vjwIFCiAkJATyz1kl3fNKBiAmHYe6B4pJmdgsClCAAhSgAAWcKJCWlqYutEtISICXlxe8vb2dWDqLyiggwV9KSor6UcGCBVGiRIksgxDd80oGICYdi7oHikmZ2CwKUIACFKAABZwocOnSJZw6dUr91r1UqVLZ/tbdidXy2KJu3LiBkydPqoCwdOnSql8yS7rnlQxATDrkdA8UkzKxWRSgAAUoQAEKOFHg2LFjuHbtGqpUqQJfX18nlsyishKQICQmJkZ9hlWmTBkGIBwquRdgAJJ7O75JAQpQgAIUoIBjBGSim5ycjGrVqjmmAOaaK4GDBw+qz+EqV67MACRXgnxJCTAA4UCgAAUoQAEKUMBoAocOHVJVkhUQJuMI5NQvuueV/ATLOH2vtSa6B4rWyjEzClCAAhSgAAU8UiCnia5Hohig0Tn1i+55JQMQA3S6I6qge6A4oo7MkwIUoAAFKEABzxLIaaLrWRrGaW1O/aJ7XskAxDh9r7UmugeK1soxMwpQgAIUoAAFPFIgp4mu2VF69eqFX3/9VR1FbKSUU7/onlcyAMmi9+VynA8++ADR0dHYsmULTp8+jZ49e2LmzJk2j5fdu3fjpZdewh9//KHeadq0KcaOHYs6dercloc9z9pSAd0DxZYy+QwFKEABClCAAhTITiCnia5Z9GQe9t1330ECjgoVKlibxQAknYIBSBYjXSLTihUromTJkoiKisKyZcvsCkD++usvNGzYEIULF8bgwYNVKRMnToScf71p0yZUrVrVWrI9z9r6LyYDEFul+BwFKEABClCAAs4S8JQAZMaMGejTpw/WrFmDe+65x8orR97K5X+BgYHOIrepnJz6Rfe8kgFIFt0iA+TcuXPqUhY5Lk7OqrZnBaRjx45YsWIF9u3bh7Jly6pSjh8/jpo1a6J169ZYvHixtWR7nrVpFPEULFuZ+BwFKEABClCAAk4UyGmi68SqOLSorAIQhxaah8xz6hcGIHnAze2r9gYg8vlWWFgYOnfujFmzZt1UrAQxCxcuVMFNcHAw7HnWnvrrHij2lM1nKUABClCAAhSgQGYCOU10zaD2xhtv4M0337ytKV988YXa/3HrHhD5REsuAJQvZZ577jls3rxZzSOHDRuGQYMG4ciRIxgyZAh+++03+Pv7o1+/fnjrrbduu0V+/fr1ePvtt7FhwwZcv34dNWrUUPnJ3DOnlFO/6J5XcgUkpx4B7F4BkY5v0qQJPvnkE/Tv3/+mEuRnAwcOVIPjjjvuUH+19Vkbqmp9RPdAsaXsc1dvYNmOk1i55wxGPlwLNUoUsOU1PkMBClCAAhSggIcI5DTRNQPDzp07VTDx2WefYfjw4errF0ky3xs1alSmAYiXl5cKGrp37672jMgvsCUQkaBFApqHHnoItWvXVl/QyGdd8uc9evSwcn3zzTfo1KkT6tevD/myRj7xkj0oP/30E8aMGaP2JGeXcuoX3fNKBiA2jHR7V0CWLFmiOl86vm3btjeVID975JFH8NVXX6FDhw6w59msqnrq1CnI/zIm+fRLBrFsopfB6Iz0Z8x5PD59oyrq+ZbVMKTFf/tcnFE+y6AABShAAQpQwNgCOU10jV1722uX1SdYmW1Cl4Dj6NGjkCBC5oiS5EsZWRVJTEzElClTMGDAAPXzhIQE9Wm/3CS/bt069TMJXMqVK6d+sS3zzHz58lkrKvPRH3/8ESdPnkTBggWzbEBO/cIAxPa+1/akvQHI7NmzVVS6cuVK3H///TfVQyLRVq1aQZ6RAMGeZ7NqUFZLffK8MwOQlNQ0NBz9My5cS0Sd0gWwbPD/aesDZkQBClCAAhSggPsLZDfRfWbBNhw+e9UwjaxcNBgfdY7MVX3sDUAk0JAgIWOSU1MPHDiAa9euwc/Pz/pHDzzwALZt24YzZ86on1l+uS2/1L777rtvyuP777/HU089pQ5TatOmDQOQXPWmi16yNwCxZ1XDnmezar5RVkCkfi99tQOLtpxQVf1j2L0oU8hYpzy4aAixWApQgAIUoAAFAGQXgDw0aS12/3PFME55+WWqvQGInLoqn+VnTHJ9g6yMyCFGGZPsMZb5Y1JSkvqxXPEg+0WyS59//jmefPJJBiCGGV02VMTeAMSefR32PGtDVa2P6F4qs7Xsn/eeQe8vt6jHX3+oFp5qWtHWV/kcBShAAQpQgAImF+AKyO0XEVo2oVvujbMMAQlATpw4cdulhRKAyKf8Mj+V9N577+GVV15Re4+rVKmS6QiqVasWSpUq5dkBiEA9/fTTKFq0qFv8a2ZvAJLTyVYLFizA+fPnbToFK+Oz9mC5KgBJSEpB/bdWIT4xBaUL5seKZ/8PIQG+9lSdz1KAAhSgAAUoYFKBnPYamKXZsgG9d+/et90DktUeENnvkdsAxPI1zfz589UJrLlJOfWL7nmlSzahy05/+ZZNNmHLKVH/93/G3iuQXQAiy1+HDx9GaGiourTQkqRtsgdk//79ahORJMs9ILIHRAZLbp61dVDpHii2livPvbt8H6b9FqNeeaxBGYztGGHP63yWAhSgAAUoQAGTCuQ00TVLs+UXyF26dMHSpUvRrl07a7McEYDIL77Lly+v5qEbN25Uv+DOmM6ePYsiRYrcdmxvxmdy6hfd80qXBSDSaMsufVkWkqNpZVN2SEiIYcbe5MmT1c3lqampGDlyJCIjI9G+fXtVv4cffhjh4eFqWUxuTL/1kkLZNNSoUSN1jrOc3SxJjmSTlQ+5Cb169erWdtrzrK04ugeKreXKc7IK8tCkP3AoNn0j2QedItAxKj0IY6IABShAAQpQwHMFcpromkXmr7/+UidVNWzYUJ1glT9/fjRu3FjdD5LVPSC5XQERs2+//VYdw1u8eHFIkCMBSWxsrNqsLpvUZSO7j49Plrw59YvueaXLAhBL8JGWlqYw5J+DgoLQrVs31VEyuXd1shyLllk95Fxm6eCsAhB5R86BlnOXLcekybd8chZzZm2z51lbXHQPFFvKzPjMnpOX0f7j9biRnAp/Hy8sHdgEtUuF2psNn6cABShAAQpQwEQCOU10TdRUdXzuuHHj1EbylJQUdadHdhcR5iUAETc5+fTdd9/F2rVrcfHiRbXVQX7J/+ijj6pf9Gc8nvdW55z6Rfe80iUBiAQZsiR148YNK0bGQERQ5CxjwZJoLuPRY2YamI5si+6Bkpu6Lt5yHC9+tVO9Wq5wIL4bdBcKBv53jFxu8uQ7FKAABShAAQq4r0BOE133bZl71zynftE9r3RJACJdJBesyAad6dOnqyvmJWVcFbH8vXzCJOcX9+3bF5UqVXLv3nVi7XUPlNxW/ZWluzB/0zH1eqMKhfHl040Q4Oud2+z4HgUoQAEKUIACbiyQ00TXjZvm1lXPqV90zytdFoBYeklWPlasWKGODlu+fLlaosrs8yz5mVzqJ6sicpFKdstIbj0CNFVe90DJbbVuJKeg66d/IvroRZVF24hS+OjxevDy+u+WztzmzfcoQAEKUIACFHAvgZwmuu7VGvPUNqd+0T2vdHkAkrHr5JSoadOmQS5LOX36tPojCTRu/TxLTpXq16+fOspXNtsw3S6ge6DkxVhuRu/wyXocOXdNZTPgnsoY9kCNvGTJdylAAQpQgAIUcEOBnCa6btgkU1Q5p37RPa80VABi6UE59lb2iMiqyG+//WYNRORvMgYjvr6+6mgzOcq3WbNmphgAuhqhe6DktV5/n7uG9p+shwQjkka3q4NujcvnNVu+TwEKUIACFKCAGwnkNNF1o6aYqqo59YvueaUhA5CMPSr3aHz88ceYPXs2Ll++nOWm9Ro1aqjTs3r06IECBQqYalDkpjG6B0pu6nDrO/IZVtdPN6qTsfLlAyZ1icRD4VnfyqmjTOZBAQpQgAIUoIBxBHKa6Bqnpp5Vk5z6Rfe80vABiKX7r1+/jrlz52Lq1KkQBEkZP8+y/HNgYCC6du2qgpF69ep51ujJ0FrdA0UX5IrdpzFwbjRS0wAfr3z4tEcD3FujmK7smQ8FKEABClCAAgYWyGmia+Cqm7pqOfWL7nml2wQgGXt98+bNmDRpEubMmXPbioglEJG/3n333eoejtatW5t60GTWON0DRSdgxuN55Y6QWU81wh2VwnQWwbwoQAEKUIACFDCgQE4TXQNW2SOqlFO/6J5Xul0AIpvTZ8yYgU8//RSyaT2r07Bkr4jlz+QCQHlebqT0lKR7oOh2+/yPIxi1bK/KNtjfB/P6NEZ4mYK6i2F+FKAABShAAQoYSCAmJgZJSUlqTsYTTY3RMampqZCb22VvdVZXXuieV7pNALJmzRq1KV2umpdN6hmTZWO63BkiNz8K5K1H+ebPnx+LFi1SR/h6QtI9UBxh9tHPf2H8zwdV1gUDfbGw752oXiLEEUUxTwpQgAIUoAAFDCBw8uRJtae3YMGCKFGiBIMQF/eJBIOxsbG4cuUKZB5drFjmn8XrnlcaOgCRATpz5ky17+PgwfSJasZTsCx/LyscgwYNQvv27XHq1Cm12iGrJLJakjEQkf0hu3fvRoUKFVzc3Y4vXvdAcUSNpf9G/7APM/5Iv4iySLAf5ve5A1WLMwhxhDfzpAAFKEABCrhaQO57O3bsGBISEuDl5QVvb15O7Ko+kV/YS39ICggIQNmyZeHj45NpdXTPKw0ZgMgeDwk6Fi5cCNl8fus9IPLPls3mEniEh4ffhiURnQQhr732mloVsXyS9cwzz2DcuHGu6munlat7oDiq4tIvw7+W29KPW4OQeX3uQDUGIY4iZ74UoAAFKEABlwrIpPfMmTMqCJFJMJNrBCQAlIBDTo8NDQ3NdjVK97zSMAGIBBrz5s1Tn1lt27ZN9URmqx2VK1dWJ1w99dRTavkupyRRdv369a1BSK1atdQqiNmT7oHiSK/U1DS8snQXFm5hEOJIZ+ZNAQpQgAIUoAAFciOge17p8gBk3759arXjyy+/VN+fZbbaIZ9RtWrVSn1mJSda2btpaeTIkXjrrbeUt6ycXL16NTf2bvWO7oHi6MZLECIrIQs2pwchYUF+mN+XKyGOdmf+FKAABShAAQpQICcB3fNKlwQgGW86//3331WbMws8ZIWjV69e+N///gdZ+chtklvVO3bsqF6X4MXyvVtu83OH93QPFGe0WYKQEd/89zkWgxBnqLMMClCAAhSgAAUokL2A7nmlSwKQkiVLqh33lsDj1hOr6tatq4KO7t27qxWLvKbVq1ejZcuW1pUTBiB5FXXc+5kFIbInhKdjOc6cOVOAAhSgAAUoQIHsBEwRgMiml4yfUcnqh2yCefTRR9VnVnKBoM5kCUC4AqJT1XF5pQchuzF/0zFViKyEMAhxnDdzpgAFKEABClCAAh4TgEjgUbx4cfTp0wf9+/dHqVKlHNL7cgqWRG6W1KJFC4eUY6RMdUeqzm4bgxBni7M8ClCAAhSgAAUokLmA7nmlSz7BkhWQO+64Q612dOrUSd28yKRXQPdA0Vs723KTIOTVb3dj3p/pKyGFZWM6P8eyDY9PUYACFKAABShAAU0CuueVLglA5JjdyMhITSTMJjMB3QPFVcoShLz27W7MzRCEzHm6MWqVKuCqKrFcClCAAhSgAAUo4FECuueVLglAPKrHXNRY3QPFRc1Qxd4ahITm98WspxqhXtmc74FxZb1ZNgUoQAEKUIACFDCDgO55pUsCELlw0JLuu+8+FCtWLNd9IzdpyiZzS+ratWuu8zLTi7oHiqttJAh54/s9+HLDUVWVID9vfN6rIRpXCnN11Vg+BShAAQpQgAIUMLWA7nmlSwKQjKdgrVq1Cs2bN891p3niEbu2YOkeKLaU6ehn5NCC91bsx7TfYlRRAb5emPZEAzSrVtTRRTN/ClCAAhSgAAUo4LECuueVLgtApAflKF5dAYglP0+448OW0a97oNhSpjOekSBk0i+HMG7VQVWcn7cXJnWNRKvaJZxRPMugAAUoQAEKUIACHiege17psgDEcg8IAxDHjGHdA8Uxtcx9rjPWxuDtH/apDLy98mHcYxF4pF7p3GfINylAAQpQgAIUoAAFMhXQPa9kAGLSgaZ7oBiRae6fR/HqN7uRliaracC77eqic6NyRqwq60QBClCAAhSgAAXcVkD3vNLtA5AffvgBbdu2VR0aGBiIq1evum3n6qy47oGis24681q69QSGLt6B1LT0XF9/qBaealpRZxHMiwIUoAAFKEABCni0gO55pdsHIJMnT8aQIUPUoJBb1U+dOuXRA8TSeN0Dxcioy3edwpAF25CUkh6FDL2/Gv53bxW1x4iJAhSgAAUoQAEKUCBvArrnlW4fgDRr1gxr165Vk80mTZqov2cCdA8Uo5uu2R+LfnOikZicqqr6dNOKGPFgTXh5MQgxet+xfhSgAAUoQAEKGFtA97zSYQHI+vXrs5Rs2rSp9bfT48ePR4MGDWxWl1OQ4uPjERMTg4ULF+K3336D/EwCkMGDB2PChAk252XmB3UPFHewWn/4HPrM2oJriSmquu0jS2NMx3D4enu5Q/VZRwpQgAIUoAAFKGBIAd3zSocFIBnv+rhVUgIGS8rrZzKWvCQfWf2QVRAmz1sBsfT5zhOX0OuLzbhwLVH9qHmNYpjStT7y+3lzWFCAAhSgAAUoQAEK5ELA7QKQjMFGLtqb7SuW4EXK6NChAxYvXqy7CLfNT/dAcSeIw2evosdnm/DPpeuq2g0rFMKMng0Rmt/XnZrBulKAAhSgAAUoQAFDCOieVzp8BcSRAYj0iJ+fH3r37o0PP/wQ/v7+hugkI1RC90AxQpvsqcOpy9fxxGebcCg2/VS0GiVC8OXTjVAsJMCebPgsBShAAQpQgAIU8HgB3fNKhwUg3bt3z/IUorlz51r/rHnz5ihRwvZbrOXTrqCgIISFhSE8PBz33nsvihQp4vED41YA3QPFHYEvXkvEkzM3Y/vxS6r65QoHYvbTjVA+LMgdm8M6U4ACFKAABShAAZcI6J5XOiwAyU4n4/6QvN6E7pJecINCdQ8UN2hyplW8diMZ/edEY+1f59SfFwn2x6ynGqJ2qVB3bRLrTQEKUIACFKAABZwqoHte6bIARNRkDwcDEMeMH90DxTG1dE6ucjTv84u2Y9nO9Dtigv19ML1HFJpU5sqZc3qApVCAAhSgAAUo4M4CuueVLglAXnvtNWsfPPXUU6hYkTdX6x6UugeK7vo5O7+U1DS8+f0efLnhqCraz9sL4x+vhzbhJZ1dFZZHAQpQgAIUoAAF3EpA97zSJQGIW4m7aWV1DxQ3Zbip2nIgwpQ1h/DBTwfVz+Wi9Dcfro0ed1YwQ/PYBgpQgAIUoAAFKOAQAd3zSgYgDukm12eqe6C4vkX6arBw8zG8snQXUv+9jmbgPZXxYqvqWR6aoK9k5kQBClCAAhSgAAXcT0D3vJIBiPuNAZtqrHug2FSoGz30894zGDR/KxKSUlWt29cvjTEdeGu6G3Uhq0oBClCAAhSggJMEdM8rGYA4qeOcXYzugeLs+jujvOijF9F71mZcjE9Sxf1f1SL4pHuU2qTORAEKUIACFKAABSiQLqB7Xqk1AFm/fv1t/dSkSZPbfpbZc7o6OLPydOXtTvnoHiju1HZ76iq3pvf8fBNOXEy/Nb1O6QL4vFdDXlhoDyKfpQAFKEABClDA1AK655VaA5CM93tIL8gxu8nJybd1yK3P6eqxrMrTlb875aN7oLhT2+2ta+yVBPT6YjP2nrqiXi1bOD9mPdkIlYoG25sVn6cABShAAQpQgAKmE9A9r3RIACKnDVkCkJSUlCwDEMtzunpJApDMytOVvzvlo3uguFPbc1PXuIQkDJizFX8cSr+wsHCQHz7r2QCR5QrlJju+QwEKUIACFKAABUwjoHteyQDENEPj5oboHigmZbqpWXJh4bAlO/H1tn/UzwN8vTCla320qFncE5rPNlKAAhSgAAUoQIFMBXTPK7UGIN27d7/tKNPZs2ff1pDMntPV35mVpytvd8pH90Bxp7bnpa6pqWkYu/IApv52WGXjlQ94p11ddG5ULi/Z8l0KUIACFKAABSjgtgK655VaAxC3VTVhxXUPFBMSZdukmeuO4M1le/Hv14R49r6qeKZFVd4V4mkDge2lAAUoQAEKUMDYp2Cxf4wjwAAk733x465TeHbBdiSmpN8V0rlhWbz9aB34eHvlPXPmQAEKUIACFKAABdxEQPe8kisgbtLx9lZT90Cxt3yzPL8x5jz6fLkFcQnpp7m1qFEMk7pGItCPd4WYpY/ZDgpQgAIUoAAFshfQPa9kAGLSEad7oJiUyaZmHTgdh15fbMKpywnq+YiyBdUJWUWC/W16nw9RgAIUoAAFKEABdxbQPa9kAOLOoyGbuuseKCZlsrlZJy9dV0HIwTNX1TvlCgdi5pMNeVeIzYJ8kAIUoAAFKEABdxXQPa90uwAkNTUVFy9eRMGCBeHt7e2u/ejweuseKA6vsBsUcPl6EvrN3oKNMRdUbQsG+mJGjwZoUKGwG9SeVaQABShAAQpQgAK5E9A9r3SbAGTlypUYM2YM1q1bp25Xl0sH69Wrh4EDB+Kpp57KnaaJ39I9UExMZVfTbiSnYNhXO/HN9pPqPT8fL4x/rB7ahJe0Kx8+TAEKUIACFKAABdxFQPe80iUByLfffotXXnnFar5kyRLUrFkzyz547733MGLECPXnt96eLoHIQw89BMnDx4cbgy2IugeKu/wL4ox6yhj88KeDmLzmkLW4EQ/WRO//q8hjep3RASyDAhSgAAUoQAGnCuieV7okAOnYsSOWLl2q4MLDw7F9+/YsEVevXo37779fBR4SbNyaLD/v06cPpk6d6tTOMHJhugeKkdvqqrot2HQMI77ZjZTUNFWFnneWx+tta8Nbbi9kogAFKEABClCAAiYR0D2vdEkAUqRIEbWPQ9KoUaOsqxuZ9VGDBg3U5ScSfFhWP8LCwnD16lXcuHHD+nMvLy9s3LgR8jwTtF8YQ9PMBX49EIv/zd2Ka4kp6oH7ahbHxC71eEwvBwwFKEABClCAAqYRcPsA5PDhw6hatarqEAkq/vzzzyyDBvmzO++807ryUatWLSxYsAC1a9dWwcf777+P119/3frnXbp0wZw5c0zT2XlpiO6Bkpe6mP3dPScv48kvNiM27oZqakSZUMzo2RBFQ3hMr9n7nu2jAAUoQAEKeIKA7nml01dAZDN569atVV/JKVbXrl2Dn59fpn33wgsvYPz48dZnd+/ejerVq9/0bM+ePTF79mz1s+DgYJw/fx6+vr6eMBaybaPugeLxoDkA/HPpOp7McExvmUL5MfPJRqhSLJh0FKAABShAAQpQwK0FdM8rnR6AzJgxA3379lWrFuXLl0dMTEyWHSKnXO3cuVM9K/tAli9fftuzO3bsQGRkpPq5PLdlyxbrP7t1T+ex8roHSh6r4xGvyzG9A+ZEY/3h86q9ofl9Mf2JKDSuFOYR7WcjKUABClCAAhQwp4DueaXTA5Bx48Zh6NChKliIiIhQexUyS5cuXYLs9bCkKVOmoH///pk+K3eCXLlyReX55Zdfolu3bubsfTtapXug2FG0Rz+amJyKl5fuxNKt/ygHX+98GNMhHO3rl/FoFzaeAhSgAAUoQAH3FdA9r3R6ADJ69Gi89tprKli444471L0emaWMn2rJs/v377fuHbn1+YwrJfLJ1pAhQ9y3hzXVXPdA0VQtj8hGDksY//NfmLj6L2t7hzSvgudaVuMxvR4xAthIClCAAhSggLkEdM8rnR6AfPDBB3jppZfUREyO4N22bVumPSSby99++231Z3JqVmxsbJY9edddd2HDhg0qz7feegvDhw83V6/nojW6B0ouquDxryyJPqFWQ82Xo48AACAASURBVJJS0o/pbRtRCu93DEeAr7fH2xCAAhSgAAUoQAH3EdA9r3R6ADJ9+nTrp1SFChVSm8YzS82aNcPatWtVUCEXDcrlhVmljCsgcmnhiy++6D496qCa6h4oDqqm6bP9M+Y8+s2JxqX4JNXW+uUK4tMeDRAWzBOyTN/5bCAFKEABClDAJAK655VOD0DkYsGWLVuq7sjq06oLFy6gRIkSSElJv1shp6CiXLlyOHHihMovu70iJhkDNjVD90CxqVA+lKnAkXPX8NTMzZC/SipbOD8+79kQVYuHUIwCFKAABShAAQoYXkD3vNLpAYiseBQrVswK/fzzz6v7PDImy2daliBFPq9q1KhRpp0TFxcH2YRuSbJSIismnp50DxRP98xr+y/FJ6Lf7Gj8eeSCyiokwAefdItC06pF8po136cABShAAQpQgAIOFdA9r3R6ACI68nnVH3/8oW429/HxwbRp09CrVy+1gvH777+jbdu26qZz+fOyZcvi6NGjWaLKZ1qSnyVY2bNnD2rUqOHQTnCHzHUPFHdos9HrKCdkvbJ0F5ZsPaGq6u2VD28/WgddGpUzetVZPwpQgAIUoAAFPFhA97zSJQHIokWL0LlzZxVwSJAhfy1QoACCgoJw6tQp9TNLQCEb0V955ZUsu1xO1JKTtSTJRYSXL1/mSUOAOt44KioK0dHRqF+/vgf/K2OspsvY/vjXw3h/5QFrxfrdXQnDHqgBL698xqosa0MBClCAAhSgAAUcMK90SQAiPSm3octRu5YgJGPvWn5WsWJFdRGhBCZZpbp160JWPeSde+65B7LHhIkBiNHHwA87T+H5RdtxIzlVVfX+WsUxoXM9BPr5GL3qrB8FKEABClCAAh4moPsX2y4LQK5du4Z27drh559/Vl0oAYQky+pHqVKlVIBSu3btLLt406ZN6i4Ry7uyEvLyyy972JDIvLm6BwpR9QtsPXYRfb/cgnNXE1XmdUuHYkbPBiheIEB/YcyRAhSgAAUoQAEK5FJA97zSZQGIpf1LlizBwoUL1UWD8fHxKF26NB544AEMHDgQoaGh2TJ17doVCxYssD6zY8cOyIoIE1dA3GUMHL8Qj6dnbcbBM1dVlUuGBuCzng1Rq1QBd2kC60kBClCAAhSggMkFTBeA5KW/Ll26ZD2qV/IJCwvLS3amelf3QDEVjsEacyUhCf+buxVr/zqnahbk541JXSPRvEZxg9WU1aEABShAAQpQwBMFdM8rXb4CYuROlHtI5IjgGTNm4Pjx4+pErt69e6uLDr29c77NeubMmZg4cSL27t2LwMBA3HvvvepOk6pVq97UbDkBbNasWZlSSLllypSxm0n3QLG7AnzBLoHklFSM/G4P5v55TL0n+9Fff6gWet1V0a58+DAFKEABClCAAhTQLaB7XskAJJseks/APvnkEzz55JNo0qQJ1q9fjy+++EJ9HiYXHmaXRo4ciVGjRqF58+Zqr4uczjV58mQkJydj8+bNqFChgvV1SwAiQYiXl9dN2bZv314FL/Ym3QPF3vL5vP0Csv/p83V/4+0f9uLfg+DQ887yeO2hWvDxvnlc2J8736AABShAAQpQgAK5E9A9r2QAkkU/7Nq1CxERERg8eDA++ugj61PPPPMMJk2ahOz2m5w+fRpyO7useKxYscK6SV7uM6lVqxYeeeQRzJs377YAJCkpSd2LoiPpHig66sQ8bBNYtfcMhszfhutJKeqFe6oXxaQukQgJ8LUtAz5FAQpQgAIUoAAFNAronlcyAMmic0aMGIF33nkHMTExkOOALenIkSOoVKkShg8fbr1/5NYsli5dig4dOmD27Nno3r37TX/cpk0b/PLLLzh37pz1eGHLCkhiYiKuX7+u7jO5dSXE3jGke6DYWz6fz5vA7n8uq83pZ67cUBnVKBGCz3o1ROmC+fOWMd+mAAUoQAEKUIACdgronlcyAMmiA1q1aqVWOWQ149ZUvHhxREZGqtWNzJKczNWlSxfICV/yCVXG1KlTJ3z11VfYuHEjGjdurP7IEoCEhIQgLi4O+fPnV/ekjB07FpUrV7ZziKQ/rnug5KoSfClPAqcvJ6ggZM/JKyqfIsH++KxnA0SULZinfPkyBShAAQpQgAIUsEdA97zSEAHI4cOH8cMPP6hbu+UzJdkvIUfyWu4EsRVI7gM5cOC/G6ZtfS+z5+Q4Xz8/P1WnW5PcLC6fS8lnWpklCVzq1aunPt+STeiWlJCQoDagnzhx4qbgRO4ukfzk5nJ/f38VnMhnXhKQbNmyBeXLl8+2KXJ7vPwvY9q3b59afeFN6HkZBa5/99qNZDyzYDt+3ndGVSbA1wsTHq+HB+qUdH3lWAMKUIACFKAABTxCwFQByF9//aUm6XIZYcZgw97Aw9LzEoDIyVU6kqw8yEqHbDy/NcmG9NjYWBw6dCjLopo1a4YNGzaoz7gsm9BlY7pcrijBRmafZ2XMbNWqVbj//vvRs2dPyGla2aU33ngDb775ZqaPMADRMRpcm0dKahre/XEfZvxxxFqRl1vXQL+7K1n3F7m2hiydAhSgAAUoQAEzC5gmAFm+fDnkcyTZ82AJOCw3muemAyUPnQFIXlZApP4SoEjwkPEzrTvvvFOdiiU3tn/zzTdqM3p2SVZEZGXj5MmT2T7HFZDcjBj3e2f2xqN447s9kIBEUqeoMni7XR34++R8JLT7tZY1pgAFKEABClDAKAKmCEDk1vMGDRqoz6wkSeCQ21WPjB2jMwDJyx6QjHWSz61kI3uxYsVQo0YNyOdWY8aMgXwiJf+cXZL9I8uWLYNsTrc36R4o9pbP5x0j8NvBsxg0dyvibiSrAqLKF8LU7lEoGuLvmAKZKwUoQAEKUIACHi+ge17pkj0gjz76KL777jvr5yMSfMjpUN26dVP7IOTTp6CgoFx1ti0XBNqSsZxy9e677+bqFKzs8pe9IbLHRYKSnFZ8wsPDcfHiRXUJor1J90Cxt3w+7ziBg2fi0HvWFhy7kB7AlwoNwPQeDVCndKjjCmXOFKAABShAAQp4rIDueaXTAxCZUBctWlSteMj/AgICsHDhQrRt29ZQnSobyeWkq6zuAdm+fTskQJD9HLKJPjQ0FCVLZr8xWG5U79Onj7rcsH///qq9165dU3d/yObzjElMOnfujH79+mHq1Kl22+geKHZXgC84VODitUQMnLsVG2LOq3Ly+3rjg04RaBPOzekOhWfmFKAABShAAQ8U0D2vdHoA8u2336pN2Zbf/n/44Yd49tlnDdmVEiRMmzZN3YR+1113Yd26deom9IxBwd9//63uCbl1s/jQoUNx9uxZtaIjQZbc/bFo0SI8/vjj6hJCS/slkJHVH1kVkhOy5OQt2bwuz5QtWxZ//vmnWhGyN+keKPaWz+cdL5CUkopR3++F7A2xpCEtquLZFlXh5ZXP8RVgCRSgAAUoQAEKeISA7nml0wOQjz/+GIMGDVKdJZPt8+fP5/pzK0f3eHJysrqLQ1YuZC9HmTJl0Lt3b7z00kvWG8uzCkAkgBg3bhwOHjyoTuaS/R59+/ZVKyAZLxmUe0ZeeOEFbN68WW02lxUVCTxkRUguQyxSpEiumql7oOSqEnzJKQJz/t2cnvzv5vRWtYtj3GP1EOTv45TyWQgFKEABClCAAuYW0D2vdHoAIsfSvvrqq2oFoFatWlnepWHubnR863QPFMfXmCXkRWBjzHkMmBONi/FJKpvqxUPwaY8GKBcWmJds+S4FKEABClCAAhTQfsG10wMQuZhPPrmSAKRhw4bq0j0m/QIMQPSbGj3H4xfi0efLLdh/Ok5VtWCgLz7uWh9NquRuFc3o7WX9KEABClCAAhRwjoDueaXTAxC5F+PBBx9UAYjsncjuMj/nkJqzFN0DxZxK5muV3Jz+wqIdWLHntGqct1c+vNamJno2qZDjqWvm02CLKEABClCAAhTQIaB7Xun0AOTq1avqFKwbN26ovRCyByK3+xx0gJo1D90DxaxOZmxXamoaJq85hHGrDlqb91iDMnjrUV5aaMb+ZpsoQAEKUIACjhbQPa90egAiQHKq1KxZs9RvZGWj9jPPPONoN4/LX/dA8ThAEzR45Z7TeH7hdlxLTFGtqV+uIKY+EYViIQEmaB2bQAEKUIACFKCAswR0zytdEoDIaU+yAT0uLg5hYWFqY4ucMMWkT0D3QNFXM+bkTIEDp+PUvhDLpYUlCgRg2hNRiChb0JnVYFkUoAAFKEABCrixgO55pUsCEPH/8ccf1X0gcuysBCNyP0jlypXduGuMVXXdA8VYrWNt7BGQSwsHzd+KdYfSLy308/HCmA510S6SQb89jnyWAhSgAAUo4KkCuueVLgtApANXrlyJrl27Qm5HDw4OVhf8devWDREREdwwm8cRrnug5LE6fN3FAskpqXjnx/34fN0Ra0363l0Jwx6ooTaqM1GAAhSgAAUoQIGsBHTPK10SgFSrVs3avgsXLkD+J8lyO7i/v7/6NCt//vx2jQR5/8CBA3a9Y9aHdQ8Uszp5WrsWbTmOV7/ejcSUVNX0ZtWKYmLnSIQG+noaBdtLAQpQgAIUoICNArrnlS4JQOT0KwkW0tLSblrpkH/OS5I85dZxJmi/MIam5hHYeuwi+s2Oxtm4G6pRFYsEqUsLqxQLNk8j2RIKUIACFKAABbQJmCoA0aYCWIMZBiDpqroHis6+Yl6uFzh9OQH9Zm/BjhOXVWVC/H0woXM9tKhZ3PWVYw0oQAEKUIACFDCUgO55pctWQByhyhWQ/1R1DxRH9BfzdK1AQlIKhi/dhaXb/lEVyZcPGHp/dQy8pzL3YLm2a1g6BShAAQpQwFACuueVLglAHLlK4e3tbagOc1VldA8UV7WD5TpWQD57/OyPI3jnx31I/fcLyIfCS+L9jhHI78d/lxyrz9wpQAEKUIAC7iGge17pkgDEPajdu5a6B4p7a7D2OQn8fvAsBs3biisJyerR2qUKYHqPBihd0L6DIHIqh39OAQpQgAIUoID7CeieVzIAcb8xYFONdQ8UmwrlQ24tcOTcNXVp4aHYq6odYUF++KR7FBpVLOzW7WLlKUABClCAAhTIm4DueSUDkLz1h2Hf1j1QDNtQVkyrQFxCEp5buB0/74tV+fp45cObj9RGt8bltZbDzChAAQpQgAIUcB8B3fNKBiDu0/d21VT3QLGrcD7s1gKpqWkYt+ogJq85ZG1H9zvKYWTb2vD19nLrtrHyFKAABShAAQrYL6B7XskAxP4+cIs3dA8Ut2g0K6lVYNnOk3hx8U5cT0q/W0c+xfqkW32EBftrLYeZUYACFKAABShgbAHd80pDBSB79uzBb7/9hr179+L8+fO4dOmSuljw/fffR0REhLF7xmC10z1QDNY8VsdJAntOXkbfL6Pxz6XrqkTZlD69RxRqlwp1Ug1YDAUoQAEKUIACrhbQPa80RAAye/ZsjB07VgUeGZPlpvRVq1ahefPmt9kPGjQIv//+u/p5xYoV8e2337q6fwxTvu6BYpiGsSJOFzh/9QYGzN2KTUcuqLLz+3rjg04RaBNe0ul1YYEUoAAFKEABCjhfQPe80qUByLlz59C1a1esXr1a3WQuSS4TlJTxn7MKQH744Qe0bdvW+t62bdsQHh7u/F4xYIm6B4oBm8gqOVEgMTkVo5btwZyNx6ylDrq3Cp5vWQ1eXun/zjJRgAIUoAAFKGBOAd3zSpcFIBcuXMCdd96JQ4cOqWBDAg9L0GEJRHJaAZHnZOXj6NGj6v2XX34Zo0ePNmfP29kq3QPFzuL5uEkF5v55FCO/3YPkf28tvK9mMYx/vB5CAnxN2mI2iwIUoAAFKEAB3fNKlwUg8knVr7/+ag088ufPj169eqF169YoV64c6tWrp3pbAousVkDkz59//nlMmDBBPdegQQP8+eefHCUAdA8UolLAIiCfYg2YE43z1xLVj6oUC8anPRqgYpEgIlGAAhSgAAUoYEIB3fNKlwQgS5cuRceOHa2fW9WtW1ft3yhf/r+7Bry8vKx/nl0AIn/WqlUr1dXe3t5q43pQECdCugeKCf9dYpPyICCb0vt+uQV7Tl5RuRQI8MHkrvVxd7WieciVr1KAAhSgAAUoYEQB3fNKlwQgjRo1wpYtW5RvqVKl1G/rixUrdpO3rQGInJZVtGj6pEdWQTZt2oSoqCgj9p1T66R7oDi18izMLQSuJ6bgxa92YNnOU6q+shVk+IM18XTTitZfHrhFQ1hJClCAAhSgAAWyFdA9r3R6AHL27FmUKFHC2shp06ahd+/etzXa1gBEXixZsiTOnDmjJj1z5sxBly5dPH4Y6R4oHg9KgEwFZJ/Wx78exgc/HcC/50igff3SeKddXQT4elONAhSgAAUoQAETCOieVzo9AJFPrdq1a6e6wsfHB7IZPTg4OE8BSJ06ddQRvhKAfPTRR5DjeT096R4onu7J9mcvsHrfGTyzYDuu3khWD0aULYhp3aNQIjSAdBSgAAUoQAEKuLmA7nml0wMQWfEYMGCAChaqVauGffv2Zdol9qyAyGlasvlc8pRTsOQ0LE9PugeKp3uy/TkLHIqNQ58vo3Hk3DX1cLEQf0x9Igr1yxXK+WU+QQEKUIACFKCAYQV0zyudHoC8++67GDFihAoWZC/Ihg0b8hyAyC3pu3btUnnKhYYvvPCCYTvQWRXTPVCcVW+W494Cl+OTMHjBNvx+8KxqiJ+3F0a3q4NODcq6d8NYewpQgAIUoIAHC+ieVzo9AJkyZQoGDx6sgoVatWqpwCGzZM8KSMY9IJ999pk6ztfTk+6B4umebL/tAimpaRi7Yj+m/R5jfenJuypgxIM14ePtZXtGfJICFKAABShAAUMI6J5XOj0AWbRoETp37qwwQ0JC1LG5ltvPMwrbGoAcO3YMFSpUUK9KPsuWLVN3iXh60j1QPN2T7bdf4OttJzBsyS7ILeqS7qoShsld6qNQkJ/9mfENClCAAhSgAAVcJqB7Xun0AGTnzp03XTK4ceNGNGzY8DZQWwMQuYRQLiOUJO/IaVhhYWEu6yCjFKx7oBilXayHewnsPHEJfb+MxukrCari5QoHYkbPBqhWPMS9GsLaUoACFKAABTxYQPe80ukBiPSdfDIVGxurulE+l5LPpm5NtgQg8fHx6jOu48ePQ44DlQsNd+zY4cHD47+m6x4oRKVAbgVi4xLQf3Y0th67pLII8vPGuMfroVXt/47jzm3efI8CFKAABShAAccL6J5XuiQAkVOw5DQsSRJo/PLLL7j77rtv0rMlAHniiScwd+5c9Z58fvXee+/hxRdfdHwvuEEJugeKGzSZVTSwwI3kFLz2zW4s2nLCWsvn7quGwc2rwEtuMGSiAAUoQAEKUMCwArrnlS4JQGTFQo7gTUxMVCsXoaGh+Prrr3HPPfdY4bMLQC5evIg+ffqodyRJHoUKFYLsBwkKCjJs5zmzYroHijPrzrLMKSD/nn654ShGLdsL2aguqXWdEvigUwSC/H3M2Wi2igIUoAAFKGACAd3zSpcEINIPb7zxBkaNGqVWLmRiIn997LHH8OSTT6rjeQsXLqy6S37+888/46677kJ0dLQKOj799FNcuXJFvWd5ZsaMGepdpnQB3QOFrhTQJbD+0DkMnLcVl+KTVJY1SoTg0x4NULZwoK4imA8FKEABClCAAhoFdM8rXRaAiEnHjh2xdOnSm4KQjFaWwERWN+S0LEvAkTHwkL/v3bs3pk+frpHZ/bPSPVDcX4QtMJLAsfPx6PPlFhw4E6eqVSjQF1O61UeTykWMVE3WhQIUoAAFKEABB/xi26UBiHyCNWjQIMjqheUoXktwkV1vW1ZN5Bm5U2T8+PFqLwnTfwIMQDgajC5w7UYynl+0HSv3nFFV9fbKh5Fta+GJO8pnejS30dvD+lGAAhSgAAXMKqB7XunSAMTSSQsWLMDrr7+OQ4cOqR9ldi+I5VlLgCJ3f4wePRpdunQxa1/nqV26B0qeKsOXKZCFQGpqGib+8hcm/PyX9YnODcti1CN14OfDXypw4FCAAhSgAAWMIKB7XmmIAERgJbCQ/R3ff/89fv31Vxw9evQ2b9kX0qJFCzz00EMq8PDx4cbVrAal7oFihMHPOphXYMXu02o1JD4xRTWyQflC+KR7FIqG+Ju30WwZBShAAQpQwE0EdM8rDROA3Op//fp1nD9/HnLiVWBgIIoWLYoCBQq4STe5vpq6B4rrW8QamF1g/+kral/I8QvXVVNLhgZg+hMNULdMqNmbzvZRgAIUoAAFDC2ge15p2ADE0L3gBpXTPVDcoMmsogkELl5LxP/mbcX6w+dVa/x9vDC2YzgeqVfaBK1jEyhAAQpQgALuKaB7XskAxD3HQY611j1QciyQD1BAk0BySire/mEfZq7/25pjv7sr4cVW1eHjzX0hmpiZDQUoQAEKUMBmAd3zSpcEIPPmzbM2+L777kOxYsVsBrj1wTNnzmD16tXWH3ft2jXXeZnpRd0DxUw2bIt7CCzafByvfrMbiSmpqsJNKodhUpdIhAVzX4h79CBrSQEKUIACZhHQPa90SQCS3S3n9naUBB8tW7a0npyVkpK+idXTk+6B4umebL9rBKKPXsSAOdGIjbuhKiD7QmRzer2yBV1TIZZKAQpQgAIU8EAB3fNKlwUg0ndy3O6qVavQvHnzXHelJQCx5McAJJ1S90DJdQfxRQrkUSA2LgGD5m7Dpr8vqJz8vL3w5iO10aVRuTzmzNcpQAEKUIACFLBFQPe80mUBiOWuDwYgtnS7/c/oHij214BvUECfQFJKKt79cT8+X3fEmunjDcqqQCTA11tfQcyJAhSgAAUoQIHbBHTPKxmAmHSQ6R4oJmVis9xM4Nvt/+DlJbtwPSn9U8vwMqH4uFt9lCkU6GYtYXUpQAEKUIAC7iOge17p9gHIDz/8gLZt26oelPtCrl696j696cCa6h4oDqwqs6aAXQJyX0i/2dE4ej5evVco0BeTutRH06pF7MqHD1OAAhSgAAUoYJuA7nml2wcgkydPxpAhQ5Re8eLFcerUKdskTf6U7oFici42z80ELl9PwvMLt2P1/lhVc698wNBW1TGgWWXrgRRu1iRWlwIUoAAFKGBYAd3zSrcPQJo1a4a1a9eqSUeTJk3U3zNxEzrHgPkFUlPTMHnNIYz/+SDS0tLb26p2cXzQKQIhAb7mB2ALKUABClCAAk4ScJsAZP369VmSNG3a1PpbyvHjx6NBgwY286WlpSE+Ph4xMTFYuHAhfvvtN8jPJAAZPHgwJkyYYHNeZn5Q90AxsxXb5t4Caw7E4tkF2yGrIpIqFQ3CtO5RqFo8xL0bxtpTgAIUoAAFDCKge17psBWQjHd93GonAYMlWU7Dyq2vJS/JR1Y/ZBWEiSsgHAOeJXDsfDz6z4nG3lNXVMOD/LzxfqcIPFi3pGdBsLUUoAAFKEABBwi4XQCSMdjQ7WEJXqSMDh06YPHixbqLcNv8dA8Ut4VgxT1G4HpiCkZ8vQtLt/1jbXO/uyvhxVbV4ePt5TEObCgFKEABClBAt4DueaXDV0AcGYAIrp+fH3r37o0PP/wQ/v7+ur3dNj/dA8VtIVhxjxKQ/97M3ngUo77fi+TU9JXWOyuFYVLXSBQJ5n8fPGowsLEUoAAFKKBNQPe80mEBSPfu3bM8jWbu3LnWP5Nb0EuUKGEzkHzaFRQUhLCwMISHh+Pee+9FkSI8fvNWQN0DxeYO4oMUMIBA9NELGDBnK2LjbqjalAwNUPeFRJYrZIDasQoUoAAFKEAB9xLQPa90WACSHWvG/SF5vQndvbrPebXVPVCcV3OWRAE9ArFxCRg0dxs2/X1BZejn7YXX29ZCt8bleFSvHmLmQgEKUIACHiKge17psgBE+kv2cDAAcczI1T1QHFNL5koBxwokpaTi3R/34/N1R6wFdahfBqPb1UGAr7djC2fuFKAABShAAZMI6J5XuiQAee2116zd8dRTT6FixYom6R7jNEP3QDFOy1gTCtgv8N2Okxj21U5cT0pRL9cqWQBTu0ehXFig/ZnxDQpQgAIUoICHCeieV7okAPGwPnNJc3UPFJc0goVSQKPAwTNx6D87GjHnrqlcQ/P7YsLj9XBvjWIaS2FWFKAABShAAfMJ6J5XMgAx3xhRLdI9UEzKxGZ5mEBcQhKGLt6BlXvOqJbnywcMaV4Vz7SoCi+vfB6mweZSgAIUoAAFbBPQPa9kAGKbu9s9pXuguB0AK0yBLATkqN5pv8dg7Ir9+PekXjSvUQzjH6+nVkWYKEABClCAAhS4WUD3vJIBiElHmO6BYlImNsuDBdYfOodB87fhwrVEpVA+LBDTnohCjRIFPFiFTacABShAAQrcLqB7XmmYAGTnzp2Ijo7G0aNHcfnyZcTHx8PeSwzlVK1p06Zx3PATLI4BCtgk8M+l6xgwJxo7T1xWz+f39cY77eugXWQZm97nQxSgAAUoQAFPEDBVAJKUlKRuMJeg4dixY3nqPwlWJABJSUk/5cbTk+6B4umebL95BRKSUvD6t7uxaMsJayPlrhC5M8Tfh0f1mrfn2TIKUIACFLBVQPe80mUrIEeOHEHbtm2xb9++m1Y6JIjITWIAcrOa7oGSmz7hOxRwFwH578eCzccx8rs9SExOVdUOLxOKKV3ro2xhHtXrLv3IelKAAhSggGMEdM8rXRKAXLx4EY0aNcLhw4eVkiXosPeTq1uJuQLyn4jugeKY4cxcKWAsgV0nLmPgvGgcv3BdVUw2pY9/PALNaxQ3VkVZGwpQgAIUoIATBXTPK10SgDz//POYMGHCTYFHkSJF0L59e0RFRaF48eIICgrKFWuLFi1y9Z7ZXtI9UMzmw/ZQICuBy/FJeGHxdvy8L9b6yKB7q+C5ltXgzaN6OXAoQAEKUMADBXTPK50egCQmJqJQoUJISEhQn17JqsWzzz6Ld955B/7+/h7YpY5psu6B4phaMlcKGFMgNTUNu0HCZwAAIABJREFUU38/jA9WHrAe1dukchgmdolEkWD+d8qYvcZaUYACFKCAowR0zyudHoD88ssvuO+++6yrH08//TSmT5/uKC+PzVf3QPFYSDbcowU2HD6PwfO34dzVG8qheAF/tS+kQYXCHu3CxlOAAhSggGcJ6J5XOj0A+fTTT9GvXz/Va7L6cfz4cZQqVcqzetEJrdU9UJxQZRZBAUMKxF5JwKB527Dp7wuqfj5e+fBy6xp4umlF6y9SDFlxVooCFKAABSigSUD3vNLpAch7772H4cOHq//jrly5Mg4ePKiJRn82cqTv+++/jxkzZqhAqWzZsujduzdefPFFeHvnfDznzJkzMXHiROzduxeBgYG49957Ie2vWrXqbZXdvXs3XnrpJfzxxx/qz5o2bYqxY8eiTp06uWqY7oGSq0rwJQqYRCA5JRXvrzygblC3pNZ1SmBsx3CEBPD2dJN0M5tBAQpQgAJZCOieVzo9AJF7P2QCLwGInIS1YcMGw3b2wIED8cknn+DJJ59EkyZNsH79enzxxReQn0+ZMiXbeo8cORKjRo1C8+bN0a5dO3W54uTJk5GcnIzNmzejQoUK1vf/+usvNGzYEIULF8bgwYPVzyVwuXTpEjZt2pRpwJITmu6BklN5/HMKeILAyj2nMXTRDsTdSFbNrVgkCJ90r8/b0z2h89lGClCAAh4soHte6fQA5KuvvsJjjz2mApDq1aur1QEjpl27diEiIkIFBB999JG1is888wwmTZqEHTt2oG7duplW/fTp0yhXrpxa8VixYoX1Mw255b1WrVp45JFHMG/ePOu7HTt2VM/JnSiyyiJJVlxq1qyJ1q1bY/HixXYT6R4odleAL1DApAJHz19D/zlbse/UFdXCAF8vvP1oXXSM4u3pJu1yNosCFKCAxwvonlc6PQA5efKkmmTLCVgBAQG4cOGC+qvR0ogRI9TJXDExMahYsaK1enKBYqVKldRnZKNHj8602kuXLkWHDh0we/ZsdO/e/aZn2rRpA9mIf+7cOXXU8NWrVxEWFobOnTtj1qxZNz3bs2dPLFy4UD0bHBxsF5HugWJX4XyYAiYXyOz29M4Ny+KNh2sjwDfnzzNNzsPmUYACFKCAyQR0zyudHoBIf8gpWDIJl1UQ+e2+3P9htNSqVSu1yiGrGbcmuackMjJSrVpklhYsWIAuXbpgyZIlt7WtU6dOkFWgjRs3onHjxuoTNPm8Sz716t+//03Zyc/kcy955o477rCLSPdAsatwPkwBDxFYtPk4Xvt2N278e3t67VIF8Em3KJQL4+3pHjIE2EwKUIACHiGge17pkgBEJt+yyVpWQeSTpOjoaPj5+RmqA+XzKqmT1O3WVL9+fSQlJUE+08osSeBSr1499fmW7OWwJLn7RDagnzhxwhqcSJAin2B99913aNu27U3Zyc/kcy0JWGRFJat06tQpyP8yJvmcS1ZfpP5SXyYKUMAxAntOXsbAuVtx9Hy8KqBAgA8+fKweWtbi7emOEWeuFKAABSjgbAFTBCCC9tZbb0E2assqiEyu58yZY6ggRE7okpUO2Xh+a5IVi9jYWBw6dCjL/m/WrJlauZDPuCyb0KW9K1euVMGL5fMs+WuPHj3Uz++///6b8vvpp58gKzGZfcqV8cE33ngDb775ZqZ1YQDi7H9FWZ4nClxJSFKb03/ae8ba/P7NKmPo/dXg4+3liSRsMwUoQAEKmEjANAGI9IlMyCUQkSBEVhxkT4XskTBCyssKiNRfAhTZw5HxM60777xTnYol7fzmm2/U6gZXQIzQ26wDBfIuICu6M9YewXsr9iMlNU1l2LhiYUzqGoliIcbb55b3FjMHClCAAhTwFAFTBCB9+/a19teaNWtw+PBh9c8SiISGhqojaYsVK4b8+fPb1a/y/rRp0+x6J6uH87IHJGOe8rmVbGSX9tSoUQMvv/wyxowZo068kn/mHhAt3cVMKGAYgU1HLmDQvK2IjUu/Pb1oiD8md4lE40phhqkjK0IBClCAAhSwR8AUAYiXl9dtNwjLbw8tSQIJe5O8L+/J5YE6kpxy9e677+bqFKzsype9IXIniAQlUt+cTsGSDe3nz5/nKVg6OpV5UMBJArFxCRgyfxs2xqTfnu7tlQ8vtaqOvndX4u3pTuoDFkMBClCAAvoETBuA5JVIdwAiG8nlpKus7gHZvn07wsPD1X4OWcGRlZuSJUtm2wy5Ub1Pnz63nXgle2BkD8j+/ftRpkz6XQKWe0BkJUY+07I36R4o9pbP5yng6QJye/q4VQfx8a/pK7ySZGP6B50iEJqft6d7+vhg+ylAAQq4k4DueaVLTsGSFRBHJJ0rIFI/ORZXPumSm9DvuusurFu3Tt2E3q9fP0ydOlU14e+//1b3hMh+j5kzZ1qbNXToUJw9exZRUVHqnhM5dnjRokV4/PHH1SWEGVd5Dhw4oG6Fl/tAhgwZovKQ07Nk5UNuQpcLG+1NugeKveXzeQpQIF1g9b4zeG7hdlxJSL89vXxYID7uVh+1S4WSiAIUoAAFKOAWArrnlS4JQFavXu0w7BYtWmjLOzk5GWPHjoWsXMheDlmd6N27N1566SX4+PhkG4BIkDFu3DgcPHhQfRYm+z1k74usgGQWgO3cuVPlK0GOJDmmWPaKyCpLbpLugZKbOvAdClAgXeD4hXgMmBuN3f+k357u5+OFtx6pjccbliMRBShAAQpQwPACuueVLglADK9sggrqHigmIGETKOBSAbk9fdSyvZj35zFrPTpFlcGoR+ogvx9vT3dp57BwClCAAhTIVkD3vJIBiEkHnO6BYlImNosCThdYuvUEhn+9CwlJqarsGiVC1CdZlYoGO70uLJACFKAABShgi4DueSUDEFvU3fAZ3QPFDQlYZQoYVmD/6SsYMGcrjpy7puoY5OeN9zqEo21EKcPWmRWjAAUoQAHPFdA9r2QAYtKxpHugmJSJzaKAywTiEpLw8pJd+GHXKWsdut9RDq+2qYUAX36S5bKOYcEUoAAFKHCbgO55JQMQkw4y3QPFpExsFgVcKiDHh8/eeBRvL9uHxJT0T7JqlyqAKV3ro0KRIJfWjYVTgAIUoAAFLAK655UMQEw6tnQPFJMysVkUMITArhOXMXBeNI5fuK7qE+zvgzEdwtEmPPu7hQxReVaCAhSgAAVML6B7XskAxKRDRvdAMSkTm0UBwwhcvp6EYV/txIo9p6116nFneYxoUxP+PvwkyzAdxYpQgAIU8EAB3fNKhwQg77zzjsu6Zvjw4S4r20gF6x4oRmob60IBswrIJ1mz1v+N0T/uQ1JKmmpm3dKhmNw1EuXD+EmWWfud7aIABShgdAHd80qHBCBy0V7Gm76diSqX/jEBugcKTSlAAecJ7Dh+Cf+btxUnLqZ/khXi74OxHcPRui4/yXJeL7AkClCAAhSwCOieVzo0AJHf5jkzSdDDACRdXPdAcWY/siwKUAC4HJ+EF7/agZ/2nrFy9GpSAcMfrKluUmeiAAUoQAEKOEtA97zSFAGIBB4S7DAA+W8Y6h4ozhrgLIcCFPhPQP679vm6v/Huj/uQnJr+C52IsgUxpWskyhQKJBUFKEABClDAKQK655UOCUCaNm3qsk+w1q5d65SOMHohugeK0dvL+lHAzALbjl3EoHnb8M+l9E+yQvP7YtxjEWhRs7iZm822UYACFKCAQQR0zysdEoAYxMqjq6F7oHg0JhtPAQMIXIpPxAuLdmD1/lhrbfrdXQlDW1WHrzc/yTJAF7EKFKAABUwroHteyQDEpENF90AxKRObRQG3EkhNTcP0tTF4f+UBpPz7SVZkuYKY1IWfZLlVR7KyFKAABdxMQPe8kgGImw0AW6ure6DYWi6fowAFHC+w6cgFDJm/DaevJKjCCgT44P1OEWhVu4TjC2cJFKAABSjgcQK655UMQEw6hHQPFJMysVkUcFuBC9fkk6ztWHPgrLUNckrWKw/W4MWFbturrDgFKEABYwronlcyADFmP+e5VroHSp4rxAwoQAHtAvJJ1md/HMGYFfutp2TVKV0Ak7vUR4UivLhQOzgzpAAFKOChArrnlQxATDqQdA8UkzKxWRQwhYCckjV4/jbrxYXB/j54p31dPBxRyhTtYyMoQAEKUMC1ArrnlQxAXNufDitd90BxWEWZMQUooEXg8vUkDPtqJ1bsOW3Nr0ujshjZtjYCfL21lMFMKEABClDAMwV0zysZgJh0HOkeKCZlYrMoYCoBubhwzsajeGvZPiSmpKq2VS8egsldI1G1eIip2srGUIACFKCA8wR0zysZgDiv75xaku6B4tTKszAKUCBPArv/uaw+yTpy7prKJ7+vN958uDY6NSjjskti89QgvkwBClCAAi4V0D2vZADi0u50XOG6B4rjasqcKUABRwhcvZGMEV/vwrfbT1qzf6ReKYxuVxeyR4SJAhSgAAUoYKuA7nklAxBb5d3sOd0Dxc2az+pSgAIA5JOsxdEnMPLbPbielKJMKoQFYnLX+qhTOpRGFKAABShAAZsEdM8rGYDYxO5+D+keKO4nwBpTgAIWgUOxcfjf3G04cCZO/cjP20vdFyL3huTLl49QFKAABShAgWwFdM8rGYCYdMDpHigmZWKzKOAxAglJKRi1bC/m/XnM2uaWtYrj/Y7hKBjo5zEObCgFKEABCtgvoHteyQDE/j5wizd0DxS3aDQrSQEK5Cjww85TeHnJTsTdSFbPli6YHxO7RCKqfKEc3+UDFKAABSjgmQK655UMQEw6jnQPFJMysVkU8EiBY+fjMWj+Vuw8cVm139srH15sVR19/68SvLz4SZZHDgo2mgIUoEA2ArrnlQxATDrcdA8UkzKxWRTwWIHE5FSMWbEfn/1xxGrQrFpRjHssAmHB/h7rwoZTgAIUoMDtArrnlQxATDrKdA8UkzKxWRTweIFVe89g6OIdkJvUJRUL8cdHnSNxZ+Uwj7chAAUoQAEKpAvonlcyADHpyNI9UEzKxGZRgAIA/rl0Hc/M34YtRy8qD/kKa3DzqhjSoqr6PIuJAhSgAAU8W0D3vNIlAcj69eutvVi3bl2EhITkulfj4uKwa9cu6/tNmjTJdV5melH3QDGTDdtCAQrcLpCckorxPx/Ex78eRlpa+p83rlhYrYaUCA0gGQUoQAEKeLCA7nmlSwIQLy8v69nzq1atQvPmzXPdpatXr8b999+v3pfz7JOT00928fSke6B4uifbTwFPEVj711k8t3A7zl1NVE0uHOSHDx+LwL3Vi3kKAdtJAQpQgAK3COieV7osALEEDDoCkJYtW1oDkJSU9Nt+PT3pHiie7sn2U8CTBGLjElQQsu7QeWuz+91dCUNbVYevt5cnUbCtFKAABShglj0guldAGIDc/u8GAxD+94ICFMiLQEpqGj759RDGrTqI1H8/yYosVxATO0eibOHAvGTNdylAAQpQwM0EdM8rXbYCIp9LSeIKiGNGoO6B4phaMlcKUMDoApuOXMCQ+dtw+kqCqmqBAB+M7RiOB+qUNHrVWT8KUIACFNAkoHte6fYByPLly9GmTRvFGxAQgPj4eE3U7p2N7oHi3hqsPQUokBeBC9cS8eLiHVi9P9aaTY87y2P4gzUR4Oudl6z5LgUoQAEKuIGA7nml2wcg06dPR//+/VXXFS1aFGfOnHGDbnR8FXUPFMfXmCVQgAJGFkhLS1OXFsrlhUkp6d9k1SpZAJO7RqJS0WAjV511owAFKECBPAronle6fQDy8MMPY9myZeoErAYNGuDPP//MI7E5Xtc9UMyhwlZQgAJ5Fdhx/BIGzd+K4xeuq6wC/bwxul0dtIssk9es+T4FKEABChhUQPe80mEByMmTJ7MkLFOmjPUY3vnz56Np06Y2c8tv4eQzq5iYGCxYsACzZs1S70oA0qdPH0ydOtXmvMz8oO6BYmYrto0CFLBP4EpCEl5Zugs/7DxlfbFTVBm8+UhtBPr52JcZn6YABShAAcML6J5XOiwAyXjS1a2qEkRYkmUzem7lLXlJPrIfxHInSG7zM8t7ugeKWVzYDgpQQI+A/Ld33qZjGPX9XtxITlWZVikWrD7JqlGigJ5CmAsFKEABChhCQPe80uEBSMZgQ7egJXiRMpo1a4Y1a9boLsJt89M9UNwWghWnAAUcKrD/9BX8b+5WHD57TZXj7+OFkW1ro0ujstaVbodWgJlTgAIUoIDDBXTPK906ALFoP/jgg+pTrLCwMId3gLsUoHuguEu7WU8KUMD5AvGJyXj92z34KvqEtfA24SXxbvu6KBDg6/wKsUQKUIACFNAqoHte6bAARPZ1ZPV51bp166x/Vrt2bYSGhtqMJJ92BQUFqWAjPDwcDzzwAOrWrWvz+57yoO6B4ilubCcFKJB7gaVbT+DVb3YjPjFFZVKucKD6JCu8TMHcZ8o3KUABClDA5QK655UOC0Cyk9J5E7rLe8SgFdA9UAzaTFaLAhQwmMDhs1fVJ1n7T8epmvl658OwB2rg6aYV+UmWwfqK1aEABShgq4DueaXLAhBpsKyQ5PUmdFvhPO053QPF0/zYXgpQIPcCCUkpGP3DPszeeNSayX01i+H9jhEoFOSX+4z5JgUoQAEKuERA97zSJQHIE088YcUbPnw4atas6RJMMxeqe6CY2YptowAFHCOwfNcpvLRkJ+ISklUBJUMDMLFLJBpWKOyYApkrBShAAQo4RED3vNIlAYhDZJjpTQK6Bwp5KUABCuRG4PiFeAyavw1ygaEkb698eL5lNQxoVhleXvlykyXfoQAFKEABJwvonlcyAHFyBzqrON0DxVn1ZjkUoID5BBKTU/HBTwcw/fcYa+OaVimCcY9HoFhIgPkazBZRgAIUMJmA7nklAxCTDRBLc3QPFJMysVkUoIATBdbsj8Xzi7bjYnySKrVIsD8mPF4PTasWcWItWBQFKEABCtgroHteyQDE3h5wk+d1DxQ3aTarSQEKGFzg9OUEDFmwDZuOXFA1zZcPGHhPZTx3XzX4eHsZvPasHgUoQAHPFNA9rzRMAHL58mVI444ePQr5+/j4eOTmFnXZ1M4EZRkVFYXo6GjUr1+fJBSgAAUMI5CckoqJvxzCpF/+QlpaerWiyhfCR53roUyhQMPUkxWhAAUoQIF0Ad3zSpcHIAsWLMDHH38MuZxQR0pJSb8Ay9OT7oHi6Z5sPwUooF9g/eFzeHbBdsTG3VCZFwjwwZgO4Whdt6T+wpgjBShAAQrkWkD3vNJlAcj58+fRpUsXrF69WmHcutqR1S3qWcnJ+/IOAxDHRKq5HrF8kQIUoEA2Aueu3sDQxTvw64Gz1qe6Ni6H1x+qhQBfb9pRgAIUoIABBEwRgFy/fh333HMPtmzZogIPS7CRm0+uMvYJA5D/NHQPFAOMfVaBAhQwqUBqaho+X3cEY1bsR1JK+jdZ1YoHY3LX+qhWPMSkrWazKEABCriPgO55pUtWQN566y2MHDnypsDDx8cHzZs3V/sWihcvjqCgoFz1ytNPP52r98z2ku6BYjYftocCFDCewK4TlzF4/lb8fT5eVc7fxwuvt62Fro3KWf//wni1Zo0oQAEKmF9A97zS6QGIfCJVpEgRXLlyxfrZVYcOHTBp0iSUKFHC/D3opBbqHihOqjaLoQAFPFzg6o1kvPbNbny97R+rROs6JfBe+3CEBvp6uA6bTwEKUMA1ArrnlU4PQNavX4+mTZtaf5v18MMP4+uvv3aNpolL1T1QTEzFplGAAgYUWBJ9Aq99uxvxiekHi5QumB8Tu9RDVPnCBqwtq0QBClDA3AK655VOD0A+++wz9OnTR/WS7NnYv38/qlatau5ec0HrdA8UFzSBRVKAAh4uEHP2KgbP34Y9J68oCW+vfHjuvqoYcE8V9fdMFKAABSjgHAHd80qnByBjxozBK6+8orTKly+PI0eOOEfOw0rRPVA8jI/NpQAFDCJwIzkFY5YfUJvULenOSmGY0LkeihcIMEgtWQ0KUIAC5hbQPa90egDy/vvvY9iwYWr1o2HDhti4caO5e8xFrdM9UFzUDBZLAQpQQAn8sv8Mhi7eiQvXEtU/Fw7ywwedwtG8RnEKUYACFKCAgwV0zyudHoDMnz8f3bp1UwFIzZo1sXv3bgeTeWb2ugeKZyqy1RSggJEEzlxJUBcXbog5b63WU3dVxLDW1eHvwztDjNRXrAsFKGAuAd3zSqcHIDExMahSpYrqFTlq9+LFi5AjeJn0CugeKHprx9woQAEK5E4gJTUNn/x6CON//gvy95JqlyqASV0iUalocO4y5VsUoAAFKJCtgO55pdMDEGld48aNsXnzZrUKsmLFCrRs2ZLdrllA90DRXD1mRwEKUCBPAtFHL2DI/O3459J1lU+gnzfeeqQOOkSVyVO+fJkCFKAABW4X0D2vdEkAsnz5crRp00YFIE2aNMHatWvZ15oFdA8UzdVjdhSgAAXyLHA5PgkvL92J5btPW/NqF1kabz1aB8H+XFnPMzAzoAAFKPCvgO55pUsCEGnLgAEDMG3aNBWEDB48GBMmTGAnaxTQPVA0Vo1ZUYACFNAmkJaWhnmbjmHU93txIzlV5VshLBCTutRH3TKh2sphRhSgAAU8WUD3vNJlAUhqaiqefPJJzJ49WwUhrVq1gpyQVbt2bU/uX21t1z1QtFWMGVGAAhRwgMCB03EYPH8rDp65qnL39c6HYQ/UgGxS9+KdIQ4QZ5YUoIAnCeieV7okAHnnnXesffbFF1/g8OHD1pvR69atq/aIFCtWDPnz57e7b4cPH273O2Z8QfdAMaMR20QBCphL4HpiCt76YS/m/XnM2rB7qhfFB50iUCTY31yNZWsoQAEKOFFA97zSJQGIl5eXNeCw2MkyuiXJikhuU0pKSm5fNdV7ugeKqXDYGApQwNQCy3edwrAlO3ElIVm1s2iIPyY8Xg93VSli6nazcRSgAAUcJaB7XunSAESCjrwEGxmRLXkxAElX0T1QHDWgmS8FKEABRwicuBiPZxZsR/TRiyp7+b3WgGaV8VzLavD19nJEkcyTAhSggGkFdM8rXRaAOKKHJJhhAMIAxBFji3lSgALuJ5CckooJP/+FKb8egmWRPbJcQUzsHImyhQPdr0GsMQUoQAEXCZgiAPnss88cxvf00087LG93ylj3QHGntrOuFKAABTIKrD90Ds8u3I7YuBvqxyEBPnivfTjahJckFAUoQAEK2CCge17pkhUQG9ppiEdkNUVO5poxYwaOHz+OsmXLonfv3njxxRfh7e2dYx3nzp2LSZMm4cCBA5BPxOQG+L59+6o8ZB+MJfXq1QuzZs3KND8pt0wZ+y/W0j1QcmwsH6AABShgYIHzV2/gxa924pf9sdZadmlUFq8/VBv5/XL+77mBm8aqUYACFHC4gO55JQOQbLps4MCB+OSTT9RxwXJh4vr16yGndsnPp0yZkm1njx49Gq+++qo6Xvjhhx9WAciSJUuwZs0avPDCC/jggw9uC0AkCMkYmMgD7du3R2Cg/Z8K6B4oDh/ZLIACFKCAgwXkv8Ofr/sb7y3fh6SU9INPqhQLxqQukahZsoCDS2f2FKAABdxXQPe8kgFIFmNh165diIiIUJckfvTRR9annnnmGbWqsWPHDsiRwVklOUa4fPny2LRpk3Wjvdx9Ur9+ffz999+4dOnSbQFIUlISfHz03N6re6C4778yrDkFKECBmwV2/3MZg+dvw5Fz19Qf+Hl74aUHqvPOEA4UClCAAlkI6J5XMgDJAnrEiBGQ+0piYmJQsWJF61NHjhxBpUqVIPeNyCpHVknuMGnRogWWLVt20yOyIiLBzcmTJ28LQBITE3H9+nUEBwffthJi778RugeKveXzeQpQgAJGFrh6Ixmvf7sbS7f+Y61mk8ph6s6QUgXtv4PKyG1l3ShAAQrkVUD3vJIBSBY9IoGCrHKcPn36tieKFy+OyMhIrFixIsv+fOihh7B8+XL1qdUjjzyiPsFavHgxJLCZPHkyBgwYcFsAEhISgri4OHUBY+vWrTF27FhUrlw5V2NG90DJVSX4EgUoQAGDC3y/4yRGfL3LemdIgQAfvN2uLh6OKGXwmrN6FKAABZwnoHteabgA5OrVq9i/fz/Onz+vPlOSjeD33XefuhndmUk+r/Lz80N0dPRtxcpnVPK5lKxkZJUkcHniiSfw888/Wx8JCAjAp59+iu7du9/02ssvv6zyi4qKgr+/PzZu3Kg+85KAZMuWLepTruzSqVOnIP/LmPbt26fKkfpLfZkoQAEKUCBzgVOXr2Po4h1Yd+i89YFH65XCqEfroECAL9koQAEKeLyAKQMQCTqmTp2KhQsXYvv27ZC9EhnTqlWr0Lx589s6f/78+daJd8mSJdGlSxdtA0RWHmSlQzae35pkQ3psbCwOHTqUZXkXL17EsGHDIJ9VPfjggyrA+PLLL/HLL79A6t2xY8ds6yptvv/++9GzZ0/MnDkz22ffeOMNvPnmm5k+wwBE25BgRhSggIkFUlPT8MX6vzFmxX4kJqf/f1Dpgvkx/vF6aFSxsIlbzqZRgAIUyFnAdAHInDlz1EbvK1euqNbLp0oZk1wumFUAIkfkyuqBJF9fXxWMFCpUKGdFG57IywqIrNo0atRIfT61aNEia2nStqZNm6pjeeV4XfnUKrskKyLSpoz7RTJ7nisgNnQoH6EABShgg8CB03F4ZsE27D8dp572khvU76mMZ+/jDeo28PERClDApAKmCkCee+45TJw40Rp0SLBhCUAsf59dAHL58mWULl0a8fHx6qSpjz/+GP369dPS9XnZAyJH7cqKjazoPPbYYzfV58MPP8TQoUOxefNmNGjQINu6yhG8soldVlHsTboHir3l83kKUIAC7iqQkJSCD1YewIw/jlibEF4mFBMer4dKRYPdtVmsNwUoQIFcC+ieV7psD4isXsgnSpIswYasXtx7771qz8P48eOtf5bVCog80KFDB3z99dcqj7Zt2+InXzb0AAAgAElEQVSbb77JNW7GF+WUq3fffTdXp2DJJ1Zdu3bFvHnzbvssbMyYMWrVZsOGDbjjjjuyrWt4+P+3dx/gVRRrH8D/phISOoTea0IvKhfw0kGkiqjIRUTlimBBpYMKSFMQvRZEBD5QiopgAbyACAjSy6VKD6GEDiFACKTyPe/Ec0io2ZPJ2T17/vM8PmqyuzPzmzkw75mdmWqQV7lktsRo0t1RjObP6ylAAQp4usCag+fR94ftOHM59QT1IH9fvNMmHHKAofydw0QBClDAWwR0jytNCUBk7UR4eLhaYC5JThUfNWoU5IwNWYQtSQ7kc/wBf68AZObMmWqdhKRcuXIhOjpay18MsgOW7HR1t3NAZK2KBAiytiMiIkLlLetQJG3btk0t/H700UfVTliOlJSUpGY95BWsc+fOqe12r169qs7+cNTbca3MnnTu3FnN6Mj6GKNJd0cxmj+vpwAFKGAHgYtXEzD4x11Y8tfNHRGbhxfE+x2rIl9I6t9XTBSgAAXsLqB7XGlKACK7Q82ePdsZaMgMhsxepE0ZDUAOHjyIihUrqlslYJHBfbly5bT0g5dffhmTJ09WJ6HXr18fa9euVSehpw0K5FBBOSfk1sXisg3vr7/+ioYNG6rTzCX4kPUuEpwMGzYMsnBckgQyrVu3RocOHVC+fHm185bMjsjsSfHixbFx40a1GN5o0t1RjObP6ylAAQrYRUBto74lCsMX/oW4hNQvzgrkCMT4TtXQqKJ7d2i0iynrQQEKeJaA7nGl2wMQmTHImzevWrchSQb5EydOvK0VMhqAyF8MOXPmVDMJEoAsWLBADeh1JAka5CyOqVOnIioqCsWKFUOPHj0wYMAA54nldwtA4uPj1ZoU2flKDjOUdRwy6yPnf8gzHEm26+3bt69aEyKLzcVHAg8JyOTMkPz587tUFd0dxaVC8CYKUIACNhI4cv4q3vh+O7Yfj3HWqnu9UhjUqhKy+fvaqKasCgUoQIH0ArrHlW4PQFavXo1GjRqpWt1rxiKjAYg8p0KFCmpLXN0L0T258+nuKJ5swbJTgAIU0CWQmJyCz1YcwucrDiLl700bKxQMwX+eronwIjl1ZcPnUIACFLCUgO5xpdsDkK+//lq90iSpSJEiambhTslIAPLggw+qA/ckAJFF3rLLlLcn3R3F2z1ZfwpQgAJpBbYejVazIcejr6kfB/j6oH/LinixQWn4yN69TBSgAAVsJKB7XOn2AOTDDz9UrzBJsCDnXGzatCnTAYicuSEnhsszR48e7TwbxEbtbrgqujuK4QLwBgpQgAI2F7hyPRHDF+zB/P/d/CKtfrl8mPBkDRTKlc3mtWf1KEABbxLQPa50ewAi52D0799fBQs1atRQMxd3SkZmQGTRuayz4CtYNyV1dxRv+pCxrhSgAAWMCPy68xSG/LQLl64lqttyBfljbMeqeKxq6s6ITBSgAAU8XUD3uNLtAYgsyu7evbtqBx2vYMni89y5cyMlJUU9U7av7dSpk6e3c6bLr7ujZLpAfAAFKEABGwucunQNfefuwLqIC85adqpdDMPbVUZIoJ+Na86qUYAC3iCge1zp9gBk+fLlaN68uWormbGQNSCO8zPSNmBGZ0Bkq1vHFr7yvF27dqndprw96e4o3u7J+lOAAhS4n0BKyg1MWxOJ8Uv3IyE59UuxEnmz4+Ona6B2yTz3u52/pwAFKGBZAd3jSrcHILL9rmzDK9vNSpIT0d96663bwDMagMhhf7/99pu6PzQ0FLKtLROgu6PQlAIUoAAFMiaw5+Rl9PluGw6ejVU3+Po8gFcal8OrjcshwM8nYw/hVRSgAAUsJKB7XOn2AEQsmzVrhhUrVijWQoUK4a+//kKePOm/HcpIAPLtt9/iX//6l/PkczkM8P/+7/8s1FzmFUV3RzGvJsyZAhSggOcJXE9MxvuL92HGuiPOwlcqlAPjO1VH1WK5PK9CLDEFKODVArrHlaYEIEuWLMFjjz3mDBwaN26MRYsWIVu2m7uG3C8A+eWXX1Twce3aNchhhHL9zp07+frV3x8P3R3Fqz91rDwFKEABFwVW7j+LAfN24tyVeOdsyFvNK6BXw7LcrtdFU95GAQq4X0D3uNKUAETY5DBCOZRQ1m1ICgsLg+yQ1bJlS/X/dwtA9uzZg7Fjx0JmPxwLz+UZXbp0wcyZM93fIhbNUXdHsWg1WSwKUIAClhe4FJeIkb/uwbytN7frfaR8frVTVrE82S1ffhaQAhSggO5xpWkByKlTpyAHCMq/JckshgQS8kqW/HzBggXq5/Kzp59+WgUkmzdvVieep71e7qtYsaL6XUhICHvI3wK6OwphKUABClAgcwIr953Fm3O3IyYudQ1k9gBfDG5VCV3rlnR+GZe5HHg3BShAgawR0D2uNC0AEZ4dO3agdevWOHnypPrDV4IJR9Dh+G/H/zuCDger4/oyZcpAXumSs0CYbgro7ii0pQAFKECBzAtEXYzDwPk7sfbQze1621UvgtGPV0GObP6Zz4BPoAAFKJAFArrHlaYGIOJz/vx5dO3a1bmTleOVrHvZOYIT2c53zpw5yJcvXxZQe/YjdXcUz9Zg6SlAAQpYR0D+Dpu75ThGLdqLK/FJqmCFc2XDqA5V0DSsoHUKypJQgAIU+FtA97jS9ADE0bKyHkS25JVzQq5fv37XBvfz88MjjzyCwYMHq920mO4soLuj0JkCFKAABfQKHD4Xi54ztzq365WnP1m7GEZ2qIJs/r56M+PTKEABCmRCQPe40jIBiMMkISEB69evhyw2v3DhAi5evIjs2bOjQIECkNetGjZsiBw5cmSC0Dtu1d1RvEONtaQABSjgXgHZrvfzFYfw5aoIJKWkvoZcuUhOjOtUDZWLcLte97YGc6MABe4moHtcabkAhE2vR0B3R9FTKj6FAhSgAAXuJCCHF/aevRVHLsSpX8vhhS83LIM+TSvw8EJ2GQpQwHQB3eNKBiCmN2nWFEB3R8maUvKpFKAABSjgEJDtet9dsBu/bD/pRKlTMg8+71ILhXLdPCeLYhSgAAXcLaB7XOn2AOSnn35C3759nW5ffvklWrRo4W5H2+enu6PYHowVpAAFKGARgRX7zmDoT7tx6lLqesgcgX4Y9FglPPNgCR5eaJE2YjEo4G0CuseVbg9Axo0bh0GDBql2CwoKQnR0NAIDA72tHbO8vro7SpYXmBlQgAIUoIBTIPpqAl779n/ptuutVzYfPn2mJvKH8O9MdhUKUMC9ArrHlW4PQGSnq4EDB6pzP8qXL499+/a5V9BLctPdUbyEjdWkAAUoYBmBlJQbmL3pGD5YvA+xabbrff+JamhYoYBlysmCUIAC9hfQPa50ewAya9YsdOvWTQUgcuL5hg0b7N9qJtRQd0cxoQrMkgIUoAAFAJy6dA19vtuOTZHRTo+ONYvi3bbhyJ09gEYUoAAFslxA97jS7QHIypUr0bRpUxWAlCpVChEREVmO5o0Z6O4o3mjIOlOAAhSwikBicgom/HYAU/48jOS/t+stnjcIX3atze16rdJILAcFbCyge1zp9gBEDhkMDQ1FbGwsfHx8cPLkSfX/THoFdHcUvaXj0yhAAQpQwBWB3Scuof+8ndh76rK63c/nAbz0zzJ4vWl5Hl7oCijvoQAFMiSge1zp9gBEatm1a1fMmTNHzYLIovS0u2JlSIEX3VdAd0e5b4a8gAIUoAAF3CIQn5SM4Qv+wrebjjvzK50/GF/8qxbCCud0SxmYCQUo4F0CuseVpgQgcsp5zZo1kZSUhFy5cmHLli3qlHMmfQK6O4q+kvFJFKAABSigQ2D53jN45+fdOPn3dr1B/r4Y/FgldH24JLfr1QHMZ1CAAk4B3eNKUwIQqc3kyZPRq1cvVbESJUpg7ty5eOihh9jUmgR0dxRNxeJjKEABClBAo8DV+CSMXLQH322+ORsihxd+0KkayhYI0ZgTH0UBCnizgO5xpWkBiDTiN998g969eyMuLg6+vr54/PHH8eyzz6Ju3booUIBbDGamo+vuKJkpC++lAAUoQIGsFfhhy3G8t3APrvy9XW/2AF988EQ1tK1eJGsz5tMpQAGvENA9rjQlAEn7upUcRHj5cupiOlkT4kjBwcHq9Sx/f/8MN6zcz121Url0d5QMNwIvpAAFKEABUwTOXL6uXsn6bc8ZZ/6tqhTCiPaVEZojmyllYqYUoIA9BHSPK00JQGT3KwkWbty4kS7okP/PTJJnJicnZ+YRtrlXd0exDQwrQgEKUMDGAvL36KyNx/Dewr+QmJz6d2rObH4Y16kaHq1S2MY1Z9UoQIGsFNA9rjQ1ANEJ5QhmGICkquruKDrbis+iAAUoQIGsFThw5goGzt+JbcdinBk994+S6P9oJYQE+mVt5nw6BShgOwHd40rTApCsaBnOgNxU1d1RsqK9+EwKUIACFMg6ATmwcPraSLy/eB+S/j68sGjuIIx+vAoaVeT5W1knzydTwH4CuseVpgQg9msW69VId0exXg1ZIgpQgAIUyIjA1qPR6Dt3B45ciHNe/mrjcnijWXn4+fpk5BG8hgIU8HIB3eNKBiA27VC6O4pNmVgtClCAAl4hcD0xGZ8sP4ivVh+GzIxIqlwkp9opq0rRXF5hwEpSgAKuC+geVzIAcb0tLH2n7o5i6cqycBSgAAUokCGBDYcv4NU5/8P52AR1va/PA3ilcTn0aVpe/TcTBShAgTsJ6B5XMgCxaT/T3VFsysRqUYACFPA6gQux8Rj16178tO2Es+7/KJNPzYaUyJfd6zxYYQpQ4P4CuseVDEDub+6RV+juKB6JwEJTgAIUoMBdBf7Yfxb9ftjhnA3J5u+Dfi0q4oX6peHD2RD2HApQII2A7nElAxCbdi/dHcWmTKwWBShAAa8WOH3putqud9WBc06HxhULYMJTNZA3OMCrbVh5ClDgpoDucaWlApAdO3Zg1apV2LNnDy5cuICYmBh1sOCECRNQs2ZN9gMDAro7ioGseSkFKEABCniQgJyjJa9jvbdoD2LiElXJJfh4t0042tcoku7AYA+qFotKAQpoFNA9rjQ9AJE/+GbMmIHx48dj//796agchwsuW7YMTZo0uY2xd+/e+OOPP9TPS5cujV9//VUjtWc/SndH8WwNlp4CFKAABe4ncDLmGnrN/h92HL95eGHTSqGY8FR15M7O2ZD7+fH3FLCzgO5xpakByNmzZ/HMM8+oIEKCDUlymKCktP9/twBk8eLFaN26tfO+rVu3okaNGnZu/wzXTXdHyXDGvJACFKAABTxWIDE5BdPWROLjZQcQn5Si6iGHF77/RFU8Ur6Ax9aLBacABTInoHtcaVoAIq9Y1a1bF4cPH1bBhgQejqDDEYjcbwZEritTpgyOHDmi7h84cCDGjBmTOWGb3K27o9iEhdWgAAUoQIEMCBy9cBWvf7c93WzIE7WK4d224cgV5J+BJ/ASClDATgK6x5WmBSCNGjXC6tWrnYFH9uzZ8eKLL6JVq1YoUaIEqlSpotpNAou7zYDI7/v27YuPP/5YXVe7dm1s2rTJTu3tcl10dxSXC8IbKUABClDAIwXik5Ix4bcDmPrnYfx9diFK5M2OiV1qoWoxHl7okY3KQlPARQHd40pTApB58+bhqaeecr5uVb16dSxYsADFihVzsvj4+Dh/f68A5Pfff0eLFi3Ufb6+vmrhenBwsIu89rlNd0exjwxrQgEKUIACRgR2RsVgwLyd2Hf6SurftT4PoEeD0nijWQUEBfgaeRSvpQAFPFRA97jSlACkTp06kIpIKlq0KLZt24b8+fOna5KMBiDR0dHOe2UWRGZAZCbE25PujuLtnqw/BShAAW8WkNmQEQv3YM7GY06GkvlSZ0OqFOVsiDf3DdbdOwR0jyvdHoDIwvPChQs7W2vKlCl44YUXbmu9jAYgcqM878yZM2rGZNasWWphu7cn3R3F2z1ZfwpQgAIUAOTwwqE/7caJmGuKI8DPB282q4B/P1Iafr4+JKIABWwqoHtc6fYA5Oeff0bHjh1V8/j7+0NmMO70ypSRAETWi8jZIRKAfPLJJ3j11Vdt2vwZr5bujpLxnHklBShAAQrYWSAuIQmjft2bbjakcpGcGN+pOsKL5LRz1Vk3CnitgO5xpdsDkMmTJ6NXr14qWKhQoQL27t17x8Y0EoD84x//wMaNG9UzR48ejUGDBnltB3FUXHdH8XpQAlCAAhSgQDqBBTtOYviCvxB9NcE5GzKiXWV0frA4Dy9kX6GAzQR0jyvdHoCMHTsWQ4cOVX84PfTQQ1i/fn2mAxBZxL5r1y71zHHjxqmdsbw96e4o3u7J+lOAAhSgwO0CF2LjMXLRHvy8/aTzl4+Uz48xj1dF8bzZSUYBCthEQPe40u0ByMSJE/Haa6+pYCEsLAy7d+/OdABSqFAhyNoSeea0adPQvXt3mzS369XQ3VFcLwnvpAAFKEABuwvM3xqFt3/ejWuJyaqqQf6+GN4uHE/V4WyI3due9fMOAd3jSrcHIHPnzkXnzp1Va4WEhODSpUt3nKrN6CtYR48eRenSpdXzJABZtGiROkvE25PujuLtnqw/BShAAQrcW+B4dByG/LQLfx4877zw0cqFMKJ9ZRTMmY18FKCABwvoHle6PQCRV6XklSlHwCCvYMmrWLemjAYgH330Efr166dul3tkJiRv3rwe3MR6iq67o+gpFZ9CAQpQgAJ2Frhx4wZ+2BKFdxfsxvXEFFXVHNn8MPSxMDzNtSF2bnrWzeYCuseVbg9ApH1k21wJFCR169YN06dPdykAuXr1KsLDwxEVFQX5Q69atWrYvn27zbtAxqqnu6NkLFdeRQEKUIACFAAOnb2C/vN2YtuxGCdH66qFMapDFeQJDiARBSjgYQK6x5WmBCCvvPIKJk2apOhl1kJOOm/cuHG6psjIDEiXLl3w3Xffqfu4AD19T9bdUTzsc8LiUoACFKCAyQLJKTcwa8NRfLBkH+ISUteG5AsOwLttw9GuehHulGVy+zB7ChgR0D2uNCUAOXHiBMqVK4eEhAQ1c5EzZ07MmzcPzZo1c1rcKwC5cOECevTogQULFqjr5Rny2tWxY8eQPTt33RAT3R3FSCfltRSgAAUoQAGHQMS5WPT5bht2n7jsRGlZuSDGPVEdubL7E4oCFPAAAd3jSlMCEHEeOXIkhg0bpr4BkQBC/v3444+rU9FlTUhoaKhqDvn577//jrp162LLli346aef1E5XsbGx6j7HNfIal7zOxZQqoLuj0JUCFKAABSjgqkBicgqmrYnEx8sOID4pdW1IwZyBGN62Mh6tUoizIa7C8j4KuElA97jStABEvJ5++mn88MMP6YKQtI6OwERmSC5fvvnNSdrAQ/67Z8+ezle63NQOls9Gd0exfIVZQApQgAIUsLxA5PmreP3bbdh14pKzrM3CQjG2YzUUyBFo+fKzgBTwVgHd40pTA5DExET06dMHX375pfPbD0dwca8GdsyayDVvvvkmPvzwQ357cguY7o7irR841psCFKAABfQKxCclY9IfEfhiZQQSklNnQ/KHBGLM41XQonIhvZnxaRSggBYB3eNKUwMQh4is/3j33Xexb98+9SMJMO6WHAFK2bJlMWbMGDz55JNaYO32EN0dxW4+rA8FKEABCpgrIGtDBs3fic1HLjoL0iK8IN5rXwWFcvHcEHNbh7lTIL2A7nGlJQIQqaIEFgsXLlT//PHHH4iIiLit7WVdSNOmTdGmTRs89dRT8PX1Zf+4i4DujkJoClCAAhSggG6BpOQUfL7yECauPITE5NR1nbmC/PF+x6poVbWw7uz4PApQwEUB3eNKywQgt3rI61my29XFixfVzlYFChTgDlcGOo3ujmIga15KAQpQgAIUMCRw8MwVDP5xF7YcvTkbIqeoD29XmbMhhiR5MQWyRkD3uNKyAUjW8HnPU3V3FO+RY00pQAEKUMAMgZSUG/i/tZHq3BDHbEhIoB/ea18ZHWsVM6NIzJMCFPhbQPe4kgGITbuW7o5iUyZWiwIUoAAFLCaw/7TMhuzE/9Keol6tMIa1CUdoTq4NsVhzsTheIqB7XMkAxKYdR3dHsSkTq0UBClCAAhYUkNmQWRuPYvSve53nhuQI9MOARyuia92S3PnSgm3GItlbQPe4kgGITfuL7o5iUyZWiwIUoAAFLCwgsyED5+/E9uMxzlLKuSFjHq/K2RALtxuLZj8B3eNKBiD26yOqRro7ik2ZWC0KUIACFLC4QHLKDXy76Rg+WLwPV+KTVGkdsyFdHi4JX5+7b91v8aqxeBTwGAHd40pLBCCyBe+2bduwdetWHD16FJcuXUJcXJzamtdIkvNDpk2bZuQW216ru6PYFooVowAFKEABjxA4Hh2HPt9tS7c25KFSefHJMzVQOFeQR9SBhaSApwroHleaGoDEx8dj/Pjx+Oqrr3DixIlMtYkEKxKAJCcnZ+o5drlZd0exiwvrQQEKUIACnisgsyFzNh7FuCX7082GvNWiAp6tWxJ+vj6eWzmWnAIWFtA9rjQtAJGDBuVAwQMHDqSb6bjXKej3ahcGIOl1dHcUC38mWDQKUIACFPAygbOXr6PvDzvw58HzzpqHF86Jj56ujkqFcnqZBqtLgawX0D2uNCUAkQMGH3roIURGRioxR9Bh9JWrW7k5A3JTRHdHyfquzRwoQAEKUIACGReQnbLmbDqGcUv24fL11LUhgX4+6NuiAp6vXxr+nA3JOCavpMB9BHSPK00JQPr06YPPPvssXeARGhqKTp06oXbt2ihYsCCCg4Nd6gwNGzZ06T673aS7o9jNh/WhAAUoQAF7CJyPjceoRXvw8/aTzgpVKBiidsqqUyqvPSrJWlDAZAHd40q3ByCy7iNPnjyQfztem+rXrx9GjhyJgIAAk3ntk73ujmIfGdaEAhSgAAXsKPDzthMYsfAvXIxLVNV74AGgV8OyeL1peWTz97VjlVknCrhNQPe40u0ByIoVK9CsWTPn7MdLL72ESZMmuQ3QWzLS3VG8xY31pAAFKEABzxW4eDUB45buV9v2OlKpfNkxskMVPFK+gOdWjCWngMkCuseVbg9ApkyZgp49eypGHx8fREVFoVChQiaz3jl72VFLdumaOnUqjh8/juLFi6NHjx7o378/fH3v/23K7Nmz1atm+/fvV7M95cqVgwRc8gype9q0e/duDBgwAGvWrFE/btCgAcaNG4cqVaq4ZKO7o7hUCN5EAQpQgAIUMEFg7aHz6Dt3B05fvu7M/cnaxTC8XWUEB/qZUCJmSQHPFtA9rnR7ADJ27FgMHTpUzYCULVtW7YJl1dS7d281O/P888+jXr16WLduHaZPnw75+cSJE+9Z7NGjR+Ptt99Gy5Yt0a5dOxWAzJ8/HytXrkTfvn3x4YcfOu8/ePAgHnzwQeTNmxevvfaa+vmnn36KmJgYbNq0CeXLlzdMpLujGC4Ab6AABShAAQqYKHDleiI+XnYQM9ZFIuXvY8WK5g7CsLbhaB5e0PkmholFZNYU8BgB3eNKtwcgEyZMUDMIEoDITljr16+3JP6uXbtQvXp1FRB88sknzjI6FtDv2LEDVatWvWvZZVF9yZIlVQDh2OUrJSUFtWrVwpEjR1Rw4Uiy+H7JkiXYu3evmmWRJDMuYWFhaNWqFX744QfDRro7iuEC8AYKUIACFKCABQR2n7ikDjCMOHfVWZomlULx/hNVEZojmwVKyCJQwPoCuseVbg9AZBbgySefVIPyChUqqEG3FZPM0owZMwaHDx9G6dKlnUWUrYPLlCmDIUOGQGY57paCgoLQtGlTLFq0KN0lMiMiwc3Jk6m7dcTGxiJfvnzo3Lkzvv7663TXPvfcc/j+++9x/vx5hISEGGLS3VEMZc6LKUABClCAAhYSiEtIwqfLD2HamsNITE6dDskfEoB321ZG22qFORtiobZiUawpoHtc6fYA5NSpUyhWrJh6JSkwMBDR0dGQwbrVkgQKMstx+vTp24om2wTXrFlTzVrcLckhi4sXL1avWrVv317VV2YyJLD5/PPP0atXL3WrzADJ613yqtfLL7+c7nHyM3ndS66pW7euISLdHcVQ5ryYAhSgAAUoYEGBQ2djMeSnXdgUGe0s3cOl82JUhyooXzCHBUvMIlHAGgK6x5VuD0CEsXnz5li+fLn6xkG+4ZdXkKyW5PUq2RZ469attxVNXqNKTExUMxl3SxK4PPvss/j999+dl2TLlg2yCL9r167On8mMkNR/wYIFaNu2bbrHyc8keJk3bx6eeOKJu+YlQZ38kzbJzJLkI+WX8jJRgAIUoAAFKAAkp9zAFysP4bOVh5CQlKJI5ADDd9qEo8tDJeDj8wCZKECBWwRsEYDIugj51l9mBSpVqgSplMyGWCnJAnmZ6ZCF57cmKfvZs2dx6NChuxb54sWLGDhwIBISEvDYY4+pgOWbb76BbEP87bffOoOumTNnolu3bli6dClatGiR7nm//fabWsQu16QNWm7NdPjw4RgxYsQdy8IAxEq9imWhAAUoQAGrCByPjsN7i/Zg2Z4zziJVL54bI9pVRo3iua1STJaDApYQsEUAIpJpd8Pq0KED5syZY6kgJDMzILJ9ryywlyBm7ty5zo4jAZdsryvb8soic3n1jDMglvhcsRAUoAAFKOClAvO3RuHdX3bjakKyEpADDF/6Zxn0a1ER/r7pt8z3UiJWmwJqsqB27dra3qwx5RUsRzvK6efy7b2k8PBwjBo1Sm1Z69g1ysz2zswaENlqt0mTJur1sqeeeipdNWQXMDn5ffPmzahTpw7XgJjZyMybAhSgAAUoAOD0pesYu3gvftmeukGMpPKhIRjWtjIalM9PIwp4vYClA5AXXnjBcAOtXr1a7TSV+q3DA8iRI4eKsGQbW6OL0+X+adOmGS7DnW6QXa5klsaVXbDkFasuXbqoWZ1nnnkm3eM/+OADDBo0yLmw/H67YH333cgzt44AACAASURBVHe4cOECd8HS0qp8CAUoQAEKUODuAusOncdbtxxg+FjVQhjdoSryBAeQjgJeK2DpAERO93Z19kJeT3IkV54h98t98vqTjiQ7YMlOV3c7B2T79u2oVq2aWtsRERGBXLlyoXDhwirrbdu2qYXfjz76qNoJy5GSkpLUrIe8gnXu3DlnUCELzGUNyL59+9QOYZIc54DITIy8pmU06e4oRvPn9RSgAAUoQAFPFIiJS8CE3w5g9sajzgMMC+XMphapSzDiyhjFEx1YZgqkFdA9rtT6ClZmApDMNrPuAETKI9viTp48WZ2EXr9+faxdu1adhN6zZ098+eWXqshyqKCcEyJndsyYMcNZDdmG99dff0XDhg3RsWNHSPAxa9YsFZwMGzbM+eqZ3CABiawZkfNAXn/9dfUMOQldZj5kwX7FihUN8+juKIYLwBsoQAEKUIACHiyw5+RlDP5xJ3ZEXXLWglv2enCDsuiZEtA9rtQegGSqdpm8WecMiBRFgoZx48Zh6tSpiIqKUrMTPXr0wIABA+Dn53fPACQ+Ph5ffPGF2vlKXuOS3bBknYuc/yHPuDXt3LlTPVeCHEmyWF1e15JZFleS7o7iShl4DwUoQAEKUMCTBWSb3s9WHMTkVYeRkJy6ZW82fx8MeSwM/3q4JHy5Za8nNy/LbkBA97hSawCyatUqA1XJmktlxoEJ2ncroCkFKEABClDAWwWOXYjDyF/Tb9kbVjgnhrUNR90y+byVhfX2IgFLByBe1A6Wr6rujmL5CrOAFKAABShAgSwUkFe9522NwoiFexAbn+TM6YX6pTHg0YrI5u+bhbnz0RQwV0D3uFLrDIi5NMw9rYDujkJdClCAAhSgAAWAc1fiMX7pPvywNQqO/XNK5M2Ooa3D0CK8IBeps5PYUkD3uJIBiC27CV/BsmmzsloUoAAFKGARga1Ho/H6t9txIuaas0SNKhbAuE7VEJojm0VKyWJQQI8AAxA9jrZ/iu6OYnswVpACFKAABShgUODK9UR8vuIQ/m9tJBKTU48TyBscoE5Rf/rB4lykbtCTl1tXQPe4kjMg1m3rTJVMd0fJVGF4MwUoQAEKUMDGApHnr2LgvJ3YdCTaWcvwwjnx/hNVUa1YbhvXnFXzFgHd40pLBiAnT55U52tER0fjypUr6nT0vHnzolSpUihSpIi3tHWm6qm7o2SqMLyZAhSgAAUoYHOB5JQb+L81kfhk+UHnInV/3wfwSuNy6PnPsggK4CJ1m3cBW1dP97jSMgHI1q1bMXHiRKxcuRLHjh27ayOWKFECTZo0Qe/evVG7dm1bN3ZmKqe7o2SmLLyXAhSgAAUo4C0Cskj9gyX71I5ZjlQ4VzYMfiwMbasV5iJ1b+kINqun7nGl6QHI8ePH8cILL2DFihWqqWSbu/slOXBQUrNmzdQhgcWLF7/fLV73e90dxesAWWEKUIACFKBAJgSW7D6Nd3/ZjbNX4p1PaVe9CN5uE8ZF6plw5a3mCOgeV5oagCxbtgydO3dGTEyMM/BwBBe3BiN3+3mePHnw/fffq2CE6aaA7o5CWwpQgAIUoAAFjAlcjU/C5FURmLz6MOKTbp6k3rtRObzcsCwC/HyMPZBXU8AkAd3jStMCkC1btqBx48a4evWqopQAI+3sh7xqVaBAAYSEhCA2Nhbnzp1L92qWIyCRe4KDgyGnsNeqVcukZrFetro7ivVqyBJRgAIUoAAFPEPg0NkrasvePacuOwtcLjQEYztWxYOl8npGJVhKrxbQPa40JQBJSEhAhQoVVECRNpBo1KgRXn75ZTRv3hwys3FrunjxImTWZPLkyWqtSNp7ZYH6/v374e/v79UdxFF53R2FqBSgAAUoQAEKuC4gi9R/3nYCYxfvw/nYm69ldX6wOAa3CkOu7By/uK7LO7NaQPe40pQA5KOPPkK/fv2csx6yy9W0adPQqVOnDPv9+OOPau2I7JIlsyASjEyYMAFvvPFGhp9h5wt1dxQ7W7FuFKAABShAAXcJXLqWqBapz9l4c8Od/CEBeKdNOGSNSNpXzt1VJuZDgfsJ6B5XmhKAVKxYEYcOHVKBQ7Zs2dTrUw8++OD96n7b7zdv3oyGDRsiPj5ePat8+fJqFoSJJ6GzD1CAAhSgAAWsLLDlSDQG/7gLB8/GOov5zwoFMKp9FZTIl93KRWfZvFDA4wMQOd+jTJkyzgh/+PDheOedd1xuypEjR2LYsGHqfvnWICIiQp0X4u1Jd0fxdk/WnwIUoAAFKKBbICEpBVP+PKzODpH/lpTN3wd9mlZAj0dKw9+Xi9R1m/N5rgnoHle6fQZk3rx5eOqpp1TtfX19ERUVhYIFC7qmAeDs2bPqcELHAva5c+fiiSeecPl5drlRd0exiwvrQQEKUIACFLCawJHzVzH0511Ye+iCs2iVCuXAmI5VUavE7WtirVZ+lsf+ArrHlW4PQCZNmoRXXnlFzVbITMjBgwcz3Wry6pXMfMgzP//8c/Tq1SvTz/T0B+juKJ7uwfJTgAIUoAAFrCwgX6T+vP0ERi7ai+irCaqocuxZ14dLov+jFZEzGxepW7n97F423eNKtwcgY8eOxdChQ1WwIOs+NmzYkOk2e/jhhyHrQeSZo0aNwuDBgzP9TE9/gO6O4ukeLD8FKEABClDAEwQuXk3A2MV7MXfLzZPUQ3MEYni7ymhVpRAXqXtCI9qwjLrHlW4PQNLOgJQuXVotRs9s4gzI7YK6O0pm24j3U4ACFKAABSiQcYENhy9gyE+7cPhc6nlpkppWCsV7HaqgaO6gjD+IV1JAg4DucaXbA5Bb14CcOHECoaGhLtNwDcid6XR3FJcbiDdSgAIUoAAFKOCSQHxSMib9EYEvVkYgITl1kXr2AF+81bwCutcrBT8uUnfJlTcZF9A9rnR7ABIZGYmyZcs6pxDfe+899UqWq2nMmDF4++231e3cBeumou6O4mr78D4KUIACFKAABTIncOhsLIb+tAsbI6OdD6pSNCfGPl4NVYvlytzDeTcFMiCge1zp9gBE6pj2HJCgoCCsWbMGNWvWzED101+yfft21K9fH9evX+c5ILfo6e4ohhuHN1CAAhSgAAUooE1AFqn/sDUKY/67FzFxieq5Pg8A3euVxlstKiAk0E9bXnwQBW4V0D2uNCUA+fDDDzFgwADnSei5c+fGzJkz0bp16wy3+OLFi/Hss8/i4sWLzpPQx40bh759+2b4GXa+UHdHsbMV60YBClCAAhTwFIELsfEY/ete/LjthLPIhXJmw8BWFdG+elH4SFTCRAHNArrHlaYEIHJyebly5XDy5EnFI1G9vD7VsmVL9OjRA82aNUPOnDlvo7t8+TKWL1+OKVOmYOnSpc775P7ixYurLX0DAgI0k3vm43R3FM9UYKkpQAEKUIAC9hRYc/C8Ojvk6IU4ZwVrFM+N0Y9XQeUifC3Lnq1uXq10jytNCUCEb+PGjWjSpIl6fSptEOKgldPMCxQogODgYFy9ehXnzp2DnKLuSI6gRf6dPXt2rFixAg899JB5LWOxnHV3FItVj8WhAAUoQAEKeL3A9cTUReqTV0fgemLqInU/nwfQs2EZvNK4HLIH8LUsr+8kmgB0jytNC0DE47///S+6dOkCmdmQGRBHIJLWSn7uOOXc8fO018pMybfffotWrVppIrbHY3R3FHuosBYUoAAFKEAB+wmcjLmG9xfvw4IdqW+WSJLXsoa0DkPbaoV5doj9mtztNdI9rjQ1ABE92RXr+eefx+rVqxWmI7i4l6wjIGnYsCGmT58OmS1hSi+gu6PQlwIUoAAFKEABawss23MGw37ZjZOXUt8ukVS3TF68174KKhTMYe3Cs3SWFtA9rjQ9AHFoyytZckihvEoVFXXz9M9bW6NYsWLq1a3evXvzlat7dFXdHcXSnwoWjgIUoAAFKEABJXAtIRmTVkVg8qoIxCelvpbl6/OAOjfkjWblkSObP6UoYFhA97jSMgFIWgkJQGRmRHa4unLlCnLkyIE8efJATk6XAITp/gK6O8r9c+QVFKAABShAAQpYReB4dBxGLNyD3/eecRapYM5AjGhXGS0rF8rQGydWqQvLYb6A7nGlJQMQ85k9vwS6O4rni7AGFKAABShAAe8TWLnvLIYv/CvdblnNwgrivfaVUSR3kPeBsMYuCegeVzIAcakZrH+T7o5i/RqzhBSgAAUoQAEK3ElAdsv64o8ITPrjEBKTb6hLsgf4ok/T8ni+fmkE+PkQjgL3FNA9rmQAYtMOp7uj2JSJ1aIABShAAQp4jcChs1cw5Kfd2BQZ7axzmQLB6rWsR8oX8BoHVtS4gO5xJQMQ423gEXfo7igeUWkWkgIUoAAFKECBewqkpNzAvK1ReH/JPkRfTXBe+2jlQhjWLhyFc/G1LHah2wV0jytNCUAuXbqETz75xFm79u3bo3r16obbe/v27ViwYIHzvn79+qlDCZkA3R2FphSgAAUoQAEK2EfgUlwiPv79AL5ZfwQpqW9lITjAF282r6B2zPLz5WtZ9mntzNdE97jSlABk4sSJeO2119QODEFBQTh69Cjy5ctnWOfChQsoUaKE8zT1r776Ci+++KLh59jxBt0dxY5GrBMFKEABClDA2wX2nLyMd3/ZjS1HLzopwgrnxKgOVVC7ZB5v52H9/xbQPa40JQB55JFHsHbtWhWAPPvss5gxY4bLDdytWzfMmjVLPatRo0ZYvny5y8+y0426O4qdbFgXClCAAhSgAAVuCshrWT9sPY6xi/chJi7R+YuONYtiYKtKKJgzG7m8XED3uNLtAUhsbCzy5s2L5ORk1ZQ///wz2rZt63KzLly4EPIKl6TAwEBER0erWRVvT7o7ird7sv4UoAAFKEABuwtciI1XQYisEXEkeS3rtabl8WKD0vDna1l27wJ3rZ/ucaXbA5B169ahQYMGqoIyayGvUeXOndvlBo2JiVEBjeN58vyHH37Y5efZ5UbdHcUuLqwHBShAAQpQgAL3Fth69CKGL/gLu05ccl5YpWhOtVtW7ZKpYy4m7xLQPa50ewAyffp05zqN4sWLq/UfmU2yDkROT5eARl7nkte6vD3p7ije7sn6U4ACFKAABbxJwLFb1gdL9uFCmt2yOtYqimFtKiNXdn9v4vD6uuoeV7o9AJkwYQL69++vgoUaNWpg69atmW7U2rVrY9u2beqZH374Id58881MP9PTH6C7o3i6B8tPAQpQgAIUoIBxgYtXEzBy0R78uO2E8+bc2f3xepPyeK5eKfj6PGD8obzD4wR0jyvdHoCMHTsWQ4cOVcFCnTp1sHHjxkw3grxytXnzZvXMkSNHYsiQIZl+pqc/QHdH8XQPlp8CFKAABShAAdcFdhyPwcD5O7Hv9BXnQ6oXz413WoehTim+luW6rGfcqXtc6fYAZPLkyejVq5cKFkqWLInDhw9nWr5s2bKIjIxUz5TzRV599dVMP9PTH6C7o3i6B8tPAQpQgAIUoEDmBBKSUjBzw1F8uvwgLl27uVtW+xpFMLR1GEJzcLeszAlb927d40q3ByDz58/Hk08+qYR9fX1x7ty5TC9CL1CgAFJSUtQzZ8+ejc6dO1u3Bd1UMt0dxU3FZjYUoAAFKEABClhc4NyVePVa1oIdJ50lzZHND/1bVsS/Hi7J17Is3n6uFE/3uNLtAcj+/fsRFhamZiskyRkezzzzjCsW6p45c+aga9eu6r/lmVu2bEHNmjVdfp5dbtTdUeziwnpQgAIUoAAFKKBHYGdUDN75eTd2RN3cLatykZwY+lgY6pXLrycTPsUSArrHlW4PQESxWLFiOHXqFG7cuIHw8HDs2rXLGZAYUZZZj2rVqmHv3r3qWfnz58fZs2eNPMK21+ruKLaFYsUoQAEKUIACFHBZIDnlBuZsPIpxS/fjyvUk53M6qNeywlEgR6DLz+aN1hHQPa40JQCRNRpffPGFM+jo3bs3PvvsM8PKtz7n+eefx9SpUw0/x4436O4odjRinShAAQpQgAIU0CNw9sp1jFuyH/P/F4UbN1KfGRLoh1cal8Pz9Ushm7+vnoz4FFMEdI8rTQlAjh8/jnLlyiEpKUnNXMirU926dcPnn3+O4ODg+8LGxcVBgpaZM2eqa+UZfn5+aiZEFqQzAbo7Ck0pQAEKUIACFKDA/QR2n7iEQT/uxO4Tl52XFs8bhCGtwvBolUIuvfFyvzz5+6wX0D2uNCUAEabBgwfjgw8+UB3REYTky5cPPXv2RJs2bSBne0hQ4UgSrMiZIQsXLsSUKVNw/vx5dZ8kecYbb7wBOWOEKVVAd0ehKwUoQAEKUIACFMiIgLyWNXfLcXy4dH+6QwwblMuP4e0qo1xoSEYew2ssJKB7XGlaACLrN1q1aoVly5Y5gxBHMCH/lh2ycufOjZCQEMTGxiImJgbJycmqKdIGHvLfTZs2xdKlS+Hj42OhpjK3KLo7irm1Ye4UoAAFKEABCniawOXriZi48hCmrzmChOTU3Ur9fR/Aiw3K4LUm5RAcePOLZk+rm7eVV/e40rQARBpOAgt59ernn392Tsk5got7NaxjBy25tn379vjmm2+QI0cOb+sL96yv7o5CXApQgAIUoAAFKOCKQOT5qxi+4C+sOnDOeXtojkC1be8TtYrBh6epu8Lq1nt0jytNDUAccv/5z38gJ6TLmSCSHAHGnWQdAYqc/SEnnvfp08etDeApmenuKJ5Sb5aTAhSgAAUoQAHrCcj47bc9Z/Dewj04EXPNWUDZtvedNuGoWyaf9QrNEjkFdI8rLRGASO2uXbuG6dOnY/Hixfjzzz9x+fLNxUuO2ufMmROPPPIIHnvsMXTv3h1BQUHsGncR0N1RCE0BClCAAhSgAAUyK3AtIRlT/jyMSX9E4Fpi6qv1kmQmZPBjlZA/hNv2ZtY4K+7XPa60TACSFkvWh5w5cwYXLlxQgYgEHrJAvWDBglznkcFepbujZDBbXkYBClCAAhSgAAXuK3Dm8nW1SH1emm17g/x90atRWbzcsCwC/Liu976IbrxA97jSkgGIGz1tm5XujmJbKFaMAhSgAAUoQAHTBOQ09YHzd2HvqZtvvpQtEIz+LSuhZeWC3LbXtJZJn7HucSUDEIs0rO5i6O4ousvH51GAAhSgAAUoQAERSEpOwY/bTmD80v04dyXeidK0UiiGtA5D2QLcttfsnqJ7XMkAxOwWzaL8dXeULComH0sBClCAAhSgAAWUwKW4RHz8+wHM2XjMuW2vbJD1XL1S6NuiojpZnckcAd3jSgYg5rRjluequ6NkeYGZAQUoQAEKUIACFAAg2/YO+XEX1h++4PTIFxyAVxqXQ9e6Jbk+xIReontcyQDEhEZ0R5a6O4o7ysw8KEABClCAAhSggAjItr1/7D+H4Qv/wtELcU6UMvmDMbR1GJpUCuX6EDd2Fd3jSksEINLJtm3bhq1bt+Lo0aO4dOkS4uLinCeeZ9RXzg+ZNm1aRi+39XW6O4qtsVg5ClCAAhSgAAUsKXA9MRlfrzuCSasiEBOX6CzjPysUwLttwlEulOtD3NFwuseVpgYg8fHxGD9+PL766iucOHEiU34SxEgAkpx8c0/pTD3Qw2/W3VE8nIPFpwAFKEABClDAgwVkfcgnyw/im/VHkJRyQ9XEz+cBdK9XCq83K4+c2fw9uHbWL7rucaVpAUhERATatGmDAwcOpJvpuNcp6PdqHgYg6XV0dxTrfzRYQgpQgAIUoAAF7C4QcS5Wnaa+6sA5Z1XzBgfg9Sbl0OVhrg/JqvbXPa40JQCRAwYfeughREZGKidH0CFBRGYSZ0Bu6unuKJlpF95LAQpQgAIUoAAFdAnIeHHFvrN4b9GedOtDSubLjn4tKqJNtcJcH6IL++/n6B5XmhKA9OnTB5999lm6wCM0NBSdOnVC7dq11YnnwcHBLtE1bNjQpfvsdpPujmI3H9aHAhSgAAUoQAHPFohPSsaMtUcwceUhXL6e5KxMrRK5MaxtZVQvntuzK2ih0useV7o9AJF1H3ny5IH82/HaVL9+/TBy5EgEBARYiNqzi6K7o3i2BktPAQpQgAIUoIBdBWLiEvDFHxGYse4IEpJSnNXsVLsYBjxaEaE5stm16m6rl+5xpdsDkBUrVqBZs2bO2Y+XXnoJkyZNchugt2Sku6N4ixvrSQEKUIACFKCAZwpEXYzDB0v2Y+GOk84KBAf4omfDsujxSGlkD+BBhq62rO5xpdsDkClTpqBnz56q/j4+PoiKikKhQoVc9eB9dxHQ3VEITQEKUIACFKAABTxBYFNkNEYs/At/nbzsLG5ojkC81bwCZFbEz9fHE6phqTLqHle6PQAZO3Yshg4dqmZAypYtq3bBYtIvoLuj6C8hn0gBClCAAhSgAAWyRiA55QbmbT2OCb8dwNkr8c5MKhbMgbfbhOGR8gWyJmObPlX3uNLtAciECRPQv39/FYDITljr16+3aVOZWy3dHcXc2jB3ClCAAhSgAAUoYFwgLiEJ0/6MxJerInA14eZZcU0rheLtNuEond+1TY+Ml8Sz79A9rnR7ADJ//nw8+eSTKgCpUKEC9u7d69ktYtHS6+4oFq0mi0UBClCAAhSgAAXuK3A+Nh4fLTuA7zYdw9/nGMLf9wG82KAMXm1SDiGBXB9yL0Td40q3ByCnTp1CsWLF1A5YgYGBiI6ORlBQ0H07Di8wJqC7oxjLnVdTgAIUoAAFKEAB6wnsO30ZIxftwdpDF5yFk/Uh/VtWxOM1i3J9yF2aTPe40u0BiNSrefPmWL58uZoF+f7779X5H0x6BXR3FL2l49MoQAEKUIACFKCAOQLyJfjSv85g1K97EHXxmrMQ5UJDMKpDFdQtk8+cglk4V93jSlMCkE2bNqFevXpqFqRSpUqQSslsiNVScnIyxo8fj6lTp+L48eMoXrw4evToodaw+Pr63rW4R44cQenSpe/6ez8/PyQmJjp/3717d3z99dd3vF7ylRkjo0l3RzGaP6+nAAUoQAEKUIACVha4npiMr1YfxqQ/InAt8eb6kGZhoRj8WBjKFgixcvHdWjbd40pTAhARS7sbVocOHTBnzhzLBSG9e/dWZ5Q8//zzKmBat24dpk+fDvn5xIkT79rwV69exU8//XTb7yWYGDJkCNq2bYsFCxbcFoBIECJbE6dNHTt2RPbs2Q13Mt0dxXABeAMFKEABClCAAhTwAIFzV+Lxn98PYM6mY7hxI7XAAX4+eLFBafT8Zxnkzs6DsnWPK00LQKRx5fTz4cOHq4YODw/HqFGj0K5dO+chhWb22V27dqF69ep47bXX8MknnziL0qdPH3z22WfYsWMHqlataqiIUr933nkHshBfAgtHcsyAyKyIzI7oSLo7io4y8RkUoAAFKEABClDAqgI7jsfgw9/248+D551FzJPdH0Nbh6v1Ib4+D1i16FleLt3jSlMCkBdeeMEJtXr1ahw+fFj9v6wJyZEjB2rXro3Q0FDDi9Pl/mnTpmlpBDmrZMyYMapsaV+nioyMRJkyZdRMxujRow3lJbt+yaL7kydPIiDgZjTtCEASEhJw7do1hISE3DYTYigjQL3WJo5bt25FrVq1jN7O6ylAAQpQgAIUoIDXCcjygP/uOo0x/92LEzE314eUDw3BiPaVUa9sfq8zkQrrHleaEoDIa0YSLKRN0uCOdOvvMtLScr/cJ+s2dKSWLVuqWY7Tp0/f9riCBQuiZs2aWLJkSYazkte36tevj1dffVXNoKRNjgBEgq8rV66owKtVq1YYN26cOqzRlaS7o7hSBt5DAQpQgAIUoAAFPFEgISkFM9ZFqq17ryemOKvQsnJB9G1RERUK5vDEarlcZt3jSssEIC6L/H2j7gBEXq+SWQqZQbg1yYyCvC4lr2llNPXs2RNfffUVNm/ejDp16qS7bdCgQep5MmMhi/E3bNigghQJSLZs2YKSJUveMxvZ2lj+SZvkfJWuXbtyBiSjDcTrKEABClCAAhSgwC0Cpy9dx+crD2L2xpvrQ+RNrH89XBJvNq+AvMHesT7ENgFIVvRwnTMgMvMgMx0yc3FrkgXpZ8+exaFDhzJUjfj4eBQqVAhFixbF7t27M3TPsmXL0KJFCzz33HOYMWPGPe+RdTQjRoy44zV8BStD3LyIAhSgAAUoQAEK3FVA1oeMW7ov3fkhObL54a3mFdC1bkn4+6bfRMhulLYIQFatWpVl7dKwYUMtz9Y5AzJ37lw8/fTT6pUq2cI3o0lmRGRmQ9aM3CtxBiSjoryOAhSgAAUoQAEKuC6w6sA5jFq0BwfPxjofUqFgCIa1rYz65ey7PsQWAYjrze6+O3WuAWndujWWLl2qzhIpXLhwhishO2UtWrQIsjjdaNLdUYzmz+spQAEKUIACFKCAHQWSklPUlr0TfjuAS9dunusm54fI+pCwwjltV23d40pT1oB4QqvILldyVklmd8GSRexykKC8TvXf//7XUNWrVauGixcvqsDFaNLdUYzmz+spQAEKUIACFKCAnQUuXk3AhGX7MWfjMaT8vZeSrA/p8nAJ9G1eEXlstD5E97iSAchdPhmyA5bsdHW3c0C2b98OCRBk8XhERARy5cp1x9mNCRMmoF+/fvj+++/x1FNP3ZabHFooZ3/cehK8XN+5c2fI4vUvv/zS8OdXd0cxXADeQAEKUIACFKAABbxA4K+Tl/D+4n3pzg/JFeSPvi0qoMtDJeBng/UhuseVDEDu8cF4+eWXMXnyZHUSumyhu3btWnUSetqg4MiRI+qckLstFpcgRWYwZCbk1iBDspZARl7RktPgy5cvr3beWr9+vToZvnjx4ti4caNaDG806e4oRvPn9RSgAAUoQAEKUMCbBNYeOo/hC/5Ktz6kTP5gDGpVCS0qF/JoCt3jSgYg9+gOSUlJauH41KlTERUVpV6l6tGjBwYMGOA8sfxeAci2bdvUIYASyEyaNOmOOUlg0rdvX7U9ryw2lxkVCTzatm0LOQwxf37XyVI4jgAAGLtJREFUFjTp7ige/alh4SlAAQpQgAIUoIAbBBKTUzBrw1F1fsiV60nOHJtWSl0fEl7EM9eH6B5XMgBxQ2c0IwvdHcWMOjBPClCAAhSgAAUo4IkCF2LjMXFlBGZuOILE5NQFInIGd9eHS+KNZuWRLyTQo6qle1zJAMSjmj/jhdXdUTKeM6+kAAUoQAEKUIACFBCBQ2evYMTCPenWh2QP8EWvhmXRs2FZBPh5xvkhuseVWRKAvPfee6b1unfffde0vK2Use6OYqW6sSwUoAAFKEABClDAkwTWRZzH2z/vxuFzV53FLp43CK83KY8nahWDj2yfZeGke1yZJQGIj48P5FRyM1JycrIZ2VouT90dxXIVZIEoQAEKUIACFKCABwnI+pAf/xeF8UsP4HxsvLPkNYrnRv+WFVGvbD7Txs/3Y9Q9rszSAOTGjb83Rb5frTT9XoIeBiCpmLo7iqYm4mMoQAEKUIACFKCAVwtcikvEpFUR+Gb9EcQl3PzivEmlULzdOgxlCoRYzkf3uJIBiOWaWE+BdHcUPaXiUyhAAQpQgAIUoAAFROD0pesYu3gvFuw4Ccd39vIm1pO1i2Ngq0rIa6GDDHWPK7MkAGnUqJFpU0grV65kr+YMCPsABShAAQpQgAIU8AiBA2euqPND1kVccJY3RzY/9GhQBi/9swyCAnxNr4dHBCCmK7EAfAWLfYACFKAABShAAQp4iIAsW1h14Jw6UX3f6SvOUhfNHYQ3m1dAhxpFTD1RnQGIh3Qks4upu6OYXR/mTwEKUIACFKAABewukJScgm83HVNniJy+fN1Z3UqFcmBEu8p4uEw+Uwh0jyuz5BUsU2SYaToB3R2FvBSgAAUoQAEKUIAC7hGIS0jCZysOYfraSFxPTHFmWr9cPhWIlAvN4Z6C/J2L7nElAxC3Np/7MtPdUdxXcuZEAQpQgAIUoAAFKCAC567EY/zSfZi7JcoJ4u/7AF5tXB59mpV3G5LucSUDELc1nXsz0t1R3Ft65kYBClCAAhSgAAUo4BD46+QlfPTbASzfd1b9qG/zCnitKQMQ9hCLCTAAsViDsDgUoAAFKEABClAgkwK//XUaMzccxdTn6iDQz327Y+keV3IGJJMdwaq36+4oVq0ny0UBClCAAhSgAAUokLUCuseVDECytr1Me7rujmJaRZgxBShAAQpQgAIUoICpArrHlQxATG3OrMtcd0fJupLyyRSgAAUoQAEKUIACVhbQPa5kAGLl1s5E2XR3lEwUhbdSgAIUoAAFKEABCniwgO5xJQMQD+4M9yq67o5iUyZWiwIUoAAFKEABClDgPgK6x5UMQGza5XR3FJsysVoUoAAFKEABClCAAgxA2Ad0CDAA0aHIZ1CAAhSgAAUoQAEK6B5XcgbEpn1Kd0exKROrRQEKUIACFKAABSjAGRD2AR0CDEB0KPIZFKAABShAAQpQgAK6x5WcAbFpn9LdUWzKxGpRgAIUoAAFKEABCnAGhH1AhwADEB2KfAYFKEABClCAAhSggO5xJWdAbNqndHcUmzKxWhSgAAUoQAEKUIACnAFhH9AhwABEhyKfQQEKUIACFKAABSige1zJGRCb9indHcWmTKwWBShAAQpQgAIUoABnQNgHdAgwANGhyGdQgAIUoAAFKEABCugeV3IGxKZ9SndHsSkTq0UBClCAAhSgAAUowBkQ9gEdAgxAdCjyGRSgAAUoQAEKUIACuseVnAGxaZ/S3VFsysRqUYACFKAABShAAQpwBoR9QIcAAxAdinwGBShAAQpQgAIUoIDucSVnQGzap9auXYsGDRpg1qxZCAsLs2ktWS0KUIACFKAABShAgawW2Lt3L7p27Yo1a9agfv36mc6OAUimCa35gNmzZ6uOwkQBClCAAhSgAAUoQAEdAvLF9r/+9a9MP4oBSKYJrfmA8+fPY+nSpShVqhSCgoLcVkhHhMyZF7eRuz0jtrHbyd2aIdvXrdymZMY2NoXdrZmyjd3K7fbMzGjfa9eu4ciRI2jZsiXy58+f6TozAMk0IR+QVkD3O4LUtZ4A29h6baKzRGxfnZrWfBbb2JrtorNUbGOdmtZ7lh3alwGI9fqVR5fIDh8Kj24ANxSebewGZBOzYPuaiO+mrNnGboI2MRu2sYn4bsjaDu3LAMQNHcWbsrDDh8Kb2suVurKNXVHznHvYvp7TVq6WlG3sqpzn3Mc29py2cqWkdmhfBiCutDzvuauAHT4UbN57C7CN7d1D2L72bl+pHduYbWx/AXvX0A6fYQYg9u6jbq/dqVOnMHnyZPTs2ROFCxd2e/7MMOsF2MZZb2xmDmxfM/Xdkzfb2D3OZubCNjZTP+vztkP7MgDJ+n7CHChAAQpQgAIUoAAFKECBvwUYgLArUIACFKAABShAAQpQgAJuE2AA4jZqZkQBClCAAhSgAAUoQAEKMABhH6AABShAAQpQgAIUoAAF3CbAAMRt1MyIAhSgAAUoQAEKUIACFGAAwj5AAQpQgAIUoAAFKEABCrhNgAGI26jtnVFycjLGjx+PqVOn4vjx4yhevDh69OiB/v37w9fX196V99DaxcbG4sMPP8TWrVuxZcsWnD59Gs899xxmzJhxW42MtK+Raz2UziOKLW06a9YsrFixApGRkQgODkblypUxePBgNGvWLF0djLSZkWs9AsqDC7l3716MGDFCfYZlW04fHx+UKVMG3bt3R69evRAYGOisnZF2M3KtB/N5ZNHl89y0aVNV9oMHD6JcuXJsY49sydRCHzlyBKVLl75jDV588UU1pnIkI59LI9eaxccAxCx5m+Xbu3dvTJo0Cc8//zzq1auHdevWYfr06ZCfT5w40Wa1tUd1HH/wyXkttWvXxqJFi+4agBhpXyPX2kPSmrXo1KkTVq1ahSeeeAK1atWCBJzymdy9eze++OILNUB1JCNtZuRaa8rYp1S//fab+hLh4YcfRrFixSCDjrVr1+Lbb79F69atsXDhQraxfZobCQkJqF69uvqS7+rVq7cFIEY+m0autRGh5ari+Hu4ffv2kD+z0yYJLuvWrWvbzzADEMt1R88r0K5du9Qfiq+99ho++eQTZwX69OmDzz77DDt27EDVqlU9r2I2L3F8fDzOnz+PokWLIikpCf7+/ncMQIy0r5Frbc5revVkIFqnTp1034Jfu3YNNWrUwLlz53D27Fn4+fnBSJsZudZ0AC8uwKuvvqq++Nm3bx8qVqzINrZJXxg7diz+85//oEuXLurfaWdAjHw2jVxrEzrLVsMRgAwdOhSjRo26azmNtJmRa82EYQBipr5N8pYPzpgxY3D48OF0U4ny2oe8DjBkyBCMHj3aJrW1ZzXuFYAYaV8j19pT0vq16tu3Lz766CMcO3ZMvSpppM2MXGt9CfuWUGZF5PXXDRs2qNkRI+1m5Fr7ClqvZkePHkV4eDg+//xzyH/Lq3dpAxAj7WbkWutJ2KtEaQMQaRdJQUFBt1XSSJsZudZMTQYgZurbJO+WLVuqWQ5ZQ3BrKliwIGrWrIklS5bYpLb2rMa9AhAj7WvkWntKWr9WzzzzDObNm4eYmBi1LsRImxm51voS9ilhXFwc5B95LWfTpk145ZVX1IzmoUOH1GDGSLsZudY+gtavibyiI7OW8nqzBB+3BiBG2s3ItdaX8ewSOgKQkJAQ9ZqsJHn16o033lCfY0cy0mZGrjVTjwGImfo2yVterwoICFALIW9N8u55YmKiegWAyboC9wpAjLSvkWutq2HfksmiZXkFq02bNpg/f76qqJE2M3KtfRWtV7Phw4erAakjyazHV199hWrVqrGNrddchksk6/MkAJHgUtbrOdo77QyIkc+mkWsNF5Y3GBKQmegXXngBHTp0QMmSJXHy5Em18Fw2EenXr5/a3Meuf04zADHUVXjxnQTKli0LmemQb2ZuTbIgXb61kW/imKwrcK8AxEj7GrnWuhr2LNmlS5fwj3/8Q+2WJDOWJUqUUBU10mZGrrWnojVrJa+/yj8XLlxQu57t3LkTsl6gUaNGbGNrNlmGSyXrtmT3uubNm2Py5MnqvjsFIEY+m0auzXBBeaE2AdlMokmTJlizZg0OHDig/ow20mZGrtVWaBcexADEBTTekl6A36Z4fo/gDIjnt+G9aiCDGJmW37x5s3odsmHDhs7LjXx+jVxrb1Fr1+7jjz/GwIEDVaAZFhbGWS5rN9c9S/f222+rXetkIJo/f/67BiBGPptGrvVgOo8uuuxg165dOxV0vvTSS7b8DDMA8eguao3Ce8r7htbQsmYpuAbEmu2io1Sydaf8RSbfjP/444/q9au0ycjn18i1OsrOZ7gmcObMGRQqVMi5AYiRdjNyrWul410ZFZDXcWQjF1kPIOdqOdKnn36qdphcvnw5SpUqpa4x0m5Grs1oWXmdXgGZxZTdRWUDH9nIx0ibGblWb6mNPY0BiDEvXn0HAflwyHQ/d8Hy3O5xrwDESPsaudZztTyn5NKusre8fJs2e/ZsdO7c+bbCG2kzI9d6jpL9SupY2Cpnvci350bazci19pOzVo22b9+uNnG5V5KNJGTxspF2M3KttUS8pzTyZZGc4STrQeRAQiNtZuRaM0UZgJipb5O8ZZpf/pC82zkg8oeoYzGkTapsu2rcKwAx0r5GrrUdosUqlJKSos4L+P7779WC5H//+993LKGRNjNyrcU4bFkcWV8XGhp6W90GDBigFq/OmDFDne1jpN2MXGtLVAtVStZtLVu27LYSzZ07Fz/88IOaBZEDKGUBs5F2M3KthThsWZTo6GjkzZs3Xd3kldn69eurQ2MjIiLUdulG2szItWaiMgAxU99Geb/88svqXUU5CV0+OHIImpy63LNnT3z55Zc2qqm9qiJ7yst2rDJYHTZsmAokO3bsqCopr+04Akcj7WvkWntpWqs2b731FmQtgKz3SPv6hqOUsqhVNo+QZKTNjFxrLRH7lebxxx9XC89lsbkMUuSzvHTpUvVqToMGDbBy5Up12CTb2F5tf6dF6Gxjz2xj+ftWttCWE88lmJTX7r7++mv1Rom8WTJo0CBnxYz82WvkWrPkGICYJW+zfOUb9HHjxqnpwqioKPVBkkGPfBPn+AvQZlW2RXXk/WE51OpOSQLI7t27q18ZaV8j19oC0aKVkEHpqlWr7lo6GZw6dkky0mZGrrUojW2KJbNbMssh74vL6faBgYGoVKkSnn76aTUjLf/vSEbazci1tsH0oIrcLQAx0m5GrvUgGo8r6rRp01TAsX//fshsiJwHIscX9OnTR30JmDYZaTMj15qFxgDELHnmSwEKUIACFKAABShAAS8UYADihY3OKlOAAhSgAAUoQAEKUMAsAQYgZskzXwpQgAIUoAAFKEABCnihAAMQL2x0VpkCFKAABShAAQpQgAJmCTAAMUue+VKAAhSgAAUoQAEKUMALBRiAeGGjs8oUoAAFKEABClCAAhQwS4ABiFnyzJcCFKAABShAAQpQgAJeKMAAxAsbnVWmAAUoQAEKUIACFKCAWQIMQMySZ74UoAAFKEABClCAAhTwQgEGIF7Y6KwyBShAAQpQgAIUoAAFzBJgAGKWPPOlAAUoQAEKUIACFKCAFwowAPHCRmeVKUABClCAAhSgAAUoYJYAAxCz5JkvBShAAQpQgAIUoAAFvFCAAYgXNjqrTAEKUIACFKAABShAAbMEGICYJc98KUABClCAAhSgAAUo4IUCDEC8sNFZZQpQgAIUoAAFKEABCpglwADELHnmSwEKUIACFKAABShAAS8UYADihY3OKlOAAhSgAAUoQAEKUMAsAQYgZskzXwpQgAIUoAAFKEABCnihAAMQL2x0VpkCFKAABShAAQpQgAJmCTAAMUue+VKAAhSgAAUoQAEKUMALBRiAeGGjs8oUoAAFvEGgd+/emDRpkrOqBQoUwLZt21C0aNEMVf9///sf6tWrh/j4eOf1H330Ed58880M3c+LKEABClDgzgIMQNgzKEABClDAlgIJCQlo0KABNm/e7KyfBBSrVq2Cn5/fPescExODWrVqITIy0nndE088gXnz5tnSipWiAAUo4E4BBiDu1GZeFKAABSjgVoGjR4+qQCI6OtqZr8xgyEzG3dKNGzfQvn17LFy40HlJ+fLlsWXLFuTMmdOt5WdmFKAABewowADEjq3KOlGAAhSggFNg8eLFaN26NSSwcCSZyZAZjTulsWPHYsiQIc5fBQUFYePGjahatSpVKUABClBAgwADEA2IfAQFKEABClhb4N1338XIkSOdhZSZDJnRkJmNtGnlypVo3rw5kpOTnT/++uuv0a1bN2tXkKWjAAUo4EECDEA8qLFYVApQgAIUcE0gJSUFjz76KJYtW+Z8QLVq1bBhwwbIDIekkydPqte1zpw547zm3//+N7766ivXMuVdFKAABShwRwEGIOwYFKAABSjgFQLnz59HzZo1ERUV5azvc889hxkzZiApKQmNGzfGmjVrnL+TYGTt2rXIli2bV/iwkhSgAAXcJcAAxF3SzIcCFKAABUwXWL9+PRo2bIjExERnWaZMmYK9e/emW5ieO3duyDa8pUuXNr3MLAAFKEABuwkwALFbi7I+FKAABShwT4FPP/0Uffr0cV4TEBAA2bLXkR544AH88ssvaNu2LSUpQAEKUCALBBiAZAEqH0kBClCAAtYWePrppzF37tw7FnLQoEGQnbCYKEABClAgawQYgGSNK59KAQpQgAIWFoiNjUWdOnWwf//+dKVs1KgRfv/9d/j6+lq49CwaBShAAc8WYADi2e3H0lOAAhSggAsCFy9eVDteHTlyJN3dw4cPx7Bhw1x4Im+hAAUoQIGMCjAAyagUr6MABShAAVsIyIGE7dq1w6JFi26rj8x8LF++XC1UZ6IABShAgawRYACSNa58KgUoQAEKWFTg1pPOby1moUKFsG3bNsi/mShAAQpQQL8AAxD9pnwiBShAAQpYVOCPP/5As2bN0p10PnXqVHzxxRdq211HkhkQmQnhWhCLNiSLRQEKeLQAAxCPbj4WngIUoAAFMipw+vRpdRCh/NuRXn75ZUyaNAmRkZFqTUhMTIzzd9wNK6OyvI4CFKCAMQEGIMa8eDUFKEABCnigQHJyMpo2bYpVq1Y5S1+7dm110nlgYKD62YIFC9ChQwfIGhFJch6I/KxNmzYeWGMWmQIUoIB1BRiAWLdtWDIKUIACFNAkMGDAAIwfP975tDx58qhXrkqVKpUuh4EDB2LcuHH3vU5TsfgYClCAAl4pwADEK5udlaYABSjgPQJGZjZkpqRJkyZYvXq1E0jOC1mzZo1zpsR75FhTClCAAlkjwAAka1z5VApQgAIUsIDAndZ2DB48GGPGjLlr6U6dOqXWipw5c8Z5Ta9evdRCdSYKUIACFMi8AAOQzBvyCRSgAAUoYEGB+Ph41KtXL93uVhk96XzlypVo3rx5ut2yZs+ejS5duliwpiwSBShAAc8SYADiWe3F0lKAAhSgQAYFXnrpJUyZMsV5deHChdX5HgULFszQE2SWZOjQoc5rg4ODsXnzZoSFhWXofl5EAQpQgAJ3FmAAwp5BAQpQgAK2E5g5cya6devmrJefn5861+Of//xnhusqu2G1bdsWv/76q/MeCT4kCJFghIkCFKAABVwTYADimhvvogAFKEABClCAAhSgAAVcEGAA4gIab6EABShAAQpQgAIUoAAFXBNgAOKaG++iAAUoQAEKUIACFKAABVwQYADiAhpvoQAFKEABClCAAhSgAAVcE2AA4pob76IABShAAQpQgAIUoAAFXBBgAOICGm+hAAUoQAEKUIACFKAABVwTYADimhvvogAFKEABClCAAhSgAAVcEGAA4gIab6EABShAAQpQgAIUoAAFXBNgAOKaG++iAAUoQAEKUIACFKAABVwQYADiAhpvoQAFKEABClCAAhSgAAVcE2AA4pob76IABShAAQpQgAIUoAAFXBBgAOICGm+hAAUoQAEKUIACFKAABVwTYADimhvvogAFKEABClCAAhSgAAVcEGAA4gIab6EABShAAQpQgAIUoAAFXBNgAOKaG++iAAUoQAEKUIACFKAABVwQYADiAhpvoQAFKEABClCAAhSgAAVcE2AA4pob76IABShAAQpQgAIUoAAFXBBgAOICGm+hAAUoQAEKUIACFKAABVwTYADimhvvogAFKEABClCAAhSgAAVcEGAA4gIab6EABShAAQpQgAIUoAAFXBNgAOKaG++iAAUoQAEKUIACFKAABVwQYADiAhpvoQAFKEABClCAAhSgAAVcE2AA4pob76IABShAAQpQgAIUoAAFXBBgAOICGm+hAAUoQAEKUIACFKAABVwTYADimhvvogAFKEABClCAAhSgAAVcEGAA4gIab6EABShAAQpQgAIUoAAFXBNgAOKaG++iAAUoQAEKUIACFKAABVwQYADiAhpvoQAFKEABClCAAhSgAAVcE/h/gpBQIYsU390AAAAASUVORK5CYII=\" width=\"640\">"
+      ],
+      "text/plain": [
+       "<IPython.core.display.HTML object>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plotting.single_curve(times, sf[0], label='time', ylabel='Incoherent Intensity')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Notice that at long times the scattering is not reliable due to low statistical sampling. A more reliable curve at times~1000 can be obtained if we use all frames, not just the first 1000:\n",
+    "\n",
+    "<code>tr = np.asarray([sel.positions for _ in u.trajectory[::1]]) # all frames in the trajectory</code>"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plot all curves sequentially"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Calculate the coherent dynamic structure factor `S(q)` in a log scale of Q-values\n",
+    "\n",
+    "We will calculate more values of Q using the first 100 oxygen atoms and all the trajectory frames. Here we calculate scattering from `Q=0.01` to `Q=10` for a total of 15 Q-values spaced in a log scale."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "common = dict(xlabel='time', ylabel='Incoherent Intensity', xscale='linear', yscale='linear')\n",
+    "for i, curve in enumerate(sf):\n",
+    "    common['label'] = 'Q={:.2f}'.format(q_mod[i])\n",
+    "    plotting.single_curve(times, curve, **common)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/nscsim/quasielastic/coherent.py b/nscsim/quasielastic/coherent.py
index f445176..f8a996b 100644
--- a/nscsim/quasielastic/coherent.py
+++ b/nscsim/quasielastic/coherent.py
@@ -1,5 +1,3 @@
-from __future__ import (print_function, absolute_import)
-
 import numpy as np
 from nscsim.utilities import (glog, map_parallel, shared_array)
 from nscsim import qvec
@@ -38,14 +36,18 @@ def intermediate_amplitudes(tr, q, bc, n_cores=None):
         ----------
         frame: numpy.ndarray
             shape = (#atoms, 3) coordinates of the system
+
+        Returns
+        -------
         """
         exponents = np.tensordot(frame, q, axes=(1, 1))  # shape=(#atoms, #q's)
         exponentials = np.exp(1j * exponents)
         return np.tensordot(bc, exponentials, axes=1)  # shape=(#q's,)
 
     glog.info('\nCalculating coherent amplitudes for one set of q vectors\n')
-    #TODO: close_pool=False is a temporary fix. See issue #31
-    amps = np.array(map_parallel(serial_worker, shared_array(tr), n_cores, close_pool=False))
+    # TODO: close_pool=False is a temporary fix. See issue #31
+    amps = np.array(map_parallel(serial_worker, shared_array(tr), n_cores,
+                                 close_pool=False))
 
     return amps
 
diff --git a/nscsim/quasielastic/incoherent.py b/nscsim/quasielastic/incoherent.py
index 77f410d..70e81e0 100644
--- a/nscsim/quasielastic/incoherent.py
+++ b/nscsim/quasielastic/incoherent.py
@@ -1,176 +1,145 @@
-from __future__ import (print_function, absolute_import)
-
 import numpy as np
-import multiprocessing
-from tqdm import tqdm as progress_bar
+import functools
+from nscsim.utilities import glog
+from nscsim import qvec
+from tqdm import tqdm
 import pathos
 
-from nscsim.utilities import namedtuplefy
-
 
-def si_vector_self_intermediate(q, tr, ns=200, nt=100, dt=1):
+def intermediate_vector_set(tr, q, s_inc, n_cores=None):
     r"""
-    Single atoms self incoherent structure factor for each vector in
-    a set of q-vectors
+    Incoherent scattering for a collection of vectors.
+
+    We sum the intensities math:`A_i` derived for each atom math:`i`.
+    The intensities math:`A_i` are self-correlations of atomic
+    amplitudes. The amplitudes math:`a_i` are complex numbers.
+    math:`a_i(\vec{q}, t) = e^{-i\vec{q}\vec(r)_i(t)}
+    math:`A_i(\vec{q}, t) = si*<a_i^*(\vec{q}, t_0)a_i^*(\vec{q}, t_0+t)>_t_0`
 
     Parameters
     ----------
+    tr: numpy.ndarrary
+        Atomic trajectory, shape=(#atoms, #frames, 3)
     q: numpy.ndarray
-        q vectors shape=(#vectors, 3) for a set of q-vectors or
-        shape=(#sets, #vectors, 3) for a list of sets of q-vectors
-    tr: numpy.ndarray
-        atomic trajectory for one atom, shape=(#frames, 3)
-    ns: int
-        Number of t0's for a given t
-    nt: int
-        Number of t's (number time lapses)
-    dt: int
-        Spacing between consecutive t's, in number of frames
+        Array of q-vectors shape=(#vectors, 3)
+    s_inc: numpy.ndarray
+        Incoherent scattering cross-section, shape=(#atoms, 1)
+    n_cores: int
+        Number of CPU's to use. All cores minus 1, or 1 if only 1 core
 
     Returns
     -------
-    numpy.ndarray
-        array of shape=(#vectors) and the structure factor is normalized
-        by the number of frames
+    np.ndarray
+        Scattering atomic intensities shape=(#q's, 2 * #frames - 1)
     """
-    pass
-
 
-def si_self_intermediate_si(q, tr, ns=200, nt=100, dt=1):
+    def serial_worker(atomic_tr, atomic_s_inc):
+        r"""
+        Calculate intensities for a single atom
+
+        Parameters
+        ----------
+        atomic_tr: numpy.ndarray
+            shape = (#frames, 3) trajectory for one atom
+        atomic_s_inc: float
+            Incoherent scattering cross-section for one atom
+        Returns
+        -------
+        numpy.ndarray
+            shape = (#q's, 2 * #frames - 1)
+        """
+        r"""Pseudocode:
+        - calculate a_i(\vec{q}, t)
+        - return self-correlation of a_i to obtain A_i(\vec{q}, t)
+        """
+        ai = np.exp(-1j * np.tensordot(q, atomic_tr, axes=(1, 1)))
+        w = np.arange(1, 1 + len(ai[0]))
+        w = 1.0 / np.concatenate((w, w[::-1][1:]))
+        result = np.apply_along_axis(
+            lambda x: w * atomic_s_inc**2 * np.correlate(x, x, "full"),
+            1,
+            ai,
+        )
+        return result[:, len(result[0])//2:]
+
+    glog.info("\nCalculating incoherent intensities for a set of q vectors\n")
+    pool = pathos.pools.ProcessPool(nodes=n_cores)
+    glog.info("\nUsing pool of {} cpus.\n".format(n_cores))
+    try:
+        pool.restart(force=True)
+        result = functools.reduce(
+            np.add,
+            pool.uimap(
+                serial_worker,
+                tr,
+                s_inc,
+                ),
+        )  # shape = (#q's, 2*#frames - 1)
+    finally:
+        pool.close()
+        pool.join()
+    return result/np.sum(s_inc*s_inc)
+
+
+def intermediate(
+    tr, q, s_inc, n_cores=None, averaging=(np.average, None, dict(axis=0))
+):
     r"""
-    Single atom self incoherent structure factor for each vector in a set
-    of q-vectors, or for each set average in a list of sets of q-vectors.
+    Incoherent scattering averaged over a collection of q-vectors
 
     Parameters
     ----------
-    q: numpy.ndarray
-        q vectors shape=(#vectors, 3) for a set of q-vectors or
-        shape=(#sets, #vectors, 3) for a list of sets of q-vectors
     tr: numpy.ndarray
-        atomic trajectory for one atom, shape=(#frames, 3)
-    ns: int
-        Number of t0's for a given t
-    nt: int
-        Number of t's (number time lapses)
-    dt: int
-        Spacing between consecutive t's, in number of frames
+        Atomic trajectory, shape=(#atoms, #frames, 3)
+    q: numpy.ndarray
+        Array of q-vectors shape=(#vectors, 3)
+    s_inc: numpy.ndarray
+        Incoherent scattering cross-section, shape=(#atoms, 1)
+    n_cores: int
+        Number of CPU's to use. All cores minus 1, or 1 if only 1 core
+    averaging: tuple
+        The Q-average to carry out over the :math:`S(\vec{q}, t)`, the
+        return value from `intermediate_vector_set`.
+        Tuple elements:
+            - averaging function
+            - positional arguments for function, excluded the first which is
+              assumed to be array :math:`S(\vec{q}, t)`
+            - keyword arguments for function
 
     Returns
     -------
-    numpy.ndarray
-        If q.shape=(#vectors, 3), returns shape=(#vectors) and the structure
-        factor is normalized by the number of frames and the sum of the
-        squares of scattering lengths. If q.shape=(#sets, #vectors, 3),
-        returns shape=(#sets) and the structure factor is normalized by
-        the number of vectors in each set and the number of frames
+    np.ndarray
+        shape = (#frames)
     """
-    kw = dict(ns=ns, nt=nt, dt=dt)
-    if q.ndim == 2:
-        return si_vector_self_intermediate(q, tr, **kw)
-    elif q.ndim == 3:
-        return np.asarray([np.mean(si_vector_self_intermediate(v, tr, **kw))
-                           for v in progress_bar(q)])
+    sf = intermediate_vector_set(tr, q, s_inc, n_cores=n_cores)
+    av_args = list() if averaging[1] is None else averaging[1]
+    return averaging[0](sf, *av_args, **averaging[2])
 
 
-@namedtuplefy
-def si_self_intermediate_spherical(q_mod, tr, ns=200, nt=100, dt=1):
+def intermediate_spherical(tr, q_mod, s_inc, n_cores=None):
     r"""
-    Single atom self incoherent structure factor for a set of q-vector moduli,
-    spherically averaged in q-vector space
-
-    The analytical spherical average is carried out.
+    Coherent scattering averaged over scattering for a collection of q-vectors
 
     Parameters
     ----------
-    q_mod: numpy.ndarray
-        Array of q-vector moduli
     tr: numpy.ndarray
-        atomic trajectory for one atom, shape=(#frames, 3)
-    ns: int
-        Number of t0's for a given time lapse `t`, i.e., number of
-        sampling starting times.
-    nt: int
-        Number of t's, i.e, number time lapses
-    dt: int
-        Spacing between consecutive time lapses, in units of number of frames
-
-    Returns
-    -------
-    namedtuple
-        Fields of the namedtuple
-        sf: numpy.ndarray, structure factor of shape (q_mod.shape, nt)
-        t: numpy.ndarray, list of time lapses where the structure factor
-            is calculated
-    """
-    co = 1e-05  # avoid later dividing by zero
-    n_q = len(q_mod)
-    nfr = len(tr)
-
-    it = 0  # current number of evaluated time lapses
-    t = 0  # first time lapse is no lapse
-    sf = np.zeros((nt, n_q))  # shape is (nt, n_q)
-    ts = []
-    while it < nt:  # cycle over all t's (+1.0 for t=0)
-        dt0 = max(1, int((nfr - t) / ns))  # separation between t0's
-        t0s = np.arange(1, nfr - t, dt0)  # from 1 to nfr-t every delta
-        dij = np.linalg.norm(tr[t0s] - tr[t0s + t], axis=1)  # distances
-        phase = co + np.outer(q_mod, dij)
-        sf[it] = (np.sin(phase) / phase).mean(axis=1)
-        ts.append(t)
-        it += 1
-        t += dt
-    return dict(sf=sf.transpose(), t=np.asarray(ts))
-
-
-@namedtuplefy
-def self_intermediate_spherical(q_mod, tr, bi, ns=200, nt=100, dt=1,
-                                n_cores=None):
-    r"""
-    Self incoherent structure factor for a set of q-vector moduli,
-    spherically averaged in q-vector space and by the sum of the
-    squares of the scattering lengths
-
-    The analytical spherical average is carried out.
-
-    Parameters
-    ----------
+        Atomic trajectory, shape=(#atoms, #frames, 3)
     q_mod: numpy.ndarray
-        Array of q-vector moduli
-    tr: numpy.ndarray
-        atomic trajectory for one atom, shape=(#atoms, #frames, 3)
-    bi: numpy.ndarray
-        incoherent scattering lengths
-    ns: int
-        Number of t0's for a given time lapse `t`, i.e., number of
-        sampling starting times.
-    nt: int
-        Number of t's, i.e, number time lapses
-    dt: int
-        Spacing between consecutive time lapses, in units of number of frames
+        Array of q-vector moduli shape=(#q-moduli, 1)
+    s_inc: numpy.ndarray
+        Incoherent scattering cross-section, shape=(#atoms, 1)
     n_cores: int
-        Number of CPU cores to use. `None` means use all cores
+        Number of CPU's to use. All cores minus 1, or 1 if only 1 core
 
     Returns
     -------
-    namedtuple
-        Fields of the namedtuple
-        sf: numpy.ndarray, structure factor of shape (q_mod.shape, nt)
-        t: numpy.ndarray, list of time lapses where the structure factor
-            is calculated
+    np.ndarray
+        shape = (#frames)
     """
-    sf = np.zeros((len(q_mod), nt))
-
-    def serial_worker(atomic_tr):
-        return si_self_intermediate_spherical(q_mod, atomic_tr, ns=ns, nt=nt,
-                                              dt=dt).sf
-    if n_cores == 1:
-        for i, atomic_tr in enumerate(tr):
-            sf += bi[i] ** 2 * serial_worker(atomic_tr)
-    else:
-        if n_cores is None:
-            n_cores = multiprocessing.cpu_count() - 1
-        pool = pathos.pools.ProcessPool(ncpus=n_cores)
-        for i, _sf in enumerate(pool.map(serial_worker, tr)):
-            sf += bi[i] ** 2 * _sf
-        pool.terminate()
-    return dict(sf=sf / sum(bi * bi), t=dt * np.arange(nt))
+    q_sets = qvec.sphere_average(q_mod)
+    glog.info("\nCalculating coherent scattering for a set of q-moduli\n")
+    sf = [
+        np.real(intermediate(tr, q_set, s_inc, n_cores=n_cores))
+        for q_set in tqdm(q_sets)
+    ]
+    return np.asarray(sf)


Test patch:
diff --git a/tests/unit/quasielastic/test_incoherent.py b/tests/unit/quasielastic/test_incoherent.py
index 23c5861..f241cb8 100644
--- a/tests/unit/quasielastic/test_incoherent.py
+++ b/tests/unit/quasielastic/test_incoherent.py
@@ -34,26 +34,18 @@ def diffusive_particles(n_atoms, jump, duration):
     return shared_array(from_array=tr)
 
 
-def test_si_self_intermediate_spherical():
-    q = np.array((0.2, 0.3, 0.4, 0.5))
-    tr = diffusive_particles(1, 1.0, 10000)[0]
-    out = incoherent.si_self_intermediate_spherical(q, tr, ns=5000, nt=100)
-    time_at_fwhm = np.abs(out.sf - 0.5).argmin(axis=1)
-    assert_almost_equal(np.log(2)/(time_at_fwhm * np.square(q)),
-                        0.5 * np.ones(len(q)), decimal=1)
-
-
-def test_self_intermediate_spherical():
+def test_intermediate_spherical():
     q = np.array((0.2, 0.3, 0.4, 0.5))
     n_atoms = 100
     tr = diffusive_particles(n_atoms, 1.0, 1000)
     bi = np.ones(n_atoms)
-    serial = incoherent.self_intermediate_spherical(q, tr, bi, ns=1000, nt=100,
-                                                    n_cores=1)
-    parall = incoherent.self_intermediate_spherical(q, tr, bi, ns=1000, nt=100,
-                                                    n_cores=2)
-    assert_almost_equal(serial.sf, parall.sf, decimal=3)
-    time_at_fwhm = np.abs(parall.sf - 0.5).argmin(axis=1)
+    seed = np.random.randint(200)
+    np.random.seed(seed)
+    serial = incoherent.intermediate_spherical(tr, q, bi, n_cores=1)
+    np.random.seed(seed)
+    parall = incoherent.intermediate_spherical(tr, q, bi, n_cores=2)
+    assert_almost_equal(serial, parall, decimal=4)
+    time_at_fwhm = np.abs(parall - 0.5).argmin(axis=1)
     assert_almost_equal(np.log(2)/(time_at_fwhm * np.square(q)),
                         0.5 * np.ones(len(q)), decimal=1)