Cortex-imputation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imputation results of scPI-based methods (Cortex)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "plt.style.use('ggplot')\n",
    "import time\n",
    "from scPI import scPI_FA\n",
    "from scPI import scPI_ZIFA\n",
    "from scPI import scPI_ZINB\n",
    "from utils import *\n",
    "from sklearn.decomposition import FactorAnalysis\n",
    "from sklearn.metrics import silhouette_score\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.metrics import normalized_mutual_info_score as NMI\n",
    "from sklearn.metrics import adjusted_rand_score as ARI\n",
    "def cluster_scores(latent_space, K, labels_true):\n",
    "    labels_pred = KMeans(K, n_jobs=8, n_init=20).fit_predict(latent_space)\n",
    "    return [silhouette_score(latent_space, labels_true), NMI(labels_true, labels_pred), ARI(labels_true, labels_pred)]\n",
    "n_celltype = 7\n",
    "n_repeat = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def imputation_error(X_mean, X, X_zero, i, j, ix):\n",
    "    \"\"\"\n",
    "    X_mean: imputed dataset\n",
    "    X: original dataset\n",
    "    X_zero: zeros dataset\n",
    "    i, j, ix: indices of where dropout was applied\n",
    "    ========\n",
    "    returns:\n",
    "    median L1 distance between datasets at indices given\n",
    "    \"\"\"\n",
    "    all_index = i[ix], j[ix]\n",
    "    x, y = X_mean[all_index], X[all_index]\n",
    "    return np.median(np.abs(x - y))\n",
    "\n",
    "X = np.loadtxt(\"/Users/jiazhao/Downloads/Results/Cortex/data_train\", dtype='float32')\n",
    "X_zero = np.loadtxt(\"/Users/jiazhao/Downloads/Results/Cortex/imputation/X_zero.txt\", dtype='float32')\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/imputation/\"\n",
    "i = np.load(data_path + \"i.npy\")\n",
    "j = np.load(data_path + \"j.npy\")\n",
    "ix = np.load(data_path + \"ix.npy\")\n",
    "n_repeat = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### example of performing scPI-FA, scPI-ZIFA, scPI-ZINB (vary n_epoch and linearity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When evaluating the impact of the number of epochs, we vary n_epoch to be $\\{250,500,10000\\}$. When evaluating the impact of the linearity and non-linearity, we vary non_linear to be False or True. We repeat the experiments ten times to draw boxplots. Note: we use learning rate 1e-3 for scPI-ZIFA-linear, and learning rate 5e-4 for all the other methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_zero = np.log(np.loadtxt(data_path + \"X_zero.txt\", dtype='float32') + 1)\n",
    "n_repeat = 10\n",
    "file_name = \"scPI_FA\"\n",
    "for k in range(n_repeat):\n",
    "    result, result_test = scPI_FA.fitModel(Y=X_zero, K=10, batch_size=128, n_epoch=500, lr=5e-4, non_linear=True, Y_test=None)\n",
    "    A = result[\"A\"]\n",
    "    Z = result[\"latent\"]\n",
    "    mu = result[\"mu\"]\n",
    "    X_impu = np.exp(Z @ A.T + mu.reshape([1, D])) - 1\n",
    "    impu_err = imputation_error(X_impu, X_true, np.exp(X_zero)-1, i, j, ix)\n",
    "    np.savetxt(str(k) + file_name + '_imperr.txt', np.array(impu_err).reshape([-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_repeat = 10\n",
    "file_name = \"scPI_ZIFA\"\n",
    "for k in range(n_repeat):\n",
    "    result, result_test = scPI_ZIFA.fitModel(Y=expression_train, K=10, batch_size=128, n_epoch=500, lr=5e-4, act_nng_exp=True, non_linear=True, Y_test=None)\n",
    "    # save training results\n",
    "    A = result[\"A\"]\n",
    "    Z = result[\"latent\"]\n",
    "    mu = result[\"mu\"]\n",
    "    X_impu = np.exp(Z @ A.T + mu.reshape([1, D])) - 1\n",
    "    impu_err = imputation_error(X_impu, X_true, np.exp(expression_train)-1, i, j, ix)\n",
    "    np.savetxt(str(k) + file_name + '_imperr.txt', np.array(impu_err).reshape([-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_repeat = 10\n",
    "file_name = \"scPI_ZINB\"\n",
    "for k in range(n_repeat):\n",
    "    result, result_test = scPI_ZINB.fitModel(Y=expression_train, K=10, batch_size=128, n_epoch=500, lr=5e-4, non_linear=True, act_nn_exp=True, Y_test=None)\n",
    "    # save training results\n",
    "    X_impu = result[\"rho\"]\n",
    "    impu_err = imputation_error(X_impu, X_true, expression_train, i, j, ix)\n",
    "    np.savetxt(str(k) + file_name + '_imperr.txt', np.array(impu_err).reshape([-1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Impact of the number of epochs on imputation error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ZIFA \n",
    "n_repeat = 10\n",
    "zifa = np.zeros([n_repeat,])\n",
    "N, D = X_zero.shape\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/ZIFA/\"\n",
    "A  = np.loadtxt(data_path + \"0impu_CORTEX_ZIFA_A.txt\")\n",
    "mu = np.loadtxt(data_path + \"0impu_CORTEX_ZIFA_mus.txt\")\n",
    "Z  = np.loadtxt(data_path + \"0impu_CORTEX_ZIFA_Z.txt\")\n",
    "Xtilde = Z @ A.T + mu.reshape([1, D])\n",
    "imp_err = imputation_error(np.exp(Xtilde)-1, X, X_zero, i, j, ix)\n",
    "for t in range(n_repeat):\n",
    "    zifa[t] = imp_err\n",
    "\n",
    "# ZINB\n",
    "zinb = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/ZINB/\"\n",
    "for t in range(n_repeat):\n",
    "    Xtilde  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_ZINB_Xtilde.txt\")\n",
    "    imp_err = imputation_error(Xtilde, X, X_zero, i, j, ix)\n",
    "    zinb[t] = imp_err\n",
    "\n",
    "# FA\n",
    "fa = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/FA/imputation/\"\n",
    "for t in range(n_repeat):\n",
    "    fa[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_FA_imperr.txt\")\n",
    "\n",
    "# zifa-nln\n",
    "zifa_nln = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/scPI_ZIFA/500-impu/\"\n",
    "for t in range(n_repeat):\n",
    "    zifa_nln[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_ZIFAVAE_nonlinear_imperr.txt\")\n",
    "\n",
    "# zifa-nln\n",
    "zifa_nln_1w = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/scPI_ZIFA/10K-impu/\"\n",
    "for t in range(n_repeat):\n",
    "    zifa_nln_1w[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_ZIFAVAE_nonlinear_imperr.txt\")\n",
    "\n",
    "# zifa-nln\n",
    "zifa_nln_250 = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/scPI_ZIFA/250-impu/\"\n",
    "for t in range(n_repeat):\n",
    "    zifa_nln_250[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_ZIFAVAE_nonlinear_imperr.txt\")\n",
    "\n",
    "# zinb-nln\n",
    "zinb_nln = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/scPI_ZINB/500-impu/\"\n",
    "for t in range(n_repeat):\n",
    "    zinb_nln[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_ZINBVAE_nonlinear_imperr.txt\")\n",
    "\n",
    "# zinb-nln\n",
    "zinb_nln_1w = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/scPI_ZINB/10K-impu/\"\n",
    "for t in range(n_repeat):\n",
    "    zinb_nln_1w[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_ZINBVAE_nonlinear_imperr.txt\")\n",
    "\n",
    "# zinb-nln\n",
    "zinb_nln_250 = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/scPI_ZINB/250-impu/\"\n",
    "for t in range(n_repeat):\n",
    "    zinb_nln_250[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_ZINBVAE_nonlinear_imperr.txt\")\n",
    "\n",
    "# fa-nln\n",
    "fa_nln = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/scPI_FA/500-impu/\"\n",
    "for t in range(n_repeat):\n",
    "    fa_nln[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_FAVAE_nonlinear_imperr.txt\")\n",
    "\n",
    "# fa-nln\n",
    "fa_nln_1w = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/scPI_FA/10K-impu/\"\n",
    "for t in range(n_repeat):\n",
    "    fa_nln_1w[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_FAVAE_nonlinear_imperr.txt\")\n",
    "\n",
    "# fa-nln\n",
    "fa_nln_250 = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/scPI_FA/250-impu/\"\n",
    "for t in range(n_repeat):\n",
    "    fa_nln_250[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_FAVAE_nonlinear_imperr.txt\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_repeat = 10\n",
    "model = list(\"FA\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-FA-10K\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-FA-500\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-FA-250\" for _ in range(n_repeat)) + \\\n",
    "        list(\"ZIFA\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-ZIFA-10K\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-ZIFA-500\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-ZIFA-250\" for _ in range(n_repeat)) + \\\n",
    "        list(\"ZINB\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-ZINB-10K\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-ZINB-500\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-ZINB-250\" for _ in range(n_repeat))\n",
    "error = list(fa) + list(fa_nln_1w) + list(fa_nln) + list(fa_nln_250) + \\\n",
    "        list(zifa) + list(zifa_nln_1w) + list(zifa_nln) + list(zifa_nln_250) + \\\n",
    "        list(zinb) + list(zinb_nln_1w) + list(zinb_nln) + list(zinb_nln_250)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Cortex')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAGPCAYAAADWY0OgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmcHVWd//9XJaQ7Gr4YCAGiLM0MMDMCItuMBAUlMAZwi8sH0EQCYaTzVQhfhCCbbCKmHZgJWxoISwIqfBxWSdAfyYBsIhLZgiKKhBBoIAQCCZruDl2/P07dTuV23b51e7u9vJ+Px33ce889dep0Ba1PnfqcU1Ecx4iIiIgUG1btDoiIiEj/pCBBREREMilIEBERkUwKEkRERCSTggQRERHJpCBBREREMilIEBERkUwKEkRERCSTggQRERHJpCBBREREMm1S7Q6IyNATRdFw4CvA54BPAFsBHwRWA88DDwI/ieN4aRX6NhWoA+6P4/j+vt6/SH8S6dkNItKXoij6BDAP2CVV3AqsAUaz8QjnbcBRcRy39GH/7gcOBM6L4/jcvtqvSH+k2w0i0meiKPo8cD8hQFgFnA7sEsdxTRzHY4AaYF/gR8C7wJcJIwwiUgW63SAifSKKop2Bm4Ba4A/AZ+M4XpGuE8fx+8DjwONRFP0YuK7POyoi7TSSICJ95QfAZsA6YFJxgFAsjuO34jj+EvBOujyKom2iKPpxFEXPRlG0Noqi95LPDVEUbZ3VVhRFdVEUxcmrLoqif4yi6Oooil6Moqg5iqJlURRNjaIoJtxqADgntU37thlt7xlF0XVRFL0QRdHfkj49FUXRD6Io2jKj/v5RFK1P2jupRH+3jaJoVVLn6s6Ok0hvUk6CiPS65OT9KuHC5No4jo/rYjsHAncQchcA/gbEwKjk+9vAF+I4fqhouzrgxeTrN4CrgE1T278JnAbMBrYARgDvAWuLurBvHMcvp9o9DzgbiFL92YRw2wSgCTg8juMnivpzNnA+0AJ8Iv17FEXDgP8lBCvPAXvHcfy3Tg+MSC/RSIKI9IXPsOH/b27vSgNRFG3HhgDhD8An4zgeFcfxpsABwJ+AzYE7oyj6SCdNXQU8SzjhF7b/9ziOb4njeBvgkaTef8ZxvE3RKx0gnAR8nxBInA6Mi+N4FCGHYh/CiX4ccFcURZsW9eFCQm5GDXBzFEWjUr+dRQgQmoEjFSBINSlIEJG+sGvq8xMla3XuDEKA8DYwIY7jhws/xHH8IHAwIdlxC8JJu5RVwMFxHD+e2v75SjqS3Ea4kDAKMSmO4x/Fcfxa0tb7cRwvAT4LLAG2BTYaOYnjuA2YnPRlF+DypN39CYEHwMw4jp+qpF8iPU1Bgoj0hTGpz29VunEURRFgydfGwgk5LclxaEy+HtlJc5fHcVx8G6FS3yCMGDwex/HirApxHK8HfpZ8/WzG768AxyZfp0ZR9H+BnwLDgQVxHF/azT6KdJtmN4hIX4jKV+nUjoQRAoBFndS7F5gJjImiaMc4jl/MqPNwRlmlPpm87xZFUYeAJeUDyfsOWT/GcXxXFEWXA98BrkiKm4CpPdBHkW5TkCAifeHN1OctCEmMldgq9fmVTuqlZ0xsxYZkxbQ3Ktx3lg8n7x9gQyDQmc7WejgFmAQU8iiOjeP4zU7qi/QZ3W4Qkb7wbOrznt1sK++UrFL13u/m/iHcEoBw6yPK8arrpK3D2RAgwIYpmCJVpyBBRPrCfUBb8nlSF7ZPX/1v10m9bVOfV3ZhP3kVbjHs3p1Gkhkbc5OvTyfvM6MoOqg77Yr0FAUJItLr4jh+Hbg1+fr1KIp26ax+WpK0+CIbEh4ndFL94OR9VYl8hDwKwUxneRSFvIZPRFGUmW9QTvKQq58Qpm3+gfCgq9sJ/798YxRFYzrZXKRPKEgQkb5yFmFNgQ8At5VZy4AoijaPouhW4ENxWPXtluSn46Mo2iaj/oeB45OvPyv+vQLvJu+jO6lzI/B3wm2HK5ITfqYoioZFUZTV1lnApwjrIRwVx/HfCVMlVxByHq7vQt9FepSCBBHpE8laBFMIqwzuCjwZRdFpURTtVKgTRdHwZJnj84G/Eh7wVPBDwqOktwAWRVE0PrXd/oRZD6MJIw4/6kZXC4+nPqxUIJNMwfxe8vVw4N5kueXhSX+iKIr+OYqik5P2PpfePunv2cnXU+M4fjpp9y3C+gltwOejKPpON/4OkW7Tsswi0qeSE+QNwE6p4hbCKEP6UdExcDNwdBzHrcm2BwJ3Ah9K6ryXvBdWLFxNWJb5waJ91rFhpsOOcRwv66R/OxPyA0YSTtYrCc+bgLDK44pU3VOBi9iQyNhCeOT1ZoSlnQsmx3H8k2Sb0cBTwPaE9RA2CiCSOucTgoh1wL/GcfxMqf6K9CYFCSLS55Ir7q8RrrD/jTBd8QOEhzk9B/wauDGO4z9lbLsNYdrgYUAdIZh4CVgAXJy10FIlQUJS/xOEVRv/jbAQVGG6eIdtk5GQ/0vIlagjBCzvAn8hJGzeATya3DIhiqKfA18lJD9+LI7jDgmWyfF5ABhPyFfYJ7kdIdKnFCSIiIhIJuUkiIiISCYFCSIiIpJJQYKIiIhkUpAgIiIimRQkiIiISCYFCSIiIpJJj4qWcsYAnwWWsWFBGRERGVhGEtbx+BWwKu9GChKknM8SHkIjIiID3zeAn+atrCBBylkG8Pbbb7N+/fqqdmTs2LGsXNmbT/8dGHQcdAwK0sehpaWFH/7whwCcccYZ1NTUVLNrfUb/LQTljsMmm2zC5ptvDsn/p+elIEHKWQewfv16Wltbq9aJ8LTg0I+hvEqojoOOQUHxcVizZg1/+ctfAFizZg2bbbZZNbvXJ/TfQlDhcajotrESF0VERCSTggQRERHJpCBBRGQQaGtry/ws0h0KEkREBoH33nsv87NIdyhxcYAxs6nA9SV+vtjdT0nVrQFeAbYEJrn7Hb3fQxERGSwUJAxcxwDPFZW9WvT9i4QAAWAaoCBBRERyU5AwcC1198fL1JkGNAMPAYea2Th3b+r9romIyGCgnIRBysy2Aw4hjB5cDAwHjq5qp0REZEDRSMLANdzMNvr3c/f0kojHEILA64BFwArgWOBHnTVqZrVAbeF7Q0PDqLq6uh7qsoiIDCQaSRi4HgVa069C0GBmESFIWA4scvc2YB6ws5kdUKbd04F3Cq85c+Y80Ev9FxGRfk5BwsD1TWDf9Cs1kjCB8LSvG5IAAcKMiJgwmtCZi4APFV7Tp08vF1SIiMggpdsNA9cfO0lcnJa8325mo5PPqwijD18zsxPd/d2sDd29mZDsWKAJ1yIiQ5SChEHGzDYHJiVfnyhR7Ujg6r7pkYiIDFQKEgafyYTEwzOA3xT9FgG3Em45KEgQEZFOKUgYfI4F3gQuSW4dbMTMbgRONLNd3f3ZPu+diIgMGEpcHETMbG/g44SExQ4BQqIwgjCtxO8iMgCNGjUq87NId2gkYYBx9xuAG0r8toRwS6Gz7Z8tV0dEBp5hw4ZlfhbpDv2XJCIiIpkUJIiIiEgmBQkiIoNAbW1t5meR7lBOgojIIFBTU0NjY2P7Z5GeoCBBRGQQiKJIIwjS4xQkiIgMYnEc09LS0ml5TU0NUZQ96amz32TwU5AgIjKItbS0UF9f3+XtGxsbNUIxhClxUURERDJpJEFEZIi44OUV1MQxAC1RxNnbbduhvPg3GdoUJIiIDBE1cUxtKhgoVy6iIEFEZBDIm4g4UPcn1aGcBBGRQaCQoFhfX585m2Gg70+qY8CPJJjZucA5qaJW4BXgLuAcd1+d1JsKXA/s6+6Pl2mzUDfLxe5+SlH9mmSfWwKT3P2OCvr/IeBMYM/kNQY4291/UKL+PsAs4N8If+v/At9192WpOjsBfwb+n7v/d6p8OHAd8E3ge+4+K28/RURk6BlMIwkTgf2Aw4E7gBOAe8ysO2NgxyRtpl+XZtT7IiFAgMofwTwW+BZQQ+h3SWa2K3A/MBz4GnAc8C/Ag2Y2psy2tcCtwGTgWwoQRESknAE/kpCyxN3fTD7fm5w0pwDjgYe72ObScqMOiWlAM/AQcKiZjXP3ppz7eMHdRwOY2TZ0HmRcAKwFPufua5NtngD+BJxMGJHowMw2Be4EPgkc5e6es28iIjKEVTVIMLOxwIXAocBWwLvA84TbBIuSOhOBU4F9gBHAS8B8d7+oTPOPEoKEHeh6kFCWmW0HHAL8HJgHTACOBn6UZ3t3z5VSnIwEHAbMLQQIyfZ/NbMHgElkBAlJsHQPsCvwBXf/VZ79icjAEqdmJzQ3N2d+zjt/IV0vvT0lymPNjBi0qj2ScCOwF+Hk9jwwOvk+BsDMpgHXAL8G6oE3gF2A3XK0vVPyvrIb/RtuZhsdI3dfX1TnGMJtm+uARcAK4FhyBgkV2AmoBZ7O+O1p4EAzG+HuranyjwAPAB8GDnH3R8rtJAlG2pdXa2hoGFVXV9edfotIH0gnD86YMSOzTmsUMTLHCb01NVOhVFvF+x45cmSOXspAU+0gYX/ClfE1qbI7oX2I/BLCKMBBqSvuxSXaKpzQNyXkJdQDLwMPdqN/jxYXJCfi9cnniBAkLAcWuXubmc0DzjSzA9z9gW7su1gh5+CtjN/eIuQpjGbjoKiQYHlQngAhcTqpRNA5c+Ywa5bSF0REhqJqBwmPAVPNbBXhKnxJ6kp4PLAZcGXOIfnXir4/TEjQW5dVOTnBD0+XZYwSfBP4Yyd1JgB1wPnu3paUXQ+cQRhNeKCCfeXV2bEo/u2XwEHAf5nZwamcjc5cRAjOAJg+ffoeJH+HiPRf6cdDz549u/15C83Nze2jASNy3hZI10u3lZZuV4+mHryqHSQcAZxFyNK/AFhrZrcDMwlZ/xCG7/M4GHiHMC1whbuvKlP/QOC+dIGZ7ZieSgj8sUziYiHJ8HYzG518XkUYgfiamZ3o7u8Sgol7i/a1nbvn/dsK7cKGEYW0LYD3CX9/2q+AKwizGhab2YRygYK7NxOSMAveq6CPIlIl6cWMamtrM0/sead6peuVaqvUvmVwqWqQkJywTgJOMrPtgS8Q7uVvxYar2bwLiD+V80q5YAmwb1HZq3k3NrPNCcmCAE+UqHYkcDXw24x9vZ53X4k/E07eu2f8tjvwfFE+AgDufreZfZkQKNxnZge5e3fyNEREZIio9khCO3dfDlxuZhMIuQqPEK6M683s5ryzACrY3xogz/TGUiYTEvzOAH5T9FtEOCkfC1zdA/vC3VvMbAHwVTP7nru/B2BmdYRRkZKJA+6+wMwmAbezIVB4ozv9ERGRwa9qQUKy0uB9wE+B54A1hKvticBt7r7WzL4LzAUWmdk1hKvvnYA93P071el5u2OBN4FLkiH6jZjZjcCJZraruz/bWUNmdhjwQULiIcCuZvbV5PMCd/978vn7hFGJu81sVrLNBYR8jP/qbB/ufo+ZfZGwYFMhUKh0NENERIaQaq64uI5wwpsC/IQwl/84whXxfwC4+7WEtQGGE4KFuwm3J5ZXob/tzGxv4OPADVkBQuLq5D3PCoxXE9ZZKMzyODL5/nNSOQhJsPEZQoLirYQkyeeAA3LkYJCskfBF4B8IgcI2OfomIiJDVKRFMKSMvYAlK1eupLW1Q8pDn4miiHHjxtHU1DSkF27RcdAxKCg+DqWeytjc3Ex9fT0As5a/3P5I6OYo4rTtt+tQXvxbY2NjZuJif3gKpP5bCPIchxEjRjB27FiAvYHf52273+QkiIhI10VRVHYWwkDen1SHggQRkSGiJXW1X+pz1ncZuhQkiIgMEWdvlz2jvFS5yGB6VLSIiIj0II0kiIgMYjU1NTQ2NnYoz5t4qCWXhzYFCSIig1hnCYZ6cqOUoyBBRGQAi+OY5uZSy7Xk276npzJWa0qk9DwFCSIiA1hLS0v7Ogj9Ram1FWTgUeKiiIiIZNJIgojIIHHByyuoqXDlwZYoap8C2ZXts9qRwUNBgojIIFETxxstr9zX28vgoyBBRGQA6A/PSujvdIx63qAIEszsXOCcVFEr8ApwF3COu69O6k0lPDlxX3d/vEybhbpZLnb3U4rq1yT73BKY5O53VND/g4F7S/zcoa9m9u/A+cAewFrgF8BMd38zo09nAkcD2wAvApe5+5V5+yYi/UM6QbGxsVHTFzMUHyMlT3bfYEtcnAjsBxwO3AGcANxjZt0JJ49J2ky/Ls2o90VCgAD5Hg+d5bSMff0xXcHMJgALCQHJ54H/B3wWWJQEBWlXATOBywjH5i7gcjOb2cX+iYjIEDIoRhJSlqSupu81szHAFGA88HAX21xabtQhMQ1oBh4CDjWzce7eVOG+nnf3R8vU+THwLGDu/j6AmS0Hfk0YMbgmKdsDmEoYYbg42fZ+MxsLnGNmV7n7OxX2T0REhpCqBwnJSetC4FBgK+Bd4HnCbYJFSZ2JwKnAPsAI4CVgvrtfVKb5RwlBwg50PUgoy8y2Aw4Bfg7MAyYQTtg/6uH97ADsCZxaCBAA3P0BM3sBmEQSJABfSt5vKGrmesLoyGcB78n+iYjI4FL1IAG4EdiLcO/8eWB08n0MgJlNI5z4fg3UA28AuwC75Wh7p+R9ZTf6N9zMNjpO7r6+qM4xhFs31wGLgBXAsVQeJFxlZv8D/A14BDjf3R9J/V74m5/O2PYZYO+iuk3uXvy3P536XUGCyAARp2YdNDc3E0UR69at22i1xWrOS0jvuzsrQFYifQyKV56MNUujR/SHIGF/YK67X5MquxPAzDYFLiGMAhzk7oV/9cUl2iqc0Dcl5CXUAy8DD3ajfx2G/81sRCFQSPIdjgGWA4vcvc3M5gFnmtkB7v5Ajn2sBv6LEAi9RQiCTgV+bWYT3b3w945J3t/KaOOt1O+Fuh3qufs7ZvZ+Ud3031YLtGf7NDQ0jKqrq8vxJ4hIbypk7QPMmDEjs05rFDGySifH1tRMglL960stLS1K7uwB/SFIeAyYamarCFfhS9y9NfltPLAZcGUqQOjMa0XfHwa+5e7rsionJ/jh6bKMUYJvUpQ8WFRnAlBHuOpvS8quB84gjCY8UG5fSc5DOu/hQTO7jZB70MDGIwRQ+oKhuLxUvaiT304nNVNkzpw5zJo1q0RVEREZzPpDkHAEcBZwHHABsNbMbidk5Y9N6qzI2dbBwDuEKZAr3H1VmfoHAvelC8xsR3dflir6Y5nExcJMhtvNbHTyeRVhBOJrZnaiu79LCCY2muZoZtu5e+bf5u5vm9kCYJqZ1bh7S9IuZI8CbMHGIwergH8urmRmHyLcGskajQC4iDB6A8D06dP3IAl0RKR60o9snj17NiNHjmSbbbbhpZde4sQTTwRgRBWH2NP7nj17dp9MP4yiiG222YbXXnut/XZDYRRDj7juGVUPEpLZCCcBJ5nZ9sAXCPfyt2LDySrvWp9PFa8VUMYSYN+islfzbmxmmxOSBQGeKFHtSOBq4LcZ+3q9zC4K43eF//U9k7zvDvx/RXV3B5amvj8DfMXMtiw6Jrsn7+m67dy9mTBLo+C9Mn0UkT6QXhiotraW2tpaRo4cudHJuJpLB6X3Xehfr+8zitqPQXEOghZS6hlVDxLS3H05YR7/BEKuwiOEkYF6M7s55y2HSva3ho2H+Ss1mXD//gzgN0W/RcCthFsOV1e6LzPbgpBX0X77xd2Xm9nvgSlm9t+pKZCfBP6RjRMl7wDOI8yyuDhVPpWQGPmrvH0REZGhqapBQjL0fR/wU+A5YA3hansicJu7rzWz7wJzCYsFXUO4+t4J2MPdv1Odnrc7FngTuCS5At+Imd0InGhmu7r7s6UaMbNbgL8SRjbeJCQunkK4rfD1ouozCSf4W8xsDmElxVnAU8D8QiV3f9rMbgB+YGZx0vahSZ+/pzUSRESknGqvuLiOMAw/BfgJcA8hN2EW8B8A7n4tcBgh6W8ucDfh9sTyKvS3nZntDXwcuCErQEhcnbyXW4HxacIJfC4hefMHSdl+7r5RzkQy0+FwYDvCsZhNyHU4JMlbSDuesPjSSYTbE18CTnD3hrJ/oIiIDHmR5pJKGXsBS1auXElra2vZyr0liiLGjRtHU1PTkJ7/rOMwdI9B8cOLhg0bxrhx41i2bBnHH388ALOWv1zxUxybo4jTtt+uy9tntdNXz00o/m9hqD7gKc//JkaMGMHYsWMhzJb7fd62+1VOgoiIZIuiSA8sKkPHqOcpSBARGSRaunDlnN6mK9v3xLbSfylIEBEZJM7eLu9s8d7ZXgafaicuioiISD+lkQQRkQGspqaGxsbGLm/fG8l+Wu1w8FCQICIygPVEsp4ehCSlKEgQEZEuS49E9LX0o6Lb2tr6/fTH/tqvzihIEBGRLmtpaaG+vr7a3RgQ+mr9iJ6kxEURERHJpJEEERHpERe8vIKaKq2C2RJF7VM4q9mPYul+DUQKEkREpEfUxHGXl3UejP0YDBQkiIgMEkP12QWDWbX/TQdFkGBm5wLnpIpagVeAu4Bz3H11Um8qcD2wr7s/XqbNQt0sF7v7KUX1a5J9bglMcvc7Kuj/wcA3gPGEpzu+BTwOnOfuTxTVvSmpW+xZd98to09nAkcTHin9InCZu1+Zt28iMnCkkwgHYpKcdFTtf9PBlrg4EdiP8CjlO4ATgHvMrDuh1zFJm+nXpRn1vkgIEKD8o6GLfZsQHPwX4ZHR/w8YB/zWzA7MqL82o09fz6h3FTATuIxwbO4CLjezmRX2T0REhqBBMZKQssTd30w+32tmY4AphCv0h7vY5tJyow6JaUAz8BBwqJmNc/emnPs43t3fSBeY2S+BvwKnA78uqv++uz/aWYNmtgcwFZjp7hcnxfeb2VjgHDO7yt3fydk/EREZgqoeJCQnrQsJV9BbAe8CzxNuEyxK6kwETgX2AUYALwHz3f2iMs0/SggSdqDrQUJZZrYdcAjwc2AeMIEwxP+jPNsXBwhJ2Roz+yNhhKErvpS831BUfj1hdOSzgHexbRERGQKqHiQANwJ7Ee6dPw+MTr6PATCzacA1hKvpeuANYBdgt6zGiuyUvK/sRv+Gm9lGx8nd1xfVOYZw6+Y6YBGwAjiWnEFCFjPbHPg48MuMnzc1s9cJtzdeJdxa+b67v52qsxvQ5O7Ff/vTqd8VJIgMInEqo7+5ublP9pnej+YTdJQ+Jl35N9no+FZhxkZ/CBL2B+a6+zWpsjsBzGxT4BLCKMBB7l44QotLtFU4oW9KyEuoB14GHuxG/zoM65vZiEKgkOQ7HAMsBxa5e5uZzQPONLMD3P2BLu73SmAk8MOi8icISY1LgQj4DHAS8Bkz+1d3/1tSbwwhAXIj7v6Omb2f/N6BmdUC7ZkxDQ0No+rq6rr4J4hIX0ovjzxjxow+339rFDFSUw830pqajdDdf5OWlpY+f85GfwgSHgOmmtkqwlX4EndvTX4bD2wGXJkKEDrzWtH3h4Fvufu6rMrJCX54uixjlOCbwB87qTMBqAPOd/e2pOx64AzCaMIDFeyr0K+LgCOB6e7+ZNE2FxdVv9fMngJuJuRFXJb6rdQxizr57XRSM0XmzJnDrFmzSlQVEZHBrD8ECUcAZwHHARcAa83sdkJW/tikzoqcbR0MvEOYArnC3VeVqX8gcF+6wMx2dPdlqaI/lklcLMxkuN3MRiefVxFGIL5mZie6+7uEYOLeon1t5+4risouAL4HnObueZ//+j/AOuATbAgSVgH/XFzRzD5EuDXSYZQhcRFh9AaA6dOn70ES6IhI/5Z+RPPs2bP7ZLpcc3Nz+xXyCI0idJA+Jl35N0kf32o8grvqQUIyG+Ek4CQz2x74AuFe/lZsOFnlXdPyqdTshjyWAPsWlb2ad+Mkb2BS8vWJEtWOBK4Gfpuxr9eL2ruAEDCd7e4NefuR0pb6/AzwFTPbsuiY7J68L81qwN2bCbM0Ct7rQj9EpArSC+3U1tb2+Zx6Ld3UUfqYdPffpBqLY1U9SEhz9+WEefwTCLkKjxBGBurN7Oactxwq2d8awv39rppMuH9/BvCbot8i4FbCLYery+3LzM4jBAjnuvsPKuzHEYT8hXT+xB3AeYRZFulbFFOBvwG/qnAfIiIyxFQ1SEiGvu8Dfgo8B6whXG1PBG5z97Vm9l1gLrDIzK4hXH3vBOzh7t+pTs/bHQu8CVySXIFvxMxuBE40s13d/dlSjZjZacD3gQXAr8zsE6mfY3f/bVLvHwlTGm8B/kIIRD4NnEiYtXBdYSN3f9rMbgB+YGYxYdTk0KTP39MaCSIiUk61RxLWEYbhpxCS/0YQZgnMAhoA3P1aM3sVOI0QLETAMsJ6BFVjZnsTpij+Z1aAkLiacAKfBpzcSXOfT94PT15p77Ph3+kdQlByCrA1IbdgGfDfwEXu/veibY8nzO44Kan/InCCu1/R2d8mIiICEFVj3qUMKHsBS1auXElra2vZyr0liiLGjRtHU1NTVeYK9xc6DjoGBVnHoRoPA2pubm5/tsCs5S9X7emLzVHEadtvV/V+FEv3qyvPXsjzb5rnfxMjRoxg7NixAHsDv8+7/2qPJIiISA+JokgPdRpkqv1vqiBBRER6REsVH02d3nc1+1GsP/WlKxQkiIhIjzh7u7yz1XtXf+nHYDDYHhUtIiIiPUQjCSIi0mU1NTU0NuZdHLZnRVHENttsw2uvvUZbW1ufJ21WqhorJnaXggQREemyaibWRVHEyJEjqa2tJY7jPn/40VCgIEF6VTWmZImISM9QToL0qpaWFurr66mvr9/oMbYiItL/KUiQXtXc3Jz5WURE+j8FCSIiIpJJQYKIiIhkUpAgvaqtrS3zs4iI9H8KEqRXvffee5mfRUSk/9MUyD5mZnkeTXaeu5+b1L3C3b+TbFtHeNxzliXuvk/G/u4iPIp6truf1MVui4jIEKQgoe/tV6J8E2A+8BFgYZk2LgN+WlS2triSmY0DDk2+Tjazme6ueYgiIpKLgoQ+5u6PZpWb2aXAjsDx7v5YmWaWl2qnyNGEf+MFwOHAlwCvoLsiIjKEKSehHzCzKcAJwLXufnUPNn0s8CowFViXfBcREclFIwlVZmZ7AlcBvwO+nXOzYWZW/G/3vru35zuY2QHAzsBF7v6mmd0Rim17d1/eSX9qgfaF2BsaGkbV1dXl7JaIiAw+sipcAAAgAElEQVQmGkmoIjPbEridkE/wFXfPuyThLKC16DWhqM605P265P1awr/31DJtnw68U3jNmTPngZx9EhGRQUZBQpWY2XDgZmBb4Ah3f7mCzWcD+xa9fptqezPgq8AD7v6XpHgx8BJwjJl19pSli4APFV7Tp08/oIJ+iYjIIKLbDdXTQLj6P8Xd76tw2xXu/ngnvx8FfBBwMxudKnfgVOAgQtDQQTKakR7R0OIGIiJDlEYSqsDMjgJOBm5x94t7YReFWw2XA2+nXqcW/S4iIlKSgoQ+ZmYfA+YCS+mFk7WZ7U64/eDAZzJe9wOTzGzznt63iIgMLrrd0IeSE/MdhNkDs4DdzSyr6kp3f6GLuykEHrPc/fcl+vBp4OvAFV3ch4iIDAEKEvrWHoQFkwBu7KTePMrPQujAzGqAycDjWQFC4hdAEyGY6PUgYdSoUZmfRUSk/1OQ0Ifc/X6gs5kFxfWjou/LOts+WXJ5yzJtrgc+nLcP3TVs2LDMzyIi0v/p/7VFREQkk4IEERERyaQgQXpVbW1t5mcREen/lJMgvaqmpobGxsb2zyIiMnAoSJBeFUWRRhBERAYoBQnSp+I4pqWlpf0dwghDFG2YtFH8XUREqkNBgvSplpYW6uvrO63T2Nio0QcRkX5AiYsiIiKSSUGCVN2fdzuOv+x2XLW7ISIiRXS7QaouHjaCuNqdEBGRDhQkSK8qTlDMo7m5ub2+EhhFRKpHQYL0qnSiYmG9hHJmzJjRXl8JjCIi1VP1IMHMzgXOSRW1Aq8AdwHnuPvqpN5U4HpgX3d/vEybeUavz3P3c1P1r3D37yTf64AXS2y3xN33ydjnXcDngdnuflKO/ae3/S7wSWBPYAdgsbsfXKLu1sCPgcOADwJPAGe5+31F9VYQngb5paLy44ErgTuBo9y9uZK+iojI0NGfEhcnAvsBhwN3ACcA95hZV8ab9yvx+hTh5N8CLMzRzmUZbUwtrmRm44BDk6+Tk0c2V2I6sD2wCFhVqpKZfQC4DzgQOBH4YlL/V2b2yXI7MbPvAY3AfOBrChBERKQzVR9JSFni7m8mn+81szHAFGA88HAlDbn7o1nlZnYpsCNwvLs/lqOp5aXaKnI04VguIAQ5XwI8Z3cBdnH3tqSPnZ3s/wP4F+Bf3f13Sf37gWeAWcD+pTY0swbgVOC/gO+6u3IFRUSkU90KEsxsLHAh4Sp6K+Bd4HnCbYJFSZ2JhJPTPsAI4CVgvrtfVKb5RwlBwg5UGCSU6OsUwujEte5+dXfbK3Is8CphlOHl5HvuIKEQIOQwCXi2ECAk27aa2U+A881sa3d/Pb2BmQ0DrgKOA77v7hfk7VdPiOMNsUghIbFExU63FRGRvtfdkYQbgb2AMwnBwejk+xgAM5sGXAP8GqgH3gB2AXbL0fZOyfvKbvYRM9uTcKL8HfDtCjYdZmbFx+j99FW4mR0A7Axc5O5vmtkdodi2d/fl3e17kd2AezPKn07edwXSQUINcDPwVeDb7n5luR2YWS3Qni3Y0NAwqq6urqv9bZ/ZABsSErNE8frMbUeOHNnlfYuISPd0N0jYH5jr7tekyu4EMLNNgUsIowAHpU6si0u0NTw5IW9KGLKvJ1yVP9idDprZlsDtwFrgKxXeh5+VvNIOIeQOFExL3q9L3q8FjiSMKpxfaX/L2AJ4K6O8UDamqLyQJ3F+ngAhcTqpRNI5c+Ywa1bxIRARkaGgu0HCY8BUM1tFOHEucffW5LfxwGbAlTnvf79W9P1h4Fvuvi6rcpLQODxd5u7ri+oMJ1xJbwsc4u4v5+hH2mzgpqKyP6Xa34xwlf6Au/8lKV5MuKVyjJldUPjby41IVKCzbYp/ewIYC5xoZnenb1N04iJCcAfA9OnT9wAeqLiXifTaCLNnzwayRxTiqON/inq0tIhIdXU3SDgCOItwv/sCYK2Z3Q7MJJycAFbkbOtg4B3CFMgV7l4yyz9xICHTv52Z7ejuy1JFDcAE4JTiKYI5rSgz3fIowjREN7PRqXIn5GEcBCxOAoTWom2n0DEAKectOo4WQBhhKPyethz4GuE43Wtm/14uYTMZaUmPtrxXYR83kl4MqdM1DzIWTdJCSiIi1dWtICGZjXAScJKZbQ98AfgRIYmxcDW6bc7mnkrNbshjCbBvUdmrhQ9mdhRwMnCLu19cQbuVKNxquDx5Zf2+2N3Xm1lxX//ahf09A+yeUV4oe7b4B3d/wcw+zYZA4bM5Z2yIiMgQ12NTIJMkvcvNbAIhV+ERwshAvZnd3NNT7tx9DZB5lW9mHwPmAkvZcCLvUWa2OyFIcWBORpVzgElmtrm7v11uAaicbgcuNbO93X1J0o9NgG8ADxfPbChw97+mAoVfmdlEd/9ND/RHREQGsS4HCWb2IcJJ56fAc8AawklzInCbu69NVhKcCywys2sImfc7AXsUVjfsaWa2OWExplpC0uHuZpZVdaW7v9CNXRWCj1nu/vsS/fg08HXgijJ93pcw1RPg/wBbmdlXk++PpWZJzCUsvHSrmZ1OmPlxAvCPhDUUSnL3FzMChUc6/QtFRGRI686Ki+uA3xLurf8EuIeQmzCL5ITl7tcSlg8eTjjB3U24PdHTUwPT9iAsmDScMEXzNyVeZ3d1B8mKipMJyx53CBASvwCayDeSMQP4efL6MOH2QeH7AYVK7v53Qp7DA4TbG78g5H5MdPeHyu0kydf4NPAm8EszK7n4Uk+pqamhsbGRxsbG3ImIs2fPrqi+iIj0jkgL1kgZewFLVq5cSWtrce5l5Zqbm9sf+FTw/MemA7DL0+GuTdaDnaIoYty4cTQ1NQ3pRZZ0HHQMCnQcdAwK8hyHESNGMHbsWIC9gVIXtx30p2WZZYiK2lrRPAYRkf5HQYJU3c5L51a7CyIikqE/PQVSRERE+hGNJEifKiQyxnHc/lyHmpqajRZOUsKiiEj/oCBB+lQURe1JiXp4k4hI/6YgQfqdwihDerShtraWdevW0dzcTBzHHUYfRESk5ylIkH6npaWlwzTJYlnTJEVEpGcpcVFEREQyKUiQAeHPux3HX3Y7rtrdEBEZUnS7QQaEeNgIhu56aiIi1aEgQfqFdJJinuVVm5ubgY7TJ0VEpOfodoP0C4Vkxfr6+vZgoTMzZszIXVdERLqm6iMJZnYucE6qqBV4BbgLOMfdVyf1pgLXA/u6++Nl2swzMn2eu5+bqn9F4fHVZlYHvFhiuyXuvk/GPu8CPg/MdveTcuy/sN0mhL85y6nu/p9F9XcC/hP4DCHI+01S76mMtr8BzAT+CVgF3AKc7e7v5e2fiIgMXf1pJGEisB9wOHAHcAJwj5l1ZSx5vxKvTxFO/i3AwhztXJbRxtTiSmY2Djg0+To5eZR0pW7J2NdPivazNfAg8A9JP44ERgEPmNnORXWPBm4iBBGHAhcQHlv98y70TUREhqCqjySkLHH3N5PP95rZGGAKMB54uJKG3P3RrHIzuxTYETje3R/L0dTyUm0VOZpwLBcQgpwvAZ6zuwWv5djXacDmhNGUFQBm9jDwV+Bc4BtJ2SZAA7DQ3QsLDtxnZn8D5pnZIe5+b4X9ExGRIaZbQYKZjQUuJFypbgW8CzxPuE2wKKkzETgV2AcYAbwEzHf3i8o0/yghSNiBCoOEEn2dQhiduNbdr+5ue0WOBV4lXN2/nHyvNEjIYxKwqBAgALj7ajO7HTAzG+bubcD+hH+P64u2vxloTNrpV0FCOlmxkJSYUanT7UREpGd1dyThRmAv4ExCcDA6+T4GwMymAdcAvwbqgTeAXYDdcrS9U/K+spt9xMz2BK4Cfgd8u4JNhyVX5Wnvu3v7mcnMDgB2Bi5y9zfN7I5QbNu7+/IK9vVNMzuecAvoGeAyd5+X2s+mQB3ws4xtnyYEJnWEUYXdUuXt3L3FzP5EJ8ffzGqB9qUMGxoaRtXV1VXwZ3RNOgFx5syZmXWieH3mdnoGhIhI7+hukLA/MNfdr0mV3QntJ7VLCKMAB6VOrItLtDU8OSFvShiyrydclT/YnQ6a2ZbA7cBa4CvuXuIyNdOs5JV2CLAo9X1a8n5d8n4tIVdgKnB+jn3EhNyDhcByYBvgOOAGM6tz9/OSelsk729ltFEoG0MIEsaUqbtNJ/05nVQi6Zw5c5g1q/gQiIjIUNDdIOExYKqZrSKcOJe4eyFTfzywGXBl+sq7E68VfX8Y+Ja7r8uqnCQ0Dk+Xufv6ojrDCUPs2wKHuPvLOfqRNpuQ/Jf2p1T7mwFfBR5w978kxYsJt1SOMbMLCn97qREJd38fmFz02/+Y2T3AGWZ2ubuvSv3W2bEs/q1U3c7auIgQ3AEwffr0PYAHOqnfI9KPh25oaMgcTYijjv+56rHSIiK9p7tBwhHAWYQr3wuAtcn98ZnA2KTOihLbFjsYeIcwHXBF0Ykxy4HAfekCM9vR3ZelihqACcAp7r5R3ZxWlJlueRTwQcDNbHSq3Al5GAcBi0tMc5xCxwAk7SbCjI+9CPkD6dGCYsWjDKtSdYuP4xZkjzCEjoeRlvRoS59Ml0wviFTywU0ZiyZpISURkd7TrSAhmY1wEnCSmW0PfAH4ESFprnA1um3O5p5KzW7IYwmwb1HZq4UPZnYUcDJwi7tfXEG7lSjcarg8eWX9vtjd15tZcV//WqbtwtmvDcDd15rZi8DuGXV3J9xOWZZ8fyZV/nyhUjI185+AG8rsW0REpOemQCZJepeb2QRCrsIjhJGBejO7Oecth0r2twbIvMo3s48Bc4GlbDiR9ygz250QpDgwJ6PKOcAkM9vc3d8utwBUhsmE9RyeSJXdDnzbzD7i7q8k/fgQYbbCncnMBgjH/g1CXsSt6W4DHwBuq7AvIiIyBHU5SEhOTvcBPwWeA9YQTpoTgduSK9/vEk7Wi8zsGuB1wqyFPQqrG/Y0M9ucsBhTLSHpcHczy6q60t1f6MauCsHHLHf/fYl+fBr4OnBFJ/39HmF2xH2EWzNbE27fHAyc5e7pWwMNhLUQFiQrVbYAZxD+HQsJjiQjFzMJyY9XEgKZfyIcj3sK01NFREQ6052RhHXAbwn31usIayAsJ5yIGgDc/Voze5WwCNBcwhD6MmBex+Z6zB6EBZMgTNEsZR4ZqyfmkQzbTwYezwoQEr8AmgjBRMkgAfgjYTbH5wkLJf2NMHpwhLtvtNaCu79uZp8iLMs8n5C4+TBwoLv/uajuPDNbTzj20wi5CdcTckj6nZqaGhobG4F8ax/Mnj2b2tpaJS6KiPSiSIvRSBl7AUtWrlxJa2upR0z0rObmZurr6zcqe/5j0wHY5elwZ6exsbF0guMgFkUR48aNo6mpacguJKVjEOg46BgU5DkOI0aMYOzYsQB7A6UubjvoT8syi5QUtbWieQwiIn1LQYIMCDsvnVvtLoiIDDn96SmQIiIi0o9oJEH6nUISYxzH7c90qK2tZdy4cbz22mvEcayERRGRPqAgQfqdKIrakxILD2+KooiRI0dSW1s7pBOURET6koIEGfLiOKa5ubl91KKmpqZ9uef0ZxGRoUZBggx5LS0tTJ8+PfO3oTrVUkQElLgoIiIiJShIECnywkePrnYXRET6BQUJIkXaho2odhdERPoFBQkypBWSFosKN/pdRGSoUpAgQ1pLSwszZszYqCyK12/0u4jIUFX12Q3JI4/PSRW1Aq8AdwHnuPvqpN5UwlMM93X3x8u0mefy7zx3PzdV/4rC46vNrA54scR2S9x9n4x93kV4kuNsdz8px/4L2/0z8C3CY6X/AXif8GTIS9z9tqK6xwHXlGhqrLu/WVT/G8BMwmOiVwG3AGe7+3t5+zfQPPnkk9x0001MnjyZj3/84yXLemM/0jkdM5GBpz+NJEwE9iM8NvkO4ATgHjPryiT1/Uq8PkU4+bcAC3O0c1lGG1OLK5nZOODQ5Ovk5FHSeU1MXj8HvkZ4BPWLwK1mdkaJbb6Z0a/VRX06GrgJ+E3StwsIj4z+eQV9G1Cam5uZP38+q1atYv78+TQ3N2eW9cZ+pHM6ZiIDU9VHElKWpK6E7zWzMcAUYDzwcCUNufujWeVmdimwI3C8uz+Wo6nlpdoqcjThWC4gBDlfAjxnd39CGH1Ij37cY2ZjgdPN7MfuXvyM5mfc/clSDZrZJkADsNDdC89cvs/M/gbMM7ND3P3enP0bMBYsWMDq1SFWWr16NQsXLiSO4w5lkyZN6vH9dLfNwU7HTGRg6laQkJzILiRcqW4FvAs8T7hNsCipMxE4FdgHGAG8BMx394vKNP8oIUjYgQqDhBJ9nUIYnbjW3a/ubntFjgVeJYwyvJx8zxUkuPvKEj89BnwWGA2UqlPK/oR/j+uLym8GGoFJwKAKEl5//fX2oABCwuHdd9/d/rnwvmDBAsaPH8/WW2+90W/d2U9xm7IxHTORgau7txtuJFw1nw/8O3AcsAgYA2Bm0wjD+sOAesI9+0uBbXO0vVPyXukJsgMz2xO4Cvgd8O0KNh1mZpsUvTa6/WFmBwA7A/OSkZA7gEPMbPtudvszwOuEXIJivzSz983sLTO71cw+WvT7bsn70+lCd28B/pT6vQMzqzWzzQqvZcuWjerG39An4jjmpptu6nDCb2tro62trdO6lSQmltpPqXLRMRMZ6Lp7u2F/YK67p5Pp7gQws02BSwijAAelhtMXl2hreDJMvilhyL6ecFX+YHc6aGZbArcDa4GvuHslN0NnJa+0QwiBUMG05P265P1a4EjCqML5lfYXwMzqgU8C33b39FnuVeAHwG+BNcDHgNOA35rZfu6+NKk3Jnl/K6P5t4BtOtn96aQSSefMmcOsWcWHoH9pampi6dKl5SsSAoelS5fS1NTEhz/84R7ZT3faHOx0zEQGtu4GCY8BU81sFeHEuSR1/3w8sBlwZdH99lJeK/r+MPAtd1+XVTm5oh+eLnP39UV1hhOG2LcFDnH3l3P0I202Ifkv7U+p9jcDvgo84O5/SYoXE26pHGNmFxT+9iQASns/67iY2ecIoy03u/uVRX/fQjZOuPy1mS0kjBicB3ylqLlSx72zf4+LCMEdANOnT98DeKCT+lU3btw4dtttN/7whz90GDkoNmzYMD760Y8ybtw4gIoeOV1qP8VtygY6ZiIDW3eDhCOAswi3GS4A1prZ7YRpd2OTOitytnUw8A5hCuQKd88aZk87ELgvXWBmO7r7slRRAzABOMXdN6qb04oy0y2PAj4IuJmNTpU7IQ/jIGBxEiAUJx9OoSgAMbPDgP8B7iHMYCjL3V8ws0eAT6SKC8duDB1vV2xB9ghDob1mID3a0u+nS0ZRxOTJkznzzDM3Kh82LNxNS5+coihiypQp7U92rOQJj6X2U9ymbKBjJjKwdStISO7BnwSclNyD/wLwI0LSXOFqNE/+AcBTxfP8y1gC7FtU9mrhg5kdBZwM3OLuF1fQbiUKtxouT15Zvy929/VmVtzXv6a/JAHCbYSRiK9lzGjoTASkL6GfSd53JySSFvZRQ1gz4YYK2h4Qtt56aw477DDuvvtu4jgmiiI+97nPtScwFsoOP/xwttpqqx7dT3fbHOx0zEQGrh6bAunuy4HLzWwCIVfhEcLIQL2Z3ZzzlkMl+1sDZF7lm9nHgLnAUjacyHuUme1OCFIcmJNR5Rxgkplt7u5vdzYiYWaHEgKE+4EvJwmGefuxE2GdhPRtiEeANwh5EbemqwMfSPY16Bx++OE89NBDvP3224wePZrDDjsMILOsN/YjpemYiQxMXQ4SzOxDhOH+nwLPERLp9iUsDHSbu681s+8STtaLzOwaQrb+TsAehdUNe5qZbU6YYVBLSDrc3cyyqq509xe6satC8DHL3X9foh+fBr4OXNFJfw8knLRfIYzC7FnU32eTgAgzuw/4X8JIwbuExMWZwHpSyYbJyMVM4AYzu5IQyPwT4XjcU5ieOtjU1tbyzW9+s31Vv9raWoDMst7Yj5SmYyYyMHVnJGEdIct+ClBHWANhOeFE1ADg7tea2auEDPy5hGHxZcC8buy3nD0ICyZBmKJZyjwyVk/MIxm2nww8nhUgJH4BNBGCiZJBAmG2xEjCksxZeROfAh5KPj9DyIM4lTAi8AZhvYMLUomTALj7PDNbTzj20wi5CdcTckgGrY9//OMdlvzNKuuN/UjndMxEBp5I85SljL2AJStXrqS1tZI0iZ4VRRHjxo2jqampR+fWx3HMmjVrNnrI0593ncbOz14LhCmgI0eO7LH9dVdvHYeBRMcg0HHQMSjIcxxGjBjB2LFjAfYGSl3cdtCfnt0g0ueiKOo49J3KuFf2vYgMZQoSRIoMa6veiImISH+iIEGkyD/+oTdTZkREBg4FCSIiIpKpPz0qWqQqampqmDNnTvvDnmpqatpzESpZtllEZLBRkCBDXhRFjBw5sl/NYhAR6Q8UJIj0A3Ec09zc3GE0Iz2qISLS1xQkiPQDLS0tTJ8+vUN5Y2OjVicUkapR4qKIiIhkUpAgIiIimRQkiIiISCblJIhUURzHtLS00NzcnPl7oVwJjCJSDQoS+oCZ5XnyyHnufm5S94rCo7TNrA54MalzlLvfXNT2uYTHRI919zeTshuAo1PV2ghPpHwEON/dl3b9r5Ge1NLSQn19fcnfCw+eUgKjiFSDbjf0jf1KvD5FCABagIU52rnQzEbk3OffU/s5kPCI6L2AR8zsIxX1XkREhiSNJPQBd380q9zMLgV2BI5398fKNHMPcChQD1yWY7dtRft9yMyWA4uBw4Grc7QhIiJDmEYSqsTMpgAnANe6e54T9v8CvwLONrP/08XdvpO86zGHIiJSlkYSqsDM9gSuAn4HfLuCTU8DngBOBb6fYz+Ff99NgJ2AHwNvAwsq6a/0njjOk66Sv56ISE/SSEIfM7MtgduBtcBX3D07rT2Duz8F/BQ42cy2KVN9FGHEoJWQn/AM8M/A5939jU76V2tmmxVey5YtG5W3f1K5wjLMPVVPRKQnKUjoQ2Y2HLgZ2BY4wt1f7kIzZwEjCDMaOvN3YN/k9W/Al4HngYVmtl8n251OuC3xDvDOnDlzHuhCH0VEZBDQ7Ya+1QBMAE5x9/u60oC7LzOzK4HvmNklnVRtc/fH0wVm9ivgZeASwqyHLBclvwMwffr0PQAFCr0k76Oo9chqEakGBQl9xMyOAk4GbnH3i7vZ3A+AY4EfAs/m3cjd/2ZmLwB7dFKnGUjfAnmvq52U8vIukKSFlESkGnS7oQ+Y2ceAucBSYFp323P3VcAs4KvAv1bQj00JCYwlcxJEREQKNJLQy8xsc+AOoJZwYt/dzLKqrnT3Fypo+r8JMyMOLfH7MDP7ROEz8BHgRGBz4NwK9iMiIkOUgoTetwdhwSSAGzupNw+YmrfR5NbBuZReFOkDwG9S398A/ghMcvc78u5HRESGrkjzr6WMvYAlK1eupLW1emswRVHEuHHjaGpqGlRrBqQf8FR4TkPa7Nmzqa2tbX/A02A9DpXQMQh0HHQMCvIchxEjRjB27FiAvYHf521bIwkiVRRFUacPbqqtrdWDnUSkapS4KCIiIpkUJIiIiEgmBQkiIiKSSTkJIv1ATU0Nc+bMaX9GQyFRUSstikg1KUgQ6QeiKGLkyJGMHDmy2l0REWmnIEEkh8JUxazvNTU11NbWaulkERl0FCSI5NDS0kJ9fX3J3xsbGzVVUUQGHSUuioiISCYFCSIVeu3Lr/Hal1+rdjdERHqdbjeIVCjeZOgu/yoiQ4uCBJEcNloPPQbK5CgWJzYqqVFEBqKqBAnJ0wvPSRW1Aq8AdwHnuPvqpN5U4HpgX3d/vEybeS7vznP3c1P1r3D37yTf64AXk3pHufvNJfo81t3fTMpuAI5OVWsDmoBHgPPdfWmOPmFmJwEHAHsCdcCv3f3TJepuBTQAnwM+CDwFnOXui4vqLQOWuvvnisqPA64CFgDm7uvy9HGoS89siN6Pyo4mpBMdldQoIgNVtXMSJgL7AYcDdwAnAPeYWVcuu/Yr8foU4eTfAizM2daFZjYiZ92/p/Z1IHAW4cmJj5jZR3K2UQ/sAPwvsLJUJTOrBRYDE4AZwBeB14FfmtmB5XZiZqcC1wA/Ab6sAAGefPJJTjnlFJ588smNPvf2vkREBoJq325YUrgqB+41szHAFGA88HAlDbn7o1nlZnYpsCNwvLs/lqOpe4BDCSfuy3LUbyva90NmtpxwMj8cuDpHGx9197akv52NPkwDdgPGu/tvkvr3EUYTGoB/K7Whmf0QOJ3wN81w9yF/Y725uZn58+fz9ttvM2/ePABWr17N/Pnz+Zd/+ZcevfpP76s32hcR6Q0VBwlmNha4kHAi3Qp4F3iecJtgUVJnInAqsA8wAngJmO/uF5Vp/lFCkLADFQYJJfo6hTA6ca275zlZQ7iaHwacbWY3uPuaLuz6neS9NU/lQoCQwyTgT4UAIdl2vZndBPzQzD7i7q+kNzCzYcAVhKDnfHc/BwFgwYIFrF69GqD9vfB54cKFTJo0qdf21dPti4j0hq6MJNxIGE4/kxAcjE6+jwEws2mEIe1fE05MbwC7EK6Ay9kpeS855J6Xme1JuPf+O+DbFW5+GvAEIdD5fo59FY7jJoS/4cfA24T7/j1pN+DBjPKnk/ddCbkdBSMItxaOIIweXNrD/RmwXn/9dRYuXLhxQmIijmMWLFjA+PHj2XrrrXO119zcXPJ78b660r6ISDV0JUjYH5jr7tekyu4EMLNNgUsIowAHpYa0F5NteHKC3ZQwNF8PvEz2iTA3M9sSuB1YC3zF3ZvLbLIRd3/KzH4KnGxmV7p7Z5PiR9FxxKAJ+Ly7v1HJfnMYA7yVUf5W6ve0f0/ef5g3QEjyHtrHwRsaGkbV1dVV2M3+LY5jbrrppswAobjOySefnGtmwowZM0r+9rOf/azDviptX0SkGroSJDwGTDWzVcAiQl5B4SQ5HtgMuDLnPe/ik+/DwLdKJdQlCY3D02Xuvr6oznDgZk5lD2sAACAASURBVGBb4BB3fzlHP7KcBXyNMKNheif1/k6YmQDhNsVHCEmFC81sYip3oPhYv9/FvIDOtin+7UlgC+A7ZvaLUnkbRU4nNfNkzpw5zJo1q/Je9mNNTU0sXdr5xJO2tjaWLl1KU1MTH/7wh7u1v+eee65X2xcR6S1dCRKOIJxAjwMuANaa2e3ATGBsUmdFzrYOJty/bwVWuPuqMvUPBO5LF5jZju6+LFXUQMj+P8X9/2/vzuOkqK4Fjv+abfBhFEWQUaPgFh/BYGJ8URL3JSyJBo0HjWyKMUOiwaBxeS7gLhhNiAoYURFx4ZigGBYXXIOKGhKNiDwTI6JxUEJAIZFhm/fHqR6Knuqe6m26ezjfz2c+PV19697bd3q6Tt2695ZulTYbqrpURCZgB9hbMiTdnDo9U0SewHpEbsFmPUDj3oYzgSlZVmsljXsLwAIBaNzL8A/gZKzNnggHLRncgNUbgBEjRvQCXsiynmWturqanj17snjxYjZvjh4O0qpVK3r06EF1dXWsPMePH7/VQMS6urqG3oUePXqwZMmSrcrKNn/nnCuFrIOEYDbC+cD5IrIncCJwIzaIMXlw2SNmdm+EZjfEsRA4JGXbR8lfROR0YBQwXVVvziLfdK4FzgKuB96Ku5Oq/kdE3gV6hTan1vs9svcmcGDE9uS2RqfHqvqeiBzF1oHCS+kKCC7NhC/P/DuHepa1RCLBoEGDuOyyyzKmGTx4cOxLAVVVVWlnK5x22mlcddVVeeXvnHOlkNcUSFVdBtwmIsdiYxVewnoGakTkoUJPswtmGkQuqiQiXwEmYwfK4QUqb6WIjMVmc3SIu18wNmNfbNBmMq+Mi0HF9AgwQUS+oaqvBGW1AQYBr6jqR1E7Bb0iR2GBwuMi0ldV8549Usl23XVX+vXrx6xZsxqNF0gkEvTv358uXboUpKwuXbpsVVah83fOuWLJKkgQkR2xA80DwBJgDXaG3AeYoaprReQC7GA9T0TuxBb72RfolVzdsNBEZCdsMaYqYCxwoIhEJV2hqu9mmf2vsNkRfdO83kpEDk3+jo1J+CmwEzAmTgEi8nVspUWwMR0JEfl+8Pw1VX0/+P3uoC4Pi8glWBDyY+BL2KWbtFT1/ZRAoZ+q5jVAtNL179+f+fPns2rVKjp27AjY9MSOHTvSr1+/opZV6Pydc64Ysl1xcR3wCraWwf3YwkNnYwfmHwKo6l1AP2yA4WRgFnZ5YllhqhypF7ZgUmtsiubLaX6uyDZjVf0PmQ/224XyfxG4DdgEDMhiyuG5wMPBzxeBHqHnR4fqUoeNt3gWWxTp90A10FdVn4/xXpYBR2GB21wROSLzHi1bVVUVQ4YMoVOnTgwdOpShQ4fSqVMnhgwZUvCFjsJlFSN/55wrhkSmaWDOYWtgLFyxYgUbNsRaG6ooEokE1dXV1NbWZpy6WCzr1q1jxAib5FJ7ai0koFq3DDpMvT9DsW7wVOp2KAfeBsbbwdsgKU47tG3bls6dOwMcDPwpbt6lXpbZuYqw1UE+xvE+kUh4b4FzruJ5kOBclhIbfUaCc27b4EGCc1nqOqNrqavgnHPNotS3inbOOedcmfKeBOdiaNeuHZMmTWp4njowsV27dqWqmnPOFY0HCc7FEDUQsX379iWqjXPONQ8PEpyrUOHejFz3y3Z6ZiGnczrnyp8HCc5VqPXr11NTU9OsZaauB+Gca9l84KJzzjnnInlPgnMtwPKTl1PfJt6Kc4mNiYZpnHH2C6d3zm1bPEhwrgWob1MfO0goxH7OuW2DBwnONYNi3cuhlFrie3LObc3HJDjXDJKDDGtqanKakVCOWuJ7cs5trSQ9CSIyBhgd2rQB+AfwGDBaVVcH6YYB9wCHqOofm8gzTp/pVao6JpT+dlU9N3jeDXgvSHe6qj6Ups6dVfWfwbYpwNBQss1ALfAScLWqLopRp0x1v1RVb0xJuzfwC+AY7O/3MnCxqja6q5eInAZcAhwA/AtQ4HJVXRunXs61JK+//jrTpk1j0KBBHHTQQaWujnMVodQ9CX2Aw4D+wKPAecBcEcml3/KwND+HYwf/9cCcmHldJyJtY6b9PFTWkcDl2O2VXxKR3WPXHn5L47pPDScQkc7AH4D9gbMAAdoDz4nIl1LSngE8CLwG9AWuAoYBM7Kok3MtQl1dHVOnTmXlypVMnTqVurq6UlfJuYpQ6jEJC5Nn5cBTItIJGAz0Bl7MJiNVXRC1XUR+DXQHfqSqr8bIai52UK0Bbo2RfnNK2fNFZBnwNBb8/CZGHgAfp3sPIT8HOgO9VfV9ABGZD7wLXA0MDLa1Bm4CnlTVHwb7Pisia4D7RaSvqs6NWS/nKt7s2bNZvXo1AKtXr2bOnDkMGDCgxLVyrvxlHSQEZ7PXYQfSLsBnwDvYZYJ5QZo+2AHt60Bb4H1gqqre0ET2C7AgYS+yDBLS1HUw1jtxl6rGPVg/g/WwXCEiU1R1TQ5Ffxo8bshh30wGAM8kAwQAVf1MRGYAg0WkjapuBA4FqoFRKfs/DNwR5ONBQjOqr99yRSmfs9hEIsG6deuoq6vbOp9iTlAI5R0uM/x7+P2Vm48//pg5c+Y01LG+vp7Zs2fTu3dvdt111xLXzrnylktPwn1Yd/plWHDQMXjeCUBEhgN3As9jZ+OfYN3jPWPkvW/wuCKHem1FRL6KHRBfA36S5e4XA3/GAp0rY5SVbMc22Hu4CVgFzM6izB8EbdcKWATcpqr3hMrYDtgHeCRi378A2wF7Y3+TnqHtDVR1g4gsIcPfQkSqgIYl9caNG9ehW7duWbwNFyU8sG/kyJEFzz+xKUF92+IcqBObtlz9S1f39evXl+W9LOrr65k2bVqjICa5fdSoUT4rw7kMcgkSvglMVtU7Q9tmAojI9sAtWC/AMaqa/M98Ok1erYMD7PZY13wN8AF23T1nIrILdjBdC5yiqlmduqnqGyLyADBKRCao6vIMyTvQuMegFviuqn4Ss8gHsIDiA6x3Zjhwt4jsrapXBGl2AhLYAMRUyW2dUh7Tpe2WoS6XEhpUOnHiRMaOHRvjLThXfmpra1m0qPH44c2bN7No0SJqa2vZbbfdSlAz5ypDLkHCq8AwEVkJzMPGFSQPkr2BHYAJoQAhk9SD74vAOaq6LipxMKCxdXhb0L0eTtMaeAjYAzheVT+IUY8olwOnYgfMERnSfQ4cEfzeCtgdGAnMEZE+qvpyUK/Utt6UbCNVPSPltd+JyO+BS0Tk16oa7lnJ1K6pr6VLmymPG7BAD4ARI0b0Al7IkN7FEL6V9Pjx43O+/0EikaBr164sX76cdevWNZzZ17cuXnd/OO9w3evq6hrKL9dbZVdXV9OzZ08WL17M5s2bG7a3atWKHj16UF1dXcLaOVf+cgkSBmIH0LOBa4C1IvIIcBE2qA7gw5h5HYddv98AfKiqK5tIfyTwbHiDiHRX1aWhTeOAY4ELVXWrtNlQ1aUiMgE4V0RuyZB0c+r0TBF5AusVuAWbpQCNexvOBKZkyHca8B1sXMdc7PJFPVt6CcJ2Dh6TPQfJduwEfByRNqqHAYCg1yXc8/LvDHV0MYW7tKuqqvIKEtq3b09VVdXWXejF7DEP5Z2u7uXaZZ9IJBg0aBCXXXZZo+2DBw8u23o7Vy6yDhKC2QjnA+eLyJ7AicCNWDd58mC6R8zs3gjNbohjIXBIyraPkr+IyOnYYL3pqnpzFvmmcy021fB64K24O6nqf0TkXaBXaHNqvd8js+S31+Ygz89F5G/AgRFpD8R6NP4ePH8ztH1xMlHQm3EANjXSuW3CrrvuSr9+/Zg1axb19fUkEgn69+9Ply5dSl0158peXlMgVXUZcJuIHIuNVXgJ6xmoEZGHYl5yyKa8NUDkokoi8hVgMjbob3iBylspImOx2Rwd4u4XjM3YFxu0mcwr42JQEQZjvQ8LQ9sewYKzLyYvo4jIF4CTgcdCl15ewcZFDAOmh/b/Pjb+w9dKcNuU/v37M3/+fFatWkXHjh3p169fqavkXEXIKkgQkR2x7v4HgCXAGuwMuQ8wQ1XXisgF2MF6nojciXV37wv0Sq5uWGgishO2GFMVMBY4UESikq5Q1XezzP5X2OyIvmlebyUihyZ/x8Yk/BQbaDimqcxF5OdAD2xw54dsGbh4AjAmpaflF1jwMFtErsQuC1yCLajUUJaqbhKRi4D7ROQOrOdgP+xSzFOq+niT79q5FqSqqoohQ4Y0rLiY6+Ue57Y12fYkrMPOUgdjI+TbAsuwA/M4AFW9S0Q+wqYRTsa6zZcC9xakxtF6YQsmgU3RTOde7Ow6tuDSwRjSL4q0HbY0ctInwNvAAFV9NEYRS7BLNv2xwOJz4HUiloZW1RUicjgWLNzLlmWZj1LVJSlpp4nIJiyIGIaNQ5iKTV11zaxdu3ZMmjSp4feWoNLe00EHHeTLMTuXpUQ5L4LiysLXgIUrVqxgw4ZCrw0VXyKRoLq6mtra2rJeuKfYwu2wbt06ampqAKiV2ti3fE5sTFCt1bH3C6efNGlSyc/C/bNgvB28DZLitEPbtm3p3LkzwMFAo3v9pFPqZZmdcwWQ2Bh/lH44bZz9ssnbOdeyeJDgXAvQdUbXZt3PObdtKPVdIJ1zzjlXprwnwbkKFR44mI36+vqGe0m0a9cuqwWFKmGAonOucDxIcK5CJRKJnAcRluPNmJxz5ceDBOecixDucWmOMrLt1UkkEtv0iH7XPDxIcM65COvXr2+YYlquZs6cWeoquBbOBy4655xzLpL3JDjnXBOWn7w89mJV2UhsTDRMQ41bRngf54rNgwTnnGtCfZv6ogQJzV2Gc9nyIME51yLkMwjQFY//XSqbj0lwzrUIyYGGNTU1RZ+V4OLzv0tlK0lPQnBXxdGhTRuAfwCPAaNVdXWQbhhwD3CIqv6xiTzj9NNdpapjQulvT96+WkS6Ae8F6RrdgTFU587J2zeLyBRgaCjZZqAWeAm4WlUXNVUhEdkfOAc4GtgnyONt4Jeq+tuUtMOw9ohSrarLU9Kfht0F8gDsLpAKXK6qa5uql3POOVfqnoQ+wGHYbZIfBc4D5opILv1Rh6X5ORw7+K8H5sTM6zoRaRsz7eehso4ELsfunPiSiOweY/8TsPf/O+BU4Azgr8DDInJlmn3OpPH7XBlOICJnAA8CrwF9gauwW0bPiPm+nHPObeNKPSZhYfKsHHhKRDoBg4HewIvZZKSqC6K2i8ivge7Aj1T11RhZzcUOqjXArTHSb04pe76ILAOexg7+v2li/4ewHo1wT8hcEdkFuFhExqpqXco+izL1rIhIa+Am4ElV/WGw+VkRWQPcLyJ9VXVujPfmnHNuG5Z1kCAinYHrsANpF+Az4B3sMsG8IE0f4OfA14G2wPvAVFW9oYnsF2BBwl5kGSSkqetgrHfiLlVt6mCd9AzWw3KFiExR1TU5FP1p8LihqYShICnVq1iQsTN2CSMbhwLVwKiU7Q8DdwADsGDIuRYjvPpgXV1qXJ29rfIop0kHobqsW7eOurq6sl55MdyO5VxPFy2XnoT7sO70y7DgoGPwvBOAiAwH7gSex87GPwH2B3rGyHvf4HFFDvXaioh8FTsgvgb8JMvdLwb+jAU66br8w2Ul27EN9h5uAlYBs7MsN+xorB0+iXhtVhCsfQo8B1yZMv4h2dZ/Ce+kqhtEZAkZ/hYiUgU03BBg3LhxHbp165ZL/Z1rVuFBcSNHjixo3olNCerblscBLrFpy9XYgQMHlrAm2Vu/fr3fN6TC5BIkfBOYrKp3hrbNBBCR7YFbsF6AY0Jd6E+nyat1cIDdHjtrrgE+AP6QQ70aBF31jwBrgVMiuuszUtU3ROQBYJSITEgdEJiiA417DGqB76pq1AG+SSJyNnAUMFJVN4VeWo714izAenAOxAYmLhCRb6rqG0G6TsHjvyKy/xfQLUPxlxIaVDpx4kTGjh2bw7twzjlX6XIJEl4FhonISmAeNq4geZDsDewATEi5xp5O6sH3ReAcVV0XlTgY0Ng6vE1VN6akaY1d598DOF5VP4hRjyiXYwMJRwMjMqT7HDgi+L0VsDswEpgjIn1U9eWgXqltvSmqjUSkL3A78FtSxkSo6uPA46FNL4jIbOBN4GrgpJTs0v0NMv1tbsACPQBGjBjRC3ghQ3rnykL4Ntbjx4/P+Q6ZSXV1dQ09EvWty6MXAbauy/Tp01m9enVZd+OH29FvNV55cgkSBmIH0LOBa4C1IvIIcBHQOUjzYcy8jsO6zDcAH6rqyibSHwk8G94gIt1VdWlo0zjgWOBCVd0qbTZUdamITADOFZFbMiTdnDqIUESewHpEbsFmHkDj3oYzgSkp+30bm33wFHBGnEArqOd8bBxCUrIdOwEfp+yyM9E9DMn86oBwz8u/m6qDc+UgvEhPVVVV3kHC1pkXLqu8herSvn17qqqqyjpICPOFlCpP1kFCMNDufOB8EdkTOBG4ERvEmDyY7hEzuzcyDNyLshA4JGXbR8lfROR0bLDedFW9OYt807kWOAu4Hngr7k6q+h8ReRfoFdqcWu/3wk+CAOFRbCzHKaqazaojCWx9haQ3g8cDgcWhMtpgayY8mEXezjnntlF5TYFU1WXAbSJyLDZW4SWsZ6BGRB6Keckhm/LWAJFT/0TkK8BkYBEwvEDlrRSRsdg4gA5x9wvGZuxLaNBhE1MWT8AChPnA97IZQyEi3bG2nxfa/Ao2LmIYMD20/fvY+A9fK8E551yTsgoSRGRHrLv/AWAJsAY7Q+4DzFDVtSJyAXawnicid2Ld3fsCvZKrGxaaiOyEHWSrgLHAgSISlXSFqr6bZfa/wmZH9E3zeisRSXb1J8ck/BTYCRjTVOYi8i2s7suxHouDUuq+WFU/C9LOw8YH/IUtAxcvwsYYXJHcQVU3ichFwH0icgfWc7AfdinmqWBsg3POOZdRtj0J67Cz1MHYCPm2wDLswDwOQFXvEpGPsGmEk7Gu8KXAvQWpcbRe2IJJYFM007kXO7uOLbh0MIb0iyJtB7wcev4JtqzyAFV9NEYRxwV5dMPWaEh1NDbNEewywkDgwmCfT4J9rlHVd1LqPU1ENmGzH4Zh4xCmYlNXnWtx2rVrx6RJkxp+d+XB/y6VLVEpA15cyXwNWLhixQo2bGhybaiiSSQSVFdXU1tbWzGDtIrB26H52qCuro6amhoAaqW2KLdxTmxMUK3VWZUR3mfmzJmsWrXKPwvb8P8DxGuHtm3b0rlzZ4CDgT/FzbvUyzI751zZS2wszqj8cL5xyyhWXZyL4kGCc841oeuMri2iDOeyVeq7QDrnnHOuTHlPgnPORQgPuCuW+vr6hntOtGvXLqvFhhKJRGEXjHIuggcJrintAdq0KY+PSrnUo9S8HZqnDcp9NH4ikfDPAv7/kJSpHUKvZXWHLZ/d4JryA+D+UlfCOedcQZyBrXUUiwcJrimdgG9ja11E3nirOSxdurTDxIkTXxgxYsQR3bp122bvJ+Ht4G2Q5O3gbZAUsx3aY+vxPMGW+/s0yYMEVxFEZAdsye8dkytQbou8HbwNkrwdvA2SitkOPrvBOeecc5E8SHDOOedcJA8SXKWoA64KHrdl3g7eBkneDt4GSUVrBx+T4JxzzrlI3pPgnHPOuUgeJDjnnHMukgcJzjnnnIvka1m6khGRY4BBQG/gi8Bq4I/A1aq6sIl9TwZOBQ4Bdgc+Bl4ExqjqX4tZ70LLsx2OAy4GvowtfPUpsAj4harOKWa9CymfNojI61rgMuAtVe1Z6LoWU56fhWHAPWlerlbV5QWsatEU4rMgIicBo4CvAq2xxeDGq+pvilHnYsjzs/AccGSGJLE/D96T4EppBLYC2HigHzAS6AIsCP5BMrkY+C/gOqAPcDn2hfAnEflysSpcJPm0QyfgLeBnwAnAj4ANwGwRGVSsChdBPm3QQEQOAi7EgsZKVIh2OBM4LOUn9gp7ZSCvNhCRS4AZWLAswInABKC8b8TRWD7t8GMafwaOxb4bFmQTMHpPgiuln6jqJ+ENIvI48Dfgf4FnMuz73Yh9n8HOGH4GnF3YqhZVzu2gqtOB6Sn7zgLeA84BphW8tsWRz2chmb4NdiZ9B9AL2KUI9Sy2vNsBWKSqfyxG5ZpJzm0gIgdjJw6Xquq40EtPF6OiRZbP98Li1G0iMhRoC0zOphLek+BKJvUfINi2FliMda9lu+9HwIdN7Vtu8mmHNPltwLomN+Zfu+ZRoDa4BNgZu9RQkQr9WahEebbBudhaAbcWoWrNqgifheHAWlJOKpriPQmurIjIjsDXiHfGlLrv3sBewKOFrldzy7YdRKQVFvR3wS457I9dkqlY2bSBiPTALjmdrKprRaTY1Ws2OfxPzBKRztj4lOeAK1V1UZGq1yyyaIMjgLeBU0TkCmBfoBbrUbtSVdcXtaJFluv3o4jsBxwOTA4Cjdi8J8GVm9uBDliXYWxBV/NdWKT8yyLUq7ll2w5zsOuN/wDOBwaq6uwi1a25xGqDIEC6G5hRSYM1sxD3s7A8SHM2cDRwBTawd4GI9CpqDYsvbhvsDuwH/Dr4OQ6Ygo1TSTeos5Lk9P2I9SKAfUdmxXsSXNkQkWuwe52fl82IdhFJYB/+w4FTVPWDIlWxWeTYDucBHYFqbET0dBEZqqoPFqmaRZVlG4zCDgwnFr1izSybdlDVx4HHQ5teEJHZwJvA1cBJRatoEWX5WWgFfAE4XVUfCrY9KyIdgPNFZLSq/q2I1S2aPL4f2wBDsdk+C7It14MEVxZEZDTWXXyZqt6WxX4JbCDOIGCoqs4sUhWbRa7tkDLt8zERmQvcLiLTVXVzoetZTNm0gYjsiR0ALwHWi0jH4KU2QKvgeZ2qfl7MOhdDrp+FMFVdKiLzgUMLWrlmkkMbrAS6Ak+kbJ+L9bB9DRv4V1Hy/Cz0w9pkbC5l++UGV3LBP8AYbI2D67PYLxkgnAmcraqVMpI/Uq7tkMarwE5A53zr1ZxyaIO9ge2waWKrQj/fBP47+P2GolS2iAr8WUgAFRUoQs5t8Jc02xPB47bSDmHDgfXAfbmU7z0JrqSCwUVjgGtV9aos9ksAd2IBwo9UtaKvN+baDmnySmALqaymgubH59gGr2PX31P9CtgR+3x8WJAKNpMCfxa6YwHTvAJUrdnk0Qa/w9YL6Qs8ENreDwsQXitUHZtDvp8FEemKvfcZqprTd4HfBdKVjIhcAPwCu47a6B8gef1MRO7Crqnto6rvB9tuxaY73Y0FC2F1qvrnIla9oPJsh5nAG9jBciWwGzAM+6L8iapOaIa3kLd82iBNfs8Bu1Tgiov5fBbmAS9gZ9OfAQcCF2HX6HtXygyHPNugLbAA2Ae4EpsueBw2cHGiqp7XHO+hEArxPyEiFwM3Aieo6lO51MN7ElwpfTd47BP8pEp2EbYOfhKh15L7nhX8hL2PrVRWKfJphxeB72MB0w5sWbr1OxU2uyGfNmhJ8mmHN4GB2AFxO+ATbKrcNar6TlFqWxw5t4GqbhCR44HrsQWHdsYWFrsEuKVYFS6SQvxPnIUtMJdzT5L3JDjnnHMukg9cdM4551wkDxKcc845F8mDBOecc85F8iDBOeecc5E8SHDOOedcJA8SnHPOORfJgwTnnHPORfIgwTnnnHORfMVF5yqUiNQDqGpLXX2w6ERkGHAPcJWqjiltbZwrPx4kOOfKgoiMAUYDZ6rqlALlOQVb1/5oVX2uEHk6ty3xIME5ty17BLsh0D9LXRHnypEHCc65bZaqfgp8Wup6OFeu/AZPzlWoqDEJItINu+vd80B/4BrgVGAX4G1gtKr+Pkh7KvBz4MvAWmA6cLGqfp5SzlJgL2yg80+BHwF7Y7em/h1wpaquTtnnOeBIoLuqLk15raGOqnpUShlRjlbV50SkPXAGcCJ2G+RqoA67NfIEVX0oqn3S6K6qSzONSRCR/wIuwO6suA+wHrstd6Oywu9BVRMicjbWVvthQchMrG1Xp+6XSdBWlwLfDt7vZ9jf9mpV/UtK2ob3AjwAXA0cDXQGTlbVR8N/F6B3UMcvA5tUtWMor8OwOyf2xu4uWgvMAa5V1Y+yLTeb9+zKi89ucK5lagc8DQwGXse61HsBj4jIcSLyM+wLfSPwJHar2fOAyRnyvBW4CfgQO+gl93leRL6QZ31/ix2AwW5/fW/oZ3mwvVtQv28Ay4I6vA4cCjwYjGkIuxd4N/j9iZQ812aqTPB+XsAOeF2AWUG9/ico61cZ9h0H3I4d0B/HbuF7DvCYiMQeZCoi38La5Jygvo8BfwVOBhaIyNFpdv0S8FpQ12eBp4ANKWkuBe7DAp9ZwKJQuYOAP2C3Kv4/YAYWjI0A/iQiB+RRrqswfrnBuZbpMOA5YH9VXQVbnfFNBHYGjlHVPwSv7Qb8GfiBiFyhqn+PyHMwcJiqLgz22R47UB+DnUWOyrWyqnphcJDvBUxOM3BxBXZGPU9VNyc3ikh34BngChGZkuy5UNVhwcDFfYAbsxy4eD1wMDAPGKCqa4OyDsDO5EeKyJOqOidi30HAN1T19WCfXYCXgcOxM+xnmipcRHYAHga2A05V1d+GXjsOmA3cJyJ7q+r6lN1PA24DzlfVTWmKGIL9/Z9PKfeLwG+AeuBEVZ0VbG8F3AycD0zFAoFUccp1FcZ7EpxrmTYBP0wGCIGp2IF2X+C2ZIAAEHQh3x88PSJNnrclA4Rgn7XAudgBZbiIVBWw/o2o6kpVfTIcIATb3wOuw77PvptvOSLSARgObAZ+nAwQgrKWANcGT3+aJosrkgFCsM8/scAM0rdtqrOArsAvwgFCkN88YAKwO/CdiH1XYJc2Mh2o70oNEAJnY4HJg8kAIShzM3b54SPgEBE5NMdyXYXxIMG5lmmpqv4tvCH4on8/ePpUxD7JrvnqNHk2ug6vqm9jXeI7AF/JrarZEZFvicjlIjJRRO4JegtODV7erwBFHIwdKF9V1b9GvH5f8PjNNJcPhMHGVAAAA91JREFUnozY9k7wmK5tUx0fPKa7nj8/eDwk4rV5qvqfJvJ/LM32w4PH+1NfUNU6rHcjnC7bcl2F8csNzrVM/0iz/d8ZXk++lq5H4P0025cCBwG7xapZjkRkR+z6+DEZkuU7NgK2vI+lUS+q6moR+RTYEQuOUmdHfBixW7I3Im5vS7fg8RURyZRul4hty2Lkny5Nxvce2h71t45TrqswHiQ41zI1NW2pkNOasl3xMdcezLFYgPACcCU22G61qm4SkROwwYmFXH0yThs1SqOqhWjb1sHjw0Cms/NXIrati5F/U2ly+fzEKddVGA8SnHNx7QW8GbF9z+AxPDUuOZhu+4j0X8yx/AHYWIsTg/UNwvbOMc8oyffRPerFoEdjR6znZU0Byw37EJstcG3qVMci+ygotztbLpGEJaep1jZbjVxJ+ZgE51xcA1M3BKP9D8IOluGDWfIgsn9EPiekyT8ZWKQ7edkJWBMRIACk65NvKs8oC4HPgf8RkagxDoOCx/kF6jWIMi94/F6R8k8nOZj1jNQXRKQdW8Z+/CH1ddcyeZDgnIvrXBH5avJJMAvgVqyL/+5gYFtScuT8BcGiRMl9jsOm0UVJnsF/Kc3r7wAdRWSrYCVY8yHdmgFN5dmIqv4buBv7frw9eJ/JsvYHLg+e3ho3zxzcgc0W+F8ROTN1gKSIdBCRISKyR4HLvQsLkE4Xkf6h8lph00J3B15T1QUFLteVKb/c4JyLaxo2kO4ZbLDeEdg0vbewGzOFPQhchK3Y97aIvIZdZvg6cAtwYUT+T2LXtX8mIj2xA3w9cJOq/h9wQ1CHh0TkJ1iXfC/gAOCXwM8i8vw9Nn7hZhE5ni33aLhYVVdmeK+XYos0HQ/8XUSeBzpgYyLaA79W1dkZ9s+Lqq4SkQHYLIS7gdEisghb1GhP4L+D+nyV6IGSuZa7TETOAaYAvxeRF4EPgK9hgdbH2BoLbhvhPQnOubjOww6eewEnYQfw24HDUy8BBEs7H4sFC18A+mHfNwODfRoJ1mo4CVsd8lvYWgHDCaYNqur92FLTC7BLHH2xQOIY0kzpC9Z1GIQFMicE+Q2niVkQqroGW754NBZYnIhN+/sj8ANVHZlp/0JQ1Rex5advxs7uj8Heww7YKokDgcVFKHcaFgDOwoKR72NTQicCBwdrRbhthN+7wTmXUfieBKWui3OueXlPgnPOOecieZDgnHPOuUgeJDjnnHMuko9JcM4551wk70lwzjnnXCQPEpxzzjkXyYME55xzzkXyIME555xzkTxIcM4551wkDxKcc845F8mDBOecc85F8iDBOeecc5H+H+s3kAawguSZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "d = {'model': model, 'error': error}\n",
    "df = pd.DataFrame(data=d)\n",
    "fig = plt.figure(dpi=100)\n",
    "grid = plt.GridSpec(1, 10, wspace=0.2, hspace=0.2)\n",
    "ax = plt.subplot(grid[:, 2:10])\n",
    "sns.boxplot(y=\"model\", x=\"error\", data=df)\n",
    "ax.artists[0].set_facecolor('tab:red')\n",
    "ax.artists[1].set_facecolor('tab:red')\n",
    "ax.artists[2].set_facecolor('tab:red')\n",
    "ax.artists[3].set_facecolor('tab:red')\n",
    "ax.artists[4].set_facecolor('tab:blue')\n",
    "ax.artists[5].set_facecolor('tab:blue')\n",
    "ax.artists[6].set_facecolor('tab:blue')\n",
    "ax.artists[7].set_facecolor('tab:blue')\n",
    "ax.artists[8].set_facecolor('tab:green')\n",
    "ax.artists[9].set_facecolor('tab:green')\n",
    "ax.artists[10].set_facecolor('tab:green')\n",
    "ax.artists[11].set_facecolor('tab:green')\n",
    "plt.xlabel('Imputation error', fontsize=15)\n",
    "plt.ylabel('', fontsize=15)\n",
    "plt.xticks(fontsize=12)\n",
    "plt.yticks(fontsize=12)\n",
    "plt.title(\"Cortex\", fontsize=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Impact of the linearity on imputation error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# zifa-ln\n",
    "zifa_ln = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/scPI_ZIFA_linear/500-impu/\"\n",
    "for t in range(n_repeat):\n",
    "    zifa_ln[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_ZIFAVAE_linear_imperr.txt\")\n",
    "\n",
    "# zinb-ln\n",
    "zinb_ln = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/scPI_ZINB_linear/500-impu/\"\n",
    "for t in range(n_repeat):\n",
    "    zinb_ln[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_ZINBVAE_linear_imperr.txt\")\n",
    "\n",
    "# fa-ln\n",
    "fa_ln = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/scPI_FA_linear/500-impu/\"\n",
    "for t in range(n_repeat):\n",
    "    fa_ln[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_FAVAE_linear_imperr.txt\")\n",
    "    \n",
    "# scvi\n",
    "scvi = np.zeros([n_repeat,])\n",
    "data_path = \"/Users/jiazhao/Downloads/Results/Cortex/scVI/imputation/\"\n",
    "for t in range(n_repeat):\n",
    "    scvi[t]  = np.loadtxt(data_path + str(t) + \"impu_CORTEX_scVI_imperr.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Cortex')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhgAAAGPCAYAAAAeD6l2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XuYXFWZ9/3vTuiuxmQgMSTQgNC8g4wXqIkcZgRHcADHJKgDHm7kkDEQR0mDTIKAolECiiFxIKKkM3IQMIJ4zwiRkeT1IXhAQESiMJIBeYYhhJAGkgwEgqa6Q+/nj7Wrs1Opqq7qru6q6v59rquurlp77bXW7jTsu9ZpR3EcIyIiIlJNo2rdABERERl+FGCIiIhI1SnAEBERkapTgCEiIiJVpwBDREREqk4BhoiIiFSdAgwRERGpOgUYIiIiUnUKMERERKTqFGCIiIhI1e1W6waIiFQiiqLRwEeBDwLvBiYBbwJeAZ4CfgXcGsfx4zVo20ygDfhFHMe/GOr6RepJpGeRiEijiKLo3cAtwCGp5G7gNWAcO/fK3gGcFsdx1xC27xfAccBlcRzPH6p6ReqRhkhEpCFEUfQh4BeE4GIzcAlwSBzHzXEcTwCagaOAK4FXgY8QejZEpAY0RCIidS+KorcC3wcywH8BH4jjeH06TxzHbwCPAI9EUfQN4LtD3lAR6aUeDBFpBF8D9gC2AafkBxf54jj+3ziOTwa2pNOjKNoniqJvRFG0JoqirVEUvZ68XxRF0d6FyoqiqC2Kojh5tUVR9JdRFF0XRdEzURRloyhaG0XRzCiKYsLwCMClqXN6zy1Q9ruiKPpuFEVPR1H0p6RNj0VR9LUoivYqkP89URRtT8qbU6S9+0dRtDnJc12p35PIYNIcDBGpa8mNfwPhC9GNcRx/qp/lHAcsJ8zVAPgTEANjks8vAx+O4/j+vPPagGeSj2cA3wHGps7fBHweuAZ4M9AEvA5szWvCUXEcP5cq9zLgy0CUas9uhKEegE7gpDiOf5/Xni8DlwNdwLvTx6MoGgX8jBDoPAkcEcfxn0r+YkQGiXowRKTe/R07/l91Z38KiKLoLewILv4L+Ns4jsfEcTwWOBb4IzAe+HEURfuVKOo7wBpCsJA7/+/jOP5hHMf7AA8m+f4ljuN98l7p4GIO8BVCEHIJ0BrH8RjCnJEjCUFCK3BXFEVj89pwBWEuSjNwexRFY1LH5hGCiyzwCQUXUksKMESk3h2Wev/7orlK+yIhuHgZOCGO4wdyB+I4/hVwImFi6JsJN/xiNgMnxnH8SOr8pyppSDL0cQWh9+OUOI6vjOP4haSsN+I4Xg18AFgN7A/s1GMTx3EPcGbSlkOAa5Ny30MIWgAujuP4sUraJVJtCjBEpN5NSL3/30pPjqIoAiz5+K+5m3laMqfjX5OPnyhR3LVxHOcPfVTqDEJPxSNxHN9bKEMcx9uBHyQfP1Dg+PPA2cnHmVEUtQO3AaOBu+M4/tYA2ygyYFpFIiL1Luo7S0kHEXomAFaVyHcPcDEwIYqig+I4fqZAngcKpFXqb5Ofb4+iaJdgJ2X35OeBhQ7GcXxXFEXXAucBS5LkTmBmFdooMmAKMESk3m1KvX8zYcJnJSal3j9fIl96ZcokdkzsTHupwroL2Tf5uTs7gohSSu3lcSFwCpCbN3J2HMebSuQXGTIaIhGRercm9f5dAyyr3GVzxfK9McD6IQxjQBiuicp4tZUo6yR2BBewY5msSM0pwBCRevdzoCd5f0o/zk/3OrylRL79U+839qOecuWGRd4xkEKSlTE3JB//M/l5cRRFxw+kXJFqUYAhInUtjuMXgR8lH0+PouiQUvnTkgmez7BjcugJJbKfmPzcXGT+RTlygVCpeSO5eRzvjqKo4PyKviQPfLuVsLT2vwgPfbuT8P/0ZVEUTShxusiQUIAhIo1gHmHPiN2BO/rYq4IoisZHUfQjYM847Cb4w+TQZ6Io2qdA/n2BzyQff5B/vAKvJj/HlcizDPgzYahkSRIsFBRF0agoigqVNQ94L2G/i9PiOP4zYTnresIcj5v60XaRqlKAISJ1L9lrYgZh98rDgEejKPp8FEUH5/JEUTQ62Xr7cuB/CA87y/k64XHubwZWRVF0TOq89xBWl4wj9HRcOYCm5h4RP71YEJQsk/1C8vEk4J5kC/DRSXuiKIreFkXRBUl5H0yfn7T3y8nHi+I4/s+k3P8l7I/RA3woiqLzBnAdIgOmrcJFpGEkN9ebgYNTyV2E3o3049pj4Hbgk3EcdyfnHgf8GNgzyfN68jO3E+YrhK3Cf5VXZxs7VpQcFMfx2hLteythPkQL4Ua/kfD8FAi7h65P5b0IWMCOSZ9dhMfO70HYbjznzDiOb03OGQc8BhxA2O9ip+AjyXM5IQDZBvx1HMd/KNZekcGkAENEGkryTf/jhG/2f0NYUro74cFmTwK/BJbFcfzHAufuQ1jaOR1oIwQizwJ3A1cV2oSrkgAjyf9uwm6gf0PYJCy3HcAu5yY9MO2EuSFthGDnVeC/CZNblwMPJcM8RFH0b8DHCBNF3xnH8S6TUZPfz33AMYT5GUcmQygiQ0oBhoiIiFSd5mCIiIhI1SnAEBERkapTgCEiIiJVpwBDREREqk4BhoiIiFSdAgwRERGpOj2uXQbTBOADwFp2bDYkIiKNpYWwT8tPgc3lnqQAQwbTBwgPZBIRkcZ3BnBbuZkVYMhgWgvw8ssvs3379kGrZOLEiWzcOJhP164vtbzeOI7p7u4GoKmpifCw0sGlf9/hTddb/3bbbTfGjx8Pyf/Tyz5vUFojEmwD2L59e+9NqdpyN7jt27czEnalrfX1ZrNZ2tvbAejo6CCTyQxqfbW+3qGm6x3ehsH1VjTUrUmeIlK2bDZb8L2ISD4FGCIiIlJ1CjBERESk6hRgiIiISNUpwBCRsvX09BR8LyKSTwGGiJTt9ddfL/heRCSflqmOIGY2E7ipyOGr3P3CVN5m4HlgL+AUd18++C0UEZHhQgHGyHQW8GRe2oa8z/9ACC4AZgEKMEREpGwKMEamx939kT7yzAKywP3ANDNrdffOwW+aiIgMB5qDIbsws7cA7yf0WlwFjAY+WdNGiYhIQ1EPxsg02sx2+rd39/TDQs4iBJ/fBVYB64GzgSuHrIUiItLQ1IMxMj0EdKdfuYDDzCJCgLEOWOXuPcAtwFvN7NhShZpZxsz2yL3Wrl07ZlCvQkRE6pYCjJHpH4Gj0q9UD8YJQBtwcxJcQFh5EhN6MUq5BNiSey1duvS+KrdbREQahIZIRqYnSkzynJX8vNPMxiXvNxN6PT5uZue7+6tFzl0AXJ37MHv27MmAggwRkRFIAYb0MrPxwCnJx98XyfYJ4LpCB9w9S1h5kqOdmERERigFGJJ2JpABvgj8Ou9YBPyIMExSMMAQERHJUYAhaWcDm4Crk96InZjZMuB8MzvM3dcMeetERKRhaJKnAGBmRwBTCJM7dwkuErmei1lFjouIiADqwRhR3P1m4OYix1YThkFKnb+mrzwiIiKgHgwRqcCYMWMKvhcRyacAQ0TKNmrUqILvRUTy6f8QIiIiUnUKMERERKTqFGCIiIhI1SnAEJGyZTKZgu9FRPJpmaqIlK25uZmOjo7e9yIixSjAEJGyRVGkngsRKYuGSERERKTq1IMhIn2K45iurq6yjjc3NxNFxTd87eu4iAwPCjBEpE9dXV20t7dXpayOjg4Ns4iMABoiERERkapTD4aIVOSyZ9fRHMc7pXVFEZceeEBZx0VkZFCAISIVaY5jMnkBRCXHRWRk0BCJiIiIVJ16MESkqNzqkGw2W9M2ZLNZ4jjWChSRBqIAQ0SKqubqkf7KZrPMnj0b0AoUkUaiIRIRERGpOvVg1Ckzmw9cmkrqBp4H7gIudfdXknwzgZuAo9z9kT7KzOUt5Cp3vzAvf3NS517AKe6+vOILERGREUkBRv2bCmwB/gKYDvwz8Ndmdoy793eq/lnAk3lpGwrk+wdCcAEwC1CAISIiZVGAUf9Wu/um5P09ZjYBmAEcAzzQzzIf76u3IzELyAL3A9PMrNXdO/tZp4iIjCAKMAaBmU0ErgCmAZOAV4GnCEMbq5I8U4GLgCOBJuBZ4HvuvqCP4h8iBBgH0v8Ao09m9hbg/cC/AbcAJwCfBK4crDql/sQF9rPoT7dZ+pxKVqREUcS2bdtKtkdE6pMCjMGxDDgc+BIhsBiXfJ4AYGazgOuBXwLnAC8BhwBvL6Psg5OfGwfQvtFmttO/vbtvz8tzFmES8HeBVcB64GxKBBhmlgF6p/gvWrRoTFtb2wCaKbVW6AFn3VFES4U3+u7U0tK5c+cOqD0tLS39Pl9Eho4CjMHxHuAGd78+lfZjADMbC1xN6H04PjWP4t4iZeWCgbHASYSA5DngVwNo30P5CWbWlAsyzCwiBBjrgFXu3mNmtwBfMrNj3f2+IuVeQmpi6tKlS1m4cOEAmikiIo1KAcbgeBiYaWabCd/+V7t7d3LsGGAPoKPMSZov5H1+APi0u28rlDkJDkan0wr0Tvwj8ESJPCcAbcDl7t6TpN0EfJHQi1EswFhACJ4AmD179uQSeaUBNDc375LW1I9hivQ5ixcvLnsviyiKGDduHKeddlrR9ohIfVKAMThOBeYBnwK+Cmw1szuBi4GJSZ71ZZZ1ImEVSTew3t0395H/OODn6QQzO8jd16aSnuhjkues5OedZjYueb+Z0PPxcTM7391fzT/J3bOESaE5r/fRVqlzhXbN7M8+mulzMplMRQFGekhEu3iKNA4FGIMgWfUxB5hjZgcAHybMXZjEjm/4+5dZ3GOpVSTlWA0clZdWaAlqQWY2Hjgl+fj7Itk+AVxXQZtERGSEUYAxyNx9HXCtmZ1AmJvxIKFH4hwzu30Ae1kUq+81oJwlqMWcSZio+UXg13nHIuBHhGESBRgiIlKUAowqM7M9CUMUtxE2s3qN0KMwFbjD3bea2eeAG4BVZnY98CJhdchkdz+vNi3vdTawCbg6GfLYiZktA843s8Pcfc2Qt05ERBqCnkVSfduA3xD2qrgVWEmYi7EQ+CcAd7+RsCvnaEKg8RPCkMq6GrS3l5kdAUwBbi4UXCRyPRezihwXEREh0sY1MogOB1Zv3LiR7u7uPjP3RxRFtLa20tnZOSI2YRrq600/rj23f8WCtc+Syas7G0Vc0nZgWccreSJqFEXss88+PPvssyPice36ex7eGvV6m5qamDhxIsARwO/KPU9DJCJSVBRFNX88eq4NjfQ/ZBHREImIiIgMAvVgiEhFugoMUaTT+jouIiODAgwRqcilBx4woOMiMjJoiERERESqTj0YItKn5uZmOjo6ih7PrTbJ5S210kPPExEZGRRgiEifyllNoseoi0iahkhERESk6tSDISJDKj2c0tfxTCbDtm3byGazu+yDMdw33RJpdAowRGRIdXV10d7ePuByKtkRVESGnoZIREREpOrUgyEiNXPZs+tozhv66Iqi3r008o+nj4lIfVOAISI10xzHuzwYrZLjIlK/NEQiIiIiVaceDBGpSCWbajWK4XhNIrWmHgwRqUhuFUh7e3vJ5aaNZDhek0itKcAQERGRqmv4IRIzmw9cmkrqBp4H7gIudfdXknwzgZuAo9z9kT7KzOUt5Cp3vzAvf3NS517AKe6+vIL2fwq4HniLu69P0r4PvNvdDy63HJFqe/TRR7nttts4/fTTmTJlSq2bIyINZjj1YEwFjgZOApYDnwVWmtlABlPPSspMv75VIN8/EIILgFkDqC/nUuBjVShHpF+y2SzLli1j8+bNLFu2jGw2W+smiUiDafgejJTV7r4peX+PmU0AZgDHAA/0s8zH++rtSMwCssD9wDQza3X3zn7Wibs/3d9zh0IStGXcfVut2yKDY8WKFWzZsgWALVu2sHLlSk4++eQat0pEGklNAwwzmwhcAUwDJgGvAk8RhjZWJXmmAhcBRwJNwLPA99x9QR/FP0QIMA6k/wFGn8zsLcD7gX8DbgFOAD4JXDmAMncaIjGz3QhDP9cAvwcuAQ4g/K4ucfeVeef/FXAZcDywJ/A08G13X5rKszvhd38C0Aa8ATwJLHD3/0jlS9f9FHA+8JfAbOCG/l6j1K8XX3yRFStW9D77I45jVqxYwdFHH83ee++90zNB+tOzkT6n0h0u0vmr2auyU5u074ZIVdS6B2MZcDjwJcLNa1zyeQKAmc0izE/4JXAO8BJwCPD2MsrOzV/YOID2jU5usL3cfXtenrMIQ03fBVYB64GzGUCAUcI/AO8G5gGvA18AlpvZIe7+LICZvZ0QUD0DXAC8CEwHrjWzN7v7FUlZuxN+34uADUCGECgtN7MZ7n5bXt0fI/z+5ydlvpjfODPLJOUAsGjRojFtbW0Dv2oZMnEcc+uttxZNnzt37k6rLObOnTug+rqjiJYKbujdqeWjA627mK6uLj16XqQKah1gvAe4wd2vT6X9GMDMxgJXE26Wx7t77v9C9xYpKxcMjCXMwzgHeA741QDa91B+gpk15YKMZKjgLGAdsMrde8zsFuBLZnasu983gLoLyQAnuPvrSf2PESaXfhz4lyTPN4FXgPe6+2tJ2j1m1pK061p33+Lu/0sIhHLXNZoQIE0A5gD5AcabgL9z9y0l2ncJqQm3S5cuZeHChf27UqmJzs5O1qxZs0t6T08Pa9asobOzk7Fjx9agZSLSaGodYDwMzDSzzYSb22p3706OHQPsAXSkgotSXsj7/ADw6WLzBJLgYHQ6rUDvxD8CT5TIkxteuNzde5K0m4AvEm7e91VQVznuzQUXSRkbzGwTYRgIMxsDvI8wEfXPeb0vKwhB118D9yT5TwX+GZhMCCBythaoe1UfwQXAAkJQCMDs2bMnk/wOpDG0trZy2GGH8cQTT9DT09ObPmrUKA499FBaW1t3Gk5YvHhxxU80zWazvb0PTRUOR6Tz96fuctrU3NxclTJFRrpaBxinErr7PwV8FdhqZncCFwMTkzzryyzrRGALYb7Aenff3Ef+44CfpxPM7CB3X5tKeqKPSZ65FSN3mtm45P1mQs/Hx83sfHd/lRCI3JNXV++y1AoUuqYsYbgDwkqW0cDc5FXIXkn9BtwO/JAwTPICsJ2w+uaMAuf1OWnV3bNJe3JeL5ZX6lMURZxxxhnMmzevYHoURTvtcpnJZAZ0k690iVc6/0DrLlqHdvEUqYqaBhjJqo85wBwzOwD4MGHuwiR2fBPev8ziHkutIinHauCovLQN5Z5sZuOBU5KPvy+S7RPAdcBvCtS1yxyGKvhfoIfQi/KvRfL8T/LzTOD/Aqele4iSeRSFaObbCLH33nszffp07r77buI4Jooipk+fzqRJk2rdNBFpILXuwejl7usIExFPIMzNeJDQI3GOmd1e5jBJJfW9BpSzBLWYMwlzIr4I/DrvWAT8iDBMcl0V6iqLu79mZvcB7yIEXN0lssdAV15wsS/wwUFupjSA6dOnc//99/PKK68wbtw4pk2bVusmiUiDqVmAYWZ7EoYobiMsj3yN8C1/KnCHu281s88RlkKuMrPrCd/6DwYmu/t5tWl5r7OBTcDVydDATsxsGXC+mR3m7rvOmhs85xMmtt5nZv8KrCXMZXkrcJK7n5Dk+wlwnZl9G7iTsOz1K4RenP9vCNsrdSiTyTBjxozenTwHYyhCRIa3Wu7kuY0wdDADuBVYSZiLsRD4JwB3v5GwxHI0IdD4CWFIZV0N2tvLzI4ApgA3FwouEtclP6uxs2fZ3P0PhKW+TxL2ubgHuBH4CDvPA7mBsDz4g4QJoBcBXwN8KNsr9WvKlCksWrRI24SLSL9E2lRGBtHhwOqNGzfS3V1qtKb/oiiitbWVzs7OEbFBUj1c70AfbZ7NZmlvbwdgwdpnyeRdRzaKuKTtwILH08c6Ojqq1rNSL49rr4d/36Gk620MTU1NTJw4EeAI4Hflnlc3czBEpDFEUTTshkyG4zWJ1NpwetiZiIiI1An1YIhIzXQVGIpIp+UfL5RfROqTAgwRqZlLDzxgQMdFpH5piERERESqTj0YIjKkmpub6ejoKHo8vaIjk8nQ2trKCy+8sMusez0zRKS+KcAQkSFVzoqN3OPSoyiipaWFTCbTUMv6RERDJCIiIjII1IMhdateNj8SEZHKqQdD6lZXVxft7e20t7f3BhoiItIYFGBI3cpmswXfi4hI/VOAISIiIlWnAENERESqTgGGiIiIVJ0CDKlbPT09Bd+LiEj9U4Ahdev1118v+F5EROqfAgwRERGpOm201UDMrJy9ki9z9/lJ3iXufl5ybhvwTJFzVrv7kQXquwv4EHCNu8/pZ7NFRGQEUoDRWI4ukr4b8D1gP2BFH2V8G7gtL21rfiYzawWmJR/PNLOL3V27XYmISFkUYDQQd3+oULqZfQs4CPiMuz/cRzHripWT55OEv4+7gZOAkwGvoLkiIjKCaQ5GgzOzGcBngRvd/boqFn02sAGYCWxLPouIiJRFPRgNzMzeBXwH+C1wbpmnjTKz/H/3N9y9d36HmR0LvBVY4O6bzGx5SLYD3H1difZkgN7ncC9atGhMW1tbmc0SEZHhRD0YDcrM9gLuJMyf+Ki7l/uwjoVAd97rhLw8s5Kf301+3kj4W5nZR9mXAFtyr6VLl95XZptERGSYUYDRgMxsNHA7sD9wqrs/V8Hp1wBH5b1+kyp7D+BjwH3u/t9J8r3As8BZZlbqmekLgD1zr9mzZx9bQbtERGQY0RBJY1pE6HW40N1/XuG56939kRLHTwPeBLiZjUulO3ARcDwh4NhF0ouS7knR7lgiIiOUejAajJmdBlwA/NDdrxqEKnLDI9cCL6deF+UdFxERKUoBRgMxs3cCNwCPMwg3ejN7B2HIxIG/K/D6BXCKmY2vdt0iIjK8aIikQSQ39eWEVRoLgXeYWaGsG9396X5WkwtaFrr774q04X3A6cCSftYhIiIjgAKMxjGZsJkWwLIS+W6h79UeuzCzZuBM4JFCwUXiP4BOQiCiAENERIpSgNEg3P0XQKkVHPn5o7zPa0udn2wDvlcfZW4H9i23DQM1ZsyYgu9FRKT+aQ6G1K1Ro0YVfC8iIvVP/9cWERGRqlOAISIiIlWnAENERESqTgGGiIiIVJ0CDKlbmUym4HsREal/WqYqdau5uZmOjo7e9yIi0jgUYEjdiqJIPRciIg1KQyQiIiJSderBkIYTxzHZbJauri6iKGLcuHFks1mampqIorI3OxURkUGkAEMaTldXF+eee+4u6R0dHRpSERGpExoiERERkapTgCEiIiJVpwBDREREqk4BhoiIiFSdAgxpKLkVJIVks1niOB7iFomISCEKMKShdHV1MXfu3ILH5s6dS1dX1xC3SERECtEy1UFgZvOBS1NJ3cDzwF3Ape7+SpJvJnATcJS7P9JHmeV8Nb/M3een8i9x9/OSz23AM0XOW+3uRxao8y7gQ8A17j6njPpFREQA9WAMtqnA0cBJwHLgs8BKM+vPblBHF3m9lxA4dAEryijn2wXKmJmfycxagWnJxzPNTA8DERGRsqkHY3CtdvdNyft7zGwCMAM4BnigkoLc/aFC6Wb2LeAg4DPu/nAZRa0rVlaeTxL+Pu4mBEgnA15mc0VEZIRTgJHHzCYCVxC+vU8CXgWeIgxtrEryTAUuAo4EmoBnge+5+4I+in+IEGAcSIUBRpG2ziD0itzo7tcNtLw8ZwMbCL0bzyWfFWCIiEhZFGDsahlwOPAlQmAxLvk8AcDMZgHXA78EzgFeAg4B3l5G2QcnPzcOtJFm9i7gO8BvgV33zS5ulJnl/7u/4e69czzM7FjgrcACd99kZstDsh3g7usG2vaB6GuViFaRiIjUBwUYu3oPcIO7X59K+zGAmY0Frib0PhyfuinfW6Ss0cnNfCxhmOEcQm/ArwbSQDPbC7gT2Ap81N0Lr9ssbGHySns/sCr1eVby87vJzxuBTxB6My4v0a4M0PswkEWLFo1pa2uroGl962uVSFdXFy0tLVWtU0REKqcAY1cPAzPNbDPhprva3buTY8cAewAd6W/8JbyQ9/kB4NPuvq1Q5mTy5+h0mrtvz8szGrgd2B94v7s/V0Y70q4Bvp+X9sdU+XsAHwPuc/f/TpLvJQwDnWVmXy1x7ZeQWj2zdOlSFi7Mj2VERGQkUICxq1OBecCngK8CW83sTuBiYGKSZ32ZZZ0IbCEsU13v7pv7yH8c8PN0gpkd5O5rU0mLgBOAC919p7xlWt/HktjTgDcBbmbjUulOmHdyPMV7bBYQengAmD179mTgvn60sajm5tKLWfo6LiIiQ0MBRp5k1cccYI6ZHQB8GLiSMOEzd/Pcv8ziHkutIinHauCovLQNuTdmdhpwAfBDd7+qgnIrkRseuTZ5FTpeMMBIhmrSwzWvV7dpEEWlV/j2dVxERIaGAowSkgmN15rZCYS5GQ8SeiTOMbPbyxwmqaS+14CCvQtm9k7gBuBxdgQBVWVm7yAEOA4sLZDlUuAUMxvv7i8PRhtERGR4UICRYmZ7EoYobgOeBF4j3HCnAne4+1Yz+xzhRr/KzK4HXiSsDpmc2zVzENo1nrBRV4YwQfMdZlYo60Z3f3oAVeUCl4Xu/rsi7XgfcDqwZAD1iIjIMKcAY2fbgN8Q9qpoI+xxsY5wU18E4O43mtkG4POEQCMC1gK3DGK7JhM204KwjLaYWyiwK2c5kp06zwQeKRRcJP4D6CQEIgowRESkqEj7BsggOhxYvXHjRrq7u/vMXI5sNkt7e3vR4x0dHWQymaLHG10URbS2ttLZ2Tki9vzQ9Q5vut7G0NTUxMSJEwGOAIp9Ad2FnkUiDaW5uZnFixcXPLZ48WKtIhERqRMKMKShRFFUtIcik8loFYmISJ1QgCEiIiJVpwBDREREqk4BhoiIiFSdAgwRERGpOu2DIQ2nubmZJUuW0NXVRRRFHHDAAbz44os0NTXVumkiIpJQgCENJ4oiWlpaaGlpIYoidt99dzKZTEOtKxcRGe40RCIiIiJVpx4MGTbiOCabzdLV1QWEoZT0vhj5n0VEZPAowJBho6uri3PPPbfo8eG+jbiISD3REImIiIhUnQIMGbb+eOjZPHXo2bVuhojIiKQhEhmFFeTVAAAgAElEQVS24lFNaF2JiEhtqAdDREREqk4BhgwLuRUkpWSzWe2VISIyRBRgyLDQ1dXF3LlzS+aZO3du7xJWEREZXDWfg2Fm84FLU0ndwPPAXcCl7v5Kkm8mcBNwlLs/0keZ5XxNvczd56fyL3H385LPbcAzRc5b7e5HFqjzLuBDwDXuPqeM+nPnHQz8X2CGu38/SfsUcD3wFndfX25ZIiIi9aKeejCmAkcDJwHLgc8CK82sPzsjHV3k9V5C4NAFrCijnG8XKGNmfiYzawWmJR/PNLPmfrQ57cdJXS8NsBwREZGaqHkPRspqd9+UvL/HzCYAM4BjgAcqKcjdHyqUbmbfAg4CPuPuD5dR1LpiZeX5JOF3eTchQDoZ8DKbuwt33whs7O/5Q8HM3uTuf6p1O0REpD4NKMAws4nAFYRv75OAV4GnCEMbq5I8U4GLgCOBJuBZ4HvuvqCP4h8iBBgHUmGAUaStMwi9Ije6+3UDLS/P2cAGQu/Gc8nnfgcYhYZIzOx+YCzwGeBq4F1AJ/Ad4BvuHqfO3xP4CvARYD9CT4gD89JBgZl9FjDgbcCbgKeB7wHfdPftqXy5uucCXwcmA3cAZ/b3GkVEZHgbaA/GMuBw4EuEwGJc8nkCgJnNItwofwmcQ7jRHQK8vYyyD05+DvibvJm9i3Aj/i1QfC/pXY0ys/zf0Rt5N/NjgbcCC9x9k5ktD8l2gLuvG2jb8+xLCACuIsxb+RiwEFgP3Ja0ZyzwK2AfQvD3OPAO4DLg7Wb2gVT7/xK4lTBs1A1MAb5I+Df6dF7d+wO3AFcClwBvVPnaBqTo6pC8dK0iEREZGgMNMN4D3ODu16fSfgy9N7qrCb0Px6duavcWKWt0cjMfSxhmOIfQG/CrgTTQzPYC7gS2Ah9199JrGXe2MHmlvR9Ylfo8K/n53eTnjcAnCL0Zl1fa3j7sBUx1998BmNm9wPuA00kCDGAOcChhMuzvk7R7zawTuB04EbgHID0Z1cxGAfcBLwPXmdmF7v5qqu4JwEfc/b5ijTOzDND7sI9FixaNaWtr6/fFVqLY6pAo3r5LvpaWlqFokojIiDbQAONhYKaZbSbcdFe7e3dy7BhgD6Aj/Y2/hBfyPj8AfNrdtxXKnEz+HJ1OS3frJ3lGE26q+wPvd/fnymhH2jXA9/PS/pgqfw9CL8J97v7fSfK9hGGgs8zsq7lr76snpEzP54ILAHePzewPhCGOnA8CjwF/yKvz/09+vo8kwDCzI4D5hH+rN+fV9VZgderzxlLBReISUiuCli5dysKF+fGZiIiMBAMNME4F5gGfAr4KbDWzO4GLgYlJnnKXWZ4IbCF01a9398195D8O+Hk6wcwOcve1qaRFwAnAhe6+U94yre9jSexphLkLbmbjUulOmHdyPKH3YDfCdaXNYNfgpS+FfidZYPfU572BtgL15ewF4XdF6LH4L+CfCcMkWUKwcU1emRDme/RlAaHXCoDZs2dPTuoYdM3NhRfuxNFuZeUTEZHqGlCAkaz6mAPMMbMDgA8TxugnseNGs3+ZxT2WWkVSjtXAUXlpG3JvzOw04ALgh+5+VQXlViI3PHJt8ip0/F53325m+W39n0Fq0ybgFeCfihzPzWk5hRAcnZLea8PMdtnjI9Fnb0sy/JQegnq9z9ZWSRQVWc2cl140n4iIVFXVlqkmExqvNbMTCHMzHiT0SJxjZrf3Yzigr/peAwr2LpjZO4EbCBMcZxXKM1Bm9g5CgOPA0gJZLgVOMbPx7v5yX5uDVdFPgAsJQxrPlsiX+/foDQiSeRifGsS2iYjICNHvACNZCvlzwuTCJ4HXCDfcqcAd7r7VzD5HuNGvMrPrgRcJq0Mm53bNrDYzG0/YqCtDmKD5DjMrlHWjuz89gKpygcvC9LyIvHa8jzABc8kA6qnU1YTeifvM7JvAHwhzVQ4A/p7Q3keA/0MYRvmBmV1FGBJpJ8ybERERGZCB7OS5DfgNYS7BrcBKwrffhSTd8+5+IzCdcIO7gfDteg5Q7eWbaZMJm2mNJiyj/XWR15f7W0GyU+eZwCOFgovEfxDmLQxKD0oxSc/O3xKu/TOEzb9+CJxH+L0/m+RbQ5igOpGwp8U1hB6h0g/0EBERKUOkfQFkEB0OrN64cSPd3cXmnA5MFEW0traydu1aZs+evdOxJ9/+GQDe9vh3etM6OjrIZDI0qtz1dnZ2jog9PXS9w5uutzE0NTUxceJEgCOAYl+qd1FPzyIR6bfm5mYWL15cMs/ixYu1ikREZIgowJBhIYqiPnsmMpmMVpGIiAwRBRgiIiJSdfX0NFWRqop6ulF/hYhIbSjAkGHrr/7ru31nEhGRQaEhEhEREak69WDIsNHc3MySJUt6n6za3Ny806ROrSARERk6CjBk2IiiiJaWFj2OXUSkDmiIRERERKpOPRgifejp6WHr1q1AGGbRfhoiIn1TgCHSh61btzJ37o5HtDT6duMiIkNBQyQiIiJSdQowREREpOoUYIiIiEjVKcAQERGRqlOAIVJCHMdks9md0rLZLHEc16hFIiKNQQGGSAldXV184Qtf2Clt7ty5vbuFiohIYQowREREpOq0D0adM7Ny+uIvc/f5Sd4l7n5ecm4b8EyS5zR3vz2v7PnApcBEd9+UpN0MfDKVrQfoBB4ELnf3x/t/NSIiMlKoB6P+HV3k9V5C8NAFrCijnCvMrKnMOv+cquc4YB5wOPCgme1XUetFRGREUg9GnXP3hwqlm9m3gIOAz7j7w30UsxKYBpwDfLuManvy6r3fzNYB9wInAdeVUYaIiIxg6sFoQGY2A/gscKO7l3Oz/xnwU+DLZvYX/ax2S/Kzu5/nN6Riq0W0ikREpDT1YDQYM3sX8B3gt8C5FZz6eeD3wEXAV8qoJ/e3sRtwMPAN4GXg7hLnZIDeh3QsWrRoTFtbWwVNrD/FVot0dXXpsfAiIiWoB6OBmNlewJ3AVuCj7p7t45Re7v4YcBtwgZnt00f2MYSeim7CfIw/AG8DPuTuL5U47xJCT8cWYMvSpUvvK7d9IiIyvCjAaBBmNhq4HdgfONXdn+tHMfOAJsLKkVL+DByVvP4G+AjwFLDCzI4ucd4CYM/ca/bs2cf2o411pbm5uaJ0EREJNETSOBYBJwAXuvvP+1OAu681sw7gPDO7ukTWHnd/JJ1gZj8FngOuJqwuKVR+Fkj3qrzen3bWkyiKKkoXEZFAPRgNwMxOAy4AfujuVw2wuK8BfwK+XslJ7v4n4Glg8gDrFxGREUABRp0zs3cCNwCPA7MGWp67bwYWAh8D/rqCdowlTPYsNQdDREQE0BBJXTOz8cBywsqMhcA7zKxQ1o3u/nQFRX+TsAJlWpHjo8zs3bn3wH7A+cB4YH4F9YiIyAilAKO+TSZspgWwrES+W4CZ5Rbq7n9KtgkvtofG7sCvU59fAp4ATnH35eXWIyIiI1ekDYNkEB0OrN64cSPd3YOzP1cURbS2ttLZ2Tkom1/FccymTZt2eqLq4sWL+Yu/+IuaTPQc7OutN7re4U3X2xiampqYOHEiwBHA78o9T3MwREqIoohMJrNTWiaT0SoSEZE+KMAQERGRqlOAISIiIlWnAENERESqTgGGiIiIVJ2WqYr0YezYsSxevBgIzyDRc0hERPqmAEOkD6NGjWKPPfaodTNERBqKhkhERESk6tSDIdJPcRzT1dVVMr25ubngnhnF0kVEhgsFGCL91NXVRXt7e7/O7ejo2GUDLxGR4URDJCIiIlJ16sEQqYINJ28g3i08WyDaHrHv8n1LpouIDHcKMESqIN4t7g0kykkXERnuNEQiIiIiVaceDJE85awCqUYd2Wx2UOsQEakl9WCI5MmtDmlvby+4DLVR6hARqSUFGCIiIlJ1GiKpEjObD1yaSuoGngfuAi5191eSfDOBm4Cj3P2RPsosZ3bgZe4+P5V/ibufl3xuA55J8p3m7rcXafNEd9+UpN0MfDKVrQfoBB4ELnf3x8tokwyxRx99lNtuu43TTz+dKVOm1Lo5IiLqwRgEU4GjgZOA5cBngZVm1p9B9qOLvN5LCBy6gBVllnWFmTWVmffPqbqOA+YBhwMPmtl+ZbdehkQ2m2XZsmVs3ryZZcuW9c7tEBGpJfVgVN/qXG8AcI+ZTQBmAMcAD1RSkLs/VCjdzL4FHAR8xt0fLqOolcA04Bzg22Xk78mr+34zWwfcSwicriujDBkiK1asYMuWLQBs2bKFlStXcvLJJ9e4VSIy0inAAMxsInAF4SY8CXgVeIowtLEqyTMVuAg4EmgCngW+5+4L+ij+IUKAcSAVBhhF2jqD0Ctyo7uXe6P/GaG36stmdrO7v9aPqrckP7v7cW5DieMdI1OlegN2OlbOYFYqT/rcdH2VevHFF1mxYkVvGXEcs2LFCo4++mj23nvvfpcrIjJQCjCCZYQhgC8RAotxyecJAGY2C7ge+CWhF+Al4BDg7WWUfXDyc+NAG2lm7wK+A/wWOLfC0z8P/J4QJH2ljLpyfxu7Ea7hG8DLwN0lzskAvQ/YWLRo0Zi2trYKm1l76VUdc+fOLeuc6I2IuKl0oBC9sWOU7Atf+MJO9bW0tFTYyhBM3HrrrUXT586dq+WvIlIzCjCC9wA3uPv1qbQfA5jZWOBqQu/D8e6eu4vcW6Ss0cnNeSxhOOEc4DngVwNpoJntBdwJbAU+6u4VDbS7+2NmdhtwgZl1uPsLJbKPYdeeik7gQ+7+UonzLiE10XXp0qUsXLiwkmZKBTo7O1mzZs0u6T09PaxZs4bOzk723Vdbk4tIbSjACB4GZprZZmAVYR5F7gZ7DLAH0JEKLkrJv3E/AHza3bcVypxM/hydTnP37Xl5RgO3A/sD73f358poRyHzgI8TgoDZJfL9GTg2eT8K2A/4Z2CFmU11918XOW8BIRgDYPbs2ZOB+/rZ1pppbm7ufb948eKiTz3NZrO9PRzx6L7/NNJ5rrzyyt5ejHR9lWhtbeWwww7jiSeeoKenpzd91KhRHHroobS2tvarXBGRalCAEZxKuPl+CvgqsNXM7gQuBiYmedaXWdaJhPkK3cB6d9/cR/7jgJ+nE8zsIHdfm0paBJwAXOjuO+WthLuvNbMO4Dwzu7pE1p78JbRm9lNCT8zVhNUlhcrPAumeldf729ZaSg8rZDKZ8h6rXs5IRCpPusz+DmNEUcQZZ5zBvHnzCqZreEREakkBBpCs+pgDzDGzA4APA1cSJnzmbsT7l1ncY6lVJOVYDRyVl7Yh98bMTgMuAH7o7ldVUG4xXwPOBr4O7Nq/XoS7/8nMngYmV6ENUiV7770306dP5+677yaOY6IoYvr06UyaNKnWTROREU4BRh53Xwdca2YnEOZmPEjokTjHzG4vc5ikkvpeAwpuuGVm7wRuAB4HZlWpvs1mtpCwamZMueclc1EOJkxwlToyffp07r//fl555RXGjRvHtGnTat0kEREFGGa2J2GI4jbgSeA1Qo/CVOAOd99qZp8j3OhXmdn1wIuEm+3k3K6Zg9Cu8YSNujLAQuAdZlYo60Z3f7rC4r9JWIVS7E40yszenXtPmINxPjAemF9hXTLIMpkMM2bM6N3Js6whHRGRQTbiAwxgG/Abwl4VbYQ9LtYRbuqLANz9RjPbQFjqeQNhNH0tcMsgtmsyYTMtCMtoi7kFmFlJwclwx3yKb5i1O5CeyPkS8ARwirsvr6QuGRpTpkzRFuEiUleigWzyI9KHw4HVGzdupLt7cPbniqKI1tZWOjs7B7RhVVq5j2vPZrO0t7cD8PzHnifeLdQfbY/Y79/3K5m+ZMmS3nIreVz7YFxvPdP1Dm+63sbQ1NTExIkTAY4AflfueerBEMkTRdGgDzMMRR0iIrWkh52JiIhI1akHQ6QKou1RRe9FRIY7BRgiVbDv8sJbchdLFxEZ7jREIiIiIlWnHgyRfmpubqajo2OX9HJWofT3+SMiIo1CAYZIP5VaCdKfx6+LiAwnGiIRERGRqlMPhsggiOOYbDZb1oZdfRnIuSIitaIAQ2QQdHV1ce6551alrI6ODm3KJSINR0MkIiIiUnXqwRAZAhtO3tD7TJJyRNsj7aEhIg1NAYbIEIh3iysKMEREGp2GSERERKTq1IMh0g/lPtJ9pLZHREQ9GCL90NXVRXt7O+3t7b03drVHRGQHBRgiIiJSdTUZIjGz+cClqaRu4HngLuBSd38lyTcTuAk4yt0f6aPMcmbQXebu81P5l7j7ecnnNuCZJN9p7n57kTZPdPdNSdrNwCdT2XqATuBB4HJ3f7yvBqXqPcvdb07SZhKu+yB3X1vGdYmIiNSVWvdgTAWOBk4ClgOfBVaaWX8GkI8u8nov4QbeBawos6wrzKypzLx/TtV1HDAPOBx40Mz2K7v1O7s7Ka+zn+eLiIjUVK0nea7O9QYA95jZBGAGcAzwQCUFuftDhdLN7FvAQcBn3P3hMopaCUwDzgG+XUb+nry67zezdcC9hMDpujLK2Im7bwQ2VnreUDKzN7n7n2rdDhERqU8VBxhmNhG4gnATngS8CjxFGNpYleSZClwEHAk0Ac8C33P3BX0U/xAhwDiQCgOMIm2dQegVudHdy73R/4zQs/NlM7vZ3V/rR9Vbkp/d/Ti34BCJmf0C2As4C7gaOAJ4gRDALHL3ntT5ewBfAT4K7EcIVv4N+JK7v57Kdy5wKvA2YAzwP8AyYLG7d6fy5epuB64EphCGsz7Rn+sbDuJ4x4hcNpvd5XihtMoqKFxWFEVs27aNbDZbtA3pdBGRWulPD8YywhDAlwiBxbjk8wQAM5sFXA/8ktAL8BJwCPD2Mso+OPk54G/vZvYu4DvAb4FKHwrxeeD3hCDpK2XUlfs97ka4hm8ALxOGOqppH+BW4CrgMuAUYAGwAfhe0pY3EX73+wNfB/4TOAy4HHiHmZ3o7rk70F8Ct7FjCGky4d/1bcDZeXW3At8HFgFfJMw32YmZZYDeh2YsWrRoTFtb2wAvuT6lV2rMnTu36uVHb+wYJay0/K6uLj0uXkRqrj8BxnuAG9z9+lTajwHMbCzh2/UDwPGpG9m9RcoandycxxKGE84BngN+1Y929TKzvYA7ga3AR929oq+T7v6Ymd0GXGBmHe7+QonsY9i1p6IT+JC7v1RJvWWYAExPDfWsMrP3AaeTBBjA+cA7gb9JTYy918yeB/6dMO9lJYC7X5Ar2MxGEX7vm4GbzOxz7v5yqu43Ax9395+VaN8lpCbvLl26lIULF/brQkVEpLH1J8B4GJhpZpuBVYR5FLkb7DHAHkBHKrgoJf/G/QDwaXffVihzMvlzdDrN3bfn5RkN3E74Bv9+d3+ujHYUMg/4OOGGObtEvj8DxybvRxGGJP4ZWGFmU93910m78n/Xb5T5O0p7ocA8kv8kDFnkfBB4HHg0r86fEjre30cSYCS9PJcRgsY355V7CPCb1OeX+wguIPSmXJ37MHv27MnAfX2c05Cam5t73y9evHiXp51ms9kB9WzEo3f8aaTLj6KIffbZhxdeeGGXIZJcfem2iYjUSn8CjFMJN99PAV8FtprZncDFwMQkz/oyyzqRMF+hG1jv7pv7yH8c8PN0gpnlL+VcBJwAXOjuO+WthLuvNbMO4Dwzu7pE1p78JbRm9lNCT8zVhNUgsGsvx1nAzRU2q9DvJwvsnvq8N2GYptj8j72SNh5A6LH4IyEgWgtsA/4aWJJXJpSxoiXpKUr3Fr1eLG+jS++Umclkqv849dQ6qnT5URTR0tJCJpMpOtdCu3iKSD2oOMBIVn3MAeYkN6kPEyb+TWLHt9f9yyzusdQqknKsBo7KS9uQe2NmpwEXAD9096sqKLeYrxHmInwdWFPuSe7+JzN7mjCnISe/3c8wODYRelXy51CkjwOcTBje+Yi7P5s7aGZTCp6107RDERGR0ga0TNXd1wHXmtkJhG72Bwk9EueY2e39GALoq77XgIIbbpnZO4EbCMMDs6pU32YzW0hYNTOm3POSuSgHEya45soquVFYFf2EMAlzs7uXCmJy/za9PQ7JENQ/DWLbRERkhKgowDCzPQlDFLcBTwKvEb6ZTwXucPetZvY5wo1+lZldD7xIuNlOzu2aWW1mNp6wUVcGWEhYLVEo60Z3f7rC4r9JWIUyrcjxUWb27tx7whyM84HxwPwK66qGbxKWp95nZosJczRGAQcAfw9c5e6/Ae4hrBz5gZktAloIc03G16DNIiIyzFS6k+c2wsS/GYTlkisJczEWknzzdfcbgemEyZg3EL5RzwHWVafJBU0mbKY1mrCM9tdFXl+utOBkM6n5JbLsnir/AeBa4A3gFHf/VqX1DVSyz8V7CfM7Pk1YKuuEoGc9Ya4F7v4kIRAZD9xB2FTs0SSfiIjIgETalEcG0eHA6o0bN9Ld3a89x/oURRGtra10dnYO6QZTfT0ePZvN0t7e3vv5+Y89T7xb+e2Ltkfs9+9hp/mOjo6dJnkWut7h+rj2Wv371oqud3hr1Ottampi4sSJEDZ4/F2559V6q3CRhhRFUfVXjgxAvbVHRKTWDzsTERGRYUgBhoiIiFSdhkhEhkC0vbI5EZXmFxGpNwowRIbAvsv3rXUTRESGlIZIREREpOrUgyEyCJqbm1myZElVlo7q4WUi0ogUYIgMgtxDyVpaWmrdFBGRmtAQiYiIiFSdejBE6lB6Z860KIrYtm0b2Wx2UHcCrJedQQfjeofTTqci9UwBhkgd6urq2mmrcame9NbrIjJ4NEQiIiIiVaceDJE6t+HkDRU9KK0aou1R794dtai/2tLXIyJDQwGGSJ2Ld4treoOvdf0i0pg0RCIiIiJVpx4MkRqpl5UaIgOlv2UpRD0YIjWSWynS3t5ecEmqSKPQ37IUoh6MEcrMfg9MAg5w9zeK5HkA+Ctg3+T1DHCRu//LkDVUREQaknowRq4bCUHDBwodNLNDgGOAZe6uryQiIg3o0Ucf5eKLL+bRRx8d8roVYIxctwLbgLOLHM+lf3domiMiItWUzWZZtmwZmzdvZtmyZWSz2SGtX0MkDcrMJgJXANMIQx2vAk8Bl7r7qiTPVOAi4EigCXgW+J67L3D3l83sTuCjZjbB3Tenyh4NzAB+6+5/GMrrEhGR6lixYgVbtmwBYMuWLaxcuZKTTz55yOpXgNG4lgGHA18iBBbjks8TAMxsFnA98EvgHOAl4BDg7akybgROA84Erkmlf4AwfHL5oF7BCJd+tkb+N4udPmsLioFL/Q6H+ltcJYbqWTPVlv6dNlK7h7MXX3yRFStW9P57xHHMihUrOProo9l7772HpA0KMBrXe4Ab3P36VNqPAcxsLHA18ABwvLvn/ou/N6+MnxEmbp7NzgHG2cCfgB9U0iAzywC9D3lYtGjRmLa2tkqKGFHSs+3nzp1bNF/0RkTcpP9pD0T0xo5lk6V+1zJwXV1dtLS01LoZI1ocx9x6661F0+fOnTskS4k1B6NxPQzMNLN5ZvZuM2tKHTsG2APoSAUXu0iO3QS808yOADCzCcCHgB+5+6sVtukSYEvutXTp0vsqPF9ERAaos7OTNWvW0NPTs1N6T08Pa9asobOzc0jaoR6MxnUqMA/4FPBVYGsyp+JiYGKSZ30Z5dwEzAfOAlYDZwDNhOGTSi0g9JwAMHv27MmAgowimpube98vXrx4pyd8ZrPZ3m/a8Wj1XgxU+neY/7uuJ1EUsc8++/DCCy801FBD+u81/XcttdHa2sphhx3GE088sVOQMWrUKA499FBaW1uHpB0KMBqUu28C5gBzzOwA4MPAlYQJn7mb/P5llLPezP4PcLqZfY4QaPw3/QgM3D0LpAe4X6+0jJEk3UWZyWSK3/S0KeLApX6HJX/XNRZFES0tLWQymYYKMNK0i2ftRVHEGWecwbx58wqmD9W/kYZIhgF3X+fu1wL3ECZ6PkgYpjjHzMr5S7oRGE+Y1DkFuKnU0IqIiNS3vffem+nTp/cGE1EUMX36dCZNmjRkbVAPRgMysz2BnwO3AU8CrwFHAVOBO9x9a9IbcQOwysyuB14EDgYmu/t5eUXeBWwiLGl9A7hlSC5EREQGzfTp07n//vt55ZVXGDduHNOmTRvS+tWD0Zi2Ab8h7FVxK7CSMBdjIfBPAO5+IzAdGE0INH5CGFJZl19YslPnMkJH8k/d/fnBvwQRERlMmUyGGTNmMGHCBM4888whHxpUD0YDSuY6zC4j30pC8FFOmRcAF5Q4vhbNBhARaShTpkxhypQpNalbAYZIjTQ3N9PR0dH7XqRR6W9ZClGAIVIjURTV7WoGkUrob1kK0RwMERERqTr1YIjUuWj70E99SddZi/qrbThcg0ijUYAhUuf2Xb7viK5fRBqThkhERESk6tSDIVKH0rPy04bqWRVxHPc+7bW5ublm2z8PxvVqlYPI0FCAIYOpBWC33Qb/z2wo6hhqxW6ELS0tjBkzZohbUzsj7XpheP49l6LrrW+p9rZUcl7UqA/UkYZwOmGnURERaXxnEB5RURYFGDKYJgAfANYStjevurVr145ZunTpfbNnzz62ra1t2D+9Vdc7vOl6h7cGvt4WoA34KbC53JMUYEhDM7M9CE+O3dPdX611ewabrnd40/UObyPterWKRERERKpOAYaIiIhUnQIMaXRZ4LLk50ig6x3edL3D24i6Xs3BEBERkapTD4aIiIhUnQIMERERqToFGCIiIlJ1jbVfqYwIZnY8cCZwDPAW+H/t3Xv4VFW9x/E3WFmZlzzkU3FSwXta3lLT1ALB8oZh8aUUb2H2ZN5QT4aikKGeo9lFQNMUMTU730pLyYwQIfMRI8rM8pgdQ0TMC4lJGijYH981Og575jcze+8fM/F5Pc/v2fz2nrX3+s7w2/s7a6+1NkuBXwPnufv8HsoeBowAdgP6A08CdwMT3P3hMuvdrjzxZuxrInA28Ad336HouhahiHjN7FDgNGBnYB1iMrdvuvuVZdQ5j7zxmtkg4CxgR5dnL+wAAA2iSURBVOCtwCPAVcAUd19ZVr3bZWY7AecD7wPeAbwIPETU9/omym8CXAQcTMT7O2Ccu99RWqVzyBNvN56vWqEWDOlEnydmjfsmcCBwCrAJMDedrBs5kzgpnQ98DBhHXIR+Y2bbl1XhnPLE+6p0ojuDOEl1slzxmtmXgJuABwADhgGXAZ36FLO24zWzIcBM4svgZ4GPA7PTvr5WWo3z2Qh4jEiKDgSOIhLA68xsXKOCZrYucAewH/E+HUr8f77dzD5cYp3zaDteuvN81TS1YEgn+oK7P1W9wsxuB/5M/BHPalD2kIyys4g/+DHAccVWtRB54q28/g3ANcAVxDfdfiXUsyhtx2tmuxIn47HuflHVpo78dpvk+XyPAV4CDnb3ytTSM81sm7TtlMJrm5O7zyaSoGrTzWwAcDwwsUHx0cAOwF7ufg+Amd1JtGJcBOxRdH3zyhlvN56vmqYWDOk4tX9wad0y4I9EE3OrZRcDi3oqu6bkibfKl4CNidsjHS1nvCcScwhMKqFqpcgZ70vACqLZvdpSSnq+T4meAV7u4TXDgYcqyQWAu78MXA/sbmb9S6xf0XqMtxvPV61QC4Z0BTPbENiFJr7NZ5QdCGwG/KjoepWllXjN7L1E0+ph7r7MzMquXuFaiHdf4EHgE2Z2DrAl8ARxATrX3VeUWtGCtBDvt4BPA5ea2QXAC8AhxIV4bKmVzMnM+hJfYt9O9DP4KJEgNrIDcFfG+vvTcnvg8aLqWKQ2483aT9edr+pRC4Z0iynAekTzeNPSrYOrgWXA10uoV1maijed1KYCN7n7bb1RsZI0+/n2B7YCLk0/Q4BpRN+Ta0qsX9Gaitfd7wUGEwnF48CzRJxnu/slZVcyp8uIFpiniL+9k939ih7K/Afwt4z1f6va3qnaifd1uvh8lUkJhnQ8M/sKcAQwppVRFWbWh/hj3Qc4yt0fK6mKhWox3tOIC+6ppVesJC3G2xdYHzjB3ae4+53uPo64ZXK4mW1ZcnVzayXe1OfkZmA+0XIxGLgQmJhacDrZBcToiIOIJHiymZ3RRLlG00t38tTT7cYLdO/5qhHdIpGOZmbjieb/s919cgvl+hBD+UYBR7v7j0uqYqFaidfMNgXOI/pfrDCzjdKmNwB90+/L3b32/n3HaOPzXQK8E/hZzfqfEknWLkTnyY7URrxTiFEUw6uGpN5pZquACWZ2g7s/UlJ1c3H3hcDC9Ott6dbdhWZ2rbs/XafYErJbKTZOy6zWjY7QZrxA956veqIWDOlY6WQ8gRgTfkEL5Sp/rMcCxzUz9r4TtBHvQOAtxJDFZ6t+PgRsl/59YSmVLUCbn+/9ddb3SctVeetVljbj3QmYnzHfxTzi/L1dcTUs3a+I5Hdgg9f8nphPolZl3QNFV6pEzcTbteerZijBkI6Umn8nABPd/cstlOsDfJv4Y/2cu3fFffk2470PGJTx8ztimNsgoOlWn97U7ucL/DAtD6hZfyCRXMzLX7vi5Yh3MfABM1unZv2eabmogOr1lkHEZ9SoxeVmYFsze3U4auqXMAq4N42w6BY9xtut56tm6Wmq0nHM7HTgq8DtxKONX8fd56bXXQ0cDWzh7o+mdZOInttTiT/casvd/bclVr0teeKts7/ZQL8Onskzz+f7RmAusAVwLjHUcwjRyfNydz+pN2JoRc54TyI6s/6UmOPkBWISqtOB2e4+tDdiaIWZXQn8nfgG/yQxJ8sIYCRwsbt/Mb0uK951if4mGxC3/p4CTiD6nwxx9zm9G03PcsbbdeerVqgPhnSiQ9LyY+mnVqU5fJ3006dqW6XsZ9JPtUeJGRU7TZ54u1Hb8br7S2Y2lOhQdxZxb/4vxMWoU2e2zBPvJDN7nJh06SriltgCIlHp1FEG9xDfyI8mZrlcRrSqHVnT/J8V73Iz24+YVGsSMcvlfcABnZhcJG3HS3eer5qmFgwREREpnPpgiIiISOGUYIiIiEjhlGCIiIhI4ZRgiIiISOGUYIiIiEjhlGCIiIhI4ZRgiIiISOGUYIiIiEjhNJOnyFrIzF4BcPdunxV0jTGzY4BrgC+7+4Q1WxuRzqMEQ0S6nplNAMYDx7r7tIL2OY2Y/nmQu88uYp8iaxMlGCIi7bmZePDaM2u6IiKdSAmGiEgb3P054Lk1XQ+RTqWHnYmshbL6YJjZ5sSTSecABwFfIR473Q94EBjv7rem144A/gvYnnh65P8CZ7r7izXHWQBsRnQoPxn4HDAQWAL8EDjX3ZfWlJkNfBgY4O4Lara9Wkd3/0jNMbIMcvfZZvZm4AhgGPA+4F3AcuB+4DJ3/17W+1PHAHdf0KgPhpm9lXik+kji0fIriCdsrnas6hjcvY+ZHUe8V1sRCcyPifd2aW25RtJ7NRb4aIr378Rne56731/z2ldjAb4LnAcMAt4BHObuP6r+XIC9Uh23B1a6+0ZV+9qTeLrtXsRj158AbgMmuvviVo/bSszSWTSKRERqvQm4AziSeFT2XGBH4GYzG2JmY4iLwcvADOIR1CcRjxOvZxJwMbCIuGBWyswxs/Vz1vcHxMUb4G7g2qqfv6b1m6f67QEsTHW4D/ggcGPqw1HtWuD/079/VrPPZY0qk+L5BXGx3ASYnuq1ezrWNxqUvQiYQiQDtxOP9j4euMXMmu6Qa2Z7E+/J8am+twAPA4cBc81sUJ2i2wDzUl3vBH4OvFTzmrHAdUTSNB14oOq4o4C7iMeQPwTcRCRynwd+Y2bb5jiudBndIhGRWnsCs4Gt3f1ZeN03zcuBjYHB7n5X2vZu4LfA4WZ2jrs/krHPI4E93X1+KvM24iI/mPj2elq7lXX3M1KCsCNwVZ1Onk8T3+RnuvuqykozGwDMAs4xs2mVFhN3PyZ18twC+O8WO3leAOwKzASGu/uydKxtiRaEU8xshrvfllF2FLCHu9+XyvQD7gH2Ib7Zz+rp4Ga2AfB94C3ACHf/QdW2IcBPgOvMbKC7r6gp/ilgMnCqu6+sc4ijiM9/Ts1x3wNcCbwCDHP36Wl9X+AS4FTgO0QSUauZ40qXUQuGiNRaCXy2klwk3yEu0lsCkyvJBUBq9r4h/bpvnX1OriQXqcwy4ETiYjTazNYtsP6rcfcl7j6jOrlI6/8CnE+cCw/JexwzWw8YDawCTqgkF+lY/wdMTL+eXGcX51SSi1TmGSKpg/rvba3PAO8EvlqdXKT9zQQuA/oDB2eUfZq4HdPoIn91bXKRHEckNTdWkot0zFXELZPFwG5m9sE2jytdRgmGiNRa4O5/rl6RLhKPpl9/nlGmcjvhXXX2uVq/A3d/kGjG3wB4f3tVbY2Z7W1m48zscjO7JrVSjEibtyrgELsSF9lfufvDGduvS8sP1bnlMSNj3Z/Sst57W2toWtbrv/DLtNwtY9tMd3+hh/3fUmf9Pml5Q+0Gd19OtKpUv67V40qX0S0SEan1eJ31/2iwvbKtXkvEo3XWLwB2At7dVM3aZGYbEv0BBjd4Wd6+IPBaHAuyNrr7UjN7DtiQSKxqR6EsyihWaQVptpVn87S818wava5fxrqFTey/3msaxl61Puuzbua40mWUYIhIrZ6GlhU59KzVmUTbbXX9HyK5+AVwLtExcam7rzSz/YmOnEXOatrMe7Taa9y9iPd2nbT8PtCoVeDejHX/bGL/Pb2mnf8/zRxXuowSDBHpDZsBv89Yv2laVg9frHQ8fFvG69/T5vGHE31LhqX5K6oNbHOfWSpxDMjamFpSNiRafJ4v8LjVFhGjMibWDkct2eJ03AG8dlunWmUo8RO9ViNZo9QHQ0R6w8jaFWlUxU7Ehbb6Qli5AG2dsZ/96+y/kpTU+9L0duD5jOQCoN59hJ72mWU+8CKwu5ll9ekYlZa/LKi1IsvMtPx4Sfuvp9Lx94jaDWb2Jl7r63JX7Xb596QEQ0R6w4lmtnPllzTaYhJxW2Jq6gRYURmhcHqasKpSZggx1DFLpeVgmzrb/wRsZGavS3TSnB715oToaZ+rcfd/AFOJc+uUFGflWFsD49Kvk5rdZxuuIEZlnGVmx9Z2JjWz9czsKDP7z4KPezWRXH3azA6qOl5fYuhuf2Ceu88t+LjSoXSLRER6w/VEp8NZRMfGfYmhlH8gHlJW7Ubgi8RMkA+a2Tzi1sgHgK8BZ2TsfwZxH3+Mme1AJAevABe7+0PAhakO3zOzLxC3EXYEtgW+DozJ2OetRH+NS8xsKK89c+RMd1/SINaxxAReQ4FHzGwOsB7RB+TNwKXu/pMG5XNx92fNbDgx2mMqMN7MHiAmvNoU2C7VZ2eyO5W2e9yFZnY8MA241czuBh4DdiGStCeJOTRkLaEWDBHpDScRF97NgEOJi/8UYJ/a2xZpuvH9iERjfeBA4lw1MpVZTZqL41Bi1tG9ibkgRpOGdrr7DcT053OJ2zIHEEnIYOoMu0zzdowikqD90/5G08NoE3d/nphSezyRlAwjhmb+Gjjc3U9pVL4I7n43MSX6JUSrwmAihg2I2TdHAn8s4bjXE8njdCKR+SQxbPdyYNc0F4isJfQsEhEpTfUzNtZ0XUSkd6kFQ0RERAqnBENEREQKpwRDRERECqc+GCIiIlI4tWCIiIhI4ZRgiIiISOGUYIiIiEjhlGCIiIhI4ZRgiIiISOGUYIiIiEjhlGCIiIhI4ZRgiIiISOH+BeG9CzjD0vRnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_repeat = 10\n",
    "model = list(\"FA\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-FA\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-FA-linear\" for _ in range(n_repeat)) + \\\n",
    "        list(\"ZIFA\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-ZIFA\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-ZIFA-linear\" for _ in range(n_repeat)) + \\\n",
    "        list(\"ZINB\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-ZINB\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scPI-ZINB-linear\" for _ in range(n_repeat)) + \\\n",
    "        list(\"scVI\" for _ in range(n_repeat)) \n",
    "error = list(fa) + list(fa_nln) + list(fa_ln) + \\\n",
    "        list(zifa) + list(zifa_nln) + list(zifa_ln) + \\\n",
    "        list(zinb) + list(zinb_nln) + list(zinb_ln) + list(scvi)\n",
    "\n",
    "d = {'model': model, 'error': error}\n",
    "df = pd.DataFrame(data=d)\n",
    "fig = plt.figure(dpi=100)\n",
    "grid = plt.GridSpec(1, 10, wspace=0.2, hspace=0.2)\n",
    "ax = plt.subplot(grid[:, 2:10])\n",
    "sns.boxplot(y=\"model\", x=\"error\", data=df)\n",
    "ax.artists[0].set_facecolor('tab:red')\n",
    "ax.artists[1].set_facecolor('tab:red')\n",
    "ax.artists[2].set_facecolor('tab:red')\n",
    "ax.artists[3].set_facecolor('tab:blue')\n",
    "ax.artists[4].set_facecolor('tab:blue')\n",
    "ax.artists[5].set_facecolor('tab:blue')\n",
    "ax.artists[6].set_facecolor('tab:green')\n",
    "ax.artists[7].set_facecolor('tab:green')\n",
    "ax.artists[8].set_facecolor('tab:green')\n",
    "ax.artists[9].set_facecolor('tab:green')\n",
    "plt.xlabel('Imputation error', fontsize=15)\n",
    "plt.ylabel('', fontsize=15)\n",
    "plt.xticks(fontsize=12)\n",
    "plt.yticks(fontsize=12)\n",
    "plt.title(\"Cortex\", fontsize=18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}