Code.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "52d56c18-bdaf-44d0-a669-fc1ba5f41091",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setting up the environment\n",
    "import folium\n",
    "import json\n",
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "from matplotlib import pyplot as plt\n",
    "import timeit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "636d4126-6ffd-48f0-8389-21acfe3bb48f",
   "metadata": {},
   "source": [
    "# Part 0: Pre-process the datasets\n",
    "In this part, we will pre-process the datasets, such as crime, commerical land use, population, street. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "390b24c7-082c-43e5-bda8-7b715883ec75",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7533119\n",
      "THEFT                                1585804\n",
      "BATTERY                              1381696\n",
      "CRIMINAL DAMAGE                       858601\n",
      "NARCOTICS                             742126\n",
      "ASSAULT                               485919\n",
      "OTHER OFFENSE                         468856\n",
      "BURGLARY                              416782\n",
      "MOTOR VEHICLE THEFT                   349047\n",
      "DECEPTIVE PRACTICE                    327095\n",
      "ROBBERY                               283143\n",
      "CRIMINAL TRESPASS                     209853\n",
      "WEAPONS VIOLATION                      97558\n",
      "PROSTITUTION                           69560\n",
      "OFFENSE INVOLVING CHILDREN             53505\n",
      "PUBLIC PEACE VIOLATION                 51559\n",
      "SEX OFFENSE                            29337\n",
      "CRIM SEXUAL ASSAULT                    27754\n",
      "INTERFERENCE WITH PUBLIC OFFICER       17923\n",
      "LIQUOR LAW VIOLATION                   14679\n",
      "GAMBLING                               14608\n",
      "ARSON                                  12806\n",
      "HOMICIDE                               11841\n",
      "KIDNAPPING                              7090\n",
      "CRIMINAL SEXUAL ASSAULT                 4880\n",
      "INTIMIDATION                            4452\n",
      "STALKING                                4311\n",
      "CONCEALED CARRY LICENSE VIOLATION        899\n",
      "OBSCENITY                                767\n",
      "PUBLIC INDECENCY                         188\n",
      "NON-CRIMINAL                             178\n",
      "OTHER NARCOTIC VIOLATION                 141\n",
      "HUMAN TRAFFICKING                         89\n",
      "NON - CRIMINAL                            38\n",
      "RITUALISM                                 24\n",
      "NON-CRIMINAL (SUBJECT SPECIFIED)           9\n",
      "DOMESTIC VIOLENCE                          1\n",
      "Name: Primary Type, dtype: int64\n"
     ]
    }
   ],
   "source": [
    "# Crime dataset\n",
    "# Read csv file\n",
    "data= pd.read_csv(\"Crimes_-_2001_to_Present.csv\")\n",
    "\n",
    "# Number of primary crime\n",
    "print(data['Primary Type'].count())\n",
    "\n",
    "# Catalog of primary crime\n",
    "print(data['Primary Type'].value_counts())\n",
    "\n",
    "# Select rows based on the condition\n",
    "options = ['THEFT', 'BATTERY', 'ASSAULT', 'ROBBERY'] \n",
    "rslt_df = data[data['Primary Type'].isin(options)] \n",
    "\n",
    "# Select columns from the rows\n",
    "new = rslt_df[['Primary Type', 'Latitude', 'Longitude', 'Location']].copy()\n",
    "\n",
    "# Remove rows with NA values\n",
    "new = new.dropna()\n",
    "\n",
    "# Save the new dataset\n",
    "new.to_csv('new_crimedata_Chicago.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ffd049ab-11c2-4d2a-9ede-7ce89478c667",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>id</th>\n",
       "      <th>geometry</th>\n",
       "      <th>addr:street</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>relation/2814726</td>\n",
       "      <td>POINT (-87.69410 41.92517)</td>\n",
       "      <td>North Talman Avenue</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>relation/2817831</td>\n",
       "      <td>POINT (-87.66430 41.94689)</td>\n",
       "      <td>North Southport Avenue</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>relation/9765115</td>\n",
       "      <td>POINT (-87.69038 42.00202)</td>\n",
       "      <td>North Western Avenue</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>way/28100638</td>\n",
       "      <td>POINT (-87.66769 41.92281)</td>\n",
       "      <td>North Ashland Avenue</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>way/28293211</td>\n",
       "      <td>POINT (-87.63548 41.88853)</td>\n",
       "      <td>Merchandise Mart Plaza</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 id                    geometry             addr:street\n",
       "1  relation/2814726  POINT (-87.69410 41.92517)     North Talman Avenue\n",
       "2  relation/2817831  POINT (-87.66430 41.94689)  North Southport Avenue\n",
       "3  relation/9765115  POINT (-87.69038 42.00202)    North Western Avenue\n",
       "4      way/28100638  POINT (-87.66769 41.92281)    North Ashland Avenue\n",
       "5      way/28293211  POINT (-87.63548 41.88853)  Merchandise Mart Plaza"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Commerical land use dataset\n",
    "# Read geojson\n",
    "gdf = gpd.read_file('osm_shops_restaurants.geojson')\n",
    "\n",
    "# Keep only certain columns\n",
    "gdf = gdf[['id', 'geometry', 'addr:street']]\n",
    "\n",
    "# Remove rows with NA values\n",
    "gdf = gdf.dropna()\n",
    "\n",
    "# Show first rows in gdf\n",
    "gdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "08b549bc-1697-418a-9a6e-204a8317bc2e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\eugen\\anaconda3\\envs\\myenv\\lib\\site-packages\\geopandas\\io\\file.py:362: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.\n",
      "  pd.Int64Index,\n"
     ]
    }
   ],
   "source": [
    "# Street dataset\n",
    "# Read geojson\n",
    "street = gpd.read_file('Street_Copy.geojson')\n",
    "\n",
    "# Keep only certain columns\n",
    "street = street[['geometry','street_nam','street_typ','class']]\n",
    "\n",
    "# Save the new dataset\n",
    "street.to_file(\"street_clean.geojson\",driver=\"GeoJSON\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a72ad461-fc97-4c90-90ee-d634cf021619",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\eugen\\anaconda3\\envs\\myenv\\lib\\site-packages\\geopandas\\io\\file.py:362: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.\n",
      "  pd.Int64Index,\n"
     ]
    }
   ],
   "source": [
    "# Boundaries dataset\n",
    "# Read geojson\n",
    "boundaries = gpd.read_file('boundaries.geojson')\n",
    "\n",
    "# Reorganize the dataset\n",
    "boundaries = boundaries[['zip','geometry','shape_area']]\n",
    "\n",
    "# Save the new dataset\n",
    "boundaries.to_file(\"boundaries_new.geojson\",driver=\"GeoJSON\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "03723a2c-8c82-4fd6-8a28-a934d4721627",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Demographic dataset (merging them into boundaries)\n",
    "# Load the boundaries dataset\n",
    "with open('boundaries_new.geojson', encoding=\"utf-8\") as json_file:\n",
    "    data = json.load(json_file)\n",
    "data = data['features']\n",
    "data_zip = pd.DataFrame.from_dict(pd.json_normalize(data), orient='columns')\n",
    "\n",
    "# Extract useful information (e.g. zip code, area, coordinates) from the dataset\n",
    "data_zip = data_zip.drop(['type','geometry.type'],axis =1)\n",
    "data_zip = data_zip.rename(columns={'properties.zip':'zip','properties.shape_area':'area','geometry.coordinates':'coordinates'})\n",
    "data_zip.dropna(inplace=True)\n",
    "data_zip.reset_index(inplace=True)\n",
    "data_zip = data_zip.drop(['index'],axis =1)\n",
    "data_zip['zip'] = data_zip['zip'].astype(int)\n",
    "data_zip['area'] = data_zip['area'].astype(float)\n",
    "data_copy = data_zip\n",
    "data_copy\n",
    "\n",
    "# Load the demographic dataset and extract useful information\n",
    "demographic = pd.read_csv(\"Chicago_Population_Counts.csv\")\n",
    "demographic = demographic[['Geography','Population - Total']]\n",
    "demographic = demographic.drop(range(0,61),axis=0)\n",
    "demographic = demographic.rename(columns={'Geography':'zip','Population - Total':'population'})\n",
    "demographic.dropna(inplace=True)\n",
    "demographic.reset_index(inplace=True)\n",
    "demographic = demographic.drop(['index'],axis =1)\n",
    "demographic['zip'] = demographic['zip'].astype(int)\n",
    "demographic['population'] = demographic['population'].astype(int)\n",
    "demographic\n",
    "\n",
    "# Merge and save the new dataset\n",
    "data_demographic = pd.merge(data_copy,demographic,on='zip')\n",
    "data_demographic.to_json('demographic.json')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a7405c6-05de-4f4b-b330-30594f205b5f",
   "metadata": {},
   "source": [
    "# Part 1: Calculate interanction points of streets in Chicago\n",
    "In this part, we will calculate the new variable 'street weight' that combines the road density and road class, which could be used to represent the condition of transportation. We mainly use the 'segment interaction' and 'spatial index'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a69cfc3a-bb3e-4ce2-8706-9cf8e1e99d69",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>name</th>\n",
       "      <th>type</th>\n",
       "      <th>class</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>YALE</td>\n",
       "      <td>AVE</td>\n",
       "      <td>4</td>\n",
       "      <td>[[[-87.63187018233509, 41.78080974758186], [-8...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>COTTAGE GROVE</td>\n",
       "      <td>AVE</td>\n",
       "      <td>2</td>\n",
       "      <td>[[[-87.60548932047928, 41.76222964352745], [-8...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>CAMPBELL</td>\n",
       "      <td>AVE</td>\n",
       "      <td>4</td>\n",
       "      <td>[[[-87.6838184819274, 41.69724384442973], [-87...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>SANGAMON</td>\n",
       "      <td>ST</td>\n",
       "      <td>4</td>\n",
       "      <td>[[[-87.64584189298007, 41.69339840044404], [-8...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>118TH</td>\n",
       "      <td>ST</td>\n",
       "      <td>4</td>\n",
       "      <td>[[[-87.67053671319235, 41.67975927489014], [-8...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55193</th>\n",
       "      <td>BRAINARD</td>\n",
       "      <td>AVE</td>\n",
       "      <td>3</td>\n",
       "      <td>[[[-87.55784956501411, 41.65791820995521], [-8...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55194</th>\n",
       "      <td>MILWAUKEE</td>\n",
       "      <td>AVE</td>\n",
       "      <td>3</td>\n",
       "      <td>[[[-87.73902939370839, 41.94834404748583], [-8...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55195</th>\n",
       "      <td>KENNETH</td>\n",
       "      <td>AVE</td>\n",
       "      <td>4</td>\n",
       "      <td>[[[-87.73879459600205, 41.94718195189966], [-8...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55196</th>\n",
       "      <td>OHIO</td>\n",
       "      <td>ST</td>\n",
       "      <td>4</td>\n",
       "      <td>[[[-87.74339277679763, 41.891415000932874], [-...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55197</th>\n",
       "      <td>LEAVITT</td>\n",
       "      <td>ST</td>\n",
       "      <td>E</td>\n",
       "      <td>[[[-87.68107407701358, 41.8639693463979], [-87...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>55198 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                name type class  \\\n",
       "0               YALE  AVE     4   \n",
       "1      COTTAGE GROVE  AVE     2   \n",
       "2           CAMPBELL  AVE     4   \n",
       "3           SANGAMON   ST     4   \n",
       "4              118TH   ST     4   \n",
       "...              ...  ...   ...   \n",
       "55193       BRAINARD  AVE     3   \n",
       "55194      MILWAUKEE  AVE     3   \n",
       "55195        KENNETH  AVE     4   \n",
       "55196           OHIO   ST     4   \n",
       "55197        LEAVITT   ST     E   \n",
       "\n",
       "                                                geometry  \n",
       "0      [[[-87.63187018233509, 41.78080974758186], [-8...  \n",
       "1      [[[-87.60548932047928, 41.76222964352745], [-8...  \n",
       "2      [[[-87.6838184819274, 41.69724384442973], [-87...  \n",
       "3      [[[-87.64584189298007, 41.69339840044404], [-8...  \n",
       "4      [[[-87.67053671319235, 41.67975927489014], [-8...  \n",
       "...                                                  ...  \n",
       "55193  [[[-87.55784956501411, 41.65791820995521], [-8...  \n",
       "55194  [[[-87.73902939370839, 41.94834404748583], [-8...  \n",
       "55195  [[[-87.73879459600205, 41.94718195189966], [-8...  \n",
       "55196  [[[-87.74339277679763, 41.891415000932874], [-...  \n",
       "55197  [[[-87.68107407701358, 41.8639693463979], [-87...  \n",
       "\n",
       "[55198 rows x 4 columns]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Based on our demands, we will read and clean the dataset that we have pre-processed.\n",
    "# Read the geojson and extract the properties that we need.\n",
    "with open('street_clean.geojson', encoding=\"utf-8\") as json_file:\n",
    "    data = json.load(json_file)\n",
    "data = data['features']\n",
    "data_street = pd.DataFrame.from_dict(pd.json_normalize(data), orient='columns')\n",
    "data_street = data_street.drop(['type','geometry.type','geometry'],axis =1)\n",
    "data_street = data_street.rename(columns={'properties.street_nam':'name','properties.street_typ':'type',\n",
    "                                          'properties.class':'class','geometry.coordinates':'geometry'})\n",
    "\n",
    "# Check and rebuild the pre-processed street dataset\n",
    "data_street.dropna(inplace=True)\n",
    "data_street.reset_index(inplace=True)\n",
    "data_street = data_street.drop(['index'],axis =1)\n",
    "data_street['name'] = data_street['name'].astype(str)\n",
    "data_street['type'] = data_street['type'].astype(str)\n",
    "data_street['class'] = data_street['class'].astype(str)\n",
    "data_copy = data_street\n",
    "data_copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d189c811-8a62-446b-a492-7c0c6f7f0dca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0,\n",
       " 1,\n",
       " 2,\n",
       " 3,\n",
       " 4,\n",
       " 5,\n",
       " 6,\n",
       " 7,\n",
       " 8,\n",
       " 9,\n",
       " 10,\n",
       " 11,\n",
       " 12,\n",
       " 13,\n",
       " 14,\n",
       " 15,\n",
       " 16,\n",
       " 17,\n",
       " 18,\n",
       " 19,\n",
       " 20,\n",
       " 21,\n",
       " 22,\n",
       " 23,\n",
       " 24,\n",
       " 25,\n",
       " 26,\n",
       " 27,\n",
       " 28,\n",
       " 29,\n",
       " 30,\n",
       " 31,\n",
       " 32,\n",
       " 33,\n",
       " 34,\n",
       " 36,\n",
       " 37,\n",
       " 38,\n",
       " 39,\n",
       " 41,\n",
       " 42,\n",
       " 43,\n",
       " 44,\n",
       " 45,\n",
       " 46,\n",
       " 48,\n",
       " 49,\n",
       " 50,\n",
       " 51,\n",
       " 52,\n",
       " 53,\n",
       " 54,\n",
       " 55,\n",
       " 56,\n",
       " 57,\n",
       " 58,\n",
       " 59,\n",
       " 61,\n",
       " 63,\n",
       " 65,\n",
       " 66,\n",
       " 67,\n",
       " 68,\n",
       " 69,\n",
       " 70,\n",
       " 71,\n",
       " 73,\n",
       " 74,\n",
       " 75,\n",
       " 76,\n",
       " 77,\n",
       " 78,\n",
       " 79,\n",
       " 80,\n",
       " 81,\n",
       " 82,\n",
       " 83,\n",
       " 84,\n",
       " 85,\n",
       " 86,\n",
       " 87,\n",
       " 88,\n",
       " 89,\n",
       " 90,\n",
       " 91,\n",
       " 99,\n",
       " 100,\n",
       " 102,\n",
       " 103,\n",
       " 104,\n",
       " 105,\n",
       " 106,\n",
       " 107,\n",
       " 108,\n",
       " 110,\n",
       " 111,\n",
       " 112,\n",
       " 115,\n",
       " 116,\n",
       " 118,\n",
       " 119,\n",
       " 121,\n",
       " 122,\n",
       " 124,\n",
       " 125,\n",
       " 127,\n",
       " 129,\n",
       " 130,\n",
       " 131,\n",
       " 132,\n",
       " 134,\n",
       " 135,\n",
       " 136,\n",
       " 138,\n",
       " 139,\n",
       " 140,\n",
       " 141,\n",
       " 142,\n",
       " 143,\n",
       " 144,\n",
       " 145,\n",
       " 146,\n",
       " 147,\n",
       " 148,\n",
       " 149,\n",
       " 150,\n",
       " 151,\n",
       " 152,\n",
       " 153,\n",
       " 154,\n",
       " 155,\n",
       " 156,\n",
       " 157,\n",
       " 158,\n",
       " 159,\n",
       " 160,\n",
       " 161,\n",
       " 162,\n",
       " 163,\n",
       " 164,\n",
       " 165,\n",
       " 166,\n",
       " 167,\n",
       " 169,\n",
       " 170,\n",
       " 171,\n",
       " 173,\n",
       " 174,\n",
       " 175,\n",
       " 176,\n",
       " 177,\n",
       " 178,\n",
       " 179,\n",
       " 180,\n",
       " 181,\n",
       " 182,\n",
       " 183,\n",
       " 186,\n",
       " 187,\n",
       " 188,\n",
       " 189,\n",
       " 190,\n",
       " 191,\n",
       " 192,\n",
       " 193,\n",
       " 194,\n",
       " 196,\n",
       " 197,\n",
       " 198,\n",
       " 200,\n",
       " 201,\n",
       " 202,\n",
       " 203,\n",
       " 204,\n",
       " 205,\n",
       " 206,\n",
       " 207,\n",
       " 211,\n",
       " 212,\n",
       " 213,\n",
       " 214,\n",
       " 215,\n",
       " 217,\n",
       " 218,\n",
       " 219,\n",
       " 221,\n",
       " 222,\n",
       " 223,\n",
       " 224,\n",
       " 225,\n",
       " 226,\n",
       " 227,\n",
       " 228,\n",
       " 229,\n",
       " 230,\n",
       " 231,\n",
       " 232,\n",
       " 233,\n",
       " 234,\n",
       " 235,\n",
       " 236,\n",
       " 237,\n",
       " 238,\n",
       " 239,\n",
       " 240,\n",
       " 241,\n",
       " 242,\n",
       " 243,\n",
       " 244,\n",
       " 245,\n",
       " 246,\n",
       " 247,\n",
       " 248,\n",
       " 249,\n",
       " 251,\n",
       " 252,\n",
       " 253,\n",
       " 254,\n",
       " 255,\n",
       " 258,\n",
       " 259,\n",
       " 261,\n",
       " 262,\n",
       " 263,\n",
       " 266,\n",
       " 268,\n",
       " 269,\n",
       " 270,\n",
       " 271,\n",
       " 272,\n",
       " 273,\n",
       " 274,\n",
       " 275,\n",
       " 276,\n",
       " 277,\n",
       " 278,\n",
       " 279,\n",
       " 280,\n",
       " 281,\n",
       " 282,\n",
       " 283,\n",
       " 284,\n",
       " 285,\n",
       " 287,\n",
       " 288,\n",
       " 289,\n",
       " 290,\n",
       " 291,\n",
       " 292,\n",
       " 293,\n",
       " 294,\n",
       " 296,\n",
       " 297,\n",
       " 298,\n",
       " 299,\n",
       " 300,\n",
       " 301,\n",
       " 302,\n",
       " 303,\n",
       " 304,\n",
       " 305,\n",
       " 306,\n",
       " 307,\n",
       " 308,\n",
       " 309,\n",
       " 311,\n",
       " 313,\n",
       " 314,\n",
       " 315,\n",
       " 316,\n",
       " 319,\n",
       " 320,\n",
       " 321,\n",
       " 322,\n",
       " 323,\n",
       " 324,\n",
       " 325,\n",
       " 326,\n",
       " 327,\n",
       " 328,\n",
       " 329,\n",
       " 330,\n",
       " 331,\n",
       " 332,\n",
       " 333,\n",
       " 334,\n",
       " 335,\n",
       " 336,\n",
       " 337,\n",
       " 338,\n",
       " 340,\n",
       " 341,\n",
       " 343,\n",
       " 345,\n",
       " 346,\n",
       " 347,\n",
       " 348,\n",
       " 349,\n",
       " 350,\n",
       " 351,\n",
       " 352,\n",
       " 353,\n",
       " 354,\n",
       " 355,\n",
       " 356,\n",
       " 357,\n",
       " 359,\n",
       " 360,\n",
       " 361,\n",
       " 362,\n",
       " 363,\n",
       " 364,\n",
       " 365,\n",
       " 367,\n",
       " 368,\n",
       " 369,\n",
       " 370,\n",
       " 371,\n",
       " 372,\n",
       " 373,\n",
       " 374,\n",
       " 375,\n",
       " 376,\n",
       " 377,\n",
       " 378,\n",
       " 379,\n",
       " 380,\n",
       " 382,\n",
       " 383,\n",
       " 384,\n",
       " 385,\n",
       " 386,\n",
       " 388,\n",
       " 389,\n",
       " 390,\n",
       " 391,\n",
       " 394,\n",
       " 395,\n",
       " 396,\n",
       " 397,\n",
       " 398,\n",
       " 399,\n",
       " 400,\n",
       " 401,\n",
       " 402,\n",
       " 403,\n",
       " 404,\n",
       " 405,\n",
       " 406,\n",
       " 407,\n",
       " 408,\n",
       " 409,\n",
       " 410,\n",
       " 411,\n",
       " 412,\n",
       " 413,\n",
       " 414,\n",
       " 415,\n",
       " 416,\n",
       " 417,\n",
       " 419,\n",
       " 420,\n",
       " 421,\n",
       " 423,\n",
       " 424,\n",
       " 425,\n",
       " 426,\n",
       " 428,\n",
       " 429,\n",
       " 431,\n",
       " 432,\n",
       " 433,\n",
       " 434,\n",
       " 435,\n",
       " 436,\n",
       " 437,\n",
       " 438,\n",
       " 439,\n",
       " 440,\n",
       " 441,\n",
       " 442,\n",
       " 443,\n",
       " 444,\n",
       " 445,\n",
       " 446,\n",
       " 447,\n",
       " 448,\n",
       " 449,\n",
       " 450,\n",
       " 451,\n",
       " 452,\n",
       " 453,\n",
       " 454,\n",
       " 455,\n",
       " 456,\n",
       " 457,\n",
       " 458,\n",
       " 459,\n",
       " 460,\n",
       " 461,\n",
       " 462,\n",
       " 463,\n",
       " 464,\n",
       " 465,\n",
       " 467,\n",
       " 468,\n",
       " 469,\n",
       " 470,\n",
       " 471,\n",
       " 472,\n",
       " 473,\n",
       " 474,\n",
       " 475,\n",
       " 476,\n",
       " 477,\n",
       " 479,\n",
       " 480,\n",
       " 481,\n",
       " 482,\n",
       " 483,\n",
       " 484,\n",
       " 485,\n",
       " 486,\n",
       " 487,\n",
       " 488,\n",
       " 489,\n",
       " 490,\n",
       " 491,\n",
       " 492,\n",
       " 493,\n",
       " 494,\n",
       " 496,\n",
       " 497,\n",
       " 498,\n",
       " 499,\n",
       " 500,\n",
       " 501,\n",
       " 503,\n",
       " 504,\n",
       " 505,\n",
       " 506,\n",
       " 507,\n",
       " 508,\n",
       " 510,\n",
       " 511,\n",
       " 513,\n",
       " 514,\n",
       " 515,\n",
       " 516,\n",
       " 517,\n",
       " 518,\n",
       " 519,\n",
       " 521,\n",
       " 522,\n",
       " 523,\n",
       " 524,\n",
       " 525,\n",
       " 526,\n",
       " 527,\n",
       " 529,\n",
       " 530,\n",
       " 531,\n",
       " 532,\n",
       " 533,\n",
       " 534,\n",
       " 535,\n",
       " 536,\n",
       " 537,\n",
       " 538,\n",
       " 539,\n",
       " 540,\n",
       " 541,\n",
       " 542,\n",
       " 543,\n",
       " 544,\n",
       " 545,\n",
       " 546,\n",
       " 547,\n",
       " 548,\n",
       " 549,\n",
       " 550,\n",
       " 551,\n",
       " 552,\n",
       " 553,\n",
       " 554,\n",
       " 555,\n",
       " 556,\n",
       " 557,\n",
       " 558,\n",
       " 559,\n",
       " 560,\n",
       " 561,\n",
       " 562,\n",
       " 563,\n",
       " 564,\n",
       " 565,\n",
       " 566,\n",
       " 567,\n",
       " 570,\n",
       " 571,\n",
       " 572,\n",
       " 573,\n",
       " 574,\n",
       " 575,\n",
       " 576,\n",
       " 577,\n",
       " 578,\n",
       " 582,\n",
       " 583,\n",
       " 584,\n",
       " 585,\n",
       " 586,\n",
       " 587,\n",
       " 588,\n",
       " 589,\n",
       " 590,\n",
       " 592,\n",
       " 593,\n",
       " 594,\n",
       " 595,\n",
       " 596,\n",
       " 597,\n",
       " 598,\n",
       " 599,\n",
       " 600,\n",
       " 601,\n",
       " 602,\n",
       " 603,\n",
       " 604,\n",
       " 605,\n",
       " 606,\n",
       " 607,\n",
       " 608,\n",
       " 609,\n",
       " 610,\n",
       " 612,\n",
       " 613,\n",
       " 614,\n",
       " 615,\n",
       " 616,\n",
       " 617,\n",
       " 618,\n",
       " 620,\n",
       " 621,\n",
       " 622,\n",
       " 623,\n",
       " 625,\n",
       " 627,\n",
       " 628,\n",
       " 630,\n",
       " 631,\n",
       " 632,\n",
       " 633,\n",
       " 634,\n",
       " 635,\n",
       " 636,\n",
       " 637,\n",
       " 638,\n",
       " 639,\n",
       " 641,\n",
       " 642,\n",
       " 643,\n",
       " 644,\n",
       " 645,\n",
       " 646,\n",
       " 647,\n",
       " 648,\n",
       " 649,\n",
       " 650,\n",
       " 651,\n",
       " 652,\n",
       " 653,\n",
       " 654,\n",
       " 655,\n",
       " 656,\n",
       " 658,\n",
       " 660,\n",
       " 661,\n",
       " 662,\n",
       " 663,\n",
       " 664,\n",
       " 666,\n",
       " 667,\n",
       " 668,\n",
       " 669,\n",
       " 670,\n",
       " 671,\n",
       " 672,\n",
       " 673,\n",
       " 674,\n",
       " 675,\n",
       " 676,\n",
       " 677,\n",
       " 678,\n",
       " 679,\n",
       " 681,\n",
       " 682,\n",
       " 683,\n",
       " 685,\n",
       " 686,\n",
       " 687,\n",
       " 689,\n",
       " 690,\n",
       " 691,\n",
       " 692,\n",
       " 693,\n",
       " 694,\n",
       " 695,\n",
       " 696,\n",
       " 697,\n",
       " 698,\n",
       " 699,\n",
       " 700,\n",
       " 701,\n",
       " 702,\n",
       " 703,\n",
       " 704,\n",
       " 706,\n",
       " 707,\n",
       " 708,\n",
       " 709,\n",
       " 710,\n",
       " 711,\n",
       " 712,\n",
       " 713,\n",
       " 715,\n",
       " 716,\n",
       " 717,\n",
       " 718,\n",
       " 720,\n",
       " 721,\n",
       " 722,\n",
       " 723,\n",
       " 724,\n",
       " 725,\n",
       " 726,\n",
       " 727,\n",
       " 728,\n",
       " 729,\n",
       " 730,\n",
       " 731,\n",
       " 733,\n",
       " 735,\n",
       " 736,\n",
       " 737,\n",
       " 739,\n",
       " 740,\n",
       " 741,\n",
       " 742,\n",
       " 743,\n",
       " 744,\n",
       " 745,\n",
       " 746,\n",
       " 747,\n",
       " 748,\n",
       " 749,\n",
       " 750,\n",
       " 751,\n",
       " 753,\n",
       " 754,\n",
       " 755,\n",
       " 756,\n",
       " 757,\n",
       " 758,\n",
       " 759,\n",
       " 760,\n",
       " 761,\n",
       " 762,\n",
       " 763,\n",
       " 764,\n",
       " 765,\n",
       " 766,\n",
       " 767,\n",
       " 768,\n",
       " 769,\n",
       " 770,\n",
       " 771,\n",
       " 772,\n",
       " 773,\n",
       " 774,\n",
       " 775,\n",
       " 776,\n",
       " 777,\n",
       " 778,\n",
       " 779,\n",
       " 780,\n",
       " 781,\n",
       " 782,\n",
       " 784,\n",
       " 785,\n",
       " 786,\n",
       " 787,\n",
       " 788,\n",
       " 789,\n",
       " 790,\n",
       " 791,\n",
       " 792,\n",
       " 793,\n",
       " 794,\n",
       " 795,\n",
       " 796,\n",
       " 797,\n",
       " 798,\n",
       " 799,\n",
       " 800,\n",
       " 801,\n",
       " 802,\n",
       " 803,\n",
       " 806,\n",
       " 807,\n",
       " 808,\n",
       " 809,\n",
       " 810,\n",
       " 811,\n",
       " 812,\n",
       " 813,\n",
       " 814,\n",
       " 815,\n",
       " 816,\n",
       " 817,\n",
       " 819,\n",
       " 820,\n",
       " 821,\n",
       " 822,\n",
       " 823,\n",
       " 824,\n",
       " 825,\n",
       " 826,\n",
       " 828,\n",
       " 829,\n",
       " 830,\n",
       " 831,\n",
       " 832,\n",
       " 833,\n",
       " 834,\n",
       " 835,\n",
       " 836,\n",
       " 838,\n",
       " 839,\n",
       " 840,\n",
       " 841,\n",
       " 842,\n",
       " 843,\n",
       " 844,\n",
       " 846,\n",
       " 847,\n",
       " 848,\n",
       " 849,\n",
       " 852,\n",
       " 854,\n",
       " 855,\n",
       " 856,\n",
       " 857,\n",
       " 858,\n",
       " 859,\n",
       " 860,\n",
       " 864,\n",
       " 865,\n",
       " 867,\n",
       " 869,\n",
       " 870,\n",
       " 871,\n",
       " 872,\n",
       " 873,\n",
       " 875,\n",
       " 876,\n",
       " 877,\n",
       " 878,\n",
       " 879,\n",
       " 880,\n",
       " 881,\n",
       " 882,\n",
       " 883,\n",
       " 884,\n",
       " 885,\n",
       " 886,\n",
       " 887,\n",
       " 888,\n",
       " 889,\n",
       " 890,\n",
       " 891,\n",
       " 892,\n",
       " 893,\n",
       " 894,\n",
       " 895,\n",
       " 897,\n",
       " 898,\n",
       " 899,\n",
       " 900,\n",
       " 901,\n",
       " 902,\n",
       " 903,\n",
       " 904,\n",
       " 905,\n",
       " 906,\n",
       " 907,\n",
       " 908,\n",
       " 909,\n",
       " 910,\n",
       " 911,\n",
       " 914,\n",
       " 915,\n",
       " 916,\n",
       " 917,\n",
       " 918,\n",
       " 919,\n",
       " 920,\n",
       " 921,\n",
       " 922,\n",
       " 923,\n",
       " 924,\n",
       " 925,\n",
       " 926,\n",
       " 928,\n",
       " 931,\n",
       " 932,\n",
       " 933,\n",
       " 934,\n",
       " 935,\n",
       " 936,\n",
       " 937,\n",
       " 938,\n",
       " 939,\n",
       " 940,\n",
       " 941,\n",
       " 942,\n",
       " 943,\n",
       " 944,\n",
       " 945,\n",
       " 946,\n",
       " 947,\n",
       " 948,\n",
       " 949,\n",
       " 951,\n",
       " 952,\n",
       " 953,\n",
       " 954,\n",
       " 955,\n",
       " 956,\n",
       " 957,\n",
       " 958,\n",
       " 959,\n",
       " 960,\n",
       " 961,\n",
       " 963,\n",
       " 964,\n",
       " 965,\n",
       " 966,\n",
       " 967,\n",
       " 968,\n",
       " 969,\n",
       " 970,\n",
       " 971,\n",
       " 972,\n",
       " 973,\n",
       " 974,\n",
       " 975,\n",
       " 976,\n",
       " 980,\n",
       " 981,\n",
       " 982,\n",
       " 983,\n",
       " 984,\n",
       " 985,\n",
       " 986,\n",
       " 987,\n",
       " 988,\n",
       " 989,\n",
       " 990,\n",
       " 991,\n",
       " 992,\n",
       " 993,\n",
       " 994,\n",
       " 995,\n",
       " 996,\n",
       " 997,\n",
       " 998,\n",
       " 1001,\n",
       " 1003,\n",
       " 1004,\n",
       " 1005,\n",
       " 1006,\n",
       " 1007,\n",
       " 1008,\n",
       " 1009,\n",
       " 1010,\n",
       " 1011,\n",
       " 1012,\n",
       " 1013,\n",
       " 1014,\n",
       " 1015,\n",
       " 1016,\n",
       " 1017,\n",
       " 1018,\n",
       " 1019,\n",
       " 1020,\n",
       " 1022,\n",
       " 1023,\n",
       " 1024,\n",
       " 1025,\n",
       " 1026,\n",
       " 1027,\n",
       " 1028,\n",
       " 1029,\n",
       " 1030,\n",
       " 1031,\n",
       " 1032,\n",
       " 1033,\n",
       " 1034,\n",
       " 1035,\n",
       " 1036,\n",
       " 1037,\n",
       " 1038,\n",
       " 1039,\n",
       " 1040,\n",
       " 1041,\n",
       " 1042,\n",
       " 1043,\n",
       " 1044,\n",
       " 1045,\n",
       " 1046,\n",
       " 1047,\n",
       " 1049,\n",
       " 1050,\n",
       " 1051,\n",
       " 1052,\n",
       " 1053,\n",
       " 1054,\n",
       " 1055,\n",
       " 1056,\n",
       " 1057,\n",
       " 1058,\n",
       " 1059,\n",
       " 1060,\n",
       " 1062,\n",
       " 1063,\n",
       " 1064,\n",
       " 1066,\n",
       " 1067,\n",
       " 1069,\n",
       " 1071,\n",
       " 1072,\n",
       " 1073,\n",
       " 1074,\n",
       " 1076,\n",
       " 1077,\n",
       " 1079,\n",
       " 1080,\n",
       " 1081,\n",
       " 1085,\n",
       " 1086,\n",
       " 1087,\n",
       " 1088,\n",
       " 1089,\n",
       " 1090,\n",
       " 1091,\n",
       " 1092,\n",
       " 1093,\n",
       " 1095,\n",
       " 1096,\n",
       " 1097,\n",
       " 1098,\n",
       " 1099,\n",
       " 1100,\n",
       " 1102,\n",
       " 1103,\n",
       " 1104,\n",
       " 1105,\n",
       " 1106,\n",
       " 1107,\n",
       " 1108,\n",
       " 1109,\n",
       " 1110,\n",
       " 1111,\n",
       " 1112,\n",
       " 1113,\n",
       " 1115,\n",
       " 1116,\n",
       " 1117,\n",
       " 1119,\n",
       " 1120,\n",
       " 1121,\n",
       " 1122,\n",
       " 1124,\n",
       " 1125,\n",
       " 1126,\n",
       " 1127,\n",
       " 1128,\n",
       " 1129,\n",
       " 1130,\n",
       " 1131,\n",
       " 1132,\n",
       " 1133,\n",
       " 1136,\n",
       " 1138,\n",
       " 1139,\n",
       " 1140,\n",
       " 1141,\n",
       " 1142,\n",
       " 1144,\n",
       " ...]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Due to the large area of Chicago, we decide to using bounding box to set up our own study area. \n",
    "# Range of study area (Boundaries of bounding box): x from -87.743 to -87.633, y from 41.736 to 41.816.\n",
    "# k is the indexes of records that we want to delete.\n",
    "k=[]\n",
    "for i in data_copy.index:\n",
    "    # Here, I delete the roads with types that we don't need, like the ramp of a expy/highway.\n",
    "    if data_copy.loc[i,'type'] == 'XR' or data_copy.loc[i,'type'] == 'TER' or data_copy.loc[i,'type'] == 'ER':\n",
    "        k.append(i)\n",
    "\n",
    "    # If any vertice of a segment is outside of the boundary of bounding box, I would delete the whole segment.\n",
    "    for j in range(0,len(data_copy.loc[i,'geometry'][0])):\n",
    "        \n",
    "        # Go through all coordinates of a record. Compare the coordinates with boundary.\n",
    "        if (data_copy.loc[i,'geometry'][0][j][0] <= -87.743 or data_copy.loc[i,'geometry'][0][j][0] >= -87.633 or data_copy.loc[i,'geometry'][0][j][1] <= 41.736 or data_copy.loc[i,'geometry'][0][j][1] >= 41.816):\n",
    "            k.append(i)\n",
    "\n",
    "# list of k have a lot of duplicate values.\n",
    "    # 1. Some records have the road type of XR, TER or ER, as well as have vertice(s) outside the bounding box.\n",
    "    # 2. Some records have more than 1 vertices outside the bounding box, so the index was saved multiple times.\n",
    "# So I have to deduplicate it.\n",
    "deduplicated_k = list()\n",
    "for item in k:\n",
    "    if item not in deduplicated_k:\n",
    "        deduplicated_k.append(item)\n",
    "deduplicated_k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "eccf6221-c55e-4b74-a8a4-0723c7f99e26",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop the records with index in list k.\n",
    "data_new=data_copy.drop(deduplicated_k)\n",
    "data_new.reset_index(inplace=True)\n",
    "data_new = data_new.drop(['index'],axis =1)\n",
    "\n",
    "# Still need some sentences to explain why we just make the class E or 9 as 5\n",
    "for i in data_new.index:\n",
    "    if data_new.loc[i,'class'] == 'E' or data_new.loc[i,'class'] == 9:\n",
    "        data_new.at[i,'class']=5\n",
    "data_new['class'] = data_new['class'].astype(int)\n",
    "data_new\n",
    "\n",
    "# Save the street dataset\n",
    "data_new.to_csv('street.csv')\n",
    "data_new.to_json('street.json')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c67bcee-c539-4dc3-aa5e-7eecff107eac",
   "metadata": {},
   "source": [
    "In the data, if a segment intersected with another, it would be cut. So, for the segment intersection, we just need to consider, if the first point or the last point of a segment is shared by at lest two of other segments. If yes, it would be the crossing point. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "645810ff-db05-4ea1-97d2-b108488b36e0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\eugen\\AppData\\Local\\Temp\\ipykernel_96452\\3097041537.py:28: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.\n",
      "  points = points.append(data,ignore_index = True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time:  1757.0008193999997\n"
     ]
    }
   ],
   "source": [
    "# Process the segment interaction\n",
    "start = timeit.default_timer()\n",
    "points = pd.DataFrame()\n",
    "\n",
    "for i in data_new.index:\n",
    "    # Only need to focus on starting and ending points here.\n",
    "    for j in (0,len(data_new.loc[i,'geometry'][0])-1):\n",
    "        # a is the point we use for comparison.\n",
    "        # y is used to cound the segments that sharing the point.\n",
    "        # z is the weight of a segment.\n",
    "        # As I suppose crimes happen more near high-class roads, I use 6-(class value) to represent the weight.\n",
    "        a = data_new.loc[i,'geometry'][0][j]\n",
    "        y = 1\n",
    "        z = 6-data_new.loc[i,'class']\n",
    "        \n",
    "        # Then I compare it with other points of segments with greater indexes.\n",
    "        for k in range(i+1,max(data_new.index)+1):\n",
    "            for l in (0,len(data_new.loc[k,'geometry'][0])-1):\n",
    "                # If find a point that have the same coordinate, I would make y add 1, means I find another segment sharing this vertice.\n",
    "                # I sum z up, which is the overall weight of points. If there are a lot of segments sharing the vertice, or the segments that sharing the vertice have larger weight, then the vertice would have a greater z.\n",
    "                if (data_new.loc[k,'geometry'][0][l]) == a:\n",
    "                    y = y+1\n",
    "                    z = z+6-data_new.loc[k,'class']\n",
    "        \n",
    "        # If y is equal or greater than 3, it would be considered as crossing points. So I add the point to point dataframe.\n",
    "        if y >= 3:\n",
    "            data = [[a,z]]\n",
    "            points = points.append(data,ignore_index = True)\n",
    "\n",
    "            \n",
    "stop = timeit.default_timer()\n",
    "print('Time: ', stop - start) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "fd288c56-6310-4745-a1a2-071bf2768323",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geometry</th>\n",
       "      <th>weight</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[-87.6540413652614, 41.76232762971525]</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[-87.65402050244748, 41.761416800680564]</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[-87.73178609194234, 41.744287454640435]</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>[-87.73295978195388, 41.744909524324726]</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>[-87.71298413347624, 41.771520084807925]</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5980</th>\n",
       "      <td>[-87.66195614621004, 41.755836444830884]</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5981</th>\n",
       "      <td>[-87.71311228988141, 41.79336771030571]</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5982</th>\n",
       "      <td>[-87.64263165885555, 41.78069141374427]</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5983</th>\n",
       "      <td>[-87.64376373665309, 41.776142196816636]</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5984</th>\n",
       "      <td>[-87.63919014217377, 41.74343406986939]</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5985 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                      geometry  weight\n",
       "0       [-87.6540413652614, 41.76232762971525]       8\n",
       "1     [-87.65402050244748, 41.761416800680564]      10\n",
       "2     [-87.73178609194234, 41.744287454640435]       6\n",
       "3     [-87.73295978195388, 41.744909524324726]       6\n",
       "4     [-87.71298413347624, 41.771520084807925]       8\n",
       "...                                        ...     ...\n",
       "5980  [-87.66195614621004, 41.755836444830884]       5\n",
       "5981   [-87.71311228988141, 41.79336771030571]      10\n",
       "5982   [-87.64263165885555, 41.78069141374427]       6\n",
       "5983  [-87.64376373665309, 41.776142196816636]       4\n",
       "5984   [-87.63919014217377, 41.74343406986939]       5\n",
       "\n",
       "[5985 rows x 2 columns]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Built the new interaction points with their geometry and weight (means the combination of the road density and road class).\n",
    "points = points.rename(columns={0:'geometry',1:'weight'})\n",
    "points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "e25a0013-909e-47c2-9519-efa1a852256e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geometry</th>\n",
       "      <th>weight</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[-87.6540413652614, 41.76232762971525]</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[-87.65402050244748, 41.761416800680564]</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[-87.73178609194234, 41.744287454640435]</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>[-87.73295978195388, 41.744909524324726]</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>[-87.71298413347624, 41.771520084807925]</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3639</th>\n",
       "      <td>[-87.64266584225966, 41.77612081106971]</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3640</th>\n",
       "      <td>[-87.67469165929282, 41.81033317369515]</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3641</th>\n",
       "      <td>[-87.66195614621004, 41.755836444830884]</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3642</th>\n",
       "      <td>[-87.71311228988141, 41.79336771030571]</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3643</th>\n",
       "      <td>[-87.63919014217377, 41.74343406986939]</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3644 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                      geometry  weight\n",
       "0       [-87.6540413652614, 41.76232762971525]       8\n",
       "1     [-87.65402050244748, 41.761416800680564]      10\n",
       "2     [-87.73178609194234, 41.744287454640435]       6\n",
       "3     [-87.73295978195388, 41.744909524324726]       6\n",
       "4     [-87.71298413347624, 41.771520084807925]       8\n",
       "...                                        ...     ...\n",
       "3639   [-87.64266584225966, 41.77612081106971]       6\n",
       "3640   [-87.67469165929282, 41.81033317369515]       7\n",
       "3641  [-87.66195614621004, 41.755836444830884]       7\n",
       "3642   [-87.71311228988141, 41.79336771030571]      10\n",
       "3643   [-87.63919014217377, 41.74343406986939]       5\n",
       "\n",
       "[3644 rows x 2 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We may have some duplicate points. For those vertices shared by more than 3 segments, it would be counted more than one time.\n",
    "k=[]\n",
    "for i in points.index:\n",
    "    x=points.loc[i,'geometry'][0]\n",
    "    y=points.loc[i,'geometry'][1]\n",
    "    for j in range(i+1,max(points.index)+1):\n",
    "        if (points.loc[j,'geometry'][0])==x and (points.loc[j,'geometry'][1])==y:\n",
    "            k.append(j)\n",
    "\n",
    "# I deduplicate in the same way as before.\n",
    "deduplicated_k = list()\n",
    "for item in k:\n",
    "    if item not in deduplicated_k:\n",
    "        deduplicated_k.append(item)\n",
    "deduplicated_k\n",
    "\n",
    "# Save the interaction points dataset\n",
    "points_new=points.drop(deduplicated_k)\n",
    "points_new.reset_index(inplace=True)\n",
    "points_new = points_new.drop(['index'],axis =1)\n",
    "points_new"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb98075f-5374-44e5-bdfc-d73c7b86ad4f",
   "metadata": {},
   "source": [
    "This step is simplification. If a point is too near to another, we have to delete the point. The algorithm to calculate distance with latitute and longitute comes from the first lecture of this course. The method is similar to the method of selecting segments inside of a bounding box. I use the k list as the indexes of records that I want to delete."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e413670b-0499-44e2-8c4e-93051224da6e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time:  456.6334698999999\n"
     ]
    }
   ],
   "source": [
    "# The method and process is same as before. \n",
    "start = timeit.default_timer()\n",
    "\n",
    "k=[]\n",
    "for i in points_new.index:\n",
    "    # I use if i not in k to avoid deleting the points that actually don't need to be deleted.\n",
    "    if i not in k:\n",
    "        for j in range(i+1,max(points_new.index)+1):\n",
    "            if j not in k:\n",
    "                d = 2 * 6371000 * np.arcsin(np.sqrt(np.sin((np.radians(points_new.loc[j,'geometry'][1]) - np.radians(points_new.loc[i,'geometry'][1]))/2)**2 + np.cos(np.radians(points_new.loc[i,'geometry'][1])) * np.cos(np.radians(points_new.loc[j,'geometry'][1])) * np.sin((np.radians(points_new.loc[j,'geometry'][0])-np.radians(points_new.loc[i,'geometry'][0]))/2)**2))\n",
    "                if d < 100:\n",
    "                    k.append(j)\n",
    "\n",
    "stop = timeit.default_timer()\n",
    "print('Time: ', stop - start)    \n",
    "\n",
    "deduplicated_k = list()\n",
    "for item in k:\n",
    "    if item not in deduplicated_k:\n",
    "        deduplicated_k.append(item)\n",
    "deduplicated_k\n",
    "\n",
    "# Sort out the final interaction points dataset\n",
    "points_final = points_new.drop(deduplicated_k)\n",
    "points_final.reset_index(inplace=True)\n",
    "points_final = points_final.drop(['index'],axis =1)\n",
    "points_final\n",
    "\n",
    "# Save interaction points dataset\n",
    "points_final.to_csv('points_final.csv')\n",
    "points_final.to_json('points_final.json')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48c281ea-df79-4173-bccb-cac399dfd88b",
   "metadata": {},
   "source": [
    "# Part 2: Calculate the number of crimes and commerical land use\n",
    "In this part, we will calculate the number of crimes and commerical land use.We mainly use the 'points in circles' and 'spatial index'."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7c977ac-aede-495b-b6cf-835f64fe7ff8",
   "metadata": {},
   "source": [
    "## 2.1 Crimes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "fad63004-1cae-483d-ba60-f12b70240192",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "544285\n"
     ]
    }
   ],
   "source": [
    "# Read the crime dataset\n",
    "data= pd.read_csv(\"new_crimedata_Chicago.csv\")\n",
    "crime_df = pd.DataFrame(data)\n",
    "\n",
    "# Built the empty column\n",
    "crime_df[\"Distance\"] = 0\n",
    "\n",
    "# Select the dataset using the bounding box (x from -87.743 to -87.633, y from 41.736 to 41.816)\n",
    "crime_selected = data[(data['Longitude'] <= -87.43) &\n",
    "          (data['Longitude'] >= -87.633) &\n",
    "          (data['Latitude'] <= 41.816) &\n",
    "          (data['Latitude'] >= 41.736)]\n",
    "\n",
    "# Print the selected dataset\n",
    "print(len(crime_selected))\n",
    "\n",
    "# Save the selected crime\n",
    "crime_selected.to_csv('crime_selected.csv')\n",
    "crime_selected.to_json('crime_selected.json')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "ce2e06c7-d368-430f-8afb-758494af5919",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BATTERY    228707\n",
      "THEFT      183179\n",
      "ASSAULT     80988\n",
      "ROBBERY     51411\n",
      "Name: Primary Type, dtype: int64\n",
      "-87.575997062\n",
      "41.766247597\n"
     ]
    }
   ],
   "source": [
    "# Read final interaction points\n",
    "data= pd.read_csv(\"crime_selected.csv\")\n",
    "crime_selected_df = pd.DataFrame(data)\n",
    "\n",
    "# Catalog of primary crime\n",
    "print(crime_selected_df['Primary Type'].value_counts())\n",
    "\n",
    "# Show the records of crime\n",
    "crime_selected_df[:5]\n",
    "\n",
    "# Built the geodataframe type of crime_selected_df\n",
    "crime_selected_df1 = gpd.GeoDataFrame(\n",
    "    crime_selected_df, geometry=gpd.points_from_xy(crime_selected_df['Longitude'], crime_selected_df['Latitude']))\n",
    "\n",
    "# Quick check one point in intersection_points_df1\n",
    "print(crime_selected_df1.loc[1,'geometry'].x)\n",
    "print(crime_selected_df1.loc[1,'geometry'].y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9c1fc9f0-0aea-40c0-9134-c57123ca833c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ID</th>\n",
       "      <th>geometry</th>\n",
       "      <th>weight</th>\n",
       "      <th>Latitude</th>\n",
       "      <th>Longitude</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>[-87.6540413652614, 41.76232762971525]</td>\n",
       "      <td>8</td>\n",
       "      <td>[-87.6540413652614</td>\n",
       "      <td>41.76232762971525]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>[-87.65402050244748, 41.761416800680564]</td>\n",
       "      <td>10</td>\n",
       "      <td>[-87.65402050244748</td>\n",
       "      <td>41.761416800680564]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>[-87.73178609194234, 41.744287454640435]</td>\n",
       "      <td>6</td>\n",
       "      <td>[-87.73178609194234</td>\n",
       "      <td>41.744287454640435]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>[-87.73295978195388, 41.744909524324726]</td>\n",
       "      <td>6</td>\n",
       "      <td>[-87.73295978195388</td>\n",
       "      <td>41.744909524324726]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>[-87.71298413347624, 41.771520084807925]</td>\n",
       "      <td>8</td>\n",
       "      <td>[-87.71298413347624</td>\n",
       "      <td>41.771520084807925]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3022</th>\n",
       "      <td>3022</td>\n",
       "      <td>[-87.65595139607863, 41.78868826260606]</td>\n",
       "      <td>8</td>\n",
       "      <td>[-87.65595139607863</td>\n",
       "      <td>41.78868826260606]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3023</th>\n",
       "      <td>3023</td>\n",
       "      <td>[-87.64266584225966, 41.77612081106971]</td>\n",
       "      <td>6</td>\n",
       "      <td>[-87.64266584225966</td>\n",
       "      <td>41.77612081106971]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3024</th>\n",
       "      <td>3024</td>\n",
       "      <td>[-87.66195614621004, 41.755836444830884]</td>\n",
       "      <td>7</td>\n",
       "      <td>[-87.66195614621004</td>\n",
       "      <td>41.755836444830884]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3025</th>\n",
       "      <td>3025</td>\n",
       "      <td>[-87.71311228988141, 41.79336771030571]</td>\n",
       "      <td>10</td>\n",
       "      <td>[-87.71311228988141</td>\n",
       "      <td>41.79336771030571]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3026</th>\n",
       "      <td>3026</td>\n",
       "      <td>[-87.63919014217377, 41.74343406986939]</td>\n",
       "      <td>5</td>\n",
       "      <td>[-87.63919014217377</td>\n",
       "      <td>41.74343406986939]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3027 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        ID                                  geometry  weight  \\\n",
       "0        0    [-87.6540413652614, 41.76232762971525]       8   \n",
       "1        1  [-87.65402050244748, 41.761416800680564]      10   \n",
       "2        2  [-87.73178609194234, 41.744287454640435]       6   \n",
       "3        3  [-87.73295978195388, 41.744909524324726]       6   \n",
       "4        4  [-87.71298413347624, 41.771520084807925]       8   \n",
       "...    ...                                       ...     ...   \n",
       "3022  3022   [-87.65595139607863, 41.78868826260606]       8   \n",
       "3023  3023   [-87.64266584225966, 41.77612081106971]       6   \n",
       "3024  3024  [-87.66195614621004, 41.755836444830884]       7   \n",
       "3025  3025   [-87.71311228988141, 41.79336771030571]      10   \n",
       "3026  3026   [-87.63919014217377, 41.74343406986939]       5   \n",
       "\n",
       "                 Latitude             Longitude  \n",
       "0      [-87.6540413652614    41.76232762971525]  \n",
       "1     [-87.65402050244748   41.761416800680564]  \n",
       "2     [-87.73178609194234   41.744287454640435]  \n",
       "3     [-87.73295978195388   41.744909524324726]  \n",
       "4     [-87.71298413347624   41.771520084807925]  \n",
       "...                   ...                   ...  \n",
       "3022  [-87.65595139607863    41.78868826260606]  \n",
       "3023  [-87.64266584225966    41.77612081106971]  \n",
       "3024  [-87.66195614621004   41.755836444830884]  \n",
       "3025  [-87.71311228988141    41.79336771030571]  \n",
       "3026  [-87.63919014217377    41.74343406986939]  \n",
       "\n",
       "[3027 rows x 5 columns]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Read final interaction points\n",
    "points = pd.read_csv('points_final.csv')\n",
    "intersection_points = pd.DataFrame(points)\n",
    "intersection_points.rename(columns={'Unnamed: 0':'ID'}, inplace=True)\n",
    "\n",
    "# Create two lists for the loop results to be placed\n",
    "lat = []\n",
    "lon = []\n",
    "\n",
    "# For each row in a varible,\n",
    "for row in intersection_points['geometry']:\n",
    "    row = \"\".join(row)\n",
    "    # Try to,\n",
    "    try:\n",
    "        # Split the row by comma and append\n",
    "        # everything before the comma to lat\n",
    "        lat.append(row.split(',')[0])\n",
    "        # Split the row by comma and append\n",
    "        # everything after the comma to lon\n",
    "        lon.append(row.split(',')[1])\n",
    "    # But if you get an error\n",
    "    except:\n",
    "        # append a missing value to lat\n",
    "        lat.append(np.NaN)\n",
    "        # append a missing value to lon\n",
    "        lon.append(np.NaN)\n",
    "\n",
    "# Create two new columns from lat and lon\n",
    "intersection_points['Latitude'] = lat\n",
    "intersection_points['Longitude'] = lon\n",
    "\n",
    "intersection_points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ff658d27-df69-4847-9ff0-ec096cf234d7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "41.761416800680564\n",
      "-87.65402050244748\n"
     ]
    }
   ],
   "source": [
    "# Clean the two new columns\n",
    "intersection_points['Latitude'] = [x.replace('[', '') for x in intersection_points['Latitude']]\n",
    "intersection_points['Longitude'] = [x.replace(']', '') for x in intersection_points['Longitude']]\n",
    "\n",
    "# Built the geodataframe type of intersection_points\n",
    "intersection_points_df1 = gpd.GeoDataFrame(\n",
    "    intersection_points, geometry=gpd.points_from_xy(intersection_points['Longitude'], intersection_points['Latitude']))\n",
    "\n",
    "# Quick check one point in intersection_points_df1\n",
    "print(intersection_points_df1.loc[1,'geometry'].x)\n",
    "print(intersection_points_df1.loc[1,'geometry'].y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "2a3f2965-4b70-448c-842c-980f7c7f1127",
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[22], line 22\u001b[0m\n\u001b[0;32m     20\u001b[0m \u001b[38;5;66;03m# the formula of Haversine distance\u001b[39;00m\n\u001b[0;32m     21\u001b[0m r \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m6371000\u001b[39m\n\u001b[1;32m---> 22\u001b[0m d \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m2\u001b[39m \u001b[38;5;241m*\u001b[39m r \u001b[38;5;241m*\u001b[39m arcsin(sqrt(sin((radians(intersection_points_df1\u001b[38;5;241m.\u001b[39mloc[\u001b[38;5;241m1\u001b[39m,\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgeometry\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mx) \u001b[38;5;241m-\u001b[39m radians(crime_selected_df1\u001b[38;5;241m.\u001b[39mloc[\u001b[38;5;241m1\u001b[39m,\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgeometry\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39my))\u001b[38;5;241m/\u001b[39m\u001b[38;5;241m2\u001b[39m)\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m \u001b[38;5;241m+\u001b[39m cos(radians(intersection_points_df1\u001b[38;5;241m.\u001b[39mloc[\u001b[38;5;241m1\u001b[39m,\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgeometry\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mx)) \u001b[38;5;241m*\u001b[39m cos(radians(crime_selected_df1\u001b[38;5;241m.\u001b[39mloc[\u001b[38;5;241m1\u001b[39m,\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgeometry\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39my)) \u001b[38;5;241m*\u001b[39m sin((radians(\u001b[43mintersection_points_df1\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mloc\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m,\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mgeometry\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[38;5;241m.\u001b[39my) \u001b[38;5;241m-\u001b[39m radians(crime_selected_df1\u001b[38;5;241m.\u001b[39mloc[\u001b[38;5;241m1\u001b[39m,\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgeometry\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mx))\u001b[38;5;241m/\u001b[39m\u001b[38;5;241m2\u001b[39m)\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m))\n\u001b[0;32m     24\u001b[0m \u001b[38;5;66;03m# Count the number of crime and Calculate the distance\u001b[39;00m\n\u001b[0;32m     25\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m d \u001b[38;5;241m<\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m50\u001b[39m:\n",
      "File \u001b[1;32m~\\anaconda3\\envs\\myenv\\lib\\site-packages\\pandas\\core\\indexing.py:960\u001b[0m, in \u001b[0;36m_LocationIndexer.__getitem__\u001b[1;34m(self, key)\u001b[0m\n\u001b[0;32m    958\u001b[0m     key \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mtuple\u001b[39m(com\u001b[38;5;241m.\u001b[39mapply_if_callable(x, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mobj) \u001b[38;5;28;01mfor\u001b[39;00m x \u001b[38;5;129;01min\u001b[39;00m key)\n\u001b[0;32m    959\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_is_scalar_access(key):\n\u001b[1;32m--> 960\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mobj\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_get_value\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkey\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtakeable\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_takeable\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    961\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_getitem_tuple(key)\n\u001b[0;32m    962\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m    963\u001b[0m     \u001b[38;5;66;03m# we by definition only have the 0th axis\u001b[39;00m\n",
      "File \u001b[1;32m~\\anaconda3\\envs\\myenv\\lib\\site-packages\\pandas\\core\\frame.py:3623\u001b[0m, in \u001b[0;36mDataFrame._get_value\u001b[1;34m(self, index, col, takeable)\u001b[0m\n\u001b[0;32m   3618\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mindex, MultiIndex):\n\u001b[0;32m   3619\u001b[0m     \u001b[38;5;66;03m# CategoricalIndex: Trying to use the engine fastpath may give incorrect\u001b[39;00m\n\u001b[0;32m   3620\u001b[0m     \u001b[38;5;66;03m#  results if our categories are integers that dont match our codes\u001b[39;00m\n\u001b[0;32m   3621\u001b[0m     \u001b[38;5;66;03m# IntervalIndex: IntervalTree has no get_loc\u001b[39;00m\n\u001b[0;32m   3622\u001b[0m     row \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mindex\u001b[38;5;241m.\u001b[39mget_loc(index)\n\u001b[1;32m-> 3623\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mseries\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_values\u001b[49m\u001b[43m[\u001b[49m\u001b[43mrow\u001b[49m\u001b[43m]\u001b[49m\n\u001b[0;32m   3625\u001b[0m \u001b[38;5;66;03m# For MultiIndex going through engine effectively restricts us to\u001b[39;00m\n\u001b[0;32m   3626\u001b[0m \u001b[38;5;66;03m#  same-length tuples; see test_get_set_value_no_partial_indexing\u001b[39;00m\n\u001b[0;32m   3627\u001b[0m loc \u001b[38;5;241m=\u001b[39m engine\u001b[38;5;241m.\u001b[39mget_loc(index)\n",
      "File \u001b[1;32m~\\anaconda3\\envs\\myenv\\lib\\site-packages\\geopandas\\array.py:360\u001b[0m, in \u001b[0;36mGeometryArray.__getitem__\u001b[1;34m(self, idx)\u001b[0m\n\u001b[0;32m    358\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__getitem__\u001b[39m(\u001b[38;5;28mself\u001b[39m, idx):\n\u001b[0;32m    359\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(idx, numbers\u001b[38;5;241m.\u001b[39mIntegral):\n\u001b[1;32m--> 360\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_geom_to_shapely\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdata\u001b[49m\u001b[43m[\u001b[49m\u001b[43midx\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    361\u001b[0m     \u001b[38;5;66;03m# array-like, slice\u001b[39;00m\n\u001b[0;32m    362\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m compat\u001b[38;5;241m.\u001b[39mPANDAS_GE_10:\n\u001b[0;32m    363\u001b[0m         \u001b[38;5;66;03m# for pandas >= 1.0, validate and convert IntegerArray/BooleanArray\u001b[39;00m\n\u001b[0;32m    364\u001b[0m         \u001b[38;5;66;03m# to numpy array, pass-through non-array-like indexers\u001b[39;00m\n",
      "File \u001b[1;32m~\\anaconda3\\envs\\myenv\\lib\\site-packages\\geopandas\\array.py:115\u001b[0m, in \u001b[0;36m_geom_to_shapely\u001b[1;34m(geom)\u001b[0m\n\u001b[0;32m    113\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m geom\n\u001b[0;32m    114\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m--> 115\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mvectorized\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_pygeos_to_shapely\u001b[49m\u001b[43m(\u001b[49m\u001b[43mgeom\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32m~\\anaconda3\\envs\\myenv\\lib\\site-packages\\geopandas\\_vectorized.py:71\u001b[0m, in \u001b[0;36m_pygeos_to_shapely\u001b[1;34m(geom)\u001b[0m\n\u001b[0;32m     69\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m compat\u001b[38;5;241m.\u001b[39mPYGEOS_SHAPELY_COMPAT:\n\u001b[0;32m     70\u001b[0m     geom \u001b[38;5;241m=\u001b[39m shapely\u001b[38;5;241m.\u001b[39mgeos\u001b[38;5;241m.\u001b[39mlgeos\u001b[38;5;241m.\u001b[39mGEOSGeom_clone(geom\u001b[38;5;241m.\u001b[39m_ptr)\n\u001b[1;32m---> 71\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mshapely\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mgeometry\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mbase\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mgeom_factory\u001b[49m\u001b[43m(\u001b[49m\u001b[43mgeom\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m     73\u001b[0m \u001b[38;5;66;03m# fallback going through WKB\u001b[39;00m\n\u001b[0;32m     74\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m pygeos\u001b[38;5;241m.\u001b[39mis_empty(geom) \u001b[38;5;129;01mand\u001b[39;00m pygeos\u001b[38;5;241m.\u001b[39mget_type_id(geom) \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[0;32m     75\u001b[0m     \u001b[38;5;66;03m# empty point does not roundtrip through WKB\u001b[39;00m\n",
      "File \u001b[1;32m~\\anaconda3\\envs\\myenv\\lib\\site-packages\\shapely\\geometry\\base.py:86\u001b[0m, in \u001b[0;36mgeom_factory\u001b[1;34m(g, parent)\u001b[0m\n\u001b[0;32m     84\u001b[0m ob\u001b[38;5;241m.\u001b[39m_set_geom(g)\n\u001b[0;32m     85\u001b[0m ob\u001b[38;5;241m.\u001b[39m__p__ \u001b[38;5;241m=\u001b[39m parent\n\u001b[1;32m---> 86\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[43mlgeos\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmethods\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mhas_z\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m(\u001b[49m\u001b[43mg\u001b[49m\u001b[43m)\u001b[49m:\n\u001b[0;32m     87\u001b[0m     ob\u001b[38;5;241m.\u001b[39m_ndim \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m3\u001b[39m\n\u001b[0;32m     88\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n",
      "File \u001b[1;32m~\\anaconda3\\envs\\myenv\\lib\\site-packages\\shapely\\geos.py:573\u001b[0m, in \u001b[0;36merrcheck_predicate\u001b[1;34m(result, func, argtuple)\u001b[0m\n\u001b[0;32m    566\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m TopologicalError(\n\u001b[0;32m    567\u001b[0m             \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mThe operation \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m could not be performed.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m    568\u001b[0m             \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mLikely cause is invalidity of the geometry.\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39mformat(\n\u001b[0;32m    569\u001b[0m                 func\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m))\n\u001b[0;32m    570\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m errcheck_just_free(result, func, argtuple)\n\u001b[1;32m--> 573\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21merrcheck_predicate\u001b[39m(result, func, argtuple):\n\u001b[0;32m    574\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Result is 2 on exception, 1 on True, 0 on False\"\"\"\u001b[39;00m\n\u001b[0;32m    575\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m result \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m2\u001b[39m:\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "from numpy import sqrt, radians, arcsin, sin, cos\n",
    "\n",
    "start = timeit.default_timer()\n",
    "\n",
    "points = pd.DataFrame()\n",
    "\n",
    "# Load the intersection points\n",
    "for a in intersection_points_df1.index:\n",
    "    # i is the interaction point we use for calculate.\n",
    "    i = intersection_points_df1.loc[a,'geometry']\n",
    "    # j is used to count the number of crime points that meets the requirements.\n",
    "    j = 1\n",
    "\n",
    "    # Load the crime points\n",
    "    for b in crime_selected_df1.index:\n",
    "        # c is the crime point we use for calculate.\n",
    "        c = crime_selected_df1.loc[b,'geometry']\n",
    "        # z is used to show the distance that is from the crime point to the interaction point.\n",
    "        z = crime_selected_df1.loc[b,'Distance']\n",
    "        # the formula of Haversine distance\n",
    "        r = 6371000\n",
    "        d = 2 * r * arcsin(sqrt(sin((radians(intersection_points_df1.loc[1,'geometry'].x) - radians(crime_selected_df1.loc[1,'geometry'].y))/2)**2 + cos(radians(intersection_points_df1.loc[1,'geometry'].x)) * cos(radians(crime_selected_df1.loc[1,'geometry'].y)) * sin((radians(intersection_points_df1.loc[1,'geometry'].y) - radians(crime_selected_df1.loc[1,'geometry'].x))/2)**2))\n",
    "        \n",
    "        # Count the number of crime and Calculate the distance\n",
    "        if d <= 50:\n",
    "            j = j + 1\n",
    "            z = z + d\n",
    "            \n",
    "data = [[i,j]]\n",
    "points = points.append(data,ignore_index = True)\n",
    "    \n",
    "stop = timeit.default_timer()\n",
    "print('Time: ', stop - start) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "755b4752-993f-4b14-a0c0-b929a55d4915",
   "metadata": {},
   "source": [
    "## 2.2 Commerical Land Use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec9a9871-6cb6-4469-89b4-736796230409",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter points after the bounding box\n",
    "corner1 = Point(-87.633, 41.736) \n",
    "corner2 = Point(-87.633, 41.816)\n",
    "corner3 = Point(-87.43, 41.816)\n",
    "corner4 = Point(-87.43, 41.736)\n",
    "\n",
    "bbox_polygon = Polygon([corner1, corner2, corner3, corner4])\n",
    "bbox_gdf = gpd.GeoDataFrame(gpd.GeoSeries(bbox_polygon), \n",
    "                            columns=['geometry'], crs='EPSG:4326')\n",
    "\n",
    "commercials_selected = gpd.overlay(gdf, bbox_gdf, how='intersection')\n",
    "\n",
    "# Save the selected commercials\n",
    "commercials_selected.to_file('commercials_selected.geojson', driver=\"GeoJSON\")\n",
    "\n",
    "commercials_selected_gdf = gpd.read_file('commercials_selected.geojson')\n",
    "\n",
    "# Read final interaction points\n",
    "points = pd.read_csv('points_final.csv')\n",
    "intersection_points = pd.DataFrame(points)\n",
    "intersection_points.rename(columns={'Unnamed: 0':'ID'}, inplace=True)\n",
    "\n",
    "# Create two lists for the loop results to be placed\n",
    "lat = []\n",
    "lon = []\n",
    "\n",
    "# For each row in a varible,\n",
    "for row in intersection_points['geometry']:\n",
    "    row = \"\".join(row)\n",
    "    # Try to,\n",
    "    try:\n",
    "        # Split the row by comma and append\n",
    "        # everything before the comma to lat\n",
    "        lat.append(row.split(',')[0])\n",
    "        # Split the row by comma and append\n",
    "        # everything after the comma to lon\n",
    "        lon.append(row.split(',')[1])\n",
    "    # But if you get an error\n",
    "    except:\n",
    "        # append a missing value to lat\n",
    "        lat.append(np.NaN)\n",
    "        # append a missing value to lon\n",
    "        lon.append(np.NaN)\n",
    "\n",
    "# Create two new columns from lat and lon\n",
    "intersection_points['Latitude'] = lat\n",
    "intersection_points['Longitude'] = lon\n",
    "\n",
    "# Clean the two new columns\n",
    "intersection_points['Latitude'] = list(set([x.replace('[', '') for x in intersection_points['Latitude']]))\n",
    "intersection_points['Longitude'] = list(set([x.replace(']', '') for x in intersection_points['Longitude']]))\n",
    "\n",
    "# Built the geodataframe type of intersection_points\n",
    "intersection_points_df1 = gpd.GeoDataFrame(\n",
    "    intersection_points, geometry=gpd.points_from_xy(intersection_points['Longitude'], intersection_points['Latitude']))\n",
    "\n",
    "# Quick check one point in intersection_points_df1\n",
    "print(intersection_points_df1.loc[1,'geometry'].x)\n",
    "print(intersection_points_df1.loc[1,'geometry'].y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e34c8a30-5d49-46c0-a899-3879083d00dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import sqrt, radians, arcsin, sin, cos\n",
    "\n",
    "start = timeit.default_timer()\n",
    "\n",
    "points = pd.DataFrame()\n",
    "\n",
    "for a in tqdm(intersection_points_df1.index):\n",
    "    # i is the interaction point we use for calculation.\n",
    "    i = intersection_points_df1.loc[a,'geometry']\n",
    "    # j is used to count the number of commercial points that meets the requirements.\n",
    "    j = 1\n",
    "\n",
    "    for b in commercials_selected_gdf.index:\n",
    "        # c is the commercial point we use for calculation.\n",
    "        c = commercials_selected_gdf.loc[b,'geometry']\n",
    "        # z is used to show the distance that is from the commercal point to the interaction point.\n",
    "        z = commercials_selected_gdf.loc[b,'Distance']\n",
    "        # the formula of Haversine distance\n",
    "        r = 6371000\n",
    "        d = 2 * r * arcsin(sqrt(sin((radians(intersection_points_df1.loc[1,'geometry'].x) \\\n",
    "            - radians(commercials_selected_gdf.loc[1,'geometry'].y))/2)**2 \\\n",
    "            + cos(radians(intersection_points_df1.loc[1,'geometry'].x)) \\\n",
    "            * cos(radians(commercials_selected_gdf.loc[1,'geometry'].y)) \\\n",
    "            * sin((radians(intersection_points_df1.loc[1,'geometry'].y) \\\n",
    "            - radians(commercials_selected_gdf.loc[1,'geometry'].x))/2)**2))\n",
    "        \n",
    "        # Count the number of commercials and Calculate the distance\n",
    "        if d <= 50:\n",
    "            j = j + 1\n",
    "            z = z + d\n",
    "            \n",
    "data = [[i,j]]\n",
    "points = points.append(data,ignore_index = True)\n",
    "    \n",
    "stop = timeit.default_timer()\n",
    "print('Time: ', stop - start) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4032ec5-6fa0-4a0e-a4ce-090acdbcef57",
   "metadata": {},
   "source": [
    "# Part 3 Calculate the number of population\n",
    "In this part, we will calculate the population number of interanction points. We mainly use the 'points in polygon'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "032152e7-e47d-4eac-9348-938b6afeba2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the environment\n",
    "import random\n",
    "from shapely.geometry import Point, Polygon\n",
    "\n",
    "# Import the module that I wrote\n",
    "from geospatial import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94049bd9-54fe-4fbd-97bb-803702919b97",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read data\n",
    "df_demog = pd.read_json(\"demographic.json\")\n",
    "df_points_orig = pd.read_json('points_weight_commercial.json')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e204404-f440-4078-936a-9651b35dd09f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get color for each demographic region\n",
    "np.random.seed(22)\n",
    "colors = [[np.random.uniform(0, 1) for rep in range(3)] for region in range(df_demog.shape[0])]\n",
    "\n",
    "# Convert coordinates of points to numbers\n",
    "df_points = pd.DataFrame(df_points_orig)\n",
    "df_points['geometry'] = df_points['geometry'].astype(str)\n",
    "\n",
    "# Create two lists for the loop results to be placed\n",
    "lat = []\n",
    "lon = []\n",
    "\n",
    "# For each row in a varible,\n",
    "for row in df_points['geometry']:\n",
    "    row = \"\".join(row)\n",
    "    # Try to,\n",
    "    try:\n",
    "        # Split the row by comma and append\n",
    "        # everything before the comma to lat\n",
    "        lat.append(row.split(',')[0])\n",
    "        # Split the row by comma and append\n",
    "        # everything after the comma to lon\n",
    "        lon.append(row.split(',')[1])\n",
    "    # But if you get an error\n",
    "    except:\n",
    "        # append a missing value to lat\n",
    "        lat.append(np.NaN)\n",
    "        # append a missing value to lon\n",
    "        lon.append(np.NaN)\n",
    "\n",
    "# Create two new columns from lat and lon\n",
    "df_points['Latitude'] = lat\n",
    "df_points['Longitude'] = lon\n",
    "\n",
    "# Clean the two new columns\n",
    "df_points['Latitude'] = list(set([x.replace('[', '') for x in df_points['Latitude']]))\n",
    "df_points['Longitude'] = list(set([x.replace(']', '') for x in df_points['Longitude']]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71fbf4ad-c77f-4e7e-982f-6ec48f750c49",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get polygons from coordinates of each demographic area\n",
    "polygons = []\n",
    "\n",
    "# Iterate over regions\n",
    "for index, row in df_demog.iterrows():\n",
    "    for i in range(len(row[\"coordinates\"])):\n",
    "        for j in range(len(row[\"coordinates\"][i])):\n",
    "            points_poly = np.array(row[\"coordinates\"][i][j])\n",
    "            polygons.append((index, Polygon(points_poly)))\n",
    "    \n",
    "# Plot polygons\n",
    "plt.figure(figsize=(10, 10), dpi=80)\n",
    "for index, poly in polygons:\n",
    "    x,y = poly.exterior.xy\n",
    "    plt.plot(x,y, color=colors[index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b9620f5-4ce5-4b77-8ad5-769efc8e3e0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "points = []\n",
    "\n",
    "# Check which demographic area corresponds to each point\n",
    "for point_idx, row in df_points.iterrows(): # Create 1000 random Points to test\n",
    "    x = float(row['Latitude'])\n",
    "    y = float(row['Longitude'])\n",
    "    p = Point(x,y)\n",
    "    \n",
    "    # Find the corresponding polygon\n",
    "    corr_poly_idx = -1\n",
    "    for poly_idx, poly in polygons:\n",
    "        if (poly.contains(p)):\n",
    "            corr_poly_idx = poly_idx\n",
    "            break;\n",
    "            \n",
    "    points.append((point_idx, p, corr_poly_idx))\n",
    "            \n",
    "# Plot points colored by their corresponding demographic region\n",
    "# plt.figure(figsize=(10, 10), dpi=80)\n",
    "points_colors = [colors[entry[2]] for entry in points]\n",
    "plt.scatter(x = [entry[1].x for entry in points],\n",
    "            y = [entry[1].y for entry in points], \n",
    "            c = points_colors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "834fd20e-dc2c-481e-98ab-7becf18c2ad9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Append demographic info to points\n",
    "df_points['area'] = 0\n",
    "df_points['population'] = 0\n",
    "\n",
    "for point_idx, _, poly_idx in points:\n",
    "    if poly_idx == -1:\n",
    "        continue\n",
    "        \n",
    "    df_points.loc[point_idx, 'area'] = df_demog.loc[poly_idx, 'area']\n",
    "    df_points.loc[point_idx, 'population'] = df_demog.loc[poly_idx, 'population']\n",
    "\n",
    "# Calculation of population density\n",
    "df_points = df_points.drop(['Latitude','Longitude'],axis =1)\n",
    "df_points['pop/km2']=(df_points['population']/df_points['area'])*100000/0.0929\n",
    "df_points\n",
    "\n",
    "# Save final dataframe\n",
    "df_points.to_json('points_all_x.json')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39c8b9e0-426e-4859-9283-ba1b47004868",
   "metadata": {},
   "source": [
    "# Part 4: Correlation Analysis and Visualization\n",
    "\n",
    "In this part, we will calculate the correlation coffincient of crime, commerical land use, population, and street weight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "21260982-121a-4eb5-9ed9-11923b16e4be",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the function for correlation analysis\n",
    "def corr(df,factors):\n",
    "    corrmat = pd.DataFrame(columns = factors, index = factors)\n",
    "    \n",
    "    for i in factors:\n",
    "        for j in factors:\n",
    "            de = 0\n",
    "            mo1 = 0\n",
    "            mo2 = 0\n",
    "            for k in df.index:\n",
    "                de = de + (df[i][k]-np.mean(df[i]))*(df[j][k]-np.mean(df[j]))\n",
    "                mo1 = mo1 + (df[i][k]-np.mean(df[i]))*(df[i][k]-np.mean(df[i]))\n",
    "                mo2 = mo2 + (df[j][k]-np.mean(df[j]))*(df[j][k]-np.mean(df[j]))\n",
    "            mo = np.sqrt(mo1)*np.sqrt(mo2)\n",
    "            correlation = round(de/mo,4)\n",
    "            corrmat[i][j]=correlation\n",
    "            \n",
    "    return corrmat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "68859778-0cad-4850-8a84-f9c1df479e08",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geometry</th>\n",
       "      <th>weight</th>\n",
       "      <th>commercial</th>\n",
       "      <th>pop/km2</th>\n",
       "      <th>Crime</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[-87.6540413653, 41.7623276297]</td>\n",
       "      <td>8</td>\n",
       "      <td>0</td>\n",
       "      <td>364.488762</td>\n",
       "      <td>3669</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[-87.6540205024, 41.7614168007]</td>\n",
       "      <td>10</td>\n",
       "      <td>0</td>\n",
       "      <td>561.012857</td>\n",
       "      <td>4236</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[-87.7317860919, 41.7442874546]</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>310.413612</td>\n",
       "      <td>913</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>[-87.732959782, 41.7449095243]</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>310.413612</td>\n",
       "      <td>2473</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>[-87.7129841335, 41.7715200848]</td>\n",
       "      <td>8</td>\n",
       "      <td>0</td>\n",
       "      <td>310.413612</td>\n",
       "      <td>646</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3016</th>\n",
       "      <td>[-87.6559513961, 41.7886882626]</td>\n",
       "      <td>8</td>\n",
       "      <td>0</td>\n",
       "      <td>456.774667</td>\n",
       "      <td>1119</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3017</th>\n",
       "      <td>[-87.6426658423, 41.7761208111]</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>310.413612</td>\n",
       "      <td>912</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3018</th>\n",
       "      <td>[-87.6619561462, 41.7558364448]</td>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>310.413612</td>\n",
       "      <td>2778</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3019</th>\n",
       "      <td>[-87.7131122899, 41.793367710300004]</td>\n",
       "      <td>10</td>\n",
       "      <td>0</td>\n",
       "      <td>307.256734</td>\n",
       "      <td>4328</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3020</th>\n",
       "      <td>[-87.6391901422, 41.7434340699]</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>344.295033</td>\n",
       "      <td>1710</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3021 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                  geometry  weight  commercial     pop/km2  \\\n",
       "0          [-87.6540413653, 41.7623276297]       8           0  364.488762   \n",
       "1          [-87.6540205024, 41.7614168007]      10           0  561.012857   \n",
       "2          [-87.7317860919, 41.7442874546]       6           0  310.413612   \n",
       "3           [-87.732959782, 41.7449095243]       6           0  310.413612   \n",
       "4          [-87.7129841335, 41.7715200848]       8           0  310.413612   \n",
       "...                                    ...     ...         ...         ...   \n",
       "3016       [-87.6559513961, 41.7886882626]       8           0  456.774667   \n",
       "3017       [-87.6426658423, 41.7761208111]       6           0  310.413612   \n",
       "3018       [-87.6619561462, 41.7558364448]       7           0  310.413612   \n",
       "3019  [-87.7131122899, 41.793367710300004]      10           0  307.256734   \n",
       "3020       [-87.6391901422, 41.7434340699]       5           0  344.295033   \n",
       "\n",
       "      Crime  \n",
       "0      3669  \n",
       "1      4236  \n",
       "2       913  \n",
       "3      2473  \n",
       "4       646  \n",
       "...     ...  \n",
       "3016   1119  \n",
       "3017    912  \n",
       "3018   2778  \n",
       "3019   4328  \n",
       "3020   1710  \n",
       "\n",
       "[3021 rows x 5 columns]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load the dataset and extract useful informations\n",
    "points = pd.read_json(\"points_all_x.json\")\n",
    "crime = pd.read_csv(\"intersection_points_crime.csv\")\n",
    "points = pd.merge(points, crime, left_index = True, right_index = True)\n",
    "\n",
    "points = points.drop(['Unnamed: 0','ID','geometry_y','weight_y','Latitude','Longitude','area','population'],axis =1)\n",
    "points = points.rename(columns={'geometry_x':'geometry','weight_x':'weight'})\n",
    "\n",
    "points.dropna(inplace=True)\n",
    "points.reset_index(inplace=True)\n",
    "points = points.drop(['index'],axis =1)\n",
    "\n",
    "points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "217e8a65-deb2-4979-8a61-d6817c5e486b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Crime', 'weight', 'commercial', 'pop/km2']"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list = ['Crime','weight','commercial','pop/km2']\n",
    "list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "43a622b6-8645-4e0b-a920-e410e7ab2942",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Crime</th>\n",
       "      <th>weight</th>\n",
       "      <th>commercial</th>\n",
       "      <th>pop/km2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Crime</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0013</td>\n",
       "      <td>-0.0054</td>\n",
       "      <td>-0.0211</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>weight</th>\n",
       "      <td>0.0013</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.2647</td>\n",
       "      <td>0.0001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>commercial</th>\n",
       "      <td>-0.0054</td>\n",
       "      <td>0.2647</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-0.0257</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pop/km2</th>\n",
       "      <td>-0.0211</td>\n",
       "      <td>0.0001</td>\n",
       "      <td>-0.0257</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             Crime  weight commercial pop/km2\n",
       "Crime          1.0  0.0013    -0.0054 -0.0211\n",
       "weight      0.0013     1.0     0.2647  0.0001\n",
       "commercial -0.0054  0.2647        1.0 -0.0257\n",
       "pop/km2    -0.0211  0.0001    -0.0257     1.0"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "corrmat = corr(points,list)\n",
    "corrmat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "24d0ac62-cd1b-4613-9238-7241f92e5419",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAALOCAYAAACgZ3uKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB4QklEQVR4nO3dd3gU5RbH8d8mkIQeSEKk9xI6hJYEpEkVAUVE6V2kEzpKL1EUDEWQjghIk6LAFZAm0qSDgiKI0oKQ0IsBkrl/RBaWJJBdNiTrfj/32eeSd9+ZPbsZJycnZ94xGYZhCAAAAHACLkkdAAAAAPCikPwCAADAaZD8AgAAwGmQ/AIAAMBpkPwCAADAaZD8AgAAwGmQ/AIAAMBpkPwCAADAaZD8AgAAwGmQ/AIAAMBpkPwCAADghfvhhx/02muvKWvWrDKZTFq1atUzt9m2bZv8/f3l4eGhvHnz6vPPP7f6dUl+AQAA8MLdvn1bJUuW1JQpUxI0//Tp06pXr54qV66sgwcPavDgwerRo4e+/vprq17XZBiGYUvAAAAAgD2YTCatXLlSjRo1infOgAED9M033+j48ePmsc6dO+vw4cPatWtXgl+Lyi8AAACeW2RkpG7cuGHxiIyMtNv+d+3apVq1almM1a5dW/v27dP9+/cTvJ8UdosIAAAAL1yq0t2SOgRJ0oCG3hoxYoTF2LBhwzR8+HC77P/ixYvy9fW1GPP19dWDBw8UHh6uLFmyJGg/JL8AAAB4boMGDVJwcLDFmLu7u11fw2QyWXz9sHv3yfGnIfkFAADAc3N3d7d7svu4l156SRcvXrQYu3TpklKkSCEvL68E74fkFwAAwJGZnOMSroCAAH377bcWYxs2bFDZsmWVMmXKBO/HOT4tAAAAJCu3bt3SoUOHdOjQIUkxS5kdOnRIZ86ckRTTRtGqVSvz/M6dO+uvv/5ScHCwjh8/rjlz5mj27Nnq27evVa9L5RcAAMCRWdHvmpzs27dP1apVM3/9sF+4devWmjdvnsLCwsyJsCTlyZNH69atU+/evfXZZ58pa9asmjRpkho3bmzV67LOLwAAgANL5d8zqUOQJN3dPzGpQ0gQ2h4AAADgNGh7AAAAcGROcsGbvfBpAQAAwGmQ/AIAAMBp0PYAAADgyBx0tYekQuUXAAAAToPkFwAAAE6DtgcAAABHxmoPVuHTAgAAgNOg8gsAAODIuODNKlR+AQAA4DRIfgEAAOA0aHsAAABwZFzwZhU+LQAAADgNkl8AAAA4DdoeAAAAHBmrPViFyi8AAACcBpVfAAAAR8YFb1bh0wIAAIDTIPkFAACA06DtAQAAwJFxwZtVqPwCAADAaZD8AgAAwGnQ9gAAAODIWO3BKnxaAAAAcBpUfgEAABwZF7xZhcovAAAAnAbJLwAAAJwGbQ8AAACOjAverMKnBQAAAKdB8gsAAACnQdsDAACAI6PtwSp8WgAAAHAaJL8AAABwGrQ9AAAAODIXbnJhDSq/AAAAcBpUfgEAABwZF7xZhU8LAAAAToPkFwAAAE6DtgcAAABHZuKCN2tQ+QUAAIDTIPkFAACA06DtAQAAwJGx2oNV+LQAAADgNKj8AgAAODIueLMKlV8AAAA4DZJfAAAAOA3aHgAAABwZF7xZhU8LAAAAToPkFwAAAE6DtgcAAABHxmoPVqHyCwAAAKdB8gsAAACnQdsDAACAI2O1B6vwaQEAAMBpUPkFAABwZFzwZhUqvwAAAHAaJL8AAABwGrQ9AAAAODIueLMKnxYAAACcBskvAAAAnEayaXtIVbpbUocAB7Hmq+FJHQIcSPk8mZI6BAD/Qenck1H9kNUerJKMvnMAAABA4ko2lV8AAADYgAverMKnBQAAAKdB8gsAAACnQdsDAACAI6PtwSp8WgAAAHAaJL8AAABwGrQ9AAAAODLW+bUKlV8AAAA4DZJfAAAAOA3aHgAAABwZqz1YhU8LAAAAToPKLwAAgCPjgjerUPkFAACA0yD5BQAAgNOg7QEAAMCRccGbVfi0AAAA4DRIfgEAAOA0aHsAAABwZKz2YBUqvwAAAHAaVH4BAAAcmInKr1Wo/AIAAMBpkPwCAADAadD2AAAA4MBoe7AOlV8AAAA4DZJfAAAAOA3aHgAAABwZXQ9WofILAAAAp0HlFwAAwIFxwZt1qPwCAADAaZD8AgAAwGnQ9gAAAODAaHuwDpVfAAAAOA2SXwAAADgN2h4AAAAcGG0P1qHyCwAAAKdB8gsAAACnQdsDAACAA6PtwTpUfgEAAOA0qPwCAAA4Mgq/VqHyCwAAAKdB8gsAAACnQdsDAACAA+OCN+tQ+QUAAIDTIPkFAACA06DtAQAAwIHR9mAdKr8AAABwGjYlvyNHjtSdO3dijd+9e1cjR4587qAAAACQMCaTKVk8HIVNye+IESN069atWON37tzRiBEjnjsoAAAAIDHYlPwahhFnhn/48GFlypTpuYMCAAAAEoNVF7xlzJjRXNouWLCgRQIcFRWlW7duqXPnznYPEgAAAHFzpJaD5MCq5Dc0NFSGYahdu3YaMWKEMmTIYH7Ozc1NuXPnVkBAgN2DBAAAAOzBquS3devWkqQ8efIoMDBQKVOmTJSgAAAAgMRg0zq/VapUUXR0tE6cOKFLly4pOjra4vmXX37ZLsEBAADgGeh6sIpNye/u3bvVrFkz/fXXXzIMw+I5k8mkqKgouwQHAAAA2JNNyW/nzp1VtmxZrV27VlmyZKHRGgAAAA7BpuT3999/1/Lly5U/f357xwMAAAArUIS0jk3r/FaoUEEnT560dywAAABAokpw5ffIkSPmf3fv3l19+vTRxYsXVbx48VirPpQoUcJ+EQIAACBeVH6tk+Dkt1SpUjKZTBYXuLVr187874fPccEbAAAAkqsEJ7+nT59OzDgAAACARJfg5DdXrlyJGQcAAABsQNuDdWxa7eGbb76Jc9xkMsnDw0P58+dXnjx5niswAAAA/LdNnTpVH3/8scLCwlS0aFGFhoaqcuXK8c5fuHChxo0bp99//10ZMmRQnTp19Mknn8jLyyvBr2lT8tuoUaNY/b+SZd9vpUqVtGrVKmXMmNGWl3BqQWXyqXerV1SmSE5l8cmgt3rP0Ldbjzx7Q/xnbFu3Qt+vXKTrVyOUJWceNWnfQ/mLlopz7sFdW7X9fyt17vRJPbh/T1ly5tGrb7dXkTIVLObduXVT3yyYoUO7t+nOrZvy8s2ixm27qVjZwBfwjmCNS3//rcmh47Xzxx/0T2SkcuXKrSEjRsuvSNE45x86sF+TQsfrr9N/6J9//tFLWbLqjSZvqXnLNs8dy8WwC/pozCjt/WmP3D3cVadeffXq008pU7pJki6cP68GdV+Jtd2kqTMUWCn+H2Cwv2WLF+nLeXMUHn5ZefPlV5/+g1Tav2y88/fv+0mffvyR/jh1Uj4+mdWybXu9+dbbFnM2bdygzz+bpHNnzyh7jpzq0r2nqtWoaX7+wL69+nLeHB0//ovCL1/WJ6GTVbW65fEwfeoUbfhunf6+eFEpU6aUX5Ei6tK9l4qVKGnfDwAOZ8mSJerVq5emTp2qoKAgTZ8+XXXr1tWxY8eUM2fOWPN//PFHtWrVSp9++qlee+01nT9/Xp07d1aHDh20cuXKBL+uTUudbdy4UeXKldPGjRt1/fp1Xb9+XRs3blT58uW1Zs0a/fDDD4qIiFDfvn1t2b3TS5PKXUdPnFfvD5cmdShIAvu2f6/lsyeqTpNWGvTpXOUvUkKfjeyrK5cvxjn/5C+HVLhUeXUZ+okGTpijgsXLaNqY/jr7xwnznAf372vSsF6KuBSmjgNGa9jUr9S86wB5evm8qLeFBLpx47rat26mFClSaOLUGVq2co169emvdOnSxbuNR6pUeuvt5pox90stW7VW7Tt11rTJk7Ri+fOdQ6KiotSza2fdvXtXs75YoLHjxmvz9xv06SfjYs2dOmOOvtv8g/lRrkKFOPaIxLLhu3UaP+5Dtev4rhYuXaHSZfzVo8u7uhh2Ic7558+dU88unVW6jL8WLl2hth066ZMPx2rTxg3mOUcOH9Tg/sGqV7+Bvlq+SvXqN9DAfsH6+chh85y7d++qQKFC6j/og3hjy5Urt/oP/kCLV6zWrC8WKEvWbOrauYOuXrlivw/A2ZmSycNKEyZMUPv27dWhQwf5+fkpNDRUOXLk0LRp0+Kcv3v3buXOnVs9evRQnjx5VKlSJb377rvat2+fVa9rU+W3Z8+emjFjhgIDH1WMatSoIQ8PD3Xq1Em//PKLQkNDLVaDQMJt2HFMG3YcS+owkEQ2r16iwFfqK6hWA0lSkw69dOzgT/rhfyvVqNV7seY36dDL4uuGLTvryJ7tOvrTj8qRt6Akaef3a3Tn1g31+2i6XFPE/GfvlfmlxH0jsMkXc2bJ1zeLho0aax7Lmi3bU7cp7FdEhf2KWMzfsmmjDh7YpzfefMs8/s2qFZo/d7YunD+nLFmz6e1mLdTk7Wbx7nf3zh06/ccpTdkwSz6ZM0uSevXprxFDBqtL915KmzateW4GT095e/PLVFJZOP8LNXz9DTVq3ESS1GfAYO3auUPLly5Wt57BseZ/vWyxXsqSRX0GDJYk5cmbT8eO/aIFX8xRjZq1JElfLZivChUD1bZDJ0lS2w6ddGDfXi1aMF9jx42XJAVVfllBlV9+amx1Xq1v8XXvfgO1euXX+v3EbypfMeD53jiSlcjISEVGRlqMubu7y93dPdbce/fuaf/+/Ro4cKDFeK1atbRz58449x8YGKj3339f69atU926dXXp0iUtX75cr776qlVx2lT5PXXqlNKnTx9rPH369Prjjz8kSQUKFFB4eLgtuwec1oP793Xm1G/yK1XeYtyvVHn98evPCdpHdHS0/rl7V6nTPfpv9OjeH5WnUDEtnj5eA1rV16juLfTdsi8UzbKEyc4PW7fIr2hRDejTSzWrBKnZW29opZUV3F+PH9ORQ4fk71/OPLZy+VJNnRyqLt17atmqterao5c+/2yS1qxeFe9+jh45pHz5C5gTX0kKCKqke/fu6ddjv1jMDe7RVTWrBKldq2b6fsN6q+LF87l//55+Pf6LKgYGWYxXDAjSkUMH49zm6OFDqhhgOT8gMEjHjv2iB/fvS5KOHD6sCoGWbVEVg4J05HDc+0xorCuXL1XadOlUsFBhm/cDSyaTKVk8QkJClCFDBotHSEhInDGHh4crKipKvr6+FuO+vr66eDHuv3QGBgZq4cKFatq0qdzc3PTSSy/J09NTkydPturzsqny6+/vr379+mn+/Pny8Yn5Tf/y5cvq37+/ypWLOdn+/vvvyp49e5zbx/WbgREdJZOLqy3hAP8Zt25cU3R0lNJ5ZrIYT++ZUTeuRiRoH5tWfaV7kXflH1TDPBZ+8YIiLh1QuSq11HXoJ7p04ZyWzBiv6Kgo1Xubv9AkJ+fPndXXSxerecs2atuhk375+ag++WisUrq5qX6DRk/dtt4rVXX16hVFRUWp03tdzVVASZo143P16jNA1V+Jqeply55df/xxSiuWL1H9hnHvNyI8XJmeuIgkffoMSpkypbm4kTp1avXuN0AlS5WRi4uLfti6WYP7B+vevRDVq9/A9g8CCXbt6jVFRUUpk5e3xXgmL694i1AREbG/t5m8vBX14IGuXbsqb5/MiggPl1cmy316ZfJWhA2Fre3btmhw/77655+78vbx0WfTZ8uTa4L+cwYNGqTgYMu/NMRV9X3ckytVPLx2LC7Hjh1Tjx49NHToUNWuXVthYWHq16+fOnfurNmzZyc4TpuS39mzZ6thw4bKnj27cuTIIZPJpDNnzihv3rxavXq1JOnWrVsaMmRInNuHhIRoxIgRFmOuvuWUMkv5OOcDzib2ySBhS9ns/WGj1i6eo86DP1Q6z0c/WAzDULoMGdW8S3+5uLoqZ/7Cun41XBtXLiL5TWaiow0VKVpUXXv2lhTT0vDHqZP6euniZya/M+ct0N07d3T0yCFNmThB2XPkUp16r+rqlSv6+2KYRg3/QGNGDDXPj4p6oLRpY3qJe7zXSQcP7JckZcmaRUtXrpEU93EXczzG/NszY0aLC+uKFC2mGzduaP7c2SS/L9iT36qnJREx82MnHbF29OQcPX2f8SlbroIWLVuha1evauWKZRrUt7fmLVwSKwGHY4uvxSEu3t7ecnV1jVXlvXTpUqxq8EMhISEKCgpSv379JMXcUThNmjSqXLmyRo8erSxZsiTotW1KfgsVKqTjx49r/fr1OnHihAzDUOHChVWzZk25uMR0UjRq1Cje7eP6zSBz5QG2hAL8p6RN7ykXF9dYVd6b16/GqgY/ad/277Vgcog6DBitwqXKWTyXPqOXXF1TyMX10V9XXsqeSzeuRujB/ftK8cQtypF0vH28lSdvPouxPHnyavP3G+LZ4pFs//61LX/Bgoq4EqEZ06aoTr1XFW1ES5I+GDZSxYpb3n7e5d+/uH0wfJQiI/+RJKV42Bfu7a2fj1quNHPjxnU9eHBfXk9UGR9XvERJrV6x/Jnxwj48M3rK1dU1VkX26pUr8S7/5OUVu4J79UqEXFOkkGcGz5g53t6KiLhsMefKlQibEtZUqVMrR85cypEzl4qXLKXX69fW6pVfm/uJ8XwccZ1fNzc3+fv7a+PGjXr99dfN4xs3blTDhg3j3ObOnTvm89NDrv/+XHtyBbKnsSn5lWI+6Dp16qhOnTpWbxvXbwa0PABSipQplTNfIR0/vFelAqqYx389tFclKlSKd7u9P2zUgslj1bbPCBWPY+myfH7FtfeHjYqOjjb/gvr3hbPKkNGLxDeZKVmqjP7680+Lsb/++lNZsmS1bkeGofv370mKSXQyZ/bV+XNnVffV1+KcnjmOSkvxEqU0Z+Z0hV++JG+fmL7f3Tt3yM3NTYXjWXZNkn47fpyL316glCndVNivqPbs2mmxDNme3TtVpVr1OLcpXrKUtm/bajG2e+cOFSlS1HxOKFGypPbs2mlR2d+zc6dKlCz93DEbRswFT3BuwcHBatmypcqWLauAgADNmDFDZ86cUefOnSXFFEvPnz+v+fPnS5Jee+01dezYUdOmTTO3PfTq1Uvly5dX1qwJP0cmOPmdNGmSOnXqJA8PD02aNOmpc3v06JHgABBbmlRuypfj0Q+O3Nm8VKJgNl29cUdnL15NwsjwIlRv2FRfhI5SrvyFladQMe1Yv1pXw/9W5Toxvxmvmj9N1yLC1aZ3TFvR3h826ovQUWrSoZfyFCqq6/9Wjd3c3JUqTczV+JXrvK6ta5Zr2axQVX31TV0KO6f1y+arav0mcQeBJNOsZWu1a9VMc2ZOV83adfTL0aNauXyZ3h/2qFVsysQJuvT33xo59iNJ0tLFC/XSS1mV+9+bCx06eEBffjFXTd9pbt6m03td9fFHY5UmTVoFVqqs+/fu69ixn3Xjxg21aNUmzlgqBgYpT958Gjp4gHoE99ON69c1cfzHatS4iXmlhzWrVylFihQq5Ocnk8lF27dt0eJFC9S9d+wVBpB4mrdqraGDB8qvaDGVKFlKK5Yv1cWwMDVu0lRS7GOmcZO3tfSrRZrw8Yd6vXETHTl8SKtXrtCYjz4x7/Pt5q3UqW1LzZszU1Wr1dDWLZu0Z88uzZ63wDznzp3bOnvmjPnr8+fP6bdfjytDhgx6KUtW3b1zR3NmTtfLVavJ28dH169d07IlX+nS3xf1Sq3aL+jTQXLVtGlTRUREaOTIkQoLC1OxYsW0bt06812Fw8LCdOax46tNmza6efOmpkyZoj59+sjT01PVq1fXRx99ZNXrmowE1onz5Mmjffv2ycvL66l3bzOZTOYVH6yRqnQ3q7f5r6rsX0AbZvWMNf7lN7vVadiCOLZwLmu+Gp7UISS6betWaOPKhbpxJUJZcuXVm+17qMC/N7mYP3G0Ii5dVO8xUyRJn77fTb//HPvq64rV66pVz0drb/7x689aPnuizp0+KU8vbwW+Ul+13mhh0QrxX1Q+z9PbRZKj7du2aMrET3X2zF/Kmi27mrdsrdcfW7Js+AeDdOHCec2YE1MNWbxogVYsW6IL58/LNYWrsmfPodcbN9EbTZqaK/2S9N3aNZr/xRydPnVSqVKlVv4CBfROi1YW1cInXQy7oA/HjNTen/bIw91dtevVV68+/eXmFnOTizWrV+mLubMUduGCXF1dlDNXbr3TohX9vklg2eJFmj9vtsIvX1a+/AUU3G+gypSNaYF68piRYm5yMWHch+abXLRq1yHWTS6+37Be06ZM1Plz55Q9Rw516d7TfNGkJO3b+5M6t28dK5b6DRpp+OgQRUZG6oOBffXz0SO6dvWqMnh6qkjR4mrfqbOKFiueSJ/Ei5HO3aYFsxJFlk5fJ3UIkqSwGY2TOoQESXDym9hIfpFQzpD8wn4cMfkFkPyR/MbmKMmvzT2/Uky/zunTp5UvX75YDcgAAABIfI54wVtSsunXljt37qh9+/ZKnTq1ihYtau7H6NGjhz788EO7BggAAADYi03J76BBg3T48GFt3bpVHh4e5vFXXnlFS5YssVtwAAAAgD3Z1KuwatUqLVmyRBUrVrQotRcpUkSnTp2yW3AAAAB4BroerGJT5ffy5cvK/Ni93h+6ffs2fScAAABItmxKfsuVK6e1a9eav36Y8M6cOVMBAQH2iQwAAACwM5vaHkJCQlSnTh0dO3ZMDx480MSJE/XLL79o165d2rZtm71jBAAAQDz4q7t1bKr8BgYGaufOnbpz547y5cunDRs2yNfXV7t27ZK/v7+9YwQAAADswqbKb/PmzVW1alW9//77KliwoL1jAgAAABKFTZXftGnTavz48fLz81PWrFn1zjvv6PPPP9evv/5q7/gAAADwFCaTKVk8HIVNye/06dP166+/6vz585owYYIyZMigiRMnqmjRosqSJYu9YwQAAADs4rnuSZwuXTplzJhRGTNmlKenp1KkSKGXXnrJXrEBAADgGRyp6poc2FT5HTBggCpWrChvb2998MEHunfvngYNGqS///5bBw8etHeMAAAAgF3YVPn9+OOP5ePjo2HDhqlhw4by8/Ozd1wAAACA3dmU/B48eFDbtm3T1q1bNX78eLm6uqpKlSqqWrWqqlatSjIMAADwotD1YBWbkt+SJUuqZMmS6tGjhyTp8OHDCg0NVY8ePRQdHa2oqCi7BgkAAADYg80XvB08eFBbt27V1q1btX37dt24cUOlSpVStWrV7BkfAAAAYDc2Jb8ZM2bUrVu3VLJkSVWtWlUdO3bUyy+/rPTp09s7PgAAADwFqz1Yx6bk98svvyTZBQAAgMOxKfmtX7++veMAAACADaj8WsemdX4BAAAAR0TyCwAAAKfxXLc3BgAAQNKi7cE6VH4BAADgNEh+AQAA4DRoewAAAHBgtD1Yh8ovAAAAnAbJLwAAAJwGbQ8AAACOjK4Hq1D5BQAAgNOg8gsAAODAuODNOlR+AQAA4DRIfgEAAOA0aHsAAABwYLQ9WIfKLwAAAJwGyS8AAACcBm0PAAAADoyuB+tQ+QUAAIDToPILAADgwLjgzTpUfgEAAOA0SH4BAADgNGh7AAAAcGB0PViHyi8AAACcBskvAAAAnAZtDwAAAA6M1R6sQ+UXAAAAToPKLwAAgAOj8GsdKr8AAABwGiS/AAAAcBq0PQAAADgwFxf6HqxB5RcAAABOg+QXAAAAToO2BwAAAAfGag/WofILAAAAp0HyCwAAAKdB2wMAAIAD4/bG1qHyCwAAAKdB5RcAAMCBUfi1DpVfAAAAOA2SXwAAADgN2h4AAAAcGBe8WYfKLwAAAJwGyS8AAACcBm0PAAAADoy2B+skm+R3zVfDkzoEOIj67wxP6hDgQK7unZLUIcBBPIgykjoEAC9Askl+AQAAYD0Kv9ah5xcAAABOg+QXAAAAToO2BwAAAAfGBW/WofILAAAAp0HyCwAAAKdB2wMAAIADo+vBOlR+AQAA4DRIfgEAAOA0aHsAAABwYKz2YB0qvwAAAHAaVH4BAAAcGIVf61D5BQAAgNMg+QUAAIDToO0BAADAgXHBm3Wo/AIAAMBpkPwCAADAadD2AAAA4MDoerAOlV8AAAA4DSq/AAAADowL3qxD5RcAAABOg+QXAAAAToO2BwAAAAdG14N1qPwCAADAaZD8AgAAwGnQ9gAAAODAWO3BOlR+AQAA4DRIfgEAAOA0aHsAAABwYHQ9WIfKLwAAAJwGlV8AAAAHxgVv1qHyCwAAAKeR4Mpv6dKlE/ybxYEDB2wOCAAAAEgsCU5+GzVqlIhhAAAAwBZ0PVgnwcnvsGHDEjMOAAAAINHR8wsAAACnYdNqD1FRUfr000+1dOlSnTlzRvfu3bN4/sqVK3YJDgAAAE/Hag/WsanyO2LECE2YMEFvvfWWrl+/ruDgYL3xxhtycXHR8OHD7RwiAAAAYB82Jb8LFy7UzJkz1bdvX6VIkULvvPOOZs2apaFDh2r37t32jhEAAADxMJlMyeLhKGxKfi9evKjixYtLktKmTavr169LkurXr6+1a9faLzoAAADAjmxKfrNnz66wsDBJUv78+bVhwwZJ0t69e+Xu7m6/6AAAAAA7sin5ff3117Vp0yZJUs+ePTVkyBAVKFBArVq1Urt27ewaIAAAAOJnMiWPh6OwabWHDz/80PzvN998U9mzZ9fOnTuVP39+NWjQwG7BAQAAAPZkU/L7pIoVK6pixYr22BUAAACQaBKc/H7zzTeqW7euUqZMqW+++eapc6n+AgAAvBiOtNJCcpDg5LdRo0a6ePGiMmfOrEaNGsU7z2QyKSoqyh6xAQAAAHaV4OQ3Ojo6zn8DAAAg6VD4tY5Nqz0AAAAAjsim5LdHjx6aNGlSrPEpU6aoV69ezxsTAAAAkChsSn6//vprBQUFxRoPDAzU8uXLnzsoAAAAJExS39bYKW5vHBERoQwZMsQaT58+vcLDw587KAAAACAx2JT85s+fX999912s8f/973/KmzfvcwcFAAAAJAabbnIRHBysbt266fLly6pevbokadOmTRo/frxCQ0PtGR8AAACewoE6DpIFm5Lfdu3aKTIyUmPGjNGoUaMkSblz59a0adPUqlUruwYIAAAA2IvVye+DBw+0cOFCvf7663rvvfd0+fJlpUqVSmnTpk2M+AAAAAC7sTr5TZEihd577z0dP35ckuTj42P3oAAAAJAwLvQ9WMWmC94qVKiggwcP2jsWAAAAIFHZ1PPbpUsX9enTR+fOnZO/v7/SpElj8XyJEiXsEhwAAACejsKvdWyq/DZt2lSnT59Wjx49FBQUpFKlSql06dLm/wcAAACeZerUqcqTJ488PDzk7++v7du3P3V+ZGSk3n//feXKlUvu7u7Kly+f5syZY9Vr2lT5PX36tC2bAQAAAJKkJUuWqFevXpo6daqCgoI0ffp01a1bV8eOHVPOnDnj3Oatt97S33//rdmzZyt//vy6dOmSHjx4YNXr2pT85sqVy5bNAAAAYGeOdGvhx02YMEHt27dXhw4dJEmhoaFav369pk2bppCQkFjzv/vuO23btk1//PGHMmXKJClmqV1r2dT2IElffvmlgoKClDVrVv3111/moFevXm3rLgEAAOCgIiMjdePGDYtHZGRknHPv3bun/fv3q1atWhbjtWrV0s6dO+Pc5ptvvlHZsmU1btw4ZcuWTQULFlTfvn119+5dq+K0KfmdNm2agoODVa9ePV27dk1RUVGSJE9PT+7wBgAA4IRCQkKUIUMGi0dcFVxJCg8PV1RUlHx9fS3GfX19dfHixTi3+eOPP/Tjjz/q559/1sqVKxUaGqrly5era9euVsVpU/I7efJkzZw5U++//75cXV3N42XLltXRo0dt2aVT2LZuhYZ0fFM93qymkOB2OvnLoXjnHty1VZOG9lT/lq8q+O2a+rh/Jx07sCfWvDu3bmrx5+M1sE0D9XizmkZ0baaf98X9GxP+e4LK5NPy0Hf1x4Yxuntwil6rykorzsAwDE37bLJeqVpJ5cuUUPs2LXXy5O/P3O77Dev1+mv1VLZUMb3+Wj1t+n6jxfOzZ05Xs7caK6BcaVWtHKBe3bvoz9N/xNrPH6dOqUfXzgqq4K+AcqXV4p23FHbhgt3eH2xnGIamT52s2jUqK7BcSXVq11KnEnBsbNq4Xm82elUV/YvrzUavavOmjbHmLF28SK/VqaGAsiXUvOkbOrh/n8Xzwz4YKP8ShS0erZs3jTfO7u91lH+Jwtqy+Xvb3izMXEzJ4zFo0CBdv37d4jFo0KCnxv5ky4ZhGPG2cURHR8tkMmnhwoUqX7686tWrpwkTJmjevHlWVX9tSn5Pnz4d56oO7u7uun37ti27/M/bt/17LZ89UXWatNKgT+cqf5ES+mxkX125HPdvNyd/OaTCpcqry9BPNHDCHBUsXkbTxvTX2T9OmOc8uH9fk4b1UsSlMHUcMFrDpn6l5l0HyNOLG484izSp3HX0xHn1/nBpUoeCF2ju7Jn68ou5Gvj+UC1cslxe3t7q3KGtbt++Fe82hw8dVP++vVW/QUMtW7Fa9Rs0VP8+vXTkyGHznH17f1LTd5rry6+WavrMuXoQFaXOHdvrzp075jlnz5xRm5bNlCdPXs2a96WWrfhGnTp3kZu7e6K+ZyTMF3NnaeGX8zRg0BDNX7RMXt4+6vJuu6ceG0cOH9Sg/sGqV7+Bvlq+WvXqN9DAfr119LFjY8N36zR+XIjadeysRUtXqnSZsurepZPCwix/6QkMqqz1m7ebH5OmTo/zNRct+MJh+1QRP3d3d6VPn97i4R7PucHb21uurq6xqryXLl2KVQ1+KEuWLMqWLZsyZMhgHvPz85NhGDp37lyC47Qp+c2TJ48OHToUa/x///ufihQpYssu//M2r16iwFfqK6hWA2XJkVtNOvSSp3dm/fC/lXHOb9Khl2q90Vy5C/gpc9YcatiyszJnya6jP/1onrPz+zW6c+uGOg/+UPn8Ssgr80vKX6Sksucp8KLeFpLYhh3HNGLqGq3efPjZk/GfYBiGFn45Xx06ddYrNWupQIGCGj32I/3zzz9at3ZNvNst+PILVQwIVPuO7ypP3nxq3/Fdla9QUQvnf2GeM23GbDV8/Q3lz19AhQoX1sjRIQoLu6Djx34xz5k86VNVevll9e7bX35+RZQ9Rw69XKWqvLy8EvV949kMw9CiBfPVrmNnVX+llvIXKKgRoz/UP//8o+/WxX9sLFowXxUqBqpdh3eVJ09etevwrsqXr6ivFjw6NhbMn6eGrzfW642bKE/efOo7YLB8X3pJy5d+ZbGvlG5u8vb2MT8yZPCM9XonfvtVC+fP09CRY+z23p2dyWRKFg9ruLm5yd/fXxs3Wv6VYePGjQoMDIxzm6CgIF24cEG3bj36Ze7EiRNycXFR9uzZE/zaNiW//fr1U9euXbVkyRIZhqGffvpJY8aM0eDBg9WvXz9bdvmf9uD+fZ059Zv8SpW3GPcrVV5//PpzgvYRHR2tf+7eVep06c1jR/f+qDyFimnx9PEa0Kq+RnVvoe+WfaHof3uwAfz3nD93TuHhlxUQVMk85ubmJv+y5XT4KXfePHLokAICK1mMBQZV1uFD8W9z6+ZNSVL6f6ss0dHR2r5tq3Llyq3OHdurauUANX+7iTZv4s/WycH58+cUEX5ZFQOCzGMxCUa5p36fjxw+pIqBQRZjAUGVdPjwIUnS/fv39OvxX2LNqRgQpCNP7Hf/vp/0SpVAvf5abY0aPkRXIiIsnr97964GD+ij/oOHyNubv1I6u+DgYM2aNUtz5szR8ePH1bt3b505c0adO3eWFNNG0apVK/P8Zs2aycvLS23bttWxY8f0ww8/qF+/fmrXrp1SpUqV4Ne1aamztm3b6sGDB+rfv7/u3LmjZs2aKVu2bJo4caLefvttW3b5n3brxjVFR0cpnWcmi/H0nhl142pEPFtZ2rTqK92LvCv/oBrmsfCLFxRx6YDKVamlrkM/0aUL57RkxnhFR0Wp3tvt7PoeACQP4eGXJSlWpdXLy1sXntJ3Gx4eHsc2Xub9PckwDH0yLkSly/irQIGCkqQrERG6c+eO5syeqW7de6lXcF/t+HG7gnt206y581W2XPk494UXIyKeYyOTl1es9gTL7cKVKdMT22TyMu/v2tWrioqKivP4iQgPN38dVOllvVKrjrJkyaoL589p2meT1LlDGy1Y8rXc3NwkSRM+DlGJkqVVtVoNAU2bNlVERIRGjhypsLAwFStWTOvWrTMvqRsWFqYzZ86Y56dNm1YbN25U9+7dVbZsWXl5eemtt97S6NGjrXpdm5JfSerYsaM6duyo8PBwRUdHK3PmzAneNjIyMtbSF/fuRcrN7b/dMxa7qTtha/Pt/WGj1i6eo86DP1Q6z4yPbW8oXYaMat6lv1xcXZUzf2FdvxqujSsXkfwC/xFr13yjUcOHmb+eMi2mhzLui0Sevq9Y2yj+C0tCRo/U7ydOaN6Xi8xj0Ua0JKlatRpq2bqNJKmwn58OHzqgZUsWk/y+YOvWfquxIx8dGxM/+zzmH3H9rNHTD47Yx0bssbh+hj3+WrXq1DP/O3+BgvIrWkz1a9fQjz9sVfVXamnbls3a+9MeLVq64llvDVZy5PbpLl26qEuXLnE+N2/evFhjhQsXjtUqYS2bk9+HvL29rd4mJCREI0aMsBhr2bWfWnfr/7zhJEtp03vKxcU1VpX35vWrsarBT9q3/XstmByiDgNGq3CpchbPpc/oJVfXFHJ5bMWNl7Ln0o2rEXpw/75SpExpvzcBIElUrVZdxYuXNH997/49STGVXB+fR0WHK1ci5OUV//nY29tb4Y9V6STpSsSVOLcJGTNKW7du1pwvFsj3pZfM4xk9MypFihTKmy+fxfw8efPp0IH91r0xPLcqVaupePFHK7zcuxdzbEQ8cWxcvRKhTE/pyfby9lZEhOWxEbNNzLHhmTGjXF1dYx8/VyKe2uvt45NZWbJm1ZkzMfcC2PvTbp07e0ZVgyx/Seof3EOly/hrxpwvn/Z2Abuxqec3IiJCXbt2VZEiReTt7a1MmTJZPJ4lrqUw3unU05ZQHEKKlCmVM18hHT+812L810N7lbdwsXi32/vDRn05aYza9hmu4mVjN3/n8yuuyxfPKTo62jz294WzypDRi8QX+I9IkyatcubKZX7ky5df3t4+2r1zh3nO/Xv3tH/fXpWMYxWeh0qUKqXdu3ZYjO3a+aNKlnq0jWEYGjt6pDZ9v0Ez53yh7NlzWMxP6eamosWK688/LW9x/9dffypL1mzP8zZhgzRp0ipHzlzmR958+eXl7aM9ux4td3n//j3t37/X4vv8pBIlS1lsI0m7d+5QyZKlJEkpU7qpsF/RWHP27N6pEk/Z77VrV/X3xTBzb2+b9h21ePlqLVq60vyQpOB+AzVsZNxrwQKJwabKb4sWLXTq1Cm1b99evr6+Vl/h5+7uHmvpCze3e7aE4jCqN2yqL0JHKVf+wspTqJh2rF+tq+F/q3Kd1yVJq+ZP07WIcLXpPURSTOL7RegoNenQS3kKFdX1f6vGbm7uSpUmrSSpcp3XtXXNci2bFaqqr76pS2HntH7ZfFWt3yRp3iReuDSp3JQvx6OLRnJn81KJgtl09cYdnb14NQkjQ2IxmUxq3rKVZs+crpy5citnrlyaPWO6PDw8VO/V+uZ57w/qr8yZfdWzdx9JUvMWrdSudQvNmTVD1arX0JbNm7Rn9y7NfaytYeyoEfrfujUKnTxVaVKnUfjlmJ7PtOnSycPDQ5LUum179e/TW/7+5VSufAXt+HG7fti6RbPmzn+BnwLiYjKZ1KxFK82ZPV05cuVSzpy5NGdWzLFRp96jY2Po4AHy8c2s7j1jjo13mrdUx7YtNW/OTFWpVkPbtmzSnj27NHveQvM2LVq10ZDBA1SkaDGVKFlKK5Yv1cWwML3ZJOY6nzt3bmv61CmqUbOWvL19dOHCeX026VN5emZUtRqvSJJ5BYgnvZQlq7JZcaU+YntWWwss2ZT8/vjjj/rxxx9VsmTJZ0+GJKls5Vd0++YNrVsyVzeuRChLrrzqMvQTeWWO+ZPijasRuhr+t3n+j+tXKzoqSkumj9eS6ePN4xWr11Wrnh9IkjL5+Kr7iFAtnz1RY3q2lqeXt6q91kS13mjxYt8ckkyZIrm0Ydajv5qM69tYkvTlN7vVadiCpAoLiaxt+46KjIzU2FEjdOPGdRUvUVLTZs5Rmn9/MZaki2FhcjE9+uNeqdJl9NHHEzRlcqg+mzxJOXLm0EeffKoSJR6dx5cuiVm2qn2blhavN3J0iBq+/oYkqcYrNfXBsOGaM3OGPgoZrdy582h86CSV8S+bmG8ZCdS6bQdF/vOPPhwzUjdvXFex4iX02eezLY+NixdkcnmULJUsVUZjPxqvqVMmatqUScqeI4c+HDdBxR87NmrVibmj68zpnyn88mXly19Akz6bbq74u7i46uTJE1r77WrdvHlT3j4+KluuvEI+/tTitYHkwGQYhmHtRuXKldPkyZNVsWJFuwWy6dfwZ08CJNV/Z3hShwAHcnXvlKQOAQ7iQZTVPw7hxNK6J59qa/3pe5896QVY8265Z09KBmzq+Z06daref/99bdu2TREREbpx44bFAwAAAC9GUt/W+OHDUdjU9uDp6anr16+revXqFuMP78ccxU0WAAAAkAzZlPw2b95cbm5uWrRokU0XvAEAAMA+yMOsY1Py+/PPP+vgwYMqVKiQveMBAAAAEo1NPb9ly5bV2bNn7R0LAAAAkKhsqvx2795dPXv2VL9+/VS8eHGlfOKGCiVKlIhnSwAAANgTXQ/WsSn5bdq0qSSpXbt25jGTycQFbwAAAEjWbEp+T58+/exJAAAAQDJjU/KbK1cue8cBAAAAG7jQ92AVm5JfSTp//rx27NihS5cuKTo62uK5Hj16PHdgAAAAgL3ZlPzOnTtXnTt3lpubm7y8vCzWlzOZTCS/AAAALwiFX+vYlPwOHTpUQ4cO1aBBg+TiYtNqaQAAAMALZ1PmeufOHb399tskvgAAAHAoNmWv7du317Jly+wdCwAAAKxkMpmSxcNR2NT2EBISovr16+u7776L8yYXEyZMsEtwAAAAgD3ZlPyOHTtW69evV6FChSQp1gVvAAAAQHJkU/I7YcIEzZkzR23atLFzOAAAALAGdUfr2NTz6+7urqCgIHvHAgAAACQqm5Lfnj17avLkyfaOBQAAAFZyMZmSxcNR2NT28NNPP2nz5s1as2aNihYtGuuCtxUrVtglOAAAAMCebEp+PT099cYbb9g7FgAAACBR2Xx7YwAAACQ9x2k4SB5sSn4funz5sn777TeZTCYVLFhQPj4+9ooLAAAAsDubLni7ffu22rVrpyxZsujll19W5cqVlTVrVrVv31537tyxd4wAAACAXdiU/AYHB2vbtm369ttvde3aNV27dk2rV6/Wtm3b1KdPH3vHCAAAgHgk9W2NneL2xl9//bWWL1+uqlWrmsfq1aunVKlS6a233tK0adPsFR8AAABgNzZVfu/cuSNfX99Y45kzZ6btAQAAAMmWTclvQECAhg0bpn/++cc8dvfuXY0YMUIBAQF2Cw4AAABP52JKHg9HYVPbQ2hoqOrWravs2bOrZMmSMplMOnTokNzd3bVhwwZ7xwgAAADYhU3Jb/HixfX7779rwYIF+vXXX2UYht5++201b95cqVKlsneMAAAAiIcjXWyWHNiU/IaEhMjX11cdO3a0GJ8zZ44uX76sAQMG2CU4AAAAwJ5s6vmdPn26ChcuHGu8aNGi+vzzz587KAAAACAx2FT5vXjxorJkyRJr3MfHR2FhYc8dFAAAABKGrgfr2FT5zZEjh3bs2BFrfMeOHcqaNetzBwUAAAAkBpsqvx06dFCvXr10//59Va9eXZK0adMm9e/fnzu8AQAAINmyKfnt37+/rly5oi5duujevXuSJA8PDw0YMECDBg2ya4AAAACIH6s9WMem5NdkMumjjz7SkCFDdPz4caVKlUoFChSQu7u7veMDAAAA7Mam5PehtGnTqly5cvaKBQAAAFZypLurJQc2XfAGAAAAOCKSXwAAADiN52p7AAAAQNLigjfrUPkFAACA0yD5BQAAgNOg7QEAAMCB0fRgHSq/AAAAcBokvwAAAHAatD0AAAA4MBdWe7AKlV8AAAA4DSq/AAAADozCr3Wo/AIAAMBpkPwCAADAadD2AAAA4MC4vbF1qPwCAADAaZD8AgAAwGnQ9gAAAODA6HqwDpVfAAAAOA0qvwAAAA6MO7xZh8ovAAAAnAbJLwAAAJwGbQ8AAAAOjK4H61D5BQAAgNMg+QUAAIDToO0BAADAgXF7Y+tQ+QUAAIDToPILAADgwKhkWifZJL/l82RK6hDgIK7unZLUIcCBZCzXLalDgIO48hPnFsAZ8MsCAAAAnEayqfwCAADAelzwZh0qvwAAAHAaJL8AAABwGrQ9AAAAODAXuh6sQuUXAAAAToPkFwAAAE6DtgcAAAAHRtuDdaj8AgAAwGlQ+QUAAHBgrPNrHSq/AAAAcBokvwAAAHAatD0AAAA4MC54sw6VXwAAADgNkl8AAAA4DdoeAAAAHBiLPViHyi8AAACcBpVfAAAAB+ZC6dcqVH4BAADgNEh+AQAA4DRoewAAAHBgVDKtw+cFAAAAp0HyCwAAAKdB2wMAAIADY7EH61D5BQAAgNMg+QUAAIDToO0BAADAgXGTC+tQ+QUAAIDToPILAADgwCj8WofKLwAAAJwGyS8AAACcBm0PAAAADsyFtgerUPkFAACA0yD5BQAAgNOg7QEAAMCBsc6vdaj8AgAAwGlQ+QUAAHBgFH6tQ+UXAAAAToPkFwAAAE6DtgcAAAAHxjq/1qHyCwAAAKdB8gsAAACnYXXyO3XqVL3yyit66623tHnzZovnwsPDlTdvXrsFBwAAgKczJZP/OQqrkt9JkyapX79+Kly4sNzd3VWvXj2FhISYn4+KitJff/1l9yABAAAAe7Dqgrfp06dr5syZatasmSSpS5cuatSoke7evauRI0cmSoAAAACAvViV/J4+fVqBgYHmrwMCArR582bVqFFD9+/fV69evewdHwAAAJ6C1R6sY1Xy6+3trbNnzyp37tzmsaJFi2rz5s2qXr26zp8/b+/4AAAAALuxque3UqVK+vrrr2ONFylSRJs2bdJ3331nt8AAAADwbC6m5PFwFFZVfgcOHKj9+/fH+VzRokW1ZcsWLV++3C6BAQAAAPZmVfJbokQJlShRIt7nixYtqqJFiz53UAAAAEBisPn2xv/884+OHDmiS5cuKTo62jxuMpn02muv2SU4AAAAPJ3J5EA9B8mATcnvd999p1atWik8PDzWcyaTSVFRUc8dGAAAAGBvNt3euFu3bmrSpInCwsIUHR1t8SDxBQAAQHJlU/J76dIlBQcHy9fX197xAAAAwApJvcrD86z2MHXqVOXJk0ceHh7y9/fX9u3bE7Tdjh07lCJFCpUqVcrq17Qp+X3zzTe1detWWzYFAAAAtGTJEvXq1Uvvv/++Dh48qMqVK6tu3bo6c+bMU7e7fv26WrVqpRo1atj0uibDMAxrN7pz546aNGkiHx8fFS9eXClTprR4vkePHlYHcjMy+tmTAEkpXW36nQ1OKmO5bkkdAhzElZ+mJHUIcCCpUj57zosy4Yc/kjoESVLwy3mtml+hQgWVKVNG06ZNM4/5+fmpUaNGCgkJiXe7t99+WwUKFJCrq6tWrVqlQ4cOWfW6Nl3wtmjRIq1fv16pUqXS1q1bLa4yNJlMNiW/AAAAcFyRkZGKjIy0GHN3d5e7u3usuffu3dP+/fs1cOBAi/FatWpp586d8b7G3LlzderUKS1YsECjR4+2KU6bSmgffPCBRo4cqevXr+vPP//U6dOnzY8//kgev30AAADgxQkJCVGGDBksHvFVcMPDwxUVFRXr+jFfX19dvHgxzm1+//13DRw4UAsXLlSKFDav1mtb5ffevXtq2rSpXFz48zMAAEBSckkm6/wOGjRIwcHBFmNxVX0f9+QaxYZhxLlucVRUlJo1a6YRI0aoYMGCzxWnTdlr69attWTJkud6YQAAAPx3uLu7K3369BaP+JJfb29vubq6xqryXrp0Kc7VxG7evKl9+/apW7duSpEihVKkSKGRI0fq8OHDSpEihTZv3pzgOG2q/EZFRWncuHFav369SpQoEeuCtwkTJtiyWwAAADgBNzc3+fv7a+PGjXr99dfN4xs3blTDhg1jzU+fPr2OHj1qMTZ16lRt3rxZy5cvV548eRL82jYlv0ePHlXp0qUlST///LMtuwAAAIAd2LrGblILDg5Wy5YtVbZsWQUEBGjGjBk6c+aMOnfuLCmmjeL8+fOaP3++XFxcVKxYMYvtM2fOLA8Pj1jjz2JT8rtlyxZbNgMAAAAkSU2bNlVERIRGjhypsLAwFStWTOvWrVOuXLkkSWFhYc9c89cWNvX8btq0Kd7npkxxnnUSL/39t4YM6q8alSsqqHxpNWvyuo4f+yXe+YcO7Fe7Vs1i5pcrpcYN6mnhl/PsEsvFsAvq3e09VSpfRjVeDtDHH47R/fv3zM9fOH9eZUv4xXrs/DFhd1JB4jIMQ9M+m6xXqlZS+TIl1L5NS508+fszt/t+w3q9/lo9lS1VTK+/Vk+bvt9o8fzsmdPV7K3GCihXWlUrB6hX9y7683TsFVn+OHVKPbp2VlAFfwWUK60W77ylsAsX7Pb+kDwElcmn5aHv6o8NY3T34BS9VrVEUoeERPbw3FKzWiVV8Lfi3LJxvd5oUE/lShfTGw3qaXNc55amjRVYvrSqvRygXj1in1s2bdyg9zq1V9VKFVSqWCH9+utxu743PGIyJY+HLbp06aI///xTkZGR2r9/v15++WXzc/PmzXvqTdWGDx9u9Rq/ko3Jb+PGjbV3795Y46GhoRo8eLAtu3Q4N25cV/vWzZQiRQpNnDpDy1auUa8+/ZUuXbp4t/FIlUpvvd1cM+Z+qWWr1qp9p86aNnmSVixf+lyxREVFqWfXzrp7965mfbFAY8eN1+bvN+jTT8bFmjt1xhx9t/kH86NchQrP9dqwj7mzZ+rLL+Zq4PtDtXDJcnl5e6tzh7a6fftWvNscPnRQ/fv2Vv0GDbVsxWrVb9BQ/fv00pEjh81z9u39SU3faa4vv1qq6TPn6kFUlDp3bK87d+6Y55w9c0ZtWjZTnjx5NWvel1q24ht16txFbs+4QheOJ00qdx09cV69P3y+cw4cx7w5M7Vg/lwNHDxUCxcvl7e3t97r+Oxzy4C+vfXqaw219OvVevW1hurft5eOPnZu2b8v5twyf9FSfT5jrqIeROm9Tu1197Fzy927d1SqdGn16NU3Ud8jYC2b7vA2d+5c9e/fX9u2bVORIkUkSZ988olGjRqlNWvWqHLlylYH4mh3eJscOl6HDx7UrC8WPNd++vXuLo9UqTRq7KNE9ZtVKzR/7mxdOH9OWbJm09vNWqjJ283i3ceO7T+od/f3tHbDFvlkzixJWv+/tRoxZLA2bN2htGnT6sL582pQ9xUtXLpChQr7PVfMSe2/doc3wzD0StXKat6yldp16CQpZjnB6i8HqmdwXzV56+04t+vXp5du37qlqdNnmcfe69Re6dNn0EefxH3R6ZUrV1StcoDmfLFA/mXLSZL69+2tFClSaOyHH9v5nSUP3OEtbncPTtFbvWfo261HkjqUZOO/doc3wzBUs1rMuaVt+8fOLVUC1at3X70Zz7mlf59eun37lj77/NG5pcu7MeeWDz+O/9xS/eUAzZ736Nzy0Pnz5/Rq7RpavHyVCjv4z5/HJac7vE3ecTqpQ5AkdQ9K+EVnScmmLKJt27YaMGCAatWqpT///FMfffSRRo0apf/97382Jb6O6IetW+RXtKgG9OmlmlWC1OytN7TSygrur8eP6cihQ/L3f3SiWLl8qaZODlWX7j21bNVade3RS59/NklrVq+Kdz9HjxxSvvwFzImvJAUEVdK9e/f06xNtGME9uqpmlSC1a9VM329Yb1W8SBznz51TePhlBQRVMo+5ubnJv2w5HT54MN7tjhw6pIDAShZjgUGVdfhQ/NvcunlTkpQ+QwZJUnR0tLZv26pcuXKrc8f2qlo5QM3fbqLNm75/nrcEIBkwn1sCLc8tZcuW06GnnCeOHD6kik+cWwKedW65FXNuyfDvuQUvlotMyeLhKGy+PUbfvn0VERGhsmXLKioqShs2bFCFBP4JPa7b391TymcuhJycnD93Vl8vXazmLduobYdO+uXno/rko7FK6eam+g0aPXXbeq9U1dWrVxQVFaVO73VVo8ZNzM/NmvG5evUZoOqv1JIkZcueXX/8cUorli9R/YZx7zciPFyZvLwsxtKnz6CUKVMqPDxckpQ6dWr17jdAJUuVkYuLi37YulmD+wfr3r0Q1avfwPYPAs8tPPyyJMnrie+hl5e3Ljyl7zY8PDyObbzM+3uSYRj6ZFyISpfxV4ECMQuEX4mI0J07dzRn9kx1695LvYL7aseP2xXcs5tmzZ2vsuXKP89bA5CEHp4Lnvz5kMnL+6k9/bacW8b/e27JX+D5bj4AvAgJTn4nTZoUayxLlixKnTq1Xn75Ze3Zs0d79uyRJPXo0eOp+woJCdGIESMsxga+P1SDhwxLaDhJLjraUJGiRdW1Z29JUmG/Ivrj1El9vXTxM5PfmfMW6O6dOzp65JCmTJyg7DlyqU69V3X1yhX9fTFMo4Z/oDEjhprnR0U9UNq0Mb3EPd7rpIMH9kuSsmTNoqUr10iKfYcUSTKMRw3onhkzqnnLNubnihQtphs3bmj+3Nkkvy/Y2jXfaNTwR8f6lGnTJcV3l5un7yvWNor7zjiSFDJ6pH4/cULzvlxkHos2YtqNqlWroZat20iSCvv56fChA1q2ZDHJL+BA1q75RqNHPDq3TJ5qx3NLPHfdkqSQMSN14sQJzZu/KM7ngeQmwcnvp59+Gue4q6urduzYoR07dkiK+Q/mWclvXLe/u6dk1DyTAN4+3sqTN5/FWJ48ebX5+w3P3DZb9uySpPwFCyriSoRmTJuiOvVeNSciHwwbqWLFLa/CdnFxjXlu+ChFRv4jSeb7Wnt5e+vno5Z9ezduXNeDB/fl5eUdbxzFS5TU6hXLnxkv7KtqteoqXryk+et7/67KER4eLh+fR60rV65EPPX75+3tba7sm7eJuBLnNiFjRmnr1s2a88UC+b70knk8o2dGpUiRQnnzPXEs582nQ//+kgXAMVStVl3FSzx2brkXc26JeOLccvVKhDJZe265Eve55cOxo7RtS+xzC16sZHJ3Y4eR4OT39Gn7NVO7u7vHanFwtAveSpYqo7/+/NNi7K+//lSWLFmt25FhmJck8/LyVubMvjp/7qzqvvpanNMzx3HLv+IlSmnOzOkKv3xJ3v+e4Hbv3CE3NzcVLlI03pf+7fhxeXv7WBcvnluaNGmVJk1a89eGYcjb20e7d+6Qn1/MBaT3793T/n171TM4/qukS5Qqpd27dpgrtpK0a+ePKlmqtMW+Q8aM0uZNGzV73pfKnj2HxT5SurmpaLHi+vNPy/++//rrT2XJmu153iaAFyy+c8uuXTtU+OG55f497du3V716P+XcUvLfc0urNuax3XGcWz4cG3NumTX3S2V74twCJGdW9fwGBgaqUaNGatCggQoXLpxYMTmEZi1bq12rZpozc7pq1q6jX44e1crly/T+sEftHFMmTtClv//WyLEfSZKWLl6ol17Kqtz/3oLv0MED+vKLuWr6TnPzNp3e66qPPxqrNGnSKrBSZd2/d1/Hjv2sGzduqMVjJ6LHVQwMUp68+TR08AD1CO6nG9eva+L4j9WocROlTRtzIlyzepVSpEihQn5+MplctH3bFi1etEDdewfHuU+8OCaTSc1bttLsmdOVM1du5cyVS7NnTJeHh4fqvVrfPO/9Qf2VObOvevbuI0lq3qKV2rVuoTmzZqha9RrasnmT9uzepbmPtTWMHTVC/1u3RqGTpypN6jQKvxzTs5c2XTp5eHhIklq3ba/+fXrL37+cypWvoB0/btcPW7do1tz5L/BTwIuQJpWb8uV49Atv7mxeKlEwm67euKOzF68mYWRIDI+fW3LljDm3zJo5Xak8PFT3sXPLB/+eW3r8e25p1qKV2rdpobmzZ6hqtRrauuXfc8tjbQ1jR/97bpk0VWnSpDH3A6dN++jccv36NYWFhenypUuSpL/+LaJ5e3tTeEGSsmqps3nz5unbb7/Vhg0blCVLFjVs2FANGjRQpUqV4u0FSihHq/xK0vZtWzRl4qc6e+YvZc2WXc1bttbrb75lfn74B4N04cJ5zZgTk0QsXrRAK5Yt0YXz5+WawlXZs+fQ642b6I0mTeXi8mjhje/WrtH8L+bo9KmTSpUqtfIXKKB3WrRStRo1443lYtgFfThmpPb+tEce7u6qXa++evXpLzc3N0kxye8Xc2cp7MIFubq6KGeu3HqnRSuH7Pf9ry11JsVUUT6fOkXLly7RjRvXVbxESQ36YKj5wjRJat+mpbJmzaZRYz80j21c/52mTA7VubPnlCNnDnXr0Vuv1Kxlfr5k0UJxvt7I0SFq+Pob5q9XrliuOTNn6O+/Lyp37jx6r1t3Vav+SiK80xePpc4eqexfQBtm9Yw1/uU3u9Vp2PMt2/hf8F9b6kx6dG75etlj55b3h1pcmNa+TUtlzZZNo8Y8dm7Z8J0+e3huyRFzbqnx2LmlVLG4zy0jRoeoYaOYc8vqVSs07INBsea8+143vde1u73eYpJJTkudfb7rz6QOQZLUOSB3UoeQIDat8xsZGalNmzZp9erV+vbbb3X//n29+uqratiwoWrXrq3UqVNbHYgjJr9IGv/F5BeJh+QXCfVfTH6ReEh+Y3OU5NemLMLd3V316tXT9OnTdeHCBa1Zs0bZsmXT0KFD5e3trfr165svgAMAAEDicTGZksXDUdilhFahQgWNGTNGR48e1dGjR1WjRg2FhYXZY9cAAACA3dh8kwtJ2rdvn44fPy6TySQ/Pz/5+/srX7586t27t73iAwAAAOzGpuT33Llzeuedd7Rjxw55enpKkq5du6bAwEB99dVXypGDJU8AAABeBAfqOEgWbGp7aNeune7fv6/jx4/rypUrunLlio4fPy7DMNS+fXt7xwgAAADYhU2V3+3bt2vnzp0qVOjRUieFChXS5MmTFRQUZLfgAAAAAHuyKfnNmTOn7t+/H2v8wYMHypaNu0IBAAC8KI600kJyYFPbw7hx49S9e3ft27dPD5cJ3rdvn3r27KlPPvnErgECAAAA9mLTTS4yZsyoO3fu6MGDB0qRIqZ4/PDfadKksZh75cqVBO2Tm1wgobjJBazBTS6QUNzkAtZITje5mLP3TFKHIElqVy5nUoeQIDa1PYSGhto5DAAAACDx2ZT8tm7d2t5xAAAAAInO5ptcREVFadWqVeabXBQpUkQNGjSQq6urPeMDAADAU9AMaB2bkt+TJ0+qXr16On/+vAoVKiTDMHTixAnlyJFDa9euVb58+ewdJwAAAPDcbPploUePHsqXL5/Onj2rAwcO6ODBgzpz5ozy5MmjHj162DtGAAAAwC5sqvxu27ZNu3fvVqZMmcxjXl5e+vDDD7nJBQAAwAtkYp1fq9hU+XV3d9fNmzdjjd+6dUtubm7PHRQAAACQGGxKfuvXr69OnTppz549MgxDhmFo9+7d6ty5sxo0aGDvGAEAAAC7sCn5nTRpkvLly6eAgAB5eHjIw8NDgYGByp8/vyZOnGjvGAEAABAPUzJ5OAqben49PT21evVqnTx5UseOHZMkFSlSRPnz57drcAAAAIA92bzO7+zZs/Xpp5/q999/lyQVKFBAvXr1UocOHewWHAAAAJ7OhQverGJT8jtkyBB9+umn6t69uwICAiRJu3btUu/evfXnn39q9OjRdg0SAAAAsAebkt9p06Zp5syZeuedd8xjDRo0UIkSJdS9e3eSXwAAACRLNiW/UVFRKlu2bKxxf39/PXjw4LmDAgAAQMLQ9GAdm1Z7aNGihaZNmxZrfMaMGWrevPlzBwUAAAAkhue64G3Dhg2qWLGiJGn37t06e/asWrVqpeDgYPO8CRMmPH+UAAAAgB3YlPz+/PPPKlOmjCTp1KlTkiQfHx/5+Pjo559/Ns/jdnsAAACJi3TLOjYlv1u2bLF3HAAAAECis7ntAQAAAEmPv7Rbx6YL3gAAAABHRPILAAAAp0HbAwAAgAOjkmkdPi8AAAA4DZJfAAAAOA3aHgAAABwYqz1Yh8ovAAAAnAaVXwAAAAdG3dc6VH4BAADgNEh+AQAA4DRoewAAAHBgXPBmHSq/AAAAcBokvwAAAHAatD0AAAA4MCqZ1uHzAgAAgNMg+QUAAIDToO0BAADAgbHag3Wo/AIAAMBpUPkFAABwYNR9rUPlFwAAAE6D5BcAAABOg7YHAAAAB8b1btah8gsAAACnQfILAAAAp0HbAwAAgANzYb0Hq1D5BQAAgNOg8gsAAODAuODNOlR+AQAA4DRIfgEAAOA0aHsAAABwYCYueLMKlV8AAAA4DZJfAAAAOA3aHgAAABwYqz1Yh8ovAAAAnAbJLwAAAJwGbQ8AAAAOjNsbW4fkFw7nQZSR1CHAgVz5aUpShwAHkal8t6QOAQ7k7kHOLY6K5BcAAMCBccGbdej5BQAAgNMg+QUAAIDToO0BAADAgdH2YB0qvwAAAHAaJL8AAABwGrQ9AAAAODAT6/xahcovAAAAnAaVXwAAAAfmQuHXKlR+AQAA4DRIfgEAAOA0aHsAAABwYFzwZh0qvwAAAHAaJL8AAABwGrQ9AAAAODBub2wdKr8AAABwGiS/AAAAcBq0PQAAADgwVnuwDpVfAAAAOA0qvwAAAA6M2xtbh8ovAAAAnAbJLwAAAJwGbQ8AAAAOjAverEPlFwAAAE6D5BcAAABOg7YHAAAAB8btja1D5RcAAABOg8ovAACAA6Pwax0qvwAAAHAaJL8AAABwGrQ9AAAAODAXrnizCpVfAAAAOA2SXwAAADgN2h4AAAAcGE0P1qHyCwAAAKdB5RcAAMCRUfq1CpVfAAAAOA2SXwAAADgN2h4AAAAcmIm+B6tQ+QUAAIDTIPkFAACA06DtAQAAwIFxd2PrUPkFAACA0yD5BQAAgNOg7QEAAMCB0fVgHSq/AAAAcBpUfgEAABwZpV+rUPkFAACA0yD5BQAAgNOg7QEAAMCBcXtj61D5BQAAQJKYOnWq8uTJIw8PD/n7+2v79u3xzl2xYoVq1qwpHx8fpU+fXgEBAVq/fr3Vr2lz8nvq1Cl98MEHeuedd3Tp0iVJ0nfffadffvnF1l0CAADASSxZskS9evXS+++/r4MHD6py5cqqW7euzpw5E+f8H374QTVr1tS6deu0f/9+VatWTa+99poOHjxo1euaDMMwrA1227Ztqlu3roKCgvTDDz/o+PHjyps3r8aNG6effvpJy5cvt3aXuhkZbfU2cE78eQfWcHXheEHCZCrfLalDgAO5e3BKUodgtv/PG0kdgiTJP3d6q+ZXqFBBZcqU0bRp08xjfn5+atSokUJCQhK0j6JFi6pp06YaOnRogl/XpsrvwIEDNXr0aG3cuFFubm7m8WrVqmnXrl227BIAAAAOLDIyUjdu3LB4REZGxjn33r172r9/v2rVqmUxXqtWLe3cuTNBrxcdHa2bN28qU6ZMVsVpU/J79OhRvf7667HGfXx8FBERYcsuAQAAYANTMnmEhIQoQ4YMFo/4Krjh4eGKioqSr6+vxbivr68uXryYoPc9fvx43b59W2+99VaC5j9k02oPnp6eCgsLU548eSzGDx48qGzZstmySwAAADiwQYMGKTg42GLM3d39qduYTJataYZhxBqLy1dffaXhw4dr9erVypw5s1Vx2lT5bdasmQYMGKCLFy/KZDIpOjpaO3bsUN++fdWqVStbdgkAAAAH5u7urvTp01s84kt+vb295erqGqvKe+nSpVjV4CctWbJE7du319KlS/XKK69YHadNye+YMWOUM2dOZcuWTbdu3VKRIkX08ssvKzAwUB988IEtuwQAAIAtkrrf4eHDCm5ubvL399fGjRstxjdu3KjAwMB4t/vqq6/Upk0bLVq0SK+++qp1L/ovm9oeUqZMqYULF2rkyJE6ePCgoqOjVbp0aRUoUMCmIAAAAOBcgoOD1bJlS5UtW1YBAQGaMWOGzpw5o86dO0uKaaM4f/685s+fLykm8W3VqpUmTpyoihUrmqvGqVKlUoYMGRL8us91h7d8+fIpX758z7MLAAAAOKGmTZsqIiJCI0eOVFhYmIoVK6Z169YpV65ckqSwsDCLNX+nT5+uBw8eqGvXruratat5vHXr1po3b16CX9emdX4Nw9Dy5cu1ZcsWXbp0SdHRlmv0rlixwtpdss4vEox1fmEN1vlFQrHOL6yRnNb5PfjXzaQOQZJUOle6pA4hQWyq/Pbs2VMzZsxQtWrV5Ovrm6Cr8gAAAICkZlPyu2DBAq1YsUL16tWzdzwAAABAorFptYcMGTIob9689o7FoS1bvEgN6ryiwLIl1aJpYx3cv++p8/fv+0ktmjZWYNmSali3ppYvXRxrzqaNG9SkUX0F+JdQk0b1tWWT5RWRB/btVe9u76lOjZdVtoSftm7+PtY+pk+dosYN6qlS+TKqFlRBXTq21c9HDj/fm8VzMwxD06dOVu0alRVYrqQ6tWupUyd/f+Z2mzau15uNXlVF/+J6s9Gr2vzEMSFJSxcv0mt1aiigbAk1b/pGrGNx2AcD5V+isMWjdfOm8cbZ/b2O8i9RWFviOL6QNAzD0LTPJqtmtUqq4F9C7du01MkEHD/fb1yvNxrUU7nSxfRGg3ra/L3l8TN75nQ1a9pYgeVLq9rLAerVo4v+PP2HxZxNGzfovU7tVbVSBZUqVki//nrcru8NyUNQmXxaHvqu/tgwRncPTtFrVUskdUh4CpMpeTwchU3J7/DhwzVixAjdvXvX3vE4pA3frdP4cR+qXcd3tXDpCpUu468eXd7VxbALcc4/f+6cenbprNJl/LVw6Qq17dBJn3w4Vps2bjDPOXL4oAb3D1a9+g301fJVqle/gQb2C7ZIXO/evasChQqp/6D4l5fLlSu3+g/+QItXrNasLxYoS9Zs6tq5g65euWK/DwBW+2LuLC38cp4GDBqi+YuWycvbR13ebafbt2/Fu82Rwwc1yHxMrP73mOito48dEzHHYojadeysRUtXqnSZsurepZPCnjgWA4Mqa/3m7ebHpKnT43zNRQu+oK0pGZo3Z6YWzJ+rgYOHauHi5fL29tZ7Hds+9fg5fOigBvTtrVdfa6ilX6/Wq681VP++vSyOn/37flLTd5pr/qKl+nzGXEU9iNJ7ndrr7p075jl3795RqdKl1aNX30R9j0haaVK56+iJ8+r94dKkDgWwO5sueLtz547eeOMN7dixQ7lz51bKlCktnj9w4IDVgTjyBW+tmzVVYT8/DRoy3Dz2ZsNXVbV6DXXrGRxr/qRPP9EPW7do+eq15rGxo4br999+1dwFMRXgQf166/at25o0bYZ5TvfOHZUufXqNHTc+1j7LlvDTJ6GTVbX60xd7vnXrlqoGltPUGXNUvmKAtW81WXD0C94Mw1DtGi+rWYtWatOuo6SYe5zXrBakHr36qHGTt+PcbmC/3rp965YmT5tpHuvWuYPSp0+vseMmSJJaNXtLhf2KaPBjx2LjhvVUtXoNde/ZR1JM5ffmzZuaMPGzp8Z54rdf1atbZ83/aplqV6+sT0KnqNozjq/k6L92wZthGKpZrbKat2yltu07SYo5fqpXCVSv3n315ltxHz/9+/TS7du39Nnns8xjXd5tr/TpM+jDjyfEuc2VK1dU/eUAzZ63QP5ly1k8d/78Ob1au4YWL1+lwoX97PTukhYXvMXt7sEpeqv3DH279UhSh5KsJKcL3g6fSR4XvJXM6RgXvNlU+W3Tpo3279+vFi1aqHHjxmrYsKHFw5ncv39Pvx7/RRUDgyzGKwYE6cihg3Fuc/TwIVUMsJwfEBikY8d+0YP79yVJRw4fVoUnFnmuGBSkI4fj3mdCY125fKnSpkungoUK27wfPJ/z588pIvyyxTEQs9h3OR2O55iRpCOHD8U6zgKCKunw4UOSrDsW9+/7Sa9UCdTrr9XWqOFDdCUiwuL5u3fvavCAPuo/eIi8vX1seZtIJOfPnVN4+GUFBFYyj7m5uals2XI69Mzjp5LFWEBQ5acec7duxfxAtWb9TABI7my64G3t2rVav369KlWq9OzJ/3HXrl5TVFSUMnl5W4xn8vJSeHh4nNtERIQrk5fXE/O9FfXgga5duypvn8yKCA+XVybLfXpl8lZEPPt8mu3btmhw/77655+78vbx0WfTZ8szY0ar9wP7iAi/LEnyinUMeMVqT7DcLlyZMj2xTSYv8/6uXb2qqKioWPv18vKyOG6CKr2sV2rVUZYsWXXh/DlN+2ySOndoowVLvpabm5skacLHISpRsrSqVqth+xtFogj/9/sd1zkk7EL8x094eHicx8bD/T3JMAyNHxei0mX8lb9AweeMGgCSD5uS3xw5cih9+vQ2v2hkZKQiIyMtxu4pZbz3f3YET7ZFGobx1F7JJ58zd588Pv7kHD19n/EpW66CFi1boWtXr2rlimUa1Le35i1cEuuHJxLHurXfauzIYeavJ372ecw/Yh0Dz27piHXcxDEW+9iyfK1adR6t0pK/QEH5FS2m+rVr6Mcftqr6K7W0bctm7f1pjxYttX69btjf2jXfaPSIR8fP5H/7s+M6hzzr9BD3NnFvFDJmpE6cOKF58xfZEDWAF+q/1d2V6GxKfsePH6/+/fvr888/V+7cua3ePiQkRCNGjLAYG/j+UA0eMiyeLZIvz4yecnV1jVWRvXrlSqwqy0NeXrEruFevRMg1RQp5ZvCMmePtrYgIy4rMlSsRNiWsqVKnVo6cuZQjZy4VL1lKr9evrdUrv1bbDp2s3hesV6VqNRUv/uhK6Xv37kmKqeT6+GQ2j199xvc35piIfdw8/KuDZ8aMcnV1jfUXhytXIuI9FiXJxyezsmTNqjNn/pIk7f1pt86dPaOqQeUt5vUP7qHSZfw1Y86XT3u7sLOq1aqreImS5q+ffvx4x9r+IW9v7ziOjSvyimObD8eO0rYtmzXniwXyfeml530LAJCs2NTz26JFC23ZskX58uVTunTplClTJovHswwaNEjXr1+3ePTpP9CWUJJcypRuKuxXVHt27bQY37N7p0qUKh3nNsVLltKe3Zbzd+/coSJFiirFvxcPlihZMvY+d+5UiZJx79MahvHoBygSX5o0ac2/fOTImUt58+WXl7ePxff3/v172r9/r0rGc8xIUomSpWIdE7t37lDJkqUk2XYsStK1a1f198Uwc29vm/YdtXj5ai1autL8kKTgfgM1bGSIVe8dzy9NmrTKmTOX+ZEvX355e/to164d5jn379/Tvn17VeoZx8/ux7aRpN07f7Q45gzDUMiYkdr0/QbNmPOFsmXPYf83BABJzKbKb2ho6HO9qLu7e6wWB0de7aF5q9YaOnig/IoWU4mSpbRi+VJdDAtT4yYxa6dOmThBl/7+WyPHfiRJatzkbS39apEmfPyhXm/cREcOH9LqlSs05qNPzPt8u3krdWrbUvPmzFTVajW0dcsm7dmzS7PnLTDPuXPnts4+ds/r8+fP6bdfjytDhgx6KUtW3b1zR3NmTtfLVavJ28dH169d07IlX+nS3xf1Sq3aL+jTwZNMJpOatWilObOnK0eumIRmzqzp8vDwUJ169c3zhg4eIB/fzOZVGt5p3lId/z0mqlSroW3mY2KheZsWrdpoyOABKvLEsfjmvytI3LlzW9OnTlGNmrXk7e2jCxfO67NJn8rTM6Oq1YhZycHb2yfOi9xeypJV2bJnT8yPBglgMpnUvGUrzZ45Xbly5lbOXLk0a+Z0pfLwUN1XHx0/Hwzqr8yZfdWjd8zx06xFK7Vv00JzZ894dE7ZvUtzH2trGDt6hP63bo1CJ01VmjRpzP3AadOmk4eHhyTp+vVrCgsL0+VLlyRJf50+LSmmsszFkf8daVK5KV+OR9/P3Nm8VKJgNl29cUdnL15NwsgQF0dfBelFsyn5bd26tb3jcGi16tTT9WvXNGv6VIVfvqx8+Qto4mefK0vWbJKk8MuXdfFimHl+tuzZNXHq55ow7kMtW7xIPj6Z1XfgYNWoWcs8p2Sp0hrz0XhNmzJRn0+ZrOw5cihk3HgVe+zPn8d++UWd2z/6Xnz6cUxyXb9BIw0fHSIXV1f9+ecfWtNnla5dvaoMnp4qUrS4Zs5boHz5CyT2x4KnaN22gyL/+Ucfjhmpmzeuq1jxEvrs89lKkyatec7FixdkemyZrpKlymjsR+M1dcpETZsySdlz5NCH4yZY/Em8Vp16unbtmmZO/8x8LE76bLr5WHRxcdXJkye09tvVunnzprx9fFS2XHmFfPypxWsjeWvTrqP++SdSY0eP0I0b11W8RElNmzHH4nsYFhYmk8ujP+6VKl1GH348QZ9NDtVnkycpR44c+ujjTy2On2VLvpIkdWjb0uL1RowOUcNGb0iStm7ZrGEfDDI/N6Bfb0nSu+9103tdu9v/zSJJlCmSSxtm9TR/Pa5vY0nSl9/sVqdhC+LbDHAICV7n98aNG+aL3G7cuPHUubZcDOfIlV+8WPyGC2v819b5ReJhnV9YIzmt83v0XPw3uHmRimd3jCJKgiu/GTNmVFhYmDJnzixPT884rxB+eOVwVFSUXYMEAAAA7CHBye/mzZvNF7Nt2bIl0QICAAAAEkuCk98qVapIkh48eKCtW7eqXbt2ypGDK4EBAACSEs1d1rF6qbMUKVLok08+obUBAAAADsemdX5r1KihrVu32jkUAAAAIHHZtNRZ3bp1NWjQIP3888/y9/dXmjRpLJ5v0KCBXYIDAADAM9D3YJUEL3X2OBeX+AvGtq72wFJnSCiWOoM1WOoMCcVSZ7BGclrq7OfzyWOps2LZ/mNLnT0uOppEFQAAIDmgKGQdq3p+N2/erCJFisR5k4vr16+raNGi2r59u92CAwAAAOzJquQ3NDRUHTt2jPMObhkyZNC7776rCRMm2C04AAAAwJ6sSn4PHz6sOnXqxPt8rVq1tH///ucOCgAAAAljMiWPh6OwKvn9+++/lTJlynifT5EihS5fvvzcQQEAAACJwarkN1u2bDp69Gi8zx85ckRZsmR57qAAAACAxGBV8luvXj0NHTpU//zzT6zn7t69q2HDhql+/fp2Cw4AAABPZ0omD0dh1Tq/f//9t8qUKSNXV1d169ZNhQoVkslk0vHjx/XZZ58pKipKBw4ckK+vr9WBsM4vEoolXWAN1vlFQrHOL6yRnNb5PX7hdlKHIEnyy5rm2ZOSAavW+fX19dXOnTv13nvvadCgQXqYN5tMJtWuXVtTp061KfEFAAAAXgSrb3KRK1curVu3TlevXtXJkydlGIYKFCigjBkzJkZ8AAAAeBr+wGUVm+7wJkkZM2ZUuXLl7BkLAAAAkKhsTn4BAACQ9LgWxjpWrfYAAAAAODKSXwAAADgN2h4AAAAcmCPdWjg5oPILAAAAp0HyCwAAAKdB2wMAAIADo+vBOlR+AQAA4DSo/AIAADgySr9WofILAAAAp0HyCwAAAKdB2wMAAIAD4/bG1qHyCwAAAKdB8gsAAACnQdsDAACAA+P2xtah8gsAAACnQfILAAAAp0HbAwAAgAOj68E6VH4BAADgNKj8AgAAODJKv1ah8gsAAACnQfILAAAAp0HbAwAAgAPj9sbWofILAAAAp0HyCwAAAKdB2wMAAIAD4/bG1qHyCwAAAKdB5RcAAMCBUfi1DpVfAAAAOA2SXwAAADgN2h4AAAAcGX0PVqHyCwAAAKdB8gsAAACnQdsDAACAA+P2xtah8gsAAACnQeUXAADAgXGHN+tQ+QUAAIDTIPkFAACA06DtAQAAwIHR9WAdKr8AAABwGiS/AAAAcBq0PQAAADgwVnuwDpVfAAAAOA2SXwAAADgN2h4AAAAcGn0P1qDyCwAAAKdB5RcAAMCBccGbdUyGYRhJHQRii4yMVEhIiAYNGiR3d/ekDgfJHMcLEopjBdbgeHEM56/dS+oQJEnZPN2SOoQEIflNpm7cuKEMGTLo+vXrSp8+fVKHg2SO4wUJxbECa3C8OAaSX+vQ9gAAAODA6HqwDhe8AQAAwGmQ/AIAAMBp0PaQTLm7u2vYsGFcYIAE4XhBQnGswBocL46B1R6swwVvAAAADizsevK44C1LBi54AwAAQCIzccmbVej5BQAAgNMg+QUAAIDTIPlNhnLnzq3Q0NAEz//zzz9lMpl06NChRIsJkCSTyaRVq1YleP7w4cNVqlSpRIsHtrP2PAPnsHXrVplMJl27di2pQ4E1TMnk4SBIfpOhvXv3qlOnTnbd57x58+Tp6WnXfcL5hIWFqW7dukkdBmw0b948VaxYMdH2f/jwYb3zzjvKkSOHUqVKJT8/P02cODHRXg9Pd/HiRXXv3l158+aVu7u7cuTIoddee02bNm2Kd5vAwECFhYUpQ4YMLzBS4MXigrdkyMfHJ6lDgBO6f/++UqZM+dQ5L7300guKBonhm2++UcOGDRNt//v375ePj48WLFigHDlyaOfOnerUqZNcXV3VrVu3RHtdxPbnn38qKChInp6eGjdunEqUKKH79+9r/fr16tq1q3799ddY29y/f19ubm78d47/PCq/dvDtt9/K09NT0dHRkqRDhw7JZDKpX79+5jnvvvuu3nnnHUnSzp079fLLLytVqlTKkSOHevToodu3b5vnPvnnyF9//VWVKlWSh4eHihQpou+//z7OPz//8ccfqlatmlKnTq2SJUtq165dkmL+jNW2bVtdv35dJpNJJpNJw4cPT5wP4z8uOjpaH330kfLnzy93d3flzJlTY8aMkSQdPXpU1atXV6pUqeTl5aVOnTrp1q1b5m3btGmjRo0aaezYsfL19ZWnp6dGjBihBw8eqF+/fsqUKZOyZ8+uOXPmmLd52NKydOlSVa5cWalSpVK5cuV04sQJ7d27V2XLllXatGlVp04dXb582SLWuXPnys/PTx4eHipcuLCmTp0a536rVq0qDw8PLViwQJI0Z84cFS1aVO7u7sqSJYtF0vLkcTdgwAAVLFhQqVOnVt68eTVkyBDdv3/frp/5f13VqlXVrVs3devWTZ6envLy8tIHH3ygh6tQXr16Va1atVLGjBmVOnVq1a1bV7///rt5+4d/1Vm1apUKFiwoDw8P1axZU2fPnrV4nX/++UcbNmxQgwYN4oxj7ty5ypAhgzZu3GiOq3v37urVq5cyZswoX19fzZgxQ7dv31bbtm2VLl065cuXT//73//M+2jXrp0mTZqkKlWqKG/evGrRooXatm2rFStW2PtjwzN06dJFJpNJP/30k958800VLFhQRYsWVXBwsHbv3i0p5r/nzz//XA0bNlSaNGk0evToWG0PD4+vNWvWqFChQkqdOrXefPNN3b59W1988YVy586tjBkzqnv37oqKijK//r1799S/f39ly5ZNadKkUYUKFbR169Yk+CScQ1J3OzhY14Nk4Lldu3bNcHFxMfbt22cYhmGEhoYa3t7eRrly5cxzChYsaEybNs04cuSIkTZtWuPTTz81Tpw4YezYscMoXbq00aZNG/PcXLlyGZ9++qlhGIYRFRVlFCpUyKhZs6Zx6NAhY/v27Ub58uUNScbKlSsNwzCM06dPG5KMwoULG2vWrDF+++0348033zRy5cpl3L9/34iMjDRCQ0ON9OnTG2FhYUZYWJhx8+bNF/b5/Jf079/fyJgxozFv3jzj5MmTxvbt242ZM2cat2/fNrJmzWq88cYbxtGjR41NmzYZefLkMVq3bm3etnXr1ka6dOmMrl27Gr/++qsxe/ZsQ5JRu3ZtY8yYMcaJEyeMUaNGGSlTpjTOnDljGIbl9/a7774zjh07ZlSsWNEoU6aMUbVqVePHH380Dhw4YOTPn9/o3Lmz+bVmzJhhZMmSxfj666+NP/74w/j666+NTJkyGfPmzbPYb+7cuc1zzp8/b0ydOtXw8PAwQkNDjd9++8346aefzMeiYRgWx51hGMaoUaOMHTt2GKdPnza++eYbw9fX1/joo4/Mzw8bNswoWbJkonwv/iuqVKlipE2b1ujZs6fx66+/GgsWLDBSp05tzJgxwzAMw2jQoIHh5+dn/PDDD8ahQ4eM2rVrG/nz5zfu3btnGIZhzJ0710iZMqVRtmxZY+fOnca+ffuM8uXLG4GBgRavs2bNGiNfvnzmrx8/z3z88cdGpkyZjF27dlnElS5dOmPUqFHmY9PFxcWoW7euMWPGDOPEiRPGe++9Z3h5eRm3b9+O9/01b97caNy4sb0+LiRARESEYTKZjLFjxz51niQjc+bMxuzZs41Tp04Zf/75p7FlyxZDknH16lXDMB4dXzVr1jQOHDhgbNu2zfDy8jJq1aplvPXWW8Yvv/xifPvtt4abm5uxePFi876bNWtmBAYGGj/88INx8uRJ4+OPPzbc3d2NEydOJOZbd1oXr99LFg9HQfJrJ2XKlDE++eQTwzAMo1GjRsaYMWMMNzc348aNG0ZYWJghyTh+/LjRsmVLo1OnThbbbt++3XBxcTHu3r1rGIblD6X//e9/RooUKYywsDDz/I0bN8aZ/M6aNcs855dffjG/pmHEnMAyZMiQSO/eOdy4ccNwd3c3Zs6cGeu5GTNmGBkzZjRu3bplHlu7dq3h4uJiXLx40TCMmOQ3V65cRlRUlHlOoUKFjMqVK5u/fvDggZEmTRrjq6++Mgwj7u/tV199ZUgyNm3aZB4LCQkxChUqZP46R44cxqJFiyxiHDVqlBEQEGCx39DQUIs5WbNmNd5///14P4Mnk98njRs3zvD39zd/TfL7bFWqVDH8/PyM6Oho89iAAQMMPz8/48SJE4YkY8eOHebnwsPDjVSpUhlLly41DCPmv21Jxu7du81zjh8/bkgy9uzZYx7r2LGjERwcbP764Xlm4MCBRpYsWYwjR47EiqtSpUrmrx8emy1btjSPPTy3PZ40P27nzp1GypQpjQ0bNlj7seA57Nmzx5BkrFix4qnzJBm9evWyGIsr+ZVknDx50jzn3XffNVKnTm1RRKldu7bx7rvvGoZhGCdPnjRMJpNx/vx5i33XqFHDGDRo0PO8NcQjqZNeR0t+6fm1k6pVq2rr1q0KDg7W9u3bNXr0aH399df68ccfde3aNfn6+qpw4cLav3+/Tp48qYULF5q3NQxD0dHROn36tPz8/Cz2+9tvvylHjhwWPVjly5ePM4YSJUqY/50lSxZJ0qVLl1S4cGF7vlWndfz4cUVGRqpGjRpxPleyZEmlSZPGPBYUFKTo6Gj99ttv8vX1lSQVLVpULi6Puo18fX1VrFgx89eurq7y8vLSpUuXLPb/+Pf24b6KFy9uMfZwm8uXL+vs2bNq3769OnbsaJ7z4MGDWBexlC1b1vzvS5cu6cKFC3G+v/gsX75coaGhOnnypG7duqUHDx4offr0Cd4eMSpWrCjTY/cnDQgI0Pjx43Xs2DGlSJFCFSpUMD/n5eWlQoUK6fjx4+axFClSWHwvCxcuLE9PTx0/flzly5eXYRj69ttvtXjxYovXHT9+vG7fvq19+/Ypb968seJ6/Lh7eGw+edxJinW8StIvv/yihg0baujQoapZs6Y1Hweek/Fvy4wpAfe8ffy4iU/q1KmVL18+89e+vr7KnTu30qZNazH28Dg4cOCADMNQwYIFLfYTGRkpLy+vBL0HWIfbG1uH5NdOqlatqtmzZ+vw4cNycXFRkSJFVKVKFW3btk1Xr15VlSpVJMX0jL777rvq0aNHrH3kzJkz1phhGAk6gUmyuFjp4TYP+5Dx/FKlShXvc0/7Pj0+/uQFZSaTKc6xJ79vcX1vnxx7uM3D/585c6ZF0iTFJDCPezxZf9r7i8vu3bv19ttva8SIEapdu7YyZMigxYsXa/z48VbtB9aL63iL6/h7OPbTTz/p3r17qlSpksXzlStX1tq1a7V06VINHDgw1vbPOl7jO88cO3ZM1atXV8eOHfXBBx9Y8c5gDwUKFJDJZNLx48fVqFGjp859/BwQH2vPW9HR0XJ1ddX+/ftjnXMeT5iBpMIFb3by8ssv6+bNmwoNDVWVKlVkMplUpUoVbd26VVu3bjUnv2XKlNEvv/yi/Pnzx3q4ucW+J3bhwoV15swZ/f333+axvXv3Wh2fm5ubxcUIsF6BAgWUKlWqOJcJKlKkiA4dOmRx4eKOHTvk4uISq/qR2Hx9fZUtWzb98ccfsY6xPHnyxLtdunTplDt37qcug/S4HTt2KFeuXHr//fdVtmxZFShQQH/99Ze93oZTeXgB0uNfFyhQQEWKFNGDBw+0Z88e83MRERE6ceKExV+JHjx4oH379pm//u2333Tt2jXzX31Wr16tV199NVYiUr58eX333XcaO3asPv74Y7u8l19++UXVqlVT69atzReD4sXKlCmTateurc8++8zinPRQYq/hW7p0aUVFRenSpUuxzkGsJJE4TMnkf46C5NdOMmTIoFKlSmnBggWqWrWqpJiE+MCBAzpx4oR5bMCAAdq1a5e6du2qQ4cO6ffff9c333yj7t27x7nfmjVrKl++fGrdurWOHDmiHTt26P3335eUsD9pPZQ7d27dunVLmzZtUnh4uO7cufNc79cZeXh4aMCAAerfv7/mz5+vU6dOaffu3Zo9e7aaN28uDw8PtW7dWj///LO2bNmi7t27q2XLluY/Db9Iw4cPV0hIiCZOnKgTJ07o6NGjmjt3riZMmPDM7caPH69Jkybp999/14EDBzR58uQ45+bPn19nzpzR4sWLderUKU2aNEkrV65MjLfzn3f27FkFBwfrt99+01dffaXJkyerZ8+eKlCggBo2bKiOHTvqxx9/1OHDh9WiRQtly5bNYsmylClTqnv37tqzZ48OHDigtm3bqmLFiuYWqactcRYQEKD//e9/GjlypD799NPneh8PE9+aNWsqODhYFy9e1MWLF2OtRILEN3XqVEVFRal8+fL6+uuv9fvvv+v48eOaNGmSAgICEvW1CxYsqObNm6tVq1ZasWKFTp8+rb179+qjjz7SunXrEvW1gYQg+bWjatWqKSoqypzoZsyYUUWKFJGPj4+5SlOiRAlt27ZNv//+uypXrqzSpUtryJAh5h7dJ7m6umrVqlW6deuWypUrpw4dOpj/jOjh4ZHg2AIDA9W5c2c1bdpUPj4+Gjdu3PO9WSc1ZMgQ9enTR0OHDpWfn5+aNm2qS5cuKXXq1Fq/fr2uXLmicuXK6c0331SNGjU0ZcqUJImzQ4cOmjVrlubNm6fixYurSpUqmjdv3lMrv5LUunVrhYaGaurUqSpatKjq169vsazW4xo2bKjevXurW7duKlWqlHbu3KkhQ4Ykxtv5z2vVqpXu3r2r8uXLq2vXrurevbv5Rjdz586Vv7+/6tevr4CAABmGoXXr1ln82Tl16tQaMGCAmjVrpoCAAKVKlcrc33vq1CmdPHlStWvXjvf1g4KCtHbtWg0ZMkSTJk2y+X0sW7ZMly9f1sKFC5UlSxbzo1y5cjbvE7bJkyePDhw4oGrVqqlPnz4qVqyYatasqU2bNmnatGmJ/vpz585Vq1at1KdPHxUqVEgNGjTQnj17lCNHjkR/beBZTMbDzng4jB07dqhSpUo6efKkxUUIABxP1apVVapUKZtvNTxv3jz16tUr3j9lT5gwQd9//z0VN+A/7PKtB0kdgiTJJ61jXErmGFE6uZUrVypt2rQqUKCATp48qZ49eyooKIjEF8AzZc+eXYMGDUrqMAAg2SD5dQA3b95U//79dfbsWXl7e+uVV17hinoACfLWW28ldQgAkKzQ9gAAAODAwpNJ24O3g7Q9cMEbAAAAnIZjpOgAAACIE3d4sw6VXwAAADgNkl8AAAA4DdoeAAAAHJgj3Vo4OaDyCwAAAKdB8gsAAACnQdsDAACAA2O1B+tQ+QUAAIDTIPkFAACA0yD5BQAAgNMg+QUAAIDT4II3AAAAB8YFb9ah8gsAAACnQfILAAAAp0HbAwAAgAPj9sbWofILAAAAp0HyCwAAAKdB2wMAAIADY7UH61D5BQAAgNOg8gsAAODAKPxah8ovAAAAksTUqVOVJ08eeXh4yN/fX9u3b3/q/G3btsnf318eHh7KmzevPv/8c6tfk+QXAAAAL9ySJUvUq1cvvf/++zp48KAqV66sunXr6syZM3HOP336tOrVq6fKlSvr4MGDGjx4sHr06KGvv/7aqtc1GYZh2OMNAAAA4MW7GRmd1CFIktK5W1dTrVChgsqUKaNp06aZx/z8/NSoUSOFhITEmj9gwAB98803On78uHmsc+fOOnz4sHbt2pXg16XyCwAAgOcWGRmpGzduWDwiIyPjnHvv3j3t379ftWrVshivVauWdu7cGec2u3btijW/du3a2rdvn+7fv5/gOEl+AQAA8NxCQkKUIUMGi0dcFVxJCg8PV1RUlHx9fS3GfX19dfHixTi3uXjxYpzzHzx4oPDw8ATHyWoPAAAADiy53N540KBBCg4Othhzd3d/6jamJxYpNgwj1tiz5sc1/jQkvwAAAHhu7u7uz0x2H/L29parq2usKu+lS5diVXcfeumll+KcnyJFCnl5eSU4TtoeAAAAHJjJlDwe1nBzc5O/v782btxoMb5x40YFBgbGuU1AQECs+Rs2bFDZsmWVMmXKBL82yS8AAABeuODgYM2aNUtz5szR8ePH1bt3b505c0adO3eWFNNG0apVK/P8zp0766+//lJwcLCOHz+uOXPmaPbs2erbt69Vr0vbAwAAAF64pk2bKiIiQiNHjlRYWJiKFSumdevWKVeuXJKksLAwizV/8+TJo3Xr1ql379767LPPlDVrVk2aNEmNGze26nVZ5xcAAMCB3bmXPFK51G7J48K7Z6HtAQAAAE6D5BcAAABOg55fAAAAR+YY3QbJBpVfAAAAOA2SXwAAADgN2h4AAAAcWHK5vbGjoPILAAAAp0HlFwAAwIFZe2thZ0flFwAAAE6D5BcAAABOg9sbAwAAwGlQ+QUAAIDTIPkFAACA0yD5BQAAgNMg+QUAAIDTIPkFAACA0yD5BQAAgNMg+QUAAIDTIPkFAACA0yD5BQAAgNP4P69pHxWlOYB6AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 900x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Visualization\n",
    "import seaborn as sns\n",
    "\n",
    "a = points.corr()\n",
    "\n",
    "plt.subplots(figsize=(9,9))\n",
    "sns.heatmap(a, annot = True, vmax = 1, square = True, cmap = \"Blues\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4155949-d4c7-4eeb-803f-1d00fd03d0b4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}