diff --git a/.gitignore b/.gitignore index 5b698d8..72b5a03 100644 --- a/.gitignore +++ b/.gitignore @@ -17,6 +17,11 @@ secret.py *.sh !/data/more_gtfs_sources.csv !/data/cleaned_sources.csv +/notebooks/cache/* +/notebooks/checkpoints/* +/notebooks/files/* +/notebooks/lightning_logs/* +*/*.zip # Byte-compiled / optimized / DLL files __pycache__/ diff --git a/README.md b/README.md index f0f64d9..4a7aa3c 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,35 @@ -# open_bus_tools +# Open Bus Tools Tools for open bus data collection, travel time prediction, energy consumption and visualization. -If not installed as package, add [export PYTHONPATH="${PYTHONPATH}:/my/other/path"] to .bashrc \ No newline at end of file +This is companion code for a dissertation in Civil and Environmental Engineering at the University of Washington. It will be distributed as a Python package when completed. + +If not installed as package, you may need to add [export PYTHONPATH="${PYTHONPATH}:/my/other/path/to/open_bus_tools"] to .bashrc to run the scripts. + +### Workflow +There are two folders in the base directory /scripts and /notebooks. These are useful starting points for understanding how different parts of the code work together. The two main modules of the package are: +* drivecycle: Used to estimate energy use based on vehicle and driving profile from a given actual or predicted GPS trajectory. + +* traveltime: Used to predict travel times using ML-based geospatial data mining. + +The /webscraper folder has code used to gather GTFS and GTFS-RT feeds. These data are used in training the travel time models. + +### Motivation +Standardized and open source bus data including static schedules and realtime positions have become widely available in public web APIs. Though these data primarily underly popular map-based trip planning applications, they also enable new analyses in understanding, forecasting and improving bus operations across agencies and cities. Due to their lower resolution and simpler, anonymized features, these data are more challenging to work with than those of the underlying sensors, making them unfavorable for direct operations analysis by transit agencies. However, the wide scale of these data and their open source nature make them a valuable resource for researchers and planners. + +This work develops and tests a set of tools for analyzing bus operations using open source data. Central to this endeavor is the ongoing collection of a multi-year dataset from the King County Metro transit network in Seattle, Washington, rapidly approaching one billion tracked bus locations. First, these data are used in basic roadway segment matching and aggregation to visualize the spatiotemporal dynamics of different types of delays in the transit system. The findings are used to identify locations in the network where targeted transit priority treatments can generate the greatest benefit. Then a set of deep learning models and baselines are developed to forecast bus travel times under different data availability scenarios. Their generalizability is tested across different cities and transit networks. Finally, these models are used to estimate the drive cycle energy demands of a system-wide bus electrification project, and are validated using results of a more traditional approach with direct sensor data. + +### Past Work +This work will replace code in these repositories. +* [ML-Based Bus Travel Time Models](https://github.com/zackAemmer/generalizable_travel_times) +* [Bus Performance Visualization Tool](https://github.com/zackAemmer/transit_vis) +* [TransitVis Web-App](https://www.transitvis.com) + +### Data Sample +The full dataset is stored on Amazon S3 and downloaded for local analysis as needed, please feel free to reach out if you would like a sample. + +### Other Libraries and Code You Might Find Useful +* https://github.com/kraina-ai/srai +* https://github.com/NREL/fastsim +* https://github.com/smohiudd/drivecycle +* https://github.com/EricaEgg/Route_Dynamics +* https://github.com/mrcagney/gtfs_kit \ No newline at end of file diff --git a/data/cleaned_sources.csv b/data/cleaned_sources.csv index 9333983..25cb99e 100644 --- a/data/cleaned_sources.csv +++ b/data/cleaned_sources.csv @@ -17,7 +17,6 @@ US,Georgia,Atlanta,Metropolitan Atlanta Rapid Transit Authority (MARTA),https:// US,Illinois,Springfield,Springfield Mass Transit District (SMTD),http://data.smtd.org/gtfs/smtd_gtfs_feed.zip,http://ride.smtd.org/gtfsrt/vehicles,-89.75837,-89.527032,39.667827,39.898621,dba09c2c-6433-4924-acdd-61745c439460,America/Chicago,32616 US,Wisconsin,Madison,Metro Transit,http://transitdata.cityofmadison.com/GTFS/mmt_gtfs.zip,http://transitdata.cityofmadison.com/Vehicle/VehiclePositions.pb,-89.563863,-89.244818,42.987718,43.176537,38a8b09e-078a-491f-b4bf-4b001ee51a61,America/Chicago,32616 US,Ohio,Cleveland,Greater Cleveland Regional Transit Authority,http://www.riderta.com/sites/default/files/gtfs/latest/google_transit.zip,http://gtfs.gcrta.org/TMGTFSRealTimeWebService/Vehicle/VehiclePositions.pb,-81.96736,-81.437811,41.227828,41.637051,1656b240-d1ee-45ac-afc4-0bf17f4b4466,America/New_York,32617 -US,Pennsylvania,Pittsburgh,Port Authority of Allegheny County,https://www.portauthority.org/developerresources/GTFS.zip,http://truetime.portauthority.org/gtfsrt-train/vehicles,-80.258865,-79.70545,40.273078,40.667597,2f18ef13-fd56-4e71-94ce-e96973a0175c,America/New_York,32617 US,Pennsylvania,Pittsburgh,Port Authority of Allegheny County,https://www.portauthority.org/developerresources/GTFS.zip,https://truetime.portauthority.org/gtfsrt-bus/vehicles,-80.258865,-79.70545,40.273078,40.667597,ead6fd88-9158-44c2-b71f-cd093a358072,America/New_York,32617 US,Massachusetts,Boston,Massachusetts Bay Transportation Authority (MBTA),https://cdn.mbta.com/MBTA_GTFS.zip,https://cdn.mbta.com/realtime/VehiclePositions.pb,-71.848488,-70.276583,41.581289,42.797837,1287a85e-ca1c-47d3-b51f-d58fa9bd6784,America/New_York,32619 US,Virginia,Arlington,Arlington Transit,https://www.arlingtontransit.com/shared/content/gtfs/art/google_transit.zip,https://realtime.arlingtontransit.com/gtfsrt/vehicles,-77.162277,-77.049105,38.839365,38.925147,1161e355-c52c-4925-a92b-1eabc12d0153,America/New_York,32618 @@ -44,7 +43,6 @@ CA,Québec,La Vallée-du-Richelieu,Exo Vallée du Richelieu,https://exo.quebec/x CA,Québec,Sainte-Julie,Exo Ville de Sainte-Julie,https://exo.quebec/xdata/omitsju/google_transit.zip,https://opendata.exo.quebec/ServiceGTFSR/VehiclePosition.pb?agency=OMITSJU,-73.56448,-73.288934,45.497418,45.63169,8f392db0-4b1b-4b7f-8036-4fb7504bb196,America/Toronto,32618 CA,Québec,L'Assomption,Exo L'Assomption,https://exo.quebec/xdata/mrclasso/google_transit.zip,https://opendata.exo.quebec/ServiceGTFSR/VehiclePosition.pb?agency=MRCLASSO,-73.539828,-73.277404,45.589483,45.901438,9523f994-7c76-4cfc-ba9b-2b5751f90f28,America/Toronto,32618 CA,British Columbia,Port Alberni,Port Alberni,https://bct.tmix.se/Tmix.Cap.TdExport.WebApi/gtfs/?operatorIds=11,https://bct.tmix.se/gtfs-realtime/vehicleupdates.pb?operatorIds=11,-124.85047,-124.78265,49.21847,49.27743,16c4e44f-108d-4e04-80ab-b305093f1328,America/Vancouver,32610 -FI,Keski-Suomi,Jyväskylä,"Jyväskylän Liikenne Oy, Koivuranta Oy, Uuraisten Liikenne Ky, Savo-Karjalan Linja Oy, Linja-Karjala Oy, Oulun taksipalvelut Oy, Oubus Oy, Suorsan Liikenne Ky, Koskilinjat Oy, Oy Pohjolan liikenne Ab, Kuopion Liikenne Oy, Autolinjat Oy, ML-Charter Oy, Etelä-Suomen Linjaliikenne Oy, Linjaliikenne Martti Laurila Oy, Kari Väisänen Ky, Matkatoimisto Matka-Majuri Ky, Linja-autoliikenne P. Puolakka Ky, Elimäen Liikenne Oy, Anjalankosken Linja Oy, Savonlinja Oy, Liikenne M. Heikura Oy, Jääskeläisen Auto Oy, Kilpailutus käynnissä, TILAUSAJOT MENNÄÄN BUSSILLA OY, Revon Turistiliikenne Oy, Matkamestarit Oy, R-Kioski, A. Valppu Oy, Invataksi Niemi Oy, Oy Wiik & Ström Ab, Oravaisten Liikenne Oy, Vaasan Paikallisliikenne Oy, Koiviston Auto Oy, Lehtimäen Liikenne Oy, Järvisen Liikenne Oy, Niemisen Linjat Oy, Bus Travel Oy Reissu Ruoti, Taksi- ja linja-autoliikenne Timo Lahtinen, SL-Autoyhtymä Oy, Linja-autoliike E Ahonen Ky, Lehdon Liikenne Oy, HämeBus Oy, Tuuri Jaakko Juhani, Pekolan Liikenne Oy, Mikkolan Liikenne Oy, Vekka Liikenne Oy, Inter Kuljetus Oy, Valkeakosken Liikenne Oy, Matkahuolto Oy, ELY, Ventoniemi Oy, Yhdysliikenne Oy, MIKA K. NISKANEN OY, Kantamatkat, LINJA-AUTOLIIKE YRJÖ MAKKONEN & KUMPP, SOISALON LIIKENNE OY, Pohjolan Matka, HENKILÖKULJETUS RISSANEN OY, V.Alamäki Oy, Kai Heikkinen Oy, Savo Kainuun Tilausmatkat Ky, Kajaanin Tila-autot Oy, A Kyllönen Oy, Puolangan Bussipalvelu Ky, Kainuun Tilausliikenne P. Jääskeläinen Ky, Matka-Kyllönen Oy, Jyrkilä Oy, Liikenne Vuorela Oy, Pyhtään Kuljetuspalvelut Oy, Erkki Itkonen Oy, Ihastjärven Linja Oy, Tilausliikenne Himanen, Mikkelin Palvelutaksit Oy, Tuplataxi, Matka Mäkelä Oy, Tilausliikenne Hokkinen Ky, S. J. TOIVONEN OY, E. RANTANEN OY, Kuljetus Mikkonen KY, Linja-autoliike Veljekset Laitinen Oy, TILAUSLIIKENNE HÄNNINEN OY, JOUTSENON AUTO- JA KULJETUSPALVELU OY, TOIMI VENTO KY, NINA´S ERIKOISTAKSIPALVELUT OY, Tarkiainen Matti Olavi, Anssin Tilausliikenne Oy, Rautalammin Auto Oy, Oy Rytkönen & Co, SAIMAAN TURISTILIIKENNE OY, Studio Ticket Pusatec-DO-NOT-USE, Studio Ti",https://tvv.fra1.digitaloceanspaces.com/209.zip,https://data.waltti.fi/jyvaskyla/api/gtfsrealtime/v1.0/feed/vehicleposition,25.3666653747545,26.4944285887841,61.7425528275114,62.62864111217454,cd173e54-d152-4aaf-9c45-ede279ffbc98,Europe/Helsinki,32635 US,California,Commerce,Commerce Municipal Bus Lines,https://citycommbus.com/gtfs,https://citycommbus.com/gtfs-rt/vehiclepositions,-118.264252,-118.127642,33.965568,34.063075,97b3ab44-c33b-4ac8-9468-cfe582f201ca,America/Los_Angeles,32611 FI,Päijät-Häme,Lahti,"Haarasilta Toivo Samuli, Järvisen Liikenne Oy, Koiviston Auto Oy, Lehtimäen Liikenne Oy, Bus Travel Oy Reissu Ruoti, Tilausliikenne Kuisma Ky",https://tvv.fra1.digitaloceanspaces.com/223.zip,https://data.waltti.fi/lahti/api/gtfsrealtime/v1.0/feed/vehicleposition,25.2492841801289,26.2609529033383,60.6772517226169,61.5800619928604,8ccc2bab-ae57-44b9-b989-e6dbae4d41ed,Europe/Helsinki,32635 FI,Pohjois-Karjala,Joensuu,"Linja-Karjala Oy, Savo-Karjalan Linja Oy",https://tvv.fra1.digitaloceanspaces.com/207.zip,https://data.waltti.fi/joensuu/api/gtfsrealtime/v1.0/feed/vehicleposition,28.924184940398234,30.936747,62.092129,63.325670884853885,ff767665-d709-43a5-a91a-f301c076b7f8,Europe/Helsinki,32635 diff --git a/notebooks/explore_inference.ipynb b/notebooks/explore_inference.ipynb index d28bbed..ec6088b 100755 --- a/notebooks/explore_inference.ipynb +++ b/notebooks/explore_inference.ipynb @@ -2,31 +2,67 @@ "cells": [ { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "import json\n", - "import os\n", + "import pickle\n", "import sys\n", + "from zoneinfo import ZoneInfo\n", "sys.path.append(\"../\")\n", "\n", - "import matplotlib.pyplot as plt\n", + "from dotenv import load_dotenv\n", + "load_dotenv()\n", "import geopandas as gpd\n", + "import importlib\n", + "import contextily as cx\n", + "import logging\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import os\n", "import pandas as pd\n", - "import plotly.express as px\n", + "from pathlib import Path\n", "import plotly.figure_factory as ff\n", - "import pyproj\n", - "import seaborn as sns\n", - "import lightning.pytorch as pl\n", + "import plotly.express as px\n", "import numpy as np\n", - "from sklearn import metrics\n", - "from torch.utils.data import DataLoader, SequentialSampler\n", - "from dotenv import load_dotenv\n", - "load_dotenv()\n", + "px.set_mapbox_access_token(os.environ[\"MAPBOX_TOKEN\"])\n", + "import plotly.express as px\n", + "import lightning.pytorch as pl\n", + "import rasterio as rio\n", + "from rasterio.plot import show\n", + "import seaborn as sns\n", + "import shapely\n", + "import statsmodels.api as sm\n", + "import torch\n", + "from torch.utils.data import DataLoader\n", + "\n", + "from openbustools import plotting, spatial, standardfeeds\n", + "from openbustools.traveltime import data_loader, model_utils\n", + "from openbustools.drivecycle import trajectory\n", + "from openbustools.drivecycle.physics import conditions, energy, vehicle" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "if torch.cuda.is_available():\n", + " num_workers=4\n", + " pin_memory=True\n", + " accelerator=\"cuda\"\n", + "else:\n", + " num_workers=0\n", + " pin_memory=False\n", + " accelerator=\"cpu\"\n", + "\n", + "logging.getLogger(\"lightning\").setLevel(logging.ERROR)\n", + "logging.getLogger(\"pytorch_lightning\").setLevel(logging.ERROR)\n", + "logging.getLogger(\"lightning.pytorch.accelerators.cuda\").setLevel(logging.ERROR)\n", "\n", - "from openbustools.traveltime import data_loader, grids\n", - "from openbustools import data_utils" + "model = model_utils.load_model(\"../logs/saved_models/\", \"kcm_1month_resampled\", \"GRU\", 0)\n", + "model.eval()" ] }, { @@ -35,9 +71,17 @@ "metadata": {}, "outputs": [], "source": [ - "run_folder = \"../results/big_run/\"\n", - "atb_network_folder = \"atb/\"\n", - "kcm_network_folder = \"kcm/\"" + "train_data_folders = [f\"../data/kcm_realtime/processed/\"]\n", + "test_data_folders = [f\"../data/atb_realtime/processed/\"]\n", + "train_dates = standardfeeds.get_date_list('2023_03_15', 30)\n", + "test_dates = standardfeeds.get_date_list('2023_04_15', 7)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Inference on Holdout Routes" ] }, { @@ -46,71 +90,58 @@ "metadata": {}, "outputs": [], "source": [ - "with open(f\"{run_folder}{atb_network_folder}deeptte_formatted/train_summary_config.json\", \"r\") as f:\n", - " atb_config = json.load(f)\n", - "atb_crs = pyproj.crs.CRS.from_epsg(atb_config['epsg'][0])\n", - "\n", - "with open(f\"{run_folder}{kcm_network_folder}deeptte_formatted/train_summary_config.json\", \"r\") as f:\n", - " kcm_config = json.load(f)\n", - "kcm_crs = pyproj.crs.CRS.from_epsg(kcm_config['epsg'][0])\n", + "holdout_data, holdout_routes, holdout_config = data_loader.load_h5(train_data_folders, test_dates, only_holdout=True, holdout_routes=model.holdout_routes, config=model.config)\n", + "holdout_dataset = data_loader.H5Dataset(holdout_data)\n", + "holdout_dataset.include_grid = model.include_grid\n", "\n", - "px.set_mapbox_access_token(os.getenv(\"MAPBOX_TOKEN\"))\n", - "default_crs = pyproj.crs.CRS.from_epsg(4326)\n", - "\n", - "fold_num = 4\n", - "grid_s_size = 500\n", - "model_type=\"GRU\"\n", - "skip_gtfs=False\n", - "num_workers=0\n", - "pin_memory=False\n", + "test_data, holdout_routes, test_config = data_loader.load_h5(train_data_folders, test_dates, only_holdout=False, holdout_routes=model.holdout_routes, config=model.config)\n", + "test_dataset = data_loader.H5Dataset(test_data)\n", + "test_dataset.include_grid = model.include_grid" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "gtfs = standardfeeds.get_gtfs_shapes_lookup(\"../data/kcm_gtfs/2023_01_23/\")\n", + "gtfs_shapes = standardfeeds.get_gtfs_shapes(\"../data/kcm_gtfs/2023_01_23/\", epsg=32148)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# # Read the pickle files, splitting holdout and non-holdout samples\n", + "# res = {}\n", + "# all_holdout_shingles = []\n", + "# all_other_shingles = []\n", + "# for day in train_dates:\n", + "# print(day)\n", + "# shingles = pd.read_pickle(Path('..', 'data', 'kcm_realtime', 'processed', day))\n", + "# holdout_shingles = shingles[shingles['route_id'].isin(model.holdout_routes)].sample(10)\n", + "# other_shingles = shingles[~shingles['route_id'].isin(model.holdout_routes)].sample(10)\n", + "# all_holdout_shingles.append(holdout_shingles)\n", + "# all_other_shingles.append(other_shingles)\n", + "# all_holdout_shingles = pd.concat(all_holdout_shingles)\n", + "# all_other_shingles = pd.concat(all_other_shingles)\n", "\n", - "# Define embedded variables for network models\n", - "embed_dict = {\n", - " 'timeID': {\n", - " 'vocab_size': 1440,\n", - " 'embed_dims': 27\n", - " },\n", - " 'weekID': {\n", - " 'vocab_size': 7,\n", - " 'embed_dims': 4\n", - " }\n", - "}\n", - "hyperparameter_dict = {\n", - " 'FF': {\n", - " 'batch_size': 512,\n", - " 'hidden_size': 128,\n", - " 'num_layers': 2,\n", - " 'dropout_rate': .2\n", - " },\n", - " 'CONV': {\n", - " 'batch_size': 512,\n", - " 'hidden_size': 64,\n", - " 'num_layers': 3,\n", - " 'dropout_rate': .1\n", - " },\n", - " 'GRU': {\n", - " 'batch_size': 512,\n", - " 'hidden_size': 64,\n", - " 'num_layers': 2,\n", - " 'dropout_rate': .05\n", - " },\n", - " 'TRSF': {\n", - " 'batch_size': 512,\n", - " 'hidden_size': 64,\n", - " 'num_layers': 3,\n", - " 'dropout_rate': .1\n", - " },\n", - " 'DEEPTTE': {\n", - " 'batch_size': 512\n", - " }\n", - "}" + "# print(pd.unique(all_holdout_shingles['route_id']))\n", + "# print([x in pd.unique(all_other_shingles['route_id']) for x in pd.unique(all_holdout_shingles['route_id'])])\n", + "# print(pd.unique(all_other_shingles['route_id']))\n", + "# holdout_routes_in_data = pd.unique(all_holdout_shingles['route_id'])" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "### Inference on Sample of Shingles" + "# all_other_shingles" ] }, { @@ -119,45 +150,15 @@ "metadata": {}, "outputs": [], "source": [ - "# Set up model\n", - "atb_base_model_list, atb_nn_model = model_utils.make_one_model(model_type, hyperparameter_dict=hyperparameter_dict, embed_dict=embed_dict, config=atb_config, skip_gtfs=skip_gtfs, load_weights=True, weight_folder=f\"{run_folder}{kcm_network_folder}models/{model_type}/logs/{model_type}/version_{fold_num}/checkpoints/\", fold_num=4)\n", - "print(f\"Evaluating: {atb_nn_model.model_name}\")\n", - "# Set up dataset\n", - "atb_dataset = data_loader.LoadSliceDataset(f\"{run_folder}{atb_network_folder}deeptte_formatted/test\", atb_config, skip_gtfs=skip_gtfs)\n", - "atb_ngrid = grids.NGridBetter(atb_config['grid_bounds'][0], grid_s_size)\n", - "atb_ngrid.add_grid_content(atb_dataset.get_all_samples(keep_cols=['shingle_id','locationtime','x','y','speed_m_s','bearing']), trace_format=True)\n", - "atb_ngrid.build_cell_lookup()\n", - "atb_dataset.grid = atb_ngrid\n", - "atb_dataset.add_grid_features = atb_nn_model.requires_grid\n", - "loader = DataLoader(atb_dataset, sampler=SequentialSampler(atb_dataset), collate_fn=atb_nn_model.collate_fn, batch_size=atb_nn_model.batch_size, pin_memory=pin_memory, num_workers=num_workers, drop_last=False)\n", - "trainer = pl.Trainer(\n", - " accelerator=\"cpu\",\n", - " logger=False\n", - ")\n", - "preds_and_labels = trainer.predict(model=atb_nn_model, dataloaders=loader)\n", + "# # Plot the holdout shapes over heatmap of all other shingles\n", + "# df = all_other_shingles\n", + "# holdout_shapes = gtfs_shapes[gtfs_shapes['route_id'].isin(holdout_routes_in_data)].groupby('route_id').nth(0)\n", "\n", - "# Extract predictions\n", - "atb_preds = np.concatenate([x['out'][x['mask']] for x in preds_and_labels])\n", - "atb_labels = np.concatenate([x['y'][x['mask']] for x in preds_and_labels])\n", - "# Extract data points and connect to predictions\n", - "atb_data = [x for x in atb_dataset]\n", - "atb_feats = np.concatenate([x['samp'] for x in atb_data])\n", - "atb_ys = np.array([x['norm_label'] for x in atb_data])\n", - "atb_res = pd.DataFrame(atb_feats, columns=atb_dataset.col_names)\n", - "atb_res['preds'] = atb_preds\n", - "atb_res['labels'] = atb_labels\n", - "# Get geometries and other features for every prediciton point\n", - "atb_res['pred_speeds'] = atb_res['dist_calc_km']*1000 / atb_res['preds']\n", - "atb_res['label_speeds'] = atb_res['dist_calc_km']*1000 / atb_res['labels']\n", - "atb_res['absolute_error'] = abs(atb_res['preds'] - atb_res['labels'])\n", - "atb_res['hour'] = atb_res['timeID']//60\n", - "points = gpd.points_from_xy(atb_res['lon'], atb_res['lat'], crs=\"EPSG:4326\")\n", - "atb_res = gpd.GeoDataFrame(atb_res, geometry=points)\n", - "atb_res = atb_res.sample(100000)\n", - "# Overall accuracy on data points\n", - "print(f\"MAE: {metrics.mean_absolute_error(atb_labels,atb_preds)}\")\n", - "print(f\"MAPE: {metrics.mean_absolute_percentage_error(atb_labels,atb_preds)}\")\n", - "atb_res" + "# fig, axes = plt.subplots(1, 1, figsize=(10,10))\n", + "# sns.kdeplot(ax=axes, x=df.x, y=df.y, cmap=\"Reds\", fill=True, bw_adjust=.4)\n", + "# holdout_shapes.plot(ax=axes, column='route_id', linewidth=3)\n", + "# cx.add_basemap(ax=axes, crs=holdout_shapes.crs.to_string(), alpha=0.3, source=cx.providers.MapBox(accessToken=os.getenv(key=\"MAPBOX_TOKEN\")))\n", + "# plt.show()" ] }, { @@ -166,45 +167,63 @@ "metadata": {}, "outputs": [], "source": [ - "# Set up model\n", - "kcm_base_model_list, kcm_nn_model = model_utils.make_one_model(model_type, hyperparameter_dict=hyperparameter_dict, embed_dict=embed_dict, config=kcm_config, skip_gtfs=skip_gtfs, load_weights=True, weight_folder=f\"{run_folder}{kcm_network_folder}models/{model_type}/logs/{model_type}/version_{fold_num}/checkpoints/\", fold_num=4)\n", - "print(f\"Evaluating: {kcm_nn_model.model_name}\")\n", - "# Set up dataset\n", - "kcm_dataset = data_loader.LoadSliceDataset(f\"{run_folder}{kcm_network_folder}deeptte_formatted/test\", kcm_config, skip_gtfs=skip_gtfs)\n", - "kcm_ngrid = grids.NGridBetter(kcm_config['grid_bounds'][0], grid_s_size)\n", - "kcm_ngrid.add_grid_content(atb_dataset.get_all_samples(keep_cols=['shingle_id','locationtime','x','y','speed_m_s','bearing']), trace_format=True)\n", - "kcm_ngrid.build_cell_lookup()\n", - "kcm_dataset.grid = kcm_ngrid\n", - "kcm_dataset.add_grid_features = kcm_nn_model.requires_grid\n", - "loader = DataLoader(kcm_dataset, sampler=SequentialSampler(kcm_dataset), collate_fn=kcm_nn_model.collate_fn, batch_size=kcm_nn_model.batch_size, pin_memory=pin_memory, num_workers=num_workers, drop_last=False)\n", + "# Make predictions for the holdout route shingles\n", + "loader = DataLoader(\n", + " holdout_dataset,\n", + " collate_fn=model.collate_fn,\n", + " batch_size=model.batch_size,\n", + " shuffle=False,\n", + " drop_last=False,\n", + " num_workers=num_workers,\n", + " pin_memory=pin_memory\n", + ")\n", "trainer = pl.Trainer(\n", - " accelerator=\"cpu\",\n", - " logger=False\n", + " accelerator=accelerator,\n", + " logger=False,\n", + " inference_mode=True,\n", + " enable_progress_bar=False,\n", + " enable_model_summary=False,\n", ")\n", - "preds_and_labels = trainer.predict(model=kcm_nn_model, dataloaders=loader)\n", + "preds_and_labels = trainer.predict(model=model, dataloaders=loader)\n", + "\n", + "# Extract predictions for full shingles, and for individual points\n", + "all_preds = np.concatenate([x['preds'] for x in preds_and_labels])\n", + "all_labels = np.concatenate([x['labels'] for x in preds_and_labels])\n", + "\n", + "all_preds_raw = np.concatenate([x['preds_raw'][x['mask']] for x in preds_and_labels])\n", + "all_labels_raw = np.concatenate([x['labels_raw'][x['mask']] for x in preds_and_labels])\n", + "\n", + "mape = np.mean(np.abs(all_preds - all_labels) / all_labels)\n", + "mape_raw = np.mean(np.abs(all_preds_raw - all_labels_raw) / all_labels_raw)\n", + "mae = np.mean(np.abs(all_preds - all_labels))\n", + "mae_raw = np.mean(np.abs(all_preds_raw - all_labels_raw))\n", + "rmse = np.sqrt(np.mean(np.square(all_preds - all_labels)))\n", + "rmse_raw = np.sqrt(np.mean(np.square(all_preds_raw - all_labels_raw)))\n", + "\n", + "print(f\"MAPE: {mape:.2f}, MAPE PT {mape_raw:.2f}\")\n", + "print(f\"MAE: {mae:.2f}, MAE PT {mae_raw:.2f}\")\n", + "print(f\"RMSE: {rmse:.2f}, RMSE PT {rmse_raw:.2f}\")\n", + "\n", + "# Extract the input features from the dataset at the point-level\n", + "shingles = []\n", + "for i in range(len(holdout_data)):\n", + " # Skip the first element of each shingle, which is the start point\n", + " shingle = holdout_data[i]['feats_n'][1:]\n", + " shingles.append(shingle)\n", + "shingles = np.concatenate(shingles)\n", "\n", - "# Extract predictions\n", - "kcm_preds = np.concatenate([x['out'][x['mask']] for x in preds_and_labels])\n", - "kcm_labels = np.concatenate([x['y'][x['mask']] for x in preds_and_labels])\n", - "# Extract data points and connect to predictions\n", - "kcm_data = [x for x in kcm_dataset]\n", - "kcm_feats = np.concatenate([x['samp'] for x in kcm_data])\n", - "kcm_ys = np.array([x['norm_label'] for x in kcm_data])\n", - "kcm_res = pd.DataFrame(kcm_feats, columns=kcm_dataset.col_names)\n", - "kcm_res['preds'] = kcm_preds\n", - "kcm_res['labels'] = kcm_labels\n", - "# Get geometries and other features for every prediciton point\n", - "kcm_res['pred_speeds'] = kcm_res['dist_calc_km']*1000 / kcm_res['preds']\n", - "kcm_res['label_speeds'] = kcm_res['dist_calc_km']*1000 / kcm_res['labels']\n", - "kcm_res['absolute_error'] = abs(kcm_res['preds'] - kcm_res['labels'])\n", - "kcm_res['hour'] = kcm_res['timeID']//60\n", - "points = gpd.points_from_xy(kcm_res['lon'], kcm_res['lat'], crs=\"EPSG:4326\")\n", - "kcm_res = gpd.GeoDataFrame(kcm_res, geometry=points)\n", - "kcm_res = kcm_res.sample(100000)\n", - "# Overall accuracy on data points\n", - "print(f\"MAE: {metrics.mean_absolute_error(kcm_labels,kcm_preds)}\")\n", - "print(f\"MAPE: {metrics.mean_absolute_percentage_error(kcm_labels,kcm_preds)}\")\n", - "kcm_res" + "# Combine the input features, labels and predictions into analysis dataframe\n", + "df = pd.DataFrame(shingles, columns=data_loader.NUM_FEAT_COLS)\n", + "df['preds'] = all_preds_raw\n", + "df['labels'] = all_labels_raw\n", + "df['residuals'] = np.abs(df['preds'] - df['labels'])\n", + "df['residuals_sq'] = (df['preds'] - df['labels']) ** 2\n", + "df = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.x, df.y), crs=\"EPSG:32148\")\n", + "\n", + "# Reproject to WGS84\n", + "df = df.to_crs(\"EPSG:4326\")\n", + "df['lon'] = df.geometry.x\n", + "df['lat'] = df.geometry.y" ] }, { @@ -213,15 +232,7 @@ "metadata": {}, "outputs": [], "source": [ - "# axes = geoplot.pointplot(atb_res, projection=geoplot.crs.AlbersEqualArea(), s=0.1)\n", - "# geoplot.kdeplot(atb_res, fill=True, cmap='coolwarm', alpha=0.5, bw_adjust=0.5, ax=axes)\n", - "# axes.set_title(f\"Point Heatmap(atb)\")\n", - "# plt.savefig(\"../plots/model_spatial_performance_atb.png\", dpi=600, bbox_inches='tight')\n", - "\n", - "# axes = geoplot.pointplot(kcm_res, projection=geoplot.crs.AlbersEqualArea(), s=0.1)\n", - "# geoplot.kdeplot(kcm_res, fill=True, cmap='coolwarm', alpha=0.5, bw_adjust=0.5, ax=axes)\n", - "# axes.set_title(f\"Point Heatmap (KCM)\")\n", - "# plt.savefig(\"../plots/model_spatial_performance_kcm.png\", dpi=600, bbox_inches='tight')" + "sns.histplot(df['residuals'])" ] }, { @@ -231,36 +242,43 @@ "outputs": [], "source": [ "fig = ff.create_hexbin_mapbox(\n", - " data_frame=kcm_res,\n", + " data_frame=df,\n", " lat=\"lat\",\n", " lon=\"lon\",\n", + " width=500,\n", + " height=700,\n", " nx_hexagon=30,\n", - " opacity=0.7,\n", - " labels={\"color\": \"SD of Speed Predictions (m/s)\"},\n", - " color=\"pred_speeds\",\n", - " agg_func=np.std,\n", + " labels={\"color\": \"MAE\"},\n", + " color=\"residuals\",\n", + " agg_func=np.mean,\n", " color_continuous_scale=\"Icefire\",\n", - " range_color=[0,10]\n", + " range_color=[2,3],\n", + " mapbox_style='open-street-map',\n", ")\n", - "fig.show()\n", - "fig.write_image(f\"../plots/within_sd_preds_hexbin_kcm.eps\")\n", - "fig.write_image(f\"../plots/within_sd_preds_hexbin_kcm.png\")\n", - "\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ "fig = ff.create_hexbin_mapbox(\n", - " data_frame=kcm_res,\n", + " data_frame=df,\n", " lat=\"lat\",\n", " lon=\"lon\",\n", - " nx_hexagon=30,\n", - " opacity=0.7,\n", - " labels={\"color\": \"SD of Speed Labels (m/s)\"},\n", - " color=\"label_speeds\",\n", + " width=500,\n", + " height=700,\n", + " nx_hexagon=20,\n", + " opacity=0.9,\n", + " labels={\"color\": \"Std. Residuals\"},\n", + " color=\"residuals\",\n", " agg_func=np.std,\n", " color_continuous_scale=\"Icefire\",\n", - " range_color=[0,10]\n", + " range_color=[2,5]\n", ")\n", - "fig.show()\n", - "fig.write_image(f\"../plots/within_sd_labels_hexbin_kcm.eps\")\n", - "fig.write_image(f\"../plots/within_sd_labels_hexbin_kcm.png\")" + "fig.show()" ] }, { @@ -270,36 +288,43 @@ "outputs": [], "source": [ "fig = ff.create_hexbin_mapbox(\n", - " data_frame=atb_res,\n", + " data_frame=df,\n", " lat=\"lat\",\n", " lon=\"lon\",\n", - " nx_hexagon=30,\n", - " opacity=0.7,\n", - " labels={\"color\": \"SD of Speed Predictions (m/s)\"},\n", - " color=\"pred_speeds\",\n", - " agg_func=np.std,\n", + " width=500,\n", + " height=700,\n", + " nx_hexagon=20,\n", + " opacity=0.9,\n", + " labels={\"color\": \"Avg. Speed (m/s)\"},\n", + " color=\"calc_speed_m_s\",\n", + " agg_func=np.mean,\n", " color_continuous_scale=\"Icefire\",\n", - " range_color=[0,10]\n", + " range_color=[0,20]\n", ")\n", - "fig.show()\n", - "fig.write_image(f\"../plots/within_sd_preds_hexbin_atb.eps\")\n", - "fig.write_image(f\"../plots/within_sd_preds_hexbin_atb.png\")\n", - "\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ "fig = ff.create_hexbin_mapbox(\n", - " data_frame=atb_res,\n", + " data_frame=df,\n", " lat=\"lat\",\n", " lon=\"lon\",\n", - " nx_hexagon=30,\n", - " opacity=0.7,\n", - " labels={\"color\": \"SD of Speed Labels (m/s)\"},\n", - " color=\"label_speeds\",\n", + " width=500,\n", + " height=700,\n", + " nx_hexagon=20,\n", + " opacity=0.9,\n", + " labels={\"color\": \"Std. Speed (m/s)\"},\n", + " color=\"calc_speed_m_s\",\n", " agg_func=np.std,\n", " color_continuous_scale=\"Icefire\",\n", - " range_color=[0,10]\n", + " range_color=[0,8]\n", ")\n", - "fig.show()\n", - "fig.write_image(f\"../plots/within_sd_labels_hexbin_atb.eps\")\n", - "fig.write_image(f\"../plots/within_sd_labels_hexbin_atb.png\")" + "fig.show()" ] }, { @@ -308,31 +333,19 @@ "metadata": {}, "outputs": [], "source": [ - "fig, axes = plt.subplots(3,2,figsize=(6,8))\n", - "fig.tight_layout()\n", - "axes = axes.flatten()\n", - "sns.histplot(kcm_res, x=\"pred_speeds\", ax=axes[0])\n", - "axes[0].set_title(f\"KCM Prediction Metrics\")\n", - "axes[0].set_xlabel(\"Predicted Speed (m/s)\")\n", - "axes[0].set_xlim(0,30)\n", - "sns.histplot(atb_res, x=\"pred_speeds\", ax=axes[1])\n", - "axes[1].set_title(f\"AtB Prediction Metrics\")\n", - "axes[1].set_xlabel(\"Predicted Speed (m/s)\")\n", - "axes[1].set_xlim(0,30)\n", - "sns.histplot(kcm_res, x=\"label_speeds\", ax=axes[2])\n", - "axes[2].set_xlabel(\"Label Speed (m/s)\")\n", - "axes[2].set_xlim(0,30)\n", - "sns.histplot(atb_res, x=\"label_speeds\", ax=axes[3])\n", - "axes[3].set_xlabel(\"Label Speed (m/s)\")\n", - "axes[3].set_xlim(0,30)\n", - "sns.histplot(kcm_res, x=\"absolute_error\", ax=axes[4])\n", - "axes[4].set_xlabel(\"Absolute Error (s)\")\n", - "axes[4].set_xlim(0,30)\n", - "sns.histplot(atb_res, x=\"absolute_error\", ax=axes[5])\n", - "axes[5].set_xlabel(\"Absolute Error (s)\")\n", - "axes[5].set_xlim(0,30)\n", - "plt.savefig(\"../plots/model_prediction_distribution_comparison.eps\", format='eps', dpi=600, bbox_inches='tight')\n", - "plt.savefig(\"../plots/model_prediction_distribution_comparison.png\", format='png', dpi=600, bbox_inches='tight')" + "# plot_df = df[df['MAPE']<1].copy()\n", + "# fig, axes = plt.subplots(1, 1, figsize=(10,10))\n", + "# plt.hexbin(plot_df.x, plot_df.y, plot_df.MAPE, cmap='plasma', gridsize=15)\n", + "# sns.kdeplot(ax=axes, x=plot_df.x, y=plot_df.y, weights=plot_df.MAPE, cmap=\"plasma\", fill=True, bw_adjust=.2)\n", + "# cx.add_basemap(ax=axes, crs=plot_df.crs.to_string(), alpha=0.3, source=cx.providers.MapBox(accessToken=os.getenv(key=\"MAPBOX_TOKEN\")))\n", + "# plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Inference on All Routes" ] }, { @@ -341,38 +354,87 @@ "metadata": {}, "outputs": [], "source": [ - "fig, axes = plt.subplots(3,2,figsize=(6,8))\n", - "fig.tight_layout()\n", - "axes = axes.flatten()\n", - "sns.lineplot(kcm_res, x=\"hour\", y=\"pred_speeds\", ax=axes[0])\n", - "axes[0].set_title(f\"KCM Prediction Metrics\")\n", - "axes[0].set_xlabel(\"Hour of Day\")\n", - "axes[0].set_xlim(0,24)\n", - "sns.lineplot(atb_res, x=\"hour\", y=\"pred_speeds\", ax=axes[1])\n", - "axes[1].set_title(f\"AtB Prediction Metrics\")\n", - "axes[1].set_xlabel(\"Hour of Day\")\n", - "axes[1].set_xlim(0,24)\n", - "sns.lineplot(kcm_res, x=\"hour\", y=\"label_speeds\", ax=axes[2])\n", - "axes[2].set_xlabel(\"Hour of Day\")\n", - "axes[2].set_xlim(0,24)\n", - "sns.lineplot(atb_res, x=\"hour\", y=\"label_speeds\", ax=axes[3])\n", - "axes[3].set_xlabel(\"Hour of Day\")\n", - "axes[3].set_xlim(0,24)\n", - "sns.lineplot(kcm_res, x=\"hour\", y=\"absolute_error\", ax=axes[4])\n", - "axes[4].set_xlabel(\"Hour of Day\")\n", - "axes[4].set_xlim(0,24)\n", - "sns.lineplot(atb_res, x=\"hour\", y=\"absolute_error\", ax=axes[5])\n", - "axes[5].set_xlabel(\"Hour of Day\")\n", - "axes[5].set_xlim(0,24)\n", - "plt.savefig(\"../plots/model_hourly_comparison.eps\", format='eps', dpi=600, bbox_inches='tight')\n", - "plt.savefig(\"../plots/model_hourly_comparison.png\", format='png', dpi=600, bbox_inches='tight')" + "# Make predictions for all route shingles\n", + "loader = DataLoader(\n", + " test_dataset,\n", + " collate_fn=model.collate_fn,\n", + " batch_size=model.batch_size,\n", + " shuffle=False,\n", + " drop_last=False,\n", + " num_workers=num_workers,\n", + " pin_memory=pin_memory\n", + ")\n", + "trainer = pl.Trainer(\n", + " accelerator=accelerator,\n", + " logger=False,\n", + " inference_mode=True,\n", + " enable_progress_bar=False,\n", + " enable_model_summary=False,\n", + ")\n", + "preds_and_labels = trainer.predict(model=model, dataloaders=loader)\n", + "\n", + "# Extract predictions for full shingles, and for individual points\n", + "all_preds = np.concatenate([x['preds'] for x in preds_and_labels])\n", + "all_labels = np.concatenate([x['labels'] for x in preds_and_labels])\n", + "\n", + "all_preds_raw = np.concatenate([x['preds_raw'][x['mask']] for x in preds_and_labels])\n", + "all_labels_raw = np.concatenate([x['labels_raw'][x['mask']] for x in preds_and_labels])\n", + "\n", + "mape = np.mean(np.abs(all_preds - all_labels) / all_labels)\n", + "mape_raw = np.mean(np.abs(all_preds_raw - all_labels_raw) / all_labels_raw)\n", + "mae = np.mean(np.abs(all_preds - all_labels))\n", + "mae_raw = np.mean(np.abs(all_preds_raw - all_labels_raw))\n", + "rmse = np.sqrt(np.mean(np.square(all_preds - all_labels)))\n", + "rmse_raw = np.sqrt(np.mean(np.square(all_preds_raw - all_labels_raw)))\n", + "\n", + "print(f\"MAPE: {mape:.2f}, MAPE PT {mape_raw:.2f}\")\n", + "print(f\"MAE: {mae:.2f}, MAE PT {mae_raw:.2f}\")\n", + "print(f\"RMSE: {rmse:.2f}, RMSE PT {rmse_raw:.2f}\")\n", + "\n", + "# Extract the input features from the dataset at the point-level\n", + "shingles = []\n", + "for i in range(len(test_data)):\n", + " # Skip the first element of each shingle, which is the start point\n", + " shingle = test_data[i]['feats_n'][1:]\n", + " shingles.append(shingle)\n", + "shingles = np.concatenate(shingles)\n", + "\n", + "# Combine the input features, labels and predictions into analysis dataframe\n", + "df = pd.DataFrame(shingles, columns=data_loader.NUM_FEAT_COLS)\n", + "df['preds'] = all_preds_raw\n", + "df['labels'] = all_labels_raw\n", + "df['residuals'] = np.abs(df['preds'] - df['labels'])\n", + "df['residuals_sq'] = (df['preds'] - df['labels']) ** 2\n", + "df['pred_speed_m_s'] = df['calc_dist_m'] / df['preds']\n", + "df = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.x, df.y), crs=\"EPSG:32148\")\n", + "\n", + "# Reproject to WGS84\n", + "df = df.to_crs(\"EPSG:4326\")\n", + "df['lon'] = df.geometry.x\n", + "df['lat'] = df.geometry.y" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "### Inference on Entire Network" + "# fig = ff.create_hexbin_mapbox(\n", + "# data_frame=df,\n", + "# lat=\"lat\",\n", + "# lon=\"lon\",\n", + "# width=1000,\n", + "# height=1200,\n", + "# nx_hexagon=200,\n", + "# labels={\"color\": \"MAE\"},\n", + "# color=\"residuals\",\n", + "# agg_func=np.mean,\n", + "# color_continuous_scale=\"Icefire\",\n", + "# range_color=[1.5,3.5],\n", + "# mapbox_style='open-street-map',\n", + "# )\n", + "# fig.show()" ] }, { @@ -381,43 +443,21 @@ "metadata": {}, "outputs": [], "source": [ - "# Create grid of regularly spaced fake shingles to feed model\n", - "inference_shingles = data_utils.create_grid_of_shingles(100, atb_config['grid_bounds'][0], atb_config['coord_ref_center'][0])\n", - "\n", - "# Make predictions for the fake shingles\n", - "atb_dataset = data_loader.ContentDataset(inference_shingles, atb_config, skip_gtfs=True)\n", - "loader = DataLoader(atb_dataset, sampler=SequentialSampler(atb_dataset), collate_fn=atb_nn_model.collate_fn, batch_size=atb_nn_model.batch_size, pin_memory=pin_memory, num_workers=num_workers, drop_last=False)\n", - "trainer = pl.Trainer(\n", - " accelerator=\"cpu\",\n", - " logger=False\n", - ")\n", - "preds_and_labels = trainer.predict(model=atb_nn_model, dataloaders=loader)\n", - "atb_preds = np.concatenate([x['out'][x['mask']] for x in preds_and_labels])\n", - "atb_labels = np.concatenate([x['y'][x['mask']] for x in preds_and_labels])\n", - "\n", - "# Extract predictions\n", - "atb_preds = np.concatenate([x['out'][x['mask']] for x in preds_and_labels])\n", - "atb_labels = np.concatenate([x['y'][x['mask']] for x in preds_and_labels])\n", - "# Extract data points and connect to predictions\n", - "atb_data = [x for x in atb_dataset]\n", - "atb_feats = np.concatenate([x['samp'] for x in atb_data])\n", - "atb_ys = np.array([x['norm_label'] for x in atb_data])\n", - "atb_res = pd.DataFrame(atb_feats, columns=atb_dataset.col_names)\n", - "atb_res['preds'] = atb_preds\n", - "atb_res['labels'] = atb_labels\n", - "# Get geometries and other features for every prediciton point\n", - "atb_res['pred_speeds'] = atb_res['dist_calc_km']*1000 / atb_res['preds']\n", - "atb_res['label_speeds'] = atb_res['dist_calc_km']*1000 / atb_res['labels']\n", - "atb_res['absolute_error'] = abs(atb_res['preds'] - atb_res['labels'])\n", - "atb_res['hour'] = atb_res['timeID']//60\n", - "points = gpd.points_from_xy(atb_res['lon'], atb_res['lat'], crs=\"EPSG:4326\")\n", - "atb_res = gpd.GeoDataFrame(atb_res, geometry=points)\n", - "# Transform x and y to replace dummy lat and lon for mapping\n", - "transformer = pyproj.Transformer.from_crs(atb_crs, default_crs)\n", - "atb_res['lat'], atb_res['lon'] = transformer.transform(atb_res['x'], atb_res['y'])\n", - "points = gpd.points_from_xy(atb_res['lon'], atb_res['lat'], crs=default_crs)\n", - "atb_res = gpd.GeoDataFrame(atb_res, geometry=points)\n", - "atb_res" + "# fig = ff.create_hexbin_mapbox(\n", + "# data_frame=df,\n", + "# lat=\"lat\",\n", + "# lon=\"lon\",\n", + "# width=1000,\n", + "# height=1200,\n", + "# nx_hexagon=100,\n", + "# labels={\"color\": \"MAE\"},\n", + "# color=\"calc_speed_m_s\",\n", + "# agg_func=np.mean,\n", + "# color_continuous_scale=\"Icefire\",\n", + "# range_color=[0,30],\n", + "# mapbox_style='open-street-map',\n", + "# )\n", + "# fig.show()" ] }, { @@ -426,43 +466,21 @@ "metadata": {}, "outputs": [], "source": [ - "# Create grid of regularly spaced fake shingles to feed model\n", - "inference_shingles = data_utils.create_grid_of_shingles(100, kcm_config['grid_bounds'][0], kcm_config['coord_ref_center'][0])\n", - "\n", - "# Make predictions for the fake shingles\n", - "kcm_dataset = data_loader.ContentDataset(inference_shingles, kcm_config, skip_gtfs=True)\n", - "loader = DataLoader(kcm_dataset, sampler=SequentialSampler(kcm_dataset), collate_fn=kcm_nn_model.collate_fn, batch_size=kcm_nn_model.batch_size, pin_memory=pin_memory, num_workers=num_workers, drop_last=False)\n", - "trainer = pl.Trainer(\n", - " accelerator=\"cpu\",\n", - " logger=False\n", - ")\n", - "preds_and_labels = trainer.predict(model=kcm_nn_model, dataloaders=loader)\n", - "kcm_preds = np.concatenate([x['out'][x['mask']] for x in preds_and_labels])\n", - "kcm_labels = np.concatenate([x['y'][x['mask']] for x in preds_and_labels])\n", - "\n", - "# Extract predictions\n", - "kcm_preds = np.concatenate([x['out'][x['mask']] for x in preds_and_labels])\n", - "kcm_labels = np.concatenate([x['y'][x['mask']] for x in preds_and_labels])\n", - "# Extract data points and connect to predictions\n", - "kcm_data = [x for x in kcm_dataset]\n", - "kcm_feats = np.concatenate([x['samp'] for x in kcm_data])\n", - "kcm_ys = np.array([x['norm_label'] for x in kcm_data])\n", - "kcm_res = pd.DataFrame(kcm_feats, columns=kcm_dataset.col_names)\n", - "kcm_res['preds'] = kcm_preds\n", - "kcm_res['labels'] = kcm_labels\n", - "# Get geometries and other features for every prediciton point\n", - "kcm_res['pred_speeds'] = kcm_res['dist_calc_km']*1000 / kcm_res['preds']\n", - "kcm_res['label_speeds'] = kcm_res['dist_calc_km']*1000 / kcm_res['labels']\n", - "kcm_res['absolute_error'] = abs(kcm_res['preds'] - kcm_res['labels'])\n", - "kcm_res['hour'] = kcm_res['timeID']//60\n", - "points = gpd.points_from_xy(kcm_res['lon'], kcm_res['lat'], crs=\"EPSG:4326\")\n", - "kcm_res = gpd.GeoDataFrame(kcm_res, geometry=points)\n", - "# Transform x and y to replace dummy lat and lon for mapping\n", - "transformer = pyproj.Transformer.from_crs(kcm_crs, default_crs)\n", - "kcm_res['lat'], kcm_res['lon'] = transformer.transform(kcm_res['x'], kcm_res['y'])\n", - "points = gpd.points_from_xy(kcm_res['lon'], kcm_res['lat'], crs=default_crs)\n", - "kcm_res = gpd.GeoDataFrame(kcm_res, geometry=points)\n", - "kcm_res" + "# fig = ff.create_hexbin_mapbox(\n", + "# data_frame=df,\n", + "# lat=\"lat\",\n", + "# lon=\"lon\",\n", + "# width=1000,\n", + "# height=1200,\n", + "# nx_hexagon=100,\n", + "# labels={\"color\": \"MAE\"},\n", + "# color=\"calc_speed_m_s\",\n", + "# agg_func=np.std,\n", + "# color_continuous_scale=\"Icefire\",\n", + "# range_color=[0,10],\n", + "# mapbox_style='open-street-map',\n", + "# )\n", + "# fig.show()" ] }, { @@ -471,36 +489,52 @@ "metadata": {}, "outputs": [], "source": [ - "fig = ff.create_hexbin_mapbox(\n", - " data_frame=atb_res,\n", - " lat=\"lat\",\n", - " lon=\"lon\",\n", - " nx_hexagon=50,\n", - " opacity=0.7,\n", - " labels={\"color\": \"Mean Speed (m/s)\"},\n", - " color=\"pred_speeds\",\n", - " agg_func=np.mean,\n", - " color_continuous_scale=\"Icefire_r\",\n", - ")\n", - "fig.show()\n", - "fig.write_image(f\"../plots/mesh_mean_pred_speeds_hexbin_atb.eps\")\n", - "fig.write_image(f\"../plots/mesh_mean_pred_speeds_hexbin_atb.png\")\n", - "\n", - "fig = ff.create_hexbin_mapbox(\n", - " data_frame=kcm_res,\n", - " lat=\"lat\",\n", - " lon=\"lon\",\n", - " nx_hexagon=50,\n", - " opacity=0.7,\n", - " labels={\"color\": \"Mean Speed (m/s)\"},\n", - " color=\"pred_speeds\",\n", - " agg_func=np.mean,\n", - " color_continuous_scale=\"Icefire_r\",\n", - ")\n", - "fig.show()\n", - "fig.write_image(f\"../plots/mesh_mean_pred_speeds_hexbin_kcm.eps\")\n", - "fig.write_image(f\"../plots/mesh_mean_pred_speeds_hexbin_kcm.png\")" + "# fig = ff.create_hexbin_mapbox(\n", + "# data_frame=df,\n", + "# lat=\"lat\",\n", + "# lon=\"lon\",\n", + "# width=1000,\n", + "# height=1200,\n", + "# nx_hexagon=100,\n", + "# labels={\"color\": \"MAE\"},\n", + "# color=\"pred_speed_m_s\",\n", + "# agg_func=np.mean,\n", + "# color_continuous_scale=\"Icefire\",\n", + "# range_color=[0,30],\n", + "# mapbox_style='open-street-map',\n", + "# )\n", + "# fig.show()" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# fig = ff.create_hexbin_mapbox(\n", + "# data_frame=df,\n", + "# lat=\"lat\",\n", + "# lon=\"lon\",\n", + "# width=1000,\n", + "# height=1200,\n", + "# nx_hexagon=100,\n", + "# labels={\"color\": \"MAE\"},\n", + "# color=\"pred_speed_m_s\",\n", + "# agg_func=np.std,\n", + "# color_continuous_scale=\"Icefire\",\n", + "# range_color=[0,100],\n", + "# mapbox_style='open-street-map',\n", + "# )\n", + "# fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { @@ -519,7 +553,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.6" + "version": "3.9.18" }, "orig_nbformat": 4, "vscode": { diff --git a/notebooks/explore_other_cities.ipynb b/notebooks/explore_other_cities.ipynb index ef1367b..80a2180 100644 --- a/notebooks/explore_other_cities.ipynb +++ b/notebooks/explore_other_cities.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 2, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -27,18 +27,231 @@ "import seaborn as sns\n", "import shapely\n", "import statsmodels.api as sm\n", + "import torch\n", "from torch.utils.data import DataLoader\n", "\n", "from openbustools import plotting, spatial, standardfeeds\n", "from openbustools.traveltime import data_loader, model_utils\n", "from openbustools.drivecycle import trajectory\n", - "from openbustools.drivecycle.physics import conditions, energy, vehicle" + "from openbustools.drivecycle.physics import conditions, energy, vehicle\n", + "\n", + "from srai.embedders import GTFS2VecEmbedder\n", + "from srai.joiners import IntersectionJoiner\n", + "from srai.loaders import GTFSLoader\n", + "from srai.loaders.osm_loaders.filters import HEX2VEC_FILTER\n", + "from srai.neighbourhoods.h3_neighbourhood import H3Neighbourhood\n", + "from srai.regionalizers import H3Regionalizer, geocode_to_region_gdf\n", + "from srai.plotting import plot_regions, plot_numeric_data" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, + "outputs": [], + "source": [ + "if torch.cuda.is_available():\n", + " num_workers=4\n", + " pin_memory=True\n", + " accelerator=\"cuda\"\n", + "else:\n", + " num_workers=0\n", + " pin_memory=False\n", + " accelerator=\"cpu\"\n", + "\n", + "logging.getLogger(\"lightning\").setLevel(logging.ERROR)\n", + "logging.getLogger(\"pytorch_lightning\").setLevel(logging.ERROR)\n", + "logging.getLogger(\"lightning.pytorch.accelerators.cuda\").setLevel(logging.ERROR)\n", + "\n", + "cleaned_sources = pd.read_csv(Path('..', 'data', 'cleaned_sources.csv'))" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "London Transit Commission\n", + "Barrie Transit\n", + "Big Blue Bus\n", + "Mountain View Transportation Management Association (MVgo)\n", + "Capital Metro\n", + "Regional Transportation District (RTD)\n", + "Metro St. Louis\n", + "Intercity Transit\n", + "Duluth Transit\n", + "Central Florida Regional Transit Authority (LYNX)\n", + "Votran\n", + "Nashville Metropolitan Transit Authority (Nashville MTA)\n", + "Transit Authority of River City (TARC)\n", + "Metropolitan Atlanta Rapid Transit Authority (MARTA)\n", + "Springfield Mass Transit District (SMTD)\n", + "Metro Transit\n", + "Port Authority of Allegheny County\n", + "Massachusetts Bay Transportation Authority (MBTA)\n", + "Arlington Transit\n", + "Rochester-Genesee Regional Transportation Authority (RGRTA)\n", + "BC Transit (Victoria Regional Transit System)\n", + "Edmonton Transit System\n", + "Saskatoon Transit\n", + "Hamilton Street Railway\n", + "MiWay\n", + "Kingston Transit\n", + "Halifax Transit\n", + "Thunder Bay Transit\n", + "Port Alberni\n", + "Commerce Municipal Bus Lines\n", + "King County Metro\n", + "Stanislaus Regional Transit Authority (StanRTA)\n" + ] + } + ], + "source": [ + "# Test model on each city\n", + "res = {}\n", + "test_sources = cleaned_sources.iloc[:]\n", + "for i, row in test_sources.iterrows():\n", + " model = model_utils.load_model(\"../logs/saved_models/\", \"kcm_1month_resampled\", \"GRU\", 0)\n", + " model.eval()\n", + " test_dates = ['2024_01_03.pkl']\n", + " train_data_folders = [Path('..', 'data', 'other_feeds', f\"{row['uuid']}_realtime\", 'processed')]\n", + " available_files = [x.name for x in Path('..', 'data', 'other_feeds', f\"{row['uuid']}_realtime\", 'processed').glob('*.pkl')]\n", + " if test_dates[0] not in available_files:\n", + " continue\n", + " print(row['provider'])\n", + " test_data, holdout_routes, test_config = data_loader.load_h5(train_data_folders, test_dates, holdout_routes=model.holdout_routes, config=model.config)\n", + " test_dataset = data_loader.H5Dataset(test_data)\n", + " test_dataset.include_grid = model.include_grid\n", + " test_loader = DataLoader(\n", + " test_dataset,\n", + " collate_fn=model.collate_fn,\n", + " batch_size=model.batch_size,\n", + " shuffle=False,\n", + " drop_last=False,\n", + " num_workers=num_workers,\n", + " pin_memory=pin_memory\n", + " )\n", + " trainer = pl.Trainer(\n", + " accelerator=accelerator,\n", + " logger=False,\n", + " inference_mode=True,\n", + " enable_progress_bar=False\n", + " )\n", + " preds_and_labels = trainer.predict(model=model, dataloaders=test_loader)\n", + " preds = np.concatenate([x['preds'] for x in preds_and_labels])\n", + " labels = np.concatenate([x['labels'] for x in preds_and_labels])\n", + " res[row['uuid']] = {'preds':preds, 'labels':labels}" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "London Transit Commission\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/zack/Desktop/open_bus_tools/.venv/lib/python3.9/site-packages/lightning/pytorch/trainer/configuration_validator.py:74: You defined a `validation_step` but have no `val_dataloader`. Skipping val loop.\n", + "/home/zack/Desktop/open_bus_tools/.venv/lib/python3.9/site-packages/lightning/pytorch/callbacks/model_checkpoint.py:639: Checkpoint directory /home/zack/Desktop/open_bus_tools/notebooks/checkpoints exists and is not empty.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Barrie Transit\n", + "Big Blue Bus\n", + "Mountain View Transportation Management Association (MVgo)\n", + "Capital Metro\n", + "Regional Transportation District (RTD)\n", + "Metro St. Louis\n", + "Intercity Transit\n", + "Duluth Transit\n", + "Central Florida Regional Transit Authority (LYNX)\n", + "Votran\n", + "Nashville Metropolitan Transit Authority (Nashville MTA)\n", + "Transit Authority of River City (TARC)\n", + "Metropolitan Atlanta Rapid Transit Authority (MARTA)\n", + "Springfield Mass Transit District (SMTD)\n", + "Metro Transit\n", + "Port Authority of Allegheny County\n", + "Massachusetts Bay Transportation Authority (MBTA)\n", + "Arlington Transit\n", + "Rochester-Genesee Regional Transportation Authority (RGRTA)\n", + "BC Transit (Victoria Regional Transit System)\n", + "Edmonton Transit System\n", + "Saskatoon Transit\n", + "Hamilton Street Railway\n", + "MiWay\n", + "Kingston Transit\n", + "Halifax Transit\n", + "Thunder Bay Transit\n", + "Port Alberni\n", + "Commerce Municipal Bus Lines\n", + "King County Metro\n", + "Stanislaus Regional Transit Authority (StanRTA)\n" + ] + } + ], + "source": [ + "# Tune, then re-test\n", + "res_tuned = {}\n", + "test_sources = cleaned_sources.iloc[:]\n", + "for i, row in test_sources.iterrows():\n", + " test_dates = ['2024_01_03.pkl']\n", + " train_data_folders = [Path('..', 'data', 'other_feeds', f\"{row['uuid']}_realtime\", 'processed')]\n", + " available_files = [x.name for x in Path('..', 'data', 'other_feeds', f\"{row['uuid']}_realtime\", 'processed').glob('*.pkl')]\n", + " if test_dates[0] not in available_files:\n", + " continue\n", + " print(row['provider'])\n", + " res_tuned[row['uuid']] = {}\n", + " for n_batches in [1, 2, 4]:\n", + " model = model_utils.load_model(\"../logs/saved_models/\", \"kcm_1month_resampled\", \"GRU\", 0)\n", + " model.eval()\n", + " test_data, holdout_routes, test_config = data_loader.load_h5(train_data_folders, test_dates, holdout_routes=model.holdout_routes, config=model.config)\n", + " test_dataset = data_loader.H5Dataset(test_data)\n", + " test_dataset.include_grid = model.include_grid\n", + " test_loader = DataLoader(\n", + " test_dataset,\n", + " collate_fn=model.collate_fn,\n", + " batch_size=model.batch_size,\n", + " shuffle=False,\n", + " drop_last=False,\n", + " num_workers=num_workers,\n", + " pin_memory=pin_memory\n", + " )\n", + " trainer = pl.Trainer(\n", + " accelerator=accelerator,\n", + " logger=False,\n", + " inference_mode=True,\n", + " enable_progress_bar=False,\n", + " enable_model_summary=False,\n", + " max_epochs=10,\n", + " limit_train_batches=n_batches\n", + " )\n", + " trainer.fit(model=model, train_dataloaders=test_loader)\n", + " preds_and_labels = trainer.predict(model=model, dataloaders=test_loader)\n", + " preds = np.concatenate([x['preds'] for x in preds_and_labels])\n", + " labels = np.concatenate([x['labels'] for x in preds_and_labels])\n", + " res_tuned[row['uuid']][f\"{n_batches}_batches\"] = {'preds':preds, 'labels':labels}" + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "metadata": {}, "outputs": [ { "data": { @@ -61,163 +274,178 @@ " <thead>\n", " <tr style=\"text-align: right;\">\n", " <th></th>\n", - " <th>country_code</th>\n", - " <th>subdivision_name</th>\n", - " <th>municipality</th>\n", " <th>provider</th>\n", - " <th>static_url</th>\n", - " <th>realtime_url</th>\n", - " <th>min_lon</th>\n", - " <th>max_lon</th>\n", - " <th>min_lat</th>\n", - " <th>max_lat</th>\n", - " <th>uuid</th>\n", - " <th>tz_str</th>\n", - " <th>epsg_code</th>\n", + " <th>n_batches</th>\n", + " <th>experiment</th>\n", + " <th>mae</th>\n", + " <th>mape</th>\n", + " <th>rmse</th>\n", + " <th>ex_var</th>\n", + " <th>r_score</th>\n", " </tr>\n", " </thead>\n", " <tbody>\n", " <tr>\n", - " <th>0</th>\n", - " <td>CA</td>\n", - " <td>Ontario</td>\n", - " <td>London</td>\n", - " <td>London Transit Commission</td>\n", - " <td>http://www.londontransit.ca/gtfsfeed/google_tr...</td>\n", - " <td>http://gtfs.ltconline.ca/Vehicle/VehiclePositi...</td>\n", - " <td>-81.363110</td>\n", - " <td>-81.137591</td>\n", - " <td>42.905244</td>\n", - " <td>43.051188</td>\n", - " <td>0fd41f26-6bc5-45f0-88ff-594bde1e8c24</td>\n", - " <td>America/Toronto</td>\n", - " <td>32617</td>\n", + " <th>18</th>\n", + " <td>Arlington Transit</td>\n", + " <td>0_batches</td>\n", + " <td>not-tuned</td>\n", + " <td>69.569825</td>\n", + " <td>0.224919</td>\n", + " <td>8.340853</td>\n", + " <td>0.775192</td>\n", + " <td>0.771584</td>\n", + " </tr>\n", + " <tr>\n", + " <th>20</th>\n", + " <td>BC Transit (Victoria Regional Transit System)</td>\n", + " <td>0_batches</td>\n", + " <td>not-tuned</td>\n", + " <td>111.589098</td>\n", + " <td>0.259579</td>\n", + " <td>10.563574</td>\n", + " <td>0.777157</td>\n", + " <td>0.728569</td>\n", " </tr>\n", " <tr>\n", " <th>1</th>\n", - " <td>CA</td>\n", - " <td>Ontario</td>\n", - " <td>Barrie</td>\n", " <td>Barrie Transit</td>\n", - " <td>http://www.myridebarrie.ca/gtfs/Google_transit...</td>\n", - " <td>http://www.myridebarrie.ca/gtfs/GTFS_VehiclePo...</td>\n", - " <td>-79.740632</td>\n", - " <td>-79.610896</td>\n", - " <td>44.321804</td>\n", - " <td>44.420207</td>\n", - " <td>3196b95d-1eb2-4e05-bb36-7a2099d2348b</td>\n", - " <td>America/Toronto</td>\n", - " <td>32617</td>\n", + " <td>0_batches</td>\n", + " <td>not-tuned</td>\n", + " <td>116.481377</td>\n", + " <td>0.237706</td>\n", + " <td>10.792654</td>\n", + " <td>0.793000</td>\n", + " <td>0.792672</td>\n", " </tr>\n", " <tr>\n", " <th>2</th>\n", - " <td>US</td>\n", - " <td>California</td>\n", - " <td>Santa Monica</td>\n", " <td>Big Blue Bus</td>\n", - " <td>http://gtfs.bigbluebus.com/current.zip</td>\n", - " <td>http://gtfs.bigbluebus.com/vehiclepositions.bin</td>\n", - " <td>-118.549205</td>\n", - " <td>-118.237266</td>\n", - " <td>33.929498</td>\n", - " <td>34.075669</td>\n", - " <td>dbeb97e3-012d-4c23-93aa-a80696884340</td>\n", - " <td>America/Los_Angeles</td>\n", - " <td>32611</td>\n", - " </tr>\n", - " <tr>\n", - " <th>3</th>\n", - " <td>US</td>\n", - " <td>California</td>\n", - " <td>Mountain View</td>\n", - " <td>Mountain View Transportation Management Associ...</td>\n", - " <td>http://data.trilliumtransit.com/gtfs/mountainv...</td>\n", - " <td>https://ridemvgo.org/gtfs-rt/vehiclepositions</td>\n", - " <td>-122.111591</td>\n", - " <td>-122.047584</td>\n", - " <td>37.387656</td>\n", - " <td>37.431429</td>\n", - " <td>6885d33d-2f4e-4ac7-84a7-2e4747c723f0</td>\n", - " <td>America/Los_Angeles</td>\n", - " <td>32610</td>\n", + " <td>0_batches</td>\n", + " <td>not-tuned</td>\n", + " <td>346.089513</td>\n", + " <td>0.434805</td>\n", + " <td>18.603481</td>\n", + " <td>0.686116</td>\n", + " <td>0.263355</td>\n", " </tr>\n", " <tr>\n", " <th>4</th>\n", - " <td>US</td>\n", - " <td>California</td>\n", - " <td>Merced</td>\n", - " <td>Merced County Transit (The Bus)</td>\n", - " <td>http://data.trilliumtransit.com/gtfs/mercedthe...</td>\n", - " <td>https://thebuslive.com/gtfs-rt/vehiclepositions</td>\n", - " <td>-121.020240</td>\n", - " <td>-120.248000</td>\n", - " <td>36.964503</td>\n", - " <td>37.521786</td>\n", - " <td>ea803659-c28f-43ac-8c01-619d767e0b1d</td>\n", - " <td>America/Los_Angeles</td>\n", - " <td>32610</td>\n", + " <td>Capital Metro</td>\n", + " <td>0_batches</td>\n", + " <td>not-tuned</td>\n", + " <td>215.147161</td>\n", + " <td>0.263385</td>\n", + " <td>14.667896</td>\n", + " <td>0.743784</td>\n", + " <td>0.720664</td>\n", + " </tr>\n", + " <tr>\n", + " <th>9</th>\n", + " <td>Central Florida Regional Transit Authority (LYNX)</td>\n", + " <td>0_batches</td>\n", + " <td>not-tuned</td>\n", + " <td>161.216318</td>\n", + " <td>0.260830</td>\n", + " <td>12.697099</td>\n", + " <td>0.550255</td>\n", + " <td>0.526094</td>\n", + " </tr>\n", + " <tr>\n", + " <th>29</th>\n", + " <td>Commerce Municipal Bus Lines</td>\n", + " <td>0_batches</td>\n", + " <td>not-tuned</td>\n", + " <td>205.992501</td>\n", + " <td>0.269350</td>\n", + " <td>14.352439</td>\n", + " <td>0.683146</td>\n", + " <td>0.656156</td>\n", + " </tr>\n", + " <tr>\n", + " <th>8</th>\n", + " <td>Duluth Transit</td>\n", + " <td>0_batches</td>\n", + " <td>not-tuned</td>\n", + " <td>126.026914</td>\n", + " <td>0.285953</td>\n", + " <td>11.226171</td>\n", + " <td>0.785196</td>\n", + " <td>0.784499</td>\n", + " </tr>\n", + " <tr>\n", + " <th>21</th>\n", + " <td>Edmonton Transit System</td>\n", + " <td>0_batches</td>\n", + " <td>not-tuned</td>\n", + " <td>115.017353</td>\n", + " <td>0.258047</td>\n", + " <td>10.724614</td>\n", + " <td>0.777081</td>\n", + " <td>0.776323</td>\n", + " </tr>\n", + " <tr>\n", + " <th>26</th>\n", + " <td>Halifax Transit</td>\n", + " <td>0_batches</td>\n", + " <td>not-tuned</td>\n", + " <td>129.650576</td>\n", + " <td>0.325215</td>\n", + " <td>11.386421</td>\n", + " <td>0.680104</td>\n", + " <td>0.676025</td>\n", " </tr>\n", " </tbody>\n", "</table>\n", "</div>" ], "text/plain": [ - " country_code subdivision_name municipality \\\n", - "0 CA Ontario London \n", - "1 CA Ontario Barrie \n", - "2 US California Santa Monica \n", - "3 US California Mountain View \n", - "4 US California Merced \n", - "\n", - " provider \\\n", - "0 London Transit Commission \n", - "1 Barrie Transit \n", - "2 Big Blue Bus \n", - "3 Mountain View Transportation Management Associ... \n", - "4 Merced County Transit (The Bus) \n", - "\n", - " static_url \\\n", - "0 http://www.londontransit.ca/gtfsfeed/google_tr... \n", - "1 http://www.myridebarrie.ca/gtfs/Google_transit... \n", - "2 http://gtfs.bigbluebus.com/current.zip \n", - "3 http://data.trilliumtransit.com/gtfs/mountainv... \n", - "4 http://data.trilliumtransit.com/gtfs/mercedthe... \n", + " provider n_batches experiment \\\n", + "18 Arlington Transit 0_batches not-tuned \n", + "20 BC Transit (Victoria Regional Transit System) 0_batches not-tuned \n", + "1 Barrie Transit 0_batches not-tuned \n", + "2 Big Blue Bus 0_batches not-tuned \n", + "4 Capital Metro 0_batches not-tuned \n", + "9 Central Florida Regional Transit Authority (LYNX) 0_batches not-tuned \n", + "29 Commerce Municipal Bus Lines 0_batches not-tuned \n", + "8 Duluth Transit 0_batches not-tuned \n", + "21 Edmonton Transit System 0_batches not-tuned \n", + "26 Halifax Transit 0_batches not-tuned \n", "\n", - " realtime_url min_lon max_lon \\\n", - "0 http://gtfs.ltconline.ca/Vehicle/VehiclePositi... -81.363110 -81.137591 \n", - "1 http://www.myridebarrie.ca/gtfs/GTFS_VehiclePo... -79.740632 -79.610896 \n", - "2 http://gtfs.bigbluebus.com/vehiclepositions.bin -118.549205 -118.237266 \n", - "3 https://ridemvgo.org/gtfs-rt/vehiclepositions -122.111591 -122.047584 \n", - "4 https://thebuslive.com/gtfs-rt/vehiclepositions -121.020240 -120.248000 \n", - "\n", - " min_lat max_lat uuid \\\n", - "0 42.905244 43.051188 0fd41f26-6bc5-45f0-88ff-594bde1e8c24 \n", - "1 44.321804 44.420207 3196b95d-1eb2-4e05-bb36-7a2099d2348b \n", - "2 33.929498 34.075669 dbeb97e3-012d-4c23-93aa-a80696884340 \n", - "3 37.387656 37.431429 6885d33d-2f4e-4ac7-84a7-2e4747c723f0 \n", - "4 36.964503 37.521786 ea803659-c28f-43ac-8c01-619d767e0b1d \n", - "\n", - " tz_str epsg_code \n", - "0 America/Toronto 32617 \n", - "1 America/Toronto 32617 \n", - "2 America/Los_Angeles 32611 \n", - "3 America/Los_Angeles 32610 \n", - "4 America/Los_Angeles 32610 " + " mae mape rmse ex_var r_score \n", + "18 69.569825 0.224919 8.340853 0.775192 0.771584 \n", + "20 111.589098 0.259579 10.563574 0.777157 0.728569 \n", + "1 116.481377 0.237706 10.792654 0.793000 0.792672 \n", + "2 346.089513 0.434805 18.603481 0.686116 0.263355 \n", + "4 215.147161 0.263385 14.667896 0.743784 0.720664 \n", + "9 161.216318 0.260830 12.697099 0.550255 0.526094 \n", + "29 205.992501 0.269350 14.352439 0.683146 0.656156 \n", + "8 126.026914 0.285953 11.226171 0.785196 0.784499 \n", + "21 115.017353 0.258047 10.724614 0.777081 0.776323 \n", + "26 129.650576 0.325215 11.386421 0.680104 0.676025 " ] }, - "execution_count": 2, + "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "cleaned_sources = pd.read_csv(Path('..', 'data', 'cleaned_sources.csv'))\n", - "cleaned_sources.head()" + "metrics = {}\n", + "for key, value in res.items():\n", + " metrics[key] = model_utils.performance_metrics(value['labels'], value['preds'])\n", + "metrics_df = pd.DataFrame(metrics).T\n", + "metrics_df.index.name = 'uuid'\n", + "metrics_df = metrics_df.reset_index()\n", + "metrics_df = pd.merge(metrics_df, cleaned_sources, on='uuid')\n", + "metrics_df['experiment'] = 'not-tuned'\n", + "metrics_df['n_batches'] = '0_batches'\n", + "metrics_df[['provider','n_batches','experiment','mae', 'mape', 'rmse', 'ex_var', 'r_score']].sort_values(['provider', 'n_batches']).head(10)" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 46, "metadata": {}, "outputs": [ { @@ -241,127 +469,178 @@ " <thead>\n", " <tr style=\"text-align: right;\">\n", " <th></th>\n", - " <th>uuid</th>\n", - " <th>municipality</th>\n", " <th>provider</th>\n", - " <th>trip_id</th>\n", - " <th>file</th>\n", - " <th>locationtime</th>\n", - " <th>lat</th>\n", - " <th>lon</th>\n", - " <th>vehicle_id</th>\n", + " <th>n_batches</th>\n", + " <th>experiment</th>\n", + " <th>mae</th>\n", + " <th>mape</th>\n", + " <th>rmse</th>\n", + " <th>ex_var</th>\n", + " <th>r_score</th>\n", " </tr>\n", " </thead>\n", " <tbody>\n", " <tr>\n", - " <th>0</th>\n", - " <td>0fd41f26-6bc5-45f0-88ff-594bde1e8c24</td>\n", - " <td>London</td>\n", - " <td>London Transit Commission</td>\n", - " <td>2007915</td>\n", - " <td>2024_01_01</td>\n", - " <td>1704116300</td>\n", - " <td>43.018391</td>\n", - " <td>-81.243179</td>\n", - " <td>3155</td>\n", - " </tr>\n", - " <tr>\n", - " <th>1</th>\n", - " <td>0fd41f26-6bc5-45f0-88ff-594bde1e8c24</td>\n", - " <td>London</td>\n", - " <td>London Transit Commission</td>\n", - " <td>2007915</td>\n", - " <td>2024_01_01</td>\n", - " <td>1704116331</td>\n", - " <td>43.019211</td>\n", - " <td>-81.240379</td>\n", - " <td>3155</td>\n", - " </tr>\n", - " <tr>\n", - " <th>2</th>\n", - " <td>0fd41f26-6bc5-45f0-88ff-594bde1e8c24</td>\n", - " <td>London</td>\n", - " <td>London Transit Commission</td>\n", - " <td>2007915</td>\n", - " <td>2024_01_01</td>\n", - " <td>1704116361</td>\n", - " <td>43.020302</td>\n", - " <td>-81.236549</td>\n", - " <td>3155</td>\n", - " </tr>\n", - " <tr>\n", - " <th>3</th>\n", - " <td>0fd41f26-6bc5-45f0-88ff-594bde1e8c24</td>\n", - " <td>London</td>\n", - " <td>London Transit Commission</td>\n", - " <td>2007915</td>\n", - " <td>2024_01_01</td>\n", - " <td>1704116393</td>\n", - " <td>43.020969</td>\n", - " <td>-81.234413</td>\n", - " <td>3155</td>\n", - " </tr>\n", - " <tr>\n", - " <th>4</th>\n", - " <td>0fd41f26-6bc5-45f0-88ff-594bde1e8c24</td>\n", - " <td>London</td>\n", - " <td>London Transit Commission</td>\n", - " <td>2007915</td>\n", - " <td>2024_01_01</td>\n", - " <td>1704116423</td>\n", - " <td>43.022369</td>\n", - " <td>-81.229729</td>\n", - " <td>3155</td>\n", + " <th>151</th>\n", + " <td>Arlington Transit</td>\n", + " <td>128_batches</td>\n", + " <td>tuned</td>\n", + " <td>46.755706</td>\n", + " <td>0.179666</td>\n", + " <td>6.837814</td>\n", + " <td>0.886084</td>\n", + " <td>0.882567</td>\n", + " </tr>\n", + " <tr>\n", + " <th>148</th>\n", + " <td>Arlington Transit</td>\n", + " <td>16_batches</td>\n", + " <td>tuned</td>\n", + " <td>47.968289</td>\n", + " <td>0.182138</td>\n", + " <td>6.925914</td>\n", + " <td>0.886586</td>\n", + " <td>0.878768</td>\n", + " </tr>\n", + " <tr>\n", + " <th>144</th>\n", + " <td>Arlington Transit</td>\n", + " <td>1_batches</td>\n", + " <td>tuned</td>\n", + " <td>47.231036</td>\n", + " <td>0.179317</td>\n", + " <td>6.872484</td>\n", + " <td>0.885291</td>\n", + " <td>0.880821</td>\n", + " </tr>\n", + " <tr>\n", + " <th>145</th>\n", + " <td>Arlington Transit</td>\n", + " <td>2_batches</td>\n", + " <td>tuned</td>\n", + " <td>47.063368</td>\n", + " <td>0.175729</td>\n", + " <td>6.860275</td>\n", + " <td>0.888068</td>\n", + " <td>0.881468</td>\n", + " </tr>\n", + " <tr>\n", + " <th>149</th>\n", + " <td>Arlington Transit</td>\n", + " <td>32_batches</td>\n", + " <td>tuned</td>\n", + " <td>47.408175</td>\n", + " <td>0.180911</td>\n", + " <td>6.885359</td>\n", + " <td>0.887721</td>\n", + " <td>0.879983</td>\n", + " </tr>\n", + " <tr>\n", + " <th>146</th>\n", + " <td>Arlington Transit</td>\n", + " <td>4_batches</td>\n", + " <td>tuned</td>\n", + " <td>47.043689</td>\n", + " <td>0.178609</td>\n", + " <td>6.858840</td>\n", + " <td>0.887190</td>\n", + " <td>0.882800</td>\n", + " </tr>\n", + " <tr>\n", + " <th>150</th>\n", + " <td>Arlington Transit</td>\n", + " <td>64_batches</td>\n", + " <td>tuned</td>\n", + " <td>47.366756</td>\n", + " <td>0.177036</td>\n", + " <td>6.882351</td>\n", + " <td>0.886998</td>\n", + " <td>0.879714</td>\n", + " </tr>\n", + " <tr>\n", + " <th>147</th>\n", + " <td>Arlington Transit</td>\n", + " <td>8_batches</td>\n", + " <td>tuned</td>\n", + " <td>47.008838</td>\n", + " <td>0.176873</td>\n", + " <td>6.856299</td>\n", + " <td>0.888629</td>\n", + " <td>0.881684</td>\n", + " </tr>\n", + " <tr>\n", + " <th>167</th>\n", + " <td>BC Transit (Victoria Regional Transit System)</td>\n", + " <td>128_batches</td>\n", + " <td>tuned</td>\n", + " <td>58.076670</td>\n", + " <td>0.151077</td>\n", + " <td>7.620805</td>\n", + " <td>0.928405</td>\n", + " <td>0.927545</td>\n", + " </tr>\n", + " <tr>\n", + " <th>164</th>\n", + " <td>BC Transit (Victoria Regional Transit System)</td>\n", + " <td>16_batches</td>\n", + " <td>tuned</td>\n", + " <td>59.399674</td>\n", + " <td>0.152677</td>\n", + " <td>7.707118</td>\n", + " <td>0.925121</td>\n", + " <td>0.922669</td>\n", " </tr>\n", " </tbody>\n", "</table>\n", "</div>" ], "text/plain": [ - " uuid municipality \\\n", - "0 0fd41f26-6bc5-45f0-88ff-594bde1e8c24 London \n", - "1 0fd41f26-6bc5-45f0-88ff-594bde1e8c24 London \n", - "2 0fd41f26-6bc5-45f0-88ff-594bde1e8c24 London \n", - "3 0fd41f26-6bc5-45f0-88ff-594bde1e8c24 London \n", - "4 0fd41f26-6bc5-45f0-88ff-594bde1e8c24 London \n", - "\n", - " provider trip_id file locationtime lat \\\n", - "0 London Transit Commission 2007915 2024_01_01 1704116300 43.018391 \n", - "1 London Transit Commission 2007915 2024_01_01 1704116331 43.019211 \n", - "2 London Transit Commission 2007915 2024_01_01 1704116361 43.020302 \n", - "3 London Transit Commission 2007915 2024_01_01 1704116393 43.020969 \n", - "4 London Transit Commission 2007915 2024_01_01 1704116423 43.022369 \n", + " provider n_batches experiment \\\n", + "151 Arlington Transit 128_batches tuned \n", + "148 Arlington Transit 16_batches tuned \n", + "144 Arlington Transit 1_batches tuned \n", + "145 Arlington Transit 2_batches tuned \n", + "149 Arlington Transit 32_batches tuned \n", + "146 Arlington Transit 4_batches tuned \n", + "150 Arlington Transit 64_batches tuned \n", + "147 Arlington Transit 8_batches tuned \n", + "167 BC Transit (Victoria Regional Transit System) 128_batches tuned \n", + "164 BC Transit (Victoria Regional Transit System) 16_batches tuned \n", "\n", - " lon vehicle_id \n", - "0 -81.243179 3155 \n", - "1 -81.240379 3155 \n", - "2 -81.236549 3155 \n", - "3 -81.234413 3155 \n", - "4 -81.229729 3155 " + " mae mape rmse ex_var r_score \n", + "151 46.755706 0.179666 6.837814 0.886084 0.882567 \n", + "148 47.968289 0.182138 6.925914 0.886586 0.878768 \n", + "144 47.231036 0.179317 6.872484 0.885291 0.880821 \n", + "145 47.063368 0.175729 6.860275 0.888068 0.881468 \n", + "149 47.408175 0.180911 6.885359 0.887721 0.879983 \n", + "146 47.043689 0.178609 6.858840 0.887190 0.882800 \n", + "150 47.366756 0.177036 6.882351 0.886998 0.879714 \n", + "147 47.008838 0.176873 6.856299 0.888629 0.881684 \n", + "167 58.076670 0.151077 7.620805 0.928405 0.927545 \n", + "164 59.399674 0.152677 7.707118 0.925121 0.922669 " ] }, - "execution_count": 3, + "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "all_data = []\n", - "for uuid in cleaned_sources.uuid:\n", - " provider_path = Path('..', 'data', 'other_feeds', f\"{uuid}_realtime\")\n", - " available_files = [x.name for x in provider_path.glob('*.pkl')]\n", - " for fname in available_files:\n", - " data = pd.read_pickle(provider_path / fname)\n", - " data['uuid'] = uuid\n", - " all_data.append(data)\n", - "all_data = pd.concat(all_data)\n", - "all_data = pd.merge(cleaned_sources[['uuid', 'municipality', 'provider']], all_data, on='uuid')\n", - "all_data.head()" + "metrics = {}\n", + "for k_uuid, res_uuid in res_tuned.items():\n", + " for k_n_batches, res_n_batches in res_uuid.items():\n", + " metrics[(k_uuid, k_n_batches)] = model_utils.performance_metrics(res_n_batches['labels'], res_n_batches['preds'])\n", + "metrics_tuned_df = pd.DataFrame(metrics).T\n", + "metrics_tuned_df.index.names = ['uuid', 'n_batches']\n", + "metrics_tuned_df = metrics_tuned_df.reset_index()\n", + "metrics_tuned_df = pd.merge(metrics_tuned_df, cleaned_sources, on='uuid')\n", + "metrics_tuned_df['experiment'] = 'tuned'\n", + "metrics_tuned_df[['provider','n_batches','experiment','mae', 'mape', 'rmse', 'ex_var', 'r_score']].sort_values(['provider', 'n_batches']).head(10)" ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 65, "metadata": {}, "outputs": [ { @@ -385,680 +664,126 @@ " <thead>\n", " <tr style=\"text-align: right;\">\n", " <th></th>\n", - " <th>uuid</th>\n", - " <th>municipality</th>\n", - " <th>trip_id</th>\n", - " <th>file</th>\n", - " <th>locationtime</th>\n", - " <th>lat</th>\n", - " <th>lon</th>\n", - " <th>vehicle_id</th>\n", + " <th>mape</th>\n", " </tr>\n", " <tr>\n", - " <th>provider</th>\n", - " <th></th>\n", - " <th></th>\n", - " <th></th>\n", - " <th></th>\n", - " <th></th>\n", - " <th></th>\n", - " <th></th>\n", + " <th>n_batches</th>\n", " <th></th>\n", " </tr>\n", " </thead>\n", " <tbody>\n", " <tr>\n", - " <th>King County Metro</th>\n", - " <td>1035847</td>\n", - " <td>1035847</td>\n", - " <td>1035847</td>\n", - " <td>1035847</td>\n", - " <td>1035847</td>\n", - " <td>1035847</td>\n", - " <td>1035847</td>\n", - " <td>1035847</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Roma Servizi per la Mobilità</th>\n", - " <td>995006</td>\n", - " <td>995006</td>\n", - " <td>995006</td>\n", - " <td>995006</td>\n", - " <td>995006</td>\n", - " <td>995006</td>\n", - " <td>995006</td>\n", - " <td>995006</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Massachusetts Bay Transportation Authority (MBTA)</th>\n", - " <td>868739</td>\n", - " <td>868739</td>\n", - " <td>868739</td>\n", - " <td>868739</td>\n", - " <td>868739</td>\n", - " <td>868739</td>\n", - " <td>868739</td>\n", - " <td>868739</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Edmonton Transit System</th>\n", - " <td>765328</td>\n", - " <td>765328</td>\n", - " <td>765328</td>\n", - " <td>765328</td>\n", - " <td>765328</td>\n", - " <td>765328</td>\n", - " <td>765328</td>\n", - " <td>765328</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Capital Metro</th>\n", - " <td>736330</td>\n", - " <td>736330</td>\n", - " <td>736330</td>\n", - " <td>736330</td>\n", - " <td>736330</td>\n", - " <td>736330</td>\n", - " <td>736330</td>\n", - " <td>736330</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Port Authority of Allegheny County</th>\n", - " <td>688369</td>\n", - " <td>688369</td>\n", - " <td>688369</td>\n", - " <td>688369</td>\n", - " <td>688369</td>\n", - " <td>688369</td>\n", - " <td>688369</td>\n", - " <td>688369</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Metropolitan Atlanta Rapid Transit Authority (MARTA)</th>\n", - " <td>412615</td>\n", - " <td>412615</td>\n", - " <td>412615</td>\n", - " <td>412615</td>\n", - " <td>412615</td>\n", - " <td>412615</td>\n", - " <td>412615</td>\n", - " <td>412615</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Adelaide Metro</th>\n", - " <td>364396</td>\n", - " <td>364396</td>\n", - " <td>364396</td>\n", - " <td>364396</td>\n", - " <td>364396</td>\n", - " <td>364396</td>\n", - " <td>364396</td>\n", - " <td>364396</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Central Florida Regional Transit Authority (LYNX)</th>\n", - " <td>353448</td>\n", - " <td>353448</td>\n", - " <td>353448</td>\n", - " <td>353448</td>\n", - " <td>353448</td>\n", - " <td>353448</td>\n", - " <td>353448</td>\n", - " <td>353448</td>\n", + " <th>0_batches</th>\n", + " <td>0.110069</td>\n", " </tr>\n", " <tr>\n", - " <th>York Region Transit</th>\n", - " <td>300899</td>\n", - " <td>300899</td>\n", - " <td>300899</td>\n", - " <td>300899</td>\n", - " <td>300899</td>\n", - " <td>300899</td>\n", - " <td>300899</td>\n", - " <td>300899</td>\n", + " <th>1_batches</th>\n", + " <td>0.070764</td>\n", " </tr>\n", " <tr>\n", - " <th>Halifax Transit</th>\n", - " <td>292752</td>\n", - " <td>292752</td>\n", - " <td>292752</td>\n", - " <td>292752</td>\n", - " <td>292752</td>\n", - " <td>292752</td>\n", - " <td>292752</td>\n", - " <td>292752</td>\n", + " <th>2_batches</th>\n", + " <td>0.060201</td>\n", " </tr>\n", " <tr>\n", - " <th>Hamilton Street Railway</th>\n", - " <td>264707</td>\n", - " <td>264707</td>\n", - " <td>264707</td>\n", - " <td>264707</td>\n", - " <td>264707</td>\n", - " <td>264707</td>\n", - " <td>264707</td>\n", - " <td>264707</td>\n", - " </tr>\n", - " <tr>\n", - " <th>BC Transit (Victoria Regional Transit System)</th>\n", - " <td>225325</td>\n", - " <td>225325</td>\n", - " <td>225325</td>\n", - " <td>225325</td>\n", - " <td>225325</td>\n", - " <td>225325</td>\n", - " <td>225325</td>\n", - " <td>225325</td>\n", - " </tr>\n", - " <tr>\n", - " <th>London Transit Commission</th>\n", - " <td>216393</td>\n", - " <td>216393</td>\n", - " <td>216393</td>\n", - " <td>216393</td>\n", - " <td>216393</td>\n", - " <td>216393</td>\n", - " <td>216393</td>\n", - " <td>216393</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Metro St. Louis</th>\n", - " <td>197283</td>\n", - " <td>197283</td>\n", - " <td>197283</td>\n", - " <td>197283</td>\n", - " <td>197283</td>\n", - " <td>197283</td>\n", - " <td>197283</td>\n", - " <td>197283</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Nashville Metropolitan Transit Authority (Nashville MTA)</th>\n", - " <td>196899</td>\n", - " <td>196899</td>\n", - " <td>196899</td>\n", - " <td>196899</td>\n", - " <td>196899</td>\n", - " <td>196899</td>\n", - " <td>196899</td>\n", - " <td>196899</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Grand River Transit</th>\n", - " <td>193908</td>\n", - " <td>193908</td>\n", - " <td>193908</td>\n", - " <td>193908</td>\n", - " <td>193908</td>\n", - " <td>193908</td>\n", - " <td>193908</td>\n", - " <td>193908</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Regional Transportation District (RTD)</th>\n", - " <td>185767</td>\n", - " <td>185767</td>\n", - " <td>185767</td>\n", - " <td>185767</td>\n", - " <td>185767</td>\n", - " <td>185767</td>\n", - " <td>185767</td>\n", - " <td>185767</td>\n", - " </tr>\n", - " <tr>\n", - " <th>MiWay</th>\n", - " <td>179608</td>\n", - " <td>179608</td>\n", - " <td>179608</td>\n", - " <td>179608</td>\n", - " <td>179608</td>\n", - " <td>179608</td>\n", - " <td>179608</td>\n", - " <td>179608</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Rochester-Genesee Regional Transportation Authority (RGRTA)</th>\n", - " <td>110889</td>\n", - " <td>110889</td>\n", - " <td>110889</td>\n", - " <td>110889</td>\n", - " <td>110889</td>\n", - " <td>110889</td>\n", - " <td>110889</td>\n", - " <td>110889</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Metro Transit</th>\n", - " <td>110668</td>\n", - " <td>110668</td>\n", - " <td>110668</td>\n", - " <td>110668</td>\n", - " <td>110668</td>\n", - " <td>110668</td>\n", - " <td>110668</td>\n", - " <td>110668</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Intercity Transit</th>\n", - " <td>68357</td>\n", - " <td>68357</td>\n", - " <td>68357</td>\n", - " <td>68357</td>\n", - " <td>68357</td>\n", - " <td>68357</td>\n", - " <td>68357</td>\n", - " <td>68357</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Kingston Transit</th>\n", - " <td>62635</td>\n", - " <td>62635</td>\n", - " <td>62635</td>\n", - " <td>62635</td>\n", - " <td>62635</td>\n", - " <td>62635</td>\n", - " <td>62635</td>\n", - " <td>62635</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Saskatoon Transit</th>\n", - " <td>57982</td>\n", - " <td>57982</td>\n", - " <td>57982</td>\n", - " <td>57982</td>\n", - " <td>57982</td>\n", - " <td>57982</td>\n", - " <td>57982</td>\n", - " <td>57982</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Arlington Transit</th>\n", - " <td>51833</td>\n", - " <td>51833</td>\n", - " <td>51833</td>\n", - " <td>51833</td>\n", - " <td>51833</td>\n", - " <td>51833</td>\n", - " <td>51833</td>\n", - " <td>51833</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Duluth Transit</th>\n", - " <td>44031</td>\n", - " <td>44031</td>\n", - " <td>44031</td>\n", - " <td>44031</td>\n", - " <td>44031</td>\n", - " <td>44031</td>\n", - " <td>44031</td>\n", - " <td>44031</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Transit Authority of River City (TARC)</th>\n", - " <td>38799</td>\n", - " <td>38799</td>\n", - " <td>38799</td>\n", - " <td>38799</td>\n", - " <td>38799</td>\n", - " <td>38799</td>\n", - " <td>38799</td>\n", - " <td>38799</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Big Blue Bus</th>\n", - " <td>35066</td>\n", - " <td>35066</td>\n", - " <td>35066</td>\n", - " <td>35066</td>\n", - " <td>35066</td>\n", - " <td>35066</td>\n", - " <td>35066</td>\n", - " <td>35066</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Thunder Bay Transit</th>\n", - " <td>34248</td>\n", - " <td>34248</td>\n", - " <td>34248</td>\n", - " <td>34248</td>\n", - " <td>34248</td>\n", - " <td>34248</td>\n", - " <td>34248</td>\n", - " <td>34248</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Barrie Transit</th>\n", - " <td>33623</td>\n", - " <td>33623</td>\n", - " <td>33623</td>\n", - " <td>33623</td>\n", - " <td>33623</td>\n", - " <td>33623</td>\n", - " <td>33623</td>\n", - " <td>33623</td>\n", - " </tr>\n", - " <tr>\n", - " <th>TransLink Sunbus Cairns</th>\n", - " <td>33248</td>\n", - " <td>33248</td>\n", - " <td>33248</td>\n", - " <td>33248</td>\n", - " <td>33248</td>\n", - " <td>33248</td>\n", - " <td>33248</td>\n", - " <td>33248</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Stanislaus Regional Transit Authority (StanRTA)</th>\n", - " <td>15619</td>\n", - " <td>15619</td>\n", - " <td>15619</td>\n", - " <td>15619</td>\n", - " <td>15619</td>\n", - " <td>15619</td>\n", - " <td>15619</td>\n", - " <td>15619</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Votran</th>\n", - " <td>12285</td>\n", - " <td>12285</td>\n", - " <td>12285</td>\n", - " <td>12285</td>\n", - " <td>12285</td>\n", - " <td>12285</td>\n", - " <td>12285</td>\n", - " <td>12285</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Merced County Transit (The Bus)</th>\n", - " <td>9455</td>\n", - " <td>9455</td>\n", - " <td>9455</td>\n", - " <td>9455</td>\n", - " <td>9455</td>\n", - " <td>9455</td>\n", - " <td>9455</td>\n", - " <td>9455</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Beach Cities Transit</th>\n", - " <td>5512</td>\n", - " <td>5512</td>\n", - " <td>5512</td>\n", - " <td>5512</td>\n", - " <td>5512</td>\n", - " <td>5512</td>\n", - " <td>5512</td>\n", - " <td>5512</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Commerce Municipal Bus Lines</th>\n", - " <td>3610</td>\n", - " <td>3610</td>\n", - " <td>3610</td>\n", - " <td>3610</td>\n", - " <td>3610</td>\n", - " <td>3610</td>\n", - " <td>3610</td>\n", - " <td>3610</td>\n", - " </tr>\n", - " <tr>\n", - " <th>Port Alberni</th>\n", - " <td>3290</td>\n", - " <td>3290</td>\n", - " <td>3290</td>\n", - " <td>3290</td>\n", - " <td>3290</td>\n", - " <td>3290</td>\n", - " <td>3290</td>\n", - " <td>3290</td>\n", + " <th>4_batches</th>\n", + " <td>0.052487</td>\n", " </tr>\n", " </tbody>\n", "</table>\n", "</div>" ], "text/plain": [ - " uuid municipality \\\n", - "provider \n", - "King County Metro 1035847 1035847 \n", - "Roma Servizi per la Mobilità 995006 995006 \n", - "Massachusetts Bay Transportation Authority (MBTA) 868739 868739 \n", - "Edmonton Transit System 765328 765328 \n", - "Capital Metro 736330 736330 \n", - "Port Authority of Allegheny County 688369 688369 \n", - "Metropolitan Atlanta Rapid Transit Authority (M... 412615 412615 \n", - "Adelaide Metro 364396 364396 \n", - "Central Florida Regional Transit Authority (LYNX) 353448 353448 \n", - "York Region Transit 300899 300899 \n", - "Halifax Transit 292752 292752 \n", - "Hamilton Street Railway 264707 264707 \n", - "BC Transit (Victoria Regional Transit System) 225325 225325 \n", - "London Transit Commission 216393 216393 \n", - "Metro St. Louis 197283 197283 \n", - "Nashville Metropolitan Transit Authority (Nashv... 196899 196899 \n", - "Grand River Transit 193908 193908 \n", - "Regional Transportation District (RTD) 185767 185767 \n", - "MiWay 179608 179608 \n", - "Rochester-Genesee Regional Transportation Autho... 110889 110889 \n", - "Metro Transit 110668 110668 \n", - "Intercity Transit 68357 68357 \n", - "Kingston Transit 62635 62635 \n", - "Saskatoon Transit 57982 57982 \n", - "Arlington Transit 51833 51833 \n", - "Duluth Transit 44031 44031 \n", - "Transit Authority of River City (TARC) 38799 38799 \n", - "Big Blue Bus 35066 35066 \n", - "Thunder Bay Transit 34248 34248 \n", - "Barrie Transit 33623 33623 \n", - "TransLink Sunbus Cairns 33248 33248 \n", - "Stanislaus Regional Transit Authority (StanRTA) 15619 15619 \n", - "Votran 12285 12285 \n", - "Merced County Transit (The Bus) 9455 9455 \n", - "Beach Cities Transit 5512 5512 \n", - "Commerce Municipal Bus Lines 3610 3610 \n", - "Port Alberni 3290 3290 \n", - "\n", - " trip_id file \\\n", - "provider \n", - "King County Metro 1035847 1035847 \n", - "Roma Servizi per la Mobilità 995006 995006 \n", - "Massachusetts Bay Transportation Authority (MBTA) 868739 868739 \n", - "Edmonton Transit System 765328 765328 \n", - "Capital Metro 736330 736330 \n", - "Port Authority of Allegheny County 688369 688369 \n", - "Metropolitan Atlanta Rapid Transit Authority (M... 412615 412615 \n", - "Adelaide Metro 364396 364396 \n", - "Central Florida Regional Transit Authority (LYNX) 353448 353448 \n", - "York Region Transit 300899 300899 \n", - "Halifax Transit 292752 292752 \n", - "Hamilton Street Railway 264707 264707 \n", - "BC Transit (Victoria Regional Transit System) 225325 225325 \n", - "London Transit Commission 216393 216393 \n", - "Metro St. Louis 197283 197283 \n", - "Nashville Metropolitan Transit Authority (Nashv... 196899 196899 \n", - "Grand River Transit 193908 193908 \n", - "Regional Transportation District (RTD) 185767 185767 \n", - "MiWay 179608 179608 \n", - "Rochester-Genesee Regional Transportation Autho... 110889 110889 \n", - "Metro Transit 110668 110668 \n", - "Intercity Transit 68357 68357 \n", - "Kingston Transit 62635 62635 \n", - "Saskatoon Transit 57982 57982 \n", - "Arlington Transit 51833 51833 \n", - "Duluth Transit 44031 44031 \n", - "Transit Authority of River City (TARC) 38799 38799 \n", - "Big Blue Bus 35066 35066 \n", - "Thunder Bay Transit 34248 34248 \n", - "Barrie Transit 33623 33623 \n", - "TransLink Sunbus Cairns 33248 33248 \n", - "Stanislaus Regional Transit Authority (StanRTA) 15619 15619 \n", - "Votran 12285 12285 \n", - "Merced County Transit (The Bus) 9455 9455 \n", - "Beach Cities Transit 5512 5512 \n", - "Commerce Municipal Bus Lines 3610 3610 \n", - "Port Alberni 3290 3290 \n", - "\n", - " locationtime lat \\\n", - "provider \n", - "King County Metro 1035847 1035847 \n", - "Roma Servizi per la Mobilità 995006 995006 \n", - "Massachusetts Bay Transportation Authority (MBTA) 868739 868739 \n", - "Edmonton Transit System 765328 765328 \n", - "Capital Metro 736330 736330 \n", - "Port Authority of Allegheny County 688369 688369 \n", - "Metropolitan Atlanta Rapid Transit Authority (M... 412615 412615 \n", - "Adelaide Metro 364396 364396 \n", - "Central Florida Regional Transit Authority (LYNX) 353448 353448 \n", - "York Region Transit 300899 300899 \n", - "Halifax Transit 292752 292752 \n", - "Hamilton Street Railway 264707 264707 \n", - "BC Transit (Victoria Regional Transit System) 225325 225325 \n", - "London Transit Commission 216393 216393 \n", - "Metro St. Louis 197283 197283 \n", - "Nashville Metropolitan Transit Authority (Nashv... 196899 196899 \n", - "Grand River Transit 193908 193908 \n", - "Regional Transportation District (RTD) 185767 185767 \n", - "MiWay 179608 179608 \n", - "Rochester-Genesee Regional Transportation Autho... 110889 110889 \n", - "Metro Transit 110668 110668 \n", - "Intercity Transit 68357 68357 \n", - "Kingston Transit 62635 62635 \n", - "Saskatoon Transit 57982 57982 \n", - "Arlington Transit 51833 51833 \n", - "Duluth Transit 44031 44031 \n", - "Transit Authority of River City (TARC) 38799 38799 \n", - "Big Blue Bus 35066 35066 \n", - "Thunder Bay Transit 34248 34248 \n", - "Barrie Transit 33623 33623 \n", - "TransLink Sunbus Cairns 33248 33248 \n", - "Stanislaus Regional Transit Authority (StanRTA) 15619 15619 \n", - "Votran 12285 12285 \n", - "Merced County Transit (The Bus) 9455 9455 \n", - "Beach Cities Transit 5512 5512 \n", - "Commerce Municipal Bus Lines 3610 3610 \n", - "Port Alberni 3290 3290 \n", - "\n", - " lon vehicle_id \n", - "provider \n", - "King County Metro 1035847 1035847 \n", - "Roma Servizi per la Mobilità 995006 995006 \n", - "Massachusetts Bay Transportation Authority (MBTA) 868739 868739 \n", - "Edmonton Transit System 765328 765328 \n", - "Capital Metro 736330 736330 \n", - "Port Authority of Allegheny County 688369 688369 \n", - "Metropolitan Atlanta Rapid Transit Authority (M... 412615 412615 \n", - "Adelaide Metro 364396 364396 \n", - "Central Florida Regional Transit Authority (LYNX) 353448 353448 \n", - "York Region Transit 300899 300899 \n", - "Halifax Transit 292752 292752 \n", - "Hamilton Street Railway 264707 264707 \n", - "BC Transit (Victoria Regional Transit System) 225325 225325 \n", - "London Transit Commission 216393 216393 \n", - "Metro St. Louis 197283 197283 \n", - "Nashville Metropolitan Transit Authority (Nashv... 196899 196899 \n", - "Grand River Transit 193908 193908 \n", - "Regional Transportation District (RTD) 185767 185767 \n", - "MiWay 179608 179608 \n", - "Rochester-Genesee Regional Transportation Autho... 110889 110889 \n", - "Metro Transit 110668 110668 \n", - "Intercity Transit 68357 68357 \n", - "Kingston Transit 62635 62635 \n", - "Saskatoon Transit 57982 57982 \n", - "Arlington Transit 51833 51833 \n", - "Duluth Transit 44031 44031 \n", - "Transit Authority of River City (TARC) 38799 38799 \n", - "Big Blue Bus 35066 35066 \n", - "Thunder Bay Transit 34248 34248 \n", - "Barrie Transit 33623 33623 \n", - "TransLink Sunbus Cairns 33248 33248 \n", - "Stanislaus Regional Transit Authority (StanRTA) 15619 15619 \n", - "Votran 12285 12285 \n", - "Merced County Transit (The Bus) 9455 9455 \n", - "Beach Cities Transit 5512 5512 \n", - "Commerce Municipal Bus Lines 3610 3610 \n", - "Port Alberni 3290 3290 " + " mape\n", + "n_batches \n", + "0_batches 0.110069\n", + "1_batches 0.070764\n", + "2_batches 0.060201\n", + "4_batches 0.052487" ] }, - "execution_count": 4, + "execution_count": 65, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "all_data.groupby('provider').count().sort_values('uuid', ascending=False)" + "plot_df.groupby('n_batches')[['mape']].std()" ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "GRU(\n", - " (loss_fn): MSELoss()\n", - " (min_em): MinuteEmbedding(\n", - " (em): Embedding(1440, 48)\n", - " )\n", - " (day_em): DayEmbedding(\n", - " (em): Embedding(7, 4)\n", - " )\n", - " (rnn): GRU(4, 64, num_layers=2, dropout=0.05)\n", - " (feature_extract): Linear(in_features=116, out_features=1, bias=True)\n", - " (feature_extract_activation): ReLU()\n", - ")" + "<Axes: xlabel='MAPE', ylabel='Count'>" ] }, - "execution_count": 4, + "execution_count": 57, "metadata": {}, "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGwCAYAAACHJU4LAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAACXyklEQVR4nOzdd3gUVdvH8e9syWZTNr2SSu8dITRRilJsiNgBxYagKBYeUBAsgAWxgx2xob6PKIigdKRJDb0IAgHS66ZunfePyEoektACu4H7c117JTt7ZuaeTWB/OXPmjKKqqooQQgghhAfSuLsAIYQQQoiqSFARQgghhMeSoCKEEEIIjyVBRQghhBAeS4KKEEIIITyWBBUhhBBCeCwJKkIIIYTwWDp3F3AhnE4nqamp+Pv7oyiKu8sRQgghxFlQVZXCwkKio6PRaKrvM6nVQSU1NZXY2Fh3lyGEEEKI83Ds2DFiYmKqbVOrg4q/vz9QfqAmk8nN1QghhBDibJjNZmJjY12f49Wp1UHl5Okek8kkQUUIIYSoZc5m2IYMphVCCCGEx5KgIoQQQgiPJUFFCCGEEB6rVo9REUIIUfMcDgc2m83dZYhaTK/Xo9Vqa2RbElSEEEIA5XNbpKenk5+f7+5SxGUgMDCQyMjIC57nTIKKEEIIAFdICQ8Px8fHRybSFOdFVVVKSkrIzMwEICoq6oK2J0FFCCEEDofDFVJCQkLcXY6o5YxGIwCZmZmEh4df0GkgGUwrhBDCNSbFx8fHzZWIy8XJ36ULHe8kQUUIIYSLnO4RNaWmfpckqAghhBDCY0lQEUIIIYTHcmtQcTgcTJgwgcTERIxGI/Xq1eOll15CVVV3liWEEMINEhISeOutt9xdhlutXLkSRVEu+iXiPXr04Iknnrio+6gpbg0qr776KjNnzuS9995j7969vPrqq7z22mu8++677ixLCCHE/1AUpdrHpEmTLngfmzZt4qGHHrrwYs9g+/bt3HjjjYSHh+Pt7U1CQgK3336763La2s7hcDBt2jQaN26M0WgkODiYjh078sknn7ja/Pjjj7z00kturPLsufXy5HXr1nHTTTfRv39/oDxNf/vtt2zcuNGdZQkhhPgfaWlpru+/++47Jk6cyP79+13L/Pz8LngfYWFhF7yNM8nKyqJnz54MGDCA3377jcDAQI4cOcL8+fMpLi6+6Pu/FCZPnsyHH37Ie++9R/v27TGbzWzevJm8vDxXm+DgYDdWeG7c2qPSuXNnli1bxoEDB4DylLtmzRr69u1baXuLxYLZbK7wEEIIcfFFRka6HgEBASiK4no+a9YsunbtWqH9W2+9RUJCguv5sGHDuPnmm3njjTeIiooiJCSEkSNHVrh09X9P/SiKwieffMItt9yCj48PDRo0YP78+RX2M3/+fBo0aIC3tzfXXHMNX3zxRbWnTtauXUtBQQGffPIJbdq0ITExkWuuuYYZM2aQmJgIlPdIDB8+3DUsoVGjRrz99tsVtnPyeKZMmUJERASBgYG8+OKL2O12nnnmGYKDg4mJieHzzz93rXPkyBEURWHu3Ll07twZb29vmjdvzqpVq6p979esWUO3bt0wGo3Exsby+OOPVxuq5s+fz6OPPsptt91GYmIirVq1Yvjw4Tz99NOuNqee+jl5uul/H8OGDXO1//nnn2nbti3e3t7UrVuXyZMnY7fbq627prg1qPznP//hjjvuoHHjxuj1etq0acMTTzzB3XffXWn7qVOnEhAQ4HrExsZe4oqFJ2vXugXREaHVPtq1buHuMoW4Yq1YsYJDhw6xYsUKvvjiC2bPns3s2bOrXWfy5MkMHjyYHTt20K9fP+6++25yc3MBOHz4MIMGDeLmm29m+/btPPzwwzz33HPVbi8yMhK73c68efOqHA/pdDqJiYnhhx9+YM+ePUycOJHx48fz/fffV2i3fPlyUlNTWb16NW+++SYvvPACAwYMICgoiD///JNHHnmEhx9+mOPHj1dY75lnnuGpp55i27ZtJCUlccMNN5CTk1NpLYcOHeL666/n1ltvZceOHXz33XesWbOGUaNGVXuMy5cvJysrq9r34qTOnTuTlpbmeixfvhxvb2+6d+8OwB9//MGQIUMYPXo0e/bs4cMPP2T27Nm88sorZ7X9C6a60bfffqvGxMSo3377rbpjxw51zpw5anBwsDp79uxK25eVlakFBQWux7Fjx1RALSgouMSVC08UFR6iqsunVPuICg9xd5lCeKTS0lJ1z549amlp6Rnbfv7552pAQIDr+QsvvKC2atWqQpsZM2ao8fHxrudDhw5V4+PjVbvd7lp22223qbfffrvreXx8vDpjxgzXc0B9/vnnXc+LiopUQF20aJGqqqo6duxYtXnz5hX2+9xzz6mAmpeXV2X948ePV3U6nRocHKxef/316muvvaamp6dXe8wjR45Ub7311tOOx+FwuJY1atRI7datm+u53W5XfX191W+//VZVVVU9fPiwCqjTpk1ztbHZbGpMTIz66quvqqqqqitWrKhQ//Dhw9WHHnqoQi1//PGHqtFoqvxZ7d69W23SpImq0WjUFi1aqA8//LD666+/Vmhz9dVXq6NHjz5t3ezsbLVu3brqo48+6lrWs2dPdcqUKRXaffnll2pUVFSl+z+put+pgoKCs/78dmuPyjPPPOPqVWnRogX33nsvTz75JFOnTq20vcFgwGQyVXgIIYSoHZo1a1ZhKvWoqKgzDmBt2bKl63tfX19MJpNrnf3799OhQ4cK7a+66qoz1vHKK6+Qnp7OrFmzaNasGbNmzaJx48bs3LnT1eb999+nXbt2hIWF4efnx0cffURKSsppx6PR/PsxGhERQYsW//baarVaQkJCTjvGpKQk1/c6nY727duzd+/eSmvdvn07s2fPxs/Pz/W47rrrcDqdHD58uNJ1mjZtyq5du9iwYQP3338/mZmZ3HDDDTzwwAPVvi82m41bb72V+Pj4Cqe6tm/fzosvvlihhgcffJC0tDRKSkqq3WZNcOtg2pKSkgo/ZCj/wTqdTjdVJIQQ4lxpNJrTTqNUNm26Xq+v8FxRlDP+f38+65yNkJAQbrvtNm677TamTJlCmzZteOONN/jiiy+YO3cuTz/9NNOnTycpKQl/f39ef/11/vzzzzPWVtP1FhUV8fDDD/P444+f9lpcXFyV62k0Gjp06ECHDh144okn+Oqrr7j33nt57rnnXGNx/teIESM4duwYGzduRKf7Nx4UFRUxefJkBg4ceNo63t7e53FU58atQeWGG27glVdeIS4ujmbNmrFt2zbefPNN7r//fneWJYQQ4hyEhYWRnp6OqqquadOTk5Mv+n4bNWrEr7/+WmHZpk2bznk7Xl5e1KtXzzVAde3atXTu3JlHH33U1ebQoUMXVuwpNmzY4Br/Ybfb2bJlS5VjTtq2bcuePXuoX7/+Be2zadOmAFUOwn3zzTf5/vvvWbdu3Wk3pWzbti379++/4BrOl1uDyrvvvsuECRN49NFHyczMJDo6mocffpiJEye6sywhhBDnoEePHmRlZfHaa68xaNAgFi9ezKJFiy766fmHH36YN998k7FjxzJ8+HCSk5Ndg3Orus/ML7/8wty5c7njjjto2LAhqqqyYMECfv31V9cVOg0aNGDOnDn89ttvJCYm8uWXX7Jp06YqeyLO1fvvv0+DBg1o0qQJM2bMIC8vr8o/0MeOHUunTp0YNWoUDzzwAL6+vuzZs4clS5bw3nvvVbrOoEGD6NKlC507dyYyMpLDhw8zbtw4GjZsSOPGjU9rv3TpUp599lnef/99QkNDSU9PB8rvgBwQEMDEiRMZMGAAcXFxDBo0CI1Gw/bt29m1axcvv/xyjbwn1XHrGBV/f3/eeustjh49SmlpKYcOHeLll1/Gy8vLnWUJIYQ4B02aNOGDDz7g/fffp1WrVmzcuLHCpbAXS2JiIv/3f//Hjz/+SMuWLZk5c6brqh+DwVDpOk2bNsXHx4ennnqK1q1b06lTJ77//ns++eQT7r33XqA8AA0cOJDbb7+djh07kpOTU6F35UJNmzaNadOm0apVK9asWcP8+fMJDQ2ttG3Lli1ZtWoVBw4coFu3brRp04aJEycSHR1d5favu+46FixYwA033EDDhg0ZOnQojRs35vfff69wSuekNWvW4HA4eOSRR4iKinI9Ro8e7dreL7/8wu+//06HDh3o1KkTM2bMID4+vmbekDNQ1P89sViLmM1mAgICKCgokIG1guiIUFLnPlV9mzumk5qRfYkqEqL2KCsr4/DhwyQmJl6ScQcXyyuvvMKsWbM4duyYu0s5zZEjR0hMTGTbtm20bt3a3eVcdNX9Tp3L57dbT/0IIYQQF+KDDz6gQ4cOhISEsHbtWl5//fVq5xgRtY8EFSGEELXWX3/9xcsvv0xubi5xcXE89dRTjBs3zt1liRokQUUIIUStNWPGDGbMmOHuMs5KQkJClbPhiqq5dTCtEEIIIUR1JKgIIYQQwmNJUBFCCCGEx5KgIoQQQgiPJUFFCCGEEB5LrvoRQghRrZSUFLKzL91EiaGhodXecO9yMmzYMPLz8/npp5/cXYrHkqAihBCiSikpKTRu0oTSkpJLtk+jjw/79u4967AybNgwvvjiC6ZOncp//vMf1/KffvqJW2655bwvCe7RowerVq2q8vWrr76alStXnte2T3r77bflkuUzkKAihBCiStnZ2ZSWlHD32NeJiKt30feXkXKIr199huzs7HPqVfH29ubVV1/l4YcfJigoqEZq+fHHH7FarQAcO3aMq666iqVLl9KsWTOAGrkvXUBAwAVv43InY1SEEEKcUURcPWIaNLvoj/MNQ7169SIyMpKpU6dW2+6///0vzZo1w2AwkJCQwPTp06tsGxwcTGRkJJGRkYSFhQEQEhJCZGQk+/btIyQkhPz8fFf75ORkFEXhyJEjAMyePZvAwEB+++03mjRpgp+fH9dffz1paWmudYYNG8bNN9/set6jRw8ef/xxnn32Wdf+J02aVKGuffv20bVrV7y9vWnatClLly5FUZTL9vSRBBUhhBC1nlarZcqUKbz77rscP3680jZbtmxh8ODB3HHHHezcuZNJkyYxYcIEZs+efdHqKikp4Y033uDLL79k9erVpKSknPHO0l988QW+vr78+eefvPbaa7z44ossWbIEAIfDwc0334yPjw9//vknH330keuO0ZcrCSpCCCEuC7fccgutW7fmhRdeqPT1N998k549ezJhwgQaNmzIsGHDGDVqFK+//vpFq8lmszFr1izat29P27ZtGTVqFMuWLat2nZYtW/LCCy/QoEEDhgwZQvv27V3rLFmyhEOHDjFnzhxatWpF165deeWVVy5a/Z5AgooQQojLxquvvsoXX3zB3r17T3tt7969dOnSpcKyLl268Ndff+FwOC5KPT4+PtSr9+/prKioKDIzM6tdp2XLlhWen7rO/v37iY2NJTIy0vX6VVddVYMVex4JKkIIIS4b3bt357rrrrvod1DWaMo/Pk+9Ysdms53WTq/XV3iuKMoZr/KpbB2n03m+pdZ6ctWPEEKIy8q0adNo3bo1jRo1qrC8SZMmrF27tsKytWvX0rBhQ7Ra7Tnt4+Tg2rS0NNdVRsnJyedf9Flq1KgRx44dIyMjg4iICAA2bdp00ffrTtKjIoQQ4rLSokUL7r77bt55550Ky5966imWLVvGSy+9xIEDB/jiiy947733zji4tTL169cnNjaWSZMm8ddff7Fw4cJqryCqKb1796ZevXoMHTqUHTt2sHbtWp5//nmgvOflciQ9KkIIIc4oI+VQrdrPiy++yHfffVdhWdu2bfn++++ZOHEiL730ElFRUbz44osMGzbsnLev1+v59ttvGTFiBC1btqRDhw68/PLL3HbbbTVSf1W0Wi0//fQTDzzwAB06dKBu3bq8/vrr3HDDDXh7e1/UfbuLotbiKfHMZjMBAQEUFBRgMpncXY5ws+iIUFLnPlV9mzumk5px6aYCF6K2KCsr4/DhwyQmJlb4wKsNM9Ne6dauXUvXrl05ePBghYG77lbV7xSc2+e39KgIIYSoUlxcHPv27pV7/XiQefPm4efnR4MGDTh48CCjR4+mS5cuHhVSapIEFSGEENWKi4uT4OBBCgsLGTt2LCkpKYSGhtKrV69LMj7GXSSoCCGEELXIkCFDGDJkiLvLuGTkqh8hhBBCeCwJKkIIIYTwWBJUhBBCCOGxJKgIIYQQwmNJUBFCCCGEx5KgIoQQQgiPJZcnCyGEqFZKSopM+HaZmjRpEj/99NMluaHi+ZKgIoQQokopKSk0adKYkpLSS7ZPHx8je/fuO+uwsnr1al5//XW2bNlCWloa8+bN4+abb67QRlVVXnjhBT7++GPy8/Pp0qULM2fOpEGDBq42ubm5PPbYYyxYsACNRsOtt97K22+/jZ+fn6vNjh07GDlyJJs2bSIsLIzHHnuMZ599ttr65s2bx6uvvsrevXtxOp3ExcXRu3dv3nrrrbN+T65kElSEEEJUKTs7m5KSUr4aP5gmcWEXfX97U7K4Z8r3ZGdnn3VQKS4uplWrVtx///0MHDiw0javvfYa77zzDl988QWJiYlMmDCB6667jj179rjuQ3P33XeTlpbGkiVLsNls3HfffTz00EN88803QPn9afr06UOvXr2YNWsWO3fu5P777ycwMJCHHnqo0v0uW7aM22+/nVdeeYUbb7wRRVHYs2cPS5YsOY9358rk1qCSkJDA0aNHT1v+6KOP8v7777uhIiGEEJVpEhdG24Z13F1Gpfr27Uvfvn2rfF1VVd566y2ef/55brrpJgDmzJlDREQEP/30E3fccQd79+5l8eLFbNq0ifbt2wPw7rvv0q9fP9544w2io6P5+uuvsVqtfPbZZ3h5edGsWTOSk5N58803qwwqCxYsoEuXLjzzzDOuZQ0bNqzQ43Po0CHGjBnDhg0bKC4upkmTJkydOpVevXq52iQkJPDAAw9w4MABfvzxR0JCQnj33XdJSkrigQceYNmyZdStW5fPPvvMVf/s2bN54oknmD17Ns888wzHjh3j6quv5pNPPiE2NrbK9+uTTz5h+vTpHD58mISEBB5//HEeffRRAKxWK2PGjOG///0veXl5RERE8MgjjzBu3Lgz/JTOn1sH027atIm0tDTX42TCvNi3yRZCCHHlOHz4MOnp6RU++AMCAujYsSPr168HYP369QQGBro+5AF69eqFRqPhzz//dLXp3r07Xl5erjbXXXcd+/fvJy8vr9J9R0ZGsnv3bnbt2lVlfUVFRfTr149ly5axbds2rr/+em644QZSUlIqtJsxYwZdunRh27Zt9O/fn3vvvZchQ4Zwzz33sHXrVurVq8eQIUNQVdW1TklJCa+88gpz5sxh7dq15Ofnc8cdd1RZy9dff83EiRN55ZVX2Lt3L1OmTGHChAl88cUXALzzzjvMnz+f77//nv379/P111+TkJBQ5fZqglt7VMLCKnYjTps2jXr16nH11Ve7qSIhhBCXm/T0dAAiIiIqLI+IiHC9lp6eTnh4eIXXdTodwcHBFdokJiaeto2TrwUFBZ2278cee4w//viDFi1aEB8fT6dOnejTpw933303BoMBgFatWtGqVSvXOi+99BLz5s1j/vz5jBo1yrW8X79+PPzwwwBMnDiRmTNn0qFDB9cf92PHjiUpKYmMjAwiIyMBsNlsvPfee3Ts2BGAL774giZNmrBx40auuuqq0+p94YUXmD59uusUWmJiInv27OHDDz9k6NChpKSk0KBBA7p27YqiKMTHx1fxrtccj7k82Wq18tVXX3H//fejKEqlbSwWC2azucJDCCGE8FS+vr4sXLiQgwcP8vzzz+Pn58dTTz3FVVddRUlJCVDeo/L000/TpEkTAgMD8fPzY+/evaf1qLRs2dL1/cmA1KJFi9OWZWZmupbpdDo6dOjget64cWMCAwPZu3fvabUWFxdz6NAhhg8fjp+fn+vx8ssvc+jQIQCGDRtGcnIyjRo14vHHH+f333+/0LfojDwmqPz000/k5+czbNiwKttMnTqVgIAA16O6c2xCCCEE4OpdyMjIqLD81J6HyMjICh/wAHa7ndzc3AptKtvGqfuoSr169XjggQf45JNP2Lp1K3v27OG7774D4Omnn2bevHlMmTKFP/74g+TkZFq0aIHVaq2wDb1e7/r+5B/0lS1zOp3V1lKVoqIiAD7++GOSk5Ndj127drFhwwYA2rZty+HDh3nppZcoLS1l8ODBDBo06Lz2d7Y8Jqh8+umn9O3bl+jo6CrbjBs3joKCAtfj2LFjl7BCIYQQtVFiYiKRkZEsW7bMtcxsNvPnn3+SlJQEQFJSEvn5+WzZssXVZvny5TidTtdpk6SkJFavXo3NZnO1WbJkCY0aNar0tE9VEhIS8PHxobi4GIC1a9cybNgwbrnlFlq0aEFkZCRHjhy5kEN2sdvtbN682fV8//795Ofn06RJk9PaRkREEB0dzd9//039+vUrPE495WUymbj99tv5+OOP+e677/jvf/9Lbm5ujdRbGY+4PPno0aMsXbqUH3/8sdp2BoPBdU5PCCGEgPKegIMHD7qeHz58mOTkZIKDg4mLi0NRFJ544glefvllGjRo4Lo8OTo62nX1TZMmTbj++ut58MEHmTVrFjabjVGjRnHHHXe4/oC+6667mDx5MsOHD2fs2LHs2rWLt99+mxkzZlRZ26RJkygpKaFfv37Ex8eTn5/PO++8g81mo3fv3gA0aNCAH3/8kRtuuAFFUZgwYcJ594r8L71ez2OPPcY777yDTqdj1KhRdOrUqdLxKQCTJ0/m8ccfJyAggOuvvx6LxcLmzZvJy8tjzJgxvPnmm0RFRdGmTRs0Gg0//PADkZGRBAYG1ki9lfGIoPL5558THh5O//793V2KEEKISuxNyfLY/WzevJlrrrnG9XzMmDEADB06lNmzZwPw7LPPUlxczEMPPUR+fj5du3Zl8eLFrjlUoPyKl1GjRtGzZ0/XhG/vvPOO6/WAgAB+//13Ro4cSbt27QgNDWXixIlVXpoMcPXVV/P+++8zZMgQMjIyCAoKok2bNvz+++80atQIgDfffJP777+fzp07ExoaytixY2tsDKaPjw9jx47lrrvu4sSJE3Tr1o1PP/20yvYPPPAAPj4+vP766zzzzDP4+vrSokULnnjiCQD8/f157bXX+Ouvv9BqtXTo0IFff/0VjebinaBR1FOvY3IDp9NJYmIid955J9OmTTundc1mMwEBARQUFGAymS5ShaK2iI4IJXXuU9W3uWM6qRmXbipwIWqLsrIyDh8+TGJiYoUP79owM62o3Ml5VPLz892y/6p+p+DcPr/d3qOydOlSUlJSuP/++91dihBCiP8RFxfH3r375F4/wm3cHlT69OmDmzt1hBBCVCMuLk6Cg3Abj7nqRwghhBA1Z9iwYW477VOTJKgIIYQQwmNJUBFCCCGEx5KgIoQQQgiPJUFFCCGEEB5LgooQQgghPJYEFSGEEEJ4LLfPoyKEEMKzpaSkyIRvl6mTlzD/9NNP7i6lShJUhBBCVCklJYXGTRpTegmn0Df6GNl3DlPoT506lR9//JF9+/ZhNBrp3Lkzr776quteOlA+nftTTz3F3LlzsVgsXHfddXzwwQdERES42qSkpDBixAhWrFiBn58fQ4cOZerUqeh0/35Urly5kjFjxrB7925iY2N5/vnnGTZsWLX1ffzxx7z33nscOnQInU5HYmIigwcPZty4cef2xlyhJKgIIYSoUnZ2NqUlpQx9ZSiRiZEXfX/ph9P54rkvyM7OPuugsmrVKkaOHEmHDh2w2+2MHz+ePn36sGfPHnx9fQF48sknWbhwIT/88AMBAQGMGjWKgQMHsnbtWgAcDgf9+/cnMjKSdevWkZaWxpAhQ9Dr9UyZMgUovytz//79eeSRR/j6669ZtmwZDzzwAFFRUVx33XWV1vbZZ5/xxBNP8M4773D11VdjsVjYsWMHu3btqoF368ogQUUIIcQZRSZGEtfEM0/HLF68uMLz2bNnEx4ezpYtW+jevTsFBQV8+umnfPPNN1x77bUAfP755zRp0oQNGzbQqVMnfv/9d/bs2cPSpUuJiIigdevWvPTSS4wdO5ZJkybh5eXFrFmzSExMZPr06QA0adKENWvWMGPGjCqDyvz58xk8eDDDhw93LWvWrFmFNps2bWL8+PFs27YNm81G69atmTFjBm3btnW1URSFWbNmsWDBApYvX058fDyfffYZYWFhPPDAA2zatIlWrVrx5ZdfUq9ePQAmTZrETz/9xIgRI3j55ZfJyclhwIABfPzxxwQEBFRar9Pp5NVXX+Wjjz4iPT2dhg0bMmHCBAYNGgRAXl4eo0aN4vfff6eoqIiYmBjGjx/Pfffddy4/snMig2mFEEJcVgoKCgAIDg4GYMuWLdhsNnr16uVq07hxY+Li4li/fj0A69evp0WLFhVOBV133XWYzWZ2797tanPqNk62ObmNykRGRrJhwwaOHj1aZZvCwkKGDh3KmjVr2LBhAw0aNKBfv34UFhZWaPfSSy8xZMgQkpOTady4MXfddRcPP/ww48aNY/PmzaiqyqhRoyqsc/DgQb7//nsWLFjA4sWL2bZtG48++miVtUydOpU5c+Ywa9Ysdu/ezZNPPsk999zDqlWrAJgwYQJ79uxh0aJF7N27l5kzZxIaGlrl9mqC9KgIIYS4bDidTp544gm6dOlC8+bNAUhPT8fLy4vAwMAKbSMiIkhPT3e1OTWknHz95GvVtTGbzZSWlmI0Gk+r54UXXmDgwIEkJCTQsGFDkpKS6NevH4MGDUKjKe8rONnLc9JHH31EYGAgq1atYsCAAa7l9913H4MHDwZg7NixJCUlMWHCBFdvzujRo0/r2SgrK2POnDnUqVMHgHfffZf+/fszffp0IiMrnsqzWCxMmTKFpUuXkpSUBEDdunVZs2YNH374IVdffTUpKSm0adOG9u3bA5CQkHDaMdc06VERQghx2Rg5ciS7du1i7ty57i4FgKioKNavX8/OnTsZPXo0drudoUOHcv311+N0OgHIyMjgwQcfpEGDBgQEBGAymSgqKiIlJaXCtlq2bOn6/mRgatGiRYVlZWVlmM1m17K4uDhXSAFISkrC6XSyf//+02o9ePAgJSUl9O7dGz8/P9djzpw5HDp0CIARI0Ywd+5cWrduzbPPPsu6detq4F2qnvSoCCGEuCyMGjWKX375hdWrVxMTE+NaHhkZidVqJT8/v0KvSkZGhqtXITIyko0bN1bYXkZGhuu1k19PLju1jclkqrQ35VTNmzenefPmPProozzyyCN069aNVatWcc011zB06FBycnJ4++23iY+Px2AwkJSUhNVqrbANvV7v+l5RlCqXnQxA56qoqAiAhQsXVgg3AAaDAYC+ffty9OhRfv31V5YsWULPnj0ZOXIkb7zxxnnt82xIj4oQQoha7eTYjHnz5rF8+XISExMrvN6uXTv0ej3Lli1zLdu/fz8pKSmuUxxJSUns3LmTzMxMV5slS5ZgMplo2rSpq82p2zjZ5uQ2ztbJ7RUXFwOwdu1aHn/8cfr160ezZs0wGAw1Nm9NSkoKqamprucbNmxAo9FUuHT71LoMBgMpKSnUr1+/wiM2NtbVLiwsjKFDh/LVV1/x1ltv8dFHH9VIrVWRHhUhhBBnlH443WP3M3LkSL755ht+/vln/P39XWNKAgICMBqNBAQEMHz4cMaMGUNwcDAmk4nHHnuMpKQkOnXqBECfPn1o2rQp9957L6+99hrp6ek8//zzjBw50tWb8Mgjj/Dee+/x7LPPcv/997N8+XK+//57Fi5cWGVtI0aMIDo6mmuvvZaYmBjS0tJ4+eWXCQsLcwWcBg0a8OWXX9K+fXvMZjPPPPPMGXtozpa3tzdDhw7ljTfewGw28/jjjzN48ODTxqcA+Pv78/TTT/Pkk0/idDrp2rUrBQUFrF27FpPJxNChQ5k4cSLt2rWjWbNmWCwWfvnlF5o0aVIjtVZFgooQQogqhYaGYvQx8sVzX1yyfRp9jOd0JcnMmTMB6NGjR4Xln3/+uWsythkzZqDRaLj11lsrTPh2klar5ZdffmHEiBEkJSXh6+vL0KFDefHFF11tEhMTWbhwIU8++SRvv/02MTExfPLJJ1VemgzQq1cvPvvsM2bOnElOTg6hoaGunpmQkBAAPv30Ux566CHatm1LbGwsU6ZM4emnnz7r469O/fr1GThwIP369SM3N5cBAwZUOO7/9dJLLxEWFsbUqVP5+++/CQwMpG3btowfPx4ALy8vxo0bx5EjRzAajXTr1u2ijwdSVFVVL+oeLiKz2UxAQAAFBQWYTCZ3lyPcLDoilNS5T1Xf5o7ppGZcuqnAhagtysrKOHz4MImJiXh7e1d4TabQr51OzqOSnJzslv1X9zt1Lp/f0qMihBCiWnFxcRIchNvIYFohhBBCeCwJKkIIIcRlaNKkSW477VOTJKgIIYQQwmNJUBFCCCGEx5KgIoQQQgiPJUFFCCGEEB5LgooQQgghPJYEFSGEEEJ4LJnwTQghRLVkZtrL17Bhw8jPz+enn35ydylVkqAihBCiSikpKTRp3JiS0tJLtk8fo5G9+/adV1iZNm0a48aNY/To0bz11luu5WVlZTz11FPMnTu3wr1+IiIiXG1SUlIYMWIEK1aswM/Pj6FDhzJ16lR0un8/KleuXMmYMWPYvXs3sbGxPP/88677CVXl448/5r333uPQoUPodDoSExMZPHgw48aNO+fjuxJJUBFCCFGl7OxsSkpL+eDeITSs5I67Ne1AejqPfjmH7Ozscw4qmzZt4sMPP6Rly5anvfbkk0+ycOFCfvjhBwICAhg1ahQDBw5k7dq1ADgcDvr3709kZCTr1q0jLS2NIUOGoNfrmTJlCgCHDx+mf//+PPLII3z99dcsW7aMBx54gKioqCpvTPjZZ5/xxBNP8M4773D11VdjsVjYsWMHu3btOsd35solQUUIIcQZNYyMpGVsrLvLqFJRURF33303H3/8MS+//HKF1woKCvj000/55ptvuPbaa4HyOys3adKEDRs20KlTJ37//Xf27NnD0qVLiYiIoHXr1rz00kuMHTuWSZMm4eXlxaxZs0hMTGT69OkANGnShDVr1jBjxowqg8r8+fMZPHgww4cPdy1r1qxZhTabNm1i/PjxbNu2DZvNRuvWrZkxYwZt27Z1tVEUhVmzZrFgwQKWL19OfHw8n332GWFhYTzwwANs2rSJVq1a8eWXX1KvXj3g35sSjhgxgpdffpmcnBwGDBjAxx9/TEBAQKX1Op1OXn31VT766CPS09Np2LAhEyZMYNCgQQDk5eUxatQofv/9d4qKioiJiWH8+PHcd9995/LjOicymFYIIUStN3LkSPr370+vXr1Oe23Lli3YbLYKrzVu3Ji4uDjWr18PwPr162nRokWFU0HXXXcdZrOZ3bt3u9r87/avu+461zYqExkZyYYNGzh69GiVbQoLCxk6dChr1qxhw4YNNGjQgH79+lFYWFih3UsvvcSQIUNITk6mcePG3HXXXTz88MOMGzeOzZs3o6oqo0aNqrDOwYMH+f7771mwYAGLFy9m27ZtPProo1XWMnXqVObMmcOsWbPYvXs3Tz75JPfccw+rVq0CYMKECezZs4dFixaxd+9eZs6cSWhoaJXbqwlu71E5ceIEY8eOZdGiRZSUlFC/fn0+//xz2rdv7+7ShBBC1AJz585l69atbNq0qdLX09PT8fLyIjAwsMLyiIgI0tPTXW1ODSknXz/5WnVtzGYzpaWlGI3G0/b9wgsvMHDgQBISEmjYsCFJSUn069ePQYMGodGU9xWc7OU56aOPPiIwMJBVq1YxYMAA1/L77ruPwYMHAzB27FiSkpKYMGGCqzdn9OjRp/VslJWVMWfOHOrUqQPAu+++S//+/Zk+fTqR/3Mqz2KxMGXKFJYuXUpSUhIAdevWZc2aNXz44YdcffXVpKSk0KZNG9dndEJCwmnHXNPc2qOSl5dHly5d0Ov1LFq0iD179jB9+nSCgoLcWZYQQoha4tixY4wePZqvv/4ab29vd5dzmqioKNavX8/OnTsZPXo0drudoUOHcv311+N0OgHIyMjgwQcfpEGDBgQEBGAymSgqKiIlJaXCtk4de3MyMLVo0aLCsrKyMsxms2tZXFycK6QAJCUl4XQ62b9//2m1Hjx4kJKSEnr37o2fn5/rMWfOHA4dOgTAiBEjmDt3Lq1bt+bZZ59l3bp1NfAuVc+tPSqvvvoqsbGxfP75565liYmJbqxICCFEbbJlyxYyMzMrjOdwOBysXr2a9957D4vFQmRkJFarlfz8/Aq9KhkZGa5ehcjISDZu3Fhh2xkZGa7XTn49uezUNiaTqdLelFM1b96c5s2b8+ijj/LII4/QrVs3Vq1axTXXXMPQoUPJycnh7bffJj4+HoPBQFJSElartcI29Hq963tFUapcdjIAnauioiIAFi5cWCHcABgMBgD69u3L0aNH+fXXX1myZAk9e/Zk5MiRvPHGG+e1z7Ph1h6V+fPn0759e2677TbCw8Np06YNH3/8cZXtLRYLZrO5wkMIIcSVq2fPnuzcuZPk5GTXo3379tx9990kJyej1Wpp164der2eZcuWudbbv38/KSkprlMcSUlJ7Ny5k8zMTFebJUuWYDKZaNq0qavNqds42ebkNs7Wye0VFxcDsHbtWh5//HH69etHs2bNMBgMNTZvTUpKCqmpqa7nGzZsQKPR0KhRo0rrMhgMpKSkUL9+/QqP2FMGUoeFhTF06FC++uor3nrrLT766KMaqbUqbu1R+fvvv5k5cyZjxoxh/PjxbNq0iccffxwvLy+GDh16WvupU6cyefJkN1QqhBBXtgP/jNPwtP34+/vTvHnzCst8fX0JCQlxLQ8ICGD48OGMGTOG4OBgTCYTjz32GElJSXTq1AmAPn360LRpU+69915ee+010tPTef755xk5cqSrN+GRRx7hvffe49lnn+X+++9n+fLlfP/99yxcuLDK+kaMGEF0dDTXXnstMTExpKWl8fLLLxMWFuYKOA0aNODLL7+kffv2mM1mnnnmmTP20Jwtb29vhg4dyhtvvIHZbObxxx9n8ODBp41POflePv300zz55JM4nU66du1KQUEBa9euxWQyMXToUCZOnEi7du1o1qwZFouFX375hSZNmtRIrVVxa1BxOp20b9/edY16mzZt2LVrF7Nmzao0qIwbN44xY8a4npvN5gopTwghRM0KDQ3Fx2jk0S/nXLJ9+hiNNX4lyYwZM9BoNNx6660VJnw7SavV8ssvvzBixAiSkpLw9fVl6NChvPjii642iYmJLFy4kCeffJK3336bmJgYPvnkkyovTQbo1asXn332GTNnziQnJ4fQ0FBXz0xISAgAn376KQ899BBt27YlNjaWKVOm8PTTT9fIcdevX5+BAwfSr18/cnNzGTBgQIXj/l8vvfQSYWFhTJ06lb///pvAwEDatm3L+PHjAfDy8mLcuHEcOXIEo9FIt27dmDt3bo3UWhVFVVX1ou6hGvHx8fTu3ZtPPvnEtWzmzJm8/PLLnDhx4ozrm81mAgICKCgowGQyXcxSRS0QHRFK6tynqm9zx3RSMy7dVOBC1BZlZWUcPnyYxMTE0walyhT6tdPJeVSSk5Pdsv/qfqfO5fPbrT0qXbp0OW3k8YEDB4iPj3dTRUIIIf5XXFycBAfhNm4dTPvkk0+yYcMGpkyZwsGDB/nmm2/46KOPGDlypDvLEkIIIYSHcGtQ6dChA/PmzePbb7+lefPmvPTSS7z11lvcfffd7ixLCCGEqPUmTZrkttM+NcntM9MOGDCgwsx7QgghhBAnyb1+hBBCuLjx+gpxmamp3yUJKkIIIVwznJaUlLi5EnG5OPm7dOrsuefD7ad+hBBCuJ9WqyUwMNA1M6uPj49rSnYhzoWqqpSUlJCZmUlgYCBarfaCtidBRQghBPDvPW1OnUZeiPMVGBhY6Qy450qCihBCCKD8pnZRUVGEh4djs9ncXY6oxfR6/QX3pJwkQUUIIUQFWq22xj5khLhQMphWCCGEEB5LgooQQgghPJYEFSGEEEJ4LAkqQgghhPBYElSEEEII4bEkqAghhBDCY0lQEUIIIYTHkqAihBBCCI8lQUUIIYQQHkuCihBCCCE8lgQVIYQQQngsCSpCCCGE8FgSVIQQQgjhsSSoCCGEEMJjSVARQgghhMeSoCKEEEIIjyVBRQghhBAeS4KKEEIIITyWBBUhhBBCeCwJKkIIIYTwWBJUhBBCCOGxJKgIIYQQwmNJUBFCCCGEx5KgIoQQQgiPJUFFCCGEEB5LgooQQgghPJYEFSGEEEJ4LLcGlUmTJqEoSoVH48aN3VmSEEIIITyIzt0FNGvWjKVLl7qe63RuL0kIIYQQHsLtqUCn0xEZGenuMoQQQgjhgdw+RuWvv/4iOjqaunXrcvfdd5OSklJlW4vFgtlsrvAQQgghxOXLrT0qHTt2ZPbs2TRq1Ii0tDQmT55Mt27d2LVrF/7+/qe1nzp1KpMnT3ZDpZeXti1bkpaWVm2bqKgotu7YcYkq8iwtWrch/QzvT2RUFDuTt12iioQQ4sqlqKqquruIk/Lz84mPj+fNN99k+PDhp71usViwWCyu52azmdjYWAoKCjCZTJey1FotKiyMHRNfqLZNyxcnk5aVdYkqqhnREaGkzn2q+jZ3TCc1I7vaNmERkYz/anW1babc052sjPRzrlEIIUT553dAQMBZfX67fYzKqQIDA2nYsCEHDx6s9HWDwYDBYLjEVQkhhBDCXdw+RuVURUVFHDp0iKioKHeXIoQQQggP4Nag8vTTT7Nq1SqOHDnCunXruOWWW9Bqtdx5553uLEsIIYQQHsKtp36OHz/OnXfeSU5ODmFhYXTt2pUNGzYQFhbmzrKEEEII4SHcGlTmzp3rzt0LIYQQwsN51BgVIYQQQohTSVARQgghhMeSoCKEEEIIjyVBRQghhBAeS4KKEEIIITyWBBUhhBBCeCwJKkIIIYTwWBJUhBBCCOGxJKgIIYQQwmNJUBFCCCGExzqvoFK3bl1ycnJOW56fn0/dunUvuCghhBBCCDjPoHLkyBEcDsdpyy0WCydOnLjgooQQQggh4BxvSjh//nzX97/99hsBAQGu5w6Hg2XLlpGQkFBjxQkhhBDiynZOQeXmm28GQFEUhg4dWuE1vV5PQkIC06dPr7HihBBCCHFlO6eg4nQ6AUhMTGTTpk2EhoZelKKEEEIIIeAcg8pJhw8fruk6hBBCCCFOc15BBWDZsmUsW7aMzMxMV0/LSZ999tkFFyaEEEIIcV5BZfLkybz44ou0b9+eqKgoFEWp6bqEEEIIIc4vqMyaNYvZs2dz77331nQ9QgghhBAu5zWPitVqpXPnzjVdixBCCCFEBecVVB544AG++eabmq5FCCGEEKKC8zr1U1ZWxkcffcTSpUtp2bIler2+wutvvvlmjRQnhBBCiCvbeQWVHTt20Lp1awB27dpV4TUZWCuEEEKImnJeQWXFihU1XYcQQgghxGnOa4yKEEIIIcSlcF49Ktdcc021p3iWL19+3gUJIYQQQpx0XkHl5PiUk2w2G8nJyezateu0mxUKIYQQQpyv8woqM2bMqHT5pEmTKCoquqCChBBCCCFOqtExKvfcc4/c50cIIYQQNaZGg8r69evx9vauyU0KIYQQ4gp2Xqd+Bg4cWOG5qqqkpaWxefNmJkyYUCOFCSGEEEKcV1AJCAio8Fyj0dCoUSNefPFF+vTpUyOFCSGEEEKcV1D5/PPPa7oOIYQQQojTXNAYlS1btvDVV1/x1VdfsW3btgsqZNq0aSiKwhNPPHFB2xFCCCHE5eO8elQyMzO54447WLlyJYGBgQDk5+dzzTXXMHfuXMLCws5pe5s2beLDDz+kZcuW51OOEEIIIS5T59Wj8thjj1FYWMju3bvJzc0lNzeXXbt2YTabefzxx89pW0VFRdx99918/PHHBAUFnU85QgghhLhMnVdQWbx4MR988AFNmjRxLWvatCnvv/8+ixYtOqdtjRw5kv79+9OrV68ztrVYLJjN5goPIYQQQly+zuvUj9PpRK/Xn7Zcr9fjdDrPejtz585l69atbNq06azaT506lcmTJ5/19kXlzOY8Wj7/n+rblJVdomqgbcuWpKWlVdumoLQIXz+/atsU5ubwxuuvVb+d/MJzrq8y+fkFhEVEVtsmMiqKnckXNnZLCCGudOcVVK699lpGjx7Nt99+S3R0NAAnTpzgySefpGfPnme1jWPHjjF69GiWLFly1pPEjRs3jjFjxriem81mYmNjz/0ArnBOp5OD43tX2yb0+Z8vUTWQlpbGjokvVNsm8qnRvLG4+pD6ZPuRjLi6R7VtJiypmeNyOp2M/2p1tW2m3NO9RvYlhBBXsvMKKu+99x433ngjCQkJrqBw7NgxmjdvzldffXVW29iyZQuZmZm0bdvWtczhcLB69Wree+89LBYLWq22wjoGgwGDwXA+JQshhBCiFjqvoBIbG8vWrVtZunQp+/btA6BJkyZnNc7kpJ49e7Jz584Ky+677z4aN27M2LFjTwspQgghhLjynFNQWb58OaNGjWLDhg2YTCZ69+5N797lpxAKCgpo1qwZs2bNolu3bmfclr+/P82bN6+wzNfXl5CQkNOWCyGEEOLKdE5X/bz11ls8+OCDmEym014LCAjg4Ycf5s0336yx4oQQQghxZTunoLJ9+3auv/76Kl/v06cPW7ZsOe9iVq5cyVtvvXXe6wshhBDi8nJOQSUjI6PSy5JP0ul0ZGVlXXBRQgghhBBwjkGlTp067Nq1q8rXd+zYQVRU1AUXJYQQQggB5xhU+vXrx4QJEyirZDKw0tJSXnjhBQYMGFBjxQkhhBDiynZOV/08//zz/PjjjzRs2JBRo0bRqFEjAPbt28f777+Pw+HgueeeuyiFCiGEEOLKc05BJSIignXr1jFixAjGjRuHqqoAKIrCddddx/vvv09ERMRFKVQIIYQQV55znvAtPj6eX3/9lby8PA4ePIiqqjRo0EDufCyEEEKIGndeM9MCBAUF0aFDh5qsRQghhBCignMaTCuEEEIIcSlJUBFCCCGEx5KgIoQQQgiPJUFFCCGEEB5LgooQQgghPJYEFSGEEEJ4LAkqQgghhPBYElSEEEII4bEkqAghhBDCY0lQEUIIIYTHkqAihBBCCI8lQUUIIYQQHkuCihBCCCE8lgQVIYQQQngsCSpCCCGE8FgSVIQQQgjhsSSoCCGEEMJjSVARQgghhMeSoCKEEEIIjyVBRQghhBAeS4KKEEIIITyWBBUhhBBCeCwJKkIIIYTwWBJUhBBCCOGxJKgIIYQQwmNJUBFCCCGEx3JrUJk5cyYtW7bEZDJhMplISkpi0aJF7ixJCCGEEB7ErUElJiaGadOmsWXLFjZv3sy1117LTTfdxO7du91ZlhBCCCE8hM6dO7/hhhsqPH/llVeYOXMmGzZsoFmzZm6qSgghhBCewq1B5VQOh4MffviB4uJikpKSKm1jsViwWCyu52az+VKVJ4QQQgg3cHtQ2blzJ0lJSZSVleHn58e8efNo2rRppW2nTp3K5MmTL3GFnqNty5akpaVV26a4qAhfP79q26iqWiP1tGxz5nqioqLYsW1HtW1yCvNpPGFctW1qqmanqhIWGVZ9PTkFvD799UtSjxBCiOq5Pag0atSI5ORkCgoK+L//+z+GDh3KqlWrKg0r48aNY8yYMa7nZrOZ2NjYS1muW6WlpbFj4gvVtol9YjSHztAm5snHaqyeyYurD44vXF99LQAOp5P/vHZrtW3GjpxzTrVV50w1j2w/mh5396i2zbdLZ9dYPUIIIarm9qDi5eVF/fr1AWjXrh2bNm3i7bff5sMPPzytrcFgwGAwXOoShRBCCOEmHjePitPprDAORQghhBBXLrf2qIwbN46+ffsSFxdHYWEh33zzDStXruS3335zZ1lCCCGE8BBuDSqZmZkMGTKEtLQ0AgICaNmyJb/99hu9e/d2Z1lCCCGE8BBuDSqffvqpO3cvhBBCCA/ncWNUhBBCCCFOkqAihBBCCI8lQUUIIYQQHkuCihBCCCE8lgQVIYQQQngsCSpCCCGE8FgSVIQQQgjhsSSoCCGEEMJjSVARQgghhMeSoCKEEEIIjyVBRQghhBAeS4KKEEIIITyWBBUhhBBCeCwJKkIIIYTwWBJUhBBCCOGxJKgIIYQQwmNJUBFCCCGEx5KgIoQQQgiPJUFFCCGEEB5LgooQQgghPJYEFSGEEEJ4LAkqQgghhPBYElSEEEII4bEkqAghhBDCY0lQEUIIIYTH0rm7AOEZnBYHlsxSbGYrjmIbk8IiSR07FjRatMFB6ELD8IqPw1C/PvqYGBSNZFwhhBAXnwSVK5i9xEbxQTMlh83Yci0VXrvRFEDBz/MrXU8TEIBP+/YM8vLG92g2xbEhoFEuRclCCCGuMBJUrkBxej3Zq1IpOWwG9d/lugAvvIIMaP30vLxqP69MexXV7sCRm4s9KxPL34ex/v03zoICipYt43EfPxj7LVZ/b3LaJpJ1VT3yWsTi9JJfKyGEEDVDPlGuIE6rlZING/gxri4lf5sBMEQY8a0fgDHWD63x31+HL35ez6zhw0/bhmq3U7Z7N8V/bmThq9No5+eDV2EZUav2ErVqL3ZvPTltEsjqWI/c1vGX7NiEEEJcniSoXCGsx49TtHw5zsJCdIqCd4wvgW3D8ArxPqftKDodxlatMLZqxdMTx/Pid08QsD+NsE2HCN14CO/cYiLW/0XE+r9w6LVM8fWnYMEC/K65Fq2fb4VtldnLOJR/CNNVJvZoM7BgL98HCga0+KhemFQDAaqxxt4HIYQQtYsElcucqqqUbt1KyYYNAGj8/Rm+bze/3denZrav05LfLIb8ZjH8NaQ7pkMZhG48RNjGQ/hkFNBVbyD1mWdRDAa8r+7KX9fWJzmilD+ztvBX/l84VScxD8awmeNV7kOnakh4JoHkzGTqBtTFZDDVSO1CCCE8nwSVy5hqt1O0YgWWAwcAMDRpgm/Xrmwa+9TF2aFGwdwgEnODSP6+qzO+KTnsGf0pfTo2YFFEGusbr8RcsAoK/l0lyBDI8R0naFY/BgM6FBWcCliwU6xYKVBKsSlOfBv5sjF1HfuPryHaO4xGwY2IDqyP3eALigzkFUKIy5UElcuUarNhXrQI27FjoCj4duuGsUWLS7d/YE9QMfOGGPi1XgYnp+wxlWlot99Oy8MqTY6phDiLWZBhIeKeILLiA3DqytsF2IuJtWYRaS0lyF6IT1EeoV4nt3ICSAbAoSgUeweQZwqn+BoDWkcODk0AKPKrLYQQlwO3/m8+depUfvzxR/bt24fRaKRz5868+uqrNGrUyJ1l1XqqzYZ54UJsJ06AToepXz+8YmMvzb5VlZTCFDalbyKrNAtDPQNaRcu1cdcysMFAOkZ2xL5nP+bFiyhctBhbaioD/E3w83504U50zZ2EhBZh0pRW3LDXP9unPJyoqooe0KoqptJ8TKX5PNPFANbtqCg4NEE4tBHYtWESWoQQohZz6//gq1atYuTIkXTo0AG73c748ePp06cPe/bswdfX98wbEKfRK4orpCh6PaYbbkAfFXVJ9p1TmsO61HUcLyofb6LT6Mhfns/699YT6Rv5b40tmmNs0Zzwp5+mbPN6Fo+9ketag4/J5mqjOiG/1Ifj3qGkBIUy4c0/eOz/XqHM4IuqaLA77ezP2sGRtM2EWoqpb7XR8FgRt/pq0KkWdM5cdM5cvGwHsGujseliUDUyKFcIIWobtwaVxYsXV3g+e/ZswsPD2bJlC927d3dTVbWXqqpMjoj4N6TceCP6yMgzr3iBbA4bG9M3sjN7JyoqGkVDi9AWtAlvw7RnplUIKS4luSh/zsL454fc0r08oNhVDVmFJsoO6LAc0aHaNeiwE+OdzY2OMLRpZaiJ/kB5CGoW0Zb6oc3YmLaR73N2QTN/Ztk1TMsLoE2ZGZ0jHY1ait5xDJ3jOHZtJLEmGc8ihBC1iUf1iRcUlI+yDA4OrvR1i8WCxfLvDKpms/mS1FVblKxfT19/E2g0+Pfte0lCSoo5hVXHV1FkKwKgbkBdkqKSqr4yx1oM6z+AtW+DtRCAAzlOjiW2YrdPHJY6XtBQJeREIVEHcog6kIuh1M7tAYEwbi75jaI4dkNbstslgqJg0BroFtONRsGN+GbV12RHe/NgaB5DioN4qDAeL2ceOvsxdM5c9I40dj/qx649v7Gl4TXYdV4X/f0RQghxYTwmqDidTp544gm6dOlC8+bNK20zdepUJk+efIkruzRatmlJWlpatW3MednEjnm80tf6+5l4OaL8FI/ftdde8JgUVVUx6KrufVC8FEIHRbLw8EIA/L386V6nO3GmuArtivPziI4IBVQGNlCZkOQg4p+zeruz4e0tWr7abMbgsxnYfNp+tEBrL2+u9/Kjj58fgfvTCNy/kP0WCzPz8lhaXPxvYz08N70D83zMfOGXxwG9hRfzI/DXhmBzFOBlP4RBl0+7AysJ3bSckYus/HrQWenxrV+37lzerkrl5+f/c+xVi4qKYkvyzmrbtGjdhvQz/G4UFxcT6Gu4JPuKjIpiZ/K2atucze9zUVEZfmc4xXs2+xJCXN48JqiMHDmSXbt2sWbNmirbjBs3jjFjxriem81mYi/RINGLLS0tjcmLqw9hT7YfybKpHU5brs12YFpQDA74MDeb52pgMLIKLJ92VaWvpeicfOZvIV1XPv9+i9AWdIzqiF6jP62t0+Ek9dOh8NdiKPhnrhTvQEjsTrOwJnx0q8IXvcfzy4RW1dbT7T9bMY+4l2a7DtBk9wEaAW9FRpIREcrGTm3Iigxj+jMvMtYcThurkVcCMllvKGF4yHFm5EVThwDKNG24e85iPhscQnxAEb/c4c2fNOInOmNR/u1dmfHSDNrHVz+r7pFqX/3n2J1OUudWfyl49B3Tz7id9LQ0xn+1uto2T/drTuqC8ZdkX1PuOfNp2bP5fX6sw5OMn7flgvclhLi8eURQGTVqFL/88gurV68mJiamynYGgwGDofq/Gq80SpkTv6UlKA6wxuqYdTCb5y7Svpyo/G60M9/XhlMBe76Nm9rcfFoviqs21cnojnrY8jmoDtDoIb4zxHQAzbn/6pX6+rC5Y2t2tGpCs537ab5jLxEZ2dzw8xL2N67Hx//c0blPmT/xdi/GBqWRorPxQPBxZuRF0djuzfwDDl5lMNezmavZQUf2U59UvlR7kaKEX9D7I4QQouZp3LlzVVUZNWoU8+bNY/ny5SQmJrqznNpHVfH9owxtkYrDpFDcw3jqPQZrVLGiMtNk4Se/8pDSxqLlyMSDVYYUU1E2N//xEa/28i4PKUF1ocMDEJd0XiHlVFZvA9s6tOT/7riBA43qAtBo3yEWxMVRcrR83Esju4FPcmJoaPMiT+tgRPAJtniVAGBT9CxQkviAG8jBnxAKGcXPdFV3gXqx3kEhhBDnw61BZeTIkXz11Vd88803+Pv7k56eTnp6OqWlpWdeWWDYZ8PrqB1VA0XX+qAaLs4VLce1TqYGlbHT4ESnwr2FXjxk9sJZ7Di9sarS7PAGbl/xDlG5Rym0qNCwL7S4DbwDarSuUl8f1vToxMIbe5EXaCJUpyN7+Qly16XjtDsJdeqYmRtDB4uRUo3KU4Fp+DTyca3/txLNdG5lO4nocDKQtdzOKvRu/VchhBDiVG79L3nmzJkUFBTQo0cPoqKiXI/vvvvOnWXVCpo8Bz4bygAo7WDAEaq9KPvZYLDzalAZ2VqVEIfCs/nedCnToXB6KNLbyuiz6Vuu3v4zeoeN46F1af9xMUS1uqjT3GdEhfPzoL58lpcHQNH+fDJ+OYq90IqvquGNvCg6WXwo06jEj44h1zvHtW6ZYuALevMTSThR6Mh+Ft1pwFeVsCyEEJ7A7ad+KnsMGzbMnWV5PqeK38rS8nEpMVrKml+Ey2w18L2vldkmKzYFmlk1jM/zJs5e+a9McEE6g1a9T/3UnTgUDWua92d+l+EcLbg0p1KcWi3Tc3II6xOLxluLLc9C+oKjlKUWY0DDq3mRdLL4oDFo2Ba+uUJYQVFYrbTkU66nDD09E7W8bP2CcGfeJaldCCFE1aSTuxby3m5Fl+PEaYDi7sYa760oVpzUeSye5T52APoX6xhZYMBXrXw/DVO2cuvqDwgqyqbIGMBP3R5iR/2uoFz6Xy9jHV8ib0zAK9Qbp8VB5u/HKD5U4AorhTuLcGqcp4cVYK8Sx3vcxDGzkxg1h6nW2TR0Vn1XZyGEEBefBJVaRpvjwLitfNK7kiQjqk/N/ghTtTYeDD6Ob0t/9Co8WODFDSVeaCo51aNVoOv2+fTa+gN6h42U8AZ83+MxMoKrv7T3YtP56gnvG4dPXROokLM6DfOuHAxoOPbeCUJKQnFqnCSHb8HsVVBh3VQlhI6fWfhbicRECROtX9PCcdhNRyKEEEKCSi2iBXxXl6I4wRqvw1qvZq8u36EvZXjwcf7WW7Hn23g630A7a+X70KoO5t/pQ8vD6wHY1KgnC5OGUWbwjHs0aXQaQrpH4d80CID8TVkUJGej2lVaZbUlqDQYh8bBtojNlOhKKqybVqQy0etetmnqYsDOf2zf0cZx0B2HIYQQVzwJKrXI4ICg8lM+XlDcxbtGT/ksNxQxKjiVPK2DhjYvUl4+RLy98gG63k4rLUuO0quuDptWz6Kr7mFTk16objjVUx1FUQi8KpyAdmEAFGzLZkRQEFpVS6vMtvhZ/bFqrWyL2IRVY6mwrkXx4jX9bWzUNMQLB8/afuAqxz53HIYQQlzRPOuTRVTJK7eIkSHl07GXdvCu0VM+/+eTz3OB6VgVla5lPszKjcGeZ6+0bYC9mJYlR/BRrRwzO/mx2yMcjm5WY7XUNEVRCGgZQuA/YWVUSAit9hxEr+ppk9Eeb7uREn0J2yK24FAqHrNd0fGmfiDrNE3Q4WSM7Uc6O3a74zCEEOKKJUGllqg/5w/8NFrs4VosjU+fqv58qKjM9MvhDVM2qgI3l5iYlh+Fj1r5r0WYrYCmZcfQ46RQ402XT4vJCYyukVouNtMpYeWq7XtpdCgFb4c3bTPao3foMRsK2BGWjJOK9/5xKFre1t/MKk0LtKiMtv1MV8cudxyCEEJckSSo1ALByUeJ2HAQu6rW2CkfOyqvmDL5wq/8EtwHC4MZaw5DV8mgWYBoaw4NLWlogCydPzuNcWQU165ZXE0tQ/j4n7lWum7aTvzxdHxtfrTObIfGqSHbJ4v9wXtPW8+paHhffwNLta3RoDLKNp+bGl6ceWuEEEJUJEHFw2msdhp+thKAbwvycIRc+Aek4qXwbGAav/gUolHhPwVhDC8OrnQSN1SVBEsGidYsAE7ogzhgiPa48Shn662cHPbXjUWjwrXrthCam0+gJYgW2a1BheOmFIKvCTxtPVVR+EjXz9Wz8t1AL+LyNlzy+oUQ4kpTOz9triDx8zZjzDRTFuzLzJzsC95egVJKzDOJrPMuwaAqTMuP4ubSyqe2V1QnDS1p1LGV90Ic9grjiCHios4yeyn80aElx6LC0Dmc9Fm9CZ+SUsJLIqif1xCAqDsjKPZKPW09VVH4QD+ADZrGGHQKN+57mjoF2y51+UIIcUWRoOLBjKl5xM3fAsBfw7pTojrPsEb1MpUiHveZj7GuDyanhndzo+luqfxyYn8vaFp2nDC7GSdwwBBFqlfIBe3fU6gaDcs6tyPP5IdvaRl9Vm9Ca3eQYK5LVFEdFK3CiaDVWLQFp63rVDS8rb+ZXw860Dst3LT3SSIKZYCtEEJcLBJUPFiDOX+gcTjJaRNPdod6F7St45oCHvf9mWPafGw5Vj7MiaGlzVhpW0W1sHSIL4GOEhxo2OsdS5a+Zm8o6G42Lz2/Xd2RUoMXYXkFdN66CwWFptnNKP6rBKfGxvHgFTgUy2nr2hUtt/6fhWOmdhgcxdyyZzRBJUcu/UEIIcQVQIKKhwredoSQ5KM4tRr+GtL9gk63/K3JYbTPz2RqiohxBHBs2mESHZXfH0hxluJt2ULrSC1WRctOYxz5Os+YxK2mFfr5sLxzW1Sg8aEUGvx9DA1aUt4/gc7ui01XyImg1aic3pNVZoefm0wnza8ZRnsBt+x5HF9L1qU/CCGEuMxJUPFAit1B/S/XAHD8+laURgWe97b2ajN40ncBeZpS6jqCeavkRuy5tsr36yzG27IVjVrG4TwnO43xFGu9z3vftUFqZBhbWjQCoOvmHQTnmXEUOojJuwaNU0eJIZ1M05ZK17XpfPmp6VvkescRYEnjlj2j8bIXXcryhRDisidBxQPV+X0nvql5WE1Gjtza4by3s017gqd9FlKoWGhqD2dG8Y0Eqz6VtlWcRRgtW9Fgwan4cO2cYso0F+GuzB5oW7MGHIssH1zba+1mfBUFb3sQUfldAMjz3Ue+sfIp9Mv0gcxr9g7F+hDCSv7ixr1Po3VaL2X5QghxWZOg4mH05lIS/rsRgL9v74TDx3Be21mnO8J/fBZRqthoa6/D6yUD8KfybWmcZoyWrSjYcCh+lBraklpYu+ZIuSCKworObSjy8SagsJiXIyJAVfG3xBFa2BKAjIA/KdVXfmrH7F2HeU3fxqL1Jda8hesPvAAXOPBZCCFEOQkqHibxhz/RF1sojA8l7Zqm57WNFbqDvGBcgk1x0NkWz5SS6zFS+Wy2Gkce3pZtKNhxaEyUGdqAcmX0pJzKYjCwrEt7HBqF6/z86LivfOK3kKKW+JXFoipOTgStwqYpqXT9LL9GLGj8Og5FR8OcpfQ4/OalLF8IIS5bElQ8iG9KNtFLy6dnPzi0O2jO/cfjf1UArxiX41Cc9LLWZ1Jpb7yo4g7Ijhy8rdtRcODQBFHm1RqUmpmevzbKDA3iz9bl4bDPps2E5eWhoBCV3wWDLRC7tpQTQatw4qh0/WOBHVjcYDIAbdK+45mkmr27tRBCXIkkqHiQ+nP+QFFVMjvWJ79pnXNef6PBTuSDMTgVlb7WRvyn7Fp0VD6TrdaRicG6AwUndk0IZV4tQZEP1t0NE1ldXIze6WDQ6tVoHQ60qp46eT3QOL0o88omPaDqGWkPhPVhZeKTAEy7Vg/pOy9V6UIIcVmSoOIhuuq8CN51HIdey6G7u5zz+hsNdj73t6JoFPpaG/FU2dVoqrhvz53N9Risu1FQsWvDsXi1AEXuXQOAovB8ZibFBgNReblcu20rAF4Of+rkdQdVwezzNyG9gqrcxLbou9hc597yJwcWQe7hS1G5EEJcliSoeACn1cqjxvK5So4NaENZuOmc1j8ZUlQFClbnVhtSvMs28tlN3iio2LSRWPTNoJbet+diyXY4+LlzeVjssmsXCenpAPhaowg3twMgcnA4Jyzbq9zGH/Gj+HqnvXxQ7Z55UJh+8QsXQojLkHxCeYC8OXOI0WqxBPmSclO7c1r3z1NCSpdSLRlzUqsMKcbSNfgX/4xGUbBpY7Dqm9T6+/ZcLPvi49nSoAEaYOAfqzFYyy85DippjKmkLopWYWXumxTY0rDanZjLbGQWlpGSW8JfmYUcyCxm5F+d+NFwMz9YOjFvawq/HdeyJsePrfk+HCgykGeVXiwhhDgTGZTgZvasLLJnzgLg0J1JOLzP/oqbPw12Zp8SUu4u8uLzyq4qVlV8SpfjW7ocgNfWWhjZs4GElDNYdFVHEtPTCS4s5NqNW/mu47WUOr2xpd9CoS2dYn0YX+zNQ1Urn+QtsN8Yxpy8XZAdqOSWQPp7P6bba8uJ8PcmJshIvTA/6oX7US/Mj/gQH7z1EmaEEFc2CSpulvnWWziLi9ljt5HRtfFZr1dZSKm0J0VV8S1ZjE9Z+Uy3xcbeTFjxX0b2kpBSJY2WIruRQkcoz1wzAkeplgzfECg+tVEYzlNuA6RRwKjX4q3X4qXToFUUDm7fQO+2ddGqNuzmDIodOko0fhTrgylyaMm36VB0XhzLLeVYbimbj+ZVLEOBBuH+eHW5j+3H84kweRPq54XuPK4GE0KI2kqCihuV7txFwY/zAHintJjbNWcXHk4NKV1LtdxVZUhx4lf8C0bLnwAU+fSn1NgZ+G9NHcJlQVWhQPUlwxFMljOQ2FFfsbfonxs2aoB/bnUUVpKH0w/0ege7f/k/bn3yYTaXvgW6PNqYbqGtaXCF7T49YRJfPlR+uTLFWZD8FdgtENoQmt6MRdWS+NBs5i9bR7q5jGO5pRzKKuJQVhEHM4soLLOzP6MQfcOurNxfPtmcVqMQFeBNTKCRmCAfIgIMElyEEJc1CSpuoqoqGVOmgKoScNON7Pnis7Na7+xDigP/4nl4W7aholDkexNl3uc/Hf/lxqkqZDsDSHWEkuYMoVT9955GGi/QKXb8dSX4aUvwU4p5dOlP1M9KY09cHHOvuZY/ty+mVcR0jCWDWJP/PslF3xGijyfe2LHyHfqGQbNbYcd3kH0ADi7FUL83FGXTPiH4tOaqqpJZaGHH8QKGjplIg553kGEuo8zm5HheKcfzSuFwLlqNQkyQEb/WfUkp8SLOR6bvF0JcXuRPMTcxL/yV0m3bUHx8CBsz5qzW8e8U4Aop3Up11YeUoh/+CSkaCv0GSUgBnCoY6jRlq7UBC8uSWGNtxd+OOpSq3mhxEKXJpqX+IKmzR9PatJ/6vseI9M7Bz1DGb0lXYddoaJqSQqu/D7m22dDnWpr69gNgVf475NlSqi4gMA4aDyj/PnUrHPuzyqaKohBh8qZ30whs237m5tZ1eKhbXe7tFM81jcJoEO6HUa/F4VQ5mlNCUM+H6P5HY3quaciU/VEk5xtRr6C7IAghLl/So+IGzpISMl9/HYDQhx5CHxFxxnX25+4ncniMK6TcWaSvNKQYtGAq/AaDbR8qWsx+t2M1NKvxY6hNjjpM/GRtyC/WekTc6c+RfyaWNWAlUptDtDaHME0eOqX8/jy27KOnjTPOCA5mZevW9Nq6lX4b/uQ97b+DXK8yDSXPdow0606W5k7jxrBXMWj8Ky8mvAlYC+HQcji8klsanP1gWUVRCPb1ItjXi5YxgaiqSk6xlaM5JSxdsRLfuGYcKvbmULE3Hx0JI8ZopX9EPgMiC2huKj2n90wIITyFBBU3yPnkE+wZGehjYgi+b9gZ2+/P3c/yY8tRNEq1IUWjOpl3u88/IUWH2f8urF6NLsIReD6bqrDSFs93liast/87y6/TUkKiTwGx2kzCNPnndOHTmuYtaJySQkx2Ni9FRFCkqqAoaBQd1wQ/xfysZyl0ZLAibwZ9gp+rekMxV0FZIZzYxIxrHPD3Sqjb45yPUVEUQv0MhPoZ+Pb7Cexf+BJ/5PjzW0YAy7L8OV7qxYdHwvnwSDhxRgvadrexN81Mk6hzm6dHCCHcSU79XGLW4yfI+bR8PEr42GfRGKq/O/LJkAKQvzK3ypCiVR00LT1Gz7o6VLwoMA29IkNKttPI+6Vt6VNwB08U92K9vQ4KKt10x3jDdzknZg6jndcBwrXnFlIAnBoNP3bthk2rpauPD7Frf3O95q3xp1fwf9Ap3qRatvNnwRnGHNW7FsIa46UFvrsX0ned+8H+jwC9kwGRBbzbKoUt1+xhZqsj9I/Mx6h1klJqQNtyAH3f/oP+7/zB7LWHySuW8SxCCM8nQeUSy3zjDVSLBZ9OnfDv1avatqeGlKYhTcn8qvLJ3HSqg2alxwhwlpJfppJvug+bvu5Fqd9Tme0Ggns/Su+C2/mgrC2Zqi8hSikPeifzm+k7Zvn/Rl+vv1HtF/bhnB0YyJJ25ZPyNZr3OcasNNdrwfp4ugc+BijsLVlMcO+qp9lHUaDxANanKmAxw9e3QcHxC6rtVEatSt9IM++3SmFLj9283+ooziOb8NJq2J1qZtKCPXScsoxHv97Cin2Z2B3OGtu3EELUJAkql1Dxxo0ULl4MGg0R48ahVPMn/f+GlO51ukMlgyP1TjvNS1Pwd5ZhQ0ufL4ux6+Mu1iF4nGyrL2vy6vJbTjP8WvXBio6W2kxe913O0oBvecK4mTrayidkO19/NmnKxtJSdNYyWnz9Djj/vZtygrETHUxDAIi4LYwlJdUMrtXouH+RFsIaQ2EqfHUrlOTWaK0APjqV/pEF2Fe8x5/jezLphqY0izZhdTj5dWc6983eRJdXl/PGb/tRfE+/AkkIIdxJgsolojocZEyZCkDQHbfj3ahhlW335e5zhZRmIc3oXqd7paHGy2mjeWkKvk4LVkXLTp84tmdcGX8ZZ1t9WZVXnxV5jUizBgIqJX9t4Ev/BXzjP59+Xn/jpVyc90JVFJ7LyMBu8Cb40B4SVv5S4fXmvjfQxOd6FI3CuJz1bLdkV7mtAqsCd/8f+EdB1j74ehBYCi9K3QBBvl4M65LIwse78evj3bivSwLBvl5kmC28t+IgxkGvsmB7KkdzilHlsiEhhAeQoHKJ5P/wf1j27UMTEEDoY49V2W5f7j5WHFsBlIeUbnW6VRpSDE4rLUpT8FGtlCk6dhrjKdVUP97lcnBqQMm0mlBQSTRmc33IHrJ/nkZbXcYluTPACbudfbcMB6DBgi/xTT/mek1RFDoG3E/h9iIsqoPHslZxzFZN+AiMhXt/AmMwnNgC394JtrKLfATQNNrECzc0Y8O4nsy8uy1d6oegaDT8nV3MT8mpzFl/lK0peZTZHGfemBBCXCRuDSqrV6/mhhtuIDo6GkVR+Omnn9xZzkVjz8sj6623AAgbNQpdUOVjF842pPg4ymhZehRv1UapomeXMZ4yzdnfI6g2MtRpclpAqWvMom/obtqbUvDXWc68kRp2vHNvspq0RWu30fLLt1Ac/36gaxQtxz9KpYk+iDynhRFZK8lxVBM+whvDPf8FL3848gf8333gsF2CowAvnYa+LaL4+oFOlPz4HK1jAvHSasgvtfHHX9l8uuYwS/ZkkGG++OFJCCH+l1uDSnFxMa1ateL99993ZxkXXeYbb+DIz8fQoAFBd9xeaRtjO+NZhRR/RwktSlPwUh0UawzsNMZj0egvav3udNQewasFdxFx59TTAko70zF8tW68ckVR2HXXKGxGXwJSDpK4pOKtCVSLyvvhPYjW+nLUXsgjmSswO6upt05buGsu6Lxh/6/w06PgvLSn8tSCdK5uFMbwrolc2zicUD8v7E6VPWlm5m46xnebjqFN7CCDb4UQl4xb51Hp27cvffv2dWcJF13J5s0U/PdHACInT0bRnx4q5v01D9Og8rktqgsp19fT0az0GFpUzBoje4wxOJTL8+66WY4A/q+kB2stLVFRUB126vnl0dg3w73h5H9YAkPYM/hhWn3xJvUXzSWreXsKY/694ipMa+Sj8GsYmrGUfbY8RmWtYlbYNfhoqvinl9AVBs+BuXfBzu/B20Slo6gvMi+dhhZ1AmgebSKtoIwdJwo4mFFEurkM7x6P0P21FQztnMAdV8URYLx8g7IQwv1q1RgVi8WC2Wyu8PBkqtVK2qRJAATedhs+bduc1mbeX/N4Yd0LKBql2pDS4Fgy/zfYiBaVXK0vu42xl2VIKXJ683Vxb57OG8kaSytUFDp57SLts1Hu70GpQlq77qS3TkLjdJSfArJVPGUTrzfxYfg1+Ct6tlmyeDJ7NVa1mnEfDa+DWz4EFNj0CeOvsl/cA6iGoihEBxq5vlkk93VJoGNiMGqpmdSCMqYu2kfS1GW88PMujmQXn3ljQghxHmrVzLRTp05l8uTJl2x/Ab5GrNbqxz6ogH9ISKWv3W0w8rDRlwKg4VOn38/nZEhRUclbkccm259s4vT7vwwJKOPR8BLQKmTqTBw0RKFWM2K08YRx1dYM8N6KZWdsczam/3MrgKpYVS1+/V6oto0NHVZVy2+lHZlf2oUStfzOxU31h7nDZxn19Kl8X5DOsWPHqt2OikLEpIVnaFNz3nzjNdf3AQ4HH2v0BKYeJXfSE3weFF/hqplGXkF8EN6DhzKXs64snf9kr+P10C5k5RcRFhFZ6faHNNUx/WobT7az88q39/N+apMqa1FVldemv1bl6wAZOVXv66ScnNwz/kxTV31DVLvr0DbvQ0lwLF+sP8rsdYdxpmzHsfs31LR95OQXMLbn+Gq341SdZ9xXfn5+ta8DtGjdhvS0tGrbREZFsTN52xm3VRNatmlJ2hnqiYqKYse2HRe8Lx8/fyyW6v+PMhgMlBRdvCvJLkY9nvYzPRtnUzNAUXExfr6+F9ympo6/NrzXtSqojBs3jjGn3MDPbDYTGxt70fZntVrIfvmmatsEP/cTkxefHp68Mwq46umvwebgneJCvgkMrPD6qSHljkZ3MPHBiUx7796KG1FVOhfuo0vRPgDe32ih9TVRnOmylv+8dmu1rz//6Bd07VX1hx4ASzdW//o/Hrn66mpfH/vd14yYMqHK11UVZn68hGfzHiLLWT7IOFabwR2+y2ilP1jhUKMDAs9QjcpD/3my2hYzXppxhm2cvTGDOld4/vdf0bSd/xN3FJ4gsX93hv9PFmxtCOOt0O6MylrFktJjPJezASdOxn+1usp9/HF8Dt2OvstzcTtJaFOfhYG9K203d+lsrrn7mmrrVRcvqHZfAGP6NmPEGX6mE5bPJvnOa1FVOxstaXxbFMBaiw/a+DZo49vQUG+hcP4s+o8ZiqaaPttv7xt2xn1NXDGn2tcB0tPSznhcU+7pfsbt1JS0tLRK/0841QvXVx/ez5bFYuGND2ZV2+bpRx+pkX2djZqqx9N+pmfjbGoGeLpfC8bP23LBbWrq+GvDe12rTv0YDAZMJlOFh0dSVRp+thKtzUFesxiW2Cr+hfG/IWV8x9P/8tSoTvoUbHOFlDX+TRjzu+WMIaU2ybcZWZXXgLCbxpLlDCJIY+Zhv5+YEvgRrb0O1rpDzWzQgBNNm6KoKi0WL8K7kgPobIzi9dAu6FBYWHKEmIeicahVX92zOWYIL6wqP911d848BuQtuWj1nytFgY7eZbwVmsEP4ce51deMQXFywGYgtO9otvxl4limAZu9lv0ghRAepVYFldoifP1fhGxPwanTsH94jwqvnRpS7mx8J+M7jj9tTIreaWNg7npalRzFCSwJaMV6/8aX7gAuMotTxxZzLEtyG5Nl88dps3CLcTVvBL1Pd+8daJTaO9HY3mt6Uubnh29eHk9UcUqwp08sb4Z2Q48GU3t/lue+gV2teuzNi3/Y+W9QPwDuyvWssHJSgt7GfwJzWBh5jJGmXOyF2djsGo5lGdl8wMTBE0aKy+S/GyHEuXPr/xxFRUUkJyeTnJwMwOHDh0lOTiYlpZppxz2cV34xDT5fBcDRm9pTGv3vnCn/G1LGXXX6NPp+jlLuzP6DREsmVkXLT8GdSPa9PO7b41QVDhSHsyi7KX+XhgEKMYY80j4fxSDflXgrl2bekIvJ7u3NzuuuB+DewECK96dX2u4anxjeDeuO0+rkmGUzy3KnYXdWfW7/v8EDPD6sAARonAzzL+DExw/TIKYYP6MdVVXIzDew/ZCJ3Ud8ySvUIZPeCiHOlluDyubNm2nTpg1t2pRfDTNmzBjatGnDxIkT3VnW+VNVGn68Aq/CMgoTQjl6S3vXSz/+9aMrpNzV+K5KQ0qorYC7s1YRYS+gWGNgbkg3DnlHXeqjuCjSLCZ+z2nC9qIYbKqOQF0JPYIOkBR4GIc5y93l1aichERSWrUCIPXzP7CbSytt18UYTcrbx9Ep3pywbOf33JexOKu+L9H/hpWb8xbhsZ/4TgdhATZaJBbRPLGQEJMVUCko1rM3xY/kg/74tboem4eWL4TwHG4NKj169EBV1dMes2fPdmdZ5y3ij32EbTmMU6th74jeqLryy4eNHY0VQsp/rvrPaSGlZ4KWu7JXY3KWkqPz4+vQq8nwqubuu7WELrgOf+TVY01+fQod3hgUG+38j9IreB9hXjV7s0BPsr97Dw5ZrdgLSjnx2R+oVUzcVrK/lOuCJ6BXjKRb9/BL9njM9sp7YaA8rPwQNACAwbkLuCvnR88NK5SPYzH5OGgUW0LbBoVEh5Sh1aiUWrWE9H6Eb4sC2GgxUuSUcSxCiMrJSeMaYsgupMHs8pHTh2/rSHF8KAA7s3YSMDAAgHua3FNpSGHb1yy83RuDaueYVwjfhF5Nga76S9M8ncWhZ012K6KGvk26NQAFJw19Mugbupu6Pjm1bqDsuXJ4efFkWhqKQUfJ/nSyftleZdsIQ2P6hb6MjyaYAvsJfskeR6Z1f5Xt5wX348uQ8iu7BhQs48Gsr9HUgvfT28tJQmQZ7RsWkBhZgi0vDQsatlu9mVscwLJSXzIdl9/cQEKICyNBpSaoKo0/XIa+xIq5XgTHbmgLQHJmMmtS1wBwX7P7eLbDsxVDitMBSybCz4+i1yrsMcbwQ0iXWn3fHqcKu8x1+epYX7abG6JodUR55XNdyF5a+Z9Ar7lypl4/ZLMRdVcnAHIW7aRo14kq24boE7ghbBoh+kTKnGYWZb/A36Vrq2y/KLAnH4bdgxOFawrX8c0tXihq7XhvtVqICrGS+tlI+hiLiNLaUFH42+7FzyUmfi7255BNj9NzO4qEEJeQBJUaEL10F8E7j+HQa9n7aC9UrYatGVtZn7YegKJlRTzZ7smKIaXMXH6X3LVvAzBlrZWFge1r9Wyzx0vD+P5Eb1Zlt6PMaSBIX0DmD5PoGvS3W24a6AkCrqpLUPdGAKTO/gNbTtWnu3y1IfQLeYlYQ3sc2FiZ9yZbzd9CFb0lq0ydeTviAexoGdxUR4Pj29E43TeL7TlTncTrbAzwKWKgj5mGOgsaVDKdOpaX+TG3OABTx0FkFV6ZvztCiHISVC5Qgt6L+l+V95r8fWdniqOD2JS+iT/Ty2eY7RDZgaLfiyqGlJxD8Ekv+Ou38hvQ3fopE1ZZa+0cKQU2XxalJ/FzWg9yrIEYNFa6hmzj9pgllB1Ndnd5bhc+qD3e8SE4iq0c/3gVTlvV0+frNUZ6Bj9LM9/ycSjJRf9HwtOx2NTKZ/Lc5NeG16NGUGxVCSjJpUnKFvS22neX4xCtg6uNJdzpW0Bbr1K8FSfFqoaArnfRedoyHvt2GxsP51aY7VcIcWWQoHIBnHYnr0ZGo7XYyWsWw7HrWrIxfSObMzYD0DGyI+0j2ldc6dBy+PhayN4P/tFw3yJoMcgN1V84RW9kfW5zvj1+HX+XxKDgpIXpIPfELqJVwEG0tXg+lJqk0WuJefBqtL5elB3NIe2r9dV+4GoULR0D7uPqwCfQKd74NfblL/VzCtW/K22/06cp13xpwab1wsdSRNOUzRjLLs2U6TXNR6PSzlDGXb4F9PAuxpK6D5tDZcH2VAZ/uJ6+b//BVxuOUmSpRT1HQogLIkHlAuRvyqSRwRurycjukb1Zn7GBrZlbAegc3Zm2EW3/bex0wqrX4cuBUJYPddrDQyugTtvKN+7BnCqsLmtJ2PBZbM1vgkPVEmPM4PaYJXQP3Ya3B9440N30IX7UGd4dNArmjX+TvWjnGdep59ONm8JepzSlDAclHFG/J925Eqd6+of05jQne+LbU+rlg5fdQpOULZiKcy7GoVwSWgUa6K1kfjueXx7ryp1XxWLUa9mXXsjzP+2i05RlTPx5FwcyamcgE0KcPQkq56n4bzNF+/IB2D2yJ78Xb2J7VvmVHV2ju9IqrJWrbYCXCt/eASteBlRoOwSGLQT/6m8Q54n+ssXwQsFwPiy6Ga1fCCZdEX0j1nJj5GpCvDz7btbu5tskmsg7OgKQvSCZgo2V95CcKkAXzd+vHCWY8rmGstjAQfUzitXTb9Bo1RvZG9ceszEQreqg4fHthOcd8+jLl89G8zoBTB3Ykg3jezJxQFPqhvpSZLEzZ/1R+sxYjXf/cexKLcBqrx2DiYUQ56ZW3ZTQU1hzy8hdW363yY8LsjkYsJ/DuYdRUOgR24PGwf9Odx+af4Jlg2z/jkfpPx3a3OOu0s9bjsOf70p6stbSEgBvxULmym94ZJgJrSIfEGcrqFtDrBlmcpftIfWLtXTz8TnjOqpdpY7mOvzUBFLV37GQy9/q1wSrrYlUeqBVvF1tHVo9B2LbkJC+l1BzOvGZB/ApK+RoRKOLeFSXRoBRz/1dE7mvSwJrD+bw5YYjLN2bCeH1WbY3k9UHsmgY4U+zaBORJu/TpwEQQtRKElTOkdPiIHv5CVS7ijPOyM9t/PApOIxG0dA7vjd1A/6Z7l5VafH3epJ2L0JnAoISYPCXENXSneWfs1KnFwtLO/NraScseKGg0t2QzGDf5QzctALtfUPdXWKtEz6wHXZzKeZNh3kzIoLth/aQX6/pGdcLUBrhRzxp6gry2E4uyZjVv4jkmgpXBqmKhsORTSk1+BGTdZAwcxpGazF1/C+PD25FUejaIJSuDULJLCyj9S0jCO92J/mlNnanmtmdaibY14umUSYaRfrjZ5D/5oSozeRf8DlQnSrZK1OxF9ooCtXxxiANPgY/dBodfRP6EuMfA4C3pYhrt/4fCRnlk3YtOqzQd+wqMAa6sfpzpNGytLQd/y25GrPqB0BDXQpD/H4jUZfm5uJqN0WjED20C45SK+w6QftZL7J5xAvk121yxnW1ijcxSl8C1WacUBdjJZfj6i/UeyGBZGcRrRTf8p4ERSE9OJ4Sgx/1UnfhV2Zm83BvVhVs5URA7RsXVZVwf29sOxcx5Nn/kJpfxu7UAv7KLCK32Mqag9msPZhNbLAPunpJFFnsElqEqIVkjMpZUlWVvD8zKEstJj1UYcIQHXsMFhxFdm6sd6MrpMRk/sXty98hIWM/do2O1S1vYMhvuloTUlQVNlkaETbsfT4v7o9Z9SNSk8MT/t8zMWC2hJQaomg1xDx4NX+WlqIrK6X9+5MIOrj7rNf3U+JooNxPhHI1Ggx4x3rzjj2Nqfbj7HeWuq4qMvuGsDv+KkoMfkT4KQza9ShXHfsMRa36EunaSFEU6gQZ6dMskge6JXJNozCiArxRgZTcEgzdH6D9y0sYPXcbK/ZnYnPI6Uohagv58+IsFe7Jo2hfPgfqwOt36ijQ24my61g/7S8iFkSgddjouOd3Wh8qn1Ml1z+cJe3vICcgCvjNvcWfpb9sMXxT3IsD9jh0wWBSirnFZxXXem9FJ+NQapzGS8fItDQW9xtI6L5k2n8wie1DnyKzVaezW1/REU4SIbRh1cIJRPcL5aBaxqv24yQqBnprgmiv8QOv8kG2e35cyr0toEvKTGILNrG4wYsUG8Iu8lFeegadlpYxgbSMCaSg1Ma+dDPrtu2hLCCSn5NT+Tk5lUAfPX2aRtC3RRRd6oXipZO/2YTwVPKv8ywUHyogf2MmfzZUePFuHQV6lcY2A5/kxmBNsxCZc5TBK95xhZRdiR35oceof0KK5ztij+SNgjuYVHA/B+xxeGGjcP1cpge9Sx/jZgkpF1GZqrL1oefIbNYerc1Km0+nEb9ywTltQ6t4k/nfLKbpE+ihMaFD4bBq4SNHOmNtR1joyCVfURn6s5Xf6k/EpvEmrmAz9yTfRUJu1dP0Xw4CjHo6JoZQ+uNz/DyyC8M6JxDq50V+iY3vNx/nvs830e7lJYz5PpmlezIoq2YyPiGEe0iPyhmUHi8i+480FnZQ+LKnFlWBLmU+vFQQiY9T5dVeBm7540MUVIq8TaxqfQtHIxufecMe4Lg9lP+W9GCjtXwgp4KT7obtDPJZyaC1y/C54So3V3hlcOq92PbgeJr+8CGxa3+jyX8/wf/EEfbc9hBOL8NZbydI0TFEF8FNaggrnQWscBSQh53/OnL4yZFD3OMx/GIK4ljLz7jxr0mEFx/glr1PkBx5G2sSRmHTnvkKpNqsVWwgrWIDmTCgKRsP57JoVxqLdqWTVWjhx60n+HHrCfwMOro1COWaxuH0aBRGuL/3mTcshLioJKhUo623kdTVJ/ion4ZVLcs7n24pMfGUOQwvRx4G2wGe6GQAVPbFtmVtiwFYvIzuLfosaAOj+KDwZtZZWqCioKCSZNjFQJ9VRGlz3V3eFUnVatl9+whKQiJouOBLYjYsxXTsENvve/qctxWg6LhJG0I/TRAbnUUsc+ZzRLXg38qPlXkzWKsY+T22AyNzIrk6azWt038gMW8tvzeYeBGOzPNoNQpJ9UJIqhfCpBuasSUlj193prF4VzppBWUs2pXOol3pALSMCeCaRuFc2zicFnUC0NSG21QLcZmRoFKFkk2bmNwwhomDNByKVtCq8FhhKHcUe2Ow7UHnyADghNlJcp/7akUvSr7Vj60FjQm7fyBrLeU3P+zgtZdbfVYSq8tyc3UCReFw71sxx9aj5RdvYjpxmM7TnmB4YCClTieq5tzO1OoVDV20JrpoTaSpVh76cQuJNzelyJHJXssfjPKDTtpIXs7OJ8KSym27HuHd6/WUqGXYlCujJ0GjUeiQEEyHhGAm9G/KzhMFLN+XyYr9mew4XuB6vL3sL0L9vEiqF0pS3RA61wshPsRH5moR4hKQoFKJki1bWDLpIV6/X0+Bn4LJoeGV/AiSSvPwsm1HwYEK2LV1aPPhXiYN8eyQkmUJZEt+Yw4VxwAKigZa6/9ikO9KuYrHA+U0bs26sTNo/s27hO3dxlOhoWR8+gnLru3Jobp1z+vmlVGKF5k/ZfPsQx+QYd3H0bI/OW7ZwgZjKjfVCeWp3HxuKyxiVHs96c7n+MaYxDGSCLSGo0N/EY7S82g0iuv00JO9G5JpLmPl/iyW78tkzcFssousLNieyoLtqQBEB3jTqV4IneuF0jExmJggowQXIS4CCSqVWHBwAVNvteHQKtSzevFejjdRlt1o1SIAHIo/Vq9GODUmCix73Vxt5VQV0spC2ZLfmJTSfwf1Jviksvnj6TwzOsCN1YkzsQSGsGXEC9T5czlxc2YQkZnJXXO/5Uh8POs7duJg/frnFVgURSHS0IRIQxM6MgyzPZ3jZVv51LiVZfmbmZiVTrSjjDHFK1hpXM+r4cEUKpEEWSMJtkYSZItA8boyPozDTd4M7hDL4A6xWO1OtqXksf7vHNYdymFbSh6pBWWusS0AYf4G2sQGom/ej5RMA5HBVrx0tfv2BUJ4Agkqlah7dX/4bR4x23L4PliPl7N83IaKDqu+LnZtnfP6kLgkFA1/F0ezLb8R6ZbQ8kWo1PdLoV3gPkK8zKxP2w/IQFmPpyic6NSTh156gjfvvIsOmzeRcPQoCUePkhMczM7mLdjdtCm5ISHnvQuTLpKmfv1o6tcPZ7CDFjOuYvawllxv2U2P0jKSUlP5JKCYzwMyOeK7C4Cm7zfkXvPHxOUoxGUpxOZqiM/TElGkRato0Xh781pEBMHLllLk54/ZZCIrLIyc4OBzPn3lKbx0GjrWDaFj3RCe6AWlVgebj+ay/lB5cNl1ooCsQgu/78nAq/1g5i4HRVEJC7ARGWwhPNBGWKCVsEAr3l4SXoQ4FxJUKtHBN5YVxhaYTPPROhVUFOzaGKz6BFA8sxu8yOnNyrI2hD3wEIsyym92qMFBE/8jtAncT4C+2M0VivNV4HSytFcvNnbowFWbNtEmeRshubn0WL2KHqtXkRsURNPwcHSb83CGGVBDvFB9tFDZwE+nA52lDH1xEd55WRjzsv75mo13bhazrWHEf11EijGUiHYF+EVaGZlfwD3ZRfygmPgq3Jcck4YTgU5OBML6ev9u2sumEpMNsVkqcXVM+GT8SeNdKkFF5TP823Q6skLDyIgIJzUqmvV6fXnXn6eG/moYvbR0axBGtwbl89CU2RzsOlHAtpR8Jn8wh6AmbSgs1ZGZ70VmvleFdU0+dgzXjmb67/upF+ZHQqgviSG+BPh45v8tQribBJXKbPuSoD0LQKNg14Rh1ddD1XjmpZsp9nB+L+3AWktLrOjRBYBBY6GZ/2FaBvyFr67M3SWKGmIOCGBpr16s7taNJvv20nTPXhKPHCY4L49bTSZYmO5qqwIYtag6BbQK2JxsqlsXn9EDq91HmMEAdivWQj0H18XgU9dBXON0ArwtPEAeg9LgzW8LaN/7Ro74lXHUu4jDXgWkeBVi0Tv5Owr+jqoYPHwsCnGZKnGZTuKy0onLSqPHvu30j4+n9L/vkBWZQFZkPOl16mMx+l2Ed+7i89ZraZ8QTPuEYMYPfp8R/5lMYYmW1BwDmXleZObrycr3wlyiw1yiQxfXhneXH6ywjSAfvSu0JIT6Eh1oJDrAm8gAb6ICjBi9tG46utrD4VSx2p2U2RxY7E4s9n++2k753u7AYnNic6o4nSoOp4pD/ed79dRluJapannu12oUFEU5/XtFQVe/C/vTC9FqFHQaBe0/j5Pf67QadBoFRe+NqqoynukcSFCpTKeRkJrMNf/5jl9GXOvuak5jVbVssTZmeVlb9tgSXcvjtOnsWPgdYx+JQaeRSdouV1aDge2tWrO9VWsMZWXEHj9O8ayZDOsQgybLilJgK79HYanj1HsV4nPKaReHTk9ZUBilQaGUBYf9830Yb73+H24f8RKlPv449OU9AVrVRnvnKro4fyPQK48X+4BVl0yJ9lpsal2wKjisTtI0hRzW5HJYk8vbO5cT3iERC7mUGFT2xcK+2IqnfYLNKvFZJcRm7SEuazdN/lDxUaLJim5MWmxDikznf0rLE/j7OGjkU0Kj2BLXsjKrQla+F7OnLGHY4+P5O6uYw9nFZBZayCuxkZeSz7aU/Eq3F+ijJyrASKTJQLCvgWBfPf5XDWSX2RujVsWoceKtdaLXqHgpKnqNWmmn2qWmquBQwaYq2FUFXVA0u04UUGZzUGpzUGot/1rm+t6JvvVN/PFXFnaHis3pxO5QsTtVbI6T3zsxDnyFzlOX/RM+ygOIzeG+02qGbvezeHf6GdvFPP4N7y4/iF6nwUurwauSrwHdh7D5SC7eei1GLy3eun++6jV467RX3GXyElQq4+UDd3zNhnu/dXclLqoK+sgGfF7Ul3WW5pSo5fO1aHDS3msffYwbaaxLoeeujeg0ckfjK4XF25uD9eszNyeHO+7uXL7QqaKUOFBKHOWfEA4VdAq3jFvBs9+sxuFtxKnTV3rKZX1pKf0DKgYEh6LnT20vtms609n5G22sS/DmCF7mz7DpYikx9sCqb0SMM4AYZwDdSOTpWV/RteMUnKodCzmUkUWZmkUZWVjIxoaZXJNCrklh2ymnjxRnBlF5GcRlriQ81YeBXSM4eGAjifXbodXU/h4Fby+V2HAL9n3LmDpwrmt5scXOkZzy0HIku5ijOSWkFZSRVlBKWkEZJVYH+SU28kts7D3lQr3AbveyrJqZBXSKSvQjn9Hj9RV467XotAp6reafx+nfK5T3xqnqya9qee+cCirlPQuOfwKD1eHEZlfLvzqcWO1OIu9/n8+Ohv7za6fgUBXsKpx6e++o+99nwLtrqn2fvNrcyNYqAttJmoBIUguq7jHWaRQMOg0Gvbb8q06DQVf+Ye+l06DTaNBqFDQaBe0/PSQaRTllmeJapijg/KdnxamW97ic+r1ThcVLlpDQomN5b8w/D7vrq7P8q6P8/VQBq738PcNyeu2mDjez9lBOlcdm0Gkw3vISd328gXB/A2H+BsL9vQk3/ft9hMmAv/flcTpRgoqHK3EYOFAYx97CRELvuY2l//y7DNEU0N2wnR7eWwnVmt1bpPAsGgXVT4fqV/Gf93G7HZv/+V/tVab4sFx7C4PfX8D+sddiLNuE3n6MgMIvsWsjKTF2w+LVHJR/96tRdBiJwEjEqZ9VONQyfpnyJIMfeZACfTb5+iwKdJlYtRZSQyA1RAHK4KoAblk/HK81Com6CBrGtKZRRHPqB9anQVADwoyXx72KfA06mkUH0Cz69J+PqqqYy+ykF5SRWlBKprmM3GIbeSVWZnzwCc2u6kKZQ6HUoaHMqcHmVHD+82bbVQWtbxBHckpO2+7FoA+KptBe9etaVGxlxUSFh2DUa109BkZ9+cP7n++//WoOHa8biE6rQX/ytIlWQa8p/6rTKHz7ymMsXbQQg/7fEFIeTMp7JnTaSztwO+zZAdw6dHW1bVRV5dmbOjDphw1Y/wl3VvvJ0OfE4ijvHfrtu8/oeMM9lNmclFodrt4ni728p9xid6IJjGZdNWEGwN+gIzrQSFSgN9GBRuoEGokO9CY6wEh0oJHIgNoxX5IEFQ9U6vDi7+I6HCqO4XhpOOo/t2RS7VY6++7nasN2mukPo1Hk6gFx6Z0oVCn2HUCJsQc+pWvwLvsTnSMdU9EPOJVFlHp3IPIMQ020ijclB0upX9zGtUxFpUxTRL4+myJNGhbrITKdx8gP12DVw341nf3HFrPg2GLXOiYvE8GPBLAu/2OC9LEE6eIJ0sdh0PherMO/5BRFIcCoJ8Cop1Gkf4XXJtz0Ljfd1KzCspOnWqyqgs2pYdrLr/DHug3/jMso/0C0/XP6xPrP97Z/ekVOHdusKAoK5c+Vk88V0ChKhdMUJ3tkvHQael57LY8//XR5D4WiolXKe3X0GhWdUn4q6ulHHyHFZq32mL8Y9TXdRj1cbRtn5kFaxNSuaRYURUG1W/A16KjuN/T7VV/QZ+zps1I7nSpl9vJTZB9PfowPP/+azMIyMs0WMgstZBaWkVVY/n1hmZ1Ci539GYXszyisoh4wDn6d7zYdw2TUYfIu/z0z/fP75mfQofWA00wSVDyExieQXea6HCqO4URpmCucAIQbcmnif5jvXniNUZObVbMVIS4dVeNHse/1lBi7Yyz7E++yP9GqhfiWriDlMSObMj5jpX8Su42NUJUz/3WroGB0+mO0+AOJQGcmThpO8bEj7F80l92bF3Go7DgpYQrHwhTSghXMVjNedb3YV7K4wrZ8NSEE6eMI0sURpI9DF63D4rBg0J79vZNqK0UBnQI6VNA6sGUdoUNC8CXZtzV1H5He1XSpiAui0Sj4eOnw8dLhTNvHzW3qVNm2xGonraCM1PxSUvNLOZH/7/ep+aWkFpRhtTv/v707j46qvhs//r4zk9mX7BuEJUFEqwKShoKy2BMV5Si21scuetyKTz3YPpbHtlrb8mtV5FFsqRarpQFsrUIfD92sVX/yEy2KIpQoBYwGA0KWIQlZZp/M3O/vjwljYiAhMckk8fPifM+dufc7937uZyYzH+6KwZFJQ3uYhpNsmNe0xFYZ8+f/YwjXqm9SqKSIUooD9T5efb+RV6qOUXjrBl5t+ngffI65hRLnEUoctaSnJS40tykipxiLkUcZ7ATtFxG0zccS3Yct/CZpHGaufxdz/btoNqaz3TWb11xfoN6c1+/5W/MLmH7Td5l+03eJHDxI23PP0f635wjUH6UuCz7K0agZb6N6sofa9A78WisBvZlApJmjkT0AZN+RQdkfypjgmsAZGWdQ5CpinHMcBY6CxNBZgM008u/TJcTpsptNlOQ4Kck5+eZNpRTNgSifK5vP1fc8hi8Uoy3cQVuog/ZQB+3hGHE9sdsRc2q3UEqhMozagh38s7qRV6saefX9Ro75Pj6KSjMYybUcp8RxlBLHUbnuiRh9NCMRy3lELOdxyarl/M9/L2Ku/22y4q0saX2RJa0vUmMuYqdzBm87ZvY9v5OwlJSQ+1//Rc53vkNoTyW5f/srGX94mvn7gkDiOAzvhHG8O+dcqqbm0GhuoSX2EXWt/wY7HGo/xKH2Qyedd6Y1kwJHATm2HLJsWWTbsru1LFsWWdYs7Gkj81IFpyUWhagfYmHoCEEsArEQdIQT4/R4Z8fO3cqqy+5lowlM1s5mAaMlMTRZwepODMWooWka2U4LevNhzsh19ZiulCIQjdMe6uDJP78EPDD8QXaSQmUIhSIGjjZaMJdey5K1r7P3aCt6l797W5qRuSVZLDgzh29eMYdlKy9NXbBCDKJKr2JjzrU8lf1lzg/sZb7vTaYH9zM5eoTJx49w7fG/cfN/WmmJ/4WD2tnUapPRtdP/OtI0Dfv5M7GfP5PzH/s1D/3gIQrffpXcvTvJ+6iWiz+q5WLg+JTPUVe6gB8+/hqvf/gO1S3VVLdWU+uvpc5fR12gjjp/Hf4OP8fDxzke7vvu4XaTnUxr5sfNlkmGJYNMayYZ1gzMZ5hpDDZiS7NhM9qG9GwlixE8sSac8TYc8Tac8Xacehv2eDsWPYRVhbjkVgs8PA3CbdAxhAfVGtKov8NKWsP9RAw2ogYrYc1O0OgkaHARNLoIGlx8YZwBjteAMzfl/1MXp6ZpGk6LCafFhGqpTWksUqgMEqWgxW+ivtnC0cZEa25PXIci7ZzLeOdIKwBT85wsmJrDgqm5fH5yBhZT4kvspvZjqQpdiCET09LY6Tyfnc7zccV9zAq8S1mgknOC73F2DqC/xAW8RAQLh7Wp1GjTOGw4gyYKTuu4FoAY0HjubBrPnY0pFCCvcgeFb28js/rfZFbvI7N6H3/2eAjeupwzy8spLb8Y8+e679tvj7ZT56+j3l9PU7iJplATzaFmmkPNNIWaki0cDxOMBQn6gxz1Hz1pPJnfzOTZD55NPjcbzNhMNmwmG1aTFZvJhvNSJ0/tf4oMa0a3oifDmoHJYAJdh2Az+OqgvR58na29rnOYeB68ywbe+3vNT2GOIfGargxdt45YIa1zaOjyk5A8fb3zpOV4DOKRxJaX2IlhNDFEgd5BrkODeBPEOaUrb7TAIzMST9Ic4MwB93jwjAPPeHCP4+IJcbIC1fgteUSMzlF59WIxeKRQGYC40jike9gfy8a18Js8szUPb6uZaEfPL9Ysd5SGt1/nVz+5g9mTsyhMl/3g4rPJZ3SxzX0B29wXYIuHCK+5nbu/dAGT1Xs48DNV7WWq2gs6RLBSq01Cn2+C91+C/HPBld/nD1bM5qB2Tjm1c8qxtjRSsPuf5O95Hc9H1YR27Sa0azfHVv0PljPOwHHhhTjnXYittBS32Y070820zFPfCV0pRTAWpCnUREu4heZwMy3hluSWmBNt++7teCZ4CMVCKBRRPUo0GqUt2gaAVdc5e56Fra/9lNx4nNxYnNx4nLxYjNy4Tr6uyI51YFKnd1ZfDBN+o6dbCxrchA12IgYbv3zkcV5+7U2wesCWDmYnGAfx+hq6ntidFGln+lkl/PCH38Oih7CoMFY9gF33Y4/7sOs+7HEfkcZDTMq2d+5+CkBLAFoOdZvl04uByq8BEDE68Jtz8Vny8Fnyabfk47PkM7cgnnide9zgro8YcaRQ6Y3BxMF4Oh+eaHpiWBP3ECLxh+EshSOdF1wyGXVy0jsYnx1mfE6EcTkR7BadFY88yZdmrk7higgxsoSMNp75d5yzr7kJlE4etRTrB5ikqhinarAQpli9x4/nmeDpaxIvsmVC3ucSLfsMyJgMmZMxGU7+gx7OyKGm/MvUlH+Z316/gH8++CC+l18muGsXkQ8+IPLBBxzfsAHNZsP++VLspZ/HPut8rOeei8Fs7jE/TdNwpDlwpDmY6J7YfaJSiR9r/zGu+H4Zy+6/CHvYhznUii3Ugj3UjjPswxMJYI/3fmpuV80GA16TkWNGI8dMJhqNRo6bLXQ4ctBc49jy4g5mlC3CiQencuFUbhw4MXb5an/lsA6FM057mf1mMCSOUbG6+Xejos5S3Gv3O//Pt4hH25P5wu+FtlpoP5oYth1l7+svMKXAjS3WhiUewBKqIStU020+i64Cfjkd0MBVAOlFiS0yni7DE+Oso+s0ZtGdFCon8Y+99Tz0UhW53/lfrmw/eYpsdDDN2Mz2t3dyzW1zycuIkunuYJivMSTE6KcZ8FKE11jEDi5BUzo51DFe/5Cmyqe5/otnQXM1hI7DoX8mWhdHlkJg1xIC5mwC5iwCaVkEzNmETR4iJidRo4Pi/BiZF08nc9HnifkDhCrfJfj2vwjs3EWsqZnwW9uI7NpGq0FhsJiwTivBOqUY88TxWMblYvKY0SI+CLUmjvUIdw4DjYkfWv+x5PEff7sKePvpXlc50AGOvBJwFyZ+ZF356K58/FY3jaY06gw6H+kRasONNAQakq0xdOIytD7wvQdzM6jkre4zV2DDgVO5cSoXuVdn88x7z1DoKKTAWUChoxCnOcX3VNI0sLgSLaukx+QvfjefHz71MqZ4GFekAVfUizPixR3xJp5HGogdfpuSLHNid5SvLtGOvHWShQEWT6Jg+WQxkz4h8R4482SrzAg2IgqVtWvX8tBDD9HQ0MD06dN59NFHKSsrS1k8BoPGh40BNKMJO1EmG9soMbZSbGil2JhoEwztGDVF5v/7M+c8OD1lsQox1ijNwDHGc8w4np/87Umu3/J24gyVxirw7oNj+6H5ILTUQMshTLEwnkgdnkjdKed51VXAE/OAxJeeC3BZgHmnekUdRP4J75Nop8vs5MNGP7bJEwlaXQQtToJWN36rm4DNTcCaaHdf+QCNDf/q9lID4O5sPX+6E6LxKN6gN1m4fPO/l3LhlRcR0Hz4tXb8+IhrMUIECGkBGqkn+7IsVr61stt8XGYXBY6CZPGSa88ly9r9TKfk8TIpFDNaabFPosU+qce0lavm09hQnygW245A6xFoO5p4fGLYeiRR4Eba4FgbHNt36oXZs8CZD668LsMuzZ4F9kywZUhRM8xSXqhs3ryZ5cuX8/jjjzN79mzWrFnDpZdeSlVVFbm5uSmJqWxSJr+7uYzLvvA5vN+/QI7jEiLV0myJ3Ref3IWh65xTUsCPVv8KR0czjmgT9mgzjmgz1ng7lpgfc9xP4Oh7lIzPBT3W2eKJYbwDVOeRn4Y0lMkCWhoqDnpMR4/EiQc7iEdAj2rEowbiHYbE4w4D8bAB3ZyB5ilAyyzClDeeX25Zx/w7F9JhttHhsBBzWhNDhyWxmwTodj+BfjAbzRS5iihyFQHQ9PfjXLj44uR0hSJMiEBn0eLX2vnHK3/mmpu/Qn2gnrpAHW2RNnxRH76oj/dbTl2FaWhkWDOShYvH4kkcy2N2d3vstriTjw12A3HiGDCgDXAd+0XTEmcPOXNh3KyT94kGuhcuXYuZ1iOJLTF6LHHwcrC592LmBIsH7BmJ4sWWCfYs7rsgyvTDjxM1OYgaHUSMia15J54nmh17GqB0OM2DxcUIKFR+/vOfs3TpUm666SYAHn/8cf7+97+zfv167rrrrpTElOEwM39qDrqvSYoUIUYygwFvUKPe3ftWzZUPzafRe4ofZaUSzWCg63kuJ35GVDxOtKaGcFUVsZpDRD/8kMihGqK1h1ChEBAGajobfMvmgLX/96SL6rCbiTmsPOlKp+bqr6BZrRis1sTQYkGzWTFYEs81oxFMRjSDEc1khBNDY2IcJiNfcbkYt/fdT5yhAyr5vWWh+e8hfnbNFxN9HBpBewSv3kZDvBWvasWrt9Ks+ziu/BzXfTTrPlp0PzoqeYBwbwVNV9MeOYMneQRNaZh1I2m6EZNuxKwSj9N0IzO/O4Flm76GSTNiwkiaZsSkmTBhIK1zePYiFx+8/3OMXQoeTWl0/XfuPAcb/7QCjc4bB3b+M/CJxyduBZB8d81gLIHMEshMvP/mWBBrLIC1w/9xi3382BILYo6FMMdDiblE2hKty0HA/3kecLSizxzd9n0rvDGbGAY6NBMxTJ1DIzEt8fjC2+y4dl5DHAO6ZkAn0eJdHj++RPHBY4tQWqKP0owoDMmz5ZSWWOvEPSE1VJdPd9dpd19sYHLlDzpvlvjx+BP9FRpfmp7aoiqlhUo0GmX37t3cfffdyXEGg4Hy8nJ27NjRo38kEiES+fgiaW1tiaPo29uH5qZ8iZuBdfTRB0L+UK99dF3vM0alFOFQ7wfZKQWBcC/n/XX2Gaz5DM6yFJHwqe9weqJP38tSRMInuc3oJwIarD4dfcSMUoRDvb/vSinaA32vezjg7zueUO/Xv1BK4Q/29Vk9vWWFw32vl6+PdT+dmE93WX397ei63ud6nc7fYK9yc9Fyc7HMm4eFxK4jpRTx1lZix47R4fUS83qJHWvkj0/8mtLzijEGw6T5o5iCYUyhzvfGFwJfiBw0mt59d+DxdLozKwv+8Xyvfe7Jzub9O7/XY3xWZzv7JK/RNfDboNUObU6NNjv4rRC0agSsELBAwApBi0bQCn4LBK0QM338P7sYn7iUvgYYgWIzr7S+0/uKzUljG73f4I9yCw96/7f3Pp+GsbNZAWyADYNSuHQdT1zHo+ukdw49uo5b17HrOg6l49AVDj0xtOsKp9Kx6wobXQ/8jne2CCa6/xhnOAB/94OHP+nsLODIG596NYuLgONbe+3jPM866L+zJ+anTufsNpVCtbW1ClBvvPFGt/Hf+973VFlZWY/+K1asOHGHbGnSpEmTJk3aKG9Hjhzps1ZI+a6f/rj77rtZvnx58rmu6xw/fpysrCw0TaO9vZ2ioiKOHDmC2+1OYaSfTZL/1JL8p5bkP7Uk/6kzkNwrpfD5fBQWFvbZN6WFSnZ2NkajEa/X22281+slPz+/R3+LxYLF0v3up+np6T36ud1u+aCmkOQ/tST/qSX5Ty3Jf+r0N/cej+e0+qX0CBmz2cysWbPYuvXj/WO6rrN161bmzJmTwsiEEEIIMRKkfNfP8uXLueGGGygtLaWsrIw1a9YQCASSZwEJIYQQ4rMr5YXKtddeS2NjIz/5yU9oaGhgxowZvPDCC+Tl5fV7XhaLhRUrVvTYPSSGh+Q/tST/qSX5Ty3Jf+oMde41pU7zzldCCCGEEMNMLo0nhBBCiBFLChUhhBBCjFhSqAghhBBixJJCRQghhBAj1qgrVNauXcukSZOwWq3Mnj2bnTt3nrLvvn37uPrqq5k0aRKaprFmzZrhC3SM6k/+161bx7x588jIyCAjI4Py8vJe+4u+9Sf/W7ZsobS0lPT0dBwOBzNmzOD3v//9MEY79vQn/11t2rQJTdO46qqrhjbAMaw/ud+4cWPiZoRdmtVqHcZox57+fvZbW1tZtmwZBQUFWCwWpk6dyvPP935fqlManLv2DI9NmzYps9ms1q9fr/bt26eWLl2q0tPTldfrPWn/nTt3qjvvvFM988wzKj8/X/3iF78Y3oDHmP7m/+tf/7pau3at2rNnjzpw4IC68cYblcfjUUePHh3myMeG/ub/lVdeUVu2bFH79+9X1dXVas2aNcpoNKoXXnhhmCMfG/qb/xNqamrUuHHj1Lx589SSJUuGJ9gxpr+537Bhg3K73aq+vj7ZGhoahjnqsaO/+Y9EIqq0tFRdfvnlavv27aqmpkZt27ZNVVZWDmj5o6pQKSsrU8uWLUs+j8fjqrCwUD3wwAN9vnbixIlSqHxKnyb/SikVi8WUy+VSTz755FCFOKZ92vwrpdTMmTPVj370o6EIb8wbSP5jsZiaO3eu+u1vf6tuuOEGKVQGqL+537Bhg/J4PMMU3djX3/z/+te/VsXFxSoajQ7K8kfNrp9oNMru3bspLy9PjjMYDJSXl7Njx44URvbZMBj5DwaDdHR0kJmZOVRhjlmfNv9KKbZu3UpVVRXz588fylDHpIHm/2c/+xm5ubnccsstwxHmmDTQ3Pv9fiZOnEhRURFLlixh3759wxHumDOQ/P/1r39lzpw5LFu2jLy8PM455xxWrlxJPB4fUAyjplBpamoiHo/3uGJtXl4eDQ0NKYrqs2Mw8v+DH/yAwsLCbh94cXoGmv+2tjacTidms5nFixfz6KOPcvHFFw91uGPOQPK/fft2KioqWLdu3XCEOGYNJPdnnnkm69ev5y9/+QtPPfUUuq4zd+5cjh49OhwhjykDyf+HH37Is88+Szwe5/nnn+fHP/4xDz/8MPfdd9+AYkj5JfTFZ8OqVavYtGkT27Ztk4PahpHL5aKyshK/38/WrVtZvnw5xcXFLFy4MNWhjWk+n4/rr7+edevWkZ2dnepwPnPmzJnT7ca2c+fO5ayzzuKJJ57g3nvvTWFknw26rpObm8tvfvMbjEYjs2bNora2loceeogVK1b0e36jplDJzs7GaDTi9Xq7jfd6veTn56coqs+OT5P/1atXs2rVKl5++WXOO++8oQxzzBpo/g0GA1OmTAFgxowZHDhwgAceeEAKlX7qb/4PHjzIoUOHuOKKK5LjdF0HwGQyUVVVRUlJydAGPUYMxnd/WloaM2fOpLq6eihCHNMGkv+CggLS0tIwGo3JcWeddRYNDQ1Eo1HMZnO/Yhg1u37MZjOzZs1i69atyXG6rrN169ZulbMYGgPN/4MPPsi9997LCy+8QGlp6XCEOiYN1udf13UikchQhDim9Tf/06ZNY+/evVRWVibblVdeyUUXXURlZSVFRUXDGf6oNhif/Xg8zt69eykoKBiqMMesgeT/ggsuoLq6OlmcA7z//vsUFBT0u0gBRt/pyRaLRW3cuFHt379f3XrrrSo9PT152tn111+v7rrrrmT/SCSi9uzZo/bs2aMKCgrUnXfeqfbs2aM++OCDVK3CqNbf/K9atUqZzWb17LPPdjtN0OfzpWoVRrX+5n/lypXqpZdeUgcPHlT79+9Xq1evViaTSa1bty5VqzCq9Tf/nyRn/Qxcf3P/05/+VL344ovq4MGDavfu3eqrX/2qslqtat++falahVGtv/n/6KOPlMvlUrfffruqqqpSzz33nMrNzVX33XffgJY/qgoVpZR69NFH1YQJE5TZbFZlZWXqzTffTE5bsGCBuuGGG5LPa2pqFNCjLViwYPgDHyP6k/+JEyeeNP8rVqwY/sDHiP7k/5577lFTpkxRVqtVZWRkqDlz5qhNmzalIOqxoz/5/yQpVD6d/uT+jjvuSPbNy8tTl19+ufrXv/6VgqjHjv5+9t944w01e/ZsZbFYVHFxsbr//vtVLBYb0LI1pZTq/3YYIYQQQoihN2qOURFCCCHEZ48UKkIIIYQYsaRQEUIIIcSIJYWKEEIIIUYsKVSEEEIIMWJJoSKEEEKIEUsKFSGEEEKMWFKoCCGEEGLEkkJFCCGEECOWFCpCiEF34403omka3/rWt3pMW7ZsGZqmceONN3Ybv2PHDoxGI4sXL+7xmkOHDqFpWrJlZWVxySWXsGfPnmSfhQsXdutzop0sBiHE6CGFihBiSBQVFbFp0yZCoVByXDgc5umnn2bChAk9+ldUVPDtb3+b1157jbq6upPO8+WXX6a+vp4XX3wRv9/PZZddRmtra3L60qVLqa+v79YefPDBQV83IcTwkUJFCDEkzj//fIqKitiyZUty3JYtW5gwYQIzZ87s1tfv97N582Zuu+02Fi9ezMaNG086z6ysLPLz8yktLWX16tV4vV7eeuut5HS73U5+fn635na7h2T9hBDDQwoVIcSQufnmm9mwYUPy+fr167npppt69PvjH//ItGnTOPPMM7nuuutYv349fd0v1WazARCNRgc3aCHEiCKFihBiyFx33XVs376dw4cPc/jwYV5//XWuu+66Hv0qKiqS4xctWkRbWxuvvvrqKefb2trKvffei9PppKysLDn+sccew+l0dmt/+MMfBn/FhBDDxpTqAIQQY1dOTk5yV45SisWLF5Odnd2tT1VVFTt37uRPf/oTACaTiWuvvZaKigoWLlzYre/cuXMxGAwEAgGKi4vZvHkzeXl5yenf+MY3uOeee7q9put0IcToI4WKEGJI3Xzzzdx+++0ArF27tsf0iooKYrEYhYWFyXFKKSwWC7/61a/weDzJ8Zs3b+bss88mKyuL9PT0HvPyeDxMmTJl8FdCCJEysutHCDGkFi1aRDQapaOjg0svvbTbtFgsxu9+9zsefvhhKisrk+2dd96hsLCQZ555plv/oqIiSkpKTlqkCCHGJtmiIoQYUkajkQMHDiQfd/Xcc8/R0tLCLbfc0m3LCcDVV19NRUVFv66DEgwGaWho6DbOYrGQkZExwOiFEKkmW1SEEEPO7Xaf9DThiooKysvLexQpkChUdu3axbvvvnvay1m3bh0FBQXd2te+9rVPFbsQIrU01dc5gEIIIYQQKSJbVIQQQggxYkmhIoQQQogRSwoVIYQQQoxYUqgIIYQQYsSSQkUIIYQQI5YUKkIIIYQYsaRQEUIIIcSIJYWKEEIIIUYsKVSEEEIIMWJJoSKEEEKIEUsKFSGEEEKMWP8fuq1T5fqxvC8AAAAASUVORK5CYII=", + "text/plain": [ + "<Figure size 640x480 with 1 Axes>" + ] + }, + "metadata": {}, + "output_type": "display_data" } ], "source": [ - "# For each city, test torch model on that city, then fine-tune and re-test\n", - "model = model_utils.load_model(\"../logs/\", \"kcm\", \"GRU\", 0)\n", + "plot_df = pd.concat([metrics_df, metrics_tuned_df], axis=0)\n", + "plot_df = plot_df[plot_df['n_batches'].isin(['0_batches', '1_batches', '2_batches', '4_batches'])]\n", + "plot_df['Tuning Sample Size'] = plot_df['n_batches'].replace({'0_batches': 'No Tuning', '1_batches': '1000 Samples', '2_batches': '2000 Samples', '4_batches': '4000 Samples'})\n", + "plot_df['MAPE'] = plot_df['mape']\n", + "sns.histplot(plot_df, x='MAPE', hue='Tuning Sample Size', bins=50, kde=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# loader = GTFSLoader()\n", + "# regionalizer = H3Regionalizer(resolution=11)\n", + "# joiner = IntersectionJoiner()\n", + "\n", + "# area = geocode_to_region_gdf(\"Seattle\")\n", + "# gtfs_file = Path(\"seattle_gtfs.zip\")\n", + "# # download_file(\"https://metro.kingcounty.gov/GTFS/google_transit.zip\", gtfs_file.as_posix())\n", + "\n", + "# features = loader.load(gtfs_file)\n", + "# # Need to name index to match, and fill stops with NAs\n", + "# features.index = pd.Index(name=\"stop_id\", data=features.index)\n", + "# for colname in features.columns:\n", + "# if colname.startswith(\"directions_at_\"):\n", + "# for row in features.loc[features[colname].isnull(), colname].index:\n", + "# features.at[row, colname] = set()\n", "\n", - "model.eval()\n", "\n", - "# gdf = self.gdf.copy()\n", - "# gdf['t_min_of_day'] = self.traj_attr['t_min_of_day']\n", - "# gdf['t_day_of_week'] = self.traj_attr['t_day_of_week']\n", - "# # Fill modeling features with -1 if not added to trajectory gdf\n", - "# for col in data_loader.NUM_FEAT_COLS:\n", - "# if col not in gdf.columns:\n", - "# gdf[col] = -1\n", - "# feats_n = gdf[data_loader.NUM_FEAT_COLS].to_numpy().astype('int32')\n", - "# return {0: {'feats_n': feats_n}}\n", + "# regions = regionalizer.transform(area)\n", "\n", - "# data_loader.normalize_samples(samples, model.config)\n", - "# dataset = data_loader.H5Dataset(samples)\n", - "# if model.is_nn:\n", - "# loader = DataLoader(\n", - "# dataset,\n", - "# collate_fn=model.collate_fn,\n", - "# batch_size=model.batch_size,\n", - "# shuffle=False,\n", - "# drop_last=False,\n", - "# num_workers=0,\n", - "# pin_memory=False\n", - "# )\n", - "# trainer = pl.Trainer(\n", - "# accelerator='cpu',\n", - "# logger=False,\n", - "# inference_mode=True,\n", - "# enable_progress_bar=False,\n", - "# enable_checkpointing=False,\n", - "# enable_model_summary=False\n", - "# )\n", - "# preds_and_labels = trainer.predict(model=model, dataloaders=loader)\n", - "# else:\n", - "# preds_and_labels = model.predict(dataset, 'h')\n", - "# preds = [x['preds_raw'].flatten() for x in preds_and_labels][0]\n", - "# preds[0] = 0\n", - "# self.pred_time_s = preds\n", - "# self.gdf['cumul_time_s'] = np.cumsum(preds)" + "# joint = joiner.transform(regions, features)\n", + "# joint.index = pd.MultiIndex.from_tuples(joint.index, names=[\"region_id\", \"stop_id\"])\n", + "\n", + "# embedder = GTFS2VecEmbedder()\n", + "# neighbourhood = H3Neighbourhood(regions_gdf=regions)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# embedder.fit(regions, features, joint)\n", + "# embeddings = embedder.transform(regions, features, joint)" ] }, { @@ -1067,26 +792,8 @@ "metadata": {}, "outputs": [], "source": [ - "test_data, holdout_routes, test_config = data_loader.load_h5(args.train_data_folders, test_dates, holdout_routes=model.holdout_routes, config=model.config)\n", - "test_dataset = data_loader.H5Dataset(test_data)\n", - "test_dataset.include_grid = model.include_grid\n", - "test_loader = DataLoader(\n", - " test_dataset,\n", - " collate_fn=model.collate_fn,\n", - " batch_size=model.batch_size,\n", - " shuffle=False,\n", - " drop_last=False,\n", - " num_workers=num_workers,\n", - " pin_memory=pin_memory\n", - ")\n", - "trainer = pl.Trainer(\n", - " accelerator=accelerator,\n", - " logger=False,\n", - " inference_mode=True\n", - ")\n", - "preds_and_labels = trainer.predict(model=model, dataloaders=test_loader)\n", - "preds = np.concatenate([x['preds'] for x in preds_and_labels])\n", - "labels = np.concatenate([x['labels'] for x in preds_and_labels])" + "# folium_map = plot_regions(area, colormap=[\"rgba(0,0,0,0.1)\"], tiles_style=\"CartoDB positron\")\n", + "# plot_numeric_data(regions, 0, embeddings, map=folium_map)" ] } ], diff --git a/openbustools/standardfeeds.py b/openbustools/standardfeeds.py index 02b540d..9172a29 100644 --- a/openbustools/standardfeeds.py +++ b/openbustools/standardfeeds.py @@ -25,6 +25,29 @@ GTFS_LOOKUP = dict(zip(GTFS_NAMES, GTFS_TYPES)) +def validate_realtime_data(row): + """ + Validates the realtime data for a given city. + + Args: + row (pandas.Series): The row containing the city info to be validated. + + Returns: + bool: True if there is loadable realtime data, False otherwise. + """ + uuid = row['uuid'] + provider_path = Path('data', 'other_feeds', f"{uuid}_realtime") + available_files = [x.name for x in provider_path.glob('*.pkl')] + if len(available_files) != 0: + for file in available_files: + data = pd.read_pickle(provider_path / file) + if len(data) == 0: + return False + return True + else: + return False + + def get_time_embeddings(timestamps, timezone='America/Los_Angeles'): """ Converts a list of timestamps to time embeddings. @@ -74,13 +97,13 @@ def get_gtfs_shapes_lookup(gtfs_folder): routes = pd.read_csv(f"{gtfs_folder}routes.txt", low_memory=False, dtype=GTFS_LOOKUP) trips = pd.read_csv(f"{gtfs_folder}trips.txt", low_memory=False, dtype=GTFS_LOOKUP) data = trips.merge(routes, on='route_id') - data = data[['service_id','route_id','direction_id','trip_id']].drop_duplicates().sort_values(['service_id','route_id','direction_id','trip_id']) + data = data[['service_id','route_id','trip_id']].drop_duplicates().sort_values(['service_id','route_id','trip_id']) return data def get_gtfs_shapes(gtfs_folder, epsg=None, stop_dist_filter=None): """ - Use stop locations to create unique shapes from GTFS files. + Use stop locations to create unique shapes for each trip in GTFS files. Args: gtfs_folder (str): The path to the folder containing the GTFS files. @@ -110,9 +133,9 @@ def get_gtfs_shapes(gtfs_folder, epsg=None, stop_dist_filter=None): data = data[~data['trip_id'].isin(drop_ids)] # Avoid making duplicate shapes; one per service/route/direction shapes_lookup = get_gtfs_shapes_lookup(gtfs_folder) - data = data.merge(shapes_lookup, on='trip_id').drop_duplicates(['service_id','route_id','direction_id','stop_id']).sort_values(['service_id','route_id','direction_id','stop_sequence']) - data = data.groupby(['service_id','route_id','direction_id'], as_index=False)['geometry'].apply(lambda x: LineString(x.tolist())) - data.crs = f"EPSG: {epsg}" + data = pd.merge(shapes_lookup, data, on=['trip_id'])[['service_id','route_id','trip_id','stop_id','stop_sequence','geometry']].sort_values(['service_id','route_id','trip_id','stop_sequence']).drop_duplicates() + data = data.groupby(['service_id','route_id','trip_id'], as_index=False)['geometry'].apply(lambda x: LineString(x.tolist())) + data = gpd.GeoDataFrame(data, geometry='geometry', crs=f"EPSG: {epsg}") return data @@ -129,7 +152,7 @@ def load_gtfs_files(gtfs_folder): # Read stop locations, trips/times, and route ids stop_times = pd.read_csv(Path(gtfs_folder, "stop_times.txt"), low_memory=False, dtype=GTFS_LOOKUP)[['trip_id','stop_id','arrival_time','stop_sequence']].dropna() stops = pd.read_csv(Path(gtfs_folder, "stops.txt"), low_memory=False, dtype=GTFS_LOOKUP)[['stop_id','stop_lon','stop_lat']].sort_values('stop_id').dropna() - trips = pd.read_csv(Path(gtfs_folder, "trips.txt"), low_memory=False, dtype=GTFS_LOOKUP)[['trip_id','service_id','route_id','direction_id']].sort_values(['service_id','route_id','direction_id','trip_id']).dropna() + trips = pd.read_csv(Path(gtfs_folder, "trips.txt"), low_memory=False, dtype=GTFS_LOOKUP)[['trip_id','service_id','route_id']].sort_values(['service_id','route_id','trip_id']).dropna() # Deal with schedule crossing midnight, get scheduled arrival stop_times['t_sch_sec_of_day'] = [int(x[0])*60*60 + int(x[1])*60 + int(x[2]) for x in stop_times['arrival_time'].str.split(":")] stop_times = stop_times.sort_values(['trip_id','t_sch_sec_of_day']) @@ -348,30 +371,4 @@ def combine_phone_sensors(phone_folder, timezone, chop_n=None): 'start_epoch': start_epoch, 'end_epoch': end_epoch } - return (metadata, all_sensors) - - -# def filter_gtfs_w_phone(phone_df, gtfs_df, route_short_name, gtfs_calendar): -# """Filter bus schedule using best guess for trip_ids that correspond to a phone trip.""" -# # Filter to the route number shown on the headsign -# filtered_df = gtfs_df[gtfs_df['route_short_name']==route_short_name] -# # # Filter by service ids that are active on day of week -# # weekdays = ("monday","tuesday","wednesday","thursday","friday","saturday","sunday") -# # phone_start_time = datetime.fromtimestamp(int(str(min(phone_df.time))[:10])) -# # phone_start_time = phone_start_time.hour*3600 + phone_start_time.second -# # valid_service_ids = gtfs_calendar[gtfs_calendar[weekdays[phone_start_time.weekday()]]==1].service_id -# # filtered_df = filtered_df[filtered_df['service_id'].isin(valid_service_ids)] -# # Filter to direction 0; will need to be improved -# filtered_df = filtered_df[filtered_df['direction_id']==0] -# # # Filter to trips scheduled active when the phone started recording -# # trip_times = filtered_df[['trip_id','arrival_s']].groupby('trip_id', as_index=False).agg(['min','max']) -# # trip_times = trip_times[phone_start_time > trip_times[('arrival_s','min')]] -# # trip_times = trip_times[phone_start_time < trip_times[('arrival_s','max')]] -# # filtered_df = filtered_df[filtered_df['trip_id'].isin(trip_times[('trip_id','')])] -# # filtered_df = filtered_df.sort_values(['trip_id','arrival_s']) -# # Return only one trip_id in the dataframe for plotting, regardless of total remaining -# remaining_trip_ids = np.unique(filtered_df.trip_id) -# keep_trip_id = remaining_trip_ids[0] -# print(f"Filtered down to {len(remaining_trip_ids)} possible trip_ids; returning the first one ({keep_trip_id}).") -# filtered_df = filtered_df[filtered_df['trip_id']==keep_trip_id] -# return filtered_df.copy(), remaining_trip_ids \ No newline at end of file + return (metadata, all_sensors) \ No newline at end of file diff --git a/scripts/process_data.py b/scripts/process_data.py index 67f114f..2d60f7b 100755 --- a/scripts/process_data.py +++ b/scripts/process_data.py @@ -15,41 +15,49 @@ def prepare_run(**kwargs): + logger = logging.getLogger('data_processing') + logger.setLevel(logging.DEBUG) + ch = logging.StreamHandler() + ch.setLevel(logging.DEBUG) + formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') + ch.setFormatter(formatter) + logger.addHandler(ch) + """Pre-process training data and save to sub-folder.""" print(f"PROCESSING: {kwargs['network_name']}") for day in kwargs['dates']: print(day) # Loading data and unifying column names/dtypes - try: - data = pd.read_pickle(Path(kwargs['realtime_folder'], day)) - num_pts_initial = len(data) - print(f"Loaded {num_pts_initial:_} points") - except: - logging.warning(f"File failed to load: {day}") - continue - if num_pts_initial == 0: - logging.warning(f"No points found: {day}") - continue - elif num_pts_initial < 100: - logging.warning(f"Too few points found: {day}") - continue + data = pd.read_pickle(Path(kwargs['realtime_folder'], day)) data['file'] = day[:10] data = data[kwargs['given_names']] data.columns = standardfeeds.GTFSRT_LOOKUP.keys() data['locationtime'] = data['locationtime'].astype(float) data = data.astype(standardfeeds.GTFSRT_LOOKUP) + logger.debug(f"Loaded and unified columns: {len(data):_} points") + if len(data) == 0: + continue # Sensors seem to hold old positions right at start/end of trip data = data.reset_index(drop=True) for _ in range(3): data = data.drop(data.groupby('trip_id', as_index=False).nth(0).index) data = data.drop(data.groupby('trip_id', as_index=False).nth(-1).index) + logger.debug(f"Removed trip start/end points: {len(data):_} points") + if len(data) == 0: + continue # Split full trip trajectories into smaller samples, resample data = spatial.shingle(data, min_break=2, max_break=5, min_len=3, max_len=200) + logger.debug(f"Shingled: {len(data):_} points") + if len(data) == 0: + continue # Project to local coordinate system, apply bounding box, center coords data = spatial.create_bounded_gdf(data, 'lon', 'lat', kwargs['epsg'], kwargs['coord_ref_center'], kwargs['grid_bounds'], kwargs['dem_file']) + logger.debug(f"Created bounded gdf: {len(data):_} points") + if len(data) == 0: + continue # Calculate geometry features data['calc_time_s'] = data['locationtime'] - data['locationtime'].shift(1) @@ -63,6 +71,9 @@ def prepare_run(**kwargs): data = data.drop(data.groupby('shingle_id', as_index=False, sort=False).nth(0).index) data['calc_speed_m_s'] = data['calc_dist_m'] / data['calc_time_s'] data['calc_dist_km'] = data['calc_dist_m'] / 1000.0 + logger.debug(f"Calculated geometry: {len(data):_} points") + if len(data) == 0: + continue # Filter trips with less than 2 points toss_ids = [] @@ -77,6 +88,9 @@ def prepare_run(**kwargs): # Filter the list of full shingles w/invalid points toss_ids = np.unique(toss_ids) data = data[~data['shingle_id'].isin(toss_ids)].copy() + logger.debug(f"Filtered invalid trips: {len(data):_} points") + if len(data) == 0: + continue # Calculate time features data['t'] = pd.to_datetime(data['locationtime'], unit='s', utc=True).dt.tz_convert(kwargs['timezone']) @@ -92,6 +106,9 @@ def prepare_run(**kwargs): # For calculating absolute time differences in trips (midnight crossover) data['t_sec_of_day'] = data['t'] - datetime.datetime(min(data['t_year']), min(data['t_month']), min(data['t_day']), 0, tzinfo=ZoneInfo(kwargs['timezone'])) data['t_sec_of_day'] = data['t_sec_of_day'].dt.total_seconds() + logger.debug(f"Calculated time features: {len(data):_} points") + if len(data) == 0: + continue # Get GTFS features best_static = standardfeeds.latest_available_static(day[:10], kwargs['static_folder']) @@ -108,12 +125,16 @@ def prepare_run(**kwargs): data = data.merge(trips, on='trip_id', how='left') data['calc_stop_dist_km'] = data['calc_stop_dist_m'] / 1000.0 else: - logging.warning(f"No data after joining static feed: {day}") + logger.warning(f"No data after joining static feed: {day}") data['stop_sequence'] = 0 data['t_sch_sec_of_day'] = 0 data['calc_stop_dist_m'] = 0 data['calc_stop_dist_km'] = 0 data['route_id'] = 'NotFound' + logger.debug(f"Joined static feed: {len(data):_} points") + if len(data) == 0: + continue + # Passed stops data['pass_stops_n'] = data.groupby('shingle_id')['stop_sequence'].diff().fillna(0) # Scheduled time @@ -138,8 +159,8 @@ def prepare_run(**kwargs): data_grid.build_cell_lookup(data[['locationtime','x','y','calc_speed_m_s','calc_bear_d']].copy()) # Save processed data - num_pts_final = len(data) - print(f"Kept {np.round(num_pts_final/num_pts_initial, 3)*100}% of original points") + logger.debug(f"Final data: {len(data):_} points") + # Full geodataframe processed_path = Path(kwargs['realtime_folder'], "processed") processed_path.mkdir(parents=True, exist_ok=True) @@ -161,7 +182,7 @@ def prepare_run(**kwargs): g.create_dataset('feats_n', data=data_n) g.create_dataset('feats_c', data=data_c) g.create_dataset('feats_g', data=data_g) - print(f"PROCESSING COMPLETED: {kwargs['network_name']}") + print(f"PROCESSING COMPLETED: {kwargs['network_name']}\n") if __name__=="__main__": @@ -210,19 +231,18 @@ def prepare_run(**kwargs): # given_names=['trip_id','file','locationtime','lat','lon','vehicle_id'], # ) cleaned_sources = pd.read_csv(Path('data', 'cleaned_sources.csv')) - cleaned_sources = cleaned_sources.iloc[35:] - cleaned_sources = cleaned_sources[cleaned_sources['provider']=='York Region Transit'] for i, row in cleaned_sources.iterrows(): - print(row['provider']) - prepare_run( - network_name=row['uuid'], - dates=['2024_01_01.pkl'], - static_folder=f"./data/other_feeds/{row['uuid']}_static/", - realtime_folder=f"./data/other_feeds/{row['uuid']}_realtime/", - timezone=row['tz_str'], - epsg=row['epsg_code'], - grid_bounds=[row['min_lon'], row['min_lat'], row['max_lon'], row['max_lat']], - coord_ref_center=[np.mean([row['min_lon'], row['max_lon']]), np.mean([row['min_lat'], row['max_lat']])], - dem_file=[x for x in Path('data', 'other_feeds', f"{row['uuid']}_spatial").glob(f"*_{row['epsg_code']}.tif")][0], - given_names=['trip_id','file','locationtime','lat','lon','vehicle_id'], - ) \ No newline at end of file + if standardfeeds.validate_realtime_data(row): + print(row['provider']) + prepare_run( + network_name=row['uuid'], + dates=['2024_01_03.pkl'], + static_folder=Path('data', 'other_feeds', f"{row['uuid']}_static"), + realtime_folder=Path('data', 'other_feeds', f"{row['uuid']}_realtime"), + timezone=row['tz_str'], + epsg=row['epsg_code'], + grid_bounds=[row['min_lon'], row['min_lat'], row['max_lon'], row['max_lat']], + coord_ref_center=[np.mean([row['min_lon'], row['max_lon']]), np.mean([row['min_lat'], row['max_lat']])], + dem_file=[x for x in Path('data', 'other_feeds', f"{row['uuid']}_spatial").glob(f"*_{row['epsg_code']}.tif")][0], + given_names=['trip_id','file','locationtime','lat','lon','vehicle_id'], + ) \ No newline at end of file