mnn_correction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "badc38dd-4b45-41ce-83c4-29369517fc9d",
   "metadata": {},
   "source": [
    "## CCA Seurat impl"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc9a59ba-6b53-46e1-9871-8c9a183be125",
   "metadata": {},
   "source": [
    "#### why does ols fail when we have large number of variables and they are highly correlated "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67f9c9ad-79dd-4ebe-89c7-935332195c6d",
   "metadata": {},
   "source": [
    "models the cross covariance of $X^TY$ using latent variables. useful when the columns of $X$ and $Y$ are collinear. We basically compute the latent space projection vectors of $X$ and $Y$ by maximizing cross covariance of $X^TY$. This way we choose latent vectors that maximize this quantity and are maximally related by the loss function. So it has greater predictive power as it will retain those eigen vectors that represent the maximal. \n",
    "\n",
    "So we basically compute the svd of $cov(X^TY)$ and then choose the maximal eigen vector for both. This is the best rank 1 approximation of $cov(X^TY)$. So we choose ξ and ω that maximize Cov(Xu , Yv) and then we find the residual matrix by saying that X can be approximated by Xu. and we say that X_new = coeff * Xu_r and find the residual matrix by estimating coef since we already know Xu_r. Then we subtract X_new from X to obtain the residual matrix and we do the similar for Y.\n",
    "\n",
    "\n",
    "So basically we approximate the X and Y matrices by the eigen vectors that best explain the cross covariance of $X^TY$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d83a394c-7887-4d58-bad1-40261c70ba58",
   "metadata": {},
   "source": [
    "### here i implement the Mutual Nearest Neighbor correction as described in laleh hagverdhi et al. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e0d62f8-244c-4a97-88ec-c92389b71c85",
   "metadata": {},
   "source": [
    "basically computes the mutual nearest neighbors between two datasets. Then says that the difference between these mutual nearest neighbors is the batch effect. So we get a reduced set of batch effect for every unique cell in the query dataset. Then we say that the batch effect for each cell in the query is a gaussian kernel weighted sum of the batch effects of the mutual nearest neighbor cells so for every query we compute a gaussian weighted G(x, B) where B is the unqiue mnn cell in the query dataset and do G(x, B) * batch_vector / sum_over_B(G(x, B)). Then we subtract off this batch effect computed for every cell to get the corrected data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 386,
   "id": "be3c7377-89af-458f-a7e9-9ee067ef65f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cross_decomposition import CCA\n",
    "from sklearn.cross_decomposition import PLSRegression\n",
    "from sklearn.cross_decomposition import PLSCanonical\n",
    "import scanpy as sc\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "np.set_printoptions(suppress = True)\n",
    "\n",
    "from sklearn.svm import SVR\n",
    "from sklearn.model_selection import StratifiedShuffleSplit\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from tqdm import tqdm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 387,
   "id": "c065922b-3b9e-441b-b321-33d65d64fffe",
   "metadata": {},
   "outputs": [],
   "source": [
    "eps = 1e-20\n",
    "nhvg = 2000\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 388,
   "id": "c999ff47-f018-4651-a1f9-1e6eee91140e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AnnData object with n_obs × n_vars = 13999 × 14053\n",
       "    obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'stim', 'seurat_annotations'\n",
       "    var: 'features'"
      ]
     },
     "execution_count": 388,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adata = sc.read_h5ad(\"ifnb.h5ad\")\n",
    "n_cells, n_genes = adata.shape\n",
    "adata ## axis 1 is genes "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 389,
   "id": "fc300d75-9fd7-40b1-be00-f00345c09874",
   "metadata": {},
   "outputs": [],
   "source": [
    "#sc.pp.normalize_total(adata, 10000)\n",
    "#sc.pp.log1p(adata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 390,
   "id": "ea9a7ac9-1a41-4aab-abd6-41590c81451d",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_ctrl = adata[adata.obs[\"stim\"]==\"CTRL\"]\n",
    "adata_stim = adata[adata.obs[\"stim\"]==\"STIM\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "803d4aa0-3637-4867-bc25-2b22769da766",
   "metadata": {},
   "source": [
    "### choose highly variable features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 391,
   "id": "d0c189c5-c8f7-4900-88c1-41cb604e0d4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def filter_count_matrix(cnt, axis, flt_val, flt_name):\n",
    "    print (f\"Executing {flt_name} on axis {axis} by val {flt_val}\")\n",
    "    mask_vec = cnt.sum(axis = axis)\n",
    "    print (f\"shape of mask vec {mask_vec.shape}\")\n",
    "    mask = mask_vec > flt_val\n",
    "    \n",
    "    if axis == 0:\n",
    "        cnt_flt = cnt[:, mask]\n",
    "    else:\n",
    "        cnt_flt = cnt[mask, :]\n",
    "    print (f\"shape before {cnt.shape}, shape after {cnt_flt.shape}\")\n",
    "    return cnt_flt, mask\n",
    "\n",
    "\n",
    "def compute_z_scores(raw_counts, std_estimated):\n",
    "\n",
    "\n",
    "    n_cells = raw_counts.shape[0]\n",
    "    raw_counts_mean = raw_counts.mean(axis = 0)[None, :]\n",
    "    raw_counts_std = np.sqrt(10**(std_estimated))[None, :]\n",
    "    z = (raw_counts - raw_counts_mean)/raw_counts_std\n",
    "    z = np.clip(z, 0, np.sqrt(n_cells))\n",
    "    return z\n",
    "\n",
    "def fit_log_mean_log_var(log_mean_genes, log_var_genes):\n",
    "    \"\"\"\n",
    "    estimating best hyperparameters over the given grid \n",
    "    \"\"\"\n",
    "    svr = GridSearchCV(SVR(kernel=\"rbf\", gamma=0.1), param_grid={\"C\": [1e-2, 1e-1, 1e0, 1e1], \"gamma\": np.logspace(-2, 2, 4)})\n",
    "    svr.fit(log_mean_genes, log_var_genes)\n",
    "    \n",
    "    svr_best = SVR(kernel = \"rbf\", **svr.best_params_)\n",
    "    svr_best.fit(log_mean_genes, log_var_genes)\n",
    "    print (svr.best_params_, svr_best.score(log_mean_genes, log_var_genes))\n",
    "    log_var_smoothed = svr_best.predict(log_mean_genes)\n",
    "    return log_var_smoothed\n",
    "\n",
    "def filter_adata_and_get_gene_stats(adata):\n",
    "    raw_counts = adata.X.toarray()\n",
    "    raw_counts, low_genes_mask = filter_count_matrix(raw_counts,  0, 3, \"filter low genes\")\n",
    "    raw_counts, low_total_count_mask = filter_count_matrix(raw_counts, 1, 3, \"filter low total counts\")\n",
    "    adata_flt = adata[:,low_genes_mask]\n",
    "    adata_flt = adata_flt[low_total_count_mask, :]\n",
    "    log_mean_genes = np.log10(raw_counts.mean(axis = 0))[:, None]\n",
    "    log_var_genes = np.log10(raw_counts.var(axis = 0))\n",
    "    \n",
    "    return log_mean_genes, log_var_genes, adata_flt, adata, raw_counts\n",
    "\n",
    "def estimate_hvg(adata):\n",
    "    log_mean_genes, log_var_genes, adata_flt, adata, raw_counts = filter_adata_and_get_gene_stats(adata)\n",
    "    log_var_smoothed = fit_log_mean_log_var(log_mean_genes, log_var_genes)\n",
    "    z = compute_z_scores(raw_counts, log_var_smoothed)\n",
    "    variance_genes = z.var(axis = 0)\n",
    "    \n",
    "    hvg_idx = np.argsort(variance_genes)[::-1]\n",
    "    hvg_vals = variance_genes[hvg_idx]\n",
    "    hvg = adata_flt.var[\"features\"].iloc[hvg_idx].to_list()\n",
    "    \n",
    "    adata_flt.var[\"hvg_sorted\"] = hvg\n",
    "    adata_flt.var[\"hvg_variance\"] = hvg_vals\n",
    "    adata_flt.layers[\"z_scores\"] = z\n",
    "    adata_flt.var[\"hvg_idx_sorted\"] = hvg_idx\n",
    "    \n",
    "    return adata_flt, hvg, log_var_genes, log_var_smoothed, log_mean_genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 392,
   "id": "768b82a0-a79e-4521-9e67-e33cb2eba14f",
   "metadata": {},
   "outputs": [],
   "source": [
    "top10_ctrl = [\"HBB\", \"HBA2\", \"HBA1\", \"CCL3\", \"CCL4\", \"CXCL10\", \"TXN\", \"CCL7\", \"CCL2\", \"GNLY\" ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 393,
   "id": "6dff6738-adac-4edb-9eff-2bdb4dbc1290",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Executing filter low genes on axis 0 by val 3\n",
      "shape of mask vec (14053,)\n",
      "shape before (6548, 14053), shape after (6548, 13110)\n",
      "Executing filter low total counts on axis 1 by val 3\n",
      "shape of mask vec (6548,)\n",
      "shape before (6548, 13110), shape after (6548, 13110)\n",
      "{'C': 10.0, 'gamma': 0.21544346900318834} 0.976500917028416\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_394191/2342671033.py:59: ImplicitModificationWarning: Trying to modify attribute `.var` of view, initializing view as actual.\n",
      "  adata_flt.var[\"hvg_sorted\"] = hvg\n",
      "/home/user-dnf/anaconda3/envs/hiplex/lib/python3.8/contextlib.py:120: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.\n",
      "  next(self.gen)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Executing filter low genes on axis 0 by val 3\n",
      "shape of mask vec (14053,)\n",
      "shape before (7451, 14053), shape after (7451, 12971)\n",
      "Executing filter low total counts on axis 1 by val 3\n",
      "shape of mask vec (7451,)\n",
      "shape before (7451, 12971), shape after (7451, 12971)\n",
      "{'C': 10.0, 'gamma': 0.21544346900318834} 0.9804551129458917\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_394191/2342671033.py:59: ImplicitModificationWarning: Trying to modify attribute `.var` of view, initializing view as actual.\n",
      "  adata_flt.var[\"hvg_sorted\"] = hvg\n",
      "/home/user-dnf/anaconda3/envs/hiplex/lib/python3.8/contextlib.py:120: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.\n",
      "  next(self.gen)\n"
     ]
    }
   ],
   "source": [
    "adata_ctrl_flt, hvg_ctrl, lvg_c, lvs_c, lmg_c = estimate_hvg(adata_ctrl)\n",
    "adata_stim_flt, hvg_stim, lvg_s, lvs_s, lmg_s = estimate_hvg(adata_stim)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cf0387e-0603-4dc5-9d40-5802a2e63035",
   "metadata": {},
   "source": [
    "### plot the mean variance relationship and the fitted curves "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 405,
   "id": "32835094-d0ed-455d-9747-b6903398790e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA64AAAHwCAYAAACmHTLcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzde5zVZbn///elKQrMMEtABSZg0lRwaRRYexD8SW01mfZu5yFT/OYpRzMTdhme2mkHd1qZklSKeWpLSSq61fGAe2+MIaYMbNQlaKmjNQMq6AdnQEHF+/fHvT6sNYu1ZtYw6zjr9Xw8fMxhfdbnc4898lrXfV/3dZtzTgAAAAAAlKpdij0AAAAAAAB6QuIKAAAAAChpJK4AAAAAgJJG4goAAAAAKGkkrgAAAACAkkbiCgAAAAAoaSSuQAkxs9vM7AdZXjvezJyZfSjf4wIAYKAys4fN7LQ83fsKM7sjH/cGKg2JK1AAZna6mS0v4vOPNLP2gfIcAAB2RrpE0jl3rHPu9mKNKdSXyetyeA6QaySuQIkws12LPQYAAACgFJG4AjlkZh82s8Vmtt7M3jCz+WY2QdINkurNbJOZbYxfe5uZ/dLMHjKzzZJm9HLvXc3sJ2a2wcxektSQ8voZZrbGzLrM7CUzOyf++yGSHpY0Ov78TWY22sw+aWYtZrbRzNbFx7p7/D1mZtea2etm9paZPW1m0fhrg+Lj+LuZvWZmN5jZnpmek9t/wwAA9M7MLjKzjnhMfN7MPmNmn5V0qaST4jHqqfi1j5vZV+Lfn25mf4jHwI3xeDo1/vt/xONixrJiM6szs9/Hn/uYpBEpr99lZq/GY+syMzs4/vtGSbMkzY2P7YH47y82sxfj91ttZl9Iutf+8We9Ff9ssCjptYPM7DEzezP+93+xp+cA5YDEFciR+Irpg5JekTRe0hhJdzrn1kg6V1KLc26oc64m6W2nSLpSUpWk3kqJz5b0OUkflzRF0gkpr78ef71a0hmSrjWzTzjnNks6VtLa+POHOufWStom6d/lg2q9pM9IOi9+r6MlHSHpAEk1kk6S9Eb8tavjv58kaf/43/mdHp4DAEDBmNmBks6XdJhzrkrSMZJeds49Iuk/JS2Kx6iPZbjFpyQ9LWm4pN9IulPSYfIx71RJ881saIb3/kbSKvnY+n1JqUnuw5I+KmlvSU9KWihJzrkF8e9/FB/bv8Svf1HSdEnDJH1X0h1mNir+2vclLZEUkVQr6fr43z9E0mPxsewt6WRJvzCzg3t4DlDySFyB3PmkpNGSvuWc2+yc2+Kc6y0Z/W/n3B+ccx8457b0cu0XJV3nnPuHc+5NST9MftE51+Sce9F5v5cPZtMz3cw5t8o590fn3PvOuZcl3Sjp/4u//J58Mn2QJHPOrXHOrTMzk0+g/90596Zzrkv+Q8CXehk7AACFsk3SIEkTzWw359zLzrkX+/D+Nufcrc65bZIWSfqwpO8557Y655ZIelc+ie3GzMbKJ7j/Eb92maRuK5rOuVucc13Oua2SrpD0MTMblmkgzrm7nHNr458TFkn6m/znDcnH6nGSRqd85vicfKJ+azzGPynpHu044Q2UFRJXIHc+LOkV59z7fXjPP/pw7eiU619JftHMjjWzP8bLgjZKmqmUEqWU6w8wswfjJUud8gnoCElyzv2fpPmSfi7pNTNbYGbVkkZKGixpVbyEaqOkR+K/BwCg6JxzL0iaI58Yvm5md/Zx68prSd+/E79n6u/SrbiOlhTEK5BC22N1fMvPVfHS305JL8df6ilWf9nMWpNibjTp+rmSTNITZvasmZ0Z//04SZ8K3xN/3yxJ+/b8ZwOljcQVyJ1/SBpr6Y+ncRnek+n36ayTT45DY8NvzGyQ/GzqTyTtEy9Hfkg+oGV6zi8lPSfpo865avl9P+H1cs79zDk3WdLB8qXB35K0QT5gH+ycq4n/M8w5Fwbwvvw9AADkhXPuN865afJJnJPf5iLlN06tkxSJl+qGxiZ9f4qkz0v6Z/nS3/Hx36eN1WY2TtJN8mXPw+OxPRZe75x71Tl3tnNutKRz5MuB95f/PPL7pDhdEy8L/mq65wDlgsQVyJ0n5IPWVWY2xMz2MLPD46+9Jqk2bH60k34n6QIzqzWziKSLk17bXb4sar2k983sWPl9qqHXJA1PKUeqktQpaZOZHSQpDGgys8PM7FNmtpukzZK2SNrmnPtAPohea2Z7x68dY2bH9PAcAAAKxswONLNPxyd1t8hPuG6Lv/yapPFmlvPPwM65VyStlPRdM9vdzKZJSt5DWiVpq3zPiMHylU7JXpP0kaSfh8gnmesl34RRfsVV8Z9PNLPa+I9B/Npt8v02DjCz/2dmu8X/Ocx8s8h0zwHKAokrkCPxvTD/Ir/v5e+S2uWbGknS/0l6VtKrZrZhJx9xk6RHJT0l39BhcdKzuyRdIJ/cBvKzuvcnvf6cpN9KeileNjRa0oXx67ri997ejVC+wdNN8Xu9Ih9kfxJ/7SJJL0j6Y7zU6X8kHdjDcwAAKKRBkq6SrxJ6Vb5B0aXx1+6Kf33DzJ7Mw7NPkW/u9KakyyX9Oum1X8vH1A5JqyX9MeW9N8vvy91oZvc551ZLukZSi3yyeYikPyRdf5ikP5nZJvmYP9s51xb/THC0fP+JtfL/Dq6W//eyw3Ny82cD+WfOUS0AAAAAAChdrLgCAAAAAEoaiSsAAAAAoKSRuAIAAAAAShqJKwAAAACgpJG4AgAAAABK2oeK+XAze1n+KI5tkt53zk3p6foRI0a48ePHF2BkAIBKsGrVqg3OuZHFHkc5IzYDAHIpU2wuauIaN8M5l9W5luPHj9fKlSvzPR4AQIUws1eKPYZyR2wGAORSpthMqTAAAAAAoKQVO3F1kpaY2SozayzyWAAAAAAAJajYpcKHO+fWmtnekh4zs+ecc8uSL4gntI2SNHbs2GKMEQAAAABQREVdcXXOrY1/fV3SvZI+meaaBc65Kc65KSNH0j8DAAAAACpN0RJXMxtiZlXh95KOlhQr1ngAAAAAAKWpmKXC+0i618zCcfzGOfdIEccDAAAAAChBRUtcnXMvSfpYsZ4PAAAAACgPxe4qDAAAAABAj0hcAQAAAAAljcQVAAAAAFDSSFwBAAAAACWNxBUAAAAAUNJIXAEAAAAAJY3EFQAAAABQ0khcAQAAAAAljcQVAAAAAFDSSFwBAAAAACWNxBUAUD6CQGpu9l8BAEDxFSg2k7gCAMpHLCYtXeq/AgCA4itQbP5QXu8OAEAuRaPdvwIAgOIqUGwmcQUAlI9IRJo+vdijAAAAoQLFZkqFAQAAAAAljcQVAAAAAFDSSFwBAAAAACWNxBUAAAAAUNJIXAEAQDdmtquZ/cXMHiz2WAAgJzgHvOyRuAIAgFSzJa0p9iAAIGc4B7zskbgCAIDtzKxWUoOkXxV7LACQM9GoNGMG54BnUgYr0iSuAAAg2XWS5kr6oNgDAYCcCc8ajUSKPZLSVAYr0h8q9gAAAEBpMLPPSXrdObfKzI7s4bpGSY2SNHbs2AKNDgCQN+FKdAmvSLPiCgAAQodL+lcze1nSnZI+bWZ3pF7knFvgnJvinJsycuTIQo8RAJBrZbAiTeIKAAAkSc65S5xztc658ZK+JOn/nHOnFnlYAACQuAIAAAAASht7XAEAwA6cc49LerzIwwAAQBIrrgAAAACAEkfiCgAAAAAoaSSuAAAAAICSRuIKAAAAAChpJK4AAAAAgJJG4goAAAAAKGkkrgAAAACAkkbiCgDInSCQmpv9VwAAgBwhcQUA5E4sJi1d6r8CAADkyIeKPQAAwAASjXb/CgAAkAMkrgCA3IlEpOnTiz0KAAAwwFAqDAAAAAAoaSSuAICyEQSBmpubFdD8CQCAvivjJookrgCAshGLxbR06VLFaP4EAEDflXETxaLvcTWzXSWtlNThnPtcsccDAChd0XjTpyjNnwAA6LsybqJY9MRV0mxJayRVF3sgAIDSFolENJ3mTwAA7JwybqJY1FJhM6uV1CDpV8UcBwAAAACgdBV7j+t1kuZK+qDI4wAAAAAAlKiiJa5m9jlJrzvnVvVyXaOZrTSzlevXry/Q6ABUrDLutgcAADBQFXPF9XBJ/2pmL0u6U9KnzeyO1Iuccwucc1Occ1NGjhxZ6DECqDRl3G0PAABgoCpacybn3CWSLpEkMztS0oXOuVOLNR4AkFTW3fYAAAAGqlLoKgwApaOMu+0BAAAMVMVuziRJcs49zhmuAAAAAFBegiBQc3Ozgjz3B2HFFQAAAABKWRD4/hvRqK8OKwFBECgWi6mrq0srV66UpLyetV4SK64AAAAAgAxKsHlkLBbT0qVL5ZzTjBkzFM1zfxBWXAEAAACglJVg88gwUY1Go4oUYBWYxBUAAAAASlkJNo+MRCJ5LQ1ORakwAAAAAKCkkbgCAAAAAPomCKTmZv+1AEhcAQAAAAA76PGomwI3jGKPKwAAAABgB2HnYCnNUTcFbhhF4goAAAAMRCV49ifKS3Ln4B0UuGEUiSsAAAAwEIWlnFLJdaRFeSh05+CekLgCAAAAA1EJnv0J7CwSVwAAAGAgKsGzP4GdRVdhAAAAAKhgGbsHF/jIm56QuAIAAABABQu7B8dSj7Yp8JE3PaFUGAAAAAAqWMbuwSW0T5rEFQCQF0EQKBaLKRqNKsIxDAAAlKyM3YNLaJ80pcIAgLzIWHYEAADQR6y4AgDyosdDywEAQGEEgd+jGo36FdQyxYorACAvwrIjyoQBABWpVDryJjVYytg9uAyQuAIAcqacAyIAADlVKh15o1FpxgwpGi3rbTyUCgMAciYMiJLSN3kAAKASBIHU2SlNmVL8jrxJDZbKeRsPiSsAIGfKOSBCMrM9JC2TNEj+M8LdzrnLizsqAChDsZi0apVf6SyhLTMZuweXARJXAEDOlHNAhCRpq6RPO+c2mdlukpab2cPOuT8We2AAUFZK6PzTgXI8HXtcAQCAJMl5m+I/7hb/xxVxSABQnsLy3BJIFMt5X2syVlwBAMB2ZrarpFWS9pf0c+fcn4o8JABAPwyUbTysuAIAgO2cc9ucc5Mk1Ur6pJnt8EnHzBrNbKWZrVy/fn3hBwkAyNpAOZ6OxBUAkB+lcn4ddopzbqOkxyV9Ns1rC5xzU5xzU0aOHFnwsQEoQfw3H3lG4goA2FEuPoCUyvl1yJqZjTSzmvj3e0r6Z0nPFXdUAMoC/81HnrHHFQCwo/ADiLT97Lc+K6GOisjaKEm3x/e57iLpd865B4s8JgDlgP/mI89IXAEAO8rFB5CkA89RHpxzT0v6eLHHAaAM8d985BmJKwBgR3wAAQAAJYQ9rgAAAACAkkbiCgAAAAAoaSSuAAAAAICSRuIKAAAAAChpJK4AAAAAgJJG4goAAAAAKGkkrgAAAACAkkbiCgBIKwgCNTc3KwiCYg8FAABUOBJXAEBasVhMS5cuVSwWK/ZQAACoaEwmk7gCANIIgkBdXV2aPHmyotFosYcDAEBFYzJZ+lCxHmxme0haJmlQfBx3O+cuL9Z4AAAJsVhMK1eu1IwZMxSJRIo9HAAAKlo4iVzJk8lFS1wlbZX0aefcJjPbTdJyM3vYOffHIo4JACACJAAApSQSiWj69OnFHkZRFa1U2Hmb4j/uFv/HFWs8AIosCKTmZv8VRReJRDQ9GlUkFuN/EwAAUHRF3eNqZruaWauk1yU95pz7UzHHA6CIYjFp6VL/FaWB/00AAECJKGapsJxz2yRNMrMaSfeaWdQ51+0Tkpk1SmqUpLFjxxZhlAAKIixJpTS1dPC/CQAAKBEl0VXYObdR0uOSPpvmtQXOuSnOuSkjR44s+NgAFEgkIk2f7r+iNPC/CQAAKBFFS1zNbGR8pVVmtqekf5b0XLHGAwAAAAAoTcUsFR4l6XYz21U+gf6dc+7BIo4HAAAAAFCCipa4OueelvTxYj0fAAAAwAATBL6pYDTav60uuboPcqYk9rgCAAAAQL/lqiM+nfVLTlG7CgMAMmCmFwCAvstVR3w665ccElcAKEXhTK/kO/sCAIDehR3xS+U+yBkSVwAoRcz0AgAAbMceVwAoRZyhCgDItyCQmpv9V6DEkbgCAAAAlYgGRCgjlAoDAAAAlYhtKSgjJK4AAABAJaIBEcoIpcIAAAAAylN/9umyx7eskLgCAAAAKE/92adbwD2+QRCoublZAUnyTqNUGAAAAEB56s8+3QLu8Y3FYloaP599OuXZO4XEFQAAAEB56s8+3dT3BoFffa2tldrbfUKbo2PpovHkOEojrJ1G4goAAACgsoRJanJyGpYO77WX9Oab/nc5Wh2NRCKstPYTiSsAAACAyhImqVIiOQ1XQ5NXXPshCALFYjFFo1FFcrRyW8lIXAEAAABUlnT7W5NLh+vq+v0I9rXmFokrAORCupIjAADQP/mKrwU4w5Z9rbnFcTgAkAsFbKkPAEDFKNP4Splw7rHiCqD/WG0saEt9AAAqRpnGV8qEc4/EFUD/pWtwUGkKUHIEAEDFKVZ87eekPGXCuUfiCqD/ynQ2FAAAlJBSquDq56Q8x9/kHokrgP5jtREAAPRXmCx2dkrV1cVNYJmULzkkrgAAAACKL0wSu7qKvwWJSfmSQ+IKAMVQSuVQAACUgjBZDAKpqio/q53E37LFcTgAUAxl2t4fAIC8CxPYfCSWOY6/QRCoublZQRD0dqHU3Oy/YqeQuAKoLKUSOKJRacYM9s6gpJjZh81sqZmtMbNnzWx2sccEADmV4/jb0tKiW2+9VS0tLT1fyIR1v1EqDKCylMrRPeydQWl6X9I3nXNPmlmVpFVm9phzbnWxBwYAOZGj+NvW1qbf/e53evPNN7V161Y553p+A82e+o3EFUBlIXAAGTnn1klaF/++y8zWSBojicQVAJI0NTXpzjvv1Ihhw/T5SZM0deLEnt/AhHW/kbgCqCwEDiArZjZe0scl/am4IwGA0tPQ0KDNmzdr/Dvv6OitWxVpb5fq6oo9rAGNPa4AAKAbMxsq6R5Jc5xznWlebzSzlWa2cv369YUfIAAUSdiMqaamRhc1NuqkCRMUqamRamuLPbQBjxVXAACwnZntJp+0LnTOLU53jXNugaQFkjRlypReNnYBwMARi8W0NN4rY7okPfqof4EV17wjcQWATDjrDRXGzEzSzZLWOOd+WuzxAECpicZ7ZIRfdcIJknP0zigASoUBIBNa16PyHC7p/0n6tJm1xv+ZWexBAUDO9fF4vLa2Ns2fP18bN27U9OnTFZH854P6eqmhgQnuAiBxBYBMOGsVFcY5t9w5Z865Q51zk+L/PFTscQGoEIU8az0Wkx56SFq4MKvnNTU1afHixWpqavK/aGmRbr3Vf0VBUCoMAJnQgRgAgMIp5Fnr0aj01FNSR4d/bobnBUGgWCymadOmSfLdhCX58uDkr8g7ElcAAAAAxVfIs9YjEWnWrEQviwxaWlp0991364QTTtD555+feGHqVKm6mqqsAiJxBYASFM7wRqNRRdg3AwCoBIWudMriea6rS3rtNf+1j+9FbrHHFQBKUNhuP0ZjKAAA+q8v+2eDQE1XXaUjDz9cm1av1hn77KOpVVV9vw9yihVXAChBO7TbBwAAO68P+2eDlhZddu21WvPmm5Kkx6+6KlESXMh9uOiGxBUASlAkEtF0AiIAALnRh/2zdzzzjN7ctk21I0boWxdc0D1BLeQ+XHRDqTAAAACA8rCzpbrhntRe+kYEQaA7779fb2zapGhVlRpGj96p+yD3WHEFiiEIEl3s+A8fAABAdnamVDfLz11BEGjhwoUaP368dpP0zZNOYmW1hBRtxdXMPmxmS81sjZk9a2azizUWoODC/+jSeAcAACB7tbXSXntJVVXpV17Trcim+9yV5rpYLKaOjg7V19fr3gcf1BEXXOBfoBlTSSjmiuv7kr7pnHvSzKokrTKzx5xzq4s4JqAw2B+BOI69AQBA2VejtbdLb74pLV/uv0rdV17TrchGo1JXl9TZ6Z8TiaS9Lrkx4vaYTDOmklG0xNU5t07Suvj3XWa2RtIYSSSuGPg4+wtx4bE3kmjGBACoXNkmiOGkf22tT2JTFwHSLQ5EIn6FdulSqbra3z/NddsbIwaB9NBDknPSxInSjBksNpSAktjjambjJX1c0p+KOxIAKJwgCNTZ2akpU6Zw7A0AoLJlW42WPPlfV9f9tSCQWlp8wrlxY/cV3NT7x+/T1tampoUL1dDQoLrwfrGYdPfd/vsw0Q1Li+lPUjRFT1zNbKikeyTNcc51pnm9UVKjJI0dO7bAowMwoKWWJYU/J8/i5ik4BUGghQsWqOOZZzTz5JMpEwYAVLZcVKMlJ5xtbd1LidPcPwgCzZs3T6tWrZIknX/++f6FaFQ65hjp5Zf9Z4Lw3pQMF1VRE1cz200+aV3onFuc7hrn3AJJCyRpypQproDDA1As/em63Jf3pgah8Oe99kq/byaHYrGYOp5+WmOCQFGzvDwDAICKEo1KJ5yQKPFNV0qcZPHixXr88cc1ceJENTQ0JF6IRKTRo6Xnn/f3qKujP0kJKFriamYm6WZJa5xzPy3WOACUoP7MavblvalBqLd9MzkUjUalU05R1EyR+vq8PQcAgIoRiUgzZyZ+Ti0lTrJs2TL94Ac/0FtvvaX6+vpEmXAoQ2kxiqeYK66HS/p/kp4xs9b47y51zj1UxDEBKAX9mdXsy3tTg1BP+2Z2VsoKcHIX4ekNDZzpCwBAEVx//fV64403tPfee+ucc87p/iKxuSQV7RxX59xy55w55w51zk2K/0PSCiCRQO5MsOjPe/MhFvOdCRculIJAsRUrtPTWWxVbsSLxOmf6AgAqRbpzVgt4jyAI1NTUpMmTJ2vGjBm65ZZbNGnSpO4XZXnuKwqr6M2ZAKBgijGDGo1KTz2l4NlnFfvud1V7yCGaISX2tbJnBgBQSXLR5Kgf94jFYrrnnnskSRdeeKGmR6M7dgtOF5tpzlR0JK4AKkcxgk4koqChQQt//3t1tLdr5ujRmn7GGeyZAYCBoNJKSnPx9/ZlwjbT83Zi0retrU233XabtmzZosMPP1yjRo1StLbWV0R1dPiLwnicLjYz0Vx0JK4AKkeRgk7LmjVqkfSxSZMUPfHE3O2fBQAUV6WtwuXi7+3LhG2m5yXfI8tk+sYbb9QNN9yg6upqzZ07VzNnzvQrrR0d0pgxfTs/FkVB4gqgchQ66ASB2pYs0f2PPKJ3JR38uc8pQtIKAANHpa3C5frv7S3pzOZ5WSbTr732mpxzGjVqVOLom+T7Jz+/0lbSywSJKwDkSdDSonk//ame2rRJnzrqKE2dOrXYQwIA5FKlrcLl+u/tLenM5nm9JLdBEGjJkiWqrq7Wcccdp9mzZyeOvsl0/0pbSS8TJK4AkCcx52Rjx+pTw4dr9uzZijBrCwBAQi5WcHtIboMg0MKFC3X//fdrw4YNOvPMM7t3EA4CacUKyUyqr/f3CgKpq0uaPLlyVtLLBIkrAORYW1ubmpqaNG3aNB13+umKmilSU1PsYQEAUFryvGIdi8XU0dGhqVOnarBzahg50iem4URyLCbFOwyrqsqPJRaTVq6UZsygTLjEkLgCQI7ddddd+u1vf6u3335bc+vrfblRGBABAKgkudwv2sd7ReMrptFoVJHwXPX166WGBmnNGr+yeswxPkanrv6y2lpySFwB5M5AbmbQh79t3LhxGjVqlMaNG+ev7+z0wTF5lhcAgErQ3/2iyfG3t3vFrw1qa7Vi9WqZmerr6/1Wnfi56urokJqapCef9O8544zMHYtRUkhcAeTOQG5m0Ie/7eijj9bo0aP9TG8kIlVXs+oKAKhMPa1gZjMpnBx/e1sNjcUUPPSQFr7zjlrWr9cgSVUvvaTps2b5+8+a5e9XWyt95COSc6yslhESVwC5M5DLa3r528J9rQ0NDaqrq9P05AR1IP97AQCgJz2duZrNpHDqkTU9TQBHo1pwxx36+YMP6tNHHaUTJk5UtKPDP2f69O7v53i6skPiCiB38l1eU8xS5F7+tt/97ne68847tXnzZl3U2Nh9nJQdAQCQSFS7unwVUm2tb4IUJqfp4nwfYmjrK6/o2v/+b73xxhv629/+ppnXXpu4H8reLsUeAABkLQx4sVixR7JdEARqbm7WiBEjNGrUKI0fP15qaZFuvdV/BQAAXjTqE1XnfDxvb0+shErd43wQSM3N/muqDK/dcMMN2rJli4YPH65LL700kfSGx9w89JDf35runih5rLgCKB8lWHK7ZMkS/frXv9bxn/2sLpk5U9FPflJavdq/6FxxBwcAQCkJE8kg8P0fUuN5cpzvqYw45bUgCBSLxTRx4kQdeOCBamxsVENDw47vuftu/311NZVQZYjEFUB5KFSZcLbPiV/3yurVWrdunTY8+aTONPN7ZqZOTR+QAQCoZMkxNl3imFwWnGmyOgh8qfHkydtf2z6JfPzx+slPfrL9GJxuolHphBNoyFTGKBUGUB4KUSYcBNLChb6UqJfntC5erDkXXqiD99xTZ555pk6cNs2/YNa9NAkAAHh9ieWZYmksJq1c6SeIIxEtuvlmzfnqV/XsM89ow4YNmj59uj/+Jt39Zs70Z7gSn8sSK64AykNPM6/JK6T9WZmNxfz5bmPG9DgbGwSBrrz/fv2prU16+WVdd8MN/rkHHMAsLgAAqcLYnNqMaWckfR4IgkCXfec7ejUINHqXXXTiiSemf/aKFX5iub6epLWMseIKoO96apiQLz3NvCbP3vY2m9vT2KNRadIk6e23pY0bMw6lpaVF70r6eH29Tj/lFH8/iVVWAADSCZsWrlnT/1gZ/zzQtnGjLr30Uo3Yd1+N32cf/fTqq1VXU7NjjI/FpHvu8ftbS6i5I/qOxBVA35VSd9+wQ2HyimxPs7k9dfyNRKT166VHHvFdBzNwzmn48OFqbGzUpK4uf78VKxIXFCOxBwCgVIXNCnPUtLC1tVWnnnqqHnjgAb1vpp/fcotOOuus9J9PolHp+OP9/laqosoapcIA+q6Uuvumnu+W3LGwuXnHkmHnpK1bfVBLVzIUdiFM7UaYZOrUqaqurvbNH8IE2CxxQTYHqgMAUCkyNS1MLiFub/DxSN8AACAASURBVO95m0/82qC2Vpdccon+8pe/qKamRg0NDaqvr/fXpPt8Eon0GNNRPkhcAfRdHw4DL5pMyePUqVJbm9/LGovt+HfU1Unnn++/T9ov2/rMM7rt+ut1+te/rklHHKHpYav+CROkM87oHiRLKbEHAKDYMn1uCGP1Xnv5uPzUU9KsWemT1/i1N3Z0aPny5drVTEeOHq05X/hC92ZMnZ1+Upn9rAMOiSuAgamn5LGuTho50ge3INgxsIUJa2entGqVJOmGa67Rff/7v9ry7ru64YgjfGnwPff48qPUmdxySOwBACiEnpomhjG6ttZv0eno8PE1XJ1Nvj5+7eI5c/T222+rtrpaV+62myLLl/v+FJJ/zm9+k9iqM3Nmnv84FBJ7XAHkRqnt6+ypmdOqVX4v66pV6ffphjPAZtKMGQpqazVs3DjVjhmjwz7zGX9NWBqcXCIMAAC666kvRhir6+r8SuvMmT6uprk+kNTU2amJhxyi0aNH66tnn626WbOkadMSnz+iUenQQ/0Kbo7206J09LriamaDJX1T0ljn3Nlm9lFJBzrnHsz76ACUj3LZ15k8u9ve7r+m7oWNXxPU1irW3q7O1au1y5AhOv3rX9dxs2b5UuNnnpGOOcaXIgEFRmwGUBaCQOrqkiZPzn77zIQJUlVV4vr4im2ss1P33HOP3nvvPV100UWaFZYUNzd3//zR2JhY4cWAks2K662StkoKP521S/pB3kYEoDz11s23VCTP7k6f7pPX1Jnd+DVLnnhCP/zhD/Xqq69q5syZmjVrlt9H09QkPfqoX7WVSmulGZWC2Ayg9MVi0sqVvvQ33X7T5GqtcAI8bNIUi23/fetvfqPfLFigj3/845o1a1YiHkv+2gMP9CXGbW3+d+E+V2LzgJLNHtf9nHMnmdnJkuSce8eM2jgAKUppX2cQ+IDlnG/GJPW+vyZNwv3yyy9r3bp12rBhg84666zEC8mdh8tlpRkDDbEZQOnrrVlhcgyNRn3C2dXlk9B4j4lW53T6Y49p/ebNGlRXp6+feuqOMX3p0sSxdFOn+h4Ukl+5JTYPGNkkru+a2Z6SnCSZ2X7ys7wAUJpaWqT58/0el+pq/7tMyWWahLutrU1NTU2qr6/XkCFD1NDQ4Gdxm5p8sprcebimxn8t9ZVmDDTEZgA966kpUqGev2JF5l4QQeAT1SlTEmOsrvbxesqU7T0mLp8zR/8IAo0aNUqnn376jsnuwoXSunXSnntK48cnzm01IzYPMNkkrpdLekTSh81soaTDJZ2ez0EBQL845wPgIYf4oLVxo09ia2t7feuyZct03nnnacuWLWpsbNTcuXP9C/PnS4sX++/DpFUqrZVmVJK8xWYzu0XS5yS97pzjUx9QropdERSL9bzyGTZLnDFjhx4TikYVSLrpppv02muvacSIEfr3f/93TZo0KZHwhiuzHR1+lfWQQxJH4HBu64DUa+LqnHvMzJ6U9E+STNJs59yGvI8MQGXKxQxx8kHnkYi/35tv+n0zNTUZ7x8EgS655BL99a9/VU1NjcaNG5d4Mbk8GCiyPMfm2yTNl/TrHN0PQDHk60zxbON0byufqeNL2uYTbNyoebffrvvuu09DBw3SWZMn67h//MNXP9XVSZs2Sbff7u8/c2bxVpVRUNl0Ff6CpP9zzjXFf64xs39zzt2X99EBqDy5mCFOXQVNDo493H/JkiXauHGjhg8fri9+8Ys6+uijEy8mlwcDRZbP2OycW2Zm4/t7HwBFlq+KoGzjdG8rn6nji8Wku++WJC3585/129/+Vps3b1b95Mk6e8MGRZ57Thoxwh+b8/jj0j/+Ib3xBlVPFSSrUmHn3L3hD865jWZ2uSQSVwC530OTjxni5OCY4f6tra362c9+pnfffVef//zndcUVVyQ6FhZ7nxCwo6LGZjNrlNQoSWPHji3EIwGUinyt5Eaj0gknqPWFF/SzhQv1/vvva/LkyZr73e8q8vzz0iuvJJoi7rGH9JnPSCeemNsxoKRlk7imOzInm/cBqAS53kOT7z2jae4fBIGuvPJKvfzyy/roRz+qiy66KJG0SsXfJwTsqKix2Tm3QNICSZoyZYor1HMBlIB8xOn4BHFTV5car75aW7Zs0SGHHKLrrrtOdXV10qRJ3a/fb79Es0RUjGyC3Eoz+6mkn8t3L/y6pFV5HRWA8pGvmde+6MeKaBAEWrhggao2bdLUKVN02Te/qbpwL2yaZhFAiSA2AyhfqXF7xQotuvJKnfnEE9rinPbdd99E0pqqvT2rvhUYeLJJXL8u6T8kLZJvALFE0tfyOSgAZaQUuuru5Ipoa2urrrzySlVt2qRPmGnWOecoYrbjvUrhbwS6IzYDKE9B4I+w6ejwP0+frrZXX9U3/vxnvb1tmwYPGqQFCxb4DsLJ7wkT1Cz7VmDgyaar8GZJFxdgLACQ0JdV1J1YEW1qalJjY6PeffddHTl1qk9a6+u735O9rShR+YzNZvZbSUdKGmFm7fL7aW/Ox7MA9FMu4lS+Yl1Sl2BNnZq4dyzmk9YxY6RoVG2trZqzYIE0ZIj2lXTdNdeoYepUqbm5++kAyQlqct+K8GicICBWD3DZdBU+QNKFksYnX++c+3T+hgWg4vVlFnUnVkT/8z//U6+//rqGDx+uy777XUWGDfMzwA0NiXs1NzOTi5KUz9jsnDu5v/cAUCA7u+KYnKzmetUyvHdX1/YuwaquztgksenKK/XCiy9qWE2NLr/6ap100kk7xt/k96Qm2tXV/tp0Z8ViQMmmVPguSTdI+pWkbfkdDoCyke/VyDztKw2CQCtWrND++++vv//97zr//PN9OdL8+dLixf6i8Ngb9raidBGbAfSc0PUkOVnNdawL7z15snTCCT6B7exMrIgmTTbfPHeufvToozqsrk5fvOCCxDF00eiO78s0qUysrhjZJK7vO+d+mfeRACgv+d5XkhykcpQkB0Ggiy++WPfcc48OPvhgXXTRRZo1a5Z/MTxrrqGh+/OYvUVpIjYD6Dmh60lyspfrPg6p9w7HlbzqKmnZsmW64Prr9faWLdq9s1P3HHecj721tb7xkiStWrXD+3ZIVOlDUTGySVwfMLPzJN0raWv4S+fcm3kbFVAM5bifsZhjLuQMZy5KoSIRrVixQosXL1YQBHrjjTc0a9asxLE3dXX+UPNYTFq92gfLvj4PKBxiM4Du+hKX85nspd47zerpsmXLdPLJJ2vre+9p8ODBuuTb307E+r328l2DJ0+WZszY8e8hUa1Y2SSup8W/fivpd07SR3I/HKCIyrEzXTHHHIkk9sbkO3HOthQq9bX4v5/Wv/5V85YvV1dXlw444AC99dZbuvrqq7uf1Sol/n1OmZI+WAKlg9gMoLtSS+iSY3JV1fZV19aqKp122mlav369Ru61lxZ9+9s64rjjEu8LV1zDGFxuiwrIm2y6CnOyLypDOe6RKPaYC5U4Z1sKlTqeaFRBV5dm/+AH+sMTT6i6ulpf/vKXdfnllyuycaN09dXS+PFSuKems9MnrfX1BEiUNGIzgJKXYR/tjZdeqg0bNmjw4MG69Etf0hEbN/prw27BQZAoFS7HRQXkTTYrrjKzqKSJkvYIf+ec+3V/H25mt0j6nKTXnXNllC1gQCq1mcpsFHvMxUice3pmmn0vi9et05PPPKMPPvhAI0aM0OzZs/1K64IF0s03S+PGSaNH++tXrfIrrSStKAP5is0ABpBCbSlKfU4Q+PLgyZMTv4tG1bp4sf62erVqa2v1r0cfrVMPPVTasMGvsoby2TgKZS2b43Aulz/PbaKkhyQdK2m5pFwEx9skzc/RvQBkI5dBrC+Jc66e29MzU15ra2vTLQsWaBfnNG70aP3q/PNVV1PjXxw/3iethx/ePSASHFEG8hybAQwU+V6xTD76ZuVKX7lUXe2/Jk0Gt7W16bvnnaf7li7VLrvvrs8cc4wunjpVkdtv97F76lTfa0LKb+MolLVdsrjmBEmfkfSqc+4MSR+TNCgXD3fOLZNEIwmgkMIgFovt3PuDwJfrBkH/npt6n529bwZtbW2aO3euOtev18eHD9e9s2b5cqSWFv+cAw+UjjpKmjDBvyEMjqy2ojzkLTYDGECi0dz0bMgUo8PY7px/jpn00EPSs88mVlsl3X777Vr42GN6a+tW7bLrrrrssssUGTrUN2I69NDu4yMeI4NsSoXfcc59YGbvm1m1pNdF8wegfPW37GZnZ29Tn5t6n97uG87qJjdtyBDUli1bpgsuuEDvvfeeJh5yiH50yimq++Qn/fs6OxNdC597zv8zejQzuig3xGYAvcu0YtnXKqhMMTp1dTQIpJdekjo6/MqrpKarrtKtv/qVBu+5p/b44AP97IIL/Pnp48b5a2i8hCxlk7iuNLMaSTdJWiVpk6Qn8jqqJGbWKKlRksaOHVuoxwIDV3/LbnY28U3XHj/1a2enLzcKDxtPltomX0r7d7S1tem8c87RCy+8oP33318/uu461dXVJZo9TJzoA2VVlbR5sy8ZpjwY5aeosRlAGegpOe3rJHSm2B/G9nBFNhr1R8u1tEidnWpbvFizr75a67q6NPbDH9ZjV12luqOPzi5xLsdjCpFX2XQVPi/+7Q1m9oikaufc0/kdVrfnL5C0QJKmTJniCvVcAHmWmshGIj6hXLrUJ5WpgTQMlqlt8pMEra2aN2eO3nnjDdXsvrvOOeoov6f1lluk//5vafhw6ROfSJzXunWrX20lIKLMFDs2AygDmZLTIEh00Q9jaaYkMbnaKZ3w9XBPa/isqio13XCDzluyRBvfe0+j9t1X8+bPV11Dg3/PwoV+VVbKfLQeHYWRImPiamYHOeeeM7NPpHntE865J/M7NABZKfSMZD4DSTaruTU1iQYOSYIg0LxLLtGfn3xSn/rIR/T5f/s3HX3aaX68v/iF9MIL0sEHS/vsk/j31duzgBJDbAaQtUxxLhbbsYt+ptge/n7QIOn556Xjj5caGnZ8/cADfUVUPMFtq6rS7McfV/sbb2i33XfXGV/5ihrC98ViPmmtqfFVVi0tvrFT6rOJ00jR04rrN+RLdK9J85qT9On+PtzMfivfFXGEmbVLutw5d3N/7wtUlELPSGYbSHYmoU5XctRbUJVPWhcsWKD7/v53de62m/71s5/VSVdc4V+sqZH++Z/991/6kjRpEp0KUc7yHpsBDBCZ4lxyBVMYazPF9vDntWt94mrWPb6Hr3d2+tfXrFHbE09ozm23aesuu2jU6NGaedRROm3ChMQ2oNpaH5vfecc/f9q09A2kiNNIkTFxdc41mtkukr7tnPtDPh7unDs5H/cFKkohZyT7koz2J6FO997UQBsvGQ5qa3XTXXfp4Ycf1pCaGv3T9On64jnnJO4ViUjnnOOvb2hIu1oLlItCxGYAA1yYEDY3d4+16WJ18oTy6NGJst7U98W7DQe//70uuOMOPfb66xq577665pprdNLo0f76oUMTR+U8/7zfrlNf74/CYcsOstDjHtd4x8KfSKov0HgA9FUhZyT7koz2NaFOncFNbdQUP7x8+76YMWOkN99Uyx576LHHHpMkff7zn9fZZ5+tiNR9xba93Td0am8ncUXZIzYDyInkON3bxHTyZ4108T3ep2LJ8uVatnat3pU0dM89dfTRRyeu6erynyGmTPElx2Y+cSVpRZay6Sq8xMyOl7TYOUdzJKCS9SUZ7SmhThcgU5PidI2awn0xY8YomDZNseXL1bXnntp777116KGH+qQ1Euk+ixwmwclNKIDyR2wG0D/JcTp19TXb94WCQG1//at+s2WLqqurtZekn5x5po/J4T3b2vxRORMmMImMnZJN4voNSUMkvW9mWySZJOecq87ryACUnlycB5faTTDTDG66JDn+fVBbq8uvvVaPPPKIzjzzTJ177rmK1tYqEo6httY3iaiqSjxr5kxmdTGQEJsBeEEgrVjRvxXMdDE3ObZLmbsOr1ihRXfcoW888IAGDxumfzrySH9+evJqqyStXi09+aRPWklcsROyOQ6nqhADAVDG+lJCnLRqukOZUfje1GCZVPYbRKO66aabtHjxYq1fv14tLS26+OKLu88WS740ePny9M8CyhyxGcB2sZh0zz3++3THyaWTzYRzcmyXMnYdDm67TRc/8IDWbt2qfYYMSZyfnsos8ZUzWrETsllxlZlFJH1U0h7h75xzy/I1KABlpi8lxMnXRiI9lw13dvrSoqSz3v7riit0/aJFqh40SOMmTdI3v/nNzGOorfUzvGGwBAYQYjMAST7uhXtGs52kTZ1wDn/u6vK/c06aONFvs+ns9N8nd/6Nx+425/S9l17SB7vsouo99tA5557rz09PaqK4Pb5PmODPUp8wgTNasVN6TVzN7CuSZkuqldQq6Z8ktYiW+wBCfWkQlXptpg7CnZ3Ss89K69ZJ++0nRaNquvFGXXHjjdr83nsaO3q0Hrz8ckWOOCL9fcPv29t33CsLlDliM4DtIhHfmTcW6/3acLK4trZ7Ipp8rE24eltd7WNnGKOrk3YirFgh3XOPbtiwQfesWaPdd91VJ59wgubMmZOI63vt5aufwvs3NUkvvui/NjSkPwIH6EE2K66zJR0m6Y/OuRlmdpCk7+Z3WAAqRqYOwtXV/uf99lPQ0KDFixfrsmuv1VtbtyoyZIh++LWvKVKf1FQ1U9kRB5hjYCI2A0iIJ5I6/nifFGYSJpUzZnSfzE0+9sbMr7gmx82wI7AkTZ+uZc88ox8++qjatm6V3ntP+40YoYtmz/bNmJKPrwtXXMNtQs75r+3tTCajz7JJXLc457aYmcxskHPuOTM7MO8jA1Cacr0vJUxSU1dFkxLOe++9V5dddpnefucd7Tt8uBYcc4yOOPzwxPMzNXwK709wxMBDbAaQkLx/tCepk7mpMT0S8Q2eUldvJ0zwMToaVVtrq8775S/1/KuvavCee+oT48bpuiOOUF1YZtzTc5OTWaCPsklc282sRtJ9kh4zs0DS2vwOC0DJyvW+lCDwK64HHrjjquv06Vq2aJG+f+GFevvddzV2/Hjd8YtfaFLyPp4waX3xxe0lxUAFIDYDSKiv355YppWcoCbH7nQrtemaMsVXaIMg0Ldmz1ZbR4cG77GHjpoxQz/+t39T3dChiRieWirc2eknqMPkmI7C2EnZdBX+QvzbK8xsqaRhkh7J66gAlK5clN4mB9BYTFq1yge455/vtura1tamC77xDb2+caNGVVX5pDXc0xoEvvlDV5dfad1vP2nWLLoToiIQmwF001t1UaZJ53QrtRmOowuCQFddcYWWPPGEtmzbpikHHaSbvvpVRVau7B7DwyPppk3zMTqlzBjYWdk0Z5onaZFzboVz7vcFGBOAUpZt6W1q+VFqshoGsQzlQ62trZozZ47eGTRIH91rL/1sxgy/0hpqaZHuvls65hh/Rist9VFBiM1ACcvXUS/9uW+mSed0K7XhPtWWFr8ndepUBZLmzZunO/7rv7Rt61aNHTZMP/7e93yviTBZDWN4S4v/5+23pbPP9vfsaTUYyFI2pcJPSvq2mR0g6V75QLkyv8MCUPYytdqXdjwSJ0yE4+VDy5Yt02mnnaYgCHTwwQfrjttuU93TT/vAF57p+uqr/qibadOYwUUlIjYDpSpfR730575hMpou8e3q8iXDU6cmntPV5SeHJam6WgtWrNB1116rD733niZVVenn9fWaVF2diOFtbdJTT/kE1jlp40b/cyzmXydOIweyKRW+XdLtZraXpOMlXW1mY51zH8376ACUr9TZ3eSV1VjMB7ZLLpEOO0w67rjtQbStrU0XXHCB1q1bp8GDB+vMM89U3aZN0pNP+tnbLVv8fTZskN59138FKgyxGShh+epmn+6+fVmFzbSfNUxQpcTZ6dOmSYcfrmD5crW0tuq2225TV1eXxg4dqjvOPVd1RxzhJ5Pnz/f3amqSFi/295g1y3/ty7myQBayWXEN7S/pIEnjJa3Oy2gA5Fauy5X6cr/UkuLw5+ZmP2P8zDP++7/8RTrgAGn6dLW2tuprX/ua3nzzTdXV1em800/XcaNGJQ5EHzdOGj3aP7+2Vho8uOe2/8DAR2wGSk2+utmnu2/qKmxPcTrTftYTTvBxNvns9KlTpZtu0o2PPKLr7rpL++y3nw458EBdecQRqjv3XF8hNX9+IlkNY3FDg38usRl5kM0e16slHSfpRUmLJH3fObcx3wMDkAO5LlcK79fVldivki6BTT7gfPVqHyTr6xOlSl1d/vshQ/y4olEta2rSmY2NWt/VpciwYbripJN0UjQqrVwpTZ4snXFGYuY23HdDMyZUKGIzUEbytedV2nEVtqdENtN+1pkz/URya6tPWhsapFhMc//3f/XjV1+VJI3YulXNP/2pb8TU3u4T1+Rkta5OOv/83P5tQIpsVlzbJNU756jHA8pNX8uVeguu4X06O3tOiJNb4T/5pP9d2C04EvHfB4F01lmJ9voXXqi/r12rfUaM0H+ccoqO3rrVJ7wzZnRv8LRwoU9cBw3y7fXZN4PKRGwGykFP54znQure1XSJ7EMP+f2ms2al3+caBNLatVJNjU9C29vVeuutumbJku2Pufr007s3Ylq0SHr5ZemLX+R4GxRMNntcbyjEQADkQV/LlXpboU1uwtDW5oNXOsn7WevqdtznkhRY29raNG/ePG3ZdVcNHjxYX6iv11nnnpvoTpicQMdiPvh/7GPSwQezdwYVi9gMlIkwbo0Z03PM6s+qbHJymlqiG43633d0+Os6O9Pvc330Uf99e7uWdXbqtPvu04ckvSfpwkmT1HDOOYkkeeFCvy+2o8PH97lz+zZeYCf1ZY8rgIEudaa2rc03XAjLgELt7f5Q8bBcKFWY4AZB90PHU7SuWKE5l12mt7Zt0ycOO0yzpk/XiXvs4e+bLnFO7UYMAEApyzZu9WdrT3Jy2tTk43N4n0jEr7SGSXFLi38t3T5X59Ta2anTzjpL7W+9pSF77KFfNDTorG9/u/txdh0d0kc+Iu22m+89ARQIiStQqopxDlzqCm1yl8Bw70oQ+BnbKVN2nD1OvXemQByLqemGG/SVBx7Qpi1bdOgBB+g73/mO6mpqEu/PdE9KgwEA5SLbuNWfTsTJyWnKmejdYqjke0wcdJA0YUL3e1RV6ea//EXfuOQSvbd1q4bsuadOO/VUHTdjhrRkib9POL6uLv/P0KGJI3SAAsimOdM+ksZIcpLWOudey/uoAOSmsVK6JLUv901uvJA8ruXLfdlTfX3PY84QiFudU+Ojj+q1ri4N3XNPnXXOOaoLV25765gIgNgMDDR9nZhNboIYJqopZ6JL6h5DpW4lwQoni9eu1c3z5qmxpUUfSBq6xx665oc/1HH77afIwoXS669LRx2V+CxRVeUbJ86YQfUTCipj4mpmkyTdIGmYpPiOctWa2UZJ5znnnizA+IDKlYtz4NIlfdneNwh8YAs79yYHyTFjEvtlkgNt6r1TA3EQqG3JEp3zox9p4zvvaPiIEbrwwgv1hVNPTTwzeWY4fF7YoAmocMRmoMJkqpJKboKYXBqcKjUuh0ffdHb6suGVK7Xoqad0TjxpNTNdN3++zjrrLL9dqLran5k+blyio//EicRlFEVPK663STrHOfen5F+a2T9JulXSx/I4LqCy5apMOF2Smu2sbmrSG/48Y0b3/TLJUrsbpoy9dfFizfne9/TKW29pt91205mnnKKLwjKj1M6LUuJ5rLQCodtEbAYqR6aqo9pan7ROmyY9/7xPKsOGiJkkH33z0ENSTY2Wbdumr9x/v7bFL7nw/PN90iolVmUPPtivst59t/89Hf1RJD0lrkNSA6MkOef+aGZD8jgmALkqj+3PntDUpDe1wUSm+2Y467WtrU1zbrlFz3V1adx+++nYf/kXnTZhQvcSpo4Of+/OTmZ0gfSIzUAlyVQlFTZJXLvWJ61/+pP//vLL06/MhveIxXxsfucdNbW26kuPPaZN27ZpF0lfrq3VJXvt5Vda16yRXn3V36uhwSewkl9xJS6jSHpKXB82syZJv5b0j/jvPizpy5IeyffAgIqWizLh/kruDNzcnPnoG8lfk1pClHTWa1ttrc477zw9/9JLGvfRj+rGH/9Yk8z8PUeP7v53rl3rW/VLO7b1B0BsBipJponi5HPVzaS99/YJ54IFUmOjfy15u01tbaKqacwYLVu/Xqf+z/9o0za/1jp+2DD9dPRoRR58UHrnHemZZ6T33pM+/OHECQIzZxbojwbSy5i4OucuMLNjJX1evgGESWqX9HPn3EMFGh9QGVJLg0ule25Yvvviiz4p3XNP//vUcuBYbMcSovhROK3Oac6pp+r5NWv0vpmOPfZYn7SmKwOePt2XL0ndW/UDkERsBipCNtuFwq05K1b4xkldXdKvfiX9z/9I48dL69f7JHXmTJ+0Xn21tG6dNGmSFj39tM763e+0+f33JUn7DBumnx11lCLvviu99JJvxhSJ+PcNH97zxDVQQD12FXbOPSzp4QKNBahcvZUG52LP687cIxbzSWt7uzRihLT//omkNbV78DHHSC+/nAhwkYiWOaczzzxTb73+uj4ydKiOPeYYnXbaaYmSo3QryvX1iRJjADsgNgNlKts4nK68N/k9YZVTLOaT0f328xVKr7ziz3N95RV/zZgx/n0LF/oEd/fd1br77jr77ru1Of6oQw89VPfdfrvqurr8BPXixdJxx0mbNvlV17BRY097Z4EC2WVn3mRmC3I9EKCiRaM97+cMg1gslvkeYUlveNZaNvfo7T3RqDRqlP9+n338geNSoilEUpKq0aOlrVt9gJPU1NSk4447Th0dHRo2fLhu/P73dcVPf+qPvUleUU59fvgaLfaBPiE2AyUum1gudf9MkO49YZXTn//sY244wXz22b554rhx0qRJvkpq3jyfhB5yiK7qpLmldwAAIABJREFU6tKUu+9WV/w2IwYN0u3XX6+6SZP8s55+2r/HzE8gb9yYSH6BEtDTcTh7ZXpJEkXuQC71Vhqcbs9r6sxtT6u2QeDLiCZP7n6P3lZ6IxHfTfC553wwW7nSBzPJlyA1NSWOy0lKZpuamnTyySerq6tLkaoq3TJ/viZNndq/M2UBEJuBcpZt/4rkzwTJ7wnjflWV33v6979Le+zhV1zD2Oqc9MAD0m67JVZMIxH97LXXdEl8YlmSxgwerN9MmqRJf/yjdMghiQor53w8T66MYiIZJaKnUuH1kl6RD4YhF/9573wOCkCKdIltunJdKX1y29WV/rDw3oJoEPiZ1+OP97Ouy5cnAlpLi/TYY9Kzz0pz527vcNi0aJFO+eEP1dXVpd0/9CFd9KlP6Yjqan/93Xf7M+RmzsycTAPoCbEZKFe9HBmX8T3JFUrh2a0vvii99prf39rQkGiQuGmT7zb8zjvS0KHSAQdo7hNP6MdJSeu4kSN13xFHaNKGDb6EeO1a6fTTfQLc0ZEoDWZCGSWmp8T1JUmfcc79PfUFM/tHmusB7Iyd3b+amnT2lNxOnpy+FLm3ld4VK3yH32OO8UlrckA7+GD/2ksv+e9nzVLb2rWaPXeuNnV2qmrIEF124YVqnDDBdz3ctMkfYv7MM34fayyWPpkG0BNiM1DO+lNpFMbw2lpp82a/PUfyR9fcfbf/+dBDpdNO2/6Wxquu0k0vvrj957rqav3vhAmqW7lSmjLFlwPfd5+fkD7tNF9JRTMmlKieEtfrJEUk7RAcJf0oP8MBKtDOBrFsOg+nnr3am9QkOuzsm9roQfLH3hx7rCSpra5O1zU2aslTT+mDd99VXU2N5l10kRouvjhx0Hkk4ps7Pf20T4inTu0+RgDZyHtsNrPPSponaVdJv3LOXZWL+wJQbo67q6nxR94MGeJXXt9+Wxo71sfphx+WjjxSmj1bc7/6Vd3U2rr9bfuPHau7vvQl1S1f7pPczZv9ymy8u/D2s2FpxoQS1dNxOD/v4bXr8zMcoALl48zW5AQ0ObntbXU3FvNJ5lNP+b2rYYff2lofyGprE+9vb5e2blWw++667LLLdM9TT+mDXXbRwRMm6N7/+A/VffKTifNfx4zxq7WRiDRokE+IS+XIH6CM5Ds2m9mukn4u6Sj5Y3b+bGb3O+dW9/feAJQ59vUWn8Pj6To6ElVMGzZIw4b5LTubN0v77utj8x/+oJvffls/XrRo+9vHDhmiJ046SZGTT5Y+8Qnp8cf9C/vsIw0eLJ14Ys8d/4ES0ONxOJJkZsel+fVbkp5xzr2e+yEBFSYfCVymVdzeVnejUZ+0dnT4/TLhsTSRiJ99DffXJF173bx5WvzCC/pA0qhRo/Szq69W3aZN0o9+5JtGHHWU70ZcV+dXadvbCYpAP+UxNn9S0gvOuZfiz7lT/sxYElcgn3qLz7GYj81jxvjJ3xtukP72Nz/JfPDB0h/+4DsJH3SQFrW365xrrtn+1klDhujWoUMV+fWvpeefl667Tjr6aGnJEl9RdeKJiRVWJpRRwnpNXCWdJaleUvz/TTpS0h8lHWBm33PO/VeexgZgZ0Wjfka2q2t7R0FJOx5jkzrDG4n4IBiL+fcvXeq/Vlf7a6JRf8/OTrW98oouePBBPfjCC5KkyKBBuv/BBzWpq8sHxb/9TZowwb/33Xd9Q6a6OsqPgNzIV2weIyl5r2y7pE+lXmRmjZIaJWns2LE7+SigQHJxFnq+n9tb9VXq6+PHS21t0uuv+2Tzz3+W9tlHF/zhD7r++kTxxbfGjdOPBg2Stm3zZcCrV/t9rOefL61fLz3yiF9xPf/8/v2tQAFkk7h+IGmCc+41STKzfST9Uj6QLZNE4gr0Jt9BMwj8vlEzX94bifiEcelSv2oazqCm7l/paYZ34kR/j64uf01Xl7+Xc2p98EGdtmSJnn7ppe2X//KK/5+9c4+vqj7T/bMCSUggIQHCLRGIoggERYkXKKjRFhyj04p4qUy91EvnzHjrmWrltGNrP3W0OnOm0HbaorZoD6hFaL0EQbyAYFAU5ZJwEUJAkkAksEM25EIg6/zx8HOvrL3Wzk7IPc/388lnZ6+91m+tjf30yfv+3vd5f46JEydSSM84g4JYUwPs2sVn0i6rEK1JW2mz5XHMDjtg2/MBzAeA7OzssM+F6FR01Og153293ITdfxtMmxaar+7+e8FZnbVxI/U+O5slwxs3AhUVmPH443h727avL3nskUfwUGoq8Le/0YSpXz9qe24uT3C/CtHJiSZwHWWE8RRfATjHtu3DlmXVt9FzCdF2dETmtbVF02uG65Il/MwEql7ZW/cx9w6s81mzs7nW2LF8LSsDFi7EB7GxuOnvf0f5kSMA6N7yy5EjcXO/fnyukhL2zNx6K4egJyWFgmkhRGvRVtpcAuAMx/sMAGWnsZ4QHU9beEk0975efwd4HTNu/lOn0nQpN7dxpVIgADzxBPDxx8AFFzAY3boVT23bhreDwa9Pu3/6dDx+xhnAvn0cm2NZwDe/Cfz7v4fWy8zUTqvoUkQTuK6xLOtNAItPvZ8F4APLsvoCqGyzJxOireiIzGtri6bXDNcZM9ir4gxCq6rYq2oCR2fGNhBguZAZcWNmswaDDFptOzRKJzkZAPBKYSFu27QJx08tnxwTgzdGjcJl55xDZ0MTTFdVNd79FUK0Nm2lzZ8AONuyrEwApQBuAXDraT2pEB1Ne5oBeu2iAtElk4GQm//atcDnn9OA6aqrQuvl5wP19Qxa/+3f8MGSJfjOSy8hYIcKH6YnJWHewIHAb38bMkQcPpx/Jxi9lzaLLkg0geu/ApgJYCpYQvQCgCW2bdsActrw2YRoGzoi89pS0fTbHfaa4Tp8OE0XjCgtXMhAND6+cbmwwWn0YLLBr77Kz+68kwFwcTHLkDZswAdxcfj+9u1fB60WgEX/9m+4bNAgYP9+Di43z1NczEB29272zEoghWht2kSbbds+YVnWfQBWgAUVf7Jtu7AVnleIrkFzqrK8zvVLjnv9HeA+FggwaTxrFnX1wAEed5YbFxZyDM7s2Xhl40Z8d968RrX8TwJ4NDGRo+eM83BSEjBsGPtZjXOwTJhEF6TJwNW2bduyrLUAjoN9LutPCaMQXZOuNIalOQLoLkkqKqJY5eSEPnOKrHNn1Fw3axZF06xx+DCQmYm8kycx66c/Re3JkwCA3paF/5oxA7mXX87rx4/nXNbUVPbmlJZyndJSrtNV/r2F6CK0pTbbtr0MwLLWWEuILkdzqrLc5wYCrFqaNKnp5LhX0GvKhG+4AZg5EzjnnMbj6J59Fvj0UyA7G/M2b8aDc+Y0WjInLQ2PjhpFQ8QTJzjbtbSUOr1rF9DQAFx5pTwnRJclmnE4NwF4BsAqMKv7G8uyHrZt+9U2fjYhRHN2h90lSWaszfDh/plgt4HT2LEsHx437ut7znvnHfzvX/wCJ0/dppdl4ZNVqzDRstj3umIFA15zD2f/rEbfCNEmSJuFaCOao7vucwsKGFjm5DS9W+sVIJtE8tGj/NwZtC5ezDnrsbGYV1iIB198sdFytyUl4YU77wS2b/96zjoqK4G4OGDwYPbMDhzI0TeqghJdlGhKhX8C4CIzF86yrDQA7wCQOAoRidYwgWrJ7rC5b25u48AxEOAOa3Y2RXDNGr46d2Tz8oClS/n7ffdh3uef48Ff/OLrpeNjYvDzu+7CxI8+ogB++SWzy1u2ePfRavSNEG2FtFmItqA5uus+N1LQ6/6bwOvcyZP5umULnYLT05mArqwETrkFPxUTgzmuoPXJjAw8eu65QG0tA+ejR6nRyckMWq+5BrjnHgWsossTE805rmHmh6K8rkksy7rasqwdlmXtsizr0dZYU4hWw1jSBwItu76ggNnRhQubv4bfvaN5JpPF3bo1/PiGDdxdLSnhOcY2v7KS6553Hn+mTkXe88/jfz/44NeXp8bHY/+99+LRQAD4y19YsrRrF4W1spLrCyHaizbTZiF6BKer8V6YQNYrQDTabLTSeW5xMY2UKiup0ZWV1NbcXL6WlQGHDuH5gwcx5+23Gy17W2wsHk1IYKVUXR0T1ACD4DlzgEceUdAqug3R7LgutyxrBYCXTr2/Ga3Q+2JZVi8AvwPwLdB+/xPLsl63bXtr5CuFaCdO133YWa7b3D5Pv3t7HffL4paVhXplcnNDPa3BIEuCx4wBXn8dSEgAhg6lsdOYMcCAAXh+0SL84Jlnvi4PTgXw+7Q0pO7Zw9msiYnAd77DIHfsWJUEC9H+tIk2C9EpaI+xde0xYcDtKwF4a+XixUwIFxYy0Jw2jVr9xht0Fh42DNeuWIG8urpGl90DYH5MDEferFjBtfv25esvf6mqJ9HtiMac6WHLsm4A8A2wj2a+bdt/a4V7Xwxgl23buwHAsqyXAXwbgAJX0Tk4Xffh1FQ66hrRao17ex33G3Bu5rmZnhmAWd3SUmZ09+6lEcSkSUD//szUjhqF2//4R7xo1gMwZeRIvJmTg9SdOxmgVlTQeKm2luVHgMRRiHamDbVZiI6nPYLK1pww4AxQgdDv7u/h911GjmTAWVFBnU1KAhYt4kicmhpMrqnBR66g9WUwW4XYWP707093/5oa4LbbpMuiWxLNjits214CYEkr3zsdwD7H+xIAl7TyPYRoOa3hPtzSNfyui8ZN+P33Q6VCs2aFembM+JuUFAa1CQncOR08GDhyBMVxcfj2T3+KLbt2fb30BSNH4ne/+x1Sjx7lga1bGbTW1VEchRAdRhtpsxAdT3uMrWvNCQPOABVonEx2vnphSpVnzKALcFkZ36enAwMH4vadO/GRSUSf4jmcCloti32sycksFf7ud6nvqoAS3RTfwNWyrCAAL2t9C3TiTz7Ne1sex8LuZ1nWvQDuBYARI0ac5i2FaCbtUa50uvfwGnAeDIY7GzpLiF99FUhLo3lDRgaKP/sMV/z1r/iyuvrrZafHxuLlb38bqW++yV7YhAR+kJzMgPfCC1v4hYUQLaUdtFmIjqerjK0z+u02OgRCmu73Pcy1ZWXAc88B1dXAoEF83bsXgTFjcM327fjIocvxAN6OicFlDQ080L8/f4JBjq8LBrvGv5sQLcQ3cLVtO6mN710C4AzH+wwAZR7PMR/AfADIzs7W/FjRvrRHuZK5RzDI8qDmBrDuwHfaNJYE795NMfWjTx9gyBB8sGEDZi1ahIP19V9/dD2ApYMHc4f288/5XP3706kwORm46SZg+vSWf2chRItoB20WQgCNg1Lj4+DWZqPfOTmN/0aI5u8Fc21DA82YysuBL74A4uNRfOQIpuzciQOO07MB/BVAZkICNbl3b+64NjRQz4cP106r6PZEVSrcRnwC4GzLsjIBlAK4BcCtHfg8or1pj93M06U9ypXM2lVV0QXJ7n83r+C6pIRBZ14e+2zNecuWsVR41ixg7FjkvfIKvjNnDk44lr8HwPzUVPbcpKfTOXj8eODWW4E33+Q4nZkzO+9/MyGEEOF0Bc3tTBhtHTCAu5lAuCliUhI/j5Qk9iMjg9cmJADnnMMe140b8XxdHe45erRRWUUugDfNm+pq4Phxnl9dzQqob36Tpk767yq6OR0WuNq2fcKyrPsArADQC8CfbNsu7KjnER1AS3Yz21t426NcyblLWljIsqHi4qYzvEDIKTg7O3xea3w8sHo1he2ee/j5unV0Op4wAbfPmYMXX3ml0dKPJSTg8bi40D2HDwf69WPZcUYGMGwYBVbiKIQQXYv2qCDqThhNde64GtxBbUmJvxmS198tgQATy6Wl/G/xz/8M5OfjkXXr8Mzx440un5qQgDdjY0PeFbZNM6b0dAavl1wC/PCHMmMSPYKO3HGFbdvLIPv+nktLdjO7s/CWlADbt/Pn4MHwDK/Bbca0di0FbOxYugQXFjLw3bOHP1u3Mnh98EHuoq5cibzf/x4vvvlmo2XvAfB4374MVsvKKJL5+SwN3rmTZUx33aVSJCGE6Iq0RwUR0HV3dv2cgVNTw4NCZ1C7bRv1MhDguX5VUVVV1FOj3UVFwJEjwPr1wHe+gwdeeAG/cQWt1wNYevw44GjlQXw8TRcnTaLXxPTpXevfWYjToEMDV9HDacluZnsJr6E9BTgriyW8tk13QL/ZqG4zpk2bKIC/+hX7Wnv3ZuBZUQEcOgScPAksWACccQZw6BDuLSjAswcONFryfgDzEhI44mbnTpYhJSfTvGnqVO6y3nEHMHFi2/4bCCGEaBvay/CotRLM7R0A+zkDN+XwX1LS+LrCwpBT8LRpjY0RFy0Czj6bieahQ+khsXo1HnnxRfymuLjRLW4D8AJADe/ViwFrYiJw5ZXAjTdybuuKFUw2d7dEvhA+KHAVXYv2dhpszx3e1NTQXFQgurIfMyt24UKWBZsSIvPT0MBMbSAAzJ+PaVu3Yq0ZbQNgqGXhOdtGbkwMR9z06UORTEwMZXOvvx74/vfb4AsLIYTodrRWgrkp/W3twNbruZv6DoEAjRUnTaJR0quvUksnTw5da/5uycsDvvoK2LGDlVDZ2cAll+CWDRvwikOXAe60vuA8kJEBnH8+8I1vsPUHYH+tbasKSvQoFLgKEYn23uH1w0+gzfGpU7nDCjAbu24dBS0/n8djY3HB2rXYaCz0AYxLTMSbloXM6mpmc03w2rcvcNFFwN13qwRJCCFE82itBHNT+tvaiWX3c0frDGxGz2VlUXcti4Grs2w4I4PHp0xhafGAAQgUF+Oc555DxYkTjZa8H8A854G0NOB3vwuVGRtNdia6heghKHAVIhKdZZacl0AXFwNPP83A9KKL6P5bXs7d0gkTKKYXXogP/vY3/OPmzTjiWC6zb1+snTABqRs20HxpwAD23xw5wnKkCRNYfqSgVQghREfQlP52hsSy8xlSU4Hc3MaJ5vx84KWXWAF19CjbgG6/HXnr1+PGX/wCNY6lEgAs7tMHubW1DHJjYhgIx8Wx/Piee6TJosejwFWIroBTHI0o5ufTmCk+nmI4ejSD2RdfBM49FzjrLDzy17/imT17Gi11FoBPrrwSqatXs4y4vp7lwaNHc7f1+9+ntX5H7zILIYQQfnR0YtmvEsqZaD56lO8rK+k/AeCRHTvwjMvRPwnApthYZBoTJtvmT1ISR9ht2sR1OkMiXYgORIGrEJ2d4mL2xuTmUhzXrKEojhkD/NM/ATU13GkFOK5m2zagrAz3rlyJZ419/ikuBbCud2/gwAFmdAGWCffqxZ3Xe+4JzX0VQgghTpeu6jLcFF6VUGYcnRlR9/nnNG861YZz7cqVyDtypNEyOQDeS05mkOto50H//tT4b3yDAaySyUIocBXCk84ktHl5wNKl/P2++8LLo/LzgTff5POecQaQkoIZ+/bhbVffzMMAngaAEyeAvXuBSy9lCdLw4RRI2eoLIYRobbrrGDvnSJy8PCaDg0Hg9dfpHGzbwEcfcRxdQwOmlZdjrWuJSwG8BzBoNfTuzRF3558PXH21elmFcKDAVQgvohXa9ghwc3P5OnUqd1uzskLP9MorwHPPsVT4qquAykpce+BAWNCag1NBq+Grrxi83nknMH48DSMUsAohhGhtOkMvamvj1P6CAmDJEu6qJiVxDntpKbByJUuEExIw4fBhFLiW+HrcDcCd1thYICGBSeX//E8Gwd3p30yIVkCBq+hetFYgGa3QNjeTHMkd2DgBu4PIzEzutC5bRqv9WbOYgQ0EgFWrgC++YOlvRgbGvPMOvvDaaY2Pp6gCNHzo1Yu7rZs2Adu3062wO2XChRBCdA46uhcVOP2/DZwaPW4cd1iLinhs5Ehgxgy26WzdClxwAT977z2guhqD6upwyLXcYwAedx7o3RsYMYL9rPfco5npQvigwFV0fdyZz9YoSYpWaJubSfZ7voICBqUGY3sPsBTYlCABFM6NG4EnnuD7jAxsLC3FN5YvR7XrdnMBPABQDA8e5LXnn8+gdfZs4KyzNAdOCCFE9+Z0k8xOjS4sZNI3OZmvmzZRV40D8MCBwGef4anqaswxCWMHTwJ41H0wKYmafMklHKUjhPBEgavo+jgFqb1LkpqbSXb2xCxbxmA0KQkYO5Y7qcEgRTEQ4E7q8uXA8eMUwhkzaMh09Cjwhz8AH38MDBuGW774Aq9UVobd6k0AuTExdB0OBhn89u/P8uCJEztH/64QQgjR1rQ0yWw0OiODGm3b1OBPPgFqazk6bvVq9rUOHw6ceSbw3nu45cMP8YoraE0AUAgg03mwd28Gu1OnAg89JE0WogkUuIquhzsT6p6j1tElSZEwz7dmDbO35eXAkCEMJq+5hsc3bqQxwyefMDi94IKQYG7eDLz1FtCnDzBsGG7fvz8saE0AkA9gIsCy4JEj2c8aH0+jh3/6J4mjEEKItqUzmRy6/zZwPhsQ/pzmeFUVA9icHAamCxYAM2dSoz/+mAnm0lJWNQ0bBrz2GiZs3x7Wz5oM4IjrGCyL7sOVlRxB19H/RkJ0ARS4iq6Hu+SnswerXmRlhXZYnTb3WVkUSssCzjsvJGyTJ1Pc4uJorV9VhSsPHcL7Tut8ANcDWOo8YNuc0TpgANd74gmJoxBCiLbHXQ3lDg7bIrCNdk3nswGsgNq0ieW65nOjy8nJ3HF95BG6BJeWAj/5CT8rLOSu6fjxQHExMrZvR6nrVtMBrHDff+hQYNAgBsC1tRxrJ4RoEgWuouvR2RwKWyK+qanhFvdmHQD49FNmeG+9leJaUACUlXEXtqYGM48cCQtaGzkUGpKSuNN6/vnA00/T6EkIIYRoDSLpn1OrvXpM22JMTrRreo2VW7cOSEsLuQKbNUyV1KBB9IUYNIhJ5wcfpK7W1CCwZw/G5efjgOs2U+ERtPbuzfuecQZntNbUhKYHCCEiosBVdD062w6rEcqqqpCpUkuyx2ad7GwGrVlZ3GUdMIDZ3nfeAUpLMbmyEh+5Lg1zKATYzzppEtd58EEFrUIIIVqXgoLGu5VO7XNqtVfCuS2S0E4fCTM+zkuP3X9HTJgA7NjBtppAgDuh7me97jrgiiv4vqwMePddYP9+PLVvH+Z8+WXYLe4HMM95ID4eSEzkjNd//3dWVHWGMmohuhAKXEXXoTP1yzgx4hYMRs70BgIhh2DjGlhQQIEtKeHu6IABNGoyQWZBAW31587FUx9/jDmuftbeAN4FcJnzYGws17/hBmaHNaNVCCFEW5CVxaC1tJR65ZdU9ko4t0US2ukj4afHXv2tY8fSa8Losfk7w3luUhLXDQQ4Rq68HFfu2IH3bTvsMcKCVoCJ7ZkzgR//WIlkIVqIAlfRdWiLsiIn0QTGXucYoQwEGverej3/kiX8PSmJr8uWsUwoIYEZ3sOHKZpG1LKygIULkffGG5hTVtZouWQAb8AjaL3ySuCllxSsCiGEaFtSU7nT6gwEo6Utk9GRdnOdu8SZmaxmqqnhfNZ164Abb+TzFBcDc+cy2QxQt3ftoqtwRQWu3L4d77uWTgPwldfzpKQA3/0u8POfS5uFOA0UuIquQ1v3tkYTGPud4xZgr/dlZRwwPnYss7rbtrF0aNcu4KKL2ONiMr0OblmxAq+4gtazAOxyP9vQoXQQvvNOCaMQQoj2oaU7p22RjHZqr99Oa0YGn9n0tNo28OGHHGsTE8Pfb72VZcMbNoRabhYuBA4cAMaPR/Ly5Qi6bu2pywDv8aMfAffcI20W4jRR4Cq6Dm3d2xpNYOx3zrp1HG8zaxZNl9yCnJ8PvPgis65XXcUAdc0a9rACQHU1sH49M7nvvsv3q1cj4/PPUVpf3+hWZwHYZVkUW4BZ4OnTgcsuC5UGCyGEEO1BS3ZOAwG215igsLXuEann1uhyTg5dgLdvB/r1owdESkpofvquXcBzzwEXXsjz7riD65x3Hp6fPx93//GPYbe9GcDLXs+TmMidVgWtQrQKClyFMEQTGPvNggueyr2aYDIjI2SqBLDUKCWF7r5GpDdtog1+ejqwcyewahWwZw9w4gRQUYHRdXVhtvpfj7sxpUu2zfLgK64IF2khhBCirWnJzmlBQcg9PxrdiuYegQBNElNTvXtu3YlnY6YI0Cm4ro5lwuvX09xwxw7usO7YAfTvj1uuuw6vVFWF3dYzaE1LAy69FPj2t9nXKm0WolVQ4CrE6eB0Ar7zzpAIlpRQOPPyGFBOnty4/7WgAJg6lb9PncpypHnzgMOH8XwwiLs9bjUVjhmtZhRO376002/KGEMIIYRoC1rSxtPca6I5v6CAWpqd3Vhvnbu1znmyWVmsltqyhUnj+Hgmg48fZ9J57Vrg0CFg1SqM+f738UV1daPbxQE4AKBRSNqrFxPJzz4rAyYh2gAFrkKcjkGEU0yd13o5LZqgcs0aYOlSBrcZGRS3994DduzAI3V1eMbjNnMBPOA+OHIk8NhjzFh79MYKIYQQbU5z2ngi9aBGc49AwH/MjVuPzbnBIHd3Af5u2nqSkti3WlrKHdaBA4HXXwfi4oDCQu7eVlZi3B/+gC9cjzMQQIXXc950E/DEEwpahWgjFLiKnoc7UD0dg4jU1MYZXKfTsHFazMjgzqsZg2NMHgoLmd2dPx/4/HNcUFeHje7lARz2uu8FFwB/+hMwcSLfSySFEEJ0dk7XkMlcHwyGdlXdDv/ucydNCs1Gz8/nZ7ZNbY6Pp9dE//7cXd2xAxg1Chg2DCguxoBjxxBwPUI6gBL3c8XHAw88AMyZo7JgIdoQBa6i5+EWTq8SJK9dWL9jJmNr1jOkplIYf/UrXhcbC1x8MT+rqOD7VauA6mpk1NSE9bN6iuPAgcDo0cAPfxgKWoUQQoiuwOlOBzDXVVU1HQC771VQwJE3prf17beBzZuBPn0YqI4axYTy9u145L338IzLGBFgy84a98HkZOCXvwTuv79l30kIETUKXEXPwy1mXmVOXllh9zETtBYV0c3XHfgap+EVKzgjLiVLYIKkAAAgAElEQVQFOHgQKC+nUNbXY+OhQ5gEoMH1iJ62+r16AffdB3z5ZWgOrBBCCNFZcSd8T3c6gLNk2GmuFOlcgCXDxlHYHNuzh/2s6ek817KAPn1wQWFhWPUTADwJ4FH3wcsuA/7lX+jsL4RocxS4is6NX//p6fSlRiOcXllhr+xtaSmD1tmzeWzZMpYgWRaD1t27GXD27k2RHT+eYmnbeH7zZk8TphwA73k901VXAbffrn5WIYQQXYNoS4Obo+kt0X/j9J+UFOqRnTyZva+zZgGffQasWoV+hYU45nG57bXmyJHAggVq1RGiHVHgKjo3fqLXWoPL/QTQGdxGciQEQteuWcNgFQBmzADOPZeOvwBnwwEManftwpWbN+P92tqwx3kTQK7zQEICcNFFwODBnCWXmSmRFEII0TWItjS4OZreEv0vKQEOH6ZTcGkpe1337OG4m1Wr8MEbb+DyUnfDjk/1EwCMGQO8/LL0WIh2RoGr6Nz4id7p9skYIg0rd57z/vvsqSkubtzPOm0ajy1cyLE2M2ZQDAEGvImJPP7uu8AHHwAVFRhTVxfmUJgG4Cv3fWNigJtvpnOwdlmFEEJ0NaItDTb6lpHR2DXYJI4zMkI6eDrjd4xZ4kcfcW7rkCF46tVXMaci3CN4IoDP3Qctiy07P/yhglYhOgAFrqJz4yd6p9snA1AQy8r4WlvrPwfVCF4wyKA1Pb2xYP71r8y8HjsGTJkCbNzIsqPNm4GtW4G//52Bb0MDYgGccH8VeAStCQnAP/0THQq1yyqEEKIzcjptO06MppteVIDvTeJ4wADumJrjzdF/9zPOnk1NPXoU4265Bds8LvHsZ01NBX7/eyaUhRAdggJX0b1xC5bzfUEBjZPq6hpncb2oqgKOHqVYTp7MYyYrnJDAwNe2eV5qKrByJU2b0tKAffvws4YG/MJj2VywPLgR8fF0IpZDoRBCiM5Ma7XtGNy7qc6d0uZWHhUXc3c1LY0J5U2bWAG1YgUwaBAS7r4b7oadOAB1Xms98ADw859r1I0QHYwCV9G9cYqqmZ9qSn2zsoAbbgjNV/UTpIICYMkS/n7nnXx99tlQefHhwwx+33mHTsJ9+3JXtqoKOHECF+3fj089lg3rZ42NZRCcm8vdViGEEKIz01ptOwZ3NZXzfUqK/+6uV0nx4sXACy8AI0YAl1zCSqhf/xo4cABxx47BPewmCUCV1zPde6+CViE6CQpcRdvRWiVEp3M/p6gaF2BT6puayiCxKTIyaLQ0cmRogPnKlUB9PbB+PUfd9O/PsuPKSh6PjwcCAQyorg4bXg4AuwE0Kv7t1YtB8ZlnAjfdJIEUQgjR+WmNtp1o8RpJ56ygev99tuysWsVgc+RIVkJ9+SVwxRXAoUOYVlSEtR5LJwM44nXP228HnnpKmixEJ0GBq2g7WruEqCX3c4qq2wW4KYwoBoMhoyWAO7SpqTy+ahWD1qFDge3buftaXw8kJsKqrg5bMgFA2FHLAn7yE+ChhySOQogOw7KsGwH8HMBYABfbtu1VLCJEx+A1ks6YK+bmckbrn/4E7N3LVp6bbgKuuw54/XWgrAxjli8PM0YEgPsBzHMfHDQI+NGPGABLl4XoNChwFW1Ha5cQne79/Ebc+IlSfj5LhGfM4C5tURFLhEeOBCZNYlnwvn0MYGtrgf37gbo6vALglqNHw5Y7B8AO98HkZOCNNzjEXAghOpYCADMB/LGjH0R0ASJ5SJxusOe1llvDjadEURF7WWfPBh58kBVPtbXAI48wmXzgAFK2bw/bUe0LIFypwST19Ok0W1TQKkSnQoGraDvas4SouffzKzly9sdYFj9PSqKIPfssS4QHDwauvJJmTRUVPO/IEaCuDjMAvO1xu9sAvOA8EBfHdX/0IwWtQohOgW3b2wDAMv/fJ0Qk3DramlVWTa3lTCwDTCSfeSYwdiwD2mXLgGAQxb164UzbDrvccwRdXByNnHJzgW9/WyPohOiEKHAVHUt798EavEqO3n+fmdodO0JieO65FMLUVGD8eI64mTCBAeuuXexv7dsXOHECKfDukTkMjrxpxLnnAk8+GXIoFkIIIboSfg7ArRHwuddy/63gTCyPH89y4dWrgf/7f4GPPwaOHsVTAOY0NIQtnQ6gxH0wMZGj7ZKT2//vESFE1ChwFeG0ZzAZKavals/h3p014vjFFxS/bdto6ABwBzYzExg3Drj0UmZj9+0Dli4Fjh8HKith1bv9CSM4FPbtC9x1F3DNNa37nYQQogksy3oHwFCPj35i2/ZrzVjnXgD3AsCIESNa6elElyKSA3Bzceo9EK797r8Vxo4FLrwQGD6cxoijRwOvvQbs3g3U1mIQgEMetwmrfgIYBM+eHZ1ZoxCiQ1HgKsJpT1MlvwxtINB4dM3pDBuPNgCuqmLgWl8PDBzI0TSDBvF4IEAH4ddeY9BaVAQUFeHeQADPemR0hwA44D4YF8fv8d3vAjNnRv99hBCilbBt+5uttM58APMBIDs7O7wWU/Q8TifZ7Py7Awj/GyQjAxgwgC7BTz/NYHXTJuDDD4GYGODss4E+fYATJ2CdOOF5C8//kQ4eTD1+5JHmPa8QokNQ4CrCaU9TJb8MrXt0jZtIAtmcvhuzTlUVsGABS4GHDWMJcGUlDZcKCjjIfN8+/r5tG1BZieHHjmG/x1fydCjs14/OwT/4gUqQhBBCdD+M1gaDLOGNNoAtLmaP6pgxjfXe+XtJCY2Wli4F3nuP+nz8OI8PHgzU1uKV4mLcUlMTtnx/AJXug716ARdfDNx9N3D99dJlIboIClxFOO1tquRFU6NrIgWjzem7MetkZ3OXtXdvliDdeCNdCouKGKi+8AKzvcEgcOwY/KxLPPtZR40C/vY3YOLEiF9ZCCE6EsuyrgfwG9C7Js+yrI22bc/o4McSnRV3AtlobFVV86q28vKA5ctDlUjOkuE1a/i72XGdNIk+FLGx9KSoqADKynBBfj42eix9FoBd7oOxscDvf8/7KWAVokuhwFV0TpoKniMFo9H23Rg7/TFjWH503XXAoUMs5c3MZM9Lfj6weDHPDQQAwDNoTQRwzOs5hwxR0CqE6BLYtv03AH/r6OcQXQR3AtlobSAQMjmKhtxcoLqajr7r1gGfOsYHm/WDQX5WWMgKppMnWSK8Zw/679/v6SfxMICn3QeNCZP6WYXoknRI4Koh5z2I1jJYcq8T7a6wV79rfj7NGGwb2LCBWdwdO1gWvGsXy5Asi6VOb77JXVjLwgO2jd943MLThKlPH+Bb3wJ+8QsFrUIIIToPraXLfgnk5s5Mz8ykw/777zORPGAA9besjO/Lylj5VFZGjT7Vy4pDh9Br/36Eu0z49LOq+kmILk9H7bhqyHlPobWMnlq6jle/60svsVfmhhtYapSQwPKjAwcoioWFnNdaW0tn4UOHMNq2UeSx/M0AXnYfHDMGePlliaMQQojOR2vpcqQEsglYg8HQDqrzXL8y42CQieS336b5UlwcsGcPj1dWskqqrg5PVVVhjsdtYwCc9Hqeiy4CXnmFQbIQosvSIYGrhpz3IFrL6Ckri8JlHH6jzRK775+RwUxuXR17Y3bs4M+ddwJTptCU6aWXKKi9egFHjsDyGF4O+GR0R48G3npL4iiEEKJz0h4GjCY4njQJyMkJv5fbyCkjg8eHD+eOa2IiA9Xqar7u3QscOQLExGDQiROeo2483fwB4JZbgP/5H/WzCtENUI+riI6WlhZFU9IbzdqpqRS3999n70y0WWJnz82yZcCWLSwzuvxy4LzzgO3bue4XXzAozszkTmtiIk2YPEbd9ALgabafkcF+WAWtQgghOittZcDo1PKsLGqqZTXWdnNOUhID1KIiJnsnTAAOHgTOPZcafOaZwL/8CyuhFi+mVts2Yk+c8NTfJwE86vVMDzwAzJ3b+t9VCNEhtFngqiHn3QxnaVFWVuv0x3it3VJDpqaC34IC4NVXudN63nk0gVi0iLNZa2uBFSuAM86giJaX48qaGrwfvgomAvjc69nuvx/44Q8VtAohhOgatFavq8Gt5cnJfJ+UFNJ2c86AAWzZ2b+fPykp3FG94AKaJP7lL2y32bQJWLsWqK31dfP3rH4aNw547DFg+vTT/15CiE5DmwWuGnLehYhGvJxBY2v1x3itHYlI7sALF3Lua1VVyM3QmeHduZO/X3kl57T++c80e+jdm7PgqqoolocPo69to9rj9p7zWQHg4Yc5EF0IIYToKrSWlpu/ITIyGpcFe2m7s22npIRBbWYm9be+Hti6lT2xu3bRa6KuDj87dgy/8Lit53xWgC07b76pRLIQ3RCVCovoxMsZNLZlf0xLMsAFBQxa09NZluT+LgUFwJIlzOpeeCHntAYCFM1TJcGoqQFsu3n9rGeeCTz4IPC97zX7awohhBDthpe2toaWOxPH11zT9Cg6Y8TkfJZgkK08/foBmzdTl+PjgZoa9Dl6FHUet80CsMXreS68EHj+eQWtQnRTOmocjoacdyaaK16t3R/jDJyB5mWAzSzWadNop19ZCezezWyus5dm9GigoYFlwiUlDFTr65nlPX4c9wJ41ucWnkHrzJnAc8/J7EEIIUTnxhlcAiFtbQ0tdyaOo/0bwt16FAyyLPi119iuU1cHHD0K64Snm4S3JvftC8ybB1x/vXRZiG5MR7kKa8h5Z6KtjBqixZg4BIPcDfVyIPSjoICzWHNy+D3WrQM++4y7ocbM6eRJYPlyBquLFgG33sprT54Ejh9HLLzNlqYDWOE+2L8/Z7N+73sSRyGEEJ2flgSX0eJMfPtpotfom7Iy4N13gXfeYbtOaiqD1hMnAMtqXtB6xhnA669rBJ0QPQCVCou2p6ny39RUbxOHaNZx9susWcNy4AMH+GrbnKn62mu8NjaWwnjbbbTWB2e+eQnhy+CM1kb07w+sWiVxFEII0XE0t6UmUnB5OgZN0V67bh3NEWfNYjkxQC196y3OUT/nHOCqq4DjxzHi0CHs81iiH4Cg19qXXsqEtEqDhegRKHAVbU80PbTROAaXldH99xvf4O5pbi6dCAEaOmzYwHJgAPjyS5o7VFZy1A0ADB3KMuK9e/EKgFt8Htczozt4MPDTnypoFUII0bE011QpUlVVSwyajCZXVVF33de6A1rjHREMMsEcDHKOumUB5eV837cvBmzciIDH7cYC2Oo+GBMD/Pd/q/pJiB6GAlfRtpge1Ozs0K6oV3bWS1id4rh2LYPQujr+vmkTB5NnZdEhODaWI25KS3l80CBmeVevpsV+QwOD2/p6ZAAo9XjUBMDTTRg33ADceScwZUqr/JMIIYQQLaY5vhRN7Yq2xKDJBLvZ2d6tPe5geNw4oLgYOHqUen3BBdTsxESOxDl8GNZLL3neyneUxH/8B2e0CiF6FApcReviFklnD2pJSfMyuwUFdBpMTeXOam0tMH48S4u++orB6ZYtwJ49HFx+8iT7Y4YOZZ/r6tXcpXXgNwfuegBL3QdjYoB//Vfg8ceV0RVCCNE+RNNeE+3uaFM7qtGs5dWjCoTeBwJAXh53UCdPbuxbYRz8Dx8G+vRh8nnZMv5dcOwYnqquxhyPW8YAOOn1LP36AXPn0oRJCNHjUOAqWhe3SEaa49YUWVncWS0tDa1lWcDw4Swz+vhjYP164JJLGKB+/DH7UMeN4+7qwYONlmvW8PK4OGD7dvXNCCGEaF9ac1Z6a4y8cfeommA3EGAVVVUVR84BIZ8K41sB0HeirIzGiIcOAR9+CFRXY3R1NYo8bjcEwAGv55gxA/j976XLQvRgFLiK6IjWhMFY21dV8Rp3Ntcpdk25EK5bx/LftDT2yFgWr6upYVD6ySfceU1IYHBbU8Ny5IoKYOdOliI1NGDeyZN40OdxPYPW1FSJoxBCiI6hNWelOzW4pUZMtk3n3y1buKNqrnWWDN9wAzU6K4v3CQaBSZP4On8+A9YtW1g5tW8fLONH4eIxAI97ffDII8Cjj6r6SYgejgJXER3RZoBTU0NjaJKTvc+NZq2CAmZ4AQ4U37iRfarbt7NUKDkZ+O53aaVfXMwAtndvBqxbtnBGa3w8Bp08iUMeyw8DUOZxHNOnAy+/LHEUQgjRMTR3RJ07IPULUJu7k2vWGTeO7r2lpXxvrs3IAOLjGZxOnx66RzAIfPopXf2XLw8lnjdvBurqfINWz0Rynz7Ab38L3HVX1P8cQojuiwJXEU6ksTPRZIAzMoABA/jaks/NfWbNouCNG8f+mddfp3CmpVEsN26kuUNsLEWxro4/p7Acvzt5EsCjXh/cfDN3WhW0CiGE6Cq4A1K/ANXde9qU1pl1cnKA2bNDfxeYvxGCQWDHDv4MH85r3n+fO605OUB+PiunamuB7GxM3rkTH/ncyjNoveUW4H/+R5oshPgaBa4iHC/Ra04G2BgxlJQ0Lrd1ip35PCWFxzIy+N4Ey6mp7KUxJcMJCXQiHDSIZcJr1/K1oYEZ2bg4GjOdoln9rGPGsJTpssui+35CCCFER+FOLruTwX6J5qZmpkda1/k3wJo1oQB1xgzORbdt4L33+PmBA8Bjj9Hg8PBhoL4e8UuW4LjHV8kG8InXd3zzTY68E0IIBwpcRThO0WtJT4yfaJqA2GRjMzKAhQu5ixoXxwB18mTgoYdC91q3DnjiCc5lzcqiY3BREY0e6uspjDExwIgRQHk5RgcCnmYPvQHUez3r/fcD8+ZF972EEEKIjsadXHYni93mSdFWT/mtu20bf2ybY+GcaxQUcMd10SJg5Urep7yc157CL5H8HADPAuC5cxW0CiE8UeAqwvHKrALR77j6iabbQn/NGgat6ek0VCoupgFEWhrLkgCOtNmxgyXAhw/TgKm+nqNvGhr4EwwCtg3ryBHPxxkFoNjrg2uv5agbIYQQoqvgDj6bShYD/tVTzuS033pVVQxMKyv5Pjc3pPGff87RNqNGUZsB4NgxAMAjAJ7x+Qqe1U/Dh7MtaOJE/+8uhOjRKHAVkTkdd0O3aLrdDauq+H7yZJYalZezX/WDD4DqamDgQBol1dQAZ5zB63fv5vUxMaH7NDT4Bq2e4piQwJ1WORQKIYToariDT79Wnkj6bQLWqioGngDX8Fo3EAAKCzme7uhRzmENBjkvfdEi9rAeOEBd7tULGDIEvQMBzzmsSQCqvL7TvfcCTz0lTRZCRESBqwjhVRbcXHdDJ+7ROJWVzKbm5rKMaMMGlgwD7FlNSeFPUREFMjkZOHKEu7ANDTxWUUFhPBW4PgDgN7ZneBoetMbFAc88A5x1FsudJJBCCCG6Kk218phKp0gOw2PGNG2WmJoK3HNPyKNi4ULq9LFjTCynpQGXX85JAHv3Ira83DNoHQWf6qeHHwbmzJEmCyGaRIGrCNGaQ8+B8NE4mzYBS5fys9xcimVSEkWwsJBZ2+pqCuyFF3LszapVwL59DF6rq1kiHBMD9OuH3rW1nuJ4KYB1Xs+zdKn6ZoQQQnRN3IFqU5odCIR8JNznOMuAd+xobKZo7mNME52v773Htp0BA4DBg5lcTk4G3ngDOHQIVpXnfqp39VNiIvDWWzJGFEJEjQJXEaI1hp47hRUIDSFPSuJu6WWXMXhcv57lRvv2sfw3EAD27KEz8Nlns0x49246B9fX89r6es5qtSxYh7ymswJzwV3YRvTqBbz2moJWIYQQXRd3oOqsaioubhxkZmWFxtGcdx41+Le/pQ66DZySk0OaHQgAzz7LRPP55zOh3KcP38fFsae1shL4h38AhgzhzmtxMRAIwHI4+zvx7Wddu7bx5AEhhGgCBa49jUilRU2VBbuDUq/fCwoYkG7aBJx5JoeQZ2cDf/gDzZjM+qtWMTAtK2PGtq6OZUcZGexBffFF7rgOGcLjffrwuuPHfeezeopj797A3/+uoFUIIUTXxLkDmpMT0l1nVVNxcWgn9PBhBrSFhTxvwgQGiabiyTmT1an7xcV09P3iC+DQIV53+DDbaw4eBPbvZ5Dc0ABs3QqMHAlMm4aUP/4RXi4TvqNucnKA559X0CqEaDYKXHsa0ZQD+wW3zmsB79+zspjhXbeOfS85ORS6/ftp6rB/P/tc+/QB+ven+NXWMli1LIpkIMAd1pQUinJtLdDQgPjjxz3nwCUDnqKJ9HTOgpNDoRBCiK6Is9z3mmvCddsEse6y3rIyYPNm7rZOngyMHcvzcnP9/w7Iy6P3xIABwOjRTECvX08NHTSIJoo1NQyK9+8HiopgmTJkF/cDCBs0FxsLvPOOSoOFEC1GgWtPo6ly4Gh6Ypy9McEgd0YPH2aQmZoKjB8PbN/O49OmMYs7bBivmTIlVCq8ZQszt5bFXpeTJxnc1tWFyoNtG4iLg3XcK2QF7gEw3+uDBx4Afv5zmT0IIYTouhQUhMbGuXXbnWQ2O5iZmax8iovjrmlqKn/uu4+fp6SE9DsQCOmkqUw67zwGvZs3syw4P58+E7GxIcNEwDdo9ax+Sk+nA7GCViHEaaDAtafhHknj3ll1i6T7HGcgm5zMrO3hwyE7/TvuYHnSjBkMUgEGqTt2cM7bjBnMBu/dSwEcPZoiWlNDkXT2yNg2UFHhO7zcUxxHjwYWL9YuqxBCiK6Pe/65k0gVVGPH0uRw7FjvNh/LYitPUhKPv/02fSZuuokaXVTEdp7+/UMz1E+17AwC4OUyEQt4VkXhttuAX/9aiWQhxGmjwLU70pRFvsFL9NwiuWYNzwkGQwJn1jTnJiXx1bKABQsYgF59NcuFbZvlwKWl/H3BAmZ6i4q429q3b2h3t6EBiI9nX2pNDWbU1+Ntn0f3DFonTqTjocRRCCFEZyNabXYSyXsiUgVVSQkDzpIS/rz/PrX3k0+ADz8ErrsuNAanoAB44QWW//bty53X3buBbduoyWPHUsdtG72PH/d08/fsZ7UsGj3ddVd031UIIZpAgWt3JNqxNhkZ4fPb3CLpLAtetIjluz/5CYNE57k/+xnvm5REEwjbBh5/HBg6FLjgAuDWW3nekSPckd28mYJYX8/S4F69+JOaChw9isT6etR4PLKv2cP99/N+ClqFEEJ0Rk7HY8ILr6DWXJ+UFNL3lBR+FgwyaN2zh+7Ahw8DK1dSn4cMoaN/WhoTwEVFrICKiWHFVHV181yDAZowqjRYCNGKKHDtjrizsH5C6MzIOt393GVFADBuHHdLP/6Yr7/+tfe9+/enhf6773Lual0dsHMnhXHcOJo/HD5MIayv/9p4CZbFwLWiwtc1+DYAL7gP9uoF/PKXwKOPRv/vI4QQQrQ30YycW7cOePVVYNYsmjG5aSqwNc7+NTXU1bw87qAGg/y57jrupI4dC/z+93QPNs7+DQ0sH+7dG/jqKwauwSAm19TgI5/H9QxaBw0ClixR0CqEaHUUuHZH3FlYvyyvl4g6Z7jNnh2y2s/J4U7rggXsY3XiNHRKT+drfDzFsqgI2LWLIlZeHiofHjKEfa5ml/XIEaCmpvn9rH/6U+RdZSGEEKIz0NTIOYD66Hx1Y/TctO+457b+x3/w98RErlFayuB13Tpq9RVXcAf20CFg+nQGqOecQ7fffftokGhZNHY6eRIJNTWo9XiM/gAqvZ7v2ms5zk7VT0KINkCBa0/AL8vrJaIFBQxaKyspjLbNOazm2htu4Ow2IJT5DQYpjqmpHGtTeUrOyspYjnTsGJ2Fv/ySolhaynPq6xm4DhmC5EOHEPR5fE/5vv129s1EylwLIYQQHUk0pb/Oc6ZMofGhn7aZ42Vl3Jk991xWLgHAM8+EjBKfeor6vXQpXYKrq6ntI0fy+LvvskIqPT1UPnzWWUw2n9ptbVYiGQCee079rEKINkWBa0/AL8sbCIR2QKdM4XlZWdxptW3gwAEaNtx+e2OjJoDrmczvpEl8v2ULe2OOH2cwW1zM4DQujj2ygQBdCQ8eZAlSQwPQ0ABr61bPxz4HwA6vD2bOBP77v5XRFUII0bmJpq/VfU6kXVmj53l5fD9yJDB8OLX7f/0vJoUfeIDn/Pa3oUD2wQc5z9U5a92ZTD50iEnmmBgFrUKITosC155MQQEztgAzvEYsjXvw/Pl0Gdyzh8fdO7dZWaFdWYCCOXo0e1oLC9nLOnIkzR62bWPQmpXF+a6FhUBZGSyfcihfcXzySeAHP1DQKoQQovMTTV9rpHP8dmzNzqwpFQYYwP7DP1Czly3jTuvKlUwib93K3dXHHwe+9a3QTNa6Ol5XUgIcPYoRJ09in89jeury0KEsUXb6ZAghRBuhwLWn4DXHLSmJZUYjRzY+bjK/N91Ea/ypU7nbmpXVOBOcmsqgddEiimefPhTBbdsY0MbHA2ecQUG0bfa07t4NnH8+ri0tRZ7Po3qKY3w8d39vvrmV/kGEEEKINiaavtZI5zS1Y7t1a2hXNSuL5cCrV/NnwgTqb58+HIMzfz7w2Wf8PSGBpcJDhtCQqa4O8SdPes5hTQPwldezXXcdMHeuglYhRLuhwLWnkJ9Pg6QZM1iqa4yUamuZbU1NZXBbVRXqaU1NBe67j5lbt8uhKTNev55r7djB8+vr+TNuHNfdvp0/lsXssGXBWrHC8xGHASjz+iAnB3j+eYmjEEKInoWfiaIxRJw2jRppNHv2bJb8btnCXtXjx5lY3rmT1w4ezGC2Vy8Gr8nJwIABsPbv97y9b/XT3LnA976n6ichRLuiwLWnYJ3qWNm7l6KXnk7XX+NGCDCzu2EDRdApRl4uh6bM+Phx7tqaweWmX+bkSa5tTCNOnmTfTNDbgslXHB97jKVNQgghRFekObNZ3RjvCef1BQWh5PPkyeFrjhpFzR84kGXCt95Kbf70U+r47t3AmWcCdXUY/d57KGpo8Ly1py5bFvDSS6p+EkJ0CApcewqTJ4db5wOh3hjAf/5rejpw4YXcRS0upinE1KncgTWB6Gef0U7/2DH2zTQ0cJbb6NFAWRmSampwtDlB64gRHF6uXVYhhBCdCa/WG9PngPgAACAASURBVBNUepkeRir3jSaodV9v7puRwWoqywoFsAUFwIoVDFS3bqX78KJFTDC//TZ7UkePBhITMeKzzzz7Wc8CsMvrOcaMAV5+GZg4sXn/XkII0UoocO1MGAFzBpetVYbj7KExwaDbJdicEwjws2CQGdoBA2i0VFLC/pmlS2mtP3kyr/30U/a17twJnH02d1ePHOF6/frBqqnxfSzPoDU7G/jrXxW0CiGE6Hw4A0kg3G3fbXoYyXwpGtdh9/VGq9esYQtQXR13UXNz2e5zwQUMMINBlgSXlzOpXFXFcydMgLV8ueetVP0khOjMKHDtTBgBM4EiEBKy0yk18sMphs71nWNuJk3iWJyTJ5m5nTqV16SlNT5nyRKWBffrB1x0EfD3vwP19bCWLPG8ta843nwzBVcIIYTojHgFos5js2Zxx9UdaEa7lhu/67OyOFu9sJDzV+fOZd9qTQ3beIYM4Q7r7t0MWk+cAOrqPIPW3gDq/e7/2GPAQw/5P58QQrQTClw7E87yH2c5LxCelXWXKjlLkyorWc6bmxt519Iphs7dV+dz/PWvLAEGgH37aLg0ezbLk8aMYTnS3r3sb01M5OfbtgEVFbDKy8NumQjgmNezyDVYCCFEV8AdSLrd9o2JIdB00rkp12FTemxabZKSQtVOpkx4+nRgwQLupo4ezQSzcfgvL2cFVFUVBlVV4ZDPbTyD1v79aQKVm+v/fEII0Y4ocO2MpKSEB5zurKy7VMlZmmTKeQG6Ajvx6r/xWh+gY/BbbwEVFcDYsXQjLCvjz+uv815VVTw3NpZ9rZ98ggFFRQj4fDXPoDU7m703cicUQgjRnYimFLip6199lQGobXMHNSmJn5mKpgsv5E7rsGGshiotBa64gknloiKgvh4pVVU44rF8OoASj+O49lpg3jy17AghOhUKXDsTRuCqqhgUOjO07qysO9B0liZlZPCYc/6qWcer/8a5vrHZLypiuRHAGXAHD3LHdft23qeigsHqFVcwqH32WQCAVVTk+dXOAbDD64Np07jTqqBVCCFEZ8btQxGNH4VfKXC07T+m9NhpblhVRbPEG27gjuvw4TRkOnyYyebYWODuu/laWQnr6FHPpX1bdl5+WdVPQohOiQLXzoQRtmCweRlad2mSmb/qNl8y93AGuW7xXLeOP8nJzOBecw0t9deuZcb34EFeU1NDl8KxY0Olwc1xDU5OBn7wA2DOHAWtQgjRCbAs6xkA1wE4DqAIwJ22bVd27FN1Itw+FF5+FG68SoEDASZ7N21i641Tv72ud37u1PXkZO6+/uEPwCefcKc1NpaO/l99hZTXX/fcZQV8dLlvX+CVV1QaLITotChw7Uw4dz2TkiKbNbTEidDcY/JkXl9ZCSxeHBLPyZM5tNyYOowdy2v27GH5cnw857WWl/MnPh6YOxd9161Dtc9j+roGf/KJ/3cTQgjREawEMMe27ROWZf0KwBwAP+7gZ+o8uH0oMjI4ciYYpG5Hm4QtKKDuVlY2no/uhRlBl5tLHQ4G6S9RWEg93r6dzv61tTRRTEkBevVC4gMPwMvPPwlAldd9Zs4E/vM/VRoshOjUKHDtjLgztE7hysykQAaDdPONFNwaKisb97U6s8abNnG9/ftpdV9Xx2xtYSGDU9PDevbZfF9RAfTuTXfCvXuRsmOHZ9AaB6DO61meew64667m/5sIIYRoU2zbftvx9iMAszrqWTolXmPlSkqop0lJ3lMAgPDfMzKAf/xH9qCOG9f4Hu4qqLw8elZUVNDh/9Ah4MwzgY8/ZlVUIMCAdfhwHh8xAsnz53sGrb4h8sMPA08/fRr/MEII0T4ocO1M+PW8GOECWAJcUMAMa05O5AxvQQGvKylhH0y/fiwtyshg0Dp1Ko0c9u4FNm8G3ngDuPxy4Oqr2XcKUAjLy0P3BOgk3KsXrNpaz9t6imNcHJ9FJUhCCNEV+D6AVzr6ITo9Rk+Nt4TxiSgtDZ3jnvk6aVKo7aakhLukJrlsWSGtnTaNOr1rF8996y1WOsXGcv3Bg1kZdckl/FvgootgXXut52P6tuzMn69+ViFEl0GBa2fCr/zXBHvmNdLcN2fwm5VF18GSEuD884HrruOxdeuAzz5jUJqUxFKjigrurvbvT2v9pCSK6CefsKfGtimYp6z1LZ+v4CmOQ4YAy5cDEye24B9FCCFEa2FZ1jsAhnp89BPbtl87dc5PAJwAsDDCOvcCuBcARowY0QZP2omIZKRUUsI+15IS7sIWFDCoTE/3nvMKsGKqtJRrVVVRk41p4owZdOh/8006+AMMcm2b56ekAN/8JnDOOayM+uILugl/+imsxx8Pe/RkwLvP9coreU95TAghuhAKXDsTfgFpZmbjsTaR5r65g98HHwwvM/74Y4rs/v0MWI8f526sEUVnH+yAAUCvXgxkq6ow6MgRzzlwvq7BQ4eyB0jiKIQQHY5t29+M9LllWbcDuBbAVbbt34Bp2/Z8APMBIDs7u4lGzS6O0dVgMOQ/4R4ll5HBJG9GBs2U3OeYwNfpYxEMcnd10qSQaeLRowxaS0uZYJ4yhTuuQ4cCt94KjBzJNRMSqOHbtsHKz/d87ET4BK1JSWzbkS4LIboYHRK4yrnwFO4sbqSA1N3n6oc7+M3MpPFSfj7df4NB/g6wV6ayErj0UuC887jbatvAnXeyBCkmhs908cVA796wTPmwC9+/WKZPp9uhxFEIITo9lmVdDZoxXW7btp/nXs/D6GlZGXcpZ8xgT6nR7mnTQm6/OTl8X1zMkuHc3FAfLMDPvIwYjU7m5bG6qU8fJpIDAQap554beo4//5k7rSkpsALeU9N9dTkjg7osEyYhRBeko3Zc5VwIhO+ORipHcva5zp7tf54RxOJi2u0PHMgAddMmfp6cDCQmsrRo+nQGsqZ8eO9eYOVK4KuvWOY0ahQD2IaG5gWtsbF0J3zggZb+ywghhGh/fgsgHsBKy7IA4CPbtv+5Yx+pE2B0ddkyvt+7F9hxqsbIJJvdva5uzQbCq6m8ktVTptC46d13mViurKRmDxgQ2p3t1w8oLYW1ebPn43rqckwM+1lljiiE6MJ0SOAq58JTuHdHI424cfa5RjMKJy+PQ8Tj4oCzzuKOqmWxtCgzkxncYDAUMK9fz/eJidx13b8fqKxE/IEDOO6xvK9r8D33UByFEEJ0KWzbHt3Rz9CpmTyZO6RJSZxtboJUoHGva0oKjQ9nzKBmp6aGEsTG3d+vEik1Fbj3XiaOt22jVhcX896WBWRkIGHBAnhZI/rusloW8PrrMkcUQnR5OkOPa89wLvTaTXVnWyOZLjn7XFNS/M8z5OYC1dXccR06lJb727YB48dzN3XtWorvmjUMWHfuBOrrKbxHjwIA4o4cQb3H0hFLgxW0CiGE6E449XvaNCaGP/uMumxKbp36XVDAHdmcnNDnBQUhAybj7u+cz5qfzwBz8mSec/Agx9NVVrIs+L33gEGDYF12mecj+uryiBHAX/4C+FwnhBBdiTYLXHuUc2GkEl9DNLukkXpczT3M4HNzr0DAO4ubmQk88kjo+mXL2G9z/vl8X1rKPpdgkGVJs2fz97lzgdJS9AeiD1oTEii6cg0WQgjR3XDrt3XKV99y+Os79dsZxAYC1McDBxhEjh3L4wsXAq+8QuOl6dM5AaCujtVPBw5w7VGj2NZTVwfs2ePpGhwD4KTfcz/5JPDoo63wDyCEEJ2DNgtce5RzYTRBaaTd1KZwzoVLT+euqLmXO4vrd3/bZuZ20yY+Q3o63X6//JKfffObwNChsNat87y8H4Cg1wd9+3L3VkGrEEKI7ohbv03JsAlMI1VTrVnDoPTAAVY/XXUVP8/NZdBqWUwajxjBndUlS/j+6quBCROAhQthLVrk+2ieQWvv3kxMq59VCNHN6ChX4e7lXBhNUBppN9Xgt3PrnAtnHAqd9zQ2+pHuP2UKX48epSgmJNBE6eBBjsTZsgWxBQWel8bAJ2gdNQr4298UtAohhOi+uPXbHZi+/z7bbTZvBu64o7EmZmUBN9xA7TXBLsCqKDOuDmAS2cxr7dOHbT3BIHr5BK2+WfxHHgH++Z/lGiyE6JZ0VI9r93IujCYo9cIdqJqd26oq7p6a484gNTW1sSClpnJmnF/JsPO83FyK5JNPMtN7/DidBvv0geUTtPqK43PPKZsrhBCi5+CVXDb6vHAhsGIFf//Zzxqf5zZFMuPt0tJYQZWWBpxxBl/r6oDCQqCiApaPM7+vLj/2GOBRTiyEEN2FjnIV7jnOhZH6X90lxkYAg8HwYedNBcZ+JcNGIKdO5Qy6Dz5gZvirr4BevWCd9O6O8e2bGT4ceOklGT0IIYTo3vgll8vKuDtqjJUA4NZbuVN6xx2Ntd2YNTn9KcyonMsuAwYN4nkbN9LQaeRIWO+8A7zzTtjjROyVUtAqhOgBdAZX4e5NpP5Xd4mxeyi5GXY+axZ7aiIZQGVl0Xp/zx4KpBHcd98FXnuNGdy6OmDfPpYJJyZyyLnH8HJfcXz4YeDpp1vyryCEEEJ0XrwMEP2Sy+vWAcuX8/fzz+c5OTnAr38dWgsIrbFsGVBTQ+0FQjuwaWkMWI8dYwVUWRksMx82WqZMAf7P/9GoGyFEj0CBa1sTqf/Vr8TYHDe9L7bdtAFUaip3Q3fs4Nib3bvZFxsXR1fC7GwOLd+2DRg3Dta8eRyX48AC0OD3PZTNFUII0V0xGjtgQMgA0eh2RgZ7WU31U0YGk7/OHVenxrsdhtetY3vORReFks/33ccA97PPmKC+5hpYJhh2EFGXs7OBDz9sjW8vhBBdAgWubU00I278dlGnTAn1uhq8AmBnpjgnhzu1q1ezBOkHP6BjsJkZV14O6w9/CFsiFsBxr2fs2xf46U9lqS+EEKL7kpHBoHXqVLbpGF2eNi1kwATwvXOuutldraz01vPUVGDkSGDVKr6az8yYnA0bgC++gPXpp2GPpNJgIYRojALX9sIrSI1mF9V53JQRm8yve52cHIrvn/7EPlaAQaxtAwsWABs3wnr9dc/H8wxax47lyBwhhBCiO1NSwp3WYLDpth4n+fkcYTNmDNtxnN4URqP79eMonPJy4KGH2Ae7YwfwzDNAaSmsqqrmPavmswoheigKXNsLryC1JbNdvdYxmeKkJGDuXJYlnTzJgHXxYqCkhGYPHvhmdHNygPfei/65hBBCiK5KS9p6AM5hBTgebvhwTgVwa7Spnlq0iL2xdXW8dMMGzyXlMyGEEN4ocG0vvESxqTE6Xru0JkjNyAidZzLFf/87jR6GDAGOHKEh06FDsA4e9FzeUxwnTgQ+/7xZX00IIYTo9ERqz2nOrHWngdPkyY13WAMBBqmmLzYjg94SwSCrmGybzsFz5njewjdovf9+Ba1CiB6PAtf2oiWzXr12V02QWlISmudqguHFi1kiPGkSEBfnu8sK+Ijj/fcD8+Y17xmFEEKIrkBT7TnRXt+nD7B9Ox3/r7nGe62tW9m/OmAADZjKy4HevanNf/yj5/KeunzPPcD8+c1/ViGE6IYocO1I/LK/GzeyJ3XmTJbsZmWFGzB5lTMlJHDEzaFDsPLzPW+pEiQhhBA9kua257g12lxXVsbA1bZD5yQlAWvXUodXrKCx4dVXs/e1uhpIS4Plo7G+ujx1qoJWIYRwoMC1IzHZ26qqUGlRSQnNld56C6itBYwDsHE1zMkJN2ky65xzDjBpEqxXX/W8na84ZmcraBVCCNG9iVT5ZBLGd9zBlhkgfIfWOWt9+PDG+nv4MLB5M3X4yy8ZrE6axGtXrYL11luet/XV5TPOoMYLIYT4GgWu7UVxMcfR5OaGSnxNv+rRo6GSosOHgcGDafRw8cWh6529rYEAsHAh57Tu3AmsXw9cfDGsu+7yvLWvMA4cCFRUtOa3FEIIIboeCxYAr73G33/9a7767dA6A+CqKhoZfu97wOjR3CU189TT0oBt2zyD1oijbqZOVdAqhBAeKHBtL/LygKVL+buZ/2b6VTMzQ6NsSkr4eu21jcVy/XruwqalMdNbVMQypddeA1aubF7PDKAZcEIIIbo/Tc1LN8ycyWTwzJmhY84dVucYOpOIPu884Cc/AXbvBvr3B370I77268dznn8eVklJ2K18dfn660N/JwghhAhDgWtbY0Rz6lS+z80NfebM5hpBNbux5tWwdy+wfz9fp08H3n0X+PBDoLoaVk2N5619xTE3V0GrEEKI7k+0hkyWBUyYEBpv47WGaevJz2cf6+rVPP/MM4GrruI5AwYAeXmwli/3vI2vLo8dq6BVCCGaQIFrW2MELycntNNqqKwENm3iDqtXJtiZKb7xRiAxMRT4lpcDBw7AKigIu8xXGC+9lDNehRBCiJ5AtIZMkc5zt/WMGcOd2fPOC807v+wyjrz54ouWBa1btzb5VYQQoqejwLU1iFSKFEkMFy8GXnqJfaZXXRW63phETJjA0qVgENi1i7usmZnA//t/sF5+2fNRfIUxJycksEIIIURPoCWj6AxG26uq2NaTlsYA9uKLqcWBAHdJN2xgRdSOHbBWrw5bxleX09PZHiSEECIqFLi2BpFKkbxE04jhwIHAsGFATQ3w5z+HZsIZk4jaWmD2bFrv/9d/McAtKYH12Weej+EpjsrkCiGEEJFxlwM7HYPHjGHACjCA3bqVAWdVFXDwIBAMwvIZW+Opy6NGsQdWCCFEs1Dg2ho0dzacEcP0dFrnZ2R8LX5YsyZkDmFs+fPymN2NifEMWn2zuampClqFEEIIILrqqGAwlIg2JcJAYyNFc05KCqxFi3xv56nNEycCn39+2l9FCCF6IgpcW4NIpUheQmnEcN8+zn0bPhy4886QGKamskS4qornHz3qWX4ERAhak5IotEIIIYQIr45y67M5lpQU2nE9fBhISOBrejoDz+efB158ERg50vM2vrqcmKigVQghTgMFrm2Bc2ZrSQmFMhgMiaEZgzN2LDBoUGi2a3ExbfWXLaNrYXk5MGQIrFdf9byNSpCEEEIIB5F2VU1CNymJ752BrAlUs7JCiWh3knnRIrbyLFsGFBfDKipqtHzE2aypqUomCyHEaaLAtTVwC6VzZut55wFbtgApKXQRBkJiaAweDCagTUvjjNavvqJYuoho9KCgVQghRE8lkufE2rXU1LVruXPqbPMx15WVsXXHJJ4PHwb69AEaGoDt24EtW2Dt2RN224jVT6Z6SgghxGmhwLU1cAulGVmTmwvMnQt8+il3Vm+9tXH5UUlJ48DVZIPXrWN58I4dYbfyFceBA+VOKIQQomcTyXPCqc1A4zYf22aSuaKCGg0wibx0KU0U9+71DFgj0qePglYhhGhFFLi2Bu7yo8zM0MzWmTPZr/rd74YE0i2s8+YBv/oV56wGArB27gy7hW/A2rs3UF/fal9FCCGE6BJ4lQVH8pxwarP7+qVLmWS+4grqdm4uvSe2bweqq4FevTyX9NVmy+LEACGEEK1GTEc/QJfkgw+AG2/kK9C4/MiNZXEe6/79wG9/y1JeI6yVlTz2xBNAWRmspUthmZ3baBg2TEGrEEKInompdjI7pKdz/R13AN/+NvDggwxuMzOBhx8GLr0U1s6dnv2snkFrbCx3bxsaWvZMQgghfNGOa0v4zW+AlSv5+2WXhZcfOTG7quvWAcuXh87Ly2NJ0gcfAGeeyR3bEyfCLo9o9lBW1uKvIIQQQnRpmjuKLtL1+fnAxo2skJo4kceTklgu7CCiJtsRPxVCCHGaaMfVTSDAWaqBgP85998PfOtbfAVC5UemX9W5htldHT+eGdjMzMbmTTNnwvroI1iuoNU3mwtQHCWQQgghejJGX93uwU4iafrevcCSJXx95hngo4+Af/1XelI89RSsyy/3dfUPQ5oshBBtjnZc3URyJDRcdhl/mrPGypUUx5UrWYoEwDKBr4OI0tenj3pmhBBCiGiJpOkLFgCvvcbfH34YeOaZ0Mz0OXManaqwVAghOh4Frm5aUnrkNohwrmE+mzmTx+64A9i6FYg2i2tQNlcIIYRoHk49fv554MknGZTedRf1GODrhg0tc+aXNgshRLuhUmE3LSk9chtEONcwn736KrB4MfDBB7CuvTaU1T1FxNJgIYQQQjQmmtaet99mW8/bbzNoLSriKwCMHAnccAMwciSsu+/2NWDy1Ga17AghRLujHdeW4C49irRLa47deCNQXs6xNw6alD0JoxBCCBFONK09Tz7J+axmp9W8uq938LXqpqZGDoqFEEK0K9pxbQlZWUBOTigo/ctfgFtu4aubU9leq7wcFgArWidgZXOFEEIIf9xaDISPq5szhyPp5swBhg4FMjL46ne9oX9/76BV2iyEEB2GZXeh/wPOzs62P/30045+DFrmL1jAvpiJE4H0dI6mGT6cVvrFxXQOzs0Frr8e2LIFlmOmm+z0hRCic2BZ1gbbtrM7+jm6Mp1GmwEGrStX0vl/8eLGn11xBZ2DL70UWLUKlmV9/VEj5TU67Phc2iyEEO2HnzZrx7UlGCfCBQv4/sc/ZtD64x/zvRl3k5cHa9OmRkErACA5uT2fVgghhOh+bNwIPPQQXw3ucXXOHdiHH2bQ+vDD4Ws557VaVuOgVQghRKdAPa7R4HYNdjoRAkB1NVBXx1eAO63m1THyxgaAfv2AqqrweyibK4QQQvz/9u4/1u76ruP48yUbP9YGCNBtruVHSeoiTqxbN2cky8imskWpgIb6D6JOxAiziWQySUYiYjSdhkimyXTLYoLTGFycjDjEbKAmzBVEYIEtzRxZy1gLboXCoNv69o/vudzTe86991Da8/nee56P5Ob8+H7vOW8+ueWV9/fz/Xy/kxu+nc0tt3SPC29Xd+ut3QwszN+T9Z57GDnbbO/e7tFZVknqLWdcJ7HwqsGbN3chuXlz9/rDH4ann+4egZx7Lrn2WnLuuYd/ThU8++zoe4ajJEmLG3cF4SuvhK1b5w8iw+ga14UzsENq/fruqsHr1zvLKkkrgDOuk1ju3q7XXdc1rdddN7KpTj+9a2rBUJQk6UiMu4Lw3EHkYUMzrC/Nvg7PwA4bvm+rM62S1Hs2rpOYuy/rYq6/vvsZ56mnukdDUZKkxS1cljNs3AHkhRdKhPmZ1cHjYRdgOvFEeOEFOPFEBhuP/n+DJOmYsXGdxPBVgjduHNl8WDBWUWvWwHPPwZo1BqMkSZNY6r6s4w4gT7LGddh3vrP4d3tAWZJ6z8Z1EnNXCQa45prl9z9wYP65M62SJC1vuWU5Cy28UCLACSfAwYNw/PHdRROHmceStKLZuE5i+CrBAK99Lezb110+f+5KhMOcZZUk6eVZblnOQmPWuObgwe7JwYPdhZeGG9SF2WzzKkkrio3rJDZuPGymNfv2dU/27RsNxoUMRknSCpDkJmArcAjYC1xZVU+0rUqSpI63w5EkSQA7qur8qtoM3AF8qHVBS1q7tptFXbt28X3mbnOTzN9+zgPKkrQi2bhO4rTTutA77bTx24eDEQxHSdKKU1XPDL1cA/Q6xPLcc2TwOKeqXvqRJK0unio8gczd8Pxb31r+1GBJklaoJDcDVwD7gQsblyNJ0kuazLgmuSnJQ0keTHJXkje0qEOSpFmS5O4kj4z52QpQVTdU1ZnAbcCil9FPclWSnUl27pu77sOxtGFDd1bThg1L7+cZUJK0aqXF7GGSk+dOSUryfuC8qrp6ud/bsmVL7dy589gW95rXdPd6O+kkeP75uXpf2lxVXlJfklaJJPdX1ZbWdfRNkrOBz1TVsvemmUo2T5q75rMkrXiLZXOTGdder6OZu0H50I3KXTMjSVrtkmwaenkx8FirWkasX3/4oyRp5jRb4zrpOpokVwFXAZx11lnHvq6h57apkqQZ8sdJ3kh3O5zHgWXPhJqW7NnTPdmzZ+ls9gCzJK1ax2zG9Wito6mqj1bVlqrasm7duqNT3LZt3elE27ZNtr9rZiRJq1xVXVZVbxrcEufnq2rPVAu48UY45ZTucRILs1mStKo1WeN6WAEt1tEssQZmZD3rMvtLklYu17i+ckctm085BZ55Bk4+GfbvP2yT2SxJs6NXa1ybr6O5/PLDH4e4nlWSpAa2b++a1u3bRzaZzZKkVlcVvh04bB3NJKckTeXKhZKkmeGM6ytnNkuSjqbFsrnJxZmq6rIW3ztn7ClH8xsZ2jiliiRJmm1msyRpKU1OFZYkSZIkaVI2rpIkSZKkXmt2H9eWllzX6ylIkiRNndksSVqKM66SJEmSpF6byRlXL/IgSVLPmM2SpCU44ypJkiRJ6jUbV0mSJElSr83mqcKegiRJUr+YzZKkJTjjKkmSJEnqNRtXSZIkSVKv2bhKkiRJknrNxlWSJEmS1Gs2rpIkSZKkXrNxlSRJkiT1mo2rJEmSJKnXbFwlSZIkSb1m4ypJkiRJ6jUbV0mSJElSr9m4SpIkSZJ6zcZVkiRJktRrNq6SJEmSpF6zcZUkSZIk9ZqNqyRJkiSp12xcJUmSJEm9ZuMqSZIkSeq1VFXrGiaWZB/weOMyzgCealxDXzgWHcdhnmPRcRzm9X0szq6qda2LWMl6ks3j9P1vrwXHZJRjMp7jMsoxGXWsxmRsNq+oxrUPkuysqi2t6+gDx6LjOMxzLDqOwzzHQq34tzfKMRnlmIznuIxyTEZNe0w8VViSJEmS1Gs2rpIkSZKkXrNxffk+2rqAHnEsOo7DPMei4zjMcyzUin97oxyTUY7JeI7LKMdk1FTHxDWukiRJkqRec8ZVkiRJktRrNq5HIMlNSR5K8mCSu5K8oXVNrSTZkeSxwXh8KsmprWtqIckvJflSkkNJZu6Kc0kuSvLlJLuSXN+6nlaSfDzJ3iSPtK6lpSRnJvlckkcH/y5+p3VNmk3m9Shze9SsZ/gw83yU2T6qVc7buB6ZHVV1flVtBu4APtS6oIb+FXhTVZ0PfAX4YON6WnkEuBS4t3Uh05bkOOAjBvllIgAABWlJREFUwHuA84BfTnJe26qa+QRwUesieuB7wO9W1Q8Dbwd+e4b/JtSWeT3K3B41sxk+zDxf1Ccw2xdqkvM2rkegqp4ZerkGmNmFwlV1V1V9b/DyPmBDy3paqapHq+rLreto5G3Arqr6alUdBP4O2Nq4piaq6l7g/1rX0VpVfaOqHhg8fxZ4FFjftirNIvN6lLk9asYzfJh5PobZPqpVzr/qWH/BapXkZuAKYD9wYeNy+uLXgL9vXYSmbj3w9aHXu4GfaFSLeibJOcCPA19oW4lmlXm9JHNbw8xzvWzTzHkb10UkuRt4/ZhNN1TVP1XVDcANST4IXAPcONUCp2i5sRjscwPdaQO3TbO2aZpkHGZUxrw387MagiRrgduB7QtmvqSjxrweZW6PMsMnYp7rZZl2ztu4LqKq3j3hrn8LfIZVHITLjUWSXwF+DnhXreL7K72Mv4lZsxs4c+j1BuCJRrWoJ5K8mi7Mbquqf2xdj1Yv83qUuT3KDJ+Iea6Jtch517gegSSbhl5eDDzWqpbWklwE/B5wcVU937oeNfFFYFOSjUmOB7YBn25ckxpKEuBjwKNV9Wet69HsMq9HmdtagnmuibTK+czIgbajKsntwBuBQ8DjwNVVtadtVW0k2QWcADw9eOu+qrq6YUlNJLkEuBVYB3wbeLCqfrZtVdOT5L3ALcBxwMer6ubGJTWR5JPAO4EzgG8CN1bVx5oW1UCSC4B/Bx6m+/8kwO9X1Z3tqtIsMq9HmdujZj3Dh5nno8z2Ua1y3sZVkiRJktRrniosSZIkSeo1G1dJkiRJUq/ZuEqSJEmSes3GVZIkSZLUazaukiRJkqRes3GVjrEkB17B716TZFeSSnLG0PtJ8ueDbQ8lefPRqVaSpNXPbJZWHhtXqd/+E3g33f0Hh70H2DT4uQr4yynXJUnSrDKbpQZsXKUpGRyJ3ZHkkSQPJ7l88P4PJPmLJF9KckeSO5P8IkBV/XdVfW3Mx20F/qY69wGnJvnBMd95IMmfJLk/yd1J3pbk80m+muTiwT7HDer64uAI8W8O3l+b5N+SPDCod+vg/XOSPJrkrwY135XkpGMzapIkHTtms7Ry2LhK03MpsBn4MbojtTsGgXYpcA7wo8D7gJ+c4LPWA18fer178N5Ca4DPV9VbgGeBPwR+GrgE+IPBPr8O7K+qtwJvBX4jyUbgBeCSqnozcCHwp0ky+J1NwEeq6keAbwOXTVCzJEl9YzZLK8SrWhcgzZALgE9W1feBbya5hy6MLgD+oaoOAU8m+dwEn5Ux79WY9w4C/zJ4/jDwYlV9N8nDdIEM8DPA+XNHkoFT6MJvN/BHSd4BHKIL39cN9vnfqnpw8Pz+oc+SJGklMZulFcLGVZqecYG21PtL2Q2cOfR6A/DEmP2+W1VzoXkIeBGgqg4lmfv3H+DaqvrsYUUlVwLrgLcMAvVrwImDzS8O7fp9wNORJEkrkdksrRCeKixNz73A5YN1K+uAdwD/BfwHcNlgPc3rgHdO8FmfBq4YrM15O93pRN84wro+C/xWklcDJPmhJGvoju7uHQTjhcDZR/j5kiT1ldksrRDOuErT8ym6NTL/Q3fq0Aeq6skktwPvAh4BvgJ8AdgPkOT9wAeA1wMPJbmzqt4H3Am8F9gFPA/86iuo66/pTid6YLBOZh/wC8BtwD8n2Qk8CDz2Cr5DkqQ+MpulFSLzZypIaiXJ2qo6kOR0uiO9P1VVT7auS5KkWWU2S/3ijKvUD3ckORU4HrjJYJQkqTmzWeoRZ1wlSZIkSb3mxZkkSZIkSb1m4ypJkiRJ6jUbV0mSJElSr9m4SpIkSZJ6zcZVkiRJktRrNq6SJEmSpF77fwEu12IkqbnaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, a = plt.subplots(1, 2, figsize =(16, 8))\n",
    "a[0].scatter(lmg_c, lvg_c, c = \"r\", s = 2, alpha = .3)\n",
    "a[0].scatter(lmg_c, lvs_c, c = \"k\", s = 2, alpha = .3)\n",
    "\n",
    "a[1].scatter(lmg_s, lvg_s, c = \"r\", s = 2, alpha = .3)\n",
    "a[1].scatter(lmg_s, lvs_s, c = \"k\", s = 2, alpha = .3)\n",
    "\n",
    "a[0].set_title(\"ctrl dataset\")\n",
    "a[1].set_title(\"stim dataset\")\n",
    "for i in range(2):\n",
    "    a[i].set_xlabel(\"log10 mean\")\n",
    "    a[i].set_ylabel(\"log10 variance\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 395,
   "id": "0af7933e-005e-4e62-9868-ec370e8d85f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nsc.pp.normalize_total(adata, target_sum=1e4)\\nsc.pp.log1p(adata)\\n\\nstim_cat = adata.obs[\"stim\"].unique()\\nadata_list = [adata[adata.obs[\"stim\"] == sc] for sc in stim_cat]\\nadata_list\\n'"
      ]
     },
     "execution_count": 395,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "sc.pp.normalize_total(adata, target_sum=1e4)\n",
    "sc.pp.log1p(adata)\n",
    "\n",
    "stim_cat = adata.obs[\"stim\"].unique()\n",
    "adata_list = [adata[adata.obs[\"stim\"] == sc] for sc in stim_cat]\n",
    "adata_list\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e37eaa0-7047-4c3c-9b74-0e4150767035",
   "metadata": {},
   "source": [
    "### Seurat CCA\n",
    "We employ a variant of CCA, diagonal CCA, to account for cases where there are more cells than genes and apply this using the single-cell RNA-seq data sets as input"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71ed0288-10ae-48f8-a0f4-1ceea7d5cbcf",
   "metadata": {},
   "source": [
    "### MNN batch effect correction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c84a04eb-4be7-406a-920a-5c49cbb08d66",
   "metadata": {},
   "source": [
    "##### implement the mnn correction for batch effects first and then the seurat cca function\n",
    "\n",
    "The paper though introduces a lot of confusing and vague terminology and its a masterpiece in obfuscation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 406,
   "id": "838db4a1-0762-4c49-aec7-aa75847fb2cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def cosine_normalize(rc):\n",
    "    return rc / np.linalg.norm(rc, axis = 1)[:, None]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 407,
   "id": "67b46fed-68ea-4f35-abb8-597ea0d8c04d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1222\n"
     ]
    }
   ],
   "source": [
    "hvg0 = adata_ctrl_flt.var[\"hvg_sorted\"][:3000]\n",
    "hvg1 = adata_stim_flt.var[\"hvg_sorted\"][:3000]\n",
    "hvg = list(np.intersect1d(hvg0, hvg1))\n",
    "print (len(hvg))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 408,
   "id": "be9058a2-0ddf-4582-857e-8c18d3f7ff03",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_ctrl_hvg = adata_ctrl_flt[:,hvg]\n",
    "adata_stim_hvg = adata_stim_flt[:, hvg]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 409,
   "id": "ba673774-a041-4229-88b3-35653e051dda",
   "metadata": {},
   "outputs": [],
   "source": [
    "rc_hvg_b0 = cosine_normalize(adata_ctrl_hvg.X.toarray())\n",
    "rc_hvg_b1 = cosine_normalize(adata_stim_hvg.X.toarray())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 410,
   "id": "a35d2bff-c493-40ed-8709-947d09c355e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "gi = np.intersect1d(adata_ctrl_flt.var[\"features\"].to_list(), adata_stim_flt.var[\"features\"].to_list())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 411,
   "id": "a07e6106-01b1-49cd-afe7-f3dd16e2a874",
   "metadata": {},
   "outputs": [],
   "source": [
    "rc_b0 = cosine_normalize(adata_ctrl_flt[:, gi].X.toarray())\n",
    "rc_b1 = cosine_normalize(adata_stim_flt[:, gi].X.toarray())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 412,
   "id": "ca918ee9-4b02-4815-8987-9bae609ac53e",
   "metadata": {},
   "outputs": [],
   "source": [
    "rc_hvg_ref = rc_hvg_b0\n",
    "rc_ref = rc_b0\n",
    "\n",
    "\n",
    "bi_hvg = rc_hvg_b1\n",
    "bi = rc_b1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b792246-8b22-49f8-a9da-6d872cd7ffb8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "485dac97-e2f6-4808-8c4f-a5a288d85da7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 415,
   "id": "5e58fe8e-ca7b-479c-bb1d-e040fa183379",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KDTree\n",
    "from scipy.spatial import cKDTree\n",
    "from scipy.spatial.distance import cdist\n",
    "### might have to swap out this thingy for the fast nearest neighbor searches "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 416,
   "id": "91df444d-d17f-4b76-8519-84373f22e10a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_mnn(rc_hvg_ref, bi_hvg, n_neighbors):\n",
    "    \n",
    "    ## for each bi, compute the mnn between ref and bi\n",
    "    index_ref = cKDTree(rc_hvg_ref)\n",
    "    index_bi = cKDTree(bi_hvg)\n",
    "    dist_bi_to_ref, indices_bi_to_ref = index_ref.query(bi_hvg,    k = n_neighbors)\n",
    "    dist_ref_to_bi, indices_ref_to_bi = index_bi.query(rc_hvg_ref, k = n_neighbors)\n",
    "    mnn_ref = []\n",
    "    mnn_bi = []\n",
    "    for bi_node, ref_neighbors in tqdm(enumerate(indices_bi_to_ref)):\n",
    "        ## check if any of ref_neighbors are in indices_ref_to_bi\n",
    "        for rn in ref_neighbors:\n",
    "            if bi_node in indices_ref_to_bi[rn]:\n",
    "                mnn_ref.append(rn)\n",
    "                mnn_bi.append(bi_node)\n",
    "                \n",
    "    return mnn_ref, mnn_bi, indices_bi_to_ref, indices_ref_to_bi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 417,
   "id": "726f845e-f97b-4a49-91af-f24ac96ae14f",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_neighbors = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 418,
   "id": "72310493-1d53-4ff2-a3f3-b8aef98004ff",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "7451it [00:00, 17082.24it/s]\n"
     ]
    }
   ],
   "source": [
    "mnn_ref, mnn_bi, ind_bi_to_ref, ind_ref_to_bi = get_mnn(rc_hvg_ref, bi_hvg, n_neighbors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 419,
   "id": "2de02a52-c285-4db2-805b-31bc58939ddb",
   "metadata": {},
   "outputs": [],
   "source": [
    "### compute the batch vectors \n",
    "batch_vectors_hvg_ref_to_bi = rc_hvg_ref[mnn_ref] - bi_hvg[mnn_bi]\n",
    "batch_vectors_all_ref_to_bi = rc_ref[mnn_ref] - bi[mnn_bi]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 420,
   "id": "bf2edbbc-1858-446a-85e3-dedf43ccc656",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((11100, 1222), (11100, 12393))"
      ]
     },
     "execution_count": 420,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "batch_vectors_hvg_ref_to_bi.shape, batch_vectors_all_ref_to_bi.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 421,
   "id": "7fa4d68b-34f7-4f03-b91f-9d88445f83de",
   "metadata": {},
   "outputs": [],
   "source": [
    "## reduce the batch_vectors_hvg_ref_to_bi\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 422,
   "id": "d268bc3d-2baf-4eb6-98c0-b574f22b9428",
   "metadata": {},
   "outputs": [],
   "source": [
    "mnn_bi_unq, mnn_bi_cnt = np.unique(mnn_bi, return_counts = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 423,
   "id": "3c522318-997b-491e-b11e-ad742bf00a29",
   "metadata": {},
   "outputs": [],
   "source": [
    "### average the batch vectors at the same index in bi\n",
    "### i.e. have a shape of mnn_bi.shape[0] x hvg \n",
    "batch_vectors_hvg_ref_to_bi_red = np.zeros((bi_hvg.shape[0], batch_vectors_hvg_ref_to_bi.shape[1]))\n",
    "for index, ve in zip(mnn_bi, batch_vectors_hvg_ref_to_bi):\n",
    "    batch_vectors_hvg_ref_to_bi_red[index] += ve\n",
    "batch_vectors_hvg_ref_to_bi_red = batch_vectors_hvg_ref_to_bi_red[mnn_bi_unq] / mnn_bi_cnt[:, None]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 462,
   "id": "9af71314-61e2-44c2-82a4-e5a57536d6c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_mat = np.zeros((len(bi_hvg), len(kernel_cells_bi)))\n",
    "for i, bh in tqdm(enumerate(bi_hvg)):\n",
    "    for j, kc in enumerate(kernel_cells_bi):\n",
    "        dist_mat[i, j] = np.dot(bh , kc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3c72187-d434-4e9a-a7fc-9b361948bd0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "## total contribution of a kernel point to the entire dataset is normed to 1 \n",
    "## so if thats too far away it doesnt adversely affect the results "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 463,
   "id": "af0a30d7-1090-43d7-bd14-de5b1bd9514e",
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel_cells_bi = bi_hvg[mnn_bi_unq]\n",
    "cdist_kernel_bi = dist_mat\n",
    "\n",
    "sigma = 1.0\n",
    "W = np.exp(-cdist_kernel_bi / sigma**2)\n",
    "norm_W = np.sum(W[mnn_bi_unq], axis=0)\n",
    "W = np.divide(W, norm_W)\n",
    "denom_Wm = W.sum(axis = 1, keepdims = True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 469,
   "id": "a7e8eed7-64dc-434a-8d80-1e50d2024d34",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((7451, 1833), (1833,), (1833,))"
      ]
     },
     "execution_count": 469,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W.shape, norm_W.shape, mnn_bi_unq.shape### across all kernel cells the probability is normalized to 1 \n",
    "### and across all cells its normalized to 1 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 492,
   "id": "26ed801d-d335-4e43-b988-405ac8c651a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "correction = np.divide(np.dot(W , batch_vectors_hvg_ref_to_bi_red), denom_Wm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 497,
   "id": "f65f4328-c60e-47c1-86e7-768eac12f966",
   "metadata": {},
   "outputs": [],
   "source": [
    "bi_corrected = bi_hvg + correction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3d39341-ff92-4d91-be71-23419a27d299",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6652b07b-c7ca-4469-b538-c4118846341a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07c773fa-93f9-48a8-a7b7-3f50beb014b2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c13a2c60-3e01-4acd-b365-d44f28040bfd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 434,
   "id": "dfe8e293-f6e9-43fb-b246-e6922a1cc53b",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Performing cosine normalization...\n",
      "Starting MNN correct iteration. Reference batch: 0\n",
      "Step 1 of 1: processing batch 1\n",
      "  Looking for MNNs...\n",
      "  Computing correction vectors...\n",
      "  Applying correction...\n",
      "MNN correction complete. Gathering output...\n",
      "Packing AnnData object...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/user-dnf/.local/lib/python3.8/site-packages/scipy/sparse/_index.py:116: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.\n",
      "  self._set_arrayXarray_sparse(i, j, x)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done.\n"
     ]
    }
   ],
   "source": [
    "corrected = mnnpy.mnn_correct(adata_ctrl_hvg, adata_stim_hvg, do_concatenate = False, var_adj=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 483,
   "id": "70aca6cd-61e5-4c9a-bf74-8895110bb1ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "corrected_mnnpy = corrected[0][1].X.toarray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 448,
   "id": "df238e2d-da1a-4c88-be91-729aad10581c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.00001207, -0.00002247,  0.00030707, ...,  0.00003268,\n",
       "        -0.00001538,  0.00010187],\n",
       "       [ 0.00000916, -0.00002378,  0.00034182, ...,  0.00003332,\n",
       "        -0.00000494,  0.00011209],\n",
       "       [ 0.00001171, -0.00002378,  0.00033741, ...,  0.0000445 ,\n",
       "        -0.00001758,  0.00012099],\n",
       "       ...,\n",
       "       [ 0.00000989, -0.00002577,  0.00035307, ...,  0.00003647,\n",
       "        -0.0000072 ,  0.00012651],\n",
       "       [ 0.00001008, -0.00002174,  0.00032928, ...,  0.00005182,\n",
       "        -0.00003217,  0.00014998],\n",
       "       [ 0.00001266, -0.00002262,  0.00035557, ...,  0.00004068,\n",
       "        -0.00001323,  0.0001155 ]], dtype=float32)"
      ]
     },
     "execution_count": 448,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7dca8b3-dd98-4ff7-ad81-9ba8a89b11b5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03c72bdc-5d1e-4637-a74c-3692f5241ebb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fdda0cb-e260-4181-9f34-87e9237187c3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70bfb348-4f6f-482a-9a2b-3a4bb100eada",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d00f6986-33ef-473a-a161-fd727bcadd3c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0195a2bd-2541-413f-90da-07537332488e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "hiplex",
   "language": "python",
   "name": "hiplex"
  },
  "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.8.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}