mixture_distribution_quantiles.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as sp\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# On Chip Generation, Estimation of Extreme Values\n",
    "\n",
    "When generating projection data on chip, some values that are required in order to allocate memory and estimate performance must be treated as random variables. For memory allocation, typically we want to calculate extreme quantiles of these random variables, i.e., values large enough that the true values exceed them with very low probability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantiles of the number of synapses per sub-row\n",
    "\n",
    "One quantity that we will want an extreme quantile of is the number of is $N_\\texttt{conns}$, the number of connections per sub-row of the connection matrix of a projection.\n",
    "\n",
    "In the remainder of the document, we use the following example values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "post_slice_size = 1024\n",
    "pre_size = 5000\n",
    "post_size = 5000\n",
    "full_matrix_size = pre_size * post_size\n",
    "\n",
    "quantile = 0.9999"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means we are dealing with a projection from a population of 5000 pre-synaptic neurons to a population of 5000 post-synaptic neurons, and that there are $5000^2$ possible connections. The full synaptic matrix contains $5000^2$ entries. Each row in the matrix specifies the post-synaptic neurons that a particular pre-synaptic neuron is connected to. The synaptic matrix has been partitioned, and we are dealing with a partition that is 1024 entries wide.\n",
    "\n",
    "We want to calculate a value for the number of synapses appearing in a row within our partition such that, with probability 0.9999, this value will not be exceeded for any row in the matrix.\n",
    "\n",
    "The probability distribution over the number of synapses per sub-row (i.e., per row within our partition) depends on the connector type that is used. For example, the $\\texttt{FixedProbabilityConnector(p)}$ connects each pre-synaptic neuron to each post-synaptic neuron with probability $p$. The number of possible connections in the sub-row is $\\texttt{post_slice_size}$, so the number of connections is distributed as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$N_{\\textrm{conns}} \\sim \\textrm{Binomial}(\\texttt{post_slice_size}, p)$.\n",
    "\n",
    "The following figure shows the probability mass function of this random variable with $p=0.05$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0xb5aaa1ec>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAF5CAYAAABEPIrHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X2UXXV97/H3R0At2sZ7ixAsta3PUFsh8aGpXrVwEakt\nVK3Sab22RVFabLmh66ql7U2C1VpdJJUqtyz6gFEZFxaLepXGgl5rJYgSgSoBsaLVekihanwAKw/f\n+8feo4fhnCQzmcn5JfN+rXVWcn77u/f57d/aM/M5+zFVhSRJUivuN+kOSJIkDTOcSJKkphhOJElS\nUwwnkiSpKYYTSZLUFMOJJElqiuFEkiQ1xXAiSZKaYjiRJElNMZxIkqSmNBNOkpyW5OYkdyS5MsmT\ndlL/giRb+/prkxw/oubwJO9J8vUk30ry8SSHLd5aSJKk3dVEOElyEnA2sAY4CrgW2JTkoDH1q4AL\ngfOBI4FLgEuSHDFU80jgo8D1wNOBnwJeA3xn8dZEkiTtrrTw4L8kVwIfr6rT+/cBvgScU1VvGFH/\nTuDAqjphqG0z8Kmq+u3+/TTw3ar69T2xDpIkaWFMfM9JkgOAlcDlM23VJabLgFVjZlvVTx+2aaa+\nDzfPAW5K8vdJtvWHik5c6P5LkqSFNfFwAhwE7Adsm9W+DVg+Zp7lO6k/GHgw8CrgA8CxwN8B707y\n30YtMMmBSVYkOXDOayBJ0hK20H9D91+IhSySAHM55jRcPxO6Lqmqc/r/X5fkZ4FT6c5Fme1I4GPA\nliTfmjXt7+n2zEiStNQdBzx7VtuDgRXAU4ErdvcDWggntwF3A4fMaj+Y++4dmXHLTupvA+4Cts6q\n2Uo3cKP8eP/vihHTng68bsx8kiSp8+PsC+Gkqu5McjVwDPBe+N45I8cA54yZbfOI6cf27TPL/ATw\n2FnzPQb44phlfgHg7W9/O4cffvjcV2QJW716NRs2bJh0N/Yqjtn8OG5z55jNj+M2N1u3buVFL3oR\n9H9Ld9fEw0lvPfDWPqRcBawGDgQuAEiyEfhyVZ3Z178J+EiSM4D3A1N0J9WeMrTMNwLvTPJR4MPA\n8cAvAM8Y04fvABx++OGsWDFq54nGWbZsmWM2R47Z/Dhuc+eYzY/jNm8LcruOJsJJVV3U39PkLLrD\nNdcAx1XVrX3JYXSHaWbqNyeZAl7bv24CTqyq64dqLklyKnAmXZi5EXheVW3eE+skSZLmp4lwAlBV\n5wLnjpl29Ii2i4GLd7LMC+j3vkiSpL1DC5cSS5IkfY/hRLttampq0l3Y6zhm8+O4zZ1jNj+O22Q1\ncfv6FiRZAVx99dVXexKUJElzsGXLFlauXAmwsqq27O7y3HMiSZKaYjiRJElNMZxIkqSmGE4kSVJT\nDCeSJKkphhNJktQUw4mkfcJgMGDt2rUMBoNJd0XSbjKcSNonDAYD1q1bZziR9gGGE0l7BfeMSEuH\n4UTSXsE9I9LSYTiRtM9zr4u0dzGcSNrnuddF2rsYTiRJUlMMJ5IkqSmGE0mS1BTDiaQmeNKqpBmG\nE0lN8KRVSTMMJ5IkqSmGE0mS1BTDiSRJaorhRJIkNcVwImnJ80ohqS2GE0lLnlcKSW0xnEiSpKYY\nTiRJUlMMJ5IkqSmGE0mS1BTDiSRJaorhRNIe4eW6knaV4UTSHuHlupJ2leFEkiQ1xXAiSZKaYjiR\nJElNMZxIkqSmGE4kSVJTDCeSJKkphhNJ2gnv0SLtWYYTSdoJ79Ei7VnNhJMkpyW5OckdSa5M8qSd\n1L8gyda+/tokx8+a/jdJ7pn1+sDiroUkSdpdTYSTJCcBZwNrgKOAa4FNSQ4aU78KuBA4HzgSuAS4\nJMkRs0ovBQ4BlvevqUVZAUmStGCaCCfAauC8qtpYVTcApwK3AyePqT8duLSq1lfVjVW1BtgCvGJW\n3X9W1a1V9e/9a/uirYEkSVoQEw8nSQ4AVgKXz7RVVQGXAavGzLaqnz5s04j6ZybZluSGJOcm+a8L\n1G1JkrRIJh5OgIOA/YBts9q30R2KGWX5LtRfCrwYOBp4JfAM4ANJsrsdliRJi2f/SXdgBwLUfOur\n6qKhaZ9J8s/AvwDPBD48biGrV69m2bJl92qbmppiasrTVSRJmp6eZnp6+l5t27cv7FkTLYST24C7\n6U5cHXYw9907MuOWOdZTVTcnuQ14FDsIJxs2bGDFihU767MkSUvSqC/sW7ZsYeXKlQv2GRM/rFNV\ndwJXA8fMtPWHXo4Brhgz2+bh+t6xfftISQ4DfhjwRgWSJDVs4uGktx54WZIXJ3kc8BfAgcAFAEk2\nJnndUP2bgOOTnJHksUnW0p1U++a+/kFJ3pDkKUl+LMkxdJcbf5buxFlJi8A7qUpaCE2Ek/78kN8D\nzgI+Bfw0cFxV3dqXHMbQya5VtZnuniUvA64BngecWFXX9yV398t4D3Aj3f1QPgE8vd9TI2kReCdV\nSQuhhXNOAKiqc4Fzx0w7ekTbxcDFY+q/Azx7QTsoSZL2iCb2nEiSJM0wnEiSpKYYTiRJUlMMJ5Ik\nqSmGE0mS1BTDiSRJaorhRJJ2kzefkxaW4USSdpM3n5MWluFEkiQ1xXAiSZKaYjiRJElNMZxIkqSm\nGE4kSVJTDCeSJKkphhNJktQUw4kkSWqK4USSJDXFcCJJkppiOJEkSU0xnEiSpKYYTiTtMp++K2lP\nMJxI2mU+fVfSnmA4kSRJTTGcSJKkphhOJElSUwwnkiSpKYYTSZLUFMOJJElqiuFEkhaR94aR5s5w\nIkmLyHvDSHNnOJEkSU0xnEiSpKYYTiRJUlMMJ5IkqSmGE0mS1BTDiSRJaorhRJIkNcVwIkmSmmI4\nkSRJTTGcSJKkpjQTTpKcluTmJHckuTLJk3ZS/4IkW/v6a5Mcv4Pa85Lck+R3F77nkiRpITURTpKc\nBJwNrAGOAq4FNiU5aEz9KuBC4HzgSOAS4JIkR4yo/SXgycC/LU7vJUnSQmoinACrgfOqamNV3QCc\nCtwOnDym/nTg0qpaX1U3VtUaYAvwiuGiJD8CnAP8KnDXovVekiQtmImHkyQHACuBy2faqqqAy4BV\nY2Zb1U8ftmm4PkmAjcAbqmrrQvZZkiQtnomHE+AgYD9g26z2bcDyMfMs34X6VwPfrao3L0QnJUnS\nntFCOBknQM2nPslK4HeB31yEfkn7rMFgwNq1axkMBpPuiqQlbP9JdwC4DbgbOGRW+8Hcd+/IjFt2\nUv804KHAl7qjO0C3d2Z9kv9ZVY8Y15nVq1ezbNmye7VNTU0xNTW1k9WQ9n6DwYB169ZxwgkncOih\nh066O5IaND09zfT09L3atm/fvqCfMfFwUlV3JrkaOAZ4L3zvfJFj6E5mHWXziOnH9u3QnWvyD7Pm\n+WDf/jc76s+GDRtYsWLFXFZBkqQlY9QX9i1btrBy5coF+4yJh5PeeuCtfUi5iu7qnQOBCwCSbAS+\nXFVn9vVvAj6S5Azg/cAU3Um1pwBU1deArw1/QJI7gVuq6qZFXxtJkjRvTYSTqrqov6fJWXSHa64B\njquqW/uSwxi6FLiqNieZAl7bv24CTqyq63f0MYvSeUmStKCaCCcAVXUucO6YaUePaLsYuHgOyx97\nnokkSWpHy1frSJKkJchwIkkT5OXb0n0ZTiRpgmYu3zacSN9nOJEkSU0xnEiSpKYYTiRJUlMMJ5Ik\nqSmGE0mS1BTDiSRJaorhRJIkNcVwIkmSmmI4kSRJTTGcSJKkphhOJElSUwwnkiSpKYYTSZLUFMOJ\nJElqiuFEkiQ1xXAiSZKaYjiRJElNMZxIS8hgMGDt2rUMBoNJd0WSxjKcSEvIYDBg3bp1hhNJTTOc\nSJKkphhOJElSUwwnkiSpKYYTSWqYJzFrKTKcSFLDPIlZS5HhRJIkNcVwIkmSmmI4kSRJTTGcSJKk\nphhOJElSUwwnkiSpKYYTSZLUFMOJJElqiuFEkiQ1xXAiSZKaYjiRJElNMZxIkqSmGE4kSVJTDCeS\nJKkpzYSTJKcluTnJHUmuTPKkndS/IMnWvv7aJMfPmr6mn/6tJF9N8g9Jnry4ayFJknZXE+EkyUnA\n2cAa4CjgWmBTkoPG1K8CLgTOB44ELgEuSXLEUNmNwGnA44GnAl8APpjkhxdpNSRJ0gJoIpwAq4Hz\nqmpjVd0AnArcDpw8pv504NKqWl9VN1bVGmAL8IqZgqp6Z1V9qKq+UFVbgTOAHwJ+elHXRJIk7Zbd\nDifp7cb8BwArgctn2qqqgMuAVWNmW9VPH7ZpXH3/GS8Hvk63V0aSJDVq3uEkyUuSfBr4DvCdJJ9O\n8tJ5LOogYD9g26z2bcDyMfMs35X6JM9J8s2+j6cDx1bVV+fRR0mStIfsP5+ZkpxFd5jkz4HNffMq\nYEOSh1fV/16AvgWo3az/EPAEugB0CvCuJE+uqtsWoH+SJGkRzCucAL8FnFJV00Nt701yHV1gmUs4\nuQ24GzhkVvvB3HfvyIxbdqW+qu4APt+/rkryWeAlwJ+O68zq1atZtmzZvdqmpqaYmpra8VpIjRgM\nBpx33nm8/OUv59BDD510dyTtY6anp5menr5X2/bt2xf0M+YbTg4APjmi/eq5LrOq7kxyNXAM8F7o\nzmPp358zZrbNI6Yfy/f34oxzP+ABOyrYsGEDK1as2IWeS20aDAasW7eOE044wXAiacGN+sK+ZcsW\nVq5cuWCfMd9zTt5Gt/dktpcB75jH8tYDL0vy4iSPA/4COBC4ACDJxiSvG6p/E3B8kjOSPDbJWrqT\nat/c1x+Y5LVJnpLk4UlWJPlr4GHAu+bRP0mStIfMd88JwEuSPAu4sn//M8CPAhuTrJ8pqqozdrag\nqrqov6fJWXSHa64BjquqW/uSw4C7huo3J5kCXtu/bgJOrKrr+5K7gccBL6Y73+Q/gE8AT+svK5ak\nfYKH8bQvmm84eTzdfUUAHtn/e2v/evxQ3S6f0FpV5wLnjpl29Ii2i4GLx9T/J/D8Xf1sSdpbeRhP\n+6J5hZOq+rmF7ogkSRK0c4dYSZIkwHAiSZIaYziRJElNMZxIkqSmGE4kSVJTDCeSJKkphhNJktQU\nw4kkSWqK4USSJDXFcCJJkppiOJEkSU0xnEiSpKYYTiRJUlMMJ5IkqSmGE0mS1BTDiSRJaorhRJIk\nNcVwIkmSmmI4kSRJTTGcSJKkphhOpL3MYDBg7dq1DAaDSXdFjXNb0d7KcCLtZQaDAevWrfMPjnbK\nbUV7K8OJJElqiuFEkiQ1xXAiSZKaYjiRJElNMZxIkqSmGE4kSVJTDCeSJKkphhNJktQUw4kkSWqK\n4USSJDXFcCJJkppiOJEkSU0xnEiSpKYYTiRJUlMMJ5IkqSmGE0mS1BTDiSRJaorhRJIkNaWZcJLk\ntCQ3J7kjyZVJnrST+hck2drXX5vk+KFp+yf50yTXJflWkn9L8tYkhy7+mkiSpN3RRDhJchJwNrAG\nOAq4FtiU5KAx9auAC4HzgSOBS4BLkhzRlxzYt6/rl/dc4LHAexZxNSRJ0gJoIpwAq4HzqmpjVd0A\nnArcDpw8pv504NKqWl9VN1bVGmAL8AqAqvpGVR1XVRdX1U1VdVU/bWWSwxZ/dSSpfYPBgLVr1zIY\nDCbdFeleJh5OkhwArAQun2mrqgIuA1aNmW1VP33Yph3UAzwEKODr8+6sJO1DBoMB69atM5yoORMP\nJ8BBwH7Atlnt24DlY+ZZPpf6JA8AXg9cWFXfmn9XpcXnt1lJS93+k+7ADoRuT8du1SfZH3hXP+23\nd7aQ1atXs2zZsnu1TU1NMTU1NYeuSPM38232hBNO4NBDPYdbUlump6eZnp6+V9v27dsX9DNaCCe3\nAXcDh8xqP5j77h2Zccuu1A8Fkx8Fjt6VvSYbNmxgxYoVu9BtSZKWnlFf2Lds2cLKlSsX7DMmflin\nqu4ErgaOmWlLkv79FWNm2zxc3zu2b59ZxkwweQRwTFV9bQG7LUmSFkkLe04A1gNvTXI1cBXd1TsH\nAhcAJNkIfLmqzuzr3wR8JMkZwPuBKbqTak/p6/cDLqa7nPgXgAOSzOxp+WofiCRJUoOaCCdVdVF/\nT5Oz6A7XXAMcV1W39iWHAXcN1W9OMgW8tn/dBJxYVdcP1f9C//9r+n9nzkn5OeAfF3F1JEnSbmgi\nnABU1bnAuWOmHT2i7WK6vSOj6r9IdwWQJEnay0z8nBNJkqRhhhNJktQUw4kkSWqK4USSJDXFcCJJ\nkppiOJEkSU0xnEiSpKYYTiRJUlMMJ5IkqSmGE0mS1BTDiSRppMFgwNq1axkMBpPuipYYw4kkaaTB\nYMC6desMJ9rjDCeSJKkphhNJktQUw4kkSWqK4UTawzzJUJJ2zHAi7WGeZChJO2Y4kSRJTTGcSJKk\nphhOJElSUwwnkiSpKYYTSZLUFMOJJElqiuFEkiQ1xXAiSZKaYjiRJElNMZxIkqSmGE4kSfPic6K0\nWAwnkqR58TlRWiyGE0mS1BTDiSRJaorhRJIkNcVwIkmSmmI4kSRJTTGcSIvASywlaf4MJ9Ii8BJL\nSZo/w4kkSWqK4USSJDXFcCJJkppiOJEkSU1pJpwkOS3JzUnuSHJlkiftpP4FSbb29dcmOX7W9Ocm\n+fsktya5J8lPL+4aSJKkhdBEOElyEnA2sAY4CrgW2JTkoDH1q4ALgfOBI4FLgEuSHDFU9iDgn4BX\nAbV4vZckSQupiXACrAbOq6qNVXUDcCpwO3DymPrTgUuran1V3VhVa4AtwCtmCqrq7VX1x8DlQBa3\n+5KkYd7rR7tj4uEkyQHASroQAUBVFXAZsGrMbKv66cM27aBekrQHea8f7Y6JhxPgIGA/YNus9m3A\n8jHzLJ9jvSRJ2ku0EE7GCXM7V2Su9ZIkqUH7T7oDwG3A3cAhs9oP5r57R2bcMsf6XbZ69WqWLVt2\nr7apqSmmpqZ2d9GSJO31pqenmZ6evlfb9u3bF/QzJh5OqurOJFcDxwDvBUiS/v05Y2bbPGL6sX37\nyI/Z1f5s2LCBFStW7Gq5JElLyqgv7Fu2bGHlypUL9hkTDye99cBb+5ByFd3VOwcCFwAk2Qh8uarO\n7OvfBHwkyRnA+4EpupNqT5lZYJL/Ajwc+BG6Qz6P60PPLVW123tYJEnS4mginFTVRf09Tc6iO1xz\nDXBcVd3alxwG3DVUvznJFPDa/nUTcGJVXT+02BOAv6Hba1LAzD6odf3nSJKkBjURTgCq6lzg3DHT\njh7RdjFw8Q6W91bgrQvWQWnIYDDgvPPO4+UvfzmHHnropLsjSfuUlq/WkZrlPRwkafEYTiRJUlMM\nJ5IkqSmGE0nSHuezd7QjhhNJ0h7neVvaEcOJJElqiuFEkiQ1xXAiSZKaYjiRJElNMZxIkqSmGE4k\nSVJTDCfSCN6DQZImx3AijeA9GCRpcgwnkqTmuPdyaTOcSJKa497Lpc1wIkmSmmI4kSRJTTGcSJKk\nphhOJElSUwwnkiSpKYYTSZLUFMOJJElqiuFES5I3eJL2bv4M79sMJ1qSvMGTtHfzZ3jfZjiRJElN\nMZxIkqSmGE4kSVJTDCeSJKkphhNJ0j7FK3n2foYT7bP8BSUtTV7Js/cznGif5S8oSdo7GU4kSVJT\nDCeSJKkphhNJ0pLi+WjtM5xIkpYUz0drn+FEey2//UjSvslwor2W334kad9kOJEkaYh7ZSfPcCJJ\n0hD3yk6e4UTN8tuLJC1NhhM1y28vklrkF6fF10w4SXJakpuT3JHkyiRP2kn9C5Js7euvTXL8iJqz\nknwlye1J/iHJoxZvDZau6enpSXdBkvYYvzgtvibCSZKTgLOBNcBRwLXApiQHjalfBVwInA8cCVwC\nXJLkiKGaVwGvAF4OPBn4dr/M+y/iqixJ8w0nfvuQtK/x99rCaCKcAKuB86pqY1XdAJwK3A6cPKb+\ndODSqlpfVTdW1RpgC10YGa55TVW9r6o+DbwYeBjwS4u2FpoTv31I2tfs7Pea4WXXTDycJDkAWAlc\nPtNWVQVcBqwaM9uqfvqwTTP1SR4BLJ+1zG8AH9/BMrUI/EGUpO8zvOyaiYcT4CBgP2DbrPZtdAFj\nlOU7qT8EqDkuU/MwGAy48cYbd/iD5t4RSdo1hpfO/pPuwA6ELmAsZP2Oah4IcMUVV4yceOutt/Lu\nd7+b5z3veTz0oQ9d0OmLuezF/uytW7fy2c9+lg996EMcfvjhI6cP/7ur0xZ7up/tZ/vZfvbe+tnr\n1q3j0Y9+9H1+507yb8lQfx84suNzlO4IyuT0h3VuB55fVe8dar8AWFZVzx0xzxeBs6vqnKG2tcCJ\nVXVUkp8A/gU4sqquG6r5f8Cnqmr1iGX+KvCOhVovSZKWoF+rqgt3dyET33NSVXcmuRo4BngvQJL0\n788ZM9vmEdOP7dupqpuT3NLXXNcv84eApwBvGbPMTcCvAV8AvjP/NZIkacl5IPDjdH9Ld9vE95wA\nJHkh8Fa6y36vort655eBx1XVrUk2Al+uqjP7+lXAR4BXA+8Hpvr/r6iq6/uaVwKvAn6DLnC8BvhJ\n4Cer6rt7bOUkSdKcTHzPCUBVXdTf0+QsupNZrwGOq6pb+5LDgLuG6jcnmQJe279uojukc/1QzRuS\nHAicBzwE+ChwvMFEkqS2NbHnRJIkaUYLlxJLkiR9j+FEkiQ1ZcmFkyS/n+SqJN9Isi3J3yV5zKya\nByR5S5Lbknwzyd8mOXhSfZ60JKf2D1fc3r+uSPLsoemO10702909SdYPtTlusyRZ04/T8Ov6oemO\n2QhJHpbkbf243N7/vK6YVeODUIf0D5qdva3dk+TP++lua7MkuV+S1yT5fL8dfS7JH46o2+1tbcmF\nE+C/AX9Od1nxfwcOAD6Y5AeGav4MeA7wfODpdM/kuXgP97MlX6K78mll//oQ8J4kM3cAcrx2oH/C\n9il0D7Qc5riN9mm6E+OX96+nDU1zzGZJ8hDgY8B/AscBhwO/B3xtqMYHod7XE/n+Nrac7nYUBVzU\nT3dbu69X021Dvw08Dngl8Mok33uu3YJta1W1pF90t8+/B3ha//6H6H7InztU89i+5smT7m8rL+A/\ngN90vHY6Tg8GbgSOBj4MrO/bHbfR47UG2DJmmmM2elxeD3xkJzVfAVbPGss7gBdOuv+tvOjCyGeH\nxsdt7b5j9D7g/FltfwtsHHq/INvaUtxzMttD6NLyV/v3K+kusR5+aOCNwL/iQwNnduv9CnAg3U3v\nHK8dewvwvqr60Kz2J+K4jfPoJP+W5F+SvD3Jj/btbmuj/SLwySQX9YeqtyR56czE/o7ZPgh1B/o7\nlf8a8Fd9kz+fo10BHJPk0QBJngA8FfhA/37BtrUm7nMyKf2daP8M+Kf6/j1SlgPf7Qd02JJ+aGCS\nx9OFkQcC36T7RnFDkqNwvEbqQ9xRdH9UZzsEx22UK+lunHgjcCiwFvjHfvvzZ3O0RwC/BZxNd9+n\npwDnJPlOVb2dbmx8EOqOPRdYRnczUPDnc5zX0+0JuSHJ3XSnhvxBVb2zn75g29qSDifAucAR3PuY\n9jhzfRDhvuYG4Al0e5qeD2xM8vQd1C/p8UpyGF3wPbaq7pzLrCzhcauq4VtffzrJVcAXgRcy/rES\nS3rM6P5AXFVVf9S/vzbJT9IFlrfvYL6lPm7DTgYurapbdlK31MfsJOBXgV8BrgeOBN6U5CtV9bYd\nzDfncVuyh3WSvBn4eeCZVfWVoUm3APfvn8Uz7GDumwaXjKq6q6o+X1VbquoP6E7uPB3Ha5yVwEOB\nq5PcmeRO4BnA6Um+Szc2D3DcdqyqtgOfBR6F29o4A2D2I2y3Ag/v/38L3R+HQ2bVLPVxAyDJw+ku\njjh/qNltbbQ3AH9SVe+qqs9U1TuADcDv99MXbFtbkuGkDyYnAj9XVf86a/LVdLfKP2ao/jF0P+ib\n91gn23c/4AE4XuNcBvwU3TeLJ/SvT9J9k535/504bjuU5MHAI+lOsnNbG+1jdCdrDnss3R4nqupm\nuj8aw+M28yDUK/ZQH1t2Mt0fzg8MtbmtjXYg990Dcg99lljIbW3JHdZJci7dgwJPAL6dZCbhba+q\n71TVN5L8FbA+ydfozq84B/hYVV01mV5PVpLXApfSXVL8g3Qnjj0DeJbjNVpVfZtut+f3JPk28B9V\ntbV/77jNkuSNdFcEfBH4EWAd3R+Jd7qtjbUB+FiS36e7DPYpwEvpLl+f8WfAHyb5HN9/EOqXgffs\n2a62pT/v8DeAC6rqnpl2t7Wx3gf8QZIvAZ8BVtA9qPcvh2oWZlub9KVJE7gU6h7g7hGvFw/VPIDu\nXii30W2U7wIOnnTfJzhmfwl8nu5ysFuADwJHO15zHscP0V9K7LiNHaPp/hfZHXRXRlwI/IRjttNx\n+3ngOuD2/o/GySNq1tLtgbqd7rH2j5p0vyf9oru3yd2jxsJtbeR4PQhYD9xMd/+Sm+i+QOw/q263\ntzUf/CdJkpqyJM85kSRJ7TKcSJKkphhOJElSUwwnkiSpKYYTSZLUFMOJJElqiuFEkiQ1xXAiSZKa\nYjiRJElNMZxIkqSmGE4kzVuSC5Lck+SVs9pPTHLPuPkkaUcMJ5J2R9E9pO9VSZaNmCZJc2Y4kbS7\nLqN7WvWZ85k5yVOTfDjJt5N8NcmlM0Enyf2TnJNkW5I7knw0yROH5n1Gv+fm6CSf6JfxsSSPGapZ\nk+RTSV6U5OYkX08yneRBQzW/nOS6JLcnuS3JB5P8wPyHRNLuMJxI2l130wWT30nysLnMmORIunDz\naeBngKcC7wP260veCDwX+B/AUcDngE1JHjJrUX8MrAZWAncBfzVr+iOBE4GfB54DPAN4dd+H5cCF\nwF8Cj+unvRvIXNZF0sLZf9IdkLT3q6r3JLkGWAecModZ/xfwiar6naG2rQBJDgROBV5cVR/s204B\njgVeApw98/HAmVX1T33N64H/m+T+VfXdvibAr1fV7X3N24BjgD8CDqULQ39XVV/q6z8zh3WQtMDc\ncyJpobwK+PUkj5vDPEcCl4+Z9ki6L1BXzDRU1V3AVcDhs2r/eej/g/7fg4favjATTIZqZqZf2/fh\n00kuSvLSEXtmJO1BhhNJC6KqPgpsAv5kDrPdsYNpM4dVZp9YmxFtdw53pf/3fmOmz9TcD6Cq7qmq\nZwHPptuvJL6tAAABRElEQVRj8jvADUl+bMddl7RYDCeSFtKZwC8CP7uL9dfRHV4Z5XN0oeJpMw1J\n9geeCFy/G30cqao2V9U6unNb7qQ710XSBHjOiaQFU1X/nOQddHsfdsWfANcleQvwF3Sh4JnARVX1\n1ST/B3hjkq8BXwJeCfwA8NdDyxh14uoun8ya5Ml0AemDwL/TnZh7EIsQgCTtGsOJpIX2R8AL2YX7\nnFTVTUmeBbwO+DjdYZ6P0109A90VNQE2Aj8IfBJ4VlVtH17MqEXPob/fAJ4OnA78EPBF4IyZk3Al\n7Xmp8j5JkiSpHZ5zIkmSmmI4kbRoknwgyTdHvL6R5NWT7p+kNnlYR9KiSXIo3Qmso3y1qr6+J/sj\nae9gOJEkSU3xsI4kSWqK4USSJDXFcCJJkppiOJEkSU0xnEiSpKYYTiRJUlMMJ5IkqSn/H7V7A+gk\nVtzuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xb5a8e02c>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = 0.05\n",
    "\n",
    "N_conns_dist = sp.binom(post_slice_size, p)\n",
    "\n",
    "# A range of values for which the distribution has non-negligible probability mass\n",
    "N_conns_range = np.arange(20,80)\n",
    "ps = N_conns_dist.pmf(N_conns_range)\n",
    "plt.bar(N_conns_range, ps, width=0.01)\n",
    "plt.xlabel(\"N_conns\"); plt.ylabel(\"p\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scipy distribution objects provide a ppf or *percent point function* method (which we will also refer to as the *quantile function* $Q(p)$). For a continuous random variable, the ppf is the inverse of the cumulative distribution function (cdf). For a discrete random variable, the ppf gives the lowest value $k$ such that the cumulative probability of $k$, i.e., $P(X\\le k)$ exceeds the specified probability.\n",
    "\n",
    "The ppf function will allow us to compute a value such that a sample from the above distribution will be less than that value with probability 0.9999. Looking at the plot of the probability mass function above, a value around 80 seems plausible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Each row will have fewer than 79 synapses with probability 0.9999\n"
     ]
    }
   ],
   "source": [
    "single_row_extreme_value = int(N_conns_dist.ppf(quantile))\n",
    "\n",
    "print(\"Each row will have fewer than {} synapses with probability {}\".format(single_row_extreme_value, quantile))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantiles of the maximum number of synapses in a sub-row\n",
    "\n",
    "Note that *each sub-row* of our synaptic matrix has a 0.0001 probability of having more than 79 synapses. Our synaptic matrix contains $\\frac{5000^2}{1024} = 24414$ such sub-rows, and the probability that at least one of them will contain more than 79 synapses is 0.91. Worse, this probability increases with the size of the synaptic matrix.\n",
    "\n",
    "What we should really be calculating is the 0.9999 quantile on the *maximum number of synapses contained within a sub-row of the whole matrix*. We want the value such that, with probability 0.9999, no sub-row within the matrix has a number of synapses that exceeds that value.\n",
    "\n",
    "There is a simple expression for the cumulative distribution function of the maximum of $n$ identically and independently distributed random variables.\n",
    "\n",
    "If the cdf of the number of synapses in a single subrow is $F_X(x) = P(X \\le x)$, then the cdf of the maximum of $n$ identically and independently distributed random variables is\n",
    "\n",
    "$F_n(x) = P[(X_1 \\le x) \\cap (X_2 \\le x) \\ldots \\cap (X_n \\le x)] = (F_X(x))^n \\quad . $\n",
    "\n",
    "The quantile function $Q(p)$ is given by the inverse of the cdf, i.e., if $F_X(x) = p$ then $Q_X(p) = x$.\n",
    "\n",
    "If the cdf of the maximum $F_n^\\texttt{max}(x) = (F_X(x))^n = p$, then $F_X(x) = \\sqrt[n]{p}$, and so the quantile function of the maximum of $n$ samples is\n",
    "\n",
    "$Q_n^\\texttt{max}(p) = Q_X(\\sqrt[n]{p}) \\quad .$\n",
    "\n",
    "### Coin flipping example\n",
    "\n",
    "The reader might be reassured by an example. If we flip a coin $100$ times, the number of heads we will observe is binomially distributed. The median observed value, calculated using the scipy ppf function and by sampling, is as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The median number of heads when flipping 100 coins, exact: 50, sample: 50\n"
     ]
    }
   ],
   "source": [
    "coin_distribution = sp.binom(n=100, p=0.5)\n",
    "\n",
    "exact_median = int(coin_distribution.ppf(0.5))\n",
    "\n",
    "samples = coin_distribution.rvs(1000)\n",
    "sample_median = int(np.percentile(samples, 50))\n",
    "\n",
    "print(\"The median number of heads when flipping 100 coins, exact: {}, sample: {}\".format(exact_median, sample_median))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "What if we flip 100 coins 10 times, and record, amongst those 10 trials, the maximum number of heads observed?\n",
    "\n",
    "The following code calculates the median value of the maximum of 10 trials by sampling, and using the formula above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The median value of the maximum number of heads observed amongst 10 trials of 100 coin flips, exact: 57, sample: 58\n"
     ]
    }
   ],
   "source": [
    "samples = coin_distribution.rvs((10, 1000))\n",
    "samples_of_maximum = samples.max(axis=0)\n",
    "sample_median_of_maximum = int(np.percentile(samples_of_maximum, 50))\n",
    "\n",
    "maximum_median_quantile = 0.5 ** (1/10.0)\n",
    "exact_median_of_maximum = int(coin_distribution.ppf(maximum_median_quantile))\n",
    "\n",
    "print(\"The median value of the maximum number of heads observed amongst 10 trials of 100 coin flips, exact: {}, sample: {}\".format(exact_median_of_maximum, sample_median_of_maximum))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have good reason to want to be able to calculate quantiles without sampling. In the above coin flipping examples we were able to obtain a good estimate of the median by sampling 1000 values from the distribution. However, in order to obtain good estimates of extreme quantiles, we would need a very large number of samples.\n",
    "\n",
    "We can use the above to calculate the 0.9999 quantile of the maximum number of synapses appearing in a sub-row in the whole matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With probability 0.9999, the maximum number of synapses within a sub-row of this matrix will be less than or equal to 96.\n",
      "With probability 0.5, the maximum number of synapses within a sub-row of this matrix will be less than or equal to 81.\n",
      "With probability 0.9999, the number of synapses within a given sub-row of this matrix will be less than or equal to 79.\n"
     ]
    }
   ],
   "source": [
    "max_full_matrix_quantile = 0.9999 ** (post_slice_size / full_matrix_size)\n",
    "extreme_max_row_synapses = int(N_conns_dist.ppf(max_full_matrix_quantile))\n",
    "\n",
    "max_full_matrix_median = 0.5 ** (post_slice_size / full_matrix_size)\n",
    "median_max_row_synapses = int(N_conns_dist.ppf(max_full_matrix_median))\n",
    "\n",
    "print(\"With probability 0.9999, the maximum number of synapses within a sub-row of this matrix will be less than or equal to {}.\".format(extreme_max_row_synapses))\n",
    "print(\"With probability 0.5, the maximum number of synapses within a sub-row of this matrix will be less than or equal to {}.\".format(median_max_row_synapses))\n",
    "print(\"With probability 0.9999, the number of synapses within a given sub-row of this matrix will be less than or equal to {}.\".format(single_row_extreme_value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from the above that if we allocate memory per sub-row such that we allow for up to 96 synapses per sub-row, we will with probability 0.9999 have allocated enough memory for each sub-row and in the median case we will be \"wasting\" 15 synapses worth of memory in the sense that, in the median case, the maximum number of synapses in a sub-row of our matrix will be 81."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantiles of mixture distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some of the quantities that we need to estimate reasonable upper bounds for are random variables whose parameters are themselves random variables.\n",
    "\n",
    "For example, suppose that the delays associated with the synapses in our projection are gamma distributed with a scale parameter of 10. (We will ignore here that delays are quantised, and are multiples of the simulation time-step.)\n",
    "\n",
    "One quantity that we need to calculate is an extreme quantile (e.g., 0.9999) of the maximum number of synapses per sub-row whose delay is less than 8 * the time-step of the simulation.\n",
    "\n",
    "If the simulation time-step is 1 ms, then what we want is a value such that with probability 0.9999 the number of synapses in a sub-row whose delay is less than or equal to 7 ms does not exceed that value for any sub-row.\n",
    "\n",
    "Below is shown our example delay distribution, and the shaded region represents the probability that a given synapse will have a delay of less than or equal to 7 ms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The probability that a given synapse has a delay less than 7 is 0.17\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhIAAAFkCAYAAAB1rtL+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl4lNX9/vH3JyFsotGKoEjrgoLyrQvEpbgLWtRWsXbR\nuKCCiLtF/Wlta7XaulQBcUFxl6KpKEpFhbApuKAoAaqAG6CiAoJL2CHMnN8fZyIhJCQzmcyZ5X5d\n13NBnjnzzJ1xJJ+c5yzmnENEREQkEXmhA4iIiEjmUiEhIiIiCVMhISIiIglTISEiIiIJUyEhIiIi\nCVMhISIiIglTISEiIiIJUyEhIiIiCVMhISIiIglTISEiIiIJS6iQMLNLzWyhma01s7fN7OCttO1s\nZs/F2kfN7Io6rn19rN2gRLKJiIhI6sRdSJjZ6cBA4EagCzAbKDWz1rU8pSUwH7gOWFzHtQ8G+sWu\nKSIiImkukR6JAcAw59xw59yHwEXAGqBPTY2dc+85565zzo0ENtR2UTNrBYwALgB+SCCXiIiIpFhc\nhYSZFQBFwKTKc85vHzoR6NbALPcDY5xzkxt4HREREUmRJnG2bw3kA0urnV8KdEo0hJmdgb9NUlTP\n9jsCPYHPgHWJvq6IiEgOag7sDpQ6575t6MXiLSRqY4BL6Ilm7YG7geOdcxX1fFpP4KlEXk9EREQA\nOAt4uqEXibeQWA5EgLbVzrdhy16K+ioCdgJmmJnFzuUDR5nZZUCz2O2Tqj4DGDFiBPvuu2+CLyvx\nGjBgAIMHDw4dI6foPU89veepp/c8tebNm8fZZ58NsZ+lDRVXIeGcqzCzGUAP4EWA2A//HsA9CWaY\nCOxX7dwTwDzg9hqKCIjdzth3333p2rVrgi8r8SosLNT7nWJ6z1NP73nq6T0PJilDAxK5tTEIeDJW\nUEzHz+Joif/hj5kNB750zv059nUB0Bl/+6MpsKuZHQCscs7Nd86tBuZWfQEzWw1865ybl9B3JSIi\nIikRdyHhnBsZWzPiZvwtjllAT+fcsliT9sDGKk9pB8xk0xiKa2LHFKB7bS8Tby4RERFJvYQGWzrn\nhgJDa3mse7WvPyfOaabVryEiIiLpSXttSL0VFxeHjpBz9J6nnt7z1NN7ntms5rGM6c3MugIzZsyY\noQE6IiIicSgrK6OoqAigyDlX1tDrqUdCREREEqZCQkRERBKmQkJEREQSpkJCREREEqZCQkRERBKm\nQkJEREQSpkJCREREEqZCQkRERBKmQkJEREQSpkJCREREEpbQpl0issk338ADD8Bnn8EPP2w6nINe\nveCcc2CvvUKnFBFpHOqREEnQ99/DX/4Ce+4JAwfCRx/B+vXQrh0cdhgccAAMHgx77w2HHw7DhvkC\nQ0Qkm6hHQiROK1fC3Xf74qGiAq64Aq65Bnbcccu2DzwAL74Iw4fDJZfA3/4G//kPHHts6nOLiDQG\n9UiIxGHxYjjkEPjnP+H882HBArjttpqLCICWLeGMM+CVV+Dzz2G//eC44/xzotHUZhcRaQzqkRCp\np6++gu7dYc0a+N//oGPH+J7fvj2UlsJNN8Gf/wxvveV7KnbYoVHiioikhHokROph0SI4+mhYtw6m\nTIm/iKiUnw+33AIvvwxvvgldu8L77yc3q4hIKqmQEKnD55/7ImLjRl9E7Llnw6950klQVgaFhdCz\np38NEZFMpEJCZCsWLYJjjvF/nzIFdt89edfefXd/q6N5czjxRD8LREQk06iQEKlFNOrXgIhEfBGx\n227Jf422bWHsWFi6FE491d86ERHJJCokRGoxZIgvIIYPh5/+tPFep1MnGDMGpk+Hc8/VbA4RySwq\nJERqMHcuXH89/PGPm25tNKbDDoOnn4Znn4Vrr2381xMRSRYVEiLVVFRA796wxx5w662pe93f/Mb3\nggwc6BetEhHJBFpHQqSaW2+FWbNg2jRo0SK1r3355X5a6CWXwFFH+eW2RUTSmXokRKp47z34xz/8\nglEHHxwmw/33Q7Nm0Lev3/hLRCSdqZAQidmwwd/S2H9/+Otfw+XYcUd49FEYNw4eeihcDhGR+lAh\nIRLzyCPw4Yfw+OPQtGnYLCedBBdeCFdfDfPnh80iIrI1KiRE8Ptn/OMfcNZZvkciHQwc6NeZOPdc\nv5aFiEg6UiEhAgwdCsuW+Q210kWrVvDEE35zr7vuCp1GRKRmKiQk561YAbffDn36QIcOodNs7sgj\n4aqrfIGj/ThEJB2pkJCcN2QIrFwZdoDl1tx0k99qXAtViUg6UiEhOe277/xtg4svbtxlsBuiVSu4\n7TYYORJefz10GhGRzSVUSJjZpWa20MzWmtnbZlbrjHsz62xmz8XaR83sihraXG9m081shZktNbMX\nzKxjItlE4nHnnX578OuvD51k6845x69r8cc/ai8OEUkvcRcSZnY6MBC4EegCzAZKzax1LU9pCcwH\nrgMW19LmSOBe4FDgOKAAGG9mKV5XUHLJ0qVwzz1wxRV+dkQ6y8uDu++GsjI/AFNEJF0k0iMxABjm\nnBvunPsQuAhYA/SpqbFz7j3n3HXOuZHAhlranOSc+7dzbp5z7n3gPOBnQFEC+UTq5fbboUkT+H//\nL3SS+jnsMCgu9qturlgROo2IiBdXIWFmBfgf7pMqzznnHDAR6JbEXNsDDvguidcU+dEPP/hVI6+8\nEn7yk9Bp6u+OO3wRkcrNxEREtibeHonWQD6wtNr5pcDOyQhkZgbcDbzhnJubjGuKVPfkk35J7Isu\nCp0kPj/9qZ+9MXiwVrwUkfSQrFkbhu9BSIahQGfgjCRdT2Qz0ahfgOq00zJzd81rr4U2beAvfwmd\nREQk/m3ElwMRoPrQtDZs2UsRNzO7DzgJONI5V9vAzB8NGDCAwsLCzc4VFxdTXFzc0CiSxSZNgo8/\n9ntrZKKWLeGGG3xvyg03wP/9X+hEIpKuSkpKKCkp2exceXl5Ul/DXJz7FJvZ28A7zrkrY18b8AVw\nj3PuzjqeuxAY7Jy7p4bH7gN6AUc75xbUcZ2uwIwZM2bQtWvXuPKL9OoFCxfC7NlgFjpNYjZsgI4d\n4ZBD/PoSIiL1VVZWRlFREUCRc66soddL5NbGIOBCM+ttZvsAD+KneD4BYGbDzezHoWBmVmBmB5jZ\ngUBTYNfY1x2qtBkKnAWcCaw2s7axo3nC35lIDT77DF56CS67LHOLCPC7k95wAzz7LLz/fug0IpLL\n4i4kYtM4rwZuBmYC+wM9nXPLYk3as/nAy3axdjNi568ByoCHq7S5CNgOeA34usrxh3jziWzNgw/C\nttv6XT4zXe/esOee8Pe/h04iIrks3jESADjnhuIHRdb0WPdqX39OHQWLc05LdUujW7fOj4s4/3zY\nZpvQaRquoMDvD9Knj79Nc8ABoROJSC7SD3DJGc88A99+C5dcEjpJ8pxzjt+xNJ22PxeR3KJCQnLG\n/fdDz56w996hkyRPkyZ+rMTo0TBzZug0IpKLVEhITpg+Hd591w+yzDZnneWLI/VKiEgIKiQkJzzy\nCPzsZ3DiiaGTJF9lr8SLL/pNvUREUkmFhGS99ev9NMmzz4b8/NBpGseZZ/oZHHdudSUXEZHkUyEh\nWW/sWL9JVzZM+axNfj5cdZVfnGrhwtBpRCSXqJCQrDdiBBx4IHTuHDpJ4zr/fNhhB7+hl4hIqqiQ\nkKxWXu5Xsjz77NBJGl/Lln4w6aOP+mmuIiKpoEJCstqoUX5fijNyZC/ZSy/dtLupiEgqqJCQrDZi\nBBx7LOy6a+gkqbHTTn6ly3vvhbVrQ6cRkVygQkKy1ldfwWuv5cZtjaquusrf2njyydBJRCQXqJCQ\nrFVS4nfJPO200ElSq0MH/z0PHAiRSOg0IpLtVEhI1nrqKTj5ZCgsDJ0k9a69Fj79FP7739BJRCTb\nqZCQrDRnDsyalXu3NSodfDAcfTT861+hk4hItlMhIVnpqaf8mgrZuCR2fV1zDbzzjj9ERBqLCgnJ\nOtEoPP00/P73foxErjrpJD9e4p57QicRkWymQkKyzjvvwOef+/0ncllenl9X4tlnYfHi0GlEJFup\nkJCs88IL0KYNHHFE6CThnX++75UZNix0EhHJViokJKs45wuJU07J3p0+47H99tC7Nzz4oF/hU0Qk\n2VRISFaZN89Pe/zNb0InSR+XXQZLl/pbHCIiyaZCQrLKCy9Aq1bQvXvoJOmjc2c47ji/bLaISLKp\nkJCsMnq0n/LZvHnoJOnl8sv9INTp00MnEZFso0JCssaiRfDee7qtUZNf/Qr22EO9EiKSfCokJGv8\n979QUODXT5DN5ef7qaDPPANLloROIyLZRIWEZI3Ro/3YiFzcW6M++vTxhdZDD4VOIiLZRIWEZIXv\nv/dbhp96augk6WuHHfzeIw89BBs3hk4jItlChYRkhZde8ltmn3JK6CTp7aKL4Kuv4OWXQycRkWyh\nQkKywujR8ItfQLt2oZOkty5d4JBD/AJVIiLJoEJCMt7atTBunG5r1NfFF0NpKSxYEDqJiGQDFRKS\n8SZMgDVrNO2zvv7wBz8gVYMuRSQZVEhIxhs9GvbdFzp2DJ0kM7RsCeeeC489BuvXh04jIplOhYRk\ntGjUD7Ts1St0kszSvz8sW+aXFBcRaQgVEpLRZszwPxB/9avQSTLLvvvC0Udr0KWINJwKCcloY8f6\n+/2/+EXoJJnnootgyhSYOzd0EhHJZAkVEmZ2qZktNLO1Zva2mR28lbadzey5WPuomV3R0GuKVBo7\nFo4/Hpo0CZ0k85x2Guy0EwwbFjqJiGSyuAsJMzsdGAjcCHQBZgOlZta6lqe0BOYD1wGLk3RNEb79\n1u9oeeKJoZNkpqZNoW9fePJJP+tFRCQRifRIDACGOeeGO+c+BC4C1gB9amrsnHvPOXedc24ksCEZ\n1xQBGD8enIMTTgidJHNdeCGsWOE38xIRSURchYSZFQBFwKTKc845B0wEuiUSoDGuKblh7Fg44ACt\nZtkQe+wBxx0HDz8cOomIZKp4eyRaA/nA0mrnlwI7J5ihMa4pWS4a9atZ6rZGw/XrB9OmwZw5oZOI\nSCZK1hA1A1ySrlXvaw4YMIDCantGFxcXU1xcnOQokm7Kyvy0TxUSDderlx90+fDDcPfdodOISDKV\nlJRQUlKy2bny8vKkvka8hcRyIAK0rXa+DVv2KDT6NQcPHkzXrl0TfFnJZGPHwnbbQTfd/Gqwpk03\nrXR5++3QvHnoRCKSLDX9cl1WVkZRUVHSXiOuWxvOuQpgBtCj8pyZWezrtxIJ0BjXlOw3dqy/t19Q\nEDpJdrjgAvjuO610KSLxS2TWxiDgQjPrbWb7AA/ip3g+AWBmw83s1srGZlZgZgeY2YFAU2DX2Ncd\n6ntNkaq++07TPpOtUyc46igNuhSR+MU9RsI5NzK2vsPN+NsRs4CezrllsSbtgY1VntIOmMmm8Q7X\nxI4pQPd6XlPkRxMm+MGWmvaZXP36wTnnwKefwl57hU4jIpkioZUtnXNDnXO7O+daOOe6Oefeq/JY\nd+dcnypff+6cy3PO5Vc7utf3miJVjR0L++0H7duHTpJdfvtb2H57eOSR0ElEJJNorw3JKJr22Xha\ntPA9Ek88ARUVodOISKZQISEZZdYsWLpUhURj6dfPv79jxoROIiKZQoWEZJSxY2GbbaIcemgkdJSs\ntN9+cOihur0hIvWnQkIyyksvrWP16pd4443JoaNkrX79/O2jL74InUREMoEKCckYq1fDu+8WABNY\ntWpV6DhZ6w9/gJYt/VgJEZG6qJCQjPHGGxCJ5OP3c5PGsu22cPrp8PjjfnCriMjWqJCQjDFhgsNs\nMfBh6ChZr29f+OwzmKw7SCJSBxUSkjFefnkdzo0PHSMndOsG++wDjz4aOomIpDsVEpIRli2DDz9s\ngW5rpIaZ75V44QW/JLmISG1USEhG2NTFnugmsxKv3r0hEoGnngqdRETSmQoJyQgTJkQx+xA4MHSU\nnNGmDZx8sl9Twrm624tIblIhIRnhlVc2xMZHHBc6Sk7p2xf+9z+YMSN0EhFJVyokJO0tWACLFzcn\nL+9toEvoODmlZ09o106DLkWkdiokJO1NnAiwkWi0CZAfOE1uadIEzjsPnn4a1qwJnUZE0pEKCUl7\n48ZVAO8Ch4WOkpP69IEVK2DUqNBJRCQdqZCQtBaNwqRJDpgA9AgdJyd16ADHHqvbGyJSMxUSktZm\nz4YVK5qSn/8BsFfoODmrb1+YMgU+/TR0EhFJNyokJK1NnAhma4hEdgAsdJycddppUFjo998QEalK\nhYSktZdfXo9zU4BjQkfJaS1awJln+h1BI5HQaUQknaiQkLS1fj1Mm1a522f30HFyXp8+8PXXUFoa\nOomIpBMVEpK2pk2DDRuakJ+/CGgbOk7OKyqC/feHxx4LnURE0okKCUlbkyc7zL4jEmkfOorgN/Lq\n0wdefNFvoiYiAiokJI2NG7ce515F4yPSx9ln+4JixIjQSUQkXaiQkLS0Zg3MnFkATAGOCh1HYnbc\nEXr18mtKaCMvEQEVEpKmpk2DjRvzyc9fCmwfOo5U0bcvzJkD774bOomIpAMVEpKW/PiIb4lEfhY6\nilRz3HHw059q0KWIeCokJC2Vlq7DuclofET6yc/3G3mVlGgjLxFRISFpaPVqmDWrKfAacETgNFKT\n88/3G3k991zoJCISmgoJSTtvvQWRSD75+cuBwtBxpAZ77AHdu2sjLxFRISFpyI+PWEYkslvoKLIV\nffvC1KnwySehk4hISCokJO1sGh9xbOgoshW/+Q1sv70GXYrkOhUSklZWrYLZs5vi14/Q+Ih01qIF\nnHUWPPkkbNwYOo2IhKJCQtLKG29ANJpPXt63wLah40gd+vaFxYth7NjQSUQklIQKCTO71MwWmtla\nM3vbzA6uo/3vzWxerP1sMzux2uPbmNl9ZrbIzNaY2Rwz659INslsr77qMFtKNLpn6ChSD126+EOD\nLkVyV9yFhJmdDgwEbgS6ALOBUjNrXUv7bsDTwMPAgcBoYLSZda7SbDDwS+BMYB/gbuA+M/t1vPkk\ns40duw7nJqH1IzJH377w0kuwZEnoJCISQiI9EgOAYc654c65D4GLgDVAn1raXwmMdc4Ncs595Jy7\nESgDLqvSphvwpHPudefcF865h/EFyiEJ5JMMtWIFzJnTDHgdODx0HKmnM8+EJk1g+PDQSUQkhLgK\nCTMrAIqASZXnnHMOmIgvBmrSLfZ4VaXV2r8FnGJm7WKvcyywd6yd5Ag/PiIvNj6iVeg4Uk877AC/\n/a028hLJVfH2SLQG8oGl1c4vBXau5Tk716P95cA84Esz2wC8AlzqnHszznySwfz6EYuJRvcOHUXi\n1LcvfPwxvKn/Y0VyTpMkXceAeH4Xqd7+CuBQ4NfAF/h9o4ea2dfOLyhQowEDBlBYuPnKh8XFxRQX\nF8cRRdLFuHGV4yO0fkSmOeYYv9rlo4/CEZq1K5I2SkpKKCkp2exceXl5Ul8j3kJiORAB2lY734Yt\nex0qLdlaezNrDvwT6OWcGxd7/AMz6wJcA9RaSAwePJiuXbvG9Q1Ieiovh7lzK8dHDAodR+KUl+f3\n37j9dhgyBLbbLnQiEYGaf7kuKyujqKgoaa8R160N51wFMAPoUXnOzCz29Vu1PG1a1fYxx8fOAxTE\njuo9GpF480nmevNNcC6PvLxyYJvQcSQB558P69bBf/4TOomIpFIiP6gHAReaWW8z2wd4EGgJPAFg\nZsPN7NYq7YcAJ5rZVWbWycxuwg/YvA/AObcSv4zhnWZ2tJntbmbnAb2B5xP7tiTTvPaaw2wJ0WiH\n0FEkQe3bwwknwCOPhE4iIqkUdyHhnBsJXA3cDMwE9gd6OueWxZq0p8pASufcNKAYuBCYBZyGv40x\nt8plTwfeBUYAc4Brgeudcw/Fm08yk99f41X88BjJVP36wbvvwuzZoZOISKokNNjSOTcUGFrLY91r\nODcKGLWV630D9E0ki2S+Vavggw8qx0fcETqONMCvfgVt2/peiXvvDZ1GRFJBYxAkuLfe8utH5Odr\nf41MV1Dgx0qMGAFr14ZOIyKpoEJCgpsyxWG2nEhkt9BRJAn69IEffoBRtfZBikg2USEhwY0fvw7n\nXgOODh1FkmDvvf26Ehp0KZIbVEhIUGvXwqxZTYGpgFYyyhb9+sGUKX61SxHJbiokJKi334aNG/PJ\nz18KFNbZXjLDaaf5PTi0vbhI9lMhIUFNmQJm3xOJtA8dRZKoeXM45xx44gnYsCF0GhFpTCokJCi/\nfsQUtH5E9rngAvjmG3jppdBJRKQxqZCQYNavhxkzmuAXNtX4iGyz335w6KHw8MOhk4hIY1IhIcG8\n+y5UVDQhP/9rYMfQcaQR9OsHpaXw2Wehk4hIY1EhIcH48REriESqbw4r2eL006FVKw26FMlmKiQk\nGL9+xFQ0PiJ7tWoFZ5/tC4mKitBpRKQxqJCQICoq4J138vHrR6iQyGb9+8PixRp0KZKtVEhIEDNm\nwPr1BeTnfwG0CR1HGtEBB8Ahh8BD2stXJCupkJAg/PiI1UQiGmSZC/r316BLkWylQkKCmDBhPc69\nARwZOoqkwOmnw7bbaiqoSDZSISEpF4nAW2/l4deP0PiIXLDNNn7Q5WOPadClSLZRISEpN2sWrF1b\nQH7+AqBd6DiSIv37w5IlMGZM6CQikkwqJCTlpk4Fs3VEItuHjiIptP/+8ItfwLBhoZOISDKpkJCU\nmzBhA85NAw4LHUVSrH9/GD8eFi4MnUREkkWFhKRUNApvvAF+fMTRgdNIqv3hD1BYqKmgItlEhYSk\n1Jw5sHJlU/LyPgR2Cx1HUqxlSzj3XL/S5fr1odOISDKokJCU8utHVBCNtgodRQK5+GJYtgxGjQqd\nRESSQYWEpNSkSRtw7h3gF6GjSCD77APdu8PQoaGTiEgyqJCQlHEOXnvN4ffX0PiIXHbJJfDmm/C/\n/4VOIiINpUJCUuajj+CHH5qRl/c+sFfoOBLQKadAu3bwwAOhk4hIQ6mQkJSZOhVgI9FoM8ACp5GQ\nCgqgXz/4979hxYrQaUSkIVRISMpMmrQRmAkcGjqKpIF+/WDdOl9MiEjmUiEhKeEcTJ68Ee2vIZV2\n3RVOPdUPunQudBoRSZQKCUmJhQth+fLm5OXNBDqHjiNp4pJLYO7cytteIpKJVEhISkyZAhAlGs1D\n4yOk0rHHQqdOmgoqkslUSEhKvPrqRuB94KDQUSSNmPleieefh8WLQ6cRkUSokJCUmDixAngNjY+Q\n6nr3hmbNtCuoSKZSISGNbtEiWLy4BXl5M4D9Q8eRNLP99n7/jQcf1P4bIpkooULCzC41s4VmttbM\n3jazg+to/3szmxdrP9vMTqyhzb5m9l8z+8HMVpnZO2bWPpF8kl4qB9JFoxEgP2gWSU+XXQZLl8Kz\nz4ZOIiLxiruQMLPTgYHAjUAXYDZQamata2nfDXgaeBg4EBgNjDazzlXadABeB+bi+773A24B1sWb\nT9LP5MkRzObi//OLbGnffeGXv4QhQzQVVCTTJNIjMQAY5pwb7pz7ELgIWAP0qaX9lcBY59wg59xH\nzrkbgTLgsipt/gG87Jy73jn3P+fcQufcS8655QnkkzQzfvwGnHsVOCZ0FEljV1wB770H77wTOomI\nxCOuQsLMCoAiYFLlOeecAyYC3Wp5WrfY41WVVrY3MwN+BXxiZuPMbGnsdkmveLJJevrqK/jyyxaY\nTcd3YInU7MQToUMHuPfe0ElEJB7x9ki0xt/kXlrt/FJg51qes3Md7dsArYDrgFeA44EXgOfN7Mg4\n80ma8etHAFQATQImkXSXlweXXw4jR8LXX4dOIyL1laxZGwbEc2ezavvKDKOdc/fEbm3cAbyEv20i\nGezVVyOYfYhzGh8hdTvvPGjeXFNBRTJJvL8iLgciQNtq59uwZa9DpSV1tF8ObATmVWszDzh8a2EG\nDBhAYWHhZueKi4spLi7e2tMkhTaNjzg6dBTJAIWFvph48EH485/9+hIikriSkhJKSko2O1deXp7U\n14irkHDOVZjZDKAH8CL8OMahB3BPLU+bVsPjx8fOV17zXaBTted1BD7fWp7BgwfTtWvXeL4FSaHF\ni+GLL1pg9jbOXRA6jmSIyy6D++7ztzjOOSd0GpHMVtMv12VlZRQVFSXtNRK5tTEIuNDMepvZPsCD\nQEvgCQAzG25mt1ZpPwQ40cyuMrNOZnYTfsDmfVXa3AmcbmYXmFkHM7sM+DVwfwL5JE1Ujo9wrgIo\nCJpFMkenTnDCCZoKKpIp4i4knHMjgauBm4GZ+KUKezrnlsWatKfKwEvn3DSgGLgQmAWcBvRyzs2t\n0mY0fjzEtcD/8FNJT4s9VzKUHx/xMVrNUuJ15ZUwYwa88UboJCJSl4SG0TvnhgI17tfnnOtew7lR\nwKg6rvkEsV4NyQ6lpRtwbjJaP0Li1bMndO4MAwfCkZq7JZLWtNeGNIolS+Dzz1tg9g7+TpZI/ZnB\nVVfBiy/CJ5+ETiMiW6NCQhpF5f4afpVzjY+Q+J11Fuy0EwweHDqJiGyNCglpFJMnRzH7FOc0PkIS\n07y5n8HxxBPw7beh04hIbVRISKMYP35dbHyE1o+QxF18sf/zgQfC5hCR2qmQkKT75htYuLAlZtOA\ng0LHkQzWujWce65fV2Kd9gIWSUsqJCTpNu2vsRZoGjCJZIMBA3xx+vTToZOISE1USEjSvfpqFLP5\nOLdf6CiSBTp2hJNPhkGDtECVSDpSISFJV1q6XutHSFJdfTXMmQOlpaGTiEh1KiQkqZYtgwULWsTG\nRxwcOo5kiSOPhIMOgjvvDJ1ERKpTISFJ9eqrlX9bjcZHSLKYwXXXweTJMH166DQiUpUKCUmqiRP9\n/hrOHRg6imSZ3/zGj5e47bbQSUSkKhUSklTjxm3AuYnAFluuiDRIfr7vlRg9GubNC51GRCqpkJCk\nWbQIFi1qQV7eNLS/hjSGs8+GXXeFO+4InUREKqmQkKSpHB8RjUZJcGNZka1q2tTP4HjqKfjii9Bp\nRARUSEgSTZgQAWaj1SylMfXrB9ttB3fdFTqJiIAKCUkS52D8+ApgEhofIY2pVSu4/HJ45BE/3VhE\nwlIhIUlfAevJAAAcWUlEQVQxfz58801z8vLeA7SipTSuyy+HvDy4557QSUREhYQkxeTJABuJRpug\nj5U0th13hAsv9Jt5rVgROo1IbtO/+JIUpaUVwAygW+gokiOuugpWr4ahQ0MnEcltKiSkwZyDSZOi\naHyEpFL79tC3rx90uWpV6DQiuUuFhDTYnDlQXt6MvLzZQMfQcSSHXH+9v7Vx//2hk4jkLhUS0mCT\nJ4PZBqLRVoCFjiM55Gc/gz591CshEpIKCWkwvyz2m8CRoaNIDrr+eigv11gJkVBUSEiDRCIwZQrA\nZODYwGkkF+22m++VuPNOP/hSRFJLhYQ0yMyZsGZNU/LzPwJ2Cx1HcpR6JUTCUSEhDeLHR6whEtkx\ndBTJYbvtBuefr14JkRBUSEiDvPLKepybChwdOorkuOuvh++/hwceCJ1EJLeokJCErVsHb72VD0xA\n4yMktN13970S//qXZnCIpJIKCUnYG29ARUUT8vM/B9qGjiPCX/7ix0oMGRI6iUjuUCEhCRs3zpGX\nt4RIpEPoKCKAHytx8cW+V+Lbb0OnEckNKiQkYWPGrCMaHQf8MnQUkR/9+c8QjcIdd4ROIpIbVEhI\nQhYvho8/boHZa8DhoeOI/KhNG7+h1733wldfhU4jkv1USEhCJk6s/Ns6oHnAJCJbuvpq2GYbuPnm\n0ElEsl9ChYSZXWpmC81srZm9bWYH19H+92Y2L9Z+tpmduJW2w8wsamZXJJJNUmPs2I3ATJw7NHQU\nkS1st52/xfHoo/Dxx6HTiGS3uAsJMzsdGAjcCHQBZgOlZta6lvbdgKeBh4EDgdHAaDPrXEPbU4FD\nAHVIprFoFMaNiwClaHyEpKtLLoFddoG//S10EpHslkiPxABgmHNuuHPuQ+AiYA3Qp5b2VwJjnXOD\nnHMfOeduBMqAy6o2MrNdgXuAM4GNCeSSFHn/ffj++2bk5ZUBW9SDImmheXO46SZ45hkoKwudRiR7\nxVVImFkBUARMqjznnHPARKBbLU/rFnu8qtKq7c3MgOHAv5xz8+LJJKk3fjyYrSUaLUTbhks6O/dc\n6NQJ/vQncC50GpHsFG+PRGsgH1ha7fxSYOdanrNzPdr/CdjgnLsvzjwSwJgx63FuMtA9dBSRrWrS\nxE8DnTABxo4NnUYkOyVr1oYB8dT7P7Y3syLgCuD8JGWRRrRmDbz9dhP8stjHh44jUqdTToFjjvEz\nOSoqQqcRyT5N4my/HIiw5XrIbdiy16HSkjraHwHsBCzydzgA3+sxyMz+6Jzbs7YwAwYMoLCwcLNz\nxcXFFBcX1/FtSKKmToWKinzy878kEqlxfK1IWjGDwYOha1d46CG49NLQiURSp6SkhJKSks3OlZeX\nJ/U14ioknHMVZjYD6AG8CD+Ob+iBHyhZk2k1PH587Dz4sRETqj1nfOz841vLM3jwYLp27RrPtyAN\nVFrqMPuaSKRj6Cgi9XbggX5DrxtvhDPPhB12CJ1IJDVq+uW6rKyMoqKipL1GIrc2BgEXmllvM9sH\neBBoCTwBYGbDzezWKu2HACea2VVm1snMbsIP2LwPwDn3vXNubtUDqACWOOc+Sfg7k0YxZsw6nNOy\n2JJ5/vEPv2PtP/4ROolIdom7kHDOjQSuBm4GZgL7Az2dc8tiTdpTZSClc24aUAxcCMwCTgN6xQqG\nWl8m3lzS+L76CubPb4HZVGqfpCOSnnbZxc/euPde+PTT0GlEske8YyQAcM4NBYbW8tgWQ/mdc6OA\nUXFcv9ZxERLO+PEAUZyrAJoFTiMSv6uv9uMkrr0Wnn8+dBqR7KC9NqTeXnyxAngP+EXoKCIJadEC\nbr8dXngBXnstdBqR7KBCQuplwwYoLQU/xvZXgdOIJO6MM6BbN7jsMk0HFUkGFRJSL1Onwtq1BeTn\nzwE6hI4jkrC8PBg6FObNgyFDQqcRyXwqJKRexoxx5OV9pWmfkhUOPNCvJ3HTTfDll6HTiGQ2FRJS\nJ+dg1Kj1RKMvAr8OHUckKW65BVq1ggEDQicRyWwqJKROH30EX33VnLy819C0T8kWhYUwcCA891zl\njCQRSYQKCanTmDGVu322IMEZwyJp6cwz/T4cl17qF6sSkfipkJA6jRq1HucmACeEjiKSVGZw//3w\n2Wdw552h04hkJhUSslXffQfvvlsAvAL0DB1HJOk6d4arroJbb9WKlyKJUCEhW1VaCtFoHmbLAe10\nJNnpb3/zS2hfcAFEo6HTiGQWFRKyVS+8sBGzMpzTapaSvbbZBh5+GKZM8Utoi0j9qZCQWm3cCK+8\nEsW5McDJoeOINKoePXyPxLXXwqJFodOIZA4VElKrt96C1aubkp8/G9BCVJL97roLtt0W+vf366eI\nSN1USEit/GqWS4lEdgcsdByRRldYCA8+CGPHwogRodOIZAYVElKr555bTzQ6Bm3SJbnk5JOhuBj+\n+EdYujR0GpH0p0JCavTJJ/DZZ83Jy3sVODJ0HJGUuuceyM+HSy7RLQ6RuqiQkBo991zlapZ5QNPQ\ncURSqnVrv0Po88/D8OGh04ikNxUSUqMRI9bHZmv0Ch1FJIjf/Q7OPRcuuwwWLAidRiR9qZCQLcyf\nD3PnNsPsReCk0HFEgrnnHthpJzjnHD8dWkS2pEJCtvDss2C2DucMaBk6jkgw220H//43vP023HZb\n6DQi6UmFhGxh022NU0JHEQnu8MPhL3+Bv/8dpk8PnUYk/aiQkM0sWABz5jTD7L/otoaId8MNUFQE\nZ50Fq1aFTiOSXlRIyGY23dYA2CZ0HJG0UFAATz0Fixf7wZeaEiqyiQoJ2Yy/rfESuq0hsrm99oJh\nw+DJJ+HRR0OnEUkfKiTkRwsWwAcfVN7W0GqWItWddRZcdJHvlZg5M3QakfSgQkJ+5BehWodzDt3W\nEKnZ4MHw85/7dSZ++CF0GpHwVEjIj0aMWIdzL6Mtw0Vq17y5H0v03Xdw3nkaLyGiQkIAWLgQ3n+/\nuW5riNTDHnv4pbP/+1+/9bhILlMhIUDV2xoRoFXoOCJp7+ST4U9/guuvh8mTQ6cRCUeFhADw73+v\n120NkTjdcgv06OHHS3zySeg0ImGokBDmzIH332+G2XPotoZI/TVpAs88A23awK9/Dd9/HzqRSOqp\nkBAee8yRl/cdzjUHtg0dRySjbL89vPQSLFsGf/iDNveS3KNCIsdVVMDjj1cQjQ4HeoeOI5KR9toL\nRo2C116DAQNCpxFJLRUSOW7sWPj++6bk548Hjg4dRyRjHXss3HefP4YODZ1GJHUSKiTM7FIzW2hm\na83sbTM7uI72vzezebH2s83sxCqPNTGzO8zsf2a2ysy+MrMnzWyXRLJJfB5+eCNms4hEfoHqSpGG\n6d8frrwSLr8cRo8OnUYkNeL+yWFmpwMDgRuBLsBsoNTMWtfSvhvwNPAwcCAwGhhtZp1jTVrGzv89\ndr3fAJ2A/8abTeLzzTcwdmwezj0KnBs6jkhWGDgQTjsNzjgDpk4NnUak8SXyK+gAYJhzbrhz7kPg\nImAN0KeW9lcCY51zg5xzHznnbgTKgMsAnHMrnHM9nXOjnHOfOOemxx4rMrP2CeSTehoxAqLRjZh9\nCewWOo5IVsjP9/9vHX44nHIKzJ4dOpFI44qrkDCzAqAImFR5zvmNGSYC3Wp5WrfY41WVbqU9wPaA\nA7SSfSNxDoYNW49zo3Hud6HjiGSVZs3ghRegQwc44QS/IZ5Itoq3R6I1kA8srXZ+KbBzLc/ZOZ72\nZtYMuB142jm3Ks58Uk8zZsDHHzcjL28kcFroOCJZZ7vt/GDmVq2gZ09YWv1fQZEs0SRJ1zF8D0KD\n2ptZE+DZ2GOX1HWRAQMGUFhYuNm54uJiiouL44iSmx59NEpe3jdEozsBLULHEclKbdrA+PH+Nsfx\nx8OkSbDTTqFTSS4pKSmhpKRks3Pl5eVJfY14C4nlQARoW+18G7bsdai0pD7tqxQRPwW616c3YvDg\nwXTt2rUesaWqdetgxIgI0ehjwHmh44hktT32gIkToXt3v5z25MnQusah6SLJV9Mv12VlZRQVFSXt\nNeK6teGcqwBmAD0qz5mZxb5+q5anTavaPub42PnKa1QWEXsCPZxzWmi2EY0eDatWFZCf/zpwSOg4\nIlmvc2dfQHzzjS8oli8PnUgkeRKZtTEIuNDMepvZPsCD+CmcTwCY2XAzu7VK+yHAiWZ2lZl1MrOb\n8AM274u1zwdGAV2Bs4ECM2sbOwoS/L5kKwYPrsDsdSKR7vi7TCLS2FRMSLaKu5Bwzo0ErgZuBmYC\n+wM9nXPLYk3aU2UgpXNuGlAMXAjMwo/s6+Wcm1ul/a9jf84CvgYWx/7c2swOScA778D06QX4Oq62\nGbsi0hiqFxPffBM6kUjDJTTY0jk3FKhxEVjnXPcazo3C9zrU1P5z/EwQSYE774yQl/c50WhrYMfQ\ncURyTmUx0aMHHHEElJb6cRQimUprIueQzz6D5583otGB+HXCRCSEzp3hzTchGoXDDtOiVZLZVEjk\nkCFDHGYrgWVAx9BxRHLannv6YqJdOzjqKJgyJXQikcSokMgR5eUwbFiEaPQ+YquTi0hgbdvCq6/C\nwQf7Rauefz50IpH4qZDIEY88AuvWOfLypgJHho4jIjHbbQcvvwy9esHvfgf/+pdfwl4kUyRrZUtJ\nYxUVcNddG3DuaZw7D035FEkvzZpBSQnstRdcdx188AE89BA0bx46mUjd1CORA0aNgiVLmpKf/zSg\nDbpE0lFeHvzzn/D00/Dss3DMMbB4cehUInVTIZHlnIPbb9+A2QQikV8CWuNLJJ0VF8PUqbBokR87\n8d57oROJbJ0KiSw3eTLMnt0UsweAfqHjiEg9VBYQ7dv7Db/uvVfjJiR9qZDIYtEoXH11BXl504lG\nOwCFdT5HRNLDLrv4KaEXXQRXXAG//S18r12IJA2pkMhizzwDs2cXADcCfwodR0Ti1KwZDBkCL7zg\np4l27eqXuRdJJyokstT69XDttRXAi0Sjx6LlsEUy16mnwsyZft2JI46AW2+FjRtDpxLxVEhkqQce\ngC+/zCc/fzBweeg4ItJAu+8Or78O11wDN9wA3brB3Ll1Pk2k0amQyEI//AA33VQBPEIkci7QInQk\nEUmCggK47TZ46y1YtQq6dPELWEUioZNJLlMhkYVuv92xcmWEvLxngHNCxxGRJDv0UCgrgyuvhD/9\nyc/smDUrdCrJVSokssyiRTB4cJRo9E6i0WvQDu0i2alFC98b8eabsHIlFBX5wqK8PHQyyTUqJLLM\nDTdEiUTKMXsHOCF0HBFpZN26+d6IO+6ARx+FTp1gxAitOyGpo0Iii7z6Kjz5ZB6RyF9x7ia0p4ZI\nbigo8IMwP/zQL619zjlw5JEwbVroZJILVEhkifJyOPvsCvLypgBR4KDQkUQkxdq3h//8ByZOhNWr\n4bDD4LTTfIEh0lhUSGSJK6+MsmTJBsz+CtwVOo6IBNSjB8yY4W9xzJwJP/859O8PX34ZOplkIxUS\nWWD0aH9LIxq9nEjkdqBV6EgiElheHpx1lu+NuOsuvwtwhw6+oFiwIHQ6ySYqJDLcN99Anz4bMRsD\n7AQcHjqSiKSRZs3gj3+EhQvhllv8Lx4dO0Lv3jBvXuh0kg1USGQw5+CCCyKsWLECs4HAzaEjiUia\n2nZbuPZaX1AMHux3Bv6//4OTT4YJEzTLQxKnQiKDPf44jBmTTyTSn2j0bqBZ6EgikuZatoTLL4f5\n8/100UWL4Je/9EXFgw/6QZoi8VAhkaEmToT+/aPAI0AX4MDAiUQkkzRrBuef7wdjTpkC++4Ll14K\n7dr5rcvffVe9FFI/KiQyUFkZnHJKhEhkAjAeuC50JBHJUGZw1FF+MOaCBXDFFfDyy3DIIbD//nD3\n3X4slkhtVEhkmPnz4Ze/3Mj69bOAu4HhaBlsEUmG3XbzAzI/+wzGjvW9FNdeC7vsAscfD488At99\nFzqlpBsVEhlk6VLo0WMjP/zwBXAVzpUAzUPHEpEsk58PJ5wAI0fC4sXwwAMQjfqpo23bwkknwcMP\n+8dEVEhkiPJyOOGEjXz55fdAb6LREmD70LFEJMvtuCNceCFMmgRffeVvdaxe7cdRtGvnb4Hccovf\n7yMaDZ1WQlAhkQE+/hgOOmgj77+/Fud+RyTyCNAudCwRyTE77+wHZE6Z4sdN/PvfsPvucOed0KWL\nvwVy5pl+RtmiRaHTSqo0CR1Atu6VV+D00zeydu1CoC/R6GBgn9CxRCTH7bgjnH22PzZsgDfe8LPJ\nJk70+304B3vt5TcPO+II/+dee/nBnZJdVEikKef8tsB//rMDxmE2hGj0GWCX0NFERDbTtCl07+6P\nW2/1AzInT4apU+H11+GJJ/y/aW3awKGH+tshBx8MBx3kCxLJbCok0tDixXDppVFeeCEPv1rlVzj3\nElpwSkQywU9+Ar/7nT/Aj/GaNs33WkyfDgMHwg8/+Mf23BMOOGDzY/fd1XORSVRIpJHVq/3mOrfd\nFmHjxhXAhUB34G9AOvxfVQIUhw6RU0pKSigu1nueSnrPk6+w0M8COeEE/7Vz8OmnftGr996D0tIS\npk4t5ttv/eOtWvmpp1WPjh190dFMv0+lnYQGW5rZpWa20MzWmtnbZnZwHe1/b2bzYu1nm9mJNbS5\n2cy+NrM1ZjbBzPZKJFsmikTgscdgjz0q+PvfN7B+/UCcOwu/0NTFpEcRAb6QkFQqKdF7nmp6zxuf\nGey9tx+YOWgQdOhQwrJlfpvzl1+GG27wS3Z/9BHcdhuceip07gwtWvjeiuOO81NRb7sNSkp8b8fX\nX2vWSChx90iY2enAQPyvy9OBAUCpmXV0zi2voX034Gn8T8WXgTOB0WbWxTk3N9bmOuAy4FxgIfCP\n2DX3dc5tSOg7ywCffAJPPQWPPrqeL79sBowiP/9BIpEriUb/H+lTQIiINC4z2HVXf5x00qbzzvki\n4dNPNz+mT/frXFTeIgEoKPBTUtu333TssoufbVL5Z9u2sMMOfpt1SY5Ebm0MAIY554YDmNlFwK+A\nPsC/amh/JTDWOTco9vWNZvZLfOFwSZU2tzjnxsSu2RtYCpwKjEwgY1pyzq9MWVrqi4eZM5uRl7eK\naHQkeXn/IRr9NZHIeKBp6KgiImmhaoFx9NFbPl5eDp9/7lfj/OILv9bFl1/6P8vK/JizVas2f05+\nPrRuDTvttOnPn/zED/z8yU/8scMOsP32m47CQthuO/9c2VxchYSZFQBFwK2V55xzzswmAt1qeVo3\nfA9GVaVAr9g19wR2BiZVueYKM3sn9tyMLCTWrfPzqOfP9/cBp05dxzvvGCtXNsNsI86Nx2wk0agD\nziAafRkoCB1bRCSjFBb6PUH237/2NqtW+ZWBlyzxx7JlWx6ffupnm3z3HaxcWfu1ttnGFxTbbuv/\nbNVq82ObbfwOqy1bbvp7ixabH82bb340a7b50bRpZg02jbdHojV+Y4el1c4vBTrV8pyda2m/c+zv\nbQFXR5vqmgM8//w83nuv7tBV1bSbXdVzzm1+RKObjg0bYONG/2dFhS8WVq50lJdHKC+v4Icfoixf\nDsuW5bNqVdWlq1cC7wPvY/YFzuUBB+HcBcC2sTbvx/eNBFEOlKVBBpg/fz5lZaGzNL7y8vKc+D7T\nid7z1EvVe96iBeyxhz+2pqLCFx8rV/pj1SpYsQLWrPF/X71607F2rX9s6VL/+Nq1/mdD5bF2bWK7\nqDZpsukoKPBH1XNNmvjekZq+zsvb8u/5+Zv+Xl4+r/JlkrPHgnOu3gd+EYMocGi18/8C3qrlOeuB\n06uduwT4Ovb3bkAEaFutzUjg6VqueSa++NChQ4cOHTp0JHacGU8NUNsRb4/EcmI/9Kudb8OWPQqV\nltTRfgl+VGHbatdoA8ys5ZqlwFnAZ8C6euQWERERrzmwO/5naYPFVUg45yrMbAbQA3gRwMws9vU9\ntTxtWg2PHx87j3NuoZktibX5X+ya2wGHAvfXkuNb/EwQERERid9bybpQIrM2BgFPxgqKyumfLYEn\nAMxsOPClc+7PsfZDgClmdhV++mcxfsBmvyrXvBv4q5l9iu9luAX4EvhvAvlEREQkReIuJJxzI82s\nNX7t5rbALKCnc25ZrEl7YGOV9tPMrBj4Z+z4BOhVuYZErM2/zKwlMAy/N/brwInZvIaEiIhINjCX\nyHBSERERERJcIltEREQEVEiIiIhIA2RkIRHvpmGSODO70cyi1Y65dT9T6svMjjSzF83sq9j7e0oN\nbXJ2U7vGUNd7bmaP1/C5fyVU3mxgZteb2XQzW2FmS83sBTPrWK1NMzO738yWm9lKM3vOzNqEypzp\n6vmev1btcx4xs6HxvE7GFRJVNg27EegCzMZv8NU6aLDs9gF+YO3OseOIsHGyzjb4QcuX4heJ2UyV\nTe36A4cAq/GfeW3KkritvucxY9n8c6+9xRvmSOBe/NT+4/B7Aow3sxZV2tyN37vpt8BRQDtgVIpz\nZpP6vOcOeIhNn/VdgGvjeZGMG2xpZm8D7zjnrox9bcAi4B7nXE2bhkkDmNmN+Fk2XUNnyQVmFgVO\ndc69WOXc18CdzrnBsa+3wy/edq5zLiP3okkntbznjwOFzrnTwiXLbrFf/r4BjnLOvRH7XC8DznDO\nvRBr0wmYB/zCOTc9XNrsUP09j517FZjpnLsq0etmVI9ElU3Dqm7w5YCtbRomDbd3rAt4vpmNMLOf\nhg6UK8xsD2rY1A6o3NROGs8xse7gD81sqJn9JHSgLLM9/rfh72JfF+GXJKj6Wf8I+AJ91pOl+nte\n6SwzW2Zm75vZrdV6LOqUyIJUISWyaZg0zNvAecBH+C6vm4CpZvZz59zqgLlyxc74//Hj2dROGm4s\nvkt9IdABuA14xcy6uUzrxk1DsZ7ku4E3qqwptDOwIVYoV6XPehLU8p4DPAV8DnwN7I/fO6sj8Lv6\nXjvTConaGLXf55QGcM5VXYv9AzObjv/Q/QF4PEwqQZ/5RlXtltEcM3sfmA8cA7waJFR2GQp0pn7j\nrfRZT47K9/zwqiedc49U+XJObMuKiWa2h3NuYX0unFG3Nkhs0zBJIudcOfAxoFkDqVF1U7uq9JlP\nodg/qMvR577BzOw+4CTgGOfc11UeWgI0jY2VqEqf9Qaq9p4vrqP5O/h/c+r9Wc+oQsI5VwFUbhoG\nbLZpWNI2IJHamVkrfFdvXR9GSYLYD7DKTe2AzTa102c+RcysPbAj+tw3SOwHWi/gWOfcF9UenoHf\nXqHqZ70j8DNimzxK/Op4z2vSBd8DVO/Peibe2tjqpmGSXGZ2JzAGfztjV+Dv+P/ZS0LmyiZmtg2+\n+rfYqT3N7ADgO+fcIrSpXdJt7T2PHTfix0gsibW7A98Tl5Rtl3NRbG2CYuAUYLWZVfaylTvn1jnn\nVpjZo8AgM/seWInfNfpNzdhITF3vuZntCZwJvAJ8CxyA/xk7xTn3Qb1fyDmXcQdwCf4f1LX4SvWg\n0Jmy9cAXDF/G3usv8Nu37xE6VzYdwNFAFH/brurxWJU2N+EHQ63B/zDbK3TuTD629p4DzYFx+CJi\nHbAAeADYKXTuTD5qeb8jQO8qbZrh1z1Yji8kngXahM6eqUdd7zl+k83X8NNu1+AH1d8GtIrndTJu\nHQkRERFJHxk1RkJERETSiwoJERERSZgKCREREUmYCgkRERFJmAoJERERSZgKCREREUmYCgkRERFJ\nmAoJERERSZgKCREREUmYCgkRERFJmAoJERERSdj/BzZTMTsTpLO9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xaded8b6c>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "delay_dist = sp.gamma(10)\n",
    "delay_threshold = 7\n",
    "probability_delay_below_threshold = delay_dist.cdf(delay_threshold)\n",
    "\n",
    "delay_range = np.linspace(0,25,100)\n",
    "pds = delay_dist.pdf(delay_range)\n",
    "plt.plot(delay_range, pds)\n",
    "plt.fill_between(delay_range[delay_range < delay_threshold],0,pds[delay_range < delay_threshold])\n",
    "print(\"The probability that a given synapse has a delay less than {} is {:.2f}\".format(delay_threshold, probability_delay_below_threshold))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the probability that a given synapse will have a delay less than our threshold is 0.17.\n",
    "\n",
    "The number of synapses with a delay less than the threshold is a random variable distributed as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$N_\\texttt{delay} \\sim \\textrm{Binomial}(N_\\texttt{conns}, p=0.17)$,\n",
    "\n",
    "where $N_\\texttt{conns}$ is the random variable representing the number of synapses in the sub-row discussed above.\n",
    "\n",
    "This distribution is a *mixture distribution*. If $f_{N_\\texttt{conns}}(k)$ is the probability mass function giving the probability that $N_\\texttt{conns} = k$, then\n",
    "- with probability $f_{N_\\texttt{conns}}(0)$, $N_\\texttt{delay}$ is distributed as $\\textrm{Binomial}(0, p=0.17)$\n",
    "- with probability $f_{N_\\texttt{conns}}(1)$, $N_\\texttt{delay}$ is distributed as $\\textrm{Binomial}(1, p=0.17)$\n",
    "- with probability $f_{N_\\texttt{conns}}(2)$, $N_\\texttt{delay}$ is distributed as $\\textrm{Binomial}(2, p=0.17)$\n",
    "- ...\n",
    "- with probability $f_{N_\\texttt{conns}}(1024)$, $N_\\texttt{delay}$ is distributed as $\\textrm{Binomial}(1024, p=0.17)$\n",
    "\n",
    "(Note that, in our example, 1024 is the maximum possible number of synapses in a sub-row.)\n",
    "The distributions listed above are the *component distributions* of the mixture distribution.\n",
    "\n",
    "We need a quantile or ppf function for our mixture distribution. Unfortunately, there is no simple expression of the quantile function of mixture distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing the CDF of a mixture distribution\n",
    "\n",
    "However, there is a simple expression for the cdf of our mixture distribution.\n",
    "\n",
    "If $f_{N_\\texttt{conns}}(k)$ is the probability that the number of synapses in a sub-row is $k$, and\n",
    "\n",
    "$F_\\texttt{delay}(k; N_\\texttt{conns}=k^\\prime)$ is the cdf of the number of synapses whose delay is less than or equal to our threshold when the number of synapses in the sub-row is equal to $k^\\prime$, then the cdf of the mixture distribution can be written as\n",
    "\n",
    "$F_\\texttt{delay}(k) = \\sum_{k^\\prime = 1}^{1024} f_{N_\\texttt{conns}}(k^\\prime) F_{delay}(k; N_\\texttt{conns} = k^\\prime) \\quad .$\n",
    "\n",
    "That is, the cdf of the mixture distribution is a weighted sum of the cdfs of the component distributions, weighted by the probability of the parameter value that would give us that cdf.\n",
    "\n",
    "The following function computes cdf(k) for a mixture distribution which is composed of the component distributions $\\texttt{dists}$ weighted by the probabilities $\\texttt{ps}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def eval_mixture_cdf(dists, ps, k):\n",
    "    return sum(dist.cdf(k) * p for dist, p in zip(dists, ps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing the quantile function of a mixture distribution\n",
    "\n",
    "Note that the quantile function $Q(p)$ is the minimum value $k$ such that $F(k) \\ge p$.\n",
    "\n",
    "We can use the above expression for the cdf of a mixture distribution along with binary search over plausible values $k$ to compute the quantile function of a mixture distribution.\n",
    "\n",
    "First we will need to import $\\texttt{lazyarray}$, for creating arrays of lazily-evaluated elements, and the $\\texttt{bisect_left}$ which **TODO**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import lazyarray as la\n",
    "from bisect import bisect_left"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function computes the quantile function of the mixture distribution $\\texttt{second_dist}$ which has its $\\texttt{param_name}$ parameter distributed according to $\\texttt{first_dist}$, where both $\\texttt{first_dist}$ and $\\texttt{second_dist}$ are discrete.\n",
    "\n",
    "It works as follows.\n",
    "- First, we restrict ourselves to a range of \"plausible\" values for the parameter that is to be drawn from $\\texttt{first_dist}$, i.e., to values whose probability is greater than $\\epsilon=10^{-10}$\n",
    "- For each of those plausible values, we compute the probability of obtaining that value\n",
    "- We create an array of scipy distributions, one for each plausible value of the parameter $\\texttt{param_name}$\n",
    "- We need to return a value such that it is the lowest value that has cumulative probability for $\\texttt{second_dist}$ that is greater than $\\texttt{quantile}$. We will perform binary search over a plausible range of values. If it is the case that compontent distributions in our mixture distribution that have higher parameter values have lower values for a given quantile (as will be the case in our example), then the lowest value we need to consider is $\\texttt{ppf(quantile)}$ for the component distribution with the lowest parameter value, and the \n",
    "- We construct a lazily evaluated array $\\texttt{cdfs}$, such that $\\texttt{cdfs}$[$\\texttt{k}$] is the cumulative probability of the value $k$. The value of $\\texttt{cdfs}$[$\\texttt{k}$] for any $k$ will remain unevaluated until we access that element.\n",
    "- $\\texttt{bisect_left(list, value, min, max}$ returns the index after the last element in $\\texttt{list}$ (between $\\texttt{min}$ and $\\texttt{max}$) whose value is less than $\\texttt{value}$. We use it here to find the first value $k$ such that $\\texttt{cdf[k]}$ is greater than or equal to $\\texttt{quantile}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "epsilon = 1e-10\n",
    "\n",
    "def eval_mixture_ppf(quantile, first_dist, second_dist, param_name, **other_params):\n",
    "    plausible_first_values = np.arange(int(first_dist.ppf(epsilon)), int(first_dist.ppf(1-epsilon)) + 1)\n",
    "    plausible_first_value_ps = first_dist.pmf(plausible_first_values)\n",
    "    component_dists = [second_dist(**{**{param_name: k}, **other_params}) for k in plausible_first_values]\n",
    "    min_plausible_second_value = int(component_dists[0].ppf(quantile))\n",
    "    max_plausible_second_value = int(component_dists[-1].ppf(1-epsilon))\n",
    "\n",
    "    cdfs = la.larray(lambda k: eval_mixture_cdf(component_dists,\n",
    "                                                plausible_first_value_ps, k),\n",
    "                     shape=(max_plausible_second_value,))\n",
    "    k = bisect_left(cdfs, quantile, min_plausible_second_value, max_plausible_second_value)\n",
    "\n",
    "    return k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following use of $\\texttt{eval_mixture_ppf}$ computes a value such that with probability 0.9 a given sub-row will not have more than that many synapses with a delay less than our delay threshold.\n",
    "\n",
    "Also shown is the same value calculated by sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With probability 0.9, no sub-row will have more than 13 synapses with a delay less than our delay threshold. By sampling, we calculate the value to be 12.\n"
     ]
    }
   ],
   "source": [
    "exact_90th_percentile_delay_num = eval_mixture_ppf(0.9, N_conns_dist, sp.binom, 'n', p = probability_delay_below_threshold)\n",
    "\n",
    "n_conns_samples = N_conns_dist.rvs(20000)\n",
    "n_delays_samples = np.random.binomial(n=n_conns_samples, p=probability_delay_below_threshold)\n",
    "estimated_90th_percentile_delay_num = int(np.percentile(n_delays_samples, 90))\n",
    "\n",
    "print(\"With probability 0.9, no sub-row will have more than {} synapses with a delay less than our delay threshold. By sampling, we calculate the value to be {}.\".format(exact_90th_percentile_delay_num, estimated_90th_percentile_delay_num))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following use of $\\texttt{eval_mixture_ppf}$ compute a value such that with probability 0.9999 there will be no sub-row within our matrix which has more than that number of synapses with a delay less than our delay threshold. We cannot compute the same thing by sampling (which is the reason for this whole exercise), as the number of samples needed is too large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With probability 0.9999, no sub-row will have more than 30 synapses with a delay less than our delay threshold.\n"
     ]
    }
   ],
   "source": [
    "exact_extreme_delay_num = eval_mixture_ppf(max_full_matrix_quantile, N_conns_dist, sp.binom, 'n', p = probability_delay_below_threshold)\n",
    "\n",
    "print(\"With probability 0.9999, no sub-row will have more than {} synapses with a delay less than our delay threshold.\".format(exact_extreme_delay_num))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum of vector of random length"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "In the above, we noted that if $Q_X(p)$ is the $p$th quantile of random variable $X$, then the $p$th quantile of the maximum of $n$ random variables distributed identically to (and independently from) $X$ is\n",
    "\n",
    "$Q_n^\\texttt{max}(p) = Q_X(\\sqrt[n]{p}) \\quad .$\n",
    "\n",
    "But what if $n$, i.e., the number of samples of which we are taking the maximum, is itself a random variable?\n",
    "This is the case in the following example.\n",
    "\n",
    "Suppose we want to compute a quantile on the maximum delay in a given sub-row. In our example, the delays are independently drawn from the gamma distribution given above. How many delays are there? There are $N_\\texttt{conns}$ of them, which is itself a random variable.\n",
    "\n",
    "As noted above, if the cdf of a single delay is $F_X(x)$, then the cdf of the maximum of $n$ delays is $(F_X(x))^n$. When the number of synapses is a random variable, then the maximum delay is a mixture distribution, with cdf given by\n",
    "\n",
    "$F_{N_\\texttt{conns}}^\\texttt{max}(x) = \\sum_{k=0}^{1024} f_{N_\\texttt{conns}}(k) (F_X(x))^k \\quad .$\n",
    "\n",
    "We can use a modified version of the above trick to compute quantiles on this mixture distribution. First, the following function computes the cdf of the maximum of $n$ random variables following the $\\texttt{dist}$ distribution, where $n$ itself is a random variable. The array $\\texttt{vals}$ contains values $n$ can take with non-negligible probability, and the array $\\texttt{ps}$ contains the associated probability of each value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def eval_max_cdf(dist, vals, ps, k):\n",
    "    return sum(dist.cdf(k)**val * p for val, p in zip(vals, ps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the component distributions in this case are continuous, we need to perform binary search over a continuous variable. For this purpose, we import here our own function $\\texttt{continuous_bisect_fun_left(f,v,lo,hi)}$, which finds the maximum value $v^\\prime$, between $\\texttt{lo}$ and $\\texttt{hi}$, such that $f(v^\\prime)$ is less than $v$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from utilities import continuous_bisect_fun_left"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function computes the $\\texttt{quantile}$ quantile of the maximum of $n$ random variables distributed as $\\texttt{second_dist}$, where $n$ is distributed as $\\texttt{length_dist}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def eval_max_ppf(quantile, length_dist, second_dist):\n",
    "    plausible_lengths = np.arange(int(length_dist.ppf(epsilon)), int(length_dist.ppf(1-epsilon) + 1))\n",
    "    plausible_length_ps = length_dist.pmf(plausible_lengths)\n",
    "    lo = second_dist.ppf(epsilon)\n",
    "    hi = second_dist.ppf(1-epsilon)\n",
    "    \n",
    "    k = continuous_bisect_fun_left(lambda k: eval_max_cdf(second_dist,\n",
    "                                                          plausible_lengths,\n",
    "                                                          plausible_length_ps, k),\n",
    "                                   quantile, lo, hi)\n",
    "\n",
    "    return k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The use of $\\texttt{eval_max_ppf}$ below calculates the 90th percentile of the maximum delay in a single sub-row.\n",
    "\n",
    "For comparison, an estimate of the 90th percentile is calculated by sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With probability 0.9, the maximum delay of a given sub-row will be less than or equal to 21.49 ms. By sampling, we calculate the value to be 21.56 ms.\n"
     ]
    }
   ],
   "source": [
    "exact_90th_percentile_of_max_sub_row_delay = eval_max_ppf(0.9, N_conns_dist, delay_dist)\n",
    "\n",
    "n_conns_samples = N_conns_dist.rvs(10000)\n",
    "max_samples = np.array([delay_dist.rvs(k).max() for k in n_conns_samples])\n",
    "\n",
    "estimated_90th_percentile_of_max_sub_row_delay = np.percentile(max_samples, 90)\n",
    "\n",
    "print(\"With probability 0.9, the maximum delay of a given sub-row will be less than or equal to {:.2f} ms. By sampling, we calculate the value to be {:.2f} ms.\".format(exact_90th_percentile_of_max_sub_row_delay, estimated_90th_percentile_of_max_sub_row_delay))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The use of $\\texttt{eval_max_ppf}$ below calculates a delay value such that with probability 0.9999, no sub-row will contain a delay with greater than that value. We cannot calculate this by sampling, as the number of samples would be too large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With probability 0.9999, no sub-row in the matrix will contain a delay of greater than 44.63 ms\n"
     ]
    }
   ],
   "source": [
    "exact_extreme_of_max_sub_row_delay = eval_max_ppf(max_full_matrix_quantile, N_conns_dist, delay_dist)\n",
    "\n",
    "print(\"With probability 0.9999, no sub-row in the matrix will contain a delay of greater than {:.2f} ms\".format(exact_extreme_of_max_sub_row_delay))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}