notebooks/2018-09-16-hdb-feature-engineering-ii.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "<script>\n",
    "  function code_toggle() {\n",
    "    if (code_shown){\n",
    "      $('div.input').hide('500');\n",
    "      $('#toggleButton').val('Show Code')\n",
    "    } else {\n",
    "      $('div.input').show('500');\n",
    "      $('#toggleButton').val('Hide Code')\n",
    "    }\n",
    "    code_shown = !code_shown\n",
    "  }\n",
    "\n",
    "  $( document ).ready(function(){\n",
    "    code_shown=false;\n",
    "    $('div.input').hide()\n",
    "  });\n",
    "</script>\n",
    "<div align=\"right\">\n",
    "    <form action=\"javascript:code_toggle()\"><input type=\"submit\" id=\"toggleButton\" value=\"Show Code\"></form>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "In my [second post on HDB resale flat prices](https://dataandstuff.wordpress.com/2017/09/09/resale-flats-and-clusters/), I attempted to price resale flats in Jurong West by creating clusters of flats. The big idea behind this methodology was the realisation that it is impossible to explicitly account for all qualitative reasons why home buyers would want to buy a flat in a specific area. It could be because of good schools, proximity to amenities like transport centers and malls, or the liveliness of the areas. Furthermore, each factor would have different importance to different people. Modeling this would be a nightmare. Hence, I chose just two features to develop a spatial representation of these preferences: (1) latitude and (2) longitude. In fact, I showed that clusters of nearby flats shared a statistically meaningful relationship with resale prices in those clusters.  \n",
    "  \n",
    "In that post, I computed clusters for Jurong West. This time, I apply that methodology to all 26 towns in the dataset. The objective is to develop a set of clusters for **each town** to be used as categorical features in our model of HDB resale flat prices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules\n",
    "import gmaps\n",
    "import json\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn.apionly as sns\n",
    "from scipy.spatial import ConvexHull\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "import tabulate\n",
    "import urllib.request as ur\n",
    "import warnings\n",
    "\n",
    "# Settings\n",
    "%matplotlib inline\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# Colours\n",
    "def get_cols():\n",
    "    \n",
    "    print('[Colours]:')\n",
    "    print('Orange:     #ff9966')\n",
    "    print('Navy Blue:  #133056')\n",
    "    print('Light Blue: #b1ceeb')\n",
    "    print('Green:      #6fceb0')\n",
    "    print('Red:        #f85b74')\n",
    "\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Modify settings\n",
    "mpl.rcParams['axes.grid'] = True\n",
    "mpl.rcParams['axes.grid.axis'] = 'y'\n",
    "mpl.rcParams['grid.color'] = '#e8e8e8'\n",
    "mpl.rcParams['axes.spines.right'] = False\n",
    "mpl.rcParams['axes.spines.top'] = False\n",
    "mpl.rcParams['xtick.color'] = '#494949'\n",
    "mpl.rcParams['xtick.labelsize'] = 12\n",
    "mpl.rcParams['ytick.color'] = '#494949'\n",
    "mpl.rcParams['ytick.labelsize'] = 12\n",
    "mpl.rcParams['axes.edgecolor'] = '#494949'\n",
    "mpl.rcParams['axes.labelsize'] = 15\n",
    "mpl.rcParams['axes.labelpad'] = 15\n",
    "mpl.rcParams['axes.labelcolor'] = '#494949'\n",
    "mpl.rcParams['axes.axisbelow'] = True\n",
    "mpl.rcParams['figure.titlesize'] = 20\n",
    "mpl.rcParams['figure.titleweight'] = 'bold'\n",
    "mpl.rcParams['font.family'] = 'sans-serif'\n",
    "mpl.rcParams['font.sans-serif'] = 'Raleway'\n",
    "mpl.rcParams['scatter.marker'] = 'h'\n",
    "\n",
    "# Colours\n",
    "def get_cols():\n",
    "    \n",
    "    print('[Colours]:')\n",
    "    print('Orange:     #ff9966')\n",
    "    print('Navy Blue:  #133056')\n",
    "    print('Light Blue: #b1ceeb')\n",
    "    print('Green:      #6fceb0')\n",
    "    print('Red:        #f85b74')\n",
    "\n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geocoding\n",
    "ArcGIS defines geocoding as \"the process of transforming a description of a location—such as a pair of coordinates, an address, or a name of a place—to a location on the earth's surface\". In the HDB resale flat dataset, we are provided with a block number and a street name. Combining these two features gives us an address that can be converted into latitude and longitude. If we think of latitude and longitude as the *y* and *x* axes on a graph, then each flat is simply a point on the graph, and groups of flats can be identified easily.  \n",
    "  \n",
    "How do we collect this data? Easy: we get it [HERE](https://developer.here.com/). HERE provides a host of location services like interactive maps, geocoding, traffic, tracking and routing. Its clients include Bing, Samsung, Audi, and Grab. Yes, Singapore's Grab. What's amazing is that it provides users with an API (see the link above) with up to **250,000 free requests**. That is an insanely large amount of requests, at least when compared to Google Places' 40,000 request limit.  \n",
    "  \n",
    "To begin, we load the HDB resale flat dataset and create two addresses:  \n",
    "  \n",
    "1. **Full Address:** To facilitating identification of the address. For example: *174 ANG MO KIO AVE 4*\n",
    "2. **Search Address:** This is the address that we are going to append to our query to HERE. For example: *174+ANG+MO+KIO+AVE+4+SINGAPORE*\n",
    "  \n",
    "Note that I extracted the unique addresses from the full list of search addresses we created. This is to save time and stick within the query limit, because we need not geocode the same address more than once. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Read data\n",
    "hdb = pd.read_csv('resale-flat-prices-based-on-registration-date-from-jan-2015-onwards.csv')\n",
    "\n",
    "# Create addresses\n",
    "hdb['full_address'] = hdb.block + ' ' + hdb.street_name\n",
    "hdb['search_address'] = hdb.block + '+' + hdb.street_name.str.replace(' ', '+') + '+SINGAPORE'\n",
    "\n",
    "# Extract search addresses\n",
    "all_adds = hdb.search_address.unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to configure two parameters: our individual app ID and app code. These can be found on your project page after you have created an account. We save these as variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set parameters\n",
    "APP_ID = '[YOUR HERE APP ID]'\n",
    "APP_CODE = '[YOUR HERE APP CODE]'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will loop through the unique list of addresses, use `urllib` to query the API, `json` to process the JSON response (into a Python dictionary), and a custom function to extract the data we need from the processed response. The function `get_loc` retrieves two sets of latitude and longitude: the display position and the navigation position. Based on some testing, I discovered that the **average** of these two sets of coordinates was more accurate than either of them individually. Hence, I chose to compute the average latitude and average longitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define function to extract average location\n",
    "# Takes a dictionary object\n",
    "# Returns the average of display position and navigation position in a dictionary\n",
    "def get_loc(result):\n",
    "    \n",
    "    # Output\n",
    "    output = dict()\n",
    "    \n",
    "    if len(result['Response']['View']) > 0:\n",
    "    \n",
    "        # Get display position lat/long\n",
    "        lat_dp = result['Response']['View'][0]['Result'][0]['Location']['DisplayPosition']['Latitude']\n",
    "        lon_dp = result['Response']['View'][0]['Result'][0]['Location']['DisplayPosition']['Longitude']\n",
    "\n",
    "        # Get navigation position lat/long\n",
    "        lat_np = result['Response']['View'][0]['Result'][0]['Location']['NavigationPosition'][0]['Latitude']\n",
    "        lon_np = result['Response']['View'][0]['Result'][0]['Location']['NavigationPosition'][0]['Longitude']\n",
    "    \n",
    "        # Configure output\n",
    "        output['lat'] = (lat_dp + lat_np) / 2\n",
    "        output['lon'] = (lon_dp + lon_np) / 2\n",
    "        \n",
    "    else:\n",
    "        \n",
    "        # Configure output\n",
    "        output['lat'] = np.nan\n",
    "        output['lon'] = np.nan\n",
    "    \n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the search addresses created, app ID and code configured, and the function defined, we are ready to geocode. Note that the output of `get_loc` for each search address was a dictionary, and all the dictionaries were saved into a list. This facilitated conversion into a Pandas dataframe. The code below will take approximately 5 minutes to run and 8.5k queries out of your HERE API limit. I've run the code and saved the results to a CSV file for quick loading."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise results\n",
    "all_latlon = []\n",
    "\n",
    "# Loop through to get lat lon\n",
    "for i in range(len(all_adds)):\n",
    "    \n",
    "    # Extract address\n",
    "    temp_add = all_adds[i]\n",
    "    \n",
    "    # Configure URL\n",
    "    temp_url = 'https://geocoder.api.here.com/6.2/geocode.json' + \\\n",
    "      '?app_id=' + APP_ID + \\\n",
    "      '&app_code=' + APP_CODE + \\\n",
    "      '&searchtext=' + temp_add\n",
    "    \n",
    "    # Pull data\n",
    "    temp_response = ur.urlopen(ur.Request(temp_url)).read()\n",
    "    temp_result = json.loads(temp_response)\n",
    "    \n",
    "    # Process data\n",
    "    temp_latlon = get_loc(temp_result)\n",
    "    \n",
    "    # Add address\n",
    "    temp_latlon['address'] = temp_add\n",
    "    \n",
    "    # Append\n",
    "    all_latlon.append(temp_latlon)\n",
    "    \n",
    "    # Update\n",
    "    print(str(i) + '. ', 'Getting data for: ' + str(temp_add))\n",
    "\n",
    "# Convert to data frame\n",
    "full_latlon = pd.DataFrame(all_latlon)\n",
    "\n",
    "# Save\n",
    "full_latlon.to_csv('latlon_data.csv', index = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a preview of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>address</th>\n",
       "      <th>lat</th>\n",
       "      <th>lon</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>174+ANG+MO+KIO+AVE+4+SINGAPORE</td>\n",
       "      <td>1.375270</td>\n",
       "      <td>103.837640</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>541+ANG+MO+KIO+AVE+10+SINGAPORE</td>\n",
       "      <td>1.374025</td>\n",
       "      <td>103.855695</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>163+ANG+MO+KIO+AVE+4+SINGAPORE</td>\n",
       "      <td>1.373885</td>\n",
       "      <td>103.838110</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>446+ANG+MO+KIO+AVE+10+SINGAPORE</td>\n",
       "      <td>1.367855</td>\n",
       "      <td>103.855395</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>557+ANG+MO+KIO+AVE+10+SINGAPORE</td>\n",
       "      <td>1.371540</td>\n",
       "      <td>103.857790</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           address       lat         lon\n",
       "0   174+ANG+MO+KIO+AVE+4+SINGAPORE  1.375270  103.837640\n",
       "1  541+ANG+MO+KIO+AVE+10+SINGAPORE  1.374025  103.855695\n",
       "2   163+ANG+MO+KIO+AVE+4+SINGAPORE  1.373885  103.838110\n",
       "3  446+ANG+MO+KIO+AVE+10+SINGAPORE  1.367855  103.855395\n",
       "4  557+ANG+MO+KIO+AVE+10+SINGAPORE  1.371540  103.857790"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load data\n",
    "map_latlon = pd.read_csv('latlon_data.csv')\n",
    "\n",
    "# View\n",
    "map_latlon.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use this data frame to map the search addresses in the main dataframe to a latitude and longitude pair. With that, we are done with geocoding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set index - I should have named the search address as search_address in the loop\n",
    "map_latlon = map_latlon.rename({'address': 'search_address'})\n",
    "map_latlon = map_latlon.set_index('address')\n",
    "\n",
    "# Separate maps\n",
    "map_lat = map_latlon['lat']\n",
    "map_lon = map_latlon['lon']\n",
    "\n",
    "# Map\n",
    "hdb['lat'] = hdb.search_address.map(map_lat)\n",
    "hdb['lon'] = hdb.search_address.map(map_lon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clustering\n",
    "Next, we develop clusters within each of the 26 towns using the resale flats' coordinates. Let's use Jurong West as an example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Jurong West\n",
    "First, we extract the coordinates of flats in Jurong West and scale the coordinates using the `MinMaxScaler`. This transforms each of the coordinates into features that range between 0 and 1 to ensure that longitude, which has a larger magnitude (100 vs. 1), does not affect the distance calculations in clustering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lat</th>\n",
       "      <th>lon</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.685821</td>\n",
       "      <td>0.624126</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.703477</td>\n",
       "      <td>0.616068</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.758060</td>\n",
       "      <td>0.652565</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.685821</td>\n",
       "      <td>0.624126</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.824356</td>\n",
       "      <td>0.927597</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        lat       lon\n",
       "0  0.685821  0.624126\n",
       "1  0.703477  0.616068\n",
       "2  0.758060  0.652565\n",
       "3  0.685821  0.624126\n",
       "4  0.824356  0.927597"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For Jurong West\n",
    "dat_jw = hdb[['lat', 'lon']][hdb.town == 'JURONG WEST']\n",
    "dat_jw = dat_jw.reset_index(drop = True)\n",
    "\n",
    "# Normalise\n",
    "mmscale = MinMaxScaler()\n",
    "mmscale.fit(dat_jw)\n",
    "dat_jw_scaled = pd.DataFrame(mmscale.transform(dat_jw), columns = ['lat', 'lon'])\n",
    "dat_jw_scaled.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we apply K-Means clustering on the dataset. Testing out different values of *k* (2 to 22) and using the elbow method to decide on an optimal *k*, we get *k* = 7."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up values of k\n",
    "all_k = np.arange(2, 23, 1)\n",
    "\n",
    "# Initialise results\n",
    "k_results = []\n",
    "\n",
    "# Loop through values of k\n",
    "for k in all_k:\n",
    "    \n",
    "    # Set up kmeans\n",
    "    km1 = KMeans(n_clusters = k, random_state = 123)\n",
    "\n",
    "    # Fit data\n",
    "    km1.fit(dat_jw_scaled)\n",
    "\n",
    "    # Score data\n",
    "    k_results.append(km1.inertia_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAngAAAIECAYAAACKd1OAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3X18lXX9x/HXxs2GjE0FN2CADNhQQQEr766Vgmj1zTK1xG40s0R/JVRWlnmHd/nzJksoNStvsvzJ1MzKb5oIars0TQUSvGELhIFsU1A2hI27/f64rgOHsY2z7XzPObt4Px+PPTjnOtd1fT7nbLrPvrdZLS0tiIiIiEh0ZKc7ARERERFJLhV4IiIiIhGjAk9EREQkYlTgiYiIiESMCjwRERGRiFGBJyIiIhIxvdOdgEiUeZ73DHB8Jy6Z7Pv+M12IMwu4Knz6rO/7J3T2Ht25l+d5vYBvAucARxD88bgaeB54FLC+72/rak6ZyPO8Y4EbgY8BDcBLwLW+77/UzvkHAN8HTgNGA5uB5cCzwIPtXSeZw/O8e4GvhU/v833/3PRlI9IxteCJSLd4npcLPAHcCRwH5AH7AWXAucBjwBue5+WlK8dk8zxvBPA08HEgFygETgGGtHN+CbAIuAw4DMgB9geOBL4HvOh53oPuMxeRfYVa8ETcep2gAIgpBErCx1uAha3Ob0hFUkl2BTA17vmbBP9vKQF6hcde9X1/Y6oTc+hyoF/4eAewBCgG/tnO+X8ARoSPtwBvAIPCa2KeTX6aIrKvUoEn4pDv+9+Kf+553oXAHeHTtb7vH5P6rJLH87ws4Py4Q9/1ff+28LWBBN1ZFxMURFFybNzj633fv9LzvCzf9/fYGsjzvPEELZsxR/m+vzh8rRT4NkGB/FuXCYvIvkUFnkiG8jzvEOASgl/+RQSte/8CZvu+/1TceW8DB8dderzneS3EjZ/zPC8b+CJBl+lHgQJgPfBieL+nu5jmQeFXTFbsge/764BbPc/7le/7ze3k+3Xf9++Ney2+QJrs+/4znuedACwIjz0CzARmAZ8OY68ArvZ9/0HP8yYQtCgeD/Qn+Lym+75f3dGb8DyvN7vGEI4D+oT3fRS4xff9D1pdEt/dvDx8v+3t+3hYq+fxn1EV8F3P83J839/aKqeDgavD9zmAoJXwemAN8O/wtPt83z/X87xzgXvCYyt93x8Zd59ZtDOm0vO8cQSf51RgGEHr4pvAXGBOq+9b7P1tJ/h8vg/MAJ72ff+88Jz9gIuAaQRd9NnAfwk+x9t8318fd78+wDcICtwygjGJC4E7fd+fu8enyM6u8bfjPsPxvu8vbee9/tD3/Vs8zzsQ+BFwFkEX+vsELa03+r4f+xy7xfO8coKf0d7h+yj3ff/VZNxbpKs0Bk8kA3me9yVgMfB1YDjQl6BL7xTgH57n3dSJe+UT/EJ7EPhUeJ8+BEXj54B5nufN7GKqHwLxhc3PPM/7k+d5X/I8b3+A+CIhCYYCLxO0Gg4jGMt2CPBHz/O+RzDR4QyC99gPmAz8qaMbep5XQPDL+Q6Clrn88NrDCMbM/ScstvE874Sw0BkZd4t7PM9rCSfUtKV11/Rznufd7nneyWGRs8dn5HneROBVghbQwjCfjwF/Bi7s6P0kyvO8HwH/AaYDowh+xvII/gC4GVjgeV5OO5f/LDxnBOHvEc/zhhJ8/jcSjC2MjcU8HLiS4HMcH3ePewg+8/Fh7ALgBOBBz/Pubiuo7/urCIr2mJNbnfLJuMcPxf3sXxLm2ofg8zwD+JfneV9v5/0lzPO8wUAFuxpMLlBxJ5lABZ5IhvE873DgXoJfehC0mLzB7oXCDz3Pi83mW0jQqhPTSNAy93rc8/7h483AWwStQfGzWm8Mu1Q7xff9D9l9HGE2wSzRB4DasNjrzCzivTmWoBXmXWBtq7i3EnxmVQSFZ8zhnucd0cE97wTK456/S9DqFDMceCxWjHXBK8CmuOcDgP8BniT4jG73PG907MWwtfWPwIGtcno7fPyNLubR2laCz20Hwc/PIoLWrZhjCVo1W+sFfDd8XAesDB8/QND6SXjPpQQ/t7E/AIoJPse+nucdB3wl7p4rgHfinrfXGgpB62LMzoIunKX8sfDpi77vryRoTYy1oMb+O4q1xmaz+38DnRa2/Fawa3LNbb7v39+de4okiwo8kczzI3YVd+8CR/i+fxhBoVEZd94VAL7vn8bu47de9X3/mNj4v7Dr8AKCrr583/cP8X3/cIJfxrFuwVx2HyfWGRfQ9uSQHIJi75mw6yxZ/kbQxVsK1LR6bbrv+2XAJ1odL23rRuEYuGlxh2YDg33fH0PQpR0rNMrC8xoIiuctcdcsZ/eCeje+79cRjENsq2g5kKDY+4/neSY89kl279a9Hijyfb8kPLej4qczbif43hX5vj/M9/1JBK268+LOmdLOtdUE3ZCDfd+/yvO8j7NrOaB3gUN83x8f/tyOY1dRNYqg9Wxo3L3+7vv+KN/3i4GxwK8Iun/b8xBBAQnwiXAWNwTdzLFJPRXhv/FxfuL7/mG+7x8Qvq+bklCM3UgwkxrgGeAH3byfSNKowBPJPPG/VG/1ff91gHAcWPwvkNHxLT978QZBi8wbnuc1e563juAXfFPcOQe2eeVe+L7/MkE33B20Pwv4Ks/zju7K/dvwLd/3N4ethy/EHV/g+/5vwpxeBerjXmtviZbJ7BrP1QBc4vv+jvAeDxO0ssV80vf9V8OJMfGth9fGF9Rt8X3/1wSthH+m7Vaj/YDfe57Xn92L0+XAlbHxfb7v39kqpy7zfb+JoAXuUc/zGj3P20TQGjso7rT2fia+4/u+H/f8pLjHBwHLwm7rFoLCd/+4148BfHb9rEzxPG+253mfAlb7vn9RG2Me4/N+h11/6PRjV4H1qfDfFoIiEILle3bm7HneFeHP4bO+7/+ovRgJmkxQuENQwJ4ZtbUepWdTgSeSeeInLbzV6rWqVs+L9nYzz/MOIuh+uxIYQ1DUbSAoOAbEnZq159WJ8X1/VVjgDAFOJ2hB2drqtM939f6tYsW32m2Oe7yq1anxr7X33uI/67fbGC8Y/3nv9bPuiO/7z4etrcUE495aL6kykOB7El9gLY0VnHHqupNHjOd5vyQoOMsJWnDrCL5/E+NOa+9z29zq+eBOhM7xfX8tQXfq7wl+TmYAfwfe8zzvQc/zDt3LPdrqpo2Nx3sh9jPi+/7fwuNPEXz/riEYw1fred7Pwm7drhoR93h/4CPduJdI0qnAE8k878Y9HtPqtda/+N5l7y5h19p7vwAG+r4/imAywZp2r0qQ53nHhQPs8X1/k+/7j/q+Pw04mt27MmOtQfFdjDvXCPQ8rx+pF//5jWhjnN2h7ZzbKZ7nnRLrSvR9v973/d/4vv8Jghmk8Q4E1sU9b6trua3VD9r8TEN7fK7h7NlY7CpghO/7Jb7vDySYPNFZ8WP31oYx47/6xz2+CMD3/WW+73+NoKA9iaBrtjdBV/i/Pc8r6yDewwRj6gA+Gb6fYeHzivgTfd9/yvf9kwmK+WnAX8LHFwN+OI6uq5bHPb7L87wB7Z4pkmIq8EQyzzNxj78X+0UXzkq9Me61mnCZjdYKwvNjv9jjJxgsj+tGOplgRmGMYU9DPM/LCgf+7yFcB+8O4C3P864Nl/aIWcXurXixYjK+++2k8D69Cca/xSvGvfjFhfcHfhp7r57nfZHdF3Ce35UAYfH7MPC653kXhrN2Y/7b6vQ17N7tfIjneZfG3evzwJlthIn/TIvC5WJis3Gnx702yPO8vuz+M7GesHj1PG8YwUzWmDLP81r/kdGWZ+IeDyFYeqU57AY+hKCI/AvBGMlt4c/Unz3PuxzI831/nu/7FxHsDgJBQfjp9oL5vl8fF3M8cF74OL57Fs/zijzP8z3POwfY5Pt+BfAFdnUPH8quiSGdtYxgxvGj4fPhdK04FnFC6+CJZJ6bCQb49yb4ZbnE87wqgokF/ePOuyHucWPc44me5/0XKAyX91gZ99psL1hsOZ+gxSO+S/KL4b/x4+jKCIqQbeHj1j7HrmLhcuByz/PqCFpxRsflu4NgDTvC+8W6AU/3PO9NgqK0dTffLwlmkzrj+/5bnuc9QjDwH4Ixjmd7nvchwYSAmJUEs0S74ocEE05KCIrhX3met4rgc44vLmJ792YRFMexLsCfhsvYfEBQLLWldaFY6XleNcHYyF5xx8cRTNSIX//taGC153nvset73BLmMYSgC7fDdQQJxrq9TFDwQPCHyEzP8zYQFFFZBBMeYi1e5wGnhl+zPM9bFh6PbzGNnwndlrnAieHjWGukH47Ri7mFYPLQccAd4X9Hgwh+/nfG8TxvCLsKtd/6vp/IotMv+L7/fvi9OYlgnOd0z/Pm+r6/YC/XijinFjyRDOP7/kKC5SlirV99CGZVxhd3t/u+f0fc8+db3WYUwS+cE4D/ZfcC8DCC4u4dglmUrT3H7l1+JUBpuN5Xa/3ZvUsRgrFOE9l9fN+lcQvS/q7V+WMJirt3CH4hx3R5TGAnnU8wCzamiN2Lu1rgVN/3W487S9R2dm/JzCZYR+8IdhVfm4Czfd/fFi54fB67F9+DCYq79bTRkuj7/pJW7yGP4HvQi2AW7ntxr2X5vv888HjcsSKC4q83QRdq6/GMHQongXyB3ccsFhP8rMW+j/PYNUmokl1FaS+Cwi6+uFtBXEtcO/7ErgkrsfX6Klqd8wi7CsX9gAns3jL8p3AR7ByCQvdodnX1JsT3/dXsWlw5C/htuOCzSFqpwBPJQL7v30ewWOz9BN12WwkKqSeAz/q+/+1W578AfIugpamJoMXlPoJB+ssBj2Cg+WaCX/b3EAwK/wPBL8XmME5sVux0gl/AzQS/7P/QTp4PELQ0TQvjvUFQTG4jKNgeBk7wff+muGv+TrAMxmqCMXorCLpnJxC0/Dwfvt/4marO+L7/PsHM1ZkERdJGgs/wLeAm4HA/3Fqsi/f/AUH33UUEkxpWEnwfthB8n+4IYzwTd83TBC1n/wjz2UDQxXkcey4NE/Pl8PzYJJoFwOm+719OsKTO+vC12G4S0whai+vDfHzgU77v3wVcSvBzsinu/L29z5UEP1NXECzSvYmg1fElgj8kjB/uR+z7/lsEXavfJSj2PgjPf5PgZ+Bo3/c37CXeOnZf0mUHwc9b/Dl/JhjHeiPwWlxOrxC0+p2VyHtLwG0EE5kg+OPgp0m6r0iXZbW0JGtJJRERcc3zvHsJdrgAqAgntIiI7EYteCIiPUv8X+Vf8DzvtLRlIiIZSwWeiEjPEj9pJptw1rSISDwVeCIiPcutwGME4+k+ZPcJNCIiQA8cg2eMOZdg1lsL8Gtr7f1xrx1NMGMsj2A9pVnW2tar6YuIiIhEWkaug2eMySKY0XY6wcyopcCPCabhH8+u/SN/bYzZYa39ozGmjGCV/tMIttw5n2CvzfNT/w5ERERE0idTu2hHE0xnP8laO5lgSvvVBMtAzLTWbrHWNhMUgTPCay4CrrTW1lprW6y1dwHjjDH7t3H/na644ooWgtZAfelLX/rSl770pa9M/0pIRrbgWWur2X3l9ByCBVVzrbWNcec1GGNiC0qOJliDK96bBGsgvRx/0BgznXD7nkGDBrFuXet1WkVEREQyz8CBAxM6LyMLvBhjzEnAbwhmjZ0JPGGM6Wut3RK+nsOulczXECy4ujruFiNpY1HQsHXvLoCZM2e2JPphiYiIiPQEmdpFC4C19imCbZIqCPa8nEuw32XM5ezaq/Ie4CpjTC7sLA4brbV1qctYREREJP0yusADsNa2AHcSjL+7BehljHnJGBPbLPsX4XkvEBR5/zTG/ItgL8fz0pCyiIiISFplZBdt2Po2HrjNWrsD+BKw1Fq7Dbgs/NqDtfZB4MGUJSoiIiKSgTKywAPmE2y0/rQxBmAJwSxZEREREdmLjCzwrLXbgRvDLxERERHphIwfgyciIiIinaMCT0RERCRiVOCJiIiIRIwKPBEREZGIUYEnIiIiEjEq8EREREQiRgWeiIiISMSowBMRERGJGBV4IiIiIhGjAk9EREQkYlTgOdbS0sK76zewbdv2dKciIiIi+wgVeI7d9MI/OPf3v6V65TvpTkVERET2ESrwHGvZrxf5RQNYtnxNulMRERGRfYQKPMemDi3j9SffZNkKFXgiIiKSGirwHPto4Qiy67eowBMREZGUUYHn2JINtRx+/FiqVOCJiIhIiqjAc2xBfTWDPjaU6pXvsGPHjnSnIyIiIvsAFXiOFebkcUCvXDY3baFm7XvpTkdERET2ASrwHJs2YiInF4wG0ExaERERSQkVeClQVjIMgGUrVqc5ExEREdkX9E53AlE3p6oSgIMOLNBMWhEREUkJteClSGlJsbpoRUREJCXUgufYmcMnALCy5E0effJ5WlpayMrKSnNWIiIiEmVqwXOsKHcARbkDKCspZkPjh9Sv25DulERERCTi1ILn2Py6agDKRhUDwUSLokH7pzMlERERiTi14Dm2tKGWpQ21lJUEBZ52tBARERHX1ILn2Lj8wQAMKTyQvP79NNFCREREnFOB59iUojE7H5eOHKqlUkRERMQ5ddE6VtfUSF1TIwBjRw1TgSciIiLOqcBzrKJmMRU1iwEoKymm9t33aWjclOasREREJMpU4KVQaWyixdtqxRMRERF3NAbPsRml5Tsfx2bSvrV8NR85vDRdKYmIiEjEqQUvhQ4uLqRvn94ahyciIiJOqcBzbO6qRcxdtQiA3r17MfrgISrwRERExCl10TpW37xxt+dlJcX8580VacpGRERE9gVqwXNscuEYJhfGrYVXUszKNfU0NW9JY1YiIiISZSrwHBtfMJjxBYN3Pi8rKWbHjhb+u3JtGrMSERGRKFOB59iSDbUs2VC783lsJq3G4YmIiIgrGoPn2IL6aoCdrXijDx5KdnaWCjwRERFxRgWeY4U5ebs975fbl4OLC1XgiYiIiDMq8BybNmLiHsdKS4qpUoEnIiIijmgMXhqUlRTz35Vr2bZte7pTERERkQhSgefYnKpK5lRV7nasrKSY5i1bWbmmPk1ZiYiISJSpwEsDzaQVERERl1TgOXbm8AmcOXzCbsfKSoYBaByeiIiIOKFJFo4V5Q7Y41j+gP0oGnSAWvBERETECbXgOTa/rpr5ddV7HC8bVawCT0RERJxQgefY0oZaljbU7nG8rCQo8FpaWtKQlYiIiESZumgdG5c/uM3jY0cVs/HDzdS++z5DCg9McVYiIiISZSrwHJtSNKbN46WxmbTLV6vAExERkaRSF61jdU2N1DU17nE8NpNW4/BEREQk2VTgOVZRs5iKmsV7HC8cWEDBgP4q8ERERCTpVOClSVZWFmUlxby1XAWeiIiIJJfG4Dk2o7S83dfKSor5xz9fTWE2IiIisi9QC14alY4q5t31G3h/w8Z0pyIiIiIRogLPsbmrFjF31aI2X9u1J+3qVKYkIiIiEacCz7H65o3UN7fdQrezwNM4PBEREUkijcFzbHJh2+vgAQwfchC5OX00k1ZERESSSgWeY+ML2t7JAqBXr2xKRxZTpQJPREREkkhdtI4t2VDLkg177kUbUxruSSsiIiKSLCrwHFtQX82C+up2Xy8rKaZm7Xts2tycwqxEREQkylTgOVaYk0dhTl67r5eVFNPS0kL1yndSmJWIiIhEWcaOwTPGZAHnA2cA/YCngOustS0dXHM0cD2QBzwNzLLWbk1Buu2aNmJih6+XjQpm0r61fDVHHFKSipREREQk4jK2wAPGAoXAKUALcB9wmjFmIHAeEOvTfM1aO8MYUwb8AjgNqCMoDm8P/81Yo0YMoVevbE20EBERkaTJ2ALPWvsmcF3suTFmDbAN+BRwqrW2vtUlFwFXWmtjMxruMsaca4zZ31r7QUqSbsOcqkqg/S3L+vbpzchhRVoLT0RERJImYwu8GGPMCOBSoBF4HLgGuNkYczCwHLjUWlsHjAbeaHX5m8AY4OVW95wOTAcYNGgQ69atc5b/tq1BD3FHMUYWH8Qb1auc5iEiIiI938CBAxM6L6MLPGPMBcAJBGPp3jLGHAasB35kra01xnwGeAj4BLAGGAHE7/s1EqhpfV9r7V3AXQAzZ85sSfTD6oov9/8oAANzB7R7zvixJTz70lLy8wvo0yejvyUiIiLSA2RsNWGMmQh83Fr7pdgxa+3rwJS4548bY75tjBkO3APMMsacaq1tMsacBDSGrXtpU9RBYRdTVlLMtm3beXt1HaXh9mUiIiIiXZWxBR7waeBjxphn4o49CvzGWrsp7lhfYKu19gVjzD3AP40x24EVBJMx0mp+XbAG3pSi9rcsKxs1DIBlK9aowBMREZFuy9gCz1p7A3BD/DFjTDnwqDFmmrX2A2PMF4EtsYkV1toHgQdTn237ljYEcz46KvBiRd2yFWv4TEqyEhERkSjrUQsdW2srgV8CDxljFgDHAmelN6uOjcsfzLj89vejBcjbL5fiwQO1ZZmIiIgkRca24LXHWvtX4K/pziNRHbXcxSsrKdZSKSIiIpIUPaoFryeqa2qkrqlxr+eVlRRTtWINO3bsSEFWIiIiEmUq8ByrqFlMRc3ivZ5XWlLMpqZm1tRqLTwRERHpHhV4GaIsbqKFiIiISHf0uDF4PU17W5S1Vlaya6mUE72JLlMSERGRiFMLXoYYdGA+B+4/QC14IiIi0m0q8Bybu2oRc1ctSujc2EQLERERke5QgedYffNG6ps3JnRu2ahiteCJiIhIt2kMnmOTCxNbBw+CFrz1HzTy3voGBh2Y7zArERERiTIVeI6NL+h4F4t4u2bSrmbQgYe5SklEREQiTl20ji3ZUMuSDbUJnRs/k1ZERESkq9SC59iC+mogsZa84sED6d8vRxMtREREpFtU4DlWmJOX8LlZWVmMGVnMW9qTVkRERLpBBZ5j00Z0btHislHF+C+/7igbERER2RdoDF6GKSsp5p26dTR+uDndqYiIiEgPpQLPsTlVlcypqkz4/NhM2qq31U0rIiIiXaMCL8PsLPA0Dk9ERES6SGPwHDtz+IROnT9yWBF9evfSUikiIiLSZSrwHCvKHdCp8/v06c2oEUNU4ImIiEiXqYvWsfl11cyvq+7UNWUl2pNWREREuk4FnmNLG2pZ2pDYThYxZSXFvL26ji1btznKSkRERKJMBZ5j4/IHMy4/8f1oAUpHFbN9+w6Wr1rrKCsRERGJMo3Bc2xK0ZhOXxObSfvW8tUcMnp4slMSERGRiFMLnmN1TY3UNTV26poxBw8lKyuLZVoqRURERLpABZ5jFTWLqahZ3Klr9uuXw/AhgzTRQkRERLpEBV6GKisppkoFnoiIiHSBxuA5NqO0vEvXlZYUU/nyUrZv30GvXqrDRUREJHGqHDJUWUkxTc1bqVn7brpTERERkR5GBZ5jc1ctYu6qRZ2+buyoYQAahyciIiKdpgLPsfrmjdQ3b+z0dWWjgqVSVOCJiIhIZ2kMnmOTCzu/Dh7A/vl5FA4soEpLpYiIiEgnqcBzbHxB53axiFeqPWlFRESkC9RF69iSDbUs2dC5vWhjykqKeWv5alpaWpKclYiIiESZCjzHFtRXs6C+ukvXlpUU07BxE3XvfZDkrERERCTKVOA5VpiTR2FOXpeuje1Jq25aERER6QyNwXNs2oiJXb62LFwqpWrFGj5x1PhkpSQiIiIRpxa8DDb4oAPI699PLXgiIiLSKSrwHJtTVcmcqsouXZuVlUWZZtKKiIhIJ6nAy3BjRxWzTGvhiYiISCeowHPszOETOHP4hC5fX1pSTN1779PQuCmJWYmIiEiUqcBzrCh3AEW5A7p8fVlJbE/a1clKSURERCJOBZ5j8+uqmV/XtXXwQEuliIiISOepwHNsaUMtSxu6tpMFwMHFheT07cNbGocnIiIiCdI6eI6Ny+/6XrQAvXplM/rgIWrBExERkYSpwHNsStGYbt+jtKSYxa8vT0I2IiIisi9QF61jdU2N1DU1duseZSXFrFxTz+amLUnKSkRERKJMBZ5jFTWLqahZ3K17lJUU09LSwn9XvpOkrERERCTKVOD1AJpJKyIiIp2hMXiOzSgt7/Y9Rh88hOzsLBV4IiIikhC14PUAuTl9Obi4iCoVeCIiIpIAFXiOzV21iLmrFnX7PmUlxWrBExERkYSowHOsvnkj9c0bu32fslHF/HflWrZt256ErERERCTKNAbPscmF3V8HD4IWvC1bt7FyTT2jDx6SlHuKiIhINKkFz7HxBYMZX9C93Sxg10zat5av7va9REREJNpU4Dm2ZEMtSzZ0fS/amNKRQYGniRYiIiKyN+qidWxBfTVAt1vx8gfsx+CDDtBECxEREdkrFXiOFebkJe1emkkrIiIiiVCB59i0EROTdq+ykmIe/NtztLS0kJWVlbT7ioiISLRoDF4PUjZqGBs/3Mza+vXpTkVEREQymAo8x+ZUVTKnqjIp9yobpT1pRUREZO9U4PUgsaVSli1XgSciIiLt0xg8x84cPiFp9zrowAL2z++vFjwRERHpkAo8x4pyByTtXllZWZpJKyIiInuVkQWeMSYLOB84A+gHPAVcB/QBrgZOAD4ELrXW/jvuunPD61qAX1tr709p4m2YXxesgzelKDlblpWWFPPks68k5V4iIiISTZk6Bm8sUAicAkwBDgFOA+4ClgPHAWcDs40xo2FncXc8MBk4ETjRGPOVlGfeytKGWpY2dH8ni5iykmLee7+B9R80Ju2eIiIiEi0JteAZY0YCU6y1d4fPvwr8nKDYOstauyKZSVlr3yRosYvFXwMMAMZYa88ND681xlwDfBu4GPgWcKK1dkt4zUzgH8Afk5lbZ43L7/4+tPHKSoYBwUzaYyYdktR7i4iISDQk2kV7E3A/gDGmD0Ex9RGCLtRfAKe6SM4YMwK4FGgEXgfebHXKUoICDyDXWruzWcta22CM2a+d+04HpgMMGjSIdevWJTv1nSb0PgAgaTEOOiB4Swtfe4vSEQcl5Z4iIiLSMwwcODCh8xIt8E4Avhw+Phb4i7V2lTFmNvD9TmeXAGPMBWHcWdbat4wxQ4ERrU4bCdSEj7cbY/rGteDlANvaure19i6C7l70qxFxAAAgAElEQVRmzpzZkuiH1RV1TUHNmazJFgcccAD9cvuypn5Dwt9kERER2bckOgZvi7U2VixNBp4EsNZu78Q9EmaMmQh83Fr7JWvtW2Gsd4Btxpip4Tm5wGXA3eFlc4HL425zOWnungWoqFlMRc3ipN0vOzubMSOHUqWZtCIiItKORFvwFhpjLgQeA06y1l4NYIwZSzAOL9k+DXzMGPNM3LFHgXOAXxljrgN6ATfHzaK9BbjaGPMSkAU8QdB9HDllJcW8tGhZutMQERGRDJVogfddgtawnwLfizt+JXBbspOy1t4A3NDOy9PauWYbQYveZcnOpztmlJYn/Z5lJcU88nefDzc30b9fbtLvLyIiIj1bQgWetfa/wDFtvHSFtdZFC550IDaTtvrtd5hw6Kg0ZyMiIiKZJuGFjo0xA4CfAGOttaeHkxjedpVYVMxdtQiAaSMmJu2eO/ekXbFGBZ6IiIjsIaEJEsaYQuB5ghmrx4WHjwZeNMYUOMotEuqbN1LfvDGp9ywZMZhevbJZtlwTLURERGRPic6AvQm41lp7O8E2YFhrnwP+ANzoKLdImFw4hsmFydmmLKZvn96UDB/MshWrk3pfERERiYZEu2g/CXw9fNwSd/x2oDqpGUXM+ILk7mQRU1ZSzDItlSIiIiJtSLQFr5e1tqWN41kEW4hJO5ZsqGXJhuTtRRtTVlLMipo6tm5tcy1nERER2YclWuC9ZIw5NnycFXf8G8CLyU0pWhbUV7OgPvmNnGUlxWzbtp0VNckvHkVERKRn68w6eBXGmIsAjDHjgLOB84GjHOUWCYU5eU7uWzYqWCpl2Yo1Ox+LiIiIQIIteNbaaoIFhn8M5AEvAGXAUeEaedKOaSMmJnWJlJjSkqEAGocnIiIie0ioBc8YU2atXQZ8znE+kqD+/XIZNniQCjwRERHZQ6Jj8J4xxiS8KLLsMqeqkjlVlU7uXTZKM2lFRERkT4kWeHcDF7pMRDqvtKSY6hXvsGPHjnSnIiIiIhkk0Va5J4HLjTFDgWdav2it/Ucyk4qSM4dPcHbvspJiNjU1s7r2PUYMLXQWR0RERHqWRAu8a8J/jw2/4rUAKvDaUZTrbpnAnXvSLl+jAk9ERER2SqjAs9ZOdp1IVM2vC9bAm1KU3O3KIK7AW7GGqeWTkn5/ERER6ZkSHYPXJmNMtjHm9GQlE0VLG2pZ2uBmMeKBB+QzcP8BmmghIiIiu0l0mZT9gP8FTgeGxL3UBDwL/Cn5qUXDuHw3e9HGlJYUU6UCT0REROIk2oJ3AzAQOAa4E5gMDCKYfHFNB9ft86YUjXHSPRszdtQwlq1YQ0tLW1sFi4iIyL4o0QLvVOB8a+1qguLueWvt+8DFwFWukouCuqZG6poand2/bFQx72/YyHvvNziLISIiIj1LogVeP2vtpvBxvrV2G4C19m3gMBeJRUVFzWIqahY7u3/8TFoRERERSLzAW2SMGRs+ftsYMwnAGFME9HGSmSQkfiatiIiICCS+Dt5PgbOAq4HZwAPGmLvDY79zlFskzCgtd3r/oUUD6b9friZaiIiIyE6JroP3LMFsWay1FcaYZmAKcKu19o8O85O9yMrKonTkULXgiYiIyE6JtuDtxlr7GPAYgDEmy1qrKZztmLtqEQDTRkx0FqOspJh//nups/uLiIhIz5LoOnj/R7AlWVvXnwBon6x21DdvdB6jrKSYisf/SePGTQzI2895PBEREclsibbgPdHGsc8COwDtkdWByYXu1sCLKRs1DICqt9/hyPHu44mIiEhmS3QM3n1tHL7PGPMz4HjggaRmFSHjC9zuZAHBbhYQLJWiAk9ERES6tRctMAv4YRLyiKwlG2pZssHNXrQxJcOK6NunN69Xr3IaR0RERHqGLk2yiPMhMCwZiUTVgvpqwG1LXu/evfjslYY127Y7iyEiIiI9R6KTLPq2OtQLGA5cAryR7KSipDAnLzVxcvN4dUk179StY2jRwJTEFBERkcyUaBdtE7A57quBYF28fOBsN6lFw7QRE50ukRLzlZIjefG+l5jnL3IeS0RERDJbopMsujtWTxwbO2oYwwYPYl7lQs45/cR0pyMiIiJplNTCzRgzOZn3i4I5VZXMqap0HueX1T4nfP8Enn3xNZq3bHUeT0RERDJXomPw1tL2Qset7/UBUNbdpKRrDijI48NNTby48E0+cfTh6U5HRERE0iTRWbRHEKx1dxvwWnjscOCbBBMtmsNjbtcD6YHOHD4hZXE2H7SFX/f5A09VLlSBJyIisg9LtMC7AbjJWvtU3LGVxphtwPnWWq2F146i3AGpi5ML3kcPY56/kGu/f05K4oqIiEjmSXQM3kmtijsArLVPAKcnN6VomV9Xzfy66pTFmVo+iaoV7/D26jrnMUVERCQzJVrg9THG5Lc+aIzJA/ZPbkrRsrShlqUN7nuuY3GmesHWwE9ruRQREZF9VqIF3gPAHcaYrFbHfw78LbkpRcu4/MGMy3e/H20szuiDh1AyfDBPVS50HlNEREQyU6Jj8H4CPAYsNMY8AuwATg3/neoot0iYUjQm5XFOKp/E7/80j81NW+iX23oTEhEREYm6hFrwrLVbrLWfBq4CDgIGAzcDx1prNzrMr8era2qkrqkxpXGmlk+kqXkr/stLnccVERGRzNOphY6ttY8B3wHuBWo6e/2+qKJmMRU1i1Ma57iPHEa/3L7MUzetiIjIPqndAs0YU2iMuaXVsd4EY+4scBfwsjGmyG2K0lm5OX35+FHjmecvoqVlb+tTi4iISNR01AJ3OVDf6thXgRHAaGvtEQQteVe7SS0aZpSWM6O0POVxpnqTeHt1HdUr1zqPLSIiIpmlowLvFIJWunhfBm6OG3d3B3Cyi8Ske6aWTwRQN62IiMg+qKMCL9da+0HsiTEmF/CAx2PHrLVbgFx36fV8c1ctYu4q92vStY4zYmghY0cNU4EnIiKyD+qowNtkjOkV9/w44E1r7brYAWNMH2eZRUR980bqm91PNG4rztTySTz/yuts3NTkPL6IiIhkjo4KvPlA/IamFwGPtjrnZODfyU4qSiYXjmFyofu18NqKM7V8Ilu3bee5F19zHl9EREQyR0cLHV8F/N0Y8xUgn2BLsvNiL4atd7OAK10m2NONL3C/i0V7cY6eeAh5/fsxz1+ImfyxlOQhIiIi6dduC561di1wNHArwUzZj8SPyQOKgcettX93m2LPtmRDLUs2uN+Ltq04ffv05oSjD2depZZLERER2Zd0uFWZtbaZYM27tl57m6AFTzqwoL4acN+S116cqeWT+Nv8l3ijuobDSkc4zUFEREQyg3aicKwwJ4/CnLy0xTnR03IpIiIi+5oOW/Ck+6aNmJjWOEMKD+TwsSN5qnIhM79+akpyERERkfTqaKuyUcaYrFQmI25MLZ/ES4vfYkPjh+lORURERFKgoy7aRwgXMTbGXJeadKJnTlUlc6oq0xpnavlEtm/fwTP/+o/zPERERCT9OirwhgDbwsfndXCeZLiPjC9l//z+GocnIiKyj+hoDJ4FnjLG/AnoZ4yZ3u6J1rbes1ZCZw6fkPY4vXv3YvKxE3jaX8SOHTvIztbcGhERkSjrqMCbDnwF+BiQAxzbznktgAq8dhTlDsiIOCeVT+LRJ5/nP2++zcTDRqUkJxEREUmPdgs8a+024D7gPmNMvrX266lLKzrm1wXr000pcrtd2d7iTD5uAllZWcyrXKgCT0REJOISWibFWnsOgDEmGzicYOzea2ERKB1Y2hDsLuG6wNtbnIMOLGDSuNHMq1zID6af4TQXERERSa+EB2OFe9K+AzwGPAqsNcac6yivyBiXP5hx+e73o00kztTyibyypJp17zc4z0dERETSJ6ECzxhzAsG2ZKdYa0daa0cCBrjMGHOSs+wiYErRGOetd4nGmVo+iZaWFua/sNh5PiIiIpI+ibbg3Qx81Vr7cuyAtfbfwNnATS4Si4q6pkbqmhozIs7EQ0cx6IB8LZciIiIScYkWeMXW2hdbH7TW/gvQDvYdqKhZTEWN+xazROJkZ2dzojeR+c8vZvv2Hc5zEhERkfRItMDbZIwpan3QGFMMbExuSuLS1PJJvL9hI68uqU53KiIiIuJIQrNogdnA/caY06y1HwIYY3KBO4E5rpIL4wwBzgKetNa+vpdzjwauB/KAp4FZ1tqtLvPbmxml5RkVZ/KxR5CdHSyX8rEJZY6zEhERkXRIdJmU2caYMuBNY8xTwHbgJOAJa+0trpIzxkwDPhk+rQFeN8acT7B1WnN4/DVr7Ywwv18ApwF1wPnA7eG/Eto/P4+jJozlKX8hl357WrrTEREREQcSbcHDWnuRMeY3wInhoTusta+6SWtnzLnAXGPMrLjDnwJOtdbWtzr9IuBKa21t+PwuY8y5xpj9rbUfxJ8Ybrs2HWDQoEGsW7fOzRsA/rauCoBTBpY6i9HZOMdNKuPWu//CG8uWUziwwGleIiIikjwDBw5M6LyECzwAa+1iIG1rbBhjegGjgZuNMQcDy4FLrbV14fE3Wl3yJjAGeDn+YLh37l0AM2fObEn0w+qKhvVBSi5jdDbO5072uPXuv7DwjZV8+dTJTvMSERGR1Otpu86PBdYDP7LWngA8AjwUvraGPWf0jiTo2k2byYVjmFzofh28zsQZV3Ywgw86QMuliIiIRFSPKvCsta9ba6fEumGttY8DG40xw4F7gKvCyR+ECzA3hq17aTO+YDDjC9zvZNGZOFlZWUwtn8SCf/2HrVu125yIiEjUZHSBZ4wZY4z5G/Bl4MfGmLuNMfu1Oq0vsNVa+wJBkfdPY8y/CCZinJfajPe0ZEMtSzbU7v3EFMc5qXwSjRs389LiZQ6zEhERkXTo1Bi8VLPWVgOnxJ4bY8qBR40x06y1HxhjvghsiWvRexB4MD3Ztm1BfbDenOtWvM7G+cRR4+nTuxfz/IV4Hz3MZWoiIiKSYgkVeOGEhluBcmAgkBW+lAW0WGt7uUlvd9baSmPML4GHjDG9gYUEa+RlrMKcvIyMMyBvP4458lDmVS7kqu98xVFWIiIikg6JtuA9DLwKfJxgMkPaBm5Za/8K/DVd8Ttr2oiJGRtnqjeRq37+B1avfY9hQwY5yEpERETSIdECrwQ4ylrb4jIZSa2Tyidx1c//wDx/Ied+4aR0pyMiIiJJkugki9cJumalk+ZUVTKnqjIj45SWFDNi6EHMq1zkKCsRERFJh0QLvIuBObElSCQaYsulPPfiazRvSeuWvSIiIpJEiXbRvgqsAt42xiwAduuqtdZ+OdmJRcWZwydkdJyp5ZO4u+IfvPDqG5xwzBFJzkpERETSIdEC727gGOCXwFrSOMmipynKHZDRcco/Oo6cvn14qnKhCjwREZGISLTA+ywwxlr7vstkomh+XbA+3ZQit9uVdTXOfv1yKP/oOOZVLuT6H3zNRWoiIiKSYomOwasFGlwmElVLG2pZ2uB+J4vuxJlaPpH/rlzL8lXu8xQRERH3Ei3wZgOXukwkqsblD2Zcvvu9aLsT50RvEgBP+wuTmZKIiIikSaJdtO8Bpxlj/gI8zp6TLO5KdmJR4bprNhlxRo0YzOiDhzCvchHnf+nTScxKRERE0iHRAu8UgskVEEy2iNcCqMBrR11TI+B+skV340wtn8S9Dz3Fps3N7NcvJ5mpiYiISIolVOBZa7/uOpGoqqhZDMCM0vKMjnNS+SR+/UdL5ctLOfnjRyYzNREREUmxhAo8Y0zfDq4/1Vr7f8lLSdLh2CMPZb/cHOZVLlSBJyIi0sMl2kXbRNAVmxV3rAVYA7wCqMBrh+uWu2TFyenbh08cfTjzKhfS0tJCVlbW3i8SERGRjJRoF+0es22NMZ8FLgC+muykJD2mlk/kiWdfpmrFGspGDUt3OiIiItJFiS6Tsgdr7V+BO4D/TV460TN31SLmrlrUI+JMDZdLeapSy6WIiIj0ZF0u8ACstY8Dn0lSLpFU37yR+uaNPSLOsCGDOHTMcOb57gtSERERcSfRMXhtMsYMB/KSlEskTS5MzTp4yYoztXwSd/7hcRo3bmJA3n5JuaeIiIikVqKzaKe3OtQLGE4w/u73yU4qSsYXuN/FIplxpnqTmHPvX3jupSV8ZspRSbmniIiIpFaiLXjHtnq+nWB3iwuBvyc1o4hZsiHY39V1oZesOEdNKGNAXj+eqlyoAk9ERKSH0kLHji2orwbcF3jJitOnT28mH3OElksRERHpwTqcZGGM2WPFW2PMWcaY14wxbxpjfugutWgozMmjMMf9MMVkxplaPonad99n6bKVSbmfiIiIpFa7BZ4x5pvAta2OfYRg39nrgK8AU4wxFzrNsIebNmIi00ZM7FFxphwX3GeelksRERHpkTpqwZsO/KjVsYuAOdbaudbaVwgWOr7AVXKSHoMPOoAjDi3RcikiIiI9VEcF3ghr7ZJWx04mblsya+0qYJCLxKJiTlUlc6oqe1yck7xJvLT4LT5ocL+Gn4iIiCRXRwVec/wTY8yhQHYbRV/fpGclaTe1fBI7drSw4IX/pDsVERER6aSOCrzlxpjxcc9PB56MP8EYcxiwykViUXHm8AmcOXxCj4tz5PgxHFCQp3F4IiIiPVBHy6TcADxgjLkJ2B+4BJjc6pxvAxWOcouEotwBPTJOr17ZTDluAk/7i9ixYwfZ2d3a1U5ERERSqN3f2tbafwAzgM8BxwNftta+GnvdGFMGnAjc7jrJnmx+XTXz66p7ZJyp5ZN47/0GFr2xPKn3FREREbc6XOjYWvss8Gw7ry0zxhxurd3qJLOIWNoQ7DAxpcjtnrQu4kw5dgJZWVnMq1zIkeNSs6euiIiIdF+3+t1U3O3duPzBjMt3vx+tizgDD8jnI+PHMK9Sy6WIiIj0JInuRStd5LrlznWcqeWTuPHOh3h3/QYOOrDASQwRERFJLo2cd6yuqZG6psYeG2dq+SRaWlpY8PzipN9bRERE3OhWgWeMyQq3NJN2VNQspqLGfXHkKs4Rh4ykcGABT2m5FBERkR6ju2PwWmi1X61ES3Z2Nid6E1nwwmK2bdue7nREREQkAe2OwTPG/B/QspfrS4D1Sc0oYmaUlvf4OFPLJ/F/f3mWV5ZUcfTEQ5zFERERkeToaJLFEwlcvwWYn6RcJEOdcMwR9OqVzbzKRSrwREREeoB2Czxr7X2pTCSq5q4KlhiZNmJij41TMKA/R00Yy7zKhVx20VlJv7+IiIgkV8LLpBhj+gPHAEOBrLjrP2OtPcNBbpFQ37wxEnFOKp/ENbMfYG39eoYUHug0loiIiHRPQpMsjDETgDeAS4HvAd8ETgduAv7sLLsImFw4hsmF7tfCcx1navkkAJ72teixiIhIpkt0Fu2twMXW2qnAGuBsa+3ngWuAga6Si4LxBYMZX+B+JwvXcQ4dM5yhRQOZp+VSREREMl6iBd6h1tqHw8fjrbUrw8e/Ac5NelYRsmRDLUs21Pb4OFlZWXzyE0cyz1/Iihr370dERES6LtECb5sxJjZer8UYkwtgrd0MuG+e6sEW1FezoL46EnEGnzqGEy4+ge9f9xtaWva2go6IiIikS6IF3mPAp8PHLwBfAzDGfAaoc5BXZBTm5FGYkxeJOMMG7E9Z4WCee2kJc//6nNNYIiIi0nWJzqK9BpgcPr4KeNoYMxMYBHzRRWJR4Xp5lFTGmTZiIjuGHcGLE1/iilt/z4nlEznowALncUVERKRzEmrBs9a+a62tCB9XA2OBrwKjrbVqytmHZGdnc+vl09n44WYuv0VLJYqIiGSiRJdJ+Vb8c2ttk7V2IZBjjHnKSWYRMaeqkjlVlZGIE4sxdvQwvnPe53nk775m1YqIiGSgRMfgXd7O8feB1PRBSkb53jdOo7RkKD/86e/YuKkp3emIiIhInHbH4Blj/g+ITZUsMMY80MZpg4GVbRyX0JnDJ0QmTnyMnL59+PkVF3DKeVdx4x0VXPv9c5zHFxERkcR0NMniibjHBniyjXO2AOqi7UBR7oDIxGkd45hJh/C1M6by6wcsp3/KY9K40c5zEBERkb3LSmQ9M2PMhdbaO1OQT8rNnDmzZfbs2c7uP78uWJtuSpHb7cpSEaetGA2Nmzj29IspHFjAP+6/nj59Et7eWERERDovK5GTEp1FeyeAMeZwY8x3jDHfM8Z8pDvZ7SuWNtSytMH9zg+piNNWjPwB+3Hjj7/Oa2+9zR1/fNxpfBEREUlMopMsMMbcRdBNOwmYAPzVGHOPq8SiYlz+YMblu9/sIxVx2otxyolHYyZ/jJvufEjbmImIiGSARLtovw2cCpxhrW0Mj+UDDwF/t9b+wmmWDrnuot1XrK1fz3FnXMykw0bzyJ2Xk5WVUAuyiIiIdE7yumiB7wBfixV3ANbaBuCb4WvSjrqmRuqaGvd+Yg+I01GMIYUHcsWML2sbMxERkQyQaIGXZ61d2/qgtbYGyE9uStFSUbOYiprFkYiztxjnfmEqR08cyxW3/p53129wmouIiIi0L9ECb60x5qOtD4bHVic3Jemp4rcxu+Jnv093OiIiIvusRNe0+BHwoDHmTGvtqxDMqAXuBy52lVwUzCgtj0ycRGLEtjG75a5HOPMzH2fKcdroREREJNUSmmQBYIz5InAb0ARsBwYAF1tr29rhosfQJIvka96ylRPOuoTm5q388+Fb6N8vN90piYiIREVSJ1lgrX0IOBj4AjANGNHTi7tUmLtqEXNXLYpEnERjxLYxW/XOu/zv7RVOcxIREZE9dbQX7Xxr7ZT4Y9barcCrzrOKkPrmjZGJ05kY2sZMREQkfTpqwTs0ZVlE2OTCMUwudLtNWaridDbGVd/5CgcduD8XX3sXW7duc5iZiIiIxOuowEtscJ50aHzBYMYXuN/JIhVxOhtD25iJiIikR4ezaI0xfUhgMJ+1dkvSMto9/hDgLOBJa+3rxpi+wNXACcCHwKXW2n/HnX8ucD5Bcfpra+39LvLqjCUbgq27XBdfqYjTlRjx25h99sSjKRnuvtgVERHZ13XUgjeYYMbs5g6+Yq8nnTFmGnA9cDhwWHj4LmA5cBxwNjDbGDM6PP9c4HhgMnAicKIx5isucuuMBfXVLKivjkScrsa48cfn0adPb35w/W9JdNa2iIiIdF1HLXh11tohKcukFWvtXGCuMWYWgDFmIDDGWntueMpaY8w1wLcJ1uL7FnBirDXRGDMT+AfwxxSnvpvCnLzIxOlqjNg2Zpfc8Dsq/vYc0z57fJIzExERkXgdFXiZ1tQyCniz1bGlBAUeQG7rvXKNMfu1dSNjzHRgOsCgQYNYt26dg3QDU/sPB3AaI1VxuhPjs5Mn8eBfRnPZLfdx5GEHc+D+A5KdnoiISOQNHDgwofMS3ckiE6wBRrQ6NhKoCR9vN8b0jWvBywHanLpprb2LoLuXmTNntiT6YUn3zLn6W5xw1iX87O6/cuf1M9KdjoiISGR1NAbP7c71nWStfQfYZoyZCmCMyQUuA+4OT5kLXB53yeWkuXsWYE5VJXOqKiMRp7sxYtuYPWwrmf+8+8WfRURE9lXttuBZaz+dykRaM8aMAX4BlAENxpiTgXOAXxljrgN6ATfHzaK9BbjaGPMSwczfJ8LrJYN87xun8dhTL/CD63+rbcxEREQcSXgv2qhyvRdtXVMwLLAo1+2Ys1TESVaMfy18k1POu4r/+epnuPb75yQjNRERkX1Fcveila4pyh3gvLhLVZxkxYjfxmzh0v8mITMRERGJpwLPsfl11cyvc78OXiriJDNG/DZm27ZtT8o9RUREJKACz7GlDbUsbaiNRJxkxtA2ZiIiIu70pGVSeqRx+anZmisVcZIdo/U2ZiOHFSX1/iIiIvsqTbJwPMlCOra2fj3HnXExR44bw8N3XEZWVkJjR0VERPZVmmSRCeqaGnfOPu3pcVzEiG1j9uyLr1Hxt+eSem8REZF9lQo8xypqFlNR437N6FTEcRXj3C9M5eiJY7ni1vt5b31D0u8vIiKyr1GBJ2mXnZ3NrZdPp3HjJi7/2X3pTkdERKTH0yQLx2aUlkcmjssYsW3MbrnrEc78zMeZctxEZ7FERESiTi14kjG+943TKC0Zyg+u/y2NH25OdzoiIiI9lgo8x+auWsTcVYsiEcd1jJy+ffj5FRewpm4dZ110Aw2Nm5zFEhERiTIVeI7VN2+kvnljJOKkIsYxkw7hf24/m4JPDuf0C6/l/Q3uPzsREZGo0Rg8xyYXjolMnFS9ly+MPZKF2wr4cfWzfP78q3n4zss56MCClMQWERGJArXgOTa+YDDjC9zvMpGKOKl8L2cf5/HAbZewvKaWz31zFmvr1zuPKyIiEhUq8BxbsqGWJRvc70Wbijipfi/HH3MEFb/6Ce/Uredz35xFzTvvOo8tIiISBSrwHFtQX82C+upIxEnHezn2yEP5052Xs/6DjXz2G7NYvsp9gSkiItLTqcBzrDAnj8KcvEjESdd7+cjhpTx61xVsamrms9+YxbLlq53nICIi0pNltbS0pDuHtJo5c2bL7Nmz052GJODN/9Zw+gXXsmNHCw/fcRnjx45Md0oiIiKplpXISWrBkx7jkNHD+cvvZtG3bx8+P/0aXl3qvrtYRESkJ1KB59icqkrmVFVGIk4mvJcxBw/lr7+bRcGA/pxx4XW8uOhN5/mIiIj0NCrwpMc5uLiQv/5uFoUD9+eL//NTnntpSbpTEhERySgag+d4DF5dUyMARbkDnMVIVZxMey/16z7gjAuvY0VNLffe8n2mlk9ympeIiEgG0Bi8TFCUO8B5QZSqOJn2XgoH7s+f77qS0pJizv7ezTw+/yXnuYmIiPQEKvAcm19Xzfw695MBUhEnE9/LwAPy+fNdVzLh0FGcd8nP+dMTvuPsREREMp8KPMeWNtSytMH94rypiJOp76VgQH8euuMyjpowlgt+MocHHlvgMDsREZHM1zvdCUTduHz3e7emKk4mv6GCd7MAACAASURBVJcB/fsx95eXcs7FtzBz1p00NW/hvDM/6SA7ERGRzKdJFlroOFKamrfwjUt+wZPPvcI1F5/Nt84+Jd0piYiIJJMmWWSCuqbGnbNCe3qcnvBecnP6cs8tF3PqScdw5a3387PfPJLk7ERERDKfumgdq6hZDMCM0vIeH6envJe+fXrz65/OJKdvH264vYKm5q385NvTyMpK6I8eERGRHk8FnkRS7969+OU13yI3ty8//92jbG7awrXfP1tFnoiI7BNU4DnmurUrlXF62nvJzs7mZ5edT27fvtz5x8dpat7CTZeeR3a2RiaIiEi0qcCTSMvKyuL6H36N3Jy+zL73MZqat3DbVRfSq5eKPBERiS4VeI7NXbUIgGkjJvb4OD31vWRlZXHFzC/RL7cvN975EE3NW7jjuovo00c//iIiEk36DedYffPGyMTpye8lKyuLH17wBXJz+3L1L/5I85at/PbG75LTt0/SY4mIiKSbCjzHJheOiUycKLyXGV/7HP1y+vLjG+/hSzP+l3tv+T75A/ZzFk9ERCQdVOA5Nr4gNbs/pCJOVN7LN8/6FAP69+M71/wa8/UreHDOpQwbMshpTBERkVTSSHPHlmyoZckG9/u3piJOlN7LuE+M5bZfzeSduvWcfM5lLHp9udN4IiIiqaQCz7EF9dUsqK+ORJyovZf6A7Zh772Gvr1787lvzOLJ515xGlNERCRVVOA5VpiTR2FOXiTiRPG9HDJ6OE/8/jpKS4o5+3s3c3fFk07jioiIpEJWS0tLunNIq5kzZ7bMnj073WlImn24uYnpP57Nk8+9wrfOPoVZ3/3/9u48vqr6zv/46yYhC9kXkiCbrMoO0lZFKwpS5Ku1Y7XaVset1XHG5ae1i3bGtjo6nWoVHRxbbbXVjrbquLR1PnUDSgXUugAqqAUFASUBsocsLMnvj3MTIwYaknNObk7ez8cjjyQ3l/M+35Bzz+d+z/f7PWdrQWQREUlEXbolk85gIkBmRjoP3PZtvnnWPO76zVN847u309i0q7d3S0REpFtU4AVs4bplLFy3LBI5UW9LcnISP/7eBdz47XN5avFfOe3iG9heWRPofoiIiARBBZ5IB7FYjEvOPplf3XIVb/1tI/PPu451Gz/q7d0SERE5KBqDF/AYvPKmOgBK0rMDywgrp7+15bU313H2/7uZPXv38sBt32bmjAmB7Y+IiEgXaQxeIihJzw68IAorp7+1ZcbksTz9wI0UFeRyxj/fxP/+KfjL0yIiIn5QgRewxeXrWVwe/NpxYeT0x7YcOrSEP/363/nMlHFc8v2FLLj3Cfp7r7eIiCQ+FXgBW1Nbxpra4O/+EEZOf21Lfm4Wj971fc5wx3LTnb/jyhvuZvfuPYHun4iISE/oXrQBm5gTzv1bw8jpz21JSx3Az268jBGHFHPrLx/nw7IKfnXLVWRnDQxoD0VERLpPkyy00LEcpAefXMLVN/2CcSOH8Nv/+h5DSot6e5dERKT/0CSLRFDeVNc+Y7Ov56gtnrP/4QR+t/AaNm/dzrxz/403393o786JiIj0kAq8gD2yeTWPbF4diRy15WPHHzWF/7vvBpKTkzjlwh/y/LKVPu6diIhIz6jAE+mmCWOH8/T9NzJ6xGDOvvJmfvXoc729SyIiIoDG4GkMnvRYfUMTF11zB8+98DqXnfdFfnDF10lK0nsnEREJhMbgiYQha2A6v7nt21zwlbncef8fueiaO2hs2tXbuyUiIv2YlkkJ2MObVgFw1vBpfT5Hbdm/lJRkbr72Gxw6tIQfLvgftm6r5DcLvkNhfo4v2xcRETkY6sEL2LbmerY110ciR205sFgsxqXnfpH7br6KN97ZwEnnXcd7H2z1NUNERKQr1IMXsBOKx0QmR23pmlPnHsXg4gLOufJm5px9Lf/0dcc/n+PIy8kKLFNERKQjTbLQJAsJyMYt5Vx/+4P8cdHLZGdlcPFX53PJOSeTn6tCT0REuk2TLBLBWzVlvFUT/P1bw8hRWw5OfXYrV1/3dZY+fDMnHDWFW3/5ONNPvoyb7vwdldXBLxgtIiL9ly7RBmzJtvUATMoN9v6qYeSoLd3LuHzcsdx3y7d4e/0mfvqLx7n9vie557d/4ptnzeNf/vEUTcQQERHfqQcvYMVpWRSnBX9JLowctaVnGePHDOfen1zJC4/ewtxjp/Nfv/4DR5x8GTfc8RA7KmsD3RcREelfNAZPY/Ckl7z73hZu/eXjPPHMCjLSUrnwrC9w6blfZFBBbm/vmoiIJC6NwRNJZIeNHso9P76CFY/dipv9We76zVPMOPlyfnDbb9hWUd3buyciIn2YCryALVy3jIXrlkUiR20JJmPsyCH8/KbLWf7YbZw853P8/MH/Y8Ypl3PdrQ9QvkOFnoiIHLxITbJwzp0MfBcYAPwvsMDM+vc1aOkzxh56CD+78TKuvujLLPjlE9z9kPGrR5/lvDPmcvl5p1I6KL+3d1FERPqIPjkGzzl3PTAb2Bt/6DlgBV5xdybQAFwPNJrZTQfaVtBj8MqbvOUwStKzA8sIK0dtCTfjvQ+2suDeJ3jUXmBASjLnfvlELj//VAYXF/i5myIi0rd0aQxeXy3wlgBzzWxPh8f+AFxtZuvi3ycDq8xs8oG2pUkWkug2bC5jwb1P8PBTfyElOZlzvzyHKy74kgo9EZH+KZoFnnNuEGDABqAEeB24DlgOTOt4SdY5txg408x27LONi4GLAYqKimYsWLAgsP19sXYLAEfnDA0sI6wctaV3MzZ/tIOf//Zpnnz2JZKSkviKO4aLv/oFXboVEelHCgsLu1Tg9cUxeNOBTcA3zKwuXqzdA5QDg4GPOjx3EFCx7wbM7J74v+GKK65oLSwsDGxnN1S+DcApAWaElaO29G5GYWEhP598GNde+jUW3PsEv/vjUh615Zxz2my+durxTB0/kqQkzZsSEZE+2IPXGefc68BtwBy8wq/FOfdNYKqZXX6gfxv0JdrF5d7dDGaXBHdz+7By1JbEytj00TZuv+9Jfvv7P7N7z15KivKZN+sI5h03g+M+N5mM9NRAckVEpFdF9hJtNtBgZns7PLYar2fvCuBsoAXv0u3VZtZwoO1pDJ70dRVVtTy/bCXP/OU1Fq1Yzc6GJjLSU5l15BTmHXcEXzhuBiVFeb29myIi4o/IXqI9FfiCc+4SoAm4CnjRzFqA2+MfCUMzTxMzJyoZAHsyYhw/9wjO+uIsmnftZsVra3l66Ws8s/Q1nl76KgBHTBrDvOOO4KRZn2HC2OHEYl16fRARkT6qzw3YMbMHgaV4Ey0WAVl4PXcJ6ZHNq3lk8+pI5KgtiZexb05a6gBOOHoqP7nmQlbanSx9+Gau/ZczAfjxXY8w66zvMt1dxjX/eR9LXlxN867dge+fiIiEry/24GFm9wH39fZ+iCSyWCzGxHEjmDhuBFdfdDpl26t47oXXeeYvr/Pg75fwy4efISszg9lHT2HecTM48djpFObn9PZui4iID/rcGDy/aQye9EcNjc288MpbPL30NZ79y+uU76giKSnG56YexrxZMzjpuBmMOfQQXcoVEUk80Zxk4TcVeNLftbS0sPrtDfFi7zXefHcjACOHlXLSrBnMmzWDo6YdTkpKcu/uqIiIgAq8rgm6wHt40yoAzho+LbCMsHLUlsTLCCJny9YdPPvCazy99DWWvbKGXbv3kJudyeyZU5l77HROmDmVQQW5vmSJiMhBi+ws2j5lW3N9ZHLUlsTLCCJn6OAiLjxzHheeOY+6nY0sfekNntv7AVW19fzLdf8NwLQJo5hzzHROPGYaR0waQ3Jyn5uvJSISaSrwAnZCcbCLAoeZo7YkXkbQOdmZGZwy50gOrRlBa0sL35pwPIuWr2TR8lUsuPdxbv3FY+TnZnH8UVM48ZhpnDBzKsWFWnNPRKS36RKtxuCJdEtVTT1LXlzNouWrWLxiNdsrawCYOn4UJx4zjROPna7ePRER/2kMXlcEXeC9VVMGwKTc0sAywspRWxIvI6ycv5fR0tLCG+9sbO/de/XNdbS0tKp3T0TEfxqDlwiWbPPuRxr0ST6MHLUl8TLCyvl7GUlJSUybMIppE0Zx9UWnU1VTz59feiPeu7eKJ55ZAXzcuzfnmGnMmDxWvXsiIgFRgRew4rSsyOSoLYmXEVbOwWbk52Zx2ryZnDZvJi0tLbz57kYWLV/F88tWsuC+J7j1l4+Tl5MZ792bzuxj1LsnIuInXaLVGDyRUO3bu7etwhu7N2X8yI979yaN1bp7IiKd0xi8rlCBJ9J7OvbuLVq+ilfe+BstLa3kZA3k+KMmM3vmNGbPnMohJYW9vasiIolCY/ASwcJ1ywC4fOyxfT5HbUm8jLBygspISkpi6vhRTB0/igGzijlxzxxGb0lj0QpvZu4fnn8ZgPFjhjF75jTmzJzKkdMPJy11gK/7ISISNSrwRCRhpKQkc+rcozh17lG0trbyznubWbRiNYuXr+IXv/0T//3AHxmYnsaxn5vInHjv3shhwU5iERHpi3SJNuBLtOVNdQCUpGcHlhFWjtqSeBlh5SRCRn1DE8tfXcPiFd7aexu3lAPePXPnzJzK7GOmccxnJpCZkR7YPoqIJACNwesKjcET6Zve31TG4hWrWLRiNctfWUNDUzOpA1I4+ojxzJ45ldkzp3L46GHEYl16LRQR6StU4HVF0AXe4nJv/bDZJcHetiqMHLUl8TLCykn0jKbmXby88h0Wv/gGi1es4u31mwEYXFzA7JlTOfGY6Rx35CRyszN93WcRkV6gSRaJYE2tdweAoE/yYeSoLYmXEVZOomekp6Uy66gpzDpqCtdfdQ4flVeweMVqFq9YxR+ff5kHn1xCcnISn5k8ltkzp/KZyWMZP3a41t4TkchSgRewiTnhDAAPI0dtSbyMsHL6WsYhJYWcc9pszjltNnv27OW1t9axePlqFq1YxY/veqT9eUX5OUwYO5zxY4czcexwxo8ZzmGjhjIwI823fRER6Q26RKsxeCL9SmV1HWvWfcDadZvaP959bwsNTc0AJCXFGDmslAljhjNh7McfI4YUk5SkW6uJSK/TJdpEoJmUiZkTlYywcqKSAbA7HcZNHsHnPzup/bG9e1vYuKWct9dvYs26Tby9zvv81OK/0vYmODMjjcPHDGP8mOFe8TfO+1yQF+z+ioh0hwq8gD2yeTUQ/GK3YeSoLYmXEVZOVDL2l5OcnMToEYMZPWIwp8w5sv3xnY1NvPv+Ftb+bVN78WdLXuF/nljc/pzSQfkf9/SN8S73jhs5RIsxi0ivUoEnIrIfmRnpHDFxDEdM/HjiR2trK+U7qj/R27d2/SbueehP7Nq9B/AKxsNGDWXahFHxj9FMGDuc9LTU3mqKiPQzGoOnMXgi4oM9e/by/qatrImP63vz3Y2sWvMeO6pqARiQksz4scPbC75pE0YxfvQwBgzQ+2wROShaB68rVOCJSFBaW1v5sKyCVWvfY9Xa91m5xvtcU7cTgLTUAUw6bER7wTdtwijGjRxKcrImc4jIfqnA64qgC7yHN60C4Kzh0wLLCCtHbUm8jLByopIRVs6BMlpbW9m4pby92Fu19n1Wv/0+OxuaABiYnsaU8SOZGi/4pk8czahhpZrBKyJtNIs2EWxrro9MjtqSeBlh5UQlI6ycA2XEYt4yLCOHlfLlk44BoKWlhfUbP2ov+FatfY8HHnueux/cBUB2VgZTx39c8E2bMJrhhwzSbdhEZL9U4AXshOJg72IQZo7akngZYeVEJSOsnIPNSEpKYtyooYwbNZQzTzkO8Mb0vfv+lvaCb+Wa97j7QWP3nr0A5OdmMeXwkYwYWsyQkkKGlhYxZHARQ0uLOKSkkFSN7RPp13SJVmPwRKSPaN61m7fXb45f3n2Pt979gA/LdrRP5GgTi8UoLsxjSGkhQ0sLGVJaxNB48TektIghpYUU5eeoB1Ckb9Il2kTwVo13f81JucHe6imMHLUl8TLCyolKRlg5QWWkpQ5on4zxVs3k9ozGpl18VF7Bh2U72BL/+HBrBVvKdrB2/WaeW7aSxqZdn9hWetoADikpYmhpIUMHxwu/kk9+rVu2ifRdKvACtmTbeiD4k1YYOWpL4mWElROVjLByws7ISE9tX6i5M62trVRW13mFX1lbIVgRLwR3sHjFasp3VLPvFZ3CvGyGDC5i+OBBjBxeyqFDSxg5tIRDh5UwpKRIs31FEpgKvIAVp2VFJkdtSbyMsHKikhFWTqJlxGIxCvNzKMzPYer4UZ0+Z9fuPZRtr2TLVq/4+7BsB1u27uDDsh288/4Wnn3h9faFnMFb12/4kGKv4IsXfSOHlTJyaAnDhxRrUWeRXqYxeBqDJyLyd+3d28LWbZW8v7mMjVvK2bi5nI1bytiwpZwNm8up39nY/txYLMbg4gJGDivp0OtXyshhJYwcWkpO9sBebIlIn6cxeCIi4o/k5CRvosbgIo773KRP/Ky1tZWKqroOBV9ZvAAs57kXXmdbRc0nnl+Ql+31+g39uNdv6CFFFORmU5DnfWgWsEjP6AgK2MJ1y4Dgb6AeRo7akngZYeVEJSOsnKhkdDUnFotRVJBDUUEOn5ky7lM/r9vZyAdbvIJvQ7z3b8PmMl594288+ewKWlo+fSUpc2B6h4Ivi/zcbApysyjIz6YgN9v7Pi/L+3luNvl52WRmpGlmsEicCjwREQlUdmYGkw47lEmHHfqpn+3avYc73vkLzbt2M7kyh8rqeiqra6msqaequo7Kmnoqq+vYuKWcyur69tu8dSZ1QAr5udkU5meTn5sVL/yyKMzLIT83i8K8bAbH1ww8pKRA4wQl0jQGL+AxeOVNdQCUpGcHlhFWjtqSeBlh5UQlI6ycqGSElXMwGXv27KW6dmd7EVhZXUdVTT0V1bVUdfi+srru469r6ti7t+VT2youzOWQfRaJ9r4uZEhJEYMKcnSLOElEuhdtV2iShYhItLW2tlJb38COylo+2lbBh1sr+LD807OFdzY2f+LfpQ5IYUjbQtHxtQG9QjDeC1haRNbA9F5qlfRjmmSRCBaXe2tVzS4J9vZIYeSoLYmXEVZOVDLCyolKRlg5QWfEYjFyszN5rWErDE/nq5+d9anntLa2Ul27c5+1Aj8u/v7y17co2175qfGCeTmZ7XcI8RaJLqS4II/8vI8vERfkZpObnal1AyVUKvACtqbWW9E+6BfhMHLUlsTLCCsnKhlh5UQlI6ycRGhLLBYjPzeL/NwsJncyVhC8y8Nl26s63C2k7c4hFWzZup2XVr6z3zGCsViMvJzM9ski+Xnxz/HJInk5bRNGPvkz3U1EuksFXsAm5gS7Kn+YOWpL4mWElROVjLByopIRVk5faUtKSnL7UjH786fNa9nZ0MzIXVlUVtdTVVNHVW19+9dt4wLLtlfy9rpNVNXUferScEfpaQM6LQrzcjLJSE8jIz2VjPQ0BqankZGR2v7YwPS0T/88PZUBWn6m39AYPI3BExGRXtS8azeV1XVUdywEO8wi9grDtsfrqK7ZSXVtPbv37D3orJSUZNLTUhmYnvqJAtArAj/9WH5uJgV5ORTl51CQ581QbluaJiUlOYDfhnSBxuAlgv460y3Rc6KSEVZOVDLCyolKRlg5/b0taakDGFxcwODigoPKKEwZSGPzLhqbmmls8j43NDXT2NjJY/Gvm5p30dD+847P2UVFdR2NTTvaH9vZ2MzOhqb97kdudiYFedmfLv7ysinMiz+Wlx2/TV42OVkDNSs5RCrwAvbI5tVA8IuRhpGjtiReRlg5UckIKycqGWHlqC3dz8hOySA7MyOQnIXrltHa0sqZBZOorK5jR1VtfPmZWiqqvMvNFdV1VFTV8lF5BW++u5GKqlqad+3udHvJyUmfuFtJW/GXn5tFXm4WedmZ5OdmkZuT2eHrLC1g3U0q8ERERKRTsaQYpYPyKR2U36Xnt7a20tDUTGXVxwVhRbworKyuay8Md1TVsm7jR7y08h2qaus7XaewzYCUZPJyssjLzSQv2ysA8+MFYV5uljd5JSf+eI73fW78c1rqAL9+FX2OxuBpDJ6IiEivaW1tpb6hieqaeqrrvPGFbeMMq2rrqand2f65usPXVTX11NY3HHDbA9PTyMvNIjd7IJkZ6QzMSGNg++c0MjPSyByYzsD0fR4fmB7/ecfnpjMwI530tAG93aOoMXgiIiKS2GKxGNmZ3qXmYQw6qH+7d28LtfUNVNXUU11X7xWJtfEiMV4QVtd6t7hriI8prKz2Zi43NDbT0NTEzoamTu+HvD9JSTEy0r3icGBGevvngRlpnHv6HL409+iD/RUEQgVewB7etAqAs4ZP6/M5akviZYSVE5WMsHKikhFWjtqSeBlh5fQ0Izk5qX39wu5mtLa20rxrt1fwNTbT0NhEfWNTh++9wrChw2M7G5vZ2fjJxxoam2netadb7QiCCryAbWuuj0yO2pJ4GWHlRCUjrJyoZISVo7YkXkZYOYmQEYvFSE9LJT0tlYK87s+wXrhuGeXd/tf+U4EXsBOKg12ZPcwctSXxMsLKiUpGWDlRyQgrR21JvIywcqKSEWZOV2mShSZZiIiISN/RpUkWWnEwYG/VlPFWTVkkctSWxMsIKycqGWHlRCUjrBy1JfEywsqJSkaYOV2lS7QBW7JtPQCTcoO912IYOWpL4mWElROVjLByopIRVo7akngZYeVEJSPMnK5SgRew4rT9z+zpazlqS+JlhJUTlYywcqKSEVaO2pJ4GWHlRCUjzJyu0hg8jcETERGRvkNj8ERERET6IxV4AVu4bhkL1y2LRI7akngZYeVEJSOsnKhkhJWjtiReRlg5UckIM6erVOCJiIiIRIzG4AU8Bq+8qQ6AkvTur46dKDlqS+JlhJUTlYywcqKSEVaO2pJ4GWHlRCUjzBy6OAZPBZ4mWYiIiEjfoUkWiWBx+XoWl6+PRI7akngZYeVEJSOsnKhkhJWjtiReRlg5UckIM6erVOAFbE1tGWtqg1/ZOowctSXxMsLKiUpGWDlRyQgrR21JvIywcqKSEWZOV2mh44BNzAlnReswctSWxMsIKycqGWHlRCUjrBy1JfEywsqJSkaYOV0VqTF4zrkc4BZgKrAduMrMDthfqjF4IiIi0of0yzF4jwG/N7OjgKuB3znn8npzh8qb6tpn1vT1HLUl8TLCyolKRlg5UckIK0dtSbyMsHKikhFmTldFpsBzzh0BbDMzAzCzvwH3Auf25n49snk1j2xeHYkctSXxMsLKiUpGWDlRyQgrR21JvIywcqKSEWZOV0VpDN5o4O19HlsDnLnvE51zFwMXAxQVFVFRURHYTu3ZvRsg0IywctSWxMsIKycqGWHlRCUjrBy1JfEywsqJSkaYOYWFhV16XmTG4DnnZgLnm9nFHR47FxhsZj85wL/bDnwQ8O4VATsCzggrR21JvIywcqKSEVZOVDLCylFbEi8jrJyoZISVs8PMTvq7z2ptbY3Ex/z585Pnz5//8vz58yfHv8+fP3/+S/Pnzx+WAPv2alRy1JbEy4hSW/T7SrwMtaX/ZkSpLVH6fXX1IzJj8MxsL/BV4CfOuRcBA64xs829u2ciIiIi4YrSGDzMbAPgens/RERERHpTZHrwEtw9EcpRWxIvI6ycqGSElROVjLBy1JbEywgrJyoZYeb8XZGZZCEiIiIiHvXgiYiIiERMpMbgiXSFcy4fuAs4FO8YMOBHZuZLd7Zz7lvAqR0eigGlwLFmtt2PjA5Zc4DrgFa8tvzCzB7wM0MkKpxzacAdeLezHAC8CHzbzJp7uN3BeJP8ngGmEV9ntYOhwJfMbE1PcvbJnIZ3a84UvLY8aWY/9Wv70vepwAuQcy4GXAScDmQAzwE3+lVI9BbnXC6w08z29Pa+dNOdwAIz+yuAc+5LQDFQ7sfGzew24La2751zZwDHB1DcpQALgc+bWYVzLhVY7Jx72cze9TFnOvADYBDwPnC1X23peGI0s7XxNlwPHA/sBK41s1cCyMmPf7/azFb0dPv7yTgDuADIAl7Hm9XvWyFhZmt7us8Hm+GcywT2mllTkDl+6STjPwAzs0viP5+F90av28eLc+4sYF78281m9hDwUIeffxb4QU+Lu07aci9whpltiJ9rHnbOzTOzZ/zKcM6NBm4AhgEVeMf++z3YfqfnRLwC1bfj/gA5efh03B8g43R8Pu67SwVesA7DKxxOwethuR84DXjcr4D4H9kVwJeBFry7d1xjZvV+ZeyTVwS8ApwAbAxg+9cDs4G98YeeM7ObfNz+IUCj96X7KdAE3GFmvhR3neQlA9fi/Q34rQWowSu8KoACvL8z35ZRd86NBe4Gvmpm7zvnPod3IjnRzFp6uO1PnBiBtXgDlJcD38fr9XzcOXeOmb3nV45zbhDwz8AWIBPocYHXScY24Ei8O+k0AP8JXEqHwr+nGcBa59w5eCeTJLwF268xszI/Mzr8LB1YAnwX+HN3M/aX45y7CLgQaDsZvmlml/uV4Zxbj1fMbXDO/Rnv+LnfzJZ2NwPAzB7GOyZ+tJ+n3Ah8rycZ+/l/2Q4MBjYA2fGPbi8L1snvqxx4FLjQzFY55w4D7nfOze/B+WV/58RT8fe47yznB8B4/DvuO8v4CvBZfDzue0IFXoDM7B28gxsA59yHgN+9XqPx/pDmmtku59wFeO+ErvY5p63H6GcEUNh1cBxwQoC9g9OBzwOXmtmPnHOlwGPOuQYzWxJA3rnA82a21e8Nm1mLc+5y4EXn3GZgFHCemfm5ivp5wH+0vWs3s78655YAs/BO9t2274nROVcIjDGz8+NP2eqcuwHvBfJbfuXET+hLnXPn460632OdZOwAvtP2c+fcB/Tw9baT31c+3snwFDNrdM6dAPw3Xg+CLxn7WIh/vdyd5ZyEdxlzW0AZY4EZwBNmdrxzLgd4wDkXM7Nf+5G5L+fc8UC1ma3qyXb2BisMbAAABvNJREFU8/u6BFjunKsGhgP/1pOe0E4yTgMeaNt3M3vXOXcPXhHzq25mdHZOzMbn434/OSvN7Hq/jvv9ZDSZma/HfU+owAuBc244Xi9OnZn9wc9tm9l6YH2Hh9Lw3p0E4ad4J5Dzg9h4vGclC3jIOVeC1719nc+9kXnAs2b2PICZlTnnrgIup4cFy77ilxuvwrvs4Lv47+su4Kj4i28J3gnrAzN71aeYFLx3px21AIfj8+8Lr0B9Z5/H1uC90PdJ8R7vq4ESvBOyb8ysCu+YbBPYse+c+w6wFBgT0PaT8d6s3uKcG4E3FOBan3vW8/B6Bf8HwMxqnXP/hNdL9Wsfczr6d7zLeL6KjyX8LXC6mb0UHzbzc+fcBh/PMQc69nuk4zkRrzcykOM+yHPvgTKCPO4PhmbRBiz+AvIT4HYz+26AOXOdcxuBrwE/DGD7FwPrzezPfm+7g+nAJuAbZjYLeBv/1xTajHdZtqOdBPNm5xLgMTOrDGDb4PVEPts23i5+MrwbOMPHjIeA6+In3bbxeKfhTRzx24d4PREdHUoPLjv1JufcqcB9eH8DF5rZroByJjvn3sEbcN/tS5oH2P58oKitMArIYUAl8D0zOx54DK/w8lNnx34DkOxzDtD+/78+3tPjt4nAe2b2EoCZ1eD9/3/dx4zfA990zk0GcM6NwuvR79Gx38k5MZDjPoxzb2cZYR33XaECL0DOm+X0eTP7mp+D3jtjZs8BI4FH+OQMzh5zzh0LTDezO/3c7r7M7FkzO93M6uLf3wMcHr807JeXgeOdcxMBnHMZeN3sD/qY0TYY/WKCHXvxOvCF+LvFtrb8I94YSV+Y2RvAlcBd8XFL5wL/i1eI+8rMPgL2OOdOhPYxX/+K92LZp8Qv/V8KnOZjb2qnzOxNMzscb4zROX5uOz7u6gK8sVGBMbO1Zja7bfygmf0fUO+cG+Zjxiag0HkTK9p6Df8Tryes25xzY5xzT+EVV9c45+5xziXh/X/8qGd7vV9/AyY750bG9yEFb/yin8f+VrwOgx8555YC/4b3u+r2sd/ZOTGI4z6Mc29nGWEe912hhY4D5Jy7Fu9yZsfxV0+Y2R0BZiYDr5rZdB+3+SDeO6rd8Ycm4V0WPid+idivnGygwbz7Crc9thqvuOzRgP59cibgFV4ZeD1395jZ/X5tP57xr0CzBbxsgXNuLt4LbwveJbr7zexunzNy4z0Ebd8/hTfWr0eTOZxzY4DbgXFALV7B+n28YQAj8HpWbjGzR3zOeQNvcPRQIBVYZ2Zf8jkjN77tDR2e9qKZXetjxutmdvE+z3nTzCb7mLESmII3MQm8S7Q7gMvMbJmPOW/Et9nQ4TnP473GdGvSyH7+vm7AWyalCG/m5hPAT83nlQ2cc2cDR5vZZT5tr7O23I1XoKYA6XjLPXV7lYb9ZHxnn2P/58DPzGx1NzM6PSfivcH27bg/QM5cfDru95OxBK/Q9u247wkVeH1c/AQ/CW8maIvzZtWdZGa+vpPfJ/PPwPlmttHn7Z4NfAHv0mYT3vi1cRZf0qAvcc79A95SA41/98kJLF50PwV83cw+jP8fzdq3sJDwxd+onIc3TnVXfED/j+KXOIPK/DXwa7+HasSvElwHnGVm1c65rwAXmFmfvLd4vDdqtfm8NFKYnLdCwyK8CWlvO+dmA1eama9XiCQ4ukTb9y3GezeyyHmzG48EfHnXeABl+D8bGDN7EG8gt+G9sGThLQHT55jZk329uAOIXy7/Ft4A7qV4f19X9u5eCXiXNYF1gDnnFuMVe4G9sYvbwafHsfVYvDfwTuDR+OvY0XjrlfVJZvZ8Xy7uAOI9gRcBN3QYnnFBr+6UHBT14ImIiIhEjHrwRERERCJGBZ6IiIhIxKjAExEREYkYFXgiIiIiEaMCT0RERCRiVOCJiIiIRIwKPBEREZGIUYEnIiIiEjEq8EREREQiRgWeiIiISMSowBMRERGJGBV4IiIiIhGjAk9EREQkYlTgiYiIiESMCjwREZ845zY65w7v7f0QEVGBJyIiIhIxKvBEREREIialt3dARCSqnHOlgAHfMbNFvb0/ItJ/qMATEQlAvLh7GrhGxZ2IhE0FnoiIz+LF3TPAv5rZ0729PyLS/2gMnoiIv0qB54HtZvZUb++MiPRPKvBERPx1H3A50OKcO7e3d0ZE+icVeCIi/rrUzJYA5wPXO+cO7d3dEZH+SAWeiIi/NgCY2UfAd4EHnHN6rRWRUOlFR0QkIGb2KLARuKaXd0VE+plYa2trb++DiIiIiPhIPXgiIiIiEaMCT0RERCRiVOCJiIiIRIwKPBEREZGIUYEnIiIiEjEq8EREREQiRgWeiIiISMSowBMRERGJGBV4IiIiIhHz/wGtwNgPIlNh6wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x242a303b390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# CODE FOR CUSTOM GRAPHICS NOT INCLUDED"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, we fit a K-Means model with *k* = 7 and assign each resale flat to a cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fitting 7 clusters:\n",
    "km_final = KMeans(n_clusters = 7, random_state = 123)\n",
    "km_final.fit(dat_jw_scaled)\n",
    "dat_jw['label'] = km_final.labels_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We find some differentiation in resale prices across the clusters. For example, Jurong West Cluster 3 contains flats that are generally more expensive than Clusters 1, 2, 5, and 6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "dat_jw_full = hdb[hdb.town == 'JURONG WEST']\n",
    "dat_jw_full['label'] = km_final.labels_ + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAo8AAAIECAYAAAB458RgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xmc3VV9//FXJBBERW3iQi0UBVFrXSj9FfQrlVX0U4sbdamiEQVb0QjYahFawKWWgoCpoqQKWooCikWqR0RoUXpEoDUuRVCpIoiAJC4syhKY3x/nXOc7w8zkJJnJ3GRez8cjj9z7Xc/93u/c+75n+X7njYyMIEmSJLV4wGwXQJIkSRsOw6MkSZKaGR4lSZLUzPAoSZKkZoZHSZIkNTM8SpIkqdn82S6AJEkzpeu6pwDPAz6ac1452+XR+tV13YOB/YGrc87/Odvl2VjM8zqPmg1d130MeM0ks28CPgUclXP++XorVKOu6y4Gnl2fHpNzPnqat38t8LurWezLOefdxi372pzzx9Zhvx9j9D35eM558WqWXwycNsGsEeDHwCeBd+Wcb2/c/27Abz7cc87zWtabLWt6vKZ531sAfwHsBzwBeCBwA/BfwIdzzpf1lr2YGTxfpyjjrOx3gnJ8EXgOcCfwtJzz9yZY5lpG/47+Muf84fVXwuHSdd03gKfVpxfknPep0z8OvLpOf03O+V+6rnsocAuwKfBL4FE557u6rjsQWNbb7JNyzlfPYJmvZYLPwa7rHg58H1gI/C/l/b9vpsoxl9hsrWH0aODNwKVd1y2a7cJojc0DtgbeBny5/vLXNOm67g+B7wDvA54B/BYlPG4PLKb83fxj13VDHb7Xh67rdqcER4A8UXDU/fR/EO5eAxiUH/QD+9b/n0sJjgDn5pzvqo9f0lv20pkMjlOplQ/n1qe/D7xyNsqxMTI8ahisBC6r/74BDH4ZPgH4wGwVakj8kNFj0//3ndks1ATuYbRs36LUPgL8AfD3jdu4lbGvUeN0Xff7wH8wtmb6B8CVwL31+Tzgr4HD12/phtLbeo//cdZKsWE5g/L3DCUYDoLiBcAv6uPndF23GfCnvfXOAui67mHAHr3pE7VOrE8n9B7/1ayVYiNjn0cNgwtzzi8fPOm67rXAqfXpS7uue2vO+YbZKdqs+8cNpAntZznnXQZPuq57M7C0Pj2w67rDc853TLWBnPPXgV2mWmYuqzWJpwMPqZN+Brw45/zlOn8r4GOUmrafAl+bhWIOja7rnkSpGQP4Xs75gtksz4Yi57yi67p/B15cJ72E0i3j7q7rPkvpqvEQYC9KX1Io5+KF9fG+jNZG/poaKmdLzvk7Xdf9J7A78NSu6/bMOV80m2XaGBgeNYw+BpwMbE6pRdkZ+MxgZtd1fwQcBvwxsAj4OXAp8MGc85f6G6pfIH8L7A1sSekX9nngvTnnn/SWewjwBkofsifVff8EuAg4Nuf8/dbCd103H3gdpX/Qk4EFwHV1vyfknH/cuq111XXdc4GDgGdS+v38ElgOnJJz/nRvufGdn1/Tdd1rWPu+fP/KaHjcHHh613WbMtqn8SLKl9NxwJ8BhwI/YpI+jzU4vRQ4kFKbuTmlX+XgmF4/7nU/kVL7tjfwKOB24L+BjwKfyjmP9Jb9HeBIypfeIuDmWr735py/2/qCu67brr6e3YFNKDWEbx9so+u6bwFPqYsfmnM+qbfuCxhtXjs75/yyCXbxfODp/ec550sHT3LON3Zd95L6Gg9r+cE17n3fPed8cZ2+LaXWe+CxOedr67wp/6YmWBfgqK7rjmJc38fW92lcf9j7nTs5549P8PL6x/CcKQ/EBMb16f1Rznnb3ryjgaPq0y/nnHer0y9mtJ/nrsBDgXcDD845P74uM6+W7UBgR2ALRs/lY8d9LvX382bgEkqN8rMp3RWuAd6Wc/78uLL/MXAEoz/Gvlqf70E5bjB1H+nTGA2Pz+m67iE559uAsxnt5/vuWgaAz+ScB7WV/Sbrz+Scbx1Xtn2AgyldLgZ9Ji8GTso5XzFu2Z0pf5u7Uv7mfwj8G3Bczvnn4/tJD8redd1p417fOZS/S4CXU84hrQObrbUh6H/RH0KpUXkZsBXlF+4jgRcAF3Rd1/9CfjpwBfAKSijYDHgs8CbgW13X/UFvuasoH6o7U74QNwO2pYTA5V3XNdWI1RD6H8CHKYHtoZQPvR0oAel/az+sGdV13Sa1luALwIsoX8zzKQFyL+BTXde9bwaLcO+45+P7321KaQY7CHj4BPN/o4bxs4EzgT3r8g8EHg8cQnl/ntlb/qXAN4EDKH0vN6N8yT2HUgvyqbpNuq7bmnKOvIHR8+l3KF+Qy2v4bvEYyg+YFwEPo9TMvAC4rOu6p9ZlPtpb/sVjV2ef3uNPMbH+Ol/pB8eBnPPtOeeXzVRNfevfVOO2mt+ncVrPnX6T6nmt5ZpGL6z7fTq1Jq7+gPo0ZTDZHpTyLwC2A5ZQPh+eNcn29qQc+5dR+oVvBvwe8Omu6x4zWKi23Pwn5ThuWf89lxIg/6ix7OdTBi5Sy/cn9fGXKD/WoQTfgUGT9YMZ7WMK45qsu647sW77Tynnz6bAbwN/Dnyt67q39pZ9LpApP5oeWsvxREp4/p/6I6XVZ3uP/9T+wOvO8KhhtB8lcA18G37T+f0ERr8sVlBq0fqX33hL13Wvr4//EXhQffxrSj/BO+vzBwG31cc/pYQqKDVz36b0I6O37D81lv2fKL+SB75Xt7eqPn8o8Jk1GAj0tq7rvjbu39+ubqWc872MHqe7gf+j9EW8s7fYobUPHZQ+hit681bUaf/XWM7xXtp7fC9lpGPfrpSgDqU2YQWTO4pyTgz8nDKCcmAh8Imu6zavNVkfp3yxwmhN60295V8CHF0f/y3lixjKcfoOMGhe3xz41RTl6tsLeARwLeV8Gngoo6HxX+s+ALqu6x7ZW24QHn8FpEn28dTe4682lmu6tfxN3UU5d27rrXdDnfZj+E2N45q8T32rPXdqOHhyfboK+HrDa5tuh1G+Y2+j/DgFeCdjfwT8Evguo/28Hw78W9d1C7m/F1LC1g2MPcc2p9b2dV33WOBDjP1uv4byN7OAUlO7WjnnVZTzdWC/Ov0eSs1f3y2M1v4Fo5/dP6L8kKaWbTHlx97ATyjv+aBm8gHA8V3X7V2fL6XU4FOXuZrRz9HBuTboJz34u4Ly2X1ZLdfg9fyYctyg/JBe3dUstBqGRw2DvXrB6ApKh+2B/8g5X1MfH8loIHof8Oic8x9QPgw+2Fvnr+v/v92btlvO+cmUL/PXAEsGTdG1mejPgZ1yzg/LOT8157wd8Nre+n/Ydd2WU72I+kt4//r0bmDnnPMTcs5PpdRiDgLpwyhNVi0eS/mi7P/brnHdt1E+9B+ac94+5/y0Wo7BF888alNO7a/Yb/r6fM55l5zzuxr39Vu99/B/KDWvA2fmnH8xbvl5lFrRx+acH5dz/txEG601GYf1Jp0NbJVz3oES2D5BGUG5Y875Tsp7P/jyOoty6ZA/oJwLf9PbzpKu6xYw9hzZv54jW1K+qN+Rc/5K4+sHeHPO+bGUGsz+j40/7Lpup3qNwUHT9AMoNZN0Xbc98Lg6PeWcJwusD+k9/skky8y01f5N5ZxvrOdTP7B9pJ5PH6nP1/R96ms5dxZRwhLALTnnuydYZjLTdf26XwOvAhbmnJ9XWyXe0pt/DuUz7ImULjiDEL4I+MtJtnkiJfhsz2gYglILD6UGffC67wH2rs3lj6H0l10T/VrD53Xl8lAw2p9x4AGMZol+k/XH+91DKD/UBt6Sc35M7/P73N68waCW/rn2+Jzzkyjh+hDgDTnnlTnnr9dz7cbesu+q59qYpnyg37VlK7RO7POoYbCQ0Zq/vhuA1wPUL5A/7s17K/DWrusm2t4O9fISX2C09uHk2g/mgpzzv0ywzjeBv61NJQspv1q/PW6ZhzP6K3kiezD6IboZpclysmXXx8CQH1Kaq/6n9se7kzIqd1Vvmd+aaMW1sCmjtUF9P2Bs+Ot7Sc7516vZ7jMofcIG3jy4HEjt9D6+79LevccvA142yXvwEEqT3xcYbZJ7V9d121DOkc8ytqlrdS7IOX+gluu+ruveRgkOg8ucdMD/UGohB7WyLwb+mbHNfJ9mcrcwGhIm+ntZH9bkb2oqa/o+jbe6c6dfMbL5pEuN6n8X3jXpUmvmAznn/g/hZ1C6Wwy8uf7gIeecu3Ld0L+o8/ah9Cns+wHwV/U6hbd1XZcZPZcGl8Pqf0Z+Mud8Yd3+r7uuewOlufhhLYWvA00upzR1bwE8t+u6RPlh1beQEi4voNQ8QgngHxss0HXd4xj9gQTw/q7r3j/JrgefjV9gtMXhk13XnQF8Kec82Xqr0z8Pxner0RoyPGqYvT7nPOh4v5A1O18XAG+nNNn8JbBT/UfXdd+l1FSenHO+t+u6XSl9qAYfLjdTPlz2HrfN1fWTefRq5o8vX4u1umBx13WbUzrX/2Gd9CvKiMinMzaMzVTfn59SahOOzDnfMtECDcERSlPwwK05559OumSxRu9BzvmDXdfdQqnN2IXS7/W4ruuupwS74xvL2a/5IOd8Z9d1/8fo8R90U7iQMnhqG2CPrlxkedBk/WtgwhrY6muUfrTUdY5uKNd0a/qbatjOOv2tNLwnt1DO+S2Ah3dd96A89Wj/fq1u00XtG4wv4/hz+cZx8/tdMR41wfauz2MvcN3f/uDvuL+Pb/ZXrgHylzSGx+o0RvtJvoTSjDw4VonRsPjKWoZBiP1K77Mb1u793p8yeOr1lB/mewB0Xfd14MSc879OvPqktu49vnYN19U4NltrGJxVR9Y+kLFNC2/tPf4lY5uTDq7LD/5t0X+ec74p53xfzvmUnPPTKc2/b6R8oD6B0p9mEMpOYDQ4/lnO+dE5520YHRnbqn83nFWUwNsv44N6j//0fmtPr9cxGlw+DSzKOT+O0iQ7E5dwuTnnPK/371E55zdMFhzXQD8sbtl13eq+hPrvwT8y9viPP08uA8g5n51zfialmWwxJXRvTemf1jrQYswPm65cA69f0/Kzuq/7GG0O3IxSizMYQHX+agLOmb3Hu/T6hvX3u0XXded1Xbfj+HkN+jUzD5xogTX4m1qdNX6f1kQ9zv/VmzTpILXa/7jfJWXQJaD/eTO+9nLC47Ma/b+FLbuue8S4+U+aZNk10T+uT5xg/ppWGJ3JaHP68yl/HwOHUPosQvk8688bf23H8XcKez6TfzY+FMoPsJzzsbXZ/fcozdk/oFxt4fSu645ofRFd1z2N0dr67zb8CNVqGB41NGoTTr9fzF5d1z2nzruDMtJw4BDgEXWd+ygB8CpKX8jBJTH26bru/DrQ5kc55w9R+gQNDJpE+gMRBpdV2ZRySYe+4P4eU5d/AOVyEwPzKTVDq2oZt6L0AfsSpRl3E6Zfvyz913Rdr6bmjyhf9AO79foy9T20bmttviSn06WMDmAB+ECtVaXrumd2XXdO13WLe1/EF/eWPQDYvh7/uyk/Rn5AGVDwhznnka7rfr/ruq905VI5t+RyyZf+oIK9au3g6ryo67pn13JtQjkf+10CLuk9PpXRARLHMFqTM9koawByuYxJv2byXwf7rPv9Hcolrf4UuLzruon6q/bPERi96DPU5vP6no8fiT9Yr/Vvarzx59PFvXmrfZ8m2ebq9PvRvWjwoOu63xt3Xvf7Id5LGVgGY4/No2oAGYw4P6g3b1H9sbA6lzG2Sfz9g/XqCOv9e/P+g7XTH4H/mq5cFoeu6+Z1Xfd26vvYqvZVHhzHLRn9DPyv2md8cD3eBzLanH079+9+cTVjB0O9A9iivufzKX2Xv05pqt+qlvnAruvOqv2Fr8o5v6+uN9B6rsHYQUrnonVmeNSwOZ2xfQ2P7V1Wod8H6PHA97ty7bybKH2FtqV8Ed1e+0h+gNK89x/AT+tAji/2tjEIJT/qTbuiK/d2vZkSZPsf9oP+Rf1+j6/ruu5K4PSc81WMvZ7cm4Gf1DJ+nxLangX88aCv0zTol+Xoruv+F3gvY1/TYV3Xfa82o36Vsc3Wz6Zc6gjGjo59YW2K/HEN0rOi/mjoB5mXADd1XfdtSs3Siym1HF+vXxb/wGifzkXAN+ox+Qnl/Hk0ZXDHoGn1ZMro3XOBn9Umsf75dw+jd9uYyhbAxV3XfYdy7hw89mXk/+k9uY7RCyoPRn3eydRN1gMHMjpQ4pF1nz+o5+C1jDaBz2e0v+WE52t93h9Nf2g9V29g9OLPAx9Yw78pGHs+Lem67ipKv09Y8/dpbXyK0abdV3SjI5jfD/yy67qr6t/Ekb11zsuj92Iff6WB/+q6bjmlKbXf9PtkJh/g8hs5518y9o5Zr6Ccy1cDX2G0ufaXlPNybSxj9JhtBpxfX+N1lGO+Nia6Q8wgNJ7B/fuInj2+Br3+AHhPb9Izgeu6rvsmpcvHCyg1r/sDP6+1wcdSPnP/u+u6n9RzrV+Wyc614+oxPRd+UxEwCPsjrPnAIU3A8KihUpub3t6b9HTKSGhyzv/O2A/6BZSm5cGX5F3Ay3POP6yDKs5mtOlpEaW5o1+LdGz9v9/8sQB4Wt3mxYy9tdXAl3uP51GaVHaqzw8ELu/Nf0Qt46Cm8ZuMXmR3OvTLMp/yRfY0Sq1N/2Lkj6c0o97K5Pd3zeOe70CpPWu9NtxMeRflvRx4KOU+tYMfFT8D/jzn/Ouc83JKH6lB4NuEckwGfcjuAw7Oo9dI/GRv2S0p167rNyf+0xSjn/sGF6l+EmMHs9zIxO/3R8c9/2IuF2GeUs75JsqgiP6ljx5LOQcH59gI5bw9tD6f6nwdX47B39OV4+bNW8O/KRh7Pj2A0oz6pK7rHrMW79Mayzmv6L2GB1IC7AMpP+A2reXpdy24mdFjRs75fxnbZP5gyufRJpQg1L9EUGvf4Xcw9kfCwyk/Kgfr3wbsN0F/yCY55ysZO1odymv8HcrnUn+EdmuN7oWM/Sy5nfr3mMfeO3pgwtsR1gFly3qTHkRpIRn0k/w58MIasn/J2AFrW1HOtUFt4n2MXuwcxp5rm1J/qNfguD+jfS4/V4+R1pHhUUMn5/wFxt414N2D5p2c83soNUWfoXzY302pcfk45XItn+1t5wjKl+RpdZm7KLWUFwJ/knP+p7rcpyjNWt/ube9oSg3LByjNKfcw2hdqaf13M6VT/lXUvnH1w3RXSlPY4Fp3t1NC49uBLud887odoTGOpPyS/llvP1+q5XgWpfnodkoT3DnA/8s5nwOcVMu+oh4Xcs5nUvr53dR7XR9i3GCQ9S3nvCqXO668lPLe/ZxSU/d9Si3SjjnnS3rLf5zyvp9O+dK7u/7/GeCZtal1sOyHKGHqg5Rrct5JOSaZMlq69V64H6CE8m/UbdxEGXCzU855omtljr9O45RN1n0550G/r4Prdn5BqYW5qpZjp5zzW/PoHT8mPV8pX+bHUvqW3kXptvFuysj5Yxm9RumNdd9Nf1PV+yjnz0rKOfgt4HjqNfnW5H1aB8cyWkN1CCWsHVPLcnud913K38NTc84/Grf+n1MG091JCTT/Sbkl5JGUlomf1Xk/aylMLpcMegHl7lNfrtu8m3JlhA8BTxmMkF5bOefjKU3IX6PUvN5Mqcn8U8b+MBrfD3Gy7d3H6I8juH/NYv9HxjU5535f0/HbegPlTk5foJwXd1H+jj8APG3wYyHnfE/O+bWUQTJnUULv4Pz4LLBr/7OeEpjPohzPX1Bqh4+jtAgMukLdx+wMMtsozRsZma5LWkmSptJ13Wa53CP4SEqNKpQv0Efmcbdx0/Soff0OowzMOblx9PwGq+u6P6MM8PnUoHtM13UPooT5Qf/U+4DfnuYfskOn9j9eTGldOj/n/MbZLdHGw/AoSetJ13Xfo9SG9Acu/EvOeTq7MqinNl1u2tj9YIPWdd1fMdqcey/l6hW3UpquH9xbdFmtBZwT5tI5sL4M7XUeI2JPSnXzCKWc/5xSWtML0UrSMNmMscHxx4ztx6tpVpvvWwY9bQw+ShlUtgulb+a2EyxzHmNHmG/05tg5sF4MZc1jRMyn9EnZNaW0MiI2o4zuex2lb8gxwG6U/iqHp5Su6K27mDJoYQQ4JaV0em/ezpSOzg+m3Jni6JTSPXXeVpRbP21b93FISunmOm+zqfYpSS26rvsCpR/XCsoo5aNyztdPvZbUrjbVHkTps/tEyqWgbqFc6uxjOWcvVaN1Nqzh8QGUDuuvTSldHRGPpnQofxGlw3UGPkIZQfUZ4FUppf+rwfHZlH4d84BTgC+llM6IiB0onX5fROlAfCDw/1JKB9ZweAnw5pTS5TVknkQJr6si4mOT7XM9HA5JkqShMZThESAi/pByQeXrKf01XkO5dMpnU0rP6i33PGDvlNJhEXE5sGdK6bY6b0vggpTSLhGxFPj3lNKXeut+lXLR0z2AnVJKR/TmHUsJjHmqfU5Q7oOo15R6zGMes9M//MPaXlpLkiRp/Vm4cGHTZaeGss9jRDyCcmmBXVJK342IRwH/QrnQ6dXjFr+S0Qvybj4IjgAppVsjYnBB5O0ol6jouxrYfpJ5V1Kuc3fjavY5RkppGfVaVkuWLBlZuHDhRItJkiRtkIb1Oo+7UmoMvwtQ+x6eQmlq3mbcstsyej/ke2sTNAARsYDRuxjcMMW6aztPkiRpThnW8Ph14DkRsQggIh5IuUr8F4FVEbFXnb455fpNg1slncXYkYtHUm6fBOWitkfVdYiIvYHbajD9HLBfRGxd521DuZDp51NKP1nNPiVJkuaMYe7zuDcl/N1HuWXcx1NKp9RA+UHKPWE3AY5LKZ1d15lPGRW9N2XAzPmUEdX31vkvB95Kuf7VD4E3pZRW1nlPpwySWUC5Y8BbUkrfqvMm3edUlixZMrJ06dJpOBqSJEkzrqnP49CGx42B4VGSJG1AmsLjsDZbS5IkaQgZHiVJktTM8ChJkqRmhkdJkiQ1MzxKkiSpmeFRkiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1MzxKkiSpmeFRkiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1MzxKkiSpmeFRkiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1MzxKkiSpmeFRkiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1MzxKkiSpmeFRkiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1MzxKkiSpmeFRkiRJzebPdgEmEhGHAfv2Js0DHg08C/gT4EBgBDglpXR6b72dgfcADwYuAo5OKd1T520FnAhsC/wQOCSldHOdtxlwDLAbcAdweErpit52F0+2T0mSpLlkKGseU0onpJR2G/wD/gn4EiU4PhvYHdgT2DMiXgkQETsAJwGvAp4B/Ag4uc7bDDgXOCGltEtd7tyIGITnZcAPgGcC+wNLI2K7uu7iyfYpSZI01wxlzWNfRGwCHA48H/gssGdK6e46bwlwAXAG8Cbg71JKN9VVl0XE4oh4GLAHcGFK6XKAlNJlEfGVsonIwPYppcV1vRsj4p3AwcBhwBun2OdE5T0IOAhg0aJFrFy5cvoOhiRJ0gxZuHBh03JDHx6BV1OC340RsXlK6bbBjJTSrRGxRX26HXDVuHWvBrafZN6VwA7AjXW58fMOro+n2uf9pJSWUWoyWbJkyUjrGyFJkrQhGMpm64Ha3HwocGyddG+dNpi/AFhVn94AbDNuE9sC16/DvNXtU5IkaU4Z6vAI/AVwTkrpZ/X5WcCRvflHMtp8fBpwVERsDhARewO31UExnwP2i4it67xtKANyPp9S+gmwKiL2qvM2B44ATm3YpyRJ0pwytOExIh5E6Tt4Qm/y8cAmEXF5RAxGQ58EkFK6lBIgL4mIrwEH1H+klG6tj0+PiEuBjwMHpJRur9t4NXBgXe8S4KO90daT7lOSJGmumTcyMjLbZZhQRBwB3JVSOn62y7K2lixZMrJ06dLZLoYkSVKLeS0LDfOAmSuBL852ISRJkjRqaMNjSunc2S6DJEmSxhraPo+SJEkaPoZHSZIkNTM8SpIkqZnhUZIkSc0Mj5IkSWpmeJQkSVIzw6MkSZKaGR4lSZLUzPAoSZKkZkN7hxlJkrR+nbd8xTpvY98dF01DSTTMDI+SJAlYffA7b/kKw6FstpYkSVI7w6MkSZKaGR4lSZLUzPAoSZKkZoZHSZIkNTM8SpIkqZnhUZIkSc0Mj5IkSWrmRcIlSRsN75AizTzDoyRpo+EdUqSZZ7O1JEmSmhkeJUmS1MzwKEmSpGaGR0mSJDUzPEqSJKmZ4VGSJEnNDI+SJElqZniUJElSM8OjJEmSmhkeJUmS1MzwKEmSpGaGR0mSJDUzPEqSJKmZ4VGSJEnNDI+SJElqZniUJElSM8OjJEmSmhkeJUmS1MzwKEmSpGaGR0mSJDUzPEqSJKmZ4VGSJEnNDI+SJElqNn+2CzCZiHg4cDKwLaWcCTg6pTQym+WSJEmay4Y2PAIfAE5MKV0OEBEvAB4ZET8F3gq8GLgHeG9K6fzBShHxJ8DbgE2BT9dtjNR5OwAnAIuA5cDbUkq31XlbAscBTwNuAQ5NKV1T582bap+SJElzxVCGx4j4beDX5WEcD9wJvD+ldHNEHA1sAuwKPAg4KyLuSCldEhG7A28Cng/8CjgGeAfwnohYCHwSeHlK6fsR8XzgbOB5dbfnUILmG2rIPDMi9kop/QI4arJ9rofDIUmSNgDnLV+xztvYd8dF01CSmTWU4RHYkRLUDk4pHR0RjwbOiYhfAfsBT00p3QfcGhGHAe8BLgEOBZb0ahOPAr5R578G+HBK6fsAKaXPRcT+EfFUynH4aUop1Xnfi4iPAq8Glq5mn2NExEHAQQCLFi1i5cqVM3F8JElryc/ldePxm1y3zbwp5+frRla7zGwe34ULFzYtN6zh8WHABSmlCwFSSjdFxKGUmsSbaogbuBp4XH38u8A1gxkppXsj4paIWARsR6l57LsS2AGYB1w1wbxFo4ORAAAgAElEQVSXRsQjVrPPMVJKy4BlAEuWLBlpfSMkSevBdSuavyA1AY/futlIjt+wjra+ntJU3XcH8AvgkeOmbw3cWB/fDGw1bv4jgJXADcA24+ZtW/c11byVq9mnJEnSnDGs4fEyYLeIeDJARDwQeDdwBvDViDigTn8A8C7gI3W9Uyn9Gx9Q578euLgOmPkE8JY6ipvaXP1E4Iq6v6dFxFPqvIcDbwQ+UWscp9qnJEnSnDGU4TGldBelj+L7IuLLwIXAuSmlz1FGPe8cEZcBlwP/nVI6p653JvBN4LI6fyfg7XXetZTBMykiLgXeC7wipXRfSule4OXAsXVeAv4mpXR9LdKk+5QkSZpL5o2MeNnEmbJkyZKRpUuXznYxJEnVectXbBCjWYeVx2/dbADHb+rRPNVQ1jxKkiRpOBkeJUmS1MzwKEmSpGaGR0mSJDUzPEqSJKmZ4VGSJEnNDI+SJElqZniUJElSM8OjJEmSmhkeJUmS1MzwKEmSpGaGR0mSJDUzPEqSJKmZ4VGSJEnNDI+SJElqZniUJElSM8OjJEmSmhkeJUmS1Gz+bBdA0sbnvOUr1nkb++64aBpKIkmaboZHSdNudcHvvOUrDIeStIGy2VqSJEnNDI+SJElqZniUJElSM8OjJEmSmhkeJUmS1MzwKEmSpGaGR0mSJDUzPEqSJKmZ4VGSJEnNDI+SJElqZniUJElSM8OjJEmSmhkeJUmS1MzwKEmSpGaGR0mSJDWbP9sF0Mw5b/mKdd7GvjsumoaSSJKkjYXhcSO2uuB33vIVhkNJkrRGbLaWJElSM8OjJEmSmhkeJUmS1MzwKEmSpGaGR0mSJDUzPEqSJKmZl+qRJuA1MiVJmpjhUZqA18iUJGlihkdJkuaIL357JXetGlmnbaxLy8yC+fPY5ykL12n/mn0bXHiMiD8B3gZsCnwaODGlNFLn7QCcACwClgNvSyndVudtCRwHPA24BTg0pXRNnTcPeCvwYuAe4L0ppfNb9ilJ0oZiXYPjhr5/TY+hDY8RcQywB3BvnfQl4KvAm4DnA78CjgHeAbwnIhYCnwRenlL6fkQ8HzgbeF5d/xxK6HtDDZlnRsReKaVfAEcBmwC7Ag8CzoqIO1JKl0TE7pPtc2aPgCRJ0vAZ2vAI/DGwe0pp1WBCRJwHLOnVJh4FfIMS5F4DfDil9H2AlNLnImL/iHgq5XX+NKWU6rzvRcRHgVcDS4H9gKemlO4Dbo2Iw+o2LwEOnWKf9xMRBwEHASxatIiVK1dO5zGZdsNevmHmsVs3Hj/NFs+92TXXj/8wv/6FC9u6FAxleIyIRwAPBj4REY8Cvg78LfC7wDWD5VJK90bELRGxCNiOUvPYdyWwAzAPuGqCeS+t+7qpBseBq4HH1ceT7jOldL+OHymlZcAygCVLloy0vhGz4roVzSeKxvHYrRuPn2bLXD/3rlv3K0msq7l+/DeG1z+s13ncEbgOeF1K6dmU4LcMuBnYatyyjwBWAjcA24ybty1w/WrmrQQeOW7e1sCN9fFU+5QkSZpThrLmMaV0AXBB7/myiPgLymCY90TE61JK90XE64GLU0ojEfEJSl/FL6SUfl6bq58IXEGpeTwxIp6SUvp2RDwceCPwZ3U7X42IA1JKp0bEA4B3AR+puz91sn2ur+MhSZJmn6PVi6GseYyIh0TEJuMmbwJ8AvgmcFlEXAbsBLwdIKV0LWUgS4qIS4H3Aq9IKd2XUroXeDlwbJ2XgL9JKV1ft/1WYOe6zcuB/04pnVO3e+Zk+5QkSXPHbI8Wn+39DwxlzSOwL/CcWtt4J2XQyqW1X+JJ9d/9pJQuAi6aZN4PgZhk3h3AGyYrTEpp0n1K0nTzDkeShtlQhseU0hkRsYBSQzgCXAwsmdVCSdJ64h2OJA2zoQyPACmlUyn9DSVJkjQkhrLPoyRJkoaT4VGSJEnNDI+SJElqZniUJElSM8OjJEmSmq1ReIyIJ0bEfjNVGEmSJA235vAYEScBHwNOqc93johDZ6hckiRJGkJN4TEi3ghsnlLaBbi7Tr6KcheY185U4SRJkjRcWmseD6HcNxrKHV9IKd0KvB44bAbKJUmSpCHUGh4fllL62fiJKaUbgEdNb5EkSZI0rFrD4y0R8YjxEyNiJ+B+oVKSJEkbp9bw+PfAcf0JEfFk4F+Bd093oSRJkjScmsJjSukM4KqI+Cbw8Ij4b+BS4LSU0r/OZAElSZI0POa3LphSOjYi/gV4JmXQTE4p3TxjJZMkSRoiBxxw8GwXgRXLz5rtIrSHx4h4Wkrpm8A59fm8iNgzpXTRjJVOkqSeL357JXetGlmnbZy3fMVar7tg/jz2ecrCddq/NlynnvrB2S7CUGi9zuPeQIqILXqTFwKnRcQ+M1IySZLGWdfguKHvXxoGrQNm3gm8MqX0q8GElNIK4JXAUTNRMEmSJA2f1vC4bUrp4vETU0qXAI+b1hJJkiRpaLWGx3sj4qHjJ0bEbwGrprdIkiRJGlat4fFfgH/oT4iIecCJwEenu1CSJEkaTq2jrY8GPh8RlwKfAzYBXgj8GDhgZoomSZKkYdMUHlNKdwN7R8R+wLMp13k8JqX02ZksnCRJkoZL83UeAVJKnwY+PUNlkSRJ0pCbNDxGxMkppTfWx5+k1DZOKKX05zNQNkmSJA2ZqWoeL+49Pn+GyyFJkqQNwKThMaV0du/pDSmlC9dDeSRJ0gzx3syaDq19Hj8SEb/Xv8OMJEnasHhvZk2H1us8ngAcPpMFkSRJ0vBrrXlcCbwwIj7B2L6QAKSUlk1noSRJkjScWsPjXsCP6uNnjJs3AhgeJUmS5oDWi4S/dqYLIkmSpOE3ZXiMiO2AQ4EdgGuA96eUvrs+CiZJkqThM+mAmYjYhdK/8VrgxPr/RRHxrPVRMEmSJA2fqWoeTwH2SyldVp9/ISK+ApwM/MGMl0ySJElDZ6pL9WzVC44ApJS+Bvz2zBZJkiRJw2qq8Lhqkun3zURBJEmShtmC+fPm9P4HWi/VI0mSNKft85SF67T+ectXsO+Oi6apNLNnqvC4ICL2BsbH3M3GT08pXTAThZMkSdJwmSo8fht4R8P0EcDwKEmSNAdMGh5TSrutx3JIkiRpAzDVgBlJkiRpDAfMSNJ69sVvr+SuVSPrtI3zlq9Y63UXzJ+3zh3/Jc1d1jxK0nq2rsFxQ9+/pA2bNY+S1pg1Z5I0d61ReIyIfYAnp5ROmKHyTLbfhwJ3pJRWRcTOwHuABwMXAUenlO6py21FuQ/3tsAPgUNSSjfXeZsBxwC7AXcAh6eUrujtYzFwIGX0+CkppdN78ybdpzQXzXbN1WzvX5LmsqbwGBGbAmcD9wC7AydExJ7Avimlt8xg+YiIRcAVwO41AJ4EvAi4mRL2TgYOrPPOBd6cUrq8Br5zI2LXlNIqYBmQKZcZejTwmYh4VUrp/2pwfHZ9bfOAUyLivpTSGRGxw2T7nMnXLUm6vwMOOHi2i8CK5WfNdhGkWdVa83g48K2U0lERcWOd9mVg/4j465TScTNRuIiYD3wIuLZOehPwdymlm+rzZRGxOCIeBuwBXJhSuhwgpXRZRHylbCYysH1KaXFd78aIeCdwMHAY8EZgz5TS3XW/SyjXrjxjqn2mlH4xQZkPAg4CWLRoEStXrpyuwzEjhr18w8xjN7s8/utmQz1+p576wdkuwgZ77IbFXD9+w/z6Fy5s6w7UGh4PAJ5cH48A1Cbkw4BLgRkJj8DxwAeBxfX5dsBV45a5Gth+knlXAjsAN9blxs8b/ITdPKV022BGSunWiNiiYZ//Pb7AKaVllFpOlixZMtL6RsyK61Y0nygbm+nos5evW/v1N/g+e9etfX/F6bJBn7sev7XnsVs3Hr/ZtZF877aOtt4ipXTHBNN/DszIUag1eNeklC7uTb4B2GbcotsC16/DPIB7a7P3YN8LgFUN+9QGarb7zM32/iVJWlut4fH6iPjd+rh/r+s9mIEQFRHPAnZMKX1g3KzTgKMiYvO63N7AbXVQzOeA/SJi6zpvG2Bf4PMppZ8AqyJirzpvc+AI4NS63bOAI3v7OZLSZL26fUqStMFYMH/e6hfaiPev6bEmfR4/GBH7UZutI+K5wEcoTdrT7S+BbSPi4vr894EnAq+ihLlLIuJeyojqA+A3Tc0HAKfXmsM7gQNSSrfXbby6voZ3A5sAx/VGWx8PHBMRl1PC8fmUQTKklC6NiAn3KUnShmRdu8uct3wF++64aJpKow1VU3hMKV0QEQsp/QQXRcQtwL3AW1JKF0x3oVJKr+w/ryFycUrpWuAa4MxJ1vsG5VI8E81bAbxsknmrKDWRR0wy/8zJ9ilJkjSXNN9hJqX0SeBJwDOAfYCtU0pn1VA5025itA+iJEmSZskaXSS8XsrmfwbPI2Ie8DXg8dNcrvH7fflMbl+SJEltJg2PEfGchvW3A+z8IEmSNEdMVfN4eMP6d1Muoi1JkqQ5YNLwmFLafX0WRJIkScNvjfo8jhcRDwBemFL6zDSVR5I2et6fWdKGrCk81lv1/QPwYmCr3qw7Kfe4NjxKUqNhuD+zJK2t1kv1vJdyG8JdgA8Du1MGynwReOfMFE2SJEnDpjU8vgA4MKX0Y0pw/GpK6efAYcBRM1U4SZIkDZfW8PjAlNKv6uMt6x1ZqHd8+b2ZKJgkSZKGT2t4/EZEPKE+vjYidgSIiEcBm85IySRJkjR0WsPj3wODu7wsBT4REX8NJOCjM1EwSZIkDZ+m0dYppS9TRlWTUjo7Iu4C9gBOSCmdMYPlkyRJ0hBprXkkIjbpPU2UkdbfmPYSSZIkaWg1hceIeAlwbG/SR4G/Ai6IiJfORMEkSZI0fFprHt9Bub4jEfEI4CEppT2A59R5kiRJmgNab0+4dUrpmvp4D+AcgJTSlTVMahZ88dsruWvVyDpt47zlK9Z63QXz57HPUxau0/61YfL2epI0d7WGx9sj4iEppduAPYG/A4iI+cCqmSqcprauwXFD379mj7fXk6S5qzU8fhI4JyL+E3hcSummOj2oo7AlSZK08Wvt83gk8G/AI4GDetNfABwz3YWSJEnScGq9zuMI8CH4zV1lBtNfN0PlkqSN1oL582a128eC+fNmbd+SNnytzdZExIuA9wEPBx4eEY8FHplSumymCidJG6N1HWh23vIV7LvjomkqjSStmabwGBF7AocCOwHf6a37sYh4XUrpqzNUPmlGOFpYkqS101rz+PfAa1JKP4+IEYCU0vcj4tXAccBuM1Q+aUY4WliSpLXTOmBm25TS1eMnppSuAJ4wvUWSJEnSsGoNj3dFxGbjJ0bEVoAX+5MkSZojWsPjx4A39idExObAKcBp01wmSZIkDanW8PguYK+IWAY8KCKOpwyceQBw1EwVTpIkScOl9TqP9wDPj4gXA3cBmwGHppQ+O5OFkyRJ0nBpvs4jQErpM8Bn+tMi4oEppV9Pa6kkSZI0lKYMjxHxEuAdwA7AD4DjU0qn9+bvC5wEPG4mCylJkqThMGl4jIjFwCHAEuDbwNOAEyJiE+Ai4J+AxwOvnfliSpLkrR2lYTBVzeMxwLNSStfX5xfXWxR+DdiEUuP4Z7U/pCRJM85bO0qzb6rwuFkvOAKQUvpRRGwK/FFK6QczWzRJkiQNm6ku1TNZu8DdBkdJkqS5qfU6j5IkSdKUzdZbRMRBLdNTSsumt1jSzLLTvSRJa2eq8PhvwDMapo8AhkdtUOx0L0nS2pk0PKaUvASPJEmSxrDPoyRJkpoZHiVJktRsje5tLUmSpImdt3zFOi+zIfSnNzxKkiRNgw0h+E0Hm60lSZLUzPAoSZKkZoZHSZIkNTM8Slpjs32HnNnevyTNZRvkgJmIWAwcSLm7zSkppdN783YG3gM8GLgIODqldE+dtxVwIrAt8EPgkJTSzXXeZsAxwG7AHcDhKaUrWvYpzTXeoUeS5q6hDY8RMQ9YArwYuA+4EvgbYD/g2cDuwDzglIi4L6V0RkTsAJwEvAi4mRL2TgYOrOHwXODNKaXLa8g8NyJ2TSmtotxiMQPvAB4NfCYiXpVS+r8aHCfc53o5GJIkSUNiaMMjsB3wK2DvlNLdEfFaSs3grsCeKaW7ASJiCXABcAbwJuDvUko31W0si4jFEfEwYA/gwpTS5QAppcsi4itlE5GB7VNKi+t6N0bEO4GDgcOAN06xzzEi4iDgIIBFixaxcuXKaT0ow2Zjf31TmcuvfTp4/CaXrxtZ7TKru1Zct41N+5Px3Fs3Hr+N18KFba1KQxseU0rXANf0Ji0AHgRsnlK6rbfcrRGxRX26HXDVuE1dDWw/ybwrgR2AG+ty4+cdXB9Ptc/x5V5GqcVkyZIlI61vxFq5bvUXI51pM/r6htl1K+bua58OHr8p7euhmTmee+vG4yeGODwORMTewD8DPwJeCpwfEZv1agEXAKvq4jcA2wA/7m1iW+D6Om/bcZvfFvhub73x866vj++dYp+SJG0U5sodUrRuhj48ppS+FBGPpTQd7wucBRwJ/F1d5EhGm49PA46OiBeklO6swfO2lNLNEfE54OKIOD2ldH1EbFO3t1tK6faIWBURe6WULoyIzYEj6rZZzT4lSdooGPzUYoO4VE9KaQT4MCVAHg9sEhGXR8RgNPRJdblLKQHykoj4GnBA/UdK6db6+PSIuBT4OHBASun2uo1XUwbWfA24BPhob7T1pPuUJEmaS4a25rHWGv4+8P6U0n3AK4Ar68joI+q/+0kpnQmcOcm8b1AuxTPRvBXAyyaZN+U+JUmS5oqhDY/AfwB/AFwUEQD/SxlNLUmSpFkytOExpXQvcGz9J0mSpCEwtOFRq3fAAQevfqEZtmL5WbNdBEmStB4ZHjdgp576wdkugiRJmmM2iNHWkiRJGg6GR0mSJDUzPEqSJKmZ4VGSJEnNDI+SJElqZniUJElSM8OjJEmSmhkeJUmS1MzwKEmSpGaGR0mSJDUzPEqSJKmZ4VGSJEnNDI+SJElqZniUJElSM8OjJEmSmhkeJUmS1MzwKEmSpGaGxw3Ygvnz5vT+JUnS+jd/tgugtbfPUxau0/rnLV/BvjsumqbSSJKkucCaR0mSJDUzPEqSJKmZzdbSBM5bvmKdl7FLgCRpY2R4lCZg8JMkaWI2W0uSJKmZ4VGSJEnNDI+SJElqZniUJElSM8OjJEmSmhkeJUmS1MzwKEmSpGaGR0mSJDUzPEqSJKmZ4VGSJEnNvD2hpGnnvcElaeNleJQ07Qx+krTxstlakiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1G9rrPEbEPOBA4CXAA4EvAe8GNgWOAXYD7gAOTyld0VtvcV1vBDglpXR6b97OwHuABwMXAUenlO6p87YCTgS2BX4IHJJSurnO22yqfUqSJM0Vw1zz+ATgkcDzgT2AJwIvApYBPwCeCewPLI2I7eA3wfHZwO7AnsCeEfHKOm8H4CTgVcAzgB8BJ9d5mwHnAieklHapy50bEYNwPek+JUmS5pKhrXlMKV1NqWkEICJuAB4CbJ9SWlwn3xgR7wQOBg4D3gjsmVK6u66zBLgAOAN4E/B3KaWb6rrLImJxRDyMEk4vTCldXvd9WUR8pWwi8mr2OUZEHAQcBLBo0SJWrlw5Lcdjpgx7+SRpuvm5J01s4cKFTcsNbXgciIhtgMOB24DvAFePW+RKSpAD2DyldNtgRkrp1ojYoj7dDrhq3LpXA9tPMu9KYAfgxtXsc4yU0jJKTSVLliwZaX0jZsV1K5pPFEnaKPi5J62zYW62JiLeABwLnJRSehtwA7DNuMW2Ba6vj++tTdCD9RcAq+rTqdZd23mSJElzytCGx4h4OrBrSukVKaXvAqSUfgKsioi96jKbA0cAp9bVzgKO7G3mSEqTNcBpwFF1HSJib+C2Oijmc8B+EbF1nbcNsC/w+YZ9SpIkzRnD3Gz9POD/RcTFvWn/Brwa+GBEvBvYBDiuN/L5eOCYiLgcmAecTxn8Qkrp0og4DbgkIu6ljKg+oM67NSIOAE6vtZV3AgeklG6v251qn5IkSXPGvJGRkdkuw0ZryZIlI0uXLp3tYkzqvOUr2HfHRbNdDElab/zck6Y0r2WhoW22liRJ0vAxPEqSJKmZ4VGSJEnNDI+SJElqZniUJElSM8OjJEmSmhkeJUmS1MzwKEmSpGbDfIcZraPzlq9Y52W8mK6kDYmfe9LMMzxuxPwAlDTX+LknzTybrSVJktTM8ChJkqRmhkdJkiQ1MzxKkiSpmeFRkiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1MzxKkiSpmeFRkiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1MzxKkiSpmeFRkiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1MzxKkiSpmeFRkiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1MzxKkiSp2fzZLsDqRMRWwMuBL6aUvhMRmwHHALsBdwCHp5Su6C2/GDgQGAFOSSmd3pu3M/Ae4MHARcDRKaV7evs5EdgW+CFwSErp5jpvyn1KkiTNFUNd8xgRL6OEvacAv1cnLwN+ADwT2B9YGhHb1eUXA88Gdgf2BPaMiFfWeTsAJwGvAp4B/Ag4uc7bDDgXOCGltEtd7tyImL+6fUqSJM0l80ZGRma7DKsVEUcD/wv8J/DZlNKzevOeB+ydUjosIi4H9kwp3VbnbQlckFLaJSKWAv+eUvpSb92vAgHsAeyUUjqiN+9YINd/k+5zgrIeBBwEsGjRop1OPPHE6ToMkiRJM2bhwoXzWpYb+mbrcR4HXD1u2pXAwfXx5oPgCJBSujUitqhPtwOuGrfu1cD2k8y7EtgBuHE1+xwjpbSMUlPJkiVLRhYuXLialyRJkrThGOpm6wncAGwzbtq2wPX18b21CRqAiFgArGpYd23nSZIkzSkbVHhMKf0EWBURewFExObAEcCpdZGzgCN7qxwJnFEfnwYcVdchIvYGbquDYj4H7BcRW9d52wD7Ap9v2KckSdKcMdR9HiNie8rglR2AW4GvA+8APgj8LrAJcFxK6ey6/HzKqOi9gXnA+ZQR1ffW+S8H3grcSxlR/aaU0so67+l1XwuAO4G3pJS+VectmmyfU1myZMnI0qVL1/1ASJIkzbymPo9DHR43dIZHSZK0AWkKjxtUs7UkSZJml+FRkiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1MzxKkiSpmeFRkiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1MzxKkiSpmeFRkiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1MzxKkiSpmeFRkiRJzQyPkiRJamZ4lCRJUjPDoyRJkpoZHiVJktTM8ChJkqRmhkdJkiQ1Mzz+//buP9iOurzj+DuShDIoAjID0qoZRelYGRS0oDhVazoyjyh2OkQKSovTSZXWQYsgBYGgFEJBQESpIDSdTC1haEto+ZRS+SFV8EeKA8ZSC1Tll3SCFX8UpEm4/WP3Mpc7N8nee8695258v2buTM7ePbvPPtnzvc9+9/vdI0mSpM4sHiVJktSZxaMkSZI6s3iUJElSZxaPkiRJ6sziUZIkSZ1ZPEqSJKkzi0dJkiR1ZvEoSZKkziweJUmS1JnFoyRJkjqzeJQkSVJnFo+SJEnqzOJRkiRJnVk8SpIkqbOFow6gL6pqF+A8YH9gA/DhJPeNNipJkqS5Zc9jd38LrE1yMHACcFVV7TrimCRJkubUgrGxsVHHMO9V1QHACUmOnrDsA8CiJBdPWnc5sLx9uS/wnTkLdPr2AB4bdRA9Ze4GY/4GY/5mztwNxvwNZr7n77Ekh25rJW9bd/My4J5Jy74NLJu8YpLLgMvmIqhBVdW6JK8ddRx9ZO4GY/4GY/5mztwNxvwNZnvJn7etu3kYePGkZUuAB+c+FEmSpNGxeOzma8D+VbUfQFXtBhwHfGGkUUmSJM0xi8cOkmwGjgTOrao7gAAnJ+l7z2Mvbq/PU+ZuMOZvMOZv5szdYMzfYLaL/DlhRpIkSZ3Z8yhJkqTOLB4lSZLUmY/q0Xarqo4DjgI2ATsB1wDnJxmrqvVJXjXAtn8f2CPJ+UOI88XAZ4HdaT6Tf5XkM4Nud1B9yd98Zg41Kp57mk0Wj/OcDcDMVNXFwBPAm5JsrqqFwDnAHwGXjDS4CapqAXAp8KEk97bL3lNVOyV5coRx9SV/i2i+8ek3gV8C/jrJ50YbVaMvOQSoqpcCRwCrkzwy6njAtm8QfTn3qurPgV+fsGgRsCNwSJKnRhOV514XFo/zWF8agPmmql4N7JvkbePLkmyqqpOYYqhGVa0CLkmyrn29pH19WLutC4EFwA7ARcB/AicDi6rqoCRHVNVi4AzgwHa9HwKnJPmvqloBfBc4HFiX5OwJuz8I+CbwwXZfPwH+bMSFY5/y90bgR8DbaBr569uH8P7b0BIyA33KYVUdT/NFCEuAm4CRF4+2fTPXp3MvyUmTYvlIu3yUhWMvz72qekGSH87V/iwe56k+NQDt/q6g+QM0bnWSKwbPxIwsBa6dvDDJGLB5mts6DTg+yd1V9RxgvyTfqqqVPPvqcQVwQ5JTAarqhcDlVfXb7e+PBo5KMvlrqQ6kKXz+MMmdVbUP8HdVtSzJf0wz1mHpTf6S3ALc0r58oqo2ABunGeNs6FMOP9Wuv2qacc2KPrV9VbVju62lNL1UtwNnJNk0xJRMV2/OvYmqahfgGJoL6pHo07k3KY59gK/TDH2aExaP81dvGoD21uvLk/zGNOOaLRsZXgFxNbCszcWXk9y1hfWOBN5QVbfz8NYAAAbbSURBVBOX7Qa8ov33X26h4dwVuCrJnQBJ7quqs2hyfdowDmAG+pQ/AKrqlcCpwB1J7h487IH1LofzSG/aPpo/+N8CzgKepikCPgB8eppxDlNfz70TgEtHedeFfp17tOvvAnwS+P404xuIxeP81acG4ABgY1WtpSmGbgbOTjKqHqAAn6uqK9oP/TOq6nlJfjpp/Y00Y23GLX5mQ8kaYE1V/SpwTFW9JMnJU+xzAXBYkp9N/kVVHcGW/y8fBJ43adn/MtrPZp/yR1WdCewJnJTk4W0c21zpVQ7nmd60fUluH/93WyAsAnYeRuAD6N25V1V7AO8CRv2dz7059+CZc+5S4E9pJl3OGR/VM38FOKrt1XuWqppcbMA2GoAkH6O5MjmmvfKZyngD8OYJP/sn+faEfUxlf+BeYBnwljaOM7dybLMqzcST24BL2lsCVNUOVXUqzVUewMJ2LAs0twUOGl8P+OD4tqrqsKpa1N5CvgI4vJpJGmM0EzTGXQVcOL7NqnpjVX2iQ7j/BBxbVb/Svm834KM0Dc9I9Cl/VfUOYHGS98+jwrFXOZyH+tT2jcf1XuBhYF9GPC6up+feKcAFI+xwGNe3c+/jwJok/76VdWaFxeM81acGIMmV7R/vp5I8DZwOvH2mxz4MSVYA3wBurapbaCYCPAGc2K7yFeBrVXUo8HngXVV1LbCq/dm7XW8R8C9VdRvN/8eFbQP3DeD3quqmqnoNzTE/Cvxru+5xNONdthXnBuB4YHVVfQm4DliZ5JuDZWAwfckfcBhQVXXrhJ93D3b0w9GXHFbVIVX1j8BbaT7/Zw167IPoU9s3IebVwIuAe4CRD9/py7kH0F44vxlYPcAhD0Wfzr2qOhL4eZLrZnzAA/DrCee5aqb1LweeorlCWQt8KsnT1UxSeTXNWK87aXqrHgd+SjM49/IkB1QzduJ4mouFxcCVSS6rZpzYWuAB4CPAeppGYCnN1c5DwB8n+Z9qBu6uT3LNFDHumuTxCa8XAV9P8pph50PSL4Y+tH1TxLwEuDjJO4eShF8AVXU58M9d8jtX+nDutZ0N0BSj0NzyX8cWhg8Mm8WjBtZele1M8wEAOA/47yRb6qaXpN6rqvcBP0tydXur8xRgp/Z2pTpo7xRcPXl8pqanqr6XZMlc7c8JMxqGc2hmyt1Ec+V0LXM8eFeSRuBvgNOreaj008AdNOPQ1FE7KUeDe2gud2bPoyRJkjpzwowkSZI6s3iUJElSZxaPkiRJ6sziUZIkSZ1ZPEqSJKkzi0dJkiR1ZvEoSZKkziweJUmS1JnFoyRJkjqzeJQkSVJnFo+SJEnqzOJRkiRJnVk8SpIkqbOFow5AkrZXVbUY+DBwDLA7sAm4ATgjySNVtQr4apK/GNL+lrTb22sY25OkqdjzKEmzoKqeC9wCvAxYmuSFwKuA7wJ3VNUvjzI+SZopex4laXacA9yfZPn4giQ/Bs6uqieAV44sMkkagMWjJA1ZVe1Kc6v616b6fZKL2vWOnvS+W4GVSW6YsGwJE25FV9VrgfOAfYEdgOuBjwJ7A18Edq+qR4GbkxzVvmdP4ALgt4AFwO3AnyS5v/39Cppb6vcDHwfOTfL5AdMgaTvlbWtJGr6DgbuSPDTMjVbVAuAfgAuS7A0sAW4D3pDkLuB1wIYke00oHF8A3Ah8AdgL2BP4DHB9W5iOWwq8FXi9haOkrbHnUZKGb3dgwyxs9/nALjSTbkjyJLBqG+9ZAZyX5PoJy26sqhOAc4F3t8t2TPIHQ41W0nbJ4lGShu8RYOgTYpI8XlWfAL5YVTcANwPrkmzeytveAvxuVZ0/afkC4AcTXn9puNFK2l5ZPErS8H0V2KeqXp7k3mm8bzPNOMYtSrKyqj4JHAS8Hfh0VX0syY1beMsY8L4k100jDknaIsc8StKQJfk5cCFwUTtO8Vmq6j1VddoUb30MeNGkZe+YYvsbk3w5yenAicDZWwnnJuDoyQur6uB2Yo8kTYs9j5I0O84FrgGurqoTk3yvqnYGlgMfAg6leYD4SwGq6jk0k1/eX1VrgJ8A76XpYVzYrvMmYBnNbOgH2oeQv5Pm2ZEA/wc8t6qeT9Pj+CRwJrC+vd29sl12JM1YyNfNagYkbZfseZSkWZBkE/A7wFdoZjY/CqyneWj465PcA6wBjq2q7wP7AVcCdwP3Ad8B9gGOBR6uqkNoisv1wN9X1Q9oisa9gOPafT4CrAUeaPf7kiQ/Ag4EXgE82P4cTvPg8h/Pdh4kbX8WjI2NjToGSZIk9YQ9j5IkSerM4lGSJEmdWTxKkiSpM4tHSZIkdWbxKEmSpM4sHiVJktSZxaMkSZI6s3iUJElSZ/8PzqyZRtMYTjYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x242a806f908>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# CODE FOR CUSTOM GRAPHICS NOT INCLUDED"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even though the clustering looks good statistically, it always makes sense to do a visual check, which can be done using the `gmaps` module. `gmaps` uses the Google Maps Javascript API (Google provides a free API key) to generate a HTML Google Map. We plot the **unique** coordinates of flats, color-coded by label, to give us a map with 7 distinct zones. The map shows that 7 clusters is a good fit (or at least isn't terrible)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3efd234e0deb49749c78aba49891a3ee",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Figure(layout=FigureLayout(height='420px'))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# GMAPS\n",
    "GOOGLE_API = '[YOUR GOOGLE API KEY HERE]'\n",
    "gmaps.configure(api_key = GOOGLE_API)\n",
    "\n",
    "# Configure labels and column names\n",
    "dat_jw_plot = dat_jw[['lat', 'lon', 'label']].copy()\n",
    "dat_jw_plot['label'] = dat_jw_plot['label'] + 1\n",
    "dat_jw_plot = dat_jw_plot.rename(columns = {'lat': 'latitude', 'lon': 'longitude'})\n",
    "\n",
    "# Remove duplicates\n",
    "dat_jw_plot = dat_jw_plot.drop_duplicates()\n",
    "latlon_jw = dat_jw_plot[['latitude', 'longitude']]\n",
    "\n",
    "# Create layers\n",
    "# Colour code:\n",
    "#  - 1: Blue\n",
    "#  - 2: Orange\n",
    "#  - 3: Green\n",
    "#  - 4: Red\n",
    "#  - 5: Purple\n",
    "#  - 6: Brown\n",
    "#  - 7: Pink\n",
    "#  - 8: Grey\n",
    "c1 = gmaps.symbol_layer(latlon_jw[dat_jw_plot.label == 1], fill_color = '#1f77b4', stroke_color='#1f77b4', scale = 2)\n",
    "c2 = gmaps.symbol_layer(latlon_jw[dat_jw_plot.label == 2], fill_color = '#ff7f0e', stroke_color='#ff7f0e', scale = 2) \n",
    "c3 = gmaps.symbol_layer(latlon_jw[dat_jw_plot.label == 3], fill_color = '#2ca02c', stroke_color='#2ca02c', scale = 2)\n",
    "c4 = gmaps.symbol_layer(latlon_jw[dat_jw_plot.label == 4], fill_color = '#d62728', stroke_color='#d62728', scale = 2)\n",
    "c5 = gmaps.symbol_layer(latlon_jw[dat_jw_plot.label == 5], fill_color = '#9467bd', stroke_color='#9467bd', scale = 2)\n",
    "c6 = gmaps.symbol_layer(latlon_jw[dat_jw_plot.label == 6], fill_color = '#8c564b', stroke_color='#8c564b', scale = 2)\n",
    "c7 = gmaps.symbol_layer(latlon_jw[dat_jw_plot.label == 7], fill_color = '#e377c2', stroke_color='#e377c2', scale = 2)\n",
    "\n",
    "# Create base map\n",
    "t1 = gmaps.figure()\n",
    "\n",
    "# Add layers\n",
    "t1.add_layer(c1)\n",
    "t1.add_layer(c2)\n",
    "t1.add_layer(c3)\n",
    "t1.add_layer(c4)\n",
    "t1.add_layer(c5)\n",
    "t1.add_layer(c6)\n",
    "t1.add_layer(c7)\n",
    "\n",
    "# Visualise\n",
    "t1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Function to Plot Elbow Graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clust_town(town):\n",
    "    # Extract town data\n",
    "    temp_dat = hdb[['lat', 'lon']][hdb.town == town]\n",
    "    temp_dat = temp_dat.reset_index(drop = True)\n",
    "\n",
    "    # Normalise\n",
    "    temp_mm = MinMaxScaler()\n",
    "    temp_mm.fit(temp_dat)\n",
    "    temp_dat_scaled = pd.DataFrame(temp_mm.transform(temp_dat), columns = ['lat', 'lon'])\n",
    "\n",
    "    # Set up values of k\n",
    "    temp_k = np.arange(2, 23, 1)\n",
    "\n",
    "    # Initialise results\n",
    "    temp_results = []\n",
    "\n",
    "    # Loop through values of k\n",
    "    for k in temp_k:\n",
    "\n",
    "        # Set up kmeans\n",
    "        temp_km = KMeans(n_clusters = k, random_state = 123)\n",
    "\n",
    "        # Fit data\n",
    "        temp_km.fit(temp_dat_scaled)\n",
    "\n",
    "        # Score data\n",
    "        temp_results.append(temp_km.inertia_)\n",
    "\n",
    "    # Plot\n",
    "    plt.plot(temp_k, temp_results)\n",
    "    plt.title(town)\n",
    "    plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOT RUN\n",
    "# Perform clustering on each town using k = 2 to 22 to choose an optimal k\n",
    "# for t in hdb.town.value_counts().index:\n",
    "    \n",
    "#     clust_town(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Function to Plot Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clust_map(town):\n",
    "    \n",
    "    # Extract town data\n",
    "    temp_dat = hdb[['lat', 'lon']][hdb.town == town]\n",
    "    temp_dat = temp_dat.reset_index(drop = True)\n",
    "    \n",
    "    # Normalise\n",
    "    temp_mm = MinMaxScaler()\n",
    "    temp_mm.fit(temp_dat)\n",
    "    temp_dat_scaled = pd.DataFrame(temp_mm.transform(temp_dat), columns = ['lat', 'lon'])\n",
    "    \n",
    "    # Get optimal clusters\n",
    "    opt_clust = disp_clust.loc[town][0]\n",
    "    \n",
    "    # Fitting 7 clusters:\n",
    "    temp_km = KMeans(n_clusters = opt_clust, random_state = 123)\n",
    "    temp_km.fit(temp_dat_scaled)\n",
    "    \n",
    "    # Configure labels and column names\n",
    "    plot_labels = temp_dat[['lat', 'lon']].copy()\n",
    "    plot_labels['label'] = temp_km.labels_ + 1\n",
    "    plot_labels = plot_labels.rename(columns = {'lat': 'latitude', 'lon': 'longitude'})\n",
    "\n",
    "    # Remove duplicates\n",
    "    plot_labels = plot_labels.drop_duplicates()\n",
    "    plotdata = plot_labels[['latitude', 'longitude']]\n",
    "    \n",
    "    # Configure colours\n",
    "    temp_colors = sns.color_palette().as_hex()\n",
    "    \n",
    "    # Create base graph\n",
    "    out_graph = gmaps.figure()\n",
    "    \n",
    "    # Add layers sequentially\n",
    "    for i in range(opt_clust):\n",
    "        \n",
    "        # Add layer\n",
    "        out_graph.add_layer(\n",
    "            gmaps.symbol_layer(\n",
    "                plotdata[plot_labels.label == i + 1],\n",
    "                fill_color = temp_colors[i], stroke_color = temp_colors[i],\n",
    "                scale = 2\n",
    "            )\n",
    "        )\n",
    "    \n",
    "    # Output\n",
    "    return out_graph\n",
    "\n",
    "# Function to plot all points\n",
    "def clust_map_all(town):\n",
    "    \n",
    "    # Extract town data\n",
    "    temp_dat = hdb[['lat', 'lon']][hdb.town == town]\n",
    "    temp_dat = temp_dat.reset_index(drop = True)\n",
    "    \n",
    "    # Remove duplicates\n",
    "    temp_dat = temp_dat.drop_duplicates()\n",
    "    \n",
    "    # Create base graph\n",
    "    out_graph = gmaps.figure()\n",
    "    \n",
    "    # Add layer\n",
    "    out_graph.add_layer(gmaps.symbol_layer(temp_dat, scale = 2))\n",
    "    \n",
    "    # Output\n",
    "    return out_graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOT RUN\n",
    "# Run clust_map for all towns to check that the k values make sense\n",
    "# clust_map('[TOWN HERE]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Decide on clusters\n",
    "clust_results = [7, 5, 6, 7, 7, 4, 5, 6, 5, 5, 3, 6, 5, 7, 7, 6, 6, 5, 5, 8, 5, 5, 4, 2, 1, 2]\n",
    "\n",
    "# Create dataframe\n",
    "disp_clust = pd.DataFrame(\n",
    "    [hdb.town.value_counts().index, clust_results], index = ['Town', 'Clusters']\n",
    ").T.set_index('Town')\n",
    "\n",
    "# Create markdown table\n",
    "# print(tabulate.tabulate(disp_clust, tablefmt=\"pipe\", headers = 'keys'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate Optimal Clusters for All Towns\n",
    "Performing the same process of selecting an optimal *k* using the elbow and graphical methods, we obtained the following results:  \n",
    "  \n",
    "| Town            |   Clusters |\n",
    "|:----------------|-----------:|\n",
    "| JURONG WEST     |          7 |\n",
    "| SENGKANG        |          5 |\n",
    "| WOODLANDS       |          6 |\n",
    "| TAMPINES        |          7 |\n",
    "| BEDOK           |          7 |\n",
    "| YISHUN          |          4 |\n",
    "| PUNGGOL         |          5 |\n",
    "| HOUGANG         |          6 |\n",
    "| ANG MO KIO      |          5 |\n",
    "| CHOA CHU KANG   |          5 |\n",
    "| BUKIT BATOK     |          3 |\n",
    "| BUKIT MERAH     |          6 |\n",
    "| BUKIT PANJANG   |          5 |\n",
    "| TOA PAYOH       |          7 |\n",
    "| KALLANG/WHAMPOA |          7 |\n",
    "| PASIR RIS       |          6 |\n",
    "| SEMBAWANG       |          6 |\n",
    "| GEYLANG         |          5 |\n",
    "| QUEENSTOWN      |          5 |\n",
    "| CLEMENTI        |          8 |\n",
    "| JURONG EAST     |          5 |\n",
    "| SERANGOON       |          5 |\n",
    "| BISHAN          |          4 |\n",
    "| CENTRAL AREA    |          2 |\n",
    "| MARINE PARADE   |          1 |\n",
    "| BUKIT TIMAH     |          2 |\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Attach Clusters to Dataset\n",
    "To obtain the cluster labels for each town, we fit a K-Means model to the scaled data and the optimal *k* value for each town. This generates labels, which we then attach to the original dataset, town by town, under the feature `label`. To distinguish between the clusters of different towns, we generate a new feature `cluster` that appends the town name to the cluster label. That gives us a total of 134 clusters across Singapore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Get list of towns\n",
    "all_towns = hdb.town.value_counts().index\n",
    "\n",
    "# Initialise label\n",
    "hdb['label'] = 0\n",
    "\n",
    "# Loop through\n",
    "for town in all_towns:\n",
    "    \n",
    "    # Extract town data\n",
    "    temp_dat = hdb[['lat', 'lon']][hdb.town == town]\n",
    "    temp_dat = temp_dat.reset_index(drop = True)\n",
    "\n",
    "    # Normalise\n",
    "    temp_mm = MinMaxScaler()\n",
    "    temp_mm.fit(temp_dat)\n",
    "    temp_dat_scaled = pd.DataFrame(temp_mm.transform(temp_dat), columns = ['lat', 'lon'])\n",
    "    \n",
    "    # Get optimal clusters\n",
    "    opt_clust = disp_clust.loc[town][0]\n",
    "    \n",
    "    # Fit optimal clusters:\n",
    "    temp_km = KMeans(n_clusters = opt_clust, random_state = 123)\n",
    "    temp_km.fit(temp_dat_scaled)\n",
    "    \n",
    "    # Attach labels\n",
    "    hdb['label'][hdb.town == town] = temp_km.labels_ + 1\n",
    "\n",
    "# Attach town name to cluster label\n",
    "hdb['clust'] = hdb.town + '_' + hdb.label.astype('str')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the full map of all clusters in Singapore:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract town data\n",
    "full_dat = hdb[['lat', 'lon', 'clust', 'town']]\n",
    "full_dat = full_dat.reset_index(drop = True)\n",
    "\n",
    "# Remove duplicates\n",
    "full_dat = full_dat.drop_duplicates()\n",
    "\n",
    "# Rename columns\n",
    "full_dat = full_dat.rename(columns = {'lat': 'latitude', 'lon': 'longitude'})\n",
    "\n",
    "# Sort\n",
    "full_dat = full_dat.sort_values(['longitude', 'latitude'])\n",
    "\n",
    "# Extract towns\n",
    "all_towns = list(full_dat.town.unique())\n",
    "\n",
    "# Extract cluster names\n",
    "all_clust = list(full_dat.clust.unique())\n",
    "\n",
    "# Configure colours\n",
    "all_colors = sns.color_palette().as_hex()\n",
    "\n",
    "# Configure town colors\n",
    "town_colors = all_colors * 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create base graph\n",
    "out_graph = gmaps.figure()\n",
    "\n",
    "# Loop through clusters\n",
    "for t in range(len(all_towns)):\n",
    "    \n",
    "    # Get clusters\n",
    "    temp_clust = full_dat.clust[full_dat.town == all_towns[t]].unique()\n",
    "    \n",
    "    # Town coordinates\n",
    "    temp_towndata = full_dat[['latitude', 'longitude']][full_dat.town == all_towns[t]].copy()\n",
    "    \n",
    "    # Get town coords\n",
    "    temp_towncoords = np.array([[n, m] for n, m in zip(temp_towndata.longitude, temp_towndata.latitude)])\n",
    "    \n",
    "    # Calculate convex hull\n",
    "    temp_townhull = ConvexHull(temp_towncoords)\n",
    "    \n",
    "    # Drawing\n",
    "    temp_towndrawing = gmaps.drawing_layer(features=[\n",
    "         gmaps.Polygon(\n",
    "             [(n, m) for n, m in zip(temp_towncoords[temp_townhull.vertices, 1], temp_towncoords[temp_townhull.vertices, 0])],\n",
    "             fill_color = town_colors[t], fill_opacity = 0.3, stroke_color = 'black'\n",
    "         )\n",
    "    ], show_controls = False)\n",
    "\n",
    "    # Add drawing\n",
    "    out_graph.add_layer(temp_towndrawing)\n",
    "    \n",
    "    # Get town center\n",
    "    \n",
    "    # Add point\n",
    "    out_graph.add_layer(\n",
    "        gmaps.symbol_layer(\n",
    "            [(temp_towndata.latitude.median(), temp_towndata.longitude.median())],\n",
    "            fill_color = town_colors[t], stroke_color = town_colors[t],\n",
    "            scale = 1\n",
    "        )\n",
    "    )\n",
    "    \n",
    "    for c in range(len(temp_clust)):\n",
    "\n",
    "        # Extract coordinates\n",
    "        temp_plotdata = full_dat[['latitude', 'longitude']][full_dat.clust == temp_clust[c]].copy()\n",
    "\n",
    "        # Get coords\n",
    "        temp_coords = np.array([[x, y] for x, y in zip(temp_plotdata.longitude, temp_plotdata.latitude)])\n",
    "\n",
    "        # Calculate convex hull\n",
    "        temp_hull = ConvexHull(temp_coords)\n",
    "\n",
    "        # Drawing\n",
    "        temp_drawing = gmaps.drawing_layer(features=[\n",
    "             gmaps.Polygon(\n",
    "                 [(x, y) for x, y in zip(temp_coords[temp_hull.vertices, 1], temp_coords[temp_hull.vertices, 0])],\n",
    "                 fill_color = all_colors[c], fill_opacity = 0.3, stroke_color = all_colors[c]\n",
    "             )\n",
    "        ], show_controls = False)\n",
    "\n",
    "        # Add drawing\n",
    "        out_graph.add_layer(temp_drawing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e43f243e39b04161882786fe0edee5d2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Figure(layout=FigureLayout(height='420px'))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "out_graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion\n",
    "In this post, I demonstrated how an address created from block numbers and streets in the HDB resale flat dataset could be used to generate new features. Geocoding was used to convert addresses into geographic coordinates, and coordinates were used to generate clusters within each town. This produced a total of 134 clusters across Singapore. Hopefully, these will be useful when we develop our machine learning model to predict resale flat prices."
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}