step_2_filter_mapmatch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "The primary goal of this notebook is to obtain the map matching for GPS data. To accomplish this goal, we are going to parse the trajectories, analyze the travel speed and duration as well as the departures of each trip"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import osmnx as ox\n",
    "import networkx as nx\n",
    "from collections import Counter\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "#import seaborn as sns\n",
    "\n",
    "import os\n",
    "import glob\n",
    "import math\n",
    "import numpy as np\n",
    "\n",
    "import calendar\n",
    "from datetime import datetime\n",
    "import time\n",
    "\n",
    "import requests\n",
    "import json "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data from the parsed trajectories\n",
    "data=pd.read_csv('data_parsed.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = #add an id\n",
    "test=data[data['activity_id']==idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This image shows an example fo why we need to parse trajectories before using the data. We took a trip registered by the app and parsed it. We can see how a single trip is divided into 6 other trips/tours in which the latitude of the person changes over the time. Also, we can see how points classified as `stays` have a similar latitutude over the time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VOXVwPHfAcIW1oSIIGggskgkBAjKoiCgIIoIKuKCIta6VH1V6kZrlVLUWrUutVWpfRGLviqlokWpIIKAKBoUkEX2CBGFmAhCANnO+8dzJwwhkwzJzNws5/v5zGdm7nruMMzJfZ57zyOqijHGGFNa1fwOwBhjTMVmicQYY0yZWCIxxhhTJpZIjDHGlIklEmOMMWViicQYY0yZWCIxFY6IXCciC6OwXRGRSSLyo4h8FuntlxcikiUi50Zx+1eLyKxobd+UP5ZIzFG8H5m9IrI76PGc33HFyFnAeUALVT3D72D8ICIvi8iEsmxDVV9V1QFljCNZRFREapRlOyY27B/JFOUiVf3A7yB8cAqQpar5fgcSTERqqOpBv+OA8hVLcSpKnJWFnZGYsInI8yLyr6D3j4nIHK9JqLGIzBCRHK9paIaItAhadp6ITBCRRd5Zzn9EJFFEXhWRn0TkcxFJDlpeReR/RGSjiPwgIo+LSJHfVxFpLyKzRSRPRNaIyOXFHENzEXnHW3a9iPzSm/4L4CWghxff74tYN0VEPhSRXC+mV0WkUTH7CnkMJW3LOzO8T0SWA/kiUsN7/62I7PKOs7+3bDURuV9ENnjbe1NEEoK2dY2IfOPN+20x8d4IXA3cG/g3KiaWwP52icgqERkWtJ2jmh6L+/cRkToi8qQX304RWSgidYD53iI7vFh6eMf5gLfsdhF5RUQaetsJnMH8QkQ2Ax+KyLsicnuhY1wuIkNDfQamlFTVHvYoeABZwLkh5tUF1gLXAWcDP+CagQASgUu9ZeoDU4HpQevOA9YDKUBDYJW3rXNxZ8avAJOClldgLpAAnOwte4M37zpgofc6HtgCjPa208WLKzXEMXwE/A2oDaQDOUD/wtsNse6puKavWkAS7sfu6WKWL+4Yit2W9++wFGgJ1AHaecfZ3JufDKR4r+8EPgVaeNt7Efg/b14HYDfQ25v3Z+BgMf/GLwMTivhOFMTiTRsONMf9MToCyAeaHe+/D/BX77txElAd6OnFmex9fjWC4rge9x1qDdQD/g38M+jzUNz3KN77zC4HFget3wnIBWr6/f+ssj18D8Ae5evh/WjsBnYEPX4ZNP8MIA/4BriymO2kAz8GvZ8H/Dbo/ZPAzKD3FwFLg94rcH7Q+18Bc7zXwT9UI4AFhfb9IvBQETG1BA4B9YOmPQq8XHi7YX5WQ4Evi5kf8hhK2pb373B90PtTge24xBtXaN3VeMnQe98MOOD9cD8IvB40Lx7Yz/EnkutDHae3zFLg4uP598Elob1ApyK2F0gMwYlkDvCroPftgo4zsHzroPm1vO9qG+/9E8Df/P4/Vhkf1kdiijJUQ/SRqOpnIrIROAF4MzBdROoCTwHnA429yfVFpLqqHvLebwva1N4i3tcrtLstQa+/wf0FXNgpwJkisiNoWg3gn0Us2xzIU9VdhbabUcSyxxCRE4BncWdj9XE/hD+WsFqRxxDmtgrWVdX1InInMA5IFZH3gTGquhX3GbwlIoeD1j0ENPX2F7ydfBHJDed4izkORORaYAzuBxzcv12TItYr7t+nCe7McEOYMTTHfYYB33jbalpUnKr6s4i8CYz0miqvBC4Lc1/mOFgfiTkuInIr7i+9rcC9QbN+jfsL8UxVbYBrSgGQMuyuZdDrk719FrYF+EhVGwU96qnqLUUsuxVIEJH6hbb7bZjxPIr7qzfNO8aRlHx8oY4hnG0dVZpbVV9T1bNwP84KPObN2gIMKvQZ1FbVb4HvgmPwEn5iMfGGKgdeMF1ETgH+DtwGJKpqI2BFEfEHYgv17/MDsA/X3BlOHIGkGXAyrpku+A+SwutNxvX79Af2qOonIY7PlIElEhM2EWkLTMD96F2D65RN92bXx51V7PA6eh+KwC7vEdeJ3xK4A3ijiGVmAG29DuU479FNRE4rvKCqbgEWAY+KSG0RSQN+AbwaZjz18Zr9ROQk4J4yHMNxbUtE2olIPxGphfvx3Ys76wB4AXjY+4FHRJJE5GJv3r+AwSJylojUBMZT/P/7bbg+iOLE436wc7z9jQZOD7FsyH8fVT0M/C/wZ3EXQVT3OtVreds+XCiW/wPuEpFWIlIPeAR4Q4u5OstLHIdxTalFnaWaCLBEYoryHzn6PpK3xF3PPwV4TFWXqeo64DfAP73/+E/jOjh/wHX8/jcCcbwNLMG1v78L/KPwAl4z1QDgCtxfrN/j/lKvFWKbV+KaY7YCb+H6UmaHGc/vcZ3FO714/l2GYzjebdUC/oj7fL/HNS3+xpv3DPAOMEtEduE+/zMBVHUlcCvwGu7s5Ecgu5j9/APoICI7RGR6UQuo6ircD/MnuMTTEfg4xLIl/fvcDXwFfI7rz3gMqKaqe4CHgY+9WLrjks4/cRcmbMIl1KOuygrhFS/GKWEsa0pBVG1gK1P+iIjiOknX+x1LaVWGYygNEbkeGKmq/fyOBQr6c270mgVNFNgZiTEm0lJxZwy+8/qEfgVM9DuWyswSiTEmYrzmsPNxTV9+xzIQ19eyDde0Z6LEmraMMcaUiZ2RGGOMKZMqcUNikyZNNDk52e8wjDGmQlmyZMkPqppU0nJVIpEkJyeTmZnpdxjGGFOhiMg3JS9lTVvGGGPKyBKJMcaYMrFEYowxpkyqRB+JMaZqO3DgANnZ2ezbt8/vUMql2rVr06JFC+Li4kq1viUSY0yll52dTf369UlOTkakLAWpKx9VJTc3l+zsbFq1alWqbVjTljGm0tu3bx+JiYmWRIogIiQmJpbpbM0SiTGmSrAkElpZPxtLJMV59VV44QW/ozDGmHLNEklxpk2DZ5/1OwpjjCm13Nxc+vbtS7169bjtttuisg9LJMU59VTYsAEOHSp5WWOMKYdq167NH/7wB5544omo7cMSSXFSUmD/ftiyxe9IjDEVWFZWFu3bt+eGG27g9NNP5+qrr+aDDz6gV69etGnThs8++wyAvLw8hg4dSlpaGt27d2f58uUcPnyY5ORkduzYUbC9U089lW3btpGTk8Oll15Kt27d6NatGx9/fOxAlfHx8Zx11lnUrl07asdnl/8Wp2dP9zx9Otx5p7+xGGMi4847YenSyG4zPR2efrrYRdavX8/UqVOZOHEi3bp147XXXmPhwoW88847PPLII0yfPp2HHnqIzp07M336dD788EOuvfZali5dysUXX8xbb73F6NGjWbx4McnJyTRt2pSrrrqKu+66i7POOovNmzczcOBAVq9eHdljC4OdkRSnY0f3mDnT70iMMRVcq1at6NixI9WqVSM1NZX+/fsjInTs2JGsrCwAFi5cyDXXXANAv379yM3NZefOnYwYMYI33ngDgNdff50RI0YA8MEHH3DbbbeRnp7OkCFD+Omnn9i1a1fMj83OSEpyxhnw1ltw+DBUs7xrTIVXwplDtNSqVavgdbVq1QreV6tWjYMHDwLu5sDCRIQePXqwfv16cnJymD59Og888AAAhw8f5pNPPqFOnToxOILQ7JexJL17Q14efPWV35EYYyq53r178+qrrwIwb948mjRpQoMGDRARhg0bxpgxYzjttNNITEwEYMCAATz33HMF6y+NdJNdmOyMpCRnnOGely2DTp38jcUYU6mNGzeO0aNHk5aWRt26dZk8eXLBvBEjRtCtWzdefvnlgmnPPvsst956K2lpaRw8eJDevXvzQhH3viUnJ/PTTz+xf/9+pk+fzqxZs+jQoUPE4q4SY7ZnZGRoqQe2OnAAGjSAW26BP/85soEZY2Ji9erVnHbaaX6HUa4V9RmJyBJVzShpXWvaKklcHHTtCosX+x2JMcaUS5ZIwtGhA6xd63cUxhhTLlkiCUerVvDDD7B7t9+RGGNMuWOJJByBGv0bN/obhzHGlEOWSMLRtq17XrfO3ziMMaYcskQSjkAi8aH0gDHGlHeWSMJRr55r3lq2zO9IjDHmuMyePZuuXbvSsWNHunbtyocffhjxfdgNieFKS7MzEmNMhdOkSRP+85//0Lx5c1asWMHAgQP59ttvI7oPOyMJV2BsksOH/Y7EGFPB+FlGvnPnzjRv3hyA1NRU9u3bx88//xzR47MzknC1bQv79sHmzZCc7Hc0xpjSqsJl5KdNm0bnzp2PKiAZCZZIwnX66e551SpLJMaY4xYoIw8UW0Z+2rRpwLFl5MePH8/o0aOPKSO/atWqgn0EysjXr1//mP2vXLmS++67j1mzZkX82CyRhOvEE91zTo6/cRhjyqYKlpHPzs5m2LBhvPLKK6SkpETqkApEvY9ERKqLyJciMiPE/MtFZJWIrBSR14Km/1dEdhReT0ReFpFNIrLUe6RH+xgAaNLEPW/fHpPdGWOqnmiUkd+xYwcXXnghjz76KL169YpK3LHobL8DKLLRTkTaAGOBXqqaCgSPZ/s4cE2Ibd6jquneIzYF+OvXh2bNYPnymOzOGFP1jBs3jszMTNLS0rj//vuPKSM/ZcqUgmYtcGXkA8t36NChyBLyzz33HOvXr+cPf/gD6enppKensz3CfxBHtYy8iLQAJgMPA2NUdXCh+X8C1qrqSyHWPwe4O3g9EXkZmKGq/wo3jjKVkQ92ySWwYAFs22ajJRpTgVgZ+ZKV5zLyTwP3AqGumW0LtBWRj0XkUxE5P8ztPiwiy0XkKREp8vIDEblRRDJFJDMnUv0ageKN554bme0ZY0wlELVEIiKDge2quqSYxWoAbYBzgCuBl0SkUQmbHgu0B7oBCcB9RS2kqhNVNUNVM5KSko43/KJddZV7njsXqsCAYMYYE45onpH0AoaISBbwOtBPRKYUWiYbeFtVD6jqJmANLrGEpKrfqfMzMAk4I/Khh9C1K0ya5F7PmROz3RpjTHkWtUSiqmNVtYWqJgNXAB+q6shCi00H+gKISBNcU1extdpFpJn3LMBQYEWEQy/elVdCzZowe3ZMd2uMMeVVzO8jEZHxQKaqvgO8DwwQkVXAIdzVWLnecgtwTVj1RCQb+IWqvg+8KiJJgABLgZtjegC1arm6W0uKa7EzxpiqIyaJRFXnAfO81w8GTVdgjPcovM7ZIbbVLypBHo927aCImjbGGFMV2TWspXHyyZCdDYcO+R2JMcYU67PPPiu4f6RTp0689dZbEd+HlUgpjZYt4eBB+P57OOkkv6MxxpiQTj/9dDIzM6lRowbfffcdnTp14qKLLqJGjcj9/NsZSWl06OCe7S53Y0wY/CwjX7du3YKksW/fPtx1SpFlZySlERh616vYaYypQKpgGfnFixdz/fXX88033/DPf/4zomcjYImkdJKSQMQ1bRljTBj8LCN/5plnsnLlSlavXs2oUaMYNGgQtWvXjtixWSIpjRo14IQT4Lvv/I7EGHO8qmAZ+YDTTjuN+Ph4VqxYQUZGiSW0wmZ9JKXVrBlEeNxjY0zVFo0y8ps2bSpIVN988w1r1qwhOcKD89kZSWm1aWM3JRpjImrcuHGMHj2atLQ06tate0wZ+W7duvHyyy8XTHv22We59dZbSUtL4+DBg/Tu3fuYUvILFy7kj3/8I3FxcVSrVo2//e1vNAmMrxQhUS0jX15ErIx8sAkT4He/g507oUGDyG7bGBNRVka+ZOW5jHzl1bmze162zN84jDHGZ5ZISiuQSIYN8zcOY4zxmSWS0mrWzD3n5rrmLWOMqaIskZSWCPz61+510HXcxhhT1VgiKYs77nDPn37qbxzGGOMjSyRl0bIlnHYavPuu35EYY4xvLJGU1cUXw7x51k9ijCnXNm/eTL169XjiiScivm1LJGXVvbsbl2TtWr8jMcaYkO666y4GDRoUlW1bIimr1q3d88Zih5o3xlRhfpaRB5g+fTqtW7cmNTU1KsdnJVLKqlUr9/z11/7GYYwJTxUrI5+fn89jjz3G7Nmzo9KsBXZGUnb16kGXLjBnjt+RGGPKsUAZ+WrVqhVbRv6aa64Bji0j/8YbbwAcU0b+tttuIz09nSFDhhSUkQ/20EMPcdddd1GvXr2oHZudkURCp04wa5bfURhjwlHFysgvXryYf/3rX9x7773s2LGDatWqUbt2bW677baIHZudkUTCCSfA9u1QBQpgGmOiJxpl5BcsWEBWVhZZWVnceeed/OY3v4loEgE7I4mM5s3hwAH44Qc3eqIxxpRCNMrIx4KVkY+Et9+GoUNh9mw499zo7ccYUypWRr5kVkbeb+3auefhw/2NwxhjfGCJJBLatYMmTWDHDvj5Z7+jMcaYmLJEEgki8Le/udcrV/obizHGxJglkkjp1s09h7iz1BhjKitLJJGSnAxt28KkSeBdE26MMVWBJZJI+t3v4MsvYe5cvyMxxpiYsUQSSeed556XLfM3DmOM8WRlZVGnTh3S09NJT0/n5ptvjvg+7IbESGra1JVLmTYN7r7b72iMMQaAlJSUIu96jxQ7I4m0IUPgs8/gu+/8jsQYU074XUY+2uyMJNKGD4fHH4feveG//4WUFL8jMsYEq2Jl5AE2bdpE586dadCgARMmTODss8+O6OFbIom0jh3hww/hoovgyivd2YkxpsoLlJEHii0jP23aNODYMvLjx49n9OjRx5SRX7VqVcE+AmXk69evXzCtWbNmbN68mcTERJYsWcLQoUNZuXIlDRo0iNixWSKJhh494Je/hCefdHe6B5WPNsb4rIqVka9Vq1bBvrp27UpKSgpr164lI6PEElphi3ofiYhUF5EvRWRGiPmXi8gqEVkpIq8FTf+viOwovJ6ItBKRxSKyTkTeEJGa0T6GUunZ01UEtjMSY0yYolFGPicnh0OHDgGwceNG1q1bR+vAEOEREovO9juAYxvtABFpA4wFeqlqKnBn0OzHgWuKWO0x4ClVbQP8CPwisuFGSJcu7vmrr/yNwxhTYYwbN47MzEzS0tK4//77jykjP2XKlIJmLXBl5APLd+jQocgS8vPnzyctLY1OnTpx2WWX8cILL5CQkBDRuKNaRl5EWgCTgYeBMao6uND8PwFrVfWlEOufA9wdWE9EBMgBTlTVgyLSAxinqgOLiyPqZeSLouoKOfbrB1OnxnbfxpijWBn5kpXnMvJPA/cCh0PMbwu0FZGPReRTETm/hO0lAjtUNVCDJBs4qagFReRGEckUkcycnJzSxF42IjBypBurZP/+2O/fGGNiJGqJREQGA9tVdUkxi9UA2gDnAFcCL4lIo+I2W8S0Ik+pVHWiqmaoakaSX6MW9u7t+knmzfNn/8YYEwPRPCPpBQwRkSzgdaCfiEwptEw28LaqHlDVTcAaXGIJ5QegkYgErjZrAWyNbNgRNHiwG889qDPMGGMqm6glElUdq6otVDUZuAL4UFVHFlpsOtAXQESa4Jq6NhazTQXmApd5k0YBb0c49MipVQsuvhgWLfI7EmOMiZqYl0gRkfEiMsR7+z6QKyKrcAniHlXN9ZZbAEwF+otItogEOtTvA8aIyHpcn8k/YnsEx6ltW8jNdaMnGmNMJRSTGxJVdR4wz3v9YNB0BcZ4j8LrFHkPv6puBM6IRpxR0cZrqVu92t2oaIwxlYwVbYy2wMiJM4q8H9MYY6Ju+fLl9OjRg9TUVDp27Mi+ffsiun1LJNHWvLmrv/XMMzZyojEm5g4ePMjIkSN54YUXWLlyJfPmzSMuLi6i+7BEEgtjxkB+Pqxd63ckxhgf+FlGftasWQV3tgMkJiZSvXr1iB6fFW2MhXPOgRo14Kmn4O9/9zsaY6q2KlZGfu3atYgIAwcOJCcnhyuuuIJ77703oodviSQWkpPhttvcF+3++22MEmOqIL/KyB88eJCFCxfy+eefU7duXfr370/Xrl3p379/xI7NEkmsXHutSyRLllgiMcZPVayMfIsWLejTpw9NmjQB4IILLuCLL76IaCIJq49EnJEi8qD3/mQRqTiX4JYH7du7+ltBfz0YY0ywaJSRHzhwIMuXL2fPnj0cPHiQjz76iA4dOkQ07nDPSP6GK7zYDxgP7AKmAd0iGk1lVqcOtG5ticQYE9K4ceMYPXo0aWlp1K1b95gy8t26dePll18umPbss89y6623kpaWxsGDB+ndu/cxpeQbN27MmDFj6NatGyLCBRdcwIUXXhjRuMMqIy8iX6hqFxH5UlU7e9OWqWqniEYTJb6UkS/KRRdBVpaNUWJMjFkZ+ZLFooz8ARGpjldpV0SSCF0a3oTSvr27BNgbrcwYYyqDcBPJs8BbwAki8jCwEHgkalFVVm3auLFJsrP9jsQYYyImrD4SVX1VRJYA/XFjggxV1SKHzzXFaN/ePS9ZAqec4m8sxhgTIcWekYhIQuABbAf+D3gN2OZNM8ejRw9o2hSmFB6WpRgzZ8KKFdGLyRhjyqikpq0lQKb3nAOsBdZ5r4sb+dAUJS4OLr8c3nsP9u4teflVq+CCC1ytrtdfd+PAG2NMOVNsIlHVVqraGjduyEWq2kRVE4HBwL9jEWClc+658PPP8MUXJS87adKR11deCW++Gb24jDGmlMLtbO+mqu8F3qjqTKBPdEKq5NLS3PPKlSUv+8EHLvEERlh84AH46afoxWaMiZmnn36aPXv2+B1GRISbSH4QkQdEJFlEThGR3wK50Qys0jr5ZKhXr+R+j/x8WL7c9av06AH//jesXw9//nNs4jTGRFVVTCRXAkm4S4CnAyd408zxqlYNUlNLvinxyy/h8GE4w6tEM2wYDBrkxjWxYXuNqVDy8/O58MIL6dSpE6effjq///3v2bp1K3379qVv374A3HLLLWRkZJCamspDDz0EwJw5cxg2bFjBdmbPns0ll1ziyzEUJ9zLf/OAO6IcS9XRsSO89ZbrPBcpeplAogk0hQHcfTf07w+ffgrnnx/9OI0xEfHf//6X5s2b8+677wKwc+dOJk2axNy5cwuKKT788MMkJCRw6NAh+vfvz/Lly+nXrx+33norOTk5JCUlMWnSJEaPHu3noRQp3KKNc0Xkw8KPaAdXaaWlQW4ubN0aepmVK6F+fWjZ8sg0b2Aabr8d3n47ujEaU8Xl5e/nxY82kJe/v8zb6tixIx988AH33XcfCxYsoGHDhscs8+abb9KlSxc6d+7MypUrWbVqFSLCNddcw5QpU9ixYweffPIJgwYNKnM8kRZu0ca7g17XBi4FbNzY0urc2T2/+y7ceGORi+zq1YdV++vQZs8BEuJruomJidCsmesrGTrU9bOkpsYoaGOqlqmZW3h05tcA3NSnbEM/tG3bliVLlvDee+8xduxYBgwYcNT8TZs28cQTT/D555/TuHFjrrvuuoJx1UePHs1FF11E7dq1GT58ODVqlL/RP8I6I1HVJUGPj1V1DHBmlGOrvHr2hO7dYcIEOHCgyEVeatKJEQnnMHlR1tEz5s6FJ55wTWKvvRb9WI2pooZntGTsoPYMz2hZ8sIl2Lp1K3Xr1mXkyJHcfffdfPHFF9SvX59du3YBbkCq+Ph4GjZsyLZt25g5c2bBus2bN6d58+ZMmDCB6667rsyxRENYqa3QXezVgK7AiVGJqCqoVg3uvRcuuQQWLIB+/YpYSAs9e9q1c4+ZM+E//4GHH452tMZUSQnxNct8JhLw1Vdfcc8991CtWjXi4uJ4/vnnC5qpmjVrxty5c+ncuTOpqam0bt2aXr16HbX+1VdfTU5OTsTHEYmUcM+RluB+0QTXpLUJ+EW0gqoSAl+UpUuPSSSBNtk7+rdhVM/kotfv1g0efxw2b3aXFBtjyq2BAwcycODAo6ZlZGRw++23F7wPHmeksIULF/LLX/4yWuGVWbiX/56mqq29O93bqOoA4PNoBlbpnXACtGgBn31WMCnQuTd5URbPzFlP3ZrVj/SPFHbDDa4c/RNPxChgY4wfunbtyvLlyxk5cqTfoYQUbiJZVMS0TyIZSJXUqxcsXFhQQ2vyoiwenfk1P+bvp2+7JM7t0DT0uikprq9l+fIYBWuM8cOSJUuYP3/+UWO+lzfFNm2JyInASUAdEemMa9oCaADUjXJsld6uIZewkfq03JEPNWuSmZUHwPx1OWTl7qF7622k9KkXegPt2rn7UfbudUP5GmNCUlUk1H1bVVw4I+UWp6Q+koHAdUALILg2xy7gN2Xas+G1k7ryaHI8faevIq1FIz7ekEtCfBxZuXvo2y6p5KtFrrnGFXacOhWuvTY2QRtTAdWuXZvc3FwSExMtmRSiquTm5lK7du1SbyPcMdsvVdVppd6Lz8rNmO2F5OXv5/bXvuDjDbnceHYr1m3fzdw1OfRtl8STl6eH7h8JUHWjLoKrJtygQfSDNqYCOnDgANnZ2QX3Zpij1a5dmxYtWhAXF3fU9HDHbC+paWukqk4BkkVkTOH5qmoVBMsgIb4mGckJfLwhl5krvqdv+xMKrtQqMYmAu5fkmWdg8GD45z/h1lujH7QxFVBcXBytWrXyO4xKq6TO9njvuR5Qv9CjmMZ7E65RPZNJTqzLlh/38son3xR/pVZRLrzwyH0lxhjjg2LPSFT1Re/lB6r6cfA8EelVxCrmOCXE12RA6olMnL+RM1sllO4u2tNOg/nziy8CaYwxURLu5b9/CXOaKYWb+6QwdlB7nh/Z9fjORgJSUyEvz85KjDG+KKmPpAfQE0gq1EfSAKgezcCqkjKXYhg7Fh55xI2keMEFkQvMGGPCUNLlvzVxfSE1cP0iAT8Bl0UrKHOc4uOha1fXvGWMMTFWUh/JR8BHIvKyqn4To5hMafTvD08+Cbt3u6F8jTEmRsLtI9kjIo+LyHvHO7CViFQXkS9FZEaI+ZeLyCoRWSkirwVNHyUi67zHqKDp80RkjYgs9R4nhHkMldt558HBgzB7tt+RGGOqmHATyavA10Ar4PdAFuEXbbwDWF3UDBFpA4wFeqlqKnCnNz0BeAg35skZwEMi0jho1atVNd17bA8zjsqtTx9o3NiVljfGmBgKN5Ekquo/gAOq+pGqXg90L2klEWkBXAi8FGKRXwJ/VdUfAYKSwkDPt8zXAAAclUlEQVRgtqrmefNmAzZIeXFq1IABA1wiCTFYljHGREO4iSTwy/SdiFzoFXBsEcZ6TwP3AodDzG8LtBWRj0XkUxEJJIuTgC1By2V70wImec1avxMrnHPEFVfADz/YyInGmJgKN5FMEJGGwK9x47e/hNcMFYqIDAa2q+qSYharAbQBzgGuBF4SkUYcqTIcLFAU7GpV7Qic7T2uCbH/G0UkU0Qyc3Jyigu18rjwQujSBf74R78jMcZUIeGO2T5DVXeq6gpV7auqXYGSbnzoBQwRkSzgdaCfiEwptEw28LaqHlDVTcAaXGLJBoJv8W4BbPVi+dZ73gW8hutDKSrmiaqaoaoZSUlJ4RxmxRcXB4MGwbp1ruPdGGNiINwzkqIcU8QxmKqOVdUWqpoMXAF8qKqFh/iaDvQFEJEmuKaujcD7wAARaex1sg8A3heRGt5yiEgcMBhYUYZjqHySk93IiRs2+B2JMaaKKEsiKVXfhIiMF5Eh3tv3gVwRWQXMBe5R1VxVzQP+gLsy7HNgvDetFi6hLAeWAt8Cfy/DMVQ+GV7F5/btYYXlWGNM9IU1HkmRK4psVtWTIxxPVJTX8Uii5u23YehQuPNOeOopv6MxxlRQkRqPZBdHOrmPmgXY2K7l1cUXw+mnw8aNfkdijKkCSiqRUr+4+aYca9sWVq70OwpjTBVQlj4SU5517Oiu3srP9zsSY0wlZ4mksurSBQ4fhuXL/Y7EGFPJWSKprDp1cs/LlvkbhzGm0rNEUlmdfDI0amSJxBgTdSUNbGUqKhF2//Fx5uyvz9n5+0s3hK8xxoTBzkgqsb+3Pps7vq3H//zfF+Tl7/c7HGNMJWWJpFJzxQcWrs/lhXnrfY7F+O5wqCLcxpSNJZJKbFTPZJIT6wKwLHsnT81ew1Oz1x57dnKcBR7z8vfz4kcbjt7O99/DQw/Zj1V59PPPcP31UL06dO8OEyfCrFnw9dd+R2YqCesjqcQS4mvyj+u68eD0Few9cIhn5rizkuXZO3hyWCoJM99xPzKjR4MI7NoF8fGASxaTF20ChFE9k4/qY5m8KItn5qxjwbofePbKziTE12T7m29xwvjxPLdhPz9ddS0390mptP0yefn7mZq5heEZLSvGMV5+Obzzjnu9bh3cdNOReV9+Cenp/sRlKg07I6nkUpLqkdq8AV9s3gFA84a1qfbuDA62aQtXXeWSCIAq9/x1NhtydgOBZLGeZ+asY2rmFvLy9/PU7DU88u5qFqxz47ssXP8Dv35zKRtydnPl/vY8MOBXPNc0g4nzNzI1c0uR8ZRngTOtL775kdGTPiv4LILnBRLsozO/9hJtObdx45EkMnmyG/hs8WJXRgegc2fYu9e/+EylYGckVcCy7J0Fr09sWIsahw6xqXZjfuzTiXYfzYSbbuL21hfwn7zqbJ2+gr9c1YXMrDwAzjo1keEZLZmauaXgjCagZeM6zF2Tw869y9jwwx62db+IfT8fosvJjdiz/xB5FexqscmLNvHMnPUkJ9YlK3cPsIpJo91wN1Mzt/DozEBTkBR6LsfuuefI66FD3ZnnGWfA9Onwu9/BhAku0YwY4V+MpsKzRFIFdGrZiMWb8jizVQKdWjZiYruevN+uJ3f0P5V2894DoNl7q2H+Rjbn7eGFjzbw8YZc+rZL4snL00mIr8nwjJbs2X+QvftdH0idmtXZe+AQE+dvZN32XQC0bVqfgaknsmf/IZ6Zs466NatzU59jxz8r3DRUXDNatAXv+8d8N6L0Ga0SaNUkngcGdyhYbnhGy6Oe69asXvC63Dp0yDVdATz+ODRocPT8ceNg2jS46y53A2v79jEP0VQOlkiqgJv7pJDoJQOAOnHVCPxoBy8zZ/U2NuTks/LbnYwd1J5zOzQ96gf/rvPaHbXdvPz9rNu2i7lrXFNXximNualPChtydrM8ewfndmh61LKBH2yAZ+asA+CmPikFzWhA2MmnLIK3FXymdWarxgA0a1ibP13W6ah1EuJrHhVX4Rj9TIYhrVwJmzbByy/DqFHHzq9eHaZOhX794IYbYOHCmIdoKgdLJFVA4R/BwgkhsMzEazOYMGMVDwzuQEpSPV78aENBc05RP+4J8TV58vL0o35AAT5YtY25a3LYnJfJxGszSEmqd9QP9rXdT6FXSiK5+fu9K7/cSAVntmocskks0Oy0Z//BIuM/HoFmqk835vLA4A4sWJfDwvW5gHBH/zZHJdhwBeKD0Mkw5hYscM9duoReJjXVdb4//DDs3AkNG8YmNlOpWCIxBVKS6hX0CcCxzTlFKepMZXhGS6Z8+g0bcvIZO205j16axpzV22jesDZbd+5j/rocsnL38PGGXBLjazKqZyvq1qxRQpNYZPolNuTsZv7aHM5slcDcNTmktfiWDs0akP3jXhZvyqNf+xNKdTax94Br8juzVUKxn9eGnN1HJeuyCJwFBZobwTU5DklvzpyvtnL98y9Qo1s36NChmK3gKkUfPgybN7vXxhwnSyQmpMJnMqUhIkyYsYrPsn4EICUpng05+fRKSSQjOaGgqeqmPink5e8P2fcwqmdyRPolJsxYxccbcumVksjYQe3Zs/8QExe4q6/6tksq9fZdcyF0b51wTCIKbkqbMGOV1xS4iicvTy+yua64ZrzgecFNgsGWZ+9g7pocaj38D0adebJrwirOCSe45++/t0RiSsUSiYm4qZlb2PLjXpIT69KpZSPOTz2R/QcPk3pSQ0Z0a8kHq7YV+SNZXOKKRFIDvA70I2cER5rWyta3ETirKioRBV/xFbz/4OmBRDo1c0vBmRngJQzXdDgkvXlBItqz/yCZXnJuWKcGO/cepGXjOlzSpQVD0pvTvfU2LspoCeEcz+mnQ82a8O67cN55pTp+U7WVesz2iqTKjdnus8I/iGMHtS8ffQY+CXWGUXh6oE/qjv6nFiSlwM2fAKck1OWbvD2cdWoTup7SmGfmrKNvuyRu79eGv3y4rmzNZQMHwnff2fg15igRGbPdmNIIp6mqKgl1NlV4enCf1JGE4/7QO3JvC3Q9pdFRTX0J8TWP6tsqlbPOciVu8vIgIaFs2zJVjp2RGFOOBc5azu3QlHeWfkvULi9euBDOPhv+9S+49NLIbttUWOGekVgiMcbAgQOQlOTucH/xRb+jMeVEuInEam0ZYyAuztXdshE1TSlYIjHGOKef7u6G37fv6Om5uTBvng0RYEKyRGKMcS65BHbvhn//+8i0nBw46STo2xeee86/2Ey5ZonEGOP06QONG8P77x+Z9uqrbswagDvugOef9yc2U65ZIjHGONWqwbBh8PrrsGMHqMIrr0BamrvHpEcP+PWvrYnLHMMSiTHmiMsvh/37YfBgmD/flaG/5RY48UQ3ENreva65a8sWV+xxS8UbwMxEniUSY8wRPXq4548/hnPOca/79nXPJ5/snj/5xM2bONH1qwSavkyVZYnEGHNEgwZu5MTAWcgpp8Cpp7p5bdq45xtucEP43nYbZGayYOxjXs2yI4KHJjaVnyUSY8zRxo+Hv/0NNmyAtWuPVA8+7TTXIZ+bCykp8Je/MOPJV7g2rjPDX1h01Bj3gYKUUzOt6asqsERijCla3bquKnCwMWPc82WXAXDaNcNoXK8WG3LymTBjVcFiwzNaMnZQ+ypfZ62qsKKNxpjwDRniRl7s2hVwo2Hm5R+gTlw1RvVILlgsUmX/TcVgZyTGmONz1llQpw7gzjxSkuLZe+Awkz/JKljE+kiqFkskxphSS4ivycRrM+jbLskbtMuxPpKqxZq2jDFlkpJU75jxUILHVjGVX9TPSESkuoh8KSIzQsy/XERWichKEXktaPooEVnnPUYFTe8qIl+JyHoReVZEJNrHYIwxJrRYNG3dAawuaoaItAHGAr1UNRW405ueADwEnAmcATwkIo291Z4HbgTaeI/zoxq9Mea4Tc3cwp/eXcmf/zzN+kmqgKgmEhFpAVwIvBRikV8Cf1XVHwFUdbs3fSAwW1XzvHmzgfNFpBnQQFU/UTci1yvA0GgegzHm+A3PaMmkT17irkdvZvpHq0pewVRo0T4jeRq4FwhV5a0t0FZEPhaRT0UkcHZxEhDcS5ftTTvJe114+jFE5EYRyRSRzJycnLIcgzHmOCXE16TThPtI3PsTI1bP8zscE2VRSyQiMhjYrqpLilmsBq556hzgSuAlEWkEFNXvocVMP3ai6kRVzVDVjKSkpOOK3RhTdg379IJ27YifNdPvUEyURfOMpBcwRESygNeBfiIypdAy2cDbqnpAVTcBa3CJJRsIvtyjBbDVm96iiOnGmPJo8GA3uuLu3SUuaiquqCUSVR2rqi1UNRm4AvhQVUcWWmw60BdARJrgmro2Au8DA0SksdfJPgB4X1W/A3aJSHfvaq1rgbejdQzGmDI65xxXln75cr8jMVEU8xsSRWS8iAzx3r4P5IrIKmAucI+q5qpqHvAH4HPvMd6bBnALrvN+PbABsPNmY8qr9u3d85o1/sZhokrcxU+VW0ZGhmZmZvodhjFVz8GDrjT9qFE2TG8FJCJLVDWjpOWsRIoxJnpq1HDD906ZAtnZJS9vKiRLJMaY6Bo/3nW2Tyl8rY2pLCyRGGOiKyXFjbS4bJnfkZgosURijIm+tm1h/Xq/ozBRYonEGBN97du7K7eqwMU9VZElEmNM9HXoALt2weLFfkdiosASiTEm+gYPdmPAP/mk35GYKLBEYoyJvhYt3F3u1k9SKVkiMcbERuvWsGEDHA5VDNxUVJZIjDGxccYZrp/EqkxUOpZIjDGx0b+/e/78c3/jMBFnicQYExvNmrkOd+snqXQskRhjYkPE9ZOsXOl3JCbCLJEYY2Ln4oth9mx47TW/IzERZInEGBM7Dz4IbdrAxIl+R2IiyBKJMSZ2ataEoUNh0SJYsMDvaEyEWCIxxsTWmDGuIvD558Mnn/gdjYkASyTGmNg68USYN89dwfXXv/odjYkASyTGmNhr2hS6drUruCoJSyTGGH907AirV7tx3U2BvPz9PPLuKka+tJgNObv9DicslkiMMf7o0gV+/hleeSXkInn5+3nxow3k5e+PYWD+2JCzm9GTPuOP761m4oJNLFz/Aw9OX1Ehjr+G3wEYY6qogQPd8yOPwMiR7oouXPJ4Yd56Vn23iw7NGzBx/kYAbuqTEpUw8vL3MzVzC8MzWpIQXzMq+wjHhBmrmLsmh9px7u/7xnXjSD2pIY/O/BqI3vFHgiUSY4w/EhJg+nR3OfAHH8AFFwAweVEWExdsAkBVGTuoPcMzWkYtjKmZW8rFj/UDgzuwOS+TDTn5pCTFM/HaDBrXrUlifM2oHn8kWCIxxvinTx/3/NVXcMEFbMjZzdtLvwUgObEu44eeTkpSvaiGEPiR9vvHOiWpHlNv7nnM2VF5PhMJsD4SY4x/GjVyg14tXQrAQ2+vJCt3D8mJdfn3r3pFPYkAJMTX5NwOTfn1m0t979xOiK/JTX1Swm5iKy99SJZIjDH+6t0b5s0jL38/+T8fcJPaNIlpf8WD01cwd00OD05fEbN9RkKgWW5q5hbAv8RiicQY468uXeD773n5P0v4cstOABrH14ppCKknNTzquaIYntHyqD6kyYs28ejMr5m8aFNMk4r1kRhj/JWaCsDPn30ONVNITqzLqJ7JMQ3h5j4pFaJTu7BAU9gRUvAcy4sILJEYY/zVsyeIcFPcdta2684DgzvE/DLcY3+QK6ZRPZOpW7P6UQkxFslRVDXqO/FbRkaGZto40caUX6mpcMIJMHeu35GYICKyRFUzSlrO+kiMMf477zz49FOoAn/YVkaWSIwx/ktJgX374Icf/I7ElIIlEmOM/9LTYcQI2F++a0qZollnuzHGf2ef7R6mQrIzEmOMMWViicQYY0yZRD2RiEh1EflSRGYUMe86EckRkaXe44ageY+JyArvMSJo+ssisilonfRoH4MxxpjQYtFHcgewGmgQYv4bqnpb8AQRuRDoAqQDtYCPRGSmqv7kLXKPqv4rWgEbY4wJX1TPSESkBXAh8NJxrtoB+EhVD6pqPrAMOD/S8RljjCm7aDdtPQ3cCxwuZplLRWS5iPxLRAL38i8DBolIXRFpAvQFgu/zf9hb5ykRKbK6m4jcKCKZIpKZk5MTiWMxxhhThKglEhEZDGxX1SXFLPYfIFlV04APgMkAqjoLeA9YBPwf8Alw0FtnLNAe6AYkAPcVtWFVnaiqGaqakZSUFIEjMsYYU5RonpH0AoaISBbwOtBPRKYEL6Cquar6s/f270DXoHkPq2q6qp6HK2m5zpv+nTo/A5OAM6J4DMYYY0oQk6KNInIOcLeqDi40vZmqfue9Hgbcp6rdRaQ60EhVc0UkDXgNSFfVg4F1RESAp4B9qnp/CfvPAb4pZfhNgIpYt8Hiji2LO7Ys7tg4RVVLbNKJ+Z3tIjIeyFTVd4D/EZEhuGarPOA6b7E4YIHLFfwEjFTVQNPWqyKShDtLWQrcXNI+w/kgiok3M5zql+WNxR1bFndsWdzlS0wSiarOA+Z5rx8Mmj4W1+dRePl9uCu3itpWv6gEaYwxplTsznZjjDFlYomkZBP9DqCULO7Ysrhjy+IuR6rECInGGGOix85IjDHGlIklEmOMMWViiaQYInK+iKwRkfUiUuy9KjGIpaWIzBWR1SKyUkTu8KaPE5Fvg6ohXxC0zlgv9jUiMjBoekyPS0SyROQrL75Mb1qCiMwWkXXec2NvuojIs15sy0WkS9B2RnnLrxORUVGOuV3QZ7pURH4SkTvL6+ctIv8rIttFZEXQtIh9xiLS1fs3XO+tK1GK+XER+dqL6y0RaeRNTxaRvUGf+wslxRbq+CMhROwR+26ISCsRWezF/oaI1IxU7FGhqvYo4gFUBzYArYGauPpfHXyMpxnQxXtdH1iLu0R6HO5mz8LLd/BirgW08o6luh/HBWQBTQpN+xNwv/f6fuAx7/UFwEzcfULdgcXe9ARgo/fc2HvdOIbfhe+BU8rr5w30xlXMXhGNzxj4DOjhrTMTGBSlmAcANbzXjwXFnBy8XKHtFBlbqOOP4ucdse8G8CZwhff6BeCWWHzXS/uwM5LQzgDWq+pGVd2PK/NysV/BqCsN84X3eheuNP9JxaxyMfC6qv6sqpuA9bhjKi/HdTFebTXveWjQ9FfU+RRoJCLNgIHAbFXNU9UfgdnEriJ0f2CDqhZXHcHXz1tV5+Nu6i0cU5k/Y29eA1X9RN0v2ytB24pozKo6S4/cfPwp0KK4bZQQW6jjL7MQn3cox/Xd8M6o+gGBoTIiGns0WCIJ7SRgS9D7bIr/4Y4ZEUkGOgOLvUm3eU0B/xt0+h4qfj+OS4FZIrJERG70pjVVrzyO93yCN708xR1wBa54aEB5/7wDIvUZn+S9Ljw92q7HnWEEtBI3SN5HIhIY4L242EIdfzRF4ruRCOwISqjl5rcnFEskoRXVBuz7tdIiUg+YBtypbqCv54EU3CBg3wFPBhYtYnUtZno09VLVLsAg4FYR6V3MsuUpbry26SHAVG9SRfi8S3K8scb8GETkt7jSSa96k74DTlbVzsAY4DURaeBHbMWI1HejPB1TWCyRhJbN0WOgtAC2+hQLACISh0sir6rqvwFUdZuqHlLVw7gKyoFqyKHij/lxqepW73k78JYX4zavWSLQPLG9vMXtGQR8oarboGJ83kEi9Rlnc3QTU1SPwevkHwxc7TVX4TUL5Xqvl+D6FtqWEFuo44+KCH43fsA1N9YoNL3cskQS2udAG+/qiZq45o13/ArGazf9B7BaVf8cNL1Z0GLDgMBVJO8AV4hILRFpBbTBdUrG9LhEJF5E6gde4zpTV3j7DFwVNAp4Oyjua70ri7oDO71mifeBASLS2GsyGOBNi7YrCWrWKu+fdyER+Yy9ebtEpLv3Pbw2aFsRJSLn48YYGqKqe4KmJ4mrCo6ItMZ9vhtLiC3U8UdFpL4bXvKcC1wWq9jLzO/e/vL8wF3dshb3189vfY7lLNzp7XJc1eOlXnz/BL7ypr8DNAta57de7GsIusomlseFuyJlmfdYGdgfrh14Dm6cmTlAgjddgL96sX0FZARt63pcR+V6YHQMPvO6QC7QMGhaufy8ccnuO+AA7i/dX0TyMwYycD+MG4Dn8KpiRCHm9bh+g8B3/AVv2Uu9788y4AvgopJiC3X8Ufy8I/bd8P7ffOZ9HlOBWtH+vpflYSVSjDHGlIk1bRljjCkTSyTGGGPKxBKJMcaYMrFEYowxpkwskRhjjCkTSySmwhGRQ1511ZUiskxExohIsd9lcdVjV3iv04Mrs4a5z2QRuSrofYaIPFu6IyhxX3eKyLXR2HZpiKve3KSY+a+LSJtYxmTKF0skpiLaq6rpqpoKnIe7Fv+h41g/3VvneCQDBYlEVTNV9X+Ocxsl8u5mvh54LdLbjqLngXv9DsL4xxKJqdDUlV25EVcsT0SkurgxLT73iufdFLy8dwfxeGCEd1YzQkTOEJFFXkHARSLSrohd/RE421vnLhE5R0RmeNscJyKTRWSW99f7JSLyJ3FjZPzXK20TGDfjI3HFK98vdCd0QD9cSZaD3jr/IyKrvGN53ZsWL64o4OdezBd706uLyBPefpeLyO3e9P7ecl9569XypmeJyO9F5AtvXntveqJ3LF+KyIt4tZ+8/b7rnQWuEJERXswLgHPlSEkPU9X4fUekPexxvA9gdxHTfgSa4pLKA960WkAmbgyIZLyxI4DrgOeC1m3AkTEwzgWmFbH9c4AZRb3HjUOxEIgDOgF7ODImxlu4EuBxwCIgyZs+AvjfIvbze+D2oPdb8e5qBhp5z48AIwPTcHdGxwO34GqxBY4lAaiNu1O8rTftFVzBT3DjxNzuvf4V8JL3+lngQe/1hbiKCk1wd5f/PSi24Dv+ZwNd/f5u2MOfh52RmMoiUDF1AK6O1FJcmf1EXG2j4jQEpnp9KE8BqaXY/0xVPYArkVEd+K83/StcEmsHnA7M9mJ7gKLH2mgG5AS9Xw68KiIjcdVwwR3j/d525uGSxcm4JPiCemczqprn7XeTqq711p2MG5Qp4N/e8xIvTrz5U7xtvItL0oFjOVdEHhORs1V1Z9B2tgPNQ3w2ppKzU1FT4XlF/A7hfswE91f2+4WWSS5mE38A5qrqMG+5eaUI42cAVT0sIgdUNVB76DDu/5kAK1W1Rwnb2YtLDAEX4n7YhwC/E5FUb1uXquqa4BW9ooWFax6VNCTuz97zIY7+PTimdpKqrhWRrrj+pUdFZJaqjvdm1/ZiN1WQnZGYCk1EknBDkT7n/Xi/D9wS1C/RVlzV4WC7cMMVBzQEvvVeXxdiV4XXOV5rgCQR6eHFFeclhcJWA6d6y1QDWqrqXFxndiOgHu4Yb/cSByLS2Vt3FnBzoK9CRBKAr4FkETnVW+Ya4KMSYp0PXO1tYxBu2F1EpDmwR1WnAE/ghpoNaIsrqmiqIEskpiKqE7j8F/gA9wP6e2/eS8Aq4AuvqepFjj3zngt0CHS248b2flREPsY1SxVlOXDQ62i+63gDVjeU6mXAYyKyDFfZtmcRi87kSNNTdWCKiHwFfAk8pao7cGdQccBy7xj/EHTsm73py4CrVHUfMBrXdPcV7gzphRLC/T3QW0S+wDWjbfamdwQ+85rUfgtMABCRprgr6b4L+wMxlYpV/zWmnBGRt4B7VXWd37GEw0usP6nqP/yOxfjDzkiMKX/ux3W6VxQ7cJ34poqyMxJjjDFlYmckxhhjysQSiTHGmDKxRGKMMaZMLJEYY4wpE0skxhhjyuT/AR79+6yWBJFGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "color=['r','y','b','g','p','v']\n",
    "for i in test['stays'].unique():\n",
    "    if i!=0:\n",
    "        plt.plot(test[test['stays']==i]['delta_seconds'],test[test['stays']==i]['lat'],'r');\n",
    "    else:\n",
    "        plt.scatter(test[test['stays']==i]['delta_seconds'],test[test['stays']==i]['lat'],1);\n",
    "plt.xlabel ('Delta time (seconds)');\n",
    "plt.ylabel ('Latitude');\n",
    "plt.legend(['move 1','move 2','move 3','move 4','move 5','move 6','stay']);\n",
    "plt.title('Example of a parsed trajectory');\n",
    "plt.savefig('parsed_example.eps', format='eps', dpi=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter the data that has stays greater than  so that we consider only moving people\n",
    "data=data[data['stays']>0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create new activity ids for the \"new\" trips\n",
    "data.loc[:,'new_activity_id']=data['activity_id'].astype('str')+'-'+data['stays'].astype('int').astype('str')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trips Description"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate instantaneous speed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to calculte distance\n",
    "def calculate_dist (lat_1,lat_2,lon_1,lon_2):\n",
    "    if (math.isnan(lat_1)|math.isnan(lat_2)|math.isnan(lon_1)|math.isnan(lon_2)):\n",
    "         return float('nan')         \n",
    "    elif ((lat_1==lat_2)&(lon_1==lon_2)):\n",
    "        dist=0.0\n",
    "    else:\n",
    "        R_aver = 6374;\n",
    "        lat_1 = math.radians(lat_1)\n",
    "        lon_1 = math.radians(lon_1)\n",
    "        lat_2 = math.radians(lat_2)\n",
    "        lon_2 = math.radians(lon_2)\n",
    "        aux=math.cos(lat_1)*math.cos(lat_2)*math.cos(lon_1-lon_2) + math.sin(lat_1)*math.sin(lat_2)\n",
    "        if aux>1:\n",
    "            aux=1\n",
    "        dist = R_aver * math.acos(aux);\n",
    "            \n",
    "    return dist*1000  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.reset_index(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.rename(columns={'lat':'y1','lon':'x1'},inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the following point to each record to avoid using loops\n",
    "data.loc[1:,'x2']=data['x1'].tolist()[0:-1]\n",
    "data.loc[1:,'y2']=data['y1'].tolist()[0:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate instantaneous distance\n",
    "data.loc[:,'distance_inst']=data.apply(lambda row: calculate_dist(row['y1'],row['y2'],\\\n",
    "                                row['x1'],row['x2']),axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the time\n",
    "data.loc[:,'date_1']=data['date'].apply(lambda x: datetime.strptime(x,'%Y-%m-%dT%H:%M:%S'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the time for the destination\n",
    "data.loc[1:,'date_2']=data['date'].tolist()[0:-1]\n",
    "data.loc[0,'date_2']=data['date'][0]\n",
    "data.loc[:,'date_2']=data['date_2'].apply(lambda x: datetime.strptime(x,'%Y-%m-%dT%H:%M:%S'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Delta time (instantaneous)\n",
    "data.loc[:,'delta_sec']=(data['date_1']-data['date_2']).apply(lambda x: x.total_seconds())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate Speeds\n",
    "data.loc[:,'speed_ms']=data['distance_inst']/data['delta_sec']\n",
    "data.loc[:,'speed_kmh']=data['speed_ms'].apply(lambda x: (x/1000)*3600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use only the speed between same activity IDs. Check that two continuous data points belong to the same\n",
    "# activity_id to calculate mean speed\n",
    "data.loc[1:,'activity_id2']=(data['new_activity_id'].tolist()[0:-1])\n",
    "data.loc[:,'same_activity']=(data['new_activity_id']==data['activity_id2'])*1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter the data so that we use only the same activity ids for the speeds\n",
    "data=data[data['same_activity']==1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set a speed limit for the speeds which is 30kmh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_speed=30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter the speeds that are reasonable\n",
    "data=data[data['speed_kmh']<=max_speed]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_speed=data.groupby(['new_activity_id'])[['speed_kmh']].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEXCAYAAAC6baP3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHGhJREFUeJzt3X20XVV57/HvjyBvvgESLCZgUFKVUoo0glarVBQFVGyrLdTWaKlp78W3aoeidogv5Q68tUUplhoFAasiFau00EpEEXuVl6DIq0oKKUQQQgGRgiDw3D/WPLI5nCR7JeecfU7y/Yxxxl5rrrnWelZ2zn7OnHPtuVJVSJI0rM1GHYAkaXYxcUiSejFxSJJ6MXFIknoxcUiSejFxSJJ6MXFIM0yS/ZKsWsv2u5I8ZTpjGkaSk5P81ajj0NQzcWikkqxMcl+SHcaVX5qkkiwYTWQzV1U9pqqu7btfkgXt33Tztp4kf5fk+0nmrW88SX6Y5JfXd3/NPiYOzQTXAYeNrST5VWDr0YWz8UsS4OPAfsALqupH63mcpwKbVdUPJzE8zXAmDs0EnwZeO7C+GDh1sEKSLZN8OMn1SW5O8g9Jtm7btkvyr0lWJ7m9Lc8f2Pe8JB9M8v+S/DTJOeNbOAN1d2j735HktiTfTLJZ27YyybuSXNXO86kkWw3s+7LWUrojybeS7Dmw7UlJzmgxXpfkzQPbtm7dPLcnuQp41tr+sVqrYbe2fHKSjyU5q13bhe3DfG3mACcDi4D9qurmdqzXtX+jY9s1XJvkN1r5DUluSbJ43LEOBs4eWN+uZyyahUwcmgkuAB6X5BlJ5gC/D/zjuDofAn4Z2AvYDZgHvLdt2wz4FPBkYBfgHuD4cfv/AfB6YEdgC+Av1hDL24FVwFzgicC7gcF5eV4DvAR4aovnLwGS7A2cBPwp8AS6v+bPbAlvM+BfgO+1uPcH3prkJe2YR7XjPbUde/yH87ocBrwf2A5YARy9jvqfAZ4OvLCq/nvctn2By9o1fBY4jS6R7Qb8IXB8kscM1D8IOGsDYtEsZOLQTDHW6ngx8H3gF10nrVvlDcCfV9VtVfVT4P8AhwJU1X9X1RlVdXfbdjTwgnHH/1RV/bCq7gFOp0tAE/k5sBPw5Kr6eVV9sx4+odvxVXVDVd3WzjPWxfYG4ONVdWFVPVBVpwD3As+m++CdW1UfqKr72vjEJ8biB34POLpd2w3AcX3+4YAvVtVFVXU/XVJY07WNOQA4varumGDbdVX1qap6APg8sDPwgaq6t6rOAe6jSyIk2aZd2zc2IBbNQpuPOgCp+TRwPrAr47qp6P763wa4pMshAISuy2XsA+xY4KV0f+kCPDbJnPYBCPDjgePdDQz+1Tzor4H3Aee0cy2tqmMGtt8wsPxfwJPa8pOBxUneNLB9i7b9AeBJSQY/qOcA32zLT5rguH0Me21jXgb8a5Lbq+qkcdtuHli+B2CsK2ugbOz4+wPfqqqfbUAsmoVMHJoRquq/klxH1/Vx+LjNt9J9YP3KGgZx3w48Ddi3qn6cZC/gu3TJpW8cP23He3uSXwG+nuTiqjq3Vdl5oPouwI1t+Qa6VsMjumaSPIfuL/mFazjtTe24Vw4cdyp9C3g5XfL4WVV9dj2PM76bSpsIu6o0kxxO1+/+P4OFVfUgXdfOsUl2BEgyb2CM4LF0ieWOJNvTjRmslzbAvVvrHruTrrXwwECVI5LMb+d5N113Di2+P0uyb7vN9dFJDk7yWOAi4M4k72wD4XOS7JFkbBD8dOBdbZB/PjDYapkSVfUN4HeApUletZ6HOZCHD4xrE2Hi0IxRVf9ZVcvXsPmddIOtFyS5E/gqXSsD4CN0t+/eSjfQ/u8bEMbCduy7gG8Df19V5w1s/yxwDnBt+/mrFvtyunGO44HbW6yva9seoPsLfy+6W49vBT4JPL4d8/103VPXtWN/egPiH1pVLaO7EeHkJC/vs2+SPYC7qur6KQlOM1p8kJM0nCQrgT+pqq+OOpZRS/IOYIeqeseoY9H0c4xD0vpYSXeLsTZBJg5JvVXV6aOOQaNjV5UkqRcHxyVJvWyUXVU77LBDLViwYNRhSNKscskll9xaVXPXVW+jTBwLFixg+fI13dUpSZpIkqFmLbCrSpLUi4lDktSLiUOS1IuJQ5LUi4lDktSLiUOS1MuUJY4kJ7VnFF8xwba/aM9N3qGtJ8lxSVYkuaw9hnOs7uIk17Sfvo/UlCRNsqlscZxM90S2h0myM93jQQenYz6QbjrrhcAS4IRWd+zZCvsC+wBHJdkOSdLITFniqKrzgdsm2HQs8A5gcJKsQ4BTq3MBsG2SnYCXAMvas5hvB5YxQTKSJE2faf3meJJXAD+qqu8NPDsaYB4Pf+byqla2pvKJjr2ErrXCLrtM9ZM3pfW34MjRPG115TEHj+S82vhM2+B4km2A9wDvnWjzBGW1lvJHFlYtrapFVbVo7tx1TrUiSVpP03lX1VOBXYHvtSepzQe+k+SX6FoSOw/UnQ/cuJZySdKITFviqKrLq2rHqlpQVQvoksLeVfVj4Ezgte3uqmcDP6mqm4CvAAck2a4Nih/QyiRJIzKVt+N+Dvg28LQkq5IcvpbqZwPXAiuATwD/G6CqbgM+CFzcfj7QyiRJIzJlg+NVddg6ti8YWC7giDXUOwk4aVKDkyStN785LknqxcQhSerFxCFJ6sXEIUnqxcQhSerFxCFJ6sXEIUnqxcQhSerFxCFJ6sXEIUnqxcQhSerFxCFJ6sXEIUnqxcQhSerFxCFJ6sXEIUnqxcQhSerFxCFJ6sXEIUnqxcQhSerFxCFJ6mXKEkeSk5LckuSKgbK/TvL9JJcl+eck2w5se1eSFUl+kOQlA+UvbWUrkhw5VfFKkoYzlS2Ok4GXjitbBuxRVXsCPwTeBZBkd+BQ4FfaPn+fZE6SOcDHgAOB3YHDWl1J0ohMWeKoqvOB28aVnVNV97fVC4D5bfkQ4LSqureqrgNWAPu0nxVVdW1V3Qec1upKkkZklGMcfwz8W1ueB9wwsG1VK1tT+SMkWZJkeZLlq1evnoJwJUkwosSR5D3A/cBnxoomqFZrKX9kYdXSqlpUVYvmzp07OYFKkh5h8+k+YZLFwMuA/atqLAmsAnYeqDYfuLEtr6lckjQC09riSPJS4J3AK6rq7oFNZwKHJtkyya7AQuAi4GJgYZJdk2xBN4B+5nTGLEl6uClrcST5HLAfsEOSVcBRdHdRbQksSwJwQVX9WVVdmeR04Cq6LqwjquqBdpw3Al8B5gAnVdWVUxWzJGndpixxVNVhExSfuJb6RwNHT1B+NnD2JIYmSdoAfnNcktSLiUOS1IuJQ5LUi4lDktSLiUOS1IuJQ5LUi4lDktSLiUOS1IuJQ5LUi4lDktSLiUOS1IuJQ5LUi4lDktSLiUOS1IuJQ5LUi4lDktSLiUOS1IuJQ5LUi4lDktSLiUOS1MuUJY4kJyW5JckVA2XbJ1mW5Jr2ul0rT5LjkqxIclmSvQf2WdzqX5Nk8VTFK0kazlS2OE4GXjqu7Ejg3KpaCJzb1gEOBBa2nyXACdAlGuAoYF9gH+CosWQjSRqNKUscVXU+cNu44kOAU9ryKcArB8pPrc4FwLZJdgJeAiyrqtuq6nZgGY9MRpKkaTTdYxxPrKqbANrrjq18HnDDQL1VrWxN5Y+QZEmS5UmWr169etIDlyR1ZsrgeCYoq7WUP7KwamlVLaqqRXPnzp3U4CRJD5nuxHFz64Kivd7SylcBOw/Umw/cuJZySdKITHfiOBMYuzNqMfDlgfLXtrurng38pHVlfQU4IMl2bVD8gFYmSRqRzafqwEk+B+wH7JBkFd3dUccApyc5HLgeeHWrfjZwELACuBt4PUBV3Zbkg8DFrd4Hqmr8gLskaRpNWeKoqsPWsGn/CeoWcMQajnMScNIkhiZJ2gAzZXBckjRLmDgkSb2YOCRJvZg4JEm9TNnguDSTLTjyrFGHIM1atjgkSb2YOCRJvZg4JEm9mDgkSb0MlTiS7DHVgUiSZodh76r6hyRb0D3V77NVdcfUhSRpKozyTrKVxxw8snNr8g3V4qiq5wGvoZvifHmSzyZ58ZRGJkmakYYe46iqa4C/BN4JvAA4Lsn3k/zOVAUnSZp5hh3j2DPJscDVwAuBl1fVM9rysVMYnyRphhl2jON44BPAu6vqnrHCqroxyV9OSWSSpBlp2MRxEHBPVT0AkGQzYKuquruqPj1l0UmSZpxhxzi+Cmw9sL5NK5MkbWKGTRxbVdVdYytteZupCUmSNJMNmzj+J8neYytJfh24Zy31JUkbqWHHON4K/FOSG9v6TsDvT01IkqSZbKjEUVUXJ3k68DQgwPer6udTGpkkaUbqM8nhs4A9gWcChyV57fqeNMmfJ7kyyRVJPpdkqyS7JrkwyTVJPt+mOCHJlm19Rdu+YH3PK0nacMN+AfDTwIeB59ElkGcBi9bnhEnmAW8GFlXVHsAc4FDgQ8CxVbUQuB04vO1yOHB7Ve1G92XDD63PeSVJk2PYMY5FwO5VVZN43q2T/Jzu7qyb6L6F/gdt+ynA+4ATgEPaMsAXgOOTZBJjkST1MGxX1RXAL03GCavqR3Stl+vpEsZPgEuAO6rq/lZtFTCvLc8Dbmj73t/qP2H8cZMsSbI8yfLVq1dPRqiSpAkM2+LYAbgqyUXAvWOFVfWKvidMsh1dK2JX4A7gn4ADJ6g61qLIWrY9VFC1FFgKsGjRIlsjkjRFhk0c75vEc74IuK6qVgMk+SLwG8C2STZvrYr5wNitv6vopnNflWRz4PHAbZMYjySph2Gfx/ENYCXwqLZ8MfCd9Tzn9cCzk2yTJMD+wFXA14FXtTqLgS+35TPbOm371xzfkKTRGfauqjfQDUx/vBXNA760Piesqgvbsb4DXN5iWEr3nI+3JVlBN4ZxYtvlROAJrfxtwJHrc15J0uQYtqvqCGAf4ELoHuqUZMf1PWlVHQUcNa742naO8XV/Brx6fc8lSZpcw95VdW9V3Te20sYa7C6SpE3QsInjG0neTffdixfT3Qn1L1MXliRppho2cRwJrKYbk/hT4Gy6549LkjYxw05y+CDdo2M/MbXhSJJmuqESR5LrmPhLd0+Z9IgkSTNan7mqxmxFd5fT9pMfjiRpphv2C4D/PfDzo6r6CN2khJKkTcywXVV7D6xuRtcCeeyURCRJmtGG7ar6m4Hl++mmH/m9SY9GkjTjDXtX1W9NdSCSpNlh2K6qt61te1X97eSEI0ma6frcVfUsuplqAV4OnE97wJIkadPR50FOe1fVTwGSvA/4p6r6k6kKTJI0Mw075cguwH0D6/cBCyY9GknSjDdsi+PTwEVJ/pnuG+S/DZw6ZVFJkmasYe+qOjrJvwG/2YpeX1XfnbqwJEkz1bBdVQDbAHdW1Ufpnv+96xTFJEmawYZ9dOxRdI92fVcrehTwj1MVlCRp5hq2xfHbwCuA/wGoqhtxyhFJ2iQNmzjuq6qiTa2e5NFTF5IkaSYbNnGcnuTjwLZJ3gB8lQ14qFOSbZN8Icn3k1yd5DlJtk+yLMk17XW7VjdJjkuyIsll4yZclCRNs2GnVf8w8AXgDOBpwHur6u824LwfBf69qp4O/BpwNd3jac+tqoXAuW0d4EBgYftZApywAeeVJG2gdd6Om2QO8JWqehGwbENPmORxwPOB1wFU1X3AfUkOAfZr1U4BzqMbkD8EOLV1lV3QWis7VdVNGxqLJKm/dbY4quoB4O4kj5+kcz4FWA18Ksl3k3yyjZk8cSwZtNcdW/15PHxOrFWt7GGSLEmyPMny1atXT1KokqTxhv3m+M+Ay5Mso91ZBVBVb17Pc+4NvKmqLkzyUR7qlppIJiib6PnnS4GlAIsWLXrEdknS5Bg2cZzVfibDKmBVVV3Y1r9AlzhuHuuCSrITcMtA/Z0H9p8P3DhJsUiSelpr4kiyS1VdX1WnTNYJq+rHSW5I8rSq+gGwP3BV+1kMHNNev9x2ORN4Y5LTgH2Bnzi+IUmjs64Wx5foupVIckZV/e4knfdNwGeSbAFcC7yebrzl9CSHA9cDr251zwYOAlYAd7e6kqQRWVfiGBxfeMpknbSqLqV7ONR4+09Qt4AjJuvckqQNs667qmoNy5KkTdS6Why/luROupbH1m2Ztl5V9bgpjU6SNOOsNXFU1ZzpCkSSNDv0eR6HJEkmDklSPyYOSVIvJg5JUi8mDklSLyYOSVIvJg5JUi8mDklSLyYOSVIvJg5JUi8mDklSLyYOSVIvJg5JUi8mDklSLyYOSVIvJg5JUi8mDklSLyYOSVIvI0scSeYk+W6Sf23ruya5MMk1ST6fZItWvmVbX9G2LxhVzJKk0bY43gJcPbD+IeDYqloI3A4c3soPB26vqt2AY1s9SdKIjCRxJJkPHAx8sq0HeCHwhVblFOCVbfmQtk7bvn+rL0kagVG1OD4CvAN4sK0/Abijqu5v66uAeW15HnADQNv+k1b/YZIsSbI8yfLVq1dPZeyStEmb9sSR5GXALVV1yWDxBFVriG0PFVQtrapFVbVo7ty5kxCpJGkim4/gnM8FXpHkIGAr4HF0LZBtk2zeWhXzgRtb/VXAzsCqJJsDjwdum/6wJUkwghZHVb2rquZX1QLgUOBrVfUa4OvAq1q1xcCX2/KZbZ22/WtV9YgWhyRpesyk73G8E3hbkhV0YxgntvITgSe08rcBR44oPkkSo+mq+oWqOg84ry1fC+wzQZ2fAa+e1sAkSWs00sQhLTjyrFGHIKmnmdRVJUmaBUwckqReTBySpF5MHJKkXkwckqReTBySpF5MHJKkXkwckqReTBySpF5MHJKkXkwckqReTBySpF5MHJKkXkwckqReTBySpF5MHJKkXkwckqReTBySpF5MHJKkXqY9cSTZOcnXk1yd5Mokb2nl2ydZluSa9rpdK0+S45KsSHJZkr2nO2ZJ0kNG0eK4H3h7VT0DeDZwRJLdgSOBc6tqIXBuWwc4EFjYfpYAJ0x/yJKkMdOeOKrqpqr6Tlv+KXA1MA84BDilVTsFeGVbPgQ4tToXANsm2Wmaw5YkNSMd40iyAHgmcCHwxKq6CbrkAuzYqs0DbhjYbVUrG3+sJUmWJ1m+evXqqQxbkjZpI0scSR4DnAG8taruXFvVCcrqEQVVS6tqUVUtmjt37mSFKUkaZySJI8mj6JLGZ6rqi6345rEuqPZ6SytfBew8sPt84MbpilWS9HCbT/cJkwQ4Ebi6qv52YNOZwGLgmPb65YHyNyY5DdgX+MlYl5Ymz4Ijzxp1CNqIjer/18pjDh7JeTd20544gOcCfwRcnuTSVvZuuoRxepLDgeuBV7dtZwMHASuAu4HXT2+4kqRB0544quo/mHjcAmD/CeoXcMSUBiVJGtooWhxaA7uLJM0GTjkiSerFxCFJ6sXEIUnqxcQhSerFxCFJ6sXEIUnqxcQhSerFxCFJ6sXEIUnqxcQhSerFxCFJ6sXEIUnqxcQhSerFxCFJ6sXEIUnqxcQhSerFxCFJ6sXEIUnqxcQhSerFZ45PwGd/SxuHUf4urzzm4JGde6rNmhZHkpcm+UGSFUmOHHU8krSpmhWJI8kc4GPAgcDuwGFJdh9tVJK0aZotXVX7ACuq6lqAJKcBhwBXjTQqSVqDUXWTTUcX2WxJHPOAGwbWVwH7DlZIsgRY0lbvSvKDDTjfDsCtG7D/TLIxXQtsXNezMV0LbFzXM2uvJR+asHjY63nyMOeYLYkjE5TVw1aqlgJLJ+VkyfKqWjQZxxq1jelaYOO6no3pWmDjup6N6Vpg8q9nVoxx0LUwdh5Ynw/cOKJYJGmTNlsSx8XAwiS7JtkCOBQ4c8QxSdImaVZ0VVXV/UneCHwFmAOcVFVXTuEpJ6XLa4bYmK4FNq7r2ZiuBTau69mYrgUm+XpSVeuuJUlSM1u6qiRJM4SJQ5LUi4ljwMY2rUmSlUkuT3JpkuWjjqevJCcluSXJFQNl2ydZluSa9rrdKGMc1hqu5X1JftTen0uTHDTKGIeVZOckX09ydZIrk7yllc/W92ZN1zPr3p8kWyW5KMn32rW8v5XvmuTC9t58vt1ktP7ncYyj06Y1+SHwYrrbfy8GDquqWfvt9CQrgUVVNTu/yJQ8H7gLOLWq9mhl/xe4raqOacl9u6p65yjjHMYaruV9wF1V9eFRxtZXkp2AnarqO0keC1wCvBJ4HbPzvVnT9fwes+z9SRLg0VV1V5JHAf8BvAV4G/DFqjotyT8A36uqE9b3PLY4HvKLaU2q6j5gbFoTjUhVnQ/cNq74EOCUtnwK3S/4jLeGa5mVquqmqvpOW/4pcDXd7A6z9b1Z0/XMOtW5q60+qv0U8ELgC618g98bE8dDJprWZFb+5xlQwDlJLmlTsmwMnlhVN0H3Cw/sOOJ4NtQbk1zWurJmRdfOoCQLgGcCF7IRvDfjrgdm4fuTZE6SS4FbgGXAfwJ3VNX9rcoGf7aZOB6yzmlNZqHnVtXedLMKH9G6SzRznAA8FdgLuAn4m9GG00+SxwBnAG+tqjtHHc+GmuB6ZuX7U1UPVNVedDNs7AM8Y6JqG3IOE8dDNrppTarqxvZ6C/DPdP+JZrubW5/0WN/0LSOOZ71V1c3tl/xB4BPMoven9Z+fAXymqr7YimftezPR9czm9wegqu4AzgOeDWybZOwL3xv82WbieMhGNa1Jkke3gT6SPBo4ALhi7XvNCmcCi9vyYuDLI4xlg4x9yDa/zSx5f9oA7InA1VX1twObZuV7s6brmY3vT5K5SbZty1sDL6Ibs/k68KpWbYPfG++qGtBut/sID01rcvSIQ1pvSZ5C18qAbmqZz86260nyOWA/uimhbwaOAr4EnA7sAlwPvLqqZvyg8xquZT+6bpACVgJ/OjZGMJMleR7wTeBy4MFW/G66cYHZ+N6s6XoOY5a9P0n2pBv8nkPXMDi9qj7QPg9OA7YHvgv8YVXdu97nMXFIkvqwq0qS1IuJQ5LUi4lDktSLiUOS1IuJQ5LUi4lDktSLiUMCkrynTUN9WZtCe98pPNd5SRatZftda9o2xLFPTvKqddeU1t+seOa4NJWSPAd4GbB3Vd2bZAdgg55XIG3MbHFIsBNw69g3aavq1qq6sT0I60PtwTgXJdkNfjGtwxlJLm4/z23lj26zqF6c5LtJDmnlWyc5rbVmPg9sPUxQSXZI8u0kByfZL8k3kpye5IdJjknymhbX5UmeOrDr85N8K8m1tj40FUwcEpwD7Nw+kP8+yQsGtt1ZVfsAx9NNRwPwUeDYqnoW8LvAJ1v5e4CvtfLfAv66zRP2v4C7q2pP4Gjg19cVUJInAmcB762qs1rxr9E9lOdXgT8CfrnF9kngTQO77wQ8j64VdUyPfwdpKHZVaZPXnpb268Bv0n3gfz4PPTr4cwOvx7blFwG7d3PjAfC4NqHkAcArkvxFK9+Kbt6m5wPHtXNdluSydYT0KOBc4Iiq+sZA+cVjcyUl+U+6hAfdHEu/NVDvS21G16taApImlYlDonuGAd0U1OcluZyHZnkdnMxtbHkz4DlVdc/gMdosq79bVT8YVz7+OOtyP93jS18CDCaOwUnpHhxYf5CH/y4P1pvoOTPSBrGrSpu8JE9LsnCgaC/gv9ry7w+8frstnwO8cWD/vdriV4A3tQRCkme28vOB17SyPYA91xFSAX8MPH2g5SPNGLY4JHgM8HftOQb3AyuAJXRjBFsmuZDuj6zDWv03Ax9rXU6b0yWGPwM+SDcOcllLHivbMU4APtXqXwpctK6AquqBJIcC/5LkTuCqSbpWaYM5rbq0BklWAouq6tZRxyLNJHZVSZJ6scUhjUjrAttyXPEfVdXlo4hHGpaJQ5LUi11VkqReTBySpF5MHJKkXkwckqRe/j8nva2wGwB4igAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(mean_speed['speed_kmh'],bins=10)\n",
    "plt.xlabel('Speed_kmh')\n",
    "plt.ylabel('Frequency')\n",
    "plt.title('Mean speed in Km/h')\n",
    "plt.savefig('speeds.eps', format='eps', dpi=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Trip Duration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "def duration(data):\n",
    "    minutes=(data['date_1'].max()-data['date_1'].min()).total_seconds()/60\n",
    "    return minutes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "durations=data.groupby('new_activity_id').apply(duration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEWCAYAAACT7WsrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAH5JJREFUeJzt3Xu8VXWd//HXW/CaF0TREDCwqJEuIp6f2XQZlS6CFZSa+jMlpWhm6GLpTHTV5jGPwrnk5Ew/G9IMrTQ0TSat5IeoPysvB0JFyTgRwgmC4wW8pYZ+fn+s79btZq+z14Gzzt7nnPfz8TiPvdZ3fddan+9ecD7n+103RQRmZmb17NTsAMzMrHU5SZiZWS4nCTMzy+UkYWZmuZwkzMwsl5OEmZnlcpKwUkg6WNKTkob00vZC0mt6Y1sF9/czSTP6YD8fkXR71fyTkg7ppW1/QdIlaXps+g6H9tK2e/X4WutykrAek7RG0ju7qxMRayNiz4h4vq/i2l6Szpf0/eqyiJgSEfP7Opb0na3uro6koyV1FtjW1yLio70RV+0x70/H13aMk4T1ut76a7U3tFIsfWmwttt6n5OE9YikK4CDgf9Jww3/WDWUMVPSWuDm2uENSbdI+rqkuyRtkXS9pOHd7OcfJG2QtF7SWTXLbpH00ar52iGbkDRb0ipgVSr7pqR1kh6XtFTS21P5ccAXgJNTe+6p3YeknSR9SdJDkjZJulzSPmlZpZ0zJK2V9LCkL3bTrv0kLUxx3AW8umb5i8NqkqZKekDSE5L+KOlcSa8AfgYclOJ9UtJBqTd0jaTvS3oc+Ei9HhJwVvpON0g6p2q/35P0z1XzL/ZWGhzzyvE9KLXrUUkdkj5Wta3zJS1I39sTku6X1Jb3HVlrcZKwHomI04G1wPvScMO/VC3+G+BQ4D05q58BnAUcBGwFLqpXKf3iPhd4FzAe6HZoK8d04M3AhDR/NzARGA78ELha0m4R8XPga8CPUnsOq7Otj6SfY4BDgD2B/6qp8zbgdcBk4CuSDs2J61vAM8BIsu/irJx6AJcCH4+IvYA3ADdHxFPAFGB9infPiFif6k8DrgGGAT/I2eYxZN/pu4E5jYYNoeExr7gS6CQ7ticCX5M0uWr5+4GrUmwL2fb7sxblJGG96fyIeCoi/pyz/IqIWJF+0X0Z+FDOic8PAZdV1T1/O2L5ekQ8WoklIr4fEY9ExNaI+HdgV7Jf6kWcBnwjIlZHxJPA54FTaoZ0vhoRf46Ie4B7gG2STWrrCcBX0ve0AujuvMdfgAmS9o6IxyJiWYM4fx0RP4mIF7o5Bl9N+74PuAw4tcE2G5I0hixJfi4inomI5cAlwOlV1W6PiBvTOYwrqPP9WGtykrDetK4Hyx8Cdgb2r1PvoDp1dygWSedIWpmGujYD++Tsu56DamJ4CBgKHFhV9qeq6afJehu1RqT1irbtBGAq8JCkWyW9pUGcjb7/2joPkbVtRx0EPBoRT9Rse1TVfO33s5vPm/QPThK2PfIeHdzokcJjqqYPJvtL+eE69TbUqVvtKWCPqvlXdhdLOv/wObIeyr4RMQzYAqhg3OuBV9XEsxXY2GC9Wl1pve7a9qKIuDsipgEHAD8BFjSIt8gjnWv3XRmqavSddrft9cBwSXvVbPuPBeKxFuckYdtjI9nYfE99WNIESXsA/wRck3MJ5QKyE6+VuufVLF8OfFDSHukk78wG+92L7JdzFzBU0leAvWvaM1ZS3v+HK4HPSBonaU9eOoextcF+Xya19Vrg/BT7BKDuvRiSdpF0mqR9IuIvwONA5bvaCOxXOXneQ19O+349cCbwo1S+HJgqabikVwJn16yXe8wjYh3wK+DrknaT9CayY5J3XsT6EScJ2x5fB74kabOkc3uw3hXA98iGHnYDPlWvUkT8DPgP4GagI31WuxB4juwX13wa/zL6BdkVQb8jGwZ5hpcPu1ydPh+RVG/c/7sp9tuAP6T1P9lgn3k+QTYU9Sey7+KybuqeDqxJVyv9LfBhgIj4LVniWp2OQU+GjG4l+04XA/8WETel8ivIzqWsAW7ipeRR0eiYnwqMJetVXAecFxGLehCXtSj5pUPWFyTdAnw/Ii5pdixmVpx7EmZmlstJwszMcnm4yczMcrknYWZmufr1zSz7779/jB07ttlhmJn1K0uXLn04IkYUqduvk8TYsWNpb29vdhhmZv2KpMJPMfBwk5mZ5XKSMDOzXE4SZmaWy0nCzMxyOUmYmVkuJwkzM8vlJGFmZrmcJMzMLJeThJmZ5Sr1jmtJnwE+Svbqw/vI3oQ1ErgKGA4sA06PiOck7QpcDhwBPAKcHBFryoyvqLFzbihUb83c40uOxMysb5XWk5A0iuzNY20R8QZgCHAKcAFwYUSMBx7jpVdPzgQei4jXkL157IKyYjMzs2LKHm4aCuwuaSjZS9Y3AMcC16Tl84HpaXpamictnyxJmJlZ05SWJCLij8C/AWvJksMWYCmwueoF8p3AqDQ9ivTe4bR8C7BfWfGZmVljZQ437UvWOxgHHAS8AphSp2rlrUf1eg3bvBFJ0ixJ7ZLau7q6eitcMzOro8wT1+8E/hARXQCSrgX+GhgmaWjqLYwG1qf6ncAYoDMNT+0DPFq70YiYB8wDaGtr2+7X6hU9GV3WNn2S28z6gzLPSawFjpK0Rzq3MBl4AFgCnJjqzACuT9ML0zxp+c3hd6uamTVVmeck7iQ7Ab2M7PLXnch6AJ8DPiupg+ycw6VplUuB/VL5Z4E5ZcVmZmbFlHqfREScB5xXU7waOLJO3WeAk8qMx8zMesZ3XJuZWS4nCTMzy+UkYWZmuZwkzMwsl5OEmZnlcpIwM7NcThJmZpbLScLMzHI5SZiZWS4nCTMzy+UkYWZmuUp9dpPl83uzzaw/cE/CzMxyOUmYmVkuJwkzM8vlJGFmZrmcJMzMLFdpSULS6yQtr/p5XNLZkoZLWiRpVfrcN9WXpIskdUi6V9KksmIzM7NiynzH9YMRMTEiJgJHAE8D15G9u3pxRIwHFvPSu6ynAOPTzyzg4rJiMzOzYvpquGky8PuIeAiYBsxP5fOB6Wl6GnB5ZO4Ahkka2UfxmZlZHX2VJE4BrkzTB0bEBoD0eUAqHwWsq1qnM5W9jKRZktoltXd1dZUYspmZlZ4kJO0CvB+4ulHVOmWxTUHEvIhoi4i2ESNG9EaIZmaWoy96ElOAZRGxMc1vrAwjpc9NqbwTGFO13mhgfR/EZ2ZmOfoiSZzKS0NNAAuBGWl6BnB9VfkZ6Sqno4AtlWEpMzNrjlIf8CdpD+BdwMeriucCCyTNBNYCJ6XyG4GpQAfZlVBnlhmbmZk1VmqSiIingf1qyh4hu9qptm4As8uMx8zMesZ3XJuZWS4nCTMzy+UkYWZmuZwkzMwsl5OEmZnlcpIwM7NcThJmZpbLScLMzHI5SZiZWS4nCTMzy+UkYWZmuZwkzMwsl5OEmZnlcpIwM7NcThJmZpbLScLMzHKVmiQkDZN0jaTfSlop6S2ShktaJGlV+tw31ZWkiyR1SLpX0qQyYzMzs8bK7kl8E/h5RPwVcBiwEpgDLI6I8cDiNA8wBRiffmYBF5ccm5mZNVBakpC0N/AO4FKAiHguIjYD04D5qdp8YHqangZcHpk7gGGSRpYVn5mZNVZmT+IQoAu4TNJvJF0i6RXAgRGxASB9HpDqjwLWVa3fmcpeRtIsSe2S2ru6ukoM38zMykwSQ4FJwMURcTjwFC8NLdWjOmWxTUHEvIhoi4i2ESNG9E6kZmZWV5lJohPojIg70/w1ZEljY2UYKX1uqqo/pmr90cD6EuMzM7MGSksSEfEnYJ2k16WiycADwEJgRiqbAVyfphcCZ6SrnI4CtlSGpczMrDmGlrz9TwI/kLQLsBo4kywxLZA0E1gLnJTq3ghMBTqAp1NdMzNrolKTREQsB9rqLJpcp24As8uMx8zMesZ3XJuZWS4nCTMzy+UkYWZmuZwkzMwsl5OEmZnlcpIwM7NcThJmZpbLScLMzHI5SZiZWS4nCTMzy+UkYWZmuZwkzMwsl5OEmZnlcpIwM7NcThJmZparYZKQNLwvAjEzs9ZTpCdxp6SrJU2VpNIjMjOzllEkSbwWmAecDnRI+pqk1xbZuKQ1ku6TtFxSeyobLmmRpFXpc99ULkkXSeqQdK+kSdvbKDMz6x0Nk0RkFkXEqcBHgRnAXZJulfSWAvs4JiImRkTlNaZzgMURMR5YnOYBpgDj088s4OIetsXMzHpZkXMS+0n6dOoJnAt8EtgfOAf44XbscxowP03PB6ZXlV+ektIdwDBJI7dj+2Zm1kuKDDf9GtgbmB4Rx0fEtRGxNSLagW83WDeAmyQtlTQrlR0YERsA0ucBqXwUsK5q3c5U9jKSZklql9Te1dVVIHwzM9teQwvUeV1ERL0FEXFBg3XfGhHrJR0ALJL0227q1jspvs1+I2Ie2TkS2tra6sZlZma9o0hP4iZJwyozkvaV9IsiG4+I9elzE3AdcCSwsTKMlD43peqdwJiq1UcD64vsx8zMylGkJzEiIjZXZiLisdQz6JakVwA7RcQTafrdwD8BC8lOfs9Nn9enVRYCn5B0FfBmYEtlWGowGzvnhsJ118w9vsRIzGwwKpIknpd0cESsBZD0KuoMA9VxIHBdurViKPDDiPi5pLuBBZJmAmuBk1L9G4GpQAfwNHBmj1piZma9rkiS+CJwu6Rb0/w7yC5R7VZErAYOq1P+CDC5TnkAswvEY2ZmfaRhkkh//U8CjiI7ufyZiHi49MjMzKzpivQkAHYFHk31J0giIm4rLywzM2sFDZOEpAuAk4H7gRdScQBOEmZmA1yRnsR0snslni07GDMzay1F7pNYDexcdiBmZtZ6ivQkngaWS1oMvNibiIhPlRaVmZm1hCJJYmH6MTOzQabIJbDzJe0OHBwRD/ZBTGZm1iKKPCr8fcBy4OdpfqIk9yzMzAaBIieuzyd7MN9mgIhYDowrMSYzM2sRRZLE1ojYUlPmR3SbmQ0CRU5cr5D0v4EhksYDnwJ+VW5YZmbWCor0JD4JvJ7s8tcrgceBs8sMyszMWkORq5ueJnsS7BfLD8fMzFpJkWc3LaH+a0SPLSUiMzNrGUXOSZxbNb0bcAKwtZxwzMyslRQZblpaU/TLqhcQNSRpCNAO/DEi3itpHHAVMBxYBpweEc9J2hW4HDgCeAQ4OSLWFN2PmZn1viI30w2v+tlf0nuAV/ZgH58GVlbNXwBcGBHjgceAmal8JvBYRLwGuDDVMzOzJipyddNSsp7AUuDXwDm89Iu9W5JGA8cDl6R5AccC16Qq88keRQ4wLc2Tlk9O9c3MrEmKDDftyN3V/wH8I7BXmt8P2BwRlXMancCoND0KWJf2uVXSllT/Za9KlTSL9I7tgw8+eAdCMzOzRopc3fTB7pZHxLU5670X2BQRSyUdXSmut4kCy6r3Nw+YB9DW1uY7v83MSlTk6qaZwF8DN6f5Y4BbgC1kv8TrJgngrcD7JU0luypqb7KexTBJQ1NvYjSwPtXvBMYAnZKGAvuQvVfbzMyapMg5iQAmRMQJEXEC2d3XRMSZEXFW7koRn4+I0RExFjgFuDkiTgOWACemajOA69P0wjRPWn5zRLinYGbWREWSxNiI2FA1vxF47Q7s83PAZyV1kJ1zuDSVXwrsl8o/C8zZgX2YmVkvKDLcdIukX5A9tynIegVLerKTiLiFbIiKiFhN9ujx2jrPACf1ZLtmZlauIlc3fULSB4B3pKJ5EXFduWGZmVkrKNKTgOzO6Cci4v9K2kPSXhHxRJmBmZlZ8xW54/pjZDe3/XcqGgX8pMygzMysNRQ5cT2b7HLWxwEiYhVwQJlBmZlZayiSJJ6NiOcqM+keBl+aamY2CBRJErdK+gKwu6R3AVcD/1NuWGZm1gqKJIk5QBdwH/Bx4EbgS2UGZWZmraHbq5vSuyDmR8SHge/0TUhmZtYquu1JRMTzwAhJu/RRPGZm1kKK3CexhuxtdAuBpyqFEfGNsoIyM7PWkNuTkHRFmjwZ+Gmqu1fVj5mZDXDd9SSOkPQqYC3wn30Uj5mZtZDuksS3gZ8D48heX1ohsvskDikxLjMzawG5w00RcVFEHApcFhGHVP2MiwgnCDOzQaDhfRIR8Xd9EYiZmbWeIjfTmZnZIFVakpC0m6S7JN0j6X5JX03l4yTdKWmVpB9V7sGQtGua70jLx5YVm5mZFVNmT+JZ4NiIOAyYCBwn6SjgAuDCiBgPPAbMTPVnAo9FxGuAC1M9MzNrotKSRGSeTLM7p58AjiV7PwXAfGB6mp6W5knLJ0tSWfGZmVljpZ6TkDRE0nJgE7AI+D2wOSK2piqdZC8xIn2uA0jLtwD7lRmfmZl1r9QkERHPR8REYDRwJHBovWrps16vYZv3VkiaJaldUntXV1fvBWtmZtvok6ubImIzcAtwFDAsvbgIsuSxPk13AmPgxRcb7QM8Wmdb8yKiLSLaRowYUXboZmaDWplXN42QNCxN7w68E1gJLAFOTNVmANen6YVpnrT85ojwG/DMzJqoyFNgt9dIYH56J8VOwIKI+KmkB4CrJP0z8Bvg0lT/UuAKSR1kPYhTSoxtQBo754ZC9dbMPb7kSMxsoCgtSUTEvcDhdcpXk52fqC1/BjiprHjMzKznfMe1mZnlcpIwM7NcThJmZpbLScLMzHI5SZiZWS4nCTMzy+UkYWZmuZwkzMwsl5OEmZnlcpIwM7NcThJmZpbLScLMzHI5SZiZWS4nCTMzy+UkYWZmuZwkzMwsl5OEmZnlKvMd12MkLZG0UtL9kj6dyodLWiRpVfrcN5VL0kWSOiTdK2lSWbGZmVkxZfYktgLnRMShwFHAbEkTgDnA4ogYDyxO8wBTgPHpZxZwcYmxmZlZAaUliYjYEBHL0vQTwEpgFDANmJ+qzQemp+lpwOWRuQMYJmlkWfGZmVljfXJOQtJY4HDgTuDAiNgAWSIBDkjVRgHrqlbrTGW125olqV1Se1dXV5lhm5kNeqUnCUl7Aj8Gzo6Ix7urWqcstimImBcRbRHRNmLEiN4K08zM6ig1SUjamSxB/CAirk3FGyvDSOlzUyrvBMZUrT4aWF9mfGZm1r0yr24ScCmwMiK+UbVoITAjTc8Arq8qPyNd5XQUsKUyLGVmZs0xtMRtvxU4HbhP0vJU9gVgLrBA0kxgLXBSWnYjMBXoAJ4GziwxNjMzK6C0JBERt1P/PAPA5Dr1A5hdVjxmZtZzvuPazMxyOUmYmVkuJwkzM8vlJGFmZrmcJMzMLJeThJmZ5XKSMDOzXE4SZmaWy0nCzMxylflYDmtRY+fcULjumrnHlxiJmbU69yTMzCyXk4SZmeVykjAzs1xOEmZmlstJwszMcjlJmJlZrjJfX/pdSZskragqGy5pkaRV6XPfVC5JF0nqkHSvpEllxWVmZsWV2ZP4HnBcTdkcYHFEjAcWp3mAKcD49DMLuLjEuMzMrKDSkkRE3AY8WlM8DZifpucD06vKL4/MHcAwSSPLis3MzIrp63MSB0bEBoD0eUAqHwWsq6rXmcq2IWmWpHZJ7V1dXaUGa2Y22LXKiWvVKYt6FSNiXkS0RUTbiBEjSg7LzGxw6+sksbEyjJQ+N6XyTmBMVb3RwPo+js3MzGr0dZJYCMxI0zOA66vKz0hXOR0FbKkMS5mZWfOU9hRYSVcCRwP7S+oEzgPmAgskzQTWAiel6jcCU4EO4GngzLLiMjOz4kpLEhFxas6iyXXqBjC7rFjMzGz7tMqJazMza0FOEmZmlstJwszMcjlJmJlZLr/j2rpV9H3Yfhe22cDknoSZmeVykjAzs1xOEmZmlstJwszMcjlJmJlZLicJMzPL5SRhZma5nCTMzCyXk4SZmeXyHdfWK4remQ2+O9usP3GSsD7nR32Y9R8ebjIzs1wt1ZOQdBzwTWAIcElEzG1ySNZEHsIya76W6UlIGgJ8C5gCTABOlTShuVGZmQ1urdSTOBLoiIjVAJKuAqYBDzQ1KusXetLrKKpo78Q9HhvIWilJjALWVc13Am+urSRpFjArzT4p6cHt3N/+wMPbuW6rGmhtamp7dEEp2/Qxan0DrU312vOqoiu3UpJQnbLYpiBiHjBvh3cmtUdE245up5UMtDYNtPbAwGvTQGsPDLw27Wh7WuacBFnPYUzV/GhgfZNiMTMzWitJ3A2MlzRO0i7AKcDCJsdkZjaotcxwU0RslfQJ4Bdkl8B+NyLuL3GXOzxk1YIGWpsGWntg4LVpoLUHBl6bdqg9ithm2N/MzAxoreEmMzNrMU4SZmaWa1AmCUnHSXpQUoekOc2OZ3tIWiPpPknLJbWnsuGSFklalT73bXac3ZH0XUmbJK2oKqvbBmUuSsfsXkmTmhd5fTntOV/SH9NxWi5patWyz6f2PCjpPc2JunuSxkhaImmlpPslfTqV98vj1E17+u1xkrSbpLsk3ZPa9NVUPk7SnekY/ShdEISkXdN8R1o+ttsdRMSg+iE7Kf574BBgF+AeYEKz49qOdqwB9q8p+xdgTpqeA1zQ7DgbtOEdwCRgRaM2AFOBn5HdT3MUcGez4y/YnvOBc+vUnZD+7e0KjEv/Joc0uw114hwJTErTewG/S7H3y+PUTXv67XFK3/WeaXpn4M703S8ATknl3wb+Lk3/PfDtNH0K8KPutj8YexIvPv4jIp4DKo//GAimAfPT9HxgehNjaSgibgMerSnOa8M04PLI3AEMkzSybyItJqc9eaYBV0XEsxHxB6CD7N9mS4mIDRGxLE0/AawkezpCvzxO3bQnT8sfp/RdP5lmd04/ARwLXJPKa49R5dhdA0yWVO9mZmBwDjfVe/xHd/9IWlUAN0lamh5VAnBgRGyA7D8DcEDTott+eW3oz8ftE2no5btVQ4D9rj1pWOJwsr9U+/1xqmkP9OPjJGmIpOXAJmARWY9nc0RsTVWq436xTWn5FmC/vG0PxiRR6PEf/cBbI2IS2VNzZ0t6R7MDKll/PW4XA68GJgIbgH9P5f2qPZL2BH4MnB0Rj3dXtU5Zy7WrTnv69XGKiOcjYiLZkyqOBA6tVy199qhNgzFJDIjHf0TE+vS5CbiO7B/GxkrXPn1ual6E2y2vDf3yuEXExvQf+AXgO7w0VNFv2iNpZ7JfqD+IiGtTcb89TvXaMxCOE0BEbAZuITsnMUxS5Ybp6rhfbFNavg/dDJMOxiTR7x//IekVkvaqTAPvBlaQtWNGqjYDuL45Ee6QvDYsBM5IV88cBWypDHe0sprx+A+QHSfI2nNKutJkHDAeuKuv42skjVVfCqyMiG9ULeqXxymvPf35OEkaIWlYmt4deCfZuZYlwImpWu0xqhy7E4GbI53FrqvZZ+ab8UN2BcbvyMbtvtjseLYj/kPIrri4B7i/0gayccXFwKr0ObzZsTZox5VkXfu/kP11MzOvDWRd5G+lY3Yf0Nbs+Au254oU773pP+fIqvpfTO15EJjS7Phz2vQ2sqGIe4Hl6Wdqfz1O3bSn3x4n4E3Ab1LsK4CvpPJDyBJaB3A1sGsq3y3Nd6Tlh3S3fT+Ww8zMcg3G4SYzMyvIScLMzHI5SZiZWS4nCTMzy+UkYWZmuZwkrM9IGibp7xvU+dUObP9oST/d3vXrbO9sSXtUzd9YuR59B7f7t5LO2NHt9PZ2JX2hN+OxgcGXwFqfSc/K+WlEvKHOsiER8fwObv9osid5vrdgfZH9H3ghZ/kasuv8H96RuPoLSU9GxJ7NjsNai3sS1pfmAq9Oz+v/1/SX/xJJPyS7kQlJT6bPoyXdJuk6SQ9I+rakbf69Kns3yG8l3Q58sKr8fEnnVs2vkDQ2/ayU9H+AZcAYSRdLaq95Fv+ngIOAJZKWpLI1kvZP059N21wh6exUVtn2d9K2bkp3wNbG/GJskm6RdIGy9wH8TtLb69Q/WtKtkhakOnMlnZbWuU/Sq4tuV9JHJP1X1bZ/mrY/F9g9HZsfpGUfTusvl/Tfyh4iN0TS91K775P0maIH3/onJwnrS3OA30fExIj4h1R2JNkd4xPq1D8SOAd4I9nD1z5YvVDSbmTP2Xkf8HbglQXjeB3Z46wPj4iH0v7byO5c/RtJb4qIi8iedXNMRBxTs98jgDOBN5M9I+djkg5Pi8cD34qI1wObgRMKxDM0Io4EzgbOy6lzGPBpsu/idOC1aZ1LgE/uwHYBiIg5wJ/TsTlN0qHAyWQPkpwIPA+cRvYAvFER8YaIeCNwWYH2WT/mJGHNdldkz+nPW7Y6DUNdSfZIhWp/BfwhIlZFNm76/YL7fCiydx1UfEjSMrJHG7ye7EUz3XkbcF1EPBXZc/yvJUtSpHiWp+mlwNgC8VQemtdd/bsjexfCs2SPiLgpld/XzTpFtptnMnAEcLeyR1BPJnvMw2rgEEn/Kek4oLsnwtoAMLRxFbNSPdXNstoTZvVOoOWdVNvKy/8I2q3ePtND284F/ldEPCbpezV168l9QQvwbNX088A2w03drPM8+f8nq7f7QtX8CwXWqd5ud99LNQHzI+Lz2yyQDgPeA8wGPgSclbMNGwDck7C+9ATZKyOLOlLZ03p3Ihv6uL1m+W+BcZUxeeDUqmVryF4lirL3LI/L2cfeZElji6QDyd7P0Sje24DpkvZQ9hTeDwD/r3CrmmsNMFHSTpLG8PK3rP1F2WO0IXto34mSDoAX32n9qnROZqeI+DHwZdJ3bAOXexLWZyLiEUm/lLSC7D3INzRY5ddkJ7vfSPaL+bqa7T2j7K18N0h6mCyJVK6c+jHZI6uXkz0e/nc5Md0j6TdkT9NdDfyyavE84GeSNlSfl4iIZanHUXlk9CUR8Rs1eqF8a/gl8AeyYaoVZCfvK+YB90pals5LfIns7Yc7kT3ZdjbwZ+CyqosItulp2MDiS2CtJfX0clYzK4eHm8zMLJd7EmZmlss9CTMzy+UkYWZmuZwkzMwsl5OEmZnlcpIwM7Nc/x8FP1AkDhsYYAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Filter the trips of more than 5 hours\n",
    "plt.hist(durations[durations<300],bins = list(range(0,300,10)));\n",
    "plt.xlabel('trip duration in minutes');\n",
    "plt.ylabel('frequency');\n",
    "plt.title('trip duration distribution');\n",
    "plt.savefig('duration.eps', format='eps', dpi=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "count    4179.000000\n",
       "mean       28.407442\n",
       "std        20.033622\n",
       "min         0.000000\n",
       "25%        13.250000\n",
       "50%        25.833333\n",
       "75%        39.958333\n",
       "max       185.933333\n",
       "dtype: float64"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "durations.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_filtered=pd.DataFrame(durations[durations<300].reset_index())['new_activity_id'].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Thanks to the parsing, we did not get any trip longer than 300 minutes, thus, buth lists have equal lenght\n",
    "len(data_filtered),len(data['new_activity_id'].unique())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Departures analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "data=data[data['new_activity_id'].isin(data_filtered)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4179, 4179)"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(data['new_activity_id'].unique()),len(data['new_activity_id'].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.loc[:,'hour']=data['date'].apply(lambda x: datetime.strptime(x[11:],'%H:%M:%S').hour+datetime.strptime(x[11:],\\\n",
    "    '%H:%M:%S').minute/60.0 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.loc[:,'day_num']=data['date_1'].apply(lambda x: x.weekday())\n",
    "data.loc[:,'day_name']=data['date_1'].apply(lambda x: calendar.day_name[x.weekday()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_departures=data.drop_duplicates(subset='new_activity_id',keep='first')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEWCAYAAACT7WsrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAG0JJREFUeJzt3XuYJVV57/HvT65ykVEHlfuIokY9UXEUEqLgJUZERc8RlZgECEpMMNHo8xyQeJAYNWhMwASjwWBEvOIdNTdEwHhBHBBBJYaJDjAOwij3mwZ4zx+1WrbNru7dw+zpPd3fz/Psp6tWrap6V1d3vXutql07VYUkScPcZ74DkCRNLpOEJKmXSUKS1MskIUnqZZKQJPUySUiSepkktNFK8i9JDpnnGPZJclmSm5O8YAzbPzTJV9b3dkfY7x8mubq164Ej1J+XONdVklVJnjnfcWwMTBIbgfYHfVuSm5Jcn+RrSV6ZZF6OX5JlSSrJphtwn8cl+eBgWVXtX1WnbqgYerwJOKmqtqmqz8xzLOtFks2AvwGe1dr102nLN/jx1/wxSWw8nldV2wK7AccDRwGnbOgg1seJIZ2F8re3G/DddVlxgk+yDwa2ZB3bpYVlofyjLhpVdUNVnQG8BDgkyWMBkmyR5B1JrmjDBO9Jct+2bL8kq5Mck+QnrWfysqltJjkgybeS3JjkyiTHDSybetd4eJIrgC8BX26Lr2/DEb82/Z3+9HebSc5J8pYkXwVuBXZPsl2SU5JcleRHSd6cZJPpbU7ybOAY4CVtf98e2ObL2/ShSb6a5ITW2/pBkl9v5VcmuWZwaGqW39fSJJ9v27k2yX8MS2pJ/hvYHfhci2uLJDsmOaOttzLJKwbqH5fkE0k+mORG4NAh23xgW//GJOcDD5u2/J2tPTcmuSDJU1r5Q5LcOjg0lOSJSda2nsH0/WyR5MQka9rrxFb2COD7A8f3S9PXZcjxH9juO5Jcl+SHSfYfKB/1WG+Zrte8tM2/IckdSe7X5t+c5MSBNgw9hm35c5NclLt73786pC0keVSL96XDli92JomNVFWdD6wGntKK3gY8Ang88HBgJ+DYgVUeAixt5YcAJyd5ZFt2C/B7wBLgAOAPc8/x9X2BXwF+C3hqK1vShiO+PmLYvwscAWwLXA6cCtzR4n0C8Czg5UPa+q/AW4GPtf09rmf7ewEXAw8EPgx8FHhS2/7vACcl2abVnen39Tq63+32dO+qjwHu8fyaqnoYcAVdL2+bqvoZ8JG27o7Ai4C3JnnGwGoHAp+g+11/aEgb3gXcDuwA/H57Dfpmi/kBrY0fT7JlVf0YOAd48UDd3wE+WlX/M2Q/fwbs3bb1OODJwBuq6r+Ax7Q6S6rq6UPW7Tv+e9ElmKXA24FTkqQtG/VY397auO/Avi4H9hmYP7dN9x7DJHsC7wP+gO7v4R+AM5JsMbi/Vu/fgT+uqo8OaauqyteEv4BVwDOHlJ9H988euhP9wwaW/Rrwwza9H90/6NYDy08H/l/P/k4ETmjTy+hOkLsPLJ8q23Sg7Djgg3116E5gbxpY/mDgZ8B9B8oOBs7uiemXtj+wzZe36UOBywaW/a+2/wcPlP2U7oQy2+/rTcBngYfP5dgAuwB3AtsOLP9L4P0DbfjyDNvaBPgf4FEDZW8FvjLDOtcBj2vTLwG+OrCtHwNP7lnvv4HnDMz/FrCq7/hOW3fY8T8UWDkwv1Wr85B1ONZ/AfwtsGlrw6vphli3BG6jS0KzHcN3A38xbbvfB/YdOG5/TpfQnzaO/9uF8prUMVGNZifgWrp3vFsBF9z9xo3QnSimXFdVtwzMX073bpcke9H9Ez4W2BzYAvj4tH1duR7iHdzGbsBmwFUDMd/nXu7n6oHp2wCqanrZNsz++/oruhP6v7flJ1fV8SPsf0fg2qq6aaDscmD5wPxM7due7sQ4WOfywQpJXkf3DnxHupPw/ehOmtAltvck2Z3uHfYN1fU4+2Id3PYv/h7uhR9PTVTVre13tw1dr2cux/pcugvnewKXAGfSXX/bmy4R/STJg5j5GO5GNxz7xwPb3ZxfbuMrgXOr6uw5t3QRcbhpI5XkSXRJ4ivAT+hOgI+pqiXttV1VbTOwyv2TbD0wvyuwpk1/GDgD2KWqtgPeQ/cPN6h6pqfcQvdPO+UhQ+oMrncl3bvLpQMx36+qHjNkvb59rqsZf19VdVNVva6qdgeeB7x22pBRnzXAA5JsO1C2K/CjgfmZ2rGWrse3y7T1AWjXH46iG1K6f1UtAW6gHavqhmpOB15GN7R32iyx7jZtP2t66k4312Mx12P9NeCRwAvpTuLfa/EdwN1DTbP9zV8JvGVg2ZKq2qqqPjKwn1cCuyY5YY7tWVRMEhuZJPdL8ly68fYPVtUlVXUX8F7ghPYOiyQ7Jfmtaav/eZLN28nmudzdW9iW7h3w7UmeDPz2LGGsBe6iu2g75SLgqUl2TbId8PqZNlBVV9GNBf91a9N9kjwsyb49q1wNLBt2AXmuZvt9tQueD2/j6TfSDSHdOcJ2r6Q7wf1luwD7q8DhDL/2MGz9O4FPAccl2SrJo+muH03Zli6JrAU2TXIsXU9i0Afohn6eD3yQfh8B3pBk+3aR+NhZ6g8advx7zfVYV9WtwAXAkdydFL5Gd33h3FZntr/59wKvTLJXOlunu0FjMIHfBDyb7u92lJ7iomSS2Hh8LslNdO+Q/oyuO37YwPKjgJXAeenunPki3buxKT+mG79eQ3fSemVV/Wdb9kfAm9r2j6V7N9qr/RO/Bfhqu3Nk76o6E/gY3YXjC4DPj9Cm36MbAvhei+0TdBdsh5lKaD9NcuEI257NTL+vPdr8zcDXgb+vqnNG3O7BdGP2a4BPA29sv5tRvYpuiObHwPuBfxpY9m/AvwD/RTc8dDvThmyq6qt0J/ALq2rVDPt5M7CC7nhdAlzYymY17PiPsNpcjjV0yWAz4PyB+W25+84qmOEYVtUK4BXASW1/KxlyN1lVXQ/8JrB/kr8YoR2LTqr80qGFLsl+dL2Onec7Fo1fu231w1X1j/MdizZ+XriWFpB2rWpPulttpXvN4SZpgUhyKt2Qy2um3WElrTOHmyRJvexJSJJ6bdTXJJYuXVrLli2b7zAkaaNywQUX/KSqth+l7kadJJYtW8aKFSvmOwxJ2qgkuXz2Wh2HmyRJvUwSkqReJglJUi+ThCSpl0lCktTLJCFJ6mWSkCT1MklIknqZJCRJvTbqT1xLmkzLjv7CnNdZdfwBY4hE95Y9CUlSL5OEJKmXSUKS1MskIUnqZZKQJPUySUiSepkkJEm9TBKSpF4mCUlSL5OEJKmXSUKS1MskIUnqZZKQJPXyKbAaiU/1lBYnexKSpF4mCUlSL5OEJKmXSUKS1MskIUnqZZKQJPUySUiSepkkJEm9/DCdxsYP4EkbP3sSkqReY08SSTZJ8q0kn2/zD03yjSSXJflYks1b+RZtfmVbvmzcsUmSZrYhehKvBi4dmH8bcEJV7QFcBxzeyg8HrquqhwMntHqSpHk01iSRZGfgAOAf23yApwOfaFVOBV7Qpg9s87Tlz2j1JUnzZNw9iROB/wvc1eYfCFxfVXe0+dXATm16J+BKgLb8hlb/lyQ5IsmKJCvWrl07ztgladEbW5JI8lzgmqq6YLB4SNUaYdndBVUnV9Xyqlq+/fbbr4dIJUl9xnkL7D7A85M8B9gSuB9dz2JJkk1bb2FnYE2rvxrYBVidZFNgO+DaMcYnSZrF2HoSVfX6qtq5qpYBLwW+VFUvA84GXtSqHQJ8tk2f0eZpy79UVffoSUiSNpz5+JzEUcBrk6yku+ZwSis/BXhgK38tcPQ8xCZJGrBBPnFdVecA57TpHwBPHlLnduCgDRGPJGk0fuJaktTLJCFJ6mWSkCT1MklIknqZJCRJvUwSkqReJglJUi+ThCSpl0lCktTLJCFJ6mWSkCT1MklIknqZJCRJvTbIU2AlaTbLjv7CnOqvOv6AMUWiQfYkJEm9TBKSpF4mCUlSL5OEJKmXSUKS1MskIUnqZZKQJPUySUiSepkkJEm9TBKSpF4mCUlSL5OEJKmXSUKS1MskIUnqZZKQJPUySUiSepkkJEm9TBKSpF4mCUlSL5OEJKmXSUKS1MskIUnqNbYkkWTLJOcn+XaS7yb581b+0CTfSHJZko8l2byVb9HmV7bly8YVmyRpNOPsSfwMeHpVPQ54PPDsJHsDbwNOqKo9gOuAw1v9w4HrqurhwAmtniRpHo0tSVTn5ja7WXsV8HTgE638VOAFbfrANk9b/owkGVd8kqTZjfWaRJJNklwEXAOcCfw3cH1V3dGqrAZ2atM7AVcCtOU3AA8css0jkqxIsmLt2rXjDF+SFr2RkkSSx67Lxqvqzqp6PLAz8GTgV4ZVm9rNDMsGt3lyVS2vquXbb7/9uoQlSRrRqD2J97SL0H+UZMlcd1JV1wPnAHsDS5Js2hbtDKxp06uBXQDa8u2Aa+e6L0nS+jNSkqiq3wBeRncSX5Hkw0l+c6Z1kmw/lVCS3Bd4JnApcDbwolbtEOCzbfqMNk9b/qWqukdPQpK04Ww6e5VOVV2W5A3ACuBvgSe0C8vHVNWnhqyyA3Bqkk3oktHpVfX5JN8DPprkzcC3gFNa/VOA05KspOtBvHSdWyVJWi9GShJJfhU4DDiA7gL086rqwiQ7Al8H7pEkqupi4AlDyn9Ad31ievntwEFzil6SNFaj9iROAt5L12u4baqwqta03oUkaQEaNUk8B7itqu4ESHIfYMuqurWqThtbdJKkeTXq3U1fBO47ML9VK5MkLWCjJoktBz49TZveajwhSZImxahJ4pYke07NJHkicNsM9SVJC8Co1yReA3w8ydQH33YAXjKekCRJk2KkJFFV30zyKOCRdI/P+M+q+p+xRiZJmncjf5gOeBKwrK3zhCRU1QfGEpUkaSKM+mG604CHARcBd7biAkwSkrSAjdqTWA482mcpSdLiMurdTd8BHjLOQCRJk2fUnsRS4HtJzqf7WlIAqur5Y4lKkjQRRk0Sx40zCEnSZBr1Fthzk+wG7FFVX0yyFbDJeEOTJM23Ub++9BXAJ4B/aEU7AZ8ZV1CSpMkw6oXrI4F9gBuh+wIi4EHjCkqSNBlGTRI/q6qfT82076D2dlhJWuBGTRLnJjkGuG/7buuPA58bX1iSpEkwapI4GlgLXAL8AfDPgN9IJ0kL3Kh3N91F9/Wl7x1vOJKkSTLqs5t+yJBrEFW1+3qPSJI0Meby7KYpWwIHAQ9Y/+FIkibJSNckquqnA68fVdWJwNPHHJskaZ6NOty058Dsfeh6FtuOJSJJ0sQYdbjprwem7wBWAS9e79FIkibKqHc3PW3cgUiSJs+ow02vnWl5Vf3N+glHkjRJ5nJ305OAM9r884AvA1eOIyhJ0mSYy5cO7VlVNwEkOQ74eFW9fFyBSZLm36iP5dgV+PnA/M+BZes9GknSRBm1J3EacH6ST9N98vqFwAfGFpUkaSKMenfTW5L8C/CUVnRYVX1rfGFJkibBqMNNAFsBN1bVO4HVSR46ppgkSRNi1K8vfSNwFPD6VrQZ8MFxBSVJmgyj9iReCDwfuAWgqtbgYzkkacEb9cL1z6uqkhRAkq1nWyHJLnQXtx8C3AWcXFXvTPIA4GN0d0etAl5cVdclCfBO4DnArcChVXXhHNsjaT1bdvQX5jsEzaNRexKnJ/kHYEmSVwBfZPYvILoDeF1V/QqwN3BkkkfTfcvdWVW1B3BWmwfYH9ijvY4A3j2nlkiS1rtR7256R/tu6xuBRwLHVtWZs6xzFXBVm74pyaXATsCBwH6t2qnAOXTXOw4EPlBVBZyXZEmSHdp2JEnzYNYkkWQT4N+q6pnAjIlhhm0sA54AfAN48NSJv6quSvKgVm0nfvkxH6tbmUlCkubJrEmiqu5McmuS7arqhrnuIMk2wCeB11TVjd2lh+FVh+1+yPaOoBuOYtddd51rOJpwcx3/XnX8AWOKRBKMfuH6duCSJGfS7nACqKo/mWmlJJvRJYgPVdWnWvHVU8NISXYArmnlq4FdBlbfGVgzfZtVdTJwMsDy5cvvkUQkSevPqEniC+01sna30inApdMeJX4GcAhwfPv52YHyVyX5KLAXcIPXIyRpfs2YJJLsWlVXVNWp67DtfYDfpeuBXNTKjqFLDqcnORy4AjioLftnuttfV9LdAnvYOuxTkrQezdaT+AywJ0CST1bV/xl1w1X1FYZfZwB4xpD6BRw56vYlSeM32+ckBk/yu48zEEnS5JktSVTPtCRpEZhtuOlxSW6k61Hct03T5quq7jfW6CRJ82rGJFFVm2yoQCRJk2cu3ychSVpkTBKSpF4mCUlSL5OEJKmXSUKS1MskIUnqZZKQJPUySUiSepkkJEm9TBKSpF4mCUlSL5OEJKmXSUKS1MskIUnqZZKQJPUySUiSepkkJEm9Zvv6UmmiLTv6C3NeZ9XxB4whEmlhsichSeplkpAk9XK4SdJGyaHGDcMksUityz+YRucJTAuFSUKaECYWTSKThDSCSe15mVg0biYJaZGZ1ISnyeTdTZKkXiYJSVIvk4QkqZfXJLToOCYvjc6ehCSpl0lCktTLJCFJ6jW2JJHkfUmuSfKdgbIHJDkzyWXt5/1beZL8bZKVSS5Osue44pIkjW6cPYn3A8+eVnY0cFZV7QGc1eYB9gf2aK8jgHePMS5J0ojGliSq6svAtdOKDwRObdOnAi8YKP9Adc4DliTZYVyxSZJGs6GvSTy4qq4CaD8f1Mp3Aq4cqLe6ld1DkiOSrEiyYu3atWMNVpIWu0n5nESGlNWwilV1MnAywPLly4fWkaRhfCDi3G3onsTVU8NI7ec1rXw1sMtAvZ2BNRs4NknSNBs6SZwBHNKmDwE+O1D+e+0up72BG6aGpSRJ82dsw01JPgLsByxNshp4I3A8cHqSw4ErgINa9X8GngOsBG4FDhtXXJKk0Y0tSVTVwT2LnjGkbgFHjisWSdK6mZQL15I0keZ6sXuhXej2sRySpF4mCUlSL5OEJKmXSUKS1MskIUnqZZKQJPUySUiSepkkJEm9TBKSpF4mCUlSL5OEJKmXSUKS1MskIUnqZZKQJPUySUiSepkkJEm9TBKSpF4mCUlSL5OEJKmX33EtSevRXL8TGyb7e7HtSUiSepkkJEm9TBKSpF4mCUlSL5OEJKmXSUKS1MskIUnqZZKQJPUySUiSepkkJEm9TBKSpF4mCUlSLx/wtwCsywPFJGkU9iQkSb3sSUjSPJvkx4tPVE8iybOTfD/JyiRHz3c8krTYTUySSLIJ8C5gf+DRwMFJHj2/UUnS4jZJw01PBlZW1Q8AknwUOBD43jh25sVeSZrdJCWJnYArB+ZXA3tNr5TkCOCINntzku+vw76WAj9Zh/UWCtu/eNu/mNsOC6j9eds6rTbV/t1GXWGSkkSGlNU9CqpOBk6+VztKVlTV8nuzjY2Z7V+87V/MbQfbvy7tn5hrEnQ9h10G5ncG1sxTLJIkJitJfBPYI8lDk2wOvBQ4Y55jkqRFbWKGm6rqjiSvAv4N2AR4X1V9d0y7u1fDVQuA7V+8FnPbwfbPuf2pusewvyRJwGQNN0mSJoxJQpLUa9ElicX+6I8kq5JckuSiJCvmO55xS/K+JNck+c5A2QOSnJnksvbz/vMZ47j0tP24JD9qx/+iJM+ZzxjHKckuSc5OcmmS7yZ5dStf8Md/hrbP+fgvqmsS7dEf/wX8Jt0tt98EDq6qsXyqexIlWQUsr6oF8YGi2SR5KnAz8IGqemwreztwbVUd394o3L+qjprPOMehp+3HATdX1TvmM7YNIckOwA5VdWGSbYELgBcAh7LAj/8MbX8xczz+i60n8YtHf1TVz4GpR39ogaqqLwPXTis+EDi1TZ9K98+z4PS0fdGoqquq6sI2fRNwKd2THRb88Z+h7XO22JLEsEd/rNMvbiNWwL8nuaA94mQxenBVXQXdPxPwoHmOZ0N7VZKL23DUghtqGSbJMuAJwDdYZMd/Wtthjsd/sSWJkR79scDtU1V70j1t98g2JKHF493Aw4DHA1cBfz2/4Yxfkm2ATwKvqaob5zueDWlI2+d8/Bdbklj0j/6oqjXt5zXAp+mG4Babq9uY7dTY7TXzHM8GU1VXV9WdVXUX8F4W+PFPshndSfJDVfWpVrwojv+wtq/L8V9sSWJRP/ojydbtIhZJtgaeBXxn5rUWpDOAQ9r0IcBn5zGWDWrq5Ni8kAV8/JMEOAW4tKr+ZmDRgj/+fW1fl+O/qO5uAmi3fJ3I3Y/+eMs8h7TBJNmdrvcA3SNZPrzQ25/kI8B+dI9Ivhp4I/AZ4HRgV+AK4KCqWnAXeHvavh/dUEMBq4A/mBqfX2iS/AbwH8AlwF2t+Bi6sfkFffxnaPvBzPH4L7okIUka3WIbbpIkzYFJQpLUyyQhSeplkpAk9TJJSJJ6mSS00Uty87T5Q5OctIFjOKg9cfPsaeXLkvz2OGJr216wn3PQZDBJSD3aU4NHdTjwR1X1tGnly4Dfvmd1aeNgktCClmS3JGe1B5qdlWTXVv7+JC8aqHdz+7lfew7/h+k+iDR9ewe37+P4TpK3tbJjgd8A3pPkr6atcjzwlPbs/j9tZTsm+df2fQZvH9j2s5J8PcmFST7enrszff9PTPLtJF8HjhwoX5bkP9q6Fyb59VZ+WpIDB+p9KMnz5/hr1GJWVb58bdQv4E7gooHXFcBJbdnngEPa9O8Dn2nT7wdeNLCNm9vP/YBbgIcO2c+Obdvb031i/UvAC9qyc+i+p2P6OvsBnx+YPxT4AbAdsCVwOd3zxJYCXwa2bvWOAo4dsr2LgX3b9F8B32nTWwFbtuk9gBVtet+BNm8H/BDYdL6Pma+N57XpuqcXaWLcVlWPn5pJciiwvM3+GvC/2/RpwNuZ3flV9cMh5U8CzqmqtW0/HwKeSveYj7k4q6puaNv4HrAbsAR4NPDV7rE7bA58fXClJNsBS6rq3IH27N+mNwNOSvJ4uqT5CICqOjfJu5I8iO738MmqumOO8WoRM0losZl6Ds0dtOHW9jC0zQfq3NKz7rBHza+Lnw1M30n3fxjgzKo6eIb1Qv+j7f+U7vlMj6Nr1+0Dy04DXkb3QMvfX8eYtUh5TUIL3dfoTo7QnSi/0qZXAU9s0wfSvROfzTeAfZMsbRe1DwbOnWWdm4BtR9j2ecA+SR4OkGSrJI8YrFBV1wM3tIe3QdeeKdsBV1X3COjfpXuA5ZT3A69p2/juCLFIv2CS0EL3J8BhSS6mO3m+upW/l+6Efz6wF/29h1+o7mmZrwfOBr4NXFhVsz1m+mLgjnax+U/7KrUhrEOBj7RYzwMeNaTqYcC72oXr2wbK/x44JMl5dENNv2hPVV1N9/WV/zRLrNI9+BRYaYFLshXdnVp7Tl0LkUZlT0JawJI8E/hP4O9MEFoX9iQkSb3sSUiSepkkJEm9TBKSpF4mCUlSL5OEJKnX/wdSxkutwbK8NgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_dw=data_departures[data_departures['day_num']<5]\n",
    "plt.hist(data_dw['hour'],bins=24);\n",
    "plt.xlabel('Hour of the day');\n",
    "plt.ylabel('Frequency');\n",
    "plt.title('Departure times for day of the week');\n",
    "plt.savefig('week.eps', format='eps', dpi=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHMlJREFUeJzt3Xm4HFWd//H3hwQIgbCZgGwhoBFERxavgINKGBbZCf5ARX9OYNCI4riMzzNExh8iwhDHBZxBxQBK2EWQzR0ji4KICSA7A0IIkJBc2RJ2E76/P+o0KS59b9e9SXUnfT6v57lPV53avlXVt759TlWfVkRgZmb5WqXTAZiZWWc5EZiZZc6JwMwsc04EZmaZcyIwM8ucE4GZWeacCKwjJP1S0qQOx7CLpPslPStpYg3rP1zSH5b3eits91OS5qf9ekOF+TsS51BJmi1pjzZs52xJJ9a9nRWBE0HN0pv2BUmLJD0t6UZJR0nqyLGXNE5SSBrexm0eL+m8cllE7BMR09sVQz9OAE6LiLUi4vIOx7JcSFoV+DawV9qvJ/pMb/v5txWfE0F7HBARo4DNganAMcBZ7Q5iefzzq9At75vNgbuGsuAKfCHdEBjBEPfL8tQt/9ArhYh4JiKuBD4ETJL0dgBJq0v6pqQ5qUp/uqQ10rQJkh6VdKykv6Uaxkcb65S0n6RbJS2U9Iik40vTGp/+jpQ0B/gdcH2a/HRqOnh330/sfT81SrpW0kmSbgCeB7aUtI6ksyTNk/SYpBMlDeu7z5L2Bo4FPpS295fSOj+ehg+XdIOkU1Kt6UFJ/5jKH5G0oNyM1OJ4jZb0s7SeJyX9vlnikvRXYEvgqhTX6pI2lnRlWu4BSZ8ozX+8pEsknSdpIXB4k3W+IS2/UNLNwJv6TP9O2p+FkmZJem8qf6Ok58vNOJLeKak3fcLvu53VJZ0qaW76OzWVvQW4r3R+f9d3WZqc/9J6vynpKUkPSdqnVF71XI9QUfsdnca/LGmxpLXT+ImSTi3tQ9NzmKbvL+k2La1Fv6PJviBp6xTvh9P4xpIuTcfuIUmfLc17vKSLJZ2jooZ+l6Se0vTtJd2Spv2YIqHmISL8V+MfMBvYo0n5HOBTafhU4EpgfWAUcBVwcpo2AVhMUd1fHdgVeA7YqjT9HyiS+juA+cDENG0cEMA5wJrAGqWy4aVYjgfOK42/Zh7g2hTv24DhwKrA5cAP0no3AG4GPtnPMXjN+kvr/HgaPjzt4xHAMODEtL3vpn3eC1gErFXheJ0MnJ5iXBV4L6Aq5wa4DvgexQVgO6AX2L20D38HJqZjvUaT9V0EXJyOyduBx4A/lKb/X+AN6Rh+EXgcGJGm/aLxfkjjpwD/00/cJwA3peM+BrgR+Fqzc9dk2Wbn//C0b59Ix/9TwNzGcRvkub4e+D9p+DfAX4F9StMOrnAOdwAWADuleCalc7V6+byl+eYA+6fyVYBZwHHAahSJ/kHg/aVz+CKwb1rvycBNadpqwMPAF9L75pB0TE7s9DWkLdepTgfQ7X99Lzal8puA/wBEcWF/U2nau4GH0vAEiovkmqXpFwP/r5/tnQqckoYb//RblqY3uxAcT+tEcEJp+obAS5QuhsBhwDX9xPSa9ZfWWU4E95em/UPa/oalsicoLs6tjtcJwBXAmwdzboDNgCXAqNL0k4GzS/tw/QDrGpYuHFuXyv6TUiJossxTwLZp+EPADaV1PQ7s2M9yfwX2LY2/H5jd3/nts2yz83848EBpfGSa541DONdfA/6bItk9DnyOojl0BPACMLrCOfw+KbGVpt8H7Fo6b18FHgV2K82zEzCnz3JfAn5UOoe/LU3bBnghDb+PUvJLZTeSSSJYUds5c7AJ8CTFJ7qRwCxJjWmiuBg0PBURz5XGHwY2BpC0E8U/2tspPtWsDvykz7YeWQ7xltexOcWnpnmlmFdZxu3MLw2/ABARfcvWovXx+gbFP/xv0vRpETG1wvY3Bp6MiEWlsoeBntL4QPs3huLiV57n4fIMkr4IfDxtK4C1KS6MUCSv0yVtCbwFeCYibh4g1vK6X30/LIPHGwMR8Xw6dmtRfGIfzLm+jqL2ugNwB3A1xf2wnSmSzd8kbcDA53BziqbTfy2tdzVeu49HAddFxDWlss2BjSU9XSobBvy+2X5SNHOOUNEEujHwWKQMkLzm/HUz3yPoAEnvokgEfwD+RnGRe1tErJv+1omItUqLrCdpzdL4WIpPLwAXUFSxN4uIdSiaRcRrRT/DDc9R/GM2vLHJPOXlHqH4lDi6FPPaEfG2ZvvbzzaHasDjFRGLIuKLEbElcADwb5J2r7DeucD6kkaVysZSNO80DLQfvRQ1t836LA9Auh9wDPBBYL2IWBd4hnSuIuJFipreR4GPAee2iHXzPtuZ28+8fQ32XAz2XN8IbAUcTHGhvjvFtx9FkoDW7/lHgJNK09aNiJERcWFpO0cBYyWd0ifWh/osNyoi9q2wn/OATVTKTJTOX7dzImgjSWtL2p+iLfm8iLgjIl4BzgBOSZ+UkLSJpPf3WfyrklZLF5T9WfqpfxTFJ9kXJe0IfKRFGL3AKxTtpw23Ae+TNFbSOhTV6X5FxDyK9t9vpX1aRdKbJO3azyLzgXHNbtoOVqvjlW4yvjn9Qy+kaO5ZUmG9j1BcxE5ONz3fARwJnF8xriXAT4HjJY2UtA1F23bDKIpE0QsMl3QcRY2g7ByKZpoDgfPo34XAlyWNSTdmj2sxf1mz89+vwZ7riHieop3+aJZe+G8EPtkYr/CePwM4StJOKqyp4qGIcpJeBOxN8b5t1PhuBhZKOkbSGpKGSXp7+uDVyh8pzs9nJQ2X9AFgxwrLdQUngva4StIiik8s/0FRdT6iNP0Y4AHgJhVPpPyW4lNVw+MU7clzKS5MR0XEvWnap4ET0vqPo/hU2a/0j3oScEN6ImPniLga+DFwO8U/8c8q7NM/U1TX706xXQJs1M+8jaT1hKRbKqy7lYGO1/g0/izFP/f3IuLaius9jKINfS5wGfCVdGyq+gxFc8rjwNnAj0rTfg38EvhfiiaHF+nTvBIRN1BcpG+JiNkDbOdEYCbF+boDuCWVtdTs/FdYbDDnGooL/qoUF+bG+CiWPrEEA5zDiJhJceP6tLS9B2jylFZEPA3sCewj6WspGR9AcS/pIYqax5nAOq12MCJeBj6QtvMUxT2bn7Zarls0ngqwFZSkCRS1h007HYvVT8UjnxdExJmdjsXy4ZvFZiuI1ISxA3BQp2OxvLhpyGwFIGk6RfPI5/s8uWRWu9qahiRtRdHu3LAlRRv2Oal8HMXzwB+MiKdqCcLMzFpqyz0CFV9Hf4ziCx9HUzzlMlXSFIpH6Y6pPQgzM2uqXYlgL4onMHaRdB8wISLmSdoIuDYithpo+dGjR8e4ceNqj9PMrJvMmjXrbxExptV87bpZ/GGKZ5+h6DZgHhTPKDeeI+5L0mRgMsDYsWOZOXNmWwI1M+sWkip9O7r2m8WSVqP4gkzfbg8GFBHTIqInInrGjGmZ0MzMbIja8dTQPhRfkGn0GzM/NQmRXhe0IQYzM+tHOxLBYSxtFoKiX5zGV+8nUXS2ZWZmHVJrIpA0kuIr4OWvak8F9pR0f5pWpWdIMzOrSa03i1O/Jm/oU/YEUKU3SDMzawN/s9jMLHNOBGZmmXMiMDPLnBOBmVnm3A212TIYN+Xng15m9tT9aojEbOhcIzAzy5wTgZlZ5tw0ZCuFdjXBDGU7Zis71wjMzDLnRGBmljknAjOzzDkRmJllzonAzCxzTgRmZplzIjAzy5wTgZlZ5pwIzMwy50RgZpY5JwIzs8y5ryHrWu43yKwa1wjMzDLnRGBmlrlaE4GkdSVdIuleSfdIerek9SVdLen+9LpenTGYmdnA6q4RfAf4VURsDWwL3ANMAWZExHhgRho3M7MOqS0RSFobeB9wFkBEvBwRTwMHAdPTbNOBiXXFYGZmrdVZI9gS6AV+JOlWSWdKWhPYMCLmAaTXDZotLGmypJmSZvb29tYYpplZ3upMBMOBHYDvR8T2wHMMohkoIqZFRE9E9IwZM6auGM3MsldnIngUeDQi/pTGL6FIDPMlbQSQXhfUGIOZmbVQWyKIiMeBRyRtlYp2B+4GrgQmpbJJwBV1xWBmZq3V/c3ifwXOl7Qa8CBwBEXyuVjSkcAc4NCaYzAzswHUmggi4jagp8mk3evcrpmZVedvFpuZZc6JwMwsc04EZmaZcyIwM8ucE4GZWeacCMzMMudEYGaWOScCM7PMORGYmWXOicDMLHNOBGZmmXMiMDPLnBOBmVnm6u6G2sz6GDfl54NeZvbU/WqIxKzgGoGZWeacCMzMMuemIWu7oTSNmFl9XCMwM8ucE4GZWeacCMzMMudEYGaWOScCM7PM1frUkKTZwCJgCbA4InokrQ/8GBgHzAY+GBFP1RmHmZn1rx01gt0iYruI6EnjU4AZETEemJHGzcysQzrRNHQQMD0NTwcmdiAGMzNL6k4EAfxG0ixJk1PZhhExDyC9btBsQUmTJc2UNLO3t7fmMM3M8lX3N4t3iYi5kjYArpZ0b9UFI2IaMA2gp6cn6grQzCx3tdYIImJuel0AXAbsCMyXtBFAel1QZwxmZjaw2hKBpDUljWoMA3sBdwJXApPSbJOAK+qKwczMWquzaWhD4DJJje1cEBG/kvRn4GJJRwJzgENrjMHMzFqoLRFExIPAtk3KnwB2r2u7ZmY2OP5msZlZ5pwIzMwy50RgZpY5JwIzs8w5EZiZZc6JwMwsc04EZmaZcyIwM8ucE4GZWeacCMzMMld3N9SWgXFTft7pEMxsGbhGYGaWOScCM7PMORGYmWXOicDMLHNOBGZmmauUCCS9ve5AzMysM6rWCE6XdLOkT0tat9aIzMysrSolgoh4D/BRYDNgpqQLJO1Za2RmZtYWle8RRMT9wJeBY4Bdgf+WdK+kD9QVnJmZ1a/qPYJ3SDoFuAf4J+CAiHhrGj6lxvjMzKxmVbuYOA04Azg2Il5oFEbEXElfriUyMzNri6qJYF/ghYhYAiBpFWBERDwfEecOtKCkYcBM4LGI2F/SFsBFwPrALcDHIuLlIe+BWQYG25/T7Kn71RSJdaOq9wh+C6xRGh+Zyqr4HEWTUsPXgVMiYjzwFHBkxfWYmVkNqiaCERHxbGMkDY9stZCkTYH9gDPTuCjuK1ySZpkOTBxMwGZmtnxVTQTPSdqhMSLpncALA8zfcCrw78ArafwNwNMRsTiNPwps0mxBSZMlzZQ0s7e3t2KYZmY2WFXvEXwe+ImkuWl8I+BDAy0gaX9gQUTMkjShUdxk1mi2fERMA6YB9PT0NJ3HzMyWXaVEEBF/lrQ1sBXFxfzeiPh7i8V2AQ6UtC8wAlibooawrqThqVawKTB3gHWYmVnNBtPp3LuAdwDbA4dJ+ueBZo6IL0XEphExDvgw8LuI+ChwDXBImm0ScMWgozYzs+WmUo1A0rnAm4DbgCWpOIBzhrDNY4CLJJ0I3AqcNYR1mJnZclL1HkEPsE1EDKmtPiKuBa5Nww8COw5lPWZmtvxVbRq6E3hjnYGYmVlnVK0RjAbulnQz8FKjMCIOrCUqMzNrm6qJ4Pg6gzAzs86p+vjodZI2B8ZHxG8ljQSG1RuamZm1Q9VuqD9B0S3ED1LRJsDldQVlZmbtU/Vm8dEUXxBbCK/+SM0GdQVlZmbtUzURvFTuKlrScPrpGsLMzFYuVRPBdZKOBdZIv1X8E+Cq+sIyM7N2qZoIpgC9wB3AJ4FfUPx+sZmZreSqPjX0CsVPVZ5RbzhmZtZuVfsaeogm9wQiYsvlHpGZLbPB/rQl+OctczaYvoYaRgCHUvzmsJmZreQq3SOIiCdKf49FxKkUPzlpZmYruapNQzuURlehqCGMqiUiMzNrq6pNQ98qDS8GZgMfXO7RmJlZ21V9ami3ugMxM7POqNo09G8DTY+Iby+fcMzMrN0G89TQu4Ar0/gBwPXAI3UEZWZm7TOYH6bZISIWAUg6HvhJRHy8rsDMzKw9qnYxMRZ4uTT+MjBuuUdjZmZtV7VGcC5ws6TLKL5hfDBwTm1RmZlZ21R9augkSb8E3puKjoiIW+sLy8zM2qVq0xDASGBhRHwHeFTSFgPNLGmEpJsl/UXSXZK+msq3kPQnSfdL+rGk1ZYhfjMzW0ZVf6ryK8AxwJdS0arAeS0Wewn4p4jYFtgO2FvSzsDXgVMiYjzwFHDkUAI3M7Plo2qN4GDgQOA5gIiYS4suJqLwbBpdNf0FRR9Fl6Ty6cDEQcZsZmbLUdWbxS9HREgKAElrVllI0jBgFvBm4LvAX4GnI2JxmuVRYJN+lp0MTAYYO3ZsxTDNbKjcdXW+qtYILpb0A2BdSZ8AfkuFH6mJiCURsR2wKbAj8NZms/Wz7LSI6ImInjFjxlQM08zMBqvqU0PfTL9VvBDYCjguIq6uupGIeFrStcDOFMlkeKoVbArMHXzYZma2vLRMBKl559cRsQdQ+eIvaQzw95QE1gD2oLhRfA1wCHARMAm4YiiBm1ke3GRVv5aJICKWSHpe0joR8cwg1r0RMD0lklWAiyPiZ5LuBi6SdCJwK3DWkCI3M7PlourN4heBOyRdTXpyCCAiPtvfAhFxO7B9k/IHKe4XmJnZCqBqIvh5+jMzsy4zYCKQNDYi5kTE9HYFZGZm7dXq8dHLGwOSLq05FjMz64BWiUCl4S3rDMTMzDqjVSKIfobNzKxLtLpZvK2khRQ1gzXSMGk8ImLtWqMzM7PaDZgIImJYuwIxM7POGMzvEZiZWRdyIjAzy5wTgZlZ5qp+s9gyMZQOvsxs5eYagZlZ5pwIzMwy56YhM+s6g23izP33C1wjMDPLnBOBmVnmnAjMzDLnRGBmljknAjOzzPmpITNrmxX1C4tDiaubnjRyjcDMLHNOBGZmmastEUjaTNI1ku6RdJekz6Xy9SVdLen+9LpeXTGYmVlrddYIFgNfjIi3AjsDR0vaBpgCzIiI8cCMNG5mZh1SWyKIiHkRcUsaXgTcA2wCHARMT7NNBybWFYOZmbXWlnsEksYB2wN/AjaMiHlQJAtgg36WmSxppqSZvb297QjTzCxLtScCSWsBlwKfj4iFVZeLiGkR0RMRPWPGjKkvQDOzzNWaCCStSpEEzo+In6bi+ZI2StM3AhbUGYOZmQ2szqeGBJwF3BMR3y5NuhKYlIYnAVfUFYOZmbVW5zeLdwE+Btwh6bZUdiwwFbhY0pHAHODQGmMwM7MWaksEEfEHQP1M3r2u7ZqZ2eC4ryEzG7IVte8gGxx3MWFmljknAjOzzDkRmJllzonAzCxzTgRmZplzIjAzy5wTgZlZ5pwIzMwy50RgZpY5JwIzs8w5EZiZZc6JwMwsc04EZmaZcyIwM8ucE4GZWeacCMzMMudEYGaWOScCM7PMORGYmWXOicDMLHNOBGZmmastEUj6oaQFku4sla0v6WpJ96fX9eravpmZVVNnjeBsYO8+ZVOAGRExHpiRxs3MrINqSwQRcT3wZJ/ig4DpaXg6MLGu7ZuZWTXtvkewYUTMA0ivG/Q3o6TJkmZKmtnb29u2AM3McrPC3iyOiGkR0RMRPWPGjOl0OGZmXavdiWC+pI0A0uuCNm/fzMz6aHciuBKYlIYnAVe0eftmZtZHnY+PXgj8EdhK0qOSjgSmAntKuh/YM42bmVkHDa9rxRFxWD+Tdq9rm2ZmNngr7M1iMzNrDycCM7PMORGYmWXOicDMLHNOBGZmmXMiMDPLXG2Pj5qZdbNxU34+6GVmT92vhkiWnWsEZmaZcyIwM8tc1zcNdVP1bbCGsu9mVp/B/k+261rkGoGZWeacCMzMMudEYGaWOScCM7PMORGYmWXOicDMLHNd//hoN/HjoGZWB9cIzMwy50RgZpY5Nw01saJ++8/MrA6uEZiZZc6JwMwscx1pGpK0N/AdYBhwZkRM7UQcy4uf5jGzlVnbawSShgHfBfYBtgEOk7RNu+MwM7NCJ5qGdgQeiIgHI+Jl4CLgoA7EYWZmdKZpaBPgkdL4o8BOfWeSNBmYnEaflXTfELc3GvjbEJftJj4OS/lYFHwcllohj4W+vsyr2LzKTJ1IBGpSFq8riJgGTFvmjUkzI6JnWdezsvNxWMrHouDjsFTux6ITTUOPApuVxjcF5nYgDjMzozOJ4M/AeElbSFoN+DBwZQfiMDMzOtA0FBGLJX0G+DXF46M/jIi7atzkMjcvdQkfh6V8LAo+DktlfSwU8brmeTMzy4i/WWxmljknAjOzzHVtIpC0t6T7JD0gaUqn4+kkSbMl3SHpNkkzOx1PO0n6oaQFku4sla0v6WpJ96fX9ToZYzv0cxyOl/RYel/cJmnfTsbYDpI2k3SNpHsk3SXpc6k8u/dEWVcmAndj0dRuEbFdhs9Knw3s3adsCjAjIsYDM9J4tzub1x8HgFPS+2K7iPhFm2PqhMXAFyPircDOwNHp2pDje+JVXZkIcDcWlkTE9cCTfYoPAqan4enAxLYG1QH9HIfsRMS8iLglDS8C7qHo7SC790RZtyaCZt1YbNKhWFYEAfxG0qzUdUfuNoyIeVBcGIANOhxPJ31G0u2p6Sir5hBJ44DtgT+R+XuiWxNBpW4sMrJLROxA0VR2tKT3dTogWyF8H3gTsB0wD/hWZ8NpH0lrAZcCn4+IhZ2Op9O6NRG4G4uSiJibXhcAl1E0neVsvqSNANLrgg7H0xERMT8ilkTEK8AZZPK+kLQqRRI4PyJ+moqzfk90ayJwNxaJpDUljWoMA3sBdw68VNe7EpiUhicBV3Qwlo5pXPiSg8ngfSFJwFnAPRHx7dKkrN8TXfvN4vQo3Kks7cbipA6H1BGStqSoBUDRpcgFOR0LSRcCEyi6GZ4PfAW4HLgYGAvMAQ6NiK6+kdrPcZhA0SwUwGzgk4128m4l6T3A74E7gFdS8bEU9wmyek+UdW0iMDOzarq1acjMzCpyIjAzy5wTgZlZ5pwIzMwy50RgZpY5JwJbKUh6ts/44ZJOa3MMh6ZeK6/pUz5O0kfqiC2tu+uf77fOciKwrKWeaqs6Evh0ROzWp3wc8JHXz262cnAisJWepM0lzUidp82QNDaVny3pkNJ8z6bXCalP+gsovljUd32Hpd9vuFPS11PZccB7gNMlfaPPIlOB96Y+/b+QyjaW9KvUv/1/lda9l6Q/SrpF0k9Snzd9t/9OSX+R9Efg6FL5OEm/T8veIukfU/m5kg4qzXe+pAMHeRgtZxHhP/+t8H/AEuC20t8c4LQ07SpgUhr+F+DyNHw2cEhpHc+m1wnAc8AWTbazcVr3GIpvYv8OmJimXQv0NFlmAvCz0vjhwIPAOsAI4GGKvq9GA9cDa6b5jgGOa7K+24Fd0/A3gDvT8EhgRBoeD8xMw7uW9nkd4CFgeKfPmf9Wnr/hQ08hZm31QkRs1xiRdDjQ+JGddwMfSMPnAv9FazdHxENNyt8FXBsRvWk75wPvo+iWYjBmRMQzaR13A5sD61L8UNINRZc3rAb8sbyQpHWAdSPiutL+7JOGVwVOk7QdRWJ8C0BEXCfpu5I2oDgOl0bE4kHGaxlzIrBu1Og3ZTGp+TN1NrZaaZ7n+lm2WRfmQ/FSaXgJxf+agKsj4rABlhP9d5n+BYp+gral2K8XS9POBT5K0cHivwwxZsuU7xFYN7iR4gIIxcXwD2l4NvDONHwQxSfqVv4E7CppdLqRfBhwXYtlFgGjKqz7JmAXSW8GkDRS0lvKM0TE08AzqXM0KPanYR1gXhTdRn+MokPFhrOBz6d13FUhFrNXORFYN/gscISk2ykukJ9L5WdQXNRvBnai/1rAq6LoffNLwDXAX4BbIqJVl8S3A4vTDd4v9DdTam46HLgwxXoTsHWTWY8AvptuFr9QKv8eMEnSTRTNQq/uT0TMp/jZxR+1iNXsddz7qFkXkDSS4gmoHRr3Jsyqco3AbCUnaQ/gXuB/nARsKFwjMDPLnGsEZmaZcyIwM8ucE4GZWeacCMzMMudEYGaWuf8PRL2wWo+dneoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_dnw=data_departures[data_departures['day_num']>=5]\n",
    "plt.hist(data_dnw['hour'],bins=24);\n",
    "plt.xlabel('Hour of the day');\n",
    "plt.ylabel('Frequency');\n",
    "plt.title('Departure times for day of the weekend');\n",
    "plt.savefig('weekend.eps', format='eps', dpi=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figures formatted for paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XuYHGWZ9/Hvz8mJg0IM0SUzExME0XCQwxAQEZUIBkXialgCrgmH3ez6GkHR5QUPiKi7IojiC7tuhHDGgFlgI0QCAuLqAmaCEQjhEENIhqDkgGCUSBLu94+qgUpPz0z3THf1dM/vc119TXfVU1139czcc89TTz2liMDMzMzMzBKvq3UAZmZmZmYDiQtkMzMzM7MMF8hmZmZmZhkukM3MzMzMMlwgm5mZmZlluEA2MzMzM8vItUCWNFnS45KWSzqryPrDJT0oaYukqZnl+0m6T9JSSQ9JOj7PuM3MGkFfc3C6boakJ9PHjPyiNjPLn/KaB1lSE/AEcCTQASwCToiIRzNtxgFvAL4AzI+IeenytwEREU9KGgMsBt4REX/MJXgzszrXzxz8RqAdaAOCJAcfGBHP53gIZma5ybMHeSKwPCJWRMTLwFxgSrZBRKyMiIeAVwqWPxERT6bP1wDPAaPzCdvMrCH0OQcDHwTujIgNaVF8JzA5j6DNzGphSI77agZWZ153AAeX+yaSJgLDgN8VWTcTmAmwww47HPj2t7+9b5GamfXT4sWL10XEQPpHvj85uNi2zcUaOg+b2UDQ3xycZ4GsIsvKGt8haVfgGmBGRBT2cBARs4HZAG1tbdHe3t6XOM3M+k3S07WOoUB/cnDJ2zoPm9lA0N8cnOcQiw6gNfO6BVhT6saS3gDcBnw5Iu6vcGxmZo2uPzm4X/nbzKze5FkgLwL2kDRe0jBgGjC/lA3T9jcDV0fEj6sYo5lZo+pzDgYWAkdJGilpJHBUuszMrCHlViBHxBZgFklSXQbcGBFLJZ0n6VgASQdJ6gCOA/5T0tJ0878DDgdOkrQkfeyXV+xmZvWuPzk4IjYAXycpshcB56XLzMwaUm7TvOXNY9+sGppbx7KmY3WPbca0tPLM6lU5RTTwbd68mY6ODjZt2lTrUKpixIgRtLS0MHTo0G2WS1ocEW01CmtAcB42qz3n4L7J8yI9s7q3pmM1My++tcc2s08/Jqdo6kNHRwevf/3rGTduHFKxa73qV0Swfv16Ojo6GD9+fK3DMTPrwjm4b3yraTOrqk2bNjFq1KiGS8wAkhg1alTD9syYWf1zDu4bF8hmVnWNmJg7NfKxmVljaOQ8Va1jc4FsRjK2WFKvDzMzM2t8HoNsRmlji8HjiyuhlAsdy1HKRZFNTU3ss88+RARNTU1ccsklHHrooaxZs4bTTjuNefPmceWVV9Le3s4ll1xSsdjMzAYa5+DSuEA2s1yV+s9IqUr5p2W77bZjyZIlACxcuJCzzz6be++9lzFjxjBv3ryKxWJmNtA5B5fGQyzMKqxpyNCShms0t46tdaiD0osvvsjIkSMBWLlyJXvvvXeXNrfddhvvete7WLduHU8//TSTJk1i3333ZdKkSaxa5Sn8zMz6ql5ysHuQzSps65bNHq4xwLz00kvst99+bNq0iWeffZa7776727Y333wzF110EQsWLGDkyJGcfPLJTJ8+nRkzZjBnzhxOO+00brnllhyjNzOrb/WYg10gm1nDy57eu++++5g+fTqPPPJIl3b33HMP7e3t3HHHHbzhDW94tf1NN90EwCc/+UnOPPPM/AI3M2sA9ZiDPcTCzAaVztN2a9eu7bJut912409/+hNPPPFEt9t7NhMzs76rlxzsAtnMBpXHHnuMrVu3MmrUqC7r3vKWt3DTTTcxffp0li5dCsChhx7K3LlzAbjuuus47LDDco3XzKyR1EsO9hALM8vVmJbWio6/HtPS2mubzvFvkNya9KqrrqKpqalo2z333JPrrruO4447jp/85Cd8//vf55RTTuGCCy5g9OjRXHHFFRWL3cwsb87BpXGBbGa56m2+zGrYunVr0eXjxo17dRzcSSedxEknnQTA/vvvz6OPPvpqu54uKDEzqyfOwaXxEAszMzMzswwXyGZmVhWLFy/2vN9mVpc8xMLMqi4iGnb2h4iodQgD1i6tu/OxL3wP8LzfZrXkHFw+9yCbWVWNGDGC9evXN2QhGRGsX7+eESNG1DoUM7OinIP7xj3IZlZVLS0tdHR0FJ3zshGMGDGClpaWWodhZlaUc3Df5FogS5oMXAw0AZdFxLcK1h8OfA/YF5gWEfMy62YAX05ffiMirsonajPrj6FDhzJ+/Phah2FmNig5B/dNbkMsJDUBlwJHAxOAEyRNKGi2CjgJuL5g2zcCXwUOBiYCX5U0stoxm5mZmdngk+cY5InA8ohYEREvA3OBKdkGEbEyIh4CXinY9oPAnRGxISKeB+4EJucRtJmZmZkNLnkWyM3A6szrjnRZxbaVNFNSu6T2Rh1rY2ZmZmbVlWeBXGx+kVIvqSxp24iYHRFtEdE2evTosoIzM7P60dw69tU5lrMPz7dsZpWQ50V6HUD2ht0twJoytn1fwbY/r0hUZmaDRAkXSg8HrgYOBNYDx0fESklDgcuAA0j+blwdEf+Wa/AF1nSsZubFt3ZZ7vmWzawS8uxBXgTsIWm8pGHANGB+idsuBI6SNDK9OO+odJlZj7rrZSp8mDW6Ei+UPhV4PiJ2B74LnJ8uPw4YHhH7kBTP/yRpXB5xm5nVQm49yBGxRdIsksK2CZgTEUslnQe0R8R8SQcBNwMjgY9I+lpE7BURGyR9naTIBjgvIjbkFbvVr+56mQq518kGgVcvlAaQ1Hmh9KOZNlOAc9Pn84BLlPwHGcAOkoYA2wEvAy/mFLeZWe5ynQc5IhYACwqWnZN5vohk+ESxbecAc6oaoJlZ4yp2sfPB3bVJOzVeAEaRFMtTgGeB7YHPuZPCzBqZbzVtZjY4lHKxc3dtJgJbgTHAeODzknYrupPMbEKbNr7Qn3jNzGrGBbKZ2eBQyoXSr7ZJh1PsBGwATgRuj4jNEfEc8CugrdhOsrMJjdhxpwofQv945gszK1WuQyzMzKxmXr1QGniG5ELpEwvazAdmAPcBU4G7IyIkrQKOkHQtyRCLQ4DvVSKo5taxrOlYXXTdmJZWnlm9qhK7ATzzhZmVzgWymdkgUMqF0sDlwDWSlpP0HE9LN78UuAJ4hGQYxhXpXU/7racLaftSuDYNGeqZacys31wgm5kNEiVcKL2JZEq3wu02Fls+EG3dsrmiBbeZDU4eg2xmZtYNj1s2G5zcg2xmZgNSd8MlxrS0FmldHR63bDY4uUA2M7MBqbvhEi5OzazaPMTCrEY6e8d6evg0rpmZWf7cg2xWIz1dTNTJPWVmZmb5cw+ymZmZmVmGC2QzMzMzswwXyGZmZmZmGS6QzczMzMwyXCBb3epuAv/sw8wGhu5mbRkIeppRxswGJ89iYXWruwn8szwLhNnAMJDnNPbtqc2skHuQzczMzMwyXCCbmZmZmWXkWiBLmizpcUnLJZ1VZP1wSTek6x+QNC5dPlTSVZIelrRM0tl5xm1mZmZmg0duBbKkJuBS4GhgAnCCpAkFzU4Fno+I3YHvAueny48DhkfEPsCBwD91Fs9mZmZmZpWUZw/yRGB5RKyIiJeBucCUgjZTgKvS5/OASUouIw5gB0lDgO2Al4EX8wnbzMzMzAaTPAvkZmB15nVHuqxom4jYArwAjCIplv8MPAusAi6MiA2FO5A0U1K7pPa1a9dW/gjMzMzMrOHlWSAXm1AySmwzEdgKjAHGA5+XtFuXhhGzI6ItItpGjx7d33jNzMzMbBDKs0DuAFozr1uANd21SYdT7ARsAE4Ebo+IzRHxHPAroK3qEZuZmZnZoFNWgZwWrX21CNhD0nhJw4BpwPyCNvOBGenzqcDdEREkwyqOUGIH4BDgsX7EYmZWd/qZg83MrETl9iA/K+lCSe8od0fpmOJZwEJgGXBjRCyVdJ6kY9NmlwOjJC0HzgA6p4K7FNgReISk0L4iIh4qNwYzszrX5xxsZmalK7c34ovAycDnJP0auAy4ISI2lrJxRCwAFhQsOyfzfBPJlG6F220sttzMbJDpVw42M7PSlNWDHBE/jIhDgb2BXwLfIOnRmCPp3dUI0MzMEs7BZmb56NNFehGxLCL+heRCuy+SXET3C0mPSfpnSb6FtZlZlTgHm5lVV5+SqKRhkqYBPyW54939wEnAFcBXgOsrFaCZmW2rrzlY0mRJj0taLumsIuuHS7ohXf9A9o6lkvaVdJ+kpZIeljSi8kdmZjYwlDUGWdIBwCnACcBm4Grg0xHxZKbNXcD/VDJIMzPrXw6W1ERywfORJFNqLpI0PyIezTQ7FXg+InZPC/DzgePT2TOuBT4ZEb+VNCrdv5lZQyr3Ir1fA3cCM4H/TmemKLSM5DbSZmZWWf3JwROB5RGxAkDSXGAKkC2QpwDnps/nAZdIEnAU8FBE/BYgItb3/1DMzAaucgvkt0bE0z01iIg/k1xlbWZmldWfHNwMrM687gAO7q5NRGyR9AIwCngbEJIWAqOBuRHx7WL7lzSTpIBnx5GNe0fTpiFDSf536GpMSyvPrF6Vc0RmVknlFsj3SDqosPdA0s7AgxHR5fbPZmZWMf3JwcWquSixzRDgMOAg4C/AXZIWR8RdXRpHzAZmA4weu0fh+zeMrVs2M/PiW4uum336MTlHY2aVVu5FeuOApiLLh5P0PJiZWfWMo+85uANozbxuAdZ01yYdd7wTsCFdfm9ErIuIv5DMZ39AucGbmdWLknqQJX0s8/LD6Wm3Tk3AJGBlBeMyM7NUhXLwImAPSeOBZ4BpJNPDZc0HZgD3AVOBuyOic2jFmZK2B14G3ksye4aZWUMqdYjFvPRrkNwOOmszSWL+fIViMjOzbfU7B6djimcBC0mK6jkRsVTSeUB7RMxP3/saSctJeo6npds+L+kikiI7gAURcVtFjszMbAAqqUCOiNcBSHoKOCgi1lU1KjMze1WlcnBELCAZHpFddk7m+SbguG62vZZkqjczs4ZX1kV6ETG+WoGYmVnPnIPNzPLRa4Es6Qzg3yNiU/q8WxFxUcUiMzMz52AzsxoopQf5M8BVwKb0eXcCcHI2M6ss5+A6090cyZ4f2ax+9FogZ0/p+fSemVm+nIPrT3dzJHt+ZLP6Ue48yF1IGlqJQMzMrHzOwWZmlVdWgSzpNEkfz7yeA7wk6XFJe1Y8OjMze5VzsJlZPsrtQT4NWAsg6XCS6YBOBJYA36lsaGZmVsA52MwsB+UWyM28dremjwA/jogbgXOBQ3rbWNLktKdjuaSziqwfLumGdP0DksZl1u0r6T5JSyU9LGlEmbGb1Z3Oi316ezS3jq11qJaPfuVgMzMrTVnzIAMvAqOBVcCRwAXp8s1AjwWrpCbg0nS7DmCRpPkR8Wim2anA8xGxu6RpwPnA8ZKGkExQ/8mI+K2kUek+zRpadxf7FPLFP4NGn3OwmZmVrtwe5DuAH0q6HNgd+Gm6fC/gqV62nQgsj4gVEfEyMBeYUtBmCsl0RpDcWnWSkrlyjgIeiojfAkTE+ojYWmbsZmb1rj852MzMSlRugfxp4FfALsDUiNiQLj8A+FEv2zYDqzOvO9JlRdtExBbgBWAU8DYgJC2U9KCkM4vtQNJMSe2S2teuXVvGYZmZ1YX+5GAzMytRubeafpEiE9VHxFdL2LzrrOnJxPaltBkCHAYcBPwFuEvS4oi4qyCO2cBsgLa2tsL3NjOra/3MwWZmVqJyxyADIGkM8CYKeqAj4sEeNusAWjOvW4A13bTpSMcd7wRsSJffGxHr0v0vIOkxuQszs0GmjznYzMxKVFaBLGl/kovl3k7X3t4AmnrYfBGwh6TxwDPANJLpibLmAzOA+4CpwN0REZIWAmdK2h54GXgv8N1yYjczq3f9zMFmZlaicnuQZ5OMEf5Hkt7fkocxRMQWSbOAhSRJfE5ELJV0HtAeEfOBy4FrJC0n6Tmelm77vKSLSIrsABZExG1lxm51orl1LGs6Vvfe0Gzw6XMONjOz0pVbIE8A9o+IJ/qys4hYACwoWHZO5vkmkonvi217LUnPiTW4NR2rPbWZWXH9ysFmZlaacmexeBj4m2oEYmZmvXIONjPLQbkF8heBb0v6gKQ3S3pj9lGNAM3M7FXOwWZmOSh3iMXP0q93sO3YN+ELRMzMqs052MwsB+UWyO+vShRmZlYK52AzsxyUe6OQe6sViJmZ9cw52MwsH+WOQUbSPpIukfRTSbumyz6azs9pZmZV5Bxcv5qGDEVS0Udz69hah2dmGeXeKOQokpt5/BQ4AtguXfVW4CTgo5UMzszMXuMcXN+2btnc7RSWnrbSbGAptwf568AZEfG3JHe06/RzYGKlgjIzs6L6lYMlTZb0uKTlks4qsn64pBvS9Q9IGlewfqykjZK+0K+jMDMb4MotkPei4EYfqQ2ApxgyM6uuPudgSU3ApcDRJDccOUHShIJmpwLPR8TuwHeB8wvWf5ek99rMrKGVWyA/DzQXWX4A0NH/cMzMrAf9ycETgeURsSIiXgbmAlMK2kwBrkqfzwMmSRIk45yBFcDSPsZuZlY3yi2QrwcukNRCMufmEEnvBS4Erq50cGZmto3+5OBmYHXmdQddi+1X20TEFuAFYJSkHYD/C3yttwAlzZTULql908YXSjgkM7OBp9wC+cvAU8DTwI7Ao8A9wC+Bb1Y2NDMzK9CfHKwiy6LENl8DvhsRG3sLMCJmR0RbRLSN2HGn3pqbmQ1I5c6DvBn4hKSvkJzSex3wm4h4shrBmZnZa/qZgzuA1szrFmBNN206JA0BdiIZ33wwMFXSt4GdgVckbYqIS/p1QGZmA1TJBbKk7YAzgY8Du5H0KqwAfizpOxHxUnVCNDOzCuTgRcAeksYDzwDTgBML2swHZgD3AVOBuyMigPdk4jgX2Oji2MwaWUkFctqTcDdJj8XtwG0kp+ImAOcAR0t6bzpmzczMKqgSOTgitkiaBSwEmoA5EbFU0nlAe0TMBy4HrpG0nKTneFo1j8vMbKAqtQd5JrA7cEBEbHMFs6S9ScbAzQT+vbLhmZkZFcrBEbGAgmniIuKczPNNwHG9vMe55QRuZlaPSr1IbyrwzcLEDBARjwD/lrYxM7PKcw42M8tRqQXyXiSn97rzM2Dv/odjZmZFOAc3uKYhQ5HU5dHcOrbWoZkNSqUOsRgJrO1h/VqSK5t7JGkycDHJ+LfLIuJbBeuHk8zleSCwHjg+IlZm1o8lmdbo3Ii4sMTYzczqXUVysA1cW7dsZubFt3ZZPvv0Y2oQjZmV2oPcBPR0Ad4raZtu+Tan1tw6tmgPSeHDzLrodw42M7PSldqDLOBaSX/tZv3wEt7j1ducAkjqvM3po5k2U4Bz0+fzgEskKSIic5vTP5cYsw0wazpWF+0hKeQek/J1np7tzZiWVp5ZvSqHiKzCKpGDzcysRKUWyFeV0KYvtzk9uLs26ZREnbc5fYnkNqdHAl/obgeSZpJcyc3YsR63ZYNHd6dnC/mfj7pViRxsZmYlKqlAjoiTK7CvitzmtKdesoiYDcwGaGtrK3xvM7O6VKEcbGZmJSrrVtP95NucmpmZmdmAl2eB7NucmpmZmdmAl1uB7NucmpmZmVk9yLMH2bc5NTMzM7MBr9R5kM3MzMzMBgUXyGZmZmZmGS6QzczMzMwyXCCbmZmZmWW4QDYzMzMzy3CBbGZmNoA1t45FUpdHc+vYWodm1rBynebNzMzMyrOmYzUzL761y/LZpx9Tg2jMBgf3IJuZmZmZZbhANjMzMzPLcIFsZmZmZpbhAtnMzMzMLMMFspmZmZlZhgtkq4jupiHKPsystiRNlvS4pOWSziqyfrikG9L1D0galy4/UtJiSQ+nX4/IO3Yzszx5mjeriO6mIcrylERmtSOpCbgUOBLoABZJmh8Rj2aanQo8HxG7S5oGnA8cD6wDPhIRayTtDSwEmvM9AivUNGRot50PY1paeWb1qpwjMmscLpDNzAaHicDyiFgBIGkuMAXIFshTgHPT5/OASyQpIn6TabMUGCFpeET8tfphW3e2btncbceEOyTM+sdDLMzMBodmYHXmdQdde4FfbRMRW4AXgFEFbT4O/Ka74ljSTEntkto3bXyhIoGbmeXNBbKZ2eBQ7Fx8lNNG0l4kwy7+qbudRMTsiGiLiLYRO+7Up0DNzGot1wLZF4iY1VbnmMWeHs2tY2sdplVHB9Caed0CrOmujaQhwE7AhvR1C3AzMD0iflf1aM3Maii3Mci+QKQ+NbeOZU3H6t4bWl3oacxiJ49dbFiLgD0kjQeeAaYBJxa0mQ/MAO4DpgJ3R0RI2hm4DTg7In6VY8xmZjWR50V6vkCkDpUyOwW4qDIb6CJii6RZJB0MTcCciFgq6TygPSLmA5cD10haTtJzPC3dfBawO/AVSV9Jlx0VEc/lexRmZvnIs0AudoHIwd21SZN55wUi6zJterxApF6V0lPraXvMrD8iYgGwoGDZOZnnm4Djimz3DeAbVQ/QzGyAyLNAruQFIkcV3YE0E5gJMHZsfY2j9DzCZmZmZgNDnhfpVf0CkezV06NHj65w+GZmZmY2GORZIL96gYikYSRj2+YXtOm8QAR8gYiZmZmZ1UBuBXI66XznBSLLgBs7LxCRdGza7HJgVHqByBlA51Rw2QtElqSPN+UVu5mZWT3pbkpHT+NoVppcbzXtC0QGFk/hZmbWmLqb0tHXspiVJtcC2QYWXxhoZmZm1pVvNV1HSrkLmk+hmZmZmfWPe5CrrJLDGEq5Cxq419fMzMysP1wgV5nvRGdmZgNF55nIYnwzKrPXuEBuQD0lQDMzG7x6OhPpjhqz17hAbkAeimFmZmbWd75Iz8zMzMwswwWymZmZmVmGC2QzMzMzswwXyGZmZmZmGS6QzWwbviGNmZkNdp7Fwsy24VlQzAan7qYI9fzINhi5QO6HSt4lz8zMrJa6++fY/wzbYOQCuR9KuUueE4uZmTWq7jqK3Ots9c4FspmZmXWrt7uzutfZGpELZDPrk1Jvae6eJLP65ttT22DkAtnM+sQX85mZWaMaVNO8NbeO9fRVZmZmOejub67/xlo9yLUHWdJk4GKgCbgsIr5VsH44cDVwILAeOD4iVqbrzgZOBbYCp0XEwnL3X8pFdQCXf/5vSzp1bGa9K2UoxtDhI9j81029vpeHa/RPrXOwDS7d/c3t6axSXy7684WCVg25FciSmoBLgSOBDmCRpPkR8Wim2anA8xGxu6RpwPnA8ZImANOAvYAxwM8kvS0itlYjVp86NqucUn6fZp9+jH/nqqyecrANXn0pqvuyjVlv8hxiMRFYHhErIuJlYC4wpaDNFOCq9Pk8YJKSrqcpwNyI+GtEPAUsT9+vW4sXL+5yWsfMbBDLNQebdae7u3X2NPSipzt89mWbYSO28/AP65EiIp8dSVOByRHxD+nrTwIHR8SsTJtH0jYd6evfAQcD5wL3R8S16fLLgZ9GxLyCfcwEZqYv9wYeqepBVd8uwLpaB9EPjr+2HH9t7RkRr691EJ3yyMHpuoGUh2v9M+T9e//ef+30KwfnOQa52L95hdV5d21K2ZaImA3MBpDUHhFt5QY5kNT7MTj+2nL8tSWpvdYxFKh6DoaBlYe9f+/f+x/c++/P9nkOsegAWjOvW4A13bWRNATYCdhQ4rZmZtY952AzsxLlWSAvAvaQNF7SMJILPuYXtJkPzEifTwXujmQMyHxgmqThksYDewC/ziluM7NG4BxsZlai3IZYRMQWSbOAhSRTDM2JiKWSzgPaI2I+cDlwjaTlJL0W09Jtl0q6EXgU2AJ8uoSrp2dX61hyVO/H4Phry/HX1oCKvwY5GGr/GXj/3r/37/33SW4X6ZmZmZmZ1YNBdSc9MzMzM7PeuEA2MzMzM8toyAJZ0mRJj0taLumsWsdTLkkrJT0sackAnCqqC0lzJD2XzqHaueyNku6U9GT6dWQtY+xNN8dwrqRn0u/DEkkfqmWM3ZHUKukeScskLZV0erq8Lr4HPcRfF58/gKQRkn4t6bfpMXwtXT5e0gPp9+CG9OK4hpd3Dh4ovwOSmiT9RtKt6evcvv+SdpY0T9Jj6efwrjyPX9Ln0s/+EUk/Sn8nqnr85fztUeL76c/kQ5IOqNL+L0i/Bw9JulnSzpl1Z6f7f1zSB6ux/8y6L0gKSbukr3M5/nT5Z9JjXCrp25nlVT9+SftJuj/9m9EuaWK6vPzjj4iGepBcfPI7YDdgGPBbYEKt4yrzGFYCu9Q6jjLiPRw4AHgks+zbwFnp87OA82sdZx+O4VzgC7WOrYTYdwUOSJ+/HngCmFAv34Me4q+Lzz+NW8CO6fOhwAPAIcCNwLR0+Q+AT9U61hw+i9xz8ED5HQDOAK4Hbk1f5/b9J7kD4j+kz4cBO+d1/EAz8BSwXea4T6r28Zfztwf4EPDT9Hf1EOCBKu3/KGBI+vz8zP4npL8Lw4Hx6e9IU6X3ny5vJbkY92nSWiLH438/8DNgePr6TXkeP3AHcHTmmH/e1+NvxB7kUm6nahUUEb8gueI9K3vL2quAj+YaVJm6OYa6EBHPRsSD6fM/ActI/mDVxfegh/jrRiQ2pi+Hpo8AjiC5ZTMM4O9BheWegwfC74CkFuDDwGXpa5HT91/SG0iKhcsBIuLliPgj+eaAIcB2SubP3h54lioff5l/e6YAV6e/q/cDO0vatdL7j4g7ImJL+vJ+kjnDO/df0du19/B367vAmWx7M59cjh/4FPCtiPhr2ua5zP7zOP4A3pA+34nX5msv+/gbsUBuBlZnXndQZ39sSb7Bd0harOS2rfXozRHxLCR/vIA31TievpqVno6ZU+3Ts5UgaRywP0kPZt19Dwrihzr6/NPT60uA54A7SXpI/pj5Y1mPuagvapqDa/g78D2SouSV9PUo8vv+7wasBa5Ih3hcJmkHcjr+iHgGuBBYRVIYvwAspjY//90dcy1+Lk8h6bXMbf+SjgWeiYjfFqzK6/jfBrwnHVpzr6SDct7/Z4ELJK0m+Zk8u6/7b8QCueRbog5g746IA4CjgU9LOryelbM4AAAgAElEQVTWAQ1S/wG8FdiPJOl/p7bh9EzSjsB/AZ+NiBdrHU+5isRfV59/RGyNiP1IeowmAu8o1izfqGqiZjm4Vr8Dko4BnouIxdnFRZpW63MYQnKq+T8iYn/gzyTDC3KR/vM6heTU+RhgB5K/X4Vq+fOf68+lpC+RzBl+XV77l7Q98CXgnGKrq73/1BBgJMkwhn8BbkzPpuS1/08Bn4uIVuBzpGdV+rL/RiyQ6/6WqBGxJv36HHAz/TwNUSN/6Dx9kX59rpf2A05E/CEtel4BfsgA/j5IGkpSGFwXETeli+vme1As/nr6/LPSU9s/J/kDsXN6yhnqMBf1UU1ycI1/B94NHCtpJcmQkiNIepTz+v53AB0R0XnmZR5JwZzX8X8AeCoi1kbEZuAm4FBq8/Pf3THn9nMpaQZwDPCJSAfA5rT/t5L8k/Lb9GexBXhQ0t/ktH/S/dyUDmX4NckZlV1y3P8Mkp8/gB/z2t+NsvffiAVyKbdTHbAk7SDp9Z3PSQb8d7lCtQ5kb1k7A/jvGsbSJwXjk/6WAfp9SP87vxxYFhEXZVbVxfegu/jr5fMHkDS682p1SduRFAzLgHtIbtkMA/h7UGG55+Ba/w5ExNkR0RIR40iO9+6I+AQ5ff8j4vfAakl7posmkdz1MK8csAo4RNL26feic/+1+Pnv7pjnA9PT2QwOAV7oHIpRSZImA/8XODYi/lIQV1Vv1x4RD0fEmyJiXPqz2EFy8ervyen4gVtI/kFE0ttILhhdR363q18DvDd9fgTwZPq8/OPv7Sq+enyQXK34BMkYwC/VOp4yY9+N5ErP3wJL6yF+4Eckp8A3k/xCnkoy/u6u9IfzLuCNtY6zD8dwDfAw8FD6y7VrrePsJvbDSE4VPQQsSR8fqpfvQQ/x18Xnnx7DvsBv0lgfAc5Jl+9G8kdgOUlvxvBax5rT55FrDh5IvwPA+3htFovcvv8kQ5Ha08/gFpLT3LkdP/A14LH05/8aktkKqnr85fztITnFfmn6M/kw0Fal/S8nGeva+XP4g0z7L6X7f5x0poVK779g/Upem8Uir+MfBlyb/hw8CByR5/GnuWAxSQ31AHBgX4/ft5o2MzMzM8toxCEWZmZmZmZ95gLZzMzMzCzDBbKZmZmZWYYLZDMzMzOzDBfIZmZmZmYZLpBtUJN0q6QrS2h3paRzMq9XSvpCVYNL9nOSpI29tJklqW7m+jYzMxvoXCBbVaU3UPj3tKD8q6Q/SLpL0pG1jq1UkvYhuY3q93La3zBJL0h6Z4mb/BBok/SeasZlZo1P0v2SLsxru0qR9M+SFlTgfX4vaVYlYqo2SZdI+nat42hULpCt2v6L5FaPpwJvI7n95k9JJnOvF58B/isiXsxpf+8Hno+I35bSOCL+ClwPnFbVqMysXyRFL48rax1jb9JCdF2RVR8iuVlI7iRtD5xbof3vA8ypwPu8StLb0+/v3pV8X+BfgU9Jau21pZXNBbJVTXrr3fcAZ0XEXRHxdEQsiogLI2Jupt1KSedKulbSxvQ/+C8UvNdOkmZLek7SnyTdK6mtoM2h6fK/SHpG0n9IekNm/fbpUImNaU/2F0s4hibg7+jlVrmS/l7Si5KOTV//PN3/dyRtkLRW0unpbTYvlfRHSaskfbLI202h4JaskiZJekTSnyXdk96qM2s+cGz6h8LMBqZdM49/LLLs9GIbSRqaS3T9EBEbIuJPNdr9NGBdRDzQ3zeKiLWx7S2iB6yIWAPcC8ysdSyNyAWyVdPG9HGspBG9tD0DWAYcAHwV+FdJHwOQJOA2oJmkB3p/4BfA3ZJ2TdvsA9xBUii+E/gYya1Xsz0BFwJHAh8HJqXvc3gvce0L7ERyC9eiJJ0G/D/gmIjIFtKfAP4EHAx8i2SIxi0kt+BtA64CLpM0JvNeAo5l2wJ5OHA2cArwLmBn4AcFYbQDQ9L1ZjYARcTvOx/AHwuXRcQLmd7G49J/+DcBMyS9WdIN6T//f0n/Yf5E53un/4CvTnMImeU3Sboh8/pjkn4jaZOkFWnnREkFuKTJwH8AozK93mel67YZYpF2dJyd6fh4Ot33GyXNS5c9Jul9BfvYR9LtmY6MayWN7iW0EynoxJA0N93PV9L3+aOk8yQ1SfpXSeskPSvpswXbvTrEQtKI9BhPlnRz+rn/TtLfZdp36R3ObHdM+rdvWbrq4XT57Zm2M9PPYVP6dVb2eyjpM5KWKxmiuFZdh5HMB07o5fOxvqjkfdH98KPwQVKMbgA2AfeRFKkHF7RZCdxZsOwy4Jfp8yNICu3tCtosAc5Mn18NXF6wfj8ggDcBOwJ/BT6RWb8jyR+pK3uI/6PAK8DrisT8BeA84A/A/gXrfw7cl3ktYC0wP7NsKPAyMDWz7KD08xqSvj4pPYY9M20+kW5XGNMG4NRaf8/98MOP3h/A1ORPcJflb09/53+X5p/xwBhgHPC5NK/tBnwa2AIclm73JmAz8L7Me+2c5t6PpK+PTXPe9PQ9PpDu5xuZbe4HLuwm5mHAmcB64G/Sxw7FtgN+D6wj6SnfnaQT4c/AApKCdg/gGuAZYFi6TSvwPPB1YM/0WG8H/qeHz7Epfd+PFSyfC7wIfD99rxnp53o7yVCMPUh67F8B9i6Ie1b6fES6zSqSXurdgYvSz3TXgu9X9j06tzsmfX1Y+vq96Wc2Ml3+GaAD+Nv0+/xRkr8T/5Cuf3ea648H3pJ+Hp8vOM7Ov3PNtf6ZbrSHe5CtqiLiv0iS+0dIxh4fCtxfZHjDfUVeT0ifHwhsD6xNexU2KpnZYW/grZk2f1+w/lfpuremj2HZ/UTERuDhXg5hO2BzRLxSZN3pJON+D4uI3xRZ/1BmXwE8l91fRGwm+WPwpsw2U4DbImJLZtlfI+LxzOs1JMX1zgX7eymN18zq30URcUtEPBURayJiZUR8NyKWRMSKiLgUuJWkcCMingN+RvIPdKepJJ0LnT2WXwa+GRFXp+/xM+BLwP8pJaCIeJmk6Ix4rdf7zz1s8pOI+GFELCcpSrcHHomI6yPiSeCbJH8f3pa2/wzwvxHxlYh4PCKWkHQSHCZp3272MTp932eLrFsLfDZ9r6uApcAbI+KrEfFkRFxM0sHxvl4O/fKImJsex9kkRfm7e9kmq3PM9vr0M3s+7SX+EvC5iLg5/T7fAnyH174fbyH5vH8SyRDFJRHxnYL3XpN+HVdGPFYCF8hWdRGxKSLujIjzIuJQ4HLgXEnDSnyL15Eksf0KHm8HvpJpc1nB+neS9BIsIenB7Yt1wLBuxvb+kuQ/9+5Ob20ueB3dLMv+Hn6UgvHHJL1EhdtA19/fN5L8QTCz+rfNsC5JQyR9VdLDSq5r2Ah8GBibaXYtMFXS8PT1J4AbImJzWpDtD5xX0JEwBxgpaWQVjiHbSbAO2Mq2nRJ/SL92dhIcCBxZEN/ydN1bKa6zU2BTkXWPFHRu/IGunSLPsW0nRTHZ4/grydm63rbpTQvwZuCqguM9l9eOdUEa30pJ1yi51mWHgvd5Kf3qzpEKG1LrAGxQepTkZ28EyekjgEMK2hzCa+O2HiRJJK9ExIpu3vNBYK/0P/wuJC0nKU4PAVaky3Yg6YX+XQ+xLkm/TqDrOOTFJKfb7pQUEfH1Ht6nV5LeSnIK7/be2naz7QiSz8HM6l9hz+yXSIZVfJakJ/TPJL2NwzNtbia5PuFDkhaRXGPRebZOJP9Uf5mu/4RD0lNZaYUdAoXLCv/Zfx3JdRrFLqD+fTf76OydLVbg96WTopietukswLOdMKWM6e7c/mSSvyVZrwBExB+VTPf5PpLhMOcA35R0UHrGAJKOEXDnSMW5QLaqkTQK+DFJD8VDJBestZGMYbsrtp027RBJZwPzSJLBdF47VfgzkuES/y3pTOAxknFck4GfRcT/AOeTDN34AfCf6b7eTjL27p8iYqOky4HzJa0lOS11Dsmpsm5FxFpJD5KMIetyoV5ELJJ0FHBHWiR/o6wPaVtTSD6XHm8M0o33ACvS05Zm1ngOA26OiOsBJL2OZGjC050NIuIvkm4hyZ27A09FxH3pulckLQHe1l1HQoleppe82Q8PkuT1pyJiaykbRMSfJP2OpBPj7irF1ZPOwnRXXuud3q+gTWdHUPZzW01S3O8WETfQjXQo3p0kHTHnkoz/nkxy3Q0knTwv8VqHklWIC2Srpo0kF26cTpKsh5NckHE9UFhIXkQyY8SXSHpGzomIeZAMdpP0oXSbH5Kc2voDSdF8ddrmIUmHp23uJUlEK0h6VDp9AdghXfYXkotGCk9XFTMb+Ae6uVFIRPw6UyTTjyL5oySnSPviBJLPxswa0xPAhyW9i+RCuzNIxu8+XdDuWpIe4r2B6wrWfQ34L0nPkMxR/wrJvL/7RUSv016mVgI7SXovaU92RLzU8yYlu5ikR/X6dEaM9SR/O44HPpWOgS5mIck/EJdUKI6SpeOJlwBnS1oN7EIya1HWsyRF8mRJzwKbIuJFSV8j6bTZSHIMw0mGmYyOiAuUzOQ0hmQ43/PAUSRnCrPF8HuAe3r4bKyv8r4q0A8/Ch+kM0LUOo4e4hsOPAW8p4r72IXkNN7f9GHbvUn+Ydip1p+VH374UdqD3mex2Ltg+S4kU3ptTH/fv0lyPcftBe2aSM6QBbBHkff/EPC/JL2OLwC/Bv4ps77bWSzS9Ur3uz7dx1nFtiMzG0Rm2RZgWub1zul7fKDg+G8h+SfgLyRnDC8GmnqI6Z0kHSuvzyybC8wraPcz4AcFy5aw7SwexWaxOKZgm22OjeSfjPvTeJeQnAXdZjuSC+9Wk4zDvj2zfHq6TefY5l8AH0/XvZ+kw2dD+t4PAX9f8L1YCXy01j/PjfhQ+iGb1YyklcAlEVGz25T2Ju2dHhkRxcbuVeL93wZ8MCL+Xx+2PQpQRCysfGRmZgNfOrTkVxFxQa1jyYukj5NcqH5AFJ9pyfrBQyzMShARv6jy+z9Bcgq1L9veUeFwzMzqzRnA0bUOImfDSea+d3FcBe5BNjMzMzPL8DzIZmZmZmYZLpDNzMzMzDJcIJuZmZmZZbhANjMzMzPLcIFsZmZmZpbhAtnMzMzMLMMFspmZmZlZhgtkMzMzM7MMF8hmZmZmZhkukM3MzMzMMlwgm5mZmZlluEA2MzMzM8twgWxmZmZmljGk1gFUyy677BLjxo2rdRhmNkgtXrx4XUSMrnUcteQ8bGa10t8c3LAF8rhx42hvb691GGY2SEl6utYx1JrzsJnVSn9zsIdYmJmZmZlluEA2MzMzM8twgWxmZmZmltGwY5DNbGDYvHkzHR0dbNq0qdahVMWIESNoaWlh6NChtQ7FzKwL5+C+cYFsZlXV0dHB61//esaNG4ekWodTURHB+vXr6ejoYPz48bUOx8ysC+fgvsl1iIWkyZIel7Rc0llF1h8u6UFJWyRNLbL+DZKekXRJPhGbWX9t2rSJUaNGNVxiBpDEqFGj6qZnxjnYbPBxDu6b3ApkSU3ApcDRwATgBEkTCpqtAk4Cru/mbb4O3FutGM2sOhoxMXeql2NzDjYbvOolT/VFtY4tzx7kicDyiFgRES8Dc4Ep2QYRsTIiHgJeKdxY0oHAm4E78gjWzKzBOAebmZUozwK5GVided2RLuuVpNcB3wH+pZd2MyW1S2pfu3ZtnwM1s+ppbh2LpIo9mlvH9rrPpqYm9ttvP975zndywAEH8L//+78ArFmzhqlTk5EEV155JbNmzarqsddY1XNw2tZ52GwAcw4uTZ4X6RXrA48St/0/wIKIWN1TV3pEzAZmA7S1tZX63maWozUdq5l58a0Ve7/Zpx/Ta5vtttuOJUuWALBw4ULOPvts7r33XsaMGcO8efMqFssAV/UcDM7DZgOdc3Bp8uxB7gBaM69bgDUlbvsuYJaklcCFwHRJ36pseGY2GLz44ouMHDkSgJUrV7L33nt3aXPbbbfxrne9i3Xr1vH0008zadIk9t13XyZNmsSqVavyDrlSnIPNrObqJQfn2YO8CNhD0njgGWAacGIpG0bEJzqfSzoJaIuILldgm5kV89JLL7HffvuxadMmnn32We6+++5u2958881cdNFFLFiwgJEjR3LyySczffp0ZsyYwZw5czjttNO45ZZbcoy+YpyDzawm6jEH59aDHBFbgFnAQmAZcGNELJV0nqRjASQdJKkDOA74T0lL84rPzBpX5+m9xx57jNtvv53p06cT0fXs/z333MP555/Pbbfd9moPx3333ceJJyZ15Cc/+Ul++ctf5hp7pTgHm1mt1GMOzvVGIRGxAFhQsOyczPNFJKf9enqPK4ErqxCemQ0Cnaftil1Atttuu7FixQqeeOIJ2traim5fz9MlOQebWa3VSw7O9UYhZma19thjj7F161ZGjRrVZd1b3vIWbrrpJqZPn87SpUnn6aGHHsrcuXMBuO666zjssMNyjdfMrJHUSw72rabNLFdjWlpLuuq5nPfrTef4N0huTXrVVVfR1NRUtO2ee+7Jddddx3HHHcdPfvITvv/973PKKadwwQUXMHr0aK644oqKxW5mljfn4NKo2BiQRtDW1hbt7e21DsMsd82tY1nTsbr3hqkxLa08s7p6VwUvW7aMd7zjHVV7/4Gg2DFKWhwRxc8RDhLOw2a15xzcN+5BNmsw5c5xWcmeBDMzs0bgMchmZmZmZhkukM2s6hp1KBc09rGZWWNo5DxVrWNzgWxmVTVixAjWr1/fkAk6Ili/fj0jRoyodShmZkU5B/eNxyCbWVW1tLTQ0dFRdM7LRjBixAhaWnqcOtjMrGacg/vGBbKZVdXQoUMZP358rcMwMxuUnIP7xkMsjObWsUgq+dHcOrbWIZuZmZlVjXuQzdOCmZmZmWW4B9nMzMzMLMMFspmZmZlZRq4FsqTJkh6XtFzSWUXWHy7pQUlbJE3NLN9P0n2Slkp6SNLxecZtZtYInIPNzEqTW4EsqQm4FDgamACcIGlCQbNVwEnA9QXL/wJMj4i9gMnA9yTtXN2Izcwah3OwmVnp8rxIbyKwPCJWAEiaC0wBHu1sEBEr03WvZDeMiCcyz9dIeg4YDfyx+mGbmTUE52AzsxLlOcSiGVided2RLiuLpInAMOB3RdbNlNQuqb1RJ8Q2M+ujqufgdL3zsJnVvTwLZBVZVtZ9DyXtClwDnBwRrxSuj4jZEdEWEW2jR4/uY5hmZg2p6jkYnIfNrDHkWSB3AK2Z1y3AmlI3lvQG4DbgyxFxf4VjMzNrdM7BZmYlyrNAXgTsIWm8pGHANGB+KRum7W8Gro6IH1cxRjOzRuUcbGZWotwK5IjYAswCFgLLgBsjYqmk8yQdCyDpIEkdwHHAf0pamm7+d8DhwEmSlqSP/fKKvR6Vc/toM2t8zsFmZqXL9VbTEbEAWFCw7JzM80Ukp/0Kt7sWuLbqATaQcm4f7VtHmw0OzsFmZqXxnfTMzMzMzDJcIJuZmZmZZbhANjMzMzPLcIFsZmZmZpbhAtnMzMzMLMMFspmZmZlZhgtkMzMzM7MMF8hmZmbWJ+XclKq5dWytwzUrWa43CjEzM7PG4ZtSWaNyD7KZmZmZWYYLZDMzMzOzDBfIZmZmZmYZuRbIkiZLelzScklnFVl/uKQHJW2RNLVg3QxJT6aPGflFbWbWGJyDzcxKk1uBLKkJuBQ4GpgAnCBpQkGzVcBJwPUF274R+CpwMDAR+KqkkdWO2WwwaBoy1FehDwLOwYOXZ5owK1+es1hMBJZHxAoASXOBKcCjnQ0iYmW67pWCbT8I3BkRG9L1dwKTgR9VP2yzxrZ1y2ZfhT44OAcPUp5pwqx8eQ6xaAZWZ153pMsqtq2kmZLaJbWvXbu2z4GamTWgqudgcB42s8aQZ4GsIsuikttGxOyIaIuIttGjR5cVnJlZg6t6DgbnYTNrDHkWyB1Aa+Z1C7Amh23NzMw52MysZHkWyIuAPSSNlzQMmAbML3HbhcBRkkamF4YclS4zM7PSOAebmZWorAJZUp8v6ouILcAskqS6DLgxIpZKOk/Ssen7HySpAzgO+E9JS9NtNwBfJ0nwi4DzOi8WMTMbLJyDzczyUW6yfVbSVcDlEbGs3J1FxAJgQcGyczLPF5Gcuiu27RxgTrn7NDNrIM7BVlWd0z6aDXblFshfBE4GPifp18BlwA0RsbHikZmZWSHnYKuqcqZ9BE8LZ42rrCEWEfHDiDgU2Bv4JfANkh6NOZLeXY0Azcws4RxsZpaPPl2kFxHLIuJfSE7FfRE4EfiFpMck/bOkXG9hbWY2mDgHm5lVV58u+EivgP4YcApwBElPxuXAGOArwPtIrpA2M7MKcw42M6uusgpkSQeQJOQTgM3A1cCnI+LJTJu7gP+pZJBmZuYcbGaWl3J7kH8N3AnMBP47nTao0DJgbn8DMzOzLpyDzcxyUG6B/NaIeLqnBhHxZ5KrrM3MrLKcg83MclDuhRz3SBpVuFDSzpJWVCgmMzMrzjnYzCwH5RbI44CmIsuHA839jsbMzHoyDudgq1OdNyEp9dHcOrbWIdsgVtIQC0kfy7z8sKQXMq+bgEnAygrGZWZmKedgawS+CYnVk1LHIM9LvwbJVEJZm0kS8+crFJOZmW3LOdjMLEclFcgR8ToASU8BB0XEuqpGZWZmr3IOtqzm1rGs6Vhd6zDMGlpZs1hExPhqBWJmZj1zDjaANR2rPVTBrMp6LZAlnQH8e0RsSp93KyIu6uW9JgMXk4yZuywivlWwfjjJxPcHAuuB4yNipaShwGXAAWnMV0fEv/UWu5lZvXMONjPLXyk9yJ8BrgI2pc+7E0C3yVlSE3ApcCTQASySND8iHs00OxV4PiJ2lzQNOB84HjgOGB4R+0jaHnhU0o8iYmUJ8ZuZ1TPnYDOznPVaIGdP6fXz9N5EYHlErACQNBeYAmST8xTg3PT5POASSSJJ/DtIGgJsB7wMvNiPWMzM6oJzsJlZ/sqdB7mL9NRbKZqB7FUFHXSdt/PVNuktVF8ARpEk6j8DzwKrgAsjYkORWGZKapfUvnbt2rKOw8ysHg2kHJzG4zxsZnWvrAJZ0mmSPp55PQd4SdLjkvbsbfMiy6LENhOBrcAYYDzweUm7dWkYMTsi2iKibfTo0b2EY2ZWXwZ6Dgbn4b5qbh1b8g00zKz6yprFAjgNOAVA0uEk49JOBD4OfAfo6VLZDqA187oFWNNNm470VN5OwIZ0H7dHxGbgOUm/AtoA31rVzAYT5+AGVc7MFJ6Vwqz6yh1i0cxrd2v6CPDjiLiRZMzaIb1suwjYQ9J4ScOAacD8gjbzgRnp86nA3RERJKf0jlBih3Rfj5UZu5lZvXMONjPLQbkF8otA5zmzI4G70uebgRE9bZiOZ5sFLASWATdGxFJJ50k6Nm12OTBK0nLgDOCsdPmlwI7AIyRJ/oqIeKjM2M3M6p1zsJlZDsodYnEH8ENJvwF2B36aLt8LeKq3jSNiAbCgYNk5meebSE4ZFm63sdhyM7NBxjnYzCwH5fYgfxr4FbALMDVzFfMBwI8qGZiZmXXhHGxmloNybzX9IkUmqo+Ir1YsIjMzK8o52MwsH+UOsQBA0hjgTRT0QEfEg5UIyszMuuccbGZWXWUVyJL2B64F3k7X+TIDaKpQXGZmVsA52MwsH+X2IM8mucvSP5LMn1k4ybyZmVWPc7CZWQ7KLZAnAPtHxBPVCMbMzHrkHGxmloNyZ7F4GPibagRiZma9cg42M8tBuQXyF4FvS/qApDdLemP2UY0AzczsVc7BZmY5KHeIxc/Sr3ew7dg34QtEzMyqzTnYrIjm1rGs6VhdcvsxLa08s3pVFSOyeldugfz+qkRhZmalcA62QaNpyFCkwslaujfz4ltLbjv79GP6EpINIuXeKOTeagViZmY9cw62wWTrls0lF70ueK3Syh2DjKR9JF0i6aeSdk2XfTSdn9PMzKrIOdjMrPrKKpAlHQUsApqBI4Dt0lVvBXq91amkyZIel7Rc0llF1g+XdEO6/gFJ4zLr9pV0n6Slkh6WNKKc2M3M6p1zsJlZPsrtQf46cEZE/C3wcmb5z4GJPW0oqQm4FDiaZC7PEyRNKGh2KvB8ROwOfBc4P912CMndo/45IvYC3gdsLjN2M7N65xxsZpaDcgvkvYAFRZZvAHqbYmgisDwiVkTEy8Dc/9/encfZUdV5H/98SSCggywSFLMQJJHdBcKmoFEeWdwiY5SAIlEUcIzzjMqoOIoRRUFwHdzQMAHElcEhA9E4EaKjIhIWDQTBAJE0YQnLQ4gmgcDv+eOcm1RX7u2+t9N9b3f19/163Vd3VZ2qOlV161fnnnOqCphaSjMVuDj/fzlwhFIP/SOBP0XEHwEi4pGIeLrFvJuZDXWOwWZmbdBqAfkxUtNe2f5AVy/zjiG9IrWmq86yNqSJiPXA48BzgRcBIWm+pJskfaTeCiSdImmRpEUrV67sdWPMzIaYQR2DwXHYzKqh1QLy94HzJI0lPXNzpKRXAecDl/Qyb71ntUSTaUYChwFvz3+PlXTEJgkjLoyIyRExefTo0b1kx2xoGDNuPJKa/lilDeoYDI7DZlYNrT4H+RPAHOCvpEC6hFTIvgw4u5d5u4BxheGxwIoGabpyn7ftSE2HXcCvIuJhAEnzSDUmv2wx/2ZDzoqu5X6+p9U4BpuZtUFLNcgR8VREvB2YBLwNOAHYIyJObKI/2g3AJEm7SdoKmA7MLaWZC5yU/58GXBMRAcwHXizpWTlov4p0YTAzGzYcg83M2qPpGmRJ2wAfAd4CvJDU7HY38BNJX4yINT3NHxHrJc0kBdoRwEURcZuks4BFETEXmA1cKmkpqdZiep73MUlfIgX4AOZFxNUtbqv1k1bfbuRXepptPsdgM7P2aaqAnGsMriE1qf0cuJrUvLc3cCZwjKRX5Zs6GoqIeZTuwI6IMwv/rwXe2mDe75EeM2Qd1srbjcBN/mabyzHYzKy9mq1BPgWYCOwfEbcVJ0jaF4BsNhwAAByVSURBVLg2p/lG/2bPzMxwDDYza6tm+yBPA84uB2aAiLgV+HxOY2Zm/c8x2MysjZotIO9Dat5rZAGw7+Znx8zM6nAMNjNro2YLyDsAPT3xfSWw/eZnx8zM6nAMNquoVp51P2bc+E5nd9hotg/yCKCnmz+eyWnMzKz/OQabVVQrz7r3Te/t02wBWcD3JK1rMH1UP+XHzMw25RhsZtZGzRaQL24iTW+vOTUzs75xDDYza6OmCsgR8a6BzoiZmdXnGGxm1l4tvWrazMzMzKzqXEA2MzMzMytwAdnMzMzMrMAFZDMzMzOzAheQzczMzMwK2lpAlnS0pDskLZX0sTrTR0n6UZ5+vaQJpenjJa2WdHq78mybb8TILf2WILNBwDHYzKw5zT4HebNJGgF8HXgt0AXcIGluRCwpJDsZeCwiJkqaDpwLHFeY/mXgZ+3Ks/WPp9c/5bcEmXWYY7CZWfPaWYN8ELA0Iu6OiCeBHwJTS2mmsvGB+JcDR0gSgKQ3A3cDt7Upv2ZmVeIY3EZjxo1vuuUs72JrI7dsWm/aVoMMjAGWF4a7gIMbpYmI9ZIeB54raQ3wUVLNR8OmPUmnAKcAjB/vL7SZWcGAx2BwHK5Z0bW86ZYzcOtZu7ll03rTzhrkej+Ro8k0nwa+HBGre1pBRFwYEZMjYvLo0aP7mE0za6SVWhfXvAw6Ax6DwXHYzKqhnTXIXcC4wvBYYEWDNF2SRgLbAY+SajmmSfoCsD3wjKS1EXHBwGfbqmLMuPGs6Free8LsBWPHcd/yewcwR0NPK7Uu4JqXQcYx2MysSe0sIN8ATJK0G3AfMB04oZRmLnAScB0wDbgmIgI4vJZA0ixgtQOzQeuF3lYKd7M/fGzTfQNdmLYhwDHYzKxJbSsg5/5sM4H5wAjgooi4TdJZwKKImAvMBi6VtJRUazG9XfmzoamVfn6t1ma2UlvaSmEaYMtRW/PUurUt5cdsczgGm5k1r501yETEPGBeadyZhf/XAm/tZRmzBiRzZpuhL10PfIOItZtjsNnQVrsPpFlu3ey7thaQre9a7UpgZmZm1eL7QNrHBeQhYrg8MqjVX8dmZmZWXyvXVNc2d+cCsg0q/nVsZmbWP/y8575r53OQzczMzMwGPReQzczMzMwKXEA2MzMzMytwAdnMBkwrr6b2a6nNbDBqJY5JYqutt2kpvQ1OvknPzAaMbxAxs6FuIJ9zX0tvg49rkM3MzMzMClxANjMzMzMrcAHZzAaFVvv5uc+ymZkNFPdBNrNBwS+JMTOzwaKtNciSjpZ0h6Slkj5WZ/ooST/K06+XNCGPf62kGyUtzn9f0858m5lVgWOwmVlz2lZAljQC+DpwDLA3cLykvUvJTgYei4iJwJeBc/P4h4E3RsR+wEnApe3JtZlZNTgGm5k1r501yAcBSyPi7oh4EvghMLWUZipwcf7/cuAISYqImyNiRR5/G7C1pFFtybWZWTU4BpuZNamdBeQxwPLCcFceVzdNRKwHHgeeW0rzFuDmiFhXXoGkUyQtkrRo5cqV/ZZxM7MKGPAYDI7DZlYN7Swg13tdTLSSRtI+pCa/U+utICIujIjJETF59OjRfc6omVkFDXgMBsdhM6uGdhaQu4BxheGxwIpGaSSNBLYDHs3DY4GfAu+MiLsGPLdmZtXiGGxm1qR2FpBvACZJ2k3SVsB0YG4pzVzSDSAA04BrIiIkbQ9cDZwREb9tW47NzKrDMdjMrEltKyDn/mwzgfnA7cCPI+I2SWdJelNONht4rqSlwIeA2mOIZgITgU9KuiV/dm5X3s3MhjrHYDOz5rX1RSERMQ+YVxp3ZuH/tcBb68z3WeCzA55BM7MKcwzuuzHjxrOia3nvCc2sEvwmPTMzs16s6FruNz2aDSNtfZOemZmZmdlg5wKymZmZmVmBC8hmNiSNGLklkpr6jBk3vtPZNTOzIcR9kM1sSHp6/VNN9wl1f1AzM2uFa5DNzMzMhrlWWuWGQ8uca5DNzMzMhrlWWuWg+i1zrkE2MzMzMytwAdnMzAbEjTfeOOyaZc2sGtzFwszMBsRO4ybyj6d/pdu4qjfLmlk1uAbZzMyGpTHjxjd9Q5KZdVf1R226BtnMzIalVl4f7Zpvs+6q/qhN1yB3kGsvzMzMzAafttYgSzoa+CowAvhuRJxTmj4KuAQ4AHgEOC4iluVpZwAnA08D/xwR89uY9QHh2gsza6fBEINrzbLNesHYcdy3/N6+rMrMrM/aVkCWNAL4OvBaoAu4QdLciFhSSHYy8FhETJQ0HTgXOE7S3sB0YB/gBcACSS+KiKfblX8zG7pcKBs8MbjVZ63O/vCxTR+7LUdtzVPr1raaJTOzTbSzBvkgYGlE3A0g6YfAVKAYnKcCs/L/lwMXKEXGqcAPI2IdcI+kpXl51zVaWe3xQmWtBtBW0js4mw1OfgA+0OYY3F9a7efo42xm/UER0Z4VSdOAoyPiPXn4RODgiJhZSHNrTtOVh+8CDiYF7N9HxPfy+NnAzyLi8tI6TgFOyYP7ArcO6EYNDjsBD3c6E23g7ayO4bCNAHtExLadzkRNO2Jwnjbc4vBw+T57O6tjOGwjbGYMbmcNcr02snLpvFGaZuYlIi4ELgSQtCgiJreayaHG21ktw2E7h8M2QtrOTuehZMBjMAy/ODwcthG8nVUyHLYRNj8Gt/MpFl3AuMLwWGBFozSSRgLbAY82Oa+ZmTXmGGxm1qR2FpBvACZJ2k3SVqQbPuaW0swFTsr/TwOuidQHZC4wXdIoSbsBk4A/tCnfZmZV4BhsZtaktnWxiIj1kmYC80mPGLooIm6TdBawKCLmArOBS/MNII+SAjg53Y9JN5OsB97fxN3TFw7Utgwy3s5qGQ7bORy2EQbZdnYgBsMg2wcDZDhsI3g7q2Q4bCNs5na27SY9MzMzM7OhwG/SMzMzMzMrcAHZzMzMzKygkgVkSUdLukPSUkkf63R+BoqkZZIWS7plED5Sqs8kXSTpofxM1tq4HSX9j6S/5L87dDKPm6vBNs6SdF8+nrdIel0n89gfJI2TdK2k2yXdJun/5vGVOZ49bGPljmezHIOHNsfg6pyzwyEGw8DE4cr1QVZ6neqdFF6nChxfep1qJUhaBkyOiEo98FvSK4HVwCURsW8e9wXg0Yg4J19wd4iIj3Yyn5ujwTbOAlZHxPmdzFt/krQLsEtE3CRpW+BG4M3ADCpyPHvYxrdRsePZDMfgoc8xuDrn7HCIwTAwcbiKNcgbXqcaEU8Ctdep2hAREb8m3UFfNBW4OP9/MemLP2Q12MbKiYj7I+Km/P8TwO3AGCp0PHvYxuHKMXiIcwyujuEQg2Fg4nAVC8hjgOWF4S6qe7EK4BeSblR6vWuVPS8i7od0IgA7dzg/A2WmpD/l5r8h3eRVJmkC8DLgeip6PEvbCBU+nj1wDK6mSp6zdVT2nB0OMRj6Lw5XsYDc9CtRK+AVEbE/cAzw/txkZEPXN4HdgZcC9wNf7Gx2+o+kfwD+E/iXiFjV6fwMhDrbWNnj2QvHYBuqKnvODocYDP0bh6tYQB42r0SNiBX570PAT0lNm1X1YO5jVOtr9FCH89PvIuLBiHg6Ip4BvkNFjqekLUkB67KIuCKPrtTxrLeNVT2eTXAMrqZKnbP1VPWcHQ4xGPo/DlexgNzM61SHPEnPzh3RkfRs4Ejg1p7nGtKKr8A9Cbiyg3kZELVglR1LBY6nJJHeznZ7RHypMKkyx7PRNlbxeDbJMbiaKnPONlLFc3Y4xGAYmDhcuadYAOTHeHyFja9TPbvDWep3kl5IqrGA9Mrw71dlOyX9AJgC7AQ8CHwK+C/gx8B44F7grRExZG+waLCNU0jNQAEsA06t9REbqiQdBvwvsBh4Jo/+OKlvWCWOZw/beDwVO57Ncgwe2hyDq3PODocYDAMThytZQDYzMzMz66sqdrEwMzMzM+szF5DNzMzMzApcQDYzMzMzK3AB2czMzMyswAVkMzMzM7MCF5ArRNIcSVe1OM/pkpYNUJasH0maJaljz+WUNEnSg5K2y8MzJK3uVH56I+kPkv6x0/mw4aEv8bcf171a0oxOrHs4kbRM0ukdXP8sSRcVhhdKuqBT+emJpP0k3ZefET4kuYDcjySdJulv+eH4tXFbSfq7pMWltJMkhaTXtD+n7SVpQt7W2me1pDskfVfSizudv6JOF0JzHmr7a3Jp0vnAqzqRp+xzwDci4vEO5qEVnwHOleQ4Nww4/jaW41ot/q6X9Kik30k6I7+ad9DodCE056HRdeBA4Bvtzg+ApJ2BDwGf7cT6WxURi4Hfk/I8JPnC0b+uAZ5F91cZHgw8DrxI0ujC+CnAOuB3bctd5x0N7ALsB3wQ2Bm4UdL0juYqy6+p7M/lbdV7quZFxOqIeKQ/l9ksSeOANwP/0Yn1lzW5b+cB2wLHDHB2bHBw/O3ZHaT4Ow44HLgYOBW4WdLzO5kx6P94KWkLSSP6c5kRsTIi/t6fy2zBe4A/RMTdHVr/Bi3s2/8A3idp5EDnaSC4gNyPIuJOYAXw6sLoVwMLgEWkoFwcf11ErIUNNR3nSurKtSA3SDqquHxJe0u6WtITkh6S9IOeApukl0i6X9LZhXEfkfRArsW9BPiH0jwHSvqFpIclrZL0G0mHFqZfVG5GzCfLvZJ6+6X4SEQ8EBH3RMS8iHgT8BPgW5K2Lyzv5ZJ+lWt+7pP0TUnPKUxfKOlbkr4q6bH8Oa9YUyjpHXkf1vbVTySNKUyfkmtTXpeb4p8kXSw+BexTqG2ZkdOHpGml7e5W05HTvF/SFZL+Rqpxbfm4AffkvzfkZS7My+lWq6HcpCvpo/mYPi7pnHw8ZuV1PSDpo6V8byfpwjz9ibyvy7XVZccBiyPi3vIESUdIujV/b6+VtFtp+qmSlkp6Mv99b2l6n/atpC0lfU3SCknrJC2XdE5tnoh4mlRIPr6XbbMKGCLx942SbpS0VtI9ks5W9xrvZZI+IenbOf52SfrX0nIn5hi4Vqkl7g1N7qL1Of7eHxG3RcS3gUOBHYFzC8uX0nXiLklrJC2W9I7C9FoL1wlK14e1kv4s6chCmhGSZudtXCPpL3mZxRhdjF9dQFeOdbsC59VicE67SXcubYzhOxXTKMX0W4Engb3ytHdJWpLzeqekD6pBy5JSzG90HagXl94n6Uql69Wdkl4taayk+fm7dIuk/Uvr6PEa18AJ1H9l+xaSPqd0zX5I0vml/byDpIuVrpNrJC2QtE9xe/u6b5W6Ufwyf1efkPRHScXz7xek79eUXrZtcIoIf/rxA3wPuKYwfC1wMnA2qXm6Nn4FcGZh+DJSc8QrgRcCM0lfwpfk6bsAD5MC2V7Ai4H/Bv4AbJHTzAGuyv8fDjwKfLCwjrflZZ4KvAj4N2AVsKyQ5jXAiXkdewIXAI8BO+XphwLrgV0K8xyVlzu6wT6ZQHrN4+Q60/bP06bl4f2A1cCHgUmkGqDrgMsL8ywEngD+PefxbaRaog8V0rwbeF3elwfl4/DrwvQpeb2LgSNzunGkbgx/Bp6fP9vk9BvyWFjGMuD0wnAAD5F+6b8Q2K2Z41ZnnxyYl3VUzsOOefws4NZCujn5+H0r74fjSa/Y/Dnw+XyMT8vLOiDPI+A3wNV5v0wkdUVYVTymdfJ0JXBhadwM4ClSAeSgvG03A/MLaY7NaWbm/HwgD7+xtN/6sm8/DCwnnTPjgZcD7yot5zRgeafjgj/t+TC44+9R+Tx7F7A7qZB+B3B+Ic0y4JG8/on5fAng0Dx9C1LM+jXwMuAVpML/U8CMHvbLLAqxozTta6T4WduOs3O+js7n2QnA34DX5+kTcp66SLF3T1IsXgOMyWm2BM4ixbIJOd3/A04urHcOKY5fBuxLiv075nP60+QYnNPOAFaX8j0l52OnQpr1pFaBV5DizbbAe4H7gWl5e94IPADMbLA/tqHxdWAZm8al+0ixdxLw/bzs+cDUnId5wJ8K8/R6jauTpx1Jsf3lpfEL87E7K6/rbXkfHF9Ic2Xellfmdc/N+3ibfti3i0nn3J6k7+ux5O9qYVm/Bz7T6djQp3jS6QxU7UMKxmuAUcDW+f/dSYWw23OaPfOX77A8vHv+8o8vLeu/yEE9nwC/LE3fIS/noDw8B7gKeAMpEL+zlP53wHdK4xZQKCDX2R7l4PKOwrhbgY8Vhn/Uy8k9gcYF5K3ztI/k4UuA2aU0tfeo75yHFwJ3kl+Vnsd9AujqIQ+1fT42D9cCwFtK6WZR50JC84W4fy+l6fW4Nbu/ynnLx3s5MKIwbhGFYFzOJ+kH0GpycCykuaV2DBrk6Rbg06VxM3I+9yiMezupYFG72P4WuKg03xzgN/2wb78G/LL4PaiT7zeRzq2RzZ7D/gzdD4M7/v4a+GRp3Jvz+ag8vAz4QSnNX4BP5P+PBJ4u5hU4LOdjRg/7pVvsKE2r/YjeGXh23meHl9J8BZiX/5+Q0/9bYfoWpJj82R7ycA6woDA8B1gJjCql63bu53EzaK4QF+TKgEK6e4ETS+P+BVjS6v4q5y2v7/OF4X3zuA/1kM9er3F11lubvltp/EJSS0hx3P8A383/T8rzvbIwfTtSofo9/bBvVwEn9XJOXgFc2tdzupOfIdkvZJC7lhSYDyUVLh+OiLskPQDsnpvkXg38Hbg+z7N/TrtEUnFZo0j96gAOAF5ZbgrJdifVZNTS/RQ4ISJ+Ukq3F/Dd0rjrSL/8gA03Anwm5/F5wAjSL+rxhXm+A/wTcI6kHUm/lI+ttzOaUNvgKOR/oqTj6qTZnVSLCPD7yGdfYTs+I+k5EbEqN2l9ihRYdiwsYzyp5qNmUR/z3Uh5ec0et75aEqkrQc2DpJoaSuN2LuTnWcDK0ndt65yfRrYB1tYZvy4i7igMryDVHm1PqkHbC7ioNM9vSAXXVpX37RzSxeBOSb8g1dT8LCKeKaRZQzr2W5MKIlZtgzn+HgAcpO5dnrYgnVvPJ1VEAPypNN8KNp6/ewH3RfeuTteTCvh9VYzBe5P2389r3RuyLUmFw6Lrav9ExDOSrs/zp4VKp5FafHYlbeOWwF9Ly7g1ItZtRt7L1pN+zNfyMJrUMvhtSd8spBvJxu3eXMXj9WD+u7jOuJ1JrRDNXuOKtsl/68Xg3r4vz9D9WD2udNPq3rSm277NvgR8V9JJpMqK/4yIP5fSrCnkf0hxAbmfRcTdkv5K+gUm0i88IuJvkm7M46eQatCeyrNtQQpOB5KayorWFNJcDdS7u/fBwv/3kE6wd0ua24fgczGpYPxBUkBcR/riF2+guJT0dIDDSM18D5P6GvVF7SSt3XiwBakQ/+U6ae9rZoFKj5WZT6odP5G0P3YC/pfu2wGp6bAZwaYBtd5NfeXlNXvc+qr8fYkG42p90rbI6z28zrJW9bCeh0k1ZmXr66yrtp7yuHrpav+3vG8j4iZJE0hNwa8hfXf/KOm1hULyjsDaiHDheBgY5PF3C1LXgXLBGVJNak1P529/FeqK9iad+4+QupdA6oZQvt+gnK+GcuHvK6T99bu8/PezaUVKs/H3GZqLEetKFQa1/XYaA3dDZnG/RA/jijG41Wvcw/nvDmz8IVVv/bX1NfN9qeWrr/uWiJgl6TLSjdBHAZ+SdFpEFCtFdmTTH1dDggvIA+NaUi2FSBftmoWkC/kU0i+vmptz2udHxLUNlnkTqX/RXwuBvZ5HSbVzvwR+KunYQpC+HTiE7jV6h5TmPwz454i4GkDS80j97zaIiEclXUHq5/syYE75xGnB6aTmngV5+CZgn4hY2st8B0tSoRb5EGBFrj0+gFQg/nhE3JO3o9nn4T5JqjUvW0lhP9TbLw00e9zKeaBBPjbXTaQfQM9Ea3dD30zrNQ6QvnOH0f07dxiwpDDc131LRDxBKnD8RNIcUn+3iaTmXkhNnjf1Id82dA3W+HsTsGcTsa0nS4AxksZFxPI87iD6eMO9pF1IfYyvyLXAS0iVIrtGxDU9z80h5Bp2par3g4DL87TDgOsjYsMzeiX11EJVVC8GrwSeVWshzONe2tuCIuJBSfcBu0fEJU2uv1Ee+kuz17iiu0g/Mvame+zszRLSd+NQUhcf8s2A+7HxiUR92rc1EfEXUjegr+Va+vfQPd7vS+pmMeT4KRYD41pS8DiYXIOR/QqYTmr+2BCII919fRkwR9I0SS+UNFnpJR61gt3XSX2HfiTp4Jzm/yg9jWDb4soj4mHgCGAscIWkUXnSV4GTJL1X6TmgZ+Q8Ft0JvEPpju0DgR+yscBW9B1Sf9OX0Pyjv54r6fmSdpN0jKS5pBsnTouNz9Y9l9QM+S1JL1O6Y/sNkr5dWtYLgK9I2kPpCQj/ysZf5PeSgvzMvJ9eT+o20oxlwK6S9pe0U2HfXQO8Px+Xl5Ga9+s1d5U1fdwKHiLVXB0l6XnKL+boJwtI/YKvzMdgN0mHSvq0pHq1yjXzgUPU+uN6zgNOVHoCxSRJHyB9b75QSNOnfSvpQ5KOl7SXpImkC/0qunehOZx006INH4M1/p4FnCDpLEn7Stozr694LvRmAemGq0skvVTpCUNfZtOWnHpG5vi7i6R9JJ1Canp/FDgj5/0J0g1q50t6d46/L1V6xvQppeW9L+d/D1Jt8a5ArRvDncD+OcZMkvRJmn+G+zLgcEljlJ+iQOpG8jfg8zlPbyF182vGLOAjSk+u2CPv+3fm619Peah3HegPzV7jNsgtYgtIPzyalguvV5K6mBwuaT/STXWrSDcUQh/3raRtJH1d6YkXEyQdTKnyQ6mFbwx9b2HurE53gq7ih9TnKSjdPU96pNpTpBrTEaVpW5JO5LtJBdIHSHebHlBIM4n0C/0xUgHqDtLdw1vl6XPId1Hn4Z1I/ZOuIt8IQQqED5H6Y34/r3NZYZ6XkE6YNaRfrSeSbsqbVcqv8vRrmtgfE/L+qH3+Tgqgs8l3iZfSTyYValaRTtzFwFmF6QtJT264gNTf9jHgi3S/We24nL+1pP6BR+V1T8nTp1C4CaEw36jCPt5w4wupQP6zvN/uAt5C/Rs2ptXZnh6PW4N99h5SQf9pYGEeN4tNb9K7qjTfVaQa/eK439P9TvltST+WukjfteWkH0K795CfETnd6wvjZtDLzR153GnAUtJ3fynw3tI8fdq3pLvTbyLdCb+KVAB6eWH6mLx9YzsdE/xp34fBHX+PJHX1+nv+zi6i8DSF8vc+j1sIXFAYflH+rq8j1dy9KZ87M3rYJ7PYGH+fzttwHfBxYNtSWpGenlGrTV5J6uv/2jx9Ql7O20ndFtbmfXFMYRlbkeL7Y6QYPRs4k+7Xmm77qzD+EOCPeblRGD+VdN1YQ/rB/g42vZFsdYPtPz7HirU5T78BpvewvxpdB7odH0pxKR/zDdeZPK52U+i+hXE9XuMa5OkoUveK4nWu23ejwfdwB1JLSu17u4BUg12cp+V9m4/x90n9yteR+j5fCDynkOYM4Oedjgl9/dTunDVriaRtSP2lPhARl7V53QtJBcWZ7VzvcKd0081bI+KITuelGZLOA7aLiHLNl5n1Ua4VvAc4MCL6+yZn64Gk60hPVrm003npTa51/wvpkXO/7XR++sJ9kK0lSg8gr93Et4b6N5xYNX0H2FHSdjE0Xjf9EKm52MysCk4l3fczFOwKnD1UC8fgArK1bjyp9qCL9FKGev2TrYIi3Yj5uU7no1kRcV6n82Bm1l8i4k9s+li3QSlS3/47e004iLmLhZmZmZlZgZ9iYWZmZmZW4AKymZmZmVmBC8hmZmZmZgUuIJuZmZmZFbiAbGZmZmZW8P8B3gn12qrOis8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nb, b = np.histogram(mean_speed['speed_kmh'], bins = np.linspace(0,30,31),density = True)\n",
    "b = b[:-1]\n",
    "frac = np.divide(nb, float(np.sum(nb)))\n",
    "\n",
    "n_dur, dur = np.histogram(durations[durations<300], bins = np.linspace(0,180,45), density = True)\n",
    "dur = dur[:-1]\n",
    "frac_dur = np.divide(n_dur, np.sum(n_dur))\n",
    "\n",
    "n_week, week = np.histogram(data_dw['hour'], bins = np.linspace(0,25,26), density = True)\n",
    "week = week[:-1]\n",
    "frac_week = np.divide(n_week, np.sum(n_week))\n",
    "\n",
    "n_weekend, weekend = np.histogram(data_dnw['hour'], bins = np.linspace(0,25,26), density = True)\n",
    "weekend = weekend[:-1]\n",
    "frac_weekend = np.divide(n_weekend, np.sum(n_weekend))\n",
    "\n",
    "\n",
    "#######\n",
    "plt.figure(figsize=(10,6))\n",
    "#00BFFF\n",
    "#6691B4 azul inicial\n",
    "plt.subplot(2,2,1)\n",
    "plt.bar(b, frac, width=1, linewidth=1, edgecolor='#000000', color='#6691B4',alpha=1, label=\"Biko\")\n",
    "plt.legend()\n",
    "xt = plt.xticks()[0]  \n",
    "xmin, xmax = min(xt), max(xt)  \n",
    "plt.xlim(0,30)\n",
    "plt.ylim(0,0.12)\n",
    "plt.xlabel(\"Speed (km/h)\",size=14)\n",
    "plt.ylabel(\"Density\",size=14)\n",
    "plt.tight_layout()\n",
    "\n",
    "######\n",
    "plt.subplot(2,2,2)\n",
    "# plt.bar(b, frac, width=1, linewidth=1, edgecolor='#000000', color='#6691B4',alpha=1, label=\"Biko\")\n",
    "plt.bar(dur, frac_dur, width=4, linewidth=1, edgecolor='#000000', color='#6691B4',alpha=1, label=\"Biko\")\n",
    "plt.legend()\n",
    "xt = plt.xticks()[0]  \n",
    "xmin, xmax = min(xt), max(xt)  \n",
    "plt.xlim(0,180)\n",
    "plt.ylim(0,0.1)\n",
    "plt.xlabel(\"Travel time (minutes)\",size=14)\n",
    "plt.ylabel(\"Density\",size=14)\n",
    "plt.tight_layout()\n",
    "\n",
    "\n",
    "######\n",
    "plt.subplot(2,2,3)\n",
    "# plt.bar(b, frac, width=1, linewidth=1, edgecolor='#000000', color='#6691B4',alpha=1, label=\"Biko\")\n",
    "plt.bar(week, frac_week, width=1, linewidth=1, edgecolor='#000000', color='#6691B4',alpha=1, label=\"Biko\")\n",
    "plt.legend()\n",
    "xt = plt.xticks()[0]  \n",
    "xmin, xmax = min(xt), max(xt)  \n",
    "plt.xlim(0,25)\n",
    "plt.ylim(0,0.15)\n",
    "plt.xlabel(\"Weekday Departure time (hours)\",size=14)\n",
    "plt.ylabel(\"Density\",size=14)\n",
    "plt.tight_layout()\n",
    "\n",
    "######\n",
    "plt.subplot(2,2,4)\n",
    "# plt.bar(b, frac, width=1, linewidth=1, edgecolor='#000000', color='#6691B4',alpha=1, label=\"Biko\")\n",
    "plt.bar(weekend, frac_weekend, width=1, linewidth=1, edgecolor='#000000', color='#6691B4',alpha=1, label=\"Biko\")\n",
    "plt.legend()\n",
    "xt = plt.xticks()[0]  \n",
    "xmin, xmax = min(xt), max(xt)  \n",
    "plt.xlim(0,25)\n",
    "plt.ylim(0,0.15)\n",
    "plt.xlabel(\"Weekend Departure time (hours)\",size=14)\n",
    "plt.ylabel(\"Density\",size=14)\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.savefig('./img/descriptive.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fourth Part: Map Matching"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_all=data[['date','activity_id','new_activity_id','user_id','altitude','stays','x1','y1']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1419217, 8), (1419217, 25))"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Delete points where the person is not moving\n",
    "data_all=data_all[data_all['stays']>0]\n",
    "data_all.shape,data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4179"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(data_all.new_activity_id.unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare the coordinates for the link\n",
    "data_all.loc[:,'coordinates_str']=(data_all['x1']).astype(str)+','+(data_all['y1']).astype(str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_all.reset_index(inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate the delta time "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_all.loc[:,'date_1']=data_all['date'].apply(lambda x: datetime.strptime(x,'%Y-%m-%dT%H:%M:%S'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sort the data to be in chronologycal order from the oldest to the newest date\n",
    "data.sort_values(['activity_id','date_1'],inplace=True)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to calculate the delta seconds from the start of an activity id to each point recoded\n",
    "def total_second(arr):\n",
    "    return (arr-arr.min())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the total_seconds function to create a list\n",
    "cum= data[['date_1']].groupby(data['activity_id']).transform(total_second)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "cum.rename(columns={'date_1':'delta_seconds'},inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add a column where we see the seconds from the start of the trip\n",
    "merge=pd.merge(data_all,cum,left_on='index',right_index=True,how='left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check that no duplicates were created\n",
    "merge.shape,data_all.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "merge.loc[:,'timestamp']=merge['delta_seconds'].apply(lambda x: str(int(x.total_seconds())))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_all=merge.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract coordinates "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_coordinates (data):\n",
    "    list_coordinates=data['coordinates_str'].tolist()\n",
    "    return ';'.join(list_coordinates)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_times (data):\n",
    "    list_times=data['timestamp'].tolist()\n",
    "    return ';'.join(list_times)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### API Query\n",
    "Note: The online OSRM service allows map matching for the auto mode only. To allow other modes, a local version of the service has to be set. We followed this post to set our local version: https://github.com/Project-OSRM/osrm-backend/wiki/Running-OSRM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "metadata": {},
   "outputs": [],
   "source": [
    "def api_query(profile,coordinates,times):\n",
    "    url_base='http://localhost:5000/match/v1/' # link to make local queries \n",
    "    url=url_base+profile+'/'+coordinates+'?steps=false&geometries=polyline&overview=full&annotations=true&tidy=true&timestamps='+times\n",
    "    # Make the query\n",
    "    response = requests.get(url)\n",
    "    result=response.text\n",
    "    dict_result=json.loads(result)\n",
    "    return dict_result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 632,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_result(result,new_activity_id,i):\n",
    "    list_df=[]\n",
    "    tracepoints=pd.DataFrame(index=[list(range(len(result['tracepoints'])))],columns=[])\n",
    "    tracepoints.loc[:,'points']=result['tracepoints']\n",
    "    tracepoints.loc[:,'new_activity_id']=new_activity_id\n",
    "    tracepoints.loc[:,'repeat_num']=i\n",
    "    for match in range(len(result['matchings'])):\n",
    "        frame=pd.DataFrame(result['matchings'][match]['legs'])\n",
    "        frame.loc[:,'match']=match\n",
    "        frame.loc[:,'repeat_num']=i\n",
    "        frame.loc[:,'new_activity_id']=new_activity_id\n",
    "        frame.loc[:,'confidence']=result['matchings'][match]['confidence']\n",
    "        frame.loc[:,'code']=result['code']\n",
    "        list_df.append(frame)\n",
    "    dfs=pd.concat(list_df)\n",
    "    return dfs,tracepoints\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table with each activity and node sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 528,
   "metadata": {},
   "outputs": [],
   "source": [
    "profile='bicycle'\n",
    "map_match=pd.DataFrame(index=[],columns=['new_activity_id','legs','node_sequence','error'])\n",
    "idx=0\n",
    "# Limit is a variable with a function to stop querying from time to time, this was necessary when doing online queries.\n",
    "limit=100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 531,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "idx=0\n",
    "dfs=[]\n",
    "tracepoints=[]\n",
    "for activity in data_all['new_activity_id'].unique():\n",
    "    # the program cannot handle more than a 100 points at the same time\n",
    "    repeat_number=np.ceil(len(data_all[data_all['new_activity_id']== activity])/100)\n",
    "    for i in range (1,int(repeat_number)+1):\n",
    "        # We divide the data in subsets of a 100 points so that the program does not crash\n",
    "        subset=data_all[data_all['new_activity_id']==activity].reset_index()[(i-1)*100:i*100]\n",
    "        list_coordinates=extract_coordinates(subset) \n",
    "        list_times=extract_times(subset) \n",
    "        result=api_query(profile,list_coordinates,list_times)        \n",
    "        if result['code']!='Ok':\n",
    "            map_match.loc[idx,'new_activity_id']=activity\n",
    "            map_match.loc[idx,'error']=1\n",
    "            map_match.loc[idx,'node_sequence']=result['code']\n",
    "            idx=idx+1\n",
    "        else:            \n",
    "            map_match=map_match.astype('object')\n",
    "            map_match.loc[idx,'new_activity_id']=activity\n",
    "            map_match.loc[idx,'error']=0\n",
    "            result_read=read_result(result,activity,i)\n",
    "            dfs.append(result_read[0])\n",
    "            tracepoints.append(result_read[1])\n",
    "            idx=idx+1\n",
    "            if idx==limit:\n",
    "                time.sleep(5)\n",
    "                limit=limit+100\n",
    "# Save the data after running\n",
    "tracepoints=pd.read_csv('tracepoints.csv')\n",
    "new=pd.concat(dfs)\n",
    "new.to_csv('matched.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reopen data\n",
    "new=pd.read_csv('matched.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare the unique activity id on the `new` dataframe so that we can know how many activities the OSRM was not able to locate\n",
    "len(new['new_activity_id'].unique())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filter data and create the pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "new.reset_index(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "new.loc[:,'nodes']=new['annotation'].apply(lambda x: ast.literal_eval(x)['nodes'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "new.loc[:,'pairs']=new['nodes'].apply(lambda x: list(zip(x[0:-1],x[1:])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "def unique_nodes(data):\n",
    "    list_list=data['pairs'].tolist()\n",
    "    flat_list=[y for x in list_list for y in x]\n",
    "    list_pairs=pd.Series(flat_list).unique()\n",
    "    return list_pairs\n",
    "# you can try both options\n",
    "# uniques=[list_nodes[0]]\n",
    "# for i in range(1,len(list_nodes)):\n",
    "#     if not list_nodes[i]==list_nodes[i-1]:\n",
    "#         uniques.append(list_nodes[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Need to keep with unique nodes since the OSRM may have repeated sequence of nodes from a trace to the other. \n",
    "df=pd.DataFrame(data=new.groupby(['new_activity_id','confidence'],sort=False).apply(unique_nodes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.reset_index(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Edge sequence is the sequence of nodes a person is travelling through translated into an edge for future network analysis\n",
    "df.rename(columns={0:'edges_sequence'},inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "list_edges=df['edges_sequence'].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [],
   "source": [
    "flat_list=[y for x in list_edges for y in x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a network to get the weight for each edge. \n",
    "g = nx.DiGraph((x, y, {'weight': v}) for (x, y), v in Counter(flat_list).items())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "840346.0"
      ]
     },
     "execution_count": 127,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g.size(weight='weight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Translate network to a df\n",
    "pairs=pd.DataFrame([[x, y, v] for (x, y), v in Counter(flat_list).items()])\n",
    "pairs.columns=['source','target','weight']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "840346"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pairs['weight'].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [],
   "source": [
    "pairs.set_index('source',inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save pairs to csv\n",
    "pairs.to_csv('pairs.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Network edges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download the network to check how many people travelled through its links\n",
    "west=-74.3068778-0.01\n",
    "east=-73.9337932+0.01\n",
    "north=4.924683900000001+0.01\n",
    "south=4.4266025-0.01\n",
    "\n",
    "# create network from that bounding box\n",
    "G1 = ox.graph_from_bbox(north, south, east, west, network_type='all_private',simplify=False)\n",
    "df_g1=ox.graph_to_gdfs(G1, nodes=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [],
   "source": [
    "nodes_bog=pd.DataFrame(index=[],columns=['node_id'])\n",
    "nodes_bog.loc[:,'node_id']=list(G1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [],
   "source": [
    "nodes_bog.loc[:,'lat']=nodes_bog['node_id'].apply(lambda x: G1.node[x]['y'])\n",
    "nodes_bog.loc[:,'lon']=nodes_bog['node_id'].apply(lambda x: G1.node[x]['x'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "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>node_id</th>\n",
       "      <th>lat</th>\n",
       "      <th>lon</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1111282148</td>\n",
       "      <td>512157.519138</td>\n",
       "      <td>604398.927445</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1111282149</td>\n",
       "      <td>512171.151775</td>\n",
       "      <td>604416.746714</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      node_id            lat            lon\n",
       "0  1111282148  512157.519138  604398.927445\n",
       "1  1111282149  512171.151775  604416.746714"
      ]
     },
     "execution_count": 136,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Nodes bog have the nodes id and coordinate for a certain bbox\n",
    "nodes_bog.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [],
   "source": [
    "nodes_bog.rename(columns={'node_id':'id'}).set_index('id').to_csv('nodes_network_bike_.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [],
   "source": [
    "lista=[]\n",
    "for line in nx.generate_edgelist(G1,data=False):\n",
    "    lista.append(line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [],
   "source": [
    "df=pd.DataFrame(index=[list(range(len(lista)))],columns=['lista','source','target','pairs'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[:,'lista']=lista"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[:,'source']=df['lista'].apply(lambda x: str(x.split()[0]))\n",
    "df.loc[:,'target']=df['lista'].apply(lambda x: str(x.split()[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [],
   "source": [
    "# invert origin and destination to represent two ways edges\n",
    "df.loc[:,'pairs']=df.apply(lambda row: ((row['source']),(row['target'])),axis=1)\n",
    "df.loc[:,'pairs_2']=df.apply(lambda row: ((row['target']),(row['source'])),axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[:,'weight']=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.drop('lista',axis=1,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {},
   "outputs": [],
   "source": [
    "pairs.reset_index(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [],
   "source": [
    "pairs.loc[:,'source']=pairs['source'].astype('str')\n",
    "pairs.loc[:,'target']=pairs['target'].astype('str')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [],
   "source": [
    "pairs.loc[:,'pairs']=pairs.apply(lambda row: ((row['source']),(row['target'])),axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [],
   "source": [
    "merge=pd.merge(df,pairs[['pairs','weight']],left_on='pairs',right_on='pairs',how='left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [],
   "source": [
    "pairs_notfound=pairs[~pairs['pairs'].isin(merge[merge['weight'].notnull()]['pairs'].tolist())]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [],
   "source": [
    "merge.drop(['lista'],axis=1,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "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>source</th>\n",
       "      <th>target</th>\n",
       "      <th>pairs</th>\n",
       "      <th>pairs_2</th>\n",
       "      <th>weight</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1111282148</td>\n",
       "      <td>1927389468</td>\n",
       "      <td>(1111282148, 1927389468)</td>\n",
       "      <td>(1927389468, 1111282148)</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1111282148</td>\n",
       "      <td>5643183061</td>\n",
       "      <td>(1111282148, 5643183061)</td>\n",
       "      <td>(5643183061, 1111282148)</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       source      target                     pairs                   pairs_2  \\\n",
       "0  1111282148  1927389468  (1111282148, 1927389468)  (1927389468, 1111282148)   \n",
       "1  1111282148  5643183061  (1111282148, 5643183061)  (5643183061, 1111282148)   \n",
       "\n",
       "   weight  \n",
       "0     NaN  \n",
       "1     NaN  "
      ]
     },
     "execution_count": 156,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "merge.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [],
   "source": [
    "swap_od=pd.merge(merge,pairs_notfound[['pairs','weight']],left_on='pairs_2',right_on='pairs',how='left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(126758, 4)"
      ]
     },
     "execution_count": 158,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pairs_notfound2=pairs[(~pairs['pairs'].isin(df['pairs'].tolist()))&(~pairs['pairs'].isin(df['pairs_2'].tolist()))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [],
   "source": [
    "swap_od['weight_x'].fillna(0,inplace=True)\n",
    "swap_od['weight_y'].fillna(0,inplace=True)\n",
    "swap_od['weight_raw']=swap_od['weight_x']+swap_od['weight_y']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "297515.0"
      ]
     },
     "execution_count": 160,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This is the amount of links we have, from now on, we are going to add one weight to be able to draw the entire network\n",
    "swap_od['weight_raw'].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "297515"
      ]
     },
     "execution_count": 161,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check we had assigned all the links to the network or the not_found group\n",
    "pairs['weight'].sum()-pairs_notfound2['weight'].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "375138.0"
      ]
     },
     "execution_count": 162,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "swap_od['weight']=swap_od['weight_raw']+1\n",
    "swap_od['weight'].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "swap_od[['source','target','weight']].set_index('source').to_csv('network_with_flow.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}