diff --git a/.gitattributes b/.gitattributes
index 4bc8eaf..f204e49 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -1,39 +1,42 @@
 data/distributions.pickle filter=lfs diff=lfs merge=lfs -text
 data/stop_times_array_version2.csv filter=lfs diff=lfs merge=lfs -text
 data/transfer_array_version2.csv filter=lfs diff=lfs merge=lfs -text
 data/routes_array_version2.csv filter=lfs diff=lfs merge=lfs -text
 data/route_stops_array_version2.csv filter=lfs diff=lfs merge=lfs -text
 data/stop_routes_array_version3.csv filter=lfs diff=lfs merge=lfs -text
 data/stops_array_version2.csv filter=lfs diff=lfs merge=lfs -text
 object.data filter=lfs diff=lfs merge=lfs -text
 *.pkl filter=lfs diff=lfs merge=lfs -text
 data/*.pkl filter=lfs diff=lfs merge=lfs -text
 data/dere.pkl filter=lfs diff=lfs merge=lfs -text
 data/transfer_array.pkl filter=lfs diff=lfs merge=lfs -text
 data/stops_array.pkl filter=lfs diff=lfs merge=lfs -text
 data/stop_times_array.pkl filter=lfs diff=lfs merge=lfs -text
 data/stop_routes_array.pkl filter=lfs diff=lfs merge=lfs -text
 data/routes_array.pkl filter=lfs diff=lfs merge=lfs -text
 data/route_stops_array.pkl filter=lfs diff=lfs merge=lfs -text
 data/route_stops_df.pkl filter=lfs diff=lfs merge=lfs -text
 data/routes_array_df.pkl filter=lfs diff=lfs merge=lfs -text
 data/stop_routes_df.pkl filter=lfs diff=lfs merge=lfs -text
 data/stop_times_df.pkl filter=lfs diff=lfs merge=lfs -text
 data/stops_df.pkl filter=lfs diff=lfs merge=lfs -text
 data/transfer_df.pkl filter=lfs diff=lfs merge=lfs -text
 data/distrib_recov_tab_stopID_hour.pkl.gz filter=lfs diff=lfs merge=lfs -text
 data/join_distribution_all.pkl.gz filter=lfs diff=lfs merge=lfs -text
 data/join_distribution_cumulative_p.pkl.gz filter=lfs diff=lfs merge=lfs -text
 data/join_distribution_cumulative_p_2.pkl.gz filter=lfs diff=lfs merge=lfs -text
 data/route_stops_array_cyril.pkl filter=lfs diff=lfs merge=lfs -text
 data/route_stops_df_cyril.pkl filter=lfs diff=lfs merge=lfs -text
 data/routes_array_cyril.pkl filter=lfs diff=lfs merge=lfs -text
 data/routes_array_df_cyril.pkl filter=lfs diff=lfs merge=lfs -text
 data/stop_routes_array_cyril.pkl filter=lfs diff=lfs merge=lfs -text
 data/stop_routes_df_cyril.pkl filter=lfs diff=lfs merge=lfs -text
 data/stop_times_array_cyril.pkl filter=lfs diff=lfs merge=lfs -text
 data/transfer_array_cyril.pkl filter=lfs diff=lfs merge=lfs -text
 data/transfer_df_cyril.pkl filter=lfs diff=lfs merge=lfs -text
 data/stops_array_cyril.pkl filter=lfs diff=lfs merge=lfs -text
 data/stop_times_df_cyril.pkl filter=lfs diff=lfs merge=lfs -text
 data/join_distribution_cumulative_p_3.pkl.gz filter=lfs diff=lfs merge=lfs -text
+data/d_all.pkl.gz filter=lfs diff=lfs merge=lfs -text
+data/d_real.pkl.gz filter=lfs diff=lfs merge=lfs -text
+data/transfer_probabilities.pkl.gz filter=lfs diff=lfs merge=lfs -text
diff --git a/notebooks/MC_RAPTOR.ipynb b/notebooks/MC_RAPTOR.ipynb
index 5f03e54..bd4d182 100644
--- a/notebooks/MC_RAPTOR.ipynb
+++ b/notebooks/MC_RAPTOR.ipynb
@@ -1,1033 +1,1033 @@
 {
  "cells": [
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "# McRAPTOR for latest departure problem \n",
     "\n",
     "*Authors: Felix Grimberg, Cyril Pulver*\n",
     "\n",
     "Multiple-criteria adaptation of RAPTOR (Round bAsed Public Transit Optimized Router) algorithm (cf. [README](../README.md)).\n",
     "\n",
     "## Left out at this stage:\n",
     "\n",
     "- Realistic time to get out of one transport and walk to the platform of the next. Instead, we just set it to 2 minutes, no matter what.\n",
     "We also reserve 2 minutes for walking from the entrance of a station to the correct platform, and vice versa, when walking between stations (total: `time_to_walk_distance + 4 minutes`)\n",
     "\n",
     "# Run the algorithm\n",
     "\n",
     "Please make sure that you follow these steps:\n",
     "1. Execute `git lfs pull` in a terminal to download our data files (if you haven't already).\n",
     "2. **Execute the code cells in the other sections of this notebook (\"Load the data and onwards\") to load the data and define the classes and functions that make up the algorithm.**\n",
     "3. Execute the two cells below\n",
     "4. Plan journeys !"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "After executing the two cells below, please enter:\n",
     "- The departure stop (only BPUIC-like stop ID are currently supported, with a minimum of 7 characters. They may also contain platform information etc.)\n",
     "- The arrival (target) stop\n",
     "- The arrival time (format HH:mm)\n",
     "- The minimum probability of success for the journey (slider)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "# Execute to set input\n",
     "grid"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "# execute to get output\n",
     "output"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "## Load the data\n",
     "### General considerations\n",
     "We adhere to the data structures proposed by Delling et al. These structures aim to minimize read times in memory by making use of consecutive in-memory adresses. Thus, structures with varying dimensions (e.g dataframes, python lists) are excluded. We illustrate the difficulty with an example. \n",
     "\n",
     "Each route has a potentially unique number of stops. Therefore, we cannot store stops in a 2D array of routes by stops, as the number of stops is not the same for each route. We adress this problem by storing stops consecutively by route, and keeping track of the index of the first stop for each route.\n",
     "\n",
     "This general strategy is applied to all the required data structures, where possible.\n",
     "\n",
     "### routes\n",
     "The `routes` array will contain arrays `[n_trips, n_stops, pt_1st_stop, pt_1st_trip]` where all four values are `int`. To avoid overcomplicating things and try to mimic pointers in python, `pt_1st_stop` and `pt_1st_trip` contain integer indices."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
     "lines_to_next_cell": 0
    },
    "outputs": [],
    "source": [
     "import numpy as np\n",
     "import pickle\n",
     "\n",
     "def pkload(path):\n",
     "    with open(path, 'rb') as f:\n",
     "        obj = pickle.load(f)\n",
     "    return obj"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
     "lines_to_next_cell": 0
    },
    "outputs": [],
    "source": [
     "routes = pkload(\"../data/routes_array_cyril.pkl\").astype(np.uint32)\n",
     "print(routes.shape)\n",
     "routes"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "### routeStops\n",
     "`routeStops` is an array that contains the ordered lists of stops for each route. `pt_1st_stop` in `routes` is required to get to the first stop of the route. is itself an array that contains the sequence of stops for route $r_i$."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "routeStops = pkload(\"../data/route_stops_array_cyril.pkl\").astype(np.uint16)\n",
     "print(routeStops.shape)\n",
     "print(routeStops.max())\n",
     "routeStops"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "### stopTimes\n",
     "\n",
     "The i-th entry in the `stopTimes` array is itself an array which contains the arrival and departure time at a particular stop for a particular trip. `stopTimes` is sorted by routes, and then by trips. We retrieve the index of the first (earliest) trip of the route with the pointer `pt_1st_trip` stored in `routes`. We may use the built-in `numpy` [date and time data structures](https://blog.finxter.com/how-to-work-with-dates-and-times-in-python/). In short, declaring dates and times is done like this: `np.datetime64('YYYY-MM-DDThh:mm')`. Entries with a `NaT` arrival or departure times correspond to beginning and end of trips respectively.\n",
     "\n",
     "Note that trips are indexed implicitely in stopTimes, but we decided to change a little bit from the paper and index them according to their parent route instead of giving them an absolute index. It makes things a bit easier when coding the algorithm."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "stopTimes = pkload(\"../data/stop_times_array_cyril.pkl\")\n",
     "print(stopTimes.shape)\n",
     "stopTimes"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "`NaT` is the `None` equivalent for `numpy datetime64`."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "print(np.isnat(stopTimes[0]))"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "### stopRoutes\n",
     "\n",
     "`stopRoutes` contains the routes (as `int`s representing an index in `routes`) associated with each stop. We need the pointer in `stops` to index `stopRoutes` correctly."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "stopRoutes = pkload(\"../data/stop_routes_array_cyril.pkl\").flatten().astype(np.uint32)\n",
     "print(stopRoutes.shape)\n",
     "stopRoutes"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "### transfers\n",
     "`transfers` is a 2D `np.ndarray` where each entry `[p_j, time]` represents (in seconds) the time it takes to walk from stop p_j to the implicitely given stop p_i.\n",
     "p_i is given implicitely by the indexing, in conjunction with `stops`. In other words:\n",
     "`transfers[stops[p_i][2]:stops[p_i][3]]` returns all the footpaths arriving at stop p_i.\n",
     "\n",
     "As we cannot store different data types in numpy arras, `time` will have to be converted to `np.timedelta64`, the format used to make differences between `np.datetime.64` variables. We will consider all `time` values as **positive values in seconds**."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "transfers = pkload(\"../data/transfer_array_cyril.pkl\").astype(np.uint16)\n",
     "print(transfers.shape)\n",
     "transfers"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "### stops\n",
     "\n",
     "`stops` stores the indices in `stopRoutes` and in `transfers` corresponding to each stop.\n",
     "\n",
     "`stopRoutes[stops[p][0]:stops[p][1]]` returns the routes serving stop p.\n",
     "\n",
     "`transfers[stops[p][2]:stops[p][3]]` returns the footpaths arriving at stop p."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "stops = pkload(\"../data/stops_array_cyril.pkl\")\n",
     "stopRoutes = pkload(\"../data/stop_routes_array_cyril.pkl\")\n",
     "print(np.isnan(stops.astype(np.float64)).sum(axis=0))\n",
     "print(np.equal(stops, None).sum(axis=0))\n",
     "print(stops.shape)\n",
     "stops = stops[:,[0,0,1,1]]\n",
     "# Make column 1 contain the start_index of the next stop in stopRoutes\n",
     "stops[:-1,1] = stops[1:,0]\n",
     "stops[-1, 1] = stopRoutes.shape[0]\n",
     "# Deal with NaN's in column 2 (for stops with 0 foot transfers within 500m)\n",
     "if np.isnan(stops[-1,2]).item():\n",
     "    stops[-1,2] = transfers.shape[0]\n",
     "for i in np.isnan(stops[:-1,2].astype(np.float64)).nonzero()[0][::-1]:\n",
     "    stops[i,2] = stops[i+1,2]\n",
     "print(np.isnan(stops.astype(np.float64)).sum(axis=0))\n",
     "# Make column 3 contain the start_index of the next stop in stopRoutes\n",
     "stops[:-1,3] = stops[1:,2]\n",
     "stops[-1, 3] = transfers.shape[0]\n",
     "# Convert to int\n",
     "stops = stops.astype(np.uint32)\n",
     "stops"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "### Example"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "p = 0\n",
     "routes_serving_p = stopRoutes[stops[p][0]:stops[p][1]]\n",
     "print(\"routes serving stop 0:\", routes_serving_p)\n",
     "for r in routes_serving_p:\n",
     "    print(\"stops of route {}:\".format(r), routeStops[routes[r][2]:routes[r][2]+routes[r][1]])\n",
     "for pPrime, walking_seconds in transfers[stops[p][2]:stops[p][3]]:\n",
     "    print(\"stop {} can be reached from stop {} by walking for {} seconds.\".format(p, pPrime, walking_seconds))"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "### Distribution of delays"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "import gzip \n",
-    "with gzip.open(\"../data/join_distribution_cumulative_p_3.pkl.gz\") as distrib_pkl:\n",
+    "with gzip.open(\"../data/transfer_probabilities.pkl.gz\") as distrib_pkl:\n",
     "    distrib_delays = pickle.load(distrib_pkl)\n",
     "    \n",
     "distrib_delays[0:2]"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "### Relate `stop_id`s and `trip_headsign`s to the integer indices used in the algorithm"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "stop_times_df = pkload(\"../data/stop_times_df_cyril.pkl\")[['route_id', 'stop_id_general', 'stop_name', 'trip_headsign', 'route_int', 'stop_int', 'route_desc']]\n",
     "stop_times_df.head()"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "stop_ids_names = stop_times_df[['stop_id_general', 'stop_int', 'stop_name']].drop_duplicates()\n",
     "assert np.all(stop_ids_names == stop_ids_names.drop_duplicates(subset='stop_int'))\n",
     "assert np.all(stop_ids_names == stop_ids_names.drop_duplicates(subset='stop_id_general'))\n",
     "assert np.all(stop_ids_names == stop_ids_names.drop_duplicates(subset='stop_name'))\n",
     "stop_ids_names = stop_ids_names.sort_values(by='stop_int')\n",
     "stop_ids = stop_ids_names['stop_id_general'].to_numpy()\n",
     "stop_names = stop_ids_names['stop_name'].to_numpy()\n",
     "print(stop_ids, stop_names)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "p = np.random.randint(stops.shape[0])\n",
     "print(p, stop_ids[p], stop_names[p])\n",
     "stop_times_df[stop_times_df['stop_int'] == p].head(2)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "# Define the classes and functions making up the journey planner\n",
     "\n",
     "Based on a modified version of RAPTOR (reversed RAPTOR), we implement a multiple criteria RAPTOR algorithm.\n",
     "The optimization criteria are:\n",
     "- Latest departure\n",
     "- Highest probability of success of the entire trip\n",
     "- Lowest number of connections (implicit with the round-based approach)\n",
     "\n",
     "We specify the time to change platforms as an absolute constant:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "tau_change_platform = np.timedelta64(2, 'm')\n",
     "np.datetime64('2020-05-11T15:30') - tau_change_platform"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "Next, we define some shorthand functions to help with the indexing in the data structures:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "def gen_random_stop_id():\n",
     "    \"\"\"Generate a random stop_id to test the journey planner.\"\"\"\n",
     "    return str(stop_ids[np.random.randint(stops.shape[0])])+':0:5DE'\n",
     "\n",
     "def stop_id_to_int(p_id):\n",
     "    \"\"\"Given a stop id, returns the corresponding stop_int\"\"\"\n",
     "    return np.asarray(stop_ids == int(p_id[:7])).nonzero()[0].item()\n",
     "\n",
     "def get_stops(r):\n",
     "    \"\"\"Returns the stops of route r\"\"\"\n",
     "    idx_first_stop = routes[r][2]\n",
     "    return routeStops[idx_first_stop:idx_first_stop+routes[r][1]] # n_stops = routes[r][1]\n",
     "\n",
     "def calc_stopTimes_idx(r, t, offset_p):\n",
     "    \"\"\"Returns the index of the entry in stopTimes\n",
     "    corresponding to the offset_p-th stop of the t-th trip\n",
     "    of route r.\n",
     "    \"\"\"\n",
     "    return (routes[r][3]       # 1st trip of route\n",
     "            + t * routes[r][1] # offset for the right trip\n",
     "            + offset_p         # offset for the right stop\n",
     "           )\n",
     "\n",
     "def get_arrival_time(r, t, offset_p):\n",
     "    \"\"\"Returns 2000 (instead of 0) if t is None.\n",
     "    Otherwise, returns the arrival time of the t-th trip of route r\n",
     "    at the offset_p-th stop of route r.\n",
     "    trips and stops of route r start at t=0, offset_p=0.\n",
     "    \"\"\"\n",
     "    if t is None:\n",
     "        return np.datetime64('2000-01-01T01:00')\n",
     "    \n",
     "    return stopTimes[calc_stopTimes_idx(r,t,offset_p)][0] # 0 for arrival time\n",
     "\n",
     "def get_departure_time(r, t, offset_p):\n",
     "    \"\"\"Throws TypeError if t is None.\n",
     "    Otherwise, returns the departure time of the t-th trip of route r\n",
     "    at the offset_p-th stop of route r.\n",
     "    trips and stops of route r start at t=0 & offset_p=0.\n",
     "    \"\"\"\n",
     "    if t is None:\n",
     "        raise TypeError(\"Requested departure time of None trip!\")\n",
     "    \n",
     "    return stopTimes[calc_stopTimes_idx(r,t,offset_p)][1] # 1 for departure time\n",
     "\n",
     "# This one isn't indexing related. Forgive us.\n",
     "def time2str(t):\n",
     "    \"\"\"Prints the hour and minutes of np.datetime64 t.\"\"\"\n",
     "    return str(t.astype('datetime64[m]')).split('T')[1]"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "In the cell below, we define the following class hierarchy of label types:\n",
     "```\n",
     "BaseLabel\n",
     " +-- RouteLabel\n",
     " +-- ImmutableLabel\n",
     "      +-- WalkLabel\n",
     "      +-- TargetLabel\n",
     "```\n",
     "Labels come together recursively to represent a (partial) solution, going from the `stop` attribute of the first label to the target stop (that stop will be referred to as \"here\").\n",
     "Each partial solution is associated with a departure time `tau_dep`, an overall probability of success (all the way to the target stop) `Pr`, and a number of trips necessary to get from here to the target stop.\n",
     "Any two labels can be compared via their `dominates` method, to retain only non-dominated labels.\n",
     "Further, all types of labels implement a `print_instructions` method, used to recursively print instructions for the entire journey.\n",
     "Each `RouteLabel` corresponds to the portion of a journey that begins with standing on the correct platform to board a vehicle, and ends with exiting from that vehicle at another stop.\n",
     "Each `WalkLabel` corresponds to the portion of a journey that begins with walking away from the main hall of one stop and ends with having walked all the way to the main hall of another stop."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "class InstantiationException(Exception):\n",
     "    pass\n",
     "\n",
     "class BaseLabel:\n",
     "    \"\"\"An abstract base class for Labels. Do not instantiate.\n",
     "    A label corresponds to a recursive (partial) solution, going\n",
     "    to the target stop from the stop currently under consideration.\n",
     "    \"\"\"\n",
     "    indent = \" \"*4\n",
     "    def __init__(self, stop, tau_dep, Pr, n_trips):\n",
     "        self.stop    = stop\n",
     "        self.tau_dep = tau_dep\n",
     "        self.Pr      = Pr\n",
     "        self.n_trips = n_trips\n",
     "    \n",
     "    def dominates(self, other):\n",
     "        \"\"\"Returns True if self dominates other, else returns False.\n",
     "        other: another Label instance.\n",
     "        \"\"\"\n",
     "        if  self.tau_dep >= other.tau_dep \\\n",
     "        and self.Pr      >= other.Pr      \\\n",
     "        and self.n_trips <= other.n_trips :\n",
     "            return True\n",
     "        return False\n",
     "    \n",
     "    def stop2str(self):\n",
     "        \"\"\"Returns a printable str representing self.stop\"\"\"\n",
     "        return \"{stopN} (stop {stopI})\".format(\n",
     "                    stopN = stop_names[self.stop],\n",
     "                    stopI = stop_ids[self.stop]\n",
     "        )\n",
     "        \n",
     "    def print_journey(self):\n",
     "        print(\"Departure stop:        \", self.stop2str())\n",
     "        print(\"Departure time:        \", time2str(self.tau_dep))\n",
     "        print(\"Number of trips used:  \", self.n_trips)\n",
     "        print(\"Probability of success: {:.4f}\".format(self.Pr))\n",
     "        \n",
     "        self.print_instructions()\n",
     "    \n",
     "    def to_str(self):\n",
     "        s = \"Departure at {0} from stop {1} (id: {2}, int: {3}).\".format(\n",
     "            self.tau_dep,\n",
     "            stop_names[self.stop],\n",
     "            stop_ids[self.stop],\n",
     "            self.stop\n",
     "        )\n",
     "        return repr(type(self)) + s\n",
     "    \n",
     "    def pprint(self, indent=''):\n",
     "        print(indent, self.to_str())\n",
     "    \n",
     "    def copy(self):\n",
     "        raise InstantiationException(\"class BaseLabel should never \"\n",
     "                                     \"be instantiated.\"\n",
     "                                    )\n",
     "\n",
     "class ImmutableLabel(BaseLabel):\n",
     "    \"\"\"Base class for immutable Labels\"\"\"\n",
     "    def copy(self):\n",
     "        return self\n",
     "\n",
     "class TargetLabel(ImmutableLabel):\n",
     "    \"\"\"A special type of label reserved for the target stop.\"\"\"\n",
     "    def __init__(self, stop, tau_dep):\n",
     "        BaseLabel.__init__(self, stop, tau_dep, 1., 0)\n",
     "        \n",
     "    def print_instructions(self):\n",
     "        \"\"\"Finish printing instructions for the journey.\"\"\"\n",
     "        print(\"Target stop:           \", self.stop2str())\n",
     "        print(\"Requested arrival time:\", time2str(self.tau_dep))\n",
     "\n",
     "class WalkLabel(ImmutableLabel):\n",
     "    \"\"\"A special type of label for walking connections.\"\"\"\n",
     "    def __init__(self, stop, tau_walk, next_label):\n",
     "        \"\"\"Create a new WalkLabel instance.\n",
     "        stop:       stop where you start walking.\n",
     "        tau_walk:   (np.timedelta64) duration of the walk.\n",
     "        next_label: label describing the rest of the trip after walking.\n",
     "        \"\"\"\n",
     "        if isinstance(next_label, WalkLabel):\n",
     "            raise ValueError(\"Cannot chain two consecutive WalkLabels!\")\n",
     "        tau_dep = next_label.tau_dep - tau_walk - tau_change_platform\n",
     "        BaseLabel.__init__(self, stop, tau_dep, next_label.Pr, next_label.n_trips)\n",
     "        self.tau_walk   = tau_walk\n",
     "        self.next_label = next_label\n",
     "        \n",
     "    def print_instructions(self):\n",
     "        \"\"\"Recursively print instructions for the whole journey.\"\"\"\n",
     "        print(self.indent + \"Walk {:.1f} minutes from\".format(\n",
     "                self.tau_walk / np.timedelta64(1,'m')\n",
     "                ),\n",
     "              self.stop2str()\n",
     "             )\n",
     "        print(self.indent*2 + \"to\", self.next_label.stop2str())\n",
     "        self.next_label.print_instructions()\n",
     "\n",
     "class RouteLabel(BaseLabel):\n",
     "    \"\"\"A type of label for regular transports.\"\"\"\n",
     "    def __init__(self,\n",
     "                 tau_dep,\n",
     "                 r,\n",
     "                 t,\n",
     "                 offset_p,\n",
     "                 next_label,\n",
     "                 Pr_connection_success\n",
     "                ):\n",
     "        \n",
     "        self.tau_dep      = tau_dep\n",
     "        self.r            = r\n",
     "        self.t            = t\n",
     "        self.offset_p_in  = offset_p\n",
     "        self.offset_p_out = offset_p\n",
     "        self.next_label   = next_label\n",
     "        # Store Pr_connection_success for self.copy()\n",
     "        self.Pr_connection_success = Pr_connection_success\n",
     "        \n",
     "        self.route_stops = get_stops(self.r)\n",
     "        self.stop        = self.route_stops[self.offset_p_in]\n",
     "        self.Pr          = self.Pr_connection_success * self.next_label.Pr\n",
     "        self.n_trips     = self.next_label.n_trips + 1\n",
     "    \n",
     "    def update_stop(self, stop):\n",
     "        self.stop     = stop\n",
     "        self.offset_p_in = self.offset_p_in - 1\n",
     "        # Sanity check:\n",
     "        assert self.offset_p_in >= 0\n",
     "        assert self.route_stops[self.offset_p_in] == stop\n",
     "        self.tau_dep  = get_departure_time(self.r, self.t, self.offset_p_in)\n",
     "        \n",
     "    def print_instructions(self):\n",
     "        \"\"\"Recursively print instructions for the whole journey.\"\"\"\n",
     "        stopTimes_idx = calc_stopTimes_idx(self.r, self.t,\n",
     "                                           self.offset_p_in)\n",
     "        \n",
     "        print(self.indent   + \"At\", self.stop2str())\n",
     "        print(self.indent*2 + \"take the\",\n",
     "              stop_times_df['route_desc'][stopTimes_idx], \"to\",\n",
     "              stop_times_df['trip_headsign'][stopTimes_idx]\n",
     "             )\n",
     "        print(self.indent*2 + \"leaving at\", time2str(self.tau_dep),\n",
     "              \"(route id: {}).\".format(\n",
     "                  stop_times_df['route_id'][stopTimes_idx]\n",
     "                  )\n",
     "             )\n",
     "        \n",
     "        tau_arr = get_arrival_time(\n",
     "                    self.r,\n",
     "                    self.t,\n",
     "                    self.offset_p_out\n",
     "                    )\n",
     "        assert self.next_label.stop == self.route_stops[self.offset_p_out]\n",
     "        \n",
     "        print(self.indent   + \"Get out at\", self.next_label.stop2str())\n",
     "        print(self.indent*2 + \"at\", time2str(tau_arr)+\".\")\n",
     "\n",
     "        self.next_label.print_instructions()\n",
     "        \n",
     "    def copy(self):\n",
     "        \"\"\"When RouteLabels are merged into the bag of a stop,\n",
     "        they must be copied (because they will subsequently\n",
     "        be changed with self.update_stop()).\n",
     "        \"\"\"\n",
     "        l = RouteLabel(self.tau_dep,\n",
     "                       self.r,\n",
     "                       self.t,\n",
     "                       self.offset_p_in,\n",
     "                       self.next_label,\n",
     "                       self.Pr_connection_success\n",
     "                      )\n",
     "        l.offset_p_out = self.offset_p_out\n",
     "        return l"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "The functions below define the bag merge operation and implement the journey planner:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "def run_mc_raptor(stop_id_start,\n",
     "                  stop_id_destination,\n",
     "                  arrival_time_max,\n",
     "                  Pr_min\n",
     "                 ):\n",
     "    \"\"\"Run MC RAPTOR, using the data defined in cells above (stopRoutes etc.).\n",
     "    Inputs:\n",
     "        stop_id_start:       source stop id\n",
     "                             str (format: '1234567')\n",
     "        stop_id_destination: target stop id\n",
     "                             str (format: '1234567')\n",
     "        arrival_time_max:    latest acceptable arrival time\n",
     "                             str (format: '14:56')\n",
     "        Pr_min:              minimum acceptable probability of success\n",
     "                             float in open interval (0,1)\n",
     "    Outputs:\n",
     "        bags_p_s: bags_p_s[k] contains the Pareto set of non-dominated journeys\n",
     "                  from stop_id_start to stop_id_destination that use at most\n",
     "                  k different trips (i.e, getting in at most k different vehicles),\n",
     "                  under the given constraints:\n",
     "                      1. Each journey must succeed with a probability\n",
     "                         greater or equal to Pr_min.\n",
     "                      2. The journey is a success if and only if all individual\n",
     "                         connections succeed, including the user's appointment\n",
     "                         in p_t at tau_0.\n",
     "                      3. A connection succeeds if, and only if, the user reaches\n",
     "                         the platform before or on the scheduled departure time\n",
     "                         (allowing some time to change platforms)\n",
     "                  Non-dominated:\n",
     "                      A journey J1 is *dominated* by another journey J2, if\n",
     "                      J2 departs no earlier than J1 AND the probability of\n",
     "                      success of J2 is no less than that of J1.\n",
     "                  Pareto set:\n",
     "                      Each bag in bags_p_s contains only journeys that are not\n",
     "                      dominated by any other possible journey. Such a collection\n",
     "                      of non-dominated solutions is called a *Pareto set*.\n",
     "        bags:     bags[k][stop_id_to_int(stop_id)] contains the Pareto set of\n",
     "                  non-dominated journeys from stop_id to stop_id_destination\n",
     "                  that use at most k different trips under the given constraints.\n",
     "                  \n",
     "    Each journey is represented as a Label that forms the start of a chain.\n",
     "    The journey can be reconstructed by calling label.print_journey().\n",
     "    \"\"\"\n",
     "# Input sanitization:\n",
     "    try:\n",
     "        tau_0 = np.datetime64('2020-05-24T'+arrival_time_max)\n",
     "        p_s   = stop_id_to_int(stop_id_start)\n",
     "        p_t   = stop_id_to_int(stop_id_destination)\n",
     "    except Exception as e:\n",
     "        print(\"ERROR parsing input! Please make sure that:\")\n",
     "        print(\"- stop_id_start and stop_id_destination are strings \"\n",
     "              \"beginning with 7 or more digits.\")\n",
     "        print(\"- arrival_time_max is in the format 'hh:mm'. Example: '09:42'.\")\n",
     "        raise e\n",
     "    \n",
     "# initialization\n",
     "    k_max = 10 # Maximum number of rounds\n",
     "    \n",
     "    # For each route and for each label at each stop p, we will look at the n latest\n",
     "    # trips until we find a trip for which the individual connection at stop p\n",
     "    # succeeds with a probability at least equal to Pr_threshold.\n",
     "    # Under some reasonable assumptions, setting Pr_threshold = Pr_min**(1/k)\n",
     "    # guarantees that we will find a solution, if a solution exists involving at\n",
     "    # most k connections (including the user's appointment in p_t at tau_0).\n",
     "    Pr_threshold = Pr_min**(0.1)\n",
     "    \n",
     "    # Initialize empty bags for each stop for round 0:\n",
     "    n_stops = stops.shape[0]\n",
     "    bags = [\n",
     "        [\n",
     "            [] # an empty bag\n",
     "        for _ in range(n_stops)] # one empty bag per stop\n",
     "    ]\n",
     "\n",
     "    # Create a TargetLabel for p_t, and mark p_t\n",
     "    bags[0][p_t].append(TargetLabel(p_t, tau_0))\n",
     "    marked = {p_t}\n",
     "\n",
     "    print(\"Searching for journeys from\",\n",
     "          stop_names[p_s],\n",
     "          \"(stop {})\".format(stop_ids[p_s]),\n",
     "          \"to\", bags[0][p_t][0].stop2str())\n",
     "\n",
     "# Define bag operations (they depend on p_s and Pr_min for target pruning):\n",
     "    def update_bag(bag, label, k, target_pruning=True):\n",
     "        \"\"\"Add label to bag and remove dominated labels.\n",
     "        bag is altered in-place.\n",
     "\n",
     "        k: Round number, used for target pruning.\n",
     "\n",
     "        returns: Boolean indicating whether bag was altered.\n",
     "        \"\"\"\n",
     "        if target_pruning:\n",
     "            # Apply the Pr_min constraint to label:\n",
     "            if label.Pr < Pr_min:\n",
     "                return False\n",
     "\n",
     "            # Prune label if it is dominated by bags[k][p_s]:\n",
     "            for L_star in bags[k][p_s]:\n",
     "                if L_star.dominates(label):\n",
     "                    return False\n",
     "\n",
     "        # Otherwise, merge label into bag1\n",
     "        changed = False\n",
     "        for L_old in bag:\n",
     "            if L_old.dominates(label):\n",
     "                return changed\n",
     "            if label.dominates(L_old):\n",
     "                bag.remove(L_old)\n",
     "                changed = True\n",
     "        bag.append(label.copy())\n",
     "        return True\n",
     "\n",
     "    def merge_bags(bag1, bag2, k, target_pruning=True):\n",
     "        \"\"\"Merge bag2 into bag1 in-place.\n",
     "        k: Round number, used for target pruning.\n",
     "        returns: Boolean indicating whether bag was altered.\n",
     "        \"\"\"\n",
     "        changed = False\n",
     "        for label in bag2:\n",
     "            changed = update_bag(bag1, label, k, target_pruning) or changed\n",
     "        return changed\n",
     "    \n",
     "    globals().update({'merge_bags': merge_bags})\n",
     "    \n",
     "# Define the footpaths-checking function (depends on update_bag)\n",
     "    def check_footpaths(bags, marked, k):\n",
     "        \"\"\"Modify bags and marked in-place to account for foot-paths.\"\"\"\n",
     "        q = []\n",
     "        for p in marked:\n",
     "            for pPrime, delta_seconds in transfers[stops[p][2]:stops[p][3]]:\n",
     "                q.append((p, pPrime, delta_seconds))\n",
     "        for p, pPrime, delta_seconds in q:\n",
     "            for L_k in bags[k][p]:\n",
     "                # We do not allow two consecutive walking trips\n",
     "                if not isinstance(L_k, WalkLabel):\n",
     "                    L_new = WalkLabel(pPrime, np.timedelta64(delta_seconds, 's'), L_k)\n",
     "                    if update_bag(bags[k][pPrime], L_new, k):\n",
     "                        marked.add(pPrime)\n",
     "\n",
     "# main loop\n",
     "    indent= ' '*4\n",
     "\n",
     "    k = 0\n",
     "    # Check footpaths leading to p_t at k=0:\n",
     "    check_footpaths(bags, marked, k)\n",
     "    while k < k_max:\n",
     "        k += 1 # k=1 at fist round, as it should.\n",
     "\n",
     "        # Instead of using best bags, carry over the bags from last round.\n",
     "    #     if len(bags <= k):\n",
     "\n",
     "        bags.append([bags[-1][p].copy() for p in range(n_stops)])\n",
     "\n",
     "        print('\\n', ' '*30, 'STARTING round k =', k)\n",
     "        if len(marked) < 50:\n",
     "            print('Marked stops at the start of the round:', [stop_ids[p] for p in marked])\n",
     "        else:\n",
     "            print('Number of marked stops at the start of the round:', len(marked))\n",
     "        \n",
     "        # accumulate routes serving marked stops from previous rounds\n",
     "        q = []\n",
     "        for p in marked:\n",
     "            for r in stopRoutes[stops[p][0]:stops[p][1]]: # foreach route r serving p\n",
     "                append_r_p = True\n",
     "                for idx, (rPrime, pPrime) in enumerate(q): # is there already a stop from the same route in q ?\n",
     "                    if rPrime == r:\n",
     "                        append_r_p = False\n",
     "                        p_pos_in_r = np.where(get_stops(r) == p)[0][-1]\n",
     "                        pPrime_pos_in_r = np.where(get_stops(r) == pPrime)[0][-1]\n",
     "                        if p_pos_in_r > pPrime_pos_in_r:\n",
     "                            q[idx] = (r, p) # substituting (rPrime, pPrime) by (r, p)\n",
     "                if append_r_p:\n",
     "                    q.append((r, p))\n",
     "        marked.clear() # unmarking all stops\n",
     "#         print(\"Queue:\", q)\n",
     "\n",
     "#         print('Queue before traversing each route: {}'.format(q))\n",
     "        # traverse each route\n",
     "        for (r, p) in q:\n",
     "#             print('\\n****TRAVERSING ROUTE r={0} from stop p={1}****'.format(r, p))\n",
     "            B_route = [] # new empty route bag\n",
     "\n",
     "            # we traverse the route backwards (starting at p, not from the end of the route)\n",
     "            stops_of_current_route = get_stops(r)\n",
     "#             print('Stops of current route:', stops_of_current_route)\n",
     "            offset_p = np.asarray(stops_of_current_route == p).nonzero()[0]\n",
     "            if offset_p.size < 1:\n",
     "                print(\"WARNING: route {r} is said to serve stop {p} in stopRoutes, but stop {p} \"\n",
     "                      \"is not included as a stop of route {r} in routeStops...\".format(p=p, r=r))\n",
     "                offset_p = -1\n",
     "            else:\n",
     "                offset_p = offset_p[-1]\n",
     "            for offset_p_i in range(offset_p, -1, -1):\n",
     "                p_i = stops_of_current_route[offset_p_i]\n",
     "#                 print('\\n\\n'+indent+\"p_i: {}\".format(p_i))\n",
     "\n",
     "                # Update the labels of the route bag:\n",
     "                for L in B_route:\n",
     "                    L.update_stop(p_i)\n",
     "\n",
     "                # Merge B_route into bags[k][p_i]\n",
     "                if merge_bags(bags[k][p_i], B_route, k):\n",
     "#                     print(\"marking stop\", p_i)\n",
     "                    marked.add(p_i)\n",
     "\n",
     "                # Can we step out of a later trip at p_i ?\n",
     "                # This is only possible if we already know a way to get from p_i to p_t in < k vehicles\n",
     "                # (i.e., if there is at least one label in bags[k-1][p_i])\n",
     "                for L_k in bags[k-1][p_i]:\n",
     "                    # Note that k starts at 1 and bags[0][p_t] contains a TargetLabel.\n",
     "#                     print('\\n'+indent+'----scanning arrival times for route r={0} at stop p_i={1}----'.format(r, p_i))\n",
     "\n",
     "                    # We check the trips from latest to earliest\n",
     "                    for t in range(routes[r][0]-1, -1, -1): # n_trips = routes[r][0]\n",
     "                        # Does t_r arrive early enough for us to make the rest \n",
     "                        # of the journey from here (tau[k-1][p_i])?\n",
     "                        tau_arr = get_arrival_time(r, t, offset_p_i)\n",
     "#                         print(indent+'arrival time: ', tau_arr)\n",
     "                        if tau_arr <= L_k.tau_dep - tau_change_platform:\n",
     "\n",
     "                            max_delay = L_k.tau_dep - tau_arr - tau_change_platform\n",
     "                            max_delay_int = min(max_delay.astype('timedelta64[m]').astype('int'), 30)\n",
     "                            \n",
     "                            Pr_connection = distrib_delays[calc_stopTimes_idx(r, t, offset_p_i),\n",
     "                                                           max_delay_int + 1]\n",
     "#                             print(Pr_connection)\n",
     "                            L_new = RouteLabel(get_departure_time(r, t, offset_p_i),\n",
     "                                               r,\n",
     "                                               t,\n",
     "                                               offset_p_i,\n",
     "                                               L_k,\n",
     "                                               Pr_connection\n",
     "                                              )\n",
     "                            update_bag(B_route, L_new, k)#:\n",
     "#                                 print(indent+\"Explored connection from\")\n",
     "#                                 L_new.pprint(indent*2)\n",
     "#                                 print(indent+\"to\")\n",
     "#                                 L_k.pprint(indent*2)\n",
     "\n",
     "                            # We don't want to add a label for every trip that's earlier than tau_dep.\n",
     "                            # Instead, we stop once we've found a trip that's safe enough.\n",
     "                            if Pr_connection > Pr_threshold:\n",
     "                                break\n",
     "        \n",
     "#         print(marked)\n",
     "        # Look at foot-paths (bags and marked are altered in-place):\n",
     "        check_footpaths(bags, marked, k)\n",
     "#         print(marked)\n",
     "        \n",
     "        # Report the number of new journeys found this round:\n",
     "        n_new_journeys = len(bags[k][p_s]) - len(bags[k-1][p_s])\n",
     "        if n_new_journeys:\n",
     "            print(\"Found\", n_new_journeys, \"non-dominated journeys with\", k, \"trips.\")\n",
     "            # If there's a journey with X trips, we won't look for a solution with X+3 trips:\n",
     "            if not bags[k-1][p_s]:\n",
     "                # The condition above means we found a journey for the first time.\n",
     "                print(\"Will search for journeys with up to\", k+1, \"trips.\")\n",
     "                k_max = k+1\n",
     "\n",
     "        # Additional stopping criterion: reached equilibrium\n",
     "        if not marked:\n",
     "            print(\"\\n\", ' '*15, \"*\"*15, \" THE END \", \"*\"*15)\n",
     "            print(\"Equilibrium reached with\" + \"out\"*(not bags[k][p_s]),\n",
     "                  \"finding a solution. The end.\")\n",
     "            return [bags[K][p_s] for K in range(len(bags))], bags\n",
     "    \n",
     "    # We have exited the while-loop because k==kmax:\n",
     "    print(\"\\n\" + \"*\"*15 + \" THE END \" + \"*\"*15)\n",
     "    if bags[k][p_s]:\n",
     "        print(\"There is a solution with\", k-1, \"trips.\")\n",
     "    else:\n",
     "        print(\"There are no solutions with up to\", k, \"trips.\")\n",
     "    print(\"We shall not search for journeys with\", k+1, \"or more trips.\")\n",
     "\n",
     "    return [bags[K][p_s] for K in range(len(bags))], bags\n",
     "\n",
     "def time_sorted_pareto(bags_p):\n",
     "    \"\"\"Input: list of bags, e.g. for one stop and various k.\n",
     "              It is assumed that Pr >= Pr_min for each label.\n",
     "    Output: A Pareto set of non-dominated labels (np.array),\n",
     "            sorted by decreasing departure time.\n",
     "    \"\"\"\n",
     "    res_bag = []\n",
     "    for bag in bags_p:\n",
     "        globals()['merge_bags'](res_bag, bag, 0, target_pruning=False)\n",
     "    res = np.array(res_bag)\n",
     "    return res[np.argsort([label.tau_dep for label in res])[::-1]]\n",
     "\n",
     "def print_solutions(bags_p):\n",
     "    print(\"Showing Pareto set of non-dominated journeys\")\n",
     "    if not bags_p[-1]:\n",
     "        print(\"There are no journeys to print.\")\n",
     "        return\n",
     "    L_s = bags_p[-1][0]\n",
     "    L_t = L_s\n",
     "    while not isinstance(L_t, TargetLabel):\n",
     "        L_t = L_t.next_label\n",
     "    print(\"from\", L_s.stop2str(), \"to\", L_t.stop2str())\n",
     "    print(\"with criteria:\")\n",
     "    print(\"  - maximize departure time\")\n",
     "    print(\"  - maximize probability of success\")\n",
     "    print(\"  - minimize number of individual trips\")\n",
     "    print(\"and constraints:\")\n",
     "    print(\"  - arrival time at least\", tau_change_platform,\n",
     "              \"minutes before\", time2str(L_t.tau_dep), \"(2 minutes to walk from the platform to the extraction point).\")\n",
     "    print(\"  - probability of success at least {Pr_min:.4f}.\".format(**globals()))\n",
     "    print(\"\\nJourneys are sorted by descending departure time.\")\n",
     "    \n",
     "    for i, label in enumerate(time_sorted_pareto(bags_p)):\n",
     "        print('\\n'*2, '-'*10, 'OPTION', i+1, '-'*10)\n",
     "        label.print_journey()"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "# basic user interface\n",
-    "Pr_min = 0.9\n",
     "\n",
     "import ipywidgets as widgets\n",
+    "\n",
     "departure_stop_input = widgets.Text(description='departure stop ID')\n",
     "arrival_stop_input = widgets.Text(description='arrival stop ID')\n",
     "arrival_time_input = widgets.Text(description = 'arrival time (HH:mm)')\n",
     "\n",
     "probability_min_input = widgets.FloatSlider(\n",
     "    value=0.9,\n",
     "    min=0,\n",
     "    max=1.0,\n",
     "    step=0.01,\n",
     "    description='Min. probability of success:',\n",
     "    disabled=False,\n",
     "    continuous_update=False,\n",
     "    orientation='horizontal',\n",
     "    readout=True,\n",
     "    readout_format='.2f',\n",
     ")\n",
     "\n",
     "from IPython.display import display\n",
     "find_journey_button = widgets.Button(description=\"Find journey\")\n",
     "output = widgets.Output()\n",
     "\n",
     "def on_find_journey_button_clicked(b):\n",
     "    with output:\n",
     "        output.clear_output()\n",
     "        bags_p_s, bags = run_mc_raptor(departure_stop_input.value,\n",
     "                               arrival_stop_input.value,\n",
     "                               arrival_time_input.value,\n",
     "                               probability_min_input.value,\n",
     "                                      )\n",
     "        print_solutions(bags_p_s)\n",
     "        \n",
     "find_journey_button.on_click(on_find_journey_button_clicked)\n",
     "\n",
     "from ipywidgets import GridspecLayout\n",
     "grid = GridspecLayout(3, 2)\n",
     "grid[0, 0] = departure_stop_input\n",
     "grid[0, 1] = arrival_stop_input\n",
     "grid[1, 0] = arrival_time_input\n",
     "grid[1, 1] = probability_min_input\n",
     "grid[2, :] = find_journey_button"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     " "
    ]
   }
  ],
  "metadata": {
   "jupytext": {
    "formats": "ipynb,md,py:percent"
   },
   "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
 }
diff --git a/notebooks/gen_transfer_proba.ipynb b/notebooks/gen_transfer_proba.ipynb
index cae65d8..67f464b 100644
--- a/notebooks/gen_transfer_proba.ipynb
+++ b/notebooks/gen_transfer_proba.ipynb
@@ -1,890 +1,921 @@
 {
  "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": [
        "<Figure size 720x360 with 1 Axes>"
       ]
      },
      "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": [
        "<Figure size 576x432 with 1 Axes>"
       ]
      },
      "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": null,
+   "execution_count": 221,
    "metadata": {},
-   "outputs": [],
+   "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": null,
+   "execution_count": 224,
    "metadata": {},
-   "outputs": [],
+   "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": [
+       "<Figure size 576x432 with 1 Axes>"
+      ]
+     },
+     "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",
-    "\n",
     "plt.show()"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 33,
+   "execution_count": 225,
    "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>[16, 497, 62, 14, 2, 3, 4, 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>[2, 512, 58, 18, 3, 2, 2, 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, 417, 69, 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, 405, 80, 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, 389, 81, 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  [16, 497, 62, 14, 2, 3, 4, 0, 0, 1, 0, 0, 0, 0...  d_all  \n",
-       "1  [2, 512, 58, 18, 3, 2, 2, 1, 0, 0, 1, 0, 0, 0,...  d_all  \n",
-       "2  [0, 417, 69, 11, 2, 0, 0, 1, 0, 0, 1, 0, 0, 0,...  d_all  \n",
-       "3  [0, 405, 80, 12, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0,...  d_all  \n",
-       "4  [0, 389, 81, 24, 4, 1, 0, 1, 0, 0, 0, 0, 1, 0,...  d_all  "
+       "                                        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": 33,
+     "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": 354,
+   "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": 355,
+   "execution_count": 227,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "array([[0.03 , 0.92 , 0.982, 0.99 , 0.992, 0.992, 0.994, 0.994, 0.994,\n",
-       "        0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996,\n",
-       "        0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 1.   , 1.   , 1.   ,\n",
-       "        1.   , 1.   , 1.   , 1.   , 1.   ],\n",
-       "       [0.002, 0.886, 0.968, 0.99 , 0.992, 0.992, 0.992, 0.994, 0.994,\n",
-       "        0.994, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996,\n",
-       "        0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 1.   , 1.   , 1.   ,\n",
-       "        1.   , 1.   , 1.   , 1.   , 1.   ],\n",
-       "       [0.   , 0.83 , 0.966, 0.988, 0.992, 0.992, 0.992, 0.994, 0.994,\n",
-       "        0.994, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996,\n",
-       "        0.996, 0.996, 0.996, 0.996, 0.996, 1.   , 1.   , 1.   , 1.   ,\n",
-       "        1.   , 1.   , 1.   , 1.   , 1.   ],\n",
-       "       [0.   , 0.806, 0.964, 0.988, 0.99 , 0.992, 0.992, 0.994, 0.994,\n",
-       "        0.994, 0.994, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996,\n",
-       "        0.996, 0.996, 0.996, 0.996, 1.   , 1.   , 1.   , 1.   , 1.   ,\n",
-       "        1.   , 1.   , 1.   , 1.   , 1.   ],\n",
-       "       [0.   , 0.776, 0.934, 0.982, 0.99 , 0.992, 0.992, 0.994, 0.994,\n",
-       "        0.994, 0.994, 0.994, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996,\n",
-       "        0.996, 0.996, 0.996, 0.996, 1.   , 1.   , 1.   , 1.   , 1.   ,\n",
-       "        1.   , 1.   , 1.   , 1.   , 1.   ]])"
+       "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": 355,
+     "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": 243,
+   "execution_count": 228,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
        "True"
       ]
      },
-     "execution_count": 243,
+     "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": 244,
+   "execution_count": 229,
    "metadata": {},
    "outputs": [],
    "source": [
     "################################## WRITE FINAL NUMPY ##################################\n",
-    "with gzip.open(\"../data/join_distribution_cumulative_p_3.pkl.gz\", \"wb\") as output_file:\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
 }