{ "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": 160, "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": 161, "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": 162, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "len dict_real : 16760\n", "[('1.TA.26-20-j19-1.1.R__8503000', array([ 0, 25, 78, 18, 5, 2, 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])), ('1.TA.26-20-j19-1.1.R__8503003', array([ 0, 83, 31, 9, 4, 1, 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])), ('1.TA.26-20-j19-1.1.R__8503101', array([ 3, 89, 22, 6, 5, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,\n", " 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])), ('1.TA.26-20-j19-1.1.R__8503104', array([ 0, 83, 21, 13, 4, 3, 2, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0,\n", " 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])), ('1.TA.80-173-Y-j19-1.1.H__8502276', array([ 0, 315, 61, 17, 2, 1, 0, 0, 0, 0, 0, 0, 0,\n", " 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,\n", " 0, 0, 0, 0, 0, 0]))]\n", "len dict_all : 240567\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/d_real.pkl.gz\", \"rb\") as input_file:\n", " d_real = pickle.load(input_file)\n", "\n", "with gzip.open(\"../data/d_all.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": 63, "metadata": {}, "outputs": [], "source": [ "def plot_data_points_hist(dico, binwidth = 100):\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", " 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": 64, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoMAAAFNCAYAAAB7ZGXXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dedwdVX348c+XAIpsCRJTJEhQUy1qRUgBKyoVhQAVqCuKEigFFbX4s7XFpaKAFrVqpSIWBQmiIOICVShGBK0LSFAWAZXIIkEggYRN3IDv749zLkxunvs8N8lznyXzeb9e93VnzpyZOXNm+94zM3ciM5EkSVI7rTPeBZAkSdL4MRiUJElqMYNBSZKkFjMYlCRJajGDQUmSpBYzGJQkSWoxg8GGiMiIeOoEKMfFEfEP4zTvDSLifyLinoj4ch/5d42IxWNRtjUREQdExLfGuxxjJSLOj4h5g5pGRMyq+8u6azKP8RIR90fEk3sMOygivt9j2IRe7uHKPlFExJNq/U8ZJk/fx+KIeF9EnD56JVw1k6HOR9tw+09bRMQeEfH18S7HcCLiMRHx84iYPlLeSREM1g2v83k4In7X6D+gxziTIkiZgF4BzAAen5mvHM0Jj+dBMzO/kJm795N3Ih/c+z3xZeaemTl/Tea1ptOIiCeu6j4YEadGxLGrO89+ZeZGmXnDoOezJhqB53ld6adHxPvGqVhrLDN/Xev/IRjfH7+rarx/DNT9449d58QpjeG71ZP/AxFxUURsPYhyrMr+M4hGloh4YZ3usV3pT46Ib0TEfRFxZ0R8uDHs9Ii4LSLujYhfdm9zEfGqiLiujnttROw3QjE+ABzXGP+YiLg6Ih4caf+MiKkRMT8iltTP+7qGXxQRS2tZr4yIfRvDnh0R19Tle3sjfb2IuDQituqkZeYfgFOAI0dYlskRDNYNb6PM3Aj4NfDSRtoXxrt8E1UUq7qOtwZ+mZkPDqJM6t/qnHBWc50Pyl7A/453IdYCO0XEX493ITRhfLh5TmwE1ZsDXwX+DdgMWAh8aRzLORARsR7wCeDSrvT1gQXAd4A/A2YCzR/O/w7MysxNgH2AYyNihzruljXv24FNgHcAX4yIJ/Qow18Bm2bmJY3kRcC/AN/sYzE+DjwOmAXsCLw+Ig5uDD8C2KKW9TDg9IjYorEc/ww8G3h3RPxZTX878JXMvKVrXl8E5kXEY4YtUWZOqg9wE/Di2v0Y4D+B39TPf9a0DYHfAQ8D99fPE2ul/wi4G7gN+CSwfmPaCTy1x3wvBo4BfgDcB3wL2LwO2xVYPEw53wd8mbKx3QdcDfw58E5gCXALsHvXvP4d+DFwL3AOsFlj+M7AD+tyXAns2jXuB2o5fzfU8gB/UfPdDVwD7FPT3w/8EfhTrbNDhhh3A+BUYDlwLWWnWdwYfiTwq7qc1wJ/15jn74GH6rTvrul7Az+ty3kL8L5h1v2uwGLgXcCdtY4PaAzfFDgNWArcDLwHWKcOOwj4fte6fiNwfa2HE4AYppx71eW5D7gV+OceZTyo1v0ngXuAnwO7dZXxZMr2dytwLDCla9yPA3cBx3ZNe27X+rmy1zqvaf/QT5mGqe/mNKYA/1Hr/QbgzbUO1x1m/K8CLxsiPeoyLqnr/WrgmZSD3p/qMt4P/M9w22sddirwacpJ4D7gu8DWfSzbI/s68Hjg3FqWH1P28+/3GG9Wc7mBl1O2w2dSflx3tv+7gLOo+y3lBPHWrmldRd0/RpjXvwIXNdJPp8d+wsrb+UeA79ftbshtD1gfWAY8qzHeE4AHgOnA5sA3av0vA/6Pul91zfv9wH/V7vWA3wIfaRw3fk8JUh6pQ8p2+1Addj/wyeH2zx7L/D7g9D6PjwdRtt/7gBupxw/KPvNdyv5xJ/ClHvP6dS1b57zy3E6dU/aP5XW6ezbGORi4rs7zBuANQxzT/omyP9wGHDzMNnEqXceFxrDDgB82+jvnwaf3yH8T5Rx0bS3354DHNoYfSglwllH2jyf22H9Orevnm3UZLwWeUod9r+b9ba2vV/e7PQ1TB0cCH+6ui7r8/9fnNJ5W6/pVtX8nYElXnqXAc3uM/17gsz2G9dw/G3nuBP6q0f+uXmWnxC2/B3as/dcBj6ndl9ThW1OOXev1mMb1wAuHLVO/K2CifFgxyDq6VsYTKAetHwLH1GG7snKAtgPlQLEu5YB0HfC2oTbwIeZ7MeUg/+eUA9vFwHHDzKtZzvfVlblHnfdplAPGuykHzUOBG7vmdSvlBLMh8BXqwQ7YknKi2Yty8nlJ7Z/eGPfXwDPqvNbrKtd6lB38XZSTwIsoO/DTGmU9fZj6P46y824GbAX8jBWDwVdSAu91KDv+bym/cKDrRNWou2fV/H8J3AHs12PeuwIPAh+jBP0vrNPvlP00SuC8cV2/v6QGtN3zruv6G8BU4EmUHX/uMOW8DXh+7Z4GbN+jjAfVMv6/WtevppxgOkHB14D/ruv1CZQd+A1d4761rrsNhpj+SutnqHXOysFgzzINs66b03gjJYjcqq77ixgmGKzzuRPYeIhhewCX17rvBOCdbeRUVjzAj7S9nlr7X1C3iU90r7se5WuezM6kBG4bUva5W3tNgxUDmYNr2TrTOYJyPJpZy/LfwBl12KuASxvTeTZlv11/mDJ25rVxLVPneDJiMEjZnz4DXAA8ro9t71PAhxrTOYJHg/F/pwTc69XP8xkiMKvr5ura/deU4+WljWFXdtdh93bWz/453D7BMMfHutz3NradLYBn1O4zKMfjdYDHAruMtP676vxPlOP4FOBNlMaJqMP3Bp5C2dZfSAmyt+86ph1d63avOnxaj/mfSgmgllH2oZc3hn0COLEr/8+aebqG3VSHd/bpH1D3vbq+7gS2p2zL/wV8r8f+c2qt4x0p+8UXgDOHyjvS9kTZDj81zD6xNeW4vhErHytOAT4PnF/LfjGNHziN6T9Qy/QTYKOaPoXyY2Cf2r0fJUjfsEc5vgy8o8ewfoPBHRv97waWd+X5BiVuSMoVlnUa834p5ThzO+XH7NcZJtijBPP/OGyZhhs4ET+sGGT9CtirMWwP4KbGTrZ4hGm9Dfhar422K+/FwHsa/YcD/9trXqwcDC5oDHsp5VdSp0Vo4zrvqY15HdfIvy2ltWQKpZXg813zugCY1xj36GGW+fl1A1qnkXZGZ+Nl5GDwBhoHZcqvsZ71DFwB7Fu7D2KEEzWldffjPYbtSjlwbthIO4tyWWRKraNtG8PeAFw81Lxrfe/SNZ0je5WTEmy9AdhkhPIfRONEUNN+DLyeci/mH2gEecBrqK0+ddxfjzD9ldbPUOuclYPBIcs0wrya0/gO8MbGsN0ZPhjcDbiwx7AXUQ7oO9PVIsDKB/iRttdTWfHEsxGlpWmrEZYtKa1BUygn8qc3hn2w13bKo8HAP1NaVGY2hl3Hiq3AW9Rpr0sJMJYDs+uw/2CYk17XvNalHG8uqekjBYOXUi4PfoUabPax7e1E2cY7J+WFPNpqcjTlR9aQx8bG9Dqtf4+ntN68i3JC3YjSanh893J1b2f97J/D7RMMc3ykBIN3U1pzN+jKcxpwUnN9jrROuup8UaP/cTXPn/WYxteBI2r3rpTWu+b0lgA79xh3+1q/61ICx/uA59VhJ9M4b9S0HwAH9ZjWTay4T+8F/KoxrQ937Vd/olxm7ayfZjD42a7p/Lx7X2v097U99SjzOcCrG/NtHiu+Vcu4J+WH4zso56v1u6YxBdiFcuVovUb6IZTz8oOUgHHvYcqxoFl3XcP6CQZPp1w52ZhyHPoV8Ich8q1Xl+ftjbStgfMowexrKAHs5yk/ms6hBLWv7JrOF4D3DlemiXJv0ep6IuVyYMfNNW1IEfHn9ebS2yPiXspBf/NVmN/tje4HKDtIv+5odP8OuDPrvR61n67pNa/730zZKDanbAivjIi7Ox/Khr1Fj3G7PRG4JTMf7pr+ln0uxxOHKNsjIuLAiLiiUbZnMkwdR8ROjZtl76G0QA23TpZn5m+75v/EOs56rLw9DLdcq7I+X045yN0cEd+NiOcOk/fWrHtgVxm3rmW8rVE//01ppekYbt0NZ6TxepWpX8Ou9yHsRTlgrSQzv0O5ZH0CsCQiToqITYab7wjb6yPlysz7Ka0m/S7bdMqJdVWWDcqJ5oTMbD4gszXwtca6vY4SmM7IzN9TArTX1Xs6X0M5gPfrs8CMiHhpH3mfCuwLvD8z/9goW89tLzMvpewDu0bE0+s0zq3jfoTSAvqtiLghIoa8GT0zf0cJIl9Iaan9LuVqzfNq2ndXYXlh9Y63PY+P9bjxasox5raI+GZdVij3egXw43pz/t+vblkz84HauRFAROwZEZdExLJanr1Y8Rh3V654j3bPZc3Mn2TmXZn5YGaeRznJv6wOvp9yv1vTJpSAsZfu7b6z36xwbq371V30Pp6uyrrqa3vqVrf9jTOz132Qv6P8iDu/bvf/QQmc/6KZKTMfyszvU1rW3lSn/WLKpeddKYHkC4HPRsR2Pea1nBLIra5/rOW9nhLAnUH54bSCzPxTZp4P7B4R+9S0mzNzr8zcvo57DOXH6X9QjjH7AB+LiM0ak9qY8kOop8keDP6GsvN3PKmmQfk10u1EyqWu2VluzHwX5QCwpn5L+TUIQH26a8RHuUewVaP7SZRfPHdSdt7PZ+bUxmfDzDyukX+oZe/4DbBV10MGT6JchurHbUOUDYD65NpngLdQnkaeSrkM0anjocr1RcpJZ6vM3JRy+WC4dTItIjbsmv9vKHXzJ1beHvpdrqaVypmZl2XmvpST59cpLRW9bBkRzWXolPEWSuvM5o11t0lmPmO4eY9Utj7H61WmfvVc7z30DAYBMvP4zNyB0ur955TgClZejn6210fKFREbUS559btsSyktAauybFBaRt8TES9vpN1CuVesuW8+NjM7ZZ0PHEBpNX0gM3/UZxmpJ7f3Uw78Ix2zrqNcwj4/Ip7WKNtI29584HWUVuyzawBLZt6Xmf+UmU+mnGjeHhG79Zj3dyktv88BLqv9e1AuIX6v1+KNsDyrYtjjY2ZekJkvofx4/jnleEVm3p6Zh2bmEylXAD7V4wnYVSprvWn/K5QT9Yx6TDyP0TnvdMrTmdY1lNsPOvPekHJ5+pphxu/e7jv7zQrn1jqtx7N6x9MVC7xq21PTbsCc2phzOyWwf1tEnFOHX8WqrZ91KfUDsB3lMvjCzHw4My+jtLC/uMe4V1GOW6slM5dl5gGZ+Wd1H1yHcrWmn7I2vRf4TGbeQbndamFm3kMJLJvb719Q7p/tabIHg2dQDsjT65NU7+XRp4fuAB4fEZs28m9MuWfk/vqL8E2jVI5fAo+NiL3rk07vodxnsSZeFxHbRsTjKM3qZ9eWxNOBl0b5j6MpEfHYKH+jM7PP6XZaAP6lPoq+K+Wy9Zl9jn8W8M6ImFbn+dbGsA0pO+NSgPp01DMbw+8AZtanvjo2BpZl5u8jYkfgtX2U4f0RsX5EPB/4W+DLtW7OAj4QERvXwPTtrPg0Wb9WKGed1wERsWlm/omyDT08zPhPAP6x1u8rKTvieZl5G+VSxkcjYpOIWCcinhIRL1zFss1ajSeGhyzTKox/Vh1/ZkRMY5i/KoiIbSg3OF/XY/hf1RbhzkMGv+fR+rwDaP5/WT/b614RsUtdX8dQLqf21cJat5uvAu+LiMdFxLaUS4ojuYbyQM8JnV/slB8yH6jbHvW49MhfQtTg72Hgo6xaq2DH5ymXm+eOlDEzz6D82P12RDylz23vdODvKAHhaZ3EiPjbiHhq/TFxD6W1s9f2/13gQODaGsBeDPwD5Z7opT3G6V7na6Ln8TEiZkTEvjWw+QOlJe3huoyvbBxDl1OOY0Mt49Ka3m9516ecC5YCD0bEnpQfEqslIl4RERvV9bc7ZV11WnC/BjwzIl4eEY+lnA+vysyfDzPJN9e62Yxy31qn1e0M4OCI2K4GtB+k3P9502oUe4X1u4rbU9O/UQKw7ernXEow33kK93Rg54h4cW2QeRulkeC6iHhCROxf625KROxBaZ2/sI57GfD8TktgRDyHcovKVT3Kch6l9fAR9fj0WEpctW7d9ob8L8267z2+lmVPyu1Wx9ZhT4/SmrxBnebreLSlvTmNbSktmSfWpBuBF0XEDGA25baPzpPSm1HuZ+5psgeDx1IuS1xFeSLxJzWNugOcAdwQ5XLBEylNqa+lNJt/hlF67L5G4odTLuXcSjnBrel/HH6eck/E7ZQTwD/Wed1CuQT0LsoB5hZKq0pf67IeoF9KuQ/hTsoNtQeOcMBoej/l8sGNlJPLIye1zLyWcqL7EeUA8CzKPSsd36GcRG+PiDtr2uHA0RFxH+XgNVyLG5T6WE755foFyn0bnbK/lVL3N1Buov8i5abiVTVUOV8P3BTl9oI3Ulp4ermUsjPeSXla8hWZeVcddiDlBNF5gu9sVrzEP5LOH4HfFRE/WYXxhitTPzoPI1xJ2c++OkzevRk+0NykTm85ZVu6i3LpCMq9StvWffbrfW6vXwSOolwe3oFyglwVb6Fc1rqdss99rp+RMvNKyo+Rz9QD+icoJ6hv1e35Esq9eE2nUfaL0wEi4tMR8ek+5/cQZR/ZbKS8Nf98yg/J70TELEbY9uqx5SeUQOj/GpOaDXybEjz9iHKv40U9ZvtDyr2DnVbAaynBfq9WQSj19oqIWB4Rx/ezbL2McHxch/ID8TeUbeWFPNog8FfApRFxP2UdHpFD/I9evQT8AeAHdRvdeYTy3Ec5dp9FqfPX8mjwtjqOoJxj7qbsM4dm5sV1Xkspt7N8oM5rJ2D/Eab3Rcpx/AbKfWud8+e3KcHXVyhXBZ7Sx7R6eR8wv9bXqxhmexpuf6gtird3PpTLrL/NzGV1+C8o+/6nKcu/L+WfB/5I2abfRDkvL6e01L4tM8+t4363lvPsuu9+BfhgZg75ooLM/AlwT0Q09+/P1DK9hhJY/45y3iAinl+3rY4dKDHLfZQHag7IzE4LbtSyLKFsw0dQ7pPsPt6fQNlOO7ebvZOyrV1Ty965dP9aYH6W/xzsqXOzsDTh1Vah0zOz31bQMRcRB1Fuht9lvMvSMdZlivInyZ+s9zQNel6nUh5ges+g5zUaIuJA4LCJtH00RcQpwG8mS31q9UXETZTjwrfHuyyTUW2ZPTwzR/pz6nFTW3WvBF6QmUuGyzshX6kkaVK7mPLXM2qIcsvH4ZTWzQmnth6+jHK/n6Rh1FbDCf2K09oa+PQRMzL5LxNLWgOx4mutmp/nr+40M/PDWZ4sHTedyzJDfcapPHtQLvncQbk0N6FExDGUh70+kpk3jnd5JI0tLxNLkiS1mC2DkiRJLWYwKEmS1GKte4Bk8803z1mzZo13MSRJkkZ0+eWX35mZa/oii2G1LhicNWsWCxcuHO9iSJIkjSgi+nlF5hrxMrEkSVKLDSwYjIinRcQVjc+9EfG2iNgsIhZExPX1e1rNHxFxfEQsioirImL7xrTm1fzXR8S8RvoOEXF1Hef4+nobSZIk9WlgwWBm/iIzt8vM7SivXnmA8u7EI4ELM3M25b2AnXec7kl5Tc1synv6TgSo70w8ivJqnR2BozoBZM1zaGO8Ed/ZKUmSpEeN1WXi3YBfZebNlPcFzq/p84HOq1z2BU7L4hJgakRsAewBLMjMZZm5HFgAzK3DNsnMS7L8WeJpjWlJkiSpD2MVDO4PnFG7Z2TmbbX7dmBG7d6S8lLxjsU1bbj0xUOkS5IkqU8DDwYjYn1gH+DL3cNqi97AX4ESEYdFxMKIWLh06dJBz06SJGnSGIuWwT2Bn2TmHbX/jnqJl/q9pKbfCmzVGG9mTRsufeYQ6SvJzJMyc05mzpk+faB/1SNJkjSpjEUw+BoevUQMcC7QeSJ4HnBOI/3A+lTxzsA99XLyBcDuETGtPjiyO3BBHXZvROxcnyI+sDEtSZIk9WGgfzodERsCLwHe0Eg+DjgrIg4BbgZeVdPPA/YCFlGePD4YIDOXRcQxwGU139GZuax2Hw6cCmwAnF8/kiRJ6lOU2/baY86cOekbSCRJ0mQQEZdn5pxBzsM3kEiSJLVY695NLJh15DdX6L/puL3HqSSSJGm82TIoSZLUYgaDkiRJLWYwKEmS1GIGg5IkSS1mMChJktRiBoOSJEktZjAoSZLUYgaDkiRJLWYwKEmS1GIGg5IkSS1mMChJktRiBoOSJEktZjAoSZLUYgaDkiRJLWYwKEmS1GIGg5IkSS1mMChJktRiBoOSJEktZjAoSZLUYgaDkiRJLWYwKEmS1GIGg5IkSS1mMChJktRiBoOSJEktZjAoSZLUYgaDkiRJLWYwKEmS1GIDDQYjYmpEnB0RP4+I6yLiuRGxWUQsiIjr6/e0mjci4viIWBQRV0XE9o3pzKv5r4+IeY30HSLi6jrO8RERg1weSZKktc2gWwY/AfxvZj4deDZwHXAkcGFmzgYurP0AewKz6+cw4ESAiNgMOArYCdgROKoTQNY8hzbGmzvg5ZEkSVqrDCwYjIhNgRcAJwNk5h8z825gX2B+zTYf2K927wuclsUlwNSI2ALYA1iQmcsyczmwAJhbh22SmZdkZgKnNaYlSZKkPgyyZXAbYCnwuYj4aUR8NiI2BGZk5m01z+3AjNq9JXBLY/zFNW249MVDpEuSJKlPgwwG1wW2B07MzOcAv+XRS8IA1Ba9HGAZAIiIwyJiYUQsXLp06aBnJ0mSNGkMMhhcDCzOzEtr/9mU4PCOeomX+r2kDr8V2Kox/syaNlz6zCHSV5KZJ2XmnMycM3369DVaKEmSpLXJwILBzLwduCUinlaTdgOuBc4FOk8EzwPOqd3nAgfWp4p3Bu6pl5MvAHaPiGn1wZHdgQvqsHsjYuf6FPGBjWlJkiSpD+sOePpvBb4QEesDNwAHUwLQsyLiEOBm4FU173nAXsAi4IGal8xcFhHHAJfVfEdn5rLafThwKrABcH79SJIkqU8DDQYz8wpgzhCDdhsibwJv7jGdU4BThkhfCDxzDYspSZLUWr6BRJIkqcUMBiVJklrMYFCSJKnFDAYlSZJazGBQkiSpxQwGJUmSWsxgUJIkqcUMBiVJklrMYFCSJKnFDAYlSZJazGBQkiSpxQwGJUmSWsxgUJIkqcUMBiVJklrMYFCSJKnFDAYlSZJazGBQkiSpxQwGJUmSWsxgUJIkqcUMBiVJklrMYFCSJKnFDAYlSZJazGBQkiSpxQwGJUmSWsxgUJIkqcUMBiVJklrMYFCSJKnFDAYlSZJazGBQkiSpxQYaDEbETRFxdURcERELa9pmEbEgIq6v39NqekTE8RGxKCKuiojtG9OZV/NfHxHzGuk71OkvquPGIJdHkiRpbTMWLYN/k5nbZeac2n8kcGFmzgYurP0AewKz6+cw4EQowSNwFLATsCNwVCeArHkObYw3d/CLI0mStPYYj8vE+wLza/d8YL9G+mlZXAJMjYgtgD2ABZm5LDOXAwuAuXXYJpl5SWYmcFpjWpIkSerDoIPBBL4VEZdHxGE1bUZm3la7bwdm1O4tgVsa4y6uacOlLx4iXZIkSX1ad8DT3yUzb42IJwALIuLnzYGZmRGRAy4DNRA9DOBJT3rSoGcnSZI0aQy0ZTAzb63fS4CvUe75u6Ne4qV+L6nZbwW2aow+s6YNlz5ziPShynFSZs7JzDnTp09f08WSJElaawwsGIyIDSNi4043sDvwM+BcoPNE8DzgnNp9LnBgfap4Z+Ceejn5AmD3iJhWHxzZHbigDrs3InauTxEf2JiWJEmS+jDIy8QzgK/Vf3tZF/hiZv5vRFwGnBURhwA3A6+q+c8D9gIWAQ8ABwNk5rKIOAa4rOY7OjOX1e7DgVOBDYDz60eSJEl9GlgwmJk3AM8eIv0uYLch0hN4c49pnQKcMkT6QuCZa1xYSZKklvINJJIkSS1mMChJktRiBoOSJEktZjAoSZLUYgaDkiRJLWYwKEmS1GIGg5IkSS1mMChJktRiBoOSJEktZjAoSZLUYgaDkiRJLWYwKEmS1GIGg5IkSS1mMChJktRiBoOSJEktZjAoSZLUYgaDkiRJLWYwKEmS1GIGg5IkSS1mMChJktRiBoOSJEktZjAoSZLUYgaDkiRJLWYwKEmS1GIGg5IkSS1mMChJktRiBoOSJEktZjAoSZLUYgaDkiRJLTbwYDAipkTETyPiG7V/m4i4NCIWRcSXImL9mv6Y2r+oDp/VmMY7a/ovImKPRvrcmrYoIo4c9LJIkiStbcaiZfAI4LpG/4eAj2fmU4HlwCE1/RBgeU3/eM1HRGwL7A88A5gLfKoGmFOAE4A9gW2B19S8kiRJ6tNAg8GImAnsDXy29gfwIuDsmmU+sF/t3rf2U4fvVvPvC5yZmX/IzBuBRcCO9bMoM2/IzD8CZ9a8kiRJ6tOgWwb/E/gX4OHa/3jg7sx8sPYvBras3VsCtwDU4ffU/I+kd43TK30lEXFYRCyMiIVLly5d02WSJElaa4wYDEbE8yJiw9r9uoj4WERs3cd4fwssyczLR6GcayQzT8rMOZk5Z/r06eNdHEmSpAmjn5bBE4EHIuLZwD8BvwJO62O85wH7RMRNlEu4LwI+AUyNiHVrnpnArbX7VmArgDp8U+CuZnrXOL3SJUmS1Kd+gsEHMzMp9+N9MjNPADYeaaTMfGdmzszMWZQHQL6TmQcAFwGvqNnmAefU7nNrP3X4d+p8zwX2r08bbwPMBn4MXAbMrk8nr1/ncW4fyyNJkqRq3ZGzcF9EvBN4HfCCiFgHWG8N5vmvwJkRcSzwU+Dkmn4y8PmIWAQsowR3ZOY1EXEWcC3wIPDmzHwIICLeAlwATAFOycxr1qBckiRJrdNPMPhq4LXAIZl5e0Q8CfjIqswkMy8GLq7dN1CeBO7O83vglT3G/wDwgSHSzwPOW5WySJIk6VH9BIOvBD6XmcsBMvPX9HfPoCRJkia4fu4ZnAFcFhFn1Td+xKALJUmSpLExYjCYme+hPLRxMnAQcH1EfDAinjLgskmSJGnA+vrT6fpU7+318yAwDTg7Ij48wLJJkiRpwEa8ZzAijgAOBO6kvFbuHZn5p/pU8fWUN4xIkiRpEurnAZLNgJdl5s3NxNyK3q4AABKJSURBVMx8uL5lRJIkSZPUiMFgZh4VEbtExIsy83MRMR3YKDNvzMzrxqCMGrBZR35zhf6bjtt7nEoiSZLGWj/vJj6K8kfR76xJ6wGnD7JQkiRJGhv9PEDyd8A+wG8BMvM39PE6OkmSJE18/QSDf6xPEydARGw42CJJkiRprPQTDJ4VEf8NTI2IQ4FvU54qliRJ0iTXz9PEHwVeDNwLPA14L/C9QRZKkiRJY6OfYPDkzPx7YAFARGwEnAfsNsiCSZIkafD6uUx8a0R8CiAipgHfwqeJJUmS1gr9vJv434D7I+LTlEDwo5n5uYGXTJIkSQPX8zJxRLys0Xsp8G/Aj4GMiJdl5lcHXThJkiQN1nD3DL60q/+nlD+cfinlb2YMBiVJkia5nsFgZh48lgWRJEnS2OvnARJJkiStpQwGJUmSWsxgUJIkqcVG/NPpiHgM8HJgVjN/Zh49uGJJkiRpLPTzBpJzgHuAy4E/DLY4kiRJGkv9BIMzM3PuwEsiSZKkMdfPPYM/jIhnDbwkkiRJGnP9tAzuAhwUETdSLhMHkJn5lwMtmSRJkgaun2Bwz4GXQpIkSeNixGAwM28ei4JIkiRp7A3sfwYj4rER8eOIuDIiromI99f0bSLi0ohYFBFfioj1a/pjav+iOnxWY1rvrOm/iIg9Gulza9qiiDhyUMsiSZK0thrkn07/AXhRZj4b2A6YGxE7Ax8CPp6ZTwWWA4fU/IcAy2v6x2s+ImJbYH/gGcBc4FMRMSUipgAnUC5jbwu8puaVJElSnwYWDGZxf+1dr34SeBFwdk2fD+xXu/et/dThu0VE1PQzM/MPmXkjsAjYsX4WZeYNmflH4MyaV5IkSX0a6OvoagveFcASYAHwK+DuzHywZlkMbFm7twRuAajD7wEe30zvGqdXuiRJkvo00GAwMx/KzO2AmZSWvKcPcn69RMRhEbEwIhYuXbp0PIogSZI0IQ00GOzIzLuBi4DnAlMjovMU80zg1tp9K7AVQB2+KXBXM71rnF7pQ83/pMyck5lzpk+fPirLJEmStDYY5NPE0yNiau3eAHgJcB0lKHxFzTaP8u5jgHNrP3X4dzIza/r+9WnjbYDZwI+By4DZ9enk9SkPmZw7qOWRJElaG/Xzp9Orawtgfn3qdx3grMz8RkRcC5wZEccCPwVOrvlPBj4fEYuAZZTgjsy8JiLOAq4FHgTenJkPAUTEW4ALgCnAKZl5zQCXR5Ikaa0zsGAwM68CnjNE+g2U+we7038PvLLHtD4AfGCI9POA89a4sJIkSS01JvcMSpIkaWIyGJQkSWoxg0FJkqQWMxiUJElqMYNBSZKkFjMYlCRJajGDQUmSpBYzGJQkSWoxg0FJkqQWMxiUJElqMYNBSZKkFjMYlCRJajGDQUmSpBYzGJQkSWoxg0FJkqQWMxiUJElqMYNBSZKkFjMYlCRJajGDQUmSpBYzGJQkSWoxg0FJkqQWMxiUJElqMYNBSZKkFjMYlCRJajGDQUmSpBYzGJQkSWoxg0FJkqQWMxiUJElqMYNBSZKkFhtYMBgRW0XERRFxbURcExFH1PTNImJBRFxfv6fV9IiI4yNiUURcFRHbN6Y1r+a/PiLmNdJ3iIir6zjHR0QMank0fmYd+c0VPpIkafQMsmXwQeCfMnNbYGfgzRGxLXAkcGFmzgYurP0AewKz6+cw4EQowSNwFLATsCNwVCeArHkObYw3d4DLI0mStNYZWDCYmbdl5k9q933AdcCWwL7A/JptPrBf7d4XOC2LS4CpEbEFsAewIDOXZeZyYAEwtw7bJDMvycwETmtMS5IkSX0Yk3sGI2IW8BzgUmBGZt5WB90OzKjdWwK3NEZbXNOGS188RLokSZL6NPBgMCI2Ar4CvC0z720Oqy16OQZlOCwiFkbEwqVLlw56dpIkSZPGQIPBiFiPEgh+ITO/WpPvqJd4qd9LavqtwFaN0WfWtOHSZw6RvpLMPCkz52TmnOnTp6/ZQkmSJK1FBvk0cQAnA9dl5scag84FOk8EzwPOaaQfWJ8q3hm4p15OvgDYPSKm1QdHdgcuqMPujYid67wObExLkiRJfVh3gNN+HvB64OqIuKKmvQs4DjgrIg4BbgZeVYedB+wFLAIeAA4GyMxlEXEMcFnNd3RmLqvdhwOnAhsA59ePJEmS+jSwYDAzvw/0+t+/3YbIn8Cbe0zrFOCUIdIXAs9cg2JKkiS1mm8gkSRJarFBXiaWBmKot5DcdNze41ASSZImP1sGJUmSWsxgUJIkqcUMBiVJklrMYFCSJKnFDAYlSZJazGBQkiSpxQwGJUmSWsxgUJIkqcUMBiVJklrMYFCSJKnFDAYlSZJazGBQkiSpxQwGJUmSWsxgUJIkqcUMBiVJklrMYFCSJKnFDAYlSZJazGBQkiSpxQwGJUmSWsxgUJIkqcUMBiVJklrMYFCSJKnFDAYlSZJazGBQkiSpxQwGJUmSWsxgUJIkqcUMBiVJklpsYMFgRJwSEUsi4meNtM0iYkFEXF+/p9X0iIjjI2JRRFwVEds3xplX818fEfMa6TtExNV1nOMjIga1LBPFrCO/ucJHkiRpTQ2yZfBUYG5X2pHAhZk5G7iw9gPsCcyun8OAE6EEj8BRwE7AjsBRnQCy5jm0MV73vCRJkjSCgQWDmfk9YFlX8r7A/No9H9ivkX5aFpcAUyNiC2APYEFmLsvM5cACYG4dtklmXpKZCZzWmJYkSZL6NNb3DM7IzNtq9+3AjNq9JXBLI9/imjZc+uIh0iVJkrQKxu0Bktqil2Mxr4g4LCIWRsTCpUuXjsUsJUmSJoWxDgbvqJd4qd9LavqtwFaNfDNr2nDpM4dIH1JmnpSZczJzzvTp09d4ISRJktYWYx0Mngt0ngieB5zTSD+wPlW8M3BPvZx8AbB7REyrD47sDlxQh90bETvXp4gPbExLkiRJfVp3UBOOiDOAXYHNI2Ix5ang44CzIuIQ4GbgVTX7ecBewCLgAeBggMxcFhHHAJfVfEdnZuehlMMpTyxvAJxfP5IkSVoFAwsGM/M1PQbtNkTeBN7cYzqnAKcMkb4QeOaalLEN/D9CSZI0HN9AIkmS1GIGg5IkSS1mMChJktRiA7tnUFod3uMoSdLYsmVQkiSpxQwGJUmSWsxgUJIkqcUMBiVJklrMYFCSJKnFDAYlSZJazL+W0UqG+nuXm47bexxKIkmSBs1gUKvFgFGSpLWDwaD64p9BS5K0dvKeQUmSpBazZVCjprv10MvGkiRNfAaDA2BQJEmSJgsvE0uSJLWYLYMTWBse2mjDMkqSNJHZMihJktRitgyuZWxpkyRJq8KWQUmSpBazZVAD41tKJEma+GwZlCRJajFbBrVW8L8dJUlaPQaDk5gPi/TWT90YMEqSZDA4YbQlsGvLckqSNFl4z6AkSVKL2TKo1vI+Q0mSDAalUWeQKUmaTCZ9MBgRc4FPAFOAz2bmceNcpJV4n9zk0M//IrouJUlrm0kdDEbEFOAE4CXAYuCyiDg3M68d35JpbTEawZ9/vi1JmsgmdTAI7AgsyswbACLiTGBfwGBQE9p4Xkr2MrYkqWmyB4NbArc0+hcDO41TWaTVNp6Xn/1PRklqt8keDPYlIg4DDqu990fELwY8y82BOwc8jzazfgdrpfqND41TSdZebsODZf0OlvU7WN31u/WgZzjZg8Fbga0a/TNr2goy8yTgpLEqVEQszMw5YzW/trF+B8v6HTzreLCs38GyfgdrPOp3sv/p9GXA7IjYJiLWB/YHzh3nMkmSJE0ak7plMDMfjIi3ABdQ/lrmlMy8ZpyLJUmSNGlM6mAQIDPPA84b73J0GbNL0i1l/Q6W9Tt41vFgWb+DZf0O1pjXb2TmWM9TkiRJE8Rkv2dQkiRJa8BgcBRFxNyI+EVELIqII8e7PJNJRJwSEUsi4meNtM0iYkFEXF+/p9X0iIjjaz1fFRHbN8aZV/NfHxHzxmNZJpqI2CoiLoqIayPimog4oqZbv6MkIh4bET+OiCtrHb+/pm8TEZfWuvxSfdCNiHhM7V9Uh89qTOudNf0XEbHH+CzRxBQRUyLipxHxjdpv/Y6SiLgpIq6OiCsiYmFN8xgxSiJiakScHRE/j4jrIuK5E6p+M9PPKHwoD7D8CngysD5wJbDteJdrsnyAFwDbAz9rpH0YOLJ2Hwl8qHbvBZwPBLAzcGlN3wy4oX5Pq93TxnvZxvsDbAFsX7s3Bn4JbGv9jmodB7BR7V4PuLTW3VnA/jX908CbavfhwKdr9/7Al2r3tvXY8Rhgm3pMmTLeyzdRPsDbgS8C36j91u/o1e1NwOZdaR4jRq9+5wP/ULvXB6ZOpPq1ZXD0PPJqvMz8I9B5NZ76kJnfA5Z1Je9L2YGo3/s10k/L4hJgakRsAewBLMjMZZm5HFgAzB186Se2zLwtM39Su+8DrqO8vcf6HSW1ru6vvevVTwIvAs6u6d113Kn7s4HdIiJq+pmZ+YfMvBFYRDm2tF5EzAT2Bj5b+wPrd9A8RoyCiNiU0uBxMkBm/jEz72YC1a/B4OgZ6tV4W45TWdYWMzLzttp9OzCjdveqa9fBCOrlsudQWq6s31FUL2FeASyhHKR/BdydmQ/WLM36eqQu6/B7gMdjHQ/nP4F/AR6u/Y/H+h1NCXwrIi6P8tYu8BgxWrYBlgKfq7c5fDYiNmQC1a/BoCaFLG3kPvq+BiJiI+ArwNsy897mMOt3zWXmQ5m5HeVNSDsCTx/nIq01IuJvgSWZefl4l2Uttktmbg/sCbw5Il7QHOgxYo2sS7kN6sTMfA7wW8pl4UeMd/0aDI6evl6Np1VyR20ap34vqem96tp10ENErEcJBL+QmV+tydbvANTLPxcBz6Vc3un8n2uzvh6pyzp8U+AurONengfsExE3UW7BeRHwCazfUZOZt9bvJcDXKD9oPEaMjsXA4sy8tPafTQkOJ0z9GgyOHl+NN/rOBTpPS80DzmmkH1ifuNoZuKc2tV8A7B4R0+pTWbvXtFar90qdDFyXmR9rDLJ+R0lETI+IqbV7A+AllHszLwJeUbN113Gn7l8BfKe2DJwL7F+fht0GmA38eGyWYuLKzHdm5szMnEU5tn4nMw/A+h0VEbFhRGzc6abs2z/DY8SoyMzbgVsi4mk1aTfgWiZS/Q7iqZm2fihPAP2Scq/Qu8e7PJPpA5wB3Ab8ifIr6hDKPT4XAtcD3wY2q3kDOKHW89XAnMZ0/p5yU/gi4ODxXq6J8AF2oVx+uAq4on72sn5HtY7/EvhpreOfAe+t6U+mBBuLgC8Dj6npj639i+rwJzem9e5a978A9hzvZZtoH2BXHn2a2PodnTp9MuUp6yuBazrnL48Ro1rH2wEL6zHi65SngSdM/foGEkmSpBbzMrEkSVKLGQxKkiS1mMGgJElSixkMSpIktZjBoCRJUosZDEpqrYi4f4ThUyPi8DEox9ER8eIR8uwaEX896LJIah+DQUnqbSow8GAwM9+bmd8eIduugMGgpFFnMChp0ouIWRFxXUR8JiKuiYhv1TeBdOfbJiJ+FBFXR8SxjfSNIuLCiPhJHbZvHXQc8JSIuCIiPjJMvu753B8RH69luTAiptf07SLikoi4KiK+Vt8iQEScGhGvqN03RcT7G/N4ekTMAt4I/L9aludHxCsj4mcRcWVEfG8061NSuxgMSlpbzAZOyMxnAHcDLx8izycoL4t/FuWNNx2/B/4uM7cH/gb4aH2N35HArzJzu8x8xzD5um0ILKxl+S5wVE0/DfjXzPxLypsFjhpiXIA76zxOBP45M28CPg18vJbl/4D3Antk5rOBfUasHUnqwWBQ0trixsy8onZfDswaIs/zKK8+BPh8Iz2AD0bEVZTXQm0JzBhi/H7zPQx8qXafDuwSEZsCUzPzuzV9PvCCHsvy1RGWA+AHwKkRcSgwpUceSRrRuuNdAEkaJX9odD8ErHSZuBrqHZwHANOBHTLzTxFxE+X9tqubr595DqezLA/R4zidmW+MiJ2AvYHLI2KHzLxrFecjSbYMSmqVHwD71+4DGumbAktqgPc3wNY1/T5g4z7ydVsHeEXtfi3w/cy8B1geEc+v6a+nXELu1wpliYinZOalmfleYCmw1SpMS5IeYcugpDY5AvhiRPwrcE4j/QvA/0TE1cBC4OcAmXlXRPwgIn4GnA98aKh8Q/gtsGNEvAdYAry6ps8DPh0RjwNuAA5ehbL/D3B2fWjlrZSHSWZTLl1fCFy5CtOSpEdE5qpevZAkDSci7s/Mjca7HJLUDy8TS5IktZgtg5IkSS1my6AkSVKLGQxKkiS1mMGgJElSixkMSpIktZjBoCRJUosZDEqSJLXY/wcLgyHrVNxg7wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plot_data_points_hist(d_all, binwidth=50)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Reformat stop_times table\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": 66, "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 key as distribution dictionnaries d_real/d_all\n", "stoptimes['key'] = stoptimes['trip_id'] + '__' + stoptimes['stop_id']\n", "stoptimes = stoptimes.set_index('key')\n", "\n", "# Subset stoptimes \n", "stoptimes = stoptimes[['trip_id','stop_id', 'route_desc', 'arrival_time', 'departure_time']]\n", "\n", "# Iterate over stoptimes to generate rounded hours\n", "list_hours = []\n", "size_stop_times = stoptimes.shape[0]\n", "\n", "for x in range(size_stop_times):\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", " # Print % progression \n", " if (x % 10000) == 0 :\n", " print('{}%'.format(round(100*x/size_stop_times,2)), end = ', ')\n", " \n", "# Add new column with rounded hour\n", "stoptimes['hour'] = list_hours\n", "stoptimes = stoptimes.drop(columns=['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": "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." ] }, { "cell_type": "code", "execution_count": 163, "metadata": {}, "outputs": [], "source": [ "with gzip.open(\"../data/stop_times_wHour.pkl\", \"rb\") as input_file:\n", " stoptimes = pickle.load(input_file)\n", "\n", "# load dictionnary df and stoptimes df\n", "distrib_df = pd.DataFrame(d_all).transpose()\n", "distrib_df_real = pd.DataFrame(d_real).transpose()\n", "\n", "stoptimes_df = pd.DataFrame(stoptimes)\n", "\n", "# Make join btw. stoptimes and all_distrib\n", "recovery_df = distrib_df.join(stoptimes_df)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We compute a script to make recovery tables. We use data from `real` distribution to make the recovery tables, as they are more likely to be precise." ] }, { "cell_type": "code", "execution_count": 164, "metadata": {}, "outputs": [], "source": [ "def make_recovery_tab(recovery_df, list_params):\n", " list_bins = [x for x in range(32)]\n", " # Aggregate counts on columns (delay in minute)\n", " recov_tab = recovery_df.groupby(list_params)[list_bins].apply(lambda x : x.astype(float).sum())\n", " recov_tab = recov_tab.astype('int')\n", " return recov_tab" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use our 4 parameters and make every possible combination of them. Note that `trip_id` does not need to be combine with `hour` and `route_desc` (transport type), as it already imply a specific time and transport type." ] }, { "cell_type": "code", "execution_count": 165, "metadata": {}, "outputs": [], "source": [ "# Combination of 3 parameters\n", "recov_stop_time_type = make_recovery_tab(recovery_df, ['stop_id','hour','route_desc'])\n", "# Combination of 2 parameters\n", "recov_time_type = make_recovery_tab(recovery_df, ['hour','route_desc'])\n", "recov_stop_type = make_recovery_tab(recovery_df, ['stop_id','route_desc'])\n", "recov_stop_time = make_recovery_tab(recovery_df, ['stop_id','hour'])\n", "# single parameter aggregations\n", "recov_trip_id = make_recovery_tab(recovery_df, ['trip_id'])\n", "recov_time = make_recovery_tab(recovery_df, ['hour'])\n", "recov_stop = make_recovery_tab(recovery_df, ['stop_id'])\n", "recov_type = make_recovery_tab(recovery_df, ['route_desc'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Overall aggregate distribution (used when everything else fails)" ] }, { "cell_type": "code", "execution_count": 166, "metadata": {}, "outputs": [], "source": [ "last_chance_distrib = np.array(recovery_df3.sum(axis=0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Evaluate precision of recovery tables - establishing hierarchy of _cascade_ \n", "\n", "We use only data with `geschaetz` or `real` status and many data points to evaluate perfomence of our recovery tables - namely probabilities coming from aggregated delay distribution from similar transfer." ] }, { "cell_type": "code", "execution_count": 167, "metadata": {}, "outputs": [], "source": [ "def cdf_distrib(distr_arr):\n", " # get total number of elements \n", " N_tot = np.sum(distr_arr)\n", " \n", " # make cumulative distribution probabilities\n", " cdf_distrib = np.empty((len(distr_arr)), dtype=float)\n", " save_x = 0\n", " for x in range(len(distrib)):\n", " cdf_distrib[x] = float(distr_arr[x])/float(N_tot) + float(save_x)/float(N_tot)\n", " save_x += distr_arr[x]\n", " \n", " return cdf_distrib\n", "\n", "def residual_in_cdf(cdf1, cdf2):\n", " res = N_tot = 0.0\n", " for i in range(len(cdf1)):\n", " res += abs(float(cdf1[i]) - float(cdf2[i]))\n", " N_tot += 1.0\n", " \n", " return res/N_tot" ] }, { "cell_type": "code", "execution_count": 168, "metadata": {}, "outputs": [], "source": [ "recovery_df_real = recovery_df.iloc[recovery_df.index.isin(distrib_df_real.index),:]" ] }, { "cell_type": "code", "execution_count": 169, "metadata": {}, "outputs": [], "source": [ "# declare empty lists\n", "res_recov_stop_time_type, res_recov_stop_time, res_recov_stop_type, res_recov_time_type, \\\n", "res_recov_trip_id, res_recov_type, res_recov_stop, res_recov_time, res_fail, \\\n", "= ([] for i in range(9))\n", "\n", "for index, row in recovery_df_real.iterrows():\n", " \n", " distrib_real = np.array(row[0:32])\n", " trip_id = row['trip_id']\n", " stop_id = str(row['stop_id'])\n", " transport_type = row['route_desc']\n", " hour = row['hour']\n", " key = str(trip_id) + '__' + str(stop_id)\n", " \n", " # control distrib with trip_id x stop_id \n", " cdf_real = cdf_distrib(distrib_real)\n", " sum_counts = np.sum(distrib_real)\n", " \n", " if sum_counts > 100 :\n", " # recov - stop_time_type\n", " if (stop_id, hour, transport_type) in recov_stop_time_type.index:\n", " cdf = cdf_distrib(np.array(recov_stop_time_type.loc[(stop_id, hour, transport_type)]))\n", " res_recov_stop_time_type.append(residual_in_cdf(cdf, cdf_real))\n", " \n", " # recov - time_type\n", " if (hour, transport_type) in recov_time_type.index:\n", " cdf = cdf_distrib(np.array(recov_time_type.loc[(hour, transport_type)]))\n", " res_recov_time_type.append(residual_in_cdf(cdf, cdf_real))\n", "\n", " # recov - stop_time\n", " if (stop_id, hour) in recov_stop_time.index:\n", " cdf = cdf_distrib(np.array(recov_stop_time.loc[(stop_id, hour)]))\n", " res_recov_stop_time.append(residual_in_cdf(cdf, cdf_real))\n", " \n", " # recov - stop_type\n", " if (stop_id, transport_type) in recov_stop_type.index:\n", " cdf = cdf_distrib(np.array(recov_stop_type.loc[(stop_id, transport_type)]))\n", " res_recov_stop_type.append(residual_in_cdf(cdf, cdf_real))\n", " \n", " # recov - type\n", " if (transport_type) in recov_type.index:\n", " cdf = cdf_distrib(np.array(recov_type.loc[(transport_type)]))\n", " res_recov_type.append(residual_in_cdf(cdf, cdf_real))\n", " \n", " # recov - trip_id \n", " if (trip_id) in recov_trip_id.index:\n", " cdf = cdf_distrib(np.array(recov_trip_id.loc[(trip_id)]))\n", " res_recov_trip_id.append(residual_in_cdf(cdf, cdf_real))\n", " \n", " # recov - stop\n", " if (stop_id) in recov_stop.index:\n", " cdf = cdf_distrib(np.array(recov_stop.loc[(stop_id)]))\n", " res_recov_stop.append(residual_in_cdf(cdf, cdf_real))\n", " \n", " # recov - time\n", " if (hour) in recov_time.index:\n", " cdf = cdf_distrib(np.array(recov_time.loc[(hour)]))\n", " res_recov_time.append(residual_in_cdf(cdf, cdf_real))\n", " \n", " # Evaluate compare with overall_average\n", " cdf_overall = cdf_distrib(last_chance_distrib)\n", " res_fail.append(residual_in_cdf(cdf_overall, cdf_real))\n", " " ] }, { "cell_type": "code", "execution_count": 217, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAADZCAYAAACNbSIWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3df7xVVZ3/8ddbIzUgSk2mHANNGykQ6qL90EycIktSm9HR1BSpyL5GU301+KalNjlqOlP4o6ycBM3CEdOUJqHkopa/QAUETXMStB8ziSV5jYzw8/1jr6ubwznnnnvv2ecX7+fjcR53n3XWXmt99jn3nHX2WmcvRQRmZmZm9bZNsxtgZmZmncmdDDMzMyuEOxlmZmZWCHcyzMzMrBDuZJiZmVkh3MkwMzOzQryk2Q2w5tl5551j9OjRdSnr2WefZejQoXUpq5k6JQ7onFg6JQ7onFg6JQ7onFiaHce99967LiJeVZruTsZWbPTo0SxbtqwuZS1ZsoSDDjqoLmU1U6fEAZ0TS6fEAZ0TS6fEAZ0TS7PjkLS2XLqHS8zMzKwQ7mSYmZlZITxcYluF8WcvYv2GjbVlvvmHhbVj+JhZPPPQeYWVv4U6xzJihyGsOHNyXcs0s87lToZtFdZv2Mia8w7tM1/R45rj5s6qqR31UEQso2cV1wEzs87j4ZIGkTRV0muqPP5FSe8aQLknSzqhTPpoSav6W95gSWp0lWZWEP8/22BtlWcylP3nKCKeb2C1U4FVwG/KtGfbiPjCQAqNiMsG2S4zM7NCbDVnMtI3+4clXUn2Yf95SUslrZR0di7fCSlthaSrcvsuTum3SHqtpBGS1kraJuUZKukJSUPK1H0kMBG4WtJySTtIWiPpfEn3AUdJmpPykR77sqQHJN0jac8qcZ0l6dS03ZXavQI4pX5Hz8zMrP+2tjMZewEnAi8HjgT2AwTcKOlA4CngDODtEbFO0o5pv4uBuRExV9I04KKIOELScuCdQDcwBVgYEVvMLoyI+ZI+AZwaEcvghdOQT0XEm9P9Q0p2Wx8R49JQyFdT+X25AvhERNwm6YJyGSRNB6YDjBw5kiVLltRQbN96enpeKKtVx+1riTUfRzPbUQ9FxdKU57fAybgN12axlHsNNeL/pFE6JZaWjSMitoobMBp4LG1fCKwBlqfbo8CHgRnAOWX2XQcMSdtDgHVp+1jgsrR9PfDuKvUvASbm7q8BRuXuzwGOzD22R66+p6qUexZwKvAK4PFc+j7AqmrHpKurK+qlu7s7Iqu4bmXW06iZC2rK1xtHUcbOGVto+XlFxFLrcaynop+TRmq3WCr9P7dbHNV0SizNjgNYFmU+Z7aa4ZLk2fRXwLkRMSHd9oyI/xhAeTcCh6QzHl3A4gG2p5yosG1mZtYWtrZORq+FwDRJwwAk7SppF7JOwlGSdkrpvcMldwDHpO3jgNsBIqIHWArMBhZExKYqdT4DDO9HG4/O/b2zr8wR8TTwtKQDcu1suKxDa2adwP/PNlhb25wMACJikaQxwJ1pbkQPcHxErJZ0DnCrpE3A/WS/CpkBXCHpNOBJ4KRccdcA1wIH9VHtHOAySRuAt9XQzFdKWgk8B3ywxtBOAr4tKYBFNe6z1ah5LkGhF+Nq8JyGAi7GZWZWq62mkxERa4Cxufuzyc5AlOabC8wtSVsLHFyh3Plkwy991X8dcF0uaXTJ41NLdrkgImbWUO5Zue17gfG5hz/b1/5bi1ovgFX8IkONuRAXNH/BJDOzrXW4xMzMzAq21ZzJaBRJlwL7lyTPjograi0jIkaXKfd04KiS5Gsj4px+N9LMzKwB3Mmos4go5CJYqTPhDoWZmbUND5eYmZlZIdzJMDMzs0K4k2FmZmaFcCfDzMzMCuGJn2ZljD97Ees3bLHWXeGGj5nFMw+dV78Cm7QY14gdhrDizMlNqdvMWoc7GWZlrN+wseYLeNXTuLmz6lZvMy/G1aor8ZpZY3XMcImkgyS9vc5lTpD0vtz9wyTNqmcdJfV9StLLiirfzMyskQbdyVCmFTorBwF17WQAE4AXOhkRcWNE1PFc9hY+BbRlJ2PSpEnNboKZNUBa78msJgPqHEgaLelhSVcCq4DPS1oqaaWks3P5TkhpKyRdldt3cUq/RdJrJY2QtLa3syJpqKQnJJVdjUnSJyU9mMqYJ2k0cDLwaUnLJb2jXD1p3zmSLpO0TNIjkqZUqOOlwBeBo1OZR0uaKumSXDlfl3SXpF+mMynflvSQpDm5ciZLulPSfZKu7V35tVxMwGuAbkndkqZJ+mru8Y9K+kqK6+eSrk51ze89+yGpS9Ktku6VtFDSq2t6Qs3MzAqggSzlmz7Uf0l25uDlwJHAx8gWCrsR+DLwFHA98PaIWCdpx4j4vaSbgPkRMVfSNOCwiDhC0g+Ar0ZEt6SjgXdHxEcq1P8bYPeIeE7SKyLiaUlnAT0RcWHKU6meOcDfkJ2heB3QDewZEX8uU89UYGJEfKL0fipne7IVUg8DriK7nPhqsuXfPwz8Cvg+8N6IeFbSTGC7iPhihbjWpPLXpc7ICmDviNgo6Y50jJ8BHgMOiIifSfo28CDZYm+3AodHxJPpGL4nIqaV1DEdmA4wcuTIrnnz5pVrSr9NmjSJUTMX1KWsVjHnkKENr3PG2hlcPOriupTV09PDsGFl+7SFm3rzs02p14q39vwp3HTTTU17bdVbM/9P6qnZcUyaNOneiJi4xQMR0e8b2Qqij6XtC4E1wPJ0e5TsA3YGcE6ZfdcBQ9L2EGBd2j4WuCxtX0/WyahU/83AfOB4YFhKOws4tYZ65gDTcvluAyZUqGcqcEm5+6mc49L2HsAvcvmuBI4ApqR29B6bB4H/qBLXGmDn3P1vAR8A9gaW5o7947k8BwM3kK0w+8dcXQ8Ai6o9j11dXVEv2Uup/XV3d0dExKiZC5pS/9g5Y+tWVm8szVDP49fMOOqtE2IBOiKOXp0SS7PjAJZFmc+Zwfy6pPerioBzI+Ib+QclzehneTcC/yppR6ALWFwl76HAgcD7gdMljetnXaWnb/p/OifzXPr7fG679/5LgE3AjyPigwMs/3Lgc8DPgfwCa+XaL2B1RLxtgHWZmZnVVT0mbC4EpvXONZC0q6RdyDoJR0naKaXvmPLfARyTto8DbgeIiB6yYYbZwIKI2FSusjRvY7eI6AZmAiOAYWTDCMNzWcvWkxwlaRtJryM7C/FwhdhKy+yvu4D9Je2Z2j5U0uur5N+svoi4G9iN7CzP93L5XiuptzNxLPBTshhe1ZsuaYikNw6i7f3S3d3dqKrMrIliAEPstvUadCcjIhYB3wXulPQA2TDG8IhYTbZq6K2SVgD/nnaZAZwkaSXwIeCfc8VdQzYEck2VKrcFvpPquh+4KCKeBm4CPtA78bOPeh4H7gF+BJwcZeZjJN3AG3onftZyPPIi4kmyIZbvpXbcSTb0Uck3gZsl5T+x/xP4WUT8IZf2MHCKpIeAVwJfj4i/kM2NOT8d7+XU/9c2ZmZmNRvQcElErCGbA9B7fzbZGYjSfHOBuSVpa8nmEZQrdz7Zaf9qdW8EDiiT/giwT0ly2XqAn0TEydXqSWX+Hti3JHlOemxqLt8aNj8e+ccWlymjUn0XA6Wz/g4AvlKS9teIOL7M/svJhpGsDppxQanhY+pcbxOv+Glm5it+tihJryA727IiIm5pdnu2Ns242memfvU284qfZmbQ4p0MSZeS/Sw0b3ZEXFEufy3yZxly9bwHOL8k+bGI+MBA66lG0vXA7iXJMyNiYe+dNAS0xfyN0rMmZmZmraqlOxkRcUqD6llINoG1IYrqvJiZmbWSVrgcuJmZmXUgdzLMzMysEO5kmJmZWSHcyTAzM7NCuJNhZmZmhWjpX5eYNdP4sxexfsPGhtY5fMwsnnnovPoV2KSLcUF2Qa4VZ05uWv1m1nzuZJhVsH7DxoZflGvc3Fl1q7PZF+NqxhVTzay1eLhkACQdJKmu64JImiDpfbn7h0maVc86ijRp0qRmN8HMGsT/71artupkKNMKbT6I+i8+NgF4oZMRETdGRB3Pm5uZmTWWWn3ZXkmjya7GeTfQRbYq6RRgO+D6iDgz5TsBOBUIYGVEfCjt+21gZ+BJ4CRgPbAS2D0inpc0FPg5sEdafK20/k8CJwN/BR4EZpEt4b4plTkDeKK0noh4XNIc4M/ARODlwGciYkGZOl4KPArsAPwaODdtT4yIT6RyNgBvAnYBpgEnAG8D7u69VLqkycDZ6dj8d2pHT0ld04HpACNHjuyaN29exWPfH5MmTWLUzC1Ca3tzDhna0PpmrJ3BxaNK18gbmJ6eHoYNG1aXsgZi6s3PNq1uK9ba86fQ3d3dd8Y20Oz/k3ppdhyTJk26NyImbvFARLT0DRgNPA+8FZhMthy6yM7CLCBbdfSNwCPAzmmfHdPfm4AT0/Y04Ia0/QNgUto+Gri8Sv2/AbZL269If88CTs3lqVTPHODm1Na9gF8B21eoZypwSbn7qZx5Ke7DgT8C41K595KdBdkZuA0YmvaZCXyh2rHt6uqKesleSu2vu7v7he1RMxc0vP6xc8bWrax8LM1Qr+PX7DjqqVNi6ZT/94jOeU6aHQewLMp8zrTC0EMt1kbEXWSdjMnA/cB9wN5kH94HA9dGxDp4YYl2yL7pfzdtX8WLS8RfQ9a5ADgm3a9kJXC1pOPJzmaUU6kegP+MiOcj4hfAL1ObB+Km9EQ+APxvRDwQEc8Dq8k6Ym8F3gD8TNJy4ERg1ADrMjMzG7R2+XVJ73lXAedGxDfyD0qa0c/ybgT+VdKOZEMwi6vkPZTsbMn7gdMljetnXaXjUQMdn3ou/X0+t917/yVkwzc/jogPDrB8MzOzumqXMxm9FgLTJA0DkLSrpF3IOglHSdoppe+Y8t9BdqYC4DjgdoDI5iksBWYDCyJiU7nK0iTT3SKim2z4YQQwDHgGGJ7LWrae5ChJ20h6HbAH8HCF2ErL7K+7gP0l7ZnaPlTSFkvFF6VTxmfNrG/+f7datcuZDAAiYpGkMcCdkgB6gOMjYrWkc4BbJW0iG06ZSjYp8wpJp/HixM9e1wDXkv1SpJJtge9IGkF2FuWiiHha0k3AfEmHpzqq1fM4cA/ZxM+TI+LPFerqBmaloY5zazogORHxpKSpwPckbZeSzyCbq2ID1OhrPQwfU+c6m3wxLjPburV8JyMi1gBjc/dnk52BKM03F5hbkraWbL5GuXLnk3UcqtW9kc3nV/SmPwLsU5Jcth7gJxFxcrV6Upm/B/YtSZ6THpuay7eGzY9H/rHFZcqwAWr0hbgy9auz2RfjMjNrt+ESMzMzaxMtfyajUSRdCuxfkjw7Iq4YaJn5swy5et4DnF+S/FhEfGCg9ZiZmbUidzKSiDilQfUsJJvAamZm1tE8XGJmZmaFcCfDzMzMCuFOhpmZmRXCnQwzMzMrhCd+mhVk/NmLWL9hi4V9G6vkYlzDx8zimYfOa1JjBqHMRcVG7DCEFWdObkJjzKxW7mQMgKSDgL9ExB11LHMC8JqI+K96lWnNtX7DxiZd0CtT7mJc4+bOamqbBqLSRcUafTVWM+u/thouUaYV2nwQ8PY6lzkBeF+dy2yYSZMmNbsJZrYVS0tNWItphQ/sqiSNlvSwpCuBVcDnJS2VtFLS2bl8J6S0FZKuyu27OKXfIum1kkZIWtvbWUkLiT0hqexCC5I+KenBVMY8SaOBk4FPS1ou6R3l6kn7zpF0maRlkh6RNKVCHS8Fvggcnco8WtIvJL0qPb6NpEclvapSmZK2lXRB7th8rD7PgJmZ2cC0y3DJXsCJZIuMHQnsR7buyI2SDgSeIlsM7O0RsS63CuvFwNyImCtpGtkCZ0ekRcjeSbYo2RRgYVqnpJxZwO4R8ZykV6QF0i4DeiLiQoC0YNpm9QBHpP1Hp/a+DuiWtGfpImkR8RdJXwAmRsQnUpl7k63o+lXgXcCKtAha2TKBE4D1EbFvWiDtZ5IWRcRj+bokTQemA4wcOZIlS5b0efBr1TGnr+u4qFg9j29/9fT0lK2/mW0aiEpxQJu+5pq4aF1dtWAcA3ltV3t9tZOWjSMiWvpG9oH6WNq+EFgDLE+3R4EPk62Cek6ZfdcBQ9L2EGBd2j4WuCxtXw+8u0r9NwPzgeOBYSntLODUGuqZA0zL5bsNmFChnqnAJbn7uwH3pe15wJRqZaY2PpI7No8Bk6sd266urqiX7KXU/rq7u+tW1qiZC+pW1kCUi2XsnLGNb8ggVXpOmn18B6Ker69masU4Bvoe1IqxDESz4wCWRZnPmZYfLkmeTX8FnBsRE9Jtz4j4jwGUdyNwSDrj0QUsrpL3UOBS4M3AUkn9PfsTfdwvv1PEE8D/SjqY7KzFj/ooU8CM3LHZPSIW9bOtZmZmddMunYxeC4FpkoYBSNpV0i5knYSjJO2U0nuHS+4AjknbxwG3A0RED7CUbMn4BRGxqVxlad7GbhHRDcwERgDDgGeA4bmsZetJjkpzKl4H7AE8XCG20jIBLge+A1xb0sZyZS4EPt47t0TS6yUNrVCXmZlZ4dplTgYAEbFI0hjgzjQ3oQc4PiJWSzoHuFXSJuB+suGHGcAVkk4DngROyhV3DXAt2S9FKtkW+I6kEWRnCi6KbE7GTcB8SYenOqrV8zhwD9l8kpOjZD5GTjcwK80XOTciriE743JFuuVtUaaky8mGlu5TdnCe5MV5IYXr7u5uVFVmZlvIzthbq2n5TkZErAHG5u7PJjsDUZpvLjC3JG0tcHCFcueTdRyq1b0ROKBM+iPAPiXJZesBfhIRJ1erJ5X5e2DfkuTxZBM+f95XmRHxPPC5dLMW0fSJiVtcjKsF2jQQFS7GZWatreU7GVsrSbOAj5MNv1gbavZFr8pfxKq9LsQFlS/GZWatz52MRNKlwP4lybMjonSoomYRMbVMPe8Bzi9JfiwiPlCy73nAFtd/LlemmZlZK3InI4mIUxpUz0KySZpmZmYdrd1+XWJmZmZtwp0MMzMzK4Q7GWZmZlYIdzLMzMysEO5kmJmZWSH86xKzAo0/exHrN1Ra4LcBylzEaviYWTzz0Ba/jm5tFVb8HLHDEFacObnBjTGzWrmT0QYkHQT8JSLuaHZbrH/Wb9jYtItyVbqI1bi5s5p+obD+qHYxrra8eqnZVsTDJVUo0wrH6CDg7c1uxGCktWbMzArh95jW1AofoC1F0mhJD0u6ElgFfF7SUkkrJZ2dy3dCSlsh6arcvotT+i2SXitphKS1vZ0VSUMlPdG7WmqZ+j8p6cFUxjxJo4GTgU9LWi7pHeXqSfvOkXSZpGWSHpE0pdijZWZmVpmHS8rbCziRbJXTI4H9yBZTu1HSgcBTwBnA2yNiXW5p+YuBuRExV9I0slVbj0grq76TbKXVKcDCtPhaObOA3SPiOUmvSKu+Xgb0RMSFAGkV2M3q4cUVV0en9r4O6Ja0Z5WVX83MzArjTkZ5ayPiLkkXApPJlo4HGEbWARkPXBsR6+CFFVQB3gb8Q9q+Cvhy2r4GOJqsk3EM8LUqda8ErpZ0A3BDhTyV6gH4z7Qi6y8k/RLYG1je+6Ck6cB0gJEjR7JkyZIqTaldT09Pn2W1zfh5hUmGA1WvY9xf1Z6TZrVpIPp6bbXN66pXnV9fTdOCcQzkdV3Le1c7aNk4IsK33I3sTMCqtP1vwMfK5JkBnFMmfR0wJG0PAdal7WHAGmBH4HFg2yr1bwtMAv4deIisI3gWcGoN9cwBTsrluw0YX6murq6uqJfu7u6qj2cvtdbXVxz9NWrmgrqW1x+VYhk7Z2xjGzJI1Z6TZh7fgaj366tZWjGOgb7HtGIsA9HsOIBlUeZzxnMyqlsITJM0DEDSrpJ2ARYDR0naKaX3DpfcQXamArIl2m8HiIgeYCkwG1gQEZvKVZbmbewWEd3ATGAEWQflGWB4LmvZepKjJG0j6XXAHsDDA4zdzMxsUDxcUkVELJI0BrgzzVzuAY6PiNWSzgFulbSJbDhlKtkZjisknQY8CZyUK+4a4FqyX4pUsi3wHUkjyOaAXBTZnIybgPmSDk91VKvnceAesvkkJ0eLzMfIOrpmZsXwe0xrciejRESsAcbm7s8mOwNRmm8uMLckbS1wcIVy55N1HKrVvRE4oEz6I8A+Jcll6wF+EhEnV6vHGqupcwbKXoyrc+YxjNih7I+0zKxFuJNhVqBmXvSq8kWs2udCXFD9Ylxm1trcyWgSSZcC+5ckz46IKwZaZkRMHVSjzMzM6sidjCaJiFOa3QYzM7Mi+dclZmZmVgh3MszMzKwQ7mSYmZlZIdzJMDMzs0K4k2FmZmaF8K9LzNrI+LMXsX5DpQV8y2jAIlbDx8zimYfOK7aSJi7GNWKHIaw4c3LT6jdrZ+5kmLWR9Rs21nyBr0ZdxGrc3FmFXnSs2Rfjaruro5q1EA+XVCFpqqTX5O5fLukNBdU1WtKxRZRtZmbWDG3TyVCm0e2dCrzQyYiIj0TEgwXVNRrY6joZaeE5M7OO4Pe0zbV0JyN9u39Y0pXAKuDzkpZKWinp7Fy+E1LaCklX5fZdnNJvkfRaSSMkre3trEgaKukJSVussiTpSGAicLWk5ZJ2kLRE0sT0eI+kCyStlvQTSfulx38p6bCUZ9uUp7fNH6sS7nnAO1Jdn5Z0m6QJufb8VNJ4SWdJukrSnZJ+IemjuTynlTs+ZmZmzdAOczL2Ak4kW7r8SGA/stVMb5R0IPAUcAbw9ohYJ2nHtN/FwNyImCtpGtmy6UdIWg68E+gGpgAL0+qnm4mI+ZI+AZwaEctgix7qUGBxRJwm6XrgS8C7gTeQrc56I/BhYH1E7CtpO+BnkhZFxGNl4pyV6pqS6vo92ZmUT0l6PbB9RKyQ9AGyFVnfmtpwv6Qfkq0cu1fp8YmI2/KVSJoOTAcYOXIkS5YsqXLoa9fT0zPgslpuzLuJkwxrUetxHsxz0l9F1tPIOCqp62u0xV9fNeuUOKDusTTj9doK/ydlRUTL3siGEB5L2xcCa4Dl6fYo2Yf4DOCcMvuuA4ak7SHAurR9LHBZ2r4eeHeV+pcAE8vdB54DlLa/CJyetrcBnk7b84FHcm1+DJhcoa6DgAW5+y9LMQ4hO8vxiZR+FvDFXL4rgSMqHZ9qx7erqyvqpbu7e0D7ZS/B1jHQOBpl1MwFNedtVCxj54wttPxmPyf9OeZ9aXYs9dIpcUTUP5Zmvac1+zkBlkWZz5l2OJPxbPor4NyI+Eb+QUkz+lnejcC/pjMeXcDiAbZrYzqwAM+TdTqIiOcl9R5XATMiYmF/C4+IP0n6MXA48E+prS88XJqdCsfHzMysWVp6TkaJhcA0ScMAJO0qaReyTsJRknZK6b3DJXcAx6Tt44DbASKiB1gKzCY7c7CpSp3PAMMH2eaP9875kPR6SUP7UdflwEXA0oj4Qy79cEnbp5gPIoun0vFpaS/208zM2p/f0zbXDmcyAIiIRZLGAHemuRE9wPERsVrSOcCtkjYB95PNZZgBXCHpNOBJ4KRccdcA15J9QFczB7hM0gbgbQNo9uVkQz73KWv0k2RDG+WsBDZJWgHMiYivRMS9kv4IXFEmbzewM/AvEfEb4Dfljg/wuwG028zMbNBaupMREWvIJjT23p9NdgaiNN9cssmW+bS1wMEVyp1PNrzQV/3XAdflkg7KPTYst31WyX7D0t/ngc+lW191bSxtb7pGxzbAopLsKyPihDJllD0+1ln6NQmxIVf8bMDk3SZf8dPMBqalOxlbM0knAOcAn0mdFbN+XVmzcVfKLO5qn9D8K36a2cC5kwFIuhTYvyR5dkSUDlPUo65xwFUlyc9FxFvyCRFxJdkvRyhJP6vebTIzMyuCOxlARJzSwLoeACb0mdHMzKzNtdOvS8zMzKyNuJNhZmZmhXAnw8zMzArhToaZmZkVwp0MMzMzK4R/XdJCJE0FFqUreJrVxfizF7F+wxYLDbeU4WNm8cxD51XO0IYrfo7YYQgrzpzc7GaYNZU7GRWky4CrwRfCmgqsAtzJsLpZv2Fjvy7i1Qzj5s6q2MZ2vRhX4VdBNWsDHi7JkTRa0sOSriT7sP+8pKWSVko6O5fvhJS2QtJVuX0Xp/RbJL1W0ghJayVtk/IMlfRE74JpJXUfCUwErpa0XNKhkm7IPf5uSden7R5JX5G0OtX1qpT+Okk3S7pX0u2S9i7yeBUtrcFiZtbxOvX9zp2MLe0FfA34NLArsB/ZxbO6JB0o6Y3AGcDBETEe+Oe038XA3IjYB7gauCgi1gPLgXemPFOAhWmdks2k9VSWAcdFxATgv4C9ezsQZAu8fTttDwWWRcQbgVuBM1P6N8mWlu8CTk1xmJmZNYW8LO2LJI0GuiNid0kXAkcCT6eHhwHnAi8D/iYiTi/Zdx3w6ojYmM5U/DYidpZ0LHBgRJyczkR8LSJ+XKH+JcCpEbEs3T8d+BPZKqz3A3tFxF/TarPbpe09gO8DB5Ct8vpwrsjtImJMSR3TgekAI0eO7Jo3b17/D1QZPT09DBs2rO+M/TBp0iRGzVxQ1zK3VnMOGdrsJlQ1Y+0MLh51cdnHinhtNcLUm59tdhOsjaw9fwrd3d0D3r/Z/yeTJk26NyImbvFARPiWbmTLsq9K2/8GfKxMnhnAOWXS1wFD0vYQYF3aHgasAXYEHge2rVL/EmBi7v5rgHuBjwNfzqVvAl6Stvcg64C8nKxjU3O8XV1dUS/d3d11K6tX9vJsrCLiaJbeWEbNXNDchtRg7JyxFR9r1+ek3HFv11hKdUocEa0Ty2Df75odB9nZ9S0+ZzxcUtlCYJqkYQCSdpW0C7AYOErSTil9x5T/DuCYtH0ccDtARPQAS8mWYF8QEZuq1PkMMLz3TmS/MvkN2fBMfrG2bcjOsgAcC/w0Iv4IPCbpqNQuSRo/kMDNzMzqwb8uqSAiFkkaA9yZJuT0AMdHxGpJ5wC3pmGL+8l+FTIDuELSaWTDFiflirsGuBY4qI9q5wCXSdoAvC0iNpDN73hVRDyUy/cssJ+kM4DfAUen9OOAr6f0IcA8YMUAwm8J4aE8M9tKdOr7nTsZORGxBhibuz+b7AxEab65wNyStLXAwRXKnQ/0OXU4Iq4DritJPgD4Vpm8nymT9hhwSNx7tzcAABUnSURBVF/1mJmZNYI7GS1M0r1kZy3+b7PbYu2t1a/ZMHxMH21s04txmW3t3MloAkmXAvuXJM+OiPy8CyL7KeoWIqL9ptpb07T6hbgyldvYrhfjMjN3MpoiIk5pdhvMzMyK5l+XmJmZWSHcyTAzM7NCuJNhZmZmhXAnw8zMzArhToaZmZkVwp0MMzMzK4R/wmrWIcafvYj1GzZunthmF7EaPmYWzzx03pYPtFkc1Yy4dRErzpzc7GaYNYQ7GWYdYv2GjZtdeKsdL2I1bu6sLS4e1o5xVLJkyRIvAW9bFQ+XtAhJUyW9Jnf/cklvaGabGiUtQGdmZgVr9PutOxllpGXSG31spgIvdDIi4iMR8WCD22BmZlY3Hi5JJI0GFgJ3A13Af0qaAmwHXB8RZ6Z8JwCnAgGsjIgPpX2/DezMi8u8rwdWArtHxPOShgI/B/aIiM0GziUdCUwEru5d5h34EXBqRCyT1AN8HXgf8Fvgc8CXgdcCn4qIGyVtC5xHtpz8dsClEfGNMnFOB6YDjBw5kiVLlgzquPXq6ekZVFkttYBXG4//55+DwT4nzVLa5naNo5yenh5ArfV6H6g2/j/ZQqfEUmMcDf1/igjfIgBGA88DbwUmA98kW559G2ABcCDwRuARYOe0z47p703AiWl7GnBD2v4BMCltHw1cXqX+JcDEcvfJOjTvTdvXA4uAIcB4YHlKnw6ckba3A5aRdXAqxtzV1RX10t3dPeB9s5dhaxhMHM02auaCze63Yyxj54zdIq0d46iku7t7i+epHXXac9IJao2jqPdbYFmU+ZzxmYzNrY2IuyRdSNbRuD+lDwP2IvtQvzYi1gFExO/T428D/iFtX0V2lgHgGrLORTdwDPC1AbbrL8DNafsB4LmI2CjpAbLOEam9+6SzIgAjUpsfG2CdZmZmg+JOxuZ6p30LODdKhhskzehneTcC/yppR7IhmMUDbNfG1FOE7GzLcwCRDcP0PocCZkTEwgHWYWZmVlee+FneQmCapGEAknaVtAtZJ+EoSTul9B1T/jvIzlQAHAfcDhARPcBSYDawICI2VanzGWD4INv8cUlDUtten+aBtLwX+09mZlakRr/f+kxGGRGxSNIY4M70c58e4PiIWC3pHOBWSZvIhlOmAjOAKySdxosTP3tdA1xLNiGzmjnAZbmJn/11OdnQyX3KGv0kcMQAyrE2tsWEwjab0DZ8TIVJwG0WRzUjdhjS7CaYNYw7GUlErAHG5u7PJjsDUZpvLjC3JG0tcHCFcueTDWX0Vf91wHW5pINyjw3LbZ9Vst+w9Pd5sl+dfK6vuqwzdcZFrA7dIqU94yivk2Ixq4WHS8zMzKwQPpPRYJIuBfYvSZ4dEVc0oz1mZmZFcSejwSLilGa3wczMrBE8XGJmZmaFcCfDzMzMCuFOhpmZmRXCnQwzMzMrhCd+mnWg8WcvYv2GjW15EavhY2bxzEPnbZ7YhnFUVBLLiB2GsOLMyU1qjFmxqnYy0hLmCyJibLV8AyXpILLlzKdIOgx4Q0Sc18duZtaH9Rs2MueQoW154adxc2dtdmGxTrqAVblYOmLZd7MKWuZMRkTcSLagWEuStG0fa480tP5a29PsdtdCktcvMTNrsEa899YyJ+Mlkq6W9JCk+ZJelhr3BUlLJa2S9M20XgaSPinpQUkrJc1LaUMlfVvSPZLul3R4aSWSpkq6JG3PkXSRpDsk/TK3fDmSTkv1rpR0drkGS/q6pGWSVvfmkXSIpGtzeQ6StCBtT5Z0p6T7JF2bWxhtjaTzJd1HtjDaR1PdKyRdlzsWr5N0l6QHJH1JUk8/21tr/aX3P5jqXCXp/Fx5PZL+TdIKBrYOipmZ2aDV0sn4O+BrETEG+CPwf1L6JRGxbxpK2QGYktJnAW+KiH2Ak1Pa6cDiiNgPmARcUMMKoa8GDkjlngfZhzGwF7AfMAHoknRgmX1Pj4iJwD7AOyXtA/wEeEuu3qOBeZJ2Bs4A3hURbwaWAZ/JlfVURLw5IuYB308xjwceAj6c8swmu2rnOOBXvTvW0t5+1v/CfeA24HyyNVMmAPtK6l0QbShwd0SMj4iflj+8ZmZmxapluOSJiPhZ2v4O8EngQmCSpM8CLwN2BFYDNwErgasl3QDckPabDBwm6dR0f3vgtX3Ue0Na9OtBSSNz5UwmW/0UYBjZh/htJfv+k6TpKb5Xk831WCnpZuD9kuaTrcT0WeCdwBuAn6WTMS8F7syVdU1ue6ykLwGvSHUvTOlv48UVT7+bjk+t7X1rP+rP398XWBIRTwJIuho4kOyYb2LzxdZekI7LdICRI0eyZMmSctn6raenZ1BltdS4dIdMMhzsc9JM+Xa3cxylKsXSUq//WnXI/wnQObEMII6i/7dq6WSUDtiEpO2BrwETI+IJSWeRdRwg+/A+EHg/cLqkcWSrkP5jRDycLyjXeSjnuXzW3N9zI+IblXaStDtwKrBvRPxB0pxc2+YBnwB+DyyLiGfSMM+PI+KDFYp8Nrc9BzgiIlZImkrfy7f32d6Up9b6y90v58+V5mFExDeBbwJMnDgx6jWhbrCT80pXEG2WTplkOHrWDxk2bFh7xjKXzdrdKc8JVIjl5h+2zOu/Vh3/nLShgcSh8yk89lqGS14rqXdc/1jgp7z4ob0uzR84EkDSNsBuEdENzARG8OI3/hm5eRtvGmB7FwLTcnMWdpW0S0mel5N9EK9PnZj35h67FXgz8FGyDgfAXcD+kvZMZQ6V9PoK9Q8HfitpCHBcLv0u4B/T9jH9bG9/6s+7h2woaGdJ2wIfTPGZmZm1hFrOZDwMnCLp28CDwNcj4k+SvgWsAv4HWJrybgt8R9IIsm/oF0XE05L+BfgqsDJ1RB7jxTkcNYuIRZLGAHem/koPcDzwu1yeFZLuB34OPAH8LPfYpjTZcypwYkp7Mp2V+J6k7VLWM4BHyjTh88DdwJPp7/CU/qkU9+nAzcD6frS3P/Xnj8VvJc0CusmO9Q8j4gfV9mlV/mWJmVnjNeK9V36DH7z0K5MNERGSjgE+GBFb/IKm1UycODGWLVtWl7K25lOOragtx/iTshfj6mDteDGuTvk/gc6JpdlxSLo3/eBiMy1znYw21wVckoaDngamNbk9tpVbc96hTX/TGbjN5ye0bxxb6qRYzGrhTkYdRMTtwPhmt8PMzKyVeIE0MzMzK4Q7GWZmZlYIT/zcikl6Elhbp+J2BtbVqaxm6pQ4oHNi6ZQ4oHNi6ZQ4oHNiaXYcoyLiVaWJ7mRYXUhaVm5mcbvplDigc2LplDigc2LplDigc2Jp1Tg8XGJmZmaFcCfDzMzMCuFOhtXLN5vdgDrplDigc2LplDigc2LplDigc2JpyTg8J8PMzMwK4TMZZmZmVgh3MmwLkg6R9LCkR9MibKWPbyfpmvT43ZJG5x77fyn9YUnvqbXMohQUyxpJD0haLqk+i78UFIeknSR1S+qRdEnJPl0pjkclXdS7SnKbxrIklbk83UpXO26lON4t6d507O+VdHBun4Y/JwXF0fDnY5Cx7Jdr6wpJH6i1zDaKo+HvW0C2CptvvvXeyFbS/W9gD+ClwArgDSV5/g9wWdo+Brgmbb8h5d8O2D2Vs20tZbZLLOmxNcDObfKcDAUOAE4GLinZ5x7grWSr+P4IeG8bx7IEmNgmz8mbgNek7bHAr5v1nBQYR0OfjzrE8jLgJWn71WQrZb+kljLbIY50fw0NfN/qvflMhpXaD3g0In4ZEX8B5gGlK8oeDsxN2/OBv0/fuA4H5kXEcxHxGPBoKq+WMtsllmYYcBwR8WxE/BT4cz6zpFcDL4+IuyJ7B7oSOKLQKDJ1j6VJBhPH/RHxm5S+GtghfTNtxnNS9zgKbm81g4nlTxHx15S+PdA7WbEZ711FxNE07mRYqV2BJ3L3f5XSyuZJL+j1wE5V9q2lzCIUEQtk/7iL0ini6QW0u9Rg4qhW5q/6KLMIRcTS64p0KvjzDRhmqFcc/wjcFxHP0ZznpIg4ejXy+disnUm/YpH0FkmrgQeAk9PjzXjvKiIOaPz7FuBVWM0G4oCI+HUaZ/6xpJ9HxG3NbtRW7rj0nAwHrgM+RHYmoGVJeiNwPjC52W0ZjApxtN3zERF3A2+UNAaYK+lHzW7TQJSLIyL+TJPet3wmw0r9Gtgtd/9vU1rZPJJeAowAnqqyby1lFqGIWIiI3r+/A66n+GGUwcRRrcy/7aPMIhQRS/45eQb4Li3+nEj6W7LXzgkR8d+5/I1+ToqIoxnPx2btTAb02oqIh4Ae0jyTGsqstyLiaMb7FuBOhm1pKbCXpN0lvZRsUtGNJXluBE5M20cCi9MY8o3AMWl8eXdgL7KJbLWU2RaxSBqavp0haSjZt7dVLRxHWRHxW+CPkt6aTmWfAPyg/k3fQt1jkfQSSTun7SHAFFr4OZH0CuCHwKyI+Flv5iY9J3WPo0nPBwwult3ThzWSRgF7k02UbMZ7V93jaNL7VqbRM019a/0b8D7gEbIZzqentC8Ch6Xt7YFrySZD3gPskdv39LTfw+Rmxpcrsx1jIZvxvSLdVjcqlkHGsQb4Pdm3ml+RZqoDE8neaP4buIR0cb52i4XsVyf3AivTczKb9EugVowDOAN4Flieu+3SrOek3nE06/kYZCwfSm1dDtwHHFGtzHaLgya9b0WEr/hpZmZmxfBwiZmZmRXCnQwzMzMrhDsZZmZmVgh3MszMzKwQ7mSYmZlZIdzJMDMzs0K4k2FmZmaF8NolZlY4ScPILiBUzVPhC/eYdRRfjMvMCifpW8BH+sj2qohY14j2mFljuJNhZoWT1AUsAdaSXTL5T2Wy+UyGWYfxnAyzFiLpjgrpZ0k6dRDl9gy8VYMXEfcC/wC8HrgU+ENErCu5bdHBkLRE0uhyZUr6mKT/kbRc0i8lTU3pO0i6NT2+PN2ez21/ZSAxSPpbSUeXpL1U0m29i1I1Wg3HYFtJH6nXcWjFY9Bf1V5TVn/uZJgVRJl+/Y9FxNuLas9AlMZQa0zl8kXEj4FpwKHA1+vQvHHAWRExgWwlyn9L6dOA70fEN9JjhwJPRMSEdPv0AOv7e+DN+YSI+AtwC3B02T2K19cx2BQRl9fxOLTiMbAW5k6GWR1JGi3pYUlXkq2muZuk4yXdk749fiN9uxwq6YeSVkha1fvtMH/GQdLpkh6R9FPg73Llr8rlOVXSWWn7Bkn3SlotaXqN7S3XttIY3lEmps+kdq+S9KlKsZfWFxHfAWYCH5X0hf4f4c3sA/w8bf8K2DZtH8fmS6SPBR6opcB0BuB9ki6W9L5c+gHAvwNHpmO1R263G1KdzVDrMYAaj0MbHgNrYW1xesuszewFnBgRd0kaQ/YNb/+I2Cjpa2Rvxs8Cv4mIQwEkjcgXoGwOwzHABLL/0/vIls+uZlpE/F7SDsBSSddFxFOVMldp220lMYwuud8FnAS8BRBwt6RbgT/k81WqNyIukPQa4GxJv4qIb/cRVyXjgIckCfgksEDSS8mWvV5Tkm9Vmf0BkLQn8F6yuSK7Ad3Af6W/vW3+qaSlwKkRUVrWKmDfXHmvjIg/DDCm/qr1GPTmLXsc6n0M+jLYY9TgY2yD4DMZZvW3Nvch+/dAF9mH/vJ0fw+yb5TvlnS+pHdExPqSMt4BXB8Rf4qIPwI31lDvJyWtAO4i+6DYq4/8ldpWGkPp/QNS256NiB7g+6m95far5HyyyZ+f6e+QEoCk3YBhwELgHuCVwCnAzsDTJdn7+gZ/HXABsAJ4V0TMiIgfRcSfS/L9HS+eNUDSvwBExCbgL5KGp4e2mOuQOgE1k/ST3Jmi/O3wXJ7+HAOofhzqfQz68sIx6i2jn/qcTyLpPZI+NICyrY58JsOs/p7NbQuYGxH/rzSTpDeTfXP8kqRbIuKLNZT9Vzb/crB9Kusg4F3A2yLiT5KW0Pd1Kcq2LZ25eLYkb+n9SvrMJ+llZJ2mZ4HDIuL5GsvOGwfcEhGHlJS9DVvGPY7qH0oTgDeRPRfzJQ0BFgHfiojHU7k7A+sj4q/p/t8AQ3JlbAf8WdIhwN6STgOuAq4nG0q4Q9JbIuJCSZcCnyMbNnoZsE1EfDLfoIh4V52PQV/HoZ7HYAhwdj621CH6ArCe7FdG+WM0JJV1DfBD4I3AHcC7yeabrEpDgq8EniLrUPXu/9V8XcCpZEM6fyQ703ZKtQNoxfOZDLNi3UI2hr0LgKQdJY1KwwV/SnMULqBkMh3ZkMURaXx8OPD+lP6/wC6SdpK0HTAlpY8g+8XGnyTtDbx1oG2rYb/bU9teJmko8IGU1qf0Afhdsg+8wyLil7XsV8Y+ZN+6N5NOoW8rqbfztQ3ZGZ2HKhUUmfsi4ksRcQDZsMFq4FW5bKOB3+TuTwCWpzp2AtZFxEZgHfCdiLgg5fleRJwP7J5r7zDgeGAHsjMOmw2V9UNNxyC1sepxqPMxmF4mtr2BvwAXAf/D5sdoefp7dUR8Oe3zLeBaYJSkXcm+ED8N7M/mx7i0ro+TdZw/R9bpebjK8bMG8JkMswJFxIOSzgAWpTf6jWTfrkYAF0h6PqV9vGS/+yRdQ/Yh8jtgaUrfKOmLZN/mfs2Lp65vBk6W9BDZG2ufQxZV2vY/fex3n6Q5qQ0Al0fE/artZ4FfAQ4D/qnGYZVKxpHNGShnEdmQzk+APYFfpV9AlCXpJ8DflHnodF6cB/NzYGdlk26nk30o3pAem0T2DRw2/+CfwIuTL8cBCyW9HAiyswanRMRzVWLsS63HAPo4DnU+BlvEFhE/lvQEcAlwN1knmlwZRwA3pLMgT0XE85LGknU2/gX4Z7IOz25sfow3qyu9Lr+ROuZlfxZtjeWLcZlZQ0j6Z7LT259N30Jr2WcJMLXMJMZq+7wZ+HREFDYeL+k/gI+mD8PvA7Mi4hFJhwH/CJxHduq+N8+Hyb7N/5VsyGA1cCzwBLA4Im6uc/uaeQzeT0lsks4n++XL9sCPya6Z8sIxIutMfJSsA/GuNKx0VUR8SNL/BV4O7ETWgV7Li8d4z3xdZEMq7yGb79MTEZ8t0+4l9PM1ZQPnToaZFS5NWPw+2STFz1TI9nhEbHYl0IF+IEiaRnbafFP/W9uvel4KHBMRVxZZz0D4GJTnTkZjuZNhZoWTtJLs9H41kyJiScl+U4EbIqLcryXM+s2vqcZyJ8PMzMwK4V+XmJmZWSHcyTAzM7NCuJNhZmZmhXAnw8zMzArhToaZmZkVwp0MMzMzK4Q7GWZmZlYIdzLMzMysEP8fSx/XXrXI93AAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Assemble residuals to make datasframe \n", "all_residuals = pd.DataFrame([res_recov_stop_time_type, res_recov_time_type, \\\n", " res_recov_type, res_recov_trip_id, \\\n", " res_recov_stop_type, res_recov_stop_time, res_recov_stop, \\\n", " res_recov_time, res_fail], \\\n", " index = ['recov_stop_time_type', 'recov_time_type', \\\n", " 'recov_type', \\\n", " 'recov_trip_id', 'recov_stop_type',\\\n", " 'recov_stop_time', \\\n", " 'recov_stop', 'recov_time', \\\n", " 'baseline average error']).transpose()\n", "# Sort before plotting\n", "meds = all_residuals.median()\n", "meds.sort_values(ascending=False, inplace=True)\n", "all_residuals = all_residuals[meds.index]\n", "\n", "# Boxplot or residuals per recovery tables\n", "plt.figure(figsize = (8,6))\n", "plt.gcf().subplots_adjust(bottom=0.5, left = 0.2)\n", "all_residuals.boxplot(vert=False, showfliers=False)\n", "plt.xlabel('residual error $\\sum $ | $P(T \\leq t_i)_{true} - P(T \\leq t_i)_{estimated} $ |')\n", "plt.savefig(\"../data/figs/recov_tab_residuals_final.pdf\", dpi=300)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The order of importance of recovery tables is therefore :\n", "1. recov_trip_id \n", "2. recov_stop_time_type\n", "3. recov_stop_time\n", "4. recov_stop\n", "5. recov_time_type\n", "6. recov_type\n", "7. recov_time" ] }, { "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 multiple 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 tables` with different combinations of parameters, are used in the order defined in previous section\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": 211, "metadata": {}, "outputs": [], "source": [ "def try_recovery(recov_df, tup, keep_trying, n_recov, all_distrib):\n", " if keep_trying and tup in recov_df.index:\n", " distrib = np.array(recov_df.loc[tup])\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_recov += 1\n", " return n_recov, keep_trying, all_distrib" ] }, { "cell_type": "code", "execution_count": 221, "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 FINAL DISTRIBUTION TABLE USING RECOVERY TABLES ################\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", "\n", "# declare empty variables \n", "size_stop_times = stoptimes.shape[0]\n", "\n", "n_fail = n_real = n_all = n_recov_stop_time_type = n_recov_stop_time =\\\n", "n_recov_stop_type = n_recov_time_type = n_recov_stop =\\\n", "n_recov_time = n_recov_type = n_recov_trip_id = i = 0\n", "\n", "all_distrib, all_transport_type, all_hours, all_keys = \\\n", " ([] for i in range(4))\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", " \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", " \n", " # 3) try recovery tables in order previously defined\n", " n_recov_trip_id, keep_trying, all_distrib =\\\n", " try_recovery(recov_trip_id, (trip_id), keep_trying, n_recov_trip_id, all_distrib)\n", " n_recov_stop_time_type, keep_trying, all_distrib =\\\n", " try_recovery(recov_stop_time_type, (stop_id, hour, transport_type), \\\n", " keep_trying, n_recov_stop_time_type, all_distrib)\n", " n_recov_stop_time, keep_trying, all_distrib =\\\n", " try_recovery(recov_stop_time, (stop_id, hour), keep_trying, n_recov_stop_time, all_distrib)\n", " n_recov_stop, keep_trying, all_distrib =\\\n", " try_recovery(recov_stop, (stop_id), keep_trying, n_recov_stop, all_distrib)\n", " n_recov_time_type, keep_trying, all_distrib =\\\n", " try_recovery(recov_time_type, (hour, transport_type), keep_trying, n_recov_time_type, all_distrib)\n", " n_recov_type, keep_trying, all_distrib =\\\n", " try_recovery(recov_type, (transport_type), keep_trying, n_recov_type, all_distrib)\n", " n_recov_time, keep_trying, all_distrib =\\\n", " try_recovery(recov_time, (hour), keep_trying, n_recov_time, all_distrib)\n", " \n", " # save parameters for summary\n", " all_keys.append(key)\n", " all_hours.append(hour)\n", " all_transport_type.append(transport_type)\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", " n_fail += 1 \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": 224, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAADTCAYAAABJG/MPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de5zVVb3/8ddbVFQgFC8dJHW8UCpeSEbLS0pmZGqpj5+khSlSeix/Zp6y6JeVWmadTpnXlCzwgolYdrycBBRFVFSGO6ioqRy8lHcENS/4+f2x1tav272HmWH2DHvm/Xw89mO+3/Vd37XWd82W+bjW2nspIjAzMzOrF2t1dgPMzMzMWsPBi5mZmdUVBy9mZmZWVxy8mJmZWV1x8GJmZmZ1xcGLmZmZ1RUHL2ZdkKSQtF1ntwNA0lBJT9ao7Ib8rGvn879JOradyv6UpMWF8yckHdAeZefyFkka2l7lNVPPGZKuqnU9aypJH5Z0p6Tlkn7d2e2x9rF2ZzfAzLoWSQEMjIhHO7ruiPh8S/K1pI0RMR34WHu0S9I44MmIOL1Q/qD2KLs9VWpnF3AC8DzwofAXm3UZHnkxs6pKIxrdTXd97i5qK+CBaoGLf9f1ycGLWZ2QdJykGwvnj0iaWDhfKmlw4ZYDcp6XJV0kSYW8oyQ9KOklSZMkbVW4FpJOkvQI8EhOO0TS3FzWPZJ2qdLGO/PhPEkrJB1ZuPYdSc9KekbScYX0npL+S9L/SvqnpEskrV+l/B457/OSHgMOLrt+h6Sv5+PtJE2TtCznn1CtjaWpLUnfl/QPYGyV6a7dJT2Q+22spPVymSMl3VXWlshtOAEYAXwv13djvv7uNFTug99Kejq/fiupZ75WalvF/qvQR1vn514uaQqwSdn1iZL+kfvlTkmDcnq1do6W9Pdc3gOSDm+m7h6S/l8h/yxJW+Rre0mameudKWmvst/bz/J7a4WkGyVtLGm8pFdy/oZC/u0lTZH0oqTFkr5UpT3jgGMLz3SA0jTadZKukvQKMFLS5pJuyOU9Kun4Qhln5D67Kj/TAkkflfSD/PtYKmlYtT6xGokIv/zyqw5ewDbAy6T/6dgcWEIa4i9dewlYK58HcBOwIbAl8BxwYL52KPAosANp6vh04J5CPQFMAfoB6wMfB54FPgH0IP0xeALoWaWdAWxXOB8KvA2cBawDHAS8BmyUr58L3JDr6wPcCJxTpewTgYeALXL+23N9a+frdwBfz8d/An6Y+2s9YJ8WtPGXQM/83ENL/ZvzPAEsLNR9N/CzfG0kcFe1fgDGlfKWlXdAPj4LuBfYDNgUuAf4aUv6r0IfzQB+k59jX2A5cFXh+qjczz2B3wJzC9cqtXM46f22FnAk8CrQv0rdpwELSNNtAnYFNs799RLwVdJ77sv5fOPC7+1RYFugL/AA8DBwQM5/BTA25+0FLAWOy9c+TpoW2rFKm973TMAZwFvAYfmZ1gfuBC7O75PBpP9e9i/k/xfwuUJbHie9t9YBjgce7+x/H7rbq9Mb4JdffrX8lf/R3g04ChgD3A9sn/8hv6GQL3j/H+trgdH5+G/A1wrX1sp/DLcq3Lt/4frvSn9IC2mLgf2qtLFSYPA6OcDIac8Cn8x/4F4Fti1c27PaHwNgKnBi4XwY1YOXK3IffaSFbXwTWK8srTx4KdZ9EPD3fDyS1Qte/g4cVLj2OeCJVfVfhefakhTo9CqkXU0heCnLv2FuZ99q7axwz1zg0CrXFle6Rgpa7i9LmwGMLPzefli49mvgb4XzL5CDLFIANb2srEuBn1Rp0/ueiRSM3Fk43wJYCfQppJ0DjCvkn1LWlhVAj3zeJ/fhhi3979iv1X952sisvkwj/THbNx/fAeyXX9PK8v6jcPwa0DsfbwWcpzQF9DLwIimIGFDIv7RwvBXwnVL+fM8WpP8bb6kXIuLtCu3ZFNgAmFUo+5acXsnmZW1b0kyd3yM91/1Kn+wZtYo2PhcR/1pFnvK6W9MHzSmNpFUru1r/VSrnpYh4taws4N1pnV/kaZ1XSAEUlE0tFUk6Ru9NGb4M7NRM/i1IgVildpX/rpbw/vfcPwvHr1c4L75/P1H2fhwB/Fu1Z6ig+HvcHHgxIpa3om3PR8TKwjlU/n1YjXihkll9mUb6P7+tgZ+TppFGkEYrLmxhGUuBsyNifDN5iosbS/nPbn1zV+l50j/+gyLiqRbkf4b0B7Jky2oZI+IfpCF9JO0D3Crpzqj+CaOWfBKlvO6n8/GrpCCMXF/5H9JVlf006Y/yogplt8YzwEaSehUCmC0L9X+FNG14AClw6Uuavimth3pfO5XWQv0e+AwwIyJWSppbyF9uKWnqZ2FZeun5irYkBaqttRSYFhGfbcO9JcXnfBroJ6lPIYDZEmjJ+9E6iUdezOrLNODTwPoR8SQwHTiQtK5gTgvLuAT4QWGhZl9Jw5vJ/3vgREmfUNJL0sGS+lTJ/0/SGpxVioh3cvnnStost2eApM9VueVa4FuSPiJpI2B0tbIlDZf0kXz6EukP1jutbWOZk3Ld/UhrHibk9HnAIEmDlRbxnlF236rq+xNwuqRNJW0C/Bho9XezRMQSoAk4U9K6OWj7QiFLH+AN4AVSsPXzVbSzF6nfnoO0aJw08lLNZcBPJQ3M75VdJG0M/A/wUUlfkbS20kLuHUnrslrrplzWVyWtk1+7S9qhDWUREUtJa4zOkbSe0mL0r9GG/reO4+DFrI5ExMOk+fbp+fwV4DHg7sIw9qrKuJ60MPWaPHWwEKj6/SgR0UQawbiQFAQ8SlrjUc0ZwOV5SL/ip0DKfD+XeW9uz61U/36V3wOTSMHCbOAvzZS7O3CfpBWkBcGnRMRjbWxjydXAZFKf/x34Gbz7ezkrt/0R4K6y+/4A7Jjr+2uFcn9GCjrmkxa8zi6V3QZfIS2ufhH4CWntT8kVpCmRp0iLYu9trp0R8QBp/ckMUmCzM2mhcjW/IQWYk4FXcnnrR8QLwCHAd0iB0/eAQyLi+dY+XB4dGUZa9/U0aXq0tNC6rb4MNOTyrietn7l1NcqzGlNecGRmZmZWFzzyYmZmZnXFwYuZmZnVFQcvZmZmVlccvJiZmVldcfBiZmZmdcVfUmc1sckmm0RDQ0NnN8PMzOrYrFmzno+ID3zjtoMXq4mGhgaampo6uxlmZlbHJFXcAsTTRmZmZlZXHLyYmZlZXXHwYmZmZnXFwYuZmZnVFQcvZmZmVlf8aSOriQVPLaNh9M0dUtcTvzi4Q+oxM7M1g0dezMzMrK44eDEzM7O64uBlNUl6QtImVa6NljSiA9syVNJN7VTWYZJ2bI+yzMzM2pODl0xJe/fH54DJ7VxmRzkMcPBiZmZrnG4dvEhqkLRY0hXAQmALSadJmilpvqQzC3n/KmmWpEWSTmhB2R8C1o2I58rSN5U0JZdzmaQlpZEbSUdLul/SXEmXSuqRX+MkLZS0QNKpOe92km6VNE/SbEnb5ip6S7pO0kOSxktSzv/j/FwLJY0ppG8r6Zb8bNMlbS9pL+CLwK9yW/bOP0uvlZK2Wv3fgJmZWet16+AlGwhcHBGDgI/l8z2AwcAQSfvmfKMiYgjQCHxL0sarKPcA4LYK6T8Bpub6rgO2BJC0A3AksHdEDAZWAiNyOwZExE4RsTMwNpczHrgoInYF9gKeyekfB75NGjXZBtg7p18YEbtHxE7A+sAhOX0McHJ+tu/mvrgHuAE4LSIGR8Td+edg4PfAnyOi4n4TZmZmteaPSsOSiLg3Hw/Lrzn5vDcpmLmTFLAcntO3yOkvNFPugbwXaBTtAxwOEBG3SHopp38GGALMzIMi6wPPAjcC20i6ALgZmCypDymguT6X8y+AfN/9EfFkPp8LNAB3AZ+W9D1gA6AfsEjS7aTAZ2K+F6BntQeStDdwfH6GStdPAE4A6PGhD2wCamZm1i4cvMCrhWMB50TEpcUMkoaSRlL2jIjXJN0BrLeKcvcAvtGKdgi4PCJ+8IEL0q6k9TMnAl8CTmmmnDcKxyuBtSWtB1wMNEbEUklnkNq/FvByHlFpvnFSf+APwBcjYkWlPBExhjSSQ8/+A2NVZZqZmbWFp43ebxIwSlJvAEkDJG0G9AVeyoHL9sAnmytE0iDgoYhYWeHy3aQABEnDgI1y+m3AEbk+JPWTtFVeD7NWRPwZOB3YLSKWA09KOizn7Slpg2aaVAq0ns/PdgRARLwCPC5peC5HOVACWA70yenrABOB70fEw809u5mZWa05eCmIiMnA1cAMSQtIa1L6ALeQRjAeBH4B3Fu9FAA+n++p5ExgmKSFwHDgH8DyiHiAFJxMljQfmAL0BwYAd+QpoKuA0sjMV0lTWfOBe4B/a+a5XiatVVlICtBmFi6PAL4maR6wCDg0p18DnCZpDmlqqRE4s7Bod/NV9IGZmVlNKMKj++1N0hTgmIh4psK1nsDKiHhb0p7A71oybVNvevYfGP2P/W2H1OXtAczMuiZJsyKisTzda15qICI+28zlLYFr83fKvElaAGtmZmYt5OClg0XEI6SPM5uZmVkbOHixmth5QF+aPJ1jZmY14AW7ZmZmVlccvJiZmVld8bSR1cSCp5bRMPrmzm6GmZl1sI74BKhHXszMzKyuOHgxMzOzulLz4EXSE/kr7itdGy1pRK3bUKhvqKSb2qmswyTt2Mp71pE0uz3qNzMz665aHLzkfW/aO9j5HDC5ncvsKIcBrQpeSLsx392WyiR5fZKZmRmrCF4kNUhaLOkK0r44W0g6TdJMSfMlnVnI+1dJsyQtknTCqiqW9CFg3Yh4rix9U0lTcjmXSVpSGrmRdLSk+/PeOpdK6pFf4yQtlLRA0qk573aSbpU0T9JsSdvmKnpLuk7SQ5LGS1LO/+P8XAsljSmkbyvplvxs0yVtL2kv4IvAr3Jb9i7s+TNX0kpJW1V47AOBv1XoixWSzs3PfJukTXP6HZJ+K6kJOEXSZyTNyc/5x7zVAJIOys8zS9L5pdElSWfkfHdIekzStwp1/kd+1oWSvp3Tekm6OffZQklH5vQhkqbl8icp7TBtZmbWKVoykjIQuDgiBgEfy+d7AIOBIZL2zflGRcQQ0gZ+35K08SrKPYC0k3K5nwBTc33Xkb5OH0k7AEcCe+e9gFaSNhUcDAyIiJ0iYmdgbC5nPHBRROxK2liwtM/Qx4Fvk0ZNtgH2zukXRsTuEbETsD5wSE4fA5ycn+27uS/uAW4ATouIwRFxd/45mLQB4p8jYkmFZ/s0cEeF9F5AU37mabkPStbN+zpcBIwDjszPuTbwDUnrAZcCn89t3LSs7O1JI1x7AD/JU1dDgOOAT5B2yD5e0sdJwdXTEbFr7odblHaUvgA4Ipf/R+DsCs9gZmbWIVoyFbEkIkq7KA/Lrzn5vDcpmLmTFLAcntO3yOkvNFPugbwXaBTtAxwOEBG3SHopp38GGALMzIMi6wPPAjcC20i6ALiZtCtzH1JAc30u518A+b77I+LJfD4XaADuAj4t6XvABkA/YJGk20mBz8R8L0DPag8kaW/SXkX7VLg2AHgxIl6rcOs7wIR8fBXwl8K1UvrHgMcj4uF8fjlwEikYeiwiHs/pfwKKI183R8QbwBuSngU+nNt3fUS8mtv2F+BTpJ2wfy3pl8BNETFd0k7ATsCU3Ac9eC8QLH/GE0p19/hQeQxlZmbWPloSvLxaOBZwTkRcWswgaShpJGXPiHhN0h3Aeqsodw/gGy1vKgIuj4gffOCCtCtpdOFE4EvAKc2U80bheCWwdh69uBhojIilks4gtX8t4OWW7Pqcp1L+AHwxIlZUyHIgMGlV5WTFrb5frZqrZT7wvFUrjXhY0m7AQcDPJN0GXA8siog9V1VRRIwhjVTRs/9Ab1duZmY10doFuJOAUZJ6QxpNkLQZ0Bd4KQcu25OmIqqSNAh4KCJWVrh8NykAQdIwYKOcfhtwRK4PSf0kbZXXw6wVEX8GTgd2i4jlwJOSDst5e0raoJkmlQKt5/OzHQEQEa8Aj0sanstRDpQAlgN9cvo6wETg+4WRkXIV17tka5XqBL5CGgkqtxhokLRdPv8qaYppMWnkqSGnH1n9Md81HThM0gaSepFGuqZL2hx4LSKuAn4F7JbL31TSnqVnzb8/MzOzTtGq4CUiJgNXAzMkLSCtSelDmm5YW9KDwC+Ae6uXAsDn8z2VnAkMk7QQGA78A1geEQ+QgpPJkuYDU4D+wADgjjwFdBVQGpn5Kmkqaz5wD/BvzTzXy6S1KgtJAdrMwuURwNckzQMWAYfm9GuA0yTNIU0tNQJnFhbtbl4qQFIPYLuIeKhKE14F9sjPvD9wVoU2/ou0TmVi7vt3gEsi4nXgm6T1KbNIQdWyas+ay5pNWj9zP3AfcFlEzAF2Bu7PffkT4GcR8SYpsPpl7oO5+XnNzMw6hSI6fnRf0hTgmIj4wNqJ/AmalRHxdv6//d+1ZNpmTSZpH+DoiDixyvUVEdF7NcrvHRErlBalXAQ8EhHntrW89tCz/8Dof+xvO7MJZmbWCdpzewBJs/KHVt6nU747JCI+28zlLYFrlb5T5k3SAti6FhF3UXkqqL0cL+lYYF3SYupLV5HfzMysbq1xX3wWEY+QPs7cbazOqEu+/1ygU0dazMzMOsoaF7xY17DzgL40dcDOomZm1v14Y0YzMzOrKw5ezMzMrK542shqYsFTy2gYfXPN62nPVe1mZlYfPPJiZmZmdcXBi5mZmdUVBy9mZmZWVxy81DlJI4tbEVS4fpakA9pQ7omSjqmQ3pC3MTAzM+sUXrDbjvLX8ysi3unAakeS9mR6ukJ7ekTEj9tSaERcsprtMjMzqwmPvKymPBKxWNIVpCDiR5JmSpov6cxCvmNy2jxJVxbunZrTb5O0paS+kpbk7RGQ1EvS0rxzdXndR5A2hByfN4NcX9ITkn4paTYwXNK4nI987T8lLZB0f2GH6krPdYak7+bjIbnd84CT2q/3zMzMWs/BS/sYCFwMnEra5XoPYDAwRNK+kgaRdsTePyJ2BU7J910AXB4RuwDjgfMjYhlp5+b9cp5DgEkR8VZ5pRFxHdAEjIiIwXmHaYAXImK3iLimQluXRcTOwIVAS3dOHAucnNtelaQTJDVJalr5WrMbW5uZmbWZg5f2sSQi7gWG5dccYDawPSmw2R+YGBHPA0TEi/m+PYGr8/GVwD75eAJwZD4+Kp+3RnP5/1T4ueeqCpK0IbBhRNxZaGdFETEmIhojorHHBn1b3FgzM7PWcPDSPl7NPwWck0dBBkfEdhHxhzaUdwNwoKR+wBBgahvbU0lUOTYzM6sLDl7a1yRglKTeAJIGSNqMFHwMl7RxTu+X899DGlkBGAFMB4iIFcBM4DzgpohY2Uydy4E+rWjjkYWfM1aVOSJeBl6WVBoVGtGKuszMzNqdP23UjiJisqQdgBnpg0esAI6OiEWSzgamSVpJmlYaCZwMjJV0GvAccFyhuAnARGDoKqodB1wi6XVaMA0EbCRpPvAG8OUWPtpxwB8lBTC5hfeYmZnVhCI8c9BdSHoCaCytvamlnv0HRv9jW7oeuO28t5GZWdclaVZENJane9rIzMzM6oqnjeqEpIuAvcuSz4uIsS0tIyIaKpT7Q2B4WfLEiDi71Y0s2HlAX5o8KmJmZjXg4KVORERNvhwuBymrFaiYmZl1JE8bmZmZWV3xyIvVxIKnltEw+uZ3z72w1szM2otHXszMzKyuOHgxMzOzuuLgxczMzOqKg5dVkDRU0l7tXOZgSQcVzr8oaXR71lFW37clbVCr8s3MzDrSGhu8KFkT2jcUaNfgBRgMvBu8RMQNEfGLdq6j6NuAgxczM+sS1oTg4F2SGiQtlnQFsBD4kaSZkuZLOrOQ75icNk/SlYV7p+b02yRtKamvpCWlIEhSL0lLJa1Tpf5vSXogl3GNpAbgROBUSXMlfapSPfnecZIukdQk6WFJh1SpY13gLODIXOaRkkZKurBQzu8k3SvpsTzy80dJD0oaVyhnmKQZkmZLmljaDLLSMwGbA7dLul3SKEm/LVw/XtK5+bkekjQ+13VdabRG0hBJ0yTNkjRJUv8W/ULNzMxqYI0KXrKBwMXAqcAAYA/SSMUQSftKGgScDuwfEbsCp+T7LgAuj4hdgPHA+RGxDJgL7JfzHAJMioi3qtQ9Gvh4LuPEiHgCuAQ4NyIGR8T0SvUU7m/I7T2YtFnieuUVRMSbwI+BCbnMCRXasRFpk8VTgRuAc4FBwM55ymmT3AcHRMRuQBPwH5UeKCLOB54GPh0RnwauBb5QCOCOA/6Yjz8GXBwROwCvAN/M+S4AjoiIITlvxS+1k3RCDt6aVr62rFIWMzOz1bYmBi9LIuJeYFh+zQFmA9uTApv9SV9f/zxARLyY79sTuDofXwnsk48nAEfm46PyeTXzgfGSjgberpKnWj0A10bEOxHxCPBYbnNb3Bhpx8wFwD8jYkFEvAMsIgVInwR2BO6WNBc4FtiqJQVHxApgKnCIpO2BdSJiQb68NCLuzsdX5Wf7GLATMCXXdTrwkSplj4mIxoho7LFB31Y/tJmZWUusiV9S92r+KeCciLi0eFHSya0s7wbg55L6AUNIf7irORjYF/gC8ENJO7eyrvItutu6Zfcb+ec7hePS+drASmBKRHy5jeVfBvw/4CGguDdSpfYLWBQRe7axLjMzs3a1Jo68lEwCRpXWckgaIGkzUvAxXNLGOb1fzn8PaWQFYAQwHd4daZgJnAfcFBErK1WW18VsERG3A98H+gK9geVAn0LWivVkwyWtJWlbYBtgcZVnKy+zte4F9pa0XW57L0kfbSb/++qLiPuALYCvAH8q5NtSUilI+QpwF+kZNi2lS1onT92ZmZl1ijU2eImIyaTpmRmSFgDXAX0iYhFpzcU0SfOA3+RbTgaOkzQf+CrvrYWBNFV0NM1PGfUArsp1zSGtmXkZuBE4vLRgdxX1/C9wP/A30pqZf1Wp63Zgx9KC3Zb0R1FEPAeMBP6U2zGD5qeoxgC3SLq9kHYtcHdEvFRIWwycJOlB0rqb3+U1OkcAv8z9PZf2//SVmZlZiyktrbDVlT8JdFNEXNfZbWkJSTeRFiLfls8bSO3fqT3K79l/YPQ/9t0PNXlvIzMzazVJsyKisTx9jR15sdqQtKGkh4HXS4GLmZlZPemWIy+SLgL2Lks+LyLGVsq/GvV8DvhlWfLjEXF4e9ZTqO96YOuy5O9HxKRa1NecxsbGaGpq6uhqzcysC6k28rImftqo5iLipA6qZxJp4XGHqFVQZGZmtibxtJGZmZnVFQcvVhMLnvI37JqZWW04eDEzM7O64uDFzMzM6oqDl24s71jtL5wzM7O64uClEyhZE/p+KP62XDMzqzNrwh/QbkFSg6TFkq4AFgI/kjRT0nxJZxbyHZPT5km6snDv1Jx+m6QtJfWVtKQUBOX9jZZKWqdK/d+S9EAu45r8jbonAqeWtj6oVE++d5ykSyQ1SXpY0iG17S0zM7PquuX3vHSigcCxwIdI+wXtQdq1+QZJ+wIvAKcDe0XE84VNJy8ALo+IyyWNIu27dJikucB+pL2SDgEmRcRbVeoeDWwdEW9I2jAiXpZ0CbAiIv4LQNKN5fUAh+X7G3J7twVul7RdM3s3mZmZ1YxHXjrWkoi4FxiWX3OA2aRNFQcC+wMTI+J5gIh4Md+3J2mTSoArgX3y8QSgtLHjUTS/8eR8YLyko4G3q+SpVg/AtRHxTkQ8AjxGhY0gJZ2QR2eaVr7mj0qbmVltOHjpWK/mnwLOiYjB+bVdRPyhDeXdAByYR2iGAFObyXswcBGwGzBTUmtH3cr3kfjAvhIRMSYiGiOisccGfVtZvJmZWcs4eOkck4BRknoDSBogaTNS8DFc0sY5vTRtdA9pZAVgBDAdICJWADOB80g7Qq+sVFleF7NFRNwOfB/oC/QGlgN9Clkr1pMNl7SWpG2BbYDFbXx2MzOz1eI1L50gIiZL2gGYIQlgBXB0RCySdDYwTdJK0rTSSOBkYKyk04DngOMKxU0AJpI+OVRND+AqSX1Joz7n5zUvNwLXSTo019FcPf8L3E9ar3Oi17uYmVln6Za7SlvrSBpHGtm5rqX39Ow/MN545pHaNcrMzLq8artKe9rIzMzM6oqnjboYSRcBe5clnxcRY9taZkSMXK1GmZmZtSMHL11MRJzU2W0A2HmAP21kZma14WkjMzMzqysOXszMzKyuOHixmljwlL9h18zMasPBi5mZmdUVBy9mZmZWVxy8mJmZWV1x8NIJJI2UtHnh/DJJO9aorgZJX6lF2WZmZp2h2wcvSjq6H0YC7wYvEfH1iHigRnU1AA5ezMysy+iWwUsejVgs6QpgIfAjSTMlzZd0ZiHfMTltnqQrC/dOzem3SdpSUl9JS0pBkKRekpZKWqdC3UcAjcB4SXMlrS/pDkmN+foKSb+StEjSrZL2yNcfk/TFnKdHzlNq878387i/AD6V6zpV0p2SBhfac5ekXSWdIelKSTMkPSLp+EKe0yr1j5mZWWfolsFLNhC4GDgVGADsAQwGhkjaV9Ig4HRg/4jYFTgl33cBcHlE7AKMJ+3QvAyYC+yX8xwCTIqIt8orzZsbNgEjImJwRLxelqUXMDUiBgHLgZ8BnwUOB87Keb4GLIuI3YHdgeMlbV3lOUcD03Nd5wJ/II38IOmjwHoRMS/n3QXYH9gT+LGkzSUNy331vv6pVJGkEyQ1SWpa+Zo/Km1mZrXRnYOXJRFxLzAsv+YAs4HtSX+s9wcmRsTzABHxYr5vT+DqfHwlsE8+ngAcmY+Pyudt8SZwSz5eAEzLQdAC0hQQub3HSJoL3AdsnNvcEhOBQ/Ko0ChgXOHaf0fE6/mZbycFLNX65wMiYkxENEZEY48NvD2AmZnVRnfe2+jV/FPAORFxafGipJNbWd4NwM8l9QOGAFPb2K63IiLy8TvAGwAR8Y6k0u9LwMkRMam1hUfEa5KmAIcCX8ptffdyeXaq9I+ZmVln6c4jLyWTgFGSegNIGiBpM1LwMVzSxjm9X85/D2lkBWAEMB0gIlYAM4HzgJsiYmUzdS4H+qxmm79RWlMj6aOSerWirsuA84GZEfFSIf1QSevlZx5Kep5q/WNmZtYpuvPICwARMVnSDgbU9w4AAAygSURBVMAMSQArgKMjYpGks4FpklaSpk1GAicDYyWdBjwHHFcobgJpWmboKqodB1wi6XXSNFRrXUaaQpqt1OjngMOq5J0PrJQ0DxgXEedGxCxJrwBjK+S9HdgE+GlEPA08Xal/gGfb0G4zM7PVpvdmKKy7yN8xcwewfUS8k9POAFZExH+1Rx09+w+MN555pD2KMjOzbkrSrIhoLE/3tFE3I+kY0iLfH5YCl1rYeYAX7JqZWW10+2mjWpJ0EbB3WfJ5EVE+XdMede1M+vRT0RsR8YliQkRcAVxRfn9EnNHebTIzM6sFBy81FBEndWBdC0jfw2JmZtaledrIzMzM6opHXqwmFjy1jIbRN797/sQvDu7E1piZWVfikRczMzOrKw5ezMzMrK44eOkGJI3M3+1iZmZW9xy8dDAlHd3vIwEHL2Zm1iU4eOkAkhokLZZ0BbAQ+JGkmZLmSzqzkO+YnDZP0pWFe6fm9NskbSmpr6QlpSBIUi9JS0t7HZXVfQTQCIyXNFfSwZL+Wrj+WUnX5+MVks6VtCjXtWlO31bSLZJmSZouafta9peZmVlzHLx0nIHAxcCpwABgD9L3sgyRtK+kQcDpwP4RsStwSr7vAuDyiNgFGA+cHxHLgLnAfjnPIcCkiHirvNKIuA5oAkZExGDgf4DtS4EJaW+mP+bjXkBTRAwCpgE/yeljSLtYDwG+m5/DzMysUzh46ThLIuJeYFh+zQFmA9uTApv9gYkR8TxARLyY79sTuDofXwnsk48nAEfm46Py+SpF2szqSuBoSRvm8v+WL79TKOcqYJ+8m/RewERJc4FLgf6VypZ0gqQmSU0rX1vWkuaYmZm1mr/npeO8mn8KOCciLi1elHRyK8u7Afi5pH7AEGBqK+4dC9wI/IsUML1dJV+QAtyX86hNsyJiDGmUhp79B3rHTzMzqwmPvHS8ScCoPKKBpAGSNiMFH8MlbZzT++X895BGVgBGANMBImIFMBM4D7gpIlY2U+dyoE/pJCKeBp4mTVMV91laCzgiH38FuCsiXgEelzQ8t0uSdm3Lg5uZmbUHj7x0sIiYLGkHYIYkgBXA0RGxSNLZwDRJK0nTSiOBk4Gxkk4DniOtUSmZAEwEhq6i2nHAJZJeB/aMiNdJ62c2jYgHC/leBfaQdDrwLO9NS40AfpfT1wGuAea14fHNzMxWm9ISCOtuJF0IzImIPxTSVkRE7/Yov2f/gdH/2N++e+7tAczMrLUkzYqIxvJ0j7x0Q5JmkUZZvtPZbTEzM2stBy9diKSLgL3Lks+LiOK6FvJHnj+gvUZdAHYe0Jcmj7aYmVkNOHjpQiLipM5ug5mZWa3500ZmZmZWVxy8mJmZWV1x8GJmZmZ1xcGLmZmZ1RUHL2ZmZlZXHLyYmZlZXXHw0g1IGilp88L5ZZJ27Mw2mZmZtZWDlw6WNzbs6H4fCbwbvETE1yPigQ5ug5mZWbtw8NIBJDVIWizpCmAh8CNJMyXNl3RmId8xOW2epCsL907N6bdJ2lJSX0lLSkGQpF6Slkpap0LdRwCNwHhJcyWtL+kOSY35+gpJv5K0SNKtkvbI1x+T9MWcp0fOU2rzv9e+18zMzCpz8NJxBgIXA6cCA4A9gMHAEEn7ShoEnA7sHxG7Aqfk+y4ALo+IXUg7QZ8fEcuAucB+Oc8hwKSIeKu80oi4DmgCRkTE4LyjdFEvYGpEDAKWAz8DPgscDpyV83wNWBYRuwO7A8dL2rq8LkknSGqS1PTcc8+1tn/MzMxaxMFLx1kSEfcCw/JrDjAb2J4U2OwPTIyI5wEi4sV8357A1fn4SmCffDwBODIfH5XP2+JN4JZ8vACYloOgBUBDTh8GHCNpLnAfsHFu8/tExJiIaIyIxk033bSNzTEzM2ue9zbqOK/mnwLOiYhLixclndzK8m4Afi6pHzAEmNrGdr0VEZGP3wHeAIiIdySV3h8CTo6ISW2sw8zMrN145KXjTQJGSeoNIGmApM1IwcdwSRvn9H45/z2kkRWAEcB0gIhYAcwEzgNuioiVzdS5HOizmm3+RmlNjaSPSuq1GuWZmZm1mUdeOlhETJa0AzBDEsAK4OiIWCTpbGCapJWkaaWRwMnAWEmnAc8BxxWKmwBMBIauotpxwCWSXidNQ7XWZaQppNlKjX4OOKwN5ZiZma02vTdjYNZ+Ghsbo6mpqbObYWZmdUzSrIhoLE/3tJGZmZnVFU8bdSGSLgL2Lks+LyLGdkZ7zMzMasHBSxcSESd1dhvMzMxqzdNGZmZmVlc88mI1seCpZTSMvvnd8yd+cXAntsbMzLoSj7yYmZlZXXHwYmZmZnXFwUs3J+lbkh6UNL7K9UZJ5+fjkZIu7NgWmpmZvZ/XvNg3gQMi4slKFyOiibQrtZmZ2RrBIy/dmKRLgG2Av0n6vqQZkuZIukfSx3KeoZJu6tyWmpmZvccjL91YRJwo6UDg08CbwK8j4m1JBwA/B/5PpzbQzMysAgcvVtIXuFzSQCCAdVpbgKQTgBMAenxo0/ZtnZmZWeZpIyv5KXB7ROwEfAFYr7UFRMSYiGiMiMYeG/Rt9waamZmBgxd7T1/gqXw8shPbYWZm1iwHL1byn8A5kubg6UQzM1uD+Y9UNxcRDfnweeCjhUun5+t3AHfk43HAuI5qm5mZWSUeeTEzM7O64pEXq4mdB/SlyZsxmplZDXjkxczMzOqKgxczMzOrKw5ezMzMrK44eDEzM7O64uDFzMzM6oqDFzMzM6srDl7MzMysrigiOrsN1gVJWg4s7ux2dFGbkL4R2dqf+7a23L+101X7dquI2LQ80V9SZ7WyOCIaO7sRXZGkJvdtbbhva8v9WzvdrW89bWRmZmZ1xcGLmZmZ1RUHL1YrYzq7AV2Y+7Z23Le15f6tnW7Vt16wa2ZmZnXFIy9mZmZWVxy8WLuSdKCkxZIelTS6s9uzJpP0hKQFkuZKaspp/SRNkfRI/rlRTpek83O/zpe0W6GcY3P+RyQdW0gfkst/NN+rjn/KjiPpj5KelbSwkFbz/qxWR1dSpW/PkPRUfv/OlXRQ4doPcj8tlvS5QnrFfx8kbS3pvpw+QdK6Ob1nPn80X2/omCfuOJK2kHS7pAckLZJ0Sk73e7c5EeGXX+3yAnoAfwe2AdYF5gE7dna71tQX8ASwSVnafwKj8/Fo4Jf5+CDgb4CATwL35fR+wGP550b5eKN87f6cV/nez3f2M9e4P/cFdgMWdmR/VqujK72q9O0ZwHcr5N0x/7ffE9g6/5vQo7l/H4BrgaPy8SXAN/LxN4FL8vFRwITO7osa9G1/YLd83Ad4OPeh37vNvDzyYu1pD+DRiHgsIt4ErgEO7eQ21ZtDgcvz8eXAYYX0KyK5F9hQUn/gc8CUiHgxIl4CpgAH5msfioh7I/3LdEWhrC4pIu4EXixL7oj+rFZHl1Glb6s5FLgmIt6IiMeBR0n/NlT89yGPAuwPXJfvL/89lfr2OuAzXW0EMSKeiYjZ+Xg58CAwAL93m+XgxdrTAGBp4fzJnGaVBTBZ0ixJJ+S0D0fEM/n4H8CH83G1vm0u/ckK6d1NR/RntTq6g/+bpy7+WJhyaG3fbgy8HBFvl6W/r6x8fVnO3yXlabGPA/fh926zHLyYdZ59ImI34PPASZL2LV7M/5fkjwO2k47oz272O/sdsC0wGHgG+HXnNqe+SeoN/Bn4dkS8Urzm9+4HOXix9vQUsEXh/CM5zSqIiKfyz2eB60nD6v/Mw7zkn8/m7NX6trn0j1RI7246oj+r1dGlRcQ/I2JlRLwD/J70/oXW9+0LpKmPtcvS31dWvt435+9SJK1DClzGR8RfcrLfu81w8GLtaSYwMH9yYF3SArsbOrlNayRJvST1KR0Dw4CFpP4qfUrgWOC/8/ENwDH5kwafBJbl4d5JwDBJG+Vh+2HApHztFUmfzGsEjimU1Z10RH9Wq6NLK/3Ryw4nvX8h9cdR+ZNCWwMDSQtGK/77kP+P/3bgiHx/+e+p1LdHAFNz/i4jv5/+ADwYEb8pXPJ7tzmdvWLYr671Iq2Ef5j0qYIfdnZ71tQX6RMX8/JrUamvSPP5twGPALcC/XK6gItyvy4AGgtljSItinwUOK6Q3kj6g/J34ELyl1J21RfwJ9L0xVukef2vdUR/VqujK72q9O2Vue/mk/4I9i/k/2Hup8UUPuVW7d+H/N/D/bnPJwI9c/p6+fzRfH2bzu6LGvTtPqTpmvnA3Pw6yO/d5l/+hl0zMzOrK542MjMzs7ri4MXMzMzqioMXMzMzqysOXszMzKyuOHgxMzOzuuLgxczMzOqKgxczMzOrKw5ezMzMrK78fwQz7/od3C45AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure(figsize = (8,6))\n", "plt.barh(np.array(['real / geschaetz','real / geschaetz / prognose','recov_trip_id',\\\n", " 'recov_stop_time_type','recov_stop', 'recov_time_type',\\\n", " 'recov_type', 'recov_time', 'fail'])[::-1],\\\n", " np.array([n_real, n_all, n_recov_trip_id, n_recov_stop_time_type, \\\n", " n_recov_stop_time, n_recov_stop, n_recov_time_type, \\\n", " n_recov_type, n_recov_time])[::-1])\n", "plt.title('where the distribution data come from')\n", "plt.gcf().subplots_adjust(bottom=0.5, left = 0.2)\n", "plt.savefig(\"../data/figs/data_final_origin.pdf\", dpi=300)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 225, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
keyhourtransport_typedistribution
02064.TA.26-13-j19-1.24.H__85762407Tram[16, 497, 62, 14, 2, 3, 4, 0, 0, 1, 0, 0, 0, 0...
12064.TA.26-13-j19-1.24.H__85913537Tram[2, 512, 58, 18, 3, 2, 2, 1, 0, 0, 1, 0, 0, 0,...
22064.TA.26-13-j19-1.24.H__85910397Tram[0, 417, 69, 11, 2, 0, 0, 1, 0, 0, 1, 0, 0, 0,...
32064.TA.26-13-j19-1.24.H__85911217Tram[0, 405, 80, 12, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0,...
42064.TA.26-13-j19-1.24.H__85914177Tram[0, 389, 81, 24, 4, 1, 0, 1, 0, 0, 0, 0, 1, 0,...
\n", "
" ], "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 \n", "0 [16, 497, 62, 14, 2, 3, 4, 0, 0, 1, 0, 0, 0, 0... \n", "1 [2, 512, 58, 18, 3, 2, 2, 1, 0, 0, 1, 0, 0, 0,... \n", "2 [0, 417, 69, 11, 2, 0, 0, 1, 0, 0, 1, 0, 0, 0,... \n", "3 [0, 405, 80, 12, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0,... \n", "4 [0, 389, 81, 24, 4, 1, 0, 1, 0, 0, 0, 0, 1, 0,... " ] }, "execution_count": 225, "metadata": {}, "output_type": "execute_result" } ], "source": [ "summary_df = pd.DataFrame([all_keys, all_hours, all_transport_type, all_distrib],\\\n", " index = ['key','hour','transport_type','distribution']).transpose()\n", "summary_df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Write down summary table used to make final distribution" ] }, { "cell_type": "code", "execution_count": 226, "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": 227, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.0266223 , 0.85357737, 0.95673877, 0.98003328, 0.98336106,\n", " 0.98835275, 0.99500832, 0.99500832, 0.99500832, 0.99667221,\n", " 0.99667221, 0.99667221, 0.99667221, 0.99667221, 0.99667221,\n", " 0.99667221, 0.99667221, 0.99667221, 0.99667221, 0.99667221,\n", " 0.99667221, 0.99667221, 0.99667221, 0.99667221, 1. ,\n", " 1. , 1. , 1. , 1. , 1. ,\n", " 1. , 1. ],\n", " [0.00332779, 0.85524126, 0.95174709, 0.98169717, 0.98668885,\n", " 0.99001664, 0.99334443, 0.99500832, 0.99500832, 0.99500832,\n", " 0.99667221, 0.99667221, 0.99667221, 0.99667221, 0.99667221,\n", " 0.99667221, 0.99667221, 0.99667221, 0.99667221, 0.99667221,\n", " 0.99667221, 0.99667221, 0.99667221, 0.99667221, 1. ,\n", " 1. , 1. , 1. , 1. , 1. ,\n", " 1. , 1. ],\n", " [0. , 0.82902584, 0.96620278, 0.98807157, 0.99204771,\n", " 0.99204771, 0.99204771, 0.99403579, 0.99403579, 0.99403579,\n", " 0.99602386, 0.99602386, 0.99602386, 0.99602386, 0.99602386,\n", " 0.99602386, 0.99602386, 0.99602386, 0.99602386, 0.99602386,\n", " 0.99602386, 0.99602386, 0.99602386, 1. , 1. ,\n", " 1. , 1. , 1. , 1. , 1. ,\n", " 1. , 1. ],\n", " [0. , 0.80516899, 0.96421471, 0.98807157, 0.99005964,\n", " 0.99204771, 0.99204771, 0.99403579, 0.99403579, 0.99403579,\n", " 0.99403579, 0.99602386, 0.99602386, 0.99602386, 0.99602386,\n", " 0.99602386, 0.99602386, 0.99602386, 0.99602386, 0.99602386,\n", " 0.99602386, 0.99602386, 1. , 1. , 1. ,\n", " 1. , 1. , 1. , 1. , 1. ,\n", " 1. , 1. ],\n", " [0. , 0.77335984, 0.93439364, 0.98210736, 0.99005964,\n", " 0.99204771, 0.99204771, 0.99403579, 0.99403579, 0.99403579,\n", " 0.99403579, 0.99403579, 0.99602386, 0.99602386, 0.99602386,\n", " 0.99602386, 0.99602386, 0.99602386, 0.99602386, 0.99602386,\n", " 0.99602386, 0.99602386, 1. , 1. , 1. ,\n", " 1. , 1. , 1. , 1. , 1. ,\n", " 1. , 1. ]])" ] }, "execution_count": 227, "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": 228, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 228, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sum(np.array(final_df.index == stoptimes.index)) == stoptimes.shape[0]" ] }, { "cell_type": "code", "execution_count": 229, "metadata": {}, "outputs": [], "source": [ "################################## WRITE FINAL NUMPY ##################################\n", "with gzip.open(\"../data/transfer_probabilities.pkl.gz\", \"wb\") as output_file:\n", " pickle.dump(final_np, output_file)" ] } ], "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 }