doc/userdoc/model_details/IAF_neurons_singularity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IAF neurons singularity\n",
    "\n",
    "This notebook describes how NEST handles the singularities appearing in the ODE's of integrate-and-fire model neurons with alpha- or exponentially-shaped current, when the membrane and the synaptic time-constants are identical.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing(use_latex=True)\n",
    "from sympy.matrices import zeros\n",
    "tau_m, tau_s, C, h = sp.symbols('tau_m, tau_s, C, h')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For alpha-shaped currents we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = sp.Matrix([[-1/tau_s,0,0],[1,-1/tau_s,0],[0,1/C,-1/tau_m]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-singular case ($\\tau_m\\neq \\tau_s$) \n",
    "The propagator is: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}e^{- \\frac{h}{\\tau_{s}}} & 0 & 0\\\\h e^{- \\frac{h}{\\tau_{s}}} & e^{- \\frac{h}{\\tau_{s}}} & 0\\\\\\frac{\\tau_{m} \\tau_{s} \\left(- h \\left(\\tau_{m} - \\tau_{s}\\right) e^{\\frac{h \\left(\\tau_{m} + \\tau_{s}\\right)}{\\tau_{m} \\tau_{s}}} + \\tau_{m} \\tau_{s} e^{\\frac{2 h}{\\tau_{s}}} - \\tau_{m} \\tau_{s} e^{\\frac{h \\left(\\tau_{m} + \\tau_{s}\\right)}{\\tau_{m} \\tau_{s}}}\\right) e^{- \\frac{h \\left(2 \\tau_{m} + \\tau_{s}\\right)}{\\tau_{m} \\tau_{s}}}}{C \\left(\\tau_{m} - \\tau_{s}\\right)^{2}} & \\frac{\\tau_{m} \\tau_{s} \\left(- e^{\\frac{h}{\\tau_{m}}} + e^{\\frac{h}{\\tau_{s}}}\\right) e^{- \\frac{h \\left(\\tau_{m} + \\tau_{s}\\right)}{\\tau_{m} \\tau_{s}}}}{C \\left(\\tau_{m} - \\tau_{s}\\right)} & e^{- \\frac{h}{\\tau_{m}}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡                                         -h                                  \n",
       "⎢                                         ───                                 \n",
       "⎢                                          τₛ                                 \n",
       "⎢                                        ℯ                                    \n",
       "⎢                                                                             \n",
       "⎢                                          -h                                 \n",
       "⎢                                          ───                                \n",
       "⎢                                           τₛ                                \n",
       "⎢                                       h⋅ℯ                                   \n",
       "⎢                                                                             \n",
       "⎢      ⎛               h⋅(τₘ + τₛ)          2⋅h          h⋅(τₘ + τₛ)⎞  -h⋅(2⋅τ\n",
       "⎢      ⎜               ───────────          ───          ───────────⎟  ───────\n",
       "⎢      ⎜                  τₘ⋅τₛ              τₛ             τₘ⋅τₛ   ⎟       τₘ\n",
       "⎢τₘ⋅τₛ⋅⎝- h⋅(τₘ - τₛ)⋅ℯ            + τₘ⋅τₛ⋅ℯ    - τₘ⋅τₛ⋅ℯ           ⎠⋅ℯ       \n",
       "⎢─────────────────────────────────────────────────────────────────────────────\n",
       "⎢                                                2                            \n",
       "⎣                                     C⋅(τₘ - τₛ)                             \n",
       "\n",
       "                                                  ⎤\n",
       "                                                  ⎥\n",
       "                                                  ⎥\n",
       "                          0                    0  ⎥\n",
       "                                                  ⎥\n",
       "                          -h                      ⎥\n",
       "                          ───                     ⎥\n",
       "                           τₛ                     ⎥\n",
       "                         ℯ                     0  ⎥\n",
       "                                                  ⎥\n",
       "ₘ + τₛ)         ⎛   h     h ⎞  -h⋅(τₘ + τₛ)       ⎥\n",
       "────────        ⎜   ──    ──⎟  ─────────────   -h ⎥\n",
       "⋅τₛ             ⎜   τₘ    τₛ⎟      τₘ⋅τₛ       ───⎥\n",
       "          τₘ⋅τₛ⋅⎝- ℯ   + ℯ  ⎠⋅ℯ                 τₘ⎥\n",
       "────────  ──────────────────────────────────  ℯ   ⎥\n",
       "                     C⋅(τₘ - τₛ)                  ⎥\n",
       "                                                  ⎦"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PA = sp.simplify(sp.exp(A*h))\n",
    "PA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the entry in the third line and the second column $A_{32}$ would also appear in the propagator matrix in case of an exponentially shaped current"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Singular case ($\\tau_m = \\tau_s$) \n",
    "We have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{1}{\\tau_{m}} & 0 & 0\\\\1 & - \\frac{1}{\\tau_{m}} & 0\\\\0 & \\frac{1}{C} & - \\frac{1}{\\tau_{m}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡-1           ⎤\n",
       "⎢───   0    0 ⎥\n",
       "⎢ τₘ          ⎥\n",
       "⎢             ⎥\n",
       "⎢     -1      ⎥\n",
       "⎢ 1   ───   0 ⎥\n",
       "⎢      τₘ     ⎥\n",
       "⎢             ⎥\n",
       "⎢      1   -1 ⎥\n",
       "⎢ 0    ─   ───⎥\n",
       "⎣      C    τₘ⎦"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "As = sp.Matrix([[-1/tau_m,0,0],[1,-1/tau_m,0],[0,1/C,-1/tau_m]])\n",
    "As"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The propagator is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}e^{- \\frac{h}{\\tau_{m}}} & 0 & 0\\\\h e^{- \\frac{h}{\\tau_{m}}} & e^{- \\frac{h}{\\tau_{m}}} & 0\\\\\\frac{h^{2} e^{- \\frac{h}{\\tau_{m}}}}{2 C} & \\frac{h e^{- \\frac{h}{\\tau_{m}}}}{C} & e^{- \\frac{h}{\\tau_{m}}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡  -h                 ⎤\n",
       "⎢  ───                ⎥\n",
       "⎢   τₘ                ⎥\n",
       "⎢ ℯ         0      0  ⎥\n",
       "⎢                     ⎥\n",
       "⎢   -h      -h        ⎥\n",
       "⎢   ───     ───       ⎥\n",
       "⎢    τₘ      τₘ       ⎥\n",
       "⎢h⋅ℯ       ℯ       0  ⎥\n",
       "⎢                     ⎥\n",
       "⎢    -h      -h       ⎥\n",
       "⎢    ───     ───   -h ⎥\n",
       "⎢ 2   τₘ      τₘ   ───⎥\n",
       "⎢h ⋅ℯ     h⋅ℯ       τₘ⎥\n",
       "⎢───────  ──────  ℯ   ⎥\n",
       "⎣  2⋅C      C         ⎦"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PAs = sp.simplify(sp.exp(As*h))\n",
    "PAs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numeric stability of propagator elements\n",
    "For the lines $\\tau_s\\rightarrow\\tau_m$ the entry $PA_{32}$ becomes numerically unstable, since denominator and enumerator go to zero.\n",
    "\n",
    "**1.** We show that $PAs_{32}$ is the limit of $PA_{32}(\\tau_s)$ for $\\tau_s\\rightarrow\\tau_m$.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAF0AAAAaCAYAAADVLFAXAAAABHNCSVQICAgIfAhkiAAABMNJREFUaIHt2VuIVlUUB/CfmjKa0kOWSZYllVmmjlJa2J2kLBnH7j2oFD0kJBJdqB6yh+iCkUpFSERmVFTeMruYRUhlNxXtImVm9GKpPXTVSrOHtQ/f6cz5ZubT0Wlg/nBY395n7b3XXmfd9v66zJw5UycOLrq2twCtwIN4vb2FaEt0BKWPwPpW8D2JRw6wLG2CjqD04VpWeldchqUHXpz9R5nSb8FeXHeQZSnDUeiHv/EGfsdmXFDgG43ueC+1l4g9lD2TD7jULaBM6SMTXXswBamC+kSnY5aw+i/wcIGvAcuxO7VvQH8ch39wZWr3x/P7IU8/7MHckncn4n6swXZhKNuxEtPQK2M8pGTwSGFRX++HcG2FEfgFV2Nr6luI+wp8Dbg71/4p0XphWKuwrQ3kaUjzLc71dcG9uAM98AFexs8YiHG4ENfibJoq/VAMxmphIW2NmbinBZ7z8W76PQKvqCgcBuGbXPskYdFvlsw1DD9oG4VDo/igq3J9T2EqvhQhuZh/eoqQfULWUQwvI1LfWpyCZ5PQv4kvOLqKMJeLsm4H/sIm3IVuBb5HMaSF5+OCPKsLc9RjXa7dgLeFdxYxTHkSXobZ+DDJOlpY73e4tcoeDxO5ZJkIMYQyp2IjxlRZa6fwzGlZR9HSRyU6EJ+mzTyDUzFeWN0J+DXxdRMf5hphfS/hT1ySFhqMKbn5d6SnNeiV1lpX6K8X7puhAU9XmWMQvirpHyqUNwPzRI4Yj74iN8wqGXOpCB+LUru/2ONuXKWik2rYmf0oWnqWRMdgLCbg9rTgIhwprC/DHKHwB4SV3pQ2MlR4xmThMfuC4YluyPUdjgEqH+IIYaXLqszRVRjQAJW99km/56X2LjwmcscuEYvL0Ci86a3UnoE6LMDnrdlQXqg8MqVP1bR62ZhoXaKjhcssxZ0qlQORuefn+PYFw4Xr58NGvfCkTJYJ+AQ/VpljThqzRZSUhEHkQ1i+PVS5AutwsQihu1LfxEQXtLyV/yIfXuqEtX6r/Ng9KNHNid4sMvcfIkEWMTTRLrUKlfBEevJYqfLRidCypJk53hGWXpQr7z3Hiz0TOeCzknnGobdK1dJbJPC9+KiZ9UuRV/rw1F5RhXekcL0tOUGIUqg5fF+rUDXgfbXX3aeJXEXE5R+E8rJ380vGNIoCYXlqH5HoL8LoakJe6VkSXVPC10cU/6uSgHVp4VU4t9ZF2xAP7cOY6bnfW3Fmrj1FU3QTYewdlXifhbye6f2eknFVkY/pWTwvU3p2yMjeZSGjby2LdVCcIxJ4/kC0TZSXPXBeC+ObnPqLSv9LeSLJvCBLrjtFXDwFk6osNlbTOr0jYpI4KBYv02Yn+jhOLhnXRZTOLxRfZOGlh0gw60XlUURZ6LlNxLiFIsFtEB/x6MTfHcc2t5sOgomi/C1WSHNF4r1eJN8V4urkbxwjjG4AnitOmCn9NKGkapdco8SpdFOubwXOEncOY4Wb/Szi5Eq8WMPG/q84XSiu7J5+r7hYW4gbcQYuEnraKk7Sy/FqcWCm9DWaL+2GVOn/BFe0LHuHRWOii5vheS09rUZH+BOjPdEoQu6WlhhrQdnVbicqqObh+4VOS28HdCq9HdCp9HZAp9LbAf8CUHr5FKFsEGIAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\frac{h e^{- \\frac{h}{\\tau_{m}}}}{C}$"
      ],
      "text/plain": [
       "   -h \n",
       "   ───\n",
       "    τₘ\n",
       "h⋅ℯ   \n",
       "──────\n",
       "  C   "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PA_32 = PA.row(2).col(1)[0]\n",
    "sp.limit(PA_32, tau_s, tau_m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2.** The Taylor-series up to the second order of the function $PA_{32}(\\tau_s)$ is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAArCAYAAACHItq+AAAABHNCSVQICAgIfAhkiAAAEfxJREFUeJztnX+0VNV1xz8PjYpgQNEoDZEfBkF8j/dALAQN0ahgalpEm/gjEZ+xXam2EpYrJBrTtZ62haSGClapK79KEtPSNohKrJEKGowpFVE0SWlDlRdjwZ+gGBXyXrV/fM/h3rnv3pk7b+bemXmzP2vNmrnn3rmz79xzzz57n332aenq6sJoSj4AfA94H9AD3AjcVVOJmo+PA3OAa2L2rQK6gP/KUyDDMBqHg2stgFEzeoHPAU8hJb4F+BHwVi2FajImo/8/yiBgPKa8DcMowqBaC2DUjF0EyuMlYA9wdO3EaUragROAnwLPA1Nd+SmoQ3UksBUp8r3u8935i2kYRj1iCnzg8lXg/pTHTgPeA/w6VHYk8CJSMM3MD4BrMzr3ZGAHMBNYBsxz5XOAdahT1QEscHJ0AOdnJIthGA2GudAHLh3AkymOGwF8F7gSeDdU/iXkUn+m+qI1FDcBDwPfAl4v87s/TyifDexGnaavu7Ie4FX3+Szg9tDxrcB/lvnbhmEMcMwCH7i0Ez++GuZQYA2wBLlxPYcDfwx8MxvRqspKFOyVFU8jK/nT/fhua8Jrp3t/PHRsG1L4R6Dnck9o3yRMgRuGEaERFHg5rmBDHAcci6y6HwFvIkv6o6FjWpDy24Ci0cOcB7wD/KRCOe5GVn3ca36F585TlnuBSyNlo4DvI0v6NWA1CgZMSzvws9D2ZLd9FvBQ5NhjUcyCZwzwILAI+HvgT5yMI8r4fcMwGpxGUOAdlLYkQdbiLRnL0ihMce8LgK8hZfELYGnomNOAi9CY6lb3anP7TkdBVGGXen+4EhiJFM47wCfc9kjgHys8d56y/AdwKjDYbY9D/89z6H88AwUArihDHq+wQc/h+1DMwWzggcix96MOwsluux34J+Bm4BjgDuRBKacDYRiVcj2wGQVYvgysRZ4lIycaQYGncQUPQnNq78lenIagAz1UFyFL7X+QAjgmdMxP0P/WEXp5hTIGuXkr5VXgBeAo91sb3fYLyDuQJ5XIshONV/+O274DjYlfD2xDnZ/FyHpOy+eQ1QzqUIxxn38XeCxy7ArUafiF225H93UYwVSzD6L73OicBwyvtRBGKs5AdXMm8u71onp5VA1lGqh8ErVBBcQp8GuR5RV1GdaCNK5ggOno4rzLt55ct7WgAymHsNt1HOkb+MOAfTHlXST/r/51Rsz3JiNF+VLK3y/Gl4DfhF6fiin7cJHv90eWt937YGA0cA6wMPKbdwO/LeOcSUwD/q/EMWOBbgqt+IPJv1NUbU4FjkdDEnF8F923IblJ1Nycgp7pKxP2z0FDOD9H9fAyZCSclot0zcXDwF9EC+MUuJ+L+kSW0qQkjSsYYC5wH+oBQnau22NR43przL7xKBhsC3In9bj3B4GrUWBYXnQA/x4pm0K6qHSAV9A0sii3ASeVeEWtR0hOWLIWTZ/aBGxHHbE1SDl9PkG2Oyj0GtwbU/Z4wneLyQIK3HsSNUarQ+XeongZ1cHX3XnCv9mGlG8eXIEa1keA77iyzpx+Ow3jgL9C/+VuYD+6pyvR/xdHC/Blgqj8KNNQIOFXUEfeyJ4tqGP6l8DQFMf7AMzdWQrVpPhcHR8JF8ZNI5uKHpBf5iBUKcKuYG9NrkaNQ5i5wA2hbT8dZwqBu7Qa1t9cd741obIWlIb0i8AhaCzyB6iRH43GNM8CLqG4ZVgtDkfu1KiynuLkSsOTxCuEV9yrXMYB/x1T3oqU+ELUcC8Ffg+NJ9+HOm1RdlPYQLzhttN6F5JkORL4M2R19FLoxm1DbvQXUcdsqPv8RsrfbBa8Er4BPQs/Rl6xN9GzPB95TD4LfDvy3YtRxH+S92Exagv+rupSG8VYgmJAFqB7UIxlaDhpU9ZCNSnfQ7EvB/RI1AIfAkxAN+GdDAToojwXbBpX8InI0o4G/kB1XbegRBuvog6B59vAnzuZOpD76CrgOqS0R6FGrdzxyU6SXdLF8BbO06GyEU6OtBb4A8iarlZmtkGoMzOKoM753rq3uPahuc973edy51xXIgsESvurSGGH3biz0BAOqHHaA9yJOkUnIJf67UiBNSqd9K++eVrQs3AT8mBMBM5Eed6vA85F/xPonk+JfP861DjFcSJwNvDPBMMZA5lOKrsX1eQxFGfxWeCgIsfdjJ6TT1B6CMjoHzvRfzvLF0QVeIcrewLNPb0TKcDfIMtyesKJL0SRsq+gccDtaFwyesPLdcGmcQXPBdYT71arput2GBp7X0tQQa9FD9s2YEbCb72NPAZXJ5y32rSjawr/H1OQG3NbynP8DPW6L66STMudDDsIAjFaKbzX4e1WkpOgZCELyJpuQ/dwFarToHHv84FvuO09wMdQB+Qh1Nm9GWWxqzRqv5G5Dj0LW5CFEOfBW48s6IOQ18UzEXmNkurnZ1AHIUnBV8rRyGApZVzsR/Eh9UpW17EKxSacnbB/KfKunEWyoTIS3edhNZC/XMYRP4QYR94yb0EBbUBfF/op7n00GkdcjwJHTkauzXvRg+ZdhwchJX8xunH/4gT9GFJaE4DLQ+cvxwWb1hU8F42txVFN1+15yC3oV+waia6xF/2hpdypeVkOd7hXmAcpv/LcCPwtanAr7VFvQHUqTCuFXoKxwLPuczg4qxSdVZAFFMOwHdX3aSjJDSieYhOFbsHH6RtI2cyMRfVlH7LA4gIgPQ8gq/xDobKPoPudVM/Odvuycs0ORZ4DzxjUbm0Bfhgqf5ni11ZrsrqOR937OfT1dN6KhjjPpHiyoSVOlitIbq/r4T4MR0M/L6J6lxRQ6clb5idRZxnoq8B9ANsMNBc4HMi2GrgAWcWPuLLlSHl/BbmRfRDZIhQ1Nx+5JPuTRSqNK/gYZD1fSDxhd+lO1FMq5rodSrLrdh6yav/NbS9EStFHYQ40HkAek1HArzI4fxvqIII6Qy8QWLBtBMFZefFlVO/fRF4fb+39lvjlPo2AzyNvxgrk2SiGz7cftsTa0GIucQxBbc42sgte66Ywm98foUZ4FfGd+Xqlm2yuY7N7nxUpX4ECC89HcSjHuXI/M6Ncuqn9fXgNXdditB7BORQfzusmX5n/FxnGBwO9SQq8k75R6N695S256cgtfA+aDxumBzXAM91x/VXgpVzBv48q14sJ51iOFOwOpJz3E++6XR76HKeMD0NjePcT9KL8ohLRLGYDibho+2qxIPR5F4UW2eXkT9JvJkVFG6KFYBGWO1Mc77PFhT1xHyAIPI3yfuTp25WwPws63PvWHH8zC6p1Ha+jdu/4SPlV7n19pPxG+qY37qR8b1mt7sMS974YGTKzkYGXhqxl3ouU93HA82EFfhgag36W+NSl49y7X9ziGvTwvkV8Lmqfkae/gT1pXMFzKb68YrVct7NRB8BHnw9FgTXvorHiSumOkdMTTasJ6hx1VuF3jeakm+rVNz89s5d0jdYM9x42EIaQbLl7hb8nYT8Uv544vk/x3Pbe+5cmA2SldJPds1/N69iNptGGyTpos9r3oZvy6gnIAP0WGhpKQ9Z1x3ckhkKhC73dba9L+OJU1BPzD9ps935JiR98rnwZU/Mo5c/r7o/rdh5ypd7ntn1Gs72oA1Mpy+ibfaoDdVC+gypemEa3DIzaUs365mcqvEEwhJZEC5pGBkEsCWh8u0+WKYePHSkWw/EM5Y0vFssy2II68jvRuGXWZPXsV/s6BpPvDIAs7kM59WQQiouB9HFbedQd/5z0QKEC9wFsW2K+dAS6mI1I2R2GlNhGIhPLc+av+/Gdcl23ByFX/QaCsRDv1h/s9lca5LUspqyTIEDv4QrPbxhhqlnffJDPcBR4WqxDeykKiN1GYdrjPcB7E77jp4AWW6ilnBS2pRjrZHm01IFVIqtnv5rXMQjd31LxDdUki/uQtp4MQmtrjEeB2WljYPKoO/452QOF08j8+HecAvcJUfw+7zqp1jzhemYWajzCyVteQj3jQyg9VzPvfPOlpjPYq3leefAcCqxpIZjnHceJKDioFwX6hPNMPEty/vNdyJqZULGk6Zjo3osFpo5C/+8nUcf+LeQynYBcro+6sk3u2FqQ5jpagC+gmTr7ULsWl+xpgjs2T89fNeUvh0Eon8EV7lyXUtqz5Mmj7gxHhuRuL6xnKnITx/24t879uNXbaBx5EopMj+N0ik/8bxQuQI1NdKEU33NeQXDjwrSg6XSrshMtlhZ7pXodj6aLbEN1+cI6kKnar7zwqwD+DcGCL2E+jhqmIShd7U8j+zcT/wyBGruNyFj4YMWSlsZbOMWClnyg0tUoYGs6wYyUxWgWzmnov1gYd4IcSHMdi5CiuhopkD8gmGUTxsctPFSBPCvRvexMeXw15S+HIShO6i40PJxWeUM+dWciobgr70I/xAn9FPELIsS51xehMeHVKLjsadQheL87/j30jVpsRM5HDU400v1WNN7xGRT4tg4lr+hBUbWnox7UP+QmqVEOvWhFsKfQMpxbULa1asQ0NBu3oIbp0yhr1z3IMvcLW0xCsSbzUP6FKBvQs/Je4hu/1aiDNYfsV1zzCWgWohz4m+kbZ9NOkOLZtwvr3PZEglS/P0ZxNrUgzXWcC/wrQUzQr4ifaz8bDRNWstqjNxbTKsRqyl8ObyBX+1uUp7whn7ozhSAr5IE/tQ0p3KQFTE5B8/q2h8rWoWliq933F6DVaE5CCv2KhHM1EqeihmVNzL53UZKP89CE/Q7gT5F7cDKaS9xJ4Zi7UT/sIogU9QsFNMOQUBa8g579P0S5z89Frs25aHrYtchCilPeIJfgXfSdZ+xZjRq7PFYSfALlBOhBDXHcAjUd6FrCnfrR6BrCefqPJ99x4zBprmON2/cgarvi6v8wZMT8kGAOf39oQ8rxvlIHOqolf394nf6t7Jd13Wlx5zzQKWjp6urqh5xNw2I0x30ctXsQjeyZhubzTyK/sWOjkPHIpZi0jPH16HmcSvqc/lnxS+SBuy1U9muUzGplqOw11MkPr2xXb5yAOlrzURDWDApT2l6DrnUWQQKvchmOOnJLUceumpSSv96opO7MQsMEB9J95x1g1WjMQ1aaKe98uB65nfaiwKW1BPkEorSj+bzPowCWHUgJn1zmb45A6VOvxJR3LdmOvCKjE/bfgtzyNyXsz4shSGmEOxFH03exoLHIeq11Z6MUz6DYhWnIwpsc2jcYPZOr6b/yBuXG73G/U22KyV9vVFp3PoWynh7AFHhxTiIIOjCy5wwUFDgT5RrvRe6xoyLHzUdj1vvR2NGJBO7VcoYsDkVuuCX0Dawy8ucGkofe9iE3/eOoIawVXkGEE3VMRXUxnHGyg8K8GfXGF9EQ3yT0/NyIgpgfDh0zBmUiTFrgKS1rUaDWCxWeJ0wa+euNSurO2SgqvmBOetx64IZRK+ZEti9DFfk0gvHTmWiaxyKC6GeQdfYIgbK/G7nW4rgcWesrUQDVQE6H20jsQ9kXhxGff3ojhUv51gKf4jmc63sKmr3TEzluK/Xr1TkULYoxGuW12ISCt8Jjs9uIz7JZD6SRv96opO48jwJEC7AxcKOeGYmyGp1OkBzhMdTQJwU8eUagwMxD0Tzji1CAFWg8bjpSBuG0upeRfhU0wzCMmmIWuFHPLEM9UT81ZCKaGZBmnXK/OIZPQrSRIKsXSJnbEJJhGA2LKXCjXrkZWdkfJkhV6+MR4rIFJjEZjb29VOpAwzCMRsIUuFGPLEVJQc6kMHHHYPdezlrDk8lnVSnDMIxcMReiUW/cSqC8o+vI+zS/SePfh8eUjUO5kg3DMAYUpsCNemIFmhpyCcpKdJx7DXX7N6PUibe548aj/NgXoXSK7fRlEIpUHUVhfV+Lxtg3ocjQ6WhKWTeVT5sxDMPIHFPgRj1xFVq6dj1K6uFfYYV6AfA1lJ5zK1LqX0BR6nFJM5ajQLYdFK453Yos/BlokYalaHrZR1FSF8MwjLrGxsCNeiLNClr70TrwadeC30Df7F5HoM7r1932PuB2lAFuKPFzkA3DMOoKs8CNZqQVzSeP226l+Hq+hmEYdYEpcKMZaaUwgctYlOwFFLVuyVwMw6h7TIEbzUgbgQIfieaJvxvaZwrcMIy6x8bAjWYkvODJLuBDoe3Lc5bFMAyjX5gFbhiGYRgNiClwwzAMw2hATIEbhmEYRgPy/07DM35F11SiAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\frac{h e^{- \\frac{h}{\\tau_{m}}}}{C} + \\frac{h^{2} \\left(- \\tau_{m} + \\tau_{s}\\right) e^{- \\frac{h}{\\tau_{m}}}}{2 C \\tau_{m}^{2}} + O\\left(\\left(- \\tau_{m} + \\tau_{s}\\right)^{2}; \\tau_{s}\\rightarrow \\tau_{m}\\right)$"
      ],
      "text/plain": [
       "   -h                   -h                           \n",
       "   ───                  ───                          \n",
       "    τₘ    2              τₘ                          \n",
       "h⋅ℯ      h ⋅(-τₘ + τₛ)⋅ℯ       ⎛          2         ⎞\n",
       "────── + ────────────────── + O⎝(-τₘ + τₛ) ; τₛ → τₘ⎠\n",
       "  C                 2                                \n",
       "              2⋅C⋅τₘ                                 "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PA_32_series = PA_32.series(x=tau_s,x0=tau_m,n=2)\n",
    "PA_32_series "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore we have \n",
    " \n",
    "$T(PA_{32}(\\tau_s,\\tau_m))=PAs_{32}+PA_{32}^{lin}+O(2)$ where $PA_{32}^{lin}=h^2(-\\tau_m + \\tau_s)*exp(-h/\\tau_m)/(2C\\tau_m^2)$\n",
    " \n",
    "**3.** We define\n",
    "\n",
    "$dev:=|PA_{32}-PAs_{32}|$\n",
    " \n",
    "We also define $PA_{32}^{real}$ which is the correct value of P32 without misscalculation (instability).\n",
    " \n",
    "In the following we assume $0<|\\tau_s-\\tau_m|<0.1$. We consider two different cases\n",
    " \n",
    "**a)** When $dev \\geq 2|PA_{32}^{lin}|$ we do not trust the numeric evaluation of $PA_{32}$, since it strongly deviates from the first order correction. In this case the error we make is\n",
    " \n",
    "$|PAs_{32}-PA_{32}^{real}|\\approx |P_{32}^{lin}|$\n",
    "        \n",
    "**b)** When $dev \\le |2PA_{32}^{lin}|$ we trust the numeric evaluation of $PA_{32}$. In this case the maximal error occurs when $dev\\approx 2 PA_{32}^{lin}$ due to numeric instabilities. The order of the error is again\n",
    "\n",
    "$|PAs_{32}-PA_{32}^{real}|\\approx |P_{32}^{lin}|$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The entry $A_{31}$ is numerically unstable, too and we treat it analogously."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tests and examples\n",
    "We will now show that the stability criterion explained above leads to a reasonable behavior for $\\tau_s\\rightarrow\\tau_m$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nest\n",
    "import numpy as np\n",
    "import pylab as pl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Neuron, simulation and plotting parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "taum = 10.\n",
    "C_m = 250.\n",
    "# array of distances between tau_m and tau_ex\n",
    "epsilon_array = np.hstack(([0.],10.**(np.arange(-6.,1.,1.))))[::-1]\n",
    "dt = 0.1\n",
    "fig = pl.figure(1)\n",
    "NUM_COLORS = len(epsilon_array)\n",
    "cmap = pl.get_cmap('gist_ncar')\n",
    "maxVs = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loop through epsilon array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'voltage V (mV)')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeVhV1dfA8e9iFlAUQQVRQMEJRVQ0fM0xLae0QS0z08zqZzbZZGVZWTZaWTll5ZCVpuZUOZSz5QA4I4aoOKCIguKAyrjfP+4VUQFRuPcw7M/z3AfY555z1rV0sffZe21RSqFpmqZpt8rG6AA0TdO00kknEE3TNO226ASiaZqm3RadQDRN07TbohOIpmmadlvsjA7Amjw8PJSfn5/RYWiappUqW7duTVJKeV7fXq4SiJ+fH5GRkUaHoWmaVqqIyOG82vUQlqZpmnZbdALRNE3TbotOIJqmadptKVfPQDRNK18yMjKIj4/n8uXLRodSKjg5OeHj44O9vX2h3q8TiKZpZVZ8fDwVK1bEz88PETE6nBJNKUVycjLx8fH4+/sX6hw9hKVpWpl1+fJlqlatqpNHIYgIVatWvaXemk4gmqaVaTp5FN6t/lkZlkBE5DkRiRGRPSLyqbltgIjsyPXKFpGQPM51F5G/RSTW/LWK9T+BVpJlZytWrNjPhx9uYPr07aSk6DFwTStuhiQQEekI9AaClVJBwDgApdTPSqkQpVQIMBA4pJTakcclXgdWKaUCgVXmnzUNgOTki3TqNJOuXX9m1KjVDBmyhDp1vmL+/GijQ9PKqSFDhlCtWjUaN26c53GlFM8//zwBAQEEBwezbds2K0d4e4zqgQwDPlZKpQEopU7m8Z7+wOx8zu8NzDR/PxO4r9gj1Eqlixcz6Nx5Fps3xzN1ak8uXHiD8PCh1K/vQd++85g2bbvRIWrl0ODBg1m+fHm+x5ctW0ZsbCyxsbFMnTqVYcOGWTG622dUAqkHtBWRLSKyTkRa5vGeh8g/gVRXSiUAmL9Ws1CcWinz8ssr2LnzBAsWPMSTT7bAxcWBli1rsnbtIO6+uy5PPfU7a9ceMjpMrZxp164d7u7u+R5fvHgxjz32GCJCWFgYKSkpJCQkWDHC22OxabwishKokcehUeb7VgHCgJbAXBGpo8z764rIHcBFpVRUMcTxFPAUQO3atYt6Oa0E27r1OFOmbOWll8Lo3j3wmmOOjnbMm9eXli2/47HHFrJ79zDc3JwMilQzwosvLmfHjhPFes2QkBqMH9+1yNc5duwYtWrVyvnZx8eHY8eO4eXlVeRrW5LFeiBKqc5KqcZ5vBYD8cACZRIOZAMeuU5/mPx7HwCJIuIFYP6a1xDYlTimKqVClVKhnp43FJPUypCRI1fi4eHM6NEdOJ4MZ85fe7xSJUdmzryPY8fO8/rrK40JUtPyYP7d+RqlYfaYUQsJFwGdgLUiUg9wAJIARMQG6Au0K+D8JcAg4GPz18UWjVYr8SIjj7NqVRzDRvXks1e6cmnnfpyquXLwjq8ZO7QLdcy/yIWF+TB8eEsmTozg2WdbERSkRz/Li+LoKViKj48PR48ezfk5Pj4eb29vAyMqHKOegUwD6ohIFDAHGKSupuB2QLxS6mDuE0TkexEJNf/4MdBFRGKBLuaftXJswoRwGoTZEDepM2O/X8sXEfF8+Od/hH08nH9+HMiOA1ffO3p0eypWdODVV/82LmBNy6VXr178+OOPKKXYvHkzbm5uJX74CgzqgSil0oFH8zm2FtOzkevbh+b6Phm4y1LxaaXL6dOX+HXJPu7iS5afOcu7d9fhmemL2TR/LN2/jMf27UR+VgPwevpnqlcBDw9n3nyzLSNHriQ8/BitWtU0+iNoZVz//v1Zu3YtSUlJ+Pj48N5775GRkQHA//73P7p3787SpUsJCAjA2dmZ6dOnGxxx4UheY29lVWhoqNIbSpU9P/ywjTt2P8vprzKJuO8cLy/8L+fYqWMbOHDH8/wvYT+t33mfSW+/iAicP5+Gr+942rf3Y+HChwyMXrOkvXv30rBhQ6PDKFXy+jMTka1KqdDr36tLmWil3oaIvwj63pmwMFtGzL924p5nzbakft2T/dmX2Pf5u8xemw5AxYqOPPtsKxYt+o+9e08ZEbamlXo6gWil2smTqVSPnMDw1H3s+bAnNrY3jsre9cD7jOzqx+pzZ9n7bWsupZnan3uuFRUq2PHll5utHLWmlQ06gWil2o9z/+KHbSeJ97hAs46j8n3fq/O20NTBhR/n7WHiYtMTdU9PF/r3b8zPP+/m7FldK0vTbpVOIFqpdm7Z6ySrDHoMa13g+5xcq/L8S3dwJDuNtJ/vJdWcL4YNa8nFixnMmrXLCtFqWtmiE4hWamVlZbNhwwnq2VXgibcX3vT9g95fzjLfIF4P9+S7FYkAhIZ6ExrqzeTJkXku5tI0LX86gWil1oxZX7Hu/Dm6B1fBzt7hpu+3tbOn7pjG2J6oQOi2gWRnm9qHDQslOvoUGzceLfgCmqZdQycQrdSqu+kXnqU2LZ7qU+hz6vSfwddV9vPxF5EsjTD1OPr2bYSzs70extIsZvny5dSvX5+AgAA+/vjGdc/r16+nefPm2NnZMX/+fAMivD06gWilVptIe8bX8eXRp78q9Dm29k6caan488IZNs4bAJim9N5/fwN+/XUPaWmZlgpXK6eysrIYPnw4y5YtIzo6mtmzZxMdfe3eNLVr12bGjBk88sgjBkV5e3QC0UqlndH/Er4zk+NtLt3yuYMn/IALthz+9S+OmMtwDhwYTErKZf74Y18xR6qVd+Hh4QQEBFCnTh0cHBx4+OGHWbz42vJ9fn5+BAcHY2NTuv5JNqqYoqYVyfIvXub1rAg+rtSRkbd4rm9gB/r6u/NrXDL1fl/DO0905K676uDl5cqsWbt48MFGFolZM9aLU7imJlpxCKkL4/9X8HvyKtW+ZcuW4g3EIKUr3WmaWezaOJyx4e4nx9/W+fc8055LZJOx8BmUAjs7Gx55pAlLl8aSlHSxmKPVyrPSWqq9MHQPRCuV1h8+T7sqlWjWNPi2zu/74i9EjNlMz/32bIiCdk1gwIAmfP75JhYt+o+hQ5sXc8Sa0W7WU7CU0lqqvTB0D0QrddaunUls5iWCgyrd9jVs7ex59YlahMV4s2nzZMC0u1ydOlWYPz/6JmdrWuG1bNmS2NhY4uLiSE9PZ86cOfTq1cvosIqFTiBaqRM+cwIAbi2K9qwi44mn2ME5vBd9RXqGaVihb99GrFoVx+nTt/5wXtPyYmdnx4QJE7jnnnto2LAh/fr1IygoiNGjR7NkyRIAIiIi8PHxYd68eTz99NMEBQUZHHXh6HLuWqlzeEAocb/ChSVf0bN7myJdq31FN5LSMvn0nwv0aCVERh6nZcvvmDatF48/3qyYItaMosu53zpdzl0r02ruduTORo50aNeyyNfq3LYa0RkX2fTnCwC0aOGFn19l5s/fW+Rra1pZpxOIVqocOxHD6D0J/Ot3DlfXm5cvuZmeb72BAKl/LiAt3TSM1adPQ/7++wApKbpCr6YVRCcQrVT5c+obfJQdxxqX4plA2Oz/htDGpRIrd53hr+2m4dw+fRqRkZHN77/HFMs9NK2s0glEK1WObtgJgE2DHsV2zS53ehCXkcaOdW8A0KpVTWrVqsRvv+lhLE0riE4gWqlyKOYcnmLPHa0fLbZr9v/wHU7SkYFx68jINA1j9epVn7//PsilSxnFdh9NK2t0AtFKlaiTl2hayZmWob7Fds3A5o/hEHyRWivt2GDeUr1Xr/pcvJjB6tVxxXYfTStrdALRSo0zp48Qn5ZOg9oVcHevUKzX/iPoPKH7d7Jm5acAtG/vi6urA0uW6OcgWtHdrJx7WloaDz30EAEBAdxxxx0cOnQIgOTkZDp27IirqyvPPvuslaO+OZ1AtFIj5t9JJNKRHj0aFPu1Xe/vwQ7Ok7V8KkqBo6MdXbsG8Pvv+8jOLj9rpbTiV5hy7j/88ANVqlRh//79jBgxgpEjTSVCnZyceP/99xk3bpwRod+UYQlERJ4TkRgR2SMin5rbBojIjlyvbBEJyePcd0XkWK73dbf+J9CsrdL2CGwQ9roVff3H9Tr3/YhAuwpERiYRdcjUdu+99UhIuMC2bQnFfj+t/ChMOffFixczaNAgAPr06cOqVatQSuHi4sKdd96Jk5OTEaHflCHFFEWkI9AbCFZKpYlINQCl1M/Az+b3NAEWK6V25HOZL5VSJTMtaxYxf8E+ku2E+lVGW+T6Xeq58X10IgvXraeJfzu6dw/ExkZYsiSG0NCyUfyuPHsxazc71NlivWaIuDHetkmB7ylMOffc77Gzs8PNzY3k5GQ8PDyKNd7iZlQPZBjwsVIqDUApdTKP9/QHZls1Kq1EW3XgHOEOKTRvWtsi1w97uAXpKLL+egkADw9n2rSppZ+DaEVSmHLupbXku1Hl3OsBbUVkLHAZeEUpFXHdex7C1EvJz7Mi8hgQCbyslDqT15tE5CngKTBtG6mVXv+lXuKeWm40bGiZ38r6jJjB0jGNCYu7wPFk8K5qGsZ67bWVHD6cgq9vZYvcV7OOm/UULKUw5dyvvMfHx4fMzEzOnj2Lu7u7tUO9ZRbrgYjIShGJyuPVG1PiqgKEAa8CcyVXuhWRO4CLSqmofC4/GagLhAAJwOf5xaGUmqqUClVKhXp6ehbTp9OsLSLiN06qDLxru+DmZpnx4AquHnzfxZeue2uyaLNpqKNXr/oAeqtb7bYVppx7r169mDlzJgDz58+nU6dO5bsHopTqnN8xERkGLFCmflu4iGQDHsAp81sepoDhK6VUYq5rfQf8USxBayXWtj9mAGDvW92i9znR3QPHZalkbx8OPX6ifn0P6tWrypIl+xg+vJVF762VTbnLuWdlZTFkyJCccu6hoaH06tWLJ554goEDBxIQEIC7uztz5szJOd/Pz49z586Rnp7OokWL+Ouvv2jUqGRsu2zUENYioBOwVkTqAQ5AEoCI2AB9gXb5nSwiXkqpK1Nj7gfy66loZYTt/sM0wIU0rzCL3se++2v4P9eR/82rQfpIcLCH7t0DmDw5kosXM3B2trfo/bWyqXv37nTvfu1k0TFjxuR87+TkxLx58/I898qakJLIqIfo04A6IhIFzAEGqatPkdoB8Uqpg7lPEJHvReRKPfpPRWS3iOwCOgIjrBW4Zoz7VQWinVpTy/8hi96ndp32hLlU5N+Yc/xrnqrfrVsgaWlZrFmjV6VrWm6GJBClVLpS6lGlVGOlVHOl1Opcx9YqpW74NVMpNVQpFWn+fqBSqolSKlgp1StXb0QroyoetCHb/xIhTWpa/F53NnVje1oqq9f9AEC7dr44O9uzbNl+i99b00oTvRJdKxXu3LqHD9Vh6tWravF73fFYFwBk9RcAODnZ0amTP0uXxuY53VLTyiudQLQSb+/e1URknudsZaFaNReL36/H4G+oLg7s3XKCePO0jm7dAoiLS2HfvmSL31/TSgudQLQSb+Pv3wFQoZa7VaY22js6MyOsDuMvNOSPcNN03m7dAgD0MJam5aITiFbindq+B4BLng2tds/A/nWoecGVCnFvAeDvX4UGDTxYujTWajFoWkmnE4hW4iUdOI0AVOpptXtWvO9VpnCE84uWk5FpauvePYB16w6TmpputTi0suF2y7kDfPTRRwQEBFC/fn1WrFiR0z5kyBCqVatG48aNrfER8qQTiFbiBaY70t+xGo3rB1ntntVqdWB6haP8tDGRTeadbbt1CyQ9PYs1aw5ZLQ6t9CtKOffo6GjmzJnDnj17WL58Oc888wxZWVkADB48mOXLl1v98+SmE4hW4j2RXYMZjesSEGDd2kDtgtyIuHSeFRsWANC2bW1cXOz1MJZ2S4pSzn3x4sU8/PDDODo64u/vT0BAAOHh4QC0a9fO8HpZRq1E17RCU0ccSO2USmCgdf+ytBpwJ9mRv2K7dgw88wCOjnbcdVcdli3bj1KqVNQq0q6ae+JF4i/ntzvE7fFxCqFfjfEFvqco5dyPHTtGWFjYNeceO3asGD9B0egeiFai7d27Guez6/jq9CmrTOHNrdfTU3AXOw5siifBPHu3W7cADh1KISZGT+fVCqco5dxLepl33QPRSrSIFT+SiSLTs5LV/+I4VqhM5xpunE2EP7elMbSLY67pvLE0aFCyN/vRrnWznoKlFKWce2HONZLugWgl2ontuwBI8www5P4fjGjGHxmhZOx/DwBf38o0auTJ0qV6PYhWOEUp596rVy/mzJlDWloacXFxxMbG0qpVyakKrROIVqIl70sCIKtSN0Pu73L/cAA6R6/BPPmFbt0CWL/+MBcu6Om82s3lLufesGFD+vXrl1POfcmSJQA88cQTJCcnExAQwBdffJEz1TcoKIh+/frRqFEjunbtysSJE7G1tQWgf//+tG7dmpiYGHx8fPjhhx+s/tmkPNX2CQ0NVZGRkUaHod2CgbU9WRV/jncmb+Tpp1sYEsOzlb2J4xLvbDlDq/qwatVBOneexZIlD3PvvfUNiUkrnL1799KwofUWoJYFef2ZichWpVTo9e/VPRCtROtiU4XhNWrg72/cdrJSS/H32bP8uXETAHfeWRtXVwdd1kQr93QC0Uq0gak+PBdaHT8/4xJI6H3BZKBg3csA5um8ujqvpukEopVYKWeOcDwpm+Tqitq13QyL474R3+KMDfH/HuDMeVNbt24BHD58lv/+SzIsLk0zmk4gWom17q/v8WEd3x8/h5OTcTPO3dz9aFelEusPnWPlDlOPo1u3QEBX59XKN51AtBIrfpvpmYPy9DQ4Erj/Li8GZtbkwN5vAKhd241GjTx1AtHKNZ1AtBLrTEw8AKmu1iuimJ8uI59mNAF0PfAbVx576Om8WnmnE4hWYiUfPY8NYOvaxehQ8G02nPPVUzn/dzLRh01t3boFmKvzxhkbnFbiWaKce37XnDBhAgEBAYgISUmWfUanE4hWYiWevIy3jSMN6vgbHQo2tna86BJDj+3/8fuWI4BpOq+Li70extIKZIly7gVds02bNqxcuRJfX1+LfzadQLQSa4BjNUZ7+xg6hTe3Fvf4c54sUv8dBnBDdV5Ny4slyrkXdM1mzZrh5+dnlc+miylqJVb31OpcaJ5KQglJID1GfMrzkzuSsm47Fy+Ds5NpGGvJkhhiYpJ1ccUS7hRjSCP65m+8BY40wpPRBb7HUuXcb3ZNayiwByIiPiLyiogsFpEIEVkvIpNEpIeIFLn3IiLPiUiMiOwRkU/NbfYiMlNEdovIXhF5I59z/UVki4jEisivIuJQ1Hi0kiP1QjKbEy9ypHKmoWtAcvMNbE9rl4psjDnHut2mtivVefUmU1p+LFHOvaSUec+3ByIi04GawB/AJ8BJwAmoB3QFRonI60qp9bdzYxHpCPQGgpVSaSJSzXyoL+ColGoiIs5AtIjMVkoduu4SnwBfKqXmiMgU4Alg8u3EopU8G9fN4G4VzptJdRlr4BqQ67Vp6sYnG4+yetNMurUclFOdd9my/bz0Umujw9MKcLOegqVYqpx7SSjzXlAv4nOl1N1Kqa+VUhuVUvuVUlFKqQVKqeeADsDxItx7GPCxUioNQCl10tyuABcRsQMqAOnAudwniinVdgLmm5tmAvcVIRathDmweQ0AqlrJGhbq+WwPdtCae0//nNOmp/NqBbFEOffCXNMaCkogXUXEJ7+DSql0pVRRpp/UA9qah6HWiUhLc/t8IBVIAI4A45RSp687tyqQopTKNP8cj6m3pJURZ2IOAXC5cgNjA7lOWJ8vCHarQLPIs8SdMLXp6bxaQSxRzj2/awJ8/fXX+Pj4EB8fT3BwMEOHDrXch1NK5fkCvgSOAusx9RY88ntvAddYCUTl8ept/vo1IEArIM78fRvgZ8AeqAbEAHWuu64nsD/Xz7WA3fnE8BQQCUTWrl1baaXDK3f4KEC98NpPRodygzWtG6nnHH3U+MXnlVJKXb6coVxcxqphw/4wODLtetHR0UaHUOrk9WcGRKo8/n3NtweilBoB1AbeBoKBXSKyTEQeE5GKhUxOnZVSjfN4LcbUa1hgji8cyAY8gEeA5UqpDGUa1voXuL4OfRJQ2TzMBeBDPsNpSqmpSqlQpVSoZwkoiaEVTuKJy9QQBwJqBRodyg0iawvfpMWTuOlZQE/n1cqvAmdSmf9xX6eUGobpt/zxwAggsRjuvQjTcwxEpB7ggCkxHAE6iYkLEAb8d31cwBqgj7lpEHDtxGqtVHu2ohdTfOqWmBlYuXV+8RUAMtatIcM8iNqtWwCHDqUQE5NsYGSaZl2FmoorIk2AMcBETA+13yyGe08D6ohIFDAHGGRODBMBV0xDXBHAdKXULnMcS0XkylSDkcBLIrIf0zMR6+/nqFlM6IXKdAmsiI9PJaNDuUFI2GCaOrgQsTuFTXtNbVem8y5bpqfzauVHQdN4A4GHgf5AFqZ/5O9WSh0sjhsrpdKBR/Nov4BpKm9e53TP9f1BTM9OtDJoYUIKfo2FWiUwgQC0q1+JKbtPsGzzOto1aY+vb2UaNvRg2bL9jBihp/Nq5UNBPZAVmNZ9PKSUaqKUGltcyUPTCrJv3zr6pO3m1+QzeHo6Gx1OnlrcH4wXjtSMei+nrVu3ANatO0xqqp7Oq5UPBT1Er6OUGqWU2g0gIpVExP3Ky3ohauXNjo2mx1kOnq6GrK4tjIdenUacYxseSbhA4hlTW/fugaSnZ7F6tZ7Oq5UPN30GIiJPi0gisAvYan5FWjowrfw6GbUTgGz3ajd5p3GcXL3JaJWK278VWL7VNPNKV+fV8jNkyBCqVatG48aNb/ncrVu30qRJEwICAnj++edzZvq9++671KxZk5CQEEJCQli6dGlxh31ThXmI/goQpJTyU0r5m191LB2YVn6djUsA4KJryZvCm9v82qnUPL6ZqM2mEhl6Oq+Wn8GDB7N8+fLbOnfYsGFMnTqV2NhYYmNjr7nOiBEj2LFjBzt27KB79+4FXMUyCpNADgAXLR2Ipl2Rcuw8AOLcwdhAbqJmn/tIJB27dbPJzja16em8Wl7atWuHu/u1I/8HDhyga9eutGjRgrZt2/Lff//dcF5CQgLnzp2jdevWiAiPPfYYixYtslbYN1WYKnVvABtFZAuQdqVRKfW8xaLSyrUnqlalh6snu31Kdg+k3b2j8bf9jN3bktl+AFoEXjudV5d3L2H2jYTzu4v3mhWbQL1PbuvUp556iilTphAYGMiWLVt45plnWL169TXvOXbsGD4+VytK5S7nDqbdB3/88UdCQ0P5/PPPqVKlyu19jttUmB7It8BqYDNXn4FstWRQWvnmd86Jtt4u1K5dMqfwXmFja0dH30qsO32eJVsOAVwznVfT8nPhwgU2btxI3759CQkJ4emnnyYhIeGG9+U1FHplYsmwYcM4cOAAO3bswMvLi5dfftnicV+vMD2QTKXUSxaPRNPM5saewae6XYlcRHi90G51mDbxJFmb/gePmsamu3ULYMKECFJT03Fx0dvUlBi32VOwhOzsbCpXrsyOHTuuac/KyqJFixaAqULvsGHDiI+Pzzmeu2x79erVc9qffPJJevbsaYXIr1WYHsgaEXlKRLz0NF7NGl46Gces86dKRQLpMWIcw6UWnROOczbV1NatW6C5Ou8hQ2PTSq5KlSrh7+/PvHnzAFNPY+fOndja2uY8FB8zZgxeXl5UrFiRzZs3o5Tixx9/pHfv3gDX9FgWLlx4WzO8iqowCeQRzM9B0NN4NQs7kRhDssqgiocDnp4uRodzU7XrtuHLEB/a7nFj5XbTcEPbtlem8+qyJppJ//79ad26NTExMfj4+PDDDz/w888/88MPP9C0aVOCgoJu2Cf9ismTJzN06FACAgKoW7cu3bp1A+C1116jSZMmBAcHs2bNGr788ktrfiSgEENYSil/awSiaQDhG2YD4Ozpgo1NyVxEeL2TbRXxX2eweftPPHjnwJzpvH/+GcuECarELobUrGf27Nl5thdmam9oaChRUVE3tM+aNavIcRVVvj0QEbmzoBPNK9Ot32fSyrTjOyIAUFWrGhxJ4UU0aEAYW3Bc/QlXnnnee289Dh8+S1TUyYJP1rRSrKAhrAdFZKOIjBaRHiLSSkTaicgQEZmFaa/0ClaKUysnUg6a9nm+7FZ61qr2eGIinmLP/vATxJifd/bsWQ+AJUtiDIxM0yzrZhtK9cC0tWxf4H3gJSAQ+FYp1U4pFWGVKLVy434PF/ZxJ7aVOhgdSqHZOzjT0asSqxPOsTjctAiyRg1X7rijJkuW7DM4Ok2znJttKHVGKfWdUmqwUuoepdR9Sqk3lFL/WCtArXzxSMoiwMWRmtWbGB3KLWnV3otTKoMz4cNz2nr1qk94+DESEs4bGJmmWU6hNpTSNGv5bXsyP1Y4Su1aJX8Kb273vDASgKz1a7lkrtfQq1d9AP74Q/dCtLJJJxCtRPn+0El+ykikZs3SlUAa3/Eo63yCGX3Jm/XmahlBQZ74+1fWw1hamVXQLCx7awaiaQDH09OpXskeLy9Xo0O5ZUH3VMA1qhLro8IBU8mJXr3qs3LlQb3JVDlniXLuAN988w3169cnKCiI1157rThDLpSCeiDHROQ7EekkeiK7ZgVpaakkZmdQpbI9NWqUvgRyuH093so6gNvGq5V/evWqz+XLmaxcqTfzLM8sUc59zZo1LF68mF27drFnzx5eeeWV4gy5UApKIA0xrTh/GzgqIuNF5A7rhKWVRzu2LyYThWtVR+ztbY0O55YF3juWrzjMf6v3ccS8/KNt29q4uTnq6bzlnCXKuU+ePJnXX38dR0dHAKpVs/4GbPmuRFdKJWOqxPutiHhjmso7XkSqAXOUUqOsFKNWThzYtg4Au6qlr/cBULFyLTq4V2TtofP8GZHOsB4O2Nvb0r17IL//vo+srGxsbfVjR8O8OAV2HCjea4bUhfH/u61Ti1rOfd++fWzYsIFRo0bh5OTEuHHjaNmy5e1/lttQqP+blVLHgR+AycB5YKglg9LKp0CHc6TSmVqNGhgdym1r3dqTuKzLxIWPzGnr1as+p05dJDz8WAFnauVJcZRzz8zM5MyZM9ceFOkAACAASURBVGzevJnPPvuMfv36WX0nzAJrYYmIE3Av0B9oAyzHVFjxL8uHppU3LidO4YwdJ21Lzyr063V5+Vne+vMFWLmA9Le+xMEeunYNwM7OhiVLYmjdupbRIZZft9lTsITiKOfu4+PDAw88gIjQqlUrbGxsSEpKwtPT02qfo6BZWL8AR4CHgF8AX6XUIKXUMqVUlrUC1MqPTZuP8hL/cS67udGh3LZWHZ+nbYVKOOzPZp15Om/lyk60b+/L4sX6OYhmUhzl3O+7776cIa99+/aRnp6Oh4d1d8EsaAhrBVBXKdVHKTVfKXW5OG8sIs+JSIyI7BGRT81t9iIyU0R2i8heEXkjn3NniEiciOwwv0KKMzbNGJtjUpjGMWp5lu4C0PMGN+T9U434Z8fSnLbeveuzd28SMTFJBkamGcUS5dyHDBnCwYMHady4MQ8//DAzZ860euXngh6iz7TUTUWkI9AbCFZKpZkfzIPpQb2jUqqJiDgD0SIyWyl1KI/LvKqUmm+pGDXrSzqbgZedPT41KxodSpGce7AT1Sdvo1PsOJTqjgg88EBDnn9+Ob/9tpc332xrdIialVminLuDgwM//fRTkWMrCqOmhAwDPlZKpQEopa7UvFaAi4jYYar0mw6cMyZEzdoSUzOpUcEeL6/SnUDqdhhDD4cIvvh+R86kn5o1KxEW5sP8+dHGBqdpxcioBFIPaCsiW0RknYhcmXs2H0jFVAH4CDBOKXU6n2uMFZFdIvKliDjmdyPzdryRIhJ56tSpYv0QWvE6kZZBNVe7UrkKPTcbWzu8atmxMuksCzZcndvfp09Dtm8/wcGDZwyMTtOKz00TiIg4i8jbIvKd+edAEbnp7u0islJEovJ49cY0dFYFCANeBeaaV7u3ArIAb8AfeFlE8pqS8wbQAGgJuAMj83gPAEqpqUqpUKVUqDVnJ2i3Jiszg4wsRdVSugr9eq0faMRFsmHD1RnvDzzQEIAFC/YaFZamFavC9ECmA2lAa/PP8cAHNztJKdVZKdU4j9di8zUWKJNwIBvwwLT/+nKlVIZ5WOtfIDSPayeYz00zx9eqEJ9DK8FOJG7nKB14rpMfjo433Wm5xHtw5DRcsSV+7T4OJ5ra/P2r0Ly5lx7G0sqMwiSQukqpT4EMAKXUJaCoj/oXAZ0ARKQe4AAkYRq26iQmLph6KDes7xcRL/NXAe4DbnzCpJUqCbGmVegn3crGJpeVq/rT2bMSfx0+y7yNF3Pa+/RpyJYtx4iP14/2tNKvMAkkXUQqYHrAjYjUxdQjKYppQB0RiQLmAIOUaQnlRMAVU0KIAKYrpXaZ77vUXFIF4GcR2Q3sxtRzuWmPSCvZ9mzYQF92EHmh7BSB7vdgXcZm18Pp8Js5bQ8+2AjQw1ha2VCYBPIOphXotUTkZ2AVUKS6wUqpdKXUo+YhreZKqdXm9gtKqb5KqSClVCOl1Ge5zuluLqmCUqqTUqqJ+fxHlVIXihKPZryE6EPMJ5FEVXaeU3V581MG2XjTO3oTKeb/Q+vVq0rjxtX0MFY5N3ToUKKji///gQ4dOhAZGVns183PTROIUupv4AFgMDAbCFVKrbVsWFp5cz7BNKST7Vp2Hmd51OrIkUbJrFyUwh/hV2sU9enTkH/+OcKJE/r3nvLq+++/p1GjRkaHQVZW0YqKFGYWVnPAF9PU2uNAbRGpa16roWnF4sypNJywwcuzhdGhFKsfa55hyKl97Fx3da+GBx9shFKwcKEexioPUlNT6dGjB02bNqVx48b8+uuv1/QUXF1dGTVqFE2bNiUsLIzERNOsiwMHDhAWFkbLli0ZPXo0rq6m2Ylr166lZ8+rE2GfffZZZsyYccN9hw0bRmhoKEFBQbzzzjs57X5+fowZM4Y777wzp5TK7SrMENYkYDMwFfgO2ITpucU+Ebm7SHfXNLPkMxl42TrgV9v95m8uRbq9YUoc2cvnceGSqS0oyJMGDTyYM2ePgZGVTx06dLjhNWnSJAAuXryY5/Er/zgnJSXdcKwwli9fjre3Nzt37iQqKoquXbteczw1NZWwsDB27txJu3bt+O677wB44YUXeOGFF4iIiMgpoHgrxo4dS2RkJLt27WLdunXs2rUr55iTkxP//PMPDz/88C1fN7fCJJBDQDPzWooWQDNMD7k7A58W6e6aZuaaYUuQc4VSvwr9eqHthxHmXJHVO0/nDGOJCI880pgNGw5z9OhZgyPULK1JkyasXLmSkSNHsmHDBtzc3K457uDgkNOjaNGiBYcOHQJg06ZN9O3bF4BHHnnklu87d+5cmjdvTrNmzdizZ881z1weeuih2/w01yrMMFQDpVTOr0pKqWgRaaaUOqh3utWKy3cuDbgcfJnEUr4KPS93h1VlzOpDRK57h4fbjwGgf/8mjB69ll9/3cMrr/yfwRGWH2vXrs33mLOzc4HHPTw8Cjyen3r16rF161aWLl3KG2+8wd13XztwY29vn1ME0dbWlszMzAKvZ2dnR3Z2ds7Ply/fWOc2Li6OcePGERERQZUqVRg8ePA173Nxcbnlz5GXwvRAYkRksoi0N78mYRq+csS8NkTTikqS7LnoocpcDwSg68jh2ADV/plNqvnvcECAOy1bejN7tl7CVNYdP34cZ2dnHn30UV555RW2bdtWqPPCwsL47bffAJgzZ05Ou6+vL9HR0aSlpXH27FlWrVp1w7nnzp3DxcUFNzc3EhMTWbZsWfF8mOsUJoEMBvYDLwIjgIPmtgygo0Wi0sqV48ejaZO6lbkpZ3ByKntzM1rf/QrH67Xi5WPV+TPXbKz+/RuzbVuCLvFexu3evZtWrVoREhLC2LFjeeuttwp13vjx4/niiy9o1aoVCQkJOUNftWrVol+/fgQHBzNgwACaNWt2w7lNmzalWbNmBAUFMWTIENq0aVOsnymHUqrcvFq0aKG0kuf338YoQI3qFmh0KBZz/OX/U4p71LuTp+a0HTt2Tom8q955Z41xgZVx0dHRRodw21JTU1V2drZSSqnZs2erXr16WeW+ef2ZAZEqj39TCzONN1BE5otItIgcvPKyTDrTyqPjUaYuvbhXMTgSyzn34EN0JZK0OR9x0TyM5e1dkQ4d/Pjll91W38taK/m2bt1KSEgIwcHBTJo0ic8//9zokG5Q2GKKk4FMTENWPwKzLBmUVr6cPXgUgAw3H4MjsZz6rZ8n0fEyayKSWJZroXD//o2JjT3Ntm0JxgWnlUht27Zl586d7Nq1i/Xr1xMQEGB0SDcoTAKpoJRaBYhS6rBS6l3MhRA1rTicP27aHyOrUtnemfieFu5svnief9Zfnf3+4IONsLe34ZdfdhsYWdmme3eFd6t/VoVJIJdFxAaIFZFnReR+oNrNTtK0wvJIU4RSiSpV7jQ6FIvqNvJJAGz/nMh5c4Fed/cKdOsWyOzZUWRlZRdwtnY7nJycSE5O1kmkEJRSJCcn4+TkVOhzCjPl5UXAGXgeeB/TMNZjtxWhpuVhcGUPnvP0ZZ5f2f69pH2vt7mjwmf8HZlEk42KQZ1Nc/8feyyYJUti+Pvvg3TtWvKGKUozHx8f4uPj0buRFo6TkxM+PoUfSi5MAvFTSkUAF4DHAUSkL7DltiLUtOs4nrIl2yMDb++ytwbkeoO61+Tcb1mcjH4HOpsWFd57b32qVq3A9Ok7dAIpZvb29vj7+xsdRplVmCGsNwrZpmm3pffOGF46d6DU74VeGA989gWv2fjzyI4VnDhtanNwsOWRR5qwaNF/nD59ydgANe0W5JtARKSbiHwD1BSRr3O9ZmCakaVpxWL7pQuctc8qk6vQr1fdvxspLc+waWEqP66+Ws798cdDSE/PYvZs/TBdKz0K6oEcB7YCl81fr7yWAPdYPjStPLhwPolTKoPKVexxdi47uxEWZHrdy/RN2cOZ9QNy2po186Jp0+rMmLHTwMg07dbkm0CUUjuVUjMw7Yk+M9drgVLqjPVC1MqyrRELUYBL1cLP/Cjt+o2dSgVsSFy6iX3xV9sffzyEyMjjREWdNC44TbsFBQ1h7RaRXcA2Edl1/cuKMWpl2KFd/wJgU7WSwZFYj49fGD1qVGbJwRRmrT6U0z5gQDD29jZMn77duOA07RYUNITVE7i3gJemFZnjmQR64YmD561vmFOadenfiGSVASse4kplbg8PZ3r2rMdPP+0mI6NoW41qmjUUNIR1+MoL03OQJubXJXObphVZsH0mi2mOY82ytZXtzQwcsxAvcSBpzWHW5Hrs8cQTzTh5MpXFi2OMC07TCqkwxRT7AeFAX6AfsEVE+lg6MK18qJhoWpKd5dDS4Eisq4KrB3/1b8TEE8Gs2jI3p71r1wB8fd2YPDmygLM1rWQozDqQUUBLpdQgpdRjQCvgbcuGpZUXY1fE0Ug24OdV/hZ7uYwYjI2yYdCeLzlz3tRma2vDU0+1YPXqOL1PiFbiFSaB2Cilck8LSS7keZp2U4kp6ShbRc2aZX8NyPX8Q1/gU6//eHRqFDNXpee0P/FEM+ztbZgyRfdCtJKtMIlguYisEJHBIjIY+BNYWtQbi8hzIhIjIntE5FNzm4OITDfPANspIh3yOdddRP4WkVjz17K7kUQZdzI1k+pO9uViFXpepI0zkekXOLH6oZy26tVdeeCBhsyYsZOLF/Wu0VrJddMEopR6FfgWCAaaAlOVUiOLclMR6Qj0BoKVUkHAOPOhJ833bAJ0AT43VwK+3uvAKqVUILDK/LNWCiWmZVDN1a5crELPyyPjZuKCLceX/MP2/Vfbhw0LJSXlMr/+qvdM10quwjxEHwGEK6VeUkqNUEotLIb7DgM+VkqlAeQaImuEKSFcaUsBQvM4vzcw0/z9TOC+YohJszKlFAlZ6bi72ePq6mB0OIao6duSB2pXYeHhM8xctj6nvV07Xxo29NAP07USrTBDWJWAFSKyQUSGi0j1YrhvPaCtiGwRkXUicmUKzk6gt4jYiYg/0AKolcf51ZVSCQDmr/nWAReRp0QkUkQidUnnkuX4iZ0MpibBdcrPIsK89BrRmQtk4fz7EM6lmtpEhGHDQomIOE54+DFjA9S0fBRmCOs98zDTcMAbWCciK292noisFJGoPF69MZWRrwKEAa8Cc0VEgGlAPBAJjAc2UsTCjUqpqUqpUKVUqKenZ1EupRWzhP1rmUgjQprVNDoUQz3w3Cze9/Dl0ejKzFh59ZnHoEEhVKrkyJdfbjYwOk3L363MpjoJnMA0C+umO/8opTorpRrn8VqMKUksUCbhQDbgoZTKNA+ThSilegOVgdg8Lp8oIl4A5q+6eFApdPlQFJlkc8LZzehQDGVja8ejrwbQ6JgHbrFPcWXzvEqVHHnyyebMm7eHI0fOGhukpuWhMM9AhonIWkzPJjyAJ5VSwUW87yLM+6qLSD3AAUgSEWcRcTG3dwEylVLReZy/BBhk/n4QsLiI8WgG+Hfldhz4m+2nKxgdiuE8nprCJudkTn63mrW5Ks0991wrAL75Ru/fppU8hemB+AIvKqWClFLv5PMP+q2aBtQRkShgDjBImTYtroapeONeYCQw8MoJIvK9iFx5oP4x0EVEYjHN1vq4GGLSrOzC8RQUYFO5fJUxyYtr5QB+9ovnzX1HmffnpJx2X9/K9OnTiO++28b582kGRqhpNyrMM5DXlVI7ivOmSql0pdSj5iGt5kqp1eb2Q0qp+kqphuYhsMO5zhmqlIo0f5+slLpLKRVo/nq6OOPTrOPsyUvYIdSu1dboUEqEfh8OIxOF44IPOZprvsdLL7Xm7Nk0pk8v1r+GmlZkekW5ZpgzZzKobmNPXb/imNhX+rXrPZq7K7vx67ZTTPz9UE57q1Y1adOmFuPHbyYzM9u4ADXtOjqBaIY5eT4DLwcHvL3L5yLCvDz0dHMSVDp2i3pwIdf26C+/3Jq4uBTmzt1jXHCadh2dQDTD9HOoxqCaHuW2jEleBn2wgo5ObnhvSufbFVfrY/Xu3YCgIE/Gjt1AdrYyMEJNu0onEM0wQzJ9eLSJJxUrOhodSolha2fP9++34pmkutSOGUimeV8pGxth1Ki2REefYtGi/4wNUtPMdALRDHHp0lmOnMnktJse07+ezzM/keVxmarT9zB3/dXeRr9+QQQGuvPBB+tRSvdCNOPpBKIZIjJyEX5qAz/EnzE6lBLHwbkanzU7yl0xe9j6x6CchYW2tja8+WZbtm8/wdKlea2v1TTr0glEM8TBHRsAsKuqH6Dnpf/EH6iILYd/+ZOV26+2DxjQBF9fN95/X/dCNOPpBKIZIuWAqXa5jYeewpsX38AOPB5cjQUnTvP7vBdzeiH29ra8+WZbtmw5xh9/7DM2SK3c0wlEM8S5I6byZeLRwOBISq6np43HGVsSZ/7Eqly9kMcfDyEgwJ0331xNVpZ+hqQZRycQzRDnT1wAoKp3e4MjKbkatejHkCbV2J5wifWr3r+mF/LBBx2JijrJL7/sNjZIrVzTCUQzRJeKroxzCKRRveZGh1KivThnMlGOYby25o9reiF9+wbRrFkNRo9eS3p6lnEBauWaTiCaIdqoioyo7kXNmvohekHqNOrNhUcuYLO8IiuWjc7phdjYCB99dBeHDqUwdepWY4PUyi2dQDRDxBxN44jbBb0KvRAuvzGapmzkwNcT+TP8avvdd9elQwc/xoxZR0rKZeMC1MotnUA0Qzy0fx+vnDqkV6EXgndgD3q3cGfhydOs+6VPzup0EeHzz+8mKekiY8asMzZIrVzSCUQzxInMdKq62RsdRqnxwrx5VBN7tkz7i++WX62R1by5F0OHNuebb8LZu/dUAVfQtOKnE4hmdXGHt3KeLNyqOhgdSqlRy681L3T3Y8OF88TPaXdNpd6xYzvh4mLPiBEr9OJCzap0AtGsbteW3wFwqupscCSlywu/rKeeXQUSFx1m3IKre6h5errw7rsdWLHigF5cqFmVTiCa1R3fHQGAraeHwZGULi6VajDr4458dyGYB7f35nDi1WPDh7ekYUMPXnhhOamp6flfRNOKkU4gmtXVt7nEzzShUkATo0MpdUJfXEzGHedwm6T46IfpOe329rZMmdKTuLgU3nlnrXEBauWKTiCa1fleSOMRvKno3cHoUEodG1s7/vvkPkIubeHcpFdZmmtab7t2vjz9dAu+/HIzkZHHjQtSKzd0AtGs7vB/Z9nokIyvTzOjQymVgtuPZFBINWYfS2btjJ5cSrt67JNPOlO9ugtDhy4hI0OvUNcsSycQzeombUtgcPYe/HyrGh1KqTXy9z+obePIsu/XMPaXQzntbm5OTJrUg507E/nkk3+NC1ArF3QC0awu8XwmXk72+PhUMjqUUquGTzPefK4FURkXSZ1wJ9ty7S91330NePjhxrz33jo9lKVZlGEJRESeE5EYEdkjIp+a2xxEZLqI7BaRnSLSIZ9z3xWRYyKyw/zqbtXgtSI5fjmd6q52VKigFxIWxdPj/6V/zapkbLdh5ZLXSM+4emzSpO54ebkyYMACPStLsxhDEoiIdAR6A8FKqSBgnPnQkwBKqSZAF+BzEckvxi+VUiHm11KLB60Vi8yMdI5npePhrhcRFocvNs7nmyoBvPztJj6Yc7W3UaVKBWbOvI/Y2GRefvkvAyPUyjKjeiDDgI+VUmkASqmT5vZGwKpcbSlAqCERahaxK2o5l8mmkoeugVUcatTuwP7P3Nkem0mFqW2IiLl6rGNHf1555f/49tutLFiw17ggtTLLqARSD2grIltEZJ2ItDS37wR6i4idiPgDLYBa+VzjWRHZJSLTRKRKfjcSkadEJFJEIk+d0rWCjJZ6NIL1tKJF8/z+s2q3KnDIfMb77Oetfw7x25fdOZt69dj773ekVauaDB68iJiYJOOC1MokiyUQEVkpIlF5vHoDdkAVIAx4FZgrIgJMA+KBSGA8sBHIzOPyk4G6QAiQAHyeXxxKqalKqVClVKinp2dxfkTtNjgkHqAtVRCf2kaHUqZ8uGEJde0qMOu7Vbz+9ZycfUMcHe2YP78vDg62PPjgXP08RCtWFksgSqnOSqnGebwWY0oSC5RJOJANeCilMpVSI8zPNXoDlYHYPK6dqJTKUkplA98BrSz1ObTidXz7fmaTwOkKAUaHUqbU9mvNuEkDSMrO4NBnw5jw+4WcY7VquTF79oNER5/iySd/1wUXtWJj1BDWIqATgIjUAxyAJBFxFhEXc3sXIFMpFX39ySLilevH+4Eoy4esFYd/tyXwCLtw82pndChlTq8nv+PtrnVYfiaFzAnNrnke0qVLXT74oBOzZ0cxduwG44LUyhSjEsg0oI6IRAFzgEHK9GtRNWCbiOwFRgIDr5wgIt+LyJUH6p+ap/ruAjoCI6wbvna7zpxKo6rY0ywoxOhQyqQ3ft/Ld03q8Ozffhz/7QGOJ+c69sadDBwYzNtvr+GXX3YbF6RWZkh56s6GhoaqyMhIo8Mo1+71qMLh8+lsOn0GFxc9ldcSLqTsp0LY4yQeFj59sRkfvfsVFcyT3tLSMrnnnp/YtCmelSsH0ratr7HBaqWCiGxVSt0wI1avRNesKiE1Ey9ne508LMi1cgDHFvyP3hk7+OOzqbz01a9kZ5uOOTrasWDBQ/j7V6Z37zns3HnC2GC1Uk0nEM2qjqal4VVZJw9Lq91oAK9+1o1jWelEvjOUVyf9kzMzy929AsuXP4qLiwNdusziv//09F7t9ugEolnNqaRYtqjWDL3T2+hQyoV+I35l/PNhbL+cyo5RPXnvp/05x/z8KrNq1WPY2Ah33fUjBw6cLuBKmpY3nUA0q4mLWowfFbD1dzc6lHLj6a/+5ZOHg1h97iyXRrfnqyVXn6rXq1eVlSsfIy0tk44dZ+qFhtot0wlEs5r969fyKXHE2FU2OpRy5eXZu/m6Vz1eP9SAJ7/uzhcLE3KONW5cjVWrHiMtLYu2baezfXtCAVfStGvpBKJZTWxELCPZx2nn+kaHUu48tziGi2+mY7OqEikjWjH2l6vDWU2b1mDDhsepUMGeDh1msn79YQMj1UoTnUA0qzlz7AI2QNPQvkaHUi7VHLuB2f1P8P7heNYMC+WNKatzZmfVq1eVf/8dgrd3Rbp0mcWMGTuMDVYrFXQC0azm5Mk0ato4EtKkgdGhlFuP/7KbT/s1Yt25c/z5wr2M+OiLnC1xfXwq8e+/Q2jbtjaPP76YV1/9i6ysbGMD1ko0nUA0qzl+NoNaTg64u1cwOpRy7dVf9zD1jU4cTk9j/ug3eO+tQZwwT8Jyd6/AsmUDGD68JePGbaJnz9mcOpVa8AW1cksnEM1q4i+l4+PmgKnwsmakxz9cyS8z/4e72DNk/H6Oft6RdbtMWxra29syYUJ3pkzpwZo1cTRtOoU1a+IMjlgriXQC0awiPe08e2jDex10GfeSosdjE1ixbx5+LbIJ/dCBbY8GMvbn7TnPRZ5+OpTNm4dSqZIjd931I6NGrSItLa/dFbTySicQzSp275yHY5Ydl31djQ5Fy8W7Tjds1v/FigeP8/Luw/w4uA1vvfwgh8wVTkJCahAZ+RSDB4fw4Yf/0Lz5VDZvjjc2aK3E0AlEs4qdf8znOfay3V6vASlp7Bwq0nX+bqaPvoe07Gw+Gb+Az7p48fXC/8jOBldXB6ZN682ffz7C+fNp/N///cALLywjJeWy0aFrBtMJRLOKIxGxTOAIF6s1NzoULR+D3lvO2n0rGejnyaSoEyx6uDXLPnuQiH2mMa3u3QPZs+cZhg0L5ZtvwgkM/IZJkyLIzNQztcornUA0qzh55DxO2HBX5yeNDkUrgF/dO5kRd5Jpb3bm5co+9Hg9lXqDW/P2B69x4jRUrOjIxIk9iIx8iqAgT4YPX0pw8GQWLtxLdnb52RpCM9EJRLOKo4np1LV3ol49XUixNHh87N90PrSG+NcvMCXyNB+9PY6Xm3ry7qejOZUCzZt7sWbNIBYs6EdmZjYPPDCXkJApzJ27R68dKUd0AtGs4uC5NOq4OWJjo6fwlhaOFTzw+egfuv7zDk81rMaC+NN8MPIDXmrqyQcfPM/RU8L99zckOno4s2bdT0ZGNg89NJ+goElMmBDOuXNpRn8EzcJ0AtEs7kLqadIzFLWqOxodinYbmrZ6lEnRJ/jn36k8Ua8a846cZs3bP1O9V2uWf96XDXsuMWBAMFFRw5gz50Hc3Jx47rll1Kz5Bc8+u5RduxKN/giahegtbTWL27D8Q9p2W8+Wt9K44/01RoejFdGerXNJHv8hbZe7k5CUTUfbcLo1daf6o4MZ+OgH+HhCePgxJk6MYM6cKNLTs2jSpBoDBjThkUeaUKuWm9EfQbtFektbzTCOOzYAEO8bbHAkWnEIatGPdrN2kH5kLltfy6S2qwNfbzvGmy+Npa9fJd7o4M/2PT/x/rjeHDv2EhMndsfV1YHXX1+Fr+94wsK+Z+zY9ezalUh5+gW2LNI9EM3i3v8/X3ZvSmPU7lU0bRxkdDiaBYSvn8zitz7h9y2niEq/SAId8PTPYmlQMnsCa9Gk71f4u3vy27woFi+OITLyOAC+vm7cc09dOnTwo317P7y9Kxr8SbS85NcD0QlEs7ieVatw8Hwae9JSdR2scmDz35/gveYPvP7K4r6tsSwlCUdsaF3Rleb1K+LdrjHenUdyMq4yq1fEsHbtoZwH7oGB7rRr50vLlt6EhnrTpEl1HBxsDf5Emk4g6ARihIz0dLydXLnHx42fjpwyOhzNyg7GrmHZ16PYszqOzQfOsyMtldZU5l/uILvWJd6w34+Njz0VG/pyzvNOIrY3Y9vGA6ScMa1yd3CwJTi4Os2a1aBhQw8aNvSkQQMPatd20zP6rCi/BGJnRDBa+fHLjNdJUhnUDdb7oJdHdQI7MvybjTk/x0QvI+bPScTHJ1FlVzYLN5wh9uAlWH8YWI8rtgyv6sVb3b047SNMSzhHursn+8/68s6n9TmfqEBl4+xsT/36ValTpwq+vm74viCDEwAADXFJREFU+VXG17dyzvdubk7GfehyxLAEIiK/Alf2Nq0MpCilQszH3gCeALKA55VSK/I43x+YA7gD24CBSql0a8SuFd6+uYsB8L67h8GRaCVB/UbdqN+oW87P+4DdW+cQ/scsjkXEciz2LA0ynXCJcCFtqfAe4cB/wAYEqGnjwKgatXi8rgeJldIZF3sC+1MOnIh15pCLG+dUJc5m1SA+2ZvkpGpUdnHDq7oz1aq5UL26C9Wru+Lh4Uzlyk43vJyc9O/Tt6pEDGGJyOfAWaXUmP9v786D5CjPO45/fzOzh3ZXWmm1QqzuA8lGHAKFAmEF7GCVMFSMEHZcSuyYlEnwAamAKhVkq0ypnLIrBDuOk2AwNhSyYwfjYIywHRfgIpBKhC0BEpIsJK1OFlZanbva+5gnf/Q78uzszIodsTNj6/lUjab7fbunH73d0+/0273vK2kB8B/AlcAU4HlgvpkNZKzzBPBjM3tc0kPAFjN7cLjteBNWYfX39bFu9gLWnzjBU6eOeJODG5G21oNs+OV32LdxA8d3NnH0wCkOt/TwkeqJLE9OYmdLLwtbN9LP4HPYwyzgr5jOJlr5IzZSG4tTG09QWxanuizGPbOnsHjyWHb0dfKtgy1UVMQor4iRqIiTqIxz9Xvqqa2r4VB7H9sPtUNZAiXKIFGGyiqYMnUy5dU1tPeIk139xBLVxMqrSVSMpaysmtpxEyivHIMQxMqoKKukqnIMlVU1VFbWUDOmhsrKchLxGIlEjFhMxOMK77FB0+l5sZiKdg+xZJuwFJXIx4DrQtJy4HEz6wH2SWokqkw2ZKxzHfBnIWkdsBYYtgLJ14HPLubOdbvZ0981KH1JRS3fqp8HwNLm12keGHwBtGxMHV+vmwvA4ubXaEsOHkthRVU9Xx4/G4BLmjfRn1GZ/3n1ZL4wbgY9lmThoVeGxPXp6gbuHjuN48k+rm55bUj+qpppfLpmCgf7u1l65PUh+feOm8Unqiazo6+D5ce2Dcm/b9xcVoyp55XeU6w88Zsh+Q/UzmNZRR0v9pzkttadQ/LXVV/IbW1zWfx3HV55uBEbVzuD62/5EtySPX8B0NHTwe6dz7Fz80sc2bWDjpbjzKpNsDfeQd9b3XxiUz3tnQO0dw3Q2pOkrXcAWuJUHaukta2Hp1qP00mSbn7b/cpzL9Tzfip4ghP8PVuGbHcDV7GY8TxKE3/N9iH521nCAmr4Bge4izeG5B/kWqYzhi9rD2ttDzEghk6/H6hawrh4GV/q3cu/9TYhIPXtEeLg+Kspi8X4QucevttzaFBepWLsqr8KBKvaGvlJz1EE3F7VwKL7lrDss98bwR44s6JXIMA1wGEz2x3mpwIvp+U3hbR0E4mavPqHWQYASbcDtwPMmJHfYEYnJ1Yya1IZ5Rm9V08dm6B7dvT0yAX9FUzoG/y0yJTaOF0zo5Uu7KugY6BsUP75E+J0TY/yL+4dQ9Lst0cKcN5E0TW1m95kksv6hw4DW3+e6Gropq9/gEVWnZFr1DVA53ldWE8fV6hqcLbE+KlGZ30Xsa4+rtiVuT7UzBigs66TRHsfVzYOza+aNUDHhE4q2/q5au/Q/PK5vRy88ART1vxsSJ5z74byimouuvRmLrr05iF5c4Crh1l3KXAsTPf0dHDkSCOHDjWigTY29Z6k7uhRHmu6gr6uTga6uxno7mKgt5euaePYlBATmyfxlT2XMNA3QLKvn4G+JMmBJO3v7aexop1ZzeWsPjgTS9rpV9Lg1EUdvJXoZv5bCT7XdD5JAzPDkpBMGh2XdGDxGHPejHPTofGkflem3k8t7CAuMefNBB9sGRflhX8SEqcujYYgnrU/zuIT1QhomCSSdVlPkWdlVJuwJD0PnJ8la42ZPR2WeRBoNLOvhfkHgA1m9u9h/hHg52b2ZNrnTgrLXBDmp4dlLhkuHm/Ccs65kStKE5aZLR0uX1KC6AL1D9KSm4DpafPTgLczVj0KjJeUCFch2ZZxzjk3iordlclS4A0zSx8jcz2wUlJFeNJqHvDr9JUsumx6AfhoSLoVeLoA8TrnnAuKXYGsJHri6jQz2w48AfwG+AVwR+oJLEk/l5QaUOIeYFW4yT4ReKRgUTvnnCuNx3gLxe+BOOfcyHlvvM45595VXoE455zLi1cgzjnn8uIViHPOubycUzfRJR0BDuS5ej3R35+UmlKNC0o3No9rZEo1Lijd2H7f4pppZpMyE8+pCuRsSNqU7SmEYivVuKB0Y/O4RqZU44LSje1cicubsJxzzuXFKxDnnHN58QrknXu42AHkUKpxQenG5nGNTKnGBaUb2zkRl98Dcc45lxe/AnHOOZcXr0Ccc87lxSuQd0DShyTtlNQoaXUR45gu6QVJOyRtl/Q3IX2tpLckbQ6vG4sQ235JW8P2N4W0OknPSdod3icUOKb3pJXJZkltku4qVnlJelRSi6RtaWlZy0iRfwnH3OuSFhU4rvslvRG2/ZSk8SF9lqSutLJ7qMBx5dx3kj4fymunpOsLHNcP02LaL2lzSC9keeU6P4zeMWZm/hrmBcSBPUQjZJYDW4AFRYqlAVgUpscCu4iGhl4L/G2Ry2k/UJ+R9o/A6jC9GrivyPvxEDCzWOUFXAssAradqYyAG4H/IhrkeDHwqwLHtQxIhOn70uKalb5cEcor674L34MtQAUwO3xn44WKKyP/a8C9RSivXOeHUTvG/ArkzK4kGnJ3r5n1Ao8Dy4sRiJk1m9mrYfoUsIMcY8GXiOXAujC9Dhg6cHXhfBDYY2b59kRw1szsJeB4RnKuMloOfNciLxONwNlQqLjM7FmLRvsEeJlo1M+CylFeuSwHHjezHjPbBzQSfXcLGpckAR8jY5yjQhjm/DBqx5hXIGc2FXgzbb6JEjhpS5oFXA78KiTdGS5DHy10U1FgwLOSXpF0e0ibbGbNEB3cwHlFiCslc/CyYpdXSq4yKqXj7lNEv1RTZkt6TdKLkq4pQjzZ9l2plNc1wGEz252WVvDyyjg/jNox5hXImSlLWlGffZZUAzwJ3GVmbcCDwFzgMqCZ6BK60JaY2SLgBuAOSdcWIYasJJUDNwE/CkmlUF5nUhLHnaQ1QD/w/ZDUDMwws8uBVcAPJI0rYEi59l1JlBfwpwz+oVLw8spyfsi5aJa0EZWZVyBn1gRMT5ufBrxdpFiQVEZ0cHzfzH4MYGaHzWzAzJLAtxmlS/fhmNnb4b0FeCrEcDh1SRzeWwodV3AD8KqZHQ4xFr280uQqo6Ifd5JuBf4Y+LiFRvPQRHQsTL9CdK9hfqFiGmbflUJ5JYBbgB+m0gpdXtnOD4ziMeYVyJltBOZJmh1+ya4E1hcjkNC++giww8z+KS09vd1yBbAtc91Rjqta0tjUNNEN2G1E5XRrWOxW4OlCxpVm0K/CYpdXhlxltB74ZHhSZjHQmmqGKARJHwLuAW4ys8609EmS4mF6DjAP2FvAuHLtu/XASkkVkmaHuH5dqLiCpcAbZtaUSihkeeU6PzCax1ghng74XX8RPa2wi+jXw5oixvGHRJeYrwObw+tG4HvA1pC+HmgocFxziJ6A2QJsT5URMBH4JbA7vNcVocyqgGNAbVpaUcqLqBJrBvqIfv3dlquMiJoXHgjH3FbgigLH1UjUPp46zh4Ky34k7OMtwKvAhwscV859B6wJ5bUTuKGQcYX0x4DPZCxbyPLKdX4YtWPMuzJxzjmXF2/Ccs45lxevQJxzzuXFKxDnnHN58QrEOedcXrwCcc45lxevQJzLQdJ4SZ9Lm58i6T9HYTsfkPS+YfJvlnTvu7Cdr0q67mw/x7kUf4zXuRxCf0I/NbOLR3k7a4F2M/tqjvz/I/qDvqNnuZ2ZwLfNbNnZfI5zKX4F4lxu/wDMDeM43B/GdtgGIOkvJP1E0jOS9km6U9Kq0Gney5LqwnJzJf0idDL5P5Lem76BUEl9Brg7bOeajPz5QE+q8pD0mKQHFY37sFfS+0OngjskPRaWiYfltikao+VuAIt6Ip4o6fzRLDR37kgUOwDnSthq4GIzuwxOn+zTXUzU42kl0V9u32Nml0v6OvBJ4J+Bh4n+Onm3pKuAbwKnm5HMbL+iQYZyXYEsIfoL5nQTwmfcBDwTlvlLYKOky4jGPpmaunJSGAwqeDUs/+QIysG5rLwCcS5/L1g07sIpSa1EJ3OIuoW4NPSK+j7gR1E3RUA04NFINABHMtKeMTOTtJWo6/CtAJK2Ew1g9CIwR9K/Aj8Dnk1btwWYMsIYnMvKKxDn8teTNp1Mm08SfbdiwMnUFUyeuoDaHNtN3+bp7ZrZCUkLgeuBO4gGOPpUWKYyfKZzZ83vgTiX2ymioUHzYtFYDPsk/QmcHoN64Qi3swO4YCTblVQPxMzsSeCLRMOvpsynuL0Pu98jXoE4l4NF4zj8b7gZfX+eH/Nx4DZJqZ6Ksw2H/AywIttNdOAl4HKltYG9A1OB/5a0maiH2M/D6bEiLgA2jey/4Fx2/hivcyVO0jeI7ns8f5afswJYZGZffHcic+c6vwJxrvR9hWhck7OVoDSH73W/o/wKxDnnXF78CsQ551xevAJxzjmXF69AnHPO5cUrEOecc3nxCsQ551xe/h+O52GEGisEwQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i,epsilon in enumerate(epsilon_array):\n",
    "    nest.ResetKernel() # reset simulation kernel \n",
    "    nest.resolution = dt\n",
    "\n",
    "    # Current based alpha neuron \n",
    "    neuron = nest.Create('iaf_psc_alpha') \n",
    "    neuron.set(C_m=C_m, tau_m=taum, t_ref=0., V_reset=-70., V_th=1e32,\n",
    "               tau_syn_ex=taum+epsilon, tau_syn_in=taum+epsilon, I_e=0.)\n",
    "   \n",
    "    # create a spike generator\n",
    "    spikegenerator_ex = nest.Create('spike_generator')\n",
    "    spikegenerator_ex.spike_times = [50.]\n",
    "    \n",
    "    # create a voltmeter\n",
    "    vm = nest.Create('voltmeter', params={'interval':dt})\n",
    "\n",
    "    ## connect spike generator and voltmeter to the neuron\n",
    "    nest.Connect(spikegenerator_ex, neuron, 'all_to_all', {'weight':100.})\n",
    "    nest.Connect(vm, neuron)\n",
    "\n",
    "    # run simulation for 200ms\n",
    "    nest.Simulate(200.) \n",
    "\n",
    "    # read out recording time and voltage from voltmeter\n",
    "    times = vm.get('events','times')\n",
    "    voltage = vm.get('events', 'V_m')\n",
    "    \n",
    "    # store maximum value of voltage trace in array\n",
    "    maxVs.append(np.max(voltage))\n",
    "\n",
    "    # plot voltage trace\n",
    "    if epsilon == 0.:\n",
    "        pl.plot(times,voltage,'--',color='black',label='singular')\n",
    "    else:\n",
    "        pl.plot(times,voltage,color = cmap(1.*i/NUM_COLORS),label=str(epsilon))\n",
    "\n",
    "pl.legend()\n",
    "pl.xlabel('time t (ms)')\n",
    "pl.ylabel('voltage V (mV)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show maximum values of voltage traces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f68764a5750>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEKCAYAAADenhiQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZRU9Zn/8ffDJpsCgiIKCCSgBmgabBRFDAjuJiQxLrihxjGD0VHPT1F/uGTyO06cOZqYGY2GuOHEiBsCjlWiEBEYRTbZ3KO22oqAREFAEOjn98e9hU3bVV3dXbdvVfXndU6dqrpL3efbRffDd7nfr7k7IiIiudQs7gBERKT4KLmIiEjOKbmIiEjOKbmIiEjOKbmIiEjOKbmIiEjOtYg7gHzQpUsX79WrV9xhiIgUlKVLl37u7vvVtE/JBejVqxdLliyJOwwRkYJiZh+m26dmMRERyTklFxERyTklFxERyTn1uaSxY8cOKioq2LZtW9yhFITWrVvTvXt3WrZsGXcoIpIHlFzSqKioYO+996ZXr16YWdzh5DV3Z8OGDVRUVNC7d++4wxGRPKBmsTS2bdtG586dlViyYGZ07txZtTwR2U3JJQMlluzpZyVSgKZPh/LySD5ayaUJKS8vp3v37lRWVu6xvbS0lEWLFsUUlYjEYvNmOOssuOuuSD5eyaUJ6dWrFz169GD+/Pm7t7311lt89dVXHHHEETFGJiKN7m9/g2++gVNOieTjlVzyWHl5OYceeiiXXHIJAwYM4Nxzz2X27NkMHz6cvn37smjRIhYtWsTRRx/N4MGDOfroo3n77bcB+N3vfsfFF18MwKpVqxgwYABbt25l3LhxTJ06dfc1pk6dyrhx42Ipn4jEKJmE9u3hmGMi+XjTMsdQVlbm1ad/efPNNznssMOCN1ddBcuX5/aipaVw550ZDykvL+f73/8+r732Gv3792fo0KEMGjSI+++/n5kzZ/Lggw/y8MMP07ZtW1q0aMHs2bO55557eOqpp6isrGTkyJFcffXV3HrrrfzhD39g+PDhfPbZZwwePJiPP/6YFi1acNhhh/HEE08wYMCABhdpj5+ZiOQvd+jVC4YMgaefrvfHmNlSdy+raZ+GIue53r17M3DgQAD69+/P6NGjMTMGDhxIeXk5GzduZPz48bz77ruYGTt27ACgWbNmPPTQQ5SUlPDLX/6S4cOHA3DAAQfQv39/5syZQ9euXWnZsmVOEouIFJA334SPPoJJkyK7hJJLNmqpYURpr7322v26WbNmu983a9aMnTt3ctNNNzFq1CiefvppysvLGTly5O7j3333Xdq3b8+nn366x2emmsa6du2qJjGRpiiZDJ5PPjmyS6jPpcBt3LiRgw46CICHHnpoj+1XXnkl8+bNY8OGDTz55JO7951++ukkEgkee+wxzj777MYOWUTilkjAgAHQo0dkl1ByKXATJ07khhtuYPjw4ezatWv39quvvprLLruMfv36cf/993P99dezbt06ADp27MiwYcPo2rWr7qgXaWq++grmz4+01gLq0Aey6NCXrOhnJlIApk+Hn/4UXnwRqjSj10emDn3VXEREmpJkEvbeG8JBPlFRchERaSrcg/6W44+HiGcwV3IREWkqXn8dKioi728BJRcRkaYjkQielVxERCRnkkkoKYHw9oUoKbmIiDQFmzbBggWRTVRZnZJLgbnkkkt44403cv65I0eOpPpwbBEpIrNnw86djdIkBpr+peDcd999cYcAwK5du2jevHncYYhIthIJ6NABjjqqUS6nmkse27JlC6eeeiqDBg1iwIABPPbYY3vUMNq3b8+kSZMYNGgQw4YNY+3atQC89957DBs2jKFDh3LzzTfTvn17AObOnctpp522+/Mvv/zyPaaMSZkwYQJlZWX079+fW265Zff2Xr168Zvf/IZjjjmGJ554IsKSi0hOuQf9LY0wBDlFNZcsXHXVVSzP8ZT7paWl3FnLhJjPPfccBx54IM8++ywQzBd2zz337N6/ZcsWhg0bxq233srEiRP585//zI033siVV17JlVdeybhx47j33nvrHNutt97Kvvvuy65duxg9ejQrV66kpKQEgNatW7NgwYI6f6aIxGjlSvj000brbwHVXPLawIEDmT17Ntdddx3z58+nQ4cOe+xv1arV7prI4YcfTnm4FvYrr7zCGWecAcA555xT5+s+/vjjDBkyhMGDB/P666/v0cdz1lln1bM0IhKb1CzIJ53UaJdUzSULtdUwotKvXz+WLl1KIpHghhtu4IQTTthjf8uWLTEzAJo3b87OnTszfl6LFi2orKzc/X7btm3fOeaDDz7g9ttvZ/HixXTq1IkLL7xwj+PatWvXkCKJSBwSCRg8GLp1a7RLquaSxz799FPatm3LeeedxzXXXMOyZcuyOm/YsGE89dRTAHssaXzwwQfzxhtvsH37djZu3MicOXO+c+6mTZto164dHTp0YO3atSRT/+MRkcL05Zfw8suNNkosRTWXPLZq1SquvfZamjVrRsuWLbnnnnu45ppraj3vzjvv5LzzzuOOO+7g1FNP3d2c1qNHD84880xKSkro27cvgwcP/s65gwYNYvDgwfTv358+ffrsXsFSRArUCy/Arl2N2t8CmnIfKL4p97du3UqbNm0wM6ZOncqjjz7KjBkzIr9uIf/MRIrWxRfD00/D+vXQIrf1iUxT7qvmUoSWLl3K5ZdfjrvTsWNHHnjggbhDEpE4VFYGnfknnpjzxFIbJZciNGLECFasWBF3GCIStxUr4LPPGr2/BdShLyJSvFKzIDfiEOQUJZcM1B+VPf2sRPJQMgmHHw5duzb6pZVc0mjdujUbNmzQH80suDsbNmygdevWcYciIin/+Ae88kosTWIQY5+LmV0BXA7sBJ5194nh9hLgT8A+QCUw1N23VTt3EHAv0B4oB851901m1hK4DxhCULaH3f239Ymve/fuVFRUsH79+vqc3uS0bt2a7t27xx2GiKS88ELQod/IQ5BTYkkuZjYKGAuUuPt2M9s/3N4C+AtwvruvMLPOwI4aPuI+4Bp3f8nMLgauBW4CzgD2cveBZtYWeMPMHnX38rrG2LJlS3r37l2v8omIxC6RgH33hSOOiOXycTWLTQBuc/ftAO6+Ltx+ArDS3VeE2ze4+64azj8EmBe+fgE4PXztQLswSbUBvgE2RVMEEZE8VVkJzz0XDEGOaWmMuJJLP2CEmb1qZi+Z2dAq293MZpnZMjObmOb81cCPw9dnAD3C108CW4A1wEfA7e7+j2iKICKSp157Ddati62/BSJsFjOz2cABNeyaFF63EzAMGAo8bmZ9wu3HhNu2AnPCO0CrT4J1MfCfZnYzMJOghgJwBLALODD8/PlmNtvd368hvkuBSwF69uzZkKKKiOSXRALMgppLTCJLLu4+Jt0+M5sATPNgKNYiM6sEugAVwEvu/nl4XIKgc36P5OLubxE0oWFm/YBTw13nAM+5+w5gnZn9L1AGfCe5uPtkYDIE0780oKgiIvklmYSyMth//9hCiKtZbDpwHOxODq2Az4FZQImZtQ37TX4IfGfB+CoDAJoBNxKMHIOgKew4C7QjqBm9FXFZRETyx4YNsHBhbKPEUuJKLg8AfcxsNTAVGO+BL4DfAYuB5cAyd38WwMzuM7PUBGnjzOwdgsTxKfBguP1uguHJq8PPeNDdVzZWoUREYvf888GyxjH2t4BmRQZqnhVZRKQgnX9+MFLss88iHymWaVZk3aEvIlIs8mAIcoqSi4hIsViyBD7/PPb+FlByEREpHslkMAT5hBPijkTJRUSkaCQScOSR0KVL3JEouYiIFIX162Hx4thHiaUouYiIFINZs4IhyHnQ3wJKLiIixSGZDO7IHzIk7kgAJRcRkcK3a1cwBPmkk6BZfvxZz48oRESk/hYvDlaezJP+FlByEREpfIlEUGPJgyHIKUouIiKFLpmEYcOClSfzRFbJxcw6mVl/M+sTzkQsIiL5YO3a4M78PGoSgwzruZhZB+BXwDiCKfHXA62Brma2EPiju7/YKFGKiEjNZs0KnvNkCHJKpsXCngQeBka4+5dVd5jZ4cD5ZtbH3e+PMkAREckgkYCuXaG0NO5I9pA2ubj78Rn2LQWWRhKRiIhkZ+fOYP2WsWPzZghyStpozOwNM5tkZt9rzIBERCRLr74KX3yRd/0tkLlDfxzBqo7Pm9mrZnaVmR3YSHGJiEhtkslg3Zbj0zY0xSZtcnH3Fe5+g7t/D7gSOBhYaGZ/M7N/arQIRUSkZokEHHUUdOoUdyTfkVUjnbsvdPergQuATsBdkUYlIiKZrVkDr72Wd6PEUjKNFgPAzIYSNJGdDpQDk4Enog1LREQyeu654DkP+1sg830u/wacBXwBTAWGu3tFYwUmIiIZJJPQrRsMGhR3JDXKVHPZDpzs7u80VjAiIpKF1BDk008PljXOQ5nuc/lXADNrDpwK9Kp6vLv/LurgRESkBq+8Ahs35m1/C2TR5wI8A2wDVgGV0YYjIiK1SiahRQsYMybuSNLKJrl0d/eSyCMREZHsJBIwfDh06BB3JGllMxQ5aWb5s0iAiEhT9sknsGJF3o4SS8mm5rIQeDqcan8HYIC7+z6RRiYiIt+VGoKcx/0tkF1yuQM4Cljl7h5xPCIikkkyCQcdBAMGxB1JRtk0i70LrFZiERGJ2Y4d8MILQa0lT4cgp2RTc1kDzDWzJMG9L4CGIouINLqXX4ZNm/K+vwWySy4fhI9W4UNEROKQSARDkEePjjuSWtWaXFI3U4qISMySSRgxAvbJ//FUmRYLm2xmA9Psa2dmF5vZudGFJiIiu338MaxaVRBNYpC55vJH4KYwwawG1gOtgb7APsADwCORRygiIgUzBDkl09xiy4Ezzaw9UAZ0A74G3nT3txt6YTO7Argc2Ak86+4Tw5rQtVUOKwGGhLFUPXdf4DGC+c7KgTPd/QszM+APwCnAVuBCd1/W0FhFRGKXSECPHvCDH8QdSVay6XPZDMzN5UXNbBQwFihx9+1mtn94rUcIa0NhjWlG9cQSuh6Y4+63mdn14fvrgJMJalZ9gSOBe8JnEZHC9c03MHs2nHtu3g9BTslqJcoITABuc/ftAO6+roZjxgGPpjl/LDAlfD0F+EmV7Q97YCHQ0cy65S5sEZEYLFgAmzcXTH8LxJdc+gEjzOxVM3spXO2yurNIn1y6uvsagPB5/3D7QcDHVY6rCLeJiBSuZBJatiyIIcgp2dznUi9mNhs4oIZdk8LrdgKGAUOBx82sT2oWADM7Etjq7qvretkattU4s4CZXQpcCtCzZ886XkZEpBElEnDssdC+fdyRZC1jcjGz7sDZwAjgQIIO/dXAs0DS3dOu7+LuaRcaMLMJwLQwmSwys0qgC8GINMJrpqu1AKw1s27uviZs9ko1q1UAPaoc1x34NE18k4HJAGVlZZraRkTy04cfwhtvwC9+EXckdZLpPpcHCYYbfwP8O0EfyGXAbOAkYIGZHVvP604Hjguv04/gzv/Pw/fNgDOAqRnOnwmMD1+PB2ZU2X6BBYYBG1PNZyIiBSmZDJ4LqL8FMtdc7kjTLLUamGZmrYD6tic9ADxgZqsJktf4KhNjHgtUuPv7VU8ws/uAe919CXAbQVPaL4CPCJIRQIJgGPLfCYYiX1TP+ERE8kMyCb16waGHxh1JnWRKLieZ2ZfuXlHTTnf/huCPeJ2F556XZt9cgr6Y6tsvqfJ6A/Cdnq0wQf2qPjGJiOSd7dthzhy44IKCGYKckmm02EHAK2Y2z8wmmFmXxgpKRESA+fNhy5aCuSu/qrTJxd2vJmj2uongTvmVZpY0swvMbO/GClBEpMlKJKBVKxg1Ku5I6izjfS7hzYgvufsEglFYdwJXA2sbIzgRkSYtmYSRI6Fdu7gjqbOsbqIMp2L5DXA3QQf8/40yKBGRJu+DD+CttwpulFhK2g59M+tLcL/JOGAXwdDgE6qP4hIRkQikhiAXYH8LZB4tNovgRsaz3H1VI8UjIiIQJJc+faBv37gjqZdMU+73acxAREQktG1bMAT54osLbghySlwTV4qISDrz5sHXXxdskxgouYiI5J9EAvbaKxgpVqCyHS3WxswOiToYEREh6G8ZNQrato07knqrNbmY2Y+A5cBz4ftSM5sZdWAiIk3Se+/BO+8U7BDklGxqLr8GjgC+BAiXHe4VXUgiIk1YgQ9BTskmuex0942RRyIiIkF/y/e/HzwKWDbJZbWZnQM0N7O+ZvZfwMsRxyUi0vR8/TW8+GLB11ogu+RyBdAf2E5wU+Um4KoogxIRaZLmzg3ucSnw/haoZZljAHffSrDu/aTowxERacKSSWjTBn74w7gjabBak4uZPQNUX2N+I7AE+JO7b4siMBGRJieRCIYgt2kTdyQNlk2z2PvAZuDP4WMTwZT7/cL3IiLSUO++GwxDLoL+Fsii5gIMdvdjq7x/xszmufuxZvZ6VIGJiDQpiUTwXAT9LZBdzWU/M+uZehO+Ti15/E0kUYmINDXJJBxySDATchHIpubyf4AFZvYeYEBv4DIzawdMiTI4EZEmYevWYKTYhAlxR5Iz2YwWS4QLhx1KkFzeqtKJf2eUwYmINAkvvgjbtxdNfwtkV3MB6AscArQGSswMd384urBERJqQRCKYpPLYY2s/tkBkMxT5FmAk8AMgAZwMLACUXEREGso9SC6jRwfT7BeJbDr0fw6MBj5z94uAQUDx/AREROL09ttQXl40o8RSskkuX7t7JbDTzPYB1gHFMZxBRCRuqVmQiyy5ZNPnssTMOhLcMLmU4IbKRZFGJSLSVCQScNhh0KtX3JHkVDajxS4LX95rZs8B+7j7ymjDEhFpAjZvhnnz4PLL444k57JZiXJO6rW7l7v7yqrbRESknv72N/jmm6IagpyStuZiZq2BtkAXM+tEcI8LwD7AgY0Qm4hIcUsmoV07OOaYuCPJuUzNYr8kWLflQGBZle2bgLujDEpEpOi5B8llzJiiGoKckrZZzN3/4O69gWvcvXeVxyB3v6sRYxQRKT5vvgkfflh0o8RSMjWL/Sx8+UmV17u5+7TIohIRKXZFOgQ5JVOz2I8y7HNAyUVEpL4SCejfH3r2rP3YApQ2uYR340fGzK4ALgd2As+6+0QzOxe4tsphJcAQd19e7dx9gceAXkA5cKa7fxGef1142GZggruviLIcIiJ19tVXMH8+XHVV3JFEJpuhyB3M7HdmtiR83GFmHRpyUTMbBYwFSty9P3A7gLs/4u6l7l4KnA+UV08soeuBOe7eF5gTvgf4APihu5cA/w+Y3JA4RUQiMWcO7NhRtE1ikN30Lw8AXwFnho9NwIMNvO4E4DZ33w7g7utqOGYc8Gia88fy7VoyU4CfhJ/zsrt/EW5fCHRvYJwiIrmXTMLee8Pw4XFHEplsksv33P0Wd38/fPwrDZ9brB8wwsxeNbOXzGxoDcecRfrk0tXd1wCEz/vXcMwvgGQD4xQRya3ULMhjxkCrVnFHE5ls5hb72syOcfcFAGY2HPi6tpPMbDZwQA27JoXX7QQMA4YCj5tZH3f38Nwjga3uvjq7Ynzn2qMIkkvaO5PM7FLgUoCeRdqhJiJ56PXXoaICbrkl7kgilU1ymQBMqdLP8gUwvraT3H1Mun1mNgGYFiaTRWZWCXQB1oeHnE36WgvAWjPr5u5rzKwbwUzNqc8uAe4DTnb3DRnim0zYJ1NWVua1lUdEJCcSieD5pJPijSNi2SSXVe4+KJxuH3fflIPrTgeOA+aaWT+gFfA5gJk1A84AMi3JNpMgwd0WPs8Iz+1JMET6fHd/JwdxiojkVjIJJSXQvbi7hLPpc/nAzCYTNF99laPrPgD0MbPVwFRgfKpJjCCpVLj7+1VPMLP7zKwsfHsbcLyZvQscH74HuBnoDPzRzJab2ZIcxSsi0nCbNsGCBUU9SizFvv2bnuYAszYEN1SeDQwB/geYmuqDKQZlZWW+ZInykIhEbNo0OP10eOklODZT40xhMLOl7l5W075aay7u/rW7P+7uPwMGE8yK/FKOYxQRKX6JBOyzDxx1VNyRRC6bZjHM7Idm9keC2ZFbE9zvIiIi2UrNgnzCCdCyZdzRRK7WDn0z+wBYDjwOXOvuWyKPSkSk2KxcCZ9+2iT6WyC70WKDcjRCTESk6UrNglzkQ5BT0jaLmdmNZrZvusRiZseZ2WnRhSYiUkQSCSgthQObxkK+mWouq4BnzGwbQV/LeoL+lr5AKTAb+LfIIxQRKXRffgkvvwwTJ8YdSaPJNOX+DGCGmfUFhgPdCCat/AtwqbvXOgWMiIgAL7wAu3bBKafEHUmjyabP5WN3f6jqBjPrQhbzi4mICEF/S8eOMGxY3JE0mmyGIi82s90/ETM7HXg5upBERIpIZeW3Q5BbZPP/+eKQTUnPAR4ws7nAgQTTqxwXZVAiIkVjxQr47LMmMwQ5pdbk4u6rzOxW4L8J5hY71t0rIo9MRKQYNJFZkKvL5ibK+4HvEaxn349gBNld7n531MGJiBS8ZBKGDIEDalreqnhl0+eyGhjl7h+4+yyCBb6GRBuWiEgR+OILeOWVJjVKLCWbZrHfV3u/kWCVRxERyeT554MO/SbW3wLZNYv1BX4L/IDgJkoA3L1PhHGJiBS+ZBL23ReOPDLuSBpdNs1iDwL3ADuBUcDDBJ37IiKSTtUhyM2bxx1No8smubRx9zkEC4t96O6/RkORRUQye+01WLeuSfa3QHb3uWwL17V/18wuBz4B9o82LBGRApcagnziifHGEZNsai5XAW2BfwEOB84HxkcZlIhIwUsmYehQ2L9p/l88m9Fii8OXm4GLog1HRKQIbNgACxfCTTfFHUlsshktVgZMAg6uery7l0QYl4hI4Xr++WBZ4yba3wLZ9bk8AlxLsL5LZbThiIgUgUQCOneGsrK4I4lNNsllvbvPjDwSEZFiUFkJzz0XzCXWBIcgp2STXG4xs/uAOcD21EZ3nxZZVCIihWrJEvj88yZ5V35V2SSXi4BDgZZ82yzmgJKLiEh1ySSYNdkhyCnZJJdB7j4w8khERIpBIgFHHAFdusQdSayyuc9loZn9IPJIREQK3fr1sHhxk28Sg+xqLscA483sA4I+FwNcQ5FFRKqZNavJD0FOySa5NK3l00RE6iuZhP32g8MPjzuS2GVzh/6HjRGIiEhB27UrGIJ86qnQLJseh+Kmn4CISC4sXgz/+If6W0JKLiIiuZBIBDWWE06IO5K8oOQiIpILyWSw4mTnznFHkheUXEREGmrt2uDOfI0S2y225GJmV5jZ22b2upn9R7jtXDNbXuVRaWalNZy7r5m9YGbvhs+dqu0fama7zOznjVUeEWnCZs0KntXfslssycXMRgFjgRJ37w/cDuDuj7h7qbuXEixKVu7uy2v4iOuBOe7el2DOs+urfHZz4N+BWREXQ0QkkEhA164weHDckeSNuGouE4Db3H07gLuvq+GYccCjac4fC0wJX08BflJl3xXAU0BNnykikls7dwbrt5x0koYgVxHXT6IfMMLMXjWzl8xsaA3HnEX65NLV3dcAhM/7A5jZQcBPgXsjiFlE5LtefRW++EL9LdVkc4d+vZjZbOCAGnZNCq/bCRgGDAUeN7M+7u7huUcCW919dR0veydwnbvvMrPa4rsUuBSgZ8+edbyMiEgomQxqLMcfH3ckeSWy5OLuY9LtM7MJwLQwmSwys0qgC7A+PORs0tdaANaaWTd3X2Nm3fi2CawMmBomli7AKWa2092n1xDfZGAyQFlZmdetdCIioWQSjj4aOnWq/dgmJK5msenAcQBm1g9oBXwevm8GnAFMzXD+TGB8+Ho8MAPA3Xu7ey937wU8CVxWU2IREcmJzz6DZcs0SqwGcSWXB4A+ZraaIImMTzWJAccCFe7+ftUTzOw+M0stSH0bcLyZvQscH74XEWlczz0XPKu/5Tvs27/pTVdZWZkvWbIk7jBEpNCceSYsWACffBKsPtnEmNlSdy+raZ/GzYmI1EdqCPLJJzfJxFIbJRcRkfp45RXYuFH9LWkouYiI1EcyCc2bw5i0A2ObNCUXEZH6SCRg+HDo2DHuSPKSkouISF198gmsWKEmsQyUXERE6kpDkGul5CIiUlfJJBx0EAwcGHckeUvJRUSkLnbsgBde0BDkWii5iIjUxcsvw6ZN6m+phZKLiEi2NmyAP/0JWrTQEORaRDYrsohIUdixI+hjmTIFnnkmeH/RRbDPPnFHlteUXEREqnOH5cuDhPLXv8L69bDffvCrX8H48VBaGneEeU/JRUQk5bPP4JFHgqSyahW0bAk/+lGQUE4+OXgvWVFyEZGmbds2mDkzSCizZsGuXXDEEXD33XDWWdC5c9wRFiQlFxFpetxh4cIgoTz2GHz5JRx4IFxzTVBLOeywuCMseEouItJ0fPQR/Pd/w8MPwzvvQJs28NOfBgll9OhgIkrJCSUXESluW7bAtGnw0EPw4otBrWXECJg4Ec44Q6O+IqLkIiLFp7IS5s0Lmr2efBI2b4beveHmm+GCC6BPn7gjLHpKLiJSPP7+96DJ6+GH4cMPYe+9g6WIx4+HY46BZrpvvLEouYhIYdu4ER5/PKil/O//BvN9jRkDt94a9Ke0bRt3hE2SkouIFJ5du4LJI6dMgenTg+HEhx4Kv/0tnHcedO8ed4RNnpKLiBSO118PEspf/gJr1kCnTnDxxUGz19ChmqU4jyi5iEh++/xzePTRIKksXRoMFz7llCChnHYa7LVX3BFKDZRcRCT/fPNNsEb9lCnw7LPBZJGlpfD738M558D++8cdodRCyUVE8oM7LFsWJJRHHw1qLF27whVXBLWUkpK4I5Q6UHIRkXitWRP0oUyZEvSptGoFY8cGCeXEE4O1U6Tg6FsTkcb39dcwY0aQUJ5/PrjpcdgwuOeeYLLITp3ijlAaSMlFRBqHe7BE8JQpwX0pGzdCjx5w/fXBXfOHHBJ3hJJDSi4iUrPKyuD+ka1bg5pGTc912ffOO/Dee8FNjaefHjR7jRqlu+aLlJKLSCFxh+3bG/ZHPttjtm2rX4wtWgQJpE2b4Dn1ul8/uPHGILHsvXdufy6Sd8zd444hdnu3bGR2qsIAAAgMSURBVOmH13dBIP38Gi6qn2FNn1t9W13fx/EZlZXBHemVlcGjvpo1C+4RadZsz9e53qYbGQtKaWkpd955Z73ONbOl7l5W0z7VXCD4xd24Me4oilOcf2hqunb1bXV9X9O26s06ufjMqu9zkQD0B18amZILcMiQIcxdsiTuMEREioZ60kREJOdiSy5mdoWZvW1mr5vZf4TbzjWz5VUelWZWWsO5+5rZC2b2bvjcqcq+keG5r5vZS41ZJhERCcSSXMxsFDAWKHH3/sDtAO7+iLuXunspcD5Q7u7La/iI64E57t4XmBO+x8w6An8Efhx+7hnRl0ZERKqLq+YyAbjN3bcDuPu6Go4ZBzya5vyxwJTw9RTgJ+Hrc4Bp7v5Rhs8VEZGIxZVc+gEjzOxVM3vJzIbWcMxZpE8uXd19DUD4nJoitR/QyczmmtlSM7sg55GLiEitIhstZmazgQNq2DUpvG4nYBgwFHjczPp4eNONmR0JbHX31XW8bAvgcGA00AZ4xcwWuvs7NcR3KXApQM+ePet4GRERySSy5OLuY9LtM7MJBM1XDiwys0qgC7A+PORs0tdaANaaWTd3X2Nm3YBU81cF8Lm7bwG2mNk8YBDwneTi7pOByQBlZWW6E1JEJIfiahabDhwHYGb9gFbA5+H7ZgQd8VMznD8TGB++Hg/MCF/PIGhua2FmbYEjgTdzHr2IiGQUy/QvZtYKeAAoBb4BrnH3v4X7RhJ09g+rds59wL3uvsTMOgOPAz2Bj4Az3P0f4XHXAhcBlcB97l7rvAZmth74sMqmDsDGNO9Tr1PPXQgTYz1Vv1Zdjqlpezaxp3vdkLI0pBzp9hViWepajurvq//7gsIpS5TfSaY4szkmn8qSD78rufr3dbC771fjHnfXo9oDmJzufep1leclubxWXY6paXs2sWcoU73L0pByFFNZ6lqO2v59FVJZovxOiqks+fC7kqt/X5keukO/Zs9keP9MmmNyda26HFPT9mxiz/S6vhpSjnT7CrEsdS1H9ff695VesZQlH35XcvWdpKVZkRvIzJZ4mllBC43Kkp+KpSzFUg5QWbKhmkvDTY47gBxSWfJTsZSlWMoBKkutVHMREZGcU81FRERyTslFRERyTslFRERyTsklQmbWzMxuNbP/MrPxtZ+Rv8J1cuab2b3hja4FzczahZObnhZ3LPVlZoeF38eT4ZRKBcvMfmJmfzazGWZ2QtzxNISZ9TGz+83sybhjqavw92JK+F2c25DPUnJJw8weMLN1Zra62vaTwkXO/m5m19fyMWOBg4AdBPOexSJHZXFgM9Cawi8LwHUEszzEIhflcPc33f2fgTOB2IbF5qgs0939n4ALCWZEj0WOyvK+u/8i2kizV8cy/Qx4MvwuftygC0dxZ2YxPIBjgSHA6irbmgPvAX0I5kNbAfwAGAj8T7XH/gSLmP0yPPfJAi9Ls/C8rsAjBV6WMQSTo14InFao5QjP+THwMnBOIX8nVc67AxhSJGWJ7Xe+AWW6ASgNj/lrQ64b2azIhc7d55lZr2qbjwD+7u7vA5jZVGCsu/8W+E7ziplVEMydBrArumgzy0VZqvgC2CuKOLORo+9lFNCO4JfpazNLuHtlpIFXk6vvxN1nAjPN7Fngr9FFnF6OvhMDbgOS7r4s2ojTy/HvSl6oS5kIWiW6A8tpYMuWkkvdHAR8XOV9BcHMy+lMA/7LzEYA86IMrB7qVBYz+xlwItARuCva0OqsTmVx90kAZnYhwRINjZpYMqjrdzKSoBljLyARaWR1V9fflSsIapQdzOz77n5vlMHVUV2/l87ArcBgM7shTEL5Jl2Z/hO4y8xOpYFTxCi51I3VsC3tXajuvhXIm7bXaupalmkEyTIf1aksuw9wfyj3oTRIXb+TucDcqIJpoLqW5T8J/rDlo7qWZQPwz9GFkxM1lsmDtbAuysUF1KFfNxVAjyrvuwOfxhRLQ6ks+adYygEqS76LvExKLnWzGOhrZr3DNWnOJli4rBCpLPmnWMoBKku+i75McY9kyNcHwTLLa/h2GPEvwu2nECyb/B4wKe44VZbCLEuxlENlyf9HXGXSxJUiIpJzahYTEZGcU3IREZGcU3IREZGcU3IREZGcU3IREZGcU3IREZGcU3IRyWNm9uPUdOhm9mszuybumESyobnFRPKYh7Mexx2HSF2p5iISITM7z8wWmdlyM/uTmTU3s81mdoeZLTOzOWa2X3jsv5jZG2a2MpwCHTO70My+Mwu1mZWa2cLw2KfNrFO4fa6Z/Xt4zXfCGblFGp2Si0hEzOwwglUVh7t7KcGaPucSrCWzzN2HAC8Bt4SnXA8MdvcSap9V92HguvDYVVU+A6CFux8BXFVtu0ijUbOYSHRGA4cDi4O1sGgDrAMqgcfCY/7Ct0sZrAQeMbPpwPR0H2pmHYCO7v5SuGkK8ESVQ1KftxTo1eBSiNSDai4i0TFgiruXho9D3P3XNRyXmuDvVOBugoS01Mzq+5+/7eHzLvQfSImJkotIdOYAPzez/QHMbF8zO5jg9+7n4THnAAvMrBnQw91fBCYSrPjZvqYPdfeNwBdV+lPOJ2heE8kb+l+NSETc/Q0zuxF4PkweO4BfAVuA/ma2FNhI0C/THPhL2ORlwO/d/cuwOa0m44F7zawt8D45Wj1QJFc05b5IIzOzze5eY61EpFioWUxERHJONRcREck51VxERCTnlFxERCTnlFxERCTnlFxERCTnlFxERCTnlFxERCTn/j+8HgDYZDVqqAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = pl.figure(2)\n",
    "pl.semilogx(epsilon_array,maxVs,color='red',label='maxV')\n",
    "#show singular solution as horizontal line\n",
    "pl.semilogx(epsilon_array,np.ones(len(epsilon_array))*maxVs[-1],color='black',label='singular')\n",
    "pl.xlabel('epsilon')\n",
    "pl.ylabel('max(voltage V) (mV)')\n",
    "pl.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The maximum of the voltage traces show that the non-singular case nicely converges to the singular one and no numeric instabilities occur. "
   ]
  },
 {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----------------------------\n",
    "### License\n",
    "\n",
    "This file is part of NEST. Copyright (C) 2004 The NEST Initiative\n",
    "\n",
    "NEST is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 2 of the License, or (at your option) any later version.\n",
    "\n",
    "NEST is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details."
    ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}