willirath/profiling_Agulhasparticles.ipynb

## profiling_Agulhasparticles.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Profiling Parcels\n",
    "\n",
    "We'll run the Agulhas example with 12-hourly output and find that almost all the time is spent in the `ParticleFile.write` routine where the `ParticleSet` is sorted back to netCDF arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from parcels import (AdvectionRK4, ErrorCode, FieldSet,\n",
    "                     JITParticle, ParticleFile, ParticleSet)\n",
    "from datetime import timedelta\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wrap experiment in a function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def DeleteParticle(particle, fieldset, time, dt):\n",
    "    particle.delete()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_experiment(runtime=timedelta(days=5),\n",
    "                   output_timedelta=timedelta(hours=12),\n",
    "                   num_of_particles=1e3):\n",
    "    \n",
    "    filenames = \"GlobCurrent_example_data/20*.nc\"\n",
    "    variables = {'U': 'eastward_eulerian_current_velocity',\n",
    "                 'V': 'northward_eulerian_current_velocity'}\n",
    "    dimensions = {'lat': 'lat',\n",
    "                  'lon': 'lon',\n",
    "                  'time': 'time'}\n",
    "    fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)\n",
    "    \n",
    "    num_of_particles_sqrt = int(np.round(num_of_particles ** 0.5))\n",
    "    lons, lats = np.meshgrid(np.linspace(15, 35, num_of_particles_sqrt),\n",
    "                             np.linspace(-40, -30, num_of_particles_sqrt))\n",
    "    pset = ParticleSet(fieldset=fieldset, pclass=JITParticle, lon=lons, lat=lats)\n",
    "\n",
    "    output_file = ParticleFile(\n",
    "                 name=(\"agulhas_profiling\"),\n",
    "                 particleset=pset,\n",
    "                 outputdt=output_timedelta)\n",
    "    \n",
    "    pset.execute(AdvectionRK4,\n",
    "                 runtime=runtime,\n",
    "                 dt=timedelta(minutes=5),\n",
    "                 output_file=output_file,\n",
    "                 recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Do one reference run to gauge timing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Casting lon data to np.float32\n",
      "WARNING: Casting lat data to np.float32\n",
      "WARNING: Casting depth data to np.float32\n",
      "INFO: Compiled JITParticleAdvectionRK4 ==> /tmp/parcels-1000/fe8d600031ee877c32716e5c5bbf6ff8.so\n",
      "100% (432000.0 of 432000.0) |############| Elapsed Time: 0:00:49 Time:  0:00:49\n"
     ]
    }
   ],
   "source": [
    "run_experiment(num_of_particles=8000);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up profiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from line_profiler import LineProfiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "profile = LineProfiler()\n",
    "profile.add_function(ParticleFile.write);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Profile experiments wiht increasing number of particles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Compiled JITParticleAdvectionRK4 ==> /tmp/parcels-1000/b32a611f2b8d83b2a4bd2de9a70fa026.so\n",
      "INFO: Compiled JITParticleAdvectionRK4 ==> /tmp/parcels-1000/e4c018e2e204b2ea3fa6ecd34c04436b.so\n",
      "INFO: Compiled JITParticleAdvectionRK4 ==> /tmp/parcels-1000/eef5ef1bbff2b34629fdd1727d91e461.so\n",
      "100% (432000.0 of 432000.0) |############| Elapsed Time: 0:00:04 Time:  0:00:04\n",
      "INFO: Compiled JITParticleAdvectionRK4 ==> /tmp/parcels-1000/e1fe1157af36081d75e5edbf452757d2.so\n",
      "100% (432000.0 of 432000.0) |############| Elapsed Time: 0:00:16 Time:  0:00:16\n",
      "INFO: Compiled JITParticleAdvectionRK4 ==> /tmp/parcels-1000/b4922612a887126cb49d8b18846f8617.so\n",
      "100% (432000.0 of 432000.0) |############| Elapsed Time: 0:00:43 Time:  0:00:43\n",
      "INFO: Compiled JITParticleAdvectionRK4 ==> /tmp/parcels-1000/57ffe0235411c193dd83640067665ef7.so\n",
      "100% (432000.0 of 432000.0) |############| Elapsed Time: 0:01:37 Time:  0:01:37\n"
     ]
    }
   ],
   "source": [
    "lstats = {}\n",
    "for nop in [100, 500, 1000, 2000, 4000, 8000]:\n",
    "    lstats[nop] = profile.run(f\"run_experiment(num_of_particles={nop})\").get_stats()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wrap results in a list and plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "particle_numbers = np.array(sorted(lstats.keys()))\n",
    "lines = np.stack([[tpl[0] for tpl in list(lstats[pn].timings.items())[0][1]] for pn in particle_numbers])\n",
    "hits = np.stack([[tpl[1] for tpl in list(lstats[pn].timings.items())[0][1]] for pn in particle_numbers])\n",
    "totals = np.stack([[tpl[2] for tpl in list(lstats[pn].timings.items())[0][1]] for pn in particle_numbers])\n",
    "particle_numbers = particle_numbers[:, np.newaxis] + 0 * lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAD8CAYAAACihcXDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAF4xJREFUeJzt3X+QXWV9x/H3hwTCD4UEEp2YhCY4gRE6NsBKUQuNQOVHHUEHNRl/RIqNUKT+6FRBZ4x2ygw6IOhQwWgi0NEAAgWGWpEiwrSjYRIIEAhpwg9hTSSEFGirhmT32z/Oc/Gy7O49uffuOc/mfl4zZ/bc5577nO9mT7777Pc85xxFBGZmVq896g7AzMycjM3MsuBkbGaWASdjM7MMOBmbmWXAydjMLANOxmZmGXAyNjPLgJOxmVkGJtYdAMDUAyfE7Fl71h2GmY0Dqx/avjUipnXSx8nv2i+e3zZQdn93RMQpneyvjCyS8exZe3LfHQfXHYaZjQMTpm/4Vad9PL9toHTOmTB9w9RO91dGFsnYzKxKAQwyWHcYr+JkbGY9Jwh2RLkyRVWcjM2sJ3lkbGZWsyAYyOz2wU7GZtaTBnEyNjOrVQADTsZmZvXLbWRc6go8ScslbZG0tqltnqRfSlojaZWkY1K7JH1L0kZJD0k6aqyCNzNrRwA7IkotVSl7OfTVwNArUL4OfDUi5gFfTq8BTgXmpmUxcGXnYZqZdU8QDJRcqlIqGUfEvcC2oc3A/mn9AGBTWj8duDYKvwQmS5rejWDNzLoiYKDkUpVOasafAe6QdAlFUn9Hap8BPNO0XX9q29z8YUmLKUbOHDzDpWszq05xBV5eOrlr27nAZyNiFvBZYFlq1zDbvub3S0QsjYi+iOibdtCEDsIwM9tVYqDkUpVOkvEi4Oa0/iPgmLTeD8xq2m4mfyhhmJnVrjiBp1JLVTpJxpuAP0/rJwAb0vptwMfSrIpjgRcjYvNwHZiZ1aGYZ5zXyLhUsVbSCmA+MFVSP7AE+Gvgm5ImAr8n1X+BHwOnARuB3wJndTlmM7OODVY46i2jVDKOiIUjvHX0MNsGcF4nQZmZjaXGyDgnnsZgZj0nEAOZPXXOydjMetK4LFOYme1OAvFy5DWl1snYzHpOcdGHyxRmZrXzCTwzs5pFiIHwyNjMrHaDHhmbmdWrOIGXV/rLKxozswr4BJ6ZWSYGPM/YzKxevgLPzCwTg55NYWZWr+JGQU7GZma1CsQOXw5tZlavCHzRh5lZ/eSLPszM6hZ4ZGxmlgWfwDMzq1kg31zezKxuAezwvSnMzOqm8Xc/Y0nLgfcAWyLij1Pb9cBhaZPJwAsRMU/SbGAdsD6998uIOKfbQZuZdSIYn1fgXQ1cAVzbaIiIDzXWJV0KvNi0/eMRMa9bAZqZjYVujYwlfRb4BEWOfxg4C5gOXAccCNwPfDQiXh6tn5a/GiLiXmDbCEEI+CCwYleCNzOrU4QYjD1KLaORNAP4W6AvVQ4mAAuArwGXRcRc4L+Bs1vF1Ok4/Tjg2YjY0NQ2R9IDku6RdFyH/ZuZdV1xAm9CqaWEicA+kiYC+wKbgROAG9P71wBnlOmkEwt59ah4M3BwRDwv6WjgFklHRMRLQz8oaTGwGODgGT6PaGZV6s4z8CLi15IuAZ4Gfgf8FFhNcR5tZ9qsH5jRqq+2o0m/Bd4PXN8U2PaIeD6trwYeBw4d4ZtYGhF9EdE37aC8bthhZru34gSeSi3AVEmrmpbFjX4kTQFOB+YAbwL2A04dYZej6mRIehLwWET0NwU2DdgWEQOSDgHmAk90sA8zszGxC1fgbY2IvhHeOwl4MiKeA5B0M/AOYLKkiWl0PBPY1GonLaORtAL4BXCYpH5JjUL0Al574u544CFJD1LUS86JiGFP/pmZ1aVxBV7JkfFongaOlbRvmtBwIvAocDdwZtpmEXBrq45ajowjYuEI7R8fpu0m4KZWfZqZ1a0bDySNiJWSbqSYvrYTeABYCvwrcJ2kf0xty1r15TNnZtZzImDHYHcu+oiIJcCSIc1PAMfsSj9OxmbWc4oyxfi7As/MbLcz7u5NYWa2u2lMbcuJk7GZ9SCXKczMsuBn4JmZ1ayYTZHXlb9OxmbWc/zYJTOzTLhMYWZWM8+mMDPLhGdTmJnVLELsdDI2M6ufyxRmZjVzzdjMLBNOxmZmNfM8YzOzTHiesZlZzSJgZ5duLt8tTsZm1pNcpjAzq5lrxmZmmYjMknHLoomk5ZK2SFrb1Ha9pDVpeUrSmqb3LpS0UdJ6SSePVeBmZp0YRKWWqpQZGV8NXAFc22iIiA811iVdCryY1g8HFgBHAG8C/l3SoREx0MWYzcw6EjEOa8YRca+k2cO9J0nAB4ETUtPpwHURsR14UtJGisdV/6Ir0ZqZdYUYyGw2RafRHAc8GxEb0usZwDNN7/enNjOzrESo1FKVTk/gLQRWNL0eLvIY7oOSFgOLAQ6e4fOIZladHO9N0fbIWNJE4P3A9U3N/cCsptczgU3DfT4ilkZEX0T0TTsor2dRmdluLoq6cZmlKp2UKU4CHouI/qa224AFkiZJmgPMBe7rJEAzs7Ew7mZTSFoBzAemSuoHlkTEMopZE80lCiLiEUk3AI8CO4HzPJPCzHITGZ7AKzObYuEI7R8fof0i4KLOwjIzG1tVliDK8JkzM+tJuV2B52RsZj2nODnnZGxmVrvcprY5GZtZT3LN2MysZoEYHG+zKczMdkeZDYw7vjeFmdn4E929N4WkyZJulPSYpHWS3i7pQEl3StqQvk4ZrY8sRsabdu7DkueOqDsMMxsXNrTepIzuDo2/CfwkIs6UtBewL/BF4K6IuFjSBcAFwBdG6sAjYzPrSd0aGUvaHzgeWFb0Gy9HxAsUtxS+Jm12DXDGaP1kMTI2M6tSAIODpae2TZW0qun10ohY2vT6EOA54PuS/gRYDXwaeGNEbAaIiM2S3jDaTpyMzaz3BFB+nvHWiOgb5f2JwFHA+RGxUtI3KUoSu8RlCjPrSV28hWY/0B8RK9PrGymS87OSpgOkr1tG68TJ2Mx6U5RcWnUT8RvgGUmHpaYTKe5ceRuwKLUtAm4drZ8syhTbfr8v1z12dN1hmNm4cEsX+uj6I5XOB36QZlI8AZxFMdi9QdLZwNPAB0brIItkbGZWuS5ObYuINcBwdeUTy/bhZGxmvScgys+mqISTsZn1KCdjM7P6ZXZzCidjM+tNTsavpd/uwZ4P7ld3GGbWK3btoo9KZJGMzcyqltvN5Vte9CFpuaQtktYOaT9f0npJj0j6emqbLel3ktak5aqxCtzMrCODKrdUpMzI+GrgCuDaRoOkd1HckeitEbF9yA0wHo+IeV2N0sysyzTeRsYRcS+wbUjzucDFEbE9bTPqNddmZlkpeyl0hQm73XtTHAocJ2mlpHskva3pvTmSHkjtx3UhRjOzLlNxAq/MUpF2T+BNBKYAxwJvo7j++hBgM3BwRDwv6WjgFklHRMRLQzuQtBhYDDBx/1GfRmJm1n2ZlSnaTcb9wM0REcB9kgaBqRHxHNAoXayW9DjFKHrV0A7SzZmXArz+gJkx7f4duxRAbvUeM6vGum51NNitjrqj3TLFLcAJAJIOBfYCtkqaJmlCaj8EmEtxByMzs3w05hmPpzKFpBXAfIpHj/QDS4DlwPI03e1lYFFEhKTjgX+QtBMYAM6JiKEn/8zMapfbX9ctk3FELBzhrY8Ms+1NwE2dBmVmNuYyS8Z+0oeZWQZ8ObSZ9aRxV6aowh7/+zv2vmdt6w27bTCzn4aZVSOo9FLnMrJIxmZmlctsLOZkbGY9yWUKM7McOBmbmWXAydjMrF4KlymGF3hmg5lVy7MpzMzq55GxmVkOnIzNzGrmmrGZWSacjM3M6qfMbi6fRTKOCGLnrj3pIzz7wsx2I1kkYzOzymU2nnMyNrPe4xN4ZmaZcDI2M8tAZsnYj10ys54jitkUZZZS/UkTJD0g6fb0eo6klZI2SLpe0l6t+shmZOzZEfaKyGzOke1+ul8z/jSwDtg/vf4acFlEXCfpKuBs4MrROvDI2Mx6U5RcWpA0E/hL4HvptYATgBvTJtcAZ7Tqp2UylrRc0hZJa4e0ny9pvaRHJH29qf1CSRvTeye3/lbMzGrQpWQMXA58Hmj8SXcQ8EJE7Eyv+4EZrTopMzK+GjiluUHSu4DTgbdGxBHAJan9cGABcET6zLclTSixDzOzSjXuadxqAaZKWtW0LH6lD+k9wJaIWN3c9TC7a5nWW9aMI+JeSbOHNJ8LXBwR29M2W1L76cB1qf1JSRuBY4BftNqPmVmlyteMt0ZE3wjvvRN4r6TTgL0pasaXA5MlTUyj45nAplY7abdmfChwXDpbeI+kt6X2GcAzTduNODyXtLjxm2YH29sMw8ysDdGd2RQRcWFEzIyI2RRVgZ9FxIeBu4Ez02aLgFtbhdRuMp4ITAGOBf4euCEVrUsPzyNiaUT0RUTfnkxqMwwzszZ1r2Y8nC8An0vVgYOAZa0+0O7Utn7g5ogI4D5Jg8DU1D6rabtSw3Mzs6p1+3LoiPg58PO0/gRFiba0dkfGt1BM3UDSocBewFbgNmCBpEmS5gBzgfva3IeZ2dgZ25HxLms5Mpa0AphPcUaxH1gCLAeWp+luLwOL0ij5EUk3AI8CO4HzImJgrII3M2tLxYm2jDKzKRaO8NZHRtj+IuCiToIyMxtLwndtMzPLgpOxmVkOnIzNzDLgZGxmVjM/6cPMLBNOxmZm9St74/iqOBmbWU9ymcLMrG7j8aIPM7PdkpOxmVm9fAWemVkmlNlDkJ2Mzaz3uGZsZpYHlylyEplNNDSz6jgZm5nVzyNjM7McOBmbmdUsfDm0mVntPM/YzCwXkVc2djI2s57kkbGZWd0yvOhjj1YbSFouaYuktU1tX5H0a0lr0nJaap8t6XdN7VeNZfBmZu3SYLmlKmVGxlcDVwDXDmm/LCIuGWb7xyNiXqeBmZmNpXE3myIi7pU0e+xDMTOrSJDdCbyWZYpRfErSQ6mMMaWpfY6kByTdI+m4kT4sabGkVZJW7WB7B2GYme06RbmlKu0m4yuBNwPzgM3Apal9M3BwRBwJfA74oaT9h+sgIpZGRF9E9O3JpDbDMDNrU5RcKtJWMo6IZyNiICIGge8Cx6T27RHxfFpfDTwOHNqtYM3MuqFx0ce4HxlLmt708n3A2tQ+TdKEtH4IMBd4otMgzcy6KgINlluq0vIEnqQVwHxgqqR+YAkwX9I8ikH8U8An0+bHA/8gaScwAJwTEdvGIG4zs87kdf6u1GyKhcM0Lxth25uAmzoNysxsrPkKPDOzugXgZ+CNwE/dMLMq5ZWLO5pnbGY2bnVrNoWkWZLulrRO0iOSPp3aD5R0p6QN6euU0fpxMjazntTF2RQ7gb+LiLcAxwLnSTocuAC4KyLmAnel1yNyMjaz3lP2go8SuTgiNkfE/Wn9f4B1wAzgdOCatNk1wBmj9ZNPzdjMrCLFRR+li8ZTJa1qer00IpYO229xH58jgZXAGyNiMxQJW9IbRtuJk7GZ9abycwa2RkRfq40kvY5iau9nIuIlSbsUjssUZtaTFFFqKdWXtCdFIv5BRNycmp9tXK2cvm4ZrQ8nYzPrPV2sGasYAi8D1kXEN5reug1YlNYXAbeO1o/LFGbWg7p634l3Ah8FHpa0JrV9EbgYuEHS2cDTwAdG68TJ2Mx6U5duLh8R/0FxTnA4J5btx8nYzHpPjMPHLpmZ7ZYye+ySk7GZ9aa8crGTsZn1Jg3mVadwMjaz3hPsykUflXAyNrOeI8pf0FEVJ2Mz601OxmZmGXAyNjOrmWvGZmZ5yG02RcsbBUlaLmmLpLVNbV+R9GtJa9JyWtN7F0raKGm9pJPHKnAzs/ZFUaYos1SkzF3brgZOGab9soiYl5YfA6RHjSwAjkif+bakCd0K1sysK4Lxl4wj4l5gW8n+Tgeui4jtEfEksBE4poP4zMzGxmDJpSKd3M/4U5IeSmWMxlNPZwDPNG3Tn9rMzLLSzZvLd0O7yfhK4M3APGAzcGlqH+42csN+N5IWS1oladUOtrcZhplZmzIrU7Q1myIinm2sS/oucHt62Q/Matp0JrBphD6WAksB9teBeU34M7PdWwQMjLPZFMNpPNcpeR/QmGlxG7BA0iRJc4C5wH2dhWhmNgbG28hY0gpgPsXjqvuBJcB8SfMoShBPAZ8EiIhHJN0APArsBM6LiIGxCd3MrAPj7Qq8iFg4TPOyUba/CLiok6DMzMZUAN17Bl5X+Ao8M+tBAZFXzdjJ2Mx6T5DdCTwnYzPrTeOtZmxmtltyMjYzq1u109bKcDI2s94TQGa30HQyNrPe5JGxmVnd8rsc2snYzHpPQHiesZlZBnwFnplZBlwzNjOrWYRnU5iZZcEjYzOzugUxkNfdfZ2Mzaz3+BaaZmaZ8NQ2M7N6BRAeGZuZ1Sx8c3kzsyzkdgJPkcH0DknPAb+qafdTga017buZ43i1HOLIIQZwHEMdFhGv76QDST+h+H7K2BoRp3SyvzKySMZ1krQqIvoch+PIMQbHkW8c3bZH3QGYmZmTsZlZFpyMYWndASSO49VyiCOHGMBxDJVLHF3V8zVjM7MceGRsZpaB3S4ZS1ouaYuktU1tB0q6U9KG9HVKapekb0naKOkhSUc1fWZR2n6DpEVtxDFL0t2S1kl6RNKn64hF0t6S7pP0YIrjq6l9jqSVqc/rJe2V2iel1xvT+7Ob+rowta+XdHIb/yYTJD0g6fa6Ykh9PCXpYUlrJK1KbXUcI5Ml3SjpsXScvL2G4+Ow9O/QWF6S9Jka4vhsOj7XSlqRjttajo/aRMRutQDHA0cBa5vavg5ckNYvAL6W1k8D/g0QcCywMrUfCDyRvk5J61N2MY7pwFFp/fXAfwGHVx1L6u91aX1PYGXq/wZgQWq/Cjg3rf8NcFVaXwBcn9YPBx4EJgFzgMeBCbv4b/I54IfA7el15TGkfp4Cpg5pq+MYuQb4RFrfC5hcRxxN8UwAfgP8UZVxADOAJ4F9mo6Lj9d1fNS11B7AmHxTMJtXJ+P1wPS0Ph1Yn9a/Aywcuh2wEPhOU/urtmszpluBv6gzFmBf4H7gTykm709M7W8H7kjrdwBvT+sT03YCLgQubOrrle1K7nsmcBdwAnB76rPSGJo+9xSvTcaV/lyA/SkSkOqMY8i+3w38Z9VxUCTjZygS+cR0fJxc1/FR17LblSlG8MaI2AyQvr4htTcOgob+1DZSe1vSn1FHUoxKK48llQfWAFuAOylGDC9ExM5h+nxlf+n9F4GDuhDH5cDngcYNAQ6qIYaGAH4qabWkxamt6p/LIcBzwPdT6eZ7kvarIY5mC4AVab2yOCLi18AlwNPAZoqf92rqOz5q0SvJeCQapi1Gad/1HUivA24CPhMRL9URS0QMRMQ8itHpMcBbRumz63FIeg+wJSJWNzdXGcMQ74yIo4BTgfMkHT/KtmMVy0SKctqVEXEk8H8U5YCq4yg6L+qx7wV+1GrTbseR6tGnU5QW3gTsR/GzGam/Mf9/W4deScbPSpoOkL5uSe39wKym7WYCm0Zp3yWS9qRIxD+IiJvrjAUgIl4Afk5R65ssqXGjqOY+X9lfev8AYFuHcbwTeK+kp4DrKEoVl1ccwysiYlP6ugX4F4pfUFX/XPqB/ohYmV7fSJGc6zo+TgXuj4hn0+sq4zgJeDIinouIHcDNwDuo6fioS68k49uAxtndRRT120b7x9IZ4mOBF9OfZHcA75Y0Jf3WfndqK02SgGXAuoj4Rl2xSJomaXJa34fiwF8H3A2cOUIcjfjOBH4WRQHuNmBBOpM9B5gL3Fcmhoi4MCJmRsRsij+FfxYRH64yhgZJ+0l6fWOd4t9zLRX/XCLiN8Azkg5LTScCj1YdR5OF/KFE0dhfVXE8DRwrad/0/6bxb1H58VGruovW3V4oDqjNwA6K35RnU9ST7gI2pK8Hpm0F/BNFDfVhoK+pn78CNqblrDbi+DOKP5EeAtak5bSqYwHeCjyQ4lgLfDm1H0JxoG6k+NN0UmrfO73emN4/pKmvL6X41gOntvnzmc8fZlNUHkPa54NpeQT4Umqv4xiZB6xKP5tbKGYh1BHHvsDzwAFNbVUfp18FHkvH6D9TzIio5Rita/EVeGZmGeiVMoWZWdacjM3MMuBkbGaWASdjM7MMOBmbmWXAydjMLANOxmZmGXAyNjPLwP8DY5dvLITeVCsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1061914898>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.pcolormesh(\n",
    "    particle_numbers,\n",
    "    lines,\n",
    "    (totals / totals.sum(axis=1)[:, np.newaxis] * 100).cumsum(axis=1))\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAHwCAYAAABZrD3mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3X20XFWZ7/vvj0ACXhISDGBCEnYAwz3kDnmLCHQHEG0RtRtfoC8eBy8H6HTbHQ9i03R4OWIPYdxoE1AIKukmoHSL2MqBNIjSnmNA9BJIIiEEpAVMTEIIhpdABBHIc/9Ysy4rRe29q/aunaq59u8zRo29as65quaTubOfWnPNWksRgZmZmeVlh053wMzMzFrnBG5mZpYhJ3AzM7MMOYGbmZllyAnczMwsQ07gZmZmGXICN8ucpC2S9u2nTY+kkLTj9uqXmQ0tJ3CzjEhaLOnscllE7BoRT3aqT2WSFkh6TNJWSWc0qD9X0tOSNktaKGlUqa5H0k8kvSzpl5Lev107b5YZJ3CzDKiQw//XFcBfA8vrKyQdD8wB3gf0APsC/1BqchPwC+DtwEXA9yTtMcT9NctWDn8QzLIlabWkCyQ9Iul5SddL2jnVjZN0u6TfprrbJU0q7btY0mWSfga8DNwIzATmp2nz+aldSNo/be8iaZ6kNeko915JuzTo126SrpO0QdJ6SZdKGpHq9pd0d9p/k6Sbm403Iq6JiP8F/L5B9enAdRGxKiKeB74InJHecxpwKHBJRLwSEd8HVgKfaPa9zYYbJ3Czofcp4HhgP2AacHEq3wG4HtgHmAK8Asyv2/dUYBYwmiLZ/RSYnabNZzd4r8uBw4CjgN2B84GtDdp9E3gd2B84BPgAUJua/yJwFzAOmARcXdspfciY01zYbzGd4gi9ZgWwl6S3p7onI+KluvrpA3wvs8rzghazoTc/ItYCSLqMIiFeHBHPAt+vNUp1P6nb94aIWFVq0+ubpCn2M4EjImJ9Kv55/X6S9gJOAMZGxCvA7yRdSfFB4VrgNYoPFRMjYh1wb23fiPhIS5Fva1dgc+l5bXt0g7pa/d6DeD+zSvMRuNnQW1vaXgNMBJD0NknXpunuF4F7gLG1qewG+/ZnPLAz8EQ/7fYBdgI2SHpB0gsUiXvPVH8+IOB+SaskndlCH/qyBRhTel7bfqlBXa3+JcysISdws6E3ubQ9BXgqbf8tcADwnogYAxydysuH2fW3C+zr9oGbKM4979dPf9YCrwLjI2JseoyJiOkAEfF0RPxFREwE/hL4Wu0c+yCtAg4qPT8I2JhmIlYB+0oaXVe/CjNryAncbOj9jaRJknYHLgRqi8JGU5z3fiHVXdLEa22kWL39FhGxFVgIXCFpoqQRko4sf1UrtdtAcY57nqQxknaQtJ+kYwAknVxaTPc8xYeGN5oJVNLItEhPwE6Sdi6tnv8WcJakAyWNo1gLcEPq038CDwKXpH0+BryL0ikGM9uWE7jZ0Ps2RcJ8Mj0uTeVfAXahOHK+D/hhE6/1VeCktGr9qgb151Gs3n4AeA74Eo3/n58GjAQeoUjS3wMmpLp3A0skbQEWAedExK8BJN0p6cI++ncXxYeSo4AFaftogIj4IfBlivP8a9Kj/KHlFGBG6s9c4KSI+G0f72U2rCmirxk5MxsMSauBsyPix53ui5lVi4/AzczMMuQEbmZmlqF+E7ikyen6xI+mr5Sck8p3l/Qfkn6Vfo4b+u6a5SUiejx9bmZDod9z4JImABMiYnn6iscy4KMUV4V6LiLmpiszjYuIvx/qDpuZmdkAFrFJuo3ico/zgWMjYkNK8osj4oC+9h0/fnz09PQMtK9mZmaVt2zZsk0R0e+NfFq6lKqkHorrJi8B9krfJyUl8T172WcWxSUamTJlCkuXLm3lLc3MzIYVSWuaadd0Ape0K8VFFT4bES/2dU3msohYQPF9UGbMmOHvrJl1UM+cO/pts3ruh4f8tYeyH2bDRVMJXNJOFMn7XyPillS8UdKE0hT6M0PVSTPb/rolyXZLP8y6TTOr0AVcBzwaEVeUqhZR3N+X9PO29nfPzMzMGmnmCPyPKO5JvFLSg6nsQopLHX5X0lnAb4CTh6aLZmZmVq/fBB4R97Lt3ZHK3tfe7piZmVkzfCU2MzOzDDmBm5mZZcgJ3MzMLENO4GZmZhlyAjczM8uQE7iZmVmGnMDNzMwy5ARuZmaWISdwMzOzDDmBm5mZZcgJ3MzMLENO4GZmZhlyAjczM8uQE7iZmVmGnMDNzMwy5ARuZmaWISdwMzOzDDmBm5mZZcgJ3MzMLENO4GZmZhlyAjczM8uQE7iZmVmGnMDNzMwy5ARuZmaWISdwMzOzDDmBm5mZZcgJ3MzMLENO4GZmZhlyAjczM8vQjoPZWdIHga8CI4B/joi5bemV2TDWM+eOftusnvvh7dATM+tmAz4ClzQCuAY4ATgQ+KSkA9vVMTMzM+vdYI7ADwcej4gnASR9BzgReKQdHRsKPrIZuGb+7cD/fr0Zyt89/15vy/8eNlwoIga2o3QS8MGIODs9PxV4T0TMrms3C5iVnh4APDaAtxsPbBpQR/NR9RirHh84xqqoeoxVjw/yj3GfiNijv0aDOQJXg7K3fBqIiAXAgkG8D5KWRsSMwbxGt6t6jFWPDxxjVVQ9xqrHB8MjRhjcKvR1wOTS80nAU4PrjpmZmTVjMAn8AeCdkqZKGgmcAixqT7fMzMysLwOeQo+I1yXNBn5E8TWyhRGxqm0929agpuAzUfUYqx4fOMaqqHqMVY8PhkeMA1/EZmZmZp3jK7GZmZllyAnczMwsQ12RwCUtlPSMpIcb1J0nKSSNT8+PlbRZ0oPp8fnt3+PWNIpP0hckrS/F8aFS3QWSHpf0mKTjO9Pr1rQSo6QeSa+Uyr/RuZ43r7ffU0mfSWO1StKXS+VZjWMr8VVpDCXdXIpjtaQHS3VZjSG0FmPFxvFgSfelOJZKOjyVS9JVaRwfknRo53reZhHR8QdwNHAo8HBd+WSKRXJrgPGp7Fjg9k73ebDxAV8AzmvQ9kBgBTAKmAo8AYzodAxtjrGnfqxzePQS43uBHwOj0vM9cx3HFuOrzBjW1c8DPp/rGA4gxsqMI3AXcELa/hCwuLR9J8W1S44AlnS6/+16dMUReETcAzzXoOpK4HwaXCAmJ33E18iJwHci4tWI+DXwOMVla7taizFmqZcYPw3MjYhXU5tnUnl249hifFnq6/dUkoA/B25KRdmNIbQcY5Z6iTGAMWl7N968LsmJwLeicB8wVtKE7dPTodUVCbwRSX8GrI+IFQ2qj5S0QtKdkqZv77610ew0pbNQ0rhUtjewttRmXSrLVaMYAaZK+oWkuyXN7FjvBm8aMFPSkhTLu1N5Vcaxt/igOmNYMxPYGBG/Ss+rMoZl9TFCdcbxs8A/SloLXA5ckMqrOI5AlyZwSW8DLgIand9eTnGd2IOAq4Fbt2ff2ujrwH7AwcAGimktaPIStZnoLcYNwJSIOAT4HPBtSWMav0TX2xEYRzE193fAd9NRTlXGsbf4qjSGNZ9k2yPTqoxhWX2MVRrHTwPnRsRk4FzgulRexXEEujSBU/zRnwqskLSa4jKtyyW9IyJejIgtABHxA2AnpQVuOYmIjRHxRkRsBf6JN6fmKnOJ2t5iTFOSz6btZRTnFqd1rqeDsg64JU3P3Q9spbiRQlXGsWF8FRtDJO0IfBy4uVRclTEEGsdYsXE8Hbglbf8bFfybWq8rE3hErIyIPSOiJyJ6KAbg0Ih4WtI70hEAaZXhDsCzHezugNSdg/kYUFtNuQg4RdIoSVOBdwL3b+/+tUNvMUraQ8X95JG0L0WMT27/HrbFrcBxAJKmASMp7oJUlXFsGF/FxhDg/cAvI2JdqawqY1jzlhgrNo5PAcek7eOA2mmCRcBpaTX6EcDmiNjQiQ62XadX0UWxSvAmiqmc1yiS9Vl19at5cxX6bGAVxerQ+4CjOt3/gcQH3AisBB6i+AWbUGp/EcUn4cdIqyq7/dFKjMAnSmO4HPjTTvd/EDGOBP6F4sPJcuC4XMexlfiqNIap/Abgrxq0z2oMW42xSuMI/DGwLMWyBDgstRVwTRrHlcCMTve/XQ9fStXMzCxDXTmFbmZmZn1zAjczM8uQE7iZmVmGnMDNzMwy5ARuZmaWISdwMzOzDDmBm5mZZcgJ3MzMLENO4GZmZhlyAjczM8uQE7iZmVmGnMDNMidpS7qTVF9teiRFuqWkmVWAE7hZRiQtlnR2uSwido2IrrgFZPqQ8Lv0oWKLpH8u1UnSlyQ9mx5frt0aONUfLGmZpJfTz4M7E4VZHpzAzTKQkl8u/18PSh8qdo2I8oeNWcBHgYOAdwEfAf4SQNJI4DaKW5eOA74J3JbKzayBXP4gmGVJ0mpJF0h6RNLzkq6XtHOqGyfpdkm/TXW3S5pU2nexpMsk/Qx4meL+6jOB+enodn5qF5L2T9u7SJonaY2kzZLulbRLg37tJuk6SRskrZd0qaQRqW5/SXen/TdJurlN/xynA/MiYl1ErAfmAWekumOBHYGvRMSrEXEVxX2cj2vTe5tVjhO42dD7FHA8sB8wDbg4le8AXA/sA0wBXgHm1+17KsWR62iKZPdTYHY6up3d4L0uBw4DjgJ2B84HtjZo903gdWB/4BDgA0DtaPmLwF0UR8KTgKtrO6UPGXP6ifceSU9LukVST6l8OrCi9HxFKqvVPRQRUap/qFRvZnWcwM2G3vyIWBsRzwGXAZ8EiIhnI+L7EfFyRLyU6o6p2/eGiFgVEa9HxGt9vUmaYj8TOCci1kfEGxHx84h4ta7dXsAJwGcj4ncR8QxwJXBKavIaxYeKiRHx+4i4t7ZvRHwkIub20Y1jgB7g/wSeAm4vLZzbFdhcarsZ2DWdB6+vq9WP7itms+HMCdxs6K0tba8BJgJIepuka9N094vAPcDY2lR2g337Mx7YGXiin3b7ADsBGyS9IOkF4Fpgz1R/PsX09f2SVkk6s9kORMQ9EfGHiHgBOAeYCvyXVL0FGFNqPgbYko666+tq9S81+95mw40TuNnQm1zankJxZArwt8ABwHsiYgxwdCpXqX15SrnR87JNwO8ppur7shZ4FRgfEWPTY0xETAeIiKcj4i8iYiLFIrOv1c6xD0DwZjyrKBaw1RyUymp17yqvSqdY6LYKM2vICdxs6P2NpEmSdgcuBGqLwkZTnPd+IdVd0sRrbQQafuc7IrYCC4ErJE2UNELSkZJG1bXbQHGOe56kMZJ2kLSfpGMAJJ1cWkz3PEUSfqO/jkmanr4KNkLSrhSL1NYDj6Ym3wI+J2lvSRMpPsDckOoWp/f475JGSaqd3//f/f6LmA1TTuBmQ+/bFAnzyfS4NJV/BdiF4sj5PuCHTbzWV4GT0qr1qxrUnwesBB4AngO+ROP/56cBI4FHKJL094AJqe7dwBJJW4BFFOfUfw0g6U5JF/bSt70oPpy8mOLsAT5SOnd/LfDvqX8PA3ekMiLiDxRfMTsNeIHiXP5HU7mZNaBtF32aWTtJWg2cHRE/7nRfzKxafARuZmaWISdwMzOzDPWbwCVNlvQTSY+mr5Sck8p3l/Qfkn6Vfo4b+u6a5SUiejx9bmZDod9z4JImABMiYrmk0cAyisUmZwDPRcTcdGWmcRHx90PdYTMzM2viCDwiNkTE8rT9EsVXQvYGTqS4HCPp50eHqpNmZma2rZZWoafrGt8D/F/AbyJibKnu+Yjocxp9/Pjx0dPTM6COmpmZDQfLli3bFBF79Ndux/4a1KQLM3yf4vrJL257waQ+95tFcTMGpkyZwtKlS5t9SzPLSM+cO/pts3ruh1tuazbcSFrTTLumEriknSiS979GxC2peKOkCRGxIZ0nf6bRvhGxAFgAMGPGDH/p3KyDckyy3dIPs27TzCp0AdcBj0bEFaWqRRT39yX9vK393TMzM7NGmjkC/yOKexKvlPRgKrsQmAt8V9JZwG+Ak4emi2ZmZlav3wSe7gXc2wnv97W3O2ZmZtYMX4nNzMwsQ07gZmZmGXICNzMzy5ATuJmZWYacwM3MzDLkBG5mZpYhJ3AzM7MMOYGbmZllyAnczMwsQ07gZmZmGXICNzMzy5ATuJmZWYacwM3MzDLkBG5mZpYhJ3AzM7MMOYGbmZllyAnczMwsQ07gZmZmGXICNzMzy5ATuJmZWYacwM3MzDLkBG5mZpYhJ3AzM7MMOYGbmZllyAnczMwsQ07gZmZmGXICNzMzy5ATuJmZWYacwM3MzDK042B2lvRB4KvACOCfI2JuW3plNoz1zLmj3zar5354O/TEzLrZgI/AJY0ArgFOAA4EPinpwHZ1zMzMzHo3mCPww4HHI+JJAEnfAU4EHmlHx6y7NHNUCD4yNDPbXhQRA9tROgn4YEScnZ6fCrwnImbXtZsFzEpPDwAeG8DbjQc2Daij+ah6jFWPDxxjVVQ9xqrHB/nHuE9E7NFfo8EcgatB2Vs+DUTEAmDBIN4HSUsjYsZgXqPbVT3GqscHjrEqqh5j1eOD4REjDG4V+jpgcun5JOCpwXXHzMzMmjGYBP4A8E5JUyWNBE4BFrWnW2ZmZtaXAU+hR8TrkmYDP6L4GtnCiFjVtp5ta1BT8JmoeoxVjw8cY1VUPcaqxwfDI8aBL2IzMzOzzvGV2MzMzDLkBG5mZpahrkjgkhZKekbSww3qzpMUksan58dK2izpwfT4/PbvcWsaxSfpC5LWl+L4UKnuAkmPS3pM0vGd6XVrWolRUo+kV0rl3+hcz5vX2++ppM+ksVol6cul8qzGsZX4qjSGkm4uxbFa0oOluqzGEFqLsWLjeLCk+1IcSyUdnsol6ao0jg9JOrRzPW+ziOj4AzgaOBR4uK58MsUiuTXA+FR2LHB7p/s82PiALwDnNWh7ILACGAVMBZ4ARnQ6hjbH2FM/1jk8eonxvcCPgVHp+Z65jmOL8VVmDOvq5wGfz3UMBxBjZcYRuAs4IW1/CFhc2r6T4tolRwBLOt3/dj264gg8Iu4BnmtQdSVwPg0uEJOTPuJr5ETgOxHxakT8Gnic4rK1Xa3FGLPUS4yfBuZGxKupzTOpPLtxbDG+LPX1eypJwJ8DN6Wi7MYQWo4xS73EGMCYtL0bb16X5ETgW1G4DxgracL26enQ6ooE3oikPwPWR8SKBtVHSloh6U5J07d339podprSWShpXCrbG1hbarMuleWqUYwAUyX9QtLdkmZ2rHeDNw2YKWlJiuXdqbwq49hbfFCdMayZCWyMiF+l51UZw7L6GKE64/hZ4B8lrQUuBy5I5VUcR6BLE7iktwEXAY3Oby+nuE7sQcDVwK3bs29t9HVgP+BgYAPFtBY0eYnaTPQW4wZgSkQcAnwO+LakMY1fouvtCIyjmJr7O+C76SinKuPYW3xVGsOaT7LtkWlVxrCsPsYqjeOngXMjYjJwLnBdKq/iOAJdmsAp/uhPBVZIWk1xmdblkt4RES9GxBaAiPgBsJPSArecRMTGiHgjIrYC/8SbU3OVuURtbzGmKcln0/YyinOL0zrX00FZB9ySpufuB7ZS3EihKuPYML6KjSGSdgQ+DtxcKq7KGAKNY6zYOJ4O3JK2/40K/k2t15UJPCJWRsSeEdETET0UA3BoRDwt6R3pCIC0ynAH4NkOdndA6s7BfAyoraZcBJwiaZSkqcA7gfu3d//aobcYJe2h4n7ySNqXIsYnt38P2+JW4DgASdOAkRR3QarKODaMr2JjCPB+4JcRsa5UVpUxrHlLjBUbx6eAY9L2cUDtNMEi4LS0Gv0IYHNEbOhEB9uu06voolgleBPFVM5rFMn6rLr61by5Cn02sIpideh9wFGd7v9A4gNuBFYCD1H8gk0otb+I4pPwY6RVld3+aCVG4BOlMVwO/Gmn+z+IGEcC/0Lx4WQ5cFyu49hKfFUaw1R+A/BXDdpnNYatxlilcQT+GFiWYlkCHJbaCrgmjeNKYEan+9+uhy+lamZmlqGunEI3MzOzvjmBm5mZZcgJ3MzMLENO4GZmZhlyAjczM8uQE7iZmVmGnMDNzMwy5ARuZmaWISdwMzOzDDmBm5mZZcgJ3MzMLENO4GZmZhlyAjfLkKQt6faPfbXpkRTpPtBmVjFO4GZdTtJiSWeXyyJi14joivs2S1og6TFJWyWd0aD+XElPS9osaaGkUaW6Hkk/kfSypF9Ken+z+5oNd07gZl1KhRz+j64A/priftLbkHQ8MAd4H9AD7Av8Q6nJTcAvgLdT3Hv7e5L2aHJfs2Ethz8OZlmRtFrSBZIekfS8pOsl7Zzqxkm6XdJvU93tkiaV9l0s6TJJPwNeBm4EZgLz07T5/NQuJO2ftneRNE/SmnSkeq+kXRr0azdJ10naIGm9pEsljUh1+0u6O+2/SdLNzcYbEddExP8Cft+g+nTguohYFRHPA18EzkjvOQ04FLgkIl6JiO8DK4FP9LevmTmBmw2VTwHHA/sB04CLU/kOwPXAPsAU4BVgft2+pwKzgNEUCeunwOw0bT67wXtdDhwGHAXsDpwPbG3Q7pvA68D+wCHAB4Da1PwXgbuAccAk4OraTulDxpzmwn6L6RRH6DUrgL0kvT3VPRkRL9XVT29iX7Nhz4tbzIbG/IhYCyDpMoqEeHFEPAt8v9Yo1f2kbt8bImJVqU2vb5Km2M8EjoiI9an45/X7SdoLOAEYGxGvAL+TdCXFB4VrgdcoPlRMjIh1wL21fSPiIy1Fvq1dgc2l57Xt0Q3qavV7N7Hvs4Pok1kl+AjcbGisLW2vASYCSHqbpGvTdPeLwD3A2NpUdoN9+zMe2Bl4op92+wA7ARskvSDpBYrEvWeqPx8QcL+kVZLObKEPfdkCjCk9r22/1KCuVl87Iu9rX7NhzwncbGhMLm1PAZ5K238LHAC8JyLGAEen8vJhdtS9Vv3zsk0U557366c/a4FXgfERMTY9xkTEdICIeDoi/iIiJgJ/CXytdo59kFYBB5WeHwRsTDMRq4B9JY2uq1/VxL5mw54TuNnQ+BtJkyTtDlwI1BaFjaY47/1CqrukidfaSLEC+y0iYiuwELhC0kRJIyQdWf91q4jYQHGOe56kMZJ2kLSfpGMAJJ1cWkz3PMWHhjeaCVTSyLRIT8BOknYurZ7/FnCWpAMljaNYC3BD6tN/Ag8Cl6R9Pga8izdPMfS6r5k5gZsNlW9TJMwn0+PSVP4VYBeKI+f7gB828VpfBU5Kq9avalB/HsXq7QeA54Av0fj/9mnASOARiiT9PWBCqns3sETSFmARcE5E/BpA0p2SLuyjf3dRfCg5CliQto8GiIgfAl+mOM+/Jj3KH1pOAWak/swFToqI3za5r9mwpoi+ZufMrFWSVgNnR8SPO90XM6suH4GbmZllyAnczMwsQ/0mcEmT07WKH01fLzknle8u6T8k/Sr9HDf03TXrfhHR4+lzMxtq/Z4DlzQBmBARy9PXPZYBH6W4QtRzETE3XaVpXET8/VB32MzMzAawiE3SbRSXfpwPHBsRG1KSXxwRB/S17/jx46Onp2egfTUzM6u8ZcuWbYqIPfpr19KlVCX1UFxDeQmwV/puKSmJ79nLPrMoLtfIlClTWLp0aStvaWZmNqxIWtNMu6YTuKRdKS6w8NmIeLGv6zOXRcQCiu+GMmPGDH9nzayDeubc0W+b1XM/3HLboXztVvthNlw0tQpd0k4UyftfI+KWVLwxTZ3XzpM/MzRdNDMzs3rNrEIXcB3waERcUapaRHG/XtLP29rfPTMzM2ukmSn0P6K4P/FKSQ+msgspLnv4XUlnAb8BTh6aLpqZmVm9fhN4RNzLtndKKntfe7tjZmZmzfCV2MzMzDLkBG5mZpYhJ3AzM7MMOYGbmZllyAnczMwsQ07gZmZmGXICNzMzy5ATuJmZWYacwM3MzDLkBG5mZpYhJ3AzM7MMOYGbmZllyAnczMwsQ07gZmZmGXICNzMzy5ATuJmZWYacwM3MzDLkBG5mZpYhJ3AzM7MMOYGbmZllyAnczMwsQ07gZmZmGXICNzMzy5ATuJmZWYacwM3MzDLkBG5mZpYhJ3AzM7MMOYGbmZllaFAJXNIHJT0m6XFJc9rVKTMzM+vbjgPdUdII4BrgT4B1wAOSFkXEI+3qnNlw1DPnjn7brJ774e3QEzPrZoM5Aj8ceDwinoyIPwDfAU5sT7fMzMysL4qIge0onQR8MCLOTs9PBd4TEbPr2s0CZqWnBwCPDeDtxgObBtTRfFQ9xqrHB46xKqoeY9Xjg/xj3Cci9uiv0YCn0AE1KHvLp4GIWAAsGMT7IGlpRMwYzGt0u6rHWPX4wDFWRdVjrHp8MDxihMFNoa8DJpeeTwKeGlx3zMzMrBmDSeAPAO+UNFXSSOAUYFF7umVmZmZ9GfAUekS8Lmk28CNgBLAwIla1rWfbGtQUfCaqHmPV4wPHWBVVj7Hq8cHwiHHgi9jMzMysc3wlNjMzsww5gZuZmWWoKxK4pIWSnpH0cIO68ySFpPHp+bGSNkt6MD0+v/173JpG8Un6gqT1pTg+VKq7IF2e9jFJx3em161pJUZJPZJeKZV/o3M9b15vv6eSPpPGapWkL5fKsxrHVuKr0hhKurkUx2pJD5bqshpDaC3Gio3jwZLuS3EslXR4Kpekq9I4PiTp0M71vM0iouMP4GjgUODhuvLJFIvk1gDjU9mxwO2d7vNg4wO+AJzXoO2BwApgFDAVeAIY0ekY2hxjT/1Y5/DoJcb3Aj8GRqXne+Y6ji3GV5kxrKufB3w+1zEcQIyVGUfgLuCEtP0hYHFp+06Ka5ccASzpdP/b9eiKI/CIuAd4rkHVlcD5NLhATE76iK+RE4HvRMSrEfFr4HGKy9Z2tRZjzFIvMX4amBsRr6Y2z6Ty7Maxxfiy1NfvqSQBfw7clIqyG0NoOcYs9RJjAGPS9m68eV2SE4FvReE+YKykCdunp0OrKxJ4I5L+DFgfESsaVB8paYWkOyVN3959a6PZaUpnoaRxqWxvYG2pzbpUlqtGMQJMlfQLSXdLmtmx3g3eNGCmpCUplnen8qqMY2/xQXXGsGYmsDEifpWeV2UMy+pjhOqM42eBf5S0FrgcuCCVV3EcgS5N4JLeBlwENDq/vZziOrEHAVcDt24fWIWNAAAgAElEQVTPvrXR14H9gIOBDRTTWtDkJWoz0VuMG4ApEXEI8Dng25LGNH6JrrcjMI5iau7vgO+mo5yqjGNv8VVpDGs+ybZHplUZw7L6GKs0jp8Gzo2IycC5wHWpvIrjCHRpAqf4oz8VWCFpNcVlWpdLekdEvBgRWwAi4gfATkoL3HISERsj4o2I2Ar8E29OzVXmErW9xZimJJ9N28sozi1O61xPB2UdcEuanrsf2EpxI4WqjGPD+Co2hkjaEfg4cHOpuCpjCDSOsWLjeDpwS9r+Nyr4N7VeVybwiFgZEXtGRE9E9FAMwKER8bSkd6QjANIqwx2AZzvY3QGpOwfzMaC2mnIRcIqkUZKmAu8E7t/e/WuH3mKUtIeK+8kjaV+KGJ/c/j1si1uB4wAkTQNGUtwFqSrj2DC+io0hwPuBX0bEulJZVcaw5i0xVmwcnwKOSdvHAbXTBIuA09Jq9COAzRGxoRMdbLtOr6KLYpXgTRRTOa9RJOuz6upX8+Yq9NnAKorVofcBR3W6/wOJD7gRWAk8RPELNqHU/iKKT8KPkVZVdvujlRiBT5TGcDnwp53u/yBiHAn8C8WHk+XAcbmOYyvxVWkMU/kNwF81aJ/VGLYaY5XGEfhjYFmKZQlwWGor4Jo0jiuBGZ3uf7sevpSqmZlZhrpyCt3MzMz65gRuZmaWISdwMzOzDDmBm5mZZcgJ3MzMLENO4GZmZhlyAjczM8uQE7iZmVmGnMDNzMwy5ARuZmaWISdwMzOzDDmBm5mZZcgJ3CxDkrak2z/21aZHUqT7QJtZxTiBm3U5SYslnV0ui4hdI6Lj922WNE3SbZJ+K+k5ST+SdEBdm3MlPS1ps6SFkkaV6nok/UTSy5J+Ken9ze5rNtw5gZt1KRW6/f/oWIp7vR8A7AXcD9xWq5R0PDAHeB/QA+wL/ENp/5uAXwBvp7j39vck7dHkvmbDWrf/cTDLjqTVki6Q9Iik5yVdL2nnVDdO0u3piPX5tD2ptO9iSZdJ+hnwMnAjMBOYn6bN56d2IWn/tL2LpHmS1qQj1Xsl7dKgX7tJuk7SBknrJV0qaUSq21/S3Wn/TZJubibWiLg/Iq6LiOci4jXgSuAASW9PTU4HrouIVRHxPPBF4Iz0ntOAQ4FLIuKViPg+sBL4RH/7mpkTuNlQ+RRwPLAfMA24OJXvAFwP7ANMAV4B5tfteyowCxhNkbB+CsxO0+azG7zX5cBhwFHA7sD5wNYG7b4JvA7sDxwCfACoTc1/EbgLGAdMAq6u7ZQ+ZMxpLmyOBp6OiGfT8+nAilL9CmCvlOCnA09GxEt19dOb2Nds2PPiFrOhMT8i1gJIuowiIV6cEtv3a41S3U/q9r0hIlaV2vT6JmmK/UzgiIhYn4p/Xr+fpL2AE4CxEfEK8DtJV1J8ULgWeI3iQ8XEiFgH3FvbNyI+0kzAaSbhGuBzpeJdgc2l57Xt0Q3qavV7N7Hvs5gNcz4CNxsaa0vba4CJAJLeJunaNN39InAPMLY2ld1g3/6MB3YGnuin3T7ATsAGSS9IeoEice+Z6s8HBNwvaZWkM1voA+m89V3A1yLiplLVFmBM6Xlt+6UGdbX62hF5X/uaDXtO4GZDY3JpewrwVNr+W4oFX++JiDEUU85QJM+aqHut+udlm4DfU0zV92Ut8CowPiLGpseYiJgOEBFPR8RfRMRE4C+Br9XOsfdH0jiK5L0oIi6rq14FHFR6fhCwMc1ErAL2lTS6rn5VE/uaDXtO4GZD428kTZK0O3AhUFsUNprivPcLqe6SJl5rI8UK7LeIiK3AQuAKSRMljZB0ZP3XrSJiA0WSnSdpjKQdJO0n6RgASSeXFtM9T/Gh4Y3+OiZpDPAj4GcR0eg8+beAsyQdmBL9xcANqU//CTwIXCJpZ0kfA97Fm6cYet3XzJzAzYbKtykS5pPpcWkq/wqwC8WR833AD5t4ra8CJ6VV61c1qD+PYvX2A8BzwJdo/H/7NGAk8AhFkv4eMCHVvRtYImkLxdfCzomIXwNIulPShb307WNp3/+WVsnXHlMAIuKHwJcpzvOvSY/yh5ZTgBmpP3OBkyLit03uazasKaKv2Tkza5Wk1cDZEfHjTvfFzKrLR+BmZmYZcgI3MzPLUL8JXNLkdK3iR9PXS85J5btL+g9Jv0o/xw19d826X0T0ePrczIZav+fAJU0AJkTE8vR1j2XARymuEPVcRMxNV2kaFxF/P9QdNjMzswEsYpN0G8WlH+cDx0bEhpTkF0fEAX3tO378+Ojp6RloX83MzCpv2bJlmyJij/7atXQpVUk9FNdQXgLslb5bSkrie/ayzyyKyzUyZcoUli5d2spbmpmZDSuS1jTTrukELmlXigssfDYiXuzr+sxlEbEAWAAwY8YMf2fNrIN65tzRb5vVcz/cctuhfO1W+2E2XDS1Cl3SThTJ+18j4pZUvDFNndfOkz8zNF00MzOzes2sQhdwHfBoRFxRqlpEcb9e0s/b2t89MzMza6SZKfQ/org/8UpJD6ayCykue/hdSWcBvwFOHpoumpmZWb1+E3hE3Mu2d0oqe197u2NmZmbN8JXYzMzMMuQEbmZmliEncDMzsww5gZuZmWXICdzMzCxDTuBmZmYZcgI3MzPLkBO4mZlZhpzAzczMMuQEbmZmliEncDMzsww5gZuZmWXICdzMzCxDTuBmZmYZcgI3MzPLkBO4mZlZhpzAzczMMuQEbmZmliEncDMzsww5gZuZmWXICdzMzCxDTuBmZmYZcgI3MzPLkBO4mZlZhpzAzczMMuQEbmZmliEncDMzsww5gZuZmWVoUAlc0gclPSbpcUlz2tUpMzMz69uOA91R0gjgGuBPgHXAA5IWRcQj7eqc2XDUM+eOftusnvvh7dATM+tmgzkCPxx4PCKejIg/AN8BTmxPt8zMzKwvioiB7SidBHwwIs5Oz08F3hMRs+vazQJmpacHAI8N4O3GA5sG1NF8VD3GqscHjrEqqh5j1eOD/GPcJyL26K/RgKfQATUoe8ungYhYACwYxPsgaWlEzBjMa3S7qsdY9fjAMVZF1WOsenwwPGKEwU2hrwMml55PAp4aXHfMzMysGYNJ4A8A75Q0VdJI4BRgUXu6ZWZmZn0Z8BR6RLwuaTbwI2AEsDAiVrWtZ9sa1BR8JqoeY9XjA8dYFVWPserxwfCIceCL2MzMzKxzfCU2MzOzDDmBm5mZZagrErikhZKekfRwg7rzJIWk8en5sZI2S3owPT6//XvcmkbxSfqCpPWlOD5UqrsgXZ72MUnHd6bXrWklRkk9kl4plX+jcz1vXm+/p5I+k8ZqlaQvl8qzGsdW4qvSGEq6uRTHakkPluqyGkNoLcaKjePBku5LcSyVdHgql6Sr0jg+JOnQzvW8zSKi4w/gaOBQ4OG68skUi+TWAONT2bHA7Z3u82DjA74AnNeg7YHACmAUMBV4AhjR6RjaHGNP/Vjn8OglxvcCPwZGped75jqOLcZXmTGsq58HfD7XMRxAjJUZR+Au4IS0/SFgcWn7ToprlxwBLOl0/9v16Ioj8Ii4B3iuQdWVwPk0uEBMTvqIr5ETge9ExKsR8WvgcYrL1na1FmPMUi8xfhqYGxGvpjbPpPLsxrHF+LLU1++pJAF/DtyUirIbQ2g5xiz1EmMAY9L2brx5XZITgW9F4T5grKQJ26enQ6srEngjkv4MWB8RKxpUHylphaQ7JU3f3n1ro9lpSmehpHGpbG9gbanNulSWq0YxAkyV9AtJd0ua2bHeDd40YKakJSmWd6fyqoxjb/FBdcawZiawMSJ+lZ5XZQzL6mOE6ozjZ4F/lLQWuBy4IJVXcRyBLk3gkt4GXAQ0Or+9nOI6sQcBVwO3bs++tdHXgf2Ag4ENFNNa0OQlajPRW4wbgCkRcQjwOeDbksY0fomutyMwjmJq7u+A76ajnKqMY2/xVWkMaz7JtkemVRnDsvoYqzSOnwbOjYjJwLnAdam8iuMIdGkCp/ijPxVYIWk1xWVal0t6R0S8GBFbACLiB8BOSgvcchIRGyPijYjYCvwTb07NVeYStb3FmKYkn03byyjOLU7rXE8HZR1wS5qeux/YSnEjhaqMY8P4KjaGSNoR+Dhwc6m4KmMINI6xYuN4OnBL2v43Kvg3tV5XJvCIWBkRe0ZET0T0UAzAoRHxtKR3pCMA0irDHYBnO9jdAak7B/MxoLaachFwiqRRkqYC7wTu3979a4feYpS0h4r7ySNpX4oYn9z+PWyLW4HjACRNA0ZS3AWpKuPYML6KjSHA+4FfRsS6UllVxrDmLTFWbByfAo5J28cBtdMEi4DT0mr0I4DNEbGhEx1su06vootileBNFFM5r1Ek67Pq6lfz5ir02cAqitWh9wFHdbr/A4kPuBFYCTxE8Qs2odT+IopPwo+RVlV2+6OVGIFPlMZwOfCnne7/IGIcCfwLxYeT5cBxuY5jK/FVaQxT+Q3AXzVon9UYthpjlcYR+GNgWYplCXBYaivgmjSOK4EZne5/ux6+lKqZmVmGunIK3czMzPrmBG5mZpYhJ3AzM7MMOYGbmZllyAnczMwsQ07gZmZmGXICNzMzy5ATuJmZWYacwM3MzDLkBG5mZpYhJ3AzM7MMOYGbmZllyAncLEOStqTbP/bVpkdSpPtAm1nFOIGbdTlJiyWdXS6LiF0joqvu2yzp9PSB4exSmSR9SdKz6fFlSSrVHyxpmaSX08+Dm93XbLhzAjfrUimBZfF/VNI44AKKe0uXzQI+ChwEvAv4CPCXaZ+RwG0U9xsfB3wTuC2V97mvmTmBm7WdpNWSLpD0iKTnJV0vaedUN07S7ZJ+m+pulzSptO9iSZdJ+hnwMnAjMBOYn6bN56d2IWn/tL2LpHmS1kjaLOleSbs06Ndukq6TtEHSekmXShqR6vaXdHfaf5Okm1sM+/8BrgI21ZWfDsyLiHURsR6YB5yR6o4FdgS+EhGvRsRVgIDjmtjXbNhzAjcbGp8Cjgf2A6YBF6fyHYDrgX2AKcArwPy6fU+lOPocTZGwfgrMTtPmsxu81+XAYcBRwO7A+cDWBu2+CbwO7A8cAnwAqE13fxG4i+JIeBJwdW2n9CFjTm+BSjocmAF8o0H1dGBF6fmKVFareygiolT/UF19b/uaDXte3GI2NOZHxFoASZdRJMSLI+JZ4Pu1RqnuJ3X73hARq0pten2TNMV+JnBEOkoF+Hn9fpL2Ak4AxkbEK8DvJF1J8UHhWuA1ig8VEyNiHXBvbd+I+Egf7z8C+BrwmYjY2qCvuwKbS883A7umc9n1dbX60f3tW5f0zYYlH4GbDY21pe01wEQASW+TdG2a7n4RuAcYW5vKbrBvf8YDOwNP9NNuH2AnYIOkFyS9QJG490z151NMX98vaZWkM5t8/7+mOIr+f3up3wKMKT0fA2xJCbi+rlb/UhP7mg17TuBmQ2NyaXsK8FTa/lvgAOA9ETEGODqVlw9d6xNUXwlrE/B7iqn6vqwFXgXGR8TY9BgTEdMBIuLpiPiLiJhIsVDsa7Vz7P14H/AxSU9LeppiGn9e7Vw9xaK2g0rtD+LNhW6rgHfVrSx/V119b/uaDXtO4GZD428kTZK0O3AhUFsUNprivPcLqe6SJl5rI9DwO98RsRVYCFwhaaKkEZKOlDSqrt0GinPc8ySNkbSDpP0kHQMg6eTSYrrnKT40vNFE384A/gtwcHosBf4BuCjVfwv4nKS9JU2k+ABzQ6pbnN7jv0saJal2fv9/N7Gv2bDnBG42NL5NkTCfTI9LU/lXgF0ojpzvA37YxGt9FTgprVq/qkH9ecBK4AHgOeBLNP6/fRowEniEIkl/D5iQ6t4NLJG0BVgEnBMRvwaQdKekCxt1LCJeSEfvT0fE08AfgBcjonbu+lrg31P/HgbuSGVExB8oviZ2GvACxbn8j6byPvc1M5BPJ5m1l6TVwNkR8eNO98XMqstH4GZmZhlyAjczM8tQvwlc0mRJP5H0aPp6yTmpfHdJ/yHpV+nnuKHvrln3i4geT5+b2VDr9xy4pAnAhIhYLmk0sIxi4ckZwHMRMTddpWlcRPz9UHfYzMzMBrCITdJtFJd+nA8cGxEbUpJfHBEH9LXv+PHjo6enZ6B9NTMzq7xly5Ztiog9+mvX0qVUJfVQXEN5CbBX+m4pKYnv2cs+sygu18iUKVNYunRpK29pZmY2rEha00y7phO4pF0pruH82Yh4sdnb8kbEAmABwIwZM/ydNbMO6plzR79tVs/98JC/9lC1NRtOmkrgknaiSN7/GhG3pOKNkiaUptCfGapOmtn258Rp1t2aWYUu4Drg0Yi4olS1iOJ+vaSft7W/e2ZmZtZIM0fgf0Rxf+KVkh5MZRcCc4HvSjoL+A1w8tB00czMzOr1m8Aj4l62vVNS2fva2x0zMzNrhq/EZmZmliEncDMzsww5gZuZmWXICdzMzCxDTuBmZmYZcgI3MzPLkBO4mZlZhpzAzczMMuQEbmZmliEncDMzsww5gZuZmWXICdzMzCxDTuBmZmYZcgI3MzPLkBO4mZlZhpzAzczMMuQEbmZmliEncDMzsww5gZuZmWXICdzMzCxDTuBmZmYZcgI3MzPLkBO4mZlZhpzAzczMMuQEbmZmliEncDMzsww5gZuZmWXICdzMzCxDg0rgkj4o6TFJj0ua065OmZmZWd92HOiOkkYA1wB/AqwDHpC0KCIeaVfnzIajnjl39Ntm9dwPb4eemFk3G8wR+OHA4xHxZET8AfgOcGJ7umVmZmZ9GfAROLA3sLb0fB3wnsF1x7pVM0eF4CPD3gzlUbWP2Lflfw8bLhQRA9tROhk4PiLOTs9PBQ6PiM/UtZsFzEpPDwAeG8DbjQc2Daij+ah6jFWPDxxjVVQ9xqrHB/nHuE9E7NFfo8Ecga8DJpeeTwKeqm8UEQuABYN4HyQtjYgZg3mNblf1GKseHzjGqqh6jFWPD4ZHjDC4c+APAO+UNFXSSOAUYFF7umVmZmZ9GfAReES8Lmk28CNgBLAwIla1rWdmZmbWq8FMoRMRPwB+0Ka+9GVQU/CZqHqMVY8PHGNVVD3GqscHwyPGgS9iMzMzs87xpVTNzMwy1BUJXNJCSc9IerhB3XmSQtL49PxYSZslPZgen9/+PW5No/gkfUHS+lIcHyrVXZAuT/uYpOM70+vWtBKjpB5Jr5TKv9G5njevt99TSZ9JY7VK0pdL5VmNYyvxVWkMJd1cimO1pAdLdVmNIbQWY8XG8WBJ96U4lko6PJVL0lVpHB+SdGjnet5mEdHxB3A0cCjwcF35ZIpFcmuA8ansWOD2Tvd5sPEBXwDOa9D2QGAFMAqYCjwBjOh0DG2Osad+rHN49BLje4EfA6PS8z1zHccW46vMGNbVzwM+n+sYDiDGyowjcBdwQtr+ELC4tH0nIOAIYEmn+9+uR1ccgUfEPcBzDaquBM4Hsj5R30d8jZwIfCciXo2IXwOPU1y2tqu1GGOWeonx08DciHg1tXkmlWc3ji3Gl6W+fk8lCfhz4KZUlN0YQssxZqmXGAMYk7Z3483rkpwIfCsK9wFjJU3YPj0dWl2RwBuR9GfA+ohY0aD6SEkrJN0pafr27lsbzU5TOgsljUtljS5Ru/f271rbNIoRYKqkX0i6W9LMjvVu8KYBMyUtSbG8O5VXZRx7iw+qM4Y1M4GNEfGr9LwqY1hWHyNUZxw/C/yjpLXA5cAFqbyK4wh0aQKX9DbgIqDR+e3lFJeZOwi4Grh1e/atjb4O7AccDGygmNaCYpqnXq4zEL3FuAGYEhGHAJ8Dvi1pTOOX6Ho7AuMopub+DvhuOsqpyjj2Fl+VxrDmk2x7ZFqVMSyrj7FK4/hp4NyImAycC1yXyqs4jkCXJnCKP/pTgRWSVlNcpnW5pHdExIsRsQX+/++h76S0wC0nEbExIt6IiK3AP/Hm1FxTl6jNQW8xpinJZ9P2Mopzi9M619NBWQfckqbn7ge2UlyHuSrj2DC+io0hknYEPg7cXCquyhgCjWOs2DieDtyStv+NCv5NrdeVCTwiVkbEnhHRExE9FANwaEQ8Lekd6QiAtMpwB+DZDnZ3QOrOwXwMqK2mXAScImmUpKnAO4H7t3f/2qG3GCXtoeJ+8kjalyLGJ7d/D9viVuA4AEnTgJEUN1Goyjg2jK9iYwjwfuCXEbGuVFaVMax5S4wVG8engGPS9nFA7TTBIuC0tBr9CGBzRGzoRAfbrtOr6KJYJXgTxVTOaxTJ+qy6+tW8uQp9NrCKYnXofcBRne7/QOIDbgRWAg9R/IJNKLW/iOKT8GOkVZXd/mglRuATpTFcDvxpp/s/iBhHAv9C8eFkOXBcruPYSnxVGsNUfgPwVw3aZzWGrcZYpXEE/hhYlmJZAhyW2gq4Jo3jSmBGp/vfroevxGZmZpahrpxCNzMzs745gZuZmWXICdzMzCxDTuBmZmYZcgI3MzPLkBO4mZlZhpzAzczMMuQEbmZmliEncDMzsww5gZuZmWXICdzMzCxDTuBmZmYZcgI3y4ykLenWj3216ZEU6R7QZlZBTuBmXUzSYklnl8siYteI6Ip7Nks6TtJySS9KelLSrLr6/yppjaTfSbpV0u6lut0l/c9Ut0bSf212XzNzAjfrSip09f9PSTsB/xO4FtgN+L+BKyQdlOqnp7pTgb2Al4GvlV7iGuAPqe5TwNfTPs3sazbsdfUfCLPcSFot6QJJj0h6XtL1knZOdeMk3S7pt6nudkmTSvsulnSZpJ9RJKwbgZnA/DRtPj+1C0n7p+1dJM1LR6qbJd0raZcG/dpN0nWSNkhaL+lSSSNS3f6S7k77b5J0c5Ph7g6MAW6MwgPAo8CBqf5TwL9HxD0RsQX4H8DHJY2W9H8AnwD+R0RsiYh7gUUUCbvPfZvsm1nlOYGbtd+ngOOB/YBpwMWpfAfgemAfYArwCjC/bt9TgVnAaOAM4KfA7DRtPrvBe10OHAYcRZFQzwe2Nmj3TeB1YH/gEOADQG1q/ovAXcA4YBJwdW2n9CFjTqMgI2IjcBPw3ySNkHRkiu3e1GQ6sKLU/gmKI+5p6fFGRPxn6SVXpH3629fMAC9wMWu/+RGxFkDSZRQJ8eKIeBb4fq1RqvtJ3b43RMSqUpte3yRNsZ8JHBER61Pxz/+/9u41Vq6qDOP4/6GlVSS1NYdCKy1TScGACYjlolhANJGr9QIGPggJIFKsARQRSoKE+KFeEEFRqdAYlavaQCNWlEQgfmihLZRSsHIrtlDaANICQYL09cNeJ50Oc87MnJlzZtY+zy+ZdM/ae89Zb9/TvrP2XrOm9jxJewLHAxMj4k3gDUnXULxRuAF4m6LwTo2IjewowETESQ1ivRW4Ebg2PZ/bHzuwO7C15vitFG9O3hlkX6NzzQyPwM2Gw4aq7eeAqQCSdpN0Q7rcvQ14AJjYfym7zrmN9AHvAZ5ucNw+wK7AJkmvSnqVonBPTvsvAQQ8KGmtpLOa+eGSPgzcDpwBjKMYNV8i6cR0yOsUl9irTQBea7Cv0blmhgu42XCYVrU9HXghbX8L2B84PCImAEel9uphdtS8Vu3zai8B/6W4VD+YDcBbQF9ETEyPCRFxIEBEvBgRX42IqcDXgJ/332Nv4CPAuoi4JyK2R8Q64G6K0T7AWuCg/oPTR9/GA/9Kj7GSZla93kHpnEbnmhku4GbD4euS9k4fe5pPMUqF4vLvm8Crad93m3itzUDdz3xHxHZgEcXM76n996Elja85bhPFPe6rJU2QtIukfSUdDSDp1KrJdP+heNPwThN9exiYmT5KJkn7Aiex4971zcDJkmanSWtXAYsj4rWIeANYDFwl6X2SjgTmUEzcG/TcJvplNiq4gJt13i0UBfOZ9Pheav8J8F6KkfMy4C9NvNa1wClp1vp1dfZfDKwBHgJeAb5P/X/X/Ze5H6co0n8ApqR9hwLLJb1OMRP8goh4FkDSUknz63UsTSw7C7gO2AbcT3GP/6a0fy1wHkUx3kLxBub8qpc4n+LvYwvFvfS5/ff/mzjXbNRTxGBX6MysFZLWA+dExL3d7ouZlZtH4GZmZhlyATczM8tQwwIuaZqkv0t6In3E5ILU/gFJf5P0ZPpz0vB316y3RUTFl8/NbCQ0vAcuaQowJSJWpWUMVwKfp1gl6pWIWJBWapoUEd8Z7g6bmZnZECaxSbqLYvnHnwHHRMSmVOTvi4j9Bzu3r68vKpXKUPtqZmZWeitXrnwpIvZodFxLS6lKqlCso7wc2DN9vpRUxCcPcioAlUqFFStWtPIjzczMRhVJzzVzXNMFXNLuFJ/xvDAitg22RnPNeedSrLnM9OnTm/1xZpaZyqV3Nzxm/YITh/VYs9GkqQKevvf3j8DNEbE4NW+WNKXqEvqWeudGxEJgIcCsWbP8oXOzLnLhNCuPZmahi2JlpSci4sdVu5YAZ6btM4G7Ot89MzMzq6eZEfiRFN9RvEbSI6ltPrAAuEPS2cC/gVOHp4tmZmZWq2EBj4h/sPO3JVX7dGe7Y2ZmZs3wSmxmZmYZcgE3MzPLkAu4mZlZhlzAzczMMuQCbmZmliEXcDMzswy5gJuZmWXIBdzMzCxDLuBmZmYZcgE3MzPLkAu4mZlZhlzAzczMMuQCbmZmliEXcDMzswy5gJuZmWXIBdzMzCxDLuBmZmYZcgE3MzPLkAu4mZlZhlzAzczMMuQCbmZmliEXcDMzswy5gJuZmWXIBdzMzCxDLuBmZmYZcgE3MzPLkAu4mZlZhlzAzczMMuQCbmZmlqGx7Zws6TjgWmAMcGNELOhIr8xGscqldzc8Zv2CE0egJ2bWy4Y8Apc0BrgeOB44ADhd0gGd6piZmZkNrJ0R+GHAUxHxDICk24A5wOOd6Jj1lmZGheCR4UCGc1TtEfvO/Pdho4UiYmgnSqcAx0XEOen5V4DDI2JezXHnAuemp/sD64bw4/qAl4bU0XyUPcayxweOsSzKHmPZ44P8Y9wnIvZodEcddUcAAASZSURBVFA7I3DVaXvXu4GIWAgsbOPnIGlFRMxq5zV6XdljLHt84BjLouwxlj0+GB0xQnuz0DcC06qe7w280F53zMzMrBntFPCHgJmSZkgaB5wGLOlMt8zMzGwwQ76EHhH/kzQPuIfiY2SLImJtx3q2s7YuwWei7DGWPT5wjGVR9hjLHh+MjhiHPonNzMzMuscrsZmZmWXIBdzMzCxDPVHAJS2StEXSY3X2XSwpJPWl58dI2irpkfS4YuR73Jp68Um6UtLzVXGcULXvMklPSVon6bPd6XVrWolRUkXSm1Xtv+xez5s30O+ppG+kXK2V9IOq9qzy2Ep8ZcqhpNur4lgv6ZGqfVnlEFqLsWR5PFjSshTHCkmHpXZJui7l8VFJh3Sv5x0WEV1/AEcBhwCP1bRPo5gk9xzQl9qOAf7U7T63Gx9wJXBxnWMPAFYD44EZwNPAmG7H0OEYK7W5zuExQIyfAu4Fxqfnk3PNY4vxlSaHNfuvBq7INYdDiLE0eQT+Chyftk8A7qvaXkqxdskRwPJu979Tj54YgUfEA8ArdXZdA1xCnQVicjJIfPXMAW6LiLci4lngKYpla3taizFmaYAY5wILIuKtdMyW1J5dHluML0uD/Z5KEvBl4NbUlF0OoeUYszRAjAFMSNvvZ8e6JHOA30RhGTBR0pSR6enw6okCXo+kzwHPR8TqOrs/Lmm1pKWSDhzpvnXQvHRJZ5GkSantg8CGqmM2prZc1YsRYIakhyXdL2l213rXvv2A2ZKWp1gOTe1lyeNA8UF5cthvNrA5Ip5Mz8uSw2q1MUJ58ngh8ENJG4AfAZel9jLmEejRAi5pN+ByoN797VUU68QeBPwUuHMk+9ZBvwD2BQ4GNlFc1oIml6jNxEAxbgKmR8RHgW8Ct0iaUP8let5YYBLFpblvA3ekUU5Z8jhQfGXKYb/T2XlkWpYcVquNsUx5nAtcFBHTgIuAm1J7GfMI9GgBp/hPfwawWtJ6imVaV0naKyK2RcTrABHxZ2BXpQluOYmIzRHxTkRsB37FjktzpVmidqAY0yXJl9P2Sop7i/t1r6dt2QgsTpfnHgS2U3yRQlnyWDe+kuUQSWOBLwK3VzWXJYdA/RhLlsczgcVp+/eU8P/UWj1ZwCNiTURMjohKRFQoEnBIRLwoaa80AiDNMtwFeLmL3R2SmnswXwD6Z1MuAU6TNF7SDGAm8OBI968TBopR0h4qvk8eSR+iiPGZke9hR9wJHAsgaT9gHMW3IJUlj3XjK1kOAT4D/DMiNla1lSWH/d4VY8ny+AJwdNo+Fui/TbAEOCPNRj8C2BoRm7rRwY7r9iy6KGYJ3kpxKedtimJ9ds3+9eyYhT4PWEsxO3QZ8Ilu938o8QG/BdYAj1L8gk2pOv5yinfC60izKnv90UqMwJeqcrgKOLnb/W8jxnHA7yjenKwCjs01j63EV6YcpvZfA+fVOT6rHLYaY5nyCHwSWJliWQ58LB0r4PqUxzXArG73v1MPL6VqZmaWoZ68hG5mZmaDcwE3MzPLkAu4mZlZhlzAzczMMuQCbmZmliEXcDMzswy5gJuZmWXo/6Uk22oF85YtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f10628c8d30>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(particle_numbers.shape[0], 1, figsize=(7, 7))\n",
    "for n, _ax in enumerate(ax):\n",
    "    _ax.bar(lines[n, :], totals[n, :] / totals[n, :].sum() * 100);\n",
    "    _ax.set_title(\"particles: \" + str(particle_numbers[n, 0]))\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAHwCAYAAABZrD3mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmYXVWZ7/Hvj4DMCJgKEgiUQBwAJa0xthcHVERQFFBAxsapg92NXC62Emm6oRG7UVttEKeoDCoEUAKCenHAizSDQ4WOEEQUMEAgkApTEhAkyXv/WKs6h6JO1alhn7N21e/zPPs55+zxreFd7x7W3kcRgZmZmdXLep0OwMzMzIbPBdzMzKyGXMDNzMxqyAXczMyshlzAzczMasgF3MzMrIZcwK1tJJ0s6RstzHe+pDPaEZOZNSfpq5L+uYX5rpX0oXbEZOu4gBdO0nGSeiQ9Len8Aaa/RdLvJT0p6f9J2rFh2oaSzpW0QtKDkk5sY9x7SVrSOC4i/i0inOQ2ruVi9pSkVXm4o9/0IyTdI+kJSVdI2rph2taSLs/T7pF0RBvjfp+k6xvHRcSHI+KT7YrBhscFvHwPAGcA5/afIGkyMB/4Z2BroAe4pGGW04DpwI7Am4CPS9q34niRtH7V2zAr3HERsVkeXtI3UtJuwNeAo4FtgCeBLzcs9yXgL3nakcBX8jKVcs7Wkwt44SJifkRcATw8wOR3A7dFxHcj4ilSwd5D0kvz9L8BPhkRj0bE7cDXgfcNtJ28932DpC9Kejwf1b+lYfr7Jd0uaaWkuyUd2zBtL0lLJJ0k6UFgHvB/gakNRyFTJZ0m6TsNy71O0o2SHpN0n6Rmse0vaWGe70ZJr2iYdpKk+3NcdzTGbFagI4GrIuK6iFhF2vl+t6TNJW0KvAf454hYFRHXA1eSiv1z5Hz6nqRL8v//zZL2aJg+R9JdedrvJB3UMK0v378g6RHSjv9XgdfmfH0sz/esy1mSDsi5uCKve8ADAkkfyO3Fo5J+3HdmUMkXJC3L7cwtknYf3a904nIBr7fdgN/2fYiIJ4C7gN0kbQVMbZye3w+2N/8a4G5gMnAqML/h9N4yYH9gC+D9wBckvbJh2ReSzgLsSNpx2A94oOEo5IHGDUnagVTkvwh0ATOAhf0Dyts4FzgWeAHp6OXKfHngJcBxwKsjYnPgbcDiQX4+s3b5d0nLc5Hcq2F8/5y9i3TE/eI8rImIPzTMP1TOHgB8l5R7FwFXSNogT7sLeD3wfOBfge9I2rZh2b58nwIcBXwYuCnn65b9NyRpFvAt4GPAlsAbGCDfJB0InEw6wOgC/ou0Uw+wT17uxXkd72XggxNrgQt4vW0GPN5v3OPA5nka/ab3TWtmGfCfEfFMRFwC3AG8AyAifhgRd0XyC+AnpMahz1rg1Ih4OiL+3ELsRwI/i4h5eXsPR8RzCjjwt8DXIuJXEbEmIi4Angb+GlgDbAjsKmmDiFicG0SzTjoJ2AnYDpgLXCVp5zxtqJxtNq2ZBRHxvYh4Bvg8sBEpN8hn5h6IiLU5n/8IzGpY9oGI+GJErG4xZz8InBsRP83rvD8ifj/AfMcC/x4Rt0fEauDfgBn5KPyZ/PO8FFCeZ2kL27YBuIDX2yrSEXGjLYCVeRr9pvdNa+b+ePa329xDOopH0n6SfinpkXx67e2kI/U+vfk0fqumkY4QhrIj8NF8+vyxvO1pwNSIuBM4gXTpYJmkiyVNHUYMZmMu72yuzDuzFwA3kPIFhs7ZZtOaua9hu2uBJazL2b9puPT0GLA7z87Z+xie4eTsWQ3bfQQQsF1E/Bw4h3St/yFJcyX1/5mtRS7g9XYb0HjNa1NgZ9J18UeBpY3T8/vbBlnfdpLU8HkH4AFJGwKXAf8BbJNPr/2IlJR9+n+t3VBfc3dfjnUo9wGfiogtG4ZNImIeQERcFBGvIzUaAXy6hXWatVOwLlf65+xOpLNIf8jD+pKmNyw7VM5Oa1jXesD2pJzdkdTn5TjgBTlnF9G+nD22X85uHBE3AkTE2RHxKtKlgReTTsnbCLiAF07S+pI2AiYBkyRt1NBj9HJgd0nvyfP8C3BLw2mtbwGnSNoqd2z7W+D8QTY3BThe0gaSDgFeRirUzyM1Mr3Aakn7ka5lDeYh4AWSnt9k+oXA3pIOzT/jCyTNGGC+rwMflvSa3AFmU0nvyJ1+XiLpzXkH4yngz6TT6mYdIWlLSW/ry1NJR5Ku+f44z3Ih8E5Jr8873KcD8/MR+xOku0pOz//ne5KucX97kE2+StK7c5twAuny0i+BTUkFuTfH9X7SEfhgHgK2l/S8JtO/Cbxf6dbV9SRtp3UdZht9FfiEcu95Sc/P7QmSXp1zeQPgCVLeOmdHyAW8fKeQCtMcUkeTP+dxREQvqdfqp4BHSZ1SDmtY9lTSKa97gF8An42IqwfZ1q9It50tz+s8OF+bXgkcD1yat3MEqXdsU3knYh5wdz6VNrXf9HtJpxU/SjrFtpBnny3om6+HtONxTt72nazrSb8hcGaO90HSDsjJg8VlVrENSLd99pL+Lz8CHBgRdwBExG2kzmIXkvqcbA78fcPyfw9snKfNA/4uL9PM90kdwR4l9VZ/d+5T8jvgc8BNpML8ctKp/MH8nHS0/6Ck5f0nRsSvyR1YSdfmf0E689V/vstJZ8IulrSCdOS/X568BWmn/FFSu/Qw6cyejYCefcnTJiqlW7g+lE9Hm1nhJJ0G7BIRR3U6FusMH4GbmZnVkAu4mZlZDfkUupmZWQ35CNzMzKyGKnuAvaRppNuYXkh6StfciDgrP5rzEqCb9Bi+Q/M9y01Nnjw5uru7qwrVbFxYsGDB8ojo6nQc4Jw1a8Voc7bKb6BZDXw0Im6WtDmwQNJPSbcAXRMRZ0qaQ7o96qTBVtTd3U1PT0+FoZrVn6R7Oh1DH+es2dBGm7OVnUKPiKURcXN+vxK4nfRs4AOAC/JsFwAHVhWDmZnZeNWWa+CSuoG/Ij0oZJu+h9fn1yntiMHMzGw8qbyAS9qM9BztEyJixTCWmy2pR1JPb29vdQGa2Zhwzpq1V6UFPD/v9jLgwoiYn0c/1PedtPl12UDLRsTciJgZETO7uorol2Nmg3DOmrVXlb3QRXr4/e0R8fmGSVcCx5CeYX0M6Vm+Nkrdc344quUXn/mOMYpkYKOND54bY+k/s1mdjXV+VZGvE70NqLIX+p6kh+vfKmlhHncyqXBfKumDwL3AIRXGYNayid4YWHv5/81Gq7ICHhHX8+zvnm30lqq2a2PDjYtZvThnJx4/ic3MzKyGXMDNzMxqyAXczMyshlzAzczMasgF3MzMrIZcwM3MzGrIBdzMzKyGXMDNzMxqyAXczMyshlzAzczMasgF3MzMrIZcwM3MzGrIBdzMzKyGXMDNzMxqyAXczMyshlzAzczMasgF3MzMrIZcwM3MzGrIBdzMzKyGXMDNzMxqqLICLulcScskLWoYd5qk+yUtzMPbq9q+mZnZeFblEfj5wL4DjP9CRMzIw48q3L6Zmdm4VVkBj4jrgEeqWr+ZmdlE1olr4MdJuiWfYt+q2UySZkvqkdTT29vbzvjMbAScs2bt1e4C/hVgZ2AGsBT4XLMZI2JuRMyMiJldXV3tis/MRsg5a9ZebS3gEfFQRKyJiLXA14FZ7dy+mZnZeNHWAi5p24aPBwGLms1rZmZmza1f1YolzQP2AiZLWgKcCuwlaQYQwGLg2Kq2b2ZmNp61VMAlbQzsEBF3tLriiDh8gNHfbHV5MyvHSNoAM6vWkKfQJb0TWAhcnT/PkHRl1YGZWRncBpiVqZVr4KeROps9BhARC4Hu6kIys8KchtsAs+K0UsBXR8TjlUdiZqVyG2BWoFaugS+SdAQwSdJ04HjgxmrDMrOCuA0wK1ArR+AfAXYDngbmASuAE6oMysyK4jbArEBDHoFHxJPAP+XBzCYYtwFmZRqygEuaCZxM6rTyP/NHxCuqC8vMSuE2wKxMrVwDvxD4GHArsLbacMysQG4DzArUSgHvjQjf82k2cbkNMCtQKwX8VEnfAK4hdWIBICLmVxaVmZXEbYBZgVop4O8HXgpswLrTZwE4ec0mBrcBZgVqpYDvEREvrzwSMyuV2wCzArVyH/gvJe1aeSRmViq3AWYFauUI/HXAMZL+RLr+JSB8C4nZhOE2wKxArRTwfSuPwsxK5jbArEBNC7ikLSJiBbCyjfGYWSHcBpiVbbAj8IuA/YEFpB6napgWwE4VxmVmnec2wKxgTQt4ROyfX1/UvnDMrBRuA8zKNmQvdEnXtDLOzMYntwFmZRrsGvhGwCbAZElbse702RbA1DbEZmYd5DbArGyDXQM/lvSdv1NJ18D6kncF8KWhVizpXNL1s2URsXsetzVwCelbjRYDh0bEoyOM3cyqNao2wMyq1fQUekScla99/WNE7BQRL8rDHhFxTgvrPp/n3n4yB7gmIqaTnqs8Z6SBm1m1xqANMLMKDXkNPCK+OJIVR8R1wCP9Rh8AXJDfXwAcOJJ1m1n7jLQNMLNqtfIo1bG0TUQsBcivU5rNKGm2pB5JPb29vW0L0MxGxjlr1l5NC7ikPfPrhu0LZ52ImBsRMyNiZldXVydCMJvQhtsGOGfN2muwI/Cz8+tNY7i9hyRtC5Bfl43hus1sbFXRBpjZGBmsF/ozks4DtpN0dv+JEXH8CLZ3JXAMcGZ+/f4I1mFm7VFFG2BmY2SwAr4/sDfwZtItJMMiaR6wF+ke0iXAqaTCfamkDwL3AocMd71m1jajagPMrFqDPUp1OXCxpNsj4rfDXXFEHN5k0luGuy4za7/RtgFmVq1WeqE/LOlyScskPSTpMknbVx6ZmZXCbYBZgVop4OeRrl1PBbYDrsrjzGxicBtgVqBWCviUiDgvIlbn4XzA94iYTRxuA8wK1EoB75V0lKRJeTgKeLjqwMysGG4DzArUSgH/AHAo8CCwFDg4jzOzicFtgFmBBruNDICIuBd4VxtiMbMCuQ0wK1O7n4VuZmZmY8AF3MzMrIZcwM3MzGqo5QIu6a8l/VzSDZL8Pd5mE4zbALOyNO3EJumFEfFgw6gTSR1ZBNwIXFFxbGbWQW4DzMo2WC/0r0paAHw2Ip4CHgOOANYCK9oRnJl1lNsAs4I1PYUeEQcCC4EfSDoaOIGUuJsAPn1mNs65DTAr26DXwCPiKuBtwJbAfOCOiDg7InrbEZyZdZbbALNyNS3gkt4l6Xrg58Ai4DDgIEnzJO3crgDNrDPcBpiVbbBr4GcArwU2Bn4UEbOAEyVNBz5FSmYzG7/cBpgVbLAC/jgpQTcGlvWNjIg/4sQ1mwjcBpgVbLBr4AeROqusJvU8NbOJxW2AWcGaHoFHxHLgi22MxcwK4jbArGxDfhtZFSQtBlYCa4DVETGzE3GYmZnVVUcKePamvIdvZmZmw+QvMzEzM6uhThXwAH4iaYGk2R2KwczMrLY6dQp9z4h4QNIU4KeSfh8R1zXOkAv7bIAddtihEzGa2TA4Z83aqyNH4BHxQH5dBlwOzBpgnrkRMTMiZnZ1dbU7RDMbJuesWXu1vYBL2lTS5n3vgX1Ij2k0MzOzFnXiFPo2wOWS+rZ/UURc3YE4zMzMaqvtBTwi7gb2aPd2zczMxhPfRmZmZlZDLuBmZmY15AJuZmZWQy7gZmZmNeQCbmZmVkMu4GZmZjXkAm5mZlZDLuBmZmY15AJuZmZWQy7gZmZmNeQCbmZmVkMu4GZmZjXkAm5mZlZDLuBmZmY15AJuZmZWQ23/PvA66p7zw1GvY/GZ7xiDSMysFaPNWeer1YGPwM3MzGrIBdzMzKyGXMDNzMxqyAXczMyshjpSwCXtK+kOSXdKmtOJGMzMzOqs7QVc0iTgS8B+wK7A4ZJ2bXccZmZmddaJ28hmAXdGxN0Aki4GDgB+N1YbqMMtJHWI0UbHf+PWlf67Kj0+Gxt1+zsrItq7QelgYN+I+FD+fDTwmog4rt98s4HZ+ePuwKK2Bjo8k4HlnQ5iCKXHWHp8UH6ML4mIzTu1cefsmCs9xtLjg/JjHFXOduIIXAOMe85eRETMBeYCSOqJiJlVBzZSpccH5cdYenxQfoySejq5fefs2Co9xtLjg/JjHG3OdqIT2xJgWsPn7YEHOhCHmZlZbXWigP8GmC7pRZKeBxwGXNmBOMzMzGqr7afQI2K1pOOAHwOTgHMj4rYhFptbfWSjUnp8UH6MpccH5cdYUnwlxTKQ0uOD8mMsPT4oP8ZRxdf2Tmw29iStiojNJE0Fzo6IgzsdUzOS3gfM7N9p0Wwicc7aWPCT2MaRiHig5IZgLOTnCJiNC85ZGw0X8HFEUrekRfn9+yTNl3S1pD9K+kzDfPtIuknSzZK+K2mzAdZ1raRPS/q1pD9Ien3Des9pmO8HkvbK71flZRZI+pmkWXk9d0t6V8Pqp+W47pB0asO6jsrbWyjpa32Jn9d7uqRfAa8d41+bWcc4Z200XMDHtxnAe4GXA++VNE3SZOAUYO+IeCXQA5zYZPn1I2IWcAJwapN5Gm0KXBsRrwJWAmcAbwUOAk5vmG8WcGSO7xBJMyW9LMe6Z0TMANbkefrWuygiXhMR17f4s5vVkXPWWtaJ+8Ctfa6JiMcBJP0O2BHYkvQI2xskATwPuKnJ8vPz6wKgu4Xt/QW4Or+/FXg6Ip6RdGu/5X8aEQ/nuOYDrwNWA68CfpPj2hhYludfA1zWwvbN6s45ay1zAR/fnm54v4b09xYpGQ8fxvJ9y0JK2sYzNxs1vH8m1vWKXNu3fESsldT4v9a/52TkuC6IiE8MEMdTEbGmhXjN6s45ay3zKfSJ55fAnpJ2AZC0iaQXD2P5xcAMSetJmkY6tTZcb5W0taSNgQOBG4BrgIMlTclxbS1pxxGs22y8cc7agHwEPsFERG++LWSepA3z6FOAP7S4ihuAP5FOty0Cbh5BGNcD3wZ2AS6KiB4ASacAP5G0HvAM8A/APSNYv9m44Zy1ZnwfuJmZWQ35FLqZmVkNuYCbmZnVkAu4mZlZDbmAm5mZ1ZALuJmZWQ25gJuZmdWQC7iZmVkNuYCbmZnVkAu4mZlZDbmAm5mZ1ZALuJmZWQ25gJuZmdWQC7h1jKRVknYaYp5uSdHvu4nNrAOcs2VxAS+MpOMk9Uh6WtL5A0x/i6TfS3pS0v9r/P5dSRtKOlfSCkkPSjqx1WWrJulaSR9qHBcRm0XE3e2KwWys5Zz7pqR7JK2U9N+S9us3j3PWKuECXp4HgDOAc/tPkDQZmA/8M7A10ANc0jDLacB0YEfgTcDHJe3b4rKVUOL/Mxuv1gfuA94IPJ+UX5dK6gbnrFUsIjwUOJCK+Pn9xs0Gbmz4vCnwZ+Cl+fP9wD4N0z8JXNzKsgNsfzHwCeB3wKPAecBGedpWwA+A3jztB8D2DcteC3wKuCFv40JgDfAUsAo4J88XwC75/cbA54B7gMeB6/O47jzf+nm+5wPfBJbmn/cMYFKetgvwi7z8cuCSTv8dPUy8AbgFeE9+75x1zlY2eC+rXnYDftv3ISKeAO4CdpO0FTC1cXp+v9tQyw6yvSOBtwE7Ay8GTsnj1yM1DjsCO5AS/px+yx5NaoA2B94H/BdwXKRTcMcNsK3/AF4F/C/S0cbHgbUDzHcBsJqU+H8F7AP0neb7JPATUmO1PfDFQX42szEnaRtSrtyWRzlnnbOVcQGvl81Ie6qNHicl3GYNn/tPG2rZZs6JiPsi4hHS3vnhABHxcERcFhFPRsTKPO2N/ZY9PyJui4jVEfHMYD9UPl33AeB/R8T9EbEmIm6MiKf7zbcNsB9wQkQ8ERHLgC8Ah+VZniE1UFMj4qmIuH6w7ZqNJUkbkI5cL4iI3+fRzlnnbGVcwOtlFbBFv3FbACvzNPpN75s21LLN3Nfw/h7S0QKSNpH0tdxxZwVwHbClpElNlh3KZGAj0tHFYHYENgCWSnpM0mPA14ApefrHAQG/lnSbpA8MIwazEcsF7dvAX4DGo1XnrHO2Mi7g9XIbsEffB0mbkk6V3RYRj5KuMe3RMP8erDuV13TZQbY3reH9DqQOdgAfBV4CvCYitgDe0Lfahvmj37r6f260nHStbedB5oHUwDwNTI6ILfOwRUTsBhARD0bE30bEVOBY4MuSdhlinWajIkmka7zbkK59Nx69Omeds5VxAS+MpPUlbQRMAiZJ2qjhfsrLgd0lvSfP8y/ALQ2n674FnCJpK0kvBf4WOL/FZQfyD5K2l7Q1cDLresBuTrqG9liedmoLP9pDwID3j0bEWlKv+89LmippkqTXStqw33xLSdfLPidpC0nrSdpZ0hsBJB0iafs8+6OkBmhNC7GZjcZXgJcB74yIP/eb5px1zlan073oPDx7IN1WEv2G0xqm7w38npSM1wLdDdM2JCXVClLyndhv3U2XHSCOxazr0foYqSPKJnna1Lz8KuAPpD3nxl6n1wIf6re+1+Z5HwXOzuP692j9T1Iv1cdJp/ia9Wj9CrAkz/ffwGF52mfy8qtIp/Zmd/rv6WF8D6RTxMG63tp9w5EN8zhnnbOVDMq/RLNnkbSYlNA/63QsZjY05+zE41PoZmZmNeQCbmZmVkM+hW5mZlZDPgI3MzOrocq+7k3SNNItEi8kPV5vbkSclW9huITUU3ExcGik+yGbmjx5cnR3d1cVqtm4sGDBguUR0dXpOMA5a9aK0eZsld/Xuhr4aETcLGlzYIGkn5KesXtNRJwpaQ4wBzhpsBV1d3fT09NTYahm9Sfpnk7H0Mc5aza00eZsZafQI2JpRNyc368Ebge2Aw4g3Z9Ifj2wqhjMzMzGq7ZcA8/fjftXwK+AbSI9nYf8OqXJMrMl9Ujq6e3tbUeYZjYKzlmz9qq8gEvaDLiM9G00K1pdLiLmRsTMiJjZ1VXEZT0zG4Rz1qy9Ki3g+ev1LgMujIj5efRDkrbN07cFllUZg5mZ2XhUZS/0vm/ouT0iPt8w6UrgGODM/Pr9qmKYSLrn/HBUyy8+8x1jFMnARhsfPDfGsf6ZS/8d2vhS+v9bHfKr9N9h1arshb4ncDRwq6SFedzJpMJ9qaQPAvcCh1QYg43QRE8Ms7pxzk48lRXwiLieZ3/XbKO3VLVdMzOzicBPYjMzM6shF3AzM7MacgE3MzOrIRdwMzOzGnIBNzMzqyEXcDMzsxpyATczM6shF3AzM7MacgE3MzOrIRdwMzOzGnIBNzMzqyEXcDMzsxpyATczM6shF3AzM7MacgE3MzOrIRdwMzOzGnIBNzMzqyEXcDMzsxpyATczM6shF3AzM7MaqqyASzpX0jJJixrGnSbpfkkL8/D2qrZvZmY2nlV5BH4+sO8A478QETPy8KMKt29mZjZuVVbAI+I64JGq1m9mZjaRdeIa+HGSbsmn2LdqNpOk2ZJ6JPX09va2Mz4zGwHnrFl7tbuAfwXYGZgBLAU+12zGiJgbETMjYmZXV1e74jOzEXLOmrVXWwt4RDwUEWsiYi3wdWBWO7dvZmY2XrS1gEvatuHjQcCiZvOamZlZc+tXtWJJ84C9gMmSlgCnAntJmgEEsBg4tqrtm5mZjWctFXBJGwM7RMQdra44Ig4fYPQ3W13ezMoxkjbAzKo15Cl0Se8EFgJX588zJF1ZdWBmVga3AWZlauUa+GmkzmaPAUTEQqC7upDMrDCn4TbArDitFPDVEfF45ZGYWancBpgVqJVr4IskHQFMkjQdOB64sdqwzKwgbgPMCtTKEfhHgN2Ap4F5wArghCqDMrOiuA0wK9CQR+AR8STwT3kwswnGbYBZmYYs4JJmAieTOq38z/wR8YrqwjKzUrgNMCtTK9fALwQ+BtwKrK02HDMrkNsAswK1UsB7I8L3fJpNXG4DzArUSgE/VdI3gGtInVgAiIj5lUVlZiVxG2BWoFYK+PuBlwIbsO70WQBOXrOJwW2AWYFaKeB7RMTLK4/EzErlNsCsQK3cB/5LSbtWHomZlcptgFmBWjkCfx1wjKQ/ka5/CQjfQmI2YbgNMCtQKwV838qjMLOSuQ0wK1DTAi5pi4hYAaxsYzxmVgi3AWZlG+wI/CJgf2ABqcepGqYFsFOFcZlZ57kNMCtY0wIeEfvn1xe1LxwzK4XbALOyDdkLXdI1rYwzs/HJbYBZmQa7Br4RsAkwWdJWrDt9tgUwdagVSzqXdPptWUTsnsdtDVxC+lKExcChEfHoKOI3s4qMtg0ws2oNdgR+LOna10vza9/wfeBLLaz7fJ7be3UOcE1ETCc9lnHOMOM1s/YZbRtgZhUa7Br4WcBZkj4SEV8c7ooj4jpJ3f1GHwDsld9fAFwLnDTcdZtZ9UbbBphZtYa8Bj7GibtNRCzN610KTBnDdZtZBVy8zcrUyqNUO0LSbEk9knp6e3s7HY6ZDcE5a9ZeTQu4pD3z64ZjuL2HJG2b17stsKzZjBExNyJmRsTMrq6uMQzBzFox3DbAOWvWXoMdgZ+dX28aw+1dCRyT3x9D6gxjZmWqog0wszEy2JPYnpF0HrCdpLP7T4yI4wdbsaR5pA5rkyUtAU4FzgQulfRB4F7gkJEGbmaVG1UbYGbVGqyA7w/sDbyZdOvIsETE4U0mvWW46zKzjhhVG2Bm1RrsNrLlwMWSbo+I37YxJjMrgNsAs7K10gv9YUmXS1om6SFJl0navvLIzKwUbgPMCtRKAT+P1PlsKrAdcFUeZ2YTg9sAswK1UsCnRMR5EbE6D+cDvkfEbOJwG2BWoFYKeK+koyRNysNRwMNVB2ZmxXAbYFagVgr4B4BDgQeBpcDBeZyZTQxuA8wKNNhtZABExL3Au9oQi5kVyG2AWZmKfRa6mZmZNecCbmZmVkMu4GZmZjXUcgGX9NeSfi7pBkkHVhmUmZXHbYBZWZp2YpP0woh4sGHUiaSOLAJuBK6oODYz6yC3AWZlG6wX+lclLQA+GxFPAY8BRwBrgRXtCM7MOsptgFnBmp5Cj4gDgYXADyQdDZxAStxNAJ8+Mxvn3AaYlW3Qa+ARcRXwNmBLYD5wR0ScHRG97QjOzDrLbYBZuZoWcEnvknQ98HNgEXAYcJCkeZJ2bleAZtYZbgPMyjbYNfAzgNcCGwM/iohZwImSpgOfIiWzmY1fbgPMCjZYAX+clKAbA8v6RkbEH3Himk0EbgO8p+ZwAAAdJElEQVTMCjbYNfCDSJ1VVpN6nprZxOI2wKxgTY/AI2I58MU2xmJmBXEbYFa2Ib+NrAqSFgMrgTXA6oiY2Yk4zMzM6qojBTx7U97DNzMzs2Hyl5mYmZnVUKcKeAA/kbRA0uyBZpA0W1KPpJ7eXj8zwqx0zlmz9upUAd8zIl4J7Af8g6Q39J8hIuZGxMyImNnV1dX+CM1sWJyzZu3VkQIeEQ/k12XA5cCsTsRhZmZWV20v4JI2lbR533tgH9JjGs3MzKxFneiFvg1wuaS+7V8UEVd3IA4zM7PaansBj4i7gT3avV0zM7PxxLeRmZmZ1ZALuJmZWQ25gJuZmdWQC7iZmVkNuYCbmZnVkAu4mZlZDbmAm5mZ1ZALuJmZWQ25gJuZmdWQC7iZmVkNuYCbmZnVkAu4mZlZDbmAm5mZ1ZALuJmZWQ25gJuZmdVQ278PvI665/xw1OtYfOY7xiASM2vFaHPW+Wp14CNwMzOzGnIBNzMzqyEXcDMzsxpyATczM6uhjhRwSftKukPSnZLmdCIGMzOzOmt7AZc0CfgSsB+wK3C4pF3bHYeZmVmddeI2slnAnRFxN4Cki4EDgN+N1QbqcAtJHWK00fHfuHWl/65Kj8/GRt3+zoqI9m5QOhjYNyI+lD8fDbwmIo7rN99sYHb+uDuwqK2BDs9kYHmngxhC6TGWHh+UH+NLImLzTm3cOTvmSo+x9Pig/BhHlbOdOALXAOOesxcREXOBuQCSeiJiZtWBjVTp8UH5MZYeH5Qfo6SeTm7fOTu2So+x9Pig/BhHm7Od6MS2BJjW8Hl74IEOxGFmZlZbnSjgvwGmS3qRpOcBhwFXdiAOMzOz2mr7KfSIWC3pOODHwCTg3Ii4bYjF5lYf2aiUHh+UH2Pp8UH5MZYUX0mxDKT0+KD8GEuPD8qPcVTxtb0Tm409SasiYjNJU4GzI+LgTsfUjKT3ATP7d1o0m0icszYW/CS2cSQiHii5IRgL+TkCZuOCc9ZGwwV8HJHULWlRfv8+SfMlXS3pj5I+0zDfPpJuknSzpO9K2myAdV0r6dOSfi3pD5Je37Decxrm+4GkvfL7VXmZBZJ+JmlWXs/dkt7VsPppOa47JJ3asK6j8vYWSvpaX+Ln9Z4u6VfAa8f412bWMc5ZGw0X8PFtBvBe4OXAeyVNkzQZOAXYOyJeCfQAJzZZfv2ImAWcAJzaZJ5GmwLXRsSrgJXAGcBbgYOA0xvmmwUcmeM7RNJMSS/Lse4ZETOANXmevvUuiojXRMT1Lf7sZnXknLWWdeI+cGufayLicQBJvwN2BLYkPcL2BkkAzwNuarL8/Py6AOhuYXt/Aa7O728Fno6IZyTd2m/5n0bEwzmu+cDrgNXAq4Df5Lg2Bpbl+dcAl7WwfbO6c85ay1zAx7enG96vIf29RUrGw4exfN+ykJK28czNRg3vn4l1vSLX9i0fEWslNf6v9e85GTmuCyLiEwPE8VRErGkhXrO6c85ay3wKfeL5JbCnpF0AJG0i6cXDWH4xMEPSepKmkU6tDddbJW0taWPgQOAG4BrgYElTclxbS9pxBOs2G2+cszYgH4FPMBHRm28LmSdpwzz6FOAPLa7iBuBPpNNti4CbRxDG9cC3gV2AiyKiB0DSKcBPJK0HPAP8A3DPCNZvNm44Z60Z3wduZmZWQz6FbmZmVkMu4GZmZjXkAm5mZlZDLuBmZmY15AJuZmZWQy7gZmZmNeQCbmZmVkMu4GZmZjXkAm5mZlZDLuBmZmY15AJuZmZWQy7gZmZmNeQCbh0haQdJqyRNGmK+vSQtaVdcZtacpNsk7dXCfNH39adWHRfwQkmaLukpSd/pN/4ISfdIekLSFZK2bpi2taTL87R7JB3R6rJt+HkWS9q773NE3BsRm0XEmnbFYFYFSd2SfiTpUUkPSjpH0voN02dIWiDpyfw6o2GaJH1a0sN5+IwktbJsG36u8yWd0TguInaLiGvbFYMNzgW8XF8CftM4QtJuwNeAo4FtgCeBL/db5i952pHAV/IyrSxbicaGzGyc+jKwDNgWmAG8Efh7AEnPA74PfAfYCrgA+H4eDzAbOBDYA3gFsD9wbIvLVmaoM2NWiIjwUNgAHAZcCpwGfKdh/L8BFzV83plUsDcHNs3vX9ww/dvAmUMt2ySGAI4H7gaWA58F1mtY9ufAw3nahcCWDcsuBk4CbgGeBuYBa4E/A6uAjwPdeRvr52W2Bs4DHgAeBa7I4/cCljSseypwGdAL/Ak4vmHaLKAHWAE8BHy+039LD+N/AG4H3t7w+bPA1/L7fYD7ATVMvxfYN7+/EZjdMO2DwC9bWXaAOK4F/h34NfA4qfhv3TD9u8CDedp1wG4N084HvgL8CHiCtGPxTG4jVgFX5fkWA3vn95OAk4G7gJXAAmBanhbALvn9hsB/5NgfAr4KbJynTQZ+ADwGPAL8V18742HowUfghZG0BXA68NEBJu8G/LbvQ0TcRS7aeVgTEX9omP+3eZmhlm3mIGAm8ErgAOADfWGSGoqpwMuAaaSdjUaHA+8gFfbDScn7zkinzT8zwLa+DWyS45wCfKH/DJLWA67KP8d2wFuAEyS9Lc9yFnBWRGxB2sm4dJCfzWysnAUcJmkTSdsB+wFX52m7AbdErlbZLTTJS56bs4MtO5C/IeXpVGA1cHbDtP8LTCfl182kHe9GRwCfIh0QfCtP/0zO2XcOsK0TSXn+dmCLvN0nB5jv06R2ZgawCyl3/yVP+yiwBOginRk8mVT8rQUu4OX5JPDNiLhvgGmbkfaeGz1OSrjBpg21bDOfjohHIuJe4D9JyUpE3BkRP42IpyOiF/g86bRho7Mj4r6I+PMg6wdA0rakRu/DEfFoRDwTEb8YYNZXA10RcXpE/CUi7ga+TjpjAemIYRdJkyNiVUT8cqhtm42BX5CK6gpSMeoBrsjThpuXjwOb5evgI8nZb0fEooh4Avhn4NC+0+ERcW5ErIyIp0k73HtIen7Dst+PiBsiYm1EPDXUDw18CDglIu6I5LcR8XDjDPnn+Fvg/+S2ZCXpbGBjzm4L7Jjz/r/67bDYIFzAC5I7qOzNAEef2SrSnm6jLUinrwabNtSyzTTuRNxD2qtH0hRJF0u6X9IK0jW6yYMsO5RpwCMR8egQ8+0ITJX0WN9A2mPfJk//IGlP//eSfiNp/2HEYDZs+azQj4H5pMtYk0nXqz+dZxluXm4BrMpFbCxydgNgsqRJks6UdFfO2cV5nslNlm3FNNLp88F0kc6sLWjI2avzeEiXG+4EfiLpbklzhhnDhOYCXpa9SNeG75X0IPCPwHsk3Zyn30bq7AKApJ1I15f+kIf1JU1vWN8eeZmhlm1mWsP7HUjXpyGdPg/gFfl09VGk0+qN+u9FD7ZXfR+wtaQtB5mnb74/RcSWDcPmEfF2gIj4Yz5dP4XUgH5P0qZDrNNsNLYm5ck5+YzUw6S+HG/P028DXtHYs5zUWW3AvOS5OTvYsgPpn7PPkPqpHEG6DLY38HxSOwPPztvh5CykfNx5iHmWk/q+7NaQs8+PiM0A8hmBj0bETsA7gRMlvWWIdVrmAl6WuaSEmJGHrwI/BPqu8V4IvFPS63NhOh2Yn5PgCdJRwOmSNpW0Jylhvz3UsoPE8zFJW0maBvxv4JI8fnPS0cFj+Zrfx1r42R4CdhpoQkQsJV2f+3Le3gaS3jDArL8GVkg6SdLG+ahid0mvBpB0lKSuiFhL6hQD4NvUrDIRsZzUmfLvJK2fd0KPYd117WtJ/4PHS9pQ0nF5/M/z67dIRWs7SVNJ14TPb3HZgRwlaVdJm5By/HuRbtXcnNSh9GHSEfG/tfDjNc3Z7BvAJ5VueZWkV0h6QeMMORe/DnxB0hSA/LO+Lb/fX9IueSdlRf55nbMtcgEvSEQ8GREP9g2kIvlUvs5MRNwGfJhUjJeRkvLvG1bx98DGedo84O/yMq0sO5Dvk3qWLiTtSHwzj/9XUse2x/P4+S38eP8OnJJPo/3jANOPJh0t/D7Hd0L/GXJD9E7Szs2fSHv33yAdUQDsC9wmaRW5Y1GL1/LMRuPdpP+9XtLp4NXA/wGIiL+QbhP7G9JO5QeAA/N4SLd2XgXcCiwi5dPXWlx2IN8m7QA8CGxEupME0o7CPaRe7b8DWukf8k1g15yzVwww/fOkjqI/IRXfb5Lan/5OIv1efplP3/8MeEmeNj1/XgXcBHw5fJ95y+T+AjYQSQFMj4g7Ox2LmQ1N0rWk206/0elYrD18BG5mZlZDLuBmZmY15FPoZmZmNeQjcDMzsxqq7Ism8q1H3wJeSHoO9tyIOEvpG7AuId2HuBg4dKgHeEyePDm6u7urCtVsXFiwYMHyiOgaes7qOWfNhjbanK3ym6JWAx+NiJslbU56Es9PgfcB10TEmfmpO3NItxk01d3dTU9PT4WhmtWfpHs6HUMf56zZ0Eabs5WdQo+IpRFxc36/kvSNPduRHi5yQZ7tAtJ9jmZmZjYMbbkGLqkb+CvgV8A2+clbfU/gmtJkmdmSeiT19Pb2tiNMMxsF56xZe1VewCVtRvr+5hMiYkWry0XE3IiYGREzu7qKuKxnZoNwzpq1V5XXwJG0Aal4XxgRfY/bfEjSthGxNH+N5LIqY7AydM/54ajXsfjMd4zpOvuvz8zWGev8qiJfJ3obUGUvdJGejXt7RHy+YdKVpIf9n5lfv19VDBPJRP9HHgv+HVo7+f/NRqvKI/A9SV9QcaukhXncyaTCfamkDwL3AodUGIONkBsXs3pxzk48lRXwiLie535HdB9/36uZmdko+ElsZmZmNeQCbmZmVkMu4GZmZjXkAm5mZlZDLuBmZmY15AJuZmZWQy7gZmZmNeQCbmZmVkMu4GZmZjXkAm5mZlZDLuBmZmY15AJuZmZWQy7gZmZmNeQCbmZmVkMu4GZmZjXkAm5mZlZDLuBmZmY15AJuZmZWQy7gZmZmNeQCbmZmVkOVFXBJ50paJmlRw7jTJN0vaWEe3l7V9s3MzMazKo/Azwf2HWD8FyJiRh5+VOH2zczMxq3KCnhEXAc8UtX6zczMJrJOXAM/TtIt+RT7Vs1mkjRbUo+knt7e3nbGZ2Yj4Jw1a692F/CvADsDM4ClwOeazRgRcyNiZkTM7Orqald8ZjZCzlmz9mprAY+IhyJiTUSsBb4OzGrn9s3MzMaLthZwSds2fDwIWNRsXjMzM2tu/apWLGkesBcwWdIS4FRgL0kzgAAWA8dWtX0zM7PxrKUCLmljYIeIuKPVFUfE4QOM/mary5tZOUbSBphZtYY8hS7pncBC4Or8eYakK6sOzMzK4DbArEytXAM/jdTZ7DGAiFgIdFcXkpkV5jTcBpgVp5UCvjoiHq88EjMrldsAswK1cg18kaQjgEmSpgPHAzdWG5aZFcRtgFmBWjkC/wiwG/A0MA9YAZxQZVBmVhS3AWYFGvIIPCKeBP4pD2Y2wbgNMCvTkAVc0kzgZFKnlf+ZPyJeUV1YZlYKtwFmZWrlGviFwMeAW4G11YZjZgVyG2BWoFYKeG9E+J5Ps4nLbYBZgVop4KdK+gZwDakTCwARMb+yqMysJG4DzArUSgF/P/BSYAPWnT4LwMlrNjG4DTArUCsFfI+IeHnlkZhZqdwGmBWolfvAfylp18ojMbNSuQ0wK1ArR+CvA46R9CfS9S8B4VtIzCYMtwFmBWqlgO9beRRmVjK3AWYFalrAJW0RESuAlW2Mx8wK4TbArGyDHYFfBOwPLCD1OFXDtAB2qjAuM+s8twFmBWtawCNi//z6ovaFY2alcBtgVrYhe6FLuqaVcWY2PrkNMCvTYNfANwI2ASZL2op1p8+2AKYOtWJJ55JOvy2LiN3zuK2BS0hfirAYODQiHh1F/GZWkdG2AWZWrcGOwI8lXft6aX7tG74PfKmFdZ/Pc3uvzgGuiYjppMcyzhlmvGbWPqNtA8ysQoNdAz8LOEvSRyLii8NdcURcJ6m73+gDgL3y+wuAa4GThrtuM6veaNsAM6vWkNfAxzhxt4mIpXm9S4EpzWaUNFtSj6Se3t7eMQzBzIaj1TbAOWvWXq08SrUjImJuRMyMiJldXV2dDsfMhuCcNWuvpgVc0p75dcMx3N5DkrbN690WWDaG6zazMVRRG2BmY2SwI/Cz8+tNY7i9K4Fj8vtjSJ1hzKxMVbQBZjZGBnsS2zOSzgO2k3R2/4kRcfxgK5Y0j9RhbbKkJcCpwJnApZI+CNwLHDLSwM2scqNqA8ysWoMV8P2BvYE3k24dGZaIOLzJpLcMd11m1hGjagPMrFqD3Ua2HLhY0u0R8ds2xmRmBXAbYFa2VnqhPyzpcknLJD0k6TJJ21cemZmVwm2AWYFaKeDnkTqfTQW2A67K48xsYnAbYFagVgr4lIg4LyJW5+F8wDd5mk0cbgPMCtRKAe+VdJSkSXk4Cni46sDMrBhuA8wK1EoB/wBwKPAgsBQ4OI8zs4nBbYBZgQa7jQyAiLgXeFcbYjGzArkNMCtTsc9CNzMzs+ZcwM3MzGrIBdzMzKyGWi7gkv5a0s8l3SDpwCqDMrPyuA0wK0vTTmySXhgRDzaMOpHUkUXAjcAVFcdmZh3kNsCsbIP1Qv+qpAXAZyPiKeAx4AhgLbCiHcGZWUe5DTArWNNT6BFxILAQ+IGko4ETSIm7CeDTZ2bjnNsAs7INeg08Iq4C3gZsCcwH7oiIsyOitx3BmVlnuQ0wK1fTAi7pXZKuB34OLAIOAw6SNE/Szu0K0Mw6w22AWdkGuwZ+BvBaYGPgRxExCzhR0nTgU6RkNrPxy22AWcEGK+CPkxJ0Y2BZ38iI+CNOXLOJwG2AWcEGuwZ+EKmzympSz1Mzm1jcBpgVrOkReEQsB77YxljMrCBuA8zKNuS3kVVB0mJgJbAGWB0RMzsRh5mZWV11pIBnb8p7+GZmZjZM/jITMzOzGupUAQ/gJ5IWSJo90AySZkvqkdTT2+tnRpiVzjlr1l6dKuB7RsQrgf2Af5D0hv4zRMTciJgZETO7urraH6GZDYtz1qy9OlLAI+KB/LoMuByY1Yk4zMzM6qrtBVzSppI273sP7EN6TKOZmZm1qBO90LcBLpfUt/2LIuLqDsRhZmZWW20v4BFxN7BHu7drZmY2nvg2MjMzsxpyATczM6shF3AzM7MacgE3MzOrIRdwMzOzGnIBNzMzqyEXcDMzsxpyATczM6shF3AzM7MacgE3MzOrIRdwMzOzGnIBNzMzqyEXcDMzsxpyATczM6shF3AzM7Maavv3gddR95wfjnodi898xxhEYmatGG3OOl+tDnwEbmZmVkMu4GZmZjXkAm5mZlZDLuBmZmY11JECLmlfSXdIulPSnE7EYGZmVmdt74UuaRLwJeCtwBLgN5KujIjfjdU26tADtQ4xmrVL6flQenw2MXXiNrJZwJ0RcTeApIuBA4AxK+BmJXCjb1YvdctZRUR7NygdDOwbER/Kn48GXhMRx/WbbzYwO3/cHVjU1kCHZzKwvNNBDKH0GEuPD8qP8SURsXmnNu6cHXOlx1h6fFB+jKPK2U4cgWuAcc/Zi4iIucBcAEk9ETGz6sBGqvT4oPwYS48Pyo9RUk8nt++cHVulx1h6fFB+jKPN2U50YlsCTGv4vD3wQAfiMDMzq61OFPDfANMlvUjS84DDgCs7EIeZmVlttf0UekSslnQc8GNgEnBuRNw2xGJzq49sVEqPD8qPsfT4oPwYS4qvpFgGUnp8UH6MpccH5cc4qvja3onNzMzMRs9PYjMzM6shF3AzM7MaKqKASzpX0jJJi/qN/0h+5Optkj7TMP4T+TGsd0h6W0nxSXqrpAWSbs2vb646vuHG2DBtB0mrJP1jafFJeoWkm/L4WyVtVFKMkjaQdEGO7XZJn+hEfJIukbQwD4slLWyYVkmelJ6vw42xEzlber6OJMZ252zp+dosxjHN2Yjo+AC8AXglsKhh3JuAnwEb5s9T8uuuwG+BDYEXAXcBkwqK76+Aqfn97sD9pf0OG6ZfBnwX+MeS4iN1rrwF2CN/fkHVf+MRxHgEcHF+vwmwGOhud3z9pn8O+Jf8vrI8KT1fRxBj23O29Hwdwe+w7Tlber42i7Hf9FHlbBFH4BFxHfBIv9F/B5wZEU/neZbl8QeQ/hBPR8SfgDtJj2ctIr6I+O+I6Luv/TZgI0kbVhnfcGMEkHQgcHeOsXLDjG8f4JaI+G0e/3BErCksxgA2lbQ+sDHwF2BFB+IDQJKAQ4F5eVRleVJ6vg43xk7kbOn5OoIY256zpefrIDECY5OzRRTwJl4MvF7SryT9QtKr8/jtgPsa5luSx7Vbs/gavQf4775/pg4YMEZJmwInAf/aobj6NPsdvhgIST+WdLOkjxcY4/eAJ4ClwL3Af0TEgInaJq8HHoqIP+bP7c6T0vMVys/Z0vMVys/ZuuQrjEHOduJRqq1aH9gK+Gvg1cClknaixUextsGA8UU+FyJpN+DTpD3TTmn2O/xX4AsRsSrtBBYX3/rA6/K4J4FrJC2IiGsKinEWsAaYmqf/l6SfRf6Sng44nHV78tD+PCk9X6H8nC09X6H8nK1LvsIY5GzJBXwJMD8n168lrSU9mL6UR7E2i69X0vbA5cDfRMRdHYhtqBhfAxycO3hsCayV9FREnFNIfEuAX0TEcgBJPyJdR+pEAW8W4xHA1RHxDLBM0g3ATNJpzrbKpwXfDbyqYXS786T0fIXyc7b0fB0sxlJytvh8hbHL2ZJPoV8BvBlA0ouB55G+VeZK4DBJG0p6ETAd+HUp8UnaEvgh8ImIuKEDcTUaMMaIeH1EdEdEN/CfwL91qDFo9jf+MfAKSZvkf/Q30rmvm20W473Am5VsStrj/32HYtwb+H1ELGkY1+48KT1fm8ZYUM6Wnq9NY6ScnK1DvsJY5WwrPemqHkinEZYCz5D2Qj5I+sV/h/SVhDcDb26Y/59IPfTuAPYrKT7gFNK1loUNw5SSYuy33Gm0pxf6cP/GR5E67CwCPlPa/yGwGalH8G2khupjnYgvjz8f+PAA81eSJ6Xn6wj+lm3P2dLzdYR/57bmbOn52izGPH5MctaPUjUzM6uhkk+hm5mZWRMu4GZmZjXkAm5mZlZDLuBmZmY15AJuZmZWQy7g44CkVfl1qqTvdTqewUh6n6RO3cNqVgTnrI0FF/BxJCIeiIiDOx1HlSRN6nQMZmPFOWuj4QI+jkjqVv7e2bzXPF/S1ZL+qGd/b+8+St/be7Ok70rabIB1XSvp05J+LekPkl7fsN5zGub7gaS98vtVeZkFkn4maVZez92S3tWw+mk5rjskndqwrqPy9hZK+lpf4uf1ni7pV8Brx/jXZtYxzlkbDRfw8W0G8F7g5cB7JU2TNJn05Km9I+KVQA9wYpPl14+IWcAJwKlN5mm0KXBtRLwKWAmcAbwVOAg4vWG+WcCROb5DJM2U9LIc654RMYP0xQNHNqx3UUS8JiKub/FnN6sj56y1rOQvM7HRuyYiHgeQ9DtgR9KXIewK3KD0zUbPA25qsvz8/LoA6G5he38Brs7vb+X/t3fHKnFFQRzGv7+NsbVIJzZinT51eos0dvZ5goAvIwGtzAOYNoKgXbo0ai3pI7qOxbkXF1G4a3fW71ctF+bcgcsw7BRn4Laq7pL8eRb/q6r+DXn9pG0xuqdd7H8+5LUGjLt8Z8DxhPdLvbNmNZkNfLnN7zSe0b53aMW4u0D8GAutaOcnNx/mft/V0928D2N8VT0MCw5Gz+/vrSGvg6r6/kIe/6tqNiFfqXfWrCZzhP7+nAGfk2wBDNuDtheIvwI+JVlJskEbrS3qS5L1JGvADnBKWzv4NcnHIa/1JJtvOFtaNtasXuQ/8Hemqm6S7AFHSVaHx/vA34lHnAKXtHHbuPFnUb+BH8AWcFhVFwBJ9oGTJCu07T3fgOs3nC8tDWtWr3EbmSRJHXKELklSh2zgkiR1yAYuSVKHbOCSJHXIBi5JUods4JIkdcgGLklShx4BLjiRjYOA+nMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f105f69af28>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(int(particle_numbers.shape[0] / 2), 2, figsize=(7, 7), sharex=True, sharey=True)\n",
    "for n, _ax in enumerate(ax.flat):\n",
    "    _ax.bar(lines[n, :], totals[n, :] / totals[n, :].sum() * 100);\n",
    "    _ax.axis([160, 170, None, None])\n",
    "    _ax.set_title(\"{} particles\".format(particle_numbers[n, 0]))\n",
    "    _ax.set_xlabel(\"line number\")\n",
    "    _ax.set_ylabel(\"% of time\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show relevant code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import inspect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   161\t                for p in pset:\n",
      "   162\t                    if p.dt*p.time <= p.dt*time:  # don't write particles if they haven't started yet\n",
      "   163\t                        i = p.fileid\n",
      "   164\t                        self.id[i, self.idx[i]] = p.id\n",
      "   165\t                        self.time[i, self.idx[i]] = time\n",
      "   166\t                        self.lat[i, self.idx[i]] = p.lat\n",
      "   167\t                        self.lon[i, self.idx[i]] = p.lon\n",
      "   168\t                        self.z[i, self.idx[i]] = p.depth\n",
      "   169\t                        for var in self.user_vars:\n",
      "   170\t                            getattr(self, var)[i, self.idx[i]] = getattr(p, var)\n"
     ]
    }
   ],
   "source": [
    "whole_file = !cat -n {str(inspect.getsourcefile(ParticleFile))}\n",
    "for line in whole_file[160:170]:\n",
    "    print(line)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:py3_parcels_profile]",
   "language": "python",
   "name": "conda-env-py3_parcels_profile-py"
  },
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}