scMODD_with_combined_feature.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import seaborn\n",
    "import sys, os, random\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import scanpy as sc\n",
    "import statsmodels.api as sm\n",
    "from scipy import sparse\n",
    "import anndata as ad\n",
    "from scipy.stats import nbinom\n",
    "from scipy.special import gamma\n",
    "from sklearn.metrics import auc\n",
    "from sklearn.metrics import precision_recall_curve\n",
    "from sklearn.metrics import roc_curve\n",
    "from scipy.special import gamma, loggamma, gammaln\n",
    "import itertools\n",
    "import warnings\n",
    "from sklearn.metrics import roc_auc_score, precision_score, accuracy_score, recall_score, f1_score\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from scipy.optimize import curve_fit\n",
    "from tqdm import tqdm\n",
    "import numpy.matlib\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "plt.style.use('seaborn-white')\n",
    "from sklearn.ensemble import GradientBoostingClassifier\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.neural_network import MLPClassifier\n",
    "\n",
    "#R interface\n",
    "import rpy2.rinterface_lib.callbacks\n",
    "import logging\n",
    "\n",
    "from rpy2.robjects import pandas2ri\n",
    "import anndata2ri\n",
    "# import pandas as pd\n",
    "import rpy2.robjects as ro\n",
    "from rpy2.robjects.packages import importr\n",
    "# from rpy2.robjects import pandas2ri\n",
    "\n",
    "from rpy2.robjects.conversion import localconverter\n",
    "#Note: this can be commented out to get more verbose R output\n",
    "rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)\n",
    "import xgboost as xgb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The rpy2.ipython extension is already loaded. To reload it, use:\n",
      "  %reload_ext rpy2.ipython\n"
     ]
    }
   ],
   "source": [
    "# Automatically convert rpy2 outputs to pandas dataframes\n",
    "pandas2ri.activate()\n",
    "anndata2ri.activate()\n",
    "%load_ext rpy2.ipython\n",
    "\n",
    "plt.rcParams['figure.figsize']=(8,8) #rescale figures\n",
    "sc.settings.verbosity = 3\n",
    "#sc.set_figure_params(dpi=200, dpi_save=300)\n",
    "# sc.logging.print_versions()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "library(Seurat)\n",
    "library(scDblFinder)\n",
    "library(SingleCellExperiment)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def roc_pr_plot(y_label, y_pre ,path=None, title=None):\n",
    "    print(type(y_label), type(y_pre))\n",
    "    precision, recall, thresholds_pr = precision_recall_curve(y_label, y_pre)\n",
    "    fpr, tpr, thersholds_roc = roc_curve(y_label, y_pre, pos_label=1)\n",
    "    pr_auc = auc(recall, precision)\n",
    "    roc_auc = auc(fpr, tpr)\n",
    "\n",
    "    fig = plt.figure(figsize=(16, 8))\n",
    "    ax = fig.add_subplot(1, 2, 1)\n",
    "    ax.plot(fpr, tpr, 'r--', label='ROC (area = {0:.2f})'.format(roc_auc), lw=2)\n",
    "     \n",
    "    ax.set_xlim([0, 1.05])  \n",
    "    ax.set_ylim([0, 1.05])\n",
    "    ax.set_xlabel('False Positive Rate')\n",
    "    ax.set_ylabel('True Positive Rate')  \n",
    "    ax.set_title('ROC Curve - {}'.format(title))\n",
    "    ax.legend(loc=\"center right\")\n",
    "    ax.plot((0, 1), (0, 1), transform=ax.transAxes, ls='-',c='b', lw=0.5)\n",
    "\n",
    "    ax2 = fig.add_subplot(1, 2, 2)\n",
    "    ax2.plot(recall, precision, 'r--', label='PR (area = {0:.2f})'.format(pr_auc), lw=2)\n",
    "     \n",
    "    ax2.set_xlim([0, 1.05])  \n",
    "    ax2.set_ylim([0, 1.05])\n",
    "    ax2.set_xlabel('Recall')\n",
    "    ax2.set_ylabel('Precision') \n",
    "    ax2.set_title('PR Curve - {}'.format(title))\n",
    "    ax2.legend(loc=\"center right\")\n",
    "    ax2.plot((0, 1), (1, 0), transform=ax2.transAxes, ls='-',c='b', lw=0.5)\n",
    "    if path:\n",
    "        fig.savefig(path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NB model\n",
    "def func_fit(x, a):\n",
    "    return a* x\n",
    "def NB_pre(this_mtx, CURVE_FIT=True):\n",
    "    all_miu = np.mean(this_mtx, axis=0)\n",
    "    all_var = np.var(this_mtx, axis=0)\n",
    "    # curve fitting\n",
    "    X1 = (all_var-all_miu)\n",
    "    Y1 = np.power(all_miu, 2)\n",
    "    popt, _ = curve_fit(func_fit, X1, Y1)\n",
    "    all_theta = np.repeat(popt, this_mtx.shape[1])\n",
    "    return all_miu, all_var, all_theta\n",
    "\n",
    "def NB_model(y, theta, miu):\n",
    "    eps = 1e-10\n",
    "    #! log(0) -> nan\n",
    "    log_t1 = theta * (np.log(theta+eps) - np.log(theta+miu+eps))\n",
    "    log_t2 = y * (np.log(miu+eps) - np.log(miu+theta+eps))\n",
    "    log_t3 = gammaln(y+theta+eps) - gammaln(y+1+eps) - gammaln(theta+eps)\n",
    "    #! log(0) -> nan\n",
    "    results =log_t1 + log_t2 + log_t3\n",
    "    return results\n",
    "\n",
    "def generate_classify_param(score_lst, rank_lst=[1, 2]):\n",
    "    rank_a, rank_b = rank_lst\n",
    "    max_idxa, max_idxb = np.argsort(score_lst)[-rank_a], np.argsort(score_lst)[-rank_b]\n",
    "    max_a, max_b = np.sort(score_lst)[-rank_a], np.sort(score_lst)[-rank_b]\n",
    "    return max_a, max_b, max_idxa, max_idxb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ZINB model\n",
    "def func_fit(x, a):\n",
    "    return a* x\n",
    "def ZINB_pre(matrix):\n",
    "    this_pi = np.mean(matrix == 0, axis=0)\n",
    "    this_pi[this_pi == 0] = 1e-10\n",
    "    this_pi[this_pi == 1] = 1 - 1e-10\n",
    "    matrix = np.where(matrix, matrix, np.nan)\n",
    "    this_mean = np.nanmean(matrix, axis=0)\n",
    "    this_var = np.nanvar(matrix, axis=0)\n",
    "    filter_flag = ~np.isnan(this_mean) & ~np.isnan(this_var)\n",
    "    Y1 = np.power(this_mean[filter_flag], 2)\n",
    "    X1 = this_var[filter_flag] - this_mean[filter_flag]\n",
    "    popt, _ = curve_fit(func_fit, X1, Y1)\n",
    "    this_theta = np.repeat(popt, matrix.shape[1])\n",
    "    this_mean[np.isnan(this_mean)] = 1e-10\n",
    "    return this_mean, this_theta, this_pi\n",
    "\n",
    "def ZINB_model(y, theta, miu, pi):\n",
    "    this_nb = NB_model(y, theta, miu)\n",
    "    pi = numpy.matlib.repmat(pi, y.shape[0], 1)\n",
    "    this_nb[y==0] = np.log(np.exp(this_nb[y==0]) * (1-pi[y==0])+pi[y==0])\n",
    "    this_nb[y!=0] = np.log(1-pi[y!=0]) + (this_nb[y!=0])\n",
    "\n",
    "    return this_nb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def umap_preproc(matrix):\n",
    "    adata = ad.AnnData(matrix)\n",
    "    adata_filter_gene = ad.AnnData(matrix)\n",
    "    print('Before filter genes: ', adata.shape)\n",
    "    sc.pp.filter_genes(adata_filter_gene, min_cells=3)\n",
    "\n",
    "    adata.raw = adata.copy()\n",
    "    sc.pp.normalize_total(adata, target_sum=1e4)\n",
    "    sc.pp.log1p(adata)\n",
    "    sc.pp.highly_variable_genes(adata, n_top_genes=2000)\n",
    "    print(adata.var.highly_variable.shape)\n",
    "    idxs = adata.var.highly_variable\n",
    "    adata = adata[:, adata.var.highly_variable]\n",
    "    print('After filter genes: ', adata.shape)\n",
    "\n",
    "    adata_hvg = adata.raw.to_adata()[:, idxs]\n",
    "    sc.tl.pca(adata, svd_solver='arpack')\n",
    "\n",
    "    sc.pp.neighbors(adata)\n",
    "    sc.tl.umap(adata)\n",
    "    print(adata.X.shape, adata_hvg.X.shape, adata_filter_gene.X.shape)\n",
    "    return adata, adata_hvg, adata_filter_gene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_paths = []\n",
    "data_files = []\n",
    "data_dir = os.walk('../datasets/real_datasets/') # The directory storing all real-world data\n",
    "for root, dirs, files in data_dir:\n",
    "    for name in dirs:\n",
    "        data_paths.append(os.path.join(root, name))\n",
    "    for file in files:\n",
    "        if file.endswith('csv') and 'meta' not in file:\n",
    "            data_files.append(file)\n",
    "data_files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# scMODO with Combined Feature from scDblFinder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "REAL\n",
      "pbmc-1C-dm.csv --- meta_pbmc-1C-dm.csv\n",
      "Before filter genes:  (5270, 15865)\n",
      "filtered out 2893 genes that are detected in less than 3 cells\n",
      "normalizing counts per cell\n",
      "    finished (0:00:00)\n",
      "If you pass `n_top_genes`, all cutoffs are ignored.\n",
      "extracting highly variable genes\n",
      "    finished (0:00:01)\n",
      "--> added\n",
      "    'highly_variable', boolean vector (adata.var)\n",
      "    'means', float vector (adata.var)\n",
      "    'dispersions', float vector (adata.var)\n",
      "    'dispersions_norm', float vector (adata.var)\n",
      "(15865,)\n",
      "After filter genes:  (5270, 2000)\n",
      "computing PCA\n",
      "    on highly variable genes\n",
      "    with n_comps=50\n",
      "    finished (0:00:00)\n",
      "computing neighbors\n",
      "    using 'X_pca' with n_pcs = 50\n",
      "    finished: added to `.uns['neighbors']`\n",
      "    `.obsp['distances']`, distances for each pair of neighbors\n",
      "    `.obsp['connectivities']`, weighted adjacency matrix (0:00:00)\n",
      "computing UMAP\n",
      "    finished: added\n",
      "    'X_umap', UMAP coordinates (adata.obsm) (0:00:07)\n",
      "(5270, 2000) (5270, 2000) (5270, 12972)\n",
      "running Leiden clustering\n",
      "    finished: found 7 clusters and added\n",
      "    'leiden', the cluster labels (adata.obs, categorical) (0:00:00)\n",
      "(7,)\n",
      "Now Running:  (5270, 2000)\n",
      "<class 'list'> <class 'numpy.ndarray'>\n",
      "('FINISH!!!---', 0.8)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6kAAAHoCAYAAACrRKlTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAChuElEQVR4nOzddXiTVxsG8Ds1SgXX4b5hw4a7DR3DZYNiw10rFCm0FNcBg8FwX2HY2HA2GAwfMmC4e5EK1Xx/PF+blra0hSQncv+uq9cba3JToG+enHOeo9FqtVoQERERERERmQAb1QGIiIiIiIiIYrBIJSIiIiIiIpPBIpWIiIiIiIhMBotUIiIiIiIiMhksUomIiIiIiMhksEglIiIiIiIik2GnOgCZr2LFiqFo0aKwsbGBRqNBaGgoXFxcMH78eJQqVQoAEBISgnnz5mH//v1wcHAAANStWxd9+/aFo6Nj7HNt2bIF69evx9u3bxEREYHy5ctj5MiRSJcuXaKvndrHm6u7d+9i6tSpmDdvnt6f+/jx45g4cSJ27Nih9+dOysOHD9GuXTv88ssvyJQpU+ztBw4cwLJly/D69WtERkaiSJEiGD16NHLmzJnscy5duhT//fcf/P39DRmdiMgs3bt3Dw0aNEDRokVjb9NqtejSpQvatGmT6P0hISHIkSMH/Pz8kCdPngTPGRUVhZUrV2L79u2IiopCREQE6tSpg8GDB8ee6y3FP//8g82bN8PHx0fvzx0QEIDffvsNP/zwg96fOymXL19Gz5498eeff8a7/WPeV/n4+CBjxowYOHCgoWKTFWKRSh9lxYoV8YqNpUuXYtKkSdiwYQMiIyPRrVs3lClTBlu3bkXatGkRGhqKGTNmoEePHlixYgXs7OywaNEiHD58GN9//z2yZMmCiIgI+Pn5oU+fPli7dm2C10zt483ZgwcPcPPmTdUx9GLr1q2YO3cunjx5Eu/27du3Y+HChVi4cCHy5csHrVaLxYsXo0uXLti5c6fFveEhIjI2R0dH/PLLL7HXHz9+jGbNmqFkyZJwcXFJcL9Wq8WkSZMwa9YszJw5M8HzjR8/Hq9evcKKFSvg6uqKkJAQjBgxAl5eXpg2bZpR/kzGcu3aNTx+/Fh1jI8WGRmJ1atXY8mSJQgJCYl3nzW9ryLzwSKV9CYyMhIPHz5E+vTpAQC7d+9GdHQ0PDw8Yh+TNm1aeHl54euvv8aePXtQq1Yt/PDDD9iyZQuyZMkCALC3t8eoUaOwZ88ehIeHxytSQkJCkn38Dz/8gMDAQIwdOxYAMG/evNjrnTt3Rvr06XHjxg20b98eCxYswB9//AEHBwdERUWhdu3aWL58ObJlywZfX19cvXoVERERqFKlCkaNGgU7u5T/lzl+/DhmzpyJnDlz4ubNm0ibNi169eqFVatW4ebNm2jYsCE8PT0BAPv378fChQsREREBR0dHjB49GqVLl8aYMWPw+PFj9OjRA0uXLsWiRYuwb98+vH37FqGhoRg9ejQaNGjw3hx169ZF06ZNceTIEbx58wbdunVDp06dYn+egwYNwu3bt5EuXTr4+PigQIECcHd3h6OjI65evYrnz5+jbt26yJAhAw4cOICnT59i0qRJqFKlCoKDgzFp0iScPn0atra2qF+/PoYOHQqNRhMvw+PHj7F3714sXboUjRo1inffrFmzMHHiROTLlw8AoNFo0KtXL+TMmTPB3z8AREREYNKkSTh69CgyZ86MzJkzw9XVFQDQuXNnlChRAmfPnsWLFy/Qrl07PHv2DH///TdCQ0Mxe/ZsFCtWLMHPKObflJ2dHfLlywd/f3+8ffsWo0ePRmBgIACgVq1aGDJkCDp06IBu3brhyy+/BIDYN2QjR45M/h8FEZEJyJ49O/Lly4dbt26hZMmSCe4PCwvDkydPYs+zcd27dw/bt2/Hn3/+CRcXFwCAk5MTJkyYgNOnTwMA3N3dUaRIEfTo0SPB9bp166J06dK4cuUKBg4ciIULF2L79u0AgNevX6NevXrYu3cv3r59Cx8fHzx8+BARERFo2rQp+vTpk6o/Z0BAAH7//XdER0fjwYMHyJ49O9q1a4fVq1fj1q1b6NatG7p37w4A2LRpE9atW4fo6GhkyJAB3t7ecHJywty5c/HmzRt4eHjA19cXfn5+OHfuHIKDg2OL+fLly783R/HixfHdd9/hjz/+QEhICIYNG4aGDRsCAJ4+fYoePXrgyZMnyJUrFyZOnIisWbOm+Hz29OlTjBs3Djdu3ICNjQ06dOiALl26JMhw6dIlXLlyBfPnz4/9MwMpe1/17nk4KCgIXl5euHz5MrJlywZbW9vYn0HdunXRrFkzHDt2DK9evULPnj1x+vRpXLx4EXZ2dli4cCGyZ88e7/kiIyMxbdo0HDx4ELa2tihbtizGjRuHu3fvwsvLC+Hh4dBqtWjTpg06dOiAunXr4vvvv4/9tztkyBBUrFgx9r0NWQauSaWP4ubmhubNm6N69eqxb9onT54MADhz5gwqVKiQ4Hs0Gg2qVKmCU6dO4caNG3B0dET+/PnjPSZt2rT46quvEvxiTO3jE5MuXTrs2rULbm5uKFKkCPbv3w8A+PPPP5E7d24UKlQIfn5+KFGiBAICArB161YEBgbip59+SumPJdb58+fRq1cv/PLLL3BxccHixYvxww8/ICAgAGvXrsXjx49x69YtzJo1C4sXL8bWrVsxceJEDBw4EGFhYZg0aRLy5s2LpUuX4v79+zh69ChWrVqF7du3Y+jQoZg7d26Kcrx69Qo///wzVq1ahblz5+LKlSsAZPpt165d8csvv6BZs2YYNWpU7PdcunQJK1aswOrVq7Fs2TI4OTlh/fr16NKlC5YsWQIAmDt3LsLCwrBr1y5s3boVp0+fxt9//53g9bNnz4758+ejQIEC8W4PDAzE/fv3Ua5cuXi3azQafPXVV7FvgOJau3Ytbt26hZ07d2LZsmV4+PBhvPvv37+P9evXY9q0aZg2bRoqVqyIgIAA1KhRA6tXr07wfPv27UNAQAA2bNiAHTt2IHfu3Fi9ejU2btyI3LlzY8uWLVizZg1u376NN2/eoG3btggICAAgU962bduGtm3bpuSvgYjIJJw5cwZ37tzB559/DgB4+/YtWrRogebNm6Nq1apo2bIlChYsiBEjRiT43osXL6Jw4cIJfj9nzZo19n1AcooUKYJff/0VjRs3RnBwMM6fPw8A2LFjB2rVqoX06dNj5MiRaN26NQICArB582YcPXoUu3btSvWf9eTJk5gwYQK2bduGR48eYefOnVixYgWWLFmC2bNnIzo6Gn///Te2bt2KNWvWYOvWrejZsycGDBiAnDlzYtCgQahQoQImT56Mc+fO4cmTJ9iwYQN27dqFli1bxp4P3ycqKgpp06ZFQEAAZs+eDU9PT7x48QIAcPPmTYwdOxbbt29H0aJF4evrG/t9KTmfTZgwAfnz58fu3buxYcMGbNy4Ebdv306QoXTp0pg8eTKyZs0a7/YPeV81d+5cODo6Yvfu3ZgzZ06CGV9hYWHYuHEjBg8ejLFjx8LNzQ3btm1Dzpw5sWXLlgTPt3btWly8eBG//PILduzYgeDgYOzatQtLly5F3bp1ERAQgMWLF+PkyZPQaDSx/y4AeX/z119/oXnz5sn+PZB54UgqfZSY6b4XL15Er169UKlSJWTOnDn2/sjIyES/Lzw8HLa2trCxsUF0dHSKXy+1j09M3MK5TZs22LJlCxo1aoSAgAC0a9cOAHDw4EGcP38emzdvBiAn8A+RO3duFC9eHACQN29euLq6wsHBAZkyZYKzszNevXqFEydO4MmTJ+jatWvs92k0Gty5cyfec+XKlQtTp07F9u3bcfv27dhPclOiU6dO0Gg0yJEjB2rUqIEjR46gRIkSKFasWGyB2LJlS4wfPx5v3rwBANSpUwf29vbImjUrnJycUKNGjdg/x8uXLwEAR48ehYeHB2xtbWFra5toEfg+NjbyOVlq/k7/+usvNGvWDA4ODnBwcEDz5s1ji24AsSPLMeuo4uZOrID+66+/0KhRo9gZADEj///88w969eqFhw8fomrVqhg+fDhcXV3RpEkTTJ06FU+fPsWlS5eQP3/+BCd3IiJTElOEAlIwZcyYEdOmTUPOnDlx7969eNN9//jjD4wcORJ16tSBs7NzgufS53k4puDYsmULSpUqhYCAAIwaNQohISE4ceIEXr16hTlz5gCQEb/Lly+jSZMmqXqtUqVKxfY3yJ07N6pXrw4bGxvkyZMHYWFhCA0NxcGDB3H79m106NAh9vtev34de66LUbZsWaRPnx7r16/H3bt3cfz48UR/Ron59ttvAQCffvopihYtihMnTgAAqlatGjuTqE2bNmjTpk3s96TkfHb06NHYmTyurq6p7jPxIX+ff/31Fzw9PaHRaJApU6YEM7piRonz5MmDLFmy4NNPP43N/erVqwTPd/ToUbRo0SK2V8ns2bMBAOnTp8fo0aPxzz//oEqVKhgzZgxsbGzQunVrtGnTBu7u7tixYwfq1q0bO6OKLAeLVNKLEiVKwMPDA+7u7vjss8+QO3dulCtXDj/++COio6NjixFACpITJ06gb9++KFy4MCIjI3Hr1q14b/TDwsIwYMAATJo0Kd60kJQ8XqPRQKvVxt4XERERL6uTk1Ps5caNG8Pf3x/Xr1/HiRMnYpvvREdHY86cOShUqBAAOVm9O4UVQOxJHwAmTZoU2zAqxrufQCY2XTg6OhpVqlSJ/aUMyAhntmzZcPLkydjbLl68iH79+qFr166oVq0avvjiC0yYMCHB8yUm7uvG/fuI+/cCyBuGmMemJLudnV28n8vDhw/h6OgYr+BO7OcSI3369MifPz/OnTuHqlWrxrtv8ODB6Nu3L1atWoULFy4AQLw3EDFsbW3jXX83t729faKvHff74/4ZXr9+jdevX6N06dLYt28f/vrrLxw7dgxt27bFkiVLULJkSXz55ZfYsWMHzpw5w1FUIjJ57645fZ8aNWqgW7duGDx4MHbu3JlgxLR06dK4ceMGgoKC4t33+PFjeHt7Y+7cuak6D7dp0wYtW7ZE27Zt8ebNG1SsWBFBQUHQarVYv3490qZNCwB48eIF0qRJE+95Hj9+jF69esVeX7x4cYKppCk9D7do0SK22IuOjsaTJ09iP7yMcfDgQfj6+qJbt26oV68eChYsiG3btiV4vsTEPVdFR0fHXn/39rj5UnI+e/c8fPfuXdja2qJv376xtyX2c4mRkvdVY8aMie0nMWjQIACI9/f7vvNwcufgmD9DXM+ePUN0dDTq1KmD3377DUePHsVff/2F77//HgEBAciVKxeKFy+OgwcPIiAgIHbpFFkWTvclvWnWrFnsdBIA+PLLL5E2bVr4+fnFjkS+ffsWEydOhLOzMxo0aAAHBwd899138PLywrNnzwDIKKufnx9CQ0MTPdkk9/iMGTPi4sWL0Gq1CAoKwoEDB5LMnCZNGjRt2hTu7u5o2LBh7MmwevXqWL58ObRaLcLDw9G3b99ERwl/+eWX2K+kCrHkVKlSBUeOHMH169cBAIcOHcJXX32Ft2/fwtbWNvbkfuLECZQsWRLdunVDxYoVsW/fPkRFRaXoNbZu3QpAGjEdOXIENWvWBABcuXIF//77LwBgw4YNKF++fOzPIKXZt2zZgujoaISHh2PQoEE4ceJEqn4uAwYMgK+vb+z0pKioKCxYsACXL19GwYIF4evrG/tcHTt2RI0aNbB161aEhYXFTjX+GFWrVsWePXsQFBQEQNYwL1++HNOnT8eCBQtQv359eHl5oXDhwvjvv/8AAO3atcOWLVtw+vTpFE9vIyIyF927d4ezs3OiS0qyZ8+O5s2bw9PTM/b3ZlBQEMaPH48MGTLA0dERGTNmjP1w8fHjx4nOYon7fKVLl8bYsWNjRxFdXFxQpkyZ2GU2r1+/RseOHbFv374E3xv3fJNUIZac6tWrY+fOnbGF2Lp16+Dm5gZACrCYWWFHjhxBnTp10KlTJ5QsWRJ79+5N9Xn44sWLuHnzJr744gsA0r/iwYMHAID169fHnp9TqkqVKvj5558BAG/evIGbmxtevHiR4p9LSt5XLVmyJPa56tWrhxo1amDz5s2Ijo7Gq1evEvy9pFaVKlWwY8cOhIeHIzo6GuPHj8fOnTsxfPhw7Nq1C02bNsW4cePg4uISO8usXbt2WLJkCUJDQ5NdE0zmiSOppFfe3t746quv8Mcff6BGjRpYtmwZFixYgFatWsHGxgZRUVGoW7culi1bFvvpWp8+fZA2bdrYBgthYWGoWLEiFixYkOhrJPf4mNdv2LAhsmfPjooVK8b7xO9dbdu2xerVqzF+/PjY27y8vODr64vmzZsjIiICVatWRc+ePfXxI0qgcOHC8PHxwbBhw6DVamMbCzg7O6Nw4cJIkyYN2rRpg0WLFuH3339H48aNYz9hfPXqVYJPsxNz7949tGrVCm/fvsWYMWNQsGBBPH36FAULFsT8+fNx9+5dZM6cOdXbuMQUmC1atEBUVBSaNGkSO80npZo3bw6tVothw4YhMjISYWFhKFGiBFasWJHoWpgOHTrgzp07aNasGTJkyBA7TSo19u3bh/Xr12PJkiWoVasWrl27ho4dOwKQv4+JEyciNDQU7u7usVOLixUrhqZNmwIASpYsCVtbWzRq1CjBJ/tERObO3t4e3t7e6NmzJ9q0aRNvexoAGDduHBYsWIAOHTrA1tYW4eHhqF+/fuwWJJ07d8aIESPw5ZdfInfu3KhcufJ7X69t27YYPHgwFi5cGHvb9OnTMXHiRDRv3hzh4eFo1qwZvvrqK/3/YSFF6nfffYfu3btDo9HAxcUF8+fPh0ajQZkyZfD9999jwIABGDp0KIYPH47mzZsjMjIS1apVi23M9O7MpHedPn0aGzduRHR0NGbNmhU7Slu0aFF4enri2bNnKFiwYKq3uhk7dizGjx8fey7t3bt3os2w3ie178MGDhyIcePGoXHjxsiUKVOCfx8psW7dOly4cAG+vr7o0KED7t+/j1atWkGr1aJixYro3Lkzbt++DS8vL2zYsCG2OWNMcV+3bl1MmDAB3333Xapfm8yDRvu+d+9EZPbq1q2LOXPmfPBILxEREX24YsWK4a+//oq3ZR8RvR9HUonM3LZt27B06dJE72O3OyIiIsP68ccfY7fReVfM6CQRpQ5HUomIiIiIiMhksHESERERERERmQyTnO779u1bXLhwAVmzZk3Q1pqIiCi1oqKi8PTpU5QsWTJ2Lz5KHZ6biYhIn953bjbJIvXChQv45ptvVMcgIiILs2bNGlSoUEF1DLPEczMRERlCYudmkyxSs2bNCkAC58iRQ3EaIiIyV8HBwJQpgJvbI/Tq9U3s+YVSL+ZnN23aGpQrx3MzERF9nEePHuGbbxI/N5tkkRozjShHjhzInTu34jRERGSO3r4F/P2BqVPlMgBOU/0IMT+77dtzoEKF3PjkE8WBiIjIIiR2bmbjJCIisjiRkYCHB+DuDmTJojqNZRk9GvD1BZ4/V52EiIgsFYtUIiKyKNHRgJcX0K8fwMk4+ufoCPj5Ad7ewJs3qtMQEZElYpFKREQWQ6sFfHyATp2AIkVUp7Fc6dMDEybISHXMVGoiIiJ9YZFKREQWY/p0oGFD4PPPVSexfFmzSpHq7i7Tq4mIiPSFRSoREVmEhQulOK1aVXUS65Enj0yr9vKSadZERET6wCKViIjM3qpVQLZsMopKxlW0KNCxo0yz1mpVpyEiIkvAIpWIiMza1q0yite6teok1qtMGaBBA2DGDNVJiIjIErBIJSIis7V3L3DvHuDmpjoJVasGlCwJLFqkOgkREZk7FqlERGSWjh8HTpwABgxQnYRiNGok+9KuXq06CRERmTMWqUREZHYuXAB27pTOsmRa2rQBIiKAX35RnYSIiMwVi1QiIjIrN24AK1YA48cDGo3qNJSYbt2A27eBfftUJyEiInPEIpWIiMzGgwfA3LmAnx9gwzOYSRs0SKZkHz+uOgkREZkbnuKJiMgsPH8uxenkyYC9veo0lBIeHsC2bTI9m4iIKKVYpBIRkcl78wYYO1aK1LRpVaehlNJoZP/U5ctlmjYREVFKsEglIiKT9vYt4OkJTJgApEunOg2llq2tjH7Png08fKg6DRERmQODFannzp1D586dE9y+f/9+tG7dGu3bt8fGjRsN9fJERGQBIiNlyujo0bK1CX0cVedme3vA3x+YNAl48ULvT09ERBbGzhBPumTJEmzbtg1p35mTFRERgcmTJ2Pz5s1ImzYtOnbsiDp16iBr1qyGiEFERGYsOhrw8gL69QNy51adxvypPjc7Ocl0bU9PYMoUwMVFr09PREQWxCAjqXnz5sW8efMS3H79+nXkzZsX6dOnh4ODA8qXL4+TJ08aIgIREZmTEiWAAgViv7T5C8BnwBN06gQUKaI6nGUwhXNz+vSydZC7u0zjJiIiSoxBRlK//PJL3Lt3L8HtQUFBcHV1jb3u7OyMoKAgQ0QgIiJjevECmDVLOhw1aQI0bCi3//MPsGxZ0t83aZIMqd2+DQQHx948HSPQoEIgPv88m4GDWw9TOTdnzSrTt93dgenTATuDvBMhIiJzZtRTg4uLC4LjvAkJDg6Od2IkIiIzcuKErrDcsAFYtEguZ8miK1KvXwfmzEn6Oby8pEi9cAHQagEAi9a4onSeSFTrkNGA4SmGinNznjwyjdvLS5oqcc9bIiKKy6hFaqFChXD79m28fPkSTk5OOHnyJHr06GHMCEREpC89egDnz8e/rVs3oEED3fVSpWSENSnOznLMnx8AsHo1kPUz4MvW+o1KSVN1bi5aFOjQQbaoGTdOtqshIiICjFSkbt++HSEhIWjfvj3c3d3Ro0cPaLVatG7dGtmzZzdGBCIi0rcKFYBMmXTXs2eXUdO4o3CFCwNDhqTo6X75Rbr5tmaBahSmcG4uW1YG46dPB0aONMpLEhGRGTBYkZo7d+7YNvbNmzePvb1u3bqoW7euoV6WiIgM6cQJYPNmGf5631rTVNq3D7h7FxgwQG9PSYkwxXNz9eqylHnRIqBPHyURiIjIxHAVCBERpVyLFsDUqYCjI3D1ql6e8vhxqX1ZoFqvxo2BzJmBNWtUJyEiIlPAIpWIiBL6fxOj2MvR0cDFi8DDh3LbokWyqPAjXbgA7Nwp3V7JurVtC4SHy7RvIiKybixSiYgoPl9fYMoU3fWAAMDWFihZUndbr14f/TI3bgArVsi+mWyaQ4D03bp9G9i/X3USIiJSibuTERFZswEDgDt34t/2xx9A1arSevX/XXdjq0hbWyliP7KqfPBAeixNn87tRyi+QYMAPz/ZmahiRdVpiIhIBRapRESWJjIS2LIFeP064X0vXsjWMTFdeQ8elGm87zp6FMiRQy63bi3TffXkxQspQqZNA+zt9fa0ZEE8PABvb8DJKf4APhERWQcWqUREliA8XFqkurgAjx4BO3YA27cDgYEJH+vnB6xcCTRvDsybJ9/3rhIlpDmSngUFSfExeTKQNq3en54shEYDTJgga5X79QMKFlSdiIiIjIlFKhGRuXvwAChTBnj6FNi0CWjTBpg9W9qlvnqV8PEZMgC1a8vlOnWMFvPtWxkhGz8eSJfOaC9LZsrWVj5PGTFC/t3kzKk6ERERGQuLVCIic7VhA/D991KkPn0qcyMdHOS+jBmBmTPV5osjMlIKjVGjgKxZVachc+HgAPj7AyNHAhMn6mapExGRZWORSkRkju7fl7WlwcFyPU0a4PRpoFgxtbkSER0NjBkj0zbz5FGdhsyNk5OMqHp4yBa9Li6qExERkaGxSCUiMiVBQbq9SBNToABgZwf4+EiB2rQp4O4O5M6t68RrQrRaGQHr2BEoUkR1GjJX6dPLNHF3d2DGDPlMhoiILBeLVCIiUxEcLB1inj5N+jGPHwPZsgHdu0uDpKlTTXL0NMaMGUD9+sDnn6tOQuYuWzaZLu7uLp2h7fgOhojIYvFXPBGRSiEhssVLnjzA4sXA7dtAr17AsWOJP97WVo6VKgG//GK8nB9g0SKgVCmgWjXVSchS5M0L9O0r08f9/LjHLhGRpWKRSkSkwsOHwJIlwM6dwN9/y22LF8u+LKtWqc2mB6tXA1myAF9+qToJWZqiRYH27WUa+dixsl0NERFZFhapRESG8uoVcOWKXM6ZU9c16OVLoFWr+KOl1asbPZ6h/PILEBEBfPut6iRkqcqWldnxM2bIFjVERGRZOFGGiMgQtFrZu7RSJflatEh337FjugK1cmVZV7p2rZKY+rZ/P3DnDtCtm+okZOmqVwdKlIj/X4uIiCwDi1QiIn3q2FFGTdOkAW7dktu++EK678ZIn15u+/prYO9e2QTSAvZm+ftv4PhxYOBA1UnIWjRuLHunrlmjOgkREekTi1QiIn0KDJSuuxERcr1tW6ne+vbVPaZKFbltyxbA2VlNTj27cAHYvl06rxIZU7t2QFiYyfcRIyKiVOCaVCKilHrzBjhyJOHtP/0ENGgAtGkj03bDwuR2jQbInt24GRW4cQNYsQKYMoVNbEiN7t2BOXNkunnduqrTEBHRx2KRSkSUUrduyfzCxGzcKPflymXUSKo9fCjFwfTp3A6E1Bo8GPD1BVxcgIoVVachIqKPwSKViCilXFyS3lOlYkWrK1BfvJCiYNo0wN5edRoiwNNT9lB1dpamSkREZJ5YpBIRxTh0CFi4MOHtu3cDo0cDbm5ymRAUBHh7A5Mny9auRKZAowF8fOS/a//+QIECqhMREdGHYJFKRNbh1Svg99+ByMiE92XPLgvZbt4ENmxI/Pu9vKRDC+HtW2mQNH48kC6d6jRE8dnaAn5+sn+qh4c02yYiIvPCIpWILFtUlLxrHTpUGhwlpl49KVJr1QLWrUv8MQUKAIUKGS6nmYiMlDf+o0cDWbOqTkOUOAcHwN9fdneaOFG2qSEiIvPBIpWILNegQUD+/FKgli4NFCwo+5d+/nn8x5UsKccCBTg/8D2io2W9X9++FrGtK1k4JycZUfXwAKZOlSXlRERkHlikEpHlmjcPKFVKmh0NGSJf9EG0WhmR6tABKFpUdRqilEmfXqalu7sDM2bIZ1RERGT6uGEAEVmmR4/keP4856XqwcyZQP36QJkyqpMQpU62bMCoUVKoJrYknYiITA+LVCKyLKGhwP37MooKyOad2bKpzWTmfvhBtvOoVk11EqIPkzcv0KeP9D+LjladhoiIksMilYgsy2+/Ablzy2I0APjmG7V5zNyaNUDmzECjRqqTEH2cYsVkuvrEiTJ9nYiITBeLVCKyLI6OwCefyFeZMtLikz7Itm1ARATQpo3qJET6UbasNPKeOVN1EiIieh82TiIi83LoEPDHH4nf5+oKDB4s033po+zfD9y+DQwcqDoJkX7VqAEEBck09t69VachIqLEsEglItNx/bqsKY3r1CkgJESGP4oVk+rJxyfx7//kEylS6aP8/Tdw7Bjg6ak6CZFhNG4MbNwIrF0LdOqkOg0REb2LRSoRqaHVAocPAzVrAhqN3NahA3DyZOKPnztXitSaNaX7SWJcXQ2T1YpcuCDTfCdOVJ2EyLDatQOWLZN/7199pToNERHFxSKViIyvRw8ZEb11S4Y0du6UQrVgwYQjqQCQMSPQpIlcrldPvkjvbt4Eli8HpkzRfW5AZMm6d5fPvw4cAOrUUZ2GiIhisEglIuN59Uq67i5bprvtyy91FdGGDWpyER4+BGbPBqZPB2xtVachMp5BgwBfX8DZGahYUXUaIiICWKQSkbFotdImdu9eue7oCFy7JutISakXL+RN+tSpgL296jRExufpKasInJ1lT2AiIlKLW9AQkXFERQEODjJ1d9o06c6TKxfnlSoWFAR4e0uR6uSkOg2RGhqNrMNetkymvRMRkVosUolIP7JkATJlSvxr2TLAzk7m0m3dCowYAZQqpTqx1QsLAzw8gPHjgfTpVachUsvWFpg8WfZQffhQdRoiIuvG6b5ElHpz5gBnzwK1agFdu8ptgYFAdHTijw8Lk+O4ccZIRykQGSkF6qhRQNasqtMQmQYHB2kcNnKkjKxmyqQ6ERGRdWKRSkQiPBw4dEj2JAVkX9KYLV1OnQLu3ZPLt28DQ4bIZQcHXZH69GnSz815pCYlOhoYMwbo0wfIk0d1GiLT4uQk09+9vKRgdXFRnYiIyPqwSCUiMWsW4O6uu37pEvDZZ7r71qxJ+D1ubrrLHHIwC1qtjBB16AAULao6DZFpypBBJn54eEjH6zRpVCciIrIuLFKJSPj76y5/9VX84YNy5YA3b3TXbW2BAQOAqlWNl4/0YsYM2Wa2TBnVSYhMW7ZsMu3X3V16vdnxHRMRkdHwVy4RSeXy8qVc/vlnoFWr+PcPGyZfZNZ++EG216heXXUSIvOQN69Mix8zRrZ4tmG7SSIio+CvWyJr9+aNDBcA0uK1SRO1ecgg1q6VGdmNG6tOQmReihUD2rWTafJareo0RETWgUUqkbVzdZWmSRs2ACdPAo6OqhORnm3bJg2W27ZVnYTIPJUrB9SpI9vTEBGR4bFIJSJZbNWuHVC4sOokpGcHDkhD5m7dVCchMm81a0ovuR9+UJ2EiMjycU0qkTU4elRaVYaHJ7wvKEh2sK9alXstWJi//wb++gvw9FSdhMgyNGkCbNwo0+c7dVKdhojIcrFIJbIk4eHAihXxO/F27ixtKr/+WpofJVaoNmsG3L/PItWCXLwo03wnTlSdhMiytGsHLF0KbN8ONG+uOg0RkWVikUpkKbRaaUP500/xb2/QAChVCnBwAEqWTPx7CxQAsmY1fEYyips35Z/BlCmARqM6DZHl6dEDmDNHptPXqaM6DRGR5WGRSmQpRo6UykSjAfr3B+zt5fYsWeSYN698kUV7+BCYPRuYPl22syUiwxg8GPD1lQkoX3yhOg0RkWVhkUpkKY4ckWPLlsC8eWqzkBIvXgCTJgHTpuk+oyAiw/H0lD1UnZ2B4sVVpyEishwsUolMxb17wPLlslcIAPTsCeTLJ5d37gSOHUv4Pf/9B3z3HVCjBuDtDbRqBcyfb7TIZDqCguSfgJ8f4OSkOg2RddBoAB8fYNQoYMAAWTlBREQfj0UqkamYPBlYsEB3vVEjXZH6229Jj45u3ixDaDVqSKEbM72XrEZYGODhIQ2c06dXnYbIutjayq/vESPk/2HOnKoTERGZPxapRKYiKEiOhQoBbm7x1482aZJ0Y6OyZYF06eSyq6thM5LJiYyUN8YjR0oTZyIyPgcHwN9f/h9OnAhkyqQ6ERGReWORSmQKfHyATZvksre3FKlxNWokX0RxREfLP5c+fdgTi0g1JydppOTpCUydyh29iIg+ho3qAEQE4NUrIDRUPo4vUUJ1GjIDWq00SWrXDihaVHUaIgKADBlk2r27u669ABERpR6LVCIV7twB9uwB9u3TDYfdvQs8eQJUqKA6HZmBmTOBunVltjcRmY7s2aWRkru7TMcnIqLUY5FKZGzBwUDJkkDDhkD9+rKe9O1bIHdudr2hFPnhB9nuonp11UmIKDF58wK9e8vnj1qt6jREROaHRSqRMQ0ZIguV3ryR6/XrAwMHAjlyKI1F5mPtWmnK0rix6iRE9D6ffgq0bSuNlFioEhGlDhsnERlT3L1O+/QBFi5Ul4XMzrZtMujeqZPqJESUEuXKSeP2WbOAYcNUpyEiMh8cSSUypg0b5PjwIQtUSpUDB4Bbt4Du3VUnIaLUqFlTRlUXL1adhIjIfLBIJTKGsDD5ypkTOHOG03spVU6cAI4eBQYNUp2EiD5EkybS+XfdOtVJiIjMA4tUIkNr1w5wdAQ++0w+Si9TRnUiMiOXLgG//CJ7LxKR+WrXDggJAbZvV52EiMj0sUglMgQvL+ls07gxsGmT3HbzpgyJEaXQzZvAsmXAhAmARqM6DRF9rB49gOvXZfo+EREljY2TiPThyBHg/Hmgdm1ZfHTihOyDGiNdOuDxYxlRJUqBhw+B2bOBadMAW1vVaYhIXwYPBnx9pdH7F1+oTkNEZJpYpBJ9iMBAWWMKAK9f6zasXLpUitRJk4ChQ3WPL16cBSql2IsX8iZ26lTAwUF1GiLSJ41GJtt4eQHOznJ6ICKi+FikEn2Ib78Fdu1KeHuxYnKsWNG4echiBAUB3t6Anx/g5KQ6DREZgkYj+6eOGiVbZefPrzoREZFp4ZpUouRotbKQqFIl4ORJuS1DBiB79vhf7u5AtWpKo5J5CwsDPDyAceOA9OlVpyEiQ7K1BSZPBmbOBB49Up2GiMi0cCSVKDFLlgDPnsnlGzekew0gU3sBYM0aNbnIYkVGSoE6ciSQLZvqNERkDA4OgL+//L+fNAnImFF1IiIi08AilayPVgv8959UBe8qUABIm1Y61ly6lPD+8uUNHo+sT3Q0MGYM0Ls3kDev6jREZExOTrIG3ctL1qE7O6tORESkHotUsj6jR0vL1MScOAFUqAD07Ak8eaK7XaMBvv6aczBJ77RaGUFp3163pJmIrEuGDMDYsbJqZPp0IE0a1YmIiNRikUqWY/x4YO3axO/LlAk4dkwuR0Xpbv/ss/iPi+nAG7czL5EBzZoF1KkDlC2rOgkRqZQ9OzBihBSq06YBdnyHRkRWjL8CyXJcuSLTeBOTNavu8tSpwNatwM6dsl0MkSKLF8vnJDVqqE5CRKYgXz6Z9j9mjDRV0mhUJyIiUoNFKpm/58+lOF23TvbtiIhI+Bhb2/iXz5wB0qUzXkaid6xbJ01SGjdWnYSITMmnnwJt2wI+PjIFmIUqEVkjFqlk3qKjgdKlgRcvgNBQaXyUEixQSaHt2+Wfa/fuqpMQkSkqXx4IDpbtaYYPV52GiMj4uE8qmadz54DChaWR0YMHwNu3qhMRpciBA8DNmyxQiej9ataUUdXFi1UnISIyPoMUqdHR0Rg7dizat2+Pzp074/bt2/Hu37ZtG1q2bInWrVtjbVKNbojep2NH4Pp1IChIrtesqTYPUQqcOAH89RcwaJDqJGSNeG42P02bymex69apTkJEZFwGKVL37t2L8PBwbNiwAcOHD4e/v3+8+6dOnYqffvoJ69atw08//YRXr14ZIgZZqjdvgH//lctDhgBPn8rwFJEJu3RJ+nV5eKhOQtaK52bz1L49EBIC7NihOgkRkfEYpEg9deoUavy/XWWZMmVw4cKFePcXK1YMb968QXh4OLRaLTTsCkCp4eQEXLggPfonTQKyZAFsOHOdTNfNm8DSpdIIhb/uSBWem81Xjx7AtWvAwYOqkxARGYdBGicFBQXBxcUl9rqtrS0iIyNh9/9Nv4oUKYLWrVsjbdq0aNCgAdKxiQ2l1J9/Art3S3FaooTqNETJevRI9kKdPj1+k2kiY+O52bwNHgz4+gIuLkCFCqrTEBEZlkGGn1xcXBAcHBx7PTo6OvYkePnyZRw8eBD79u3D/v378eLFC/z666+GiEGWaNo02d902TLVSYiSFRgon6f4+wMODqrTkLXjudm8aTSAlxcQECDLB4iILJlBitRy5crh8OHDAICzZ8+iaNGisfe5urrC0dERadKkga2tLTJlyoTXr18bIgZZiogIYNs2YOVKOZ49Kx1oiExYcDAwZowUqU5OqtMQ8dxsCTQaYOJE4McfgVu3VKchIjIcg0z3bdCgAY4cOYIOHTpAq9XCz88P27dvR0hICNq3b4/27dujU6dOsLe3R968edGyZUtDxCBLsXEj8O238W/76is1WYhSICxMGiSNGwdkyKA6DZHgudky2NrK7Izhw2VkNUcO1YmIiPTPIEWqjY0NfHx84t1WqFCh2MsdO3ZEx44dDfHSZInSptVd7twZ+OQToG5ddXmI3iMyUgrUESOAbNlUpyHS4bnZcjg4SKE6cqSsU82YUXUiIiL9YktUMn2tWgFlywK//ipTfv39gTRpVKciSkCrBby9gd69gbx5VachIkvm7Az4+cloapylxkREFsEgI6lEqfbTT8C765+uXZOz8NdfA6dPK4lFlFJarawVa9cOKFZMdRoisgYZMsiyAnd36SDOz2+JyFKwSCXTMHGibCaZmKAgoHJl4+YhSqVZs4A6dWTQn4jIWLJnl+UFHh7SAJ9bXRGRJWCRSmoMHAgsXQrMnQv07Al07Qo8f57wcY6OQN++Ro9HlBqLFwOffgrUqKE6CRFZo3z5gF69pKO4n590ASYiMmcsUsl4Ll0C+vSRab137wKhodJlBgDGjlWbjegDrVsnU+6aNFGdhIis2aefAm3ayMQkb28WqkRk3likknG8fCnDTC9exL+9cGElcYj0YccO+ayle3fVSYiIgPLlpYnSrFnAsGGq0xARfTh29yXjuHwZqFdPujq0bw+cOSM7kdevrzoZ0Qc5eBC4fp0FKhGZlpo1pXnb4sWqkxARfTgWqWQclSsDU6fKhm4LFwJlysgiGiIzdPIkcOQIMHiw6iRERAk1bQqkSyfLEYiIzBGLVDK8I0eA1q2B/PmB4cO56ziZtUuXgC1bAE9P1UmIiJLWoYNM/d2xQ3USIqLU45pU0q+rV2WOUXi47rbDh4E3b4CLF4ESJdRlI/pIt24By5YBU6awKQkRmb6ePWV96sGDQO3aqtMQEaUci1TSLz8/YMWKxO/Tao2bhUiPHj2SN3vch5CIzMmQIcCkSYCLC1Chguo0REQpwyKV9KtvXyB3biAsDMiTR3d7oUJAyZLqchF9hMBA2dZh2jTAwUF1GiKilNNoZP9ULy/A2Rn47DPViYiIkscilT7c7dvSsffxY7mePz9w/jxQqZLSWET6FBwsb/B8fQEnJ9VpiIhST6ORD9pGjgQGDZLTNRGRKWPjJPpwvr6yB0dQkHwFB6tORKRXYWGAuzswdiyQIYPqNEREH87WFvD3B2bOlOULRESmjEUqfZioKGDJErlcvz7w+jXwzz9qMxHpUVQU4OEhIw/Zs6tOQ0T08RwcgMmTZVQ1MFB1GiKipLFIpQ+zaZPuso8P4OoqXRmILIBWC3h7A716AXnzqk5DRKQ/zs4yEWrMGE6AIiLTxSKVPkyjRtLTvk8foHJl1WmI9EarlU6YbdoAn36qOg0Rkf5lyCDLGNzdZVkDEZGpYZFKHyZDBqBWLWDhQm4YSRZl9mz5p12unOokRESGkz07MGKELGuIilKdhogoPhap9GEePQIePFCdgkivliwBihUDatZUnYSIyPDy5QO++06m/nIrcyIyJSxSKeWePweePJGv9u2B/v2BtWtVpyLSi3XrgPTpgSZNVCchIjKezz6T5Q2TJrFQJSLTwSKVUsbbG8iSReYHZc8OHD4MbN3Kqb5kEXbsAEJDgXbtVCchIjK+8uVlBsns2aqTEBEJFqmUMnfu6C5nzSpf5cpx2InM3sGDwLVrQPfuqpMQEalTqxZQtCjw44+qkxARsUil9wkPB5YulSGmOXPktn//1U35PXVK5kcSmamTJ4E//wQGD1adhIhIvaZNZTe59etVJyEia8cilRJ69Ai4cgWYMQPo2RO4elX2Qb15k3tykMX4918gIADw8uKsdSKiGB06AEFBsgyCiEgVFqkU3759wCefSDHq6Sm3TZkC2NoC+fMrjUakL7duyZS2iRNZoBIRvatnT+C//4BDh1QnISJrxSKV4ps8Wdfer2hRoFQpoGtXpZGI9OnRI2DmTPmnbmurOg0RkWkaMkR6JJ48qToJEVkjFqkUX1iYHAcMkCm///wDNGyoNhORngQGyjYL/v6Ag4PqNEREpkujkf1Tf/5ZlkcQERkTi1TSeflStpU5flym+BJZkOBgecM1aRLg5KQ6DRGR6dNo5HfmkiWyTIKIyFhYpJI4fFhGT11cgIoV+S6eLEpYGODhAYwdC2TIoDoNEZH5sLWV5REzZgCPH6tOQ0TWgkUqiVu3pItvTLMkIgsRFSUF6ogRQPbsqtMQEZmfNGlkmYSPj0y6IiIyNBapBBw7BvTpA5w4ATx9qjoNkd5otTLFt1cvIG9e1WmIiMyXs7NM/fXykuUTRESGxCLVWty7B+zaBezeHf/2w4eBKlWA0FC5njmz8bMRGYBWK1vMtGnD7X2JiPQhY0ZZNuHuDoSHq05DRJbMTnUAMoKoKKBcORkldXDQdfAFgH794l8eN874+YgMYPZsoHZtoHx51UmIiCxH9uyyfMLdHZg2jVt5EZFhcCTVkkVEyDt0OzvdNN4vv4z/mBo1gMaNgb59gTlz2FWGLMKPP8o2vzVrqk5CRGR58uUDvvsO8PbWba1ORKRPHEm1FBs2AJs3y+VNm+R49y5w+rTuMc2aAdu2xf++hQuNk4/ISNavB1xdgaZNVSchIrJcn30GtG4t61THjJHtaoiI9IVFqjm7fh346y+53Llzwvtz5ABKlAAcHeVxdvzrJsu2Y4c09OjRQ3USIiLLV7488OaNLK8YOlR1GiKyJKxazFmDBsDNm/FvW7RId9nJSc4cJUsC9vZGjUZkbIcOAdeuAUOGqE5CRGQ9atcGgoKAJUtkCjARkT5wTao5y5dPjl9/DXzzDTBrFtC7d/zH1K8vI6pEFuzkSeCPP4DBg1UnISKyPs2aAS4ustyCiEgfOJJqbu7ckTk1gwYBBw6oTkOk3L//Alu2yLoorokiIlKjY0cZTd25kz0BiOjjcSTVXPz2GzBvHlC3LhAQII2SiKzcrVvSydfHhwUqEZFq330HXLkiyy+IiD4GR1LNwe3bQKNG8W+LiFCThchEPHoEzJwJTJ/OffqIiEzF0KEys8XFhftUE9GHY5FqDuLOmxkwQBoiDRqkLg+RYi9fAhMnAlOnAg4OqtMQEVEMjUa2pPH0BJydgU8/VZ2IiMwRp/uag4ED5bf82LEy5XfKFCBXLtWpiJQIDpY3QL6+8gaIiIhMi0Yjo6mLF8tkMCKi1GKRag5695buMBMmqE5CpFR4OODhAXh7AxkyqE5DRERJsbUFJk+WJRmPH6tOQ0TmhkWqKbt1C9izB7hwQYaPiKxYVJQUqMOHA9mzq05DRETJSZMG8PeXz9hfvlSdhojMCYtUU/XLL0CBAkDDhkCpUtLTnchKabUyevrdd7rtgYmIyPQ5O8vyDE9Pft5ORCnHItUUPX0KfP217nqDBkCdOsriEKmk1craptat2YCDiMgcZcwobTXc3WXZBhFRclikmiJvb93lFSuA338HsmZVl4dIodmzgVq1uJUBEZE5y5FDlmt4eMjyDSKi92GRaoo0Gjk2agR06aI2C5FCP/4IFC0K1KypOgkREX2s/PmBnj3ls3itVnUaIjJlLFJNkZsbMH8+8NNPqpMQKbN+vWwGH3ebYCIiMm+ffQa0aiXrVFmoElFS7FQHoDgWLgQcHIAePYDKlVWnIVJmxw4gKEg+cSciIstSoYL8jp8zBxgyRHUaIjJFHElV7dUrWaTh5gb06yfvyosVU52KSJlDh4D//mOBSkRkyWrXBgoXlmUdRETvYpGq2ubNwMyZwMqVuts2bFCXh0ihU6eAP/7gJ+tERNagWTNZ1sG3PUT0Lhapqm3erLu8fDmwdy/w+efK4hCpcvky8PPPgJeXrncYERFZtg4dZFIZt4MnorhYpKpWqJAcvbxkym+9enyHTlbn9m1gyRJg4kT+8ycisja9eskHlYcOqU5CRKaCRapq8+dLe7tJk1QnIVLi8WNgxgxg8mTA1lZ1GiIiUmHYMClST51SnYSITAGLVJX27gWCg1WnIFLm5UsZPZ08WRpbExGRddJoZP/UTZtkVJWIrBuLVFVOnAAaNAAyZQKOHFGdhsjogoOBMWNkEoGzs+o0RESkmkYj+6cuXizLQIjIerFIVeHePaBiRbkcHq67TGQlwsMBDw/51DxDBtVpiIjIVNjayuya6dNlOQgRWScWqcbm5wfkyaO7vmABYG+vLg+RkUVFAe7usj1w9uyq0xARkalJkwbw9wd8fGRZCBFZHxapxjZunO5yjx5A377qshAZmVYro6fffQfky6c6DRERmSpnZ1kO4uUFhISoTkNExmanOoDVePgQcHSU1nUnTsgU3woVVKciMhqtVtYatW4NfPaZ6jRERGTqMmaUDzbd3WX6LxvsEVkPjqQaQ0AA8MknwNKlQNWqwODBQJUqnOZLVmXOHKBGDaB8edVJiIjIXOTIIdvTuLvLchEisg4sUg3tyBEZOgKA5cuVRiFS5ccfgSJFgFq1VCchIiJzkz8/0LOndITXalWnISJjYJFqaBMn6i67u6vLQaTIhg2AiwvQtKnqJEREZK6KFwdatZJ1qixUiSwfi1RD++8/OfbuDXz7rdosREa2cyfw5g3QoYPqJEREZO6++EKWjcyZozoJERkai1RDiooCbtyQy82aqc1CZGSHDgFXr8oULSIiIn2oXRsoVEjafBCR5WKRakjh4dIgCQDq1FGbhciITp0CDh8GhgxRnYSIiCxN8+aAk5MsJyEiy8Qi1ZDSpgWOHpXFE87OqtMQGcXly8DmzdLgQqNRnYaIiCxRx47Aq1fArl2qkxCRIbBINRStFrh0SXUKIqO6fRtYvFgaW7BAJSIiQ+rVSz4YPXxYdRIi0rdki9SgoCDMmjULnp6e+P3333H79m1j5DJ/K1YAJUrITtREVuDxY2DGDGDyZMDWVnUaIiKyBkOHAgcPAqdPq05CRPqUbJHq6emJPHny4NatW8iSJQu8vLyMkcu8HT8OdOsml+3s1GYhMoKXLwEfHylQ06RRnYaIiKyFRgN4ewObNsmoKhFZhmSL1JcvX6JNmzaws7NDuXLloOXmVMnz8dFd/uEHdTmIjCAkRNafTprEpddERGR8Go2cg374QZadEJH5S9Ga1OvXrwMAHj16BBsbLmNNVswq/j59ZOdpIgsVHg64u8un2JzZTkREqtjaAv7+wLRpsvyEiMxbshXnmDFj4OnpiUuXLmHQoEHw8PAwRi7L0Lmz6gREBhMVBXh4AMOHA9mzq05DRETWLk0aKVQnTJBlKERkvpJdMHn//n1siLMR1a5du1C8ePH3fk90dDTGjx+PK1euwMHBAZMmTUK+fPli7//nn3/g7+8PrVaLrFmzYtq0aUhjSQvZhg0DHByASpVUJyEyCK1WRk979gTi/NcmIhNm9edmsgouLoCvL+DlJaOqTk6qExHRh0iySD1w4ABOnz6NnTt34syZMwDkBLdv3z40adLkvU+6d+9ehIeHY8OGDTh79iz8/f2xcOFCAIBWq4W3tzfmzp2LfPnyYdOmTbh//z4KFiyoxz+WAkeOAPPmAVWrSotTIgul1cobgNatgc8+U52GiFLKKs/NZJUyZpQPUt3dgenTZdyAiMxLkkXqp59+ipcvXyJNmjQoUKAAAECj0aBp06bJPumpU6dQo0YNAECZMmVw4cKF2Ptu3ryJDBkyYMWKFbh69Spq1aplGSfBjRuBDRtk9f6gQarTEBnMnDlA9epA+fKqkxBRaljluZmsVo4cMrHNwwOYOpVboxGZmySL1Jw5c6Jly5Zo0aJFvGZJT548SfZJg4KC4OLiEnvd1tYWkZGRsLOzQ2BgIM6cOQNvb2/ky5cPffr0QcmSJVGlSpWP/KMoFBYGzJ0rl9OmVZuFyIB+/BEoXBioXVt1EiJKLas7N5PVy58f6NFDRlV9fWUcgYjMQ7KNk+bPn4/KlSujfPnyKFGiBLrF7P/5Hi4uLggODo69Hh0dDbv/7xeaIUMG5MuXD4ULF4a9vT1q1KgR79Ncs/Toke6yu7u6HEQGtGGDrPVp1kx1EiL6EFZ3biYCULw40LIl4OenOgkRpUayRerhw4dx+PBhNG/eHLt27UL2FLTxLFeuHA4fPgwAOHv2LIoWLRp7X548eRAcHIzb/9/I6uTJkyhSpMiH5jcNX3yhuxznz0pkKXbtAl6/Bjp0UJ2EiD6U1Z2bif7viy+AatVkuQoRmYdku/tmyJABDg4OCA4ORr58+RAaGprskzZo0ABHjhxBhw4doNVq4efnh+3btyMkJATt27eHr68vhg8fDq1Wi7Jly6K2uc8djJk/0r692hxEBnD4MHDlCjB0qOokRPQxrO7cTBRH7drAmzfA0qUyBZiITFuyRWqOHDmwefNmpE2bFjNmzEBQUFCyT2pjYwMfH594txUqVCj2cpUqVbB58+YPiGtiHj0CnJ2B33+XTSNLlVKdiEivTp8GDh0CxoxRnYSIPpbVnJuJktC8ObB2rSxf4bgCkWlLdrqvj48PqlSpglGjRiFbtmyYPXu2EWKZgbNnpb3pzz8DGTIA5coB9vaqUxHpzeXLwObNUqCy2QQREVmCTp2Aly9lGQsRma4ki9TIyEj8/vvv+Pvvv5ErVy64uLigUaNGmDdvnjHzma7Jk4EHD4Bu3QC7ZAekiczK7dvAkiXAxIksUImIyLL07g38+y/wxx+qkxBRUpKsrkaMGAFbW1s8ffoU165dQ+7cueHl5YUuXboYM59pWrlS9kUFpNVprlxq8xDp0ePHsvn59OncV46IiCzTsGHyQayzs0yGIyLTkmSReufOHQQEBCA8PBytW7eGvb09Vq5cGW/9ilV68gRwc9Nd9/JSl4VIz16+BHx8ZOPzNGlUpyEiIjIMjUb2T/XwkEK1WDHViYgoriSn+8Zs+O3g4IDo6GgsW7aMBSoAXLqku7xuHVCpkrosRHoUEiKfuUyaJCdsIiIiS6bRyDlv0SLgzh3VaYgormQbJwFA5syZkSFDBgNHMRMFCwIzZwIDB8qmkVywRxYgPBxwd5dPlTNmVJ2GiIjIOOzsAH9/YNo0mSxHRKYhyem+165di90vLeZyjBkzZhglnEnKm5cbRpJFiYoCPD1lfU6OHKrTEBERGVeaNNIPc/RowNdXNm0gIrWSLFLjbjXToUMHY2QhIiPTaoGxY2Vj8/z5VachIiJSw8VFpv56ecmoqpOT6kRE1i3JIrVixYrGzGE+tm4F/vwTaNECqFFDdRqij+LnB7RsCXz2meokRJRa//77LzZs2ICwsLDY2yZPnqwwEZF5y5hR9gYfPRqYMQNwcFCdiMh6cYPPlAoLkx7lMY2T8uVjkUpmbc4coFo1oEIF1UmI6EO4u7vj22+/RQ7O0yfSm5w5ZfmLh4d0uudWbERqsEhNichI4OuvdQVqxozy7p7ITC1dChQqBNSurToJEX2oLFmyoG3btqpjEFmcAgWA7t2lmaCvL3tkEqmQbJH6+PFjTJs2DYGBgfjyyy9RrFgxfP7558bIZjoOHQJ275bLNWvKdSIztWGDbDHTrJnqJET0MXLlyoXFixfjs88+g+b/76KrV6+uOBWRZShRQpbD+PnJOlUiMq5kt6Dx9vZG69atER4ejgoVKsDX19cYuUzLb7/pPkb7+We1WYg+wq5dwOvXsnsSEZm3iIgI3Lx5E7t27cLOnTuxc+dO1ZGILMoXXwBVq8ryGCIyrmRHUsPCwlClShUsXLgQBQsWRJo0aYyRy7RMnQp89x3w/DmQJYvqNEQf5PBh4PJlWWtDROZv8uTJuHr1Kq5du4YCBQrgM3ZAI9K7OnWAoCBg2TKZAkxExpHsSKqDgwP++OMPREdH4+zZs3Cw1lZnRYoAlSurTkH0QU6fBg4e5Ba/RJZk1apV8Pb2xpkzZ+Dt7Y2lS5eqjkRkkZo3BxwdgY0bVSchsh7JFqkTJ05EQEAAAgMDsWzZMowfP94IsUzMr78C27cDERGqkxCl2pUrwKZN0gCCzR+ILMeOHTuwZs0aeHl5Yd26ddi1a5fqSEQWq1MnIDBQls0QkeElO933t99+w/jx45E+fXpj5DFNTZrI8dUrwN5ebRaiVLhzB/jhB9mYnAUqkWXRarWws5PTuL29Pex5fiIyqN69Zf9UV1fuQkhkaMkWqZGRkejWrRsKFCiAdu3aoVKlSsbIZTru3NFd5rt8MiNPngDTp0uByn3eiCxP+fLlMWjQIJQvXx6nTp1C2bJlVUcisnjDhgE+PtIlv1w51WmILFey03179OiBgIAAuLm5Ye3atWjYsKExcpmOvXt1l11d1eUgSoWXL4EJE4DJkwFr7HVGZA1Gjx6NVq1aITIyEq1bt8bo0aNVRyKyeBqNLJ/ZuFGW0xCRYSRbpL59+xa//PILZs2ahVevXmHQoEHGyGU6goPlmCOH2hxEKRQSAowZA0yaJJ/0EpFlOXDgAABgw4YNePz4MVxcXPDo0SNs2LBBcTIi62BjI+fYRYviT7gjIv1JdrrvV199hS+//BLjx49Hvnz5jJHJtMTMk2zRQm0OohQIDwfc3aVIzZhRdRoiMoSXL18CAJ4+fao2CJEVs7MD/P2B4cOBsWOBbNlUJyKyLEkWqZGRkbCzs8OWLVtimzGEh4cDgHVtQ3P2rOoERCkSFQV4eMh6GQ78E1muli1bAgAGDBiAN2/eQKPRYO/evahTp47iZETWJU0aKVRHjwb8/ABr7jFKpG9JFqmjR4/GjBkz0Lx5c2g0Gmi1WgCARqPBvn37jBZQubZtgdBQoHBh1UmIkqTVyhqZHj2A/PlVp6GUOn78OIYMGYLC///9EhwcjNy5c2P69OlwcHDAixcvMGXKFDx48ABRUVHImTMn3N3dkTVrVgDAyZMn8f333yMyMhIhISFo1aoVvvnmmwSvc/LkSVy8eBFubm5G/fO96/bt23B3d4dGo0GRIkUwbtw42NjoVp1ERETA3d0d9+/fh42NDSZOnIhChQol+n0ajQbu7u6YMGECHB0dFf6p1Bk1ahSqVauGM2fOIDo6Gnv27MH333+vOhaRVXFxASZOBDw9pVGhk5PqRGRo7567w8LC0Lx5c3Tu3Bnu7u64ePEiMmTIAK1Wi5cvX6Jbt25o3bp1gudZtGgRqlevjpIlSxr7jxDP/v378f3338POzg6tW7dGu3bt4t3/4MEDeHp6IioqClqtFj4+PihYsCB++uknbN68GZkyZQIATJgwAa6urli4cCHGjh378cG0yTh37ly868eOHUvuWz7a3bt3tUWLFtXevXvX4K9FZAkmTdJq//5bdQpKrWPHjmmHDBkS77Zhw4Zpf/31V210dLS2Y8eO2j179sTed+TIEW3Lli21kZGR2jt37mhbtGihffr0qVar1WpDQ0O1bdu21R46dCje80VHR2u7dOmiDQsLM/wfKBm9e/eOPYd4e3trf//993j379mzRzto0CCtVqvV/vnnn9oBAwa89/sOHz6snTdvXope2xLPK506ddJqtVrtt99+q9VqtdouXboY9PUs8WdIpC8PHmi1AwdqtSbwq5YM7N1zd1hYmLZOnTraV69eaUePHh3vPBwYGKitVq2aNjo6Ot5zPHjwQDt06FCjZU5KeHi4tn79+tqXL19qw8LCtK1atdI+efIk3mNGjRoV+17k8OHD2v79+2u1Wq12+PDh2vPnzyd4zunTp2uPHz+eotd/33klyZHUkydP4tq1a1i+fDm6desGAIiOjsaaNWuwY8eOj6+OzcHNm8DDh0Dx4kCGDKrTECVq7lygWjXgiy9UJ7EA79tm6ocfgF695PLixbJhXlL+P/MktcLDw/HkyROkT58eFy5cgKurK+rXrx97f9WqVZE3b16cOHECJ0+exNdff40sWbIAABwdHbF06VI4vfMx/pEjR1C4cGE4ODggKioKY8eOxaNHjxAYGIiaNWtiyJAhcHd3x8uXL/Hy5Uv88MMP+PHHH3HixAlotVp07doVjRs3xt9//4358+cDkIZ6U6ZMQYECBWJfZ/fu3VizZk281x45ciRKly4de/3ixYuoWLEiAKBmzZo4cuQIGjRoEHt/gQIFEBUVhejoaAQFBcXuAZrU91WtWhX+/v7o169fvBFZaxEREYFdu3ahcOHCePHiRexaVSIyvpw5gaFDZdnN1Knc+s2oFJ+7g4KCYGNjA9tE/tKfPXsGBwcHaN7JuG7dOnz55ZcAgEePHmH8+PEICwvDy5cv0b9/f9SvXx/NmjVD/vz54eDggAkTJsDLywuBgYEAgDFjxqBYsWJYvXo1fv/9d0RGRsLV1RXz5s2Ltyxz1qxZOH36dLzXXrp0aexjrl+/jrx58yL9/+eqly9fHidPnkTjxo1jHz969Gi4/n+Hk6ioKKT5/7YRFy9exOLFi/H06VPUrl0bvf//s23WrBnmzZsXe97+UEkWqenSpcOzZ88QHh4e25xBo9Fg5MiRH/WCZqVgQTkOGgTMmaM2C1Eili0DChQAatdWnYQ+1LFjx9C5c2c8f/4cNjY2aNeuHapUqYJdu3YhT548CR6fJ08ePHjwAE+ePMGnn34a7z7XRLbJ+vvvv1GsWDEAwMOHD1GmTBm0bdsWYWFhsUUqAFSuXBldu3bFoUOHcO/ePaxfvx5hYWFo164dqlWrhv/++w/Tpk1D9uzZsWjRIuzevRt9+/aNfZ1GjRqhUaNG7/2zarXa2BO1s7Mz3rx5E+9+Jycn3L9/H40bN0ZgYCAWLVr03u+ztbVFpkyZcPXq1QQ/C2vQs2dP7Ny5Ex4eHli1alXs3yURqVGgANC9uzRSmjTp/bUTmbeYc7dGo4G9vT28vb3h/P8tFaZNm4ZFixbhwYMHKFSoEOYkUkP8/fffaNWqFQDgxo0b6NatGypVqoTTp09j3rx5qF+/PkJCQtCvXz8UL14c06ZNQ+XKldGpUyfcunULHh4eWLNmDV6+fInly5fDxsYGPXr0wPnz51G+fPnY1xk6dOh7/xxBQUHx3js4OzsjKCgo3mNipvPeuHEDU6ZMiV1W0rRpU3Tq1AkuLi4YMGAADhw4gDp16qBw4cIJCuMPkWSRWrRoURQtWhTt2rVDNmtsWRb305T/jyQTmZKNG4G0aYHmzVUnsSAp/RS1Vy/dJ7MfqXLlypg1axYCAwPRvXt35M6dGwCQPXt23L9/P8Hjb9++japVq+LJkyd49OhRvPsuX74MrVaLzz77LPa2wMBAfP755wCADBky4Pz58zh27BhcXFxim+EBiB0VvXr1Ki5evIjOnTsDkCZ6Dx48QPbs2eHr6wsnJyc8fvwY5d7ZxT4lI6lxRzuDg4ORLl26eI9fvnw5qlevjuHDh+Phw4dwc3PD9u3b3/t92bJls7oRxJjGhrVr10bt/39CFfcDAyJSp0QJ2RDCzw/w8lKdxkooPHcnZuTIkahZsyYOHTqE6dOnI2/evAkeExgYGDsTKmvWrFi4cCE2b94MjUaDyMjI2MfFPTcfO3YMv/76KwDg9evXsLGxgb29PYYNGwYnJyc8evQo3vcCyY+kuri4IDhmu03IOTaxD7yPHTuGCRMmYOrUqShYsCC0Wi3c3NxiH1urVi1cunQJderUga2tLWxtbREdHf1Rs5ySLFIHDRqEuXPnxlb5cf35558f/IJmY+1aOTo4AGXKKI1C9K5du4CXL/X2u5ZMQMaMGTFt2jR06dIFW7duRbly5fDs2TPs378fdevWBQAcPnwYt2/fRsWKFZEnTx70798fTZo0QaZMmRAcHIyxY8eif//+8YrUTJkyxY48BgQEwNXVFT4+Prh9+zY2btwYrykeABQsWBCVKlXCxIkTER0djQULFiB37tzo2rUr9u7dCxcXF4wePTr2+2KkZCS1ePHiOH78OCpVqoTDhw+jcuXK8e5Ply5dbDf59OnTIzIyElFRUe/9vlevXiFz5swf8iM3WzGNDRs1apRgCplVNTYkMlEVKwJBQTIJb/Bg1WlIlVq1auHMmTPw9vbG3Llz492XKVMmvH79Gi4uLpgzZw7atm2LWrVq4eeff8aWLVtiHxdT5BUsWBBfffUVmjdvjufPn2PTpk24fPky9u7di02bNiE0NBStWrVKcG5ObiQ1pjnhy5cv4eTkhJMnT6JHjx7xHnPs2DH4+vrixx9/RK5cuQDICGyzZs2wa9cuODk54fjx47HNobRaLezs7D56GU6SRWrMD9MqCtLEfPutHBP5NIFIpT/+AC5flq1myLIULlwYnTt3xqRJkzB37lwsWrQIfn5++OGHHwAAOXLkwOLFi2Fra4vcuXNj5MiRGDBgAGxtbREcHIw2bdqgVq1a8Z6zUqVK2LNnD77++mtUqVIFw4YNw6lTp5A2bVrky5cPT548iff4unXr4u+//0anTp0QEhKC+vXrw8XFBS1atEC7du2QLl06ZMmSJcH3pcTo0aPh7e2NmTNnomDBgrHrcUaNGoUhQ4aga9eu8PT0RKdOnRAREYGhQ4fCyckpye+Ljo7G48ePYzssWosZM2YAkIL00aNHyJkzJ/755594o9ZEpFbdulKoLlsmU4DJOvXr1w+tWrXCwYMHY2e+AEDFihVx7tw5fPLJJ2jUqBF8fX3xww8/IGfOnLHrTuPq06cPvLy8sHHjRgQFBWHAgAHIly8f0qZNi1atWsHBwQFZs2ZN9bnZ3t4e7u7u6NGjB7RaLVq3bo3s2bPj5cuXGDNmDObPnw8/P7/Y7vuAjO76+Phg6NCh6NKlCxwcHFClSpXY9x9XrlxBGT0M8Gm075bc7zhx4gRCQ0Oh1WoxceJEDB48GM0NPL/w3r17qFevHvbt2xc79c2orl/XbTlz8iQQZ243kUqnTwM7dwJjxnCtC6VMdHQ03Nzc4k3vsRSHDh3CxYsX0a9fv2Qfq/y8YgBjx45Fjhw50K9fP0yaNAkajQZeBpxfaIk/QyJDW7NGJuW1bas6CZmS+/fvY8qUKQlGWC3B1KlTUbduXVSoUCHZx77vvJLsOOy0adOQP39+rFy5EuvWrcP69es/PLW5sLEBevYEypZlgUom48oVYNMmFqiUOjY2Nujfvz/WxixhsBBarRbbt29H165dVUdR5t9//40t0MeMGYNLly4pTkRE7/rmG+DFC+D/SwmJAAC5cuVCsWLFcP78edVR9Orp06cICgpKUYGanCSn+8ZIkyYNMmfODDs7O2TNmjVeow2LVaAAsGSJ6hREse7cARYtko3CWaBSalWuXDnB+k9zp9FoMH36dNUxlNJqtQgMDETGjBnx+vVrREVFqY5ERIno3RuYPh1wcQFq1FCdhkxF//79VUfQu6xZs8LHx0cvz5Vskeri4oJu3bqhU6dOWLNmDXLmzKmXFyailHnyRIrT6dMBu2T/xxKRtejfvz9at26N9OnT482bNxg3bpzqSESUhOHDAR8fKVTLllWdhsj0JfuWd86cObhz5w4KFy6M//77D22tYVL9pEnSmWbwYOCLL1SnISv26hUwYQIwZQrw/72TiYgAAHXq1EHNmjXx7NkzZMuWLUGnXyIyHRoN4O0NeHoCTk7A/7evJqIkJLsm9cWLF5g7dy6aNm2K2bNnf1BHR7Oydq38FlmzRuZYEikSEiJ7rE2aJJ+8EhHF9ffff6NFixbo2rUr5s6di02bNqmORETvYWMj5/RFi/gWkyg5yRapY8aMQYsWLbBu3Tq0bNnSoJ0DTULcLlvVq6vLQVYtPBzw8JAmSRkzqk5DRKZozpw5WL16NbJkyYI+ffpg3bp1qiMRUTLs7AB/f2DqVFnOQ0SJS7ZIDQsLQ7169ZAuXTrUr18fkZGRxsilzokTchwzBsieXW0WskpRUTIdaOhQIEcO1WmIyFTZ2NggQ4YM0Gg0SJMmDZydnVVHIqIUSJNGCtXx42VZDxEllGyRGhUVhStXrgCQzVktes3Lq1dAdLRcrldPbRaySlotMHasbPydP7/qNERkyvLmzYsZM2bg5cuXWLx4MT755BPVkYgohVxcZOqvp6cs7yGi+FI03dfT0xM1a9aEl5cXxowZY4xcakRGApkzy2X2CCcFJk8Gvv4aKF5cdRIiMnXjxo3DJ598gvLlyyNt2rSYOHGi6khElAqZMsnEPXd3WeZDRDrv7e4bFBSEAgUK4OeffzZWHrUyZwZ69AAyZABsbVWnISszZw5QpQobShNRyvTp0wfLli1THYOIPkLOnLK8x9NTOvnz7SeRSHIkdfXq1fjqq6/QokUL/PHHH8bMpNaUKdKxhsiIli0DChYE6tRRnYSIzIWrqyv27t2L69ev4+bNm7h586bqSET0AQoUALp1k+U+Wq3qNESmIcmR1B07dmD37t0ICgrCqFGjUMPSp79euwb4+gI1a8pvCiIj2bgRcHQEmjdXnYSIzEVQUBDu3buHFStWxN6m0WiwcuVKhamI6EOVKAG0aCHLfjw9VachUi/JItXBwQEODg7IlCkTIiIijJlJDT8/YPly4MYNFqlkNL/+CgQGAr17q05CROZi9erVWLZsGWxtbTF48GDUrFlTdSQi0oOKFYGgINkNcdAg1WmI1Eq2cRIAaK1h7sGqVXJ0cVGbg6zGH38Aly6xQCWi1ImZ6bRhwwaOnBJZmLp1pbs/l5uTtUtyJPXatWsYPnw4tFpt7OUYM2bMMEo4o5kzRzr7AoCbm9osZBXOnAEOHAC8vVUnISJzY9Yznf76C6haFahUCThyhF1iiBLx1VfAmjXApk1A27aq0xCpkWSROnv27NjLHTp0MEYWNe7dA4YM0V1v2FBZFLIOV67IOlQ/P8CStx0mIsMzu5lOp07J8fhx+XCYRSpRor75Bli0SJYFNW6sOg2R8SVZpFasWNGYOdSxtwe6dpX1qEePyvYzRAZy5w7www/A1KksUInow5j1TKdLl1QnIDIbffoA06bJSjRL719K9K737pNqFbJnl3kVjRvLJpVEBvLkCTB9upxw7Pg/j4g+kFnPdLp4UXUCIrMyYgQwYYIUqmXLqk5DZDzW/Vb5xQtZFPj996qTkIV79UpOMlOmAGnSqE5DRObMrGc6cSSVKFU0Gtk/1cMDcHYGihZVnYjIOJLt7vv48WOMGDECPXr0wMaNG3Hu3Dlj5DKOnTuBH3+U5g1EBhISAnh5AZMmsXk0EVmxJ0+AZ89UpyAyOzY2gK8vsHAhcPeu6jRExpFskert7Y3WrVsjPDwcFSpUgK+vrzFyGUdEBBAeDixdqjoJWajwcMDdXYrUjBlVpyEiUoijqEQfzM4OmDxZZmQ9eaI6DZHhJVukhoWFoUqVKtBoNChYsCDSWNJcRZ4wyYCiogBPT2DoUCBnTtVpiIgUy5gR+O67+Lc9fQrUri0tz4novRwdAX9/WT706pXqNESGlWyR6uDggD/++APR0dE4e/YsHBwcjJHLODJlkuO9e2pzkMXRaoFx44Bu3YACBVSnISIyAZ9/DixeDLx+LV979gDZsgGHDgHt26tOR2QWXFyAiRPlQ/CQENVpiAwn2SJ14sSJCAgIQGBgIJYtW4bx48cbIZaRxGyCXrmy2hxkcSZPBlq0AEqUUJ2EiMjEuLrKV/PmqpMQmaVMmYAxY2Q5UXi46jREhpFsd98cOXJg1qxZxshifDFFqr292hxkUebOld2MvvhCdRIiIoWiooDdu4H8+eUTu19+AYoVk693N4rOmBGIjpYOMUSUrJw5ZTmRp6esU7W1VZ2ISL+SLVKrV68ee/nly5fIkycPfv31V4OGMppt2+TI/9mkJz/9JNN769RRnYSISLH794FmzeTyhQvA11/LCGqlSkBkJFC1KnD0qNwfGCjTT7ZvVxaXyNwUKCDLisaOlR0E3v3sh8icJVuk/vnnn7GX79+/j/nz5xs0kFH9/jvw119A3ryqk5AF2LRJ9kDlDDYiIgC3b+suHz8ux+LFgcOHZY7i0KG6IhUALGmLOyIjKVEC+OorWWbk6ak6DZH+pGpeTa5cuXDjxg1DZTG+HDmAli2B8uVVJyEz9+uvMhDQqZPqJEREJiJukRrTTf+zz+I/Jl06IEsWudy4sXFyEVmYSpXka9481UmI9CfZkdRhw4ZB8//5A0+ePEHmzJkNHsooXr8G/vxT1soUL646DZmxP/6Q91/Dh6tOQkRkQu7c0V2+elWOxYrpbvvtNzkXFyoEPHvGJoZEH6FePSAoSJYddeumOg3Rx0u2SG3SpAnSpUsHAEiTJg1Klixp8FBGMXOmbDRVtSpw5IjqNGSmzpwB9u+X9SBERBRH3JHU//6TY5EiutsuXQL27dONoK5dy3fXRB+hRQtg9WpZftS2reo0RB8n2em+S5cuRcWKFVGxYkV8/vnnsLWUJkN//y1HbjJFH+jKFWDDBsDbm80KiIgSSKpIjYzU3X7smG4PjVOnjJeNyEJ9+y3w/Lk01iYyZ8mOpKZPnx4rVqxAgQIFYPP/1vBxO/6arZgOxRUrqs1BZunuXWDRImDaNO6YQESUqLhFalSUHAsXlq1mAOCTT+KvUeUvUyK96NMHmD4dcHEBLOEtO1mnZIvUjBkz4vLly7h8+XLsbWZfpD59qrsc0x6fKIWePAGmTpUTgF2y/4OIiKzUyZOyzjRNGpl2EhEBODnp7i9fXtqSzpkDDB6sLieRBRo+XFa1ubgAZcqoTkOUekm+xR4yZAhmz56NyZMnGzOPccycqbvMboKUCq9eyS/9KVPkfRcRESXB2Rk4fz7h7UOHArNmyUwmW1tpi84ilUivNBrpl+HhIZ8NFS2qOhFR6iQ5t+bFixfGzGFcadPKsU4dDoVRioWEAF5ewMSJ8skkEREl4eFD3bTed505I8eyZY2Xh8gK2dgAvr7AggWyTInInCRZod29excz4444xjFs2DCDBTKKsWPZjpVSJTxcPo308gIyZVKdhojIhB08KB8CA7I2IjAQOHoU8PQEGjaUdaiPH7NIJTICOzvA3x8YMULe+mbLpjoRUcokWaQ6OjqiQIECxsxCZJKiouS91ZAhQM6cqtMQEZm4GTN0l0eN0l0eOlSOCxbEf7ydHVClCpA+veGzEVkhR0dg8mRg9Gg58r8amYMki9QsWbKgZcuWxsxiPPv3yxyI6tU53ZfeS6sFxo2Trfv4mQ0RUQp06QLs2JHw9rh7pMaVIYOMtBKRwbi6ApMmyYfu06frVr4Rmaok16SWLFnSmDmM59EjoF49mYoUGqo6DZm4yZOl+WSJEqqTEBGZgaVLgXbtEr+vYEHjZiGieDJlkmVL7u7SbJvIlCVZpI4ePdqYOYznm290l9n9ht5j7lzZPYFb6RIRpVDPnknf5+iY+O1aLRAWJl9EZFCffCLNtD08dNsXE5ki69o5++5dmeoLyHQkjUZtHjJZP/0E5M8P1K2rOgkRkQV433qJ58+lgM2Vy3h5iKxYwYJA166ynEmrVZ2GKHHWVaR6eekuv9u4gej/Nm8GHBxkmi8REaVC3EZJcSW1HpWIlChZEmjeXDr/Epki6ypSV62S47ffyibjRO/YvVs+1I87K5yIiFIo7kK3zp11l+vXN34WInqvSpVkSdO8eaqTECVkXUXqn3/KanHukUqJ+PNP4OJFoHdv1UmIiMxQaChw5ozu+sqVwM8/AytWAH36GOY137wBvL2Ba9cSvz8ykvMZid6jXj0gb15Z5kRkSqxr/5Vq1eSL6B1nzwL79vHzCyKiD3b0KHDwYPzbWrUy3OvduqVb63rsGLBnT/z7X78GypYFihYFfv3VcDmIzFyLFjLZcNMmoG1b1WmIhEFGUqOjozF27Fi0b98enTt3xu3btxN9nLe3N6ZPn26ICEQpdvUqsH69fBjPXlpEZKkMfm6+eVN32RiL+uM2Y3r7NuH9s2cDN27IOg4ieq/OnYFnz4DfflOdhEgYpEjdu3cvwsPDsWHDBgwfPhz+iazKXr9+Pa5evWqIl0/c4cNA5sxA69bGe00yeXfvAgsXygbXNtY1+Z2IrIzBz81xi9QP6dT7/Dlw6lTyj7t0Cfjrr/i3FS8e/3pUlBSpRJRiffsC//wDHDmiOgmRgYrUU6dOoUaNGgCAMmXK4MKFC/HuP3PmDM6dO4f27dsb4uUTt3Mn8OIFEBxsvNckk/bkCTBlCjB5MmBnXRPficgKGfzcHLdITc3SGhcX4Pvv5XKFCjKck5Tbt4ESJYCqVRPed+MG8O+/ctnWFli3Tvf8RJQiI0bIzPmzZ1UnIWtnkCI1KCgILnFOCra2toiMjAQAPHnyBPPnz8dYYy/+e/pUjs+fG/d1ySS9egVMmCCt15PaX56IyJIY/NwcU6QePJi6FumOjkC/frrrgYFJP3blyvjXGzSQY3CwFK6DBwPR0XJbpUpytLVNeRYiK6fRSH+OdetkORSRKgYZP3JxcUFwnBHL6Oho2P1/qGr37t0IDAxEr1698PTpU7x9+xYFCxZEK0M2VwCkKgHYupUQGipb5k6cyA/Yich6GPzcfOuWHOOuFU2NwoWT7tIbY+NG3eX06YGmTWXYZ8MG6eT78KF0gGnfXja8/u47IG3aD8tDZKVsbABfX2DkSGDYMCBPHtWJyBoZZCS1XLlyOHz4MADg7NmzKFq0aOx9Xbp0QUBAAFatWoVevXqhWbNmhi9QASAgQI78RNWqRUTILkReXkCmTKrTEBEZj0HPzaGhwKNHgL39h61HfZ/bt4HNm6UAjTtFuW1b+UX+ySdSoAKAkxPQtavc5+QELF4MpEkD1KwZfw9XInovOztZDjVlim4yIpExGWQktUGDBjhy5Ag6dOgArVYLPz8/bN++HSEhIcZdhxojbiOGRo2M//pkEqKiAE9PYMgQIGdO1WmIiIzLoOdmOztpUPj4sf4/DM6fX47Vq0tB+sUXwNatQFiYFKBlygClS8tj/vlHjt9+K8dLl4Bp0+Ty1auynpWIUsTRUQpVd3fAz08mLxAZi0GKVBsbG/j4+MS7rVChQgkeZ5QRVABwdQWWLAGWL2d1YqW0WmD8ePmA/UNnohERmTODnpvt7YH/N2X6aMHB0jDgm2/keWP8+afc3qePXE+TRo5Lluge8/YtkC4d0LKljK527Ki7b+xYWdPq7KyfnERWwNVVlkd5egLTp3P2PBmPdWy6UbQo0KkTcOiQ6iSkiL8/0Lw5P0QnIjJ5bm6AhweQN6+8K47r4MH4wzkhIbIGNa7evWU96tu3upFVQJb9vPtYIkpWpkyyTGr0aM6aJ+Ox/CI1NFT6aN+5w/WoVmrePGnyWLGi6iRERBZq9Wp5B3v69Ic/x759wPXr8QvLGTPiP6Zhw/jX584FXr6UyzGNBgYOlGPM1jYZMugeH7N2lYhS5ZNPpHm2h4eugTaRIVl+kbptG1C2LNCmjeokpMBPPwH58gF166pOQkRkwbZuBaZOBa5c+fDnyJtX3gm/T/Pm8a937Qq0agUcOwacOCFrUOO2Iu3XT7a06dZNrjs5fXg+IitXqJD8lxs3TpZRERmSQdakmpRly+R4/braHGR0mzbJcqavvlKdhIjIwsXskfqxi/6PHk36vlWrZJuauHLkAH7++eNek4hSrGRJ3bJxDw/VaciSWfZI6p07wO+/y+UePdRmIaPavRt4/lzX4JGIiAxIH0XqwIGyr2m1avFv/+474MAB/fxCf/YMePPm45+HyIrFLKGaN091ErJklj2S2qKF7vLw4epykFEdOSJb6Y0YoToJEZEVCAuTKbW2tkC2bB/+PLt3AzduyDFzZunaW7GibDnj4vJxGWNmU40eDfz3X/yOwESUavXqyec9y5fLFGAifbPskdT79+U4YgT3HbESZ88Ce/fyMwkiIqMJDJRjpkyARqOf58yUSQrKOnU+vkAFgOLFdZe5mI5IL77+WrZI3rxZdRKyRJZdpI4eDXh7A+/sC0eW6epVYP16+SvX1/skIiJKxosXcozprvuhrl2L/3z65OMD1K8vl5cuBX75BfD1lXmLr17p//WIrMS33wJPnwK//aY6CVkay57uy+E0q3H3LrBwITBtGmBj2R+9EBGZluho6aKvrxlL8+dL8ahPWbNKF729e+X611/r7jt4ULoG8+RB9EH69pX3Xy4uCZeUE30o/kYms/f0qex8MHmyTDshIiIjKllS9kfVV5fduP0k9CkqKvHbv/5a1tPOnm2Y1yWyAiNGSK/Ss2dVJyFLYblFakgIsHixzP8ki/XqFTB+vBSojo6q0xAR0QfbsAHw9ARatzbM87u5Af37A+fPy4aP7xo27ONf488/gV27Pv55iMyMRiP7p65dK8uviD6W5Rapt28DvXtzyq8FCw0FvLyAiRP101eDiIg+QESEfpoRtWsn60QN1VQgY0aZSlyypGxpc/Uq0KmT7v5SpaTz74cICwOKFQNq1ACaNpUPyomsjI0N4OcHLFgA3LunOg2ZO8stUi9elOODB2pzkEFERADu7lKkfmyvDiIi+gjjxwP29rIozVzkyQMUKQK0aiXvrN3cgH/+0TVXSq1Hj+IPHz1/DmzaFH9PVnYVJitgZwf4+8vX06eq05A5s9wi9fJlOZYqpTYH6V1UFODhAQweDOTMqToNEZGVe/FCfjE7OalOknqtWwOvX+s+2L5z5/2Pv3VLNohcvVquv3oFbNkCpE8f/4SUN6+MDKdLB8yaJS1QP/0UCA42yB+DyJQ4OsoyrPHj5b8X0Yew3CLV1VWOefOqzUF6pdXKmoeuXYGCBVWnISIivW1Bo4qzs5xYkvPjj9LBeP9+WUq0dSuQIYOMxu7dK8uM6tRJ+H3DhgFr1shI66pVSTdwIrIgrq6y85OnpyzPIkotyy1SN22SY2LNEchs+fvLTgElS6pOQkREAMy/SAWAmjXl6OoKHDkia1WfPAE2bgQ2bwb27QO++073+CdPgJYtddfv3JEpz3v3Am/fAh076u5bu1Z3uW9fYPt2eY3bt3W3x50WTGQhMmeWItXdXZZpEaWG5RaptrZyjI5Wm4P0Zv58oGJF/W+fR0REHyGmSM2cWW0OfXjzBqheHVi3DqhcGWjfHmjbVpohTZyY+PfEdA0GZH1rmjRA9+7SqGnRovgFKyDFbfXqutvd3WW6cNxidtMmoG5d4MYN/f8ZiYzok0+AQYNkmRbfklNqWG6ROmoUMHOm/K8gs7d8ufS5qFdPdRIiIorHEkZSNRogX774t928qbscFgYMHChrTu7eBSpUAAoXBq5flxNUmjTxv7d+fWme1Lu3XO/eXYr4mOsA8NdfQMOGwJQp8rzffCMjths3Ah06SAfi334zyB/3g9y8yU0w6YMUKiTLtMaNY/8wSjnLLVKbNgWGDpWPcMisbd4s3eIMtb87ERF9hOfP5WjORaqrqzRAsonztmjsWN3lW7ekCRIA5M4t03UPHHh/c4S4W+n8+CNw/z4wZEj8x7zb3LF+fRm9jRlyMqV39CNHynQmdsKhD1CyJNCsmSzbIkoJyyxSv/wSmD6d289YgN9+A549k8aIRERkgr7/XjrYxhRx5qpIEWlyBEgn3lKlgJ07gStXZJQ1btHp4CDFakppNDLa+umn8d+b+PomPuOrSBEpgGN+prt3A+vXp/z1Xr4EDh2SYjo8XHd7ZKTsfpDS4jduN+Ldu2VhYfr08mcxpVFeMguVKgFffAHMm6c6CZkDyytStVrg99/lE7+Ydalklo4cAc6fB/r0UZ2EiIiS1LmzjBDamPlbChcX2es1Ohr491+gTRugSROgaFH9vk7OnNIZ+OxZ2atjyBDZq+PmTWDOHPm6ckUaMXTuLAVu48ayhjWpUczjx4GyZYFdu4AffpD1sLVry9rX+fPlMY8fS+H92Wcyanz7NuDtLZ1tYoSFSeFcooS8rouLTD8GgJAQ3ePCw4FTp6SIDQrS78+HLFr9+vL5zvLlqpOQqTPzM0oi/vtPdzlbNnU56KOcPQvs2SNd/omIiIwm7oipobRoAXz+uVzOlk0W6+XPLx1mBg2SDL//nvD7tFrgn3+ALl2AbdukWFy+XJo8nT0rS51On47/PXfvAmfOyBBWzB7y7u5SsE6aJBtaurkBgYHA4sXAmDHApUu6748ZMY0ZZQbkcf/+K0VsuXLJ/3kDAqTInjAhhT8gsmQtW8o40s8/q05CpszyitQff9RdNsaJhvTu6lVprDh2LP8KiYhM0qxZMsq3ZQswdy7w66+qE1mer7/WXf7pJ9lftX9/KW5XrZItcdKkAbp10z2uYEEpGitX1t22bRtQrZoUqzFy55ZpuzFWrgTatQOWLImf4aefpEMxIMuotFr5KlMGWL1abo87OJCYFy9kS5/162XE2JTW2ZIynTvLTk6cNU5JsbwiNUb16qoT0Ae4exdYuFBmG5n7zDEiIou1f78UKb/9BgweLFNUSb9q1ZLuwZGR0pjBzQ1Ys0Z3f9x94AcMAK5dk2nCvXtL5+Bbt+TvqXNnIDRU2quGhMhj9u2TKcIZMsj3f/65FKPbt8s04KdPpZjs2lX2f31X2bKyrCrGkCEygluvHpA3r3zCrNHopgJnzarXHw1Zhr59gXPngKNHVSchU2SnOoDBNG+uOgGl0tOn0ol/+nTp5ktERCYqNFSOYWFyjCl2SL9iugdv2CAjkba2wIgRMtXIyQk4dkwKwfr1E35vvnzyVaeO7Llao4YUjjFrbD/5BFiwQIrW0aOBtGlTnsvOTk7Y06bJdRcX2U4nptNzDG9vGXW/c4efPFOiRo6UAXZnZ90MeCLAkotUMiuvX8svKX9/6SNBREQmbN8+Oe7dK0cWqYa1eLF0+t2+HahaVXd73Gm971OzZuK3d+z4cbmKF5duxZUrS3OnuLPY2rcHJk6Uy1y7Q0nQaGRJtoeHfO5SpIjqRGQqLK9IHTxYNsHm/qhmIzRUmgtOnChb1RERkZm4d0+Ocdc3kv799ptM+zWlT3E1Gunwa2urmxIcGSkjpkkVpc7O8phbt4DChVm8EgD5J+PrK5MERoxI3e5OZLksb+5FrlzSNCBHDtVJKAUiIqTJoKenee8DT0Rk1TiSalh2dqZVoMZwdIy/ZtXWNunCs2RJ2RandGmZchxTzH73nVGikmmzs5PZdP7+svyLyPKK1J9+it86nUxWdLQUp4MHc+CbiMissUil5Jw4IdOLY7bBibFsmZo8ZHIcHWVHpHHjkt4SmKyH5RWp3btLN7v791UnoffQauWXkJubri8EERGZIVtbFqmUPEdHGUXt2hX49FO5rWxZ2U4HkC12EpPU7WSRXF1l+ZeHh64/G1knyypS37yR48CB3IfLxPn7A82ayewfIiIyMzENd9askXUb7dqpzUPmoVQpmfH277/yPu30aeDwYZn2a2cn3YJr1ADmzQNmztTd3qCBbHPD/XitQubMMtPO3V1+vZB1sqwiNe5HLlx1bbLmzQMqVgQqVVKdhIiIPppGI6OpRB8ipkM0AIwaBVy9KsNpw4fHf8ycOdKoiaxCrlzAoEFSrEZHq05DKlhWkRqDm0abrOXLgTx5ZL9vIiIyQ1qt/CKvVYtNCunj+fvHvz5qFNCpk7xhiOvzz4HgYKBtW1YtVqJQIaBLF1kexgmS1seyitSXL1UnoPf4+WeZtfP116qTEBHRBwsJAaZOBU6elPl4ZcsCz5+rTkXmqk0b4L//5N/Uo0cygurgIE0rtFrg4UOZKXf2rHxt3gzcuaM6NRlJqVJA06bAlCmqk5CxWVaRunmzHAMD1eagBH77DXjyBPj2W9VJiIjoo7x9K8e0aYF//pHCwclJaSQyc4ULAyNHAtmzJ7wvRw7d9jv//itHb29dg8ywMGDBAmDYMCA8PHWvy+E5s1C5MlChAjB/vuokZEyWVaTGrIn56iu1OSieo0flfUzfvqqTEBHRR4spUl+8kALB3t409/AkyxMzY271auk9snAhkC8f0L8/MGsWUL48cPeu7vEhIbrCNoZWKyO3AwfK2tfmzbnW1QzUry/rVN+dBU6Wy7KK1NGj5ZfPzz+rTkL/d+4c8PvvwIgRqpMQEZFexBSpMesCM2SQ5klEhta0qe5ys2bAd99JN+AYFy4AGzZIcTpzJpAlC1C8ePzmTJs3A0WLyrBccDCwY4cM05HJa9lSxqP4Nt86WFaRSiblv/+AtWuBsWP5/oWIyGLEFKkxuEcqGcvq1cCtW8Ds2cDSpdLoYtMmWc8aI00auW3CBN2uDw0ayEgroBtZLVw4/nNPny7NwLje1aR17gw8fizLyMiyWU6RGh0tv5Ri1qWSUvfuyRIRX1/AxnL+lREREYtUUkWjkem9gwcD2bLpbh85EjhwQEZNS5cGfvhBNmL/5BPdY2Km9Hp5Se+S//6TgjQwUNaynjole7YeOmTcPxOlWr9+shT+6FHVSciQ7FQH0JvISNlM3N5eOsWRMk+fShe2adPkQ04iIrIgLFLJFNWuDVy8KJdjqpebN2VKl4ODvEcEZL5ozL/ZPHl03x/zhoXNlMzCqFGyNY2zs+xORJaHJQTp1evXwPjxsu0Z+2gQEVmgokVlREurlZGqVq1UJyJKXIECMnKaGm5uMhpbpowMgJw7B5w+DXToII2WyCRoNPJ+091dmosXKaI6Eemb5UzEjIhQncDqhYYCnp6Ajw9/jxMRWaysWeUNOyBTZvr0UZuHSB/iNs/w8QGePQM+/VSaKvXqBSxbpi4bJcrGRpaVff+9LDMjy2I5RerBg3JksapERIR8muXpCWTOrDoNERERUSq0b6+7XKQI4OcHXL+uu23RIl0lxCnBJsPeXmbv+fvL5wpkOSynSF20SI6JbQRNBhUdLcXpoEHxexQQEZEFunBBPpmsX1+mQz54oDoR0cdr2lSKz3v3pEB1cwP27ZOmTABw+bLsp7dsGZA3r2xfQybB0RGYPFnWqL5+rToN6YvlFKl79sixdWu1OayMVitrAtzcgEKFVKchIiKDO35cOunv3Su//DduVJ2ISH9y5ZLmSp9/DtStq9tDtXZt2cYmIEAK2fbtZWub9u2ls/CFC3JfzP7BZFSurjJL29NTt/MQmTfLKVIPHgQ8PLg2xsimTJEPH0uWVJ2EiIiMIiws/vV06dTkIDKGdu2ko/WBA0DOnMD583L7zp3yQc3GjcD9+3Js3Ro4cQJ4/pzFqgKZM0sp4O7O1X+WwHKK1MqVZXpGqVKqk1iN+fOBL74AKlVSnYSIiIzm3S1oWKSSpUuTRnf53feZ330HtG0rU4MBeT+aJQswbFjiz6XV6grYy5eB5s2B4cNl+5ygIP1ntzK5csnyM09Pfk5g7iynSH33k10yqBUrgNy5gXr1VCchIiKjerdIZTt3siabN8tIqVYrX4sXSxH77v+DEyfkGBUlW9hMmCDtaG1spNCNjgYePgR27ABmzpQpablyAYGBxv8zWZhChYAuXWSNKntcmS/LKFIfP5aNmb/5RnUSq/Dzz7Jc4+uvVSchIiKj40gqWTNHRyBTpoS3e3gAffsCnTrJ9Xr1gJ49ATs7oHx5aeARUzF9+ikQHg7kyRP/OV6/ZiMyPSlVCmjSRJalkXmyUx1AL+7elZPm2rXAmjWq01i0334DnjyR38NERGSFWKQSJVSrlnwBwIIFgLOzDKDEaN8eyJ8f+OoroEoV2Ze1cGHp8nPzJlC8ONdP6VmVKjKD+vvvgf79Vaeh1LKMIjUyUo4VK6rNYeGOHgXOnQNGjVKdhIiIlLF7560Di1Si+NKnl+MffwC3b0thapPE5EVHR+Czz4AvvwS6dUu4VUJ4OPDmDTeh/0ANGkihumKFNCMn82EZ031jitR3T5ykN+fOyShqzHZhRERkpfz8gI4d5fKYMUCOHGrzEJmqsmVlbVRSBWpcmzdLw6WVK4GGDaVY1WhkvWvOnMC//xo8rqVq2VJ+lAEBqpNQalhGkRoVJUcWqQbx338yk3rcOPlPTkREBEBGgOztVacgMn8uLsC6dUDv3sCePcCNG7r7fvwRyJtXXTYL0KUL8OgR8PvvqpNQSllGkbplixxZpOrdvXsyl9/XN2UfBBIRkRVwcZHph3G35iCij1O7tu5ykSLSbOn4cen827SpjBRoNNJRmFKtXz/gzBlZvkamzzLKjrt35RgSojaHhXn2TLqi+fuz/iciov/r0UOmJmbIABw+rDoNkeX49lvd1jZXr8oUtooVgYMHgUOHdI/r3Vu+uBFoqo0aJcvX/vlHdRJKjmUUqbNmAfv3A7t2qU5iMV6/lt+Nfn6ypp+IiAiAfDAcGAhcvy7v9ojIsHr3BubNk31VYuzeHX9KMKWIRiPvb1evluVsZLoso0jNmxeoUwfImFF1EosQGgp4egI+PtyjnYiI3hF3CxqeJIgMz9kZGDBAhv9Wr5Y1WH/9JVvYUKrZ2MiPcP58WdZGpsn8i9QbN2STZNKLiAjZj9rTk93OiYgoEeHhussnT6rLQWSNvvlG3qRt3y5dgLt1k1HWCxdUJzMr9vaypG3yZFneRqbH/IvUQ4eAs2elZTd9lOho+b03cCDwySeq0xARkUkKC1OdgIhOnZIuwMuXA4MGyVTgc+fkviNH+P80BRwdpe/KuHGyzI1Mi/kXqa9eSXXFVl0fRauV/6RduiTcR5qIiChW3JFUW1t1OYisWc2aCW9bvx54/FgaiixbZvxMZsjVFZgwQQZpQkNVp6G4zL9IPXBAjlqt2hxmbsoU6W4ed00+ERFRAnGL1HTp1OUgsmbffiuDNFqtLLD09JSpwIsXA8eOAQ8eqE5oNrJkkaVu7u6y7I1Mg3kXqS9fAtu2yWVnZ6VRzNn33wMVKgCVK6tOQkREJu+773SXS5ZUl4PI2mk0cvT0lEI15v/jixfApElA0aLSaImSlSuXLHfz8uLOPqbCvIvU8+d1l93c1OUwYytWyPrT+vVVJyEiIrMwahTQsaNc7tNHbRYiii9tWt3l//6TrsBRUTIVeMwY4No1ddlMXOHCMkA9fjwnaJoC8y5SY8bka9UCPv9cbRYzFBAgbbhbtlSdhIiIiIg+WteuMhwY0wGze3d5s9exo4y2FikCtGkj0+dGjOD81neULg00bgxMnao6CZl3kfrZZ9LVbNQo1UnMzu+/y9r6zp1VJyEiIrNy+DCQNass4ipeXHUaIoorSxaZ6nv/vgwHzpwp04LLltU95uefgePHZTrwy5fKopqqKlWAcuVkORypY6c6wEfJmZPTfD/A0aOyaw9reyIiSjU3NyBHDuDuXdVJiCil9u8H+vcH1q6V699+C5QvLwVscDB7u7yjQQMgKEiWxbHUUMO8R1JXrgR++011CrPyzz8yijpypOokRERkttKkUZ2AiFIjQwZgzRoZXdVqgVWrgN27ZVaEn5/qdCapZUup4QMCVCexTuZdpLq5Ab/+KlMaKFn//Se/n8aO1TWEIyIiSjUWqUSWw89P3hhyil0CXboADx8Ce/aoTmJ9zLtIBYA5c9grOgXu3ZO59ZMmyfp5IiKiD8Yilcj8pU8f//rs2UpimLr+/YHTp4G//lKdxLpYRrmSJ4/qBCbt2TPA31++7O1VpyEiIrPHIpXI/Lm5yQhGz55yPTJSvthMKYFRo2Ty5j//qE5iPcy3SH3xQnUCs/D6NTBuHDB5MuDoqDoNERFZBBapRObP1RXo1w9YvFjeVz9/LpuEZswo03/37AFu3ZL5rlZOo5EfzerV3GrWWMy3SF23TnUCkxcaCnh6Aj4+8nuIiIhIL1ikElkOjUYK04wZZfuHGA0bAgUKyJ6rYWHK4pkKGxvZanbePLbDMQbzLVJjVjBXqqQ2h4mKiJAt7Dw8gMyZVachIiKLcfgwsHSp6hREZAg//QR06BD/tuzZAQcHNXlMjL09MGWKzFB89kx1GstmvkVqw4ZA6dLAoEGqk5ic6GjAywsYOBDIlUt1GiIisih58gC5c6tOQUSGkDWr7KX655/AxYtAVBTw6BGwcaN0ENq6VaozrVZ1UmUcHaVIHTdOltWRYZhvkdqvH3DuHNCpk+okJkWrlTnznTsDhQqpTkNEREREZkWjAapVA4oXlzmuM2fK6OqCBbJ5aNasQI0aqlMq5eoKTJggy+pCQ1WnsUwGKVKjo6MxduxYtG/fHp07d8bt27fj3b9jxw60bdsWHTp0wNixYxGd2i1koqM51SgJU6cCTZoApUqpTkJERKZEb+fmvn2l0QoRWYfq1YGSJePfduSImiwmJEsWWVbn7i7L7Ei/DFKk7t27F+Hh4diwYQOGDx8Of3//2Pvevn2L2bNnY+XKlVi/fj2CgoJw4MCB1L3AiROy8fCFC3pObt6+/x4oXx6oXFl1EiIiMjV6Ozfv3h2/uQoRWbaKFYHz52W63pMngIuLvOEk5Moly+s8PWUMjfTHIEXqqVOnUOP/0wDKlCmDC3GKSQcHB6xfvx5p06YFAERGRiJNarsEhoYCd+8CgwfrLbO5W7FCmq/Vr686CRERmSK9npvZ3ZfIOmXNCly/LnuxDBokU4PjfpUqZXXrVQsXlmV248db3R/doAxSpAYFBcHFxSX2uq2tLSIjI+UFbWyQJUsWAMCqVasQEhKCatWqpe4FTpyQcfWoKL1lNmcBAfJ7oWVL1UmIiMhU6fXczI23iaxXtmzAp58m3p3zwgXZc9XKlC4NNG4MTJumOonlsDPEk7q4uCA4ODj2enR0NOzs7OJdnzZtGm7evIl58+ZBo9Gk7gX++EOOd+/qI65Z+/13abrWr5/qJEREZMr0em7mSCoRDRkiDZb27wcePpQpwQDw8iXw4IE0XrK1VZnQqKpUAYKCpL8U35d/PIOMpJYrVw6HDx8GAJw9exZFixaNd//YsWMRFhaGBQsWxE4tSpVz5+TYpMnHRjVrf/0FnDnD/whERJQ8vZ6bWaQSUZo00lRp7Fhg4ULZtubPP6W5aenSUsS+fas6pVE1aADkyAGsXKk6ifkzyEhqgwYNcOTIEXTo0AFarRZ+fn7Yvn07QkJCULJkSWzevBkVKlSAm5sbAKBLly5o0KBByl/gzh05NmpkgPTm4Z9/pHfF+PGqkxARkTnQ67mZRSoRJeXSJTnOny9fK1fKok0r0aqV9IoJCJDL9GEMUqTa2NjAx8cn3m2F4mzaefny5Y97gbRppXmSle7RdO0asGaNbCSc2pnSRERknfR2bq5WDShQQJ/RiMiSjBgB/PKL7nqfPlZVpAKAm5vsurFnj4yuUuoZZLqvwYWESNMkV1fVSYzu/n35UGrSJNlfmYiIyKhWr2anPiJKWvXq8j495kOxkBC1eRTp3x84fVqW51HqmW+ZY2NjdcOIz54B/v7yZW+vOg0RERERUSJsbIAxY2SJXswyPSs0ahTw66+yTI9Sx3yLVCvz+jUwbhzg58fO/0REpNCrV7INHBHR+2g0QJ480u332DFpfWtlNBrpH7NqlSzXo5QzvyL1zh0gd26rWo8aGgp4esqsCSuc4UxERKakTBlg1y7VKYjIXHTvLvuz/Puv6iRK2NjIINO8ebJsj1LG/IrUt2/lb/jRI9VJjCIiAnB3Bzw8gMyZVachIiICu/sSUcpptXI8cACIjFSbRRF7e1muN3ky8Py56jTmwfyK1AcP5OjgoDaHEURHA15ewMCBQK5cqtMQERH9H4tUIkqpmB4yo0cDP/2kNotCadPKiOrYscCbN6rTmD7zK1LXrJFjcLDaHAam1coc9m+/BQoXVp2GiIgoDhapRJRScbef6dVLhhStVLp0wIQJMkPy7VvVaUyb+RWpMZ/GNG6sNoeBTZsmf8TSpVUnISIiegeLVCJKqUGDgOXL5fLIkUCXLkrjqJYliyzlGz2aPejex/yK1Ji/zS++UJvDgBYsAMqWlTXmREREJodFKhGlhpubzIKcMgX45BPVaZTLnVuW83l5yfI+Ssh8i1QL3Sh05UogZ06gQQPVSYiIiJLAIpWIUsvJSTcjklC4sCzrGz9e11uKdMyvSJ09G7h+HWjRQnUSvQsIkGPLlmpzEBERJWnePPk0lYgoNY4eBcqVk0J1/nxg/35g926rrtBKl5blfdOmqU5iesyvSM2SBShYUFYeW5A9e2RXHSufpk9ERKauWTPAxUV1CiIyN4GBwJkzcnngQKBePanQxoxRm0uxKlVk++kFC1QnMS3mV6RaoL/+Ak6fBvr1U52EiIiIiMgAvvxSZmIA8beStLWVhZnOzjLKunKlmnwKNWwIZM8OrFqlOonpMK8iNSJC/vHa2QGhoarT6MU//wC//gqMGqU6CRERUQosWgRERalOQUTmxs4OGDBAfn8EB0thevu2bBwKAGFhcnRzk/f7XbsC164pi2tsrVvLj2TLFtVJTIN5FalHjsgxKsoiGidduwasXi0LprmOnIiIzMKUKTxpEdGHs7GRglWjAfLmlcs2NtJ3Jq4VKwA/P7ms1QJ37gDHjskb6Dt3dEWtBXFzAx48APbuVZ1EPfMqUmP+MebIIf+gzdj9+7Jm3NdX/l8SERGZBXt7nriISP8GDADu3QP8/WWt6po1wMKFQOfOQPr0QL58soCzSBG5XKaM6sQG0b8/cPKk1OPWzLzOMuPHy7F0aaUxPtbz58DkyfJ/0AIGhImIyJrEXUtGRKRPuXIBo0cDu3YBnTrJdlfh4cCbNwkf++yZLP/75Rd5jAWJ+RGcP686iTrmU6RqtbqPFHLnVpvlI7x5I1PvJ08GHB1VpyEiIkolfrpKRMY0dSrw+LHUAtHRuq+nT6VQbdsW2LBBdUq90mhkbG7lSqtalhuP+cyZvXtXd9nXV12Oj/D2LeDpCUyYALi6qk5DRET0ATiSSkTGlC+f7vK76+FnzZLGqlOnyn2OjsDnn8uUYDNnYyNLckeMkAaruXKpTmRc5lOkZs0K7NsnC6Vz5FCdJtUiIgB3d/nKkkV1GiIiog/EkVQiMhW2tnK8cEHWrgJS1J4/D7x+DVy9CpQtC2TIoCzix7C3l+WBI0fKIFfmzKoTGY/5TPdNmxaoW1faUZuZ6GjAy0vWg1vbpyBERGRhzPTNHhFZIDc3oGbN+M3cWraUPVcHDZLaoWdPdfn0IG1aGVH19k58aa6lMp8i9dgx4NEj1SlSTauVOeXffgsULqw6DRER0UfasUN1AiIiUbIkcOiQbE8Zs2Y1ZusMFxd5zL17ajPqQbp0gI8P4OEhywetgfkUqb17AzlzApMmqU6SKtOmSRdtM29ITERERERk2jQawMlJLvftK0etVl0ePcqSRbd0MDJSdRrDM58i9dYtOX7xhdIYqbFggUyDr1JFdRIiIiIiIiv0999SvMYtVs10y5rcuWUfVU9PGTS2ZOZRpL55I4ufASBPHrVZUmjlSunv1KCB6iRERER61L+/6gRERMkrWFA3qgpIkXrlijRjTZNGCtdCheSo0QB79qjLmgpFigDffCONlCxkkDhR5lGk7t+vu1y0qLocKbRlixxbtVKbg4iISO+ePVOdgIgoedmyyf6qx47Jm3OtVubMxv0dduOG7nLM7S1aSCOZly+NGjc1Pv8c+PJLWVZoqcxjC5qYJg0aDWBn2pH37gUePOAHzUREZKG4TyoRmQsXF6BSJd31TJmkkdKLF8D9+zJD89w5YP58oGNHecy2bXI8fx6oUcP4mVOoalUgKEiWF/brpzqN/pnHSOqlS3Ls3l1tjmQcOwacOsUClYiILJiJf1hMRJQkjUb2gyxVCmjUCChRAujUCTh6VPeYMmXk+OqVkoip0bAhkD07sGqV6iT6Zx5F6qRJwIgRwMCBqpMk6fx54NdfgVGjVCchIiIyIBapRGTJYrauad4c2L1bbZYUaN1amijFLDe0FOZRpNapI5OuP/9cdZJEXbsmn2CMGycf0BAREVkse3vVCYiIDKdePd3lf/5RlyMV3NxkueHevaqT6I95FKkm7P59mcYes28wERGRRWORSkSWbPx4GXkCgMaNgcBAICpKaaSU6N8fOHlSlh9aAtMvq+7dA4oXl38kJub5c2DyZPniOZuIiKwC91YjIkvXr5+8wc+bVzoT2dkB33+vOlWyRo8Gdu6UZYjmzvSL1NBQ4N9/gf/+U50knjdvgLFjAT8/IG1a1WmIiIiMpFkz1QmIiAwrWzbA3V0uz58vxwEDgMGDgbAwdbmSodHI/qkrVwLXr6tO83FMv0i9cEGOJtSo4e1bwMND/hGkS6c6DRERERER6V369MDChbrrc+cCjo7AL7/I9efPpQvw27dq8iXCxkYG0ebOlXWq5sr0i9QnT+R47ZraHP8XGSkfrLi7y37AREREVuXmTdUJiIiM5+uvgYMH4992+bIc27QBMmSQaZUaDVCypNy2aZORQ8Znbw/4+0ux+vy50igfzPSL1Bcv5GgCe6RGRwNeXjLanzu36jREREQKrFunOgERkXHVqiWFwP37QEiILP4EgDNn4j/u4kXg559lNEtxdZg2rRSpY8fKMkVzY/pF6oIFcgwPVxpDq5Xpvd98AxQurDQKERGROia0/IaIyGg0GuCTT+I3o3n5Uqb6Hj0K9OgBlC8vBerFi8CJEzKqOnu2qsRIl07qFw8Pk5qRnCKmX6QWKCDHhg2Vxpg+HWjUCChdWmkMIiIitdjOnohIJ00aoEoV4McfZQ+YyZNl3er16zKqOnSoFLjduwNTpxp9O5ssWXRLFSMjjfrSH8X0Pw49fFh1AixcCJQpI//+iIiIrBqLVCKi5FWsGP/6Tz/JsWZNoHJlo0bJnVv2UfX0lLWqNqY/TGkGI6mKrVoFZM/ObeGIiIgAcLovEVFKfPGFjJpeuABs3Ki7/dNP5VitGpAvn4yyFigge7Lu2mWwOEWKyLLFCRNkGaOpM+0zjVYLRETICVGjMfrLb9kia6RbtTL6SxMREZkmjqQSEaWMjQ1QooR8PX8O/PGHbv/K27elERMA3LolxzZtpDGTgXz+ORAUJMsYR4402MvohWmPpIaGAg4OgIuL0V96717ZW8jNzegvTUREZLpYpBIRpV6mTECLFrq5tr//LqOsf/4pnY0AXaPYJ0+Ax48NEqNaNemxE3f7V1Nk2iOpERFyNPIJ8dgxWffs7m7UlyUiIjJ9LVuqTkBEZP6KF9ddLl1augY7O8v1OXNk/5gYWbIAc+cCHTvq5aW//BLYvFmWNXburJen1DvTHkmNaUFlxCL1/HmZDh6z/RERERHF4eqqOgERkWVxdQW8vYFhw+T65cvx73/2DOjUCZg0SW8LStu0kSWzW7fq5en0zrRHUq9dk6ORdqC9fh1YuRKYMkXJElgiIiIiIrJ2P/8sU39PnZLpnYMGATlzSp+e4GBg3DjZ+iZbNuD1a+DRI6BDB6Bq1VQ1t+vaFZg/X5Y51q9vuD/OhzDtIjU6Wo6OjgZ/qQcPgHnzgGnTzKMtMxERkRKnT8t+BkREZDgODrL/ZZUqQNmyMrr67bfScOnIEeD48fiPX7hQ9mhN5XrFAQPk21xcjL4zznuZdpF6964cy5Uz6Ms8fy7TvqdNYz8IIiKi93r6VHUCIiLrUr26fAGyPjV/ftnKJn164PBh4OxZuS9fPjnu3g00bgzMmgXY2soi1CJFkpwq6u4OjB0rS2JLlTL4nyZFTLtIfftWjv/+a7CXePNG/lImT/5fe/cfVFWd/3H86QVBDEEWkGkH9ZugNiurpO5oMjZtamK1lfJrRTRzd8zd1Vk1d8BfhBVIgTmlrYVNSqgp2zJmmu6uurWbpYZjmrvrsuKPlVrFMXRBVn7ds3+cuBtfEbh44Rzj9Zhp7N5zzue+7wfGt+/zOZ/Px5yvLCIiIi3Q3VwREesEB8PWrU3fq6iAjRvNUVcwF14CWLCg6XmFhZCQcEOT3bqZ+6empsKcORAR4fmw3WXvIjUx0RzWDgrqkOavX4clS8wfSuOWRSIiItICN+Y7iYhIJwgKalqQ/vKXEBYGp0+be7M2Ki83V0s6cwbefhtqasxFmX7wAxzjxpGV9X8sWmQWq9/9bud/jW+yd6bp0aPD5r3U15tD26mp5qi5iIiItIFGUkVE7C021vyv0VtvQV6eOen044/N/Wc2bjQXXQJ4/XUAuvfsSfbuD1iU+QOefdYctLWKvZcIWrrU3MDHw5xOcwT1F7/Q2g8iIiJu0UiqiMjtZcYM+OgjGDHC3J81JATmzYOxY/+3NytAXBx+Vy+w8rl60kNepbLv1+c6HOYzwevXd1rI9s40W7aYK1glJ5uTfj3AMMzHe6dNM+cPi4iIiBs0kioicvsKDjb3ZP2mmhrzseD+/WHgQAIOHCDDN5vFZanksogefL036+zZ5l41d93V4WHaeyQV4MMPPVagAuTmmgtcDRvmsSZFRES6jsBAqyMQERFP8vU1i8/GEbyRIwk9tpfUTyaTNruC+h3v/+/cr77qlJDsPZIK8PjjHmtq3TqzOB0zxmNNioiIdC2dcAddREQs5OsLgwfTF/hFMCx9YxIrr1bi+LLM3P6mE9h/JLVxT6BbVFAAffrAgw96pDkREREREZFvtYEDzZmXz77kjzH4bnNh205g/yLVAxNHt283F0uKi7v1cERERLo0p9PqCEREpBMNGwYTJpjTJjuL/YtUH59bunzvXvjiC3jiCQ/FIyIi0pVdvGh1BCIi0sliYmDoUHP6ZGewd5E6bx5ERrb78oMHobjY3GpGREREPEBb0IiIdEkTJ0JoaIfsEHoDexepCxe2u0j9/HN4/31ITfVwTCIiIl3ZLT7hJCIit6/4eGhoMKdTdiR7F6ntVFoKb70FGRnmvrMiIiLiIRpJFRHp0mbOhPPnYd++jvsMexepzzzj9iVffgmvvAJZWeCw97cTERG5/ahIFRHp8ubNg8OH4dChjmnf3mXcb3/r1umXL5vFaXY2dO/eQTGJiIh0ZUqwIiICpKXBzp1w4oTn27Z3kWoYbT61shLS080i1c+vA2MSERHpyvSYkoiIYE6rXLECNm6E06c927a9M011dZtOu34dliwxOykgoINjEhERERERERwOWLkSXn7ZnHbpsXY911QHePDBVk+pr4fFi81VfENCOiEmERERERERAcxZINnZ5hOtly97pk17F6mtzHtxOmHpUvj5zyE8vJNiEhERERERERc/P7NITU83p2HeKnsXqS3MezEMePZZSE6GgQM7MSYRERERERFpIiDA3AJ08WJzOuatsHeROmTITQ/l5ppPAw8b1onxiIiIiIiISLNCQ81pmGlp5rTM9rJ3kfrUU82+vW4dDB0KY8Z0cjwiIiIiIiJyU337mtMxly41p2e2h72L1GYUFJgV+sSJVkciIiIiIiIi/9+gQTB1qjk9041dRV3sXaSeOtXk5fbt0NAA8fHWhCMiIiIiIiKti46GCRNg1Sr3r7V3kfree67/3bcPyspg5kzrwhEREREREZG2iYmBqCh47TX3rrN3kerlBcChQ/DppzB3rsXxiIiIiIiISJvFxkJICGza1PZr7F2kAidOwK5d5ipRIiIiIiIicnuJj4e6Onj33badb+si9fyQSeTnm/vtdOtmdTQiIiIiIiLSHk8+Cf/8pzmNszUdUqQ6nU7S09NJSkpi+vTpnDt3rsnx/fv3ExcXR1JSEoWFhTdtZ8OuULKywGHrUlpERMT+PJWbRURE2mvePDh82JzO2RLvjvjwvXv3Ultby7Zt2/jss8/Izs5m3bp1ANTV1bFy5Ureeecd/Pz8mDp1Kj/84Q8JDQ11Xd/Q0ADAE4+d5OLF3h0RooiIdCEXLlwA/pdfuiJP5ebGvhQREWmPlBRYvRruvffmublDitQjR44wduxYAKKjozlx4oTrWGlpKf369SMwMBCAESNGUFxczKRJk1znXLp0CYBZzy7tiPBERKSLunTpEv3797c6DEt4KjdPmzatE6MWEZFvq127zD+by80dUqRWVVXh7+/veu3l5UV9fT3e3t5UVVXRq1cv17E77riDqqqqJtdHRUWxefNmQkND8fp6hV8REZH2amho4NKlS0RFRVkdimWUm0VExE5ays0dUqT6+/tz7do112un04m3t3ezx65du9YkMQL06NGDkSNHdkRoIiLSRXXVEdRGys0iImI3N8vNHbIk0fDhw/nTn/4EwGeffcagQYNcxyIiIjh37hxXrlyhtraW4uJi7rnnno4IQ0RERL6m3CwiIreLboZhGJ5u1Ol0kpGRQUlJCYZhkJWVxV//+leqq6tJSkpi//79vPrqqxiGQVxcnOa3iIiIdDDlZhERuV10SJHaVo0J8+9//zs+Pj48//zzTYZ8GxOmt7c3cXFxJCYmWhWqbbTWZzt37iQ/Px8vLy8GDRpERkYGDu3h02q/NVq+fDmBgYEsWrTIgijtpbU+O378ONnZ2RiGQWhoKDk5Ofj6+loYsfVa67MdO3awYcMGHA4HcXFxJCcnWxitvRw7dozc3FwKCgqavK880PmUm92n3Nw+ys3uU252n3Jz+1mamw0L/e53vzNSU1MNwzCMo0ePGnPmzHEdq62tNcaPH29cuXLFqKmpMaZMmWKUl5dbFapttNRn//nPf4xx48YZ1dXVhmEYxoIFC4y9e/daEqfdtNRvjd5++20jMTHRyMnJ6ezwbKmlPnM6ncajjz5qnD171jAMwygsLDRKS0stidNOWvs9i4mJMSoqKoyamhrX329iGHl5ecYjjzxiJCQkNHlfecAays3uU25uH+Vm9yk3u0+5uX2szs2W3sZr63L4Pj4+ruXwu7qW+szHx4etW7fi5+cHQH19fZe/e9aopX4DOHr0KMeOHSMpKcmK8GyppT47c+YMvXv3Jj8/n5SUFK5cucKAAQOsCtU2Wvs9Gzx4MJWVldTW1mIYBt26dbMiTNvp168fa9asueF95QFrKDe7T7m5fZSb3afc7D7l5vaxOjdbWqTebDn8xmOtLYffFbXUZw6Hg5CQEAAKCgqorq4mJibGkjjtpqV+Ky8vZ+3ataSnp1sVni211GcVFRUcPXqU5ORkNmzYwMGDB/nkk0+sCtU2WuozgIEDBxIXF8fDDz/M/fffT0BAgBVh2s7EiRNdq8x+k/KANZSb3afc3D7Kze5TbnafcnP7WJ2bLS1Sb3U5/K6opT5rfP3CCy9w4MAB1qxZo7tBX2up3/bs2UNFRQWzZ88mLy+PnTt3UlRUZFWottFSn/Xu3Zv+/fsTGRlJ9+7dGTt27A13Jruilvrs5MmTfPDBB+zbt4/9+/fz1VdfsXv3bqtCvS0oD1hDudl9ys3to9zsPuVm9yk3e1Zn5QFLi1Qth+++lvoMID09nZqaGn7961+7Hi2SlvttxowZFBUVUVBQwOzZs3nkkUeYMmWKVaHaRkt91rdvX65du8a5c+cAKC4uZuDAgZbEaSct9VmvXr3o0aMHvr6+eHl58Z3vfId///vfVoV6W1AesIZys/uUm9tHudl9ys3uU272rM7KAzeO4XaiCRMmcODAAX784x+7lsN/7733XMvhp6Wl8ZOf/MS1HH5YWJiV4dpCS30WFRXFO++8w8iRI3niiScA8y/5CRMmWBy19Vr7XZMbtdZnmZmZPP300xiGwT333MP9999vdciWa63PkpKSSE5Opnv37vTr14/JkydbHbItKQ9YS7nZfcrN7aPc7D7lZvcpN3tGZ+cBS7egEREREREREfkmbdIlIiIiIiIitqEiVURERERERGxDRaqIiIiIiIjYhopUERERERERsQ0VqSIiIiIiImIblm5BI2JHZWVlPProowwZMsT13qhRo5g7d26z56elpfHQQw9x3333tevzHnjgAe68804cDgeGYdC7d2+ys7Px9/dvcxt5eXmMHj2awYMHs2PHDhISEigqKiIwMJBx48bdclwNDQ1UV1fz3HPP8f3vf/+m12zatImUlJR2fZ6IiIjdHDp0iPnz5xMZGQnAtWvXCA8PJzc3Fx8fn3a3W1ZWxsKFCyksLOSBBx5g9+7d+Pr6eipskdueilSRZkRGRlJQUNBpn/fmm2+6klNOTg5FRUXMmDGjzdfPnj0bMJPeb37zGxISEjyy6fk34/rzn//M2rVref311296/rp161SkiojIt8ro0aNZvXq16/XTTz/N/v37iY2NtTAqkW83FakibdTQ0EB6ejoXLlygoqKC++67j/nz57uOnzlzhsWLF+Pt7Y2XlxcvvvgiYWFhrFq1ik8//RTDMJg5cyaTJk266Wc4nU4qKyu56667qKurY8mSJZw/f56GhgaefPJJHnroITZv3sz27dtxOBwMHz6c1NRU12ju73//e06dOsXatWsxDIOQkBDOnj3L3XffzeTJk7l06RJPPfUURUVFbsUF8OWXXxIQEADAnj172Lx5s+vYyy+/zLZt27h69SoZGRksXbqUZ555hnPnzuF0Opk/fz6jRo26tR+AiIiIxWpraykvLycwMLDZPHrs2DEyMzMxDIOwsDByc3M5fvw4a9euBeD69eu88MILdO/e3eJvImJvKlJFmnHq1CmmT5/uep2bm0tdXR3R0dEkJCRQU1NzQ5H68ccfM2TIENLS0iguLubq1aucPHmSsrIytm7dSk1NDYmJicTExLiKvUazZs3C4XDQrVs3hg4dyuOPP87WrVsJCgoiJyeHqqoqpkyZwujRoykqKmL58uVER0ezZcsW6uvrXe3MmTOHkpIS5s6dy5o1awBITExkxYoVTJ48mXfffZcpU6bw4YcftjmumpoaysvLGTt2LKmpqQCcPXuWvLw8/Pz8SE9P56OPPuJnP/sZmzZtIiMjgy1bthAUFERWVhYVFRWkpKSwa9cuT/+YREREOtzBgweZPn06ly9fxuFwkJiYSG1tbbN5dPny5axevZqIiAg2b95MaWkp//jHP8jJySEsLIzXXnuNPXv28KMf/cjqryViaypSRZrR3OO+VVVVfP755xw8eBB/f39qa2ubHI+Pj2f9+vX89Kc/pVevXixYsICSkhL+8pe/uAre+vr6JiOSjb75WG2j0tJSxowZA4C/vz8RERGcP3+elStX8uabb5Kbm0t0dDSGYbT4XSIiImhoaOCLL77g/fffZ+PGjWzbts2tuF566SXKysoIDg4GIDg4mNTUVO644w5Onz5NdHR0k+tKSko4cuQIx48fd7VfUVFBUFBQi7GKiIjYTePjvhUVFcyaNYvw8PCb5vfLly8TEREBwLRp0wD417/+RWZmJj179uTixYsMHz7csu8icrvQ6r4ibVRUVESvXr1YtWoVs2bN4vr1600KxH379jFixAjy8/OJjY3ljTfeYMCAAYwaNYqCggLy8/OZNGkS4eHhbfq8iIgIiouLAbNALikpITw8nMLCQlasWMGmTZv429/+xtGjR13XOBwOnE7nDW3Fx8eTk5NDZGQkAQEBbsc1f/58ysvL2bJlC5WVlbzyyiusXr2a559/Hl9fX1c/NP45YMAAHn74YQoKCli/fj2xsbEEBga26XuLiIjYUePTTcuWLSMkJKTZPNqnTx/Onj0LmIsa/uEPf2DZsmVkZWWRnZ1Nnz59Wr25LCIaSRVps3vvvZeFCxdy5MgR/Pz86N+/P+Xl5a7jUVFR/OpXv2LNmjU4HA4WL17M9773PQ4fPkxycjLV1dWMHz++zav2JiYmsnz5cqZOnUpNTQ1z584lODiYwYMHEx8fT1BQEGFhYQwbNoyioiLAHOGsq6sjJyeHHj16uNqKjY0lMzOTdevWAebKve7E5XA4yMzMZNq0aYwfP57hw4czefJkevbsSUBAgKsfIiIiWLRoEVlZWSxbtoyUlBSqqqpITk7G4dA9MRERub1FRkYyffp0/vjHP3LnnXfekEdXrFjBkiVLcDgchIaGMnPmTB577DESExMJCAggJCSkyb8dRKR53QzdzhERERERERGb0NCGiIiIiIiI2IaKVBEREREREbENFakiIiIiIiJiGypSRURERERExDZUpIqIiIiIiIhtqEgVERERERER21CRKiIiIiIiIrbxX6ZGq8kOYmbAAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1152x576 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import rpy2.robjects as robjects\n",
    "\n",
    "root_dir = './demo_data/'\n",
    "#! Use raw count data to build simulated doublets\n",
    "def generate_dbl_vector(ct1_idx, ct2_idx, num, mtx):\n",
    "    a = np.random.choice(ct1_idx, num, replace=True)\n",
    "    b = np.random.choice(ct2_idx, num, replace=True)\n",
    "    rand_mtx1 = mtx[a]\n",
    "    rand_mtx2 = mtx[b]\n",
    "    sampled_dbl_mtx = (rand_mtx1+rand_mtx2)\n",
    "    rand = np.random.rand(sampled_dbl_mtx.shape[0], sampled_dbl_mtx.shape[1]) - 0.5\n",
    "    sampled_dbl_mtx = np.round(sampled_dbl_mtx + rand)\n",
    "    return sampled_dbl_mtx\n",
    "roc_collect = {}\n",
    "prc_collect = {}\n",
    "this_resolution = 0.8\n",
    "\n",
    "file = 'example_data.npz'\n",
    "meta_name = 'ground_truth_for_example_data.csv'\n",
    "print(file,'---' ,meta_name)\n",
    "simulated_data = sparse.load_npz(os.path.join(root_dir, file.split('.')[0]+'.npz')).todense()\n",
    "meta_data = pd.read_csv(os.path.join(root_dir, meta_name)).reset_index(drop=True)\n",
    "filtered, adata_hvg, adata_filter_gene = umap_preproc(simulated_data)\n",
    "sc.tl.leiden(filtered, resolution=this_resolution)\n",
    "leiden_lst = filtered.obs['leiden']\n",
    "print(np.unique(leiden_lst).shape)\n",
    "\n",
    "for name, this_mtx in {'hvg_gene': adata_hvg}.items():\n",
    "np_array_transform = np.array(this_mtx.X)\n",
    "print('Now Running: ', np_array_transform.shape)\n",
    "search_lst = leiden_lst\n",
    "\n",
    "matrix = np_array_transform\n",
    "ct_idx = {}\n",
    "# create index dict for each cell type\n",
    "for ct in np.unique(search_lst):\n",
    "    ct_idx[ct] = search_lst[search_lst == ct].index\n",
    "\n",
    "# pre_compute sgl theta and miu\n",
    "nb_sgl_param_dict = {}\n",
    "for this_ct, this_idx in ct_idx.items():\n",
    "    this_mtx = np.array(matrix)[[int(i) for i in this_idx]]\n",
    "    all_miu, all_var, all_theta = NB_pre(this_mtx)\n",
    "    nb_sgl_param_dict[this_ct] = [all_miu, all_var, all_theta]\n",
    "\n",
    "# pre_compute dbl theta and miu\n",
    "doublet_ct_lst = None\n",
    "nb_dbl_param_dict = {}\n",
    "artificial_dbl_clt = np.array([])\n",
    "doublet_ct_lst = list(itertools.combinations(np.unique(search_lst), 2))\n",
    "for ct1, ct2 in doublet_ct_lst:\n",
    "    ct1_dx, ct2_idx = [int(i) for i in ct_idx[ct1]], [int(i) for i in ct_idx[ct2]]\n",
    "    artifitial_cmb_mtx = generate_dbl_vector(ct1_dx, ct2_idx, num=500, mtx = matrix)\n",
    "    if artificial_dbl_clt.shape[0] == 0:\n",
    "        artificial_dbl_clt = artifitial_cmb_mtx\n",
    "    else:\n",
    "        artificial_dbl_clt = np.vstack((artificial_dbl_clt, artifitial_cmb_mtx))\n",
    "    cmb_miu, cmb_var, cmb_theta = NB_pre(artifitial_cmb_mtx)\n",
    "    dbl_type = (ct1+'+'+ct2) if ((ct1) < (ct2)) else (ct2+'+'+ct1)\n",
    "    nb_dbl_param_dict[dbl_type] = [cmb_miu, cmb_var, cmb_theta] \n",
    "\n",
    "whole_lst = list(nb_sgl_param_dict.keys()) + list(nb_dbl_param_dict.keys())\n",
    "scores_df_fast = pd.DataFrame(index=range(matrix.shape[0]))\n",
    "scores_artificial_dbl = pd.DataFrame(index=range(artificial_dbl_clt.shape[0]))\n",
    "\n",
    "for sgl, sgl_param in nb_sgl_param_dict.items():\n",
    "    miu, var, theta = sgl_param\n",
    "    this_score_all = NB_model(matrix, theta, miu)\n",
    "    this_score =  np.mean(this_score_all, axis = 1)\n",
    "    scores_df_fast[sgl] = this_score\n",
    "\n",
    "    this_score_all_artificial = NB_model(artificial_dbl_clt, theta, miu)\n",
    "    this_score_artificial =  np.mean(this_score_all_artificial, axis = 1)\n",
    "    scores_artificial_dbl[sgl] = this_score_artificial\n",
    "\n",
    "for dbl, dbl_param in nb_dbl_param_dict.items():\n",
    "    dbl_miu, dbl_var, dbl_theta = dbl_param\n",
    "    this_score_all = (NB_model(matrix, dbl_theta, dbl_miu))\n",
    "    this_score = np.mean(this_score_all, axis = 1)\n",
    "    scores_df_fast[dbl] = this_score\n",
    "\n",
    "    this_score_all_artificial = (NB_model(artificial_dbl_clt, dbl_theta, dbl_miu))\n",
    "    this_score_artificial = np.mean(this_score_all_artificial, axis = 1)\n",
    "    scores_artificial_dbl[dbl] = this_score_artificial\n",
    "\n",
    "#! New Approach\n",
    "combined_mtx = pd.concat((scores_df_fast, scores_artificial_dbl), axis=0)\n",
    "y_lables = [0]*scores_df_fast.shape[0]+[1]*scores_artificial_dbl.shape[0]\n",
    "real_arti_mtx = np.concatenate((matrix, artificial_dbl_clt), axis=0).T\n",
    "\n",
    "%R -i real_arti_mtx\n",
    "%R sce <- SingleCellExperiment(list(counts = real_arti_mtx))\n",
    "%R sce1 <- scDblFinder(sce, returnType = 'table')\n",
    "return_r_df = robjects.r(\"sce1[which(sce1$type == 'real'), ]\")\n",
    "cy_return_r_df = return_r_df.copy()\n",
    "cy_return_r_df = cy_return_r_df.drop(['type', 'src', 'class', 'include.in.training', 'nearestClass', 'score'], axis=1)\n",
    "cy_return_r_df = (cy_return_r_df-cy_return_r_df.min())/(cy_return_r_df.max() - cy_return_r_df.min())\n",
    "cy_return_r_df = cy_return_r_df.reset_index().drop(['index'], axis=1)\n",
    "cmb_feature_mtx = pd.concat([cy_return_r_df.reset_index(), combined_mtx.reset_index()], ignore_index=True, axis=1).drop([0], axis=1)\n",
    "clf = GradientBoostingClassifier(n_estimators=100, learning_rate=1, max_depth=1, random_state=0).fit(cmb_feature_mtx, y_lables)\n",
    "\n",
    "pred_mtx = cmb_feature_mtx.iloc[:scores_df_fast.shape[0], :]\n",
    "rltss = clf.predict_proba(pred_mtx)\n",
    "y_label = meta_data['x'].copy()\n",
    "y_label[y_label == 'doublet'] = 1\n",
    "y_label[y_label == 'singlet'] = 0\n",
    "y_label = [int(i) for i in y_label]\n",
    "y_pre = rltss[:, 1]\n",
    "\n",
    "pd.DataFrame({'score':y_pre}).to_csv('./demo_data/scMODD_predicted_score_combined_feature.csv')\n",
    "\n",
    "roc_pr_plot(y_label=y_label, y_pre=rltss[:, 1], title=meta_name)\n",
    "print(('FINISH!!!---', this_resolution))\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.12 ('base')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "e88ed867cd18119a9b2c02d720f6306e9a53bfdbbd7b446bd3d11a1671cb0173"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}