{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute probability of transfer success from delays distributions\n",
    "\n",
    "To be able to compute the probability of success of a given transfert, we use the arrival delay distribution compared with the next trip departure. To be able to do that, we need delay distributions for each trip arrival to a given station. Whenever we have a clear match, we can use an __cumulative distribution function__ to compute $P(X \\leq x)$ :\n",
    "\n",
    "$${\\displaystyle F_{X}(x)=\\operatorname {P} (T\\leq t)=\\sum _{t_{i}\\leq t}\\operatorname {P} (T=t_{i})=\\sum _{t_{i}\\leq t}p(t_{i}).}$$\n",
    "\n",
    "The strategy was to rely entirely on past data we have to compute $p(t_i)$, without the need of building a model which imply making additionnal assumptions. If we have enough data for a given transfer with known trip_id x stop_id, we use the the abovementionned formula to compute each $p(t_i)$ by simply using :\n",
    "\n",
    "$$p(t_i) = \\frac{x_i}{\\sum x_i}$$\n",
    "\n",
    "with $x_i$ being the number of delays at time $t_i$ from SBB dataset.\n",
    "\n",
    "### Recover missing data \n",
    "\n",
    "As we are using SBB data to compute delays from timetable trip_id, we may encounter problems with the translation between the two datasets (certain trip_id/stop_id have no correspondance datasets!). We may also encounter To recover missing or faulty data, the strategy is the following :\n",
    "\n",
    "1. If we have more than 100 data points in `real` group, we rely exclusively on its delay distribution to compute probabilities for a given transfer on a `trip_id x stop_id`.\n",
    "\n",
    "_Note : `real` group corresponds to arrival time with status `geschaetz` or `real`, meaning it comes from actual measurments._\n",
    "\n",
    "2. If we do not find enough data within `real` group, we use delay distributions in `all` group (contains all delays including `prognose` status), if there is more than 100 data points for a given `trip_id x stop_id`.\n",
    "\n",
    "3. If `all` group still does not have more than 100 data points, we rely on `recovery tables` to estimate delay distributions. The strategy is the following :\n",
    "    - As we will always know the `stop_id`, the `time` and the `transport_type`, we rely on arrival delays from aggregated values of similar transfer. \n",
    "        - First, we compute a table of distribution with all possible combination of `stop_id`, `time` (round to hours) and `transport_type`, and aggregate all the counts we have to compute cumulative distribution probabilities. \n",
    "        - Is there is less than 100 data points in one of these intersections, we use the last possibilities : a table with `transport_type` x `time` aggregate counts.\n",
    "        - The last values with no match are given the overall average of cumulative distribution probabilities for each `transport_type` with no limit for the minimum number of data points.\n",
    "\n",
    "Following this approach, we can find cumulative distribution probabilities for every combination of `trip_id x stop_id` as defined in `stop_times_df`. We will make a table with the same row order so that McRaptor can easily find their indexes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 245,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle \n",
    "import gzip\n",
    "from itertools import islice\n",
    "import matplotlib as mlt \n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np \n",
    "import pandas as pd \n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functon to take a slice from a dictionnary - head equivalent\n",
    "def take(n, iterable):\n",
    "    \"Return first n items of the iterable as a list\"\n",
    "    return list(islice(iterable, n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load dictionnaries of distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "len dict_real :  6417\n",
      "[('10.TA.1-11-B-j19-1.1.R__8590314', array([ 0, 84, 54, 35, 24,  5,  3,  0,  0,  1,  1,  0,  0,  0,  0,  0,  0,\n",
      "        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0])), ('10.TA.1-11-B-j19-1.1.R__8590317', array([ 0, 96, 69, 36, 22,  5,  4,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,\n",
      "        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0])), ('10.TA.1-11-B-j19-1.1.R__8594304', array([ 0, 70, 78, 45, 22,  9,  1,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,\n",
      "        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0])), ('10.TA.1-11-B-j19-1.1.R__8594307', array([ 0, 87, 70, 43, 16,  5,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,\n",
      "        1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0])), ('10.TA.1-11-B-j19-1.1.R__8594310', array([ 0, 77, 65, 33, 13,  5,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,\n",
      "        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0]))]\n",
      "len dict_all :  221945\n",
      "[('1.TA.26-161-j19-1.1.H__8587347', array([  0, 214, 191,  71,  23,   5,   0,   1,   0,   0,   1,   1,   0,\n",
      "         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,\n",
      "         0,   0,   0,   0,   0,   0])), ('1.TA.26-161-j19-1.1.H__8590671', array([  0,  82, 155, 135,  86,  28,  12,   5,   1,   0,   2,   1,   0,\n",
      "         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,\n",
      "         0,   0,   0,   0,   0,   0])), ('1.TA.26-161-j19-1.1.H__8590677', array([  0, 118, 197, 113,  57,  16,   3,   1,   0,   0,   1,   1,   0,\n",
      "         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,\n",
      "         0,   0,   0,   0,   0,   0])), ('1.TA.26-161-j19-1.1.H__8590678', array([  0,  71, 140, 150,  88,  36,  12,   6,   1,   0,   2,   1,   0,\n",
      "         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,\n",
      "         0,   0,   0,   0,   0,   0])), ('1.TA.26-161-j19-1.1.H__8590680', array([  1, 216, 163,  74,  37,   9,   4,   0,   1,   0,   1,   1,   0,\n",
      "         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,\n",
      "         0,   0,   0,   0,   0,   0]))]\n"
     ]
    }
   ],
   "source": [
    "with gzip.open(\"../data/distributions_geschaetzAndReal.pkl.gz\", \"rb\") as input_file:\n",
    "    d_real = pickle.load(input_file)\n",
    "\n",
    "with gzip.open(\"../data/distributions_allDelays.pkl.gz\", \"rb\") as input_file:\n",
    "    d_all = pickle.load(input_file)\n",
    "\n",
    "# display a slice of it\n",
    "print('len dict_real : ', len(d_real))\n",
    "print(take(5, d_real.items()))\n",
    "\n",
    "# display a slice of it\n",
    "print('len dict_all : ', len(d_all))\n",
    "print(take(5, d_all.items()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probability using cumulative distribution based on frequency of delays \n",
    "\n",
    "When we have __enough data__ and no ambiguity about `trip_id` and `stop_id` for a given distribution, then we can compute the probability $P(T \\leq t)$ for every $t$ (delay in minute). \n",
    "\n",
    "Let's take a __threshold of 100__ sample points (=number of time we could measure a delay) as a minimum number of points to use this approach. \n",
    "\n",
    "_How many keys in our distionnary of distribution have at least this number of samples ?_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_data_points_hist(dico):\n",
    "    list_tot_points = []\n",
    "    for key in dico:\n",
    "        distrib = dico[key]\n",
    "        list_tot_points.append(np.sum(distrib))\n",
    "\n",
    "    tot_per_key = np.array(list_tot_points)\n",
    "    binwidth = 100\n",
    "    n_keys_less_than_binwidth = np.sum(np.array(tot_per_key < binwidth))\n",
    "    perc_key_to_recover = round(100 * ( n_keys_less_than_binwidth / len(tot_per_key) ), 2)\n",
    "    plt.figure(figsize = (10,5))\n",
    "    plt.hist(tot_per_key, bins = range(min(tot_per_key), max(tot_per_key) + binwidth, binwidth))\n",
    "    plt.title(\"Total number of data points per trip_id / stop_id key. N keys with less than {0} points: {1} ({2}%)\"\\\n",
    "             .format(binwidth, n_keys_less_than_binwidth, perc_key_to_recover))\n",
    "    plt.xlabel('n data points')\n",
    "    plt.ylabel('n keys')\n",
    "    return plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAokAAAFNCAYAAABsRvUeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd7gtVX3/8fdHioUiIFceKXJRMflhV35IYgmRSI1ibMGgICESY4kmpmCiYk0wxlh+YosiIDYkJhDFICKWqKAXC1XlCigg5UrHCvL9/bHWgc2ZU/blcso9vF/Ps589e60pa2bNrPnuNTN7p6qQJEmSRt1toQsgSZKkxccgUZIkSQMGiZIkSRowSJQkSdKAQaIkSZIGDBIlSZI0YJA4hiSV5EGLoBxfTPJnC7Tseyb57yTXJfnkGOPvkuSS+SjbmkiyX5LPLXQ55kuSzyY5YK7mkWR5P17WXZNlLJQkNyZ5wDR5z0/yv9PkLer1nqnsi0WS+/ftv84M44zdFid5bZJj7rwSrp61YZvPp7taWzudJF9N8qgFLsMWSc5LcvfZxl2rg8TeoEy8bknyi5HP+00zzVoRvCxCzwS2AO5TVc+6M2e8kI1pVX2kqnYbZ9zF3OiPe0Ksqj2r6qg1WdaaziPJlqt7DCY5Mskb7+gyx1VVG1bVBXO9nDUxEpCeOCn9mCSvXaBirbGq+nHf/r+Bhf1SvLoW+ktCkmcn+VqSnyf54hT5j0xyRs8/I8kjR/KS5M1JruqvNyfJnV3GhWpr+zn/lknxwgEj+f8nyRd6B8jKJH80afo/6+k3JvmfJFuO5P1VkguSXJ/kJ0neNtM+kOQpwA1V9e3++YBeH9cnuSTJv0w3fZInTFqHG/s+94yev2uSC5NcnmTfkek2SfKtJBtNpFXVFcCpwMGzbb+1OkjsDcqGVbUh8GPgKSNpH1no8i1WvVFY3brfFvhBVd08F2XS+O7IiegO1vlc2Qv4n4UuxBLw2CS/u9CF0KJwNfB24LDJGUnWB44HjgE2BY4Cju/p0AKFpwGPAB4OPAX483ko83z6yWi8MPElt7elxwOfBjajbYtjkjy45+8C/BOwT8+/EPjYyHxPAB5dVRsDD6Vtw7+coRwvBD488vlewMuBzYHHArsCfzPVhFX1lUkxzx8CN3JbW/p2Wt3tDrx7pEf+n4HDquqGSbP8COPUc1UtiRdwEfAHffjufYP9pL/e3tM2AH4B3NI37o3AlsBOwNeBa4HLgHcB64/Mu4AHTbPcLwJvAL4K3AB8Dti85+0CXDJDOV8LfJJ28N4AnAU8GHglcCVwMbDbpGX9M/AN4Hrazr3ZSP7OwNf6enwX2GXStG/q5fzFVOsD/J8+3rXAOcBTe/rrgF8DN/VtdtAU094TOBK4BjgX+NvRdQcOAX7Y1/Nc4I9GlvlL4Dd93tf29L2Bb/f1vBh47Qx1vwtwCfAPwE/7Nt5vJP/ewNHAKuBHwKuAu/W85wP/O6muXwic37fD4UBmKOdefX1uAC4F/maaMj6/b/t3AdcB3wN2nVTGD9L2v0uBNwLrTJr2bcBVwBsnzXuPSfXz3enqvKf92ThlmmF7j85jHeBf+3a/AHhx34brzjD9p4CnT5Gevo5X9no/i9bwHtzX7dd9/f57pv215x0JvBc4udfNl4Btx1i3W4914D60k8D1tGPuDaP7yqTplo+uN/AM2n74UNqX8Yn9/yrgWPpxC3wGeOmkeZ1JPz5mWdbfA6eOpB/DNMcJw/38LcD/9v1uyn0PWJ8WfDxsZLr7Aj8HltFObJ/u2/9q4Cv042rSsl8H/L8+vB7wM+AtI+3GL2kn4Fu3IW2//U3PuxF410zH5zTr/FrgmDHbx+fT9t8baIHAfj39QX3fuY62j39immX9uJdt4rzyOxPbnHZ8XNPnu+fINAcC5/VlXgD8+RRt2itox8NlwIFj7L9/BnxxUtpuvV4zqbx79OGvAQeP5B0EnLa2trXTlXmavIf2ZYxum88Bb+jD/wocPpK3ZS/3A6eY132AzwPvnmZZ69Pa4a1nKOtf09u3MdbrQ8CHRj5fMDJ8Oe1Y3Qn4n2mmX5d2LG8743LGKcza8OL2wdfrgdP6RlrWD4KJSh/sMMBjaA3IurSG6jzg5ZN25pmCxB/Sgrt79s+HzbCs0XK+th8Mu/dlH01rSP6R1pi+ALhw0rIu7Tv2BsB/0BtBYCvaCWgv2knpyf3zspFpfww8pC9rvUnlWg9YSTv41weeRDsYf2ukrMfMsP0Po50kNgO2Ac7m9kHis/oBdjfgj2knivv1vOcz6eTbt93D+vgPB64AnjbNsncBbgb+jfZl4Pf6/CfKfjQtoN6o1+8P6IHu5GX3uv40sAlwf1pjt8cM5bwMeEIf3pT2rXKqMj6/l/Gv+rb+Y9qJZyJY+E/gfb1e70sLSv580rQv7XV3zynmP6ifqeqcYZA4bZlmqOvRebyQFlxu0+v+VGYIEvtyfgpsNEXe7sAZfdtPnCwm9pEjGQmOmX1/PbJ/fmLfJ94xue6mKd9okPhxWkC3Ae2Yu3S6eXD7AOfAXraJ+byM1h5t3cvyPuBjPe/ZwOkj83kE7bhdf4YyTixro16mifZk1iCRdjz9O3AScK8x9r13A28emc/LuC1I/2daIL5efz2BKQK2Xjdn9eHfpbWXp4/kfXfyNpy8n41zfM50TDBD+9jX+/qRfed+wEP68Mdo7fHdgHsAj5+t/idt85to7fg6wF/QOi3S8/cGHkjb13+PdsJ+9KQ27fV92+7V8zedZf+dKkj8K+Czk9I+DbyiD18HPHYkb0faJdG1pq2lBZnT1c0utC+YV9DOr28DNuh5UwWJJwP/2Yf/lZGgr+9HBewzkvYnff+pvg6PmKYcDwF+Nkv9/Rc9fphlvA1o7dsuI2mn0dqPR/T9bD1a59eDZ5jPmYx8uZ7qtVguP93Z9gNeX1VXVtUq2jfZ5003clWdUVWnVdXNVXURrcH8vdVY3oeq6gdV9QvaSeWRs00w4itVdVK1y7ifpDVah1XVTbST1PIkm4yM/+GqOruqfga8Gnh271Z+LnBiVZ1YVbdU1cnAClrjMuHIqjqnr+dNk8qxM7BhX/avq+oLtAP4OWOux7OBN1XV1VV1MfDO0cyq+mRV/aSX7RO0b487TTezqvpiVZ3Vxz+T1ljPVievrqpfVdWXaD00E9tmX+CVVXVDr9+3MsP+QNsG11bVj2lBz0z1eROwQ5KNq+qaqvrWDONeCby9qm7q2+D7wN5JtqDV08ur6mdVdSWtIdt3ZNqfVNX/63X3ixm3wu3NVOfTlmk15v/sPv3FVXU1LXCYyRNpQcHkSx/QtuVGwG/TGu3zquqyaeYzzv76mar6clX9inai/50k24yzUn2/eQbwml4nZ9Mu083m5bRe9F2qamVPeyHwj1V1SS/La4Fn9ktdJwAPTrJ9H/d5tN6qX4+xrF/QetzGvVdzPdpxtBnt1pyfj7HvHQU8Z+Qetedx2+Wym2gB1bZ9//lK9TPPJF8Htk9yH1r9fxDYKsmGtGP6S2OWf8LqHJ8TZmsfbwEemuSeVXVZVZ0zso7bAltW1S+ranXvk/tRVf17tfssj6Jtry0AquozVfXDar5E68F6wsi0N9HOYzdV1Ym0YOa3VnP50I6T6yalXUc71qbKvw7YcJb7EhdVW1tVm8xQN9/r870f7UvJY2hBLrT27krgb5Osl2Q32j55r57/P33dHp7knsBraMHgRD5V9dFql5sfTPvSdMU05diEFthNKcmf0gL0f51unBFPp33ZHj12Xkj7Mvx+2jb/C1rP5j2SnJTk1CSTz6E39HJNa6kGiVvSuron/KinTSnJg5N8ut/weT3tHoTNV2N5l48M/5x20I1rdIf6BfDT3qBMfGbS/C4eGf4RreHfnNaQPSvJtRMv4PG0A2OqaSfbEri4qm6ZNP+txlyPLaco262S7J/kOyNleygzbOMkj+079aok19EOgJnq5JoeOI8uf8s+zXoM94eZ1mt16vMZtBPNj5J8KcnvzDDupZNOohNl3LaX8bKR7fM+Wq/OhJnqbiazTTddmcY1Y71PYS/gxKkyeqD3LtplpyuTvD/JxjMtd5b99dZyVdWNtEui467bMlqv4OqsG7QA8fCqGn0wZ1vgP0fq9jzapbQtquqXwCeA5/Z7Rp/D7e9Zms0HgC36DfGzeRDt3qrXjQShM+57VXU67RjYJclv93mc0Kd9C63H9HP95v1Dplpo/1KzgnbyfSLtxPY14HHcsSDxjrS307aPvd34Y1obc1mSz/R1Bfg7Wk/fN5Kc00/kd6isVfXzPrghQJI9k5yW5Openr24fRt3Vd3+HvDVPbdMuBGYfBxtzG0By+T8jYEbpwn4Ye1oa29VVZdX1bn9y8GFtDp9Rs+7iXY/5t69LK+gdfRc0vM/DxxKu2p3UX/dMJE/aTnn0257efc0RbmG2wLz20nyNNoX7D2r6qdjrNYBwNGjdVRV36mqXarqsbTL8n9Ki2U+QOsoOxD48KTgfyNaL+y0lmqQ+BNaozDh/j0N2reAyd5D+7axff9G8A+0hmFN/YyRbxz9m9ayNZznaE/I/Wnfrn5KO5l9uH+jmnhtUFWjNzJPd9BD2z7bTHq44f60y1njuGyKsgGQZFvaJa6X0J6O3oR2OXpiG09Vro/STkbbVNW9ad/QZqqTTZNsMGn5P6Ftm4negNG8cddr1KCcVfXNqtqHdlL9L1oDM52tJh2gE2W8GPgV7V7WibrbuKoeMtOyZyvbmNNNV6ZxTVvv05g2SASoqndW1WOAHWjfzP92ImvSqOPsr7eWq/dabcb467aKdlltddYN2v1fr5p44rC7mNb4jx6b96iqibIeRbv6sSvw86r6+phlpAd7r6PdLzlbm3Ue7UTx2SQTPVLj7HtH0Xringcc1wNbem/RK6rqAcBTgb9Osus0y/4SrRfnUcA3++fdaVcTvjzd6s2yPqtjxvax2tWcJ9O+VH+P1l5NBBgvqKotaTf5vztT/wTPapU17adH/oPWa7RFbxNP5M4570x2DvDwScf5w3v6RP4jRvIeMZI3lbWhrZ1t3re2G1V1ZlX9XlXdp6p2Bx5Au+ViIv/wqtq+qrag1dm6tPPXVNal3UIwlZW0ZwhvFzQn2YO2vz2lqs6arfD9asgutEv703kb8Kr+Be1hwIres7sePQbpVzIeRLs/d1pLNUj8GK2hXpZkc1oX8cTPg1wB3CfJvUfG34h2T8GN/RvkX9xJ5fgBrat37yTr0W7infV3iWbx3CQ7JLkX7X6V43rP4zHAU5LsnmSdJPdIe/R/6zHnO9Fj8He9230X2pNSHx9z+mOBVybZtC/zpSN5G3Db/RokOZDWkzjhCmDrkaftoNXJ1VX1yyQ70e77mM3rkqyf5Am0J78+2bfNscCbkmzUA9a/5rb9YXXcrpx9WfsluXf/Rno97bLVdO4L/GXfvs+i3XN3YrVLqp8D3ppk4yR3S/LAKS4NzFa25Vn9J5inLNNqTH9sn37rJJvSHtCYUpLtgLtX1XnT5P/f3oM88XDDL7lte15Ba7wnjLO/7pXk8b2+3kC7GX+sHtm+33wKeG2SeyXZgfbtfTbn0B4kOjzJU3vae2n737Z9PZcl2WdkWV/v6/lWVq8XccKHaffL7THbiFX1MdqX4M8neeCY+94xwB/RAsVbT0xJ/jDJg3rwcR2td3S6/f9LwP7AuT2w/SLt/rkLq90SNJXJdb4mpm0f034zbp8e+PyK1rN2S1/HZ420odfQ2rGp1nFVTx+3vOvTzgWrgJuT7En7gnGHTKwTLUi5W1+/9Xr2F2l185dJ7p7kJT39C/39aFqAv1Xaz7u8gnZP70wWe1t7qyS/n2TbNNvQ7p8/fiT/4X173SvJ39C+KBzZ8+6R5KF92vvTLuW+o6qu6fl/luS+fXgH2kOnp0xVjr7ff56R26aSPIn2lPEzquobU003hecBX6uqH06zvk8G7lFVn+5JFwJPSvIQ2j53VU/fCbioqma8QrJUg8Q30i5vnEl7QvJbPY2q+h4tiLwg7bLDlrRHzv+E1o3877TLP2usqq4DXkTr7r2UduJb099o/DBtB76cdmL4y76si2mXkv6B1vBcTOuFGauO+w78FGBP2jfCdwP79+01jtfRLi1cSDvp3Hqyq6pzaSfAr9MO/ofRnqqd8AXayfXyJBNd7S8CXp/kBlqQP9u3xstpjfhPaAfdC0fK/lLatr+AdvP+R4EjxlyvUVOV83nARWm3KbyQ1iM0ndOB7Wnb903AM6tq4oDdn3biOLevx3Hc/laB2Uz8wPlVSWa6L3J1yjSOiYcgvks7zj41w7h7M3MAunGf3zW0fekq2iVNaPex7dCP2f8ac3/9KO1S0dW0+5CeuxrrBa3ne0PavnUk7WnCWVXVd2knzn/vJ/930HrFP9f359NoP3cx6mjacXEMQJL3JnnvmMv7De0Y2WzM8Y+ifcH8QpLlzLLv9bblW7QA6Ssjs9qedtK7kXZsv7uqTp1msV+jPdg30Wt4Lu1LwHS9iNC22zOTXJPknTOMN6tZ2se70YKZn9D2ld/jto6C/wucnuRGWh2+rKb4Hc1+KflNwFf7PrrzLOW5gdZ2H0vb5n/CbZfx74jn0W5Peg/tvsZfcFtv6K9pl1T3p11a/FPaQ4ATtxy8D/hv2rnybNo9hu+bYVmLrq1N+83AJwxnA7Te66/1cn2Ntp6jP1PzPNoVkStpvflPrnbvMLRz7Edp+/g3aPv5q0emfRxwVpKf0dq2E2n72HTex+3v0Xw17YnwE3Pbbx9+dmS9Pptk8vz2Z5r7o9N6qN9Ce8BswktpX1Q/D7yobrudbb+ePqOJp6yktVbvRTqmqsbtNZ13SZ5Pe1Lz8QtdlgnzXaa0H39+V7Wb8Od6WUfSnq5/1Vwv686QZH/az5Asmv1jVJIjaA9PrRXbU3NjbWhrF7skXwVeUv0HtReoDPel9e4/auL2keksyr+QkrQkfZH2BKNGpN068iKmv+F9QfXexqfTemQkrYGqetwiKMOVtFuLZrVULzdLWgMZ/v3TjbNc0plVVf1Lrd7P99zpMvVfW93YLycuRHl2p13+vIJ2WWtRSfIG2iXIt/QnQyXdhXi5WZIkSQP2JEqSJGnAIFGSJEkDPrjSbb755rV8+fKFLoYkSdKszjjjjJ9W1Zr+QceMDBK75cuXs2LFioUuhiRJ0qySjPNXoWvEy82SJEkaMEiUJEnSgEGiJEmSBgwSJUmSNGCQKEmSpAGDREmSJA0YJEqSJGnAIFGSJEkDBomSJEkaMEiUJEnSgEGiJEmSBvzvZg0sP+QzY4130WF7z3FJJEnSQrEnUZIkSQMGiZIkSRowSJQkSdKAQaIkSZIG5ixITHJEkiuTnD2StlmSk5Oc39837elJ8s4kK5OcmeTRI9Mc0Mc/P8kBI+mPSXJWn+adSTLTMiRJkjS+uexJPBLYY1LaIcApVbU9cEr/DLAnsH1/HQy8B1rABxwKPBbYCTh0JOh7D/CCken2mGUZkiRJGtOcBYlV9WXg6knJ+wBH9eGjgKeNpB9dzWnAJknuB+wOnFxVV1fVNcDJwB49b+OqOq2qCjh60rymWoYkSZLGNN/3JG5RVZf14cuBLfrwVsDFI+Nd0tNmSr9kivSZljGQ5OAkK5KsWLVq1R1YHUmSpKVpwR5c6T2AtZDLqKr3V9WOVbXjsmXL5rIokiRJa5X5DhKv6JeK6e9X9vRLgW1Gxtu6p82UvvUU6TMtQ5IkSWOa7yDxBGDiCeUDgONH0vfvTznvDFzXLxmfBOyWZNP+wMpuwEk97/okO/enmvefNK+pliFJkqQxzdl/Nyf5GLALsHmSS2hPKR8GHJvkIOBHwLP76CcCewErgZ8DBwJU1dVJ3gB8s4/3+qqaeBjmRbQnqO8JfLa/mGEZkiRJGtOcBYlV9ZxpsnadYtwCXjzNfI4AjpgifQXw0CnSr5pqGZIkSRqf/7giSZKkAYNESZIkDRgkSpIkacAgUZIkSQMGiZIkSRowSJQkSdKAQaIkSZIGDBIlSZI0YJAoSZKkAYNESZIkDRgkSpIkacAgUZIkSQMGiZIkSRowSJQkSdKAQaIkSZIGDBIlSZI0YJAoSZKkAYNESZIkDRgkSpIkacAgUZIkSQMGiZIkSRowSJQkSdKAQaIkSZIGDBIlSZI0YJAoSZKkAYNESZIkDRgkSpIkacAgUZIkSQMGiZIkSRowSJQkSdKAQaIkSZIGDBIlSZI0YJAoSZKkAYNESZIkDRgkSpIkacAgUZIkSQMGiZIkSRowSJQkSdKAQaIkSZIGDBIlSZI0YJAoSZKkgQUJEpP8VZJzkpyd5GNJ7pFkuySnJ1mZ5BNJ1u/j3r1/Xtnzl4/M55U9/ftJdh9J36OnrUxyyPyvoSRJ0tpt3oPEJFsBfwnsWFUPBdYB9gXeDLytqh4EXAMc1Cc5CLimp7+tj0eSHfp0DwH2AN6dZJ0k6wCHA3sCOwDP6eNKkiRpTAt1uXld4J5J1gXuBVwGPAk4rucfBTytD+/TP9Pzd02Snv7xqvpVVV0IrAR26q+VVXVBVf0a+HgfV5IkSWOa9yCxqi4F/hX4MS04vA44A7i2qm7uo10CbNWHtwIu7tPe3Me/z2j6pGmmSx9IcnCSFUlWrFq1as1XTpIkaYlYiMvNm9J69rYDtgQ2oF0unndV9f6q2rGqdly2bNlCFEGSJGlRWojLzX8AXFhVq6rqJuBTwOOATfrlZ4CtgUv78KXANgA9/97AVaPpk6aZLl2SJEljWogg8cfAzknu1e8t3BU4FzgVeGYf5wDg+D58Qv9Mz/9CVVVP37c//bwdsD3wDeCbwPb9aen1aQ+3nDAP6yVJkrRkrDv7KHeuqjo9yXHAt4CbgW8D7wc+A3w8yRt72gf7JB8EPpxkJXA1Leijqs5JciwtwLwZeHFV/QYgyUuAk2hPTh9RVefM1/pJkiQtBfMeJAJU1aHAoZOSL6A9mTx53F8Cz5pmPm8C3jRF+onAiWteUkmSpLsm/3FFkiRJAwaJkiRJGjBIlCRJ0oBBoiRJkgYMEiVJkjRgkChJkqQBg0RJkiQNGCRKkiRpwCBRkiRJAwaJkiRJGjBIlCRJ0oBBoiRJkgYMEiVJkjRgkChJkqQBg0RJkiQNGCRKkiRpwCBRkiRJAwaJkiRJGjBIlCRJ0oBBoiRJkgYMEiVJkjRgkChJkqQBg0RJkiQNGCRKkiRpwCBRkiRJAwaJkiRJGjBIlCRJ0oBBoiRJkgYMEiVJkjRgkChJkqQBg0RJkiQNGCRKkiRpwCBRkiRJAwaJkiRJGjBIlCRJ0oBBoiRJkgYMEiVJkjRgkChJkqQBg0RJkiQNGCRKkiRpwCBRkiRJAwsSJCbZJMlxSb6X5Lwkv5NksyQnJzm/v2/ax02SdyZZmeTMJI8emc8Bffzzkxwwkv6YJGf1ad6ZJAuxnpIkSWurhepJfAfwP1X128AjgPOAQ4BTqmp74JT+GWBPYPv+Ohh4D0CSzYBDgccCOwGHTgSWfZwXjEy3xzyskyRJ0pIx70FiknsDTwQ+CFBVv66qa4F9gKP6aEcBT+vD+wBHV3MasEmS+wG7AydX1dVVdQ1wMrBHz9u4qk6rqgKOHpmXJEmSxjBrkJjkcUk26MPPTfJvSbZdg2VuB6wCPpTk20k+0Oe/RVVd1se5HNiiD28FXDwy/SU9bab0S6ZIlyRJ0pjG6Ul8D/DzJI8AXgH8kNY7d0etCzwaeE9VPQr4GbddWgag9wDWGixjLEkOTrIiyYpVq1bN9eIkSZLWGuMEiTf3oG0f4F1VdTiw0Ros8xLgkqo6vX8+jhY0XtEvFdPfr+z5lwLbjEy/dU+bKX3rKdIHqur9VbVjVe24bNmyNVglSZKkpWWcIPGGJK8Engt8JsndgPXu6AKr6nLg4iS/1ZN2Bc4FTgAmnlA+ADi+D58A7N+fct4ZuK5flj4J2C3Jpv2Bld2Ak3re9Ul27k817z8yL0mSJI1h3THG+WPgT4CDquryJPcH3rKGy30p8JEk6wMXAAfSAtZjkxwE/Ah4dh/3RGAvYCXw8z4uVXV1kjcA3+zjvb6qru7DLwKOBO4JfLa/JEmSNKZxgsRnAR/qTxBTVT9mze5JpKq+A+w4RdauU4xbwIunmc8RwBFTpK8AHromZZQkSborG+dy8xbAN5Mcm2QPf5hakiRp6Zs1SKyqV9F+kPqDwPOB85P8U5IHznHZJEmStEDG+jHtfsn38v66GdgUOC7Jv8xh2SRJkrRAZr0nMcnLaE8I/xT4APC3VXVTf8r5fODv5raIkiRJmm/jPLiyGfD0qvrRaGJV3ZLkD+emWJIkSVpI49yTeCiwTZIDAZIsS7JdzztvjssnSZKkBTDOfzcfCvw98MqetB5wzFwWSpIkSQtrnAdX/gh4Ku0/lqmqn7Bmf8snSZKkRW6cIPHX/enmAkiywdwWSZIkSQttnCDx2CTvAzZJ8gLg87SnnCVJkrREjfN081uBPwCuB34LeA3w5bkslCRJkhbWOEHiB6vqT4GTAZJsCJzIFP+zLEmSpKVhnMvNlyZ5N0CSTYHP4dPNkiRJS9o4v5P4auDGJO+lBYhvraoPzXnJJEmStGCmvdyc5OkjH08HXg18A6gkT6+qT8114SRJkrQwZron8SmTPn+b9kPaT6H9HI5BoiRJ0hI1bZBYVQfOZ0EkSZK0eIzz4IokSZLuYgwSJUmSNGCQKEmSpIFZf0w7yd2BZwDLR8evqtfPXbEkSZK0kMb5x5XjgeuAM4BfzW1xJEmStBiMEyRuXVV7zHlJJEmStGiMc0/i15I8bM5LIkmSpEVjnJ7ExwPPT3Ih7XJzgKqqh89pySRJkrRgxgkS95zzUkiSJGlRmTVIrKofzUdBJEmStHj4O4mSJEkaMEiUJEnSgEGiJEmSBgwSJUmSNGCQKEmSpAGDREmSJA0YJEqSJGnAIFGSJEkDBomSJEkaMEiUJEnSgEGiJEmSBgwSJUmSNGCQKEmSpAGDREmSJA0sWJCYZJ0k307y6f55uySnJ1mZ5BNJ1u/pd++fV/b85SPzeGVP/36S3UfS9+hpK5McMt/rJkmStLZbyJ7ElwHnjXx+M/C2qq9tS9UAAA3PSURBVHoQcA1wUE8/CLimp7+tj0eSHYB9gYcAewDv7oHnOsDhwJ7ADsBz+riSJEka04IEiUm2BvYGPtA/B3gScFwf5SjgaX14n/6Znr9rH38f4ONV9auquhBYCezUXyur6oKq+jXw8T6uJEmSxrRQPYlvB/4OuKV/vg9wbVXd3D9fAmzVh7cCLgbo+df18W9NnzTNdOmSJEka07wHiUn+ELiyqs6Y72VPUZaDk6xIsmLVqlULXRxJkqRFYyF6Eh8HPDXJRbRLwU8C3gFskmTdPs7WwKV9+FJgG4Cef2/gqtH0SdNMlz5QVe+vqh2rasdly5at+ZpJkiQtEfMeJFbVK6tq66paTnvw5AtVtR9wKvDMPtoBwPF9+IT+mZ7/haqqnr5vf/p5O2B74BvAN4Ht+9PS6/dlnDAPqyZJkrRkrDv7KPPm74GPJ3kj8G3ggz39g8CHk6wErqYFfVTVOUmOBc4FbgZeXFW/AUjyEuAkYB3giKo6Z17XRJIkaS23oEFiVX0R+GIfvoD2ZPLkcX4JPGua6d8EvGmK9BOBE+/EokqSJN2l+I8rkiRJGjBIlCRJ0oBBoiRJkgYMEiVJkjRgkChJkqQBg0RJkiQNGCRKkiRpwCBRkiRJAwaJkiRJGjBIlCRJ0oBBoiRJkgYMEiVJkjRgkChJkqQBg0RJkiQNGCRKkiRpwCBRkiRJAwaJkiRJGjBIlCRJ0oBBoiRJkgYMEiVJkjSw7kIXQGuv5Yd8ZqzxLjps7zkuiSRJurPZkyhJkqQBg0RJkiQNGCRKkiRpwCBRkiRJAwaJkiRJGjBIlCRJ0oBBoiRJkgYMEiVJkjRgkChJkqQB/3FFS9q4/woD/jOMJEmj7EmUJEnSgEGiJEmSBgwSJUmSNGCQKEmSpAGDREmSJA0YJEqSJGnAIFGSJEkDBomSJEkaMEiUJEnSgEGiJEmSBuY9SEyyTZJTk5yb5JwkL+vpmyU5Ocn5/X3Tnp4k70yyMsmZSR49Mq8D+vjnJzlgJP0xSc7q07wzSeZ7PSVJktZmC9GTeDPwiqraAdgZeHGSHYBDgFOqanvglP4ZYE9g+/46GHgPtKASOBR4LLATcOhEYNnHecHIdHvMw3pJkiQtGfMeJFbVZVX1rT58A3AesBWwD3BUH+0o4Gl9eB/g6GpOAzZJcj9gd+Dkqrq6qq4BTgb26HkbV9VpVVXA0SPzkiRJ0hgW9J7EJMuBRwGnA1tU1WU963Jgiz68FXDxyGSX9LSZ0i+ZIl2SJEljWrAgMcmGwH8AL6+q60fzeg9gzUMZDk6yIsmKVatWzfXiJEmS1hoLEiQmWY8WIH6kqj7Vk6/ol4rp71f29EuBbUYm37qnzZS+9RTpA1X1/qrasap2XLZs2ZqtlCRJ0hKyEE83B/ggcF5V/dtI1gnAxBPKBwDHj6Tv359y3hm4rl+WPgnYLcmm/YGV3YCTet71SXbuy9p/ZF6SJEkaw7oLsMzHAc8DzkrynZ72D8BhwLFJDgJ+BDy7550I7AWsBH4OHAhQVVcneQPwzT7e66vq6j78IuBI4J7AZ/trSVp+yGfGGu+iw/ae45JIkqSlZN6DxKr6X2C63y3cdYrxC3jxNPM6AjhiivQVwEPXoJiSJEl3af7jiiRJkgYMEiVJkjRgkChJkqQBg0RJkiQNGCRKkiRpwCBRkiRJAwaJkiRJGjBIlCRJ0oBBoiRJkgYMEiVJkjRgkChJkqQBg0RJkiQNGCRKkiRpwCBRkiRJAwaJkiRJGjBIlCRJ0oBBoiRJkgYMEiVJkjRgkChJkqSBdRe6ANJisfyQz4w13kWH7T3HJZEkaeHZkyhJkqQBg0RJkiQNGCRKkiRpwCBRkiRJAwaJkiRJGjBIlCRJ0oBBoiRJkgb8ncS7iHF/A1CSJAnsSZQkSdIUDBIlSZI04OXmebQ6l3z96zdJkrSQ7EmUJEnSgD2JmnP2oEqStPYxSNSi4lPYkiQtDl5uliRJ0oA9iYuUPWqSJGkh2ZMoSZKkAYNESZIkDRgkSpIkacAgUZIkSQMGiZIkSRpYskFikj2SfD/JyiSHLHR5JEmS1iZLMkhMsg5wOLAnsAPwnCQ7LGypJEmS1h5LMkgEdgJWVtUFVfVr4OPAPgtcJkmSpLXGUv0x7a2Ai0c+XwI8doHKoiVmbfih87XlP7DnYluuLesuSYvdUg0Sx5LkYODg/vHGJN+f40VuDvx0jpehO27J1E/evNAluNONXTdLcN3XBkvm2FmCrJvFbU3qZ9s7syBTWapB4qXANiOft+5pt1NV7wfeP1+FSrKiqnacr+Vp9Vg/i5d1s7hZP4uXdbO4Lfb6War3JH4T2D7JdknWB/YFTljgMkmSJK01lmRPYlXdnOQlwEnAOsARVXXOAhdLkiRprbEkg0SAqjoROHGhyzHJvF3a1h1i/Sxe1s3iZv0sXtbN4rao6ydVtdBlkCRJ0iKzVO9JlCRJ0howSJwn/k3g/EtyRJIrk5w9krZZkpOTnN/fN+3pSfLOXj9nJnn0yDQH9PHPT3LAQqzLUpNkmySnJjk3yTlJXtbTrZ9FIMk9knwjyXd7/byup2+X5PReD5/oDwaS5O7988qev3xkXq/s6d9PsvvCrNHSk2SdJN9O8un+2bpZJJJclOSsJN9JsqKnrZ1tW1X5muMX7eGZHwIPANYHvgvssNDlWuov4InAo4GzR9L+BTikDx8CvLkP7wV8FgiwM3B6T98MuKC/b9qHN13odVvbX8D9gEf34Y2AH9D+QtP6WQSvvp037MPrAaf37X4ssG9Pfy/wF334RcB7+/C+wCf68A69vbs7sF1vB9dZ6PVbCi/gr4GPAp/un62bRfICLgI2n5S2VrZt9iTOD/8mcAFU1ZeBqycl7wMc1YePAp42kn50NacBmyS5H7A7cHJVXV1V1wAnA3vMfemXtqq6rKq+1YdvAM6j/VOS9bMI9O18Y/+4Xn8V8CTguJ4+uX4m6u04YNck6ekfr6pfVdWFwEpae6g1kGRrYG/gA/1zsG4Wu7WybTNInB9T/U3gVgtUlru6Larqsj58ObBFH56ujqy7OdYvfz2K1ltl/SwS/XLmd4AraSeoHwLXVtXNfZTRbX1rPfT864D7YP3MlbcDfwfc0j/fB+tmMSngc0nOSPtnN1hL27Yl+xM40myqqpL4eP8CSrIh8B/Ay6vq+tbB0Vg/C6uqfgM8MskmwH8Cv73ARRKQ5A+BK6vqjCS7LHR5NKXHV9WlSe4LnJzke6OZa1PbZk/i/BjrbwI1L67oXfn09yt7+nR1ZN3NkSTr0QLEj1TVp3qy9bPIVNW1wKnA79AuhU10Loxu61vroeffG7gK62cuPA54apKLaLcuPQl4B9bNolFVl/b3K2lfsHZiLW3bDBLnh38TuHicAEw8JXYAcPxI+v79SbOdgev6pYGTgN2SbNqfRtutp2kN9HuiPgicV1X/NpJl/SwCSZb1HkSS3BN4Mu2+0VOBZ/bRJtfPRL09E/hCtbvvTwD27U/YbgdsD3xjftZiaaqqV1bV1lW1nHYu+UJV7Yd1sygk2SDJRhPDtDbpbNbWtm2+n5S5q75oTzD9gHZfzz8udHnuCi/gY8BlwE20+zkOot2LcwpwPvB5YLM+boDDe/2cBew4Mp8/pd3UvRI4cKHXaym8gMfT7ts5E/hOf+1l/SyOF/Bw4Nu9fs4GXtPTH0ALJFYCnwTu3tPv0T+v7PkPGJnXP/Z6+z6w50Kv21J6Abtw29PN1s0iePV6+G5/nTNxvl9b2zb/cUWSJEkDXm6WJEnSgEGiJEmSBgwSJUmSNGCQKEmSpAGDREmSJA0YJErSJElunCV/kyQvmodyvD7JH8wyzi5JfneuyyLprscgUZJW3ybAnAeJVfWaqvr8LKPtAhgkSrrTGSRKWrKSLE9yXpJ/T3JOks/1fxCZPN52Sb6e5KwkbxxJ3zDJKUm+1fP26VmHAQ9M8p0kb5lhvMnLuTHJ23pZTkmyrKc/MslpSc5M8p/9HxZIcmSSZ/bhi5K8bmQZv51kOfBC4K96WZ6Q5FlJzk7y3SRfvjO3p6S7FoNESUvd9sDhVfUQ4FrgGVOM8w7gPVX1MNq/9Ez4JfBHVfVo4PeBt/a/FDwE+GFVPbKq/naG8SbbAFjRy/Il4NCefjTw91X1cNq/Lhw6xbQAP+3LeA/wN1V1EfBe4G29LF8BXgPsXlWPAJ4669aRpGkYJEpa6i6squ/04TOA5VOM8zja3zgCfHgkPcA/JTmT9ldaWwFbTDH9uOPdAnyiDx8DPD7JvYFNqupLPf0o4InTrMunZlkPgK8CRyZ5AbDONONI0qzWXegCSNIc+9XI8G+AweXmbqr/KN0PWAY8pqpuSnIR7b9w7+h44yxzJhPr8humab+r6oVJHgvsDZyR5DFVddVqLkeS7EmUJFrv2759eL+R9HsDV/bA7/eBbXv6DcBGY4w32d2AZ/bhPwH+t6quA65J8oSe/jzapehx3a4sSR5YVadX1WuAVcA2qzEvSbqVPYmSBC8DPprk74HjR9I/Avx3krOAFcD3AKrqqiRfTXI28FngzVONN4WfATsleRVwJfDHPf0A4L1J7gVcABy4GmX/b+C4/rDMS2kPsWxPuwR+CvDd1ZiXJN0qVat7tUOSdEckubGqNlzockjSOLzcLEmSpAF7EiVJkjRgT6IkSZIGDBIlSZI0YJAoSZKkAYNESZIkDRgkSpIkacAgUZIkSQP/H3sCKcfZuYy+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_data_points_hist(d_all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 250,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoEAAAFNCAYAAAB/kbXqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZgmVXn///dHFhdQFpnwk3WIEv3iGkOQxI2ERFnUMa4YFSQYYtyNiaIx4v7FGOPyVTEoCKgBCTERBaO4axQENCKLxpF1EGRYBUEFvX9/nNPw0PQ6M9093fV+XVdfXXXqVNVd+/2cqnqeVBWSJEkalrssdACSJEmafyaBkiRJA2QSKEmSNEAmgZIkSQNkEihJkjRAJoGSJEkDNIgkMEklud96EMdXkjx/geZ99ySfTnJ9kn+bQf09kqyaj9jWRpJnJ/n8QscxX5J8NskBczWNJMv78bLh2sxjoSS5MclvTzLseUm+Mcmw9Xq5p4p9fZFkh77+N5iizozPxUnekORj6y7C2VkM63w+DeVcm+Svkrx7oeNYG0nemeSvZ1J3QZPAfsIY+/tNkptH+p89yTiLIjlZDz0N2Bq4d1U9fV1OeCFPllX18ap63Ezqrs8n9Zle8Kpq76o6Zm3mtbbTSLLNbI/BJEcnecuaznOmqmrTqrpgruezNkYSzlPGlX8syRsWKKy1VlWX9PX/a1jYD72ztdAfApI8I8k3k9yU5CsTDH9YkrP68LOSPGxkWJK8PcnV/e/tSbKuY1yoc22S147LFW7u+cJWffhdkxyV5GdJrkjyNyPjbpzkxCQX9e27xzTz2hh4HfCOkbIjkvywz/N5E4zzij7fn/U47jrF9J+fZGVfjv9Kss3IsL9Lck6SG5JcmOTvRoZtmOT4JNf18e41bv38zbhZ/RPw2r48U1rQJLCfMDatqk2BS4AnjpR9fCFjW5/1g362225H4H+r6ta5iEkztyYXmjXc5nNlH+C/FjqIJeARSf5woYPQeuEa4N3AYeMH9Av5p4CPAVsAxwCfGrnAHww8GXgo8BDgicBfzUPM86Kq3jYuV3g78JWquqpXeQOwM+0a90fAq5LsNTKJbwDPAa6YwexWAD+oqstGyr4HvBD4zvjKSR4PHALs2ef/28AbJ5pwT0Df1uexJXAhcNxoFWB/2jbeC3hxkv36sKcABWwFXE/b5iTZCXgS8N7ReVXV5cAP+rCpVdV68QdcBPxJ774r7YD4Sf97dy/bBLgZ+A1wY//bBtgN+BZwHXA58D5g45FpF3C/Seb7FeDNwH8DNwCfB7bqw/YAVk0R5xuAf6MdnDcA3wd+B3gNcCVwKfC4cfP6v8C3gZ/RDuwtR4bvDnyzL8f3gD3GjfvWHufNEy0P8H96veuAc4En9fI3Ar8Cbunr7KAJxr07cDRwLXAe8Hejy07b0X/cl/M84M9G5vkL4Nd92tf18n2B7/blvBR4wxTbfg9gFfBa4Kq+jp89Mnwz4FhgNXAx7ZPaXfqw5wHfGLetXwD8qK+H99MOrsni3Kcvzw3AZcDfThLj8/q6fx/tIPwBsOe4GI+k7X+XAW8BNhg37ruAq4G3jJv2XuO2z/cm2+a97PkziWmK9T06jQ1onxqvAi4AXtTX4YZTjP9J4CkTlKcv45V9u38feBDthHVLX8YbgU9Ptb/2YUcDHwRO7dvmq8COM1i224514N7AST2Wb9OO829MMt7y0eUGnkrbDx9E+7A8tv9fDZxAP26Bk4GXjJvW2fTjY5p5vRr48kj5x5jkOOHO+/k7aBe3zSbb94CNacnFg0fG+y3gJmAZ7YLymb7+rwG+Tj+uxs37jcD/690bAT8H3jFy3vgF7aJ22zqk7be/7sNuBN431fE5yTK/AfjYDM+Pz6PtvzfQLq7P7uX36/vO9bR9/BOTzOuSHtvYdeUPxtY57fi4tk9375FxDgTO7/O8APirCc5pr6QdD5cDB85g/30+LcEZLXtc364ZF+9evfubwMEjww4CTlus59pp1k/6uj5gpOwn3PE6+2bg+AnGXTW6z0wy/aOA100y7BvA88aV/SvwtpH+PYErJhn/n4D3j/Rv09fhfSep/15uP+5ePbZ/9XX+gd79aeCRk4z/98BHpl2ns90Ic/XHHZOrNwGn0U5Yy/pO/ubRnXjcuL9HO0FsSDsRnQ+8fNzOOlUS+GNa8nb33n/YFPMajfMNfWd/fJ/3sbQTxd/TTpZ/CVw4bl6X0S4smwD/Tj/JAdvSLjD70C46f9r7l42MewnwwD6vjcbFtRGwknZwbwz8Me1gu/9IrB+bYv0fRrsIbAlsD5zDHZPAp/ed9i7AM2kXgvv0Yc9j3MW1r7sH9/oPAX4KPHmSee8B3Ar8My3Zf2yf/ljsx9IS5nv27fu/9ER2/Lz7tv4MsDmwA+1kttcUcV4OPLp3bwE8fJIYn9djfEVf18+kXVjGkoH/AP6lb9ffoiUdfzVu3Jf0bXf3CaZ/p+0z0TbnzkngpDFNsa1Hp/ECWvK4fd/2X2aKJLDP5yrgnhMMezxwVl/3YxeDsX3kaEaSX6bfX4/u/Y/p+8R7xm+7SeIbTQKPpyVsm9COucsmmwZ3TGAO7LGNTedltPPRdj2WfwGO68OeAZw+Mp2H0o7bjaeIcWxe9+wxjZ1Ppk0CacfTh4DPAfeYwb73AeDtI9N5Gbcn4f+Xlmhv1P8ezQQJWd823+/df0g7X54+Mux749fh+P1sJsfnVMcEU5wf+3L/bGTfuQ/wwN59HO18fBfgbsCjptv+49b5LbTz+AbAX9MSjvTh+wL3pe3rj6Ul1w8fd057U1+3+/ThW0yz/06UBL4C+Oy4ss8Ar+zd1wOPGBm2K3DDJNMfi2u9OtfSksgJt824aTyGllhuOjKdArYeqfM0+v46btyZJIFnAE+fZNhESeD3gGeO9G/V47n3BOP/Ez15G9mnC1gxQd3QGlFeMLKvfaJvs0/QPqz/GVMkebTWw+9Mt07Xl9tL4z0beFNVXVlVq2mfRJ87WeWqOquqTquqW6vqItoJ8bGzmN9Hqup/q+pm2kXjYdONMOLrVfW5ardZ/412Ujqsqm6hXYSWJ9l8pP5Hq+qcqvo58A/AM9IepH4OcEpVnVJVv6mqU4EzaSePMUdX1bl9OW8ZF8fuwKZ93r+qqi/RDtBnzXA5ngG8taquqapLuXPz8r9V1U96bJ+gffrbbbKJVdVXqur7vf7ZtJPxdNvkH6rql1X1VVoLy9i62Q94TVXd0LfvO5lif6Ctg+uq6hJaUjPV9rwF2CXJvarq2qq6U5P/iCuBd1fVLX0d/BDYN8nWtO308qr6eVVdSWsR229k3J9U1f/r2+7mKdfCHU21zSeNaRbTf0Yf/9KquoaWGEzlMbSL/g0TDLuFdvF4AO1CeX612xITmcn+enJVfa2qfkm7kP9Bku1nslB9v3kq8Pq+Tc6h3UabzstpreB7VNXKXvYC4O+ralWP5Q3A0/pt/ZOA30myc6/7XFpr069mMK+baS1mM31WciPacbQl7dGZm2aw7x0DPGvkGbHnAh/t3bfQEqYd+/7z9epXj3G+Beyc5N607X8ksG2STWnH9FdnGP+Y2RyfY6Y7P/4GeFCSu1fV5VV17sgy7ghsU1W/qKrZPqd2cVV9qNpzjsfQ1tfWAFV1clX9uJqv0u4iPXpk3Fto17FbquoUWvJy/1nOH9pxcv24sutpx9pEw68HNp3mucD16lxbVZvPcNscAJxYVTf2/k37//HLf0/WzOa0D58zNdG6Z5L5/xdtPT8kyd2B19OSwHtMUPcNtA8uH+n9p9AamM7o8zgeOJR26/utSb6W5APjngG8oS/PlNbXJHAbWlP0mIt72YSS/E6Sz4w9nEm7777VLOY3+qzATdy+Y83ET0e6bwau6ieMsX7GTe/Ske6LaSf2rWgnqqf3Bz+vS3Id8CjaSWeiccfbBri0qn4zbvrbznA5tpkgttsk2T/J/4zE9iCmWMdJHpHky0lWJ7mediGdaptc2xPj0flv08fZiDvvD1Mt12y251NpF5KLk3w1yR9MUfeycRfJsRh37DFePrJ+/oXWKjNmqm03lenGmyymmZpyu09gH9oJ6U56Ivc+2m2hK/sD1feaqC4z219vi6uf9K9h5su2jNaqN5tlg5YAvr+qRl982RH4j5Ftez7tVtfWVfUL2ifz5/RnNp/F7UnWTHwY2DrJE2dQ936054neOJJkTrnvVdXptGNgjyQP6NM4qY/7DlqL5+eTXJDkkIlm2j+0nElL+B5DS/q+CTySNUsC1+R8O+n5sZ83nkk7x1ye5OS+rACvorWqfDvJuUn+Yk1jraqbeuemAEn2TnJakmt6PPtwx3Pc1XXHZ7Bne20ZcyMw/ji6F7cnK+OH3wu4cZKEHhbHufZOktyDdkdq9MPcWDI4fvlnk8iNupbZJZATrXsmmn9VfYGWuP077Y7iRb3eHV6yS/Ji2rOB+/YPnfQPGodU1UOq6mDa4ykfBH6f1vL7WNodldH9+560FtYpra9J4E9oB/2YHXoZtMx5vMNpt7R2rqp70W4xrYu3o37OSJbePyktW8tpjrZk7ED7dHQV7WL10f6JaOxvk6oafVB4soMa2vrZftzLAzvQbjfNxOUTxAZAkh1pt6BeTGvm3px2u3hsHU8U17/SLjbbV9VmtB12qm2yRZJNxs3/J7R1M/ZpfnTYTJdr1J3irKozqmoF7aL5n7SW4MlsO+7T9ViMlwK/pD1LOrbt7lVVD5xq3tPFNsPxJotppibd7pOYNAkEqKr3VtXvAbvQHrEYe8Nt/HLMZH+9La7e6rQlM1+21bTbXrNZNmjPX70uyVNHyi6lPQs2emzerW5/ePwY2t2LPYGbqupbM4yRnsy9kfYc03TnrPNpt6o/m2SsRWkm+94xtJa059JaUX7R531DVb2yqn6b9gD53yTZc5J5f5V26/d3aa0RX6Xd/t8N+NpkizfN8szGlOfHandj/pT2ofkHtPMVVXVFVf1lVW1De1niA5n4K2pmFWt/A/Tfabf4tu7nxFNYN9ed8c4FHjLuOH9ILx8b/tCRYQ8dGTaRxXCuncif0T4IfmVkmtfSzmGzWf6pnE07b83UROv+p1V19USVq+r9VbVzVW1N2382pF1LAegfUg6hPds94TcwJHkw7bGMI2iPXJ3VE/4zaPvFmP9Du109pfU1CTyOdiJelvYa+Otpz8tAa3m7d5LNRurfk/ZMyI39E+CMvh9nBv4XuFuSfZNsRHtIdtLXv2foOUl26Z9q3kQ7Kf+atnxPTPL4JBskuVva1+FsN8Ppjn3if1WSjfqbSE+kNRvPxAnAa5Js0ef5kpFhm9AO6tUASQ6ktQSO+Smw3bim6HsC11TVL5LsBvz5DGJ4Y9or/Y8GngD8W183JwBvTXLPnpD+DbfvD7Nxhzj7vJ6dZLNqt1p/RrutNJnfAl7a1+/TaQfZKdVueX4eeGeSeyW5S5L7JpnNIwk/pT06MNtjcsKYZjH+CX387ZJsQTsBTSjtTbS7VtX5kwz//d4CPPbywC+4fX3+lPbm3JiZ7K/7JHlU315vpj3sPqMW1b7ffBJ4Q5J7JNmFditpOufSXtR5f5KxN+s+SNv/duzLuSzJipF5fasv5zuZXSvgmI/Snlfba7qKVXUc7UPuF5Lcd4b73sdoF9Dn0J75oi/HE5LcrycX19NaNyfb/79Ka504ryeuX6E9v3ZhtUd2JjJ+m6+NSc+PSbZOsqInNr+ktc78pi/j00fOodfSzmMTLePqXj7TeDemXQtWA7cm2Zv2AWKNjC0TLSm4S1++jfrgr9C2zUvTvg7lxb38S/3/sbQEftu0rxx5Je2Z2qms7+faiRwAHDtBC+extHxhi379/0tGlr+vs7v13o37up0sWT+FcY8t9djvRkvwN+rjj52njwUO6tf0zWk5wtFMoI/3oDQ70JK49/RElrSvxXsb8Kc1yddc9bjfB7y030W5EBg7Rz6W9tLMmMcCn51kOW+zviaBb6Hdfjib9obhd3oZVfUDWpJ4QdptgW2Av6UlGTfQPgF+Yl0EUVXX014N/zDt09DPGdd0uwY+SttJrqCd+F/a53Up7VbPa2knlktprSgz2kb9xPxEYG/aJ7oPAPv39TUTb6Q1/V9Iu6jcdjGrqvNoF7hv0Q7uB9PeSh3zJdrF84okY6/tvxB4U5IbaEn8dJ/6rqCdpH8CfJz2QOxY7C+hrfsLaA/n/ivtLa7ZmijO5wIXpT1G8AJai85kTqd9FcFVtGe5njbyiW9/2oXhvL4cJ3LHW/nTGfsC76uTTPVc4mximomxlwy+RzvOPjlF3X2ZOsG8V5/etbR96Wpu/76tI2nPA12X5D9nuL/+K+32yTW0l7+eM4vlgtZyvSlt3zqa25+vmVJVfY92YfxQv7i/h9aq/fm+P58GPGLcaMfSjouPAST5YJIPznB+v6YdI1vOsP4xtA+QX0qynGn2vX5u+Q4tAfr6yKR2Br5AS5q+RXto/cuTzPabtBfnxlr9zqMl+ZO1AkJbb09Lcm2S905Rb1rTnB/vQktWfkLbVx7L7Q0Bvw+cnuRG2jZ82UQX2H6r963Af/d9dPdp4rmBdu4+gbbO/5zbb7OviefSHh86nPZc4c3c3pr5K9pXwOxPu733F7SX7MYeCfgX2lui36e1Kp3cyyaz3p1r074379F3nsxtw7eltUQfO8HgQ2kvK11M+7Dyjqoa/QqrH9LW57a0c93N3LG1c9SngQdk5Pv7aNfDm7m99e1m2mMR9Pn8I+15yEt6DIeOxH1ubv/O47vR1ueNtJe3vkV7L2DMW2jfaHBGbv9OxPHnkAOBc6rqrN7/Sdp2XN3HPaLP9z60uzH/Ocly3iZ3Tqql+dVbgT5WVTNt9Zx3aV8S+vyqetRCxzJmvmNK+3Lj91V7yH2u53U07e301831vNaFJPvTvqZjvdk/RiU5ivZy0qJYn5obi+Fcu9CSHAzsUlUvX+hY1lSSdwI/rqoPTFd3vfyJJEnrpa/QPvFqRNqjHS+ktWaud3pr4VNoz/NJmkJVHbHQMaytqnrlTOuur7eDJa2F3PFnlkb/Jr3lMp2q+sea3dfbrHNJHj3Zsi1QPI+n3Yr5Ke1Wz3olyZtptwjfUVUXLnQ8ktYv3g6WJEkaIFsCJUmSBsgkUJIkaYCW5IshW221VS1fvnyhw5AkSZrWWWeddVVVre2PUczakkwCly9fzplnnrnQYUiSJE0ryUx+1nKd83awJEnSAJkESpIkDZBJoCRJ0gCZBEqSJA2QSaAkSdIAmQRKkiQNkEmgJEnSAJkESpIkDZBJoCRJ0gCZBEqSJA2QSaAkSdIALcnfDpZmavkhJy90CHPqosP2XegQJEnrKVsCJUmSBsgkUJIkaYBMAiVJkgbIJFCSJGmATAIlSZIGyCRQkiRpgEwCJUmSBsgkUJIkaYBMAiVJkgbIJFCSJGmA5iwJTHJUkiuTnDOu/CVJfpDk3CT/OFL+miQrk/wwyeNHyvfqZSuTHDJX8UqSJA3JXP528NHA+4BjxwqS/BGwAnhoVf0yyW/18l2A/YAHAtsAX0jyO3209wN/CqwCzkhyUlWdN4dxS5IkLXlzlgRW1deSLB9X/NfAYVX1y17nyl6+Aji+l1+YZCWwWx+2sqouAEhyfK9rEihJkrQW5vuZwN8BHp3k9CRfTfL7vXxb4NKReqt62WTlkiRJWgtzeTt4svltCewO/D5wQpLfXhcTTnIwcDDADjvssC4mKUmStGTNd0vgKuCT1Xwb+A2wFXAZsP1Ive162WTld1JVR1TVrlW167Jly+YkeEmSpKVivpPA/wT+CKC/+LExcBVwErBfkrsm2QnYGfg2cAawc5KdkmxMe3nkpHmOWZIkacmZs9vBSY4D9gC2SrIKOBQ4Cjiqf23Mr4ADqqqAc5OcQHvh41bgRVX16z6dFwOfAzYAjqqqc+cqZkmSpKGYy7eDnzXJoOdMUv+twFsnKD8FOGUdhiZJkjR4/mKIJEnSAJkESpIkDZBJoCRJ0gCZBEqSJA2QSaAkSdIAmQRKkiQNkEmgJEnSAJkESpIkDZBJoCRJ0gCZBEqSJA2QSaAkSdIAmQRKkiQNkEmgJEnSAJkESpIkDZBJoCRJ0gCZBEqSJA2QSaAkSdIAmQRKkiQNkEmgJEnSAJkESpIkDZBJoCRJ0gCZBEqSJA3QnCWBSY5KcmWScyYY9soklWSr3p8k702yMsnZSR4+UveAJD/qfwfMVbySJElDMpctgUcDe40vTLI98DjgkpHivYGd+9/BwOG97pbAocAjgN2AQ5NsMYcxS5IkDcKcJYFV9TXgmgkGvQt4FVAjZSuAY6s5Ddg8yX2AxwOnVtU1VXUtcCoTJJaSJEmanXl9JjDJCuCyqvreuEHbApeO9K/qZZOVS5IkaS1sOF8zSnIP4LW0W8FzMf2DabeS2WGHHeZiFpIkSUvGfLYE3hfYCfhekouA7YDvJPn/gMuA7UfqbtfLJiu/k6o6oqp2rapdly1bNgfhS5IkLR3zlgRW1fer6reqanlVLafd2n14VV0BnATs398S3h24vqouBz4HPC7JFv2FkMf1MkmSJK2FufyKmOOAbwH3T7IqyUFTVD8FuABYCXwIeCFAVV0DvBk4o/+9qZdJkiRpLczZM4FV9axphi8f6S7gRZPUOwo4ap0GJ0mSNHD+YogkSdIAmQRKkiQNkEmgJEnSAJkESpIkDZBJoCRJ0gCZBEqSJA2QSaAkSdIAmQRKkiQNkEmgJEnSAJkESpIkDZBJoCRJ0gCZBEqSJA2QSaAkSdIAmQRKkiQNkEmgJEnSAJkESpIkDZBJoCRJ0gBtuNABSNKaWn7IyQsdwpy76LB9FzoESUuULYGSJEkDZBIoSZI0QCaBkiRJAzRnSWCSo5JcmeSckbJ3JPlBkrOT/EeSzUeGvSbJyiQ/TPL4kfK9etnKJIfMVbySJElDMpctgUcDe40rOxV4UFU9BPhf4DUASXYB9gMe2Mf5QJINkmwAvB/YG9gFeFavK0mSpLUwZ0lgVX0NuGZc2eer6tbeexqwXe9eARxfVb+sqguBlcBu/W9lVV1QVb8Cju91JUmStBYW8pnAvwA+27u3BS4dGbaql01WLkmSpLWwIElgkr8HbgU+vg6neXCSM5OcuXr16nU1WUmSpCVp3pPAJM8DngA8u6qqF18GbD9SbbteNln5nVTVEVW1a1XtumzZsnUetyRJ0lIyr0lgkr2AVwFPqqqbRgadBOyX5K5JdgJ2Br4NnAHsnGSnJBvTXh45aT5jliRJWorm7GfjkhwH7AFslWQVcCjtbeC7AqcmATitql5QVecmOQE4j3ab+EVV9es+nRcDnwM2AI6qqnPnKmZJkqShmLMksKqeNUHxkVPUfyvw1gnKTwFOWYehSZIkDZ6/GCJJkjRAJoGSJEkDZBIoSZI0QCaBkiRJA2QSKEmSNEBz9nawJEla+pYfcvJChzCnLjps34UOYc7YEihJkjRAJoGSJEkDZBIoSZI0QCaBkiRJA2QSKEmSNEAmgZIkSQNkEihJkjRAJoGSJEkDZBIoSZI0QCaBkiRJA2QSKEmSNEAmgZIkSQNkEihJkjRAJoGSJEkDZBIoSZI0QCaBkiRJAzRnSWCSo5JcmeSckbItk5ya5Ef9/xa9PEnem2RlkrOTPHxknAN6/R8lOWCu4pUkSRqSuWwJPBrYa1zZIcAXq2pn4Iu9H2BvYOf+dzBwOLSkETgUeASwG3DoWOIoSZKkNTdnSWBVfQ24ZlzxCuCY3n0M8OSR8mOrOQ3YPMl9gMcDp1bVNVV1LXAqd04sJUmSNEvz/Uzg1lV1ee++Ati6d28LXDpSb1Uvm6xckiRJa2HBXgypqgJqXU0vycFJzkxy5urVq9fVZCVJkpakaZPAJI9Msknvfk6Sf06y4xrO76f9Ni/9/5W9/DJg+5F62/WyycrvpKqOqKpdq2rXZcuWrWF4kiRJwzCTlsDDgZuSPBR4JfBj4Ng1nN9JwNgbvgcAnxop37+/Jbw7cH2/bfw54HFJtugvhDyul0mSJGktbDiDOrdWVSVZAbyvqo5MctB0IyU5DtgD2CrJKtpbvocBJ/TxLwae0aufAuwDrARuAg4EqKprkrwZOKPXe1NVjX/ZRJIkSbM0kyTwhiSvAZ4DPCbJXYCNphupqp41yaA9J6hbwIsmmc5RwFEziFOSJEkzNJPbwc8EfgkcVFVX0J7Le8ecRiVJkqQ5NZOWwKcDH+nf00dVXcKaPxMoSZKk9cBMWgK3Bs5IckKSvZJkroOSJEnS3Jo2Cayq19F+zu1I4HnAj5K8Lcl95zg2SZIkzZEZfVl0f3Hjiv53K7AFcGKSf5zD2CRJkjRHpn0mMMnLgP2Bq4APA39XVbf0t4R/BLxqbkOUJEnSujaTF0O2BJ5SVRePFlbVb5I8YW7CkiRJ0lyayTOBhwLbJzkQIMmyJDv1YefPcXySJEmaAzP57eBDgVcDr+lFGwEfm8ugJEmSNLdm8mLInwFPAn4OUFU/Ae45l0FJkiRpbs0kCfxVfzu4AJJsMrchSZIkaa7NJAk8Icm/AJsn+UvgC7S3hCVJkrRIzeTt4HcCfwL8DLg/8Hrga3MZlCRJkubWTJLAI6vqL4BTAZJsCpwC7DmXgUmSJGnuzCQJvCzJB6rqhUm2AE4GPjTHcUlaB5YfcvJChyBJWk/N5HsC/wG4MckHgc8D76yqj8x5ZJIkSZozk7YEJnnKSO/pwD8A3wYqyVOq6pNzHZwkSZLmxlS3g584rv+7tC+KfiLt62JMAiVJkhapSZPAqjpwPgORJEnS/JnJ9wRKkiRpiTEJlCRJGqAFSQKTvCLJuUnOSXJckrsl2SnJ6UlWJvlEko173bv2/pV9+PKFiFmSJGkpmfZ7ApPcFXgqsHy0flW9aU1mmGRb4KXALlV1c5ITgP2AfYB3VdXx/etoDgIO7/+vrar7JdkPeDvwzDWZtyRJkpqZtAR+ClgB3Ar8fORvbWwI3D3JhsA9gMuBPwZO7MOPAZ7cu1f0fvrwPZNkLecvSZI0aDP5xZDtqmqvdTXDqrosyT8BlwA3076A+izguqq6tVdbBWzbu7cFLu3j3prkeuDewFXrKiZJkqShmUlL4K4Mh1QAAA60SURBVDeTPHhdzbD/9NwKYCdgG2ATYK2TzCQHJzkzyZmrV69e28lJkiQtaTNJAh8FnJXkh0nOTvL9JGevxTz/BLiwqlZX1S20L51+JLB5vz0MsB1wWe++DNgeoA/fDLh6/ESr6oiq2rWqdl22bNlahCdJkrT0zeR28N7reJ6XALsnuQftdvCewJnAl4GnAccDB9CeRQQ4qfd/qw//UlXVOo5JkiRpUKZNAqvq4nU5w6o6PcmJwHdoL5t8FzgCOBk4PslbetmRfZQjgY8mWQlcQ3uTWJIkSWthJi2B61xVHQocOq74AmC3Cer+Anj6fMQlSZI0FP5iiCRJ0gCZBEqSJA2QSaAkSdIAmQRKkiQNkEmgJEnSAJkESpIkDZBJoCRJ0gCZBEqSJA2QSaAkSdIAmQRKkiQN0IL8bNxSsfyQkxc6hDl10WH7LnQIkiRpjtgSKEmSNEAmgZIkSQNkEihJkjRAJoGSJEkDZBIoSZI0QCaBkiRJA2QSKEmSNEAmgZIkSQNkEihJkjRAJoGSJEkDZBIoSZI0QAuSBCbZPMmJSX6Q5Pwkf5BkyySnJvlR/79Fr5sk702yMsnZSR6+EDFLkiQtJQvVEvge4L+q6gHAQ4HzgUOAL1bVzsAXez/A3sDO/e9g4PD5D1eSJGlpmfckMMlmwGOAIwGq6ldVdR2wAjimVzsGeHLvXgEcW81pwOZJ7jPPYUuSJC0pC9ESuBOwGvhIku8m+XCSTYCtq+ryXucKYOvevS1w6cj4q3rZHSQ5OMmZSc5cvXr1HIYvSZK0+C1EErgh8HDg8Kr6XeDn3H7rF4CqKqBmM9GqOqKqdq2qXZctW7bOgpUkSVqKFiIJXAWsqqrTe/+JtKTwp2O3efv/K/vwy4DtR8bfrpdJkiRpDc17ElhVVwCXJrl/L9oTOA84CTiglx0AfKp3nwTs398S3h24fuS2sSRJktbAhgs035cAH0+yMXABcCAtIT0hyUHAxcAzet1TgH2AlcBNva4kSZLWwoIkgVX1P8CuEwzac4K6BbxozoOSJEkaEH8xRJIkaYBMAiVJkgbIJFCSJGmATAIlSZIGyCRQkiRpgEwCJUmSBsgkUJIkaYBMAiVJkgbIJFCSJGmATAIlSZIGyCRQkiRpgEwCJUmSBsgkUJIkaYBMAiVJkgbIJFCSJGmATAIlSZIGyCRQkiRpgEwCJUmSBsgkUJIkaYBMAiVJkgbIJFCSJGmAFiwJTLJBku8m+Uzv3ynJ6UlWJvlEko17+V17/8o+fPlCxSxJkrRULGRL4MuA80f63w68q6ruB1wLHNTLDwKu7eXv6vUkSZK0FhYkCUyyHbAv8OHeH+CPgRN7lWOAJ/fuFb2fPnzPXl+SJElraKFaAt8NvAr4Te+/N3BdVd3a+1cB2/bubYFLAfrw63t9SZIkraF5TwKTPAG4sqrOWsfTPTjJmUnOXL169bqctCRJ0pKzEC2BjwSelOQi4HjabeD3AJsn2bDX2Q64rHdfBmwP0IdvBlw9fqJVdURV7VpVuy5btmxul0CSJGmRm/cksKpeU1XbVdVyYD/gS1X1bODLwNN6tQOAT/Xuk3o/ffiXqqrmMWRJkqQlZ8Ppq8ybVwPHJ3kL8F3gyF5+JPDRJCuBa2iJoyRJ673lh5y80CFIk1rQJLCqvgJ8pXdfAOw2QZ1fAE+f18AkSZKWOH8xRJIkaYBMAiVJkgbIJFCSJGmATAIlSZIGyCRQkiRpgEwCJUmSBsgkUJIkaYBMAiVJkgbIJFCSJGmATAIlSZIGyCRQkiRpgBb0t4MlSVNbfsjJCx3CnLrosH0XOgRpsGwJlCRJGiCTQEmSpAEyCZQkSRogk0BJkqQBMgmUJEkaIJNASZKkATIJlCRJGiCTQEmSpAEyCZQkSRogk0BJkqQBmvckMMn2Sb6c5Lwk5yZ5WS/fMsmpSX7U/2/Ry5PkvUlWJjk7ycPnO2ZJkqSlZiFaAm8FXllVuwC7Ay9KsgtwCPDFqtoZ+GLvB9gb2Ln/HQwcPv8hS5IkLS3zngRW1eVV9Z3efQNwPrAtsAI4plc7Bnhy714BHFvNacDmSe4zz2FLkiQtKQv6TGCS5cDvAqcDW1fV5X3QFcDWvXtb4NKR0Vb1svHTOjjJmUnOXL169ZzFLEmStBQsWBKYZFPg34GXV9XPRodVVQE1m+lV1RFVtWtV7bps2bJ1GKkkSdLSsyBJYJKNaAngx6vqk734p2O3efv/K3v5ZcD2I6Nv18skSZK0hhbi7eAARwLnV9U/jww6CTigdx8AfGqkfP/+lvDuwPUjt40lSZK0BjZcgHk+Engu8P0k/9PLXgscBpyQ5CDgYuAZfdgpwD7ASuAm4MD5DVeSJGnpmfcksKq+AWSSwXtOUL+AF81pUJIkSQPjL4ZIkiQNkEmgJEnSAJkESpIkDZBJoCRJ0gCZBEqSJA2QSaAkSdIAmQRKkiQNkEmgJEnSAJkESpIkDZBJoCRJ0gCZBEqSJA2QSaAkSdIAmQRKkiQN0IYLHYDWX8sPOXmhQ5AkSXPElkBJkqQBMgmUJEkaIJNASZKkATIJlCRJGiCTQEmSpAEyCZQkSRogk0BJkqQBWjRJYJK9kvwwycokhyx0PJIkSYvZokgCk2wAvB/YG9gFeFaSXRY2KkmSpMVrUSSBwG7Ayqq6oKp+BRwPrFjgmCRJkhatxfKzcdsCl470rwIesUCxSJLWEX+eUlo4iyUJnFaSg4GDe++NSX44x7PcCrhqjuehueU2XPzchoub22/xW/LbMG+fl9nsOC9zGWexJIGXAduP9G/Xy25TVUcAR8xXQEnOrKpd52t+Wvfchouf23Bxc/stfm7DxW2xPBN4BrBzkp2SbAzsB5y0wDFJkiQtWouiJbCqbk3yYuBzwAbAUVV17gKHJUmStGgtiiQQoKpOAU5Z6DhGzNutZ80Zt+Hi5zZc3Nx+i5/bcBFLVS10DJIkSZpni+WZQEmSJK1DJoFrwJ+wW/8l2T7Jl5Ocl+TcJC/r5VsmOTXJj/r/LXp5kry3b9Ozkzx8YZdAY5JskOS7ST7T+3dKcnrfVp/oL4uR5K69f2Ufvnwh41aTZPMkJyb5QZLzk/yBx+HikeQV/Rx6TpLjktzNY3DpMAmcJX/CbtG4FXhlVe0C7A68qG+nQ4AvVtXOwBd7P7TtuXP/Oxg4fP5D1iReBpw/0v924F1VdT/gWuCgXn4QcG0vf1evp4X3HuC/quoBwENp29LjcBFIsi3wUmDXqnoQ7cXM/fAYXDJMAmfPn7BbBKrq8qr6Tu++gXbh2Za2rY7p1Y4Bnty7VwDHVnMasHmS+8xz2BonyXbAvsCHe3+APwZO7FXGb8OxbXsisGevrwWSZDPgMcCRAFX1q6q6Do/DxWRD4O5JNgTuAVyOx+CSYRI4exP9hN22CxSLZqDfkvhd4HRg66q6vA+6Ati6d7td10/vBl4F/Kb33xu4rqpu7f2j2+m2bdiHX9/ra+HsBKwGPtJv6X84ySZ4HC4KVXUZ8E/AJbTk73rgLDwGlwyTQC1pSTYF/h14eVX9bHRYtVfjfT1+PZXkCcCVVXXWQseiNbYh8HDg8Kr6XeDn3H7rF/A4XJ/1ZzVX0JL5bYBNgL0WNCitUyaBszftT9hp/ZBkI1oC+PGq+mQv/unY7aX+/8pe7nZd/zwSeFKSi2iPXfwx7fmyzfutKbjjdrptG/bhmwFXz2fAupNVwKqqOr33n0hLCj0OF4c/AS6sqtVVdQvwSdpx6TG4RJgEzp4/YbcI9OdQjgTOr6p/Hhl0EnBA7z4A+NRI+f797cTdgetHbldpAVTVa6pqu6paTjvOvlRVzwa+DDytVxu/Dce27dN6fVuYFlBVXQFcmuT+vWhP4Dw8DheLS4Ddk9yjn1PHtp/H4BLhl0WvgST70J5VGvsJu7cucEgaJ8mjgK8D3+f258leS3su8ARgB+Bi4BlVdU0/wb2PdqvjJuDAqjpz3gPXhJLsAfxtVT0hyW/TWga3BL4LPKeqfpnkbsBHac9/XgPsV1UXLFTMapI8jPZiz8bABcCBtAYIj8NFIMkbgWfSvnHhu8Dzac/+eQwuASaBkiRJA+TtYEmSpAEyCZQkSRogk0BJkqQBMgmUJEkaIJNASZKkATIJlLTkJblxmuGbJ3nhPMTxpiR/Mk2dPZL84VzHIkkmgZIEmwNzngRW1eur6gvTVNsDMAmUNOdMAiUtGkmWJzk/yYeSnJvk80nuPkG9nZJ8K8n3k7xlpHzTJF9M8p0+bEUfdBhw3yT/k+QdU9QbP58bk7yrx/LFJMt6+cOSnJbk7CT/0X+DlSRHJ3la774oyRtH5vGAJMuBFwCv6LE8OsnTk5yT5HtJvrYu16ekYTMJlLTY7Ay8v6oeCFwHPHWCOu8BDq+qBwOjPzv2C+DPqurhwB8B7+y/UnEI8OOqelhV/d0U9cbbBDizx/JV4NBefizw6qp6CO1Xaw6dYFyAq/o8Dqf9IspFwAeBd/VYvg68Hnh8VT0UeNK0a0eSZsgkUNJic2FV/U/vPgtYPkGdRwLH9e6PjpQHeFuSs4Ev0H7+ausJxp9pvd8An+jdHwMelWQzYPOq+movPwZ4zCTL8slplgPgv4Gjk/wl7acqJWmd2HChA5CkWfrlSPevgTvdDu4m+k3MZwPLgN+rqluSXATcbS3qzWSeUxlbll8zyfm4ql6Q5BHAvsBZSX6vqq6e5Xwk6U5sCZS0FP03sF/vfvZI+WbAlT2x+yNgx15+A3DPGdQb7y7A03r3nwPfqKrrgWuTPLqXP5d2q3im7hBLkvtW1elV9XpgNbD9LKYlSZOyJVDSUvQy4F+TvBr41Ej5x4FPJ/k+cCbwA4CqujrJfyc5B/gs8PaJ6k3g58BuSV4HXAk8s5cfAHwwyT2AC4ADZxH7p4ET+8soL6G9JLIz7Rb1F4HvzWJakjSpVM327oUkCdrbwVW16ULHIUlrwtvBkiRJA2RLoCRJ0gDZEihJkjRAJoGSJEkDZBIoSZI0QCaBkiRJA2QSKEmSNEAmgZIkSQP0/wORGbRLsZ6DeAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_data_points_hist(d_real)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we generate a dictionnary with cumulative probability based on frequency of delays, for each keys in our reference dictionnary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 251,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cumul_distri_probas_dict(dico):\n",
    "    list_tot_points = []\n",
    "    for key in dico:\n",
    "        distrib = dico[key]\n",
    "\n",
    "        # get total number of elements \n",
    "        N = np.sum(distrib)\n",
    "\n",
    "        # make cumulative distribution probabilities\n",
    "        cdf_distrib = np.empty((len(distrib)), dtype=float)\n",
    "        save_x = 0\n",
    "        for x in range(len(distrib)):\n",
    "            cdf_distrib[x] = float(distrib[x])/float(N) + float(save_x)/float(N)\n",
    "            save_x += distrib[x]\n",
    "\n",
    "        dico[key] = cdf_distrib\n",
    "    return dico"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Construct recovery tables \n",
    "\n",
    "First approach is to simple sum up similar distribution to get a new distribution we can use. For that, we need to have transport type (`route_desc`), `time` (rounded to hour) and `stop_id` which are valid. We then make all combination of these tree parameters and get the associate distributions.\n",
    "\n",
    "First we need to reformat stoptimes table in order to get time rounded to the hour, correct stop_id format and type, generate `key` column using `trip_id` and `stop_id`, and droping unnecessary columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0%, 3.84%, 7.68%, 11.52%, 15.36%, 19.2%, 23.04%, 26.88%, 30.72%, 34.55%, 38.39%, 42.23%, 46.07%, 49.91%, 53.75%, 57.59%, 61.43%, 65.27%, 69.11%, 72.95%, 76.79%, 80.63%, 84.47%, 88.31%, 92.15%, 95.98%, 99.82%, "
     ]
    }
   ],
   "source": [
    "with open(\"../data/stop_times_df_cyril.pkl\", \"rb\") as input_file:\n",
    "    stoptimes = pickle.load(input_file)\n",
    "\n",
    "# Use stop_id_general as stop_id  \n",
    "stoptimes['stop_id'] = stoptimes['stop_id_general'].apply(str)\n",
    "\n",
    "# Set same stoptimes index as distribution dict \n",
    "stoptimes['key'] = stoptimes['trip_id'] + '__' + stoptimes['stop_id']\n",
    "stoptimes = stoptimes.set_index('key')\n",
    "\n",
    "stoptimes = stoptimes[['trip_id','stop_id', 'route_desc', 'arrival_time', 'departure_time']]\n",
    "\n",
    "list_hours = []\n",
    "size_stop_times = stoptimes.shape[0]\n",
    "for x in range(size_stop_times):\n",
    "    if (x % 10000) == 0 :\n",
    "        print('{}%'.format(round(100*x/size_stop_times,2)), end = ', ')\n",
    "    \n",
    "    arr_time_hour = pd.to_datetime(stoptimes.iloc[x,:]['arrival_time']).hour\n",
    "    if math.isnan(arr_time_hour): # if arrival is NaT, use departure time\n",
    "        arr_time_hour = pd.to_datetime(stoptimes.iloc[x,:]['departure_time']).hour\n",
    "    list_hours.append(int(arr_time_hour))\n",
    "    \n",
    "stoptimes['hour'] = list_hours\n",
    "stoptimes = stoptimes.drop(columns=['trip_id', 'arrival_time', 'departure_time'])\n",
    "\n",
    "# Write this pickle to avoid re-running this above code all the time\n",
    "with gzip.open(\"../data/stop_times_wHour.pkl\", \"wb\") as output_file:\n",
    "    pickle.dump(stoptimes, output_file)         \n",
    "              "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(16895, 32)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>...</th>\n",
       "      <th>22</th>\n",
       "      <th>23</th>\n",
       "      <th>24</th>\n",
       "      <th>25</th>\n",
       "      <th>26</th>\n",
       "      <th>27</th>\n",
       "      <th>28</th>\n",
       "      <th>29</th>\n",
       "      <th>30</th>\n",
       "      <th>31</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>stop_id</th>\n",
       "      <th>hour</th>\n",
       "      <th>route_desc</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">8500926</th>\n",
       "      <th>7</th>\n",
       "      <th>Bus</th>\n",
       "      <td>49</td>\n",
       "      <td>681</td>\n",
       "      <td>108</td>\n",
       "      <td>56</td>\n",
       "      <td>38</td>\n",
       "      <td>26</td>\n",
       "      <td>18</td>\n",
       "      <td>16</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <th>Bus</th>\n",
       "      <td>63</td>\n",
       "      <td>700</td>\n",
       "      <td>115</td>\n",
       "      <td>61</td>\n",
       "      <td>22</td>\n",
       "      <td>10</td>\n",
       "      <td>9</td>\n",
       "      <td>5</td>\n",
       "      <td>3</td>\n",
       "      <td>7</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <th>Bus</th>\n",
       "      <td>2</td>\n",
       "      <td>407</td>\n",
       "      <td>210</td>\n",
       "      <td>73</td>\n",
       "      <td>31</td>\n",
       "      <td>14</td>\n",
       "      <td>5</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <th>Bus</th>\n",
       "      <td>0</td>\n",
       "      <td>204</td>\n",
       "      <td>121</td>\n",
       "      <td>99</td>\n",
       "      <td>47</td>\n",
       "      <td>17</td>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>4 rows × 32 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                         0    1    2   3   4   5   6   7   8   9   ...  22  \\\n",
       "stop_id hour route_desc                                            ...       \n",
       "8500926 7    Bus         49  681  108  56  38  26  18  16   4   2  ...   0   \n",
       "        8    Bus         63  700  115  61  22  10   9   5   3   7  ...   0   \n",
       "        9    Bus          2  407  210  73  31  14   5   2   0   1  ...   0   \n",
       "        10   Bus          0  204  121  99  47  17   5   4   2   2  ...   0   \n",
       "\n",
       "                         23  24  25  26  27  28  29  30  31  \n",
       "stop_id hour route_desc                                      \n",
       "8500926 7    Bus          0   0   0   0   0   0   0   0   0  \n",
       "        8    Bus          0   0   0   1   0   0   0   0   0  \n",
       "        9    Bus          0   0   0   0   0   0   0   0   2  \n",
       "        10   Bus          0   0   0   0   0   0   0   0   0  \n",
       "\n",
       "[4 rows x 32 columns]"
      ]
     },
     "execution_count": 270,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with gzip.open(\"../data/stop_times_wHour.pkl\", \"rb\") as input_file:\n",
    "    stoptimes = pickle.load(input_file)\n",
    "\n",
    "distrib_df = pd.DataFrame(d_all).transpose()\n",
    "stoptimes_df = pd.DataFrame(stoptimes)\n",
    "\n",
    "recovery_df = distrib_df.join(stoptimes_df)\n",
    "list_bins = [x for x in range(32)]\n",
    "\n",
    "recovery_df = recovery_df.groupby(['stop_id','hour', 'route_desc'])[list_bins].apply(lambda x : x.astype(float).sum())\n",
    "recovery_df = recovery_df.astype('int')\n",
    "\n",
    "print(recovery_df.shape)\n",
    "recovery_df.head(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 271,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_df_missing(df, max_bin = 10000):\n",
    "    tot_per_key = np.array(df.sum(axis=1)).astype('int')\n",
    "    binwidth = 100\n",
    "    n_keys_less_than_binwidth = np.sum(np.array(tot_per_key < binwidth))\n",
    "    perc_key_to_recover = round(100 * ( n_keys_less_than_binwidth / len(tot_per_key) ), 2)\n",
    "    plt.figure(figsize = (10,5))\n",
    "    plt.hist(tot_per_key, bins = range(min(tot_per_key), max_bin + binwidth, binwidth))\n",
    "    plt.title(\"Total number of data points per stop_id / hour key. N keys with less than {0} points: {1} ({2}%)\"\\\n",
    "             .format(binwidth, n_keys_less_than_binwidth, perc_key_to_recover))\n",
    "    plt.xlabel('n data points')\n",
    "    plt.ylabel('n keys')\n",
    "    return plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAFNCAYAAABv3TlzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3debgcVZn48e8r+x6WDAOBIYiog7siMoMLij8XEOMgKopsg8O478OiKKLo4L4LooyAoIAMCioOIgjuaEBkEYEAwRDABITIJoK8vz/OuaTT6e7bl9y+3ZV8P89zn1td66k6dareOudUd2QmkiRJapZHDDsBkiRJmjiDOEmSpAYyiJMkSWoggzhJkqQGMoiTJElqIIM4SZKkBlqhgriIyIh41Aik4/yIeN2Qtr1GRHw3IhZFxLf6mH/HiLhxKtK2LCJiz4j44bDT0XQR8Z6I+GqP6XMj4vkPc93DPO/3jYifDWPb/YiI4yLiiGGno5fxythErxXDPB/q9kf+mE+liDg6It437HQMWkT8d0S8fdjp6CUinhgRv+hn3pEI4iLirpa/ByPi3pbPe3ZZphHBxQjaHdgY2DAzXzGZKx7mjTIzT8rMF/Qz76jf0DuJiJn1IWTlQW4nMz+Smct0Y42IqyLi0ZOVpqap51dGxIFt42+MiB2HlKxl1l7GRuWhuB/DLvMR8aGIuCwiHoiID3SY/pqIuCEi7o6I70TEBi3TNoiIb9dpN0TEawaRxsx8fWZ+qJ95JzMAjojnRsSPa8XC3A7T/zUifh0Rd0bEpRHxzLbp0yPiG3X52yPipB7bmg7sDXy5fl41Ik6rD6c5XvmMiDdHxOyIuC8ijuswfc2I+FJE3FrT85OWaa+JiJvrtp7bMn6riPhFRKw0Ni4zLwXuiIhde6UHRiSIy8y1x/6APwK7tozrmiEruigmmodbAFdn5gODSJP6N+iAbBgiYitgpcy8eohpGIXj+mfgwIhYZ9gJ0UiYAxwIfL99QkQ8jhJU7EV5wL4H+FLLLF8E/lan7QkcVZdZXtwN/A/wX+0TajD7XeDjwDTgY8B3I2L9ltlOB24B/gn4B+ATPba1L3BWZt7bMu5nwGvrOsZzE3BETW8nxwAbAP9c/7+j7sfKwJHAU4E3A59vWeZzwDsy8+9t6zoJ+M9xU5SZI/UHzAWeX4dXAz5TD9xNdXg1YC3gXuBB4K76tymwHfBL4A7gZuALwKot607gUV22ez7wIeDnwJ3AD4GN6rQdgRt7pPMDwLeAE+uylwGPBg4BFgDzgBe0beu/gV8DfwHOADZomb498Iu6H78Ddmxb9sM1nfd22h/KCXR+Xf4K4KV1/OGUi8H99Zjt32HZNYDjgNuB31MK1o0t0w8Grq37+Xvg31q2+Vfg73Xdd9TxuwC/rfs5D/hAj7zfEbgReA9waz3Ge7ZMXw84AVgI3AAcCjyiTtsX+FlbXr8euKYehy8C0SOdO9f9uROYD7y7Sxr3rcf+C8Ai4A/ATm1pPJZy/s2nFPiV2pb9NHAbcESH9W8HzK7H60/Ap+r4P9Z9Gjvf/4XyEHZoPRYL6rFZr84/s85/AKXs3Nxtn9q2/wHgxJbPe9X13wa8l5bzvsvybwU+N9EyVqe/lHK+3lHn/eduZZdyjh7Rdt4cRLkQf71LvrWeHx+nXLzX65ZnwKqUYOwJLcv9A+UmO73HMdi3rvu7wGEt42+kpSy3LdO6P+sAP6Zc3AN4LHBOTctVwCvrfE+v58hKLevZDfhdr3Opw7YvAF5eh3eox3qX+nkn4JL2Ywj8pM53N+V8fFVLPryLcj7eDOzX4zidD7yu5fO/A1dSrj1nA1vU8UEpMwvqvlwGPL7fckv3Mn8c5brw/br8hcBWLct9lnLN+gtwEfCstnJyKqXM3Uk5b7fto3ydSNs1EPgI8I2Wz1tRrtPrUO51fwMe3TL968CRPcrvacApNV0XA09qOxbn03Zv6FGmlspLyjXl/pquu4Dv1vEH1Ty4k3Ke7jTe8WhL+/OBuW3jXgJc0Tbuauq9C3gB5Zq0Up/bOA94bZdpXctnh3mPAI5rG/fYeq6s22H+jYFf1uHVgXvq8O7AMV22MYNyj1+tV1pGoiauh/dSAponA0+iXJQOzcy7gRcDN+XiGrubKIX0HcBGlJvcTsAbJ7C91wD7US7UqwLvnsCyu1IK1/qUoOVsyk12BvBBavVti70pF61NgAcoF2wiYgblonIEJZJ/N/C/tRp4zF6UgrQO5Qb7kIhYhXLz+GHdj7cAJ0XEYzLzMMoF45R6zI7tsB+HUS4iWwEvBPZpm34t8CzKje9w4MSI2CQzr6QETb+s655W57+77us0SkD3hoh4WacDWP0jJf9m1G0fExGPqdM+X7f7SOA5db379VjXSyg3uicCrwRe2COdxwL/mZnrAI+nFPZunlGPw0aU43V6S/PHcZT8fBTwFMpF5nVty15HKdQf7rDuzwKfzcx1KXlwah3/7Pp/Wk33Lyk31X2B51KOydqU4LLVc4GtazoOmkh/tojYBjiKcr5tCmwIbDbOYjvTobahRccyVptfvwm8HZgOnEV54l61z+T+I6W8bEEpGx1FxCMi4iuUc+IFmbmILnmWmX8DTqY8pY95NXBuZi7sI03vA97e2jQ2nojYEDgX+HlmvhVYkxLAfYNyzPYAvhQR22TmbyjBdWs3gr0ogQV0P5faXUC5aUMpV9ex+Hx7Tp2+hMwcm/6kej6eUj//I6WMzgD2B77YVmvSbb9nUR7edqPk/08p5wN1/55NeTBej1KWb6vTxi23Pco8lON5OOW6PYcly+RvKPeeDSjH/1sRsXrL9JdSzo9pwJksXfb69TjKw/pYeq+lBm7174Fcsmb7d3WZbmZRKhXG0v2diFil172hy3o65mVmHkOpJfpYPZ671nW8GXh6zYsXUoIrIuKZEXFH30djadHh8+Pr8PaUgPH4iLgtIn4TEc/psa4n1PkHYTvK/fjw2px6WUS8vE5bCGwYEZsB/w+4otbSH0qp7FlKZs6nBMvd8gcYkebUHvYEPpiZC+pF83DKRaqjzLwoM3+VmQ9k5lxK4NQrQ9t9LTOvzlLVeiqlAPfrp5l5dpZmym9RLkRHZub9lII+MyJaLx5fz8zLa0D6PuCVtU38tZTq3rMy88HMPIfyNL1zy7LHZeYVdT/vb0vH9pSb+ZGZ+bfMPA/4HuXm049XAh/OzD9n5jxqcDkmM7+VmTfVtJ1CqenartvKMvP8zLyszn8p5cI8Xp68LzPvy8wLKAHB2LHZAzgkM++s+ftJepwPlGNwR2b+kVKz0Ss/7we2iYh1M/P2zLy4x7wLgM9k5v31GFwF7BIRG1Py6e2ZeXdmLqDUIOzRsuxNmfn5mnf3Lr1q7gceFREbZeZdmfmrHunYk1K7cl1m3kW5GOzR1px4eE3LZcDX6P88gPKU+L3M/Elm3kc5Tx/sNnNErEkJms/vsc5uZexVwPcz85x6Tn+CUiv8r32m9UFKrdd9XY4rwCqU828DSpeNe/rIs+OBV0fE2I1kL8rD2rgy8xJKAHZQn/uwKSVg+lZmHlrHvYRSO/G1es78FvhfYKw/6/HUILMGiy+k3Lih/3PpAhaXyWdTWgnGPncM4nq4n3LNvj8zz6LU1PS8CVWvB/47M6+s19CPAE+OiC3qOteh1HREnefmlu31W247+XZm/rpu8yRarhGZeWJm3laP+ycprUCt+/Kzep3+O+WceNIEtz1mbUqtfqtFlH1em1K702laNxdl5mm1HH2KUvOzPRO/N0wkL/9OOT7bRMQqmTm3BqNk5s/aAueJ+CWwaUS8ugai+1AeSNas0zejBPk/pgSdnwTOiIiNuqxvGqWmcBA2owSXiyhl+c2U4PKfM/NB4A2UWtJ3A/9BiWc+Dzyx9gk8OyIe37bOO2uauxr1IG5TlqxpuqGO6ygiHh0R34uIWyLiL5QLQbfM7KS1Tfweygnfrz+1DN8L3JqL27jHbiqt65vXMnwD5QazEaUm4RURccfYH/BMSo1dp2XbbQrMqydN6/pn9Lkfm3ZI20MiYu+IuKQlbY+nxzGOiGfUE3RhRCyiXKx75cntNbBt3f6mdZlVWPp86LVfE8nPl1Nu5jdExAUR8S895p2fWeq729K4RU3jzS3H58uUp94xvfIOyhPvo4E/1KfKl/SYt1P5WJlSy9dpez3LT5f1P7R8zZfbus/OTsAvasDXTbc8WWJf6vk7j/7P24WZ+ddx5nkUpZbi8Cy1bDBOnmXmhTWdO0bEY+s6zuwzTQDvp9Q+bzzunKWmeg3g6JZxWwDPaLse7Em5YUFpnts1ItaiPID9tCXA6fdc+iXw6JrGJ1Nq8javN8LtKE2n/botl+xv2+91dAvgsy37+GdKjcuMGmx8gdL0uSAijomIdetyEym3nXS9RkTEuyPiyigd1O+g1Ept1GPZ1R9mf8y7gHXbxq1LuYH3mtZNa5l9kNJMuCkTvzf0nZeZOYdSi/4BSh6dHBETudZ0lJm3UcrsOyn32BcBP6LsE5R769zMPLYGmydT9n+HLqu8nd4B8LK4lxL4HlGD5AsoweUL6r6cm5nbZ+ZzKF0RtqW0ApxAaVH5END+zQDrUJq+uxr1IO4mSuEe8091HJSD0O4oSh+lrbM0IbyHpatiH467WRz5U2uFpnefvS+btwz/EyXzb6WcgF/PzGktf2tl5pEt83fa9zE3US7ArXn7T5S+Cv24uUPaAKhPxV+hPGFsWJ+uLmfxMe6Urm9QbnqbZ+Z6lBtUrzxZv96QWrd/E+XY3M/S50O/+9VqqXRm5m8ycxbl5v0dujc9AcxoqZlpTeM84D5KP6+xvFs3M1ubPnrlHZl5TWa+uqbjo8Bp9Xh0Wq5T+XiAJR8o2vPyJvq3xLlQa9o27DH/zpRm0IdjiX2px3dzFufvPbSUQRYHMWN6HtfqSkpT7g9ampD6ybOx2q69gNP6CBYXJyrzD5SO1+/tY/avAP8HnNVSBuYBF7RdD9bOzDfU9c+nBGG70VZL2ONcak/jPZQ+X28DLq8B7i8oN85rM/PWfvd3GcyjNIu27ucamfmLmsbPZebTgG0ogel/1fH9ltt+zo+HRMSzKC8ivBJYv17rFjE595N2V9BSixcRj6TUal1d/1aOiK1b5n9SXaab1jL7CEoN0Vi/8mW5N7TqdA39RmY+k1KOk3LOLbPMvCAzn56ZG1DO8cdS+pMDXNohLb3y+lLK+TMIl3YYt1Ra6rXtC5T+wxtR+vPdQGm+f2LLfDMoXU56Nv+OehD3TeDQKK8Qb0R5qj2xTvsTpY15vZb516FUPd9Vn5rfMEnpuJrylLVL7VdwKKWQLYvXRsQ29cb4QcrN4e8sfrJ+YUSsFBGrR/k6lfH6Io0Zqzk4sFY/70jpr3dyn8ufChwSEevXbb6lZdpYMLEQICL2Y3HfBCh5slks2Y9pHeDPmfnXiNiO0idqPIdHefX7WZTmpG/VY3Mq8OGIWKcGlO9k8fkwEUuks25rz4hYrzZB/IUezYaUG8Zb6/F9BaWz8Fm1BuSHwCcjYt0o/a+2it59NJYQEa+NiOn1aXnsCexByjF/kNL3bcw3gXdExJYRsTaL+zu2Pj2/L8pr74+jBDCn0L/TgJdE6dOyKuU87XXNeDG9+8P1ciqlSXqnWsbeRQmuxr4r6RLgNbVMvIiJdZN4SGZ+k/Jw96OI2KrPPDsR+DdKIHfC0msd1+GUY99Pk9KbKRft70bEGpTmrkdHxF71fFslIp4eEf/csswJlIDjCZSAEeh5LnVyQd32WNPp+W2fO/kTS56Py+JoynXncQARsV4tW9T9fUY9L+6mvKTw4ATLbadrUy/rUB6IFlKCqPezdI1Y32q+rU4pPyvX6/rYV0qcRLnmP6sG2R8ETs/SbeRuSp5+MCLWiogdKDVTvZr0nxYRu9VawbdTytGvWPZ7Q6sl8j4iHhMRz4uI1Sj5M/bi4bhqmVudUiMe9dis2jL9KTW961K6WczLzLPr5G9THvz3qdeG3SlB68+7bO4s2q4dEbFaLO7ruGrdfsdgPSJWrvOuBIzdn8dqX39CeQHtkDrfDpQ+yWe3reZ1wMW1u8VtwBpR+h8/l9IfdcxzgPPGadkY+SDuCEp/sEspbyRdXMeNPeF+E7guShX8ppS25tdQqpq/wsRuWF1l6fz8RkpV53zKhWRZv6Pu65Sq1FsofRbeWrc1j1JI30O5gMyjPHX2lVf1KXpXyg31Vsqr6nvX49WPwylV7NdTbm6tT/a/p/Q5+CWlED+BJQvLeZQnxFsiYuzp/Y2UC9CdlCC8Vw0XlONxO+Wp8STg9S1pfwvl2F9HefvvG3R/1buXTuncC5gbpRn+9ZQmq24upLwscCulI/TutdofyssWq1LemLudEght0mklXbyI0un1LkrH9D0y895aW/Jh4Of1fN+esu9fp1w8rqdcPN/Str4LKB22zwU+kZl9fyFyZl4BvIlynG+u+9PxvI/Sl+OuLP0PJywzr6IESZ+nHNddKf3Wxpo931bHjTUnfufhbKdu63jKjfK8iJjJOHlWy+TFlAeYnwJExBXR5TssO2zveko+LVUL1mHesTeKb6S8tX4/pTlmD0qZuIVSw9H6EPltSu3Ht+t5MqbjudRl0xdQApefdPncyQcofX7uiIhXjrdvvWTmtyn7dXItg5dTrmFQgqevUPJm7E3pj9dp/ZbbTmW+l7MptaJX123+lfG7QvTyFUpg82pKrey9Ne1j5ez1lOvdAspxb30h742UZvYFlHveG+oy3ZxB6WN6e93GbrWpcVnvDa2OpfR/uyMivkM5H4+s672F8qB7CJRazXoOdvNsyvE4i1IzeC/l3jPmQBa3Um1CeaACIDP/THnB5N2UmtKDgVk9ao9PAHauD0hjrqrbnEHJ93uprQJRvvz8By3zHlqnH0y5Xt1bx1EfJGZRWiQWUfJ8ieMbpTLqbZT+xdQH7jdTzs+jWfL6vSdLdq3oKJbs2iMNT30yPDEz+611nHIRsS/lzcVnjjfvMNXg5HpglZyC7wSM8sW2G2XmgePO3EAR8T+Ul1IOHXfmIYiIaynNkT8adlo0PFG+SPhRmfna8eZdUUXER4AFmfmZYaelm4h4IvDlzBy3j+cofCmmpOabS/n6guVODYh3o3z9yMiJ8jUGSe+vxZEEZOZ7hp2G8WT5Joe+XtIZ9eZUSZMsIn4QS/7U3djfw764ZeapWb6Pa7kSER+iNO19vDaLjpSIOJ/yQteb2t46lLQCsDlVkiSpgayJkyRJaiCDOEmSpAZq9IsNG220Uc6cOXPYyZAkSRrXRRdddGtmLuuPBTyk0UHczJkzmT179rCTIUmSNK6IuGH8ufpnc6okSVIDGcRJkiQ1kEGcJElSAxnESZIkNZBBnCRJUgMZxEmSJDWQQZwkSVIDGcRJkiQ1kEGcJElSAxnESZIkNZBBnCRJUgM1+rdTtWKYefD3lxo398hdhpASSZJGhzVxkiRJDWQQJ0mS1EAGcZIkSQ1kECdJktRABnGSJEkNZBAnSZLUQAZxkiRJDWQQJ0mS1EAGcZIkSQ1kECdJktRABnGSJEkNZBAnSZLUQAZxkiRJDWQQJ0mS1EAGcZIkSQ1kECdJktRABnGSJEkNZBAnSZLUQAZxkiRJDWQQJ0mS1EAGcZIkSQ1kECdJktRABnGSJEkNZBAnSZLUQAZxkiRJDWQQJ0mS1EADDeIi4h0RcUVEXB4R34yI1SNiy4i4MCLmRMQpEbFqnXe1+nlOnT5zkGmTJElqsoEFcRExA3grsG1mPh5YCdgD+Cjw6cx8FHA7sH9dZH/g9jr+03U+SZIkdTDo5tSVgTUiYmVgTeBm4HnAaXX68cDL6vCs+pk6faeIiAGnT5IkqZEGFsRl5nzgE8AfKcHbIuAi4I7MfKDOdiMwow7PAObVZR+o8284qPRJkiQ12SCbU9en1K5tCWwKrAW8aBLWe0BEzI6I2QsXLlzW1UmSJDXSIJtTnw9cn5kLM/N+4HRgB2BabV4F2AyYX4fnA5sD1OnrAbe1rzQzj8nMbTNz2+nTpw8w+ZIkSaNrkEHcH4HtI2LN2rdtJ+D3wI+B3es8+wBn1OEz62fq9PMyMweYPkmSpMYaZJ+4CykvKFwMXFa3dQxwEPDOiJhD6fN2bF3kWGDDOv6dwMGDSpskSVLTrTz+LA9fZh4GHNY2+jpguw7z/hV4xSDTI0mStLzwFxskSZIayCBOkiSpgQziJEmSGsggTpIkqYEM4iRJkhrIIE6SJKmBDOIkSZIayCBOkiSpgQziJEmSGsggTpIkqYEM4iRJkhrIIE6SJKmBDOIkSZIayCBOkiSpgQziJEmSGsggTpIkqYEM4iRJkhrIIE6SJKmBDOIkSZIayCBOkiSpgQziJEmSGsggTpIkqYEM4iRJkhrIIE6SJKmBDOIkSZIayCBOkiSpgQziJEmSGsggTpIkqYEM4iRJkhrIIE6SJKmBDOIkSZIayCBOkiSpgQziJEmSGsggTpIkqYEM4iRJkhrIIE6SJKmBDOIkSZIayCBOkiSpgQziJEmSGsggTpIkqYEM4iRJkhrIIE6SJKmBDOIkSZIayCBOkiSpgQziJEmSGsggTpIkqYEM4iRJkhrIIE6SJKmBBhrERcS0iDgtIv4QEVdGxL9ExAYRcU5EXFP/r1/njYj4XETMiYhLI+Kpg0ybJElSkw26Ju6zwP9l5mOBJwFXAgcD52bm1sC59TPAi4Gt698BwFEDTpskSVJjDSyIi4j1gGcDxwJk5t8y8w5gFnB8ne144GV1eBZwQha/AqZFxCaDSp8kSVKTDbImbktgIfC1iPhtRHw1ItYCNs7Mm+s8twAb1+EZwLyW5W+s4yRJktRmkEHcysBTgaMy8ynA3SxuOgUgMxPIiaw0Ig6IiNkRMXvhwoWTllhJkqQmGWQQdyNwY2ZeWD+fRgnq/jTWTFr/L6jT5wObtyy/WR23hMw8JjO3zcxtp0+fPrDES5IkjbKBBXGZeQswLyIeU0ftBPweOBPYp47bBzijDp8J7F3fUt0eWNTS7CpJkqQWKw94/W8BToqIVYHrgP0ogeOpEbE/cAPwyjrvWcDOwBzgnjqvJEmSOhhoEJeZlwDbdpi0U4d5E3jTINMjSZK0vPAXGyRJkhrIIE6SJKmBDOIkSZIayCBOkiSpgQziJEmSGsggTpIkqYEM4iRJkhrIIE6SJKmBDOIkSZIayCBOkiSpgQb926kakpkHf3+pcXOP3GUIKZEkSYNgTZwkSVIDGcRJkiQ1kEGcJElSAxnESZIkNdC4QVxE7BARa9Xh10bEpyJii8EnTZIkSd30UxN3FHBPRDwJeBdwLXDCQFMlSZKknvoJ4h7IzARmAV/IzC8C6ww2WZIkSeqln++JuzMiDgFeCzw7Ih4BrDLYZEkaZX4PoSQNXz81ca8C7gP2z8xbgM2Ajw80VZIkSeqpn5q4VwBfy8zbATLzj9gnTpIkaaj6qYnbGPhNRJwaES+KiBh0oiRJktTbuEFcZh4KbA0cC+wLXBMRH4mIrQacNkmSJHXR15f91rdTb6l/DwDrA6dFxMcGmDZJkiR1MW6fuIh4G7A3cCvwVeC/MvP++pbqNcCBg02iJEmS2vXzYsMGwG6ZeUPryMx8MCJeMphkSZIkqZd++sQdBmweEfsBRMT0iNiyTrtywOmTJElSB/38duphwEHAIXXUKsCJg0yUJEmSeuvnxYZ/A14K3A2QmTfhz25JkiQNVT9B3N/q26kJEBFrDTZJkiRJGk8/QdypEfFlYFpE/AfwI8pbqpIkSRqSft5O/STwfOAvwGOA9wM/GWSiJEmS1Fs/QdyxmfnvwDkAEbE2cBaw0yATJkmSpO76aU6dHxFfAoiI9YEf4tupkiRJQ9XP98S9D7grIo6mBHCfzMyvDTxlkiRJ6qprc2pE7Nby8ULgfcCvgYyI3TLz9EEnTpIkSZ316hO3a9vn31K+6HdXyteNGMRJkiQNSdcgLjP3m8qESJIkqX/9vNggSZKkEWMQJ0mS1EAGcZIkSQ007pf9RsRqwMuBma3zZ+YHB5csSZIk9dLPLzacASwCLgLuG2xyJEmS1I9+grjNMvNFA0+JJEmS+tZPn7hfRMQTBp4SSZIk9a2fmrhnAvtGxPWU5tQAMjOfONCUSZIkqat+grgXDzwVkiRJmpBxg7jMvGEqEiJJkqT+Dfx74iJipYj4bUR8r37eMiIujIg5EXFKRKxax69WP8+p02cOOm2SJElN1U9z6rJ6G3AlsG79/FHg05l5ckQcDewPHFX/356Zj4qIPep8r5qC9GmEzDz4+8NOgiRJjTDQmriI2AzYBfhq/RzA84DT6izHAy+rw7PqZ+r0ner8kiRJajPo5tTPAAcCD9bPGwJ3ZOYD9fONwIw6PAOYB1CnL6rzS5Ikqc3AgriIeAmwIDMvmuT1HhARsyNi9sKFCydz1ZIkSY0xyJq4HYCXRsRc4GRKM+pngWkRMdYXbzNgfh2eD2wOUKevB9zWvtLMPCYzt83MbadPnz7A5EuSJI2ugQVxmXlIZm6WmTOBPYDzMnNP4MfA7nW2fSi/zQpwZv1MnX5eZuag0idJktRkA/+KkQ4OAt4ZEXMofd6OreOPBTas498JHDyEtEmSJDXCVHzFCJl5PnB+Hb4O2K7DPH8FXjEV6ZEkSWq6YdTESZIkaRkZxEmSJDWQQZwkSVIDGcRJkiQ10JS82KDlU7ffOZ175C5D2fZUbFeSpFFhTZwkSVIDWRMnTSFrECVJk8WaOEmSpAYyiJMkSWoggzhJkqQGMoiTJElqIIM4SZKkBjKIkyRJaiCDOEmSpAbye+IkSZoEfg+kppo1cZIkSQ1kECdJktRANqdKmhQ2JUnS1LImTpIkqYEM4iRJkhrIIE6SJKmB7BO3ArHPkiRJyw9r4iRJkhrIIE6SJKmBDOIkSZIayD5xy4FOfd0kSdLyzZo4SZKkBjKIkyRJaiCDOEmSpAayT9yI8DvcJEnSRFgTJ0mS1EAGcZIkSQ1kECdJktRABnGSJEkNZBAnSZLUQAZxkiRJDWQQJ0mS1EAGcZIkSQ3kl/1KkjQgfpG7BsmaOEmSpAayJk5aQVlDIEnNZk2cJElSAxnESZIkNSMRplUAAAwISURBVJBBnCRJUgMZxEmSJDWQQZwkSVIDGcRJkiQ1kEGcJElSAw0siIuIzSPixxHx+4i4IiLeVsdvEBHnRMQ19f/6dXxExOciYk5EXBoRTx1U2iRJkppukDVxDwDvysxtgO2BN0XENsDBwLmZuTVwbv0M8GJg6/p3AHDUANMmSZLUaAML4jLz5sy8uA7fCVwJzABmAcfX2Y4HXlaHZwEnZPErYFpEbDKo9EmSJDXZlPSJi4iZwFOAC4GNM/PmOukWYOM6PAOY17LYjXWcJEmS2gw8iIuItYH/Bd6emX9pnZaZCeQE13dARMyOiNkLFy6cxJRKkiQ1x8qDXHlErEIJ4E7KzNPr6D9FxCaZeXNtLl1Qx88HNm9ZfLM6bgmZeQxwDMC22247oQDw4ej0I+HgD4VLkqThGuTbqQEcC1yZmZ9qmXQmsE8d3gc4o2X83vUt1e2BRS3NrpIkSWoxyJq4HYC9gMsi4pI67j3AkcCpEbE/cAPwyjrtLGBnYA5wD7DfANMmqU/daqMlScM1sCAuM38GRJfJO3WYP4E3DSo9o2R5vyl22j+bnyVJmlz+YoMkSVIDDfTFBknjs+ZSkvRwWBMnSZLUQNbEScsZa/ZWDOazJGviJEmSGsiaOElaTlg7J61YrImTJElqIGviNCWsIZCGw7InLb+siZMkSWoga+IkrK2QJDWPQZwkSVPIh0ZNFoM49WV5/71XSZKaxj5xkiRJDWQQJ0mS1EAGcZIkSQ1kECdJktRAvtjQML5gIEnD57VYo8CaOEmSpAayJk5D45OsVhR+L5ikQTCIk1YA/QbMoxZYG/ysmLqdh+a9tCSDOK1wRi1QkSTp4TCIk5aRtUXNZv5JaiqDOC3XhlXrZm3fxE3FMRulfLHJUNKy8u1USZKkBrImTpKkIbNZf2otL8fbIE5So0zFxXd5ucBLGq1uFJPN5lRJkqQGsiZOaogmPk1aoyVJg2MQp6U0MViQRoXlR8PgA9OKySBOkh4mA7ZmM/DpzmPTDAZxKzhvQt15EdMw+OKGpH4ZxEmShsYHyYnxeKmVQZykxrNmSdKKyK8YkSRJaiBr4iRJqqzVVZMYxEkTYH8UaWmjFvhYTgdjWfK532X7zTsD68IgThpBy/NNaHnet0EZ1jEbteBM0pIM4iSpDwaf0mhbEcuoQZwkqRGGVTM4SjWhy7r8KH3n4LLs34oYsHViECdJ0gpssgMiA6ypYxA3wiwIE+PxUivPh4np93h5XJvN/Fu++D1xkiRJDWRNnCRp0k1VjY81S1qRWRMnSZLUQNbEPUx+IaEkSRoma+IkSZIayJq4AbO/hiRJGgRr4iRJkhpopIK4iHhRRFwVEXMi4uBhp0eSJGlUjUwQFxErAV8EXgxsA7w6IrYZbqokSZJG08gEccB2wJzMvC4z/wacDMwacpokSZJG0igFcTOAeS2fb6zjJEmS1KZxb6dGxAHAAfXjXRFx1YA3uRFw64C3oYkzX0aPeTKazJfRY56MoPjolOTLFpO5slEK4uYDm7d83qyOW0JmHgMcM1WJiojZmbntVG1P/TFfRo95MprMl9FjnoymJubLKDWn/gbYOiK2jIhVgT2AM4ecJkmSpJE0MjVxmflARLwZOBtYCfifzLxiyMmSJEkaSSMTxAFk5lnAWcNOR5spa7rVhJgvo8c8GU3my+gxT0ZT4/IlMnPYaZAkSdIEjVKfOEmSJPXJIK4HfwZs6kTE5hHx44j4fURcERFvq+M3iIhzIuKa+n/9Oj4i4nM1by6NiKe2rGufOv81EbHPsPZpeRERK0XEbyPie/XzlhFxYT32p9QXkYiI1ernOXX6zJZ1HFLHXxURLxzOniw/ImJaRJwWEX+IiCsj4l8sK8MVEe+o167LI+KbEbG6ZWXqRcT/RMSCiLi8ZdyklY2IeFpEXFaX+VxExNTuYZvM9K/DH+XlimuBRwKrAr8Dthl2upbXP2AT4Kl1eB3gasrPr30MOLiOPxj4aB3eGfgBEMD2wIV1/AbAdfX/+nV4/WHvX5P/gHcC3wC+Vz+fCuxRh48G3lCH3wgcXYf3AE6pw9vU8rMasGUtVysNe7+a/AccD7yuDq8KTLOsDDU/ZgDXA2vUz6cC+1pWhpIXzwaeClzeMm7Sygbw6zpv1GVfPMz9tSauO38GbApl5s2ZeXEdvhO4knJhnEW5YVH/v6wOzwJOyOJXwLSI2AR4IXBOZv45M28HzgFeNIW7slyJiM2AXYCv1s8BPA84rc7SnidjeXUasFOdfxZwcmbel5nXA3Mo5UsPQ0SsR7lRHQuQmX/LzDuwrAzbysAaEbEysCZwM5aVKZeZPwH+3DZ6UspGnbZuZv4qS0R3Qsu6hsIgrjt/BmxIatPCU4ALgY0z8+Y66RZg4zrcLX/Mt8n1GeBA4MH6eUPgjsx8oH5uPb4PHfs6fVGd3zyZXFsCC4Gv1Wbur0bEWlhWhiYz5wOfAP5ICd4WARdhWRkVk1U2ZtTh9vFDYxCnkRIRawP/C7w9M//SOq0++fg69RSJiJcACzLzomGnRUtYmdJcdFRmPgW4m9JE9BDLytSqfaxmUQLsTYG1sFZzJC1vZcMgrru+fgZMkyciVqEEcCdl5ul19J9qFTb1/4I6vlv+mG+TZwfgpRExl9Kd4HnAZylNDmPfMdl6fB869nX6esBtmCeT7Ubgxsy8sH4+jRLUWVaG5/nA9Zm5MDPvB06nlB/LymiYrLIxvw63jx8ag7ju/BmwKVT7gxwLXJmZn2qZdCYw9mbQPsAZLeP3rm8XbQ8sqtXlZwMviIj169PxC+o4TVBmHpKZm2XmTMr5f15m7gn8GNi9ztaeJ2N5tXudP+v4PeobeVsCW1M6B+thyMxbgHkR8Zg6aifg91hWhumPwPYRsWa9lo3liWVlNExK2ajT/hIR29d83rtlXcMxzLcqRv2P8ubK1ZQ3hN477PQsz3/AMylV3JcCl9S/nSn9RM4FrgF+BGxQ5w/gizVvLgO2bVnXv1M6BM8B9hv2vi0Pf8COLH479ZGUG8sc4FvAanX86vXznDr9kS3Lv7fm1VUM+W2u5eEPeDIwu5aX71DeoLOsDDdPDgf+AFwOfJ3yhqllZerz4ZuUfon3U2qt95/MsgFsW/P4WuAL1B9NGNafv9ggSZLUQDanSpIkNZBBnCRJUgMZxEmSJDWQQZwkSVIDGcRJkiQ1kEGcpOVCRNw1zvRpEfHGKUjHByPi+ePMs2NE/Oug0yJp+WYQJ2lFMQ0YeBCXme/PzB+NM9uOgEGcpGViECdppETEzIi4MiK+EhFXRMQPI2KNDvNtGRG/jIjLIuKIlvFrR8S5EXFxnTarTjoS2CoiLomIj/eYr307d0XEp2tazo2I6XX8kyPiVxFxaUR8u36zOxFxXETsXofnRsThLdt4bETMBF4PvKOm5VkR8YqIuDwifhcRP5nM4ylp+WUQJ2kUbQ18MTMfB9wBvLzDPJ+l/Aj8Eyjf0D7mr8C/ZeZTgecCn6w/kXMwcG1mPjkz/6vHfO3WAmbXtFwAHFbHnwAclJlPpHzb+2EdlgW4tW7jKODdmTkXOBr4dE3LT4H3Ay/MzCcBLx336EgSBnGSRtP1mXlJHb4ImNlhnh0oP7ED5WeOxgTwkYi4lPITOzOAjTss3+98DwKn1OETgWdGxHrAtMy8oI4/Hnh2l305fZz9APg5cFxE/AewUpd5JGkJKw87AZLUwX0tw38HlmpOrTr9buCewHTgaZl5f0TMpfxW5cOdr59t9jK2L3+nyzU3M18fEc8AdgEuioinZeZtE9yOpBWMNXGSmurnwB51eM+W8esBC2pg9lxgizr+TmCdPuZr9whg9zr8GuBnmbkIuD0inlXH70Vpau3XEmmJiK0y88LMfD+wENh8AuuStIKyJk5SU70N+EZEHASc0TL+JOC7EXEZMBv4A0Bm3hYRP4+Iy4EfAB/tNF8HdwPbRcShwALgVXX8PsDREbEmcB2w3wTS/l3gtPoyxVsoLzlsTWniPRf43QTWJWkFFZkTbRmQpBVHRNyVmWsPOx2S1M7mVEmSpAayJk6SJKmBrImTJElqIIM4SZKkBjKIkyRJaiCDOEmSpAYyiJMkSWoggzhJkqQG+v8otgFpcFKAUQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_df_missing(recovery_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Make second recovery table\n",
    "\n",
    "Here only taking combination of `transport_type x hour`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(39, 32)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>...</th>\n",
       "      <th>22</th>\n",
       "      <th>23</th>\n",
       "      <th>24</th>\n",
       "      <th>25</th>\n",
       "      <th>26</th>\n",
       "      <th>27</th>\n",
       "      <th>28</th>\n",
       "      <th>29</th>\n",
       "      <th>30</th>\n",
       "      <th>31</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hour</th>\n",
       "      <th>route_desc</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">7</th>\n",
       "      <th>Bus</th>\n",
       "      <td>38120</td>\n",
       "      <td>4107332</td>\n",
       "      <td>2539805</td>\n",
       "      <td>1021932</td>\n",
       "      <td>387799</td>\n",
       "      <td>157305</td>\n",
       "      <td>73998</td>\n",
       "      <td>39418</td>\n",
       "      <td>22943</td>\n",
       "      <td>13943</td>\n",
       "      <td>...</td>\n",
       "      <td>342</td>\n",
       "      <td>277</td>\n",
       "      <td>282</td>\n",
       "      <td>303</td>\n",
       "      <td>237</td>\n",
       "      <td>332</td>\n",
       "      <td>228</td>\n",
       "      <td>272</td>\n",
       "      <td>190</td>\n",
       "      <td>6931</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S-Bahn</th>\n",
       "      <td>97805</td>\n",
       "      <td>105802</td>\n",
       "      <td>91405</td>\n",
       "      <td>37178</td>\n",
       "      <td>13751</td>\n",
       "      <td>5179</td>\n",
       "      <td>2250</td>\n",
       "      <td>1176</td>\n",
       "      <td>812</td>\n",
       "      <td>640</td>\n",
       "      <td>...</td>\n",
       "      <td>26</td>\n",
       "      <td>34</td>\n",
       "      <td>30</td>\n",
       "      <td>11</td>\n",
       "      <td>9</td>\n",
       "      <td>9</td>\n",
       "      <td>2</td>\n",
       "      <td>5</td>\n",
       "      <td>7</td>\n",
       "      <td>33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Tram</th>\n",
       "      <td>28177</td>\n",
       "      <td>2773595</td>\n",
       "      <td>941928</td>\n",
       "      <td>315930</td>\n",
       "      <td>114498</td>\n",
       "      <td>47608</td>\n",
       "      <td>22713</td>\n",
       "      <td>11922</td>\n",
       "      <td>6706</td>\n",
       "      <td>3954</td>\n",
       "      <td>...</td>\n",
       "      <td>190</td>\n",
       "      <td>176</td>\n",
       "      <td>209</td>\n",
       "      <td>94</td>\n",
       "      <td>158</td>\n",
       "      <td>119</td>\n",
       "      <td>74</td>\n",
       "      <td>61</td>\n",
       "      <td>73</td>\n",
       "      <td>3004</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <th>Bus</th>\n",
       "      <td>36788</td>\n",
       "      <td>3507612</td>\n",
       "      <td>2269474</td>\n",
       "      <td>1040781</td>\n",
       "      <td>451125</td>\n",
       "      <td>208801</td>\n",
       "      <td>106285</td>\n",
       "      <td>57440</td>\n",
       "      <td>33605</td>\n",
       "      <td>20332</td>\n",
       "      <td>...</td>\n",
       "      <td>355</td>\n",
       "      <td>280</td>\n",
       "      <td>277</td>\n",
       "      <td>200</td>\n",
       "      <td>181</td>\n",
       "      <td>246</td>\n",
       "      <td>273</td>\n",
       "      <td>336</td>\n",
       "      <td>220</td>\n",
       "      <td>6591</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>4 rows × 32 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                    0        1        2        3       4       5       6   \\\n",
       "hour route_desc                                                             \n",
       "7    Bus         38120  4107332  2539805  1021932  387799  157305   73998   \n",
       "     S-Bahn      97805   105802    91405    37178   13751    5179    2250   \n",
       "     Tram        28177  2773595   941928   315930  114498   47608   22713   \n",
       "8    Bus         36788  3507612  2269474  1040781  451125  208801  106285   \n",
       "\n",
       "                    7      8      9   ...   22   23   24   25   26   27   28  \\\n",
       "hour route_desc                       ...                                      \n",
       "7    Bus         39418  22943  13943  ...  342  277  282  303  237  332  228   \n",
       "     S-Bahn       1176    812    640  ...   26   34   30   11    9    9    2   \n",
       "     Tram        11922   6706   3954  ...  190  176  209   94  158  119   74   \n",
       "8    Bus         57440  33605  20332  ...  355  280  277  200  181  246  273   \n",
       "\n",
       "                  29   30    31  \n",
       "hour route_desc                  \n",
       "7    Bus         272  190  6931  \n",
       "     S-Bahn        5    7    33  \n",
       "     Tram         61   73  3004  \n",
       "8    Bus         336  220  6591  \n",
       "\n",
       "[4 rows x 32 columns]"
      ]
     },
     "execution_count": 259,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with gzip.open(\"../data/stop_times_wHour.pkl\", \"rb\") as input_file:\n",
    "    stoptimes = pickle.load(input_file)\n",
    "    \n",
    "distrib_df = pd.DataFrame(d_all).transpose()\n",
    "distrib_to_rm = np.array(distrib_df.iloc[:,range(11)].sum(axis=1) == 11) # missing trips\n",
    "distrib_df = distrib_df.iloc[~distrib_to_rm,:]\n",
    "\n",
    "stoptimes_df = pd.DataFrame(stoptimes)\n",
    "\n",
    "recovery_df2 = distrib_df.join(stoptimes_df)\n",
    "list_bins = [x for x in range(32)]\n",
    "\n",
    "recovery_df2 = recovery_df2.groupby(['hour', 'route_desc'])[list_bins].apply(lambda x : x.astype(float).sum())\n",
    "recovery_df2 = recovery_df2.astype('int')\n",
    "print(recovery_df2.shape)\n",
    "recovery_df2.head(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Third recovery table \n",
    "\n",
    "Takes aggregated transport type distributions to make cumulative probabilities over minutes of delays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 260,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3, 32)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>...</th>\n",
       "      <th>22</th>\n",
       "      <th>23</th>\n",
       "      <th>24</th>\n",
       "      <th>25</th>\n",
       "      <th>26</th>\n",
       "      <th>27</th>\n",
       "      <th>28</th>\n",
       "      <th>29</th>\n",
       "      <th>30</th>\n",
       "      <th>31</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>route_desc</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bus</th>\n",
       "      <td>363268</td>\n",
       "      <td>36683792</td>\n",
       "      <td>23736240</td>\n",
       "      <td>10576138</td>\n",
       "      <td>4388686</td>\n",
       "      <td>1928373</td>\n",
       "      <td>951454</td>\n",
       "      <td>516108</td>\n",
       "      <td>303852</td>\n",
       "      <td>193705</td>\n",
       "      <td>...</td>\n",
       "      <td>5812</td>\n",
       "      <td>4984</td>\n",
       "      <td>4531</td>\n",
       "      <td>4406</td>\n",
       "      <td>3942</td>\n",
       "      <td>4194</td>\n",
       "      <td>3517</td>\n",
       "      <td>3662</td>\n",
       "      <td>3192</td>\n",
       "      <td>63646</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S-Bahn</th>\n",
       "      <td>978390</td>\n",
       "      <td>1105562</td>\n",
       "      <td>826122</td>\n",
       "      <td>300372</td>\n",
       "      <td>108761</td>\n",
       "      <td>41031</td>\n",
       "      <td>17883</td>\n",
       "      <td>9170</td>\n",
       "      <td>5731</td>\n",
       "      <td>4066</td>\n",
       "      <td>...</td>\n",
       "      <td>183</td>\n",
       "      <td>180</td>\n",
       "      <td>115</td>\n",
       "      <td>87</td>\n",
       "      <td>73</td>\n",
       "      <td>56</td>\n",
       "      <td>49</td>\n",
       "      <td>36</td>\n",
       "      <td>42</td>\n",
       "      <td>246</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Tram</th>\n",
       "      <td>247447</td>\n",
       "      <td>25497756</td>\n",
       "      <td>9524532</td>\n",
       "      <td>3371924</td>\n",
       "      <td>1289451</td>\n",
       "      <td>553734</td>\n",
       "      <td>262423</td>\n",
       "      <td>135400</td>\n",
       "      <td>77638</td>\n",
       "      <td>47291</td>\n",
       "      <td>...</td>\n",
       "      <td>1986</td>\n",
       "      <td>1953</td>\n",
       "      <td>1794</td>\n",
       "      <td>1508</td>\n",
       "      <td>1278</td>\n",
       "      <td>1321</td>\n",
       "      <td>1018</td>\n",
       "      <td>1021</td>\n",
       "      <td>867</td>\n",
       "      <td>33074</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3 rows × 32 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                0         1         2         3        4        5       6   \\\n",
       "route_desc                                                                   \n",
       "Bus         363268  36683792  23736240  10576138  4388686  1928373  951454   \n",
       "S-Bahn      978390   1105562    826122    300372   108761    41031   17883   \n",
       "Tram        247447  25497756   9524532   3371924  1289451   553734  262423   \n",
       "\n",
       "                7       8       9   ...    22    23    24    25    26    27  \\\n",
       "route_desc                          ...                                       \n",
       "Bus         516108  303852  193705  ...  5812  4984  4531  4406  3942  4194   \n",
       "S-Bahn        9170    5731    4066  ...   183   180   115    87    73    56   \n",
       "Tram        135400   77638   47291  ...  1986  1953  1794  1508  1278  1321   \n",
       "\n",
       "              28    29    30     31  \n",
       "route_desc                           \n",
       "Bus         3517  3662  3192  63646  \n",
       "S-Bahn        49    36    42    246  \n",
       "Tram        1018  1021   867  33074  \n",
       "\n",
       "[3 rows x 32 columns]"
      ]
     },
     "execution_count": 260,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with gzip.open(\"../data/stop_times_wHour.pkl\", \"rb\") as input_file:\n",
    "    stoptimes = pickle.load(input_file)\n",
    "    \n",
    "distrib_df = pd.DataFrame(d_all).transpose()\n",
    "\n",
    "stoptimes_df = pd.DataFrame(stoptimes)\n",
    "\n",
    "recovery_df3 = distrib_df.join(stoptimes_df)\n",
    "list_bins = [x for x in range(32)]\n",
    "\n",
    "recovery_df3 = recovery_df3.groupby(['route_desc'])[list_bins].apply(lambda x : x.astype(float).sum())\n",
    "recovery_df3 = recovery_df3.astype('int')\n",
    "print(recovery_df3.shape)\n",
    "recovery_df3.head(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Overall aggregate distribution (used when it fails)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 261,
   "metadata": {},
   "outputs": [],
   "source": [
    "last_chance_distrib = np.array(recovery_df3.sum(axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reconstruct cumulative distribution probabilities from multiple distributions to recover data with few/missing points \n",
    "\n",
    "At this point, we have 2 dictionnaries of distributions and 3 recovery dataframes :\n",
    "\n",
    " - `d_real` : contains delay distribution for each keys in form `trip_id + __ + stop_id` calculated from delays with status `geschaetz` or `real` in sbb datasets.\n",
    " - `d_all` : contains delay distributions for each keys in form `trip_id + __ + stop_id`. No filter was applied on status (contains `geschaetz`, `real` __and__ `prognose` = evaluated delay).\n",
    " - `recovery_df` : contains aggregated delay distributions for each combination of `stop_id`, `route_desc` (transport type) and `hour` (time rounded to hour). \n",
    " - `recovery_df2` : contains aggregated delay distributions for each combination of `route_desc` (transport type) and `hour` (time rounded to hour). \n",
    " - `recovery_df3` : contains aggregated delay distributions for `route_desc` (transport type) \n",
    " - `last_chance_distrib` : contains aggregated delay from all lines\n",
    " \n",
    "We will now use these in order to reconstruct the final table with $P(T\\leq t_i)$ for each time points between -1 and +30, using a cumulative probability function as mentionned above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 276,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0%, 3.84%, 7.68%, 11.52%, 15.36%, 19.2%, 23.04%, 26.88%, 30.72%, 34.55%, 38.39%, 42.23%, 46.07%, 49.91%, 53.75%, 57.59%, 61.43%, 65.27%, 69.11%, 72.95%, 76.79%, 80.63%, 84.47%, 88.31%, 92.15%, 95.98%, 99.82%, "
     ]
    }
   ],
   "source": [
    "###################### MAKE CUMULATIVE PROBABILITY TABLE #######################\n",
    "\n",
    "# Load stop_time table, to use its order as a template for our final table \n",
    "with open(\"../data/stop_times_df_cyril.pkl\", \"rb\") as input_file:\n",
    "    stoptimes = pickle.load(input_file)\n",
    "    \n",
    "# declare empty variables \n",
    "size_stop_times = stoptimes.shape[0]\n",
    "n_fail = n_real = n_all = n_recov1 = n_recov2 = n_recov3 = i = 0\n",
    "all_distrib, all_transport_type, all_hours, all_keys, all_failure_keys, all_data_origin = \\\n",
    "    ([] for i in range(6))\n",
    "\n",
    "# Get useful indexes\n",
    "stop_id_idx = stoptimes.columns.get_loc(\"stop_id_general\")\n",
    "trip_id_idx = stoptimes.columns.get_loc(\"trip_id\")\n",
    "transport_type_idx = stoptimes.columns.get_loc(\"route_desc\")\n",
    "\n",
    "for index, row in stoptimes.iterrows():\n",
    "    \n",
    "    trip_id = row[trip_id_idx]\n",
    "    stop_id = str(row[stop_id_idx])\n",
    "    transport_type = row[transport_type_idx]\n",
    "    key = str(trip_id) + '__' + str(stop_id)\n",
    "\n",
    "    # Compute rounded hour using arrival if possible - recover with departure\n",
    "    hour = pd.to_datetime(stoptimes.loc[index]['arrival_time']).hour\n",
    "    if math.isnan(hour): # if arrival is NaT, use departure time\n",
    "        hour = pd.to_datetime(stoptimes.loc[index]['departure_time']).hour\n",
    "    \n",
    "    distrib = np.zeros(31)\n",
    "    keep_trying = True\n",
    "    \n",
    "    # 1) try d_real to get distribution from measured delays\n",
    "    if key in d_real:\n",
    "        distrib = d_real[key]\n",
    "        sum_distrib = np.sum(distrib)\n",
    "        if sum_distrib > 100 :\n",
    "            #summary_df.loc[index, 'distribution'] = distrib\n",
    "            all_distrib.append(distrib)\n",
    "            keep_trying = False \n",
    "            n_real += 1\n",
    "            all_data_origin.append('d_real')\n",
    "    \n",
    "    # 2) try d_all to get distribution from measured + estimated delays\n",
    "    if keep_trying and key in d_all:\n",
    "        distrib = d_all[key]\n",
    "        sum_distrib = np.sum(distrib)\n",
    "        if sum_distrib > 100 :\n",
    "            #summary_df.loc[index, 'distribution'] = distrib\n",
    "            all_distrib.append(distrib)\n",
    "            keep_trying = False\n",
    "            n_all += 1\n",
    "            all_data_origin.append('d_all')\n",
    "\n",
    "            \n",
    "    # 3) try first recovery table using stop_id, transport_type and hour\n",
    "    if keep_trying and (stop_id, hour, transport_type) in recovery_df.index:\n",
    "        distrib = np.array(recovery_df.loc[(stop_id, hour, transport_type)])\n",
    "        sum_distrib = np.sum(distrib)\n",
    "        if sum_distrib > 100 :\n",
    "            #summary_df.loc[index, 'distribution'] = distrib\n",
    "            all_distrib.append(distrib)\n",
    "            keep_trying = False \n",
    "            n_recov1 += 1\n",
    "            all_data_origin.append('recov_tab1')\n",
    "            \n",
    "    # 4) use second recovery table using transport_type and hour \n",
    "    if keep_trying and (hour, transport_type) in recovery_df2.index:\n",
    "        distrib = np.array(recovery_df2.loc[(hour, transport_type)])\n",
    "        sum_distrib = np.sum(distrib)\n",
    "        if sum_distrib > 100 :\n",
    "            #summary_df.loc[index, 'distribution'] = distrib\n",
    "            all_distrib.append(distrib)\n",
    "            keep_trying = False \n",
    "            n_recov2 += 1\n",
    "            all_data_origin.append('recov_tab2')\n",
    "    \n",
    "    # 5) use third recovery table using transport_type only \n",
    "    if keep_trying and (transport_type) in recovery_df3.index:\n",
    "        distrib = np.array(recovery_df3.loc[(transport_type)])\n",
    "        sum_distrib = np.sum(distrib)\n",
    "        #summary_df.loc[index, 'distribution'] = distrib\n",
    "        all_distrib.append(distrib)\n",
    "        keep_trying = False \n",
    "        n_recov3 += 1\n",
    "        all_data_origin.append('recov_tab3')\n",
    "\n",
    "    # Record results in summary\n",
    "    all_keys.append(key)\n",
    "    all_transport_type.append(transport_type)\n",
    "    all_hours.append(hour)\n",
    "\n",
    "    # save number of failure for recovery\n",
    "    if keep_trying:\n",
    "        #print('fail{}'.format(index), end = ', ')\n",
    "        all_distrib.append(last_chance_distrib)\n",
    "        all_failure_keys.append(key)\n",
    "        n_fail += 1 \n",
    "        all_data_origin.append('fail')\n",
    "        \n",
    "    # print progression \n",
    "    if (index % 10000) == 0 :\n",
    "        print('{}%'.format(round(100*index/size_stop_times,2)), end = ', ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 277,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg4AAAEICAYAAAA3Cny5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZhU1bnv8e9PVMAJRUmCOBCVOKCRaGtinI0Sp0Q9wWBCVCQnHpMch5xzTMiN1yGTmskhahSN4pg4xcToVcEBnJVGUUDFEQ7O4oDgECO+94+1ymzLqurdTU80v8/z1MOutdde692rit5vrb26WhGBmZmZWRnLdHUAZmZmtuRw4mBmZmalOXEwMzOz0pw4mJmZWWlOHMzMzKw0Jw5mZmZWmhMHsx5KUkjaoKvjAJC0k6RnO6jtwflcl83Pb5B0cDu1vb2kWYXnsyXt2h5t5/ZmStqpvdpr0M/xki7p6H66K0mflHS7pAWSftvV8Szplu3qAMys55EUwJCIeLKz+46IPcrUKxNjRNwBbNgecUkaDzwbEccU2h/aHm23p1px9gCHAvOAVcJfXrTYPONgZg1VPskvbZbW8+6h1gUeqZc0+LVuHScOZksQSYdI+nvh+ROSriw8nytpWOGQXXOdNySdKUmFumMkPSrpdUk3SVq3sC8kfV/SE8ATuWxvSdNyW3dL+mydGG/Pmw9JWihpZGHff0t6WdILkg4plPeW9BtJ/yvpJUlnS+pbp/1eue48SU8De1XtnyTp3/P2BpImS5qf619eL8bK7RRJP5L0InBBnVssW0l6JI/bBZL65DZHS7qzKpbIMRwKjAJ+mPv7e97/4a2PPAanSno+P06V1Dvvq8RWc/xqjNGn83kvkDQRWKNq/5WSXszjcrukobm8XpxjJT2V23tE0n4N+u4l6f8U6k+VtHbe90VJU3K/UyR9sep1+3l+by2U9HdJq0u6VNKbuf7gQv2NJE2U9JqkWZK+Xiee8cDBhXPaVenWzVWSLpH0JjBa0pqSrs3tPSnpO4U2js9jdkk+p+mSPiPpx/n1mCtpeL0x6XEiwg8//FhCHsB6wBukpH9NYA5pWrmy73Vgmfw8gOuAVYF1gFeA3fO+fYAngY1JtyyPAe4u9BPARKA/0Bf4HPAy8HmgF+kH8Wygd504A9ig8Hwn4H3gp8BywJ7A28Bqef8pwLW5v5WBvwMn1mn7MOAxYO1c/7bc37J5/yTg3/P2n4Cf5PHqA2xXIsaTgd75vHeqjG+uMxuYUej7LuDned9o4M564wCMr9Stam/XvP1T4F7gE8AA4G7gZ2XGr8YY3QP8Lp/HDsAC4JLC/jF5nHsDpwLTCvtqxbk/6f22DDASeAsYWKfvo4HppFs8AjYHVs/j9TpwIOk99438fPXC6/YksD7QD3gEeBzYNde/CLgg110RmAsckvd9jnQrYpM6MX3knIDjgX8C++Zz6gvcDpyV3yfDSP9fdinUfxf4ciGWZ0jvreWA7wDPdPXPh077OdTVAfjhhx+te+QfmFsABwDjgPuBjfIP0WsL9YKPXiivAMbm7RuAbxf2LZMvROsWjt2lsP8PlYtYoWwWsGOdGGtdlN8hX9xz2cvAF/LF5S1g/cK+ber9IAZuBQ4rPB9O/cThojxGa5WM8T2gT1VZdeJQ7HtP4Km8PZrFSxyeAvYs7PsyMLul8atxXuuQkowVC2WXUUgcquqvmuPsVy/OGsdMA/aps29WrX2khOH+qrJ7gNGF1+0nhX2/BW4oPP8KOcEhJS93VLV1DnBcnZg+ck6kROD2wvO1gUXAyoWyE4HxhfoTq2JZCPTKz1fOY7hq2f/HS/LDtyrMljyTSReSHfL2JGDH/JhcVffFwvbbwEp5e13gNKXbDm8Ar5Eu4IMK9ecWttcF/rtSPx+zNulTaFmvRsT7NeIZAKwATC20fWMur2XNqtjmNOjzh6Tzul/pNxjGtBDjKxHxbgt1qvtuzRg0UplBqtd2vfGr1c7rEfFWVVvAh7cSTsq3Et4kJS9QdTujSNJB+tdtqjeATRvUX5uUBNWKq/q1msNH33MvFbbfqfG8+P79fNX7cRTwqXrnUEPxdVwTeC0iFrQitnkRsajwHGq/Hj2OF4SYLXkmkz7xfBr4JenWxSjSp/QzSrYxF/hFRFzaoE5xIVml/i9aH26L5pF+8A6NiOdK1H+BdHGqWKdexYh4kTSNjKTtgJsl3R71f5OizIr76r6fz9tvkRIgcn/VF7GW2n6edEGcWaPt1ngBWE3SioXkYZ1C/98k3aralZQ09CPdMqisf/lInEprX84FvgTcExGLJE0r1K82l3S7YUZVeeX8itYhJYmtNReYHBG7teHYiuJ5Pg/0l7RyIXlYByjzflzqeMbBbMkzGdgZ6BsRzwJ3ALuT7iM/WLKNs4EfFxbF9ZO0f4P65wKHSfq8khUl7SVp5Tr1XyKtuWhRRHyQ2z9F0idyPIMkfbnOIVcAR0haS9JqwNh6bUvaX9Ja+enrpIvFB62Nscr3c9/9Sfe4L8/lDwFDJQ1TWjB5fNVxLfX3J+AYSQMkrQEcC7T6uxciYg7QDJwgafmcMH2lUGVl4B/Aq6RE55ctxLkiadxegbRAlzTjUM95wM8kDcnvlc9KWh34f8BnJH1T0rJKi2Y3Ia3Daa3rclsHSlouP7aStHEb2iIi5pLWlJwoqY/Swt9v04bxXxo4cTBbwkTE46T7q3fk528CTwN3FaZOW2rjGtIiwD/n6eoZQN3vP4iIZtIn9zNIF+AnSff06zkeuDBPI9dc7V7lR7nNe3M8N1P/+xPOBW4iXagfAP7SoN2tgPskLSQtvjwyIp5uY4wVlwETSGP+FPBz+PB1+WmO/Qngzqrj/ghskvv7a412f0664D9MWlz4QKXtNvgmaSHra8BxpLUeFReRpuGfIy1AvLdRnBHxCGm9wT2kpGIz0qLQen5HSu4mAG/m9vpGxKvA3sB/k5KWHwJ7R8S81p5cnhUYTlrn8zzpllxlUWtbfQMYnNu7hrRe4ubFaK/HUl7YYWZmZtYizziYmZlZaU4czMzMrDQnDmZmZlaaEwczMzMrzd/jYD3eGmusEYMHD+7qMMzMlihTp06dFxEf+yI2Jw7W4w0ePJjm5uauDsPMbIkiqea3svpWhZmZmZXmxMHMzMxKc+JgZmZmpTlxMDMzs9KcOJiZmVlpThzMzMysNCcOZmZmVpoTBzMzMyvNXwBlPd705+YzeOz1DevMPmmvTorGzGzJ5hkHMzMzK82Jg5mZmZXmxMHMzMxKc+JgZmZmpTlxMDMzs9KcOJiZmVlpThzMzMysNCcOZmZmVpoTBzMzMyvNiYOZmZmV5sShB5E0W9IadfaNlTSqE2PZSdJ17dTWvpI2aY+2zMxs8Thx6IaUtPdr82VgQju32Vn2BZw4mJl1A04cuglJgyXNknQRMANYW9LRkqZIeljSCYW6f5U0VdJMSYeWaHsVYPmIeKWqfICkibmd8yTNqcxYSPqWpPslTZN0jqRe+TFe0gxJ0yX9INfdQNLNkh6S9ICk9XMXK0m6StJjki6VpFz/2HxeMySNK5SvL+nGfG53SNpI0heBrwK/zrFsm/+tPBZJWnfxXwEzMyvDiUP3MgQ4KyKGAhvm51sDw4AtJe2Q642JiC2BJuAISau30O6uwC01yo8Dbs39XQWsAyBpY2AksG1EDAMWAaNyHIMiYtOI2Ay4ILdzKXBmRGwOfBF4IZd/DjiKNFuwHrBtLj8jIraKiE2BvsDeuXwccHg+t//JY3E3cC1wdEQMi4i78r/DgHOBqyNiTvWJSTpUUrOk5kVvz29heMzMrCz/We3uZU5E3Ju3h+fHg/n5SqRE4nZSsrBfLl87l7/aoN3d+ddFvmg7YD+AiLhR0uu5/EvAlsCUPBnQF3gZ+DuwnqTfA9cDEyStTEomrsntvAuQj7s/Ip7Nz6cBg4E7gZ0l/RBYAegPzJR0GynpuDIfC9C73glJ2hb4Tj6Hj4mIcaREhN4Dh0T9oTEzs9Zw4tC9vFXYFnBiRJxTrCBpJ9IMwjYR8bakSUCfFtrdGvhuK+IQcGFE/PhjO6TNSeslDgO+DhzZoJ1/FLYXActK6gOcBTRFxFxJx5PiXwZ4I88kNA5OGgj8EfhqRCwsd0pmZtYefKui+7oJGCNpJQBJgyR9AugHvJ6Tho2ALzRqRNJQ4LGIWFRj912kiz+ShgOr5fJbgBG5PyT1l7RuXv+wTERcDRwDbBERC4BnJe2b6/aWtEKDkCpJzrx8biMAIuJN4BlJ++d2lJMUgAXAyrl8OeBK4EcR8Xijczczs/bnxKGbiogJwGXAPZKmk9YgrAzcSPrk/ihwEnBv/VYA2CMfU8sJwHBJM4D9gReBBRHxCCkxmCDpYWAiMBAYBEzKtx0uASozEgeSbp88DNwNfKrBeb1BWpswg5QcTSnsHgV8W9JDwExgn1z+Z+BoSQ+Sbmc0AScUFkiu2cIYmJlZO1GEb//2ZJImAgdFxAs19vUGFkXE+5K2Af5Q5lbBkqb3wCEx8OBTG9aZfdJenRSNmdmSQdLUiGiqLvcahx4uInZrsHsd4Ir8nRHvkRYbmpmZ1eXEYSkWEU+QfmXSzMysFK9xMDMzs9KcOJiZmVlpThzMzMysNCcOZmZmVpoTBzMzMyvNiYOZmZmV5l/HtB5vs0H9aPYXPJmZtQvPOJiZmVlpThzMzMysNCcOZmZmVpoTBzMzMyvNiYOZmZmV5t+qsB5v+nPzGTz2+q4Ow8ysU83uoN8m84yDmZmZlebEwczMzEpz4mBmZmalOXEwMzOz0pw4mJmZWWlOHMzMzKw0Jw5mZmZWmhMHMzMzK82Jg5mZmZXmxMHMzMxKW+oSB0mzJa1RZ99YSaM6MZadJF3XTm3tK2mTVh6znKQH2qN/MzNbOiyxiYOS9o7/y8CEdm6zs+wLtCpxALYD7mpLZ5L8d07MzJZCS1TiIGmwpFmSLgJmAGtLOlrSFEkPSzqhUPevkqZKminp0BJtrwIsHxGvVJUPkDQxt3OepDmVGQtJ35J0v6Rpks6R1Cs/xkuaIWm6pB/kuhtIulnSQ5IekLR+7mIlSVdJekzSpZKU6x+bz2uGpHGF8vUl3ZjP7Q5JG0n6IvBV4Nc5lm3zv5XHIknr1jjt3YEbaozFQkmn5HO+RdKAXD5J0qmSmoEjJX1J0oP5PM+X1DvX2zOfz1RJp1dmVSQdn+tNkvS0pCMKff5XPtcZko7KZStKuj6P2QxJI3P5lpIm5/ZvkjSwpdfXzMzaxxKVOGRDgLMiYiiwYX6+NTAM2FLSDrnemIjYEmgCjpC0egvt7grcUqP8OODW3N9VwDoAkjYGRgLbRsQwYBEwKscxKCI2jYjNgAtyO5cCZ0bE5sAXgRdy+eeAo0izBesB2+byMyJiq4jYFOgL7J3LxwGH53P7nzwWdwPXAkdHxLCIuCv/Oww4F7g6IubUOLedgUk1ylcEmvM5T85jULF8RDQBZwLjgZH5PJcFviupD3AOsEeOcUBV2xuRZna2Bo7Lt0u2BA4BPg98AfiOpM+REpvnI2LzPA43SloO+D0wIrd/PvCL6hOQdKikZknNi96eX+MUzcysLZbExGFORNybt4fnx4PAA6SL0pC87whJDwH3AmsXyuup+embNJ3/Z4CIuBF4PZd/CdgSmCJpWn6+HvA0sJ6k30vaHXhT0sqkZOKa3M67EfF2buf+iHg2Ij4ApgGDc/nOku6TNB3YBRgqaSVS0nFl7vMcoO6nbUnbAt8BxtTYNwh4rRBH0QfA5Xn7kjwGFZXyDYFnIuLx/PxCYAfSa/B0RDyTy/9U1fb1EfGPiJgHvAx8Mrd/TUS8FRELgb8A2wPTgd0knSxp+4iYn/vdFJiYx+AYYK3qE4iIcRHRFBFNvVboV3N8zMys9ZbE+9RvFbYFnBgR5xQrSNqJNIOwTUS8LWkS0KeFdrcGvtuKOARcGBE//tgOaXPSp+rDgK8DRzZo5x+F7UXAsvlT+1lAU0TMlXQ8Kf5lgDfyTELj4NL0/R+Br+aLcbXdgZtaaieLwvZbdWuV87HzrdtpxOOStgD2BH4u6RbgGmBmRGyzmHGYmVkbLIkzDkU3AWPyJ3EkDZL0CaAf8HpOGjYiTX/XJWko8FhELKqx+y7SxR9Jw4HVcvktwIjcH5L6S1o3r39YJiKuJn0a3iIiFgDPSto31+0taYUGIVWSnHn53EYARMSbwDOS9s/tKCcpAAuAlXP5csCVwI8KMwLV6s2wQHpfjMjb3wTurFFnFjBY0gb5+YGk2xqzSDMug3P5yPqn+aE7gH0lrSBpRWA/4A5JawJvR8QlwK+BLXL7AyRtUznX/PqZmVknWKITh4iYAFwG3JOn9K8iXTxvJH1yfxQ4iXS7opE98jG1nAAMlzQD2B94EVgQEY+QEoMJkh4GJpJuGwwCJuVp9EuAyozEgaTbJw8DdwOfanBeb5DWJswgJUdTCrtHAd/Ot2FmAvvk8j8DR0t6kHQ7owk4obBAcs1KA5J6ARtExGN1QngL2Dqf8y7AT2vE+C5pXcKVeew/AM6OiHeA75HWI0wlJTQNFxlExAOk9RL3A/cB50XEg8BmwP15LI8Dfh4R75GSmpPzGEzL52tmZp1AEdFyrR5O0kTgoIh4oca+3sCiiHg/f8r9Q5lbBd2ZpO2Ab0XEYXX2L4yIlRaj/ZUiYqEkkRZRPhERp7S1vcXVe+CQGHjwqV3VvZlZl5h90l6LdbykqXkx/EcsiWsc2l1E7NZg9zrAFUrfGfEeabHhEi0i7qT27Yf28h1JBwPLkxauntNCfTMzW0I4cWhBRDxB+pXJpcbizDbk408BumyGwczMOs4SvcbBzMzMOpcTBzMzMyvNiYOZmZmV5sTBzMzMSnPiYGZmZqU5cTAzM7PS/OuY1uNtNqgfzYv5RShmZpZ4xsHMzMxKc+JgZmZmpTlxMDMzs9KcOJiZmVlpThzMzMysNP9WhfV405+bz+Cx13d1GB9a3D91a2bWlTzjYGZmZqU5cTAzM7PSnDiYmZlZaU4czMzMrDQnDmZmZlaaEwczMzMrzYmDmZmZlebEwczMzEpz4mBmZmalOXEwMzOz0pw4mJmZWWndPnGQNFrSmiXqTZLUVKN8I0n3SPqHpP9pcLwk3SpplQZ1jpK0QvnoWybpp5J2rVG+k6Tr2tjmaElnLH50i0/SMEl7dkC7AyTd2N7tmplZY6UTh3xh7YpEYzTQYuLQwGvAEcBvWqi3J/BQRLzZoM5RQLsmDhFxbETc3J5tdhRJbfmjaMNIY9uuIuIV4AVJ27Z322ZmVl/DREDSYEmzJF0EzADWlnS0pCmSHpZ0QqHuQbnsIUkXF46/NZffImkdSf0kzakkIZJWlDRX0nI1+h8BNAGXSpomqa+kY3P/MySNk6TCIQfmejMkbQ0QES9HxBTgny2MxSjgb4WYrs/nMkPSSElHkBKY2yTdlut9Q9L0XOfkQtwLJZ0iaWY+7wENxnh8Pk8k7S7pMUkPAP9Wp/4PJJ2ftzfLfddKZtaUdKOkJyT9qnB83ZgL2yMkjS/Ed7ak+4BfFTtoKRZJywM/BUbm12VkjmdA3r+MpCfz7EGln2ZJj0vaO9fpJenXhffcfxRC+Cvpdas1TofmtpoXvT2/VhUzM2uDMjMIQ4CzImIosGF+vjXpk+SWknaQNBQ4BtglIjYHjszH/h64MCI+C1wKnB4R84FpwI65zt7ATRHxsQt7RFwFNAOjImJYRLwDnBERW0XEpkDffHzFChExDPgecH75YQBgW2Bq3t4deD4iNs/93BgRpwPPAztHxM5Kt09OBnbJY7GVpH3z8SsCzXnMJgPHtdS5pD7AucBXgC2BT9WpehqwgaT9gAuA/4iIt2vUGwaMBDYjXbjXbiHmRtYCvhgR/9WaWCLiPeBY4PL8+l0OXMK/Lva7kmZ5XsnPB5PeW3sBZ+cx+TYwPyK2ArYCviPp07l+M7B9rYAjYlxENEVEU68V+pU4RTMzK6NM4jAnIu7N28Pz40HgAWAjUiKxC3BlRMwDiIjXcv1tgMvy9sXAdnn7ctJFDeCA/LysnSXdJ2l67ndoYd+fcv+3A6tIWrUV7faPiAV5ezqwm6STJW2fk51qWwGTIuKViHiflBjtkPd9UDinS/jXeTeyEfBMRDwREZGP+5iI+IB0++ZiYHJE3FWnvVsiYn5EvAs8AqzbQsyNXBkRixYjlqLzgYPy9hhSwlFxRUR8EBFPAE+TxmQ4cJCkacB9wOqk9xzAyyzebSwzM2ulMves3ypsCzgxIs4pVpB0eCv7vRb4paT+pE/Xt5Y5KH8CPQtoioi5ko4H+hSqRNUh1c8beV/SMvnC9bikLUj35n8u6ZaI+Gkr2qrWmjjKGAIspPFF8x+F7UW0/FoXY+xTte8t6isTy786Sa/bS5J2Ic0uFG811Hr9BBweETfVaK4P8E6Zfs3MrH20drHjTcAYSSsBSBok6ROkC//+klbP5f1z/btJMwqQLhB3AETEQmAKaar7ulqfZgsWACvn7coFbV6OYURV3ZG5/+1I09utubk9C1gvH78m8HZEXAL8GtiiRiz3AztKWkNSL+AbpNsSkMa1Ets3gTtL9P8YMFjS+vn5N2pVktQPOJ00U7B6ZX1ESY1ifknSxnntyX5lGisZS3HMKs4jzahUz2Tsn9c9rE96LWaR3nPfrayBkfQZSSvm+p8hrb0xM7NO0qpV8hExQdLGwD15TeJC4FsRMVPSL4DJkhaRbmWMBg4HLpB0NPAKcEihucuBK4GdWuh2POl+9zukWx/nki4WL5KSj6J3JT0ILEeaBkfSp0j3wlcBPpB0FLBJjd+euD7H8iRpXcCvJX1AWlT53VxnHHCjpOfzOoexwG2kT8XXR8Tfcr23gK0lHUOaTh9JCyLiXUmHAtdLepuUZFVfcAFOAc7MsyLfJi3WvD0iXi7RxwsNYh4LXEd6nZqBlVpqr2QstwFj862GE/M6h2tJtyguqGrvf0nJzSrAYXlMziOtfXhA6U33ClBZl7Ez6XUzM7NOonQ73SQNBC6KiN3aoa2FEVHmwrtUUvq+jVMiYvtC2XjS7NNVrWjndmCfiHi9Ub3eA4fEwINPbWu47W72SXt1dQhmZi2SNDUiPvb9SN3+C6A6S0S8AJyrBl8AZYsvz3hcDfx4MdsZAPyupaTBzMzaV1u+0KdDSDqT9CuRRadFRPV0doeJiCvaqZ2PzTZ0h/PrDiLiJOCkGuWjW9nOK6TvcTAzs07UbRKHiPh+V8fQkXr6+ZmZ2dLBtyrMzMysNCcOZmZmVpoTBzMzMyvNiYOZmZmV1m0WR5p1lM0G9aPZ351gZtYuPONgZmZmpTlxMDMzs9KcOJiZmVlpThzMzMysNCcOZmZmVpp/q8J6vOnPzWfw2I7769v+a5dmtjTxjIOZmZmV5sTBzMzMSnPiYGZmZqU5cTAzM7PSnDiYmZlZaU4czMzMrDQnDmZmZlaaEwczMzMrzYmDmZmZlebEwczMzEpz4mBmZmalOXHoYJJGS1qzRL1JkppqlI+S9LCk6ZLulrR5neMl6VZJq0gaLGlGe8S/uCStKul7HdT2zZJW64i2zcystqUmccgX1q4439FAi4lDA88AO0bEZsDPgHF16u0JPBQRby5GXw1JassfRVsV6JDEAbi4A9s2M7MaenTikD95z5J0ETADWFvS0ZKm5E/xJxTqHpTLHpJ0ceH4W3P5LZLWkdRP0pxKEiJpRUlzJS1Xo/8RQBNwqaRpkvpKOjb3P0PSOEkqHHJgrjdD0tYAEXF3RLye998LrFXndEcBfys87yXpXEkzJU2Q1DfHNEzSvfmcrql8Yi/OeEhaQ9LsvD1a0rWSbgVuqTq/rXI7ffI4zJS0aVVcJwHr5/P6taSLJO1baONSSfvkfv6W43hC0nGFOt+SdH9u4xxJvfKua4Fv1BoMSYdKapbUvOjt+XWGzMzMWqtHJw7ZEOCsiBgKbJifbw0MA7aUtIOkocAxwC4RsTlwZD7298CFEfFZ4FLg9IiYD0wDdsx19gZuioh/VnccEVcBzcCoiBgWEe8AZ0TEVhGxKdA3H1+xQkQMI32KPr/GuXwbuKHOeW4LTK067zPzeb8BfC2XXwT8KJ/TdOA4WrYFMCIidiwWRsQU0sX758CvgEsiovoWyVjgqXz+RwN/JM3CIKkf8EWg8jevt85xfhbYX1KTpI2BkcC2eWwWkZIkckLVW9Lq1QFHxLiIaIqIpl4r9CtximZmVkZbpp6XNHMi4t68PTw/HszPVyJdYDcHroyIeQAR8Vrevw3wb3n7YtLFEeBy0sXsNuAA4KxWxLOzpB8CKwD9gZnA3/O+P+X+b89rFVaNiDcAJO1MShy2q9Nu/4hYUHj+TERMy9tTgcH5Qr1qREzO5RcCV5aIeWJhTKr9FJgCvAsc0VJDETFZ0lmSBpCShKsj4v088TIxIl4FkPQX0rm+D2wJTMl1+gIvF5p8mXQr6NUS52FmZotpaUgc3ipsCzgxIs4pVpB0eCvbvBb4paT+pIvarWUOktSHlGQ0RcRcSccDfQpVouqQyMd9FjgP2KNyYa3hfUnLRMQH+fk/CvsWkS64jbzPv2ag+lTte4v6ViclYMvl4xrVrbgI+BYp6TqkUF7r/EWa9flxnbb6AO+U6NPMzNrB0nCrougmYIyklQAkDZL0CdKFf//KlHdOCADuJl3cIE2P3wEQEQtJn7JPA66LiEUN+lwArJy3KxfkeTmGEVV1R+b+twPmR8R8SesAfwEOjIjHG/QzC1ivwX7ybZbXJW2fiw4EKrMPs0lJEDXiauQc4P+SbuWcXGN/8fwrxgNH5ZgeKZTvJql/Xo+xL3AXaV3FiPw6kfevm7cFfCrHbmZmnWBpmHH4UERMyPfM78nT3guBb0XETEm/ACZLWkS6lTEaOBy4QNLRwCt89NPx5aRp/p1a6HY8cLakd0i3Ps4lLdR8kZR8FL0r6UHSp/cxuexY0qf6s3LM70fEx35tk7ROYCfgyRbiOTjHswLwdOGcfgNcIelQ/rXmoCFJBwH/jIjL8oLFuyXtEhEfzsBExKuS7lL69dAbIuLoiHhJ0qPAX6uavB+4mrQA9JKIaM79HANMyAtS/wl8Hx9q4pkAAA29SURBVJhDSnTujYj3y8RrZmaLTxHVs8O2JJI0ELgoInbr6lhakpOW6cAWeRYESaNJt3D+sxXtnAZcGxG3NKrXe+CQGHjwqYsRcWOzT9qrw9o2M+sqkqbW+qC6tN2q6LEi4gXgXEmrdHUsjUjaFXgU+H0laVgMM1pKGszMrH0tVbcqOpKkM0m/Ell0WkRc0FkxRMQVndVXW0XEzcC6NcrHk27rtKatc9snKjMzK8uJQzuJiO93dQxmZmYdzbcqzMzMrDQnDmZmZlaaEwczMzMrzYmDmZmZlebFkdbjbTaoH83+rgUzs3bhGQczMzMrzYmDmZmZlebEwczMzEpz4mBmZmalOXEwMzOz0pw4mJmZWWlOHMzMzKw0Jw5mZmZWmhMHMzMzK82Jg5mZmZXmxMHMzMxKc+JgZmZmpTlxMDMzs9KcOJiZmVlpThzMzMysNCcOZmZmVpoTBzMzMyvNicNSRNJoSWuWqDdJUlON8n0kPSxpmqRmSdvVOb6vpMmSekkaLOmb7RF/VR/LS7pd0rLt3baZmdXnxKELKOmKsR8NtJg4NHALsHlEDAPGAOfVqTcG+EtELAIGA+2eOETEezmeke3dtpmZ1efEoZPkT96zJF0EzADWlnS0pCn5U/wJhboH5bKHJF1cOP7WXH6LpHUk9ZM0p5KESFpR0lxJy9XofwTQBFyaZwz6Sjo29z9D0jhJKhxyYK43Q9LWABGxMCIi718RCGobBfwtb58EbJ/b+kGeJRhWiOtOSZtLOl7SxZLukfSEpO8U6tQcJ+Cvua9a431onhVpfuWVV+qEaWZmreXEoXMNAc6KiKHAhvn51sAwYEtJO0gaChwD7BIRmwNH5mN/D1wYEZ8FLgVOj4j5wDRgx1xnb+CmiPhndccRcRXQDIyKiGER8Q5wRkRsFRGbAn3z8RUr5JmF7wHnVwol7SfpMeB60szCR0haHlgvImbnorHAHbnPU4A/kmY+kPQZoE9EPJTrfhbYBdgGOFbSmpKG1xqnXH8GsFWtgY6IcRHRFBFNAwYMqFXFzMzawIlD55oTEffm7eH58SDwALAR6QK5C3BlRMwDiIjXcv1tgMvy9sVAZX3B5fxruv6A/LysnSXdJ2l67ndoYd+fcv+3A6tIWjU/vyYiNgL2BX5Wo801gDca9HklsHeeFRkDjC/s+1tEvJPP/TZSslBvnMi3Qt6TtHIrztnMzBaDF5Z1rrcK2wJOjIhzihUkHd7KNq8FfimpP7AlcGuZgyT1Ac4CmiJirqTjgT6FKtW3IT7yPCJul7SepDUqSU72TlU7VB33tqSJwD7A13PMjfqsOU4FvYF36/VnZmbtyzMOXecmYIyklQAkDZL0CdKFf39Jq+fy/rn+3aQZBUj39e+AtO4AmAKcBlyXP4XXswCofDqvXNzn5RhGVNUdmfvfDpgfEfMlbVBZByFpC9JF+9XiQRHxOtArJybVfVacB5wOTMn1K/aR1Cef+075vOqNE7nevFq3ZszMrGN4xqGLRMQESRsD9+Rr8ULgWxExU9IvgMmSFpGm6EcDhwMXSDoaeAU4pNDc5aRbADu10O144GxJ75BufZxLWifwIukiXfSupAeByi0FgK8BB0n6J2lmYWRhsWTRBNKtlJuBh4FFkh4CxkfEKRExVdKbwAVVxz1MukWxBvCziHgeeL7WOAEvAzuT1lqYmVknUe2f+2Ztl2cjfhARB9bZvyYwCdgoIj7IZccDCyPiN63o5y/A2Ih4vFG9pqamaG5uLtusmZkBkqZGxMe+08e3KqzdRcQDwG2SelXvk3QQcB/wk0rS0Bb5tzf+2lLSYGZm7cszDj2QpDOBbauKT4uI6lsDSwXPOJiZtV69GQevceiBIuL7XR2DmZn1TL5VYWZmZqU5cTAzM7PSnDiYmZlZaU4czMzMrDQnDmZmZlaaEwczMzMrzYmDmZmZlebEwczMzEpz4mBmZmalOXEwMzOz0pw4mJmZWWn+WxXW401/bj6Dx17/sfLZJ+3VBdGYmS3ZPONgZmZmpTlxMDMzs9KcOJiZmVlpThzMzMysNCcOZmZmVpoTBzMzMyvNiYOZmZmV5sTBzMzMSnPiYGZmZqU5cbBuRdIRkh6VdGmd/U2STs/boyWd0bkRmpkt3fyV09bdfA/YNSKerbUzIpqB5s4NyczMKjzjYN2GpLOB9YAbJP1I0j2SHpR0t6QNc52dJF3XtZGamS29PONg3UZEHCZpd2Bn4D3gtxHxvqRdgV8CXyvblqRDgUMBeq0yoCPCNTNbKjlxsO6qH3ChpCFAAMu15uCIGAeMA+g9cEi0f3hmZksn36qw7upnwG0RsSnwFaBPF8djZmY4cbDuqx/wXN4e3YVxmJlZgRMH665+BZwo6UF8S83MrNvwD2TrViJicN6cB3ymsOuYvH8SMClvjwfGd1ZsZmbmGQczMzNrBScOZmZmVpoTBzMzMyvNiYOZmZmV5sTBzMzMSnPiYGZmZqU5cTAzM7PSnDiYmZlZaf4CKOvxNhvUj+aT9urqMMzMegTPOJiZmVlpThzMzMysNCcOZmZmVpoTBzMzMyvNiYOZmZmV5sTBzMzMSnPiYGZmZqU5cTAzM7PSnDiYmZlZaYqIro7BrENJWgDM6uo4algDmNfVQVTpjjGB42otx1Ved4wJukdc60bEgOpCf+W0LQ1mRURTVwdRTVJzd4urO8YEjqu1HFd53TEm6L5xgW9VmJmZWSs4cTAzM7PSnDjY0mBcVwdQR3eMqzvGBI6rtRxXed0xJui+cXlxpJmZmZXnGQczMzMrzYmDmZmZlebEwXosSbtLmiXpSUljO6iPtSXdJukRSTMlHZnLj5f0nKRp+bFn4Zgf55hmSfpyS/FK+rSk+3L55ZKWLxnbbEnTc//Nuay/pImSnsj/rpbLJen03MfDkrYotHNwrv+EpIML5Vvm9p/Mx6qFeDYsjMc0SW9KOqorxkrS+ZJeljSjUNbhY1Ovjxbi+rWkx3Lf10haNZcPlvROYdzObmv/jc6xQVwd/rpJ6p2fP5n3Dy4R1+WFmGZLmtaZ46X6PxO6/P3VbiLCDz963APoBTwFrAcsDzwEbNIB/QwEtsjbKwOPA5sAxwP/U6P+JjmW3sCnc4y9GsULXAEckLfPBr5bMrbZwBpVZb8CxubtscDJeXtP4AZAwBeA+3J5f+Dp/O9qeXu1vO/+XFf52D1a+fq8CKzbFWMF7ABsAczozLGp10cLcQ0Hls3bJxfiGlysV9VOq/qvd44txNXhrxvwPeDsvH0AcHlLcVXt/y1wbGeOF/V/JnT5+6u9Hp5xsJ5qa+DJiHg6It4D/gzs096dRMQLEfFA3l4APAoManDIPsCfI+IfEfEM8GSOtWa8+ZPELsBV+fgLgX0XI+R9chvVbe0DXBTJvcCqkgYCXwYmRsRrEfE6MBHYPe9bJSLujfRT6qJWxvUl4KmImNNCrB0yVhFxO/Bajf46emzq9VE3roiYEBHv56f3AmvVGS8A2th/vXNsNF71tOfrVoz3KuBLlU/XLcWV630d+FOjYNt7vBr8TOjy91d7ceJgPdUgYG7h+bM0vqAvtjyN+jngvlz0n3nq8fzClGG9uOqVrw68UbhwtOY8ApggaaqkQ3PZJyPihbz9IvDJNsY1KG9Xl5d1AB/9gd7VYwWdMzb1+ihrDOkTZsWnJT0oabKk7Qvxtrb/tv5/6ejX7cNj8v75uX4Z2wMvRcQThbJOHa+qnwlLwvurFCcOZu1A0krA1cBREfEm8AdgfWAY8AJpyrSzbRcRWwB7AN+XtENxZ/600um/j53vX38VuDIXdYex+ojOGJvW9iHpJ8D7wKW56AVgnYj4HPBfwGWSVumo/mvodq9blW/w0eS0U8erxs+ENrfVFh3ZhxMH66meA9YuPF8rl7U7ScuRfkBcGhF/AYiIlyJiUUR8AJxLmqZtFFe98ldJU5fLVpW3KCKey/++DFyTY3ipMqWa/325jXE9x0enzFszvnsAD0TESzm+Lh+rrDPGpl4fDUkaDewNjMoXBPKtgFfz9lTS+oHPtLH/Vv9/6aTX7cNj8v5+uX5Due6/AZcX4u208ar1M6ENbXXa+6u1nDhYTzUFGKK0Wnt50tT4te3dSb6P+kfg0Yj4XaF8YKHafkBl1fe1wAFKq8U/DQwhLXSqGW++SNwGjMjHHwz8rURcK0paubJNWmA3I/dfWZ1dbOta4KC8wvsLwPw85XkTMFzSankqejhwU973pqQv5DE4qExc2Uc+CXb1WBV0xtjU66MuSbsDPwS+GhFvF8oHSOqVt9fL4/N0G/uvd46N4uqM160Y7wjg1kri1IJdgcci4sMp/c4ar3o/E9rQVqe8v9okOmDFpR9+dIcHabXy46RPFj/poD62I00HPgxMy489gYuB6bn8WmBg4Zif5JhmUfhNhHrxklah309aZHYl0LtEXOuRVq0/BMystEe6P3wL8ARwM9A/lws4M/c9HWgqtDUm9/0kcEihvIl0sXgKOIP8TbQtxLUi6RNjv0JZp48VKXF5Afgn6R7xtztjbOr10UJcT5LudVfeX5XfMvhafm2nAQ8AX2lr/43OsUFcHf66AX3y8yfz/vVaiiuXjwcOq6rbKeNF/Z8JXf7+aq+Hv3LazMzMSvOtCjMzMyvNiYOZmZmV5sTBzMzMSnPiYGZmZqU5cTAzM7PSnDiYmZlZaU4czMzMrLT/D6m/8q787sFwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.barh(np.array(['real / geschaetz','real / geschaetz / prognose','recov_tab1 (stop_id x hour x type)','recov_tab2 (hour x type)','recov_tab3 (type)', 'fail'])[::-1],\\\n",
    "         np.array([n_real,n_all,n_recov1,n_recov2, n_recov3, n_fail])[::-1])\n",
    "plt.title('where the distribution data come from')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 279,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>hour</th>\n",
       "      <th>transport_type</th>\n",
       "      <th>distribution</th>\n",
       "      <th>origin</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2064.TA.26-13-j19-1.24.H__8576240</td>\n",
       "      <td>7</td>\n",
       "      <td>Tram</td>\n",
       "      <td>[15, 445, 31, 4, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0,...</td>\n",
       "      <td>d_all</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2064.TA.26-13-j19-1.24.H__8591353</td>\n",
       "      <td>7</td>\n",
       "      <td>Tram</td>\n",
       "      <td>[1, 442, 41, 11, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0,...</td>\n",
       "      <td>d_all</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2064.TA.26-13-j19-1.24.H__8591039</td>\n",
       "      <td>7</td>\n",
       "      <td>Tram</td>\n",
       "      <td>[0, 415, 68, 11, 2, 0, 0, 1, 0, 0, 1, 0, 0, 0,...</td>\n",
       "      <td>d_all</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2064.TA.26-13-j19-1.24.H__8591121</td>\n",
       "      <td>7</td>\n",
       "      <td>Tram</td>\n",
       "      <td>[0, 403, 79, 12, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0,...</td>\n",
       "      <td>d_all</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2064.TA.26-13-j19-1.24.H__8591417</td>\n",
       "      <td>7</td>\n",
       "      <td>Tram</td>\n",
       "      <td>[0, 388, 79, 24, 4, 1, 0, 1, 0, 0, 0, 0, 1, 0,...</td>\n",
       "      <td>d_all</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                 key hour transport_type  \\\n",
       "0  2064.TA.26-13-j19-1.24.H__8576240    7           Tram   \n",
       "1  2064.TA.26-13-j19-1.24.H__8591353    7           Tram   \n",
       "2  2064.TA.26-13-j19-1.24.H__8591039    7           Tram   \n",
       "3  2064.TA.26-13-j19-1.24.H__8591121    7           Tram   \n",
       "4  2064.TA.26-13-j19-1.24.H__8591417    7           Tram   \n",
       "\n",
       "                                        distribution origin  \n",
       "0  [15, 445, 31, 4, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0,...  d_all  \n",
       "1  [1, 442, 41, 11, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0,...  d_all  \n",
       "2  [0, 415, 68, 11, 2, 0, 0, 1, 0, 0, 1, 0, 0, 0,...  d_all  \n",
       "3  [0, 403, 79, 12, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0,...  d_all  \n",
       "4  [0, 388, 79, 24, 4, 1, 0, 1, 0, 0, 0, 0, 1, 0,...  d_all  "
      ]
     },
     "execution_count": 279,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summary_df = pd.DataFrame([all_keys, all_hours, all_transport_type, all_distrib, all_data_origin],\\\n",
    "                          index = ['key','hour','transport_type','distribution', 'origin']).transpose()\n",
    "summary_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 309,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                        trip_id  stop_id_general\n",
      "167848  35.TA.40-5-Y-j19-1.33.H          8503000\n",
      "167849  35.TA.40-5-Y-j19-1.33.H          8503016\n",
      "167886    8.TA.40-5-Y-j19-1.8.H          8503000\n",
      "167887    8.TA.40-5-Y-j19-1.8.H          8503016\n",
      "167920  16.TA.40-5-Y-j19-1.16.H          8503000\n",
      "167921  16.TA.40-5-Y-j19-1.16.H          8503016\n",
      "251450  18.TA.40-4-Y-j19-1.16.H          8503016\n",
      "251451  18.TA.40-4-Y-j19-1.16.H          8503000\n"
     ]
    }
   ],
   "source": [
    "print(stoptimes.iloc[np.array(stoptimes['route_desc'] == 'Eurocity'),:][['trip_id','stop_id_general']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 311,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>hour</th>\n",
       "      <th>transport_type</th>\n",
       "      <th>distribution</th>\n",
       "      <th>origin</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [key, hour, transport_type, distribution, origin]\n",
       "Index: []"
      ]
     },
     "execution_count": 311,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summary_df.loc[np.array(summary_df['key'] == '1326.TA.26-8-C-j19-1.8.R'),:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write down summary table used to make final distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 240,
   "metadata": {},
   "outputs": [],
   "source": [
    "with gzip.open(\"../data/join_distribution_all.pkl.gz\", \"wb\") as out_file:\n",
    "    pickle.dump(summary_df, out_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 241,
   "metadata": {},
   "outputs": [],
   "source": [
    "list_all_rows = []\n",
    "for index, row in summary_df.iterrows():\n",
    "    distrib = np.array(row['distribution'])\n",
    "    \n",
    "    # get total number of elements \n",
    "    N = np.sum(distrib)\n",
    "        \n",
    "    # make cumulative distribution probabilities\n",
    "    cdf_distrib = np.empty((len(distrib)), dtype=float)\n",
    "    save_x = 0\n",
    "    for x in range(len(distrib)):\n",
    "        cdf_distrib[x] = float(distrib[x])/float(N) + float(save_x)/float(N)\n",
    "        save_x += distrib[x]\n",
    "        \n",
    "    list_all_rows.append(cdf_distrib)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 242,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.03314917, 0.88950276, 0.97237569, 0.98895028, 0.98895028,\n",
       "        0.98895028, 0.99447514, 0.99447514, 0.99447514, 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        ],\n",
       "       [0.        , 0.85082873, 0.95027624, 0.98895028, 0.98895028,\n",
       "        0.98895028, 0.98895028, 0.99447514, 0.99447514, 0.99447514,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        ],\n",
       "       [0.        , 0.77348066, 0.94475138, 0.98895028, 0.98895028,\n",
       "        0.98895028, 0.98895028, 0.99447514, 0.99447514, 0.99447514,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        ],\n",
       "       [0.        , 0.76243094, 0.94475138, 0.98895028, 0.98895028,\n",
       "        0.98895028, 0.98895028, 0.99447514, 0.99447514, 0.99447514,\n",
       "        0.99447514, 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        ],\n",
       "       [0.        , 0.74033149, 0.90055249, 0.97237569, 0.98895028,\n",
       "        0.98895028, 0.98895028, 0.99447514, 0.99447514, 0.99447514,\n",
       "        0.99447514, 0.99447514, 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        , 1.        , 1.        , 1.        ,\n",
       "        1.        , 1.        ]])"
      ]
     },
     "execution_count": 242,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "final_df = pd.DataFrame(list_all_rows)\n",
    "final_df.index = summary_df.index\n",
    "final_np = final_df.to_numpy()\n",
    "final_np[0:5,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last check if all indexes corresponds to `stoptimes` indexes :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 243,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 243,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(np.array(final_df.index == stoptimes.index)) == stoptimes.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write recovery table \n",
    "with gzip.open(\"../data/join_distribution_cumulative_p_3.pkl.gz\", \"wb\") as output_file:\n",
    "    pickle.dump(final_np, output_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Poisson cumulative distribution\n",
    "\n",
    "The Poisson distribution is popular for modeling the number of times an event occurs in an interval of time or space. We modeled a poisson distribution for delays assuming parameter $k$ is the time in minutes (as it was done [here](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0126137), formulas $(4),(5),(6)$).\n",
    "\n",
    "A discrete random variable X is said to have a Poisson distribution with parameter λ > 0, if, for k = 0, 1, 2, ..., the probability mass function of X is given by:\n",
    "\n",
    "$${\\displaystyle \\!f(k;\\lambda )=\\Pr(X=k)={\\frac {\\lambda ^{k}e^{-\\lambda }}{k!}},}$$\n",
    "where\n",
    "\n",
    "e is Euler's number (e = 2.71828...)\n",
    "k! is the factorial of k.\n",
    "The positive real number λ is equal to the expected value of X __and__ to its variance.\n",
    "\n",
    "$${\\displaystyle \\lambda =\\operatorname {E} (X)=\\operatorname {Var} (X)}$$\n",
    "\n",
    "We can approximate E[𝑋]∼$\\mu_i$ for our data $X_i$, if we assume the sample $X_i$ of size N follow the distribution of $X$ meaning $X_i$∼$X$.\n",
    "\n",
    "Poisson-related __assumptions__ :\n",
    "- $k$ is the __delay time in minutes__ and can take values 0, 1, 2, ... (strictly positive and discrete)\n",
    "- We assume our sampling $X_i$ of $X$ is good enough to approximate E[X] ~ $\\mu_i$\n",
    "- The occurrence of one event does not affect probability of others. That is, events occur independently.\n",
    "    - __We assume being late one day is not affecting the delay of the day after__ \n",
    "- The average rate at which events occur is independent of any occurrences. For simplicity, this is usually assumed to be constant, but may in practice vary with time.\n",
    "    - __we assumes delays occurs with a constant rate over time__\n",
    "- Two events cannot occur at exactly the same instant\n",
    "\n",
    "We made a function `poisson_proba` that takes a `trip_id`, a `stop_id`, an `arrival time` and a `departure time` and a dictionnary {key : distribution} to compute a __probability to be at least 2 minutes before departure of next trip__. \n",
    "\n",
    "We make a few __assumptions__ on our side :\n",
    "- We assume that if we have less than 2 minutes for the transfer, we miss it.\n",
    "- We assume the next train is on time.\n",
    "- As for poisson distribution $k$ is strictly positive, we assume trains ahead of schedule were on time ($k=0$)\n",
    "\n",
    "\n",
    "_Question we should address :_\n",
    "- _Is the poisson a reasonable approximation of the binomial distribution in our case ?_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first test the poisson distribution and compare it with our distribution to see how well it fits the data. We will compute $Pr(X = k)$ for each values of k and look at the shape of the poisson distribution compared to the shape of our scaled data. Then, we will compare $\\sum_{k=0}^T Pr(X = k)$ with the cumulative distribution function which directly gives $Pr(k \\leq X)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "An error was encountered:\n",
      "Invalid status code '404' from http://iccluster044.iccluster.epfl.ch:8998/sessions/6821 with error payload: \"Session '6821' not found.\"\n"
     ]
    }
   ],
   "source": [
    "################################# POISSON FIT TEST #########################################\n",
    "\n",
    "# to do .. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are all the functions needed to calculate probability of success for a given transfer. We need the `trip_id`, `stop_id`, `departure_time`, `arrival_time` and dictionnary `d` (pickled load at the beginning of the cell) to be able to compute a probability of success with following function : \n",
    "\n",
    "`poisson_proba(trip_id, stop_id, arrival_time, departure_time, d)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lambda (expectation given distribution):  1.0194769059543685 \n",
      "\n",
      "Probability of success for transfer time = 13.0 minutes :  0.999999999994185\n"
     ]
    }
   ],
   "source": [
    "%local\n",
    "################################# POISSON FUNCTIONS ########################################\n",
    "\n",
    "import pickle \n",
    "import gzip\n",
    "import time\n",
    "import math \n",
    "import datetime\n",
    "import time\n",
    "from scipy.stats import poisson\n",
    "\n",
    "# Load dictionnary\n",
    "with gzip.open(\"../data/distributions.pickle\", \"rb\") as input_file:\n",
    "    d = pickle.load(input_file)\n",
    "\n",
    "# Load dictionnary\n",
    "with open(\"../data/stop_times_array.pkl\", \"rb\") as input_file:\n",
    "    times = pickle.load(input_file)\n",
    "\n",
    "# we take two exemple time in format numpy.datetime64\n",
    "arr_time = times[4][1]\n",
    "dep_time = times[0][1]\n",
    "\n",
    "# Load distribution in dictinonary given a key\n",
    "def get_distrib(key, dico):\n",
    "    if key in dico:\n",
    "        return dico[key]\n",
    "    else:\n",
    "        raise ValueError(\"KEY ERROR: {} not found un distribution dictionnary\".format(key))\n",
    "    \n",
    "# Evaluate lambda parameter assuming it is equal to average \n",
    "def evaluate_lamda(distrib):\n",
    "    # First calculate total number of measures N\n",
    "    N = 0 # by starting at -1 we ignore trains ahead of schedule\n",
    "    for x in distrib:\n",
    "        N += x\n",
    "\n",
    "    lambda_p = 0 # expectation - we want to calculate it\n",
    "    t        = -1 # time = index - 1\n",
    "\n",
    "    for x in distrib:\n",
    "        if t>0:\n",
    "            lambda_p += t*x\n",
    "        t += 1\n",
    "\n",
    "    # calculate lambda - the expectation of x\n",
    "    if N > 0:\n",
    "        lambda_p /= N \n",
    "        print('lambda (expectation given distribution): ',lambda_p, '\\n')\n",
    "        return lambda_p\n",
    "    else : \n",
    "        raise ValueError(\"ERROR : {} distribution has 0 counts\".format(key))\n",
    "        #print('Returning 1 to avoid later problem... \\n')\n",
    "        return 1\n",
    "\n",
    "# process time given as string in format 'hh:mm' - not needed\n",
    "def process_time_str(str_time):\n",
    "    x = time.strptime(str_time,'%H:%M')\n",
    "    return datetime.timedelta(hours=x.tm_hour,minutes=x.tm_min,seconds=x.tm_sec).total_seconds()\n",
    "\n",
    "# Calculate transfer time given two times in string format 'hh:mm'\n",
    "def get_transfer_time(arr_time, dep_time, delta=2.0):\n",
    "    diff_time_min = (arr_time - dep_time).astype('timedelta64[m]') / np.timedelta64(1, 'm')\n",
    "    return diff_time_min - delta\n",
    "\n",
    "# Calculate poisson probability of success for a given transfert \n",
    "# for a given trip_id, stop_id, arrival/departure times and dict\n",
    "def poisson_proba(trip_id, stop_id, arr_time, dep_time, dico):\n",
    "    # Generate key from trip_id / stop_id \n",
    "    key = str(trip_id) + '__' + str(stop_id[0:7]) # 7 first char to be sbb-compatible\n",
    "\n",
    "    # Get distribution from dictionnary\n",
    "    distrib = get_distrib(key, dico)\n",
    "    \n",
    "    # Calculate transfer time at disposal \n",
    "    T = get_transfer_time(arr_time, dep_time)\n",
    "    \n",
    "    # Get lambda value to calculate proba\n",
    "    lambda_p = evaluate_lamda(distrib)\n",
    "\n",
    "    # Get proba\n",
    "    if T > 2:\n",
    "        poisson_p = poisson.cdf(T, lambda_p)\n",
    "    else : \n",
    "        poisson_p = 0.0 # if we have less than 2 minutes, we miss it\n",
    "        \n",
    "    print('Probability of success for transfer time = {} minutes : '.format(T),poisson_p)\n",
    "    return poisson_p\n",
    "\n",
    "# Mock exemple of probability calculations with given inputs\n",
    "trip_id = '1286.TA.26-32-j19-1.12.H'\n",
    "stop_id = '8591184'\n",
    "\n",
    "# we take two exemple time from stop_times_array in format numpy.datetime64\n",
    "arr_time = times[3][1]\n",
    "dep_time = times[0][1]\n",
    "\n",
    "Pr = poisson_proba(trip_id, stop_id, arr_time, dep_time, d)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}