{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<h1 style=\"text-align: center\">Extracting Features From the History of Crimes Recently Committed and Within the Area of Influence: A Feature Engineering Strategy</h1>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This document presents a feature engineering strategy that supports the predictive power of an ensemble-based approach to crime classification. First, to each crime in the dataset, several features are extracted from the history of recent incidents that are within its area of influence. Next, a methodology to develop a stacked generalization ensemble is proposed. In particular, several base models are trained on stratified subsets of the whole dataset that are different from each other. Then, the first-level predictions, i.e., the outputs from the base models, are combined according to the stacking technique the meta-model implements to make the second-level predictions. Experimental results show that training a classifier to combine the predictions from the base models outperforms the straightforward stacking technique of soft voting. Nevertheless, the lowest multi-class logarithmic loss obtained, i.e., 2.56044, is the result of combining the second-level predictions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<span style=\"font-size: 128.6%;\"><b>Contents</b></span>\n",
"\n",
"<ul style=\"list-style-type: none; padding-left: 0 !important;\">\n",
" <li><a href=\"#Project-Description\">Project Description</a></li>\n",
" <li><a href=\"#0.-Requirements\">0. Requirements</a></li>\n",
" <li><a href=\"#1.-Data-Preprocessing\">1. Data Preprocessing</a></li>\n",
" <li><a href=\"#2.-Feature-Engineering\">2. Feature Engineering</a></li>\n",
" <li>\n",
" <a href=\"#3.-Crime-Classification\">3. Crime Classification</a>\n",
" <ul style=\"list-style-type: none; padding-left: 1em !important;\">\n",
" <li><a href=\"#3.1-Train-Test-Split\">3.1 Train-Test Split</a></li>\n",
" <li><a href=\"#3.2-Baseline-Model\">3.2 Baseline Model</a></li>\n",
" <li><a href=\"#3.3-Discriminative-Features\">3.3 Discriminative Features</a></li>\n",
" <li>\n",
" <a href=\"#3.4-Stacking-Ensemble\">3.4 Stacking Ensemble</a>\n",
" <ul style=\"list-style-type: none; padding-left: 1em !important;\">\n",
" <li><a href=\"#3.4.1-Results\">3.4.1 Results</a></li>\n",
" <li><a href=\"#3.4.2-Late-Submission\">3.4.2 Late Submission</a></li>\n",
" </ul>\n",
" </li>\n",
" </ul>\n",
" </li>\n",
" <li><a href=\"#4.-Conclusion\">4. Conclusion</a></li>\n",
"</ul>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Project Description"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[San Francisco Crime Classification](https://www.kaggle.com/c/sf-crime/) is a Kaggle challenge related to multi-class classification. In particular, each record in the dataset must be classified into one out of 39 categories of crime. There are 12 years of criminal records on which predictive models must be trained and evaluated. To get a sense of the dataset, an overview of its attributes is presented below:\n",
"\n",
"* **Dates**: the DateTime the crime was committed at\n",
"* **Category**: the category of the committed crime. (This is the target variable.)\n",
"* **Descript**: a further description of the crime\n",
"* **DayOfWeek**: day of the week on which the crime was committed\n",
"* **PdDistrict**: the Police Department District that attended the crime incident\n",
"* **Resolution**: an (optional) explanation on how the crime was resolved\n",
"* **Address**: \"the approximate street address of the crime incident\"\n",
"* **X**: the longitude of the location the crime was committed\n",
"* **Y**: the latitude of the location the crime was committed\n",
"\n",
"In total, the dataset comprises 1,762,311 crime records, 878,049 of which have a category (i.e., the training set). These data range from January 2003 to May 2015.\n",
"\n",
"Lastly, the performance of predictive models is evaluated on the test set using the multi-class logarithmic loss metric."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 0. Requirements"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In order to encourage reproducibility, the following is a list of technologies used, as well as their respective version:\n",
"\n",
"1. Python **3.7.2**.\n",
"2. NumPy **1.17.2**.\n",
"3. SciPy **1.3.1**.\n",
"4. Scikit-learn **0.21.3**.\n",
"5. pandas **0.25.1**.\n",
"6. Matplotlib **3.1.1**.\n",
"7. Numba **0.48.0**.\n",
"8. tabulate **0.8.6**."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import copy\n",
"import os\n",
"import re\n",
"\n",
"import joblib\n",
"import matplotlib.pyplot as plt\n",
"import numba\n",
"import numpy as np\n",
"import pandas as pd\n",
"import tabulate\n",
"from IPython.display import display\n",
"from IPython.display import HTML\n",
"from sklearn.ensemble import GradientBoostingClassifier\n",
"from sklearn.ensemble import RandomForestClassifier\n",
"from sklearn.linear_model import LogisticRegression\n",
"from sklearn.metrics import log_loss\n",
"from sklearn.model_selection import cross_val_predict\n",
"from sklearn.model_selection import GridSearchCV\n",
"from sklearn.model_selection import ParameterGrid\n",
"from sklearn.model_selection import StratifiedKFold\n",
"from sklearn.preprocessing import StandardScaler"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"CURRENT_PATH = os.path.abspath(os.getcwd())\n",
"DATA_PATH = os.path.join(CURRENT_PATH, 'data')\n",
"\n",
"DATE_FORMAT = '%Y-%m-%d %H:%M:%S'\n",
"RANDOM_STATE = 91\n",
"\n",
"ALGORITHMS = {\n",
" 'logit': {\n",
" 'estimator': LogisticRegression(\n",
" penalty='l2',\n",
" solver='liblinear',\n",
" multi_class='ovr'\n",
" ),\n",
" 'predict_proba': True,\n",
" 'param_grid': None\n",
" },\n",
" 'gb': {\n",
" 'estimator': GradientBoostingClassifier(),\n",
" 'predict_proba': True,\n",
" 'param_grid': None\n",
" },\n",
" 'rf': {\n",
" 'estimator': RandomForestClassifier(),\n",
" 'predict_proba': True,\n",
" 'param_grid': None\n",
" }\n",
" }"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"FNAMES = {\n",
" 'original-train': 'train.csv',\n",
" 'train': 'projected-train.csv',\n",
" 'original-test': 'test.csv',\n",
" 'test': 'projected-test.csv',\n",
" 'sample-submission': 'sampleSubmission.csv',\n",
" 'crime-dataset': 'crime-dataset.csv',\n",
" 'baseline': 'baseline-models.csv'\n",
" }\n",
"\n",
"FNAMES = {key: os.path.join(DATA_PATH, fname) for key, fname in FNAMES.items()}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Data Preprocessing"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's map each crime category to its numerical representation."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"with open(FNAMES['sample-submission']) as f:\n",
" for i, row in enumerate(f):\n",
" row = row.rstrip('\\n')\n",
" CRIME_CATEGORY = {category: j-1 for j, category in enumerate(row.split(',')) if j > 0}\n",
" \n",
" break"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"PROJECTION = ['Id', 'Dates', 'Category', 'DayOfWeek', 'PdDistrict', 'X', 'Y']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As specified by the variable `PROJECTION`, let's project (or filter) the datasets."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def dataset_projection(in_fname, out_fname, projection=PROJECTION):\n",
" valid_columns = []\n",
" \n",
" insert_id = False\n",
" header, data = [], []\n",
" with open(in_fname) as f:\n",
" for i, row in enumerate(f):\n",
" row = row.rstrip('\\n')\n",
" \n",
" if i == 0:\n",
" for j, col in enumerate(row.split(',')):\n",
" if col in projection:\n",
" header.append(col)\n",
" valid_columns.append(j)\n",
" \n",
" insert_id = True if 'Id' not in header else False\n",
"\n",
" continue\n",
" \n",
" for old in re.findall(r'\"[^\"]+\"', row):\n",
" new = re.sub(r',', '|', old) \n",
" row = row.replace(old, new, 1)\n",
" \n",
" record = [\n",
" re.sub(r'\\|', ',', col).strip()\n",
" for j, col in enumerate(row.split(',')) if j in valid_columns\n",
" ]\n",
" \n",
" if len(record) != len(valid_columns):\n",
" print('({}), Malformed columns at line {}'.format(in_fname, i+1))\n",
" continue\n",
" elif insert_id:\n",
" record.insert(0, i-1)\n",
" \n",
" data.append(record)\n",
" \n",
" if insert_id:\n",
" header.insert(0, 'Id')\n",
" \n",
" return header, data"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"datasets = [\n",
" [FNAMES['original-train'], FNAMES['train']],\n",
" [FNAMES['original-test'], FNAMES['test']]\n",
" ]\n",
"for in_fname, out_fname in datasets:\n",
" if os.path.isfile(out_fname):\n",
" continue\n",
" \n",
" header, data = dataset_projection(in_fname, out_fname)\n",
" df = pd.DataFrame(data, columns=header)\n",
" \n",
" df['Dates'] = pd.to_datetime(df['Dates'], format=DATE_FORMAT)\n",
" df = df.sort_values(by=['Dates'])\n",
" df['Dates'] = df['Dates'].dt.strftime(DATE_FORMAT)\n",
" \n",
" df.to_csv(out_fname, index=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, let's append the test set to the training one."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"if not os.path.isfile(FNAMES['crime-dataset']):\n",
" columns = copy.deepcopy(PROJECTION)\n",
" columns.remove('Category')\n",
" columns.insert(0, 'Dataset')\n",
" \n",
" df = None\n",
" for in_fname, dataset_type in [[FNAMES['train'], 'train'], [FNAMES['test'], 'test']]:\n",
" dataset = pd.read_csv(in_fname)\n",
" dataset['Dataset'] = dataset_type\n",
" dataset = dataset[columns]\n",
" \n",
" df = dataset.copy(deep=True) if df is None else df.append(dataset)\n",
" \n",
" df['Dates'] = pd.to_datetime(df['Dates'], format=DATE_FORMAT)\n",
" df = df.sort_values(by=['Dates', 'Dataset', 'Id'])\n",
" df['Dates'] = df['Dates'].dt.strftime(DATE_FORMAT)\n",
" \n",
" df = df[columns]\n",
" \n",
" df.to_csv(FNAMES['crime-dataset'], index=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. Feature Engineering"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This process explores the ability to extract valuable information from the history of crimes recently committed and relatively close to each other. In other words, to each crime $c$, whose position in the map is depicted by the white filled circle and area of influence corresponds to the red filled circle, the history of crimes comprises all those from the previous hours before $c$ and within the area of influence."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![San Francisco Map](images/map.png \"San Francisco Map\")\n",
"<div align=\"center\" style=\"margin-top: 10px;\"><b>Figure 1</b>: A crime instance in San Francisco, which is is depicted by the white filled circle and whose area of influence corresponds to the red filled circle. Credits: <a href=\"https://www.google.com/maps\">Google Maps</a> and <a href=\"https://www.mapdevelopers.com/draw-circle-tool.php\">Map Developers</a></div>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Formally, let $x_{i}$ be each crime in the dataset and $S_{i}$ be the set comprising all crimes from the previous $t_{w}$ hours before $x_{i}$, then:\n",
"\n",
"\\begin{equation*}\n",
"S_{i} = \\{ x_{j} | (x_{j}^{Dates} < x_{i}^{Dates}) \\wedge (hours(x_{j}^{Dates}, x_{i}^{Dates}) < t_{w}) \\}_{_{j \\neq i}^{j = 1}}^{N}\n",
"\\end{equation*}\n",
"\n",
"\\begin{equation*}\n",
"S_{i}^{PdDistrict} = \\{ x_{j} | (x_{j}^{Dates} < x_{i}^{Dates}) \\wedge (hours(x_{j}^{Dates}, x_{i}^{Dates}) < t_{w}) \\wedge (x_{j}^{PdDistrict} = x_{i}^{PdDistrict}) \\}_{_{j \\neq i}^{j = 1}}^{N}\n",
"\\end{equation*}\n",
"\n",
"\\begin{equation*}\n",
"S_{i}^{AreaOfInfluence} = \\{ x_{j} | (x_{j}^{Dates} < x_{i}^{Dates}) \\wedge (hours(x_{j}^{Dates}, x_{i}^{Dates}) < t_{w}) \\wedge (distance(x_{i}^{Y}, x_{i}^{X}, x_{j}^{Y}, x_{j}^{X}) \\leq r) \\}_{_{j \\neq i}^{j = 1}}^{N}\n",
"\\end{equation*}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Afterward, several features are created by applying the aggregation function of count to each set of crimes as follows:\n",
"\n",
"\\begin{equation*}\n",
"x_{i}^{Agg_{1}} = |S_{i}|\n",
"\\end{equation*}\n",
"\n",
"\\begin{equation*}\n",
"x_{i}^{Agg_{2}} = |S_{i}^{PdDistrict}|\n",
"\\end{equation*}\n",
"\n",
"\\begin{equation*}\n",
"x_{i}^{Agg_{3}} = |S_{i}^{AreaOfInfluence}|\n",
"\\end{equation*}\n",
"\n",
"where,\n",
"\n",
"* $N$ is the number of total crimes in the dataset,\n",
"* $| \\cdot |$ is the cardinality of a set,\n",
"* $x_{i}^{Dates}$ is the DateTime when crime $x_{i}$ was committed,\n",
"* $hours(x_{j}^{Dates}, x_{i}^{Dates})$ is a function calculating the number of hours between the DateTimes $x_{j}^{Dates}$ and $x_{i}^{Dates}$,\n",
"* $t_{w}$ is the maximum number of hours between previous crime $x_{j}$ and crime $x_{i}$ in order to include the former in the set $S_{i}$,\n",
"* $x_{i}^{PdDistrict}$ is the Police Department District that attended the crime incident,\n",
"* $S_{i}^{PdDistrict}$ is the set of crimes attended by the same Police Department District,\n",
"* $x_{i}^{Y}$ and $x_{i}^{X}$ correspond to latitude and longitude, respectively,\n",
"* $distance(x_{i}^{Y}, x_{i}^{X}, x_{j}^{Y}, x_{j}^{X})$ is a function calculating the distance, in kilometers, between previous crime $x_{j}$ and crime $x_{i}$,\n",
"* $r$ is the maximum distance in kilometers between previous crime $x_{j}$ and crime $x_{i}$ in order to include the former in the set $S_{i}$,\n",
"* Note that $r$ can also be seen as the radius of the circle representing the area of influence of $x_{i}$,\n",
"* $S_{i}^{AreaOfInfluence}$ is the set of crimes within the area of influence of $x_{i}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Before continuing, some clarifications must be provided, namely:\n",
"\n",
"* $t_{w}$ takes a value from the set $\\{ 12, 24, 72, 168, 336 \\}$.\n",
"* $r$ takes a value from the set $\\{ 1, 2, 4, 8, 16 \\}$.\n",
"* The three aggregated features described above result from assigning one value to $t_{w}$ and one to $r$ from their corresponding sets."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On the other hand, the following list comprises all attributes derived from the DateTime $x_{i}^{Dates}$:\n",
"\n",
"1. Year,\n",
"2. Month,\n",
"3. Quarter,\n",
"4. Triannual,\n",
"5. Semester,\n",
"6. Day,\n",
"7. Day of week,\n",
"8. Whether or not the day of the week is a working day,\n",
"9. Fortnight,\n",
"10. Hour,\n",
"11. Four-hour period the hour of the crime belongs to,\n",
"12. Six-hour period the hour of the crime belongs to,\n",
"13. Twelve-hour period the hour of the crime belongs to."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In summary, the feature engineering process creates 13 date-derived features plus three aggregated features per each combination of values $t_{w}$ and $r$ take. To create these features, four raw features are used, namely: $Dates$, $PdDistrict$, $X$, and $Y$. Lastly, bear in mind that terms *feature* and *attribute* are used interchangeably throughout this document."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"@numba.jit(nopython=True)\n",
"def compute_distance(\n",
" lat_1, lon_1,\n",
" lat_2, lon_2):\n",
" \"\"\"Compute distance between two locations.\n",
" \n",
" Returns\n",
" -------\n",
" float\n",
" Distance in KM.\n",
" \n",
" Source: <https://stackoverflow.com/questions/19412462/>\n",
" \"\"\"\n",
" # Approximate radius of earth in KM\n",
" earth_radius = 6373.0\n",
" \n",
" lat_1 = np.radians(lat_1)\n",
" lon_1 = np.radians(lon_1)\n",
" \n",
" lat_2 = np.radians(lat_2)\n",
" lon_2 = np.radians(lon_2)\n",
" \n",
" lon_dist = lon_2 - lon_1\n",
" lat_dist = lat_2 - lat_1\n",
" \n",
" a = (np.square(np.sin(lat_dist/2))\n",
" + np.cos(lat_1)\n",
" * np.cos(lat_2)\n",
" * np.square(np.sin(lon_dist/2)))\n",
" \n",
" c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a))\n",
" \n",
" return earth_radius * c"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"@numba.jit(nopython=True)\n",
"def compute_aggregated_features(data, time_window, crime_radius):\n",
" \"\"\"Compute aggregated features.\n",
" \n",
" Parameters\n",
" ----------\n",
" data : np.ndarray, dtype('int64')\n",
" A Numpy-like array of shape \"(n, m)\", where \"n\" is the number\n",
" of records and \"m\" is the number of columns (or attributes).\n",
" The strict order of the columns is presented below:\n",
" Dataset,\n",
" Id,\n",
" Dates,\n",
" PdDistrict,\n",
" X - Longitude,\n",
" Y - Latitude\n",
" time_window : int\n",
" Time window (in hours).\n",
" crime_radius : list\n",
" List of integers, each of which representing a radius in kilometers.\n",
" \"\"\"\n",
" n = len(data)\n",
" \n",
" # Let's transform the time window into seconds\n",
" time_window = time_window * 60 * 60\n",
" \n",
" aggregated_features = []\n",
" for i in range(n):\n",
" ts = data[i,2] \n",
" \n",
" lower_ts = ts - time_window\n",
" \n",
" mask = ((lower_ts < data[:,2])\n",
" & (data[:,2] < ts))\n",
" \n",
" historical_data = data[mask]\n",
" m = len(historical_data)\n",
" \n",
" police_district = data[i,3]\n",
" \n",
" feature_vector = [\n",
" int(data[i,0]),\n",
" int(data[i,1]),\n",
" m, # number of crimes within the time window\n",
" 0 # number of crimes attended by the same police department district\n",
" ]\n",
" feature_vector = feature_vector + [0 for j in crime_radius]\n",
" \n",
" lat_1 = data[i,5]\n",
" lon_1 = data[i,4]\n",
" \n",
" for j in range(m):\n",
" feature_vector[3] += 1 if police_district == historical_data[j,3] else 0\n",
" \n",
" lat_2 = historical_data[j,5]\n",
" lon_2 = historical_data[j,4]\n",
" \n",
" # Let's compute the number of crimes within each given radius\n",
" distance = compute_distance(lat_1, lon_1, lat_2, lon_2)\n",
" \n",
" for k, rad in enumerate(crime_radius):\n",
" feature_vector[4+k] += 1 if distance <= rad else 0\n",
" \n",
" aggregated_features.append(feature_vector)\n",
" \n",
" return aggregated_features"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"def derive_date_attributes(df, column, drop_column=True):\n",
" df['Year'] = df[column].dt.year\n",
" \n",
" df['Month'] = df[column].dt.month\n",
" df['Quarter'] = df[column].dt.quarter\n",
" \n",
" df['Triannual'] = 0\n",
" for i, (min_m, max_m) in enumerate([[1, 4], [5, 8], [9, 12]]):\n",
" df.loc[((min_m <= df['Month']) & (df['Month'] <= max_m)), 'Triannual'] = i + 1\n",
" \n",
" df['Semester'] = 1\n",
" df.loc[(df['Quarter'] > 2), 'Semester'] = 2\n",
" \n",
" df['Day'] = df[column].dt.day\n",
" df['DayOfWeek'] = df[column].dt.dayofweek\n",
" \n",
" df['WorkingDay'] = 1\n",
" df.loc[(df['DayOfWeek'] > 4), 'WorkingDay'] = 0\n",
" \n",
" df['Fortnight'] = 1\n",
" df.loc[(df['Day'] > 15), 'Fortnight'] = 2\n",
" \n",
" df['Hour'] = df[column].dt.hour\n",
" \n",
" hourly_periods = {\n",
" 'four': [[i, i+3] for i in range(0, 24, 4)],\n",
" 'six': [[i, i+5] for i in range(0, 24, 6)],\n",
" 'twelve': [[i, i+11] for i in range(0, 24, 12)],\n",
" }\n",
" \n",
" for str_period, period in hourly_periods.items():\n",
" period_column = '{}HourPeriod'.format(str_period.title())\n",
" df[period_column] = 0\n",
" \n",
" for i, (min_hr, max_hr) in enumerate(period):\n",
" df.loc[((min_hr <= df['Hour']) & (df['Hour'] <= max_hr)), period_column] = i + 1\n",
" \n",
" if drop_column:\n",
" df = df.drop(columns=[column])\n",
" \n",
" return df"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def feature_engineering(\n",
" df, data, time_windows,\n",
" idx_to_dataset, idx_to_district,\n",
" crime_radius=[1, 2, 4, 8, 16]):\n",
" \"\"\"Compute the process of feature engineering.\"\"\"\n",
" df = df[['Dataset', 'Id', 'PdDistrict', 'Dates']]\n",
" df = df.replace({'Dataset': idx_to_dataset, 'PdDistrict': idx_to_district})\n",
" \n",
" df = derive_date_attributes(df, 'Dates')\n",
" \n",
" crime_radius = np.array(crime_radius, dtype=int)\n",
" \n",
" agg_ds_fname = os.path.join(DATA_PATH, 'agg-dataset-{}H.csv')\n",
" cum_ds_fname = os.path.join(DATA_PATH, 'agg-dataset-{}H-cumulative.csv') \n",
" \n",
" for i, time_window in enumerate(time_windows):\n",
" agg_ds_fname_1 = agg_ds_fname.format(time_window)\n",
" cum_ds_fname_1 = cum_ds_fname.format(time_window)\n",
" \n",
" if (os.path.isfile(cum_ds_fname_1)\n",
" or (i == 0 and os.path.isfile(agg_ds_fname_1))):\n",
" continue\n",
" elif os.path.isfile(agg_ds_fname_1):\n",
" agg_ds = pd.read_csv(agg_ds_fname_1)\n",
" \n",
" for col in agg_ds.columns:\n",
" if col in ['Dataset']:\n",
" continue\n",
" \n",
" agg_ds[col] = pd.to_numeric(agg_ds[col])\n",
" else:\n",
" agg_ds = compute_aggregated_features(data, time_window, crime_radius)\n",
" \n",
" prefix = '{}H_'.format(time_window)\n",
" agg_ds_columns = (['Dataset', 'Id']\n",
" + [(prefix + col) for col in ['Crimes', 'CrimesAttendedByPdDistrict']]\n",
" + [(prefix + 'CrimesWithin{}KMRad'.format(rad)) for rad in crime_radius])\n",
" \n",
" agg_ds = pd.DataFrame(agg_ds, columns=agg_ds_columns)\n",
" agg_ds = agg_ds.astype({col: 'int32' for col in agg_ds_columns})\n",
" \n",
" agg_ds['Dataset'] = agg_ds['Dataset'].map(idx_to_dataset)\n",
" \n",
" agg_ds.to_csv(agg_ds_fname_1, index=False)\n",
" \n",
" if i == 0:\n",
" continue\n",
" \n",
" cum_ds = agg_ds.copy(deep=True)\n",
" agg_ds = None\n",
" \n",
" for j in range(i):\n",
" agg_ds_fname_2 = agg_ds_fname.format(time_windows[j])\n",
" \n",
" agg_ds = pd.read_csv(agg_ds_fname_2)\n",
" \n",
" for col in agg_ds.columns:\n",
" if col in ['Dataset']:\n",
" continue\n",
" \n",
" agg_ds[col] = pd.to_numeric(agg_ds[col])\n",
" \n",
" cum_ds = pd.merge(agg_ds, cum_ds, on=['Dataset', 'Id'], how='inner')\n",
" \n",
" cum_ds = pd.merge(df, cum_ds, on=['Dataset', 'Id'], how='inner')\n",
" \n",
" cum_ds.to_csv(cum_ds_fname_1, index=False)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"df = pd.read_csv(FNAMES['crime-dataset'])\n",
"\n",
"df['Dates'] = pd.to_datetime(df['Dates'], format=DATE_FORMAT)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"DATASET_TO_IDX = {dataset: i for i, dataset in enumerate(df['Dataset'].unique())}\n",
"IDX_TO_DATASET = {i: dataset for dataset, i in DATASET_TO_IDX.items()}\n",
"\n",
"df['Dataset'] = df['Dataset'].map(DATASET_TO_IDX)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"DISTRICT_TO_IDX = {district: i for i, district in enumerate(df['PdDistrict'].unique())}\n",
"IDX_TO_DISTRICT = {i: district for district, i in DISTRICT_TO_IDX.items()}\n",
"\n",
"df['PdDistrict'] = df['PdDistrict'].map(DISTRICT_TO_IDX)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"PROJECTION = (['Dataset']\n",
" + [col for col in PROJECTION if col not in ['Category', 'DayOfWeek']])\n",
"\n",
"df = df[PROJECTION]\n",
"df = df.sort_values(by=['Dates', 'Dataset', 'Id'])\n",
"\n",
"crimes = df.copy(deep=True)\n",
"crimes['ts'] = crimes['Dates'].values.astype(np.int64) // 10 ** 9\n",
"crimes = crimes[[('ts' if col == 'Dates' else col) for col in PROJECTION]].to_numpy().astype(float)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Bear in mind that running the feature engineering process takes a long time; approximately 7.3 hours."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"TIME_WINDOWS = [12, 24, 72, 168, 336]\n",
"\n",
"feature_engineering(df, crimes, TIME_WINDOWS, IDX_TO_DATASET, IDX_TO_DISTRICT)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Crime Classification"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first subsection presents how the whole dataset is split into training, validation, and test sets. Then, the baseline model and the discriminative power of each set of features are discussed in subsections 3.2 and 3.3, respectively. Finally, the methodology to build stacked generalization ensembles is developed in subsection 3.4."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3.1 Train-Test Split"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's split the whole dataset of crimes into training, validation and test datasets. Please recall that the original training and test datasets, as downloaded from Kaggle, were merged to run the feature engineering process. Lastly, the training dataset will be split into two folds: the first one, approximately 80% of the data, is intended to train prediction models; the second one is aimed at validating each prediction model on an independent dataset."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"VALIDATION_SIZE = 0.2\n",
"\n",
"for time_window in TIME_WINDOWS:\n",
" in_fname = os.path.join(DATA_PATH, 'agg-dataset-{}H-cumulative.csv'.format(time_window))\n",
" if not os.path.isfile(in_fname):\n",
" continue\n",
" \n",
" window_path = {\n",
" '/': os.path.join(DATA_PATH, '{}H'.format(time_window))\n",
" }\n",
" window_path['data'] = os.path.join(window_path['/'], 'data')\n",
" \n",
" for path in window_path.values():\n",
" if not os.path.isdir(path):\n",
" os.makedirs(path)\n",
" \n",
" window_fnames = {\n",
" ds: os.path.join(window_path['data'], '{}.csv'.format(ds))\n",
" for ds in ['train', 'validation', 'test']\n",
" }\n",
" \n",
" split_data = False\n",
" for fname in window_fnames.values():\n",
" if not os.path.isfile(fname):\n",
" split_data = True\n",
" break\n",
" \n",
" if not split_data:\n",
" continue\n",
" \n",
" df = pd.read_csv(in_fname) \n",
" \n",
" df['Id'] = pd.to_numeric(df['Id'])\n",
" \n",
" columns = df.columns \n",
" \n",
" numerical_columns = [col for col in columns if re.match(r'[0-9]{2,3}H_', col)] \n",
" \n",
" categorical_columns = ['PdDistrict']\n",
" for col in columns:\n",
" if (col not in numerical_columns\n",
" and col not in ['Dataset', 'Id', 'PdDistrict']):\n",
" categorical_columns.append(col)\n",
" \n",
" if not os.path.isfile(window_fnames['test']):\n",
" test = df.loc[df['Dataset']=='test'].drop(columns=['Dataset'])\n",
" \n",
" test = test.sort_values(by=['Id'])\n",
" \n",
" columns = ['Id'] + categorical_columns + numerical_columns\n",
" test = test[columns]\n",
" \n",
" test.to_csv(window_fnames['test'], index=False)\n",
" test = None\n",
" \n",
" df = df.loc[df['Dataset']=='train'].drop(columns=['Dataset'])\n",
" \n",
" train = pd.read_csv(FNAMES['train']) \n",
" \n",
" train['Id'] = pd.to_numeric(train['Id'])\n",
" train = train[['Id', 'Category']]\n",
" \n",
" assert len(df) == len(train)\n",
" \n",
" df = pd.merge(df, train, on=['Id'], how='inner')\n",
" train = None\n",
" \n",
" columns = ['Id'] + categorical_columns + numerical_columns + ['Category']\n",
" df = df[columns]\n",
" \n",
" # The validation dataset is a stratified random sample\n",
" \n",
" validation = None\n",
" for category in CRIME_CATEGORY.keys():\n",
" category_samples = df.loc[df['Category']==category]\n",
" \n",
" sample_size = len(category_samples) * VALIDATION_SIZE\n",
" sample_size = int(np.round(sample_size, 0))\n",
" \n",
" sample = category_samples.sample(n=sample_size, replace=False, random_state=RANDOM_STATE)\n",
" \n",
" validation = (\n",
" validation.append(sample).reset_index(drop=True)\n",
" if validation is not None\n",
" else sample.copy(deep=True)\n",
" )\n",
" \n",
" validation.to_csv(window_fnames['validation'], index=False)\n",
" \n",
" df = df.loc[~df['Id'].isin(validation['Id'].values)]\n",
" df.to_csv(window_fnames['train'], index=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3.2 Baseline Model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's build a strong baseline model according to the following criteria:\n",
"\n",
"1. The training dataset will be used to learn prediction models.\n",
"2. The validation dataset will be used to rank prediction models, as they haven't seen these data during their training process.\n",
"3. The entire set of features will be used.\n",
"4. A set of several machine learning algorithms will be used to learn prediction models. Thus, a prediction model will be learned per each combination of algorithm and time window.\n",
"5. There will not be hyperparameter optimization. Machine learning algorithms will be used with their hyperparameters default values."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"def identify_attributes(\n",
" df, attributes, time_window, window_type,\n",
" attributes_to_exclude=['Id', 'PdDistrict', 'Category']\n",
" ):\n",
" \"\"\"Identify the attribute names the sets of numerical and categorical attributes consist of.\n",
" \n",
" Returns\n",
" -------\n",
" list\n",
" List of attribute names the set of numerical features consists of.\n",
" list\n",
" List of attribute names the set of categorical features consists of.\n",
" \"\"\"\n",
" columns = df.columns\n",
" \n",
" num_attr_re = re.compile(\n",
" r'{}H_'.format('[0-9]{2,3}' if window_type == 'cumulative' else time_window)\n",
" )\n",
" num_attr = (\n",
" [col for col in columns if num_attr_re.match(col)]\n",
" if attributes in ['num', 'all']\n",
" else []\n",
" )\n",
" \n",
" cat_attr = []\n",
" for col in columns:\n",
" if attributes == 'num':\n",
" break\n",
" \n",
" if (col not in attributes_to_exclude\n",
" and not re.match(r'[0-9]{2,3}H_', col)):\n",
" cat_attr.append(col)\n",
" \n",
" return num_attr, cat_attr"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"def build_model(\n",
" train, validation,\n",
" attributes, y_column,\n",
" time_window, window_type,\n",
" estimator, param_grid,\n",
" predict_proba,\n",
" return_pred=False,\n",
" return_clf=False\n",
" ):\n",
" \"\"\"Build a prediction model.\n",
" \n",
" Parameters\n",
" ----------\n",
" train : pd.DataFrame\n",
" \n",
" validation : pd.DataFrame\n",
" \n",
" attributes : str\n",
" The set of attributes to be used as input by the estimator.\n",
" \n",
" - If \"cat\", then the date-derived attributes are used.\n",
" - If \"num\", then the aggregated attributes computed using the time window are used.\n",
" - If \"all\", then the union of the attributes \"cat\" and \"raw\" is used.\n",
" \n",
" Bear in mind that, whatever the set of attributes,\n",
" the attribute \"PdDistrict\" is always used.\n",
" \n",
" y_column : str\n",
" Which is the target variable? This variable must\n",
" be in both training and validation datasets.\n",
" \n",
" time_window : int\n",
" Time window (in hours).\n",
" \n",
" window_type : str\n",
" Whether to use other windows whose time in hours is less than the specified time window value.\n",
" \n",
" - If \"exact\", then aggregated attributes that were computed using only the specified time\n",
" window, are used.\n",
" - If \"cumulative\", then aggregated attributes that were computed using windows whose time\n",
" in hours is less than the specified time window, are also used.\n",
" \n",
" estimator\n",
" A scikit-learn classification estimator.\n",
" \n",
" param_grid : dict\n",
" Dictionary of parameters, as well as their corresponding values, for the estimator.\n",
" This enables an exhaustive search over the specified parameter values through cross-validation.\n",
" \n",
" predict_proba : bool\n",
" Whether or not the estimator supports the \"predict_proba()\" method.\n",
" \n",
" return_pred : bool\n",
" \n",
" return_clf : bool\n",
" \"\"\"\n",
" num_attr, cat_attr = identify_attributes(train, attributes, time_window, window_type)\n",
" \n",
" if ('PdDistrict' in train.columns\n",
" and 'PdDistrict' not in cat_attr):\n",
" cat_attr.insert(0, 'PdDistrict')\n",
" \n",
" X_train_cat = train[cat_attr].to_numpy()\n",
" X_valid_cat = validation[cat_attr].to_numpy()\n",
" \n",
" X_train_num = None\n",
" X_valid_num = None\n",
" \n",
" scaler = StandardScaler()\n",
" if attributes in ['num', 'all']:\n",
" X_train_num = scaler.fit_transform(train[num_attr].to_numpy().astype(float))\n",
" X_valid_num = scaler.transform(validation[num_attr].to_numpy().astype(float))\n",
" \n",
" X_train = (\n",
" np.hstack([X_train_cat, X_train_num])\n",
" if X_train_num is not None\n",
" else X_train_cat\n",
" )\n",
" \n",
" X_valid = (\n",
" np.hstack([X_valid_cat, X_valid_num])\n",
" if X_valid_num is not None\n",
" else X_valid_cat\n",
" )\n",
" \n",
" y_train = train[y_column].values.astype(int)\n",
" y_valid = validation[y_column].values.astype(int)\n",
" \n",
" clf = estimator\n",
" \n",
" if isinstance(param_grid, dict):\n",
" cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)\n",
" \n",
" scoring = 'neg_log_loss' if predict_proba else 'f1_macro'\n",
" clf = GridSearchCV(\n",
" estimator=estimator, param_grid=param_grid, cv=cv,\n",
" n_jobs=5, scoring=scoring, iid=False, refit=True\n",
" )\n",
" \n",
" clf.fit(X_train, y_train)\n",
" \n",
" y_pred = clf.predict_proba(X_valid) if predict_proba else clf.predict(X_valid) \n",
" \n",
" score = log_loss(y_valid, y_pred)\n",
" \n",
" if not return_pred and not return_clf:\n",
" return score\n",
" elif return_pred and not return_clf:\n",
" return score, y_pred\n",
" elif not return_pred and return_clf:\n",
" return score, clf, scaler\n",
" else:\n",
" return score, y_pred, clf, scaler"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"def find_baseline_model(\n",
" time_windows, out_fname,\n",
" district_to_idx=DISTRICT_TO_IDX,\n",
" crime_category=CRIME_CATEGORY,\n",
" algorithms=ALGORITHMS):\n",
" \"\"\"Find the best baseline model, as specified above.\"\"\"\n",
" results = []\n",
" results_header = ['TimeWindow', 'Algorithm', 'LogLoss']\n",
" \n",
" for time_window in time_windows:\n",
" window_path = {\n",
" '/': os.path.join(DATA_PATH, '{}H'.format(time_window))\n",
" }\n",
" \n",
" window_path['data'] = os.path.join(window_path['/'], 'data')\n",
" \n",
" train = pd.read_csv(os.path.join(window_path['data'], 'train.csv'))\n",
" train = train.replace({'PdDistrict': district_to_idx, 'Category': crime_category})\n",
" \n",
" validation = pd.read_csv(os.path.join(window_path['data'], 'validation.csv'))\n",
" validation = validation.replace({'PdDistrict': district_to_idx, 'Category': crime_category})\n",
" \n",
" for algo, settings in algorithms.items():\n",
" score = build_model(\n",
" train, validation,\n",
" 'all', 'Category',\n",
" time_window, 'cumulative',\n",
" settings['estimator'], None,\n",
" settings['predict_proba']\n",
" )\n",
" results.append([str(time_window), algo, '{:.5f}'.format(score)])\n",
" \n",
" pd.DataFrame(results, columns=results_header).to_csv(out_fname, index=False, mode='w')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finding a strong baseline model takes a long time. Therefore, the set of time windows will be reduced, namely: 24, 72, and 168 hours."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"TIME_WINDOWS.pop(0)\n",
"\n",
"if not os.path.isfile(FNAMES['baseline']):\n",
" find_baseline_model(TIME_WINDOWS[:-1], FNAMES['baseline'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's analyze these results."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"baseline_results = pd.read_csv(FNAMES['baseline'])\n",
"\n",
"for col in baseline_results.columns:\n",
" if col in ['Algorithm']:\n",
" continue\n",
" \n",
" baseline_results[col] = pd.to_numeric(baseline_results[col])"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Algorithm gb logit rf\n",
"TimeWindow \n",
"24 8.13907 2.60024 14.50282\n",
"72 2.58249 2.59720 14.36502\n",
"168 10.42041 2.59499 14.36123\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEPCAYAAABhkeIdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAbjklEQVR4nO3de5zVdb3v8deb4TImNBqikEOCl5SRm4BYoYBt62AJbCQJNl4ShUdptU+WR7bHyo4dw8tp706KPdCQttF4ASUxH5U7QzZENpBDIMT2smEz3cDxiCEQw/A5f6wljOMaZs3wW/ObWfN+Ph6/B+t3Wd/fZ7Fg3vP9fX8XRQRmZmZJ6ZJ2AWZmVlwcLGZmligHi5mZJcrBYmZmiXKwmJlZorqmXUDaTjjhhBgwYEDaZZiZdSjr1q17LSL65FrX6YNlwIABrF27Nu0yzMw6FEnbmlrnQ2FmZpYoB4uZmSXKwWJmZonq9GMsudTV1VFTU8O+ffvSLqVgSktLKS8vp1u3bmmXYmZFxsGSQ01NDb169WLAgAFISrucxEUEtbW11NTUMHDgwLTLMbMi40NhOezbt4/evXsXZagASKJ3795F3SMzs/Q4WJpQrKHytmL/fGaWHgeLmZklymMsLbRs2TKmTJnC5s2bOeuss9i6dSuXXHIJGzduTKT9a6+9lhtuuIGKigpuv/12br75ZoDE99PQkB8MSaytDVdtSKwtM+uYHCwtVFlZyfnnn09lZSXf+MY3Em27vr6eBx544NB8w2Axa0pSvxj4l4K2V6y/1PlQWAvs3r2bVatW8f3vf5+HH374Xev37NnDtGnTqKioYMqUKZx33nmHbhdTWVnJkCFDGDx4MDfddNOh9/Ts2ZMvf/nLDBs2jDVr1jB+/HjWrl3L3Llz2bt3L8OHD2fmzJlAJnhmz57N2Wefzcc//nH27t0LwPjx4/nSl77EqFGjGDRoEFVVVVx66aWcccYZ3HLLLW3wN2NmdpiDpQV+/OMfM2HCBD74wQ/Su3dv1q1b94718+fP5/jjj2fTpk3cdttth9b/8Y9/5KabbuLZZ5+lurqaqqoqli1bBsBbb73Feeedx/r16zn//PMPtTVv3jyOOeYYqqurWbx4MQAvvfQS119/PS+++CLHHXccS5cuPbR99+7dWbt2LZ/97GeZPHky9957Lxs3bmTRokXU1tYW+q/GzOwQB0sLVFZWMn36dACmT59OZWXlO9avWrXq0PrBgwczdOhQAKqqqhg/fjx9+vSha9euzJw5k5UrVwJQUlLC1KlT89r/wIEDGT58OAAjR45k69ath9ZNmjQJgCFDhnD22WfTr18/evTowamnnsr27dtb/6HNzFrIYyx5ev3113n22WfZsGEDkqivr0cS119//VG1W1paSklJSV7b9ujR49DrkpKSQ4fCGq7r0qXLO7br0qULBw4cOKoazcxawj2WPC1ZsoQrrriCbdu2sXXrVrZv387AgQPf0RsYM2YMjz76KACbNm1iw4bMYNro0aN57rnneO2116ivr6eyspJx48Y1u89u3bpRV1dXmA9kZlYgDpY8VVZWMmXKlHcsmzp1Kt/61rcOzV933XXs3LmTiooKbrnlFs4++2zKysro168f8+bN48ILL2TYsGGMHDmSyZMnN7vPOXPmMHTo0EOD92ZmHYEiIu0aUjVq1Kho/KCvzZs3M2jQoBa3VV9fT11dHaWlpbzyyitcdNFFbNmyhe7duydVbqLe/pzFespjZ+HTjTuujvx/T9K6iBiVa53HWBK0Z88eLrzwQurq6ogI5s+f325DxcysUBwsCerVq5cfc2xmnZ7HWMzMLFEOFjMzS5SDxczMEuVgMTOzRHnwPg8D5v4k0fa2zvtkq987fvx47r77bkaNynmWn5lZ6txjMTOzRLnH0o7ddttt/PCHP6RPnz7079+fkSNHAvDQQw9x7bXXcuDAARYuXMjo0aNTrtTM7DAHSztVVVXF0qVLWb9+PXV1dYwYMeJQsOzZs4fq6mpWrlzJrFmzCvJUSTOz1vKhsHZq9erVTJ48mdLSUnr16sXEiRMPrZsxYwYAY8eO5c033+SNN95Iq0wzs3dxsHRAko44b2aWJgdLOzVmzBiWL1/Ovn372L17N0899dShdY888giQebBYWVkZZWVlaZVpZvYuHmPJw9GcHtxa5557LpMmTWLo0KGcdNJJDBky5FCAlJaWcs4551BXV8fChQvbvDYzsyNxsLRjX/nKV7j11lvZs2cPY8eOZeTIkcyePTvtsszMjsjB0o7NmTOHTZs2sW/fPq666ipGjBiRdklmZs1ysLRjP/rRj9IuwcysxTx4b2ZmiXKwmJlZohwsZmaWKAeLmZklyoP3+bg14QsQb93V7CY9e/Zk9+7drWr+2muv5YYbbqCiooLbb7+dm2++uVXtmJm1hnssReiBBx6goqICgNtvvz3lasyss3GwtHMRwY033sjgwYMZMmTIodu5HDx4kOuuu46zzjqLj33sY3ziE59gyZIlQOZhYGvXrmXu3Lns3buX4cOHM3PmzDQ/hpl1Ij4U1s49/vjjVFdXs379el577TXOPfdcxo4dy+rVq9m6dSubNm1ix44dDBo0iFmzZr3jvfPmzeOee+6huro6perNrDNyj6WdW7VqFTNmzKCkpISTTjqJcePGUVVVxapVq7jsssvo0qULffv25cILL0y7VDMzwMFiZmYJc7C0cxdccAGPPPII9fX17Ny5k5UrVzJ69GjGjBnD0qVLOXjwIH/5y19YsWJFzvd369aNurq6ti3azDo1j7HkI4/TgwtlypQprFmzhmHDhiGJO++8k759+zJ16lR+8YtfUFFRQf/+/RkxYkTO57LMmTOHoUOHMmLECBYvXpzCJzCzzsbB0k69fQ2LJO666y7uuuuud6zv0qULd999Nz179qS2tpbRo0czZMgQgHf0Xu644w7uuOOONqvbrL0aMPcnibWVxjOaOhIHSwd2ySWX8MYbb7B//36++tWv0rdv37RLMjMrrmCRdCrwP4GyiPhU2vUUWlPjKmZmaWr3g/eSFkraIWljo+UTJG2R9LKkuQAR8WpEXJNOpWZmBh0gWIBFwISGCySVAPcCFwMVwAxJFW1fmpmZNdbugyUiVgKvN1o8Gng520PZDzwMTM63TUlzJK2VtHbnzp0JVmtmZu0+WJpwMrC9wXwNcLKk3pK+B5wj6Z+aenNELIiIURExqk+fPoWu1cysUymqwfuIqAU+m3S7Q34wJNH2Nly1ocXveeyxx/ja175G3759+eUvf5loPWZmSeqowfIHoH+D+fLssqIUEdx///3cf//9nH/++WmXY2Z2RB31UFgVcIakgZK6A9OBJ1OuKVFbt27lzDPP5Morr6RLly4888wzXHPNNdx4441pl2ZmdkTtPlgkVQJrgDMl1Ui6JiIOAJ8HfgZsBh6NiBfTrLMQXnrpJa677joignHjxrF48eJ3XYFvZtbetPtDYRExo4nlTwNPt3E5beqUU07hQx/6UNplmJm1SLvvsXRmxx57bNolmJm1mIPFzMwS1e4PhbUHrTk92Myss3KwtFMDBgxg48bDt0fzDSfNrKPwoTAzM0tUpw0WSRMlLdi1K72nQ5qZFaNOGywRsTwi5uR6nG92fRtX1LaK/fOZWXo6bbAcSWlpKbW1tUX7wzciqK2tpbS0NO1SzKwIefA+h/LycmpqaijmW+qXlpZSXl6edhlmVoQcLDl069aNgQMHpl2GmVmH5ENhZmaWKAeLmZklysFiZmaJcrCYmVmiHCxmZpYoB4uZmSXKpxu3sQFzf5JYW1vnfTKxtszMkuIei5mZJarTBotvQmlmVhidNliauwmlmZm1TqcNFjMzKwwHi5mZJcrBYmZmiXKwmJlZohwsZmaWKAeLmZklysFiZmaJcrCYmVmiHCxmZpaoZoNF0mmSemRfj5f0RUnHFb40MzPriPLpsSwF6iWdDiwA+gM/KmhVZmbWYeUTLAcj4gAwBfhuRNwI9CtsWWZm1lHlEyx1kmYAVwFPZZd1K1xJZmbWkeUTLFcDHwb+d0T8p6SBwEOFLcvMzDqqZp8gGRGbgC8CSDoe6BURdxS6sEKTNBGYePrpp6ddiplZUcnnrLAVkt4r6X3Ab4H7JX278KUVlp/HYmZWGPkcCiuLiDeBS4F/jYjzgIsKW5aZmXVU+QRLV0n9gGkcHrw3MzPLKZ9g+V/Az4BXIqJK0qnAS4Uty8zMOqp8Bu8fAx5rMP8qMLWQRZmZWceVz+B9uaQnJO3ITksllbdFcWZm1vHkcyjsQeBJ4P3ZaXl2mZmZ2bvkEyx9IuLBiDiQnRYBfQpcl5mZdVD5BEutpMsllWSny4HaQhdmZmYdUz7BMovMqcZ/Bv4EfAr4TAFrMjOzDqzZYImIbRExKSL6RMSJEfH3+KwwMzNrQmufIHlDolWYmVnRaG2wKNEqzMysaDR7gWQTItEqEibpWGA+sB9YERGLUy7JzKzTaLLHIumvkt7MMf2VzPUszZJ0nKQlkn4vabOkD7emSEkLsxdnbsyxboKkLZJeljQ3u/hSYElEzAYmtWafZmbWOk0GS0T0ioj35ph6RUS+PZ3vAD+NiLOAYcDmhislnSipV6NluR6QsgiY0HihpBLgXuBioAKYIakCKAe2Zzerz7NWMzNLQGvHWJolqQwYC3wfICL2R8QbjTYbByyT1CP7ntnAdxu3FRErgddz7GY08HJEvBoR+4GHgclADZlwgSY+o6SJkhbs2rWrxZ/NzMyaVrBgAQYCO4EHJb0g6YHs2Mch2Rtc/gx4RNJMMtfMXNaCfZzM4Z4JZALlZOBxYKqk+8jcguZd/KAvM7PCKGSwdAVGAPdFxDnAW8DcxhtFxJ3APuA+YFJE7D7aHUfEWxFxdUR8zgP3ZmZtq5DBUgPURMTz2fklZILmHSRdAAwGngC+3sJ9/AHo32C+PLvMzMxSks9t83OdHbY9eyv9U5t6X0T8Gdgu6czsor8DNjVq+xxgAZlxkauB3pK+2YL6q4AzJA2U1B2YTuZOzGZmlpJ8zu76FzK9jx+RuTByOnAa8FtgITD+CO/9ArA4+0P/VTLh0dB7gGkR8QqApCvJcR8ySZXZ/ZwgqQb4ekR8PyIOSPo8mXGaEmBhRLyYx2cyM7MCySdYJkXEsAbzCyRVR8RNkm4+0hsjohoYdYT1qxvN1wH359huxhHaeBp4+kh1mJlZ28lnjGWPpGmSumSnaWQG26GdX4FvZmZtL59gmQlcAezITlcAl0s6Bvh8AWszM7MOqNlDYRHxKjCxidWrki3HzMw6unzOCivPngG2IzstlVTe3PvMzKxzyudQ2INkTuF9f3Zanl1mZmb2LvkES5+IeDAiDmSnRUCfAtdlZmYdVD7BUivpckkl2elyoLbQhZmZWceUT7DMAqYBfwb+BHyKHBcxmpmZQR7BEhHbImJSRPSJiBMj4u+BqW1Qm5mZdUCtvQnlDYlWYWZmRaO1waJEqzAzs6LR2mDxrVzMzCynJq+8l/RXcgeIgGMKVpGZmXVoTQZLRPRqy0LMzKw4FPIJku2apImSFuzatSvtUszMikqnDZaIWB4Rc8rKytIuxcysqHTaYDEzs8JwsJiZWaIcLGZmligHi5mZJcrBYmZmiXKwmJlZohwsZmaWKAeLmZklysFiZmaJcrCYmVmiHCxmZpYoB4uZmSXKwWJmZolysJiZWaIcLGZmligHi5mZJcrBYmZmiXKwmJlZohwsZmaWKAeLmZklysFiZmaJ6pp2AYUg6VhgPrAfWBERi1Muycys0yh4j0VSiaQXJD11FG0slLRD0sYc6yZI2iLpZUlzs4svBZZExGxgUmv3a2ZmLdcWh8L+Edica4WkEyX1arTs9BybLgIm5Hh/CXAvcDFQAcyQVAGUA9uzm9W3unIzM2uxggaLpHLgk8ADTWwyDlgmqUd2+9nAdxtvFBErgddzvH808HJEvBoR+4GHgclADZlwgSY+o6SJkhbs2rWrBZ/IzMyaU+gey78A/wM4mGtlRDwG/Ax4RNJMYBZwWQvaP5nDPRPIBMrJwOPAVEn3Acub2PfyiJhTVlbWgt2ZmVlzCjZ4L+kSYEdErJM0vqntIuJOSQ8D9wGnRcTuo913RLwFXH207ZiZWcsVsscyBpgkaSuZQ1QflfTDxhtJugAYDDwBfL2F+/gD0L/BfHl2mZmZpaRgwRIR/xQR5RExAJgOPBsRlzfcRtI5wAIy4yJXA70lfbMFu6kCzpA0UFL37H6eTOQDmJlZq6R9geR7gGkR8UpEHASuBLY13khSJbAGOFNSjaRrACLiAPB5MuM0m4FHI+LFNqvezMzepU0ukIyIFcCKHMtXN5qvA+7Psd2MI7T9NPD0URdpZmaJSLvHYmZmRcbBYmZmiXKwmJlZohwsZmaWKAeLmZklysFiZmaJcrCYmVmiHCxmZpYoB4uZmSXKwWJmZolysJiZWaIcLGZmligHi5mZJcrBYmZmiXKwmJlZohwsZmaWKAeLmZklysFiZmaJcrCYmVmiHCxmZpYoB4uZmSXKwWJmZonqmnYBhSDpWGA+sB9YERGLUy7JzKzTKFiPRVKppN9IWi/pRUnfOIq2FkraIWljjnUTJG2R9LKkudnFlwJLImI2MKm1+zUzs5Yr5KGwvwEfjYhhwHBggqQPNdxA0omSejVadnqOthYBExovlFQC3AtcDFQAMyRVAOXA9uxm9Uf5OczMrAUKFiyRsTs72y07RaPNxgHLJPUAkDQb+G6OtlYCr+fYzWjg5Yh4NSL2Aw8Dk4EaMuECHkcyM2tTBf2hK6lEUjWwA3gmIp5vuD4iHgN+BjwiaSYwC7isBbs4mcM9E8gEysnA48BUSfcBy5uobaKkBbt27WrB7szMrDkFDZaIqI+I4WR6D6MlDc6xzZ3APuA+YFKDXs7R7PetiLg6Ij7X1MB9RCyPiDllZWVHuzszM2ugTQ4TRcQbwC/JPU5yATAYeAL4egub/gPQv8F8eXaZmZmlpJBnhfWRdFz29THAx4DfN9rmHGABmXGRq4Hekr7Zgt1UAWdIGiipOzAdeDKJ+s3MrHUK2WPpB/xS0u/IBMAzEfFUo23eA0yLiFci4iBwJbCtcUOSKoE1wJmSaiRdAxARB4DPkxmn2Qw8GhEvFuwTmZlZswp2gWRE/A44p5ltVjearwPuz7HdjCO08TTwdCvLNDOzhPlUXDMzS5SDxczMEuVgMTOzRDlYzMwsUQ4WMzNLlIPFzMwS5WAxM7NEOVjMzCxRDhYzM0uUg8XMzBJVlM+87zRuTeiW/wM/kEw7lr+kvjvw95cG/987IgeLWQsMmPuTRNrZWppIM2btkg+FmZlZohwsZmaWKAeLmZklysFiZmaJcrCYmVmiHCxmZpYoB4uZmSXKwWJmZolysJiZWaIUEWnXkCpJO4FtaddRQCcAr6VdhLWKv7uOrdi/v1Miok+uFZ0+WIqdpLURMSrtOqzl/N11bJ35+/OhMDMzS5SDxczMEuVgKX4L0i7AWs3fXcfWab8/j7GYmVmi3GMxM7NEOVjMzCxRDhYzM0uUg6XISTox7RrMrHNxsBQRSe9rNPUGfiPpeEnvS7s+a5qk90r6lqSHJP1Do3Xz06rLrDV8VlgRkXSQd9+ephyoASIiTm37qiwfkpYCLwG/BmYBdcA/RMTfJP02IkakWqA1S9IHgDcj4g1JA4BRwO8jYmOqhaXAPZbiciOwBZgUEQMjYiBQk33tUGnfTouIuRGxLCImAb8Fns32Oq2dkzQXeA74taRrgZ8CFwOPSLoh1eJS4B5LkZFUDvwzsB34OrDeodL+SdoMnB0RBxss+wyZXxZ6RsQpadVmzZP0IpkeynuArcCpEbFT0rHA8xExOM362pp7LEUmImoi4jJgBfAMmX/o1v4tBz7acEFELAK+DOxPoyBrkfqI2Au8AewFagEi4q1Uq0qJeyxFRtJZwMnA80A9mUMsGyVNiIifpludNUXSF4EnImJ72rVYy0laBHQHjgX2AAfIHA77KNArIqalV13bc7AUkewPp+uBzcBw4B8j4sfZdR4Absck7QLeAl4BKoHHImJnulVZviR1BS4DAlgCnAfMAP4LuLez9VwcLEVE0gbgwxGxO3tWyhLgoYj4jqQXIuKcVAu0Jkl6ARgJXAR8GpgErCMTMo9HxF9TLM9aQVLviKhNu440eIyluHSJiN0AEbEVGA9cLOnbgFKsy5oXEXEwIn4eEdcA7wfmAxOAV9MtzZojaZ6kE7KvR0l6lcwZYtskjUu5vDbnYCkuf5E0/O2ZbMhcQuYRqUNSq8ry8Y7gj4i6iHgyImYAPiOs/ftkRLz9GOK7gE9HxBnAx4D/k15Z6XCwFJcrgT83XBARByLiSmBsOiVZnj7d1IqI2NOWhVirdM2OswAcExFVABHxH0CP9MpKh8dYzMyOkqQvABOBeWR+iTseeJzMWWGnRsQVKZbX5hwsZmYJkDQe+BzwQaArmYuUlwELI+JAiqW1OQeLmVkBSbo6Ih5Mu4625GAxMysgSf8VER9Iu4621LX5TczM7Egk/a6pVcBJbVlLe+BgMTM7eicB/w34f42WC/hV25eTLgeLmdnRe4rMXairG6+QtKLty0mXx1jMzCxRvkDSzMwS5WAxM7NEOVjMcpDUW1J1dvqzpD80mE98MFbSC2/f501SV0m7JV3eYP06SSMkTco+BrclbS+S9KmkazZrigfvzXLI3u787R/0twK7I+LuAu5yNfARoBoYBvxHdv6H2cfbnkbmMdO/BZ4sYB1mR809FrMWkrQ7++d4Sc9J+rGkV7O3Tp8p6TeSNkg6LbtdH0lLJVVlpzE5mv0VmSAh++f3yAYbMBpYFxH1kj4j6Z5su4sk/V9Jv8ru/1PZ5ZJ0j6Qtkv4NOLFB7X+X7R1tkLRQUg9J50p6PLt+sqS9krpLKs3e/t2sRRwsZkdnGPBZYBBwBfDBiBgNPAB8IbvNd4B/johzganZdY293WMh++dK4G+SemXnmzr81g84n8zjEeZll00BzgQqyNzx+iMAkkqBRWRu6T6EzBGLzwEvcDjELgA2AueSeQri8/n9NZgd5mAxOzpVEfGniPgbmccK/zy7fAMwIPv6IuAeSdVkDmO9V1LPho1ExDagu6S+wFnAFqCKzA/3j5AJnlyWZR8QtonDV3iPBSojoj4i/gg8m11+JvCf2Vu5A/wAGJu9QeIrkgaR6R19O9vGBcC/t/hvxDo9j7GYHZ2/NXh9sMH8QQ7//+oCfCgi9jXT1q/IPDf9TxERkn4NjCHzw35NHvs/mqeErgQuBuqAfyPTsykBbjyKNq2Tco/FrPB+zuHDYjR8ymcjvwL+O4dDZA3Zh7dFxK4W7G8l8GlJJZL6ARdml28BBkg6PTt/BfBc9vW/v73viNgJ9CbTw9nYgv2aAQ4Ws7bwRWCUpN9J2kRmTCaX1cCpZIMlIv5EptfQ0tObnwBeAjYB/9qgvX3A1cBjkjaQ6VV9L/ue58kcSluZnf8dsCF8aw5rBd/SxczMEuUei5mZJcrBYmZmiXKwmJlZohwsZmaWKAeLmZklysFiZmaJcrCYmVmi/j+cICSSMoh5XQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"chart_data = baseline_results.pivot(index='TimeWindow', columns='Algorithm', values='LogLoss')\n",
"\n",
"print(chart_data)\n",
"\n",
"chart_data.plot(kind='bar', grid=False, legend=True, x=None)\n",
"\n",
"plt.xlabel('Time Window')\n",
"plt.ylabel('Log Loss')\n",
"\n",
"plt.yscale(\"log\")\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notably, the Logistic Regression algorithm outperforms all of the machine learning algorithms, even though the best result, i.e., **2.58249**, was achieved by the Gradient Boosting algorithm on data aggregated using a time window of 72 hours. This conclusion is drawn as Logistic Regression is the most computationally efficient algorithm, and the difference between its best result and the overall best result is negligible."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On the other hand, let's analyze how having a larger time window contributes to crime classification. The Logistic Regression algorithm will be used to conduct this analysis."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"FNAMES['logit-baseline'] = os.path.join(DATA_PATH, 'baseline-logit-models.csv')\n",
"\n",
"if not os.path.isfile(FNAMES['logit-baseline']):\n",
" logit_baseline = baseline_results.loc[baseline_results['Algorithm']=='logit']\n",
" logit_baseline = logit_baseline.drop(columns=['Algorithm'])\n",
" \n",
" time_window = 336\n",
" \n",
" window_data_path = os.path.join(DATA_PATH, '{}H'.format(time_window), 'data')\n",
" \n",
" train = pd.read_csv(os.path.join(window_data_path, 'train.csv'))\n",
" train = train.replace({'PdDistrict': DISTRICT_TO_IDX, 'Category': CRIME_CATEGORY})\n",
" \n",
" validation = pd.read_csv(os.path.join(window_path['data'], 'validation.csv'))\n",
" validation = validation.replace({'PdDistrict': DISTRICT_TO_IDX, 'Category': CRIME_CATEGORY})\n",
" \n",
" algo = ALGORITHMS['logit']\n",
" \n",
" score = build_model(\n",
" train, validation,\n",
" 'all', 'Category',\n",
" time_window, 'cumulative',\n",
" algo['estimator'], None,\n",
" algo['predict_proba']\n",
" )\n",
" \n",
" logit_baseline = logit_baseline.append({'TimeWindow': time_window, 'LogLoss': score}, ignore_index=True)\n",
" \n",
" logit_baseline.to_csv(FNAMES['logit-baseline'], index=False, mode='w')"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" TimeWindow LogLoss\n",
"0 24 2.600240\n",
"1 72 2.597200\n",
"2 168 2.594990\n",
"3 336 2.594697\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAAEPCAYAAAAeQPDsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAe2UlEQVR4nO3df7RV5X3n8fcnoiagISkRJ4IORK2RxQwSKTUxaGLSFEnENkYj448ULa52aaITV2cxa80sZ7rSTrQdO45VZjAq/YE/oqCVylLSFZVqLKKGRAS1ajVCcTQ/RMUkCHzmj71vOZ6ewz3ncu89PPJ5rXXWPWc/z7P395wF97OevZ+7j2wTERGxp3tPrwuIiIjoRAIrIiKKkMCKiIgiJLAiIqIICayIiCjCiF4X8G71oQ99yBMmTOh1GRERRXnsscd+bPugVm0JrCEyYcIEHn300V6XERFRFEkvtmvLKcGIiChCAisiIoqQwIqIiCIksCIioggJrIiIKEICKyIiipDAioiIIiSwIiKiCPnD4YJMmH93r0voyAvf/HyvS4iId6HMsCIioggJrIiIKEICKyIiipDAioiIIiSwIiKiCFkl2CVJo4Brga3A/bYX97ikiIi9wh41w5J0qKT7JK2T9KSki9v0+4Ck2yU9JWm9pI/vxjFvkPSKpLVN22dKelrSs5LmNzR9Ebjd9jxg9kCPGxER3dmjAgvYBlxqexJwHHChpEkt+l0F3GP7o8AUYH1jo6Sxkg5s2nZEm2MuAmY29d0HuAY4GZgEzGmoYzzwUv18e4fvKyIidtMeFVi2N9l+vH7+BlUQjWvsI2k0cAJwfd1vq+3XmnZ1InCnpP3rMfOAq9sccyXw06bN04FnbT9veytwC3Bq3baBKrSgxecn6RRJCzdv3tzBO46IiE7tUYHVSNIEYCqwqqlpIvAqcKOk70v6Vn1d6V/Yvg24F7hV0lnAecDpXRx+HDtnUVCFVF9wLgVOk7QAWNY80PYy2xeMHj26i8NFRER/9sjAknQAsAS4xPbrTc0jgI8BC2xPBbYA85v6YPsK4BfAAmC27TcHozbbW2zPtf37WXARETF89rjAkrQvVVgttr20RZcNwAbbfTOv26kCrHk/M4DJwB3AZV2WsRE4tOH1+HpbRET0yB4VWJJEdW1qve0rW/Wx/TLwkqSj6k2fAdY17WcqsJDqutNcYIykb3RRymrgSEkTJe0HnAnc1dWbiYiIQbVHBRZwPHAOcJKkNfVjFoCk5ZIOqft9FVgs6YfAMcAfN+1nJHCG7eds7wDOBV5sdUBJNwMPA0dJ2iDpfNvbgIuoroOtB75t+8nBfasREdGNPeoPh20/CKhN26yG52uAabvYz0NNr98GrmvTd06b7cuB5f1XHRERw2FPm2FFRES0lMCKiIgiJLAiIqIICayIiChCAisiIoqQwIqIiCIksCIioggJrIiIKEICKyIiipDAioiIIiSwIiKiCAmsiIgoQgIrIiKKkMCKiIgiJLAiIqIICayIiChCAisiIoqQwIqIiCIksCIioggJrIiIKEICKyIiipDAioiIIiSwIiKiCAmsiIgoQgIrIiKKMKLXBZRE0ijgWmArcL/txT0uKSJirzFkMyxJh0q6T9I6SU9KurhNvxckPSFpjaRHG7ZfLGltPfaSpjFt27qs8QZJr0ha26JtpqSnJT0raX69+YvA7bbnAbMHetyIiOjeUJ4S3AZcansScBxwoaRJbfp+2vYxtqcBSJoMzAOmA1OAL0g6or+2PpLGSjqwads7+tQWATObN0raB7gGOBmYBMypax8PvFR3277rtx8REYNpyALL9ibbj9fP3wDWA+M6HH40sMr2W7a3AQ9QzW76a+tzInCnpP0BJM0Drm5R40rgpy2OPx141vbztrcCtwCnAhuoQgvafHaSTpG0cPPmzR2+1YiI6MSwLLqQNAGYCqxq0WxghaTHJF1Qb1sLzJA0RtJIYBZwaAdt1Q7t24B7gVslnQWcB5zeRcnj2DmTgiqoxgFLgdMkLQCWtRpoe5ntC0aPHt3F4SIioj9DvuhC0gHAEuAS26+36PJJ2xsljQW+I+kp2yslXQ6sALYAa6hPwdle366tke0rJN0CLAAOt/3m7r4X21uAubu7n4iI6N6QzrAk7UsVVottL23Vx/bG+ucrwB1Up+Owfb3tY22fAPwMeKZhTNu2hmPPACbX+7ysy9I38s5Z2/h6W0RE9MhQrhIUcD2w3vaVbfqM6lscUS8Z/xzVKT/qGReSDqO6RnVTw7i2bfX2qcBCqutOc4Exkr7RRfmrgSMlTZS0H3AmcFcX4yMiYpAN5QzreOAc4KR6yfoaSbMAJC2XdAhwMPCgpB8AjwB3276nHr9E0jqqa0UX2n6tYd+7agMYCZxh+znbO4BzgRebC5R0M/AwcJSkDZLOB6gXc1xEdR1sPfBt20/u/kcSEREDNWTXsGw/CKhN26yGl1Pa9Jmxi323bavbH2p6/TZwXYt+c3axj+XA8l0dJyIihk9uzRQREUVIYEVERBESWBERUYQEVkREFCGBFRERRUhgRUREERJYERFRhARWREQUIYEVERFFSGBFREQRElgREVGEBFZERBQhgRUREUVIYEVERBESWBERUYQEVkREFCGBFRERRUhgRUREERJYERFRhARWREQUIYEVERFFSGBFREQRElgREVGEBFZERBQhgRUREUVIYEVERBFG9LqAkkgaBVwLbAXut724xyVFROw1hmyGJelQSfdJWifpSUkXt+n3gqQnJK2R9GjD9oslra3HXtI05j/W29dKulnSewdY4w2SXpG0tkXbTElPS3pW0vx68xeB223PA2YP5JgRETEw/QaWpMMl7V8//5Skr0n6QAf73gZcansScBxwoaRJbfp+2vYxtqfVx5kMzAOmA1OAL0g6om4bB3wNmGZ7MrAPcGZTzWMlHdi07YgWx10EzGzxnvcBrgFOBiYBc+raxwMv1d227/rtR0TEYOpkhrUE2F7/wl8IHArc1N8g25tsP14/fwNYD4zrsK6jgVW237K9DXiAanbTZwTwPkkjgJHAPzeNPxG4syFo5wFXt6hxJfDTFsefDjxr+3nbW4FbgFOBDVShBW0+O0mnSFq4efPmzt5pRER0pJPA2lGHxm8DV9v+A+DD3RxE0gRgKrCqRbOBFZIek3RBvW0tMEPSGEkjgVlUQYntjcCfAj8CNgGbba94xw7t24B7gVslnQWcB5zeRcnj2DmTgiqoxgFLgdMkLQCWtRpoe5ntC0aPHt3F4SIioj+dLLp4W9Ic4CvAKfW2fTs9gKQDqGZpl9h+vUWXT9reKGks8B1JT9leKelyYAWwBVhDfQpO0gepZjsTgdeA2ySdbfuvG3dq+wpJtwALgMNtv9lpze3Y3gLM3d39RERE9zqZYc0FPg78ke1/kjQR+KtOdi5pX6qwWmx7aas+9YwJ268Ad1CdjsP29baPtX0C8DPgmXrIZ4F/sv2q7bepZj2faHHsGcDkep+XdVJvg43UM7ra+HpbRET0SL+BZXud7a/Zvrme3Rxo+/L+xkkScD2w3vaVbfqM6lscUS8Z/xzV6UDqGReSDqO6ftV33exHwHGSRtbH+AzV9bHG/U6lut52KlXgjpH0jf5qbrAaOFLSREn7US3quKuL8RERMcg6WSV4v6T3S/oV4HHgOkktA6jJ8cA5wEn1kvU1kmbV+1wu6RDgYOBBST8AHgHutn1PPX6JpHVU14outP0agO1VwO11LU/U72Fh07FHAmfYfs72DuBc4MUW7+1m4GHgKEkbJJ1fH2MbcBHVdbD1wLdtP9nBe46IiCHSyTWs0bZfl/S7wF/avkzSD/sbZPtBQG3aZjW8nNKmz4xd7PsydnGaz/ZDTa/fBq5r0W/OLvaxHFjerj0iIoZXJ9ewRkj6MHAG8LdDXE9ERERLnQTWH1KdGnvO9mpJHwH+cWjLioiIeKd+TwnWf9N0W8Pr54HThrKoiIiIZp0suhgv6Y76nnuvSFoiaXx/4yIiIgZTJ6cEb6Ra0n1I/VhWb4uIiBg2nQTWQbZvtL2tfiwCDhriuiIiIt6hk8D6iaSzJe1TP84GfjLUhUVERDTqJLDOo1rS/jLVzWa/BPzOENYUERHxr3Rya6YXbc+2fZDtsbZ/i6wSjIiIYTbQbxz++qBWERER0Y+BBlbLWy5FREQMlYEGlge1ioiIiH60vdOFpDdoHUwC3jdkFUVERLTQNrBsHzichUREROzKQE8JRkREDKsEVkREFCGBFRERRUhgRUREEfr9Pqw2qwU3A48Cl9bfjxURETGk+g0s4H8BG4CbqJa0nwkcDjwO3AB8aqiKi4iI6NPJKcHZtv+v7Tdsv257IfCbtm8FPjjE9UVERACdBdZbks6Q9J76cQbwi7otd7yIiIhh0ckpwbOAq4Br69cPA2dLeh9w0VAVFjHUJsy/u9cldOSFb36+1yVE7BH6Dax6UcUpbZofHNxyIiIiWuv3lKCk8ZLukPRK/VgiafxwFBcREdGnk2tYNwJ3AYfUj2X1tr2OpFGS/kLSdZLO6nU9ERF7k04C6yDbN9reVj8WAQf1N0jSoZLuk7RO0pOSLm7T7wVJT0haI+nRhu0XS1pbj72kYftRdd++x+uN7d2QdEM9a1zbom2mpKclPStpfr35i8DttucBswdyzIiIGJhOAusnks6WtE/9OBv4SQfjtlH9YfEk4DjgQkmT2vT9tO1jbE8DkDQZmAdMB6YAX5B0BIDtp+u+xwDHAm8BdzTuTNJYSQc2bTuixXEXATObN0raB7gGOBmYBMypax8PvFR3297P+4+IiEHUSWCdB5wBvAxsAr4E/E5/g2xvsv14/fwNYD0wrsO6jgZW2X7L9jbgAarZTbPPAM/ZfrFp+4nAnZL2B5A0D7i6RY0rgZ+22O904Fnbz9veCtwCnEr1B9R91+9yW6uIiGHU7y9d2y/anm37INtjbf8WcFo3B5E0AZgKrGp1CGCFpMckXVBvWwvMkDRG0khgFnBoi7FnAje3qPk24F7g1vpa03nA6V2UPI6dMymogmocsBQ4TdICqmt5/4qkUyQt3Lx5cxeHi4iI/nTyd1itfJ3qlk39knQAsAS4xPbrLbp80vZGSWOB70h6yvZKSZcDK4AtwBqaTsFJ2o/qOtJ/bnVc21dIugVYABxu+80O31tbtrcAc/vpswxYNm3atHm7e7yIiNhpoKe11FEnaV+qsFpse2mrPrY31j9foboWNb1+fb3tY22fAPwMeKZp6MnA47b/X5tjzwAm1/u8rJN6G2zknTO68fW2iIjokYEGVr+3ZJIk4Hpgve0r2/QZ1bc4QtIo4HNUpwOpZ1xIOozq+tVNTcPn0OJ0YD1mKrCQ6rrTXGCMpG/0/7b+xWrgSEkT65ncmVRL+yMiokfaBpakN+ol482PN6j+Hqs/xwPnACc1LEGfVe97uaRDgIOBByX9AHgEuNv2PfX4JZLWUV0rutD2aw21jQJ+g+qaUisjgTNsP2d7B3Au0LwwA0k3U91q6ihJGySdD1Av9LiI6jrYeuDbtp/s4D1HRMQQaXsNy/aB7do6YftB2pw6tD2r4eWUNn1m7GLfW4Axu2h/qOn128B1LfrN2cU+lgPL27VHRMTwytLsiIgoQgIrIiKKkMCKiIgiJLAiIqIICayIiChCAisiIoqQwIqIiCIksCIioggJrIiIKEICKyIiipDAioiIIiSwIiKiCAmsiIgoQgIrIiKKkMCKiIgiJLAiIqIICayIiChCAisiIoqQwIqIiCIksCIioggJrIiIKEICKyIiipDAioiIIiSwIiKiCAmsiIgoQgIrIiKKMKLXBZRE0ijgWmArcL/txT0uKSJirzFkMyxJh0q6T9I6SU9KurhNvxckPSFpjaRHG7ZfLGltPfaSpjEfkHS7pKckrZf08QHWeIOkVyStbdE2U9LTkp6VNL/e/EXgdtvzgNkDOWZERAzMUJ4S3AZcansScBxwoaRJbfp+2vYxtqcBSJoMzAOmA1OAL0g6oqH/VcA9tj9at69v3JmksZIObNrWOL7PImBm80ZJ+wDXACcDk4A5de3jgZfqbtvbvfGIiBh8QxZYtjfZfrx+/gZVqIzrcPjRwCrbb9neBjxANbtB0mjgBOD6et9bbb/WNP5E4E5J+9dj5gFXt6hxJfDTFsefDjxr+3nbW4FbgFOBDVShBW0+O0mnSFq4efPmDt9qRER0YlgWXUiaAEwFVrVoNrBC0mOSLqi3rQVmSBojaSQwCzi0bpsIvArcKOn7kr5VX1vauUP7NuBe4FZJZwHnAad3UfI4ds6koAqqccBS4DRJC4BlrQbaXmb7gtGjR3dxuIiI6M+QL7qQdACwBLjE9ustunzS9kZJY4HvSHrK9kpJlwMrgC3AGnaeghsBfAz4qu1Vkq4C5gP/tXGntq+QdAuwADjc9pu7+15sbwHm7u5+IiKie0M6w5K0L1VYLba9tFUf2xvrn68Ad1CdjsP29baPtX0C8DPgmXrIBmCD7b7Z2u1UAdZ87BnA5Hqfl3VZ+kZ2zuigOg24sct9RETEIBrKVYKius603vaVbfqM6lscUZ/W+xzV6UDqGReSDqO6fnUTgO2XgZckHVXv5jPAuqb9TgUWUl13mguMkfSNLspfDRwpaaKk/YAzgbu6GB8REYNsKGdYxwPnACfVS9bXSJoFIGm5pEOAg4EHJf0AeAS42/Y99fglktZRXSu6sGlhxVeBxZJ+CBwD/HHTsUcCZ9h+zvYO4FzgxeYCJd0MPAwcJWmDpPMB6oUeF1FdB1sPfNv2k7v9iURExIAN2TUs2w8CatM2q+HllDZ9Zuxi32uAabtof6jp9dvAdS36zdnFPpYDy9u1R0TE8MqtmSIioggJrIiIKEICKyIiipDAioiIIiSwIiKiCAmsiIgoQgIrIiKKkMCKiIgiJLAiIqIICayIiChCAisiIoqQwIqIiCIksCIioggJrIiIKEICKyIiipDAioiIIiSwIiKiCAmsiIgoQgIrIiKKkMCKiIgiJLAiIqIICayIiChCAisiIoqQwIqIiCIksCIioggJrIiIKMKIXhdQEkmjgGuBrcD9thf3uKSIiL3GkM2wJB0q6T5J6yQ9KeniNv1ekPSEpDWSHm3YfrGktfXYSzoZM4Aab5D0iqS1LdpmSnpa0rOS5tebvwjcbnseMHugx42IiO4N5SnBbcClticBxwEXSprUpu+nbR9jexqApMnAPGA6MAX4gqQjdjWmkaSxkg5s2tY8HmARMLPF+H2Aa4CTgUnAnLr28cBLdbftbd5LREQMgSE7JWh7E7Cpfv6GpPXAOGBdB8OPBlbZfgtA0gNUs5srOjz8icDvSZpl+5eS5tXjT26qcaWkCS3GTweetf18ffxbgFOBDVShtYY2YS/pFOCUI45olY8R714T5t/d6xI68sI3P9/rEmKAhmXRRR0KU4FVLZoNrJD0mKQL6m1rgRmSxkgaCcwCDu1nzM5G+zbgXuBWSWcB5wGnd1HyOHbOpKAKqnHAUuA0SQuAZa0G2l5m+4LRo0d3cbiIiOjPkC+6kHQAsAS4xPbrLbp80vZGSWOB70h6qp75XA6sALZQzWi29zemcae2r6hnRguAw22/ubvvxfYWYO7u7iciIro3pDMsSftShdVi20tb9bG9sf75CnAH1ek4bF9v+1jbJwA/A57pb0zTsWcAk+v2y7osfSPvnNGNr7dFRESPDOUqQQHXA+ttX9mmz6i+xRH1kvHPUZ0OpJ49IekwqutPN/U3pmG/U4GFVNed5gJjJH2ji/JXA0dKmihpP+BM4K4uxkdExCAbyhnW8cA5wEn18vM1kmYBSFou6RDgYOBBST8AHgHutn1PPX6JpHVU14outP1avX1XY/qMBM6w/ZztHcC5wIvNBUq6GXgYOErSBknnA9jeBlxEdR1sPfBt208OyqcSEREDMpSrBB8E1KZtVsPLKW36zGiz/fl2Yxr6PNT0+m3guhb95uxiH8uB5bs6TkREDJ/cmikiIoqQwIqIiCIksCIioggJrIiIKELu1h4RsYfJba5aywwrIiKKkMCKiIgiJLAiIqIICayIiChCAisiIoqQwIqIiCIksCIioggJrIiIKEICKyIiiiDbva7hXUnSq7T4Dq490IeAH/e6iHeRfJ6DK5/n4Cnls/y3tg9q1ZDA2stJetT2tF7X8W6Rz3Nw5fMcPO+GzzKnBCMioggJrIiIKEICKxb2uoB3mXyegyuf5+Ap/rPMNayIiChCZlgREVGEBFZERBQhgRUREUVIYO3lJI3tdQ0REZ1IYO1FJP1K02MM8IikD0r6lV7XVxpJ75f0PyT9laT/0NR2ba/qimhW/1s9VtIHe13L7khg7V1+DDzW8HgUGAc8Xj+P7twICFgCnClpiaT967bjeldWuSQdJukD9fMJkr4kaXKv6yqNpL+W9KH6+W8Ca4HLgTWSTu9pcbshgbV3+QPgaWC27Ym2JwIb6ucf6XFtJTrc9nzbd9qeTRX8361nrtElSfOBB4B/kPS7wD3AycCtkr7e0+LKM8V2330DLwNOsP1Z4Fjgv/SurN0zotcFxPCx/T8l3Qr8maSXqP4h5w/xBm5/Se+xvQPA9h9J2gisBA7obWlFOgeYBIwEXgA+YvtVSaOAVcCVPaytNO+R9H7brwM7gB8B2P6xpGJ/72eGtZexvcH26cD9wHeofjnEwCwDTmrcYHsRcCmwtRcFFW677Z8DrwE/B34CYHtLT6sq038H7pN0HvAQcJukr0haRDVzLVLudLGXkfRRqutWq4DtVKe11kqaabvYf8i9IOlrwB22X+p1Le8G9S/T/YBRwFvANqpfricBB9o+o3fVlUfSkcDvAr9KdTZtA3Cn7Xt7WthuSGDtRepfsBcC64FjgItt/03d9rjtj/WyvtJI2gxsAZ4DbgZus/1qb6sqV32q6nSq09S3A78OzKE6nXVNZlqRwNqLSHoC+LjtNyVNoPql8Fe2r5L0fdtTe1pgYSR9n+oi9meBLwOzqVZf3gwstf1GD8t7V5A0xvZPel1HaSSNBC6iCv+rqf59ngY8Bfyh7Td7WN6A5RrW3uU9ff9Qbb8AfAo4WdKVVMuzozu2vcP2CtvnA4cA1wIzged7W1p5JH2zYSn2NEnPU60YfFHSiT0urzSLgIOBicDdwK8Bf0L1/3xB78raPZlh7UUkfRf4uu01DdtGADcAZ9nep2fFFWhXs1JJI22/Ndw1lUzSE7b/Xf38PuA/2V4t6VeBm0r/ttzhJGmN7WMkCdgEfNi269c/sP3ve1zigGSGtXc5F3i5cYPtbbbPBU7oTUlF+3K7hoTVgIxoWHL9PturAWw/A+zffli042pGsrz+2fe62FlKsevxo3u2N+yi7aHhrOXdoP5FGoPnWmC5pG8C90i6ClhKtUpwzS5HRrNHJR1g+03b5/VtlHQ4UOy11ZwSjIg9hqRPAb/PzqXYLwF3AjfY3tbD0oojaTrVpGq1pElU11afpmHGVZoEVkTs8STNtX1jr+sohaTLqG5rNYLqBgG/DtwH/AZwr+0/6mF5A5bAiog9nqQf2T6s13WUov4TlmOorv29DIy3/bqk9wGrSl10kWtYEbFHkPTDdk1US7Sjc9tsbwfekvRcfU9BbP9c0o4e1zZgCayI2FMcDPwm8LOm7QK+N/zlFG1rw59WHNu3UdJoqpvhFimBFRF7ir8FDmj8O8E+ku4f/nKKdoLtXwL0fZtAbV/gK70pafflGlZERBQhfzgcERFFSGBFREQRElgRw0jSGElr6sfLkjY2vB70hQWSvi/pmPr5CElvSjq7of0xSR+TNLv+ivpu9r1I0pcGu+aIdrLoImIY1V+V0Rcg/w140/afDuEhHwI+QXVroynAM/Xrv66/ev5wqpuhPg7cNYR1ROy2zLAi9hCS3qx/fkrSA5L+RtLz9ddunCXpEUlP1PeDQ9JBkpZIWl0/jm+x2+9RBRT1z/9DHZjAdOAx29sl/Y6kP6/3u0jS/5b0vfr4X6q3S9KfS3pa0t8BYxtq/0w9m3tC0g2S9pf0a5KW1u2nSvq5pP0kvbf+6pCIriSwIvZMU4DfA44GzgF+1fZ04FvAV+s+VwF/ZvvXqL6c71st9tM3w6L+uRL4paQD69ftTkN+GPgk8AXgm/W23waOAiZR3fn/EwCS3kv1/Utfrr8eZATV/QC/z85wnAGspfpepl8HVnX2MUTslMCK2DOttr2p/lua54AV9fYngAn1888Cfy5pDdXpvPdLOqBxJ7ZfBPaT9G+Aj1Ld/HQ1VWh8girQWrmz/nLKdey8y8QJwM22t9v+Z+C79fajgH9quHv9X1D9HdA24DlJR1PN5q6s9zED+PuuP5HY6+UaVsSe6ZcNz3c0vN7Bzv+37wGOs/2Lfvb1PeB0YFP9JX7/ABxPFSIPd3D83fk26pVUN2F9G/g7qpnYPsAf7MY+Yy+VGVZEuVaw8/QgfasBW/gecAk7w+lh6i/ztL25i+OtBL4saR9JHwY+XW9/Gpgg6Yj69TnAA/Xzv+87tu1XgTFUM7K1XRw3AkhgRZTsa8A0ST+UtI7qmlcrDwEfoQ4s25uoZjndLqO/A/hHYB3wlw37+wUwF7itvkv4DqrFHVBdqzqYKuwAfgg8Uer3MUVv5dZMERFRhMywIiKiCAmsiIgoQgIrIiKKkMCKiIgiJLAiIqIICayIiChCAisiIorw/wFjRSvdnfFVQAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"logit_baseline = pd.read_csv(FNAMES['logit-baseline'])\n",
"\n",
"logit_baseline = logit_baseline.astype({'TimeWindow': 'int32', 'LogLoss': 'float64'})\n",
"\n",
"print(logit_baseline)\n",
"\n",
"logit_baseline.plot(x='TimeWindow', y='LogLoss', kind='bar', grid=False, legend=False)\n",
"\n",
"plt.xlabel('Time Window')\n",
"plt.ylabel('Log Loss')\n",
"\n",
"plt.yscale(\"log\")\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To sum up, the **baseline result** is **2.5947**."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3.3 Discriminative Features"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The goal is to quantify how discriminative each set of features is, as well as each window type."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"def discriminative_feature(\n",
" time_window, algorithm, out_fname,\n",
" district_to_idx=DISTRICT_TO_IDX,\n",
" crime_category=CRIME_CATEGORY):\n",
" \"\"\"Quantify how discriminative each set of features is, as well as each window type.\"\"\"\n",
" window_data_path = os.path.join(DATA_PATH, '{}H'.format(time_window), 'data')\n",
" \n",
" train = pd.read_csv(os.path.join(window_data_path, 'train.csv'))\n",
" train = train.replace({'PdDistrict': district_to_idx, 'Category': crime_category})\n",
" \n",
" validation = pd.read_csv(os.path.join(window_data_path, 'validation.csv'))\n",
" validation = validation.replace({'PdDistrict': district_to_idx, 'Category': crime_category})\n",
" \n",
" grid = {\n",
" 'attributes': ['all', 'num', 'cat'],\n",
" 'window_type': ['exact', 'cumulative']\n",
" }\n",
" \n",
" results = []\n",
" results_header = ['Attributes', 'WindowType', 'LogLoss']\n",
" \n",
" for params in ParameterGrid(grid):\n",
" score = build_model(\n",
" train, validation,\n",
" params['attributes'], 'Category',\n",
" time_window, params['window_type'],\n",
" algorithm['estimator'], None,\n",
" algorithm['predict_proba']\n",
" )\n",
" \n",
" results.append([params['attributes'], params['window_type'], '{:.5f}'.format(score)])\n",
" \n",
" pd.DataFrame(results, columns=results_header).to_csv(out_fname, index=False, mode='w')"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"FNAMES['discriminative-features'] = os.path.join(DATA_PATH, 'discriminative-features.csv')\n",
"\n",
"if not os.path.isfile(FNAMES['discriminative-features']):\n",
" discriminative_feature(\n",
" 336, ALGORITHMS['logit'], FNAMES['discriminative-features']\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"discriminative_features = pd.read_csv(FNAMES['discriminative-features'])\n",
"\n",
"discriminative_features['LogLoss'] = pd.to_numeric(discriminative_features['LogLoss'])"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WindowType cumulative exact\n",
"Attributes \n",
"all 2.59470 2.60675\n",
"cat 2.65614 2.65614\n",
"num 2.60562 2.61889\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAETCAYAAACWQRnGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZAV5b3/8feHASQKdxQC+SlQGYyEiCwDArJc0MiNEhJX3AhuaDDBuJtKMKZKfpbx5pd7K0a9ikEjGiURwaXEEI1JQGRRAR0XGBdUkMnViCCIRmT7/v44PePheIY5MzBzWubzqjrlOU8/3c+3Z4r5+HT36VZEYGZmlgYtil2AmZlZNYeSmZmlhkPJzMxSw6FkZmap4VAyM7PUcCiZmVlqtCx2AV90X/7yl6OsrKzYZZiZfaEsW7bs/YjomNvuUNpNZWVlLF26tNhlmJl9oUhana/dh+/MzCw1HEpmZpYaDiUzM0sNn1Mys2Zn69atVFVVsXnz5mKXstdr06YNXbp0oVWrVgX1dyiZWbNTVVVFu3btKCsrQ1Kxy9lrRQTr1q2jqqqKbt26FbSOD9+ZWbOzefNmOnTo4EBqZJLo0KFDvWakDiUza5YcSE2jvj9nh5KZWYEuv/xyfvOb39R8PvbYY/n+979f8/nKK6/k+uuv55RTTqnXdu+66y4uuuiiBtc1bdo0ysvLKS8vp3Xr1vTu3Zvy8nImTZrU4G0Wi88pmTWRskl/atLxVrX5XpOOx+SNTTteEQwbNoz777+fyy67jB07dvD+++/z4Ycf1ixftGgRN9xwAz/72c+atK7x48czfvx4IPOF/rlz5/LlL3+5SWvYUzxTMjMr0NChQ1m8eDEAy5cvp1evXrRr144PPviATz/9lMrKStq3b0+vXr2AzAzo5JNPZtSoUXTv3p2f/OQnNduaNm0aX//61xk0aBALFy6saV+1ahVHH300ffr0YeTIkbz99tts376dbt26ERFs2LCBkpIS5s+fD8CIESN4/fXX89a7Y8cODjnkENavXw/A9u3bOfjgg1m/fj1nnnkmEydO5PDDD+frX/86f/7znwHYtm0bV1xxBYMGDaJPnz7ccccde/4HuQsOJTOzAh100EG0bNmSt99+m0WLFjFkyBCOOOIIFi9ezNKlS+nduzetW7feaZ2KigpmzJjBSy+9xIwZM1izZg3vvPMO11xzDQsXLmTBggWsWLGipv/FF1/MOeecw4svvsi4ceO45JJLKCkpoUePHqxYsYIFCxbQv39/nnrqKT799FPWrFlD9+7d89bbokULxo4dyx/+8AcAHn/8cQYOHEj79u0BWLNmDUuWLGH27NlccMEFfPrpp0ydOpVOnTrx7LPPsmTJEm655RbefvvtRvqJ5qm5yUYyM9sLDB06lEWLFtWE0pAhQ2o+Dxs27HP9R44cSWlpKW3atKFnz56sXr2aZ555hqOOOoqOHTvSunVrTj/99Jr+ixcv5nvfyxx6Peuss1iwYAEAw4cPZ/78+cyfP5+rrrqKBQsWsGTJEgYOHLjLes8//3zuvvtuAO68886aw3wAp512Gi1atKBHjx507dqV119/nb/85S8156iOOOIINmzYUOtMrDE4lMzM6mHYsGEsWrSIl156iV69ejF48GAWL17MokWLGDp06Of677PPPjXvS0pK2LZtW4PGHTFiBE899RTPPvsso0ePZsOGDcybN4/hw4fvcr2ysjIOOOAA5s6dy/PPP88xxxxTsyz3yjhJRAS33norFRUVVFRU8NZbbzFy5MgG1dwQDiUzs3oYOnQojz76KO3bt6ekpIT27duzYcMGFi9enDeU8jniiCN48sknWbduHVu3bmXmzJk7bf++++4DYPr06TWhM2jQIBYtWkSLFi1o06YN5eXl/Pa3v2XEiBF1jnf++eczbtw4zjjjDFq0+OzP/syZM4kIXnvttZrDgMceeyy33nprTXi++uqrfPLJJwX/fHaXQ8nMrB569+7N+++/z+DBg3dqKy0tLfiKtwMPPJDJkyczZMgQhg0bxqGHHlqz7Oabb2batGn06dOHe+65hxtvvBHIzLi6du1aM+7w4cPZtGkTvXv3rnO8k046iY0bN3Luuefu1N65c2cGDBjAcccdx9SpU2ndujU/+MEP6N69O+Xl5fTq1YuJEyc2eHbXEIqIJhtsbzRgwIDw85SsEL4kPD0qKyt3CoK93dNPP81VV13F3Llza9rOPPNMTjnlFE488cRGHz/fz1vSsogYkNvX31MyM9uL/eIXv2Dq1Kk1hwTTzqFkZrYXu/rqq7n66qs/137vvfcWoZq6+ZySmZmlhkPJzMxSw6FkZmap4VAyM7PUcCiZme3lVq1aVXOT2F31qb5HHsDSpUu55JJLGru0z/HVd2bW7O3p75Ct+uV39uj2mkJ1KFXfd2/AgAEMGPC5rxE1Os+U8pC0n6S7Jd0uaVyx6zGzvdPvf/97+vTpQ9++fTnrrLM499xzmTVrVs3ytm3bAjBv3jyOPPJITjjhBA4++GAmTZrE9OnTGTRoEL179+aNN94AqHX9bKtWrWL48OH079+f/v37s2jRIgAmTZrEU089RXl5OTfccAPz5s3ju9/9Ljt27KCsrIwNGzbUbKN79+7885//ZO3atYwZM4aBAwcycODAnR7B0VCNFkqSukqaK2mFpOWSLq2l3/6SZkl6RVKlpCFJ+ypJL0mqkLRbt0yQdKek9yS9nNM+StKrklZKyn5E48nArIiYABy/O2ObmeWzfPlyrrvuOv7+97/zwgsv1NxOqDYvvPACt912G5WVldxzzz289tprPPvss3z/+9/n5ptvLnjcTp068cQTT/Dcc88xY8aMmkN0v/zlLxk+fDgVFRVcfvnlNf1btGjBCSecwEMPPQTAM888w1e/+lW+8pWvcOmll3L55ZezZMkSHnjggZ2ewttQjXn4bhtwZUQ8J6kdsEzSExGxIqffjcBjEXGKpNbAvlnLvhkR79c2gKROwCcRsSmr7ZCIWJnT9S7gf4DfZ/UrAW4BvgVUAUskPZLU1wV4Kem6vfBdNjMrzN///ndOPfXUmvvlVT/jqDYDBw7kwAMPBOBrX/tazd2+e/fuvdPtg+qydetWLrroIioqKigpKeG1116rc53TTz+da6+9lvHjx3PffffVPGrjr3/9607Pgvrwww/56KOP8s7QCtVooRQR7wDvJO83SaoEOgM1eyCpFBgBnJv02wJsqccwRwI/lDQ6Ij6VNIHMLOfbObXMl1SWs+4gYGVEvJnUch9wQlJfFZlgqqCW2aSk44DjDjnkkHqUa2ZWu5YtW7Jjxw4g89TYLVs++3OY/QiMFi1a1Hxu0aJFzQ1Td7V+tRtuuIGvfOUrvPDCC+zYsYM2bdrUWdeQIUNYuXIla9eu5eGHH+bnP/95zRhPP/10QdsoVJOcU0oCoR/wTM6ibsBaYJqk5yXdIWm/ZFkAf5G0TNIF+bYbETOBx4EZybmf84BTCyyrM7Am63NV0gbwIDBG0hRgdi1jz46IC0pLSwsczszsM0cffTQzZ85k3bp1AKxfv56ysjKWLVsGwCOPPMLWrVvrtc1C1t+4cSMHHnggLVq04J577mH79szBoHbt2rFp06bP9YfMc5ZOOukkrrjiCg499FA6dOgAwDHHHLPTocOKiop61ZtPo4eSpLbAA8BlEfFhzuKWQH9gSkT0Az4Gqs/t/HtE9Ccz6/mRpLwPDYmIXwGbgSnA8RHx0e7WHBEfR8T4iJgYEdN3d3tmZrkOO+wwrr76ao488kj69u3LFVdcwYQJE3jyySfp27cvixcvZr/99qt7Q1kKWf/CCy/k7rvvpm/fvrzyyis1ffr06UNJSQl9+/blhhtu+Nx6p59+Ovfee+9OT8m96aabWLp0KX369KFnz57cdttt9fwpfF6jPrpCUivgUeDxiPh1nuX/B3g6IsqSz8OBSRHxnZx+k4GPIuK/82xjOJlAWgZsioiLaqmlDHg0Inoln4cAkyPi2OTzVQAR8Z/12Uc/usIK5UdXpEdze3RFsdXn0RWNefWdgN8BlfkCCSAi3gXWSOqRNI0EViSXZLdLtrMfcAzwcu76kvoBU8mcCxoPdJB0XYElLgG6S+qWXGBxBvBIwTtoZmZ7XGMevhsGnAUcnVzWXSFpNICkOZIOSvpdDEyX9CJQDlwPfAVYIOkF4FngTxHxWJ4x9gVOi4g3ImIHcDawOreTpD8Ci4EekqoknR8R24CLyJyTqgTuj4jle273zcysvhrz6rsFgGpZNjrrfQWQO4X7AOhbwBgLcz5vBW7P029sLevPAebUNY6ZmTUN39HBzMxSw6FkZmap4VAyM7PUcCiZme1l5s2bV3Oj1S8aP7rCzGzyHr4zS5G/szVv3jzatm3L0KFDi1pHQ3imZGZWJPfeey+DBg2ivLycH/zgB6xevZru3bvz/vvvs2PHDoYPH85f/vIXAE488UQOP/xwDjvsMKZOnVqzjccee4z+/fvTt29fRo4cyapVq7jtttu44YYbKC8v56mnnirW7jWIZ0pmZkVQWVnJjBkzWLhwIa1ateLCCy/kySef5Kc//SkTJ05k0KBB9OzZs+Zu4HfeeSft27fnk08+YeDAgYwZM4YdO3YwYcIE5s+fT7du3Vi/fj3t27fnhz/8IW3btuXHP/5xkfey/hxKZmZF8Le//Y1ly5YxcOBAAD755BM6derE5MmTmTlzJrfddttONzi96aabap5ptGbNGl5//XXWrl3LiBEj6NatG1D34y++CBxKZmZFEBGcc845/Od/7ny7zX/9619UVVUB8NFHH9GuXTvmzZvHX//6VxYvXsy+++7LUUcdxebNm4tRdqPzOSUzsyIYOXIks2bN4r333gMyj65YvXo1P/3pTxk3bhzXXnstEyZMADKPmzjggAPYd999eeWVV3j66acBGDx4MPPnz+ett96q2Qbs+jEUaedQMjMrgp49e3LddddxzDHH0KdPH771rW+xatUqlixZUhNMrVu3Ztq0aYwaNYpt27Zx6KGHMmnSJAYPHgxAx44dmTp1KieffDJ9+/ateazEcccdx0MPPfSFvNChUR9d0Rz40RVWKD+6Ij386IqmlYpHV5iZmdWXQ8nMzFLDoWRmZqnhUDKzZsnn05tGfX/ODiUza3batGnDunXrHEyNLCJYt24dbdq0KXgdf3nWzJqdLl26UFVVxdq1a4tdyl6vTZs2dOnSpeD+DiUza3ZatWpVc2seSxcfvjMzs9RwKJmZWWr48F0ekvYDbgW2APMiYnqRSzIzaxYabaYkqaukuZJWSFou6dJa+u0vaZakVyRVShqStaxE0vOSHt3NWu6U9J6kl3PaR0l6VdJKSZOyFp0MzIqICcDxuzO2mZkVrjEP320DroyInsBg4EeSeubpdyPwWER8A+gLVGYtuzTn804kdZLULqftkDxd7wJG5fQrAW4Bvg30BMZm1dcFWJO8317b+GZmtmc1WihFxDsR8VzyfhOZcOmc3UdSKTAC+F3Sb0tEbEiWdQG+A9yxi2GOBB6WtE+yzgTg5jy1zAfW5zQPAlZGxJsRsQW4DzghWVZFJpjA593MzJpMk/zBlVQG9AOeyVnUDVgLTEsO092RnM8B+A3wE2BHbduNiJnA48AMSeOA84BTCyyrM5/NhiATRNWh+SAwRtIUYHYt+3ScpKkbN35x7oxsZpZ2jR5KktoCDwCXRcSHOYtbAv2BKRHRD/gYmCTpu8B7EbGsru1HxK+AzcAU4PiI+Gh3a46IjyNifERMrO0ih4iYHREXlJaW7u5wZmaWaNRQktSKTCBNj4gH83SpAqoionoGNYtMSA0Djpe0isxhtaMl3VvLGMOBXsBDwDX1KO8fQNesz12SNjMzK5LGvPpOZM4VVUbEr/P1iYh3gTWSeiRNI4EVEXFVRHSJiDLgDODvEXFmnjH6AVPJnAsaD3SQdF2BJS4BukvqJql1Ms4jhe+hmZntaY05UxoGnEVmllORvEYDSJoj6aCk38XAdEkvAuXA9fUYY1/gtIh4IyJ2AGcDq3M7SfojsBjoIalK0vkRsQ24iMw5qUrg/ohY3rBdNTOzPaHRvjwbEQsA1bJsdNb7CuBzj8TNWj4PmFfLsoU5n7cCt+fpN7aW9ecAc2ob28zMmpYvdzYzs9RwKJmZWWo4lMzMLDUcSmZmlhoOJTMzSw2HkpmZpYZDyczMUsOhZGZmqeFQMjOz1HAomZlZajiUzMwsNRxKZmaWGg4lMzNLDYeSmZmlhkPJzMxSw6FkZmap4VAyM7PUcCiZmVlqOJTMzCw1HEpmZpYaDiUzM0uNOkNJ0tck7ZO8P0rSJZL2b/zSzMysuSlkpvQAsF3SIcBUoCvwh0atyszMmqWWBfTZERHbJJ0E3BwRN0t6vrELKyZJ+wG3AluAeRExvcglmZk1C4XMlLZKGgucAzyatLWqayVJXSXNlbRC0nJJl9bSb39JsyS9IqlS0hBJbSQ9K+mFZN3/W/gu5R3jTknvSXo5p32UpFclrZQ0KWvRycCsiJgAHL87Y5uZWeEKCaXxwBDgFxHxlqRuwD0FrLcNuDIiegKDgR9J6pmn343AYxHxDaAvUAl8ChwdEX2BcmCUpMG5K0rqJKldTtsheca4CxiV068EuAX4NtATGJtVXxdgTfJ+ewH7amZme0CdoRQRKyLikoj4o6QDgHYR8f8KWO+diHgueb+JTNh0zu4jqRQYAfwu6bclIjZExkdJt1bJK/IMcyTwcNaFGBOAm/PUMh9Yn9M8CFgZEW9GxBbgPuCEZFkVmWCCWn5Gko6TNHXjxo21/QjMzKyeCrn6bp6kf5PUHngOuF3Sr+sziKQyoB/wTM6ibsBaYJqk5yXdkZzPQVKJpArgPeCJiMhdl4iYCTwOzJA0DjgPOLXAsjrz2WwIMkFUHZoPAmMkTQFm51s5ImZHxAWlpaUFDmdmZnUp5PBdaUR8SOY8y+8j4gjgPwodQFJbMlfwXZZsJ1tLoD8wJSL6AR8DkwAiYntElJOZsQyS1Cvf9iPiV8BmYApwfNYMq8Ei4uOIGB8RE32Rg5lZ0ykklFpKOhA4jc8udCiIpFZkAml6RDyYp0sVUJU1C5pFJqRqRMQGYC4554SyxhgO9AIeAq6pR3n/IHN5e7UuSZuZmRVJIaF0LZlDZG9ExBJJBwOv17WSJJE5V1QZEXkP90XEu8AaST2SppHACkkdq7+gK+lLwLeAV/KM0Y/Md6dOIHNBRgdJ1xWwTwBLgO6SuklqDZwBPFLgumZm1ggKudBhZkT0iYiJyec3I2JMAdseBpwFHC2pInmNBpA0R9JBSb+LgemSXiRzpd31wIHA3KRtCZlzSvlmafsCp0XEGxGxAzgbWJ3bSdIfgcVAD0lVks6PiG3ARWQCtxK4PyKWF7BfZmbWSOr88qykLmSuaBuWND0FXBoRVbtaLyIWAKpl2eis9xXAgJwuH5C5MGKXImJhzuetwO15+o2tZf05wJy6xjEzs6ZRyOG7aWQOax2UvGYnbWZmZntUIaHUMSKmRcS25HUX0LGR6zIzs2aokFBaJ+nM5HtDJZLOBNY1dmFmZtb8FBJK55G5HPxd4B3gFODcRqzJzMyaqUKuvlsdEcdHRMeI6BQRJwKFXH1nZmZWLw198uwVe7QKMzMzGh5KeS/1NjMz2x0NDaV8d+w2MzPbLbV+eVbSJvKHj4AvNVpFZmbWbNUaShHRrrZlZmZmjaGhh+/MzMz2OIeSmZmlhkPJzMxSw6FkZmapUcijK/JdhbcRWApcGRFvNkZhZmbW/NQZSsBvyDy2/A9kLgc/A/ga8BxwJ3BUYxVnZmbNSyGH746PiN9GxKaI+DAipgLHRsQM4IBGrs/MzJqRQkLpX5JOk9QieZ0GbE6W+c4OZma2xxQSSuOAs4D3ktdZwJmSvgRc1Ii1mZlZM1PnOaXkQobjalm8YM+WY2ZmzVmdMyVJXSQ9JOm95PWApC5NUZyZmTUvhRy+mwY8AhyUvGYnbWZmZntUIaHUMSKmRcS25HUX0LGR6yoqSftJulvS7ZLGFbseM7PmopBQWifpTEklyetMYF1dK0nqKmmupBWSlku6tJZ++0uaJekVSZWShhS6bqEk3Zkcenw5p32UpFclrZQ0KWvRycCsiJgAHL87Y5uZWeEKCaXzgNOAd4F3gFOAcwtYbxuZOz70BAYDP5LUM0+/G4HHIuIbQF+gstB1JXWS1C6n7ZA8Y9wFjMrpVwLcAnwb6AmMzRqjC7Ameb+9gH01M7M9oM5QiojVEXF8RHSMiE4RcSIwpoD13omI55L3m8iETefsPpJKgRHA75J+WyJiQyHrJo4EHpa0T7K9CcDNeWqZD6zPaR4ErIyINyNiC3AfcEKyrIpMMIHvD2hm1mQa+gf3ivp0llQG9AOeyVnUDVgLTJP0vKQ7JO1X4LpExEzgcWBGcu7nPODUAsvqzGezIcgEUXXwPQiMkTSFzIUd+fbpOElTN27cWOBwZmZWl4aGkgruKLUFHgAui4gPcxa3BPoDUyKiH/AxMKnAdQGIiF+RucPEFDK3RPqoPjtSyzY/jojxETExIqbX0md2RFxQWlq6u8OZmVmioaFU0O2FJLUiEyrTI+LBPF2qgKqIqJ4FzSITUoWsWz3GcKAX8BBwTcF7AP8AumZ97pK0mZlZkdR6R4daHlkBmVnSl+rasCSROVdUGRG/ztcnIt6VtEZSj4h4FRgJrChk3WSMfsBU4LvAW8B0SddFxM/rqg9YAnSX1I1MGJ0BfK+A9czMGt/kJj4KMzkdpyJqnSlFRLuI+Lc8r3YRUcgjL4aRuU/e0ZIqktdoAElzJB2U9LuYTJi8CJQD1+9q3Rz7AqdFxBsRsQM4G1id20nSH4HFQA9JVZLOj4htZO7d9ziZCynuj4jlBeyXmZk1kkLCpUEiYgG1nHuKiNFZ7yuAATldal03ZzsLcz5vBW7P029sLevPAebUNY6ZmTUNX+5sZmap4VAyM7PUcCiZmVlqOJTMzCw1HEpmZpYaDiUzM0sNh5KZmaWGQ8nMzFLDoWRmZqnhUDIzs9RwKJmZWWo4lMzMLDUcSmZmlhoOJTMzSw2HkpmZpYZDyczMUsOhZGZmqeFQMjOz1HAomZlZajiUzMwsNRxKZmaWGg4lMzNLjZbFLiCNJO0H3ApsAeZFxPQil2Rm1iw02kxJUldJcyWtkLRc0qW19Ntf0ixJr0iqlDQkab9T0nuSXt4DteTdlqRRkl6VtFLSpKxFJwOzImICcPzujm9mZoVpzMN324ArI6InMBj4kaSeefrdCDwWEd8A+gKVSftdwKhdDSCpk6R2OW2H5On6uW1JKgFuAb4N9ATGZtXXBViTvN++qxrMzGzPabRQioh3IuK55P0mMmHTObuPpFJgBPC7pN+WiNiQvJ8PrK9jmCOBhyXtk2xvAnBznlrybWsQsDIi3oyILcB9wAnJsioywQQ+72Zm1mSa5A+upDKgH/BMzqJuwFpgmqTnJd2RnM8pSETMBB4HZkgaB5wHnFrg6p35bDYEmSCqDs0HgTGSpgCz860s6ThJUzdu3FhouWZmVodGDyVJbYEHgMsi4sOcxS2B/sCUiOgHfAxMoh4i4lfAZmAKcHxEfLS7NUfExxExPiIm1naRQ0TMjogLSktLd3c4MzNLNGooSWpFJpCmR8SDebpUAVURUT2DmkUmpOozxnCgF/AQcE09Vv0H0DXrc5ekzczMiqQxr74TmXNFlRHx63x9IuJdYI2kHknTSGBFPcboB0wlcy5oPNBB0nUFrr4E6C6pm6TWwBnAI4WObWZme15jzpSGAWcBR0uqSF6jASTNkXRQ0u9iYLqkF4Fy4Pqkzx+BxUAPSVWSzs8zxr7AaRHxRkTsAM4GVud2yretiNgGXETmnFQlcH9ELN9zu29mZvXVaF+ejYgFgGpZNjrrfQUwIE+fsQWMsTDn81bg9kK3FRFzgDl1jWNmZk3DlzubmVlqOJTMzCw1HEpmZpYaDiUzM0sN3yXcdm1yE385eLLvkGHWnHmmZGZmqeFQMjOz1HAomZlZavickplZAcom/alJx1vVpkmHSw3PlMzMLDUcSmZmlhoOJTMzSw2HkpmZpYZDyczMUsOhZGZmqeFQMjOz1HAomZlZajiUzMwsNRxKZmaWGg4lMzNLDYeSmZmlhkPJzMxSw3cJr4Wk/YBbgS3AvIiYXuSSzMz2ekWZKUnqKmmupBWSlku6tJZ++0uaJekVSZWShuzGmHdKek/SyzntoyS9KmmlpElZi04GZkXEBOD4ho5rZmaFK9bhu23AlRHRExgM/EhSzzz9bgQei4hvAH2ByuyFkjpJapfTdkgtY94FjMrpWwLcAnwb6AmMzaqjC7Ameb+9wP0yM7PdUJRQioh3IuK55P0mMmHTObuPpFJgBPC7pN+WiNiQs6kjgYcl7ZOsMwG4uZYx5wPrc5oHASsj4s2I2ALcB5yQLKsiE0zgc29mZk2i6H9sJZUB/YBnchZ1A9YC0yQ9L+mO5DxPjYiYCTwOzJA0DjgPOLUew3fms9kQZIKoOhwfBMZImgLMzlP3cZKmbty4sR7DmZnZrhQ1lCS1BR4ALouID3MWtwT6A1Mioh/wMTAppw8R8StgMzAFOD4iPtoTtUXExxExPiIm5rvIISJmR8QFpaWle2I4MzOjiKEkqRWZQJoeEQ/m6VIFVEVE9QxqFpmQyt3OcKAX8BBwTT3L+AfQNetzl6TNzMyKoFhX34nMuaLKiPh1vj4R8S6wRlKPpGkksCJnO/2AqWTOA40HOki6rh6lLAG6S+omqTVwBvBIvXbGzMz2mGLNlIYBZwFHS6pIXqMBJM2RdFDS72JguqQXgXLg+pzt7AucFhFvRMQO4Gxgdb4BJf0RWAz0kFQl6fyI2AZcROa8VCVwf0Qs37O7amZmhSrKl2cjYgGgWpaNznpfAQzYxXYW5nzeCtxeS9+xtbTPAebUXbWZmTW2ol99Z2ZmVs2hZGZmqeFQMjOz1HAomZlZajiUzMwsNRxKZmaWGg4lMzNLDYeSmZmlhkPJzMxSw6FkZmap4VAyM7PUcCiZmVlqFOWGrNZwZZP+1KTjrWrTpMOZWTPnmZKZmaWGQ8nMzFLDoWRmZkhXEn4AAAbbSURBVKnhUDIzs9RwKJmZWWo4lMzMLDUcSmZmlhoOJTMzSw1FRLFr+EKTtBZYXew6GtGXgfeLXYQ1iH93X2x7++/vqxHRMbfRoWS7JGlpRAwodh1Wf/7dfbE119+fD9+ZmVlqOJTMzCw1HEpWl6nFLsAazL+7L7Zm+fvzOSUzM0sNz5TMzCw1HEpmZpYaDiUzM0sNh5LZXkTSsELazNLKFzqY7UUkPRcR/etqs3SStD9wNlAGtKxuj4hLilVTU2tZdxdrDiRtAvL9H4qAiIh/a+KSrB4kDQGGAh0lXZG16N+AkuJUZQ0wB3gaeAnYUeRaisKhZABERLti12C7pTXQlsy/6ezf5YfAKUWpyBqiTURcUXe3vZcP3xkAktrvanlErG+qWqzhJH01IvbmGwTv1SRdDnwEPAp8Wt3enP79eaZk1ZaROXynrLbqzwEcXIyirN7+Jem/gMOANtWNEXF08UqyetgC/BdwNZ8dTm9W//4cSgZARHSrfp/MmrqT9UfNvjCmAzOA7wI/BM4B1ha1IquPK4FDImJvfmTFLjmUbCeSvg9cCnQBKoDBwCJgZDHrsoJ1iIjfSbo0Ip4EnpS0pNhFWcFWAv8qdhHF5FCyXJcCA4GnI+Kbkr4BXF/kmqxwW5P/viPpO8D/Ars8X2ip8jFQIWkuO59T8iXh1mxtjojNkpC0T0S8IqlHsYuygl0nqZTMYaCbyVwSfllxS7J6eDh5NVsOJctVlXyB72HgCUkfsHc/7n1vcyqwICJeBr6ZnB/8b2B2ccuyQkTE3cWuodh8SbjVStKRQCnwWERsKXY9VjdJz0dEv7raLJ0kvUWeL7FHhK++M0tOlNsXSwtJB0TEB1BzJaX/nX9xDMh634bMzLdZnRP0TMlsLyLpbOBnwMyk6VTgFxFxT/Gqst0haVlEHF7sOpqK/w/KbC8SEb+XtBSo/rLsyRGxopg1WeEkZd84twWZmVOz+jvtmZKZWUokl4JX/1HeBqwC/jsiXitaUU3MoWRmlhKS2gBj2PnRFRER1xatqCbWrKaFZmYp9zCwAXgO2FzkWorCMyUzs5SQ9HJE9Cp2HcXkx6GbmaXHIkm9i11EMXmmZGaWEpJWAIcAb5G59131k5/7FLWwJuRQMjNLCUlfzdfenB7c6FAyM7PU8DklMzNLDYeSmZmlhkPJbDdIulrSckkvSqqQdEQd/c+VdFA9x/hGsu3nJX0tz/JySSFpVFbb/pIuzPpcJul7uxjjIEmzsmr8n3rWWO/9MsvHoWTWQJKGAN8F+idXR/0HsKaO1c4F6vvH+0RgVkT0i4g38iwfCyxI/lttf+DCrM9lQN5QktQyIv43Ik6pZ13ZzqX++2X2Ob6jg1nDHQi8HxGfAkTE+9ULJB0O/BpoC7xP5o/2MDI32Jwu6RNgSER8krVOOXAbsC/wBnAeMITMk2O3SxoZEd/MLkCSyNwJ/FvAU5LaRMRm4JfA1yRVAE8Aw4FDk893Ax8AJyf1lUg6B3g064ubXSXNAzoD90bE/5VUlt1H0o+T9V/O3S+gZ+7+R8Q7ki4Bfkjmvm4rIuKMev/Ube8WEX755VcDXmT+4FYArwG3Akcm7a2ARUDH5PPpwJ3J+3nAgFq292LWNq4FfpO8nwz8uJZ1hgF/S97/ARiTvC8DXs7qdxSZQKn+fC5QBbTP7Z8sewfoAHyJz0Ind5s/Bibn7lcd+/+/wD7J+/2L/Tv0K30vz5TMGigiPkpmRMOBbwIzJE0ClgK9yDxOHqCEzB/5WkkqJfNHuvrBinfz2TORdmUscF/y/j7gbOCBAnfhiYhYv4tl65LaHgT+ncx92QrRg9r3/0UyM6qH67E9a0YcSma7ISK2k5klzJP0EnAOsAxYHhFDGnNsSSVk7ih9gqSryXz7v4OkdgVu4uNdLMv9AmOQOeSWfR66TW2lUfv+fwcYARwHXC2pd0RsK7BeawZ8oYNZA0nqIal7VlM5sBp4FeiYXAiBpFaSDkv6bAI+FxoRsRH4QNLwpOksoK7H0Y8EXoyIrhFRFhFfJTNLOinPOHnH3YVvSWov6UtkLrRYCPwT6CSpg6R9yFzkkW/7efdfUguga0TMBX4KlJI5BGpWwzMls4ZrC9wsaX8ys4iVwAURsUXSKcBNyWG5lsBvgOXAXcBt+S50IDPLuk3SvsCbwPg6xh8LPJTT9gAwMTJPoF0o6WXgz2Qekb5d0gtJDR/Use1nk211IXOhw1IASdcmy/4BvJLVf6f9AvLt/2vAvUmbgJsiYkMddVgz49sMmZlZavjwnZmZpYZDyczMUsOhZGZmqeFQMjOz1HAomZlZajiUzMwsNRxKZmaWGg4lMzNLjf8PPcZiwTR1l9QAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"chart_data = discriminative_features.pivot(index='Attributes', columns='WindowType', values='LogLoss')\n",
"\n",
"print(chart_data)\n",
"\n",
"chart_data.plot(kind='bar', grid=False, legend=True)\n",
"\n",
"plt.xlabel('Set of Attributes')\n",
"plt.ylabel('Log Loss')\n",
"\n",
"plt.yscale(\"log\")\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These results show that the union of the numerical and categorical attributes contribute the most to discriminative power. In a like manner, aggregated attributes computed using cumulative time windows rather than an exact time window are more discriminative."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3.4 Stacking Ensemble"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The methodology to develop a stacked generalization ensemble goes from splitting the initial dataset to producing the second-level predictions, as depicted in Figure 2."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![Stacking Ensemble](images/ensemble.svg \"Stacking Ensemble\")\n",
"<div align=\"center\" style=\"margin-top: 10px;\"><b>Figure 2</b>: Methodology to develop a stacked generalization ensemble</div>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Firstly, the initial dataset is split into $k$ folds through a stratified sampling method. In particular, the resulting folds are different from each other, i.e., a crime instance belongs to only one fold. Secondly, $k$ base models are trained, one per each fold. Thirdly, the first-level predictions, i.e., the outputs from the base models, are combined depending on the stacking technique the meta-model implements. Finally, the meta-model outputs the second-level predictions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Furthermore, there are other considerations regarding the methodology described above, namely:\n",
"\n",
"1. The way the first-level predictions are stacked depends on the stacking technique. If the stacking technique is model averaging, then such predictions are summed element-wise. Otherwise, matrices representing first-level predictions are stacked horizontally.\n",
"2. Accordingly, there are two stacking techniques a meta-model might implement, namely: training a classifier or model averaging.\n",
"3. Potentially, combining the predictions from several meta-models might be used to produce third-level predictions."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"def write_in_file(fname, content, mode='w', insert_new_line=True):\n",
" with open(fname, mode) as f:\n",
" content = (content\n",
" + ('\\n' if insert_new_line else ''))\n",
" \n",
" f.write(content)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"@numba.jit(nopython=True)\n",
"def is_leap_year(year):\n",
" return (True\n",
" if (year % 4 == 0 and (year % 100 != 0 or year % 400 == 0))\n",
" else False)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"@numba.jit(nopython=True)\n",
"def get_number_of_days_in_month(data):\n",
" n = len(data)\n",
" \n",
" days_in_month = [31, None, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]\n",
" \n",
" result = []\n",
" for i in range(n):\n",
" year = data[i,0]\n",
" month = data[i,1]\n",
" \n",
" result.append(\n",
" (29 if is_leap_year(year) else 28)\n",
" if month == 2\n",
" else days_in_month[month-1]\n",
" )\n",
" \n",
" return result"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"def encode_cyclical_attributes(\n",
" df,\n",
" attributes=['Month', 'Day', 'DayOfWeek', 'Hour']):\n",
" \"\"\"Encoding of cyclical attributes.\n",
" \n",
" Source: <http://blog.davidkaleko.com/feature-engineering-cyclical-features.html>\n",
" \"\"\"\n",
" df['days_month'] = get_number_of_days_in_month(df[['Year', 'Month']].to_numpy().astype(int))\n",
" \n",
" for attr in attributes:\n",
" max_val = df[attr].max() if attr != 'Day' else df['days_month']\n",
" min_val = df[attr].min()\n",
" \n",
" if min_val == 1:\n",
" df[attr] -= 1\n",
" elif min_val == 0:\n",
" max_val += 1\n",
" \n",
" df['{}_Sin'.format(attr)] = np.sin(df[attr] * 2 * np.pi / max_val)\n",
" df['{}_Cos'.format(attr)] = np.cos(df[attr] * 2 * np.pi / max_val)\n",
" \n",
" attributes.append('days_month')\n",
" df = df.drop(columns=attributes)\n",
" \n",
" return df"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"def encode_categorical_attributes(df, attributes):\n",
" \"\"\"One-Hot encoding of categorical attributes.\"\"\"\n",
" return pd.get_dummies(df, prefix=attributes, columns=attributes)"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"def split_training_data(df, n_splits, split_criterion='Category'):\n",
" \"\"\"Split training data in a stratified way without replacement.\n",
" \n",
" Please note that if the number of class instances is less than the\n",
" number of splits, a stratified random sample with replacement will\n",
" be performed.\n",
" \"\"\"\n",
" criterion_values = {\n",
" value: df.loc[df[split_criterion]==value].shape[0]\n",
" for value in df[split_criterion].unique()\n",
" }\n",
" \n",
" sampled_instances = []\n",
" \n",
" for i in range(n_splits):\n",
" last_split = True if i == (n_splits-1) else False\n",
" \n",
" sample = None\n",
" for value, size in criterion_values.items():\n",
" mask = (df[split_criterion] == value) & (~df['Id'].isin(sampled_instances))\n",
" criterion_sample = df.loc[mask]\n",
" \n",
" sample_size = size * (1/n_splits)\n",
" sample_size = int(np.round(sample_size, 0))\n",
" \n",
" if sample_size < n_splits:\n",
" criterion_sample = df.loc[(df[split_criterion] == value)]\n",
" sample_size = np.min([n_splits, criterion_sample.shape[0]])\n",
" \n",
" if not last_split:\n",
" criterion_sample = criterion_sample.sample(\n",
" n=sample_size, replace=False, random_state=RANDOM_STATE\n",
" )\n",
" \n",
" sample = (\n",
" sample.append(criterion_sample).reset_index(drop=True)\n",
" if sample is not None\n",
" else criterion_sample.copy(deep=True)\n",
" )\n",
" \n",
" sampled_instances += criterion_sample['Id'].values.tolist()\n",
" \n",
" yield sample"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"def build_meta_classifier(\n",
" prediction_fnames,\n",
" y_true,\n",
" stack_method,\n",
" estimator):\n",
" \"\"\"Build a meta classifier on the basis of several learners and a stacking method.\"\"\"\n",
" X = None\n",
" \n",
" for pred_fname in prediction_fnames:\n",
" clf_pred = np.loadtxt(pred_fname, dtype=float, delimiter=',', skiprows=1)\n",
" clf_pred = clf_pred[:,1:]\n",
" \n",
" if X is None:\n",
" X = copy.deepcopy(clf_pred)\n",
" continue\n",
" \n",
" X = np.add(X, clf_pred) if stack_method == 'soft_voting' else np.hstack([X, clf_pred])\n",
" \n",
" if stack_method == 'soft_voting':\n",
" y_pred = X / len(prediction_fnames)\n",
" else:\n",
" cv = StratifiedKFold(n_splits=5, random_state=RANDOM_STATE)\n",
" \n",
" y_pred = cross_val_predict(\n",
" estimator=estimator, X=X, y=y_true, cv=cv, n_jobs=5, method='predict_proba')\n",
" \n",
" return log_loss(y_true, y_pred)"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
"def build_classifier_ensembles(\n",
" out_fname,\n",
" base_algorithm,\n",
" meta_algorithm,\n",
" crime_category=CRIME_CATEGORY,\n",
" time_windows=TIME_WINDOWS):\n",
" \"\"\"Build classifier ensembles.\"\"\"\n",
" if not os.path.isfile(out_fname):\n",
" content = [\n",
" 'Id',\n",
" 'TimeWindow',\n",
" 'NumberOfEstimators',\n",
" 'StackingMethod',\n",
" 'EncodingOfCyclicalAttributes',\n",
" 'LogLoss'\n",
" ] \n",
" write_in_file(out_fname, ','.join(content))\n",
" \n",
" prediction_header = ['' for j in range(len(crime_category))]\n",
" for crime, idx in crime_category.items():\n",
" prediction_header[idx] = crime\n",
" \n",
" ensemble_grid = {\n",
" 'encode_cyclical_attr': ['cyclical', 'onehot'],\n",
" 'n_estimators': [3, 5, 7],\n",
" 'stack_method': ['soft_voting', 'clf'],\n",
" 'time_window': time_windows\n",
" }\n",
" ensemble_grid = ParameterGrid(ensemble_grid)\n",
" \n",
" for setting_id, settings in enumerate(ensemble_grid):\n",
" time_window = settings['time_window']\n",
" \n",
" window_path = {\n",
" '/': os.path.join(DATA_PATH, '{}H'.format(time_window))\n",
" } \n",
" window_path['/data'] = os.path.join(window_path['/'], 'data')\n",
" \n",
" train = pd.read_csv(os.path.join(window_path['/data'], 'train.csv'))\n",
" validation = pd.read_csv(os.path.join(window_path['/data'], 'validation.csv'))\n",
" \n",
" if settings['encode_cyclical_attr'] == 'cyclical':\n",
" train = encode_cyclical_attributes(train)\n",
" validation = encode_cyclical_attributes(validation)\n",
" else:\n",
" cyclical_attr = ['Month', 'Day', 'DayOfWeek', 'Hour']\n",
" train = encode_categorical_attributes(train, cyclical_attr)\n",
" validation = encode_categorical_attributes(validation, cyclical_attr)\n",
" \n",
" date_attr = [\n",
" 'Quarter', 'Triannual', 'Semester', 'Fortnight',\n",
" 'FourHourPeriod', 'SixHourPeriod', 'TwelveHourPeriod'\n",
" ] \n",
" train = encode_categorical_attributes(train, date_attr)\n",
" validation = encode_categorical_attributes(validation, date_attr)\n",
" \n",
" train = encode_categorical_attributes(train, ['PdDistrict'])\n",
" validation = encode_categorical_attributes(validation, ['PdDistrict'])\n",
" \n",
" train = train.replace({'Category': crime_category})\n",
" validation = validation.replace({'Category': crime_category})\n",
" \n",
" window_path['/prediction'] = os.path.join(\n",
" window_path['/'],\n",
" 'prediction',\n",
" 'validation',\n",
" 'n_estimators={}'.format(settings['n_estimators']),\n",
" 'cyclical_attr={}'.format(settings['encode_cyclical_attr'])\n",
" )\n",
" \n",
" if not os.path.isdir(window_path['/prediction']):\n",
" os.makedirs(window_path['/prediction'])\n",
" \n",
" learners_result_fname = os.path.join(window_path['/prediction'], 'learners-result.csv')\n",
" if not os.path.isfile(learners_result_fname):\n",
" write_in_file(learners_result_fname, ','.join(['Clf', 'LogLoss']))\n",
" \n",
" prediction_fnames = []\n",
" for i, train_sample in enumerate(split_training_data(train, settings['n_estimators'])):\n",
" pred_fname = os.path.join(window_path['/prediction'], 'clf-{}-pred.csv'.format(i))\n",
" prediction_fnames.append(pred_fname)\n",
" \n",
" if os.path.isfile(pred_fname):\n",
" continue\n",
" \n",
" score, predictions = build_model(\n",
" train_sample, validation,\n",
" 'all', 'Category',\n",
" time_window, 'cumulative',\n",
" base_algorithm['estimator'], None,\n",
" base_algorithm['predict_proba'],\n",
" return_pred=True\n",
" )\n",
" \n",
" write_in_file(\n",
" learners_result_fname, ','.join([str(i), '{:.5f}'.format(score)]), mode='a'\n",
" )\n",
" \n",
" predictions = pd.DataFrame(predictions, columns=prediction_header)\n",
" \n",
" predictions['Id'] = validation['Id'].values\n",
" predictions = predictions[['Id']+prediction_header]\n",
" \n",
" predictions.to_csv(pred_fname, float_format='%.5f', index=False)\n",
" \n",
" score = build_meta_classifier(\n",
" prediction_fnames,\n",
" validation['Category'].values.astype(int),\n",
" settings['stack_method'],\n",
" meta_algorithm['estimator']\n",
" )\n",
" \n",
" content = [\n",
" str(setting_id),\n",
" str(time_window),\n",
" str(settings['n_estimators']),\n",
" settings['stack_method'],\n",
" settings['encode_cyclical_attr'],\n",
" '{:.5f}'.format(score)\n",
" ]\n",
" write_in_file(out_fname, ','.join(content), mode='a')"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"FNAMES['classifier-ensembles'] = os.path.join(DATA_PATH, 'classifier-ensembles.csv')\n",
"\n",
"if not os.path.isfile(FNAMES['classifier-ensembles']):\n",
" build_classifier_ensembles(\n",
" FNAMES['classifier-ensembles'], ALGORITHMS['logit'], ALGORITHMS['logit']\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 3.4.1 Results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Before starting the analysis of results, some clarifications must be provided, namely:\n",
"\n",
"1. The Logistic Regression algorithm was used to train base (or first-level) models and the meta-model (or stacking model).\n",
"2. There was no hyperparameter optimization. Hence, hyperparameters default values were used.\n",
"3. The meta-model learns how to best combine the predictions from the base models.\n",
"4. However, training a meta-model is not the only technique to combine predictions. Model averaging, or soft voting, is also other technique used."
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [],
"source": [
"stacking_results = pd.read_csv(FNAMES['classifier-ensembles']).drop(columns='Id')\n",
"\n",
"stacking_results = stacking_results.astype(\n",
" {'TimeWindow': 'int32', 'NumberOfEstimators': 'int32', 'LogLoss': 'float64'}\n",
" )\n",
"\n",
"stacking_results = stacking_results.sort_values(by='LogLoss')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First of all, let's catch a glimpse of the top ten meta-models, as shown by the following table:"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<center><table>\n",
"<thead>\n",
"<tr><th style=\"text-align: right;\"> Rank</th><th style=\"text-align: right;\"> Time Window</th><th style=\"text-align: right;\"> Number of Estimators</th><th>Stacking Technique </th><th>Cyclical Attributes Encoding Technique </th><th style=\"text-align: right;\"> Log Loss</th></tr>\n",
"</thead>\n",
"<tbody>\n",
"<tr><td style=\"text-align: right;\"> 1</td><td style=\"text-align: right;\"> 336</td><td style=\"text-align: right;\"> 7</td><td>Classifier </td><td>One-hot </td><td style=\"text-align: right;\"> 2.56077</td></tr>\n",
"<tr><td style=\"text-align: right;\"> 2</td><td style=\"text-align: right;\"> 336</td><td style=\"text-align: right;\"> 7</td><td>Soft voting </td><td>One-hot </td><td style=\"text-align: right;\"> 2.56126</td></tr>\n",
"<tr><td style=\"text-align: right;\"> 3</td><td style=\"text-align: right;\"> 168</td><td style=\"text-align: right;\"> 7</td><td>Classifier </td><td>One-hot </td><td style=\"text-align: right;\"> 2.56178</td></tr>\n",
"<tr><td style=\"text-align: right;\"> 4</td><td style=\"text-align: right;\"> 168</td><td style=\"text-align: right;\"> 7</td><td>Soft voting </td><td>One-hot </td><td style=\"text-align: right;\"> 2.56212</td></tr>\n",
"<tr><td style=\"text-align: right;\"> 5</td><td style=\"text-align: right;\"> 168</td><td style=\"text-align: right;\"> 5</td><td>Soft voting </td><td>One-hot </td><td style=\"text-align: right;\"> 2.56214</td></tr>\n",
"<tr><td style=\"text-align: right;\"> 6</td><td style=\"text-align: right;\"> 72</td><td style=\"text-align: right;\"> 5</td><td>Soft voting </td><td>One-hot </td><td style=\"text-align: right;\"> 2.56239</td></tr>\n",
"<tr><td style=\"text-align: right;\"> 7</td><td style=\"text-align: right;\"> 336</td><td style=\"text-align: right;\"> 5</td><td>Classifier </td><td>One-hot </td><td style=\"text-align: right;\"> 2.56241</td></tr>\n",
"<tr><td style=\"text-align: right;\"> 8</td><td style=\"text-align: right;\"> 336</td><td style=\"text-align: right;\"> 5</td><td>Soft voting </td><td>One-hot </td><td style=\"text-align: right;\"> 2.56258</td></tr>\n",
"<tr><td style=\"text-align: right;\"> 9</td><td style=\"text-align: right;\"> 72</td><td style=\"text-align: right;\"> 7</td><td>Soft voting </td><td>One-hot </td><td style=\"text-align: right;\"> 2.56261</td></tr>\n",
"<tr><td style=\"text-align: right;\"> 10</td><td style=\"text-align: right;\"> 168</td><td style=\"text-align: right;\"> 3</td><td>Soft voting </td><td>One-hot </td><td style=\"text-align: right;\"> 2.56267</td></tr>\n",
"</tbody>\n",
"</table></center>"
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"header = [\n",
" 'Rank',\n",
" 'Time Window',\n",
" 'Number of Estimators',\n",
" 'Stacking Technique',\n",
" 'Cyclical Attributes Encoding Technique',\n",
" 'Log Loss'\n",
" ]\n",
"\n",
"table = []\n",
"\n",
"for i, (idx, row) in enumerate(stacking_results.iloc[:10].iterrows()):\n",
" stack_technique = row['StackingMethod']\n",
" stack_technique = (\n",
" 'Classifier' if stack_technique == 'clf' else stack_technique.replace('_', ' ')\n",
" )\n",
" stack_technique = stack_technique[0].upper() + stack_technique[1:]\n",
" \n",
" encoding_technique = row['EncodingOfCyclicalAttributes']\n",
" encoding_technique = (\n",
" 'one-hot' if encoding_technique == 'onehot' else encoding_technique\n",
" )\n",
" encoding_technique = encoding_technique[0].upper() + encoding_technique[1:]\n",
" \n",
" table.append([\n",
" str(i+1),\n",
" str(row['TimeWindow']),\n",
" str(row['NumberOfEstimators']),\n",
" stack_technique,\n",
" encoding_technique,\n",
" '{:.5f}'.format(row['LogLoss'])\n",
" ])\n",
"\n",
"table = ('<center>'\n",
" + tabulate.tabulate(table, header, tablefmt='html')\n",
" + '</center>')\n",
"\n",
"display(HTML(table))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The above results let us draw the following conclusions:\n",
"\n",
"1. Surprisingly, or at least for me, the technique that best handles cyclical attributes is one-hot encoding. Mapping each cyclical attribute onto a circle, through the trigonometric functions sine and cosine, doesn't outperform the basic one-hot encoding.\n",
"2. Model averaging, or soft voting, is a simple but powerful technique to combine predictions, as seven out of the top ten meta-models are built on it.\n",
"3. The larger the number of estimators, the better the discriminative power of the meta-model.\n",
"4. In a like manner, the larger the time window, the better the discriminative power of the meta-model, as eight out of the top ten meta-models use a time window of at least 168 hours to aggregate features.\n",
"5. Finally, it is worth mentioning that the best two meta-models combine the predictions from the same set of seven base models, but differ from each other by the stacking technique."
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"In total, <b>120</b> different base models were trained."
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"total_estimators = np.sum(stacking_results['NumberOfEstimators'].unique()\n",
" * stacking_results['TimeWindow'].unique().shape[0]\n",
" * stacking_results['EncodingOfCyclicalAttributes'].unique().shape[0])\n",
"\n",
"display(HTML('In total, <b>{}</b> different base models were trained.'.format(total_estimators)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Second of all, let's plot the performance of the <b>48</b> meta-models having as splitting criterion the cyclical attributes encoding technique."
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABDsAAAIvCAYAAAB6JNCVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXwV1f3/8deHgBB2FakGlSAiko0LJIilGIIYVgWpiAjK0gYXFEHgC9aNn/i1qNQgrRW/IioUZVOsilaqRAWqlQCRsAsYdhCQLZIACef3x7253uwhiUbD+/l45GHuzJkzZyaJ8+bMmTPmnENEREREREREpLKoUtENEBEREREREREpT+rsEBEREREREZFKRZ0dIiIiIiIiIlKpqLNDRERERERERCoVdXaIiIiIiIiISKWizg4RERERERERqVTU2SEivxpmFmpmzsyq+j5/aGaDKrpdJWFmNXxtv7SU2y8xs37l3S4REZHSMrM0M+vs+/5PZja9jPXlus7/EpnZp2b2R9/3A8xscUW3qaTMbI6ZPVLKbf9gZu+Vd5tEfkrq7BD5iZjZYDNLNbMTZrbPzF40s/o/0b4KDAdm9pqZPVnCOvyBpZD1Hc3sjJml5/m6tqztLy3nXDfn3OvlWacvuOQcW0beYy7PfZ0N51wn59zcitq/iIj8upnZ7WaW7Lue7fXdMPhdedXvnHvKOffH8qqvIL6skpEnh/ztp9xnUZxzs51z8eVdb57jO5PnmAeU9/5Kwjn3inPuxorYt0hpqbND5CdgZqOBp4GxQD2gHdAY+LeZnVeRbSujPc652nm+vqjoRpUnX3Cp7ZyrDXQjzzFXdPtERETOlpk9CEwBngJ+A1wO/B3oVZHtKqUb8+SQ+yq6QeUtT+7YQe5jnl3R7RP5tVBnh0g5M7O6wP8D7nfO/cs5d9o5lwbcCoQCA33lJpjZPDObaWbHzWydmUUH1BNiZm+Z2QEz+9bMRpRD227y7eeIbxhmC9/yWXiDz3u+uwb/U4q6PzWziWa23Hc8i82sQcD635nZf3z73mlmg33L6/nOwQEz225mj5hZFd+6IDObbGYHzWwb0KOAfeYMJR1sZst85Q/7zlm3gLJNzOxzX9s+NrMXzOwfZ38WwcwuM7N/5rTLzO4OWFfVzB73LT9mZivM7OKAzbub2VZfGxMDtrvbzD4xs6m+c7Q1cKSNmX1pZgMD9vG8mR0ysy1mdr+ZZQWU3Rd4t87MJlnA0GIz62Bm//XtZ5WZtS/NeRARkV8+M6sHPAEMd8697Zz7wZdN3nPOjTWzi807CvXCgG1a+67L1XyfE8xsg+8aut7MWhewnwmB19Uirvs9zGy17xq508wmlNNxFpcDLjCzV81sj2/9OwHrEnzX0+/N7F0zCwlYd4OZbTSzo+YdRWJ59xnw2fmu59/4jvsFMzPfuiAz+4svO3xrZvdZKR/Z8dX1qC9rHDSz2RYweti8o3G/9LV5h5ndHrB5AzP7yPezXG5mjX3b5DxuO6yInPJxwOceAcf5XJ6ckjd3XJ0np1xg3uy3z/c78Lj5sp9IedIvlUj5+y1QA3g7cKFzLh34ALghYPFNwBygPvAu8DcA3//w3wO+BhoB1wMjzaxLaRtlZlcBbwIjgYt8bXnPzM5zzt1B7jsHz5RyN7cDQ4CGwHnAGN++GwMfAn/17dsDpPi2+Sve0S9XALHAnb46ABKAnkArIBq4pZj9XwNsAhoAzwCv5IQM4A3gK+BCYAJwR2kO0MyC8J67/wAhQFfgT2YW6yvyENAbiMf7cx0GZAZU0dV3PK2BIWbWMWDddUCyr41/Awp79vk+oBMQCVwLlHguDzMLBd4BHgYuAB4B3jGz80tah4iI/KpcizeXLCxopXNuH/Ap3psyOe4A5jjnTptZX7zXzTuBunizy6GidljMdf8HX1318d7EuMfMepfiuApSVA6YBdQEwvHmlERfWzsBf8Z7/JcA2/FmM8x70+ZtvNfKBsBWoLgbBD2BGCDKV2dOdkvAO2LUgzcDlOWYx+DNGb8DLgVOBxzPlcD7wLN480QbYF3AtrfjzSoXAHvx3qALVFROwbePS4B5wGi8P98DeHNaSc0GjuLNfm3xnotS5TKRoqizQ6T8NQAOOueyCli317c+xzLn3AfOuWy8F+GWvuUxwEXOuSecc6ecc9uAl4Hbitn3QV8P+xEzO4L3gpajH7DIOfdv59xpYDIQjLdzpqRCAuv3fdUKWP+qc26zcy4D70XQ41t+O/Cxc+5N392kQ865FF/HwW3AQ865474RMH/hxwvercAU59xO59z3eMNIUbY75172nc/X8YaW35jZ5XjP6WO+87kMb+dSafwOqOGce9pX12bgVX782fwRGO+c2+KcO+OcW+2cOxKw/VPOuWPOuW+BzwPOEcAm59zMgPY3toLnebkV+Itzbo9z7gDeQFdSg4C3nXMf+9r3AbAeb2gSEZHK50IKzyU5XufHkadBQH+8uQS817VnnHMrnNcW59z2YvZZ4HUfwDn3qXMu1XcNWoP3RkxskbXl9k6eHJIQsK6wHHAJ3o6Gu51zh31t+sy3zQBghnNulXPuJN6OgGt9Nwe6A+uccwt82WkKsK+Y9k1yzh1xzu0AkvjxOn8r8Lxzbpdz7jAw6SyOOa+78WaNPc65TLwdFv18HTt3AO85595yzmU55w44574O2Hae71hP470R5MlTd1E5JceNwArn3Lu+ep4Bvi9Jw30dYdcBDzrnTjjn9gJTKT7jipy1X+xMxyK/YgfxDhGsWkCwuMS3PkfgBfMEUMM3nLExvo6FgPVBwFLwTlwVsDws4PsGgfs0s9cC1oXgvVsBgHPujJntxDtypKT2OOeKeptI3uPJmePiMrx3Q/JqAFQLbJfv+5w2hQA786wrin//zrkTvps5tX37+d45dyKg7E5fu85WYyC0gJ/Nx76Q0YiCjzVfG8l9jgpah2994L7g7M9LoMZAf9+duhzVfHWKiEjlc4jCc0mOfwLTzKwJ0Bw46pz7yreusGt4UQrdxsyuwfsP/Qi8o0CrA/PPou7ezrmPC1lXWA64AG8OOFzANiHAqoDt0s3sEN7rea7rrXPO+bJTUQq7zue9dhdXT4F8WeMy4AMzcwGrquDt2Cru51VUDinJesh/XrLNbHfxrQe8OaQGcODHQTdUAbaUcHuREtPIDpHy9wVwEugTuNDMcia8/KQEdewEvnXO1Q/4quOc6w65J67y3TkoiT14LzA57cm5WOZcnFxBG5WTnUDTApYfxDv0snHAsssD2rSX3B0Sl5dy/3uBC8ysZsCy0nR0gPdYNhbws7nZOefwtr2gYy1PxZ2XH/AO1c0ROGfITmB6nvbXcs4lIiIilVFOLin0sQnf6IB5eEd33MGPozqg8Gt4UYra5g28oysvc87VA6YRMA/GT2Qn3hxQ0GjJvPmoFt5Og93kud4GZKfS2Iv3kZMcpaonIGt0ynMtr+GcO0jpfl5nK+95qULum2fF5ZB04PyAttd1zuWbB0akrNTZIVLOnHNH8Q4n/KuZdTWzar6hkPOAXeQOEIX5CjhuZuPMLNg3EVWEmcWUoWnzgB5mdr15JxwbjTf8/Me3fj/eZyd/CrOBzmZ2q3kn17zQzDy+YabzgP81szq+oY0PAjkTnM0DRpjZpb45JcaXZue+4bbJwAQzO8+8r8st7evTlgGY2UjfZF5VzSzKfpysbTrwlJldYV6tCglXZTEPGGVml/ieJ847oWwK3tEbVc2sHbln238d6Ov7PQjy/X5db7knURURkUrCl0seA14ws95mVtOXTbqZWeBjkDOBwXjn5AjMKtOBMWbWxnddu9J3vS5Kgdd937o6eEdZZJpZW3I/cvuT8D0q8SHwdzM733f81/lWv4l3bgqPmVXH+8aa//oerV0EhJtZH9/I2xHk/of72ZgHPGBmjXy5YFwZDmkaMMnMLgMws4ZmlpNrZgE9zexm37m/yMyiyrCvgrwLxJhZT1+mHIt39EyOFCDOd6znE3CsvsdjvgSe8WW/KmbWzMrxNcgiOdTZIfITcN4JPv+Ed16MY8B/8fZkX+97HrS47bPxTnDlAb7FOwJiOt6JPEvbpk1479j81VffjXgnJD3lK/Jn4BHf869jCqkmxHK/+z3dzH5fgn3vwPvc62i8z3Sm8OP8JPfjvQOwDW9HwhvADN+6l4GP8E7Uuoo8k76epQF4J2k7BDwJzMXb2XNWfM+mdsc718l2vJNyvciPwzwn4Q1HS/D+7KfhHaJbnv6G95GmdXh/t+blWf8nvJOXHsH77PGcgPZvA36Pt0PuoO8YHkDXAxGRSss59xe8NxMewXvd2ol3sut3AsosB84AqwLn5HDOzQf+F+/1+bhvm8B/2Ba0v6Ku+/cCT5jZcbydMHmvYcV5L08OKXDi1QLcgXc06UbgO7wTtuN7JOZR4C28Ixaa4ps/wjdSoi/ea/shoBmw/Czbm+NlYDGwBliNd7LzLCC7FHU9A3wMLPGdx//gnVAU59wWvDc5/oT33CfjnZS13Pg6j27DO4fJAbyvM04OKLII7ySp6/F2bLyTp4r+eCeo3ehr41xfHSLlyrwjoUREzi1mNhfv4yiPV3RbysrMrgbWOuc0D5OIiJSamS0B3nDOFfY2MCkn5n0t7jTnXHGjZH4VzOxL4G/OuX8UW1jkZ6I7eSJyTjCzGDNr6hsu2RXvXY+8dxpERETOSb5HZVvjvcsu5cz32Gh336MljYDHKeR1wCJSPtTZISLniouBT/FOijUVuMc5t7pCWyQiIvILYGav430sYqRz7nhFt6eSMryPkB7G+xjLBryP8YjIT0SPsYiIiIiIiIhIpaKRHSIiIiIiIiJSqaizQ0REREREREQqFc3cX4wGDRq40NDQim6GiIjIL87KlSsPOucuquh2nAuUR0RERApWWB5RZ0cxQkNDSU5OLr6giIjIOcbMtld0G84VyiMiIiIFKyyP6DEWEREREREREalU1NkhIiIiIiIiIpWKOjtEREREREREpFLRnB0iIoU4ffo0u3btIjMzs6KbIlKhatSowaWXXkq1atUquikiIucUZRGRH51tHlFnh4hIIXbt2kWdOnUIDQ3FzCq6OSIVwjnHoUOH2LVrF02aNKno5oiInFOURUS8SpNH9BiLiEghMjMzufDCCxUu5JxmZlx44YW6qygiUgGURUS8SpNH1NkhIlIEhQsR/R2IiFQk/T9YxOts/xbU2SEichbMjNGjR/s/T548mQkTJpRL3YMHD2bBggVlqmPXrl306tWLZs2a0bRpUx544AFOnTrlX9+/f3+ioqJITExk8ODBNGnSBI/Hg8fj4be//W2h9aalpfHGG2/4PycnJzNixIgytTXHa6+9xp49e8qlrkAnTpxgwIABREZGEhERwe9+9zvS09M5cuQIf//730tdb2E/pz/+8Y+sX7++VHWmpaVhZjzyyCP+ZQcPHqRatWrcd999RW776aef8p///KfY9pVUefweiojIT0dZxEtZRFmkOOrsEBE5C9WrV+ftt9/m4MGDFd2UXLKysnDO0adPH3r37s0333zD5s2bSU9P5+GHHwZg3759rFixgjVr1jBq1CgAnn32WVJSUkhJScl1kcorb8CIjo5m6tSp5dL20gSM7OzsYss8//zz/OY3vyE1NZW1a9fyyiuvUK1atTIHjMJMnz6dsLCwUm/fpEkTFi1a5P88f/58wsPDi90ub8AQEZHKTVnES1kkP2WR3NTZISJyFqpWrcqwYcNITEzMty5vL3Tt2rUB7wUgNjaWXr16ccUVVzB+/Hhmz55N27ZtiYyMZOvWrf5tPv74Y6Kjo7nqqqt4//33Ae/FdOzYscTExBAVFcVLL73kr7dDhw7cdNNNhIWFsWTJEmrUqMGQIUMACAoKIjExkRkzZnDixAni4+PZvXs3Ho+HpUuXFnqMn332mf8OS6tWrTh+/Djjx49n6dKleDweEhMT+fTTT+nZsycAEyZMYNCgQXTo0IHGjRvz9ttv8z//8z9ERkbStWtXTp8+DcATTzxBTEwMERERDBs2DOccCxYsIDk5mQEDBuDxeMjIyOCTTz6hVatWREZGMnToUE6ePAlAaGgo48aNo3Xr1syfP5+pU6cSFhZGVFQUt912W77j2Lt3L40aNfJ/bt68OdWrV2f8+PFs3boVj8fD2LFjSU9P5/rrr6d169ZERkbyz3/+07/NzJkziYqKomXLltxxxx359vHoo48yePBgsrOz6dixI8nJyf6f/cMPP0zLli1p164d+/fvB2Dr1q20a9eOyMhIHnnkEf/vCEDNmjVp0aKFv465c+dy6623+tcfOHCA3//+98TExBATE8Py5ctJS0tj2rRpJCYm5vq5fv755/z2t7/liiuu8P9OOucYO3YsERERREZGMnfuXP/y++67j+bNm9O5c2e+++67Qn83RESk4imLKIsEUhYpgnNOX0V8tWnTxonIuWn9+vX5ltWqVcsdPXrUNW7c2B05csQ9++yz7vHHH3fOOTdo0CA3f/78XGWdcy4pKcnVq1fP7dmzx2VmZrqQkBD32GOPOeecmzJlinvggQf823fp0sVlZ2e7zZs3u0aNGrmMjAz30ksvuYkTJzrnnMvMzHRt2rRx27Ztc0lJSa5mzZpu27Ztzjnnnn/+eTdy5Mh8bfZ4PO7rr7923377rQsPD/cvHzRokAsNDXUtW7Z0LVu2dLfffrtzzrmePXu6ZcuWOeecO378uDt9+rRLSkpyPXr08G8b+Pnxxx937du3d6dOnXIpKSkuODjYffDBB84553r37u0WLlzonHPu0KFD/u0HDhzo3n33Xeecc7GxsW7FihXOOecyMjLcpZde6jZt2uScc+6OO+5wiYmJzjnnGjdu7J5++ml/HZdcconLzMx0zjl3+PDhfMe9evVqd9FFF7l27dq5hx9+2G3evNk55/Kdh9OnT7ujR48655w7cOCAa9q0qTtz5oxbu3ata9asmTtw4ECu9uf8nMeMGePuuusud+bMmXzHAfiPb+zYsf6fX48ePdwbb7zhnHPuxRdf9P+O5LTpn//8pxs9erTbsWOH69Spk3v11Vfd8OHDnXPO9e/f3y1dutQ559z27dvd1Vdf7T//zz77bK6f6y233OKys7PdunXrXNOmTZ1zzi1YsMB17tzZZWVluX379rnLLrvM7dmzx7311lv+5bt373b16tXL9Xuco6C/ByDZ/QKu1efCl/KIyLlJWURZRFkkt7PJIxrZISJylurWrcudd955VkMnY2JiuOSSS6hevTpNmzYlPj4egMjISNLS0vzlbr31VqpUqUKzZs244oor2LhxI4sXL2bmzJl4PB6uueYaDh06xDfffANA27Zty/Q60MCho7Nnzwagffv2PPjgg0ydOpUjR45QtWrxbynv1q0b1apVIzIykuzsbLp27Zrv+JKSkrjmmmuIjIxkyZIlrFu3Ll89mzZtokmTJlx11VUADBo0iM8//9y/vl+/fv7vo6KiGDBgAP/4xz8KbKPH42Hbtm2MHTuW77//npiYGDZs2JCvnHOOP/3pT0RFRdG5c2d2797N/v37WbJkCX379qVBgwYAXHDBBf5tJk6cyNGjR5k2bVqBk2Wdd955/rtNbdq08Z+DL774gr59+wJw++2359uua9eu/Pvf/2bOnDm5jhW8d9ruu+8+PB4PN910E8eOHSM9PT1fHQC9e/emSpUqhIWF+e/kLFu2jP79+xMUFMRvfvMbYmNjWbFiBZ9//rl/eUhICJ06dSqwThER+eVQFslPWSQ3ZRE9xiIiUiojR47klVde4YcffvAvq1q1KmfOnAHgzJkzuSbjql69uv/7KlWq+D9XqVKFrKws/7q8FyszwznHX//6V38Q+Pbbb/0BpVatWv6yYWFhrFy5Mtf2x44dY8eOHVx55ZUlPrbx48czffp0MjIyaN++PRs3bix2m8DjqVatmv84co4vMzOTe++9lwULFpCamkpCQkKpXmUaeLyLFi1i+PDhrFq1ipiYmFznMUft2rXp06cPf//73xk4cCAffPBBvjKzZ8/mwIEDrFy5kpSUFH7zm98U27aYmBhWrlzJ999/X+D6wHMQFBRUYNsKct5559GmTRv+8pe/cMstt+Rad+bMGb788kv/78Hu3btzDT0NFPj75r3hISIilY2ySG7KIrkpi6izQ0SkVC644AJuvfVWXnnlFf+y0NBQ/wX+3Xff9T8fejbmz5/PmTNn2Lp1K9u2baN58+Z06dKFF1980V/f5s2bcwWbHNdffz0nTpxg5syZgPf52tGjRzN48GBq1qxZ4jZs3bqVyMhIxo0bR0xMDBs3bqROnTocP378rI8nR84Fu0GDBqSnp+d6njiw7ubNm5OWlsaWLVsAmDVrFrGxsfnqO3PmDDt37iQuLo6nn36ao0ePkp6ezldffcWdd94JwPLlyzl8+DAAp06dYv369TRu3DjfsRw9epSGDRtSrVo1kpKS2L59OwCdOnVi/vz5HDp0CCBXmOjatSvjx4+nR48eZ3Ve2rVrx1tvvQXAnDlzCiwzevRonn766Vx3bwDi4+P561//6v+ckpICUOKfTYcOHZg7dy7Z2dkcOHCAzz//nLZt23Ldddf5l+/du5ekpKQSH4+IiFQcZZGzoyzidS5lEXV2iIiU0ujRo3PNhJ6QkMBnn31Gy5Yt+eKLL3L1/JfU5ZdfTtu2benWrRvTpk2jRo0a/PGPfyQsLIzWrVsTERHBXXfdVWDvvJmxcOFC5s+fT7NmzbjqqquoUaMGTz31VKH7Gzt2rH8CMI/Hw6lTp5gyZQoRERFERUVRrVo1unXrRlRUFEFBQbRs2bLACdGKU79+fRISEoiIiKBLly7ExMT41w0ePJi7774bj8eDc45XX32Vvn37EhkZSZUqVbj77rvz1Zednc3AgQOJjIykVatWjBgxgvr167Njxw6Cg4MBb1CKjY31l4mOjub3v/89F154Ie3btyciIoKxY8cyYMAAkpOTiYyMZObMmVx99dUAhIeH8/DDDxMbG0vLli158MEHc7Whb9++JCQkcNNNN5GRkVGi8zBlyhSee+45oqKi2LJlC/Xq1ctXJjw8nEGDBuVbPnXqVJKTk4mKiiIsLIxp06YBcOONN7Jw4cJiJ3u7+eab/ROcderUiWeeeYaLL76Ym2++mWbNmhEWFsadd97JtddeW6JjERGRiqcsUnLKIl7nUhYxDW8tWnR0tMuZjVZEzi0bNmygRYsWFd0MOQtjx47ljjvuICoqqqKbUqATJ04QHByMmTFnzhzefPPNXDOu/5IV9PdgZiudc9EV1KRzivKIyLlJWeTXR1nkp3U2eaT4mV5ERER+JZ599tmKbkKRVq5cyX333Ydzjvr16zNjxoyKbpKIiIiUI2WRXw51doiIiPxMOnTowNdff13RzRAREZFz1LmURTRnh4iIiIiIiIhUKursEBEREREREZFKRZ0dIiIiIiIiIlKpqLNDRERERERERCoVdXaIiIiIiIiISKWit7GIiFQyoeMXlWt9aZN6FFsmIyODrl27smTJEnbt2sXNN9/MmTNnOH36NPfffz933313sXVMnjyZ6dOnU6NGDapVq8b999/PnXfeSceOHZk8eTLR0flen37WkpOTmTlzJlOnTuXkyZP06NGDgwcP8tBDD/Hvf/+bBx98kLCwsBLXN2bMGLp3706nTp3K3DYREZHKRHmkcMojPw91doiISJnNmDGDPn36EBQUxCWXXMIXX3xB9erVSU9PJyIigptuuomQkJBCt582bRr//ve/+eqrr6hbty7Hjh1j4cKF5d7O6Ohof0hZvXo1ACkpKQD069fvrOrKzs7m/vvvJyEhQeFCRETkF0B5RHkkkB5jERGRMps9eza9evUC4LzzzqN69eoAnDx5kjNnzhS7/VNPPcWLL75I3bp1Aahbty6DBg3KV+6ee+4hOjqa8PBwHn/8cf/y8ePHExYWRlRUFGPGjAFg/vz5RERE0LJlS6677joAPv30U3r27Ml3333HwIEDWbFiBR6Ph61bt9KxY0eSk5MBWLx4Mddeey2tW7emb9++pKenAxAaGsq4ceNo3bo18+fPp3Hjxhw6dIh9+/aV9tSJ/GKFjl9U7JeIyC+J8ojySCCN7BARkTI5deoU27ZtIzQ01L9s586d9OjRgy1btvDss88WeRfl2LFjHD9+nCuuuKLYff3v//4vF1xwAdnZ2Vx//fWsWbOGRo0asXDhQjZu3IiZceTIEQCeeOIJPvroIxo1auRflqNhw4ZMnz6dyZMn8/777+dad/DgQZ588kk+/vhjatWqxdNPP81zzz3HY489BsCFF17IqlWr/OVbt27N8uXL+f3vf19s+0VEROSnoTyiPJKXRnaIiEiZHDx4kPr16+dadtlll7FmzRq2bNnC66+/zv79+8tlX/PmzaN169a0atWKdevWsX79eurVq0eNGjX4wx/+wNtvv03NmjUBaN++PYMHD+bll18mOzu7xPv48ssvWb9+Pe3bt8fj8fD666+zfft2//q8w0sbNmzInj17yuX4REREpHSUR5RH8lJnh4iIlElwcDCZmZkFrgsJCSEiIoKlS5cWun3dunWpXbs227ZtK3I/3377LZMnT+aTTz5hzZo19OjRg8zMTKpWrcpXX33FLbfcwvvvv0/Xrl0B73O3Tz75JDt37qRNmzYcOnSoRMfjnOOGG24gJSWFlJQU1q9fzyuvvOJfX6tWrVzlMzMzCQ4OLlHdIiIi8tNQHlEeyUudHSIiUibnn38+2dnZ/oCxa9cuMjIyADh8+DDLli2jefPmANx555189dVX+ep46KGHGD58OMeOHQMgPT2dmTNn5ipz7NgxatWqRb169di/fz8ffvihv+zRo0fp3r07iYmJfP311wBs3bqVa665hieeeIKLLrqInTt3luh42rVrx/Lly9myZQsAP/zwA5s3by60/ObNm4mIiChR3SIiIvLTUB5RHslLc3aIiFQyJXk1W3mLj49n2bJldO7cmQ0bNjB69GjMDOccY8aMITIyEoA1a9YU+LzsPffcQ3p6OjExMVSrVo1q1aoxevToXGVatmxJq1atuPrqq7nsssto3749AMePH+Ota+8AACAASURBVKdXr15kZmbinOO5554DYOzYsXzzzTc457j++utp2bIln332WbHHctFFF/Haa6/Rv39/Tp48CcCTTz7JVVddla/s6dOn2bJlS7m8hk5ERKQyUR5RHqlo5pyr6Db8okVHR7uc2XBF5NyyYcMGWrRoUdHN+FVYtWoViYmJzJo1q9Ayx44d4w9/+APz58//GVv201q4cCGrVq1i4sSJFd2Un1xBfw9mttI5p2T1M6iIPFKSt61UxD9mRM4lyiJnR3lEeSSQHmMREZEya926NXFxcUVOvFW3bt1KFSwAsrKy8t3xERERkYqhPCKB9BiLiIiUi6FDh1Z0E352ffv2regmiIiISADlEcmhkR0iIiIiIiIiUqmos0NEREREREREKhV1doiIiIiIiIhIpaLODhERKbOMjAxiY2PJzs5m+/bttG7dGo/HQ3h4ONOmTSt2+y+//JJrrrkGj8dDixYtmDBhQpHlT548SefOnfF4PMydO5cpU6Zw4sSJcjmWp556Ktfn3/72t6Wu6/333+exxx4ra5NERESkBJRHCnau5hG9erYYevWsyLnrV/u6twn1yrm+o8UWeeGFF8jKyuKBBx7g1KlTOOeoXr066enpRERE8J///KfA99nnaN68OfPmzaNly5ZkZ2ezadMmwsLCCi3/5Zdf8sgjj/Dxxx8DEBoaSnJyMg0aNDj748ujdu3apKenl7keAOccrVu3Zvny5dSsWbNc6qwoevVsxdKrZ0XOTb/aLALKI2WkPFIwvXq2GGZWy8xeN7OXzWxARbdHROTXbvbs2fTq1QuA8847j+rVqwPeOx5nzpwpdvvvvvuOSy65BICgoCB/sPj+++/p3bs3UVFRtGvXjjVr1vDdd98xcOBAVqxYgcfj4fnnn2fPnj3ExcURFxeXq95//etfuWYo//TTT+nZsycAb775JpGRkURERDBu3DgAxo8fT0ZGBh6PhwEDvJeH2rVr+7ft2LEjt9xyC1dffTUDBgwg54bBBx98wNVXX02bNm0YMWKEfx9mRseOHXn//fdLcValslMeEREpX8ojyiOBiu3sMLPLzCzJzNab2Toze6CQcmlmlmpmKWaW7FvW3Pc55+uYmY0M2Ka+mS0ws41mtsHMri3NQZjZDDP7zszWFrCuq5ltMrMtZjbet7gPsMA5lwDcVJp9ioiI16lTp9i2bRuhoaH+ZTt37iQqKorLLruMcePGFXkXBWDUqFE0b96cm2++mZdeeonMzEwAHn/8cVq1asWaNWt46qmnuPPOO2nYsCHTp0+nQ4cOpKSk8MADDxASEkJSUhJJSUm56u3cuTP//e9/+eGHHwCYO3cut912G3v27GHcuHEsWbKElJQUVqxYwTvvvMOkSZMIDg4mJSWF2bNn52vn6tWrmTJlCuvXr2fbtm0sX76czMxM7rrrLj788ENWrlzJgQMHcm0THR3N0qVLS3NqJYDyiIiIFEV5RHkkr5KM7MgCRjvnwoB2wHAzK2wsT5xzzpMzhMQ5t8n32QO0AU4ACwPKPw/8yzl3NdAS2BBYmZk1NLM6eZZdWcB+XwO65l1oZkHAC0A3IAzo72v7pcBOX7HsQo9cRESKdfDgQerXr59r2WWXXcaaNWvYsmULr7/+Ovv37y+yjscee4zk5GTi4+N544036NrV+7/0ZcuWcccddwDQqVMnDh06xLFjx0rctqpVq9K1a1fee+89srKyWLRoEb169WLFihV07NiRiy66iKpVqzJgwAA+//zzYutr27Ytl156KVWqVMHj8ZCWlsbGjRu54ooraNKkCQD9+/fPtU3Dhg3Zs2dPidsshVIeERGRQimPKI/kVWxnh3Nur3Nule/743gDQKNS7Ot6YKtzbjuAmdUDrgNe8dV9yjl3JM82scA7Zlbdt00C8NcC2vg58H0B+2wLbHHObXPOnQLmAL2AXXgDBpyjj/KIiJSX4OBg/52PvEJCQoiIiCjRnYSmTZtyzz338Mknn/D1119z6NChcmnfbbfdxrx581iyZAnR0dHUqVOn+I0KkTMcFrzDW7OysordJjMzk+Dg4FLvU7yUR0REpCjKI0U7F/PIWV1YzSwUaAX8t4DVDlhsZivNbFgB628D3gz43AQ4ALxqZqvNbLqZ1cpVoXPzgY+Aub5nWYcCfSm5Rvx4xwS8oaIR8DbwezN7EXivoA3N7EYz+7+jR4ufCEdE5Fx2/vnnk52d7Q8Yu3btIiMjA4DDhw+zbNkymjdvDsCdd97JV199la+ORYsW+Z83/eabbwgKCqJ+/fp06NDBP3zz008/pUGDBtStWzff9nXq1OH48eMFti82NpZVq1bx8ssvc9tttwHeOyKfffYZBw8eJDs7mzfffJPY2FgAqlWrxunTp0t8/M2bN2fbtm2kpaUB3qGpgTZv3kxERESJ65PiKY+IiEheyiPKI3mVuLPDzGoDbwEjnXMFjdn5nXOuNd4hmsPN7LqAbc/D+yzq/IDyVYHWwIvOuVbAD8B48nDOPQNkAi8CNznnyjwlrXPuB+fcEOfcPc65/A9Becu855wbVq9eOc8iLCJSCcXHx7Ns2TLAO0v2NddcQ8uWLYmNjWXMmDFERkYCsGbNmgKfl501axbNmzfH4/Fwxx13MHv2bIKCgpgwYQIrV64kKiqK8ePH8/rrrxe4/2HDhtG1a9d8E4KB945Hz549+fDDD/0TdV1yySVMmjSJuLg4WrZsSZs2bfwTmg0bNoyoqCj/hGDFCQ4O5u9//ztdu3alTZs21KlTh8BrR1JSEj166I0V5UV5RERECqM8ojwSqESvnjWzasD7wEfOuedKUH4CkO6cm+z73AsY7pyLDyhzMfClcy7U97kDMN451yNPXR3wBouVwHHn3H2F7DMUeN85FxGw7FpggnOui+/zQwDOuT8Xe9A+evWsyLnrV/26t5/ZqlWrSExMZNasWYWWOXbsGH/4wx+YP39+oWV+rdLT06lduzbOOYYPH06zZs0YNWoU+/fv5/bbb+eTTz6p6CaW2S/h1bPKI3r1rMi5Rlnk7CiPKI8EKsnbWAzvc6wbCgsWvlen1cn5HogHAmci70/uIaM45/YBO82suW/R9cD6PPW2Av4P73OtQ4ALzezJ4tocYAXQzMya+O7m3Aa8exbbi4hICbRu3Zq4uDiyswufY7Fu3bqVMlgAvPzyy3g8HsLDwzl69Ch33XUXADt27OAvf/lLBbeuclAeERGR4iiPKI8EKnZkh5n9DlgKpAI5Lyf+k3PuAzP7APgjUIMfZzWvCrzhnPtf3/a1gB3AFc65o3nq9gDTgfOAbcAQ59zhgPXtgWPOuVTf52rAYOfcy3nqeRPoCDQA9gOPO+de8a3rDkwBgoAZOe0qKY3sEDl36W6KyI8qemSH8ohGdoici5RFRHI7mzxStbjKnHPLACtkXfeAjy0LKfMDcGEh61KAQkOSc255ns+ngZcLKNc/77KAdR8AHxS2XkRERH75lEdERETkbOg1ZyIiIiIiIiJSqaizQ0REREREREQqFXV2iIiIiIiIiEilos4OEREREREREalUip2gVEREfl0iX48s1/pSB6UWWyYjI4OuXbuyZMkSgoKCCAoKIjLS247LL7+cd98t/i2bM2fO5JlnnsHMqFq1KgMGDGDMmDEMHjyYnj17csstt5T5WPbs2cOIESNYsGABAP3792fdunUMGTKEw4cPc91119G5c+cS1/e3v/2NmjVrMnTo0DK3TUREpDJRHimc8sjPQ50dIiJSZjNmzKBPnz4EBQUBEBwcTEpKSom3//DDD5kyZQqLFy8mJCSEkydPMnPmzHJvZ0hIiD9Y7Nu3jxUrVrBly5ZS1ZWVlcXQoUNp3769woWIiMgvgPKI8kggPcYiIiJlNnv2bHr16lXq7f/85z8zefJkQkJCAKhevToJCQn5yj3xxBPExMQQERHBsGHDcM4BMHXqVMLCwoiKiuK2224D4LPPPsPj8eDxeGjVqhXHjx8nLS2NiIgIAOLj49m9ezcej4elS5cyePBgf/BYuXIlsbGxtGnThi5durB3714AOnbsyMiRI4mOjub555+nZs2ahIaG8tVXX5X62EVERKR8KI8ojwRSZ4eIiJTJqVOn2LZtG6Ghof5lmZmZREdH065dO955551i61i7di1t2rQpttx9993HihUrWLt2LRkZGbz//vsATJo0idWrV7NmzRqmTZsGwOTJk3nhhRdISUlh6dKlBAcH56rr3XffpWnTpqSkpNChQwf/8tOnT3P//fezYMECVq5cydChQ3n44YdzHW9ycjKjR48GIDo6mqVLlxbbdhEREfnpKI8oj+Slx1hERKRMDh48SP369XMt2759O40aNWLbtm106tSJyMhImjZtWuZ9JSUl8cwzz3DixAm+//57wsPDufHGG4mKimLAgAH07t2b3r17A9C+fXsefPBBBgwYQJ8+fbj00ktLtI9Nmzaxdu1abrjhBgCys7O55JJL/Ov79euXq3zDhg3ZuHFjmY9NRERESk95RHkkL43sEBGRMgkODiYzMzPXskaNGgFwxRVX0LFjR1avXl1kHeHh4axcubLIMpmZmdx7770sWLCA1NRUEhIS/PtdtGgRw4cPZ9WqVcTExJCVlcX48eOZPn06GRkZtG/fvsQBwDlHeHg4KSkppKSkkJqayuLFi/3ra9Wqla9dee/SiIiIyM9LeUR5JC91doiISJmcf/75ZGdn+y/0hw8f5uTJk4D3Lsvy5csJCwsD4KGHHmLhwoX56njooYcYO3Ys+/btA7xDM6dPn56rTE79DRo0ID093f8865kzZ9i5cydxcXE8/fTTHD16lPT0dLZu3UpkZCTjxo0jJiamxOGiefPmHDhwgC+++ALwDiNdt25doeU3b97sf+5WREREKobyiPJIXnqMRUSkkinJq9nKW3x8PMuWLaNz585s2LCBu+66iypVqnDmzBnGjx/vDxepqancdNNN+bbv3r07+/fvp3PnzjjnMLN8M4rXr1+fhIQEIiIiuPjii4mJiQG8wzoHDhzI0aNHcc4xYsQI6tevz6OPPkpSUhJVqlQhPDycbt26+Sf2Ksp5553HggULGDFiBEePHiUrK4uRI0cSHh5eYPnly5czYcKEszxjIiIilZvyiPJIRbOcmWOlYNHR0S45ObmimyEiFWDDhg20aNGiopvxq7Bq1SoSExOZNWtWkeW6dOnCRx999DO16qe3evVqnnvuuWKPuzIo6O/BzFY656IrqEnnlIrII6HjFxVbJm1Sj5+hJSLnLmWRs6M8ojwSSI+xiIhImbVu3Zq4uDiys7OLLFeZggV4h8VOnDixopshIiIiKI9IbnqMRUREykXeYZ7ngpwZ0kVEROSXQXlEcmhkh4iIiIiIiIhUKursEBEREREREZFKRZ0dIiIiIiIiIlKpqLNDRETKLCMjg9jYWP+EYEFBQXg8HjweT4Gvdstr06ZNdOzYEY/HQ4sWLRg2bFix2/Tv35+oqCgSExN57bXX2LNnT5mPA2DKlCmcOHHC/7l79+4cOXKkVHWlpqYyePDgcmmXiIiIFE15pGDnah7RBKUiIpXMhqvL9xV1LTZuKLbMjBkz6NOnD0FBQQAEBweTkpJS4n2MGDGCUaNG0atXL8B7US7Kvn37WLFiBVu2bAGgY8eOREREEBISUuJ9FmbKlCkMHDiQmjVrAvDBBx+Uuq7IyEh27drFjh07uPzyy8vcNhERkV8L5ZGyUR4pO43sEBGRMps9e7Y/GJTG3r17ufTSS/2fIyMjAcjMzGTIkCFERkbSqlUrkpKSAIiPj2f37t14PB4mTpxIcnIyAwYMwOPxkJGR4a9n48aNtG3b1v85LS3NX/cnn3xCq1atiIyMZOjQoZw8eZKpU6eyZ88e4uLiiIuLAyA0NJSDBw+SlpZGixYtSEhIIDw8nPj4eP++VqxYQVRUFB6Ph7FjxxIREeHf54033sicOXNKfW5ERESkZJRHlEcCqbNDRETK5NSpU2zbto3Q0FD/sszMTKKjo2nXrh3vvPNOsXWMGjWKTp060a1bNxITE/3DNF944QXMjNTUVN58800GDRpEZmYm7777Lk2bNiUlJYVHH32U6OhoZs+eTUpKCsHBwf56r776ak6dOsW3334LwNy5c+nXrx+ZmZkMHjyYuXPnkpqaSlZWFi+++CIjRowgJCSEpKQkf5AJ9M033zB8+HDWrVtH/fr1eeuttwAYMmQIL730EikpKf67STmio6NZunTpWZ9XERERKTnlEeWRvNTZISIiZXLw4EHq16+fa9n27dtJTk7mjTfeYOTIkWzdurXIOoYMGcKGDRvo27cvn376Ke3atePkyZMsW7aMgQMHAt6g0LhxYzZv3nxW7bv11luZO3cu8GO42LRpE02aNOGqq64CYNCgQXz++efF1tWkSRM8Hg8Abdq0IS0tjSNHjnD8+HGuvfZaAG6//fZc2zRs2LDcnt8VERGRgimPKI/kpc4OEREpk+DgYDIzM3Mta9SoEQBXXHEFHTt2ZPXq1cXWExISwtChQ/nnP/9J1apVWbt2bbm0r1+/fsybN4/NmzdjZjRr1qzUdVWvXt3/fVBQEFlZWcVuk5mZmevujoiIiJQ/5ZGinYt5RJ0dIiJSJueffz7Z2dn+gHH48GFOnjwJeO+yLF++nLCwMAAeeughFi5cmK+Of/3rX5w+fRrwTvZ16NAhGjVqRIcOHZg9ezYAmzdvZseOHTRv3jzf9nXq1OH48eMFtq9p06YEBQUxceJE+vXrB0Dz5s1JS0vzTyg2a9YsYmNji62rIPXr16dOnTr897//Bcj3POzmzZtzPTMrIiIi5U95RHkkL3V2iIhImcXHx7Ns2TIANmzYQHR0NC1btiQuLo7x48f7w0VqaioXX3xxvu0XL15MREQELVu2pEuXLjz77LNcfPHF3HvvvZw5c4bIyEj69evHa6+9lutuRo7Bgwdz991355sQLEe/fv34xz/+wa233gpAjRo1ePXVV+nbty+RkZFUqVKFu+++G4Bhw4bRtWtX/4RgJfHKK6+QkJCAx+Phhx9+oF69ev51SUlJ9OjRo8R1iYiISOkojyiPBDLnXEW34RctOjraJScnV3QzRKQCbNiwgRYtyve1aZXVqlWrSExMZNasWUWW69KlCx999NHP1KqfT3p6OrVr1wZg0qRJ7N27l+eff56TJ08SGxvLsmXLqFr11/2294L+HsxspXMuuoKadE6piDwSOn5RsWXSJp1bwVnk56YscnaUR5RHAv26j1RERH4RWrduTVxcHNnZ2flm/w5UGYMFwKJFi/jzn/9MVlYWjRs35rXXXgNgx44dTJo06VcfLERERH4NlEeURwKdW0crIiI/maFDh1Z0EypMv379/M/fBmrWrFmZJiATERGRs6M8ojySQ3N2iIiIiIiIiEilos4OEREREREREalU1NkhIiIiIiIiIpWKOjtEREREREREpFJRZ4eIiIiIiIiIVCp6G4uISCXzwt1LyrW+4dM6FVsmIyODrl27smTJEj7//HNGjRrlX7dx40bmzJlD7969i6xj5syZPPPMM5gZVatWZcCAAYwZM4bBgwfTs2dPbrnlljIfy549exgxYgQLFiwAoH///qxbt44hQ4Zw+PBhrrvuOjp37lzi+v72t79Rs2bNc3rmdxERkYIojxROeeTnoc4OEREpsxkzZtCnTx+CgoKIi4sjJSUFgO+//54rr7yS+Pj4Irf/8MMPmTJlCosXLyYkJISTJ08yc+bMcm9nSEiIP1js27ePFStWsGXLllLVlZWVxdChQ2nfvr3ChYiIyC+A8ojySCA9xiIiImU2e/ZsevXqlW/5ggUL6NatGzVr1ixy+z//+c9MnjyZkJAQAKpXr05CQkK+ck888QQxMTFEREQwbNgwnHMATJ06lbCwMKKiorjtttsA+Oyzz/B4PHg8Hlq1asXx48dJS0sjIiICgPj4eHbv3o3H42Hp0qUMHjzYHzxWrlxJbGwsbdq0oUuXLuzduxeAjh07MnLkSKKjo3n++eepWbMmoaGhfPXVV6U8cyIiIlJelEeURwKps0NERMrk1KlTbNu2jdDQ0Hzr5syZQ//+/YutY+3atbRp06bYcvfddx8rVqxg7dq1ZGRk8P777wMwadIkVq9ezZo1a5g2bRoAkydP5oUXXiAlJYWlS5cSHBycq653332Xpk2bkpKSQocOHfzLT58+zf3338+CBQtYuXIlQ4cO5eGHH851vMnJyYwePRqA6Oholi5dWmzbRURE5KejPKI8kpc6O0REpEwOHjxI/fr18y3fu3cvqampdOnSpdz2lZSUxDXXXENkZCRLlixh3bp1AERFRTFgwAD+8Y9/ULWq9wnN9u3b8+CDDzJ16lSOHDniX16cTZs2sXbtWm644QY8Hg9PPvkku3bt8q/v169frvINGzZkz5495XSEIiIiUhrKI8ojeamzQ0REyiQ4OJjMzMx8y+fNm8fNN99MtWrViq0jPDyclStXFlkmMzOTe++9lwULFpCamkpCQoJ/v4sWLWL48OGsWrWKmJgYsrKyGD9+PNOnTycjI4P27duzcePGEh2Pc47w8HBSUlJISUkhNTWVxYsX+9fXqlUrX7vy3qURERGRn5fyiPJIXursEBGRMjn//PPJzs7OFzDefPPNfENGH3roIRYuXJivjoceeoixY8eyb98+wDs0c/r06bnK5NTfoEED0tPT/c+znjlzhp07dxIXF8fTTz/N0aNHSU9PZ+vWrURGRjJu3DhiYmJKHC6aN2/OgQMH+OKLLwDvMNKcOzYF2bx5s/+5WxEREakYyiPKI3npbSwiIpVMSV7NVt7i4+NZtmyZ/zVpaWlp7Ny5k9jY2FzlUlNTuemmm/Jt3717d/bv30/nzp1xzmFm+WYUr1+/PgkJCURERHDxxRcTExMDQHZ2NgMHDuTo0aM45xgxYgT169fn0UcfJSkpiSpVqhAeHk63bt38E3sV5bzzzmPBggWMGDGCo0ePkpWVxciRIwkPDy+w/PLly5kwYUJJTpOIiMg5Q3lEeaSiWc7MsVKw6Ohol5ycXNHNEJEKsGHDBlq0aFHRzfhVWLVqFYmJicyaNavIcl26dOGjjz76mVr101u9ejXPPfdcscddGRT092BmK51z0RXUpHNKReSR0PGLii2TNqnHz9ASkXOXssjZUR5RHgmkx1hERKTMWrduTVxcHNnZ2UWWq0zBAryToU2cOLGimyEiIiIoj0hueoxFRETKRd5hnueCG264oaKbICIiIgGURyTHOdnZYWa1gL8Dp4BPnXOzK7hJIiIico5RHhEREfnpFPsYi5ldZmZJZrbezNaZ2QOFlEszs1QzSzGzZN+y5r7POV/HzGxknu2CzGy1mb1f2oMwsxlm9p2ZrS1gXVcz22RmW8xsvG9xH2CBcy4ByD8zjYiIiPyiKI+IiIjI2SjJnB1ZwGjnXBjQDhhuZmGFlI1zznlyJgdxzm3yffYAbYATQN53/DwAbCioMjNraGZ18iy7soCirwFdC9g+CHgB6AaEAf19bb8U2OkrVvQDXSIiIvJLoDwiIiIiJVZsZ4dzbq9zbpXv++N4g0CjUuzremCrc257zgIzuxToAUwvZJtY4B0zq+4rnwD8tYA2fg58X8D2bYEtzrltzrlTwBygF7ALb8CAQs6Bmd1oZv939OjRkhybiMg5LSMjg9jYWLKzs0lKSsLj8fi/atSowTvvvFPk9ps2baJjx454PB5atGjBsGHDit1n//79iYqKIjExkddee409e/aUy7FMmTKFEydO+D93796dI0eOlKqu1NRUBg8eXC7tOtcpjyiPiIgUR3mkYOdqHjmrOTvMLBRoBfy3gNUOWGxmDnjJOfd/edbfBryZZ9kU4H+AOhTAOTffzJoAc81sPjAUOJvZVxrx4x0T8IaKa4CpwN/MrAfwXiH7fg94Lzo6OuEs9iciUuH+0q9nudY3em7xo/pnzJhBnz59CAoKIi4ujpSUFAC+//57rrzySuLj44vcfsSIEYwaNYpevXoB3otyUfbt28eKFSvYsmULAB07diQiIoKQkJCSHFKRpkyZwsCBA6lZsyYAH3zwQanrioyMZNeuXezYsYPLL7+8zG0TL+UREZFfPuWRslEeKbsSv3rWzGoDbwEjnXPHCijyO+dca7xDNIeb2XUB256H91nU+QHLegLfOedWFrVf59wzQCbwInCTcy69pG0uos4fnHNDnHP3aDIwEZGymz17tj8YBFqwYAHdunXzX6gLs3fvXi699FL/58jISAAyMzMZMmQIkZGRtGrViqSkJADi4+PZvXs3Ho+HiRMnkpyczIABA/B4PGRkZPjr2bhxI23btvV/TktL89f9ySef0KpVKyIjIxk6dCgnT55k6tSp7Nmzh7i4OOLi4gAIDQ3l4MGDpKWl0aJFCxISEggPDyc+Pt6/rxUrVhAVFYXH42Hs2LFERET493njjTcyZ86cszqfUjjlkV+nyNcji/0SESkr5RHlkUAl6uwws2p4g8Vs59zbBZVxzu32/fc7vM/Btg1Y3Q1Y5ZzbH7CsPXCTmaXhHc7Zycz+UcC+OwARvjofL0l7A+wGLgv4fKlvmYiIlJNTp06xbds2QkND862bM2cO/fv3L7aOUaNG0alTJ7p160ZiYqJ/mOYLL7yAmZGamsqbb77JoEGDyMzM5N1336Vp06akpKTw6KOPEh0dzezZs0lJSSE4ONhf79VXX82pU6f49ttvAZg7dy79+vUjMzOTwYMHM3fuXFJTU8nKyuLFF19kxIgRhISEkJSU5A8ygb755huG/3/27j7as7q+D/37I4MQHjqktGhkiEwVeaiUB4nV+kDERhEVbzQwEqMGvHBJ0epd3HbhvavL9MYk1FRzvQZJURSrXAQRDQgLSK9QgrGIIDrAAAWLMhSLD1cEosLA5/7x+40cD+fM+Z0ZZn5n9rxea82atb0FpQAAIABJREFUs/f+7r0/v3PWcD6893fvfcopueWWW7Lbbrvl85//fJLk+OOPz3/4D/8hN910U7bbbrtf2uewww7L3/zN30z8/WR++hEA5qMf0Y/MNsnbWCrJ2UnWdPeH5hmz8/oHd41fo/aqJDOfRH5cZk0Z7e73dveK7t47oymlX+7u35t13EOSnJXRfa3HJ9m9qt4/4WdLkuuT7FNVK8dXc96c5OJF7A/AAn7wgx9kt912e9L6++67L6tXr86rX/3qBY9x/PHHZ82aNTnmmGNy9dVX50UvelF+/vOf59prr83v/d7oV8N+++2XZz/72bnjjjsWVd+xxx6b888/P8kTzcXtt9+elStX5nnPe16S5O1vf3uuueaaBY+1cuXKHHzwwUmSF7zgBbn77rvz4x//OA8++GBe/OIXJ0l+93d/95f22WOPPZ6y+3e3ZfoRADZEP6IfmW2SmR0vSfLWjK50rH9l21FJUlWXVdWzkjwjybVV9c0kX0tyaXdfPh6zc0b3tc55BWYBOyU5trvv6u7Hk7wtyXdmD6qq85J8Ncm+VbW2qt6RJN29Lsk7k1yR0YPMLujuWzaiDgDm8Su/8iv52c9+9qT1F1xwQX77t38722+//UTHedaznpUTTjghf/VXf5Vly5bl5puf9PbOjbJq1apccMEFueOOO1JV2WeffTb6WDvssMMvvt5uu+2ybt26Bff52c9+9ktXd9ho+hEA5qUf2bBtsR+Z5G0s13Z3dfc/Wf/atu6+bLztqO7+7+Onix80/vOPu/uPZ+z/cHfv3t3zPka8u6/u7ic9waa7v9Ldq2csP9rdH5tj3HHd/Wvdvf346szZM7Zd1t3P6+7nzKwLgKfGr/7qr+axxx57UoNx3nnnPWnK6Hvf+9584Quz3/iZXH755Xn00UeTjB729cMf/jB77rlnXvayl+Xcc0ePMrjjjjvy3e9+N/vuu++T9t91113z4IMPzlnfc57znGy33Xb5oz/6o6xatSpJsu++++buu+/+xQPFPv3pT+fwww9f8Fhz2W233bLrrrvmuutGz8qcfT/sHXfc8Uv3zLJx9CMAbIh+RD8y28QPKAWA+bzqVa/Ktdde+4vlu+++O/fcc88vfmGvt3r16jzzmc980v5XXnllnv/85+eggw7Kq1/96vzZn/1ZnvnMZ+Zf/It/kccffzwHHnhgVq1alXPOOeeXrmas9/u///s5+eSTn/RAsPVWrVqVz3zmMzn22GOTJDvuuGM++clP5phjjsmBBx6Ypz3taTn55JOTJCeddFKOPPLIXzwQbBJnn312TjzxxBx88MF5+OGHs3z58l9su+qqq/La17524mMBABtHP6Ifmam6e9o1LGmHHXZYf/3rX592GcAUrFmzJvvvv/+0y9gq3HjjjfnzP//zfPrTn97guFe/+tW54oortlBVW85DDz2UXXbZJUly+umn57777suHP/zh/PznP8/hhx+ea6+9NsuWLept70vOXP8equqG7j5sSiVtU6bRj+x92qULjrn79Mka50netrL67Rt+xSNsi/Qii6Mf0Y/MtHV/UgCWhEMPPTSveMUr8thjjz3p6d8zDbGxSJJLL700f/qnf5p169bl2c9+ds4555wkyXe/+92cfvrpW31jAQBbA/2IfmSmbevTArDZnHDCCdMuYWpWrVr1i/tvZ9pnn3026QFkAMDi6Ef0I+t5ZgcAAAAwKMIOgA3wXCPw7wBgmvw3GEYW+29B2AEwjx133DE//OEPNRls07o7P/zhD7PjjjtOuxSAbY5eBEY2ph/xzA6AeaxYsSJr167N97///WmXAlO14447ZsWKFdMuA2CboxeBJyy2HxF2AMxj++23z8qVK6ddBgCwjdKLwMZzGwsAAAAwKMIOAAAAYFCEHQAAAMCgCDsAAACAQRF2AAAAAIMi7AAAAAAGRdgBAAAADIqwAwAAABgUYQcAAAAwKMumXQAAALB57H3apQuOufv0126BSgC2LDM7AAAAgEERdgAAAACDIuwAAAAABkXYAQAAAAyKsAMAAAAYFGEHAAAAMCjCDgAAAGBQhB0AAADAoAg7AAAAgEERdgAAAACDIuwAAAAABkXYAQAAAAyKsAMAAAAYFGEHAAAAMCjLpl0AAAAA89v7tEsnGnf36a/dzJXA1sPMDgAAAGBQzOwAAGCrccbJX15wzCl/ecQWqASApczMDgAAAGBQhB0AAADAoAg7AAAAgEERdgAAAACDIuwAAAAABkXYAQAAAAyKV88uQXufdumCY+4+/bVboBIAAADY+gg7AADYOH+4fLJxK39989YxywdXvW6icaee/6XNXAkA0+I2FgAAAGBQtsmZHVW1c5KPJnkkydXdfe6USwIAtjH6EQDYfBac2VFVe1XVVVV1a1XdUlXvnmfc3VW1uqpuqqqvj9ftO15e/+cnVfWexRx3ElX1iaq6v6punmPbkVV1e1XdWVWnjVe/McmF3X1ikqM39rwAwJahHwEAFmOSmR3rkpza3TdW1a5Jbqiqv+7uW+cY+4ru/sH6he6+PcnBSVJV2yW5N8kXJj1uVe2R5Kfd/eCMdc/t7jtnnfecJH+R5D/OXDk+5xlJfivJ2iTXV9XFSVYkWT0e9tgE3wMAYLr0IwDAxBac2dHd93X3jeOvH0yyJsmeG3GuVya5q7u/s4jjHp7ki1W1Q5JU1YlJPjJHjdck+dEc53xhkju7+9vd/UiSzyZ5Q0aNxorxmDm/B1X1+qo664EHHljcpwQAnnL6Ef0IACzGoh5QWlV7JzkkyXVzbO4kV1bVDVV10hzb35zkvMUct7s/l+SKJOdX1VuSnJDkmEWUvGeSe2Ysrx2vuyjJm6rqzCSXzLVjd1/S3SctXz7hU8YBgC1CPwIALGTiB5RW1S5JPp/kPd39kzmGvLS77x1P9fzrqrptfIUjVfX0jO5Ffe9ij9vdH6iqzyY5M8lzuvuhSWueT3c/nOT4TT0OALBl6UcA2FrtfdqlE427+/TXbuZKtg0Tzeyoqu0zagDO7e6L5hrT3feO/74/o/tgXzhj82uS3Njd/2Oxx62qlyV5/viY75uk3hnuTbLXjOUV43UAwFZGPwIATGqSt7FUkrOTrOnuD80zZufxQ73Wv0btVUlmPon8uMyaMjrhcQ9JclZG97Uen2T3qnr/QjXPcH2Sfapq5fhqzpuTXLyI/QGAJUA/AgAsxiQzO16S5K1JjpjxyrajkqSqLquqZyV5RpJrq+qbSb6W5NLuvnw8ZueMnj4++0rJvMedYackx3b3Xd39eJK3JfnO7AKr6rwkX02yb1Wtrap3JEl3r0vyzozus12T5ILuvmWCzwwALC36EQBgYgs+s6O7r01S82yb2QwcNM+Yh5PsvpjjzhjzlVnLjyb52BzjjtvAMS5LctmGzgMALG36EQBgMRb1NhYAAACApU7YAQAAAAyKsAMAAAAYFGEHAAAAMCjCDgAAAGBQhB0AAADAoAg7AAAAgEERdgAAAACDIuwAAAAABkXYAQAAAAyKsAMAAAAYFGEHAAAAMCjCDgAAAGBQhB0AAADAoAg7AAAAgEERdgAAAACDsmzaBQAAwJr99p9s4G+esXkLAWAQzOwAAAAABkXYAQAAAAyKsAMAAAAYFGEHAAAAMCjCDgAAAGBQhB0AAADAoAg7AAAAgEERdgAAAACDIuwAAAAABkXYAQAAAAyKsAMAAAAYFGEHAAAAMCjCDgAAAGBQhB0AAADAoCybdgEAAMDSd+CnDlxwzOq3r94ClQAszMwOAAAAYFCEHQAAAMCgCDsAAACAQRF2AAAAAIMi7AAAAAAGRdgBAAAADIqwAwAAABgUYQcAAAAwKMIOAAAAYFCEHQAAAMCgCDsAAACAQRF2AAAAAIOybNoFTEtV7Zzko0keSXJ1d5875ZIAgG2IXgQANp9NmtlRVXtV1VVVdWtV3VJV755n3N1Vtbqqbqqqr4/X7TteXv/nJ1X1nk2o5RNVdX9V3Txr/ZFVdXtV3VlVp83Y9MYkF3b3iUmO3tjzAgDToxcBAOayqbexrEtyancfkORFSU6pqgPmGfuK7j64uw9Lku6+fbx8cJIXJPm7JF+YvVNV7VFVu85a99w5jn9OkiNnjdsuyRlJXpPkgCTHzahvRZJ7xl8/tuAnBQCWIr0IAPAkmxR2dPd93X3j+OsHk6xJsudGHOqVSe7q7u/Mse3wJF+sqh2SpKpOTPKROWq5JsmPZq1+YZI7u/vb3f1Iks8mecN429qMmozEs0sAYKukFwEA5vKU/WKtqr2THJLkujk2d5Irq+qGqjppju1vTnLeXMft7s8luSLJ+VX1liQnJDlmwrL2zBNXTJJRU7G+AbooyZuq6swkl8zesapeX1VnPfDAAxOeCgCYpqH1Iol+BAA21lPygNKq2iXJ55O8p7t/MseQl3b3vVW1R5K/rqrbxlc/UlVPz+g+1ffOd/zu/kBVfTbJmUme090PbWrN3f1wkuM3sP2SJJccdthhJ27quQCAzWuIvch4jH4EADbCJs/sqKrtM2ouzu3ui+Ya0933jv++P6N7YV84Y/NrktzY3f9jA+d4WZLnj/d93yLKuzfJXjOWV4zXAQADoRcBAGbb1LexVJKzk6zp7g/NM2bn9Q/1Gr9i7VVJZj6l/LjMM210vM8hSc7K6P7W45PsXlXvn7DE65PsU1Urx1dt3pzk4gn3BQCWOL0IADCXTZ3Z8ZIkb01yxIzXth2VJFV1WVU9K8kzklxbVd9M8rUkl3b35eMxOyf5rYzuWZ3PTkmO7e67uvvxJG9L8qSHh1XVeUm+mmTfqlpbVe/o7nVJ3pnRfbZrklzQ3bds4mcGAJYOvQgA8CSb9MyO7r42Sc2z7agZiwfNM+bhJLsvcI6vzFp+NMnH5hh33Dz7X5bksg2dAwDYOulFAIC5eM0ZAAAAMCjCDgAAAGBQhB0AAADAoAg7AAAAgEERdgAAAACDIuwAAAAABkXYAQAAAAyKsAMAAAAYFGEHAAAAMCjCDgAAAGBQhB0AAADAoAg7AAAAgEERdgAAAACDIuwAAAAABkXYAQAAAAyKsAMAAAAYFGEHAAAAMCjCDgAAAGBQhB0AAADAoAg7AAAAgEERdgAAAACDIuwAAAAABmXZtAsAAABgyzjwUwdONG7121dv5kpg8zKzAwAAABgUYQcAAAAwKMIOAAAAYFCEHQAAAMCgCDsAAACAQfE2FgAAANiKeKvOwszsAAAAAAbFzA4AAOApsWa//Scat/9tazZzJcC2TtgBAAAsOR9c9bqJxp16/pc2cyXA1kjYAQAAbFFnnPzlaZcADJxndgAAAACDIuwAAAAABsVtLAAAAEPwh8sXHrPy1zd/HbAEmNkBAAAADIqZHQAAsC2bZDZAYkYAsFUxswMAAAAYFDM7tlaTJvB/+MDmrQMAAICnjmevPCXM7AAAAAAGRdgBAAAADIqwAwAAABgUYQcAAAAwKMIOAAAAYFC2ybexVNXOST6a5JEkV3f3uVMuCQDYxuhHAGDzWXBmR1XtVVVXVdWtVXVLVb17nnF3V9Xqqrqpqr4+Y/1uVXVhVd1WVWuq6sUztv2v42PeXFXnVdWOG/MhquoTVXV/Vd08x7Yjq+r2qrqzqk4br35jkgu7+8QkR2/MOQGALUc/AgAsxiQzO9YlObW7b6yqXZPcUFV/3d23zjH2Fd39g1nrPpzk8u7+nap6epKdkqSq9kzyL5Mc0N0/raoLkrw5yTnrd6yqPZL8tLsfnLHuud1956xznJPkL5L8x5krq2q7JGck+a0ka5NcX1UXJ1mRZPV42GMTfA8AgOnSjwDAIq3Zb/8Fx+x/25otUMmWt+DMju6+r7tvHH/9YJI1Sfac5OBVtTzJy5OcPd7/ke7+8Ywhy5L8SlUty6jp+O+zDnF4ki9W1Q7j452Y5CNz1HhNkh/NUcILk9zZ3d/u7keSfDbJGzJqNFaMx8z5Paiq11fVWQ888MAkHxUA2Iz0I/oRAFiMRT2gtKr2TnJIkuvm2NxJrqyqG6rqpPG6lUm+n+STVfWNqvr4+P7UdPe9Sf59ku8muS/JA9195S8dsPtzSa5Icn5VvSXJCUmOWUTJeya5Z8by2vG6i5K8qarOTHLJXDt29yXdfdLy5csXcToAYHPTjwAAC5k47KiqXZJ8Psl7uvsncwx5aXcfmuQ1SU6pqpdndKXk0CRndvchSR5Octr4eL+a0VWNlUmelWTnqvq92Qft7g8k+VmSM5Mc3d0PLeLzzam7H+7u47v7DzwMDAC2HvoRAGASE72Npaq2z6ixOLe7L5przPjKSLr7/qr6QkZTNj+TZG13r7/ycmHGzUWSf57kv3X398fnuCjJPxvvM/PcL0vy/CRfSPK+JO+c+NMl9ybZa8byivG6bcaBnzpwwTGr3756wTEAMG36EQBgUpO8jaUyusd1TXd/aJ4xO48fFrb+NWqvSnJzd38vyT1Vte946CuTrH+Q2HeTvKiqdhqf45UZ3X8787iHJDkroysuxyfZvarev4jPd32Sfapq5fhhZG9OcvEi9gcAlgD9CACwGJPcxvKSJG9NcsT4NW43VdVRSVJVl1XVs5I8I8m1VfXNJF9Lcml3Xz7e/11Jzq2qbyU5OMmfJMn46sqFSW7M6EnkT8uokZhppyTHdvdd3f14krcl+c7sAqvqvCRfTbJvVa2tqneMz7EuoysvV2TUuFzQ3bdM8o0BAJYU/QgAMLEFb2Pp7muT1DzbjpqxeNA8Y25Kctg8296X0VTQ+c79lVnLjyb52BzjjtvAMS5Lctl82wGApU8/AgAsxkTP7AAAAAC2XR9c9boFx5x6/pe2QCWTEXYAAACw2Wxt/5PMMEz86lkAAACArYGwAwAAABgUt7HAZrL3aZdONO7u01+7mSsBAADYtgg7AAAA+CVr9tt/wTH737ZmC1QCG8dtLAAAAMCgCDsAAACAQRF2AAAAAIMi7AAAAAAGRdgBAAAADIqwAwAAABgUYQcAAAAwKMIOAAAAYFCEHQAAAMCgCDsAAACAQRF2AAAAAIMi7AAAAAAGRdgBAAAADIqwAwAAABiUZdMuAAAAgK3PGSd/edolwLzM7AAAAAAGRdgBAAAADIqwAwAAABgUYQcAAAAwKMIOAAAAYFC8jQUAAAC2UUN9q46ZHQAAAMCgmNnBxCZJ/E75yyO2QCUAAAAwPzM7AAAAgEERdgAAAACDIuwAAAAABkXYAQAAAAyKsAMAAAAYFGEHAAAAMCjCDgAAAGBQhB0AAADAoAg7AAAAgEERdgAAAACDIuwAAAAABmXZtAtgWD646nUTjTv1/C9t5koAAADYVpnZAQAAAAyKsAMAAAAYFGEHAAAAMCie2QGwAXufdulE4+4+/bWbuRIAAGBS22TYUVU7J/lokkeSXN3d5065JAA2I6EVS5F+BAA2nwVvY6mqvarqqqq6tapuqap3zzPu7qpaXVU3VdXXZ6zfraourKrbqmpNVb14km2LUVWfqKr7q+rmObYdWVW3V9WdVXXaePUbk1zY3ScmOXpjzgkAbDn6EQBgMSZ5Zse6JKd29wFJXpTklKo6YJ6xr+jug7v7sBnrPpzk8u7eL8lBSdZMuC1VtUdV7Tpr3XPnOO85SY6cvbKqtktyRpLXJDkgyXHj2lckuWc87LF5PgsAsHToRwCAiS0YdnT3fd194/jrBzNqAPac5OBVtTzJy5OcPd7/ke7+8ULbZjg8yReraofxPicm+cgcNV6T5EdzlPDCJHd297e7+5Ekn03yhiRrM2owknm+B1X1+qo664EHHpjkowIAm5F+RD8CAIuxqLexVNXeSQ5Jct0cmzvJlVV1Q1WdNF63Msn3k3yyqr5RVR8f35+60LbRAbs/l+SKJOdX1VuSnJDkmEWUvGeeuGKSjJqKPZNclORNVXVmkkvm2rG7L+nuk5YvX76I0wEAm5t+BABYyMRhR1XtkuTzSd7T3T+ZY8hLu/vQjKZonlJVL8/oAaiHJjmzuw9J8nCS9fepbmjbL3T3B5L8LMmZSY7u7ocmrXk+3f1wdx/f3X/gYWAAsPXQjwAAk5go7Kiq7TNqLM7t7ovmGtPd947/vj/JFzKasrk2ydruXn/l5cKMGoossG3muV+W5PnjY75vknpnuDfJXjOWV4zXAQBbGf0IADCpSd7GUhndx7qmuz80z5id1z+4azz181VJbu7u7yW5p6r2HQ99ZZJbk2RD22Yc95AkZ2V0X+vxSXavqvcv4vNdn2SfqlpZVU9P8uYkFy9ifwBgCdCPAACLMcnMjpckeWuSI8avcbupqo5Kkqq6rKqeleQZSa6tqm8m+VqSS7v78vH+70pyblV9K8nBSf5kxrE3tC1JdkpybHff1d2PJ3lbku/MLrCqzkvy1ST7VtXaqnpHknT3uiTvzOg+2zVJLujuWyb4zADA0qIfAQAmtmyhAd19bZKaZ9tRMxYPmmfMTUkOW+y28favzFp+NMnH5hh33AaOcVmSy+bbDgAsffoRAGAxFvU2FgAAAIClTtgBAAAADIqwAwAAABgUYQcAAAAwKMIOAAAAYFCEHQAAAMCgCDsAAACAQVk27QKYvjX77T/ZwN88Y/MWAgAAAE8BMzsAAACAQRF2AAAAAIMi7AAAAAAGRdgBAAAADIqwAwAAABgUYQcAAAAwKMIOAAAAYFCEHQAAAMCgCDsAAACAQRF2AAAAAIMi7AAAAAAGRdgBAAAADIqwAwAAABgUYQcAAAAwKMIOAAAAYFCEHQAAAMCgCDsAAACAQRF2AAAAAIMi7AAAAAAGRdgBAAAADIqwAwAAABgUYQcAAAAwKMumXQCwsAM/deBE41a/ffVmrgQAAGDpM7MDAAAAGBRhBwAAADAowg4AAABgUIQdAAAAwKAIOwAAAIBBEXYAAAAAgyLsAAAAAAZF2AEAAAAMirADAAAAGBRhBwAAADAowg4AAABgUIQdAAAAwKAIOwAAAIBBWTbtAqahqnZO8tEkjyS5urvPnXJJAMA2Rj8CAJvPgjM7qmqvqrqqqm6tqluq6t3zjLu7qlZX1U1V9fUZ63erqgur6raqWlNVL56133ZV9Y2q+tLGfoiq+kRV3V9VN8+x7ciqur2q7qyq08ar35jkwu4+McnRG3teAGDL0I8AAIsxyW0s65Kc2t0HJHlRklOq6oB5xr6iuw/u7sNmrPtwksu7e78kByVZM2ufd8+xLklSVXtU1a6z1j13jqHnJDlyjv23S3JGktckOSDJcePaVyS5ZzzssXk+CwCwdOhHAICJLRh2dPd93X3j+OsHM2oE9pzk4FW1PMnLk5w93v+R7v7xjO0rkrw2ycfnOcThSb5YVTuMx5+Y5CNz1HhNkh/Nsf8Lk9zZ3d/u7keSfDbJG5KszajBSDy3BACWPP0IALAYi/rFWlV7JzkkyXVzbO4kV1bVDVV10njdyiTfT/LJ8dTQj4/vT13v/0ryr5M8Ptf5uvtzSa5Icn5VvSXJCUmOWUTJe+aJKybJqKnYM8lFSd5UVWcmuWSuHavq9VV11gMPPLCI0wEAm5t+BABYyMRhR1XtkuTzSd7T3T+ZY8hLu/vQjKZonlJVL8/oAaiHJjmzuw9J8nCS08bHe12S+7v7hg2dt7s/kORnSc5McnR3PzRpzRs45sPdfXx3/8F8DwPr7ku6+6Tly5dv6ukAgKeIfgQAmMREYUdVbZ9RY3Fud18015juvnf89/1JvpDRlM21SdZ29/orLxdm1GwkyUuSHF1Vd2c0nfOIqvrMHOd+WZLnj4/5vsk+1i/cm2SvGcsrxusAgK2MfgQAmNQkb2OpjO5xXdPdH5pnzM7rH9w1nhb6qiQ3d/f3ktxTVfuOh74yya1J0t3v7e4V3b13kjcn+XJ3/96s4x6S5KyM7ms9PsnuVfX+RXy+65PsU1Urq+rp4/NcvIj9AYAlQD8CACzGJDM7XpLkrRld6bhp/OeoJKmqy6rqWUmekeTaqvpmkq8lubS7Lx/v/64k51bVt5IcnORPFlHfTkmO7e67uvvxJG9L8p3Zg6rqvCRfTbJvVa2tqnckSXevS/LOjO6zXZPkgu6+ZRHnBwCWBv0IADCxZQsN6O5rk9Q8246asXjQPGNuSnLYXNtmjLk6ydVzrP/KrOVHk3xsjnHHbeDYlyW5bEPnBwCWNv0IALAYXnMGAAAADIqwAwAAABgUYQcAAAAwKMIOAAAAYFCEHQAAAMCgCDsAAACAQRF2AAAAAIMi7AAAAAAGRdgBAAAADIqwAwAAABgUYQcAAAAwKMIOAAAAYFCWTbsAYMs64+QvTzTulL88YjNXAgAAsHmY2QEAAAAMirADAAAAGBRhBwAAADAowg4AAABgUIQdAAAAwKAIOwAAAIBBEXYAAAAAgyLsAAAAAAZF2AEAAAAMirADAAAAGBRhBwAAADAoy6ZdAGzz/nD5wmNW/vrmrwMAAGAgzOwAAAAABsXMDoAt5MBPHTjRuNVvX72ZK2Fek8y0+sMHNn8dAABsEjM7AAAAgEERdgAAAACDIuwAAAAABsUzOwCeCt6qAwAAS4aZHQAAAMCgmNkBzOmDq1634JhTz//SFqgEAABgcYQdALAZrNlv/wXH7H/bmi1QCQDAtsdtLAAAAMCgCDsAAACAQRF2AAAAAIPimR0AsAgHfurAicZdsJnrAABgfsIOGJBJHoiY3zxj8xcCAAAwRW5jAQAAAAZF2AEAAAAMirADAAAAGBRhBwAAADAowg4AAABgUIQdAAAAwKBsk6+eraqdk3w0ySNJru7uc6dcEsAvTPIK4f1vWzPRsT646nULjjn1/C9NdCzgqaUfAYDNZ8GZHVW1V1VdVVW3VtUtVfXuecbdXVWrq+qmqvr6jPW7VdWFVXVbVa2pqhcv5riTqKpPVNX9VXXzHNuOrKrbq+rOqjptvPqNSS7s7hOTHL2x5wUAtgz9CACwGJPcxrIuyandfUCSFyU5paoOmGfsK7raTsKyAAARG0lEQVT74O4+bMa6Dye5vLv3S3JQkvWXIxc8blXtUVW7zlr33DnOe06SI2evrKrtkpyR5DVJDkhy3PgcK5LcMx722DyfBQBYOvQjAMDEFgw7uvu+7r5x/PWDGTUHe05y8KpanuTlSc4e7/9Id/94Ecc9PMkXq2qH8fFOTPKROWq8JsmP5ijhhUnu7O5vd/cjST6b5A1J1mbUYCSeWwIAS55+BABYjEX9Yq2qvZMckuS6OTZ3kiur6oaqOmm8bmWS7yf5ZFV9o6o+Pr4/daLjdvfnklyR5PyqekuSE5Ics4iS98wTV0ySUVOxZ5KLkrypqs5McslcO1bV66vqrAceeGARpwMANjf9CACwkInDjqraJcnnk7ynu38yx5CXdvehGU3RPKWqXp7RA1APTXJmdx+S5OEkp83caaHjdvcHkvwsyZlJju7uhyateT7d/XB3H9/dfzDfw8C6+5LuPmn58uWbejoA4CmiHwEAJjFR2FFV22fUAJzb3RfNNaa77x3/fX+SL2Q0ZXNtkrXdvf4KyYUZNRsTH7eqXpbk+eNjvm+Seme4N8leM5ZXjNcBAFsZ/QgAMKlJ3sZSGd3juqa7PzTPmJ3XP7hrPC30VUlu7u7vJbmnqvYdD31lklsXcdxDkpyV0X2txyfZvarev4jPd32SfapqZVU9Pcmbk1y8iP0BgCVAPwIALMYkMztekuStSY4Yv8btpqo6Kkmq6rKqelaSZyS5tqq+meRrSS7t7svH+78ryblV9a0kByf5k4WOO8NOSY7t7ru6+/Ekb0vyndkFVtV5Sb6aZN+qWltV70iS7l6X5J0Z3We7JskF3X3LJN8YAGBJ0Y8AABNbttCA7r42Sc2zbWYzcNA8Y25Kctgc6+c97owxX5m1/GiSj80x7rgNHOOyJJdt6DwAwNKmHwEAFmPBsAOApeeMk7887RIAAGDJ8k53AAAAYFCEHQAAAMCgCDsAAACAQRF2AAAAAIMi7AAAAAAGRdgBAAAADIqwAwAAABgUYQcAAAAwKMIOAAAAYFCEHQAAAMCgCDsAAACAQRF2AAAAAIMi7AAAAAAGZdm0CwAANuyDq1634JhTz//SFqgEAGDrYGYHAAAAMCjCDgAAAGBQhB0AAADAoAg7AAAAgEERdgAAAACDIuwAAAAABkXYAQAAAAyKsAMAAAAYFGEHAAAAMCjCDgAAAGBQhB0AAADAoAg7AAAAgEFZNu0CAGBbdcbJX552CQAAg2RmBwAAADAowg4AAABgUIQdAAAAwKAIOwAAAIBBEXYAAAAAgyLsAAAAAAZF2AEAAAAMirADAAAAGBRhBwAAADAowg4AAABgUIQdAAAAwKAIOwAAAIBBEXYAAAAAgyLsAAAAAAZF2AEAAAAMirADAAAAGBRhBwAAADAoy6ZdwLRU1c5JPprkkSRXd/e5Uy4JANiG6EUAYPPZpJkdVbVXVV1VVbdW1S1V9e55xt1dVaur6qaq+vqM9btV1YVVdVtVramqF29CLZ+oqvur6uZZ64+sqtur6s6qOm3GpjcmubC7T0xy9MaeFwCYHr0IADCXTb2NZV2SU7v7gCQvSnJKVR0wz9hXdPfB3X3YjHUfTnJ5d++X5KAka2bvVFV7VNWus9Y9d47jn5PkyFnjtktyRpLXJDkgyXEz6luR5J7x14/N/xEBgCVMLwIAPMkmhR3dfV933zj++sGMGoQ9J9m3qpYneXmSs8f7P9LdP55j6OFJvlhVO4z3OzHJR+ao5ZokP5q1+oVJ7uzub3f3I0k+m+QN421rM2oykjm+D1X1+qo664EHHpjk4wAAUzDkXmR8Lv0IAGyEp+wBpVW1d5JDklw3x+ZOcmVV3VBVJ43XrUzy/SSfrKpvVNXHx/eu/vKO3Z9LckWS86vqLUlOSHLMhGXtmSeumCSjpmJ9A3RRkjdV1ZlJLpnjvJd090nLly+f8FQAwDQNrRcZn1s/AgAb4Sl5QGlV7ZLk80ne090/mWPIS7v73qraI8lfV9VtSf4uyaFJ3tXd11XVh5OcluTfzN65uz9QVZ9NcmaS53T3Q5tac3c/nOT4TT0OADB9ehEAYKZNntlRVdtn1Fyc290XzTWmu+8d/31/ki9kNKVzbZK13b3+6suFGTUcc53jZUmeP973fYso794ke81YXjFeBwAMhF4EAJhtU9/GUhnd57qmuz80z5id1z/Uazw19FVJbu7u7yW5p6r2HQ99ZZJb59j/kCRnZXR/6/FJdq+q909Y4vVJ9qmqlVX19CRvTnLxxB8QAFjS9CIAwFw2dWbHS5K8NckR41e53VRVRyVJVV1WVc9K8owk11bVN5N8Lcml3X35eP93JTm3qr6V5OAkfzLHOXZKcmx339Xdjyd5W5LvzB5UVecl+WqSfatqbVW9o7vXJXlnRvfZrklyQXffsomfGQBYOvQiAMCTbNIzO7r72iQ1z7ajZiweNM+Ym5IcNte2GWO+Mmv50SQfm2PccfPsf1mSyzZ0DgBg66QXAQDmUt097RqWtKr6fua4ejMw/yDJD6ZdBE8JP8th8HMchm3h5/js7v6H0y5iW6AfYSvi5zgcfpbDsC38HOfsR4QdpKq+3t0bvKrF1sHPchj8HIfBzxEWx7+ZYfBzHA4/y2HYln+Om/w2FgAAAIClRNgBAAAADIqwg2T0Oj2Gwc9yGPwch8HPERbHv5lh8HMcDj/LYdhmf46e2QEAAAAMipkdAAAAwKAIOwAAAIBBEXYAAAAAgyLsIFW1x7RrABiSqvp7VfWCqvrVadcCWwv9CMBTRy8i7NjmVNXfn/Vn9yRfq6pfraq/P+36mMz4P15/WlWfrqrfnbXto9Oqi41TVb9eVbuNv967qn6nqp4/7bqYXFV9pqr+wfjrVye5Ocm/S3JTVR0z1eJgCdKPDIN+ZFj0I1s3vciTeRvLNqaqHk/ynVmrVyRZm6S7+x9t+apYrKr6fJL/muS/JDkhyaNJfre7f15VN3b3oVMtkIlV1WlJ/pckP0/y75P8b0m+kuRFSc7u7g9NsTwmVFWru/vA8dd/m9G/x7vHTcf/290HTbdCWFr0I8OgHxkO/cjWTy/yZMumXQBb3L9K8ltJ/lV3r06Sqvpv3b1yumWxSM/p7jeNv/5iVf0fSb5cVUdPsyg2yluTHJBkpyR3J/lH3f39qto5yXVJNBdbh6dV1d/r7p8keTzJd5Oku39QVX7XwpPpR4ZBPzIc+pGtn15klm3yQ2/LuvuDVXV+kj+vqnuSvC+J6T1bnx2q6mnd/XiSdPcfV9W9Sa5Jsst0S2ORHuvun1bVI0l+muSHSdLdD1fVdCtjMf5tkquq6oyMroR9rqouTvKKJJdPtTJYgvQjg6EfGQ79yNZPLzKL21i2YePU/X9Psnd3P3Pa9TC5qvpAkiu7+z/NWn9kko909z7TqYzFqqpzkjw9yc5J/i7Juox+IR2RZNfuPnZ61bEYVbVPkv85yfMyupiwNskXu/uKqRYGS5x+ZOulHxkO/cgw6EV+mbBjG1RV+yXZM6MpaY9lNAXx5qo6sru3ydRva1NV/zLJF7r7nmnXwqYZTys8JqMrmhcm+adJjsto6uEZ3f3wFMsD2Gz0I1s//chw6EcYImHHNmb8S+mUJGuSHJzk3d39V+NtHiS1laiqB5I8nOSuJOcl+Vx3f3+6VfFUqardu/uH066DyVXVTknemVGT+JEkq5K8KcltSf7P7n5oiuXBkqMfGQb9yLDpR7YuepEn8+rZbc+JSV7Q3f9Tkt9M8m+q6t3jbW7I23p8O6On1v9RkhckubWqLq+qt1fVrtMtjcWoqtNnvCbssKr6dpL/UlXfqarDp1wekzsnyTOSrExyaZLfSPJnGf139czplQVLln5kGPQjA6EfGYRzohf5JWZ2bGOq6pbu/sczlnfJaKrarUmO6O6Dp1YcE5t91auqtk/ymoymG/7z7v6HUyuORZn1mrCrkvzr7r6+qp6X5P/p7sOmWyGTqKqbuvvgGj3F7b4kv9bdPV7+Znf/kymXCEuKfmQY9CPDoR/Z+ulFnszMjm3P/6iqXzQQ4+lMr0vyD5IcOLWqWKxfuurV3Y9298XdfVySZ0+pJjbOshmvA/uV7r4+Sbr7jiQ7TK8sNkaPriBcNv57/bKrCvBk+pFh0I8Mh35kIPQiTxB2bHveluR7M1d097rufluSl0+nJDbCqvk2dPffbclC2GQfTXJZVR2R5PKq+nBVHV5V/zbJTVOujcl9fXxlOt19wvqVVfWcJA9OrSpYuvQjw6AfGQ79yNZPLzKL21gApqyqfjPJH+SJ14Tdk+SLST7R3eumWBqLUFUvzOgCyvVVdUCSI5PcnhlXVwBgqdKPbP30Ir9M2AGwRFXV8d39yWnXwcKq6n0Z3ae+LMlfZ/TKvquS/FaSK7r7j6dYHgBsNP3I1kEv8mTCDoAlqqq+292/Pu06WFhVrc7o9Zk7ZDQ1f0V3/6SqfiXJddviQ8EAGAb9yNZBL/JkyxYeAsDmUlXfmm9TRq8PY+uwrrsfS/J3VXVXd/8kSbr7p1X1+JRrA4AN0o8Mgl5kFmEHwHQ9I8mrk/x/s9ZXkr/d8uWwkR6pqp3GD+R7wfqVVbU8yTbZYACwVdGPbP30IrMIOwCm60tJdunuJz3pvKqu3vLlsJFe3t0/T5LuntlQbJ/k7dMpCQAmph/Z+ulFZvHMDgAAAGBQnjbtAgAAAACeSsIOAAAAYFCEHcBEqmr3qrpp/Od7VXXvjOWn/MFVVfWNqjp4/PWyqnqoqn5vxvYbqurQqjq6qk5b5LHPqarfeaprBgA2L/0IMCkPKAUm0t0/zOjd3amqP0zyUHf/+814yq8k+WdJbkpyUJI7xsufqaqdkzwnyTe7+8YkF2/GOgCAJUI/AkzKzA5gk1XVQ+O/f7Oq/nNV/VVVfbuqTq+qt1TV16pqdVU9ZzzuH1bV56vq+vGfl8xx2L/NqJnI+O+/zLi5SfLCJDd092NV9ftV9Rfj455TVf93Vf3t+Py/M15fVfUXVXV7Vf2nJHvMqP2V46s2q6vqE1W1Q1X9RlVdNN7+hqr6aVU9vap2rKpvP/XfQQBgU+lHgJmEHcBT7aAkJyfZP8lbkzyvu1+Y5ONJ3jUe8+Ekf97dv5HkTeNts62/kpLx39ck+XlV7Tpenm+q6q8leWmS1yU5fbzut5Psm+SAJG9bf9yq2jHJOUlWdfeBGc12+4Mk38gTjczLktyc5DeS/NMk1032bQAApkg/Ats4YQfwVLu+u+8bv+f7riRXjtevTrL3+Ot/nuQvquqmjKZ8/r2q2mXmQbr7O0meXlXPTLJfktuTXJ/RL/h/llHzMZcvdvfj3X1rkmeM1708yXnd/Vh3//ckXx6v3zfJf+vuO8bLn8roHeXrktxVVftndNXmQ+NjvCzJ3yz6OwIAbGn6EdjGeWYH8FT7+YyvH5+x/Hie+G/O05K8qLt/tsCx/jbJMUnu6+6uqv+S5CUZ/cL/6gTnr8UUPss1SV6T5NEk/ymjKy7bJflXm3BMAGDL0I/ANs7MDmAarswTU0iz/innc/jbJO/JE43EVzOa9vm97n5gEee7Jsmqqtquqn4tySvG629PsndVPXe8/NYk/3n89d+sP3d3fz/J7hldebl5EecFAJYu/QgMmLADmIZ/meSwqvpWVd2a0T21c/lKkn+UcXPR3fdldDVjsa+W+0KS/5rk1iT/ccbxfpbk+CSfq6rVGV3t+cvxPtdlNO30mvHyt5Ks7u5e5LkBgKVJPwIDVv6dAP9/u3ZAAgAAACDo/+t+hM6IAAAATpwdAAAAwIrYAQAAAKyIHQAAAMCK2AEAAACsiB0AAADAitgBAAAArIgdAAAAwIrYAQAAAKwEef/jUhq4uZIAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1296x648 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"stacking_results = stacking_results.replace(\n",
" {'StackingMethod': {'clf': 'Classifier', 'soft_voting': 'Soft voting'}}\n",
" )\n",
"\n",
"fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(18,9))\n",
"\n",
"alpha = 0.001\n",
"y_range = [\n",
" stacking_results['LogLoss'].min() - alpha,\n",
" stacking_results['LogLoss'].max() + alpha\n",
" ]\n",
"\n",
"for i, technique in enumerate(stacking_results['EncodingOfCyclicalAttributes'].unique()):\n",
" mask = stacking_results['EncodingOfCyclicalAttributes'] == technique\n",
" technique_results = stacking_results.loc[mask].drop(columns='EncodingOfCyclicalAttributes')\n",
" \n",
" technique_results = technique_results.pivot_table(\n",
" index='TimeWindow', columns=['NumberOfEstimators', 'StackingMethod'], values='LogLoss'\n",
" )\n",
" \n",
" technique_results.plot(kind='bar', ax=ax[i])\n",
" \n",
" technique = 'one-hot' if technique == 'onehot' else technique\n",
" ax[i].set_title('{} Encoding Technique'.format(technique.title()))\n",
" ax[i].set_xlabel('Time Window')\n",
"\n",
" ax[i].set_ylim(*y_range)\n",
" ax[i].set_yscale(\"log\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Overall, the above figures support the finding that the technique that best handles cyclical attributes is one-hot encoding. Likewise, findings with respect to the time window and the stacking technique might also be generalized."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Third of all, let's get insights into how much the meta-model strengthens the discriminative power of the base models. To accomplish this, a comparison between the two best meta-models and the set of seven base models they are built on is conducted."
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<center><table>\n",
"<thead>\n",
"<tr><th style=\"text-align: right;\"> Base Model</th><th style=\"text-align: right;\"> Log Loss</th></tr>\n",
"</thead>\n",
"<tbody>\n",
"<tr><td style=\"text-align: right;\"> 1</td><td style=\"text-align: right;\"> 2.5726 </td></tr>\n",
"<tr><td style=\"text-align: right;\"> 2</td><td style=\"text-align: right;\"> 2.56885</td></tr>\n",
"<tr><td style=\"text-align: right;\"> 3</td><td style=\"text-align: right;\"> 2.57124</td></tr>\n",
"<tr><td style=\"text-align: right;\"> 4</td><td style=\"text-align: right;\"> 2.57141</td></tr>\n",
"<tr><td style=\"text-align: right;\"> 5</td><td style=\"text-align: right;\"> 2.5707 </td></tr>\n",
"<tr><td style=\"text-align: right;\"> 6</td><td style=\"text-align: right;\"> 2.56979</td></tr>\n",
"<tr><td style=\"text-align: right;\"> 7</td><td style=\"text-align: right;\"> 2.57093</td></tr>\n",
"</tbody>\n",
"</table></center>"
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"best_setting = stacking_results.iloc[0]\n",
"\n",
"predictions_path = os.path.join(\n",
" DATA_PATH,\n",
" '{}H'.format(best_setting['TimeWindow']),\n",
" 'prediction',\n",
" 'validation',\n",
" 'n_estimators={}'.format(best_setting['NumberOfEstimators']),\n",
" 'cyclical_attr={}'.format(best_setting['EncodingOfCyclicalAttributes'])\n",
" )\n",
"\n",
"base_models_performance = pd.read_csv(os.path.join(predictions_path, 'learners-result.csv'))\n",
"\n",
"header = ['Base Model', 'Log Loss']\n",
"\n",
"table = [\n",
" [str(i+1), '{:.5f}'.format(row['LogLoss'])]\n",
" for i, (idx, row) in enumerate(base_models_performance.iterrows())\n",
" ]\n",
"\n",
"table = ('<center>'\n",
" + tabulate.tabulate(table, header, tablefmt='html')\n",
" + '</center>')\n",
"\n",
"display(HTML(table))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With this in mind, let's quantify how much the two best meta-models strengthen the discriminative power w.r.t. the best, the worst, and the average base model."
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<center><table>\n",
"<thead>\n",
"<tr><th>Stacking Technique </th><th>Best Base Model </th><th>Worst Base Model </th><th>Average Base Model </th></tr>\n",
"</thead>\n",
"<tbody>\n",
"<tr><td>Classifier </td><td>0.315 % </td><td>0.460 % </td><td>0.390 % </td></tr>\n",
"<tr><td>Soft voting </td><td>0.295 % </td><td>0.441 % </td><td>0.371 % </td></tr>\n",
"</tbody>\n",
"</table></center>"
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"header = [\n",
" 'Stacking Technique',\n",
" 'Best Base Model',\n",
" 'Worst Base Model',\n",
" 'Average Base Model'\n",
" ]\n",
"\n",
"table = []\n",
"\n",
"comparison_values = [\n",
" base_models_performance['LogLoss'].min(),\n",
" base_models_performance['LogLoss'].max(),\n",
" base_models_performance['LogLoss'].mean()\n",
" ]\n",
"\n",
"for i in range(2):\n",
" stack_technique = stacking_results.iloc[i]['StackingMethod']\n",
" \n",
" score = stacking_results.iloc[i]['LogLoss']\n",
" \n",
" comparison = ['{:.3f} %'.format((np.abs(score-val)/val)*100) for val in comparison_values]\n",
" comparison.insert(0, stack_technique)\n",
" \n",
" table.append(comparison)\n",
"\n",
"table = ('<center>'\n",
" + tabulate.tabulate(table, header, tablefmt='html')\n",
" + '</center>')\n",
"\n",
"display(HTML(table))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To conclude, let's make the following final remarks:\n",
"\n",
"1. The contributions of the two best meta-models to the discriminative power is almost negligible; all of them are below 0.5%. Even more, deploying any of these meta-models and its respective base models to a production environment might be unfeasible.\n",
"2. However, any gain, no matter how small, is worth in the context of a Kaggle competition.\n",
"3. The prediction from the two best meta-models on the test set will be used as late submissions to participate in this Kaggle competition.\n",
"4. Regarding the baseline, i.e., 2.5947 multi-class logarithmic loss, the best meta-model has outperformed it. The decrease (or gain) is around <b>1.31%</b>."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 3.4.2 Late Submission"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As pointed out above, the final predictions to submit are those from the two best meta-models. These meta-models are built on the predictions from the same set of seven base models, but differ from each other by the stacking technique. Therefore, let's recall how the two best meta-models are built, as shown by the following table:\n",
"\n",
"| Rank | Time Window | Number of Estimators | Stacking Technique | Cyclical Attributes Encoding Technique |\n",
"|------|-------------|----------------------|--------------------|----------------------------------------|\n",
"| 1 | 336 | 7 | Classifier | One-hot |\n",
"| 2 | 336 | 7 | Soft voting | One-hot |"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [],
"source": [
"def train_and_persist_base_models(\n",
" train_fname,\n",
" validation_fname,\n",
" time_window,\n",
" n_estimators,\n",
" base_algorithm,\n",
" meta_algorithm,\n",
" models_path\n",
" ):\n",
" \"\"\"Fit the base models on the training dataset and persist them.\"\"\"\n",
" if not os.path.isdir(models_path):\n",
" os.makedirs(models_path)\n",
" \n",
" results_fname = os.path.join(models_path, 'base-learners-result.csv')\n",
" if not os.path.isfile(results_fname):\n",
" write_in_file(results_fname, ','.join(['Clf', 'LogLoss']))\n",
" \n",
" train = pd.read_csv(train_fname)\n",
" validation = pd.read_csv(validation_fname)\n",
" \n",
" cyclical_attr = [\n",
" 'Month', 'Day', 'DayOfWeek', 'Hour'\n",
" ]\n",
" train = encode_categorical_attributes(train, cyclical_attr)\n",
" validation = encode_categorical_attributes(validation, cyclical_attr)\n",
" \n",
" date_attr = [\n",
" 'Quarter', 'Triannual', 'Semester', 'Fortnight',\n",
" 'FourHourPeriod', 'SixHourPeriod', 'TwelveHourPeriod'\n",
" ]\n",
" train = encode_categorical_attributes(train, date_attr)\n",
" validation = encode_categorical_attributes(validation, date_attr)\n",
" \n",
" train = encode_categorical_attributes(train, ['PdDistrict'])\n",
" validation = encode_categorical_attributes(validation, ['PdDistrict'])\n",
" \n",
" train = train.replace({'Category': CRIME_CATEGORY})\n",
" validation = validation.replace({'Category': CRIME_CATEGORY})\n",
" \n",
" X = None\n",
" \n",
" for i, train_sample in enumerate(split_training_data(train, n_estimators)):\n",
" clf_fname = os.path.join(models_path, 'clf-{}.joblib'.format(i))\n",
" scaler_fname = os.path.join(models_path, 'clf-{}-scaler.joblib'.format(i))\n",
" \n",
" if os.path.isfile(clf_fname):\n",
" continue\n",
" \n",
" score, predictions, clf, scaler = build_model(\n",
" train_sample, validation,\n",
" 'all', 'Category',\n",
" time_window, 'cumulative',\n",
" base_algorithm['estimator'], None,\n",
" base_algorithm['predict_proba'],\n",
" return_pred=True,\n",
" return_clf=True\n",
" )\n",
" \n",
" X = copy.deepcopy(predictions) if X is None else np.hstack([X, predictions])\n",
" \n",
" write_in_file(results_fname, ','.join([str(i), '{:.5f}'.format(score)]), mode='a')\n",
" \n",
" joblib.dump(clf, clf_fname)\n",
" joblib.dump(scaler, scaler_fname)\n",
" \n",
" clf_ensemble_fname = os.path.join(models_path, 'clf-ensemble.joblib')\n",
" if os.path.isfile(clf_ensemble_fname):\n",
" return\n",
" \n",
" n_cols = len(CRIME_CATEGORY) * n_estimators\n",
" \n",
" assert n_cols == X.shape[1]\n",
" \n",
" clf_ensemble = meta_algorithm['estimator']\n",
" clf_ensemble.fit(X, validation['Category'].values.astype(int))\n",
" \n",
" joblib.dump(clf_ensemble, clf_ensemble_fname)"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [],
"source": [
"def first_level_predictions(\n",
" test_fname,\n",
" time_window,\n",
" n_estimators,\n",
" models_path,\n",
" predictions_path):\n",
" \"\"\"Make predictions from the base models on the test set.\"\"\"\n",
" prediction_header = ['' for j in range(len(CRIME_CATEGORY))]\n",
" for crime, idx in CRIME_CATEGORY.items():\n",
" prediction_header[idx] = crime\n",
" \n",
" if not os.path.isdir(predictions_path):\n",
" os.makedirs(predictions_path)\n",
" \n",
" test = pd.read_csv(test_fname)\n",
" \n",
" cyclical_attr = [\n",
" 'Month', 'Day', 'DayOfWeek', 'Hour'\n",
" ]\n",
" test = encode_categorical_attributes(test, cyclical_attr)\n",
" \n",
" date_attr = [\n",
" 'Quarter', 'Triannual', 'Semester', 'Fortnight',\n",
" 'FourHourPeriod', 'SixHourPeriod', 'TwelveHourPeriod'\n",
" ]\n",
" test = encode_categorical_attributes(test, date_attr)\n",
" \n",
" test = encode_categorical_attributes(test, ['PdDistrict'])\n",
" \n",
" num_attr, cat_attr = identify_attributes(test, 'all', time_window, 'cumulative')\n",
" \n",
" X_cat = test[cat_attr].to_numpy()\n",
" \n",
" for i in range(n_estimators):\n",
" clf_pred_fname = os.path.join(predictions_path, 'clf-{}-pred.csv'.format(i))\n",
" if os.path.isfile(clf_pred_fname):\n",
" continue\n",
" \n",
" scaler_fname = os.path.join(models_path, 'clf-{}-scaler.joblib'.format(i))\n",
" scaler = joblib.load(scaler_fname)\n",
" \n",
" X_num = scaler.transform(test[num_attr].to_numpy().astype(float))\n",
" \n",
" X = np.hstack([X_cat, X_num])\n",
" \n",
" clf_fname = os.path.join(models_path, 'clf-{}.joblib'.format(i))\n",
" clf = joblib.load(clf_fname)\n",
" \n",
" predictions = pd.DataFrame(clf.predict_proba(X), columns=prediction_header)\n",
" \n",
" predictions['Id'] = test['Id'].values\n",
" predictions = predictions[['Id']+prediction_header]\n",
" \n",
" predictions.to_csv(clf_pred_fname, float_format='%.5f', index=False)"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [],
"source": [
"def second_level_predictions(\n",
" n_estimators, stack_method,\n",
" models_path, predictions_path,\n",
" out_fname\n",
" ):\n",
" \"\"\"Combine the predictions from the base models.\"\"\"\n",
" prediction_header = ['' for j in range(len(CRIME_CATEGORY))]\n",
" for crime, idx in CRIME_CATEGORY.items():\n",
" prediction_header[idx] = crime\n",
" \n",
" X = None\n",
" \n",
" identifiers = None\n",
" \n",
" for i in range(n_estimators):\n",
" pred_fname = os.path.join(predictions_path, 'clf-{}-pred.csv'.format(i))\n",
" clf_pred = np.loadtxt(pred_fname, dtype=float, delimiter=',', skiprows=1)\n",
" \n",
" if identifiers is None:\n",
" identifiers = clf_pred[:,0].astype(int).tolist()\n",
" \n",
" clf_pred = clf_pred[:,1:]\n",
" \n",
" if X is None:\n",
" X = copy.deepcopy(clf_pred)\n",
" continue\n",
" \n",
" X = np.add(X, clf_pred) if stack_method == 'soft_voting' else np.hstack([X, clf_pred])\n",
" \n",
" predictions = None\n",
" \n",
" if stack_method == 'soft_voting':\n",
" predictions = X / n_estimators\n",
" else:\n",
" clf_ensemble_fname = os.path.join(models_path, 'clf-ensemble.joblib')\n",
" clf_ensemble = joblib.load(clf_ensemble_fname)\n",
" \n",
" predictions = clf_ensemble.predict_proba(X)\n",
" \n",
" predictions = pd.DataFrame(predictions, columns=prediction_header)\n",
" \n",
" predictions['Id'] = identifiers\n",
" predictions = predictions[['Id']+prediction_header]\n",
" \n",
" predictions.to_csv(out_fname, float_format='%.5f', index=False)"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [],
"source": [
"def late_submission(\n",
" time_window,\n",
" n_estimators,\n",
" stack_method,\n",
" base_algorithm=ALGORITHMS['logit'],\n",
" meta_algorithm=ALGORITHMS['logit']\n",
" ):\n",
" \"\"\"Make final predictions.\n",
" \n",
" These predictions are made on the test set.\n",
" \"\"\"\n",
" window_path = {'/': os.path.join(DATA_PATH, '{}H'.format(time_window))}\n",
" \n",
" window_path['/data'] = os.path.join(window_path['/'], 'data')\n",
" \n",
" window_path['/model'] = os.path.join(\n",
" window_path['/'],\n",
" 'model',\n",
" 'n_estimators={}'.format(n_estimators),\n",
" 'cyclical_attr=onehot'\n",
" )\n",
" \n",
" window_path['/prediction'] = os.path.join(\n",
" window_path['/'],\n",
" 'prediction',\n",
" 'test',\n",
" 'n_estimators={}'.format(n_estimators),\n",
" 'cyclical_attr=onehot'\n",
" )\n",
" \n",
" train_base_models = False\n",
" make_first_level_pred = False\n",
" \n",
" for i in range(n_estimators):\n",
" clf_fname = os.path.join(window_path['/model'], 'clf-{}.joblib'.format(i))\n",
" clf_pred_fname = os.path.join(window_path['/prediction'], 'clf-{}-pred.csv'.format(i))\n",
" \n",
" if not os.path.isfile(clf_fname):\n",
" train_base_models = True\n",
" make_first_level_pred = True\n",
" elif not os.path.isfile(clf_pred_fname):\n",
" make_first_level_pred = True\n",
" \n",
" if train_base_models:\n",
" break\n",
" \n",
" if train_base_models:\n",
" train_and_persist_base_models(\n",
" os.path.join(window_path['/data'], 'train.csv'),\n",
" os.path.join(window_path['/data'], 'validation.csv'),\n",
" time_window,\n",
" n_estimators,\n",
" base_algorithm,\n",
" meta_algorithm,\n",
" window_path['/model']\n",
" )\n",
" \n",
" if make_first_level_pred:\n",
" first_level_predictions(\n",
" os.path.join(window_path['/data'], 'test.csv'),\n",
" time_window,\n",
" n_estimators,\n",
" window_path['/model'],\n",
" window_path['/prediction']\n",
" )\n",
" \n",
" out_fname = os.path.join(window_path['/prediction'], '{}-ensemble-pred.csv'.format(stack_method))\n",
" if not os.path.isfile(out_fname):\n",
" second_level_predictions(\n",
" n_estimators, stack_method,\n",
" window_path['/model'], window_path['/prediction'],\n",
" out_fname\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [],
"source": [
"time_window = 336\n",
"n_estimators = 7\n",
"\n",
"late_submission(time_window, n_estimators, 'clf')\n",
"late_submission(time_window, n_estimators, 'soft_voting')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, let's combine the predictions from the two best meta-models and produce third-level ones. To accomplish this, second-level predictions are combined through soft voting."
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [],
"source": [
"def third_level_predictions(\n",
" time_window,\n",
" n_estimators,\n",
" stack_methods=['clf', 'soft_voting']\n",
" ):\n",
" \"\"\"Produce third-level predictions by combining those from the two best meta-models.\"\"\"\n",
" window_path = {'/': os.path.join(DATA_PATH, '{}H'.format(time_window))}\n",
" \n",
" window_path['/prediction'] = os.path.join(\n",
" window_path['/'],\n",
" 'prediction',\n",
" 'test',\n",
" 'n_estimators={}'.format(n_estimators),\n",
" 'cyclical_attr=onehot'\n",
" )\n",
" \n",
" X = None\n",
" \n",
" identifiers = None\n",
"\n",
" for stack_method in ['clf', 'soft_voting']:\n",
" pred_fname = os.path.join(window_path['/prediction'], '{}-ensemble-pred.csv'.format(stack_method))\n",
" clf_pred = np.loadtxt(pred_fname, dtype=float, delimiter=',', skiprows=1)\n",
" \n",
" if identifiers is None:\n",
" identifiers = clf_pred[:,0].astype(int).tolist()\n",
" \n",
" clf_pred = clf_pred[:,1:]\n",
" \n",
" X = copy.deepcopy(clf_pred) if X is None else np.add(X, clf_pred)\n",
" \n",
" predictions = X / 2\n",
" \n",
" prediction_header = ['' for j in range(len(CRIME_CATEGORY))]\n",
" for crime, idx in CRIME_CATEGORY.items():\n",
" prediction_header[idx] = crime\n",
"\n",
" predictions = pd.DataFrame(predictions, columns=prediction_header)\n",
" \n",
" predictions['Id'] = identifiers\n",
" predictions = predictions[['Id']+prediction_header]\n",
"\n",
" out_fname = os.path.join(window_path['/prediction'], 'third-level-pred.csv')\n",
" predictions.to_csv(out_fname, float_format='%.5f', index=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"| Meta-model Description | Log Loss |\n",
"|------------------------------|----------|\n",
"| The first-ranked meta-model | 2.56294 |\n",
"| The second-ranked meta-model | 2.56539 |\n",
"| Third-level predictions | 2.56044 |"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The above table shows the multi-class logarithmic loss of the predictions from the three meta-models on the test set. In the first place, results show that the two best meta-models learned to accurately generalize on unseen data (i.e., the test set), as the overall difference between validation scores and test scores is really small.\n",
"\n",
"In the second place, training a classifier to combine predictions from base models outperforms the straightforward technique of soft voting. The final decrease in *log loss* is around 0.1% when a classifier is used to produce second-level predictions.\n",
"\n",
"More importantly, producing third-level predictions remarkably outperforms those from the two best meta-models. In particular, the decreases in *log loss* are 0.19% and 0.098% when compared to the performance of the second- and first-ranked meta-models, respectively.\n",
"\n",
"Overall, these results are satisfactory, although the best result achieved is relatively larger than the best-reported result for the competition, i.e., 1.95936."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4. Conclusion"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. Throughout this document, a feature engineering strategy has been developed to support the predictive power of an ensemble-based approach to crime classification.\n",
"2. To each record in the dataset, such a strategy created several features from the history of crimes recently committed and within the area of influence.\n",
"3. Then, several base models were trained to make first-level predictions and these latter combined to produce second-level predictions.\n",
"4. Model averaging, or soft voting, is a simple but powerful stacking technique, as most of the top meta-models were built on it.\n",
"5. However, training a classifier to combine the predictions from the base models outperformed the straightforward stacking technique of soft voting.\n",
"6. On the other hand, the larger the number of base models, the better the predictive power of the meta-model.\n",
"7. In a like manner, the larger the time window, the better the predictive power of the meta-model, as most of the top meta-models used a time window of at least 168 hours to set the recency, i.e., the history of recently committed crimes.\n",
"8. The best setting to build meta-models was to train seven base models, use a time window of 336 hours, and one-hot encode cyclical attributes.\n",
"9. Finally, the best score obtained, i.e., 2.56044 of multi-class logarithmic loss, was the result of combining the second-level predictions."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}