Calibration.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches\n",
    "import seaborn as sns\n",
    "from lifelines import KaplanMeierFitter\n",
    "from lifelines.utils import concordance_index\n",
    "from scipy import stats\n",
    "from scipy.stats import linregress\n",
    "\n",
    "fontsize = 14\n",
    "SMALL_SIZE = 8\n",
    "MEDIUM_SIZE = 10\n",
    "BIGGER_SIZE = 12\n",
    "legend_size =14\n",
    "plt.rc('font', size=MEDIUM_SIZE)  # controls default text sizes\n",
    "plt.rc('axes', titlesize=MEDIUM_SIZE)  # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)  # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=MEDIUM_SIZE)  # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=MEDIUM_SIZE)  # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=MEDIUM_SIZE)  # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "plt.rc('xtick', labelsize=20)\n",
    "plt.rc('ytick', labelsize=20)\n",
    "\n",
    "font = {'family': 'normal',\n",
    "        'weight': 'bold',\n",
    "        'size': 24}\n",
    "\n",
    "plt.rc('font', **font)\n",
    "params = {'legend.fontsize': 'x-large',\n",
    "          # 'figure.figsize': (15, 5),\n",
    "          'axes.labelsize': 'x-large',\n",
    "          'axes.titlesize': 'x-large',\n",
    "          'xtick.labelsize': 'x-large',\n",
    "          'ytick.labelsize': 'x-large'}\n",
    "plt.rcParams.update(params)\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "\n",
    "sns.set_style('white')\n",
    "sns.set_context('paper')\n",
    "sns.set()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data_path:  analysis/sleep/\n",
      "event:  (503,)\n",
      "empirical:  (503,)\n",
      "EVENT RATE:  0.24254473161033796\n"
     ]
    }
   ],
   "source": [
    "dataset = 'sleep'\n",
    "gamma = 3\n",
    "\n",
    "if dataset == 'seer':\n",
    "    time = 'months'    \n",
    "else:  \n",
    "    time = 'days'\n",
    "\n",
    "data_path = 'analysis/{}/'.format(dataset)\n",
    "empirical = np.load(data_path + \"Test_empirical_t.npy\")\n",
    "event  = np.load(data_path + \"Test_empirical_e.npy\")\n",
    "print(\"data_path: \", data_path)\n",
    "print(\"event: \", event.shape)\n",
    "print(\"empirical: \", empirical.shape)\n",
    "print(\"EVENT RATE: \", np.sum(event)/len(event))\n",
    "assert(len(event) == len(empirical))\n",
    "\n",
    "t_empirical_range = np.unique(np.sort(np.append(empirical, [0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(200, 503)\n"
     ]
    }
   ],
   "source": [
    "# load observed & censored\n",
    "sca = np.load(data_path + 'sca/gamma_{}/Test_predicted_t_samples.npy'.format(gamma))\n",
    "\n",
    "print (sca.shape)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "predicted_samples_sca:  (503, 200)\n"
     ]
    }
   ],
   "source": [
    "predicted_samples_sca  = np.transpose(sca)\n",
    "print(\"predicted_samples_sca: \", predicted_samples_sca.shape)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CI Index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SCA CI INDEX:  0.7498264668653791\n"
     ]
    }
   ],
   "source": [
    "def ci_index(predicted_event_times, event, event_times):\n",
    "    return concordance_index(event_times=event_times, predicted_event_times=predicted_event_times.tolist(),\n",
    "                                     event_observed=event)\n",
    "print(\"SCA CI INDEX: \", ci_index(np.median(predicted_samples_sca, axis=1), \n",
    "                                 event=event,\n",
    "                                 event_times=empirical))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RAE "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Relative absolute Error from MTLR\n",
    "def MTLR_RAE(pred_samp, emp_t_i, cens_t):\n",
    "    constant = 1e-8\n",
    "    pred_t_i = np.median(pred_samp, axis=1)\n",
    "    abs_error_i = np.abs(pred_t_i -emp_t_i)\n",
    "    pred_great_empirical = pred_t_i > emp_t_i\n",
    "    min_rea_i = np.minimum(np.divide(abs_error_i, emp_t_i + 1e-8), 1.0)\n",
    "    idx_cond = np.logical_and(cens_t, pred_great_empirical)\n",
    "    min_rea_i[idx_cond] = 0.0\n",
    "    \n",
    "    return np.sum(min_rea_i)/ len(emp_t_i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_rae(samples, method):\n",
    "    observed = event ==1\n",
    "    rae = MTLR_RAE(pred_samp=samples[observed], emp_t_i=empirical[observed], cens_t=1-event[observed])\n",
    "    print(method, rae)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SCA RAE Observed 0.5423789693409443\n"
     ]
    }
   ],
   "source": [
    "print_rae(predicted_samples_sca, \"SCA RAE Observed\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_rae(samples, method):\n",
    "    censored = event == 0\n",
    "    rae = MTLR_RAE(pred_samp=samples[censored], emp_t_i=empirical[censored], cens_t=1-event[censored])\n",
    "    print(method, rae)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SCA RAE Censored 0.2244014265485046\n"
     ]
    }
   ],
   "source": [
    "print_rae(predicted_samples_sca, \"SCA RAE Censored\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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>removed</th>\n",
       "      <th>observed</th>\n",
       "      <th>censored</th>\n",
       "      <th>entrance</th>\n",
       "      <th>at_risk</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>event_at</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.0</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>503</td>\n",
       "      <td>503</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9.0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>503</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23.0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>502</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>64.0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>501</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>74.0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5351.0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5366.0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5467.0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5490.0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5543.0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>452 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "          removed  observed  censored  entrance  at_risk\n",
       "event_at                                                \n",
       "0.0             0         0         0       503      503\n",
       "9.0             1         1         0         0      503\n",
       "23.0            1         1         0         0      502\n",
       "64.0            1         1         0         0      501\n",
       "74.0            1         1         0         0      500\n",
       "...           ...       ...       ...       ...      ...\n",
       "5351.0          1         0         1         0        5\n",
       "5366.0          1         0         1         0        4\n",
       "5467.0          1         0         1         0        3\n",
       "5490.0          1         0         1         0        2\n",
       "5543.0          1         0         1         0        1\n",
       "\n",
       "[452 rows x 5 columns]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEJCAYAAACXCJy4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3deXxU9b34/9c5Z/bJNpMFAgQwQFiKbMoismgX0KrXq2ipVvtVe7GL1P6oC4K7V6stVqtyqxf1cqvXpcC12Cu3bu21gmjdKqCgYCpbIAtJJslk9pnz+2OYIQlJZrJMkpm8n4/HPDBz5sx8Psn4ns+8z+fz/ii6rusIIYTIeGp/N0AIIUTfkIAvhBCDhAR8IYQYJCTgCyHEICEBXwghBgkJ+EIIMUhIwBdCiEHC0N8N6Ex9fTORSNeXCeTnZ1Fb605Bi/qf9C09Sd/STzr2S1UVHA57h8cHdMCPRPRuBfzYuZlK+paepG/pJ9P6JSkdIYQYJCTgCyHEIDGgUzpCiL6l6zr19TUEAj4guXRGdbVKJBJJbcP6wUDul6YZyMrKw2rtOF/fnqQD/vvvv8+VV17JE088wdlnn93q2DvvvMM999xDfX09F198MStXrkRRFHbv3s2qVas4fPgw3/zmN7n77ruxWCxdaqAQou+43Q0oisKQISNQlOQSAAaDSig0MANjTwzUfum6TjAYwOWqAehS0E/4F9V1nVdffZXly5e3e9zv97NixQrGjBnDT37yE9avX8+f/vQnAG688UasViurV69my5YtPPvss0k3TAjR97xeN9nZeUkHe9H3FEXBZDKTl1eI2+3q0rkJ/6pffPEFN954I/Pnz2/3+I4dO2hoaODKK6/kqquuoqioiK1bt1JRUUF5eTlLlixhyZIlTJ48mbfffrtLjRNC9K1IJIymSaY3HRiNJsLhUJfOSRjwi4uLefPNN/nOd77T7vGqqioA8vLyAHA6nVRXV8fvdzgcre5PtXd2HWXVb7fxfx8fJiKl/oXoMkVR+rsJA1Yqtw/p6nN35++U8KM8NzeX3NxcDhw40O7xWCPbvnhv/GLy87O6fE5OTj2fltfyaXktde4gXz+9BJNJZdTQnIx5IxcWZvd3E1JG+ta/qqtVDIaup3O6c04yjhw5wsUXn88NN9zMpZd+F5/Py09/+mN2797N/fev4eabVwDw4IO/Yd68Bezd+wXf//5lANx2212cf/4/9ej1W/brk08+5pln/pOHHnr0pHb1hM/n44kn/o0JEyZyzjnf7tK5qqp26X3V4+9uhYWFANTX1wPgcrmYNGkSRUVF8Z9j/8buS1ZtrbvLCx9OHeXgRxdP4YmXdrJl+1f8/YsqRg3J5hunlzDEacVuMXbp+QaawsJsamqa+rsZKSF963+RSKTLFypTeXEzHI4+bySi4/MFWL16JZ999im3334PFosViAa9Dz74gDlz5vHBB++jqtHZNZGI3qN2te3Xiy8+zz/+UU4oFCE318HDD/8bo0aN7nHfDx48yIsvPsfq1Xd2+bkikUir95WqKp0OlLsV8P1+Px999BGlpaVMnTqVrKwsnnnmGfbt20dlZSXz58+npKSEkSNHsmnTJlRVZdeuXVx//fXdebkuWzR7FJ99WcP+yiaq670oisKh6iaONXixWU50WUFh1NBsrGYNTZWLVEIMZL/61X1s376NlStv41vfOoePP/4QgDFjxvL3v0f/++OPP2Ts2HHs3ftFwud7//33eOSRBzl2rIYZM2Zy882rcTicPP/8M/z+98/hdrspLR3D6tV3sWfPZ7z11l8AuOSSC3jssX9nxYrrWLHiJk45ZQzXX/8jrrrqX/jf//0fjEYjP/jBj9i8eRN7937OBRdcxM9+dgPNzW7uued2PvroA1RV46yzvs4tt9we/5byi1/cDcAZZ8zjV7+6l48++pCioiEsX/7/MWfO3F75HXYr4NfU1HD11Vdz++23c8UVV/Doo49yzz338NFHH3HVVVdx7rnnAvDII4+wevVq7rvvPs4991y+//3v90qjkzFptJMzvjaUF/68j+p6L3/56DBlI/M4tTQ//hiPL8inX9WSbTVRUpRFltWIqmZG2keI3vDOrqNs23m008coCnQngztvSjFnnlqc1GN///vnOXKkAovFwuzZZ7Q6NnXqdF56aSMul4sdOz7hnHPOSxjwXS4Xt922kq9//ZtMn34azz//DI899jArVtzMb3/7KP/8z5cwa9ZsXnjhv3jrrT9zwQUXMXXqdA4dOsgdd/xru8/5/vvv8eMf/5Rf/eo+7r33Dn72sxvIzs5m48YXuPzyK9mx4+8cPHiAlStv4/PPd/P73z/PpZdexrXX/oR77rmdyy+/klmzzuCRRx7k8OFD/PznN7Nz5yfceecqNmx4mdzcvOR+sZ1IOuDPnj2bL7448Uts+d9nnnkmr7322knnTJo0ic2bN/ewiT0zcVT0onF1vReUBk4rO5FWMmVpRCI6jZ5o4B9WYMfYQS6yIMeK2aT1SZuFEK0dOVLBRRddyh//+BK//vUD/PKXD8ePTZ48hT/8YRMbNjyP293EtGnT2bTpxU6f77PPduHxNPPKKy/zyisvA+By1ZOVlcW4cWW88cafcLnqOOec81i06FysVisOh4OqqkqmTJnG0aNHTnrOpUu/xze+8S1eeOFZVFVjyZKlmM1m3nlnK01NTXzzm4txOvP58MP3+fTTXQC43U2MHTsOgNGjSykoKOCDD96joaGBe++9M/7cn3++56QPuu7I+PlX08YWMG1sQXyk/8Kf98WPTRzlYNrYAvKyTPgDYWrqve0+hy8Qpq7Rj8mgMnJINlZzxv/ahADgzFMTj8L7YoHSRRddyg03rMRgMLBx4wv8+c9vxGcAWq02ysrGs2nT73E4nIwcOSrh84VC0emMN998K6WlY/F4mrHbs1AUhd/+9mm2b9/Grl2f8B//8e/8z/9s5qmnnkn4nLFFpYqixv9bVWODRJ0NG17gyScf52c/+zljxlzGnXfuandySygUYv78hXzve1cRCPhRVZXRo0uT+TUlNGgS1xNHOShyWOM/H6p28/oHh3jhz/v45MtjmE0a2TZTuzdnjoVIRKe63ktlnYdgKNzObeCtyBMiU4wePRqAq69eRm5uLo888iBNTScuVk6ZMh2Pp5mpU6cl9XynnjoFk8nMX//6f1RVHeW+++7iued+R0XFYc477xu89dafmTVrDsXFw6mujk4xNxiMuN1NbNvWvfVEH3/8AQAmk4U333wdiF50NRiiE0l27/6UQ4cOctpps9ix4xOOHDnMf//377nhhp8SiYS79ZptZXTAb/npOW1sAZd9Y1z8tmhmCSVFWfHA/8mXxzp8Hk1VMBs1cmwmjtY28/HemnZu1dS4vDL3X4gUysnJ4Qc/+BF1dbWsXfub+P3Tpk0HYOrUGUk9j9OZz/33P0hNTRW/+MXdjBhRwvXX38Dw4SNYseJmPv98D7fdthKPp5nbb49eTP3GN74FKDzxxNputX3p0u+Rm5vLmjW/wGg0omkaBw7sZ8SIEqZPP4033niVXbt2cNNNq5k6dRoPPfRLPv98DytX3o7TmZ/4BZKg6KlcSdBD3ZmWCeDMz+Kt9/fj8YXJtRs7nX//yZfHeP2DQwCUFLU/nSmW+umM2xskEAozYaQDZ07q6gWly/S+7pC+9b/KygMMHZo4JdLSQK05A+DxeNq932w2o2mdX5MbyP2Kafv3Ssm0zIFOUxXGDs9l9/56GpsDhCI6+R0E4Vgg33Ogvt3jh6rdHKp2s+dAfaeBP8tqpK4pTEVNM1azQfL8QgwAixYtaPf+Rx99ghkzTu/j1vS/jI1KNouR0ycUEQpH+GBPFYFQGJOh/U/02IXd9nzy5TH2HKiPB/7Y49uTazNR7fJi0BTGj3KgZsjKXiHS1RNPrG/3/lNOOaWPWzIwZGzAj9FUhUKHjboGX4cBvzOxD4NY6mfPgfoOA76mqRTmWTnW6EM91ICmwciibExGmc4pRH+YPPnU/m7CgJLxAV9RFHLtJo65vDR4/OTazN16nmljC9hzoP6kqZ0xsXSPQVPJtZnw+oM0+0J4fKFWi7mGF9hxZMueAGLg0nU9Y+pOZTJdjwBd+ztlfMAHKMi1YDJqfPZVLW5vsNUxm9mQ9Ora2CKuttqme2IjeqNBI3y8pgeAxxfiy4oGNPXExTuTQWXCKAcGLaMnTIk0YTCYaG5uxG7PnGKDmUbXdcLhEE1N9ZhMXRs8DoqArygKVpNGSWHrq9cNngDeQCjpgmod5fpj6Z7YbJ/YYzRVQVNPpHOMdvWkWUcNngAf7T1RNlpToh8AWdb0LvIm0pPDUUh9fU2XNtaIFSvLNAO5X6qqYbVmkZWV26XzBkXABzAZNUqGtC4jqlc1UdHsJhSOkGvvXqoHTgT49oJ+S4qioGmtR02OrNav2+gJ4A+EsVkMctFX9DlNM1BQkFx9m5h0mXLaVZnYr0GdRxheaOfUMQXoOnj9oR7V8J82toBFM0uAjqd4tkdRlJNunx+s52BlE/5gOH4LBMMp3XxBCJH5Bs0Ivz2aqmIzKwxx2nC5/fgC4R7Nn297YTeZBVtt5dhM+INhjtZ5qGpR20fXdSaOdtK1HQWEEOKEQR3wIboy7ZTiHP5xpIHKOg/+YPs1KywmDYsp8a+rVXVOOp6z3xmzUcPcZipnfZOfw9VNjBrR/oVjIYRIZNAH/JiSomyGOu3tHqtxeag41gxEg3Fnsxc6q84JyZVpaE+2zUh9k5+DlY1oegSDqkq5ZiFEl0jAP85oUDushZ+bZcYfjODxh/D6Q9iSmNXT3hTOlmUaWj4umQ8Ag6aSbTNxuNpNY6MXq8XA1DFd/+AQQgxeEvCTkJdlJi/LzP6jDRyp9SQV8Nubwhkr0xDT1Q8As1EjL8cCoTANngANbj+qqmDQVKndI4RISKJEF+RmWag45sHjD510TCGa508m3RPT9gMg2by/qiqoCuw5WI+uR193zPBcsiyyRaMQomMS8LvAZjZQXGBvdwPPRk+wy7N82n4AtMz7J0r1xNYNRCI6DZ4An/6jjrKSXApyrR2eI4QY3CTgd4HZpFFanNPusa+ONHCw2o2mKd0q0gbdm+GjqgqOLDMNzX4OVbtxZltklC+EaFdSC69efvllFi5cyNy5c3nqqadOOr5+/XrmzZvH/PnzWx3fuHEj48ePj99++MMf9l7LB5iRQ7MpzLPi83d/K7LYrlxFDmt8pB+7dbYjF0C21YTbE6SpTa0gIYSISTjCr6ys5NZbb+Wiiy7C4XCwZs0aZs6cydSpUwH48MMPeeCBB7j22msxm82sWbOG0tJSvv71r7Nz507GjRvHqlWrUBQFp9OZ8g71F01Vycsy09Ac6PFztZ3h07I421mnj2z3HFWNrtL1+kPk2k09boMQIvMkHOG/++67BINBli1bxnXXXYeqqmzdujV+fMeOHQBcc801LF++HIfDwRtvvAHAzp07OXLkCNdeey2//vWvsdvbn+eeKXLsJjRFocnbs6Df3v67EK3V8/7uyg7PM2oKbk9ANlQXQrQrYcCvqoru2J6Xl4fZbMZut1NdfaK6Y3FxtNDStm3b2L17N01NTVRXVxOJRLDb7Zx//vmsXbuWqqoq7rjjjhR1Y2DIsho5ZVgObk+wV+vetKzTs/mv5R2meKwWA5V1XirrPDR5AgRDvbPTvRAiMyRM6cQCV0fTDRcvXsy8efO48cYbyc3NJTs7WpFSVVWef/75+OPOO+88nnvuOSKRCKqaXM22zjbjTaSwMDvxg1IgO9dKnSd4vM6+itZLF1DPOn0kNpuJHftqOHqsGYNBbTe9k50doikQwlXTTMmQHIY7LEmtGxgo+uvv1hekb+kn0/qVMOAXFhYC4HK5MJvNeDweiopOlPDSNI1HH32UiooKhg0bxrnnnsvw4cOpr69n48aNLFiwgAkTJhAOhzEYDEkHe4DaWvdJ9eOT0d9lTZVwmOpjPgKhCDaToddKIJQNy2HWpKE8/t87OFLTzOP/HU2ntZ3CqQLBYJjPy2soP1DLhFEObGbDgN/Qor//bqkkfUs/6dgvVVU6HSgnDPhz5sxB0zTWrVuH0+kkHA4ze/Zstm/fTmlpKcFgkG9+85ssXbqU4cOHU11dzeLFi7FarTz99NO89dZbXHHFFWzZsoXFixf3aucGqrKS6EXXLytcVNd58QVPLNQyGbQer4pteVG3oymcsQJsLrefXeW1jB2Ri91ilBW5Qgxi2l133XVXZw/IyclhxIgRbNiwgZ07d7J8+XKmT5/OxRdfTElJCfPnz8doNLJhwwb27t3LihUrOO+88zAYDEyZMoXXX3+dl156iZkzZ3L33XdjsSS/JZfXG2hvjVNCdrsZj6fns2V6ymo2kJ9joSDXSkGuFYvZyLEGL6qqEI7ohCM6kYiO1oXtDS0WI3k2I6eW5nNqaT5fHW3kULUbu9XIUKft5MebDARCEWobfNQ2+MjLMhEK64QjkQG3reJA+bulgvQt/aRjvxRFwWbreJaeog/gXTXSNaXTEbc3yMGq1u1q9ASxmw0dFm5rKy/Phsvlif8c214RYNHMkk4Xa7nc/nhaR1Vgcmn+gBrxD9S/W2+QvqWfdOxXj1M6ovdkWY1MGt16LcKe/XU0eAIo/mgg1nUdq9lwUj38jrTcXnHPgfpOA35ei+0UXW4/ew+5mDDKgcmgDvj8vhCi5yTg97OxI3IJhU98i6msa6a2wZ90wIcTO211Ra7dRF2Tn4/3VjNpdL4s1hJiEBhYSdxByHj8Im7sZujCLKa2DlW7E5ZgiFEUhfwcC5EI+AIhwpFI0rcBnAUUQnRCRvgDjKopBMNhXM0R8uzmxCccN3GUI15bvys7apmNKvuPNrL/aGNSj9eBYqedkUOyJA0kRJqRgD/ADMu3k201sesftSfV5cmydPzniqV1YqP8ZIN+lrVrqRx/IEzFMTe1jT5OLXVi7GZlUCFE35OAP8AoioLNYjipgNrhGje+QOelEmKj/NgF3O7un9sZs0nDbNKoa/LxyZfHsFtOvhAthBiYJOAPQAZNxZnTer1Ce1M624oF99hIP1Zhs7eDPoAjy4yuQ32TH68/NKCmdwoh2if/l6aJwjwrR2ubcSeoxBnbRSs2Pz82R7+3g76iKCgKqIrC3kMuVEUhL8tEyZDMqj0iRCaRgJ8mLCaNU4blUunyEfQFMWhqpzV6Ws7PT2WKJ9duIhzRCYUjHKn14PaFKMi1UJgnWy0KMdBIwE8TiqKQYzNisZo4Vuumut6bsChbyxRPV7ZN7ApVVVBVBaNBRVXDNDUH8AVCBEJhhhd0v9qpEKL3ScBPIxaTgZLCbKyaQlW9h0hET7h/bSzF03KDdDi5wmZvMBs1TAaVZl+II8c8EvCFGGAk4Kcho0HFqGk0ePyoioKmKgmnV7ac9ZPKC7qKomAzG2jyBjjW4KUgV1I7QgwUEvDTVNnIPMJhnXA4wheHXJgMYYyd1MSJjfSBlF/QVVUFBSU6q6idRbkWs4Esa/psyiJEppCAn6ZyjpdADUciOHPMNHtDRHQdiynxn7QrBde6K9tmpNET5MuKhlb3h8M6WTYDo4fmkJuV/EpiIUTPSS2dNKepKhNHOXFkm3F7Q4lPOG7a2AJKilKXY1cUhVy7ibwsc5ubCa8/zOcH6wlHZLN1IfqSjPAzxNB8G0drm5O6kNtSywu5yerJBV9NU8nLMtPQ7MfrD2O3KFKTR4g+IgE/Q1hN0by4PxhOetVr2/INyeitC7468NlXtUwY6cBiNiAhX4jUk4CfIdTjM3WONXiTDvgtL+Qmq+UF354s5sqzm2n0BPj8oCvafgUcTnuXn0cIkTwJ+BkkN8tEVX1zSl+jNxdz5bTYe7Pe7eejz6sY4bTJDB4hUkQCfgbRVAUFlXBER+tCHr+rUrGYy5FlJhCMUNvglYAvRIrILJ0MYrMYyMs2JSyw1lsmjnJQ5IgurKqu93Z5m8W2rGYDjc3B3miaEKIdSQX8l19+mYULFzJ37lyeeuqpk46vX7+eefPmMX/+/FbHd+/ezYUXXshpp53GypUr8fl8vddycRKLycBQp42IDg3Nfpq8qQ2e08YWcNk3xnHZN8ZR5LDGR/vJbrPYlsVkwO0LsOdAHe4Ut12IwShhwK+srOTWW29lwYIFXHLJJaxZs4YdO3bEj3/44Yc88MADXHTRRSxdupQ1a9bwl7/8BYAbb7wRq9XK6tWr2bJlC88++2zqeiIAcGSbOX18IWOH5+ELhGjy9O1oP7YBS3eCvtGgYjMbOebyUlHjTkErhRjcEgb8d999l2AwyLJly7juuutQVZWtW7fGj8eC/zXXXMPy5ctxOBy88cYbVFRUUF5ezpIlS1iyZAmTJ0/m7bffTl1PBBBd8GQ0aOTYjUwc6SCiR0sXp1pstL9oZglAt9M7ZqOGI9tCXZOPimNuquo9smm6EL0k4UXbqqoqAPLy8jCbzdjtdqqrq+PHi4uLAdi2bRtjxoyhqamJ6urq+HkOR3Sut9PppLy8vEuNy8/v/krQwsLM3Ygj2b4V6zoRTaO2wYvFbEiq7EJPnXX6SPZVNPDVkUb2Hmlk1qShXTo/L88GgNkTwO2PEI6EOKXEiT0DLuTKezL9ZFq/EkaA2Oiqo9WQixcvZt68edx4443k5uaSnZ3d6ryeqK11E4l0/XkKC7Opqel8O8B01dW+FWWbcLmaOXy0IbrK1W5K+crWccNz+epIIx/tqaJsWE7S5+Xl2XC5PK3uO1bv4bN91Ywakt2lFcQDjbwn00869ktVlU4HygkDfmFhIQAulwuz2YzH46GoqCh+XNM0Hn30USoqKhg2bBjnnnsuw4cPjz/G5XLF/215nug7pxTnMKIwm8/21x6/GBrdntBuMaQk+E8bWxCfp99e2YauTN/MzTJzpLYZR3a0Fo8QovsSBvw5c+agaRrr1q3D6XQSDoeZPXs227dvp7S0lGAwyDe/+U2WLl3K8OHDqa6uZvHixZSUlDBy5Eg2bdqEqqrs2rWL66+/vi/6JNowGjSMBhg1JBt/MJrPr6rzEArrGA2pGTV3VLahq6UZzEYNtzdIXaOP3D74diJEJtPuuuuuuzp7QE5ODiNGjGDDhg3s3LmT5cuXM336dC6++GJKSkqYP38+RqORDRs2sHfvXlasWMF5550HwOmnn85rr73GK6+8wqJFi1ixYgUGQ/J5ZK83QHcyQ3a7GU8fzU7paz3pm91qJNduItduotkbpLbRj82Smrz+UKeNU0vzT7rZrUbKjzQSCEU4tTS/1TkWixGf7+TpmKoCdU1+cu0mzMYT2zqmU/CX92T6Scd+KYqCzdbxZkiKPoCnQEgO/2S91bdIROfDL6oJR3TsZgMmY+f74/amF/68j0PVbhbNLGk1ym8vhx/T2BxAb7GbioLChNGOVuUZBjJ5T6afdOxXj3P4IjOpqsKk0U4OV7tp8gT6NOBPHOWIz9ePTd+cOMrBWaeP7PCcHHvrwO5y+ymvaKCsJA+7Jf1n8AjRF6S0wiCWZTWSZTXSjS9RPTJtbAGLZpbEN2DpTlmGXLsJnz/M4erufQsUYjCSEf4gZzZp6LpOQ7O/S+fl2ns2Y6ZlaeZYiuf93ZVJT+NUFIUcu4naRj/Go40MybfJSF+IBCTgD3IFuRacOckH72ZviF3/qKWhOUC21dgrc+NjKZ7Nfy0/Ka/fGU1VsJsNHKltJstmlIAvRAKS0hnkFEVBU9WkbzaLgUmjnVjMGt5AiGCo52UbYike6HpJBrNJw2jQ8AXC1Df5W938gXCP2yZEJpERvugSg6biyDbj8QVxuf24mvzk51p7/LzTxhbESzJ88uWxLtXVt5g0jtY2c7T2xOYvHn+IscNyceRYUrbATIh0IyN80S3DC7OYNNqJyWQgEOydkfTUcdFV3V0e5Rs1cu3mVjdHlpmK2mY++6qWJm8QfyAcvQXDUoxNDFoywhfdpigKeVkmquu9GANhDAYVW5L76bZn1qShfLSnikPV7i6P8tuymKLF4hqaA+zZX3eizSh8rdQp+X4xKMkIX/TIEIeVMcNyGV5oJxAKEw5HejSCjpVkeP2DQz3aTCUmurL4xMg/rOt8+lUtjc3ptYJSiN4gI3zRIzl2Mzl2CIYi1DX58QUjBEJhHN0sdNZyk/Su1t1JhiPLTL3bz8HqJiafkp/4BCEyiIzwRa8wGlROLc1n/Mg8gB6NoHtrM5WOZFuNNHkC7DlQz+79dbiaZOtNMThIwBe9ymYxMGmUg2Av7LI1bWxBfDVubzJo0a0UA4Ew9Y1+jtR68LRTtE2ITCMpHdGrVEUhy2ZE7cWhRG9cxG0rVnXToCnUNfqwWQxkWUJAtM6QI9ssUzlFxpERvuh1iqJgMxupbex5qiR2Ebe30zoxmqaSYzdRVeuhvKKB8ooGvjhYj8sdwOsPpeQ1hegvMsIXvU5VFCaNdvDR3mM0tJPLN2hK0tMiO9o9qyu7ZiViMmiYsk5UC21sDrDnYB05NhNlI/L6tJKoEKkkAV+khEFTmTTKcdIUTX8owt6D9Rg0tdVmJp1pu3tWKmbvtJRjNxEOR6Jz+A/UM7nUGT+m9WauSog+JgFfpESsmmVbwVCYIU4bx1w+jJqaVPG1lpU1AT758li8ln4qAj5wfMN3MzUNXj78vBo4/s3lFFm0JdKXDFdEnzIaNMaNyMNk7P5bLzZ7J3YxN1VUVWGIwxZftBUIRaisbX9HLiHSgQR80U8UGj0Bwt3cvKS3V+Qmw24xUNfkY//RRqnHI9KSpHREvxg7PIeKmmY8/hAmo4pC16ZAtlyRW13vbXVfqsSuORyta8aZYyanh5vACNHXkhrhv/zyyyxcuJC5c+fy1FNPnXR88+bNLFy4kBkzZnDDDTcQCERnZmzcuJHx48fHbz/84Q97t/UibeVmmcnLNmM2aigoNPuCXa66GVuRW+SwxmfxtLz19qhfURQsJgPhsE51vZe6Rh91jT5c7q7tFiZEf0k4wq+srOTWW2/loosuwuFwsGbNGmbOnMnUqVMB8Hq93HbbbXz7299m5syZ3H333cyePZvvfOc77Ny5k3HjxrFq1SoURcHpdCZ4NTGYDCuwM6zADkRH6nzj7A0AACAASURBVI2eAOFgOOnZOzFtZ/EAKR3159rNHGvwcazh+DqD49NQtTYXoI2aKlM6xYCSMOC/++67BINBli1bxpAhQ3jyySfZunVrPOCHw9FR2fjx45kxYwZGoxGTKTo7Y+fOnRw5coRrr72WCRMm8Jvf/CaFXRHpbOzwXBp8IcoP1HU54LedxQMn9snt7RW6EK0blNeiOFyDO8Ce/a0XhkUiOkUOK2OG5/bqawvREwlTOlVVVQDk5eVhNpux2+1UV1fHj2dlZbFq1Sp+9atf8e1vf5sJEyZwwQUXEIlEsNvtnH/++axdu5aqqiruuOOO1PVEpDWjIbritecVeKLaXtRNVZoHIDfLdLwM84mbzWyQVI8YcBKO8GOzETqqK1JXV8fatWv5zne+w5lnnsnKlSt5/PHHWb58Oc8//3z8ceeddx7PPfcckUgENcnFK/n53S+cVViY3e1zB7pM7VuD24/FaiIvz9bj5zrr9JHYbCZ27KuJ3/fVkUYOVbux2UzMmjS0x6/RmUhEp8blYd/RJoD4vxCd7jlxtJPcbpaQHogy9T2Zaf1KGPALC6PbzrlcLsxmMx6Ph6Kiovjxv/3tb9TV1XHppZcyZcoUnn76ad577z2+973vsXHjRhYsWMCECRMIh8MYDIakgz1Aba2bSDem7RUWZlNT05T4gWkok/uWm2cj6A+y/3A9FpOGxdSzSWRlw3IoG5YT/zm2YGvzX8vxeAIpn9VjBILHq3Dm5tlocEXn8Dd6Avg8fjRNJdtqZHhh71cE7UuZ+p5Mx36pqtLpQDnh/1Fz5sxB0zTWrVuH0+kkHA4ze/Zstm/fTmlpKWPHjkVRFNavX8/ChQvZs2cPl112GVarlaeffpq33nqLK664gi1btrB48eJe7ZzILCajxuRTnFTVeTlc4+5xwG8rFuBTvUo3puUqYk1V4j9n20yEQhF8gSDN3iAGTWWIs+ffaoRIRLvrrrvu6uwBOTk5jBgxgg0bNrBz506WL1/O9OnTufjiiykpKeHss89m6NChvPLKK7z66qucddZZrFq1CrvdzpQpU3j99dd56aWX4jN4LBZL0o3zegN0Z32L3W7G48nMLewyvW/BQAijQeVYgxeLUev1EsVDnTYOVruprvfy1dFGPv2qjnBEZ2iKA67FYsR3fLSvKgoGTcWoqQRCkXh5Zq8/hNcfQtd1jIb0md2Tqe/JdOyXoijYbCeXNIkf1wfwkkFJ6ZxssPRt1z+O0ewNEQxHKMix9Grg/+TLY/Fyy7EibItmlqR0xJ+XZ8Plar8sQ2Ozn9jbPBzRycsyUTosF2sPNoTvS5n6nkzHfvU4pSNEf4jtN/vhFzX4AiGs5t4rWNZyGmcsrx9L80Dvll5ORssVu8FQhIbmIJ9+VcuUMQUYtOgHnaoosiGL6DEJ+GJAigW3Yfk2DlQ1YTEZUhLwWpZogBOll/ccqO/zwA/R6anObDN1TX4+OT7DSNfhlGE5DHFInl/0jAR8MaA5si1Uu7z4AuGUpTjajvj7sj5PR5zZJ0b9TZ4gh6vdeHwhTinO6eQsITon1TLFgGazGMixGfEFulZnp7ta1udJdfnlZGVZDRhUlRqXl/qmnm8bKQYvGeGLAc+ZY+VQdTPBUASjoW/GKBNHOThU7e7X3H6MoiiYTRreQIi6Rj+O7ORnugnRkgR8MeBlWY0Mddpwuf2t0jomg5qyC5kd5fZbHutrJoNKszeI2xuM32ezGFDlYq5IkgR8MeAZDSqFDiuRFjOI/cEwzb4QWdbUbTfY0Wye2LG+ZjJqNHmDfLa/FoBQSGfauALZclEkTQK+SAv5ORbyc06kMipq3ByobMKnKb2+Irc9fb1Ktz0GTcXRov5OTYMXjy8kAV8kTS7airRUnB+tpd/sDfXZdoN9tZdusqIrdfvmYrbIDDLCF2lJVRWKHFZqG30EQpGU5vNbGkgXczVVIRTqrYLSYjCQgC/Sls1iJC/bTKM7QKMnSK694xoivaWzhVptpf6DQMHrlxG+SJ4EfJHWxgzL5ViDl/KKhj57zfYWarXVFwu3rGYNl9tPsy+IgoLFpLWq0ClEWxLwRUYIBCM0talsaDZpmFJcdbK97RUhusViqhk0FV2HXeW1RHSdCSMdOHNkjr7omAR8kfbsFiNjhufQ8tKt2xOkvsmPwab226i3ut4bD/wGgxrPt/dmqsdxvASDy+3ncI1bAr7olAR8kfasZgNWc+uSsI3NAZp9IVzNAbKtxj5boRsT21O3rbY5/94K/llWI03eALv312G3GBk1NLO25hO9QwK+yEg5dhMzygr5/EA9Lrc/vho1outkWY0YtNR+ALRN9cTq4betxd9bq3cNmordHK051OwLScAX7ZKALzLaKcNyCLaYunioqomG5gDZNhNaP6R6EtXih+6P+k1GDUNEx+0L0uQJoCgKdktqykqL9CQBX2Q0s1HDbDxx4XaI04YONHtDZNv6d4Vq2yme0PPZPYoS3TVr9/56IrrOjLLCVv0Xg5sEfDGoOHMs+ANhqus9/R7w4eTUzwt/3hdfydudoK8oSrz8QrXLy8GqJsaNyOu19or0JgFfDDoFedFa9w3NAXQdzEZ1wOwf295K3tj9Xf0AyLWZcDUF+Gx/HacMzcYmNXcGvaSuXL388sssXLiQuXPn8tRTT510fPPmzSxcuJAZM2Zwww03EAhE50Pv3r2bCy+8kNNOO42VK1fi88nmDaL/GQ0qE0c7mDjKwbACGx5/kHB4YJQomDa2gEUzSygpOjHrKPYB0NX6PSajismoUt/o52idB68/1NvNFWkmYcCvrKzk1ltvZcGCBVxyySWsWbOGHTt2xI97vV5uu+02Zs+ezapVq3jttdfYvHkzADfeeCNWq5XVq1ezZcsWnn322dT1RIguyLaZyLGbKMyzUpBjpdkXTHxSH4ntuhW7LZpZAtDuit7OKIqC2aiRbTNysLKJ2gYfHl/wpFtogHzYidRLGPDfffddgsEgy5Yt47rrrkNVVbZu3Ro/Hg5Ha3mMHz+eGTNmYDQaMZlMVFRUUF5ezpIlS1iyZAmTJ0/m7bffTl1PhOgGq9kQrbUPNDT7aWj2U1XvIRLpmwqcyehplU6zUcORbaaippld/6hrdfv7vmOybeIgkjBxWVVVBUBeXh5msxm73U51dXX8eFZWFqtWreKee+7hV7/6FTNmzOCCCy6IfwtwOKILUJxOJ+Xl5anogxA94syxMGvCkPjPn3xZQ6MniEGDLGvqC7IlI5bb724tfovJgKWdrjR5Anx1pIkDle5W948qzqYw19rd5ooBKmHAj9Ua72gub11dHWvXruU73/kOZ555JitXruTxxx/njDPO6HHj8vOzEj+oA4WFmbvwRPqWWnPsZrz+EJ8fqMNmNWHqpWmNeXm2bp971ukj2VfRwNFjzWz8a3TgNHVcIbMmDe1xm9p+m2lsDhDSFVy+5HL+rkMnUk0GVWVkcU6/rHFIhYHwfuxNCQN+YWEhAC6XC7PZjMfjoaioKH78b3/7G3V1dVx66aVMmTKFp59+mvfee48LL7wwfl7s35bnJaO21t2tr9aFhdnU1DR1+bx0IH3rG3o4glmBw5UNOHth0/DYStueGDc8l1AoQigUobreSygUoWxYTo/b1lY4HOFIVfJpntwcKw2N0fUDXn+YoD+IyRjNFhs0NW3XAQyk92OyVFXpdKCcMODPmTMHTdNYt24dTqeTcDjM7Nmz2b59O6WlpYwdOxZFUVi/fj0LFy5kz549XHbZZZSUlDBy5Eg2bdqEqqrs2rWL66+/vlc7J0SqGDSVU4bl0NAcoKHZD0BEB5vZ0G8BrOWc/VRW49Q0FXsXSk9k2UyEAtFvA6FIhL0HXaBAKBTBkWNm0mhnqpoquki766677ursATk5OYwYMYINGzawc+dOli9fzvTp07n44ospKSnh7LPPZujQobzyyiu8+uqrnHXWWaxatQqz2czpp5/Oa6+9xiuvvMKiRYtYsWIFBkPy85293ug86a6y28142pTKzRTSt76jKtFdtYrzbRTn2wmGI7jc/m7N2bdYjPh6cSbQp1/Vcajajd1qZKiz+6mi3tCybxajAYtJw2LSMGgKzd4gHn+Q3Cxz2tXqH2jvx2QoioLN1vF1J0Xvqw1Bu0FSOieTvvWfYw1e9h1ykW0zdbn4Wm+kdFqK1eEpKcrism+M67Xn7Y6O+qbrOoFQdJ+Cr52Sj9moAgo2y8BY5JbIQH8/tqfHKR0hRFS2zURelpnaJj/52eZ+LUo2bWxBl+fl97XYOgC/UeOLg/UoioICTC8rwJjijWlE+/q2SLgQacxs1Jg42onVpOEeQAu1Brqc4x+UuXYToXAkWha6qoljDd7+btqgIwFfiC4qGZKNPxCmyROgyRPA7e2/4B/bVas7C7L6g81ioK7Rz+GaZo4ea+7v5gw6ktIRooty7SZKh53YUvFITTPBUOSkueeK0vH6ld4Q21WrNzdSSTWLKRpygqEIbl+AneW1AGgqTBzlTLsLu+lGAr4QXWQ2agwrOHFhrNkbpMnTepSv6xAKR8g7Xqo4FWLTNNtupNKbe+amitGgkm01xWfhNTQHCEd0CfgpJgFfiB4qKzl5/1q3N8jf99X0yeu33Eilpxuo9KWWM53CYZ2/743+vkYOzWKo095fzcpoEvCFSIH4PPTjF3cNnkCripwWk6FXyw/ERvupXJCVSgW50dXMbm+QYFCqd6aKBHwhUkBTFUYUZhEKRXMWjjwr2vH8hccfxOsPkWVNzYYksQu56ZDaiYld61AVhcp6Lw3tLHiK6JBlMVA6LLevm5cxJOALkQKKolBSdKLwVmFhNjXW6P9ulbUe9lc14moOEInoOLN7L8+fjhdyW7JZDATDkXYXXEYiOnVNfkr7oV2ZQqZlCtHHhubbmDVxCNPHFaCq0QuWvSW2eUp3N03pb4qiYDJo7d6Mmoraj4vdMoEEfCH6gaoomAwqE0c5AR1fIBS/9cYOVD3dNGUgUpRobZ6/7ankb3sqeX93FS63v7+blVYk4AvRTxRFIdduYojDis1siFbiNGk09dJCrlh65/UPDqXV4qyOqKpCQa6FHKuJHKuJiK4TDg/YUmADkuTwhehno4tPXIT0+EI0NddS1+gjr4cVJtN1umZnWi5kk+xO10nAF2IAsVkMTC8r5NOvamloDhCL95Zu1uFvOV0zNnunrXSazdOSAnx1tJFD1Z1XtAxHdMYOzyU3hYvg0oUEfCEGGKNBpbQ4h2Aomsv3BsIcrnFHL1p2c8QfS++0FZvN0/bibjp8CGTbTISTKJ/e7AvEf5eDnQR8IQaglqPRJk+A2kYDLrcfZ073tltsuVtWS598eeykYN/Rh0BHDAaV0PGA2pcfFIqiYNASfwAapFxDnAR8IQa4bJuJr4128Pd9x3B7g726YKu9D4L2PgSSkSnXCTKZBHwh0oDRoFFanMPewy4gNSt0Yzr6NtCR2I5X7V0nGCipoSZPkFAXZ/T4dYX6ut7bpawrLGYtJYX3JOALkSZMRo3onlEDU9vrBANlxG+zGKlxdX2zFZcvRGM/bNISCIVxZFsk4AsxmBkNKkaDiuv4ytxQKEKOzYipG7N3UqHtN4MX/rwvvvCrP4N+9PfW8cbeHcnLMkMonIIWdc4fDKesTHRSC69efvllFi5cyNy5c3nqqadaHXvppZcYP358q9uVV14JwMaNG1vd/8Mf/rD3eyDEIGE1G5gxvpBZE4uYNbGIglwL/gFcWbLtwq9MWPyV7hKO8CsrK7n11lu56KKLcDgcrFmzhpkzZzJ16lQA5s2bx/r16wE4dOgQd9xxB0uXLgVg586djBs3jlWrVqEoCk6nM4VdESLztawlY7caB3RpgZYLv2DgpHgGs4QB/9133yUYDLJs2TKGDBnCk08+ydatW+MBv6ioiKKiIgD+5V/+hbPPPpvzzz8fiAb8I0eOcO211zJhwgR+85vfpLArQgwuZqMa3zFqoGqZ5ml7UXegXNAdTBKmdKqqqgDIy8vDbDZjt9uprq4+6XG7du1i69atXH/99QBEIhHsdjvnn38+a9eupaqqijvuuKOXmy/E4GXQVAZ4vG9l4igHRQ4rEB3tp1slz0yQcISvHx9CJNqM+cUXX2TKlClMmjQJAFVVef755+PHzzvvPJ577jkikQiqmlzNtvz8rMQP6kBhYXbiB6Up6Vt66u2+aWYjB495yMuz9erzdkcybTjr9JGcdfpIAJ58eRdHjzWz8a/lAEwdV8isSUNT2sbu6I/frS8QwmzSUvL/QsKAX1hYCIDL5cJsNuPxeOIpnJa2bt3Kd7/73fjP9fX1bNy4kQULFjBhwgTC4TAGgyHpYA9QW+tudyOExG3Opqam8/oa6Ur6lp5S0Td/MIzPG6C+PvGALJVi8/C7YtzwXEKhCKFQhOp6L6FQhLJhOSlqYfd0p1+9wR8MYzJq3Xq/qKrS6UA5YcCfM2cOmqaxbt06nE4n4XCY2bNns337dkpLSxk6dCiVlZVUVVUxZcqU+HlWq5Wnn36at956iyuuuIItW7awePHiLndACNE+s1HDYtJwNQdQAKOmYk/Rtom9TXL7/SPhcHvEiBHcf//9bNu2jY0bN3LTTTdRXFzM1VdfzZtvvgkQz+kXFxfHz7NYLDz66KM0NDRwyy23MHXqVFatWpWibggxOE0dW8DsSUMYOSQLj7936uj3Ncnt9x1F1wfudX5J6ZxM+paeUt03tzfIji+PkW0zoioKBq3v9jbqzdRHbJR/2TfG9crz9UR/p3Q6qnDamR6ndIQQA5/VrOHINhMMRahz+ynKs/Z3k7qto7r9IOmenpKAL0QG0FSVSaOdhMIRPthT1d/N6bbORrWycKvnJOALkUFURQEFPP4QCmAxaf06g6erOqvU2dGoXyRPAr4QGURRYHhBFqFwhGMuHwZNxWhIn4CfiMzm6RkJ+EJkEEVRGDkkumCn2RvC7Q3iyM6MvVxbpnskvdM9EvCFyFCji7PZWV5LkyeA1Wzo05k7qdB27r7ouvR+BwghOmQ1Gygdlo3VbMAX6Pu67mLgkRG+EBnKoKkUH5+Tvb/STZMnQLat6xuBDFSdTd9sSXL9J0jAFyLDDXXaMRo0du+vJxw5uX6+zWwYMLtmJSvZRUmS629NAr4QGU5RFLJtJia1EyQ9/hBHjjVjNKgnnTOQJbvRuuT6W5OAL8QgYDZqmNtZfVvf5EdVPLi9J+rw6AC6To49M2b3iBMk4AsxiDmyzcyaNKTVfXWNPj79qg69ObntE3VVpfH4YxUg22Ya8N8QBisJ+EKIVrKsRr42Ovn9p535dupqmwHYf7SRUFjPqMVemUQCvhCiFZNRIz83+Yu4hQ4bSig67bOyrpmG5iCaqhCO6CiQMQu/MoEEfCFEr5k4yknkeMV1XyDMzi+PEQpH+nXRV7LTN9syGFRCoUgKWtS5iaMcTE3RrCIJ+EKIXqOqCirRdI7ZqJNjN9HkCZBlNZ00E6gvdKemfCaTgC+ESAmjQWNyaT47yo/h8QUxGlQURcFq7ruwk+z0zfb05wYoqSIBXwiRUkOdVtyeEIoCNS4vRoOa9nV90pUEfCFESg1x2BlyPLPS4A7g8YUwaAqapmJOsxW+6U4CvhCizxQX2PD4woRCEY41eiXg9zEJ+EKIPjPUaQfAFwhxrNFL8PgsGEVB0jx9IKnf8Msvv8zChQuZO3cuTz31VKtjL730EuPHj291u/LKKwHYvXs3F154IaeddhorV67E5/P1fg+EEGlHU1VybCZ0oqUcGpoDhCN6fzcr4yUM+JWVldx6660sWLCASy65hDVr1rBjx4748Xnz5rF+/XrWr1/PPffcA8DSpUsBuPHGG7FaraxevZotW7bw7LPPpqgbQoh0YjSoTC7NZ8qY6M1qNqDrEvBTLWHAf/fddwkGgyxbtozrrrsOVVXZunVr/HhRURFz585l7ty5vPHGG5x99tmcf/75VFRUUF5ezpIlS1iyZAmTJ0/m7bffTmlnhBDpSVGgyROkoTlAsy+Y+ATRLQlz+FVVVQDk5eVhNpux2+1UV1ef9Lhdu3axdetW/vCHP7Q6z+GIXp53Op2Ul5d3qXH5xzdv6I7CwuxunzvQSd/Sk/StY/OyLYQjOvWNPvYfbSQv29JLLeuZvDxbn7+mLxDCbNJS8n5JGPBjX7MSVb978cUXmTJlCpMmTWp1Xk/U1rqJdCOvV1iYTU1NU49ffyCSvqUn6VtyAt4AjY0+lHDflzRoqz8XXpmMWrd+p6qqdDpQThjwCwsLAXC5XJjNZjweD0VFRSc9buvWrXz3u9+N/xx7jMvliv/b3nlCCBGnKET0CA3tlGbWdYjo4JRibN2WMODPmTMHTdNYt24dTqeTcDjM7Nmz2b59O6WlpQwdOpTKykqqqqqYMmVK/LySkhJGjhzJpk2bUFWVXbt2cf3116e0M0KI9JZjM3JaWfsDw3BE55Mvj/VxizJLwou2I0aM4P7772fbtm1s3LiRm266ieLiYq6++mrefPNNgHhOv7i4uNW5jzzyCD6fj/vuu49zzz2X73//+ynoghAiUyiKgsmotXszGlRU2VilRxR9AM+Fkhz+yaRv6Un61nOhcISPvqgh125K+WtB/+fwu1Pps8c5fCGEGCjCkUir/XctJk1W6HaBBHwhRFpQVYVRQ7Lj3/obPAE83hA5fTTizwQS8IUQaUFVFIYXnkhXmOo9HPBmZposVeS7kBBCDBIS8IUQYpCQlI4QIm35AqGk9so1GVQ0ubgrAV8IkZ7sFiPFBfaEZV88viBuX6jPpnMOZBLwhRBpKctqZNyIvISPO+byUl7R0ActGvjkO44QQgwSEvCFEGKQkIAvhBCDhAR8IYQYJCTgCyHEICEBXwghBgkJ+EIIMUhIwBdCiEFCFl4JITKeP9T+Prmd0TWVxi6e0xsCoQj5OZaUPLcEfCFERsuxmzi1NL/L5+Xn26mtbU5BixIzGVOTfJGAL4TIaLE9cbsqP9dKJBBKQYv6j+TwhRBikEgq4L/88sssXLiQuXPn8tRTT510vLy8nMsvv5zp06fzgx/8gKam6C407777LuPHj4/fzj///N5tvRBCiKQlDPiVlZXceuutLFiwgEsuuYQ1a9awY8eO+PFwOMyyZcvQNI0777yTjz/+mHXr1gGwY8cO8vLy+I//+A/Wr1/Pvffem7qeCCGE6FTCHP67775LMBhk2bJlDBkyhCeffJKtW7cydepUAHbu3ElFRQX33nsvZ5xxBlOnTsXhcMSP+f1+fvSjHzF8+HAefPDB1PZGCCFEhxKO8KuqqgDIy8vDbDZjt9uprq6OH6+oqADgmWeeYerUqdx88800NjYCYDAYOOuss/i3f/s3NE3jhhtuQNf1VPRDCCFEAglH+LEA3dGuMrHjVquVxx57jNWrV3Pvvfeybt06Hn300fjj9u/fz3333Ud1dTVDhgxJqnH5+VmJH9SBwsLsbp870Enf0pP0Lf1kWr8SBvzCwkIAXC4XZrMZj8dDUVFR/PjQoUMBuOCCC1i4cCGzZs1ix44dBAIB/uu//ovJkycza9YswuFw9AUNyc8Era9vJhLp+jeC/PwsamvdXT4vHUjf0pP0Lf2kY79UVcHhsHd4PGH0nTNnDpqmsW7dOpxOJ+FwmNmzZ7N9+3ZKS0uZOnUqBQUF/O53vyMYDPLhhx8ybdo0TCYTGzduZOPGjfzsZz/jxRdfZObMmeTnJ78AorOGJ9KTbwcDnfQtPUnf0k+m9SthDn/EiBHcf//9bNu2jY0bN3LTTTdRXFzM1VdfzZtvvonJZOLJJ5/E7XZzyy23UFZWxm233QbAww8/TFZWFjfffDP5+fk88MADKe+QEEKI9im6XEUVQohBQVbaCiHEICEBXwghBgkJ+EIIMUhIwBdCiEFCAr4QQgwSEvCFEGKQkIAvhBCDRMYF/ES1+wey999/n/Hjx/N///d/ALzzzjssXryYWbNm8cADD8TrFu3evZsLL7yQ0047jZUrV+Lz+QA4fPhwfF+CH/3oR9TX1/dbXyBaZ+mXv/wlc+bMYfbs2dx7771EIpG07xdAJBJhzZo1nHnmmcybN4///M//BNL/b9bWDTfcwPjx44HM6dtNN93Uap+Op59+OmP6lpCeQY4ePap/7Wtf02+77Tb917/+tV5WVqZ/8skn/d2shCKRiP6nP/1Jnzlzpl5WVqb/5S9/0X0+nz5z5kz9xz/+sb5+/Xq9rKxM37Jli67run7uuefqS5cu1Tdt2qR/7Wtf09etW6fruq7/4Ac/0BcvXqz/8Y9/1GfNmqXfeeed/dgrXX/llVf0srIy/emnn9afeOIJvaysTP/DH/6Q9v3SdV3fsmWLXlZWpv/ud7/Tb775Zr2srEzfs2dPRvQtJtbHsrKyjHg/xixatEj/+c9/rm/btk1/55139IqKiozpWyIZNcJvWbv/uuuuQ1VVtm7d2t/NSuiLL77gxhtvZP78+fH7duzYQUNDA1deeSVXXXUVRUVFbN26lYqKCsrLy1myZAlLlixh8uTJvP322wSDQd577z3OOeccLrjgAubNm8fbb7/dj72KFtZbsWIF11xzDd/73veA6Mgo3fsF8K1vfYvXXnuN7373u5SUlKAoCo2NjRnRN4iWRb/77rsZN24ckBnvR4CGhgYOHDjA22+/zQ9/+EM2btyYMe/JZGTUJuaJavcPVMXFxbz55pscOHCAV155BWjdFwCn00l1dXX8/tgmM06nk/Lycurq6ggGgyc9vj+ddtppnHbaaQD89re/BUDToptJp3O/AIxGI6NHj+bee+/l2Wef5YorrsiIvxlEU3GrV69m9uzZjBs3jn379mVM36qqqigrK+Of//mfyc7O5vbbb6esrAxI/74lI6MCvp6gdv9AlZubS25uLgcOHIjf11Ff9ASljwZi39euXcvTTz/N0qVLGT58OJAZ/QK49NJLGTVq5FnIzAAABRpJREFUFPfeey+TJ08G0r9vzz33HJ999hmbN2/mhRdeACAUCgHp37eysjL++Mc/xn9+/vnnWbt2LZD+fUtGRqV0WtbuDwQCJ9XuTyexvsQuBrlcLoqKiuL9cblcre53OBxomnbS/f1t3bp1PPbYY1x88cXcddddGdOv8vJyXn/9dcaPH8+VV16J1WrlySefBNK/b6+99hr19fUsXLiQJ554AoBHHnkESP++ff7556xduxaPxwNEP8hie3Ske9+SkVEj/PZq9y9YsKC/m9UtU6dOJSsri2eeeYZ9+/ZRWVnJ/PnzKSkpYeTIkWzatAlVVdm1axfXX389JpOJmTNn8r//+7+MGTOGrVu3cs455/RrHz744AMeeughSktLOf/883nvvfdwOp1p3y+I7td8yy23sGrVKtxuN16vlx/96Efcfffdad+3O++8E7c7uvHHxo0b2bRpE7/+9a+59tpr075vPp+Pxx57jJqaGsaMGcPevXu54447eOihh9K+b0npv+vFqbF582b9rLPO0s844wz9ySef7O/mdMl7770Xn6Wj67q+bds2fdGiRfrMmTP1X/ziF3okEtF1Xdc/++wz/cILL9RnzJih33TTTbrP59N1XdcPHjyoX3755fq0adP0a6+9Vq+rq+u3vui6rv/sZz+Lz/KI3R544IG075euR2dWPfDAA/qsWbP0M844Q3/sscd0XU//v1lbjz76qF5WVqbreub07ZlnntHPPPNMffr06foDDzygRyKRjOlbIlIPXwghBomMyuELIYTomAR8IYQYJCTgCyHEICEBXwghBgkJ+EIIMUhIwBeDwjXXXENdXR3Lli3jyy+/7JXnfOyxx7jnnnsAevV5hUiVjFp4JURH3nnnHYD4atjelqrnFaI3ScAXGW/VqlUA/L//9//48ssv2bBhAx6Ph4ceeoiioiL27duH1Wrlpz/9Kc8++yxfffUVixYtYvXq1QD85S9/4fHHHycYDGKxWFi5ciXTp09v9Rpf//rXeeSRR/B4PDz88MOUlJSwb98+AoEAd9xxB3PmzCEQCPDggw/ywQcfEA6HmTRpErfddhtZWVl9/jsRg5OkdETGu//++wH43e9+R3Fxcfz+Xbt28eMf/5hXX32V/Px81q1bx7//+7/z0ksv8fzzz1NVVcX+/ft5+OGHWbduHZs3b+Zf//Vf+elPfxqvxdKenTt3cs0117B582YuueSSeHGudevWoWkaL730En/84x8pKiriwQcfTG3nhWhBRvhi0BoxYgSTJk0CYOTIkWRnZ2MymXA6ndjtdhoaGvjggw+orq7mqquuip+nKAoHDx7s8HmHDRvGxIkTAZg0aRJ/+MMfAHjrrbdoampi+/btAASDQfLz81PUOyFOJgFfDFomk6nVz7GqiS1FIhHOOOMMfvOb38TvO3r0KEVFRbzxxhvtPq/FYon/t6Io8TK7kUiE1atXs3DhQgCam5vx+/097ocQyZKUjhgUNE2L13Tvijlz5vDOO+9QXl4OwF//+lf+6Z/+qVuBet68eTz33HMEAgEikQi33347Dz30UJefR4jukhG+GBS+9a1vcfnll9Pc3Nyl88aNG8c999zDz3/+c3Rdx2Aw8Pjjj2Oz2brchp/85Cf88pe/5KKLLiIcDjNx4kRuueWWLj+PEN0l1TKFEGKQkJSOEEIMEhLwhRBikJCAL4QQg4QEfCGEGCQk4AshxCAhAV8IIQYJCfhCCDFISMAXQohB4v8H/Nl29twj8T0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# \n",
    "kmf = KaplanMeierFitter()\n",
    "kmf.fit(empirical, event_observed=event)  \n",
    "kmf.median_\n",
    "kmf.plot()\n",
    "kmf.event_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "range_quant = np.arange(start=0, stop=1.010, step=0.010)\n",
    "t_empirical_range = np.unique(np.sort(np.append(empirical, [0])))\n",
    "km_pred_alive_prob = [kmf.predict(i) for i in t_empirical_range]\n",
    "empirical_dead = 1 - np.array(km_pred_alive_prob)\n",
    "#print(\"empirical_dead: \", empirical_dead, \"size: \", len(empirical_dead))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bounds\n",
    "def ci_bounds(surv_t,cumulative_sq_,  alpha=0.95):\n",
    "    #print(\"surv_t: \", surv_t, \"cumulative_sq_: \", cumulative_sq_)\n",
    "    # This method calculates confidence intervals using the exponential Greenwood formula.\n",
    "    # See https://www.math.wustl.edu/%7Esawyer/handouts/greenwood.pdf\n",
    "    #alpha = 0.95 \n",
    "    if surv_t > 0.999:\n",
    "        surv_t = 1\n",
    "        cumulative_sq_ = 0\n",
    "    alpha = 0.95\n",
    "    constant = 1e-8\n",
    "    alpha2 = stats.norm.ppf((1. + alpha) / 2.)\n",
    "    v = np.log(surv_t)\n",
    "    left_ci = np.log(-v)\n",
    "    right_ci = alpha2 * np.sqrt(cumulative_sq_) * 1/v\n",
    "\n",
    "    c_plus = left_ci + right_ci\n",
    "    c_neg = left_ci - right_ci\n",
    "    \n",
    "    ci_lower = np.exp(-np.exp(c_plus))\n",
    "    ci_upper = np.exp(-np.exp(c_neg))\n",
    "       \n",
    "    return [ci_lower, ci_upper]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    " #Population wise cdf \n",
    "def compute_km_dist(predicted_samples, t_empirical_range, observed=None):\n",
    "    km_dead = []\n",
    "    km_surv = 1\n",
    "    \n",
    "    km_var = []\n",
    "    km_ci = []\n",
    "    km_sum = 0\n",
    "    \n",
    "    kernel = []\n",
    "    e_event = event\n",
    "    \n",
    "    if observed:\n",
    "        predicted_samples = predicted_samples[event==1]\n",
    "        e_event = e_event[event_test==1]\n",
    "        \n",
    "            \n",
    "    for j in np.arange(len(t_empirical_range)):\n",
    "        r = t_empirical_range[j]\n",
    "        low = 0 if j ==0 else t_empirical_range[j-1]\n",
    "        area = 0\n",
    "        censored = 0\n",
    "        dead = 0\n",
    "        at_risk = len(predicted_samples)\n",
    "        count_death = 0\n",
    "        for i in np.arange(len(predicted_samples)):\n",
    "            e = e_event[i]\n",
    "            if len(kernel)!= len(predicted_samples):\n",
    "                kernel_i = stats.gaussian_kde(predicted_samples[i]) \n",
    "                kernel.append(kernel_i)\n",
    "            else:\n",
    "                kernel_i = kernel[i]\n",
    "            at_risk = at_risk - kernel_i.integrate_box_1d(low=0, high=low)\n",
    "            \n",
    "            if e == 1:\n",
    "                count_death+= kernel_i.integrate_box_1d(low=low, high=r)\n",
    "        if at_risk == 0:\n",
    "            break\n",
    "        km_int_surv = 1 - count_death/at_risk\n",
    "        km_int_sum = count_death / (at_risk * (at_risk - count_death))\n",
    "        \n",
    "        km_surv = km_surv * km_int_surv\n",
    "        km_sum = km_sum + km_int_sum\n",
    "        \n",
    "        km_ci.append(ci_bounds(cumulative_sq_=km_sum, surv_t=km_surv))\n",
    "    \n",
    "    \n",
    "        km_dead.append(1-km_surv)\n",
    "        km_var.append(km_surv * km_surv * km_sum)\n",
    "      \n",
    "    return np.array(km_dead), np.array(km_var), np.array(km_ci)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "km_dead_dist_sca, km_var_dist_sca, km_dist_ci_sca = compute_km_dist(predicted_samples_sca, \n",
    "                                                                    t_empirical_range=t_empirical_range)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1111ff518>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAEUCAYAAABnHY4tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOzdd1yV5fvA8Q8HBBkagoAmSqKCFblwhXuF29ypOFBxa5krcaVpmiscWeLWcKEl5t65za+amqvAkQtwgIIg6zy/P/x5ElmHcQ4cuN6vl6/Oecb9XOcJL2+ucz/3baQoioIQQgidU+V2AEIIUVBIwhVCCD2RhCuEEHoiCVcIIfREEq4QQuiJJFwhhNATSbhCCKEnJrkdgD5ERLxArdZ+uLGtrRVPnkTrMKKcZUjxGlKsYFjxGlKsYFjxahOrSmVEsWKW6R5TIBKuWq1kKuG+PseQGFK8hhQrGFa8hhQrGFa8ORGrlBSEEEJPJOEKIYSeSMIVQgg9kYQrhBB6UiC+NEtPbOwLoqMjSUpK1GwLD1ehVqtzMarMMaR4czJWY2MTrKysMTdP/5thIfKKAp1wY2NfEBUVgbW1HYUKmWJkZASAiYmKxETDSGBgWPHmVKyKopCQEE9k5CMASbrCIBTohBsdHYm1tR2mpma5HYrIJCMjI0xNzbC2tuPZs8eScAUP9+8h5t7dNPfHxMRw/PjvNGzYBFNT0wzbs3AsTclmzXMyxIKdcJOSEilUKOMbL/KuQoVMk5WDRMGVXrINfxTOtm1befbsGYoCnp4t9BjZfwp0wgU0ZQRhmOT/n3hbOW+fZO/VajX9m9TjSnAw1aq587XfD5QoUTJXYpNRCkKIfE2lUvHTTyvo1asv27btzrVkC5JwDVLdutWT/aldu5rm9ebN63V23RUrllK3bnXWrVud5jEPHz6gbt3qdO/eMUeuOWPG19StW50DB/bmSHuiYIiJiWHjxgDNe1fXisyd60fhwoVzMSopKRgslUpFnTr1gFe/Vr9eC/Tddx11dk0np/eoV68BZcqUSfOYwoULU69eA+zs7HUWhxDpuX//Hr16dePy5YskJibi5dU7t0PSkIRroExMTJg5c97/v9bPsLCmTT1p2tQz3WOKFbPRxCWEvt1/cJ/BnzTk0aNw3nuvLDVq1MrtkJKRkkI+NWzYAOrWrc6pU8fp1asrTZrUYdq0STx6FM7YsV/QuHEd+vTpTnDwP8B/pYAvvxzOpk0BtGzZhNatm/Hzz6s1bb5dUnh9jXXrVtGqVRMGD+6Xaknhzp3bfPnlMJo0qUPz5o2ZPn0K0dH/TXW3Zs0KOnZsTaNGH9OuXXN++GGBwTzIIfKOv/66zKZN63n0KJx69RqwZ88hXF0r5nZYyUgPNxU2NlZp7ps7dwG9enkDsHbtKkaP/jzNY8PDn2teN21an0uX/szwOG0lJiYyfvwoIHlJYfLk6Zibm2uOmzzZl4oV3+fBg/vs27ebkyePUby4Hfb29gQH/828eTP58ceVmuP/+usSV65cwsWlIn/9dYmfflpMqVKONGrUNM1YVq70p1KlKri4pPzhjot7yZdfDiMsLJSKFT8gNjaGPXt2Ehsbw4wZc9i3bw/Llv2ItXUxqlZ158qVy2zYsA43t49o0KBxpu+LKHgSExOZOnUSj/bsBKB//4FMnfothQoVyuXIUpKEa6DUajXHjv2eYntiYvIxqZ991oN+/QayevVyli//CWtrG1av3sDTp0/o0KEVN2+GJDs+NjaGZcvWUrHi++zbt5tp0yaxdevmdBNujx696d9/EPCqp/ymffv2EBYWSrNmzZkyZTqKksSoUV9QtKg1iYmJlC5dGm9vH1q0aM2775Zi48afWbzYj7t3/83qrREFwJsPOcTFx5F05AAqlYpmTT+h/bdzcjm6tEnCTcXTp9Fa1UR79fLW9HYzcuDA0eyGlYypqSmHDp0E0q/hVq3qDkDx4nYAfPihGyYmJpr3CQkJyY63s7OnYsX3Aahbtz4Ad+/eSTeWDz/8KM19t2/fBMDdvQYAhQoVYv78RZr977//IS9fvmTr1k1cuvQnN25cByA+Pj7da4qC7c2HHMxMzejQvhPRL17gWtsjF6PKmCTcfO71Y8sq1aty/ethMa/fv+3NBPy6TGFklH6p39Iy48dqk5KSNK/j4l5iZvYqjtc9Wnf3GrRv35nw8DCWL/9Jc20hUnPr9i1u37pJvzUbDOrhF/nSTCQTEfGUixcvAHD69CkAypRxSvcclco4zX1OTmUBOHv2NACJiQl4eXWlXbvmPH/+jJ07twPw+eejaNmyDY8fP872ZxD5l6IoLF36A1u3buZ/586yd+/u3A4pU6SHa6DS+tKsWrUadO78WZbbNTY2ZvToEbz//of89ddlALp06Zbl9po1a87Klf4cOXKIfv16Ehf3kocP71OnTj2KFn2HEiVKcuvWTcaM+YLixe24cuXVNWNjY7N8TZE/xcXFMXbsSDZs+JkOdnZ8XNuDTz7J2clldE0SroFK60szK6si2Wq3VClH2rRpz7p1q7CwsMDHZzB16zbIcnvm5ubMn7+IBQvmcfnyJSwtLWjRojXDh38JwIgRo4iOjubvv69jamrKwIHDWLp0MVeuXMrW5xD5S3h4ON7ePTh79gwNixenbZt2uLq+n2ZpLK8yUgpAsezJk+hUV9wMDb1DiRIpf102pPllIWfiffjwAZ07t6VMGSfWr9+aQ5GlpIt7m9b/x5xgZ1eER4+idNJ2TjOkWEH7eIOD/6FTp7Y8eHCfd98txZLPuuPgUEIn0yemRZtYVSojbG3THlIKUsMVQuRxJUqU5J13rKlevSZ79x7BwaEEgN6SbU6SkoIQIs9Rq9UkJCRgZmaGlZUVmzf/irV1MczMzIjO+PQ8SxKuAKBkyXc5fvx/uR2GEERHRzNs2EDMzc1ZsmQZoQf2EnPvrkEn2tck4Qoh8ox//71Dr17duHr1L4oWfYc7d26TlMpKDhaOpXMhuuyThCuEyBNOnz6Jt3cPnjx5Qrly5Vm3bhPvvVeW1w+fv72SgyGSL82EELlu3brVdOjQmidPntCoURN27z5I+fIVeLh/T26HlqP0mnCDgoJo0KABHh4eLF++PMX+VatWUbduXdzd3Rk9erRm8HtgYCCurq6aPwMHDtRn2EIIHdq0aT2jRo0gMTGRQYOGERAQiLV1MeC/ORMMtYTwNr2VFEJDQ5kwYQLt27enWLFizJkzhxo1alC5cmUALly4wKxZs+jbty/Ozs5MnjyZsmXLMnToUC5dukSFChUYP348RkZG2NjY6CtsIYSOtWnzKWvXrsLLqzfdunmleowhDgFLjd4S7qlTp0hISMDHxwcHBweWLVvGsWPHNAnXwsKC4cOH069fP8zNzZk7dy73798H4NKlSzx48IABAwZQsWJF/Pz89BW2EEIHgoP/wcLCFXj1d/+33/Yme2rszekX8xO9lRTCwsIAsLa2xszMDEtLS8LDwzX7XV1dGTZsGObm5gQEBBAZGYmHhwdqtRpLS0tat27N4sWLCQsLY/LkyfoKO0+KjY1l0aLv6dChFY0be9ChQxvmzp1FZGRksmP8/ZfQteunmpUU5sz5lmfPIlO0pygKHTq0om7d6nz33Qx9fhRRAO3fv4dPPmlI7969NSt7vP2I7pvJNr+UE0CPPdz/pvpLfyq1rVu3Mn36dOrWrUurVq0wMjJi/fr/VqJt1aoVAQEBqNVqrZ+jTutxu/BwFSYmqbeR1va84LvvvuHAgX0UL16catXcuXnzJtu2bSEk5B+WLl1BYmIiI0cO4a+/LmNtbU3lylUJDv6boKBfuHz5IitXrqVw4f9WhTh37n+Eh7/6B/HIkQOMGTMOU1NTncWf0/dWpVJhZ5e9OSTSo8u2c1pejlVRFObMmcNXX32FoiioVCqsrQtjZmamOebWjh1E3bmDheWrbZWGDs2tcFPIiXurt4RrZ/dqwuvIyEjMzMyIiYnB3j75yq7bt29n4sSJeHh48MMPP2BkZERERASBgYHUr1+fihUrkpSUhImJSaYmrUhrLgW1Wp3qc/15eS6FFy+iOXhwP9bW1mzatA0zs8IkJsbTrVtnLl++yD//BHP8+O/89ddl3NwqMXfuQqysrIiKimLQIG9u3gwhKCiIjh27aNrcs2cX8GpS86ioKI4cOULjxmmv8JAduri3arVaZ3MIGNL8BHk51pcvXzJq1AgCAzcCMG7cBGbO/IbHj6OB/yabD7v6t+a1hWPpPPN5cmouBb0l3Nq1a2NsbIy/vz82NjYkJSVRq1YtTp48ibOzMy9fvmTixInY2NjQs2dPzp8/j62tLU5OTqxYsYIjR47g5eXFzp078fRMf+XY/MzY2AQjIyOioqLYtWsHzZu3okgRS+bNW0BERCR2dvbs3r0DAB+fwVhZvfoBKFKkCKNHjycqKopq1apr2ouPj+fw4YOoVCp69OjNqlXL2Lt3p84Srih4QkMf0qdPd86fP4eFhQWLF/vTunXbFL/tvjkELD+MuU2N3hKuo6MjM2fOxM/Pj7i4OMaMGUPJkiVp0qQJkyZN4uHDh8TFxREXF6cZ9uXp6cnChQtZuHAh06ZN46uvvsLDw4Px48frLM6H+/fw8sG9VHvEupKZWY8KFy5M/fqNOHLkIPPmzWLRovlUq1adRo2a0qxZcxITE7j3//UvV9f3k537ermdN506dYLo6CiqVKlGu3YdWL16OWfOnCIiIoJixYpl/8OJAu/HHxdz/vw5Spcuw5o1G3BzS31Jpvw2BCw1en3SrF27drRr1y7Zths3bmhejxkzJtXzatWqxc6dO3UamyH56qtJWFlZsW/fbuLj4zl9+iSnT59k27YtfPvtXM1xFhYWGba1f/+rGfPr129I8eJ2fPCBG1euXObgwb106pT1icyFeM3XdzKJiQmMHDmW4sWLp3rMm73b/DIELDXyaO9bSjZrnqdruABWVlZ89dUkhg0byalTxzl16gSHDu3n2rWrHD58UHPcixcvKFIk7UJ/dHQ0J08eB6BevYYANGjQiCtXLrNnzy5JuCJL1Go1/v5L6NGjF0WKFMXMzIwZM2ane05B6N2CJFyDc/nyRTZtCsDV9X169vSmWbPmtGjRkvfeK4u//xL+/fc2Dg4lCAsL5caNa1SvXlNz7s8/r+bs2T/o2bMP1avX5MiRg5rVcTt3bpvsOtevX+XOnds4Ob2nz48nDFxU1HOGDPFh797dnDp1kjVr1qd7/NvjbfNz7xZkLgWDY2FhyZEjh9i48WcePfpvHPPDhw+BV8uhN2/eCoDly3/ixYtXk9o9fvyIwMCNnDv3BwkJr5Ls/v//Nc7BoQRlyzpr/hQpUhSAPXukjCO0d/v2LVq1asbevbuxtramX78BGZ6TX8fbpkV6uAamXLnyNGjQiN9/P4yXV2c++MCNJ08ec/NmCEWKFKVlyzYULfoOJ08e46+/LtG1a3vKlavAjRvXiI6Ook6detSq5cHjx4+4cOEcxsbGLF26OlltbffuHcyY8TX79u1mwIAhBrUMtcgdx48fpV+/nkRERODi4soPg4djFfIPISH/pHnOQ8v/xt/m11EJb5MergGaOHEa3t4+2NoW5+LFP3n8+DF16tRj8WJ/7O0dKFy4MIsX+/PZZ16Ymppy+fKfWFsXw9vbh2nTZqJSqdi/fy9qtZpatT5O8UVG48ZNKVKkKGFhoVy4cC6XPqUwFKtWLadLl0+JiIigWTNPdu8+iFV8nNbnF4Se7WuyiKQsIqlXsoik7uRWrGPHjmT16hV83akrzarXQPXGb0Tp9Vzz273NUw8+CCHypxkzZuPp2YL3/r2TbHtB6rlqSxKuECJTrl69wtdfT8DffxWxZ88Qc+8u772xv6DUY7NCarhCCK3t3r2TVq2aceTIIdaN+zLFFIrSq02f9HCFEBlSFAU/v7nMnPkNAB06dKaZ+6s5OTLzaHpBV+ATrqIoMuzJgBWA73xzXUxMDCNHDuXXX7diZGTEhAlT6FTxA2Lvv+rdSrLVXoFOuMbGJiQkxGNqapbxwSJPSkiIx9i4QP8Y61RsbCzt2rXg4sULWFpa8dNPK/D0bEHIqmWAlBAyq0D/pFpZWRMZ+QhrazsKFTKVnq4BURSFhIR4IiMfUaSIzGqmK+bm5tSu/TGRkRGsW7eJihWTz0AnvdvMKdAJ19zcEoBnzx6TlJSo2a5SqTRLfxgCQ4o3J2M1NjahSJFimv+PIuc8f/6MokXfAWDKlOmMGjVOs5Juflu6XJ8KdMKFV0n37b+whjQgGwwrXkOKtSBKSkrim2+msGNHEHv2HKZ48eKYmJhoki0UnJm9dEGGhQkhgFe92h49OrNkyUIePLjP//73R7rHSzkh8wp8D1cIATdvBuPl1ZXg4H+wsbFh5cqf8fCom9th5TuScIUo4I4cOYSPTx+ePYvk/fc/ZO3aDanOg/z23LUi86SkIEQBdvNmMN26deTZs0iaN2/Fzp37tEq2Ur/NGunhClGAOTuXZ9iwL1CpjBg3biIqVep9sDe/KJPabdZJwhWigAkPD+fx40d88MGHwKtFHtMbg15QFnjUBykpCFGAXL58CU/PhnTv3omwsDCADB/4kWFgOUcSrhAFxG+/BdGmzSfcv3+PkiXfzfT50rvNPkm4QuRzarWa2bO/pV+/nsTExNC1a3d+/XUnDg4OuR1agSM1XCHysRcvXjB8+CB27AhCpVIxZcp0Bg0aKvOG5BJJuELkY2fOnGTHjiCKFCmKv/9KmjT5JLdDKtAk4QqRjzVu3IzvvptP3br1qVDBJdPny0Q1OUuvNdygoCAaNGiAh4cHy5cvT7F/1apV1K1bF3d3d0aPHk1sbCwAV69epV27dri7uzNu3Dhevnypz7CFMCgBAWs5e/aM5r23d/8sJVuQEQo5Lc0e7oMHD7Ru5N13M/7GMzQ0lAkTJtC+fXuKFSvGnDlzqFGjBpUrVwbgwoULzJo1i759++Ls7MzkyZMpW7YsQ4cOZfTo0RQtWhRfX1+mTJlC+fLl8fGRheqEeFNiYiITJoxl2bKfsLOz59Spc5opFjPr7SfLZIRCzkgz4TZu3DjVwvrrJU3e3Hft2rUML3Tq1CkSEhLw8fHBwcGBZcuWcezYMU3CtbCwYPjw4fTr1w9zc3Pmzp3L/fv3uX//PiEhIUyfPp2OHTsSGBjI0aNHJeEK8YaIiKd0796fAwcOUKhQISZMmJLlZAvIY7w6kmbC3bdvn+b1sWPHWLduHV999RVubm6YmJhw9epVvvvuO7p27arVhV4Psra2tsbMzAxLS0vCw8M1+11dXXF1dQUgICCAyMhIPDw8NOcVK/ZqPk4bGxtCQkIy+TGFyL/+/vsGPXt25datmxQvbseqVQHUqlU7R9qWJc9zVpoJt0yZMprXy5cvx8/PT9MbBfDw8OCbb75hyJAhdO/ePcMLpdYzTs3WrVuZPn06devWpVWrVpw/fz7DtjNia2uV6XPs7Ipk+7r6ZEjxGlKskLfj3bNnD127duX58+dUqVKFoKCgZH93s+LWjh1YWL5a50/Xnz0v39u35USsWo1SiI6OJjExMdXtCQkJWl3Izs4OgMjISMzMzIiJicHe3j7ZMdu3b2fixIl4eHjwww8/YGRkpDkmMjJS89+3z8vIkyfRqNXar+5qaKsSGFK8hhQr5P14nz6NIioqijZtPmXDhp+JiVFnO96wq38Dr0oJuvzsef3evkmbWFUqoww7d1ol3DZt2jB27FiGDx9OxYoVURSFy5cvs2jRIj777DOtAq5duzbGxsb4+/tjY2NDUlIStWrV4uTJkzg7O/Py5UsmTpyIjY0NPXv25Pz589ja2uLq6kqZMmXYsmULKpWKy5cvM2LECK2uKUR+pCiK5jfFZs2as337XmrUqImlpSUxMVlPYPJFme5plXB9fX2xsrJi9uzZPH36FIDixYvTvXt3Bg8erNWFHB0dmTlzJn5+fsTFxTFmzBhKlixJkyZNmDRpEg8fPiQuLo64uDgGDhwIgKenJwsXLmTBggX4+voyY8YMWrRoQa9evbL4cYUwbGFhofj49OGrryZqVmTIbr02tYnF5Ysy3TBSXhdXtfQ64drY2OgkIF2QkkLeYUixQt6K988/z9O7d3cePnxA5cpV2bfvSLLvRLIaa8iqZZrX+pzvNi/d24zotaQAcPnyZYKDg0lKSgJe/VoTHx/PtWvXmD59urbNCCGy4Ndft/D550N4+fIltWt7sGLFuhyZD+HNJ8lkRILuaZVwFy5cyJIlSyhevDhPnjzBwcGBx48fk5SUxCefyLPZQuiKWq1m1qzp+PnNBcDLqzezZs3D1NQ0y21KCSH3aPVo7+bNm5k6dSrHjx+nZMmSrFu3jpMnT+Lh4UGpUqV0HaMQBdYXXwzFz28uxsbGfPvtbObNW5itZAukmmzlCzL90CrhRkZGUq9ePQDef/99zp8/T9GiRRk5ciR79sjkFkLoSseOXShevDgbN/5C//6DslVGeLh/T7J6bTlvH8p5+0iy1SOtEm6JEiW4e/fVv4rlypXj6tWrAFhaWhIREaG76IQogMLCQjWvGzRoxB9/XKJBg0bZblce1819WtVwu3TpwsiRI5k5cyZNmzald+/e2Nracvr0aSpWrKjrGIUoMFatWs7kyeNZu3YjjRo1AcDKKvNPSqZWp31NvhzLPVol3AEDBlCiRAkKFy5MpUqV8PX1ZePGjVhbWzNjxgxdxyhEvpeQkICv71jWrFkBwLlzZzUJNyvSSrbSs81dWg8La9u2LfBqCrhOnTrRuXNnnQUlREHy5MkT+vfvxYkTxzA1NWXevIV07Zrx/CTakN5s3qL1BOQBAQE0a9aMKlWqcO/ePSZNmsT8+fNRq9W6jE+IfO3atat4ejbixIlj2Ns7sG3brmwnW1mlIe/Sqoe7atUq1q5dy+eff87XX38NQJ06dfjmm28A+PLLL3UWoBD5VWJiIr16fca//96mcuWqrFmznnff1W6YZWo12oeWZsS8iNO8l/JB3qNVD3fjxo1MmzaNTz/9VDMspXnz5syePZtt27bpNEAh8isTExMWLfqJLl26sX37Hq2TLaRdo31NxtbmTVr1cENDQ3nvvfdSbC9RogTPnz/P6ZiEyLdiY2M5evQInp4tAKhd24PatT2y3N6bNVpDmpugoNKqh1ulShW2bt2abJtarWb58uVUqlRJJ4EJkd88fPiAtm2b06vXZ+zduztrbbz18IIwLFr1cCdOnIiPjw9HjhwhPj6eKVOmcPv2beLj41NdfVcIkdy5c2fp3bs74eFhlCnzHmXKOGWpHXl4wbBplXArVKjA3r172b59Ozdv3iQpKYlmzZrRtm1bLC0tdR2jEAZt8+YNjBo1gri4OOrUqcfy5WuxtbXNVpsy3MswaT0O18zMTMbeCpEJSUlJTJ/+NT/8sACAPn36MWPGbAoVKpTpttJ7ckwYjkwvk56agwcP5lhAQuQXT548ITBwIyYmJsyYMRtv7/5ZbktKCflDmgn3zaVz7t+/z9q1a+natSsfffQRJiYmXLlyhY0bN9KzZ0+9BCqEobG3t2fNmvXExsZSt279LLXxds9WSgmGLc2E+2b5oFOnTpr1xF775JNP+Oijj1iwYAHDhg3TbZRCGIjffz/M5cuXGDbscwDc3Wtkqz3p2eYvWtVwQ0JCcHFxSbHd0dGRe/fu5XhQQhgaRVFYvvwnJk/2JSkpCXf36nz8cZ0styc92/xJq3G47u7ufPvttzx8+FCz7fbt20ybNk0zMbkQBVV8fDyjRo1gwoRxJCUl8cUXo6lV6+NstSk92/xJqx7ut99+y4gRI2jcuDFWVlYoisKLFy+oXbu2Zj4FIQqiR48e0bevF2fOnKJw4cL4+f1Ahw45N5pHerb5i1YJ197eno0bNxIcHExwcDDwamxuuXLldBqcEHnZjRvX6datI/fu3aVkyXdZs2Y9VapUy1JbMuyrYNB6HC5A+fLlKV++vK5iEcKgWFsXQ1EU3N2rs3r1ehwcSmS5rdSSrZQS8p9MJVwhCjpFUVAUBZVKhYODA1u3/sa775aicOHCOdK+lBDyN0m4QmjpxYsXfP75EJyc3mPSpKkAODtnvawmZYSCR+sVH4QoyO7du0ubNp5s3/4rq1evSLayblaklmylhJD/pdnDPXXqlNaNfPyxdkNggoKCmD9/PgkJCfTt25f+/VM+6nj79m3atGnDuHHj8PLyAiAwMJCJEydqjmnYsCFLly7VOj4hsuPMmdN4e/fg8eNHlC3rzLp1m7JVr4X/arYyUXjBkmbC9fb21qoBIyMjrl27luFxoaGhTJgwgfbt21OsWDHmzJlDjRo1qFy5suaYU6dOMX78eOLj45Ode+nSJSpUqMD48eMxMjLCxsZGq9iEyK6VK1cyaNAgEhISqF+/EcuXr8bauli22nxzzTFJtgVLmgn3+vXrOXqhU6dOkZCQgI+PDw4ODixbtoxjx45pEm5UVBTe3t60bduWoKCgZOdeunSJBw8eMGDAACpWrIifn1+OxiZEalavXsHYsSMB8PEZxNSp32JikvmvPdKq1UoJoeDRuoabkJDA/fv3uXv3Lnfv3uXff/8lODiY7du3a3V+WFgYANbW1piZmWFpaUl4eLhmv5mZGbt27WLEiBHJzlOr1VhaWtK6dWsWL15MWFgYkydP1jZsIbKsTZtPqVixIvPnL2LGjNlZSraQ9pAv6d0WPFr9BO3bt49Jkyalun6Zg4MDbdu2zbANRVEA0pzy0dTUFGdn5xRzM6hUKtavX69536pVKwICAlCr1ahU2v17YWtrpdVxb7KzK5Lpc3KTIcWbl2O9desWpUuXxsTEBDu7Ily8eBFTU9Ost7djBxaWZgBUGjo0p8JMU16+t6kxpHhzIlatEu78+fPx9PSkT58+dOvWDX9/fyIjI/nmm2+STeOYHjs7OwAiIyMxMzMjJiYGe3v7DM+LiIggMDCQ+vXrU7FiRZKSkjAxMdE62QI8eRKNWq1ofbyhLcZnSPHm5VgPHdrPgAF9+eyz7kyf/h2gfbwZDfGycCyt88+dl+9tagwpXm1iVamMMuzcaZVw7927h7+/P2XKlMHNzY1Hjx7RtGlTjI2NmSKz53cAACAASURBVDVrFh07dsywjdq1a2NsbIy/vz82NjYkJSVRq1YtTp48ibOzMyVKpP6tr7m5OStWrODIkSN4eXmxc+dOPD09tQlbCK0oisJPP/3A1KkTUavV3L9/n8TExEyVEDJKtlI+EKBlwi1atCixsbEAlC1blmvXrtG0aVPKli2r9fSMjo6OzJw5Ez8/P+Li4hgzZgwlS5akSZMmTJo0STME7G2FCxdm4cKFTJs2ja+++goPDw/Gjx+v5ccTIn0vX75kzJgv2LTpVdlq9OivGD36q0z9BvUmeVJMpEerhNuoUSOmTp3K1KlTqVWrFt999x0NGjRg7969ODg4aH2xdu3a0a5du2Tbbty4key9o6Njim21atVi586dWl9HCG2EhYXRp093zp07i4WFBYsW/USbNp9mup03h3kJkR6t/hn39fXF2dmZK1eu0LRpU6pUqULXrl3ZsGED48aN03WMQujEjBlfc+7cWUqVcuS33/ZlKdlC8ocYhEiPVj1cS0tLpk+frnk/d+5cRo0aRfHixbO0AqkQecE338xErVYzefI3Wn2Bmxp5iEFkhlYJ99GjR3z11Ve4ubkxcuSrgeCdOnXCzc2NmTNnypNfwiCo1WoCAtbStWt3TE1NeecdaxYvztwj4vIQg8gOrUoKU6ZMAV4l2dcCAgJITExM1vMVIq+Kjo7G29tLsxROVslDDCI7tOrhnjlzhi1btlC69H//ir/33nv4+vry2Wef6Sw4IXLCnTu36dWrG9euXeGdd6xp2bJ1ttuU0QgiK7Tq4VpaWqY6/CssLExquCJPO3nyOJ6eDbl27Qrly1dgz56DNGrUJEttyWgEkV1a9XA7duzIhAkT+Pzzz/nwww8BuHr1KosWLaJ9+/Y6DVCIrFqzZiXjx48mMTGRxo2bsnTpSt55xzrL7cloBJFdWiXc4cOHoygK8+bN4+nTpwDY2NjQq1cvfHzkVyuR9yiKwsGD+0lMTGTw4OFMnjwNY2PjLLcnoxFETtAq4apUKr744gu++OILnj59iqmpKVZWmZ8QRgh9MTIyYskSfw4fPpjl8bVvkt6tyAlpJtwtW7bQtm1bTE1N2bJlS7qNvDl6QYjccv36Nfz85uDnt4TChQtjZVUkR5Kt9G5FTkkz4S5ZsoQmTZpgamrKkiVL0mzAyMhIEq7Idfv27WbQoP5ER0dRrlwFxozJ/nwbD/fv4WFEODEv4gDp3YrsSzPhHjp0SPN66tSpVKtWDUtLS70EJYS2FEVh0SI/Zsz4GkVRaN++I0OHfp4jbcfcu6uZy1bG2oqcoNWwsNGjR/PgwQNdxyJEpsTGxjJ4cH+mT5+CoiiMHz+Jn35aiYWFRbbbfrOMUM7bR5KtyBFafWnm6urK2bNnqVChgq7jEUIrz58/o3Pndly4cB4LC0t+/HE5LVq0yrH2X39JVsTJKcfaFELryWu++eYbFi1aRKlSpTAzM0u2PyAgQCfBCZGWIkWKUrq0E0+ePGHt2o188MGHOrlO2datDWZVApH3aZVwP/zwQ80DD0LkppcvX1K4cGGMjIxYuPBHYmJiKF68eI61n9FSOUJkh1YJd9iwYbqOQ4h0JSUlMWPGVI4ePcL27XuwsLDQ/MlJbyZbGZUgclqaCXfs2LFMnjwZKysrxo4dm24js2fPzvHAhHgtKuo5gwb1Y//+vRgbG3P69EkaN26q02vK5DRCF9JMuG8+BpmdRyKFyI6bN0Po1esz/v77BsWKFWP58rXUq9dAJ9eSyWmErqWZcGfOnJnqayH05dix3+nXryeRkZG4ulZk7dqNlC3rrJNrvVm7lVKC0BWt14E+evQomzdvJiQkBJVKhaurKz169MDd3V2X8YkC6vLli3Tp8ilJSUl4erZgyZJlFClSVGfXezPZyphboStaJdxNmzbxzTff0Lp1a3r06EFSUhJ//fUXffr0Ye7cuXh6euo6TlHAuLlVomPHLpQs+S5ffTUxx8taaY1GkGQrdEmrhPvTTz8xbdo0OnTokGx7jRo1mD9/viRckSMeP35MbGwMpUuX0Qz7Uqm0ehgyU2RdMpFbtEq4z58/p1KlSim2V69eXeq7IkdcufIXvXp9hpWVFTt37sfKqohOki1I+UDkHq1+or28vPjuu+80k48DvHjxgh9++IFevXrpLDhRMOzc+RutWjXj7t1/MTc3JyYmVi/XlWQr9E2rHu7p06e5cuUKDRs2pHTp0hgbG3P37l1iY2NxcHAgKChIc+yRI0d0FavIZxRFYf782Xz33QwAOnXqyvz5iyhcuHCOX0ueIBN5gVYJV1bmFTktJiaGESMGs337rxgZGTFx4lSGDfscIyMj3VzvrWQr9VqRG7RKuDm1UGRQUBDz588nISGBvn370r9//xTH3L59mzZt2jBu3Di8vLyAVwtWjh8/nnv37tG0aVOmTp2qk16Q0J9du35j+/ZfsbIqwtKlK2iWQ7/eZ9STlSfIRG7SKuHevHkTPz8/bt26RXx8fIr9e/fuzbCN0NBQJkyYQPv27SlWrBhz5syhRo0aVK5cWXPMqVOnGD9+fIprjB49mqJFi+Lr68uUKVMoX768LF5p4Dp27MKtWzdp27Y9rq4Vc6zd9JKt9GpFbtMq4Y4aNQqVSkXHjh2z3LM8deoUCQkJ+Pj44ODgwLJlyzh27Jgm4UZFReHt7U3btm2T1YTv379PSEgI06dPp2PHjgQGBnL06FFJuAZo06b1fPJJI4oVK4mRkVGOLIOTFunJirxIq4R769YttmzZQvny5bN8obCwMACsra0xMzPD0tKS8PBwzX4zMzN27dqFqalpsoT7+rxixYoBr5ZnDwkJydS1bW0zv8KwnV2RTJ+Tm/JyvImJiYwdO5bvv/8eV1dXLly4gLm5uU6u9fD/l8TJyfuRl+/t2wwpVjCseHMiVq0Sbr169bhw4UK2Eq6iKABpfiliamqKs7Mz9+7dS/W87HjyJBq1Wvt27OyKGNSk03k53mfPIhkwwJvDhw9SqFAhRo0aRXR0ItHRORdvanXbnLofefnevs2QYgXDilebWFUqoww7d1ol3PHjx9O+fXu2b99OqVKlUiRNbR5+sLOzAyAyMhIzMzNiYmKwt7fP8LzXx0RGRmr+q815IvcFB/9Dz55dCQkJxtbWllWrAmjTxjPH/5LJCARhKLRKuJMnTwZe/VqflJSUpQvVrl0bY2Nj/P39sbGxISkpiVq1anHy5EmcnZ0pUaJEqueVLl2aMmXKsGXLFlQqFZcvX2bEiBFZikHoz6FDBxgwwJvnz5/xwQdurF27gTJlcnZ9sLd7tlK3FXmdVgn37NmzrF+/PlvL7Dg6OjJz5kz8/PyIi4tjzJgxlCxZkiZNmjBp0iTNELDULFiwAF9fX2bMmEGLFi3k6TYD8ODBfZ4/f0bLlm1YvHgpVlaZr6OnJbUSgvRqhSHQKuFWqFCB58+fZ/ti7dq1o127dsm23bhxI9l7R0fHFNs++OADtm3blu3rC/3x8upNiRIlaNy4WY7PifD2MjjyiK4wFFol3C5dujBmzBjat2+Po6NjiqnyOnXqpJPghOEICwtjxIhBTJ36LRUrvg9A06Y5O4uclBCEodN6ekZTU1N27tyZYp+RkZEk3ALu0qU/6dWrGw8e3Cc+fjS//pry5yQ7pIQg8gutEu6hQ4d0HYcwUEFBvzBixGBiY2OpUaMWS5euyvFrSAlB5BdpJtyDBw9Sv359ChUqlObJr6dozGhVX5H/qNVqZs/+lvnzX63Y3K2bF7Nnf4+ZmVm2205rPgQpIQhDl2bCHTZsGMePH8fW1lazrWHDhgQEBFCqVCkAYmNjWbVqlSTcAmjIkP788suroXpffz2dgQOHpvpQy9vJ86GlGTEv4jJ9PSkhiPwgzYSb2hNez549Q61W6zQgYRiaNPmEgwcPsHTpSho3bprmcVmdg1ZKByI/0nrVXiGePYvknXesAejc+TOaNGmGjY1tBme98rocYEiPcwqR03SzaJTId9atW427+0dcvnxRsy2tZPtw/x5CVi0jZNUyfYUnhEGQhCvSlZCQwPjxoxk1agTPnz/j0KEDGZ4jQ7iESF26JYUdO3ZgaWmpea9Wq9m9ezc2NjYAREdH6zY6kasiIp7Sv39vjh37HVNTU+bM8aNbt5SPYMuoAiG0k2bCfffdd1mzZk2ybba2tmzcuDHZtpIlS+omMpGrbty4Ts+eXbl9+xZ2dvasWhVAzZq1NPszWspGerVCpJRmwpWHHQqu2NhYOnZsQ3h4GJUqVWHNmvWUKuWY7JjUygYyqkCI9MkoBZGCubk5M2fO4bfftuHntwQLCwvNPpnPQIisk4QrgFe92osXL1C7tgcAbdp8Sps2n2r2y3wGQmSfJFxBaOhDevfuxtWrV9i2bRfu7jWAtOu0Uj4QImsk4RZwFy6co3fv7oSGPqR06TIULvzf4o5SpxUiZ0nCLcC2bt3MF18MJS4ujo8/rsOKFesoXrw48Kp3+5rUaYXIGfLgQwGkVquZPv1rBg/uT1xcHD17ehMYGKRJtvBf71bqtELkHOnhFkD//nuHFSv8MTY2Zvr07+jb1yfN5eulhCBEzpGEWwC9915Z/P1XYmZWmPr1G6bY/2Y5QQiRcyThFhDHjv3Ow4cP6NKlGwDN0um5SjlBCN2QhJvPKYrCypXLmDhxHCqVig8+cMPN7SOtzpVyghA5SxJuPhYfH4+v71jWrl0JwJAhI3j//Q/SPUfKCULojiTcfOrJkyf069eTkyePY2ZmxvffL6ZTp65pHv/2Qw5SThAi50nCzYeuX7+Gl1cX/v33Dg4OJVizZj3VqlVP9xxZGVcI3ZOEmw8VKmRCZGQkVatWY/Xq9ZQs+W6ax8pkNELojyTcfOL1op9GRkaUK1eBX3/dQfnyLpibm6c4Nr05EoQQuqPXJ82CgoJo0KABHh4eLF++PMX+EydO4OnpSc2aNZk1a5YmiQQGBuLq6qr5M3DgQH2GnefFxMQwaFBfVq7012z76KPKqSZbSH2OhHLePlJGEELH9NbDDQ0NZcKECbRv355ixYoxZ84catSoQeXKlQGIi4tj5MiRVK9enZo1azJz5kwqVapEy5YtuXTpEhUqVGD8+PEYGRlplvgR8ODBffr18+LcuXMcPnyQTp26albWfVNqvVopHwihX3rr4Z46dYqEhAR8fHwYOnQoKpWKY8eOafZfvHiRZ8+e0bNnT/r06YO9vb1m/6VLl3jw4AEDBgxg3rx5ydZZK8jOnj1Ds2YNOHfuHE5O77F9+16tk62UD4TQP731cMPCwgCwtrbGzMwMS0tLwsPDU90PYGNjQ3h4OGq1GktLS1q3bk2jRo2YNGkSkydPZtWqVVpf29bWKtPx2tkVyfQ5+rRmzRoGDBhAfHw8jRo1IjAwEFvbNJYtjwjHwtKMIk5OlG3dWs+RppTX7+3bDCleQ4oVDCvenIhVbwn3zS91MrNfpVKxfv16zftWrVoREBCAWq1GpdKug/7kSTRqtaJ1rHZ2RXj0KErr4/Vt6dIfmDRpPAB9+/rw008/EBn5MtWYH+7fQ8yLOABK1mqQ658rr9/btxlSvIYUKxhWvNrEqlIZZdi501tJwc7ODoDIyEji4+OJiYnB3t4+xf6IiAjNcfb29kRERODv78/169cBSEpKwsTEROtkmx81b94Ke3sH5szxY9aseRQqVCjNY2VeBCHyDr31cGvXro2xsTH+/v7Y2NiQlJRErVq1OHnyJM7OzlSuXBkrKyvWrl3LP//8Q2hoKPXq1cPc3JwVK1Zw5MgRvLy82LlzJ56envoKO88IDX2Ig0MJjIyMcHJ6jzNn/sxULVtGIAiR+4y//vrrr/VxoaJFi+Lo6MjmzZu5dOkSw4YNo2rVqnTo0IHSpUvj7u6Om5sb27Zt4+jRo3Tr1g1vb28KFSpEpUqV2LdvH7/88gs1atRg6tSpFC5cWOtrx8bGo2hfUcDS0oyYmPgsfErdOHz4IB06tMbY2ISaNWsBYGpqysP9ewg/epinf54n/I8/iPjzfIo/r9lUdc+t8JPJa/c2I4YUryHFCoYVrzaxGhkZYWFhmv4xipKZVGSYDLWGqygK633HcPy3bSiKgksFF9q268DbVXALSzNNnTY1eelR3bxyb7VlSPEaUqxgWPHmVA1XnjTLo+Li4hg37kte/v/sXR/X/pg6deqnTLaOpanUvbPB/OAKUZBJws2DwsPD6dvXiz/+OE3XkiVp3rwVrb6bn9thCSGySRJuHvTD4H443grhfWdnPv20IyUcSuR2SEKIHCAJNxelNYnMpw0acigxgabNmmNlaSlDuoTIJyTh5qLXyVZRFP7++zouru9jBBQtUpTuw77IM190CSFyhiRcHUurF/tafEICc06fYseOIMaNm8CoUeP0GJ0QQp8k4epQRsn22fNnrPx1KzuuXqFo0XeoWrWaHqMTQuibJFwdevOx2rfLA6dPn6Jv3x48fvwYZ+dy/PzzZsqXr5AbYQoh9KTgTkigR28n24CAtXTs2JrHjx/ToEEj9uw5JMlWiAJAEq6OpLXceGJiIj//vIaEhAQGDhzKhg1bsbYupufohBC5QUoKOpLWLF0mJiasXr2eo0cP07nzZ7kRmhAil0gPVwfe7N2WbNacv/++ga/vGNRqNQAODg6SbIUogKSHmwUZjT54zcKxNAcO7GXgwH5ERT2nbFlnfHwG6yFCIUReJD3cTNI22ZqXKs0vf/9Njx5diIp6Tps2n9K9ey89RCiEyKukh5tJ6Q31eu3ly5eMHv05mzdvAGDsWF++/HJsgV6lQgghCTdT3q7Npubp0yf06NGFc+fOYmFhwaJFS2nTpp2+QhRC5GGScNOQXukgvclkrKyKYGpqiqNjadau3Yib20e6ClEIYWAk4aYhvWSbWu82MTERExMTTE1NWbnyZ9RqtWZhTCGEAEm4GSrn7ZPufrVazaxZ07l06U9+/nkzJiYm2Nra6ik6IYQhkYSbils7dmh1XHR0FEOG+LBnzy6MjY353//+oHZtDx1HJ4QwVPK1eSqi7twB0q/V3r59i5Ytm7Jnzy6sra3ZuPEXSbZCiHRJDzcdaY1EOHHiGP369eTp06dUqODCunUbcXYur+fohBCGRhLuW9KadOa106dP0blzOxITE2na9BN++mkFRYu+o6fohBCGTBLuW2Lu3cXC0izNckL16jWoU6cebm6VmDjxa4yNjfUcoRDCUEnCTcOb5YSnT58AYGNji4mJCevXb6FQoUK5FZoQwkDJl2YZuHbtKp980oh+/XqRkJAAIMlWCJElknDTsWfPLlq2bMq//97mxYtonj9/ntshCSEMmF4TblBQEA0aNMDDw4Ply5en2H/ixAk8PT2pWbMms2bNQlEUAK5evUq7du1wd3dn3LhxvHz5UreBKgoLFsyjd+9uvHgRTYcOnQgK2iMPNAghskVvCTc0NJQJEyZQv359OnXqxJw5c7h48aJmf1xcHCNHjqRcuXIMGTKEVatWsXv3bgBGjx6Nubk5vr6+7Ny5k3Xr1ukszoTERH755RdmzJgKwIQJU/jxxxWYm5vr7JpCiIJBbwn31KlTJCQk4OPjw9ChQ1GpVBw7dkyz/+LFizx79oyePXvSp08f7O3tOXbsGPfv3yckJISOHTvSsWNH3NzcOHr0qM7ivHzpTy7/9ReWllasWbOBzz8fhZGRkc6uJ4QoOPQ2SiEsLAwAa2trzMzMsLS0JDw8PNX9ADY2NoSHh2u2FytWTLM9JCQkU9e2tbXS+ljP7p8RnhjP1K+/xs3NLVPXyU12dkVyOwStGVKsYFjxGlKsYFjx5kSseku4r+uxafUW09r/ent2PHkSjVqtXTtFP27E123b8uhRFI8eRWX72vpgZ1dEYtURQ4rXkGIFw4pXm1hVKqMMO3d6Kym8nqowMjKS+Ph4YmJisLe3T7E/IiJCc5y9vb3mmMjIyGTbhRDC0Oith1u7dm2MjY3x9/fHxsaGpKQkatWqxcmTJ3F2dqZy5cpYWVmxdu1a/vnnH0JDQ6lXrx6lS5emTJkybNmyBZVKxeXLlxkxYoS+whZCiByjtx6uo6MjM2fO5Pjx4wQGBjJmzBhKliyJt7c3Bw4cwMLCgoULF3Lz5k0WL15Mnz59aNGiBQALFizg5cuXzJgxgxYtWtCrlyzGKIQwPEZKThRJ87jM1HDBsGpLYFjxGlKsYFjxGlKsYFjxGlwNVwghCjpJuEIIoSeScIUQQk8KxPSMKlXmnxTLyjm5yZDiNaRYwbDiNaRYwbDizShWbT5LgfjSTAgh8gIpKQghhJ5IwhVCCD2RhCuEEHoiCVcIIfREEq4QQuiJJFwhhNATSbhCCKEnknCFEEJPJOEKIYSeFLiEazBLtWcj1sDAQFxdXTV/Bg4cqPNYtYkX4Pbt23z00Uf8/PPPmm25cW+zE29u3N+MYl21ahV169bF3d2d0aNHExsbC+Tde5tWvHnx3q5cuZKGDRtSq1YtZs+enb2coBQgDx8+VD788ENl4sSJyrx58xQXFxflzz//1Ox/+fKlUqNGDWXw4MHKqlWrFBcXF2Xnzp2KoihKixYtlK5duypbtmxRPvzwQ8Xf3z/Pxjpx4kSlVatWyvHjx5UTJ04o165d02ms2sSrKIpy8uRJpUGDBoqLi4uybt06zXZ939vsxqvv+5tRrOfPn1dcXFyUWbNmKZs3b1YqVqyoLF68WFGUvHlv04s3r93bixcvKi4uLsr8+fOV7777TnFxcVEOHjyoKErW7m2B6uEaylLt2YkV4NKlSzx48IABAwYwb948LC0tdRqrNvFGRUXh7e1NzZo1k52XG/c2O/GC/u9vRrFaWFgwfPhwRowYQefOnSlatCj379/Ps/c2rXgh793bDz/8kH379jFkyBDKlSsHQKFChbJ8bwtUws3JpdrfPC8vxapWq7G0tKR169YsXryYsLAwJk+erNNYtYnXzMyMXbt2pViPLjfubXbizY37m1Gsrq6uDBs2DHNzcwICAoiMjMTDwyPP3tu04s2L99bY2BgnJyfWr1+Pr68vjRo1ol69elm+twViesbXlFxcqj2zshqrSqVi/fr1mvetWrUiICAAtVqNSqW7f18zitfU1BRnZ2fu3buX6nn6ltV4c+P+ZhTra1u3bmX69OnUrVuXVq1acf78eZ3Ek5GsxmtkZJRn723Dhg0pXrw448ePZ8WKFVSpUiVL1ytQPVxDWqo9q7FGRETg7+/P9evXAUhKSsLExESnyVabeNOSG/cWsh5vbtxfbWLdvn07EydOxMPDgx9++AEjI6M8fW9Tizcv3tsHDx6wY8cOnJycaNOmDWXKlOGPP/7I8r0tUD1cQ1qqPauxmpubs2LFCo4cOYKXlxc7d+7E09NTp7FqE2+JEiVSPS837m124s2N+5tRrC9fvmTixInY2NjQs2dPzp8/j62tLa6urnny3qYVr5OTU567t3fv3mXUqFFcv36dkiVLcuvWLT799NOs/9zm7Hd+ed+2bduUhg0bKh9//LGybNky5e7du8m+hT5+/LjyySefKDVq1FC+/fZbRa1WK4qiKFeuXFHatWunVKtWTRkzZozy8uXLPBvr6dOnlZYtWyofffSRMnDgQCUiIkLnsWoTr6IoqW7LjXubnXhz4/6mF+vs2bMVFxeXZH+GDx+uKErevLfpxZvX7q2iKMqKFSuUOnXqKDVr1lSmTZumxMfHK4qStXsrKz4IIYSeFKgarhBC5CZJuEIIoSeScIUQQk8k4QohhJ5IwhVCCD2RhCs0Xs/QdPfu3RT7NmzYgKurK99//32W2r5z5w6urq4pntxKzZkzZ3B1dSUxMTHdOF//cXNzo0WLFmzZsiVLsaWlfv36/PLLLwD07NlTq88eHx/Pxo0bs3zNwMBAGjdurPXx3bp1Y9GiRVm+Xnq+//57evbsqZO2C6oC9eCDyFihQoU4dOgQvXv3Trb9wIEDGT7+qE9+fn5Ur14deJXkdu3axYQJE3BycqJGjRo5fr1FixZRqFChDI/buXMnS5Ys4bPPPsvxGIThkx6uSKZ69eocOnQo2bbo6GguXLjABx98kEtRpfTOO+9gZ2eHnZ0dpUqVwsfHh7Jly7Jv3z6dXM/a2lqrmatkWLtIjyRckUyTJk04d+4cUVFRmm2///471atXT5FwDh8+TPv27alUqRItWrRg9+7dmn0JCQlMmzaN6tWr06BBg2RT3sGr6Q/HjRuHu7s7derUYdKkSURHR2crdhMTE00vtHHjxsyePZu6devSsmVLEhMT+eeff+jVqxeVKlWiWbNmrFy5MlmC3LhxIw0aNMDd3Z2lS5cma/vtksK6deto0qQJVapUoVevXoSEhHDmzBnGjx9PWFiYpnyiKApLliyhXr16uLu7069fP27fvq1pJywsjP79+1OlShU6dOiQYcll//79eHp6UqVKFaZPn54iwW/atIkmTZpQtWpVunXrxqVLlzT7wsPDGTFiBDVq1MDNzY1PP/2Us2fPavYHBwfTrVs3KleujLe3t2aeAJFzJOGKZMqVK0epUqWSze158OBBmjZtmuy4U6dOMXz4cNq1a0dQUBBdu3Zl9OjRmr/gixYt4siRI/z444/4+fmxbt26ZOf7+voSERFBQEAAS5cu5datW4wfPz5LMcfFxREQEEBwcHCy+uf27dtZvnw58+bNIzExUZPYXk+csmbNGs1KDseOHWPGjBmMHDmSjRs38ueff2qm4HtbYGAg33//PSNHjiQoKIgSJUowZMgQqlatiq+vL3Z2dhw/fpySJUvy888/ExQUxOzZs9m8eTNOTk707t1bs8LB559/jlqtJjAwkP79+7N27do0P2dwcDBffPEF3bp1Y+vWrcTHx3PhwgXN/kOHDrFgwQLGjx/Pr7/+Sv369endu7dm2sCxY8eSmJjIxo0b2bZtGyVKlGDKlCnAq7LMgAEDcHR05JdffqFp06YEBgZm6f+H4prMxgAABt5JREFUSIeOHk8WBsjFxUU5ceKEMmvWLOXLL79UFEVR4uPjlRo1aiiPHz9WvLy8lPnz5yuKoihDhw5VPv/882Tnf/HFF8rw4cMVtVqt1K5dW9m6datm36FDhxQXFxfl7t27yp07dxRXV9dkz8nfu3dPcXFxUR48eKCcPn1acXFxURISEtKM86OPPlKqVKmiVKlSRXF1dVVq1qyprFq1SnNMo0aNlFmzZmneb968WWnTpk2ydn755RelSZMmiqIoyvDhw5Vx48Zp9j158kRxc3PTfIY3P3uHDh00rxVFUZ4/f67MmjVLef78ubJ161alXr16mn3169dX9u3bp3mvVquVxo0bK7/++qvy999/a+7Ja7NmzVIaNWqU6ueeNWuW4uXlpXkfHx+v1KlTR1m4cKGiKIrSrVu3ZPfgddyvV1NYvXq18uDBA82+o0ePKi4uLoparVYOHz6sVKlSRYmOjtbsHz58eLLrieyTL81ECk2aNGHIkCEkJiZy+vRpypcvj62tbbJjQkJC6NKlS7JtVatWZfPmzURERPD06VMqVqyo2efm5pbsXEVRaNSoUYpr3759W6vp+KZNm0bVqlUxMjLCzMwMe3v7FF/qlSpVSvP65s2bBAcHU7VqVc02tVpNfHw88fHxhISE0LlzZ80+GxubZOe//dkHDRqkeV+kSBHGjRuX4rgXL14QGhrK6NGjk32muLg4bt++jZmZGVZWVjg6Omr2ubm5sXfv3jSv6+rqqnlfqFChZO9DQkKYP38+CxYs0GyLj4/XzHzWrVs3du3axfnz57l16xZ//fUX8GoaxODgYEqXLp2sbOTm5paiFCSyRxKuSKFatWoYGxtz7tw5Dh48SLNmzVIcU7hw4RTb1Go1SUlJmvfKG/VFE5P/ftSSkpKwsLBg27ZtKdqws7Pj8uXLGcZob2+Pk5NTuseYmZlpXicmJlKzZk2mTp2a4rjXsSlv1UPTGpWgzWgFQHMv5s+fT/ny5ZPtK1KkCGfOnElxzTfvkzbejCUpKYlx48ZRt27dZMdYWFigVqvp27cvz549o2XLljRu3JiEhASGDRumOS67sYiMSQ1XpKBSqWjYsCGHDh3i8OHDKeq3AM7Ozly8eDHZtgsXLlC2bFmKFStG8eLFkyXOa9euaV6XLVuWmJgYkpKScHJy0iTOmTNnZvuLs7SULVuW27dvU6pUKc01r127xrJly1CpVFSoUCFZvNHR0amORwZwcnLi6tWrmvcxMTHUqVOHv//+O1kvu2jRotja2vLo0SPNNR0dHZk/fz43btzAxcWFFy9ecPPmTc05b7b7trdjTEpK4saNG8k+Y2hoqOZaTk5OrFy5kj/++IPg4GDOnj3LihUrGDx4MA0bNtTUdhVFoUKFCvz77788e/ZMq1hE1kjCFalq0qQJgYGBWFtbU7p06RT7+/Tpw/79+1m9ejW3b99m9erV7N+/nx49emBkZET37t1ZvHgxJ06c4NKlS8yaNUtzbrly5ahXrx5jx47l4sWLXL9+nXHjxvHkyROdrUjQtm1b4uPjmThxIiEhIZw4cYJp06bxzjvvANCjRw/27dvHxo0bCQkJYeLEicTFxaXaVq9evVi3bh179+7l9u3bTJkyBWtra8qXL4+FhQVRUVHcunWLxMRE+vTpw4IFCzhw4AB37txh6tSpmsmty5UrR+3atfH19eX69escOHCADRs2pPkZOnfuzNWrV1m8eDE3b95k5syZhIaGavZ7e3uzbt06fv31V/79918WL17M1q1bcXZ2pmjRoqhUKnbt2sX9+/fZs2eP5oGJ+Ph4PDw8ePfdd/H19SU4OJgtW7akWdoQWScJV6SqTp06JCUlpdq7Bfjoo4+YO3cumzZtonXr1mzduhU/Pz/q1KkDwODBg2nfvj1ffvklgwYNomvXrsnOnz17Nk5OTvTt2xcvLy/s7e1ZsmSJzj6PlZUVy5f/Xzt3jKIwEIZh+EXIDYRYBouABxBscoEEBidYBFuL1GJhY5EUwQsoaqLYewLBgwTsbU3jAbZY2Gphm2Wq72kHpvnhZZp/zrxeL6y1rNdrrLUsl0sAxuMx2+2WpmmYzWb4vk8Yhr/eZYwhz3OqqsJay/v95nA40Ov1mEwmDIdDjDG0bctisSDLMsqyxBjD8/nkcrng+z7wvcDR7/fJsuzPza4gCDgej9zvd6bTKV3XEUXRz3kcx6xWK3a7HUmS8Hg82O/3jEYjBoMBRVFwvV5JkoTT6cRms8HzPNq2xfM86rrm8/mQpim32435fP6PExAAfUAuIuKIXrgiIo4ouCIijii4IiKOKLgiIo4ouCIijii4IiKOKLgiIo4ouCIijii4IiKOfAGt4AzAntQbCQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(5,4))\n",
    "plt.plot(np.array(empirical_dead),np.array(empirical_dead), '--', color='black', lw=2, label='Empirical')\n",
    "plt.step(x=np.array(km_dead_dist_sca), y=np.array(empirical_dead), label='SCA', color='brown', lw=2, alpha=0.5)\n",
    "plt.xlabel('Model Predicted dead',  fontsize=fontsize)\n",
    "plt.ylabel('Empirical dead',  fontsize=fontsize)\n",
    "plt.tight_layout()\n",
    "plt.legend(loc='best', fontsize=legend_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAAENCAYAAABU5k23AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOydd3zTdf7Hn9lJR5quFCizwFGgUvYeUvDQgj/kUHBwKii4UM/BISggCIcn6jlQTpBTz9PzhDvkTlDZtAyLILMISkEEpHvPNOP3R0lI27RN0iRN4PN8PPqAfOc77bevvvP+vIfEYrFYEAgEAoHHkLa0AQKBQHCtIYRVIBAIPIwQVoFAIPAwQlgFAoHAwwhhFQgEAg8jhFUgEAg8jBBWgUAg8DDyljbAFxQUlGE2O5+uGxkZQl5eqRct8g6BaHcg2gzCbl/jb3ZLpRLCw4Mb3H9dCKvZbHFJWK3nBCKBaHcg2gzCbl8TSHaLUIBAIBB4GCGsAoFA4GGEsAoEAoGH8bmwHjhwgG7durFz5856+/bu3cu4ceMYOHAgL7/8Mtb+MCdPnmTixIn069ePuXPnUllZ6WuzBQKBwGl8tnhlsVj45ptvWLhwocP9VVVVPPXUU/Tv35+BAweyfPlyevXqRXJyMs8++yxarZb58+ezaNEiunTpwsyZM31lukDgFhaLhYKCHAyGSsA/Fl6ys6WYzeaWNsNlfG23TCYnJESHRtPwyn9j+ExYT58+zbPPPsu4ceP48ssv6+0/evQoRUVF/P73v2fIkCGsXbuW1NRUEhMTycjIYOnSpUyePJl169aRkpIihFXg95SWFiGRSIiJaYtE4h9RN7lcitEYeMLqS7stFgvV1QYKC3MA3BJXnwlr69at2bZtG+fPn3corFlZWQDodDoAIiIiyM7Otm0PDw+3bc/IyPCanRP69sRQWMCc9u0bPCZy4GB6vfgnpHI5EoUCiUTiNXsEgUtFRSkRETF+I6oC55BIJCiVKnS6aIqKcv1bWMPCwggLC+P8+fMO91vjqXVFyhN9uCMjQ5w+9sDFC00ek3fgW3YmJwEQ1rUrvR57jOC2bdFERaEMC0OmVLpta3OJjg5tsXu7SyDaDE3bnZ1tQaVS+t0fXrk8MIXe13bLZGoKC81uPZ9+UyAQHR0NQEFBAQCFhYX06NEDvV5ve23917rNWfLySl1OLv6gTVs+/XR9ve2Hn3ua3G/32V4X/fQTBXnF5F74DsuVeyjDdGjatiU4th0yjcal+zaH6OhQcnJKfHY/TxCINoNzdpvNZkwmC/4SXwURCnAVs9ns8OcslUoaddhaVFirqqo4dOgQcXFxJCYmEhISwt///nd++uknMjMzGTFiBO3ataN9+/asX78eqVTK8ePHeeKJJ7xm09ixv2Xbti1s27YFvV7L2LG/rSWwfV5+ncxd28Fi5vjiBQDINRrkVwTUYrFgrqqk6GQ6ReknCOncGW3XbshUaq/ZLBAI/IsWFdacnBymT5/OggULmDZtGm+99RZLlizh0KFD3H///dxyyy0AvPnmm8yfP59ly5Zxyy23cO+993rNpk8/Xc/999/J5s2bAdi2bUuT5xyZP6feNm23eDr9fjqlZ89Sdv484Qm9CGrXHok0MD+GCa4fKioqeP/9v7Jz5zYKCwuIjIxm0KAhPPjgw7Y1kIqKCj7++AO2b99CdnYWWm0Yw4ePZNasRwkL09W6nsViYfLkCWRnZ3HrrZOYO/f5lnhbPkVyPQwTdDUUYP2Yp9dr6+0bO/a3vD7rEbCYOf/ZJxSfPtXgdbTd4om77wHMBgNVBQUodeGE9+qNKjLSrffhrN2BRCDaDM7ZnZl5nlatOvjIIudw5iP1okXz2L59K5GRUXTu3JWffz5LdnYWN9zQi3feeR+j0cjjjz9EevpxdDodnTv/hoyMHyksLKRTpzjWrPk7avXVT2jff3+QJ554GICQkFD++99vULq4DtFSoYCGfoZ+HQrwd6xhAXu2bdsCsx4BIO6+Bxyed/ajtRSfPmUTXalSiSYmBmNpKVkpOwnu0AldzwRkKpV334BA4CJlZaXs2LENnU7H559/gUqlpqqqkmnTpnL8+DHOnTvLnj27SU8/TkJCL1599S1CQkIoKSnh4Yenc+7cWTZt+i+TJ0+xXXPr1q8BUCqVlJaWsGdPCklJY1vqLfoE8bm0ET79dD3Z2cW2LyuPLV/a6HkNCa48JAR1dDQVly6QtWs7Vfl5HrVXIGguMpkciURCSUkJmzd/SUVFBSqVmtdee5N33nmf6Gg9X31Vky45c+YjhITUeG2hoaE8++w8/vSnVxk3Ltl2PYPBwM6d25FKpdxzz30AfPPNJt+/MR8jhNUFxo79LQCphw/x6J9ewlhehrm62qVrSKRSVJGRSGQysnfvpOTcWY+klAkEnkCtVjNy5GhMJhOvvfYy48eP4ZlnnuD48WN0794DuVzGxSspid26da91bp8+/Rg58kab2ALs37+X0tISevXqzcSJv0MikZCWtt+W/XOtIkIBLvDpp+u5++7b2bZtC3uOHiHx7ikMT+zNm48+gTIiAqlC4fS15EFBSJVKCg4fwlhcTFjPBKRy8eO4lqmbqucLogYPpc/Lr7t0znPPLSAkJIQtW77CYDCQlraPtLR9fPHFev70p1dtxwUFBTV5ra1bvwJg5MgbiYqKpkePBNLTj7N9+zfcfvudrr2ZAEL8JruIvbgC7Dl6hH4Pzah3XO+QEOa0b99gxkDcfQ8glctR6/WUnMvAWFZCRL+BIu4qaHFCQkJ47rkFzJ79FPv372Hfvj3s3LmNH344yc6d223HlZWVERracPJ8aWkp+/btAWDEiBsBGDVqNOnpx/n6681CWAW1sea12gtsXY6UlrLil18clsYWnz7F2Y/WEnffA0ikUjR6PVX5+WSn7iJ6yHDkwe41fhD4N656ji3B8eNH+de/PqFbt+78/vfTuemmm7npppvp1CmO1avf5ZdffiYmphVZWZmcPv0D/fsPtJ37j398yHffHeD3v7+f/v0HsmvXdgwGAwB33PF/te5z6tRJzp//mQ4dOvry7fkMEWNtBnUXt7Kyijh36ATDeyUCNeLa+08ran1pu8UDV8XViioiAouxmqzdOzBcqTITCHxNUFAwu3bt4LPP/kFOTrZt++XLlwGIiorm5pvHA/D++3+lrKxmDlVubg7r1n3GoUMHqK6uEVNrNkBMTCs6dYqzfYWG1qQxfv31tbuIJTxWDyKRSAhu157P//s1bTq2BqDnbeMdHts7JIS6QQKFNgxjeRnZqbuIGjIMdVS0ly0WCGrTuXMXRo0aze7dO5k27Q569EggNzeHc+fOEhqqJTn5VrTaMPbtS+XEiWNMnTqJzp27cvr0D5SWljBs2AgGDRpKbm4Ohw8fQiaT8d57HxIVFWW7x1dffcmyZS+yZctXzJr1qN/1UvAEwmP1AvKgYMY0kad3pNTxxEl5UDDy4GBy9qZQcflXb5gnEDTKCy8sYfr0mURGRnH06BHy8vIYNmwEK1euRq+PQa1Ws3Llau68cxpKpZLjx4+g04UzffpMlixZjlQqZevWbzCbzQwaNKSWqAIkJY0lNFRLVlYmhw8faqF36V1E5ZUDPFUNZKqoICtlFxJJTQ6rFasX+0mPHsDVxaxa5xoMGAryiew/iOB2Dbcw9IbdviQQbYZru/LKHwm0yivhsXoRmUZD9JChmCorMVU1PE7GUVmsTKlEFRFJ3sE0Ss6d9aaZAoHAwwhh9TIKbRiRg4diKCrCbDTW2vee5moe4JH5c2otZgFIFQpUkVEUHD5EccZPPrFXIBA0HyGsPkCjjyE8sS9VublYzGZG9usPQMqhg7YsAXDsuUrlctTR0RQePUzxj6dElZZAEACIrAAfEdopDlNZGSVnTvPuCy+SMGkCALdu+I/tmN4hIcyZP6dezFUik6GO1lN44jgWkxltfPdrciVVILhWEB6rDwnr0RNNm7ZU5eXZvFZ7jpSWcs/Jk9y64T/0vG08j7y0yLZPIpOh1uspOpVO0ckTwnMVCPwYIaw+RCKVEtGnH4pQLW89+QzpX2yyfTkS2pRDB2sJrEQqRR2tp/j0KYrST2AJwDHGAsH1gBBWHyNVKIjsPxBzdTWmK+V+AKsWLCb9i0180qMHn/ToQW+79KyUQwdt/5dIpaj1eop/PEVh+jEhrgKBHyKEtQVQhIYS0W8Ahvz8esJoXcya0769Lc+1LlZxLfnpJwqPC3EVCPwNsXjVQgS1iSW0S1dKz51FHX21dNW6aGWdQmDFWlQwsl9/Vi1YfFVcz57BYjYRntjHt29AIBA0iE891o0bNzJq1CiGDh3K+++/X2//Bx98wPDhwxkxYkSt/evWraNbt262r4ceesiXZnuNsB4JyENDqS6pX8ETd98DNY1bQmpXd1jjrj1vG0/C726l/0MzeOCZJyk4chizyeQr0wXXKMOH92/w6/PPP/XafdeufY/hw/vz8ccfNnjM5cu/Mnx4f+6+e7JH7rls2YsMH96fbdu+8cj17PGZx5qZmcnzzz/PpEmTCA8PZ8WKFQwYMIDExJpOUAcPHuTll19m1qxZqFQqVqxYQVxcHElJSRw7doyuXbsyb948JBIJERERvjLbq0jlcqL6DyRz53ZkarXDRtnWtoO9/7SCR15aVCveamXviWOUnT9H9ndqpB3jkchkXrddcO0ilUoZNmxEve1t2rT12j07dOjIiBGjaO+gzaYVtVrNiBGjiI7We80OT+EzYd2/fz/V1dXMnDmTmJgY1qxZQ2pqqk1Yjx49CsCMGTMIDw/nH//4B1u3brUJ66+//sqsWbOIj4/njTfe8JXZXkehDSM8sQ/5hw+h1usbzU9dtWBxvW3WEIFKr6f43DnMBWVE9OknxFXgNnK5nOXLX/PpPceOHcfYseMaPSY8PMLndrmLz0IBWVlZAOh0OlQqFcHBwWRnX+332Lp1TZu9PXv2cPLkSUpKSsjOzsZsNhMcHMyECRNYuXIlWVlZLFy40Fdm+4TgDh0Jim2LoZE5QEfmz3FY9mrl0aUvoomJoeziL+R9f7Be+axA4Clmz57F8OH92b9/D/feO5UxY4axZMkCcnKy+eMf/0BS0jDuv/9uzpypKcO2foR/+unH+de/PiE5eQwTJtzEP/7xoe2adUMB1nt8/PEHjB8/hlmzZjgMBZw//zNPPz2bMWOGMWHCWJYuXUSpXee4jz5ay+TJExg9eggTJ97MO++8idkHi70+81itCe0NeWTjxo1j+PDhPPvss4SFhdlGPkilUj799GpsZ/z48XzyySeYzWakUuf+LjTWhaYhoqMbHjnhDSKShnP+66+RKUCu0Vzd3rMn+enpttfFp0/xyycfkvBIzQjupEGD2JGWRsqhg7QfM4akQYNY9eSTmDLSiR48OCDmaPn6e+0pmrI7O1uKXO5/iTdN2WQ0Gpk//9l62xcvXoZGo7H9Di9aNJ/4+B78+usltmz5in37UomKikav13PmzI+89trLrFnzATJZzf1OnDhGevpxunWL5/jxo/z1rytp3749SUljkUprrimTSZDLpbZ7/O1vq0lM7E23bvG260gkNcdUVlbyzDOzyczMpHv3HpSXl/P115uorKzg5Zdf5ZtvvmLNmlWEh4fTt28/Tpw4wT//+TG9evVi9OgxtnvIZA3/nKRSqVvPp89+66KvrHwXFhaiUqkoLy9Hr78aK5HJZLz11ltcunSJNm3acMsttxAbG0tBQQHr1q1j5MiRxMfHYzKZkMvlTosqtFzbQFdR/OYGcvamoo6ORnLl/bW/536sUSdrpkB+ejrFxRUAvD1vYa3Y6460NB587XXeeOgxCgvKiew/wK/F9VpuG2g2m2u1umtslI+3GDv2t7ZRQuBc+z2z2UxKyq562ysrDSgUKpuTNHXqPTzwwEN8+OH7vP/+X9HpIvjww3+Sn5/H7343nrNnz2A0mjGZau5XUVHOmjV/Jz6+O1u2fMWSJQv4/PPPGDkyyfb7aTJZMBrNtnvcc899PPjgw8jlUi5cuAjUOGlGo5mvvtpMZmYmN910M4sWLaW6upq5c58mNDSMykoDbdq0Zfr0mdxyywTatInls8/+wcqVb3D+/Pla9zCZzA1+T8xms8Ofc1NtA332Gzd48GBkMhmrV68mIiICk8nEoEGD2LdvH3FxcVRXVzN27FimTp1KbGws2dnZjBs3Do1Gw9q1a9m1axfTpk1j06ZNjBvXeCwmUNHoYwjt2pXSjIxaKVhW4u57wOFwQmvs9fHlS9iRlkbqoYOo9Xoqsy6T9923RA4IDM9V4B8olUp27Gh6mmyfPv2AmnEtAD17JiCXy22vq+uMho+O1hMfXzMye/jwkQBcuHC+0Xv07HlDg/t+/rmmnWa/fgMAUCgUvP7627b93bv3pLKykn//+18cO3aE01fSFw12hTnewme/bW3btmX58uW88cYbVFVVMWfOHFq3bs2YMWNYsGAB06ZN46mnnmLt2rWEhISwcOFChg0bBsBbb73FkiVLeO655xg6dCjz5s3zldk+Jyy+J5WZmRhLS2s1x3aGj5Yvp11SElCzqDWyX3/+8shs8g7srxFXF8ZzCzyPved4LaBU1kwUtn56VKvVtV7XxV5or4YGG//kGezEYE2TXZphVVUlKlWNHVYPtV+/AUyadAfZ2Vm8//5ffdJnw6duzMSJE5k4cWKtbadPn7b9/+GHH+bhhx+ud96gQYPYtOnaHTxmj1QuJ7L/QLJ2bUem0TS4um/vudp3wxrZr78tLJBy6CBPrVrJXx55nNy0/UQNHIxUqfT+mxAIHFBQkM/Ro4dJTOzDt9/uB6B9+8YnLEilDWe3dOjQCYDvvvuWiRN/h9FoZNq0qRgMVXz88b/YtOm/ADz55DPExXXh1Vdf9tA7aRrx+dAPUerCCUvoRdGJY6j1MbX2abvF1+vbav/aGhawxl1TDh1EHRVFZV4eOWn7iBo4BJlK5f03IQhYjEYj8+Y9U297374DuOOOO92+rkwm49lnn6B7956cOHEcgClT7nL7ejfddDN/+9tqdu3awQMP/J6qqkouX77EsGEj0GrDaNWqNefOnWXOnD8QFRVNenrNPSsqKty+p7M0KKwXLlxw+iLt2rXziDGCq4R26kzl5csYigpRhuls2+vOxrJ6rmc/Wkvvx2fbtq9asNiW4/rIS4tYtWAxVXn55Hy7l+jBQ5Fd+bgkENTFbDaTmrq73vaQkOZlb8TGtuXWWyfx8ccfEBQUxMyZjzB8+Ci3r6fRaHj99bd5883XOH78GEFBQdxyywQef/xpAJ544hlKS0v58cdTKJVKHnpoNu+9t5L09GPNeh/O0OAwwfj4+CabKVssFiQSCT/88INXjPMUgZIVUJfq0lKydmxFERbWYHzUvqfAyJUrbdkCQK1sgfQvakIpVQUFyDVBRA8Zhswuraul8JfvtauIYYLOc/nyr9xxx//Rvn0HPv30325dI9CGCTbosW7fvt0zlgncRhESQniffuR9l4Y6JsbhHzr7TIGU2bNrxVvtvVYrqvBwDAUFZO9NIXroCORBQfWuKRAImkeDwhobG+twu7VqwWKxYDAYOHnyZIPHCppPUNt2VGRepjIrE1UDPRLs466O5mbVRRkejqGokOw9u9EPG4nciZVXgUDgPE5l2R88eJAJEybQvXt3evbsSc+ePUlISKBv377cf//9Xjbx+kYikRB+QyISmQxTpeMR2tZOWFYclb3aj3kBauK2ZjPZe3ZTXVLsWaMFAjtat27Dnj0H3Q4DBCJOCeuyZcto3749q1evRqPR8Pbbb/PCCy+g0+n485//7G0br3tkanVNY+yiokabWkf07AnUeK3WvgL2E2HriqtCqwUgK3UXhsJCL1kvEFx/OCWsZ86c4ZlnnmHEiBH07NkThULBPffcw6JFixz2VRV4HmtVVlVeXoPHJDzySL1x2nPat29cXENDkSmUZO/ZRVV+w9cWuIcY+hi4NOdn55SwajQaWzVFXFwcp07VxPF69erFuXPn3L65wDXC4nsiDw3FaNe9py7WsIBVYB2Ja90JsPLgYGQaDdl7UqjMyXZ4XYHrSKUyTCbRZSxQqa42IJO5l+rvlLAOHjyY1157jaysLPr06cOmTZvIy8tj69at6HS6pi8g8AhSuZzIvgMwVpQ32RYw7r4HaonrqgWLa02Creu9yjVBKEJCyNmbQvmvl7zzBq4zNJoQSkoKsVjETLJAomZhvorCwhxCQtzTtwbzWO3Jycnhj3/8I0lJSdx5551Mnz6dgwcPIpPJePHFF7njjjvcurmvCNQ81oYoyThDwbHDaGJa1dqu1Wpq5bHC1QIC+zQsR5MIrLO0TAYDhoJ8IvsOILhDR++9iSv4+/e6IZyx22KxUFCQg8FQCfhHSEAqlfqkH6mn8bXdMpmckBAdGo3jjJmm8lidElZrIYD96zNnzqDVaomJiWnkTP/gWhNWi9lMbto+qvLzUYWH27Y7Elb7AgL7zAFH4motIjBXV1OVl4uuV29CO3dtslCkOfj797ohhN2+xd/sbkpYnQoFjBgxgqVLl3L48GGgJgWoa9euASGq1yISqZSI3v1AAqYqxylYVuxLYO3TsFYtWEz6F5tsYgpXU7KkCgXqaD2Fx45QfOqkWIARCFzEKWGdP38+2dnZTJ8+naSkJF555RXS7braC3yPTKMhsu8ADAWFjaZgAbVirY5yXO0XtqxIZDLU+hiKT/1A4dHDWMQEWIHAaZwKBVipqKhg586dfPPNN+zZs4eoqCjGjx/PE0884U0bm821Fgqwp+D4EUrPnkUdHe0wFGDFPiQAtWOuQL3SV0eMSbqJf37m2STvQPpe2yPs9i3+ZrdHQgFWNBoNycnJPProo9x7773k5OTw0UcfNdtIgfuEdU9AHhpKdUnjD519lgDUL321zxhoiO07tqLXa9Hrtdx99+3uGSwQXAc47bEeO3aMLVu2sHXrVrKyshg5ciTJycmMHj0alZ/397yWPVaA6uIiMnduJ6p9a0orms6bdJQp0BSOFrvqzlNyh0D7XlsRdvsWf7PbIzOvbrzxRnJzcxkyZAiPPPIIY8eOJcTFsSEC76HQhhGe2IfKH49jCdY1uYpvbdriTMMWK9YG2tXFRcx+9RX2HDvi88F4AkGg4FQo4OGHHyY1NZU1a9Zw2223CVH1Q4I7dCS0UycM+flNHuusl+oIhTaMd56bb3stwgICQX0a9Fj379/PgAEDkMvldOjQwVbG6oghQ4Z4xTiB80gkEvT9+pF99iLG8nKn+6y6ExaQBwUzok9fUg9/D8C2bVu4++7br7lheQKBuzQ6QWDv3r1ERkYSHx/v6JCaC7gwQWDjxo28/vrrVFdXM2PGDB588MFa+z/44APWrl2LRCLhvvvus+0/efIk8+bN4+LFi4wdO5bFixfbJkI6w7UeY7USHR3KxdM/k7V7J6rIyEZHXtfNEgDXxBXAbDTy8IsvsPfK/CJ3Yq6B/L0WdvsOf7PbI5VXniAzM5OxY8cyadIkwsPDee+99/j8889JTEwEanq+3nPPPcyaNQuVSsXbb7/NqlWrSEpKIjk5Ga1Wyx133MGiRYt48sknmTlzptP3vp6ENSenhOIzP1J4/Gi9kteGsBdZV8XVYjbz0ML5NnEF1wQ20L/XgYaw2zN4JN3q9ttv58MPPyQrK8ttQ/bv3091dTUzZ87kscceQyqVkpqaatt/9OhRAGbMmMHs2bMJDw9n69atXLp0iYyMDCZPnszkyZNJSEggJSXFbTuuB0LjuqBpHUtVQYFTx9dt2OIKEqmU915azrAbetm2iUUtwfWOU1kBv/3tb/nf//7HK6+8Qp8+fUhOTuaWW24hooFRIY6wirJOp0OlUhEcHEx29tUWda1btwZgz549dO7cmZKSErKzs23nhV+piY+IiCAjI8Pp+wKN/mVpiOjo5k2kbCmsdkeMHckv33yDRG5B4US8tffjs0mZXTPlVat1fcjgZ2++QVVBAV0mTwbg/vvvZNOmTU2cVdvmQEPY7VsCyW6nhHXWrFnMmjWLCxcu8PXXX7NhwwaWL1/OwIEDGT9+PJOv/DI1hjXi0FAq0Lhx4xg+fDjPPvssYWFhhIaG1jqvOVxvoQAryvjeZKXsQB0VjUQmc/o6dYcSOo1MzYjefUg9cpjNmzc79T28Vr7XgYKw2zN4tPKqXbt2zJw5k/fff58nn3ySI0eOsGDBAqfOjY6OBqCwsBCDwUB5eTl6vd62XyaT8dZbb/G///2PHTt2oFAoiI2NtR1TeGV0SGFhYa3zBA2jioxEd0Milbk5Tv2BqluZZR3v4gp/fXGp7f/WKi2RkiW43nC6PXZubi5bt25l69atHDhwgG7duvHYY48xfnzTNeZQ0yxbJpOxevVqIiIiMJlMDBo0iH379hEXF0d1dTVjx45l6tSpxMbGkp2dzbhx42jXrh3t27dn/fr1SKVSjh8/7ve9CfyJ0LguVOXlUZWdjSqy8dCN1UO1X8xyNeYKNeWxdau0REqW4HrCqayAu+66i6NHj9KpUyeSk5MZP348HTt2dPlmGzdu5I033qCqqooZM2Zw8803M2bMGBYsWMC0adP461//ytq1awkJCeHBBx/knnvuAWrSrebPn8+FCxcYM2YML730kktltNdrKMCKqaqSrF07kcgkyIOcH3VtzXG17+PqChaLharcXP7w7tukHj4EQHZ27Ymw19r32t8RdnsGj6RbPfHEEzz66KON5rP6M9e7sAIYCgvI2r0DpS4cqULh1PWaK6xX711In/tr/kjWTcW6Fr/X/oyw2zN4JMb63XffYWxixpLAv1HqwglP7ENVXp7LC4JH5s9xK9569d46RvTuA9SEBETcVXCt45SwtmrVqlk5rAL/ILhDJ4I7dqIqL9ep4+0Xs8C9eKuVv764lBF9+tXaJvJdBdcqToUCnnvuOf73v//RvXt3YmNj68U3X3nlFa8Z6AlEKOAq5upqslN3YzZUodBqnb62NSxgxa10LGoqtSpzc+k/azoAycnJfPjhZy5fp6W5lp8Rf8Tf7PZIKEAikfB///d/dO3alaCgIGQyWa0vQeAgVSiIHDAIs7EaU1WV0+d5yowjvvcAACAASURBVHuVSKVo9HqG96opZd68ebMICQiuOXzWK6AlER5rfSoyL5Ozf4/LxQPgXkcsRzz84gukHqkZUJmVVeTVabCe5np4RvwJf7PbI1kBb775ZqP7n3zySdct8yFCWB1TdPoHitJPoI6JcUnUmpqf5QrWWVujBg7iXxs2O52x0NJcL8+Iv+BvdntkgsDBg7WTvU0mExcvXqS4uJhbbrmleRYKWgxt124YCgupzMpCHRXp9Hlx9z3Q7CICK0mDBrEjLY3dB9LITt1N5IBBKEIDpyZcIHBEs0IBf/7znzEajTz//POetMnjCI+1YUxVVWSn7sJiMrklaPaLWu54rlqthnZJSTXX+vs/MRsMRPQfSFCbWJdt8SXX0zPiD/ib3R7tFVCXu+++mw0bNjTnEoIWRqZSETVoCOaqKpcWs6w0NvnVVZ548zXkWi253+6j4MQxzCJ3WhCgNEtYd+7c6VInf4F/ogjVEjlwMIaCAiwmk0vnxt33QK3KLHeKCKyjt1MOHaTXlEk8tXoVJWd+ImdvCtWlpS5fTyBoaZyKsY4aNare4kZZWRmlpaXMnTvXK4YJfIumVWt0NyRSdOIYKr3e5RV6+8mvR+bPcSkssGrB4lrjtVO/P4hGr8dQVETWzm1E9O1PUGxbl9+TQNBSOBVjrftxXyKRoFAoSEhIoEOHDl4zzlOIGKtzWMxm8g8fovziBdRX2jy6Qt1sAWd6DGi1GoqLK2yvrVkCI/v1Z9WCxZgMBgz5+YR27kpYz4RG53j5kuv1GWkp/M1uj2QFTJo0qdZrg8HAqVOnXJogIPB/JFIp4Yl9MJaWYCgsRKnTuXS+1UOtW6XlCtaWg1bvVaZUotbrKf35LJV5OUT2G4gyLMzt6wsEvsCpGOu5c+e48847+f777yktLeV3v/sdU6ZM4cYbb6yXiiUIbKRyOZEDBoNEgrG8zOf3X7Vgcb1tEqkUdXQ0FqORrF3bKTmbgcVs9rltAoGzOCWsL730Enq9no4dO/Lvf/+b/Px8du/ezYwZM3j55Ze9baPAx8iDgogeMgxjWTkmg8Ht61i7YrnbGeuRlxbVeq0ICUEZHk7h0cPkHtiPqaKigTMFgpbFKWH9/vvvmTNnDhEREWzfvp0xY8YQExPDbbfdxk8//eRtGwUtgFIXTuSAQRgK8l3OFKjbVwBqUrGcFVf7LIGet42n523jbSIrlctRx8RQlZfH5R1bKb/8q0u2CQS+wClhDQoKori4mPz8fL7//ntGjRoFwM8//4zOxTicIHAIahOLLuHKzCwXPnpbU7CsX/ajtZ0R11ULFtvE1UrKoYO1PFhVeDjyoCBy9++l4Mj3mJvhWQsEnsYpYb3pppt46qmnuO+++9DpdIwcOZLNmzfzxz/+kYkTJ3rbRkELEtqlK6FxXajMzXH7GnH3PeCWuKZ/sYn0LzbV8mDtkalUNQtbv5wnc+d2p/vMCgTexilhXbBgAXfffTeDBw/mo48+QqlUYjAYeOSRR3jqqae8baOgBZFIJIT1vAFNTCuq8vLcvk5dcXUFRwtaNvukUtRRUUjkUrJTdlJ44rio2BK0OD5tG7hx40Zef/11qqurmTFjBg8++GCt/V988QV/+ctfKCkpYfTo0SxfvhylUsm6det44YUXbMfdeOONvPfee07fV+SxNh9TVRU5e1MxVVU2K92pbipWRM+etL/n/ibPs+a3pn+xqcFjLGYzVXl5yIODiew3AGW499IBxTPiW/zNbo/ksXqCzMxMnn/+eSZNmkR4eDgrVqxgwIABJCbWNDyuqKjghRdeIDk5mQEDBrB48WIGDRrElClTOHbsGF27dmXevHlIJBKRP9sCyFQqogYPIWv3ToxlZciDnZ/2ao+1QstKfno6xo/WOl2lZRVYK9ZCArialmUsKyNz9w7CunUntGs3vykqEFw/NKtXgCvs37+f6upqZs6cyWOPPYZUKiU1NdW233Rl5blbt2707dsXhUKBUqkE4NixY/z666/MmjWL1157jWA3f6kFzUMeFEz0kOGYKiowVVa6dQ37hS1X4q51F7Os1F3UApAHB6OOjKL49Cmydm2nKt/9EIZA4A4NCuuZM2cwezAJ2zqMUKfToVKpCA4OJjs727Y/JCSEefPm8corr5CcnEx8fDy33norZrOZ4OBgJkyYwMqVK8nKymLhwoUes0vgGkqdjsjBQzEUFWGurm7WteLue4CInj2BpsXVfjHL0aJW3bQsiUyGWq8HLGTt3lETe22mvQKBszT4GWnKlCls2rSJ1q1bc++997Jy5Uq0Lgyfq4s1lNtQc4/8/HxWrlzJlClTGDZsGHPnzmXVqlXMnj2bTz/91Hbc+PHj+eSTTzCbzUilzjncjcVCGiI6OjCbLfvE7uhQdCEKLu/Zg0YXg7QZc88SHnmEE6tWkZ+eTvHpU2i1GqfP/WTFK9w3bx470tJs21IOHeTx5Uv4aPnymg1aDRZ9BJXZFygvyUU/cCDBrVq5ba894hnxLYFkd4PCGhwczDvvvENiYiIHDhxgw4YNDX4Ev/32pofBRV9p6lFYWIhKpaK8vBy9Xm/bn5aWRn5+PnfccQe9evVi7dq1fPvtt9xzzz2sW7eOkSNHEh8fj8lkQi6XOy2qIBavvEJQBPK4eHKOHUWt1yNx4edhj1arof0995N/ZVErZfbsmu1Odsd6e97VTy/WDlk70tK4Z84fa2cTqEOpKK8gf+NmQuM6o+3eA5nK/ZaX4hnxLf5mt9uLVwsXLuSdd95h//79SCQSPvjgA4diJpFInBLWwYMHI5PJWL16NREREZhMJgYNGsS+ffuIi4ujS5cutvuMGjWKH374gbvuuguNRsPatWvZtWsX06ZNY9OmTYwbN87Jty/wJqGdu2KurKL4p1Oo9a7NzapL3UUta2jAlYkE9u0H6+a8Asg1GmQqFWUXfqHs14tEJPZF0yY2oIYYCgIDp9Kt4uPj2bt3L5GRzs9FcsTGjRt54403qKqqYsaMGdx8882MGTOGBQsWMG3aNNavX8+7775LQUEBI0aMYOnSpWi1WtLS0liyZAkXLlxg6NChvPzyyy5VfAmP1XtYzGYKjh6m9Pw51NFu9HGt0zYQarcfdKb1YF0ayxywYqqqwlBYgKZ1G3QJiShCXAsXiWfEt/ib3R6Z0molIyODjIwMTCYTcXFxdOvWzSNGehshrN7FYjKRe/AAlZmXXe7j6khY4Wq+qzvCat8020pD+a+GwgLMRhNhPXoS0jHO6dQs8Yz4Fn+z2yPCWlxczJw5c9i9ezdhYWGYTCZKS0sZMGAA7777LqF+PlVTCKv3MRuN5Kbtx5Cfh8qFTzZNCavtODdHbNf1XqG+B2s2GjHk5yMPDSWid19UkVFNXlc8I77F3+z2yDDBJUuWkJuby+bNm0lLS+PgwYN8+eWXVFRUsNy6+iq4rpHK5UQNGIRCq6UqP7/Z16vbIcuV7lj2OMp/revNSuXymtQsi5mslF3kfX9ItCQUNAunPNZ+/frx0UcfkZCQUGv7sWPHmDlzJml26S7+iPBYfYepspLsfamYK6tQ6poufW3IY7Wn7sgXd71XaLo01mI2YygsBECXcAPB7TogcZBO5g/fa3cQdnsGj3isDU1ilUgktoopgQBAplYTPXgYUqWC6uIij1zTvoELNH/MNtRvom1FIpWiiohAERpKwZHvydy9Q3TNEriMU8KalJTEkiVLOHfunG3b2bNneemllxg9erTXjBMEJjUTCIaDREp1SbFHrll3zLa71K3WakhgpQoFan0MmE014YFDB1pkVI0gMHF68eqxxx7ju+++I+RKWkpZWRmjRo3ilVdeaVZFli8QoYCWobq0lOw9u5FIpQ2mMzkTCrDH0aBCV0MDjrIGHKVkWbFYLBgKCsBsRtu9JyEdOxHTJsKvvtfO4m/PiLP4m90e6W6l1Wr5+OOPOX36NBkZGajVajp16kSnTp08Zqjg2kMREoJ+2AiyU3djLC1F7mKuqCPqFhKA+/1d7QXWUUGBFYlEgioiArPRSFH6cUozfiJo5BAsap3bFWeCaxuf9mNtKYTH2rIYiorI3rMbmUqJPKh2WbSrHmtdrB6sJxa0oHHP1YqpqhKlsRKDIpjwGxJdSi9rafz1GWkKf7PbI4tXAkFzUIaFoR82ElNlFcbyco9e2771oLvTYO1TshrzXK3IVGqCW7fGYqgiK2Unud+lUV3iP7/0gpZHCKvAJyh1OvTDR2KqrPCouDrKGHBVXK0tCV1FHhKCWq+nKieLzB1bKTxx3O0+tYJrC6dCASaTCVkzWsO1NCIU4D8YCgtqwgJqDfKgoGaHAuxpbr6royote+zDBHXttphMVBUUgATCuvckpEMnpAqFi+/A+wTCM+IIf7PbI6GAoUOHsmjRIg4ebPpjkkDQGEpdOPrhozzuuULz810bmlJgxdG0AisSmQx1VBTKMB1F6Se4vO0byn45j0XkeV+XOOWxpqam8tVXX7Ft2zY0Gg3JycmMHz++XiWWvyI8Vv+jxnNNQRcVRoXZ85+GPNVrwIp9BkH6F5ua9LRNBgPVhQXIQ0JtU279IYMgkJ4Re/zNbo92t6qurmbPnj188803pKamEhISwoQJExg/fjxxcXEeMdgbCGH1TwyFBVQcPUC5weL2cMKGqBsWAPc6ZdljXw7rbAjDVFlBdXExcq0OXc+EZvetbS6B9oxY8Te7PSqsAEajkb1797J9+3Y2btxIeHg4xcXFdO/enRdffJGuXbs222hPI4TVf9EqTJz8YhMyhcIjea6O8GanLHAuRctYXk51STGqiEjCeiSgiopuEYENxGcE/M9uj8RYjUYju3fvZt68eQwdOpR58+Yhk8lsnf337t1LbGwsjz76qMcMF1wfqHQ6YkaMwmw0ei1lyVGnLHdobFJsY+WxUFPmq4lphbmqiuw9KeTsSaEyN4frII38usQpj3XAgAGYzWaSkpKYMGECw4cPr5cl8PXXX/OXv/yFb775xmvGuovwWP0Xq83VJcXk7NsDFgsKL5ZI23uvzYm7WkMBjspjrTTmyVaXlmIsK0UdpUfbvQeqyCifeLCB+IyA/9ntkVDA5s2bSUpKctjlKj8/n4iIiOZZ6WWEsPov9jYby8rI3puCxWhE6cLoHVfwVAvCujHWhgS2qTBBjcCWoYqIIKx7Aqpo74YIAvEZAf+z2yPC2r17d/bu3VtPQC9evMitt97K4cOHm2+pFxHC6r/UtdlYXk7O/j2YKipQefEPdnMFtqnFq7pZBE1hLCvDWFaKIkyHNr6H17IIAvEZAf+z221h3bBhA+vXrwfg0KFDJCYmIq8zDygnJweJROKXH//tEcLqvziy2VRZSU7aPozFxV6tw3eUOeCswDqTFeBqDwKo+cNiLC1BHhKKtlt3NK3bOD2HyxkC8RkB/7Pb7cWrcePGMWTIEIYMGQJA3759ba+tX3fddRd/+9vfnDZm48aNjBo1iqFDh/L+++/X2//FF18watQo+vbtyzPPPIPBYADg5MmTTJw4kX79+jF37lwqRdngNY21WbYqMpLKHO8t8Fh7vHq6ibaVuj0IGlvcsiIPCqpJyZJKyD90gMtbv6b03FnMV34XBIGBU6GADRs2MH78eJRKpds3yszMZOzYsUyaNInw8HDee+89Pv/8cxITEwGoqKhgwIABJCcnM2DAABYvXszChQuZMmUKycnJaLVa7rjjDhYtWsSTTz7JzJkznb638Fj9l8ZsNhuN5B8+RMXFCzWxRy8n2LuSluVKKa59WMBZz9VKTaFBIRKZjNAuvyG4Q0fkQUFOn1+XQHxGwP/sdttjXb9+vc1jNJlM/Pe//2X9+vUOv5xh//79VFdXM3PmTB577DGkUimpqam2/dYRL926daNv374oFAqUSiWXLl0iIyODyZMnM3nyZBISEkhJSXHqnoLARiqXE9m3PyFxnanKyfF6eain0rLqsmrB4lqTC5zxXK3IlErUej2KsDCKf/qRy1u/Jv/wIdtcLoF/0mDw5t1332XMmDEolUrefffdBi8gkUi4/fbbm7xRVlYWADqdDpVKRXBwMNnZ2bb9ISEhzJs3jyVLlvDKK6/Qt29fbr31Vo4ePQpAeHg4ABEREWRkZDj37q7Q2F+WhoiO9u+R3g0RiHY3ZXN0zAjyYsLJO3YMjV7v0ZijPb0fn237f8rsmv9rtZoGj29sX10+WfEK982bx460NFIOHXTpXBsRoVjMZqoK8yg7eBmzXk9Ejx4ExcS45M0H4jMCgWV3g0/ojh07HP7fXawRh4ZSSfLz81m5ciVTpkxh2LBhzJ07l1WrVtlivM1BhAL8F6dtbtURRbmJnCPfo4qIQNqMsJQrWAUWaocG3OnK9fa8hbYFrXZJSYDroQEA5BosMjUVmXlkZ3yNPDgY7W+6oWkdi0ylavTUQHxGwP/sdns0y/79+526gUQiYfDgwU0eFx0dDUBhYSEqlYry8nL0er1tf1paGvn5+dxxxx306tWLtWvX8u233zJx4kTbedZ/7c8TXD+ExnVGplaT910a8pAQ5Bo3vD4naWgEzNmP1jarmcvIfv1r5bs601jbERKJBEVICIqQEExVVRQcOUzB8aOEdupMcPsOKLRNjx4XeI8GhXX69OlOXUAikfDDDz80edzgwYORyWSsXr2aiIgITCYTgwYNYt++fcTFxdGlSxckEgkffPABo0aN4ocffuCuu+6iXbt2tG/fnvXr1yOVSjl+/DhPPPGE8+9QcE0R1CYW2YhR5Ozfi8VYjSLUO1VadcXTmprV3LirvXfaVP9XZ5GpVMj0eiwmE6U/n6PkzE+ooqII7doNdVQ0kgDupRyoNCisp055Lu0EoG3btixfvpw33niDqqoq5syZQ+vWrRkzZgwLFixg2rRpLF26lHfffZddu3aRlJTE448/DsCbb77J/PnzWbZsGbfccgv33nuvR20TBBaqiEhiRo4mN20fVfn5Xi0ksBJ33wMOJ8T6ExKZDFVEBBaLBVNZGbn79yJVqQjt8huCYts2K5tA4BoNpltduHCBtm3bIpFIuHDhQqMXadeunVeM8xQixuq/NMdmU2Ulud+lYcjPQxXl/Vr7hoS1ud2y3BkL4yzm6moMRUVgNqNp04YO/RMpRe0XvWFdwd+ebbcrr+Lj49m7dy+RkZHEx8cjkUgcJmo7GwpoSYSw+i/NtdlsNFJw9DBlv/zs9Y+9jiq17HF3FIxbC1guYjGbMZaWopFbqDBKCOnSlaA2gePF+tuz7bawXrp0iTZt2iCRSLh06VKjN4mNjW2elV5GCKv/4gmbLWYzxT+epujkCVSRkT6ZNWXNCmhOWayr/QQ8gVaroSC3iOri4hovtlVrQuI613TX8uNYrL892x5tdH3mzBkyMjJQKBR06dKF9u3be8RIbyOE1X/xpM1lFy+Qf+jAlYwB73piddOt3BVYd/oJNAd7uy1mM8ayUkwVlUiVSkI6dSaobVuvLQg2B397tt1Ot7Ln0qVL/PGPf+TQoUOEhYVhsVgoKSlh9OjR/OlPf0LnpRZvAoErBLdthzwoiJxv92I2VKMM813KkVVA7QXWmfQs+/Qrd1Ov3EUilaII1aII1WKurqbkp9MUn/4BpS6ckM5dUOtjmsyLFTjGKY91xowZmM1mli1bZvvY//PPPzN//nyioqJ46623vG5ocxAeq//iDZuN5WXkpn2LsaQYZWSkVxa1mioQcLUtoa88V2cKG4zl5VSXliKRSgmKbUtw+w6oIiJbNFTgb8+2R0azHDx4kBdeeKFWLLVjx44sXLhQ1O0L/A55UDD64SPRtImlMjurRUZQOxrFffajtQ0eX7cTVksiDwpCo9ejioigMjuLnL2p/PrNZgpPnsBQWCDGyTiBU8LapUsXfvzxx3rbf/nlF79PtRJcn0gVCiL69kfX8wYqc3MwVVX53Ia6bQmLT5/iyPw5HJk/p57IrlqwuNYCliuNWryFRCpFGRaGWq9HHhJC6dkzZO3aQeb2LRSf+ZHq0tKWNtFvaTDGat+1auDAgTz//POcOHGChIQEpFIpP/74I3//+9954AH3y/sEAm8ikUrR/iYehTaMvAPfYjYYUIT6vpFH3H0P1AsNNJS2ZY25WrtgeXsxy1mkcjmqiJqm46aqKorST1B44nhNPLZjJ9T6mIBJ3fIFDcZYk640iWjyAhIJ27dv96hRnkbEWP0XX9lcXVJMbtp+jFdGvjQ37upOExYrTfV99WYaVnPsdoSxohxjaRlQUxEX3LET6mg9Mgfz8ZqDvz3bHk23ClSEsPovvrTZbDCQf/Qw5Rd/aXYxQXMEqqFCA3uB9VbxgKeF1YrFYsFUUY6xrBwsFlTReoI7dEQdFe0RkfW3Z9tjwpqTk8O5c+dsDakBDAYDJ0+e5JFHHmm+pV5ECKv/4mubLWYzJRlnKDpxFEWYzu1fek8JVEP5ryt++cXtqQON4S1htcfaq8BYUQFYUEVFE9z+isi62ZHM355tjwjrP//5T5YtW4bRaKxV2iqRSEhMTOSzzz7znMVeQAir/9JSNldkZ5F34FskMinKMNfzsD0tUI4E9i95eRy80iDeSnNF1hfCao/FYsFUXo6xrCZcoIyMJLhdB9TReuTBwU5fx9+ebY+kW61Zs4aHH36YY8eOERkZyc6dO/nyyy/p3r07Y8aM8ZixAoGv0OhjaDV6LPKgYCqzs7GYzS1qj6PBhk9FRtI7pPYvb0unYrmKRCJBHhyMWq9HFR2NuaqKgqOHubz1azJ3bqM44yeqi4uuuRQupzzWhIQEvv76a9q2bctDDz3ExIkTSU5O5uDBg8yfP58tW7b4wla3ER6r/9LSNpuNRop+OFHTwzTc+ckE3vb8HBUY3LrhPw6PdcWL9bXH2himygqqS8vAbEYWFFTjycbEoAzT1Yt/t/RzUhePlLRGRkaSn59P27ZtiYuL44cffiA5OZmYmJhac6sEgkBDKpcTfkNvlOGRFHx/EKlK6Re18nVLZItPn6J3SAhHHOSOphw66FbTbF/0JmgMmVqDTF0TczUbDJRk/Ejx6R+QKOQExbZF0zrWpT92/oRTwpqcnMzcuXNZunQpI0aM4Nlnn6V79+7s3LmTDh06eNtGgcDrBLdthzIsjNzv0qjMyUEVGekXPUvtc2Dn1Gl6VHeRy1X8KawgVSptebIWk4mKy5cp+/lnkEhQRUWj7BVPtUSDPDjY6313PYFToQCj0ciaNWv4zW9+w5gxY/jLX/7Cv/71L3Q6HcuXL6dPnz6+sNVtRCjAf/E3m81GI0XpJyg523hooCU+UjuTptUUVrvrergt7b02hMVsxlRZgUZqpqS4AllQEEFt26OJiUERpvPaxN6mEHmsCGH1Z/zV5rKLFyj4/iAShcJhl6yWjlW626bQard9EYIVfxVXuGq3yWDAWFqKxWREIpGi1segiW2LKiLSpSyD5uIxYd2/fz+fffYZZ8+eRaFQ0LlzZ6ZPn06PHj08Zqy3EMLqv/izzdWlpeQfOoChoKAmNGC3oNLSwmrF1S5ajuxuiYbbruLIbovZjLG8HFNFBVgsNd5sbFvU+poFMG/GZj2SbvXRRx8xa9YsgoKCmDp1KrfddhsWi4WpU6eyaZPzP4iNGzcyatQohg4dyvvvv19r33/+8x+6detW6+v3v/89AOvWrau1/aGHHnL6ngKBuyhCQogeNhJtt3gqc3IwVba8kNbFUaMXV6k7OdYfGsA4g0QqRRESgjo6GrVej0yppOznc+Ts28Olr74kK3UXxRk/1XTk8nGHM6cCFGvXruWll17itttuq7W9f//+vPbaa4wf3/SKZGZmJs8//zyTJk0iPDycFStWMGDAABITEwEYPnw4H3zwAVAzyHDhwoVMnToVgGPHjtG1a1fmzZuHRCIhwgdTOQUCqMkaCOveE1VUNHkHD2Asq0AZEd7SZtXDfoqsfS8CZ+OvLdlw21NIlUqUV7xUi8WCqbKCovQTFFnMSGRy1NHRaFrHotTpkIeEenVx0ilhLS8v54Ybbqi3vV+/fhQUFDh1o/3791NdXc3MmTOJiYlhzZo1pKam2oRVr9ej1+sBePDBBxk9ejQTJkwAaoT1119/ZdasWcTHx/PGG284dU+BwFOoo/W0ShpL4bFjlF08T4jG/+a8abvF1/NYra0Krft7Pz7b4blWr9W6qOXLQYfeQCKRINcE2Ub0WEwmDEVFVGRmgsWCVKFA3ao16phWKHXhHs82cEpYp02bxooVK/jzn/9M2JVAfkVFBW+//TZ33nmnUzfKulKap9PpUKlUBAcHO8yBPX78OKmpqWzYsAEAs9lMcHAwEyZMYPTo0SxYsICFCxfavFuBwFfIVGoi+g9A3aoV1WfSMVRUo9T5j/da1zN11KowZXaNsDbkydp7rnA1RzZQBdaKRCZDERKC4kolm9lopConm/ILv4BEgkypJLzfADT6GI/cr0FhHTVqlE3BLRYLWVlZjBgxgrZt2yKVSrl48SIGg4Hu3bs7dSP7/gKN8dlnn9GrVy/bophUKuXTTz+17R8/fjyffPIJZrMZqZOufGNB5oaIjvZ9305PEIh2B5zN+p5U/6YDWd99R3lmJmofTYZ1FXvv9MSqVeSnp9teWz3ZiJ49SbBrovTJilds/79v3jx2pKUBNQKr1brXQMVTePz+EVefu/KsLMJUEOahZ7FBYf3DH/7gkRtYiY6OBqCwsBCVSkV5ebnto789qamptbzggoIC1q1bx8iRI4mPj8dkMiGXy50WVRBZAf5MINoMNXYre/TFoDlL7oljyNQqv6jYaoj299xPe2rE6cjbK22ebH56OimzZzv0YN+etxC4GhZoySwIb2dhVJVVkV9QhsHJZ9HtktZJkybV21ZaWsr58+cxmUx06NDBFhZwhsGDByOTyVi9ejURERGYTCYGDRrEvn37iIuLo1WrVmRmZpKVlUWvXr1s52k0GtauXcuuXbuYNm0amzZtng1N6QAAGmdJREFUYty4cU7fVyDwFhKplNDOXVBHR5N3+BCV2dn10rL8kYYmyh6ZP6fRxa5Aj7v6EqfcPoPBwEsvvcTgwYO5/fbbmTJlCsOGDWPu3LkYDAanbtS2bVuWL1/Onj17WLduHXPmzKF169ZMnz6dbdu2Adhirq1bt7adp1areeuttygqKuK5554jMTGRefPmufo+BQKvodCGETN8FGE9elKVlxsws6AcddRyNJfLftAhBG7WgC9xqkDgpZdeIiUlhYULF9KnTx/MZjOHDx9m2bJljBkzhrlz5/rCVrcRoQD/JRBthobtNhQWkn/4INVFRX7pvTb2kdqVyQaOaOmx3c2hKjcHXa/ehHSMc+p4j1ReDR48mLfeeouBAwfW2p6WlsbTTz/N3r17nTKmpRDC6r8Eos3QuN1mo5HSs2coOnkCmSaoRQYYNoSzAtWQyDpqvm2Ptyq3Ak1YnUq3slgshIfXTyvR6XSUl5c7ZYhAcL0glcvR/iYedUxr8o8cojI7C2VEZIs1DHEHR3FYqGm+TWRkvVisO20Lr2Wc+kkPHjyYV199lVdffZXQK399i4uLef311xk0aJBXDRQIAhVlWE3steTnsxSdOI5UqXBrDExLYi+ejha7gFox2msh59UTOCWs8+fP595772XkyJG0v9IT8pdffqFjx4688847XjVQIAhkJDIZ2s5d0ehbUXj8CBVZWSjDw5EFYPPmhrzYuk24xeKWk8KqUqn48ssvSUlJ4ezZs6jVauLi4hg6dGhANJ0VCFoaRWgoUYOHUX7pIgVHD2O0gDJc5xfNtF3FkRdrbcJ9z8mTtu3O9oi9FnFKWCdOnMi7777LmDFjxPBAgcBNJFIpwe3ao47WU3gynbLz51CEhiIPCmpp09ymMS+2qbzYaxmnPVZn81UFAkHjyNRqIvv2I6RDB/KvFBYoIyICanGrLjbxvLKIZfVcOXkSHAxBvNbjsE79JEeMGMGMGTMYOXIksbGxqFSqWvuffPJJrxgnEFzLqCKjiLlxDKU/n6Xo5AmkMjkKnS6gw2t1m7g0xLUeh3VKWH/88UcSEhLIz88nPz+/1r5AfggEgpZGKpej7fIbglrHUph+nPKLF1CEaW3t7gKNhrxQ+1CB1Zu9lkMFTgnrxx9/7G07BILrGnlwMFEDB1PRsROFx45cE+EBe2qJp13Oa920rWtFZBv9qW3cuJGtW7eiUCgYO3asU5MCBAKB+2j0MahuHEPp+XMUn6xp86cMDw/I7IGmqNuY+1oS2QaFdfXq1bz55psMGTIEo9HI3LlzOX36NE8//bQv7RMIrjukcjnazl0JbtOW4tM/UHruLLIgjV+3JXSHhooPwPnJB/5Kg8L6+eefs2zZMtucqy1btjBv3jyeeuopEVcVCHyATKMhvHdfgjt2ovD4sZrigrAwZGp1S5vmcZoS2aYmH/gbDQprZmYmQ4YMsb1OSkqioqKC7OxsYmI8M75AIBA0jVIXTvTwkVRkXqbw2FGqi7NRhof75dQCV3Crv0AD6Vue4MZhw/l8w2aPXKvBwI3RaERuFziXy+Uin1UgaCEkEglBrdvQasxN6Hr1prqkhMrcXJ+PdfYEdfu7XotcG0uOAsF1glQuJzSuM0GxsZRknKHkpx+RyGQodYFTHutOYYB928CGWhrajnUjXGBtG+gpGhXWL7/8kuDgYNtrs9nMV199RURERK3jbr/9do8ZJBAImkamUqPrkUBIh04UnT5F2S/nkKnUKLTaa34NxJlptC1Ng42uk5KSnLuARML27ds9apSnEY2u/ZdAtBn8z25DYSFFp05ScfkS8qDgBptre7thtLdw1m5rJkG985vwYn3W6HrHjh1O3UAgELQ8Sp2O6MFDqcrLo+jkCSqys1CEhCAPCm765GuIurmxVuzTt+yP9VaGgVOjWTzFxo0bef3116murmbGjBk8+OCDtn3/+c9/6g0JHDhwIB9//DEnT55k3rx5XLx4kbFjx7J48WLULqScCI/VfwlEm8G/7bZYLFTl5lB44jiGwgIU2lBbiey17rE6oqmYrBVdr94MeOuvTl3TI6NZPEFmZibPP/88kyZNIjw8nBUrVjBgwAASExMBGD58OB988AEAFy5cYOHChUydOhWAZ599Fq1Wy/z581m0aBFdunRh5syZvjJdIAgoJBIJ6mg9MaNGU5GVSdHJE1RmZaEI04JW09Lm+RxHXqmzYusuPhPW/fv3U11dzcyZM4mJiWHNmjWkpqbahFWv16PX6wF48MEHGT16NBMmTODSpUtkZGSwdOlSJk+ezLp160hJSRHCKhA0gUQqJah1GzQxrajIvEzRyROUZ2ZikiqRqa8/gbWnrtj6NCvAk2Rdmeyo0+lQqVQEBweTnZ1d77jjx4+TmprKhg0bap1nHWYYERFBRkaGj6wWCAIfiVRKUJtYNDGt0FSXcG5PGpXZ2Si0ode9wHoLnwmrNZTbVCrIZ599Rq9evejRo0et85pDY7GQhoiO9p+Rxa4QiHYHos0QqHbrSJjShrJLl8g9dozq4kKUYWEBMcVA68UwRkWlmojwYMI89DP1mbBGR0cDUFhYiEqlory83PbR357U1FTuvPNO22vrMYWFhbZ/HZ3XGGLxyn8JRJshsO3Oyy8HTTjBA0bYQgTVl7L9ekyMtxfdqkorkReUYXDyZ+o3i1eDBw9GJpOxevVqIiL+v717D4ryPPs4/t1lOe8KLIoiIAZlmdQDKgKiQU06Ri2ebUajtdWaGK1JNLFYo2mtU1IxVhMPYxKsTdSx1qxtTTQGOsakQIy12EaJRuUgqFRYPALiCrL3+weybzRaMV2WXbw+M8ywz/HHDnPxcD3P3reRhoYGEhMTOXDgAFFRUXTq1Iny8nIqKiro3bu3fb+IiAi6dOnCzp070Wq15Ofn8+KLLzorthBtlr1F0CkUq6WCq18fw2qpQOfnj07/4P/lif/ntM/AhYeHs3z5cnJzczGbzaSmphIaGsqMGTPYt28fgL3nGhoaetu+a9aswWq18tprrzFy5Eh+/OMfOyu2EG2eRqvFt1MoHYd+nw4DB+Ph64u1opz6qiqHtOIeRk59jrW1SCvAdbljZmjbuZVS1F26RFXBCaznz6P19sKzXUCrjkXQ4q0AZ33ySgjxcNJoNHgHB9MheBB1V69SXVRA7dkzaLTaxsFePDxaO6LLk8IqhLgnr4AAgvv1JyDmUWpKiqkuKgSl8AoMdPvxYFuSFFYhxH3p/P0J7NELQ/dorp07S/Wpk9hu3EBnMKDzlWdh7ySFVQjRbB7ePrTrFo0+8hGuV5RTdeI41y0V6Hx80RkMbX7IwuaSwiqEeGBanQ7/sHD8Oodx4+IFqgtPcb28HG3ToNsPeR9WCqsQ4jvTaDT4tO+AT/sO1FdXc630NNWni6DBhmdAAB7e3q0dsVVIYRVCOISnwUBgz94YomO4fr6MqlMnqbt6FZ2vDzr9w9UmkMIqhHAoD29v9F2j8O/SlRsXLlBVeIobFgsaDy2eAQFodW2/7LT9n1AI0So0Wi0+ISH4hIQ0tgnOllJTXIRquIlOr7cPvt0WSWEVQrQ4T4OBwO/1pF10DNfLz1NdVMB1SwUenl6NEyC2sZtdUliFEE6j9fTEP6ILfuER1F+9wrWSEmrOlIBSeBr0bWZ8WCmsQgin02g0eAUG4dUniHaPfq/xKrbwFFaLBa2XJ54G976KlcIqhGhVHt7e6CO74t8lkvorl6kpLeHamVKUzYanwT17sVJYhRAuQaPR4BVkxBhkJODRHlgryqkuLsRqqcCzTo9N6+02TxS4R0ohxEPFw9sb/y6R+EV0ob6qCq+aC5z78hi2+np0vr7o9HqXfi5WCqsQwmVpNBq8AgLo0D0cW8dIrJUWrpWcxmqpAI2m8YaXt09rx/wWKaxCCLeg1enwC+2MX2hnbtbWcv38f6g5XYTVYkHjqcPT0M5lWgWukUIIIR6Azs8PQ7fu6KO6UX/lMrXnzlFzpgRVX4+Hrw86f32rzngghVUI4baabnh5BRlp9+j3uHHxAtdKT3P9/HmUUnjq/fHw9XN6P1YKqxCiTdDqdPh27IRvx0403LBitTT2Y29UVoJWi6dej4ePc/qxTr1W/uCDDxgyZAgDBw7k97///bfWFxUVMWXKFPr27cvMmTOprm6c9OyLL74gJibG/jVq1ChnxhZCuBkPbx/8I7oQkjyE0GEjCOrZGwVYKy1YL1yg4caNFj2/065Yy8vLWbJkCePHjycoKIiVK1cSHx9PbGwsAA0NDTz77LOEhYWxdOlSli1bRkZGBgsWLODIkSMEBgayevVqNBoNfn7u98CwEKJ16Pz90Ud1a+zHVldxvaKcayWnuW6xoPXwwNNgcPw5HX7Ee/jiiy+or6/n2WefpWPHjmzcuJGcnBx7YT169ChlZWWkpaWRlJREbGwsQUFB9nU3btxg9uzZhIWF8bvf/c5ZsYUQbYinoR2ehnYYukVTX1XF9fL/cK3kNA1Wx17BOq2wVlRUABAYGIi3tzf+/v5YLBb7+rKyMgC2bNnC7NmziYmJYdWqVQQGBqLT6Rg6dCgTJkxgxYoVLFiwgMzMzGY3pP/b/N/30qGD4/+KOYM75nbHzCC5nc3huUPaQfdw1KB4bly+jM7HB52D/ht2WmFVSgHcsxg2rff19WXdunUsXryYtLQ0MjIyWLt2rX27kpISXnvtNSwWCx07dmzWuS9erMFmU83O2qGDgcrK6mZv7yrcMbc7ZgbJ7Wwtn9sTrjXAteadQ6vV/NcLNqfdvOrQoQMAV65coa6ujtraWkJCQuzrO3XqBMDo0aMZMmQICQkJFBYWUldXxx/+8AcOHToENPZiAXQu8iCwEELcyWnVacCAAXh4eJCRkYHRaKShoYHExEQOHDhAVFQUsbGxtG/fns2bN1NfX09eXh59+vTBy8sLs9mM2Wxm3rx5/OlPfyI+Pp7g4GBnRRdCiAfitCvW8PBwli9fTm5uLmazmdTUVEJDQ5kxYwb79u3Dy8uLjRs3UlNTw6JFizCZTLz66qsAvPHGG+j1ehYuXEhwcDDp6enOii2EEA9Mo5qam22Y9FhdlztmBsntbK6W22V6rEII8bCQwiqEEA72UNxa12offACG77KPK3DH3O6YGSS3s7lS7vtleSh6rEII4UzSChBCCAeTwiqEEA4mhVUIIRxMCqsQQjiYFFYhhHAwKaxCCOFgUliFEMLBpLAKIYSDSWEVQggHk8L6DfebRba1HDp0iJiYGD799FMAPv/8c4YPH05CQgLp6en22ReOHz/O2LFjiYuL4xe/+AVWqxWAc+fO2We/nT17NpcvX26xrEopVqxYwYABA0hMTCQtLQ2bzebSmQFsNhsrV65k0KBBPPbYY7z33nuAa7/Xd1qwYAExMTFukzs1NfW22Zc3bdrkFrmbRQmllFLnz59XPXr0UK+++qpatWqVMplM6ssvv2zVTDabTX388ccqPj5emUwmtX//fmW1WlV8fLyaM2eOevfdd5XJZFIfffSRUkqpkSNHqkmTJqmdO3eqHj16qIyMDKWUUjNnzlTDhw9XH374oUpISFBLly5tscx79uxRJpNJbdq0Sb399tvKZDKpv/71ry6dWSmlPvroI2UymdTmzZvVwoULlclkUl9//bXL574zv8lkcvnfkSZPPvmkevnll1Vubq76/PPPVVlZmVvkbg65Yr3lm7PIzp07F61WS05OTqtmOnnyJD//+c9JTk62Lzty5AhXr15l2rRpTJ8+nZCQEHJycigrK6OoqIiJEycyceJEevbsSXZ2NvX19Rw8eJARI0YwevRoHnvsMbKzs1ssc6dOnXjppZf46U9/ytSpU4HGqwpXzgwwbNgwsrKymDx5MhEREWg0Gqqqqlw+NzRO1Lls2TKio6MB1/8dAbh69SqlpaVkZ2fz3HPPYTab3eL3pLkeitGtmuN+s8i2htDQUPbt20dpaSl79uwBbs8JYDQasVgs9uVNU4YbjUaKioq4dOkS9fX139q+pcTFxREXFwfAhg0bAPDw8HDpzACenp507dqVtLQ0tm7dyo9+9COXf6+hsfWyePFiEhMTiY6OpqCgwC1yV1RUYDKZGDduHAaDgV/+8peYTCaXz91cUlhvUfeZRbY1BAQEEBAQQGlpqX3ZvXKq+wxS5uyfa/369WzatIlJkyYRFhZ21wyulhngqaeeIjIykrS0NHr27HnXHK6Ue9u2bRw7doxdu3axfft2AG7evHnXHK6U22Qy8eGHH9pf//GPf2T9+vV3zeFKuZtLWgG33G8WWVfRlLOpSX/lyhVCQkLsWa9cuXLb8qCgIDw8PL61vCVlZGSwbt06JkyYwK9//Wu3yFxUVMTf/vY3YmJimDZtGr6+vmzcuNHlc2dlZXH58mWGDBnC22+/DcCaNWtcPveJEydYv349tbW1QOMfg6aZl105d3PJFestd5tFdvDgwa0d61tiY2PR6/Vs2bKFgoICysvLSU5OJiIigi5durBz5060Wi35+fm8+OKLeHl5ER8fz969e+nWrRs5OTmMGDGixfL985//ZPXq1URFRTFq1CgOHjyI0Wh06cwAR48eZdGiRbzyyivU1NRw/fp1Zs+ezbJly1w699KlS6mpqQHAbDazc+dOVq1axaxZs1w6t9VqZd26dVRWVtKtWzdOnTrFr371K1avXu3SuZut9e6buZ5du3apoUOHqqSkJLVx48bWjmN38OBB+1MBSimVm5urnnzySRUfH69++9vfKpvNppRS6tixY2rs2LGqX79+KjU1VVmtVqWUUmfOnFFTpkxRffr0UbNmzVKXLl1qsazz5s2z351u+kpPT3fpzEo1PoGRnp6uEhISVFJSklq3bp1SyrXf6zutXbtWmUwmt8m9ZcsWNWjQINW3b1+Vnp6ubDabW+RuDplBQAghHEx6rEII4WBSWIUQwsGksAohhINJYRVCCAeTwiqEEA4mhVW0mkWLFt02utGdX3/5y1+IiYnhwIEDTsu0cOFC+yhidzKbzTzxxBMtev7s7GwWLFjQoucQLU8+ICBazZIlS+xFJC8vj/nz55Obm2tfbzAYSE5OJiAgwCl58vLyOHPmDI8//rhTznc3gwcP5p133uEf//gHiYmJrZZD/G+ksIpWYzAYMBgMAPbi2fTx1yY+Pj5Oy7Nhwwaefvppp53vXqZOncqGDRuksLoxaQUIl/bNVsATTzzBjh07mDhxIr1792bmzJmUlZXx/PPPExsby7hx4ygqKrLvm5eXxw9/+EN69+5NSkoKu3btuud5SktLOXToEEOHDrUvq6io4JlnnqFPnz5MmDCBc+fO3bbPp59+yvjx4+nVqxdxcXHMnz+fmpoa6urq6N+/P3v37rVva7PZSE5OJisri+rqaubPn09CQgL9+vXj+eefp7Ky0r7t4MGDOXz4MMXFxf/r2ydaiRRW4VbWrl3Lyy+/zLZt2/jqq68YP348ycnJmM1mtFotb775JgCVlZXMmjWL0aNHs3v3bubOnUtaWhr79++/63FzcnLo0aPHbW2HefPmYbPZMJvNPPPMM2zZssW+7uzZs7zwwgtMnjyZjz/+mDVr1nDw4EG2b9+Ol5cXw4YNIzMz07794cOHqa2tZejQoaxZs4aysjK2bt3K+++/z8WLF1m+fLl9W71eT69evVp9PGDx3UkrQLiVsWPHMmjQIAASEhK4fPkykyZNAmDMmDGYzWagcTi9xMREfvKTnwAQGRlJcXExmzdvvusNqK+++oqoqCj764KCAv7973/zySefEB4eTnR0NPn5+WRlZQHQ0NDAkiVL7OcODw9n4MCBFBYWAjB69GjmzJlDbW0tfn5+7N27l2HDhuHt7U1ZWRl+fn6Eh4fj7+/P66+/TlVV1W15unfvTn5+viPfOuFEUliFW4mIiLB/7+3tTefOnW97XVdXB0BxcTE5OTn07dvXvv7mzZsYjca7HvfixYv2EfgBCgsL0ev1hIeH25f17NnTXli7du2Kl5cXb731FgUFBRQUFFBYWEhKSgrQOFqawWDgs88+Y/jw4WRlZbFixQoApk+fzpw5c0hKSiIxMZFhw4Yxbty42/IEBgZSVlb2nd4j0fqksAq30jRmZxOt9u7drJs3b5KSksLPfvazZm2v1Wqx2Wy3LbtzfKJvnvvEiRM8/fTTPP7448TFxTF9+nQ2b9582/FGjhxJZmYmRqMRpRRJSUkAJCYmkp2dzf79+/n73/9Oeno6u3fvZuvWrfb9bTabSw7gLJpHCqtokx555BEOHz5MZGSkfdm2bduwWCy89NJL39o+ODjYPmAyNI5wf+3aNYqLi+0tguPHj9vXf/DBB/Tr14/Vq1fbl5WWlt52vlGjRjFjxgxCQkIYMWKEvTC/9957REdHM2bMGMaMGUNeXh5Tp07lwoULtG/fHmgc7PnOJySE+5CbV6JNmjJlCsePH2fVqlWUlJSQmZnJypUr6dix412379GjBydPnrS/7tatGwMGDGDx4sWcOHGCffv22ac+gcZ/1U+dOsWRI0coKSkhPT2d/Px86uvr7dvExsYSHBzM+++/b28RAJSXl/Ob3/yGf/3rX5w9e5bdu3fTuXNn+5xO0DiRZNPUMML9SGEVbVJYWBjvvPMOBw4cYNSoUaxYsYIXXniBKVOm3HX7wYMHc+rUKaqrq+3L3nzzTdq3b8/kyZN54403mDZtmn3dtGnT6NevHzNmzGDy5Mn2x76+eVUL8IMf/ACj0WifYBEanzbo378/c+fOJSUlheLiYt566y37pIs1NTWcPHmSIUOGOPItEU4kA10Lccv06dNJSUnhqaeectgxX3nlFYxGI6mpqc3e589//jN79uzh3XffdVgO4VxyxSrELc899xw7duxwyLGOHj3K9u3byczMZOLEiQ+0744dO5g1a5ZDcojWIYVViFuSkpKIjIzkk08++Z+PlZOTw+uvv87cuXNvez72fj777DPCwsLsTxAI9yStACGEcDC5YhVCCAeTwiqEEA4mhVUIIRxMCqsQQjiYFFYhhHAwKaxCCOFg/wfI63DckysxVQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(5,4))\n",
    "plt.xlabel(\"Time ({})\".format(time),  fontsize=fontsize)\n",
    "plt.ylabel(\"Probability of survival\",  fontsize=fontsize)\n",
    "\n",
    "\n",
    "plt.step(x = t_empirical_range, y = 1-km_dead_dist_sca, label='SCA', color='brown', lw=2)\n",
    "plt.fill_between( t_empirical_range, km_dist_ci_sca[:, 0], km_dist_ci_sca[:, 1], \n",
    "                 alpha=0.3, color='brown', linestyle=\"-\", linewidth=1)\n",
    "\n",
    "    \n",
    "plt.step(x = t_empirical_range, y = 1-empirical_dead, label='Empirical', color='black', lw=2)\n",
    "plt.legend(loc='best', fontsize=legend_size)\n",
    "plt.show()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calibration slope\n",
    "def print_cal(predicted, model):\n",
    "    slope, intercept, r_value, p_value, std_err = linregress(x=predicted, y=empirical_dead)\n",
    "    print( model, \"slope: \", slope, \"intercept: \", intercept, \"p_value: \", p_value, \"std_err: \", std_err)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SCA Calibration slope:  1.0223155101291994 intercept:  -0.013886760654121394 p_value:  0.0 std_err:  0.006076941106039274\n"
     ]
    }
   ],
   "source": [
    "print_cal(km_dead_dist_sca, \"SCA Calibration\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}