{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Analysis of Gas Prices in the GTA\n",
"We obtain gas station price data for the Greater Toronto Area with data from yellowpages.ca and gasbuddy.com and cluster them based on location and price per litre of regular gas. We also test the price effects of having another gas station or Starbucks location within a 100m radius of a gas station."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from bs4 import BeautifulSoup\n",
"import urllib\n",
"import requests\n",
"import math\n",
"import re\n",
"import pandas as pd\n",
"import numpy as np\n",
"import datetime\n",
"import sys\n",
"import time\n",
"import json\n",
"import ast\n",
"from geopy.geocoders import GoogleV3\n",
"import seaborn as sns\n",
"from time import sleep\n",
"geolocator = GoogleV3()\n",
"%matplotlib inline\n",
"today = datetime.datetime.now()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Get the names and addressess of gas stations listed in Yellow Pages in the Toronto region\n",
"\n",
"def get_addresses (page_num):\n",
" with open('gas_dict.json', 'r') as infile:\n",
" gas_dict = json.load(infile)\n",
" \n",
" url = 'https://www.yellowpages.ca/search/si/'+str(page_num)+'/Gas%20Stations/Toronto+ON'\n",
" try:\n",
" html = urllib.request.urlopen(url).read()\n",
" except:\n",
" return page_num\n",
" soup = BeautifulSoup(html, 'html.parser')\n",
" \n",
" content = soup.findAll('div',{'class':'listing__content__wrap'})\n",
"\n",
" for n in range(len(content)):\n",
" try:\n",
" name = content[n].findAll('a')[0].get('title').split('-')[0]\n",
" except:\n",
" name = '#NO NAME'\n",
" \n",
" try:\n",
" addr_ele = content[n].findAll('span',{'class':'listing__address--full'})[0].findAll('span')\n",
" addr = ','.join([addr_ele[i].string for i,e in enumerate(addr_ele)])\n",
" except:\n",
" addr = '#NO ADDRESS#'\n",
" \n",
" gas_dict[addr] = {'name':name}\n",
" \n",
" with open('gas_dict.json', 'w') as outfile:\n",
" json.dump(gas_dict,outfile)\n",
" \n",
" return page_num + 1\n",
"\n",
"def address_main():\n",
" # If file doesn't exist, create it\n",
" try:\n",
" with open('gas_dict.json', 'r') as infile:\n",
" gas_dict = json.load(infile)\n",
" except:\n",
" with open('gas_dict.json', 'w') as outfile:\n",
" gas_dict={}\n",
" json.dump(gas_dict,outfile)\n",
" \n",
" page_num = 1\n",
" while page_num < 22:\n",
" print('Fetching page ',page_num)\n",
" page_num = get_addresses(page_num)\n",
" sleep(np.random.random(1)*2+2)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Get rid of incomplete or empty addresses\n",
"def clean_addresses():\n",
" with open('gas_dict.json', 'r') as infile:\n",
" gas_dict = json.load(infile)\n",
" \n",
" gas_dict.pop('#NO ADDRESS#',None)\n",
" gas_dict.pop('ON',None)\n",
" gas_dict.pop('12001 Hwy 400,Maple,ON,L7B 1A8',None) # This is actually both redundant and the wrong city ...\n",
" \n",
" with open('gas_dict.json', 'w') as outfile:\n",
" json.dump(gas_dict,outfile)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# for ontariogasprices.com scraping, we need a list of postal codes to search with\n",
"# the site focuses on finding the lowest price in an area and does not prioritize towards specific addresses\n",
"# but this works to our advantage since we can potentially find stations not included in the yellowpages data\n",
"# gas types are\n",
"# A: regular, B: mid, C: premium, D: diesel\n",
"\n",
"def get_area_prices(gas_type, postal_code):\n",
" with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json', 'r') as infile:\n",
" gas_stations = json.load(infile)\n",
" \n",
" url = 'http://www.ontariogasprices.com/GasPriceSearch.aspx?fuel='+gas_type+'&qsrch='+postal_code[0]+'%20'+postal_code[1]\n",
" print ('Fetching from ', url) \n",
" html = requests.get(url).text\n",
" soup = BeautifulSoup(html, 'html.parser') \n",
"\n",
" if len(soup.findAll('tbody')) == 0:\n",
" pass\n",
" else:\n",
" data_rows = soup.findAll('tbody')[0].findAll('tr')\n",
"\n",
" for i in range(len(data_rows)):\n",
" try:\n",
" name = data_rows[i].findAll('td')[0].find('a').string\n",
" st_addr = data_rows[i].findAll('td')[0].find('dd').string.split('&')[0].split('near')[0]\n",
" city = data_rows[i].findAll('td')[2].find('a').string\n",
" addr = st_addr + ', '+ city + ', Ontario'\n",
" price = float(data_rows[i].find('a').string)\n",
" if addr not in gas_stations:\n",
" gas_stations[addr] = {'name':name}\n",
" else:\n",
" pass\n",
" gas_stations[addr]['type_'+gas_type] = price\n",
" except:\n",
" # The \"other stations in the area\" row can just be ignored\n",
" pass\n",
" \n",
" with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json', 'w') as outfile:\n",
" json.dump(gas_stations,outfile)\n",
" \n",
"def price_main():\n",
" # first get postal codes\n",
" with open('gas_dict.json', 'r') as infile:\n",
" gas_dict = json.load(infile)\n",
" \n",
" try:\n",
" with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json', 'r') as infile:\n",
" gas_stations = json.load(infile)\n",
" except:\n",
" with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json', 'w') as outfile:\n",
" gas_stations={}\n",
" json.dump(gas_stations,outfile)\n",
" \n",
" postal_list=[]\n",
" for k,v in gas_dict.items():\n",
" postal_list.append(k.split(',')[-1].split())\n",
" \n",
" # If no postal code, just ignore for now...\n",
" postal_list = [postal_list[i] for i,p in enumerate(postal_list) if len(p) == 2]\n",
"\n",
" break_point = 0\n",
" while (break_point < len(postal_list)):\n",
" postal_code = postal_list[break_point]\n",
" try:\n",
" for gas_type in list('ABCD'):\n",
" get_area_prices(gas_type, postal_code)\n",
" sleep(np.random.random(1)*2+0.5) \n",
" break_point += 1\n",
" except:\n",
" # Try again\n",
" print('Trying again ...')\n",
" sleep(np.random.random(1)*2+2)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Get not only the coordinates but also the Google Maps address, which will likely clean up any messes left over\n",
"# from the ontariogasprices.com splits\n",
"\n",
"def get_coord(addr):\n",
" # splitting Toronto addresses into \"Toronto - [area]\" messes up google maps, so get rid of that when querying\n",
" if 'Toronto' in addr.split(',')[1]:\n",
" addr = addr.split(',')[0] + ', Toronto, '+addr.split(',')[2]\n",
" print ('Getting' , addr)\n",
" geoloc = geolocator.geocode(addr)\n",
" sleep(np.random.random(1)*2+2)\n",
" print ('Got ', geoloc, '\\n')\n",
" return geoloc\n",
" \n",
"def coords_main():\n",
" # Start at the beginning of gas_dict\n",
" break_point = 0\n",
"\n",
" # open gas_stations, and get the list of all the addresses\n",
" with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json', 'r') as infile:\n",
" gas_dict = json.load(infile)\n",
" addr_list = [k for k,v in gas_dict.items()]\n",
" \n",
" # do while there are still addresses to process\n",
" while(break_point < len(addr_list)):\n",
" # run collect_coords starting with the current breakpoint, and gets the\n",
" # address that it's stuck on\n",
" addr = addr_list[break_point]\n",
" try: \n",
" geoloc = get_coord(addr)\n",
" if geoloc == None:\n",
" # GoogleV3 can't get a hold of it. Need to look at it later. Most likely it's a highway. \n",
" break_point += 1\n",
" sleep(np.random.random(1)*2+2) \n",
" else:\n",
" gas_dict[addr]['address'] = geoloc.address\n",
" gas_dict[addr]['longitude'] = geoloc.longitude\n",
" gas_dict[addr]['latitude'] = geoloc.latitude\n",
" # If it works, get to the next one. Once break_point gets to the length of the list, the conditions of\n",
" # the loop will no longer be satisfied, and the loop ends.\n",
" break_point +=1\n",
" with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json','w') as outfile:\n",
" json.dump(gas_dict,outfile) \n",
" sleep(np.random.random(1)*2+2) \n",
" except:\n",
" # If it doesn't work, don't increase the break_point and try again\n",
" print ('Could not get ', addr)\n",
" sleep(np.random.random(1)*2+2)\n",
"\n",
" with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json','w') as outfile:\n",
" json.dump(gas_dict,outfile) \n",
" #return gas_dict"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def xy_distance(xlat,xlong,ylat,ylong):\n",
" # The radius of the planet is approx. 6371.01 km\n",
" xlat = math.radians(xlat)\n",
" xlong = math.radians(xlong)\n",
" ylat = math.radians(ylat)\n",
" ylong = math.radians(ylong)\n",
" dist = 6371.01 * math.acos(math.sin(xlat)*math.sin(ylat) + math.cos(xlat)*math.cos(ylat)*math.cos(xlong - ylong))\n",
" return dist\n",
"\n",
"def xy_closest(coords, x):\n",
" return (min(coords, key=lambda y: xy_distance(x[1],x[0],y[1],y[0])))"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def enter_address(df, addr, max_dist):\n",
" standard_addr = get_coord(addr)\n",
" within_range = pd.DataFrame()\n",
" for i in range(len(df)):\n",
" dist = xy_distance(standard_addr.latitude,standard_addr.longitude,df.iloc[i]['latitude'],df.iloc[i]['longitude'])\n",
" if dist <= max_dist:\n",
" within_range = pd.concat((within_range,df.iloc[i]),axis=1)\n",
"\n",
" within_range = within_range.transpose().sort_values('type_A')\n",
" \n",
" # Hide the entered address in the dataframe too so we can show it on the map\n",
" print ('Gas stations within range: ')\n",
" print (within_range[['name','type_A']])\n",
" folmap(within_range,15,custom=1,extralat=standard_addr.latitude,extralong=standard_addr.longitude)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def folmap (df,zoom=10,custom=None,extralat=None,extralong=None):\n",
" k=7\n",
" colorscale = branca.colormap.linear.YlGnBu.scale(0, 60)\n",
" cluster_colors = ['red', 'green', 'blue', 'yellow', 'cadetblue', 'purple', 'black']\n",
"\n",
" if custom == 1:\n",
" m = folium.Map(\n",
" location=[extralat, extralong],\n",
" zoom_start=zoom\n",
" )\n",
" else:\n",
" # center the map on the median coordinates\n",
" m = folium.Map(\n",
" location=[np.median(df.latitude), np.median(df.longitude)],\n",
" zoom_start=zoom\n",
" )\n",
"\n",
" marker_cluster = MarkerCluster().add_to(m)\n",
"\n",
" for j in range(k):\n",
" for r in df[:-1].iterrows():\n",
" if r[1]['label'] == j:\n",
" pop = [\"Price: \",str(r[1]['type_A'])]\n",
" folium.Marker(\n",
" location = [r[1]['latitude'],r[1]['longitude']],\n",
" popup=''.join(pop),\n",
" icon=folium.Icon(color=cluster_colors[j], icon='ok-sign'),\n",
" ).add_to(marker_cluster)\n",
" \n",
" if custom == 1:\n",
" folium.Marker(location = [extralat,extralong],\n",
" popup = 'You are here',\n",
" icon=folium.Icon(color='pink', icon='cloud')).add_to(m)\n",
" \n",
" display(m)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Run the 4 main functions to get a json file of addresses and prices, and a list of coordinates."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"#address_main()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"#clean_addresses()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"#price_main()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"#coords_main()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"712"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"with open('gas_stations_2017-12-13.json', 'r') as infile:\n",
" gas_stations = json.load(infile)\n",
"len(gas_stations)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Clean up\n",
"\n",
"For some reason, there's a station near Vaughn that only registers in Google Maps if we input it as King City. There are also a few duplicate stations that are listed under different names or address variants, so we get rid of those."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'address': '12001 ON-400 #1S1, King City, ON L7B 1A8, Canada',\n",
" 'latitude': 43.8950319,\n",
" 'longitude': -79.5574951,\n",
" 'name': 'Canadian Tire',\n",
" 'type_A': 125.3,\n",
" 'type_D': 119.7}"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Special treatment for 12001 hwy 400 nb, since gasbuddy has it as in Vaughan but Google recognizes it as being in King City\n",
"gas_stations['12001 Hwy 400 NB , Vaughan, Ontario']"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"geoloc = geolocator.geocode('12001 Hwy 400 NB , King City, Ontario')"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"gas_stations['12001 Hwy 400 NB , Vaughan, Ontario']['address'] = geoloc.address\n",
"gas_stations['12001 Hwy 400 NB , Vaughan, Ontario']['longitude'] = geoloc.longitude\n",
"gas_stations['12001 Hwy 400 NB , Vaughan, Ontario']['latitude'] =geoloc.latitude"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"gas_df = pd.DataFrame.from_dict(gas_stations, orient='index')"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" type_B | \n",
" type_C | \n",
" type_D | \n",
" address | \n",
" longitude | \n",
" latitude | \n",
"
\n",
" \n",
" \n",
" \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario | \n",
" Ultramar | \n",
" 115.6 | \n",
" 129.6 | \n",
" 132.6 | \n",
" 115.9 | \n",
" 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... | \n",
" -79.282378 | \n",
" 43.976346 | \n",
"
\n",
" \n",
" 1610 Keele St , Toronto - Central, Ontario | \n",
" Shell | \n",
" 122.9 | \n",
" 136.9 | \n",
" 144.9 | \n",
" NaN | \n",
" 1610 Keele St, York, ON M6M 3V9, Canada | \n",
" -79.471988 | \n",
" 43.682320 | \n",
"
\n",
" \n",
" 6897 Finch Ave W , Toronto - West, Ontario | \n",
" Esso | \n",
" 121.9 | \n",
" 132.9 | \n",
" 142.9 | \n",
" 118.9 | \n",
" 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada | \n",
" -79.617364 | \n",
" 43.734632 | \n",
"
\n",
" \n",
" 1 Harwood Ave S , Ajax, Ontario | \n",
" Pioneer | \n",
" 122.6 | \n",
" 123.4 | \n",
" 130.9 | \n",
" 118.6 | \n",
" 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada | \n",
" -79.025100 | \n",
" 43.861027 | \n",
"
\n",
" \n",
" 1 Thornhill Woods Dr , Vaughan, Ontario | \n",
" Esso | \n",
" 122.9 | \n",
" NaN | \n",
" NaN | \n",
" 118.9 | \n",
" 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... | \n",
" -79.463886 | \n",
" 43.826660 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A type_B \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 129.6 \n",
" 1610 Keele St , Toronto - Central, Ontario Shell 122.9 136.9 \n",
" 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 132.9 \n",
"1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 123.4 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 NaN \n",
"\n",
" type_C type_D \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 132.6 115.9 \n",
" 1610 Keele St , Toronto - Central, Ontario 144.9 NaN \n",
" 6897 Finch Ave W , Toronto - West, Ontario 142.9 118.9 \n",
"1 Harwood Ave S , Ajax, Ontario 130.9 118.6 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario NaN 118.9 \n",
"\n",
" address \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... \n",
" 1610 Keele St , Toronto - Central, Ontario 1610 Keele St, York, ON M6M 3V9, Canada \n",
" 6897 Finch Ave W , Toronto - West, Ontario 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada \n",
"1 Harwood Ave S , Ajax, Ontario 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... \n",
"\n",
" longitude latitude \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n",
" 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n",
"1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 "
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gas_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" type_B | \n",
" type_C | \n",
" type_D | \n",
" address | \n",
" longitude | \n",
" latitude | \n",
"
\n",
" \n",
" \n",
" \n",
" 995 Pape Ave , East York, Ontario | \n",
" Cocomile Service Centre | \n",
" 115.9 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" 995 Pape Ave, East York, ON M4K 3V8, Canada | \n",
" -79.347831 | \n",
" 43.687761 | \n",
"
\n",
" \n",
" 995 Pape Ave , Toronto - Central, Ontario | \n",
" National Gas & Variety | \n",
" 115.9 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" 995 Pape Ave, East York, ON M4K 3V8, Canada | \n",
" -79.347831 | \n",
" 43.687761 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A \\\n",
"995 Pape Ave , East York, Ontario Cocomile Service Centre 115.9 \n",
"995 Pape Ave , Toronto - Central, Ontario National Gas & Variety 115.9 \n",
"\n",
" type_B type_C type_D \\\n",
"995 Pape Ave , East York, Ontario NaN NaN NaN \n",
"995 Pape Ave , Toronto - Central, Ontario NaN NaN NaN \n",
"\n",
" address \\\n",
"995 Pape Ave , East York, Ontario 995 Pape Ave, East York, ON M4K 3V8, Canada \n",
"995 Pape Ave , Toronto - Central, Ontario 995 Pape Ave, East York, ON M4K 3V8, Canada \n",
"\n",
" longitude latitude \n",
"995 Pape Ave , East York, Ontario -79.347831 43.687761 \n",
"995 Pape Ave , Toronto - Central, Ontario -79.347831 43.687761 "
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gas_df[gas_df['longitude']==-79.347831]"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" type_B | \n",
" type_C | \n",
" type_D | \n",
" address | \n",
" longitude | \n",
" latitude | \n",
"
\n",
" \n",
" \n",
" \n",
" 490 Great Lakes Dr , Brampton, ON, Ontario | \n",
" Shell | \n",
" 118.9 | \n",
" 132.9 | \n",
" 140.9 | \n",
" 118.9 | \n",
" 490 Great Lakes Dr, Brampton, ON L6R, Canada | \n",
" -79.773645 | \n",
" 43.740888 | \n",
"
\n",
" \n",
" 490 Great Lakes Dr , Brampton, Ontario | \n",
" Shell | \n",
" 118.9 | \n",
" 132.9 | \n",
" 140.9 | \n",
" 118.9 | \n",
" 490 Great Lakes Dr, Brampton, ON L6R, Canada | \n",
" -79.773645 | \n",
" 43.740888 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A type_B type_C \\\n",
"490 Great Lakes Dr , Brampton, ON, Ontario Shell 118.9 132.9 140.9 \n",
"490 Great Lakes Dr , Brampton, Ontario Shell 118.9 132.9 140.9 \n",
"\n",
" type_D \\\n",
"490 Great Lakes Dr , Brampton, ON, Ontario 118.9 \n",
"490 Great Lakes Dr , Brampton, Ontario 118.9 \n",
"\n",
" address \\\n",
"490 Great Lakes Dr , Brampton, ON, Ontario 490 Great Lakes Dr, Brampton, ON L6R, Canada \n",
"490 Great Lakes Dr , Brampton, Ontario 490 Great Lakes Dr, Brampton, ON L6R, Canada \n",
"\n",
" longitude latitude \n",
"490 Great Lakes Dr , Brampton, ON, Ontario -79.773645 43.740888 \n",
"490 Great Lakes Dr , Brampton, Ontario -79.773645 43.740888 "
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gas_df[gas_df['longitude']==-79.7736453]"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Get rid of redundancies\n",
"gas_df.drop_duplicates(subset=['longitude','latitude'], inplace=True)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" type_B | \n",
" type_C | \n",
" type_D | \n",
" address | \n",
" longitude | \n",
" latitude | \n",
"
\n",
" \n",
" \n",
" \n",
" 995 Pape Ave , East York, Ontario | \n",
" Cocomile Service Centre | \n",
" 115.9 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" 995 Pape Ave, East York, ON M4K 3V8, Canada | \n",
" -79.347831 | \n",
" 43.687761 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A type_B \\\n",
"995 Pape Ave , East York, Ontario Cocomile Service Centre 115.9 NaN \n",
"\n",
" type_C type_D \\\n",
"995 Pape Ave , East York, Ontario NaN NaN \n",
"\n",
" address \\\n",
"995 Pape Ave , East York, Ontario 995 Pape Ave, East York, ON M4K 3V8, Canada \n",
"\n",
" longitude latitude \n",
"995 Pape Ave , East York, Ontario -79.347831 43.687761 "
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gas_df[gas_df['longitude']==-79.347831]"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" type_B | \n",
" type_C | \n",
" type_D | \n",
" address | \n",
" longitude | \n",
" latitude | \n",
"
\n",
" \n",
" \n",
" \n",
" 490 Great Lakes Dr , Brampton, ON, Ontario | \n",
" Shell | \n",
" 118.9 | \n",
" 132.9 | \n",
" 140.9 | \n",
" 118.9 | \n",
" 490 Great Lakes Dr, Brampton, ON L6R, Canada | \n",
" -79.773645 | \n",
" 43.740888 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A type_B type_C \\\n",
"490 Great Lakes Dr , Brampton, ON, Ontario Shell 118.9 132.9 140.9 \n",
"\n",
" type_D \\\n",
"490 Great Lakes Dr , Brampton, ON, Ontario 118.9 \n",
"\n",
" address \\\n",
"490 Great Lakes Dr , Brampton, ON, Ontario 490 Great Lakes Dr, Brampton, ON L6R, Canada \n",
"\n",
" longitude latitude \n",
"490 Great Lakes Dr , Brampton, ON, Ontario -79.773645 43.740888 "
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gas_df[gas_df['longitude']==-79.7736453]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Mean prices by brand\n",
"The top prices may be unreliable, since they are not major chains or might be out of the way, and so user data on their prices is scarce. Not surprisingly, the cheapest brand is Costco. It's a members only service, and in some sense doesn't really count. It's prices are so low that in the K Means clustering below, the price discrepancy largely overpowers any distance effects."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" count | \n",
" mean | \n",
" std | \n",
" min | \n",
" 25% | \n",
" 50% | \n",
" 75% | \n",
" max | \n",
"
\n",
" \n",
" name | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" Speedy Gas | \n",
" 1.0 | \n",
" 125.900000 | \n",
" NaN | \n",
" 125.9 | \n",
" 125.9 | \n",
" 125.9 | \n",
" 125.9 | \n",
" 125.9 | \n",
"
\n",
" \n",
" Good Guys Gas Bar | \n",
" 1.0 | \n",
" 122.900000 | \n",
" NaN | \n",
" 122.9 | \n",
" 122.9 | \n",
" 122.9 | \n",
" 122.9 | \n",
" 122.9 | \n",
"
\n",
" \n",
" Falcon | \n",
" 1.0 | \n",
" 122.900000 | \n",
" NaN | \n",
" 122.9 | \n",
" 122.9 | \n",
" 122.9 | \n",
" 122.9 | \n",
" 122.9 | \n",
"
\n",
" \n",
" Superstore | \n",
" 1.0 | \n",
" 122.600000 | \n",
" NaN | \n",
" 122.6 | \n",
" 122.6 | \n",
" 122.6 | \n",
" 122.6 | \n",
" 122.6 | \n",
"
\n",
" \n",
" Fas Gas | \n",
" 2.0 | \n",
" 121.900000 | \n",
" 0.000000 | \n",
" 121.9 | \n",
" 121.9 | \n",
" 121.9 | \n",
" 121.9 | \n",
" 121.9 | \n",
"
\n",
" \n",
" Petro-Canada | \n",
" 188.0 | \n",
" 121.862766 | \n",
" 2.145533 | \n",
" 113.6 | \n",
" 122.7 | \n",
" 122.9 | \n",
" 122.9 | \n",
" 123.3 | \n",
"
\n",
" \n",
" Esso | \n",
" 173.0 | \n",
" 121.712717 | \n",
" 2.215133 | \n",
" 113.2 | \n",
" 121.3 | \n",
" 122.9 | \n",
" 122.9 | \n",
" 123.9 | \n",
"
\n",
" \n",
" Shell | \n",
" 117.0 | \n",
" 121.341026 | \n",
" 2.771970 | \n",
" 112.4 | \n",
" 119.9 | \n",
" 122.9 | \n",
" 122.9 | \n",
" 123.9 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" count mean std min 25% 50% 75% \\\n",
"name \n",
"Speedy Gas 1.0 125.900000 NaN 125.9 125.9 125.9 125.9 \n",
"Good Guys Gas Bar 1.0 122.900000 NaN 122.9 122.9 122.9 122.9 \n",
"Falcon 1.0 122.900000 NaN 122.9 122.9 122.9 122.9 \n",
"Superstore 1.0 122.600000 NaN 122.6 122.6 122.6 122.6 \n",
"Fas Gas 2.0 121.900000 0.000000 121.9 121.9 121.9 121.9 \n",
"Petro-Canada 188.0 121.862766 2.145533 113.6 122.7 122.9 122.9 \n",
"Esso 173.0 121.712717 2.215133 113.2 121.3 122.9 122.9 \n",
"Shell 117.0 121.341026 2.771970 112.4 119.9 122.9 122.9 \n",
"\n",
" max \n",
"name \n",
"Speedy Gas 125.9 \n",
"Good Guys Gas Bar 122.9 \n",
"Falcon 122.9 \n",
"Superstore 122.6 \n",
"Fas Gas 121.9 \n",
"Petro-Canada 123.3 \n",
"Esso 123.9 \n",
"Shell 123.9 "
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gas_df.groupby('name')['type_A'].describe().sort_values(by='mean', ascending=False)[:8]"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" count | \n",
" mean | \n",
" std | \n",
" min | \n",
" 25% | \n",
" 50% | \n",
" 75% | \n",
" max | \n",
"
\n",
" \n",
" name | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" Costco | \n",
" 7.0 | \n",
" 111.471429 | \n",
" 1.511858 | \n",
" 110.9 | \n",
" 110.90 | \n",
" 110.9 | \n",
" 110.9 | \n",
" 114.9 | \n",
"
\n",
" \n",
" Danforth Gas & Wash | \n",
" 1.0 | \n",
" 112.900000 | \n",
" NaN | \n",
" 112.9 | \n",
" 112.90 | \n",
" 112.9 | \n",
" 112.9 | \n",
" 112.9 | \n",
"
\n",
" \n",
" Star Gas | \n",
" 1.0 | \n",
" 113.000000 | \n",
" NaN | \n",
" 113.0 | \n",
" 113.00 | \n",
" 113.0 | \n",
" 113.0 | \n",
" 113.0 | \n",
"
\n",
" \n",
" Mobil | \n",
" 1.0 | \n",
" 113.900000 | \n",
" NaN | \n",
" 113.9 | \n",
" 113.90 | \n",
" 113.9 | \n",
" 113.9 | \n",
" 113.9 | \n",
"
\n",
" \n",
" Cocomile Service Centre | \n",
" 1.0 | \n",
" 115.900000 | \n",
" NaN | \n",
" 115.9 | \n",
" 115.90 | \n",
" 115.9 | \n",
" 115.9 | \n",
" 115.9 | \n",
"
\n",
" \n",
" OP | \n",
" 1.0 | \n",
" 116.500000 | \n",
" NaN | \n",
" 116.5 | \n",
" 116.50 | \n",
" 116.5 | \n",
" 116.5 | \n",
" 116.5 | \n",
"
\n",
" \n",
" Santak Pump | \n",
" 1.0 | \n",
" 116.500000 | \n",
" NaN | \n",
" 116.5 | \n",
" 116.50 | \n",
" 116.5 | \n",
" 116.5 | \n",
" 116.5 | \n",
"
\n",
" \n",
" Race Trac | \n",
" 3.0 | \n",
" 116.666667 | \n",
" 2.926317 | \n",
" 114.2 | \n",
" 115.05 | \n",
" 115.9 | \n",
" 117.9 | \n",
" 119.9 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" count mean std min 25% 50% \\\n",
"name \n",
"Costco 7.0 111.471429 1.511858 110.9 110.90 110.9 \n",
"Danforth Gas & Wash 1.0 112.900000 NaN 112.9 112.90 112.9 \n",
"Star Gas 1.0 113.000000 NaN 113.0 113.00 113.0 \n",
"Mobil 1.0 113.900000 NaN 113.9 113.90 113.9 \n",
"Cocomile Service Centre 1.0 115.900000 NaN 115.9 115.90 115.9 \n",
"OP 1.0 116.500000 NaN 116.5 116.50 116.5 \n",
"Santak Pump 1.0 116.500000 NaN 116.5 116.50 116.5 \n",
"Race Trac 3.0 116.666667 2.926317 114.2 115.05 115.9 \n",
"\n",
" 75% max \n",
"name \n",
"Costco 110.9 114.9 \n",
"Danforth Gas & Wash 112.9 112.9 \n",
"Star Gas 113.0 113.0 \n",
"Mobil 113.9 113.9 \n",
"Cocomile Service Centre 115.9 115.9 \n",
"OP 116.5 116.5 \n",
"Santak Pump 116.5 116.5 \n",
"Race Trac 117.9 119.9 "
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gas_df.groupby('name')['type_A'].describe().sort_values(by='mean', ascending=True)[:8]"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" type_B | \n",
" type_C | \n",
" type_D | \n",
" address | \n",
" longitude | \n",
" latitude | \n",
"
\n",
" \n",
" \n",
" \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario | \n",
" Ultramar | \n",
" 115.6 | \n",
" 129.6 | \n",
" 132.6 | \n",
" 115.9 | \n",
" 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... | \n",
" -79.282378 | \n",
" 43.976346 | \n",
"
\n",
" \n",
" 1610 Keele St , Toronto - Central, Ontario | \n",
" Shell | \n",
" 122.9 | \n",
" 136.9 | \n",
" 144.9 | \n",
" NaN | \n",
" 1610 Keele St, York, ON M6M 3V9, Canada | \n",
" -79.471988 | \n",
" 43.682320 | \n",
"
\n",
" \n",
" 6897 Finch Ave W , Toronto - West, Ontario | \n",
" Esso | \n",
" 121.9 | \n",
" 132.9 | \n",
" 142.9 | \n",
" 118.9 | \n",
" 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada | \n",
" -79.617364 | \n",
" 43.734632 | \n",
"
\n",
" \n",
" 1 Harwood Ave S , Ajax, Ontario | \n",
" Pioneer | \n",
" 122.6 | \n",
" 123.4 | \n",
" 130.9 | \n",
" 118.6 | \n",
" 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada | \n",
" -79.025100 | \n",
" 43.861027 | \n",
"
\n",
" \n",
" 1 Thornhill Woods Dr , Vaughan, Ontario | \n",
" Esso | \n",
" 122.9 | \n",
" NaN | \n",
" NaN | \n",
" 118.9 | \n",
" 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... | \n",
" -79.463886 | \n",
" 43.826660 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A type_B \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 129.6 \n",
" 1610 Keele St , Toronto - Central, Ontario Shell 122.9 136.9 \n",
" 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 132.9 \n",
"1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 123.4 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 NaN \n",
"\n",
" type_C type_D \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 132.6 115.9 \n",
" 1610 Keele St , Toronto - Central, Ontario 144.9 NaN \n",
" 6897 Finch Ave W , Toronto - West, Ontario 142.9 118.9 \n",
"1 Harwood Ave S , Ajax, Ontario 130.9 118.6 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario NaN 118.9 \n",
"\n",
" address \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... \n",
" 1610 Keele St , Toronto - Central, Ontario 1610 Keele St, York, ON M6M 3V9, Canada \n",
" 6897 Finch Ave W , Toronto - West, Ontario 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada \n",
"1 Harwood Ave S , Ajax, Ontario 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... \n",
"\n",
" longitude latitude \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n",
" 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n",
"1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 "
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gas_df.iloc[:5]"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"gas_df.to_csv('gas_stations_2017-12-13.csv')"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"with open('gas_stations_2017-12-13.json','w') as outfile:\n",
" json.dump(gas_stations,outfile) "
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"coords_df = gas_df[['longitude','latitude']].copy()"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" longitude | \n",
" latitude | \n",
"
\n",
" \n",
" \n",
" \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario | \n",
" -79.282378 | \n",
" 43.976346 | \n",
"
\n",
" \n",
" 1610 Keele St , Toronto - Central, Ontario | \n",
" -79.471988 | \n",
" 43.682320 | \n",
"
\n",
" \n",
" 6897 Finch Ave W , Toronto - West, Ontario | \n",
" -79.617364 | \n",
" 43.734632 | \n",
"
\n",
" \n",
" 1 Harwood Ave S , Ajax, Ontario | \n",
" -79.025100 | \n",
" 43.861027 | \n",
"
\n",
" \n",
" 1 Thornhill Woods Dr , Vaughan, Ontario | \n",
" -79.463886 | \n",
" 43.826660 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" longitude latitude\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346\n",
" 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320\n",
" 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632\n",
"1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027\n",
"1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"coords_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"coords_df['long_rad'] = coords_df['longitude'].apply(lambda l:math.radians(l))\n",
"coords_df['lat_rad'] = coords_df['latitude'].apply(lambda l:math.radians(l))"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" longitude | \n",
" latitude | \n",
" long_rad | \n",
" lat_rad | \n",
"
\n",
" \n",
" \n",
" \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario | \n",
" -79.282378 | \n",
" 43.976346 | \n",
" -1.383739 | \n",
" 0.767532 | \n",
"
\n",
" \n",
" 1610 Keele St , Toronto - Central, Ontario | \n",
" -79.471988 | \n",
" 43.682320 | \n",
" -1.387048 | \n",
" 0.762400 | \n",
"
\n",
" \n",
" 6897 Finch Ave W , Toronto - West, Ontario | \n",
" -79.617364 | \n",
" 43.734632 | \n",
" -1.389585 | \n",
" 0.763313 | \n",
"
\n",
" \n",
" 1 Harwood Ave S , Ajax, Ontario | \n",
" -79.025100 | \n",
" 43.861027 | \n",
" -1.379248 | \n",
" 0.765519 | \n",
"
\n",
" \n",
" 1 Thornhill Woods Dr , Vaughan, Ontario | \n",
" -79.463886 | \n",
" 43.826660 | \n",
" -1.386906 | \n",
" 0.764920 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" longitude latitude \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n",
" 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n",
"1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n",
"\n",
" long_rad lat_rad \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -1.383739 0.767532 \n",
" 1610 Keele St , Toronto - Central, Ontario -1.387048 0.762400 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -1.389585 0.763313 \n",
"1 Harwood Ave S , Ajax, Ontario -1.379248 0.765519 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -1.386906 0.764920 "
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"coords_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"coords_df['combined'] = list(zip(coords_df['longitude'],coords_df['latitude']))"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" longitude | \n",
" latitude | \n",
" long_rad | \n",
" lat_rad | \n",
"
\n",
" \n",
" combined | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" (-79.877743, 43.6495942) | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
"
\n",
" \n",
" (-79.4275737, 43.6709778) | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
"
\n",
" \n",
" (-79.41025359999999, 43.8575578) | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
"
\n",
" \n",
" (-79.4099732, 43.75732610000001) | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
"
\n",
" \n",
" (-79.4098422, 43.8571536) | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" longitude latitude long_rad lat_rad\n",
"combined \n",
"(-79.877743, 43.6495942) 1 1 1 1\n",
"(-79.4275737, 43.6709778) 1 1 1 1\n",
"(-79.41025359999999, 43.8575578) 1 1 1 1\n",
"(-79.4099732, 43.75732610000001) 1 1 1 1\n",
"(-79.4098422, 43.8571536) 1 1 1 1"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"coords_df.groupby('combined').count().sort_values('longitude',ascending=False)[:5]"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" type_B | \n",
" type_C | \n",
" type_D | \n",
" address | \n",
" longitude | \n",
" latitude | \n",
"
\n",
" \n",
" \n",
" \n",
" 490 Great Lakes Dr , Brampton, ON, Ontario | \n",
" Shell | \n",
" 118.9 | \n",
" 132.9 | \n",
" 140.9 | \n",
" 118.9 | \n",
" 490 Great Lakes Dr, Brampton, ON L6R, Canada | \n",
" -79.773645 | \n",
" 43.740888 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A type_B type_C \\\n",
"490 Great Lakes Dr , Brampton, ON, Ontario Shell 118.9 132.9 140.9 \n",
"\n",
" type_D \\\n",
"490 Great Lakes Dr , Brampton, ON, Ontario 118.9 \n",
"\n",
" address \\\n",
"490 Great Lakes Dr , Brampton, ON, Ontario 490 Great Lakes Dr, Brampton, ON L6R, Canada \n",
"\n",
" longitude latitude \n",
"490 Great Lakes Dr , Brampton, ON, Ontario -79.773645 43.740888 "
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gas_df[coords_df['combined']==(-79.7736453, 43.7408885)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Checking for gas stations with competitors within 100 m\n",
"\n",
"While competition pushes prices down, especially when there is almost 100% price transparency (all gas stations advertise their regular gas prices prominently), multiple gas stations are usually at very busy intersections. In that case, can the extra demand counteract the competitive effect? We can check."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"within_100 = []\n",
"for i in range(len(coords_df)):\n",
" closest_coord = xy_closest(coords_df['combined'][coords_df.index != coords_df.index[i]], coords_df.iloc[i]['combined'])\n",
" within_100.append(xy_distance(coords_df.iloc[i]['combined'][0],coords_df.iloc[i]['combined'][1], closest_coord[0], closest_coord[1]) <= 0.1)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"133"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sum(within_100)"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"coords_df['within_100'] = np.asarray(within_100).astype(int)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" longitude | \n",
" latitude | \n",
" long_rad | \n",
" lat_rad | \n",
" combined | \n",
" within_100 | \n",
"
\n",
" \n",
" \n",
" \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario | \n",
" -79.282378 | \n",
" 43.976346 | \n",
" -1.383739 | \n",
" 0.767532 | \n",
" (-79.282378, 43.9763464) | \n",
" 0 | \n",
"
\n",
" \n",
" 1610 Keele St , Toronto - Central, Ontario | \n",
" -79.471988 | \n",
" 43.682320 | \n",
" -1.387048 | \n",
" 0.762400 | \n",
" (-79.4719879, 43.6823199) | \n",
" 0 | \n",
"
\n",
" \n",
" 6897 Finch Ave W , Toronto - West, Ontario | \n",
" -79.617364 | \n",
" 43.734632 | \n",
" -1.389585 | \n",
" 0.763313 | \n",
" (-79.6173637, 43.7346325) | \n",
" 0 | \n",
"
\n",
" \n",
" 1 Harwood Ave S , Ajax, Ontario | \n",
" -79.025100 | \n",
" 43.861027 | \n",
" -1.379248 | \n",
" 0.765519 | \n",
" (-79.0251004, 43.8610274) | \n",
" 1 | \n",
"
\n",
" \n",
" 1 Thornhill Woods Dr , Vaughan, Ontario | \n",
" -79.463886 | \n",
" 43.826660 | \n",
" -1.386906 | \n",
" 0.764920 | \n",
" (-79.46388639999999, 43.8266602) | \n",
" 0 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" longitude latitude \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n",
" 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n",
"1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n",
"\n",
" long_rad lat_rad \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -1.383739 0.767532 \n",
" 1610 Keele St , Toronto - Central, Ontario -1.387048 0.762400 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -1.389585 0.763313 \n",
"1 Harwood Ave S , Ajax, Ontario -1.379248 0.765519 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -1.386906 0.764920 \n",
"\n",
" combined \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario (-79.282378, 43.9763464) \n",
" 1610 Keele St , Toronto - Central, Ontario (-79.4719879, 43.6823199) \n",
" 6897 Finch Ave W , Toronto - West, Ontario (-79.6173637, 43.7346325) \n",
"1 Harwood Ave S , Ajax, Ontario (-79.0251004, 43.8610274) \n",
"1 Thornhill Woods Dr , Vaughan, Ontario (-79.46388639999999, 43.8266602) \n",
"\n",
" within_100 \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n",
" 1610 Keele St , Toronto - Central, Ontario 0 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 0 \n",
"1 Harwood Ave S , Ajax, Ontario 1 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 0 "
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"coords_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"gas_df['within_100'] = coords_df['within_100']"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" type_B | \n",
" type_C | \n",
" type_D | \n",
" address | \n",
" longitude | \n",
" latitude | \n",
" within_100 | \n",
"
\n",
" \n",
" \n",
" \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario | \n",
" Ultramar | \n",
" 115.6 | \n",
" 129.6 | \n",
" 132.6 | \n",
" 115.9 | \n",
" 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... | \n",
" -79.282378 | \n",
" 43.976346 | \n",
" 0 | \n",
"
\n",
" \n",
" 1610 Keele St , Toronto - Central, Ontario | \n",
" Shell | \n",
" 122.9 | \n",
" 136.9 | \n",
" 144.9 | \n",
" NaN | \n",
" 1610 Keele St, York, ON M6M 3V9, Canada | \n",
" -79.471988 | \n",
" 43.682320 | \n",
" 0 | \n",
"
\n",
" \n",
" 6897 Finch Ave W , Toronto - West, Ontario | \n",
" Esso | \n",
" 121.9 | \n",
" 132.9 | \n",
" 142.9 | \n",
" 118.9 | \n",
" 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada | \n",
" -79.617364 | \n",
" 43.734632 | \n",
" 0 | \n",
"
\n",
" \n",
" 1 Harwood Ave S , Ajax, Ontario | \n",
" Pioneer | \n",
" 122.6 | \n",
" 123.4 | \n",
" 130.9 | \n",
" 118.6 | \n",
" 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada | \n",
" -79.025100 | \n",
" 43.861027 | \n",
" 1 | \n",
"
\n",
" \n",
" 1 Thornhill Woods Dr , Vaughan, Ontario | \n",
" Esso | \n",
" 122.9 | \n",
" NaN | \n",
" NaN | \n",
" 118.9 | \n",
" 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... | \n",
" -79.463886 | \n",
" 43.826660 | \n",
" 0 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A type_B \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 129.6 \n",
" 1610 Keele St , Toronto - Central, Ontario Shell 122.9 136.9 \n",
" 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 132.9 \n",
"1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 123.4 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 NaN \n",
"\n",
" type_C type_D \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 132.6 115.9 \n",
" 1610 Keele St , Toronto - Central, Ontario 144.9 NaN \n",
" 6897 Finch Ave W , Toronto - West, Ontario 142.9 118.9 \n",
"1 Harwood Ave S , Ajax, Ontario 130.9 118.6 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario NaN 118.9 \n",
"\n",
" address \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... \n",
" 1610 Keele St , Toronto - Central, Ontario 1610 Keele St, York, ON M6M 3V9, Canada \n",
" 6897 Finch Ave W , Toronto - West, Ontario 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada \n",
"1 Harwood Ave S , Ajax, Ontario 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... \n",
"\n",
" longitude latitude \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n",
" 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n",
"1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n",
"\n",
" within_100 \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n",
" 1610 Keele St , Toronto - Central, Ontario 0 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 0 \n",
"1 Harwood Ave S , Ajax, Ontario 1 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 0 "
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gas_df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For general clustering, use longitude and latitude. For gasoline-specific clustering, also use the price of the specific type of gas. We can just use k means for general clustering, but when looking at only certain types of gas (such as regular type_A)."
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"cluster_df = gas_df[['name','type_A','longitude','latitude']].copy()"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" longitude | \n",
" latitude | \n",
"
\n",
" \n",
" \n",
" \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario | \n",
" Ultramar | \n",
" 115.6 | \n",
" -79.282378 | \n",
" 43.976346 | \n",
"
\n",
" \n",
" 1610 Keele St , Toronto - Central, Ontario | \n",
" Shell | \n",
" 122.9 | \n",
" -79.471988 | \n",
" 43.682320 | \n",
"
\n",
" \n",
" 6897 Finch Ave W , Toronto - West, Ontario | \n",
" Esso | \n",
" 121.9 | \n",
" -79.617364 | \n",
" 43.734632 | \n",
"
\n",
" \n",
" 1 Harwood Ave S , Ajax, Ontario | \n",
" Pioneer | \n",
" 122.6 | \n",
" -79.025100 | \n",
" 43.861027 | \n",
"
\n",
" \n",
" 1 Thornhill Woods Dr , Vaughan, Ontario | \n",
" Esso | \n",
" 122.9 | \n",
" -79.463886 | \n",
" 43.826660 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 \n",
" 1610 Keele St , Toronto - Central, Ontario Shell 122.9 \n",
" 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 \n",
"1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 \n",
"\n",
" longitude latitude \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n",
" 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n",
"1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 "
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cluster_df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"### Exploring the Starbucks Metric\n",
"Use Starbucks' location research to our benefit. As with station proximities to other stations, we will use the \"as the crow flies\" distance.\n",
"\n",
"Starbucks location data from https://opendata.socrata.com/Business/All-Starbucks-Locations-in-the-World/xy4y-c4mk"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" index | \n",
" Store ID | \n",
" Name | \n",
" Brand | \n",
" Store Number | \n",
" Phone Number | \n",
" Ownership Type | \n",
" Street Combined | \n",
" Street 1 | \n",
" Street 2 | \n",
" ... | \n",
" Country Subdivision | \n",
" Country | \n",
" Postal Code | \n",
" Coordinates | \n",
" Latitude | \n",
" Longitude | \n",
" Timezone | \n",
" Current Timezone Offset | \n",
" Olson Timezone | \n",
" First Seen | \n",
"
\n",
" \n",
" \n",
" \n",
" 0 | \n",
" 108 | \n",
" 323 | \n",
" King & Brant | \n",
" Starbucks | \n",
" 4636-99652 | \n",
" 416-596-0101 | \n",
" CO | \n",
" 527-529 King Street West | \n",
" 527-529 King Street West | \n",
" NaN | \n",
" ... | \n",
" ON | \n",
" CA | \n",
" M5V 1K4 | \n",
" (43.6446495056152, -79.3982620239258) | \n",
" 43.644650 | \n",
" -79.398262 | \n",
" Eastern Standard Time | \n",
" -300 | \n",
" GMT-05:00 America/Toronto | \n",
" 12/08/2013 10:41:59 PM | \n",
"
\n",
" \n",
" 1 | \n",
" 122 | \n",
" 335 | \n",
" 3035 Argentia Rd | \n",
" Starbucks | \n",
" 4665-100880 | \n",
" 905-785-7667 | \n",
" CO | \n",
" 3035 Argentia Road, unIT 7 | \n",
" 3035 Argentia Road | \n",
" unIT 7 | \n",
" ... | \n",
" ON | \n",
" CA | \n",
" L5N 8P7 | \n",
" (43.5956649780273, -79.7863235473633) | \n",
" 43.595665 | \n",
" -79.786324 | \n",
" Eastern Standard Time | \n",
" -300 | \n",
" GMT-05:00 America/Toronto | \n",
" 12/08/2013 10:41:59 PM | \n",
"
\n",
" \n",
" 2 | \n",
" 124 | \n",
" 337 | \n",
" 631 Commissioners Rd E | \n",
" Starbucks | \n",
" 4674-100799 | \n",
" 519-649-1444 | \n",
" CO | \n",
" 631 Commissioners Road East | \n",
" 631 Commissioners Road East | \n",
" NaN | \n",
" ... | \n",
" ON | \n",
" CA | \n",
" N6C 2V1 | \n",
" (42.957763671875, -81.2328872680664) | \n",
" 42.957764 | \n",
" -81.232887 | \n",
" Eastern Standard Time | \n",
" -300 | \n",
" GMT-05:00 America/Toronto | \n",
" 12/08/2013 10:41:59 PM | \n",
"
\n",
" \n",
" 3 | \n",
" 133 | \n",
" 354 | \n",
" 599 Taylor Kidd Blvd | \n",
" Starbucks | \n",
" 4694-102072 | \n",
" 613-634-1509 | \n",
" CO | \n",
" 599 Taylor Kidd Blvd, Unit 6 | \n",
" 599 Taylor Kidd Blvd | \n",
" Unit 6 | \n",
" ... | \n",
" ON | \n",
" CA | \n",
" K7M 0A2 | \n",
" (44.2504081726074, -76.5673141479492) | \n",
" 44.250408 | \n",
" -76.567314 | \n",
" Eastern Standard Time | \n",
" -300 | \n",
" GMT-05:00 America/Toronto | \n",
" 12/08/2013 10:41:59 PM | \n",
"
\n",
" \n",
" 4 | \n",
" 141 | \n",
" 364 | \n",
" King & Shaw | \n",
" Starbucks | \n",
" 4646-98483 | \n",
" 416-979-1953 | \n",
" CO | \n",
" 1005 King St W, Unit 7 | \n",
" 1005 King St W | \n",
" Unit 7 | \n",
" ... | \n",
" ON | \n",
" CA | \n",
" M6G 1B9 | \n",
" (43.6412811279297, -79.4147644042969) | \n",
" 43.641281 | \n",
" -79.414764 | \n",
" Eastern Standard Time | \n",
" -300 | \n",
" GMT-05:00 America/Toronto | \n",
" 12/08/2013 10:41:59 PM | \n",
"
\n",
" \n",
"
\n",
"
5 rows × 22 columns
\n",
"
"
],
"text/plain": [
" index Store ID Name Brand Store Number \\\n",
"0 108 323 King & Brant Starbucks 4636-99652 \n",
"1 122 335 3035 Argentia Rd Starbucks 4665-100880 \n",
"2 124 337 631 Commissioners Rd E Starbucks 4674-100799 \n",
"3 133 354 599 Taylor Kidd Blvd Starbucks 4694-102072 \n",
"4 141 364 King & Shaw Starbucks 4646-98483 \n",
"\n",
" Phone Number Ownership Type Street Combined \\\n",
"0 416-596-0101 CO 527-529 King Street West \n",
"1 905-785-7667 CO 3035 Argentia Road, unIT 7 \n",
"2 519-649-1444 CO 631 Commissioners Road East \n",
"3 613-634-1509 CO 599 Taylor Kidd Blvd, Unit 6 \n",
"4 416-979-1953 CO 1005 King St W, Unit 7 \n",
"\n",
" Street 1 Street 2 ... \\\n",
"0 527-529 King Street West NaN ... \n",
"1 3035 Argentia Road unIT 7 ... \n",
"2 631 Commissioners Road East NaN ... \n",
"3 599 Taylor Kidd Blvd Unit 6 ... \n",
"4 1005 King St W Unit 7 ... \n",
"\n",
" Country Subdivision Country Postal Code \\\n",
"0 ON CA M5V 1K4 \n",
"1 ON CA L5N 8P7 \n",
"2 ON CA N6C 2V1 \n",
"3 ON CA K7M 0A2 \n",
"4 ON CA M6G 1B9 \n",
"\n",
" Coordinates Latitude Longitude \\\n",
"0 (43.6446495056152, -79.3982620239258) 43.644650 -79.398262 \n",
"1 (43.5956649780273, -79.7863235473633) 43.595665 -79.786324 \n",
"2 (42.957763671875, -81.2328872680664) 42.957764 -81.232887 \n",
"3 (44.2504081726074, -76.5673141479492) 44.250408 -76.567314 \n",
"4 (43.6412811279297, -79.4147644042969) 43.641281 -79.414764 \n",
"\n",
" Timezone Current Timezone Offset Olson Timezone \\\n",
"0 Eastern Standard Time -300 GMT-05:00 America/Toronto \n",
"1 Eastern Standard Time -300 GMT-05:00 America/Toronto \n",
"2 Eastern Standard Time -300 GMT-05:00 America/Toronto \n",
"3 Eastern Standard Time -300 GMT-05:00 America/Toronto \n",
"4 Eastern Standard Time -300 GMT-05:00 America/Toronto \n",
"\n",
" First Seen \n",
"0 12/08/2013 10:41:59 PM \n",
"1 12/08/2013 10:41:59 PM \n",
"2 12/08/2013 10:41:59 PM \n",
"3 12/08/2013 10:41:59 PM \n",
"4 12/08/2013 10:41:59 PM \n",
"\n",
"[5 rows x 22 columns]"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"starbucks_df = pd.read_csv('All_Starbucks_Locations_in_the_World.csv', delimiter=';')\n",
"starbucks_df = starbucks_df[(starbucks_df['Country']=='CA') & (starbucks_df['Country Subdivision']=='ON')]\n",
"starbucks_df.reset_index(inplace=True)\n",
"starbucks_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"RangeIndex: 530 entries, 0 to 529\n",
"Data columns (total 22 columns):\n",
"index 530 non-null int64\n",
"Store ID 530 non-null object\n",
"Name 530 non-null object\n",
"Brand 530 non-null object\n",
"Store Number 530 non-null object\n",
"Phone Number 486 non-null object\n",
"Ownership Type 530 non-null object\n",
"Street Combined 530 non-null object\n",
"Street 1 530 non-null object\n",
"Street 2 191 non-null object\n",
"Street 3 110 non-null object\n",
"City 530 non-null object\n",
"Country Subdivision 530 non-null object\n",
"Country 530 non-null object\n",
"Postal Code 530 non-null object\n",
"Coordinates 530 non-null object\n",
"Latitude 530 non-null float64\n",
"Longitude 530 non-null float64\n",
"Timezone 530 non-null object\n",
"Current Timezone Offset 530 non-null int64\n",
"Olson Timezone 530 non-null object\n",
"First Seen 530 non-null object\n",
"dtypes: float64(2), int64(2), object(18)\n",
"memory usage: 91.2+ KB\n"
]
}
],
"source": [
"starbucks_df.info()"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Actually we only need the coordinates.\n",
"starbucks_on_coords = starbucks_df['Coordinates'].tolist()"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Convert from a list of strings that look like Python lists to an actual list of lists\n",
"# Also need to change it from (lat,long) to (long,lat)\n",
"starbucks_on_coords = [ast.literal_eval(starbucks_on_coords[i]) for i in range(len(starbucks_on_coords))]\n",
"starbucks_on_coords = [(starbucks_on_coords[i][1], starbucks_on_coords[i][0]) for i in range(len(starbucks_on_coords))]"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"nearest_starbucks = []\n",
"for i in range(len(coords_df)):\n",
" closest_coord = xy_closest(starbucks_on_coords, coords_df.iloc[i]['combined'])\n",
" nearest_starbucks.append(xy_distance(coords_df.iloc[i]['combined'][1],coords_df.iloc[i]['combined'][0], closest_coord[1], closest_coord[0]))"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"starbucks_within_100m = []\n",
"for i in range(len(coords_df)):\n",
" closest_coord = xy_closest(starbucks_on_coords, coords_df.iloc[i]['combined'])\n",
" starbucks_within_100m.append(xy_distance(coords_df.iloc[i]['combined'][1],coords_df.iloc[i]['combined'][0], closest_coord[1], closest_coord[0])<=1)"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"gas_df['nearest_starbucks'] = np.asarray(nearest_starbucks)\n",
"gas_df['starbucks_within_100m'] = np.asarray(starbucks_within_100m).astype(int)"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" type_B | \n",
" type_C | \n",
" type_D | \n",
" address | \n",
" longitude | \n",
" latitude | \n",
" within_100 | \n",
" nearest_starbucks | \n",
" starbucks_within_100m | \n",
"
\n",
" \n",
" \n",
" \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario | \n",
" Ultramar | \n",
" 115.6 | \n",
" 129.6 | \n",
" 132.6 | \n",
" 115.9 | \n",
" 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... | \n",
" -79.282378 | \n",
" 43.976346 | \n",
" 0 | \n",
" 11.208629 | \n",
" 0 | \n",
"
\n",
" \n",
" 1610 Keele St , Toronto - Central, Ontario | \n",
" Shell | \n",
" 122.9 | \n",
" 136.9 | \n",
" 144.9 | \n",
" NaN | \n",
" 1610 Keele St, York, ON M6M 3V9, Canada | \n",
" -79.471988 | \n",
" 43.682320 | \n",
" 0 | \n",
" 1.178252 | \n",
" 0 | \n",
"
\n",
" \n",
" 6897 Finch Ave W , Toronto - West, Ontario | \n",
" Esso | \n",
" 121.9 | \n",
" 132.9 | \n",
" 142.9 | \n",
" 118.9 | \n",
" 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada | \n",
" -79.617364 | \n",
" 43.734632 | \n",
" 0 | \n",
" 1.031603 | \n",
" 0 | \n",
"
\n",
" \n",
" 1 Harwood Ave S , Ajax, Ontario | \n",
" Pioneer | \n",
" 122.6 | \n",
" 123.4 | \n",
" 130.9 | \n",
" 118.6 | \n",
" 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada | \n",
" -79.025100 | \n",
" 43.861027 | \n",
" 1 | \n",
" 0.381652 | \n",
" 1 | \n",
"
\n",
" \n",
" 1 Thornhill Woods Dr , Vaughan, Ontario | \n",
" Esso | \n",
" 122.9 | \n",
" NaN | \n",
" NaN | \n",
" 118.9 | \n",
" 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... | \n",
" -79.463886 | \n",
" 43.826660 | \n",
" 0 | \n",
" 1.444566 | \n",
" 0 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A type_B \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 129.6 \n",
" 1610 Keele St , Toronto - Central, Ontario Shell 122.9 136.9 \n",
" 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 132.9 \n",
"1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 123.4 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 NaN \n",
"\n",
" type_C type_D \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 132.6 115.9 \n",
" 1610 Keele St , Toronto - Central, Ontario 144.9 NaN \n",
" 6897 Finch Ave W , Toronto - West, Ontario 142.9 118.9 \n",
"1 Harwood Ave S , Ajax, Ontario 130.9 118.6 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario NaN 118.9 \n",
"\n",
" address \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... \n",
" 1610 Keele St , Toronto - Central, Ontario 1610 Keele St, York, ON M6M 3V9, Canada \n",
" 6897 Finch Ave W , Toronto - West, Ontario 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada \n",
"1 Harwood Ave S , Ajax, Ontario 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... \n",
"\n",
" longitude latitude \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n",
" 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n",
"1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n",
"\n",
" within_100 \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n",
" 1610 Keele St , Toronto - Central, Ontario 0 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 0 \n",
"1 Harwood Ave S , Ajax, Ontario 1 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 0 \n",
"\n",
" nearest_starbucks \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 11.208629 \n",
" 1610 Keele St , Toronto - Central, Ontario 1.178252 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 1.031603 \n",
"1 Harwood Ave S , Ajax, Ontario 0.381652 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 1.444566 \n",
"\n",
" starbucks_within_100m \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n",
" 1610 Keele St , Toronto - Central, Ontario 0 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 0 \n",
"1 Harwood Ave S , Ajax, Ontario 1 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 0 "
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gas_df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's invent a Starbucks proximity metric, a sigmoid with the log of the inverse of the nearest starbucks in km:
\n",
"$\\frac{1}{(1+exp(log(-\\frac{1}{d})^2))}$\n",
"(This is not actually used, but keep this here just in case)"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"gas_df['starbucks_metric'] = gas_df['nearest_starbucks'].apply(lambda d:1/(1+np.exp(-np.log((2/d)**2))))"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" type_B | \n",
" type_C | \n",
" type_D | \n",
" address | \n",
" longitude | \n",
" latitude | \n",
" within_100 | \n",
" nearest_starbucks | \n",
" starbucks_within_100m | \n",
" starbucks_metric | \n",
"
\n",
" \n",
" \n",
" \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario | \n",
" Ultramar | \n",
" 115.6 | \n",
" 129.6 | \n",
" 132.6 | \n",
" 115.9 | \n",
" 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... | \n",
" -79.282378 | \n",
" 43.976346 | \n",
" 0 | \n",
" 11.208629 | \n",
" 0 | \n",
" 0.030856 | \n",
"
\n",
" \n",
" 1610 Keele St , Toronto - Central, Ontario | \n",
" Shell | \n",
" 122.9 | \n",
" 136.9 | \n",
" 144.9 | \n",
" NaN | \n",
" 1610 Keele St, York, ON M6M 3V9, Canada | \n",
" -79.471988 | \n",
" 43.682320 | \n",
" 0 | \n",
" 1.178252 | \n",
" 0 | \n",
" 0.742352 | \n",
"
\n",
" \n",
" 6897 Finch Ave W , Toronto - West, Ontario | \n",
" Esso | \n",
" 121.9 | \n",
" 132.9 | \n",
" 142.9 | \n",
" 118.9 | \n",
" 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada | \n",
" -79.617364 | \n",
" 43.734632 | \n",
" 0 | \n",
" 1.031603 | \n",
" 0 | \n",
" 0.789857 | \n",
"
\n",
" \n",
" 1 Harwood Ave S , Ajax, Ontario | \n",
" Pioneer | \n",
" 122.6 | \n",
" 123.4 | \n",
" 130.9 | \n",
" 118.6 | \n",
" 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada | \n",
" -79.025100 | \n",
" 43.861027 | \n",
" 1 | \n",
" 0.381652 | \n",
" 1 | \n",
" 0.964865 | \n",
"
\n",
" \n",
" 1 Thornhill Woods Dr , Vaughan, Ontario | \n",
" Esso | \n",
" 122.9 | \n",
" NaN | \n",
" NaN | \n",
" 118.9 | \n",
" 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... | \n",
" -79.463886 | \n",
" 43.826660 | \n",
" 0 | \n",
" 1.444566 | \n",
" 0 | \n",
" 0.657163 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A type_B \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 129.6 \n",
" 1610 Keele St , Toronto - Central, Ontario Shell 122.9 136.9 \n",
" 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 132.9 \n",
"1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 123.4 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 NaN \n",
"\n",
" type_C type_D \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 132.6 115.9 \n",
" 1610 Keele St , Toronto - Central, Ontario 144.9 NaN \n",
" 6897 Finch Ave W , Toronto - West, Ontario 142.9 118.9 \n",
"1 Harwood Ave S , Ajax, Ontario 130.9 118.6 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario NaN 118.9 \n",
"\n",
" address \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... \n",
" 1610 Keele St , Toronto - Central, Ontario 1610 Keele St, York, ON M6M 3V9, Canada \n",
" 6897 Finch Ave W , Toronto - West, Ontario 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada \n",
"1 Harwood Ave S , Ajax, Ontario 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... \n",
"\n",
" longitude latitude \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n",
" 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n",
"1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n",
"\n",
" within_100 \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n",
" 1610 Keele St , Toronto - Central, Ontario 0 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 0 \n",
"1 Harwood Ave S , Ajax, Ontario 1 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 0 \n",
"\n",
" nearest_starbucks \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 11.208629 \n",
" 1610 Keele St , Toronto - Central, Ontario 1.178252 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 1.031603 \n",
"1 Harwood Ave S , Ajax, Ontario 0.381652 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 1.444566 \n",
"\n",
" starbucks_within_100m \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n",
" 1610 Keele St , Toronto - Central, Ontario 0 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 0 \n",
"1 Harwood Ave S , Ajax, Ontario 1 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 0 \n",
"\n",
" starbucks_metric \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0.030856 \n",
" 1610 Keele St , Toronto - Central, Ontario 0.742352 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 0.789857 \n",
"1 Harwood Ave S , Ajax, Ontario 0.964865 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 0.657163 "
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gas_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"gas_df.to_csv('gas_stations_2017-12-13.csv')"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"regular_gas = gas_df[['type_A','within_100','starbucks_within_100m']].dropna().copy()"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"regular_gas.columns=['Price','Competition_Nearby','Starbucks_Nearby']"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"regular_gas['Relative_Price'] = regular_gas['Price'].apply(lambda p: p/np.mean(regular_gas['Price']))\n",
"regular_gas['Price_Difference'] = regular_gas['Price'].apply(lambda p: p-np.mean(regular_gas['Price']))"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" Price | \n",
" Competition_Nearby | \n",
" Starbucks_Nearby | \n",
" Relative_Price | \n",
" Price_Difference | \n",
"
\n",
" \n",
" \n",
" \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario | \n",
" 115.6 | \n",
" 0 | \n",
" 0 | \n",
" 0.956242 | \n",
" -5.289926 | \n",
"
\n",
" \n",
" 1610 Keele St , Toronto - Central, Ontario | \n",
" 122.9 | \n",
" 0 | \n",
" 0 | \n",
" 1.016627 | \n",
" 2.010074 | \n",
"
\n",
" \n",
" 6897 Finch Ave W , Toronto - West, Ontario | \n",
" 121.9 | \n",
" 0 | \n",
" 0 | \n",
" 1.008355 | \n",
" 1.010074 | \n",
"
\n",
" \n",
" 1 Harwood Ave S , Ajax, Ontario | \n",
" 122.6 | \n",
" 1 | \n",
" 1 | \n",
" 1.014146 | \n",
" 1.710074 | \n",
"
\n",
" \n",
" 1 Thornhill Woods Dr , Vaughan, Ontario | \n",
" 122.9 | \n",
" 0 | \n",
" 0 | \n",
" 1.016627 | \n",
" 2.010074 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Price Competition_Nearby \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 115.6 0 \n",
" 1610 Keele St , Toronto - Central, Ontario 122.9 0 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 121.9 0 \n",
"1 Harwood Ave S , Ajax, Ontario 122.6 1 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 122.9 0 \n",
"\n",
" Starbucks_Nearby \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n",
" 1610 Keele St , Toronto - Central, Ontario 0 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 0 \n",
"1 Harwood Ave S , Ajax, Ontario 1 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 0 \n",
"\n",
" Relative_Price \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0.956242 \n",
" 1610 Keele St , Toronto - Central, Ontario 1.016627 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 1.008355 \n",
"1 Harwood Ave S , Ajax, Ontario 1.014146 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 1.016627 \n",
"\n",
" Price_Difference \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -5.289926 \n",
" 1610 Keele St , Toronto - Central, Ontario 2.010074 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 1.010074 \n",
"1 Harwood Ave S , Ajax, Ontario 1.710074 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 2.010074 "
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"regular_gas.head()"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\ProgramData\\Anaconda3\\lib\\site-packages\\statsmodels\\compat\\pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n",
" from pandas.core import datetools\n"
]
}
],
"source": [
"import statsmodels.api as sm"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"Y = regular_gas['Price']\n",
"X = sm.add_constant(regular_gas[['Competition_Nearby','Starbucks_Nearby']])\n",
"lm = sm.OLS(Y,X)"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"results = lm.fit()"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"const 120.606227\n",
"Competition_Nearby 0.182607\n",
"Starbucks_Nearby 0.621326\n",
"dtype: float64"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"results.params"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"OLS Regression Results\n",
"\n",
" Dep. Variable: | Price | R-squared: | 0.011 | \n",
"
\n",
"\n",
" Model: | OLS | Adj. R-squared: | 0.008 | \n",
"
\n",
"\n",
" Method: | Least Squares | F-statistic: | 3.873 | \n",
"
\n",
"\n",
" Date: | Sat, 06 Jan 2018 | Prob (F-statistic): | 0.0213 | \n",
"
\n",
"\n",
" Time: | 21:37:31 | Log-Likelihood: | -1681.8 | \n",
"
\n",
"\n",
" No. Observations: | 675 | AIC: | 3370. | \n",
"
\n",
"\n",
" Df Residuals: | 672 | BIC: | 3383. | \n",
"
\n",
"\n",
" Df Model: | 2 | | | \n",
"
\n",
"\n",
" Covariance Type: | nonrobust | | | \n",
"
\n",
"
\n",
"\n",
"\n",
" | coef | std err | t | P>|t| | [0.025 | 0.975] | \n",
"
\n",
"\n",
" const | 120.6062 | 0.155 | 777.316 | 0.000 | 120.302 | 120.911 | \n",
"
\n",
"\n",
" Competition_Nearby | 0.1826 | 0.286 | 0.638 | 0.523 | -0.379 | 0.744 | \n",
"
\n",
"\n",
" Starbucks_Nearby | 0.6213 | 0.230 | 2.699 | 0.007 | 0.169 | 1.073 | \n",
"
\n",
"
\n",
"\n",
"\n",
" Omnibus: | 132.421 | Durbin-Watson: | 1.985 | \n",
"
\n",
"\n",
" Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 211.951 | \n",
"
\n",
"\n",
" Skew: | -1.304 | Prob(JB): | 9.45e-47 | \n",
"
\n",
"\n",
" Kurtosis: | 3.855 | Cond. No. | 2.90 | \n",
"
\n",
"
"
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Price R-squared: 0.011\n",
"Model: OLS Adj. R-squared: 0.008\n",
"Method: Least Squares F-statistic: 3.873\n",
"Date: Sat, 06 Jan 2018 Prob (F-statistic): 0.0213\n",
"Time: 21:37:31 Log-Likelihood: -1681.8\n",
"No. Observations: 675 AIC: 3370.\n",
"Df Residuals: 672 BIC: 3383.\n",
"Df Model: 2 \n",
"Covariance Type: nonrobust \n",
"======================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"--------------------------------------------------------------------------------------\n",
"const 120.6062 0.155 777.316 0.000 120.302 120.911\n",
"Competition_Nearby 0.1826 0.286 0.638 0.523 -0.379 0.744\n",
"Starbucks_Nearby 0.6213 0.230 2.699 0.007 0.169 1.073\n",
"==============================================================================\n",
"Omnibus: 132.421 Durbin-Watson: 1.985\n",
"Prob(Omnibus): 0.000 Jarque-Bera (JB): 211.951\n",
"Skew: -1.304 Prob(JB): 9.45e-47\n",
"Kurtosis: 3.855 Cond. No. 2.90\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"results.summary()"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"Gas stations within 100m of a Starbucks location charge 0.62 cents more for regular gas, and this is statistically significant with over 99% confidence.\n",
"\n",
"Proximity to competition also seems to increase prices, likely due to increased demand at busy intersections (as mentioned above). However, the p value is far too high so we can not be certain of its effect with any confidence."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Clustering and Mapping\n",
"\n",
"Let's see if we can find any meaningful clusters."
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"cluster_df.dropna(inplace=True)"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Index: 675 entries, 12731 Hwy 48 , Whitchurch-Stouffville, Ontario to Taymall Ave , Toronto - West, Ontario\n",
"Data columns (total 4 columns):\n",
"name 675 non-null object\n",
"type_A 675 non-null float64\n",
"longitude 675 non-null float64\n",
"latitude 675 non-null float64\n",
"dtypes: float64(3), object(1)\n",
"memory usage: 26.4+ KB\n"
]
}
],
"source": [
"cluster_df.info()"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"cluster_df['type_A_norm'] = cluster_df['type_A']/np.mean(cluster_df['type_A'])\n",
"cluster_df['longitude_norm'] = -cluster_df['longitude']/np.mean(cluster_df['longitude'])\n",
"cluster_df['latitude_norm'] = cluster_df['latitude']/np.mean(cluster_df['latitude'])"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" longitude | \n",
" latitude | \n",
" type_A_norm | \n",
" longitude_norm | \n",
" latitude_norm | \n",
"
\n",
" \n",
" \n",
" \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario | \n",
" Ultramar | \n",
" 115.6 | \n",
" -79.282378 | \n",
" 43.976346 | \n",
" 0.956242 | \n",
" -0.997435 | \n",
" 1.005572 | \n",
"
\n",
" \n",
" 1610 Keele St , Toronto - Central, Ontario | \n",
" Shell | \n",
" 122.9 | \n",
" -79.471988 | \n",
" 43.682320 | \n",
" 1.016627 | \n",
" -0.999820 | \n",
" 0.998849 | \n",
"
\n",
" \n",
" 6897 Finch Ave W , Toronto - West, Ontario | \n",
" Esso | \n",
" 121.9 | \n",
" -79.617364 | \n",
" 43.734632 | \n",
" 1.008355 | \n",
" -1.001649 | \n",
" 1.000045 | \n",
"
\n",
" \n",
" 1 Harwood Ave S , Ajax, Ontario | \n",
" Pioneer | \n",
" 122.6 | \n",
" -79.025100 | \n",
" 43.861027 | \n",
" 1.014146 | \n",
" -0.994198 | \n",
" 1.002935 | \n",
"
\n",
" \n",
" 1 Thornhill Woods Dr , Vaughan, Ontario | \n",
" Esso | \n",
" 122.9 | \n",
" -79.463886 | \n",
" 43.826660 | \n",
" 1.016627 | \n",
" -0.999718 | \n",
" 1.002149 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 \n",
" 1610 Keele St , Toronto - Central, Ontario Shell 122.9 \n",
" 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 \n",
"1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 \n",
"\n",
" longitude latitude \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n",
" 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n",
"1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n",
"\n",
" type_A_norm \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0.956242 \n",
" 1610 Keele St , Toronto - Central, Ontario 1.016627 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 1.008355 \n",
"1 Harwood Ave S , Ajax, Ontario 1.014146 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 1.016627 \n",
"\n",
" longitude_norm \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -0.997435 \n",
" 1610 Keele St , Toronto - Central, Ontario -0.999820 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -1.001649 \n",
"1 Harwood Ave S , Ajax, Ontario -0.994198 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -0.999718 \n",
"\n",
" latitude_norm \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 1.005572 \n",
" 1610 Keele St , Toronto - Central, Ontario 0.998849 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 1.000045 \n",
"1 Harwood Ave S , Ajax, Ontario 1.002935 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 1.002149 "
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cluster_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from scipy import linalg\n",
"from sklearn.cluster import KMeans\n",
"from sklearn import svm\n",
"from scipy.spatial.distance import cdist\n",
"from scipy.spatial import ConvexHull\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0.5,1,'SSE vs k')"
]
},
"execution_count": 65,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEWCAYAAABbgYH9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XmcHVWd9/HPNyuEJYQkaCCQZoli\nFAzShs4DuEUQFAgIjsEo+IgEFWbcUEFkRniGGcANGdCZAGqEQKIoEmdEdhhFCDQ7IQJNCBACodlC\nQoBsv+ePU03fdG6vudXV9/b3/Xrd1606dercX41M/3JOVZ2jiMDMzKzSBhQdgJmZ1SYnGDMzy4UT\njJmZ5cIJxszMcuEEY2ZmuXCCMTOzXDjBmNU4SXWSQtKgomOx/sUJxqwHJO0n6W+Slkt6SdJtkt6f\nHRsi6UeSlkhaKekJST8pOXexpNezYy2fC4q7GrN8+F80Zt0kaWvgv4EvA78BhgD7A29mVU4F6oFJ\nwLPAOOADbZo5NCJu6JWAzQriHoxZ970DICKuiIh1EfF6RFwXEQ9kx98PXBURSyNZHBG/7u6PSNo+\n6+lsW1K2l6QXJA2WtJukW7Ne1AuS5nax3SOzXtR7uhuTWXc4wZh136PAOkmzJB0saUSb43cA35D0\nFUl7SFJPfiQilgK3A0eWFH8GuDIi1gD/D7gOGAGMBf6jszYl/V/gHOCjEfFQT+Iy6yonGLNuiohX\ngf2AAC4CmiXNk/S2rMq/k/6ITwcagWckHdummT9IeqXkc3w7P3c5cDRAlqimZWUAa0jDb9tHxBsR\n8ddOQv8a8C3gQxHR1NXrNespebJLs00jaXfgMuCxiDi6zbHNgS8A5wPviYiFkhYDX+zKPRhJ25Du\n4+wKjAcuBcZFREh6O6kX8wngZeBHEfGLMm3UAU8AzwNnRsSFPbxUs25xD8ZsE0XE34FfARvd08ju\nz1xISgATetD2K6RhsH8gDY9dEdm/CiPiuYg4PiK2B04AfiZptw6aOxD4nqQjO6hjVjFOMGbdJGl3\nSd+UNDbb35E0jHVHtv81SR+StLmkQdnw2FbAvT38ycuBY0j3YlqGx5D0qZYYSAksgHUdtLMAOAi4\nUNJhPYzFrMucYMy6bwWwDzBf0mukxPIQ8M3s+OvAj4DngBeAE4EjI2JRSRt/bPMezFUd/N480vDY\nsoi4v6T8/VkMK7M6X42IJzoKPDv/EOAiSQd38XrNesT3YMzMLBfuwZiZWS6cYMzMLBdOMGZmlgsn\nGDMzy0W/nuxy1KhRUVdXV3QYZmZV5e67734hIkZ3Vq9fJ5i6ujoaGxuLDsPMrKpIerIr9TxEZmZm\nuXCCMTOzXDjBmJlZLpxgzMwsF04wZmaWCyeYbpo9G+rqYMCA9D17dtERmZn1Tf36MeXumj0bZsyA\nVavS/pNPpn2A6dOLi8vMrC9yD6YbTjutNbm0WLUqlZuZ2YacYLrhqae6V25m1p85wXTDTjt1r9zM\nrD9zgumGs86CYcM2LBs2LJWbmdmGnGC6Yfp0mDkTxo1L+0OGpH3f4Dcz25gTTDdNnw6LF8PJJ4ME\nRx1VdERmZn1TrglG0kGSHpHUJOmUMseHSpqbHZ8vqS4rP0DS3ZIezL4/UnLO3ll5k6TzJSkr31bS\n9ZIey75H5HltH/0oTJ0Ky5fn+StmZtUrtwQjaSBwIXAwMAE4WtKENtWOA16OiN2AnwDnZOUvAIdG\nxB7AscClJef8HJgBjM8+B2XlpwA3RsR44MZsPzcf+xjMnQvbbZfnr5iZVa88ezCTgKaIWBQRq4E5\nwNQ2daYCs7LtK4EpkhQR90bE0qx8AbBZ1tsZA2wdEbdHRAC/Bg4v09askvJcuQdjZlZenglmB+Dp\nkv0lWVnZOhGxFlgOjGxT50jg3oh4M6u/pJ023xYRz2ZtPQuU7VtImiGpUVJjc3Nzty+q1AknwMSJ\nm9SEmVnNyjPBqExZdKeOpHeThs1O6EabHYqImRFRHxH1o0d3uuJnh975znTD/7nnNqkZM7OalGeC\nWQLsWLI/FljaXh1Jg4DhwEvZ/ljgKuCYiHi8pP7Ydtpclg2hkX0/X7Eracfkyen7jjvy/iUzs+qT\nZ4K5CxgvaWdJQ4BpwLw2deaRbuIDHAXcFBEhaRvgf4BTI+K2lsrZ0NcKSQ3Z02PHAFeXaevYkvLc\n7LUXDB7sBGNmVk5uCSa7p3IScC2wEPhNRCyQdKakw7JqlwAjJTUB36D1ya+TgN2A0yXdl31a7ql8\nGbgYaAIeB67Jys8GDpD0GHBAtp+rzTZLScYJxsxsY0oPY/VP9fX10djYuEltXHklRMCnPlWhoMzM\n+jhJd0dEfWf1vB7MJvKb/GZm5XmqmE0UAfffDwsWFB2JmVnf4gRTAQceCOeeW3QUZmZ9ixPMJpKg\nocE3+s3M2nKCqYDJk+HRR+HFF4uOxMys73CCqYCGhvR9553FxmFm1pc4wVRAfT0MGAC33150JGZm\nfYcfU66ALbeEW2+Fd7+76EjMzPoOJ5gK2W+/oiMwM+tbPERWIUuXwr/+KzzxRNGRmJn1DU4wFbJi\nBZx+Otx0U9GRmJn1DU4wFfKOd8CIEX4fxsyshRNMhfiFSzOzDTnBVFBDQ5qTbPnyoiMxMyueE0wF\nNTTA0KHprX4zs/7OjylX0Ic/DK++mla5NDPr75xgKsiJxcyslYfIKmzOnDR9fz9eKNTMDHCCqbgV\nK+D666GpqehIzMyKlWuCkXSQpEckNUk6pczxoZLmZsfnS6rLykdKulnSSkkXlNTfStJ9JZ8XJJ2X\nHfu8pOaSY1/M89ra0zKzsh9XNrP+LrcEI2kgcCFwMDABOFrShDbVjgNejojdgJ8A52TlbwCnAyeX\nVo6IFRExseUDPAn8vqTK3JLjF1f+qjo3YUKa/NIJxsz6uzx7MJOApohYFBGrgTnA1DZ1pgKzsu0r\ngSmSFBGvRcRfSYmmLEnjge2Av1Q+9J4bOBAmTXKCMTPLM8HsADxdsr8kKytbJyLWAsuBkV1s/2hS\nj6X0dvqRkh6QdKWkHcudJGmGpEZJjc3NzV38qe75+Mdht918o9/M+rc8E4zKlLX9k9uVOu2ZBlxR\nsv9HoC4i9gRuoLVntGHjETMjoj4i6kePHt3Fn+qeb34T5s5N08eYmfVXeSaYJUBpL2IssLS9OpIG\nAcOBlzprWNJ7gUERcXdLWUS8GBFvZrsXAXv3PPTKWL266AjMzIqTZ4K5CxgvaWdJQ0g9jnlt6swD\njs22jwJuajPk1Z6j2bD3gqQxJbuHAQt7FHWFfPSjMG1akRGYmRUrtzf5I2KtpJOAa4GBwC8iYoGk\nM4HGiJgHXAJcKqmJ1HN560+ypMXA1sAQSYcDB0bEw9nhfwA+3uYn/0nSYcDarK3P53VtXTFmDNxw\nQ7oP46EyM+uP1LUOQ22qr6+PxsbGXNq+8EI46SRYvBjGjcvlJ8zMCiHp7oio76ye3+TPiV+4NLP+\nzgkmJ3vuCZtv7gRjZv2XZ1POyeDB8P3vw7vfXXQkZmbFcILJ0be/XXQEZmbF8RBZjtavh4cfhmXL\nio7EzKz3OcHk6Lnn0hDZnDlFR2Jm1vucYHK0/faw446+0W9m/ZMTTM4aGuD224uOwsys9znB5Gzy\nZHjySXj22aIjMTPrXU4wOWt54XL+/GLjMDPrbX5MOWd77QV//CPst1/RkZiZ9S4nmJxtthkcckjR\nUZiZ9T4PkfWCxx6Dc8+FNWuKjsTMrPc4wfSCu++G73wHHnqo6EjMzHqPE0wvaLnR78eVzaw/cYLp\nBePGwdvf7hcuzax/cYLpBVLqxTjBmFl/4gTTSxoa0guXy5cXHYmZWe/INcFIOkjSI5KaJJ1S5vhQ\nSXOz4/Ml1WXlIyXdLGmlpAvanHNL1uZ92We7jtrqK77ylZRchg8vOhIzs96RW4KRNBC4EDgYmAAc\nLWlCm2rHAS9HxG7AT4BzsvI3gNOBk9tpfnpETMw+z3fSVp+w1VbpnRgzs/4izx7MJKApIhZFxGpg\nDjC1TZ2pwKxs+0pgiiRFxGsR8VdSoumqsm31PPzKO/98+OY3i47CzKx35JlgdgCeLtlfkpWVrRMR\na4HlwMgutP3LbHjs9JIk0tO2es3ChXDxxWkhMjOzWpdnginXe4ge1GlrekTsAeyffT7XnbYkzZDU\nKKmxubm5k5+qrMmT4dVXU6IxM6t1eSaYJcCOJftjgaXt1ZE0CBgOvNRRoxHxTPa9AricNBTX5bYi\nYmZE1EdE/ejRo7t5SZum5YVLP65sZv1BngnmLmC8pJ0lDQGmAfPa1JkHHJttHwXcFBHt9mAkDZI0\nKtseDBwCtEzA0q22ijB+PIwY4QRjZv1DbrMpR8RaSScB1wIDgV9ExAJJZwKNETEPuAS4VFITqbcx\nreV8SYuBrYEhkg4HDgSeBK7NkstA4AbgouyUdtvqKyT4xCdg882LjsTMLH/qY//I71X19fXR2NhY\ndBhmZlVF0t0RUd9ZPb/JX5B+nNfNrJ9wgullq1fDhAlw1llFR2Jmli8nmF42ZEi6F+Op+82s1jnB\nFGDy5PQkmYfJzKyWOcEUoKEBXnoJmpqKjsTMLD9OMAXwCpdm1h84wRRgwgQ4/njYeeeiIzEzy09u\nL1pa+wYMgJkzi47CzCxf7sEUJCLdg3nzzaIjMTPLhxNMQa65Js1NNn9+0ZGYmeXDCaYgk7I5oD3x\npZnVKieYgowaBbvt5gRjZrXLCaZADQ3pUWW/cGlmtcgJpkANDfDcc/DUU0VHYmZWeX5MuUCHHgpj\nxsDIkUVHYmZWeU4wBdppp/QxM6tFHiIr2EMPwezZRUdhZlZ5TjAFu/RS+MIX/MKlmdUeJ5iCNTSk\nRcjuvbfoSMzMKivXBCPpIEmPSGqSdEqZ40Mlzc2Oz5dUl5WPlHSzpJWSLiipP0zS/0j6u6QFks4u\nOfZ5Sc2S7ss+X8zz2ipln33St2dWNrNak1uCkTQQuBA4GJgAHC1pQptqxwEvR8RuwE+Ac7LyN4DT\ngZPLNP3DiNgd2AvYV9LBJcfmRsTE7HNxBS8nN9tvn270+4VLM6s1efZgJgFNEbEoIlYDc4CpbepM\nBWZl21cCUyQpIl6LiL+SEs1bImJVRNycba8G7gHG5ngNvWLyZLjrrqKjMDOrrDwTzA7A0yX7S7Ky\nsnUiYi2wHOjSWyGStgEOBW4sKT5S0gOSrpS0YzvnzZDUKKmxubm5a1eSsx//GB58sOgozMwqq8ME\nI2nrDo519gaHypS1nRSlK3XK/fYg4Arg/IhYlBX/EaiLiD2BG2jtGW3YeMTMiKiPiPrRo0d39lO9\nYvvtYYstio7CzKyyOuvB3NKyIenGNsf+0Mm5S4DSXsRYYGl7dbKkMRx4qZN2AWYCj0XEeS0FEfFi\nRLQ87HsRsHcX2ukzzjgDLrqo6CjMzCqnswRT2sPYtoNj5dwFjJe0s6QhwDRgXps684Bjs+2jgJsi\nOp76UdK/khLR19qUjynZPQxY2El8fco118BllxUdhZlZ5XQ2VUy0s11uf8ODEWslnQRcCwwEfhER\nCySdCTRGxDzgEuBSSU2knsu0lvMlLQa2BoZIOhw4EHgVOA34O3CPJIALsifG/knSYcDarK3Pd3Jt\nfUpDQ1pGec0aGDy46GjMzDZdZwlmO0nfIPVWWrbJ9ju9gRERfwL+1Kbsn0u23wA+1c65de00W7bn\nFBGnAqd2FlNf1dAAP/1putn/vvcVHY2Z2abrbIjsImArYMuS7Zb9qnjPpFpMnpy+/T6MmdWKDnsw\nEXFGbwXS3+20E7zrXbByZdGRmJlVRocJRtLxwC0R8ZjSDY9LgCOBJ4FjI8IzaFWIBAsWpG8zs1rQ\n2RDZV4HF2fbRwHuBXYBvAOfnF1b/5ORiZrWkswSzNiLWZNuHAL/O3je5AfCrgRXW1AR77JEeWTYz\nq3adJZj1ksZI2gyYQnpDvsXm+YXVP40ZAwsXwt/+VnQkZmabrrPHlP8ZaCS9xzIvIhYASPogsKij\nE637ttgC9tzTU/ebWW3oLMEsAyYDKyLiZUnHkG7yLwNm5B1cf9TQkN7oX7cOBg4sOhozs57rbIjs\nv4CVWXL5AHA28GtSgvlp3sH1Rw0NsGJFGiozM6tmnfVgBkZEy+STnwZmRsTvgN9Jui/f0Pqn/faD\n6dP9RJmZVb9OE4ykQdlaLVPYcFiss3OtB3bZxZNemllt6CxJXAHcKukF4HXgLwCSdiMtDmY5iICl\nS2GHtsuzmZlVkQ7vwUTEWcA3gV8B+5VMpT8A+Md8Q+u/fvhD2HFHWO4UbmZVrNNhrojYaPrFiHg0\nn3AMYOLE1Iu580444ICiozEz65nOniKzAkyalG7ye2ZlM6tmTjB90PDhMGGCE4yZVTcnmD6qoSEl\nmI4XkDYz67v8qHEfddxx6f7LunUwyP8rmVkV8p+uPmry5NZVLs3MqlGuQ2SSDpL0iKQmSaeUOT5U\n0tzs+HxJdVn5SEk3S1op6YI25+wt6cHsnPOzhdCQtK2k6yU9ln2PyPPaesN998EttxQdhZlZz+SW\nYCQNBC4EDgYmAEdLmtCm2nHAyxGxG/AT4Jys/A3gdODkMk3/nDSjwPjsc1BWfgpwY0SMB27M9qva\nt74FX/960VGYmfVMnj2YSUBTRCyKiNXAHGBqmzpTgVnZ9pXAFEmKiNci4q+kRPMWSWOArSPi9uyl\nz18Dh5dpa1ZJedVqaIAHHoDXXis6EjOz7sszwewAPF2yvyQrK1snm+9sOTCykzaXtNPm2yLi2ayt\nZ4HtyjUgaYakRkmNzc3NXbyUYjQ0wPr10NhYdCRmZt2XZ4IpNx9w24duu1JnU+pvXDliZkTUR0T9\n6NGju3Nqr9tnn/Tt92HMrBrlmWCWADuW7I8FlrZXR9IgYDjwEu1bkrVTrs1l2RBay1Da8z2OvI8Y\nNQrGj/cKl2ZWnfJMMHcB4yXtLGkIMA2Y16bOPODYbPso4KaSCTU3kg19rZDUkD09dgxwdZm2ji0p\nr2p/+AP8+tdFR2Fm1n25JZjsnspJwLXAQuA3EbFA0pmSDsuqXQKMlNQEfIOSJ78kLQZ+DHxe0pKS\nJ9C+DFwMNAGPA9dk5WcDB0h6DDgg2696994Le+4JAwZAXR3Mnl10RGZmXaMOOgw1r76+Phr78B30\n2bPh+OPh9ddby4YNg5kz06qXZmZFkHR3RNR3Vs9zkfVhp522YXIBWLUqlZuZ9XVOMH3YU091r9zM\nrC9xgunDdtqpe+VmZn2JE0wfdtZZ6Z5LqYEDU7mZWV/nBNOHTZ+ebuiPG5dWuNxuOzjvPN/gN7Pq\n4On6+7jp0zdOKOvWwd/+BvvvX0xMZmZd4R5MFfrRj+BDH4I//7noSMzM2ucEU4VOPBH22AM+/Wn4\n+9+LjsbMrDwnmCq0xRZw9dUwdCgcdhi8/HLREZmZbcwJpkqNGwdXXQWLF8MxxxQdjZnZxnyTv4rt\nuy/MmgW77FJ0JGZmG3OCqXJHH926/fjjsOuuxcViZlbKQ2Q1YtYseNe74NZbi47EzCxxgqkRhx+e\nei9HHgmLFhUdjZmZE0zNGD4c5s2D9eth6lRYsaLoiMysv3OCqSHjx8NvfwsLF6a3//vxUj9m1gf4\nJn+NmTIFLrgANt88zV9mZlYUJ5ga9KUvtW4vX56Gz8zMepuHyGrYdddBXR3Mn190JGbWH+WaYCQd\nJOkRSU2STilzfKikudnx+ZLqSo6dmpU/IuljWdk7Jd1X8nlV0teyY9+X9EzJsY/neW3V4H3vgxEj\n0hNmS5YUHY2Z9Te5JRhJA4ELgYOBCcDRkia0qXYc8HJE7Ab8BDgnO3cCMA14N3AQ8DNJAyPikYiY\nGBETgb2BVcBVJe39pOV4RPwpr2urFqNGwR//CCtXpiSzalXREZlZf5JnD2YS0BQRiyJiNTAHmNqm\nzlRgVrZ9JTBFkrLyORHxZkQ8ATRl7ZWaAjweEU/mdgU14N3vhssvh3vugS98wU+WmVnvyTPB7AA8\nXbK/JCsrWyci1gLLgZFdPHcacEWbspMkPSDpF5JGlAtK0gxJjZIam5ubu3M9VevQQ+Hss2GHHdJ7\nMmZmvSHPBFPuIdm2/35ur06H50oaAhwG/Lbk+M+BXYGJwLPAj8oFFREzI6I+IupHjx7dfvQ15tvf\nTguVDRyYVsQ0M8tbnglmCbBjyf5YYGl7dSQNAoYDL3Xh3IOBeyJiWUtBRCyLiHURsR64iI2H1Ay4\n/36YMCF9m5nlKc8EcxcwXtLOWY9jGjCvTZ15wLHZ9lHATRERWfm07CmznYHxwJ0l5x1Nm+ExSWNK\ndo8AHqrYldSQ7baD115LC5U9/3zR0ZhZLcstwWT3VE4CrgUWAr+JiAWSzpR0WFbtEmCkpCbgG8Ap\n2bkLgN8ADwN/Bk6MiHUAkoYBBwC/b/OT50p6UNIDwIeBr+d1bdVszJi0GmZzc5oY8803i47IzGqV\noh8/VlRfXx+NjY1Fh1GIuXNh2rT0ZNnFF3taGTPrOkl3R0R9Z/U8VUw/9elPw0MPwV13pV7MZpsV\nHZGZ1RonmH7sjDPSY8uD/F+BmeXAc5H1YwMGpOTy7LNw0EHwyCNFR2RmtcQJxli9Or3pf+ih8PLL\nRUdjZrXCCcYYNw6uugoWL4b990/7AwakmZhnzy46OjOrVk4wBsC++8LnPw8LFsBTT6U5y558EmbM\ncJIxs55xgrG3XHfdxmWrVsFpp/V+LGZW/Zxg7C1PPdW9cjOzjjjB2Ft22qn9Y6ecAi+91HuxmFn1\nc4Kxt5x1FgwbtmHZZpvB5Mlw7rmwyy7wb//mNWXMrGucYOwt06fDzJnpKTIpfV98Mdx2W5p9+YMf\nhDvvbJ1WxmvLmFlHPBdZP52LrKfefBOGDoVHH4VPfAK+9z347GfTOjNm1j90dS4y92CsW4YOTd8r\nV8Lw4enR5j33TO/R9ON/q5hZGU4w1iPve1+aKPO3v00rZH7yk2kIzatlmlkLJxjrMQmOOirNynzJ\nJWk+s5ahsoULi43NzIrnBGObbNCgtK7Md7+b9m+5JS3LfMQRaWYAM+ufnGCs4vbeG848E266CfbY\nA449Ns1zZmb9ixOMVdxWW8Hpp8OiRXDyyfCb36S5ztasKToyM+tNuSYYSQdJekRSk6RTyhwfKmlu\ndny+pLqSY6dm5Y9I+lhJ+WJJD0q6T1JjSfm2kq6X9Fj2PSLPa7POjRyZXtBsaoJZs2Dw4PQQwA9/\nCK+8kibRrKvzzM1mtSq392AkDQQeBQ4AlgB3AUdHxMMldb4C7BkRX5I0DTgiIj4taQJwBTAJ2B64\nAXhHRKyTtBioj4gX2vzeucBLEXF2lsxGRMR3OorR78H0vltugQ9/GDbfHNau3bBXM2xYetFz+vTC\nwjOzLugL78FMApoiYlFErAbmAFPb1JkKzMq2rwSmSFJWPici3oyIJ4CmrL2OlLY1Czi8AtdgFfah\nD8G996Z3ZtoOmXnmZrPakudq7DsAT5fsLwH2aa9ORKyVtBwYmZXf0ebcHbLtAK6TFMB/RcTMrPxt\nEfFs1tazkrar5MVY5UycmGYEKOfJJ+GQQ+A970mfPfaA3XdvfcHTzKpHnglGZcrajse1V6ejc/eN\niKVZArle0t8j4n+7HJQ0A5gBsFNH0wdbrnbaKSWTtrbaKi0PcN11rT2cKVPghhvS9s9+BmPGpOSz\nyy6eosasL8tziGwJsGPJ/lhgaXt1JA0ChgMvdXRuRLR8Pw9cRevQ2TJJY7K2xgDPlwsqImZGRH1E\n1I8ePbrHF2ebptzMzcOGwc9/Dg88AK+9lt6hmTMHvv71dHzNmrT9yU/CO94BW26ZHomemfVhI+CZ\nZ9qfssYPFZj1rjx7MHcB4yXtDDwDTAM+06bOPOBY4HbgKOCmiAhJ84DLJf2YdJN/PHCnpC2AARGx\nIts+EDizTVtnZ99X53httolabuSfdlrqsey0U0o6LeWDB6eXNSdMaD1n8OC0Js3DD6fZAx58MH23\nJJSlS2HsWBgxYsMhto99DG6/PS3/vGpVqtuyHHRpLGZWYRGR2wf4OOlJsseB07KyM4HDsu3NgN+S\nbuLfCexScu5p2XmPAAdnZbsA92efBS1tZsdGAjcCj2Xf23YW39577x1WO158MeKCCyJOOCFi330j\nhg+PgIjLLosYNy5tt/2MGhWxbFk6/403Il5/vXLxtPyulL4vu6xybZsVCWiMLuQAT9fvx5RrVsuQ\n2dZbwzbbtD909pe/wH77wWWXwec+l4beRo2C0aPT9wUXpPs999+f1sNpKW/5HjEiDbuVmj17wx4T\n+DFsqx1dfUw5zyEys0JJacgM2n+oYMyYNDM0wHvfm4bpXngBmptbvwdl/1/y5z+npaPbWrIEdtgB\n/vM/4fLLU9K5/voNkwu0PobtBGP9hROM9QtnnVW+R/GDH7Q+bLDHHunTnq9+FT7zmZR4SpNQy7Mi\nQ4aknsyjj6b1csp58sl0P6i+Pt1TMqtlHiLzEFm/MXt2+w8VVFpdXfkeU4sttkjDclOnwpe/nE8M\nZnnpC2/ym/Up06enWZ3Xr0/feQ5VdfQY9m9/m1YCffppuPHG1uNf/GLqUTU2euE2qw0eIjPLQWeP\nYR91VPpueZl05Uq47ba0cBukBxM+8AE48cS0kJtZNXIPxiwnXekxtdyH2XLLtAro0qVwxRUwbVq6\nl7NsWTr+6KNpOO288+C++1KbpfwSqfVF7sGY9SFjxqTkMm1a2i99iXTBApg3L+1vuy188INp6QO/\nRGp9lW/y+ya/VZGnn05LHtx8M9x6a3ovZ++9yz9QsOOOaXjOrNK6epPfCcYJxqrcgAHtv0S6fn16\nH+g//iMlobq61s+4cWlyUbPu8ouWZv1Eey+RbrttSi4A8+fD734Hb7zRenz33dN9H4B/+Zc0wWhp\nAqqrS/eGyunNR76tejnBmFW59l4iPf/81v3LLoNLL4Xnn08PHCxe3Jp8IN3H+ctfNkxABx4I116b\ntk88Ma1CWleXzr3wwta6vudj7XGCMatynT0S3UKCt70tffZps/TfddelYbbSBDR8eDoWAX/9a3qS\nrTQBlVq1Kr3H86c/palyRo5MS2Pvv396p+f++1vLhw3bMLl1hXtM1ckJxqwGTJ++6X9w20tAUkoQ\nLQlozJjy93zeeAPuuCNNn/MZJya5AAAJ+0lEQVTqq6nO/vunc/beu7Xe0KEp0ZxxRkpKy5alIbqR\nI1uT0MiRaY64MWNS7+uEE/yUXDVygjGzLmlJQO3d8xk3Dh5/PG2vWdM6G8Hw4XDVVSnxvPhi66dl\nQdnmZvj979NaP6UzGMyaBcccAyefXH7i0JNPTgnmiSfgmmvSPaeRI9P3ttumCUiHDOn+dfZ2b6mW\ne2dOMGbWLe3d8znrrNb9wYNbXyIdNgwOP7z99t7zntTLWb8+9XxaElBdXTr+fNm1aVtfQm1sTPeI\n2mpZhuF3v4Pvfa818bQkoW99K/WQFi1Kn223Teeceiq8/npqI+/eUttlHWqtd+bHlP2Yslm39YWJ\nQ8eNS/eK1qxJvZ8XX0zfLZ9DDklDbjfdlJZSKK3z4otp2G/XXeHcc+E73+k4hoED03tFAwak7YED\n4d57YbPN0suuc+a0lg8cmHpOLfPMnXdeWr6h5diAAWkqoF/+sv1rGzky9fgg9dAi0ppGw4enNjZF\nJf6382PKZpabStzz6arOekyDB7feOyrnIx9Jn/YccwxMnpySzhFHlK+zbl2aG27duvRZv771D/1W\nW6Xfbjm2bt2GDzG8+mrqbbWct25d6/tH7b0I++KLrdvHH7/hpKjDh8OkSenBDIDvfheeey4loBEj\n0vduu8HBB6fjjz+eEuE228Af/tC7PSb3YNyDMevzeqvH1Flvqbd+b+zYNGsDpBkbnngCXnkFXn45\nfY8aBaefno4fcUQaJnzlldZ1iKZMgRtuSNu77pqGADvS3evzm/xd4ARjZqV6e6nrSv/emjWwfDms\nXQtvf3squ/rq1IN6+eXyK7JC6nG1nUC1I31iPRhJB0l6RFKTpI0uTdJQSXOz4/Ml1ZUcOzUrf0TS\nx7KyHSXdLGmhpAWSvlpS//uSnpF0X/b5eJ7XZma1Z/r09Md93Lj0R3fcuPySSx6/N3hw6t20JBdI\ns3DPmJHuM40bV/68lif6Ki23HoykgcCjwAHAEuAu4OiIeLikzleAPSPiS5KmAUdExKclTQCuACYB\n2wM3AO8AtgPGRMQ9krYC7gYOj4iHJX0fWBkRP+xqjO7BmFl/UqkeU1/owUwCmiJiUUSsBuYAU9vU\nmQrMyravBKZIUlY+JyLejIgngCZgUkQ8GxH3AETECmAhsEOO12BmVjN6u4eWZ4LZAXi6ZH8JGyeD\nt+pExFpgOTCyK+dmw2l7AfNLik+S9ICkX0gaUS4oSTMkNUpqbG5u7u41mZlVtd5cOjzPBFNutqG2\n43Ht1enwXElbAr8DvhYRr2bFPwd2BSYCzwI/KhdURMyMiPqIqB89enTHV2BmZj2WZ4JZAuxYsj8W\nWNpeHUmDgOHASx2dK2kwKbnMjojft1SIiGURsS4i1gMXkYbozMysIHkmmLuA8ZJ2ljQEmAbMa1Nn\nHnBstn0UcFOkpw7mAdOyp8x2BsYDd2b3Zy4BFkbEj0sbkjSmZPcI4KGKX5GZmXVZbm/yR8RaSScB\n1wIDgV9ExAJJZwKNETGPlCwuldRE6rlMy85dIOk3wMPAWuDEiFgnaT/gc8CDku7Lfuq7EfEn4FxJ\nE0lDaYuBE/K6NjMz65xftPRjymZm3eI3+btAUjNQZqKGPmkU8ELRQeSklq8Navv6fG3Va1Oub1xE\ndPqUVL9OMNVEUmNX/sVQjWr52qC2r8/XVr164/pynSrGzMz6LycYMzPLhRNM9ZhZdAA5quVrg9q+\nPl9b9cr9+nwPxszMcuEejJmZ5cIJxszMcuEE08d1tMharZA0UNK9kv676FgqSdI2kq6U9Pfsf7/J\nRcdUKZK+nv33+JCkKyRtVnRMmyKbgf15SQ+VlG0r6XpJj2XfZWdo7+vaubYfZP9dPiDpKknb5PHb\nTjB931rgmxHxLqABODFbkK2WfJW0tk+t+Snw54jYHXgvNXKNknYA/gmoj4j3kKaCmlZsVJvsV8BB\nbcpOAW6MiPHAjdl+NfoVG1/b9cB7ImJP0sKQp+bxw04wfVytL7ImaSzwCeDiomOpJElbAx8gzbdH\nRKyOiFeKjaqiBgGbZ7OgD2PjmdKrSkT8L2k+xFKlCyLOAg7v1aAqpNy1RcR12RpcAHeQZqyvOCeY\nKtLOImvV7jzg28D6ogOpsF2AZuCX2fDfxZK2KDqoSoiIZ4AfAk+R1l5aHhHXFRtVLt4WEc9C+oce\nacn2WvQF4Jo8GnaCqRLtLLJW1SQdAjwfEXcXHUsOBgHvA34eEXsBr1G9QywbyO5FTAV2BrYHtpD0\n2WKjsp6QdBppGH52Hu07wVSB9hZZqwH7AodJWgzMAT4i6bJiQ6qYJcCSiGjpbV5JSji14KPAExHR\nHBFrgN8D/6fgmPKwrGWdqez7+YLjqShJxwKHANMjpxcinWD6uI4WWat2EXFqRIyNiDrSTeKbIqIm\n/iUcEc8BT0t6Z1Y0hbS+US14CmiQNCz773MKNfIAQxulCyIeC1xdYCwVJekg4DvAYRGxKq/fcYLp\n+/YlLbL2EUn3ZZ+PFx2Udck/ArMlPQBMBP6t4HgqIuuVXQncAzxI+jtS1dOqSLoCuB14p6Qlko4D\nzgYOkPQYcEC2X3XaubYLgK2A67O/Kf+Zy297qhgzM8uDezBmZpYLJxgzM8uFE4yZmeXCCcbMzHLh\nBGNmZrlwgjHrAUl1pbPTdvPcX0k6qoe/+Zme/KZZEZxgzKpHHdCtBCNpYD6hmHXOCcZsE0naJZvQ\n8v1ljn1b0oOS7pe00Yt6khZLGpVt10u6Jdv+YMmLtfdK2or0ot/+WdnXs3V0fiDprmxdjxOycz+U\nrSF0OfCgpC0k/U8Ww0OSPp3n/z3MWgwqOgCzapZNBTMH+L8RcV+bYweTpnjfJyJWSdq2G02fDJwY\nEbdlE52+QZos8+SIOCRrfwZpJuP3SxoK3CapZVbjSaT1Pp6QdCSwNCI+kZ03vOdXbNZ17sGY9dxo\n0vxUn22bXDIfBX7ZMtdTRLRdb6QjtwE/lvRPwDYla3eUOhA4RtJ9pCUcRgLjs2N3RsQT2faDwEcl\nnSNp/4hY3o04zHrMCcas55YDT5PmiytHQGdzMa2l9f8P31p2OCLOBr4IbA7cIWn3dtr/x4iYmH12\nLlmX5bWSth4F9iYlmn+X9M+dxGRWEU4wZj23mjQEdkw7T3ddB3xB0jBIa7yXqbOY9Mcf4MiWQkm7\nRsSDEXEO0AjsDqwgTVDY4lrgy9lyDkh6R7lFzSRtD6yKiMtIC4XVyrIB1sf5HozZJoiI17KF066X\n9FpEXF1y7M+SJgKNklYDfwK+26aJM4BLJH2XDVcq/ZqkDwPrSNP8X0Na9XOtpPtJ66z/lPRk2T3Z\ntPnNlF/Wdw/gB5LWA2uAL2/iZZt1iWdTNjOzXHiIzMzMcuEEY2ZmuXCCMTOzXDjBmJlZLpxgzMws\nF04wZmaWCycYMzPLxf8H9mSydvgA/5cAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"SSE = []\n",
"coords = list(cluster_df[['type_A_norm','longitude_norm','latitude_norm']].values)\n",
"for k in range(1,13):\n",
" km = KMeans(n_clusters=k)\n",
" km.fit(coords)\n",
" labels = km.labels_\n",
" centroids = km.cluster_centers_\n",
" # Get the SSE\n",
" SSE.append(sum(np.min(cdist(coords, centroids, 'euclidean'),axis=1))/len(coords))\n",
" \n",
"plt.plot(range(1,13),SSE, 'bo--')\n",
"plt.xlabel('k clusters')\n",
"plt.ylabel('SSE')\n",
"plt.title('SSE vs k')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"4 seems like a decent choice, but interestingly so is 7. Let's try 7."
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAAHiCAYAAAD8hSV1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzsvX+cHFWZ7/8+3TPToQkE6JCBGKfa\nuUqIiCJEIVcur+ioCyirq16FbRAC2JuJXsN63b177b0L6Pb+8LoLcTVgrxtUur7hu351V1GCwjAJ\nSMbF+Iu4hgEcp4cQpyETCcKYyWT6fP+o7kz/qOqu7q7+NfO8X696zfSp6lOnqrvrc85znvM8SmuN\nIAiCIAidi6/VDRAEQRAEoT5EzAVBEAShwxExFwRBEIQOR8RcEARBEDocEXNBEARB6HBEzAVBEASh\nwxExF4QGopS6RSmVbHU7qkUpdZ1S6getbocTSqlPKaW+3Op2CEK7IGIuCHWilPpjpdQepdRLSqnf\nKKV2KKUu9rD+sFJKK6W6vKqzHprdHqXUeqXU/vwyrfXfaK1vbMb5BaETEDEXhDpQSn0CuB34G6AX\n6AO2Au9pZbvyaZdOgCAIjUPEXBBqRCm1DPg08FGt9Te11i9rrWe11vdqrf/M5viSEaZSalwp9fbs\n/2/OjvBfVEqllVL/mD3s4ezfF7Kj/3XZ469XSu1TSv1WKfU9pZSRV69WSn1UKfUU8FS27Gyl1ANK\nqUNKqVGl1Afzjg8ppb6dPfdjwH+p8Z4ElFK3K6UOZLfblVKBvP3vUUr9LHueXymlLs2Wb8hey++U\nUmNKqT/Jlp8I7ABWZq/9JaXUyuLpC6XUHyql/lMp9YJSaqdSak3RPf6kUupxpdRhpdT/q5Rakt23\nXCn1nez7DimlHlFKyXNR6DjkSysItbMOWAL8m0f1bQG2aK1PxhLTf82WX5L9e4rWeqnWekQp9V7g\nU8D7gNOBR4DtRfW9F7gQeG1WFB8A/h9gBXAVsFUpdU722C8CR4AzgeuzWy3EgIuA84A3AG8G/hKs\nzgrwNeDPgFOy1zWefd9zwLuBk4ENwG1KqfO11i8DlwEHste+VGt9IP+ESqmzstd+U/Ze3Afcq5Tq\nyTvsg8ClwKuA1wPXZcv/J7A/+75erHsqMa6FjkPEXBBqJwQc1Fof86i+WeDVSqnlWuuXtNY/LHPs\nnwB/q7Xelz3/3wDn5Y/Os/sPaa1/jyWU41rru7TWx7TWPwG+AXxAKeUH3g/8Vda68AvgqzVeQwT4\ntNb6Oa3188CtwDXZfTcA27TWD2itM1rrZ7XWTwBorb+rtf6VttgFfB/4by7P+SHgu9l6Z4HPAScA\n/zXvmM9rrQ9orQ8B92J1NsC652cCRtaq8oiWhBVCByJiLgi1MwUs93BO+gbgLOAJpdSPlFLvLnOs\nAWzJmodfAA4BCnhF3jHPFB1/Ye747HsiwBlYo9KuouNTNV7DyqL3prJlAK8EfmX3JqXUZUqpH2ZN\n3S8AlwPLazmn1jqDdS3592Iy7/9pYGn2//8LPA18P2ve/wuX5xSEtkLEXBBqZwTLNP1el8e/DARz\nL7Ij4tNzr7XWT2mtr8Iyg/898P9lzeN2I8VngD/RWp+St52gtd6dd4wuOn5X0fFLtdaDwPPAMSyx\nzdHn8pqKOYDVccivJ2cWfwabufjsnPo3sEbUvVrrU7BM5crmOiqeUymlsK7l2UqN1Vr/Tmv9P7XW\n/cAVwCeUUgOV3icI7YaIuSDUiNb6MPBXwBeVUu9VSgWVUt3ZUeZnbd7yJLBEKfUupVQ31lxyvnPY\n1Uqp07MjyxeyxXNYYpsB+vPquhP437k5b6XUMqXUfy/T3O8AZymlrsm2sVsp9Sal1Bqt9RzwTeCW\n7DW8FrjWxS0IKKWW5G0+rLnrv1RKna6UWp69PzlHtX8BNiilBpRSPqXUK5RSZwM92fvwPHBMKXUZ\n8M6886SBkLIcDu34V+Bd2Xq7sebBZ4DdDscfRyn1bqXUq7MdgBex7veci2sXhLZCxFwQ6kBr/Y/A\nJ7CE+Xms0efHgH+3OfYwsAn4Mtao8WUs56sclwL/qZR6CcsZ7kqt9RGt9TQQBx7Nmsgv0lr/G9bo\n/R6l1IvAL7AcxZza+TssgbwSayQ7mX1/rjPxMSzT8yTwFeAuF5f/EvD7vO1twF8De4DHgb3AT7Jl\naK0fI+vcBhwGdmHNVf8O+DiWKP8W+GPg23ltfwKrkzCWvf6c2T63fxS4Gvgn4CDWCPsKrfVRF9fw\nGuDB7LWMAFu11jtdvE8Q2golvh6CIAiC0NnIyFwQBEEQOhwRc0EQBEHocETMBUEQBKHDETEXBEEQ\nhA5HxFwQBEEQOpyOyqa0fPlyHQ6HW90MQRAEQWgKP/7xjw9qrU+vdFxHiXk4HGbPnj2tboYgCIIg\nNAWllKvQymJmFwRBEIQOR8RcEARBEDocEXNBEARB6HA6as5cEARBEKpldnaW/fv3c+TIkVY3xZEl\nS5awatUquru7a3q/iLkgCIKwoNm/fz8nnXQS4XAYK0Fee6G1Zmpqiv379/OqV72qpjrEzC4IgiAs\naI4cOUIoFGpLIQdQShEKheqyHIiYC4IgCAuedhXyHPW2T8RcEARBEBrM/fffz+rVq3n1q1/N3/3d\n33levysxV0ptU0o9p5T6hcN+pZT6vFLqaaXU40qp8/P2XauUeiq7XZtXfpVSam/2+PuVUsvrvxxB\nEARBqBPThHAYfD7rr2nWVd3c3Bwf/ehH2bFjB7/85S/Zvn07v/zlLz1pag63I/OvAJeW2X8Z8Jrs\nFgXuAFBKnQbcDFwIvBm4WSl1qlKqC9gCvFVr/XrgceBjtVyAIAiCIHiGaUI0CqkUaG39jUbrEvTH\nHnuMV7/61fT399PT08OVV17Jt771LQ8b7VLMtdYPA4fKHPIe4Gva4ofAKUqpM4E/AB7QWh/SWv8W\neACrU6Cy24nKmig4GThQx3UIgiAIQv3EYjA9XVg2PW2V18izzz7LK1/5yuOvV61axbPPPltzfXZ4\ntTTtFcAzea/3Z8tsy7XWs0qpQWAv8DLwFPBRj9oiCIIgCLUxMVFduQu01iVlXjvkeeUAZ9cq7VSu\nlOoGBoE3AiuxzOz/27ZipaJKqT1KqT3PP/+8R80VBEEQBBv6+qord8GqVat45pn5ce3+/ftZuXJl\nzfXZ4ZWY7wdemfd6FZbZ3Kn8PACt9a+01WX5V+C/2lWstU5orddqrdeefnrFLHCCIAiCUDvxOASD\nhWXBoFVeI29605t46qmn+PWvf83Ro0e55557+MM//MM6G1qIV2L+beDDWa/2i4DDWuvfAN8D3pl1\nejsVeGe27FngtUqpnDq/A9jnUVsEQRAEoTYiEUgkwDBAKetvImGV10hXVxdf+MIX+IM/+APWrFnD\nBz/4Qc455xwPG+1yzlwptR1YDyxXSu3H8lDvBtBa3wncB1wOPA1MAxuy+w4ppT4D/Chb1ae11oey\ndd4KPKyUmgVSwHXeXJIgLAxMEzZvhqkp63UoBFu21PVMEQTBDZGI5z+0yy+/nMsvv9zTOvNxJeZa\n66sq7Nc4OLBprbcB22zK7wTudHN+QVhMmOk0m/eOMXXmDHwhAF/uh6FepqZgwwbrGBF0QRDykQhw\ngtBGmOk00dFRprpmrF/nGTPwyVEYSAMwO1vXChlBEBYoIuaC0EbExsaYzmQKC5dk4C/2HRf0OlbI\nCIKwQJEUqILQRkzMzNjv6MIaoQN9T/c2r0GCIHQEIuaC0Eb0BQKknAR9SQY+Mkb8mIi5IAiFiJld\nENqIeH8/QV+Zn+WKGXF+EwShBBFzQWgjIr29JFavxu+w31gSaGp7BEHwhuuvv54VK1bwute9riH1\ni5gLQpsR6e3lq2vWlIzQFXB5KNSaRgnCIsLcaxK+PYzvVh/h28OYe+tLgQpw3XXXcf/993vQOntE\nzAWhDYn09nLtGWcUJDfQwFcnJzHT6VY1SxAWPOZek+i9UVKHU2g0qcMpovdG6xb0Sy65hNNOO82j\nVpYiYi4Ibcp9U1MU51qazmSIjY21pD2CsBiIDcWYni1MgTo9O01sqL0DPIiYC0Kb4rRMzXH5miAI\ndTNx2D6Qg1N5uyBiLghtSl/A3tnNqVwQhPrpW2af6tSpvF0QMReENsVumVrQ5yPe39+iFgnCwic+\nECfYXZgCNdgdJD5QewrUZiBiLghtSm6ZmhEIoAAjECCxejWR3sYHjTFNCIfB57P+mvU78wpCRxA5\nN0LiigTGMgOFwlhmkLgiQeTc+gI8XHXVVaxbt47R0VFWrVrFv/zLv3jUYgtlJTzrDNauXav37NnT\n6mYIwoLGNCEahek8H6BgsO6Uzm2PudckNhRj4vAEfcv6iA/E636AC+3Bvn37WLNmTaubURG7diql\nfqy1XlvpvTIyF4QOIZ02GRkJs3Onj5GRMOl0Y4bLsVihkIP1eiFna2vUciRBaBYi5oLQAaTTJk88\ncT0zMylAMzOT4oknrm+IoDtlZVvI2do6dTmSIOQQMReEDuCppzaj9dGCMq2P8tRTmz0/V5+D065T\n+UKgU5cjCUIOEXNB6ACOHZuqqjxHLY5s8bg1R55PMGiVL1Q6dTmSIOQQMReEBUrOkS2VAq2tv9Fo\nZUGPRCxnN8MApay/C935rVOXIwlCDhFzQehwnEbf9TiyRSIwPg6ZjPV3IQs5zC9Hes8rTmT7hTB0\nCdx1wRFOOPJoq5smCK4QMReEDsDvt8+WNjcXchx9L0ZHtno44cijbHzVy5yxBHwKVgQyBF+8g2/+\naFOrmyYsAJ555hne+ta3smbNGs455xy2bNniaf0i5oLQAZx11hagu6i0m3/+5y2Oo+/F6MhWD5nf\nJlhSlEh+id8qFxYXjQia1NXVxT/8wz+wb98+fvjDH/LFL36RX/7yl/VXnEXEXBA6gN7eCGvW3EVX\n1/wI3e8/mYMH7Y+fmFicjmz1cFr3XFXlwsKkVl+TSpx55pmcf/75AJx00kmsWbOGZ5991oMWW4iY\nC0IHkcn8/vj/c3NT/NmfRRkYKH3K9PUtTke2ejg066+qXFiYNCNo0vj4OD/96U+58MILPatTxFwQ\nOoSxsRiZTOFTJhCY5iMfKXzK5I++F5sjWz34To1ypGgQfmTOKhcWD432NXnppZd4//vfz+23387J\nJ5/sTaWImAtCxzAzY/80WbFiQkbfHvC+N21l+uRBDh71k9Fw8Kif6ZMHed+btra6aUITaaSvyezs\nLO9///uJRCK8733vq7/CPLo8rU0QhIYRCPRlw7kWsmRJH+PjzW/PQsQSbhHvxUw8bp9oqF5fE601\nN9xwA2vWrOETn/hEfZXZICNzQegQ+vvj+HyFHm0+X5D+fvFoEwSvaJSvyaOPPsrdd9/NQw89xHnn\nncd5553Hfffd502jkZG5IHQMvb3W02RsLMbMzASBQB/9/fHj5YIgeEMk4v1U1cUXX0wjU46LmAtC\nB9HbGxHxFgShBDGzC4IgCEKHI2IuCIIgCB2OiLkgCIIgdDgi5oIgeIJpmoTDYXw+H+FwGNOLgNaC\nILhCHOAEQagb0zSJRqNMZxfnplIpolErclpEItgIQsORkblQFhltLV5M02T58uUopVBKsXz5csfP\nPxaLHRfyHNPT08S8DGgtCB3MkSNHePOb38wb3vAGzjnnHG6++WZP65eRueCIjLYWL6ZpsmHDBmZn\nZ4+XTU1Ncd111wGln/+EQ+Bqp3JBaGfSadPzeA6BQICHHnqIpUuXMjs7y8UXX8xll13GRRdd5Emb\nZWQuOCKjrcVLLBYrEPIcx44dY/PmzSXlfQ6Bq53KBaFdSadNRkej2dDJmpmZFKOjUdLp+qySSimW\nLl0KWDHaZ2dnUUp50GILEXPBERltLV7KfcZTU1MlZfF4nGBR8vRgMEhckqcLHYZddsJMZpqxsfoH\nMXNzc5x33nmsWLGCd7zjHZICVWgOTqOqV7zitCa3RLDDNCEcBp/P+uulO0O1I+pIJEIikcAwDJRS\nGIZBIpGQ6Rih43DKTuhUXg1+v5+f/exn7N+/n8cee4xf/OIXddeZQ8RccCQej3PCCT0FZYEAXHfd\ni3WbnIT6ME0rs1MqBVpbf6NR7wS93Ig6FArZlkciEcbHx8lkMoyPj4uQCx1JIGDfkXUqr4VTTjmF\n9evXc//993tWp4i54EgkEuHP//wkenut7EG9vfDJT8LAwKwnJiehdmKxwhSNYL32yp0hEokwODhY\nUt7d3c2WLVu8OYkgtCGNyk74/PPP88ILLwDw+9//ngcffJCzzz67rjrzEW92oSzr1x9i/frSci9M\nTkLtOE1pe+nOsHXrVt7ylrcQi8WYmJigr6+PeDy+aEbcjfBoFtqfRmUn/M1vfsO1117L3NwcmUyG\nD37wg7z73e/2osmAiLlQgUCgL+vVWVoutI6+Psu0blfuJZFIZNGIdz45j+acI1TOoxkQQV8ENCI7\n4etf/3p++tOfelpnPhXN7EqpbUqp55RStjP1yuLzSqmnlVKPK6XOz9t3rVLqqex2bV55j1IqoZR6\nUin1hFLq/d5cjuA1jTI5CfURj0OR8zjBoFUu1E8jPZoFoRG4mTP/CnBpmf2XAa/JblHgDgCl1GnA\nzcCFwJuBm5VSp2bfEwOe01qfBbwW2FVL44XG09sbYfXqBIGAASgCAYPVqxMyOmkxkQgkEmAYlj+D\nYVivF+EguiE00qNZEBpBRTO71vphpVS4zCHvAb6mtdbAD5VSpyilzgTWAw9orQ8BKKUewOoUbAeu\nB87O1p8BDtZxDUKDaYTJSaifSETEu1HI9JLQaXjhzf4K4Jm81/uzZbblSqlTsq8/o5T6iVLq60qp\nXg/aIQiLkkauN1+syPSS0Gl4IeZ28eh0mfIuYBXwqNb6fGAE+Jxj5UpFlVJ7lFJ7nn/+eQ+aKwjN\nwUynCY+M4Nu5k/DICGY67f05GrzefLEi00tCp+GFN/t+4JV5r1cBB7Ll64vKdwJTwDTwb9nyrwM3\nOFWutU4ACYC1a9dqD9orCA3HTKeJjo4ynckAkJqZITo6CkCk1ztDVLn15mKCrw+ZXhI6CS9G5t8G\nPpz1ar8IOKy1/g3wPeCdSqlTs45v7wS+l51bv5d5oR8AfulBOwShbYiNjR0X8hzTmQyxsTFPz9OM\n9eaCIHjD3Nwcb3zjGz1dX56j4shcKbUdS3iXK6X2Y3modwNore8E7gMuB57GGnFvyO47pJT6DPCj\nbFWfzjnDAf8LuFspdTvwfO49grBQmJiZqaq8Vk47DWzynnCahM8XhJox02liY2NMzMzQFwgQ7+/3\nxKK2ZcsW1qxZw4svvuhBKwtx481+VYX9Gviow75twDab8hRwics2CkLH0RcIkLIR7r5AoAWtEQTB\nLY2aItu/fz/f/e53icVi/OM//qMnbc1HYrMLQgOI9/cT9BX+vII+H/H+fk/Pc+hQdeWCIJSnUVNk\nN910E5/97Gfx+RojuyLmgtAAIr29JFavxggEUIARCJBYvdpT5zdwDt/qdVhXQVgsNGKK7Dvf+Q4r\nVqzgggsuqLmOSkhsdkFoEJHeXs/Fu5h43FqKlu/RLmFdBaF2GjFF9uijj/Ltb3+b++67jyNHjvDi\niy9y9dVXk0wm62lqATIyF4QORsK6Ng9zr0n49jC+W32Ebw9j7pXF/AuRRkyR/e3f/i379+9nfHyc\ne+65h7e97W2eCjmImAtCxxOJwPg4ZDLWXxFy7zH3mkTvjZI6nEKjSR1OEb036rmgSzS/1tOsKTKv\nUZYzemewdu1avWfPnlY3QxCERUb49jCpw6Wx2o1lBuM3jXtyjlw0v+IpE7G01M++fftYs2ZNq5tR\nEbt2KqV+rLVeW+m9MjIXBEGowMRh+yg8TuW1UC6anyBUQsRcEAShAn3L7JcHOJVD9XPsEs1PqAcR\nc0EQhArEB+IEuwuzqAW7g8QH7JcNbPruJq755jVVzbHLMkOhHkTMBUEQKhA5N0LiigTGMgOFwlhm\nkLgiQeTc0slsc6/JnXvuRFPojzQ9O01syNlmHo9bc+T5yDJDwS2yzlwQBMEFkXMjtuJdTGwoViLk\nOcrNseec3GIxy7Te12cJuTi/CW4QMRcEQfCQcoJdbo4dLOEW8RZqQcRcEATBQ/qW9dkuY1Moxzl2\nYXEQDoc56aST8Pv9dHV14eVSa5kzFwRB8BA7ZzmFYuPaja7M9ELrSZtpRsIj7PTtZCQ8QtpMe1b3\n8PAwP/vZzzwVcpCRuSAIgqfkBDs2FGPi8AR9y/qID8RFyDuEtJlmNDpKZtrKnDaTmmE0aqVA7Y20\nbxQ4GZkLQpsiscA7l8i5EcZvGidzc4bxm8ZFyDuIsdjYcSHPkZnOMBarLwUqgFKKd77znVxwwQUk\nEom668tHRuaC0IbkYoFPz1ohwXLrlAERBkFoIDMT9qlOncqr4dFHH2XlypU899xzvOMd7+Dss8/m\nkksuqbtekJG5ILQlsaHYcSHPUWmdcrtjmibhcBifz0c4HMaULCJCGxLos0916lReDStXrgRgxYoV\n/NEf/RGPPfZY3XXmEDEXhDakGbHAm4lpmkSjUVKpFFprUqkU0WhUBF1oO/rj/fiChdLoC/roj9ee\nAhXg5Zdf5ne/+93x/7///e/zute9rq468xExF4Q2pJZY4O1MLBZjuiiLyPT0NDHJIiK0Gb2RXlYn\nVhMwAqAgYARYnVhdt/NbOp3m4osv5g1veANvfvObede73sWll17qUatlzlwQ2pL4QLxgzhzKxwJv\ndyYcsoU4lQtCK+mN9Hruud7f38/Pf/5zT+vMR0bmgtCGVBML3A4znSY8MoJv507CIyOYae/WydZC\nn0O2EKdyQRCqQ0bmgtCmuI0FXoyZThMdHWU6Yy2vSc3MEB211slGeluzTjYejxONRgtM7cFgkLhk\nEREET5CRuSAsMGJjY8eFPMd0JkNsrP51srUSiURIJBIYhoFSCsMwSCQSRCQQuSB4gozMBWGBMTFj\nvx42NXOkyS0p5E1vepSvfGU/oIH9rFz5KCBiLgheICNzQVhg9AUc1sNqzaYffbO5jcny5JObOHDg\nDmAuWzLHgQN38OSTm1rSHkFYaIiYC0Kbk06bjIyE2bnTx8hImHS6/NrseH8/aJt82spH4reZ0vIm\ncOCAfehKp3JBEKpDxFwQ2ph02mR0NMrMTArQzMykGB2NlhV0y8nNRsyBue7TGtPQisxVWS4IC48X\nXniBD3zgA5x99tmsWbOGkZERz+oWMReENmZsLEYmUxhsJZOZZmysfLAV/+yhqsobj7/KckFoHY0K\nPbx582YuvfRSnnjiCX7+85+zZs0aT+oFEXNBaCqmabJ8eRilfCgVZvlyk3LPiZkZ+6AqTuU5oqf6\nYK7I4W3uiFXeAlaujFZVLgitolGhh1988UUefvhhbrjhBgB6eno45ZRTvGgyIGIuCE3DNE2uvz7K\n1JRlMocUU1NRNmxwFvRAwD6oilN5jq1veh+DJ0/jP3oQdAb/0YMMnjzN1je9r76LqJGzztrKypWD\nzI/E/axcOchZZ21tSXsEwYlGhR4eGxvj9NNPZ8OGDbzxjW/kxhtv5OWXX66rznyUtnOUaVPWrl2r\n9+zZ0+pmCEJNhMNhUqmUzR4DwxhnfLx0T27OPN/U7vMFWb06QW+vLOsSBDfs27fPtUnb5/Nhp4tK\nKTKZ2h1I9+zZw0UXXcSjjz7KhRdeyObNmzn55JP5zGc+U7adSqkfa63XVmx3zS0TBKEqnOOQT+C0\nq7c3wurVCQIBA1AEAoZrIW+3kK6C0Ak0KvTwqlWrWLVqFRdeeCEAH/jAB/jJT35SV535iJgLQpNw\nfhj0Ue450dsbYd26cdavz7Bu3bhrIY+OjpKambEM+tmQriLoglCeeDxOMBgsKPMi9PAZZ5zBK1/5\nSkazoZWHhoZ47WtfW1ed+YiYC0KTiMfj9PQEi0qDdHfH8TpEeTuGdBWETqCRoYf/6Z/+iUgkwutf\n/3p+9rOf8alPfcqDFltIOFdBaBK5h8HmzTGmpiaAPkKhOFu2RPA6RLlTSFenckEQ5olEIg3JG3De\neefRKL8vEXNBaCKNekgU0xcIkLIRbsdQr4IgdDRiZheEBUi8v5+gr/DnHfT5rFCvgiAsOETMBWEB\nEuntJbF6NUYggAKMQIDE6tUty2cuCEJjETO7ICxQIr29It6CkEVrjVKq1c1wpN6YLzIyFwRBEBY0\nS5YsYWpqqm7BbBRaa6ampliyZEnNdcjIXBAEQVjQrFq1iv379/P888+3uimOLFmyhFWrVtX8fhFz\nQRCajrnXJDYUY+LwBH3L+ogPxImcu/DC05rpNLGxMSZmZugLBIj398vURwvo7u7mVa96Vaub0VBE\nzAVBaCrmXpPovVGmZ61486nDKaL3WtnTFpKg56Lw5YL35KLwASLogue4mjNXSm1TSj2nlPqFw36l\nlPq8UupppdTjSqnz8/Zdq5R6Krtda/PebzvVKwjCwiM2FDsu5DmmZ6eJDdWXlaqV2MXBlyh8QjNx\n6wD3FeDSMvsvA16T3aLAHQBKqdOAm4ELgTcDNyulTs29SSn1PuClqlstCELHMnHYPquMU3m74xQH\n3y5oD0gUPqExuBJzrfXDwKEyh7wH+Jq2+CFwilLqTOAPgAe01oe01r8FHiDbKVBKLQU+Afx1PRcg\nCEJn0bfMISuVQ3m74zQC9zscL1H4hEbg1dK0VwDP5L3eny1zKgf4DPAPQKG9TRCEBU18IE6wuygr\nVXeQ+IDH2WbqwDQhHAafz/prms7HOo2050Ci8AlNwysxt1uJr53KlVLnAa/WWv9bxYqViiql9iil\n9rTzsgJBECpjmhC7IsL0X76Ef8sz8PgfYywzSFyRaBvnN9OEaBRSKdDa+huNOgu600g7F3VPovAJ\nzUC5XUSvlAoD39Fav85m35eAnVrr7dnXo8D63Ka1/pP844BTgP8DHMXyqF8B7NZary/XhrVr1+pG\nZZwRBKGx5ERyOs8WFwxCIoHnWePqIRy2BLwYw4Dx8dLyYq91sEbgItyCFyilfqy1XlvpOK9G5t8G\nPpz1ar8IOKy1/g3wPeCdSqlTs45v7wS+p7W+Q2u9UmsdBi4Gnqwk5ELnUI2JUlg8xGKFQg7W61ib\nObFPOPjhOZVLHHyhHXC1zlyePfYbAAAgAElEQVQptR1rlL1cKbUfy0O9G0BrfSdwH3A58DTWHPiG\n7L5DSqnPAD/KVvVprXU5RzqhwykefeVMlNBeoy/Be0zTEuaJCejrg3i88DOvViRbRV+f/ci8r4x/\nnsTBF1qNazN7OyBm9vanWhNlO1FJjARn3JjQO+W7YZqwYQPMzs6XdXfDXXfJ90FoPs02swsC0Dmj\nr2KqdXoSCnFjQo/HLYHPJxi0ytuN4uRabpJtmXtNwreH8d3qI3x7GHOvfHmE5iFiLniKkymynImy\nHfByPncx+gy46cRFItZI3TAscTSM9nN+A+szP3q0sOzo0fLfhVyI2tThFBpN6nCKa755DZu+u6mx\njRWqZqH+PsXMLnhKp3gsF+PzWSPyYpSConggZenU66+XTjGhu6GW70L49jCpw6U3QKG4+313t82y\nu8VOJ/4+xcwutIROGX0V45VFoVM8tr2mk0zolajlu+AUilajOzrm/EJjIf8+RczbENM0CYfD+Hw+\nwuEwZofZgSIRazSWyVh/213IwTsx6lSfgXrp1E6cHbV8F8qFou3UmPMLkYX8+xQxbzNM0yQajZJK\npdBak0qliEajHSfonYZXYtQJPgONmjPsxE6cHbnvQig0X3bCCeXfEx+Io2wDXnZuzPmFSCf8PmtF\nxLzNiMViTBfZgaanp4ktBDtQm+OFGDXT3GzuNVn+2eWoWxXqVsXyzy6v6EHdTl777e6I9Pvfz/8/\nNVX+PkXOjbBx7cYSQW+3mPOLnYU0HVSC1rpjtgsuuEAvdJRSGiuufcGmlGp10wSXJJNaG4bWSll/\nk8kGnOPxpO7+dLfmFgq2ns/06OTjzic0DK0tGS/cDMP7NpYjmdQ6GCxsQzDYmHtVC7Xep+TjSW3c\nZmh1i9LGbUbZz0JoDfX+PpOTk9rYvVur4WFt7N6tk5OTjWjmcYA92oU+ijd7mxEOh0nZuAUbhsF4\np7kFCw3DyXsawFhmMH7TuO0+r7z266Xdvd/b5T4J7UUr4vCLN3uHEo/HCRbZgYLBIPEFYQcSvKKc\nU1W5fe0yZ9jujkjtcp+E1mOm04RHRvDt3Mm1+/bZ5q6PjY21qHXziJi3GZFIhEQigWEYKKUwDINE\nIkGkU72JhIZQzqmq3L52mTNsd7Fsl/sktJbcSDw1M4PGylFvh1NO+2YiYt6GRCIRxsfHyWQyjI+P\ni5A3iXZ3yMonPhCn29ddUt7j7ynrcFWL1/6mTdY9UcraTjqp/nvT7mK5kJbaCeXJH3mHR0Yw0+nj\n+2JjYyUjcTuccto3ExFzQaC9vLzdEDk3wl3vvYvQCfPrp0InhNj2nm0Vo41V47W/aRPccUfh/PFL\nL8F119V3bzpBLOtZ3bBp0ya6urpQStHV1cWmTRLWtR0pHnmnZmaIjo4eF3Q3I+6gz0e8v7/BLa2M\nOMAJAu3vkNUqurpgzsG2uJjuTaWMeqZpEovFmJiYIBgM8vLLL5fUMTg4yNatW5vYaqES4ZERUjaC\nbQQCjK9b57jfD2SwRuTx/v6Gpr916wAnYi4IiPeyE+WyhS2We1Mpnncu0FNxfIhi/H4/x44da3Br\nhWrw7dyJnQIqILN+fUu810vaIt7sguCednfIahV+v/O+RtybdvRbqBTP2y7Qkx1zTiYOoWU4zXXn\nyiO9vSRWr8YIBFBYI/ZmCnk1iJgLAu4csjo9Zn4tRKP25V1d3jir5Yv38uVw/fXt57dQaRndhMv1\ndP5yPSOhJcT7+wn6CmWweA480tvL+Lp1ZNavZ3zdurYUckAiwAlCjnKRoZLJpA4GgwVR+YLBoE62\nS8iyBjI4aN2TXBS0pUu9idRmFwWuHaLTFVMpGpxhGLZRG4u3wcHBVl6G4ECzI7pVCy4jwLVcoKvZ\nRMyFVuH0wDZarTR10Iyws+VwEsnirdGRjCvdh0qhZ5PJpGMYZkD7/X4RcqFmRMwFwUOaGTO/KbHd\nmxAbvdJ15I/2WzUyL7Y6ON2HStcyODhY8h1ZLJabxUD+5/+h0KT+ztLd+iGG9XZ26w+FJhvaERYx\nFwQPadbIvFkJSBqdcMXNdbgZmTcy+Uoy6dyhqOU+JJNJbRiGVkppwzBEyBcIyaTWPT3W92KASb2D\nXXqY4ePbjsCwHvjYw3rwG40xz7sVc1maJggusFt+FAwGPQ+126z17sqnQdusO1ManSmzHs0lbq7D\nbslXdzecfDIcOmS/nttLnNoIi2fZnWBhptPExsaYmJkpWTu+fLmVAhdgOyOcQem688le+OOv+rj7\nPO893WVpmnCcxeiF7TWNiplfvBTLSVy8TkDiP+XZqspzbNpkebIrZf11CmzmJpGKXRS4u+6Cgwfr\nyynvlnL3dLEvSVxMbHrySa7Zt88xClxOyAFW2Ag5wIrnQAdam3BFxHyBkxtRplIptNakUimi0agI\neg14HTPfLoSsU5AWr8Vl7q3/C7qLopR1v2yVO5AL7ZpbLj03Z722E3S36/brCZlaL6ed5ryvXWLE\nC43FTKe588CBksAxTpnQnsN+XfpzK6y/rUy4ImK+wLELaDE9PU0sF/FCaBl2wUi0LhX0RiQgMf7b\no3DFR2DZOJCx/l7xEavcgUTCfXm7J1Ipx9Kl7RUjXmgcsbEx2whwMC/Mofn0B3yZfo4UyeaRAHz5\nRut/PRloWbAjEfMFjlNAC7eBLgR31BK5zOkj0LrxCUjiA3GCF3wL/vRVcIsf/vRVBC/4VtmMa04B\nzOzKOyGRyqFD9uU2YdVrxjStOddcxrnly1sfBEeYp9xIOhcFbssWy5eDgTRD28f43P/OMHk6ZJQ1\nV/65T8LQ24EjPvhyf8uCHYkD3AInHA6TspmINQyD8cWSJcNjTBM2b56fSzvxRDh6FGZn54/Jj93t\nRKuTu5h7TWJDMSYOT9C3rI/4QLxsxjWnpCt+P3gZcrxSUhOvaPT9N00rot3Ro4Xl3d2Wb0A7dWwW\nK06JVBRw95o1x53ZNn0zzZ0njqID816R3cDJXV1MzR6DdAC+3A9D885vXn2PxAFOACAejxMssncG\ng0HinWDvbENyD+h8p5iXXy4UciiM3e1Eq03RkXMjjN80TubmDOM3jVdMneoU2tWpvBaamYq2lvtf\njQUmFisVcrC+KzLL1R7YhXNVwMaVKwu80u87c6xAyAFmgaV+P2pgPVy1rkDIwXun1Yq4Wb/WLpus\nM68NWf/qHW6jlrmNXNbqKGzVMjiotd9vXZ/fb732kkavfy+mmvtfbQyAckFxGh3VbrFRz+/ITThX\nHhrWDJduanh4/js7MKnZvlszNKzZvluHPuTNunNknbkgVEfOvJtKWabjuTnLVJZv5nVKlWpHo83l\nzTJHe0mlNpe7v8lka6+vWrN8uaWGiykXfKOplKLWi/qv6RpB99rnPY+PrWODmWb246OwZH703pPx\nse2c+tedi5ldEKog37wL83PDxWZet0vEGm0ub6Y52u7ctaQpddPmcve3EddXzbW4WTufTzwOPT2l\n5X5/Z3j1dwqVUtR6Ub/+537LwS0PNWNlV4tE4ORPjBUIOcBRX5PXnbsZvrfLJmZ2oVFUMp/nzLz5\noR3zN79f61CoeeZyr83RyceT2rjN0OoWpY3bDJ183P4C3JianUyebtpcKZOal+b2as3mtdzzwcHS\n47u72386pZMoN53hxX0+Xn+RGZ2BeTO6sjHB58zw9YKY2QXBPZXM5/nhPYu92UMha/lKM03ATu2t\nJQypudckem+U6dm8ULXdQRJXJEqc4iqZmsuZPK+5xl2bTROuvtq+rV6GWa3WbF6LObfVKxYWA/nh\nVoup1txuptNs3jvGlH8G0gFC/94PD/ba1p//GTp5xRuBAOPr1rk7uQNiZheEKqhkPs/fH4lYIUdz\n/f+DB909LMy9JuHbw/hu9RG+PYy5t3absdsIa26IDcWY/vF74LZfwy1zcNuvmf7xe9j89/9RYoKu\nZGouZ/KsJiqcYbg7th6qNZvXsna+2nMI3lKNud1Mp7n+P0eZ6pqxXNrPmGHqulFeWJu21pnnUTyN\nZucVH/RZZvhmIWIuCNgvU8rhxfy3udfk+m9dT+pwCo0mdTjF1d+8GnWrqknYvVzWlnrkLXDvP8Ph\nMOCz/n7rLqa2f7ZkftspBGpOZMuJVzVtbsayvVo6RNWGn/Wy0yXY4xT8J4fbjlNsbIyjviKzz5IM\nc9eNcfLJ5Ttxkd5eEqtXYwQCKKwReWK190lXyiFiLhynVsemhUD+qAssJyXwLnLZ5h2bOTpns+gY\nSB1OEb03WpWgexlhzT/89zB7YmHhXAAySwqKciPuciJbTryqaXMzIsg1o8PQ6lgCi4FqrGrlcIwG\nt2KGQ4cqd+Iivb2Mr1tHZv16xteta6qQg0SAE7I0ennHYiKdNhkbizEzM0Eg0Ed/f5wz7nSYBM7D\nWGYwftN44xtYhGM6VLtjFdx9t/Pysk77HjVjeV8nLiFsJyr5qNh95xhIw41jsGKGUCbAlnP7K4qr\n07w3kwGMv1jXMh8Ht3PmIuYCII46XvHkk5s4cOBOyEvf4PMF+fQvphl6vvx7FYrMzc1Pol1uPXQx\nbr4PXoqXCOHixjRhw4bSCIs9PbBtW6Gg52JEsPlJ+MMDBXbnoM9X0eydmzMvMLUf8dH9+dXcFelt\n2fdOHOCEqhBHHfek0yYjI2F27vQxMhImnTaPlxcLOUAmM030v1T+qfUta81Eqp0puKeHik4/TniV\n1rSVa+kr4Tavu1AfsVipkIMVJjffsS33nUtOplHvPVCibE4pTfOJ9Pay7ZzVhI4FrJ/wZIDQV1or\n5NUgYi4A4qjjlnTa5IknrmdmJgVoZmZSPPHE9cdN68VCnuP0QIZuX7ftPrCWgpXLWNZIfwa7+elt\n26xkIK3MetboYCC1Uk1ed6E+8gcTA6TZzghD7GQ7I7w6lS453k1K03JEens5+PZ16LeuR1+5joP3\ndIaQg4i5kEUcdQoxTZNwOIzP5yMcDmNm1fOppzajdaEjm9ZHeeqxDzNzxNlWvSRgcNd778JYZqBQ\nhE4IETohhEJhLDNs13TPt6X8CNULobcbTXs1wq6VVluLzL0my6/+OOqUcZTKsHzlS5hmdXndhfrI\nDSYGSPNJRjmDGXzAGczwZ2qUtFko6G5Smi5UZM5cOI7MT1ps2rSJO++8k/zfRjAYJJFI8IpXODiy\naQg8p5jptfs9KdasuZve3tpuZjl/hni8sxzOqqGVfhzmXpMNn3mQ2X//QoGnf8+SYxw90uX4vg56\nnHYEuTnzr82OcAalQh0wAqwbnw/K4jalaSchc+ZC1bR6JNYOmKZZIuQA09PTxCrYd/v/WeObKfYK\nV6xcubFmIYfyI9R2NUV7QSutRbGhGLPfv7lkyV45Ic8tZ6yEk8/FYqWcZSkSsaZ7VtgIOcDMRGG5\n25SmCxERc0HIIxaLlQh5jomJCfz+kO0+/2HoHYLV/1cTCBiAIhAwWLPmbs46a2tdbSrnz9BqU3Qj\nacZacycmDk/A4eocRtzkdd+6dRNveMM1vOUtKa68UvPd76YYHY0uWkF34+QYicAJxryJ/EEe5Equ\n5G28jSt9Vx6fAgP74C13r1nD1rPOauJVtYhKwduBbcBzwC8c9ivg88DTwOPA+Xn7rgWeym7XZsuC\nwHeBJ4D/BP7OTRB5LYlWhCaglNJYXmwlm2EYenIyqYeHu/XwMPPb99CTA95kArHLPV8uIUiz838v\nFozbDM2yXzve21ryuieTSR0IFH6/AgF0LIbevdto9CW1JW6/v5PJSb0ruEvHiOkAgYJ7GAwGdXIB\nZ67BZaIVN2J+CXB+GTG/HNiRFfWLgP/Ilp8GjGX/npr9/9SsmL81e0wP8AhwmZvGipgLjcYwDFsh\nV0odf2BMTib17gdCengIvXt7npCXS7lVxORkUu/ebejhYaV377Y6CclkUgeDQdsHlVMmsmozfy0G\nahHaYpKPJ3X3f79O0/1Swb3tWTJb8711+m719qKHh1VtlXYgo6ODenjYr4eH0Q884Ncf//hgiZgr\nm9sxmZzUvf5ex472QsUzMbfqIlxGzL8EXJX3ehQ4E7gK+JLTcXnlW4CPuGmHiLnQaOwEVSmlB+0U\nwUlhKzA5mdS7dgULRve7dgX1qlWhmh5UNTZjQWKXchScBb3cvUs+ntShyP/IjtDndOjM39V1b52s\nPkotnpG5JeQUbA89RImgO33lne/hwu0MNVPMvwNcnPd6CFgLfBL4y7zy/wN8sui9p2RH7P1lzh0F\n9gB7+vr6GnvXBEHbm7q9xBqRU7IpZW/eb/aDqpM7B7kRefHm95ce22yrhvPIXOnJyQ66yXWQG5EX\nbw884Hf1GTjdQxmZa08c4OyCOusy5dablOoCtgOf11o7hubRWie01mu11mtPP/30uhsrtA9Oa7lb\nTSQSYXx8nEwmw/j4OBGPPa5mZuy901assD++r4mRe9o16prbiGu5QC5uypu9EiAejxMscs8PBBR/\n9Vf1rXboLOw/IL9/zpWTo909DAaDxBdrQIw8vBDz/cAr816vAg6UKc+RAJ7SWt/uQRuEDsM0TaLR\nKKlUCq01qVSKaDTaNoLeSAIBe3HeuDHU8gdVuyx1y1+utHSp+4hrTsvD7MqbvRIgEomQSCQwDAOl\nFKGQwdKld/Oxj21d0FkK8z/LuTn7D0gpv6slscX30DAMEomE5x3ujsTN8J3yZvZ3UegA91i2/DTg\n11hOb6dm/z8tu++vgW8APjfnz20yZ75wWIzmshxOc+Y5J7hGmvgroZS9mbqZln4787fd5vOVTgdU\nM2feypUAi8Vxsfg6P/7xQf3QQ6Vm9tHRGrwUFwl46M2+HfgNMIs12r4B2AhszO5XwBeBXwF7gbV5\n770ea8na08CGbNmq7IN7H/Cz7Hajm8aKmC8c2t2RpdHzxnbe7F5Qb7vbYambUxsqbTkxdOvN3kpB\nbYf73EhynVJQGgwNyQJBf/DB3Ny5X4S8Ap6JeTttIuYLh3YemXfqqMmLdtvVoVRty7tqxck64Gar\n9uvTKme/drCA5OOlRchuRQgECwS9TfrsHYGIudDWlFtT3WraadRkJzZOAuRVuwcHS8Wm1s5MLWJZ\n68i8k0Sivb5j3v4WnTrq1gi9ddfZqYiYC21Pq+eHnWiXUZPdKLmnR+vubnuh9ardtQpNsXAPDtZm\nKbC7bp/P2nKm8xNPtG9jKFTdtbaKdrL+eG0lc46iqFp6nZ2KiLnQEDp5DbJb2mXUVM0I1TCcjw+F\nqvvMaukUOJnna72Plb5nyaTVsSmuu7u7c76Tjf4tua3fa/+VciPzhfrMaCQi5oLntNNoopG0y3VW\nM3eslH27u7tLRa/StdTSmamm4+GVhSMU8r7T5UUo2FaSH6/f7VRJPSPzasMSC9UjYi54TruMWJtB\nO1ggqh2Z27W7FsFzWt41MOD8nmo6Hvnnruc+ez0dUm0o2HbDzZI+u8+9FvGdnEzqhx8OlSwxa5cl\nlgsJEXPBc9plLnmxUO2cuR21fGZOnYjc6L/a9zi1tV4LiNedy0qhYNuhg1cON50/p8+9GvG14qsr\n27Csw8OLJ858s3Ar5pLPXHBNubzagj31hKy1y+e9bRvcdZf7HN+1fGZOEdC0do4EF49DUfA6gkHY\nuNG5rfVGm3M6Z60B88qFgm1mmNv8iGnVRIZzE7nO6XN3G8I4nTY5cOBOmI/MXYJTuGKhwbhR/HbZ\nZGTeWtplLrldKfXmbv3cYS2fWbkRXiUnuEY72pU7ZyhkbbWOnMuNzBs1xZRMFk6FnHhi9T4OOSqN\nzL34rTolCZKReeNAzOxCI2iUqbFREdGahb03t1Eg5Lmt2YFx8p2icoJV7rMrt8yt2uVp9TjaVVNX\nMlk6/VCtZ3u5OfNGTDE5eeTX2mkot6Kgnt9q/ufw0EPO5vX8OXPBO0TMO4DBwUHt9/s1oP1+v33O\n7EWA3Rxcpz0U7IWpfULW2glHT4/zA76WwDHVWgHKHV9tXU6OftWuO3fyZvdqZJ4vjE6WgHo6DV53\ntos/h+3bnUfmjzwS6qjfbKcgYt7mDA4O2j7oF5ugT04mHZ1p2sFcNzmZ1I88Mu+1+/DD9g8s+5Gb\n0RYjc61rE7tqhaEWwfMqml05IfSCWjoqdpH73CSQ8brTUI+oF38OAwNJvWNHsOi3qiS+egNxK+bK\nOrYzWLt2rd6zZ0+rm+EJXV1dzNl43Pj9fo4dO9aCFrWGkZEwMzMph72K9eszTW1PPum0yRNPXI/W\nR4v2dLNmzV0FOajDYcspqhATpaJoPe/lFQwGW5KyUSnnfV49Anw++7qUgkyVH6NTXeB8Die8uj7T\ntJzzJiYsR7J43N7xMOcsl+/cpxSceCK89FL15w0Gyzs5ujl/tXXksPscBgZMbrzxU/T2TjDnD3Hu\n6i2LKB9781FK/VhrvbbSceLN3iLshLxc+UKlnOdrV9dpDT9/Oc/hsbFYoZA/OABXboe33c++NwRI\nm+nju+w9qyNs3NiY3Mu1ejw3Ei9XOzi9Ryn7aw2F7I93Kq+FSMTKt10p77adl77W1Qu5m9UKbs9f\na056u89haCjCVRt38baH4T0/+D0PPld9vYL3iJi3CL/fX1X5QiUQcH7SN9JoZJqwfDlcfbXzcqOC\njsaDA/C5T0L6DNA+SC9nNDp6XNDtlpElErB1q7slP9W2vdplUs0QOy+XisXj9qNtre1FacsW6Okp\nLOvpscqbjZslYm6o1Gmo9vy1tMvuM6X7ZRj4FADTs9PEhmroJQieI2LeIqLRaFXlC5X+fucn/dzc\noYacMyeGU1Ol+/JHMAUdjS/fCDNLCo7NTGcYi40df+125OaGdNpkZCTMzp0+RkbCpNPzSl3LyGvL\nFujuLizr7vZW7Jw6NLXch0jEuTNnJ0qRiLUGv3hNfr19p1osIF7EXSg3bVCJ0xwMWk7l5l6T8O1h\nfLf6CN8extw7f5H5nylkYNk4XPEReP3248dMHJZ15e2AiHmL2Lp1K4ODg8dH4n6/n8HBQbZu3dri\nljWX3t4Ifr/98LDcqL0e7MQwn5xY9PfHUSo73Htuhe2xMxMzHrfOEvLR0WjWl0AzM5NidDR6XNBr\nGXlFInDjjZAz/Pj91muvp+697NBYAlKKc+CT2s9tmibLl4dRyodSYZYvN9m0qbZAMU5WBbAsIfkd\njoEB++P8/uZMnZh7TaL3RkkdTqHRpA6niN4bLRH08XEwbuuHP31VgZAD9C2TqFHtgIh5C9m6dSvH\njh1Da82xY8cWnZDnOOusLfh8hbY8ny9YdtReD5XMjTmx6O2NcPbZ2+jqCsEK+4nBQF/A49ZZc/WZ\nTGFvI5OZZmwsVtC+YsqNCE0TvvrV+Shnc3PWa7eC0eg5erv6vY7wVny+pUstUVXK5Oqro0xNWZ0n\nSDE1FeWOO8ya5p4jESvyXbGgB4OWJSS/w/Hgg1Y7ijl2rLY5boBDDgYtu/LYUIzp2cKLdDKdxwfi\nBLsLP5Bgd5D4QGN+p0J1iJgL3lHjE7+3N8Lq1QkCAQNQBAIGq1cnXHnI1hIutZzoFYtFb2+Eiy8+\nyJp/eCu+YOHPxRf00R/vr3i+anFyCsyV1yJy9ThFNTqUqVP94J3Zvvh8H/4wvPxyriQGFJtqprPl\npbiZe966Fe6+213b59tRSCpV2z2uprPnZCK3K4+cGyFxRQJjmYFCYSwzSFyRIHKueLK3BW7Wr7XL\ntpDWmS84Ghzr1ctUi07rfUOh8s2dTE7q3cZuPayG9W5jt55MTnpybcU4hczMX3ffzNCpjc6W57b+\nctdcX+Q5++A+VnnjrrvS9UP5wD5OVPNT/NBdIb19B3roIfT2HeiBL6K5BW3cZnhxaYIHIEFjhKbS\ngCd+cnJSG7t3azX8kO4dvkfHhgcKIsStWhWyfQi7CcrSzhmwJieTeteuoKcR8er5eBqdLc9N/dVG\ni+vudo7TXno+w0HMjaqj4LmlOKZ8ubCu1UaxK67f6fs9OZnUQ8M9Bd+zHQ+iL/tSt04+3kY/iEWO\nWzGXoDGCN9QRMcQuGAdvTxMdHWU6770BjvBJPsfbGQLgbW9zOqUiU22UkjYjnTYZG4sxMzNBINBH\nf3+8rsAcToFMNm60TMLlsA+IY5mOx8drblJV9Zc7Buz35ZMfNKW0LhOIUmhqD9LdneDGGyPcd1/l\nQDFuyH3PUynr3ud/d7u7YXbW+b2NeEw7BWw65gvx9ksOen9CoSYkaEwn0Y4RQKqlxogh8/OlJlqH\nSaV8XHNNmD+5/Y4CIQeYYQlf5sbjr1fYO5jTtwBysvb2Rli3bpz16zOsWzded4StSASuvbbQKUtr\nd05wjXREc1t/OQ9+N3PY+f4B8bj1U5snAiQAy2cDDEKhBHfdFWHrVm+88zdtgmuume9EFItzOSFv\nFE6+GV2ZxiwJFRqLiHmraWai5Ebi4omcNtOMhEfY6dvJSHiEtJnOOmblRkaWN7HWKV6+7e8sV98i\nnmNewTduDBEsOmcwGCTulcosMO67r1RE3Hpne+mIlgvYo5SJUmGuvtqHzxcmFDId6y/XV3Tbd8uJ\nfiQCX/uaFV41h88XYXBwHK0zaD3OwYORkutz2+cuPm7TJrjzztpH1/UG9nGKWeC09LNRS0KFBuPG\nFt8u24KcM69xMjOZTGrDMLRSShuG0dQc2WUa5ThRN5mc1LuCu/Qww8e3XcFdeoBJ5znL3l7N8HDB\n1ju8vWAOuS3vQxOpZu6/0XPfbpjP3pbUUOi82NPj7LxY7Zy5145rbp3KyqUhrbSFQvWncS2mnP9F\nI3wzBO9BHOA6hBqesLV6cbeS3cbuAiHPbV/379aO3sRKFQh5YPh+HRse6Mh8542g2rSmjfZKL8Yu\nneh8G+w7cOWcF916s9s5lFWb4awYt/eunGd6uS2/Y+KlY+YDu0JlV0bYrRIR2gsR806hhiesYVT/\nIPSKWkfCw6pUyIcZ1g8xrJWyv57QqlVZb/ZhbezerZOTtS8Fa2fv9VqpNq1pg1cPFjA4WEnAGpvr\n3e3n7faelLuWfNyOwtYD4p0AACAASURBVIt/6o34DJKPJ/XQQ/a5x4eHm2iOEepCxLzNcOwB1/CE\nVaqxD0In6rEIOI3Mdxu79eBgUivVOEtDM0WsmbgVmHya1anJjcidt+Z1SMtds9u+tNP1+P2FxznV\nVyzySlkdnkZi3Gbo7TvsxTw/ZoHQ3oiYtxEV56aqfMKWG5k3cg65HouA05x5LvBKY9vt7oHdCuoR\n11rEvFmUa1stc+a1Uqkj53aWy+29djpfboqhmZYhdYvSA1+01o4XryUXc3rnIGLeRriJ6FUNTiPk\nwcHBhs6l12sRaFYEtWLawfHLjnotBtWa2b1qsxtRKjeSTSZzbU9mR+hKh0KNcV6s1JHL3z/ApN7O\nbj2U9eXI/35W0yH0yvpRbz3GbYbmFiuqW36Utw/d1cAviOA5IuZtxPCw8nzeanBwUPv9fg1ov9+v\nBwcHGz6X3sq5+npo1cjczgHMTbtCIffzvV57P5cjmdS6q6vwfF1d9udzmjNvtGm5mHJz2Pme8ANM\n6h2Usxw1dqqmWLgHB+s/X/LxpA7Gg5pbOL4F40GJ7tZhiJi3Ec0amdt6hFcxcq71vO3sRa91a+bM\n3YiZW2epej2xvWLpUvv2LV1qf3ylzowTXl5TOe/yfA/yr/udfToa0a58qlnOVm0HNPl4Uhu3GVrd\norRxmyFC3oGImLcRlebMq10e4jRCzo3Ui7dVq0JV1V+uPZ26rrvZ3uxuHKaqWcbktRWhlvtRrn1e\ntsvLjlelNei5++q02mJYDXt1aQVtyr/3TtMldlurp4aE5uNWzCUCXBMol+IznTYZHY1mYyRrZmZS\njI5Gj0dpsmPCIX7l3NxcSUS0E07o4brrXnRdfzptsm/fhoLj9+3bcPz4SCTC+Pg4mUyG8fFxIvXm\no2wSkUj9YTmrSbeayxtertwuaJ4TTiFLa4kE3M5BB+tJ1WpHLnqdE7n76pSX3ut89Xb3fmrK/fsX\nQKRioUGImDcJp1jbY2MxMpnCp1cmM83YmPPTyyn2uGEYJBIJDMNAKYVhGPz5n5/EwEBh4Ody9T/5\n5GagOFD0bLZ88WKaJtFolFQqhdaaVCpFNBp1FHS/376e/HK7MKlOoTvtPnInUX7726Gry6qzq8sK\nJ5pPrYKZH9fdTXktOCVMqZRIpRyRyHxClmJy97U/3u9JvvpKnSu7e+9E8X31Mh6+sPAQMW8xTskO\nnMoB4vG4Y0zy4pHz+vX2SROc6p+bsx8mOJUvFmKxGNNFT+Hp6WliDgoYjdrXU1xebDHYssV9UhMn\nUR4amrcAzM3BHXcUCnq5pCXl2LixuvJacNMJqoVKqQN6I72sTqwmYARAQcAIsDqxmt5Ir2198/Hl\nrW35cuseV7J4uEkKk2vbxo3excMXFgFubPHtsnXqnHk5nJzjHnmk/Dz38blr0Ibfr5O5CcCiycVq\nne/sve6tbTFTy7K8RjuAVRNtzM1cvZt5+VqvyS2NnJf3cslY8SoCt/4O9a5gaNQ1Ce0Lks+8M8jN\nmeeb2pXqwfpc5s3dPl/w+Dz7ceySVOcnbnao37auLD/4wXKOHSsdhXd1hbj44sWb4zgcDpOysfUa\nhsG4F0m9a8Apx7cTuZ+6i69Ny2h07nQvqPa+K2VZXsDdvTf3msSGYkwcnqBvWR/xgTiRc0s/mHb+\nHAXvkHzmHYKdc5zPdxLF89a289xOdtarrz4+YVfO+c6O17xmC0r1FJQp1cNrXrOlruusmxbnfC83\ntdEI3FxuNQ50lebqqxWARn0cjc6dXg5zr0n49jC+W32Ebw9j7rW/KLem8hz5/g6V7r251yR6b5TU\n4RQaTepwiui9Udu2eO0sKHQ4bobv7bItRDO7Ha6DzFSys5ZZ0zM5mdQPPxwqMOvXulSu4VRar9Qk\nW2OzluVVszyr+NIHBuy/Cl6aw5sdQKUhSUiKg7T87SOuA6yUW1JY/JOs9r7korYVb8ZtRsmx7RrZ\nUPAWxMzePqTTJmNjMWZmJggE+ujvjzuOjAFGRsLZpWGFBAIG69aNzxe4sffZ2Cdzy8+KR/9K9XD2\n2dvKtq0llLO9xuMLztZYr6l50ybr8ufmrBF5NApbt7ZP+5qBaVoj1IkJOO00q+zQIWuUfPnl8NWv\nFn5lVPc0+oob4fXbC+oxlhmM3zReUveGDTBbtOijpwduuAHuu886b1+f9fWs5mvou9WHpvSZrFBk\nbs4UlHXC5yDUj5jZ24Ry68jTaZORkTA7d/oYGQkfX8vd3x/H5yu0Nfp8Qfr7i2yNbuysNjZBy1xf\nvPwMtD5adklcyyjnfl3J1thi83yOappRq7d5jq1b4dgxa5x27Ji3Ql6uHdWanxtF8ZK9qSlr01kP\n8zvvLP3K6NkgDP1NSV0Th0svKhKBu+4qXEYYCsG2bda9rieeQd8y+2WnduWtnJIQ2g8R8wbjtI78\nySc3O4q863nu/Ak4J2wWKJdb9lZuX8twipTR11deWdokOkq1zSh3uY1oW7V9nWa2rxYqreV2NEYe\nLr0AJ3GNRODgwXnj9sGDtQYiKrz/lx9JEuwu8s3oDhIfKFVoL3wfhIWDiHkDSKdNfvCD5ezcqWzN\n5WCt2y4XLMYpyEwJuYXKyaTrbnog4PzULbevZZQbgpRTlhZ4CNlFiau2Gc0acdXa12n1iLBSB8S9\nhcAEwliPwTCc8IWCvU4iWgt2bba7/1/9zMVc6/sexjIDhcJYZpC4ImHrzQ7eRDYUFghuJtaBbcBz\nwC8c9ivg88DTwOPA+Xn7rgWeym7X5pVfAOzNvufzYM3fl9s8d4BrgKfN5GRS79zZU3a9dvmtDu8V\nl9czOZnUw8PdJefeubOn9c5uOezSSNldWzlvrCo8hLxwbnNOgJOs2lGpYU5geRUb/mfKrosu14ZW\nOKnlZzkr52TmLuZ9aT71ru6ADkVCniclGfzGpFb37NYMDWu279YMTGrQ2uezb1ubJyEUmgxeJloB\nLgHOLyPmlwM7sqJ+EfAf2fLTgLHs31Oz/5+a3fcYsC77nh3AZZXa4amYJ5Na9/To5Llo4ya0uhlt\n/Ck6ubU+t1+nIC3FSVYeeSRku6/WTGrVUs6bveVU6y7tpCwuo6N4lQ3OOQGOUfNDu9pORnLwEW34\nn9GKOW34n9HJwUfyKyu4r4o5R7FrRaa5gutwOL9TUpL8e1kpuUowqPXSpfafVa3pfJ2ykyUnJ7W6\nf5dmeHh+27HruKCLN7pQCU/F3KqPcBkx/xJwVd7rUeBM4CrgS8XHZfc9kVdecJzT5qmYh0I6eS46\n+KnCJSDBGHX1yJ2Xlc2L9eRk0sqkNlQ4gt811EYj41biVQJyl4rkVZ52pyhxoGoSxmo7GcnBR3QP\nRwrO08OReUEvuq8Gv3YUEyfRDIWquiU1U01GOTsBzO/fhULWlt/XcxvRz01n6huPDep7dig99BB6\n+w70wBfnl7UZu3cXCnlu277b8VpyX7uBe0c1Dw5rHhrWPDisB+4dbcStFtoct2Lu1Zz5K4Bn8l7v\nz5aVK99vU948pqaIDcB0YXwUprshNlT7nGr5+Wjj+Px374Ow+nOawCSQgcCk9br3QdrGA7tleOUu\n7dJDyCkLnVO5E84JcPpqclRyEw8+/6vy4TvWcZTCLF9HCbD5S6/NXVDBvjifQlG43AksSXHK5DU1\n1ZyvYz2BWaBwLvngQWvLn1d2+qzyy90k10mnTU783Z30LtH4FJyxBD65GtadMm1FcZuZsW/wCvvy\nnO/B27/zJEMnHgA/lu3SD0MnHuDt33nS/U0RFhVeibld3iRdQ3lpxUpFlVJ7lFJ7nn/++TqaWMrE\nModym+UodtgtLevvj5dEULPoLlxaFovRu2OWdVfB+gFYdxX07piFzZu988Du1E6Bl+7SLjyE3DzY\n3VA+AU71jkqVOhnFDlQZ7LORTGVOtf4pup4I2+1/dBVoRoQxp1sfCnnjfOcmop+bztTYWIyAr/Au\nLvHDjf3Wc6Qv4JBC9bn5cr+/tJM3dMKB0qekypYLgg1eifl+4JV5r1cBByqUr7IpL0FrndBar9Va\nrz399NM9ai4QCtF32H6XT/kqhnR0Wj8OcPbZ2+jqml+E6veHWLPmrkKPdKehx9SUNx7YbbIsqyaa\n7C7tVajWSCRSkoI2kUjUnPPdqTNx4olWeTXpNAHb+2qo/baHOqViheasJ3f6CmzZUvtyrPy+bSwW\n4dprSz8riBw/JpWqbLFxWsq5ImAta4v39xP0FT1mj/jgy/3Hr+mrX7Xp5Dk9mWX9keCEG1u8rjxn\n/i4KHeAey5afBvway/nt1Oz/p2X3/Sh7bM4B7vJKbfDaAS75Bl/JnHnx5hTSsdpsZCVUOymYm0xz\n64Hk1bxzq2hyOqhmhWqttk1dXYVz5pYHdlIPDrrPmhY68ff5lRbc1+TgI47z+W6czRp7/d59Bdy4\nT5QeU9mXwuk5cM8OVeAEZ+zerdXwsA49sFuHPjRZ8Zp40GaefdiaOxcWF3jszb4d+A1W2LD9wA3A\nRmBjdr8Cvgj8Cmu52dq8916PtfzsaWBDXvla4BfZ93yBFi1NS64PHfdm99/icx0X2XX89DLnrspd\n16XnVEFqVLBSo5bzFKoFybtYE7XcNqWSWVFR2b/WMje/321/MFMxLrtTu1rt0V6Jau6nm75t6TGl\ny9eKHRAnJ5N6165gwTPg/iGlv/FYfatiBu4dtRzf8oX8IXGCW4x4KubtsjU60Yq6RdmKubqlVADr\nHplr7X4hrcuhka33c7Gg1zusqvUJ36AOQDuOqO3IroTMEwlLoEOh8m0u9zWw+yi6uqru/1Vsdzv2\n26r9GroJOWB/jPVZlft+NSoxkXizC1qLmNdENRmL7Hrku3YF9eRksmTN6eB3Bm3XoDqSe4I6PcUd\nRteOS6y8HFbVYr73eoiXvT9J0MGiJUa1rA+v9fzVKNy8waXyaC8fv9/+dvv99k2p1TRebFYPhdpH\nuO2o9mvodPyHPjQvxF//uqEHBkoD/HTKzJSwMBExr4Hk40nXaRC1tu+R29Xhdh6+hCqfWI5rZ3Pv\ncfl0LjvarSXvopfz93kdA8PmWovnND2nxo7J/OEOHS6HNg8O2t86J9N5LR9PodVgfuvubl9Br/Y6\n7T62yy5L6qGhwg75jh3BAkFXyrrX1fbfnILICEK1iJjXiOOP0OWv2Wl072a0X9qY6oTDi+AnFQOV\n1CLMlZ68g4PzQ1C/v3zy7bzzKwcxLw78URWVPucaOybzh7oLVpJPjbfHdfPKGYHadVRaq4Eo/6N9\n4AH7qbLt242COnt6rI6N2/5btYMCQSiHiLmXVCGqTvPububhHc/tckjgRVjSih2CWkam5Z68dQw9\nDQcxz7W1nIXBdp+ba6tl6KvzTdj1d7jKYXcJSmn92tc6dwjKecW3a2hRL2ZunJxYh4aU4/2w6ziU\nTHVE/kftHXhBKELE3EuqGAZ4OjKvgXodwpxDkjJfl0MHo6xVw+nJW2lSuJi8zyKJ5eBn13lJJpP/\nf3vvHubIWd/5fn7SdGuQLwPWeNrh0iU6hFmctQ8bJuTxc06Obdoh2ByH6xoTGYyNo3XPyXpM1jlL\nopNgONHB7PoJ9nIYE2WwcVAxDixnwT6PTbDbMxN2p1kwC/YAToPTaTUOGdluhwl2Z3pmWu/5o0o9\nulSVSlLp1v37PI8etd4qVb31qrq+7+V3MWNjY3XbxsbG1reNj483Xd+M11xz4+/c4ci8WKwm1mhe\nMx8fj3adP6zbWlXQ+zkyj9rdrJtjPXzIOz9C48g8qKPjaa869oLhne/tvAOvKDWomEdJO9m3olwz\nHwB+I/NWo/yWU4t+T96gp6XnieqfnkUwlojjilfTeUmlUp7XkEqlfLeBhytf4+/cxZCw3gjOMlV3\ns1Qq2nshbAiDWiO6fqyZe51nfHww6/LFJ4rm8j8bMw89Ui/kD33tJZ5GcH4dHd+23vZ3OjJXIkHF\nPEr8/mNTqVAj1Lat2QeI11R9mOngdjwB6mh3ZO5UsuWQzK/+rV4pv6d24/kbzcZDmH93OEPfNmGD\nydT2l/phzd7L5C3tjtKr9+v0p53kKNUkKW/7ozc39dWC1sz923rNs2M7Kq6UyvCgYh4lXqOxdq1i\nRohisegrdn6GWmF89D39cdtdMw9Jp2JOYz28ftMOh7JRGvUH0e7IvF+E6VR0QieTJU336zvf64ym\nWfPMsubXWfDt5//Cz5s68FGl2lU2FyrmURPWoTceH74IGx3QrmV8q5F5kF9+W+baIel0mp2qugb9\nhh0uMkftbh90Hr/l/wj7S23TKzHvyIK/9n5953udde4Ofpewv2mxWDTxeLyt/ylFMUbFvPeEmcvs\nw0g9smm7hs5KcWamvVzaLdbMI4mY19blNBu5jY+PrxvA+U6zn3lm64N3Yf4dpQFYEEFRgSPqL0VW\npzDT7EHt1pFvfe39uu3v2u4MhK2bsz146aorV0plw6Ni3mvCzmWGCb3V6unuZz0e1bSdz/CiODPT\nVkchKFBG17HsOyCoozMzM9P0UB3fsiVc242AY3a/1ufboVhsXpkKY2TXavTb6fJF9X6FtUjayn8q\n3vIVch2ZK61QMe81nj4pbT4Rwpj3BjzJ2p0K96UPi7n9Hpl74jH70NGsRr/Mv7sg6Cft1+yA13k6\nOXer27Pb5Ysobv+gOgS5e+qaudIKFfN+UH0yBYl50Bp6mHnHgCeNb/jWhg5Ey6n4PgzjAtfM+0HU\nC9Zdmn/3KjlHbfW8LndmxiOozNiLhnf+dqSeFlE2d5jbs5sOShR1DeoQ+HW64/G4CrnSEhXzsEQx\nTAk75d74hAjat0rAkyzMyDzUVHyfzKx7LWCB9MuUPASeHZuviTn65WgXsr1u7VZ+0WFjIMzMzKwb\ndMXjcTPTsAgf5M3ZLu38dJ3aUnb7GAjqcKgVu9INKuZh6KRL7jUi8xru+L1qn0BhxDzgSRbmIRFq\nKr5fZtZebdmP+V5jupt9iLievksO90nP2zyMX3Sr2ABe9gZAnaAH2QhGZNLRdJweeTmGovVSgPqX\nK52hYh6GdkdrXhY8tUOA2ge+35OsVjzCTLN7PcmqT8oQ676BmdRq695PYfW7rl52IDq2koq+nr7G\ngLMh6tMlrUbmYcKO+rlYxWsc16O2EQxze3YSf6gbikePGuvwYSMHDpjUw4fN2OVH+94fVjY+KuZh\naHe0FvSEapw/DCMeYc17a+dHG+vc4onhOzIfxPClvmLRPenD0Kko96CeviPz/f73XlQju2LRmPGt\nJ+svpyGWeKuReaCPfs15wvRno8TvfDXViozi0aMmeeiQ4cCB9df47CGTes/RjRBmQhkiwop5jM3M\n5GR75UtL/sdaXq7/nM9DMllflkw65VUyGbjnHrAsEHHe77nHKa8lk4HFRWe7MfXbVlYgl/OtVj6f\nJ9lQjySQb9zxz/7M9xg9wa8tg9q4GzIZKBTq27pQaG7rsPXpop5TU3liq1JXFjsOU/vwvPds2yab\nzVIqlTDGUCqVyGaz2Lbd/skvtDFX/g5sWwQqzvuVvwMX7gcgOZYkP910d9QRj8dblmcykEp5f9/v\n36tbfKrlWx4W24btV5eR++aQAwfZ/sgce378Y1Yqlbr9TsQqnLlngUrF+XdtdWspSpRsbjEPI7i1\ntPMUCiserlDbn6+QZpHY+zKk084DpIkOhCWTyVAoFLAsCxHBAgpA03OmUvE5aY9otyPVLbYNe/ZA\nqeR0iF54Idz3elDPiYkMO1+8kURZoAKJo7DzdpiY8773crkcKysrdWUrKyvkAjpxfuRmc5z85c/B\nh14Nt8add1fIrW0WhSsLZC4IVqFsNhuq/M472/v3CqJctpmbS3PwYIy5uTTlcvO96lMt3/Iw2DZc\nZ5dZ/sA8nLcKAstbVlk+ecpz/6XV1c5PpijdEGb4PiyvgVuzB80ddpEtIvQMcBRTvn4Li72c4vai\nn2vm3fiF97KeIe+9sC6IYWgVQ7828M977km5aUKbvQ9aWbO3eYmBeFn/f21WzJe/1XzOqCMDW5Yx\n7D9cN50e9LIOH+7uhIrSACGn2cU0TtsOMbt27TKPPfZYf09q28409tKSMxp7zWtgdrZ+n/FxuPvu\njufV0mlnwNiIZTnTdXV1yWadqfUqyWS46eIqu3fDXXd5bxNxRuj9orFt8/ng69i927nWtTVn7jSb\nhb17W5/Hr4HBo5EjqGeE2Edsrr34Wtb+ca1pm2VZLLaqewPpO9KUjjW3hbXNIj+dJ/tAlpWTK0yf\nC7fshK01U9SxWJKdOwtMTPTp2t1mv+22NOed11zn8nHh1Cs+z2U7YGEhx+rqEonEJFNT+ZZ1tMtl\ncgsLLK2uMplIkJ+aIjMx0bRfLAZm9iBI8zEaScZiFHbu9DyOonSKiHzHGLOr5X4q5gHYNlx/PZw4\ncbpsfBw++EF48MHIHu6xWPNSOPhoaxTCcuaZ8OKLzeVhhG1Q+HVCZmZaC7pfA0P/OzBtYB+xHXH9\nzgo8AJw8vS2ZTFIoFMi0+duvH/Pk6Q5hcixJ4coCudncutDv/zU4b2vz9xMJi4suWuzgatqjtt86\nOxsjFmv+/SoGCqUUu6f+mUrl9PW06nTY5TLZ+fm6NW8/Id5+dZnlf/Okp5gLMJlItOwQKEo3qJh3\nSblss3Do/axur5B4xjFOmqgOyFMpeO65yM4VemQeFVGM8PvNli3OiLyReBxOea9frtPtyHxA1I2i\nnwBmgWMQf1mcez91b9tCXsU+YpObzbF0bInJbZPkp/NkLsgQ+2gMg/M8mP1fIeY5GhUuuaT3nZ/a\nn2z/fmdk/sgjsG8fPPMM7NgB734//C9var/TkZ6bo+Sxtm0lEixedFFd2fZH5lje4r8Obi65JOQV\nKUpnhBXzzW0A50O5bDM/n2V1RwVisHoezN8C5Wl3h+XlSI3F2rXD65pOLbsHiZeQB5XXks87MyqN\njI31sJG7Z+lYjWHjhcCHgFuhctMamWuuARHsM88kvX07sViMdDodysI9c0GGxZsXqXykwuLNi+sG\nb5PbThv1PeOjX4lEjwwUG6i16dy3L8+DD45x++1QLjuTLOUy3HUHHPmG9/ePH1/y/Rf1M1LzKn8+\nQMjjOKN8RRkGVMw9WFjI1U3bAVS2wsINNQXZbGSC3rG22rYzhInF8DeBDzjp4iJt+9F0c85u6Mbv\nKJNxbBpqfaVSKW83wCGiVlzryo857zaQffFFSsvL3busAfnpPMkxp1e5bwGON/STYrEkU1PhOj/d\n3ia1zgKzsxk+9amzadTayknY91nvR1i5POn7LzqZSHif06Pcb1+ANSA7P6+CrgwFKuYerK56u3qt\n7qj50MK/u10atRVs0um0/4irOlVedbUqlSLtYHgyiHNW6dbvKJNxlkaqNunPPTfUQg714loleQLy\n7nJPDlhp+E6nLmvgjNgLVxawtlk8+qzwuZ+kOBVLAUIiYYU2foviNmmcrTp+/HnP/Z59pkIsVt9G\nx48n2bcv7/svmp+aIhmrf/QlYzHyU1Oh9q1lpVIht7AQcCWK0ifCmLwPy6tfWdMCI3R1EcoqbBSv\nYUqOUsegk5VE7Xc0AtTliL8ZU7zgdLuLTxS2TlzWOqqbb/7uaG6T2uPH45bntVqW4zK3f79lZmfF\n7N9vmenpYst/0dpQrNbhw6Z49Kh/Pdx9/dzR5MCB9i5MUdoADefaOZ5ZrR7CHJ3u/OnUTuakUMlR\ngjJZ9CKW5CDic/abfsenb5cGlbR8xNzqQ+cqOH939LdJq/+ffvQz/QRdfcuVXqJi3iV16TofTpmj\nlzfEUG8zcEgogXYJFSTEfXoV3Ye6uO/FDusXiNeTexAj814yqMxx7dAQ/KYIJtlwj/QqtWZjP8cv\nR5CTv7v6uWjAMiAGLJNKdVevoJmtfvx8XvHYk4cOBY7qFaVbVMyjpstRWztRvMKmLS2OjTU/zGsF\nPSqRDUowM2yC1ymDXkIIS0MK3uIZZxgrlYoktaafWLbqyzWOvp38QUUD9SPp8fHe5vDux8RKO9Pz\nihIFKuZDRjsj87BT8lYq5X3MKOY1a4kyOfWw0ou54X4QkYIF3XNBfTm/vk8qFf5+VxTFHxXzIaOd\nNfPq/nWjpMZ86cVi61zlvR6Zb6QH8xBdY+jRX5tzy3VLRw2x1oM6m0F9Ob9TRxlPvp3rUJSNhop5\nEAMydOo4J7XPQztwZN7rNfONMr1eZUiu0WtdVma/Zma+9eXmnf06IB5JfzyNOg8l14UwSHyDTuP3\nb+TXOUilrLrvePRRA2l1HYqy0VAx92NIHtpt4fM0LYJJNjyEk2CKqVRvrNmDnrrDbgkehiG4Bj+L\naXnoPlN8oqE+bSx/+LpbHrac8waMzDv5l/GaiRofT7pr6f7VbnXcVtehKBuNsGK++YLG5HL1Mckh\n8gAwkeORr9zGDRpiDNUYaJZlUSgWyfQiIEpQxDjbxrZt0rfdRuyRR0jfdpsT5Kaf+dGjoNOoeBHi\nF2rUJM4ld/+e+sKgnOoN97NvICS3PJ/Pk2yIKZxMJsnn8x1FKMxkMhQKBSzLQkSwLIuzzipw8mRw\nm7b6V2x1HYqyWdl8Yu4hjIHlw0DDQ9sGskA1dcgakATyV1zRcfKNbrAfeIDsTTdROu88TCxG6bzz\nyN50E/YDD/S9LqOOb/jQ1WdYOrlcXxYUV77hfvaLqV4t9xLf2qxsnfRzMpkMi4uLVCoVFhcXef75\ncPdm0L9iq+tQlM3K5hNzv9FM0Chn0DTEtvQM4wnkPvOZgYyGc29/Oytb61NXrWzdSu7tb+97XYaS\nNgKV56emkFPH6wvXjsPCvvWY7OtkMvXx5mtpuJ+npvJNYU8bY603im/UHcOw/2JB+1Wv45FH4Oqr\n4U1vgquvFr773SsGlzdAUYaAzSfmfU9RFgG185yA38BlyZiBLBcs+eRw9ivfVLQZqDwzMcGNX/0M\n8s9HwVTg+FGYv53k38+S/56HcN95Z6j7eWIiw86dBRIJi3ZjrUeF179eI63+FScmMnz/+9dy++1S\nk0HN8Psf+iz2IWtM6QAAIABJREFUddcNJm+AogwDYRbWh+U16tbskVAsGsvH8tiC9vyiI2oH6+GH\nPY22Jh79S3Uf6sTlrVg0xTeMGetmjHwEJyb7G8b8f58Rup8bw+tPT7dfdV9jvSFxLVSUKEGt2Tcu\nxZkZ/8hvYR9eDaFBDTifOxCC4tGjJjk7WyfkiUcfMrkD0+vWxl+bFfPlb238xChNdBqMplOBHmJh\nb2UVH7bqfm50IKbIe9trZ0UZclTMNzjFmRljidTHZG/Hxc4vuLaHj3Ko+tQEOpk48KU6Ic8dmDYT\nB/YbHp01sa/9pbfPdNsnHF7RqqOfwWiG3O0yKJJcKtXct/Srut/IHCyT5IV6QdeRuTLiqJhvBroR\ntCBn3y45cEDqhDxx4KH6KfhHHupO0IdVtLx+j37WNaDjMAxR08JGkmvV5/HyYXfiwDs+7BZ/Nzz3\nhKJ0SVgx33wGcBuJMP5CA7DwrXUT2scNrFJv6U58K4V/rHR+gmGMFeBn6AbtO2l3io9PV/k1Jebn\ns6yulgDD6mqJJ598HwcPCnNzacrl/hiJdeIw4nVJVTc6cIz5nPcC4LTpEpO9bWdFGUJUzDcyQZbU\nfi5NfuUhKZdtTp16AWOcz8+ww3O/tbFzOj/JMMYKCOpgdBqMpt2OmI9aLtwYp1JpdGZ0fqDVVUfo\n+yHoYazZG/HrAGQyGSxrEagAi1SFHGDSig0s6I+iDIpQYi4ibxGReRF5SkQ+7LHdEpFZEXlCRA6K\nyCtrtn1CRL7vvt5TUz4tIv9DRL4nIv9VRF4TzSWNKL0YQfsJzLXXwlVXwdhY/baxMbjzTuxymfTc\nHLGDB0nPzWGXy6FOVy7bzM9nWVtbRsQpO9c847lv/OTz7V7NaYYxVkDUHYw2XdoAX7fL1dRa4Kkq\nlRUWFno/q9HgYdmSVm5qo+hlqii9oqWYi0gc+DRwOXA+8F4ROb9ht9uBvzDGXAh8DPi4+923Ar8C\nvB74NeD3ReRs9zt3ARljzOuBLwD/Z/eXM6J08uAOg5+QrK3Bvn1www3107/33IN92WVk5+cpra5i\ngNLqKtn5+VCCvrCQaxoB/o7sY9w0B0HJvqyLSaFhfIpH3cHoZCnBJ+5qYmtr9QwbDtW2bdLpNLFY\njHQ67YTtrd8hsFNanaQoFluP0luHjO3fCoaiDDthnqhvBJ4yxiwYY04A9wFva9jnfGDW/ftAzfbz\ngUPGmFPGmBeBx4G3uNsMUBX2bcBPO7uEDUCv1oCDhOTkSfjiF5umf3MLC6xU6tezVyoVcgsLLU/n\nJQiXMcstcjvxE8+BqRA/8RwzZ6+w91ff2ebF1NDJU7zXtgNRdzA6Hel7TOl7RX9rJEw4VNu2yWaz\nlEoljDGUSiWy2expQbdtuP76+k7p9dd7tnX1J4zHmzYBzk8aLmTswMPpK8pQEEbMXwH8pObz025Z\nLY8D73L/fgdwloik3PLLRSQpItuBS4FXufvdADwoIk8D7wNu6+wSNgC9WgNutUi5vNxU5Jfow6+8\nFj9B+N8ST3Hqze/GXPomTr353d0JeZV2nuK9mvlorE+Uw8QIR/r10d/AMRo7TWNYVz9ye/aw0tDp\nXFlZIVftdO7ZAydO1H/pxAmn3INMBu69d/gmWRRlFAkj5uJRZho+3wJcLCLfBS4G/h44ZYz5OvAg\ncBjYD8wBp9zvfAi4whjzSuAe4E89Ty6SFZHHROSxZ599NkR1R5BerQFXBaadqvgk+vBNAFJDmPjf\nfnS6Th+Kflm/RzlMjHikPzGR4aKLFrnkEsPrXvf59sO62jZLHp0/gKVqp9Nnu285OlWuKFEhxjTq\ncsMOIhcBtxpjftP9/AcAxpiP++x/JvA3rkg3bvsCUAS+DXzTGPOLbvkk8DVjTONafB27du0yjz32\nWMuLGjmqI8dawUkmo3uqbd/u/UBNpeC55+qrUi6TnZ+vm2pPxmIUdu4kEyLWerlss7CQY3V1iURi\nkqmpfEuh6Pac3ge1HbFeWgK/e1zEEd5hpfYaJicdIR+UyqXTpEul9Ux9tViWxeLiIutWj160eM4o\niuKNiHzHGLOr1X5hRubfBn5JRF4tIuPA1cD9DSfbLiLVY/0BcLdbHnen2xGRC4ELga8D/whsE5HX\nut/5DeDJEHXZmEQ4PLGP2KTz24ndKqQ/JNiXbncs18fH63ccH3eSdDRWZWKCws6dWImE48GbSLQl\nqqdHgBUuumgxVCIPv3X6a598srMReuO0uh+dzHz002/fa6RfPb8IbNnivPcjfsDSEnmcVLu1VFPv\nkk77f7dLd0dFUUIQJrIMcAXwI+BvgZxb9jHgt9y/3w382N1nH5Bwy7cCP3Rf3wReX3PMdwBHcNbV\nDwJTreqhEeCCKT5RNMmPjhtuZf2V/EM3ScfMzNCGPxWPJC3VV3J21hSPHm3vgEFxQ7uJwuYRza04\nNmasVMqIiLEsyxR72a5e0eT6FQHPbdOiGz54PYzwGWf41wmMGQtIEKMoSksIGQGu5TT7MLFhp9kj\nIn1HmtKx5olQ62ew+F8sZ3Q3hKTn5igFGNhZp06xeNll4Q8YiwVPrXc6ZZ1OO6N9FxvIUp9bPplM\nUigUIs8F7nX+Jqwe/sZ+S0EveYn/mrhlDXZpQFE2AFFOsytR0eMp2qVj3tbvS9sYbHS0FuSnpkjG\n/G/FJT//JT/8ps8tqzPjtOrv1iCkOeqFHBqsu6Om1W/Yy9/YbynoeZ/gPyLqK6YofUTFvF/0wT1q\ncpu3iE0eY7DR0VpQXaePnzrluX2y3XXzqCzBbdsxHrzmGs8RsZ90LkUhql4dv1a/Ya9/Y681/GGM\nxqcomxAV815TfShfc03P3aPy03mSUm/oljwB+W+MDb3jbmZignv37SN5vD5aXPL4cfJf+UqbB4vA\noLDa+Qpwq/KTq3POceLO20ds0nekiX00RvqONPaRkB03v47fFVf4xw0YlHP2MEbjU5TNSJiF9WF5\njZwBXJDBUvUlEu0pnyga609SRj6CsW7GFC9JdWyAVCwWjWVZ4Y272k3J2rj/zIwpXn65sfbvNzI7\na6z9+03x8ssHY0AVwoiumEqZ8S1bGlJxYuJxzJ6PTZtkPllvjJhPmuITDdfi1WZBOdBrt8fj9eWD\nYlRyyyvKCILmMx8CwlhVeyVsjoouHrLFYtEkx8frRCo5Pu4v6O3m7fbbv49W98Unisb6pGXkVjHW\nJ616oW2VfNv93VKpVJOYA2bHDsz0p6kTc27FWJ+0WrdBnzp+iqIMPyrmw0ArQeilO1G74tqA5SNS\nVirl8wWrvc5Ku/s3XpuP4B89WjSHD1vmwAExhw9b5uhR7+stPlEMHjkHdcRq2lFEPNtJBLP/IUfQ\n9z+EmX3Ueb/p2mnzpZd9yTwqj5rD8S+Zo0w3H7864g4anSuKsikIK+bqmtZLWrkSpVJO4JZ+ujGF\ndF+KiTTF7AUntm/F657xcwfzi7LW7v5Vqsk8amOAj4/D3XdTvgzm57N1mdtisaRnuFJfN75tFos3\nL3q7YkHTb5ZOpyl5tPPEBOzfD6sV2Fo1xn9kGnP7Lcjq1tP14zg7uZ2J9TxFLslk87lrt2nMU0XZ\nFKhr2jAQJtFJ1Ak/qnSZvMXPuGsSvN3r2rVq7tQKOiCZh1cKVr9c3b5ufNVyLyO6YtEJf1sjovl8\nnkSiPoxpIuFkl62YGiEH2HdDnZADVNjKAjfUV6JqsOeX+LsXceUVRRlpVMx7SaMgePlL9+rB3KXL\nUD6V8g7dCd7ude1aNXdqBR2QzMMvJ7dXua8bX215iMQpmUyGP/3TG5mYEEScEfktt8Cb35ykyXX+\nmR3e9aOmvNoG1XP7xTsf4rgBiqL0HxXzXlMrCH7Tx714MPuIpX3FFaTTaWKxGOl0ej0XtW3bdeVc\ndRWFsTEsnKl1CygATXJW7YwEuYNV/bVFnNf27c53I4pHv/ty2PJHcPS495KRV2rW/HSe5Fh9+yTH\nkuSnPToTtk356u3M3SccPCDMPbKdcvn0bMru3Xt5/PHP89/+m8V99wlvfauTiWxromFkveMZ7/rF\nn/dvA/XjVhQlBLpm3k+6XMdum4asW/YVV5C99966nNTJZJJrr72Wez3KC9deS+bBB09n7fJb/w9a\n57ZtuO46OHmyvtxd525bvBsywO2+HO56IyAwfS7csrN+attvzRwcP/DcbI6lY0tMbpskP50nc0HD\nfrZN2b6O+ZtOUqmZIY9Vxtn5y3c3H7emzctXncN89udUYu6ygNeaeTLGzsJOJjI+iWx6nVFPUZSh\nJuyauYp5Pxnwg9nPWCsej7O2ttZUvp7a8vQB/AXdLw53q++024lp6Bxs+SNYqxHv6XPhhinYkYCX\nbLVCpWANJJ1m7rYSq+c1b0okLC66qKb+Hr9v+fIxFn7vbI5vWWb5RJxvfOESfv3+G0n9LMXWya1M\n5af8hbz2uMOSClVRlL6iYj6sDPDBHIvFaOf3FhEqtSNuPwvvKl4dk1ZJTzrJJ17ThvLHxlkH8MB8\nJIJ7Oxbj4CMGYvDII7BvHzzzDOzY4Ri5/cmf1Jyj3zMviqJseNSafVgJYVTVKyZ91ln90pg07V+7\nLu6FlzFf0Nqu17YwyWhq2jAe8659XOrLy3aZufQcB2MHmUvPUbZDxHu3bYjFSDzjCPntt0O57PRN\nymW4/XZZtzkAuvYgUBRF6RQV801EPp8n2WAUl8RJ49lkuZ5MkveyLG/Xyjqfh7Gx5v3Gx5st1ztI\nRpN9Q9a/3O0YlOUy5t/3PVZLq2BgtbTKfHY+WNCrdVlbY2of7PtzaMzSurpq6jOkqbGaoigDQsV8\nE5HJZCgUCliWVWehvtd9X7dct6zWObn9BOqcc+pH1gD33OMEW6mSSnkbv+VybSej2fvWvczsmlkf\nicclzsyuGf7nJUh/533EPlDiG2fdQMUk6r5XWamwkFvwv76aukzMOlPrXtRlSPPyIBBxEqQoiqL0\nEF0z36x0u77rtX4+NuaIV21Ql3YM/DqNCldbrSM2ex7aw/LK8vpa+uyts8S8+q0Cl1QuCVWXNOBl\nxlc1ErRtm1wux1KpxCSOP/76Fav1uaIoHaJr5kow3aau9PIrP/vs5uhs7QTF6XKa2j5ik30gy/I/\nL9cZxT2zzce/ezLhWe51zjz+SxG2bZPNZimVShgc0c8C64sDGrFNUZQeo2K+WYki53ejMd/zz3vv\nF9YArMsORm42x8rJZkv7fdP7OD5Wnyc9lowxlZ8KXZcMDUsRqdT6UkQul6vz0QdYAerk26sNwhj7\nKYqihECn2ZXoiMI1qwvXvdhHYxjP9DAw/cQ0Nz58A9t/PkHCasO/+/3v957iT6WcOO34u/wJsP7N\nxjbQYDCKooRAp9kHwbCNtPpdn26n7qEr1z2/eOsAc+c/yuon/ppLzKVctHhRayGv1sVvrb4mCp2f\ny996qVcbdGDspyiK4oeKeVR4uVW9732we3fwd3olttVUobX1uf76SM9hH7FJ35Em9tEY6TvS2BcS\nWbz1TvCKtw6QekmKwlWfJzOzt+Nj2zhGcDH3vbYVPV3+RJykNH5toD7piqJESZik58PyesMb3tBt\nnvfeYVnGOLJZ/xIxplhs3r9YNCaZrN83mfTetxNSKe/6pFKRHL74RNEk80nDray/kvmkKT4RUf27\nqJf1ScvIrWKsT1rd1yeVMkUwSTDUvJJgijW/VbFYNJZlGRExlmXVbfPE736xrO7qqyjKhgJ4zITQ\nR10zj4qgsKVea8a9Dv3pF9QF/OvZBuk70pSONdff2maxePNi8JdHKda4bZO+5ppAt7ROj6tr5oqi\ntELXzPtNkPuU19TpqE6zuulMl37mnTxl6ViL+ncQ5W2gZDL4XdFSN79VFN4EiqIoLirmUZHP+4+G\nvYQ+wtCf5bLN3FyagwdjzM2lnVzbtRHXavErD0N1HX55mclj3rsEGaEBgzH86tI2YdInFr2f4Vto\nBhinX1GUjYWKeVRkMnDjjc2C7mfNHYXlN46Qz//gelZXS4BhdbXE/A+up1y4qjkm+tgY3HlnW8ev\nI5dbDwqTn4VkQ3yY5FiS/HSL+vuNZkulrowAm4zxjrjHiWAmwNPAzS92vaIoygBQMY+SvXvh858P\nN3Ua0TTrwpE9VGL1qlqJnWDh7C86MdFrj3/PPd2N/mqEOHMECg+A9TMQ47wXriyQuaDF8YNGs62E\n1meEXY38VjpWwmAoHSuRfSDrCHoEMwF1Me1FwsWuVxRF6SNqADfs2Dbs2XParzmVckbXrpAcPCDe\n+bwNXHJpxL+tn9EedBfTPcyxAgzG0s/m/I3xfm+p43jv6/HWl5aYnJwkn8+rgCuK0lfUAG4jULNG\nvc7yMlx33fqoNOGTxdOvvCvyeSd1aSNjY53FdPfDayo+YITtZ3S3dGypY9uEunjrxlAqlchms/X5\nyxVFUYYEFfNBEcYoq2aNuo6TJ9eniae+kiJWH3ac2HGnPHIyGSd1aWM603an76uGX36C7iW0Adb/\nfkZ3k9smO7ZN8Iy3vrJSn79cURRlSFAxHwRhjbKCXJ/cbRNX3snO/zRG4ihQgcRR2Pmfxpi4sgtD\ntyAyGScmeTXMyXPPdb4O347QBoywvSK/rRvjdWib4Od21pU7mqIoSq8IE1lmWF5DHQGuHfyif8Xj\n9RHg/PZrjBRWLDqfRZz3qKLIdUKxWB99LpUKrk/YureImBd15DfLsuoivlVflkZoUxSljxAyAtzA\nBbqd1yiI+dGjRXP4sGUOHBBz+LBljh71EBURf5GuDelaLBozPt68TyzmiGSvxLvTzoFffcfGoqlj\nHzstxWLRJJPJ+hCuyWTrMK2KoigRomI+AI4eLZpDh5LmwAHWX4cOJZsFPWjE7TXqrh3pnnFGs2CG\njekeRgy7iRkf5rpGSAzbjreuKIoSMWHFXF3TImRuLu0Gb6knkbC46KLF0wWt3LOC3KY6jekeNhZ4\nNzHjg+LTB51TURRF8URd0wbA6qq3cVRTedUoKx73PlC7cd6DyquEDZ7STcz4MOFNow7dOmw55BVF\nUQaAinmEJBLeYuZZnsnAvfe27zblJ5ixWLCghRXpdv2ya8X0hRdgyxbv/cLUpV1GLWmLoihKj1Ax\nj5CpqTyxWL04x2JJpqZ8xDmM21TjyPOKK5o7AABra8GCFlak23EXaxTT5WXnOs480/tcrerSLoNI\n2qIoijKMhFlYH5bXsBvAGRPSmj0sfsZoMzOnDdni8dZGdEHH8jOCC2M17mfwVj13N8Z0YfDzChCJ\n5viKoigDBrVm3wC0Ektj2hO0qF27wpzb75xR1CVM+3igVuqKoowKKuYbgTBi2aGgRUIn5250tetm\nxN7ByF/9xxVFGSXCirmumQ8zYda5I8qL3hHtnru6xl6bOKZKJ2vdHYRq1ZjriqJsRNTPfJjx8A23\n3zDGA78zztutF5lIwFo8xQXPXcXE7z3oWIlPuslF+uXHbbs5w8OcOyiFKoRKS9otsVgMr3teRKj0\n+NyKoijtEqmfuYi8RUTmReQpEfmwx3ZLRGZF5AkROSgir6zZ9gkR+b77ek9NuYhIXkR+JCJPishN\nYS9u09Aw8rQvSWF/sMIHXvMi5211tG9LZZkfnPNZyv897wjh4mJ/A7JUM6CFOXcrl7SorNwDfM8n\nfc7hV64oijIKtBRzEYkDnwYuB84H3isi5zfsdjvwF8aYC4GPAR93v/tW4FeA1wO/Bvy+iJztfucD\nwKuAf2GMeR1wX9dXsxGpEcvc287k/a9eY2tDrJkYJ1hYCDlNPMggK0GCGdXSQAvf83w+T7JhaSCZ\nTJLvx7KEoihKjwgzMn8j8JQxZsEYcwJHdN/WsM/5wKz794Ga7ecDh4wxp4wxLwKPA29xt80AHzPG\nVACMMc90fhmbg6VjS+xIeG/ziz5Xh21Ttq9j7rYSBx8xzN1Womxf11tBbwwqMz7evE8qFV2I1wbf\n8/I0zH12hYMvv4a5uTSXXQaFQgHLshARLMuiUCiQ0fCyiqKMMGHE/BXAT2o+P+2W1fI48C7373cA\nZ4lIyi2/XESSIrIduBRnNA7wi8B7ROQxEXlIRH6p04vYLExum+SZVe9tftHnaik/sIf5m06yeh4Q\ng9XzYP6mk5Qf2BNtRat4BZUxxhHvqsFasdhdTvRGaqbyy9Mwf4tznQisrpaYn89y2WWwuLhIpVJh\ncXFRhVxRlJEnjJiLR1mjBdEtwMUi8l3gYuDvgVPGmK8DDwKHgf3AHHDK/U4COO4u7P85cLfnyUWy\nruA/9uyzz4ao7sYlP53nL5bGOL5WX15hvCnKnH3EJn1HmthHY6TvSGMfsVl4+zKVrQ3f3QoLb/ew\nLo8CrwhtJ086EeJ6tb5fM5W/cAPN11tZCb8koSiKMiKEEfOnOT2aBngl8NPaHYwxPzXGvNMY86+A\nnFt2zH3PG2Neb4z5DZyOwY9rjvtl9+//AlzodXJjTMEYs8sYs+vcc88NeVkbk8wFGTIX3cPnfpLi\n6HFnkHsqluKXX3c3ExOnRdE+YpN9IEvpWAmDoXSsRPaBLMcnvI+7uoPerKN3k7Sl07X9Gne51R3e\nu4RaklAURRkhQmTF4NvAL4nIq3FG3FcDv127gzuF/ry7/v0HuKNs13jupcaYZRG5EEewv+5+7SvA\nm9x9LwZ+1P3ljD72EZvcbI6lY0tMbpskP50nc8Fpoc5ckKn77EVuNsfKyQZf6pMrLK/G2L612f0q\n8SynXcaqBmPQ/ah5ctLbFa2V5XijS147dapuz+VIPFNyptgbCLMkoSiKMkq0HJkbY04Bvwv8FfAk\n8EVjzA9E5GMi8lvubpcA8yLyI2ACqM75jgHfEJEfAgXgGvd4ALcB7xKRIzjW7zdEdE0ji9+I2j7S\n3kh56Zj3yPMzCxVilXoDtNhxmPrzhh1XVmDPnu5H650GtOk2gYrrATB1abG9xDeKoigjigaNGSLS\nd6QpHWseyVrbLBZvXgx9nO3/YTvL/9y8Dp56SYofXHsnCws5VleXSCQmmfqjEhOzHgdpJJnszOK8\nnaAyVWIxZw2hkQ6CypTLdv31TuXrliQURVGGmbBBY1TMh4jYR2OYJttCEITKR8KL2Jn/95m8ePLF\npvIzxs7ghT98ob6wVVS2WizLMVrrNX516tf5FUVRhoRII8Ap/WFym090Mp9yP7yE3LfcayrcjzCG\na1EwyHjziqIoI4iK+RCRn86THGuITjaWJD/dQxHzSlaSSjXtZgPpWIxYLEY6ncbuZaCZDhKoKIqi\nbGZ0mn3IaGXNHga/NfOYOElGWh63wZrcBrJArUlaMpnUyGmKoig9RtfMNzH2EZvrv3o9J9ZO+O6T\nHEtSuLLQLOhVg7VSCeJxWFsjHY9TWltrOoZlWSzqGraiKErP0DXzAWLbNul0uj9T0h5kLshw99vu\nxtpmIQhxiTfts3Jyhdxsg6tXbfhVgLU1SCZZ8hBygKVSCbZsgd27o74ERVEUpQ10ZB4xtm2TzWZZ\nqfGTHvSUdGgreR8rct+RObBY/TAzA3v3RlFdRVEUxUVH5gMil8vVCTnAysoKubABT3pAaCt5H2v1\n/Npac9pQTkcGAhwDtR5SLtvMzaU5eDDG3Fyacrm/sx2KoijDjIp5xCz5CKJfeT9oZSW/npTljw3p\nm8G+oP77GTdNqGVZCM6IvADUzTP4TMVHQblsMz+fZXW1BJj17Gcq6IqiKA4q5hEz6RN33K+8H2Qu\nyFC4srC+hm5ts9aN3+pCyAqUXgrZK2sE3fXvzmQyTtrQeJxFGoQcHGO5HrGwkKNSqZ/t0OxniqIo\np1Exj5h8Pt88JZ1Mkh9wwJPMBRkWb16k8pEKizcvrluxeyZlGYfcNN7+3dWEJ41Uy20btm93/MNF\nnL+7NAD0y3Km2c8URVEcVMwjJpPJnJ6SFsFyp6iH1R/bLylL6aUQu26J9LO5+kQve/c6xm7VkXg8\nftr4zbbh+uthucbHfXkZrruuK0H3y3Km2c8URVEc1Jp9k+OX3KUWX5/0poOl/eO8dxFXvbpmXjvV\nHosl2bmzoElTFEXZ0Kg1uxIKL+O4RlZOrpC7f0/rgwUZ+XVhADgxkWHnzgKJhAUIiYSlQq4oilLD\nlkFXQBkstWvnQSP0pZPLzlR50HLB5KT/yLxLA8CJiYyKt6Ioig86MlfWjeME8d1n8hhOmNcg8nkY\nH28uHxtzttm2MxUfiznvfY6MpyiKslFRMVfWOecl53iWi4H8LK2nyjMZuPvu+qxrqRTcc4/zdzVU\nrDHOezZbL+gq9oqiKB2hYq4ATuCYn5/4efMGAzd+CzJHaDlVbts26VyO2PPPk7Ys7GIRnnvOEflc\nbj0L2zorK6dH+7Vx4f3EXlEURfFExVwBnDVzryxrqRXY+xDrwWP8qMakL5VKGGMolUpks9nTSWb8\nRvWl0ulMbUFiryiKoviirmkKEJCMxUDlc5Yj5AHGb+l0mpKH8dt6mtQgt7VkslnI1ysgUKl4b1MU\nRdngqGua0ha+yVhe6vqHtwh60zImfT7viLYXKyv+4WAHGAZXURRlVFAxV4DWyVha0TImfSYTnFnN\nzZ1eX4HgqX1FURTFQcVcAYKTsQTiWqDnSyWSUu/a1hSTPpNxIsF5UY0Db1nO1LpXXHhFURTFE10z\nVzqnaoHurnfbQE6EJWOYtCzybra1oO8AzghchVtRFKWJsGvmGgFO6ZwGC/QMkDEmOA57VbBzOcfC\nfXKypXGdoiiKEoyOzJXOicUcn/BG1AJdURQlEtSaXek9fpbmaoGuKIrSV1TMlc7xcjdTC3RFUZS+\no2KudE7V3Uwt0BVFUQaKGsAp3ZHJqHgriqIMGB2ZK4qiKMqIo2KuKIqiKCOOivlGZPdu2LLFWcfe\nssX5rCiKomxYVMw3ErYNZ54Jd93lxDoH5/2uu/or6G6IV2Ix511zkiuKovQUFfONgm3D9dfDiy96\nbw9KchJ1PbJZJ92pMc57NquCriiK0kM0AtxGYft2WF4O3qcfv7Vf3vKgEK+KoiiKJxoBbrPRSsj9\n8oVHjU8kTXwJAAAMgklEQVRec99yRVEUpWtUzDcL2Wx/zqMhXhVFUfqOivlGIZXy3zYzA3v3NpdH\naahWPVap5FjR16IhXhVFUXqKivlG4c47YWysvmxsDIpFfyGPylCt9ljgHK8q6BriVVEUpeeoAdxG\nwrbD5wmP0lBNjd4URVF6QlgDOBXzzUqUucg1r7miKEpPUGt2xRPbtkmn08SMIQ00Tap3YqimRm+K\noigDRcV8E2HbNtlsllKphAFKQJYaQe/UUE3zmiuKogyUUGIuIm8RkXkReUpEPuyx3RKRWRF5QkQO\nisgra7Z9QkS+777e4/HdT4nIC91dhhKGXC7HyspKXdkKkIPuDNU0r7miKMpAaZnPXETiwKeB3wCe\nBr4tIvcbY35Ys9vtwF8YY+4VkTcBHwfeJyJvBX4FeD2QAA6JyEPGmH9yj70LeGmkV6T4suQTuGVJ\npHtDNc1rriiKMjDCjMzfCDxljFkwxpwA7gPe1rDP+cCs+/eBmu3nA4eMMaeMMS8CjwNvgfVOwn8E\n/o/uLkEJy6TPGrZfuaIoijIahBHzVwA/qfn8tFtWy+PAu9y/3wGcJSIpt/xyEUmKyHbgUuBV7n6/\nC9xvjPmHoJOLSFZEHhORx5599tkQ1VX8yOfzJBvWtpPJJHld21YURRlpwoi5eJQ1+iHdAlwsIt8F\nLgb+HjhljPk68CBwGNgPzAGnROTlwL8GPtXq5MaYgjFmlzFm17nnnhuiuhuIiPOSZzIZCoUClmUh\nIliWRaFQIFMzPb5u7R6LkU6nsTXbmaIoytDT0s9cRC4CbjXG/Kb7+Q8AjDEf99n/TOBvjDGv9Nj2\nBaCI00H4LHDc3TQJLBhjXhNUl03lZ757t5OHvBG/0KwRULV2rzWSSyaTTYKvKIqi9IfIgsaIyBbg\nR8A0zoj728BvG2N+ULPPduB5Y0xFRPLAmjHmj9118ZcaY5ZF5ELgC8DrjTGnGs7xgjHmzFaV3VRi\nvmULrK01l8fjcOpUc3kEpNNpSh6R3CzLYlEjuSmKovSdsGLe0prdGHNKRH4X+CsgDtxtjPmBiHwM\neMwYcz9wCfBxETHAXwP/u/v1MeAb4sTp/ifgmkYhV3zwEvKg8gjwtXbX9KWKoihDjYZzHVaGaWQe\nj7N4773qeqYoitJnNJzrqOOXf7yHeck9rd2B/Npa5xnVFEVRlJ6jYj6s7N3rGLvF487neLynxm9Q\nY+0ejyOABRSADMDKipORTVEURRk6dJpdaUazoCmKogwFOs2udI5mQVMURRkpVMwHjW1DOu2MhtPp\n3q5L155r+3bn5XVezYKmKIoyUrR0TVN6iG07hmXVIC2l0mkDt6gtxxvPtbx8elvjeavnzuVgackZ\nkefzas2uKIoypOia+SBJpx0hbcSyus9iFvZcvT6voiiK0jG6Zj4K+AVj6UWQljDH1OAwiqIoI4mK\n+SDpp6FZmGOqgZuiKMpIomI+SPppaOZ1rn6cV1EURek5KuZ9oFy2mZtLc/BgjLm5NOWyazmeyUCh\n4KxVizjvhUJvDM0az5VKOa9en1dRFEXpOWoA12PKZZv5+SyVyum0orFYkp07C0xMhBBP21arckVR\nlE2KGsANCQsLuTohB6hUVlhYCBEatepOVio5EdmqLmQaI11RFEWpQcW8x6yueluI+5XXkcud9guv\nojHSFUVRlAZUzHtMIuFtIe5XXkc/XdcURVGUkUXFvMdMTeWJxeqtyGOxJFNTISzHNUa6oiiKEgIV\n8x4zMZFh584CiYQFCImEFd74TWOkK4qiKCHQ2Ox9YGIiE068G9EY6YqiKEoIdGQ+7GQyTrz0SsV5\nrxXyfmZcUxRFUYYWHZmPKv3MuKYoiqIMNToyH1XUbU1RFEVxUTEfVdRtTVEURXFRMR9VfNzTyled\n4x0HXlEURdmwqJiPKh5ua+XLx/ibf3OM1dUSYFhdLfHkk9fwve9dNpg6KoqiKH1BxXxU8ci49uN/\nl8DIqaZdf/azWX70o90DqKSiKIrSD1TMR5kGt7VT8Rd8d/3pTwv9q5eiKIrSV1TMNw1rg66AoiiK\n0iNUzDcQ8XgqaGvf6qEoiqL0FxXzDcRrX3snIJ7bXv7ybH8royiKovQNFfMNxMREhte97vOInFFT\nGuPlL5/hta/dO7B6KYqiKL1Fw7luMDpO6qIoiqKMLDoyVxRFUZQRR8VcURRFUUYcFXNFURRFGXFU\nzBVFURRlxFExVxRFUZQRR8VcURRFUUYcFXNFURRFGXFUzBVFURRlxFExVxRFUZQRR8VcURRFUUYc\nFXNFURRFGXFUzBVFURRlxFExVxRFUZQRR8VcURRFUUYcFXNFURRFGXHEGDPoOoRGRJ4FSgOuxnbg\nuQHXYRTRdusMbbfO0bbrDG23zulF21nGmHNb7TRSYj4MiMhjxphdg67HqKHt1hnabp2jbdcZ2m6d\nM8i202l2RVEURRlxVMwVRVEUZcRRMW+fwqArMKJou3WGtlvnaNt1hrZb5wys7XTNXFEURVFGHB2Z\nK4qiKMqIo2IOiMi/FpEfiEhFRHwtEUXkLSIyLyJPiciHa8pfLSL/XUR+LCJ/KSLjDd97t4iYoGOP\nKr1qOxH5PRH5oYg8ISKzImL143r6RQ/bLeF+fsrdnu791fQPETlHRB52r/thEXmZz36fEJHvu6/3\n1JS/SUT+h1t+r4hsccu3icgDIvK4+7tc169r6he9ajt32yUi8j237Q7143r6RS/bzd3+qyKyJiLv\n7qqixphN/wJeB+wEDgK7fPaJA38LTAHjwOPA+e62LwJXu39/Bpip+d5ZwF8D3/Q79ii/etV2wKVA\n0v17BvjLQV/riLTbbuAz7t9Xb8B2+w/Ah92/Pwx8wmOftwIPA1uAM4DHgLNxBi8/AV7r7vcx4IPu\n339YPRZwLvA8MD7o6x2Rtnsp8ENg0v28Y9DXOgrt5n6OA48CDwLv7qaeOjIHjDFPGmPmW+z2RuAp\nY8yCMeYEcB/wNhER4E3Af3b3uxd4e833/i+cm+F4xNUeCnrVdsaYA8aYFbf8m8Aro6/94OjhPfc2\n9zPu9ml3/41C7fU1/q9VOR84ZIw5ZYx5EacT9BYgBawaY37k7vcw8C73bwOc5bbVmThifqo3lzAw\netV2vw38v8aYJQBjzDM9qv+g6FW7Afxb4MtA122mYh6eV+D0sKo87ZalgJ8ZY041lCMi/wp4lTHm\n/+tnRYeQttuugQ8CD/W0hsNJJ+22/h13+zF3/43ChDHmHwDc9x0e+zwOXC4iSRHZjjPL8yqcyFxj\nNcsa73bLAf4fnNmSnwJHgD3GmErvLmMg9KrtXgu8TEQOish3ROT9Pb2K/tOTdhORVwDvwJlZ65ot\nrXfZGIjII8B5HptyxpivhjmER5nxKxeRGPBJ4AOhKzmk9LvtGs59DbALuDjEeYaKAbVbyzYddoLa\nLcz3jTFfF5FfBQ4DzwJzwCljjBGRq4FPikgC+DqnR9+/CXwPZ8bjF4GHReQbxph/6u5q+suA2m4L\n8AZgGngJMCci36wZjQ49A2q3O4B/b4xZi2LybNOIuTHmsi4P8TSne6LgTPv+FKfn9VIR2eKOhKrl\nZwH/Ejjo/lDnAfeLyG8ZYx7rsi59ZQBtB4CIXIbzz3SxMWa1yzr0nQG1W/U7T7uGNttwpoxHhqB2\nE5GyiPyCMeYfROQX8JmeNMbkgbz7nS8AP3bL54Bfd8vfjDOqBLgOuM04C5lPicjfAf8C+FY0V9Uf\nBtR2TwPPudPLL4rIXwP/EzAyYj6gdtsF3Ofqw3bgChE5ZYz5SifXoNPs4fk28EviWBGP4xgX3e/+\n8x/AmT4BuBb4qjHmmDFmuzEmbYxJ46z7jpyQR0RbbQfrSxR/htNmG20NLixttxtwv/sZd/uj7v4b\nhdrrq73udUQkLiIp9+8LgQtxRkSIyA73PQH8e05PcS7hjCwRkQkc48SFnl3FYOhV230V+HUR2SIi\nSeDXgCd7eB39piftZox5dY0+/Gdgd6dCjnvATf/CWbd4GlgFysBfueUvBx6s2e8KnN7m3+JMlVbL\np3B68E8BXwISHuc4yMa0Zu9J2wGPuMf7nvu6f9DXOiLtttX9/JS7fWrQ1xpxu6WAWZxRzyxwjlu+\nC9hX0wY/dF/fBF5f8/3/iCM088DNNeUvx3n4HgG+D1wz6GsdlbZzt/2++53vN24b9Vcv261mn8/R\npTW7RoBTFEVRlBFHp9kVRVEUZcRRMVcURVGUEUfFXFEURVFGHBVzRVEURRlxVMwVRVEUZcRRMVcU\nRVGUEUfFXFEURVFGHBVzRVEURRlx/n/JfDP+qCmi+AAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"k = 7\n",
"km = KMeans(n_clusters=k)\n",
"km.fit(coords)\n",
"labels = km.labels_\n",
"centroids = km.cluster_centers_\n",
"colors = ['r', 'g', 'b', 'y', 'c', 'm', 'k']\n",
"plt.figure(figsize=(8,8))\n",
"plt.title('Clustered Locations')\n",
"for j in range(k):\n",
" current_coords = np.asarray(coords)[np.where(labels==j)]\n",
" j_longs = current_coords[:,1]\n",
" j_lats = current_coords[:,2]\n",
" plt.plot(j_longs, j_lats, '%so' % colors[j], label=j)\n",
"plt.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"7 is interesting indeed. Mississauga and Toronto are separated, and all but one of the Costco stations are grouped together regardless of distance."
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"cluster_df['label'] = km.labels_"
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"label\n",
"2 122.913025\n",
"0 122.780663\n",
"6 120.150000\n",
"3 118.368056\n",
"1 116.464407\n",
"4 113.891667\n",
"5 111.114286\n",
"Name: mean, dtype: float64"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cluster_df.groupby('label')['type_A'].describe().sort_values('mean',ascending=False)['mean']"
]
},
{
"cell_type": "code",
"execution_count": 74,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" longitude | \n",
" latitude | \n",
" type_A_norm | \n",
" longitude_norm | \n",
" latitude_norm | \n",
" label | \n",
"
\n",
" \n",
" \n",
" \n",
" 1411 Warden Ave , Toronto - East, Ontario | \n",
" Costco | \n",
" 110.9 | \n",
" -79.297736 | \n",
" 43.759429 | \n",
" 0.917363 | \n",
" -0.997628 | \n",
" 1.000612 | \n",
" 5 | \n",
"
\n",
" \n",
" 150 Kingston Rd E , Ajax, Ontario | \n",
" Costco | \n",
" 110.9 | \n",
" -79.018454 | \n",
" 43.862950 | \n",
" 0.917363 | \n",
" -0.994114 | \n",
" 1.002979 | \n",
" 5 | \n",
"
\n",
" \n",
" 1570 Dundas St E , Mississauga, Ontario | \n",
" Costco | \n",
" 110.9 | \n",
" -79.575776 | \n",
" 43.610304 | \n",
" 0.917363 | \n",
" -1.001126 | \n",
" 0.997202 | \n",
" 5 | \n",
"
\n",
" \n",
" 55 New Huntington Rd , Vaughan, Ontario | \n",
" Costco | \n",
" 110.9 | \n",
" -79.639189 | \n",
" 43.769937 | \n",
" 0.917363 | \n",
" -1.001924 | \n",
" 1.000852 | \n",
" 5 | \n",
"
\n",
" \n",
" 71 Colossus Dr , Vaughan, Ontario | \n",
" Costco | \n",
" 110.9 | \n",
" -79.541590 | \n",
" 43.785631 | \n",
" 0.917363 | \n",
" -1.000696 | \n",
" 1.001211 | \n",
" 5 | \n",
"
\n",
" \n",
" 935 Liverpool Rd , Pickering, Ontario | \n",
" Shell | \n",
" 112.4 | \n",
" -79.086912 | \n",
" 43.828548 | \n",
" 0.929771 | \n",
" -0.994976 | \n",
" 1.002192 | \n",
" 5 | \n",
"
\n",
" \n",
" Taymall Ave , Toronto - West, Ontario | \n",
" Costco | \n",
" 110.9 | \n",
" -79.506616 | \n",
" 43.623887 | \n",
" 0.917363 | \n",
" -1.000256 | \n",
" 0.997512 | \n",
" 5 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A longitude \\\n",
"1411 Warden Ave , Toronto - East, Ontario Costco 110.9 -79.297736 \n",
"150 Kingston Rd E , Ajax, Ontario Costco 110.9 -79.018454 \n",
"1570 Dundas St E , Mississauga, Ontario Costco 110.9 -79.575776 \n",
"55 New Huntington Rd , Vaughan, Ontario Costco 110.9 -79.639189 \n",
"71 Colossus Dr , Vaughan, Ontario Costco 110.9 -79.541590 \n",
"935 Liverpool Rd , Pickering, Ontario Shell 112.4 -79.086912 \n",
"Taymall Ave , Toronto - West, Ontario Costco 110.9 -79.506616 \n",
"\n",
" latitude type_A_norm \\\n",
"1411 Warden Ave , Toronto - East, Ontario 43.759429 0.917363 \n",
"150 Kingston Rd E , Ajax, Ontario 43.862950 0.917363 \n",
"1570 Dundas St E , Mississauga, Ontario 43.610304 0.917363 \n",
"55 New Huntington Rd , Vaughan, Ontario 43.769937 0.917363 \n",
"71 Colossus Dr , Vaughan, Ontario 43.785631 0.917363 \n",
"935 Liverpool Rd , Pickering, Ontario 43.828548 0.929771 \n",
"Taymall Ave , Toronto - West, Ontario 43.623887 0.917363 \n",
"\n",
" longitude_norm latitude_norm \\\n",
"1411 Warden Ave , Toronto - East, Ontario -0.997628 1.000612 \n",
"150 Kingston Rd E , Ajax, Ontario -0.994114 1.002979 \n",
"1570 Dundas St E , Mississauga, Ontario -1.001126 0.997202 \n",
"55 New Huntington Rd , Vaughan, Ontario -1.001924 1.000852 \n",
"71 Colossus Dr , Vaughan, Ontario -1.000696 1.001211 \n",
"935 Liverpool Rd , Pickering, Ontario -0.994976 1.002192 \n",
"Taymall Ave , Toronto - West, Ontario -1.000256 0.997512 \n",
"\n",
" label \n",
"1411 Warden Ave , Toronto - East, Ontario 5 \n",
"150 Kingston Rd E , Ajax, Ontario 5 \n",
"1570 Dundas St E , Mississauga, Ontario 5 \n",
"55 New Huntington Rd , Vaughan, Ontario 5 \n",
"71 Colossus Dr , Vaughan, Ontario 5 \n",
"935 Liverpool Rd , Pickering, Ontario 5 \n",
"Taymall Ave , Toronto - West, Ontario 5 "
]
},
"execution_count": 74,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cluster_df[cluster_df['label']==5]"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"### Folium\n",
"\n",
"Let's put it on a Toronto map with Folium.\n",
"\n",
"Then run our function for getting the nearest gas stations within a specified distance in kilometers."
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import branca\n",
"import folium\n",
"from folium.plugins import MarkerCluster"
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" name | \n",
" type_A | \n",
" longitude | \n",
" latitude | \n",
" type_A_norm | \n",
" longitude_norm | \n",
" latitude_norm | \n",
" label | \n",
"
\n",
" \n",
" \n",
" \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario | \n",
" Ultramar | \n",
" 115.6 | \n",
" -79.282378 | \n",
" 43.976346 | \n",
" 0.956242 | \n",
" -0.997435 | \n",
" 1.005572 | \n",
" 1 | \n",
"
\n",
" \n",
" 1610 Keele St , Toronto - Central, Ontario | \n",
" Shell | \n",
" 122.9 | \n",
" -79.471988 | \n",
" 43.682320 | \n",
" 1.016627 | \n",
" -0.999820 | \n",
" 0.998849 | \n",
" 0 | \n",
"
\n",
" \n",
" 6897 Finch Ave W , Toronto - West, Ontario | \n",
" Esso | \n",
" 121.9 | \n",
" -79.617364 | \n",
" 43.734632 | \n",
" 1.008355 | \n",
" -1.001649 | \n",
" 1.000045 | \n",
" 0 | \n",
"
\n",
" \n",
" 1 Harwood Ave S , Ajax, Ontario | \n",
" Pioneer | \n",
" 122.6 | \n",
" -79.025100 | \n",
" 43.861027 | \n",
" 1.014146 | \n",
" -0.994198 | \n",
" 1.002935 | \n",
" 2 | \n",
"
\n",
" \n",
" 1 Thornhill Woods Dr , Vaughan, Ontario | \n",
" Esso | \n",
" 122.9 | \n",
" -79.463886 | \n",
" 43.826660 | \n",
" 1.016627 | \n",
" -0.999718 | \n",
" 1.002149 | \n",
" 2 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name type_A \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 \n",
" 1610 Keele St , Toronto - Central, Ontario Shell 122.9 \n",
" 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 \n",
"1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 \n",
"\n",
" longitude latitude \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n",
" 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n",
"1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n",
"\n",
" type_A_norm \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0.956242 \n",
" 1610 Keele St , Toronto - Central, Ontario 1.016627 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 1.008355 \n",
"1 Harwood Ave S , Ajax, Ontario 1.014146 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 1.016627 \n",
"\n",
" longitude_norm \\\n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -0.997435 \n",
" 1610 Keele St , Toronto - Central, Ontario -0.999820 \n",
" 6897 Finch Ave W , Toronto - West, Ontario -1.001649 \n",
"1 Harwood Ave S , Ajax, Ontario -0.994198 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario -0.999718 \n",
"\n",
" latitude_norm label \n",
" 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 1.005572 1 \n",
" 1610 Keele St , Toronto - Central, Ontario 0.998849 0 \n",
" 6897 Finch Ave W , Toronto - West, Ontario 1.000045 0 \n",
"1 Harwood Ave S , Ajax, Ontario 1.002935 2 \n",
"1 Thornhill Woods Dr , Vaughan, Ontario 1.002149 2 "
]
},
"execution_count": 71,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cluster_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"folmap(cluster_df,zoom=10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Enter the address and distance threshold to see nearby stations."
]
},
{
"cell_type": "code",
"execution_count": 73,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Getting 80 Bloor St W, Toronto, ON M5S 2V1\n",
"Got 80 Bloor St W, Toronto, ON M5S 2V1, Canada \n",
"\n",
"Gas stations within range: \n",
" name type_A\n",
"835 Yonge St , Toronto - Central, Ontario Canadian Tire 119.3\n",
"1077 Yonge St , Toronto - Central, Ontario Shell 122.9\n",
"150 Dupont St , Toronto - Central, Ontario Esso 122.9\n",
"241 Church St , Toronto - South, Ontario Esso 122.9\n",
"505 Jarvis St , Toronto - South, Ontario Petro-Canada 122.9\n",
"581 Parliament St , Toronto - South, Ontario Esso 122.9\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"enter_address(cluster_df,'80 Bloor St W, Toronto, ON M5S 2V1',2)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}