src/Q5_norm_train.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "590cf583-4662-43fa-853d-ae9fbdbfdfca",
   "metadata": {},
   "source": [
    "# Predicting target variables using different models\n",
    "\n",
    "### Prepare the training set (2017) and test set (2019) using \n",
    "https://github.com/sagerpascal/uzh-data-science-project/blob/main/WZ/Q5_datacleaning.ipynb\n",
    "\n",
    "### Data downloaded from:\n",
    "https://data.stadt-zuerich.ch/dataset/vbz_fahrgastzahlen_ogd <br>\n",
    "and <br>\n",
    "https://data.stadt-zuerich.ch/dataset/vbz_fahrzeiten_ogd\n",
    "\n",
    "This script documents trials to predict VBZ public transportation vehicle occupancies (0 to 1) using data on year 2017 and validation on year 2019. Those years are chosen for the similarity in traffic volumes associated with different vehicle types.\n",
    "\n",
    "In the output file from <br>\n",
    "https://github.com/sagerpascal/uzh-data-science-project/blob/main/src/Q5_datacleaning.ipynb \n",
    "\n",
    "We have feature variables in:\n",
    "- Time\n",
    "- GPS_Latitude\n",
    "- GPS_Longitude\n",
    "- Weekday\n",
    "- Vehicle_type_ <br>\n",
    "and additionally information about the maximum capacity of the vehicle, which is reserved for calculating the number of predicted empty seat ((1- Occupany) * Capacity)\n",
    "\n",
    "The predictor variable is set to y = Occupancy.\n",
    "\n",
    "Some trials were done using all the feature variables, while especially for svr because the training takes too long, the dataset was reduced to tram only plus PCA reducing the dimension.\n",
    "\n",
    "As result of this training, we obtained models with prediction precision at:\n",
    "\n",
    "> Linear Regression MAPE on test set: 2144734221.059043 <br>\n",
    "> Linear Regression MAE on test set: 0.4593566554308329 <br>\n",
    "> Linear Regression MSE on test set: 0.4039692366610636 <br>\n",
    "> \n",
    "> Random Forest MAPE on test set: 734290880.1453922 <br>\n",
    "> Random Forest MAE on test set: 0.36878384335316794 <br>\n",
    "> Random Forest MSE on test set: 0.2893897954651777 <br>\n",
    "> \n",
    "> SVR MAPE on test set: inf <br>\n",
    "> SVR MAE on test set: 0.4264129866667982 <br>\n",
    "> SVR MSE on test set: 0.40956721848620736 <br>\n",
    "\n",
    "with the mean absolute error of the median predictor being: 0.687056532327459\n",
    "\n",
    "For discussions, please take a look at out presentation (link in readme).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "db560b56-5ef5-4070-82b3-bed0982dbb30",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "# import umap\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.svm import SVR\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.metrics import mean_absolute_error, mean_squared_error\n",
    "from datetime import datetime\n",
    "from joblib import dump"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "cea3d823-bd93-4d2a-8df2-89f0f0e69a93",
   "metadata": {},
   "outputs": [],
   "source": [
    "trainset = pd.read_csv('../data/fahrgastzahlen_2017/cleaned.csv', sep=',')\n",
    "testset = pd.read_csv('../data/fahrgastzahlen_2019/cleaned.csv', sep=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "fd14b62f-f1e3-4e70-932e-4fcb2ff5caf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Too large (particularly for svr), keep only tram\n",
    "trainset = trainset[trainset['Vehicle_type_T'] == True]\n",
    "trainset = trainset.drop(columns = ['Nachtnetz', 'Vehicle_type_B', 'Vehicle_type_FB', 'Vehicle_type_N', \n",
    "                                    'Vehicle_type_SB', 'Vehicle_type_TR', 'Vehicle_type_T'])\n",
    "testset = testset[testset['Vehicle_type_T'] == True]\n",
    "testset = testset.drop(columns = ['Nachtnetz', 'Vehicle_type_B', 'Vehicle_type_FB', 'Vehicle_type_N', \n",
    "                                    'Vehicle_type_SB', 'Vehicle_type_TR', 'Vehicle_type_T'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "d48957a9-4395-4419-a133-e6ad7e4bb807",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Index: 303949 entries, 619 to 576462\n",
      "Data columns (total 6 columns):\n",
      " #   Column         Non-Null Count   Dtype  \n",
      "---  ------         --------------   -----  \n",
      " 0   Time           303949 non-null  float64\n",
      " 1   Capacity       303949 non-null  float64\n",
      " 2   Occupancy      303949 non-null  float64\n",
      " 3   GPS_Latitude   303949 non-null  float64\n",
      " 4   GPS_Longitude  303949 non-null  float64\n",
      " 5   Weekday        303949 non-null  int64  \n",
      "dtypes: float64(5), int64(1)\n",
      "memory usage: 16.2 MB\n"
     ]
    }
   ],
   "source": [
    "testset.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "62991ede-31c5-4e75-8032-c53b9e6cbcee",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Index: 199606 entries, 6680 to 374981\n",
      "Data columns (total 6 columns):\n",
      " #   Column         Non-Null Count   Dtype  \n",
      "---  ------         --------------   -----  \n",
      " 0   Time           199606 non-null  float64\n",
      " 1   Capacity       199606 non-null  float64\n",
      " 2   Occupancy      199606 non-null  float64\n",
      " 3   GPS_Latitude   199606 non-null  float64\n",
      " 4   GPS_Longitude  199606 non-null  float64\n",
      " 5   Weekday        199606 non-null  int64  \n",
      "dtypes: float64(5), int64(1)\n",
      "memory usage: 10.7 MB\n"
     ]
    }
   ],
   "source": [
    "trainset.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "f2f55fcb-b22d-4229-9f79-ff3af947aa94",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit the scaler on the training set only\n",
    "scaler = StandardScaler()\n",
    "columns_to_normalize = ['Time', 'GPS_Latitude', 'GPS_Longitude',\n",
    "                        'Occupancy' \n",
    "                    # 'Capacity' only to use to restore the number of empty seats\n",
    "                       ]\n",
    "scaler.fit(trainset[columns_to_normalize])\n",
    "\n",
    "# Transform both training and test sets\n",
    "trainset[columns_to_normalize] = scaler.transform(trainset[columns_to_normalize])\n",
    "testset[columns_to_normalize] = scaler.transform(testset[columns_to_normalize])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "8d8906a8-0c11-4cf1-beae-882c4884f479",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Time</th>\n",
       "      <th>Capacity</th>\n",
       "      <th>Occupancy</th>\n",
       "      <th>GPS_Latitude</th>\n",
       "      <th>GPS_Longitude</th>\n",
       "      <th>Weekday</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>619</th>\n",
       "      <td>-1.544443</td>\n",
       "      <td>96.0</td>\n",
       "      <td>-1.095871</td>\n",
       "      <td>-0.626557</td>\n",
       "      <td>-0.924501</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>620</th>\n",
       "      <td>-1.544443</td>\n",
       "      <td>32.0</td>\n",
       "      <td>-0.259517</td>\n",
       "      <td>-0.626557</td>\n",
       "      <td>-0.924501</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>621</th>\n",
       "      <td>-1.544443</td>\n",
       "      <td>96.0</td>\n",
       "      <td>-0.688532</td>\n",
       "      <td>-0.626557</td>\n",
       "      <td>-0.924501</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>622</th>\n",
       "      <td>-1.365261</td>\n",
       "      <td>32.0</td>\n",
       "      <td>0.228457</td>\n",
       "      <td>-0.626557</td>\n",
       "      <td>-0.924501</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>623</th>\n",
       "      <td>-0.917307</td>\n",
       "      <td>96.0</td>\n",
       "      <td>-0.296973</td>\n",
       "      <td>-0.626557</td>\n",
       "      <td>-0.924501</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>624</th>\n",
       "      <td>-0.738125</td>\n",
       "      <td>32.0</td>\n",
       "      <td>0.920881</td>\n",
       "      <td>-0.626557</td>\n",
       "      <td>-0.924501</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>625</th>\n",
       "      <td>-0.738125</td>\n",
       "      <td>96.0</td>\n",
       "      <td>-0.537839</td>\n",
       "      <td>-0.626557</td>\n",
       "      <td>-0.924501</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>626</th>\n",
       "      <td>-0.738125</td>\n",
       "      <td>48.0</td>\n",
       "      <td>0.227590</td>\n",
       "      <td>-0.626557</td>\n",
       "      <td>-0.924501</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>627</th>\n",
       "      <td>-0.648534</td>\n",
       "      <td>44.0</td>\n",
       "      <td>0.018900</td>\n",
       "      <td>-0.626557</td>\n",
       "      <td>-0.924501</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>628</th>\n",
       "      <td>-0.469352</td>\n",
       "      <td>32.0</td>\n",
       "      <td>0.453716</td>\n",
       "      <td>-0.626557</td>\n",
       "      <td>-0.924501</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Time  Capacity  Occupancy  GPS_Latitude  GPS_Longitude  Weekday\n",
       "619 -1.544443      96.0  -1.095871     -0.626557      -0.924501        1\n",
       "620 -1.544443      32.0  -0.259517     -0.626557      -0.924501        1\n",
       "621 -1.544443      96.0  -0.688532     -0.626557      -0.924501        1\n",
       "622 -1.365261      32.0   0.228457     -0.626557      -0.924501        1\n",
       "623 -0.917307      96.0  -0.296973     -0.626557      -0.924501        1\n",
       "624 -0.738125      32.0   0.920881     -0.626557      -0.924501        1\n",
       "625 -0.738125      96.0  -0.537839     -0.626557      -0.924501        1\n",
       "626 -0.738125      48.0   0.227590     -0.626557      -0.924501        1\n",
       "627 -0.648534      44.0   0.018900     -0.626557      -0.924501        1\n",
       "628 -0.469352      32.0   0.453716     -0.626557      -0.924501        1"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "testset.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "bbfa220e-b496-43bb-8a6b-0fdc59ad8e45",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Time</th>\n",
       "      <th>GPS_Latitude</th>\n",
       "      <th>GPS_Longitude</th>\n",
       "      <th>Weekday</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>6680</th>\n",
       "      <td>-1.634034</td>\n",
       "      <td>0.467447</td>\n",
       "      <td>-2.765563</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6681</th>\n",
       "      <td>-1.454852</td>\n",
       "      <td>0.467447</td>\n",
       "      <td>-2.765563</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6682</th>\n",
       "      <td>-1.454852</td>\n",
       "      <td>0.467447</td>\n",
       "      <td>-2.765563</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6683</th>\n",
       "      <td>-1.365261</td>\n",
       "      <td>0.467447</td>\n",
       "      <td>-2.765563</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6684</th>\n",
       "      <td>-1.365261</td>\n",
       "      <td>0.467447</td>\n",
       "      <td>-2.765563</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          Time  GPS_Latitude  GPS_Longitude  Weekday\n",
       "6680 -1.634034      0.467447      -2.765563        1\n",
       "6681 -1.454852      0.467447      -2.765563        1\n",
       "6682 -1.454852      0.467447      -2.765563        1\n",
       "6683 -1.365261      0.467447      -2.765563        1\n",
       "6684 -1.365261      0.467447      -2.765563        1"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train = trainset.drop(columns=['Capacity', 'Occupancy'])\n",
    "y_train = trainset['Occupancy']\n",
    "X_test = testset.drop(columns=['Capacity', 'Occupancy'])\n",
    "y_test = testset['Occupancy']\n",
    "X_train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "39906e86-0d01-43ba-af19-2e0a7e02b259",
   "metadata": {},
   "outputs": [],
   "source": [
    "# If neccesary, take a smaller training set\n",
    "X_train_s, X_val, y_train_s, y_val = train_test_split(X_train, y_train, test_size=0.8, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d80acf60-aa76-46d9-80a4-fce11913fc66",
   "metadata": {},
   "source": [
    "## Training models: linear regression and random forest with some GridSearchCV for best parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "29a5ca6e-a6af-42be-b276-221dfd091d26",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting 5 folds for each of 2 candidates, totalling 10 fits\n",
      "Fitting 5 folds for each of 2 candidates, totalling 10 fits\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 100 out of 100 | elapsed:    7.7s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 100 out of 100 | elapsed:    0.1s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 100 out of 100 | elapsed:    8.0s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 100 out of 100 | elapsed:    0.1s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 100 out of 100 | elapsed:    7.8s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 100 out of 100 | elapsed:    0.1s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 100 out of 100 | elapsed:    7.8s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 100 out of 100 | elapsed:    0.1s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 100 out of 100 | elapsed:    7.6s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 100 out of 100 | elapsed:    0.1s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:   15.6s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:    0.2s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:   16.0s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:    0.2s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:   15.7s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:    0.2s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:   15.3s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:    0.2s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:   15.2s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:    0.2s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:   19.7s finished\n",
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best parameters for Linear Regression: {'fit_intercept': False}\n",
      "Best parameters for Random Forest: {'max_depth': 10, 'n_estimators': 200}\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:   19.6s finished\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-1 {\n",
       "  /* Definition of color scheme common for light and dark mode */\n",
       "  --sklearn-color-text: black;\n",
       "  --sklearn-color-line: gray;\n",
       "  /* Definition of color scheme for unfitted estimators */\n",
       "  --sklearn-color-unfitted-level-0: #fff5e6;\n",
       "  --sklearn-color-unfitted-level-1: #f6e4d2;\n",
       "  --sklearn-color-unfitted-level-2: #ffe0b3;\n",
       "  --sklearn-color-unfitted-level-3: chocolate;\n",
       "  /* Definition of color scheme for fitted estimators */\n",
       "  --sklearn-color-fitted-level-0: #f0f8ff;\n",
       "  --sklearn-color-fitted-level-1: #d4ebff;\n",
       "  --sklearn-color-fitted-level-2: #b3dbfd;\n",
       "  --sklearn-color-fitted-level-3: cornflowerblue;\n",
       "\n",
       "  /* Specific color for light theme */\n",
       "  --sklearn-color-text-on-default-background: var(--sg-text-color, var(--theme-code-foreground, var(--jp-content-font-color1, black)));\n",
       "  --sklearn-color-background: var(--sg-background-color, var(--theme-background, var(--jp-layout-color0, white)));\n",
       "  --sklearn-color-border-box: var(--sg-text-color, var(--theme-code-foreground, var(--jp-content-font-color1, black)));\n",
       "  --sklearn-color-icon: #696969;\n",
       "\n",
       "  @media (prefers-color-scheme: dark) {\n",
       "    /* Redefinition of color scheme for dark theme */\n",
       "    --sklearn-color-text-on-default-background: var(--sg-text-color, var(--theme-code-foreground, var(--jp-content-font-color1, white)));\n",
       "    --sklearn-color-background: var(--sg-background-color, var(--theme-background, var(--jp-layout-color0, #111)));\n",
       "    --sklearn-color-border-box: var(--sg-text-color, var(--theme-code-foreground, var(--jp-content-font-color1, white)));\n",
       "    --sklearn-color-icon: #878787;\n",
       "  }\n",
       "}\n",
       "\n",
       "#sk-container-id-1 {\n",
       "  color: var(--sklearn-color-text);\n",
       "}\n",
       "\n",
       "#sk-container-id-1 pre {\n",
       "  padding: 0;\n",
       "}\n",
       "\n",
       "#sk-container-id-1 input.sk-hidden--visually {\n",
       "  border: 0;\n",
       "  clip: rect(1px 1px 1px 1px);\n",
       "  clip: rect(1px, 1px, 1px, 1px);\n",
       "  height: 1px;\n",
       "  margin: -1px;\n",
       "  overflow: hidden;\n",
       "  padding: 0;\n",
       "  position: absolute;\n",
       "  width: 1px;\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-dashed-wrapped {\n",
       "  border: 1px dashed var(--sklearn-color-line);\n",
       "  margin: 0 0.4em 0.5em 0.4em;\n",
       "  box-sizing: border-box;\n",
       "  padding-bottom: 0.4em;\n",
       "  background-color: var(--sklearn-color-background);\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-container {\n",
       "  /* jupyter's `normalize.less` sets `[hidden] { display: none; }`\n",
       "     but bootstrap.min.css set `[hidden] { display: none !important; }`\n",
       "     so we also need the `!important` here to be able to override the\n",
       "     default hidden behavior on the sphinx rendered scikit-learn.org.\n",
       "     See: https://github.com/scikit-learn/scikit-learn/issues/21755 */\n",
       "  display: inline-block !important;\n",
       "  position: relative;\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-text-repr-fallback {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       "div.sk-parallel-item,\n",
       "div.sk-serial,\n",
       "div.sk-item {\n",
       "  /* draw centered vertical line to link estimators */\n",
       "  background-image: linear-gradient(var(--sklearn-color-text-on-default-background), var(--sklearn-color-text-on-default-background));\n",
       "  background-size: 2px 100%;\n",
       "  background-repeat: no-repeat;\n",
       "  background-position: center center;\n",
       "}\n",
       "\n",
       "/* Parallel-specific style estimator block */\n",
       "\n",
       "#sk-container-id-1 div.sk-parallel-item::after {\n",
       "  content: \"\";\n",
       "  width: 100%;\n",
       "  border-bottom: 2px solid var(--sklearn-color-text-on-default-background);\n",
       "  flex-grow: 1;\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-parallel {\n",
       "  display: flex;\n",
       "  align-items: stretch;\n",
       "  justify-content: center;\n",
       "  background-color: var(--sklearn-color-background);\n",
       "  position: relative;\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-parallel-item {\n",
       "  display: flex;\n",
       "  flex-direction: column;\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-parallel-item:first-child::after {\n",
       "  align-self: flex-end;\n",
       "  width: 50%;\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-parallel-item:last-child::after {\n",
       "  align-self: flex-start;\n",
       "  width: 50%;\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-parallel-item:only-child::after {\n",
       "  width: 0;\n",
       "}\n",
       "\n",
       "/* Serial-specific style estimator block */\n",
       "\n",
       "#sk-container-id-1 div.sk-serial {\n",
       "  display: flex;\n",
       "  flex-direction: column;\n",
       "  align-items: center;\n",
       "  background-color: var(--sklearn-color-background);\n",
       "  padding-right: 1em;\n",
       "  padding-left: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "/* Toggleable style: style used for estimator/Pipeline/ColumnTransformer box that is\n",
       "clickable and can be expanded/collapsed.\n",
       "- Pipeline and ColumnTransformer use this feature and define the default style\n",
       "- Estimators will overwrite some part of the style using the `sk-estimator` class\n",
       "*/\n",
       "\n",
       "/* Pipeline and ColumnTransformer style (default) */\n",
       "\n",
       "#sk-container-id-1 div.sk-toggleable {\n",
       "  /* Default theme specific background. It is overwritten whether we have a\n",
       "  specific estimator or a Pipeline/ColumnTransformer */\n",
       "  background-color: var(--sklearn-color-background);\n",
       "}\n",
       "\n",
       "/* Toggleable label */\n",
       "#sk-container-id-1 label.sk-toggleable__label {\n",
       "  cursor: pointer;\n",
       "  display: block;\n",
       "  width: 100%;\n",
       "  margin-bottom: 0;\n",
       "  padding: 0.5em;\n",
       "  box-sizing: border-box;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       "#sk-container-id-1 label.sk-toggleable__label-arrow:before {\n",
       "  /* Arrow on the left of the label */\n",
       "  content: \"▸\";\n",
       "  float: left;\n",
       "  margin-right: 0.25em;\n",
       "  color: var(--sklearn-color-icon);\n",
       "}\n",
       "\n",
       "#sk-container-id-1 label.sk-toggleable__label-arrow:hover:before {\n",
       "  color: var(--sklearn-color-text);\n",
       "}\n",
       "\n",
       "/* Toggleable content - dropdown */\n",
       "\n",
       "#sk-container-id-1 div.sk-toggleable__content {\n",
       "  max-height: 0;\n",
       "  max-width: 0;\n",
       "  overflow: hidden;\n",
       "  text-align: left;\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-0);\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-toggleable__content.fitted {\n",
       "  /* fitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-0);\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-toggleable__content pre {\n",
       "  margin: 0.2em;\n",
       "  border-radius: 0.25em;\n",
       "  color: var(--sklearn-color-text);\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-0);\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-toggleable__content.fitted pre {\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-0);\n",
       "}\n",
       "\n",
       "#sk-container-id-1 input.sk-toggleable__control:checked~div.sk-toggleable__content {\n",
       "  /* Expand drop-down */\n",
       "  max-height: 200px;\n",
       "  max-width: 100%;\n",
       "  overflow: auto;\n",
       "}\n",
       "\n",
       "#sk-container-id-1 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {\n",
       "  content: \"▾\";\n",
       "}\n",
       "\n",
       "/* Pipeline/ColumnTransformer-specific style */\n",
       "\n",
       "#sk-container-id-1 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {\n",
       "  color: var(--sklearn-color-text);\n",
       "  background-color: var(--sklearn-color-unfitted-level-2);\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-label.fitted input.sk-toggleable__control:checked~label.sk-toggleable__label {\n",
       "  background-color: var(--sklearn-color-fitted-level-2);\n",
       "}\n",
       "\n",
       "/* Estimator-specific style */\n",
       "\n",
       "/* Colorize estimator box */\n",
       "#sk-container-id-1 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-2);\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-estimator.fitted input.sk-toggleable__control:checked~label.sk-toggleable__label {\n",
       "  /* fitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-2);\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-label label.sk-toggleable__label,\n",
       "#sk-container-id-1 div.sk-label label {\n",
       "  /* The background is the default theme color */\n",
       "  color: var(--sklearn-color-text-on-default-background);\n",
       "}\n",
       "\n",
       "/* On hover, darken the color of the background */\n",
       "#sk-container-id-1 div.sk-label:hover label.sk-toggleable__label {\n",
       "  color: var(--sklearn-color-text);\n",
       "  background-color: var(--sklearn-color-unfitted-level-2);\n",
       "}\n",
       "\n",
       "/* Label box, darken color on hover, fitted */\n",
       "#sk-container-id-1 div.sk-label.fitted:hover label.sk-toggleable__label.fitted {\n",
       "  color: var(--sklearn-color-text);\n",
       "  background-color: var(--sklearn-color-fitted-level-2);\n",
       "}\n",
       "\n",
       "/* Estimator label */\n",
       "\n",
       "#sk-container-id-1 div.sk-label label {\n",
       "  font-family: monospace;\n",
       "  font-weight: bold;\n",
       "  display: inline-block;\n",
       "  line-height: 1.2em;\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-label-container {\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       "/* Estimator-specific */\n",
       "#sk-container-id-1 div.sk-estimator {\n",
       "  font-family: monospace;\n",
       "  border: 1px dotted var(--sklearn-color-border-box);\n",
       "  border-radius: 0.25em;\n",
       "  box-sizing: border-box;\n",
       "  margin-bottom: 0.5em;\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-0);\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-estimator.fitted {\n",
       "  /* fitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-0);\n",
       "}\n",
       "\n",
       "/* on hover */\n",
       "#sk-container-id-1 div.sk-estimator:hover {\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-2);\n",
       "}\n",
       "\n",
       "#sk-container-id-1 div.sk-estimator.fitted:hover {\n",
       "  /* fitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-2);\n",
       "}\n",
       "\n",
       "/* Specification for estimator info (e.g. \"i\" and \"?\") */\n",
       "\n",
       "/* Common style for \"i\" and \"?\" */\n",
       "\n",
       ".sk-estimator-doc-link,\n",
       "a:link.sk-estimator-doc-link,\n",
       "a:visited.sk-estimator-doc-link {\n",
       "  float: right;\n",
       "  font-size: smaller;\n",
       "  line-height: 1em;\n",
       "  font-family: monospace;\n",
       "  background-color: var(--sklearn-color-background);\n",
       "  border-radius: 1em;\n",
       "  height: 1em;\n",
       "  width: 1em;\n",
       "  text-decoration: none !important;\n",
       "  margin-left: 1ex;\n",
       "  /* unfitted */\n",
       "  border: var(--sklearn-color-unfitted-level-1) 1pt solid;\n",
       "  color: var(--sklearn-color-unfitted-level-1);\n",
       "}\n",
       "\n",
       ".sk-estimator-doc-link.fitted,\n",
       "a:link.sk-estimator-doc-link.fitted,\n",
       "a:visited.sk-estimator-doc-link.fitted {\n",
       "  /* fitted */\n",
       "  border: var(--sklearn-color-fitted-level-1) 1pt solid;\n",
       "  color: var(--sklearn-color-fitted-level-1);\n",
       "}\n",
       "\n",
       "/* On hover */\n",
       "div.sk-estimator:hover .sk-estimator-doc-link:hover,\n",
       ".sk-estimator-doc-link:hover,\n",
       "div.sk-label-container:hover .sk-estimator-doc-link:hover,\n",
       ".sk-estimator-doc-link:hover {\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-3);\n",
       "  color: var(--sklearn-color-background);\n",
       "  text-decoration: none;\n",
       "}\n",
       "\n",
       "div.sk-estimator.fitted:hover .sk-estimator-doc-link.fitted:hover,\n",
       ".sk-estimator-doc-link.fitted:hover,\n",
       "div.sk-label-container:hover .sk-estimator-doc-link.fitted:hover,\n",
       ".sk-estimator-doc-link.fitted:hover {\n",
       "  /* fitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-3);\n",
       "  color: var(--sklearn-color-background);\n",
       "  text-decoration: none;\n",
       "}\n",
       "\n",
       "/* Span, style for the box shown on hovering the info icon */\n",
       ".sk-estimator-doc-link span {\n",
       "  display: none;\n",
       "  z-index: 9999;\n",
       "  position: relative;\n",
       "  font-weight: normal;\n",
       "  right: .2ex;\n",
       "  padding: .5ex;\n",
       "  margin: .5ex;\n",
       "  width: min-content;\n",
       "  min-width: 20ex;\n",
       "  max-width: 50ex;\n",
       "  color: var(--sklearn-color-text);\n",
       "  box-shadow: 2pt 2pt 4pt #999;\n",
       "  /* unfitted */\n",
       "  background: var(--sklearn-color-unfitted-level-0);\n",
       "  border: .5pt solid var(--sklearn-color-unfitted-level-3);\n",
       "}\n",
       "\n",
       ".sk-estimator-doc-link.fitted span {\n",
       "  /* fitted */\n",
       "  background: var(--sklearn-color-fitted-level-0);\n",
       "  border: var(--sklearn-color-fitted-level-3);\n",
       "}\n",
       "\n",
       ".sk-estimator-doc-link:hover span {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       "/* \"?\"-specific style due to the `<a>` HTML tag */\n",
       "\n",
       "#sk-container-id-1 a.estimator_doc_link {\n",
       "  float: right;\n",
       "  font-size: 1rem;\n",
       "  line-height: 1em;\n",
       "  font-family: monospace;\n",
       "  background-color: var(--sklearn-color-background);\n",
       "  border-radius: 1rem;\n",
       "  height: 1rem;\n",
       "  width: 1rem;\n",
       "  text-decoration: none;\n",
       "  /* unfitted */\n",
       "  color: var(--sklearn-color-unfitted-level-1);\n",
       "  border: var(--sklearn-color-unfitted-level-1) 1pt solid;\n",
       "}\n",
       "\n",
       "#sk-container-id-1 a.estimator_doc_link.fitted {\n",
       "  /* fitted */\n",
       "  border: var(--sklearn-color-fitted-level-1) 1pt solid;\n",
       "  color: var(--sklearn-color-fitted-level-1);\n",
       "}\n",
       "\n",
       "/* On hover */\n",
       "#sk-container-id-1 a.estimator_doc_link:hover {\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-3);\n",
       "  color: var(--sklearn-color-background);\n",
       "  text-decoration: none;\n",
       "}\n",
       "\n",
       "#sk-container-id-1 a.estimator_doc_link.fitted:hover {\n",
       "  /* fitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-3);\n",
       "}\n",
       "</style><div id=\"sk-container-id-1\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>RandomForestRegressor(max_depth=10, n_estimators=200, verbose=1)</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator fitted sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-1\" type=\"checkbox\" checked><label for=\"sk-estimator-id-1\" class=\"sk-toggleable__label fitted sk-toggleable__label-arrow fitted\">&nbsp;&nbsp;RandomForestRegressor<a class=\"sk-estimator-doc-link fitted\" rel=\"noreferrer\" target=\"_blank\" href=\"https://scikit-learn.org/1.4/modules/generated/sklearn.ensemble.RandomForestRegressor.html\">?<span>Documentation for RandomForestRegressor</span></a><span class=\"sk-estimator-doc-link fitted\">i<span>Fitted</span></span></label><div class=\"sk-toggleable__content fitted\"><pre>RandomForestRegressor(max_depth=10, n_estimators=200, verbose=1)</pre></div> </div></div></div></div>"
      ],
      "text/plain": [
       "RandomForestRegressor(max_depth=10, n_estimators=200, verbose=1)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Linear regression and random forest\n",
    "# Define parameter grids for GridSearchCV, linear regression and random forest\n",
    "param_grid_lr = {'fit_intercept': [True, False]}\n",
    "# param_grid_rf = {'n_estimators': [100, 200, 300], 'max_depth': [None, 10, 20]}\n",
    "param_grid_rf = {'n_estimators': [100, 200], 'max_depth': [10]}\n",
    "\n",
    "# Initialize models, verbose=1 to show progress\n",
    "lr = LinearRegression()\n",
    "rf = RandomForestRegressor(verbose=1)\n",
    "\n",
    "# GridSearchCV for Linear Regression\n",
    "grid_lr = GridSearchCV(lr, param_grid_lr, cv=5, verbose=1)\n",
    "grid_lr.fit(X_train, y_train)\n",
    "\n",
    "# GridSearchCV for Random Forest\n",
    "grid_rf = GridSearchCV(rf, param_grid_rf, cv=5, verbose=1)\n",
    "grid_rf.fit(X_train, y_train)\n",
    "\n",
    "print(\"Best parameters for Linear Regression:\", grid_lr.best_params_)\n",
    "print(\"Best parameters for Random Forest:\", grid_rf.best_params_)\n",
    "\n",
    "# Fit the training data with best parameters\n",
    "best_lr = grid_lr.best_estimator_\n",
    "best_lr.fit(X_train, y_train)\n",
    "\n",
    "best_rf = grid_rf.best_estimator_\n",
    "best_rf.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1d8c6caf-c097-4bba-86b0-d14a23a93666",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['models/random_forest_model_2024-05-14_22-18-31.joblib']"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Dump save trained model\n",
    "# Get the current date and time\n",
    "timestamp = datetime.now().strftime(\"%Y-%m-%d_%H-%M-%S\")\n",
    "\n",
    "directory_path = 'models/'\n",
    "\n",
    "lr_model_filename = f'{directory_path}linear_regression_model_{timestamp}.joblib'\n",
    "rf_model_filename = f'{directory_path}random_forest_model_{timestamp}.joblib'\n",
    "\n",
    "dump(best_lr, lr_model_filename)\n",
    "dump(best_rf, rf_model_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "19bb0f41-825e-416b-aaf3-319fbedde6f7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean absolute error of the median predictor is: 35.21378499078745\n"
     ]
    }
   ],
   "source": [
    "# Calculate the median of the predictions\n",
    "median_y_train = np.median(trainset[\"Occupancy\"].to_numpy() * trainset[\"Capacity\"].to_numpy())\n",
    "\n",
    "# Calculate the mean absolute error (MAE) of the median predictor\n",
    "mae_med = np.mean(np.abs(trainset[\"Occupancy\"].to_numpy() * trainset[\"Capacity\"].to_numpy() - median_y_train))\n",
    "print(f\"The mean absolute error of the median predictor is: {mae_med}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "59ebd7ad-864d-4b1f-87e4-2d0312f7c800",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Linear Regression MAPE on test set: 108200656574.4099\n",
      "Linear Regression MAE on test set: 23.120253718431574\n",
      "Linear Regression MSE on test set: 880.9890201909976\n",
      "Random Forest MAPE on test set: 36956074996.45427\n",
      "Random Forest MAE on test set: 18.853875447532324\n",
      "Random Forest MSE on test set: 680.6717909617971\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:    1.7s finished\n"
     ]
    }
   ],
   "source": [
    "# Check the metrics for predicting on the test set\n",
    "# Predict on the test set\n",
    "y_pred_lr = best_lr.predict(X_test)\n",
    "y_pred_rf = best_rf.predict(X_test)\n",
    "\n",
    "# Have to do this becasue the scaler was fitted on 4d data\n",
    "testset_lr = testset.copy()\n",
    "testset_rf = testset.copy()\n",
    "testset_lr['Occupancy'] = y_pred_lr\n",
    "testset_rf['Occupancy'] = y_pred_rf\n",
    "\n",
    "# Reverse scale the y\n",
    "testset_re = scaler.inverse_transform(testset[columns_to_normalize])\n",
    "testset_lr_re = scaler.inverse_transform(testset_lr[columns_to_normalize])\n",
    "testset_rf_re = scaler.inverse_transform(testset_rf[columns_to_normalize])\n",
    "\n",
    "y_test_re = testset_re[:, 3]\n",
    "y_test_re = y_test_re * testset[\"Capacity\"].to_numpy()\n",
    "y_pred_lr_re = testset_lr_re[:, 3]\n",
    "y_pred_lr_re = y_pred_lr_re * testset[\"Capacity\"].to_numpy()\n",
    "y_pred_rf_re = testset_rf_re[:, 3]\n",
    "y_pred_rf_re = y_pred_rf_re * testset[\"Capacity\"].to_numpy()\n",
    "\n",
    "# Calculate MAPE with a small denominator to prevent inf\n",
    "mape_lr = np.mean(np.abs((y_test_re - y_pred_lr_re) / (y_test_re + 1e-10))) * 100\n",
    "mape_rf = np.mean(np.abs((y_test_re - y_pred_rf_re) / (y_test_re + 1e-10))) * 100\n",
    "\n",
    "# Calculate MAE \n",
    "mae_lr = mean_absolute_error(y_test_re, y_pred_lr_re)\n",
    "mae_rf = mean_absolute_error(y_test_re, y_pred_rf_re)\n",
    "\n",
    "# Calculate MSE\n",
    "mse_lr = mean_squared_error(y_test_re, y_pred_lr_re)\n",
    "mse_rf = mean_squared_error(y_test_re, y_pred_rf_re)\n",
    "\n",
    "print(\"Linear Regression MAPE on test set:\", mape_lr)\n",
    "print(\"Linear Regression MAE on test set:\", mae_lr)\n",
    "print(\"Linear Regression MSE on test set:\", mse_lr)\n",
    "\n",
    "print(\"Random Forest MAPE on test set:\", mape_rf)\n",
    "print(\"Random Forest MAE on test set:\", mae_rf)\n",
    "print(\"Random Forest MSE on test set:\", mse_rf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "263571f2-d89b-40d2-a084-719cdb5d313f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([1.27798e+05, 9.96630e+04, 4.91590e+04, 1.81450e+04, 6.26400e+03,\n",
       "        2.02300e+03, 6.57000e+02, 1.72000e+02, 6.00000e+01, 8.00000e+00]),\n",
       " array([  0. ,  22.9,  45.8,  68.7,  91.6, 114.5, 137.4, 160.3, 183.2,\n",
       "        206.1, 229. ]),\n",
       " <BarContainer object of 10 artists>)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGdCAYAAAD+JxxnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAsiElEQVR4nO3dfXAUdYL/8U/MwxhSSRsIyTAaNXvFZcmGczW4IeAu7AIJHiFneXfgRuegjop4QXI5wgqUtz/Rqk0QEbeOrK561uIiXvwD41kHZBOfojkIxEhuCeJDlWDCkiGsDJMQ4ySG/v1h0bVDeNQJIfm+X1Vd5XR/pufb047z8Ts9kwjbtm0BAAAY6JrhHgAAAMBwoQgBAABjUYQAAICxKEIAAMBYFCEAAGAsihAAADAWRQgAABiLIgQAAIwVNdwDuNqdPn1aR48eVXx8vCIiIoZ7OAAA4BLYtq3u7m55PB5dc835530oQhdx9OhRpaamDvcwAADAt9De3q4bbrjhvNspQhcRHx8v6ZsnMiEhYZhHAwAALkVXV5dSU1Od9/HzoQhdxJmPwxISEihCAACMMBe7rIWLpQEAgLEoQgAAwFgUIQAAYCyKEAAAMBZFCAAAGIsiBAAAjEURAgAAxqIIAQAAY1GEAACAsShCAADAWBQhAABgLIoQAAAwFkUIAAAYiyIEAACMFTXcAzDdzau3D/cQLtvhdfOGewgAAIQFM0IAAMBYFCEAAGAsihAAADAWRQgAABiLIgQAAIxFEQIAAMaiCAEAAGNRhAAAgLEoQgAAwFgUIQAAYCyKEAAAMBZFCAAAGIsiBAAAjEURAgAAxqIIAQAAY1GEAACAsShCAADAWJddhN59913Nnz9fHo9HEREReu2115xt/f39WrVqlSZPnqy4uDh5PB790z/9k44ePRqyj2AwqOXLlyspKUlxcXEqKCjQkSNHQjJ+v19er1eWZcmyLHm9Xp08eTIk09bWpvnz5ysuLk5JSUkqKSlRX19fSGb//v2aMWOGYmNjdf311+uxxx6TbduXe9gAAGAUuuwi1NPTo1tuuUWVlZWDtn355Zf64IMP9Mtf/lIffPCBXn31VX3yyScqKCgIyZWWlqq6ulpVVVVqaGjQqVOnlJ+fr4GBASdTWFiolpYW1dTUqKamRi0tLfJ6vc72gYEBzZs3Tz09PWpoaFBVVZW2bdumsrIyJ9PV1aU5c+bI4/GoqalJmzZt0oYNG7Rx48bLPWwAADAKRdjfYXokIiJC1dXVuuuuu86baWpq0o9+9CN9/vnnuvHGGxUIBDR+/Hht2bJFCxculCQdPXpUqamp2rFjh/Ly8nTw4EFlZGSosbFR2dnZkqTGxkbl5OToo48+Unp6unbu3Kn8/Hy1t7fL4/FIkqqqqrR48WJ1dnYqISFBzzzzjNasWaNjx47J5XJJktatW6dNmzbpyJEjioiIuOgxdnV1ybIsBQIBJSQkfNun6rxuXr097PscaofXzRvuIQAAcEGX+v495NcIBQIBRURE6LrrrpMkNTc3q7+/X7m5uU7G4/EoMzNTu3btkiTt3r1blmU5JUiSpk6dKsuyQjKZmZlOCZKkvLw8BYNBNTc3O5kZM2Y4JehM5ujRozp8+PA5xxsMBtXV1RWyAACA0WlIi9BXX32l1atXq7Cw0GljPp9PMTExSkxMDMmmpKTI5/M5meTk5EH7S05ODsmkpKSEbE9MTFRMTMwFM2dun8mcraKiwrkuybIspaamXu5hAwCAEWLIilB/f7/uuecenT59Wk8//fRF87Zth3xUda6PrcKROfNJ4Pk+FluzZo0CgYCztLe3X3TsAABgZBqSItTf368FCxbo0KFDqqurC/lszu12q6+vT36/P+Q+nZ2dzmyN2+3WsWPHBu33+PHjIZmzZ3X8fr/6+/svmOns7JSkQTNFZ7hcLiUkJIQsAABgdAp7ETpTgj799FO98cYbGjduXMj2rKwsRUdHq66uzlnX0dGh1tZWTZs2TZKUk5OjQCCgvXv3Opk9e/YoEAiEZFpbW9XR0eFkamtr5XK5lJWV5WTefffdkK/U19bWyuPx6Oabbw73oQMAgBHmsovQqVOn1NLSopaWFknSoUOH1NLSora2Nn399df6h3/4B73//vvaunWrBgYG5PP55PP5nDJiWZaWLFmisrIyvfnmm9q3b5/uu+8+TZ48WbNnz5YkTZo0SXPnzlVRUZEaGxvV2NiooqIi5efnKz09XZKUm5urjIwMeb1e7du3T2+++aZWrlypoqIiZxansLBQLpdLixcvVmtrq6qrq1VeXq4VK1Zc0jfGAADA6BZ1uXd4//339dOf/tS5vWLFCknSokWLtHbtWr3++uuSpB/+8Ich93v77bc1c+ZMSdJTTz2lqKgoLViwQL29vZo1a5Y2b96syMhIJ79161aVlJQ43y4rKCgI+e2iyMhIbd++XcXFxZo+fbpiY2NVWFioDRs2OBnLslRXV6dly5ZpypQpSkxM1IoVK5wxAwAAs32n3xEyAb8jNBi/IwQAuNpdNb8jBAAAcLWiCAEAAGNRhAAAgLEoQgAAwFgUIQAAYCyKEAAAMBZFCAAAGIsiBAAAjEURAgAAxqIIAQAAY1GEAACAsShCAADAWBQhAABgLIoQAAAwFkUIAAAYiyIEAACMRRECAADGoggBAABjUYQAAICxKEIAAMBYFCEAAGAsihAAADAWRQgAABiLIgQAAIxFEQIAAMaiCAEAAGNRhAAAgLEoQgAAwFgUIQAAYCyKEAAAMBZFCAAAGIsiBAAAjEURAgAAxqIIAQAAY1GEAACAsShCAADAWBQhAABgLIoQAAAwFkUIAAAYiyIEAACMRRECAADGoggBAABjUYQAAICxKEIAAMBYl12E3n33Xc2fP18ej0cRERF67bXXQrbbtq21a9fK4/EoNjZWM2fO1IEDB0IywWBQy5cvV1JSkuLi4lRQUKAjR46EZPx+v7xeryzLkmVZ8nq9OnnyZEimra1N8+fPV1xcnJKSklRSUqK+vr6QzP79+zVjxgzFxsbq+uuv12OPPSbbti/3sAEAwCh02UWop6dHt9xyiyorK8+5ff369dq4caMqKyvV1NQkt9utOXPmqLu728mUlpaqurpaVVVVamho0KlTp5Sfn6+BgQEnU1hYqJaWFtXU1KimpkYtLS3yer3O9oGBAc2bN089PT1qaGhQVVWVtm3bprKyMifT1dWlOXPmyOPxqKmpSZs2bdKGDRu0cePGyz1sAAAwCkXY32F6JCIiQtXV1brrrrskfTMb5PF4VFpaqlWrVkn6ZvYnJSVFjz/+uJYuXapAIKDx48dry5YtWrhwoSTp6NGjSk1N1Y4dO5SXl6eDBw8qIyNDjY2Nys7OliQ1NjYqJydHH330kdLT07Vz507l5+ervb1dHo9HklRVVaXFixers7NTCQkJeuaZZ7RmzRodO3ZMLpdLkrRu3Tpt2rRJR44cUURExEWPsaurS5ZlKRAIKCEh4ds+Ved18+rtYd/nUDu8bt5wDwEAgAu61PfvsF4jdOjQIfl8PuXm5jrrXC6XZsyYoV27dkmSmpub1d/fH5LxeDzKzMx0Mrt375ZlWU4JkqSpU6fKsqyQTGZmplOCJCkvL0/BYFDNzc1OZsaMGU4JOpM5evSoDh8+fM5jCAaD6urqClkAAMDoFNYi5PP5JEkpKSkh61NSUpxtPp9PMTExSkxMvGAmOTl50P6Tk5NDMmc/TmJiomJiYi6YOXP7TOZsFRUVznVJlmUpNTX14gcOAABGpKih2OnZHznZtn3Rj6HOzpwrH47MmU8CzzeeNWvWaMWKFc7trq4uytBZ+DgPADBahHVGyO12Sxo829LZ2enMxLjdbvX19cnv918wc+zYsUH7P378eEjm7Mfx+/3q7++/YKazs1PS4FmrM1wulxISEkIWAAAwOoW1CKWlpcntdquurs5Z19fXp/r6ek2bNk2SlJWVpejo6JBMR0eHWltbnUxOTo4CgYD27t3rZPbs2aNAIBCSaW1tVUdHh5Opra2Vy+VSVlaWk3n33XdDvlJfW1srj8ejm2++OZyHDgAARqDLLkKnTp1SS0uLWlpaJH1zgXRLS4va2toUERGh0tJSlZeXq7q6Wq2trVq8eLHGjBmjwsJCSZJlWVqyZInKysr05ptvat++fbrvvvs0efJkzZ49W5I0adIkzZ07V0VFRWpsbFRjY6OKioqUn5+v9PR0SVJubq4yMjLk9Xq1b98+vfnmm1q5cqWKioqcWZzCwkK5XC4tXrxYra2tqq6uVnl5uVasWHFJ3xgDAACj22VfI/T+++/rpz/9qXP7zPU0ixYt0ubNm/XQQw+pt7dXxcXF8vv9ys7OVm1treLj4537PPXUU4qKitKCBQvU29urWbNmafPmzYqMjHQyW7duVUlJifPtsoKCgpDfLoqMjNT27dtVXFys6dOnKzY2VoWFhdqwYYOTsSxLdXV1WrZsmaZMmaLExEStWLEi5BogAABgru/0O0Im4HeERgculgYAswzL7wgBAACMJBQhAABgLIoQAAAwFkUIAAAYiyIEAACMRRECAADGoggBAABjUYQAAICxKEIAAMBYFCEAAGAsihAAADAWRQgAABiLIgQAAIxFEQIAAMaiCAEAAGNRhAAAgLEoQgAAwFgUIQAAYCyKEAAAMBZFCAAAGIsiBAAAjEURAgAAxqIIAQAAY1GEAACAsShCAADAWBQhAABgLIoQAAAwFkUIAAAYiyIEAACMRRECAADGoggBAABjUYQAAICxKEIAAMBYFCEAAGAsihAAADAWRQgAABiLIgQAAIxFEQIAAMaiCAEAAGNRhAAAgLEoQgAAwFgUIQAAYCyKEAAAMBZFCAAAGCvsRejrr7/Wv//7vystLU2xsbH63ve+p8cee0ynT592MrZta+3atfJ4PIqNjdXMmTN14MCBkP0Eg0EtX75cSUlJiouLU0FBgY4cORKS8fv98nq9sixLlmXJ6/Xq5MmTIZm2tjbNnz9fcXFxSkpKUklJifr6+sJ92AAAYAQKexF6/PHH9dvf/laVlZU6ePCg1q9fryeeeEKbNm1yMuvXr9fGjRtVWVmppqYmud1uzZkzR93d3U6mtLRU1dXVqqqqUkNDg06dOqX8/HwNDAw4mcLCQrW0tKimpkY1NTVqaWmR1+t1tg8MDGjevHnq6elRQ0ODqqqqtG3bNpWVlYX7sAEAwAgUYdu2Hc4d5ufnKyUlRS+88IKz7u///u81ZswYbdmyRbZty+PxqLS0VKtWrZL0zexPSkqKHn/8cS1dulSBQEDjx4/Xli1btHDhQknS0aNHlZqaqh07digvL08HDx5URkaGGhsblZ2dLUlqbGxUTk6OPvroI6Wnp2vnzp3Kz89Xe3u7PB6PJKmqqkqLFy9WZ2enEhISLno8XV1dsixLgUDgkvKX6+bV28O+Twx2eN284R4CAOAKutT377DPCN1xxx1688039cknn0iS/u///k8NDQ3627/9W0nSoUOH5PP5lJub69zH5XJpxowZ2rVrlySpublZ/f39IRmPx6PMzEwns3v3blmW5ZQgSZo6daosywrJZGZmOiVIkvLy8hQMBtXc3HzO8QeDQXV1dYUsAABgdIoK9w5XrVqlQCCg73//+4qMjNTAwIB+9atf6ec//7kkyefzSZJSUlJC7peSkqLPP//cycTExCgxMXFQ5sz9fT6fkpOTBz1+cnJySObsx0lMTFRMTIyTOVtFRYUeffTRyz1sAAAwAoV9RuiVV17RSy+9pJdfflkffPCBXnzxRW3YsEEvvvhiSC4iIiLktm3bg9ad7ezMufLfJvOX1qxZo0Ag4Czt7e0XHBMAABi5wj4j9Itf/EKrV6/WPffcI0maPHmyPv/8c1VUVGjRokVyu92SvpmtmTBhgnO/zs5OZ/bG7Xarr69Pfr8/ZFaos7NT06ZNczLHjh0b9PjHjx8P2c+ePXtCtvv9fvX39w+aKTrD5XLJ5XJ928MHAAAjSNhnhL788ktdc03obiMjI52vz6elpcntdquurs7Z3tfXp/r6eqfkZGVlKTo6OiTT0dGh1tZWJ5OTk6NAIKC9e/c6mT179igQCIRkWltb1dHR4WRqa2vlcrmUlZUV5iMHAAAjTdhnhObPn69f/epXuvHGG/WDH/xA+/bt08aNG/XP//zPkr75qKq0tFTl5eWaOHGiJk6cqPLyco0ZM0aFhYWSJMuytGTJEpWVlWncuHEaO3asVq5cqcmTJ2v27NmSpEmTJmnu3LkqKirSs88+K0m6//77lZ+fr/T0dElSbm6uMjIy5PV69cQTT+jEiRNauXKlioqKhuQbYAAAYGQJexHatGmTfvnLX6q4uFidnZ3yeDxaunSp/t//+39O5qGHHlJvb6+Ki4vl9/uVnZ2t2tpaxcfHO5mnnnpKUVFRWrBggXp7ezVr1ixt3rxZkZGRTmbr1q0qKSlxvl1WUFCgyspKZ3tkZKS2b9+u4uJiTZ8+XbGxsSosLNSGDRvCfdgAAGAECvvvCI02/I7Q6MDvCAGAWYbtd4QAAABGCooQAAAwFkUIAAAYiyIEAACMRRECAADGoggBAABjUYQAAICxKEIAAMBYFCEAAGAsihAAADAWRQgAABiLIgQAAIxFEQIAAMaiCAEAAGNRhAAAgLEoQgAAwFgUIQAAYCyKEAAAMBZFCAAAGIsiBAAAjEURAgAAxqIIAQAAY1GEAACAsShCAADAWBQhAABgLIoQAAAwFkUIAAAYiyIEAACMRRECAADGoggBAABjUYQAAICxKEIAAMBYFCEAAGAsihAAADAWRQgAABiLIgQAAIxFEQIAAMaiCAEAAGNRhAAAgLEoQgAAwFgUIQAAYCyKEAAAMBZFCAAAGIsiBAAAjDUkRehPf/qT7rvvPo0bN05jxozRD3/4QzU3NzvbbdvW2rVr5fF4FBsbq5kzZ+rAgQMh+wgGg1q+fLmSkpIUFxengoICHTlyJCTj9/vl9XplWZYsy5LX69XJkydDMm1tbZo/f77i4uKUlJSkkpIS9fX1DcVhAwCAESbsRcjv92v69OmKjo7Wzp079eGHH+rJJ5/Udddd52TWr1+vjRs3qrKyUk1NTXK73ZozZ466u7udTGlpqaqrq1VVVaWGhgadOnVK+fn5GhgYcDKFhYVqaWlRTU2Nampq1NLSIq/X62wfGBjQvHnz1NPTo4aGBlVVVWnbtm0qKysL92EDAIARKMK2bTucO1y9erX+93//V++99945t9u2LY/Ho9LSUq1atUrSN7M/KSkpevzxx7V06VIFAgGNHz9eW7Zs0cKFCyVJR48eVWpqqnbs2KG8vDwdPHhQGRkZamxsVHZ2tiSpsbFROTk5+uijj5Senq6dO3cqPz9f7e3t8ng8kqSqqiotXrxYnZ2dSkhIuOjxdHV1ybIsBQKBS8pfrptXbw/7PjHY4XXzhnsIAIAr6FLfv8M+I/T6669rypQp+sd//EclJyfr1ltv1fPPP+9sP3TokHw+n3Jzc511LpdLM2bM0K5duyRJzc3N6u/vD8l4PB5lZmY6md27d8uyLKcESdLUqVNlWVZIJjMz0ylBkpSXl6dgMBjyUd1fCgaD6urqClkAAMDoFPYi9Nlnn+mZZ57RxIkT9Yc//EEPPPCASkpK9Pvf/16S5PP5JEkpKSkh90tJSXG2+Xw+xcTEKDEx8YKZ5OTkQY+fnJwckjn7cRITExUTE+NkzlZRUeFcc2RZllJTUy/3KQAAACNE2IvQ6dOnddttt6m8vFy33nqrli5dqqKiIj3zzDMhuYiIiJDbtm0PWne2szPnyn+bzF9as2aNAoGAs7S3t19wTAAAYOQKexGaMGGCMjIyQtZNmjRJbW1tkiS32y1Jg2ZkOjs7ndkbt9utvr4++f3+C2aOHTs26PGPHz8ekjn7cfx+v/r7+wfNFJ3hcrmUkJAQsgAAgNEp7EVo+vTp+vjjj0PWffLJJ7rpppskSWlpaXK73aqrq3O29/X1qb6+XtOmTZMkZWVlKTo6OiTT0dGh1tZWJ5OTk6NAIKC9e/c6mT179igQCIRkWltb1dHR4WRqa2vlcrmUlZUV5iMHAAAjTVS4d/hv//ZvmjZtmsrLy7VgwQLt3btXzz33nJ577jlJ33xUVVpaqvLyck2cOFETJ05UeXm5xowZo8LCQkmSZVlasmSJysrKNG7cOI0dO1YrV67U5MmTNXv2bEnfzDLNnTtXRUVFevbZZyVJ999/v/Lz85Weni5Jys3NVUZGhrxer5544gmdOHFCK1euVFFRETM9AAAg/EXo9ttvV3V1tdasWaPHHntMaWlp+vWvf617773XyTz00EPq7e1VcXGx/H6/srOzVVtbq/j4eCfz1FNPKSoqSgsWLFBvb69mzZqlzZs3KzIy0sls3bpVJSUlzrfLCgoKVFlZ6WyPjIzU9u3bVVxcrOnTpys2NlaFhYXasGFDuA8bAACMQGH/HaHRht8RGh34HSEAMMuw/Y4QAADASEERAgAAxqIIAQAAY1GEAACAsShCAADAWBQhAABgLIoQAAAwFkUIAAAYiyIEAACMRRECAADGoggBAABjUYQAAICxKEIAAMBYFCEAAGAsihAAADAWRQgAABiLIgQAAIxFEQIAAMaiCAEAAGNRhAAAgLEoQgAAwFgUIQAAYCyKEAAAMFbUcA8AuBJuXr19uIdw2Q6vmzfcQwCAUY8ZIQAAYCyKEAAAMBZFCAAAGIsiBAAAjEURAgAAxqIIAQAAY1GEAACAsShCAADAWBQhAABgLIoQAAAwFkUIAAAYiyIEAACMRRECAADGoggBAABjUYQAAICxKEIAAMBYFCEAAGAsihAAADAWRQgAABiLIgQAAIw15EWooqJCERERKi0tddbZtq21a9fK4/EoNjZWM2fO1IEDB0LuFwwGtXz5ciUlJSkuLk4FBQU6cuRISMbv98vr9cqyLFmWJa/Xq5MnT4Zk2traNH/+fMXFxSkpKUklJSXq6+sbqsMFAAAjyJAWoaamJj333HP6m7/5m5D169ev18aNG1VZWammpia53W7NmTNH3d3dTqa0tFTV1dWqqqpSQ0ODTp06pfz8fA0MDDiZwsJCtbS0qKamRjU1NWppaZHX63W2DwwMaN68eerp6VFDQ4Oqqqq0bds2lZWVDeVhAwCAEWLIitCpU6d077336vnnn1diYqKz3rZt/frXv9bDDz+su+++W5mZmXrxxRf15Zdf6uWXX5YkBQIBvfDCC3ryySc1e/Zs3XrrrXrppZe0f/9+vfHGG5KkgwcPqqamRv/5n/+pnJwc5eTk6Pnnn9f//M//6OOPP5Yk1dbW6sMPP9RLL72kW2+9VbNnz9aTTz6p559/Xl1dXUN16AAAYIQYsiK0bNkyzZs3T7Nnzw5Zf+jQIfl8PuXm5jrrXC6XZsyYoV27dkmSmpub1d/fH5LxeDzKzMx0Mrt375ZlWcrOznYyU6dOlWVZIZnMzEx5PB4nk5eXp2AwqObm5nOOOxgMqqurK2QBAACjU9RQ7LSqqkoffPCBmpqaBm3z+XySpJSUlJD1KSkp+vzzz51MTExMyEzSmcyZ+/t8PiUnJw/af3Jyckjm7MdJTExUTEyMkzlbRUWFHn300Us5TAAAMMKFfUaovb1d//qv/6qXXnpJ11577XlzERERIbdt2x607mxnZ86V/zaZv7RmzRoFAgFnaW9vv+CYAADAyBX2ItTc3KzOzk5lZWUpKipKUVFRqq+v13/8x38oKirKmaE5e0ams7PT2eZ2u9XX1ye/33/BzLFjxwY9/vHjx0MyZz+O3+9Xf3//oJmiM1wulxISEkIWAAAwOoW9CM2aNUv79+9XS0uLs0yZMkX33nuvWlpa9L3vfU9ut1t1dXXOffr6+lRfX69p06ZJkrKyshQdHR2S6ejoUGtrq5PJyclRIBDQ3r17ncyePXsUCARCMq2trero6HAytbW1crlcysrKCvehAwCAESbs1wjFx8crMzMzZF1cXJzGjRvnrC8tLVV5ebkmTpyoiRMnqry8XGPGjFFhYaEkybIsLVmyRGVlZRo3bpzGjh2rlStXavLkyc7F15MmTdLcuXNVVFSkZ599VpJ0//33Kz8/X+np6ZKk3NxcZWRkyOv16oknntCJEye0cuVKFRUVMdMDAACG5mLpi3nooYfU29ur4uJi+f1+ZWdnq7a2VvHx8U7mqaeeUlRUlBYsWKDe3l7NmjVLmzdvVmRkpJPZunWrSkpKnG+XFRQUqLKy0tkeGRmp7du3q7i4WNOnT1dsbKwKCwu1YcOGK3ewAADgqhVh27Y93IO4mnV1dcmyLAUCgSGZRbp59faw7xOjw+F184Z7CAAwYl3q+zd/awwAABiLIgQAAIxFEQIAAMaiCAEAAGNRhAAAgLEoQgAAwFgUIQAAYCyKEAAAMBZFCAAAGIsiBAAAjEURAgAAxqIIAQAAY1GEAACAsShCAADAWBQhAABgLIoQAAAwFkUIAAAYiyIEAACMRRECAADGoggBAABjUYQAAICxKEIAAMBYFCEAAGAsihAAADAWRQgAABiLIgQAAIxFEQIAAMaiCAEAAGNRhAAAgLEoQgAAwFgUIQAAYCyKEAAAMBZFCAAAGIsiBAAAjEURAgAAxqIIAQAAY1GEAACAsShCAADAWBQhAABgLIoQAAAwFkUIAAAYiyIEAACMRRECAADGoggBAABjhb0IVVRU6Pbbb1d8fLySk5N111136eOPPw7J2LattWvXyuPxKDY2VjNnztSBAwdCMsFgUMuXL1dSUpLi4uJUUFCgI0eOhGT8fr+8Xq8sy5JlWfJ6vTp58mRIpq2tTfPnz1dcXJySkpJUUlKivr6+cB82AAAYgcJehOrr67Vs2TI1Njaqrq5OX3/9tXJzc9XT0+Nk1q9fr40bN6qyslJNTU1yu92aM2eOuru7nUxpaamqq6tVVVWlhoYGnTp1Svn5+RoYGHAyhYWFamlpUU1NjWpqatTS0iKv1+tsHxgY0Lx589TT06OGhgZVVVVp27ZtKisrC/dhAwCAESjCtm17KB/g+PHjSk5OVn19vX7yk5/Itm15PB6VlpZq1apVkr6Z/UlJSdHjjz+upUuXKhAIaPz48dqyZYsWLlwoSTp69KhSU1O1Y8cO5eXl6eDBg8rIyFBjY6Oys7MlSY2NjcrJydFHH32k9PR07dy5U/n5+Wpvb5fH45EkVVVVafHixers7FRCQsJFx9/V1SXLshQIBC4pf7luXr097PvE6HB43bzhHgIAjFiX+v495NcIBQIBSdLYsWMlSYcOHZLP51Nubq6TcblcmjFjhnbt2iVJam5uVn9/f0jG4/EoMzPTyezevVuWZTklSJKmTp0qy7JCMpmZmU4JkqS8vDwFg0E1Nzefc7zBYFBdXV0hCwAAGJ2GtAjZtq0VK1bojjvuUGZmpiTJ5/NJklJSUkKyKSkpzjafz6eYmBglJiZeMJOcnDzoMZOTk0MyZz9OYmKiYmJinMzZKioqnGuOLMtSamrq5R42AAAYIYa0CD344IP64x//qP/6r/8atC0iIiLktm3bg9ad7ezMufLfJvOX1qxZo0Ag4Czt7e0XHBMAABi5hqwILV++XK+//rrefvtt3XDDDc56t9stSYNmZDo7O53ZG7fbrb6+Pvn9/gtmjh07Nuhxjx8/HpI5+3H8fr/6+/sHzRSd4XK5lJCQELIAAIDRKexFyLZtPfjgg3r11Vf11ltvKS0tLWR7Wlqa3G636urqnHV9fX2qr6/XtGnTJElZWVmKjo4OyXR0dKi1tdXJ5OTkKBAIaO/evU5mz549CgQCIZnW1lZ1dHQ4mdraWrlcLmVlZYX70AEAwAgTFe4dLlu2TC+//LL++7//W/Hx8c6MjGVZio2NVUREhEpLS1VeXq6JEydq4sSJKi8v15gxY1RYWOhklyxZorKyMo0bN05jx47VypUrNXnyZM2ePVuSNGnSJM2dO1dFRUV69tlnJUn333+/8vPzlZ6eLknKzc1VRkaGvF6vnnjiCZ04cUIrV65UUVERMz0AACD8X58/37U3v/vd77R48WJJ38waPfroo3r22Wfl9/uVnZ2t3/zmN84F1ZL01Vdf6Re/+IVefvll9fb2atasWXr66adDLl4+ceKESkpK9Prrr0uSCgoKVFlZqeuuu87JtLW1qbi4WG+99ZZiY2NVWFioDRs2yOVyXdLx8PV54NLxlX8AV4tLff8e8t8RGukoQsClowgBuFpcNb8jBAAAcLWiCAEAAGNRhAAAgLEoQgAAwFgUIQAAYCyKEAAAMBZFCAAAGIsiBAAAjEURAgAAxqIIAQAAY1GEAACAsShCAADAWBQhAABgLIoQAAAwFkUIAAAYiyIEAACMRRECAADGoggBAABjUYQAAICxKEIAAMBYFCEAAGAsihAAADAWRQgAABiLIgQAAIxFEQIAAMaiCAEAAGNRhAAAgLEoQgAAwFgUIQAAYCyKEAAAMBZFCAAAGCtquAcAYPS4efX24R7CZTu8bt5wDwHAMGJGCAAAGIsiBAAAjEURAgAAxqIIAQAAY1GEAACAsShCAADAWBQhAABgLIoQAAAwFkUIAAAYiyIEAACMRRECAADG4m+NATAafx8NMJsRM0JPP/200tLSdO211yorK0vvvffecA8JAABcBUZ9EXrllVdUWlqqhx9+WPv27dOPf/xj3XnnnWpraxvuoQEAgGEWYdu2PdyDGErZ2dm67bbb9MwzzzjrJk2apLvuuksVFRUXvX9XV5csy1IgEFBCQkLYxzcSp+UB4HLxcR6utEt9/x7V1wj19fWpublZq1evDlmfm5urXbt2nfM+wWBQwWDQuR0IBCR984QOhdPBL4dkvwBwNRmq/4YC53Pm37mLzfeM6iL05z//WQMDA0pJSQlZn5KSIp/Pd877VFRU6NFHHx20PjU1dUjGCAAmsH493COAqbq7u2VZ1nm3j+oidEZERETIbdu2B607Y82aNVqxYoVz+/Tp0zpx4oTGjRt33vt8W11dXUpNTVV7e/uQfOyGS8N5uDpwHq4enIurA+fhu7FtW93d3fJ4PBfMjeoilJSUpMjIyEGzP52dnYNmic5wuVxyuVwh66677rqhGqIkKSEhgX/JrwKch6sD5+Hqwbm4OnAevr0LzQSdMaq/NRYTE6OsrCzV1dWFrK+rq9O0adOGaVQAAOBqMapnhCRpxYoV8nq9mjJlinJycvTcc8+pra1NDzzwwHAPDQAADLNRX4QWLlyoL774Qo899pg6OjqUmZmpHTt26KabbhruocnlcumRRx4Z9FEcrizOw9WB83D14FxcHTgPV8ao/x0hAACA8xnV1wgBAABcCEUIAAAYiyIEAACMRRECAADGoggNk6efflppaWm69tprlZWVpffee2+4hzSqrV27VhERESGL2+12ttu2rbVr18rj8Sg2NlYzZ87UgQMHhnHEo8e7776r+fPny+PxKCIiQq+99lrI9kt57oPBoJYvX66kpCTFxcWpoKBAR44cuYJHMfJd7DwsXrx40Gtk6tSpIRnOw3dXUVGh22+/XfHx8UpOTtZdd92ljz/+OCTDa+LKoggNg1deeUWlpaV6+OGHtW/fPv34xz/WnXfeqba2tuEe2qj2gx/8QB0dHc6yf/9+Z9v69eu1ceNGVVZWqqmpSW63W3PmzFF3d/cwjnh06Onp0S233KLKyspzbr+U5760tFTV1dWqqqpSQ0ODTp06pfz8fA0MDFypwxjxLnYeJGnu3Lkhr5EdO3aEbOc8fHf19fVatmyZGhsbVVdXp6+//lq5ubnq6elxMrwmrjAbV9yPfvQj+4EHHghZ9/3vf99evXr1MI1o9HvkkUfsW2655ZzbTp8+bbvdbnvdunXOuq+++sq2LMv+7W9/e4VGaAZJdnV1tXP7Up77kydP2tHR0XZVVZWT+dOf/mRfc801dk1NzRUb+2hy9nmwbdtetGiR/Xd/93fnvQ/nYWh0dnbakuz6+nrbtnlNDAdmhK6wvr4+NTc3Kzc3N2R9bm6udu3aNUyjMsOnn34qj8ejtLQ03XPPPfrss88kSYcOHZLP5ws5Jy6XSzNmzOCcDLFLee6bm5vV398fkvF4PMrMzOT8hNk777yj5ORk/fVf/7WKiorU2dnpbOM8DI1AICBJGjt2rCReE8OBInSF/fnPf9bAwMCgP/qakpIy6I/DInyys7P1+9//Xn/4wx/0/PPPy+fzadq0afriiy+c551zcuVdynPv8/kUExOjxMTE82bw3d15553aunWr3nrrLT355JNqamrSz372MwWDQUmch6Fg27ZWrFihO+64Q5mZmZJ4TQyHUf8nNq5WERERIbdt2x60DuFz5513Ov88efJk5eTk6K/+6q/04osvOheEck6Gz7d57jk/4bVw4ULnnzMzMzVlyhTddNNN2r59u+6+++7z3o/z8O09+OCD+uMf/6iGhoZB23hNXDnMCF1hSUlJioyMHNTaOzs7B/0fAIZOXFycJk+erE8//dT59hjn5Mq7lOfe7Xarr69Pfr//vBmE34QJE3TTTTfp008/lcR5CLfly5fr9ddf19tvv60bbrjBWc9r4sqjCF1hMTExysrKUl1dXcj6uro6TZs2bZhGZZ5gMKiDBw9qwoQJSktLk9vtDjknfX19qq+v55wMsUt57rOyshQdHR2S6ejoUGtrK+dnCH3xxRdqb2/XhAkTJHEewsW2bT344IN69dVX9dZbbyktLS1kO6+JYTBsl2kbrKqqyo6OjrZfeOEF+8MPP7RLS0vtuLg4+/Dhw8M9tFGrrKzMfuedd+zPPvvMbmxstPPz8+34+HjnOV+3bp1tWZb96quv2vv377d//vOf2xMmTLC7urqGeeQjX3d3t71v3z573759tiR748aN9r59++zPP//ctu1Le+4feOAB+4YbbrDfeOMN+4MPPrB/9rOf2bfccov99ddfD9dhjTgXOg/d3d12WVmZvWvXLvvQoUP222+/befk5NjXX3895yHM/uVf/sW2LMt+55137I6ODmf58ssvnQyviSuLIjRMfvOb39g33XSTHRMTY992223OVycxNBYuXGhPmDDBjo6Otj0ej3333XfbBw4ccLafPn3afuSRR2y32227XC77Jz/5ib1///5hHPHo8fbbb9uSBi2LFi2ybfvSnvve3l77wQcftMeOHWvHxsba+fn5dltb2zAczch1ofPw5Zdf2rm5ufb48ePt6Oho+8Ybb7QXLVo06DnmPHx35zoHkuzf/e53TobXxJUVYdu2faVnoQAAAK4GXCMEAACMRRECAADGoggBAABjUYQAAICxKEIAAMBYFCEAAGAsihAAADAWRQgAABiLIgQAAIxFEQIAAMaiCAEAAGNRhAAAgLH+P7OjW3JyEpO9AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(y_test_re)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "bc0b9776-d618-4eab-92a2-f71cc6d6095e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Linear Regression MAPE on training set: 968611893.8950824\n",
      "Linear Regression MAE on training set: 0.42944108921576596\n",
      "Linear Regression MSE on training set: 0.3530610111465693\n",
      "Random Forest MAPE on training set: 384081718.01632553\n",
      "Random Forest MAE on training set: 0.3410029712046899\n",
      "Random Forest MSE on training set: 0.2490981985761064\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:    1.2s finished\n"
     ]
    }
   ],
   "source": [
    "# Check the metrics for predicting on the training set\n",
    "y_pred_lr = best_lr.predict(X_train)\n",
    "y_pred_rf = best_rf.predict(X_train)\n",
    "\n",
    "trainset_lr = trainset.copy()\n",
    "trainset_rf = trainset.copy()\n",
    "trainset_lr['Occupancy'] = y_pred_lr\n",
    "trainset_rf['Occupancy'] = y_pred_rf\n",
    "\n",
    "# Reverse scale the y\n",
    "# all scores looks quite differently if reverse scaling is not applied\n",
    "trainset_re = scaler.inverse_transform(trainset[columns_to_normalize])\n",
    "trainset_lr_re = scaler.inverse_transform(trainset_lr[columns_to_normalize])\n",
    "trainset_rf_re = scaler.inverse_transform(trainset_rf[columns_to_normalize])\n",
    "\n",
    "y_train_re = trainset_re[:, 3]\n",
    "y_pred_lr_re = trainset_lr_re[:, 3]\n",
    "y_pred_rf_re = trainset_rf_re[:, 3]\n",
    "\n",
    "# Calculate MAPE with a small denominator to prevent inf\n",
    "mape_lr = np.mean(np.abs((y_train_re - y_pred_lr_re) / (y_train_re + 1e-10))) * 100\n",
    "mape_rf = np.mean(np.abs((y_train_re - y_pred_rf_re) / (y_train_re + 1e-10))) * 100\n",
    "\n",
    "# Calculate MAE \n",
    "mae_lr = mean_absolute_error(y_train_re, y_pred_lr_re)\n",
    "mae_rf = mean_absolute_error(y_train_re, y_pred_rf_re)\n",
    "\n",
    "# Calculate MSE\n",
    "mse_lr = mean_squared_error(y_train_re, y_pred_lr_re)\n",
    "mse_rf = mean_squared_error(y_train_re, y_pred_rf_re)\n",
    "\n",
    "print(\"Linear Regression MAPE on training set:\", mape_lr)\n",
    "print(\"Linear Regression MAE on training set:\", mae_lr)\n",
    "print(\"Linear Regression MSE on training set:\", mse_lr)\n",
    "\n",
    "print(\"Random Forest MAPE on training set:\", mape_rf)\n",
    "print(\"Random Forest MAE on training set:\", mae_rf)\n",
    "print(\"Random Forest MSE on training set:\", mse_rf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb3789d5-9a2a-42a6-bbed-b285e8a811cd",
   "metadata": {},
   "source": [
    "## Training model: SVR with some GridSearchCV for best parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "d9b772c4-e8f3-4e7b-9248-33f87cdb6ecf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1cAAAIhCAYAAACizkCYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACN9UlEQVR4nOzdeVhUZf8G8HsGhhn2nQGVHQXGFcEFzF1QM7P6lWbl7ptmm1pvyVtm2mKrW+XWotmiZmnmUoK7JrngQoqoIIoLOwICgizP7w9icgR1BgdmGO7Pdc2Vc+bMme/MM4e4Oed8H4kQQoCIiIiIiIjui9TQBRAREREREZkChisiIiIiIiI9YLgiIiIiIiLSA4YrIiIiIiIiPWC4IiIiIiIi0gOGKyIiIiIiIj1guCIiIiIiItIDhisiIiIiIiI9YLgiIiIiIiLSA4YroiZo5cqVkEgkd7zt3r27wV/7woULjfpcfZBIJHj77bfv+PjChQshkUjwxx9/3HGdL7/8EhKJBOvXr9dLTT4+Phg7dqxetmVIffr00fgOKhQKqFQqvPvuu7h582a9tpmYmIi33367zu/L2LFj4ePjc39FG1jN/qBQKHDx4sVaj/fp0wft2rUzQGXA7t27IZFI8PPPPxvk9XV14cIFDBkyBE5OTpBIJJg6depd1y8rK8Pnn3+OBx54AI6OjrCwsEDLli0xfPhw7Nmzp3GKNnHvv/8+fv31V0OXQdTozA1dABHV34oVKxAUFFRruUqlMkA19zZkyBDExcXBw8PD0KXU6ZlnnsHrr7+Ob775BoMGDapznRUrVsDV1RVDhw7Vy2tu2LABdnZ2etmWofn5+eGHH34AAGRnZ+Orr77CzJkzkZaWhuXLl+u8vcTERMyePRt9+vSpFaRmzpyJl19+WR9lG1xZWRnefPNNfPfdd4YupcmaNm0aDh48iG+++Qbu7u53/RmTk5ODQYMGISEhAePHj8d///tfODk54cqVK9i4cSP69++P+Ph4dOzYsRHfgel5//338fjjj+ORRx4xdClEjYrhiqgJa9euHcLCwgxdhtZcXV3h6upq6DLuyNnZGcOGDcOvv/6K3NxcODs7azyelJSEuLg4vPLKK5DJZPf1Wjdu3IClpSVCQkLuazvGxNLSEt27d1ffHzx4MFQqFb799lssWrQICoVCb6/l7++vt20Z2qBBg/Djjz/i1VdfbXa/0N+4cQMKhQISieS+tnPy5El07dpVq1/kR48ejRMnTmDbtm3o16+fxmNPPvkkpk+fDkdHx/uqh4iaL54WSGTC1qxZA4lEgs8//1xj+axZs2BmZobY2FgA1afUSCQSfPTRR3jvvffg5eUFhUKBsLAw7Nix456vExsbi2HDhqFVq1ZQKBQICAjApEmTkJOTo7FeXacF1pz6dPjwYfTs2RNWVlbw8/PDBx98gKqqKo3nFxYW4tVXX4Wvr6/6NJ6pU6eiuLi41nr/+c9/4OzsDBsbGwwaNAhnz57V6jObMGECbt68iR9//LHWYytWrAAAjB8/HgAwe/ZsdOvWDU5OTrCzs0Pnzp3x9ddfQwih8TwfHx889NBDWL9+PUJCQqBQKDB79mz1Y7eeFlhaWopXXnkFnTp1gr29PZycnBAeHo6NGzfWqkcikeCFF17Ad999h+DgYFhZWaFjx47YvHlzrXWTkpIwcuRIKJVKyOVyeHl5YfTo0SgrK1Ovk5GRgUmTJqFVq1awsLCAr68vZs+ejYqKCq0+u9uZm5ujU6dOuHnzJvLz89XLjxw5gieffBI+Pj6wtLSEj48PRo4cqXFq3MqVK/HEE08AAPr27as+3XDlypUA6j4tsLS0FNHR0Rrfj+eff17jteuyYMECSCQSJCcn13rs9ddfh4WFhfq7fOzYMTz00ENwc3ODXC5HixYtMGTIEFy+fFn3D+gfr732GpydnfH666/fdb2a/bTmM7jV7ae8vv3225BIJEhISMATTzyh/i5Nnz4dFRUVOHPmDAYNGgRbW1v4+Pjgo48+qvM1S0tLMX36dLi7u8PS0hK9e/fGsWPHaq135MgRPPzww3BycoJCoUBISAh++uknjXVq9v+YmBiMHz8erq6usLKy0vgO3i4tLQ3PPPOM+vMODg7Gp59+qv7ZUHP6YnJyMn7//Xf19+ROpx7Hx8fj999/x4QJE2oFqxpdunSBl5eX+v7JkycxbNgwODo6QqFQoFOnTvj22281nlNTx48//ojXX38dHh4esLGxwdChQ5GZmYnr16/j2WefhYuLC1xcXDBu3DgUFRVpbKNmf162bBnatGkDuVwOlUqFNWvW1KpRl5pWr16NN954Ay1atICdnR0GDBiAM2fO1Nrm9u3b0b9/f9jZ2cHKygo9evSo9fO/5nt16tQpjBw5Evb29lAqlRg/fjwKCgo03ktxcTG+/fZb9Zj06dMHAFBSUqL+Oa5QKODk5ISwsDCsXr26zvEganIEETU5K1asEADEX3/9JcrLyzVuFRUVGutOnjxZWFhYiMOHDwshhNixY4eQSqXizTffVK+TmpoqAAhPT0/xwAMPiF9++UWsW7dOdOnSRchkMnHgwIFar52amqpetmTJEjF37lzx22+/iT179ohvv/1WdOzYUQQGBoqbN2/e9bm9e/cWzs7OonXr1mLp0qUiNjZWTJkyRQAQ3377rXq94uJi0alTJ+Hi4iLmzZsntm/fLhYuXCjs7e1Fv379RFVVlRBCiKqqKtG3b18hl8vFe++9J2JiYsSsWbOEn5+fACBmzZp118+2srJSeHt7i06dOmksr6ioEB4eHqJ79+7qZWPHjhVff/21iI2NFbGxseKdd94RlpaWYvbs2RrP9fb2Fh4eHsLPz0988803YteuXeLQoUPqx8aMGaNeNz8/X4wdO1Z89913YufOneKPP/4Qr776qpBKpRqfhxBCABA+Pj6ia9eu4qeffhJbt24Vffr0Eebm5iIlJUW93vHjx4WNjY3w8fERS5cuFTt27BDff/+9GD58uCgsLBRCCJGeni48PT2Ft7e3WLZsmdi+fbt45513hFwuF2PHjr3rZyZE9Ti2bdu21vKwsDDh4OCg8b1ct26deOutt8SGDRvEnj17xJo1a0Tv3r2Fq6uryM7OFkIIkZWVJd5//30BQHzxxRciLi5OxMXFiaysLCGEEGPGjBHe3t7qbVZVVYmBAwcKc3NzMXPmTBETEyM++eQTYW1tLUJCQkRpaekda8/OzhYWFhbijTfe0FheUVEhWrRoIR577DEhhBBFRUXC2dlZhIWFiZ9++kns2bNHrF27VkyePFkkJibe8zO6Xc3+cPjwYbFw4UIBQOzYseOOn2nNfrpixYpa27r9uz1r1iwBQAQGBop33nlHxMbGitdee00AEC+88IIICgoSixYtErGxsWLcuHECgPjll1/Uz9+1a5f6Z8KwYcPEpk2bxPfffy8CAgKEnZ2dxvdr586dwsLCQvTs2VOsXbtW/PHHH2Ls2LG1aq15vy1bthTPPvus+P3338XPP/9c62dWjaysLNGyZUvh6uoqli5dKv744w/xwgsvCADiueeeE0IIUVBQIOLi4oS7u7vo0aOH+ntyp/Gu+U79/vvvdx2bGklJScLW1lb4+/uLVatWiS1btoiRI0cKAOLDDz+s9Xl5e3uLsWPHij/++EMsXbpU2NjYiL59+4rIyEjx6quvipiYGPHhhx8KMzMz8eKLL9YaQ09PT6FSqcTq1avFb7/9JgYNGiQAiHXr1tW7Jh8fH/H000+LLVu2iNWrVwsvLy/RunVrjc/9u+++ExKJRDzyyCNi/fr1YtOmTeKhhx4SZmZmYvv27er1bv1evfXWWyI2NlbMmzdPyOVyMW7cOPV6cXFxwtLSUjz44IPqMTl16pQQQohJkyYJKysrMW/ePLFr1y6xefNm8cEHH4jPPvtMqzEhMnYMV0RNUM0vKXXdzMzMNNYtLS0VISEhwtfXVyQmJgqlUil69+6t8T/Wml/aWrRoIW7cuKFeXlhYKJycnMSAAQNqvfatAelWVVVVory8XFy8eFEAEBs3brzrc3v37i0AiIMHD2psR6VSiYEDB6rvz507V0ilUnVIrPHzzz8LAGLr1q1CCCF+//13AUAsXLhQY7333ntPq3AlxL+/QBw9elS9bNOmTQKA+PLLL+t8TmVlpSgvLxdz5swRzs7O6rAnRHWAMjMzE2fOnKn1vNvD1e0qKipEeXm5mDBhgggJCdF4DIBQKpXqgCSEEBkZGUIqlYq5c+eql/Xr1084ODiog0ldJk2aJGxsbMTFixc1ln/yyScCgPoXozupCQI1IT89PV289dZbAoBYunTpXZ9bUVEhioqKhLW1tca4rVu3TgAQu3btqvWc28PVH3/8IQCIjz76SGO9tWvXCgBi+fLld63hscceE61atRKVlZXqZVu3bhUAxKZNm4QQQhw5ckQAEL/++utdt6WtW8NVWVmZ8PPzE2FhYervjj7C1aeffqqxXqdOnQQAsX79evWy8vJy4erqqg6RQvz7i3nnzp01vssXLlwQMplMTJw4Ub0sKChIhISEiPLyco3Xeuihh4SHh4f6M615v6NHj9bq85kxY0adPxuee+45IZFINPYnb29vMWTIkHtuc/LkyQKASEpK0qqGJ598UsjlcpGWlqaxfPDgwcLKykrk5+cLIf79vIYOHaqx3tSpUwUA8dJLL2ksf+SRR4STk5PGMgDC0tJSZGRkqJdVVFSIoKAgERAQUO+aHnzwQY31fvrpJwFAxMXFCSGq/3Dl5ORUq/bKykrRsWNH0bVrV/Wymu/V7fvZlClThEKh0PiuWFtb1/mzrV27duKRRx6ptZzIVPC0QKImbNWqVTh8+LDG7eDBgxrryOVy/PTTT8jNzUXnzp0hhMDq1athZmZWa3uPPfaYxnUxtra2GDp0KPbu3YvKyso71pGVlYXJkyfD09MT5ubmkMlk8Pb2BgCcPn36nu/D3d0dXbt21VjWoUMHjdPENm/ejHbt2qFTp06oqKhQ3wYOHKjRIXHXrl0AgKefflpje0899dQ966gxbtw4SKVSfPPNN+plK1asgLW1NUaMGKFetnPnTgwYMAD29vYwMzODTCbDW2+9hdzcXGRlZdV6P23atNHq9detW4cePXrAxsZG/Xl+/fXXdX6Wffv2ha2trfq+UqmEm5ub+rMrKSnBnj17MHz48Lte77Z582b07dsXLVq00Ph8Bw8eDABadVA7deoUZDIZZDIZPDw8MGfOHERHR2PSpEka6xUVFeH1119HQEAAzM3NYW5uDhsbGxQXF2v1fanLzp07AaBW58UnnngC1tbW9zy9ddy4cbh8+TK2b9+uXrZixQq4u7urP4OAgAA4Ojri9ddfx9KlS5GYmFivWutiYWGBd999F0eOHKl1Ot39eOihhzTuBwcHQyKRqN8TUH36ZkBAQJ0dC5966imN66G8vb0RERGh3s+Sk5ORlJSk3t9u/e48+OCDSE9Pr3UK2v/93/9pVfvOnTuhUqlq/WwYO3YshBDqMW9IO3fuRP/+/eHp6VmrhpKSEsTFxWksr+vzBqqb+dy+PC8vr9apgf3794dSqVTfNzMzw4gRI5CcnKw+9VTXmh5++GGN+x06dAAA9XgfOHAAeXl5GDNmjMb4VVVVYdCgQTh8+HCtU6/r2mZpaWmtn3t16dq1K37//XfMmDEDu3fvxo0bN+75HKKmhOGKqAkLDg5GWFiYxi00NLTWegEBAejZsydKS0vx9NNP37GTlru7e53Lbt68WeuXgBpVVVWIiorC+vXr8dprr2HHjh04dOgQ/vrrLwDQ6n+ctzeOAKpD4a3PzczMREJCgvqX95qbra0thBDqa2Jyc3Nhbm5ea5t1vbc78fb2Rv/+/fHjjz+irKwMOTk52Lx5M5544gl1kDl06BCioqIAVLdn//PPP3H48GG88cYbdb5vbTskrl+/HsOHD0fLli3x/fffIy4uDocPH8b48eNRWlpaa/17fXbXrl1DZWUlWrVqddfXzczMxKZNm2p9vm3btgWAWtfP1cXf3x+HDx/GoUOHsG7dOnTs2BFz586tdc3IU089hc8//xwTJ07Etm3bcOjQIRw+fBiurq71/kWrZtxvD5ASiQTu7u7Izc296/MHDx4MDw8P9XV1165dw2+//YbRo0er/xBhb2+PPXv2oFOnTvjf//6Htm3bokWLFpg1axbKy8vrVfetnnzySXTu3BlvvPGGXrYHAE5OThr3LSwsYGVlVau5iIWFRZ3frzv9TKj5PDMzMwEAr776aq3vzpQpUwDU/u5ouy/k5ubWuW6LFi3Uj+uq5lqq1NTUBqmhrs/7bstv/8zv9Hnf+lq61nT7zwi5XA7g359RNWP4+OOP1xrDDz/8EEII5OXl6bTNu1m0aBFef/11/Prrr+jbty+cnJzwyCOP4Ny5c/d8LlFTwG6BRM3AV199hS1btqBr1674/PPPMWLECHTr1q3WehkZGXUus7CwgI2NTZ3bPnnyJE6cOIGVK1dizJgx6uV1NQe4Hy4uLrC0tNQ4mnT740D1//QrKipqdfur673dzYQJExAbG4uNGzfi6tWruHnzJiZMmKB+fM2aNZDJZNi8ebPGL6p3mtdF225o33//PXx9fbF27VqN59ztov+7cXJygpmZ2T0bLri4uKBDhw5477336ny85he3u6lpggJUNwXo27cv2rZti6lTp+Khhx6CjY0NCgoKsHnzZsyaNQszZsxQP7esrKzWL3C6qBn37OxsjYAlhEBGRga6dOly1+ebmZlh1KhRWLRoEfLz89XBety4cRrrtW/fHmvWrIEQAgkJCVi5ciXmzJkDS0tLjfdTHxKJBB9++CEiIyPrbF1f8z27/btQn5ChrTv9TKjZt2r2u+joaDz22GN1biMwMFDjvrb7grOzM9LT02stv3r1qsZr62LgwIH43//+h19//fWO0y00dA13c6fPu6aWhqipZv3PPvtMo9vnrW49mna/rK2tMXv2bMyePRuZmZnqo1hDhw5FUlKS3l6HyFB45IrIxP3999946aWXMHr0aOzbtw8dOnTAiBEjcO3atVrrrl+/XuMvqdevX8emTZvQs2fPOk8jBP79RanmL5c1li1bpsd3UX26TUpKCpydnWsdrQsLC1N3juvbty8AqOdbqlFX97+7eeSRR+Ds7IxvvvkGK1asQJs2bfDAAw+oH5dIJDA3N9f4XG7cuHHfcxVJJBJYWFho/AKakZFRZ7dAbdR0eFu3bt1djz499NBDOHnyJPz9/ev8fLUJV7dzdnbGBx98gMzMTHz22WcAqt+fEKLW9+Wrr76qdeqpLn8N79+/P4DqcHqrX375BcXFxerH72bcuHEoLS3F6tWrsXLlSoSHh9c5j1zN++jYsSPmz58PBwcHHD169J7b18aAAQMQGRmJOXPm1DparFQqoVAokJCQoLG8vt8NbaxevVqj++XFixdx4MABdee3wMBAtG7dGidOnKjzexMWFqZx2qou+vfvj8TExFqf7apVqyCRSNT7ui46d+6MwYMH4+uvv77jaYVHjhxBWlqauoadO3eqg8utNVhZWd0xjNTXjh071EeSAKCyshJr166Fv7+/+uizvmvq0aMHHBwckJiYeMcxrDnSpovbzz6oi1KpxNixYzFy5EicOXMGJSUlOr8OkbHhkSuiJuzkyZN1tsn29/eHq6sriouLMXz4cPj6+mLx4sWwsLDATz/9hM6dO2PcuHG1jrKYmZkhMjIS06dPR1VVFT788EMUFhaq24bXJSgoCP7+/pgxYwaEEHBycsKmTZvUbd71ZerUqfjll1/Qq1cvTJs2DR06dEBVVRXS0tIQExODV155Bd26dUNUVBR69eqF1157DcXFxQgLC8Off/6pc+iRy+V4+umn8dlnn0EIgQ8++EDj8SFDhmDevHl46qmn8OyzzyI3NxeffPJJrdCgq5qW7VOmTMHjjz+OS5cu4Z133oGHh0e9T5uZN28eHnjgAXTr1g0zZsxAQEAAMjMz8dtvv2HZsmWwtbXFnDlzEBsbi4iICLz00ksIDAxEaWkpLly4gK1bt2Lp0qX3PLWwLqNHj8a8efPwySef4Pnnn4ednR169eqFjz/+GC4uLvDx8cGePXvw9ddfw8HBQeO57dq1AwAsX74ctra2UCgU8PX1rfNUyMjISAwcOBCvv/46CgsL0aNHDyQkJGDWrFkICQnBqFGj7llrUFAQwsPDMXfuXFy6dKnW0aPNmzdj8eLFeOSRR+Dn5wchBNavX4/8/HxERkaq1+vfvz/27NlT7xb2H374IUJDQ5GVlaU+LROoDnTPPPMMvvnmG/j7+6Njx444dOiQzn840EVWVhYeffRR/Oc//0FBQQFmzZoFhUKB6Oho9TrLli3D4MGDMXDgQIwdOxYtW7ZEXl4eTp8+jaNHj2LdunX1eu1p06Zh1apVGDJkCObMmQNvb29s2bIFixcvxnPPPaf1NYy3W7VqFQYNGoTBgwdj/PjxGDx4MBwdHZGeno5NmzZh9erViI+Ph5eXF2bNmqW+HvGtt96Ck5MTfvjhB2zZsgUfffQR7O3t61XDnbi4uKBfv36YOXMmrK2tsXjxYiQlJWmcWqvvmmxsbPDZZ59hzJgxyMvLw+OPPw43NzdkZ2fjxIkTyM7OxpIlS3R+L+3bt8fu3buxadMmeHh4wNbWFoGBgejWrRseeughdOjQAY6Ojjh9+jS+++47hIeHw8rKSufXITI6BmqkQUT34W7dAnFLR7tnnnlGWFlZ1er0VtOFbf78+UKIf7uQffjhh2L27NmiVatWwsLCQoSEhIht27bV+dq3dvxLTEwUkZGRwtbWVjg6OoonnnhCpKWl1epgdqdugXW18L69G5wQ1a2w33zzTREYGCgsLCyEvb29aN++vZg2bZpGh638/Hwxfvx44eDgIKysrERkZKRISkrSultgjRMnTqg7MF69erXW4998840IDAwUcrlc+Pn5iblz54qvv/661nu8WyezuroFfvDBB8LHx0fI5XIRHBwsvvzyS3WXrlsBEM8//7xW20xMTBRPPPGEcHZ2FhYWFsLLy0uMHTtWo2V1dna2eOmll4Svr6+QyWTCyclJhIaGijfeeEMUFRXd9bO60zgKIcSWLVsEAHWL+suXL4v/+7//E46OjsLW1lYMGjRInDx5ss66FyxYIHx9fYWZmZlGp7y6vh83btwQr7/+uvD29hYymUx4eHiI5557Tly7du2utd9q+fLl6q5tBQUFGo8lJSWJkSNHCn9/f2FpaSns7e1F165dxcqVK2t9Ftr87/XWboG3e+qppwSAWp9pQUGBmDhxolAqlcLa2loMHTpUXLhw4Y7dAmta29cYM2aMsLa2rvV6t49fTae57777Trz00kvC1dVVyOVy0bNnT3HkyJFazz9x4oQYPny4cHNzEzKZTLi7u4t+/fppdIq82/u9k4sXL4qnnnpKODs7C5lMJgIDA8XHH3+s0dVRCO27Bda4ceOGWLRokQgPDxd2dnbC3Nxc3XZ/y5YtGuv+/fffYujQocLe3l5YWFiIjh071urYWPN53doy/W7vua7xqdmfFy9eLPz9/YVMJhNBQUHihx9+qFX//dR0p66Te/bsEUOGDBFOTk5CJpOJli1biiFDhmg8/07fq7p+th8/flz06NFDWFlZCQCid+/eQojqLpBhYWHC0dFR/bNz2rRpIicnp9b7JGqKJELcNtslETU7Fy5cgK+vLz7++GO8+uqrhi6HiKjZkUgkeP7552tN+k5ETQuvuSIiIiIiItIDhisiIiIiIiI94GmBREREREREesAjV0RERERERHrAcEVERERERKQHDFdERERERER6wEmE61BVVYWrV6/C1tYWEonE0OUQEREREZGBCCFw/fp1tGjRAlLp3Y9NMVzV4erVq/D09DR0GUREREREZCQuXbqEVq1a3XUdhqs62NraAqj+AO3s7AxcDVBeXo6YmBhERUVBJpMZuhzSA46p6eGYmiaOq+nhmJomjqvpMaYxLSwshKenpzoj3A3DVR1qTgW0s7MzmnBlZWUFOzs7g3+5SD84pqaHY2qaOK6mh2NqmjiupscYx1Sby4XY0IKIiIiIiEgPGK6IiIiIiIj0gOGKiIiIiIhIDxiuiIiIiIiI9IDhioiIiIiISA8YroiIiIiIiPSA4YqIiIiIiEgPGK6IiIiIiIj0gOGKiIiIiIhIDxiuiIiIiIiI9IDhioiIiIiISA8YroiIiIiIiPSA4YqIiIiIiEgPGK6MXGWVwMHUPMTnSHAwNQ+VVcLQJRERERERUR3MDV0A3dkfJ9Mxe1Mi0gtKAZhh1bkj8LBXYNZQFQa18zB0eUREREREdAseuTJSf5xMx3PfH/0nWP0ro6AUz31/FH+cTDdQZUREREREVBeGKyNUWSUwe1Mi6joBsGbZ7E2JPEWQiIiIiMiIMFwZoUOpebWOWN1KAEgvKMWh1LzGK4qIiIiIiO6K4coIZV2/c7C61co/U3HySgGE4BEsIiIiIiJDY0MLI+Rmq9BqvW2JmdiWmIkW9gpEqpSIauuOrr5OkJkxMxMRERERNTaGKyPU1dcJHvYKZBSU1nndlQSAvaUM3XydsPdcDq4WlOLbuIv4Nu4i7BTm6Bfkhqi27ujVxhU2cg4xEREREVFj4G/eRshMKsGsoSo89/1RSACNgCX5578f/F97DGrngdLySvyZnIOYU5nYfjoTucU38evxq/j1+FVYmEnRI8AZkSp3DFC5aX1EjIiIiIiIdMdwZaQGtfPAkmc63zLPVTX32+a5UsjM0D9Yif7BSlRWCRxLu4aYxEzEnMrAhdwS7DqTjV1nsvHGr0AnTwdEqdwR1VYJf1cbA70zIiIiIiLTxHBlxAa180Ckyh1xyVmI2XcQUT27ITzADWZSSZ3rm0klCPNxQpiPE6IHByE5q6g6aCVm4sSlfBxLq759+EcS/Fyt1UGrUysHSO+wTSIiIiIi0g7DlZEzk0rQzdcJuacFuvk63TFY3U4ikaC10hatlbZ4vm8AMgpKEXs6E7GJmYhLycH57GIs3ZOCpXtS4Gorx4BgJaLaKhHh7wy5uVkDvysiIiIiItPDcNVMuNsrMKq7N0Z190ZhaTn2nMlGTGImdidlIft6GVYfSsPqQ2mwtjBDn0A3RKqU6BvoBnsrmaFLJyIiIiJqEhiumiE7hQxDO7bA0I4tcLOiCn+dz0VMYgZiEzORWViGLX+nY8vf6TCXStDdzxmRKiUiVUq0cLA0dOlEREREREbL4BMiLV68GL6+vlAoFAgNDcW+ffvuuO7+/fvRo0cPODs7w9LSEkFBQZg/f77GOitXroREIql1Ky3VbmLe5sbCXIpebVzx7iPtETejPzY+3wPP9/VHG6UNKqoE9ifnYNZvpxDxwU489Nk+LNpxDkkZhZy4mIiIiIjoNgY9crV27VpMnToVixcvRo8ePbBs2TIMHjwYiYmJ8PLyqrW+tbU1XnjhBXTo0AHW1tbYv38/Jk2aBGtrazz77LPq9ezs7HDmzBmN5yoUbEN+L1KpBB09HdDR0wH/HRiECznFiE3MRExiBo5cvIaTVwpx8koh5sWehaeTJaJU7ohUKRHm7QhzTlxMRERERM2cQcPVvHnzMGHCBEycOBEAsGDBAmzbtg1LlizB3Llza60fEhKCkJAQ9X0fHx+sX78e+/bt0whXEokE7u7uDf8GTJyPizX+08sP/+nlh5yiMuw8nYWYxAzsO5eDS3k38PX+VHy9PxWOVjL0D64+dbBXa1dYWrAhBhERERE1PwYLVzdv3kR8fDxmzJihsTwqKgoHDhzQahvHjh3DgQMH8O6772osLyoqgre3NyorK9GpUye88847GqHsdmVlZSgrK1PfLywsBACUl5ejvLxc27fUYGpqMGQt9nIpHu3kjkc7uaPkZgX2J+di++ks7DqTg2sl5fg5/jJ+jr8MhUyKB/yd0T/YDX0DXeFsbWGwmo2ZMYwp6RfH1DRxXE0Px9Q0cVxNjzGNqS41SISBLp65evUqWrZsiT///BMRERHq5e+//z6+/fbbWqf13apVq1bIzs5GRUUF3n77bcycOVP92F9//YXk5GS0b98ehYWFWLhwIbZu3YoTJ06gdevWdW7v7bffxuzZs2st//HHH2FlZXUf79L0VQrgfKEEf+dJ8Pc1CfLK/m0VL4GAny3Q3qkK7Z0EXHhmJhERERE1MSUlJXjqqadQUFAAOzu7u65r8HB14MABhIeHq5e/9957+O6775CUlHTH56ampqKoqAh//fUXZsyYgc8//xwjR46sc92qqip07twZvXr1wqJFi+pcp64jV56ensjJybnnB9gYysvLERsbi8jISMhkxtsaXQiBpIwibD+dhe1JWUhMv67xeBs3GwwIdsOAYFe0a2EHiaT5TlzcVMaUtMcxNU0cV9PDMTVNHFfTY0xjWlhYCBcXF63ClcFOC3RxcYGZmRkyMjI0lmdlZUGpVN71ub6+vgCA9u3bIzMzE2+//fYdw5VUKkWXLl1w7ty5O25PLpdDLpfXWi6TyQw+mLcytnrq0sHLCR28nDB9YBAuXyvB9sRMxCRm4mBqHs5mFeFsVhEW7zkPD3uFusV7N19nWJg3z4YYTWFMSTccU9PEcTU9HFPTxHE1PcYwprq8vsHClYWFBUJDQxEbG4tHH31UvTw2NhbDhg3TejtCCI2jTnU9fvz4cbRv3/6+6iXdtXK0wtgevhjbwxf5JTex60wWYk5lYs/ZbKQXlGJV3EWsirsIW4U5+ga6IaqtEr3buMJWwR+KRERERNT0GLRb4PTp0zFq1CiEhYUhPDwcy5cvR1paGiZPngwAiI6OxpUrV7Bq1SoAwBdffAEvLy8EBQUBqJ736pNPPsGLL76o3ubs2bPRvXt3tG7dGoWFhVi0aBGOHz+OL774ovHfIKk5WFng0ZBWeDSkFUrLK3EgJQcxpzKx/XQmcopu4rcTV/HbiauwMJMi3N8ZUW2ViAxWws2OF2oRERERUdNg0HA1YsQI5ObmYs6cOUhPT0e7du2wdetWeHt7AwDS09ORlpamXr+qqgrR0dFITU2Fubk5/P398cEHH2DSpEnqdfLz8/Hss88iIyMD9vb2CAkJwd69e9G1a9dGf39UN4XMDP2ClOgXpERllcDxS9cQk5iJmFOZSM0pxp6z2dhzNhtvbDiJTp4OiGqrRJTKHQFuNoYunYiIiIjojgwargBgypQpmDJlSp2PrVy5UuP+iy++qHGUqi7z58/H/Pnz9VUeNTAzqQSh3k4I9XbCjEFBSMkuUget45fy1beP/jgDPxdrRLZVIkqlRIinI6TS5tsQg4iIiIiMj8HDFVENiUSCADdbBLjZYkqfAGQWlmL76eqgdSAlB+dzirFsz3ks23MeLjZyDAiuvk4rwt8FChknLiYiIiIiw2K4IqOltFPg6W7eeLqbN66XlmPP2WzEnMrErjNZyCkqw5rDl7Dm8CVYWZihdxtXRLVVol+gEvZWbIhBRERERI2P4YqaBFuFDA91aIGHOrTAzYoqHEzNRcypTMQmZiKjsBS/n8zA7yczYCaVoJuvE6JUSkS2dUdLB0tDl05EREREzQTDFTU5FuZS9Gztip6tXTFnWFv8faUAsf9cp3Um8zoOpOTiQEou3t6UiLYt7BClckekSolgD9tmPXExERERETUshitq0iQSCTq0ckCHVg54JSoQF3OL1UHryMU8nLpaiFNXCzF/+1m0crREpKq682AXH0eYmzXPiYuJiIiIqGEwXJFJ8Xa2xsSefpjY0w+5RWXYkVQ9cfG+c9m4fO0GVvx5ASv+vAAHKxn6BbkhSuWOXm1cYGXBXYGIiIiI7g9/oyST5Wwjx/AwTwwP80TJzQrsO1c9cfHOpExcKynH+qNXsP7oFcjNpejZ2gVRKnf0C3aDi43c0KUTERERURPEcEXNgpWFOQa2dcfAtu6oqKzCkYvXqk8fTMzApbwb2H46C9tPZ0EiAcK8HdXXafm4WBu6dCIiIiJqIhiuqNkxN5Oiu58zuvs5480hwUjKuK4OWievFOLwhWs4fOEa3tt6Gm2UNurrtNq3tOfExURERER0RwxX1KxJJBIEe9gh2MMOL/VvjSv5N7D9n6B18HwezmYW4WxmEb7YlQJ3OwUGqKqv0+ru5wwLczbEICIiIqJ/MVwR3aKlgyXGRPhgTIQPCkrKsetMFmITM7H7TBYyCkvx/V9p+P6vNNjKzdEnyA1RKiX6BLrCVsGJi4mIiIiaO4Yrojuwt5LhkZCWeCSkJUrLKxGXkouYxOqJi3OKyrDpxFVsOnEVMjMJwv1dEKlSIjJYCXd7haFLJyIiIiIDYLgi0oJCZoa+QW7oG+SG9x5ph2OX8tXXaZ3PLsbes9nYezYbM389iY6eDohSKRGlUiLAzYYTFxMRERE1EwxXRDqSSiUI9XZEqLcjZgwOQnJWkTpoHUvLx4lL1bePt52Br4v1Pw0xlAjxcoQZG2IQERERmSyGK6L7FOBmgwA3GzzXxx9ZhaXYfjoLMYkZOJCci9ScYizfex7L956Hs7UFBgQrEalSopuPvaHLJiIiIiI9Y7gi0iM3OwWe6uaFp7p5oaisAnvOZCMmMQM7k7KQW3wTa49cwtojl2Apk6K1rRRlHlcR2dYDjtYWhi6diIiIiO4TwxVRA7GRm2NIBw8M6eCB8soqHDyfh9jEDMQkZiK9oBQJeVK8tv4kzH49ha4+TtUNMVRKeDpZGbp0IiIiIqoHhiuiRiAzk+KB1i54oLUL3n64LU6k5WHJbwdwocIOZzKLEHc+F3HnczFncyKCPeyqG2K0VULlYceGGERERERNBMMVUSOTSCRo28IOD3pV4cEHI5BeWI6Yf45oHbmQh9PphTidXoiFO86hpYNldUOMtkp09XGCuRknLiYiIiIyVgxXRAbm5WyFiT39MLGnH/KKb2LH6eq5tPaey8aV/BtYeeACVh64AHtLGfoHuSFSpUSvNq6wlnP3JSIiIjIm/O2MyIg4WVvgiTBPPBHmiRs3K7HvXDZiEzOx/XQmrpWUY/2xK1h/7AoszKV4IMAFUSol+gcr4WorN3TpRERERM0ewxWRkbK0MENUW3dEtXVHRWUV4i9e+2c+rUyk5ZVgZ1IWdiZlQSL5G529HP+5Tssdvi7Whi6diIiIqFliuCJqAszNpOjm54xufs54Y0gwzmYWIeZU9XVaf18pQPzFa4i/eA1zf09CgJuNOmh1aGkPKScuJiIiImoUDFdETYxEIkGguy0C3W3xYv/WuJp/A9v/uU4rLiUXyVlFSM4qwuLdKXCzlatbvIf7O0Nubmbo8omIiIhMFsMVURPXwsESo8N9MDrcBwU3yrH7TBZiEjOxOykLWdfL8MPBNPxwMA02cnP0CXRFpEqJvkFusFPIDF06ERERkUlhuCIyIfaWMgzr1BLDOrVEWUUl4lJyEZNYfVQr+3oZNiekY3NCOmRmEnT3c0aUSokBKiU87C0NXToRERFRk8dwRWSi5OZm6BPohj6Bbnh3WDucuJyPmMRMxJzKQEp2Mfady8G+czmYufEUOrSyR5RKiUiVO9oobThxMREREVE9MFwRNQNSqQQhXo4I8XLE64OCkJJdVN158FQGjl3KR8LlAiRcLsAnMWfh7WylDlqh3o4wY0MMIiIiIq0wXBE1Q/6uNvDvbYPJvf2Rdb0UO05nITYxE/uTc3AxtwRf7kvFl/tS4Wxtgf7BbohUuaNnaxcoZGyIQURERHQnDFdEzZybrQIju3phZFcvFJVVYO/Z6omLd5zORG7xTfx05DJ+OnIZljIz9Gztgqi27ugf5AZHawtDl05ERERkVBiuiEjNRm6OB9t74MH2HiivrMLh1Dz1dVpXC0qr/52YCakE6OLjVD3JsUoJTycrQ5dOREREZHAMV0RUJ5mZFBEBLogIcMGsoSqculqo7jx4Or0QB1PzcDA1D+9sTkSQu606aLVtYceGGERERNQsMVwR0T1JJBK0a2mPdi3tMT2yDS7llfwTtDJwKDUPSRnXkZRxHYt2nEMLewUiVUpEtXVHV18nyMykhi6fiIiIqFEwXBGRzjydrDDhAV9MeMAXecU3sTMpC7GJGdhzNhtXC0rxbdxFfBt3EXYKc/QLckNUW3f0auMKGzl/5BAREZHp4m86RHRfnKwt8HhoKzwe2gql5ZXYfy4HMYkZ2H46C3nFN/Hr8av49fhVWJhL0cPfubohRrAb3GwVhi6diIiISK8YrohIbxQyMwxQKTFApURllcDRtGuIOZWBmMRMXMwtwa4z2dh1JhsSCRDi6YCotu6IVCnh72pj6NKJiIiI7hvDFRE1CDOpBF18nNDFxwn/ezAY57KKEHMqA7GJmThxuQBH0/JxNC0fH/yeBH9Xa0Sq3BHVVolOrRwg5cTFRERE1AQxXBFRg5NIJGijtEUbpS1e6NcaGQWliD1d3eI9LiUXKdnFSNmTgqV7UuBqK8eAYCWi2ioR4e8MuTknLiYiIqKmgeGKiBqdu70Co7p7Y1R3bxSWlmP3mWzEnMrA7jPZyL5ehtWH0rD6UBqsLczQJ9ANUW2V6BPoBntLmaFLJyIiIrojhisiMig7hQwPd2yBhzu2QFlFJf46n6c+fTDrehm2/J2OLX+nw1wqQXc/Z0S1VWJAsBItHCwNXToRERGRBoYrIjIacnMz9G7jit5tXPHOsHZIuFKgDlrnsoqwPzkH+5Nz8NbGU2jf0v6f+bSUCFTacuJiIiIiMjiGKyIySlKpBJ08HdDJ0wGvDQrC+ewixCZmIjYxE/Fp1/D3lQL8faUA82LPwtPJElEqd0SplAj1doQ5Jy4mIiIiA2C4IqImwc/VBpN622BSb39kXy/DjtPVQWtfcg4u5d3A1/tT8fX+VDhaydA/WIkolRI9W7vC0oINMYiIiKhxMFwRUZPjaivHk1298GRXLxSXVWDfuWzEnMrEjqQsXCspx8/xl/Fz/GUoZFL0bO2KKJUS/YOVcLK2MHTpREREZMIYroioSbOWm2NQOw8MaueB8soqHL6Qh5hT1Ue1ruTfUJ9KKJUAYT5OiFIpEaVyh5ezlaFLJyIiIhPDcEVEJkNmJkWEvwsi/F0wa6gKiemF6qCVmF6IQ6l5OJSah3e3nEaQu211QwyVO9q1tGNDDCIiIrpvDFdEZJIkEgnatrBH2xb2mBbZBpfySrD9dCZiTmXi0IU8JGVcR1LGdXy2Mxke9gp10Orm5wQZG2IQERFRPTBcEVGz4OlkhXE9fDGuhy/yS25iZ1IWYk5lYs/ZbKQXlGJV3EWsirsIW4U5+gW5IVKlRO82rrBVcOJiIiIi0g7DFRE1Ow5WFniscys81rkVSssr8WdyDmJOZWL76UzkFt/ExuNXsfH4VViYSRER4IxIlRKRwUq42Snq3F5llcDB1DzE50jgnJqH8AA3mEl5miEREVFzw3BFRM2aQmaG/sHV3QQrqwSOpV1DTGImYk5l4EJuCXafycbuM9l4Y8NJhHg5qE8fDHCzAQD8cTIdszclIr2gFIAZVp07Ag97BWYNVWFQOw/DvjkiIiJqVAxXRET/MJNKEObjhDAfJ0QPDkJyVlF10ErMxIlL+TiWVn376I8z8HOxhr+bDWITM2ttJ6OgFM99fxRLnunMgEVERNSMMFwREdVBIpGgtdIWrZW2eL5vADILSxH7T9CKS8nB+ZxinM8prvO5AoAEwOxNiYhUufMUQSIiomaCLbGIiLSgtFPgme7eWDW+K+JnRuKlfgF3XV8ASC8oxaHUvMYpkIiIiAyO4YqISEd2Chn8/7nm6l7e2ngSK/9MxZX8Gw1cFRERERkaTwskIqoHN9u6Owfe7lxWEd7elIi3NyWibQs7dUOMYA9bTlxMRERkYhiuiIjqoauvEzzsFcgoKIWo43EJAFdbOSY84IvtpzNx5OI1nLpaiFNXC7Fg+zm0crRUB60uPo4w58TFRERETR7DFRFRPZhJJZg1VIXnvj8KCaARsGqOR80Z1haD2nlgUm9/5BSVYefpLMQkZmLfuWxcvnYDK/68gBV/XoCDlQz9gtwQpVKiVxtXWFnwRzMREVFTxP+DExHV06B2HljyTOdb5rmq5l7HPFcuNnIM7+KJ4V08UXKzAvvOVU9cvCMpE/kl5Vh/9ArWH70CubkUPVu7IFJVPfeWi43cEG+NiIiI6oHhiojoPgxq54FIlTvikrMQs+8gonp2Q3iA213br1tZmGNgW3cMbOuOisoqHLl4DTGnMhF7OgOX8m5g++ksbD+dBYnkb4R6OSKqrRKRKnf4ulg34jsjIiIiXTFcERHdJzOpBN18nZB7WqCbr5NO81qZm0nR3c8Z3f2cMfOhYCRlXP9nPq0MnLxSiCMXr+HIxWt4f2sSWrvZVF+n1dYdHVraQ8r5s4iIiIwKwxURkZGQSCQI9rBDsIcdXurfGlfyb2B7YiZiEzPx1/lcnMsqwrmsIizenQKlnRwDgquDVrifMyzM2RCDiIjI0BiuiIiMVEsHS4yJ8MGYCB8UlJRj15ksxCZmYveZLGQWluGHg2n44WAabOTm6BPoikiVEn2D3GCnkBm6dCIiombJ4H/qXLx4MXx9faFQKBAaGop9+/bdcd39+/ejR48ecHZ2hqWlJYKCgjB//vxa6/3yyy9QqVSQy+VQqVTYsGFDQ74FIqIGZ28lwyMhLfHF050RPzMSK8Z1wciuXnC1laOorAKbE9Lx8prjCH0nFqO+Pojv4i4gvYATFxMRETUmgx65Wrt2LaZOnYrFixejR48eWLZsGQYPHozExER4eXnVWt/a2hovvPACOnToAGtra+zfvx+TJk2CtbU1nn32WQBAXFwcRowYgXfeeQePPvooNmzYgOHDh2P//v3o1q1bY79FIiK9U8jM0DfQDX0D3fBeVTscv5xf3RAjMQMp2cXYdy4H+87lYObGU+jQyh5RquqGGG2UNpy4mIiIqAEZNFzNmzcPEyZMwMSJEwEACxYswLZt27BkyRLMnTu31vohISEICQlR3/fx8cH69euxb98+dbhasGABIiMjER0dDQCIjo7Gnj17sGDBAqxevboR3hURUeORSiXo7OWIzl6OmDE4CCnZRdUNMU5l4NilfCRcLkDC5QJ8EnMW3s5WiPznOq1Qb0edGm8QERHRvRksXN28eRPx8fGYMWOGxvKoqCgcOHBAq20cO3YMBw4cwLvvvqteFhcXh2nTpmmsN3DgQCxYsOCO2ykrK0NZWZn6fmFhIQCgvLwc5eXlWtXSkGpqMIZaSD84pqbHWMbUy0GOCRFemBDhhezrZdh5JhvbT2fhwPk8XMwtwVf7U/HV/lQ4WsnQL8gVkUFuiPB3hqWFmUHrNlbGMq6kPxxT08RxNT3GNKa61GCwcJWTk4PKykoolUqN5UqlEhkZGXd9bqtWrZCdnY2Kigq8/fbb6iNfAJCRkaHzNufOnYvZs2fXWh4TEwMrKytt3k6jiI2NNXQJpGccU9NjbGNqC+BRZ2CwA5CUL8HfeRKcuibBtZJy/HL0Kn45ehUyqUCQvUB7J4G2jgI27IdRi7GNK90/jqlp4riaHmMY05KSEq3XNXi3wNvP/xdC3POagH379qGoqAh//fUXZsyYgYCAAIwcObLe24yOjsb06dPV9wsLC+Hp6YmoqCjY2dnp8nYaRHl5OWJjYxEZGQmZjL/1mAKOqelpSmNa/s/ExbGns7HjdBauFpTi72sS/H0NkEqAMG9HDAh2Q/8gV3g5Gc8fmAyhKY0raYdjapo4rqbHmMa05qw2bRgsXLm4uMDMzKzWEaWsrKxaR55u5+vrCwBo3749MjMz8fbbb6vDlbu7u87blMvlkMvltZbLZDKDD+atjK0eun8cU9PTFMZUJgN6BbqjV6A75gwTOHW18J+JizNxOr0Qhy5cw6EL1/D+72cQ5G6rbojRrqVds22I0RTGlXTDMTVNHFfTYwxjqsvrGyxcWVhYIDQ0FLGxsXj00UfVy2NjYzFs2DCttyOE0LheKjw8HLGxsRrXXcXExCAiIkI/hRMRmRCJRIJ2Le3RrqU9pkW2waW8kn+CVgYOX7iGpIzrSMq4jkU7k9HCXoEBKiWiVO7o5ucEmZnBZ/MgIiIyKgY9LXD69OkYNWoUwsLCEB4ejuXLlyMtLQ2TJ08GUH263pUrV7Bq1SoAwBdffAEvLy8EBQUBqJ736pNPPsGLL76o3ubLL7+MXr164cMPP8SwYcOwceNGbN++Hfv372/8N0hE1MR4Ollh/AO+GP+AL64V38TOpOqJi/eczcbVglKsiruIVXEXYaswR78gN0Sp3NE70BU2coOfZU5ERGRwBv2/4YgRI5Cbm4s5c+YgPT0d7dq1w9atW+Ht7Q0ASE9PR1pamnr9qqoqREdHIzU1Febm5vD398cHH3yASZMmqdeJiIjAmjVr8Oabb2LmzJnw9/fH2rVrOccVEZGOHK0t8H+hrfB/oa1QWl6JP5NzEHMqE9tPZyK3+CY2Hr+KjcevwsJMiogAZ0SqlIgMVsLNTmHo0omIiAzC4H9qnDJlCqZMmVLnYytXrtS4/+KLL2ocpbqTxx9/HI8//rg+yiMiIlRPXNw/WIn+wUpUVgkcS7uGmH/m07qQW4LdZ7Kx+0w23thwEp08HRDVtvr0wQA3G0OXTkRE1GgMHq6IiKhpMZNKEObjhDAfJ0QPDkJyVlF10ErMxIlL+Tj+z+2jP87Az8UakW2ViFIpEeLpCCknLiYiIhPGcEVERPUmkUjQWmmL1kpbPN83AJmFperOg3EpOTifU4xle85j2Z7zcLGRY0CwG6LaKhHh7wKFjBMXExGRaWG4IiIivVHaKfBMd288090b10vLsftMNmITM7ErKQs5RWVYc/gS1hy+BCsLM/Ru44qotkr0C1TC3oqtk4mIqOljuCIiogZhq5BhaMcWGNqxBW5WVOFgai5iTmUiNjETGYWl+P1kBn4/mQEzqQTdfJ2qG2KolGjl2LwnLiYioqaL4YqIiBqchbkUPVu7omdrV8wZ1hZ/XymoPn3wVCbOZF7HgZRcHEjJxexNiWjbwg6R/8ynFexh22wnLiYioqaH4YqIiBqVRCJBh1YO6NDKAa9EBeJibrE6aB25mIdTVwtx6mohFmw/h5YOlohqW31Eq6uPE8w5cTERERkxhisiIjIob2drTOzph4k9/ZBbVIYdSVmIOZWJfeeycSX/Blb8eQEr/rwABysZ+gVWN8To1cYVVhb8XxgRERkX/p+JiIiMhrONHMPDPDE8zBMlNyuw71wOYhMzseN0Jq6VlGP9sStYf+xK9WmGAS6Ials995aLjdzQpRMRETFcERGRcbKyMMfAtu4Y2NYdFZVViL9YPXFxbGIm0vJKsCMpCzuSsiCR/I1QL8fq67TausPXxdrQpRMRUTPFcEVEREbP3EyKbn7O6ObnjDeHBONM5nXEnqqeT+vvKwU4cvEajly8hrm/JyHAzQZR/wStDi3tOXExERE1GoYrIiJqUiQSCYLc7RDkbocX+7fG1fwb2H66uiHGX+dzkZxVhOSsIizenQI3W7m6xXu4vzPk5py4mIiIGg7DFRERNWktHCwxOtwHo8N9UHCjHLvPZCEmMRO7k7KQdb0MPxxMww8H02AjN0efQFdEqpToG+QGOwUnLiYiIv1iuCIiIpNhbynDsE4tMaxTS5RVVCIuJVd9nVb29TJsTkjH5oR0yMwk6O7njCiVEgNUSnjYWxq6dCIiMgEMV0REZJLk5mboE+iGPoFueHdYO5y4nK8OWslZRdh3Lgf7zuVg5sZT6NDKHpHB1ddptVHacOJiIiKqF4YrIiIyeVKpBCFejgjxcsTrg4JwPruoeuLixEwcTbuGhMsFSLhcgE9jz8LLyerfhhgtbAxdOhERNSEMV0RE1Oz4udpgUm8bTOrtj+zrZdhxujpo7U/OQVpeCb7an4qv9qfC0UqGNtZSWJzOQp8gd1hasCEGERHdGcMVERE1a662cjzZ1QtPdvVCcVkF9p7Nrp64OCkL10rKcbBEioM/HodCJkXP1q6IUlVPXOxkbWHo0omIyMgwXBEREf3DWm6Owe09MLi9B8orq/BXcja+/P0QUkqtcCW/FLH/XLMllQBhPk7Vpw+q3OHlbGXo0omIyAgwXBEREdVBZiZFdz8n5PlWYfDgnjiXcwMxp6rDVWJ6IQ6l5uFQah7e3XIaQe62iPwnaLVraceGGEREzRTDFRER0T1IJBK0bWGPti3sMS2yDS7llagnLj50IQ9JGdeRlHEdn+1Mhoe9Qh20uvk5QWYmNXT5RETUSBiuiIiIdOTpZIVxPXwxrocv8ktuYmdSFmITM7HnbDbSC0qxKu4iVsVdhK3CHP2C3BCpUqJ3G1fYcuJiIiKTxnBFRER0HxysLPBY51Z4rHMrlJZX4s/kHMQmZmL76UzkFN3ExuNXsfH4VViYSRER4IxIlRKRwUq42SkMXToREekZwxUREZGeKGRm6B9c3U2wskrg+KVriDlV3eY9NacYu89kY/eZbLyx4SQ6eTogqq0SUSol/F05cTERkSlguCIiImoAZlIJQr2dEOrthBmDg5CSXYRt/zTEOH4pX3376I8z8HOxrr5Oq60SnTwdYSZl0CIiaooYroiIiBqYRCJBgJstAtxs8XzfAGQWlqobYsSl5OJ8TjGW7T2PZXvPw8XGAgOCq4NWhL8LFDJOXExE1FQwXBERETUypZ0CT3fzxtPdvHG9tBx7/pm4eGdSFnKKbmLN4UtYc/gSrCzM0LuNKyJVSvQLcoODFScuJiIyZgxXREREBmSrkOGhDi3wUIcWuFlRhUOpeYhJzEDMqUxkFJbi95MZ+P1kBsykEnT1cUJUWyUiVUq0cuTExURExobhioiIyEhYmEvxQGsXPNDaBbMfbouTVwoRk5iB2MRMJGVcR9z5XMSdz8XsTYlQedipg5bKgxMXExEZA4YrIiIiIySRSNC+lT3at7LHK1GBSMstqT6ilZiJIxfykJheiMT0QizYfg4tHSzVDTG6+jjBnBMXExEZBMMVERFRE+DlbIWJPf0wsacf8opvYsfp6hbv+85l40r+Daw8cAErD1yAvaUM/YPcENVWiZ6tXWEt5//qiYgaC3/iEhERNTFO1hZ4IswTT4R54sbNSuw7l62euPhaSTnWH7uC9ceuwMJcip4BLohUVc+95WorN3TpREQmjeGKiIioCbO0MENUW3dEtXVHRWUV4i9eQ2xi9VGttLwS7EjKwo6kLEgkf6OzlyOiVEpEtXWHr4u1oUsnIjI59Q5X8fHxOH36NCQSCYKDg9G5c2d91kVEREQ6MjeTopufM7r5OeONIcE4m1mEmFMZiD2diYTLBYi/eA3xF69h7u9JCHCzQZSquiFGx1YOkHLiYiKi+6ZzuMrKysKTTz6J3bt3w8HBAUIIFBQUoG/fvlizZg1cXV0bok4iIiLSgUQiQaC7LQLdbfFi/9ZIL7iB7f8c0YpLyUVyVhGSs4qweHcK3GzlGKBSIkqlRLi/M+TmnLiYiKg+dA5XL774IgoLC3Hq1CkEBwcDABITEzFmzBi89NJLWL16td6LJCIiovvjYW+JUeE+GBXug4Ib5dh9JguxiZnYfSYbWdfL8OPBNPx4MA02cnP0DnRFlEqJPoFusLeUGbp0IqImQ+dw9ccff2D79u3qYAUAKpUKX3zxBaKiovRaHBEREemfvaUMwzq1xLBOLVFWUYm4lFzEJmYiNjETWdfLsCUhHVsS0mEulSDc3xmR/5w+6GFvaejSiYiMms7hqqqqCjJZ7b9iyWQyVFVV6aUoIiIiahxyczP0CXRDn0A3vDOsHRKuFFRfp5WYiXNZRdh3Lgf7zuXgrY2n0L6lffV1Wm2VCFTacuJiIqLb6Byu+vXrh5dffhmrV69GixYtAABXrlzBtGnT0L9/f70XSERERI1DKpWgk6cDOnk64LVBQTifXaQ+ohWfdg1/XynA31cK8GnsWXg5WakbYoT5OMGMDTGIiHQPV59//jmGDRsGHx8feHp6QiKRIC0tDe3bt8f333/fEDUSERGRAfi52mBSbxtM6u2P7Otl2JmUiZhTmdiXnIO0vBJ8tT8VX+1PhZO1BfoFuSFKVT1xsaUFG2IQUfOkc7jy9PTE0aNHERsbi6SkJAghoFKpMGDAgIaoj4iIiIyAq60cI7p4YUQXLxSXVWDfuWzEJGZix+ks5BXfxM/xl/Fz/GUoZFL0bF3dEKN/sBJO1haGLp2IqNHUe56ryMhIREZG6rMWIiIiagKs5eYY1M4Dg9p5oKKyCocu5FVPXHwqE1fyb6hPJZRKgDAfJ/Xpg97OnLiYiEybVuFq0aJFePbZZ6FQKLBo0aK7rvvSSy/ppTAiIiIyfuZmUkT4uyDC3wVvPaTC6fTriEmsbohx6mohDqXm4VBqHt7dchqBSltEtVUiSuWOdi3t2BCDiEyOVuFq/vz5ePrpp6FQKDB//vw7rieRSBiuiIiImimJRAJVCzuoWthh6oA2uHytRD1x8cHUPJzJvI4zmdfx2c5keNgr1C3eu/k6w8JcaujyiYjum1bhKjU1tc5/ExEREd1JK0crjO3hi7E9fJFfchO7zmQh5lQm9pzNRnpBKVbFXcSquIuwVZijb6Abotoq0buNK2wVnLiYiJomna+5mjNnDl599VVYWVlpLL9x4wY+/vhjvPXWW3orjoiIiEyDg5UFHg1phUdDWqG0vBIHUnLU12blFN3Ebyeu4rcTV2FhJkW4vzOi2ioRGayEm53C0KUTEWlN52Pws2fPRlFRUa3lJSUlmD17tl6KIiIiItOlkJmhX5AScx/rgIP/G4BfngvHpN5+8HOxxs3KKuw5m403NpxE1/d34JEv/sQXu5KRnHUdQghDl05EdFc6H7kSQtR5AeqJEyfg5OSkl6KIiIioeTCTShDq7YRQbydEDw5GclaRuiHGsbR8HL9Ufft42xn4ulirOw+GeDly4mIiMjpahytHR0dIJBJIJBK0adNGI2BVVlaiqKgIkydPbpAiiYiIqHkIcLNBgFsApvQJQFZhKbafzkJMYgYOJOciNacYy/aex7K95+FiY4EBwdVBq0eACxQyTlxMRIandbhasGABhBAYP348Zs+eDXt7e/VjFhYW8PHxQXh4eIMUSURERM2Pm50CT3XzwlPdvFBUVoE9Z7IRk5iBnUlZyCm6iTWHL2HN4UuwsjBDr9auiGqrRL8gNzhY1T1xcWWVwMHUPMTnSOCcmofwADce/SIivdI6XI0ZMwYA4Ovri4iICMhk7ORDREREjcNGbo4hHTwwpIMHyiurcPB8HmITMxCTmIn0glL8cSoDf5zKgJlUgq4+TtUNMVRKtHKsbsD1x8l0zN6UiPSCUgBmWHXuCDzsFZg1VIVB7TwM++aIyGTofM1V79691f++ceMGysvLNR63s7O7/6qIiIiI7kBmJsUDrV3wQGsXvP1wW5y6WoiYU9VBKynjOuLO5yLufC5mb0qEysMOvi7W2PJ3eq3tZBSU4rnvj2LJM50ZsIhIL3QOVyUlJXjttdfw008/ITc3t9bjlZWVeimMiIiI6F4kEgnatbRHu5b2mB4ViLTcEnVDjMMX8pCYXojE9MI6nysASADM3pSISJU7TxEkovumcyv2//73v9i5cycWL14MuVyOr776CrNnz0aLFi2watWqhqiRiIiISCtezlaY2NMPayeF48ibkZjc2/+u6wsA6QWlOJSa1zgFEpFJ0/nI1aZNm7Bq1Sr06dMH48ePR8+ePREQEABvb2/88MMPePrppxuiTiIiIiKdOFlbINjDVqt1s66XNnA1RNQc6HzkKi8vD76+vgCqr6/Ky6v+S88DDzyAvXv36rc6IiIiovvgZqvQ63pERHejc7jy8/PDhQsXAAAqlQo//fQTgOojWg4ODvqsjYiIiOi+dPV1goe9Ane7mkppJ0dXX6dGq4mITJfO4WrcuHE4ceIEACA6Olp97dW0adPw3//+V+8FEhEREdWXmVSCWUNVAHDHgGUjN0d5ZVXjFUVEJkvna66mTZum/nffvn2RlJSEI0eOwN/fHx07dtRrcURERET3a1A7Dyx5pvMt81xVc7WV43ppOVKyi/HazwlY+GQnSCTsGEhE9adzuLqdl5cXvLy8AAA///wzHn/88fsuioiIiEifBrXzQKTKHXHJWYjZdxBRPbshPMANf53PxZhvDuG3E1fh42yF6VGBhi6ViJownU4LrKiowKlTp3D27FmN5Rs3bkTHjh3ZKZCIiIiMlplUgm6+Tgh1Eejm6wQzqQQ9Alzw/qPtAQCLdibjl/jLBq6SiJoyrcNVYmIi2rRpgw4dOiA4OBiPPfYYMjMz0bt3b4wZMwaRkZFITk5uyFqJiIiI9G54F09M6VM9H9aM9Qn463yugSsioqZK63A1Y8YM+Pr6YuPGjRg+fDh+/fVX9OzZE/3798elS5fwySefwNPTsyFrJSIiImoQr0YFYkh7D5RXCkz6Lh4p2UWGLomImiCtw9WhQ4fw8ccf46GHHsKSJUsAAP/973/x1ltvwdZWuwn6iIiIiIyRVCrBp8M7IsTLAQU3yjF+5WHkFd80dFlE1MRoHa6ysrLQsmVLAICDgwOsrKzQu3fv+y5g8eLF8PX1hUKhQGhoKPbt23fHddevX4/IyEi4urrCzs4O4eHh2LZtm8Y6K1euhEQiqXUrLeXM60RERHRnCpkZvhwdBk8nS1zMLcGzq46gtLzS0GURUROidbiSSCSQSv9dXSqVQiaT3deLr127FlOnTsUbb7yBY8eOoWfPnhg8eDDS0tLqXH/v3r2IjIzE1q1bER8fj759+2Lo0KE4duyYxnp2dnZIT0/XuCkUnHmdiIiI7s7FRo4VY7vAVmGOIxev4bWfEyCEMHRZRNREaN2KXQiBNm3aqOd/KCoqQkhIiEbgAoC8vDytX3zevHmYMGECJk6cCABYsGABtm3bhiVLlmDu3Lm11l+wYIHG/ffffx8bN27Epk2bEBISol4ukUjg7u6udR1lZWUoKytT3y8sLAQAlJeXo7y8XOvtNJSaGoyhFtIPjqnp4ZiaJo6r6dFmTL0dFfj8yY6YsOoofjtxFZ6OCkztH9BYJVI9cF81PcY0prrUoHW4WrFiRb2KuZObN28iPj4eM2bM0FgeFRWFAwcOaLWNqqoqXL9+HU5OThrLi4qK4O3tjcrKSnTq1AnvvPOORvi63dy5czF79uxay2NiYmBlZaVVLY0hNjbW0CWQnnFMTQ/H1DRxXE2PNmP6hK8Eq1PM8MXu88i/fA5dXXkEy9hxXzU9xjCmJSUlWq+rdbgaM2ZMvYq5k5ycHFRWVkKpVGosVyqVyMjI0Gobn376KYqLizF8+HD1sqCgIKxcuRLt27dHYWEhFi5ciB49euDEiRNo3bp1nduJjo7G9OnT1fcLCwvh6emJqKgo2NnZ1ePd6Vd5eTliY2MRGRl536diknHgmJoejqlp4riaHl3G9EEA9rHnsHRvKn5KNcfgXqHo5ut01+eQYXBfNT3GNKY1Z7VpQ+tw1VBqTjOsIYSotawuq1evxttvv42NGzfCzc1Nvbx79+7o3r27+n6PHj3QuXNnfPbZZ1i0aFGd25LL5ZDL5bWWy2Qygw/mrYytHrp/HFPTwzE1TRxX06PtmL42KBiXrpViy9/peH71CayfEgF/V5tGqJDqg/uq6TGGMdXl9bVuaKFvLi4uMDMzq3WUKisrq9bRrNutXbsWEyZMwE8//YQBAwbcdV2pVIouXbrg3Llz910zERERNS9s0U5EujBYuLKwsEBoaGit8yhjY2MRERFxx+etXr0aY8eOxY8//oghQ4bc83WEEDh+/Dg8PDzuu2YiIiJqftiinYi0ZbBwBQDTp0/HV199hW+++QanT5/GtGnTkJaWhsmTJwOovhZq9OjR6vVXr16N0aNH49NPP0X37t2RkZGBjIwMFBQUqNeZPXs2tm3bhvPnz+P48eOYMGECjh8/rt4mERERka7Yop2ItFHvcHXz5k2cOXMGFRUV9X7xESNGYMGCBZgzZw46deqEvXv3YuvWrfD29gYApKena8x5tWzZMlRUVOD555+Hh4eH+vbyyy+r18nPz8ezzz6L4OBgREVF4cqVK9i7dy+6du1a7zqJiIiIAtxsseyZUJhLJfjtxFXMjz1r6JKIyMjo3NCipKQEL774Ir799lsAwNmzZ+Hn54eXXnoJLVq0qNVa/V6mTJmCKVOm1PnYypUrNe7v3r37ntubP38+5s+fr1MNRERERNqICHDB+4+2x2u/JGDRzmR4O1vj/0JbGbosIjISOh+5io6OxokTJ7B7924oFAr18gEDBmDt2rV6LY6IiIjI2Azv4okpffwBADPWJyAuJdfAFRGRsdA5XP3666/4/PPP8cADD2i0TFepVEhJSdFrcURERETG6NWoQAxp74HySoHJ38cjJbvI0CURkRHQOVxlZ2drzCtVo7i4WKv5qYiIiIiaOrZoJ6K66ByuunTpgi1btqjv1wSqL7/8EuHh4fqrjIiIiMiIsUU7Ed1O54YWc+fOxaBBg5CYmIiKigosXLgQp06dQlxcHPbs2dMQNRIREREZpZoW7Y8uPqBu0b7wyU48m4eomdL5yFVERAT+/PNPlJSUwN/fHzExMVAqlYiLi0NoaGhD1EhERERktNiinYhq6HzkCgDat2+vbsVORERE1NyxRTsRAfU4crV161Zs27at1vJt27bh999/10tRRERERE0NW7QTkc7hasaMGaisrH2xphBC5wmEiYiIiEwJW7QTNW86h6tz585BpVLVWh4UFITk5GS9FEVERETUFLFFO1HzpnO4sre3x/nz52stT05OhrW1tV6KIiIiImqq2KKdqPnSOVw9/PDDmDp1KlJSUtTLkpOT8corr+Dhhx/Wa3FERERETVFNi3Zbhbm6RbsQwtBlEVED0zlcffzxx7C2tkZQUBB8fX3h6+uL4OBgODs745NPPmmIGomIiIiaHLZoJ2p+dG7Fbm9vjwMHDiA2NhYnTpyApaUlOnTogF69ejVEfURERERNFlu0EzUv9ZrnSiKRICoqClFRUfquh4iIiMikDO/iiQu5xVi8OwUz1ieghYMlwv2dDV0WETWAeoWrHTt2YMeOHcjKykJVVZXGY998841eCiMiIiIyFa9GBeJibgm2/J2Oyd/HY/2UCPi72hi6LCLSM52vuZo9ezaioqKwY8cO5OTk4Nq1axo3IiIiItLEFu1EzYPOR66WLl2KlStXYtSoUQ1RDxEREZFJqmnR/ujiP9Ut2r+f2A0KmZmhSyMiPdH5yNXNmzcRERHRELUQERERmTS2aCcybTqHq4kTJ+LHH39siFqIiIiITB5btBOZLp1PCywtLcXy5cuxfft2dOjQATKZTOPxefPm6a04IiIiIlPEFu1EpknncJWQkIBOnToBAE6ePKnxmEQi0UtRRERERKaOLdqJTI/O4WrXrl0NUQcRERFRs8MW7USmRedrroiIiIhIP9iinci01GsS4cOHD2PdunVIS0vDzZuaPwDWr1+vl8KIiIiImgO2aCcyHTofuVqzZg169OiBxMREbNiwAeXl5UhMTMTOnTthb2/fEDUSERERmTS2aCcyDTqHq/fffx/z58/H5s2bYWFhgYULF+L06dMYPnw4vLy8GqJGIiIiIpPHFu1ETZ/O4SolJQVDhgwBAMjlchQXF0MikWDatGlYvny53gskIiIiai5qWrQDwKKdyfgl/rKBKyIiXegcrpycnHD9+nUAQMuWLdXt2PPz81FSUqLf6oiIiIiameFdPDGljz8AYMb6BMSl5Bq4IiLSls7hqmfPnoiNjQUADB8+HC+//DL+85//YOTIkejfv7/eCyQiIiJqbl6NCsSQDh4orxSY/H08UrKLDF0SEWlB526Bn3/+OUpLSwEA0dHRkMlk2L9/Px577DHMnDlT7wUSERERNTdSqQSfPtERV/Nv4FhaPsavPIwNU3rAydrC0KUR0V3U67TAFi1aVD9ZKsVrr72G3377DfPmzYOjo6PeCyQiIiJqjmpatHs6WapbtJeWVxq6LCK6C63CVWFhoca/73YjIiIiIv1gi3aipkWrcOXo6IisrCwAgIODAxwdHWvdapYTERERkf6wRTtR06HVNVc7d+6Ek5MTAGDXrl0NWhARERERaapp0f7aLwlYtDMZ3s7W+L/QVoYui4huo1W46t27NwCgoqICu3fvxvjx4+Hp6dmghRERERHRv4Z38cSF3GIs3p2CGesT0MLBEuH+zoYui4huoVNDC3Nzc3zyySeorOTFlERERESNjS3aiYybzt0C+/fvj927dzdAKURERER0NzUt2kO8HFBwoxzjVx5GXvFNQ5dFRP/QeZ6rwYMHIzo6GidPnkRoaCisra01Hn/44Yf1VhwRERERaapp0f7o4j/VLdq/n9gNCpmZoUsjavZ0DlfPPfccAGDevHm1HpNIJDxlkIiIiKiB1bRof3TxAXWL9oVPdoJEIjF0aUTNms6nBVZVVd3xxmBFRERE1DjYop3I+OgcroiIiIjIONS0aAeARTuT8XP8ZQNXRNS86XxaIAAUFxdjz549SEtLw82bmhdRvvTSS3opjIiIiIju7dYW7dHrE9CSLdqJDEbncHXs2DE8+OCDKCkpQXFxMZycnJCTkwMrKyu4ubkxXBERERE1slejAnExrwRbEtIx+ft4rJ8SAX9XG0OXRdTs6Hxa4LRp0zB06FDk5eXB0tISf/31Fy5evIjQ0FB88sknDVEjEREREd0FW7QTGQedw9Xx48fxyiuvwMzMDGZmZigrK4Onpyc++ugj/O9//2uIGomIiIjoHmpatHs6WapbtJeWs9kYUWPSOVzJZDJ1m0+lUom0tDQAgL29vfrfRERERNT4alq02yrM1S3ahRCGLouo2dA5XIWEhODIkSMAgL59++Ktt97CDz/8gKlTp6J9+/Z6L5CIiIiItMcW7USGo3W4qqioAAC8//778PDwAAC88847cHZ2xnPPPYesrCwsX768YaokIiIiIq2xRTuRYWjdLdDDwwNjxozB+PHjERYWBgBwdXXF1q1bG6w4IiIiIqoftmgnanxaH7maPn06Nm3ahPbt2yM8PBxff/01ioqKGrI2IiIiIroPr0YFYkgHD5RXCkz+Ph4p2fzdjaghaR2uoqOjcebMGezevRtBQUGYOnUqPDw8MG7cOPz5558NWSMRERER1QNbtBM1Lp0bWvTs2RMrVqxARkYGFixYgOTkZPTs2ROBgYH46KOPGqJGIiIiIqontmgnajw6h6sa1tbWmDBhAvbt24dNmzYhJycH0dHR+qyNiIiIiPSALdqJGke9w1VJSQlWrFiBXr164eGHH4azszPee+89fdZGRERERHrCFu1EDU/ncLVv3z6MHz8e7u7ueOGFF+Dr64tdu3bh7NmzmDFjRkPUSERERER6wBbtRA1L61bs77//PlauXImUlBSEhYXh448/xsiRI2FnZ9eQ9RERERGRHrFFO1HD0frI1fz58zFkyBCcOHECBw8exKRJkxisiIiIiJogtmgnahhaH7m6evUqZDJZQ9ZCRERERI2gpkX71fwbOJaWj/ErD2PDlB5wsrYwdGlETZrWR64YrIiIiIhMB1u0E+lfvbsFEhEREVHTxhbtRPrFcEVERETUjLFFO5H+GDxcLV68GL6+vlAoFAgNDcW+ffvuuO769esRGRkJV1dX2NnZITw8HNu2bau13i+//AKVSgW5XA6VSoUNGzY05FsgIiIiatIiAlzw/mNs0U50v7QKV4WFhVrfdLF27VpMnToVb7zxBo4dO4aePXti8ODBSEtLq3P9vXv3IjIyElu3bkV8fDz69u2LoUOH4tixY+p14uLiMGLECIwaNQonTpzAqFGjMHz4cBw8eFCn2oiIiIiak+FhnpjSxx8AEL0+AXEpuQauiKjp0apboIODAyQSiVYbrKzU/kLIefPmYcKECZg4cSIAYMGCBdi2bRuWLFmCuXPn1lp/wYIFGvfff/99bNy4EZs2bUJISIh6ncjISERHRwMAoqOjsWfPHixYsACrV6/WujYiIiKi5ubVqEBczCvBloR0TP4+HuunRMDf1cbQZRE1GVqFq127dqn/feHCBcyYMQNjx45FeHg4gOqjRd9++22dgehObt68ifj4eMyYMUNjeVRUFA4cOKDVNqqqqnD9+nU4OTmpl8XFxWHatGka6w0cOLBWMLtVWVkZysrK1PdrjsCVl5ejvLxcq1oaUk0NxlAL6QfH1PRwTE0Tx9X0cEzv7YNHVLhyrQTHLxVg7DeHsG5SNzgbeYt2jqvpMaYx1aUGrcJV79691f+eM2cO5s2bh5EjR6qXPfzww2jfvj2WL1+OMWPGaPXCOTk5qKyshFKp1FiuVCqRkZGh1TY+/fRTFBcXY/jw4eplGRkZOm9z7ty5mD17dq3lMTExsLKy0qqWxhAbG2voEkjPOKamh2NqmjiupodjenePK4FLWWa4dO0GRn6xC8+rKiEz+JX698ZxNT3GMKYlJSVar6v1JMI14uLisHTp0lrLw8LC1Kf36eL20w2FEFqdgrh69Wq8/fbb2LhxI9zc3O5rm9HR0Zg+fbr6fmFhITw9PREVFQU7Oztt3kaDKi8vR2xsLCIjIznfmIngmJoejqlp4riaHo6p9kLDizD8y0NIvV6B3SWtMO+J9lpfJtLYOK6mx5jGVJe+EjqHK09PTyxduhSffvqpxvJly5bB09NT6+24uLjAzMys1hGlrKysWkeebrd27VpMmDAB69atw4ABAzQec3d313mbcrkccrm81nKZTGbwwbyVsdVD949jano4pqaJ42p6OKb3FtzSEcueCcXobw5h898Z8HO1wfSoQEOXdVccV9NjDGOqy+vrfIB3/vz5WLx4Mdq1a4eJEydi4sSJaNeuHRYvXoz58+drvR0LCwuEhobWOtQXGxuLiIiIOz5v9erVGDt2LH788UcMGTKk1uPh4eG1thkTE3PXbRIRERFRbWzRTqQbnY9cPfjggzh79iyWLFmCpKQkCCEwbNgwTJ48WacjVwAwffp0jBo1CmFhYQgPD8fy5cuRlpaGyZMnA6g+Xe/KlStYtWoVgOpgNXr0aCxcuBDdu3dXH6GytLSEvb09AODll19Gr1698OGHH2LYsGHYuHEjtm/fjv379+v6VomIiIiaveFhnriQU4zFu1MQvT4BLR0sEe7vbOiyiIySzuEKqD418P3337/vFx8xYgRyc3MxZ84cpKeno127dti6dSu8vb0BAOnp6RpzXi1btgwVFRV4/vnn8fzzz6uXjxkzBitXrgQAREREYM2aNXjzzTcxc+ZM+Pv7Y+3atejWrdt910tERETUHLFFO5F26hWu9u3bh2XLluH8+fNYt24dWrZsie+++w6+vr544IEHdNrWlClTMGXKlDofqwlMNXbv3q3VNh9//HE8/vjjOtVBRERERHWTSiX49ImOuJp/A8fS8jFuxWFsmBIBZ5va16wTNWc6X3P1yy+/YODAgbC0tMTRo0fV80Ndv35dL0eziIiIiMj4KGRm+HJ0GDydLJGWV4Jnv4tHaXmlocsiMio6h6t3330XS5cuxZdffqnROSMiIgJHjx7Va3FEREREZDxcbORYMbYLbBXmiL94Da/9nAAhhKHLIjIaOoerM2fOoFevXrWW29nZIT8/Xx81EREREZGRCnCzxbJnQmEuleC3E1cxP/asoUsiMho6hysPDw8kJyfXWr5//374+fnppSgiIiIiMl5s0U5UN53D1aRJk/Dyyy/j4MGDkEgkuHr1Kn744Qe8+uqrd2xMQURERESmZXiYJ57v6w8AiF6fgLiUXANXRGR4OncLfO2111BQUIC+ffuitLQUvXr1glwux6uvvooXXnihIWokIiIiIiP0SmQgLuaWYDNbtBMBqMeRKwB47733kJOTg0OHDuGvv/5CdnY23nnnHX3XRkRERERGTCqV4JMnOiLEywEFN8oxbsVh5BaVGbosIoOpV7gCACsrK4SFhaFr166wseFfKIiIiIiaI7ZoJ/qXzuGquLgYM2fOREREBAICAuDn56dxIyIiIqLm5fYW7f/9OQFVVWzRTs2PztdcTZw4EXv27MGoUaPg4eEBiUTSEHURERERURNS06J99DeHsOnEVfg4W+GVqEBDl0XUqHQOV7///ju2bNmCHj16NEQ9RERERNRE1bRof+3nBHy2MxneztZ4PLSVocsiajQ6nxbo6OgIJyenhqiFiIiIiJo4tmin5kzncPXOO+/grbfeQklJSUPUQ0RERERN3CuRgXiogwfKKwUmfx+PlOwiQ5dE1Ch0Pi3w008/RUpKCpRKJXx8fCCTyTQeP3r0qN6KIyIiIqKmp6ZF+5X8GziWlo9xKw5jw5QIONvIDV0aUYPSOVw98sgjDVAGEREREZmSmhbtjy7+U92i/YeJ3aCQmRm6NKIGo3O4mjVrVkPUQUREREQmpqZF+6OLD6hbtC8c0QlSKbtNk2mq9yTCRERERET3UtOi3VwqwaYTVzF/+1lDl0TUYLQKV05OTsjJyQHwb7fAO92IiIiIiG5V06IdAD7bmYyf4y8buCKihqHVaYHz58+Hra0tAGDBggUNWQ8RERERmaDhYZ64mFuML3alIHp9Alo6WCLc39nQZRHplVbhasyYMXX+m4iIiIhIW69EBuJibgk2J6Rj8vfxWD8lAv6uNoYui0hv7uuaqxs3bqCwsFDjRkRERERUl5oW7Z29HFBwoxzjVhxGblGZocsi0hudw1VxcTFeeOEFuLm5wcbGBo6Ojho3IiIiIqI7qWnR7ulkqW7RXlpeaeiyiPRC53D12muvYefOnVi8eDHkcjm++uorzJ49Gy1atMCqVasaokYiIiIiMiHO/7Rot1WYq1u0V1UJQ5dFdN90DlebNm3C4sWL8fjjj8Pc3Bw9e/bEm2++iffffx8//PBDQ9RIRERERCaGLdrJFOkcrvLy8uDr6wsAsLOzQ15eHgDggQcewN69e/VbHRERERGZLLZoJ1Ojc7jy8/PDhQsXAAAqlQo//fQTgOojWg4ODvqsjYiIiIhM3PAwTzzf1x8AEL0+AXEpuQauiKj+dA5X48aNw4kTJwAA0dHR6muvpk2bhv/+9796L5CIiIiITNsrkYF4qIMHyisFJn8fj5TsIkOXRFQvWs1zdatp06ap/923b18kJSXhyJEj8Pf3R8eOHfVaHBERERGZvpoW7Vfzb+BoWj7GrTiMDVMi4GwjN3RpRDq5r3muAMDLywuPPfYYgxURERER1RtbtJMp0OrI1aJFi7Te4EsvvVTvYoiIiIio+app0f7o4gPqFu0LR3SCVCoxdGlEWtEqXM2fP1+rjUkkEoYrIiIiIqq3mhbto785hE0nrsLH2QqvRAUauiwirWgVrlJTUxu6DiIiIiIiAP+2aH/t5wR8tjMZ3s7WeDy0laHLIrqn+7rmSggBITibNhERERHpF1u0U1NUr3D19ddfo127dlAoFFAoFGjXrh2++uorfddGRERERM0YW7RTU6NzuJo5cyZefvllDB06FOvWrcO6deswdOhQTJs2DW+++WZD1EhEREREzVBNi/bOXg4ouFGOcSsOI7eozNBlEd2RzvNcLVmyBF9++SVGjhypXvbwww+jQ4cOePHFF/Huu+/qtUAiIiIiar5qWrQ/svhPdYv2HyZ2g0JmZujSiGrR+chVZWUlwsLCai0PDQ1FRUWFXooiIiIiIqpR06LdTmGubtFeVcXr/sn46ByunnnmGSxZsqTW8uXLl+Ppp5/WS1FERERERLcKcLPF0mdCYS6VYNOJq5i//ayhSyKqRefTAoHqhhYxMTHo3r07AOCvv/7CpUuXMHr0aEyfPl293rx58/RTJRERERE1e2zRTsZO53B18uRJdO7cGQCQkpICAHB1dYWrqytOnjypXk8i4UzaRERERKRfw8M8cTG3GF/sSkH0+gS0dLBEmJedocsiAlCPcLVr166GqIOIiIiISCuvRAbiYm4JNiekY9J3R/DTs90MXRIRgHpcc5WZmXnHxxISEu6rGCIiIiKie7m1RXthaQX+891RFJUbuiqieoSr9u3b47fffqu1/JNPPkG3bvyrARERERE1vJoW7Z5Olrh07Qa+OmOGsvJKQ5dFzZzO4er111/HiBEjMHnyZNy4cQNXrlxBv3798PHHH2Pt2rUNUSMRERERUS23tmhPvS7B6xtOsUU7GZTO4eqVV17BX3/9hT///BMdOnRAhw4dYGlpiYSEBDz88MMNUSMRERERUZ0C3Gzx+ciOkEoEtvydwRbtZFA6hysA8PPzQ9u2bXHhwgUUFhZi+PDhUCqV+q6NiIiIiOiewv2cMcKvCgDw2c5k/Bx/2cAVUXOlc7iqOWKVnJyMhIQELFmyBC+++CKGDx+Oa9euNUSNRERERER31d1N4LlevgCA6PUJiEvJNXBF1BzpHK769euHESNGIC4uDsHBwZg4cSKOHTuGy5cvo3379g1RIxERERHRPU3tH4CHOnigvFJg0ndHkJxVZOiSqJnROVzFxMTggw8+gEwmUy/z9/fH/v37MWnSJL0WR0RERESkrdtbtI9feRi5RWWGLouaEZ3DVe/evevekFSKmTNn3ndBRERERET1dWuL9rS8Ejz7XTxK2aKdGonW4erBBx9EQUGB+v57772H/Px89f3c3FyoVCq9FkdEREREpKtbW7THX7yG//6cwBbt1Ci0Dlfbtm1DWdm/h1U//PBD5OXlqe9XVFTgzJkz+q2OiIiIiKgeAtxssfSZUJhLJdh04ipbtFOj0DpcCSHuep+IiIiIyJhEBLjg/ceqG66xRTs1hnrNc0VERERE1BQMD/PE8339AbBFOzU8rcOVRCKBRCKptYyIiIiIyJi9EhnIFu3UKMy1XVEIgbFjx0IulwMASktLMXnyZFhbWwOAxvVYRERERETGoqZF+9X8Gzialo/xKw9jw5QIONvIDV0amRitj1yNGTMGbm5usLe3h729PZ555hm0aNFCfd/NzQ2jR49uyFqJiIiIiOqFLdqpMWh95GrFihUNWQcRERERUYOqadH+2OID6hbtC0d0glTKS11IP9jQgoiIiIiaDbZop4bEcEVEREREzQpbtFNDYbgiIiIiomaHLdqpITBcEREREVGzxBbtpG8GD1eLFy+Gr68vFAoFQkNDsW/fvjuum56ejqeeegqBgYGQSqWYOnVqrXVWrlypnpPr1ltpaWkDvgsiIiIiampqWrR39nJAYWkFxq88jNwiTi9E9WfQcLV27VpMnToVb7zxBo4dO4aePXti8ODBSEtLq3P9srIyuLq64o033kDHjh3vuF07Ozukp6dr3BQKRUO9DSIiIiJqotiinfTJoOFq3rx5mDBhAiZOnIjg4GAsWLAAnp6eWLJkSZ3r+/j4YOHChRg9ejTs7e3vuF2JRAJ3d3eNGxERERFRXWpatNspzNUt2quqhKHLoiZI63mu9O3mzZuIj4/HjBkzNJZHRUXhwIED97XtoqIieHt7o7KyEp06dcI777yDkJCQO65fVlaGsrJ/DwEXFhYCAMrLy1FeXn5ftehDTQ3GUAvpB8fU9HBMTRPH1fRwTE2TPsbV21GBz0d2xPhvj2LTiavwdFBg2oAAfZVIOjKmfVWXGgwWrnJyclBZWQmlUqmxXKlUIiMjo97bDQoKwsqVK9G+fXsUFhZi4cKF6NGjB06cOIHWrVvX+Zy5c+di9uzZtZbHxMTAysqq3rXoW2xsrKFLID3jmJoejqlp4riaHo6padLHuD7hK8HqFDMs3nMe+ZfPoZsbj2AZkjHsqyUlJVqva7BwVUMi0ZwRWwhRa5kuunfvju7du6vv9+jRA507d8Znn32GRYsW1fmc6OhoTJ8+XX2/sLAQnp6eiIqKgp2dXb1r0Zfy8nLExsYiMjISMpnM0OWQHnBMTQ/H1DRxXE0Px9Q06XNcHwTgEHsOS/amYt0FczzYOxTdfJ30UyhpzZj21Zqz2rRhsHDl4uICMzOzWkepsrKyah3Nuh9SqRRdunTBuXPn7riOXC6HXC6vtVwmkxl8MG9lbPXQ/eOYmh6OqWniuJoejqlp0te4/ndQMC7ll2JzQjqm/Hgc66f0QICbjR4qJF0Zw76qy+sbrKGFhYUFQkNDax3qi42NRUREhN5eRwiB48ePw8PDQ2/bJCIiIiLTxRbtVF8G7RY4ffp0fPXVV/jmm29w+vRpTJs2DWlpaZg8eTKA6tP1Ro8erfGc48eP4/jx4ygqKkJ2djaOHz+OxMRE9eOzZ8/Gtm3bcP78eRw/fhwTJkzA8ePH1dskIiIiIroXtmin+jDoNVcjRoxAbm4u5syZg/T0dLRr1w5bt26Ft7c3gOpJg2+f8+rWrn/x8fH48ccf4e3tjQsXLgAA8vPz8eyzzyIjIwP29vYICQnB3r170bVr10Z7X0RERETU9NW0aH9s8QF1i/aFIzpBKq1/fwAybQZvaDFlyhRMmTKlzsdWrlxZa5kQd+/YMn/+fMyfP18fpRERERFRMxfgZoulz4Ri9DeHsOnEVfg4W+GVqEBDl0VGyqCnBRIRERERGbuIABe8/1h7AMBnO5Ox7sglA1dExorhioiIiIjoHoaHeeL5vv4AgP9t+BtxKbkGroiMEcMVEREREZEWXokMxEMdPFBeKTDpuyNIzioydElkZBiuiIiIiIi0wBbtdC8MV0REREREWmKLdrobhisiIiIiIh3UtGi3U5irW7RXVd29ozU1DwxXREREREQ6qmnRbi6VYNOJq5i//ayhSyIjwHBFRERERFQPbNFOt2O4IiIiIiKqp9tbtB9IyTFwRWRIDFdERERERPfh1hbtk7+LZ4v2ZozhioiIiIjoPrBFO9VguCIiIiIiuk81Ldq9nKzYor0ZY7giIiIiItIDZxs5vmGL9maN4YqIiIiISE8C3GzYor0ZY7giIiIiItIjtmhvvhiuiIiIiIj0jC3amyeGKyIiIiKiBsAW7c0PwxURERERUQNgi/bmh+GKiIiIiKiBsEV788JwRURERETUgNiivflguCIiIiIiamBs0d48MFwRERERETUCtmg3fQxXRERERESNhC3aTRvDFRERERFRI2KLdtPFcEVERERE1IjYot10MVwRERERETUytmg3TQxXREREREQGwBbtpofhioiIiIjIQALcbLB0FFu0mwqGKyIiIiIiA4rwd8Fctmg3CQxXREREREQG9gRbtJsEhisiIiIiIiPAFu1NH8MVEREREZERYIv2po/hioiIiIjISLBFe9PGcEVEREREZETYor3pYrgiIiIiIjIybNHeNDFcEREREREZIbZob3oYroiIiIiIjBRbtDctDFdEREREREaMLdqbDoYrIiIiIiIjxhbtTQfDFRERERGRkWOL9qaB4YqIiIiIqAm4vUX7q+tOsEW7kWG4IiIiIiJqIm5t0b45IR3zYtmi3ZgwXBERERERNSG3tmj/fBdbtBsThisiIiIioibmiTBPvNA3AABbtBsThisiIiIioiZoemQbtmg3MgxXRERERERNEFu0Gx+GKyIiIiKiJoot2o0LwxURERERURPGFu3Gg+GKiIiIiKiJY4t248BwRURERERkAtii3fAYroiIiIiITMStLdqj17NFe2NjuCIiIiIiMiE1LdorqtiivbExXBERERERmRC2aDcchisiIiIiIhPDFu2GwXBFRERERGSC2KK98TFcERERERGZKLZob1wMV0REREREJowt2hsPwxURERERkYlji/bGwXBFRERERNQMsEV7w2O4IiIiIiJqBmpatId6O7JFewNhuCIiIiIiaiYUMjMsHxXKFu0NhOGKiIiIiKgZYYv2hsNwRURERETUzLBFe8MweLhavHgxfH19oVAoEBoain379t1x3fT0dDz11FMIDAyEVCrF1KlT61zvl19+gUqlglwuh0qlwoYNGxqoeiIiIiKipokt2vXPoOFq7dq1mDp1Kt544w0cO3YMPXv2xODBg5GWllbn+mVlZXB1dcUbb7yBjh071rlOXFwcRowYgVGjRuHEiRMYNWoUhg8fjoMHDzbkWyEiIiIianLYol2/DBqu5s2bhwkTJmDixIkIDg7GggUL4OnpiSVLltS5vo+PDxYuXIjRo0fD3t6+znUWLFiAyMhIREdHIygoCNHR0ejfvz8WLFjQgO+EiIiIiKhpYot2/TE31AvfvHkT8fHxmDFjhsbyqKgoHDhwoN7bjYuLw7Rp0zSWDRw48K7hqqysDGVl/7ahLCwsBACUl5ejvLy83rXoS00NxlAL6QfH1PRwTE0Tx9X0cExNE8f1/n3wiApX82/gaFo+xq44hJ8ndYOztYXB6jGmMdWlBoOFq5ycHFRWVkKpVGosVyqVyMjIqPd2MzIydN7m3LlzMXv27FrLY2JiYGVlVe9a9C02NtbQJZCecUxND8fUNHFcTQ/H1DRxXO/P/7kBFzPNcPnaDYz8YheeV1VCZuAODcYwpiUlJVqva7BwVUMikWjcF0LUWtbQ24yOjsb06dPV9wsLC+Hp6YmoqCjY2dndVy36UF5ejtjYWERGRkImkxm6HNIDjqnp4ZiaJo6r6eGYmiaOq/6ERhRj+PKDSL1egV3FrTDvifaQSu/vd/P6MKYxrTmrTRsGC1cuLi4wMzOrdUQpKyur1pEnXbi7u+u8TblcDrlcXmu5TCYz+GDeytjqofvHMTU9HFPTxHE1PRxT08RxvX9BLRywdFQoRn99CFtOZsDX1QavDgw0WD3GMKa6vL7BDvRZWFggNDS01qG+2NhYRERE1Hu74eHhtbYZExNzX9skIiIiImou2KK9/gx6WuD06dMxatQohIWFITw8HMuXL0daWhomT54MoPp0vStXrmDVqlXq5xw/fhwAUFRUhOzsbBw/fhwWFhZQqVQAgJdffhm9evXChx9+iGHDhmHjxo3Yvn079u/f3+jvj4iIiIioKXoizBMXc0vw+a5kRK//Gy0dLRHh72LosoyeQcPViBEjkJubizlz5iA9PR3t2rXD1q1b4e3tDaB60uDb57wKCQlR/zs+Ph4//vgjvL29ceHCBQBAREQE1qxZgzfffBMzZ86Ev78/1q5di27dujXa+yIiIiIiauqmR7bBhdxibE5Ix+Tv4rF+Sg8EuNkYuiyjZvCGFlOmTMGUKVPqfGzlypW1lgkh7rnNxx9/HI8//vj9lkZERERE1GxJpRJ88kRHpBeUIv7iNYxbeQi/TukBZ5vavQqomoGbKxIRERERkbFSyMywfFQovJyscCnvBp79Lh6l5ZWGLstoMVwREREREdEdOdvI8c3YLrBTmCP+4jW8uu4EqqrufTZZc8RwRUREREREdxXgZoOlo0JhLpVgc0I65sWeNXRJRonhioiIiIiI7okt2u+N4YqIiIiIiLTyRJgnXugbAACIXv83DqTkGLgi48JwRUREREREWpse2QYPdfBARZXA5O/ikZxVZOiSjAbDFRERERERaa2mRXuotyMKSyswbuUh5BaVGboso8BwRUREREREOrm9Rft/Vh1hi3YwXBERERERUT3c2qL9aFo+W7SD4YqIiIiIiOqJLdo1MVwREREREVG9sUX7vxiuiIiIiIjovrBFezWGKyIiIiIium9s0c5wRUREREREesAW7QxXRERERESkJ829RTvDFRERERER6U1zbtHOcEVERERERHrVXFu0M1wREREREZHeNccW7QxXRERERETUIJpbi3aGKyIiIiIiajDNqUU7wxURERERETWY5tSineGKiIiIiIgaVHNp0c5wRUREREREDa45tGhnuCIiIiIiokZh6i3aGa6IiIiIiKjR3N6i/ScTatHOcEVERERERI3q1hbt/zOhFu0MV0RERERE1OhMsUU7wxURERERETU6U2zRznBFREREREQGUVeL9uKyChxMzUN8jgQHU/NQ2YQ6CpobugAiIiIiImq+alq0P7b4TxxNy0fou7EoLa8CYIZV547Aw16BWUNVGNTOw9Cl3hOPXBERERERkUEFuNlg/AO+APBPsPpXRkEpnvv+KP44mW6I0nTCcEVERERERAZVWSWw9nDdLdlrTgqcvSnR6E8RZLgiIiIiIiKDOpSah/SC0js+LgCkF5TiUGpe4xVVDwxXRERERERkUFnX7xys6rOeoTBcERERERGRQbnZKvS6nqEwXBERERERkUF19XWCh70Ckjs8LgHgYa9AV1+nxixLZwxXRERERERkUGZSCWYNVQFArYBVc3/WUBXMpHeKX8aB4YqIiIiIiAxuUDsPLHmmM9ztNU/9c7dXYMkznZvEPFecRJiIiIiIiIzCoHYeiFS5Iy45CzH7DiKqZzeEB7gZ/RGrGgxXRERERERkNMykEnTzdULuaYFuvk5NJlgBPC2QiIiIiIhILxiuiIiIiIiI9IDhioiIiIiISA8YroiIiIiIiPSA4YqIiIiIiEgPGK6IiIiIiIj0gOGKiIiIiIhIDxiuiIiIiIiI9IDhioiIiIiISA8YroiIiIiIiPSA4YqIiIiIiEgPGK6IiIiIiIj0gOGKiIiIiIhID8wNXYAxEkIAAAoLCw1cSbXy8nKUlJSgsLAQMpnM0OWQHnBMTQ/H1DRxXE0Px9Q0cVxNjzGNaU0mqMkId8NwVYfr168DADw9PQ1cCRERERERGYPr16/D3t7+rutIhDYRrJmpqqrC1atXYWtrC4lEYuhyUFhYCE9PT1y6dAl2dnaGLof0gGNqejimponjano4pqaJ42p6jGlMhRC4fv06WrRoAan07ldV8chVHaRSKVq1amXoMmqxs7Mz+JeL9Itjano4pqaJ42p6OKamieNqeoxlTO91xKoGG1oQERERERHpAcMVERERERGRHjBcNQFyuRyzZs2CXC43dCmkJxxT08MxNU0cV9PDMTVNHFfT01THlA0tiIiIiIiI9IBHroiIiIiIiPSA4YqIiIiIiEgPGK6IiIiIiIj0gOGKiIiIiIhIDxiuDGzv3r0YOnQoWrRoAYlEgl9//fWez9mzZw9CQ0OhUCjg5+eHpUuXNnyhpDVdx3T37t2QSCS1bklJSY1TMN3T3Llz0aVLF9ja2sLNzQ2PPPIIzpw5c8/ncV81bvUZV+6vxm3JkiXo0KGDetLR8PBw/P7773d9DvdT46fruHI/bXrmzp0LiUSCqVOn3nW9prC/MlwZWHFxMTp27IjPP/9cq/VTU1Px4IMPomfPnjh27Bj+97//4aWXXsIvv/zSwJWStnQd0xpnzpxBenq6+ta6desGqpB0tWfPHjz//PP466+/EBsbi4qKCkRFRaG4uPiOz+G+avzqM641uL8ap1atWuGDDz7AkSNHcOTIEfTr1w/Dhg3DqVOn6lyf+2nToOu41uB+2jQcPnwYy5cvR4cOHe66XpPZXwUZDQBiw4YNd13ntddeE0FBQRrLJk2aJLp3796AlVF9aTOmu3btEgDEtWvXGqUmun9ZWVkCgNizZ88d1+G+2vRoM67cX5seR0dH8dVXX9X5GPfTputu48r9tOm4fv26aN26tYiNjRW9e/cWL7/88h3XbSr7K49cNTFxcXGIiorSWDZw4EAcOXIE5eXlBqqK9CEkJAQeHh7o378/du3aZehy6C4KCgoAAE5OTndch/tq06PNuNbg/mr8KisrsWbNGhQXFyM8PLzOdbifNj3ajGsN7qfG7/nnn8eQIUMwYMCAe67bVPZXc0MXQLrJyMiAUqnUWKZUKlFRUYGcnBx4eHgYqDKqLw8PDyxfvhyhoaEoKyvDd999h/79+2P37t3o1auXocuj2wghMH36dDzwwANo167dHdfjvtq0aDuu3F+N399//43w8HCUlpbCxsYGGzZsgEqlqnNd7qdNhy7jyv20aVizZg2OHj2Kw4cPa7V+U9lfGa6aIIlEonFfCFHncmoaAgMDERgYqL4fHh6OS5cu4ZNPPuH/BIzQCy+8gISEBOzfv/+e63JfbTq0HVfur8YvMDAQx48fR35+Pn755ReMGTMGe/bsueMv4txPmwZdxpX7qfG7dOkSXn75ZcTExEChUGj9vKawv/K0wCbG3d0dGRkZGsuysrJgbm4OZ2dnA1VF+ta9e3ecO3fO0GXQbV588UX89ttv2LVrF1q1anXXdbmvNh26jGtduL8aFwsLCwQEBCAsLAxz585Fx44dsXDhwjrX5X7adOgyrnXhfmpc4uPjkZWVhdDQUJibm8Pc3Bx79uzBokWLYG5ujsrKylrPaSr7K49cNTHh4eHYtGmTxrKYmBiEhYVBJpMZqCrSt2PHjhnN4W2q/svYiy++iA0bNmD37t3w9fW953O4rxq/+oxrXbi/GjchBMrKyup8jPtp03W3ca0L91Pj0r9/f/z9998ay8aNG4egoCC8/vrrMDMzq/WcJrO/GqqTBlW7fv26OHbsmDh27JgAIObNmyeOHTsmLl68KIQQYsaMGWLUqFHq9c+fPy+srKzEtGnTRGJiovj666+FTCYTP//8s6HeAt1G1zGdP3++2LBhgzh79qw4efKkmDFjhgAgfvnlF0O9BbrNc889J+zt7cXu3btFenq6+lZSUqJeh/tq01OfceX+atyio6PF3r17RWpqqkhISBD/+9//hFQqFTExMUII7qdNla7jyv20abq9W2BT3V8Zrgyspl3o7bcxY8YIIYQYM2aM6N27t8Zzdu/eLUJCQoSFxf+3d/dBUVV9HMC/F0xX2QhEWiCIVRMCRMggDUxUMOhFdCaFcUgkZmpQQGdcaarJAW1UlqEsDZlqhpepRogUtVFWUhELR9YgwgRWJRg0QTMsx/oD2D3PHz7cvPLusz0qfj8zzHDOPfd3zj13r+7Pe891rNBqtSI3N/f/P3Aa0EjPqV6vF1OnThUqlUo4OjqKOXPmiAMHDtydwVO/+jufAER+fr7chtfq/edOziuv13tbYmKi8PT0FGPHjhXOzs4iPDxc/gIuBK/T+9VIzyuv0/vT7cnV/Xq9SkL8dyUYERERERER3TG+0IKIiIiIiMgKmFwRERERERFZAZMrIiIiIiIiK2ByRUREREREZAVMroiIiIiIiKyAyRUREREREZEVMLkiIiIiIiKyAiZXREREREREVsDkioiI/hWtra2QJAl1dXV3eyiypqYmzJ49GyqVCoGBgXd7OERENMowuSIiGqUSEhIgSRIyMzMV9Xv37oUkSXdpVHdXeno67OzsYDKZcOTIkQHbdXR0IDU1FVOmTMG4cePg4eGBRYsWDbrPgyghIQFLliy528MgIrpnMLkiIhrFVCoV9Ho9rl27dreHYjVdXV13vG9zczPmzJkDT09PODk59dumtbUVTz/9NI4ePYqsrCycPn0aBoMB8+fPR3Jy8h33TUREox+TKyKiUSwiIgIuLi7YunXrgG0yMjL6PCL34YcfQqvVyuXeOxRbtmyBRqOBg4MDNm7ciJ6eHqSlpWHixIlwd3dHXl5en/hNTU0ICQmBSqWCn58fjh07ptje0NCAF198EWq1GhqNBitWrMDVq1fl7fPmzUNKSgrWrVuHSZMmYeHChf0eh8ViwaZNm+Du7o5x48YhMDAQBoNB3i5JEmpqarBp0yZIkoSMjIx+46xevRqSJMFoNGLp0qXw8vKCn58f1q1bh5MnT8rt2trasHjxYqjVatjb2yMmJgaXL1/uM695eXl4/PHHoVarsWrVKpjNZmRlZcHFxQWPPvooNm/erOhfkiTk5ubihRdewPjx4zF58mSUlJQo2pw+fRoLFizA+PHj4eTkhDfeeAM3btzoc76ys7Ph6uoKJycnJCcno7u7W27T1dWFN998E4899hjs7Owwa9YsxbkpKCiAg4MDDh06BB8fH6jVakRFRaG9vV0+vsLCQuzbtw+SJEGSJBw7dgxdXV1ISUmBq6srVCoVtFrtoJ8/IqLRhMkVEdEoZmtriy1btmDHjh24ePHi/xTr6NGjuHTpEo4fP44PPvgAGRkZePnll+Ho6Ijq6mokJSUhKSkJFy5cUOyXlpYGnU6HH3/8ESEhIYiOjsbvv/8OAGhvb0dYWBgCAwPxww8/wGAw4PLly4iJiVHEKCwsxJgxY1BVVYVPPvmk3/F99NFHeP/995GdnY36+npERkYiOjoa586dk/vy8/ODTqdDe3s71q9f3ydGZ2cnDAYDkpOTYWdn12e7g4MDAEAIgSVLlqCzsxOVlZX49ttv0dzcjNjYWEX75uZmlJWVwWAwYNeuXcjLy8NLL72EixcvorKyEnq9Hu+++64iaQOADRs24JVXXsFPP/2EV199FcuXL0djYyMA4O+//0ZUVBQcHR1x6tQplJSU4PDhw0hJSVHEqKioQHNzMyoqKlBYWIiCggIUFBTI21977TVUVVWhqKgI9fX1WLZsGaKiouT56u0rOzsbn3/+OY4fP462tjZ53tavX4+YmBg54Wpvb0dISAi2b9+O/fv346uvvoLJZMIXX3yhSNSJiEY1QUREo9LKlSvF4sWLhRBCzJ49WyQmJgohhCgtLRW3/vGfnp4uAgICFPtu27ZNeHp6KmJ5enoKs9ks13l7e4vnnntOLvf09Ag7Ozuxa9cuIYQQLS0tAoDIzMyU23R3dwt3d3eh1+uFEEJs2LBBPP/884q+L1y4IAAIk8kkhBAiLCxMBAYGDnm8bm5uYvPmzYq64OBgsXr1arkcEBAg0tPTB4xRXV0tAIg9e/YM2ld5ebmwtbUVbW1tct2ZM2cEAGE0GoUQN+d1woQJ4vr163KbyMhIodVq+8zj1q1b5TIAkZSUpOhv1qxZYtWqVUIIIT799FPh6Ogobty4IW8/cOCAsLGxER0dHUKIf85XT0+P3GbZsmUiNjZWCCHE+fPnhSRJ4tdff1X0Ex4eLt5++20hhBD5+fkCgDh//ry8PScnR2g0Grl862esV2pqqliwYIGwWCwDzh8R0WjFO1dERA8AvV6PwsJCNDQ03HEMPz8/2Nj889eGRqOBv7+/XLa1tYWTkxOuXLmi2O/ZZ5+Vfx8zZgyCgoLkuzA1NTWoqKiAWq2Wf5588kkAN+/69AoKChp0bNevX8elS5cQGhqqqA8NDZX7Gg4hBAAM+cKPxsZGeHh4wMPDQ67z9fWFg4ODoj+tVouHH35YLms0Gvj6+vaZx8HmrLfcG7exsREBAQGKO2uhoaGwWCwwmUxynZ+fH2xtbeWyq6ur3E9tbS2EEPDy8lLMfWVlpWLeJ0yYgKlTp/YbYyAJCQmoq6uDt7c31qxZg/Ly8kHbExGNJmPu9gCIiOjfN3fuXERGRuKdd95BQkKCYpuNjY2cVPS6dW1Or4ceekhRliSp3zqLxTLkeHqTF4vFgkWLFkGv1/dp4+rqKv/e3yN6g8XtJYQY0ZsRp02bBkmS0NjYOOhb8AaKe3v9vzFngx3TUH339mOxWGBra4uamhpFAgYAarV60Bi3f1ZuN3PmTLS0tKCsrAyHDx9GTEwMIiIi8PXXXw9xhERE9z/euSIiekBkZmbim2++wYkTJxT1zs7O6OjoUHxptub/TXXreqKenh7U1NTId6dmzpyJM2fOQKvV4oknnlD8DDehAgB7e3u4ubnh+++/V9SfOHECPj4+w44zceJEREZGIicnB3/99Vef7X/88QeAm3ep2traFOvLGhoa8Oeff46ov4Hcvgbr5MmT8pz5+vqirq5OMb6qqirY2NjAy8trWPGfeuopmM1mXLlypc+8u7i4DHucY8eOhdls7lNvb2+P2NhYfPbZZyguLsbu3bvR2dk57LhERPcrJldERA8If39/xMXFYceOHYr6efPm4bfffkNWVhaam5uRk5ODsrIyq/Wbk5OD0tJSNDU1ITk5GdeuXUNiYiIAIDk5GZ2dnVi+fDmMRiN++eUXlJeXIzExsd8v7YNJS0uDXq9HcXExTCYT3nrrLdTV1WHt2rUjirNz506YzWY888wz2L17N86dO4fGxkZs375dflwvIiICM2bMQFxcHGpra2E0GhEfH4+wsLAhH2EcjpKSEuTl5eHs2bNIT0+H0WiUX1gRFxcHlUqFlStX4ueff0ZFRQVSU1OxYsUKaDSaYcX38vJCXFwc4uPjsWfPHrS0tODUqVPQ6/U4ePDgsMep1WpRX18Pk8mEq1evoru7G9u2bUNRURGamppw9uxZlJSUwMXFRX4ZCBHRaMbkiojoAfLee+/1eazLx8cHO3fuRE5ODgICAmA0Gvt9k96dyszMhF6vR0BAAL777jvs27cPkyZNAgC4ubmhqqoKZrMZkZGRmD59OtauXYtHHnlEsS5pONasWQOdTgedTgd/f38YDAbs378f06ZNG1GcyZMno7a2FvPnz4dOp8P06dOxcOFCHDlyBLm5uQBuPh63d+9eODo6Yu7cuYiIiMCUKVNQXFw8or4GsnHjRhQVFWHGjBkoLCzEl19+CV9fXwA310EdOnQInZ2dCA4OxtKlSxEeHo6PP/54RH3k5+cjPj4eOp0O3t7eiI6ORnV1tWId2VBef/11eHt7IygoCM7OzqiqqoJarYZer0dQUBCCg4PR2tqKgwcPjvh8EhHdjyQx1MPTRERE9H8jSRJKS0sHXfNFRET3Jv4zEhERERERkRUwuSIiIiIiIrICvoqdiIjoHsKn9YmI7l+8c0VERERERGQFTK6IiIiIiIisgMkVERERERGRFTC5IiIiIiIisgImV0RERERERFbA5IqIiIiIiMgKmFwRERERERFZAZMrIiIiIiIiK/gPOqH8z/837s0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# See how many dimensions to set by PCA\n",
    "pca = PCA()\n",
    "pca.fit(X_train)\n",
    "\n",
    "# Plot explained variance ratio to figure out how many dimensions to keep\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(range(1, len(pca.explained_variance_ratio_) + 1), pca.explained_variance_ratio_, marker='o', linestyle='-')\n",
    "plt.xlabel('Number of Components')\n",
    "plt.ylabel('Explained Variance Ratio')\n",
    "plt.title('Explained Variance Ratio vs. Number of Components')\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "b33a7d00-c1be-4b3d-879a-7815a0d9e5e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# PCA-reduce dimension\n",
    "pca = PCA(n_components=3)\n",
    "X_train_pca = pca.fit_transform(X_train)\n",
    "X_test_pca = pca.transform(X_test)\n",
    "\n",
    "# X_train_s and y_train_s are 0.2 sized training set from 2017\n",
    "# X_val and y_val are 0.8 sized validation set also 2017\n",
    "# X_test and y_test are full-sized test set 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7be9908d-ab21-4dc9-a552-81244e691c91",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simply fit into model\n",
    "svr = SVR()\n",
    "print(\"fitting\")\n",
    "svr.fit(X_train_pca, y_train)\n",
    "print(\"fitting done!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "5f68756e-c577-40e4-99b9-4888a3a15227",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "score on training set acquired!\n",
      "score on validation set acquired!\n",
      "SVR Train Score: 0.17670973763354014\n",
      "SVR Validation Score: 0.15826277503819675\n"
     ]
    }
   ],
   "source": [
    "# Evaluate SVR model on r squared value\n",
    "train_score = svr.score(X_train_pca, y_train)\n",
    "print('score on training set acquired!')\n",
    "test_score = svr.score(X_test_pca, y_test)\n",
    "print('score on test set acquired!')\n",
    "\n",
    "print(\"SVR Train Score:\", train_score)\n",
    "print(\"SVR Test Score:\", test_score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1b862a7-a88b-40ba-b925-e5bb782c3d40",
   "metadata": {},
   "outputs": [],
   "source": [
    "###############################################################\n",
    "## Option: If try grid search best params\n",
    "# Probably should to use a smaller training set here\n",
    "param_grid = {'C': [0.1, 1, 10]}\n",
    "\n",
    "svr_model = SVR(kernel='rbf')\n",
    "\n",
    "grid_svr = GridSearchCV(estimator=svr_model, param_grid=param_grid, cv=3, verbose=1)\n",
    "grid_svr.fit(X_train_pca, y_train_s)\n",
    "\n",
    "best_params = grid_svr.best_params_\n",
    "best_score = grid_svr.best_score_\n",
    "\n",
    "print(\"Best Parameters:\", best_params)\n",
    "print(\"Best Score:\", best_score)\n",
    "\n",
    "# Fit the training data with best parameters\n",
    "best_svr = grid_svr.best_estimator_\n",
    "best_svr.fit(X_train_pca, y_train_s)\n",
    "\n",
    "\n",
    "# Evaluate SVR model on r squared value\n",
    "train_score = svr.score(X_train_pca, y_train_s)\n",
    "print('score on training set acquired!')\n",
    "test_score = svr.score(X_val_pca, y_val)\n",
    "print('score on validation set acquired!')\n",
    "\n",
    "print(\"SVR Train Score:\", train_score)\n",
    "print(\"SVR Validation Score:\", test_score)\n",
    "\n",
    "# takes around 10 mins for 0.2-sized tram only 2017"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "4f64e0e5-d9ed-4276-87dd-d043356952d4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['models/svr_model_2024-05-15_10-55-34.joblib']"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Dump save trained model\n",
    "timestamp = datetime.now().strftime(\"%Y-%m-%d_%H-%M-%S\")\n",
    "directory_path = 'models/'\n",
    "\n",
    "svr_model_filename = f'{directory_path}svr_model_{timestamp}.joblib'\n",
    "dump(svr, svr_model_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4660f508-83c2-44bc-9f98-d0f4e811fdd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Check the metrics of svr on test set\n",
    "# Predict on the test set\n",
    "X_test_pca = pca.fit_transform(X_test)\n",
    "y_pred_svr = svr.predict(X_test_pca)\n",
    "\n",
    "# Reverse scaling\n",
    "testset_svr = testset.copy()\n",
    "testset_svr['Occupancy'] = y_pred_svr\n",
    "\n",
    "# Reverse scale the y\n",
    "testset_re = scaler.inverse_transform(testset[columns_to_normalize])\n",
    "testset_svr = scaler.inverse_transform(testset_svr[columns_to_normalize])\n",
    "\n",
    "y_test_re = testset_re[:, 3]\n",
    "y_pred_svr = testset_svr[:, 3]\n",
    "\n",
    "mape_svr_te = np.mean(np.abs((y_test_re - y_pred_svr) / y_test_re + 1e-10)) * 100\n",
    "mae_svr_te = mean_absolute_error(y_test_re, y_pred_svr)\n",
    "mse_svr_te = mean_squared_error(y_test_re, y_pred_svr)\n",
    "\n",
    "print(\"SVR MAPE on test:\", mape_svr_te)\n",
    "print(\"SVR MAE on test:\", mae_svr_te)\n",
    "print(\"SVR MSE on test:\", mse_svr_te)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09508036-6289-483f-b6bb-74ccf3127481",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Check the metrics of svr on training set\n",
    "# Predict on the training set\n",
    "X_train_pca = pca.fit_transform(X_train) # fullset\n",
    "y_pred_svr = svr.predict(X_train_pca)\n",
    "\n",
    "mape_svr_tr = np.mean(np.abs((y_train - y_pred_svr) / y_test + 1e-10)) * 100\n",
    "mae_svr_tr = mean_absolute_error(y_train, y_pred_svr)\n",
    "mse_svr_tr = mean_squared_error(y_train, y_pred_svr)\n",
    "\n",
    "print(\"SVR MAPE on train:\", mape_svr_tr)\n",
    "print(\"SVR MAE on train:\", mae_svr_tr)\n",
    "print(\"SVR MSE on train:\", mse_svr_tr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16b1a1dc-07f7-43f0-bc46-7107fa727cf7",
   "metadata": {},
   "source": [
    "### Discarded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "7791d549-7c4d-47fa-92f6-d73060d15768",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Linear Regression MAPE on test set: 201.64348264511864\n",
      "Linear Regression MAE on test set: 0.7647051788651441\n",
      "Linear Regression MSE on test set: 1.1195315325049027\n",
      "Random Forest MAPE on test set: 477.58859257554855\n",
      "Random Forest MAE on test set: 0.613925827698001\n",
      "Random Forest MSE on test set: 0.8019942406659941\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:    1.8s finished\n"
     ]
    }
   ],
   "source": [
    "# Check the metrics for predicting on the test set without reverse scaling, for lr and rf\n",
    "# Predict on the test set\n",
    "y_pred_lr = best_lr.predict(X_test)\n",
    "y_pred_rf = best_rf.predict(X_test)\n",
    "\n",
    "# Calculate MAPE with a small denominator to prevent inf\n",
    "mape_lr = np.mean(np.abs((y_test - y_pred_lr) / (y_test + 1e-10))) * 100\n",
    "mape_rf = np.mean(np.abs((y_test - y_pred_rf) / (y_test + 1e-10))) * 100\n",
    "\n",
    "# Calculate MAE \n",
    "mae_lr = mean_absolute_error(y_test, y_pred_lr)\n",
    "mae_rf = mean_absolute_error(y_test, y_pred_rf)\n",
    "\n",
    "# Calculate MSE\n",
    "mse_lr = mean_squared_error(y_test, y_pred_lr)\n",
    "mse_rf = mean_squared_error(y_test, y_pred_rf)\n",
    "\n",
    "print(\"Linear Regression MAPE on test set:\", mape_lr)\n",
    "print(\"Linear Regression MAE on test set:\", mae_lr)\n",
    "print(\"Linear Regression MSE on test set:\", mse_lr)\n",
    "\n",
    "print(\"Random Forest MAPE on test set:\", mape_rf)\n",
    "print(\"Random Forest MAE on test set:\", mae_rf)\n",
    "print(\"Random Forest MSE on test set:\", mse_rf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "cb1b8279-0ab5-4949-9859-e807510366a2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Linear Regression MAPE on training set: 183.02961228110286\n",
      "Linear Regression MAE on training set: 0.7149038139716958\n",
      "Linear Regression MSE on training set: 0.9784481069489978\n",
      "Random Forest MAPE on training set: 384.440135763696\n",
      "Random Forest MAE on training set: 0.5676781537954502\n",
      "Random Forest MSE on training set: 0.6903329825337616\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done 200 out of 200 | elapsed:    1.2s finished\n"
     ]
    }
   ],
   "source": [
    "# Check the metrics for predicting on the training set without reverse scaling, for lr and rf\n",
    "y_pred_lr = best_lr.predict(X_train)\n",
    "y_pred_rf = best_rf.predict(X_train)\n",
    "\n",
    "# Calculate MAPE with a small denominator to prevent inf\n",
    "mape_lr = np.mean(np.abs((y_train - y_pred_lr) / (y_train + 1e-10))) * 100\n",
    "mape_rf = np.mean(np.abs((y_train - y_pred_rf) / (y_train + 1e-10))) * 100\n",
    "\n",
    "# Calculate MAE \n",
    "mae_lr = mean_absolute_error(y_train, y_pred_lr)\n",
    "mae_rf = mean_absolute_error(y_train, y_pred_rf)\n",
    "\n",
    "# Calculate MSE\n",
    "mse_lr = mean_squared_error(y_train, y_pred_lr)\n",
    "mse_rf = mean_squared_error(y_train, y_pred_rf)\n",
    "\n",
    "print(\"Linear Regression MAPE on training set:\", mape_lr)\n",
    "print(\"Linear Regression MAE on training set:\", mae_lr)\n",
    "print(\"Linear Regression MSE on training set:\", mse_lr)\n",
    "\n",
    "print(\"Random Forest MAPE on training set:\", mape_rf)\n",
    "print(\"Random Forest MAE on training set:\", mae_rf)\n",
    "print(\"Random Forest MSE on training set:\", mse_rf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d84e9e22-13d9-4533-90a0-6263b309f7e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Findout a good dimension number for UMAP <- takes forever, discard\n",
    "dimensions_to_try = [2]\n",
    "\n",
    "# Plot each dimensionality reduction\n",
    " for n_components in dimensions_to_try:\n",
    "     Reduce dimensionality using UMAP\n",
    "     reducer = umap.UMAP(n_components=n_components)\n",
    "     X_umap = reducer.fit_transform(X_train)\n",
    "    \n",
    "     plt.figure(figsize=(8, 6))\n",
    "     if n_components == 2:\n",
    "        plt.scatter(X_umap[:, 0], X_umap[:, 1], alpha=0.5)\n",
    "        plt.title(\"UMAP (2D)\")\n",
    "     elif n_components == 3:\n",
    "        ax = plt.axes(projection='3d')\n",
    "        ax.scatter3D(X_umap[:, 0], X_umap[:, 1], X_umap[:, 2], alpha=0.5)\n",
    "        ax.set_title(\"UMAP (3D)\")\n",
    "     elif n_components == 4:\n",
    "        ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "197b06cc-e4aa-4d53-b433-55edb61ef73a",
   "metadata": {},
   "outputs": [],
   "source": [
    "## SVR\n",
    "# Try with SVR with UMAP dimension reduction <- Always fails at spectral initialization\n",
    "# Reduce dimensionality using UMAP\n",
    "reducer = umap.UMAP(n_components=2)\n",
    "X_train_umap = reducer.fit_transform(X_train)\n",
    "X_test_umap = reducer.transform(X_test)\n",
    "\n",
    "\n",
    "svr = SVR()\n",
    "svr.fit(X_train_umap, y_train)\n",
    "\n",
    "# Evaluate SVR model on r squared value\n",
    "train_score = svr.score(X_train_umap, y_train)\n",
    "test_score = svr.score(X_test_umap, y_test)\n",
    "\n",
    "print(\"SVR Train Score:\", train_score)\n",
    "print(\"SVR Test Score:\", test_score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb47b995-b29e-40b5-97ee-3cdc4367bfc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Check the metrics of svr on test set\n",
    "X_test_pca = pca.fit_transform(X_test)\n",
    "# Predict on the test set\n",
    "y_pred_svr = svr.predict(X_test_pca)\n",
    "\n",
    "mape_svr_te = np.mean(np.abs((y_test - y_pred_svr) / y_test + 1e-10)) * 100\n",
    "mae_svr_te = mean_absolute_error(y_test, y_pred_svr)\n",
    "mse_svr_te = mean_squared_error(y_test, y_pred_svr)\n",
    "\n",
    "print(\"SVR MAPE on test:\", mape_svr_te)\n",
    "print(\"SVR MAE on test:\", mae_svr_te)\n",
    "print(\"SVR MSE on test:\", mse_svr_te)\n",
    "\n",
    "## Check the metrics of svr on training set\n",
    "# Predict on the training set\n",
    "X_train_pca = pca.fit_transform(X_train) # fullset\n",
    "y_pred_svr = svr.predict(X_train_pca)\n",
    "\n",
    "mape_svr_tr = np.mean(np.abs((y_train - y_pred_svr) / y_test + 1e-10)) * 100\n",
    "mae_svr_tr = mean_absolute_error(y_train, y_pred_svr)\n",
    "mse_svr_tr = mean_squared_error(y_train, y_pred_svr)\n",
    "\n",
    "print(\"SVR MAPE on train:\", mape_svr_tr)\n",
    "print(\"SVR MAE on train:\", mae_svr_tr)\n",
    "print(\"SVR MSE on train:\", mse_svr_tr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "4af7737a-9229-4f9d-840c-e896b4944797",
   "metadata": {},
   "outputs": [],
   "source": [
    "# If reload model\n",
    "import joblib\n",
    "best_lr = joblib.load('../models/linear_regression_model_2024-05-14_22-18-31.joblib')\n",
    "best_rf = joblib.load('../models/random_forest_model_2024-05-14_22-18-31.joblib')\n",
    "svr = joblib.load('../models/svr_model_2024-05-15_10-55-34.joblib')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b4fff99-52c9-4083-a2de-5f7c1d57307b",
   "metadata": {},
   "source": [
    "## Remarks:\n",
    "Change to train with 2017 data to predict 2019 data for better similarity in volume with vehicle types (planned to use 2022/2023)\n",
    "see https://github.com/Chocobosaurus/uzh-data-science-project/blob/main/WZ/Q1Q2.md"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f664f70-9dd0-4cfc-b04e-0e9125b27e2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Models:\n",
    "# 2024-05-14_15-59-58: Done with unormalized Occupancy, with vehicle type. lr + rf\n",
    "# 2024-05-14_16-49-07: Done with normalized Occupancy, with vehicle type. lr + rf\n",
    "# 2024-05-14_17-16-45: Done with normalized Occupanct, 2017 as training set. lr + rf\n",
    "# 2024-05-14_20-29-31: 0.2-sized tram only 2017, svr.\n",
    "# 2024-05-14_21-57-38: 0.2-sized tram only 2017, svr with grid search. C=10.\n",
    "# 2024-05-14_22-18-31: tram only 2017, lr + rf.\n",
    "# model_2024-05-15_10-55-34: tram only 2017, svr."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2231af8a-e03a-4053-a258-cb9943948a6f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}