Code/Quick_MSD/notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute MSD in Python quickly with FFT\n",
    "\n",
    "The implementation is entirely copied from this SO answer : http://stackoverflow.com/questions/34222272/computing-mean-square-displacement-using-python-and-fft#34222273\n",
    "\n",
    "Thanks to [thomasfermi](http://stackoverflow.com/users/2609987/thomasfermi) for the implementation. The original algorithm comes from : \n",
    "\n",
    "```\n",
    "nMoldyn - Interfacing spectroscopic experiments, molecular dynamics simulations and models for time correlation functions\n",
    "V. Calandrini, E. Pellegrini, P. Calligari, K. Hinsen and G.R. Kneller\n",
    "Publié en ligne: 9 juin 2011\n",
    "DOI: 10.1051/sfn/201112010\n",
    "```\n",
    "\n",
    "This notebook compares the speed of calculation of MSD with FFT and without with trackpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%load_ext snakeviz\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "\n",
    "import timeit\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "matplotlib.style.use('ggplot')\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import trackpy as tp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def autocorr_fft(x):\n",
    "    N = len(x)\n",
    "    F = np.fft.fft(x, n=2 * N)  # 2*N because of zero-padding\n",
    "    PSD = F * F.conjugate()\n",
    "    res = np.fft.ifft(PSD)\n",
    "    res= (res[:N]).real  # now we have the autocorrelation in convention B\n",
    "    n = N * np.ones(N) - np.arange(0, N) #divide res(m) by (N-m)\n",
    "    return res / n  # this is the autocorrelation in convention A\n",
    "\n",
    "def msd_fft(traj, mpp, fps, max_lagtime=100, detail=False, pos_columns=['x', 'y']):\n",
    "    \"\"\"Implementation comes from here : \n",
    "    http://stackoverflow.com/questions/34222272/computing-mean-square-displacement-using-python-and-fft#34222273\n",
    "    \"\"\"\n",
    "    \n",
    "    r = traj[pos_columns].values\n",
    "    r *= mpp\n",
    "\n",
    "    t = traj['frame']\n",
    "\n",
    "    max_lagtime = min(max_lagtime, len(t))  # checking to be safe\n",
    "    lagtimes = 1 + np.arange(max_lagtime - 1)    \n",
    "\n",
    "    N = len(r)\n",
    "\n",
    "    D = np.square(r).sum(axis=1) \n",
    "    D = np.append(D, 0)\n",
    "    S2 = sum([autocorr_fft(r[:, i]) for i in range(len(pos_columns))])\n",
    "\n",
    "    Q = 2 * D.sum()\n",
    "    S1 = np.zeros(max_lagtime)\n",
    "\n",
    "    for m in range(max_lagtime):\n",
    "        Q = Q - D[m - 1] - D[N - m]\n",
    "        S1[m] = Q / (N - m)\n",
    "\n",
    "    msd = S1 - 2 * S2[:max_lagtime]\n",
    "    msd = msd[1:]\n",
    "\n",
    "    lagt = lagtimes / fps\n",
    "\n",
    "    results = pd.DataFrame(np.array([msd, lagt]).T, columns=['msd', 'lagt'])\n",
    "    results.index = 1 + np.arange(max_lagtime - 1)\n",
    "    results.index.name = 'lagt'\n",
    "    \n",
    "    return results\n",
    "\n",
    "def get(N, pos_columns):\n",
    "    traj = np.cumsum(np.random.choice([-1., 0., 1.], size=(N, len(pos_columns))), axis=0)\n",
    "    traj = pd.DataFrame(traj, columns=pos_columns)\n",
    "    traj['frame'] = np.arange(N)\n",
    "    return traj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Generate a random trajectory\n",
    "\n",
    "pos_columns = ['x', 'y', 'z']\n",
    "mpp = 0.3\n",
    "fps = 1\n",
    "max_lagtime = 1000\n",
    "\n",
    "traj = get(N=10000, pos_columns=pos_columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>msd</th>\n",
       "      <th>lagt</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>lagt</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.179613</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.359829</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.539037</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.721459</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.903548</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           msd  lagt\n",
       "lagt                \n",
       "1     0.179613     1\n",
       "2     0.359829     2\n",
       "3     0.539037     3\n",
       "4     0.721459     4\n",
       "5     0.903548     5"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get MSD with trackpy\n",
    "msd1 = tp.msd(traj, mpp=mpp, fps=fps, max_lagtime=max_lagtime, pos_columns=pos_columns)\n",
    "msd1[['msd', 'lagt']].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>msd</th>\n",
       "      <th>lagt</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>lagt</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.179613</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.359829</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.539037</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.721459</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.903548</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           msd  lagt\n",
       "lagt                \n",
       "1     0.179613     1\n",
       "2     0.359829     2\n",
       "3     0.539037     3\n",
       "4     0.721459     4\n",
       "5     0.903548     5"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get MSD with msd_fft\n",
    "msd2 = msd_fft(traj, mpp=mpp, fps=fps, max_lagtime=max_lagtime, pos_columns=pos_columns)\n",
    "msd2[['msd', 'lagt']].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f49cc62a358>"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzYAAAJlCAYAAAAIImesAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VPW9//H3LJmsk5WwBGQLSyAhEAj7vmm1iKiIdf3Z\n1t5Hi9d66cb19v66uXQTQb0ut9X2ttcWQRRERSpLErYAWVgTwyZ7CBACIQlJJjNnfn+08qsVCiST\nObO8nv8kc5iZvCePDzPzzjnzPRav1+sVAAAAAAQxq9kBAAAAAKCtKDYAAAAAgh7FBgAAAEDQo9gA\nAAAACHoUGwAAAABBj2IDAAAAIOjZr3WFV199VaWlpUpISNBzzz0nSXrzzTdVUlIiu92uTp06ae7c\nuYqJiZEkLV++XHl5ebLZbHrkkUc0ePDg9n0EAAAAAMLeNffYTJ48WT/84Q8/ty07O1sLFizQr3/9\na3Xp0kUrVqyQJJ04cUKFhYVauHChnnzySb3++uviNDkIRGVlZWZHQJhh5uBPzBv8jZlDILhmscnI\nyFBsbOzntmVnZ8tq/etN+/btq3PnzkmSiouLNWbMGNlsNnXs2FFdunTRwYMH2yE20DY8AcPfmDn4\nE/MGf2PmEAja/BmbvLw85eTkSJJqamrUoUOHy/+WnJysmpqatv4IAAAAAPin2lRs3n33XdlsNo0b\nN85XeQAAAADghl1z8YCryc/P144dO/SjH/3o8rbk5GRVV1dfvnzu3DklJydf8fZlZWWf2205Z86c\n1kYBbhjzBn9j5uBPzBv8jZmDvy1duvTy95mZmcrMzLy+YuP1ej+3CMDOnTu1cuVK/fSnP1VERMTl\n7bm5uXrxxRc1Y8YM1dTUqKqqSn369LnifX4W4O9VVlbe0AMCWsvpdKqurs7sGAgjzBz8iXmDvzFz\n8Ke0tLQrlmmL9xrLlr3wwgsqLy9XXV2dEhISNGfOHC1fvlxut1tOp1PSXxcQePTRRyX9dbnn9evX\ny2633/ByzxQb+AtPwPA3Zg7+xLzB35g5+FNaWtoVt1+z2PgTxQb+whMw/I2Zgz8xb/A3Zg7+dLVi\n0+ZV0QAAAADAbBQbAAAAAEGv1aui+UtcXJwsFovZMQKK1+tVfX292TEAAACAgBHwxcZisXDM5j/4\nbNEGAAAAAH/FoWgAAAAAgh7FBgAAAEDQo9gAAAAACHoUmwD0/PPP6/HHHzc7BgAAABA0KDYBipXg\nAAAAgOtHsQEAAAAQ9Cg2bTRq1Ci99tprmjZtmvr166fvf//7qq6u1kMPPaT+/fvrvvvu08WLF9Xc\n3KzHH39cWVlZGjhwoGbMmKFz585Jko4fP67Zs2crIyND999/v2pqakx+VAAAAEBwCfjz2ASDVatW\naenSpWppadHNN9+svXv3asGCBerTp48efPBBvfHGG0pNTVV9fb1KSkrkcDhUVlamqKgoSdJjjz2m\n4cOHa/HixSotLdXDDz+sW265xeRHBQAAAASPkCg2nm/M9Mn92H67slW3+9rXvqbk5GRJ0ogRI5Sa\nmqqBAwdKkm699VZt3rxZ06dP1/nz5/Xpp59qwIABysrKkiSdPHlSu3fv1pIlSxQREaGRI0dq+vTp\nPnk8AAAAQLgIiWLT2kLiKx06dLj8fVRUlFJTUz93uaGhQbNnz1ZlZaXmzp2rixcv6u6779b8+fN1\n+vRpJSQkKDo6+vJtunbtqlOnTvn1MQAAAADBjM/Y+InNZtO8efOUl5enlStXas2aNVq2bJk6deqk\n2tpaNTY2Xr5uZWWliUkBAACA4EOx8ZPCwkJVVFTIMAzFxMTIbrfLarWqa9euys7O1oIFC9TS0qLt\n27drzZo1ZscFAAAAgkpIHIpmpn8838zVzj9z5swZzZ8/X1VVVYqNjdXMmTN19913S5JefvllPfHE\nE8rKytKwYcN0zz33qLa2tt2zAwAAAKHC4vV6vWaH+MyVDsFyOp2qq6szIU3g4nfSdvwO4W/MHPyJ\neYO/MXPwp7S0tCtu51A0AAAAAEGPYgMAAAAg6FFsAAAAAAQ9ig0AAACAoEexAQAAABD0KDYAAAAA\ngh7FBgAAAEDQo9gAAAAACHoUmwD0/PPP6/HHHzc7BgAAABA0KDYBymKxmB0BAAAACBoUGwAAAABB\nj2LTRqNGjdJrr72madOmqV+/fvr+97+v6upqPfTQQ+rfv7/uu+8+Xbx4Uc3NzXr88ceVlZWlgQMH\nasaMGTp37pwk6fjx45o9e7YyMjJ0//33q6amxuRHBQAAAASWC9U1Wrlyw1X/3e7HLCFr1apVWrp0\nqVpaWnTzzTdr7969WrBggfr06aMHH3xQb7zxhlJTU1VfX6+SkhI5HA6VlZUpKipKkvTYY49p+PDh\nWrx4sUpLS/Xwww/rlltuMflRAQAAAOY78elxvbdlvzZ7kjXG6rrq9UKi2Nzxpwqf3M97D2S06nZf\n+9rXlJycLEkaMWKEUlNTNXDgQEnSrbfeqs2bN2v69Ok6f/68Pv30Uw0YMEBZWVmSpJMnT2r37t1a\nsmSJIiIiNHLkSE2fPt0njwcAAAAIRoZh6JOdFVqxu0oVStCXYqx6eXJ3JaXmXPU2IVFsWltIfKVD\nhw6Xv4+KilJqaurnLjc0NGj27NmqrKzU3LlzdfHiRd19992aP3++Tp8+rYSEBEVHR1++TdeuXXXq\n1Cm/PgYAAADAbG63W9s27tCKw426KIfu6Bih704eqKjY6GveNiSKTTCw2WyaN2+e5s2bp5MnT+rB\nBx9Uenq6xo0bp9raWjU2Nl4uN5WVlayKBgAAgLDR1Nio9etLtPK0TU616M7e0RoxLkd2+/XXFRYP\n8JPCwkJVVFTIMAzFxMTIbrfLarWqa9euys7O1oIFC9TS0qLt27drzZo1ZscFAAAA2l1tzQW9tSxP\n/7KkTCVnXHp8SIJ++cg4jZk0/IZKjcQemzb7xz0rV9vTcubMGc2fP19VVVWKjY3VzJkzdffdd0uS\nXn75ZT3xxBPKysrSsGHDdM8996i2trbdswMAAABmqDp6Uis3VajAnayRVo+eHpeq7n1y23SfFq/X\n6/VRvjarrKz8wjan06m6ujoT0gQufidtx+8Q/sbMwZ+YN/gbM4frdbDsoJYXH9VOJWl65AXdPilL\nKZ073tB9pKWlXXE7e2wAAAAAtBvDMLRz2x4t/6RGJyxxmplk19wp/RQbH+fTn0OxAQAAAOBzhmGo\nePNOLTlQr2bZdGfXKI2fNFiOSEe7/DyKDQAAAACf8Xg82r5pp5YcuiRDFt3bK1qjJgyVzWZr159L\nsQEAAADQZh6PR4UbSrT0sEs2eXVveoxGjBvS7oXmMxQbAAAAAK3mdru1Ob9Ebx9zK0oePdQvTsPG\nDJHV6t8zy1BsAAAAANywFleLNuaXaNkJQ3Fq0dcyEjRkVLbfC81nAr7YeL1eOZ1Os2MElABaoRsA\nAABh5lJdvT5eX6r3ayLVydukf8lMVvbwLNMKzWcCvtjU19ebHQEAAAAIe9WnzuiDgr1a25SobLVo\n/rDO6pc90uxYlwV8sQEAAABgnkPlB7Wy+KiKvcmaaPfqucmd1LnHULNjfQHFBgAAAMDneDwelW7d\nrff2XdBJS5y+nGDXo5PT5UzKMTvaVVFsAAAAAEiSmhoblZ9XqverLIqQRzNvita4Ce13Uk1fotgA\nAAAAYe5CdY1W5e3S6gan0r3N+pesDhqUa/6CADeCYgMAAACEqTMnq7SioEz5LckabTX09LhUde+T\na3asVqHYAAAAAGHmxKFjemfLfm03UjTVYdFLt3RVSufA/fzM9aDYAAAAAGHiUPlBLdt+VHuVqNti\nbXp1Wm/FJyWYHcsnKDYAAABAiCsrLdeyXVU6IqdmJtn1+LQMxcTFmh3Lpyg2AAAAQAgyDEM7tu7W\nsooLOmeJ0l2dHHpyymA5ogJ/hbPWoNgAAAAAIcTj8WjrhlItO9ykFlk1u3u0xk0cJntEaL/1D+1H\nBwAAAISJFleLNuQV652TXkXLrXv7xGj42BzZbDazo/kFxQYAAAAIYk2NjVq3rkTLz0aok7dZ/5KZ\nrOzhwXUOGl+g2AAAAABBqLb6vFbl7dRHDU7187r03ZwOGjBkpNmxTEOxAQAAAILIqaMn9d6mCm1w\nB/9JNX2JYgMAAAAEgQN79uvdkuPaoyTdHG3RS1O7Bf1JNX2JYgMAAAAEsPId5Vqy47ROWGI1MzlC\nj0/ppxhnnNmxAg7FBgAAAAgwhmFob0m5luyp1llLlO7q7NCUqUPkiAzNc9D4AsUGAAAACBCGYWjX\n9j1aUn5BF+TQ7LRITZwyTBGOCLOjBTyKDQAAAGAywzBUUrhLS/bVqVF2zb4pUuMn5Yb8STV9id8U\nAAAAYBLDMLR9U6mWHmyUWxbd0zNKoycMld3O2/QbxW8MAAAA8DOPx6PCDSVaetglq7y6t3e0RozP\nkc1mMzta0KLYAAAAAH7idru1Ob9Ebx9zK0oePdg3Trljh8hqtZodLehRbAAAAIB25m5xq2B9kZad\nlOLl0tcGJGrIyEEUGh+i2AAAAADtxO12q2BdkZZUWpTqbda3BqUoa1gmhaYdUGwAAAAAH/tsUYA3\nDzYrRm49np2iQbmjzI4V0ig2AAAAgA/tLtqr/91To2ZZ9XDfGD5D4ycUGwAAAMAHDlcc0h8Lj+qk\nYnTfTQ5NmDKcVc78iGIDAAAAtMGZE1X60/q92mEkanaKXU/eMlSOSIfZscIOxQYAAABohYvna7Vs\ndbHWuZJ0a7RVr34pQ7HxcWbHClsUGwAAAOAGNDZc0srV2/V+XbxG26UXb7lJKZ2Hmh0r7FFsAAAA\ngOvganZpzZrters6UpkWj345MVVde1FoAgXFBgAAAPgnDMPQxvVF+tNxr9Lk0v8d2VnpA4ebHQv/\ngGIDAAAAXMWBPfv126JTcsuqfx2UpOzhnIsmUFFsAAAAgH9w/my1/nf1TpV4EvVAWoSmTh/J0s0B\njmIDAAAA/I2r2aUPPyrUuxfiNDnSopdn9ldcgtPsWLgOFBsAAABAUmnhLr3+SYM6Wjx6dnyqbuoz\nzOxIuAEUGwAAAIS1yiMn9Pv8/TpmxOjrfaOUO3aUrFar2bFwgyg2AAAACEsXz9dq6UfFyncl6Y54\nq75/21A5ohxmx0IrUWwAAAAQVlxNLn24+q+foxlt9+qlL3dXUmoHs2Ohja5ZbF599VWVlpYqISFB\nzz33nCSpvr5eixYt0tmzZ9WxY0fNmzdPMTExkqTly5crLy9PNptNjzzyiAYPHty+jwAAAAC4Dh6P\nRxvzivWn49JN8uiZcSnq3pfP0YSKax48OHnyZP3whz/83LYVK1Zo0KBBeuGFF5SZmanly5dLkk6c\nOKHCwkItXLhQTz75pF5//XV5vd72SQ4AAABcp51bd+t7f9ii948169uDYvWjr05R9749zY4FH7pm\nscnIyFBsbOznthUXF2vixImSpEmTJqmoqOjy9jFjxshms6ljx47q0qWLDh482A6xAQAAgGs7VH5Q\nP/5dnl4tv6S7ekbq14+M06DcLLNjoR206jM2tbW1SkxMlCQlJiaqtrZWklRTU6N+/fpdvl5ycrJq\namp8EBMAAAC4fqdPnNKf15drp5Gge1LtuvnmoXJEsjBAKPPJ4gEWi8UXdwMAAAC0ycXztVq2uljr\nXEm6NdqiV76Uodj4OLNjwQ9aVWwSExN14cKFy18TEhIk/XUPTXV19eXrnTt3TsnJyVe8j7KyMpWV\nlV2+PGfOHDmdnNUV/uFwOJg3+BUzB39i3uBvgTBzTZca9fbyfL1TE60xDqveuCtDqWmdTM2E9rN0\n6dLL32dmZiozM/P6io3X6/3cIgDDhg1Tfn6+Zs2apfz8fOXm5kqScnNz9eKLL2rGjBmqqalRVVWV\n+vTpc8X7/CzA36urq7vhBwW0htPpZN7gV8wc/Il5g7+ZOXNut1v564q0uNKm3nLp2fHddVOfv650\nxv+D0OR0OjVnzpwvbLd4r7Fs2QsvvKDy8nLV1dUpISFBc+bM0fDhw7Vw4UJVV1crNTVV8+bNu7zA\nwPLly7V+/XrZ7fYbXu65srLyBh8W0Dq86MPfmDn4E/MGfzNj5gzDUGnhbv1xX72i5NEjOakamDPQ\nrxlgjrS0tCtuv2ax8SeKDfyFF334GzMHf2Le4G/+nrkDe/brf4oqdV4OPdTboZEThspqveZivwgR\nVys2Plk8AAAAAGhvp46e1Jv5FSrzxuveTnZNv3mE7BG8ncVfMQkAAAAIaBeqa7T04x0qaEnS7bFW\nPfalgYqJi732DRFWKDYAAAAISE2NjVq5aqtW1iVoXIT0X7f1UFJqitmxEKAoNgAAAAg42zaU6PVD\nbvW2GPrlxFR17TXU7EgIcBQbAAAABIzKIyf0et5+nfJGa+7AGOWMGml2JAQJig0AAABM19TYqHc/\n2KoPLyVoVoJV/37bMDmiHGbHQhCh2AAAAMBU2zeW6vWDLvWyGFo4LU0du3HYGW4cxQYAAACmqDp6\nUr9dv08nvdH6Zka0ho4ZZXYkBDGKDQAAAPyqualJ735QqA8aEnRHvEXzv8xhZ2g7ig0AAAD8pmjz\nDr2+r1k9rIaen9ZFnTjsDD5CsQEAAEC7qzpeqTfWVei4Ea1v9I9W7lgOO4NvUWwAAADQbpqbmrTi\ng0KtbEjQTKdV3/tyjiKjosyOhRBEsQEAAEC7KNmyU7+taNRNFq8WTOuszjdx2BnaD8UGAAAAPnX6\nxCm9sbZcR41YPdo/SsPHjjY7EsIAxQYAAAA+4WpyacWHm/VefYJuj7PquzOGcNgZ/IZiAwAAgDYr\n3bJTv6loVDeLVwumdFLnHhx2Bv+i2AAAAKDVzpyo0i/X5emQEa1H+0ZqxHgOO4M5KDYAAAC4Ya5m\nl977cItW1MXrjnibnvjSYEVFR5sdC2GMYgMAAIAbsmPrLv22vEGdLYZ+Pbmj+g+aqLq6OrNjIcxR\nbAAAAHBdzp48rd+t2auD3jg92jdKI8ePMTsScBnFBgAAAP+Uq9ml91dt0fLaeN0aY9UTt2dz2BkC\nDsUGAAAAV7Vz2279pqxenSyGfjW5o9J6stoZAhPFBgAAAF9QfeqMfveXPdrvjdOjfRwaMW6UrFar\n2bGAq6LYAAAA4LKmxkZ98NF2rbjo1JdiLPr2l7MVFcthZwh8FBsAAADI7XYrf22RFp+yqbc8+uWk\nVHXtxWFnCB4UGwAAgDBmGIaKN+/UmwcuKUoefTcnVQNzRpgdC7hhFBsAAIAwZBiGSrbs1OL99WqR\nTfenR2vk+Bw+R4OgRbEBAAAII5/toXnrQL3csureXtEaNWGobDab2dGANqHYAAAAhIH/X2ga5JaF\nQoOQQ7EBAAAIYf9YaL7SO0ojx1NoEHooNgAAACHIMAxt31SqJQcb5aHQIAxQbAAAAEKIx+PR1g2l\nWnq4WRZ5dW96jEaMG0KhQcij2AAAAIQAt9utwg2lWnqkRRHy6IG+ccodO4RVzhA2KDYAAABBzO12\na1N+id4+5lGM3HokI145o7IpNAg7FBsAAIAg5G5xq2B9kZadlOLl0qMDEzV4xCAKDcIWxQYAACCI\ntLhalL++SMsqrUpRs741KEVZwzIpNAh7FBsAAIAg4Gp2af26Ir1z2q7O3mb9a3aqBuWOMjsWEDAo\nNgAAAAHM1eTS2rXb9M5Zh7rJpXlDkjQwZ6TZsYCAQ7EBAAAIQM1NTfr44yK9ey5KPdWiHwxNVf/B\nFBrgaig2AAAAAaSpoVF/WVukFTXRSleL/iO3k/oOGm52LCDgUWwAAAACQGPDJX30cZHeq41Vhlr0\nnyPTlD6QQgNcL4oNAACAiS7V1evDNcV6/2KcMuXWT0YlqVdGrtmxgKBDsQEAADBBfW2dPlxbog/q\nnBps8eipMR3Uox+FBmgtig0AAIAf1V24qA/WlGhVQ7yGWj16dlyqbuozzOxYQNCj2AAAAPjBxfO1\nWrmmRKsvJWq41dAvJnZU114UGsBXKDYAAADt6MzJKq0sKFNeS7JGWb369eRO6tJjqNmxgJBDsQEA\nAGgHxw4c0buFh1TkTdbkCGnhtC7q2DXH7FhAyKLYAAAA+NAnOz7ROztP6YAlXrfFWfXq1N6KT6LQ\nAO2NYgMAANBGLa4WFW/Zpfc+bVCNJVJ3pDr0vamZioqONjsaEDYoNgAAAK107OBRrS06pILmeKUZ\n9bq1R7TGThgmewRvsQB/438dAADADWi4WK9Nm3dp7SmPzlqjNTlaemZCJ3XrzTloADNRbAAAAK7B\nMAyV7/hE68qqtM2brCxvk+b0S1LOyGz2zgABgv+JAAAAV3Gu6ozWbynXuguRssmracl2PTy2u5JS\nWQwACDQUGwAAgL/janappHC31h6+qAprkkZb3Zo3rLP6Duonq9VqdjwAV0GxAQAA0F/PO7O2+JDy\nmxPU1ajX1K7R+t64/oqOjTE7GoDrQLEBAABh61J9gzZu3Km1pzyqtkRrcoxFP5/YSV17sRAAEGwo\nNgAAIOwcKDugj3cc02ZPsjKNZs3pl8hCAECQ438vAAAICw0X61WwcafWnPaq3uLQtHibXhrbVSmd\nWQgACAUUGwAAENIOVxzSquIj2uxJ1iBvsx4a0EGDR2TJZrOZHQ2AD1FsAABAyHE1u1S4cYc+Otak\n05YY3ey06qVx7J0BQhnFBgAAhIwzJ6q0enO51jXGq7txSTN7x2n4mCGKcESYHQ1AO6PYAACAoGYY\nhnZu26OP9p1TuSVREyO8enpcR93Uh5XNgHBCsQEAAEHJ7XZrc0GJ3jnaIkm6rXOkvjMxg/POAGGK\nYgMAAIKKq9ml/PVFeveUVfFy6aH+8Ro2erCsVqvZ0QCYiGIDAACCQlNDo9asK9aKc5HqKpcey05V\n5tCBFBoAkig2AAAgwNVduKiP1pXowzqn+qlF84d1Ur/sEWbHAhBgKDYAACAgnT9brZXrd+vjpkTl\nWgz9bEwH9ejHggAAroxiAwAAAkrV8Uqt2PCJNriTNd7u1YJJndS5x1CzYwEIcBQbAAAQEA5XfKoV\n2z9VsTdZ0yMt+q9pNym5EyfUBHB9KDYAAMA0Ho9HxVt26YMDtTpujdNtTpsendpHzsR4s6MBCDIU\nGwAA4HeX6uq1Lm+HPqy2K0Ytuv2maI2dMFiOSIfZ0QAEKYoNAADwm6qjJ/XB5grltSQp2+vSt3NS\nlDE4gyWbAbQZxQYAALQrwzBUvqNcK/ecVrmSNDVKen5qF3XqxoIAAHyHYgMAANqFq8mlDQUl+uCk\nWy7ZNKOjQ/MmZyg6NsbsaABCEMUGAAD41Pmz1Vqdv0erG5zqaTTpwf5Jyhk1SDabzexoAEIYxQYA\nAPjE/t37tWrncW33JmuszdBTYzuoe19OqAnAPyg2AACg1RobLmnjhh1afcpQvSVCNyfY9dVJvZSQ\nzPlnAPgXxQYAANywo/uP6C/Fh1TQkqQMo1n390/mcDMApqLYAACA61JbfV4FW/Yor9qiC5YoTY6x\n6PkJrG4GIDBQbAAAwFW5ml0qKdyt9Ucuaq8lWblq0UMDOih7eKbsdt5GAAgcPCMBAIDPMQxDh8oO\nav3u49rUkqhuRr0mp0Xp38b0UWz8ELPjAcAVUWwAAIAkqfrUGRVsLVfe+Qi5LDZNjrPqV+M6q0uP\nYWZHA4BrotgAABDGGhsuqXDLbhUcb9QBW6JGWz365uBUDcwZIKvVanY8ALhuFBsAAMKM2+3W3uJy\n5e0/oyJvivoblzSlR6yeHJWhqNhos+MBQKu0qdh88MEHysvLk8ViUffu3TV37lw1NTVp0aJFOnv2\nrDp27Kh58+YpJibGV3kBAEArHd1/RHmln6qg0alEb5MmpTj0yJiblJTKOWcABL9WF5uamhqtXr1a\nixYtkt1u18KFC7Vp0yadOHFCgwYN0h133KEVK1Zo+fLleuCBB3yZGQAAXKf62jpt2LhTa85IFyxR\nmhgt/WRsqnr062l2NADwqTYdPGsYhpqamuTxeORyuZScnKzi4mJNnDhRkjRp0iQVFRX5JCgAALg+\nhmHok50VeuGP6/QvKw5qd7VLDw2I128fytEjX5lCqQEQklq9xyY5OVkzZszQ3LlzFRkZqezsbGVn\nZ6u2tlaJiYmSpMTERNXW1vosLAAAuLqL52uVt3GX1lbb1CKrpifZ9fC47hxqBiAstLrYNDQ0qLi4\nWK+88opiYmL0/PPPa+PGjV+4nsViaVNAAABwdYZhaG9JudaUn1axUpSrFn17eBelZ6azqhmAsNLq\nYrNnzx517NhRcXFxkqQRI0Zo3759SkxM1IULFy5/TUhIuOLty8rKVFZWdvnynDlz5HQ6WxsHuCEO\nh4N5g18xc/C189U1+mjNdq06a5Xda+i2ztH6t2lZSkxJksPhkMvlMjsiwgjPcfC3pUuXXv4+MzNT\nmZmZrS82HTp00IEDB+RyuRQREaE9e/YoPT1dUVFRys/P16xZs5Sfn6/c3Nwr3v6zAH+vrq6utXGA\nG+J0Opk3+BUzB18wDEMVuyq0es8pFXlTNMLSqG/ndFP/7H6X987U1dUxb/A7Zg7+5HQ6NWfOnC9s\nb3Wx6dOnj0aNGqX58+fLZrOpZ8+emjZtmpqamrRw4ULl5eUpNTVV8+bNa1NwAADCXX1tnQo27tJf\nzkgtsumWlAg9OqGX4pOGmB0NAAKGxev1es0O8ZnKykqzIyBM8Jcl+BszhxtlGIYO7Nmvv+w+oa1G\nirK9Nbp1YKoG5WZd87MzzBv8jZmDP6WlpV1xe5tO0AkAAHyr7sJFFWzcpY/PWtQsu25OsuvlcZxE\nEwCuhWIDAIDJPvvszMd7TmmbN0VD5dLXs1KVNSxTNpvN7HgAEBQoNgAAmOTo/iPasONTbbwUI7vX\n0LQkux4Z312JHfjsDADcKIoNAAB+dPrEKW3cvl8bLthUZ4nU+CivfjCig3oP6M15ZwCgDSg2AAC0\ns/Nnz2nztnJtPGPopNWp0Xa3vpGdooE5AzjUDAB8hGIDAEA7aLhYr62Fe7Shskn7rUnKlUuz+yVp\n8PAsOSIdZscDgJBDsQEAwEeam5pUunWvCo7UapclRZlGo6Z0j9W/j+qn6NgYs+MBQEij2AAA0AYe\nj0d7S8r+qDggAAAgAElEQVSUv++Mthsp6ump14QuUZo7ihNoAoA/UWwAAGiFquOVWr91n/LqYxTj\nbdGk5Ag9MLKrOnThfDMAYAaKDQAA16mpoVFbtuzS+uONOmJN0PgIQ/8+MkXpA/uYHQ0Awh7FBgCA\nf8Lj8ahi1z7ll5/SFk+y+hqXdEvPOI0YM1CRUVFmxwMA/A3FBgCAf2AYhg5XfKoNu49qU2Ocor1u\nTUiwatGoNKV25VAzAAhEFBsAAP7mxKfHtbHkoDbWRarFYtO4GIv+c3QH9ejXi5NnAkCAo9gAAMJa\n9akz2rjtE22sseqcNVpjIww9PjRZ/bP7UWYAIIhQbAAAYae5qUlbN+/S+qOXdNCWqJHWFj08sKOy\nhg2UPYKXRgAIRjx7AwDCgmEYOrD3gNbtPq7NnmSlG5c0pXuMnhyboajoaLPjAQDaiGIDAAhptTUX\nlL9pt9ZW2+SSTVPibVo4Kk0duw01OxoAwIcoNgCAkOPxeLSnuExrKs6qVCnKVYu+MShZmUMHyGaz\nmR0PANAOKDYAgJBRfeqM1m0p17qLUYryunVzqkPfHNdbzqQhZkcDALQzig0AIKi1uFpUsnWX1hy6\nqAprksbYPPp+brLSM/uwqhkAhBGKDQAgKFUeOaE12/ZrfWO8OhsNmp4Wre+N76/o2BizowEATECx\nAQAEjabGRhVu2qW1x5t0zOrU5EivnhqXqu59cs2OBgAwGcUGABDwjuw7rI9LPtUGV5LSvY26rVec\nho8ZJEekw+xoAIAAQbEBAASkpoZGbdq0Ux+fdOmsJUZTYy1aMKGLOrFMMwDgCig2AICA4Xa7tbek\nXPn7zmq7N1n9jUbd1SdeuaOHyB7BSxYA4Op4lQAAmO5wxSHl7TiqjU1OJXibNCklQg+N7KaUzjlm\nRwMABAmKDQDAFGdPntaG7RUqOG9TgyVCE2Okn4xNVY9+Pc2OBgAIQhQbAIDf1NfWqbBwjwoqm3XY\nlqhRVre+kZ2igTkDZLPZzI4HAAhiFBsAQLtyNblUsnW3Co7UapclRVneJt3ay6nc0QMUGRVldjwA\nQIig2AAAfM7tdqt8xycqqDijrUaSenrqNb5zlB4b3VvOpCFmxwMAhCCKDQDAZ86ePK21W8q1tj5W\nTm+zJiTZtWhEV6V2ZYlmAED7otgAANrE7XZr57a9Wr2/Rp9YkzTObuiHo1LUe0C62dEAAGGEYgMA\naJVzVWe1dnOZ1lyMVry3Wbd0idR3xvdTTFys2dEAAGGIYgMAuG4ej0e7i/ZqdUW19liSNcbq1vwR\nyeqb2dfsaACAMEexAQBcU1Njo9atK9H7Z2yKlEe3dIrUt8f3VWw8J9AEAAQGig0A4KoaLtbr/Y+L\ntKreqX5elx4f0lkDhmTIarWaHQ0AgM+h2AAAvsDV5NJHf9mqd87HKMdq6OmxHdS9b67ZsQAAuCqK\nDQDgMrfbrYJ1RVpcaVUPufXT0cnqlUGhAQAEPooNAECGYaikcJf+d1+DouTRvCGpyhw60uxYAABc\nN4oNAIS5stJyvbnzjOoUoYf6RGv4uBw+QwMACDoUGwAIUxU7K7S49JROWaI1Jy1Sk6YNl93OywIA\nIDjxCgYAYebAnv36c/EJHVOs7uns0NRpwxThiDA7FgAAbUKxAYAw8eknh7R421Ed9Do1O9WhJ6fl\nyBHlMDsWAAA+QbEBgBB3dP8RLd5ySJ8oQXcl2/W9mwcpMirK7FgAAPgUxQYAQtSJQ8e0eNMB7fEm\naFaCXf92c6aiYqPNjgUAQLug2ABAiKk8ckJLNlSo1EjSTKdNj92coRhnnNmxAABoVxQbAAgRNaer\ntWTtLm1uSdKXY2169eZ+iktwmh0LAAC/oNgAQJBruFiv5au366PGRE2OkF6e0UsJHZLMjgUAgF9R\nbAAgSLmaXPro4616pyZGOVZDz0/tok7dhpodCwAAU1BsACDIuN1ubVhfrMUnLbpJbv10dLJ6ZeSa\nHQsAAFNRbAAgSBiGodLC3frjvnpFyaMnBndQ1rCRZscCACAgUGwAIAjs21WhP5RUqVYOPZgerZHj\nc2S1Ws2OBQBAwKDYAEAAO37wmN7cdEAHvE59pUuEpkwbIXsET90AAPwjXh0BIADVXbiot1YVaYMr\nUbMS7Jr3pUGKiubkmgAAXA3FBgACiNvt1tqPt2nxmUiNsEv/NaM3SzcDAHAdKDYAECA+2Vmh35Se\nlUOGfjSqk9IHstIZAADXi2IDACarO39Rf/ywSMUt8Xr4pghNnDqChQEAALhBFBsAMIlhGNqUX6zf\nHbNoZIT04p395UyMNzsWAABBiWIDACaoOnpSr63fp3PeSM0fmqIBQ4abHQkAgKBGsQEAP2pxtWjl\nh5u1vDZesxJsumPGCEU4IsyOBQBA0KPYAICfVOys0Cul1UqSV7+e2kldegw1OxIAACGDYgMA7ay+\ntk6/WZynLc1x+mqPCI2fPIbFAQAA8DGKDQC0E8MwVFhQotePSCMipZfu6C9nEosDAADQHig2ANAO\nTp84pf9e84nOeCP1vSFJGjVxiurq6syOBQBAyKLYAIAPuVvcev/DzXrnglMz4y369xnD5Yh0mB0L\nAICQR7EBAB/Zt3ufXi0+I6cM/WpKR6X1ZHEAAAD8hWIDAG3UcLFef3p/mza3JOiRmyI0cepYFgcA\nAMDPKDYA0EqGYWjbhlL99rChHLv04sw+SkhONDsWAABhiWIDAK1w5kSVfrOmXJXeKH1ncLyyho0w\nOxIAAGGNYgMAN8Dd4taHq7bo7fNxmuG06AczcuWIYnEAAADMRrEBgOt0oOyAXtlWpRgZ+sWkVHXr\nzeIAAAAECooNAFzDpbp6/fn9bdroStDDXe2aPJ3FAQAACDQUGwD4J7ZuKNZvDxnKtkkvzkhXQock\nsyMBAIAroNgAwBWcPXlav12zV8eNGD0xKE7Zw1kcAACAQEaxAYC/4/F4tPqjQr1VHa1bY6367u05\nioyKMjsWAAC4BooNAPxN1fFK/deafWqSTc9M6KDufYaZHQkAAFwnig2AsGcYhv7y0Rb9+Wy0ZiVZ\ndMftY2S38/QIAEAw4ZUbQFg7c6JKL338iRpl0zPjU9S9L3tpAAAIRhQbAGHJMAyt/Uuh/vd0lGYm\nWHTnzDGyR/CUCABAsOJVHEDYqT51Ri+v3qvz3gg9NTZZPfuzlwYAgGBHsQEQNjwej9Z+vFVvno7S\nrXEW3XPHKEU4IsyOBQAAfIBiAyAsHD94TK9s+FQtsuqnoxPVewB7aQAACCUUGwAhzdXk0rL3N2tV\nQ4Lu7WDTrbeNZsUzAABCEK/uAELW3pK9emXXRaVZvHr+5jR17NrZ7EgAAKCdtKnYXLp0Sa+99pqO\nHz8ui8Wib33rW+rSpYsWLVqks2fPqmPHjpo3b55iYmJ8lRcArqnu/EX94YPtKnEn6Bu9IzRqwihZ\nrVazYwEAgHbUpmLz+9//Xjk5OfrOd74jj8ej5uZmvfvuuxo0aJDuuOMOrVixQsuXL9cDDzzgq7wA\n8E9t21iq1w56NMJh0Ut3ZSguwWl2JAAA4Aet/hPmpUuXVFFRocmTJ0uSbDabYmJiVFxcrIkTJ0qS\nJk2apKKiIt8kBYB/or62Ti/8cZ3eOOjSd7Oj9a0HplJqAAAII63eY3PmzBk5nU698sorOnr0qHr3\n7q1HHnlEtbW1SkxMlCQlJiaqtrbWZ2EB4Ep2bN2l//qkSbkOadFdgxQTF2t2JAAA4GetLjaGYejw\n4cP6+te/rvT0dP3P//yPVqxY8YXrWSyWNgUEgKu5VN+gP7y3VcUupx7LiNLQ0SPNjgQAAEzS6mKT\nnJyslJQUpaenS5JGjRqlFStWKDExURcuXLj8NSEh4Yq3LysrU1lZ2eXLc+bMkdPJYSPwD4fDwbwF\nudLCHfp1UbUGRVj1+kO5ik+68nNNoGDm4E/MG/yNmYO/LV269PL3mZmZyszMbH2xSUxMVEpKiior\nK5WWlqY9e/aoW7du6tatm/Lz8zVr1izl5+crNzf3irf/LMDfq6ura20c4IY4nU7mLUg1NTbqzeVb\ntMkVr2+lR2jkhFGSAv/5g5mDPzFv8DdmDv7kdDo1Z86cL2xv06poX/3qV/XSSy/J7XarU6dOmjt3\nrgzD0MKFC5WXl6fU1FTNmzevLT8CAC7bt3ufXig+p1426YUZfZTQIcnsSAAAIEBYvF6v1+wQn6ms\nrDQ7AsIEf1kKLq4ml956b5PWNsbr0e6GJkwZYXakG8bMwZ+YN/gbMwd/SktLu+L2Nu2xAYD29ukn\nh7So8JQ6WqRFt/VUcqcOZkcCAAABiGIDICB5PB6tWLlJK2rj9H+62jRl+iRZra0+9RYAAAhxFBsA\nAafmdLUWrdqjJq9Vz92cpk7dupgdCQAABDiKDYCAUrJlp17a36JpsV595c6xskfwNAUAAK6NdwwA\nAoKr2aU3392oTU1OfScrRtnDh5sdCQAABBGKDQDTVR45oefWH1ayxaKFt/dlGWcAAHDDKDYATFWw\ndpteP2HXnA7Sl2ewQAAAAGgdig0AUzQ2XNJv392i8pZY/WRUgtIHDjM7EgAACGIUGwB+d3T/Ef1q\nU6X62C16/p5BinHGmR0JAAAEOYoNAL/auH67fnPMpv+TZtG0L001Ow4AAAgRFBsAfuFuceuPyzao\nsClWPxmZoPSBfcyOBAAAQgjFBkC7u1Bdo+dW7pJNFj03q78SkhPNjgQAAEIMxQZAuzqwZ79+UXxe\nE2O9uv/uCbLbedoBAAC+xzsMAO3m49Vb9L9VkfpWL6vGTJpidhwAABDCKDYAfM7d4tbvlhZoZ3O0\nnp2Qqpv6dDc7EgAACHEUGwA+dam+Qc+9vU1uWfWruwcpLsFpdiQAABAGKDYAfKb61Bk9/dE+pUd4\n9c17xyvCEWF2JAAAECYoNgB84lD5QT2zrUZfTvLozjsmy2q1mh0JAACEEYoNgDbbvrFULx3y6ps9\nLRo7eZLZcQAAQBii2ABokw/e36i3z0Xrh0PilDEkw+w4AAAgTFFsALSK2+3W75cWaEdTtH45tYs6\n9+hqdiQAABDGKDYAbtil+gYteHubmr1W/fLOLDmT4s2OBAAAwhzFBsANOVd1Rk+v2qeeEYa+NWeC\nHJEOsyMBAABQbABcv8MVh/R04TndkuDW7DunsPIZAAAIGBQbANelZMtOLdrv0Td6SBOmTDY7DgAA\nwOdQbABc06oPNmpJdbT+Y3CsBuQMMDsOAADAF1BsAFyV2+3WH94uUEljtH4+pbPSenYzOxIAAMAV\nUWwAXFFjwyU9//ZWNRg2/eLOTMUnJZgdCQAA4KooNgC+4PzZc3r6/TJ1sxv6/n3jWPkMAAAEPIoN\ngM+pOnpSP1l3TOPi3Lp/NiufAQCA4ECxAXDZp58c0lNba3RXB7dunznF7DgAAADXjWIDQJK0p3iv\nfrW3Wd/o4dWEKePNjgMAAHBDKDYAVJhfrFeOWPTdjAgNGTXM7DgAAAA3jGIDhLnVqzbrrTMO/Wh4\ngvpm9jU7DgAAQKtQbIAwtnxFgT46H6lnJ3fhHDUAACCoUWyAMLXs3QKtveDQ07f0VMeunc2OAwAA\n0CYUGyAMLVmWp4KLkXr61nR16NLR7DgAAABtRrEBwohhGHrrnQJtrovUU7f1VUrnVLMjAQAA+ATF\nBggThmHozbfzVdQQqadv76ek1A5mRwIAAPAZig0QBgzD0B+X5qv0kkNPzRygxA7JZkcCAADwKYoN\nEOIMw9Dv3srT3iaHnpqVpYTkRLMjAQAA+BzFBghhhmHo9cV52tccoafuHCxnUrzZkQAAANoFxQYI\nUR6PR/+9OE+HXQ799O4cxSU4zY4EAADQbig2QAjyeDx65c/5OtESoZ/MHqrY+DizIwEAALQrig0Q\nYtxut17+c4Gq3Db9eE6uYuJizY4EAADQ7ig2QAgxDEOvLc7XabddP7p3hKJjY8yOBAAA4BdWswMA\n8A3DMPT7JXk60uLQf96TS6kBAABhhT02QIh4650C7W6M1NN3ZivGyWdqAABAeKHYACHg3eUF2lQX\nqWdu78+SzgAAICxRbIAgt+rDTfroQpSevaWnklJTzI4DAABgCooNEMTWf7xVy8449PTkNKV27WR2\nHAAAANNQbIAgtSW/SH+otOup8alK69nN7DgAAACmYlU0IAiVbNmp145Y9KMRierep4fZcQAAAExH\nsQGCzJ7ivVq036Mnh8QpfWAfs+MAAAAEBIoNEET27d6nX+1t1vcGRmjAkAyz4wAAAAQMig0QJA5X\nHNKzpXV6PN2iwSOyzY4DAAAQUCg2QBA48elx/aywRo92NzRi/FCz4wAAAAQcig0Q4KqOV+rHBad0\nf5cWjZ8ywuw4AAAAAYliAwSwc1Vn9eM1xzQrpVnTvzTG7DgAAAABi2IDBKja6vP68ar9mpbQpNtn\njjc7DgAAQECj2AABqL62Tj9ZuVcjYpt1z92TzI4DAAAQ8OxmBwDweY0Nl/TUO6UaEOXWg/dMNjsO\nAABAUKDYAAGkualJzy7ZprQIj75+72RZrexUBQAAuB4UGyBAtLha9KvFWxRv8+qx+ybJZrOZHQkA\nACBo8OdgIAC43W4tWrxBkvRv902Q3c7fHAAAAG4ExQYwmcfj0auLC1TrseoH941RhCPC7EgAAABB\nh2IDmMgwDL3xVr6Ot9j1H3NGKDIqyuxIAAAAQYnjXQCTGIah/1mSp4rmCP3s7hzFxMWaHQkAACBo\nUWwAExiGoTffztfuRod+dme24hKcZkcCAAAIahQbwARL3ilQUUOknp45UPFJCWbHAQAACHoUG8DP\n3n4nXxvrIvXM7f2U0CHJ7DgAAAAhgWID+NGK9zZoXW2knrktXUmpHcyOAwAAEDIoNoCfrPpgoz6s\nidQzt/RQSueOZscBAAAIKRQbwA/Wri7UO2cj9fSUrurYtbPZcQAAAEIOxQZoZwXrtunNKruentBJ\nXXp0NTsOAABASOIEnUA7Kswv1u+O2/TjUSnqlt7d7DgAAAAhi2IDtJPizTv1yhGL/u/wRPXK6G12\nHAAAgJBGsQHawa7tu/XCAY/+Y0is+mT2MTsOAABAyKPYAD5WvqNcz5W36AdZkRowJMPsOAAAAGGB\nYgP40P7d+/XzXY2a19+uQblZZscBAAAIGxQbwEcOVxzS06UX9a/pFg0dPdjsOAAAAGGFYgP4QNXR\nk/pZYY2+0d3QyPFDzY4DAAAQdig2QBudP3tOP1l3TLNTmzV+ygiz4wAAAIQlig3QBpfq6vXU+2Ua\nF9eoL88YZ3YcAACAsGVv6x0YhqEnn3xSycnJmj9/vurr67Vo0SKdPXtWHTt21Lx58xQTE+OLrEBA\ncTW79ItlReod4dH9s6eYHQcAACCstXmPzapVq9S1a9fLl1esWKFBgwbphRdeUGZmppYvX97WHwEE\nHI/Hoxff2qhIi1ff/MokWa3s/AQAADBTm96NnTt3Tjt27NDUqVMvbysuLtbEiRMlSZMmTVJRUVHb\nEgIBxjAM/X5Jvqo9dn333tGyR7R5xycAAADaqE3F5g9/+IMeeughWSyWy9tqa2uVmJgoSUpMTFRt\nbW3bEgIBZvl7G7SrKUo/vDNHUdHRZscBAACA2lBsSktLlZCQoJ49e8rr9V71en9feoBgl792m1Zd\niNaPb+0rZ1K82XEAAADwN60+hqaiokLFxcXasWOHXC6XGhsb9dJLLykxMVEXLly4/DUhIeGKty8r\nK1NZWdnly3PmzJHT6WxtHOCGOByOG563os0leuOEXc9N66pe/dLbKRlCVWtmDmgt5g3+xszB35Yu\nXXr5+8zMTGVmZsri/We7W65TeXm53n//fc2fP19vvvmm4uLiNGvWLK1YsUINDQ164IEHrut+Kisr\n2xoFuC5Op1N1dXXXff0j+w7rR4Xn9L3MSGUPH9SOyRCqbnTmgLZg3uBvzBz8KS0t7Yrbfb6U06xZ\ns7Rnzx498cQT2rt3r2bNmuXrHwH41bmqM3pqyxl9rZuHUgMAABCgfLLHxlfYYwN/ud6/LDVcrNd/\nLNupcYlu3XPXpPYPhpDFXzPhT8wb/I2Zgz9dbY8N69QCV9HiatEv3ylSRqShu2dNNjsOAAAA/gnO\nKghcgWEYevmtDYqQV9+YM5ETcAIAAAQ43q0BV7B4WYGOuyP0vTmjOAEnAABAEOAdG/APPv5oswoa\novXL2/soOjbG7DgAAAC4DhQb4O+UbNmpP52O1DMTOioptYPZcQAAAHCdOBQN+JuDZQe1aL9H/z44\nVt3Su5sdBwAAADeAPTaApKrjlXpm+3l9q5c0IGeA2XEAAABwg9hjg7BXd/6inlpzRLNSGjVm0nCz\n4wAAAKAV2GODsOZqcunZ5aXKifLojplTzY4DAACAVqLYIGx5PB69sGSjEqzSV++dZHYcAAAAtAGH\noiFs/fHtfFV77Jr3lbGy2WxmxwEAAEAbsMcGYentpR+r6FK0fj5roCKjosyOAwAAgDai2CDsbFy/\nXUtO2/XzqWlKSE40Ow4AAAB8gGKDsLJz62795phNv5rUWZ17dDU7DgAAAHyEz9ggbBwoO6AFFS36\nQVak+g3KMDsOAAAAfIhig7Bw8vBxPVN0QXN7ejUoN8vsOAAAAPAxig1C3rmqs/pp/kl9paNLoyfl\nmh0HAAAA7YBig5BWX1unn63ap6nORn3ptrFmxwEAAEA7YfEAhKzmpiY9+06JBkR6dM9dk82OAwAA\ngHZEsUFIcrvdev6tzUq0So/eO1FWKzsnAQAAQhnv9hByDMPQf79VoEuGVf/2lfGy2+nvAAAAoY53\nfAg5i5cV6KDLoWfuyZEjymF2HAAAAPgBxQYh5cMPNmljfZR+fns/xTjjzI4DAAAAP6HYIGRsytuu\nZWcj9ezUNCWlppgd5/+1d+9RVZf5Hsc/mw2biyAIihcYREW8oGhOmmYoaOc4h5qyKW+1Gj2nsWmc\nM1PUOnmczjEvoTl56eKxaSbHbmtlWso6NY7WSfBKJmomGnlHRQ0BxQtsYF/OHzWs6TKlyN7Pvrxf\n/8je7A0f1vr+8fvs5/k9AgAAwIu4xwYBYe/Hn+qlcqv+e2i8OndNMh0HAAAAXkaxgd87tP+QFh5o\n0n9khKt7nx6m4wAAAMAAig38WsWxk8rfeUHTUt3KHNzPdBwAAAAYQrGB36o+W6nZRRWa1LFRw7Jv\nNB0HAAAABlFs4JcuXbio2esO6ta29RrzL8NNxwEAAIBhnIoGv2Ovr9e8NbvVL8Kpe+7KMR0HAAAA\nPoBiA7/iaHJo4cpiJVjdemBCtkJCWHQEAAAAW9HgR1wul5atLFKTLPrtxCxZrVbTkQAAAOAjKDbw\nG6+vLtKJJpumjxsiW7jNdBwAAAD4ELaiwS+sWbtJO65EaP4dfRQV3cZ0HAAAAPgYig183rq/bNVf\nL0Ro3phUxbZvZzoOAAAAfBDFBj5t4wcf6e1Km57K6aIOSR1NxwEAAICPotjAZxUXlejVilDNvaW9\nuqQmm44DAAAAH8bhAfBJu7d/omXHLfrvwXFK6ZlqOg4AAAB8HMUGPqd0134tOejUjAGRSstIMx0H\nAAAAfoBiA59yaN9BLdhn16O9Q9X3hr6m4wAAAMBPcI8NfMbxz4/pqV0X9evu0g1DB5iOAwAAAD/C\nig18wunjpzR7e5X+NdmpoSNuNB0HAAAAfoZiA+MqK85qZmGFJnRsUPatN5mOAwAAAD9EsYFR589V\n6ckNx3RbO7t+knuL6TgAAADwU9xjA2Nqay7oyXfLlBXTpLvG5piOAwAAAD9GsYERVy5e1uyCTzUw\n0qGJd2ebjgMAAAA/x1Y0eF39lTrNfXuXeoY3acqEbIWEMIYAAAC4PqzYwKsa7Y2av2qHOoY69eDE\nHEoNAAAAWgXFBl7T1Nik37+5VVEWt34zKVtWq9V0JAAAAAQIPi6HVzgcDj375ma5JD06KUuhYXRq\nAAAAtB6KDTzO6XRq2ZtFqnWG6PGJw2QLt5mOBAAAgABDsYFHuVwuLV9ZpFNNYfrd+CGKiIw0HQkA\nAAABiP1A8BiXy6XXVxfpswab5vxsgKKi25iOBAAAgABFsYHHrF6zSTuvhCv/jr6KiWtrOg4AAAAC\nGMUGHrG2YJMKL0Zo3m09Fdu+nek4AAAACHAUG7S6de9t0brzEcof01XxHdubjgMAAIAgQLFBq/q/\nDcV6+1y48kclKTGpk+k4AAAACBIUG7SazRs/1htnQjU3q6M6d00yHQcAAABBhOOe0Sp2bN6l5SdC\n9OTQeP0oLcV0HAAAAAQZig2u2+7ivVp6VHpiUFt1693DdBwAAAAEIYoNrkvprlIt+dyh/8yMVHpm\nuuk4AAAACFIUG7TY53vLtGBfgx7rE6aMQX1NxwEAAEAQo9igRY4cOKz8PZf12zSLBt6UaToOAAAA\nghzFBtfsxKHjmvPxBf2yq1uDbxlkOg4AAABAscG1OX38lJ7cek6TkxwanjPYdBwAAABAEsUG1+CL\nU2c0s7BCEzo2aNQ/DTUdBwAAAGhGscFVqT5bqZnvl+unCQ36Se4tpuMAAAAAX0OxwQ+6UFWjmesO\n69ZYu+68Y4TpOAAAAMC3UGzwvS6dv6gn//eAhkXbNe7ubNNxAAAAgO8UajoAfNeVi5c1e+0n6h/R\npHvvyTEdBwAAAPiHKDb4Tvb6euW/XaJuNqf+bWKOQkJY3AMAAIDvotjgW5oam/TMymLFW916aFI2\npQYAAAA+jytWfI3T6dTzb26WW9LDk0bIarWajgQAAAD8IIoNmrlcLv1xZZGqnFY9PnGYwmxhpiMB\nAAAAV4WtaGj2xuoiHWqwae49NygiMtJ0HAAAAOCqUWwgSXpn7SbtuBKheXf0UZu20abjAAAAANeE\nYgOtX7dN6y9EaN6YVMW2b2c6DgAAAHDNKDZBbsvGj/XWFzbl53RWh6SOpuMAAAAALcLhAUGsZNsn\n+tOJEM0cFq8uqcmm4wAAAAAtRrEJUgf2HNBzh5yaMTBa3Xr3MB0HAAAAuC4UmyB0rOyInt5bp7xe\noYsZTogAABN1SURBVOozsLfpOAAAAMB1o9gEmTPlFZpTXKOpKS4NGjbAdBwAAACgVVBsgkjNF1Wa\n9eFJjUtsUNaoIabjAAAAAK2GYhMkLl24qFl/KdOoWLtyb7vFdBwAAACgVbX4uOfq6motXbpUtbW1\nslgsGj16tHJzc3X58mU9++yzOnfunBITE5WXl6eoqKjWzIxrZL9Sr/w1e5QZ4dC4u3JMxwEAAABa\nXYuLjdVq1eTJk5Wamiq73a7p06drwIABKiwsVP/+/XXnnXeqoKBAa9eu1X333deamXENmhqb9PtV\nxeoY6ta/TshWSAiLdAAAAAg8Lb7KjYuLU2pqqiQpIiJCSUlJqq6uVklJiUaOHClJys7O1s6dO1sl\nKK6d0+nUCys3S5L+feIIWa1Ww4kAAAAAz2iVj+8rKytVXl6u9PR01dbWKi4uTtKX5ae2trY1fgWu\nkcvl0vK3ilTpCNXj44cpzBZmOhIAAADgMS3eivY3drtdixcv1pQpUxQREfGt71ssluv9FWiB1Ws2\nqdQervyfZSqiTaTpOAAAAIBHXVexcTqdWrRokUaMGKHBgwdL+nKV5sKFC83/xsbGfud79+/fr/37\n9zc/Hj9+vGJiYq4nDr7yztsfaOOlSD0/foDad040Hccn2Ww25g1exczBm5g3eBszB29btWpV89cZ\nGRnKyMiQxe12u1v6A5cuXaqYmBhNnjy5+bk33nhD0dHRGjt2rAoKCnTlypWrPjzg9OnTLY2Cr2wt\n/FjLj1s0b1QXde6aZDqOz4qJidGlS5dMx0AQYebgTcwbvI2Zgzd16dLlO59v8YpNWVmZtmzZopSU\nFD3++OOyWCyaNGmSxo4dqyVLlqiwsFAdOnRQXl5ei0Pj2ny6s1QvlVs1e2gcpQYAAABB5bpWbFob\nKzYtd6zsiGZ+dEH/kWFT5uD+puP4PD5Zgrcxc/Am5g3exszBm1p9xQa+44tTZzS3uFpTu4pSAwAA\ngKDE/9bo52prLmj2+8d0Z4JdI0YNMR0HAAAAMIIVGz9mv1Kv/IK9GhLl1J13jDIdBwAAADCGYuOn\nHE0OLVxVrE6hbt0/Ltt0HAAAAMAotqL5IZfLpT+sLFKT26J/n5Alq9VqOhIAAABgFMXGD618Z5OO\nNoVr+rjBsoXbTMcBAAAAjGMrmp9Zv26rNl2K1NM/TVNUTLTpOAAAAIBPoNj4kY82l+jNygjNH9VJ\n7Tq0Nx0HAAAA8BlsRfMTB/Yc0P8cteiJQW3VJTXZdBwAAADAp1Bs/MCxsqN6em+d8tKtSs9MNx0H\nAAAA8DkUGx93trxCc4qr9UCyU4NuHmg6DgAAAOCTKDY+7Py5Ks368ITubm/XyFtvMh0HAAAA8FkU\nGx91ufaSZr/7mbKi7br9p1mm4wAAAAA+jVPRfFCD3a757+xSr3CnJt2TYzoOAAAA4PMoNj7G0eTQ\nopXbFBsiTZ0wUiEhLKoBAAAAP4SrZh/icrn04spNqneF6JGJWQoNpXcCAAAAV4Ni40NeX12k401h\nmjF+iGwRNtNxAAAAAL/BkoCPWFuwSTuuRGj+HX0UFd3GdBwAAADAr1BsfMCHG4r13vkIzf/nropt\n3850HAAAAMDvUGwM+3jLbr12JkxzszooMbmT6TgAAACAX+IeG4P27z6gF4649bsbopWS1tV0HAAA\nAMBvUWwMOVZ2RE9/Wq+8XqHqNaC36TgAAACAX6PYGHCmvEJzims0NcWpQcMGmI4DAAAA+D2KjZed\nP1elWR+e1D0dGjRi1BDTcQAAAICAQLHxosu1lzT73c+UHVOv226/xXQcAAAAIGBwKpqXNNjtmvfO\nLvUOd2rC3Tmm4wAAAAABhWLjBY4mhxau3K52VremTsxWSAgLZQAAAEBr4grbw1wul5atLFKjy6KH\nJ2bJarWajgQAAAAEHIqNh722qkgnmmyaPn6IbOE203EAAACAgMRWNA9aW7BJO+siNG9sX0VFtzEd\nBwAAAAhYFBsP+b8NxXrvfITmj+mq2Pg403EAAACAgEax8YAdW3br9TNheiorUYlJnUzHAQAAAAIe\n99i0stJdpVp6xK3fDYrRj9JSTMcBAAAAggLFphUd/eyIFuxr0GO9w9Qrs5fpOAAAAEDQoNi0kjPl\nFZr7UY0eTHFq4NBM03EAAACAoEKxaQU1X1Rp1ocnNa6DXVmjhpiOAwAAAAQdis11unThomb/pUzZ\nbeuVe3uW6TgAAABAUOJUtOtgr6/XvDW71TfcqQk/yzEdBwAAAAhaFJsWcjQ5tHBlsRKsbv1iYrZC\nQlj8AgAAAEzharwFXC6X/mdlkZpk0W8nZslqtZqOBAAAAAQ1ik0LvLqqUCebbPrP8TfJFm4zHQcA\nAAAIemxFu0bvrN2kkrpIzRvbV5FtokzHAQAAACCKzTX5YP12/fVChOaPSVVsfJzpOAAAAAC+QrG5\nSjs279IbZ23KH5GoDkkdTccBAAAA8He4x+YqlO4q1dKj0n/9uK2Se6SYjgMAAADgGyg2P+BY2REt\n2Negx3qHqWf/dNNxAAAAAHwHis33OFteoTnFNZqa4tTAoZmm4wAAAAD4Byg2/8D5c9Wa9eEJ3d3e\nrhGjhpiOAwAAAOB7UGy+Q92ly5r77gHdEl2v23+aZToOAAAAgB/AqWjf0Ghv1PzVO9XD5tK99+SY\njgMAAADgKlBs/o7T6dRzb21RZIj0y4kjFRLCghYAAADgD7hy/4rL5dLylUWqcVr12MThCg2l8wEA\nAAD+gqv3r6xeu1mlDeHK/1mmwiMiTMcBAAAAcA0oNpLWr9umDy9GaH5uD8XEtTUdBwAAAMA1Cvpi\nU1xUopWVNuVnd1JCpw6m4wAAAABogaC+x6Z0V6mWHbfoicFxSur2I9NxAAAAALRQ0BabY2VH9Pt9\ndj3WJ0w9M3qajgMAAADgOgRlsTlTXqE5xTWamuLSwJsyTccBAAAAcJ2CrthUnz2nJzee1LjEBmWN\nGmI6DgAAAIBWEFTF5uL5Ws1ad1D/1Nau3NtuMR0HAAAAQCsJmlPR6i5f0dy1e3VDpEN335VtOg4A\nAACAVhQUxabR3qinV3+slDCXpkzIUUhIUC1UAQAAAAEv4IuNw+HQopVbFGmRfjVpJKUGAAAACEAB\nfZXvcrn04pubVOcK0WMTsxQaGvA9DgAAAAhKAVtsXC6XXnmrSOVNYZoxfohsETbTkQAAAAB4SMAu\nYbxTsFl76sOVf1c/RUW3MR0HAAAAgAcFZLH561+26oPaCM3P7a627WJNxwEAAADgYQFXbDZv/Fir\nKsOVn9NZCZ0STccBAAAA4AUBdY/Nru2f6OUTIZo5LEFdUpNNxwEAAADgJQFTbA7sOaBnDzo1Y2C0\nuvXubjoOAAAAAC8KiGJz9LMjenpvnfJ6harPwN6m4wAAAADwMr+/x6bi2EnN2VGjX3aVBg0bYDoO\nAAAAAAP8esXmXMUXmlV0WpM6Nmp4zmDTcQAAAAAY4rfFprbqvGZtOKKfxNk15l+Gm44DAAAAwCC/\n3IpWd+myZv/vPg1p49Ddd40yHQcAAACAYX5XbBrsduWvLlEPm1P3j8sxHQcAAACAD/CrYuNocmjh\nyu2KDXHrlxNHKiTEb3fSAQAAAGhFftMMnE6nlr65SU1u6ZFJWQoN9atOBgAAAMCD/KLYuFwu/fmt\nIp12hGr6+KGyhdtMRwIAAADgQ/xi2WPVmk3aZw9X/l2ZimwTZToOAAAAAB/j88XmvXe3qPBipObf\nlqaYdm1NxwEAAADgg3y62BR+8JHWVIVr3ugkxXdsbzoOAAAAAB/ls8Vmx5bdeqUiVHOGt1enrkmm\n4wAAAADwYT55eMCnO0u19IhbvxsUo67pqabjAAAAAPBxHlux+eSTT/TKK6/I7XYrJydHY8eOvar3\nHdp/SM/sb9BjfcLUK7OXp+IBAAAACCAeWbFxuVxavny5nnjiCS1atEjbtm1TRUXFD77v5OETempn\nraalujXwpkxPRAMAAAAQgDxSbA4fPqzOnTurQ4cOCg0N1fDhw7Vz584ffN+sLWd0f6cmDcu+0ROx\nAAAAAAQojxSbmpoaJSQkND+Oj49XTU3ND77vjvgG3fqTYZ6IBAAAACCA+dThAXfeOcJ0BAAAAAB+\nyCOHB8THx6uqqqr5cU1NjeLj47/2mv3792v//v3Nj8ePH68uXbp4Ig7wnWJiYkxHQJBh5uBNzBu8\njZmDN61atar564yMDGVkZHim2KSlpens2bM6d+6c2rVrp23btunhhx/+2mv+FgAwYdWqVRo/frzp\nGAgizBy8iXmDtzFz8LbvmjePFJuQkBA98MADeuqpp+R2uzVq1CglJyd74lcBAAAAgOf+H5uBAwfq\nueee89SPBwAAAIBmPnV4AOAtbIOEtzFz8CbmDd7GzMEXWNxut9t0CAAAAAC4HqzYAAAAAPB7FBsA\nAAAAfs9jhwcAplRXV2vp0qWqra2VxWLR6NGjlZubq8uXL+vZZ5/VuXPnlJiYqLy8PEVFRUmS1q5d\nq8LCQlmtVk2ZMkUDBgww/FfAH7lcLs2YMUPx8fGaPn06MwePqaur0x/+8AedPHlSFotFv/rVr9S5\nc2fmDR7z3nvvqbCwUBaLRSkpKZo2bZrsdjszB5/Cig0CjtVq1eTJk7V48WLl5+drw4YNqqioUEFB\ngfr376/nnntOGRkZWrt2rSTp1KlTKi4u1pIlSzRjxgy9/PLL4tYztMS6deuUlJTU/JiZg6esWLFC\nN9xwg5YsWaJnnnlGSUlJzBs8pqamRuvXr9eCBQu0cOFCOZ1Obd26lZmDz6HYIODExcUpNTVVkhQR\nEaGkpCRVV1erpKREI0eOlCRlZ2dr586dkqSSkhLdfPPNslqtSkxMVOfOnXX48GFT8eGnqqurtWfP\nHo0ePbr5OWYOnlBXV6eysjLl5ORI+vLDnKioKOYNHuVyuWS32+V0OtXY2Kj4+HhmDj6HrWgIaJWV\nlSovL1d6erpqa2sVFxcn6cvyU1tbK+nLT6LS09Ob3xMfH6+amhojeeG/Xn31Vd1///2qq6trfo6Z\ngydUVlYqJiZGy5YtU3l5ubp3764pU6Ywb/CY+Ph43X777Zo2bZrCw8OVmZmpzMxMZg4+hxUbBCy7\n3a7FixdrypQpioiI+Nb3LRaLgVQIRLt371ZsbKxSU1O/d7sFM4fW4HK5dOzYMY0ZM0YLFixQeHi4\nCgoKvvU65g2t5cqVKyopKdGyZcv00ksvqaGhQVu2bPnW65g5mMaKDQKS0+nUokWLNGLECA0ePFjS\nl58mXbhwofnf2NhYSV9+klRVVdX83urqasXHxxvJDf9UVlamkpIS7dmzR42Njaqvr9cLL7zAzMEj\n4uPjlZCQoB49ekiShg4dqoKCAuYNHrNv3z4lJiYqOjpakjRkyBB9/vnnzBx8Dis2CEgvvviikpOT\nlZub2/zcj3/8YxUVFUmSioqKdOONN0qSbrzxRm3fvl0Oh0OVlZU6e/as0tLSTMSGn7r33nv14osv\naunSpXrkkUfUr18//eY3v2Hm4BFxcXFKSEjQ6dOnJX150ZmcnMy8wWPat2+vQ4cOqbGxUW63m5mD\nz7K4OaYCAaasrExPPvmkUlJSZLFYZLFYNGnSJKWlpWnJkiWqqqpShw4dlJeXpzZt2kj68ljKjRs3\nKjQ0lGMpcV0OHDigd999t/m4Z2YOnnD8+HG99NJLcjgc6tixo6ZNmyaXy8W8wWNWr16t7du3y2q1\nKjU1VQ899JDsdjszB59CsQEAAADg99iKBgAAAMDvUWwAAAAA+D2KDQAAAAC/R7EBAAAA4PcoNgAA\nAAD8HsUGAAAAgN+j2AAAvOrXv/61SktLTccAAAQYig0AwO8UFRVp5syZpmMAAHwIxQYA4JcsFovp\nCAAAHxJqOgAAIDgdPnxYr7zyiioqKmSz2XTTTTdp8uTJslqtkqS9e/dqxYoVqq2t1fDhw3Xq1CmN\nGDFCvXr10p/+9Ce5XC79/Oc/l9Vq1YoVKwz/NQAA01ixAQAYYbVaNWXKFP35z39Wfn6+SktLtWHD\nBknSpUuXtHjxYt13331avny5unTpooMHD0qSkpKSNHXqVKWnp+u1116j1AAAJFFsAACGdOvWTWlp\nabJYLGrfvr1uvfVWHThwQJK0Z88epaSkaPDgwQoJCVFubq5iY2MNJwYA+DK2ogEAjDhz5oxeffVV\nHT16VI2NjXI6nerevbskqaamRgkJCV97/TcfAwDw9yg2AAAjXn75ZXXr1k15eXkKDw/XunXrtGPH\nDklSu3bttGvXrq+9vrq6uvlrDg4AAHwTW9EAAEbY7XZFRkYqPDxcFRUVev/995u/N2jQIJ04cUIl\nJSVyuVxav369amtrm78fGxur6upqORwOE9EBAD6IYgMA8Kq/rbbcf//92rp1qyZPnqw//vGPuvnm\nm5tfExMTo0cffVSvv/66HnjgAVVUVKhHjx4KCwuTJPXr10/Jycl68MEH9Ytf/MLI3wEA8C0Wt9vt\nNh0CAIDv43a79dBDD+nhhx9W3759TccBAPggVmwAAD5p7969qqurU1NTk9asWSNJ6tmzp+FUAABf\nxeEBAACfdPDgQT3//PNyOBxKTk7W448/3rwVDQCAb2IrGgAAAAC/x1Y0AAAAAH6PYgMAAADA71Fs\nAAAAAPg9ig0AAAAAv0exAQAAAOD3KDYAAAAA/N7/A5bR7YTDN9sdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f49cc70f898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compate MSD with plots\n",
    "fig, ax = plt.subplots(figsize=(14, 10))\n",
    "\n",
    "msd1.plot(x='lagt', y='msd', ax=ax)\n",
    "msd2.plot(x='lagt', y='msd', ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 1.55 s per loop\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.554775274998974"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t1 = %timeit -o tp.msd(traj, mpp=mpp, fps=fps, max_lagtime=max_lagtime, pos_columns=pos_columns)\n",
    "t1.best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 5.4 ms per loop\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.005402266660094028"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t2 = %timeit -o msd_fft(traj, mpp=mpp, fps=fps, max_lagtime=max_lagtime, pos_columns=pos_columns)\n",
    "t2.best"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}