2_CNN_Synthetic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Machine Learning Applied to Synthetic Sitelle Data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import fits\n",
    "import tensorflow as tf\n",
    "from keras.backend import clear_session\n",
    "from keras.models import Sequential\n",
    "from keras.layers import Dense, InputLayer, Flatten, Dropout\n",
    "from keras.layers.convolutional import Conv1D\n",
    "from keras.layers.convolutional import MaxPooling1D\n",
    "from keras.optimizers import Adam\n",
    "from keras.callbacks import EarlyStopping, ReduceLROnPlateau\n",
    "from pickle import dump\n",
    "from keras.backend import clear_session\n",
    "from sklearn.preprocessing import StandardScaler, MinMaxScaler, RobustScaler\n",
    "import numpy as np\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set Input Parameters\n",
    "output_dir = '/home/carterrhea/Desktop/Test/'  # Include trailing /\n",
    "syn_num = 1000  # Number of Synthetic Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Spectra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# IF RESOLUTION IS SAMPLED USE THIS \n",
    "# This gets us the wavenumbers we are going to interpolate on\n",
    "# We use this spectra because it is a nice looking (can easily see the lines -- not relevant here)\n",
    "# spectra at a resolution of 5000 exactly. We want to keep that sampling :)\n",
    "spec98_ = fits.open('Reference-Spectrum.fits')[1].data\n",
    "channel = []\n",
    "counts = []\n",
    "for chan in spec98_:  \n",
    "    channel.append(chan[0])\n",
    "    counts.append(chan[1])\n",
    "min_ = np.argmin(np.abs(np.array(channel)-14500))  # note that we cut the values a bit more than we did during the creation of the synthetic data\n",
    "max_ = np.argmin(np.abs(np.array(channel)-15600))  # We do this to make sure all values are within the range and to avoid any extrapolation!\n",
    "wavenumbers_syn = channel[min_:max_] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in Fits -- This is not particularly fast :S\n",
    "Counts = []  # List containing the spectra\n",
    "Param_dict = {}  # {spectrum_ct: [velocity, broadening]}\n",
    "for spec_ct in range(syn_num):\n",
    "    spectrum = fits.open(output_dir+'Spectrum_%i.fits'%spec_ct)\n",
    "    header = spectrum[0].header\n",
    "    spec = spectrum[1].data\n",
    "    channel = []\n",
    "    counts = []\n",
    "    for chan in spec:  # Only want SN3 region\n",
    "        channel.append(chan[0])\n",
    "        counts.append(chan[1])\n",
    "    # We need to find the min and max wavenumber that match our min/max from our synthetic spectra    \n",
    "    # We go a bit further on either end for safety during the interpolation (so we don't need any extrapolation)\n",
    "    min_ = np.argmin(np.abs(np.array(channel)-14400))\n",
    "    max_ = np.argmin(np.abs(np.array(channel)-15700))\n",
    "    f = interpolate.interp1d(channel[min_:max_], counts[min_:max_], kind='slinear')  # interpolate using a linear spline\n",
    "    # Get fluxes of interest                                                                     \n",
    "    coun = f(wavenumbers_syn)\n",
    "    Counts.append(coun) \n",
    "    Param_dict[spec_ct] = [header['VELOCITY'], header['BROADEN']]  # Save velocity and broadening\n",
    "vel_spec = [val[0] for val in list(Param_dict.values())]  # Read velocity into list\n",
    "broad_spec = [val[1] for val in list(Param_dict.values())]  # Read broadening into list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup Machine Learning Algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/carterrhea/anaconda3/lib/python3.7/site-packages/keras/callbacks/callbacks.py:998: UserWarning: `epsilon` argument is deprecated and will be removed, use `min_delta` instead.\n",
      "  warnings.warn('`epsilon` argument is deprecated and '\n"
     ]
    }
   ],
   "source": [
    "activation = 'relu'  # activation function\n",
    "initializer = 'he_normal'  # model initializer\n",
    "input_shape = (None, len(Counts[0]), 1)  # shape of input spectra for the input layer\n",
    "num_filters = [4,16]  # number of filters in the convolutional layers\n",
    "filter_length = [8,4]  # length of the filters \n",
    "pool_length = 2  # length of the maxpooling\n",
    "num_hidden = [256,128]  # number of nodes in the hidden layers\n",
    "batch_size = 2  # number of data fed into model at once\n",
    "max_epochs = 25  # maximum number of interations\n",
    "lr = 0.0007  # initial learning rate\n",
    "beta_1 = 0.9  # exponential decay rate  - 1st\n",
    "beta_2 = 0.999  # exponential decay rate  - 2nd\n",
    "optimizer_epsilon = 1e-08  # For the numerical stability\n",
    "early_stopping_min_delta = 0.0001\n",
    "early_stopping_patience = 4\n",
    "reduce_lr_factor = 0.5\n",
    "reuce_lr_epsilon = 0.009\n",
    "reduce_lr_patience = 2\n",
    "reduce_lr_min = 0.00008\n",
    "loss_function = 'mean_squared_error'\n",
    "metrics_ = ['accuracy', 'mae']\n",
    "# Clear session for sanity :)\n",
    "clear_session()\n",
    "model = Sequential([\n",
    "    InputLayer(batch_input_shape=input_shape),\n",
    "    Conv1D(kernel_initializer=initializer, activation=activation, padding=\"same\", filters=num_filters[0], kernel_size=filter_length[0]),\n",
    "    Conv1D(kernel_initializer=initializer, activation=activation, padding=\"same\", filters=num_filters[1], kernel_size=filter_length[1]),\n",
    "    MaxPooling1D(pool_size=pool_length),\n",
    "    Flatten(),\n",
    "    Dropout(0.2),\n",
    "    Dense(units=num_hidden[0], kernel_initializer=initializer, activation=activation),\n",
    "    Dense(units=num_hidden[1], kernel_initializer=initializer, activation=activation),\n",
    "    Dense(2),\n",
    "])\n",
    "# Set optimizer\n",
    "optimizer = Adam(lr=lr, beta_1=beta_1, beta_2=beta_2, epsilon=optimizer_epsilon, decay=0.0)\n",
    "# Set early stopping conditions\n",
    "early_stopping = EarlyStopping(monitor='val_loss', min_delta=early_stopping_min_delta, \n",
    "                                       patience=early_stopping_patience, verbose=2, mode='min')\n",
    "# Set learn rate reduction conditions\n",
    "reduce_lr = ReduceLROnPlateau(monitor='loss', factor=0.5, epsilon=reuce_lr_epsilon, \n",
    "                                  patience=reduce_lr_patience, min_lr=reduce_lr_min, mode='min', verbose=2)\n",
    "# Compile CNN\n",
    "model.compile(optimizer=optimizer, loss=loss_function, metrics=metrics_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To test no scaling (other than the normal scaling), we can just use the data we read in since our synthetic spectra have already been scaled by the maximum value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train on 700 samples, validate on 200 samples\n",
      "Epoch 1/10\n",
      " - 3s - loss: 6926.8499 - accuracy: 0.5971 - mae: 61.3165 - val_loss: 6321.1572 - val_accuracy: 0.5750 - val_mae: 56.1029\n",
      "Epoch 2/10\n",
      " - 1s - loss: 5147.5604 - accuracy: 0.6329 - mae: 49.2847 - val_loss: 3530.2183 - val_accuracy: 0.8350 - val_mae: 40.9384\n",
      "Epoch 3/10\n",
      " - 1s - loss: 1794.5521 - accuracy: 0.8971 - mae: 28.3811 - val_loss: 576.1165 - val_accuracy: 0.9650 - val_mae: 16.9546\n",
      "Epoch 4/10\n",
      " - 1s - loss: 603.7346 - accuracy: 0.9343 - mae: 17.4610 - val_loss: 514.4747 - val_accuracy: 0.9650 - val_mae: 16.5065\n",
      "Epoch 5/10\n",
      " - 1s - loss: 528.3834 - accuracy: 0.9557 - mae: 16.4571 - val_loss: 430.6654 - val_accuracy: 0.9650 - val_mae: 15.2681\n",
      "Epoch 6/10\n",
      " - 1s - loss: 441.0608 - accuracy: 0.9571 - mae: 15.3234 - val_loss: 385.7453 - val_accuracy: 0.9650 - val_mae: 14.5634\n",
      "Epoch 7/10\n",
      " - 1s - loss: 426.7855 - accuracy: 0.9557 - mae: 15.1410 - val_loss: 352.6820 - val_accuracy: 0.9650 - val_mae: 13.8803\n",
      "Epoch 8/10\n",
      " - 1s - loss: 387.4838 - accuracy: 0.9443 - mae: 14.6729 - val_loss: 369.2998 - val_accuracy: 0.9600 - val_mae: 14.6632\n",
      "Epoch 9/10\n",
      " - 1s - loss: 352.5138 - accuracy: 0.9571 - mae: 14.0588 - val_loss: 323.5518 - val_accuracy: 0.9450 - val_mae: 13.8583\n",
      "Epoch 10/10\n",
      " - 1s - loss: 354.7077 - accuracy: 0.9486 - mae: 14.0730 - val_loss: 298.3474 - val_accuracy: 0.9500 - val_mae: 13.2974\n",
      "Model: \"sequential_1\"\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "conv1d_1 (Conv1D)            (None, 343, 4)            36        \n",
      "_________________________________________________________________\n",
      "conv1d_2 (Conv1D)            (None, 343, 16)           272       \n",
      "_________________________________________________________________\n",
      "max_pooling1d_1 (MaxPooling1 (None, 171, 16)           0         \n",
      "_________________________________________________________________\n",
      "flatten_1 (Flatten)          (None, 2736)              0         \n",
      "_________________________________________________________________\n",
      "dropout_1 (Dropout)          (None, 2736)              0         \n",
      "_________________________________________________________________\n",
      "dense_1 (Dense)              (None, 256)               700672    \n",
      "_________________________________________________________________\n",
      "dense_2 (Dense)              (None, 128)               32896     \n",
      "_________________________________________________________________\n",
      "dense_3 (Dense)              (None, 2)                 258       \n",
      "=================================================================\n",
      "Total params: 734,134\n",
      "Trainable params: 734,134\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "syn_num_pass = len(Counts)\n",
    "# Create divisions\n",
    "train_div = int(0.7*syn_num_pass)  # Percent of synthetic data to use as training\n",
    "valid_div = int(0.9*syn_num_pass)  # Percent of synthetic data to use as validation\n",
    "test_div = int(1.0*syn_num_pass)  # Percent of synthetic data to use as testing\n",
    "# Set training set \n",
    "TraningSet = np.array(Counts[:train_div])\n",
    "TraningSet = TraningSet.reshape(TraningSet.shape[0], TraningSet.shape[1], 1)\n",
    "Traning_init_labels = np.array((vel_spec[0:train_div], broad_spec[0:train_div])).T\n",
    "# Set validation set\n",
    "ValidSet = np.array(Counts[train_div:valid_div])\n",
    "ValidSet = ValidSet.reshape(ValidSet.shape[0], ValidSet.shape[1], 1)\n",
    "Valid_init_labels = np.array((vel_spec[train_div:valid_div], broad_spec[train_div:valid_div])).T\n",
    "model.fit(TraningSet, Traning_init_labels, validation_data=(ValidSet, Valid_init_labels), \n",
    "      epochs=10, verbose=2, callbacks=[reduce_lr,early_stopping])\n",
    "model.save('CNN_trained.h5')\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply on test test\n",
    "TestSet = np.array(Counts[valid_div:test_div])\n",
    "TestSet = TestSet.reshape(TestSet.shape[0], TestSet.shape[1], 1)\n",
    "TestSetLabels = np.array((vel_spec[valid_div:test_div], broad_spec[valid_div:test_div])).T\n",
    "test_predictions = model.predict(TestSet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEPCAYAAACZcRnqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd5wV1fXAv2c7fekgXUGaFTfYYgV7FDUWNBoxJiaWRJNobGn6U6MpErsxNmxBo1HRaFTU2IJIU6SIrDSXIr3Dtnd+f8x9Mjxemd19++r5+pnPm7nlzHmz+M7ce889R1QVwzAMw0gFBelWwDAMw8gfzOgYhmEYKcOMjmEYhpEyzOgYhmEYKcOMjmEYhpEyitKtQHPRqVMn7du3b7rVMAwjC5g2bdpqVe2cLHkiIicf1yp01SXtOXzUV5IsubmA5KrLdEVFhU6dOjXdahiGkQWIyDRVrUiWvAlP7KZ3PbQeQXjz3S1mdHzY9JphGEYSERG57a51/P3PXenQvoB3X+yVm2/2jcSMjmEYRhJ56fHuob0GldCnVzG//nkHbh67Nt0qZRRmdAzDMJJEeJRz3c86ADB0YKmNdiIwo2MYhpEk/KOcMDba2RkzOoZhGEkgcpQTZujAUtqX22gnjBmdFDDjvblUzlySbjUMw2hGoo1ywthoZwc5u08nU1j65ddcd/pYylqV8sKiuxAx70nDyDVERA6uKOOp+7pFrd9r0I7RzhGn5ve+HRvpNDMdu5dTMWIoR50+3AyOYeQuJ8Ua5YSx0Y6HjXSambKWpdz8zBXpVsMwjObl9B+c2y5ug70GlbJ1WwgRaaeqG1KkV8ZhIx3DMIymI2WliWcyXJu8nvKwkY5hGEYSqNN6arUubpsQ5sBmRscwDCMJhFAzKgGw6bUMZPG85Tzyh5dZv3pTulUxDCMgSijhf2pGyUY6mchTf32d9//9CaUtivnelcenWx3DMAJQqyFqNRS3TY4G9W8QZnQykNN+eARlrUoY+d1vpVuVhHz4ynQWzq7inF+eRGFRYbrVMYy0EQLqE4xkzOaY0clIBh/Qj8EH9Eu3GoG4+5dPsn7VJvY+dCD7fntgutUxjLQRbE3HzI4ZHaNJXHzzWXz52VcMGb5HulUxjLRSp0ptgvmz+JNv+UFGGh0RKQSmAktV9Tsi0g8YD3QApgPnq2pNOnU0PI4+8yCOPvOgdKthGGknhNr0WgAy1XvtCmCu7/p2YKyqDgDWARelRSvDyANqqutYsXx9utXIOuqBeo1/GBlodESkJ3AS8JC7FuBo4DnXZBxwanq0M4zc56YbX+C8c+9j2tSF6VYlq1C86bN4h9mdzJxe+yvwK6CNu+4IrFf9ZqtvFdAjWkcRuRi4GKB3797NrKZh5Cbl5S0pLi6iVevSdKuSVdSpUKvxI9yYy3SGjXRE5DvASlWd5i+O0jTqn05VH1TVClWt6Ny5c7PoaBi5zlW/OolXXr2KQYN2S7cqWUU9kvAwm5N5I51DgVNE5ESgDGiLN/IpF5EiN9rpCSxLo46GkfMUFmXU+2hWEHKGJR6a37E+gQwb6ajqdaraU1X7AqOBt1X1e8A7wBmu2QXAS2lS0TAMIyre9FpB3MNGOhlmdOJwDfALEanEW+N5OM36GIZh7EQ9BQmPICMdETleROaJSKWIXBulvlREnnH1k0Wkr6/uOlc+T0SOSyRTRPo5GfOdzBJXPlZEPnHHFyKSNHfGjDU6qvpfVf2OO1+gqsNVtb+qnqmq1enWzzCM5uHlVz/h4sseY8HCVelWpUGEVBIeiXB7FO8FTgCGAOeIyJCIZhcB61S1PzAWb0sJrt1oYChwPHCfiBQmkBl1O4qq/lxV91PV/YC7gX81+sFEkLFGxzCM/OSjyV8yv/JrZs9dmm5VGkQoOY4Ew4FK96Jdg7cpflREm1F4W0fA20oywm0tGQWMV9VqVV0IVDp5UWU2YDvKOcA/Aj+IBGSaI4FhGHnOL684jmNH7sW3DxmQblUaRJ0WUqvxf1IDTK/1AL7yXVcBB8Zqo6p1IrIBb9mhB/BRRN/w9pJoMhNuRxGRPkA/4O1EigfFjI5hGBlFhw6tOeKw7AseWx/Iew2AC0TkfF/xg6r6oDsPskUkVptY5dFmtOK19zMaeE5V66O0bRRmdAzDyAnq60NUfbWWPn07peX+IRXqNf6KhXrrOuNU9c4YTaqAXr7raFtEwm2qRKQIaAesTdA3WvlqEm9HGQ1cFvdLNRBb0zEMIyf4231vcdH3/8aEF6clbtwM1FFIbYIjwPTaFGCA8yorwfvRnxDRZgLe1hHwtpK8rarqykc777Z+wADg41gyXZ+Y21FEZCDQHpjUqAcSAxvpGIaRE3Ts1JqCQqF9+1ZpuX89AUY6CWS4NZrLgdeBQuARVZ0tIjcBU1V1At6WkSfcFpK1eEYE1+5ZYA5QB1wWnhaLJtPd8hpgvIjcDMxg5+0o5+A5JiR1e5EkWV7GUFFRoVOnTk23GoZhpJBQSCkoaPiufxGZpqoVjb2viDz613/vOabvoBZx2/3mvC/5bNLm9qqat2G8bXrNMIwG8+rkuZzy60eY/PmSdKuyE40xOMmiVguo0cK4h4XBMaNjGEYjmFG5lKrVG5i9aEW6VckYQlqQ8MjNeaWGYWs6hmE0mJ9/93AO23t3Dh7aJ92qZAwhF+omPjbSMaNjGEaDaVlWwuH77J5uNTKKehXqLZ9OQmx6zTCMndhWU8uNz77JP/83M92qZBXhiATxjpCNdMzoGEa+s2V7Dfe/MYnZVV8DMHvJ1zz/0Szufz2p2zNyHkUIJTgMm14zjLznlelzue+Nj5jyZRWPXHImw3bvwdWjjmBA947pVi2rCLZPxwyPGR3DyHNG7NWfmUtWcNy+ewKe2/H5RwxLs1bZhze9Vhi3jS3pmNExjLynU9tW3DL6uMQNjbh4qQ1spJMIW9MxjBTy2pTP+cP4t9lWXZtuVYwkE2SfDgESueU6ZnQMI4XcO+FDnn3vU6bOr0rpfResXktNXV3ihkajCU+vxTtC6VYyA2jQ9JqI9AS6usuvVTW1/+cYAKgqS1aso3e39njJ/4xs4YZzRjJr8QoOGtw7Zfd8c24lP332ZU7cayB3fPfElN0336iHAPl07P/XhCMdETlGRJ4SkRXAYrxQ2R8Di0VkhYg8KSIjm1tRYwePvfIxZ177GI++/HG6VTEayMFD+vCjEw6kuDD+gnMyad+yjNKiQrq1aZ2ye+YjGmR6zYg90hGRk4FbgSHhoijNuuCFvz5HRGYBN6jqK0nX0tiJ9m1aIgId2saPaGsYABV9evLJ9T+1UXEzU68F5jIdgHjTa+FkPvXA+3i5t2cBa/AMUEdgL+Ag4NvA3sCLCWQajrq6ev736SKGDe5J65alDep76pF7c/LhQykssDcnIxhmcJqfOgoSu0ybz3RcA/EVMBZ4UlVXxxMiIp2B84GfJVG3nGb8f6Zzz/j3Oemwofzmxw13VzWDYxiZhaoQShR7zUY6cY3O7uGsc4lQ1VXAHSISK++3EcFe/bvTu1t7DhjSK3FjwzAynvpAUaaNmEYnqMFpap98Zb9BPXn2zxemWw3DMJJEvRZQlzAigY10ApllERksIqeLyAHuepiIvC4i00XkFhEx824YRl4TTm0Q77AlneCbQ38P/BM4SEQKgQnASGA/4Fp3GEmitq6eiTO+YMXajelWxTCyhtlzlnLrn15h6bJ1abl/eE0n3hHEkUBEjheReSJSKSK7/LaKSKmIPOPqJ4tIX1/dda58nogcl0imiPRzMuY7mSW+urNEZI6IzBaRpxv/ZHYmqNE5wH2+CVQAuwErgCl4nmznJEshAx59cwpXPfxvTrrxEZau2ZBudQwjK3j+pWm88dZs3nhrVlruH0IC7NOJP73mXurvBU7A265yjogMiWh2EbBOVfvjOXvd7voOAUYDQ4HjgftEpDCBzNuBsao6AFjnZCMiA4DrgENVdShwZRMezU4ENTrd3OdiYF93fivwHXeeuu3VecDefbtTXFRIcWEhRealZhiBGHPeoZw3+mBOPTk9EbLrtIDaBEeAnDrDgUpVXaCqNcB4YFREm1HAOHf+HDBCPJ/4UcB4Va1W1YVApZMXVabrc7STgZN5qjv/EXCvqq4DUNWVjXooUQi6pybsINAOz+goMBfY5MrtlzGJHDy4Dx/dcTm1dSHKSmzbk2EEoXevjvxwzOFpu3+gqAPe9NoFInK+r/RBVX3QnffA264Spgo4MELKN21UtU5ENuDtm+yBt5/S37eHO48msyOwXlXrorTfE0BEPgQKgd+r6n/if7lgBDUWC9znh8CFeI/uE6CnK/86GcoYOygsKEi6wXln7pcsWp2e+e5848nPZ/Dtf/6NKV9beMJ8QSFh5lC3pDNOVSt8x4M+MdGGQpErQbHaJKscvAHJAOBIvOWTh0SkPEr7BhPU6PwdT8E9gDLgFVVdizc0A29tx8hgJn/5FZc9MYHLHn8pcWOjyUxfuYyqzRuYsyb2rMTm2mqumPQCD8+bnELNjOaiLlRIbYIjgMt0FeDfvNcTWBarjYgU4c1ArY3TN1b5aqDcyYi8VxXwkqrWuqm6eXhGqMkEMjqqeh+etbsHuJodjgNrgBtdudEEQiFl/vLVaDPFydijSwcO6NuD4/bes1nkGztz40EjeeLYszh/8P4x28xat4JXlszhkTQYneVbNrKxpjrl981lEnmuJYpW4JgCDHBeZSV4jgETItpMAC5w52cAb6v3wzEBGO282/rhGYmPY8l0fd5xMnAyw2+lLwJHAYhIJ7zptvCMV5OIF/Dzl8DLqvoFgKo+Azzjb6OqzwPPJ0ORfOf+1yfxtzcmc+V3vs0PRnwr6fI7tWnFExef1ai+y9dv4qwHnmafnt2497zINU0jGm1KSjmsR9+4bQ7s3JtbK06kf7vOqVHKsWjjOo596WF6tm7H26f9KKX3TjaTpnzJrX99jYu/fxgnH7dv4g7NSAgSOgokGum4NZrLgdfx1lIeUdXZInITMFVVJwAPA0+ISCXeCGe06ztbRJ4F5gB1wGXhDfvRZLpbXgOMF5GbgRlONq7tsSIyB29N/2pVXRP4YcQh3qLBn4A/isgC4GXgFeA936JT0hGRXsDjeN5yIbwFtjtFpAOewesLLALOCntV5Apd27WmqKCAzm1bpVuVXdhcXc36rdtYsnZ9ulXJKUSEs/eIPRJqLloUFdG6uITOLTLv31pDWVS1hg2btvHlolXpVoV6LaQu1PSAn6r6KvBqRNlvfefbgTNj9L0FuCWITFe+AM+7LbJcgV+4I6lIrOkc58d9MnAScDDeVNwmvL06LwOvuZhryVNGpDvQXVWni0gbYBqeC98YYK2q3uY2NrVX1WviyaqoqNCpU6cmU71mR1UzNhrw4jXrKW9ZRrsWZelWxUgCIVUKMvTfWkMIhZRZny9l0IBulBQ33vFGRKapakUT+j86ctyZY8r7d4zb7t2fvcyqaUvbq2revsHFi702B2+YdruItAdOxNuXcyzwXSAkIlPwDNC/VfXTpiqjqsuB5e58k4jMxXPhG4XnRQGeL/l/8YaFOUWmGhyAPh2T4riSdjZUbweBdiX5bTxzweAAFBQI+wzpmbhhClDnoWbEJ6gjwTpVfUpVzwE64y0wjcXzmrgZmC4iS5KpmAvtsD8wGejqDFLYMHVJ5r2M/GBzbTVHvvQ3Rrz0INvqatOtTl4RCimTKpewtSZ3n3tIEzsTWOy1RmzqVNWQqr6rqler6hCgP96839xkKSUirfEcFK5U1cAByETkYhGZKiJTV61K/xyvkVkUSgEti4ppUVScM2/62cJTk2Zw0SPPc+vL76RblWajXguoC8U/NJgHW07T5N2HbiHqTnc0GREpxjM4T6nqv1zx1yLSXVWXu3WfqJsf3CarB8Fb00mGPkbu0KKomHdP/QmAhRdKMXt260zXtq3Zu2e3xI2zlAa4Rec1gYyOS13wQzx/7l5AZH5lVdU9mqqMiwX0MDBXVe/wVYX90m9jZ19yw2gQqTA22+q3MX9TJXu1G0qBZf0A4MA9evHOtdntnp2IIGs69iYcfKRzG/BLdx4kTENjORQv7fVnIvKJK7ve3f9ZEbkIWEIMd0HDyASeWPQ0H675H9/rPZpjux2TbnWMFFEXEupC8V8ybHotuNE5D8/YbAY+A2qaQxlV/YDYsb9HNMc9DSPZ7NmmP/M2fUGfVn3Scv+tdduoCdVSXtI2LffPV9Sm1wIR1OgU4Y1m9lLVpHqpGUaucWSXIziyyxFpu/91n93K2pr1jN3vRjqVdkibHvlGkDUdm14L7r32IN4IxAJ3GVnL9voafvjxbVw69S/Uh+oTd8hSWhW2pLSglCKJvzveSC6KoBr/SJTELR8IOtL5LbA38JrLr7AEL7ZPGFXVi5KtXK5QXVvHJX99nhalxdx9+WkUFNg/vHRQE6rl6+3rKCoopE7rKSQ3f5Rv3vsaQqoUFeTm98tU6lWoS5BPp5ni+WYVQY3O+ezIEnpYjDZmdGKweVs1sxatoLiokNr6ekoLLDFbOmhb3IqHh19LgQilhSWJO2QpBVKAvdeknh2jGSMeQX/9bsLGhY2mY9tWPHnduZQUFVLahPhQucAVk15g9fYtPHTY2bQoKk7ZfbfXV/OneXfRrrgNPxvwk5Td18gfgq3p2M9o0F/A9nhrYKcAE1XVEnE0kD17pjZ8fSZSHwrx9rL5bK+vY031FnoWpS6e26a6zczbNJ/SghJCGkrJ/plV2z5mW91yerexdBD5QEgLqE/oMp0iZTKYoEbnWbw01XPM4BiNpbCggOdGjGFzXQ09W6U2gGjn0o7cOPQ6ygrLUrZhc/KKK6jXbbQtGUi7kj3ZWP0pbUqHUiC5O7WXz6iaUQlCUKMTTjHwXxF5lF0dCVDVx5Osm5EGFm9ax8KN6ziyx+7NIn9gefpite7Rul9K7zeg/AdsqV1Mm+J+LN7wEAvW/4VebS9kQIdrU6qHkRpCFmU6EEGNzr1402vtgd9EqVe85GtGlnPhW/9kwca1PHPcuRzYtXe61Wk062rWMnfjLL7V4WCKC1K3duRnYPsdYV9aFe9OUUFbWhU3OVpUTJZueZe60Db6tDm+2e5hxCYUkgDTa2aUGrKqbU8rD9itVRs211bTu3V2588Z/9XjfLp+GtWh7RzV5dh0q0PnViPp3Gpks8mvC23nw+W/QgnRqWw/WhXnbmDNTMWm14IR1Ogc1axaGBnDtNVV3kL/9q10b5W9YVQOaD+czbWb2LPN4HSrkhKKCsoY1P771Ia20KLInFbSQRCX6SA2SUSOx4vaXwg8pKq3RdSX4s0sHQCsAc5W1UWu7jq87Sv1wM9U9fV4MkWkHzAe6ABMB85X1RoRGQP8CVjqbnuPqj4UQP2EBDU6n6vq17EqRcTcc3KEsYeewtItGxjaoWu6VWkSwzscwvAOh6RbjZSyT8fL0q1CntP0fToiUoi3nHEMUAVMEZEJLpNzmIuAdaraX0RGA7cDZ4vIEGA0MBTYDZgoIuEoMrFk3g6MVdXxIvKAk32/6/OMql7epC8UhaBuPO+6PDa7ICKXAs8lTyUjnRzfeyAXDR6e0amzDSMTqQ9JwiOAURoOVKrqAlWtwRuFRL7UjwLGufPngBEuLcwoYLyqVqvqQqDSyYsq0/U5mh2/3+PwHMaalaBGZ0/gPRHp5S8UkT8A9zRAjmGklZrQNl5b+numrH6iQf2WrL+XGctOZ3vd0sSNjbzEW9NJFHstIT2Ar3zXVa4sahtVrQM2AB3j9I1V3hFY72REu9d3RWSmiDwX+dvfFIIai9XAHniGp5+IFInIE8CvXP30ZClkJI/nvpjFMc8/yvSVy9KtSsawtnoRX256j5nrXmhQvzVb32JzzUy21MxrJs2Sx/baL1m/9bV0q5F3JDI4PqNzgYhM9R0X+8QEyVcWq02yygFeBvqq6j7ARHaMrJpM0DWdbwNvAH2A9/GGbYe7ugnAuclSyEge7y1dxBfrVjPt66UM67JbutXJCLq1GMwx3a+nXUnDnsfgLnezpWYeHVse3UyaJY9Fq35Edd0XbG/7U7qVX5NudfKGIPl0nHfbOFW9M0aTKrzszGF6ApFvjeE2VSJSBLQD1iboG618NVAuIkVutPNNe1Vd42v/d7y1n6QQaKSjql8ABwOz8BaowgZnrKqeqqpbk6WQkTxuOngED4wYxZihw9Ktyi58tv5L/vL5eFZtX5/yew9sN5JuLYY0qE9ZUY+sMDgA5S1PpEBg1aa7qan7KnEHI3logiMxU4ABbkapBM8xYEJEmwnABe78DOBtVVVXPlpESp1X2gDg41gyXZ93nAyczJcAItbwTwHmBtI+ADFHOiLy/SjFjwA3Ay3wRjyfhttZRILMo7ysBSf0y8wUSM9+9TYfrZlDr5ZdOKt3dvyYZwKh2vlIQXuksFPMNt3KrwJqqA9tpLgwqv+P0QwEW7dJNBLSOhG5HHgdz735EVWdLSI3AVNVdQLwMPCEiFTijXBGu76zReRZYA5exJjLVLUeIJpMd8trgPEicjMww8kG+JmInOLkrAXGBHwMCRGNsZtJREIET3SnqppR4ZMrKip06tSp6VYjJazauoV3Fi3g5AGDaFGcnt33DaVy01L+u3I6Z/U+mrbFrdKtTlYQqv2c6tUnIYV9KOvydtLl14a28MGK31Fesgf7d7ok6fIzGRGZpqoVTej/aM8/XjamtE/8TbnL/u9Rts1a0F5VUz/EzxASGQrzm80CbvngXV76Yi6rt27l0ooD061OIPq36UH/NpFOOUY8vBFON6R4ULPI31jzFVVbPmDV9ll5Z3SSgYYKUIsynZB4RseiEGQJJ+wxgKqNGzisd990q2I0EdUQuuHnoLVI+V1468Rh6igpGoiUNs90ZMeyQRzR/Q+0KrIpucZiRiUxMY2Oqr6bSkWMxnPcHgM4bo8B6VbDSAa6FbZPBEKgG0E67Kir+RhqPkB1G9LyjJgimkLv1vau2WiCOwvkFCJSC2wn+rdvpao75U2P50jQTlU3NPDmDe5jNI05q1by6KfT+cmwb7FHh47pVsdoIlLQGjo8DdQhBR12riz7DqK1UPKttOjWWNbXbOGuLyZwSKfBjOy2X7rVaTY0BBpKsCKRm0bpM1WN6iIrIrvs4Yw3AblERP4iIvsmuqOI7CciY4HFwfU0ksHjMz/hubmz+cfsz3Yqf33RfC5+4wW+2mTvAImYsuoOJiw+ly21K5Mqt76uik0bbqW+rqpB/aRkX6TkgF3LpRhpeRZSlNq8QE1lytoveGPFDJ5enNuTJ8E2h+bkMnm8rIS71MVb02kDXAlcKSKLgEnAbLyopuCFUNgLb/9OH7ynGWq4vkZTuKRiOOVlZXx/n53fIJ+c8wnvL13EoT36cEEG7tPJJFZsncaGmgVsrltGq+LkJZnbsvk+tm55DNXNtC2/NWlys40juuzNJf03sH/75kkMmDHk6fQaUCsinVV1lb9QRDoBtZGN4xmdA4A/AMcC/YC+MdqFTffrwPUN1dZoGn3alXPtoYfvUv67Q47mnSULOGvg3mnQKrs4erexbK5bRtcWyZ36adHqfEKhzbRodX5S5WYbJQVFfK/vkelWIwUIJNqnk5tG6XHgSRG5RFUXwDcpEx4kSnLPeI4EM4DjRWQwcCEwAtgHb3MRePkaZgJv4W02+jyZ38JoGv3LO9K/3NZ4Qhpi+rp32K1FP3ZrEf1Nu1Vxl6SOcMIUFw+mvMNdSZdrZCh5OtJR1bEi0gaY7FIzhJ/EPcBfI9sn3NCpqnNxgT1FpABvWg1gjaradJqR0XyxaTovVN1L59IeXDnw7nSr0yTqt79FqGYaRW2uwMvjZWQUGmCkk6Oo6k3ATSLS2V2vitW2QVEEnJGJKcwwMo1eLfdkcNvh7N46+6cZazfejNYvoqBkfwrLjkm3OkYEQdJV5+I+HpeX54fAcXgjnDdE5OFYgxLLg2M0iimrFzFzXcO8stJBq6K2nNf3Wg7pdFK6VWkyxW2vo6jVDyko/XaTZVXXfM6GLc8TKwyW0QhUIJTgyE3vtZuBE4G/AXsDnYHbYjXOqHhpRnawcttGLvxgHIUFBXx80rWUFmZHvLdsp7DsWArLjk2KrGVrLqGmbh6q9bRteTIFBS2SIjefEfWOuG1So0qqORkY5oKVblPVW0VkcqzGNtIxGkx5aUuO7LYnx+82lJICe2/JRtq1OpvSkmEsXnsVC1f9IOny60J1iRvlGonSGuSuo4H4so/i0ifEXHQ0o2M0mJKCIu456BxurzgdbzrXCLOuei7/XjSCWWsy22mhQ9uf0LX9nyiQMgoLy5Mqe/KaD7jikwt56+vkZy9dt30bk5dl6LRuiADTaznJShEJx+FqC3wI3Bursb2mGkYS2Vq7jOr6tWysmZ9uVRLSomQwe/WcjUhyp0e31m9BUb7YVMkbX1/PeX3OZe92eyVF9hVvvsp7SxZx33Enc2L/DMsVlbsjmUScireFBuDHwHxVjRmdJl7std4NuauqLmlIe8PIJmpC1WyuXU+H0q5x2/VoPYKjip6kdUmfFGnWNJJtcACO6nIcQ9vuy39WvMWK7Z8wf1Nl0ozOsK7dWbh+Hbu3b58UeUknkdHJTaM0DAjPetQCfUWkb6yg0fFGOotoQBK3BLKSgogcD9yJt0H1IVWN6SFh7GDpxo10btWKksLCxI2NqIxb+EcqN3/GT/a4kX6tB8dt275saNLuW1P3NZ+v+gntyg6hT/urkya3uelS1o2ze53BXu2GNMrg1IbqeX3JPA7s2pvOLVoDEFKlXkOM2Wd/BnXsnGyVm4yEBMndKbR4/NJ3XgoMBz4hRnqcRGs60oCjWXE7Xe8FTgCGAOeISMMS3ech7y9exOGPPsSV/3k13apkNW2KyykqKKG0MLVeXtvqFrC55lPWbpuY0vsmg9LCUoa135/igoaPpv4x/xN++sFL/Prj178pW7JxPXdP+4ibP3yH2vr6OL3TRH46EaCqp/iO4/B+n1fHah9vdHJj0rVrGsOBSl9sn/HAKLx84EYMWhaXUFxYSHlZWdJlb6mtoWVRcSBngn8t/oRpq7/iun2OpXVx9u2mH937p5yll1IgqVzL79UAACAASURBVB0ttis7mMFdHqVFEyJLr9s+g8UbxtG//eW0LumfRO12Zem2KmZvmMWRXY6mpCBe8OH4fKtLL/bp0I2RPXfkierbrj2/PfQo2pe1oNhG7ZnMCrxg0FGJF3st04xOD+Ar33UVsFNuZhG5GLgYoHfvBi1JpZTn3pvJa5M/5/djjqVX5+R6DkVywG67MfOSy5M+tfb+8oVc+PYznDtgf24aflzC9vd9/j5fbVnHsT0GcUS37Ew4l2qDE6Z9iyOa1L9q0z/5eusbtCzuw54dfp4kraLz1OIn+GLzPMoKW3B458brPbh9F1468cJdyn+w764pHzKFIPt0gox2Ei0jiBcD6XG8oMxrgLNVdZGruw64CG9h/2eq+no8mS4w53igAzAdOF9Va3z3OgP4J/AtVZ0aQ9+72DHbVQjsB0RtC41YhxGRLsAur80pcCSI9jq9059QVR/Ei2xKRUVFxg5m35gyjxmVS/n0y2XNbnSAZlnLqamvI6TK9vpg+zFuGXYyn61bxqFd9ki6LkZ89ii/lJZFvejV9uxmv9fRXUbQuqg1e7XL/rBDDUab7hbtW0Y4Bu/FeoqITFBV/4zORcA6Ve0vIqOB24Gz3XLDaGAosBswUUTCLn6xZN4OjFXV8SLygJN9v9OlDfAzIOZGT4ffwNQBT6rq/2I1DmR03IO4Gc8drl2UJqlwJKgCevmuewLLmvmezcLvLjiWTyqXcdy3BqZblUYzoucAPjr9p3Qsaxmo/YGd+3Jg577Nq5QRlZbFPdmj/SUpuVdFh+FUdBieVJkL1q/lH3NmcuE+w9itddukyk4qyVm3CbKMMAr4vTt/DrjHxT8bBYxX1WpgoYhUOnlEkykic4GjgXNdm3FO7v3u+v+APwJXJdD5DWCjqm51tiLuHyno5tCrgWuActLkSABMAQaISD+343U0MCEF9006PTq146SDBlNUmN17c7u0bE1hQeq/w1db1vHG0s8tbliecP+Mj/n7p1N5bOYumY8zC90xxRbrcEbpAhGZ6jsu9kmJtozQI+JO37RxkQA24EX/j9U3VnlHYL0vmsA39xKR/YFeqvpKgG/+EtBKPP/7qcAkEYnpWRz0F+M8vMf1b3etwF+AlUAlKXA6cA/mcrxkcXOBZ1V1dnPf18g8fj75X1z+0T+ZuGxeulVJOzPXv8W4hb9i1fbc3SY3Zu9hnDVoL84Zsm+6VYlPKMDhMU5VK3zHgz4pCZcR4rRJSrlLYTOWnV2h41HmUhmMAKaq6iC8eGxRCTolFs5+dRGeZwKqerWIPAN8DGwKKKdJqOqrgPn+ZiH3zpzEB8sXcefhJ9PF7btoLMf2GISIMLi8W3KUy2LmbviAqq1zWLz1MzqXZa7zTFMY2qkLfzzq+HSrkRAhKQE/gywjhNtUiUgR3pLH2gR9o5WvBspFpMi91IfL2+B5n/3XeaZ2AyaIyClxnAk6AOcAL7qiXdJUhwk60gk/ytVhYSLSEQhnC70yoBwjT3lxwWwmrVjC52ubno7p4oGH8s+jfkDPVs3vhJHpnLDbpZzc4+fsX56c6NNGEwgncUt0xCfIMsIE4AJ3fgbwtnpzzROA0SJS6rzSBuANCqLKdH3ecTJwMl9S1Q2q2klV+6pqX+AjIKbBAf4EfIFntF4RkbZ46zxRCTrSWe0ElgNLgT7A00C1q8/q//vr60KsWrWRbt2z+mtkNPcfdRpfrFvF4T0av9/E2JW2xZ3Zp3xEutUwICmOBC49QHgZoRB4RFVni8hNeFNXE4CHgSeco8BaPCOCa/csntNBHXCZqtYDRJPpbnkNMF5EbgZmONkN1flJ4ElfUS0u23Q0ghqdWXhGZzDe9NalwMjwPYGoMXayhXv+8hr/fnE61914Gkcdm5wYUQZsqqmmpLCQ0sIi+rfrSP92HRN3MowsRULeEZcARinaMoKq/tZ3vh04M0bfW4Bbgsh05QvY4eEWS58j49WLyO9i9LtRRH6sqn/zlwc1OjfjDd3WATcAvfEyxRUA7wOp8cdsJlq3KaOgUGjZKvt2ymcqyzZt5NinH6NX23JeO+f76VbHMJqdIGs6OUq0Nf3wPOKWyIpARkdVJwGTfEWniEgZUKyqKXEiaE4uunQE5//wCEpKLNNDMhEES7eT2SzasoLHFv6HM3sdwdB2NvXZJHI4vlo8VPWOyDIR+a6rezKyrtG/sm6It72x/TMNMzjJZbc2bZl04Y8tsnWGM3HFNN5fNZMWhaVmdJpIsqbXsg0ROR0Yg+f1FqbCrSM9pqrj/O2DRiRIFNJVVdV+tY2daF3S+ICPRvOzvb6WY7oNp6ywhJFdMzemmZHx3Ar8BNjorhXP0ewqPMeznQhqKGySxDByjO//70GWbF7Dc0dcTrcWHdKtTvaTp9NrwFZV/a+/QES2qeq0aI2DGp1xEdeFQF/gUGAb8GzDdDQMI92UFhRRVFBAQeDtekZcNPH0Wo6+vR8cpezXsRoHdSTYNcY4ICJHAxPxssQZhtEE5qxdye8mv8mYwQdwUt9BzX6/Rw/+IXUaoqww+Smr85IgI50cGgmJyFQXxqfaXXfHW9s5Hy9wQNToMU16xVHVt4HNeOGvDSMwIQvWuQsfLFvIlJVVvLro88SNk0BRQaEZnCSSKNhnDrpTbxaRu0XkbBGZALwMbAUOU9XTY3UK6kgQbaNFGXAs0LRAWkbesWjNOs78+z84aPde3H1WzLiAecf3Bw2jvLQFR1jUhuwk/9Z0jsELBn0FXmy3W4B/qOqGeJ2Cruk8RuzHqXixeQwjEJura9hSU8PKjZvTrUpGUVZUzFkD9km3GkYjkQBrOrlklFS1FngUeFREBuFNrU0XkY+AR1V1YrR+DXFzjrUG9jEuRbRhBGGv3bryxk8vpLxli3SrYuQYW7ZW87s/v0y/3p24bMyRqb15/o10vkFVPweuFZHrgRPwotQ0yegcFaWsGqhS1apGaZmFzJz5Fdu21XDggZZyuan0bB8tAa2RCUytWsqzs2Zx5SEHs1vbDM7UGYWlK9YzecYiPv/y65QbnRxdt2kQqhrCy7v271htgnqvZXVAz2RQU1PH1b98mvr6EE/941K6drUfzVygrj6U9Rlck81D06byZuWX9G1fzqUHHphudRrEnrt35bbrT6NLpzaJGyebnRO1RSfPjRLEMToxnAdioqqPN12dzKWkpIgTTtyXzZu307Gj+U7kAjMWLuWiB57n9AP34tenH51udZpESJUZVcvYu3tXSoqaFhzk54ccwh4dOnD23nsnSbvUcui30jMTEWSkk6P7dBpEvH+djxHcLiuQ00YH4MqfZ372QiM4m7bXUFtfz9rNW9OtSpN5fPIM/vDmu5xbsS+/O6FpBnRg585c3blzkjTLM2wkk5BEr0RmmB2qyv/97gU2b97OzX84i5JSCzWX7Rw+uB//+Ok57NE9+0PA9OvYnvIWZezZxXIWpQtb0wlGvF9Ov/NAO+ABYCXwV7w82rvhpaneDfhxcymYKdTW1jP5o0pq6+rZuHEbnTqnYc7YSCofzV3MZfe8wPEVA7nlwhPSrU6TOGJAPyZfldVprbIfW9MJREyj43ceEJHHgK7AcL+3moi8BSwGTgdeaD41009JSRF33z+G7dtrzeDkCEWFBYhAcZGlXzCSg00NJSao285p7jNyY0VL93lKctTJbHbfowtDhvZItxpGkqjYsxfv/vlSfnfeMelWxcgFNMARABE5XkTmiUiliFwbpb5URJ5x9ZNFpK+v7jpXPk9EjkskU0T6ORnzncwSV/4TEflMRD4RkQ9EZEjDH0h0ghqdavf5hojcICJjROQG4D8R9UYesWrtJt79eD6hUPbOGbQqK0EsvWnG8dZ7c3nynx9l1b+tcESCuEci7zaRQuBevA2WQ4BzovzgXwSsU9X+wFjgdtd3CDAaGAocD9wnIoUJZN4OjFXVAcA6JxvgaVXdW1X3A/4I7JIdtLEENTqP4I0cewM3AQ+7zz6u/uFkKWRkDzfd8x+u+/MEJv4vNQEqjfzh9jv/w98ff5/KBSvTrUpwkjPSGQ5UquoCVa0BxgOjItqMYke6meeAEeK9OY0CxqtqtaouBCqdvKgyXZ+jnQyczFMBVHUjO2gVWPsABHXBuh6owXMc8G9S2YxnaW9MlkJG9vDtA3Zn85btDNy9a7pVMXKMK348gmUrNrBHv+xx3W6A99oFInK+7/pBVX3QnfcAvvLVVQGRO3S/aaOqdSKyAejoyj+K6BteD4gmsyOwXlXrorRHRC4DfgGU4BmnpBA0IkEI+K2I/AXYB+gOLAdmJoooauQuZ590AGefZGmOjeRz0rHZF/g0PIUWF88ojVPVO2OJidkrcZtY5dFmtOK1905U7wXuFZFz8ZKyXRBN4YbSoM0mzsC8n4wbG4Zh5BxNn4SqwksTEKYn3haVaG2qRKQIb0vL2gR9o5WvBspFpMiNdqLdC7zpuPsb9W2iEC8MztuAquoIdx4PVdURyVLKMAwj20jS5tApwAAR6QcsxXMMODeizQS8Ucck4AzgbVVVl0jtaRG5A2//5AC8LAASTabr846TMd7JfAlARAao6nx3v5OA+SSJeCOdI9lht/3nkUicOsMwjPwgCemq3RrN5cDrQCHwiKrOFpGbgKmqOgHPcesJEanEG+GMdn1ni8izwBygDrhMVesBosl0t7wGGC8iNwMz2OEUdrmIjARq8bzakjK1BvGNzhJ27K9dghkWo5mpqanjxZemse8+vRk4sHu61TGMhhFgTSeIc76qvgq8GlH2W9/5duDMGH1vwcvgmVCmK1+A590WWX5FAFUbRbyIBH2jnRtGc/H+B1/wwIPvsOeAbjxw35h0q2MYDSOPk7g1hEZHrRSRAqCNea8ZyeKAYX05+qghHHyQJckzsg9BEU04f5YaZTKYQEZHRE4BRgIfqOqzInIRcCfQwuXDPk1Vs2gXl5GJlJe35NfX50VEJSMHCeQybQSOSHAJcBkgItIauAsv7poABwG/bxbtDMMwsoUkxV7LdYIanXAKwQ+Bg/ECf07C890W4LgY/QIjIn8Skc9FZKaIvCAi5b66qEHsDMMwMoWwy3S8wwxPcKPTyX2uwAsmp3j5da525clwNXoT2EtV9wG+AK6D2EHsknA/wzCM5GIjnYQENTrh4G/7AIe58/nsWBPa3lRFVPUNXwygj/B2x0LsIHaGYRgZQ8II0yHLtwPBjc4s9zkFLwrpNuAToJ8rr4rWqQn8AHjNnUcLgGdJbQzDyChsei0YQV2mb8VzGChz13eo6nYROdldfxhEiIhMBLpFqbpBVcPhF27A2037VLhblPZR/3QicjFwMUDv3r2DqBSVNSs3UtaihFZtyhI3NgzDAM8d2lyiExI0yvREERkEVACLVXWaq3oRLwBoZUA5I+PVi8gFwHeAEarf/PWCBMALy38QeBCgoqKiUX/9FUvX8aPT7qZT17Y8+vKVjRFhGEY+ooGjTOc1gTeHquoSYImIlIlId1VdrqqzEnYMiIgcjxcH6AhV3eqrihXErlkoKS2iZatSytu3aq5bGIaRg9g+nWAENjoicjDwF3Ys4heJyF14Sd3uSIIBugcoBd506YM/UtWfxAti1xx06NSGpydeTWFh0OWuhjNtUiXV22s55KjBzXYPwzBSjK3ZBCJoRIJhwNt4GeT8UaXr8aKPrgWuaooiLt93rLqoQeyai+Y0ONXba/n1FU8Rqg/x5Ku/oHPXds12L8MwUoeoIqH4VicJqQ+ynqAjnZvwRiFzAf/r+WPAFXghcowAlJYVc8pZw9m2tZoOHVsn7mAYRlaQpHw6OU9Qo3Mo3ujmRGChr3yu+2y8q1gecslVJ6RbBcMwko1NrwUiqNEJ+w4vjyjvFFFvGIaRlwQa6ZhRCrw5NDy6OStcICLtgDvc5YJkKmUYhpFteN5rGvcwghud8XgOBOPYYatX42WvU1dvGEaWUVdbT1Xl1+lWI3ew2GsJCWp0bgP+i2d4wkeh+3wP+GNzKGcYRvNy7/XP8KMjbuat55pt69tObFy7mZ+f9CfuvTYH31MDhMExR4OARkdVa4BjgO8D/wAm4o1uxgDHuHrDMLKMDl3bUVRcSLsUeVKurFrL59MW8tEbM1Nyv5RSr4mPAMMdETnepXGpFJFro9SXisgzrn6yiPT11UVNAxNLpoj0czLmO5klrvwXIjLHpZp5S0T6NO3h7KAhEQnqgSfdsRMi0s7SVhtG9nH+VSfxvV+cQEFB8+1N89N/n97c/vyVdOxenrhxliE03ZHApW25F+8lvwqYIiITVHWOr9lFwDpV7S8io4HbgbMj0sDsBkwUkT1dn1gybwfGqup4EXnAyb4fmAFUqOpWEbkEbzbr7CDPIRFN+pcmIm1E5PfAomQoYxhG6kmVwQmzz6F70mP3Lim9Z0oIB/yMdyRmOFCpqgvcDNJ4vPQufkbhra8DPAeMEC+MS6w0MFFluj5HOxk4mad6X0Xf8YUj86eaaTJx/7WJyJki8qbL6PmeiPzUV/dzPGPzG6BtshQyDMPIRoLk03FcICJTfcfFPjFBUrl808blINsAdIzTN1Z5R2C9L49ZrLQxF7Ej1UyTiTm95oZtT/mKBgCHurm91sCPwk3ZkeTNMAwjLxFVJNFoxqsep6p3xhITs1fiNrHKow0u4rXfcSOR8/CyCxwRpW2jiLem89MYSoXLBdgE3I0XCNQwDCN/UaDpUaaDpHIJt6kSkSKgHV78y3h9o5WvBspFpMiNdna6l4iMBG7Ai/xf3cTv9Q3xptf2wnuM1+FNn7UDrgeK8dylXwT6qeqvVXVdshQyDMPIRsIjnURHAqYAA5xXWQmeY8CEiDYT8AItA5wBvO3yj00ARjvvtn7sSAMTVabr846TgZMZTqa5P/A34BRVXdnohxKFeEanjfu8Q1U3q+omdkQgAPihqq5NpjKGYRhZSwgIafwjAW7EcTnwOl5sy2ddepebROQU1+xhoKOIVAK/AK51fWcD4TQw/8GlgYkl08m6BviFk9XRyQb4E94yyj9F5BMRiTR8jSaIy/RBzsshkqH+clV9L1lK5TKL56+gvGMb2nWwJHGGkVMkafOnqr4KvBpR9lvf+Xa8aDDR+kZNAxNNpitfwI4caf7yZsscEMTo/DfiWqOUa0BZeU3l7KX87LS76N2/Cw+8+st0q2MYRlIJ4BZtEQkSGopoIxyjkbTr0IpO3drRd89u6VbFMIwkI/WK1JtVSUQ8ozMuTl3OsnrFejp1a57d0p27l/P4e9c3i2zDMNJMkCg3ZpNiGx1VvTCVimQCz97/Fo/+8d9c/JtRnPaDpLmlG4aRBwTxThOzOrYO46e41HscJaXFadbEMIzsI3Com7zGjI6P035wBMeddSAtW1siVMMwGkagNR2zSWZ0IjGDYxhGo1BspBMAMzqGYRjJIFAkaTNKZnQMwzCSgIRsei0IZnQMwzCSQfCcOXmNGR3DMIxkEMTomE0yo2MYhpEUQkDC6TWzOmZ0DMMwkkDA1AV5jxkdwzCMZGBrOoEwo2MYhpEMNEDOHLNJZnQMwzCSQijkHXExq2NGxzAMIxlYRIJAmNExDMNIBja9FoiCdCtgGIaRE4TqEx8BrI6IHC8i80SkUkSujVJfKiLPuPrJItLXV3edK58nIsclkiki/ZyM+U5miSs/XESmi0idiJzRtAezM2Z0DMMwkkFIEx+J8u2IFAL3AicAQ4BzRGRIRLOLgHWq2h8YC9zu+g4BRgNDgeOB+0SkMIHM24GxqjoAWOdkAywBxgBPN/p5xCDjjI6IXCUiKiKd3LWIyF3OQs8UkWHp1tEwDGMXwms68Y7EDAcqVXWBqtYA44FREW1GsSOz83PACBERVz5eVatVdSFQ6eRFlen6HO1k4GSeCqCqi1R1Jt6W16SSUUZHRHoBx+BZ2TAnAAPccTFwfxpUMwzDSEACg7PD6FwgIlN9x8U+IT2Ar3zXVa6MaG1UtQ7YAHSM0zdWeUdgvZMR615JJ9McCcYCvwJe8pWNAh5XVQU+EpFyEemuqsvToqFhGEY06uuhoD5+G8/wjFPVO2O0kGi9AraJVR5tcBGvfbOSMSMdETkFWKqqn0ZUBbH8YRkXh98eVq1a1UyaGoZhRCE502tVQC/fdU9gWaw2IlIEtAPWxukbq3w1UO5kxLpX0kmp0RGRiSIyK8oxCrgB+G20blHKov71VPVBVa1Q1YrOnTsnU3XDMIwEBJheS2x3pgADnFdZCZ5jwISINhOAC9z5GcDbbiZoAjDaebf1w1uS+DiWTNfnHScDJ9M/y9QspHR6TVVHRisXkb2BfsCn3toWPYHpIjKcYJbfMAwjrWh9PSqBptfiVGudiFwOvA4UAo+o6mwRuQmYqqoTgIeBJ0SkEm+EM9r1nS0izwJzgDrgMlWtB4gm093yGmC8iNwMzHCyEZFvAS8A7YGTReRGVR3akOcRi4xY01HVz4Au4WsRWQRUqOpqEZkAXC4i44EDgQ22nmMYRsYRdotuIqr6KvBqRNlvfefbgTNj9L0FuCWITFe+AM+7LbJ8Ct4LftLJCKOTgFeBE/Hc/7YCF6ZXHcMwjBgkXLexkAQZaXRUta/vXIHL0qeNYRhGAOrrgUTTaynRJKPJSKNjGIaRbagqmijKtAUENaNjGIaRFCyJWyDM6BiGYSSDkGLzZ4kxo2MYhpEENFSP81CO08iMkhkdwzCMZKBK4viYZnTM6BiGYSSBUCiESnyjErKRTubEXssEJr86nelvfZZuNQzDyD7eW6lfgYZiHttDmwlRB15U6LzFjI5jzfJ1/Obk27ju+JvZtmV7utUxDCO7eHI1y6nR6pgNFjGPvgwK7z3MW2x6zVHeuS0jzz+c0hYllLUsTbc6hmFkEapaO1S+xRLm05+9dqnfrlvZyFqq+LIwDeplFGZ0HIVFhfzqscvTrYZhGFnKHKaWtKZdTW8dQIns/OIaHuV8qv9LeibObMOm1wzDMJKAqtb2ZgBLmL9TeXiUM5NJeT/KATM6hmEYSWMOU0si13Z8azl5P8oBMzqGYRhJI3K0Y6OcXTGjYxiGkUT8ox0b5eyKORIYhmEkkbAn23xmsoWN5rEWgY10DMMwkswcppZsYI2NcqIgubpPSURWAYsjijsBq9OgTiZhz8DDnoM9A9jxDPqoaud0K5MP5KzRiYaITFXVinTrkU7sGXjYc7BnAPYM0oFNrxmGYRgpw4yOYRiGkTLyzeg8mG4FMgB7Bh72HOwZgD2DlJNXazqGYRhGesm3kY5hGIaRRszoGIZhGCkjZ42OiPxJRD4XkZki8oKIlPvqrhORShGZJyLH+cqPd2WVInJtejRPHiJypojMFpGQiFRE1OXFM4gk179fGBF5RERWisgsX1kHEXlTROa7z/auXETkLvdMZorIsPRpnjxEpJeIvCMic93/B1e48rx6DhmHqubkARwLFLnz24Hb3fkQ4FOgFOgHfAkUuuNLYHegxLUZku7v0cRnMBgYCPwXqPCV580ziHgeOf39Ir7r4cAwYJav7I/Ate78Wt//EycCrwECHARMTrf+SXoG3YFh7rwN8IX7t59XzyHTjpwd6ajqG6pa5y4/Anq681HAeFWtVtWFQCUw3B2VqrpAVWuA8a5t1qKqc1V1XpSqvHkGEeT69/sGVX0PWBtRPAoY587HAaf6yh9Xj4+AchHpnhpNmw9VXa6q0935JmAu0IM8ew6ZRs4anQh+gPcGA94/uq98dVWuLFZ5LpKvzyDXv18iuqrqcvB+kIEurjznn4uI9AX2ByaTx88hE8jqKNMiMhHoFqXqBlV9ybW5AagDngp3i9JeiW6AM96fPMgziNYtSlnWPoMGEOt75zs5/VxEpDXwPHClqm4UifZ1vaZRynLmOWQKWW10VHVkvHoRuQD4DjBC3aQt3ttLL1+znsAydx6rPGNJ9AxikFPPoAHE+975wNci0l1Vl7tpo5WuPGefi4gU4xmcp1T1X644755DJpGz02sicjxwDXCKqm71VU0ARotIqYj0AwYAHwNTgAEi0k9ESoDRrm0ukq/PINe/XyImABe48wuAl3zl33feWwcBG8LTT9mMeEOah4G5qnqHryqvnkPGkW5PhuY68BbHvwI+cccDvrob8LyY5gEn+MpPxPNw+RJveirt36OJz+A0vLe3auBr4PV8ewZRnklOfz/f9/wHsByodf8GLgI6Am8B891nB9dWgHvdM/kMn6djNh/At/Gmx2b6fgdOzLfnkGmHhcExDMMwUkbOTq8ZhmEYmYcZHcMwDCNlmNExDMMwUoYZHcMwDCNlmNExDMMwUoYZHWMXRGSkiKjv+HYadKjy3f+hVN8/HiLyik+3I3zlaX9uzYGIfOS+T72I7J1ufYzsxoxOBiMir/l+wGpEpHOctvdE/ODtk0pdU4mIfOD7nhNTfO+jgZPc5SRVfTeV928oIlIgIsvdszq6kWL+4D4LgD8nSTUjTzGjk9k87DsvBr4XrZHbXX+Or2iqqs5sTsVSwP8BV7vj2TTr4ucPvvOxadMiOIfgxeZbAzTWQE7A2zAJcKx/dGcYDcWMTmYzAVjtux4To93JQAff9SPNpVCqUNW/qeqf3fFGuvUBEC8R3nB3uQl4OY3qBOV09/mSqtY3RoB6O8jH+4oubbJWRt5iRieDUS/nyxO+on1FZL8oTcf4zrfhhUD5BhFpLSJXuWmpdSJSKyIrxMuoOqIxuonId0TkRRFZ6qb+NorIDBG5RUS6xOhTKCLnici/RWSZiFQ7fWa66cEOvra7rOmIyM0iosChPrEjIqYVzxOR533XU6Lo0T+iz5kBv/bFvvOXVXV7wH6IyDEistV3z3fc36UoQpdfi5fd9EPXfpmIjBWRlk7OGSIyVUS2ub/hAyLSJs6tT3Ofz/t0GSxeZtFKJ6faTcF9LCL3ichRUeT4R5uniUinoN/dMHYi3XF47Ih/AEPx4keFj79G1HfFi68Vrn8yJI150AAAB3ZJREFUor4/Xhw6jXPcEtFnZET9t311BcDjCeStBIZHyCwHPkjQby9f+ypf+UOu7OYE/RU4DzgsomxYhC6/9tWtAkoC/i2W+PpdGqU+6nMDjgG2+spfAcpcXVFEn6lAKMr3eh24KsZ3fjOGvge4+g3h74iXOXNzgmf4UBRZAqz3tTk73f9v2JGdR1anNsgHVHW2iEwGDnRF54rI1apa666/x84pKr5ZBxKRQuBFYA9XtB5vFLQMLx1veEH8ehH5VFWDrJ1cB5zvu/4ML0pvd7yIvUVAZ2CCiOypqhtduyfZeYSyGG/6cB0wCG+KMBGvue9wGdDXlc0HHvS1maaqc0VkOl66ZoAfuyPMaN/5E+qNKOMiXhIwf9j7jwPoi4gcg/d8WriiZ4DzfX+/SA7Ae6Yv4gWnPMCVH+uOT/Ge2ynAvq5upIgcoKrTImSFp9Ze8X3HHwCt3PlavMyZq/HWfXbHS3O9C6qqIjINCDsjHOG+i2E0jHRbPTsSH3jTOv430VN9dZ/6yr8EL4irqzs1ot8BEXKf99VN9ZXHemMvxPuhCpd/jntjd/U/iuj3U1e+b0T5JKBlhC6dgba+611GOr46/4hpYoxn9n1fm01AG1c+JEKXIQH/BpHPpEuANjey8wjnIaAgok/kSGdlHF1XAK1c3d4RdZdE0Weuq/uur+weX597ovQpBnrHeAaP+fpGHV3ZYUeiw9Z0soN/4P14hRkDICLDAL9r9KOqqr7rwyLkTPWvH7DjTRhgWHjdIA5DgPZ+vXTndY1xeFNDYcIjm8i355t15xxHqOoq3TEqSgbj8dI5ALQGznXnfi+/Sao6J6C8SHf1dQH6/JYdI5y/Aj9S1VCc9gATVHWTO18UpW6LO58fUef/uyAig/FGkNvYkaod4D3f+aUiMk1EnhKRG0XkNLyXgSUxdFvjO4/pvm8Y8TCjkwW4H6F/+opOdHt2LvCVhfDeRP10IDiCl2ckHpHyVvgv1JvCWRulfWS/hQ3Qq1E4XR7wFYWn1872lTVk02nMHMcBqYx4IYjFUt955LSfP4tlXURd5P/L33Wf//EbePWmUG/DM0aCNwV5Lp6B/BewXER+EkM3iXFuGIGxNZ3s4WF2GJlivNHOub76N1S1KqKP3wCEgOuBeG6zGxLosDbiupv/wu0X8huYtRGfYfoBQUcYTeF+vDWoEmB/92M6wNVtpmH7f1ZFXLdnR5rjWMwFBrvzu0Vku6o+HK8DnlNILCINTTzCo9jnIytU9ToRuQVvXW8QnrPJ0XhTdi2crq+p6uKIrv6/baLvbhhRMaOTJajq+yLyBbCnK/o94J8Oi/Zj9gHwC3deACxT1SciG4nI7kD/ANNbc/CmlcJTOeeIyG2qWu2uL2DnN+4P3ad/SgfgBhF5W1W3+XToCNT4ppbi4f9hjjklqKpfi8h4vPUd2Hkz53hV3RzgXmEqI657kfiH93LgWjzvNQEeFJFqVX2yAfdtMM7pYX+8kdIrEXX9gPWqug6Y6A6cm3t4OrII2A/P2cOP35Ei8nkYRiDM6GQXj+BNjcDOP7ar8TyaInmZnd+2HxWR7wIz8EY8vfDedvfBM1pxN2Gqar2IjAVuckUDgY9FZAI7vNfCrMRb40FVPxWRf7PDW+5gYI6IvIRnxPbAc3o4BJgVTweHf0Q3XETuxEtNDnC3zwgC3MkOo1PmK29QPDdVXSgiVUBPV1QBRHqLRVKDt09mIt5zLgAecyOe5xpy/wYSHuW8paqRo9ezgZtE5H28VOUr8AziiRHtdhqdioiww5MOGh/dwMh30u3JYEfwA286y78nJ3yMjdMnyD6dnTzEiL9PpxB4OoGs1cDBEXq0p4n7dHx1p8SRUR7lGbwX0eazRj7/v/tkjItSH8vrrz2eG3S4/P/buXvWqIIojOP/U4j4AQJibRnUb2Aq00RTBFJYWNlYBixU8jnEWuwCIhGEdEkhsTZFmmBrIImC2iQQH4szy97c3SSX63I3uzw/uLB7X4ZZ2N3DmTkzJ8DDcq1evbZaaa/ttc/l3NMhfXzR4LuwxWCV3Z1a/2fG/XvwMZmHCwkmiKR94NOQS+dueyNpjyxZXiH/fI/IuYEDstz6LbBMfxjusj6cSnpMZibrwHcyEP4BvpJ7k81K2q4995Nc2/GErKbaL8/9IrOxN5ydKL+oD+tkYcAOcHzJ7ZDZTlXbXaur64EWI+J6k4fKZ38AfCunrgFrETHfsh/nioibZFZ1Sq4PqvtADs1ukCX2v8u9P8hg9RyY12CV3XK1DUn1OS6zRkJqUlBjNrlK+XCvcOEYuCWpXtzQtK0v9BfqLkl6P4IujkxEPANeA5uShm1n06bNIEu0e4uM53TFd9e2q8uZjk2liLgREXMRscjZzOZd24BTvKy8bpQddqw3nzPKYPiIfsDZcMCx/+FMx6ZSRNxmcAHlIXBXUqNhvAva/ggslLf3JdWr86ZKRGyTQ3Z/gXuSdsbcJZtgDjo2lWpB55AsYnglaXd8vTIzBx0zM+uM53TMzKwzDjpmZtYZBx0zM+uMg46ZmXXGQcfMzDrzD8qvTuKGS6I9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate the errors\n",
    "vel_error = [(TestSetLabels[i,0]-test_predictions[i,0]) for i in range(len(test_predictions))]\n",
    "# Calculate the point density\n",
    "x = TestSetLabels[:,0]\n",
    "y = vel_error\n",
    "xy = np.vstack([x,y])\n",
    "z = gaussian_kde(xy)(xy)\n",
    "# Plots\n",
    "plt.scatter(x, y, c=z, s=5, edgecolor='')\n",
    "plt.xlabel('Velocity (km/s)', fontsize = 20, fontweight='bold')\n",
    "plt.ylabel('Residual (km/s)', fontsize = 15, fontweight='bold')\n",
    "cbar = plt.colorbar(extend='both')\n",
    "cbar.set_label('Density', rotation=270, labelpad=10)\n",
    "plt.ylim(-50,50)\n",
    "plt.savefig('Vel_Residual.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEPCAYAAACOU4kjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeVhc1dnAf+8MDBAgQIBA9n2PW4Ixi82mcde41rhVrftS7eJaW6vWtNq61GqtX2rVuDWm1iVqXKMmLonZjNmJ2UNWAkmAELaZ9/vjDgQIM3OBgZlhzu957sOce885950B5r3nvJuoKgaDwWAwBBNHqAUwGAwGQ9vDKBeDwWAwBB2jXAwGg8EQdIxyMRgMBkPQMcrFYDAYDEEnJtQCtBQZGRnas2fPUIthMBgigCVLluxV1cxgzScicvapiZ47bkpj7ORtEqx5Iwlpq67IOTk5unjx4lCLYTAYIgARWaKqOcGab9YrnfXvz+9HED6dezAqlYvZFjMYDIYgIiLyyN/38a/HsuiQ5mDuO93a5hN8AIxyMRgMhiDy7sudPEMHuujRLZbf/aoDDz9ZGGqRQoJRLgaDwRAkqlct997WAYAhA+KidvVilIvBYDAEidqrlmqidfVilIvBYDAEgfqrlmqGDIgjLTX6Vi9GuRgMBkMQaGjVUk00rl6McmkCr/75HS7qfgtXH3MXq+avC7U4BoMhxPhatVQzdGD0rV6Mcmkky+au4ZWpb1NUWMKODbt56NKnQy2SwWAIPWf6WrVUE22rF6NcGsneHXX/OA7sLaaivDJE0hgMhjDh/J9fmuK3w9CBcZQe8iAi/ju2EYxyaSTDJg6hQ9bhv42x54/AFef7acVgMEQFEh8XOBDf2ycqIvbbbG6xlqJDVipPzf0D895aSHJaIidfMibUIhkMhjCgSt1UapXfPh6ixuRilEtT6Ng1nQtvOz3UYhgMhjDCg0aV8giEUS4Gg8EQBBQPHjwB+kSP8jHKxWAwGIJApXqo1ADKJXp0i1EuBoPBEAw8gDvAyiSKdItRLgaDwRAM7Nlcoke9GOViMBgMQaBKlcoA+17+N83aFka5GAwGQxDwoGZbrBZGuRiimuJ9B3nhdzPYs20vEy4ezcmX/STUIhkiFDfgjibtEQCjXAxRzV+ufpaFHy0DYPEny0nrmMLwSUeHWCpDJKIE3vaKJt1j0r8YoprcxRv8tg0Gu1SpUBngsOOKLCKniUiuiKwXkXsauB4nIm94r38nIj1rXbvXez5XRE4NNKdYTBWRdSKyRkRua+7nUI1RLoaoZvCo/jWvRYTBI/v76W0w+MaNBDwC6RYRcQL/AE4HBgOXiMjget2uAfapal/gSeBR79jBwBRgCHAa8KyIOAPMeRXQDRioqoOAGc37FA5jtsUMUc1dL9zEyw+9yZ5tBYy/aBTHThgSapEMEYrHq0D8oYFzVo4A1qvqRgARmQFMBlbX6jMZeMD7+k3gGRER7/kZqloObBKR9d758DPnTcClqlb0p6ruCfxO7WGUiyGqaZecwI1/vSLUYhjaANa2mP/NIO/K5UoRqf1HN01Vp3lfdwG21bqWB5xQb5qaPqpaJSIHgHTv+QX1xnbxvvY1Zx/gYhE5D8gHblPVH/2+CZsY5WIwGAKyNXcH389dQ9e+2QyfaFZ3DeHGgTuApcG7cpmuqk/56NLQ0qb+bpqvPr7ONyRU9ZxxQJmq5ojI+cALQFBcJo1yMRgMfvlx2WbuOOMvlB+qAOD6qT/l/JtPCbFU4YdHBY82u1RLHpYNpJquwA4fffJEJAZIAQoDjPV1Pg/4n/f128CLzZS/BmPQNxgMfvnize9qFAvAx698HUJpfLP8q7U8ddtLvPbIu1SUVQQeEGQ8QTDoA4uAfiLSS0RcWAb6WfX6zAKu9L6+EPhcVdV7forXm6wX0A9YGGDOd4CJ3tfjgHVNfPtHYFYuBoPBL2kdU+q124dIEt+sW7qJe8/5C1WVbgA2rdzG7169tVVlqFInler/KzWQQd9rQ7kV+BhwAi+o6ioReQhYrKqzgH8Dr3gN9oVYygJvv5lYhvoq4BZVdQM0NKf3lo8Ar4nIr4AS4NrGv/OGMcrFYDD4ZfINJ5G7ZBPfffwDXftm84snws8B4od5a2oUC8CSz1e2ugxuW95igVHV2cDseufur/W6DLjIx9ipwFQ7c3rP7wfOtCFWozHKxWAw+MUVF8vvpt8UajH80nNwF7/t1sCjgjuQt1jzbTIRQ1gqF2/Qz2Jgu6qe5d0/nAF0AJYCV6hqq2yq7t93kPfeWoKqcvb5OaR1SGyN2xoMhkZw/CnHcPNjl/PZf74ho1MaNz/e+qurKpxU4vTbx0acS5shLJULcDuwBqje3H0UeFJVZ4jIc1gRqv9saSEqyqv49c0vs21LAQBffLaKf750HfHxsS19a4PB0Egm3ziJyTdOCtn93dhYubSSLOFA2HmLiUhXrD3A571twfJmeNPbZTpwbmvIsnXz3hrFApC3tZAtm/Jb49YGgyHCUBx4Ahxm5RJa/gbcBSR72+nAflWt8rZrR53WQUSuB64H6N69e7MFSc9MJi4uhvJy69YuVwyZYegpYzAYQk+lOqhQsy1WTVitXETkLGCPqi6pfbqBrg2uLlV1mqrmqGpOZmZms+VJ65DI/X+6kJ69M+nZK5PfT72ADulJzZ7XYDC0PTzqCHhE07ZYuK1cxgDniMgZQDyWzeVvQKqIxHhXLw1FrLYYI0b1ZcSovq11O4PBEKF4bKR/afhZuW0SVisXVb1XVbuqak+swKDPVfUy4AusSFSwIlPfbU25KiqqeP75udz3uzd5550lgQcYDIaow60S8LBTz6WtEG4rF1/cDcwQkYeB77EiVFuN5/7vc955ZykA8+evJz4+ltNOM9UKDQbDYexE6HvMyiX0qOqXqnqW9/VGVR2hqn1V9SJvvYJWY9Wq7XXbq7f76GkwGKIVRfAEOKKJsFUu4cTgQXWd0waHIPrXYDCEN9VxLv6OaPIWi5RtsVZj5co8Nm/Zy9FHdaN793QAbrppIgntXGzelE/O8b043WyJGQyGeljbYoFckaMHo1xq8eFHP/DYEx+iCnFxMTzx10sZNKgzLlcM1183PtTiGQyGMMZKuW+rWFhUYLbFavHeB8tqvDnKy6v46JMVoRXIYDBEDHbiXDCJK6OTlPbt6rRTUxIa7FdysJyFizeSnJzA8cN6toJkBoMh3LGzLeZpJVnCgUYpF2/eryxvc7eq5gVfpNBx6y0ns3vPAbZuLWDYcT2ZcvHII/oUl5Rx069eYdv2fQBccM5wbrvxpNYW1WAwhBlusFHPxaxcahCRScBVwElAZr1r+cBnwEuq+llLCNiadOmcxgv/uhaPR3E4Gv4jWLBoQ41iAXj7/aXcdM14YmP9P7EYDIa2jVZvfRkAPzYXETlbRFYAH2FFy3fEyl1Q++gIXAJ8LCI/eHODRTy+FAtAUmJ8nXZCfCwxMeYPymCIdgK5IRtX5MNUp1hxA18BC4CVQAGWYkkHhgIjgROBo4B3AswZ8Yw8vjdnnXY073+0nISEWH77mzOxqgJEFnPnr2NrXiEjjuvJgL7ZoRbHYIh4qnAEdkWOIl9kf4pgG/Ak8Kqq7vU3iYhkAlcAtwVRtrBERLjzttO47YaTiIlx4nRG3qrl1TcXMO2VrwB4aca3/O3hizlqUHgFhla5PXy2MJdD5ZVMzOlHSlLDzhUGQ7igKngCeIOZlYtFb1V125lEVfOBJ0TkqeCIFf7ExUVuNcpPvlxd87qyys2X3+SGnXK579n3+WLJegBe/XAx0/9wGUnt4kIslcHgG7etrMjRg89Pwq5iae4YQ+uTlVm34FnHzGQfPUPD/pJDNYoFYNvu/SzNbVOOiYY2iFsdVKnT7xFNKxdbalZEBonI+SIy3NseJiIfi8hSEZkqIkZdRxB33HwKRw3qQkpyAqdNHMIFZw4LtUh1aBcXS0K9lWGHejFIBkO4YSvlfqiFbEXsGt8fwKqncpuILANmAZ2wDPvHAAeBP7WEgKGmvKKKvz/7KatWb2fQwM7cfvMk4uMjd0sMrJXLPx65NNRi+MQVG8PDN57Bwy98Qml5JVedOYKhfTqFWiyDwS+2bC42tIuInAY8BTiB51X1kXrX44CXgeFYDlYXq+pm77V7gWuwHLFuU9WP/c0pIi8B44AD3umvUtVlgaUMjF3lMtz781MgB+gM7MSqZ388ljtym1QuL736NR98tByAzVsLSEqM45YbTNBkS7H/4CHeXrSKGKeDtx+7hnZxrlCLZDDYwoPYiHPxr3xExAn8A5iE9f26SERmqerqWt2uAfapal8RmQI8ClwsIoOxwkaGYH1HfyYi/b1j/M15p6q+afuN2sSucqn2Vd0CjPe+/hMwE9gNdA+uWOHD1m2Fddt5hT56GppLaUUlV/zjDTblW0Gqs7/P5ZVbLiYmRB55s75dxfrtexk5uAejh/QMiQyGyKFKHVQGUC42arqMANar6kYAEZkBTAZqK5fJWLtJAG8Cz4gVDzEZmOGtd7VJRNZ758PGnEHH7n9ttaE+BWsbTIE1QHEj54k4Ro3o47dtCB6r83bXKBaAFdt2kVd4wM+IluP5D77jwZc+4bVPl3Lb399m3vKNIZHDEDnYS1wJwJUisrjWcX2tabpghYFUk+c9R0N9VLUKa0sr3c/YQHNOFZHlIvKkd8stKNhduWwEjga+8QqlwDKgq/f67mAJFG6cdfoxJCa6WLlqO4MHduakCYNDLVKbpWP7JJwOwe2x/gPjY2NISwxNfMvcZRtqXqvCvB82MPbo3iGRxRAZKIFXJl6Ty3RV9RW20dAE9S01vvr4Ot/Qw3/1nPcCuwAXMA2rpPxDPmRrFHaVy7+AZ4Dqx/ZZqlooIhd624uCIUy4MmHsICaMHRRqMdo83TNSeeiiU3j642+JcTi4Z/I4UtrFBx7YAnTLSmX1lsPPTN2z0kIiR7jx2dw1PPfilyDCLddMYMKJA0ItUthQ5XFS6QlULCzgtlge0K1Wuyuww0efPBGJwdpRKgwwtsHzqrrTe65cRF4E7ggkoF1sKRdVfVZECoAxWHaXZ72XCoAHsZJXGgzNZnLOYCbnhH51ePclE6msdLNhRwEjB/fg0pPCy107FOzJL+LPf5tNVZWVOH7q4x9w9JAupKclhViy8MBjw1vMBouAfiLSC9iOZaCv79o5C7gSmI/lxfu5qqqIzAJeF5EnsAz6/YCFWCuaBucUkU6qutNrszkXK8VXUPCpXETkN8B7qroOQFXfAN6o3UdV/wf8L1jCGAzhQkpiPH+96exQixFWFOw7WKNYwMrusG9/qVEuXjzY2RYL5KqsVSJyK/AxltvwC6q6SkQeAhar6izg38ArXoN9IZaywNtvJpahvgq4pTqwvaE5vbd8zZu+S7BMHTc2+o37wN/K5a/AX0RkI/Ae8D4wz2tAMhgMUUbvnpn06ZnJhs35APTr3ZEeXdNDLFX44FYnVYG2xWzEuajqbGB2vXP313pdBlzkY+xUYKqdOb3nJwaWqGn4Uy5DgbOBM4FfALcDxSLyKZay+dCbU8wQpewsKiYhNpbUhNDYRQyNY09BMa+9v4gqt4cpZwynW3bj7Ehxrhj+/sglzP50BQ6HcPrJQ4+oY1RSWs4Pa7eTkZrIgN5ZPmZqm9jZFoum9C8+lYs3wGY18KiIpAFnAGcBpwAXAB4RWYSlaD5Q1R9aQV5DGKCq/HrWh7y/OhenCA+cOpEpxx0darEMfiirqOTmh95g+x7LtXveovW8/thVJCc27sEgKTGOn56b0+C1/UWlXHfff9i+ez8At/1sPFPOGt5g37aIInbiWKIGW/EpqrpPVV9T1UuwqlFOwErHnwI8DCwVka0tJ6YhnPhq4xbeX50LgFuVhz75gvIqs1sazuzYfaBGsQDs3X+QjXkFQb3HJ1+vrVEsAC+9tSCo84c7Hj28evF1mNxiflBVDzDXe9wpIr2xts/OCLJshjClwl03+XWlx4MnBFWQdhwsQgQ6tWsfuHOU0zE9meTEeIoPlgGQEBdLl44pQb1HfFyM33Zbx60Oqjz+n9e1+d5kEUOzI+tVdaOqPqWqpwZDIEP485PePRjW5XAiyRtHHU9CbOsm83xoyaec+O4zjHnnGR774csjrhdXlh+hBKOZpHZxPHnP+Qwf0p1jB3blr3eeS4YPL6+S0nIen/Elv/2/D5i/crPte5w+bggjj+0JWMrr7utPCYLkkUOgVUsQ3JQjCluPFt6U+tdi+VR3A+qnCFBVNXlRooS4mBheuewilubtIMnlYmin1jXc5u7fw0u5h+N2n131LRf3OZZuSamoKnd99wH/27SceGcMj486h9O7DWxV+cKVIX078czvGnQyqsPd/3yf71ZvAWDOkh958bdTGNwzcCns2BgnT/z2Agr2HySpnYs4V2RnD28sdmwuZlvsSB4BfuN9bSc9gSHIzFm8ju17DzDmqF706ZLRavf9ctMmthUdYGyPnvRITa0573I6Gdmjm5+RLUelx9PAOWuV8vmO9fxvk5XFusxdxV0L3ufUrgNwSMP/9GVVlfxjxXy2Hyzm7J4DmdDVPCMtXXe4MJvb7WH5+h22lEs16amJLSFWDeUVVbz63wVs37WfsSP7MX5MeGQJqPKI2RarhV3lcjmWUikBVgAVLSaR4Qiee+cbnn/vOwD+7935vPjbS+jfLbPF7/v0ggU8Of9bAJJcLt6acgl900Mf1zAkLYszuw/ig61rALig99H0bm/JVVJZ90/zkLuSKo8Hl7Ph+IPffDObDzavBeCdjat449RLOT6ra4N9o4VBPbNYvt7KGiICA3uEl0vxE//8lA/nWIHkn81dw1//cCEnDO8VYqns1XOJJuwqlxis1clQVTVeYa3MB9+uqXldXlHFnMXrWkW5zFixvOZ1SUUFH6xbx+2jRrX4fQMhIvx9zLn8fOAIHAjHZHSuuTaxS1/6tc/gx6K9AFzd/3ifigVg/s4tNa89qizYtTXqlctfbj6bp2bOpaColHPGDOHYfvWT8oaW71dsq9teuTUslIu9OJfowa5ymQb8FugPtJhyEZFuWBXWsrGyKUxT1adEpANW6pmewGbgp6q6z9c8bY3MtCR2FhTVtDumtU7N+/R2iewsKalpZ7QLn1LDIsJxGUd+6SXHxvG/U67km12bSY2L54SOPfzOM7hDR76upWAGd+gYdFkjjYyURP54Xfg6f/bv05Fdtdyq+/cJj5WVIja2vaJnZWNXudwPHAV8KCLfYCmY2oENqqrXBEGeKuA3qrpURJKBJd6MAFcBc1T1ERG5B7gHKzV0syncW8Kjv3+LTRt2M/yEPvzq9+fgcoWXC+Ufrj6V+6Z9wPb8A5w0vB/njh3aKvf9y6mncMv777O9qIjT+/Xj4qOOapX7Npek2DhO7WZvH/6psefwx0Vz2FFSxNm9BnNSt74AlFSWU1xZRnZCe8SHvcYQGu7+xWkkJ8WzY9d+xo7qz8QTw8Nhw61CVYBiYSHw2A8ZojberYhcCbzgr4+q+k+q0wRE5F2sVP/PAOO92Ts7AV+qqt9vj5ycHF28eHHAe/zx7pl8/fnhbacrb5zApdeMbZ7ghohmzo51/HLBW5S5qxjdsRfTTryYOGd4PXAYgouILFHVhlMP2Bv/4jHPXXVVYm//K99Vd83gwLKtaaq632/HNoDdOJeHsNZz/o6gIiI9geOA74Cs6roD3p8N/gZF5Prq6m75+fbSnu3ZdcBv2xB93L9kNmVua2H+7Z5NzNq6IsQSGSIBO3Eu0ZRbzK5yScOyRZ0NJKiqo/4RTKFEJAkrlf8vVbUoUP9qVHWaquaoak5mpj2D94RTD2/1OJzC2JOHNFZcQ5hxsKqUr/K/Y+m+pimFcnfdVDZlbpPaxhAYjzpwe/wf0bQtZnetPxO4GlitquUtKA8iEoulWF5T1be8p3fXKmrTCdgTrPudf+lIsruksnnDHo7N6cXgo0MTu2EIDgerSrlvxaPsLLOqSJ6c9ROu631Zo+b4xZCxTF32CQr0Tk7n7O6tY+MyRDaq0WVTCYRd5bIEq0rZl95SmPUN+qjqy80VxlsN7d/AGlV9otal6sprj3h/vtvce9Vm9LiBjB4XHkbBSGZRwQZ2HNrHCel9yU5IDTygBVi2f1WNYgGYs/trruz5U1wO+9HiV/YbwciOPdlzqJjj0ruSFFs/IYXBcCQekxW5DnaVyz+wtsXSgN83cF2xXIibyxjgCmCFiCzznvstllKZKSLXYCm2wDksDK3K9I1z+ce6jwFIjInjjkGnM6nTMFyO1jWEJ8bUdZeOc7iIkcb7mgxI6ciAFOOWbLCPxyO4TYR+DY35z2/xT0VVv/Zzn5Na+v6GpvPfrYfTqx+sKudPq9/kw50LeGLYzc1WMKVVRWwrXUOqqyNZ8f6D5Y5NHcKp2eP5ZNdc4p1x3Nr3ahwSVJOgwdAgZlusLnb/6ye0qBSGiCc1th17yg572gnKqqItLClcx6iMwU2e90BFPv/eeCclVYUIDs7qcgvHpU3yO+bnvabwsx4X4hSniVExtBqqgYMoo0n32FUua1V1t6+LIjI5SPIYIpTfDT2fO79/ld1l+3GKhxiHlVwy3ulq1rzf7/uUkqpCABQPX+f/t0a5FFXsYvPBBSTFZNA7+cQ642JaeTvO0DIUlpTy0Q/rSIyL5czjBhHjDOdVqJ0I/ejB7n/gXBGZUB1rUhsRuRl4Coiu/NoNUFZWyd6CYjpmtg+7KP+WZmBKF94bfzfvbf+Wp9a9jVvhzM4jOS6tb7PmdTni67UTADhQsYOZm2+m3GN5qg/rMIXRHa9v1r0M4cWB0jIufeY/5BVav+M5qzbw9yvPCbFUvnF7BLcnwMolipSP3W/A/sA8EZmoqjVZ40Tkz1hpWKJptdcgmzbnc8e9b1BYeJBO2Sk88eglZGcHt9JfJHB2l9GcnDWcCk8lKa6Gi1E1hpz0M1hXvJCtpatJcCZzeqcbANhQ/FWNYgFYtX922CkXt3pwGntPk1m0Ma9GsQB8vmoDB0rLSGkX72dU6LBsLtGjPAJhV7nsBfrgVTDANuBF4FLv9aUtIFtE8fyL8ygsPAjAzl0HeOU/33Lnr04PsVShISEmjoQj6sk1TGnVPkqrCkh1dSfGceQWmssRz5W9/kyp+wDxjkScXpfixJgOdfq1i0lrvuBBQlWZuuID3tyyhPax8Tw6/EJGZZo6MY0lM7luXZjEOBftwrgAmR2bSzRh97HqRCwX4B7AV8BnwGVYnl2zgHEtIl0EUVVVr658Zd325m0FfLXgR/YWlBBJHKwqZ9a2H/h4+yqqPMEtG7ylZAGvbryEmVuu479brqfM3XDqHREhMSa1RrEA9G9/EoNTzsApsaTEdmZS53uDKltz+HJ3Lm9sXoRbPeyrKOWeJf8LtUgRyTE9OnH7aWNIineRlZLE45efSWxM0FMYBg21k/7Fxh6PiJwmIrkist6bqLf+9TgRecN7/Ttvqqzqa/d6z+eKyKmNmPNpEQnql5OtlYuqrhORUcDHwFCguoDGk6r6G98jo4fLLxnNilXbOXSogvbJ8Uy56ISaa19+k8uDj7+P2+2hfXI8z/z5Unp2C33RrUCUuSv52dcvsvbALgAmZA/g6RFTguaB9W3+c1R5Ez7sq9jCyv2zyEm/wtZYEQcTO93BxE53BEWWYLK/orROu6jyEB71tJpL9Mw1K1i4M4+jM7O5YuixEe0xd93EEVw3cUSoxbBPMw0EIuLEiiucBOQBi0RklqqurtXtGmCfqvYVkSnAo8DFIjIYmAIMwfqO/kxE+nvH+JxTRHKAoEc9+1QuIvKzBk6/ADwMJGCtYH6o7heMCP1I5qihXXn539eydVshvXpmkFar1Ovrby/E7ba8p4qKy3j3w2Xcfn34h+18X7C1RrEAfLErl91lRWQntIwtSfXI8sWRyLisAXRKSGHnIWsldlHPnFZTLK+v+oHfzv0UgDfXrqK4ooJbhp8QYJQhGNjbFguo6EcA61V1I4CIzAAmA7WVy2TgAe/rN4FnvNlNJgMzvCm6NonIeu98+JrTq8z+imXiOC/gm2wE/lYuL+FfD//Ee0DwIvQjmoz0ZDLSjyzklRBXd584ISF8941rkxpXN9o91uEkMSZ4qVBGZlzLJzsfwq2VpMZ2Y2hq2/Bo7xCXyIyxNzB39zo6uNoxLrv1arx/k7fliLZRLq1DI4IorxSR2kv0aao6zfu6C5ZNu5o8oP4vsKaPqlaJyAEg3Xt+Qb2x1RX1fM15KzDLm7fRlvB2CbQtFrnr6TDi1msmcOdD/6Nw30EG9MliyrnHh1okWwxK6cRtAyfybO6XuBwx/OGYs0iODZ6nTq/kE7k84XVKKvPpENeLWEd4egE1hQ5xiZzX/bhWv++g9I58sGEd1c+FgzJavhx2qKhwu/2WsG5t1ONAA6Z/AWC6qj7lo0tD37n1VZavPr7ONySUikhnrFRa433I0iz8KRcTlR8k+vXO4s1/30hx8SFSU9pF1B74DQPGcl3/ExGkReROjMkgMSajSWMrKn9EcBAbazyxqrn+uON5a/NKNhYX4hBhcMe2p1x2FRVzw+vvkLt7L0M7Z/F/l55LemJ4lOAOQvqXPKB2avauwA4fffJEJAZIAQoDjG3o/HFAX2C993+7nYisV9XmBad58alcVHVuMG5gsIhxOurYYSKJcMzNtXffHZQcfA2A9knX0SH1oRBLFB58nreBjcVWRgOPKvcv+JTz+w6JiAeagtJSpn41lx3FRZw7cBBThh7dYL8n5nxN7u69AKzcsZunv5jPA2eFgQ1TCUbE3yKgn4j0ArZjGegvrdenOkv8fOBC4HNVVRGZBbwuIk9gGfT7AQuxVjRHzKmqq4Ds6klFpMSXYhGRSqDMxztMbKgSsT+DfoqqNqosY1PGGAyNpaJiVY1iASgq+RfJST8nNqan33GVVTsAJTami99+kUxZVWWddrnbjUcVZwQol19+NJtvtm0FYOH27XRKSmZczyMTlR44VLek1IGyMgAe+3Ae736/ho7tk3jkolPpl9W0FXFTUQ9ogAj9QMrHa0O5Fcsz1wm8oKqrROQhYLGqzsIqS/KK12BfiNCqev8AACAASURBVKUs8PabiWX8rwJuUVU3QENzNvLtrVDVYQ1dEJEG4xz9PZJuFZHHReSYQHcVkWNF5ElgS6C+BkNz0Yb+QwN4muXv/xMbd+awcefx7Nl3fwtJFnomde/HwLTDW2E3HzMSpyP8Vp4NsSp/T712w6XKLx5+FDHe9xTrdPLTYUfx4fJcXvhqCQUlpazZsYc73/iwxeWtT7W3mL/DjhlbVWeran9V7aOqU73n7vcqFlS1TFUvUtW+qjqi2gvMe22qd9wAVf3Q35wN3NdfSg1/SQIbvObP5pIM/BL4pYhsxlqCrQIKvNfTsWJeRmEFVwrQNnxJDSHH7d5J2aEPcDg6EJ9wLlJray7ONZTEdj/lYOlMAJITryI2trfPuSqqtlBY/ExNe1/J86QkXUZcrH8vLo9WUeEuIs6ZFhHbSgCJsS7ePutyvtu1jbS4BI7J7BRqkWwzsms3Plr/IwAOEUZ06cKbq1YxbckiEmNdPDBhIsdkZzNxQB/evO5SVu/aw1Gds+jXMYMXvlpcZ66d+4tb/w0EZ1ssXKkUkUxVraPxRSQDqGxogD/lMhz4M3AK0Avo6aNf9X/dx1iFvQyGZuF276Zgz+l4PNaTbHn5PFLT/lanT2aHp0hJvhFw4AqgJNAj//ZVK/wO2V++nnk7buOQO5/UuAGM7/wMcc7QVNdsLAkxsYzv6lvZhiuPn3IavVLT2FlSzFn9B5AY6+Kezz7B47WSXzvrHeZfez0xDgcDszMZmH14hTZhYG/++fl3HCy3fq9nHN167t+HEQgU5xK5yudl4FURualWvEwvYBo+wlD8GfS/B04TkUHA1VjFuo7G2rMDcAPLgTlYe3hrg/UuDNFNedkXNYoFoKz0TTT1caReRUlX7CBb87li+9K+3U8pql7pJJxNXOxQv2OW7X2SQ27rIW1/eS5r973MMRm3saH4G7YdXEpGfF+GpkZn7jiAOdt/ZPbWNXRNTOGmwaOJj2l+7FZCbCx3jjlcOuGj9T/WKBawDP7F5eWkJSQcMbZXZgdm3DSFz1ZvIKt9Iucc2/QaQk2mDa9cVPVJEUkGvvMGXla/22eAvzU0JmD6F1VdA9wFINbeRHXekgJtKyHVhrDC6cyq03Y4Mo5QLI2lU/rfSE26EvAQ7xoWcJurynOoTrvSU8q6oi/5cMfDNedKq/YxIqO+I0/bZ8HuLVw/778136NbSvbzt9HBD4Ad3qkz6QkJFByyfhc5nbs0qFiq6dMxnT4dQ5hWSW2sXCIYVX0IeEhEMr3tho1iXhpVdMSrTPxOaDA0l7j4CbRLuJCqsveAeBJTnwzKvAlx9oMaB6Zdwbe77kVxE+tIpm/KhSwsqJuAcvPB7yJKuczZkcv721bROaE9Nw8eS2JM0wq5Ldizpc4D+vzdm4MiX30yExP578VTmLlyJYkuF1cd2/pBqY3BToR+pJZB9qaXuRY4FWvF8omI/NvfAiO6KloZIgJ35Xoo/4gY3MBBKg/cBQc8SEw3YlP/hsT0RCuWoCXPAA4k+XYktuGYiKbSNWkCp3WfQXHFFjrEDyYhJpO0uO51+nRwdfcxOvxYlL+Fm7+deXi1cXAfz4y6qElzDU7LqtfO9tGzYT7cksvc7Rvpl5rB1YNycPhZRfZMTeOuE3/i83pYoQKBXJEjN+nJw8Bg4FngaeB74BG8u1oNYZSLIexwVywGLatpezy7UZxQuY/KA3cTm/osuu9aUKt+jhaugMw5iOPIvG7Nob2rJ+1dPWvawzv8lNKqfV6bSx9+0vHGoN6vJVm8d2ud1cbC/KZHDZzSdQAPDD+FD7w2l/uGnWx77Edbcrlp7js17b1lpdw9rG1U7BC1Dr99WkeUluBsYJg3DueQqv5JRL7zN8AoF0PY4YgdgPVvaP2nCtTYSNS9E9xbaxSLdXI/uLeDY2DLyiVOxmXd3KL3aCmGpNV1SR6a1jwX5Z/1z+Fn/XMaPW7ezs112l/t2NRmlIstg36EbosBoqpVNQ0RF/ivCGiUiyHsiHEdR3zqE1QefBUkBmfFcqpd6Z0JkyGmDzgyweM1/zk6g7Nlt6jW7J/J1pJ5tI/tzvCMW3A5IyuVz9jsvkwdfhbvbV1J53Yp3HW0/dVGMBmQWjdqvn9qG8p95sHGtljEskdE+qnqj0B74BusGjE+McrFEJa42l2Aq90FAHgq1+Ap/wJx9kA8+bjzJwAxOFxjwNkVSbwecbRc4sKNxZ+wMN/yttx1aCkVnmLGdfpji92vpbio13Fc1Cu0RvErBgwj/9BB5u3YRL/UDP5wfGiUXIvQhl2RgXOxwk8AbgB+VFW/e6v+cos16lFQVbc2pr/BYBdH7CAcsYPQyrW4C35J9X+wp3IxztR/tKhiASgoy63bLs/10dMQCIcIdxw3ljuOGxtqUVqGtrstNgxqtqcrgZ4i0tNfgmN/K5fN2P8oNMBchiDi9nh48K05fLRiHV3S2vPYpWeE1r+/lVDPHur8Seoh0CLAv3IpKZlGSfFzOBzJpKY+hiuucfV0shOOY/X+/9RpGwz1EY8gbXdbrHY5+zisCpfL8FOaxRQLi0DeWryK/y1eCcC6XXu577+fMOOWS0IsVcsjscPA2Qvcm6y2aww4svyOqShfTNGBBwDweHZRWHg1WdnL6+QqC0S3pBMZm/0gW0u+or2rG0enXdnk9xBJ7Cw9wB+//4Q9ZSVc2OsYpvRuMCmuoZo2vC2mqufUbotIJ+Dv/sb4Uy4PBkMoQ/DJLyrx224Kb27+nllbl5Od0J67jzqFb/eu4du9ufRMzOTqPhNxOUK/MBVHEs70meih90DikIRzA0bau915ddoeTyGqhxBpnEG+V/IkeiVParTMrcWhqkqcIricwfs93frt/1i+z6o19UPhdronpjE668gU+IaoZBdW4mKf+MstZpRLmHLKUf14Yd4SDlVYHlTnDGteHqWvdq/nd0vfq2mvL8ljR8Xh4ncHKku5a/C5zbpHsBBHGpL4szrnyiu+p6RkGiLxtE/+DTExXWuuueLG4HB0rMlVFhd3Mg5HZHl6BeLxFZ8zLfcbYsTB7487PWgrjNwDdVPgryvKN8rFD3biXCJ1ZSMif+fwTpYTOBZY7HtEE+wkItIROKLYuTHotx59szJ449ZL+Cp3M13S2jNpaL9mzbdm/6467e2HCqidyuuHfZubNb8/3J4yqjz7cTk7Nmqrqpoq9w7y9/4UVWv1Vl6+kOysuVjVX8HpzCQj830OHXobhyTTLjFy0rXYYeW+nfxf7jcAVKqHB7+fzWldB5Hq8p2Dyy4nZvdmzo51AMQ6nIzIjJyMBCHBVoR+xFJbkVQBr6rqt/4G2FIu3iyYD2O5oKU00MUY9FuZYCbpy8nogQPB432s6p2UxaZDh58VBqd08zW0WewvW8SqPTdR5SkiOe5Yjs56gRiHv3pFUFaZy77iFxCJJ6P9rVRWrKpRLABV7o24PfnEOA8HCcbEdCUp6WYOVm6iwn2AuJi2E1tRUlm3KqNblbX7drP9YBGD07IYlObfJuWPJ084j2m535JfVsI53YcyOLVxaV7q41HlL4vmMTdvE/3SMvjj6JNJiTviOTVyacM2F+AToEhVS736oH2gAXYVwp3A3c2RLBiIyGnAU1jLsudV9ZEQi9QmGJbejX+OnsL721aSnZDMTQPH8t72RXybb9lcbup/aovcd33BH6nyFAFQXL6MHcWv0z3lep/9q9z5bNlzPm7PfgAOln1F9wxL0ag3XYzT2RWno26gnttTzpLd17OvbBFCDEMyHqJLcnhs8zWXYRndGJbelaUFlm1pZGZPrvriDSo8bmLEwbM/OZ+Tu/Zv0twJMbHcPiR40fPTV3/PP5cvBGB1YT4CPDXhrKDNH3La8LYY8C5wlohUAguBBBF5R1Xv8TXArnK5HOtjmQ2c6X39BHAFUAS82hyp7eDVlv8AJgF5wCIRmaWqq1v63tHAuOx+jMs+vL12cY8xXNxjTLPnLSjfRqm7iE7x/Yhx1M3C66mVPwzAUy/NfX0OVayoUSwA5VW5ILFkpL9GcfE/EUcCKe3vQaRubZFdBz9kX9kiAJQq1hb86QjloqoRU22yNi6Hk+ljr2DervW4HDG8uX45FR4r1q1KPbyUu7jJyiXY/Lhvb532/B3bcHs8QSvDPP2rpcxaupqO7ZP4/XkT6Zwa8OE6uHhoy7V441U13/uAv1hVrxORVYBP5WL3t1pd1u6a6hOqeidwFtAXaI2aoiOA9aq6Ua0ygjOA4BeRMASNRQXvMG3Djby6+U5e2XwHFZ66yqRH6s1U155zObPITv6p3/niYvsgtdIZxTg64nRmEB83msyMV8joMI3YmCMrMGpNYHF1u+43wLw9r/HomvN5fO0UcosWNOYthgVxzhgmdRnIuE592VZyoM61HQeLQiTVkYzvVtcZYHdJCY8umBeUueeu3chfPpjL2p35zMvdxG9e+yAo8zYG4bBR3+fR6lIFDxHpAFyCtcgAH+WNq7GrXKoXc3urJxSRdKC6+uQvGydmk+gCbKvVzvOeq0FErheRxSKyOD/flJ0JJarK3D3Tqf7T2VW2nrVFdb9IspLOJafzewzt+C9yOr9HfIz/ZIqumB50zfgXCa4RJMaNpXvmqzjEb+48ALITT6e96yhvy0H/DofjwfJK1/BV/n9wayVl7hLezfsrVR7/JZDDmYy45JqaIarQITZ8PONO6dGPsZ16WUlE3IAHvtkeHD+gDbsL67Q37in00bMFqS4WFuiITP4KrAO6Au+LSHssO4xP7G6L7fVOmgpsB3oArwPV1sTWKC7e0G+lzg6mqk7DqulMTk5O5O5utgFEBIc46/yGHA38uSW6+pLo6mt73uSEk0hOOMlvn8V7X2b5/v8R50hmYvbddGp3FCd0foWi8jW4nGm0iz3s9VRaVffJvlLLqfCUHbGFFymc0q0fn+WtR1FAmNQtPLbEqjmlR1++2rq5pj0oPTjOFaP6difW6aTSba1SfzKgZ1DmbRRt2KCvqq9S1/xRiZ9aLmBfuazEUi6DsJZENwPVGecU8JlfJojkAbXdlroCO3z0NQSJOTtXM+3HL4lzxnDH4NM5Os2+59ik7JuYveMpPFTRM/FYBqW0fNGn7aXfs7DgRQDK3EV8tOMPXN33LRziIjX+mCP690g8inRXVwoqLIP4oPYn0i6mlffqg8hP+x1DnDOG7/Zs4+j0Tkzpd+R7DiWXDzmGwrJS5m3bTP+0DO4bPT4o8w7q0pEXrruA95auoaKiivNyhgRl3sYgHuvwS4QqHxH5Q0PnVfVBEblBVf+v/jW7yuVhYBawD7gP6A6cgbWt9hVwU5MkbhyLgH4i0gtr9TQFaFtBC2HGtoOF3LV0JlVqPQ3euvAVPj75DhKc9p7qj0o9iT5JOZS5S0hzdWowjuWQ+yBrDnxHnLMdg9qPwNGEWJfalFTW3Q495N6P21OB08dKJM7Zjqt6P8baom9xOeIZ2L75TgyhZnLvIUzu3fpfrnYQEW7PGc3tOaODPvfATpk88uMX5G7P5/1v13DNpOO57ewTg34fX1TbXJo9TwCvWBGJA14GhgMFwMWqutl77V4s27gbuE1VP/Y3p4j8G8jxir8OuEpr+/YfpiG7evVu0sEGrtlTLqo6H5hf69Q5IhIPxKpqaxjz8VZAuxX4GOsDekFVV7XGvaOVbaWFNYoF4EDlIfaVl5LQzv6WUbuYFNrFNBQaBWXugzy3/h72lm8H4NjUcVzU/fZmydwtMYd2zg6Uuq099z7J43wqlmrinUkcm3ZKs+5rCD1zV24kd/vhh4sXP1vMTWeMItbp9DMqiARhW8ymV+w1wD5V7SsiU4BHgYtFZDDWQ/cQoDPwmYhU74v6mvNXqlrkvfcTwK1Y5YvrvjXVJxqQ9QLvtQa9hZsc+KhWYEFZwI5BRFVnc9hTwdDCDE7pTEZcEnvLrQeZAe2z6RgfvFLC60uW1ygWgGX753J2l+uIdzY9hX67mA5c0OOfbCj+kjhnMgPaG6URLSS46rqgu2KdOJu5Em4MQdoWq/GKBRCRaq/Y2splMvCA9/WbwDNi+dFPBmaoajmwSUTWe+fD15y1FIsACb4kFJHzgauA2l8AOd4H/pdUdXr9MXYj9N0Buqiqmgj9Nkaqqx3TR1/HG1sWEueM4fJeo4lxBO8pMNFZ17bhcsQTGwRDenJsR47t4N+tOdxYk5/Prz6czc6SEs4fNIj7x0+IyLibUDJ2SG9OHz6AD5fk4opx8sAlk3A4Iu4zbMgr9gRffbw7OgeAdO/5BfXGVnvU+pxTRF7EMnOspm5q/dr8CbgRK64RLCX0OnAHlpniCOwqhIj7DbVlqtwefj/zEz5ebtVzeeJnZ9EvOyPwwCbQNbEDvxl8WovM3StpCGMzz+Pr/Fm4HHFc2O02nBKdzyi/+nA26woKAJi+bBnHdurE5IGDQixVZOFwCI9ceQZ3nT+eeFcs7eJiAw8KJva3xa4UkStqtad5PV3Bhlesnz6+zje0fKuZU1Wv9m7HPQ1cDLzYQP9SVf2yjhAih1R1SQN9AfvKpf6Sxwn0BMYAh4CZNucxBIG3Fq7kvaVrANiUv4/fvfEJb9wemb4Np3a6gknZl1puy1HM7pISv21/rNqxmxfnL8XldHLzuBPomtawjSta6JDcspVJfaKBt8W83/7TVfUpH13seMVW98kTK0NrClAYYKzfOVXVLSJvYKX6aki5jGrg3O98vAfAvkH/6obOi8hE4DOsimSGVmJv8UG/7Ugj2hULwHmDB/PS998DkORyMamPvdif/OKDXDX9fxSXWyFni7bk8cEtV+KKMZ9pq2Nn5RJ4ZWPHK3YWcCWWk9WFwOeqqiIyC3jda5jvDPTDygMmDc3ptbP0UdX13tdnczgwHgARWayqOV47TnWRsKuwUn+txY8NvFl7EKr6uYiUALcRoCqZIXicekx/ps9bysFyK5L8ghF+a/YYIoDfjxvPsE6d2VVSzEm9+9ArLc3WuNzde2sUC8C2fQfYU1wS8auXotIyvt+wnay0ZAZ27RhqcWxhq55LAHx5xYrIQ1g5vWYB/wZe8RrsC7GUBd5+M7FsJ1XALaqWu6ePOR3AdG+0vQA/cGRYSYmIPA18DVyGpbReAX6iqgX+3otdg/7PGjgdD5wC+M+Rbgg6fbLSmXn7pXyzbgud09ozfvCR+bQMkYWIcNaAAY0e1zezA+1iYymttNI8ZSUn0TE5fFK+NIW9RQe54rEZ7CwsQgTuvnACU8YdG2qxAhOkCP2GvGJV9f5ar8uAi3yMnQpMtTmnB8u04Y9JWImLb8faWpsK/EdVD/gdhf2Vy0v4/tiUuh4KhlagR2YaPTLtPd0a2i7ZKclMu/xc/vX1YlwxTn45cTSumMh1iig6WMarc5ays9BySlKF5z/+LiKUi9iwuURahL6qVmLZYF4UkYFYW2JLRWQB8KKqfuZrbGP+Cn15jC0EfBfhMBgMLUpOj67k9OgauGOYM3f5Bu7912wOVVVVJ8sGICEuQvK8teHcYgCquha4R0R+C5yOtYXWbOUyoYFz5UCequY1WkqDIQCbigp5Z9MqOsS145J+x+JqrShrQ8j4y4wvKausshoewAFJ8S7uu3hiSOWySzBsLpGAdzvtA+/hE7veYq2RmNJgAGD7wQOc++HLHKiwEkDM372F58adH2KpDC1NdUZjAfDADaeP5OrTjicuNkK2+ewUC4sC5VONz9+aDyO+T1T15eaLY4g2dh8s4Z0fV5PkiuOiAUNxOZ18s3NzjWIB+GTbj1R5PMQEqWKhITy5+ZzRPPzqZ3hU6ZWVxpQJx0aOYsHeyiWaotH9/eZewr6eVawsnYYoxe3xsL6ggPbx8XRKtpd/rOBQKWf+92X2HiwFgTmbN/DCGefTNalueaBO7ZKNYokCzh0zlGF9u5B/4CCDe2SR0NoR9sEgilYmgQj0WBBNitbQRMqrqvj5m2/x3bY8nCI8OOkkphxzdMBxM9esYG9xKeL9M/t800aKyssZnd2De4aN55XcpaTFJfDIqNNb+i0YwoTuWWl0z4pML8hosbnYxZ9yqW3ETwGeA/YAf8NKHdAZq7xxZ+CGlhKwreL2eHhs+ud8uWQ9XbNSefDG0+mcGZmBbx+t+5Hvtll+HW5Vpn7xJRcffVTAxIsrd++pUSwAThwkxlpPqzcOGcmNQ0baur+qsj6/gPjYGLqltUZRVIOhAYzNpQ4+lUttI76IvARkASNqe4eJyBxgC3A+8HbLidn2eOfzFbz1+XIACg+UMvX5T/jHvQ3GRYU9qnX/Yzxq7z9oQHoGs9evq2kPyeiIs5HbXx5Vfvnf9/lkzXoAfjF+FLeMs6eUDIZgY7Z6DmP3P/k878+EeuerM8SdExxxWp99hSXcdfurXHTWkzzy0DtUVFQFHFNV5WZN7k62bS9s8n13FdSt3b5rb5GPnuHPaf37MaxzZwAcItw9bqytdPHXHZfD+B69cIjQN60DT516ZqPvvXBzXo1iAXjmy/kcONSqZYYMBgu1cUQRdl0xyrHSvHwiIs9jJT/rAlxb63pE8vQTH/P94s0AzPl4Jd16ZHDZlb5Lo1ZUVvHr+2ayYlUeInDTz8dz8fkjfPb3xficvvzno6VUVlnul5NGDmyK+GFBfGwsr025iLX5+aTEx9M91d7WVEJsLC+ecz4eVRxNrl2iR7Si7H/YECbYidCPJpuMXeXyAnAX0B14qNb56m+EfwdTqNZkz64Dftv1+XrBelassnYGVWHa9HlcMDmHGGfjtnOG9OnE83+YwjfLNtE1K5VTR0WucgGIdTo5Kju7SWObrlhgRM9unDSgD3NyNwBw89gTSE2Ib/J8BkOTMU82dbCrXH4LVGAZ8GsnqiwBngQeDLJcrcaEk4eQu8YqbeBwCmMn+C/QVP+L0CHS5H3WgT2zGNgzq4mjDWB9/s9cfDZrd+cTHxtLr/TI9DQyRD7GW6wudiP0PcD9IvI4cDTQCdgJLLeTHTOcuWDKCWR3TmXThj0cO7wnQ4/u5rf/iSP7cvywnixauhmHQ7jluok4G7lqqWbmh0uZu/BHumalcusV40hONE/cTUFEGJQdGWnZ2zKqytPPfsZHn64kIz2J391zNv37NW01G4mIp+0lrmwOjQp/9SqSr1pIlpAxZuwAxoy1l+48JsbJXx68iK15hSQlushItxcwWJ/Pvl3L3176AoDvV+dRXFrOn34dsX4RBgNfzF3D27OWArAtr5Cpj77H9OevC7FUrUwUKY9A+Ev/8jmgqnqS97U/VFVPCq5o4YvDIfTsnt6sOXI37anTXlevbTBEGnsL6pZmLiiwX6q5LWC2xerib+UynsN6uPbr+oifawYfHDe4K6/NWlSnbTBEMieO7s+rr8+nuMRyBT/1lKNCLFErE5wyx20Gf8plK4fjTbcSVR9LyzP6uN48/MuzmLdoPV2zU/nZeSeEWiSDoVl07pTKc89cybfzfyQ9PYnxYyPbA7LR2LC5RFOQpb8I/Z4NvTYEj4mjBjBxVONL2xoM4UrnTqlceP7xoRYjNBhX5Do0OZ+1iDiA5Ej3Fgt33G4Pc75aS+mhcsaPGUBq+3aBBxkMhlZHUCRQ6iObqZHaAraUi4icA5wMfK2qM0XkGuApIMFbS/k8VTUW6Rbgwcfe48tvrfxbM95ZxL8e/xnJScZl2WAIN2y5IkcRdgM0bgJuAUREkoC/Y+UVE2Ak8ECLSBflFJeU1SgWgB27DvD9iq0hlMhgMPjE5Barg13lUu328Q0wCiuB5Xzgn1gK5tTgi2ZIiI8lIb5uwaS01MQQSWMwGPxR7Yrs74gmBWNXuWR4f+4ChmB9RM8Bd3rPdwqyXAasgM0H7zqHDqntiI+L5eeXjOGoQV1CLZbBYPCFWbnUYNegXwSkY6V++Yn33I+1xpsc5y3EyOG9eWf6LaEWw2AwBMCOzcWOK7KInIZl03YCz6vqI/Wux2GVlR8OFAAXq+pm77V7gWsAN3Cbqn7sb04ReQ3IASqBhcANqlppQ8yA2F25rPT+XAScCxwClgG9vOfzGhrUGETkryKyVkSWi8jbIpJa69q9IrJeRHJFxGzBGQyGsCMY22Ii4gT+AZwODAYuEZHB9bpdA+xT1b5YiYMf9Y4dDEzB2l06DXhWRJwB5nwNGIhl+kjgcBmVZmNXufwJq2aLeI8nVLUMONt7/ZsgyPIpMFRVjwbWAfeC7w8sCPczGAyG4KEa+AjMCGC9qm5U1QpgBjC5Xp/JwHTv6zeBk8SqzjcZmKGq5aq6CVjvnc/nnKo6W71grVyClirEblbkz0RkINbyaYuqLvFeegcrkeV6n4Ntoqqf1GouAC70vq75wIBNIlL9gc1v7j0NBoMhaNgoFuZduVwpIlfUOjtNVad5X3cBttW6lgfUT99R00dVq0TkAJbZogvWd2ftsdVGWr9zikgscAVwe4B3YBvbQZSquhXYKiLxItJJVXeq6sqAA5vGz4E3vK/9fWB1EJHrgesBunfv3kKiGQzRgcfjYfoTH7Nw7lq69+7IrQ+eR3KqCeL1RSPiXKar6lO+pmngXP0lj68+vs43tENVf85ngXmqGrSs97aVi4iMAh7HWjUAxIjI37GKhz1hR9GIyGdAQwUe7lPVd7197gOqsPYCwd6HbZ20tP80gJycnCjzzTAYgsvsGd8xc9qXAGzO3QXAvU9dFkKJwpzgeITlAbWLSnUFdvjokyciMUAKUBhgrM85ReQPQCZwQ7Olr4XdCP1hwOeAi7pZkN3AlVhv7I5A86jqyQHucyVwFnCSdw8Q7H3YLYbH4+GHBRtRlGNO6NPkwmAGQ6SxdX3dpBtbN5okHP4QVcTjX7vYSMm/COgnIr2A7Vj25kvr9ZmF9b07H8t88LmqqojMAl4XkSeAzkA/LDuK+JpTRK7FilM8yVsUMmjY/aZ8CIgD1tY7/xKW4H6Vhh28rnJ3A+eoammtS7OAKSIS5/1whb6itgAAGANJREFUqj+wFkdVefj21/nttS9w37Uv8sfbXsXjMfkdDNFBztgBSK2y3jk2C+pFK7a8xQKgqlXArcDHwBpgpqquEpGHvGm4AP4NpHvtz78G7vGOXQXMBFYDHwG3qKrb15zeuZ4DsoD5IrJMRO4PyoeB/W2xMVirlTOATbXOr/H+DIaB4xksBfap9w96gare6P1gqz+wKrwfWBDuF5At63czf87qmvZ3X6xlU+4u+gzq3Bq3NxhCyojxA3nguStZPC+Xbn06cualI0MtUngTpEBJVZ0NzK537v5ar8uAi3yMnQpMtTOn93yTkxcHwu7E1ZkSd9Y7n1HvepPx+mz7utbgB9bSxMW7jjgXn3DkOYOhrTJiwiBGTBgUajEiAlurkyiyBNvdFqterfy0+oSIpABPeJsbgylUuNCpWwd+dtvJiAgiwuW3nESXnhmBBxoMhqjD8hZTv0c0YXflMgMr8/F0DuvevVjKSb3X2ySX3DiRcy4bjaqS1D4h1OIYDIZwJrr0h1/sKpdHgPHeo5rqKPm5wF+CJ1L4kZhs6qe0FXZsyudAYQl9hnbDFddi281Ry8ZVeXz13vekZ6dw+uVjcMZEUTING9tidoz6bQW7EfoVIjIJuAQrP00m1srlI+B1rzeCwRDWvP/y1/zz92/i8Sj9jurGX978BfHt4kItVpthS+5Ofn3245QfqgBg3bIt/PpvVwQY1YZwq3X4JXq0i+2gDa9L26uqepmqnqKql6rqy970AyktKaTBEAxefOQ9PN597x9XbOPLd5eGWKK2xeIvVtUoFoBvZi8LoTStj2DqudSmWRGBIpIsIg8Am4MijcHQgjgc4rdtaB6demTWaWf3iDLnl+Akrmwz+FUuInKRiHzqTYU/T0R+Uevar7CUyu+B9i0rpsHQfG588AJiYi0bwNAT+jB+8vAQS9S2GH36MVx+xxlkd09n8PG9ufe5n4dapFalOreYvyOa8GlzEZEpHM7vBVZk/BgR6YGVT+y66q5YxcQMhrDmpAuOZ9jYARTvL6VL744mlU8LcNlvzuSy35wZajFCgqgigVYnUbR48WfQ/wUNJ42sPi9AMfA0VkJLgyHsSctsT1qmWWgbWgAFomx14g9/j25DsT6ue7G2vVKA3wKxWG7I7wC9VPV3qrqvpQU1GAyGcKZ65RLoiBb8KZdk788nVLVEVYs5HJEPcK2qFracaAaDwRBBeACP+j+iCDtxLiOldmrUwwypfV5V5wVPLENzWfbVWooKD3LcuIEkpyaGWhyDoe1jM/NxtGBHuXxZr60NnFebcxlagRf++Db/fcaqGp3dI4OnPrqb9h2SQiyVwdDWseFuHEXKJ5C7jDTiMIQBqsrb0z6vae/aspf5Hy0PoUQGQ3Qgbg14RBP+VhvTW00KQ9AQEdolx1NUUFJzLrF928mNtm9PEXP+u4C4eBen/H97Zx6uV1Xd4ffHJRNDCbNAAqGSMlpRQ1AGBURMAIkDKAiYACWWBrUiUpQ6UWPFqih9rC0lTCqGyBgtoAjYQERmIYRIEiCEEDQJZIBApnuXf+z9cXdOzjfc5Nz7Tet9nvPcc/Z01j73fGfttfb0yUMY4FsgOI1CLTPw20i/lFUuZnZGXwriFMcXLhvLtz89iTdeW8X7TzqYQ449sN4iFcLKFW9w3rHf4c/PLwHgvl8+yrdv/jz5XYKO07fUMhpMbaRdvJ+kBRl59AH8Yvb3WLt6HQO3aJ2W/ayHnn1TsQA8MX02L/95GTvssm0dpXKcEu23xEslXLm0KB0dm9HRQooFYIddByMJiz/gQVsOYKttfCSc0xjU1K/SRrrH179wmoZh++7GhO+cwnY7b8Nbhu3ARVeObynLzGlyDF+4MsEtF6epOG7cezlu3HvrLYbjbEhNCqR9FIwrF8dxnAJQl7vFUly5OI7jFEEbur4q4X0ujuM4RVDLZmE16B5JoyQ9LWmupAtz4gdIuj7GPyBpWBL3pRj+tKQPVitT0rkxzCQVurubKxfHcZwi6AI6rfJRbR6M1AH8CBgN7AecImm/TLKzgKVmthdwKXBJzLsfcDKwPzAK+C9JHVXKnA4cDTy/yfXP4MrFcRynAApabn8kMNfMnjWzNcBkYEwmzRi6V1C5AXh/XER4DDDZzFab2XPA3Fhe2TLN7DEzm7dpNc/HlYvjOE4R1OIWC4yV9HByjE9K2Q14IbleEMPIS2Nm64DlwPYV8tZSZuF4h77jOE4RWA17toToa8zsh2VS5K1llC20XJpy4XlGRK+PPHDl4jiOUwRdXeGoSNVv+gJgaHI9BFhYJs0CSZsTdgl+pUreamUWjrvFHMdxiqCYGfoPAcMl7SmpP6GDfmomzVRgbDw/EbjbwppIU4GT42iyPYHhwIM1llk4rlwcx3GKoOQWq3RU85qFPpRzgV8Ds4ApZjZT0sWSTojJJgHbS5oLnAdcGPPOBKYATwF3ABPMrLNcmQCSPitpAcGaeULSFUU9DneLOY7jFEFXZzgqUt16MbPbgNsyYV9NzlcBJ5XJOxGYWEuZMfwy4LKqQm0EDWe5SDo/ndCjwGVxos8Tkt5Zbxkdx3E2oJrV0tVeM/gbynKRNBT4ADA/CR5N8B0OBw4Gfhz/Oo7jNA6lPhcHaDzL5VLgAta3HccA11rgD8BgSbvURTrHcZyy1DzPpS1oGOUSO6teNLPHM1E1TwCSNL40MWnx4sW9JKnjOE4OnZ3VjzZSMH3qFpP0W+AtOVEXAV8GjsnLlhOW+x8ys8uBywFGjBjRPv9Fx3Hqj7vF1qNPlYuZHZ0XLultwJ7A42GJHIYAj0oaSW2TihzHcepMDa6vNtI9DdGhb2YzgJ1K15LmASPMbImkqcC5kiYTOvKXm9lL9ZHUcRwnH+vsxFRlKHIbWTYNoVyqcBtwLGGFz9eBM+orjuM4Tg6l4cYO0KDKxcyGJecGTKifNI7jODVS1TJpH+XTkMrFcRyn6ejsBKq5xfpEkobAlYvjOE4BmBlWbVVk73NxHMdxekQbTpSshCsXx3GcIugy2srvVQVXLn3A66++wYxpT7HNjn/DPiOH11scx3F6AevqxMyHIpdw5dLLvLZsJZ879CLmz3oRgLHf+ASnfeXEOkvlOE7hmAGbvBNly9Awa4u1Kvfd9MCbigXgum/dhLVR68Vx2oWuri6syyoeXW3023fl0ssM3HLA+tdb9CcuceM4TuswbZG9ANZV9ljV9RpdrANYXm9h+wJXLr3M4Se+m8M+GrafGTCoP+ddcU6dJXIcpxf46RJeYo2tLptgHk8zjH1KE8NbHu9z6WU6Ojr42g3ns3TRcgZtNZCBWwyonslxnKbCzNbur4OYzxz24oAN4lfZ66zgFRbwTEcdxKsLbrn0EdvutI0rFsdpYZ7i4f7lrJfEaqnW498yuHJxHMcpADNbuzvDmc+c9cJLVssT3N82Vgu4cnEcxymMPOulHa0WcOXiOI5TGFnrpV2tFnDl4jiOUyip9dKuVgv4aDHHcZxCKY0cm8MTrGRFW40QS3HLxXEcp2Ce4uH+y3m5ba0WALXqfB5Ji4HnqyTbAVjSB+L0Ba1SF69H49EqdalUjz3MbMe+FKbVaVnlUguSHjazEfWWowhapS5ej8ajVerSKvVoFtwt5jiO4xSOKxfHcRyncNpduVxebwEKpFXq4vVoPFqlLq1Sj6agrftcHMdxnN6h3S0Xx3Ecpxdw5eI4juMUTtsqF0mjJD0taa6kC+stT0+QdKWkRZKeTMK2k3SnpDnx77b1lLEakoZKukfSLEkzJX0uhjdVPQAkDZT0oKTHY12+EcP3lPRArMv1kvrXW9ZakNQh6TFJv4rXTVcPSfMkzZD0R0kPx7Cme7eambZULpI6gB8Bo4H9gFMk7VdfqXrE1cCoTNiFwF1mNhy4K143MuuAL5jZvsC7gQnxf9Bs9QBYDRxlZm8HDgRGSXo3cAlwaazLUuCsOsrYEz4HzEqum7UeR5rZgcnclmZ8t5qWtlQuwEhgrpk9a2ZrgMnAmDrLVDNmNg14JRM8Brgmnl8DfLhPheohZvaSmT0az18lfMx2o8nqAWCB1+Jlv3gYcBRwQwxvirpIGgIcB1wRr0UT1qMMTfduNTPtqlx2A15IrhfEsGZmZzN7CcKHG9ipzvLUjKRhwDuAB2jSekRX0h+BRcCdwDPAMjNbF5M0yzv2A+ACoLQe1vY0Zz0M+I2kRySNj2FN+W41K+26KrJywnxMdh2QtBVwI/DPZrYiNJSbDzPrBA6UNBi4Gdg3L1nfStUzJB0PLDKzRyQdUQrOSdrQ9YgcamYLJe0E3CnpT/UWqN1oV8tlATA0uR4CLKyTLEXxF0m7AMS/i+osT1Uk9SMolp+Z2U0xuOnqkWJmy4DfEfqRBksqNeCa4R07FDhB0jyCq/gogiXTbPXAzBbGv4sIyn4kTf5uNRvtqlweAobHUTD9gZOBqXWWaVOZCoyN52OBW+soS1WiL38SMMvMvp9ENVU9ACTtGC0WJA0Cjib0Id0DnBiTNXxdzOxLZjbEzIYRfhN3m9mpNFk9JG0paevSOXAM8CRN+G41M207Q1/SsYRWWQdwpZlNrLNINSPp58ARhCXE/wJ8DbgFmALsDswHTjKzbKd/wyDpMOBeYAbd/v0vE/pdmqYeAJL+ntBB3EFosE0xs4sl/S3BAtgOeAw4zSzZXL2BiW6x883s+GarR5T35ni5OXCdmU2UtD1N9m41M22rXBzHcZzeo13dYo7jOE4v4srFcRzHKRxXLo7jOE7huHJxHMdxCseVi+M4jlM4rlycTUbS1yVZ6aizLPMSWa6upyy1oMAjUd5OScOTuHHpc43L5DQ1kjaX9Fysz2uSdq23TE7v4MqlF5A0LPNRKB1dklYqLPN/g6SmWSzT6TVOB94Zz280szn1FKYakraWtCq+z3/X0/xxjbLvxsstgX8rVECnYXDl0rcI2AJ4K/Ax4BZJ36qvSC3HROCL8ZhcZ1kqErd+SD+ul9ZLlh5wHDAAeNLMZm9kGVcCy+P5uI1RUk7j064LV/Y1DwPXE5TLMMJ+GANi3AWSvmdmL9dSkKSt4xL1Tg5m9r/1lqEHfIgwWxzgOTO7v57C1MhH49+bKqaqgJm9IekWwhIsmwGfBr5QgGxOA+GWS98w08y+a2b/YWYTgPQD2AGkfvasS22cpI9I+r2kV4EVacGSDpB0uaTZkl6X9EZ0u02S9PasILH8SyVNk/S8pFclrVHY2fJuSeNji3oDJH1Q0r3RtbdU0tS49ElFJPWT9A+SfitpcbzfEkm/kXRSGRmzz+D9ku6StCLe/x6FDbmyeXP7XDalzJj/XZJui2lfjc/qiE3sFxmfnP+iB/mQdKqkdcl9p8TnnFfPUyQ9Gt+N5yV9TXEhSklnS3oyurpekHSJyuw0KWkgYYM9CAuOlsJHSpocy14V77NA0nRJP5A0Iqe4Kcn5WIVFTJ1Wwsz8KPggWCeWHFfHcAF7ECyZNH6PCnmnZa4tSTseWJONT461wDkZ2Y6vkL503AZslsk3jrAGWDbt64T9SzaQL+bbDniwyv2uS++X8wzurXDvfTL3m5d97gWUeRSwKidtJ/DLTNiwGt+R/sAbSb4TctKMyysbOJWwk2cp/Eqgo0w9HyrzzK8krK2XF3dVGZnHxPi5Sdj7CO9Zpf/v13PK2ibz/A+t9+/Wj2IPd4v1DWMljS0T9zMze75C3sOBlwlutUXAQQCSDgF+TLf1uQT4CeGD9ynCRkibAz+SNMPM7ovp1gGPExTcYoLvexBhs67jCQpwNMH9cUO81y6EbaFLe3usJWy1vBj4OGEV4HJcW5KZ8IGeDMwF9gc+EeU/hbBqbbn+p8OAPxFcMQcCx8bwQYQtec+pcP9y1FSmpAGE5zogyTsFmA2cQHhmG8NIYGBy/WAtmSSdSvcimQCXEfbCKTdKbwRwP6EB8Alg7xh+Rvw7Hbgb+CShLxDgU5K+bHFjrYSSS+zGJOwcut3rLwI/BV4FdiVY5IfnCWVmyyXNAUr9Le+LsjitQr21WysebNh6LHfcB2xdJe9yYPece9yYpFlH0tom/KjT1uStOfnfCpwETCD4u88n7HNTyjMpSfsvGZnOSOIGE5TfBpYLcEAm38czMlySxL1M+db3/PQ5AY8mcY9kypyXxF1d4bnWVCZBeab5/j2JG0hQMhtjuYxL8qypIY0BF7G+xfLNGt69mUC/GHdMJm5GEjc6E/ehTLmbE7bWNuDgJPyWJM+FOfIMAnYrU7/fJnmvqvfv1o9iD7dc+oZShz7AzgS3xi6EzZmmSzrSynfoX2tm83PCD0vO7zezN3faM7M5ku4jLMtPvA8AkvYgtMRzW5QJQ5Lzg5LzTkLrtHSvZZJupbslnJK9x/WSrs9JB8F9ti/BgsnyE1t/EMNsgqUFsG2Z8qpRa5lp3SHsQQOAma2SdB1hy4OesmNyXuuy799Mzr9oZt8tm7KbX5jZ2ng+r0Jcdgh09rkeGcMWsL6VNY3gLgP4psLw+tkE6/QR4Hdm9mIZ2dJ3fscyaZwmxZVL3zAz/RBI+m/Cj1nA2wj7mJQbLVNue9btkvM/58SnYemH4ma6P6KVSN1Ag5PzpckHqcRfypSxXZnwcpT7wMzLXKd7iWzsoJRayxycSZd1FeU9+1rYlP2cO4HnakybftjXZOLSHSXXZeKyz/Vj8e/NZpa64H5IaBSMBfoRduBMB0Usl3Smde80mqIy504L4MqlDpjZM5KW0P0xPbJC8pVlwl8h9KsAvCUnPg1bCqAwnyBVLJMJ80EWmlmXpAfZsKUOsCw531ZSv4yC2bmCjCmXEPqGyvFMmfCsMivXv9ATai1zWeZ6J9b/sOc9+1pYnJzXqoRnET7kHcDPJX3UzH5VJU+2nilZhZKLpM3otk7S/hbMrBM4W9IFBKWyN7AXMIrget0GuEbSHWb2eqbotN6+5XCL4UOR64DCTnk7JEEbMwwz7fx8j6R9kvKHs77brJQ2vScEt8iCqFj2BTYYuhx5KDnvAE5L7jWY7g9Plvsy16stDMle7yC42Z4p4/6rN9mO9lNLJ3Fo7ic3sty5yXk/SeUUdMrJhJ0gIbwzN0j6wEbevyccQlCii8n8TyXtLWlLM1tqZreb2Q/M7FxCX1WJrQhKMcvQ5HxuTrzTxLjl0jfsL+n8eL4T4YOUugGyH+Fa+B7w4VhOB3CvpGsJwzs/Rff/1mJaCD/gLrobFT+U9A7Cj38cYXhsHj8h9CsMitf/I+k9dI8Wy215m9kMSbfTPTfiqwrbG/+eMAx3V8JopoMIQ4NvziunzkwluI9Ka2BdLGlvgvUyhmSOUg95kDB6rjRibATwf1XyrAA+SHhWexNcl7dKGm1m/7+RctRCaZTYLdFSSfkMcJake4BnCW7CAUmeEkvTi9goSZ9db8rv1AFXLn3DiHjkMY+NWF/JzKZLmkAYiro5wSo5L5OsE/i8mU2LeRbF/p5/ivFDgH+N508Q3FLvyrnXQkmfIUz+FKHVfHaMXkOwjA7N5oucDtxOt7vtqHg0BWa2WtLphLk/Awj1L1luRqjb6CRLVw/K/R3BfQTBOqimXDCzxZKOITRIhhIU/q8kHWO9N8P/I/FvuVn56eTKPCab2bOZsEPobmC9Qo1DsZ3mwd1ifc86wiiZ3xM+7O8ws4WVs+RjZj8mKINJBMWwitAx/RxhHspBZvafmWyfJQwgeI7gj19ImC/zPuC1CveaRFhXajrB6lhB+LAeRhhSWi7fy4QPyZnArwmd/2sJLdmZhLk0ZxCGRTckZnY3oZ53EJ7RSsIoqWMIVkTKUmrn8uT842VTbSjP/HjvUv/VVsDtkjZoGGwqscxhhL6nu3KSXEWYn3Q3oaG0kvCOLwbuAf6RxI2akNb3GjPLDjZwmhytP/DDcZwssW9ldWaUFHEJlT/Qbe39yczy+hbKldtBcCWV1hc72MwaqgUvaSKhMfJTMzu9oDIHEUbdlWbp72sbvwim06C45eI41TkCmCtpoqTTJB0v6RzCzPfUWvh+TwqN/RdfSYKybs1GYJMXqszhTIJigTB50hVLC+KWi+NUQdIogguwEt83sx6v7CtJhEm27yS04vexBt/TZVOI1t4cgqttJTDcNlxmxmkBXLk4ThUkDQUuIKw4MITQ6l4NvEDoO7uiFzvTHacpceXiOI7jFI73uTiO4ziF48rFcRzHKRxXLo7jOE7huHJxHMdxCseVi+M4jlM4fwU+t5jBSShzlAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "broad_x = [TestSetLabels[i,1] for i in range(len(test_predictions))]\n",
    "vel_error = [(TestSetLabels[i,0]-test_predictions[i,0]) for i in range(len(test_predictions))]\n",
    "# Calculate the point density\n",
    "x = broad_x\n",
    "y = vel_error\n",
    "xy = np.vstack([x,y])\n",
    "z = gaussian_kde(xy)(xy)\n",
    "# Plot\n",
    "plt.scatter(x, y, c=z, s=20, edgecolor='')\n",
    "plt.xlabel('Broadening (km/s)', fontsize = 20, fontweight='bold')\n",
    "plt.ylabel('Residual (km/s)', fontsize = 15, fontweight='bold')\n",
    "cbar = plt.colorbar(extend='both')\n",
    "cbar.set_label('Density', rotation=270, labelpad=10)\n",
    "plt.savefig('Broad_Residual.png')"
   ]
  }
 ],
 "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}