diff --git a/a_ddm_feat_importance.ipynb b/a_ddm_feat_importance.ipynb new file mode 100644 index 0000000..8f59762 --- /dev/null +++ b/a_ddm_feat_importance.ipynb @@ -0,0 +1,319 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "irish-antibody", + "metadata": {}, + "source": [ + "# DDM: FEATURE IMPORTANCE" + ] + }, + { + "cell_type": "markdown", + "id": "renewable-disorder", + "metadata": {}, + "source": [ + "# 1. Dataset Definition" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "congressional-screw", + "metadata": {}, + "outputs": [], + "source": [ + "import yaml\n", + "\n", + "# SELECT CONFIG FILE TO READ DATA FROM\n", + "config_file = \"conf//data//cartpole_st_at.yaml\"" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "partial-insertion", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "DATA STRUCTURE SELECTED:\n", + " - input_cols: ['state_x_position', 'state_x_velocity', 'state_angle_position', 'state_angle_velocity']\n", + " - augmented_cols: ['config_length', 'config_masspole', 'action_command']\n", + " - output_cols: ['state_x_position', 'state_x_velocity', 'state_angle_position', 'state_angle_velocity']\n" + ] + } + ], + "source": [ + "# For quick prototyping we will test on the fitted model. If model is not even good, then we need to reiterate\n", + "##from sklearn.model_selection import train_test_split\n", + "##x_train, x_test, y_train, y_test = train_test_split(x_set, y_set, test_size=0.33, random_state=42)\n", + "\n", + "# Open and extract sim config from YAML file\n", + "with open(config_file, 'r') as file:\n", + " data_config = yaml.load(file, Loader=yaml.FullLoader)\n", + "\n", + "# Extract features from yaml file\n", + "input_cols = data_config['data']['inputs']\n", + "output_cols = data_config['data']['outputs']\n", + "augmented_cols = data_config['data']['augmented_cols']\n", + "dataset_path = data_config['data']['path']\n", + "iteration_order = data_config['data']['iteration_order']\n", + "episode_col = data_config['data']['episode_col']\n", + "iteration_col = data_config['data']['iteration_col']\n", + "max_rows = data_config['data']['max_rows']\n", + "\n", + "\n", + "config_file = \"conf//model//xgboost.yaml\"\n", + "# Open and extract sim config from YAML file\n", + "with open(config_file, 'r') as file:\n", + " model_config = yaml.load(file, Loader=yaml.FullLoader)\n", + "\n", + "print(\"DATA STRUCTURE SELECTED:\")\n", + "print(\" - input_cols:\", input_cols)\n", + "print(\" - augmented_cols:\", augmented_cols)\n", + "print(\" - output_cols:\", output_cols)" + ] + }, + { + "cell_type": "markdown", + "id": "usual-relation", + "metadata": {}, + "source": [ + "# 2. Model Definition" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "hydraulic-albany", + "metadata": {}, + "outputs": [], + "source": [ + "import gboost_models\n", + "import skmodels\n", + "\n", + "GBoostModel = gboost_models.GBoostModel\n", + "SKModel = skmodels.SKModel\n", + "\n", + "# SELECT MODEL TO BE USED\n", + "model_to_be_used = 3" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "canadian-supplement", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Model selected 'SVR' --> scale_data == False\n" + ] + } + ], + "source": [ + "if model_to_be_used == 1:\n", + " # XGBOOST\n", + " save_path = \"models/boost/moab/xgboost_model\"\n", + " model_name = \"xgboost\"\n", + " scale_data = False\n", + " print(f\"Model selected '{model_name}' --> scale_data == {scale_data}\")\n", + "\n", + " model = GBoostModel()\n", + " \n", + "elif model_to_be_used == 2:\n", + " # LIGHTGBM\n", + " save_path = \"models/boost/moab/lightgbm\"\n", + " model_name = \"lightgbm\"\n", + " scale_data = False\n", + " print(f\"Model selected '{model_name}' --> scale_data == {scale_data}\")\n", + "\n", + " model = GBoostModel()\n", + " \n", + "elif model_to_be_used == 3:\n", + " # SVR\n", + " save_path = \"models/sklearn/moab/SVR\"\n", + " model_name = \"SVR\"\n", + " scale_data = False\n", + " print(f\"Model selected '{model_name}' --> scale_data == {scale_data}\")\n", + "\n", + " model = SKModel()\n", + " \n", + "elif model_to_be_used == 4:\n", + " # GRADIENT BOOSTING\n", + " save_path = \"models/sklearn/moab/GradientBoostingRegressor\"\n", + " model_name = \"GradientBoostingRegressor\"\n", + " scale_data = False\n", + " print(f\"Model selected '{model_name}' --> scale_data == {scale_data}\")\n", + "\n", + " model = SKModel()\n", + " \n", + "elif model_to_be_used == 5:\n", + " # LINEAR MODEL\n", + " save_path = \"models/sklearn/moab/linear_model\"\n", + " model_name = \"linear_model\"\n", + " scale_data = True\n", + " print(f\"Model selected '{model_name}' --> scale_data == {scale_data}\")\n", + "\n", + " model = SKModel()\n" + ] + }, + { + "cell_type": "markdown", + "id": "circular-navigation", + "metadata": {}, + "source": [ + "# 3. Model Definition, Fitting, & Saving" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "analyzed-prospect", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:data_loaders:Iteration order set to 1 so using outputs from next 1 row\n", + "INFO:data_loaders:Next states are being added to same row with prefix next\n", + "INFO:data_loaders:Feature columns are: ['next_state_x_position', 'next_state_x_velocity', 'next_state_angle_position', 'next_state_angle_velocity', 'next_config_length', 'next_config_masspole', 'next_action_command']\n", + "INFO:data_loaders:Label columns are: ['state_x_position', 'state_x_velocity', 'state_angle_position', 'state_angle_velocity']\n", + "INFO:skmodels:Note: fit_separate should be True for SVR. Changing to True ..\n", + "INFO:skmodels:Fitting model 1 of 4\n", + "INFO:skmodels:Fitting model 2 of 4\n", + "INFO:skmodels:Fitting model 3 of 4\n", + "INFO:skmodels:Fitting model 4 of 4\n" + ] + } + ], + "source": [ + "\n", + "X, y = model.load_csv(\n", + " input_cols=input_cols,\n", + " output_cols=output_cols,\n", + " augm_cols=augmented_cols,\n", + " dataset_path=dataset_path,\n", + " iteration_order=iteration_order,\n", + " episode_col=episode_col,\n", + " iteration_col=iteration_col,\n", + " max_rows=max_rows,\n", + ")\n", + "#logger.info(\"Building model...\")\n", + "model.build_model(model_type=model_name, scale_data=scale_data)\n", + "#logger.info(\"Fitting model...\")\n", + "model.fit(X, y, fit_separate=False)\n", + "\n", + "#logger.info(f\"Saving model to {save_path}\")\n", + "model.save_model(filename=save_path)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "republican-wound", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "DATA STRUCTURE FROM UPLOADED MODEL:\n", + " - input_cols: ['state_x_position', 'state_x_velocity', 'state_angle_position', 'state_angle_velocity', 'config_length', 'config_masspole', 'action_command']\n", + " - output_cols: ['state_x_position', 'state_x_velocity', 'state_angle_position', 'state_angle_velocity']\n" + ] + } + ], + "source": [ + "print(\"DATA STRUCTURE FROM UPLOADED MODEL:\")\n", + "print(\" - input_cols:\", model.features)\n", + "print(\" - output_cols:\", model.labels)" + ] + }, + { + "cell_type": "markdown", + "id": "imported-concentrate", + "metadata": {}, + "source": [ + "# 4. Feature Importance" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "everyday-attack", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAIiCAYAAADVboEmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xt8FNX9//HXR6Rys4iCVrSIN+QW5CbilfJFbbVVi1XRghJEKFZstdVCf71oK7aoiGhrUbQCiiKISmlpsYpU8YLc5FJA8a4UqwhKFUSBfH5/nElclk2yk2wyyeb9fDzyyGbmzMxnd3bz2TPnzDnm7oiIiEj+2CPpAERERCS3lNxFRETyjJK7iIhInlFyFxERyTNK7iIiInlGyV1ERCTPKLmLiNRxZtbAzNzM7oyxzX/NbE4ljvmt6JgXpCxrGy0bWdH9SqDkXsXM7BvRm7W0n57VEMNPzOziqj5OZZnZntFrMjPpWHLFzLqa2XVm1irpWKqamV0Ynb/DyyhzRDmfh3OrIc5LzWx4VR8nF6IE6mb2npntWUqZx1Nev6/l8Nh7Ru/d7+Rqn1J9Mr5ZpEpMBf6eYflr1XDsnwAvA/dVw7FkV12Ba4EngXcSjqWqfRdY6e6vZ1H2cWBKhuULcxtSRpcC+wB/rIZj5cI24GvA6cBfU1eY2deBU6IyDXJ83D0J7927gL9lWH8IUJTjY74CNAS253i/dY6Se/VZ6u6Z/pnVamb2FcDc/fOkY6lJzGxvd/8k6Tiqi5ntRUg+47Lc5OU8/TzUA77i7p/lcLergL2BQaQld6CQkAj/DpyTw2OWqyo+8x6GTN2W6/3Wtc8j6LJ8jRNd2nzOzD4xs61mtsDMdvvQmtn3zeyvZvaOmX1uZhvM7FEz65hSZk8zc+AgoE/a5c+DUy6D35Nh/5dG605MWTYqWtbWzMaZ2X+Az4BjUsp808yeMLPNZrbNzJab2dBKvibrzOxJM+tiZk+Z2adm9oGZ3WRm9cysoZndambrzewzM/uXmR1VyvPpbWbXR69bcXznlXLc75nZC2a2JTrm/EyXKFPi61b83IGlZjYKuDsqNj/ltb8n2q6pmd1gZgvN7MPoPL5qZr8zs4Zpxzgl2naAmQ02s9VR+bfM7KelxN/NzGaY2ftR2XfM7EEzOzStXC7O2SmEBPRYzO3KZGanm9nclNiWmdngDOXOMLOHzezNqNxHZjbHzI5PK/df4FjgKMvQPGaltCNb5vbhYdGyk83sN2b2JvA5cHZKmZ5mNsvMNkbn4GUzG2HhS0Ack4DvmFmLlH0bIbnPBD7KEPNDZrZborQs2tfNrC3hsw3wg5TXaVtKmd1eq+JlZtbDzJ6OPjsfmtmfzWy/8p6kldHmbmb9o8/jpxb+Nz5vZmenlSl5btE5e97MtgAPR+tbmNkfzOyN6H3yoZktNrMflxdbbaOae/VpZGbN05Z9nvpt0sxGAyMI38J/Rbjk9T3gETO7zN1TP4zDgfcJl8zeB44AhgLPm1mX6NLoTuAi4HZgPTA6ZftNlXguDwFbgDGAR8fHzC4D7gCeB64HtgLfBO4ys0Pd/eeVOGYr4J/Ag8B04FvANYRaS1fCe/n3wP7AT4HHzKyju6dfNhxDuOz3R8KX20uA6WZ2UWpN0syuILxua4DfRGUHAX81s8Hufm/aflsDc4FphH8kjYGnCJdTB0evx9qobHFTzNej4z8CPADsAHoDI4GjgW9neB2GR8/xz8Bm4GJgjJm96+7TU+I/O4rjE+Ae4PUolm8B7YE3o3K5Omd9gbfd/aUsyzfI8HnY7u6bU55D8TmYH8X2WRT/PWbW2t1/lbLtYKAJIQn+h/DaXgr8y8xOdPfiy/2XAzcCjYCfpWxfmeax2wnvjzuBTwmvNWbWl/BeXQ3cDHwMnEh4n3YkfDazNZnwGgwAbo2W9QIOAy4Dzq9E/JmsJ7w37yW8jydGy3dmse2hwBOE/xMPAT2ifXU1s54VqfGb2RjC5/pvwC+ixecCM81siLunV1BOAL4PTIhiL457JtCdcK5WEj6n7YFvALfFjatGc3f9VOEP4U3jpfw8lFKuR7Tstxn28TfCP4bGKcsaZyjXEfgCuD1t+TrgyQzl94yOeU+GdZdG605MWTYqWjYXqJdW/mBCreW+DPu6g5C4DinntSqOZ2aG+B3om7Z8OeEL0COEpoHi5T+JyvfJ8HzeAL6asnyfaP8fAg2iZfsRktwrwN4pZZsCbxGS6lczxFeYzeuYsu4rwJ4Zlv8+2qZryrJTomXvph27CbARmJ9h2X+BAzPsf49cnbPi/QEfAOOyKHtEGZ+HZ1PKHRK9l+/NsI+7CF/qDi7n83AQ4XPzaNryBYRmgUzx/ReYk2H5t6IYL0hZNixatrL4vZPhHDzB7p+Vn0fb9czi9fovsDh6PBtYkbLuPkI/jj0IX+Ac+FrK+oeAbRn22SAqe2fcZeW9VtEyB4aV8pyvLOc1bRstG5my7Pho2a/T9mnAHEJFpWFazA6clFZ+/2j52PJe93z40WX56jMBODXtZ1TK+v6EN959ZtY89QeYRUgsxxYXdvctEC7NmdlXo3L/JdRAjqVq3eru6d/gzyMkq3szxP9XoB7QpxLHfMfd0y/5Pkv4gP/Bo09vZH70+8gM+/mTu/+v+A93/5iQLPYDTooWf5NQu7/NU66seKhV/hH4KvB/afv9gJgdFt39C3ffAWBm9c2sWfR6PREVyXQe702L/1NCJ7TU53o6sC9ws7u/l+G4xVczcnXOTgRaEGpF2XqU3T8PV6asPx+oX0Zse5JyDoo/DwBm1iS6BPw5sJiq/zzc4e7pl7+Lz8FEoFla/LOjMqfFPM69QIGF5pa9CVf1JvnuV6eStpHwZSPVbYQrL30rsL/+hC/xU9Jex/2AvwDNSGkajLzo7vPTlm0hfGE93urA3Su6LF99XnX3J8tY346QqF4to8wBxQ/MrBvhMt3JhEtLuxyrokFmaW2GZe2i3/PK2O6AMtaV540My4rbGd8sZXmmNr41GZatjn4fFv0ubpNelaHsv9PKFnu9Iv9kLdyS9QPCpcH0L9vNMmyS6XXYyK7PtTjRl3eJPFfn7LtRDOn/TMvybhafB8rZZ+rnoQ1wA+FLQtO0crns3JZJWZ+HB8rYLu7n4a+EK0yDCFetGhKaIWqaV4u/tBZz961m9ja7f26y0Y7w2SjrLoz013K3c+LuW8zsakITyVtmtorQ5PCouz9dgbhqNCX3msMI7UJnUPrtJf8GMLPWwDOEy1G/JbyRtxBq/n8g1Hiy4WWsK+u9sTXDMot+9yfUYjOpTLtmWW19pa2zDMsyPef0cpm2K0+m16RMZvYzQvvvHEIv8/cIl6JbEdrUM11Zy+a5Fj8u6/ymlqvsOfsu8NcMV3Mqozi2CwhfHDJ5FcDM9iF8CagPjCV8Tj4lfI6uZfdaXVlKe80q+nn4MV9+eUy3LtugIFzpMbMHCP0segJPu3umL3slm5SyvKr/75d23Ip8roq3207og1Lavlem/Z3x8+jut5nZI9G+Tia8v35kZve5+8AKxlcjKbnXHK8S2lXfdPfyat7fI3QI+lbqpaeo92xzQptwqowfCHffaaFn974ZVsf9hl0c84ZyamRJa8/u4w0U17KK/1EW1xA6AOnf6NunlS1PWQn2ouhYZ6Q2K1jlBw15JfrdhbJr5ZU+Z2bWmXCl48ryysZUHNsH7l7Wc4DQjLI/8H13n5oW35gM5cs6J5vI7efhkxx/Hu4lfGHoRkjyZdkE7GVmjdw9Ndll+1zK+3JYmjZmtmdq7d3MGhG+tC6qwP5eJfRdes3d06/Sxebu6whNcXdZGBhoGnCxmY1x9/QvCbWW2txrjvuj37/PdJuMmaVediquIaV/Ex5GSO7pPiXzPywIH5wTzKxkAIyovTLut9hphFrnb1P3lbLPfSzcE5+0H0btlUBJre8H7HpZ+Z+ES7k/MrPGKWW/Suit/j9Cp8JsfBr9zvT67yT8Ay05j2ZWn9BbvjKKOxldYxlGLIu+BEJuzllfwlWjJ8opF9dDhNra9RbuoU+PrVn0WkEpnwczO5Nw10G6sj4Pawnt2vun7KchoUd6HH8jNA/9wszSmwkws0Zm1iTmPnH3FYRe49cROpKWpfjS9Clpy3+S5bE+J5yD0l6r0uxH6Eia6seEZoSKjD5Z/L9xtJntlrPS/jeWyswaW9otptEXkOKEHvd51miqudcQ7v6CmV1PuAVuqZnNIFymPZBw68ZphA8HhA45vwMeMLM7CDX1Ewi9TzN9s10ADDSz3xBqdUXAXzwMtPFHQrvdU9Elv2aEW+reJNSGso3/7aj9+E5gtZlNIfTkbQF0Itz724aYlyKrwEfAi2Y2iZAMBhN6VRcWd4py940W7rO9LSo7mfBFuJBQSx3s2Q+IsZCQwH9p4R7lLYT2+UXADEK/ib9bGHK3KeESeaUGB3H3T83sUsJtWCvN7M+EKwT7E94jo4HZOTpnfYHHPbeDtuDub5jZj4A/pcT2bvQcimM7lNCJ9F+Etug/RG3v6wk12wsJ/SaOSNv9AuAUMxtHqEnuBJ5w942Ez8N3CZ+HCYTP3EB2vxpWXvz/M7OBhHO81swmEs5BM8KVonMIVxwWxNlvtO+xWRa9n9BsN8nMxhLuHPgO4Q6RbL0InG5m1xBe/53u/nA526wFboyu6iwn3AlUSEii42McGwB3n29mvyf0uD8quqz+HtCS8L+xN2GMhfIUAHPM7FHC++Jjwh1Gw6KYY5+LGi3p7vr5/sOXt8JdnWX5swi1oI8I/+TfIVxGHpJhv88R7mP+iFBTaE/oQf5aWtmvEXonf0RI7E50GxEhwY2MjvMFoX1wIGXfCndwGfGfSPh2viHa33pCp5WrgL3Kee5l3QqX6Va+jPHw5e1Wv0xZVvx8ehMS6rvR67sC6FdKPOcCLxDa77ZEr+2ZGcpljC9l/SWE4X+/IOXWQ0Jv9F8S/ulvI9xmN5rwDyc9/uJb4QZk2P8UYEeG5T0JvYk3pryXpgCtc3HOCJd3M8ZUxjbF56bc2+ai8icT7hYpju0/hKsmVxJGgisu15XwufmY8Jl4CjiODLeDERLB5GifxZ+HninrLyVc0foiOjdXEXq/l3YrXKm3tBGuHEwlJKPthC8jzwH/D2iaxfMvuRWunHK73QqXcm4XROf/A8Itji3I8rY3wheRudFr6qmvJaXfCjeHcIfC04TPzUbCXQPN08pmdStcyrqzCcM4p/5vnE34sl3m84jW7U8Yk2BF9D7ZGp3nW4D9s30P15Yfi560SF6LarJ3E+59fTbpePKBhZHxRhP+Me42QprUPRZGAFzm7t9KOpa6Tm3uIlJR/wGuUGIXqXnU5i4iFeLuDyUdg4hkppq7iIhInlGbu4iISJ5RzV1ERCTP1No29+bNm3vr1q2TDkNERKTaLFmy5EN3b1FeuVqb3Fu3bs3ixYuTDkNERKTaRBPwlEuX5UVERPKMkruIiEieUXIXERHJM7W2zV1EpDbavn0769atY9u2bUmHIjVYgwYNOPjgg6lfv375hTNQchcRqUbr1q1j7733pnXr1nw5+67Il9ydjRs3sm7dOg499NAK7UOX5UVEqtG2bdvYb7/9lNilVGbGfvvtV6mrO0ruIiLVTIldylPZ94iSu4iISJ5Rm7uISIJaj5yd0/29NfrbFdpu3LhxDB06lEaNGuWkXE0wa9YsVq9ezciRI5k5cyZt2rShffv2APz617/m5JNP5pRTTkk4yqqhmruIiDBu3Di2bt2as3I1wVlnncXIkSMBmDlzJqtXry5Z99vf/jZvEzsouYuI1Dlbtmzh29/+NkcffTQdO3bkN7/5DevXr6d379707t0bgMsuu4zu3bvToUMHrr32WgBuv/323cr985//5LjjjqNr166cd955fPrppxmPuXnzZo466iheeeUVAC688ELuvvvuUmNs0qQJP/3pT+natSt9+vRhw4YNACxbtoyePXvSqVMn+vbty0cffVQSW/v27enUqRMXXHABAJMmTWL48OE8//zzzJo1i2uuuYbOnTvz+uuvU1hYyIwZMwCYO3cuXbp0oaCggEsuuYTPP/8cCMOcX3vttXTt2pWCggJefvnlSr3u1UnJXUSkjpkzZw4tW7Zk+fLl/Pvf/+bKK6+kZcuWzJs3j3nz5gFwww03sHjxYlasWMHTTz/NihUr+NGPfrRLuQ8//JBRo0bx5JNPsnTpUrp3787YsWMzHrNp06b88Y9/pLCwkIceeoiPPvqIIUOGlBrjli1b6Nq1K0uXLqVXr1785je/AeDiiy/mxhtvZMWKFRQUFJQsHz16NC+99BIrVqzgzjvv3GVfxx9/PGeddRY333wzy5Yt4/DDDy9Zt23bNgoLC5k2bRorV65kx44djB8/vmR98+bNWbp0KZdddhljxoyp2AueACV3EZE6pqCggCeffJIRI0Ywf/58mjZtuluZ6dOn07VrV7p06cKqVat2uaRdbMGCBaxevZoTTjiBzp07M3nyZN5+u/R5TU499VQKCgq4/PLLueeee8qMcY899qBfv34ADBgwgGeffZbNmzfz8ccf06tXLwAGDhzIM888A0CnTp3o378/U6ZMYc89s+9O9sorr3DooYfSpk2b3fYJcM455wDQrVs33nrrraz3mzR1qBMRqWPatGnDkiVL+Pvf/87Pf/5zTjvttF3Wv/nmm4wZM4ZFixbRrFkzCgsLM95z7e6ceuqpTJ06NavjFhUVsWbNGho2bMimTZs4+OCDs465vFvDZs+ezTPPPMOsWbO4/vrrWbVqVVb7dfcy1++1114A1KtXjx07dmQXbA2gmruISB2zfv16GjVqxIABA7j66qtZunQpe++9N5988gkA//vf/2jcuDFNmzbl/fff5x//+EfJtqnlevbsyXPPPcdrr70GwNatW1m7dm2px7311ltp164dU6dO5ZJLLmH79u2lli0qKippE3/wwQc58cQTadq0Kc2aNWP+/PkA3H///fTq1YuioiLeffddevfuzU033cTHH3+8W9t/atyp2rZty1tvvVXyHIr3Wdup5i4ikqCK3rpWGStXruSaa65hjz32oH79+owfP54XXniB008/nQMPPJB58+bRpUsXOnTowGGHHcYJJ5xQsu3QoUN3KTdp0iQuvPDCkk5oo0aNKrnEnWrt2rXcc889LFy4kL333puTTz6ZUaNGlbSZp2vcuDGrVq2iW7duNG3alGnTpgEwefJkhg0bxtatWznssMOYOHEiO3fuZMCAAWzevBl356qrrmKfffbZZX8XXHABQ4YM4fbbby/50gBhDPeJEydy3nnnsWPHDo455hiGDRtW6dc4aVbeJYkqO7DZnu6+I3rcBGgLvOHum7LZvnv37r548eKqDFFSXbd7m1zZ5TdXTRwitdyaNWto165d0mHUeE2aNCm1531dkem9YmZL3L17edsmclnezAqB981srZmdDqwAbgSWm9mFScQkIiKSL5K6LP9T4Chgb2A50MXdXzezA4AngIy9M8xsKDAUoFWrVtUUqoiIxNG3b1/efPPNXZbdeOONfPOb39yt7LHHHltySb/Y/fffX+dr7ZWVVHLf6e4fAh+a2afu/jqAu79fVo9Id58ATIBwWb5aIhURkVgee+yxrMu++OKLVRhJ3ZVUcn/HzH5PqLm/bGa3AI8CpwDvJRSTiIhIXkjqVrgBwP+AdcBZwPPAz4H9gcKEYhIREckLidTc3f1/wO9TFj0S/YiIiEgl1bhBbKJOcyIiIlJBNXEQm7LHGBQRySdxx5Aod38VG2Oits/n3rp1axYvXkzz5s1jbXfnnXfSqFEjLr74YiZNmsRpp51Gy5YtqyjK6pNYzd3M2prZCDO73cxuix63c/e7kopJRKSuysf53LMxbNgwLr74YiBMEbt+/fqEI8qNpAaxGQE8RKilLwQWRY+nmtnIJGISEakravp87uPHj+dnP/tZyd+TJk3iiiuuAGDKlCn06NGDzp0784Mf/ICdO3futv3YsWPp2LEjHTt2ZNy4cSXL77vvPjp16sTRRx/NRRddBMB1113HmDFjmDFjBosXL6Z///507tyZ2bNn07dv35Jtn3jiiZIZ4mqDpGrug4Fj3H20u0+JfkYDPaJ1IiJSRWr6fO7nnnsujz76aMnf06ZNo1+/fqxZs4Zp06bx3HPPsWzZMurVq8cDDzywy7ZLlixh4sSJvPjiiyxYsIC7776bl156iVWrVnHDDTfw1FNPsXz5cm677bbdjtm9e3ceeOABli1bxhlnnMGaNWvYsGEDABMnTmTQoEEVfs2rW1Jt7kVASyB94t8Do3UiIlJFCgoKuPrqqxkxYgTf+c53OOmkk3YrM336dCZMmMCOHTt47733WL16NZ06ddqlTOp87gBffPEFxx13XKnHPfXUU3n44Ye5/PLLWb58eanlWrRowWGHHcaCBQs48sgjeeWVVzjhhBO44447WLJkCccccwwAn332Gfvvv/8u2z777LP07duXxo0bA2E+9vnz52NmnHvuuSVt8vvuu2+Zr5GZcdFFFzFlyhQGDRrECy+8wH333VfmNjVJUsn9SmCumb0KvBstawUcAQxPKCYRkTqhNszn3q9fP6ZPn07btm3p27cvZoa7M3DgQH7/+9+Xul1pk6G5e7lzwqcbNGgQZ555Jg0aNOC8885jzz1rYh/0zBK5LO/uc4A2wG+Ax4F/AtcBR0XrRESkitSG+dzPOeccZs6cydSpU+nXrx8Affr0YcaMGXzwwQcAbNq0ibff3vUC8Mknn8zMmTPZunUrW7Zs4bHHHuOkk06iT58+TJ8+nY0bN5Zsmy59zveWLVvSsmVLRo0aRWFhYamx1kSJfQ1x9yJgQVLHFxGpERKYHrk2zOferFkz2rdvz+rVq+nRowcA7du3Z9SoUZx22mkUFRVRv3597rjjDg455JCS7bp27UphYWHJNpdeeildunQB4Be/+AW9evWiXr16dOnShUmTJu1yzMLCQoYNG0bDhg154YUXaNiwIf3792fDhg20b9++4i94AhKbz72yNJ97NdN87iI5ofnca5fhw4fTpUsXBg+u/r7elZnPvfY0IIiIiFSjbt260bhxY2655ZakQ4lNyV1ERHIqF/O5FxQUVGmM2ViyZEnSIVSYkruIiOSU5nNPXo2bOEZEREQqR8ldREQkzyi5i4iI5Bm1uYuIJKhgcm47jq0cuDKn+5PaSTV3ERGpdVO+tm7dmg8//DCRY8+aNYvRo0cDMHPmTFavXl2y7te//jVPPvlkInGlUnIXEZFal9yTdNZZZzFyZJidPD25//a3v+WUU05JKrQSSu4iInVMEvO5Q0h8xxxzDB07dmTo0KElk7x84xvfYMSIEfTo0YM2bdowf/58IIxVf/7559OpUyf69evHscceS6aRSbOZ471YkyZN+OlPf0rXrl3p06dPyZSuy5Yto2fPnnTq1Im+ffvy0UcflTzn9u3b06lTJy644AIgzC8/fPhwnn/+eWbNmsU111xD586def311yksLGTGjBkAzJ07ly5dulBQUMAll1xScj9/69atufbaa+natSsFBQW8/PLLWZ657Cm5i4jUMUnM5w5hKNdFixbx73//m88++4y//e1vJet27NjBwoULGTduXMl483/6059o1qwZK1as4Fe/+lXGQWWymeM91ZYtW+jatStLly6lV69eJce6+OKLufHGG1mxYgUFBQUly0ePHs1LL73EihUruPPOO3fZ1/HHH89ZZ53FzTffzLJlyzj88MNL1m3bto3CwkKmTZvGypUr2bFjB+PHjy9Z37x5c5YuXcpll13GmDFjSo23opTcRUTqmIKCAp588klGjBjB/Pnzadp097kjpk+fTteuXenSpQurVq3a5dJzsdT53Dt37szkyZN3m6Ut1bx58zj22GMpKCjgqaeeYtWqVSXrzjnnHCAM+frWW28BYW724tpyx44dd5tPHkLtuHiO986dOzN37lzeeOONUmPYY489SmaZGzBgAM8++yybN2/m448/plevXgAMHDiQZ555BoBOnTrRv39/pkyZEmvK11deeYVDDz20ZBKd1H2W9nxzSb3lRUTqmCTmc9+2bRs//OEPWbx4MV//+te57rrrdtnnXnvtBUC9evXYsWNHyf7Lk80c72Upb4732bNn88wzzzBr1iyuv/76Xb6QlBdXWTI931xSchcRSVASt66tX7+efffdlwEDBtCkSRMmTZpUMpd58+bNM87n/o1vfANgl3I9e/bk8ssv57XXXuOII45g69atrFu3LuOUr8WJvHnz5nz66afMmDGDc889t8w4TzzxRKZPn07v3r1ZvXo1K1fu/lr16dOHs88+m6uuuor999+fTZs28cknn+wyDWyqoqIiZsyYwQUXXMCDDz7IiSeeSNOmTWnWrBnz58/npJNO4v7776dXr14UFRXx7rvv0rt3b0488UQefPDB3foUpM8BX6xt27a89dZbJa9N8T6ri5K7iEgdk8R87vvssw9DhgyhoKCA1q1bc8wxx5Qb5w9/+EMGDhxIp06d6NKlC506ddqtCSGbOd5TNW7cmFWrVtGtWzeaNm3KtGnTAJg8eTLDhg1j69atHHbYYUycOJGdO3cyYMAANm/ejLtz1VVXsc8+++yyvwsuuIAhQ4Zw++23l3SkA2jQoAETJ07kvPPOY8eOHRxzzDEMGzas3OecK5rPXbKj+dxFckLzuWdv586dbN++nQYNGvD666/Tp08f1q5dy1e+8pUK77NJkyZl9uivSTSfu4iI5J2tW7fSu3dvtm/fjrszfvz4SiX2ukTJXUREcirOfO5l2XvvvTPe156N0uaJry219spSchcRkZyKM597Vanr88TrPncREZE8o+QuIiKSZ5TcRURE8oza3EVEErSmbW5vi2v38pqc7k9qJ9XcRUSk1k35Wl3zuVf0OHfeeSf33XcfEGaRW79+fa5DK5OSu4iI1LrkXtMNGzaMiy++GFByFxGRalBX53MfP348P/vZz0r+njRpEldccUXW+xg7diwdO3akY8eOjBs3rmT5fffdR6dOnTj66KO56KKLALjuuusYM2YMM2bMYPHixfTv35/OnTsze/Zs+vbJIZhMAAAgAElEQVTtW7LtE088UTJDXC4puYuI1DF1dT73c889l0cffbTk72nTptGvX7+s9rFkyRImTpzIiy++yIIFC7j77rt56aWXWLVqFTfccANPPfUUy5cv57bbbtvtmN27d+eBBx5g2bJlnHHGGaxZs4YNGzYAMHHiRAYNGlTW6aoQdagTEaljCgoKuPrqqxkxYgTf+c53OOmkk3YrM336dCZMmMCOHTt47733WL169W7zqafO5w7wxRdfcNxxx5V63Hnz5nHTTTexdetWNm3aRIcOHTjzzDOB0udz//GPfwxkN587wGeffcb++++f8fgtWrTgsMMOY8GCBRx55JG88sornHDCCdxxxx3l7uPZZ5+lb9++NG7cuCTe+fPnY2ace+65NG/eHIB999231OcPYYrZiy66iClTpjBo0CBeeOGFkrb5XEokuZvZPu7+cRLHFhGp6+ryfO79+vVj+vTptG3blr59+2JmWe2jtFjcvdw54dMNGjSIM888kwYNGnDeeeex5565T8VJXZb/0MyeNLPBZrZP+cVFRPJTu5fX5PQnG+vXr6dRo0YMGDCAq6++mqVLl+4yL3mm+dyLpZbr2bMnzz33HK+99hoQ2sjXrl2b8ZiZ5nMvT/F87kCZ87nPmDGDDz74AIBNmzbx9ttvl7rPc845h5kzZzJ16lT69euX9T5OPvlkZs6cydatW9myZQuPPfYYJ510En369GH69Ols3LixZNt06XO+t2zZkpYtWzJq1CgKCwvLfR0qIqnL8muAccCFwE1m9iwwFfiLu39W2kZmNhQYCtCqVavqiFNEJO/U5fncmzVrRvv27Vm9ejU9evTIeh9du3alsLCwZJtLL72ULl26APCLX/yCXr16Ua9ePbp06cKkSZN2OWZhYSHDhg2jYcOGvPDCCzRs2JD+/fuzYcMG2rdvX+7rUBGJzOduZkvdvWv0uCFwJnAB0At43N2/X94+NJ97NdN87iI5ofncs1cV87nXFMOHD6dLly4MHjy41DK1cT73kgaKqKY+HZhuZk2B7yYUk4iI1CD5Op97t27daNy4MbfcckuVHSOp5J7xPgV33wxMruZYREQkh2ryfO4FBQUV2l8uZbqlL9cSSe7uPiaJ44qI1AQV6WFdm2g+98qrbJN5jRvEJuo0JyKSlxo0aMDGjRsr/c9b8pe7s3HjRho0aFDhfdTEQWzy9+usiNR5Bx98MOvWrSsZoUwkkwYNGnDwwQdXePvEkruZtQXOBg4CHFgPzHL3u5KKSUSkqtWvX59DDz006TAkzyVyWd7MRgAPEWrpC4FF0eOpZjYyiZhERETyRVI198FAB3ffnrrQzMYCq4DRiUQlIiKSB5LqUFcEtMyw/MBonYiIiFRQUjX3K4G5ZvYq8G60rBVwBDA8oZhERETyQlL3uc8xszZAD0KHOgPWAYvcfWcSMYmIiOSLxHrLu3sRsCCp44uIiOSrGjeIjYiIiFSOkruIiEieUXIXERHJM0ruIiIieUbJXUREJM8ouYuIiOQZJXcREZE8o+QuIiKSZ5TcRURE8oySu4iISJ5RchcREckzSu4iIiJ5RsldREQkzyi5i4iI5BkldxERkTyj5C4iIpJnlNxFRETyjJK7iIhInlFyFxERyTNK7iIiInlGyV1ERCTPKLmLiIjkGSV3ERGRPKPkLiIikmf2jLuBmX0d6AMcAEx193fMrD7QAtjg7ttzHKOIiIjEEKvmbmY3AK8D9wK/A46IVjUC1gKX5TQ6ERERiS3r5G5mQ4CfAxOAMwArXufum4G/AmflOkARERGJJ07N/XLgL+4+HFiUYf1yoG1OohIREZEKi5PcjwIeL2P9BqB55cIRERGRyoqT3D8ntK2XphWwuSJBmNkRZvY9M2tfke1FRETkS3GS+0Lgu5lWmNlewADguWx2ZGbzzKx59Pgi4O/A6cA0M7siRkwiIiKSJk5yvwU4wcwmAh2iZS3MrA/wFKHmPibLfbVw9w+jxz8CjnP3S4FjgSGlbWRmQ81ssZkt3rBhQ4zQRURE6o6sk7u7Pw4MBy4E5kWLHwT+CXQDLnP357Pc3XYzOyh6/CmwJXr8OVCvjBgmuHt3d+/eokWLbEMXERGpU2INYuPu481sFnA+oWe8Aa8C09z9nRi7ugr4p5k9AqwCnjKzOcBJwMQ4MYmIiMiuYo9Q5+7/AW6tzEHd/V9mdjzwfWBvYAmh1n6Fu79cmX2LiIjUdVkndzM7BGjv7v8oZf3pwKpsa/DRwDfjsz2+iIiIZCdOzf13wCFAxuROGL3uTWBgZQIys6HuPqEy+xCRmmNN23axyrd7eU0VRSJSd8TpLX8iMKeM9Y8DvSoXDpAyrK2IiIjEF6fmfgDw3zLWvx+VyYqZtQXOBg4CHFgPzHL3u2LEJCLVrGByQazy06soDhEpXZya+8fAYWWsP5xwW1u5zGwE8BChlr6QMFa9AVPNbGSMmERERCRNnJr7s8AQMxvn7h+krjCzA4DBUZlsDAY6pM/9bmZjCbfGjY4Rl4iIiKSIU3P/PfBV4CUzu9LMvmFmvczsSmAp0DQqk40ioGWG5QdG60RERKSCsq65u/sSMzsfuBcYS2gnh3A5fRPQz90XZrm7K4G5ZvYq8G60rBVwBGEUPBEREamguCPU/cXMWhEmeTmSkNhfAea4+5YyN951P3PMrA3Qg9ChzoB1wCJ33xknJhEREdlVRUao2wLMqOyB3b0IWFDZ/YiIiMiu4rS5i4iISC0QK7mb2blm9rSZrTezz83si7Sfz6sqUBEREclOnLHlfwLcDHxEuDd9Y1UFJSIiIhUXp839CsJgM33idJ4TERGR6hXnsvyBwH1K7CIiIjVbnOT+OmGgGhEREanB4iT3scBgM2tcVcGIiIhI5cVpc/8M2ACsNrM/E+Zu323AGXd/MEexiYiISAXESe5TUh5fV0oZB5TcRUREEhQnuZ9aZVGIiIhIzsSZOGZuVQYiIiIiuaHhZ0VERPJM7IljzKwzcCzQjN2/HLi7Zzunu4iIiFSBOMPPNgAeBs4gTNHq0W9SHjug5C4iIpKgOJflf0VI7DcCpxCS+SXAWcALhKFpC3IdoIiIiMQTJ7mfBzzi7v8PWB4te8fd/wb0BhoA/XMcn4iIiMQUJ7m3AuZFj4sHr/kKgLtvJ9zffmHuQhMREZGKiJPcP+HLNvpPgSLCZDLFPga+lqO4REREpILiJPc3gCMB3H0HsBr4Xsr67wL/yV1oIiIiUhFxkvuTwPfMrHibCcAZZrbWzF4BTgMm5jpAERERiSfOfe43EtrV9wCK3P2PZtYIGEBog/81MDr3IYqIiEgccYaf/R+wKm3ZTcBNuQ5KREREKi7ry/JmNsHMepSxvruZTchNWCIiIlJRcdrcLwWOKGP94cDgyoUjIiIilZXLiWMaA9tzuD8RERGpgDLb3M3sYMLgNcWONLPjMxTdFxhKuF1OREREElReh7rBwLWECWGc0CP+1xnKFU8ao8vyIiIiCSsvuc8C1hGS9wTgz8CCtDJOGLFuobu/lesARUREJJ4yk7u7vwS8BGBmhwDT3X1ldQQmIiIiFZPVfe5m1gQ4H/gQUHIXERGpwbLqLe/unwIHAFuqKhAz27eq9i0iIlKXxLkVbiHQLRcHNbNfpjxub2ZrgSVm9paZHVvGdkPNbLGZLd6wYUMuQhEREck7cZL7SKCfmV2Ug+Oek/L4ZuDH7n4o4dL/raVt5O4T3L27u3dv0aJFDsIQERHJP3EmjhkNbAQmmdlNwOvA1rQy7u7fjBlDS3f/R7TxQjNrGHN7ERERSREnubcn3Pa2Pvr7kAxlPMt9HWZmswi32B1sZo3cvfiLQv0YMYmIiEiaOLPCHZzD456d9vceAGZ2ADA+h8cRERGpc+LU3HPG3Z8uZfn7wB3VHI6IiEheiZ3co3ve/w84LFr0BjDX3XNym5yZDXV3TR0rIiJSQbGSu5kVEnqzf5XQXg6hnf1/Znalu0/OQUxWfhEREREpTdbJ3cy+DdwLvA2MAv4dreoADAfuNbMN7v73LPfXltD2fhBfdtSb5e53ZR++iIiIpItzn/sI4BWgk7vf4u6PRz9jgc7RupHZ7MjMRgAPEWrpC4FF0eOpZpbVPkRERCSzOJfluwDXu/sn6Svc/X9mNhn45e6bZTQY6ODu21MXmtlYYBXhnnoRERGpgDg1dwOKyli/k+zby4uAlhmWH1jOMURERKQccWruK4BCM/tTyoAzAJhZI2AgsDzLfV0JzDWzV4F3o2WtgCMI7fciIiJSQXGS+y3Aw8BiM7sNWB0t7wD8CDgKOC+bHbn7HDNrA/QgdKgzYB2wyN13xohJRERE0sQZoe4RM/sxoT18PF8ONWvAZ8CV7v5ojP0VAQtixCoiIiJZiHWfu7v/wcweAL4JHEpI7K8D/3T3TVUQn4iIiMQUe4S6KIlPrYJYREREJAcqMvzsV4CT2HX42fnu/kUuAxMREZGKiTv87PeBccB+7Dr87EYz+4m7T8lxfCIiIhJTnOFnzwWmAP8BriX0ljfCPO8/ACab2TZ3n1EVgYqIiEh24tTcfwmsBY51980pyx8xsz8QhpH9JaDkLiIikqA4I9QdBdybltgBcPePgT9HZURERCRBcZL7B5Q9NGxRVEZEREQSFCe530cYfrZx+gozawIUApNyE5aIiIhUVJw297nA6cByM7sDeJnQU7498EPgI+ApMzs+dSN3fz5HsYqIiEgW4iT3p1Ie38Kuw89mKmNRmXoVC01EREQqIk5yH1JlUYiIiEjOxJk45s9VGYiIiIjkRpwOdSIiIlILVGRs+d7Akew6BG0xd/ff5yIwERERqZg4w8+2AR4F2rF7Ui/mgJK7iIhIguLU3O8izOF+NTCfcOubiIiI1DBxkvuxwM3ufmtVBSMiIiKVF6dD3UdoeFkREZEaL05ynw6cVVWBiIiISG7ESe4/B3aa2XQzO8nMvm5mLdN/qipQERERyU6cNvfPgaWEJP+9MsppuFkREZEExUnuowk95VcCz6Le8iIiIjVSnOReCPzF3c+polhEREQkB+K0uTcC5lRVICIiIpIbcZL7QsLodCIiIlKDxUnuPwW+b2a6HU5ERKQGi9PmfhPwMfCYmb0DvAHsTCvj7v7NXAUnIiIi8cVJ7u0JE8Osj7Zrk6GM5yIoERERqbisk7u7H1yVgZjZWe4+qyqPISIiUhfEns89F8ws/XY6A+4wsz0B3P3RUrYbCgwFaNWqVZXGKCIiUlslktwJ49TPIUxEUzw3fGPgTMKl/YzJ3d0nABMAunfvriYAERGRDMpM7maWMcmWwd29rKFpix1HGPFuEXCnu7uZfcPdB8U8noiIiKQpr+b+3Zj7y6o27e6LzOxU4ArgKTMbke22IiIiUrbyknv9qjqwuxcBt5nZw8C4qjqOiIhIXVNmcnf39PvYc87d1wPnV/VxRERE6oo4I9RVi6hHvIiIiFRQjUvufNl7XkRERCogqVvhMLO2wNnAQXw58t0sd78rqZhERETyQSI196h3/EOEWvpCwi1xBkw1s5FJxCQiIpIvkqq5DwY6uPv21IVmNhZYRbgHXkRERCogqTb3IqBlhuUHRutERESkgmLX3M3s60Af4ABgqru/Y2b1gRbAhvTaeCmuBOaa2avAu9GyVsARwPC4MYmIiMiXYiV3M7sBuCbazglt5e8AjYC1wP8Dbi9vP+4+x8zaAD0IHeoMWAcsqo5760VERPJZ1sndzIYAPwf+BPwN+HvxOnffbGZ/Bc4ii+QebVMELIgVrYiIiJQrTpv75cBf3H04ocaebjnQNidRiYiISIXFSe5HAY+XsX4D0Lxy4YiIiEhlxUnunxPa1kvTCthcuXBERESksuIk94WUMgWsme0FDACey0VQIiIiUnFxkvstwAlmNhHoEC1rYWZ9gKcINfcxOY5PREREYsq6t7y7P25mw4FbgYujxQ9Gv7cDl7n78zmOT0RERGKKdZ+7u483s1mE+dfbEu5PfxWY5u7vVEF8IiIiElNWyd3MvgJ0B95399cJtXcRERGpgbJtc3fgaeDbVRiLiIiI5EBWyT0aL/59wmV4ERERqcHi9JafAZxrZkrwIiIiNVic5D4eaAI8bmanm9kRZtYy/aeK4hQREZEsxektv4bQ9m6EKV9LU69SEYmIiEilxEnuvyMkdxEREanB4gxi88uqDERERERyI06bu4iIiNQCWdfczez4bMppCFoREZFkxWlzf5bs2tzVoU5ERCRBcZL7kFK2P5wwkcwbwJ9zEZSIiIhUXJwOdaUmbjO7EVgKfJ6LoERERKTictKhzt03AncDI3OxPxEREam4XPaW3wQckcP9iYiISAXkJLlHU8L2B/6bi/2JiIhIxcW5FW5CKav2BU4ADkCX5UVERBIXp7f8paUs3wy8Coxw9/sqH5KIiIhURpzkXj/DMnf3olwFI9Wn9cjZscq/1aCKAhERkZyL0+Z+AFDf3Xem/JQkdjPbS1O+ioiIJC9Ocn8XOKeM9d+NyoiIiEiC4iR3y2JfmhJWREQkYXFvhSsreR9F6FxXIWa2b0W3FRERkS+V2aHOzC4CLkpZ9HMzG5Sh6L5AZ2BmNgc1sxOAe4Ai4BJgFHC4mdUHznf3F0rZbigwFKBVq1bZHEpERKTOKa+3fHOgXfTYgYOAZmllHPgUuA/4eZbHvRU4H2gCzAa+6+7PmllX4A+E++Z34+4TgAkA3bt3VxOAiIhIBmUmd3e/lZCIMbMi4Ap3fzAHx63v7iuj/W5w92ej4y01s4Y52L+IiEidFes+d3ffmaPjprb1p9f2v5KjY4iIiNRJWXeoy2FiB/iVmTWK9lvSTm9mhxMu74uIiEgFxam5Y2aHAj8GjiW0vad/OXB3P6q8/bj7rFKWvw7cFCcmERER2VXWNXcz6wC8BFwGfBVoA+wA9iFM9VoP+KCyAUU94kVERKSC4tzn/ltCMu8C9IqWDXf3/YHLgb2BITmIqbzBckRERKQMcS7LnwRMcPfVZrZftMwA3H28mR0P3Aicnc3OzKxtVPYgwu1064FZ7n5XjJhEREQkTZya+1eB16LHX0S/G6esfw44MZsdmdkI4CHCl4OFwKLo8VQz05zwIiIilRCn5v4+YWY43P0TM9sCHJmyvmmM/Q0GOrj79tSFZjYWWAWMjhGXiIiIpIiT3JcD3VP+ng/8yMxeIFwBuBxYkeW+ioCWwNtpyw+M1omIiEgFxUnuDwGXm1lDd/8M+BXwNCHJA2wDBmS5ryuBuWb2Kl9OE9uK0Ot+eIyYREREJE3WyT0advbBlL+XRLfHfQ/YCcx299dK2z5tX3PMrA3Qg9ChzoB1wKIcD5YjIiJS58QaxCadu78NjK3gtkXAgsocX0RERHYXO7lHE7scS+hcN8/dKz1wjYiIiOROnFvhikeP+w8wl3CJvmO0fH8z22Jml+Q+RBEREYkjzvCzfYE7CfezX0bKSHJR7f1J4JxcBygiIiLxxKm5/wx42t3PBB7JsH4RUU1eREREkhMnuXcic1Iv9h7RIDciIiKSnDjJfWc55Q8EtlYuHBEREamsOMl9BXBqphVmtgdwLrA4F0GJiIhIxcVJ7ncA3zazawnjyAOYmR1OGL2uAPhDjuMTERGRmOKMUDfVzI4GriUMPQvwD6Aeoef89e7+t9yHKCIiInHEGsTG3Uea2aNAf6AtIam/Ctzv7hptTkREpAaIPUKduy8kzMEuIiIiNVCZbe5m1sPM9q2uYERERKTyyutQ9wLwreI/zKyJmT1oZu2rNiwRERGpqPKSu6X9vRdwAfC1qglHREREKivWxDEiIiJS8ym5i4iI5BkldxERkTyTza1wZ5hZcRt7I8CB88ysc4ay7u635iw6ERERiS2b5P796CfVD0op64CSu4iISILKS+69qyUKERERyZkyk7u7P11dgYiIiEhuxB5+Nl+1Hjk7Vvm3Rn+7iiIRERGpHPWWFxERyTNK7iIiInlGyV1ERCTPKLmLiIjkGSV3ERGRPKPkLiIikmeU3EVERPKMkruIiEieSSS5m1mBmS0ws3fNbIKZNUtZt7CM7Yaa2WIzW7xhw4bqCVZERKSWSarmPh64DigA1gLPmtnh0br6pW3k7hPcvbu7d2/RokXVRykiIlILJTX8bBN3nxM9HmNmS4A5ZnYRYWY5ERERqaCkkruZWVN33wzg7vPM7HvAI8C+CcUkkjfizJWgeRJE8k9Sl+VvBNqlLnD3FUAf4NFEIhIREckTidTc3f3BUpa/Awyp5nBERETySo27Fc7MhiYdg4iISG1W45I7YEkHICIiUpsl1aEOM2sLnA0cROghvx6Y5e53JRWTiIhIPkhqEJsRwEOEWvpCYFH0eKqZjUwiJhERkXyRVM19MNDB3benLjSzscAqYHQiUYmIiOSBpJJ7EdASeDtt+YHROpGcinPfN+jebxGp3ZJK7lcCc83sVeDdaFkr4AhgeEIxiYiI5IWk7nOfY2ZtgB6EDnUGrAMWufvOJGISqU5r2rYrv1CKdi+vqaJIRCQfJdZb3t2LgAVJHV9ERCRf1cT73EVERKQSEqu5S34rmFwQq/zKgSurKBIRkbpHNXcREZE8o+QuIiKSZ5TcRURE8oySu4iISJ5RchcREckzSu4iIiJ5RsldREQkzyi5i4iI5BkNYiOSA3EH7ZleRXGIiIBq7iIiInlHyV1ERCTPKLmLiIjkGSV3ERGRPKPkLiIikmfUW15qhDVt28Uq3+7lNVUUiYhI7aeau4iISJ5RchcREckzSu4iIiJ5RsldREQkzyi5i4iI5BkldxERkTyj5C4iIpJnlNxFRETyjJK7iIhInlFyFxERyTNK7iIiInlGyV1ERCTP1LjkbmZNko5BRESkNqtxyR1YXdoKMxtqZovNbPGGDRuqMyYREZFaI5EpX83sJ6WtAkqtubv7BGACQPfu3b0KQhMRkRqk9cjZscq/NfrbVRRJ7ZJUzf13QDNg77SfJgnGJCIikhcSqbkDS4GZ7r4kfYWZXZpAPCIiInkjqeQ+CNhYyrru1RmIiIhIvkkkubv7K2Wse786YxEREck3Na5928yGJh2DiIhIbVbjkjuhx7yIiIhUUFJt7phZW+Bs4CDAgfXALHe/K6mYRERE8kFS97mPAC4EHgIWRosPBqaa2UPuPjqJuETqpOuaxit/aKuqiUNEciapmvtgoIO7b09daGZjgVWAkruIiEgFJdXmXgS0zLD8wGidiIiIVFBSNfcrgblm9irwbrSsFXAEMDyhmERERPJCUve5zzGzNkAPQoc6A9YBi9x9ZxIxiYiI5IvEesu7exGwIKnji4iI5KuaeJ+7iIiIVIKSu4iISJ5J7LK8lG1N23axyrd7eU0VRSIiIrWNau4iIiJ5RsldREQkzyi5i4iI5Bm1uYuISJ1UMLkgVvmVA1dWUSS5p+ReTeK+iaZXURwiIpL/dFleREQkzyi5i4iI5BldlhcREclCbRp/RDV3ERGRPKOae0Vd1zRe+UNbVU0cUjV0fkWkFlPNXUREJM8ouYuIiOQZXZYXEZH8EadJLY+b01RzFxERyTNK7iIiInlGyV1ERCTPKLmLiIjkGSV3ERGRPKPkLiIikmeU3EVERPKMkruIiEieUXIXERHJM0ruIiIieUbJXUREJM8ouYuIiOQZJXcREZE8k+iscGZ2AHAQ4MB6d38/yXhERETygbl79R/UrDNwJ9AU+E+0+GDgY+CH7r60lO2GAkOjP48CXqniUKtDc+DDpIOoRnq++asuPVfQ881nNfm5HuLuLcorlFRyXwb8wN1fTFveE7jL3Y+u9qASYmaL3b170nFUFz3f/FWXnivo+eazfHiuSbW5N05P7ADuvgBonEA8IiIieSOpNvd/mNls4D7g3WjZ14GLgTkJxSQiIpIXEknu7v4jMzsdOJvQoc6AdcAd7v73JGJK0ISkA6hmer75qy49V9DzzWe1/rkm0uYuIiIiVafG3ece9YgXERGRCqpxyZ1wiV5EREQqKLHL8mbWli/b3B1YD8xy9zWJBCQiIpInEqm5m9kI4CFCLX0hsCh6PNXMRiYRk4iISL5IahCbtUAHd9+etvwrwCp3P7Lag6pGZtYGuAY4hJQ7Ftz9/xILqgqZ2SPAvcA/3L0o6XiqWl06v3Xw3J4AXMeX59YAd/fDkoyrqphZI+CnQCt3H2JmRwJHufvfEg4tZ8zsJ2Wtd/ex1RVLLiV1n3sR0BJ4O235gdG6fPcwYfjdu4GdCcdSHcYDg4DbzexhYJK7v5xwTFWpLp3funZu/wxcBSwh/88twETCcz0u+nsd4f2dN8kd2Dv6fRRwDDAr+vtM4JlEIsqBpGru3wL+CLzKl4PYtAKOAIa7e14PZGNmS9y9W9JxVDczawpcCPyCcN7vBqakX8Gp7eri+a1D5/ZFdz826TiqS/EwrGb2krt3iZYtz8chws3sn8D33P2T6O+9gYfd/VvJRlYxSQ1iMye6dNmDXQexWeTudeHb8F/N7IfAY8DnxQvdfVNyIVUtM9sPGABcBLwEPACcCAwEvpFcZFWiTp3funBuzaxr9HCemd0MPMqu5zbjZFd54Asza0jo9IyZHU7K884zrYAvUv7+AmidTCiVp0FsEmBmb2ZYnM/tdo8CbYH7CZdt30tZV+snaEhXl85vXTm3ZjavjNWej/0pAMzsVOCXQHvgn8AJQKG7/yvJuKqCmf0COJ/wpdyBvsB0d/9dooFVkJK7VDkzOyN9WGEz28vd87UGUGfUtXNrZoe5+xvlLcsn0ZWZnoQrrAvcvaZOhVppZtaNcNUJ4Bl3fynJeCpDyT0BZlYfuAw4OVr0L8JUt3nVPlnMzJa6e9f/396Zh8tVVen7/QhgSDCggFMzyGCDgkwKiNCKDIr4axshiAioKNCICjY0zhrGdhYVmVFIxAE10ig4IIPBMIgBQhjEpkWkFQERSBjCEPh+f+xTSaWoOySpU/vW2ZYKldoAAB5LSURBVOt9nnrqnqGqvlX73Fpn7732WiPtawoltW+BbdvN3sbFWLRNQ3SlqdMQksYBL2TxVS535VO09OSKli+dU4EVgFOq7f2rfQdmU1QDkl5EiqlYSdIWLMo+OAmYkE1Y/TS+fUtr2yrp1sbAKpL2aDs0CRifR1WtfHmYYwYaNw0h6UPAFOBe0koIkWzdNKeupSWcex626og2vUzSjdnU1MebgPcAawLta0UfBj6RQ1CfKKF9S2vbDYH/B6xKWiLV4mHgoCyKasT2G3JryMDhpDX8/8gtpBeEc8/D05LWt/1HSHN2NHDNrO2pwFRJe9qenltPH2l8+5bWtrYvAC6QtK3tq3Pr6RclTTGRlnDOzS2iV8ScewYk7URKDnEHaehnHeAA28NF5A4ckvazfa6kI6mW0rQzqJmfRqKE9i24bU/i2fbOBWZVNwCNQtJZpCmmqdWu/YGnbTdmiqmFpG+SRmguYvFljgN5LUfPPQO2L22lcST9+N/W0OjiidXzyllV9JlC2rfItgWeQ1r698Nqe0/gFuB9kt5g+8PZlNVDCVNMLe6qHitWj4Emeu59RNKOti/rCMhZiO0f91tT0DuifZuPpMuAN9peUG0vT1r/vQtwk+1X5NTXayRdD+zVMcX0o6auhmgS0XPvL68HLmPxgJwWJmW9ahySpgKH236o2n4e8GXb782rrOcU174FtW2LfyKNWrTmZicCL7H9tKSmjc5AKoB0uaTFppjySqoHSWsAHyGtili4AmJQExSFc+8jtqdUfx5re7EsZpLWzSCpX2za+vEHsP1gtXyqURTavkW0bRtfAGZL+jXJ2b0O+C9JE4FLcgqrg0KmmFp8BziPtCriEFL65L9nVbQMxLB8BkpJhNGimqPbwfaD1fbzgRm2X5lXWT2U1L6ltS2ApBeT6mIIuNb23Zkl1Yak8cChpKxtBn4DnGb78azCaqD1Pyppju1Nq30zbL8+t7alIXrufaTARBgtvgxcJelH1fZewAkZ9dRCoe1bRNt2sBypR7c8sIGkDWwPbGnQEZhGWst/UrW9D6mOwF7ZFNVHa3nf3yS9BbiblMdhIAnn3l+KSoTRwvY0SbNYlNVqD9u35tRUE8W1b0FtC4CkzwN7kyLkn6l2mwGu+z0CG3ZEy1/e4Gj546vSxUeSbmYmAf+RV9LSE869j5SaCKNiBRalc1whs5ZaKLh9G9+2bexOcnhNnXfu5AZJr7F9DYCkbYArM2uqBdsXVn/OBQY+Q1/MufcRSR+x/YUhEmFg+7AMsmpH0uGknut0khN4G3CG7ZOGfeGAUWL7ltK2LST9nLQ07JHcWvqBpN+TRqRaxVPWBn5PGrVwa266CVRBrx8i1XBvLxzz1lyaloXoufeX31fPs7Kq6D/vA7ax/SgsHNq8mkXzeE2hxPYtpW1bPEaKlr+UxbOYNe7GrWLX3AL6yH8D3wR+yqIpl4ElnHsfsf3T6rmVyhFJywEr256XTVj9iMVzq7cqLjWKQtu3iLZt4yfVoxSWB/5i+wlJO5AqpE1rX/7YIB63/fXcInpFOPcMSPouaR3l08B1pOjqr9j+Yl5ltXE28FtJ51fbu5PukBtJYe1bVNvanippJWBt23/IracPTAdeLWkDUrv+BPgusFtWVfXwNUlTSBkH20dlBrJ2fcy5Z0DSbNubS9oXeBXwUeC6Js1fdSJpS9JaWQFX2L4hs6TaKK19C2vbfwW+BKxoe11Jm5OSFg3kvOxItHI2SPoIMN/2SZJusN24REWSPksqjPNH2lZCRIa6YElYoSqluDvwDdtPSWrcXVaV0KTFndVj4THbD/RbU59ofPsW3LZHkxLY/BrA9uwGZx8EeErSPsC7WLS8s6krIt4GrGf7ydxCekE49zycTvoxvBG4QtI6QBPnZK8jRY23z8G2tg2sl0NUHyihfUtt2wW250qLhRU06satgwNIU0wn2P5TdSNzbmZNdXEjKUfFfbmF9IIYlh8jSFq+VWkqaB7Rvs2gqvl9KfAxUrnXw4AVbB+SVVgfqIoCrWV7Tm4tdVDVC9gU+B2Lz7kP5JRLOPcMVFmQppCKTgDMIM3bzR36VYOLUjdnX2Bd28dJWht4ke1rM0urhZLat8C2nQB8EngjaZTil8BxTcy1Dgsd3ltJo7yzSWl3Z9g+IqeuOpDUNYe87Rn91tILwrlnQNJ04GagtWRqf2Az213rgA86kk4lBajsaPvlVQ/gYttbZZZWCyW1b2ltWxqt4DlJB5J67VPaC6s0EUmTWDyJzUDGj8Scex7Wt71n2/YxkmZnU1M/21QRtzfAwrKgK+YWVSMltW8RbSvppwwztz6oQ7ejYPmqCt7bSSMWjUXSwcBxwHzSDetAx4+Ec8/DfEnb254JIGk70gXVVJ6SNI7qx1HSGjQgA9QwlNS+pbTtl3ILyMSxpKmHmbZ/J2k94PbMmuriKGBj2/fnFtILYlg+A9Xa2KnAKqS7wweAdzc4UGVfUiWtLUl2TwY+ZfuHWYXVREntW1rbjoSk6R2jNsGAIOkXpKqGj+XW0gvCuWekmtuhwalJF1LVOt+J5Owutf37EV4y8JTSviW27VA0LcGLpPGk+gEbA+Nb+22/N5uompC0BVXGRRpQNyCG5TMgaTVSNPX2gCXNJEVT/yOvsnqQ9DXgPNsn59bSD0pq39LadhQ0rbf0beA24E2kIfp9WVQgqWmcDlwG3EQDppai554BSb8CrmBRMoh9gR1s75xPVX1Iejdp6PafgfNJzqCxldNKat/S2nYkWulac+voFW3R8nNsb1plXvzloKZkHQ5JV9l+bW4dvSKcewYkXWf7VR37Ztl+dS5N/aBKWbon8A5S4Y2XZZZUCyW2byltOxINHJa/1vbWkq4ADgXuAa61PZAR5MMh6QTgz6SSr+3D8rEULhg1l0t6B/CDansycFFGPf1iA2Aj4KXArXml1EqJ7VtK247ER3ML6DFnVLkLPk2qCLcy8Jm8kmrjndXzx9v2DexSuOi5Z0DSw8BEFtXBHgc8Wv1t25OyCKsJSZ8H9iBVWzoPOL+h9aCBstq3wLa9iWfPq88FZgHHNzGuIhhMoueeAdvPHe64pI1t39IvPX3gT8C2Q60fbZq9hbVvUW0L/Jx00/bdavsd1fM84BwWVU5rBJJWJVWEeymLZ20byAjy4ajiCd7PorTRvwZOt/1UNlHLQPTcxyBNC8oZibC3uTTNVklX2t6u2z5JN9l+ZS5tdSDpKuAaOiLIbU8d8kUDiqSzSOVs29NGP237wHyqlp7ouY9NNPIpjSLsbS5Ns3VlSdvY/i2ApK1J89AATaz6N76JRWKGYCvbm7VtXybpxmxqlpFw7mOT0oZTwt7m0jRbDwS+JWll0o3LPOBASROBz2ZVVg/flnQQcCENiCAfgaclrW/7jwBVqt2nR3jNmCWcexAEwSix/TvglVVZX3UED/5giJcNMk8CXyQVjWndqA1sBPkIHEVa6XIH6cZtHeCAvJKWnnDuY5MncwvoM2Fvc2mErZL2s32upCM69gNg+ytZhNXPEcAGTSmmMhy2L5X0MmBDknO/zfYTI7xszLJcbgElIul9HdvjJE1pbdt+Tf9V1UeB9krSfpI+U22vXc3NAs2ytyBbJ1TPzx3i0VRuARpRSGUkJH0AWMn2HNs3AhMkHZpb19ISPfc87CRpT1JBhtVIxQpm5JVUK6XZewopsnhHUj7uh4HpwFY5RdVEKbauXz3fWljFu6eB2ZIupwHFVEbgoPYaCbYfrOINTsmoaakJ554B2++UtDdpecljwD62r8wsqzZKsxfYxvaWkm6AhT8SK+YWVROl2LqbpE+RspeV5Nz/u3qUwHKS5Gp9uKRxwMBey+HcM1DN6xxO6uG8HNi/ykndyOGv0uwFnqp+GFo/EmvQgCpTQ1CKrb8A7gcmSmov4SsalnWwnZHWszesfv0vgR9IOo10PR9CaveBJJLYZEDSbcAHqgAOkYJW3mt748zSaqFAe/clVUrbkpQQYzLwqSYO55ZkK4CkC2z/W24dY4UmFcqRtBxwMLAz6abtYuAs2wO5HC6cewYkTbI9r2Pfy2zfXv29i+1f5VHXe0qzF0DSRsBOpB+JS203tQZ2UbYGi9O0DITDMWijFOHcxyAl/cNAc+ytyp4OSZMSf5RkazuS9gA+D7yAdDPT6GH5kWjK/+5oGLRRiphzH5s0LWXnSDTF3utIc3Xt9rS2m5b4oyRb2/kC8K8xOrGQpvzvjoaB6gmHcx+bDNRF1AMaYa/tdXNr6Bcl2drBveHYF6Np9esbQzj3IOgxkroNU84F/my7UcVFSrK1Ypak80jLw9rXff84n6T6GEX9+ov7ryobAzVKEc59bHJnbgF95s7cAnrMKaTo8TmkH4RXAjcCq0k6pGE/iCXZCjCJlKvhjW37DDTSuVNY/foRGKhRigioy4CkCcCRwNq2D2rlM7Z9YWZptVCgvd8HjrN9S7X9ClJRiuOAH9vePKe+XlKSrSVSUv16SdsBR5MKxizPomDJgYwfidzyeTibNKS3bbX9F+D4fHJqpzR7N2o5OwDbtwJb2L4jo6a6KMlWJK0p6XxJ90m6V9J0SWvm1lUjK0vaprXR8Pr13wS+AmxPSp/8agY4jXIMy+dhfdt7S9oHwPZ8tcpLNZPS7P2DpFOB71fbewP/I+k5wFP5ZNVCSbZCulH9LrBXtb1ftW+XbIrqpaT69XNt/zy3iF4Rw/IZkHQVKenHlVVe7vWB79neeoSXDiQF2rsScCipByBgJmlu+nFggu1HMsrrKSXZCiBpdudUQ7d9TUPd69c3CkmfA8aR4ifagyWvzyZqGQjnngFJbwQ+CbyClOJwO+AA25dnFVYTpdkbNBdJl5ACyb5X7dqHdC3vlE1UDWiI+vUtmli/vqp814lt79h3MT0gnHsmJK0GvIbU27nG9v2ZJdVKSfZ2CcwBYFADc4ajJFsh1asHvkGKHzFwFXCY7buyCusxkg62fYakKd2O2z6m35qCJSPm3DMg6dLqTv+iLvsaR2n2kgJz/oOUxW0gi04sASXZCmkVwLttPwgL0/B+CXhvVlW9p7j69dXUwxTgddWuGcCxtufmU7X0hHPvI5LGAxOA1SU9j0VJESYBL8kmrCZKs7eNRgXmjEBJtgJs2nLskHLoSxqYfONLQIn1678F3Ay8vdrenxQsuUc2RctAOPf+8u/Ah0mO7ToWObt5wMm5RNVIafa2uFzSF2lIYM4IlGQrwHKSntfRc2/i72iJ9evX76j6doyk2dnULCMx554BSR+yfVJuHf2iQHsbFZgzHCXZCiDpXaTe7I9Ic+5vB06w/e2swmqipPr1kq4GjrI9s9reDviS7W2Hf+XYJJx7JiRtQooeH9/aZ3taPkX1Upq9QXOpsvDtyKL69bdmlhT0AEmbA1OBVUht+wDwHts3ZhW2lIRzz0AVgboDydn9DHgzMNP25Jy66qI0ewEkvQXYmMVvZo7Np6g+SrK1NEqsXy9pEoDteSOdO5Zp4lzRIDAZ2Ay4wfYBkl4InJVZU50UZa+k00iBhG8g2TkZuDarqJooydZCaXz9+qHW9LeSaA7qmv7ILZ+H+bafARZUd4n3AY1cF1xRmr2vtf0u4MFqPfC2wFqZNdVFSbaWSAn16ydWz8/t8lh5qBeNdaLnnodZklYFziRFkT9Cs3s7pdk7v3p+TNJLgH8A62bUUycl2Voija9fb/v06s9LbF/ZfqwKqhtIYs49M5JeCkyyPSezlL5Qgr2SPg2cRMqnfzIpqvpM25/JKqwGSrK1RCSd3WW3bTctaQ+Srre95Uj7BoVw7hnolp2tyRnbSrO3nao62vj2LFeSdrH9q4yyaqEkW4PmIGlb4LWknBwnth2aBLzN9mZZhC0jMefeRySNr5JerC7peZKeXz1eSgMztpVmbzdsP9ElfeXns4ipmZJsLYVC6tevSJpbX57F59vnkQJEB5KYc+8v3TK2GXiYVIyiaZRm72hpci37TkqytYk0vn697RnADEnn2P5zbj29InrufcT212yvC5wAbF79fTZwB3B1VnE1UJq9S0BJc2El2dpE1rB9tu0F1eMcYI3comrirCrwF4BqtPGXOQUtC+Hc8zDZ9jxJ25PugM8BTs0rqVZKszcImsL9kvaTNK567EdaEdFEVrf9UGujqh/wgox6lolw7nlolcZ8C3Ca7QtI8z5NpTR7R+LO3AL6yJ25BQTLxHtJ+fPvAf5GmoM+IKui+nhG0tqtjSo2aGBHniJaPgOSLgT+CuwMvIq0VvjaQY3KHIkC7Z0AHAmsbfsgSS8DNrR9YWZpPackW0tE0lTgw5316xu6FG5X4AxSHXdIdd0Ptj2QQ/PRc8/D24FfArtWw0DPB47KK6lWSrP3bFLCj1Y1qb8Ax+eTUysl2Voiz6pfDzSxfj22fwG8GvgDcB7ppnX+sC8aw0S0fAZsP0aqf93a/htpyKuRlGYvqS703pL2AbA9X61E1c2jJFtLpJT69Ug6EDgcWBOYDbyGFPg7kOWLG9lIQZCZJyWtRDVfJ2l92lJ3NoySbC2RLwNXSVqsfn1eSbVxOLAVcI3tN0jaCDgms6alJpx7EPSeo4FfAGtJ+g6wHc0NQjqacmwtDtvTJM1iUf36PRpcv/5x249LQtJzbN8macPcopaWCKgLghqQtBppWE+knsD9mSXVRkm2Bs1F0vmkG9MPk25mHgRWsL1bVmFLSTj3IOgxJeXSL8nWoBwkvR5YBfiF7Sdz61kaYlg+CHqEpPHABKpc+ixKvTqJhuXSL8nWoDyqlLQDTTj3IOgd3XLpQypAcXIuUTVRkq1BMHDEsHwQ9BhJH7J9Um4d/aAkW4NgkAjnHgQ1IGkT4BXA+NY+29PyKaqPkmwNgkEhnHsQ9BhJU4AdSA7vZ8CbgZm2B7Y29FCUZGsQDBKRfjYIes9kYCfgHtsHAJsBz8krqTZKsjUIBoZw7kHQe+bbfgZYIGkScB+wXmZNdVGSrUEwMES0fBD0nlmSVgXOJEWSPwJcm1dSbZRkaxAMDDHnHgQ1UtWEnmR7TmYptVOSrUEw1olh+SDoMZIubf1t+07bc9r3NYmSbA2CQSKG5YOgR5SUta0kW4NgEAnnHgS9o1vWNgMPA9/IqKsOSrI1CAaOGJYPgh5h+2u21yXVu968+vts4A7g6qziekxJtgbBIBLOPQh6z2Tb8yRtD+wCnAOcmldSbZRkaxAMDOHcg6D3PF09vwU4zfYFwIoZ9dRJSbYGwcAQzj0Ies9fJZ0OvB34maTn0Nz/tZJsDYKBIda5B0GPkTQB2BW4yfbtkl4MvNL2xZml9ZySbA2CQSKcexAEQRA0jBg+C4IgCIKGEc49CIIgCBpGOPcgCLoi6XhJlrRmbi1BECwZ4dyDMYWkHSqHMtTjNX3QcISkd9X9Ob1CicmSLpR0j6QnJT0k6UpJH6vSwxaNpAPbrqEduhzfoDr21QzygqDnRPrZYKzyPeBnXfb/bx8++wjgNmBaHz5rmZA0EfgBsBtwM3A68GdgZWBb4GjgrcBrM0kci3wB2Dq3iCCok3DuwVjletvn5hbRayStSFql8kSP3vJMkmP/HPAJL7785euSXgJ8oEef1QRmAVtJ2sv2D3OLCYK6iGH5YKCRtE81/PywpMckXSNpjy7nvVPSTyXdJekJSX+X9GNJm7Sds7wkA/8E7NQxHbBm67iks7q8f2vYd/u2fa05640kfVXSX4H5wFZt57xJ0q8kzZX0uKQbJR08Stu3BPYBZvJsxw6A7bttf7LjdZtLukDSA9Vn3iLpSEkj/h5IOlfSgi77n/XdtA11f0rSOyTNkTRf0u2taQ9J61Tt8KCkeZKmSVq522dKWlXS6VXbPS5ppqStOrWMwInAPcAJkkbs3Izmumk79y+SLpG0haTLJD0i6T5JX5A0TtJKkk6UdHf1Pfxa0oZd3md89Z3dWtn5oKSfSNqs47zlqna7qbr+50m6TdJZksYt4fcSNIzouQdjlQmSVu/Y94Tth1sbkj4HfJQ0fP9p4BlgT2C6pPfbPq3ttR8E7iUNW98LbAAcDFwlaQvbfySlUt0f+DpwN6k33OKBZbDl+8CjwJdIldPurfS/HzgZuAo4DngMeBNwuqR1bX98hPfds3o+s5tj74akbYDLgSeqz74X+LdK26bAu0dv1qjZnXTDdCrwIHAgMFXSU6Qh8ouBjwPbAO8hfQ+HdEoHfkVql6OBNUjTJxdJWs/2I6PU8hhwTKXlYOCUEc4fzXXTztqVPd8lTZfsChwFPAVsSfrN/SzwAuBI4HxJm9h+BhaO7FxM+i6mka7F5wEHVZ+5ve0bqs+aAnwGuKCyw8C6pPZcgUWpgYMSsR2PeIyZB7AD6Ueq2+P7bedtXe07tst7XAg8BExs2zexy3mbAE8CX+/Y/xfgki7nL1995lldjh1YHdu+bd/x1b5LgXEd569JcrDTurzXycACYJ0RvqsLqvffdAm+39+SHM0mbfsETK/e6/Vd9K/Ztu9cYMFovhuSIzTwCLBW2/4XVbY/AxzW8T4/qY6t1PGZ7tJO+1T73zcKu1vts3ul9Q+kHvzEDq1f7Xjdkl43Bt7Wsf/GytbpVInDqv1HVOfv1LbvqOrcnTveY9XO6xKYA8yp+38yHoP5iGH5YKxyBqnKWPvj+Lbj+5J+GKdJWr39QXIQq5B6PwDYfhQWRpZPqs67hxSgtw31cqLtzl7UXqQCK9/qov+nwDhgpxHed1L1PG80IpRSw24NnG/75tZ+2yb1JgHeNpr3WkKm2/6/ts9rfe8LeHYFud+Qvpd1urzPiR3bl1XPL1sSMbYXAJ8AXkjqPQ937pJeN3fZPr9j30zSDdRJ1Xfd4jdd9O8H3ALM7rgmlifdJL5eKX8/wFxgLUkRLBk8ixiWD8Yqt9u+ZJjjLyf9YN4+zDkvbP0h6VWkoe/XARM7P2tpRY6S/+my7+XV8+XDvO6FwxyDRU79uaPUsV71fEuXY7d0nNNL7uiy70Hgr7af6rIfYLWO/c+QVgG0848hzh0R29Ml/Rb4T0mnDXXeUlw3Q9kK8Kch9rfr34h0c/P3oTQBzwf+BnwMOB+4sorn+DVp1Gp6l+81KIxw7sGgItKc4m6kH/5u3Awg6aXAFaR582NJzvZRUs//JNL85GgYbl57uP+lx7rsU/W8L3DfEK8badnfzaRlblsAN41wbvtndmO0RSaGOm84+4ea+x1uTrhTq13NS4/i3NHyEWAGKV7ja89606W7boazaahj7fqXA2aThueH4gEA21dKWp8Up/EGYEfS9XRLNTf/0DDvETSccO7BoHI7sDPwJ9sj9bz3BCYAu9puDYUiScDqpOHNdro6MNtPS5pL6jl1sqQ93pbmv48wQjEc00nDywcyujX5reCvjbsce0X13K3n2c4DwDhJk2y3TwfU0eOvFdtXSLoI+HdSj7eTJb1uesHtpGDBSzuG8LviFGD6o+qBpMNINyoH8OxpjKAgYs49GFS+XT1/ttuyH0ntQ9qtHlNnD+8Q0o90J4/Q3YFD+vHdTtL4ts9ajSWPMj+PFJR1bPt7tb3nqlXk9JDYvp6U7OdflJbdPasHK+nFko6vzv8bcC2wu6SXt50jUrQ6pGHe4WhNMezcsf+IEV43VvkYKb7huC7HlvS66QXTSCsLDu92sP267rKaBOD66nmo6zcohOi5BwOJ7aslHUcaUr1e0o9I85AvBl4NvBFYqTr9IuC/gO9IOpnU49qOtEypcx4U4Brg3ZKOIUVVPwNcYHs+8A3gHOAySd8hLVM6uHqfFyyB/j9L+iBwGnCrpHOBu0i9tk1Jy5n+mRQhPRwHkSKpPwm8VdL06n1WJgV87cGiH3yAw0jz/DMlnUJa3vVWUsDiNNszRvi875ACG78paWPSvPFuDKgzsX2zpG/T/eZsSa+bXvAV0o3TiZJ2Js2jP0xaYrdT9fcu1bm3S/oN6YbtbuAlpGvxCdIyvKBgwrkHA4vtz0iaBXyI1HOcQHJWN5OcWOu82yXtBpxAcoILgCtJQVJnkJZmtfNxUrT9YdWzgLVIjnYa6QbiUNKw5/+SbjBWYAmj7m2fKen3wH8C768+635S6ttPMnxQVes9HpX0FmAyaY34oSRH+2j1PUypbGyd/1tJ25HWen+Q9J39kTTHO+Iwru2Hqu/yy5XGh0nTA/tX2geRTwN7A4uNoCzFdbPM2H5S0q6kttmP1E6QnPc1wNS2078IvJnUy1+FdO1fBXzW9mhiMIIGo1FM6wRBEARBMEDEnHsQBEEQNIxw7kEQBEHQMMK5B0EQBEHDCOceBEEQBA0jnHsQBEEQNIxw7kEQBEHQMMK5B0EQBEHDCOceBEEQBA0jnHsQBEEQNIxw7kEQBEHQMP4/Q1uC69t1UFkAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Note, \"get_feat_importance\" method is very slow. Can take up to ~6mins per 200 rows\n", + "# - Library used, LIME, uses local fitting to find key features behind each prediction (~2secs/pred) -\n", + "feat_importance = model.get_feat_importance(X, num_times = 200)\n", + "model.plot_feature_importance(feat_importance)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "reduced-freeware", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "logical-campaign", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "narrow-third", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/base.py b/base.py index 6a19857..fef18c4 100644 --- a/base.py +++ b/base.py @@ -3,9 +3,13 @@ import os import pathlib import pickle +import random import sys +from collections import OrderedDict, defaultdict from typing import List, Tuple, Union +import lime +import lime.lime_tabular import matplotlib import matplotlib.pyplot as plt import numpy as np @@ -14,8 +18,6 @@ from sklearn.metrics import auc, roc_curve from sklearn.preprocessing import StandardScaler -from collections import OrderedDict - from loaders import CsvReader matplotlib.rcParams["figure.figsize"] = [12, 10] @@ -219,14 +221,20 @@ def predict_sequentially_all(self, X, label_col_names: List[str] = None): if label_col_names is None: # If None provided, and None stored in self.labels, we ask user to provide as input # - Currently needed to match outputs to inputs when running the model forward - - raise ValueError("Please provide a list of predicted output labels ('label_col_names')") - + raise ValueError( + "Please provide a list of predicted output labels ('label_col_names')" + ) + # prepare features & a list of predictions that are feats too (often all) feats = self.features - preds_that_are_feats = [f_name for f_name in feats if f_name in label_col_names] + preds_that_are_feats = [ + f_name for f_name in feats if f_name in label_col_names + ] # initialize feat_dict to first row & pred_dict to match first row too feat_dict = OrderedDict(list(zip(feats, X[0]))) - pred_dict = dict([(k,v) for (k,v) in feat_dict.items() if k in preds_that_are_feats]) + pred_dict = dict( + [(k, v) for (k, v) in feat_dict.items() if k in preds_that_are_feats] + ) # sequentially iterate retriving next prediction based on previous prediction preds = [] @@ -242,31 +250,34 @@ def predict_sequentially_all(self, X, label_col_names: List[str] = None): # update prediction dictionary (for next iteration) pred_dict = OrderedDict(list(zip(label_col_names, pred.tolist()[0]))) - preds = np.array(preds) #.transpose() + preds = np.array(preds) # .transpose() if self.scale_data: preds = self.yscalar.inverse_transform(preds) - #preds_df = pd.DataFrame(preds) - #preds_df.columns = label_col_names + # preds_df = pd.DataFrame(preds) + # preds_df.columns = label_col_names - return preds #preds_df + return preds # preds_df - - def predict_sequentially(self, X, label_col_names: List[str] = None, it_per_episode: int = None): + def predict_sequentially( + self, X, label_col_names: List[str] = None, it_per_episode: int = None + ): if not self.model: raise ValueError("Please build or load the model first") else: if self.scale_data: X = self.xscalar.transform(X) - + if label_col_names is None: label_col_names = self.labels if label_col_names is None: # If None provided, and None stored in self.labels, we ask user to provide as input # - Currently needed to match outputs to inputs when running the model forward - - raise ValueError("Please provide a list of predicted output labels ('label_col_names')") + raise ValueError( + "Please provide a list of predicted output labels ('label_col_names')" + ) # initialize predictions preds = [] @@ -274,22 +285,22 @@ def predict_sequentially(self, X, label_col_names: List[str] = None, it_per_epis if not it_per_episode: it_per_episode = np.shape(X)[0] - num_of_episodes = int(np.shape(X)[0]/it_per_episode) + num_of_episodes = int(np.shape(X)[0] / it_per_episode) # iterate per as many episodes as selected for i in range(num_of_episodes): - X_aux = X[i*it_per_episode:(i+1)*it_per_episode] + X_aux = X[i * it_per_episode : (i + 1) * it_per_episode] preds_aux = self.predict_sequentially_all(X_aux, label_col_names) preds.extend(preds_aux) preds = np.array(preds) - #preds_df = pd.DataFrame(preds) - #preds_df.columns = label_col_names + # preds_df = pd.DataFrame(preds) + # preds_df.columns = label_col_names - return preds #preds_df + return preds # preds_df def predict_halt_classifier(self, X): @@ -385,22 +396,28 @@ def evaluate( return results_df def evaluate_sequentially( - self, X_test: np.ndarray, y_test: np.ndarray, metric, marginal: bool = False, it_per_episode = 100 + self, + X_test: np.ndarray, + y_test: np.ndarray, + metric, + marginal: bool = False, + it_per_episode=100, ): if not self.model: raise Exception("No model found, please run fit first") else: - + if not marginal: - y_hat = self.predict_sequentially(X_test, it_per_episode = it_per_episode) + y_hat = self.predict_sequentially(X_test, it_per_episode=it_per_episode) y_hat_len = np.shape(y_hat)[0] y_test = y_test[:y_hat_len] return metric(y_test, y_hat) else: - results_df = self.evaluate_margins_sequentially(X_test, y_test, metric, False, it_per_episode=it_per_episode) + results_df = self.evaluate_margins_sequentially( + X_test, y_test, metric, False, it_per_episode=it_per_episode + ) return results_df - def evaluate_margins( self, X_test: np.ndarray, y_test: np.ndarray, metric, verbose: bool = False @@ -418,11 +435,16 @@ def evaluate_margins( return pd.DataFrame(results.items(), columns=["var", "score"]) def evaluate_margins_sequentially( - self, X_test: np.ndarray, y_test: np.ndarray, metric, verbose: bool = False, it_per_episode: int = 100 + self, + X_test: np.ndarray, + y_test: np.ndarray, + metric, + verbose: bool = False, + it_per_episode: int = 100, ): # Extract prediction and remove any tail reminder from int(len(X_test)/it_per_episode) - y_pred = self.predict_sequentially(X_test, it_per_episode = it_per_episode) + y_pred = self.predict_sequentially(X_test, it_per_episode=it_per_episode) y_pred_len = np.shape(y_pred)[0] y_test = y_test[:y_pred_len] @@ -436,6 +458,122 @@ def evaluate_margins_sequentially( idx += 1 return pd.DataFrame(results.items(), columns=["var", "score"]) + def get_feat_importance(self, X=None, num_times=200): + + if not self.model: + raise Exception("No model found, please run fit first") + + if X is None: + raise Exception( + "Training data for feature importance has not been provided, but it should" + ) + + if not self.separate_models: + raise NotImplementedError( + "Lime doesn't support regression analysis with multiple outputs" + ) + + feats = self.features + outputs = self.labels + + # Initialize feature importance final dictionary + feat_importance_d = dict() + + # Identify as categorical all inputs with 10 or less unique values + categorical_features = np.argwhere( + np.array([len(set(X[:, i])) for i in range(X.shape[1])]) <= 10 + ).flatten() + + # Prepare lime predictor + explainer = lime.lime_tabular.LimeTabularExplainer( + X, + feature_names=feats, + categorical_features=categorical_features, + verbose=False, + mode="regression", + ) + + aux_X = X.copy() + random.shuffle(aux_X) + + for output, inner_model in zip(outputs, self.models): + + feats_coef_d = defaultdict(list) + + for i in range(100): + # Get explainer for single instance + exp = explainer.explain_instance( + aux_X[i], inner_model.predict, num_features=7 + ) + # exp_data is a list of tuples of type (range_str, importance) --> e.g.: ('state_x_position > 0.12', 0.671470753024958) + exp_data = exp.as_list() + + # Iterate extracting feature importance (order is given as per feat importance) + aux_idx_list = list(range(len(feats))) + feat_coefs = [] + for i, f in enumerate(feats): + + for idx in aux_idx_list: + feat_range_str, feat_coef = exp_data[idx] + + # Extract state name & append to list + if f in feat_range_str: + feat_coefs.append(abs(feat_coef)) + break + + # Remove state before running over again + aux_idx_list.remove(idx) + + # Normalize vector of importances + coef_sum = sum(feat_coefs) + feat_coefs = [coef / coef_sum for coef in feat_coefs] + + # Append to list + for f, c in zip(feats, feat_coefs): + feats_coef_d[f].append(c) + + # Get mean across all importances per example + feat_importance_d[output] = [np.mean(feats_coef_d[f]) for f in feats] + + return feat_importance_d + + def plot_feature_importance( + self, feature_data=None, in_feats_dim=None, out_labels=None, total_width=1.0 + ): + + if not self.model: + raise Exception("No model found, please run fit first") + + if feature_data is None: + feature_data = self.get_feat_importance() + if in_feats_dim is None: + in_feats_dim = self.input_dim + if out_labels is None: + out_labels = self.features + + fig, ax = plt.subplots(figsize=(8, 7)) + + colors = plt.rcParams["axes.prop_cycle"].by_key()["color"] + n_bars = len(feature_data) + 1 # plus 1 to add a space in between predictors + bar_width = total_width / n_bars # width of single bar + bars = [] + for i, (name, values) in enumerate(feature_data.items()): + x_offset = (i - n_bars / 2) * bar_width + bar_width / 2 + for x, y in enumerate(values): + bar = ax.bar( + x + x_offset, y, width=bar_width, color=colors[i % len(colors)] + ) + bars.append(bar[0]) + + ax.legend(bars, feature_data.keys()) + plt.title("Feature Importance / Feature Multipliers", fontsize=18) + plt.xlabel("Feature Column Names", fontsize=18) + plt.ylabel("Feature Importance", fontsize=18) + + plt.xticks(ticks=range(in_feats_dim), labels=out_labels) + ax.tick_params(labelrotation=90) + plt.show() + def plot_roc_auc(self, halt_x: np.ndarray, halt_y: np.ndarray): test_halt_preds = self.predict_halt_classifier(halt_x) diff --git a/environment.yml b/environment.yml index 67609b5..e97cf7d 100644 --- a/environment.yml +++ b/environment.yml @@ -1,39 +1,40 @@ name: ddm channels: - - pytorch - - defaults + - pytorch + - defaults dependencies: - - python=3.7.7 - - pip=19.1.1 - - pytorch=1.7.0 - - torchvision=0.8 - - cryptography=3.1.1 - - pip: - - ray==1.0.0 - - dataclasses==0.6 - - black==19.10b0 - - joblib==0.17.0 - - keras==2.4.3 - - scikit-learn==0.23.2 - - scikit-optimize==0.8.1 - - skorch==0.9.0 - - tensorboard==2.3.0 - - tensorflow==2.3.1 - - tensorflow-estimator==2.3.0 - - microsoft-bonsai-api==0.1.2 - - bonsai-cli==1.0.4 - - pandas==1.1.4 - - pyyaml==5.3.1 - - pytest==6.2.1 - - scipy==1.5.4 - - h5py==2.10.0 - - nbgrader==0.6.1 - - tune-sklearn==0.1.0 - - hyperopt==0.2.5 - - lightgbm==3.1.1 - - xgboost==1.3.0.post0 - - python-dotenv==0.15.0 - - hydra-core==1.0.5 - - natsort==7.1.0 - - seaborn==0.11.1 - - rich==9.13.0 + - python=3.7.7 + - pip=19.1.1 + - pytorch=1.7.0 + - torchvision=0.8 + - cryptography=3.1.1 + - pip: + - ray==1.0.0 + - dataclasses==0.6 + - black==19.10b0 + - joblib==0.17.0 + - keras==2.4.3 + - scikit-learn==0.23.2 + - scikit-optimize==0.8.1 + - skorch==0.9.0 + - tensorboard==2.3.0 + - tensorflow==2.3.1 + - tensorflow-estimator==2.3.0 + - microsoft-bonsai-api==0.1.2 + - bonsai-cli==1.0.4 + - pandas==1.1.4 + - pyyaml==5.3.1 + - pytest==6.2.1 + - scipy==1.5.4 + - h5py==2.10.0 + - nbgrader==0.6.1 + - tune-sklearn==0.1.0 + - hyperopt==0.2.5 + - lightgbm==3.1.1 + - xgboost==1.3.0.post0 + - python-dotenv==0.15.0 + - hydra-core==1.0.5 + - natsort==7.1.0 + - seaborn==0.11.1 + - lime==0.2.0.1 + - rich==9.13.0 diff --git a/gboost_models.py b/gboost_models.py index ec18305..dc89b75 100644 --- a/gboost_models.py +++ b/gboost_models.py @@ -195,6 +195,47 @@ def sweep(self, params: Dict, X, y): tune_search.fit(X, y) return tune_search + + + + def get_feat_importance(self, X = None): + + if not self.model: + raise Exception("No model found, please run fit first") + + feats = self.features + outputs = self.labels + + + feat_importance_d = dict() + + if not self.separate_models: + + if self.model_type in ["xgboost", "lightgbm"]: + for i,(feat_name,estimator) in enumerate(zip(outputs, self.model.estimators_)): + feat_importance_d[feat_name] = estimator.feature_importances_ + + else: + raise ValueError("Unknown model type") + + else: + + if self.model_type in ["xgboost", "lightgbm"]: + for i,feat_name in enumerate(outputs): + feat_importance_d[feat_name] = self.models[i].feature_importances_ + else: + raise ValueError("Unknown model type") + + + # Normalize vector of importances + for f_name in feat_importance_d.keys(): + feat_importances = np.abs(feat_importance_d[f_name]).tolist() + importance_sum = sum(feat_importances) + feat_importances_norm = [imp/importance_sum for imp in feat_importances] + feat_importance_d[f_name] = feat_importances_norm + + return feat_importance_d + if __name__ == "__main__": diff --git a/skmodels.py b/skmodels.py index 07b788c..f219684 100644 --- a/skmodels.py +++ b/skmodels.py @@ -183,6 +183,61 @@ def sweep(self, X, y, params: Dict = None): tune_search.fit(X, y) return tune_search + + + + def get_feat_importance(self, X = None): + + if not self.model: + raise Exception("No model found, please run fit first") + + if self.model_type == "SVR": + # Note, .coef_ is only available for linear models + return BaseModel.get_feat_importance(self, X) + + if self.model_type in ["linear_model"] and not self.scale_data: + # Note, coefficients do NOT match feature importance for NON scaled data (self.scale_data == False) + return BaseModel.get_feat_importance(self, X) + + feats = self.features + outputs = self.labels + + feat_importance_d = dict() + + if not self.separate_models: + + if self.model_type in ["GradientBoostingRegressor"]: + for i,(feat_name,estimator) in enumerate(zip(outputs, self.model.estimators_)): + feat_importance_d[feat_name] = estimator.feature_importances_ + + elif self.model_type in ["linear_model"]: + print("self.model.coef_", self.model.coef_) + for i,(feat_name,coefficients) in enumerate(zip(outputs, self.model.coef_)): + feat_importance_d[feat_name] = coefficients + + else: + raise ValueError("Unknown model type") + + else: + + if self.model_type in ["GradientBoostingRegressor"]: + for i,feat_name in enumerate(outputs): + feat_importance_d[feat_name] = self.models[i].feature_importances_ + elif self.model_type in ["linear_model"]: + for i,feat_name in enumerate(outputs): + feat_importance_d[feat_name] = self.models[i].coef_ + else: + raise ValueError("Unknown model type") + + # Normalize vector of importances + for f_name in feat_importance_d.keys(): + feat_importances = np.abs(feat_importance_d[f_name]).tolist() + importance_sum = sum(feat_importances) + feat_importances_norm = [imp/importance_sum for imp in feat_importances] + feat_importance_d[f_name] = feat_importances_norm + + return feat_importance_d + if __name__ == "__main__":