step5.1_make_and_test_classifiers.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "57c88209-c077-4a49-b748-d93db0bd14c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pickle\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "from scipy.stats import gaussian_kde\n",
    "import pcntoolkit as ptk\n",
    "from pcntoolkit.normative import estimate, predict, evaluate\n",
    "from pcntoolkit.util.utils import compute_MSLL, create_design_matrix\n",
    "from pcntoolkit.util.utils import calibration_descriptives\n",
    "from pcntoolkit.model.hbr import bspline_fit, bspline_transform\n",
    "from sklearn.preprocessing import OneHotEncoder\n",
    "from sklearn.svm import SVC\n",
    "from sklearn.model_selection import train_test_split\n",
    "from pcntoolkit.util.hbr_utils import *\n",
    "import sklearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "8b46a1f2-bcbd-40a0-b4ef-963c1545a326",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('/project_cephfs/3022017.02/projects/stijdboe/Data/10_folds_sexcov/fold_0/Z_train.pkl', 'rb') as file:\n",
    "    Z_all = pickle.load(file)\n",
    "all_sites = np.unique(Z_all['site_id'])\n",
    "\n",
    "with open('/project_cephfs/3022017.02/projects/stijdboe/Data/10_folds_sexcov/patient_levels.txt','rb') as file:\n",
    "    levels = list(pickle.load(file).to_numpy())\n",
    "levels.insert(0,'Healthy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e3ff8ef0-77f9-479a-9fd6-e2eba2409645",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract some patient data\n",
    "with open('/project_cephfs/3022017.02/projects/stijdboe/Data/10_folds_sexcov/X_patient.pkl', 'rb') as file:\n",
    "    X_patient = pickle.load(file)\n",
    "with open('/project_cephfs/3022017.02/projects/stijdboe/Data/10_folds_sexcov/Y_patient.pkl', 'rb') as file:\n",
    "    Y_patient = pickle.load(file)\n",
    "with open('/project_cephfs/3022017.02/projects/stijdboe/Data/10_folds_sexcov/Z_patient.pkl', 'rb') as file:\n",
    "    Z_patient = pickle.load(file)\n",
    "labels_patient = np.loadtxt('/project_cephfs/3022017.02/projects/stijdboe/Data/10_folds_sexcov/patient_labels.txt')[:,None]\n",
    "with open('/project_cephfs/3022017.02/projects/stijdboe/Data/10_folds_sexcov/fold_0/X_test.pkl', 'rb') as file:\n",
    "    X_control = pickle.load(file)\n",
    "with open('/project_cephfs/3022017.02/projects/stijdboe/Data/10_folds_sexcov/fold_0/Y_test.pkl', 'rb') as file:\n",
    "    Y_control = pickle.load(file)\n",
    "with open('/project_cephfs/3022017.02/projects/stijdboe/Data/10_folds_sexcov/fold_0/Z_test.pkl', 'rb') as file:\n",
    "    Z_control = pickle.load(file)\n",
    "labels_control = np.zeros_like(Z_control)\n",
    "with open('/project_cephfs/3022017.02/projects/stijdboe/Data/10_folds_sexcov/fold_0/X_train.pkl', 'rb') as file:\n",
    "    X_train = pickle.load(file)\n",
    "with open('/project_cephfs/3022017.02/projects/stijdboe/Data/10_folds_sexcov/fold_0/Y_train.pkl', 'rb') as file:\n",
    "    Y_train = pickle.load(file)\n",
    "with open('/project_cephfs/3022017.02/projects/stijdboe/Data/10_folds_sexcov/fold_0/Z_train.pkl', 'rb') as file:\n",
    "    Z_train = pickle.load(file)\n",
    "labels_train = np.zeros_like(Z_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "2bbbaa7d-0a56-4c4d-8be9-6c443a9056cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(70)\n",
    "idx = np.random.choice(X_control.shape[0], 500, replace=False)\n",
    "all_Xs = np.concatenate((X_patient.to_numpy(), X_control.to_numpy()[idx]), axis = 0)\n",
    "all_Zs = np.concatenate((Z_patient.to_numpy(), Z_control.to_numpy()[idx]), axis = 0)\n",
    "all_Ys = np.concatenate((Y_patient.to_numpy(), Y_control.to_numpy()[idx]), axis = 0)\n",
    "all_labels = np.concatenate((labels_patient, labels_control[idx]), axis = 0)\n",
    "all_features_baseline = np.concatenate(( all_Xs, all_Ys, all_Zs), axis = 1)\n",
    "# all_features_baseline = all_Ys.copy()\n",
    "\n",
    "all_features_baseline = (all_features_baseline-all_features_baseline.mean(axis=0, keepdims = True))/all_features_baseline.std(axis = 0, keepdims = True)\n",
    "X_base_train, X_base_test, Y_base_train, Y_base_test = train_test_split(all_features_baseline, all_labels, stratify = all_labels, random_state = 42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "7b190236-0e75-41fc-8371-cbab0370a619",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5721518987341773\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAF7CAYAAAAzPisLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABJvUlEQVR4nO2debgcVbW+3y8hgUAIiCKCIhFkjhAgzCCDOKCgICAgKFGviAOIXOHnLA4oiMoFHJCLMsmFCAIGZEYCARlDQsKMMsgoIDNkzvf7Y+9O+nR6qK4zp9f7PP2cU7v2qr27q3r1rlV7f0u2CYIgCJZshvR3B4IgCILeJ5x9EARBBxDOPgiCoAMIZx8EQdABhLMPgiDoAMLZB0EQdADh7IMgCDqAcPZBEAQdwFL93YEgGExIGgp8BBhN1ffH9i/7q09BUIRw9kHQHpcAs4AZwIJ+7ksQFCacfRC0xztsb9TfnQiCdomYfRC0x+WSPtDfnQiCdomRfRC0xy3ARZKGAHMBAbY9qn+7FQTNUaheBkFxJD0M7AHMcHx5gkFEhHGCoD0eAu4ORx8MNiKMEwTt8TQwSdLlwOxKYUy9DAY64eyDoD0eya/h+RUEg4Jw9kFQkLygam3bB/Z3X4KgXSJmHwQFsT0fWFlSjOiDQUeM7INBjaRtgWm2X5d0ILApcKLtxwrYrkEaqV8jaQSwlO1XW5g9CtwkaSLweqUwYvbBQCdG9sFg57fAG5I2Bo4CHgPOamUk6fPABcDvctE7gIsLtPcUcCnpu7N81SsIBjQxzz4Y1Ei60/amkr4HPGn795WyFnbTgC2AW21vkstm2H5PwXaXJy2meq2bbyEI+oQY2QeDnVclfRM4EPhrfog6rIDdbNtzKhuSlgJajnwkjZE0FbgbuEfSFEkblux7EPQZ4eyDwc6+pPnun7P9DPB24PgCdtdL+hYwQtL7gfNJipatOBU4wvYattcA/hv433JdD4K+I8I4QUeStW0+B3yApG9zJXBaq5Wxku6yvXGrsiDoSfId62G2Tyh9jHD2wWBE0o22t5P0Kl3DL4WEySTtCVxme3azenXsLgLuBM7ORQcC42zv0c5xgqBdJE2yvWNp+3D2QSci6XRgZ+AG4DzgStvzCti9CfgBsF0uugH4ge0Xe6uvQQAg6RhgBWACXaf93lnIPpx9MJiRtBbwhO3ZknYENgLOsv1SAdthwK6kuP92wNW2/6tB3bNtf0rSV22f2FP9D4KiSLquTrFt71zIPpx9MJjJUyjHkXLCXglMBNa1/eGC9sOADwGfAba3vXKDeveSfhgmAjuSwkULsf1CqTcQBH1ErKANBjsLbM/LMfj/sX1ynhrZFEkfAvYDdgImAacBn2hicgpwBbAmMIWuzt65PAh6FUkfATYElqmU2f5hEdtw9sFgZ66k/YGDgN1zWZF59uNJsfovFHlIa/sk4CRJv7X9xbKdDYKySDoFWJY0QDkN2Bu4rbB9hHGCwYykDYBDgJttnyvpXcC+to/t4XZG2X5F0kr19kcYJ+htJE23vVHV35HAhbYL5UQOZx8MerIK5Tp58wHbc5vULTVlU9KltneT9Ei26xLGsR1hnKBXkXSr7S0l3QJ8HPgPKWva2kXsI4wTDGryDJwzSWqUAlaXdJDtG+rVt71d/tuWeJnt3fLfd3Wju0HQHS6VtCJphfidpEHHaUWNY2QfDGokTQE+afuBvL0OcK7tzVrYnW37U63KGti+HViDqsFSox+XIOgpJC1deb4kaWnSQ9pZRRcGxsg+GOwMqzh6ANsP5umUregiXpaF0Jr+QOR6x5Hm5d8LzK80S1pcFQQNyZIHV9repeQhbiblayA7+NmS7qyUtSKcfTDYuUPS71kkX3AAaWpkXbJCZkUA7ZVKMTCHJHLWij1I8/jbklkIAtvzJb0haQXbLxe1k/Q2ksDfCEmbsOh50SjS7Jxix4kwTjCYybezXyatgBVphP2bVs5Y0k9tf7NEe5cD+4SOfVAGSX8CtgKupqvkwWFNbA4iTRUeB9xRtetV4AzbFxZqO5x90InkRVh/q4yw8oOvHW1f3MLuz8DGwLUkaWWg+Zc1CCpkx70Yts8sYLuX7T+XbjucfTCYyTloj2bxB6ZNp0JKmmZ7bE3Z1ErWqiZ2pb+sQdAd8oDke8B7c9H1wA+LhoQiZh8Mdn4PfI0Up5/fom419RL3tPw+2D4zJyd/Z/WD4SAogqS1gZ8CG9BV8qDIOo3fkzKkVWQ9PgWcTppz35LIVBUMdl62fbntZ23/p/IqYHeHpF9KWkvSmpJOoMmD3QqSdgemkXRykDRW0sRuvYOgkzgd+C0wjyR7cBaLJhe0Yi3b37f9cH79gDY0mSKMEwxKJFWmm30CGApcSNcYelONb0nLAd8FKtPgrgKOsf16Y6uF8/p3BiaVSVQe9D+Sjmi23/Yve7HtKbY3q75mJE22vX0B25uBI23fmLe3BX5ue+sibUcYJxis/KJme1zV/yY55IZkp/4NSSPbnFkzz/bLUleF4zbsg/6nsnp6XWBzkmw1JCG9huslJM2gybm2vVGBtmfllJgPSfoK8CTw1iKdBr4InClpBdLMsxdIs3QKESP7oCORtA1pqflI2++UtDFJAfNLLex+T5qJ8w1gL+Aw0sKuQ3q7z0HPIukqYC/br+bt5YHzbX+oQf018r9fzn+r13a8UURqWNLmwH3AisCPSHPlf2b71jb6PQrA9iut6naxC2cfDGYk/YT0ZXkpb78J+G/b32lhdytJInZiVTjmbttjWtgtC3yblKgcUsKUH9ue1a03EvQ5ku4HNq6RILjL9not7G6yvW2rsga2+9g+v1VZA9sVgU+TEvVUzzwrNO03HtAGg51dq1MQOuWCLZSlyvbjNUVNZ/Pk5e4TbX/b9ub59Z1w9IOWs4HbJB0t6fvAraQHpq1YTlIlB3HlLnG5gm3WW8hXdHHfZSRHP4M0maDyKkTE7IPBztAagagRwNIF7B7PX1JnieTDSLfXDSm73D0YmNg+Jq+Irjwc/YztllnOgM8Bf8ixcwMvA59tZiBpV9Ig5O2STqraNYo0M6cIy9hu+nC5GeHsg8HOH4FrJZ1O+uJ9liR53IpDgBNJmiNPkGbjfLmpRWIWMENS4eXuwYBmWeAV26dLWlnSu2w/0szA9hRg4xw7V8Ef/qdIUgcfpeto/FXSOpEinC3p88CldJ15VihxTsTsg0FPzie7C2mGwlW2r2xS9zjb/69onLSOfaygXULIoZtxJGG7dSStRnpA2zT2LmkV4CfAarZ3zdnStrb9+wJtDiNdp4WS7dTYfhk4BniJRbOCCifOCWcfDGry9LVzcqy+SP0ZJEnYW20XkoYNlkwkTQM2Ae6sekg/vdUUyhz6OR34tu2Nszz21CJrLSTtQHou8Cg52Q5wUJF8CJL+CWxp+/lWdesRYZxgsPM24Pas6/0Hkl54sxHMFcDzpIdsr5DTEVb+ukFawgp1tHgqdpGWcPAxx7YlGRYutCvCW2z/KctlY3uepKJSHb8EPuCaZDsUyKUA3AO8UbCdxQhnHwxqbH9H0ndJUyE/A/wqy8j+3vY/65h8x/aRkv5i+2MlmiyrxRMMPP4k6XfAijkW/lmKpfl7XdKbyaEUSVuRHtIWoWyyHUjX2zRJ11FCcTWcfTDoyaOzZ4BnSDMb3gRcIOlq20fVVK9k+2lrQUoVL9u+vHxvg4GC7Z9Lej/pWlgX+J7tqwuYHkFadbuWpJuAlYF9CjbbVrKdGi7Or1JEzD4Y1Eg6DDiIFJo5DbjY9tzKknTba9XUv5uUsPl7wJG1x3ODRBDd1eIJBh6Vh/WtyurYLU0aZa9LCuM9AAxxgexlKplspwiS/mx7r4b7w9kPfCTtBlxme0F/92WgIemHpJDNY3X2rW/7vpqy7UijqU+wSBOlgm3XnS+db50bYdtNtXiCgYekO2sf0hd8QFvPbrGyvkYt8jFEGGdwsB9wolKWpNNrHViHcyYpfIOkHYGNgLNsv1Tvc8qKgTdKuqPIVLkqu516prtBfyPpi8CXgDUlTa/atTxwUxO7bueCzQO3H7H4A/6mEwMK0nTkHiP7QUJewLE/6SGkSVO/zq2IOHUqefrcONIy8itJo/V1bbeUTJA0hsWTSDRdLp8fzH2fdBtu4EZStqAiGvrBACCvfH0TKYnIN6p2vdpsgZK65oK9nUXO/hXgzEYhwJpj/IOUbGRGi1ljbdPq7iKc/SBC0luAA4HDSUv73w2cZPvk/uxXf1K5wCUdCcyyfXKr29ls931gR5KzvwzYFbjR9t4t7K4mxVn/mIsOIOWu3aWxVTBQyXpHq9BVWOxfTeoPAfa3fU7J9q4D3tcbIdlW130IoQ0CJO0u6SLgb8AwYAvbu5ISX3+9XzvX/8yVtD/pIe2luazIVLa9gfcBz9j+DOmzLKKps5LtH9l+JL9+TJKrHZBI+mqRsk4kL8j7N3A18Nf8urSZTXbSX+hGs0cBl0n6pqQjKq+C/W11Lps+WA5nPzjYBzjB9ka2j7f9LIDtN2ghwNQBfAbYmpRl6hFJ72LRqLsZM/MXd14OkT1LsRRv10naT9KQ/PoEyUkMVOrJO4zv604MUA4nhfw2tP2e/CqSgORqSV+XtLqklSqvgm0eQ1oYtQzpGUHlVYSm59L2Vc2MI4wTDHqyauV6pBj6A7bnFLD5DfAt0sPv/wZeA6blUX4zu1dJcraV2/AhLBJE66kHbd0m3+18kvRsYXLVruWB+RF2WhhSeb/toqqTFbt6QmmFVlHniQHjWtWrsemRcxnOfhAg6ePAcaT0ZaJnn+APaiR9BDgF+Cfpc3kXKeNU4YVPkkYDo2xPb1V3sKCUVeld1HkICUxv18EtieTFTeuS7syq10z0Zg7aY4G/tRqF19j0yLkMZz8IyE/wd48pl4ujlG1oN9v/yNtrAX91i2xDue5HgffmzettX1KwzWq7SbabxnmD3iXry/wWWMX2GEkbAR/Nz1Oa2X2/XrntH7SwW5a0ivadtg+WtDYpHNTyOqi6M5wNzKUPB27h7AcBKpjyrBORdIPt91Zti+S439vErDLC2hyozKrYH7jDdtOsQQ3sptj+RmOr/qPMXaFS+ru18+aDHuCJWiRdT1oN/Tu3kWKyyn45pwT0RdubQJI4+HT+cRkB3Gx7bPu9X+zYG9q+p8G+bt3hh7MfwOSTC7ADSd3xYrrebrac17ukI+m3pAUqfyLF7PchLV+/CZrKH0wHxlamwOUpeFMLrJ4sZddftHNXmJ99nArsATxCciZrABcBhxR5FlKyj0OArWz/vaT97bY3r556KGlaK+craWuSsF27SefvsD2upr27bG9cpv81x244V767d/ixgnZgs3vV/2+wKMk1JMfW8c6eNKvh36QfRIDngJVIn12rz2hFoLKIZoU22ixr1x/8uw3n8B3StNXVK4v1JC0P/Br4bn71OLYXSPoFaVZVGZ7P4buKCuXewNMF7P4H+CBZNsP2XZKa3hFm5uTRfKW9tagahHUTNdnXzrlcjHD2A5jKzBBJ29rusoxbSVe9o8mj6um2Tyhh/lNgap6RIVIM/lu9aNdf3JHDDhfT+q7w46Q1HG9U1XtV0peAW+glZ5+5StJewIUlVpZ+mXRHsp6kJ0l3JQcUMbT9eIr8LaSIbPXRpLwIq0s6B9iWnpvO2uy9t3MuFyPCOIOAerd2rZZGdwqSrnNJ3RpJq5Li7yJlrnqmN+36A6XcvLXYdQTf1EQETNIMF8jEVJaqB5fzgZm0EY+WtJntKUrJR4bkH6jdWz1wl3QBKZnIr4CtSEnnx9ner0Cbb842Am5xyexRdY7bLIxT+FzWI0b2A5gcU9wGWLlmld0oksxuAH+X9CtgAl0TgDeVHJZ0LfAL2xOryk61fXBv2PUXrdYN1FaX9CbqhxJ6VXHVdtGFRfX4X0kH2Z4BIGk/UoKZVrOrSiWdlzSRlF1qYjsPdgvS8LlIm+dyMWJkP4BRyle5I+miPKVq16vAJbYf6oU2lyLpxFSmLt4HXDFQ52WrvvSw3UJyWNLDwOOkOc8/yGUt75bK2vUX7UxLlPQoyanXc/aFFg11o58ihV7eZftHklYHVrV9WwHbNYELsv12wKdJ03F7ZRZR/l7uC3wEuI000LjU9qwWdt36bpWdYrrQPpz9wEfSGq6j194L7awGXEd6uDWV9KXfhDQTaCfbT/V2H/oKpZy1WwAnkZI+HwhcV8DZl7LrL7o7LbGvyLOqFgA7214/32FcZXvzgvbrkGLZjwN72J5ZwOZdwKEkxdRqIbSPFmxzKLAz8HngQy2ms3b7u9XdcxlhnAGMpEtY9MR/sf1FL8o2+AnwW9v/U9OPw0gPJutpc/QrklYh9Xs127tK2gDY2q216pVHVF+SNJ4kVfymIk2WtOsvlrV9W831U3ckqUXZuOrSKjTWTbZ0Ui+dmtt6MU8FbYikGXR9oLkSKbx5qyQKTIe9mDT18hLaDFPl2Ti7k0b4m5LyKjSjJ75bhc9lPcLZD2x+3sftbWV7fG2h7ZMkPVCn/kDgDJK2/7fz9oOk2+pWzn5hWMz2GdlxtIzXdsOuv2hnWuIvqv7fjK65UU0axTakm4ux5uaRcqWfK9PaAe/WxvHrMcv2Se0a5RkxW5Jm5PyatIq6VV974rtVdopp6neEcZZc2l1xpyZ62M329SftLqiRNMr2K2qgUugGySvK2vU3OZ59KulB/4ukaYkH2n60hV3h860eWIwl6QC6jpL3Br5j+/yCfdgY2D5vTrZ9VwGbT5J+nK6ijXzCkj4EXG27yDTNik23v1tlz2WFGNkPApS0N37K4lmVWj0w+xntrbhbQYtW7XbpAmkG0EDk9TwNrjLa2QpoNqL8P9KIcEq2qb4nNo1ljsva9Su2HwZ2qZ6WWNS0jWa6vRjL9jmSppByDIgUdy903Sppun+eRQvo/phnSLVK6vMe4FOkO5bKyLzlHYztKyRtoySgVx3rb5blrNvfrW6cy9RQjOwHPpJuJKXCO4EUJ/wM6dzVFXKqsmtLU6fBPN6FdHfqV2+Q48wnA2OAu4GVgb29BChYqkVSCzdRZ5R0oO0/NjpGM9tsX3iGkaS7qVmMlctHkuagF9WoaStrVJXddNJzmtfz9nIkrZpW0hf3AxsVufOosTsbWAuYxqJFWLZ9WBOb0t+t7p7LCjGyHxyMsH2tJOVZOUdLmkz6AWhGWyvuBqIzb4XtO/NUuHVJo6QHbM8tYivp7SxK/Fw53g29ZVeC7sw9X67JMeqO8CSdXLXvHZK6xLObOLMFtY4+139NUqHRpKRDSdfzv0kOVLkvRTSHRNeVrxX7VtxFkr54tkgfqxgHbOA2Rsrd/G61fS7rEc5+cDBLSSzqIaVUak+S4vCtGEWbmjqS1gUOputc4FNtP1im433EFiyaPrdpnonRKnH4caQY8b1Ujc5I+WV73K4MbiG128L2d/nfa1xcauOOqv+nNKhTtzl1fzHWV0kywWUSt/+BNAPnory9B60f0EO6i7hf0u10HQy1muV2N2nKZOGHo1D+u1XyXC7efoRxBj6SNiddGCsCPyI58eNt39LD7WxN+iE4FbiTRXOBPw98vKfb6wnK3FJnuwdIt/BtCViVtSuDpJ8BD9s+pab8a8DbbDfNOZrr9rrUhnpgMZbKZ40aQpItmEVaUCXgBttTC9juUK/c9vUF+jqWtKCq0I9ET3y3unsuw9kPIlRQd1vSUbZ/VnNbvpBGjlDS5cBxtifVlO8AfMMpyXmPI+kg8sguF90HnNRqdJ5t76PNW+psdzmwj+3Xetuu7PuTdC8wpnZaX3Zw05vFwrVIauNw0rOeCqOAPd1Ajrc756IMVXHoDSmZNUrSzbZLKWYqZYFa2/Y1SklJhrZ68FnmR6I7362y57KWCOMMAlSluw0U0d2uzGK4o8H+RqxVezFCuoglndqkfz8lJURfLPYp6bhmI1BJnyZdxEewaMSzKXB8kXAMbd5SV/0AvgFMU9K6qXYujX4Iy9p15/251tHnwgVSnVV2XRlOul6Womus9xXStMYe62t2mC85z6uXtBMplPIo8OsWD0ArfftXfg3PLygejy6lmCnp86Swykqku8O3k9ZRvK+ZXf4+rEISwwO4rd61X0Op71am7XNZjxjZDwIk3Uo6qRPdi0veJU2xvVmDfc3U+F4GXgK+ZPuvRe3y/luA/VwzVzhPazvP9lYN7Cqri5enjVvqPHJtiO26KyG7YVfq/eU6twOfdI0GktJU3HNdIHG12pDa6Ma5uJU0wnxK0ljgGtJU4Y2Aubb/q0Db+7hmTn29sga2FcXMeaRwTiHFTEnTSM97bq36XrVU95T0CeB4YFJua3vgSNsXNLEp9d2qqVf4XNYjRvaDBJfQ3VZahfj/WHx+fqN5xKvXzsCoHIo06mnEw6T5yudI2hX4uheJQrUagY6qdS65j49KavZlLbW6uOKU8/S8Wc4LY5Sm/S3d03aUf38A3wMul/RjFj0wHQd8kzQCL8Ibko4nhUlaXQNl+zrCi7RdDgT+YPsXOdw0rWA/vwnUOvZ6ZYvh8oqZs23PqXyvlITKiox+vw1sXhnN5+/ZNSQxtkaU/W5V0865XIxw9oODxyVtQ5r1MJyku11kwck5JOmAj5CUMw8iZXJqxJFN9jULCdn23ZK2JC3kukPSJ11srnszwaqG+yrx0XphojxjpulDNuBaYBegEnsfQVpJuU0P25V6fwC2L5e0B+m8HJqL7wH2cpbzLUDlGtiN1tdA2b5W/6DvTHLShcJNeXDwYeDtNc5wFAV1XyRda/t9rcrqcL2kbwEjJL0f+BKtZZEhLWiqDtv8BxjSwqbsd6uads7lYkQYZxAg6S0k3e1dSF+sq4CvusU0tcqto6qSUki63nbdB0zd6F+X29D8BT6FpAx5oJssBZf0BvCPeruANW0vV2dfw7ZzWcMkHFV1FpNUqFfWXbvuvr/u0s41ULavkk4EViU9N/kosI7tuUpJXia6iXJlfv40liTr8WPSyHo+ab79JNsvNrFdBliWpCa5I4t+dEYBl9tev5Ftth8CfI40NVnAlcBpreL+eXS9EUnTHtJU3OnNnk31BN39PsfIfhDglAWnUJq1GiqLi56W9BHgKeAdjSorJWVo1o9GU8tUU+9ySZuR5j+3WhTT9AvZCElfJI3E1lJaQVlheaBI4urXJW3qrIMiaRwtRtol7eq9P5HOQ9N0ht04H9W0cw2U7evhJIe3KrCdFy1qW5v08LMZ95Ku7eHAZ3N7q5PE7S5tYfuF3PZqdF0X8CpJqqEp+c7jYuBi2y1HyJLeTdKSP1JJ+qAy1fNm0qi7mW1fn8vF+xAj+4GLGkydrNBoBkiV/W7AZNKX52TSiOcHrsqyVFP/OZIe+LnArSzuxOuGRpSSND/onJ4vz+rYC3gMmOCahSCtyHcy/2k2wpK0Akla+KfAsaRcsAA3utgc63GkW+KnSJ/xasC+tpsuJiprl23HAp8EPkESsfqz7V81qV/qfNQco61roGxfm9hd6CYaNZJOIM00OcKLdHVGkZ7JzLT91Sa2m5MyTO1t++T8EH0v0iygo91Y1E6k1bpfgYUCgfOBk23/sEl7lwLfqg1P5mvi+7Z3b2Lbb+eyupF4DdAXKSZXeT1as31QE7tlSCOeX5FGP0sVbG8o8CGS6uBU0m31hgXs7gRWyv+/l+QI9yItALughe1WpFkNF5IWmdwNPENawv6hAm1/FZgB/AD4ITAdOLSA3T75yzKGJNT1V2DTnrYD1iE9aL2PpH1/KPBYb56PstdA2b528z0+RB501nnvD/XGdUdKWXg1KStWpWxNUhjna03s7m6yb8ZAOpd1j1PGKF59/wKmtlF3AvDHfGFcDJxYor2lgfGkB0BNnSdwV9X/vyaNqirb01rY3kGKme5Dkm3dKpevV+Q9Z+e+XNX2cqT4aUu7/Hc7ktTBx0hT8HrUjrSy9Hrg3VVlD/fm+Sh7DZTta3feI+mOsO193bnusrN9S53ylZtdc8A/yuzrj3NZ79XqCXIwcGgn3raB7QOdNDX2ZpHOd0skLZ3jkX8kJeU4iSZaOpmhedoapAUpf6va1+q50FK2r3KaT/2M87Jx2/cX7TLlRLAqNh8BTrH9FxYt5ulJu71IdyrXSfpfSRUJ30KUPB9Q7hoo29fuvMd7c9ivC5IOBFpdA2Wvu2FOz8G64BS3H9bE7nalhVi1ff0cBbSE+vhcLkY8oF0yWaj6aHtei9lvC5F0Jik8cTkpFnh3wfbOJU1je570sHJyPt67aa4tD12FsmofdBb5gTudciJYT0r6HWmG03GSlqb19Lm27WxfBFyU5+fvQQohrKKUc/Ui21c1su3G+YAS10DZvnbnPZKc3oWSPsuiXAGbk6a07tmiy2Wvu2YrepvtO5z0Pg+g67qH4a362tfnsm4f8m1CMABRWhlYOUHLkpbqQ4sVgpLmA69X1R2RbVvZLaiyq74wWq5IVEoasiopSXRFV3wdYKSbZP6p6mt1PyttLmO72UircoxNaV8Ea1lSDHWG7YfyNMH3tHBMpe1qjrESKWy1r5ssiOnm+Sh1DZTta3ftJO1MWiwk4B7b1xZsp+3rruaz6bKLAteckhxEZfX6Pbb/1qx+tun/cxnOPgiCYMknYvZBEAQdQDj7IAiCDiCc/SBD0sGDwa4/2oz3OLDa7IS+Dqb3GM5+8FH24upru/5oM97jwGqzE/o6aN5jOPsgCIIOIGbjDDCWWXEZj1x1ZMP9s16axTIrLrNY+Zz7m+d1nstshjWQXdfQxr/5cxbMYviQxdsDYEHza2cOsxneVOq9fbtm12uz99iMsnZNbVvMhZ7rWQxTg8+17Hss22YLH9D02hnWeKnOnAUzGT5kRN19nttcvbivz+WAuna6YfcqLz5ve+V6+2JR1QBj5Koj+dhZu7Vt96+tZ7Wu1IChI8up7HpOs/UnvcOCOXNbV6pr2DLXS4+iYUUW49bH88v1VUOH9ml7AEutXNevtGTe08+UbrM0Q8p9PqXp42sO4Bpf0DCTVYRxgiAIOoBw9kEQBB3AEu/sJb1Wsz1eUktd7gbH2jFrWlf+36Zq3xmSCmd6D4Ig6EuWeGffi+xI63ylQRAEA4KOdvaSVpb0Z0m359e2uXwLSX+XNDX/XbfGbjQp4e/XJE2TVJEcfW+u/3BllC/pbEkfq7I9R1KRFGRBEAQ9RifMxhkhaVrV9kpAJY3XicAJtm+U9E5Sppr1STra781yorsAPyFpdgNg+1FJpwCv2f45LNS0XpWkvrhebuMC4DSS5OtflFLpbUPKNLWQvCLuYIDl3tar+aeDIOhQOsHZz7Q9trIhaTxJgxqSJvkGVfrQoyQtD6wAnClpbZIcaUuZ3czFtheQEjKsAim3pKRfS3or8HFSLs8uk4xtnwqcCvCW9d8SCx+CIOhxOsHZN2MIsLXtLkkzlBJ9X2d7zxyymVTweLOrD1P1/9nAAcB+wGdL9zYIgqAkHR2zB64iZZgHQNLY/O8KwJP5//ENbF8Fli/YzhmkLDfYvqe9LgZBEHSfTnf2hwHjJE2XdC/poSvAz4CfSrqJlBW+HpcAe9Y8oK2L7X8D95FS6AVBEPQ5S3wYx/bImu0zSCNtctLhfevY3AysU1X03Vw+iRzSsf0gsFFVncmN2s2p7NYm5cwMgiDoczp9ZN/r5Nk89wMn226VfDsIgqBXWOJH9v2N7WuAdxatP+f+Bfxry3q5kJtz5VPT2rap8MHVxpa2DerjuX0vEud+EN7qF0GzsvTD5zOQiJF9EARBBxDOPgiCoAMIZx8EQdABhLMviKT5eZrl3ZLOzzNsGtX9qKRv9GX/giAImhHOvjgzbY+1PQaYw6I5+Ythe6LtY/uua0EQBM0JZ1+OycC7Je0u6dasjnlNRQ+nWjM/69yfVKuGGQRB0JeEs28TSUsBuwIzgBuBrWxvApwHHNXArKKGuRuw2Ihf0sGS7pB0x9wu8jpBEAQ9Q8yzL061VPJk4PfAusAESasCw4FHGtgupoZZTbXq5SitFKqXQRD0OOHsi9NFKhkWqmP+0vZESTsCRzewbaSGGQRB0CdEGKd7VKtjHtSsYhAEQX8Szr57HA2cL2ky8Hw/9yUIgqAhEcYpSK16Zi77C/CXOuVnsEhZc3yr4wRBEPQ2MbIPgiDoAMLZB0EQdAARxllC6I5M8UNnblrKbr1jysnz+19Ptq7UgAVz5pY0LClvO6RRorLmaGg5OwDPL9fXsm1qaPkxX9nzMXTkcqXs5r/ySik7oPS5LHvtLLXG6uXaAzyr5HqbJorTMbIPgiDoAMLZB0EQdADh7BsgyZLOrtpeStJzki7N2+Pz9lRJD0m6UtI2VfXPkPSIpLskPSjpLElv74/3EgRBEM6+Ma8DYySNyNvvZ9ECqgoTbG9ie22S5s2Fktav2n+k7Y1JsgpTgeskDe/tjgdBENQSzr45lwMfyf/vD5zbqKLt60j6NgfX2WfbJ5Aen+zaC/0MgiBoSjj75pwH7CdpGWAj4NYW9e8E1uvG/iAIgl4hnH0TbE8HRpNG9ZcVMGklclZ3f0gcB0HQ24Szb81E4Oc0CeFUsQlwX7v7bZ9qe5ztccNYulwvgyAImhCLqlrzB+Bl2zOyjHFdJO1AitfvVGefgENJSUyu6J1uBkEQNCacfQtsPwGc2GD3vpK2A5YlJS7Zy3b1yP14Sd/N+28BdrI9p1c7HARBUIdw9g1ooHI5CZiU/z+DrGzZwH58r3QsCIKgBBGzD4Ig6ADC2QdBEHQAEcYZiJRR5yur6gi887xyaoAvj125lN3yyy1Tyg5g6Ouzyhm+8FIpMy1V7isy79luJC4reS7LZqr33PKPkYa+5c3l2pxZ8jx2h258R8ow77HH+7S9VsTIPgiCoAMIZx8EQdABhLMPgiDoAMLZ9xCSvi3pHknTJU2TtEP+W/16RdJx/d3XIAg6j3hA2wNI2hrYDdjU9mxJbwGG2x5bVec9JH2d/+mXTgZB0NGEs+8ZVgWetz0bwHaXqRhZNfMc4Mu2n+6H/gVB0OFEGKdnuApYPWek+k3WyanmZ8BNtifWMw7VyyAIeptw9j2A7deAzUhCaM8BEySNB5C0K7AL8N9N7EP1MgiCXiXCOD2E7fkk3ZxJkmYAB0n6K/A74GO23+jP/gVB0NnEyL4HkLSupLWrisYCj5HkkU+2PbVfOhYEQZCJkX3PMBI4WdKKwDzgH8AJwA3AOyUdUFX3attH9n0XgyDoZMLZ9wC2pwDb1NnVKk1hEARBnxBhnCAIgg4gRvYDDA0ZwpAR7atCLnj99dJtLjf1X6XsPHLZUnb/ft/bStkBrPiPcoqZw5cZVsrO/yypXNjHCovdaXOp1d9Rusl5jz9R2jboW2JkHwRB0AGEsw+CIOgAwtkHQRB0AOHsu4mk+VnR8i5Jd0raJpePljRT0lRJ90m6TdJB/d3fIAg6k3hA231mVtQtJX0Q+ClQ0cb5p+1N8r41gQslDbF9er/0NAiCjiVG9j3LKODFejtsPwwcARzWpz0KgiAgRvY9wQhJ04BlSFLHOzepeyewXm2hpINJImoso+V6oYtBEHQ6MbLvPjNtj7W9HvAh4CxJjVbO1i2vVr0crnLzyIMgCJoRzr4HsX0z8BZg5QZVNgHu67seBUEQJMLZ9yCS1gOGAv+ps2808HPg5D7uVhAEQcTse4BKzB5SmOYg2/NzJGctSVNJ8fxXSXLHMRMnCII+J5x9N7E9tEH5o8CIvu1NEARBfSKMEwRB0AGEsw+CIOgAIowzwPCCBd2SKy7DvGf+XcpuwQ6blLIb/qpL2QE8vU25hOzv+uPzpezm9fG56A9CprgziJF9EARBBxDOPgiCoAPodWcvaUVJX6ra3lHSpb3dblkkTZI0rgeO8/ee6E8QBEFP0Bcj+xWBL7WqVBRJg+I5g+16CciDIAj6hR539pKOkHR3fh0OHEtaXDRN0vG52khJF0i6X9I5FS0ZSZtJul7SFElXSlo1l0+S9BNJ1wNflXSGpFMkTZb0oKTdcr0Ns278NEnTJa0t6UeSvlrVv2MkHZb/P0rSjKxFf2zV29gnH+dBSdvnustIOj3Xnyppp0Zt5vLX8t9VJd2Q999dOV4QBEFf0qOjZEmbAZ8BtiStJr0VOBAYU6X5viNJI2ZD4CngJmBbSbeSpAQ+Zvs5SfsCxwCfzYdf0fYO+RhnAKNJuvFrAddJejdwCHCi7XMkDSdJF/weuBA4UdIQYD9gC0m7AnsAW9p+Q9JK1Z+L7S0kfRj4PrAL8GUA2+/JsghXSVqnQZvVfBK40vYxkoYC5bJ0B0EQdIOeDolsB1xk+3UASRcC9Uayt9l+IteZRnLcLwFjgKvzQH8o8HSVzYSaY/zJ9gLgIUkPk6SDbwa+LekdwIW2HwIelfQfSZsAqwBTbf9H0i7A6bbfALD9QtWxL8x/p+S+Vd7bybnu/ZIeA9Zp0GY1twN/kDQMuNj2tNoPo4vEcfwWBEHQC/R0GKeRtG8ts6v+n0/60RFwT5YLHmv7PbY/UFWvdsJz7WRt2/4/4KPATOBKSRVt+dOA8aS7jj9U9bXRhO9K/yp9q9RfjCZtVvbfALwXeBI4W9Kn6xxjocTxMMrNIw+CIGhGTzv7G4A9JC0raTlgT1KYZvkCtg8AK0vaGkDSMEkbNqm/j6QhktYC1gQeyKn/HrZ9EjAR2CjXvYikNb85cGUuuwr4rKRlc3sr0ZwbgANy3XWAd7Zok1x3DeBZ2/9LCilt2qKdIAiCHqdHwzi278zx9Nty0Wm2p0i6SdLdwOXAXxvYzpG0N3CSpBVy3/4HuKdBcw8A15NCM4fYnpXj/AdKmgs8A/yw6tjXAS/Znp/LrpA0FrhD0hzgMuBbTd7eb4BTJM0A5gHjbc9u1GYVOwJH5v2vAYuN7IMgCHob2eWXrvcX+QflUtsXFKw/hJQScJ86MfUBxSit5C31vv7uRiHKyiW89vbyoaqX1yp3M/quPz5Zym7eI4+VsguC/uAaXzDFdt11Qkv8ClpJGwD/AK4d6I4+CIKgtxgUC5RqsT2+jbr3kmL6QRAEHcugdPbBwGDY82+UsnvlsDml21x9r0aPcJrjDdctZTe0pN38ex4oZRcEvcUSH8YJgiAIwtkHQRB0BP3i7CWNl/SrBvsuk7RiD7a1vaR7sjbNCEnH5+3jW1svdqxv1WyHsmUQBIOCATeyt/1h2y/14CEPAH6eV+XOBL4AbGr7yBLH6uLsQ9kyCILBQo86e0nLSfprVpG8W9K+kjaX9Pdcdpukymra1SRdIekhST+rOsajkt4i6ZA8Gp8m6ZG8KApJ+2flybslHVdl91tJd+RR+w9y2X8BnwC+p6SuORFYDrg1921lSX+WdHt+bZvtRmqRwuV0SXspqWKOyP05J9erKFtOyKJplb6ckW2G5juJ2/NxvtCTn3cQBEFReno2zoeAp2x/BCCvhJ0K7Gv7dkmjSBoyAGNJ6pezSbIDJ9t+vHIg26eQVqwOA/4G/FLSasBxwGbAiyTlyT1sXwx82/YLWVnyWkkb2T5N0nZULcCS9FqVAuf/ASfYvlHSO0lSCusD3wVetv2eXO9Ntv8s6SsV2xrOA/YFLlNSvnwf8EXgc/k4m0taGrhJ0lW2H+nWpxwEQdAmPe3sZwA/zyPuS0lKlk/bvh3A9isASqqW19p+OW/fC6wBPF7nmCcCf7N9iaSPAZNsP5ftziGJjF0MfCKrRy4FrApsAExv0d9dgA1yfwBG5TuPXUhSyOR+v9jiOJeTZB6WJv3g3WB7pqQPABtlGQiAFYC1gS7OXqF6GQRBL9PT2jgPKmnafxj4KUlsrJWyJHRVl1yIpPGkH4GvVIrqHUjSu4CvA5vbfjHLKSxToMtDgK1zLL/6eM0UMRcj6/JMAj5IGuGfW9XfQ21f2cg2258KnApJLqFou0EQBEXp6Zj9asAbtv8I/BzYihSb3zzvX14F0wrmH42vAwdm3XpIyVB2yDH9ocD+JDG0USQJ5JclrQLsWrDLV7Hoh4QsjFav/E3537k5rFSP80gSytuzSFnzSuCLFRtJ6yipgQZBEPQpPR3GeQ9wvKQFwFxS3FrAyZJGkOL1uxQ81leAlUhZqADusP1fkr4JXJePe5ntvwBImkpSyHyYJKtchMOAX0uaTvosbiBlnvpxLr+bdNfxA1JCk1OB6ZLutH1AzbGuAs4CJtquLBE9jZT85M58t/AcKTtWEARBnzIoVS+XZAaT6mVZKYHHflh+jPGOknIJZftalpBLCPqDjla9DIIgCMLZB0EQdAShehmUZv7y5ZKQzPpXN6aXbrVR6zp1mF+2vVtazd4NgsFBjOyDIAg6gHD2QRAEHUA4+yAIgg6gx5y9+la2eJKkutOLeuDYC/sq6TBJ91WEz4IgCAYrffKA1vaHW9fqWSQtZXteu3Y1ff0SsGtR4bKybQZBEPQ2LUf26l/Z4tck/ULSnZKulbRyVdf2yW0/KGn7XH+8pPMlXUJSxFxO0h+yxPDULKRWqXdhi76eQkpUPlHS1yStJOliJaniWyRtlOsfLelUSVcBZ+XtMyVdlY/1cUk/y+/viiZyC0EQBL1GkTBORbZ4Y9tjgCuACcBXbW9Mkj+oli3elySbsK+k1asPZPuULBG8OfAEXWWLd872m0vaI5ssB9xpe1OSBs73qw63lO0tgMNryrcGDrK9M/BtkmLm5sBOJCmHijZNq74eAjwF7GT7BJJkwlTbG5GSmJxVVX0z4GO2P5m31wI+AnwM+CNwXZZLnpnLuyDpYCUt/jvmdtGHC4Ig6BmKOPsZwC6Sjssj6HdSI1tcFbq41vbLtmcBFdnieiyULSY5/km2n8vHqcgWAywg/bBAcprbVR3jwvx3Ckl/psLVtl/I/38A+IakacAkkhLmO9vsa4XtgLPze/4b8GYlvX5IejjVypmX255L+uyGkn4gydvVfSUf71Tb42yPG0a5uetBEATNaBmz7y/Z4kbdqdNWbTuvVzcH7GW7i1CJpC2L9LWGev2s9Of1mvLZALYXSJrrRQJECwq0EwRB0OMUidn3l2xxpX+VxB+fBG4s9rYWciVwaFacRNImbdpXcwMpny2SdgSeryRjCYIgGOgUcdL9JltMGjFvKGkK8DIpxt4OPwL+hyRLLOBRYLc2j1HhaOD0LIf8BnBQyeMEQRD0OQNa4lgpX+zI/u5HXzKYJI7L6tT8Y9/y2jjvnvBGadtShDZOMIgIieMgCIIOZ0A/LOy0UT0AEho2vG0zz53TulKjJku0B+CSo95331LKbFDx6IRydz0Aq51ebkbW8CtuL91mJzB0/bVL2S34x2Ol7LrznewNYmQfBEHQAYSzD4Ig6ADC2QdBEHQA4ex7AEl7SrKk9fL2aEkzsx7PfVnDJ6ZqBkHQb4Sz7xn2Jy342q+q7J+2N7G9fi7/mqTP9EvvgiDoeMLZdxNJI4Ftgc/R1dkvxPbDwBHAYX3YtSAIgoWEs+8+ewBX2H4QeEHSpg3q3Qms12e9CoIgqCKcfffZHzgv/39e3q5HQ8G3LhLHntXT/QuCIBjYi6oGOpLeTNLhHyPJJDljA7+pU30T4L56x7F9KnAqwKghbx64+hVBEAxaYmTfPfYGzrK9hu3RtlcHHgHeUV1J0miSYujJfd/FIAiCGNl3l/2BY2vK/kzKZLWWpKmkhCmvAifbPr2P+xcEQQCEs+8WtnesU3YScFLf9yYIgqAxEcYJgiDoAMLZB0EQdAARxhlo2H0ujTrQpFiXBOa8VD5x/IvrDCtlt8oVret0Mpo5u3WlOiwp348Y2QdBEHQA4eyDIAg6gHD2bZLVLc+u2l5K0nOSLq0q2zWviL1P0v2Sfp7Lj5b09f7odxAEnU04+/Z5nbRidkTefj/wZGWnpDHAr4ADs+LlGODhPu9lEARBFeHsy3E58JH8//7AuVX7jgKOsX0/gO15tuvJJwRBEPQZ4ezLcR6wn6RlgI2AW6v2jQGm9EuvgiAIGhDOvgS2pwOjSaP6y7p7vC6ql5SbHhYEQdCMcPblmUgSNzu3pvweYLN2DmT7VNvjbI8bRvn52UEQBI0IZ1+ePwA/tD2jpvx44FuS1gGQNETSEX3euyAIgipiBW1JbD8BnFinfLqkw4FzJS1L0rf/ax93LwiCoAvh7NvE9sg6ZZOASVXblwKX1ql3dC92LQiCoCERxgmCIOgAwtkHQRB0ABHGGWhIaOn2Z+R4dvkpm0u9bZVyhsuUmzk079F/lWtvELHOF27v8zYf+8E2pezW+P7fS7dZ5lqF7l2vZenr667sZwPw/Kc2LWf4vxc03BUj+yAIgg4gnH0QBEEHEM4+CIKgA+g4Zy9pvqRpku6WdH6eC9+o7o6StqnaPkTSp0u2O1rSJ8vYBkEQdJeOc/bATNtjbY8B5gCHNKm7I7DQ2ds+xfZZJdsdDYSzD4KgX+hEZ1/NZODdknaXdKukqZKukbSKpNGkH4Kv5TuB7auTj0haS9IVkqZImixpvVx+hqSTJP1d0sOS9s5tHQtsn4/1tf54s0EQdC4d6+wlLQXsCswAbgS2sr0JSb74KNuPAqcAJ+Q7gck1hzgVONT2ZsDXgWrN+lWB7YDdSE4e4BvA5HysE2r6skj10rN69H0GQRBAZ86zHyFpWv5/MvB7YF1ggqRVgeHAI80OIGkkKbxzvqRKcfWk2ottLwDuldRyErvtU0k/Howa8mYXfytBEATF6ERnP9P22OoCSScDv7Q9UdKOwNEtjjEEeKn2OFVUrxhRgzpBEAR9RseGcWpYgUV5ZA+qKn8VWL62su1XgEck7QOgxMYt2qh7rCAIgr4gnH3iaFJIZjLwfFX5JcCelQe0NTYHAJ+TdBcpYcnHWrQxHZgn6a54QBsEQV/TcWGcBhLFfwH+Uqf8QVKO2QqTq/Y9Anyojs34eu3Zngu8r2y/gyAIukOM7IMgCDoA2TH5YyCxwvC3epuV923bbt7Tz/RCb5qz1KpvK2W3YJWVSre5YNq9pW2D+rx+xZqlbZff57lSdgtefbV0m0FjrvEFU2yPq7cvRvZBEAQdQDj7IAiCDiCcfRAEQQcQzj4jaU9JrtK4GS1pZtbLuU/SbZIOqqo/XtKvao4xSdK4/P+jkmbk172SfiypfOqaIAiCbhDOfhH7kzRy9qsq+6ftTWyvn8u/JukzbRxzJ9vvAbYA1iRLIgRBEPQ14exZqHWzLfA5ujr7hdh+GDgCOKzd49t+jaSguYek8lNRgiAIStJxi6oasAdwhe0HJb0gaVPghTr17gTWq9reV9J2VdvvbtSA7VckPQKsDdzaA30OgiAoTIzsE/uTpI3Jf/dvUK9W1GxCliwem0XR7mjRTl1RtGqJ4zkLZhbtcxAEQWE6fmQv6c3AzsAYSQaGAqarPn2FTYD7SrazPClb1YO1+6oljlcY/tZY5RYEQY8TI3vYGzjL9hq2R9tenaRn/47qSjlz1c+Bk9ttID8T+A1J5/7F7nc5CIKgPTp+ZE8K2RxbU/Zn4FvAWpKmAsuQJIpPtn16G8e+Tim7yRDgIuBHPdDfIAiCtul4Z297xzplJwEntbA7Azij0bFsj+6B7gVBEPQIEcYJgiDoAMLZB0EQdAAdH8YZaHjuvFJyxUOWWaZ8o8OGlTLznDml7EKmuJcYMrSU2cjdnyjd5AenPVvK7srNSspjz5pVyi6IkX0QBEFHEM4+CIKgAwhn3wZ1lDGHSDpJ0t1Z3fJ2Se/K+0L1MgiCAUM4+/aoVcbcF1gN2CirW+4JvFRVP1QvgyAYEISzL0gDZcxVgadtLwCw/US9FbKhehkEQX8Tzr44e5CVMYGKMuafgN0lTZP0C0mbNDK2/QpJhmHtPultEARBFeHsi7OYMqbtJ4B1gW8CC4BrJb2vyTFaql7OZXZP9jkIggCIefaFaKSMKeko27OBy4HLJf2bdAdwbZ1jFFK9HKWVQvUyCIIeJ0b2xWikjPleSatBmpkDbAQ8VmscqpdBEPQ3MbIvRiNlzDNI8fvKlMrbgOok5KF6GQTBgCCcfQHKKGOG6mUQBAOJCOMEQRB0AOHsgyAIOoAI4ywhdEcN8KV9Gi4PaMrIJ8upXi71txdK2Q0mllpzdGnbeQ8/Ws5wwfxSZi5pB3D5hiuWspu384al7Jb625RSdv1Bd66B0vyz8a4Y2QdBEHQA4eyDIAg6gHD2QRAEHUA4+yokvUPSXyQ9JOmfkk6UNFzSWEkfrqp3tKSv92dfgyAI2iGcfSYvfrqQtMp1bWAdYCRwDDAW+HBj67bbKpc/LgiCoCTh7BexMzDL9ukAtucDXwP+C/gZsG9Wt9w3199A0iRJD0s6rHIQSQdKui3X/V3FsUt6TdIPJd0KbN2n7ywIgo4nnP0iNgS6zOvKssSPAj8GJtgea3tC3r0e8EFSYpLvSxomaX1SQpNtbY8F5gMH5PrLAXfb3tL2jdXthOplEAS9TcyzX4SAeoqTjcr/mhUvZ0t6FlgFeB+wGXB7igoxAng2159P0tNZjFC9DIKgtwlnv4h7gL2qCySNAlYnOepaqofg80mfpYAzbX+zTv1ZOTQUBEHQ50QYZxHXAstK+jQsfIj6C5Ky5b+B5QseY29Jb83HWEnSGr3T3SAIguKEs8/YNilh+D6SHiIlGZkFfAu4jvRAtvoBbb1j3At8B7hK0nTgalKe2iAIgn4lwjhV2H4c2L3OrtnA5k3sxlT9PwGYUKfOyJ7oYxAEQRliZB8EQdABxMh+CeHKp6aVtv3gaj3XjyBRWrmyQxhM6pVlGWjXQIzsgyAIOoBw9kEQBB1AOPsgCIIOoOOdvaRvS7pH0vQ8tXJLSYdLWrbN44yWdHdv9TMIgqA7dPQDWklbA7sBm9qeLektwHDS1Mk/Am/UsRkaK2GDIBhsdPrIflXg+axxg+3ngb2B1YDrJF0HiytWSjpC0t35dXjtQSWtKWmqpM0lrSXpCklTJE2WtF7fvb0gCIJEpzv7q4DVJT0o6TeSdrB9EvAUsJPtnXK9hYqVwEzgM8CWwFbA5yUtzNgtaV2S4NlnbN9OEjg71PZmwNeB39R2IlQvgyDobTo6jGP7NUmbAdsDOwETJH2jTtVqxcrtgItsvw4g6cJsPxFYGfgLsJfteySNBLYBzs8qmABL1+lHqF4GQdCrdLSzh4VJSiYBkyTNAA6qU61asVJ19ld4GXgc2JakojkEeClr2wdBEPQbHR3GkbSupLWrisYCjwGv0ljl8gZgD0nLSlqOJJ42Oe+bA+wBfFrSJ3Pyk0ck7ZPbk6SNe/6dBEEQNKfTR/YjgZMlrQjMA/4BHAzsD1wu6emquD0Atu+UdAZwWy46zfZUSaPz/tcl7QZcLel1Uqaq30r6DjAMOA+4q9ffWRAEQRVKyr7BQGGUVvKWel/bdt3Txhlb2jYIgoHDNb5giu1x9fZ1dBgnCIKgUwhnHwRB0AF0esx+iaE7oZilVn9HKbt5jz9Rus0lnSEbr1/e9oVXS9nF+WjO0xeXOyfv+K9nS9nNf/4/pex6ixjZB0EQdADh7IMgCDqAQePsJc3PqpT3SLor69P0S/8lrSjpSyXsjpb09d7oUxAEQTMGjbMHZtoea3tD4P3Ah4Hv91NfVgTadvZBEAT9xWBy9gux/Sxp8dNX8qrUoZKOl3R71qX/AoCkHSVdL+lPWezsWEkHSLpN0gxJa+V6K0v6c7a/XdK2ufxoSX+QNEnSw5IOy104Flgr32kcn+seWdX+Dyp9zXr5D0i6Bli3Dz+mIAiChQza2Ti2H85hnLcCHwNetr25pKWBmyRdlatuDKwPvAA8TFrxuoWkrwKHAocDJwIn2L5R0juBK7MNwHokkbTlgQck/Rb4BjCmonkj6QPA2sAWJO2ciZLeC7wO7AdsQvqs7wSW/EzLQRAMOAats89URMk+AGwkae+8vQLJ+c4Bbrf9NICkf5JkjQFmkJw4wC7ABlXKlKMkVbRx/pr17mdLehZYpU4/PpBfU/P2yNz+8iSFzDdy+xPrvgnpYNKdCsvQVoKsIAiCQgxaZy9pTZL08LMkp3+o7Str6uwIXQTiF1RtL2DR+x8CbG17Zo09Nfbzqf+ZCfip7d/V2B8OtNSjCInjIAh6m0EZs5e0MnAK8CsncZ8rgS9KGpb3r5MVKYtyFfCVquOPbVG/VhXzSuCzWb8eSW+X9FaSQuaekkbkO4Xd2+hTEARBjzGYRvYjJE0jKUfOA84Gfpn3nQaMBu5UGo4/R5IaLsphwK8lTSd9JjcAhzSqbPs/km5SSjB+ue0jJa0P3JzvBl4DDswKmROAaSTp5MmNjhkEQdCbhOrlAKOs6mV3CLmEnifkEgYenSCXEKqXQRAEHU44+yAIgg5gMMXsg14ibv97ngV33Vfetgf7ESzizaeMLGW3zpUPlbK7b7NSZr1GjOyDIAg6gHD2QRAEHcCgcvaS3ibpPEn/lHSvpMskrdMP/RgvabU2bUbnqZpBEAR9zqBx9nn+/EXAJNtr2d4A+Bb15Qt6sx9DgfFAW84+CIKgPxk0zp6kYzPX9imVAtvTgBuz4uXdWclyX1ioeDlJ0gWS7pd0TlbI3FXSnyrHyPUuyf9/QNLNku6UdH7VithHJX1P0o3A/sA44JysejlC0mZZXXOKpCslrZrtNlPS3r8Z+HIffU5BEASLMZic/RjqK0Z+HBhLUrfcBTi+4mxJapOHAxsAawLbAlcDW1XJKewLTJD0FuA7wC62NwXuAI6oameW7e1s/zHvOyCrXs4DTgb2tr0Z8AfgmGxzOnCY7a2799aDIAi6x5Iw9XI74Fzb84F/S7oe2Bx4BbjN9hMAWWphdJYxvgLYXdIFwEeAo4AdSD8KN2XJg+HAzVXtTGjQ/rqkH6Krs91Q4GlJKwAr2r4+1zsb2LXeAUL1MgiC3mYwOft7gL3rlKtOWYVGipUTSGGVF0gSyK/mZwJX296/wbFeb1Au4J7a0bukFSmgeAmhehkEQe8zmMI4fwOWlvT5SoGkzYEXgX2VslWtDLwXuK3FsSYBmwKfZ9GI/RZgW0nvzsdetslMn2rVyweAlSVtne2GSdrQ9kvAy5K2y/UOKPxOgyAIephB4+yzlPGewPvz1Mt7gKOB/wOmA3eRfhCOsv1Mi2PNBy4lhVUuzWXPkWbZnJvVL28hZamqxxnAKTk0NJR0x3GcpLtICpfb5HqfIalp3gzMXPwwQRAEfUOoXg4w+kP1Mgg6gTkf2ryU3cbHTG1dqQ73bTavlF13CNXLIAiCDiecfRAEQQcwmGbjBL2Exo0pZec7Qv1hILHU28st6p735FN93ub8VVcqZdeda274FbeXsrvvinLtPf3f27Su1IDVrn+lnOHtFzTcFSP7IAiCDiCcfRAEQQcQzj4IgqADCGcfBEHQAQw4Zy/Jks6u2l5K0nOSLs3b4/P2VEkPZZXJbarqnyHpkaw2+aCksyS9vRf7O1bSh3vr+EEQBD3BgHP2JA2aMZJG5O33A0/W1JlgexPbawPHAhdKWr9q/5G2NyaJlE0FrpM0vJf6OxZoy9lLillQQRD0KQPR2QNcTlKjhKQff26jiravI4mIHVxnn22fADxDA8VJAEmvSTou69FfI2mLrIX/sKSP5jrLSDo9a+ZPlbRT/gH5IUmbZ5qkfbPt33Odv0taN9uPzxr5lwBX1bR/sKQ7JN0xt4t2WxAEQc8wUJ39ecB+kpYBNgJubVH/Thrr2BTZvxwpA9ZmJJGzH5PuKPYkOXPIyUdsv4f0A3Qm6fP7HulOY6ztCcD9wHttb5L3/aSqna2Bg2zvXN247VNtj7M9bhhLt3irQRAE7TMgwwm2p0saTXKqlxUwaSZzXGT/HKCydGIGMNv2XEkzgNG5fDtSkhJs3y/pMaCeKuYKwJmS1iZJHA+r2ne17Rda9CUIgqDHGagje4CJwM9pEsKpYhPgvm7sn+tFinALyDr4thew6Aex1Q9GhR8B19keA+wOLFO1r5EmfhAEQa8ykJ39H4Af2p7RrJKkHUjx+v+ts0+SDgNWZdHIvSw3kDXps879O0la9tXa9pBG9pUHyuO72WYQBEGPMGCdve0nbJ/YYHflgeiDwLeAvWxXj9yPz9ryD5JSFO5ke043u/QbYGgO7UwAxtueDVwHbFB5QAv8DPippJtIWvdBEAT9zoCL2dseWadsEim7FLbPICUPaWQ/vjtt2j663j7bs6gzUs8x+Fqh7OpY/ndzvTNo0u8gCILeZMCO7IMgCIKeY8CN7HsTSbfCYnMbP9XqucCSTlnZ2KErr1zKbv5zz5WyC5rz1B6jS9m99dflJY49Z245uw6Qx37bLW+Utv33lqPKGTZRce4oZ297y/7uQxAEQX8QYZwgCIIOIJx9EARBB9Anzl7S/Dw18e6sD7Nsk7o71qhYHiLp0yXbHS3pk2Vs22hjvKRyudmCIAj6iL4a2c/M2jFjSNIEhzSpuyOw0NnbPsX2WSXbHQ30qrMnTcdsy9mH6mUQBH1Nf4RxJgPvlrS7pFuzOuQ1klbJejiHAF/LdwLbSzpa0tcBJK0l6YqsTjlZ0nq5/AxJJ2WVyYcl7Z3bOhbYPh/ra/U6k0fmF0u6JOvgf0XSEblft0haKdcbm7enS7pI0ptyO+OAc3IbIyR9T9Lt+S7mVEnK9pMk/UTS9cBXe/HzDYIgWIw+dfZ5RLsrSWzsRmCrrA55HnCU7UeBU4AT8p3A5JpDnAocmtUpv05a1VphVZJY2W4kJw/wDWByPtYJTbo2hnQHsAVwDPBG7tfNQCWEdBbw/2xvlPv/fdsXAHcAB+Q2ZgK/sr15vosZkftTYUXbO9j+Rc3nEhLHQRD0Kn0VThghaVr+fzLwe1JikQmSVgWGA480O4CkkaTwzvl5sAxd58xfnIXL7pW0Spv9u872q8Crkl4GLsnlM4CNJK1ActTX5/IzgfMbHGsnSUcBywIrAfdUHW9CPQPbp5J+yBillVyvThAEQXfoK2c/0/bY6gJJJwO/tD1R0o7A0S2OMQR4qfY4VVQPiYsqVNazXVC1Xa162ZKsv/8bYJztxyUdTaheBkEwAOjPqZfV6pAHVZXXqkgCYPsV4BFJ+8BCRcuNW7RR91jtYvtl4EVJ2+eiTwGVUX51GxXH/ny+E9mbIAiCAUB/OvujSSGZycDzVeWXAHtWHtDW2BwAfC4rWt4DfKxFG9OBeUrJx+s+oG2Dg0hqmtNJeWcrGazOAE7JYarZJKnlGcDFNF28HARB0HdoUc6OYCAwSit5S72vv7tRiNDGGVg8++VtWleqw1t//ffSbcY10BhvO7a07XNjGy5Fasr0Xx0xxfa4evtiBW0QBEEH0DEje0kfBI6rKX7E9p790Z9GSHoOeKy/+xEEwaBkDdt1b7c6xtkHQRB0MhHGCYIg6ADC2QdBEHQA4eyDIAg6gHD2QRAEHUA4+yAIgg7g/wPerqAUlwrPpAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "classy = SVC(C = 50, kernel='rbf',class_weight='balanced')\n",
    "classy.fit(X_base_train, np.squeeze(Y_base_train))\n",
    "Y_base_pred = classy.predict(X_base_test)\n",
    "print(np.mean(Y_base_pred[:,None] == Y_base_test))\n",
    "fig, ax = plt.subplots(1,1, figsize=(5,5))\n",
    "plt.imshow(sklearn.metrics.confusion_matrix(Y_base_test, Y_base_pred, normalize='true'))\n",
    "plt.yticks([i for i in range(len(levels))], levels)\n",
    "plt.xticks([i for i in range(len(levels))], levels, rotation=90)\n",
    "ax.xaxis.tick_top()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "47153ed2-8d42-407f-821d-1cfb50a9df7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For BLR:\n",
    "Y_blr = np.zeros_like(all_Ys)\n",
    "all_Phis = create_design_matrix(all_Xs, site_ids = np.squeeze(all_Zs), all_sites=all_sites,basis='bspline')\n",
    "# For each feature:\n",
    "for b in range(6):\n",
    "    model_path = f'10_folds_results/fold_0_blr_sexcov/batch_{b+1}/Models/NM_0_0_estimate.pkl'\n",
    "    with open(model_path,'rb') as file:\n",
    "        model = pickle.load(file)\n",
    "    Yhat, s2 = model.predict(all_Phis)\n",
    "    this_Y = all_Ys[:,b]\n",
    "    this_Y_train = Y_train[Y_patient.columns[b]].to_numpy()\n",
    "    this_Y_s = (this_Y - this_Y_train.mean(axis = 0))/this_Y_train.std(axis=0)\n",
    "    warped_this_Y_patient = model.warp.f(this_Y_s,model.blr.hyp[1:model.blr.warp.get_n_params()+1])\n",
    "    this_Z = (warped_this_Y_patient - Yhat)/np.sqrt(s2)\n",
    "    Y_blr[:,b] = this_Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "id": "a77c2721-4a19-49f9-8662-a1711d2972eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_features_blr = np.concatenate((all_Xs, all_Ys, Y_blr, all_Zs), axis = 1)\n",
    "# all_features_blr = Y_blr.copy()\n",
    "all_features_blr = (all_features_blr-all_features_blr.mean(axis=0, keepdims = True))/all_features_blr.std(axis = 0, keepdims = True)\n",
    "X_blr_train, X_blr_test, Y_blr_train, Y_blr_test = train_test_split(all_features_blr, all_labels, stratify = all_labels, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "id": "ebfd5fed-9bfe-4a8a-86e8-5aa6c57bee7f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6594936708860759\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAF7CAYAAAAzPisLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABI2UlEQVR4nO2debhcRbW+3y8hQAIERJALikTmSQgQZpBB9IqCgoCAoES9Ig4gcoGf01WcQVQu4IBclEkuRBAQkBkJBGQmIWFGGS4KKCgzIWT4fn9UddKn08PuffoMnV7v8/RzumvXqqo+u3t17bWrviXbBEEQBIs2I4Z6AEEQBMHAE84+CIKgBwhnHwRB0AOEsw+CIOgBwtkHQRD0AOHsgyAIeoBw9kEQBD1AOPsgCIIeYLGhHkAQdBOSRgIfAMZR9f2x/ZOhGlMQFCGcfRC0x6XA68AMYN4QjyUIChPOPgja4222NxrqQQRBu0TMPgja4wpJ7x3qQQRBu8TMPgja41bgIkkjgNmAANseO7TDCoLmKFQvg6A4kh4F9gBmOL48QRcRYZwgaI9HgHvD0QfdRoRxgqA9ngYmS7oCmFUpjKWXwXAnnH0QtMdj+bF4fgRBVxDOPggKkjdUrWX7wKEeSxC0S8Tsg6AgtucCK0qKGX3QdcTMPuhqJG0LTLP9qqQDgU2BE20/UcB2NdJM/VpJo4HFbL/cwuxx4GZJlwCvVgojZh8Md2JmH3Q7vwBek7QxcDTwBHBWKyNJnwYuAH6Zi94GXFygv6eAy0jfnWWqHkEwrIl19kFXI+lu25tK+gbwN9u/qpS1sJsGbAHcZnuTXDbD9jsL9rsMaTPVK/18C0EwKMTMPuh2Xpb0FeBA4A/5JuqoAnazbL9ReSFpMaDlzEfShpKmAvcC90m6S9IGJcceBINGOPug29mXtN79U7afAd4KHF/A7gZJXwVGS3oPcD5J0bIVpwJH2F7N9mrAfwL/U27oQTB4RBgn6Emyts2ngPeS9G2uAk5rtTNW0j22N25VFgSdJF+xHmb7hNJthLMPuhFJN9neTtLL9A2/FBImk7QncLntWc3q1bG7CLgbODsXHQhMsL1HO+0EQbtImmx7x9L24eyDXkTS6cDOwI3AecBVtucUsHsT8C1gu1x0I/At288P1FiDAEDS94BlgUn0XfZ7dyH7cPZBNyNpDeCvtmdJ2hHYCDjL9gsFbEcBu5Li/tsB19j+jwZ1z7b9MUlftH1ip8YfBEWRdH2dYtveuZB9OPugm8lLKCeQcsJeBVwCrGP7/QXtRwHvAz4BbG97xQb17if9MFwC7EgKF83H9r9KvYEgGCRiB23Q7cyzPSfH4P/b9sl5aWRTJL0P2A/YCZgMnAZ8pInJKcCVwOrAXfR19s7lQTCgSPoAsAGwZKXM9reL2IazD7qd2ZL2Bw4Cds9lRdbZTyTF6j9T5Cat7ZOAkyT9wvZnyw42CMoi6RRgDGmCchqwN3B7YfsI4wTdjKT1gUOAW2yfK+kdwL62j+1wP2NtvyRp+XrHI4wTDDSSptveqOrv0sCFtgvlRA5nH3Q9WYVy7fzyIduzm9QttWRT0mW2d5P0WLbrE8axHWGcYECRdJvtLSXdCnwY+Ccpa9paRewjjBN0NXkFzpkkNUoBq0o6yPaN9erb3i7/bUu8zPZu+e87+jHcIOgPl0lajrRD/G7SpOO0osYxsw+6Gkl3AR+1/VB+vTZwru3NWtidbftjrcoa2L4VWI2qyVKjH5cg6BSSlqjcX5K0BOkm7etFNwbGzD7odkZVHD2A7YfzcspW9BEvy0JoTX8gcr3jSOvy7wfmVrolba4KgoZkyYOrbO9SsolbSPkayA5+lqS7K2WtCGcfdDt3SvoVC+QLDiAtjaxLVsisCKC9VCkG3iCJnLViD9I6/rZkFoLA9lxJr0la1vaLRe0k/RtJ4G+0pE1YcL9oLGl1TrF2IowTdDP5cvbzpB2wIs2wf97KGUv6ge2vlOjvCmCf0LEPyiDpt8BWwDX0lTw4rInNQaSlwhOAO6sOvQycYfvCQn2Hsw96kbwJ64+VGVa+8bWj7Ytb2P0O2Bi4jiStDDT/sgZBhey4F8L2mQVs97L9u9J9h7MPupmcg/YYFr5h2nQppKRptsfXlE2tZK1qYlf6yxoE/SFPSL4BvCsX3QB8u2hIKGL2QbfzK+BLpDj93BZ1q6mXuKfl98H2mTk5+durbwwHQREkrQX8AFifvpIHRfZp/IqUIa0i6/Ex4HTSmvuWRKaqoNt50fYVtv9h+5+VRwG7OyX9RNIaklaXdAJNbuxWkLQ7MI2kk4Ok8ZIu6dc7CHqJ04FfAHNIsgdnsWBxQSvWsP1N24/mx7doQ5MpwjhBVyKpstzsI8BI4EL6xtCbanxLWgr4L6CyDO5q4Hu2X21sNX9d/87A5DKJyoOhR9IRzY7b/skA9n2X7c2qPzOSptjevoDtLcBRtm/Kr7cFfmR76yJ9Rxgn6FZ+XPN6QtVzkxxyQ7JT/7KkpdtcWTPH9otSX4XjNuyDoaeye3odYHOSbDUkIb2G+yUkzaDJuba9UYG+X88pMR+R9AXgb8Bbigwa+CxwpqRlSSvP/kVapVOImNkHPYmkbUhbzZe2/XZJG5MUMD/Xwu5XpJU4Xwb2Ag4jbew6ZKDHHHQWSVcDe9l+Ob9eBjjf9vsa1F8tP/18/lu9t+O1IlLDkjYHHgCWA75DWiv/Q9u3tTHusQC2X2pVt49dOPugm5H0fdKX5YX8+k3Af9r+egu720gSsZdUhWPutb1hC7sxwNdIicohJUz5ru3X+/VGgkFH0oPAxjUSBPfYXreF3c22t21V1sB2H9vntyprYLsc8HFSop7qlWeFlv3GDdqg29m1OgWhUy7YQlmqbD9ZU9R0NU/e7n6J7a/Z3jw/vh6Ovms5G7hd0jGSvgncRrph2oqlJFVyEFeuEpcq2Ge9jXxFN/ddTnL0M0iLCSqPQkTMPuh2RtYIRI0Glihg92T+kjpLJB9GurxuSNnt7sHwxPb38o7oys3RT9humeUM+BTw6xw7N/Ai8MlmBpJ2JU1C3irppKpDY0krc4qwpO2mN5ebEc4+6HZ+A1wn6XTSF++TJMnjVhwCnEjSHPkraTXO55taJF4HZkgqvN09GNaMAV6yfbqkFSW9w/ZjzQxs3wVsnGPnKvjD/xRJ6uCD9J2Nv0zaJ1KEsyV9GriMvivPCiXOiZh90PXkfLK7kFYoXG37qiZ1j7P9/4rGSevYxw7aRYQcuplAErZbW9IqpBu0TWPvklYCvg+sYnvXnC1ta9u/KtDnKNLntFCynRrbzwPfA15gwaqgwolzwtkHXU1evnZOjtUXqT+DJAl7m+1C0rDBoomkacAmwN1VN+mnt1pCmUM/pwNfs71xlseeWmSvhaQdSPcFHicn2wEOKpIPQdJfgC1tP9eqbj0ijBN0O/8G3JF1vX9N0gtvNoO5EniOdJPtJXI6wspfN0hLWKGOFk/FLtISdh9v2LYkw/yNdkVYwfZvs1w2tudIKirV8RPgva5JtkOBXArAfcBrBftZiHD2QVdj++uS/ou0FPITwE+zjOyvbP+ljsnXbR8l6fe2P1Siy7JaPMHw47eSfgksl2Phn6RYmr9XJb2ZHEqRtBXpJm0RyibbgfR5mybpekooroazD7qePDt7BniGtLLhTcAFkq6xfXRN9Uq2n7Y2pFTxou0ryo82GC7Y/pGk95A+C+sA37B9TQHTI0i7bteQdDOwIrBPwW7bSrZTw8X5UYqI2QddjaTDgINIoZnTgIttz65sSbe9Rk39e0kJm78BHFXbnhskguivFk8w/KjcrG9VVsduCdIsex1SGO8hYIQLZC9TyWQ7RZD0O9t7NTwezn74I2k34HLb84Z6LMMNSd8mhWyeqHNsPdsP1JRtR5pNfYQFmigVbLvueul86dwI226qxRMMPyTdXXuTvuAN2np2C5UNNmqRjyHCON3BfsCJSlmSTq91YD3OmaTwDZJ2BDYCzrL9Qr3/U1YMvEnSnUWWylXZ7dSZ4QZDjaTPAp8DVpc0verQMsDNTez6nQs2T9y+w8I3+JsuDChI05l7zOy7hLyBY3/STUiTln6dWxFx6lXy8rkJpG3kV5Fm6+vYbimZIGlDFk4i0XS7fL4x903SZbiBm0jZgopo6AfDgLzz9U2kJCJfrjr0crMNSuqbC/YOFjj7l4AzG4UAa9r4MynZyIwWq8baptXVRTj7LkLSCsCBwOGkrf1rAifZPnkoxzWUVD7gko4CXrd9cqvL2Wz3TWBHkrO/HNgVuMn23i3sriHFWX+Tiw4g5a7dpbFVMFzJekcr0VdY7P+a1B8B7G/7nJL9XQ+8eyBCsq0+9yGE1gVI2l3SRcAfgVHAFrZ3JSW+PnJIBzf0zJa0P+km7WW5rMhStr2BdwPP2P4E6X9ZRFNnedvfsf1YfnyXJFc7LJH0xSJlvUjekPd34BrgD/lxWTOb7KQ/049ujwYul/QVSUdUHgXH2+pcNr2xHM6+O9gHOMH2RraPt/0PANuv0UKAqQf4BLA1KcvUY5LewYJZdzNm5i/unBwi+wfFUrxdL2k/SSPy4yMkJzFcqSfvMHGwBzFMOZwU8tvA9jvzo0gCkmskHSlpVUnLVx4F+/weaWPUkqR7BJVHEZqeS9tXNzOOME7Q9WTVynVJMfSHbL9RwObnwFdJN7//E3gFmJZn+c3sXibJ2VYuw0ewQBCtUzfa+k2+2vko6d7ClKpDywBzI+w0P6TyHttFVScrdvWE0grtos4LAya0qldj05FzGc6+C5D0YeA4Uvoy0dk7+F2NpA8ApwB/If1f3kHKOFV445OkccBY29Nb1e0WlLIqvYM6NyGB6e06uEWRvLlpHdKVWfWeiYHMQXss8MdWs/Aam46cy3D2XUC+g797LLlcGKVsQ7vZ/nN+vQbwB7fINpTrfhB4V355g+1LC/ZZbTfZdtM4bzCwZH2ZXwAr2d5Q0kbAB/P9lGZ236xXbvtbLezGkHbRvt32wZLWIoWDWn4Oqq4MZwGzGcSJWzj7LkAFU571IpJutP2uqtciOe53NTGrzLA2ByqrKvYH7rTdNGtQA7u7bH+5sdXQUeaqUCn93Vr55cMe5olaJN1A2g39S7eRYrLKfimnBPRF+5tEkjj4eP5xGQ3cYnt8+6NfqO0NbN/X4Fi/rvDD2Q9j8skF2IGk7ngxfS83W67rXdSR9AvSBpXfkmL2+5C2r98MTeUPpgPjK0vg8hK8qQV2T5ayGyrauSrM9z5OBfYAHiM5k9WAi4BDitwLKTnGEcBWtv9U0v4O25tXLz2UNK2V85W0NUnYrt2k83fanlDT3z22Ny4z/pq2G66V7+8VfuygHd7sXvX8NRYkuYbk2Hre2ZNWNfyd9IMI8CywPOl/1+p/tBxQ2USzbBt9lrUbCv7ehnP4OmnZ6qqVzXqSlgF+BvxXfnQc2/Mk/Zi0qqoMz+XwXUWFcm/g6QJ2/w38O1k2w/Y9kppeEWbeyLP5Sn9rUDUJ6ydqcqydc7kQ4eyHMZWVIZK2td1nG7eSrnpPk2fV022fUML8B8DUvCJDpBj8VwfQbqi4M4cdLqb1VeGHSXs4Xquq97KkzwG3MkDOPnO1pL2AC0vsLP086YpkXUl/I12VHFDE0PaTKfI3nyKy1ceQ8iKsKukcYFs6t5y12Xtv51wuRIRxuoB6l3attkb3CpKud0ndGkkrk+LvImWuemYg7YYCpdy8tdh1BN/URARM0gwXyMRUlqobl3OBmbQRj5a0me27lJKPjMg/ULu3uuEu6QJSMpGfAluRks5PsL1fgT7fnG0E3OqS2aPqtNssjFP4XNYjZvbDmBxT3AZYsWaX3ViSzG4Af5L0U2ASfROAN5UclnQd8GPbl1SVnWr74IGwGypa7RuorS7pTdQPJQyo4qrtohuL6vE/kg6yPQNA0n6kBDOtVleVSjov6RJSdqlL2rmxW5CG90XaPJcLETP7YYxSvsodSR/KU6oOvQxcavuRAehzMZJOTGXp4gPAlcN1XbbqSw/bLSSHJT0KPEla8/ytXNbyaqms3VDRzrJESY+TnHo9Z19o01A/xilS6OUdtr8jaVVgZdu3F7BdHbgg228HfJy0HHdAVhHl7+W+wAeA20kTjctsv97Crl/frbJLTOfbh7Mf/khazXX02gegn1WA60k3t6aSvvSbkFYC7WT7qYEew2ChlLN2C+AkUtLnA4HrCzj7UnZDRX+XJQ4WeVXVPGBn2+vlK4yrbW9e0H5tUiz7SWAP2zML2LwDOJSkmFothPbBgn2OBHYGPg28r8Vy1n5/t/p7LiOMM4yRdCkL7vgvdLzoh7INvg/8wvZ/14zjMNKNyXraHEOKpJVI417F9q6S1ge2dmuteuUZ1eckTSRJFb+pSJcl7YaKMbZvr/n81J1JakE2rrq0Co31ky2d1Eun5r6ez0tBGyJpBn1vaC5PCm/eJokCy2EvJi29vJQ2w1R5Nc7upBn+pqS8Cs3oxHer8LmsRzj74c2PBrm/rWxPrC20fZKkh+rUHw6cQdL2/1p+/TDpsrqVs58fFrN9RnYcLeO1/bAbKtpZlvjjqueb0Tc3qkmz2Ib0czPW7DxTroxzRVo74N3aaL8er9s+qV2jvCJmS9KKnJ+RdlG3Gmsnvltll5imcUcYZ9Gl3R13aqKH3ezYUNLuhhpJY22/pAYqhW6QvKKs3VCT49mnkm70P09alnig7cdb2BU+3+rAZixJB9B3lrw38HXb5xccw8bA9vnlFNv3FLD5KOnH6WrayCcs6X3ANbaLLNOs2PT7u1X2XFaImX0XoKS98QMWzqrU6obZD2lvx92yWrBrt88QSCuAhiOv5mVwldnOVkCzGeX/kmaEd2Wb6mti01jmuKzdkGL7UWCX6mWJRU3b6Kbfm7FsnyPpLlKOAZHi7oU+t0qa7p9mwQa63+QVUq2S+rwT+BjpiqUyM295BWP7SknbKAnoVcf6m2U56/d3qx/nMnUUM/vhj6SbSKnwTiDFCT9BOnd1hZyq7NrS1Gmwjnc+/V36NRDkOPPJwIbAvcCKwN5eBBQs1SKphZuoM0o60PZvGrXRzDbbF15hJOleajZj5fKlSWvQi2rUtJU1qspuOuk+zav59VIkrZpW0hcPAhsVufKosTsbWAOYxoJNWLZ9WBOb0t+t/p7LCjGz7w5G275OkvKqnGMkTSH9ADSjrR13w9GZt8L23Xkp3DqkWdJDtmcXsZX0VhYkfq60d+NA2ZWgP2vPl2rSRt0ZnqSTq469TVKfeHYTZzav1tHn+q9IKjSblHQo6fP8d5IDVR5LEc0h0Xfna8W+FfeQpC/+UWSMVUwA1ncbM+V+frfaPpf1CGffHbyuJBb1iFIqtb+R4vCtGEubmjqS1gEOpu9a4FNtP1xm4IPEFixYPrdpXonRKnH4caQY8f1Uzc5I+WU7blcGt5DabWH7y/z0WheX2riz6vldDerU7U7934z1RZJMcJnE7b8mrcC5KL/eg9Y36CFdRTwo6Q76ToZarXK7l7RksvDNUSj/3Sp5LhfuP8I4wx9Jm5M+GMsB3yE58eNt39rhfrYm/RCcCtzNgrXAnwY+3On+OkGZS+ps9xDpEr4tAauydmWQ9EPgUdun1JR/Cfg3201zjua6Ay61oQ5sxlL5rFEjSLIFr5M2VAm40fbUArY71Cu3fUOBsY4nbagq9CPRie9Wf89lOPsuQgV1tyUdbfuHNZfl82nkCCVdARxne3JN+Q7Al52SnHccSQeRZ3a56AHgpFaz82z7AG1eUme7K4B9bL8y0HZl35+k+4ENa5f1ZQc3vVksXAukNg4n3eupMBbY0w3kePtzLspQFYfegJJZoyTdYruUYqZSFqi1bF+rlJRkZKsbn2V+JPrz3Sp7LmuJME4XoCrdbaCI7nZlFcOdDY43Yo3aDyOkD7GkU5uM7wekhOgLxT4lHddsBirp46QP8REsmPFsChxfJBxDm5fUVT+ArwHTlLRuqp1Lox/Csnb9eX+udfS5cJ5UZ5ddXxYnfV4Wo2+s9yXSssaOjTU7zBec19VL2okUSnkc+FmLG6CVsf1ffiyeH1A8Hl1KMVPSp0lhleVJV4dvJe2jeHczu/x9WIkkhgdwe73Pfg2lvluZts9lPWJm3wVIuo10Ui/xAG55l3SX7c0aHGumxvci8ALwOdt/KGqXj98K7OeatcJ5Wdt5trdqYFfZXbwMbVxS55lrQ2zX3QnZD7tS7y/XuQP4qGs0kJSW4p7rAomr1YbURj/OxW2kGeZTksYD15KWCm8EzLb9HwX63sc1a+rrlTWwrShmziGFcwopZkqaRrrfc1vV96qluqekjwDHA5NzX9sDR9m+oIlNqe9WTb3C57IeMbPvElxCd1tpF+L/Y+H1+Y3WEa9auwKj0hRp1tOIR0nrlc+RtCtwpBeIQrWagY6tdS55jI9LavZlLbW7uOKU8/K81503xigt+1ui03aUf38A3wCukPRdFtwwnQB8hTQDL8Jrko4nhUlafQbKjnW0F2i7HAj82vaPc7hpWsFxfgWodez1yhbC5RUzZ9l+o/K9UhIqKzL7/RqweWU2n79n15LE2BpR9rtVTTvnciHC2XcHT0rahrTqYXGS7naRDSfnkKQDPkBSzjyIlMmpEUc1OdYsJGTb90rakrSR605JH3Wxte7NBKsaHqvER+uFifKKmaY32YDrgF2ASux9NGkn5TYdtiv1/gBsXyFpD9J5OTQX3wfs5SznW4DKZ2A3Wn8Gyo61+gd9Z5KTLhRuypOD9wNvrXGGYymo+yLpOtvvblVWhxskfRUYLek9wOdoLYsMaUNTddjmn8CIFjZlv1vVtHMuFyLCOF2ApBVIutu7kL5YVwNfdItlapVLR1UlpZB0g+26N5j6Mb4+l6H5C3wKSRnyQDfZCi7pNeDP9Q4Bq9teqs6xhn3nsoZJOKrqLCSpUK+sv3b9fX/9pZ3PQNmxSjoRWJl03+SDwNq2ZyslebnETZQr8/2n8SRZj++SZtZzSevtJ9t+vontksAYkprkjiz40RkLXGF7vUa22X4E8CnS0mQBVwGntYr759n1RiRNe0hLcac3uzfVCfr7fY6ZfRfglAWnUJq1Giqbi56W9AHgKeBtjSorJWVoNo5GS8tUU+8KSZuR1j+32hTT9AvZCEmfJc3E1lDaQVlhGaBI4upXJW3qrIMiaQItZtol7eq9P5HOQ9N0hv04H9W08xkoO9bDSQ5vZWA7L9jUthbp5mcz7id9thcHPpn7W5UkbndZC9vP5L5Xoe++gJdJUg1NyVceFwMX2245Q5a0JklL/igl6YPKUs9bSLPuZraDfS4XHkPM7IcvarB0skKjFSBV9rsBU0hfnpNJM55vuSrLUk39Z0l64OcCt7GwE68bGlFK0vywc3q+vKpjL+AJYJJrNoK0Il/J/LPZDEvSsiRp4R8Ax5JywQLc5GJrrCeQLomfIv2PVwH2td10M1FZu2w7Hvgo8BGSiNXvbP+0Sf1S56OmjbY+A2XH2sTuQjfRqJF0AmmlyRFeoKszlnRPZqbtLzax3ZyUYWpv2yfnm+h7kVYBHePGonYi7db9AswXCJwLnGz72036uwz4am14Mn8mvml79ya2Q3YuqzuJxzB9kGJylcfjNa8PamK3JGnG81PS7Gexgv2NBN5HUh2cSrqs3qCA3d3A8vn5u0iOcC/SBrALWthuRVrVcCFpk8m9wDOkLezvK9D3F4EZwLeAbwPTgUML2O2TvywbkoS6/gBs2mk7YG3SjdYHSNr3hwJPDOT5KPsZKDvWfr7HR8iTzjrv/ZGB+NyRUhZeQ8qKVSlbnRTG+VITu3ubHJsxnM5l3XbKGMVj8B/A1DbqTgJ+kz8YFwMnluhvCWAi6QZQU+cJ3FP1/GekWVXl9bQWtneSYqb7kGRbt8rl6xZ5z9m5L1X1eilS/LSlXf67HUnq4EOkJXgdtSPtLL0BWLOq7NGBPB9lPwNlx9qf90i6Imz7WH8+d9nZrlCnfMVmnzngz2WODcW5rPdodQc5GD60E29b3/aBTpoae7NA57slkpbI8cjfkJJynEQTLZ3MyLxsDdKGlD9WHWt1X2gx21c7rad+xnnbuO0Hiw6ZciJYFZsPAKfY/j0LNvN00m4v0pXK9ZL+R1JFwrcQJc8HlPsMlB1rf97j/Tns1wdJBwKtPgNlP3ejnO6D9cEpbj+qid0dShuxasf6KQpoCQ3yuVyIuEG7aDJf9dH2nBar3+Yj6UxSeOIKUizw3oL9nUtaxvYc6WbllNzemjTXloe+Qlm1NzqL/MCdTjkRrL9J+iVphdNxkpag9fK5tu1sXwRclNfn70EKIayklHP1IttXN7Ltx/mAEp+BsmPtz3skOb0LJX2SBbkCNictad2zxZDLfu6a7ehtduxw0vs8gL77HhZvNdbBPpd1x5AvE4JhiNLOwMoJGkPaqg8tdghKmgu8WlV3dLZtZTevyq76g9FyR6JS0pCVSUmiK7riawNLu0nmn6qxVo+z0ueStpvNtCptbEr7IlhjSDHUGbYfycsE39nCMZW2q2ljeVLYal832RDTz/NR6jNQdqz9tZO0M2mzkID7bF9XsJ+2P3c1/5s+hyjwmVOSg6jsXr/P9h+b1c82Q38uw9kHQRAs+kTMPgiCoAcIZx8EQdADhLPvMiQd3A12Q9FnvMfh1WcvjLWb3mM4++6j7IdrsO2Gos94j8Orz14Ya9e8x3D2QRAEPUCsxhlmjHnT4l5ulTENj7/2/BuMedPCe3heub/57/ZsZjGqgey6RjS2fcOvs7iWrHvM85rnkm7W56JgN2B9NllHPduvM6rB+WjZZyPbFj6g6Wdn8carFN+YO5PFR46ue8xvzK5bXqTPRcFuoPp8meefs71ivWOxqWqYsdwqY/j0ee0rEP9p4yKbP+szYkw5ld15M19vXanTzGuZs2V4MGJkaVONLG9bBs9uto+oOYutsmopuzlPPFm6z9L045yUYgg+q9f6goaZrCKMEwRB0AOEsw+CIOgBFnlnL+mVmtcTJbXU5W7Q1o5Z07ryfJuqY2dIKpzpPQiCYDBZ5J39ALIjrfOVBkEQDAt62tlLWlHS7yTdkR/b5vItJP1J0tT8d50au3GkhL9fkjRNUkVy9F25/qOVWb6ksyV9qMr2HElFUpAFQRB0jF5YjTNa0rSq18sDlTReJwIn2L5J0ttJmWrWI+lovyvLie4CfJ+k2Q2A7cclnQK8YvtHMF/TemWS+uK6uY8LgNNIkq+/V0qltw0p09R88o64gwGWXbn+UrUgCIL+0AvOfqbt8ZUXkiaSNKghaZKvX6UPPVbSMsCywJmS1iLJkbaU2c1cbHseKSHDSpByS0r6maS3AB8m5fKcU21k+1TgVIBVNlguNj4EQdBxesHZN2MEsLXtPkkzlBJ9X297zxyymVywvVnVzVQ9Pxs4ANgP+GTp0QZBEJSkp2P2wNWkDPMASBqfny4L/C0/n9jA9mVgmYL9nEHKcoPt+9obYhAEQf/pdWd/GDBB0nRJ95NuugL8EPiBpJtJWeHrcSmwZ80N2rrY/jvwACmFXhAEwaCzyIdxbC9d8/oM0kybnHR43zo2twBrVxX9Vy6fTA7p2H4Y2KiqzpRG/eZUdmuRcmYGQRAMOr0+sx9w8mqeB4GTbbdKvh0EQTAgLPIz+6HG9rXA24vWf/XhxbntvW9ru5+Tn7iobZsKh662bWnboAH9EMFyt4i9MUSCZmXpov/rQBAz+yAIgh4gnH0QBEEPEM4+CIKgBwhnXxBJc/Myy3slnZ9X2DSq+0FJXx7M8QVBEDQjnH1xZtoeb3tD4A0WrMlfCNuX2D528IYWBEHQnHD25ZgCrClpd0m3ZXXMayt6ONWa+Vnn/qRaNcwgCILBJJx9m0haDNgVmAHcBGxlexPgPODoBmYVNczdgIVm/JIOlnSnpDvfmDdzIeMgCIL+Euvsi1MtlTwF+BWwDjBJ0srA4sBjDWwXUsOsplr1ctlRbwnVyyAIOk44++L0kUqG+eqYP7F9iaQdgWMa2DZSwwyCIBgUIozTP6rVMQ9qVjEIgmAoCWffP44Bzpc0BXhuiMcSBEHQkAjjFKRWPTOX/R74fZ3yM1igrDmxVTtBEAQDTczsgyAIeoBw9kEQBD1AhHGGGZ4zh7nP/rNtu/7IFC8zZYVSdjMPHF3Kbs6TT5WyA8rL1I5olHBsGNJNUrwl/68jRi9Zym7eq6+WskudlvwMlDwfi73treX6A+b+/dlyhm80PhQz+yAIgh4gnH0QBEEPEM6+AZIs6eyq14tJelbSZfn1xPx6qqRHJF0laZuq+mdIekzSPZIelnSWpPLXdUEQBP0gnH1jXgU2lFQJTL+HBRuoKkyyvYnttUiaNxdKWq/q+FG2NybJKkwFrpe0+EAPPAiCoJZw9s25AvhAfr4/cG6jiravJ+nbHFznmG2fADxDElELgiAYVMLZN+c8YD9JSwIbAbe1qH83sG4/jgdBEAwI4eybYHs6MI40q7+8gEkrkbO6x6sljmf30UwLgiDoDOHsW3MJ8COahHCq2AR4oN3jtk+1PcH2hFEsUW6UQRAETYhNVa35NfCi7RlZxrguknYgxet3qnNMwKGkJCZXDswwgyAIGhPOvgW2/wqc2ODwvpK2A8aQEpfsZbt65n68pP/Kx28FdrLdZI9bEATBwBDOvgENVC4nA5Pz8zPIypYN7CcOyMCCIAhKEDH7IAiCHiCcfRAEQQ8QYZzhhoRGtX9aPKu8UuJre5bLcf7wV8qpP7zjojeXsgMY9Y+XS9nNe+Kvpfssg2f1xhLaxd5STjF17vMvdHYgRRhkNdE5f63dcD+0xMw+CIKgBwhnHwRB0AOEsw+CIOgBwtl3CElfk3SfpOmSpknaIf+tfrwk6bihHmsQBL1H3KDtAJK2BnYDNrU9S9IKwOK2x1fVeSdJX+e/h2SQQRD0NOHsO8PKwHO2ZwHYfq76YFbNPAf4vO2nh2B8QRD0OBHG6QxXA6vmjFQ/zzo51fwQuNn2JfWM+6he+vUBH2wQBL1HOPsOYPsVYDOSENqzwCRJEwEk7QrsAvxnE/sFqpdachBGHARBrxFhnA5hey5JN2eypBnAQZL+APwS+JDt14ZyfEEQ9DYxs+8AktaRtFZV0XjgCZI88sm2pw7JwIIgCDIxs+8MSwMnS1oOmAP8GTgBuBF4u6QDqupeY/uowR9iEAS9TDj7DmD7LmCbOodapSkMgiAYFCKMEwRB0APEzH64YQ+6YqJnllvuufoFM0vZPb3tUqXsAJZ7tFyO3rEvvVLKbu4/ny9l100s9m8rlbad88zfOziSYCCJmX0QBEEPEM4+CIKgBwhnHwRB0AOEs+8nkuZmRct7JN0taZtcPk7STElTJT0g6XZJBw31eIMg6E3iBm3/mVlRt5T078APgIo2zl9sb5KPrQ5cKGmE7dOHZKRBEPQsMbPvLGOBuss3bD8KHAEcNqgjCoIgIGb2nWC0pGnAkiSp452b1L0bWLe2UNLBJBE1lmTMAAwxCIJeJ2b2/Wem7fG21wXeB5wlqdHO2brlfVQvKbeOPAiCoBnh7DuI7VuAFYAVG1TZBHhg8EYUBEGQCGffQSStC4wE/lnn2DjgR8DJgzysIAiCiNl3gErMHlKY5iDbc3MkZw1JU0nx/JdJcsexEicIgkEnnH0/sT2yQfnjwOjBHU0QBEF9IowTBEHQA4SzD4Ig6AEijBMw79VXS9nNenO5ZaJLPzWvlB3Ai+PqRs1aMnZquUTunv1GKbtuImSKe4OY2QdBEPQA4eyDIAh6gAF39pKWk/S5qtc7SrpsoPsti6TJkiZ0oJ0/dWI8QRAEnWAwZvbLAZ9rVakokrriPoPtegnIgyAIhoSOO3tJR0i6Nz8OB44lbS6aJun4XG1pSRdIelDSORUtGUmbSbpB0l2SrpK0ci6fLOn7km4AvijpDEmnSJoi6WFJu+V6G2Td+GmSpktaS9J3JH2xanzfk3RYfn60pBlZi/7YqrexT27nYUnb57pLSjo9158qaadGfebyV/LflSXdmI/fW2kvCIJgMOnoLFnSZsAngC1Ju0lvAw4ENqzSfN+RpBGzAfAUcDOwraTbSFICH7L9rKR9ge8Bn8zNL2d7h9zGGcA4km78GsD1ktYEDgFOtH2OpMVJ0gW/Ai4ETpQ0AtgP2ELSrsAewJa2X5O0fPX/xfYWkt4PfBPYBfg8gO13ZlmEqyWt3aDPaj4KXGX7e5JGQshaBkEw+HQ6JLIdcJHtVwEkXQjUm8nebvuvuc40kuN+AdgQuCZP9EcCT1fZTKpp47e25wGPSHqUJB18C/A1SW8DLrT9CPC4pH9K2gRYCZhq+5+SdgFOt/0agO1/VbV9Yf57Vx5b5b2dnOs+KOkJYO0GfVZzB/BrSaOAi21Pq/1nhMRxEAQDTafDOI2kfWuZVfV8LulHR8B9WS54vO132n5vVb3axeCufW37f4EPAjOBqyRVtOVPAyaSrjp+XTXW2jZqx1cZW6X+QjTps3L8RuBdwN+AsyV9vE4bIXEcBMGA0mlnfyOwh6QxkpYC9iSFaZYpYPsQsKKkrQEkjZK0QZP6+0gaIWkNYHXgoZz671HbJwGXABvluheRtOY3B67KZVcDn5Q0Jve3PM25ETgg110beHuLPsl1VwP+Yft/SCGlTVv0EwRB0HE6GsaxfXeOp9+ei06zfZekmyXdC1wB/KGB7RuS9gZOkrRsHtt/A/c16O4h4AZSaOYQ26/nOP+BkmYDzwDfrmr7euAF23Nz2ZWSxgN3SnoDuBz4apO393PgFEkzgDnARNuzGvVZxY7AUfn4K8BCM/sgCIKBRnajSMbwJf+gXGb7goL1R5BSAu5TJ6Y+rBir5b2l3j3UwyjE67ttUcrujbHlLyhf/bdytqte9LdSdnMee6KUXRAMBdf6grts190ntMjvoJW0PvBn4Lrh7uiDIAgGiq7YoFSL7Ylt1L2fFNMPgiDoWbrS2QfDg3mjii6+6strK5W/oFz2sbml7MqGY0ZusE4pu7n3PVTKLggGikU+jBMEQRCEsw+CIOgJhsTZS5oo6acNjl0uabkO9rW9pPuyNs1oScfn18e3tl6ora/WvA5lyyAIuoJhN7O3/X7bL3SwyQOAH+VduTOBzwCb2j6qRFt9nH0oWwZB0C101NlLWkrSH7KK5L2S9pW0uaQ/5bLbJVV2064i6UpJj0j6YVUbj0taQdIheTY+TdJjeVMUkvbPypP3Sjquyu4Xku7Ms/Zv5bL/AD4CfENJXfMSYCngtjy2FSX9TtId+bFttltaCxQup0vaS0kVc3Qezzm5XkXZclIWTauM5YxsMzJfSdyR2/lMJ//fQRAERen0apz3AU/Z/gBA3gk7FdjX9h2SxpI0ZADGk9QvZ5FkB062/WSlIdunkHasjgL+CPxE0irAccBmwPMk5ck9bF8MfM32v7Ky5HWSNrJ9mqTtqNqAJemVKgXO/wVOsH2TpLeTpBTWA/4LeNH2O3O9N9n+naQvVGxrOA/YF7hcSfny3cBngU/ldjaXtARws6SrbT/Wr/9yEARBm3Ta2c8AfpRn3JeRlCyftn0HgO2XAJRULa+z/WJ+fT+wGvBknTZPBP5o+1JJHwIm2342251DEhm7GPhIVo9cDFgZWB+Y3mK8uwDr5/EAjM1XHruQpJDJ436+RTtXkGQeliD94N1oe6ak9wIbZRkIgGWBtYA+zl6hehkEwQDTaW2ch5U07d8P/IAkNtZKWRL6qkvOR9JE0o/AFypF9RqS9A7gSGBz289nOYUlCwx5BLB1juVXt9dMEXMhsi7PZODfSTP8c6vGe6jtqxrZZvtTgVMhySUU7TcIgqAonY7ZrwK8Zvs3wI+ArUix+c3z8WVUMK1g/tE4Ejgw69ZDSoayQ47pjwT2J4mhjSVJIL8oaSVg14JDvpoFPyRkYbR65W/KT2fnsFI9ziNJKG/PAmXNq4DPVmwkra2kBhoEQTCodDqM807geEnzgNmkuLWAkyWNJsXrdynY1heA5UlZqADutP0fkr4CXJ/bvdz27wEkTSUpZD5KklUuwmHAzyRNJ/0vbiRlnvpuLr+XdNXxLVJCk1OB6ZLutn1ATVtXA2cBl9h+I5edRkp+cne+WniWlB0rCIJgUOlK1ctFmW5SvXxtzy1L2b00rjZzY3HGPl5OLmHMRbeVsgu5hKCb6GnVyyAIgiCcfRAEQU8QqpdBaZ56VznVy+UeKB86HPP0zNaVOkiEY4JFhZjZB0EQ9ADh7IMgCHqAcPZBEAQ9QMecvQZXtniypLrLizrQ9vyxSjpM0gMV4bMgCIJuZVBu0Np+f+tanUXSYrbntGtXM9bPAbsWFS4r22cQBMFA03Jmr6GVLX5F0o8l3S3pOkkrVg1tn9z3w5K2z/UnSjpf0qUkRcylJP06SwxPzUJqlXoXthjrKaRE5ZdI+pKk5SVdrCRVfKukjXL9YySdKulq4Kz8+kxJV+e2Pizph/n9XdlEbiEIgmDAKBLGqcgWb2x7Q+BKYBLwRdsbk+QPqmWL9yXJJuwradXqhmyfkiWCNwf+Sl/Z4p2z/eaS9sgmSwF3296UpIHzzarmFrO9BXB4TfnWwEG2dwa+RlLM3BzYiSTlUNGmaTXWQ4CngJ1sn0CSTJhqeyNSEpOzqqpvBnzI9kfz6zWADwAfAn4DXJ/lkmfm8j5IOlhJi//O2X304YIgCDpDEWc/A9hF0nF5Bv12amSLq0IX19l+0fbrQEW2uB7zZYtJjn+y7WdzOxXZYoB5pB8WSE5zu6o2Lsx/7yLpz1S4xva/8vP3Al+WNA2YTFLCfHubY62wHXB2fs9/BN6spNcPSQ+negH4FbZnk/53I0k/kOTX1WMlt3eq7Qm2J4xiiRbDCIIgaJ+WMfuhki1uNJw6fdX282p1d8BetvvsjJG0ZZGx1lBvnJXxvFpTPgvA9jxJs71AgGhegX6CIAg6TpGY/VDJFlfGV0n88VHgpmJvaz5XAYdmxUkkbdKmfTU3kvLZImlH4LlKMpYgCILhThEnPWSyxaQZ8waS7gJeJMXY2+E7wH+TZIkFPA7s1mYbFY4BTs9yyK8BB5VsJwiCYNAZ1hLHSvlilx7qcQwm3SRx/OcTtiplt9wD5TR1AFaYXhsxK8itrTJUBkH3ExLHQRAEPc6wvlnYa7N6ACQ0avG2zTz7jdaVGjBiqXKZEtf80q2l+1zUWe+u8l+tBz+zbik733lv6T6DxpT5PkL/vpMDQczsgyAIeoBw9kEQBD1AOPsgCIIeIJx9B5C0pyRLWje/HidpZtbjeSBr+MRSzSAIhoxw9p1hf9KGr/2qyv5iexPb6+XyL0n6xJCMLgiCniecfT+RtDSwLfAp+jr7+dh+FDgCOGwQhxYEQTCfcPb9Zw/gStsPA/+StGmDencD5dbUBUEQ9JNw9v1nf+C8/Py8/LoeDbeN9pE49uudHl8QBMHw3lQ13JH0ZpIO/4aSTJIzNvDzOtU3AR6o147tU4FTAcaOePPw1a8IgqBriZl9/9gbOMv2arbH2V4VeAx4W3UlSeNIiqEnD/4QgyAIYmbfX/YHjq0p+x0pk9UakqaSEqa8DJxs+/RBHl8QBAEQzr5f2N6xTtlJwEmDP5ogCILGRBgnCIKgBwhnHwRB0ANEGGe4YQ+6NOqIFd9cym7eqyUTifQA1/3f2qVtx65eTtl76TtLdxk0YbhJFZclZvZBEAQ9QDj7IAiCHiCcfZtkdcuzq14vJulZSZdVle2ad8Q+IOlBST/K5cdIOnIoxh0EQW8Tzr59XiXtmB2dX78H+FvloKQNgZ8CB2bFyw2BRwd9lEEQBFWEsy/HFcAH8vP9gXOrjh0NfM/2gwC259iuJ58QBEEwaISzL8d5wH6SlgQ2Am6rOrYhcNeQjCoIgqAB4exLYHs6MI40q7+8v+31Ub1kVn+bC4IgWIhw9uW5hCRudm5N+X3AZu00ZPtU2xNsTxjFEp0aXxAEwXzC2Zfn18C3bc+oKT8e+KqktQEkjZB0xKCPLgiCoIrYQVsS238FTqxTPl3S4cC5ksaQ9O3/MMjDC4Ig6EM4+zaxvdBedtuTgclVry8DLqtT75gBHFoQBEFDIowTBEHQA4SzD4Ig6AEijDPM0KjFWGyFldq2m/PM30v3Oefx/yttG9RnlT3vH/Q+Hzlr01J2a3387g6PJADQEuVX1j33sXLnkv+5oOGhmNkHQRD0AOHsgyAIeoBw9kEQBD1Azzl7SXMlTZN0r6Tz81r4RnV3lLRN1etDJH28ZL/jJH20jG0QBEF/6TlnD8y0Pd72hsAbwCFN6u4IzHf2tk+xfVbJfscB4eyDIBgSetHZVzMFWFPS7pJukzRV0rWSVpI0jvRD8KV8JbB9dfIRSWtIulLSXZKmSFo3l58h6SRJf5L0qKS9c1/HAtvntr40FG82CILepWedvaTFgF2BGcBNwFa2NyHJFx9t+3HgFOCEfCUwpaaJU4FDbW8GHAlUa9avDGwH7EZy8gBfBqbktk6oGct81cs35s3s6PsMgiCA3lxnP1rStPx8CvArYB1gkqSVgcWBx5o1IGlpUnjnfEmV4upFtRfbngfcL6nlonnbp5J+PFh28be4+FsJgiAoRi86+5m2x1cXSDoZ+IntSyTtCBzToo0RwAu17VRRLUqvBnWCIAgGjZ4N49SwLAvyyB5UVf4ysExtZdsvAY9J2gdAiY1b9FG3rSAIgsEgnH3iGFJIZgrwXFX5pcCelRu0NTYHAJ+SdA8pYcmHWvQxHZgj6Z64QRsEwWDTc2GcBhLFvwd+X6f8YVKO2QpTqo49Bryvjs3Eev3Zng28u+y4gyAI+kPM7IMgCHqAnpvZdwXqjnu6s987oZTd6Bl/Ld3nnKefKW27qFNWvfKRM9pKmdyH9b7yt9aV6tAL59GzZrWu1IA3n3ZLB0eSiJl9EARBDxDOPgiCoAcIZx8EQdADhLPPSNpTkqs0bsZJmpn1ch6QdLukg6rqT5T005o2JkuakJ8/LmlGftwv6buSyqeuCYIg6Afh7BewP0kjZ7+qsr/Y3sT2ern8S5I+0UabO9l+J7AFsDpZEiEIgmCwCWfPfK2bbYFP0dfZz8f2o8ARwGHttm/7FZKC5h6Slu/HUIMgCEoRSy8TewBX2n5Y0r8kbQr8q069u4F1q17vK2m7qtdrNurA9kuSHgPWAm7rwJiDIAgKEzP7xP4kaWPy3/0b1KtdAD8pSxaPz6Jod7bop+4C+pA4DoJgoOn5mb2kNwM7AxtKMjASMH316StsAjxQsp9lSNmqHq49FhLHQRAMNDGzh72Bs2yvZnuc7VVJevZvq66UM1f9CDi53Q7yPYGfk3Tun+//kIMgCNqj52f2pJDNsTVlvwO+CqwhaSqwJEmi+GTbp7fR9vVK2U1GABcB3+nAeIMgCNqm55297R3rlJ0EnNTC7gzgjEZt2R7XgeEFQRB0hAjjBEEQ9ADh7IMgCHqAng/jDDvmzGXe8y+0baZRi5fuUkuWU3FY/PrppezmzH6jlF0wMKz96RmlbT99//2l7H6xzjrlOpw3t5xdEDP7IAiCXiCcfRAEQQ8Qzr4N6ihjjpB0kqR7s7rlHZLekY+F6mUQBMOGcPbtUauMuS+wCrBRVrfcE3ihqn6oXgZBMCwIZ1+QBsqYKwNP254HYPuv9XbIhuplEARDTTj74uxBVsYEKsqYvwV2lzRN0o8lbdLI2PZLJBmGtQZltEEQBFWEsy/OQsqYtv8KrAN8BZgHXCfp3U3aaK16SfmM9EEQBI2IdfYFaKSMKelo27OAK4ArJP2ddAVwXZ02iqlejnhzqF4GQdBxYmZfjEbKmO+StAqklTnARsATtcahehkEwVATM/tiNFLGPIMUv68sqbwdqE5CHqqXQRAMC8LZF6CMMmaoXgZBMJyIME4QBEEPEM4+CIKgB4gwzjDDNvNef31Q+3x+v81K2b26St2VpC1563F/KmXXTSy2+rjStnMefbxj4yiC+6FC+ou11ixlN2fn8aXsFvvjXaXsgpjZB0EQ9ATh7IMgCHqAcPZBEAQ9QDj7KiS9TdLvJT0i6S+STpS0uKTxkt5fVe8YSUcO5ViDIAjaIZx9Jm9+upC0y3UtYG1gaeB7wHjg/Y2t2+5rZKfaCoIgKEI4+wXsDLxu+3QA23OBLwH/AfwQ2DerW+6b668vabKkRyUdVmlE0oGSbs91f1lx7JJekfRtSbcBWw/qOwuCoOcJZ7+ADYA+67qyLPHjwHeBSbbH256UD68L/DspMck3JY2StB4pocm2tscDc4EDcv2lgHttb2n7pup+qlUvZ4fqZRAEA0Css1+AgHqKk43K/5AVL2dJ+gewEvBuYDPgjhQVYjTwj1x/LklPZyGqVS/HavlQvQyCoOOEs1/AfcBe1QWSxgKrkhx1LdVT8Lmk/6WAM21/pU7913NoKAiCYNCJMM4CrgPGSPo4zL+J+mOSsuXfgWUKtrG3pLfkNpaXtNrADDcIgqA44ewztk1KGL6PpEdISUZeB74KXE+6IVt9g7ZeG/cDXweuljQduIaUpzYIgmBIiTBOFbafBHavc2gWsHkTuw2rnk8CJtWps3QnxhgEQVCGmNkHQRD0ADGzX0S46qlppW3/fZVydsuV7nHRZ7CVK7uNUK8cfGJmHwRB0AOEsw+CIOgBwtkHQRD0AD3v7CV9TdJ9kqbnpZVbSjpc0pg22xkn6d6BGmcQBEF/6OkbtJK2BnYDNrU9S9IKwOKkpZO/AV6rYzMydsIGQdBt9PrMfmXguaxxg+3ngL2BVYDrJV0PCytWSjpC0r35cXhto5JWlzRV0uaS1pB0paS7JE2RtO7gvb0gCIJErzv7q4FVJT0s6eeSdrB9EvAUsJPtnXK9+YqVwEzgE8CWwFbApyVtUmlQ0jokwbNP2L6DJHB2qO3NgCOBn9cOIlQvgyAYaHo6jGP7FUmbAdsDOwGTJH25TtVqxcrtgItsvwog6cJsfwmwIvB7YC/b90laGtgGOD+rYAIsUWccoXoZBMGA0tPOHuYnKZkMTJY0AzioTrVqxUrVOV7hReBJYFuSiuYI4IWsbR8EQTBk9HQYR9I6ktaqKhoPPAG8TGOVyxuBPSSNkbQUSTxtSj72BrAH8HFJH83JTx6TtE/uT5I27vw7CYIgaE6vz+yXBk6WtBwwB/gzcDCwP3CFpKer4vYA2L5b0hnA7bnoNNtTJY3Lx1+VtBtwjaRXSZmqfiHp68Ao4DzgngF/Z0EQBFUoKfsGw4WxWt5b6t1t2/VPG2d8adsgCIYP1/qCu2xPqHesp8M4QRAEvUI4+yAIgh6g12P2iwwRihlejBjTltpGR5j32kIbvoMqvvPYHaXsvrnBDqXshtv5iJl9EARBDxDOPgiCoAfoGmcvaW5WpbxP0j1Zn2ZIxi9pOUmfK2F3jKQjB2JMQRAEzegaZw/MtD3e9gbAe4D3A98corEsB7Tt7IMgCIaKbnL287H9D9Lmpy/kXakjJR0v6Y6sS/8ZAEk7SrpB0m+z2Nmxkg6QdLukGZLWyPVWlPS7bH+HpG1z+TGSfi1psqRHJR2Wh3AssEa+0jg+1z2qqv9vVcaa9fIfknQtsM4g/puCIAjm07WrcWw/msM4bwE+BLxoe3NJSwA3S7o6V90YWA/4F/AoacfrFpK+CBwKHA6cCJxg+yZJbweuyjYA65JE0pYBHpL0C+DLwIYVzRtJ7wXWArYgaedcIuldwKvAfsAmpP/13UBkWg6CYNDpWmefqYiSvRfYSNLe+fWyJOf7BnCH7acBJP2FJGsMMIPkxAF2AdavUqYcK6mijfOHrHc/S9I/gJXqjOO9+TE1v146978MSSHztdz/JXXfhHQw6UqFJRn8JXtBECz6dK2zl7Q6SXr4HySnf6jtq2rq7Ah9BOLnVb2ex4L3PwLY2vbMGntq7OdS/38m4Ae2f1ljfzjQUo8iJI6DIBhoujJmL2lF4BTgp07iPlcBn5U0Kh9fOytSFuVq4AtV7Y9vUb9WFfMq4JNZvx5Jb5X0FpJC5p6SRucrhd3bGFMQBEHH6KaZ/WhJ00jKkXOAs4Gf5GOnAeOAu5Wm48+SpIaLchjwM0nTSf+TG4FDGlW2/U9JNyslGL/C9lGS1gNuyVcDrwAHZoXMScA0knTylEZtBkEQDCShejnMKKt6GQwvQi5h+NELcgmhehkEQdDjhLMPgiDoAbopZh8EXUOEVIYfn/3BYa0r1WHUBc+Wslv2/X8uZTdQxMw+CIKgBwhnHwRB0AN0lbOX9G+SzpP0F0n3S7pc0tpDMI6JklZp02ZcXqoZBEEw6HSNs8/r5y8CJttew/b6wFepL18wkOMYCUwE2nL2QRAEQ0nXOHuSjs1s26dUCmxPA27Kipf3ZiXLfWG+4uVkSRdIelDSOVkhc1dJv620ketdmp+/V9Itku6WdH7VjtjHJX1D0k3A/sAE4Jysejla0mZZXfMuSVdJWjnbbaakvX8L8PlB+j8FQRAsRDc5+w2prxj5YWA8Sd1yF+D4irMlqU0eDqwPrA5sC1wDbFUlp7AvMEnSCsDXgV1sbwrcCRxR1c/rtrez/Zt87ICsejkHOBnY2/ZmwK+B72Wb04HDbG/dv7ceBEHQPxaFpZfbAefangv8XdINwObAS8Dttv8KkKUWxmUZ4yuB3SVdAHwAOBrYgfSjcHOWPFgcuKWqn0kN+l+H9EN0TbYbCTwtaVlgOds35HpnA7vWayBUL4MgGGi6ydnfB+xdp1x1yio0UqycRAqr/Iskgfxyvidwje39G7T1aoNyAffVzt4lLUcBxUsI1csgCAaebgrj/BFYQtKnKwWSNgeeB/ZVyla1IvAu4PYWbU0GNgU+zYIZ+63AtpLWzG2PabLSp1r18iFgRUlbZ7tRkjaw/QLwoqTtcr0DCr/TIAiCDtM1zj5LGe8JvCcvvbwPOAb4X2A6cA/pB+Fo28+0aGsucBkprHJZLnuWtMrm3Kx+eSspS1U9zgBOyaGhkaQrjuMk3UNSuNwm1/sESU3zFmDmws0EQRAMDqF6OcwI1csgGBieO7jcOolRe3SPXEKoXgZBEPQ44eyDIAh6gG5ajRMMEE//5zatK9Vh5R//qcMjCYKBY4VTb2ldqR6nljN7/qAh2F5zxgUND8XMPgiCoAcIZx8EQdADhLMPgiDoAcLZB0EQ9ADDztlLsqSzq14vJulZSZfl1xPz66mSHskqk9tU1T9D0mNZbfJhSWdJeusAjne8pPcPVPtBEASdYNg5e5IGzYaSRufX7wH+VlNnku1NbK8FHAtcKGm9quNH2d6YJFI2Fbhe0uIDNN7xQFvOXlKsggqCYFAZjs4e4AqSGiUk/fhzG1W0fT1pcdTBdY7Z9gnAMzRQnASQ9Iqk47Ie/bWStsha+I9K+mCus6Sk07Nm/lRJO+UfkG+TtHmmSdo32/4p1/mTpHWy/cSskX8pcHVN/wdLulPSnbP7aLcFQRB0huHq7M8D9pO0JLARcFuL+nfTWMemyPGlSBmwNiOJnH2XdEWxJ8mZQ04+YvudpB+gM0n/v2+QrjTG254EPAi8y/Ym+dj3q/rZGjjI9s7Vnds+1fYE2xNGsUSLtxoEQdA+wzKcYHu6pHEkp3p5AZNmMsdFjr8BXJmfzwBm2Z4taQYwLpdvR0pSgu0HJT0B1FPFXBY4U9JaJInjUVXHrrH9rxZjCYIg6DjDdWYPcAnwI5qEcKrYBHigH8dne4Ei3DyyDr7teSz4QWz1g1HhO8D1tjcEdgeWrDrWSBM/CIJgQBnOzv7XwLdtz2hWSdIOpHj9/9Q5JkmHASuzYOZelhvJmvRZ5/7tJC37am17SDP7yg3lif3sMwiCoCMMW2dv+6+2T2xwuHJD9GHgq8Betqtn7sdnbfmHSSkKd7L9Rj+H9HNgZA7tTAIm2p4FXA+sX7lBC/wQ+IGkm0la90EQBENO6NkPM4ZCzz6E0IKg8wyFENpdZ/xn6NkHQRD0MsNyNc5AIek2WGht48da3RdY1IkZ+qLBn3+zSSm7NQ+c2uGRBAAr3P5cadt/bfbmDo4k0VPO3vaWQz2GIAiCoSDCOEEQBD1AOPsgCIIeYFCcvaS5eWnivVkfZkyTujvWqFgeIunjJfsdJ+mjZWzb6GOipFUGso8gCIL+Mlgz+5lZO2ZDkjTBIU3q7gjMd/a2T7F9Vsl+xwED6uxJG6facvahehkEwWAzFGGcKcCaknaXdFtWh7xW0kpZD+cQ4Ev5SmB7ScdIOhJA0hqSrszqlFMkrZvLz5B0UlaZfFTS3rmvY4Htc1tfqjeYPDO/WNKlWQf/C5KOyOO6VdLyud74/Hq6pIskvSn3MwE4J/cxWtI3JN2Rr2JOlaRsP1nS9yXdAHxxAP+/QRAECzGozj7PaHcliY3dBGyV1SHPA462/ThwCnBCvhKYUtPEqcChWZ3ySNKu1gork8TKdiM5eYAvA1NyWyc0GdqGpCuALYDvAa/lcd0CVEJIZwH/z/ZGefzftH0BcCdwQO5jJvBT25vnq5jReTwVlrO9g+0f1/xfQuI4CIIBZbDCCaMlTcvPpwC/IiUWmSRpZWBx4LFmDUhamhTeOT9PlqHvmvmLs3DZ/ZJWanN819t+GXhZ0ovApbl8BrCRpGVJjvqGXH4mcH6DtnaSdDQwBlgeuK+qvUn1DGyfSvohY6yWjy3NQRB0nMFy9jNtj68ukHQy8BPbl0jaETimRRsjgBdq26miekpcVKGynu28qtfVqpctyfr7Pwcm2H5S0jGE6mUQBMOAoVx6Wa0OeVBVea2KJAC2XwIek7QPzFe03LhFH3XbahfbLwLPS9o+F30MqMzyq/uoOPbn8pXI3gRBEAwDhtLZH0MKyUwBqvcVXwrsWblBW2NzAPCprGh5H/ChFn1MB+YoJR+ve4O2DQ4iqWlOJ+WdrWSwOgM4JYepZpGklmcAFwN39LPPIAiCjhCql8OMoVC9DBYNQhtneDFyvbVK25bVxrnj7CND9TIIgqCX6ZmZvaR/B46rKX7M9p5DMZ5GSHoWeGKoxxEEQVeymu0V6x3oGWcfBEHQy0QYJwiCoAcIZx8EQdADhLMPgiDoAcLZB0EQ9ADh7IMgCHqA/w9EHTViSoA4iwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "classy = SVC(C = 50, kernel='rbf', class_weight='balanced')\n",
    "classy.fit(X_blr_train, np.squeeze(Y_blr_train))\n",
    "Y_blr_pred = classy.predict(X_blr_test)\n",
    "print(np.mean(Y_blr_pred[:,None] == Y_blr_test))\n",
    "fig, ax = plt.subplots(1,1, figsize=(5,5))\n",
    "plt.imshow(sklearn.metrics.confusion_matrix(Y_base_test, Y_blr_pred, normalize='true'))\n",
    "plt.yticks([i for i in range(len(levels))], levels)\n",
    "plt.xticks([i for i in range(len(levels))], levels, rotation=90)\n",
    "ax.xaxis.tick_top()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "id": "508a0bed-7d81-45c8-81fe-68cded0d13fa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Right-Cerebellum-White-Matter\n",
      "EstimatedTotalIntraCranialVol\n",
      "Right-Lateral-Ventricle\n",
      "WM-hypointensities\n",
      "rh_S_interm_prim-Jensen_thickness\n",
      "Brain-Stem\n"
     ]
    }
   ],
   "source": [
    "# For SHASHb:\n",
    "Y_shashb = np.zeros_like(all_Ys)\n",
    "this_X_S = (all_Xs - X_train.to_numpy().mean(axis = 0))/X_train.to_numpy().std(axis=0)\n",
    "this_Y_S = (all_Ys - Y_train.to_numpy().mean(axis = 0))/Y_train.to_numpy().std(axis=0)\n",
    "for b in range(6):\n",
    "    print(Y_train.columns[b])\n",
    "    model_path = f'10_folds_results/fold_1_SHASHb_sexcov/batch_{b+1}/Models/NM_0_0_fold1SHASHbsexcov.pkl'        \n",
    "    with open(model_path,'rb') as file:\n",
    "        model = pickle.load(file)\n",
    "    this_X_t = bspline_transform(this_X_S, model.hbr.bsp)\n",
    "    MAP_path = f'10_folds_results/MAPS/MAP_{b}_SHASHb_fold1_sexcov.pkl'\n",
    "    with open(MAP_path,'rb') as file:\n",
    "        MAP = pickle.load(file)\n",
    "    this_Z = get_single_zscores(this_X_t,this_Y_S[:,b][:,None],all_Zs, model,MAP)\n",
    "    Y_shashb[:,b] = np.squeeze(this_Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "id": "c3f6ff4e-adc5-4389-b939-59edcc097a15",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_features_hbr = np.concatenate((all_Xs, all_Ys, Y_shashb, all_Zs), axis = 1)\n",
    "# all_features_hbr = Y_shashb.copy()\n",
    "all_features_hbr = (all_features_hbr-all_features_hbr.mean(axis=0, keepdims = True))/all_features_hbr.std(axis = 0, keepdims = True)\n",
    "X_hbr_train, X_hbr_test, Y_hbr_train, Y_hbr_test = train_test_split(all_features_hbr, all_labels, stratify = all_labels, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "id": "76f235ed-4ec4-43fe-b9cc-6e20420ae2b8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6430379746835443\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD4CAYAAADl7fPiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAARGklEQVR4nO3dfYxVdX7H8c+H4UFEFFkEBXygirZoVmop1ZhutVYXqZHdzdpCu661m+Capa1Nm9S2qW7SP2rTWJstRstWozauun1gJZGqhDZRs6KgQRRFmUVcxqHgIyggOMy3f8yBzG+8F865556Ze5n3KyFz7znfe87vzh0/nnPvub+vI0IAcMiIoR4AgNZCKABIEAoAEoQCgAShACAxcqgHUMvYCWPixKnjctXue6PYtj0ifw5Gb2+xjQ8HdnXbLvhJmMeMyb/p/fuLjuaY9pn26EDsr/litmQonDh1nH7v4Stz1b72K8X+wx1xfL6wkaTefZ8V2nYhvQer23ZRIzpyl7ojf21R8fmBQvUdZ52du/bgWz8rOpz8Cvz+Cqvo7+SFWF13HacPABKlQsH2PNtv2u60fWuN9bb9g2z9BtsXldkfgOo1HAq2OyTdLelqSbMkLbI9a0DZ1ZJmZv8WS7qn0f0BGBxljhTmSuqMiC0RcUDSo5IWDKhZIOmh6LNG0gTbp5XYJ4CKlQmFaZK29bvflS0rWiNJsr3Y9jrb6/Z9xDvFwFApEwq1Ps4Y+JlSnpq+hRHLImJORMwZe3L+j5oANFeZUOiSdHq/+9MldTdQA6CFlAmFtZJm2p5he7SkhZJWDKhZIenb2acQF0vaFRHbS+wTQMUavngpInpsL5H0lKQOSfdHxEbb383W3ytppaT5kjol7ZV0Y/khA6iSW3GSlZNGTY5LJl2Xq/afXlxeaNt/dOaljQwJOKa8EKu1Oz6seZkzVzQCSBAKABKEAoAEoQAgQSgASBAKABKEAoAEoQAgQSgASBAKABKEAoBES87mHD09OrhjZ67aot9lGP/spNy1+741ttC2e97ZdvSiQ4rOAFxkVt9W2naVCox7xLj8s3hLUu+ePfmLK/x9jzzz9KMX9XOwe0e+ws/rT9XPkQKABKEAIEEoAEgQCgAShAKABKEAIEEoAEgQCgASZXpJnm77f22/YXuj7T+pUXOZ7V2212f/bis3XABVK3NFY4+kP4uIl22Pl/SS7VUR8fqAumcj4poS+wEwiMr0fdguaXt2+xPbb6ivT+TAUGgp+xblb0m3ecn0Qts+6/GJuWtH7tpXaNvKe/mqpN59nxXadOwvcJlzkUuipWKXABfc9shTp+Su7dn5fqFtF1L0d1JAoUvnizhCa4emvKdg+yxJvyzphRqrL7H9iu3/tn3+EbZxuMHs56LBLDBUSn8hyvYJkv5T0i0RsXvA6pclnRkRn9qeL+knkmbW2k5ELJO0TJJO9MTW61ADDBOljhRsj1JfIDwcEf81cH1E7I6IT7PbKyWNsp3/a4oABl2ZTx8s6T5Jb0TEP9apOTWrk+252f4+aHSfAKpX5vThUknXS3rV9vps2V9JOkM63GD2m5Jutt0jaZ+khdGKzSsBHFbm04fnJNWfqaGvZqmkpY3uA8Dg44pGAAlCAUCCUACQIBQAJAgFAImWnOK9Sr0ffZy79pSXi02v/eGs43PXjn93dKFtH/fOu7lrPbLYyxr7K7ysvMh05gW+yyBJBz/4qJJxDHccKQBIEAoAEoQCgAShACBBKABIEAoAEoQCgAShACBBKABIEAoAEsPvMuc9e3LXfjK9WGae/FZP7trdZxT71Y8dnf+y6IMfFbj8t4X0/F/+aexRHY4UACQIBQCJslO8b7X9atYncl2N9bb9A9udtjfYvqjM/gBUrxnvKVweEfV6cl2tvuYvMyX9mqR7sp8AWlTVpw8LJD0UfdZImmD7tIr3CaCEsqEQkp62/ZLtxTXWT5PUv0NmV7bsC+glCbSGsqcPl0ZEt+3JklbZ3hQRz/RbX6svRM1mMPSSBFpDqSOFiOjOfu6UtFzS3AElXZL6z2k2XVJ3mX0CqFaZXpLjbI8/dFvSVZJeG1C2QtK3s08hLpa0KyK2NzxaAJUrc/owRdLyrH/sSEk/iognbX9XOtxLcqWk+ZI6Je2VdGO54QKoWpleklskXVhj+b39boek7zW6DwCDb9h996GIk94uNi141xX5z8amPtNbaNvt+n0GtB8ucwaQIBQAJAgFAAlCAUCCUACQIBQAJAgFAAlCAUCCUACQIBQAJLjM+Qi6v1JrOoj6znlsb0UjKejiLxerX7OhmnGgLXGkACBBKABIEAoAEoQCgAShACBBKABIEAoAEoQCgESZKd7PyxrLHvq32/YtA2ous72rX81tpUcMoFJlZnN+U9JsSbLdIeld9TWEGejZiLim0f0AGFzNOn24QtLPIuKdJm0PwBBp1ncfFkp6pM66S2y/or52cX8eERtrFWUNahdL0nE6Xh41OteO4/MDhQaad7uSdM6frim07ZZR8LsMp/x0Qu7a7X9zdqFtj1z9UqH64aDI32DRv+9mKH2kYHu0pGsl/XuN1S9LOjMiLpT0z5J+Um87EbEsIuZExJxRPq7ssAA0qBmnD1dLejkidgxcERG7I+LT7PZKSaNsT2rCPgFUpBmhsEh1Th1sn+qs2aTtudn+PmjCPgFUpNR7CraPl3SlpJv6LevfYPabkm623SNpn6SFWX9JAC2qVChExF5JXxqwrH+D2aWSlpbZB4DBxRWNABKEAoAEoQAgQSgASBAKABKtOcV7RGWXd3ZMOzV3bc/Wn1cyhlaz9n9+KXft2POLTXs/ZXXR0bSfkWedUai+1f+uOFIAkCAUACQIBQAJQgFAglAAkCAUACQIBQAJQgFAglAAkCAUACQIBQAJt+LsaCeNPCUuOXFBrtqDH++qeDQo4+1HLsxdO2PRKxWOpD15zJhC9e9ff1Guuk3L79Le97bV/CILRwoAEkcNBdv3295p+7V+yybaXmV7c/bz5DqPnWf7Tdudtm9t5sABVCPPkcIDkuYNWHarpNURMVPS6ux+Iusvebf6+kLMkrTI9qxSowVQuaOGQkQ8I+nDAYsXSHowu/2gpK/VeOhcSZ0RsSUiDkh6NHscgBbW6HsKUyJiuyRlPyfXqJkmaVu/+13ZMgAtrMqZl2q9s1n3o46kweyIcVWNCcBRNHqksMP2aZKU/dxZo6ZL0un97k9XX+fpmvo3mB3tsQ0OC0BZjYbCCkk3ZLdvkPR4jZq1kmbanpF1pl6YPQ5AC8vzkeQjkp6XdJ7tLtvfkXSHpCttb1ZfL8k7stqptldKUkT0SFoi6SlJb0j6cURsrOZpAGiWo76nEBGL6qy6okZtt6T5/e6vlLSy4dEBGHSteZlzx6S4+IRrc9X2fvJJZePomPkLheo/mzExd+2op9cVHc4xr/OuiwvVn/vg7ty1vetfLzqcY9oLsVq740MucwZwdIQCgAShACBBKABIEAoAEoQCgAShACBBKABIEAoAEoQCgAShACBR5SQrg2NER7Hy0aNy1/Zu7Sq07VGbtxSqR+rc24p9ifbGlzbkrr3vF88pNpjeg8XqjyEcKQBIEAoAEoQCgAShACBBKABIEAoAEoQCgESjDWb/wfYm2xtsL7c9oc5jt9p+1fZ620xKCLSBRhvMrpJ0QUR8WdJbkv7yCI+/PCJmR8ScxoYIYDA11GA2Ip7O+jpI0hr1dX8CcAxoxmXOfyjpsTrrQtLTtkPSv0TEsnobSXpJ6vjKpm7/8Lq5uWv3TK05A3Zd0/7+p0WHg36Kvub3nTsjd+3H1+d/3SVpwr89X6j+WFIqFGz/taQeSQ/XKbk0IrptT5a0yvam7MjjC7LAWCZJJ3pi6zWjAIaJhj99sH2DpGsk/X7U6SiTdYxSROyUtFxSsbgGMOgaCgXb8yT9haRrI2JvnZpxtscfui3pKkmv1aoF0DoabTC7VNJ49Z0SrLd9b1Z7uMGspCmSnrP9iqQXJT0REU9W8iwANE2jDWbvq1N7uMFsRGyRdGGp0QEYdFzRCCBBKABIEAoAEoQCgAShACBBKABItP0U7091ry9U/9Wp+WsnFNoyWtlw/i5DURwpAEgQCgAShAKABKEAIEEoAEgQCgAShAKABKEAIEEoAEgQCgASbX+Z81enzi5U3zHpS7lrD77/QcHRYDj627fXFqq//fzfyF3bu7fmFKiV4kgBQIJQAJBotMHs922/m83kvN72/DqPnWf7Tdudtm9t5sABVKPRBrOSdFfWOHZ2RKwcuNJ2h6S7JV0taZakRbZnlRksgOo11GA2p7mSOiNiS0QckPSopAUNbAfAICrznsIS2xuy04uTa6yfJmlbv/td2bKabC+2vc72us+1v8SwAJTRaCjcI+lsSbMlbZd0Z42aWi2b6zaOjYhlETEnIuaM0pgGhwWgrIZCISJ2RMTBiOiV9EPVbhzbJen0fvenS+puZH8ABk+jDWZP63f366rdOHatpJm2Z9geLWmhpBWN7A/A4DnqFY1Zg9nLJE2y3SXpdkmX2Z6tvtOBrZJuymqnSvrXiJgfET22l0h6SlKHpPsjYmMVTwJA81TWYDa7v1LSFz6uBNC62v67D0XxfQY0281/98eF6kf9x3u5a0+a31l0OKVxmTOABKEAIEEoAEgQCgAShAKABKEAIEEoAEgQCgAShAKABKEAIDHsLnPe/9u/mrt2zBPFpu7GF3nOBblrY12tL9u2vknLni/2gGX5S4v8/iRpyzfG56o7sHRN3XUcKQBIEAoAEoQCgAShACBBKABIEAoAEoQCgESeiVvvl3SNpJ0RcUG27DFJ52UlEyR9HBGzazx2q6RPJB2U1BMRc5oyagCVyXPx0gOSlkp66NCCiPjdQ7dt3ylp1xEef3lEvN/oAAEMrjyzOT9j+6xa62xb0u9I+s0mjwvAECn7nsKvS9oREZvrrA9JT9t+yfbiI22IXpJAayj73YdFkh45wvpLI6Lb9mRJq2xvyrpYf0FELFN2VfiJnli352RZfJ9hcG2+ZVTu2nO+VeFA2tSIPcX+Bznx9RNy1W3fd4R9FtpjP7ZHSvqGpMfq1WTNYRQROyUtV+2ekwBaSJnTh9+StCkiumqttD3O9vhDtyVdpdo9JwG0kKOGQtZL8nlJ59nusv2dbNVCDTh1sD3V9qE2cVMkPWf7FUkvSnoiIp5s3tABVKHRXpKKiD+osexwL8mI2CLpwpLjAzDIuKIRQIJQAJAgFAAkCAUACUIBQIJQAJBwRGVXFDfM9nuS3hmweJKk4fBty+HwPHmOQ+/MiDil1oqWDIVabK8bDvMxDIfnyXNsbZw+AEgQCgAS7RQKBZpttbXh8Dx5ji2sbd5TADA42ulIAcAgIBQAJNoiFGzPs/2m7U7btw71eKpge6vtV22vt71uqMfTDLbvt73T9mv9lk20vcr25uznyUM5xmao8zy/b/vd7PVcb3v+UI6xiJYPBdsdku6WdLWkWZIW2Z41tKOqzOURMbtdP9+u4QFJ8wYsu1XS6oiYKWl1dr/dPaAvPk9Juit7PWdHxMoa61tSy4eC+uZ17IyILRFxQNKjkhYM8ZiQQzZJ74cDFi+Q9GB2+0FJXxvMMVWhzvNsW+0QCtMkbet3vytbdqzJPR1+m5sSEdslKfs5eYjHU6Ultjdkpxdtc5rUDqHgGsuOxc9RL42Ii9R3mvQ9218Z6gGhlHsknS1ptqTtku4c0tEU0A6h0CXp9H73p0vqHqKxVGYYTYe/w/ZpkpT93DnE46lEROyIiIMR0Svph2qj17MdQmGtpJm2Z9gerb5ZpFcM8ZiaaphNh79C0g3Z7RskPT6EY6nMoeDLfF1t9HqW7RBVuYjosb1E0lOSOiTdHxEbh3hYzTZF0vK+1pwaKelHx8J0+Fl7gMskTbLdJel2SXdI+nHWKuDnkq4buhE2R53neZnt2eo71d0q6aahGl9RXOYMINEOpw8ABhGhACBBKABIEAoAEoQCgAShACBBKABI/D/U1ce+UxOTRQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "classy = SVC(C = 50, kernel='rbf', class_weight='balanced')\n",
    "classy.fit(X_hbr_train, np.squeeze(Y_hbr_train))\n",
    "Y_hbr_pred = classy.predict(X_hbr_test)\n",
    "print(np.mean(Y_hbr_pred[:,None] == Y_hbr_test))\n",
    "plt.imshow(sklearn.metrics.confusion_matrix(Y_base_test, Y_hbr_pred, normalize='true'))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98708f2d-7bb7-4622-bf8e-901818fb8aae",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "stijn1",
   "language": "python",
   "name": "stijn1"
  },
  "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}