6_Reproducible_results_primer.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Intro\n",
    "This notebook is used to showcases which precautions to take to ensure you can iteratively develop Neural Networks and get reproducible results. Specifically for the case where you go back and forth between preprocessing and training cells, experimenting with differnt settings. For the complete write-up and explanation, please read my blogpost [here](https://towardsdatascience.com/primer-on-developing-reproducible-neural-networks-in-jupyter-notebook-2f88d23b7c8d?sk=dbb8d4efdc36aa676cb97857b946c5ca)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "seed = 42\n",
    "os.environ['PYTHONHASHSEED']=str(seed)\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns; sns.set()\n",
    "\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "from sklearn.model_selection import GroupShuffleSplit\n",
    "from sklearn.metrics import mean_squared_error, r2_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "from tensorflow.keras import Sequential\n",
    "from tensorflow.keras.layers import Dense"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## step 1: set random seeds to get reproducible results\n",
    "\n",
    "Initial weights of neural networks are generated randomly. Due to this randomness, results may vary across training cycles, making it very difficult to judge whether the changes in model performance are because of this randomness or actuall improvements. Therefore we wan't to control the random initialization to allow comparing between iterations.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.seed(seed)\n",
    "np.random.seed(seed)\n",
    "tf.random.set_seed(seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path = './CMAPSSData/'\n",
    "train_file = 'train_FD001.txt'\n",
    "test_file = 'test_FD001.txt'\n",
    "index_names = ['unit_nr', 'time_cycles']\n",
    "setting_names = ['setting_1', 'setting_2', 'setting_3']\n",
    "sensor_names = ['s_{}'.format(i+1) for i in range(0,21)]\n",
    "col_names = index_names + setting_names + sensor_names\n",
    "train = pd.read_csv((dir_path+train_file), sep='\\s+', header=None, \n",
    "                 names=col_names)\n",
    "test = pd.read_csv((dir_path+test_file), sep='\\s+', header=None, \n",
    "                 names=col_names)\n",
    "y_test = pd.read_csv((dir_path+'RUL_FD001.txt'), sep='\\s+', header=None, \n",
    "                 names=['RemainingUsefulLife'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(20631, 26)\n"
     ]
    },
    {
     "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>unit_nr</th>\n",
       "      <th>time_cycles</th>\n",
       "      <th>setting_1</th>\n",
       "      <th>setting_2</th>\n",
       "      <th>setting_3</th>\n",
       "      <th>s_1</th>\n",
       "      <th>s_2</th>\n",
       "      <th>s_3</th>\n",
       "      <th>s_4</th>\n",
       "      <th>s_5</th>\n",
       "      <th>...</th>\n",
       "      <th>s_12</th>\n",
       "      <th>s_13</th>\n",
       "      <th>s_14</th>\n",
       "      <th>s_15</th>\n",
       "      <th>s_16</th>\n",
       "      <th>s_17</th>\n",
       "      <th>s_18</th>\n",
       "      <th>s_19</th>\n",
       "      <th>s_20</th>\n",
       "      <th>s_21</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.0007</td>\n",
       "      <td>-0.0004</td>\n",
       "      <td>100.0</td>\n",
       "      <td>518.67</td>\n",
       "      <td>641.82</td>\n",
       "      <td>1589.70</td>\n",
       "      <td>1400.60</td>\n",
       "      <td>14.62</td>\n",
       "      <td>...</td>\n",
       "      <td>521.66</td>\n",
       "      <td>2388.02</td>\n",
       "      <td>8138.62</td>\n",
       "      <td>8.4195</td>\n",
       "      <td>0.03</td>\n",
       "      <td>392</td>\n",
       "      <td>2388</td>\n",
       "      <td>100.0</td>\n",
       "      <td>39.06</td>\n",
       "      <td>23.4190</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>0.0019</td>\n",
       "      <td>-0.0003</td>\n",
       "      <td>100.0</td>\n",
       "      <td>518.67</td>\n",
       "      <td>642.15</td>\n",
       "      <td>1591.82</td>\n",
       "      <td>1403.14</td>\n",
       "      <td>14.62</td>\n",
       "      <td>...</td>\n",
       "      <td>522.28</td>\n",
       "      <td>2388.07</td>\n",
       "      <td>8131.49</td>\n",
       "      <td>8.4318</td>\n",
       "      <td>0.03</td>\n",
       "      <td>392</td>\n",
       "      <td>2388</td>\n",
       "      <td>100.0</td>\n",
       "      <td>39.00</td>\n",
       "      <td>23.4236</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>-0.0043</td>\n",
       "      <td>0.0003</td>\n",
       "      <td>100.0</td>\n",
       "      <td>518.67</td>\n",
       "      <td>642.35</td>\n",
       "      <td>1587.99</td>\n",
       "      <td>1404.20</td>\n",
       "      <td>14.62</td>\n",
       "      <td>...</td>\n",
       "      <td>522.42</td>\n",
       "      <td>2388.03</td>\n",
       "      <td>8133.23</td>\n",
       "      <td>8.4178</td>\n",
       "      <td>0.03</td>\n",
       "      <td>390</td>\n",
       "      <td>2388</td>\n",
       "      <td>100.0</td>\n",
       "      <td>38.95</td>\n",
       "      <td>23.3442</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>0.0007</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>100.0</td>\n",
       "      <td>518.67</td>\n",
       "      <td>642.35</td>\n",
       "      <td>1582.79</td>\n",
       "      <td>1401.87</td>\n",
       "      <td>14.62</td>\n",
       "      <td>...</td>\n",
       "      <td>522.86</td>\n",
       "      <td>2388.08</td>\n",
       "      <td>8133.83</td>\n",
       "      <td>8.3682</td>\n",
       "      <td>0.03</td>\n",
       "      <td>392</td>\n",
       "      <td>2388</td>\n",
       "      <td>100.0</td>\n",
       "      <td>38.88</td>\n",
       "      <td>23.3739</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>-0.0019</td>\n",
       "      <td>-0.0002</td>\n",
       "      <td>100.0</td>\n",
       "      <td>518.67</td>\n",
       "      <td>642.37</td>\n",
       "      <td>1582.85</td>\n",
       "      <td>1406.22</td>\n",
       "      <td>14.62</td>\n",
       "      <td>...</td>\n",
       "      <td>522.19</td>\n",
       "      <td>2388.04</td>\n",
       "      <td>8133.80</td>\n",
       "      <td>8.4294</td>\n",
       "      <td>0.03</td>\n",
       "      <td>393</td>\n",
       "      <td>2388</td>\n",
       "      <td>100.0</td>\n",
       "      <td>38.90</td>\n",
       "      <td>23.4044</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 26 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   unit_nr  time_cycles  setting_1  setting_2  setting_3     s_1     s_2  \\\n",
       "0        1            1    -0.0007    -0.0004      100.0  518.67  641.82   \n",
       "1        1            2     0.0019    -0.0003      100.0  518.67  642.15   \n",
       "2        1            3    -0.0043     0.0003      100.0  518.67  642.35   \n",
       "3        1            4     0.0007     0.0000      100.0  518.67  642.35   \n",
       "4        1            5    -0.0019    -0.0002      100.0  518.67  642.37   \n",
       "\n",
       "       s_3      s_4    s_5  ...    s_12     s_13     s_14    s_15  s_16  s_17  \\\n",
       "0  1589.70  1400.60  14.62  ...  521.66  2388.02  8138.62  8.4195  0.03   392   \n",
       "1  1591.82  1403.14  14.62  ...  522.28  2388.07  8131.49  8.4318  0.03   392   \n",
       "2  1587.99  1404.20  14.62  ...  522.42  2388.03  8133.23  8.4178  0.03   390   \n",
       "3  1582.79  1401.87  14.62  ...  522.86  2388.08  8133.83  8.3682  0.03   392   \n",
       "4  1582.85  1406.22  14.62  ...  522.19  2388.04  8133.80  8.4294  0.03   393   \n",
       "\n",
       "   s_18   s_19   s_20     s_21  \n",
       "0  2388  100.0  39.06  23.4190  \n",
       "1  2388  100.0  39.00  23.4236  \n",
       "2  2388  100.0  38.95  23.3442  \n",
       "3  2388  100.0  38.88  23.3739  \n",
       "4  2388  100.0  38.90  23.4044  \n",
       "\n",
       "[5 rows x 26 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(train.shape)\n",
    "train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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>unit_nr</th>\n",
       "      <th>time_cycles</th>\n",
       "      <th>RUL</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>191</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>190</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>189</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>188</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>187</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   unit_nr  time_cycles  RUL\n",
       "0        1            1  191\n",
       "1        1            2  190\n",
       "2        1            3  189\n",
       "3        1            4  188\n",
       "4        1            5  187"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def add_remaining_useful_life(df):\n",
    "    # Get the total number of cycles for each unit\n",
    "    grouped_by_unit = df.groupby(by=\"unit_nr\")\n",
    "    max_cycle = grouped_by_unit[\"time_cycles\"].max()\n",
    "    \n",
    "    # Merge the max cycle back into the original frame\n",
    "    result_frame = df.merge(max_cycle.to_frame(name='max_cycle'), left_on='unit_nr', right_index=True)\n",
    "    \n",
    "    # Calculate remaining useful life for each row\n",
    "    remaining_useful_life = result_frame[\"max_cycle\"] - result_frame[\"time_cycles\"]\n",
    "    result_frame[\"RUL\"] = remaining_useful_life\n",
    "    \n",
    "    # drop max_cycle as it's no longer needed\n",
    "    result_frame = result_frame.drop(\"max_cycle\", axis=1)\n",
    "    return result_frame\n",
    "\n",
    "train = add_remaining_useful_life(train)\n",
    "train[index_names+['RUL']].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modelling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = train[sensor_names].copy()\n",
    "y_train = train['RUL'].copy()\n",
    "y_train_clipped = y_train.clip(upper=125)\n",
    "\n",
    "X_test = test.drop(['time_cycles']+setting_names, axis=1).groupby('unit_nr').last().copy()  # get the last row of each engine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = MinMaxScaler()\n",
    "scaler.fit(X_train[sensor_names])\n",
    "X_train_scaled = X_train.copy()\n",
    "X_train_scaled[sensor_names] = pd.DataFrame(scaler.transform(X_train[sensor_names]), columns=sensor_names)\n",
    "\n",
    "\n",
    "X_test_scaled = X_test.copy()\n",
    "X_test_scaled[sensor_names] = pd.DataFrame(scaler.transform(X_test[sensor_names]), columns=sensor_names, index=X_test.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train_split_engines [  2   3   4   6   7   8   9  10  12  14  15  16  17  18  20  21  22  24\n",
      "  25  26  27  28  29  30  33  35  36  37  38  39  41  42  43  44  47  48\n",
      "  49  50  51  52  53  55  56  57  58  59  60  61  62  63  64  65  66  67\n",
      "  68  69  70  72  73  75  76  79  80  82  83  85  86  87  88  89  90  92\n",
      "  93  94  95  96  97  98  99 100] \n",
      "\n",
      "validate_split_engines [ 1  5 11 13 19 23 31 32 34 40 45 46 54 71 74 77 78 81 84 91] \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# eventhough we set np and tf seeds, gss requires its own seed\n",
    "gss = GroupShuffleSplit(n_splits=1, train_size=0.80, random_state=42)  \n",
    "\n",
    "def train_val_group_split(X, y, gss, groups, print_groups=True):\n",
    "    \n",
    "    for idx_train, idx_val in gss.split(X, y, groups=groups):\n",
    "        if print_groups:\n",
    "            print('train_split_engines', train.iloc[idx_train]['unit_nr'].unique(), '\\n')\n",
    "            print('validate_split_engines', train.iloc[idx_val]['unit_nr'].unique(), '\\n')\n",
    "\n",
    "        X_train_split = X.iloc[idx_train].copy()\n",
    "        y_train_split = y.iloc[idx_train].copy()\n",
    "        X_val_split = X.iloc[idx_val].copy()\n",
    "        y_val_split = y.iloc[idx_val].copy()\n",
    "    return X_train_split, y_train_split, X_val_split, y_val_split\n",
    "\n",
    "split_result = train_val_group_split(X_train_scaled, y_train_clipped, gss, train['unit_nr'], print_groups=True)\n",
    "X_train_split_scaled, y_train_clipped_split_scaled, X_val_split_scaled, y_val_clipped_split_scaled = split_result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_cols = sensor_names  # + setting_names\n",
    "input_dim = len(train_cols)\n",
    "\n",
    "model = Sequential()\n",
    "model.add(Dense(16, input_dim=input_dim, activation='relu'))\n",
    "model.add(Dense(32, activation='relu'))\n",
    "model.add(Dense(64, activation='relu'))\n",
    "model.add(Dense(1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.compile(loss='mean_squared_error', optimizer='adam')\n",
    "model.save_weights('reproducible_primer_weights.h5')  # important, weights are saved after compiling the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: save model weights for later reuse\n",
    "After training a model the first time, you want to be able to reset it's state to the original values. This includes the optimizer parameters and the model weights.  \n",
    "The optimizer is reset by recompining the model, but recompiling generates a new set of random weights (from the seed you've set). Therefore, you need to load the original weights to get reproducible results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train on 16561 samples, validate on 4070 samples\n",
      "Epoch 1/20\n",
      "16561/16561 [==============================] - 2s 132us/sample - loss: 2301.6502 - val_loss: 421.9053\n",
      "Epoch 2/20\n",
      "16561/16561 [==============================] - 2s 92us/sample - loss: 432.1362 - val_loss: 364.5692\n",
      "Epoch 3/20\n",
      "16561/16561 [==============================] - 1s 63us/sample - loss: 397.5807 - val_loss: 329.9561\n",
      "Epoch 4/20\n",
      "16561/16561 [==============================] - 1s 57us/sample - loss: 379.3664 - val_loss: 317.8431\n",
      "Epoch 5/20\n",
      "16561/16561 [==============================] - 1s 64us/sample - loss: 370.3313 - val_loss: 309.5288\n",
      "Epoch 6/20\n",
      "16561/16561 [==============================] - 1s 61us/sample - loss: 364.8850 - val_loss: 303.8538\n",
      "Epoch 7/20\n",
      "16561/16561 [==============================] - 1s 58us/sample - loss: 359.3875 - val_loss: 298.7764\n",
      "Epoch 8/20\n",
      "16561/16561 [==============================] - 1s 59us/sample - loss: 356.1165 - val_loss: 299.1210\n",
      "Epoch 9/20\n",
      "16561/16561 [==============================] - 1s 56us/sample - loss: 353.3764 - val_loss: 292.9368\n",
      "Epoch 10/20\n",
      "16561/16561 [==============================] - 1s 56us/sample - loss: 352.1041 - val_loss: 301.1336\n",
      "Epoch 11/20\n",
      "16561/16561 [==============================] - 1s 58us/sample - loss: 351.0232 - val_loss: 291.9363\n",
      "Epoch 12/20\n",
      "16561/16561 [==============================] - 1s 56us/sample - loss: 347.4852 - val_loss: 289.1906\n",
      "Epoch 13/20\n",
      "16561/16561 [==============================] - 1s 57us/sample - loss: 344.7926 - val_loss: 288.5135\n",
      "Epoch 14/20\n",
      "16561/16561 [==============================] - 1s 55us/sample - loss: 343.5889 - val_loss: 292.0521\n",
      "Epoch 15/20\n",
      "16561/16561 [==============================] - 1s 55us/sample - loss: 345.5568 - val_loss: 280.9507\n",
      "Epoch 16/20\n",
      "16561/16561 [==============================] - 1s 57us/sample - loss: 341.0334 - val_loss: 281.2983\n",
      "Epoch 17/20\n",
      "16561/16561 [==============================] - 1s 62us/sample - loss: 340.5049 - val_loss: 294.6393\n",
      "Epoch 18/20\n",
      "16561/16561 [==============================] - 1s 68us/sample - loss: 341.5728 - val_loss: 278.9326\n",
      "Epoch 19/20\n",
      "16561/16561 [==============================] - 1s 73us/sample - loss: 337.1285 - val_loss: 280.6086\n",
      "Epoch 20/20\n",
      "16561/16561 [==============================] - 1s 59us/sample - loss: 336.5685 - val_loss: 283.4204\n"
     ]
    }
   ],
   "source": [
    "epochs = 20\n",
    "\n",
    "# reset initial optimizer parameters by recompiling model\n",
    "# and reset initial model weights to allow repeated running of cell with consistent results\n",
    "model.compile(loss='mean_squared_error', optimizer='adam')\n",
    "model.load_weights('reproducible_primer_weights.h5')\n",
    "\n",
    "history = model.fit(X_train_split_scaled[train_cols], y_train_clipped_split_scaled,\n",
    "                    validation_data=(X_val_split_scaled[train_cols], y_val_clipped_split_scaled),\n",
    "                    epochs=epochs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxMAAAE/CAYAAAA9q540AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdfXxU9Z33//eZm9wMCTcTJgylLHTdurjYgjdbtXaDtRcEhZSVqguyoKKXLEvVuhYviFYqrcKFrF7rIl3ruvaxSiuULmD5KdTWlmppV6WtFkupawkiwTBJuElCMpmZc35/zE1mkgET5MyZGV7PxwPnnPM958x3PpnEec8553sMy7IsAQAAAMAAuZzuAAAAAIDCRJgAAAAAcFoIEwAAAABOC2ECAAAAwGkhTAAAAAA4LYQJAAAAAKeFMAEAAADgtHic7oDdjhzpkGlyK40zpaqqQi0t7U53oyhRW/tQW/tQW3tQV/tQW/tQW/s4XVuXy9CwYYOythV9mDBNizBxhlFP+1Bb+1Bb+1Bbe1BX+1Bb+1Bb++RrbTnNCQAAAMBpIUwAAAAAOC1Ff5oTAAAAzj6WZam9/Zg6O9tlmjGnu/ORHD7skmmatj+Px1OiYcMCcrv7HxEIEwAAACg6R46EZBiG/P4Rcrs9MgzD6S6dNo/HpWjU3jBhWZY6Oo7ryJGQhg8f2e/tOM0JAAAARae7u0tDh1bJ4/EWdJDIFcMwNGjQYEWj3QPajjABAACAImTJMPioOxCnE7qoMAAAAIDTQpgAAAAAbNbe3q6lS7/a7/X/8Iffa+XKb9jYozODMGGjN/5wWE9t/b3T3QAAAIDD2tqO65139vZ7/XHj/kpLlnzNxh6dGYzmZKPG5g79YvcHmlv7lyrxup3uDgAAABzy//7fw2puDmnp0q9q//59GjJkqEpLS/Xgg6u0YsU3FAodVnNzSBdf/BktWfI1/eY3u/Qf//FtrVnzbS1c+L913nl/pTff/K2OHj2ir3xlsS677HKnX5IkwoStglU+SVLTkU6Nrq5wuDcAAABnr1/87pBefeuQLfv+3KdH6vJPnXo41a98ZbFuv32B7rjjn3TddV/U97//rxo58mN66aVt+uQnz9U3v/l/FYlE9Pd/f5327v1Dn+0jkaieeOJpvfrqz/Xkk98iTJwNgv54mPig9QRhAgAAAJKkYcP8GjnyY5KkyZOn6ve/360NG76rhoZ9OnbsmDo7T/TZ5pJLLpMk/fmfn6O2tuM57e+p2Bom1qxZoxdffFGSNGnSJN1zzz1av369nnnmGRmGofPPP18PPPCASkpKtGbNGv3gBz/Q4MGDJUnXX3+95syZoz179ujee+9VR0eHLr74Yj3wwAPyeAojA40YlggTLR0O9wQAAODsdvmnPvzoQa6UlpampjdufE4/+9nL+uIXr9G1135G+/a9K8uy+mxTUlIiKT58a7Z2p9h2AfbOnTv16quvatOmTdq8ebPefvttffvb39ZTTz2l5557Ts8//7xM09R3v/tdSdLu3bv1yCOPaMuWLdqyZYvmzJkjSVq8eLHuv/9+bd++XZZlacOGDXZ1+YwrLXFrWGWpPmjtmy4BAABw9nC73YrFYn2Wv/76f+uLX5ypKVOuUnd3t955548yTXvvdn0m2RYmAoGAlixZopKSEnm9Xp1zzjnq7u7WsmXLVFFRIcMwdO6556qxsVFSPEw88cQTqqur0/LlyxUOh3Xw4EF1dXVp4sSJkqSZM2dq27ZtdnXZFkG/jzABAABwlvP7qzRiRFAPPfRAxvLrr79BTz/9bc2b93f6l3/5Z51//qd16FCjQ70cOMPKwXGShoYGzZ49W9/73vc0duxYSVJra6uuvfZarVixQueff76+8pWvaMmSJRozZoyWLFmiUaNG6YorrtCqVav0ve99T5K0f/9+3Xbbbdq+fbvdXT5jvvWDN/WzX7+v5755NbdyBwAAyJG33/69PvaxMU53o+A0Nu7X+PF/1e/1bb/44J133tGCBQt0zz33pIJEU1OTbr31Vn3pS1/SJZdcIkl68sknU9vMnz9f9fX1qqmpyfgAblnWgD+Qt7S0yzSdO69siM+rE11RvdvQoiEVpR++QZ4LBCoVCrU53Y2iRG3tQ23tQ23tQV3tQ23tk2+1NU1T0WjhnC50Kh6PK2evxTTNPj9Hl8tQVVX2wYRsvWndrl27dNNNN+nuu+/WNddcI0l69913NWvWLF1zzTVatGiRJKmxsVEbN25MbWdZljwej4LBoEKhUGp5c3Ozqqur7ezyGTcybUQnAAAAoJjYFiYOHTqkRYsWafXq1Zo2bZqk+G3Eb7nlFt15552aP39+at2ysjI9/PDDOnDggCzL0rp16zR58mSNGjVKpaWl2rVrlyRpy5YtqqmpsavLtkgOD3uIMAEAAIAiY9tpTk899ZTC4bBWrlyZWnb11VerublZTz/9tJ5++mlJ0pVXXqk777xTy5cv18KFCxWJRHThhRfq5ptvliStXr1a9913n9rb2zV+/HjNmzfPri7bwj+kTF6PS02ECQAAABSZnFyA7SSnr5mQpPuf+m9VDS7TnddNcLQfZ0K+nQ9ZTKitfaitfaitPairfaitffKtth98sF/BYHFcgJ3Layay1c2xayYQx/CwAAAAKEaEiRwIVvkUOtqlaKw4RhQAAAAAJMJETowY5pNpWQod7XS6KwAAAMhzDz74db3wwg/V3BzSV796R9Z1Pve5i0+5j8bGg1qxYrkd3ctAmMiBYFVieNgWTnUCAABA/wwfHtDq1Y+d1rYffHBIBw++f4Z71JftN60D95oAAABwWuSPv1Bk789t2bf3L2vkPffyU65TX79YU6ZM1RVXfEGSNH/+3+v22+/St7+9VuFwl9ra2nXHHXfpb/7mitQ2hw416vbbF2jz5v9Phw41avnyr6mzs1Pjx5+fWicUOqwVK76h9vY2NTeHdPXVdbr11n/Qv/zLajU2HtQ///P/1d13/x8988x39NOfvqRYzNQll1yqhQvvGPDNoLPhyEQO+Mq8Guzzcq8JAACAs1Rt7dX68Y+3S5IOHHhP3d3d+sEP1mvJkq/pP/5jnZYsuU9PPvmtk27/6KOrdPXVdfrOd76rT32qZ4TQl17arsmTa/Xtb39H//mf67Vhw/d09OhR3XnnV/WXf3me7r77/+hXv9qpvXv36Mkn/1NPP71OoVBIP/rRi2fkdXFkIkcY0QkAAMA53nMv/9CjB3b67Gc/p0cfXaUTJzr04x9vV23tVbr++hu0c+cr+ulPf6y33/6dOjtPfn3tb36zS1//+oOSpClTrtLKld+QJN1ww1z9+tdv6LvffUb79r2raDSirq7M/bzxxmv6/e9365Zb5kqSwuEujRgRPCOvizCRI8Eqn379x2anuwEAAAAHeL1eXX753+jVV3+ul19+SQ8//C9atOh/68ILL9IFF1ykiy76az3wwH2n2IORuneaYRhyudySpH/910fV2HhQkydPVU3NFXrjjdfU+zZyphnT9dfP1qxZfy9Jamtrk9vtPiOvi9OcciToH6T2zojaOyNOdwUAAAAOqK29Ws8996yGDBkqn8+nAwf265Zb/kGXXnq5Xnllh0zz5LcRuPjiz2j79hckSTt2vKzu7rAk6Y03/ls33DBXV175v/Tee/sVCh2WaZpyuz2KxWKSpAsv/Gtt3/6CTpw4oWg0qqVL79bPfvaTM/KaODKRI8HERdhNrSdUMWqIw70BAABArn360xPV3t6uv/3bazV48BBNnz5Dc+deL4/Howsv/Gt1dXWd9FSnf/qne/SNb9yv55/fpHHjzpPPN0iS9Pd/f5O+8Y37VVpaqurqoMaN+ys1Nh7Uuef+pdrb2/SNb3xNX/vaN/Q///NH3XbbTTLNmC655LO66qrpZ+Q1GVbv4yBFpqWlPXVIyEkftJ5Q/bd/pVumnafLPzXS6e6ctkCgUqFQm9PdKErU1j7U1j7U1h7U1T7U1j75VtsPPtivYHCM0904Izwel6LR3Nz8OFvdXC5DVVUVWdfnNKccGT6kTG6XwUXYAAAAKBqEiRzxuF0KDC3nxnUAAAAoGoSJHGJ4WAAAgFwxZFm5OTWoWJzO1Q+EiRwK+n1qOtKZF9dwAAAAFLOSkjIdPdqsaDRyWh+SzzaWZamj47g8npIBbcdoTjkUrPIpGjPVfLxL1UPLne4OAABA0Ro2LKD29mNqbW2Sacac7s5H4nK5Tjls7Jni8ZRo2LDAwLaxqS/IIjk87ActJwgTAAAANjIMQ5WVQ1VZOdTprnxk+TZSVjpOc8qhVJjgugkAAAAUAcJEDlX6vPKVeggTAAAAKAqEiRwyDEPBKp+aCBMAAAAoAoSJHGN4WAAAABQLW8PEmjVrNG3aNE2bNk2rVq2SJO3cuVN1dXWaMmWKHn300dS6e/bs0cyZM1VbW6t7771X0WhUktTY2Kg5c+Zo6tSpWrhwoTo6Ouzssu2Cfp+OtIXV1R11uisAAADAR2JbmNi5c6deffVVbdq0SZs3b9bbb7+trVu3qr6+XmvXrtULL7yg3bt3a8eOHZKkxYsX6/7779f27dtlWZY2bNggSXrggQd0ww03aNu2bTr//PO1du1au7qcE8mLsJtaOx3uCQAAAPDR2BYmAoGAlixZopKSEnm9Xp1zzjlqaGjQmDFjNHr0aHk8HtXV1Wnbtm06ePCgurq6NHHiREnSzJkztW3bNkUiEb3++uuqra3NWF7IglXxMHGotbCPsAAAAAC2hYlPfvKTqXDQ0NCgF198UYZhKBDouRFGdXW1mpqadPjw4YzlgUBATU1NOnLkiCoqKuTxeDKWF7LqoeUyFL/XBAAAAFDIbL9p3TvvvKMFCxbonnvukdvtVkNDQ6rNsiwZhiHTNGUYRp/lycd0vec/TFVVxUfqvx0Cfp+OdkQUCFQ63ZXTUqj9LgTU1j7U1j7U1h7U1T7U1j7U1j75Wltbw8SuXbt0xx13qL6+XtOmTdNrr72mUCiUag+FQqqurlYwGMxY3tzcrOrqavn9frW1tSkWi8ntdqfWH4iWlnaZpnXGXtOZUD2kTA2HjuXtnQxPJZ/vwFjoqK19qK19qK09qKt9qK19qK19nK6ty2Wc9At6205zOnTokBYtWqTVq1dr2rRpkqQJEyZo37592r9/v2KxmLZu3aqamhqNGjVKpaWl2rVrlyRpy5Ytqqmpkdfr1cUXX6wXXnhBkrR582bV1NTY1eWcCfp9amrtlGXlV8gBAAAABsK2IxNPPfWUwuGwVq5cmVo2a9YsrVy5UrfffrvC4bAmTZqkqVOnSpJWr16t++67T+3t7Ro/frzmzZsnSVq2bJmWLFmib33rWxo5cqQeeeQRu7qcM8Eqn8KRmI60heUfXOZ0dwAAAIDTYlhF/vV4Pp7m9PuGVq1+7rf66qyJ+quxfqe7MyBOH2YrZtTWPtTWPtTWHtTVPtTWPtTWPk7X1pHTnHByPfeaYEQnAAAAFC7ChAOGVZaq1OvWIcIEAAAAChhhwgGGYWiEv1wfECYAAABQwAgTDgn6fdy4DgAAAAWNMOGQoN+nlmNdikRjTncFAAAAOC2ECYcE/T5ZkpqOdDrdFQAAAOC0ECYcEqyKj+jEqU4AAAAoVIQJh4wYlggTXIQNAACAAkWYcEh5qUdDK0oIEwAAAChYhAkHBf0+wgQAAAAKFmHCQcGqQfqg5YQsy3K6KwAAAMCAESYcFPT7dCIcVVtnxOmuAAAAAANGmHBQ0M+ITgAAAChchAkHpYaH5boJAAAAFCDChIOGDy6Tx+0iTAAAAKAgESYc5HIZGjGsnNOcAAAAUJAIEw4bwfCwAAAAKFCECYcF/T6FjnYqGjOd7goAAAAwIIQJhwX9PsVMS83HupzuCgAAADAghAmHpUZ04roJAAAAFBjChMNS95rgugkAAAAUGMKEwyrKvaoo9xImAAAAUHA8dj9Be3u7Zs2apX/7t3/Tu+++q0ceeSTV1tTUpAkTJuiJJ57QmjVr9IMf/ECDBw+WJF1//fWaM2eO9uzZo3vvvVcdHR26+OKL9cADD8jjsb3bORWsYkQnAAAAFB5bP5W/+eabuu+++9TQ0CBJmjRpkiZNmiRJCoVCmj17tpYuXSpJ2r17tx555BFdcMEFGftYvHixvvnNb2rixImqr6/Xhg0bdMMNN9jZ7ZwL+n16690Wp7sBAAAADIitpzlt2LBBy5YtU3V1dZ+2VatWadasWRo7dqykeJh44oknVFdXp+XLlyscDuvgwYPq6urSxIkTJUkzZ87Utm3b7OyyI0b6fTre0a0TXVGnuwIAAAD0m61HJh588MGsyxsaGvTaa6+l2js6OnTeeedp8eLFGjNmjJYsWaK1a9fqiiuuUCAQSG0XCATU1NQ0oD5UVVWc/gvIkXM/USX97F2FLWlMoNLp7nyoQAH0sVBRW/tQW/tQW3tQV/tQW/tQW/vka20dufhg/fr1uuGGG1RSUiJJGjRokJ588slU+/z581VfX6+amhoZhpFabllWxnx/tLS0yzStM9Nxm/g88de0592QhpXn9/UggUClQqE2p7tRlKitfaitfaitPairfaitfaitfZyurctlnPQLekdGc/rJT36iq6++OjXf2NiojRs3puYty5LH41EwGFQoFEotb25uznrKVKELDC2XyzC4CBsAAAAFJedhorW1VV1dXRo9enRqWVlZmR5++GEdOHBAlmVp3bp1mjx5skaNGqXS0lLt2rVLkrRlyxbV1NTkusu287hdGj60jBvXAQAAoKDk/Jya999/X8FgMGOZ3+/X8uXLtXDhQkUiEV144YW6+eabJUmrV6/Wfffdp/b2do0fP17z5s3LdZdzIuhneFgAAAAUlpyEiZdffjk1/elPf1obNmzos05tba1qa2v7LB83blzGKVDFKuj3ac/+IzItS64BXhcCAAAAOIE7YOeJYJVPkaip1uNdTncFAAAA6BfCRJ4Y6fdJEqc6AQAAoGAQJvJEMBkmuAgbAAAABYIwkScGDypReambIxMAAAAoGISJPGEYBiM6AQAAoKAQJvIIYQIAAACFhDCRR4J+n1qPhxXujjndFQAAAOBDESbyyIjERdhNRzg6AQAAgPxHmMgjQYaHBQAAQAEhTOSREQwPCwAAgAJCmMgjpV63qgaXcmQCAAAABYEwkWcY0QkAAACFgjCRZ4L+Qfqg9YQsy3K6KwAAAMApESbyTLDKp67umI51dDvdFQAAAOCUCBN5JshF2AAAACgQhIk8w/CwAAAAKBSEiTwzbHCpSjwuwgQAAADyHmEiz7gMQyMY0QkAAAAFgDCRh0b4fVwzAQAAgLxHmMhDQb9PoWOdikRNp7sCAAAAnBRhIg+N9PtkWdLho51OdwUAAAA4KdvDRHt7u6ZPn673339fkrR06VJNmTJFM2bM0IwZM/TSSy9Jknbu3Km6ujpNmTJFjz76aGr7PXv2aObMmaqtrdW9996raDRqd5cdF6yKj+jUxHUTAAAAyGO2hok333xTs2fPVkNDQ2rZ7t279eyzz2rLli3asmWLJk+erK6uLtXX12vt2rV64YUXtHv3bu3YsUOStHjxYt1///3avn27LMvShg0b7OxyXmB4WAAAABQCW8PEhg0btGzZMlVXV0uSOjs71djYqPr6etXV1emxxx6TaZp66623NGbMGI0ePVoej0d1dXXatm2bDh48qK6uLk2cOFGSNHPmTG3bts3OLueF8lKPhgwq4SJsAAAA5DWPnTt/8MEHM+abm5t16aWXatmyZaqsrNSCBQu0ceNG+Xw+BQKB1HrV1dVqamrS4cOHM5YHAgE1NTXZ2eW8EWR4WAAAAOQ5W8NEb6NHj9bjjz+emp87d642b96s2tpaGYaRWm5ZlgzDkGmaWZcPRFVVxUfvuAPGjhqinW8dUiBQ6XRX+sjHPhULamsfamsfamsP6mofamsfamuffK1tTsPE3r171dDQoNraWknxcODxeBQMBhUKhVLrhUIhVVdX91ne3NycOmWqv1pa2mWa1pl5ATk01OdV24lu7XuvVRXlXqe7kxIIVCoUanO6G0WJ2tqH2tqH2tqDutqH2tqH2trH6dq6XMZJv6DP6dCwlmXpoYce0rFjxxSJRLR+/XpNnjxZEyZM0L59+7R//37FYjFt3bpVNTU1GjVqlEpLS7Vr1y5J0pYtW1RTU5PLLjsmdRE2100AAAAgT+X0yMS4ceN02223afbs2YpGo5oyZYqmT58uSVq5cqVuv/12hcNhTZo0SVOnTpUkrV69Wvfdd5/a29s1fvx4zZs3L5dddkxyeNhDrR36i48Pcbg3AAAAQF+GZVmFdw7QABTqaU4x09Q/rN6hKZ8Zreuu+Aunu5Pi9GG2YkZt7UNt7UNt7UFd7UNt7UNt7eN0bfPmNCf0n9vlUvWwck5zAgAAQN4iTOQxhocFAABAPiNM5LGg36fDRzoL8jQtAAAAFD/CRB4L+n2KmZaaj3U63RUAAACgD8JEHkuO6MSpTgAAAMhHhIk8xr0mAAAAkM8IE3ms0leiQWUejkwAAAAgLxEm8lywihGdAAAAkJ8IE3ku6PfpEGECAAAAeYgwkeeCfp+OtXerMxx1uisAAABABsJEngv6B0liRCcAAADkn36FiebmZv3kJz+RJD388MO68cYb9Yc//MHWjiEu6C+XRJgAAABA/ulXmFiyZIkOHDigX/7yl3rllVc0Y8YMffOb37S7b5BUPcwnw5CaCBMAAADIM/0KE0ePHtVNN92kn//855o+fbpmzpypzk7uypwLXo9Lw4eUcWQCAAAAeadfYSISiSgSieiVV17RZz/7WXV2durECT7c5krQP4gb1wEAACDv9CtMfOELX9Bll12mYcOG6fzzz9d1112n6dOn2903JAT9Pn1w5IRMy3K6KwAAAECKpz8r3XHHHbr++us1YsQISdLq1as1btw4WzuGHsEqn7ojpo62heUfXOZ0dwAAAABJAxjN6e2335ZhGHr44Ye1YsUKRnPKoaDfJ0ncvA4AAAB5hdGcCkAyTHDdBAAAAPIJozkVgKEVJSotcTOiEwAAAPIKozkVAMMw4hdhEyYAAACQRxjNqUCM9Ps4zQkAAAB5pV9h4o477tDWrVv1zDPPSIqP5rRo0aJ+PUF7e7umT5+u999/X5K0fv16TZ8+XXV1dVq6dKm6u7slSWvWrNHnP/95zZgxQzNmzNC6deskSXv27NHMmTNVW1ure++9V9FodMAvshiM8PvUerxL3ZGY010BAAAAJPUzTJimqR/+8IeaO3euZs+erR//+Mf9+lD/5ptvavbs2WpoaJAk7du3T0899ZSee+45Pf/88zJNU9/97nclSbt379YjjzyiLVu2aMuWLZozZ44kafHixbr//vu1fft2WZalDRs2nOZLLWxBv0+WpMNHuFYFAAAA+aFfYeKf//mf9atf/Uo33nijbr75Zv3mN7/RqlWrPnS7DRs2aNmyZaqurpYklZSUaNmyZaqoqJBhGDr33HPV2NgoKR4mnnjiCdXV1Wn58uUKh8M6ePCgurq6NHHiREnSzJkztW3bttN9rQUtNaIT100AAAAgT/TrpnWvvPKKfvCDH8jr9UqSrrjiCn3xi19UfX39Kbd78MEHM+ZHjRqlUaNGSZJaW1u1bt06rVixQh0dHTrvvPO0ePFijRkzRkuWLNHatWt1xRVXKBAIpLYPBAJqamoa0AusqqoY0Pr5qnJwuSTpeDiqQKDS0b44/fzFjNrah9rah9rag7rah9rah9raJ19r268wYVlWKkhI8SMM6fMD1dTUpFtvvVVf+tKXdMkll0iSnnzyyVT7/PnzVV9fr5qaGhmGkdGP9Pn+aGlpl2lap93XfDKsslTvvndUoVCbY30IBCodff5iRm3tQ23tQ23tQV3tQ23tQ23t43RtXS7jpF/Q9+s0p3Hjxumhhx7Se++9pwMHDmjFihU699xzT6sz7777rmbNmqVrrrkmdRF3Y2OjNm7cmFrHsix5PB4Fg0GFQqHU8ubm5tQpU2cjhocFAABAPulXmFi2bJmOHz+uWbNm6frrr1dLS4tmz5494Cdrb2/XLbfcojvvvFPz589PLS8rK9PDDz+sAwcOyLIsrVu3TpMnT9aoUaNUWlqqXbt2SZK2bNmimpqaAT9vsQhWxcOEZRXHkRYAAAAUtn6d5lRRUaGVK1dmLLvwwgv161//ekBPtnHjRjU3N+vpp5/W008/LUm68sordeedd2r58uVauHChIpGILrzwQt18882S4sPQ3nfffWpvb9f48eM1b968AT1nMQn6feoMR3X8RERDBpU43R0AAACc5foVJrIZyLfjL7/8siTppptu0k033ZR1ndraWtXW1vZZPm7cuIxToM5mI5MjOrV0ECYAAADguH6d5pTNQC+ExkfH8LAAAADIJ6cdJpB7/iFl8rhdhAkAAADkhVOe5nTBBRdkPQJhWZa6urps6xSycxmGRvjL1dTKXbABAADgvFOGia1bt+aqH+inoN+n90MdTncDAAAAOHWYSN6tGvkj6Pfpt+80Kxoz5XFzlhoAAACcw6fRAhP0+xQzLYWOcqoTAAAAnEWYKDDBKkZ0AgAAQH4gTBSYkQwPCwAAgDxBmCgwvjKvBvu8+qCFMAEAAABnESYKUNDv48gEAAAAHEeYKEDBKsIEAAAAnEeYKEBB/yC1nYiooyvidFcAAABwFiNMFKBg8iJsrpsAAACAgwgTBWiEv1wSIzoBAADAWYSJAhQYWi63yyBMAAAAwFGEiQLkcbs0fGg5YQIAAACOIkwUqJEMDwsAAACHESYKVNDvU1Nrp0zTcrorAAAAOEsRJgpUsMqnaMxUy/Eup7sCAACAsxRhokClhoflVCcAAAA4hDBRoLjXBAAAAJxGmChQlT6vfKUejkwAAADAMbaGifb2dk2fPl3vv/++JGnnzp2qq6vTlClT9Oijj6bW27Nnj2bOnKna2lrde++9ikajkqTGxkbNmTNHU6dO1cKFC9XR0WFndwuKYRgKVjGiEwAAAJxjW5h48803NXv2bDU0NEiSurq6VF9fr7Vr1+qFF17Q7t27tWPHDknS4sWLdf/992v79u2yLEsbNmyQJD3wwAO64YYbtG3bNi3r878AACAASURBVJ1//vlau3atXd0tSEGGhwUAAICDbAsTGzZs0LJly1RdXS1JeuuttzRmzBiNHj1aHo9HdXV12rZtmw4ePKiuri5NnDhRkjRz5kxt27ZNkUhEr7/+umprazOWo8cIv09H2sIKd8ec7goAAADOQh67dvzggw9mzB8+fFiBQCA1X11draampj7LA4GAmpqadOTIEVVUVMjj8WQsH6iqqorTfAX57y/HVkk//5PClvTxQGXOnjeQw+c621Bb+1Bb+1Bbe1BX+1Bb+1Bb++RrbW0LE72ZpinDMFLzlmXJMIyTLk8+pus93x8tLe1Fe2M3nydejz3vNmtwqTsnzxkIVCoUasvJc51tqK19qK19qK09qKt9qK19qK19nK6ty2Wc9Av6nI3mFAwGFQqFUvOhUEjV1dV9ljc3N6u6ulp+v19tbW2KxWIZ66NH9bByGeJeEwAAAHBGzsLEhAkTtG/fPu3fv1+xWExbt25VTU2NRo0apdLSUu3atUuStGXLFtXU1Mjr9eriiy/WCy+8IEnavHmzampqctXdglDidatqSBlhAgAAAI7I2WlOpaWlWrlypW6//XaFw2FNmjRJU6dOlSStXr1a9913n9rb2zV+/HjNmzdPkrRs2TItWbJE3/rWtzRy5Eg98sgjuepuwQj6fdy4DgAAAI6wPUy8/PLLqenLLrtMzz//fJ91xo0bp40bN/ZZPmrUKD3zzDO29q/QBf0+vXPwUNZrTAAAAAA7cQfsAhes8incHdPR9m6nuwIAAICzDGGiwAX9PknSBy3cHRwAAAC5RZgocKkwwUXYAAAAyDHCRIEbVlmqEq9LH7R2Ot0VAAAAnGUIEwXOMAwFh/k4MgEAAICcI0wUgWCVTx+0cs0EAAAAcoswUQSCfp+aj3UpEjWd7goAAADOIoSJIhD0+2RZ0uEjnOoEAACA3CFMFIFgFSM6AQAAIPcIE0VgxDDCBAAAAHKPMFEEyks9GlpRog9aCBMAAADIHcJEkQj6GR4WAAAAuUWYKBLBqkH6oPWELMtyuisAAAA4SxAmikTQ71NHV1TtnRGnuwIAAICzBGGiSAT9XIQNAACA3CJMFImgv1ySuAgbAAAAOUOYKBLDh5TL4zY4MgEAAICcIUwUCZfLUPUwRnQCAABA7hAmigjDwwIAACCXCBNFJOj36fCRTsVM0+muAAAA4CxAmCgiQb9PMdNS89Eup7sCAACAs4An10/4/e9/X88++2xq/v3339eMGTPU2dmpXbt2qbw8PirRl7/8ZU2ePFk7d+7UihUrFA6HddVVV+muu+7KdZcLRrAqPjzsodYTGpEYKhYAAACwS87DxHXXXafrrrtOkvTOO+9o0aJF+vKXv6wbb7xRzz77rKqrq1PrdnV1qb6+Xs8884xGjhypBQsWaMeOHZo0aVKuu10QUveaaDkh/YXDnQEAAEDRc/Q0p69//eu66667VF5ersbGRtXX16uurk6PPfaYTNPUW2+9pTFjxmj06NHyeDyqq6vTtm3bnOxyXqso96qi3MtF2AAAAMiJnB+ZSNq5c6e6urp01VVX6cCBA7r00ku1bNkyVVZWasGCBdq4caN8Pp8CgUBqm+rqajU1NQ3oeaqqKs501/Pa6BGVam0PKxCotO057Nz32Y7a2ofa2ofa2oO62ofa2ofa2idfa+tYmHjuued08803S5JGjx6txx9/PNU2d+5cbd68WbW1tTIMI7XcsqyM+f5oaWmXaVpnptMFoGpwqX73botCoTZb9h8IVNq277MdtbUPtbUPtbUHdbUPtbUPtbWP07V1uYyTfkHvyGlO3d3dev3113XllVdKkvbu3avt27en2i3LksfjUTAYVCgUSi0PhUIZ11Sgr6Dfp2Md3eoMR53uCgAAAIqcI2Fi7969Gjt2rHy++AXDlmXpoYce0rFjxxSJRLR+/XpNnjxZEyZM0L59+7R//37FYjFt3bpVNTU1TnS5YKQuwua6CQAAANjMkdOcDhw4oGAwmJofN26cbrvtNs2ePVvRaFRTpkzR9OnTJUkrV67U7bffrnA4rEmTJmnq1KlOdLlgpI/o9ImRgx3uDQAAAIqZI2Hi6quv1tVXX52xbM6cOZozZ06fdS+77DI9//zzuepawaseVi6XYegQRyYAAABgM+6AXWQ8bpeGDy3jNCcAAADYjjBRhIJ+X/zGdQAAAICNCBNFKOj36fCREzKts2dIXAAAAOQeYaIIBat86o6aaj3e5XRXAAAAUMQIE0VoZGJEp6bWTod7AgAAgGJGmChC3GsCAAAAuUCYKEKDB5WovNTNRdgAAACwFWGiCBmGoRHDfPqgtcPprgAAAKCIESaKVLDKx2lOAAAAsBVhokgF/T61HA8rHIk53RUAAAAUKcJEkQqmRnTi6AQAAADsQZgoUozoBAAAALsRJorUCMIEAAAAbEaYKFKlXreqBpcSJgAAAGAbwkQRC/p93GsCAAAAtiFMFLGgf5CajpyQZVlOdwUAAABFiDBRxIJVPnWGYzre0e10VwAAAFCECBNFjBGdAAAAYCfCRBEb4S+XJB0iTAAAAMAGhIki5h9cphKPi4uwAQAAYAvCRBFzGYaqh/k4zQkAAAC2IEwUuWAVYQIAAAD28DjxpHPnzlVra6s8nvjTL1++XO+9956+9a1vKRqN6sYbb9ScOXMkSTt37tSKFSsUDod11VVX6a677nKiywUr6Pfp13tDisZMedxkRwAAAJw5OQ8TlmWpoaFBP/3pT1NhoqmpSXfddZf+67/+SyUlJZo1a5YuueQSffzjH1d9fb2eeeYZjRw5UgsWLNCOHTs0adKkXHe7YI30+2Ralg4f6dTHhg9yujsAAAAoIjkPE3/6058kSfPnz9fRo0d1/fXXa9CgQbr00ks1dOhQSVJtba22bdumz3zmMxozZoxGjx4tSaqrq9O2bdsIEwMQrOoZHpYwAQAAgDMp52Hi+PHjuuyyy/S1r31NkUhE8+bN01VXXaVAIJBap7q6Wm+99ZYOHz7cZ3lTU9OAnq+qquKM9b0QDaoskyS1h2MKBCrPyD7P1H7QF7W1D7W1D7W1B3W1D7W1D7W1T77WNudh4oILLtAFF1yQmr/22mu1YsUKLVy4MLXMsiwZhiHTNGUYRp/lA9HS0i7TtD56xwvYkEEl+p8DRxQKtX3kfQUClWdkP+iL2tqH2tqH2tqDutqH2tqH2trH6dq6XMZJv6DP+RW5b7zxhn75y1+m5i3L0qhRoxQKhVLLQqGQqqurFQwGsy7HwAT9jOgEAACAMy/nYaKtrU2rVq1SOBxWe3u7Nm3apIcffli//OUv1draqs7OTv3oRz9STU2NJkyYoH379mn//v2KxWLaunWrampqct3lghes8nHjOgAAAJxxOT/N6fOf/7zefPNN/e3f/q1M09QNN9ygiy66SHfddZfmzZunSCSia6+9Vp/+9KclSStXrtTtt9+ucDisSZMmaerUqbnucsEbMcyn9s6I2jsjqij3Ot0dAAAAFAnDsqyivqCAayak3/5Psx7b+Jbq516kvxg15CPty+lz9ooZtbUPtbUPtbUHdbUPtbUPtbWP07XNq2smkHsj/YnhYTnVCQAAAGeQI3fARm4NH1omt8vQG3sPq7TErarBZaoaUqbBPu+AR8cCAAAAkggTNrKi3TKPHpJr6EgZnhLH+uF2uXTemGF6690WvfVuS2p5icclfyJYJAPG8LT5YZWlcrkIGwAAAMiOMGGjyO9fVvhXz0mGS65ho+QOjJVr+Fi5A5+Qy//xnAaMf/q7iTrRFVHzsS61HO+KPyamW4516b2mNrWdiGRs4zIM+QeXpoJG1eAyjf34UJW4pOGDy+QfXCavhzPlAAAAzlaECRt5x/8vGYMDMkMNijU3KLr/t7L2vhJvNNxy+bMEDLd9oy35yrz6szKv/mxE9jsohiMxtSbCRXNa0Gg+3qU9+4/oaHtYvS/XH1JRknE0o/djeSlvMQAAgGLFJz0bGW6PvGMvksZeJCl+gz6ro1Wx0L5UwIjs2yX94efxDVxuufwfl3v4J+QKjI0HjWEfl+HOzY+p1OvWyKpBGlk1KGt7NGbK8Hr0zr6WjLDRcrxLDYfatGtvSLFeI2cNKvP0OY3KP7hMlT6vyks9GlTmla/Mo7ISN9dvAAAAFBjCRA4ZhiGjokquiirpExdLSgSM9mbFQg0yQ/sUa96vyJ9ek/7ws/hGLo9cVaPlHj42ETA+Idewj8lw5f5H53G7FKgaJLdpZm03LUvH2rsTp1F1JoJGWC3HutR0pFO/bziicCSWdVvDkHyJcFFe5klMe+Qr88hX5pWvNDntka/U26eN060AAAByjzDhMMMwZFQG5KoMSH/+15ISAaMtpFioIX4Uo7lBkf/5lbTnp/GN3B65qv5M7uTpUYGxcg39mAyX28FXEr/GYlhlqYZVlma9n4VlWeroiqrlWJc6uiI60RXViXA08RhRR1dUnV1RdSTmjzaHU+tEotkDTJLX40oEjZ6jHb0DiC8jpMTnyxNBxOtxycWREQAAgAEhTOQhwzBkDK6Wa3C1vOd8RpJkWaas44fjAaM5fhQj8s5ORX7/cnwjd4lcw5MBY6xcwz8RH0XKlT/f2BuGoYpy72ndhTsSjelEKmgkAkhXRCfC0VQISQaSE11RHWvv1qGWjlQY6c+tGd0uQ16PSyWJcOHxuOV1u1JhI7nc63GllntSbe6M5X3+uV3yJtdJ7MuTtr6bUbMAAEABIkwUCMNwyRgSlGtIUN6/uFRSImAca1IstC9+mlRzgyJ7X1Hk7R/HN/KUyF01puf0qKEfk1FZJaO0ouCuT/B63BpS4daQitIBb2talsLdsdTRkM5EAElOR2KmItG0fzFTkWgsc1nUjB8h6UjOx9LWNRWNfbS7rBuG5HW75HIZcrsMud3xgJGc9iSmXS5DbrchtyvR7jbkSZuOb+NKWy/Rnmrrte8sz9UzferncSWeK97WM+12GQX3/gIAAKeHMFHADMMlY+hIuYaOlPeTn5UkWaYp89ih1AXesdA+RfbsUGT3Sz0bekrkqhguo3K4XBVVMiqr5KpITg+X4Rsiw8ifIxoflcswVF4aP6VJfc++OiNMy1I0aqo7I5CYikZ7wkfvtt7hpaTEq7b2sGKmqZhpKRazek1bipqmzMR8JGqqq9tMa7MUiyXWN634emZPe++L4+3kMoyTB5i0cNQ7sKSm00OLYchIPLpchlyG0qYNGa7EvGHIMJRanno0pMrBZTrR0f2h27tciX302t6VCEgZ7UbPtJHYr2GkP/b052TLey8jhAEACg1hosgYLpfcw0bJPWyUvOdeLkmyzJjMo4dkHmuS1d4ss60l9Rg5/K4U7sjcicsjo8IfDxcVw+WqrJKrcriMiipFvGNkmSWOXACez1yGoRKvWyXe079uJRCoVCjUdgZ7lcmyekJFMqiYZpYgkgofZtbpaNq0mdqfmTXQ9A1EfcNRcjociSkWi/Zsm9iXZVkyLck0LZmWlXiMBzgrtSw+X+gMKSOwxAOMZCgzmBhGfN1k+IjPx5enzysRUIzEMqlnXmnbpM8rbd/Z2lyJJ0jus6TEo0hiYIWe9fs+R/L1pfe5b3/Stj3VOln6nREkk2HQpbSQGA9vyaNqPfXtWW5kbN8rVLpOvTwZSlMBOBkUM36eaT9XQ2nP11PvZHssZsq0LK7lApD3+ER4FjBcbrn9H5fb//Gs7VZ3p8z21kTAaJbV3iKzrVlme4vM93+n6ImjqXUPSPH/2fuGpQJG8oiGq6IndDh5x29kZxiGPG5DHrck+25n4qh48OgJF6ZpybIsDfNXKBRqSwsj8UBi9Q4opiVLfbdPtqf2n9q2Z5mVDDiJ5++7XD37y7qvLPsz44/ZnsOyJEuSLMmK/0dmYkHPcvVaN3tbsnZWqo5p84nHjPWS25uWorH40bX0fvSsn72fGe29+qP05+rT5/T+9mzfU6P0GicDr1I/w0KNm9lDSM+yjLDiMjLDZ2IHfQLlKYJZT6CML3QZWdo+bL9ZAqWZeDOk/w5ISk1baT/j5DoZ66fe51bGe9Syem2jXuv02n96kEyGv+TRUJfLkDu9Le3xVG0Zp5i6+u67z37SQm3Gz7rXRE+07qlz73mjZ+VT7qPv9kbG+r1/B9N/X7P9/mW0J7atqChTW1tX1r8dyXXSf+f7/H2RJUNZavwhdTzVz6U/7cn3BE4PYQIySsrl9o+S/KOytlvRblkdrTLbWjTIaNfxQwdTRzdiH/xR0Y4jkpU52pJRPrhX0Igf4TAqhstV4ZdKyovqVCrkB8OI/0/D3eutNXhQicInCLh2sPuI2pnUEzbTAmOvMNkTNk++3EpsGzOT8z0hNpYIM5Yyg00qHCr9qFrmh9709Xy+ErW1h/uETyut3+kfoFOBqXdQzfhwqIxwZqa1JeuTESgzwl9m8Ezfp5V4XfEPmycPlOlHYHqOmCUDUvyXtmedzCNwRq9tXRn7SA8ymesoGbKSR8gMqby8RB0d4dSR1dQR1rTpnjCaOIoaMU/SHj/CG02+R7LsE4XBlSVwGH0CmtFrPrO9dyDps32WYPfhodBQqcelJTd9RmV5+rGJMIEPZXhKUhd/Dw5UKjwq84ODZcZknTiacVQjeRpVrPV9We/9VopF++7YXSLDWyp5ko+lMjylkjfx6CmR4SntafOW9F3HWxrvnye5TuIxRzf6A1A4ThY281EhhbRCk+vaJoNpLEvYSAYzKfPb+cRE+kPmdDLQ9WpIP0qQzuq9r2zb9zqq9GGnImY7lXH48Aq1tHRkPUKVfMg8PTO+tOe5eo409qnXKeqYrb3POtn2mXa6b+/2mJl55KV3zaxexe9b4+RD9p9nMmxn3ya5TnyqxOOWr8wjszvLZ6k8wCcufGSGy91zM74sLMuU1Xk8ETRaZHW0yOrukhUNS9FuWZGwFA3H5yNhWSeOyoqE40dEEstkDvAXyHBL3pJTh5Nkm9sjub2S2yPD7ZXc3sSjJzHtkdwl8XaXR/J4Zbi8kscbv74k8Si3l8OkAIAMLsOQK3mKaZGrGlKetx94C13VkPK8/YKBMAHbGYZLhm+o5Bsqd/U5p7UPy4zFg0cycKQ99kx397T1WjcVWCKdsk4cS6yT2CYWkczsd+YeMFd6AMkMJpkBxdsrwHjUUjlI4a6Y5HJLLnf8JoSpf560eU+Wdnf8OdPnXaeYN1wEHwAA8JERJlAQDJc7fp1FSbkt+7csM34qViwiKxZJm05fFl+enM62XvryzPa07bs7ZaXtJ9l+3IzJMqOZV6vZ6aTBIxlcXIng0Te4yHD1WpaYN9L32bNu9oDkjl83kzHfs7/U8p4f0slfS5+2zPnOjnJFj3UOcD+91zMSR6cSgc7du3aJNrenJ/DlaWCzLCt+tK/P+7D3siztZlRWNCKZ8fYjlYPUHbbi4djTE5QNjzdxRC85nd5W0hOmuXYKAAoaYQJQ/OiJPCXxU6Ec6kPyPF7LMuNHShL/rOR0LJo2H+3bbkbTppNt0bT5bO3Z9h2TrN77jkmWGT+ik9rWjO8zW38tM9FHs8/F+U7odOqJk+HI3SuwJU+Z6xNCegWU9PCSNp+s/anD70nCQiIEnKEXqO6POkZSr6N1PUEjscxTknYkr1cw8XhluHuCSfxiW1c8wBpG4jF9Otlu9FmeXNc4advJtu2Z7mlPXY6ZVqpey4xe7al14susWDT++9tzdnrGVZ8DCarxi6LNnn+J30srbTprW2reSvxNSFue3Gfyb0PGur2WJ0J66ueZOu00eR1cSdp8Sfw9jwGzYlEp0pU4Gt8VPzKffIx2xY/QR+KPVqQrcbQ+0R7rTpz6Wy6jpExG8ss7b/zR8JYnvtArk1Hik+Eti7e5+DIAhAkg7xiGS3K74h+Q1OfjRsHp+VCSDBxmIuSYqeBiJZanQk/6+lk/kGVxirahQwfp6NETp7GftGmrp3+KZQa11HQqlEUzp9MDXSx9Ohr/kJZcNxHW4uskw1+iLRbteX7DlXGNT/qpdUbyA3mJL0u7J+spdn1Px+t1PVDyOqHe27vcGl7lU+iD1nhoiUakWHdqOhVkUtPd8aMaycCTOM0wdbQvmhZ8km1dbbJ6tyW3zYOgapf2AW9hZISN1KNlqu9Rtjzn8pw6bGSEkQ+bT27vTS03w25Z3cmvGDLGLc3en1MczbQGvH2vdaIRWdGeD/6pD/zRXh/8E4GgZ7206cS6Azpd1+1NBIJSGZ4yyeOVFU1cz9jdGd9ff943ntJE2CiTSsoVrahUxPKmziQwvPFgokQAiU/Hw0kymMhblncBMn7hs9UTjrM89gnWadPWSbbJ2NZM/G5mPJqSaclSr3aPV5b/c47W5FQIEwBslfqmOO1Gh7kOSOWBSrX78vPCtUJnuNzxDwnespw/t2XG0oJHVPH/+Q/wf+h9lqdvf/Jt++4/fVpKG4+lz1AtVu/21ENq7FUNGlSqjo5wr/bk/rKc1td7WXI+eQTFlTyK4op/m2z0+udypf2uZrb1Wd+VPIrTd90+6yenpcTPKnlNW3fPdWv9nu+OD9CRPp/Y30B0fPgq+cVbljhqUCbDWxr/fSurkFFRlbEsPqJhIiAklstT1jOdWlb6oR/eLcuMh5XuTlmRTqm7Mx5cupPTnfG27k4p0hkPIZFOmZGwzI6WtO36G0pKeo5+eLw9vw8ZH+qt1HDFqS8SUo+Z6/XdJm1ZcjtLkhKPGe35Gb67RlRL5aOd7kZWjoSJNWvW6MUXX5QkTZo0Sffcc4+WLl2qXbt2qbw8fk78l7/8ZU2ePFk7d+7UihUrFA6HddVVV+muu+5yossAgDyTvAbHiSBjt2GBSkXzdOSWfGNZVuIIWHdGWOkbSOJtg8o9iaCWZGSd7PdR0Y+wferoibcsezDweB25rsgwXKd1nWLvYXcty4r/DNJDSSqEdKWWp0JKd2f8Zyn1hNDUqYNG3+m0RyM11qwrcXZg2va918/YJr299/Omn8Z48tMjT/8Uy2xtWdq9pSr/s7Fqz9O/CTkPEzt37tSrr76qTZs2yTAM3XrrrXrppZe0e/duPfvss6qurk6t29XVpfr6ej3zzDMaOXKkFixYoB07dmjSpEm57jYAAMhDhmGkXfNW8aHrDw1UKpKnH8qKjWEYiZBUJmmY092BTXIedwOBgJYsWaKSkhJ5vV6dc845amxsVGNjo+rr61VXV6fHHntMpmnqrbfe0pgxYzR69Gh5PB7V1dVp27Ztue4yAAAAgCxyfmTik5/8ZGq6oaFBL774otatW6fXXntNy5YtU2VlpRYsWKCNGzfK5/MpEAik1q+urlZTU9OAnq+q6sO/pcDABAKVTnehaFFb+1Bb+1Bbe1BX+1Bb+1Bb++RrbR27APudd97RggULdM899+jP//zP9fjjj6fa5s6dq82bN6u2tjZj+DvLsgY8bntLS7tMMz8vpilEvc+HxJlDbe1Dbe1Dbe1BXe1Dbe1Dbe3jdG1dLuOkX9A7MkDwrl27dNNNN+nuu+/WNddco71792r79u2pdsuy5PF4FAwGFQqFUstDoVDGNRUAAAAAnJPzMHHo0CEtWrRIq1ev1rRp0yTFw8NDDz2kY8eOKRKJaP369Zo8ebImTJigffv2af/+/YrFYtq6datqampy3WUAAAAAWeT8NKennnpK4XBYK1euTC2bNWuWbrvtNs2ePVvRaFRTpkzR9OnTJUkrV67U7bffrnA4rEmTJmnq1Km57jIAAACALAzLOsWtHYsA10ycWU6fs1fMqK19qK19qK09qKt9qK19qK19nK5t3l0zAQAAAKDwESYAAAAAnBbHhobNFZdrYEPJ4sNRU/tQW/tQW/tQW3tQV/tQW/tQW/s4WdtTPXfRXzMBAAAAwB6c5gQAAADgtBAmAAAAAJwWwgQAAACA00KYAAAAAHBaCBMAAAAATgthAgAAAMBpIUwAAAAAOC2ECQAAAACnhTABAAAA4LQQJgAAAACcFsIEslqzZo2mTZumadOmadWqVVnbP//5z2vGjBmaMWOG1q1b50AvC9PcuXM1bdq0VO3efPPNjPadO3eqrq5OU6ZM0aOPPupQLwvP97///VRNZ8yYoYsuukjLly/PWIf37cC0t7dr+vTpev/99yX1773Z2NioOXPmaOrUqVq4cKE6Ojpy2eWC0bu269ev1/Tp01VXV6elS5equ7u7zzabNm3S5z73udT7l78P2fWu7dKlSzVlypRU3V566aU+2+zZs0czZ85UbW2t7r33XkWj0Vx3O++l13XHjh0Zf28vvfRSLViwoM82vGc/XLbPWwX3t9YCevnFL35h/d3f/Z0VDoet7u5ua968edaPfvSjjHUWLFhg/frXv3aoh4XLNE3rc5/7nBWJRLK2d3Z2WpMmTbLee+89KxKJWPPnz7d+9rOf5biXhe+Pf/yjNXnyZKulpSVjOe/b/vvtb39rTZ8+3Ro/frx14MCBfr83b7vtNmvr1q2WZVnWmjVrrFWrVuW663mvd23/9Kc/WZMnT7ba2tos0zSte+65x3r66af7bLd8+XLrhz/8Ye47XEB619ayLGv69OlWU1PTKbebNm2a9Zvf/MayLMtaunSptW7dOtv7Wkiy1TXp8OHD1he+8AVr3759fbbjPXtq2T5v/fCHPyy4v7UcmUAfgUBAS5YsUUlJibxer8455xw1NjZmrLN792498cQTqqur0/LlyxUOhx3qbWH505/+JEmaP3++vvjFL+rZZ5/NaH/rrbc0ZswYjR49Wh6PR3V1ddq2bZsTXS1oX//613XXXXfJ7/dnLOd9238bNmzQsmXLVF1dLal/781IJKLXX39dtbW1kqSZM2fy/s2id21LSkq0bNkyVVRUyDAMnXvuuX3+5krS7373O23atEl1dXX66le/qmPHjuW663mvd207OzvV2Nio+vp61dXV6bHHHpNpmhnbHDx4UF1dXZo4caIk3rfZX8JihAAAB/9JREFU9K5rulWrVmnWrFkaO3Zsnzbes6eW7fNWQ0NDwf2tJUygj09+8pOpP6oNDQ168cUXNWnSpFR7R0eHzjvvPC1evFibNm3S8ePHtXbtWqe6W1COHz+uyy67TI8//ri+853v6LnnntMvfvGLVPvhw4cVCARS89XV1WpqanKiqwVr586d6urq0lVXXZWxnPftwDz44IO6+OKLU/P9eW8eOXJEFRUV8ng8kuL/o+T921fv2o4aNUqXX365JKm1tVXr1q3TF77whT7bBQIB/eM//qOef/55jRw5ss9pfOhb2+bmZl166aV66KGHtGHDBr3xxhvauHFjxja939u8b/vqXdekhoYGvfbaa5o3b17W7XjPnlq2z1uGYRTc31rCBE7qnXfe0fz583XPPfdkfOMwaNAgPfnkkzrnnHPk8Xg0f/587dixw7mOFpALLrhAq1atUmVlpfx+v6699tqM2pmmKcMwUvOWZWXM48M999xzuvnmm/ss53370fTnvZltGe/f/mtqatKNN96oL33pS7rkkkv6tD/++OO66KKLZBiGbr31Vr3yyisO9LKwjB49Wo8//riqq6tVXl6uuXPn9vm95+/u6Vu/fr1uuOEGlZSUZG3nPds/6Z+3Ro8eXXB/awkTyGrXrl266aabdPfdd+uaa67JaGtsbMz4ZseyrFQ6xqm98cYb+uUvf5ma7127YDCoUCiUmg+FQlkPKyO77u5uvf7667ryyiv7tPG+/Wj68970+/1qa2tTLBY76TrI7t1339WsWbN0zTXXaNGiRX3a29ra9J3vfCc1b1mW3G53DntYmPbu3avt27en5rP93vd+bzc3N/O+7aef/OQnuvrqq7O28Z7tn96ftwrxby1hAn0cOnRIixYt0urVqzVt2rQ+7WVlZXr44Yd14MABWZaldevWafLkyQ70tPC0tbVp1apVCofDam9v16ZNmzJqN2HCBO3bt0/79+9XLBbT1q1bVVNT42CPC8vevXs1duxY+Xy+Pm28bz+a/rw3vV6vLr74Yr3wwguSpM2bN/P+7Yf29nbdcsstuvPOOzV//vys6/h8Pv37v/97avS3Z599lvdvP1iWpYceekjHjh1TJBLR+vXr+9Rt1KhRKi0t1a5duyRJW7Zs4X3bD62trerq6tLo0aOztvOe/XDZPm8V4t9avpZDH0899ZTC4bBWrlyZWjZr1iy9/PLLuuOOO/SpT31Ky5cv///buZfQppYADuNf+ojdKNKNDyKIio0oQqUbq6FaRVwUKor4QroQCxJEqdgGKmhtLFkEEXyBC0tBEYuPRbp1I4h0GayQKgiiRUQhQouPJiZ3IfZSrBfMrbe3+P0gi0zCOTPDcDj/M3OGI0eOkMvlWLdu3ZTLSvSjzZs3k06n2bFjB4VCgf3791NbW0tzczPXrl1jwYIFJBIJjh49ypcvX2hoaGD79u0zXe1Z49WrVyxcuHBS2eHDhx2302DOnDk/HZudnZ00NjayZcsWTp8+TSwW4+rVqyxatIjz58/PcM3//+7cucP79+/p7e2lt7cXgMbGRo4dOzapby9cuMCZM2f4/PkzS5cunXLbbk0WDodpbW1l37595PN5tm3bRlNTEzD52pBMJjl16hRjY2OsXr36p+8A6G+vX7/+4XoLOGZ/wc/ut2bbtTZQLBaLM3Z2SZIkSbOWy5wkSZIklcQwIUmSJKkkhglJkiRJJTFMSJIkSSqJYUKSJElSSdwaVpL0y2pqali5ciVlZZOfSV2+fJlQKDTt53r8+DHV1dXTelxJ0r9nmJAklaSvr88bfEn6wxkmJEnTanBwkGQyyeLFi3nx4gVVVVUkEgmWL1/O6OgoXV1dZDIZAoEAkUiEtrY2KioqSKfTxONxPn36RGVlJe3t7axfvx6Aixcvkk6n+fDhA4cOHeLAgQO8e/eOjo4OstksAA0NDRw/fnwmmy5JfxzfmZAklaSlpYXm5uaJTzQanfhtaGiIgwcPkkql2LlzJydPngQgHo8zf/58UqkUd+/eZXh4mOvXr5PL5YhGo0SjUQYGBuju7qanp4dCoQDAkiVLuHfvHpcuXSKRSJDL5ejv7ycUCnH//n1u3rzJy5cvGR0dnZG+kKQ/lTMTkqSS/NMyp3A4TF1dHQC7du3i7NmzZLNZHj58yK1btwgEAgSDQfbu3UtfXx8bNmygrKyMTZs2AbBmzRpSqdTE8ZqamgBYtWoV4+PjjI2NEYlEaG1t5c2bN9TX13PixAnmzp37exstSZrEmQlJ0rQrLy+fsqxQKBAIBCbKCoUC+Xye8vLySeUAz549I5/PA1BR8e3Z1/f/FItF1q5dy4MHD9izZw8jIyPs3r2boaGh39UkSdIUDBOSpGmXyWTIZDIA3L59m9raWubNm8fGjRu5ceMGxWKR8fFx+vv7qa+vZ9myZQQCAR49egTA06dPaWlpmVjmNJVkMsmVK1fYunUrnZ2drFixgufPn/8n7ZMkfRMoFovFma6EJGl2+dnWsG1tbVRVVdHR0UE4HGZkZITq6mrOnTtHKBQim80Sj8cZHh4ml8sRiURob28nGAzy5MkTenp6+PjxI5WVlcRiMerq6n7YGvb7969fvxKLxXj79i3BYJCamhq6uroIBoMz0SWS9EcyTEiSptXg4CDd3d0MDAzMdFUkSb+Zy5wkSZIklcSZCUmSJEklcWZCkiRJUkkME5IkSZJKYpiQJEmSVBLDhCRJkqSSGCYkSZIkleQvpyq2T59K9S0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 936x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_loss(fit_history):\n",
    "    plt.figure(figsize=(13,5))\n",
    "    plt.plot(range(1, len(fit_history.history['loss'])+1), fit_history.history['loss'], label='train')\n",
    "    plt.plot(range(1, len(fit_history.history['val_loss'])+1), fit_history.history['val_loss'], label='validate')\n",
    "    plt.xlabel('Epochs')\n",
    "    plt.ylabel('Loss')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "\n",
    "plot_loss(history)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evaluate(y_true, y_hat, label='test'):\n",
    "    mse = mean_squared_error(y_true, y_hat)\n",
    "    rmse = np.sqrt(mse)\n",
    "    variance = r2_score(y_true, y_hat)\n",
    "    print('{} set RMSE:{}, R2:{}'.format(label, rmse, variance))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train set RMSE:18.004364211864736, R2:0.813339453755806\n"
     ]
    }
   ],
   "source": [
    "# predict and evaluate\n",
    "y_hat_train = model.predict(X_train_scaled[train_cols])\n",
    "evaluate(y_train_clipped, y_hat_train, 'train')\n",
    "\n",
    "# y_hat_test = model.predict(X_test_scaled[train_cols])\n",
    "# evaluate(y_test, y_hat_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# add smoothing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exponential_smoothing(df, sensors, n_samples, alpha=0.4):\n",
    "    df = df.copy()\n",
    "    # first, calculate the exponential weighted mean of desired sensors\n",
    "    df[sensors] = df.groupby('unit_nr')[sensors].apply(lambda x: x.ewm(alpha=alpha).mean())\n",
    "    \n",
    "    # second, drop first n_samples of each unit_nr to reduce filter delay\n",
    "    def create_mask(data, samples):\n",
    "        result = np.ones_like(data)\n",
    "        result[0:samples] = 0\n",
    "        return result\n",
    "    \n",
    "    mask = df.groupby('unit_nr')['unit_nr'].transform(create_mask, samples=n_samples).astype(bool)\n",
    "    df = df[mask]\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3:\n",
    "You can run the cells below with different values of alpha. When you decide to go back and use the initial value of alpha = 0.4, you get the same RMSE from the NN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prep data\n",
    "alpha = 0.4\n",
    "\n",
    "X_train_interim = train[index_names + sensor_names].copy()\n",
    "X_test_interim = test[index_names + sensor_names].copy()\n",
    "\n",
    "X_train_interim[sensor_names] = pd.DataFrame(scaler.transform(X_train_interim[sensor_names]), columns=sensor_names)\n",
    "X_test_interim[sensor_names] = pd.DataFrame(scaler.transform(X_test_interim[sensor_names]), columns=sensor_names)\n",
    "\n",
    "X_train_interim = exponential_smoothing(X_train_interim, sensor_names, 0, alpha)\n",
    "X_test_interim = exponential_smoothing(X_test_interim, sensor_names, 0, alpha)\n",
    "\n",
    "X_train_smooth = X_train_interim.drop(index_names, axis=1)\n",
    "X_test_smooth = X_test_interim.drop('time_cycles', axis=1).groupby('unit_nr').last().copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "split_result = train_val_group_split(X_train_smooth, y_train_clipped, gss, train['unit_nr'], print_groups=False)\n",
    "X_train_split_smooth, y_train_clipped_split_smooth, X_val_split_smooth, y_val_clipped_split_smooth = split_result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train on 16561 samples, validate on 4070 samples\n",
      "Epoch 1/20\n",
      "16561/16561 [==============================] - 2s 104us/sample - loss: 2329.9827 - val_loss: 398.0225\n",
      "Epoch 2/20\n",
      "16561/16561 [==============================] - 1s 61us/sample - loss: 406.2162 - val_loss: 337.6243\n",
      "Epoch 3/20\n",
      "16561/16561 [==============================] - 1s 65us/sample - loss: 380.8759 - val_loss: 315.8910\n",
      "Epoch 4/20\n",
      "16561/16561 [==============================] - 1s 63us/sample - loss: 372.8627 - val_loss: 311.3165\n",
      "Epoch 5/20\n",
      "16561/16561 [==============================] - 1s 62us/sample - loss: 370.6699 - val_loss: 308.2596\n",
      "Epoch 6/20\n",
      "16561/16561 [==============================] - 1s 67us/sample - loss: 368.6668 - val_loss: 302.1523\n",
      "Epoch 7/20\n",
      "16561/16561 [==============================] - 1s 68us/sample - loss: 364.0241 - val_loss: 299.3743\n",
      "Epoch 8/20\n",
      "16561/16561 [==============================] - 1s 61us/sample - loss: 359.5961 - val_loss: 295.8907\n",
      "Epoch 9/20\n",
      "16561/16561 [==============================] - 1s 61us/sample - loss: 354.4997 - val_loss: 287.9995\n",
      "Epoch 10/20\n",
      "16561/16561 [==============================] - 1s 59us/sample - loss: 348.8014 - val_loss: 286.1727\n",
      "Epoch 11/20\n",
      "16561/16561 [==============================] - 1s 67us/sample - loss: 344.7535 - val_loss: 279.6822\n",
      "Epoch 12/20\n",
      "16561/16561 [==============================] - 1s 59us/sample - loss: 339.1559 - val_loss: 276.1613\n",
      "Epoch 13/20\n",
      "16561/16561 [==============================] - 1s 62us/sample - loss: 333.1796 - val_loss: 271.6287\n",
      "Epoch 14/20\n",
      "16561/16561 [==============================] - 1s 58us/sample - loss: 329.0121 - val_loss: 276.0835\n",
      "Epoch 15/20\n",
      "16561/16561 [==============================] - 1s 58us/sample - loss: 328.3134 - val_loss: 261.8976\n",
      "Epoch 16/20\n",
      "16561/16561 [==============================] - 1s 60us/sample - loss: 325.2139 - val_loss: 271.7145\n",
      "Epoch 17/20\n",
      "16561/16561 [==============================] - 1s 64us/sample - loss: 323.7677 - val_loss: 274.5235\n",
      "Epoch 18/20\n",
      "16561/16561 [==============================] - 1s 57us/sample - loss: 323.9837 - val_loss: 258.4319\n",
      "Epoch 19/20\n",
      "16561/16561 [==============================] - 1s 58us/sample - loss: 319.4420 - val_loss: 257.8707\n",
      "Epoch 20/20\n",
      "16561/16561 [==============================] - 1s 59us/sample - loss: 319.0945 - val_loss: 256.8070\n"
     ]
    }
   ],
   "source": [
    "# reset initial optimizer parameters by recompiling model\n",
    "# and reload initial model weights to allow repeated running of cell with consistent results\n",
    "model.compile(loss='mean_squared_error', optimizer='adam')\n",
    "model.load_weights('reproducible_primer_weights.h5')\n",
    "\n",
    "history = model.fit(X_train_split_smooth[train_cols], y_train_clipped_split_smooth,\n",
    "                    validation_data=(X_val_split_smooth[train_cols], y_val_clipped_split_smooth),\n",
    "                    epochs=epochs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxMAAAE/CAYAAAA9q540AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de3wU5aH/8e/sLckmAbJhk0Xk560qVquinHoHq4eLhUil1iIWtNrqz+MLW2u1iBcqKvJDWntaxV9rW+2pnlZqi1Z+FGu1xQvWKlosVtFTAREwJBDIhWyyuzO/Pza72VswxMzO7vJ5v15xZ555ZvbZJ2vY7z4zzxiWZVkCAAAAgP3kcroBAAAAAIoTYQIAAADAgBAmAAAAAAwIYQIAAADAgBAmAAAAAAwIYQIAAADAgBAmAAAAAAyIx+kG2K2lpUOmya00BkttbZV27mx3uhklib61D31rH/rWHvSrfehb+9C39nG6b10uQzU1lTm3lXyYME2LMDHI6E/70Lf2oW/tQ9/ag361D31rH/rWPoXat5zmBAAAAGBACBMAAAAABqTkT3MCAADAgceyLLW371FnZ7tMM+Z0cz6RHTtcMk3T9ufxeHyqqQnK7e5/RCBMAAAAoOS0tDTJMAwFAvVyuz0yDMPpJg2Yx+NSNGpvmLAsSx0drWppadLw4SP6vR+nOQEAAKDkdHeHNWxYrTweb1EHiXwxDEOVlUMUjXbv136ECQAAAJQgS4bBR939MZDQRQ8DAAAAGBDCBAAAAGCz9vZ23XTTt/td/513/qlFi+6wsUWDgzBho9fe2aGfrfin080AAACAw9raWvXeexv6XX/06E9r7txbbWzR4GA2Jxtta+7QmvUfadako+Xzup1uDgAAABzygx/co+bmJt1007e1efNGDR06TGVlZbrrrsW6++471NS0Q83NTRo79rOaO/dWvfHGWv385z/Rfff9RFdf/XUdc8yntW7d37V7d4u++c0bdNppZzj9kiQRJmxVH/DLktTY0qlRdVVONwcAAOCA9dI/tuvFN7fbcuwzjx+hMz6z7+lUv/nNGzRnzlW69tpv6UtfOl+/+c2PNGLEQXrmmVU68sijdOed/0eRSERf+cqXtGHDO1n7RyJR/fjHD+nFF5/Xgw8+QJg4EIQCfknSR7v2EiYAAAAgSaqpCWjEiIMkSRMmTNY//7ley5b9tzZt2qg9e/aos3Nv1j6nnHKaJOnww49QW1trXtu7L4QJGyXDxM4Oh1sCAABwYDvjMx8/epAvZWVlyeXHH/+1/vKX53T++Rfowgs/q40b/yXLsrL28fl8kuLTt+ba7hQuwLZRmc+tmuoyfbSr0+mmAAAAwEFut1uxWCyr/NVXX9H550/XxInnqbu7W++9965M0967XQ8mRiZsFgr49dGu7KEqAAAAHDgCgVrV14e0cOHtaeUXXTRTS5bcrUceeUiVlVU67rjjtX37No0cebBDLd0/hAmbhQJ+/fWfjbIsi1u5AwAAHKA8Ho/+7//9eVb5ySf/m371q9/l3Oekk8ZKkh544EFFo/HRihEjDtLjjz9lX0P3E6c52SwU8KuzK6rWvRGnmwIAAAAMKsKEzUK18YuwGznVCQAAACWGMGGz1OlhAQAAgFJCmLBZ7ZByedwufbSTMAEAAIDSQpiwmctlqL6mgpEJAAAAlBzCRB4wPSwAAABKEWEiD0K1fjXt7lQ0Vjw3IAEAAAA+DmEiD+pr/IqZlpr3hJ1uCgAAAArcXXd9VytXPqXm5iZ9+9vX5qxz5plj93mMbdu26u67F9jRvDSEiTxITA/LRdgAAADor+HDg1qy5IcD2vejj7Zr69YPB7lF2bgDdh4wPSwAAICzIu++pMiG5205tvfocfIedcY+68ybd4MmTpyss88+V5J0+eVf0Zw51+knP1mqrq6w2trade211+mss85O7rN9+zbNmXOVnnji/2n79m1asOBWdXZ26thjj0vWaWraobvvvkPt7W1qbm7S5z/foK997X/rP/9zibZt26rvfe//6Prrv6Nf/vJh/fnPzygWM3XKKafq6quvlWEYn/i1MzKRB1UVXlVVeAkTAAAAB6hJkz6vP/3paUnSli0fqLu7W7/97WOaO/dW/fznj2ru3Fv04IMP9Ln/vfcu1uc/36CHH/5vfeYzJyTLn3nmaU2YMEk/+cnD+q//ekzLlv1Ku3fv1je+8W0dffQxuv767+ivf12jDRve1oMP/pceeuhRNTU16Y9//MOgvC5GJvKEGZ0AAACc4z3qjI8dPbDT6aefqXvvXay9ezv0pz89rUmTztNFF83UmjUv6M9//pPeeusf6uzs7HP/N95Yq+9+9y5J0sSJ52nRojskSTNnztLrr7+m//7vX2rjxn8pGo0oHE4/zmuv/U3//Od6XXHFLElSV1dY9fWhQXldhIk8CQX8evP9nU43AwAAAA7wer0644yz9OKLz+u5557RPff8p6655us66aSTNWbMyTr55H/T7bffso8jGDJNK75kGHK53JKkH/3oXm3btlUTJkzWuHFn67XX/ibLstL2NM2YLrroYs2Y8RVJUltbm9xu96C8Lk5zypNQrV+tHd3aG4463RQAAAA4YNKkz+vXv35EQ4cOk9/v15Ytm3XFFf9bp556hl54YbVMs+/bCIwd+1k9/fRKSdLq1c+pu7tLkvTaa69o5sxZOuecf9cHH2xWU9MOmaYpt9ujWCwmSTrppH/T00+v1N69exWNRnXTTdfrL395dlBeEyMTeZK4CLuxZa8OGzHE4dYAAAAg344//kS1t7frC1+4UEOGDNXUqdM0a9ZF8ng8Oumkf1M4HO7zVKdvfetG3XHHbfr975dr9Ohj5PdXSpK+8pXLdMcdt6msrEx1dSGNHv1pbdu2VUcddbTa29t0xx236tZb79D//M+7uvLKy2SaMZ1yyuk677ypg/KaDCtzHKTE7NzZnhwSctLW5g7d+tNX9PWpn9Zpxw3OOWpOCAar1dTU5nQzShJ9ax/61j70rT3oV/vQt/YptL796KPNCoUOcboZg8LjcSkazc/Nj3P1m8tlqLa2Kmd9TnPKk7phFTIMaTsXYQMAAKBEECbyxOtxKTi0Qo2ECQAAAJQIwkQehWqZHhYAACA/DFlWfk4NKhUDufqBMJFH9TV+Ne7aK7O0L1MBAABwnM9Xrt27mxWNRgb0IflAY1mWOjpa5fH49ms/ZnPKo1CtX91RUy2tXaodWu50cwAAAEpWTU1Q7e17tGtXo0wz5nRzPhGXy7XPaWMHi8fjU01NcP/2saktyCExPexHLXsJEwAAADYyDEPV1cNUXT3M6aZ8YoU2U1YqTnPKo2SY2Ml1EwAAACh+hIk8GlblU5nPzUXYAAAAKAmEiTwyDEOhGmZ0AgAAQGkgTORZqNbPvSYAAABQEggTeRYK+LVzT1jdkeKeVQAAAAAgTORZfaBClqQdLZ1ONwUAAAD4RAgTeTYiUClJXDcBAACAomdrmLjvvvs0ZcoUTZkyRYsXL5YkrVmzRg0NDZo4caLuvffeZN23335b06dP16RJk3TzzTcrGo1KkrZt26ZLLrlEkydP1tVXX62Ojg47m2y7+kCFJMIEAAAAip9tYWLNmjV68cUXtXz5cj3xxBN66623tGLFCs2bN09Lly7VypUrtX79eq1evVqSdMMNN+i2227T008/LcuytGzZMknS7bffrpkzZ2rVqlU67rjjtHTpUruanBflPo9qqssIEwAAACh6toWJYDCouXPnyufzyev16ogjjtCmTZt0yCGHaNSoUfJ4PGpoaNCqVau0detWhcNhnXjiiZKk6dOna9WqVYpEInr11Vc1adKktPJiV19TQZgAAABA0fPYdeAjjzwyubxp0yb94Q9/0Fe+8hUFg8FkeV1dnRobG7Vjx4608mAwqMbGRrW0tKiqqkoejyetfH/U1lZ9wlcy+A4bOUzP/32rhg+vkmEYTjdnvwWD1U43oWTRt/ahb+1D39qDfrUPfWsf+tY+hdq3toWJhPfee09XXXWVbrzxRrndbm3atCm5zbIsGYYh0zTTPlQnyhOPqfb3w/fOne0yTesTvYbBNtTvVUdnRO9/sEtD/D6nm7NfgsFqNTW1Od2MkkTf2oe+tQ99aw/61T70rX3oW/s43bcul9HnF/S2XoC9du1aXXbZZbr++ut1wQUXKBQKqampKbm9qalJdXV1WeXNzc2qq6tTIBBQW1ubYrFYWv1iFwr4JUkf7eRUJwAAABQv28LE9u3bdc0112jJkiWaMmWKJOmEE07Qxo0btXnzZsViMa1YsULjxo3TyJEjVVZWprVr10qSnnzySY0bN05er1djx47VypUrJUlPPPGExo0bZ1eT8ybEjE4AAAAoAbad5vSzn/1MXV1dWrRoUbJsxowZWrRokebMmaOuri6NHz9ekydPliQtWbJEt9xyi9rb23Xsscdq9uzZkqT58+dr7ty5euCBBzRixAh9//vft6vJeTN8aIU8boMwAQAAgKJmWJZVWBcUDLJCvGZCkm756Suqr6nQnC8e73RT9ovT5+yVMvrWPvStfehbe9Cv9qFv7UPf2sfpvnXsmgn0LRTwMzIBAACAokaYcEh9oEI7WjoVM02nmwIAAAAMCGHCIaGAXzHTUvPusNNNAQAAAAaEMOGQEYFKSczoBAAAgOJFmHBIqLbnXhOECQAAABQpwoRDqiq8qqrwEiYAAABQtAgTDqoPVHAXbAAAABQtwoSDQgG/PmohTAAAAKA4ESYcFAr4tae9W51dUaebAgAAAOw3woSDQszoBAAAgCJGmHBQKFAhiTABAACA4kSYcFBdjV+GITUSJgAAAFCECBMO8npcGj60nJEJAAAAFCXChMNCgUqmhwUAAEBRIkw4rD5QoY9a9sq0LKebAgAAAOwXwoTDRgT86o6Y2t3W5XRTAAAAgP1CmHBYKOCXxIxOAAAAKD6ECYeFarnXBAAAAIoTYcJhw6p8KvO6uQgbAAAARYcw4TDDMJIXYQMAAADFhDBRAEIBPyMTAAAAKDqEiQIQCvi1c09YkWjM6aYAAAAA/UaYKAChgF+WpMaWTqebAgAAAPQbYaIAhGp7poflVCcAAAAUEcJEAaiviYeJRi7CBgAAQBEhTBSAijKPhlX5GJkAAABAUSFMFIhQwM+N6wAAAFBUCBMFIhEmLMtyuikAAABAvxAmCkQo4FdHOKr2zojTTQEAAAD6hTBRIJIzOnGqEwAAAIoEYaJA1AeYHhYAAADFhTBRIIYPLZfbZTAyAQAAgKJBmCgQbpdLdTUVhAkAAAAUDcJEAWF6WAAAABQTwkQBCQX82tHSqZhpOt0UAAAA4GMRJgpIKOBXzLTUvCfsdFMAAACAj0WYKCCJ6WEbOdUJAAAARYAwUUBCTA8LAACAIkKYKCBVFV5Vlnu4CBsAAABFgTBRQAzDYEYnAAAAFA3CRIEhTAAAAKBYECYKTKjWr93t3ersijrdFAAAAGCfCBMFpr6mZ0anFkYnAAAAUNgIEwUmMT0sMzoBAACg0BEmCkx9TYUMiesmAAAAUPAIEwXG63Grdmg5YQIAAAAFjzBRgJjRCQAAAMWAMFGAQgG/Gnd1yrIsp5sCAAAA9IkwUYBCtX51RWLa3d7tdFMAAACAPhEmClAokJjRqcPhlgAAAAB9sz1MtLe3a+rUqfrwww8lSTfddJMmTpyoadOmadq0aXrmmWckSWvWrFFDQ4MmTpyoe++9N7n/22+/renTp2vSpEm6+eabFY2W/s3ckmGC6yYAAABQwGwNE+vWrdPFF1+sTZs2JcvWr1+vRx55RE8++aSefPJJTZgwQeFwWPPmzdPSpUu1cuVKrV+/XqtXr5Yk3XDDDbrtttv09NNPy7IsLVu2zM4mF4Rh1WXyeV3aTpgAAABAAbM1TCxbtkzz589XXV2dJKmzs1Pbtm3TvHnz1NDQoB/+8IcyTVNvvvmmDjnkEI0aNUoej0cNDQ1atWqVtm7dqnA4rBNPPFGSNH36dK1atcrOJhcEl2EoVBO/CBsAAAAoVB47D37XXXelrTc3N+vUU0/V/PnzVV1drauuukqPP/64/H6/gsFgsl5dXZ0aGxu1Y8eOtPJgMKjGxkY7m1wwQrV+bdze6nQzAAAAgD7ZGiYyjRo1Svfff39yfdasWXriiSc0adIkGYaRLLcsS4ZhyDTNnOX7o7a26pM33AGHH1yj197ZoWE1fnk9bqebkyYYrHa6CSWLvrUPfWsf+tYe9Kt96Fv70Lf2KdS+zWuY2LBhgzZt2qRJkyZJiocDj8ejUCikpqamZL2mpibV1dVllTc3NydPmeqvnTvbZZrFd7+G6nK3TEt6690dGhksnEAUDFarqanN6WaUJPrWPvStfehbe9Cv9qFv7UPf2sfpvnW5jD6/oM/r1LCWZWnhwoXas2ePIpGIHnvsMU2YMEEnnHCCNm7cqM2bNysWi2nFihUaN26cRo4cqbKyMq1du1aS9OSTT2rcuHH5bLJjemd04roJAAAAFKa8jkyMHj1aV155pS6++GJFo1FNnDhRU6dOlSQtWrRIc+bMUVdXl8aPH6/JkydLkpYsWaJbbrlF7e3tOvbYYzV79ux8NtkxvWGiQ1Jw35UBAAAABxiWZRXfOUD7oVhPc5Kk6+57UccdFtAVUz7tdFOSnB5mK2X0rX3oW/vQt/agX+1D39qHvrWP031bMKc5Yf+EavzcuA4AAAAFizBRwEK13GsCAAAAhYswUcBCAb/aOyNq74w43RQAAAAgC2GigCUvwt7JqU4AAAAoPISJApYIE9t3dTjcEgAAACAbYaKADR9WLrfL4LoJAAAAFCTCRAFzu1yqq6lgRicAAAAUJMJEgQsFmB4WAAAAhYkwUeDqA37taNlbtDfeAwAAQOkiTBS4UMCvaMxSc2vY6aYAAAAAafoVJpqbm/Xss89Kku655x5deumleuedd2xtGOKYHhYAAACFql9hYu7cudqyZYtefvllvfDCC5o2bZruvPNOu9sGxe+CLYnrJgAAAFBw+hUmdu/ercsuu0zPP/+8pk6dqunTp6uzk+lK86G6wit/mYcwAQAAgILTrzARiUQUiUT0wgsv6PTTT1dnZ6f27uXDbT4YhqFQrV+NhAkAAAAUmH6FiXPPPVennXaaampqdNxxx+lLX/qSpk6danfb0IPpYQEAAFCIPP2pdO211+qiiy5SfX29JGnJkiUaPXq0rQ1Dr1DArzXrP1K4O6pyX79+ZQAAAIDt+j2b01tvvSXDMHTPPffo7rvvZjanPErM6NS4i+tUAAAAUDiYzakIJKeH5VQnAAAAFBBmcyoCdTUVMkSYAAAAQGFhNqci4PO6VTu0nDABAACAgsJsTkWiPuDnLtgAAAAoKPs1m1MoFJLEbE5OCAX8enHrdlmWJcMwnG4OAAAA0L8wYZqmnnrqKT3//POKRqM644wz9KlPfUoeD9OU5kso4FdXd0y727tVU13mdHMAAACA/p3m9L3vfU9//etfdemll+qrX/2q3njjDS1evNjutiFFqJYZnQAAAFBY+jW08MILL+i3v/2tvF6vJOnss8/W+eefr3nz5tnaOPQK1fSGiWMOqXG4NQAAAEA/RyYsy0oGCUny+Xxp67BfzZAy+TwuLsIGAABAwehXmBg9erQWLlyoDz74QFu2bNHdd9+to446yu62IYXLMFQf8KuxhTABAACAwtCvMDF//ny1trZqxowZuuiii7Rz505dfPHFdrcNGUJMDwsAAIAC0q9rJqqqqrRo0aK0spNOOkmvv/66LY1CbvUBv17bsEORqCmvp185EAAAALDNgD+RWpY1mO1AP4wI+GVZ0o7dnU43BQAAABh4mODGafmXmB62kelhAQAAUAA4V6aIhALcawIAAACFY5/XTIwZMybnCIRlWQqHw7Y1CrlVlHk0tNLHRdgAAAAoCPsMEytWrMhXO9BP9QE/IxMAAAAoCPsMEyNHjsxXO9BPoYBfr7/b5HQzAAAAAK6ZKDahgF/tnRG1d0acbgoAAAAOcISJIpOY0YlTnQAAAOA0wkSRSc7oxEXYAAAAcBhhosgMH1out8tQYwthAgAAAM4iTBQZj9ul4LAKRiYAAADgOMJEEQoxPSwAAAAKAGGiCIUCfjW2dMo0LaebAgAAgAMYYaIIhWr9isZM7WzlLuQAAABwDmGiCCVndOJUJwAAADiIMFGEmB4WAAAAhYAwUYSq/V5VlHkYmQAAAICjCBNFyDAMZnQCAACA4wgTRYowAQAAAKcRJopUqNavlrYudXXHnG4KAAAADlCEiSLFjE4AAABwGmGiSCXCRGMLYQIAAADOsDVMtLe3a+rUqfrwww8lSWvWrFFDQ4MmTpyoe++9N1nv7bff1vTp0zVp0iTdfPPNikajkqRt27bpkksu0eTJk3X11Vero6PDzuYWlfqaChlielgAAAA4x7YwsW7dOl188cXatGmTJCkcDmvevHlaunSpVq5cqfXr12v16tWSpBtuuEG33Xabnn76aVmWpWXLlkmSbr/9ds2cOVOrVq3Scccdp6VLl9rV3KLj87oVGFLOaU4AAABwjG1hYtmyZZo/f77q6uokSW+++aYOOeQQjRo1Sh6PRw0NDVq1apW2bt2qcDisE088UZI0ffp0rVq1SpFIRK+++qomTZqUVo5eoUCFthMmAAAA4BCPXQe+66670tZ37NihYDCYXK+rq1NjY2NWeTAYVGNjo1paWlRVVSWPx5NWvr9qa6sG+AoK32Ejh+nZ17Zo+PAqGYaRt+cNBqvz9lwHGvrWPvStfehbe9Cv9qFv7UPf2qdQ+9a2MJHJNM20D7yWZckwjD7LE4+pBvKBeefOdpmmNfCGF7Chfq86u6L6n007NayqLC/PGQxWq6mpLS/PdaChb+1D39qHvrUH/Wof+tY+9K19nO5bl8vo8wv6vM3mFAqF1NTUlFxvampSXV1dVnlzc7Pq6uoUCATU1tamWCyWVh+9ktPDchE2AAAAHJC3MHHCCSdo48aN2rx5s2KxmFasWKFx48Zp5MiRKisr09q1ayVJTz75pMaNGyev16uxY8dq5cqVkqQnnnhC48aNy1dzi0J9oEIS95oAAACAM/J2mlNZWZkWLVqkOXPmqKurS+PHj9fkyZMlSUuWLNEtt9yi9vZ2HXvssZo9e7Ykaf78+Zo7d64eeOABjRgxQt///vfz1dyiEBhSLp/HRZgAAACAI2wPE88991xy+bTTTtPvf//7rDqjR4/W448/nlU+cuRI/fKXv7S1fcXMZRiqq/ETJgAAAOAI7oBd5EK1hAkAAAA4gzBR5EKBCjXvDisaM51uCgAAAA4whIkiFwr4ZVqWmnZ3Ot0UAAAAHGAIE0UuFKiUxPSwAAAAyD/CRJELMT0sAAAAHEKYKHL+cq+G+L3aTpgAAABAnhEmSkAo4FcjYQIAAAB5RpgoAUwPCwAAACcQJkpAKFCptr0RdYQjTjcFAAAABxDCRAmoT1yEzYxOAAAAyCPCRAkIBfySmNEJAAAA+UWYKAHBYRVyuwzCBAAAAPKKMFECPG6Xhg+rIEwAAAAgrwgTJWJEgBmdAAAAkF+EiRJRH6hQ465OmZbldFMAAABwgCBMlIhQwK9ozNSuPWGnmwIAAIADBGGiRDCjEwAAAPKNMFEiQrWVkqTthAkAAADkCWGiRAzxe1VR5lYjYQIAAAB5QpgoEYZhKMSMTgAAAMgjwkQJIUwAAAAgnwgTJSQU8GtXa5e6umNONwUAAAAHAMJECanvmdGpsYXRCQAAANiPMFFCmB4WAAAA+USYKCH1hAkAAADkEWGihJR53aodUkaYAAAAQF4QJkpMfcDPvSYAAACQF4SJEpOYHtayLKebAgAAgBJHmCgxoYBfnV0xtXZ0O90UAAAAlDjCRIkJ1XIRNgAAAPKDMFFiQjXxMLGdMAEAAACbESZKTGBoubweFxdhAwAAwHaEiRLjMgzV11Too52ECQAAANiLMFGCEjM6AQAAAHYiTJSg+oBfTbvDisZMp5sCAACAEkaYKEGhgF+mZalpd6fTTQEAAEAJI0yUIKaHBQAAQD4QJkpQKECYAAAAgP0IEyWostyrar+XGZ0AAABgK8JEiQoF/NxrAgAAALYiTJQopocFAACA3QgTJSpU61fr3oj2hiNONwUAAAAlijBRokI18YuwtzM6AQAAAJsQJkpUYnpYrpsAAACAXQgTJSo4rEIuw+C6CQAAANiGMFGiPG6XgsPKmR4WAAAAtiFMlLB6ZnQCAACAjQgTJSwU8KuxpVOmZTndFAAAAJQgwkQJC9X6FYma2tUadropAAAAKEGEiRI2IhCf0YlTnQAAAGAHwkQJq0+ECS7CBgAAgA08TjzprFmztGvXLnk88adfsGCBPvjgAz3wwAOKRqO69NJLdckll0iS1qxZo7vvvltdXV0677zzdN111znR5KI0tNKncp9bjbs6nW4KAAAASlDew4RlWdq0aZP+/Oc/J8NEY2OjrrvuOv3ud7+Tz+fTjBkzdMopp+jggw/WvHnz9Mtf/lIjRozQVVddpdWrV2v8+PH5bnZRMgxDoYBfH+3qcLopAAAAKEF5DxPvv/++JOnyyy/X7t27ddFFF6myslKnnnqqhg0bJkmaNGmSVq1apc9+9rM65JBDNGrUKElSQ0ODVq1aRZjYD6Fav97bstvpZgAAAKAE5T1MtLa26rTTTtOtt96qSCSi2bNn67zzzlMwGEzWqaur05tvvqkdO3ZklTc2Nu7X89XWVg1a24vR4QfX6K9vNap6aIXKfYPz6w4GqwflOMhG39qHvrUPfWsP+tU+9K196Fv7FGrf5j1MjBkzRmPGjEmuX3jhhbr77rt19dVXJ8ssy5JhGDJNU4ZhZJXvj50722WaB+59FoaUuyVJ/3yvSaPqPnmwCgar1dTU9omPg2z0rX3oW/vQt/agX+1D39qHvrWP033rchl9fkGf99mcXnvtNb388svJdcuyNHLkSDU1NSXLmpqaVFdXp1AolLMc/RdielgAAADYJO9hoq2tTYsXL1ZXV5fa29u1fPly3XPPPXr55Ze1a9cudXZ26o9//KPGjRunE044QRs3btTmzZsVi8W0YsUKjRs3Lt9NLmr1NYnpYRn9JqcAABtDSURBVLkIGwAAAIMr76c5fe5zn9O6dev0hS98QaZpaubMmTr55JN13XXXafbs2YpEIrrwwgt1/PHHS5IWLVqkOXPmqKurS+PHj9fkyZPz3eSiVuZzKzCkjJEJAAAADDrDsqySvqDgQL9mQpLu+dUbCnfHdOulYz/xsZw+Z6+U0bf2oW/tQ9/ag361D31rH/rWPk73bUFdM4H8C9X69dGuvSrx3AgAAIA8I0wcAEIBvzq7omrdG3G6KQAAACghhIkDwIgAF2EDAABg8BEmDgD1PWGisaXT4ZYAAACglOR9NifkX+2QcnncLr3w5ja1dnRrSKVP1X6vhlT6NMQf/ynzuZ1uJgAAAIoMYcJGsebNirz7kty1o+QKHibXsINkuPI/GORyGfrsMXV6/d0m/Wtra846ZV53esCo9Kran1j2aYjfq+pKn3wVPpmmJZdr/+5EDgAAgNJDmLCR2b5TkQ0vKBLpOb3IUyb38EPkCh4md/AwuesOl1EdlGHY/8H8a1M/LUnqjsTUtjei1r3dau3oVuve7vh6YrmjWztbw9q4vVVteyMyc8wAZRhSdUU8XCTCRrXfq6GVvmQAqa70aqjfp+pKn8q8jHoAAACUIsKEjbyHniTPZffL3PORzB0bFWuK/0T++awisWi8UlllPFgED5M7eLhcdYfJ5R9mW5t8Xrdqh7pVO7T8Y+ualqW94Wg8aPSEDcvl0tbGNrX1hJG2vRFt3Naq1r3dCnfHch6nzOvWkEpvPGT448HD53XL63bJ4zHkcbviy26XPB6XPG4jbd3rSdnuNuT1JJZdPcvxY7hdRl6CGQAAAOIIEzYzDJfcww6Se9hB8h51hiTJMqMyd21VrGmjzKb3FWvaqO6//z/JMuP7VNbIHTysdwQjeJiMssq8t91lGKqq8KqqwquDhseff183Ten3qMdHrYpETEVjpiIxU4N1+wtD6gkjLnndRspyIpikhJSU4OIyDBmGZBiG4mdvxR+NlHLDUFq9tHX1rivlOEbK9szjZR1f0pAhFWpvDyePq+RzxLcbOfZJPZaM+IwKme02crQzkbkyX1Pi0W0YcrkMuV25HxPtAgAABzbChAMMl0fu4YfIPfwQ6ZizJUlWtEux5g+S4SLWtFHRTa/37jO0PhksXMHD5R7+v2R4yhx6Bbntz6hHqphpKhq1FInFA0Y0avYsW/HAEe0pj5mKRNPL0vexksu926yU48W3dUSiyeNFe8KMZVkyex4T65bU57ZcdQ80LsOQ290TMj4mfLj7rOvqDSiJsp66qftmbc8MN32UJ+rXNLarvS3ce5x+hKWc7XL1vA7CFAAAkggTBcPwlMkTOlIKHZkss7o6FGva1DOCsVGx7e8q+j9/7dnBJVfNyN4RjLrD5QqMlOEqvl+p2+WS2yeVqXivregNH4mA0btsJsqUsZ4WTizV1FRq566O9JCScZz4c2UfM9c+piRlhh/lqKveNiT3NePPETMtmWbqo5lcz9pmpa/nrmcm1yMxU6YZTdue/Xzpx06UFYLUAOJOCyS9AcntzghSPduy9nW7skNU8scll8uQJ+uYrpT9U46dWHbH90kspz5XYtnjTt3mSgYlAAD6q/g+eR5AjLJKeQ4+Vp6Dj02WmXt391x/ER/BiGxaK214Pr7R7ZWr9n/1jmDUHSbX0JAMg9uJ2C1x2pE+wQexYG2l3KY5aG0qZabVR+gwLVk5QtCQoRXauasjPej0EX6SjzmeIzNMJbfF4seLxfrYnrJvJGaqK9Kzz76OmRK8YrH8jX4ZhjLCSGZwSQ8x5WVexaKx5IhO6khO6qiRy5UIYC65DGWPAPUsG5n7G0qrkzlalDYy1Y8Ales1MNIEAANHmCgyLv8wuQ4dI8+hYyT1fLvc1pQ8Ncps2hifQeqtP8V38FbIHTw0ft3FkDoZ5VUyyqrij4llN28DFBeXYcjlNuTp52BWMFitIWXFO/Il9YwUxXrDR9TsXY+Z6eEksRyNpW7rXY7uY5+05Z799vVcbpehiGkpGksPYWbPcixlOS20WUord3rEKdcIT2ow6R0t6m/QcvWM/mQcM6OOJ+N5E6NFtS1htbV1powepbYrfUTJ4yYUAXAOnyKLnGEYMobUyTWkTt4jTpEkWaYpc/f2tOsvuv/xtGTmnm1J3vL0kFFWJaO8Mjt0lFcp4q2X1W3E9+EfLSBvXIYhl8eQV4U10rivSRn2V2bgSI4y9QSPmGn2hpCMwNI7ihMPPKkBKG3kJ9Z3OIqHtNwBLJZy7MRPJBJTLBxNO3bmcjRxnDyNLiUDSWLUJWMUJ3OUJ2uboeS1TEbqqFDKtUa5rlPKPkbPKXqJUShXz+itJBmpy0oup/6bYmQsJCaQSJWYWCK9Xvo+iXq9+xgauq1NbW1huVyJSSziI2eJMJY6cmakbEtMWJHc5jLkUnw9ORlGoi8y6qY+F1BqCBMlyHC55A6MlDswUt6jz5IkWbGorM49ssLtsro6eh7b44/hdlnhtmS52bpDVrhd6t6bdexkicudI2xUpgUSZQQRo6xShqu4vx0GYJ/EiFMRXz61T2bqaFEi2KSEltTgUT2kQs0729NHiFK2RzPCSzSWPrqUHCkyTZmmskaNMk/ls6z4MeIBrfdaJtNS2nHSRpt6The0UsJcYVzRVNhSw0XqTHtGTxJKnb1PKcuJoORKLVdGvUToypwFsGff5HMl1lPqZco102LO329KocfrUjRiKuc7oZ/HM6S04Jsr0OaePMOVtT3XhBx9ThSSY3vWl6ZGzsW05dTONHIXZ7/mPjYmin0et4YPr+r7AA4jTBwgDLdHRlWtVFXb730sMxYPGF3tssIdUrhdVb6oWpuaUsJIfLu55yNZjfGyPkdAJMnnj4cKt1fxv6bu+KPLHb+2w+WWUh4NlzutnuFyp+2TrJPYp6eukXHsfdcz4uXJn951I22b0bNffNnIUT+xbKSWu1yK/wvQU0ecigAciFwuQ75+fqESDFarqdpnc4sGX+ZpbqmBRJKUMvtd6odVK2UlsZisafV+6LRSKvUu9z5YGZ+ALSt9H0kaVlOpnTvbk5NZJNpqpZ56l7lu9k5UYVqWrJ5glazXE6oSI2eJiS8Sx0+umxn79Cwn+i75Wi3JVOp6fH8lJ93o6Z+USTQS+6VNtJHoq4xypeyTOJah7H+X+vtPVaKar8yj7q5oz7792zmzWqJvUq85i0bMtPXe7b0jj9nbEkG5f6+h0C28+gyFhhbWLJ4JhAn0yXC5ZVQMkSqGJMuqg9UK7+OUBsuypEg4LWhY4bbeEZCu+MiIzJhkmrLMWPz+GimPltmdsmxKVu9j/KuxeN20fRPbi0UygPQGkQ63R5bhklyensDjjs/O5XJL7pTlzG0ud/y6l+Q2Tzw4uXMcJ+2Ynp51d+9zuj29oSvZPqPvAJUIST2Phoyc4ap3f0IUUMqKYXQpGKxWpYe/RXYYzNMeB0uuEbk+J95IrROz0kZY+golmQE2V/2sGqnhWX3US1nxelw67ohaNTe37+ulOoYwgUFlGIbkq5Dhq5Cqg3l9biv5dU1KuOgzsCS2x3r2MXv2N3t+epbN3nUrc1vKspW5j6zkvjmPm7mvaaqi3K3Ojs54+2LRnpAUk2X2LsuMyYp29dSJSWY0/vrM7GXFYupjUNpZaSEj/dHILEsZ0TESozxGxo8rcz9X+j6GS40VPnV1m2n7Ga7MEanUY2U/V29A88TDm9sbD19uj+RKLMcf46EspU5i2eUmUAFAHhVDwO2PQv63gzCBkmEkTzp1pf3RKNz//dINt+EbHcs048EiRzBJDSxpQSQWzQ5QiXCU8WjlKrcyA5iVEaSs3CEsR53eAJca7MyeMJgjoMWiafsk6nW3SbFoNC0cph7XyvEciTvSD7qe8JEdSDy9gSR1Oa1+T1nqqJEyHjPKjMyyfSwb/aiTXt+lvXuqFG3r7h3VSo54ebJHvFweyR0/VbGQ/2EEAPQfYQIoYYbLJbl6z7s+UD++DXToPS1kJMJYLP5jmZHe5VjPshmJh7RYVIr1LJspdcz0+lbPPr3bY/H9ImEpHH+O7ONFFD95uieUOTz61DnQHXMFjeQpe/vallj3yHBnnNrXx+iV4Uoffco6bS9xjVXPdVFGjhGrfY+OpR6r99qs9OvAep+jVIPUPkeHZcnwlEneMu59BJQYwgQA9KH3A2LPurPN6VPyQ5x6R316NqQs95ZZOcrS6iWv0DTjh1RiNEeS4mWJesOGlqllZ2va6XW9I2DRlBGwaPLUvP3eljy1L5q9LWtErQhu/NiPySHCXo9ipvEx9TLKEiEmGVxcio8appzamRKOs649S6ynXqOWODU041TR+Ihejrr94SmT4S2PTzHuLZfhLcu5rJ51w1sRDyF9lbu9tv66crEsS4pFpGh3z5cDEVnR7rRHRSOyYt29daKRnrrdvXWSp9r2NZLbx+m0OR6ztsvs+eLB7HskOTmibMX71Vcuw+fvWa6Q4e05bTmx7vNLvvJkueGtiK/7KiR38U0YgMFBmACAItd7KlM/6w/ic5cHq+XxFdYFl+mnrfV+KEt+cM66JirzVLeP3yf3KXhW7wfvjMkh+juJRKKez+tSV2dX7g/80a6s+lbKcVOPn2t0ZF+z3RmesuygkhJgDKM38OSeia+P55CkSJesSDg+8pa2HI5P1NG6Q1a0S1Z3WIqGe0Ptx3G5cwcTT2Ywif/sGeJX1+42KZby4T/rg340JSh0y4qmBoae7Z+E2yd5vPE+zpjIQqkTWrhcknJfZ5Y+a2DGes4ZBzOvUUuZMMNST993SpFOWR27Ze7eLnV3xsvMfrxew6295RWyPOU5QkiFlAggvor47yQlrBi+iuS6PL6SHb0rVYQJAEBJyRxRSpY705wBKcRZcfIp/s1/t6xIVzxsRMI9y50ZZYlgkl1udbb1hpVIOD4SIKkr9Yl6JkcwPN74h1i3t/eDvtsro8zfU+aV4fGlPcbr+Hr369lHWWW9x5PHFz8tr8g+LFuxSE/QCMtKBIxIZ3LZinRK3WGVu6Pa29oaDyGRsKzOVpmtjb2hpOd3sP+MlAcjR9k+1vusl/k7yPEcafeMSPnSxkg5dlpZ6o0ocpQnrztLHjFje+I/vXVkSPKUadj0b0iqViEiTAAAgIJiGEb8dChPWdr05J+EZcakaJdqayq0c3d3PBi4uH6jPwy3V0aF92N/Fx83kYgVi/YGkkiOUNIdlqJdmXtlj1JljVpZOct7p239mP1z1Usts1K2JU8pTZT1LvceNqNOrmNlPqeV+RpS6iRCbIHOgE+YAAAAJc9wuSWfX25/tYyOA3fUx0nxWemqZJQX7t2cC5VnSLVUoKOVRHIAAAAAA0KYAAAAADAghAkAAAAAA0KYAAAAADAghAkAAAAAA0KYAAAAADAghAkAAAAAA0KYAAAAADAghAkAAAAAA1Lyd8B2uQynm1By6FP70Lf2oW/tQ9/ag361D31rH/rWPk727b6e27Asy8pjWwAAAACUCE5zAgAAADAghAkAAAAAA0KYAAAAADAghAkAAAAAA0KYAAAAADAghAkAAAAAA0KYAAAAADAghAkAAAAAA0KYAAAAADAghAkAAAAAA0KYQE733XefpkyZoilTpmjx4sU5t3/uc5/TtGnTNG3aND366KMOtLI4zZo1S1OmTEn23bp169K2r1mzRg0NDZo4caLuvfdeh1pZfH7zm98k+3TatGk6+eSTtWDBgrQ6vG/3T3t7u6ZOnaoPP/xQUv/em9u2bdMll1yiyZMn6+qrr1ZHR0c+m1w0Mvv2scce09SpU9XQ0KCbbrpJ3d3dWfssX75cZ555ZvL9y9+H3DL79qabbtLEiROT/fbMM89k7fP2229r+vTpmjRpkm6++WZFo9F8N7vgpfbr6tWr0/7ennrqqbrqqquy9uE9+/Fyfd4qur+1FpDhpZdesr785S9bXV1dVnd3tzV79mzrj3/8Y1qdq666ynr99dcdamHxMk3TOvPMM61IJJJze2dnpzV+/Hjrgw8+sCKRiHX55Zdbf/nLX/LcyuL37rvvWhMmTLB27tyZVs77tv/+/ve/W1OnTrWOPfZYa8uWLf1+b1555ZXWihUrLMuyrPvuu89avHhxvpte8DL79v3337cmTJhgtbW1WaZpWjfeeKP10EMPZe23YMEC66mnnsp/g4tIZt9almVNnTrVamxs3Od+U6ZMsd544w3Lsizrpptush599FHb21pMcvVrwo4dO6xzzz3X2rhxY9Z+vGf3Ldfnraeeeqro/tYyMoEswWBQc+fOlc/nk9fr1RFHHKFt27al1Vm/fr1+/OMfq6GhQQsWLFBXV5dDrS0u77//viTp8ssv1/nnn69HHnkkbfubb76pQw45RKNGjZLH41FDQ4NWrVrlRFOL2ne/+11dd911CgQCaeW8b/tv2bJlmj9/vurq6iT1770ZiUT06quvatKkSZKk6dOn8/7NIbNvfT6f5s+fr6qqKhmGoaOOOirrb64k/eMf/9Dy5cvV0NCgb3/729qzZ0++m17wMvu2s7NT27Zt07x589TQ0KAf/vCHMk0zbZ+tW7cqHA7rxBNPlMT7NpfMfk21ePFizZgxQ4ceemjWNt6z+5br89amTZuK7m8tYQJZjjzyyOQf1U2bNukPf/iDxo8fn9ze0dGhY445RjfccIOWL1+u1tZWLV261KnmFpXW1laddtppuv/++/Xwww/r17/+tV566aXk9h07digYDCbX6+rq1NjY6ERTi9aaNWsUDod13nnnpZXzvt0/d911l8aOHZtc7897s6WlRVVVVfJ4PJLi/1Dy/s2W2bcjR47UGWecIUnatWuXHn30UZ177rlZ+wWDQf3Hf/yHfv/732vEiBFZp/Ehu2+bm5t16qmnauHChVq2bJlee+01Pf7442n7ZL63ed9my+zXhE2bNulvf/ubZs+enXM/3rP7luvzlmEYRfe3ljCBPr333nu6/PLLdeONN6Z941BZWakHH3xQRxxxhDwejy6//HKtXr3auYYWkTFjxmjx4sWqrq5WIBDQhRdemNZ3pmnKMIzkumVZaev4eL/+9a/11a9+Nauc9+0n05/3Zq4y3r/919jYqEsvvVRf/OIXdcopp2Rtv//++3XyySfLMAx97Wtf0wsvvOBAK4vLqFGjdP/996uurk4VFRWaNWtW1v/3/N0duMcee0wzZ86Uz+fLuZ33bP+kft4aNWpU0f2tJUwgp7Vr1+qyyy7T9ddfrwsuuCBt27Zt29K+2bEsK5mOsW+vvfaaXn755eR6Zt+FQiE1NTUl15uamnIOKyO37u5uvfrqqzrnnHOytvG+/WT6894MBAJqa2tTLBbrsw5y+9e//qUZM2boggsu0DXXXJO1va2tTQ8//HBy3bIsud3uPLawOG3YsEFPP/10cj3X//eZ7+3m5mbet/307LPP6vOf/3zObbxn+yfz81Yx/q0lTCDL9u3bdc0112jJkiWaMmVK1vby8nLdc8892rJliyzL0qOPPqoJEyY40NLi09bWpsWLF6urq0vt7e1avnx5Wt+dcMIJ2rhxozZv3qxYLKYVK1Zo3LhxDra4uGzYsEGHHnqo/H5/1jbet59Mf96bXq9XY8eO1cqVKyVJTzzxBO/ffmhvb9cVV1yhb3zjG7r88stz1vH7/frpT3+anP3tkUce4f3bD5ZlaeHChdqzZ48ikYgee+yxrH4bOXKkysrKtHbtWknSk08+yfu2H3bt2qVwOKxRo0bl3M579uPl+rxVjH9r+VoOWX72s5+pq6tLixYtSpbNmDFDzz33nK699lp95jOf0YIFC3T11VcrEonopJNOynlaCbJ97nOf07p16/SFL3xBpmlq5syZGjNmjKZNm6af/OQnqq+v16JFizRnzhx1dXVp/Pjxmjx5stPNLhpbtmxRKBRKK/v617/O+3YQlJWV9fnevPnmm3XOOefo3HPP1fz58zV37lw98MADGjFihL7//e873PLC9/jjj6u5uVkPPfSQHnroIUnSOeeco2984xtpffuDH/xA3/3udxUOh3XooYfmnLYb6UaPHq0rr7xSF198saLRqCZOnKipU6dKSv/bsGTJEt1yyy1qb2/Xscce2+c1AOj14YcfZv29lcR7dj/09Xmr2P7WGpZlWY49OwAAAICixWlOAAAAAAaEMAEAAABgQAgTAAAAAAaEMAEAAABgQAgTAAAAAAaEqWEBAPvt6KOP1lFHHSWXK/07qfvvv18HH3zwoD/Xyy+/rEAgMKjHBQB8coQJAMCA/OIXv+ADPgAc4AgTAIBB9corr2jJkiU66KCD9P7776u8vFyLFi3SEUccoba2Nt1+++165513ZBiGzjrrLH3rW9+Sx+PRunXrdOedd6qzs1Ner1c33nijTjvtNEnSj370I61bt067d+/WFVdcoUsuuURNTU36zne+o5aWFknS+PHj9c1vftPJlw4ABxyumQAADMill16qadOmJX+uueaa5Lb169dr1qxZeuqppzR9+nTdcMMNkqQ777xTw4YN01NPPaXf/va32rBhg37+858rEonommuu0TXXXKMVK1bojjvu0MKFC2WapiRp1KhR+t3vfqf77rtPixYtUiQS0bJly3TwwQdr+fLlevTRR7V582a1tbU50hcAcKBiZAIAMCD7Os1p9OjRGjt2rCTpi1/8ohYsWKCWlhY9//zz+tWvfiXDMOTz+TRjxgz94he/0BlnnCGXy6Wzzz5bknTcccfpqaeeSh5v6tSpkqRjjjlG3d3dam9v11lnnaUrr7xS27dv1+mnn67rr79e1dXV9r5oAEAaRiYAAIPO7XbnLDNNU4ZhJMtM01Q0GpXb7U4rl6R3331X0WhUkuTxxL/7StSxLEvHH3+8nn32WX35y1/W1q1b9aUvfUnr16+36yUBAHIgTAAABt0777yjd955R5L02GOPacyYMRoyZIjOPPNMPfLII7IsS93d3Vq2bJlOP/10HX744TIMQy+99JIk6a233tKll16aPM0plyVLlmjp0qX693//d91888361Kc+pffeey8vrw8AEGdYlmU53QgAQHHpa2rYb33rWyovL9d3vvMdjR49Wlu3blUgENBdd92lgw8+WC0tLbrzzju1YcMGRSIRnXXWWbrxxhvl8/n0j3/8QwsXLtTevXvl9Xo1d+5cjR07Nmtq2MR6LBbT3Llz1djYKJ/Pp6OPPlq33367fD6fE10CAAckwgQAYFC98soruuOOO7RixQqnmwIAsBmnOQEAAAAYEEYmAAAAAAwIIxMAAAAABoQwAQAAAGBACBMAAAAABoQwAQAAAGBACBMAAAAABuT/A1TkI989xglcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 936x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_loss(history)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train set RMSE:17.433590186466215, R2:0.8249868754417847\n"
     ]
    }
   ],
   "source": [
    "# predict and evaluate\n",
    "y_hat_train = model.predict(X_train_smooth[train_cols])\n",
    "evaluate(y_train_clipped, y_hat_train, 'train')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test set RMSE:18.681100442307972, R2:0.7979097473627429\n"
     ]
    }
   ],
   "source": [
    "# only run the test set once you're happy with the preprocessing\n",
    "y_hat_test = model.predict(X_test_smooth[train_cols])\n",
    "evaluate(y_test, y_hat_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately smoothing the data for FD001 makes the predictions worse, but that's besides the point. The point is, you can go back and forth between preprocessing and training, trying different settings and when you want to use a previous setting you can be assured the results are the same."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial NN developed here (without smoothing or feature selection) already has an RMSE of 18.50 on the FD001 dataset. A whopping 42% improvement over the baseline model and a big improvement over the previous best, which was a SVR with an RMSE of 20.54 (and an overall improvent of 35.7% on the baseline model). In my next article, we'll delve into FD002 in which the turbofan engines run on different operating conditions and the exponential smoothing does have a beneficial effect."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}