fig_6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## S conditioned on H distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "pylab.rcParams['figure.figsize'] = (10, 6)\n",
    "\n",
    "#%matplotlib inline\n",
    "\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy.random as rand\n",
    "from sklearn.kernel_ridge import KernelRidge\n",
    "from sklearn import linear_model\n",
    "from pandas.plotting import autocorrelation_plot as autocorr\n",
    "import pandas as pd\n",
    "\n",
    "import statsmodels as sm\n",
    "import seaborn as sns\n",
    "\n",
    "from lib.lif import LIF, ParamsLIF, LSM, ParamsLSM\n",
    "from lib.causal import causaleffect, causaleffect_linear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "q = 3                      #Dimension of learnt vector\n",
    "dt = 0.001                 #Simulation timestep\n",
    "DeltaT = 50               #Number of timebins over which learning rule is applied\n",
    "tsim = 1000                #Total simulation time\n",
    "T = int((tsim/dt)/DeltaT)  #Number of learning blocks\n",
    "Nsims = 3                  #Number of repeated simulations\n",
    "x_input = 0                #Input\n",
    "n = 2                      #Number of neurons\n",
    "sigma = 10                 #Their noise level\n",
    "sigma_perturb = 20        #Perturbation noise level\n",
    "mu = 1                     #Spiking threshold\n",
    "tau = 1                    #Neuron timescale\n",
    "eta = 1                    #Cost gradient learning rate (RDD)\n",
    "epsilon = 2e5              #Weight learning rate (RDD)\n",
    "epsilon_np = 1e-1          #Weight learning rate (NP)\n",
    "epsilon_np_baseline = 1e-1 #Weight learning rate with baseline (NP)\n",
    "tau_s = 0.20               #Output filter timescale\n",
    "\n",
    "c = 0.99                   #Correlation coefficient\n",
    "p = 0.1                    #Learning window\n",
    "\n",
    "N = 19\n",
    "wmax = 20\n",
    "wmin = 2\n",
    "\n",
    "sigma = 10\n",
    "mu = 1\n",
    "tau = 1\n",
    "params = ParamsLIF(sigma = sigma, mu = mu, tau = tau)\n",
    "t = 500\n",
    "lif = LIF(params, t = t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00221199216929\n",
      "0.0025150289023\n"
     ]
    }
   ],
   "source": [
    "print 1./DeltaT*(1 - np.exp(-DeltaT*dt/(tau_s)))/2\n",
    "print exp_filter[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Simulate for 500s, plot ISI distribution for a set of weights\n",
    "w = np.array([10, 12])\n",
    "\n",
    "lif.W = w\n",
    "(v, h, _, _, u) = lif.simulate()\n",
    "\n",
    "#Compute s\n",
    "len_h = h.shape[1]\n",
    "t_filter = np.linspace(0, 1, 2000)\n",
    "exp_filter = np.exp(-t_filter/tau_s)\n",
    "exp_filter = exp_filter/np.sum(exp_filter)\n",
    "ds = exp_filter[0]\n",
    "s1 = np.convolve(h[0,:], exp_filter)[0:len_h]\n",
    "s2 = np.convolve(h[1,:], exp_filter)[0:len_h]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Reshape s1 and h\n",
    "s_r = s1[:].reshape((-1, DeltaT))\n",
    "S = np.squeeze(s_r[:,-1])\n",
    "\n",
    "s_r2 = s2[:].reshape((-1, DeltaT))\n",
    "S2 = np.squeeze(s_r2[:,-1])\n",
    "\n",
    "H_r = h[0,:].reshape((-1, DeltaT))\n",
    "H = np.max(H_r, 1)\n",
    "\n",
    "H_r2 = h[1,:].reshape((-1, DeltaT))\n",
    "H2 = np.max(H_r2, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Is stationary distribution of S1 shifted by a fixed amount in spiking bins?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASkAAAEKCAYAAABOoaMaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl4VdXV/z9fIqg4FETqgAPYRluoFiWvQ2sVBRVQRK0D1llbBNQ6ddBaX7Wt/bW1WkVRREWgBHGoAypWJoG3KmpAREFRQFAQBQdwAJmyfn/sk3q5JPeem9yTe5Osz/OcJ/fs+eSGL3vvs/ZaMjMcx3GKlWaFHoDjOE4mXKQcxylqXKQcxylqXKQcxylqXKQcxylqXKQcxylqXKQcxylqXKQcxylqEhUpST0kzZM0X9JV1eRL0qAof7akA1LyhklaLumNtDoPSpoVXYskzYrS20tak5I3JMlncxynftgiqYYllQCDgaOAJcArksaa2dyUYj2B0ug6CLgr+gkwHLgDGJnarpmdltLHzcCqlOwFZtY57hh33HFHa9++fdzijuPUFTNYuBBWrmQGfGxmbbNVSUykgAOB+Wa2EEDSGKAPkCpSfYCRFs7mTJfUStIuZrbMzKZJal9T45IEnAocWdsBtm/fnoqKitpWdxwnF9auhVNOgZkz4dZb0WWXLY5TLcnlXjvg/ZT7JVFarmVq4ifAR2b2Tkpah2ipN1XST6qrJKmfpApJFStWrIjZleM4dWLNGujTB558Eu68Ey69NHbVhrxxfjrwQMr9MmCPaLl3BTBa0vbplcxsqJmVmVlZ27ZZZ5qO49SVr76C446D8ePh3nthwICcqie53FsK7J5yv1uUlmuZzZC0BXAS0KUqzczWAmujzzMkLQD2Bnw95ziF4osv4Nhj4fnnYcQIOOusnJtIcib1ClAqqYOkFkBfYGxambHA2dFbvoOBVWa2LEbb3YG3zGxJVYKkttFmPZL2ImzGL8zHgziOUwtWrYJjjoEXXoDRo2slUJDgTMrMNki6GHgWKAGGmdkcSf2j/CHAOKAXMB9YDZxXVV/SA0BXYEdJS4DrzOy+KLsvmy71AA4D/iBpPVAJ9DezT5N6PsdxMvDZZ0GgZs2Chx6Ck06qdVNqyk7vysrKzN/uOU6e+fhjOOoomDsXHnkEeveutpikGWZWlq25JPekHMdpanz0EXTvDvPnwxNPQI8edW7SRcpxnPywbBkceSQsXgxPPQXduuWlWRcpx3HqzpIlQaCWLYN//xsOOyxvTbtIOY5TNxYtCgL1ySfw7LPwox/ltXkXKcdxas+CBUGgPv8cJk6E//mfvHfhIuU4Tu2YNy8I1Nq1MHky7L9/It24SDmOkztz5waBMoPnnoN9902sq4Z8ds9xnEIwezZ07QrNmsGUKYkKFLhIOY6TCzNnwhFHwJZbwtSp8P3vJ96li5TjOPF46aWwxNtuuyBQpaX10q2LlOM42Xn++XDUpU2bIFB77VVvXbtIOY6TmSlTwmHhXXaBadNgzz3rtXsXKcdxambiROjVKwjT1KnQLq7j3PzhIuU4TvWMGxc8apaWhtnUzjsXZBguUo7jbM4TT8AJJ0CnTsFQs4Cutl2kHMfZlEcegZNPDhbkkyaFzfIC4iLlOM43jB4NffvCQQfBhAnQqlWhR+Qi5ThOxIgRcOaZcOihwd3K9psFWyoILlKO48A998B55wVHdePGwbbbFnpE/8VFynGaOoMHQ79+wdXvk09Cy5aFHtEmJCpSknpImidpvqSrqsmXpEFR/mxJB6TkDZO0XNIbaXWul7Q0ilQ8S1KvlLyro7bmSTomyWdznEbBP/4BF18cogs/9hhstVWhR7QZiYlUFANvMNAT6AicLqljWrGehPh4pUA/4K6UvOFATV7c/2FmnaNrXNRfR0Koq05RvTur4vA5jlMNf/kLXHFFeJP38MPh0HARkuRM6kBgvpktNLN1wBigT1qZPsBIC0wHWknaBcDMpgG5xM3rA4wxs7Vm9i4hlt+BdX4Kx2lsmMEf/gBXXw0/+xk88AA0b17oUdVIkiLVDng/5X5JlJZrmeq4JFoeDpPUuo5tOU7TwQyuvRauuw7OOQdGjoQtitv3ZUPcOL8L2AvoDCwDbs6lsqR+kiokVaxYsSKJ8TlOcWIGv/kN3Hgj/OIXMGwYlBT/jkiSIrUU2D3lfrcoLdcym2BmH5nZRjOrBO7hmyVdrLbMbKiZlZlZWdsCmvo7Tr1iBpddBn//O1x0EQwZEjxrNgCSHOUrQKmkDpJaEDa1x6aVGQucHb3lOxhYZWbLMjVatWcVcSJQ9fZvLNBX0paSOhA241/Ox4M4ToOmshIGDIBBg+Dyy+H22xuMQEGCgRjMbIOki4FngRJgmJnNkdQ/yh8CjAN6ETa5VwPnVdWX9ADQFdhR0hLgOjO7D/ibpM6AAYuAC6P25kh6CJgLbAAuMrONST2f4zQINm4MS7v774erroI//xmkQo8qJ2RmhR5DwSgrK7OKiopCD8NxkmHDBjj3XCgvDxvl111XVAIlaYaZlWUrV9zb+o7j1I7168M5vIceChvlv/tdoUdUa1ykHKexsW5d8GTw2GNho/zKKws9ojrhIuU4jYmvvw4W5E8/HTbKL7mk0COqMy5SjtNYWL0aTjwRxo8PJgYXXljoEeUFFynHaQx89RX07h18kQ8bFtyuNBJcpBynofPFF3DssSE23siRYcO8EdFwLLocx9mclSvh6KPhhRfCQeEiFKjycmjfPtiPtm8f7nPBZ1KO01D59NMgULNnB1crJ55Y6BFtRnl58Ke3enW4X7w43OeCi5TjNEQ+/hi6d4c334RHHw3x8YqQa675RqCqWL06pMfFl3uO09D46CPo2hXmzQvufgsgUOXlwQ26FK6SEhg4cPNyixdXX7+m9OpwkXKchsQHHwSBevfdYAt19NGJd5m+pzRwIJx9dnihWEVlJdx11+ZCVZMnmFw8xPjZPT+75zQU3n8fjjwSPvwwRHT5yU8S7a68HC69FD75ZNN0KXh+qY6SknBkMLVszcQ7u+czKcdpCLz7Lhx2GCxfHow160Ggzj13c4GCmgUKgtOFVPbcs/pyNaVXh4uU4xQ78+fD4YfDqlUh7PkhhyTeZf/+m86I4pK+jLvxxs0jZLVsGdLjklWkJO0k6T5Jz0T3HSVdEL8Lx3FqzVtvBYFaswYmT4ayrKujOlG1//Tll7Wrn25ecMYZMHRomDlJ4efQoSE9LnFmUsMJjut2je7fBi6L34XjOLXijTfCJvnGjfDcc9C5c6LdlZfD+ednf/PWsmVw9LnNNt+kNWsW0u68c/PyZ5wBixaFzfVFi3ITKIgnUjua2UNAJQSPm4B7vHScJHntNTjiiPCvf8oU+MEPEuuqvBx23DEYq69bl7380KFBjL78MuxPmQUdrU6g8kEcY86vJLUhuOulyhd5MsNxHIcZM+Coo8JUZfJkKC1NrKvy8hDZKn3DuyYGDMh9JlRX4ojUFYQgB9+R9DzQFjg50VE5TlNl+nTo0QNatw4C1aFDYl2Vl8NZZ2V+W1dFSUnYb0pqtpSJrCJlZjMlHQ7sAwiYZ2brEx+Z4zQ1/vMf6NkTdtopCNQeeyTWVadOMHduvLJt2oRTOIWiRpGSdFINWXtLwsweTWhMjtP0eO65cLxl992DmUG75IJvd+8eX6CaN4fbbktsKLHINJPqnSHPgKwiJakHcBshpNW9ZvaXtHxF+b0IIa3ONbOZUd4w4DhguZn9IKXOTdHY1gELgPPMbKWk9sCbwLyo6HQz659tjI5TcMaPhz59YK+9gkDtvHOi3U2aFK9cmzZBoOp7DyqdGkXKzOrk2k9SCTAYOApYArwiaayZpWp4T0IQz1LgIEII9YOivOHAHcDItKYnAFdHcf3+ClwN/DbKW2Bmyb6ndZx88vTTcNJJ8P3vw4QJUCRRtbfZprBLvFTiGHO2kTRI0kxJMyTdFr3ty8aBwHwzW2hm64AxQJ+0Mn2AkRaYDrSqilBsZtOAT9MbNbPxkRkEwHRCOHXHaXg8/njwAbXvvmEPKgGBGjgQttgiGFJusUVY6mWjpATuvjvvQ6k1ceykxgArgJ8S3uqtAB6MUa8d8H7K/ZIoLdcymTgfeCblvoOkWZKmSkr2cJPj1IWHH4ZTToEuXWDiRNhhh7x30b178ExQZV6wcWP2pV6bNjBiROGXeKnEMUHYxcz+mHL/J0mnJTWguEi6hhBOvcoZ6TJgDzP7RFIX4HFJnczs87R6/YB+AHsk+PbEcWqkvDz4OvnRj8Jyb/vt89p8y5bhFE1tKJYlXipxZlLjJfWV1Cy6TiUck8nGUmD3lPvdorRcy2yGpHMJm+pnWORrxszWmtkn0ecZhE31vdPrmtlQMyszs7K2RbL+d5oQw4cH46TDD4dnnsmrQHXvHpZ1tRWoIorAvglxROoXwGjC27SqvaULJX0h6fMM9V4BSiV1kNQC6EswCk1lLHC2AgcDq8xsWabBRG8MfwMcb2arU9LbRpv1SNqLsBm/MMbzOU79MHRoCDXVvTs89VRwbZknuneP/9auJvoX67twM0vsIpgWvE2Y1VwTpfUH+kefRXgDuAB4HShLqfsAYQm3nrBXdUGUPp+wjzUruoZE6T8F5kRpM4He2cbXpUsXc5x64fbbwzG3Y481W7Mmb82OGmW2zTZVJ+jiXdLmad265W1IsQEqLIaOxPLMKWk/oD0pe1jWCIw53TOnUy/cfDP86ldwwgnw4IPQokWdm8zFYjydUaNCIIT33gtG7TfeWJiNcimeZ86sG+eRUeV+hFlKZZQcy5jTcZo8f/5zUIRTTgkb5s2b17nJughU1QHhYnp7l404b/cONrOOiY/EcRoTZnDDDeE644ywYb5FfiLI1UWgCnFAuK7E2Th/UZKLlOPExSzMnm64ITgKHzEibwLVunXudQYMCENqiAIF8WZSIwlC9SGwlrDZbWa2X6Ijc5yGiFnYf7rlFrjwwqAMzfIXSmDlytzKN9TZUypxROo+4CzC27fKLGUdp+lSWRliQN1xB1xySTidm0fjo1z226VgUtDQBQriidQKM0u3b3IcJ5XKyqAK99wDV14JN92UV4Fq2RLWZ/HittVWcO+9DWtTPA5xROpVSaOBJwnLPaBxmCA4Tl7YuBF+/vOwOf6738Gf/pRXgYrbVG0tzYudOCK1NUGcUuM5uwmC40AITnfOOTB6dNgov/bagghUYyaO++A6+ZVynEbL+vVhbfXww8Ee6uqr89b0wIHBg0FcunXLW9dFRxxjzq2AC4BOwFZV6WZ2foLjcpziZu1aOO00eOKJYFF+xRV5a7pdO/jgg/jld901eHtprMR5N/pPYGfgGGAqwVPBF0kOynGKmq+/Dt40n3gCbr89rwLVsmVuAtW8OSzN6jekYRNHpL5rZtcCX5nZCOBYvnHx6zhNi9Wr4fjjg5uVu++Giy/OW9OdOuW++R0nmGdDJ87GedWLz5WSfgB8CHw7uSE5TpHy5ZfQuzdMnQrDhgVr8jyS63GXOPHyGgNxRGqopNbAtQT/T9tGnx2n6fD559CrF7z4YnAj8LOf5bX5XN7itWoFn32W1+6LmqzLPTO718w+M7OpZraXmX3bzIrITbvjJMzKlXD00fDSSzBmTEEFasCApiVQkDk4aG9gtpktju7/l+BYbjFwqZm9Wz9DdJwC8umncNRR8Prr8MgjIT5eHikpiV+2qSzv0sk0k7qREBkGSccBZxKis4wFhiQ/NMcpMCtWwBFHwJw5IfxUngUKwmmaODRVgYLMImX2jQ/xk4D7zGyGmd0LeAQDp3Hz4YfQtSu88w48+WTYjyoQW29dsK6LgkwiJUnbSmoGdANS3bxvVUMdx2n4LF0aorksXgzjxoXlXp6piuwSh9Wrs5dpzGR6u3crIajB58CbZlYBIGl/QoAEx2l8vPceHHkkLF8O//43HHpo3rto3Tq+X6imvMyrokaRMrNhkp4l2ES9lpL1IeDn+ZzGx7vvBoH67DOYMAEOyr/NcvfuLlC5ktEEwcyWmtmrZlaZkrbMzN6L07ikHpLmSZov6apq8iVpUJQ/W9IBKXnDJC2X9EZanR0kTZD0TvSzdUre1VFb8yQdE2eMjgOEvafDDgv2UJMmJSJQUPfYeE2R/Pk1TSMK1DkY6Al0BE6vxld6T0IQz1JC6PPUc9/DgR7VNH0VMMnMSgn7ZFdF/XUkBCDtFNW7sypYqONk5M03wx7U11/D5MnQpUsi3eRiD9WYvRrkSmIiBRwIzDezhWZWFfk4/R1uH2BkFCtwOtBK0i4AZjYN+LSadvsAI6LPI4ATUtLHWAi3/i4hiOiBeX0ip/HxxhvhLV5lJUyZAj/8YSLd5CJQzZo1bq8GuZJVpCT9M05aNbQjRBquYkmUlmuZdHayb0KxfwjslEtbkvpJqpBUsWLFiixdOY2aWbOCQG2xRTiP16lTIt3keuRl48ZEhtFgiTOT2uSbi5ZQycyHcyQK1ZzT9qKZDTWzMjMra9vWzb2aLBUVYZO8ZcsgUPvsk/cucjEzgLBR3tSOvMQh07GYq4HfAVtL+rwqGVgHDI3R9lJg95T73aK0XMuk85GkXcxsWbQ0XF6HtpymyIsvQo8e0KZN2INq3z7vXbRsmZvblQED8j6ERkONMykz+39mth1wk5ltH13bmVkbM4vjJ/UVoFRSB0ktCJva6VFnxgJnR2/5DgZWpSzlamIscE70+RzgiZT0vpK2lNSBsBn/coxxOk2JadPCYeFvfzvMoBIQqNr4hWoMoaeSIs5y7ylJ2wBIOlPSLZL2zFbJzDYAFwPPAm8CD5nZHEn9JfWPio0DFhI2ue8BBlbVl/QA8CKwj6Qlki6Isv4CHCXpHaB7dI+ZzQEeAuYC/wYuMjNf3TvfMHky9OwJu+0WBGr33bPXqQXuFyq/yLL8hiTNBn4I7EcwC7gXONXMDk98dAlTVlZmFRUVhR6GUx88+yyccAJ897vh1dlOO2WvUwtyje7SlAVK0gwzK8tWLs5MakO0Qd0HuMPMBgPb1XWAjlNvPPVUcPn7ve/Bc88lJlDl5bmVb8oClQtxPHN+EW2inwkcFh04bp7ssBwnTzz2WIjq8sMfhtnUDjsk0k0u5/GaNXMzg1yIM5M6jRAc9AIz+5Dw1uymREflOPngwQfhlFOCBfnEiYkJlJTbeTwXqNyIExz0Q+CWlPv3gJFJDspx6syoUSGy8I9/DE8/Ddsls0ORyx7UqFGJDKHRU6tjMZJez/dAHCdvDBsGZ58drMmfeaYoBApCsGMndzIZc55UUxYhWKjjFB9DhgTLyGOOCftRCbm19Ld49Uem5d6DQDnVHztxz5xO8TFoEFx6KRx3HDz8MGyVzJ9prgLly7y6kUmkZgN/N7M30jMkdU9uSI5TC266CX7zGzjxxBB2qkWLRLrJJboLBIHyZV7dyCRSlxFcB1fHiQmMxXFqx403wu9/H0wN/vlPaJ6MhUx5efzoLuBLvHyRyX3w/2XIczNtp/CYwfXXwx/+AGedFTbMt4hj+pc7nTrldtzFBSp/JPONOk7SmMHVV8Nf/wrnnw9Dh+a+FouJb5IXFhcpp+FhBldcAbfeCv37w+DBwYw7AVygCk+N36ykS6OfP66/4ThOFior4eKLg0D98pfBx0kCAiW5QBULmb7dqrBVt9fHQBwnK5WVcOGFQZh+/esgVLkqSQxq06QLVHJkWu69Gfls2jVy11KFCJ5790t2aI6TwsaNcMEFMGJEeJP3hz8UjUA19TDoSZPp7d7pknYmOK07vv6G5DhpbNgQjrk88EAQp2uvLfSI/svWW3sY9KTJFhz0QzP7ISGs+nbR9YGZLa6PwTkO69ZB375BoP7yl8QEqry8dntQLlDJk/XtnqTDCV4PFhGWertLOieKi+c4ybF2LZx6KowdC7fcApdfnkg35eVw5pm51fE9qPojjgnCLcDRZjYPQNLewAMUSVgrp5GyZg389KfBi8HgwTBwYPY6tSBXgXJxqn/iiFTzKoECMLO3JblnTic5Vq+GPn1g0iS45x74+c8T6cYFqmEQx8CkQtK9krpG1z1ArGMxknpImidpvqSrqsmXpEFR/mxJB2SrK+lBSbOia5GkWVF6e0lrUvKGxBmjU2R8+SX06hUiuwwfnphAQW4C5Z4MCkecmdQA4CLgl9H9/wFZo4RFkY4HA0cRQp6/ImmsmaWegOpJiI9XChwE3AUclKmumZ2W0sfNwKqU9haYWecYz+QUI59/HkJOvfRSUIXTT0+km1z8kUN4g+eeDApHHPfBawn7UrdkK5vGgcB8M1sIIGkMIeJMqkj1AUZG0WimS2oVRSVun62uJAGnAkfmOC6nGPnssxBVeObM4Jv8pz9NpJtc3+C1auWhzwtNMgeeAu2A91Pul0RpccrEqfsT4CMzeyclrUO01Jsq6Sd1GbxTj3zyCXTrBrNmwb/+VTQCtfXWLlDFQEM+YHw64S1jFcuAPczsE0ldgMcldTKzTXxiSeoH9APYY4896m2wTg0sXw7du8Pbb8Pjj4flXp4pKcnNDxS4kWYxkdNMSlIzSdvHLL4USI1jvVuUFqdMxrqStgBOIrg4BsKy1Mw+iT7PABYAe6cPysyGmlmZmZW1bds25qM4ibBsWQiWMH9+COCZgEBJuQsUuEAVE1lFStJoSdtL2gZ4A5gr6dcx2n4FKJXUQVILoC8wNq3MWODs6C3fwcAqM1sWo2534C0zW5IyzrbRhjuS9iJsxi+MMU6nECxZAocfDu+9F2yhuuffI3Vtj/a5qUFxEWcm1TFaMp0APAN0AM7KVsnMNgAXE87+vQk8ZGZzJPWX1D8qNo4gJPOBe4CBmeqmNN+XTZd6AIcBsyOThEeA/mb2aYznc+qbxYuDQH34YYgqfPjhee+iNgLVrZsLVDEiy/KtSJoDdAZGA3eY2VRJr0Vn+ho0ZWVlVlHhnpDrlYUL4YgjgrnBs8/CgQfmvQt3tdIwkDTDzMqylYszk7qbcG5vG2CapD2pOUCD49TM22/DYYcFg81Jk1ygnFjEsZMaBAxKSVos6YjkhuQ0St58E448MviFeu452C//7shat869jgtU8RNn43wnSfdJeia67wick/jInMbD669/s+80ZUoiAiXlZkVu5gLVUIiz3BtO2MDeNbp/mxCTz3Gy8+qrYQ+qRQuYOhU6dsx7F+6LvHETR6R2NLOHgEr475u3jYmOymkcvPxyWOJts00QqL03M1urMy5QjZ84IvWVpDaAAVTZMyU6Kqfh88ILwfapdWuYNg2+8528d+EC1TSIcyzmCoIh5XckPQ+0BU5OdFROw2batOBuZdddg8uV3XbLexcuUE2HOG/3ZkYuhPchuA+eZ2brEx+Z0zCZNAl694b27cPnXXbJexe5ClQC22BOPVKjSEk6qYasvSVhZo8mNCanofLvf8OJJ0JpKUycCN/+dl6br81B4Y4dYc6c7OWc4iXTTKp3hjwDXKScb3jySTj5ZOjUCSZMgDZt8tq8G2k2XTLF3TuvpjzH2YR//SuEndp//3DUpTZWlRkoKcm9jgtU4yHTcu9MMxsl6Yrq8s0sV0+dTmNkzJjgLPygg2DcOPjWt/LafIsWuS/xXKAaF5mWe9tEP7erJs//DBwYORLOOw8OPTT4g9quuj+V2uNLPAcyL/fujj5ONLPnU/Mk/TjRUTnFz333wS9+EYw1n3giGGzmERcop4o4xpy3x0xzmgp33hlCTR1zTNgwL7BAtWrlAtWYybQndQjwI6Bt2r7U9kAttjKdRsGtt4Zw5717w8MPw5Zb5rV5N9J00sm0J9UC2DYqk7rZ8Dlucd40+dvf4Le/DdFcRo8Ou9p5pLbufp3GTaY9qanAVEnDzWxxPY7JKUb++Ef43/8Npgb//Cdskd9AQ74H5dREnL+0LSUNJQTs/G95M/OgnE0BsyBOf/oTnHUW3H9/7QyXMuBLPCcTcUTqYWAIcC/uoqVpYRaWdzfdBBdcAHff7QLl1DtxRGqDmd2V+Eic4sIsbJDfdhsMGAB33AHN8hvw2gXKiUOcv7onJQ2UtIukHaquOI1L6iFpnqT5kq6qJl+SBkX5syUdkK2upOslLY3Cqc+S1Csl7+qo/DxJx8QZo1MNlZVw0UVBoC67DAYPLrhAjRqV1+6dhoSZZbyAd6u5FsaoV0KIIrwX4U3ha4QYfqllehFi+Qk4GHgpW13geuBX1fTXMSq3JSE24AKgJNMYu3TpYk4aGzaYXXBBcAH+29+aVVbmvYtvPIzHu0aNyvsQnCIAqLAsOmJmsfxJdail/h0IzDezhQCSxgB9gLkpZfoAI6MBT5fUStIuhE36bHXT6QOMMbO1wLuS5kdjeLGW4296bNgA558f3t5dey3ccEPe7QJ8iefkSqz3yJJ+QJipbFWVZmYjs1RrB7yfcr8EOChGmXYx6l4i6WygArjSzD6L6kyvpq30Z+kH9APYY489sjxCE2L9ejj77HBg+I9/hN//Pu9duEA5tSFOSKvrCMdgbgeOAP4GHJ/wuDJxF2EZ2BlYBtycS2UzG2pmZWZW1rZt2yTG1/BYty7YP40ZEww2XaCcIiLObujJQDfgQws+pn4IxPHHsRTYPeV+tygtTpka65rZR2a20cwqgXsIS7q4/TnprF0bnNU9+mg48vLrX+e9i1wEauutXaCcTYkjUmsiQdggaXtgOZuKQU28ApRK6iCpBdCXENAhlbHA2dFbvoOBVWa2LFPdaM+qihOBN1La6itpS0kdgFLg5RjjbLqsWQN9+oRDwnfeCZdemvcucp1BrV6d9yE4DZw4e1IVkloRZi0zgC+JsRltZhskXUwILFoCDDOzOZL6R/lDgHGEN3zzgdXAeZnqRk3/TVJngk+rRcCFUZ05kh4ibK5vAC4yMzc+rYmvvoLjjw8hz++9Nxhr5hlf4jn5QJbDX4ak9sD2ZjY7qQHVJ2VlZVZRUVHoYdQ/X3wBxx4Lzz8Pw4eH4y55JtegCS5QTQ9JM8ysLFu5rDMpSYdVl2Zm02o7OKeArFoFPXuG6MKjR8NppyXSjQuUky/iLPdSd1K3ImxUzwD8gHFD47PPgqO6WbPgoYfgpJqiltWNXJZ5LlBONuIYc24S2kp9znTkAAARYElEQVTS7sCtiY3ISYaPP4ajjoK5c0N0l96ZIpbVHhcoJ9/UxinQEuD7+R6IkyDLl0O3bjB/fvBH3qNHIt24QDlJEGdP6na+iQ7TjGBEOTPJQTl5ZNmyIFCLFoWILt26JdKNC5STFLFMEFI+bwAesLToMU6RsmRJiOaybFkIgX7YZu9A8kK7zQ4f1UxCGuk0YuI6vftu9HledIDXKXYWLQoC9cknMH48HHJIIt20bg0rV8YvP3FiIsNwGjE1WpxLai7pVsJB3/uB4cDCKt9OkUGlU4wsWACHHx7e5k2cmJhASbkJlC/znNqQaSZ1M9ASaG9mXwBEx2L+LukuoAfBb5NTTMybF2ZQa9fC5Mmw//6JdOPW5E59kUmkegGllmKSbmafSxoAfAz0THpwTo7MnRsEyiwcd9l330S6cYFy6pNMB4wrrZozM9F5uBVmNr2aOk6hmD0bunYNbn6nTHGBchoNmURqbuRYbhMknQm8mdyQnJyZOROOOCJEE546Fb6ffzO2du1coJzCkGm5dxHwqKTzCcdgAMqArQkuUpxi4KWXwlGXVq3CHtRee+W1+e7dYdKk3Ou5QDn5IlME46XAQZKOBDpFyePMrBZ/sk4iPP98OCzctm0QqD33zGvzLlBOMRDn7N5kYHI9jMXJhSlT4Ljjwjps8uTcLCpj4gLlFAP5Dabm1A8TJ0KvXmHmNHVqIgLVokVu5Xfd1QXKSQYXqYbGuHFhBlVaGmZTO++c9y5atAjBY+LSsSMsdW/yTkK4SDUknngCTjgBOnUKS7wEot1IuQnUqFEwZ072co5TW2rjqsUpBI88AqefDgccAM8+G97m5ZF27eCDD+KXb9UqnLpxnKTxmVRDYPToEBfvoINgwoS8C5SUm0CNGuUC5dQfiYqUpB6S5kmaX3UwOS1fkgZF+bMlHZCtrqSbJL0VlX8simSDpPaS1kiaFV1Dkny2emPECDjzTDj00OBuZfvt89p8baKon3FGXofgOBlJTKQklQCDCWf8OgKnS+qYVqwnIT5eKSH0+V0x6k4AfmBm+wFvA1entLfAzDpHV/9knqweueceOO+84IRp3DjYdtu8Nl8bgRo1Kq9DcJysJDmTOhCYb2YLzWwdMAbok1amDzDSAtOBVlHwzxrrmtl4M9sQ1Z9OiFTc+Bg8GPr1C65+n3wSWrbMa/O1FSifRTn1TZIi1Y7gi6qKJVFanDJx6gKcDzyTct8hWupNlfST2g684PzjH3DxxSG68GOPwVZb5bX5XM2qOnYMNlAuUE4haLAb55KuIbgzLo+SlgF7mFln4ApgdOT/Kr1eP0kVkipWrFhRfwOOy1/+AldcASefDA8/HA4N55lcNskHDHATA6ewJGmCsBTYPeV+tygtTpnmmepKOhc4DuhW5U4mcmu8Nvo8Q9ICYG829dGOmQ0FhkKIYFy7R0sAM/jjH+G66+BnPwsb5lvk/+uJO4tyEwOnWEhyJvUKUCqpg6QWQF9gbFqZscDZ0Vu+g4FVZrYsU11JPYDfAMeb2eqqhiS1jTbckbQXYTN+YYLPlz/M4Nprg0Cdcw6MHJmIQEG8WdSAAS5QTvGQ2EzKzDZIuhh4FigBhpnZHEn9o/whwDiCB9D5wGrgvEx1o6bvALYEJijs/k6P3uQdBvxB0nqgEuhvZp8m9Xx5wwx+8xv4+9/hF7+AIUOC47oC0bEj3Hlnwbp3nM1QNc43mwxlZWVWUVGRvWBSmMFll8GgQXDRReFnggLVqVPwMJxtSI5TH0iaYWZl2cr5sZhCUVkJAwfC3XfD5ZfDzTfXzi4gJgMHZheoAQMS695xao2LVCHYuDEs7e6/H666Cv7850QFKs4MCnyZ5xQnDdYEocGyYUPYHL///rBRXiQC5bMop1jxmVR9sn59OIf30ENw443wu98l2l15eTyBAp9FOcWLi1R9sW5d8GTw2GPhTd6VVybe5TXXxCvXrVuy43CcuuAiVR98/XWwIH/66fAG75JL6qXbxYuzl9l11+CN2HGKFReppFmzJnjTHD8+2EBdeGG9dNu6dfYyu+7qbn+d4sc3zpPkq6/g2GODo7phw+pNoMrLYeXKzGXcL7nTUPCZVFJ88UUQqOefD8dczjwz710MHAhDhwaLhiqk7AaZbdr4oWGn4eAilQQrV4agna+8Ag88AKeemvcuajItiGMxfttteR+O4ySGi1S++fRTOPpomD07uFo5Mf8R6eNYj9fE1lu7XyinYeEilU8+/jjEJn/zTXj00RAfLwGG1MF7++rV2cs4TjHhIpUvPvooGBwtWBDc/R59dCLdDBxY+0PAbg/lNERcpPLBBx8EBXjvvWALdeSRiXV1993xyqVvoHfr5vZQTsPETRDqyvvvw+GHw5IlIeRUggIFwXlCJkpKwjm8ysogUlWXC5TTUPGZVF14990gSp9+Gow1DzmkoMNxX1BOY8RnUrVl/vwwg1q1CiZNqjeB2mab3NIdp6HjIlUb3norCNSaNTB5MpRldS6YN+6+e3Pnnc2axd+rcpyGhotUrrzxBnTtGsy8n3sOOneu1+7POCMYsO+5Z9gc33PPcO+2T05jxfekcuG114IdVPPmYQb1ve8VZBhnnOGi5DQdfCYVlxkz4IgjQjThqVMLJlCO09RIVKQk9ZA0T9J8SVdVky9Jg6L82ZIOyFZX0g6SJkh6J/rZOiXv6qj8PEnH5O1Bpk8Phkbf+hZMmwalpXlr2nGczCQmUlGgzsFAT6AjcLqkjmnFehKCeJYC/YC7YtS9CphkZqXApOieKL8v0AnoAdxZFSy0TvznP3DUUbDjjmEG1aFDnZtMZ+DAEAtUCj8HDty8THk5tG8fNsnbtw/3jtMUSHImdSAw38wWmtk6YAzQJ61MH2CkBaYDrSTtkqVuH2BE9HkEcEJK+hgzW2tm7xICjh5Ypyd47jk45pgQm3zqVNhjjzo1Vx0DB8Jdd33jbmXjxnCfKlTl5dCvX/C0aRZ+9uvnQuU0DZIUqXbA+yn3S6K0OGUy1d0pCsUO8CGwUw79xWf8eOjVK0xbpkwJQpUAQ4dmT7/mms0PBq9eHd+HueM0ZBr0xrmF8Ms52VlL6iepQlLFihUrqi/09NPQuzfss08QqJ13rvtgayDVYV1N6e+9V32ZmtIdpzGRpEgtBXZPud8tSotTJlPdj6IlIdHP5Tn0h5kNNbMyMytr27bt5qN+/PHgA2rffYOZQXVl8khJDbtmqek1rTITWH06TtGRpEi9ApRK6iCpBWFTe2xambHA2dFbvoOBVdFSLlPdscA50edzgCdS0vtK2lJSB8Jm/Ms5jfjhh+GUU6BLl3Aid4cdcqpeG/r1y55+443QsuWm+S1bhnTHafSYWWIX0At4G1gAXBOl9Qf6R59FeIu3AHgdKMtUN0pvQ3ir9w4wEdghJe+aqPw8oGe28XXp0sX+y6hRZs2amR16qNmqVVafDBhgVlIS/BWUlIT7dEaNMttzTzMp/Bw1ql6H6Dh5B6iwGDoia8JH58vKyqyiogKGD4fzzw/HXcaOhW23LfTQHKfRI2mGmWU9+NqgN87zwtChcN554bjLU0+5QDlOkdG0RWr58hAL79hjwwwqfePHcZyC07RF6v33Q3ThRx8NZ/Icxyk6mvSelKQVwOIsxXYEPq6H4RQLTe15wZ+5UOxpZlltfJq0SMVBUkWczb3GQlN7XvBnLnaa9nLPcZyix0XKcZyixkUqOzUcAW60NLXnBX/mosb3pBzHKWp8JuU4TlHTpESqvt0ZFwMJPfP1kpZKmhVdverreeJQx2ceJmm5pDfS6hTt95zQ8xbPdxzngF9juIASwuHjvYAWwGtAx7QyvYBnCAefDwZeylYX+BtwVfT5KuCvhX7Wenjm64FfFfr58v3MUd5hwAHAG2l1ivJ7TvB5i+Y7bkozqfp2Z1wMJPXMxUxdnhkzmwZ8Wk27xfo9J/W8RUNTEqn6dmdcDCT1zACXREuHYcW09KFuz5yJYv2ek3peKJLvuCmJVOJYmCc3hdeldxGWF52BZcDNhR1O/dJEvuei+Y6bkkjVtzvjYiCRZzazj8xso5lVAvdQ16g8+aUuz5yJYv2eE3neYvqOm5JI1bc742IgkWeu+scacSLwBsVDXZ45E8X6PSfyvEX1HRd6574+L+rZnXExXAk98z+jsrMJ/wB2KfRz5vGZHyAsb9YT9m4uKPbvOaHnLZrv2C3OHccpaprScs9xnAaIi5TjOEWNi5TjOEWNi5TjOEWNi5TjOEWNi1QjQdJukp6ITukvlHSHpC0T6KerpB+l3PeXdHb0ebikk/PdZ0pfD0THNC5PS79C0twob5KkPXNst6ukp6pJbympXNLrkt6Q9B9JiQZmlPRlku03RFykGgGSBDwKPG5mpUApsDXh5H6+6Qr8V6TMbIiZjUygn02QtDPwP2a2n5n9Iy37VYLtz37AI+TvuS8FPjKzfc3sB8AFBHsipx5xkWocHAl8bWb3A5jZRuBygpXxtpLOlXRHVWFJT0nqGn2+S1KFpDmSbkgps0jSDZJmRjOJ70lqTzASvDzyMfSTyO/Qr9IHJKmLpKmSZkh6NuVIyS9TZj1jqqm3laT7oz5flXRElDUeaFfVb2odM3vOzFZHt9MJxz6qZkhTJD0i6a1oVqQor0eUNhM4qYbf6y6kHB8xs3lmtjaq/3j0bHMk9UsZ/5eSborSJ0o6MBrDQknHR2XOjWa9U6KZ73XVdS7p15JeiX5XN1RXpklQaGtZv+p+Ab8E/lFN+quEA6LnAnekpD8FdI0+7xD9LAGmAPtF94uAS6LPA4F7o8/Xk+JnKPUeGA6cDDQHXgDaRumnAcOizx8AW0afW1Uz5itTyn4PeA/YCmhPms+jGn4XdwC/jz53BVYRRKsZ8CJwaNTe+4QZp4CHgKeqaasz4Yzei8CfgNKUvKrf29aEIyNtonsDekafHyOIa3Pgh8CsKP1cgpV3m5T6ZVHel9HPowl+yBWN/SngsEL/rRXi8pmUc2o0m3gV6AR0TMl7NPo5gyAScdkH+AEwQdIs4PdEsxvCMYtySWcCG6qpeygwCsDM3iIEb907TqdRm2XATSnJL5vZEgsHZWdFz/E94F0ze8eCIoyqrj0zm0XwBHATsAPwiqTvR9m/lPQaYea2O0HwANYB/44+vw5MNbP10ef2Kc1PMLNPzGwN4fd8aFr3R0fXq8DMaMylNEG2KPQAnLwwlzCD+S+Stgd2BuYRBCP1P6StojIdgF8R9no+kzS8Ki9ibfRzI7n9rQiYY2aHVJN3LMEbZG/gGkn7mll1YpUTkroD1wCHW7Qki0j9nOtzYGZfEkTkUUmVQC9JOwHdgUPMbLWkKXzze1sfCR9AZVX/ZlYpKbXv9PNo6fcC/p+Z3Z3LeBsjPpNqHEwCWqa8ZSsh+P+5I/qfehHQWVIzSbvzjduN7YGvgFXRP7yeMfr6AtguS5l5QFtJh0TjaS6pk6RmwO5m9hzwW+BbQPrbsv8Dzojq7Q3sEbVXI5L2B+4GjjezOC5U3gLaS/pOdH96De3+WJGzNwUPAx0JM7tvAZ9FAvU9gkveXDlKwW/61gQvn8+n5T8LnF/1NlFSO0nfrkU/DR6fSTUCzMwknQgMlnQt0BZ40MxujIo8D7xLmHG9SVg+YGavSXqV8I/2fTb/h1IdTwKPSOoDXFLDeNYpmCIMkvQtwt/ZrYST+qOiNAGDzGxlWvU7gbskvU5YDp5rZmuj/e6auIkgdg9H5d4zs+NrKmxmX0eb3U9LWk0QxuqE9zvRWKr2hZ4G/kXwJd5f0psEAZ2eaXA18HLU1m7AKDOrSBvj+Ghp+WL0TF8CZ1I8fqzqDfeC0AhRsGN6ADjRzGYWejzOpkg6l7BRfnGhx9IQ8JlUI8TMXgByMmh0nGLFZ1KO4xQ1vnHuOE5R4yLlOE5R4yLlOE5R4yLlOE5R4yLlOE5R4yLlOE5R8/8BWjyNC9xzZ44AAAAASUVORK5CYII=",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6594f34810>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nsamps = min(np.sum(H == 0), np.sum(H == 1))\n",
    "fig,ax1 = plt.subplots(1,1,figsize = (4,4))\n",
    "sm.graphics.gofplots.qqplot_2samples(np.random.choice(S[H == 0], nsamps), np.random.choice(S[H == 1], nsamps), line='45', ax = ax1)\n",
    "plt.savefig('./fig_s1_a.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAFpCAYAAADZUI/8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAFaVJREFUeJzt3X2sZVd5H+DfG0/gjyQVYN9aLjAaj2SQcNROyMi11YBwKYlDSQypRGxFERSUwVJKk7RVxYAUUCrkBOK4qpqABmFhKjBGdUhQRT4myAJVMh9j5Dq2wcYejLE12JNxVWiLaGze/nH3ZI6HO3PH99wzd809zyMd3X3W3nuddZb2LP1m7X32ru4OAABj+JGtbgAAACcIZwAAAxHOAAAGIpwBAAxEOAMAGIhwBgAwEOEMAGAgwhkAwECEMwCAgQhnAAAD2bHVDUiSCy64oHft2rXVzQAAWNedd975N929sqj6hwhnu3btyqFDh7a6GQAA66qqby6yfqc1AQAGIpwBAAxEOAMAGIhwBgAwEOEMAGAgwhkAwECEMwCAgQhnAAADEc4AAAYinAEADEQ4AwAYiHAGADAQ4QwAYCA7troBwDZ1+/WLqffK/YupF2AQZs4AAAYinAEADGTdcFZVN1XVE1V1z0zZrVV11/R6uKrumsp3VdX3ZtZ9cJGNBwDYbs7kmrOPJPnPST56vKC7f/n4clXdkOR/zWz/UHfv2awGAgAsk3XDWXd/vqp2rbWuqirJG5P8081tFgDAcpr3mrNXJHm8u78+U3bxdErzc1X1ijnrBwBYKvPeSuPaJLfMvD+SZGd3H6uqn07yJ1V1aXd/5+Qdq2pfkn1JsnPnzjmbAQCwPWx45qyqdiT5pSS3Hi/r7u9397Fp+c4kDyV5yVr7d/eB7t7b3XtXVlY22gwAgG1lntOa/yzJ17r70eMFVbVSVedNy7uTXJLk8HxNBABYHmdyK41bktyR5KVV9WhVvXVadU2eeUozSV6Z5O7p1hr/Ncl13f3kZjYYAGA7O5Nfa157ivI3r1F2W5Lb5m8WAMBy8oQAAICBCGcAAAMRzgAABiKcAQAMRDgDABjIvE8IAE7hxoMPLKzu33rNmvd2BmAbMHMGADAQM2fAueX26xdT75X7F1MvwLNk5gwAYCDCGQDAQIQzAICBuOYMzkGL+iWoX4ECbD0zZwAAAxHOAAAGIpwBAAxEOAMAGIhwBgAwEOEMAGAgwhkAwECEMwCAgQhnAAADEc4AAAYinAEADEQ4AwAYiHAGADAQ4QwAYCDCGQDAQIQzAICBCGcAAAMRzgAABiKcAQAMZMdWNwAYx40HH9i0ui5/5Ngz3l+x+/xNqxtgOzNzBgAwEOEMAGAg657WrKqbkrwuyRPd/ZNT2XuS/FqSo9Nm7+zuz0zr9id5a5Knk/zr7v6LBbQbYHPdfv1i6r1y/2LqBbatM5k5+0iSq9Yov7G790yv48HsZUmuSXLptM8fVdV5m9VYAIDtbt1w1t2fT/LkGdZ3dZJPdPf3u/sbSR5Mctkc7QMAWCrzXHP29qq6u6puqqrnT2UvTPKtmW0encoAADgDGw1nH0iyO8meJEeS3PBsK6iqfVV1qKoOHT16dP0dAACWwIbCWXc/3t1Pd/cPknwoJ05dPpbkxTObvmgqW6uOA929t7v3rqysbKQZAADbzoZuQltVF3X3kentG5LcMy1/OsnHq+oPkvyDJJck+dLcrYQF2swbrwLAvM7kVhq3JHlVkguq6tEk707yqqrak6STPJzkbUnS3fdW1SeT3JfkqSS/3t1PL6bpAADbz7rhrLuvXaP4w6fZ/r1J3jtPowAAlpUnBAAADEQ4AwAYiHAGADCQDf1aE9g+Ln/kwFY3AYAZZs4AAAYinAEADEQ4AwAYiHAGADAQ4QwAYCDCGQDAQIQzAICBCGcAAAMRzgAABiKcAQAMRDgDABiIcAYAMBDhDABgIMIZAMBAhDMAgIEIZwAAAxHOAAAGIpwBAAxEOAMAGIhwBgAwEOEMAGAgwhkAwECEMwCAgQhnAAADEc4AAAYinAEADEQ4AwAYiHAGADAQ4QwAYCDCGQDAQIQzAICBrBvOquqmqnqiqu6ZKXt/VX2tqu6uqk9V1fOm8l1V9b2qumt6fXCRjQcA2G7OZObsI0muOqnsYJKf7O5/mOSBJPtn1j3U3Xum13Wb00wAgOWwbjjr7s8nefKksr/s7qemt19I8qIFtA0AYOlsxjVnb0nyZzPvL55OaX6uql6xCfUDACyNHfPsXFXvSvJUko9NRUeS7OzuY1X100n+pKou7e7vrLHvviT7kmTnzp3zNAMAYNvY8MxZVb05yeuS/Ep3d5J09/e7+9i0fGeSh5K8ZK39u/tAd+/t7r0rKysbbQYAwLayoXBWVVcl+fdJfrG7/+9M+UpVnTct705ySZLDm9FQAIBlsO5pzaq6JcmrklxQVY8meXdWf5353CQHqypJvjD9MvOVSX6nqv42yQ+SXNfdT65ZMQAAP2TdcNbd165R/OFTbHtbktvmbRQAwLLyhAAAgIEIZwAAAxHOAAAGIpwBAAxEOAMAGMhcTwgAzp7LHzmw1U2Yyx2Hjy2k3it2n7+QegG2ipkzAICBCGcAAAMRzgAABiKcAQAMRDgDABiIcAYAMBDhDABgIMIZAMBAhDMAgIF4QgDnjBsPPrDVTYBn7/brF1f3lfsXVzewZcycAQAMRDgDABiIcAYAMBDhDABgIMIZAMBAhDMAgIEIZwAAAxHOAAAGIpwBAAxEOAMAGIhwBgAwEOEMAGAgwhkAwECEMwCAgQhnAAADEc4AAAYinAEADEQ4AwAYyLrhrKpuqqonquqembIXVNXBqvr69Pf5M+v2V9WDVXV/Vf3cohoOALAdncnM2UeSXHVS2TuSfLa7L0ny2el9quplSa5Jcum0zx9V1Xmb1loAgG1u3XDW3Z9P8uRJxVcnuXlavjnJ62fKP9Hd3+/ubyR5MMllm9RWAIBtb6PXnF3Y3Uem5W8nuXBafmGSb81s9+hUBgDAGZj7BwHd3Un62e5XVfuq6lBVHTp69Oi8zQAA2BY2Gs4er6qLkmT6+8RU/liSF89s96Kp7Id094Hu3tvde1dWVjbYDACA7WWj4ezTSd40Lb8pyZ/OlF9TVc+tqouTXJLkS/M1EQBgeexYb4OquiXJq5JcUFWPJnl3kt9N8smqemuSbyZ5Y5J0971V9ckk9yV5Ksmvd/fTC2o7AMC2s2446+5rT7Hq1afY/r1J3jtPowAAlpUnBAAADEQ4AwAYiHAGADAQ4QwAYCDCGQDAQIQzAICBCGcAAAMRzgAABiKcAQAMRDgDABiIcAYAMBDhDABgIMIZAMBAhDMAgIEIZwAAAxHOAAAGIpwBAAxEOAMAGIhwBgAwEOEMAGAgwhkAwECEMwCAgQhnAAADEc4AAAYinAEADEQ4AwAYyI6tbgDAPO44fGwh9V6x+/yF1AuwHjNnAAADEc4AAAYinAEADEQ4AwAYiHAGADAQ4QwAYCDCGQDAQIQzAICBbPgmtFX10iS3zhTtTvLbSZ6X5NeSHJ3K39ndn9lwCwEAlsiGw1l3359kT5JU1XlJHkvyqST/MsmN3f37m9JCAIAlslmnNV+d5KHu/uYm1QcAsJQ2K5xdk+SWmfdvr6q7q+qmqnr+WjtU1b6qOlRVh44ePbrWJgAAS2fuB59X1XOS/GKS/VPRB5L8hyQ9/b0hyVtO3q+7DyQ5kCR79+7tedvBGG48+MBWNwEAzmlzh7MkP5/kK939eJIc/5skVfWhJP9tEz4DgJPdfv1i6r1y//rbAAuzGeHs2syc0qyqi7r7yPT2DUnu2YTPgHPG5Y8c2OomAHAOmyucVdWPJXlNkrfNFL+vqvZk9bTmwyetAwDgNOYKZ939f5Kcf1LZr87VIgCAJeYJAQAAAxHOAAAGIpwBAAxEOAMAGIhwBgAwEOEMAGAgwhkAwECEMwCAgQhnAAADEc4AAAYinAEADEQ4AwAYiHAGADAQ4QwAYCDCGQDAQIQzAICBCGcAAAMRzgAABiKcAQAMRDgDABiIcAYAMBDhDABgIDu2ugEADOb26xdT75X7F1MvbDNmzgAABiKcAQAMRDgDABiIcAYAMBDhDABgIMIZAMBAhDMAgIEIZwAAAxHOAAAGIpwBAAxEOAMAGMhcz9asqoeTfDfJ00me6u69VfWCJLcm2ZXk4SRv7O7/OV8zAQCWw2bMnF3Z3Xu6e+/0/h1JPtvdlyT57PQeAIAzsIjTmlcnuXlavjnJ6xfwGQAA29K84ayT/FVV3VlV+6ayC7v7yLT87SQXzvkZAABLY65rzpL8THc/VlV/P8nBqvra7Mru7qrqtXacwty+JNm5c+eczQAA2B7mmjnr7semv08k+VSSy5I8XlUXJcn094lT7Hugu/d2996VlZV5mgEAsG1seOasqn4syY9093en5Z9N8jtJPp3kTUl+d/r7p5vRUICz6Y7DxxZW9xW7z19Y3cC5b57Tmhcm+VRVHa/n493951X15SSfrKq3JvlmkjfO30wAgOWw4XDW3YeT/KM1yo8lefU8jQIAWFaeEAAAMBDhDABgIPPeSgPOWZc/cmCrmwAAP8TMGQDAQIQzAICBCGcAAAMRzgAABiKcAQAMRDgDABiIcAYAMBDhDABgIMIZAMBAhDMAgIEIZwAAAxHOAAAGIpwBAAxkx1Y3gK1x48EHtroJAMAazJwBAAxEOAMAGIhwBgAwEOEMAGAgwhkAwECEMwCAgQhnAAADEc4AAAbiJrQAZ9kdh48tpN4rdp+/kHo3ze3XL6beK/cvpl7YImbOAAAGIpwBAAxEOAMAGIhwBgAwEOEMAGAgwhkAwECEMwCAgQhnAAAD2XA4q6oXV9XtVXVfVd1bVb8xlb+nqh6rqrum12s3r7kAANvbPE8IeCrJv+3ur1TVTyS5s6oOTutu7O7fn795AADLZcPhrLuPJDkyLX+3qr6a5IWb1TAAgGW0KdecVdWuJD+V5ItT0dur6u6quqmqnr8ZnwEAsAzmDmdV9eNJbkvym939nSQfSLI7yZ6szqzdcIr99lXVoao6dPTo0XmbAQCwLcwVzqrqR7MazD7W3X+cJN39eHc/3d0/SPKhJJettW93H+juvd29d2VlZZ5mAABsG/P8WrOSfDjJV7v7D2bKL5rZ7A1J7tl48wAAlss8v9b8J0l+NclfV9VdU9k7k1xbVXuSdJKHk7xtrhYCACyReX6t+d+T1BqrPrPx5gAALDdPCAAAGIhwBgAwEOEMAGAgwhkAwECEMwCAgQhnAAADEc4AAAYinAEADEQ4AwAYiHAGADAQ4QwAYCDzPPgczorLHzmw1U2Ac8Idh48tpN4rdp+/kHqBtZk5AwAYiJmzgd148IGtbgIAcJaZOQMAGIhwBgAwEOEMAGAgwhkAwECEMwCAgQhnAAADcSsNAM5tt1+/mHqv3L+YemEdZs4AAAYinAEADEQ4AwAYiHAGADAQ4QwAYCDCGQDAQIQzAICBCGcAAANxE1o2zeWPHNjqJgBsHje3ZYsIZ5vgxoMPbHUTAIBtQjgD4LTuOHxsYXVfsfv8hdUN5yrhDADOJqdLWYcfBAAADEQ4AwAYyMLCWVVdVVX3V9WDVfWORX0OAMB2spBrzqrqvCR/mOQ1SR5N8uWq+nR337eIzwPg3LTIHxssyrA/YnAt27axqJmzy5I82N2Hu/v/JflEkqsX9FkAANvGon6t+cIk35p5/2iSf7ygzzpj7ke2ys1iAcazqFlEM33nni27lUZV7Uuyb3r7v6vq/rPwsRck+Zuz8Dmj0w8n6ItV+uEEfbFKP5yw5H3xzuML50g/vHP9Teb30kVWvqhw9liSF8+8f9FU9ne6+0CSszqFU1WHunvv2fzMEemHE/TFKv1wgr5YpR9O0Ber9MMJVXVokfUv6pqzLye5pKourqrnJLkmyacX9FkAANvGQmbOuvupqvpXSf4iyXlJburuexfxWQAA28nCrjnr7s8k+cyi6t8gV8Kv0g8n6ItV+uEEfbFKP5ygL1bphxMW2hfV3YusHwCAZ8HjmwAABnJOhbP1HglVq/7TtP7uqnr5evtW1Quq6mBVfX36+/yZdfun7e+vqp9b/Dc8Mwvqh/dX1dem7T9VVc+byndV1feq6q7p9cGz8y3PzIL64j1V9djMd37tzLplOiZunemDh6vqrql8Ox8TN1XVE1V1z0n7LNs4cap+WMZx4lR9sWzjxKn6YanGiap6cVXdXlX3VdW9VfUbM/ts3jjR3efEK6s/LHgoye4kz0nyP5K87KRtXpvkz5JUksuTfHG9fZO8L8k7puV3JPm9afll03bPTXLxtP9527gffjbJjmn592b6YVeSe7b6e5/lvnhPkn+3xuct1TFx0v43JPnt7XxMTOtemeTlJ3+/ZRon1umHpRon1umLpRknTtcPJ+2/7ceJJBclefm0/BNJHsgC8sS5NHN2Jo+EujrJR3vVF5I8r6ouWmffq5PcPC3fnOT1M+Wf6O7vd/c3kjw41bPVFtIP3f2X3f3UtP8XsnpvutEt6pg4laU6Jo6rqkryxiS3LPqLbIJ5+iLd/fkkT65R7zKNE6fshyUcJ053TJzKUh0Txy3LONHdR7r7K0nS3d9N8tWsPhXp+D6bMk6cS+FsrUdCvfAMtzndvhd295Fp+dtJLnwWn7cVFtUPs96S1f8xHHfxNC39uap6xUYbvgCL7Iu3T1PZN81MTS/rMfGKJI9399dnyrbjMXE6yzROnKllGCfWsyzjxJlYunGiqnYl+akkX5yKNm2cOJfC2cL16vzjUv98tareleSpJB+bio4k2dnde5L8myQfr6q/t1XtO0s+kNXp7j1Z/f43bG1ztty1eeb/hpfxmPg7xgnjxMQ48UxLNU5U1Y8nuS3Jb3b3d05eP+84cS6Fs3UfCXWabU637+PHp22nv088i8/bCovqh1TVm5O8LsmvTAdWpmnYY9PynVk9V/6Szfgim2AhfdHdj3f30939gyQfyonp52U8JnYk+aUktx4v28bHxOks0zhxWks2TpzSko0Tp7Vs40RV/WhWg9nHuvuPZ7bZvHGiB7g470xeWb1h7uGsXkx3/AK+S0/a5p/nmRfwfWm9fZO8P8+8gO990/KleeYFfIczxkWdi+qHq5Lcl2TlpLpWjn/vrP4v8bEkL9jqflhwX1w0s/9vZfVagaU7JmaOi88twzExs35Xfvji76UZJ9bph6UaJ9bpi6UZJ07XDzPHxVKME9P7jyb5j2vUu2njxJZ30rPs0Ndm9ZcRDyV511R2XZLrZjrtD6f1f51k7+n2ncrPT/LZJF9P8lezB0+Sd03b35/k57f6+y+4Hx7M6jnxu6bXB6fyf5Hk3qnsK0l+Yau//1noi/8ybXt3Vp8JOzsIL80xMa37yPE6Zsq28zFxS1ZPx/xtVq8LeetUvmzjxKn6YRnHiVP1xbKNE2v2w7RuacaJJD+T1dOVd8/8O3jttG7TxglPCAAAGMi5dM0ZAMC2J5wBAAxEOAMAGIhwBgAwEOEMAGAgwhkAwECEMwCAgQhnAAAD+f+ywW/Mq9C71wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6594f0be50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax1 = plt.subplots(1,1)\n",
    "ax1.hist(S[H == 0], 20, normed=1, alpha = 0.5);\n",
    "ax1.hist(S[H == 1], 20, normed=1, alpha = 0.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Is offset robust to changes in firing rate?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n"
     ]
    }
   ],
   "source": [
    "## Is the offset robust to changes in firing rate?\n",
    "N = 20\n",
    "wvals = np.linspace(2, 20, N)\n",
    "\n",
    "c = 0.99\n",
    "diffsc_0_99 = np.zeros(N)\n",
    "ratesc_0_99 = np.zeros(N)\n",
    "\n",
    "lif.params.c = c\n",
    "\n",
    "for idx, wval in enumerate(wvals):\n",
    "    print idx\n",
    "    lif.W[0] = wval\n",
    "    (v, h, _, _, u) = lif.simulate()\n",
    "    s1 = np.convolve(h[0,:], exp_filter)[0:len_h]\n",
    "    s_r = s1[:].reshape((-1, DeltaT))\n",
    "    S = np.squeeze(s_r[:,-1])\n",
    "    H_r = h[0,:].reshape((-1, DeltaT))\n",
    "    H = np.sum(H_r, 1)\n",
    "    diffsc_0_99[idx] = np.mean(S[H == 1]) - np.mean(S[H == 0])\n",
    "    ratesc_0_99[idx] = np.sum(h)\n",
    "    \n",
    "diffsc_0_01 = np.zeros(N)\n",
    "ratesc_0_01 = np.zeros(N)\n",
    "\n",
    "c = 0.01\n",
    "lif.params.c = c\n",
    "\n",
    "for idx, wval in enumerate(wvals):\n",
    "    print idx\n",
    "    lif.W[0] = wval\n",
    "    (v, h, _, _, u) = lif.simulate()\n",
    "    s1 = np.convolve(h[0,:], exp_filter)[0:len_h]\n",
    "    s_r = s1[:].reshape((-1, DeltaT))\n",
    "    S = np.squeeze(s_r[:,-1])\n",
    "    H_r = h[0,:].reshape((-1, DeltaT))\n",
    "    H = np.sum(H_r, 1)\n",
    "    diffsc_0_01[idx] = np.mean(S[H == 1]) - np.mean(S[H == 0])\n",
    "    ratesc_0_01[idx] = np.sum(h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7f6594e13f90>"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAAGFCAYAAABNKN+1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcXnV99//XZyYbZCXJZCELJCQsAcI2QsqiKFUWl4jSitUbgbaUKvbuvdii1taftupdf/fdlqqg7Y+qva3UDY2CouCOICQsWQgJkwTIvkwgeyaZmc/vj+skTIYkcyWZZM6VvJ6Px3nMdZ3zPed8z3euXPPO93uWyEwkSZJUPnU9XQFJkiTtnUFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJR4mI+HRE/HmVZT8eER/fy/ylEXF+t1du3/V4LCLOPFL7K/Z5SO10pNuo2OcRbydJ5WBQk44CEdEA3AB8sdP890TErIh4OSJeioifR8SYfWzjBOBEYH431+22iJgZES0R8eVOi/9f4BPdub8u6vKqdipJGw2NiHsjYktEvBARf9CpyBFtp65UUd8DKt/FZ0Q6pvXq6QpIql5EfA04CXi+w+wNxfv7M3Nbh7I3An8FvAt4AhgGvB14aR+bPxt4vuM2uskK4G+BK4HjOi2bAdwVEaMyc1V37GxfbZSZHwBupEM7laiNPg/sAEYC5wL3RcTTmTmvWN7t7XSIuqrvgZbf32dEOqbZoybViIgYBjwCfC0z39th+gBwNfCLTqv8IfDFzJyVFesy818zc+s+djEVaIqIOyJibUSsiIg3Hmq9M/M7mfldoHkvy7YDs6j8gT5kXbQRvLqderyNIqI/8E7gY5m5OTN/DXwP+C+7yhxoO0XEuIj4TlHH5oj43KHU8UDre6Dl9/cZkY51BjWpdlwB/ATYW4g4G1jQad424OaI+P2IGF7F9s8GGoH7qPR8fBH4y44FIuIHxRDh3qYfHOgBFeYD5xzkup3tr43g1e3U7W0EB9xOpwKtmbmww7yngc7npFXVThFRD/wAeAE4GRgD3LOf8gf6O622vgdbXlIHDn1KtWMisBB4b4denO9l5jeBIcCmTuVvAG6ncn7TmIi4H/jDzFyzj+1PBf4uMx8AiIhngMs6FsjMt3TLkexpEzC6m7a1vzaCV7dTt7cRHHA7DQA2dpq3ERjYaV617XQhlfPoPpSZrcW8X++r8EH8Tqut78GWl9SBPWpS7WgFLgKu7DCktyuAvESnP3yZuSoz/zwzx1P54z2VvfT+AEREAGcB3+8w+yzgmW4+hr0ZCLy8j3r9PCJyH9Pewsf+2gg6tVNJ2mgzMKjTvMG8Onjvs506GQe80CGkdbdq63uw5SV1YFCTakBE7Oop6pOZ7XspMpvKENNeZeYsYA7Qfx9FJhQ/mzrMOw94qlM9fhgRm/cx/bDa4+nkDCpDYXur9+WZGfuYLu1Ut67aCPbTTt3VRkVdDqSdFgK9ImJyh3nnAJ1PzN9nO3WyFBgfEVWNmBzE77Ta+h5seUkdGNSk2nAZ8BiVc6r25n7gdbveRMTtEXFJRPQtphuBy4F/28f6U4E5mZkd5p1Hp2CQmVdn5oB9TFfvbcMR0Ssi+gH1QH1E9NsVIor5F1A5r+xQddVG0KGdDlcbwYG1U2ZuAb4DfCIi+kfEpcDbgH/fVeYA2+kxYCXwmWJ7/SLikn0VPtDfaTX1PYjj2+dnRDrW+Q9Bqg2XANcDdRFxeYf5d2fmT4GvAk9FxHHFrSMGUQkcJwJbqNx64orM/O0+tn82HQJHcWL9KGBuN9T9r4C/6fD+vcD/A3wceCvw88xc0Q376aqNoEM7Ua42ej9wN7CGypWPf9rpVhevaqeit+tXmfmpjhvKzLaIeCtwB/AikMB/AA93Qz2rqu9e6tbV8e3vMyId02LP/xxKqlUR8SlgTWb+YxVlPw6QmR8/zNXqqh6/pXLyfneEnWr3aTtJqhn2qElHicz8SE/X4UBl5kU9sE/bSVLNqOoctYi4KiIWRERTRNy+l+URlRtANkXE7OjwHLx9rRsRnyzKPhURP46IEzss+3BRfkFEdMuNMCXt4efFpP37ObaTpB7U5dBncfPEhcAbgWXA48C7M/OZDmWuAT4IXEPl0vh/ysyL9rduRAzKzI3F+n8GTMnMWyNiCvB1XrkX0IPAqZnZ1o3HLUmSVHrV9KhdCDRl5uLM3EHlDtfTO5WZDny1eATLo8CQiBi9v3V3hbRCfyonvO7a1j2Z2ZKZS6hcCn/hQR6fJElSzaomqI2hcl+eXZYV86ops991I+LvImIp8B7grw9gf5IkSUe9Hr2YIDM/Cnw0Ij4M3Mael2fvV0TcAtwC0L9//wtOP/30w1NJSZKkbjRr1qx1mdlQTdlqgtpyKo8k2WVsMa+aMr2rWBfga1RuRPk3Ve6PzPwS8CWAxsbGnDlzZhWHIkmS1LMi4oVqy1Yz9Pk4MDkiJkREHyo3lJzRqcwM4Ibi6s9pwIbMXLm/dTs9TmQ68GyHbV1f3Cl8AjCZyp22JUmSjild9qhlZmtE3AY8QOXxHndn5ryIuLVYfheV3rBrqJz4vxW4aX/rFpv+TEScBrQDLwC7tjcvIr5B5UHHrcAHvOJTkiQdi46KJxM49ClJkmpFRMzKzMZqyvpQdkmSpJIyqEmSJJWUQU2SJKmkDGqSJEklZVCTJEkqKYOaJElSSRnUJEmSSsqgJkmSVFIGNUmSpJIyqEmSJJWUQU2SJKmkDGqSJEklZVCTJEkqKYOaJElSSRnUJEmSSsqgJkmSVFIGNUmSpJIyqEmSJJWUQU2SJKmkDGqSJEklZVCTJEkqKYOaJElSSRnUJEmSSsqgJkmSVFIGNUmSpJIyqEmSJJWUQU2SJKmkDGqSJEklZVCTJEkqKYOaJElSSRnUJEmSSsqgJkmSVFIGNUmSpJIyqEmSJJWUQU2SJKmkDGqSJEklZVCTJEkqKYOaJElSSRnUJEmSSsqgJkmSVFIGNUmSpJIyqEmSJJWUQU2SJKmkDGqSJEklZVCTJEkqqaqCWkRcFRELIqIpIm7fy/KIiDuK5bMj4vyu1o2Iz0bEs0X5eyNiSDH/5IjYFhFPFdNd3XGgkiRJtabLoBYR9cDngauBKcC7I2JKp2JXA5OL6RbgzirW/QlwVmZOBRYCH+6wvUWZeW4x3XqwBydJklTLqulRuxBoyszFmbkDuAeY3qnMdOCrWfEoMCQiRu9v3cz8cWa2Fus/CozthuORJEk6alQT1MYASzu8X1bMq6ZMNesC3Az8sMP7CcWw5y8i4rIq6ihJknTU6dXTFYiIjwKtwNeKWSuB8ZnZHBEXAN+NiDMzc2On9W6hMszK+PHjj2SVJUmSjohqetSWA+M6vB9bzKumzH7XjYgbgbcA78nMBMjMlsxsLl7PAhYBp3auVGZ+KTMbM7OxoaGhisOQJEmqLdUEtceByRExISL6ANcDMzqVmQHcUFz9OQ3YkJkr97duRFwF/AXwtszcumtDEdFQXIRAREykcoHC4kM6SkmSpBrU5dBnZrZGxG3AA0A9cHdmzouIW4vldwH3A9cATcBW4Kb9rVts+nNAX+AnEQHwaHGF52uBT0TETqAduDUz13fXAUuSJNWKKEYca1pjY2POnDmzp6shSZJqWGZSdB4dVhExKzMbqynrkwkkSdIxLzP52Pfm8r9/vIAydWIZ1CRJ0jHvCz9fxP999EV2th2ZXrVqGdQkSdIx7Zszl/LZBxZw7Xlj+IsrT+vp6uzBoCZJko5ZP1+whtu/M4dLJw3nf71zKnV15elNA4OaJEk6Rs1ZtoH3f+0JThs5kDvfez59epUvFpWvRpIkSYfZi81buenLj3HC8X348k2vYWC/3j1dpb3q8UdISZIkHUnNm1t43789Rmt7cs/NFzJiUL+ertI+2aMmSZKOGVt3tHLzV2ay4uVt/H/va2TSiAE9XaX9MqhJkqRjQmtbOx/8jyeZs+xl7nj3eVxw0tCerlKXHPqUJElHvV03tH3o2TX87dvP4sozR/V0lapij5okSTrq3fFQE19/bCm3vX4S7512Uk9Xp2oGNUmSdFT7z8df5B8eXMg7zx/L/3jTqT1dnQNiUJMkSUetnz67mo/cO5fXntrAZ955dqkeD1UNg5pUIus2t/DAvFWs3LCtp6siSTXvqaUv84GvPcmU0YO48z3n07u+9mKPFxNIJdDenvznzKV8+v75bNzeCsBJw47ndyYOY1oxjRpc3vv8SFLZPL9uCzd/+XGGD+zD3Te+hv59azPy1GatpaPIc6s38ZF75/D48y9x0YShfOD1k1i4ehOPLl7PfXNWcs/jSwGYMLw/0yYOZdrEYfzOxGGlvkGjdLhkJi9t3cmaTdtZs7GFNZtaWLNpO+s37+Dc8UP43TNG0q93fU9XUz1s7aYWbrj7MQC+evNFNAzs28M1OniRmT1dh0PW2NiYM2fO7OlqSAdk+842Pv+zJu76xSL69+3FR685g+suGLvH+RNt7cn8lRt5ZFEzjy5u5rEl69nUUulxm9jQf3dv27SJQxkxsLzBLTOZs3wDM55awW8WNZNAfR3U19VRH9Crro66ul0/g151QV0E9XuZ16sudr+v3zWvfu/L6uuC+njl9e5lHebtmkYN7se5Y4eU7oHMx4rWtnaat+wowtf2SgDr+HpTC2s3bmft5hZ2tr3671bv+mBnWzKoXy/edu6JXHfBOM4ZO7jmzkfSodvS0sq7/+VRFq7exNf/eBrnjT+hp6v0KhExKzMbqyprUJOOvIeb1vHRe+fwfPNW3nHeGD765jMYNqDr//G1tSfzVmzg0cXNPLp4PY8tWc/mIrid0tCf3zmlEtwumjCsFP+DfG71JmY8vYLvP72C55u30rs+uGjCMI7rU09be+45Zb563j7mt7Yn7Zm0trXTntDa3k57e/HzEL7SRgzsy5vOHMlVZ47moolDa/J8lrJpaW3b3fO1dm8BrFi2fkvLXn93Q47vzYiBfRkxsB8jBvalYdArr0cM7MuIQZXX/XrX85tF6/jWrGX8aO4qWlrbmTRiANddMJZrzxvDSHugjwk729r546/O5JcL1/IvNzRyxRkje7pKe2VQOwzu/vUSzh0/hKljBtPLL28dpObNLfzdffP5zpPLOXnY8fzdtWdzyaThB7291rZ25q3YyCOLKz1ujy9Zz5YdbQBMHjGgMkx6yjAumjC0qiDYHZau38r3Z69gxlMreHbVJuoCLj5lOG8750SuPHMUg48/vA8+zsxXh7fi577CYGt7snD1Jn40dxU/X7CWbTvbGHxcb373jJFcddYoLps83OG0Tja3tLJm4yu9XWs2bmftpleGIncFsA3bdr5q3bqA4QP6MqJT6GoY1I+G3fP70jCwL317HXi7b9y+k/tmr+Tbs5Yx84WXqAu4bHID110wljdOcWj0aJWZ/MW3ZvPNWcv49DvO5t0Xju/pKu2TQa2brdqwnWmffgiAgX17cdHEYbxpykh+r3Gs3eqqSmbyzVnL+NT989nS0sqtrzuFD7x+Urf/wWhta2fO8g08ung9jyxuZubz69laBLfTRg7cfY7bRROHMbR/n27b75pN27lv9kpmPL2CJ198GYALTjqBt51zIlefParUw7KdbdvRxi+fW8sDc1fx4PzVbNzeyvF96nn9aSO48qxRvP60Bgb2O7xhswxe2rKDBas38dzqTSxZt5XVm7aztkNP2K7PVUd96utoGPhK0NodwopAtmvZsP59qT9CQ8yL127mO08s59tPLGPlhu0M6teLt55zItddMJZzxw2p2e/wjdt38uiiZh5Z3Eyf+jqmjh3C1LGDGXvCcTV7TIfq//x4AXf8tIn/esVk/tsby32vNIPaYdC8uYVHFjfzcFMzv25ay9L127j5kgl87C1nHLP/KFSdRWs385HvzOG3S9bTeNIJfPodZzN55MAjsu+dRXDbdY7bzOdfYtvOyh/Y00cN3H2O20UThnLCAQa3DVt38sO5K/n+7BU8sqiZ9oQzRg/ibeecyFumjmbc0OMPxyEdUTta23l0cTM/mreKH89bzbrNLfSpr+PSycO56sxR/O6Ukd0aeHvCxu07eW71Jhau3syCVZt4bs0mFqzazLrNLbvLHN+nnpGDiqC1K4ANevXrwcf1Lu33YVt78siiZr41ayk/mreK7TvbOaWhP9ddMI5rzxtT+quqW1rbeOKFl/nNonX8umkdTy99mfaE43pXTiPY0dYOwLD+fTh77GDOGTuEc8YNZurYIQw/Qr3pPelrv32Bj947l3c1jquJe6UZ1A6zzOQTP3iGf3v4ef74sgl85BrDml6tpbWNL/xsEXf+fBH9etfx4WvO4F2N43r0ZPUdre3MWf5ypcdtUTMzX1jP9p3tRMDpowYxbeJQfqc4x21vQ5RbWlp5cP5qvv/0Cn6xcC0725IJw/vz1nNO5G3njGbSiCMTQHtCW3vy5Isv8aO5q/jRvFUse2kbdQEXThjK1WeN5k1njmT04ON6upr7tHVHK8+t3ry7l2zh6s0sXL2JlRu27y5zfJ96Jo8YwOSRAzlt5EAmjxzAaaMGMmpQv6PqO27j9p3cP3sl335iGY8/X86h0fb2ZP6qjTzctI5fNzXz2JJmtu9sp74uOHfcEC6ZNJxLThnGeeNPIEkWrNrE08s2MHvpy8xetoHn1mzafc7fmCHHMXVsJbSdM24wZ48ZfFT1Cv/kmdX8yb/P5HWnNvAvNzTWxOlJBrUjIDP5+Ix5fOWRF/iT103k9qtOP6q+yHRoHlnUzEfvncPidVuYfu6J/NWbp5Ti5P7OdrS2M3vZy5UetyWVHreW1kpwmzJ60O4et8xkxtMreGj+GrbtbGP04H68Zepo3nbOGM4aM+iY++xnJvNWbOSBeav40dxVPLdmMwDnjBvCVWeO4qqzRjFheP8eqdv2nW0sWlsJYQtXb2bhqk0sXLOJpetfuYlyn151TGqohLDJIwdw2siBnDpyIGOGHHfMXfW6ZN0WvvPEMr49axkrNmxnYIeh0fOO8NDo0vVb+XVTpcfskUXNrN+yA6icb3rJpOFcOmk4F00cWlXI2tLSytzlG5i9bANPL6uEtxfXbwUgAiYO7885xXDp1HFDmDJ6UCkC6oGa9cJLvOdfH+W0kQP5+i3TOL5Pbdx1zKB2hGQmH/veXP7voy/y/stP4UNXnnbM/cHSnl7asoNP3T+fb85axrihx/G3bz+b153a0NPVqlpLaxtPL61cVfrIomZmvfgSO1orQypD+/fhmrNH8bZzxtB40gnH3B/0/Wlas5kH5q3igXmrmL1sA1A5J/DKs0Zx1ZmjOGP0wG7/btjR2s6SdVuKQLZpdzB7oXnL7p6U3vXBxOEDdoexySMHctqogYwfevwRO0esVrS3J48sbuZbs5bxw7kr2b6znYkN/bnugrG847yxh2VodP2WHfxm0Toebmrm4aZ1u4PUyEF9uXRSA5dMGsYlk4Z32xWr67fsYHYR2mYve5mnl21g7abKEHevuuC0UQM5Z9wQzil63yaPGFDq3qlFazdz3Z2/YfBxvfnWn15cU0O8BrUjqL09+eh35/L1x17kg2+YxH9/46mGtWNQZnLvk8v52/vms3HbTv74tRP5szdM5rg+tfc/1I6272zjqaUv09qW3q6iSstf3sYDxfDo48+vJxPGDz2eq84axZVnjuK8cQd2r7bWtnZeWL+10jO2ejML12xi4apNLFm3hdYikdXXBScNO/6VMDZyIKeOHMDJw/v7OzsIm7bv5IdzVvGtWct47Pn11AVcOrmBd54/hivPHHXQPU/bdrTx+PPri+HMdcxbsRGoXKQ27ZRhXDppOJdMGs4pDf2PyN+RzGTVxu08vXRXcKuEuE3F01GO613PmScO2j1kes7YIZw07PhS/I1bs2k77/jCb9i2o43vvP9iThrWMz3YB8ugdoS1tycfuXcO9zy+tCauNlH3WrJuC3/13Tk83NTMeeOH8Ol3nM3powb1dLVUAms3tfDg/NX8aO4qfrNoHTvbkhED+3JlMTx64YRXwm97e7LspW0s6NRDtmjt5t29mhGV0Dd5xEBOGzWAU0cOZPKIgUxs6F+Tw1a14PldQ6NPLGf5y9sY2K8Xb5laGRo9f/z+h0Z3XYW9K5g98cLL7Ghrp3d9cMFJJ+wOZmeX6LZP7e3J881b9hgynbt8Ay3FZ3Dwcb2L892Kc97GDjniF2JsbmnlXV98hCXrtnDPLdOYOnbIEd1/dzCo9YD29uQvv125f8t/f+Op/NkVk3u0Pjr8drS288VfLOKff9ZE3/o6/uLq03nPheMdEtRebdy+k589u2aPe7UNOb43F548lFUbt/Pc6s27r8iFygngp44swljRSzZpxICa76WtVe3tyaNLiqHROavYtrONicP7884LxvKO88cwevBxZCaL1m7h4aZ1PNy0jkcWN+/unZoyehCXTq4Es9ecfELNnEsFlavHF67e9MqQ6dINLFi9ibaiR3fEwL5MHTuEc4urTKeOHcyQ4w/P1dA729q5+cuP85tFzfzr+xp5/WkjDst+DjeDWg9pa08+9K2n+c4Ty/nQlafxgddP6ukq6TBYu6mFh5vW8bmfNdG0ZjNvPns0f/PWKT57U1XreK+2J158iXEdeskmjxzI5BEDjqqr8o42m1tauX/OysrQ6JL1REDjSSewdP02Vm2sXEU79oTjuGzycC4+ZTgXnzLsiN1w+kjZvrONeSs28vTSl3ef97Z43Zbdy08adnzR41YJb2eNGXTI4TQz+R/frPyN/ft3TuX3XzPuUA+jxxjUelBbe/I/v/k09z65nL+86nT+9PJTerpKOkTbdrTx2PPr+fVza/nVc+t4dtUmoNLj8cm3n8kbTi/nI0okHX4vNG/h208s58FnVjNheP/dV2eOH1b79xE8UBu27WTu8mLIdGnl567bv9QFnDpy4B5DpqeNGkifXtUP+X72gWf5/M8WHRWjVga1HtbWnvy3/3yKGU+v4CPXnM4trzWs1ZJdz9P81XOV4YuZz7/EjrZ2+tTXVc4rmVz5Ij5rzGCvnJOk/VizaTuzl75ylensZS/z0tbKY8X69KrjjNGDOLfDPd4mDh+w19NH/v2R5/nY9+bx7gvH86lrzyrFBQ2HwqBWAq1t7fz5fz7FD2av5K/efAZ/dNnEbt9HZvLrpnX804PPseLlbbz21AauOGMkl04a7nksB2j3/YueW8fDi9bxcvFFcvqogVw6aTiXTh7OhROG1tR5JZJUNpnJ0vXbigsVKuFt7vINux9JNqBvL84aM6i4x1slvM1dvpE//dosrjh9BHe994LSXHhxKAxqJdHa1s6f3fMk989ZxV+/ZQo3Xzqh27b9m0Xr+IefLOTx51/ixMH9OHvsYB5uamZzSyt9e9Vx8SnDuOKMkbzh9BGcOKS8d0vvKRu27uSRxev41XOVq7FeaN7z/kWXTR7OxZOG1dQzKiWpFrW1J4vWbubppa/cImT+yo3sbHsln5w7bghf/+NpR00nxIEENbsHDqNe9XX80/Xn0d7+JJ/4wTPU1wXvu/jkQ9rmzOfX879/vJBHFjczclBfPjn9TH7/NePo26ueHa3tPLZkPQ89u5qH5q/hZwvmApWrja44YwRXnDGSqWMGH5NXJe5obeeJF1/i18+t41dN65izrPKcvP596pk2cRg3Xnwyl00ezikNA2q+S12Sakl9XXBq8XSM32usXCDQ0trGsys3MXvZy6zd1MKNl0w4akLagbJH7QjY0drObf/xBD9+ZjUfuvI03n/5KQccBp588SX+z08W8qvn1jF8QF/ef/kp/MFF4/d576TKZeKbeXD+Gn46fw0zX1hPe8LwAX15w+kNvOH0kVw2eTj9+x6dWT0zWbh6M796bi0PN63jt0vWs3VHG/V1wTljB3Pp5Eqv2bnjhnhDUEnSEeXQZwntbGvnQ998mu8+tYKbLjmZj715SlU9W3OWbeAfHlzIT59dw9D+fbj1dRP5L9NOPuD/Wby0ZQe/WLiWB+ev5hcL17Jpeyt96uuYdsowfveMEbzh9BGMPeHwX6WUmexoa2f7znZadraxfWc723a2sX3X1NrOth1ttLTumtdxefsr5Xa9b21j247Kei3F/G0729jS0sbmlsr9iyYO77/7AoBppwxjkLc9kCT1IINaSbW3J39733zufngJ0889kc9ed84+L02ev3Ij//CThfz4mdUMPq43t7x2IjdefHK39IDtbGvn8efX89D8NTw0fzXPF+dnnT5qIG84fQRjTjiO1raktT1pbWsvfiat7e17mVd539aetLR2Cl0dg1XrK6/bD/Ij16e+jr696ziudz39etfTr3dd8bOYetXtnl959MlgLpk8nDGeoydJKhGDWollJnf+YhF//6MFvPbUBu567/l7XEm4cPUm/unB57hvzkoG9uvFH106kZsvPfmw3vxy0drN/HT+Gh6cv5qZL7y0+27Te9O7PuhVV0evuqBXfVBfV0fv+qC+LuhbBKWOQarv7vd19OtVmX9cn/rdZft1XL6P8FUpX++tMCRJRwWDWg2457EX+ci9c5g6dgj/duNrWL91B3c89Bwznl7B8b3rufnSCfzRpRMZfPyRHabbtH0nW3e0VYJYXR296iuBrFddnUFJkqRu4FWfNeD6C8dzQv8+fPDrT/Kmf/wlzZtb6Nurnj957Snc8tqJDO1/eJ6T1pWB/Xr76BpJkkrCoNaDrjxzFF+9+UI+cu8cpp9zIn/yulNoGHh0PQ9OkiQdPINaD5s2cRg//R+X93Q1JElSCXkDKUmSpJIyqEmSJJWUQU2SJKmkDGqSJEklVVVQi4irImJBRDRFxO17WR4RcUexfHZEnN/VuhHx2Yh4tih/b0QM6bDsw0X5BRFx5aEepCRJUi3qMqhFRD3weeBqYArw7oiY0qnY1cDkYroFuLOKdX8CnJWZU4GFwIeLdaYA1wNnAlcBXyi2I0mSdEyppkftQqApMxdn5g7gHmB6pzLTga9mxaPAkIgYvb91M/PHmdlarP8oMLbDtu7JzJbMXAI0FduRJEk6plQT1MYASzu8X1bMq6ZMNesC3Azb/KH4AAASnklEQVT88AD2J0mSdNTr8YsJIuKjQCvwtQNc75aImBkRM9euXXt4KidJktSDqglqy4FxHd6PLeZVU2a/60bEjcBbgPfkK0+Hr2Z/ZOaXMrMxMxsbGhqqOAxJkqTaUk1QexyYHBETIqIPlRP9Z3QqMwO4obj6cxqwITNX7m/diLgK+AvgbZm5tdO2ro+IvhExgcoFCo8dwjFKkiTVpC6f9ZmZrRFxG/AAUA/cnZnzIuLWYvldwP3ANVRO/N8K3LS/dYtNfw7oC/wkIgAezcxbi21/A3iGypDoBzKzrduOWJIkqUbEKyOOtauxsTFnzpzZ09WQJEnqUkTMyszGasr2+MUEkiRJ2juDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVVFVBLSKuiogFEdEUEbfvZXlExB3F8tkRcX5X60bE70XEvIhoj4jGDvNPjohtEfFUMd11qAcpSZJUi3p1VSAi6oHPA28ElgGPR8SMzHymQ7GrgcnFdBFwJ3BRF+vOBd4BfHEvu12Umece/GFJkiTVvmp61C4EmjJzcWbuAO4BpncqMx34alY8CgyJiNH7Wzcz52fmgm47EkmSpKNMNUFtDLC0w/tlxbxqylSz7t5MKIY9fxERl1VRXpIk6ajT5dBnD1gJjM/M5oi4APhuRJyZmRs7FoqIW4BbAMaPH98D1ZQkSTq8qulRWw6M6/B+bDGvmjLVrLuHzGzJzObi9SxgEXDqXsp9KTMbM7OxoaGhisOQJEmqLdUEtceByRExISL6ANcDMzqVmQHcUFz9OQ3YkJkrq1x3DxHRUFyEQERMpHKBwuIDOipJkqSjQJdDn5nZGhG3AQ8A9cDdmTkvIm4tlt8F3A9cAzQBW4Gb9rcuQERcC/wz0ADcFxFPZeaVwGuBT0TETqAduDUz13fnQUuSJNWCyMyersMha2xszJkzZ/Z0NSRJkroUEbMys7Hrkj6ZQJIkqbQMapIkSSVlUJMkSSopg5okSVJJGdQkSZJKyqAmSZJUUgY1SZKkkjKoSZIklZRBTZIkqaQMapIkSSVlUJMkSSopg5okSVJJGdQkSZJKyqAmSZJUUgY1SZKkkjKoSZIklZRBTZIkqaQMapIkSSVlUJMkSSopg5okSVJJGdQkSZJKyqAmSZJUUgY1SZKkkjKoSZIklZRBTZIkqaQMapIkSSVlUJMkSSopg5okSVJJGdQkSZJKyqAmSZJUUgY1SZKkkjKoSZIklZRBTZIkqaQMapIkSSVlUJMkSSopg5okSVJJGdQkSZJKyqAmSZJUUgY1SZKkkjKoSZIklZRBTZIkqaQMapIkSSVlUJMkSSopg5okSVJJGdQkSZJKqqqgFhFXRcSCiGiKiNv3sjwi4o5i+eyIOL+rdSPi9yJiXkS0R0Rjp+19uCi/ICKuPJQDlCRJqlVdBrWIqAc+D1wNTAHeHRFTOhW7GphcTLcAd1ax7lzgHcAvO+1vCnA9cCZwFfCFYjuSJEnHlGp61C4EmjJzcWbuAO4BpncqMx34alY8CgyJiNH7Wzcz52fmgr3sbzpwT2a2ZOYSoKnYjiRJ0jGlmqA2Blja4f2yYl41ZapZ92D2J0mSdNSr2YsJIuKWiJgZETPXrl3b09WRJEnqdtUEteXAuA7vxxbzqilTzboHsz8y80uZ2ZiZjQ0NDV1sUpIkqfZUE9QeByZHxISI6EPlRP8ZncrMAG4orv6cBmzIzJVVrtvZDOD6iOgbEROoXKDw2AEckyRJ0lGhV1cFMrM1Im4DHgDqgbszc15E3Fosvwu4H7iGyon/W4Gb9rcuQERcC/wz0ADcFxFPZeaVxba/ATwDtAIfyMy2bj1qSZKkGhCZ2dN1OGSNjY05c+bMnq6GJElSlyJiVmY2dl2yhi8mkCRJOtoZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVVFVBLSKuiogFEdEUEbfvZXlExB3F8tkRcX5X60bE0Ij4SUQ8V/w8oZh/ckRsi4iniumu7jhQSZKkWtNlUIuIeuDzwNXAFODdETGlU7GrgcnFdAtwZxXr3g48lJmTgYeK97ssysxzi+nWgz04SZKkWlZNj9qFQFNmLs7MHcA9wPROZaYDX82KR4EhETG6i3WnA18pXn8FePshHoskSdJRpZqgNgZY2uH9smJeNWX2t+7IzFxZvF4FjOxQbkIx7PmLiLisijpKkiQddXr1dAUAMjMjIou3K4HxmdkcERcA342IMzNzY8d1IuIWKsOsjB8//shWWJIk6QiopkdtOTCuw/uxxbxqyuxv3dXF8CjFzzUAmdmSmc3F61nAIuDUzpXKzC9lZmNmNjY0NFRxGJIkSbWlmqD2ODA5IiZERB/gemBGpzIzgBuKqz+nARuKYc39rTsDeF/x+n3A9wAioqG4CIGImEjlAoXFB32EkiRJNarLoc/MbI2I24AHgHrg7sycFxG3FsvvAu4HrgGagK3ATftbt9j0Z4BvRMQfAi8Av1/Mfy3wiYjYCbQDt2bm+m45WkmSpBoSmdl1qZJrbGzMmTNn9nQ1JEmSuhQRszKzsZqyPplAkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSVlEFNkiSppAxqkiRJJWVQkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSVQW1iLgqIhZERFNE3L6X5RERdxTLZ0fE+V2tGxFDI+InEfFc8fOEDss+XJRfEBFXHupBSpIk1aIug1pE1AOfB64GpgDvjogpnYpdDUwupluAO6tY93bgocycDDxUvKdYfj1wJnAV8IViO5IkSceUanrULgSaMnNxZu4A7gGmdyozHfhqVjwKDImI0V2sOx34SvH6K8DbO8y/JzNbMnMJ0FRsR5Ik6ZhSTVAbAyzt8H5ZMa+aMvtbd2RmrixerwJGHsD+JEmSjnq9eroCAJmZEZEHsk5E3EJlmBVgc0Qs6P6a9bjhwLqersRRwHbsPrZl97Adu4ft2H1sy+5RbTueVO0Gqwlqy4FxHd6PLeZVU6b3ftZdHRGjM3NlMUy65gD2R2Z+CfhSFfWvWRExMzMbe7oetc527D62ZfewHbuH7dh9bMvucTjasZqhz8eByRExISL6UDnRf0anMjOAG4qrP6cBG4phzf2tOwN4X/H6fcD3Osy/PiL6RsQEKhcoPHaQxydJklSzuuxRy8zWiLgNeACoB+7OzHkRcWux/C7gfuAaKif+bwVu2t+6xaY/A3wjIv4QeAH4/WKdeRHxDeAZoBX4QGa2ddcBS5Ik1YrIPKBTw3QERcQtxRCvDoHt2H1sy+5hO3YP27H72Jbd43C0o0FNkiSppHyElCRJUkkZ1HpYFY/nml48luupiJgZEZf2RD1rQVdt2aHcayKiNSKuO5L1qxVVfCYvj4gNxWfyqYj4656oZ9lV83ks2vKpiJgXEb840nWsFVV8Jj/U4fM4NyLaImJoT9S1zKpox8ER8f2IeLr4TN7UE/WsBVW05QkRcW/x9/uxiDjroHeWmU49NFG5wGIRMBHoAzwNTOlUZgCvDFFPBZ7t6XqXcaqmLTuU+ymVC2Cu6+l6l22q8jN5OfCDnq5rmacq23EIlYumxhfvR/R0vcs4Vftvu0P5twI/7el6l22q8jP5EeB/Fa8bgPVAn56ue9mmKtvys8DfFK9Pp/LIzIPanz1qPavLx3Nl5uYsftNAf8CTCveumkedAXwQ+Dav3LdPe6q2HbV/1bTjHwDfycwXATLTz+TeHehn8t3A149IzWpLNe2YwMCICCqdBOup3H1Be6qmLadQ6RQgM58FTo6IkRwEg1rPqupxWRFxbUQ8C9wH3HyE6lZrumzLiBgDXAvceQTrVWuqfYTbxUWX/g8j4swjU7WaUk07ngqcEBE/j4hZEXHDEatdban6sYIRcTxwFZX/jGlP1bTj54AzgBXAHOC/Zmb7kaleTammLZ8G3gEQERdSeRLB2IPZmUGtBmTmvZl5OpUH13+yp+tTw/4R+Eu/eA7ZE1SG66YC/wx8t4frU6t6ARcAbwauBD4WEaf2bJVq3luBhzNzfU9XpEZdCTwFnAicC3wuIgb1bJVq1meAIRHxFJWRnCeBg7onbCme9XkMq+pxWbtk5i8jYmJEDM9Mn8m2p2rashG4p9Krz3DgmohozUyDxiu6bMfM3Njh9f0R8QU/k69SzedxGdCcmVuALRHxS+AcYOGRqWLNOJDvyetx2HNfqmnHm4DPFKfbNEXEEirnV/l0oD1V+z15E0AxlLwEWHwwO7NHrWd1+XiuiJhU/JKJiPOBvkDzEa9p+XXZlpk5ITNPzsyTgW8B7zekvUo1n8lRHT6TF1L5HvEzuadqHr33PeDSiOhVDNldBMw/wvWsBdW0JRExGHgdrzyOUHuqph1fBK4AKM6nOo2DDBdHuWq+J4cUywD+CPhlx//kHgh71HpQVvd4rndSeY7qTmAb8K4OFxeoUGVbqgtVtuN1wJ9GRCuVz+T1fib3VE07Zub8iPgRMBtoB/41M+f2XK3L6QD+bV8L/LjooVQnVbbjJ4EvR8QcIKicKmJPeSdVtuUZwFciIoF5wB8e7P58MoEkSVJJOfQpSZJUUgY1SZKkkjKoSZIklZRBTZIkqaQMapIkSSVlUJNUkyKiLSKe6jCdHBGNEXHHPsrvc9lhrOO5EXHNkdynpKOLt+eQVJMiYnNmDqiybK/MPCwPl97ftiPiRqAxM287HPuWdPSzR03SUSMiLo+IHxSvPx4R/x4RDwP/vpdldxcPRF8cEX/WYRsfi4gFEfHriPh6RPzPveznyxFxV0T8Fvj7iLgwIh6JiCcj4jcRcVpxV/JPAO8qevzeFRH9i/0+VpSdfmRaRlKt8skEkmrVccUDjwGWZOa1eykzBbg0M7dFxOWdlp0OvB4YCCyIiDupPIj6nVSeudmbygPoZ+1j/2OBizOzrXhw9WXFHct/F/hUZr4zIv6aDj1qEfEp4KeZeXNEDAEei4gHvZu+pH0xqEmqVdsy89wuyszIzG37WHZfZrYALRGxBhgJXAJ8LzO3A9sj4vv72fY3M7OteD2YyuNiJgNJJeTtzZuAt3XopesHjMdnfEraB4OapKPZ/nqqWjq8buPAvw87bvuTwM8y89qIOBn4+T7WCeCdmbngAPcl6RjlOWqS9IqHgbdGRL+IGAC8pcr1BgPLi9c3dpi/icrQ6i4PAB+MiACIiPMOrbqSjnYGNUkqZObjwAxgNvBDYA6woYpV/x74dEQ8yZ49cz8Dpuy6mIBKz1tvYHZEzCveS9I+eXsOSeogIgZk5uaIOB74JXBLZj7R0/WSdGzyHDVJ2tOXImIKlRP9v2JIk9ST7FGTJEkqKc9RkyRJKimDmiRJUkkZ1CRJkkrKoCZJklRSBjVJkqSSMqhJkiSV1P8PmtcFWTm37EsAAAAASUVORK5CYII=",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6594fe6d90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(rates/tsim, diffs)\n",
    "plt.ylim([0, 0.003])\n",
    "plt.title(\"$\\mathbb{E}(S|h= 1) - \\mathbb{E}(S|h=0)$. c = 0.01\")\n",
    "plt.xlabel(\"Firing rate\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7f658d358450>"
      ]
     },
     "execution_count": 149,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARYAAAEYCAYAAACZR9k/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAGYxJREFUeJzt3Xu4VXW97/H3h4su7yhg4WKjeJQKFLyslE22RfSEuFXyUbdaiaT7oWPmrux46mSZmeWx2qd0742GlxD3fjJNUzIvsVW2p8jLwlQEQkxMQA1C8goq8D1/jLFoMtdtLviNOdZafl7PMx/WHGPMMb7rx1qfNS5zjq8iAjOzlPqUXYCZ9T4OFjNLzsFiZsk5WMwsOQeLmSXnYDGz5BwsZpacg8XMknOwdCOSLpf0hRqXvUTSJW1MXy7pkOTFtV/Ho5JG1Wt7+Ta3aZzqPUb5Nus+TmVysHQTkgYDU4AfVU3/pKT5kv4iaa2kuZIa21nH7sBewOLEtX1OUrOktyXNrJr9feDSlNvrpJZW49RNxmgPST+X9KakP0r6RNUidR2nsvUru4D3Gkn/AewNPF8x+dX8+d0Rsa5i2anA14DTgMeBgcDHgbXtrP5A4PnKdSTyInAZMBHYoWrebOAaSe+PiJdTbKy9MYqI84CpVIxTNxqjfwPeAd4HHAT8UtKTEbEwn598nLozB0sdSRoI/Bb4dURcXTXvAeCGqpecA/woIubnz/8MXJcv39YmRgPPSroKOAN4FzgrIuZsS90RcXu+zSZgaNW89ZLmk4XOjduynXwb7Y5RbhJbjlO7Y9SO5GMkaSfgZOCAiHgD+LWkO4Ezga9A+nHq7nwoVF9HA3OAt9qYdyCwpGraOuBsSf8gaVAN6z8QaAJ+SfaX80fAlysXkHRXfsjQ1uOurn5DucXAmK18bbWOxghaj1PyMYIuj9MIYENEPFMx7Umg+pxKynHq1hws9bUv2aHEpyT9e/44NZ83AHi9avkpwH1kx+d/kvQLSXt2sP7RwLcj4r6I2AQsql4gIo6PiAHtPI7fyu/r9bz+FDoaI2g9TsnHCLo8TjsDr1VNew3YpWpaynHq1hws9bUBOByYGBGfyh+35vPWUvWDGBEvR8QXImIYcBjZL0Wrv64Ayo6NDgB+UTH5ANr5xUlsF+Av7dQ1V1K08/h1Gy/paIygapy6yRi9AexaNW03Wv+haHecehsHS51I2hd4Btgu/0tZ7SmyXeo25ecQFgA7tbPI8PzfZyumHQw8UVXHPZLeaOdxT63fT5UPke36t1X3+IhQO48jqmrrbIygg3FKNUZ5LV0Zp2eAfpL2r5g2BlhYtVy749TbOFjq56PAo2TnBNpyN3BkyxNJX5H0EUnb54+pwHjgx+28fjSwILa8c9fBVP0gR8SkiNi5ncektlYsqZ+kBqAv0FdSg6R++bwG4FCy8yLbqrMxgopxKmqMoGvjFBFvArcDl0raSdIRwInATS3LJB6nbs9XhernI8DpQB9J4yum3xARDwCzgCck7ZBfCt2V7BdkL+BNskupR0fEI+2s/0AqfkHyE5nvB55OUPvXgG9UPP8U8E3gEuAEYG5EvJhgO52NEVSME91rjD5LdrVqFbAGOLfiUjOkHaduT741Zfch6TvAqoj4YQ3LXgIQEZcUXFZndTwCnBMRKX45a92mx6mb8x5LNxIRXy27hq6KiMNL2KbHqZtzsPRcc8suoIeYW3YB70U+FDKz5HxVyMyS64mHQt7FMitHmx9Qa4v3WMwsOQeLmSXnYDGz5BwsZpZcTzx5a1azd999lxUrVrB+/fqyS+kxGhoaGDp0KP3799/qdfTE97H0uIKtPMuWLWOXXXZh4MCB7d11zypEBGvWrOH1119n+PDh1bN9VcgMYP369Q6VLpDEwIEDt3kPz8FivZ5DpWtSjJeDxcyS88lbe09ZsOLVpOs7cOhuSddXi/nz5zN16lTWrVvHcccdx5VXXtnmXsbll1/O9ddfT9++fbnqqquYOHEiABdddBGzZs1i7dq1vPHGG4XUWNgei6QbJK2S1OH9JyR9WNIGSacUVYtZb3Luuedy7bXXsnTpUpYuXcq9997baplFixZx8803s3DhQu69914++9nPsnHjRgBOOOEEHn300UJrLPJQaCZwbEcLSOoLXAH8qsA6zEo1a9YsRo8ezZgxYzjzzDO3aV0vvfQSr732GmPHjkUSU6ZM4Y477mi13J133snpp5/O9ttvz/Dhw9lvv/02h8nYsWMZMmTINtXRmcIOhSLiIUn7dLLY+cBtwIeLqsOsTAsXLuSyyy5j3rx5DBo0iFdeeaXVMg8++CBf/OIXW03fcccdmTdv3hbTVq5cydChf+0ZN3ToUFauXNnqtStXrmTs2LGdLleU0s6xKOutexJwFA4W66UeeOABTj31VAYNynqp7bHHHq2WOeqoo3jiiVaNAnq0Mk/e/hD4ckRs6uzylqRpwDSAkSNHsnBhdVcFs56rK3ssjY2NrFixYvPzFStW0NjY2Oq1jY2NLF++vNPlilLm5eYm4GZJzwOnANMlfbytBSNiRkQ0RUTTDjtU9yQ3674mTJjArbfeypo1awDaPBRq2WOpflSHCsCQIUPYddddefjhh4kIZs2axeTJk1std+KJJ3LzzTfz9ttvs2zZMpYuXcphhx2W/htsT0QU9gD2AZ6uYbmZwCm1rPPQQw8Ns1otWrSo7BJi5syZMWrUqBg9enScddZZ27y+xx57LEaNGhX77rtvnHfeebFp06aIiLjzzjvj61//+ublLrvssth3331jxIgRcffdd2+efuGFF0ZjY2NIisbGxvjGN77RahvtjFvNv/uFfVZI0k/ImkcNAv5E1pemfx5m11QtOxO4KyJ+1tl6m5qaorm5OXW51kstXryYD33oQ2WX0eO0M241vyW3yKtCZ3Rh2alF1WFm9ee39JtZcg4WM0vOwWJmyTlYzCw5B4uZJefbJth7y4u/S7u+vQ5Ou74a1HLbhHfeeYfPfOYzNDc306dPH6688krGjx8PwE9/+lO+/e1vs3HjRo4//niuuOKK5DV6j8Wsh6nltgnXXnstAAsWLGDOnDl86UtfYtOmTaxZs4YLL7yQ+++/n4ULF/Lyyy9z//33J6/RwWJWsDJum7Bo0SImTJgAwJ577smAAQNobm7mueeeY//992fw4MEAHHPMMdx2223bVFNbfChkVqCybpswZswYZs+ezRlnnMHy5cuZP38+y5cvZ8KECSxZsoTnn3+eoUOHcscdd/DOO+8k+E635GAxK1BZt004++yzWbx4MU1NTey9996MGzeOvn37svvuu3P11Vdz2mmn0adPH8aNG8cf/vCHpNsGB4tZ6Yq4bUK/fv34wQ9+sPn5uHHjGDFiBJDdmvKEE04AYMaMGfTt2zfJ97HF9pOv0cw2mzBhAieddBIXXHABAwcO5JVXXmm119KVPZbK2yYcfvjhzJo1i/PPP7/Vcm+99RYRwU477cScOXPo168fI0eOBGDVqlXsueeerF27lunTp3PLLbds+zdaxcFi7y11vjw8atQoLrroIo488kj69u3LwQcfzMyZM7dpndOnT998uXnSpElMmjQJgNmzZ9Pc3Myll17KqlWrmDhxIn369KGxsZGbbrpp8+s///nP8+STTwJw8cUXb96TSanHtVj1bROsK3zbhK2zrbdN8OVmM0vOwWJmyTlYrNfraYf7ZUsxXg4W69UaGhpYs2aNw6VGEcGaNWtoaGjYpvX4qpD1akOHDmXFihWsXr267FJ6jIaGhi3e3bs1HCzWq/Xv35/hw4eXXcZ7TmlN4SV9UtJTkhZImidpTFG1mFl9ldkUfhlwZEQcCHwLmFFgLWZWR6U1hY+Iyg9BPAxs20GdmXUb3eWq0DnAPWUXYWZplH7yVtJRZMFyRAfLbG4KP2zYsDpVZmZbq9Q9FkmjgeuAyRGxpr3loqIpfMudr8ys+yotWCQNA24HzoyIZ8qqw8zSK+xQqLIpvKQVtG4KfzEwEJie32F8Q0Q0FVWPmdVPaU3hI+IfgX8savtmVp7uclXIzHoRB4uZJedgMbPkHCxmlpyDxcySc7CYWXIOFjNLzsFiZsk5WMwsOQeLmSXnYDGz5BwsZpacg8XMknOwmFlyDhYzS87BYmbJOVjMLDkHi5kl52Axs+QcLGaWXJlN4SXpKknP5s3hDymqFjOrrzKbwk8C9s8f04CrC6zFzOqosGCJiIeAVzpYZDIwKzIPAwMkDSmqHjOrnzJ7NzcCyyuer8invdTRi5YsWcL48eMLLMvM2jJ37tyaly29KXwtKpvCb7/99iVXY2adKTNYVgJ/U/F8aD6tlYiYAcwAaGpqiq4kp5nVX5mXm2cDU/KrQ2OBVyOiw8MgM+sZymwKfzdwHPAs8Bbw6aJqMbP6KrMpfADnFbV9MyuP33lrZsk5WMwsOQeLmSXnYDGz5BwsZpacg8XMknOwmFlyDhYzS87BYmbJOVjMLDkHi5kl52Axs+QcLGaWnIPFzJJzsJhZcg4WM0vOwWJmyTlYzCw5B4uZJedgMbPkCg0WScdKWpI3fv9KG/N3k/QLSU9KWijJd+o36wUKCxZJfYF/I2v+PhI4Q9LIqsXOAxZFxBiyViH/LGm7omoys/ooco/lMODZiHguIt4BbiZrBF8pgF0kCdiZrIn8hgJrMrM6qLmvkKS/B0YBDS3TIuLSDl7SVtP3w6uW+VeyjogvArsAp0XEplprMrPuqaY9FknXAKcB5wMCTgX2TrD9icATwF7AQcC/Stq1je1Pk9QsqXn16tUJNmtmRar1UGhcREwB1kbEN4G/BUZ08ppamr5/Grg9Ms8Cy4APVq8oImZERFNENA0ePLjGks2sLLUGy7r837ck7QW8Cwzp5DWPAftLGp6fkD2d7LCn0gvA0QCS3gd8AHiuxprMrJuq9RzLXZIGAN8DHic76XpdRy+IiA2SPgfcB/QFboiIhZL+Rz7/GuBbwExJC8gOsb4cEX/eum/FzLoLZb3ZO1lI2j4i3m75muwE7vqWafXU1NQUzc3N9d6smWV//GtS66HQb1u+iIi3I+LVymlmZpU6PBSS9H6yy8Y7SDqYvybWrsCOBddmZj1UZ+dYJgJTya7o/N+K6a8BXy2oJjPr4ToMloi4EbhR0skRcVudajKzHq7Wcyy/kXS9pHsAJI2UdE6BdZlZD1ZrsPyY7LLxXvnzZ4AvFFKRmfV4tQbLoIi4BdgE2XtUgI2FVWVmPVqtwfKmpIFkb4xD0ljg1cKqMrMerdZ33l5A9nb8fSX9BhgMnFJYVWbWo9UaLIuAnwNvAa8Dd5CdZzEza6XWQ6FZZJ86/g7wL2SfbL6pqKLMrGerdY/lgIiovK3kg5IWFVGQmfV8te6xPJ6fsAVA0uGAPwloZm3q7LNCC8iuBPUH5kl6IX++N/D74sszs56os0Oh4+tShZn1Kp19VuiP9SrEzHoPd0I0s+QcLGaWnIPFzJJzsJhZcqU2hc+XGS/pibwp/H8VWY+Z1UfNLVa7qqIp/H8na6/6mKTZEbGoYpkBwHTg2Ih4QdKeRdVjZvVTdlP4T5B1QnwBICJWFViPmdVJkcHSVlP4xqplRgC7S5orab6kKQXWY2Z1UtihUBe2fyhZm9UdgN9Kejgitrglg6RpwDSAYcOG1b1IM+uaIvdYamkKvwK4LyLezFurPgSMqV6Rm8Kb9SxFBkstTeHvBI6Q1E/SjsDhwOICazKzOijsUKiWpvARsVjSvcBTZDfqvi4ini6qJjOrj5qawncnbgpvVprkTeHNzGrmYDGz5BwsZpacg8XMknOwmFlyDhYzS87BYmbJOVjMLDkHi5kl52Axs+QcLGaWnIPFzJJzsJhZcg4WM0vOwWJmyTlYzCw5B4uZJedgMbPkHCxmlpyDxcySK70pfL7chyVtkHRKkfWYWX0UFiwVTeEnASOBMySNbGe5K4BfFVWLmdVX2U3hAc4HbgPcEN6slyi1KbykRuAk4OoC6zCzOiv75O0PgS9HxKaOFpI0TVKzpObVq1fXqTQz21qFtViltqbwTcDNkgAGAcdJ2hARd1QuFBEzgBmQdUIsrGIzS6LIYNncFJ4sUE4HPlG5QEQMb/la0kzgrupQMbOep9Sm8EVt28zK5abwZlYrN4U3s/I4WMwsOQeLmSXnYDGz5BwsZpacg8XMknOwmFlyDhYzS87BYmbJOVjMLDkHi5kl52Axs+QcLGaWnIPFzJJzsJhZcg4WM0vOwWJmyTlYzCw5B4uZJedgMbPkSm0KL+mTkp6StEDSPEljiqzHzOqj7Kbwy4AjI+JA4FvkTcnMrGcrtSl8RMyLiLX504fJuiWaWQ9XalP4KucA9xRYj5nVSZEtVmsm6SiyYDminfnTgGkAw4YNq2NlZrY1itxjqaUpPJJGA9cBkyNiTVsriogZEdEUEU2DBw8upFgzS6fIYNncFF7SdmRN4WdXLiBpGHA7cGZEPFNgLWZWR2U3hb8YGAhMlwSwISKaiqrJzOrDTeHNrFZuCm9m5XGwmFlyDhYzS87BYmbJOVjMLDkHi5kl52Axs+QcLGaWnIPFzJJzsJhZcg4WM0vOwWJmyTlYzCw5B4uZJedgMbPkHCxmlpyDxcySc7CYWXIOFjNLzsFiZsmV3RRekq7K5z8l6ZAi6zGz+ii7KfwkYP/8MQ24uqh6zKx+Sm0Knz+fFZmHgQGShhRYk5nVQdlN4bvaON7MeoBu0RS+M5VN4YGFwAFdeO3TwPoi6tpKg4A/l11EBdfTse5WT0NE1PzzX5Yig6WWpvA1NY6PiBnAjK2sY313atsqqdn1tM/1dExSj2gDWmpT+Pz5lPzq0Fjg1Yh4qcCazKwOym4KfzdwHPAs8Bbw6aLqMbP6KfQcS0TcTRYeldOuqfg6gPOKrIGtP4QqiuvpmOvpWHerp03KfrfNzNLxW/rNLLleHSySBkj6maTfS1os6W9LrueLkhZKelrSTyQ11Hn7N0halV+Cb5m2h6Q5kpbm/+5ecj3fy/+/npL0c0kDyqynYt6XJIWkQWXXI+n8fIwWSvpuverpil4dLMCVwL0R8UFgDLC4rEIkNQL/BDTl70PoS3alrJ5mAsdWTfsKcH9E7A/cnz8vs545wAERMRp4BvjfJdeDpL8BPga8UMda2qxH0lFk71gfExGjgO/Xuaaa9NpgkbQb8HfA9QAR8U5E/KXcqugH7CCpH7Aj8GI9Nx4RDwGvVE2eDNyYf30j8PEy64mIX0XEhvzpw2TvbSqtntwPgP8F1PWEZDv1nAv8n4h4O19mVT1rqlWvDRZgOLAa+LGk30m6TtJOZRUTESvJ/rq8ALxE9p6dX5VVT4X3Vbx36GXgfWUWU+Vs4J4yC5A0GVgZEU+WWUeFEcBHJT0i6b8kfbjsgtrSm4OlH3AIcHVEHAy8SX1387eQn7uYTBZ4ewE7SfpUWfW0Jb/83y0uE0q6CNgA/EeJNewIfBW4uKwa2tAP2AMYC1wI3CJJ5ZbUWm8OlhXAioh4JH/+M7KgKcsxwLKIWB0R7wK3A+NKrKfFn1o+UZ7/W/qutaSpwPHAJ6Pc90P8N7I/BE9Kep7ssOxxSe8vsaYVwO35HQEeBTaRfZ6pW+m1wRIRLwPLJX0gn3Q0sKjEkl4AxkraMf8LczQlnkyuMBs4K//6LODOEmtB0rFk5zNOjIi3yqwlIhZExJ4RsU9E7EP2S31I/rNVljuAowAkjQC2o3t9SDITEb32ARwENANPkf2H7F5yPd8Efg88DdwEbF/n7f+E7PzOu2S/JOcAA8muBi0F/hPYo+R6niW7lcYT+eOaMuupmv88MKjk8dkO+Pf8Z+hxYEI9f4Zqffidt2aWXK89FDKz8jhYzCw5B4uZJedgMbPkHCxmlpyDpReTtFHSExWPfSQ1SbqqneXbnVdgjQdJOm4rXjdE0l351+Nbvq6YP1PSKR28/vuSJnS9YqtFj7hLv221dRFxUNW058ne27MFSf0iormtedsqX/eGdmYfBDRRdafBGlwAXLsNZf1L/voHtmEd1g7vsbzHVP51l3SJpJsk/Qa4qY15N0iaK+k5Sf9UsY6vK2ud++v8vjL/s43tzJR0jaRHgO9KOkzSb/MPhM6T9IH8JuuXAqfle1SnSdop3+6j+bLVTe5anAzcW8P321Sxx7ZAUgBExB+BgSW/Pb/X8h5L77aDpCfyr5dFxEltLDMSOCIi1kkaXzXvg2RvH98FWCLparI9jJPJ7m/Tn+zdn/Pb2f5QYFxEbJS0K/DRyG6yfgzwnYg4WdLFZPeo+RyApO8AD0TE2flNnh6V9J8R8WbLSiUNB9ZGfuuA3EcrvleAYcBd+V7YQfnrvseWYfQ48BHgtnbqt63kYOnd2joUqjY7Ita1M++X+S/v25JWkd1S4SPAnRGxHlgv6RcdrPvWiNiYf70bcKOk/ck+Qd2/ndd8DDixYi+ogSwkKj9XNYTslhiV/l9EHN/yRNLMypmSTiP7EOrHKiavIvukuSXmYLE3O5hXuUewka7/vFSu+1vAgxFxkqR9gLntvEbAyRGxpIP1riMLnJpIOgC4BPi7iqAjX0d7oWrbwOdYrKt+A5wgqUHSzmS3N6jFbvy1y+XUiumvkx1qtbgPOL/lHiOSDm5jXc8A+9Sy0fxw6ifAlIio3ssZQfZhPkvMwWJdEhGPkd1q4Smyu7stAF6t4aXfBS6X9Du23PN5EBjZcvKWbM+mP/CUpIX58+oa3gT+IGm/GrY7GdgbuLblJC6ApP7AfhRwFczcV8i2gqSdI+KN/A5rDwHTIuLxOtdwEnBoRHxtG15/SER8PW1lBj7HYltnhqSRZOcobqx3qABExM8lDdyGVfQD/jlVPbYl77GYWXI+x2JmyTlYzCw5B4uZJedgMbPkHCxmlpyDxcyS+/8ya6zzUE1ZWgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f658d408690>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax1 = plt.subplots(1,1,figsize = (4,4))\n",
    "ds1 = exp_filter[0]\n",
    "ds = 1./DeltaT*(1 - np.exp(-DeltaT*dt/(tau_s)))/2\n",
    "#ds = 1/T*(1 - np.exp(-T/tauS))\n",
    "sns.tsplot(data = diffsc_0_99/ds, time = ratesc_0_99/tsim, ax = ax1)\n",
    "sns.tsplot(data = diffsc_0_01/ds, time = ratesc_0_01/tsim, ax = ax1, color = 'C1')\n",
    "ax1.set_xlabel(\"Firing rate (Hz)\")\n",
    "ax1.set_ylim([0, 1.5]);\n",
    "ax1.set_ylabel('beta')\n",
    "ax1.set_title(\"$\\mathbb{E}(S|h= 1) - \\mathbb{E}(S|h=0)$\")\n",
    "sns.despine(ax = ax1, trim = True)\n",
    "ax1.legend({'c = 0.01', 'c = 0.99'})\n",
    "plt.savefig('./fig_s1_c.pdf')\n",
    "ax1.axhline(y = 1, color = (0,0,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is the theoretical term $\\Delta_s$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Compute causal effect for range of DeltaT and tau_s values\n",
    "N = 51\n",
    "deltas = np.array([10, 20, 50, 100, 200])\n",
    "M = deltas.shape[0]\n",
    "tauss = np.logspace(-2, -0.5, N)\n",
    "\n",
    "#Run a long simulation\n",
    "lif.W = w\n",
    "(v, h, _, _, u) = lif.simulate()\n",
    "\n",
    "#Compute s\n",
    "len_h = h.shape[1]\n",
    "\n",
    "ce = np.zeros((M,N))\n",
    "ce_theory = np.zeros((M,N))\n",
    "\n",
    "#Convolve according to taus, bin according to deltas, compute causal effect\n",
    "for idx, DeltaT in enumerate(deltas):\n",
    "    for idx2, tau_s in enumerate(tauss):\n",
    "        t_filter = np.linspace(0, 5, 5000)\n",
    "        exp_filter = np.exp(-t_filter/tau_s)\n",
    "        exp_filter = exp_filter/np.sum(exp_filter)\n",
    "        ds = exp_filter[0]\n",
    "        s1 = np.convolve(h[0,:], exp_filter)[0:len_h]\n",
    "        s_r = s1[:].reshape((-1, DeltaT))\n",
    "        S = np.squeeze(s_r[:,-1])\n",
    "        H_r = h[0,:].reshape((-1, DeltaT))\n",
    "        H = np.sum(H_r, 1)\n",
    "        ce[idx,idx2] = np.mean(S[H == 1]) - np.mean(S[H == 0])\n",
    "        #Compare this to the theoretcal effect\n",
    "        ce_theory[idx,idx2] = 1./DeltaT*(1 - np.exp(-DeltaT*dt/(tau_s)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVgAAAFECAYAAABieCagAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4lFX6xvHvk9AFJdKlKYIFLAhRYBVWpBcFZFFALPxkERdZu4uKqKisritrWRQREVmkSy/SpK7SAoogskY0CCLNgPS08/tjJjiMCcxApiX357pyJfO+78w8CXBzct5TzDmHiIjkvbhIFyAikl8pYEVEQkQBKyISIgpYEZEQUcCKiISIAlZEJEQUsCIiIaKAFREJEQWsiEiIFIp0ASGkKWoiEgoW6IVqwYqIhIgCVkQkRBSwIiIhooAVEQkRBayISIgoYEVEQkQBKyISIgpYEZEQUcCKiISIAlZECqyklFSGLk4mKSU1JK+fn6fKiojkKikllTtGrCQtI4siheL4qFdD6ldPyNP3UAtWRAqklVv3kZaRRZaD9IwsVm7dl+fvoYAVkQKpYY0yFCkUR7xB4UJxNKxRJs/fw5zLt4tO5dtvTETyRlJKKiu37qNhjTLBdA8EvJqWAlZEJDharlBEJNIUsCIiIaKAFREJEQWsiEiIKGBFJF8J9eysYGgml4jkG+GYnRWMsLZgzay1mW0xs2Qz65/DeTOzN73nN5hZPZ9zpc1sspl9Y2abzaxROGsXkegXjtlZwQhbwJpZPDAUaAPUBrqZWW2/y9oAtbwfvYF3fM69AXzinLsMuBrYHPKiRSSmnG52VlpaWljrCWcL9jog2Tm31TmXBowHOvhd0wEY7TxWAqXNrJKZnQc0Ad4HcM6lOef2h7F2EYkB9asn8FGvhjzS8tLfdQ/88MMPXHnllUycODFs9YSzD7Yy8KPP4+1AgwCuqQxkAHuAD8zsaiAJeNA5d9j3yWbWG0/Ll9q1a7Np06Y8/QZEJPrVr57wu37Xb775hhYtWnDo0CGqVasWtlpiZRRBIaAe8I5z7hrgMPC7Plzn3HDnXKJzLrF48eLhrlFEotC6deto3Lgx6enpLF26lIYNG4btvcMZsDuAqj6Pq3iPBXLNdmC7c26V9/hkPIErIpKrFStW0LRpU0qUKMHy5cu56qqrwvr+4QzYNUAtM7vIzIoAXYEZftfMAO7yjiZoCBxwzu10zv0M/Ghml3qvawZ8HbbKRSTmzJs3j5YtW1KpUiVWrFhBrVq1wl5D2PpgnXMZZvYAMA+IB0Y65zaZWR/v+WHAHKAtkAwcAXr6vEQ/4CNvOG/1OycicsLkyZPp3r07derUYd68eZQvXz4ideTb5QoTExPd2rVrI12GiOSBYNZt/eCDD+jVqxeNGjVi1qxZlC5dOq/LCXi5Qs3kEpGoFszsrDfeeIOHHnqIli1bMmXKFM4555wwV3uyWBlFICIFVCCzs5xzPP/88zz00EN07tyZGTNmRDxcQS1YEYly2bOz0jOycpyd5Zzj0Ucf5V//+hf33HMP7733HoUKRUe0RUcVIiK5yJ6dlVMfbGZmJr1792bkyJE8+OCDDBkyhLi46PnFXAErIlEvp9lZaWlp9OjRg0mTJvHss8/y7LPPYhbw/aewUMCKSMw5cuQInTt35pNPPmHIkCE8/PDDkS4pRwpYEYkpBw4coH379nz22WeMGDGCe++9N9Il5UoBKyIxY8+ePbRu3ZqvvvqK8ePH06VLl0iXdEoKWBGJCTt27KB58+akpKQwffp02rRpE+mSTksBKyIREczsrOTkZFq0aMG+ffuYN28ejRs3DlOVZ0cBKyJhF8zsrI0bN9KiRQvS09NZvHgx9evXD3O1Zy56BoyJSIER6N5Zq1evpkmTJsTFxbFs2bKYCldQwIpIBJxu7yyAxYsX06xZMxISElixYgW1a/tv4Rf91EUgImF3qtlZADNnzqRLly7UrFmTBQsWUKlSpQhVenYUsCISETnNzgIYO3Ysd911F/Xq1WPu3LmUKfP71m2sUBeBiESNYcOG0aNHDxo3bsyiRYtiOlxBASsiUeKVV17h/vvvp127dsyZM4dSpUpFuqSzpoAVkYhyzvHUU0/Rv39/unXrxpQpU8gvu0KrD1ZEIiYrK4t+/frx9ttvc9999zF06FDi4+MjXVaeUQtWRPJMUkoqQxcnk5SSetpr09PTufvuu3n77bd54okneOedd/JVuIJasCKSR4KZnXXs2DG6du3K9OnTGTx4ME8++WSYqw0PtWBFJE8EOjvr0KFDtGvXjunTpzN06NB8G66gFqyI5JHT7Z0F8Msvv9C2bVvWrl3Lf/7zH3r06BGBSsNHASsieeJ0s7N+/vlnWrZsyZYtW5g8eTIdO3aMUKXho4AVkTyT2+yslJQUmjdvzs6dO5kzZw7NmjWLQHXhp4AVkZDasmULzZs359ChQyxYsIBGjRpFuqSwUcCKSMisX7+eVq1aERcXx9KlS7nqqqsiXVJYaRSBiITEf//7X5o2bUrx4sVZvnx5gQtXUMCKSACCmUAAMH/+fFq2bEmFChVYvnw5tWrVCnGF0UldBCJySsFMIAD4+OOP6datG3Xq1GHevHmUL18+jNVGl7C2YM2stZltMbNkM+ufw3kzsze95zeYWT2fcz+Y2Vdm9oWZrQ1n3SIFWaATCABGjRrFbbfdxrXXXsvixYsLdLhCGAPWzOKBoUAboDbQzcz894BoA9TyfvQG3vE739Q5V9c5lxjqekXEI5DtXQDefPNNevbsSbNmzZg/fz6lS5cOc6XRJ5xdBNcByc65rQBmNh7oAHztc00HYLRzzgErzay0mVVyzu0MY50i4uN0Ewicc7z44osMHDiQW2+9lbFjx1K0aNEIVRtdwhmwlYEffR5vBxoEcE1lYCfggIVmlgm865wb7v8GZtYbT8uXatWq5V3lIgVcbhMInHM89thjDBkyhLvvvpsRI0ZQqJBu7WSLpVEENzjn6uLpRuhrZk38L3DODXfOJTrnEsuVKxf+CkUKkMzMTP785z8zZMgQ+vXrx8iRIxWufsIZsDuAqj6Pq3iPBXSNcy77825gKp4uBxGJgLS0NLp3787777/PM888wxtvvEFcXCy118IjnD+RNUAtM7vIzIoAXYEZftfMAO7yjiZoCBxwzu00s3PMrBSAmZ0DtAQ2hrF2EfE6cuQIHTt2ZOLEifzzn/9k0KBBmFmky4pKYWvPO+cyzOwBYB4QD4x0zm0ysz7e88OAOUBbIBk4AvT0Pr0CMNX7h1gIGOuc+yRctYvkV0kpqbnevMrJr7/+Svv27VmxYgXvvfcevXr1CkOVscs8N+zzn8TERLd2rYbLiuQm2AkEe/fupXXr1nz55ZeMGTOG22+/PYzVRpWAm+vqNBEpoIKZQLBjxw6aNGnCpk2bmD59ekEO16Dolp9IARXIDgQA3333Hc2bN2ffvn3MmzePJk1+N4BHcqGAFSmgTjeBAGDjxo20bNmStLQ0Pv30UxITNYkyGApYkQIstwkEAKtXr6ZNmzYUK1aMZcuWUbu2/8x2OR31wYrI7yxZsoRmzZpRunRpVqxYoXA9QwpYETnJrFmzaN26NdWqVWP58uVcdNFFkS4pZilgReSEcePG0alTJ6688kqWLVvGBRdcEOmSYpoCVkQAGD58OHfccQfXX389ixYtokyZnEcVSOAUsCL5ULBbvLz66qvcd999tG3blrlz53LuueeGuMKCQaMIRPKZYGZoOecYMGAAgwcPpmvXrowePZrChQuHueL8Sy1YkXwm0BlaWVlZ9OvXj8GDB9O7d2/GjBmjcM1jCliRfCaQLV4yMjK45557GDp0KI8//jjDhg0jPj4+AtXmb+oiEMlnTjdD69ixY3Tr1o1p06bx0ksv8eSTT2q5wRBRwIrkQ7nN0Dp06BAdO3Zk0aJFvPXWWzzwwAMRqK7gUMCKFBCpqam0bduWNWvW8OGHH3LXXXdFuqR8TwErUgDs2rWLli1b8s033zBp0iQ6deoU6ZIKBAWsSIwJdheClJQUWrRowY4dO5g9ezbNmzcPQ5UCCliRmBLsLgRbtmyhRYsWHDx4kIULF9KoUaMwVisapiUSQ3Ib45rTzK0vvviCxo0bc/z4cZYsWaJwjQC1YEViSE67EOTUqj2+YzNt27bl3HPPZeHChVxyySWRLr1AUsCKxJCcxrgOXZx8Uqv2g4nT+eC5vlSpUoUFCxZQrVq1SJddYClgRWKM/xhX31btseTPGT7jVWpffjnz5s2jQoUKEaxU1AcrEmP8+1uzW7V/cJv4ecrfSaxfn8WLFytco4BasCIxJLdRBJ/NGMOYV/5GixYtmDp1Kuecc06kSxUUsCIxxX8Uweff7WXuf4byzDPP0KlTJ8aNG0fRokUjXaZ4KWBFokQgEwh8+1sLxRurJ7zJf4b/m7vuuov333+fQoX0Tzqa6E9DJAqcagKBf/B+1Kshn327m8UjB/OfcaN54IEHeOONN4iL0y2VaKOAFYkCOU0gqF89IcfgvbLSObz6t2eZOmECAwYMYNCgQVpuMErpvzyRKJDbItn+wbts8w46derEhAkTePXVV3nhhRcUrlFMLViRKJDbItm+fa5xmccY8+x9rF/zOcOHD+fPf/5zhKuW0wlrC9bMWpvZFjNLNrP+OZw3M3vTe36DmdXzOx9vZuvNbFb4qhYJj/rVE2hYowwrt+773RjX3g3KUWTeS2xYt5qxY8cqXGNE2FqwZhYPDAVaANuBNWY2wzn3tc9lbYBa3o8GwDvez9keBDYD2lNY8p3cbnRVLHSE0U/3ZOvWrUybNo127dpFulQJUDhbsNcByc65rc65NGA80MHvmg7AaOexEihtZpUAzKwK0A4YEcaaRcImpxtdW7dupXHjxmzbto1PPvlE4RpjwhmwlYEffR5v9x4L9JrXgSeArNzewMx6m9laM1u7Z8+es69YJA/ltKSgL/8bXeUydnPDDTdw4MABPv30U/74xz+GuWI5WzFxk8vM2gO7nXNJZnZjbtc554YDwwESExNdmMoTOa1AFsr2vdFV6mAK93e9maJFi7Js2TLq1KkTocrlbIQzYHcAVX0eV/EeC+SazsAtZtYWKAaca2ZjnHM9QlivSJ7JbaFs/1ED9asncOiHDdzcoxNly5Zl4cKF1KhRI5Kly1kIZxfBGqCWmV1kZkWArsAMv2tmAHd5RxM0BA4453Y65550zlVxzl3ofd6nCleJJf6//ieUKMIdI1by2vwt3DFi5Ylug9mzZ9O6dWuqVq3K8uXLFa4x7rQtWDN75FTnnXNDAnkj51yGmT0AzAPigZHOuU1m1sd7fhgwB2gLJANHgJ6BvLZItPMf55pTi/bbz+dx5513UrduXebOnUvZsmUjXbacpUC6CEp5P18KXMtvrc6bgdXBvJlzbg6eEPU9Nsznawf0Pc1rLAGWBPO+ItHAf6Fs361fflo5k35PP0rjxo2ZOXMm556rkYj5gXkyLYALzZYB7ZxzB72PSwGznXNNQljfGUtMTHRr166NdBkiQM4rZWUf+3bhWN4Y/Cxt27Zl8uTJFC9ePMLVymkEPDc5mJtcFYA0n8dp3mMicgq5jSCoV60000YM4Y3BL3L77bczevRoihQpEulyJQ8FE7CjgdVmNtX7uCMwKs8rEolBp1rLNaf+1muqnsdDDz3EW2+9Ra9evRg2bBjx8fERql5CJeCAdc69ZGZzgcbeQz2dc+tDU5ZI7DjdGFf/rbYTq51Hz549GT16NI8++iivvvqqVsTKp4IaB+ucWwesC1EtIjEpt7Vcs/mOIKhXuSSvPH4fU6dO5YUXXuDpp59WuOZjAQesef4W3AHUcM4NMrNqQEXnXFAjCUTyG/8WakKJIgxdnPy7CQSXlS1Cp06dWLBgAW+++Sb9+vWLcOUSasGMIngHzzoANznnLjezBGC+c+7aUBZ4pjSKQMIpuw82oUQRBs3a9Lvugv3799O2bVtWrVrFyJEjufvuuyNdspy5kIwiaOCcq2dm6wGcc6neGVkiBV72GNehi5N/111QpVgarVq1YvPmzUyaNIlbb7010uVKmAQTsOneNV0dgJmV4xQrW4kURP7dBRcWPUqTJu3Yvn07s2bNokWLFpEuUcIomIB9E5gKVDCzl4A/AQNCUpVIjPK9oVWJVO7v2o5ff/2VBQsW8Ic//CHS5UmYBTvR4BWgLp4+iI7Ouc0hqUokyuQ2zjWn4/WrJ1Bo/zZatmyPc44lS5ZQt27dSJUuERRMwJbCs9bqL8AE72eRfC+3ca65Hf/ss89o164dpUqVYsGCBVx66aWR/hYkQgJertA597xzrg6exVgqAUvNbGHIKhOJEqday9X/+MKFC2nRogXlypVjxYoVCtcC7kzWg90N/AzsA8rnbTki0cd/LdeGNcrkePzo/z6nXbt21KxZk+XLl1OtWrUIVy6RFsw42L8AtwHlgEnARL8dYaOKxsFKXjpdH+ze9Qt44Yl+XHvttcyZM4eEhIRTvJrEuJCMg60KPOSc+yL4ekTyp/rVE/h85kc892g/mjVrxrRp0yhZsmSky5IoEVDAeqfJvu2c+/G0F4vkM7ndzHLOMXjwYAYMGEDHjh0ZN24cxYoVi3S5EkUC6oP17jQw57QXiuRDOd3Mcs7xt7/9jQEDBnDnnXcyadIkhav8TjA3udaZWVSuOyBytpJSUhm6OPnE5oO+/G9mXVu9NH369OHVV1+lb9++jBo1ikKFwrlBs8SKoNYiAHqY2Q/AYTwdvc45d1UoChMJl9Ot5+o7Oyux6rn86+l+jB8/nqeeeooXX3xRyw1KroIJ2FYhq0Ikgvy7AD5etz3H2Vm1yxejS5cuzJ49m1deeYUnnngiwpVLtAsmYLeRw3qwQEpIKhMJE98FWuLjjMlJ28nIPLk1++uvv3LLLbewbNkyhg0bxn333RfpsiUGBBOwb+NdDxYYBBwEPsazlbdIzPLtAvhp/1HGrd520g2tC0tm0aZNG9avX89HH31Et27dIl2yxAitByvCb+u5JqWk8vG67SeWG7z4nDSaNGnCd999x9SpU2nfvn2kS5UYovVgRXz4tmarFj5E3243s3v3bj755BNuvPHGSJcnMeZM1oMtr/VgJRb5T3fNbfpr/eoJFD+8kxYt2nHs2DEWLVrEddddF8HKJVYFs233R2aWBDTzHurgnPsmNGWJ5C3/oVgD29fJce8sgLVr19K6dWuKFCnC0qVLueKKKyJcvcSqgCcamFkXYIdzbihwPjDYzOqFrDKRPOQ/FGvuxp05LkG4bNkybrrpJkqVKsXy5csVrnJWgpnJ9Yxz7qCZ3YBnJMH7wDuhKUskb/nPxmpzRaXfLUE4Z84cWrVqRZUqVVixYgUXX3xxpMuWGBdMH2ym93M74D3n3GwzezEENYnkOd+bV9l9rpdWLHXicfLK+fTo0YOrr76aTz75hLJly0a6ZMkHgmnB7jCzd4GuwBwzKxrk8zGz1ma2xcySzax/DufNzN70nt+Q3QVhZsXMbLWZfWlmm8zs+WDeVyS3vbP6Nq3J+gUf061bNxo1asSiRYsUrpJngmnB3ga0Bl51zu03s0rA44E+2TvEayjQAtgOrDGzGX6LdrcBank/GuDpgmgAHAducs4dMrPCwAozm+ucWxlE/VJAnWqtgddee43HHnuMNm3aMHnyZEqUKBHhaiU/CaYFmgVcBPzDzD7G05JdFsTzrwOSnXNbnXNpwHigg981HYDRzmMlUNrMKnkfH/JeU9j7EdhWDFKg5LQqVm7LDQ4cOJDHHnuM2267jWnTpilcJc8F04IdDfwKvOV93B34D9AlwOdXBnwX7N6Op3V6umsqAzu9LeAkoCYw1Dm3yv8NzKw30BvQfkgFUG4tVd+1BgoXiuO6CxN48MEHeeutt+jVqxfDhg0jPj4+0uVLPhRMwF7hnKvt83ixmYVtTy7nXCZQ18xKA1PN7Arn3Ea/a4bj2VqcxMREtXALmJxaqtlTYE8sN1jtPN4e9BijRo3ikUce4Z///KeWG5SQCSZg15lZw+x+TzNrAASzq+AOPPt6ZaviPRbUNd7+38V4+oNPClgp2LJbqmnpWZgZCSV+WyqjfvUErqhYgu7duzNlyhQGDRrEgAEDFK4SUqftgzWzr8xsA1Af+MzMfvAuuv05kBjEe60BapnZRd5FYroCM/yumQHc5R1N0BA44JzbaWblvC1XzKw4nhtlmkUmJ6lfPYGB7esQF2dkOcegWZtO9MUePnyYW265hSlTpvD666/zzDPPKFwl5AJpwebJ8kHOuQwzewCYB8QDI51zm8ysj/f8MDz7frUFkoEjQE/v0ysBH3r7YePwbBk+Ky/qkvwl9UgaWc6d1E1w8XlGu3btWLlyJR988AH33HNPpMuUAuK0AeucO7GgtpldDTT2PlzunPsymDdzzs3Bb/NEb7Bmf+2Avjk8bwNwTTDvJflXbou0AL+7oXXJuVk0bdqUTZs2MXHiRDp37hyhqqUgCrgP1sweBP4MTPEeGmNmw51zb53iaSJ5Kpj9sy4sdpQHut/Mjz/+yMyZM2nVSrseSXgFc5PrXjyLbh8GMLNX8PTDKmAlbALdP+vctL00b96e/fv3s2DBAq6//voIVy4FUTABa/y2HgHer3WXQMIqkP2zNmzYQMuWLcnKymLJkiVcc416lyQyggnYD4BVZjbV+7gjnhW1RMLmdPtnpe/cQps2bShZsiQLFizgsssui3TJUoAFs+D2EDNbAtzgPdTTObc+JFWJnEJu+2fZT1/R/L47qVSpEgsXLqR69eqRLlUKOPPcuM9/EhMT3dq1wcyDkGiW28iB7OPHklfx1AP/x6WXXsr8+fOpWLFiBKuVfC7grtFgughEwsY3UIFcRw7Ur57A18tm8eD9Pbn22muZPXs2559/fiRLFzlBAStRx38oVud6VXJcYwDg7bffpm/fvtx0001Mnz6dkiVLRrh6kd8EtWC2SDj4D8Vy8LvtXZxz/P3vf6dv37506NCB2bNnK1wl6pxRC9bMmjrnFnsX3d7tXelKJE/4z8bqXK8KnetVOdFlUK9aafr3788//vEPevTowciRIylcuHCkyxb5nTPtImhlZt/i2aEgBfhr3pUkAp3rVcF5P/v2t2ZmZnL//ffz7rvv8pe//IW33nqLuDj9IibR6UwD9nzgb8ATeGZ4ieSJnPpfs6Wnp3P33Xczbtw4nnzySV566SWtiCVR7Uz/6x8E7MHTclX3gJwx/y1eclo0G+Do0aPceuutjBs3jpdffpnBgwcrXCXqBdWCNbNr8GwVcxvwM3CZc+68UBQm+V9OC7f49782rFGGgwcPcsstt7B06VLeeecd+vTpE+nSRQJy2oA1s0uAbniC9SAwCbjROfe9mX0f4vokH8uptdq3ac0TU2Eb1ijDhSWzaNasGevWrWPMmDF079490mWLBCyQFuw3eHYj+JNz7iu/c/lzGpiERcMaZSgUZ6RnOuLj7MSkguypsDt37uTGG1vy7bffMnXqVG6++eYIVywSnED6YG8Fvgfmm9l/zOxmM9OYGMkb2f2ofv2p33//PTfccAM//PADc+fOVbhKTDptwDrnpjnnuuLZLnsunm2xt5vZB8C5Ia5P8pGcbmhlZHomEmRm/nZDa/PmzTRu3JjU1FQWLlxI06ZNI1i1yJkLZjWtw8BYYKyZJQBdAC1XJAEJ9IZWUlISrVq1onDhwixbtowrrrgi0qWLnLEzGgfrnEsFhns/RE7L94bW8XTPTgSDO1150g2tI9s20r59exISEli4cCE1a9aMdNkiZ0VTYCQssm9ogefO6OSk7SSlpFK/egJ9m9Zk99cradmyJRdccAErVqxQuEq+oICVsKhfPYEuiVVPLKTp2+c6adIkOnToQO3atVm2bBlVqlTJ/YVEYogCVsLm1npVKFr45FWx3n//fbp27UqDBg349NNPKVeuXKTLFMkzWg9W8lxuuw/47qfVsEYZlk0ZxSOPPELr1q35+OOPKVGiRASrFsl7CljJU76jBQrFGV0Sq3Kr34pY9aqV5rnnnmPQoEF06dKFMWPGUKRIkQhXLpL31EUgecp3tEBapmPsqm3cMWLlibGvWVlZPPzwwwwaNIj/+7//Y9y4cQpXybcUsJKnsse2Zt/Mcvy2zkBGRgb33nsvb7zxBg8//DAjRowgPj4+kuWKhJQCVvJUdj9rtwbVTtrmpV7lknTt2pVRo0bx/PPP89prr2m5Qcn31AcrZ8V/99fsrwd3uvLENi9XVyzGc/3uZv78+bz++us8+OCDEa5aJDwUsHLG/G9oYUZG5slba198ntG+fXs+//xzRo4cSc+ePSNdtkjYhLWLwMxam9kWM0s2s/45nDcze9N7foOZ1fMer2pmi83sazPbZGZqAkWBk9ZzzXSk+63tunv3bpo2bcrq1auZMGGCwlUKnLC1YM0sHs8miS2A7cAaM5vhnPva57I2QC3vRwPgHe/nDOBR59w6MysFJJnZAr/nSpj5LtYS723BZmZ6Fm65qPgxmjRpz7Zt25g5cyatWrWKdLkiYRfOLoLrgGTn3FYAMxsPdAB8Q7IDMNo554CVZlbazCo553YCOwGccwfNbDNQ2e+5EkbZfa8D29ch9Uga3+46yMqt+6h6fgm6X1aEv3RtT2pqKvPnz+eGG26IdLkiERHOgK0M/OjzeDue1unprqmMN1wBzOxC4Bpglf8bmFlvPOvVUq1atTwoWXxlh2pCiSIMmrXpxNKDretUZNoXPwGw7bstfPLksxQrZCxevJh69epFuGqRyImpm1xmVhL4GHjIOfer/3nn3IklFBMTE7WdTR7yvaEVZ0aWcyf6W2dt8ITr8R3fsHvyc8QXKcaa5cu57LLLIly1SGSF8ybXDqCqz+Mq3mMBXePdpuZj4CPn3JQQ1ik58L2hlZXliDMj3juMNSMLjqZ8ya4JA4grVoqr+ryucBUhvAG7BqhlZheZWRGgKzDD75oZwF3e0QQNgQPOuZ3mGZH+PrDZOTckjDWLV0KJIsSZEQcUKRzHoA5XcNPlFch0cOTbleye9ByFzqtAhTte4ZFbG0e6XJGoELYuAudchpk9AMwD4oGRzrlNZtbHe34YMAfMmTHMAAATQUlEQVRoCyQDR4DscT3XA3cCX5nZF95jTznn5oSr/oLIv881yzni4oyB7esAsGjzLg5tWsy+2f+iSMWalO/yPJ0bXUr3Bur/FgEwzw37/CcxMdGtXbs20mXErNz6XOMNbrq8Aos27+JA0mx+WTCMYtWvpHynAdzf4gr6t7080qWLhFrAc7xj6iaXhI9vnyvelqvhMIOFm3ex//NJ7F/6IcVrXkf5jv0Z/Kf6armK+NFiL5Ij/z7XXjdcxBWVzyMzy/HL4lHsX/oh59S+kQqdnlK4iuRCLVj5naSU1JP6XFvXqciIFd+TnpnJL/Pf4dAXcylZtw1lW93PS52uVriK5EIBK8DJq2L5dw/M+PInMjMy2DvnXxz5einnNfwTZW68hxc6XqlwFTkFBaycdEOrSKE4BravQ5FCcaSlZ4FBZnoae6a/zNHk1Zx/493c99dH6eyzDYyI5Ex9sMLKrfs4nu7d5iU9i9QjaQxsX4e4OCPz2BF2TXqOo8lrKNvqL7z18iAGd7pS4SoSALVghYQSRcgerJflfbzxpwOkHf6VXZOeJe3nZK7v9SxvDPirglUkCArYAsq3zzX1SBpxBlnOM8BvwpptfPntD/w8dgDpqT9RqfPTCleRM6CALYBO1eeaBSRt+h+7xg8g83AqFbo8R8/bbla4ipwB9cEWQL6jBI6nZ7HxpwMMbF+HamVKkL7vR34e8wRZxw5SsetLlK55DbfWqxLpkkViklqwBVDDGmUoFGekZTocMHHtj0y27Rze8T9+njAQi4ujyp2vcEebxhotIHIWFLAFhG+fa/3qCXRJrMrYVdtwQEam49j2jeye9DzxxUrS7m9DeeaOZgpWkbOkgC0Asvtcj6d79s4a1OEKbq1XhUlrfyQt03F0axJ7pg6m0LnlqHrHSwpXkTyiPtgCIHucqwMyshwDp29ky88HubzSuRz+ZgW7P36BwmWq0OGZ95j4qG5oieQVtWDzMd/1XOPjjIwsz2jXTG/IHvhiHns/+TdFL7iMqt2eo3+nBgpXkTykgM2n/LsF2l9ViVkbdpKZ5Vly8JdV00j99D2KX1SPDo8N4bF2VylcRfKYugjyKf9ugZlf/kS9aqUxHKnLx5L66XuUuOQPVLn9WYWrSIioBZtPNaxR5uRuAQervv+F1E9HcHDtdEpe2ZxbHniOR1rVVriKhIgCNh/x7XNNPZJGrxsu4r3lW8l04LIy2ffJvzn81QJK1b+FSq17K1xFQkwBm0/49rk6IM6gUJwRF2dkpKWxd+Y/ObLlv5zfuDu9H/obf6pfVeEqEmIK2Hwie/rriVWxHKRlOrLSj7Fn6t859n0S9W97kHf/8ayCVSRMFLAxLrtb4ODRdOLMcM79FrLHD7N78vMc3/ENFds/xLv/GKhwFQkjBWwM8+8WMDxdA5XOK0bKT7vYPXEgaXtSaHzfC/yr//0KV5EwU8DGMP9uAYdntMAP27aza8IAMn/dQ5XbBypcRSJEARuDTtUtkJ76E7vGDyDr+CFaPPIGg/9ym8JVJEIUsDEmt24BzDi663t2T3gGl5VJtR4vK1xFIkwzuWJMbt0CR7d/w66x/SEujnZPvsvHz/RQuIpEmAI2hiSlpPLlj/sxv+PHUjbw84QBxBUrRfW7/8mzd7VSuIpEAXURxIiklFS6Df+ctEx30vEjyavYO+1lylxQjYeHjKbVtZcrXEWihAI2Rqzcuo90v3A9/PUS9s4aQrFKNRk/Yw7N6l4coepEJCdh7SIws9ZmtsXMks2sfw7nzcze9J7fYGb1fM6NNLPdZrYxnDVHWlJKKkMXJ3PwaPpJXQMH189h78zXqHBJXebNX6BwFYlCYWvBmlk8MBRoAWwH1pjZDOfc1z6XtQFqeT8aAO94PwOMAv4NjA5XzZGWlJLK7e9+RkbWyccPrJzM/qWjOKdWAz6ePoXrL70gMgWKyCmFs4vgOiDZObcVwMzGAx0A34DtAIx2zjlgpZmVNrNKzrmdzrllZnZhGOuNmOxxrrO+/OmkcHXOsX/Zh/y6cjIXXteS8R/9hwY1y0euUBE5pXAGbGXgR5/H2/mtdXqqayoDOwN5AzPrDfQGqFat2hkXGkn+41yzOZfFLwuGcWj9HErWbcOEcaO5rkbZiNUpIqeXr4ZpOeeGO+cSnXOJ5cqVi3Q5Z8R3J4JsLjODfbOGcGj9HM5t0JmKbfoqXEViQDhbsDuAqj6Pq3iPBXtNvpU9zvWkcM1IY8/0VziavIrSf7yb8xp2odHFCleRWBDOgF0D1DKzi/CEZlegu981M4AHvP2zDYADzrmAugdiXU7jXLPSjrJnygscS9nA+S3up1S9dtQqdw6j7/XvWRGRaBS2gHXOZZjZA8A8IB4Y6ZzbZGZ9vOeHAXOAtkAycATomf18MxsH3AiUNbPtwLPOuffDVX8oJaWkMmjmppPCNfPoQXZPeo60n7+lTPtHKVmnKReWKcGCR2+MXKEiEhTz3LDPfxITE93atWsjXcZp5XRTK/NQKrsmPkP6L9sp16E/JWo1JM5gUp8/aJaWSOT5z1bPlWZyRVBSSiqvL/wfx9J/G4uVcWA3uyY8TeahVMr/6TmKX1iXyyuW4sVOVypcRWKMAjZCcupzTd/3I7smPINLO0qF21+kaOXL6Fj3Al7vek0EKxWRM6WAjYCc+lzTdn3HrokDwYwK3V+maIWLeKnjlXRvEJvjeUVEARt2ObVcj23fxO5JzxNX7Bwq3P4iRcpUVriK5AMK2DDKqeV6dGsSe6YOJv7cclS4/QUKn1dO4SqSTyhgwySnluvhLf9l74xXKVy2GhVuG0ShkqUVriL5iAI2TIYt/e6kcD301UL2zX2TohdcSvk/PUt8sZIKV5F8RgEbBkkpqSzavOvE41/XTid10XsUu/AaynV6moaXVOJvbbQTgUh+o4ANsbGrtvHmov+R5TzLDR74bDwHVnxEiUv+QNmbH+f+my6lf9vLI12miISAAjaEXp6zmWHLtgKecE1d/D4H10zjnCuaU6ZNPzrVq6pwFcnHFLAhkpSSyrvZ4ZqVyS/zhnJow3xK1b+ZhGZ/5v4/1lS4iuRzCtgQyO4WcIDLTGfvzNc4smUF5/2hG+fd0J37/3ixwlWkAFDA5rGHxq9n2hc/AZCVfow90/7Osa1JJDS9l3Ov60THuhcoXEUKCAVsHnp5zubfwvX4YXZPHsTx7V9zfut+lLq6ldYVEClgFLB5ZOyqbSduaGUeOcDuiQNJ2/MDZW95gnMub0yfJjXUchUpYBSweWDsqm08NfUrADIO7mX3hGfIOLCL8rcOoPjF1wIoXEUKIAXsWUpKSWWAN1zTU3eya8IAso7+SvnbBlGs6hUA9GlSI5IlikiEKGDP0stzN5MFpO35gd0TnsFlZVKh298pWrEmgG5qiRRgCtgzlJSSyiMTviDllyMc/2kLuyc9hxUqTIXuL1OkrGc9AfW7ihRsCtgz4NvneixlA7unvEB8ifMof/uLFC5dEVDLVUQUsEHz7XM9kryKPdNepnBCJcrf9gKFSpUB1HIVEQ8FbJBWbt1HFnD466XsnT2EIhVqUL7L88QXP5ci8ca43o20KpaIAArYoK3auo+DX8zll3lvU7RqHcp3Hkhc0RIAClcROYkCNghjV21j5phh7F8yiuIXX0vZDv2JK1wUgMHaVltE/ChgA+Sc4++DBrJ/yYeUuLwJZds9gsV7fnyDO2knAhH5PQVsALKysvjrX//KxjkfUvLq1pzf8n4sLp7LK5biRbVcRSQXCtjTyMjIoGfPnowZM4bHH3+cq2/tyyebfqbNFZXUahWRU1LAnsKxY8fo2rUr06dP56WXXuLJJ5/EzLijYfVIlyYiMUABm4tDhw7RsWNHFi1axL///W/69u0b6ZJEJMYoYHOQmppK27ZtWbNmDaNHj+bOO++MdEkiEoPiwvlmZtbazLaYWbKZ9c/hvJnZm97zG8ysXqDPzSs///wzN954I+vWrWPy5MkKVxE5Y2ELWDOLB4YCbYDaQDczq+13WRuglvejN/BOEM89aykpKTRu3Jjk5GRmz55Nx44d8/otRKQACWcL9jog2Tm31TmXBowHOvhd0wEY7TxWAqXNrFKAzz0rW7ZsoXHjxuzdu5eFCxfSvHnzvHx5ESmAwhmwlYEffR5v9x4L5JpAnouZ9TaztWa29ujRowEX5pyjd+/eHD9+nCVLltCoUaOAnysikpt8dZPLOTccGJ79MNDnmRkfffQRR44c4ZJLLglNcSJS4IQzYHcAVX0eV/EeC+SawgE896xUqVIlL19ORCSsXQRrgFpmdpGZFQG6AjP8rpkB3OUdTdAQOOCc2xngc0VEokrYWrDOuQwzewCYB8QDI51zm8ysj/f8MGAO0BZIBo4APU/13HDVLiJyJsy5gLsqY02+/cZEJKIs0AvDOtFARKQgUcCKiISIAlZEJEQUsCIiIaKAFREJEQWsiEiIKGBFREIkX61F4CfgsWonnmC2ETgWglryWllgb6SLOI1YqBFio85YqBFip85izrkrwvFG+Tlgz8Qx51xipIs4HTNbG+11xkKNEBt1xkKNEFt1huu91EUgIhIiClgRkRBRwJ5s+OkviQqxUGcs1AixUWcs1Aiq83fy82IvIiIRpRasiEiIKGBFREKkQASsmbU2sy1mlmxm/XM4b2b2pvf8BjOrF+hzo6jOkWa22zuWN6TOtE4zq2pmi83sazPbZGYPRmGNxcxstZl96a3x+VDVeDZ1+pyPN7P1ZjYrWus0sx/M7Csz+yKUQ6TOssbSZjbZzL4xs81mljc7nzrn8vUHnh0QvgNqAEWAL4Hafte0BebimZzQEFgV6HOjoU7vuSZAPWBjFP88KwH1vF+XAv4Xip/nWdZoQEnv14WBVUDDaPtZ+px/BBgLzIrGP3PvuR+AstH699J77kOgl/frIkDpvKirILRgrwOSnXNbnXNpwHigg981HYDRzmMlUNrMKgX43GioE+fcMuCXENWWJ3U653Y659Z56z0IbCaH7dcjXKNzzh3yXlPY+xGqO8Fn9WduZlWAdsCIENWXJ3WGyRnXaGbn4WmgvA/gnEtzzu3Pi6IKQsBWBn70ebyd3/+jzu2aQJ6bV86mznDKkzrN7ELgGjwtxLx2VjV6f+3+AtgNLHDOhaLGs64TeB14AsgKUX2B1BDINQ5YaGZJZtY7Cmu8CNgDfODtbhlhZufkRVEFIWAlyphZSeBj4CHn3K+Rrsefcy7TOVcXz/bw15lZWOatB8PM2gO7nXNJka4lADd4f55tgL5m1iTSBfkphKd77R3n3DXAYSBP7rcUhIDdAVT1eVzFeyyQawJ5bl45mzrD6azqNLPCeML1I+fclGisMZv318TFQOsQ1BhQDae45nrgFjP7Ac+vwzeZ2ZgorBPnXPbn3cBUPL/OR1ON24HtPr+pTMYTuGcvlB3P0fCB53+nrXh+Dcju/K7jd007Tu78Xh3oc6OhTp/zFxL6m1xn8/M0YDTwehTXWA7vDQ6gOLAcaB9tdfpdcyOhvcl1Nj/Pc4BSPl9/BrSOphq955YDl3q/fg54NU/qCuVf9Gj5wHP38H947jI+7T3WB+jj/dqAod7zXwGJp3pulNY5DtgJpOP5H/neaKsTuAFPf9wG4AvvR9soq/EqYL23xo3AwGj9M/d5jRsJYcCe5c+zBp6w+xLYFMp/Q2f576cusNb75z4NSMiLmjRVVkQkRApCH6yISEQoYEVEQkQBKyISIgpYEZEQUcCKiISIAlZEJEQUsCIiIaKAlZhjZh3NzJnZZQFeX9rM/uL9+sJwrJkrAgpYiU3d8My66Rbg9aWBv4SunBOLOevfk5xEfyEkpnhX4roR6EXgAfsycLF3CcJXgXgze8+7Y8F8Myvu8/o9vDsafGFm75pZvPf4I2a20fvxkPfYhd4V9EfjmVb7fvY57/mXQrlrg0Q/BazEmg7AQufcl8AhM6sfwHP6A985z5J5jwO1gKHOuTrAfqAzgJldDtwOXO+9NhO4w/sePYEGeBYJ+bOZXeN97VrA297Xeh64y/tacUBXIFQrXEkMKBTpAkSC1A14z/v1RO/jYNdE/d4594X36yQ8q5ABNAPqA2vMDDyrae0GEoCpzrnDAGY2BWgMzABSnGd1fJxzP5jZPm/4VgDWO+f2Bf0dSr6hgJWYYWbn42lFdvYemggsNbM3gIFAKjDPObfoNC913OfrTDxBCp7Vlj50zj3p976n+jX/sN/jEcA9QEVg5GnqkHxOXQQSS/4EzHHOHQdwzm3Fs0Tj/wFpwL9zCdeDeDZZPJ1FwJ/MrDx4At3MquNZK7SjmZXwbiXSyXssJ1PxLNB9LTAv4O9M8iW1YCWWdAOu9q7in60Mnl/z3wD+bWb3O+8K+tmcc/vM7L/e4Vmbc3tx59zXZjYAmO/tQ00H+jrnVprZKGC199IRzrn13n3F/F8jzcwWA/udc5ln+o1K/qD1YCXmmdkreLZtLgY87JxLj2AtccA6oItz7ttI1SHRQQErkkfMrDYwC88NsUcjXY9EngJWRCREdJNLRCREFLAiIiGigBURCREFrIhIiChgRURCRAErIhIiClgRkRD5f1AI1uqp320DAAAAAElFTkSuQmCC",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f658d174110>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,1,figsize = (5,5))\n",
    "ax.plot(ce.reshape((-1,1)), ce_theory.reshape((-1,1)), '.')\n",
    "ax.plot([0,0.06], [0, 0.06], 'k')\n",
    "ax.set_xlabel('$\\Delta_s$ theory')\n",
    "ax.set_ylabel('$\\Delta_s$ observed')\n",
    "sns.despine(trim = True, ax = ax)\n",
    "plt.savefig('./fig_s1_d.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot the deltas surface as a function of DeltaT and tau_s\n",
    "N = 51\n",
    "M = 51\n",
    "deltas = np.linspace(10, 100, M)\n",
    "#tauss = np.logspace(-2, -1, N)\n",
    "tauss = np.linspace(0.01, 0.1, N)\n",
    "ce_theory = np.zeros((M,N))\n",
    "\n",
    "for idx, DeltaT in enumerate(deltas):\n",
    "    for idx2, tau_s in enumerate(tauss):\n",
    "        ce_theory[idx,idx2] = 1./DeltaT*(1 - np.exp(-DeltaT*dt/(tau_s)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVkAAAFCCAYAAABb466DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJztnX20XWV95z/f+5Lc8BpowEaSKVGjrowjECPE0lGR2kmQMe3q1IFWQWqboqDAaC3Yqa2z2hmrjlWmLNIUsFAoLIpoUxpFitouXQ0S3gIhICmKBKMhIK95vff+5o+9T7Jzcl72Pmfvc/Z+9u+z1l6cs/fz7PNscu/3fs/v+T2/R2aG4ziOUwwjwx6A4zhOyLjIOo7jFIiLrOM4ToG4yDqO4xSIi6zjOE6BuMg6juMUSOEiK2mZpEclbZZ0aYvrr5f0b5J2S/pYlr6O4zhlR0XmyUoaBb4PvBPYAtwNnG1mDyfaHAv8AvCrwM/M7HNp+zqO45Sdop3sycBmM3vczPYANwErkg3MbJuZ3Q3szdrXcRyn7BQtsscBTybeb4nPFd3XcRynFIwNewD9ImklsBJglLE3HTp21JBHFBjT08MeQWnxJekH8qI9u93Mjsna752nzbJnns32c3bfhj23m9myrJ81DIoW2aeA+Yn38+JzufU1s9XAaoAjx4+1X5zzG72N1GmL7dg57CGUkumdu4Y9hFJxx94bn+il3zPPTvPdr78yU59DXvnDOb181jAoOlxwN7BQ0gJJM4CzgDUD6OvkiA6ZNewhlJKRWRPDHoJTAQp1smY2KelC4HZgFLjGzDZKOj++vkrSzwPrgSOAaUkXA4vM7IVWfYscr9MeHTLLHW0LRmZNuKN1OlJ4TNbM1gJrm86tSrz+CVEoIFVfZ3i40DpOdnzFl5MJDx0cjIcNnE64yDqZcaE9GBdapx0usk5PuNA6TjoqnyfrDA+P0R6IT4L1xhTG89O7hz2MwnAn6zg54mEDpxkXWacvPGzgOJ1xkXX6xoX2QNzNOklcZJ1ccKE9EBdap4GLrJMbLrSOczBhZReMCCZmZuuzK9xZzWHgGQf78WwDB0IT2V7oJsouwplxod2PC213pgyendawh1EYLrLdaBZhF91UuNA6ToTHZLMyMfPAw2mLx2gjfBJsOKTYxFWSLo+vb5C0OHFttqRbJD0iaZOkt8Tnj5Z0h6TH4v923SXARbZfXHCdFLjQDpZ4I9YrgOXAIuBsSYuami0HFsbHSuDKxLUvAl83s9cDJwCb4vOXAnea2ULgzvh9R1xk88QF9yDczTpDIs1GrCuA6yxiHTBb0lxJRwJvBa4GMLM9ZvZcos+18etriXbZ7oiLbFG42O7DhTbC3exASbMRa7s2C4CngS9Juk/SVZIOjdu8wsy2xq9/Aryi20B84qtoGkJb8wkznwiL8GyDg5lEPDed2ZDMkbQ+8X51vN9fHowBi4EPm9ldkr5IFBb4o2QjMzNJXXfTdJEdFC62LrROnmw3syUdrqfZiLVdGwO2mNld8flb2B97/amkuWa2VdJcYFu3gXq4YNB4GKH2eNhgIKTZiHUNcE6cZbAUeN7MtsZbYj0p6XVxu9OBhxN9zo1fnwv8Q7eBuJMdFjV1tu5mnUGQZhNXov0DzwA2AzuA8xK3+DBwQyzQjyeufRq4WdIHgCeA93Qbi4vssKmh2LrQemx2EKTYxNWAC9r0vR84KBxhZs8QOdvUeLigLNQshOAZBx42qAsusmWiZvFaF1qnDni4oIxMzKxV+KDOeNgApmyUZ6YOG/YwCsOdbFmpiat1N+thg9BxkS07LrSOU2lcZKuAC23wuJsNFxfZqlCD8EHdhdYJExfZqhG40NYZd7Nh4iJbRQIW2rq7WRfa8PAUrqoScJqXrwirF5OM8KyncDmlJOA4bZ0drbvZsHCRDYFAhdZxQsBFNhQCFFp3s04IuMiGhAttULjQhoGLbGi40DpOqXCRDZEAhbauuJutPp7CFSqBpXh5Wle4TNkIz00dMuxhFIY72ZAJzNHWNWzgbrbauMiGTmBCW1dcaKuLi2wdCEho6+pmnexIWibpUUmbJV3a4rokXR5f3yBpceLaDyU9KOl+Setb9P2oJJM0p9s4XGTrggtt5XE3mx5Jo8AVwHJgEXC2pEVNzZYDC+NjJXBl0/XTzOxEMztgQ0VJ84FfAX6UZiwusnXChdapDycDm83scTPbA9wErGhqswK4ziLWAbMlzU1x778APg5YmoG4yNaNgIS2jribTc1xwJOJ91vic2nbGPDPku6RtLLRQNIK4CkzeyDtQDyFq44Ekt5V17Su0DZfnLQRfjZ5aNZuc5pipavNbHWOw/olM3tK0rHAHZIeAdYDnyAKFaQmLJEdGYE0f+kD+gHtGRdap9psb46VNvEUMD/xfl58LlUbM2v8d5ukrxCFH34GLAAekNRof6+kk83sJ+0GUs9wwayJ1kfd8NBBZfGwQVfuBhZKWiBpBnAWsKapzRrgnDjLYCnwvJltlXSopMMBJB1K5FwfMrMHzexYMzvezI4nCi8s7iSwEJqT7ZdWP7ihu94AHK27WacZM5uUdCFwOzAKXGNmGyWdH19fBawFzgA2AzuA8+LurwC+ErvVMeDvzOzrvY7FRbYbSeENVXBdaCtJaLHZvDGztURCmjy3KvHagAta9HscOCHF/Y9PMw4X2SyELLgBCK3jlJF6xmTzIMQ4bsVjtHXMnfXYbPlxJ9svjR/y0JxtRfGwQfWYshGenwz3D2ThTrbP9cOXSNoo6SFJN0oq75/tUJxtxd2s45SNQkW2n/XDko4DPgIsMbM3EM0QnlXkeHMhBLGtuNB62MApE0U72X7XD48BsySNAYcAPy54vPlRdbF1oXWcXChaZHtePxyvuPgcUaWbrUSJwt9o/gBJKyWtl7R+z+SOXAefCy60zoBwN1tOSptdIOkoIpe7AHglcKik9za3M7PVZrbEzJbMGCvpFhZVdrUVFlp3s04ZKFpk+1k//MvAD8zsaTPbC9wK/GKBYy2eKottRamb0LqbLR9Fp3DtWz9MJJxnAb/Z1GYNcKGkm4BT2L9++EfAUkmHADuB04mq4FSfWRPVSvnyhQqVomopXVM2wnN7w/1jWKjI9rN+2MzuknQLcC8wCdwH5FnKbLi40A6MOubOOuWh8MUIva4fjq/9MfDHhQ5wmFRtIYMLbWWompsNmdJOfNWKKsXRKjwR5jjDwEW2LFRJaCuKT4I5w8BFtkxUJfvA3azjpMZFtoy40BaGu1ln0HgVrrJSheyDik6E+SRYuZi2EV7aG+4fA3eyZaYKLqSijtZxBoWLbNmpgtBWEA8bOIPCRbYKlP0XpKJutm5CWzf6qWUdXx+VdJ+k2xLnTpS0TtL9cWGqk7uNw0W2KrjQOn1SJzfbTy3rBBcBm5rOfQb4lJmdCHwyft8RF9kqUfZfkgoKrbvZYOmrlrWkecC7gKua+hhwRPz6SFLUuPbsgqpRhawDp7SUPdMgR1rVqT4lRZvjiOpXfwH4OHB4U5+LgdslfY7IpHatDOgiW0XKLLQVTOuqW0pX2Zgy8cLezN+C5khKVuVbbWa5FJCSdCawzczukfT2pssfBC4xsy9Leg9wNVFZ1raEJbIjYnrWeOtLO/cOeDAF40KbK3US2kDc7HYzW9Lhej+1rH8deLekM4AJ4AhJ15vZe4FziWK1AH/PweGEg6hNTHZ61njLo9KUPUbrlJYaTILtq2UtaQZRLes1TW3WAOfEWQZLiWtZm9llZjbPzI6P+30zFliIYrBvi1+/A3is20DCcrI90Cy0lXO8ZXW07madIdJPLesu/C7wxXhz111EWQkdqb3INpMU3coIrgttbtRJaAMJG7Sln1rWiTbfBr6deP8d4E1ZxlGbcEEvVCqkUNavfxVM63KcPHGRTUFlYrhlFdqKUafc2RrEZoeOhwsy0hDa0oYSyhg6qGDYwBkcUyZe2hPuNx53sj1SamdbRndSsbCBu1knL1xk+8SFNlzqJLROcbjI5kBpXW3ZhLZibrZOuJstDhfZHCml0JaNigmtu1mnX1xkc6Z0rtYdipMSd7PF4CJbEC60HXA369QIT+EqkFKle5UttcvTukrJMFaBTZvYubdEpiRn3MkOgNK42rI52gpRJzfrYYN8cZEdEC60LfCwgVMDXGQHSGmEtkxUTGjrgrvZ/HCRHTClyD7wX6CecTfrZMVFdki40CZwN1tK3M3mg4vsEHGhTVAhoXU362TBU7iGzPSs8XKkeDmZqEtx70GkdJmJnbvDna9wJ1sChupo3c06TqG4yJYEF9qYCgltXcIGHpvtj6BE1iSmZ463PcqOC63j5IekZZIelbRZ0qUtrkvS5fH1DZIWN10flXSfpNsS5z4r6ZG4/Vckze42jqBEthtVEF0XWtzNlpCquVlJo8AVwHJgEXC2pEVNzZYDC+NjJXBl0/WLgE1N5+4A3mBmbwS+D1zWbSy1Etlmyiq4Q886KAMVElqnlJwMbDazx81sD3ATsKKpzQrgOotYB8yWNBdA0jzgXcBVyQ5m9g0zm4zfrgPmdRtIrUU2SdkEd2hCWzHHUgbczZaS44AnE++3xOfStvkC8HFgusNn/DbwtW4D8RSuFjSEdmT3cFOrhpbeVZaKXRWq1FWXlK4iMBN7945m7TZH0vrE+9VmtjqP8Ug6E9hmZvdIenubNn8ITAI3dLufi2wHyiC2nkfrlIlhlEJsw3YzW9Lh+lPA/MT7efG5NG1+HXi3pDOACeAISdeb2XsBJL0fOBM43cys20A9XJCCYYcRhhI6KMtXwwrFZusSNqgIdwMLJS2QNAM4C1jT1GYNcE6cZbAUeN7MtprZZWY2z8yOj/t9MyGwy4jCCO82sx1pBuJONgPDdLZDcbQeNnBaUCI32xYzm5R0IXA7MApcY2YbJZ0fX18FrAXOADYDO4DzUtz6L4GZwB2SANaZ2fmdOrjI9sD0zHEXWqclHpstD2a2lkhIk+dWJV4bcEGXe3wb+Hbi/WuyjsPDBT0yrBBCbdO7KhQ2qAMVyzQYKqmcrKQR4ATglcBO4CEz21bkwKrCsFztQHE3mwl3s06SjiIr6dXAHwC/DDwGPE002/ZaSTuAvwKuNbNOuWTBM+hYbW3DBh6bLRV5xWZtWkzuDjdy2e3J/pRoqdnvNacqSDoW+E3gfcC1xQyvWgzS1dY2tasiQlsXN1uFSbBh01FkzezsDte2Ea2KcBIELbRlcLMVoi5C63Qm1cSXpN+QdHj8+o8k3dpcscbZzyAnxAY+EVaGCQ+fBCsVPgnWmbTZBX9kZi9K+iXgdOBqDq5Y4yQIWmid1PgCBSetyE7F/30X0RrhfwJmpOnYT01HSbMl3RLXb9wk6S0px1sKhr1SrDDK4FzczZYKd7PtSSuyT0n6K+C/A2slzUzTN4eajl8Evm5mrydKIWuu7VgJBiG0HjYoL+5m603avIn3AMuAz5nZc3HNxd9P0W9fTUcASY2ajg8n2uyr6Qisi93rXKJlbm8F3g8Q14Tc0/HTRsT0xP5qPiO7pjo0HiyDmBCrbcaBU20MbG+466I6PpmkwwDMbIeZ3Wpmj8Xvt5rZN5Jt2tBPTccFRHm5X4q3gLhK0qEpnmkf0xOjBx3DJDhH6242Ne5m60u3Px//IOn/SnprUuAkvUrSByTdTuRwi2AMWAxcaWYnAS8DrWK6KyWtl7R+z56Xu9502ILrQus49aKjyJrZ6cCdwO8BGyU9L+kZ4Hrg54FzzeyWDrfop6bjFmCLmd0Vn7+FSHSbx7jazJaY2ZIZMzIZ3aGJbZCTYcPE3axTYrrGZFtVssnAvpqORMJ5FtEqsSRrgAvjeO0pxDUdASQ9Kel1ZvYoUerYwxRAQ2gHGcMtOkY70PhsGRYpVGQlmFM/Cl0wnENNxw8DN8RFdx8nXb3Hnhm02LrQ1g9fBVY/Cq/K0E9NRzO7H+i0xUQhDFJsgxLaYeNu1ikh4Za+yYFBiW0w5RLdzabC3WwTJrS7pilcDSS9Ol6AgKS3S/qIpNnFDq08DGJyrMjJsFplG/gkmFMy0v75+DIwJek1wGqibIC/K2xUJaQMebb94PUNHGc4pBXZaTObBH4N+H9m9vvA3OKGVV6KFNpgUrvczabC3Wyx9Fo3RdKEpO9JekDSRkmfaur34bieykZJn+k2jrQx2b2SzgbOBf5rfC4QRchOkbHaIuOzPgnm1IVE3ZR3EuXc3y1pjZkl00CTdVNOIaqbcgqwG3iHmb0kaRz4jqSvmdk6SacRlQI4wcx2x5sXdCStkz0PeAvwZ2b2gzjv9W9T9g2WolxtEPHZYbvZiuButjD21U2J65406qYk2Vc3xczWAbMlzY3fvxS3GY+Pxs4wHwQ+bWa7Yd/mBR1JJbJm9rCZfcTMbozf/8DM/rxxXdKX09wnRFxoOzBsoa1I2MAphH7qpiBpVNL9wDbgjsTK09cC/1nSXZL+RdKbuw0krxSuV+V0n0oyPTFaudCBUx5qn9JlMLJXWXvNkbQ+8X61ma3ObUhmU8CJcRbVVyS9wcweItLMo4GlwJuBmyW9qnkPxCR5iWzbDxgklvnfKT+qJrQDi88OO3fWY7Ohst3MOi1U6qduyj7i0q7fIiqE9RCR2701FtXvSZoG5hBVDGxJcBnAUzNHDzoGRdVSvGoTNqgAHpvNnX11U+Jl+WcR1UlJsgY4J84yWEpcN0XSMY11AJJmEU2ePRL3+SpwWnzttUQ7xGzvNJC8nOwQPWR3kkI7urvg1VsFOFoPG/SJu9na0WfdlLnAtXGGwghws5ndFl+7BrhG0kNEmwic2ylUAPmJ7B/kdJ/CaQhukWJbJaH1sEF5qH1sNmd6rZtiZhuAk9rccw/w3izjSCWykhYC/4don6593/3M7FXxf7+R5UPLQNFiW0QubeWF1nFqSNqY7JeIEnUnieIR1xEV7q48g47b1pphx2YrkNLlsdnwSCuys8zsTkBm9oSZ/QnR9uDBUJTY5j0ZVlT+rE+COcNC0zCyW5mOKpFWZHdLGgEek3ShpF8DOm2gWFmKEFsX2hLhbtYZMGlF9iLgEOAjwJuIAr/nFDWoMlBXoR0I7madGpFWZI83s5fMbIuZnWdmvw78hyIHVgbydrVVyKN1N1sO3M2GQ1qRvSzluSApq9BWOmzgbtapCR1TuCQtJ0rWPU7S5YlLRxBlGtSGqZmjuaV75ZlHW+mFCsPMnfW8WWdAdHOyPwbuAXbF/20ca4D/UuzQykedHG0twgaOMwA6OlkzewB4QNL18c4ItSfPRQxFFZXJi4EsUnA325FauFmjcmlZWegWLniQuMKWdPD/BDN7YzHD6pERMTWx35yP7pou7KPyCh/kJbSVDhsMkwoIrVNtui2rPXMgoyiIogU3zzhtHhQhtMG72QpQCzcbMB1jsvHqrifM7In41ML49Tbg2cJHlyNTEyMHiG5u980hTluF1K6gqUBKl1NdUqmOpN8FbgH+Kj41j6iuYuUoQmzLJLSVnQTzlK6OeN5sdUmrNhcApwIvAJjZY0DXXRrLTN5i60KbA8MU2gq4WRfaapK6dkFcRxEASWOUZMuZfslTbEMXWsdxspNWXf5F0ieAWZLeCfw98I/FDWvwlEloy4q72eETpJs1GNmb7agSaZXlUqKNwh4Efo+o2vj/LGpQwyIvV9uv0LqbdZxwSKUoZjZNNNH1ITP7b2b21932takyLrQd7ududugE6WYLQNIySY9K2izp0hbXJeny+PoGSYvj8xOSvifpAUkbJX0q0edoSXdIeiz+71HdxtFRTeJB/Imk7cCjwKOSnpb0yeyPXC1CEtq88SW3TtmJN0G8AlhOtG3W2ZIWNTVbDiyMj5VEu78A7AbeYWYnACcCy+LdbCH6Vn+nmS0E7ozfd6SbklxClFXwZjM72syOBk4BTpV0SbebV52icmsHTSXDBu5mO+JutisnA5vN7PF40v4mYEVTmxXAdRaxDpgtaW78/qW4zXh8WKLPtfHra4Ff7TaQbgryPuBsM/tB44SZPU5Ji3abYGpG/mug+xFad7OOMxSOA55MvN8Sn0vVRtKopPuJFl7dYWZ3xW1eYWZb49c/AV7RbSDd1GPczLY3nzSzp4nUvZRMzdC+I7d7Vlxo3c1mxN1s2ZkjaX3iWJnnzc1sysxOJFp4dbKkN7RoY6RIZe1Wu2BPj9dKQ0NoR/f0P083NTHScw2EMtQ5yLu2gdc1cPJABqPZa/RsN7MlHa4/BcxPvJ8Xn8vUxsyek/QtYBnwEPDTOKSwVdJcIqfbkW727ARJL7Q4XgT+U7ebl4m8nO2wHK2HDYaAu9kqczewUNICSTOAs4jqYCdZA5wTT/AvBZ6PxfMYSbMBJM0C3gk8kuhzbvz6XOAfug2kWz3Zcv5m90EeznZYjjaPsoiVLInobtbJiJlNSroQuB0YBa4xs42Szo+vryLK9z8D2AzsAM6Lu88Fro0zFEaAm83stvjap4GbJX0AeAJ4T7exdAsXBEu/YtuP0PZDGYV2IGGDYeH1ZiuLma0lEtLkuVWJ10ZUl6W53wbgpDb3fAY4Pcs4qp+f1Cf9hBB6DR2EvPS2MHwSrC0eMig3tRdZ6C9eOwyhLWO2QdCxWcfpAxfZBC60JcfdbFvczZaX2sZk2zE1Qz3FaYcVoy0TQcdmneKYhpFKJIT2RlBO1kZgciKHNK0BOtrQ3GzhYQN3s21xN1tOghLZBpMT2nf0Sp2E1nGc4ghSZJP0I7ZF1EFo+1lDzDhwN5sBd7NORoIX2Qb9CG1WsR105a4yhg0cx4mojcjCYF1t3cMG7maHh7vZclErkW3Qq9iWXWj7pXJu1rcRdypArVO4JifE2K5id9EZZGpXHktu8yTolK6SL7fVIbOwHTuHPYxUyGDUU7jCJaujHcRkWEhuNuiwgeOkoHCR7XUzs8T1UUn3SbqtuW9eFC20gwwblC02GzQem3VSUKjI9rmZWYOLgE3pPhCmZvbmNLPGaUMWWnezjpMfRTvZnjczA5A0D3gXcFWWD52aqX1HVsomtMOicpNgw8LdrNOFon/r+9rMDPgC8HGg55mjsgltVkIJG7ibdepKaa2VpDOBbWZ2T5d2KxubqU3uerllm15cbVFCW+ewQbC4m3U6ULTI9rOZ2anAuyX9kCjM8A5J1zd/gJmtNrMlZrZkbOLQjoOpstCGgLtZpxWNFK4sR5Uo+re9583MzOwyM5tnZsfH/b5pZu/td0BZXW0eVb1ajiOj0LqbLTnuZp02FCqyZjYJNDYz20S0IdlGSec3NjQj2oPncaLNzP4a+FCRY2pQhNAOsqBMVXE36wyKXtNHJc2X9C1JD0vaKOmiFn0/Kskkzek2jsJXfPW6mVlT+28D3857bFMzxejudCu+0q4Oy1L0O+tqsF53uu13JVgld7gdBr4KrDQk0kffSTSZfrekNWb2cKJZMn30FKL00VOASeCjZnavpMOBeyTd0egraT7wK8CP0owlqOCgCaZmZOuTJXxQhKMdVNigX/IMG7ibdQZAz+mjcbjyXgAze5HoW3gyK+oviLKeUrmpoES2wdSM3sQ23zGUK3RQtpSuYPHYbFnoN30UAEnHE20Pflf8fgXwlJk9kHYgQYpsg6xim0Zoi5gIq+MkmLtZp0/mNFI342Nl3h8g6TDgy8DFZvaCpEOATwCfzHKfWlThmpqRPu0jTZy2yvHZMhFslS6PzWbDSD03kmC7mS3pcL2f9FEkjRMJ7A1mdmt8/dXAAuABSY3290o62cx+0m4gQTvZJFlcbZ6ONrSwQaVSutzN1pme00cVKejVwCYz+3yjsZk9aGbHmtnxcWrpFmBxJ4GFGolsg2EIbVqqMgmWJ4WHDZxa0mf66KnA+4gWQN0fH2f0OpZahAuayRI+6Eaa0EGRYYNeqFVK16wJ2Llr8J/rIYOh02v6qJl9B+jqoGI325WgnGyWFK407fLMOCgqbBBCARl3s07IBCWyDdLGX9O0G0Z8tgphA4/NpsDTuRwCFdkGebnaKsRne6FMk2DuZp1QCT4m2xDQbjHYbnHavFK7ssRnsxBCSleheGy2JWWIzWqawjc0HSZBO9kkacMHHa/nFKMtKmzQC+5mHadYaiOykH2pbct7dBHaokojpiGElK5C8dhsSzw2WyxhiWyKp+nqVnMQ4jRCW6ZJMHezjlMcYYksMD0jOjrRLaugjmGD2uButiXuZosj2N/ebkILxQpt1cIG7mYdpxiCFVlI72p7uQb5CG2Z3GyZFigUirvZlribLYbgU7ggEtqRTulZHdK38lyCO0iGkdKV53LbYCt0OQcRbaToKVyVJ034oB0d3a672WribrYl7mbzJyiRNcFUh5/hTkLbT1ZBWYW26sttPTbrhEBQItugCKEdVMZBGXA3WzAld7NOvgQpsjAcoe2Eu9ke7+VuduB4yCBfghVZiIS2ndgWIbRlTesaRkpXZXA36xRM0CLbIG+h7fxZ/QlpmSbB+sHdbLVxN5sfYaVwddCnqZkw2qIYUqf0rnbpW/2kdQ2rUlcvKV397qBQGbxC13CZtsJ3Axkm5bZDPTDdwfT04mjb3qsEYQN3s06RuJvNh3L/lvbI9Hh7sc0qtEXFZ7tRlkkwj806VUXSMkmPStos6dIW1yXp8vj6BkmLE9eukbRN0kNNfU6UtC7eXHG9pJO7jSNIkW0wCKHtlWHWNhgk7ma7UPIJsKq6WUmjwBXAcmARcLakRU3NlgML42MlcGXi2t8Ay1rc+jPAp8zsROCT8fuOBC2yULzQupsNCHezIXEysNnMHjezPcBNwIqmNiuA6yxiHTBb0lwAM/tX4NkW9zXgiPj1kcCPuw0kLJFto0lZhTYrvTrdPHNny4y72S64my2C44AnE++3xOeytmnmYuCzkp4EPgdc1m0gYYksMD2z9ax8FqHNM2xQxkkwd7NOAMyJY6KNY+WAPveDwCVmNh+4BLi6W4ewUrhiGkI7svtAAZseh5GUhZ3apXb1ktaVZhPGThS1+eIgqUyFLk/nGjgyeqkYt93MlnS4/hQwP/F+Xnwua5tmzgUuil//PXBVt4EG52STtHK1rRxtXqldRYYNUn2+u1mnICoYMrgbWChpgaQZwFnAmqY2a4Bz4iwvXtnjAAAPJ0lEQVSDpcDzZra1y31/DLwtfv0O4LFuAwnSySaZnmmpHG27xQqt6GUxgrtZd7NdqbGbzRszm5R0IXA7MApcY2YbJZ0fX18FrAXOADYDO4DzGv0l3Qi8nSgssQX4YzO7Gvhd4IuSxoBdRFkJHQlKZG2kjaj2IbSDChukWQmWhqmJkdSrZ3wVmJMFHTIL27Fz2MNIjZmtJRLS5LlVidcGXNCm79ltzn8HeFOWcQQZLmgZJugjdNBPwe8sdAsbeKZB072KzDTwwjFOTgQpspBeaPthGLmzafDYrFMUFYzNDp1gRRbSCW2/brYXoe1EXm62zDVn83SzheJu1smBoGKyrRYjtIrHHtSmz/hsVkKYBCtLbNY3XBw8ucdmp60UP0tFEZyTnR5P4V5TxmfTMmg3mxZ3szngS22dPglOZKF3oW2myNVgg6prUCRlic36UtvB47HZ9AQpstCb0BaRbdDrJJi72RLhbtbpg7BEVs2i2d2tphHatOQdNuj+ee5mG7ibdcpKWCIL2MwDE/GbhTZNmKBZaIvMnXU3WxFxdDd7EB4ySEdwIgvZhTbP+GxZU7qKxN1sgbibrTxhpXAlsJnTaPf+vyHT48bI3v2C1JzaddD7DBW7msla26DflK5Un1Hwctt+yLOmgTNY8kjnklnQ//5hOVmBxvcLSb+hg7zDBu5mi8eX2jplIyyRjckitAdcG0DYoP19q7/cth8qE5t1DsJjs50JUmSbaRbaJL2kdfWDu9nicTfrlImgRFYj+wUy6WbhQKHNOhHWa9gg6ySYu1l3s1XF3Wx7ghJZgLGZk/teNwttkjQ5tJ3IaxPGg+/be0qXu9n9BOlmnUoSnMhCe6HtFJ/NI2wwiJSuPMjiZp0S4SGDSlL4b5ukZZIelbRZ0qUtrkvS5fH1DZIWx+fnS/qWpIclbZR00cF3b09aoU2Std5sWjeb5yTYoN3soOvNelHvGjJtjOzcm+moEoWKrKRR4ApgObAIOFvSoqZmy4GF8bESuDI+Pwl81MwWAUuBC1r0bfo8Y3x8f35nGqHNkm3Qq5tth7tZJzPuZitH0b9pJwObzexxM9sD3ASsaGqzArjOItYBsyXNNbOtZnYvgJm9CGwCjitikP2EDYaR0uVutiS4my01vX6Ljq9dI2mbpIea+nxW0iNx+69Imt1tHEWL7HHAk4n3WzhYKLu2kXQ8cBJwV5oPzepmm8l7m5pOuJstFl9qW0/6/BYN8DfAsha3vgN4g5m9Efg+cFm3sZT+t0zSYcCXgYvN7IUW11dKWi9p/dQLLzNrZhSvGWTYwN1sa2rhZp2y0vO3aAAz+1fg2eabmtk3zKwhKOuAed0GUrTIPgXMT7yfF59L1UbSOJHA3mBmt7b6ADNbbWZLzGzJ2JGHHHAtjdAm6RQ26EY/KV3uZoslyAkwd7PdyOVbdBd+G/hat0ZF/4bdDSyUtEDSDOAsYE1TmzXAOXF8ZCnwvJltlSTgamCTmX0+y4c23CwcKLStSBs2KDqlqx3uZh2nJXMa32DjY+UgP1zSHxJNzt/QrW2hVbjMbFLShcDtwChwjZltlHR+fH0VsBY4A9gM7ADOi7ufCrwPeFDS/fG5T5jZ2nafNyJj1vhedu4dZ9bMvezcfeAv6djMSSZ3R4+s8Wlsb/Q3Jlmxq7laV5JumzK22nwxLVkrd+VJlgpdVaXQDRdnTcDOXcXcuw5MT/fy/2+7mS3pcL2vb9GdkPR+4EzgdDPr+pW38FKHsSiubTq3KvHagAta9PsOLfefTU9DaMfHp9i7N3JUSaFNQydhTVMOsdUOt3mWQpycEGO72v87572zrZdBLBETM2FXj3/Vw2fft2gi4TwL+M2mNmuACyXdBJxC/C26000lLQM+DrzNzHakGUhQAbnRePuZWeP7fyk7TYQNahIsLaHGZn2prTNo4smpxrfoTcDNjW/RjW/SRObvcaJv0X8NfKjRX9KNwL8Br5O0RdIH4kt/CRwO3CHpfkn7DGM7givafdiM3by0Z+a+sEGSpKNtkCZs0C1M0A13s/3hbrYD7mbb0uu36Pja2W3OvybrOIJysg0OmxH90DUcbXIirEEy2yArobrZInE369SVoER2VO3dWquwQYM0YYN+FyhUIdOgzGUQnQ54OlepCUpkAY4Yj1xss5tN0hDatGURG2RJ6XI3ezCezuW0ZNqikEeWo0IEJ7LQXmjThg16nQTrhrvZchBkyMDdbGkJSmRHNM1h453z7XoJGyTpJKzuZrvjbtapG0GJbJI0YYMG3SbB3M12Gou72YPwCTAnQVAiO6rIfTbcbLewQRGTYO5mu+NutiA8ZFBKghJZgNnjOwG6hg0a9DoJlqTubraquJt1BkFQv0kNJ9sQWihuEszd7H58V9sS4W62dAQlsgBHju0X2LzCBmmomptNdX9fatsf7mbTMT2N7diZ6agSwYlsg6xhgwZ1cbNlKIPYD+5mO+ButlQEJbJjmuaosZf3udksYYMQ3Wwe1MHNOk6RBCWySTqFDdrRahKsQRXdbBWKepcBDxk4RRKUyDYmvo4ae3nfuaSbbZBmEqxBFdxslfF0roLwkEFpCEpkAWaP7q+jm2YSrJkqutkypXO5m22Bu9laE5TIjhGncI3u2Odmjxzb2XESzN3s8HE3WxDuZktBUCILcPToS/tetwsbpJ0EczebaOOLE/rD3WxbzIzpnbsyHVUiqN+cUUXiePToS13DBllJUzwmSZ3drO9qWyLczQ6doEQW4OcSTjYZNoD83GwSd7OO43QiqN+aMYzZI7v5udGXDggbQG8pXc14bDY9PgHWAg8ZDBRJyyQ9KmmzpEtbXJeky+PrGyQtTtNX0oclPSJpo6TPdBtHUCI72sKMdXKzDTw2m+FzS7o4wUMGHahhyEDSKHAFsBxYBJwtaVFTs+XAwvhYCVzZra+k04AVwAlm9h+Bz3UbS1AiC3D0SGs3m1wJBu5mW+GLEyLczQbBycBmM3vczPYANxGJY5IVwHUWsQ6YLWlul74fBD5tZrsBzGxbt4EEJ7LNNE+CQeu6BlV0s+0oeqltkbibLYj6udnjgCcT77fE59K06dT3tcB/lnSXpH+R9OZuAxnLOPBSs2HD3u1z5219ouCPmQNsL/gzBkEIzxHCM0A4z/G6Xjq9aM/efsfeG+dk7DYhaX3i/WozW93L52dkDDgaWAq8GbhZ0qvMrK2zCkpkzeyYoj9D0nozW1L05xRNCM8RwjNAWM/RSz8zW5b3WICngPmJ9/Pic2najHfouwW4NRbV70maJvoj+XS7gQQfLnAcp5bcDSyUtEDSDOAsYE1TmzXAOXGWwVLgeTPb2qXvV4HTACS9FphBl28hQTlZx3EcADOblHQhcDswClxjZhslnR9fXwWsBc4ANgM7gPM69Y1vfQ1wjaSHgD3AuZ1CBQDqct1pQtLKAcV+CiWE5wjhGcCfI3RcZB3HcQrEY7KO4zgF4iKboNdleJLmS/qWpIfjpXYXDX70+8bY81LC+PqopPsk3Ta4UR9Mn0siZ0u6JV76uEnSWwY7+gPG2c9zXBL/PD0k6UZJQ1nJkOIZXi/p3yTtlvSxLH1rgZn5EYVMRoF/B15FNGP4ALCoqc0ZwNcAEeXJ3RWfnwssjl8fDny/uW/ZnyFx/X8AfwfcVsV/i/jatcDvxK9nALOr9hxEye8/AGbF728G3l/SZziWKGf0z4CPZelbh8Od7H56XoZnZlvN7F4AM3sR2MTBq0sGQT9LCZE0D3gXcNUgB92Cnp9D0pHAW4GrAcxsj5k9N8jBJ+jr34Mo+2eWpDHgEODHgxp4gq7PYGbbzOxuoLnyfZrnDx4X2f30swxvH5KOB04C7sp9hN3p9xm+AHwcyFagIX/6eY4FRInhX4rDHldJOrTIwXag5+cws6eIio/8CNhKlMP5jQLH2o40z1BE32Bwkc0RSYcBXwYuNrMXhj2eLEg6E9hmZvcMeyx9MgYsBq40s5OAl4HKxQIlHUXk+hYArwQOlfTe4Y7K6QUX2f30swwPSeNEAnuDmd1a4Dg70c8znAq8W9IPib7WvUPS9cUNtSP9PMcWYIuZNb5J3EIkusOgn+f4ZeAHZva0me0FbgV+scCxtiPNMxTRNxhcZPfT8zI8SSKKAW4ys88PdtgH0PMzmNllZjbPzI6P+33TzIblnPp5jp8AT0pqFCs5HXh4YCM/kH6Wdv4IWCrpkPjn63SiWP+gSfMMRfQNh2HPvJXpIJrp/T7RjOgfxufOB86PX4uomO+/Aw8CS+LzvwQYsAG4Pz7OqNIzNN3j7Qwxu6Df5wBOBNbH/x5fBY6q6HN8CngEeAj4W2BmSZ/h54m+QbwAPBe/PqJd37odvuLLcRynQDxc4DiOUyAuso7jOAXiIus4jlMgLrKO4zgF4iLrOI5TIC6yjuM4BeIi6/REXE7wQ0P67JMkXd3h+jGSvj7IMTlOO1xknV6ZDQxFZIFPAJe3u2hmTwNbJZ06uCE5TmtcZJ1e+TTwakn3S/qspMMk3SnpXkkPSloBUVUyRZvOEb//mKQ/aXdTSR9RVPx8g6SbWlw/HHijmT0Qv39bPIb746pbh8dNvwr8Vo7P6zg94bvVOr1yKfAGMzsRIK55+mtm9oKkOcA6Sb2sU78UWGBmuyXNbnF9CdEy0wYfAy4ws+/GVdB2xefXA3/aw+c7Tq64k3XyQsD/lrQB+GeiuqGv6OE+G4Ab4rJ+ky2uzyWqF9vgu8DnJX2EaAeERp9tRCUCHWeouMg6efFbwDHAm2J3+1Nggkgokz9n3fapehdRwZTFwN2xQ06yM3kPM/s08DvALOC7kl6f+JydvT2K4+SHi6zTKy8S7WfW4Eiiot97JZ0G/EJ8/qfAsZJ+TtJM4Mx2N5Q0Asw3s28BfxDf87CmZpuA1yT6vNrMHjSzPycqrdcQ2ddyYFjBcYaCx2SdnjCzZyR9N57U+hrw58A/SnqQKB76SNxur6T/BXyPqGDzIx1uOwpcH+/TJeBya9qfy8wekXSkpMMt2k/t4ljUp4GN8VgATgP+Ka/ndZxe8VKHTuWQdAnwopm13fBR0r8CK8zsZ4MbmeMcjIcLnCpyJbC73UVJxwCfd4F1yoA7WWcoSLqCaF+xJF80sy8NYzyOUxQuso7jOAXi4QLHcZwCcZF1HMcpEBdZx3GcAnGRdRzHKRAXWcdxnAL5/yNKZvrrTYUSAAAAAElFTkSuQmCC",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f65872ce250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Plot it\n",
    "fig, ax = plt.subplots(1, 1, figsize = (5,5))\n",
    "a = ax.contourf(ce_theory, 40, extent = [0.01, .1, 10*dt, 100*dt])\n",
    "plt.colorbar(a)\n",
    "ax.set_xlabel('tau_s (s)')\n",
    "ax.set_ylabel('Delta_s (s)')\n",
    "ax.set_title('')\n",
    "plt.savefig('./fig_s1_e.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}