Pure Python and Cython code of an algorithm computing all pairwise correlograms between all neurons, given a set of spike trains.

correlograms.ipynb

{
 "metadata": {
  "name": "correlograms"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Computing all pairwise correlograms within a pool of neurons\n",
      "\n",
      "Pure Python and Cython code of an algorithm computing all pairwise correlograms between all neurons, given a set of spike trains."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Imports"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from itertools import product\n",
      "import numpy as np"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext cythonmagic"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Function definitions"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Pure Python"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def compute_correlograms(spiketimes, neurons, neurons_to_update=None,\n",
      "    ncorrbins=100, corrbin=.001):\n",
      "    \"\"\"Compute all pairwise cross-correlograms (and auto-correlograms)\n",
      "    between all neurons within a pool of neurons.\n",
      "    \n",
      "    Arguments:\n",
      "\n",
      "      * spiketimes: a sorted 1D array with all spiketimes,\n",
      "      * neurons: the neuron index for each spike (same shape as spiketimes),\n",
      "      * neurons_to_update=None: a 1D array with the neurons for which the\n",
      "        correlograms need to be computed. The correlograms between these neurons\n",
      "        and all neurons are computed.\n",
      "      * ncorrbins=100: the total number of bins in the correlograms. Need to be\n",
      "        an even number.\n",
      "      * corrbin=.001: the bin size, in seconds.\n",
      "\n",
      "    Returns:\n",
      "\n",
      "      * correlograms: a dictionary `(neuron0, neuron1): correlogram` where\n",
      "        correlogram is a 1D ncorrbins-long array with spike count values in each\n",
      "        bin.\n",
      "    \n",
      "    \"\"\"\n",
      "    # Ensure ncorrbins is an even number.\n",
      "    assert ncorrbins % 2 == 0\n",
      "    \n",
      "    # Compute the histogram corrbins.\n",
      "    n = ncorrbins // 2\n",
      "    halfwidth = corrbin * n\n",
      "    \n",
      "    # size of the histograms\n",
      "    nspikes = len(spiketimes)\n",
      "\n",
      "    # unique neurons\n",
      "    neurons_unique = np.unique(neurons)\n",
      "    nneurons = len(neurons_unique)\n",
      "    neuron_max = neurons_unique[-1]\n",
      "    \n",
      "    # neurons to update\n",
      "    if neurons_to_update is None:\n",
      "        neurons_to_update = neurons_unique\n",
      "    neurons_mask = np.zeros(neuron_max + 1, dtype=np.bool)\n",
      "    neurons_mask[neurons_to_update] = True\n",
      "    \n",
      "    # initialize the correlograms\n",
      "    correlograms = np.zeros(\n",
      "        ((neuron_max + 1) ** 2, ncorrbins), dtype=np.int32)\n",
      "\n",
      "    # loop through all spikes, across all neurons, all sorted\n",
      "    for i in range(nspikes):\n",
      "        t0, neuron0 = spiketimes[i], neurons[i]\n",
      "        # pass neurons that do not need to be processed\n",
      "        if neurons_mask[neuron0]:\n",
      "            # i, t0, c0: current spike index, spike time, and neuron\n",
      "            # boundaries of the second loop\n",
      "            t0min, t0max = t0 - halfwidth, t0 + halfwidth\n",
      "            j = i + 1\n",
      "            # go forward in time up to the correlogram half-width\n",
      "            while j < nspikes:\n",
      "                t1, neuron1 = spiketimes[j], neurons[j]\n",
      "                if t1 <= t0max:\n",
      "                    d = t1 - t0\n",
      "                    k = int(d / corrbin) + n\n",
      "                    ind = (neuron_max + 1) * neuron0 + neuron1\n",
      "                    correlograms[ind, k] += 1\n",
      "                else:\n",
      "                    break\n",
      "                j += 1\n",
      "            j = i - 1\n",
      "            # go backward in time up to the correlogram half-width\n",
      "            while j >= 0:\n",
      "                t1, neuron1 = spiketimes[j], neurons[j]\n",
      "                if t0min <= t1:\n",
      "                    d = t1 - t0\n",
      "                    k = int(d / corrbin) + n - 1\n",
      "                    #ind = pairs[(neuron0, neuron1)]\n",
      "                    ind = (neuron_max + 1) * neuron0 + neuron1\n",
      "                    correlograms[ind, k] += 1\n",
      "                else:\n",
      "                    break\n",
      "                j -= 1\n",
      "    return {(neuron0, neuron1): correlograms[(neuron_max + 1) * neuron0 + neuron1,:] for neuron0 in neurons_to_update for neuron1 in neurons_unique}"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Cython version"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "\n",
      "import numpy as np\n",
      "cimport numpy as np\n",
      "DTYPE = np.double\n",
      "ctypedef np.double_t DTYPE_t\n",
      "DTYPEI = np.int\n",
      "ctypedef np.int_t DTYPEI_t\n",
      "\n",
      "def compute_correlograms_cython(\n",
      "     np.ndarray[DTYPE_t, ndim=1] spiketimes,\n",
      "     np.ndarray[DTYPEI_t, ndim=1] neurons,\n",
      "     np.ndarray[DTYPEI_t, ndim=1] neurons_to_update=None,\n",
      "     int ncorrbins=100,\n",
      "     double corrbin=.001):\n",
      "    \n",
      "    # Ensure ncorrbins is an even number.\n",
      "    assert ncorrbins % 2 == 0\n",
      "    \n",
      "    # Compute the histogram corrbins.\n",
      "    cdef int n = ncorrbins // 2\n",
      "    cdef double halfwidth = corrbin * n\n",
      "    \n",
      "    # size of the histograms\n",
      "    cdef int nspikes = len(spiketimes)\n",
      "    \n",
      "    cdef int i, j, neuron0, neuron1, k, ind\n",
      "    cdef float t0, t1, t0min, t0max, d\n",
      "\n",
      "    # Unique neurons\n",
      "    cdef np.ndarray[DTYPEI_t, ndim=1] neurons_unique = np.unique(neurons)\n",
      "    cdef int nneurons = len(neurons_unique)\n",
      "    cdef int neuron_max = neurons_unique[-1]\n",
      "    \n",
      "    # neurons to update\n",
      "    if neurons_to_update is None:\n",
      "        neurons_to_update = neurons_unique\n",
      "    cdef np.ndarray[DTYPEI_t, ndim=1] neurons_mask = np.zeros(neuron_max + 1, dtype=DTYPEI)\n",
      "    neurons_mask[neurons_to_update] = 1\n",
      "    \n",
      "    # initialize the correlograms\n",
      "    cdef np.ndarray[DTYPEI_t, ndim=2] correlograms = np.zeros(\n",
      "        ((neuron_max + 1) ** 2, ncorrbins), dtype=DTYPEI)\n",
      "\n",
      "    # loop through all spikes, across all neurons, all sorted\n",
      "    for i in xrange(nspikes):\n",
      "        t0, neuron0 = spiketimes[i], neurons[i]\n",
      "        # pass neurons that do not need to be processed\n",
      "        if neurons_mask[neuron0]:\n",
      "            # i, t0, c0: current spike index, spike time, and neuron\n",
      "            # boundaries of the second loop\n",
      "            t0min, t0max = t0 - halfwidth, t0 + halfwidth\n",
      "            j = i + 1\n",
      "            # go forward in time up to the correlogram half-width\n",
      "            while j < nspikes:\n",
      "                t1, neuron1 = spiketimes[j], neurons[j]\n",
      "                if t1 < t0max:\n",
      "                    d = t1 - t0\n",
      "                    k = int(d / corrbin) + n\n",
      "                    ind = (neuron_max + 1) * neuron0 + neuron1\n",
      "                    correlograms[ind, k] += 1\n",
      "                else:\n",
      "                    break\n",
      "                j += 1\n",
      "            j = i - 1\n",
      "            # go backward in time up to the correlogram half-width\n",
      "            while j >= 0:\n",
      "                t1, neuron1 = spiketimes[j], neurons[j]\n",
      "                if t0min < t1:\n",
      "                    d = t1 - t0\n",
      "                    k = int(d / corrbin) + n - 1\n",
      "                    ind = (neuron_max + 1) * neuron0 + neuron1\n",
      "                    correlograms[ind, k] += 1\n",
      "                else:\n",
      "                    break\n",
      "                j -= 1\n",
      "    return {(neuron0, neuron1): correlograms[(neuron_max + 1) * neuron0 + neuron1,:] for neuron0 in neurons_to_update for neuron1 in neurons_unique}"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Tests\n",
      "\n",
      "Computing all pairwise correlograms within a pool of 100 neurons, with 100 time steps (+/- 50 ms, 1 ms bin), for a total of 100000 spikes."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nspikes = 100000\n",
      "nneurons = 100\n",
      "ncorrbins = 100\n",
      "corrbin = .001"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "spiketimes = np.cumsum(np.random.exponential(scale=.005, size=nspikes))\n",
      "neurons = np.random.randint(low=0, high=nneurons - 1, size=nspikes)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plotting a few correlograms."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "correlograms = compute_correlograms_cython(spiketimes, neurons,\n",
      "                                     ncorrbins=ncorrbins, corrbin=corrbin)\n",
      "for i in xrange(4):\n",
      "    for j in xrange(4):\n",
      "        subplot(4, 4, i * 4 + j + 1);\n",
      "        bar(arange(100) - 50, correlograms[i,j], width=1, ec='none');\n",
      "        xlim(-50,50);\n",
      "        grid();"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAWYAAAD9CAYAAACP8N0iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztfX1wFVWa/tMKjLIKODV8lcEJBAgEITeCkylLikTEzDpQ\nyMJUCQgJomNpOSPU6jKONZUw8zMJJZpkdWq2amtxmLVK/9jSAixAYL13ZEUHBxJnB3Rcx6QIED6k\n+IhBCEne3x+dvrc/zld33743l5yn6lbSfc77cZ5+73u7T7992iAigoaGhobGgMEN2XZAQ0NDQ8MJ\nnZg1NDQ0Bhh0YtbQ0NAYYNCJWUNDQ2OAQSdmDQ0NjQEGnZg1NDQ0BhiUEnN+fj5mzZqFkpIS/OAH\nP4jap0GFCxcuYNmyZZg+fTqKiorw8ccfZ9ul6waa22ih80J0GKLSyTAMJBIJfPe7343an0GHZ555\nBg8++CD+67/+Cz09Pejq6sq2S9cNNLfRQueF6KCUmAFAP4eSfly8eBH79+/H1q1bAQBDhgzByJEj\ns+zV9QHNbWag80I0UJrKMAwD999/P+bMmYN///d/j9qnQYPW1laMHj0aa9aswV133YXHH38cly9f\nzrZb1wU0t9FD54UIQQo4efIkERGdOXOGiouL6YMPPki2AdAfhQ8Ln3zyCQ0ZMoQOHjxIRETPPPMM\n/epXv3L0ybbfufLR3GaOW5W8oPkNx63SGfP48eMBAKNHj8aSJUtw8OBBRzsRJT8AObarq6sd+93t\nvI8lF+QTVBZQk2ONxW7TPUYe8vLykJeXh7vvvhsAsGzZMhw+fNjTL2qOVI9Jto6LTG4gc+sv/qLj\nKKhsmLwQBb9iX7MTf1FwK03Mly9fRmdnJwCgq6sLe/bswcyZM2ViGgoYN24cJkyYgC+++AIAsG/f\nPsyYMSPLXl0f0NxGC50XooX05t/p06exZMkSAEBPTw9WrlyJBx54QNlAW1tbIMeCyoWTzbzNV199\nFStXrkR3dzcKCgrw+uuvR24zjGwu2cw1bsPIZtrmYMoL2bApTcwTJ05ES0tLIOUAEIvFMioXTjbz\nNhcvXowRI0bgxhtvxPHjx31VDmSDo1yymWvchpHNtM3BlBeyYdMg2WSHTIFhOOZLDANgabT289oH\nAlR9k43Fvd/NkR0TJ07EoUOHuLWgItl0YSAfExXwOBoI3PrBQDwOYTjKNL8DkT8RRPwo3fzr7e1F\nSUkJFi1alFbHNEwMpORwvUFzGx10XogOSom5qakJRUVFMAzDt4FEIuFbJoxcONnM2wxTC5oNjnLJ\nZq5xG0Y2GzYHS17Ihk3pHPPx48exc+dOvPDCC3jllVcCGdHg48MPP8T48eNx9uxZLFiwANOmTcPc\nuXMdfaqqqpCfnw8AGDVqFGKxGMrKygCkDrx9u7wcIPK2GwYQjyccuq12wygDEVsfAJSXm+3WvKLI\nfrq3W1paPNsXLlwAIL65EgW3KtsWVPtb3ALesQaxH8ZfVW51XogW0jnmn/zkJ/jlL3+JS5cuYfPm\nzdixY4dTgWGgsrIyGdwbN45CPJ4KbsNIIB5PBZ+1zQuWsO0q27LExWt3y/P6l5e3oLo6Fdxbt25V\nuqTeuHEjbrnlFvzzP/+zg1+/l+Oqc9/2/YB4fm4gz9+pcJQubqPAQL7/wuNIlhcsWXtuYP3wyb5r\nvG37SYY9NwDB9KnoD/tD2NjYiJaWFluu3MiPPxJgx44d9NRTTxERUTwep4ULF3r6uFW4NVrb7r88\nhG1XAU+HXx95/b3bbIVdXV106dIlIiL65ptv6J577qH33ntPSVbFPz/7M8F7VGBxFBW3UUA17rIB\nFkcqeYEn6+0T1C/2vnRxmIljIeJHOMd84MABbN++HRMnTsTy5cvx/vvvY/Xq1SIRBhI++/dLDYI5\n5tOnT2Pu3LmIxWIoLS3FwoULfdWCDpb5tsHCbRjZTNpMR17IhXFm06YwMdfW1qK9vR2tra146623\ncN999+EPf/hDIEMaXli1oIcOHcLQoUNx4MCBbLt03UBzGx10XogeynXMf/zjH/Hyyy9j+/btTgWS\nOmb3/JlsHi1suwrC1lrLxuSnjhkAXnnlFRw6dAidnZ1SflWg55hTSDe3USAX55gt8PKCiqzZJ9iY\nWXIqcRxGf7oRuo4ZAObNm8ckXyMcrLvbjz322IBIEtcTNLfRQ+eFaJCBd/4lgkkNgjlmAFi/fj1e\neukl3HCD/0MxWObbBgu3YWSz5W+mbQ4WbqV1zFeuXMG8efNw9epVdHd3Y/HixairqwtkTMOJd999\nF2PGjEFJSYnwABpGFaqr8wF4yxH5tall/bKpckVnO5Lt9v688kSrnVXHrFI3HaZEKUgdsyq3ojpm\nldJM1tgtqJZ2WtxadcyqpaXWdnl5GeJx/6WiqWOLZLtqHbPOCxFDpayjq6uLiIiuXbtGpaWltH//\nfm7Jhy6XUy+Xe/755ykvL4/y8/Np3LhxNHz4cFq1apVH1i6uMn5Wf14pkbU/DCd+y/P86pHLeQVV\nuQ3rTxBeeP38flf82nH3VRsfv5MoL8hk/figKnc9lcv5Mt/V1UVz5syhI0eOcJXrxKyemO1IJBLc\nOnGdmGVyYkERt2H9GcyJ2QIrL6jK6sTMhtLkW19fH2KxGMaOHYvy8nIUFRU52quqqlBTUwPDqAHQ\nCMNIJO+QAo2uS7yEYzuRMLdT/RMwjFQfq50lbxheedN2IwyD3W7Keu3bbVvtdnl7f7e806cEgEbU\n1NSgpqYGVVVVUIXfNQcGy3xbOuZAc4HbMLLZsCnLC4CZGwyjpj8/pHJD6nuVss36rrm/242NjY72\n1FQc+nU5v+v2XNTYmMpFZn5g27Ngyjq/27L+7na7v42NjclcWVNT4+HKAT8Z/sKFC1RaWkrxeJyZ\n9e1nX6lf/Xiyzf7X++vh1GG3werH+9+Uj3t+PVXOFNxyqmeRdn9Zv9o8mr/99lv6wQ9+QMXFxTR9\n+nT6xS9+wbDFP2P2w5HbL7e/qmNl2VQ9Y+b5K9Mjl/MKqnLLQzweD3zGbOdWBe44CnLGLOPIbcv+\nv0hWJUWw8oJd1h2D7u+Lqs8sOfd4WN971jhlZ9dBfBXFArs/X6Hv9Zh/85vf4Oabb8azzz7b/yuR\nqsVzn5QQqdcx29stWRbs8qz/WT6I5Ny6Wb6LfHX3Y/kvqle8fPkyhg8fjp6eHtx7773YvHkz7r33\nXocsQFLfWf65fXX7xRqzbKwyeyr7Vfz2Ax6/KtyKwt8v10Hk7f381vz7tePuqzY+tVpvd16wy7K+\nm0HimRfX7u+eW799266Ttd/d7sdX//Eeoo7566+/Tt6l/fbbb7F3716UlJSoW9cQYvjw4QCA7u5u\n9Pb2chd11/APzW100HkhWkgTc0dHB+67777kmgOLFi3C/PnzfZhIBHIs3LxiUNngNqOcp0u3zTCy\nuWQz17gNI5tpm2HzQq6MM1s2pXXMM2fOZL72XSM9uOGGG9DS0oKLFy+ioqICiUQiWXeaQhVqavL7\n/x+FRMJbx2zVvqYCwdKRgLnL2T8FZ3+rPeWDU768vIVZW8uTt/tn1drKlk215K32oOsxq3ArqmO2\nxs6qsbZv28deXg5bbXBKHmDXJZeXp+StOmb3sbNq0Hk14W77smUvvbGBZLvF7caNQGVlG3jQeSFa\nSOeY29vbsXr1apw5cwaGYeCnP/0pfv7zn6cU6DnmUHPMdvDm6VTmmFV5CTvHzFujIKhPIt3qc6xy\nfkVzoHy9/JhR8VllPPbjIuLF6uPXz6B9U7bZHMnygqlDzzGLEGqOeejQoWhoaMCRI0fw8ccf47e/\n/S0+++wzdesaXOh5uuiguY0WOi9EC2liHjduXPIV3LfccgumT5+OkydP+jCRCOTYYJhjztY8XTjZ\nzNscPNzmzjxo2LyQK+PMlk3pHLMdbW1taG5uRmlpaSBjGk7oebrooLnNHHReiADikukUOjs7afbs\n2fTOO+94iqQrKyupurqagGoCGsh6UMNsj7uK5uNkPXRi33a3m7M17P72gm9Tt1jebt/qz7Jvt8ey\nz2pn62ug6upqqq6upsrKSuLRfOzYMSorK6OioiKaMWMGNTU1efoAao9k8/qkeGHLpsYp1uv+y7Pt\n1yeRbtXoDMqvLPxFvKv4rDIeO/eivyJXVeLDb9+UbbFCXl6wZM34r+7PD87cYP/uWt8l50Mg3tzh\n57vPa7fk7frs/tjbZbmC5a/dP+eDXA0uPvjcKoV+d3c3PfDAA9TQ0OBVYFNu/5KrBhmrP+tj7+/+\nXybPkxPZ4+lQ0e39QrJp7ujooObmZiIyA3zq1Kl09OhRj6xOzGIE5VcnZnk/EUeivGCXZX1/7HZE\nscn6y9Lp1iVqd//Psynylf09Z//Psm3+5XMrnWMmIqxduxZFRUVYt25dgHPyRACZMHLZsZlL83Th\nZDNvM6hcGH71HLMYYfNCNsaZjftdkc0xf/jhh3jjjTcwa9as5F3turo6/OhHPwpkUIMN8TxdFQwj\nv///VB2zWQPrXC/ZWtTFMMr6+yf6S3/M7VRdMZLtJpy1uVbtrLuOGWhx1OZa+gyjDETe2ln7es5m\nja/TP7u8vb9dvqGhxeFPQ4NaHbMdPH6tOuaNGwHAuda1u46YVxdsjd3i3s6tWULl1Cdbj9m+zTp2\ndm4seTu3rP6sumt3HbPFtcVtVVUbeNB5IWJwz6UVYVchusznXUKoTBf4uUxXnW7g+Syyr6rbexkj\nplk2TyfigsevCo+q7W79Xh/l+3k8snwRybPtBOOXFbtu//20+Yld3hhZumX88mKR5y/PH/Z3K3iK\nsGRFx1A2PpnvsjhS+Z9nU+U7w5J1/8+ybf7lc5uBV0tpiHDt2jUsXboUjzzyCB566KFsu3PdQfOr\nkYuQJuZHH30UY8eOxcyZMwOaSGRYLjs2szFPl2scZXqOLwy/2amjz525/7B5Qc8xiyFNzGvWrMHu\n3bsDKdcQw5qni8fjKCkpQUlJieY6jdD8RgedFyIGd5LDhtbWVrrzzjuZbXYVojkv3tyObL5RNE+n\nOs/HkxPZE81XyXR755f4NK9Zs4bGjBkj5FfEhWxeUiar0u7W7/VRvp/HI8sXkTzbDptfFW5FYxDN\nR6rGnmjekSXP0i3jlxeLPH95/rC/W/zYFeUFu6zoGMrGJ/NdFkcq//NsqnxnWLLu/1m2zb98bn09\n+ceDYVQByO/fGgUgBusuv/3OsrVtypT190/0/xVvu6sM7G/6tfrb9fPubLOqFvzat7ez9bXAMC70\nb7dBhDVr1uBnP/sZVq9eLeyn4R+a2+zDyg3mm5S8ucFEWX9f57a3Ioj9Xebps1e9iHKDtdIh77tu\nVd2oyNtzk3f1wEZUVbXY+BCAm7JtkJ0xi8/S4tKzOPavnliO9wsmkhXL+PPVn6yYZr/8ptrizLMK\nGW92WX4771c+zvCRf/bh9CvO1S3ywXr6i2+Hz6/fqz0L7lc88c6OVGLXLePWY5dn2ZXx6+aWPVav\nD9b/rFcnpfoH49aS5R9nte+a2/fg+cT7fRHFII9bFTtOn72xoMJtWs6YNaJGFexXJIlEzFHXy7oi\n4V0BpOqKwW13yzv1tzjOEuy1sfb1lWW1uLKzJPv43LW9jY0tWL/+AqqrgY0b2xAG9vWY3Wtdu8/Y\n3L7yz8gseNvtY2PVFZeXtzD0se3zzhDdddXuq0v31Wc8DtsZXwuAC/392qCRJXBTtuIvIyA7Y1b/\nyH+5gvUNI5Oej5jmIGfMmeJZpsOtS6ZbxbZMp9c2n9+gZ8y88QXh0a1f9n86jqsf3eL2YNxa/IYd\nn8px8Dsu934V3X7iQOSzs43PrbQqY/ny5bjnnnvwxRdfYMKECXj99dej/J3Q0NDIAei8EDG4KVsR\ngJ5jlsuKaQ5+xhzMX1M+6FjjwrMC2XGRHwOvDtYcs7Mfn990zjGrxxPbX7X/gx8XFTusj5tfP7Er\ngj12vT6ojdOSC8IPz6Z7v4yfbMwxZ+DJv5YMy+WWzXBnHoODI3OO2T/CcBvUZr90FmSD2ww31sBW\nMyyXWzaliXn37t2YNm0apkyZgk2bNgUwcUHeJa1yuWWzsrISI0aMwB133IGnn34aa9asidxmONnM\n27QW1fGLMNwGtdkvnQXZ4DaDjjVcbsid+MuGTWFi7u3txdNPP43du3fj6NGjePPNN/V7vdIIzW90\n0NxGC81vtBAm5oMHD2Ly5MnIz8/H0KFD8fDDD2Pbtm0+TbQFdC2oXO7YDM+vf5vhZTNvU3VpTzvC\nchvEpk06C7LBbWaD31yKv2zYFNYxnzhxAhMmTEhu5+Xl4U9/+hOjp8HYZ8dWJWcMjxq+nLevXFYu\nI7aZbtmg/DrH4d/flHzQsW5N6nBzKuPYMOQ23Tq2bjVtuttFtlS5NQwxtyIb4rF6/VX9P+hxsXPL\n1+2Fm18V+I1dtg/qsRCEH55NvzFr9lHPRV6f/cUuIEnM7qBlwby5qBEEmt/ooLmNFprfaCGcyrj9\n9tvR3t6e3G5vb0deXl7kTg0WaH6jg+Y2Wmh+I4aoDvHatWs0adIkam1tpatXr1JxcbHnZaEawaH5\njQ6a22ih+Y0WwqmMIUOG4LXXXkNFRQV6e3uxdu1aTJ8+PVO/Gdc9NL/RQXMbLTS/ESPKrL9582Yy\nDIPOnTuX3FdbW0uTJ0+mwsJCeu+99zwyzz77LE2bNo1mzZpFS5YsoQsXLijL7tq1iwoLC2ny5MlU\nX18v9O3YsWNUVlZGRUVFNGPGDGpqaiIionPnztH9999PU6ZMoQULFtD58+eZ8j09PRSLxWjhwoW+\n5NKFTHNLpM5vWG6Jco9fza06dF6QQ5qYz58/T0uXLqVp06bR9OnT6aOPPlJSfOzYMaqoqKD8/Pzk\nAThy5AgVFxdTd3c3tba2UkFBAfX29jrk9uzZk9y3YcMG2rBhg5JsT08PFRQUUGtrK3V3d0svrTo6\nOqi5uZmIzJd1Tp06lY4ePUrPPfccbdq0iYiI6uvrk/bdePnll2nFihW0aNEiIiJlOTe+//3v08yZ\nMykWi9Hdd9+tJJNpbon88RuWW6L08BuEW6Jg/A42bj///HOKxWLJz4gRI5JJTASdF9TygjQxr169\nmv7jP/6DiMx5JfsvlQjLli2jTz/91HEAamtrHb9YFRUVwkT/9ttv08qVK5VkDxw4QBUVFcnturo6\nqqurU/KViGjx4sW0d+9eKiwspFOnThGReZAKCws9fdvb22n+/Pn0/vvvJ38ZVeRYsPOjikxzSxSO\nXz/cEqWP3yDcEoXndzBwa0dvby+NGzeOjh07Ju2r84Iat8KqjIsXL2L//v149NFHAZjzSiNHjpRO\nj2zbtg15eXmYNWuWY//Jkycdd27z8vJw4sQJrp4tW7bgwQcfVJJl1VWKdNvR1taG5uZmlJaW4vTp\n0xg7diwAYOzYsTh9+rSn//r16/HSSy/hhhtS9KnI8UA+yoqywS0QnF+/3ALp5dcPt0B6+B0s3FrY\nt28fCgoKHGNgQecFdW6FN/9aW1sxevRorFmzBp9++ilmz56NpqYmDB8+HAsWLMCpU6c8Mi+++CLq\n6uqwZ8+e5D73l8OSPXHiBPbu3Yua/ves1NbWYtGiRUk9w4YNw4oVK7j+2WspVeoqWfjmm2+wdOlS\nNDU14dZbb/Xod+t99913MWbMGJSUlHDfgMuS48EwDNx///248cYb8cQTT+Dxxx8fcNyytlXgl1sg\nvfyyuAUQmN/6+nrU19cDgIPfwcitHW+99VZyvAMtdnM1LwinMj755BMaMmQIHTx4kIiInnnmGfrV\nr37l6ANAfxQ+PJw8eZKIiM6cOUPFxcX0wQcfaH7TxK/mNjpuLVy9epW+973v0ZkzZzxt2fY7Fz48\nCKcy8vLykJeXh7vvvhsAsGzZMhw+fNjTj8y5auanurpa2J5uOT+yFj9R2xRh/PjxAIDRo0djyZIl\nOHjwoDK/dpv2cbDG6R135o6LZd9uU+SvX5tRcxuEIz/js/e1OGLJi3TKOUp/7ALArl27MHv2bIwe\nPZrZzrMfdV5gjXeg5SIRhIl53LhxmDBhAr744gsA5lzSjBkzpAdLQw2XL19GZ2cnAKCrqwt79uzB\nzJkzs+zV9QHNbWbw5ptvYvny5dl247qD9GWsr776KlauXInu7m4UFBT4foVM0FW6wqzulSs2T58+\njSVLlgAAenp6sHLlSjzwwAOR2rRJB5PKEZvZ5DZXOAore+LECbzzzjs4fPgwNm3ahC1btuCHP/xh\npDYHQ14AFBJzcXExPvnkk0DKASAWi2VULpdsTpw4MdTbI8L4C2Seo0zazCa3sVisf8W2QNKBbQa2\nGFD2l7/8JX73u9/h0UcfRU9PD7q6uiK3ORjyAgAYpDCRlJ+fjxEjRuDGG2/E0KFDHXN1hmEozUUN\nRFg3SKN2PwxHqrKGwR+Huy1T43bbt/sh8te//mD8Rhm7fsbH4oUlH4azoLI8ji5evIiSkhJ89dVX\nSrLpPN4qyLS9IBDFn9I7/wzDQCKRQHNzM/MGikY49Pb2oqSkJFkSpJE+aG6jgb2U9q677sLjjz+O\ny5cvZ9ut6wbKL2MNembBq+mLSi7XbAJAU1MTioqKfNdchrEJBJPNNZvZ4DbXOAoi29PTg8OHD+Op\np57C4cOH8Q//8A/JGm87qqqq+uuRa9DY2Ji0lUgkkh+7H7LtxsZGX/3t23b7fuTdPrvbrdCS+dvY\n2Jjkw6rR5kE6xwzwC/UtVFVVIT8/HwAwatQoxGIxlJWVAUi9gdfathyXbdsH51e+paUF5eVlIOL3\nLy8vsywgkXC2l5cDRF79hgHE43J/W1paki+4lE3+Hz9+HDt37sQLL7yAV155RdhXwx80t9GBVUrL\nSsy///3vAQAbNwLr1jnbrO9QprbtecmPvDuHBLW/zkXAxo0bwQUpQFSor6gi45C5Zc5AsfvxZIMO\nVcTRsmXL6PDhw5RIJJLP16vKqvrmbuONOypYtuw202mfx1G6uA3mU7C+LK6C6EyXrIijuXPn0t/+\n9jciIqqurqZ/+Zd/4cpmOk1kI779y/EFlc6YWYX6c+fOVRHVEEDlMU5AfEWSknNuy9pl/e1XHKwr\nCOuKIx7nX2Hw7JWXy+2LrlBUrkhUuTWMKlRX52PjRqChIcWt27boakm2LeNS9diwru78+GMYCcTj\n8qtNi9uNG9sYjKXw5ZdfIhaLgYhw0003hSwT1HBAltW7urro0qVLRET0zTff0D333ONY81SmIh6P\nq/18pEnOkg1yxmzZDHLGLPKXx9Hzzz9PeXl5lJ+fT+PGjaPhw4fTqlWrlGTdNv2fMfP9dcvwOBKd\n3bF0WTZVz9itPrJYYHGkyi1vDG5+/JwVueNPdqXgbI8r9WPZFEHEOU+2v66Bq1O2ep9d1q0m6rzA\ncjsqm+nOC0SSR7IBs1B/7ty5iMViKC0txcKFC30V6mvwUVtbi/b2drS2tuKtt97Cfffdhz/84Q/Z\nduu6gOY2M6CBXpOWo1CqY+7t7cWcOXOQl5eHHTt2OBUM0DpmWR2j/Sa9ar1oumtB7fjjH/+Il19+\nGdu3b/ctK/MtSB2zrObYrkOFa3s/1TpqVb5lHIm4NddUYHMk2pb7xOZPVp8cZR0z4H8MAJ/bSZMm\nYeTIkdyigMFSxxxFXlCaY7ZKjqy1BzTSj3nz5mHevHnZduO6hOY2Gnz44YcYP348zp49iwULFmDa\ntGmee0/2+yONjez7I1HOmWdiW+V+iTn+RrS0tCT5EEI4eULsVflV50mI9ByzKRf8FrFIVs8xW/2C\n8Qs9xyyUlc0x21FTU0ObN292yUc7x+wn5tNlU9WWiqyIW+kZs7Uq/6VLl7h9oq5j9ltXbNqU6S+z\nLDDvdBtG6k66H38NI4GGBvPO9saNQGVlG3i4cuUK5s2bh6tXr6K7uxuLFy9GXV0dt7+GOjS30eLy\n5cvo7e3Frbfemly9z1zuVCMtEP0S7Nixg5566qlk5s90LWjKhr/9sjarnXcWIdovA+vsSMRRV1cX\nEZnvUywtLaX9+/e79Knx6/+MWV2fjCP1M2Z1+yp6U/3YHVW45Z8xB/OF1d/fGbNaP79Q5dwrwxb6\n6quvqLi4mG6++Wa69dZbqba2liHPP2NOB8J8/zPlh1gu4BnzgQMHsH37duzcuRNXrlzBpUuXsHr1\nan13O40YPnw4AKC7uxu9vb347ne/m2WPrh9obqPDxIkTsXr1ahw6dAidnZ14/vnns+3SdQVhuVw6\nSo5Exf1RyOWazb6+PsRiMYwdOxbl5eUoKiry9LGerzcMc70BwzBtGUYChpGwTb/As23+n/IttZ3S\nYfedtW2Xt9tntScSXn8Aaz2BhMMfyz+ZfdNmSp9hNCbXG6iqqvLwZUGFWyC1loNhmGspWL66fZH5\nasrax+9t53GVamevBWH3yb02g51Hlr+pCiRee0pHY2OjbS2HKgZfJqzH3R977DGlqiE3cuk7mg2b\nyifhiUSCFi1a5Ot0nCjaSX7RzYywN/+C3nzxO5Vh4cKFC1RaWurhC4zLQfvNNPf0gMplsf3mH++S\nnTftYHHr/rDg7Rf37GfJuOVZNp0yYn5F3LLGafnqZwrCvs/uL2tMIj2i48Ifv+wmk+w4sb8voqkM\n2ePupjx/KkPf/AsxlQF4b6I8//zzvm6iuBfyiFou12xaGDlyJH784x/jz3/+sw9dYWwGkw03zmzY\nzDy3uceRP1nVx92B1CPvgFkut359DERlKCsrY97IV3u834K/R9StfUHK7cQLmKn5m9ZyOSLxTRRF\nFaHg94xZ1ma1i87agt588XPGfPbsWTp//jwREV2+fJnmzp1L+/btc+kTnTF7/1c/Y+a38Wy45f2f\nMfN9ddsXjU/ljFmVW79c8nx17wt+xszuJztjFkH1OLH3extUHnc35b03V2XfZT/pJMz3P9O22HJ8\nQaX1mMPcRMmleZ1M2+zo6MB9992XfNx90aJFmD9/vh+rvm2GlQ3DUSZtZpPbXOEoqE197yl6WaUn\n//r6+nDXXXfh73//O5588knPTRRZHbOsDtlqt18umPv8P11j1RGHrWP2s/qaEwk0NqZW6KqqagMP\nM2fOxOHRzBfuAAAgAElEQVTDh7ntGsGhuY0W1hTnuXPncPr0ad9TnBpiKK2VYeHixYuoqKhAfX29\nLRHK13KQrT0heyecn7UrRGsNuPtZYOng7ffjU8qX9LzzT7TmBGs8PC5E61yorGvBehkIjwtWP9V1\nSnjj8/oT/J1/sK2VIfJPde0K3rFR1SM7Zm6/eDrd+lly7nb294DP7eXLlzF8+HD09PTg3nvvxebN\nm3Hvvffa5L1rkajmAhX4WR8maJ902WLLhXznnwX7TRSN8Ghvb0d5eTlmzJiBO++8E//6r/+abZeu\nK2h+o4WuE48O0sT89ddfJy/Lv/32W+zduxclJSXKBoLPzwSVy525pKFDh6KhoQFHjhzBxx9/jN/+\n9rf47LPP/Fj1bTOsbC7Nn4bjN5hNILc4CmNTrU6cYzFHvqPZsimdY+7o6EBlZSX6+vrQ19eHVatW\n+byJosHDuHHjMG7cOADALbfcgunTp+PkyZOYPn16lj27PqD5jRY33HADWlpaklOc9lK0FKpgGPkA\nzHI5IAb3/RvrbTj2+zvWAzVWOVvqfpR5/8icAUjdH7L6A+z7WYB33R6VcrvyciAeB7M99aCVs7/9\nfll5eUvSn7SWyx07dozKysqoqKiIZsyYQU1NTcolH6k+4v2ysiKZvEinyCdxmZDYnopPKV/kHLW2\nttIdd9xBnZ2dLn0goJKqq6sJqKaGhgayHggwP3HPtvNBGW+788GUuOOBHFa788EQtj6rPR5Xt+/2\n122fL99AQDVVV1dTZWVlYH4tboHq/k+D62EW71jcvrK4c/e32s3/U/otruzbLK5ZXFr6RMfOrV/m\nDxCnhoYGqq6u7o83NW6JiH7961/TSy+95Nhn8uv8Psm+b/a/7O8SP1/I7LDsqvYR+c3zU2ZDxK3U\ntY6ODmpubiYios7OTpo6dSodPXpUSbnbWd7+wZ6YOzs7afbs2fTOO+8w9LFrQUUf9zhE/US67dsy\nfSxOZPZZPKvKOMcZjF934uBx4ObBvo/FHet/1vhYumXHQeaHjFORP2w5Nrd+6sRFMSIbP2+MonHw\n7LDsqvYR+c3zU2ZDFLfSOeZx48YhFosBcF4OqkLPMYtx7do1LF26FI888ggeeughv1YD2Qwnm1s2\ng/Mb3OZgmGNubm7GhAkTcNNNN+G2227Dbbfd5muKc7DkhaD++qrKaGtrQ3NzM0pLSx37rUV2zIV2\n7AvBWPM6KeeshV2sduciN+hvS/U3SXFusxZySSGRnEvi9bfrN2072y17qUVznP7a7VntTn2NMIwa\nAOJFdogIa9euRVFREdatW8ftpxEMmt/oMHPmTOzfvx9XrlzB119/jc8//9znjWsNIeQn8yZEl4Pu\n03v3KbzokkR0KcmSd8u6IbrMYcmLLi1lPqroZnFkYf/+/WQYBhUXF1MsFqNYLEa7du1yjUdPZcjH\nGYxfi1vemEUxJOKO9T9rfCzdsuMg80PGqcgfthybWzcWL16spzIY/UU2RNwqPfkX7nJbg4d7770X\nfX192XbjuoXmNzPgXUmbqAKQDyBVlWG9HYhV1QAAhlHWv99ZNeHun7rSTm3b5d1PGVtVHtZTu3Z5\nImdVhV3e0m/Zt6ouUlfy3m2Wv4bRiMpKsypj40aIIfvF6Ovro1WrVtG6deuUsr77l4L3/jPZmZf9\nnWuyMwB3G8um6CzB+riX/ZT5aO/Lk2VxZMeaNWtozJgxdOeddzLbAdEZc5zrn3wccY9OFres4yLj\nhG9f7C/fZ9k42fyqcsuLP1EMibiTLfvJO07m/3FmPxln7qVR3e1if+ICOX7sEqnduGaNg7eUqzyW\n49y+8riKO8bPs+OWY3Ers+v+nrl1y7iVzjF/+OGHeOONNxCPx1FSUoKSkhLs3r1bJqahiDVr1mg+\nI4LmNlroK+kIwU3Z/VA563D/wrB+HVjtsl8blrxb1usP2ybvrEX86yz2UUU3iyM3WltbA54xp3cc\nbhn3+FR1qdrnn9n4kxPxq8Kt6GzI7ZfbVxZ3rP9Z42Pplh0HmR8yTkX+sOXY3MqupE15/hkzz7Zq\nLLP2qcSVjC9VblW+CyLdIm6JFM6Y9VmHhoaGGwsXLsR//ud/4t/+7d/0lXQEkN78mzt3Ltra2gIb\nMCfQy4JIBpQLbpP9SGn0snKkHms1DPtjrQlbH8u2ecPEeszVfUPCWbpo6mDdQHHf0EjdkGmE/bFa\n1g2X8nK2Pp6/5upcIn9bAKxzbV9ATQ2wcWMbwqEK1s0pYBQMw8mtmxuLW6KU716ukNRhl+fd7HK2\nt8Aw7GN13gyz+tuXzDW3rTa+Pb4/dn9bUF19oX9fG3j4xS9+gf/3//4fVq9ejebmZm4/HrKRF4LK\nBvc1uE3+ubQNsstB87FY6zHO1CPDZrt7O07ux0Ddj6na97EeYxU9Epx6ZJetX7Tttu1PPm4bn/nI\nsPmpFF6yqPDLv2yKM/enuBB9vLJB5dR18GVT4w0yTj6/UXHrd5ys/uwxi/nl64xL7fmRDcsti1+7\nTX3zT8wtv8UGWXCLButnQGF0qOj2azt98mKagyeP9I4jmxymxhvEFp/fqLgN7mu4MatyLdMrkw3L\nbYrfSrJOUuzrvJjtcZKd9LjX+vArz9Mnkw+qX7yuTAPZ173RiTkNgR5Onk/zww8/TOPHj6dhw4ZR\nXl4ebdmyxcNvpsaRTQ5T4w1ii81vlNwG9zXcmFW5lumVycq4JQp3xqw6bl6/oLypHr+ojotK3BJR\nJhJz3NdgU33iAp0yMuJSO6wP6/JK/RMXjE+JZi6/IpvBg8orG1ROXQdfNjXeIOMMxq+eyuDLqnCb\n7qkMtXiK+4hVr6za8RPzE+a4qMattCpj+fLluOeee/DFF19gwoQJeP31133OYrfIu6RVLrisfY2N\nTNkMh2z4m2s2gyJbsZBLxyWExcDftcHBrTQxV1ZWYsSIEbjjjjvw9NNPY82aNT5NXJB3SatccFnr\nTS2ZtLl7925MmzYNU6ZMwaZNmzJiM5xs7tjMPW7DyGbW5vLlyzF79mwcOXIEQ4cOxU9+8hN/FgN/\n165/bgFJYu7t7cXTTz+N3bt34+jRo3jzzTf1ClJphOY3Omhuo8Ubb7yB2267DV999RUuX76M//u/\n/9P8phHCxHzw4EFMnjwZ+fn5GDp0KB5++GFs27bNp4m2gK4FlQsuG6ZeO4jN8Pz6txleNjds5ia3\nYWQzazMsv8G/a0Hlwshm3qYwMZ84cQITJkxIbufl5eHEiROefuZryo3kX+tjbm91bYs/qT5bBTpF\nck6bfj5btwaTs2zyxxeOX5FNORfqskHl1HXwZQ1DJCsbZ+a59TtOVn/2mNVjUBTzsmMpkxVxG5Rf\nu03V7xqPH7VY9cqqHb9w+UQmK+MWkDz5ZxhiYQAwb75qBIHmNzpobqOF5jdaCM+Yb7/9drS3tye3\n29vbkZeXF7lTgwWa3+iguY0Wmt+IwS2kI6Jr167RpEmTqLW1la5evUrFxcWOF7FqhIPmNzpobqOF\n5jdaCKcyhgwZgtdeew0VFRXo7e3F2rVrMX369Ez9Zlz30PxGB81ttND8Rowos/7mzZvJMAw6d+5c\ncl9tbS1NnjyZCgsL6b333vPIPPvsszRt2jSaNWsWLVmyhC5cuKAsu2vXLiosLKTJkydTfX290Ldj\nx45RWVkZFRUV0YwZM6ipqYmIiM6dO0f3338/TZkyhRYsWJB8RbsbPT09FIvFaOHChb7k0oVMc0uk\nzm9Ybolyj1/NrTp0XpAjssR87NgxqqiooPz8/OQBOHLkCBUXF1N3dze1trZSQUEB9fb2OuT27NmT\n3LdhwwbasGGDkmxPTw8VFBRQa2srdXd3Sy+tOjo6qLm5mYjM1+NMnTqVjh49Ss899xxt2rSJiIjq\n6+uT9t14+eWXacWKFbRo0SIiImW5dCDT3BL54zcst0S5x6/mVg06L6hxK03M58+fp6VLl9K0adNo\n+vTp9NFHHykpXrZsGX366aeOA1BbW+v4xaqoqBDqe/vtt2nlypVKsgcOHKCKiorkdl1dHdXV1Sn5\nSmS+5Xfv3r1UWFhIp06dIiLzIBUWFnr6tre30/z58+n9999P/jKqyLHw/e9/n2bOnEmxWIzuvvtu\nJZlMc0sUjl8/3BKlj99sxe5g4JYoN2I3V/OC9JHsZ555Bg8++CA+++wz/OUvf1GaR9q2bRvy8vIw\na9Ysx/6TJ0867tzyax9NbNmyBQ8++KCSrHpdpRf2t/yePn0aY8eOBQCMHTsWp0+f9vRfv349Xnrp\nJdxwQ4o+FTkWDMNAIpFAc3MzDh48KO2fDW6B4Pz65RZIH7/Zit3BwC2QG7Gbq3lBePPv4sWL2L9/\nf38xuDnhP3LkSADAggULcOrUKY/Miy++iLq6OuzZsye5j1z1jJbsiRMnsHfvXtTU1AAAamtrsWjR\noqSeYcOGYcWKFVz/7LWUKnWVLHzzzTdYunQpmpqacOutt3r0u/W+++67GDNmDEpKSmxvoPD65ccf\nHj9uZItb1rYK/HILpI/fKGK3vr4e9fX1AODgd7Bxa8dAj91czQsGudmxoaWlBU888QSKiorw6aef\nYvbs2WhqasLw4cMdxjTk4NE8adIkjBw5EjfeeCOeeOIJPP744452za8a3Pzq2E0fdOxGB276Fc1z\nfPLJJzRkyBA6ePAgERE988wz9Ktf/crRR6KCKisrleZU0iVHRARk3qZIVsTRyZMniYjozJkzVFxc\nTB988IGSLMC3qVJrw5J1r8frx6YKoooFFkdhY9eyGaRuaSBy5IZ9XJmOXbtNv/xGNU6ZzaD1a0G5\nFc4x5+XlIS8vD3fffTcAYNmyZTh8+LCvX4T8/Hxf/cPK9Utn3GZQ2fHjxwMARo8ejSVLlijN1YW1\nGUY2V2yGjd1sjDOM7GCJ3cHCrTAxjxs3DhMmTMAXX3wBANi3bx9mzJgRyJCGF5cvX0ZnZycAoKur\nC3v27MHMmTOz7NX1AR270ULHbrQQ3vwDgFdffRUrV65Ed3c3CgoKfL/BZNSoUYEcCyrXL51xm0Fk\nT58+jSVLlgAAenp6sHLlSjzwwAOR2gwrm0s2w8RuNsYZRnawxO5g4BZQSMzFxcU4e/YsRowYgdbW\nVixYsMDXJUssFgvkWFC5fumM2wwiO3HixFCvs8q0v7lmc/HixRgxYgRuvPFGHD9+PFmVEaXNbMkO\nltgdDNwCkqoMCxMnTsShQ4fw3e9+16vAMAbc8n6G0f+qwwECGUe9vb2YM2cO8vLysGPHDiVZ0RiD\njt+6iZ5uvVGDx5EobkVyzj4Dc8xhoTouEUeiuJXJ+vUjKNKlPwo/RfxIHzCxMNCS7/WEpqYmFBUV\n6fKiCKDjNjrouI0O0qkMwMzs999/P7desaqqKnn3cdSoUYjFYigrKwMANDY2Orat4mvZtrXP2i4v\nLwORqnwLDGOdj/5OW6L+hpFAPC72t6WlJfmySdkrdI4fP46dO3fihRdewCuvvCLs60YikUj64RdB\nZXPJpixuAX7sOh8SKEv6AfiPXT+x39LSgnXr1vmyZ9kqLwczNv34qxq7YeLWsp1L8WfFQKZsKlXn\nieoVZSri8biKCamcnzpCIB6o7lDFV55ekayIo2XLltHhw4cpkUgkn69XkQX4NlXGzpJVqWMOejx5\nNtMhx+MoTJ2tZTOqOEq3bDzuL+btfYPErixuRbJ2m3759cuP6jhlNoPWMQfNC0pnzKx6xblz5yol\n/qC/bkHl+qUzbjOIrMpjnAD7rA4oY56VGYalx9z2e9YGJJBIONvLywGi1PjsZwF2ecMA4nG+ft5Z\npF0/a9sOs70F1dUXsHEjUFnZxmDMRDbiNluymbSpGreAN3bXr4+ByBkLrFgVxRIAaaw5+5fZ5BKe\n2JJdDfPaZf3d3DQ2NqKlpUWttlmW8bu6uujSpUtERPTNN9/QPffc41jzVEFFWuD3bCAqt4Lo5XH0\n/PPPU15eHuXn59O4ceNo+PDhtGrVKiVZnh9hxs6TtfbJ9Iax60eP2x8WR7K45cnJfBvICHrGLO7n\n7agStzxZ1WOdrljzE1uyNr9xKoMo/qSqv/rqKyouLqbi4mKaMWMG1dbWKisn0lMZpozcmXRNZagm\nZj+y9kQoHqd/myw59zYvFkSJWRa3PDm3TT2VIVbqdyrD2iXjV/Y9C5KYZRyJbAZNzEG5lVZlWCVH\nhmFg4sSJeP755+Wn4RqBoO9upw9Wne2hQ4cwdOhQHDhwINsuXbfQcZt+KNUxv/LKKzh06BA6Ozux\nfft2p4IM1TH7qSOU1eNmyo+UTHCO/NYxhxk7T9ayJRt7mPppu5yqndRfPr9hYzeX6pj9fkdU+qY7\ndlWPdbpizU9sydoAf3Eq9y1EHbNVFvPYY49lJAFraKQLOnY1chXSqgxrVf5Lly5x+6jWMfu5k2rd\nmRfdueVvtwAIVgtqgXfnlVW1YJfxUwt65coVzJs3D1evXkV3dzcWL16Muro6bn83cqmmOBs2w8Su\n8456WdIPwFtRAqRq7O21xLwKE2+sp+RbWlqwfr28Bt/9XbKqcRIJUx/gp2ohd2I3bPwFqdgacHXM\nO3bsoKeeeio5ie3n5pQF++S33xt4zm1/sum6aaM64R90kr+rq4uIiK5du0alpaW0f/9+JVl98098\n8y9s7KreZHJz5uZI5Vj4uUHF02uPeZUYSMfNvyCxq2/+2eUDVmWEKediO6LcNRAJ9r5hS1nS4UdK\nRi7U1dVFc+bMoSNHjijJigI66NhVErNMPqhdP3pUEnO6YjdIYvYjn06ZoIlZ3C+9sas61vT5ry4n\na0tHLnDKB0zMdvgt52I7otx10CTm3t5eKi4upltuuYWee+45pmxlZSVVV1cTUE0NDQ2Os6N4PO44\nizDb4kl5d7t9295uccaSt7bt8jz77nZAbF9kzy5v1wc0JPkAKqUxGCZ2dWJOT+xWV6di14L7WIti\n1dq2x6Io1nn6RbFsj1WWP257sv7u7YaGBgcfaUvMixYt8iqQHGE9laGWAC5cuEClpaWMy3bY/k/9\ntQ682zfRl1J0+WqX5Y3ZblN0Zu2WSe1z2rTbFelxfnm8PKgkZr+xq6cyrH7piV23Xb9TGfbkaW0H\nGSc7Ltk23XKyH2GWXfH3ha9AWpVx5coVlJaW4plnnsGXX36p65gjwsiRI/HjH/8Yf/7zn7PtynUB\nK25jsRiefPJJ/faSCKFjNwLwc34Kokl+RRXcXw3VvtfjVMbZs2fp/PnzRER0+fJlmjt3Lu3bt48r\nKzsj8nPGLJIVnTGLdMnPmPl2VfX4OWMOemOV5wuvXTTGYGe/wWSCnjGL+6UndkV+q/jHiwO/4wxy\nxsyzp3LGLOorij+lRYysV753d3ejt7eXu/C4hj90dHSgsrISfX196Ovrw6pVqzB//vxsu3XdQMdt\ndNCxGy2UEnNfXx/uuusu/P3vf8eTTz6JoqIiR3tVVRW2bs1HdTW/jtmqrRTVcpow64ZNeOuYrbpi\nd62mc71ms47ZfFon1d9qt2pL7bWgpr5Ev05n3bKpx5RPZx3zzJkzfb913I5cqinul0Ym60hlcQsA\nhlGF6up8AM7YteLM7GOPrVQdsRUv7jpnq5bYii137PJWPzOMBBoa/NXgp+qkU7bt+gCr3ftdY9Vd\nD/zYTSBMTfF1UcfsBmuS31Ihmji32v1cKvBu/vEu15z/eyfrVS6zWTcNRXpY42SPLfi8il3WPXZ9\n8088lWFBdHOKz1Wc65/Tttxf1rjculiyIni/D16OZL7bbfLtpCd23X7rm39ibpXWyrDjN7/5DW6+\n+WY8++yzAFLPe6s8227+EIjbrT68Z9ztf+367P3da6q49bHWfnDbV/VLBbxn4tvb27F69WqcOXMG\nhmHgpz/9KX7+859zZXlj5/nq9UONf7c8jy+WfRnHPJ9Z9lh67P1TPsvXc3DHranbAEBCrlS4trez\n/BWNXybDgxovYt9VbKUrdt1+87Zl/dzjZY3Pjx7ePpYetz0/626wYz/EWhlff/118tLm22+/xd69\ne1FSUiIT01DA0KFD0dDQgCNHjuDjjz/Gb3/7W3z22WfZduu6gI7baKFjN1pIE3NHRwfuu+8+xGIx\nlJaWYtGiRb4m+Z1rDvhBULkwssFtBhnnuHHjkq83v+WWWzB9+nScPHkyUpthZcPYDMpvEJth4zbT\nsRBeNrM2sxe7QeWyE/NBZaU3/0aNGoWRI0cmL1m+853vBDKkIUZbWxuam5tRWlrqabMvtGMYowCY\nr5Yyb9wkADhvnpr9ykDkXGjH3G+/ueptt9/QYt3gamlpcfRn3VCy97dunlo3oFiLQtn7u1+NBZiv\nknL614Kamgv9221gYebMmdi2bVvycnvr1q34zne+47nc1ggPUezab65u3OiMXVmsum92mrHRAsMo\n69ee6P/L7m9tW/rNOEJS3h6rVqzbF4FKFQ4guW3XZ+9vToV4Xz0FtNhivxGVlS3YujXfw5MH/Olr\nEx0dHdTc3ExERJ2dnTR16lQ6evSoZwJbZRJeNlnuvinhbnf/5U3Es25+yCb7eTcTZH6pQEZzZ2cn\nzZ49m9555x2hrHtcojHL+HDa4HPGGruIJ/u2e7+qnOi4sP1j86sSu7xDI/LZ7ZcotngxyrMl6ieT\nYR03ke8qtsLGLu/YqsQqy2dRnLP9929TtM37Ptj/inQ42/jOS6cywl6yaIhx7do1LF26FI888gge\neuihbLtzXUHHbrTQsRsdlOqYLfAuWaqqqgDko6bG+SZcIFXHLLvcsF+qpqZl7PWk7nZz291uXT5Y\ntaCs/nZ7zu1E/z5zmwhwXy65L3dSUwmAtX609RZnQFwLSkRYu3YtioqKsG7dOm4/PhI2v3xKhqhj\nDmozuL9hbJrgX25XwTDy+/8fhXg8lqxpT92JN22Lpn2AVCxZ8Sa6VCZi62tsTNXgs2r02W8QL0v6\nYo91+6W3te39riA51ZT67oiniYBsxm5QuRyzKb4QSIF3yWKp4F2yha1jFl028C+z40qX4d5LFLac\nTI/dX/f47Ry5sX//fjIMg4qLiykWi1EsFqNdu3ZxZVn+8tuI20bkXWiHd6nm3rbXdDr9ZMeA+5iq\nyHm34xx98ktCInHs8sftjQUeV25/WPWyrBhi7Wfxy+PV+dfLkcpxtccRWyZc7PLHKufXqYt9XCSH\nPpBN7zbbJu+v21+vH/K4VapjvnbtGhYuXIh//Md/9Pw6suqYVetVWe1WH15tJ6+OkbVf1M7ry5KT\n6XGP2VvDmJ73pvHG5reNd2zcfVTHbtfh5oDFt0yO5zv/+PH5lcUu4JQTxZDbH56cKF7t+9zj9RNb\nfuqYVcfHlgsXu7DViavUWvPqk2Vj4vvg36Zo272f9VekUzVupXPMRGEvWTREePTRRzF27FjMnDkz\n265cd9CxGx103EYLaWL+8MMP8cYbbyAej6OkpAQlJSXYvXu3soHBUsccVHbNmjW++EyHTSC36piD\nyoWL3WA2gdyrYw4iGy5ug9kMJ5dbNqU3/7Zs2YLRo0ejt7cX//u//xvIiAYfc+fOFd4g1AgOHbvR\nQcdtxBBPnRN98MEHdPjwYbrzzjuZ7ZYK98S3t594ol51Yp0/Qe/dL2rn9ZXfCGCPRzR+Gc2tra1S\nfkVj89vmtRFu7Kx9Ir5lcjzfefp4/KrErmzcKjzzfBPxwBqvn9gKGtei8bFl+LEriltTn65jDhK3\nRArrMetfxuzD/uQfkHp6CmA/KWdts57cS5WBmdvupSlZJVb2bbc9/pOFTvt2fXb7rKepWPad+lpg\nGPKSLh27AwHOckQgJjm25jYgjm37trsc0P4knkg/v1Q3tc2KVVa7pV/cvxFVVS0A8iGF+PfGhOyM\nrrKykswXYzpfuGj+KqReHgqIX7ZpfuKOPtYLEO3t8u0GZrtdX8o/u3yc259nL/VCR7u/qZeFmtyI\naQ5+xhxn7FP92MfoX453diDWKfY3qJyIXxm35stcq/s/9lh1xoJ67KX28V4e6v4u8GLX3c56mWg6\n/U3Zt/ioDMxtit+gxzRcDKl8X1jxqyLn73sTd/R1yvG55bcoHgBLOZ8YLxluUtjkxD37wx48PkFi\nX2WBwPLXOTYxzToxh5MLl5iD2lTnljcu2TF1t7v1OdvD+5tObqPjVy2GZNzy4ldFThyrXll7X6cc\nn1t+i+IBkCdmPnG8vjIdQT98guT9RX1Y/jr38Wl++OGHafz48TRs2DDKy8ujLVu2MPlNNxf+A0w8\n9nToDO4Ln9/giTl93PLGJRuvu92tL6qYUOFWFreZ4DcI5+6+PI5V7fjt65Rjc0uk+M4/jejw5ptv\nZtsFDQ3f0HEbLaR1zMuXL8c999yDL774AhMmTMDrr7/u00QimGeDpF4xHMLYDCqbOzbDxW4wm9mT\nzYbNMAhqM6hcbtmUJubKykqMGDECd9xxB55++mmsWbPGp4mWQI4Fl8stm7t378a0adMwZcoUbNq0\nKSM2w8nmjs1wsZuNcYaRHSyxOzi4FSbm3t5ePP3009i9ezeOHj2KN998M8DrYy7Iu6RVLndshuf3\n+ucoqFxuchtGdrDE7vXPLSBJzAcPHsTkyZORn5+PoUOH4uGHH8a2bdsCGdLwQvMbHTS30ULzGy2E\nifnEiROYMGFCcjsvLw8nTpzwaaItgFth5HLHZnh+/dsML5sbNnOT2zCymbWZPX6DyuWWTWFVhiFb\nF5PRz/rXKbpV2O42k+qzFcHhlRXZZMmpDN+pUzxOr6x/fr0Ix5GiC0ybPFmxTr6/QeXYulQHFj23\n/nli8+vW55VP7/dFhOzyK4+hoN9vFTlxP5ZOfzkFkCTm22+/He3t7cnt9vZ25OXlOfqY5YoaQaD5\njQ6a22ih+Y0Y3ApnIrp27RpNmjSJWltb6erVq1RcXOx4maVGOGh+o4PmNlpofqOF8Ix5yJAheO21\n11BRUYHe3l6sXbsW06dPz9RvxnUPzW900NxGC81vxIgy62/evJkMw6Bz584l99XW1tLkyZOpsLCQ\n3jDPxWMAABmLSURBVHvvPY/Ms88+S9OmTaNZs2bRkiVL6MKFC8qyu3btosLCQpo8eTLV19cLfTt2\n7BiVlZVRUVERzZgxg5qamoiI6Ny5c3T//ffTlClTaMGCBXT+/HmmfE9PD8ViMVq4cKEvuXQh09wS\nqfMbllui3ONXc6sOnRfkiCwxHzt2jCoqKig/Pz95AI4cOULFxcXU3d1Nra2tVFBQQL29vQ65PXv2\nJPdt2LCBNmzYoCTb09NDBQUF1NraSt3d3dJLq46ODmpubiYi82WdU6dOpaNHj9Jzzz1HmzZtIiKi\n+vr6pH03Xn75ZVqxYgUtWrSIiEhZLh3INLdE/vgNyy1R7vGruVWDzgtq3Col5u9///s0c+ZMisVi\ndPfddyspXrZsGX366aeOA1BbW+v4xaqoqKCPPvqIq+Ptt9+mlStXKskeOHCAKioqktt1dXVUV1en\n5CsR0eLFi2nv3r1UWFhIp06dIiLzIBUWFnr6tre30/z58+n9999P/jKqyLnx+eefJ98wHIvFaMSI\nEclfaBEyzS1ROH79cEuUPn6DxC1ReH4HA7dEOi+4kU5upY9kA2ZpTCKRQHNzMw4ePCjtv23bNuTl\n5WHWrFmO/SdPnnTcuZXVPm7ZsgUPPvigkmyYusq2tjY0NzejtLQUp0+fxtixYwEAY8eOxenTpz39\n169fj5deegk33JCiT0XOjcLCQjQ3N6O5uRmHDh3C8OHDsWTJEqFMNrgFgvPrl1sgffz6jVsgPfwO\nBm4BnRfcSCe3yqvLkav0ZcGCBTh16pSn34svvoi6ujrs2bNHKnvixAns3bsXNTU1AIDa2losWrQo\nqWfYsGFYsWIF1ydn/XSgglx88803WLp0KZqamnDrrbd69Lv1vvvuuxgzZgxKSkq4L85kycmwb98+\nFBQUYMKECQOOW9a2CvxyC6SfXzc/QPDYra+vR319PQA4+B2s3AI6L1hIN7cGsSLXhUmTJmHkyJG4\n8cYb8cQTT+Dxxx93GNOQQ0bzo48+ijlz5uCpp55y7Nf8qoHFryhuAc2tKnixq/kND25eUJnvOHny\nJBERnTlzhoqLi+mDDz5ItslUVFZWKs2ppEtOJityNzqbYo6uXr1K3/ve9+jMmTO+ZKPyN51ydvfT\nYZNFB48jUdyK5Nw2/YIlq1r/NJC+L7LF3NPBr9+6sOsp5kX8KM0xjx8/HgAwevRoLFmyRHm+DoDt\nJaL+EFQu12wCwK5duzB79myMHj06YzZziaOgctmI22zJZsOmzgvRyUoT8+XLl9HZ2QkA6Orqwp49\nezBz5sxAxjTYePPNN7F8+fJsu3FdQcdttND8Rgvpzb/Tp08nKwV6enqwcuVKPPDAA8oGRo0aFcix\noHK5ZvPEiRN45513cPjwYWzatAlbtmzBD3/4w0hthpHNFZvZittsyWbaps4L0cpKE/PEiRPR0hJ8\nBf9YLJZRuVyz+ctf/hK/+93v8Oijj6KnpwddXV2R2wwjmys2sxW32ZLNtE2dF6KVVarKECowjJxa\nRcow+m9pZNQmm6OLFy+ipKQEX331lW/ZXEG6+WbpC8pRJrnNRtyFhVlUEZwjFX5zkRcZVMck4kep\njrm3txdz5sxBXl4eduzY4ctJDT5aW1sxevRorFmzBp9++ilmz56NpqYmDB8+3NGvqqoqeRNh1KhR\niMViKCsrA4BkzaR9u7wcIOK3Z3IbSCCRSK++xsYWXLhgvrKnra0NIujYjQ6a2wihUv7hfv5bteSD\niCgej6uYSJucTFbkbnQ22UY/+eQTGjJkCB08eJCIiJ555hn61a9+pSQrsqlyVDN1XOy+pMOmn3I5\nouCxm+5YUC0LG0jfF1m5nIhbU17Or99yuai+o+mUU415ET/Sqozjx49j586deOyxx3L6knogIi8v\nD3l5ebj77rsBAMuWLcPhw4ez7NX1Ax270UFzGy2kUxnW89+XLl3i9rEutTduBBoanJfagHk5GuZS\nNsilud22vd26tC4vLwORU39ZWZkv/wwDiMe97S0tapfa48aNw4QJE/DFF19g6tSp2LdvH2bMmMHt\n74adY78IKptLNv3ELuCcJpLFgjn/mkA8nmo3DOe2W97dbm3bY9EOWSwahrktil2V747dvj12Re+r\nU+HW9LEK1dX5AFL8mvbKbOMV++fetuZwWf7b+7v5tnSIvrvubb/8ucdjypt7Gxsb0dLSgq1b81Fd\nLaRNfCGxY8cOeuqpp5Kn5NaKSbzT8ahWd06nXkuXXWdQ/apyIprHjx9PN998M9100000atQoxzqz\nMtmwfmUC6fZFdSrDb+wG8cMtLlIn6h8m/mSyKrp5fXhTGSrcmvLgHC9//onkRTpkxyOILb/+iMYv\nij/hVMaBAwewfft2TJw4EcuXL8f777+P1atXS1K9E7wFPaKSyzWb3/nOd3D8+HF8++23OH/+PEaO\nHBm5zTCyuWIzbOxmY5xhZMPYBPzJ6rzgSzqQlDAx19bWor29Ha2trXjrrbdw33334Q9/+EMgQxp8\nkJ6jSzt07EYHzW30UK5j/uMf/4iXX34Z27dvdyqw1eJFVZOYTr2WLrvOoPrTUa+oskJXZWVlTpbL\nseZhw+o3jAQaGpxz+Fu3bhX+uKnErl9YC6fZxUXxIOofJv7cOll9ZLp5fVTqmHncmvIGAGLUnYf7\n7rllRP6Ljofad9c/f7LxpY47n9u0PmCiEzOvH/8AdHR0YPz48Th79iwWLFiAV199FXPnzlWSDetX\n1FBJHEF0DoQHTHRiVvFPJ+agiVlaLnflyhWUlpYiFouhqKgIzz//vEzEgcEylxRUNswKXYOFoyBy\n2YrbbMlmeh5U5wVl6UBS0sR80003IR6Po6WlBX/5y18Qj8fxP//zP4GMaTihV+iKDjpuo4XmN1r4\nmsq4fPky5s2bh61bt6KoqMhUoKcyFPqxL1laW1uTK3T97W9/Q35+Pj777DMl2XT4FTUGylQGK25V\n5GR+AHoqAxDzq6cyIprKAIC+vj7EYjGMHTsW5eXlDvI1gmNi/wpdq1evxj/90z9hypQp2XbpuoKO\n22ih+Y0OSosY3XDDDWhpacHFixdRUVGRfHrGgv3pqcZGZ9VAY2OjY1v0dJT9Lj5gPlGTmqNhP61j\nf/rJ0tfS0oJ169Y5+ltP51hP/rn1AWXJioYgVQRWJQSrakAE69HWF154Aa+88oqwrxvu45AJ2Vyy\nKYtbgP/kn/u4Wn4A/KdIU7EF5pN4QKL/bMm5bX9SzIpdlj3WdwUADKMM9u+Nu7871llPGhpGGeJx\n9Sf/ADV+gSrU1OQDSPFrfddScI7H4sf+ZCSApH/AOs931z5eq7/7+LBykSVv58PdnkhYnKf8c9t3\n5i67P6a/QOrJPyAf/e+Z5UP+XIsTv/71r+mll16yPcUC2//e/u5FPFSfjorH48lt0RNSrKefWAuk\nuHWw9Fo2VeHU77WZauMrXbZsGR0+fJgSiQT36bTKykqqrq6m6upqamhoSI7P9Dfu4NjaZ992LgIU\nd+xzt9u3AW+73b6KvN0/q43VP3UsvO0NDQ2OdnPb5AOoFPJrwR23pi6+nGyRHVb82OOI198dk+59\nqgvm2GUtnnn93ONwjwlgx33/RbiSPzx+eXrd/PL4sbcROfMCfzz87zerj+wYO+PZa1OUp1jHJdWP\nz62U9bNnz9L58+eJiOjy5cs0d+5c2rdvn82IODHznOK1qQ6YJ8OzpZKYZf750e/sx1aajseGeZz4\n6Z+Ovjx5VR0i/ljtzuPoFZLFrSkrdy5IYhb1lyVmVXgTs9h/lUTGtsFWrMqvLD79JmbeNms8KnlD\nJTGzttUTM1uv+Zd/4KVTGR0dHaisrERfXx/6+vqwatUqzJ8/XyamoQDr0dadO3fiypUruHTpElav\nXq2fokoDdNxGC81vxOCmbEXYVbC06akMtTMz0VQGDzx/Vc6YVS6ZVY6nTN6uQ7w2rdim3zNmNf/E\n3PL8Sdll+zQYpzLYtvVURtCpDKWbfxrR4MqVK5g3bx6uXr2K8+fPY9iwYdl2SUNDYyCAm7L7cezY\nMSorK6OioiKaMWMGNTU1ubK/+IyZ92vBa1P9JeLJ8GypnDHL/POj39mPr7Srq4uIiK5du0alpaW0\nf/9+ZVm+PX/909GXJ6+qQ8Qfq13ljNlP7Mr8Yu3n+STqLztjVoX3jFnsv8oZJttGMG5NeT3HLBq/\nKP6kIdHR0UHNzc1ERNTZ2UlTp06lo0eP2ozoxMzS7+wnV9rV1UVz5syhI0eO+JLViZnPkZ/YlfnF\n2j+YE7OMW1NeJ+agiVn6gMm4ceOSr+C+5ZZbMH36dJw8eVL5jFw/Ey+GSpF+VVUVampqUFNTg8bG\nRiQSZv2r6W8iuQ2k9hkGkn2c47L6J5L97e2ybcu+u92ybxhee4aR8tfebvrI7m/X39jY6GpvBFAD\nw6gBUOXhy0KY2M21tTKCxl9QWZ0XlKWDiXFTNgOtra10xx13UGdnpy37m3W2QHX/p4HsE95mzWnc\n1t9bV2vdfDA/Vl1rPLltv7HGqsu1t1t1rvb+9psbKf3OdiJ2XbCoTtjuT2qfad+qOza5kdN84cIF\nKi0tZdwsZcumfE9t29tEZ0GWnzKw5Hk3p2R+uP21t7vl3Xbd4/R+5PzyYheo7K+HNmM3daMnTqzY\nYscqOWSsfe5Yscs7bdi3G2z+8WvKvf7EPf44v0febfsY7d+dFB9qscvi1s6vOze4v+fu8bjHb/fX\nqml3fvfY43Xza9Xg89p5uUWWi8Tx0GBrb0jmSpNjPrfKibmzs5Nmz55N77zzjof8FLHsL5bqZYf7\ni+rnksT9v72/LAHI/GPZ5Mmx+6vR7OchCNZ47G3yxCz3x8/Ptkpi5rW75WV++03MotgVxY4oHlVj\nShbPKvHr/l/VHz/j4I8/GLemDkht8bgQtcl4Eo2L14eFMByL7KT+8rlVWivj2rVrWLp0KR555BE8\n9NBDwU7NNTz4+uuvk4+/fvvtt9i7dy9KSkqy7NX1BR270UFzGx2kiZmIsHbtWhQVFSWf4feDbMyZ\n5cpcUnNzMyZMmICbbroJt912G2677TZfRfrZmPvKpTm+cLEbzCaQpTnJDM8xh80LuRR/meYWUEjM\nH374Id544w3E43GUlJSgpKQEu3fvDmRMw4mZM2di//79uHLlCr7++mt8/vnnnmU/NYJDx2500NxG\ni7S9Wsq6K2/BrlW0PqlbzpJl6eOtbeqWcdty9+Ht4/nHssmTC/Pqo4ceegg/+9nPHGfNPFm3Pbcv\nFlR858HPWrkyPyxfWO2848LTz7CuxK/X59R6wazYcfvNgiymeH6794ni1/2/2z+eP37GYe/jlAv/\naimRLRY/vO+vXYbVz+m311bKL3G7qJ/bb3c/0TH3todcj1kjerS1taG5uRmlpaWeNqtczjBqHOVq\n1lKD7vIz++WTVY5mlc/Z+wPecjerpM1e/maVr7HK2+yld5Z+u7zpS6q/2769PC619KOz3M/bbpbL\nmZ8qD18aGjkP7m3BfqxZs4bGjBlDd955J+fOJZh3IC2wSqSc8rw7l3HhnVC3vP1/3nP4sn2stSd4\nNt1tYdbKkN3ZZvli+h/3jEd0J9jNL2tMPF12m14f1T68cjm573FJPza/KrHLige3TZGvdu5Y4+TJ\nu/ex+PX3vfHa5Ntgy7LlgnFr6hBVZcS5/LjHzIohHjei4+Re84Jl0+m/lx9VjlP7+bKivCA9Y16z\nZo2eO4oQ+s52dNCxGx00txGDm7JtaG1tDXzGzPp1cbfJziZYv1Jueff/9v6isxsV/1g2eXLs/mya\n+/r6qKCggG6++WYpvyx7Kmd0orMz3jhlOrw+qn1U7KjEgvfDD2NZ7IpiR8UH+5h4+1jyIht2H1j/\nq/rjZxz88Qfj1tSh65jFY+Zzm5bV5aqqqgDk92+NAhDrf90N+/Ur9tfp2F8dZW1br36yt1v63K93\nYcm7X9cjsp+av3Tat78uiOWv2x+7/sbGFqxfb9YnV1a2gYcPP/wQX331FaZMmYIvv/wSJSUlqKur\nw49+9COujEa6kXr1ETAKiUTM8WoiwBkbJlLbrFh1x5Yo1lnb7lhkvQqJJ2/FKuu7oSbfgupqtVdL\nqaEK7twg+u65t938ur/b5eXmNi8XWNt+263vPq+/3X9V/oFGGIb5ainZTVh+yrYhyBmz9ZE/Tsv7\nVYkr/dqyt+PcfuKP09cgsqxfdRnNKvzafUrpV+OWzVXcZcetm6fDKeeHo0zPMatwyz5eKjad/e06\nnI/pssfKPiZeu+pcqfmrIuu0F4xbO79h/PWOPS5oSx9HTt3y4xnMJp9bXZWhoaGhMcAQ+UL5Qd+o\nnDr9z6RsNmzKUVVVha1b8wEAhmG/HCyD7HKaNRVjv9yKx9UvBy2Yl3B8e7xt663n7jdFi+TtK82l\n2lsApPNym4UyaQ+uZOCYD2M3GzbDIKjNoHI5ZpN7Lt2Phx9+mMaPH0/Dhg2jvLw82rJli+dyRXRq\nb7X5Of1X6evuw5Pxf7kTzGeeTTtHPISZKgrjpx/ZMDai18fmVyV22ccrXKy4dbB0qvSxDyud3Ps7\nNsG4tfMbVZxEyUtm7LK5JSL11eW4CvqV84zrOWbxASAKk5j9cxtUNtVPfFyCcORHzm+AiyBOzGrj\ntPrbdQyWOWZVfsP46x170Bjyx5FTt/x4BrPJ51Y6x7x7925MmzYNU6ZMwaZNm3yfkbe0tPiW6ZcM\nKBdGNvM2y8rKMHnyZPz1r3/FqFGj8Prrr0duM5xs7tgMF7vBxxk85sPYzbzN7PA7OLgVJube3l48\n/fTT2L17N44ePYo333zT9yI71rKW/hFULoxsZm329vbi+PHj+PLLL9Hd3Y38/Hz88Ic/jNRmeNnc\nsBk+doOPM3jMh7Gb+djNDr/XP7eAJDEfPHgQkydPRn5+PoYOHYqHH34Y27ZtC2RIwwvNb3TQ3EYL\nzW+0ECbmEydOYMKECcntvLw8nDhxwpeBtra2QI6Fu9ueGzbD8+vfZnjZ3LCZTW6Dx3wYu5m1mT1+\ng8rllk1huZwhfTxF3m/rVgDYKn/SxaEvKa3Qh71tyfqxa8mlVlbzL+uWE+lIB79+uQ0q6+zHPy4q\nNr361OVUocotYAiOl9wmK1asmGfrZO/j2fXHU/jjomrTD78qNvl2+HLB4l6NI14+CWbX/3ERJubb\nb78d7e3tye329nbk5eU5+pg3XzWCQPMbHTS30ULzGzFE5S7Xrl2jSZMmUWtrK129epWKi4vp6NGj\ngctnNJzQ/EYHzW200PxGC+EZ85AhQ/Daa6+hoqICvb29WLt2LaZPn56p34zrHprf6KC5jRaa34gR\nZdbfvHkzGYZB586dS+6rra2lyZMnU2FhIb333nsemWeffZamTZtGs2bNoiVLltCFCxeUZXft2kWF\nhYU0efJkqq+vF/p27NgxKisro6KiIpoxYwY1NTUREdG5c+fo/vvvpylTptCCBQvo/PnzTPmenh6K\nxWK0cOFCX3LpQqa5JVLnNyy3RLnHr+ZWHTovyBFZYj527BhVVFRQfn5+8gAcOXKEiouLqbu7m1pb\nW6mgoIB6e3sdcnv27Enu27BhA23YsEFJtqenhwoKCqi1tZW6u7ull1YdHR3U3NxMROYbRKZOnUpH\njx6l5557jjZt2kRERPX19Un7brz88su0YsUKWrRoERGRslw6kGluifzxG5ZbotzjV3OrBp0X1LiN\nLDEvW7aMPv30U8cBqK2tdfxiVVRU0EcffcTV8fbbb9PKlSuVZA8cOEAVFRXJ7bq6Oqqrq1P2d/Hi\nxbR3714qLCykU6dOEZF5kAoLCz1929vbaf78+fT+++8nfxlV5NKFTHNLFI5fP9wS5T6/mls+dF5Q\n4zaSZT+3bduGvLw8zJo1y7H/5MmTjju3strHLVu24MEHH1SSDVNXaX8R6unTpzF27FgAwNixY3H6\n9GlP//Xr1+Oll17CDTek6FORSweywS0QnF+/3AK5z6/mlg2dF9S5Dbzs54IFC3Dq1CnP/hdffBF1\ndXXYs2dPch+5ymYs2RMnTmDv3r2oqakBANTW1mLRokVJPcOGDcOKFSu4PthrKdXrKp345ptvsHTp\nUjQ1NeHWW2/16HfrfffddzFmzBiUlJS43vYslvODgcYta1sFfrkFBja/9fX1qK+vBwAHv5rbFAZa\n7OZqXgicmPfu3cvc/9e//hWtra0oLi4GABw/fhyzZ8/Gn/70p2TtoyX7ox/9CBs3bkRpaalDx+9/\n/3vs3LkT//3f/53c566bPH78OG6//XZuO6uu0g3rRairVq1Kvgh17NixOHXqFMaNG4eOjg6MGTPG\nIXPgwAFs374dO3fuxJUrV3Dp0iWsWrVKKucHA41bVh8Zv0G4BQY2vytXrsQvfvELAGx+NbcDL3Zz\nNi8oT7YEBGuS/+rVq/TVV1/RpEmTqK+vz9F/165dVFRURGfPnnXsl8n6ravs6+ujVatW0bp16xz7\nn3vuueScVV1dnXCyPpFIJOeS/MilC5nilsgfv+nglii3+NXc+oPOC2JEnpgnTpzoKIt58cUXqaCg\ngAoLC2n37t2e/pMnT6Y77riDYrEYxWIxevLJJ5Vld+7cSVOnTqWCggKqra0V+rV//34yDIOKi4uT\ntnbt2kXnzp2j+fPnK5W3JBKJ5N1XP3LpQia5JVLnNx3cEuUWv5pbf9B5QQyDSD83qaGhoTGQoF/G\nqqGhoTHAoBOzhoaGxgCDTswaGhoaAww6MWtoaGgMMOjErKGhoTHAoBOzhoaGxgDD/weGiuGgcrm5\n8QAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Benchmarks"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Pure Python version."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%timeit -n 2 -r 2\n",
      "correlograms = compute_correlograms(spiketimes, neurons,\n",
      "                                     ncorrbins=ncorrbins, corrbin=corrbin)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2 loops, best of 2: 21.8 s per loop\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Cython version."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%timeit -n 2 -r 2\n",
      "correlograms = compute_correlograms_cython(spiketimes, neurons,\n",
      "                                    ncorrbins=ncorrbins, corrbin=corrbin)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2 loops, best of 2: 788 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 9
    }
   ],
   "metadata": {}
  }
 ]
}