042_compute_PMBC_Python.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import numpy as np\n",
    "from scipy import sparse\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from tools import  log_transform, normalize_and_scale, pearson_residuals, sqrt_lazy, PCA_sklearn, kobak_tsne, run_glmpca"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Prepare Rpy2-bridge needed to run some of the GLM PCAs in R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rpy2 import robjects\n",
    "import rpy2.robjects.packages as rpackages\n",
    "from rpy2.robjects import numpy2ri\n",
    "\n",
    "numpy2ri.activate()\n",
    "\n",
    "utils = rpackages.importr('utils')\n",
    "utils.chooseCRANmirror(ind=1) \n",
    "if not rpackages.isinstalled('glmpca'):\n",
    "    utils.install_packages(\"glmpca\")\n",
    "glmnet = rpackages.importr('glmpca')\n",
    "glmpca = robjects.r['glmpca']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(33148, 16809)\n"
     ]
    }
   ],
   "source": [
    "folder = 'datasets/33k_pbmc/'\n",
    "picklefile = folder+'preprocessed.pickle'\n",
    "dataset = pickle.load(open(picklefile, \"rb\"))\n",
    "print(dataset['counts'].shape)\n",
    "dataset['labelshort'] = '10X v1 33k PBMC'\n",
    "counts = dataset['counts'].toarray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load sctransform results and convert to fast-loading pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sct_residuals_df = pd.read_csv('datasets/33k_pbmc/sctransform_pearson_residuals.csv')\n",
    "with open('datasets/33k_pbmc/sctransform_pearson_residuals.pickle','wb') as f:\n",
    "    pickle.dump(sct_residuals_df, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analytical Pearson residuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 10.3 s, sys: 7.18 s, total: 17.5 s\n",
      "Wall time: 17.5 s\n",
      "CPU times: user 11.7 s, sys: 2.52 s, total: 14.2 s\n",
      "Wall time: 14.2 s\n"
     ]
    }
   ],
   "source": [
    "%time dataset['residuals_theta100'] = pearson_residuals(counts,100)\n",
    "%time dataset['residuals_theta10'] = pearson_residuals(counts,10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sqrt transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4.14 s, sys: 237 ms, total: 4.38 s\n",
      "Wall time: 4.41 s\n",
      "CPU times: user 7.38 s, sys: 4.25 s, total: 11.6 s\n",
      "Wall time: 11.6 s\n"
     ]
    }
   ],
   "source": [
    "%time dataset['sqrt_lazy'] = sqrt_lazy(counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Log transform and bare normalization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.34 s, sys: 1.09 s, total: 7.43 s\n",
      "Wall time: 7.43 s\n",
      "CPU times: user 6.65 s, sys: 657 ms, total: 7.3 s\n",
      "Wall time: 7.3 s\n",
      "CPU times: user 6.46 s, sys: 3.32 s, total: 9.78 s\n",
      "Wall time: 9.78 s\n"
     ]
    }
   ],
   "source": [
    "%time dataset['normalized'] = normalize_and_scale(counts,scale=1)\n",
    "%time dataset['log1p_medianscale'] = log_transform(counts,scale_mode='median')\n",
    "%time dataset['log1p_10kscale'] = log_transform(counts,scale=10000,scale_mode='user_provided_scale')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run GLM PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def glmpca_R(data,fam,tol=10**-4,dataset_label=''):\n",
    "    \n",
    "    basepath = 'glmpca_results/glmpca_%s_%s_tol%s' % (dataset_label,fam,tol)\n",
    "    print('will save at', basepath)\n",
    "    \n",
    "    #create ctl list R object for controlling glmpca\n",
    "    list_r = robjects.r('list')\n",
    "    ctl = list_r(tol=tol,verbose=True)\n",
    "\n",
    "    #make use of sparsity\n",
    "    Y_sparse = data.T.tocoo() ## transpose needed because glmpca wants rows as features\n",
    "    r_Matrix = rpackages.importr(\"Matrix\")\n",
    "    Y_r = r_Matrix.sparseMatrix(i=robjects.IntVector(Y_sparse.row + 1),\n",
    "                                j=robjects.IntVector(Y_sparse.col + 1),\n",
    "                                x=robjects.IntVector(Y_sparse.data),\n",
    "                                dims=robjects.IntVector(Y_sparse.shape))\n",
    "    #use minibatches where possible\n",
    "    if fam=='nb2':\n",
    "        minibatch = 'none'\n",
    "    else:\n",
    "        minibatch = 'stochastic'\n",
    "    \n",
    "\n",
    "    ### Run GLM PCA in R\n",
    "    robjects.r(\"set.seed(202)\")\n",
    "    \n",
    "    res = glmpca(Y_r,50,fam=fam,minibatch=minibatch,ctl=ctl)\n",
    "\n",
    "    ### Save raw\n",
    "    with open(basepath + '_raw.pickle' ,\"wb\") as f:\n",
    "        pickle.dump(res,f)\n",
    "\n",
    "    ### Convert to python friendly shape and save\n",
    "    res_dict = {name:data for name,data in res.items()}\n",
    "    res_dict['ctl'] = {name:data for name,data in res_dict['ctl'].items()}\n",
    "    gf = {name:data for name,data in res_dict['glmpca_family'].items()}\n",
    "    res_dict.pop('glmpca_family')\n",
    "    factors = res_dict['factors'].values\n",
    "    with open(basepath + '_pyOnly.pickle',\"wb\") as f:\n",
    "        pickle.dump(res_dict,f)\n",
    "    with open(basepath + '_familyRattribute.pickle',\"wb\") as f:\n",
    "        pickle.dump(gf,f)\n",
    "    with open(basepath + '_factors.pickle',\"wb\") as f:\n",
    "        pickle.dump(factors,f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "### GLMPCA R: Poisson, higher convergence threshold\n",
    "glmpca_R(data=dataset['counts'],fam='poi',tol=10**-8,dataset_label='pbmc')\n",
    "### GLMPCA R: NB, fit shared theta\n",
    "glmpca_R(data=dataset['counts'],fam='nb',dataset_label='pbmc')\n",
    "### GLMPCA R: NB, fit theta per gene\n",
    "glmpca_R(data=dataset['counts'],fam='nb2',dataset_label='pbmc')\n",
    "### GLMPCA Python: NB, fixed theta=100\n",
    "run_glmpca(counts,fam='nb',optimize_nb_theta=False,theta=100,dataset_label='pbmc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "code_folding": []
   },
   "source": [
    "### Load GLM PCA results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "### GLMPCA R: Poisson, higher convergence threshold\n",
    "with open('glmpca_results/glmpca_pbmc_poi_tol1e-08_factors.pickle',\"rb\") as f:\n",
    "    dataset['glmpca_poi_tol10-8'] = pickle.load(f)    \n",
    "### GLMPCA R: NB, fit shared theta\n",
    "with open('glmpca_results/glmpca_pbmc_nb_tol0.0001_factors.pickle',\"rb\") as f:\n",
    "    dataset['glmpca_nb'] = pickle.load(f)    \n",
    "### GLMPCA R: NB, fit theta per gene\n",
    "with open('glmpca_results/glmpca_pbmc_nb2_tol0.0001_factors.pickle',\"rb\") as f:\n",
    "    dataset['glmpca_nb2'] = pickle.load(f)\n",
    "### GLMPCA Python: NB, fixed theta=100\n",
    "with open('glmpca_results/pbmcglmpca-py_nb_fixedTheta100_penalty1.pickle',\"rb\") as f:\n",
    "    dataset['glmpca_nb_thetaFixed100']= pickle.load(f)['factors']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute tSNEs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "code_folding": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "variance explained by top50 PCs: 21 %\n"
     ]
    }
   ],
   "source": [
    "## Prepare init\n",
    "_,sqrt_lazy_init = PCA_sklearn(dataset['sqrt_lazy'],50,42)\n",
    "## Compute \n",
    "_ = kobak_tsne(dataset['sqrt_lazy'],name='pbmc_sqrt_lazy_initsqrtLazy',do_pca=True,init=sqrt_lazy_init)\n",
    "\n",
    "_ = kobak_tsne(dataset['residuals_theta100'],name='pbmc_residuals_theta100_initsqrtLazy',do_pca=True,init=sqrt_lazy_init)\n",
    "_ = kobak_tsne(dataset['residuals_theta10'],name='pbmc_residuals_theta10_initsqrtLazy',do_pca=True,init=sqrt_lazy_init)\n",
    "\n",
    "_ = kobak_tsne(dataset['glmpca_poi_tol10-8'],name='pbmc_glmpca_poi_tol10-8_initsqrtLazy',do_pca=False,init=sqrt_lazy_init)\n",
    "_ = kobak_tsne(dataset['glmpca_nb'],name='pbmc_glmpca_nb_initsqrtLazy',do_pca=False,init=sqrt_lazy_init)\n",
    "_ = kobak_tsne(dataset['glmpca_nb2'],name='pbmc_glmpca_nb2_initsqrtLazy',do_pca=False,init=sqrt_lazy_init)\n",
    "_ = kobak_tsne(dataset['glmpca_nb_thetaFixed100'],name='glmpca_nb_thetaFixed100_initsqrtLazy',do_pca=False,init=sqrt_lazy_init)\n",
    "\n",
    "_ = kobak_tsne(dataset['normalized'],name='pbmc_normalized_initsqrtLazy',do_pca=True,init=sqrt_lazy_init)\n",
    "_ = kobak_tsne(dataset['log1p_medianscale'],name='pbmc_log1p_medianscale_initsqrtLazy',do_pca=True,init=sqrt_lazy_init)\n",
    "_ = kobak_tsne(dataset['log1p_10kscale'],name='pbmc_log1p_10kscale_initsqrtLazy',do_pca=True,init=sqrt_lazy_init)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Downsampling Experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "removing 0 cells and 1258 genes with zero counts\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWUAAAFaCAYAAADGqUzDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9eZwU9Z3///xUVXdPg+MAkpGB4VBABwwjoKIIEhHNIsZERnez3kdMTDzWNX5XgWR/azYG0HyTGK9vdL0SjbqJDMmqyCrxQECDYYSZcAiKwgyiIwLjCD3T3VWf3x/V1dNHVZ/Vc0A9H495oNXVVZ++3vWu9/F6CyklHh4eHh69A6WnF+Dh4eHh0YVnlD08PDx6EZ5R9vDw8OhFeEbZw8PDoxfhGWUPDw+PXsQha5SFEFOFEIYQYl5Pr8XDw8MjVw5JoyyEUIBfAe/09Fo8PDw88kHr6QWUiO8BfwUqenohHh4eHvnQo56yEKJKCLFYCPGaEKJdCCGFEGc67BsQQtwlhPhYCBESQrwthJhls99RwL8C/1Hi5Xt4eHi4Tk+HL44HbgeqgcYs+z4B3AI8BdwMGMBLQoipKfv9DLhHSrnf3aV6eHh4lJ6eNsrrgMFSyrHAz512EkJMAf4ZuE1KeZuU8mHgLGAncFfCfpOAU4D/KumqPTw8PEpEj8aUpZTtOe56ERABHkl4bocQ4lHgZ0KIKinlbuBrmN73LiEEmDHlqBBilJTy++6u3sPDw8N9+kqibxKwRUr5Zcr2tYAAJgK7gYeBZxMe/zWwjQxeuIeHh0dvoq8Y5Spgl8323bF/hwJIKQ8CB60HhRAh4EspZZvdQYUQ2eLORwLNUsqRea/Yw8PDowD6ilEOAp022zsSHk9DSnlVkecVFRUVIwBP39TDw8NthN3GvmKUQ0DAZntZwuN5I6UckOnxmCft1Tp7eHh0Gz1dfZEruzFDGKlY2z7uxrV4eHh4lIy+YpTXAzVCiCNStp8a+3dDN6/Hw8PDoyT0FaP8HOADrrU2CCECwNXAaiml5yl7eHgcEvR4TFkI8ePYf46L/Xu5EGI6sF9KeT+AlPKvQog/AncLIaqAD4ArgZHAVd28ZA8PD4+SIXp6cKoQwmkBO6SUoxL2KwN+ClwGDMRsy14gpVzh0jpet9k8vaKiQt2/3+vY9vDwcB3b6oseN8q9Bc8oe3h4dDN9uiSu5Egpz0zd5pXEeXh4dDd9JdHn4dFnaG3vYN2OfbS2d2Tf2cMjBc9T9vBwkfqGFhYsbcKnKEQMg4VzJ1A3ubqnl+XRh/A8ZQ8Pl2ht72DB0iY6IgbtnVE6IgYLljZ5HrNHXnhG2cPDJZr3hvApyT8pn6LQvLcgFQCPwxTPKHt4uMTwQUEihpG0LWIYDB9kq5dVUry4dt/FiynHcCiJS23r9vBwpLK8jIVzJ6TFlCvLy7I/2UW8uHbfxqtTjuHVKXu4RWt7B817QwwfFOx2g9za3sGMu1+jI9LlsZf5FFbeNrPb1+KRFa9OORNenbKHW1SWl/WYAbTi2h10GWUrru0Z5b6BF1P28DiE6E1xbY/C8Iyyh8chhBXXLvMplAc0ynxKyePaTklFL9lYGF5MOQNCiP0VFRUVXkzZo6/RXXFtp6TioZZsLNH76QkS5YtnlD16O70xqVh//enUPbjmkEk2lvAC4yX6PDxKTXcayZ72Rp2Siut3th0yycbELk3r9SxY2sT0sYNL9lo8o+zh4RLdaSR7wlik4pRUnDii4pBJNvZENYuX6IshhHg99Q+vecQjR0qle+GULHO7pbuQpJxTUnF8VUW3JxtLxfBBQcJ6915gPE/Zw8MFSuFRZfK83Sx9K8bDr5tczfSxg9NCNk7b86G7QkGZzrNq2x70hPfZp4qSX2A8oxzDax7xKAa364OzhSfcaunOJwziZLycmmVyaaJxOmZ3hIJa2zt48q0dPLxyO341/TzWexNN+FgVATVV5azbsa9kFwvPKHt4uIDbuhe5eN5ueKO5evilMJJOx+yOeHl9Qwvz6xvpjJrVZ53R9PPYvTdIuOD+NQS00l0sPKPs4eESbhhJi1w972JbunM5TymMpN0x59c3xt+/UibXrHNbBjmRxPPYvTedugRkPM5ciuSql+jz8HCRyvIyTho5sOgfaXd15qWeJ6AJvnfGsUn7lEIn2u6YnVHJk2/tKCoUlEvC0u7cduex3puAJghqCqoCAa30etmep+zh0QPkksRy0/POhHUeK776+OqPePjN7fFb81LoaZhVDXra9odXbufyqSMLCgXlGmKxez1gGly780R1iR5rstO7odTP6+jLgNfR51EKUo3HvNk1TKge0CNdeRbZJD/tDF6xF4xfvvwe9776ftK28oDGE9dM4aSRA3Ouvmht76CxpY0bft8Qjw2nrj+V+oYWbntuQzyJpyrw7+eN56ppxyQd94y7XrUNc/QPqOiGLDam7HX0eXj0NHax1Due3+TWjzzvtVhGzymO29jSBrRxZNBH/fWnEwqbnuGqbXuYcfdrBSX+rPPOnjCEh1ZuTzKkqeGDXL1jBZF0HGv9Vnw41cBPHzsYVVGIxjxf3YDFy7cwp7Yqfs7mvSFUoQDJHn3Qp/CjOeM5e3ylV33h4dHXaWxpQxHpDtKBTvOHX0ziKJ+6XjtvPbVJIhSJct3v/oYecxQ1RXD3RbUMHxQsOPGXet4LJw+j/t1dBVWsJF7g7LAMvFlp0YSqCHRDsqhuAiOP6o9fVZIMeWoycfigILpMP7YhKZlBBs8oe3i4jmUcg34l7lkCPPnWDh564wPCunPIsNAqg3xK1uy89Ttf3ERiJNPKZyUuNWpIbnuukd9cflJB1RGt7R3Mr2+iM9p13vp3dyV54Pm8btuSNaCfX8WQ5l0HwO1LGokkvJDblzTy5xunpVdWRHWC/q5EXmV5GYvqapPCHD5VsKjOax7pFrwZfR5uYBlHaUg6dUlAE0RjBiGDLY4T1nX2HQzT2t6R1cBZXjGQl+dqZ8yiNs6mogjTLUwgakjaQpG8En/WWpc1fmwbYrAMslXFkKvBs0vYBTTBvRdPora6gsryMlZs/jTJIANEdMnH+zviyUTrs1IUQd2Da5Ji/FYS1AzjED9uKfGMsoeHS9jdTtsliSz6+VQunjKc36/daRqnSBRDwi3Prs/o7aZ6xaklbBZOnqtT9UEiUQOEYb/2iqAv5+oIa62qEBwIp1dbhHWDppb9XPrI23nHp50ads4ed3TW54JZdVJTVc4F968BZPxzs4vxnz2u+xKwnlGO4bVZexSL0+20EwaS684czXVnjk6qHojoUcDe27ULPfzmjfdJtXcdESPpVjyRVGMW1nV0Q6Z5ywJINcua0uUt1lSVs35nGxNHVDC+Kv1nki3mC3DpqSNYvHxLwY0p2coGa6sr0JTkOwHrNQCEwgYBTUmLp+ca4y+FPodnlD08XCIXD9QitSZ2YL9Q1sQT2Bt+TVGRGEm36QFVEAobjkYj1Zg99daOtPK0/gGNfzq5mt+9tQNVFUgpWVRX61gil+rdZrtIBTTBCcOOTEt85htXT6zSSH29leVl3H3Ricyvb0QVCro04q8Bsn9mmSo4Eu8CIrrBgjnjkkrqCsUzyh4eLpHogeq6gZ2D6FcF3//aaC6bOjLJ6OTaoGG3ny6NmGHrMsoSWPPBZzzwyAeOhjPRmF02dSQPvPZ+Utw7FIly4cnVTB0zGOjykHNtu3YyeFZooG7SMBbUp7c7hyJRV9XunLxpy8jeMHM097yyzTbmn1jBkVqtYnn4Fnc8vwmgaMPsGWUPDxepm1zNF6EIC5dtwacaRHSrNVdy3Yx0Y2yRi6CRZUQsg2B5aPPPHceRQV9S0koI+MXL2wByDgsoikBPsEy6Ad+6fzV+TUmKr+ajTfHdM45NUmGzkmhBv0Ldg2tsY+7CpmQwG7mo6tkJLFnvlx0+VTBvdk38WInHXrhsM6qavs6Fy5JrnQvBM8oeHi7S2t7B4uVbkmKUUkr+dOM027hrInYenWWIm1r2s3j5lrjBPveEIbzY9Ak+TWHx8i0snDuB+utPjyetsontpNK8N0SZpsbj2WB62xFdEtGT46uZvHq79YLk6mmjki5I63bscwxtlGlq3mWB+Vwocol1A6iKYPHyLew9EE4/dkqoKb5dE0ULJ3lG2cPDReyMQ0BTCYVzizUnenSWNyeAUMyAWMdduv5jgHiCb8HSJu6/ZLJt0soirDuXreUSD1cV0+CcNHKgrVe/atuetEoLa70Pv7mdy6aOzOl8hehJ5KPPkWtC1jLaD638gNSOaF1Kbp41ll+8sjVpe2dEL1oLw1OJ8/BwkeGDgoQi0aRthcRIE725UBaPDoirnmUyrN+bcayjB5eoFtfPp9ruE4nK+Ouom1zNyttm8sQ1U1h520xqqsqZt8Rcr13pW6qaWmV5GeeeMCR5H1WkqeHlOqYqH1W9fBKyAJqi8L0Zx6Yd+9tThqMpyca6kNBL2vmKPoKHh0cSIiXpZvdDzVZKZeou5P4DjxgGtdXmbLxE8XaLgCY4d8KQjBMzEhslUsV9AM6bMMR24kh9QwvzljRm7FRM9dKfWP1h3Nu3EMAjV55M0KfR2t4R97xzrV9ODf8A8dcLJG23Yt1R3UhK8Nn0y3AgrHOwM5rWebhuxz6CPpX2zq6LcCGhl1Q8o+zh4SJ2sdnUH2ouKnHDBwWJOIQh/KrCRSfZa0akynBaCba6ScOoe3BNRgNnXShqqyuYf25NvJrA4qWNnzA/pdPQ8ugzGWQwvXQwjWTQr7Bw2ea0fSK65OrH36FMUwnrOoY0t+VTv5x4obDe41AkihCCMk1N+m8pDfOCmdBfrgj45olV/Gn97qTjPrr6I576604W1U3gpJEDAfdHgFl4RtnDwwUS9S4y/VCdVOKCfgVdl/Fa18ryMhbMGZdmGP2q4E83ns74qgr+9ZzjHGfm3fr147l86sj4muoeXJOxhM2uS7B/QI03UYB9ze6+g+lJMDv8muCMu15FFQoRQ7fVrOxKLEZtHk1P3Dndbdi9x5B43OT/TiRqwAuNn9ievzNqlGROYiqeUfbwKJJUg1Y3yVn5zCnJZCUC73h+E+0dUU4fM5g5tVWAWWbl00yFs3mza+JNIdmkLa3H7SodFARvbtvDqKP6E/QraUbMLrkVMcwuwV+8/F7cC7c82kz4FRLqgNPjzbmSqAuSKbSRb2dlKtEML8iSMx3YL5SkjeFmV58nch/DQZBoekVFheqJ3Hs44SQO76R8Zre/HVaDxbzZNVQP6gdAy96DSWVxmWKsm3a3xVugBx8RcDynTzXj36qiJD1eHtC4atoo/uvN7UkXmyUNLWnxalWYNc5+TSGqG/GwQ/wcCraNNLng1xQCanIIIjG0YZEoaJ/re1wImgKKEKhCEJUGP5ozvphmEU/k3sPDbZzqY0NhIx57tLBr/rCrVIAu7QVLHCcS1TGkIGpkj7EuqG/k6bXN8f+/ZMpwxwSgZdhS49edUZ1zJwxJC4HY1T/rEnzCLAe7edZxDBsYTPJiL50ygkdXf5T2vFOPGUjDjv1EHDxTnyp48NLJAGm6IGn7poQ2LpkygsdszpmIKnJT7ktENyCKOTwV3OviS8QzyjE8QSKPQrBL9oR1PSmGbNf8YXnA1/72b1nP0RXXTbYgds0Rm3a3JRlkgKfXNptt1JeexA+eXJfRCIJlqAVzH1jDojrTG8/U7AHQETPWv3hlK3ecP56Vt81Mqnb47VsfpQkerW9u4883TePXK7byvxtbbY9bW11B8950XZBUUtuhdYckaSqaItLCFXYVGBZ2m3+2bHPRXXxJ53flKB4ehylWsidxyLEhYdW2PdQ3tDDj7te44tG/csfzm+iIGLR3RumIGCxctpnmvQeKOrddpn/9zjbbfX++/D3KfIqjQQYwDBn3nMO6QWfU4LbnNrBi86e2CUwn7nzRrKywpnpXlpdx/Zlj0vbzq+YdxTFH2cuWj/lKfyrLy2wvfJpihjb6+dW4uBN0tUPnErnQpX38WBHmsXNFizXVuIVnlD08isSa92YR0SXz6xuZV+/cTBHWJXe+kF4WBqbH2t9v38AR30cRtjrKE0fY39i99t5nXP34O9jINcSxu5WPGnDT0w3UPbiGuknDki4+irA3IFFD8mLjbp7+60427TYvEpdNHRnTAOnCuqOYftxg2/V88NmBeELTagzp71fxq4Lza4cikDHX1Vy4FUoqlqBP47oZxxLQcqsTNySuTrT2jLKHR5FYt9eJGIYknOF2G+yNoKbA4gtr+fXFk7jj/PE42WYh4PHVHzHj7teob2iJd74NPiLAOeMqbZ8T0WXeMVQwW7w7IgZLGlpQEoxepqqLnzy/iQVLm5jz61UsqG+MjVayv6OYPuYrnDC0PO0YqhAsbdhFa3sHdZOrmTe7hoguUVXB0vUf0xmVHIzodEbNi2A+3nwmwrrO5VNH8ubtZ3HTWWMIaF2dfJdMGZ50YSvFeCiv+iIDQoj9FRUVFV71hUcmisn2B1SBgZnRj+oGimJWGFjVFUcGfVz/VIOjngWYCStVVdBig0Hnn1vDwpe22F4U/LHzaQp0RPL77ffzqSDgYILn39+vEopkL4t76topHHd0OWfc9VpSbNiqVPl4fwc/eGpd2ugmYmv98Xnj06QyU7nprDF8/mVnWkw9G5pi9o90DYiFuy86MV7ZkloP3dre4dZ4KFtX3DPKGfCMskeuJNYqd0Z1FEXkZKStcjJNCEIpRjSgCRbW1cY0h3M3+D5V8KM541j00mbbaomgXyESNVCVzMmzVMzbeZFmVOfNrmHRS5tRhUKnrmN3/dAUwfVnjubx1R8ltSUHVIEhLdU1HSEc5gUCwZRmllT8qkAIkfU1KYCqirgk6bzZNfz0hU1JdxGaAmvmz8p5TmKBhtkzyvniGWWPVDL9EBO7+qwOumIJxGp/nYyVE49ceTK11RXxdmtNcS6/s/CpZsigw8aQBzSFRXVdybTUWmnrtYciUS57ZK3Da0k36qnYjaCySB3rlIqdJ2+/DgUpDb41cRhXTz+GzbvbufUPG9L2+8U/nRhvrvl4vymIZHnG+UwPz4BnlPPFM8oeieT6Q2xt70gyhhHd4LwJVby08ZOYJx2N1xzniirApyk5G/pHrjw5PkC0tb2DFZta+ekLmwhFugxWQBVIBAGtq0xv4bLNaToWZT6F+y+ZnHS8TBemBfVNrNicXuKW2pDSEYnm1VQybfRRrP7gc8fH/SpIhG0IxAkldhWwW4ZPFSiQJIKvAN8+pZr6dz9Ou2OwmlfywGseOdwoxVDH3nzeUpLLCCTLGD/0xgdoiunh6oYg4FN4aeMnceGhfQfD3PLs+qTb+GzoEvQcLVjiYFAwy/be3NqaZJDB9Egfu9pUZbOqB3buPZjWdGEYkqED0tXhUkm8aPkUM5GXaB9TG1IaW/bxk+ftK1DsaGjeZ7u9f0ClM6IjEfHKhYAmMk4Sj7+2DLvYGXcDeOadlrTt+c4VzIRnlA9RXLq96jPnLTXZJlvUN7QkdcyFY9M6kJJoLA66ePkWVt42k+GDgl2Pp+BXRVbFtUR8sYGm5nPVtMGgAPf9ZRsvbfzU5tmCa3/7NxbOncCOzw/EPzdVmHKjijBL9xRFUPfgmoyfpd1FS1MgoJpZtM7YceY+sJrrZoxmUH+/Y0mgU6dd4qCA1Db0VKlRKeHiKdU8szbdgGZDVUwXNp9wkRvqcBaeUT4EyXWw5aFy3u4g2wikBUvTB4CmoiBobGnj7HFHc92M0WnTo7WYEcyHiC7jlQ1Oimn3rNhq+9ywboAO8+sbsWK91ufmVyQydndthUzm1zdSU1VOKGwKE4XCBqFIlJ2fhwj4RNpFK+jTuONbJzB/SRMg48dJfd2pZHsL+vtVfjRnPBNHVBAKG+w7GCZVelpRBLPGDaG+4eO8kpl+VfDLb5/IjU+vz2n/fn4VQ0pX1OEsPKPcQ5TyFj+feWWHwnm7g0wyjdlakC0ORnS+/+Tf+PF547ls6kgeWrk9yWA42XQF+5inxaqtezju6HJOGjkwXq9sfa+a94ZQFXBwzAFQhYJMSa9pihqLeXdt74xKzr9vNZow46ypSbmUIRxEDIOKoA+fJsiSe8sLXUrCUT2uD90ZjaYdvyNiMHSAWRttp/lhh19TWFw3gaqKflmTimAa8HsvnlRsWVwanlHuJhKNcL4TFfKlVOLbvfW83YWTTGM+44WiRpeIzaK6CWkz7exQVTAyGLXfvb2Dx1ZvZ+LwAaxv+YKA2vW9qqkqty1RS6Qzqqd5p+GobnuR0A2JvRJHeny2btIwWvYezFjGli/WhGmrZtnpQhhQBaGwEf/M7l6+hefW7cp4XBF7RUG/ktUg+1TB4gtr49ocgGuG2au+iFFK6c7EOGs22UG3KHVs18nTP1RjytmwXrehG+QyI9WnClbPOwsgVhmx0XEWn09Ann0egBmb1VQFadivqX9AtZXatMhUnpYLudYNWyix9RqG4WgU/SosurCWf//Txoylb9ZvCoiXKV5w/+qs4SG/Klh04QQW1P89bd0+xQytfPvk4Vx++kheavokqcLGGlCQB15JXCZKZZRz6fYqD2g8cc2UNKnHYilViCSb4T0Uqy9ywXrd//XmByz/u11irYugT+Gpa0/jpJEDeWL1h2kTRlLxKSKjmFC+zJ04lMumjiqoEiRXVCUmOpRTE43g26dUs6ShBaSgM4N7bxfOUYTZgOJPuEv4IhRh4bItqCroumTSiAH89UP7Co5EcrkYOSUj7zg/L31lryQuE6WS7sxlCkKht/jZjF+2yRSFkEsyrxTn7QtYr/n19z7Lum842pUkXLx8S8Z9A5pAynRTYSc7mSsvbfyE+eeNy3uycz7oBoRyPLYiJEsadsViv5lfk90RFQFLb+gaLLCscXfXhS7mUOdikMl69tghHXZauGxL0TKeniBRiXGSHUwUOSkkc2vJQl712Nq4KE13YKfElTo+/nAmH6WyPV920tjSRja7ZUhYMKcmSSXtjvPHc/dFtY6CRdnQdYMn39oBwMK5E/D1sCXwqSpqaqYwD4I+LWmwwE9fzHznYXHSiAH4M0nn5YlPK17G0/OUS4xT1r6YuV49WXrWXcm8vhgCaW3vYN/BcMZbbwtdwvn3ror/dyaklMyprWJObVXaezJ97GDHDjqLoC89hBAx4L5X3+fB1z7gmydWZU0GlhpdGiALN46J38HGlracX8+7zfu584ITuON/NmcUfcoV3ZDeNOu+gFPWvlBj05OlZ6Wa4JtITyULi7kQWGuWhsy51jjXkmQ1Nqzz7HFHp62rsryMR648hU2723hs1UcsbWhJOm5AExklRHUpWbr+49wWkie5hFcS63wBbl+ygUgBxRrzZtcAsG7HPl7Y4FxlkYohYcHSjZx27CDe3bnPtnROFWYzit27KDBrost8priRN826D+FmnLWnS89KMcHXoqfuAoq5ECSuuRR0Rg1u+H1DfDST3fk/3t/B/6zflWboc6nPdYP48FRViWt9PN+Y2dj7VcFPL/gqFUFfvNa3pqqc8+9bjZ5HrDzoU/j8QJgZd7+Wc/VLKm9v38tT107hr9v3xisqwlGdK6aO4sKTqx0rNzTV7Hz8zrRjuGzqSK/N+nDmu2ccGx/zXgpvNRulSuY17w0hU36Q0pAlvQso9kJQ7Ej7XOiMGsyvb7TV2nh45XYiUaOEZ8/O18cfzU8u+GqSSl62SrivDjuSH8XqtBNLyn5+US3zljTmfMcRihj8v9c/KDjpabFx1xeceXwl504YkqYKt2DOONsKGauU8OE3t3PZ1JFFnd/CM8p9jESPDiRXTxvl2hXaDYqNBQf9SpIqF5jdY0F/6TJRxYaDMmlZZCOficqdUcmTb+3g1q8fn6a14TbZughTeWnjp9x09liGDwqyYlMrRg4GsiFlnqBl9ObUVnHr14/nnhVbcyqnA/tZe/ny8/99j6BP5WA4CkLgU8zjXjl1FMMH9cv4nnS10Hue8mGFnUfn5hW6WNyIBYfCBmW+ZInKMp+SJEbjNsWGgyrLy2y1LHLhn06u5g/rWnJOTD30xgeMGtw/JnxfutCEVoAM5r/9YQNbP/0SIWTBbdX/+fwmFr20BUWQs0F2i6ghu+q1pYx/Jo+mqObZcTCiZwwx5YNXEtdLsTQMWts74ttKVY5md65M2+0eT7xgWBObFyxtcnyuE06GsJTx8sTBnIWWKV42dWRBZWXPvJO7QQZTsCjfSSSFENA0Lj8tv4v9xt3tRIzCDTKYnmhn1MhqkIspnysVndHCvvOpeJ5yL8TJ4yxFgs/pXNm83tTHv3fGsa5UhLhZ3ZFPKKXY5GVleRkXTBrGHzPoK7hFqQ0ymNKj/QNazrrE3U4v7UR2owrKa7POQE9MHrFry07UxnCzXMzpXPXXn542zihxDXbPc5rfVqimR7Gx6Z4oq9u0u405v15V0nM4oSowY+xXeOO9z3o04Xe4k+d33muz7gtkSzoV6tHZGTlT1jH5e+FTFNbvbMu4Brs1+lU1PurHyqbPm11TsMdQTHVHT5XVPRXrkOsJdANey6G9uzdRrOBRb8OnCq9O+VAklxBFvgbLyWtsatmfJqsYMQwmjqjIuAanNV4+dSRH9fezcNkWfJrC4uVbODLoc81DzdV77onmmlXvf5b3aPt8OdSM2KH0WgAev/oUpo/5StHH8RJ9vQw3kk6JOCXgNu1usxXCmTe7hvFVFRnX4LRGMEcehXWDA516wck+O/LR+uju5pr6hhaucJjg7AbWzYyLEg0eLuNTBUGfOz6u5yn3QuomV1NTVc76nW1MHFHB+KrCheqcvEa7EEV/v8qE6gHxNWQKk9g9bjeBww0PNd9wRHe0gieubX59Y0njuFYJbiH5tvKAoL3zUPNJeyfejL5DGDeTVE5eo12IQpcyrzBJ6uPZ5tgVmrgrJByRS+zdDdGjxpa23loIAOAZ5G7iH04Y4tpF3wtf9DLcqvdN5LtnHJsmFZotRJHvmtftMLVq7Y65atueomRG7Yx9Z1RP6/JLrauuLC+LSzkm1lOv27GPJ1Z/WLT0aX1DCzf8fl3ew049Dj2OH3KEa8fySuIy0BMlcet27OOqx9YmTYKwm0ySiyeyK9sAACAASURBVJeXOobquhmj01qy8/EW7fa18+oTPVQgY4lfriSqsHXqkrJYp0a2uurE7aFIFCFM0ZzUmXj5rimXiTIehw8Fln96JXGZcBgH5d7lL0eCfoWOaHpFhGXgEkVoEsWIUsMbdnHYh1Z+kNaSnWslh5PxtYv1rrxtZpKHmin0YGfo7bZZcfZv3b8a6Bp7v2BpEzVV5bbrsNsOkoiNTkViZ2ShFR4ehy9C4Fp1j2eUexH1DS3c9tyGJHWtxNrHVBEaq1HDLullZzQSBW3ywSnRdv8lk7PGeocPCqZdZDqiOsMHBXli9Yex8jkR16K1jm0XT1/e9EmaFkOmuuo3t+5BEbmVLEQMg6aW/Vz6yNs5xfJLOUbpUOHWr4/lly9vO+RK3+wIhc3vjxtzNr2Ycgwp5Zmpf8CX3XV+K4uf2kErMKdLWIbRruXVTv/CSbns4ZXb845PO2luADmVnqWGyKSUPLZqO3c8vympfG5+fSPz6+3j6a3tHTy08oO0tYV1+6RlKBLlFy+/l3HicSI3nDkmPrY+l1h+ZXkZ82bX4OsmDQZVwHGV/bvlXG7xq1cOD4Ns8dMXNrlS/ukZ5V5C894Qqkj/ODTVNLiZZr+FdT1e4WAltCzlslT8aroBt563aXcbKzZ/yorNnyZ9uZyqKmqrsycLm/eG0uo3FeA3b3yYtjYhhG2HofX6/Wr6QLrvzTg2LWkZ0MzR9qkJOFWY7cip9PerHHVEIC+xp/qGFhYv34IQ+ZmdQmuNdQlbWw8U9uQewsXh230CXZrVOMXihS96CcMHBc05ZSkkzvxyul02JNzzylbq392VdOt92dSRPLRye5IeRaonm5pAs/CpgrsurKVucjWV5WWce8KQpLFBdZOGUVlexvSxg7n/kslAlyA4dMWFg34lbd1OKpy6LpEpRjNxvanHCWiCy2Mx8sQSuH0Hw9zy7HoieleytJ9P5d5LJjF0QBkX3L8maR6bLmXWLsZEipk04hVqHNq0hSJFH8PzlEtINunLxP2a94aYf+44tIRPxKcKFtWZcdbmvSHmze6aaJxIRJc8vbY57dYbYFGdsyebaFxSheUjumR+vXn7/sTqD9PmuNW/uyteVnbLs+u58ekGVm3bYz6W0H1X9+Aa6iYNi6/Br4qYeFE6/zJrLN+bkV6+ZyUjU73yRXW1SV65VQJXW51uYKOGwdABZYyvqmDxhenvST4lgvlMrPY4vKgI+oo+hlcSl4FiSuJSqxXmza5hQvWAtKy+3X7Vg/oBpue5atuetMf9msrPlm1K061IpL9f5XffOZWTRg50LHuzK79LpJ9f5d6LJ3H9U+m1uEG/gq6T5HFmUpirv/50QmEjPioo1cusHXYkW1u/zFi+B7mX8GUroXM6Trbjt7Z30NjSxg2/b+gWCU2PvoOmCNbMPyufCgxb78S73JcAuwaQO57fxGWPvs30xa/yxOoPHfdbvHwLtdUVnD3uaADbxyeOqMg6WPJAWKepxbyYWB5k6pclWwWBdQ6fTSA2HDHwac4Kc6nbQ2EzFBAKG9xw5pi04zXu+iL+Ojujkoff3J70eGKDit1rSaVucjX115+OjH3vOyJGUvLO6T1x2g5ddwC3PLse3TCS7mo8PI4/+ghXSuK8r1UJaGxpsy3FCoUNwrrkjuc38cTqDx2TSIn1sk4GLvVWe+7EoWnH+dmyTWzabZ94sDxCKyTiT8lAqQrMP7eGoQPKbEfUSwkdKZUNoUjUMTbb1LI/HtK479VtBLJYtMQkWyYxokwholDYSDtPoZNazOqYrgtk1DCnX3xn2ih8nlKQB+bkFaffWz54iT6XyXWg5cJlW/jlt2vTbuM7Ika8fTjoV9Juka3k00kjBzJ97OCkbO/Lmz9NCmlEdDj/vtX8/KLajFNDzj1hCC82fULQL+gMGwhFUOZTuPPFTQhhLxgZ9Kt0RI20CRAf7+9g3uwaFi/fkhRyscrNEps4MhExzPdhxeZPmR8bf5QqRrRq2x7m1zehKmadszUfLVOSsRC1uNb2Du5eviXts/CrKnNqh3LhydV8877Vrgzv9OjbrN/ZVpSAGHhG2VUy1RKnI7n5mfVpWwOqIBQ24oZTxIyXFROdN7sm7uklxpvDum4b0tANM2FnNZds2t3GvCVNhPUuIxdP4ln23JAJxt3+tUR0c8Bp4kUgasC/PPMuhpTxGHrQb4Y0cpmppimCoE8lYhjUTRpG3YNrUFKmmYDp7Ta2tHH7ksakZpLblzTyRSjCope2xA31hZOHpVWl5Nve7XSR7Yjpb4yvquDH59mPoPc4vJg4ojiDDJ5RdhW7Lrp+PpVvTazimXeSBW+cRGwMTA85teTKMCQ3zRoTNzhR3UBiVklY53O6i1YVQfPeEKu27WHeksaiBHT6+dW40bXTY7aaNRYv3xLfRxUiTWvCDgH86p8nMnRAmW0y0CJiGLSFwmndfRFd8tMXNyddnP64roU/3zgtHtPOZJATk3xghqHmZ5gaLRDUPbiGhXMnMKe2ijtf3Ox5y4cxcycOLdpLBs8ou4pdS3FnVGfW+CEMHRjknle2Za1TlVKyefcXaTFpn6pwz4r3Myb4ynwqnVE9rStQN2Tc0BdjkAOa4N6LJ8XrkY8M+liwtAkFwcFI8utWhWDhsi1J1RnZ8PsUBvbzEwobtroS1gVh4dwJHOlQepT6/kR0yebd7Yw6KnM3nJ1wkU9JDx8lEtYN0LtazoM+1bGSxePQZ9nfP+GMhpaiJ+14RtllUksMdQk3/L4BkGZ8NksJoioE855rIpJiXMzwROZz61Ly4/PGc+eLm+KG2ap1djJ0+XDh5Op4VQh0NWy8GfPAEz3XzqiOX1PzGjefqVHGr5J0QWht70BTSLoAKWD76uYtaaRMUzMq5eUqXGSHKgQftH5JZx4XII9Dj86o4cosSK/6wkWa94bw2VQVdEYNOqMyp1vbTl2mGeSApnCRw9VXFSKp2eGqacewZv4sHrnyZB658mRWzzuLusnVOQno+FVB0Of8lXh2bXNaC/aqbXu4/bkNaaEEJFl1J1RhxZEVAlpyo0jdpGFJ++qG4ItQJGkk1d0XnUhAE/TzqQQ0wbdsKlDA9JatUrt7X32fM+5KruDIpRmkn988x3emjUqr6DgQ1vnFy+/ZVql4HF4oQhRU3ZOI5ym7SFPLfkJOPcQ5EtCSb5mt9mAgLS4NIISMx2FDYSNeg3v2uPQr9XfPODYu+dkZjWJIkXShCOsSLcN1wwB+8NQ6VEXEpTvtRJQgeXRRQFOQ0kg7nxDm/DmBINHHbW3vYEmK6LwukxOWkNxabTWl5ILl0dRUlccbWjKFWfwq3DxrLGccN5jxVRWcMKwiloSFUCzu7QndewB0RoqfBekZZZdobe+wTXwlYiXinH6/ZT5BqjNrIKmtNpMHqkJaCCPgU9nQvJ8bn97uKDmZGC8FyWnHDmLltj0oNlobIAloCgLMkrcUIrokost4HNUUUcrsEYejBj5NQaa8OCkFYUPG1eysWz+zzC+9FM9KWKa2VjvNB8yErku+df9qyjSVUCTqGBpSBUgE97/6Pr9asTX+3n4RinDni5tzOpfH4UTxF2cvfOESudwCC0Gs7tceXYcFc2ps9Rcqy8v49/PGpz0nqhs8vHK7o+RkatdgZ1TyyuZWOqOGrTBQ0KfxwKWT+VndBDJVsVmv1U5EKRWJaZhTL0Z6Snzdpyg8+dYObvj9OtsEW1Q32HcwbNsoYkqV5n6XEjFkPKwRNex/SkFNQVFEfL/USeBepYVHKv38WtHhC88ou4RdzDb1zY0aZPwhS+DIoI+Vt83kiWumsPK2mXGPt7W9gwnVA7j162Pxqwr9AyplPoXrZozGrzp3reUrntMR1WnZe5AfLW3CQTcI6JLuXFRX61q7cVjXeXjldscSNN2AW55dbztTb9W2Pegui85HDCOtxdypldzDA+BgOOqFL3oLqWPtw7puGuA87ETUMMMC9defzr6DYfYdDFNbTVrn2oI5XeJGAA+tTNaJSBwqmu+EDF2X/GzZ5vQJH6rpMSYK+1SWl8XjulZn4ebdbfzq5W05vWxNAVVR0BRBOKpzzrijeW3rZzhVlelSxkvOErPc1t2A23k2VRFpycqwrhPwZY5BexzOFN9y76nEZaAQlTirAWHfwTA3P/NuWtNEaqTUMrQWAdVMhlk2UVPM2Gvirb5PFayed1Z8RNStf9iQdExNEWixMVJ2Q0XrJpldbnb1xXb4BPzkgglMHFFh24Rhveamlv3xZpGOiB6/9U8loCkIYRr2L0IRfvrCprx1hhOHyWZTu8uGKkBRBGWaapsABegfUOkI6wghCPgUOiO64xBWj8KwLvx9mTKfwu+vPS3XsVDe4NTuwIr/trZ3pMVMwTTCiT/41GaHVF1j0/tLb4hobGmjttqswU09S9Qwy+/m1zdyZNDH9LGDqb/+dNbvbGPiiAoGHxHgrHFH0xaKsCCmK5GJiIS1H+6h8shAkpA9dCUR7br2dJsfWNCnmK3f545j+tjBTF/8qq1BLvMJOiLOP9Cw3pXlLmReXj+/im4Y8bplMDv4Pmj9kntf3UY0oX28v19l5vGVvNC4G6SMPxbQ4NcXT2Lz7rbDZhZdKcnXIKtC2P7GepJoQq19oXhGuUQkhjNUIYjoBpefNpI//K0lyaNLbXiwq7Bw4qm3dmQsxeqMSm56uoFw1EAIU1fC6lYr00yNCUsbIpvX/Kf1u/nT+t1JE0ks5bR8dIWtErLFy7fQvPeg7fp9CkR1s27a6fV9b8axSTXLC+dOyEkICtI7EyHzxSWi6/zvxk/SjqOpChua9/Pwyg/QhHnx8ug+Th41gPXNbYSjRq+5IN48a0zR8p1e+CIDxYjcW6TqKcy4+7WCxgglogp4/l+mM/eB1TmKH2U+1r+ePZYjynz8LKETMBMBTeHN22fy1Fs7uPfV9ws6b3+/SjhqpDXK5EJAE7x5+1nxO5LE93dBfRMrNrfG9z1nXCXTxgxOUq1LLRlsbe/gjLtec7y4qIrplaVeIHyKwJDSG/HUg8w+4Whe2fxpzo5MKZk7cSi/+udJ+TzFC1/0BNZV0zIc2Tw6v2rGmK0fuk8V/ONJ1fz3O83xbYoiWN70CX7VjIEWgy7hF69sA3JPUSgCljbs4jdvpE+XzpWwbuD3KUQ6UzUzwKcpSReugCqQCAJastKb3dSWVe/vSTrem+/v4Wd1pmCQlYy06r4tnnprR0Zv35S4SP+8LjppmG1Dj0f3sXzjpz29BABuPec4bpo11pVjeZ5yBgr1lBO9t1R5zetmjOazLzt4Zq3zj/nWc45j3NAjgS4DkurJBTTTUPVUa69fE4SL8NKvmTaKp9fuTDK+PhUev3oK1/72b47jpKwko51361cFPi1ZTtRKCO74/ECSAU8cCzV98at5d+R9f8axtOw7yAtN6WENj8MPTRHcnaJbngOep9wdJHpvYV3HkMnymrnc7j/w+vusvG1m3MtesfnTdD1iiet1uflQjEH2q4Lvnzmar8balRMN5fQxX0kqLbS2p0oi2nm3qki/SIV1PUkKNVUoP1tc3mn9j676sKDQi8ehSdRIlwEoFM8oxxBCvG6z+Yh8jmGvNpY/VvOHdYtuF+5IrdIoNUFN4aKTh/PcuuZ4sq5Qrvva6KQa59RBpXbbU2PHD61MD52EooZZ3ibAspeGhOVNn6S1YKtCFByC8XQuPOxQBGkyAIXgGWUXsRO5L4SwrjN8UDBhSkj3GAFVmDHmgCbSLgJh3aB/QCnaIANJHYiJMffE/0/cvqxxd5K4/7cmDkVTFTqj6dUiqW9VRJcxA558p3EgrPOrFVvzem/7B1Qi0a7hAh7pJF4QDzfcKIcDL6ackXxjyq3tHa5VV0wbcxSrP/i827LKfhUWXXhi/P//T0pDSiFoCvzjSdW2ybArThvBtyZVxxtOUmO9mUrU8iWoKcwafzSvbPoUTcn/eAFNYf65NVQP6sfmj7/gF69sLWo9hzKqMDVeDkcVU7eqLzyjnIFCEn1OMWW3UWPiRqqSPsMun2P082vxLr8lDbtQYy3Pbvyogn6FK08bxW/f+sjWw7Y880QCmmDpDdOY+8Cagl+XE4qAOV8dwl+2fEYopSa7TBMYEiYMq2Ddzq7PWwH+v/PHc2TQx23PbSiZsbF7L3ozIwcFad4XOmy9YjsSO21zxB2jLISYCtwIjAWOsjmwlFKOzuugvZRiqy+CfoWXmj7h4ZXbMQwDF+78bSn0ljGgCR649CTKfApXPvZOxlFThRL0K3lrTJ9YfSQbWr4o6Hx+TeGiycN4bl1zXlNPgLRJJvFjxjRXuyOMZD87vGfwKTCgn4/PvozEtw0+wsc9/zyJ7zzxN9cvmqWkuy56j1x5ctJ0niwUX30hhLgCeByIAFuBnfk8/3ChsrwsJiLUiJDZk3LFfmEMWZhugCIEG5r388Cr7xcUBc/lYlCI6H+hBhnMcsKAppCDpEcamWyMKbla+l91bzHIYH4nOyISTcCJwyu4evoxfKN2GL98+b2sBrlME3x3xmgaduxj9Qefd9OK7QloCpedOoLfrvmIInutuoV8E30/At4DzpZSflyC9RwSbNrdxr/9cUPOhvb6mWPY095RVCPC+KryvI1ZKGJw/6vvF2QIgj6F+y6ZHGsz3g5S0qnLHr8N37X/IL97y11fQTckh6NSpyGJSwJs3N3OlGOOYtPuNv7f69krVnQJ504YwsM2VTLdhSLgxpljGNTfz6KXNpfcIGtKemNSQcfJc/+RwL95BtmZ+oYWbn8ud4MM8Mib24kWmdEr1Lss9HsaMcyJKGePO5rLp46keW+IUCTKNY//Db0bspNBv8K00YOTWqonjxjgukEG08DYiSuVgt4UvkjEGkDw0Bsf5FSf/Y0JVYTCBtEeukIrwMNXnExtdUWsyai06/CrgsUX1hZdDgf5G+UWIFD0WQ9RrDrlfGPHbpSZdTdnj6tMKl+rLC9jxeZP0VSRdyy3EKJRg2ljjmLamMGEowYnDDuSax5/p/QnLiEqcNqxg1i9fW9PLyWNzmiUh1duzzmu/tLGT5g1vrLH7pr8PoXa6oqsLfRucNoxA7n3ksmuGGTI3yj/BrhUCPErKaUnIptCrnXKAjMG3JebEF5/77P4WCZLS3nRS5tL7pFYRAz4yfPmjDyfKrjitJFm12Mffk91YE0vNMhg6n8oIvf31qcoPPN296ec/KqCEJLvnXEsWz9tzynUUixvf7iPZY27uWraMa4cL2P1hRBiRsomBVgE+IEHgA+xmZoppVzpyup6mFLWKZsjlNLF1EtNrrfHqZKiqZT5lPh0bEX0Xm8/dYiAR/eQ7ftjkcv3MZ88xSmjBtDY8kU8x9Fd+FWFVfNm5ust518SJ4QwSH/PEg9k95iUUqr5rKy3Umid8u3PbcgphNETSTG/aspNZrqj86sw47ivsGLzZxmPlcngaYI+ken2KA1udfYpwL2XTOTmZ9f3CnlOJwKawtPfzXniiEVBJXFX53MGD1O3Yde+UE5dXz5NwU/3eZk+VXDehCr+vD5znlYiGPOV8qxGOZMHKjErNHqrB91d5OoxlhpVgKYKLpxczR/+1lzyjju3bk40VeGH/73Bhcl3paUzasTnYhZLRqMspfytK2c5jFj1/mf8+i/bcto3GjVQ1O6ptVKFmURamsUgg9mB+JuV24syKFLSp5oLSoElRXrN4+/0eP5Al6AiOO7ochQXa66dGm7s980/XNdXBtQGVFFQTb4deVkEIcRjQohTMzw+RQjxWPHL6pssqG/kskfW5vzFi0qIdJPh0iW8uS2/Iv5iVmZw+ArTdCFo/aKT73+tdzS4hqMGd764ydVjThye2+26gNjFoOdx8oNUh+UpgnTp3BSEIlwRI4I8jTJwFZDpG3YMcGXBq+nDbNrdxtNrm/N+3mFvtw5hIrpkwdImZk8YQkBz745IAc46fjBBX/7HjBrutos3tuzHn8Nrk8A3T6xy7bzFoDnYV4m9YVaE2XJuR9Cv4FcF82bXuFYS5/a9c3/MFuzDjvU72xwfU4DZJ1R221pURZDlwu46vcMH6n0oCD7e38GiugmU+RTKXDDOPk3wf2bX9IoLuioEpx6Tm7f8XMOubv9e2tHpUMxrSPvEe9SADoesdThi4NMUFi/fQn2DO6PBsn5DhBAjhBAzEsrjaqz/T/m7APgBUNgkzT7O5192Oj5mAMs3tjo+7ja6Ibs1dKDEkkhBn+L6Vb6vczCi84On1gEwb3aNO+V5Ej7e38GlU0Y43op3F6GokVdYzM3vZW8w8LqEA506HRGDBUub4rX7xZBVJU4I8R/Af5D9Tltg2p+rpZRPFr2yXkCuJXFu6Sh7HNr4NcW1mYo9rTHikY41DzKPsriCVeL+BHwUO8BjwMPAWyn7SOBL4B0pZf6B1T5O896Qd/vukRU3h9zaGWTPUPcsEcNwJdmX1ShLKTcAGwCEECOBJVLKvxd95kOI4YOCRPpI6Y7HoYtnkPPDp2Db5FVoweC5JwxxJdmXV0RKSvkTzyCnU1lexs2zjuvpZXh49DiC3JK+qjC1r50qIbqDiGE2VIE58MGvCm4957iCY9UvNO12Jaacd5pACFEmhLhNCPGWEOLT2N9bsW3uFOr1QW6aNZZzT8h54kBB9HTveqAnf0EerpHPp+hTBfkUjEicvcwyn0AVcPEpw3lrwSxumjWWc8aX9jeTDWswhJTw2NWn4NeUgu84NEXEBwAXQ76TR74CvAqcAHwBbI89NA44FbhCCDFTSpm5P/cQ5doZo3l506clu43sSVk+Ad2mAOdRGk4ZOZCaqnKeWduMpjiXeSUS1aUrpXeqgPNrh3L19GMYX2UKwbe2d/Dy5k9dOLo7XPP434qq6DAkrsSU8/WUfw6MB34IVEopJ0spJwOVwK2YxvnnRa+qjxL0F36V7e0coi/rsGLdjn08+fZOoobMySC7iS7hj+t2MfeBNfF63qfe2lESkaEyX/6WNaxLwrpBR4HJWE0RLKqb0CMi9+cDj0op70ncKKUMA78SQpwAzC16VX2UUNggoAnPo/TolRRibtz+JndGDW5fsoHKIwM8VIJRUZOHV/Cjb5wQ0/feQkQ3Sl6zf/Epw7nl68f1WEefH2jI8PjfYvsclgwfFMTwBB88XEDrDZ0RJSKiwxWPri34+ZmMVkNzG8saTdEtKUtvkG895zgWuTQGyiJfT/kdYHKGx08CCn+3+zj/8sy7eY+C8ji0UYXpbeZrHLp7+EF3Y8jCcxQii3zho6s/Kui4hTB0YJB1O/YxfFCwxzzlW4GLhBA3CSHiBl0IoQkhbgbqYvscdqx6/zPe7qWjfDx6Br8q+O13pvQadbRDAZ8q+O4Z7oxdcoN5Sxq56rG1zLj7te7TvkjhF8DnwD3AZ0KIdUKIdcBnwC9jj/1SCPFqwt9fXFlpL2fV1j09vQSPXkRAU1h8YS0dEaNor1dVzPFb2UiNeAjg1GMGFtxtevEp1QU+s3R8Y0IVT6zZgd9JZ7ObieiS9s6oq9oX+RrlYzFDHjuB/cCg2N/+2DYfpnxn4t+xRa+yDzD9uME9vQSPXsK00Ufx5u0zqZvsjlFThVlHm41U2y+Bv364r+Bk3XknDuXWr4/Nup8i4JxxlXkbynzNaplP8GLTJ3REjB4fGqBAmhyrT1FcqVPOt6NvlJTymHz/il5lH2D6mK9QdWSgp5fh0QsoL9Pi8cXa6oq8mi/sCGgaZ4zNftHPJzmoAPdfMjHe0WbHs2ubufcv71OmiZhIvf1+hoQ3tn7GY1efwtnjcpeoTTWriiDje2UYpmxpIgFV4FcV+vszt1apwjp+cR62wJTGDfiVtMk6bmlfeEqLLvHE6g/Z/YWzfKdH/vhVwQUTh/b0MvJmxebW+G1sZXkZd190YlE/tM6ozhvvZZZ+PfeEo/OS8RRC8O6O/fgyPOmFxt1EdLOm2UpWOtm0sC65+rF3eGNr/n1jPgX+4/xx/PwfT+Sui07k4inD086jKnDTWWPTpU8F/OnG0/nPC76asfFDl+b6o4ZEFYLygIZfFTmFhRKRmNK4iaOf+gdUynwKC+e6U6ecVbrT9klCjALOBo4Gfi+l/EgI4QeGAJ/E6pb7PPlId05b/BciPdlydwjiVwQGfa8SoZ9f5cnvnBqXcNy0u41v3b863tLrxDFH9WP3Fx1pErAXnTSM59fvptOm00IAT147heljvsJ9f9mW08De3oZfNWf3ZfuYA5rChZOHJQ1+9amCb0yo4vnG3Xl9T37xT7WMqzqSugfXFCW5G9QU/vWc45g7eVghBtn2MlKI9sVdwDZMCc//pCtmXAZsAq7P95h9nea9Ia/lrQSEDdnnDDKAIWX8Nra+oYULcjDIAB9+fjDNQAQ0wTXTj0E61ICV+RQ27vqCnz6/kXv+4o5B7u7b57Ce21CGzqjBkoYWFKVrhRFdsnT9x3l/TzojkvFVFSyca06EKWCyFmCK/P/yla2s2uZeoj/fwanXAf8GPAB8nQRLL6X8AvgfzK6/w4qgX/Hqkz0AM3ZpzWvbtLuNeUuaCk5KBTSFRXW1jK+qYPGFJ9ruE4oYLHppC4+u/si1luXe/FVWheJKY82Io4Ks2PwpRwZ9PHLlyUXJI3RG3au8gPybR64Hlkop/1UIcZTN443AjcUvq28RChuewLgHYH4HFr20ha2ftvPHdS22HrImzEnmmejnU7n3kkmcPc5UUaubXM30sYP51StbeXZtc4/fmClAmV/lYDi/mF0gFqoo9LeiSwNk8Ub5ykfXxtdQqH5yIqowFeK6XU8ZOA54JcPjnwGHXW3Y8EHBHp+V5tF76IwaPL222dYg+1XBE9+ZktXbi+g6tdUVSdsqy8u4pQi9X7cpZN6gRKLm8GMJaAoBTTBj7FH4VBFPpi2qq2X+uTWFLDeJxI/GjQvcgbBOU0vm3FOu5GtKOjAnVjsxErNm+bCisryMxRee6I2EdQLIZQAAIABJREFU8siIpsD3vzaa1za3Zo2BRgy4xyZp19jS5uodWc2Q8oIaMXyawmWnjsj7uboEw3AOkFji909/9zTmnzuOtR/tw68qRKIG82bXUDe5mgnVA+gf6Fl1cbtXvXj5lh5pHlmLgwqcEKIMuBxYXeyi+iJ1k6t58top/OPkYfzDCZVcdNIw7r9kYlqBucfhixCCR1d9mLM2w9Nrm1n1fnKJWVvIvcImBfjlt09EKcD17owaPPXXndx01lj8edwm6gY4qWMK4Pl/mc5Ns8YyfFCQxcu30BExOBDWCesybvSGDwq6MxU8A9muNXZn75HmEUyt5KlCiCeB2ti2IUKIfwBeB6qB/1v0qopACPH/hBC7hRBfCCGahBDf6I7zLqhv5LJH1vLHhl3878ZW/rz+Y3743xswCig59Einr13a7H7TEV1yIM8Y7FWPvRPXVKhvaGFBfVNB61GwrzH+eH8H82bXUOZTCGSwRHaPdEYNHnj9fRbMMZ8fLLSEIYYE1sb0Y5r3hvApycezjF5leVm8aiKT8Sx0Uk5AEzz/L9PzaoSBHmoekVKuAH4AXASsiG1+ElgGnAh8V0qZOum6u/k1cIyU8kjgauApIUTOM78LYdPuNp5emzzEO6JLwrrMqRTKIzu9uSLAjjPGJufBM3XOZSJqSBYsbWLV+58xb0lTwcpqU44ZmFZ2ZgA3Pd3A4uVbmDe7hqe/N5Wnrp1i6/leftoILj9tZNqdn4KgPOjj/ksmc+0ZxSsqLFzW5Q1HDPuOudb2DkYe1Z8bZo5GVZzL2Qotp7xwcjWDjwiw6v3MZW6aAn5NoZ9fJaC51zySb/UFUsqHhRD/A/wjUIN5Ed0G/EFKuavoFRWJlHJLyqYAUAXsK9U51+9sK/oY1veqrxkfj3RUYWpOJFLMxVkakqsee6eomu23P7T/+oditZyLl29h5W0zqSwv47wJQ1i6/uOk/X739k7sFDMPRnRu/cMGynwKeoaavFwrHHyaWcVw0siBzJtdw8JlW/BpAt2QLJw7gVXb9rBgaROqEBnvOjQFbp51XFozjSZAKCLj51H/7i7OGnc0PkWhw+EXaTWyLGloiVWDuOd8FXS/IaX8REp5n5TyBinl9VLKXxVikIUQVUKIxUKI14QQ7UIIKYQ402HfgBDiLiHEx0KIkBDibSHELId9HxRChDD1n1cAm/NdWz5MHFGRfacseMb40GHSiIFoKZ5xQBWOFRfZfOhOvfRNNFZooLW9gxebdtvuk+k72hExMtbq53qnENXNxpv6hhYWvbQZVUA4qjNvdg3Txw5mwdKmeJzZ6Tx+Fe6+6ERumjWWuSlt+v90ynB+NGccflVxTFJaYZNUT91CFfDoVSdT/+4uOqOSgxGdzqjsMZU4tzkeuB0zFt2YZd8ngFuAp4CbMb8jLwkhpqbuKKW8HjgCOAdYIQvpJc+D8VUVXDJleNHH8QzzocG6Hfs40JlsNIQieOKaU9JCA341u8HKtcKhGDnLjqhO0K+wtGFXSUru1BwP+r0ZZgjktuc20BmVhKIGER3ufHETjS3Z70g1RSCE+R63tnfw0sZPkh7/w9+aWfTSFlQFx6aeiGFQW13BDTNH2z7ez6+xcdcXaTrZbiX6MoYvhBCvFnBMKaW09WBtWAcMllJ+LoS4AFjqsI4pwD8Dt1jzAYUQvwP+DtwFzLBZhA6sEELcLITYIqX83wJeS84srKtFInlmrTtC192FuzdeHmD/fl46ZQTHHV3O4gsnsGBpEz5FIWIYnDFmMK9sziw2lGtHYDFyllFdcv69q0rWAJWLp+9X4cThA3hz2560Co2oAc1709vQU7HCMfPrG/nhOcenXQyiBkQzlOQBnDFmMIte3JwWwrE4EI7yf19+Ly0E4laiL1tM+VjSv2P96WoQsWqSB8T+3QN8mevJpZTtOe56ERABHkl4bocQ4lHgZ0KIKiml/T2X+RrtL3ku0trewZJ1PR5SzxvPIHcPz7zTzO/X7mTh3AmsvG0mzXtDhCJRrn7snW5fi11sWFJYR6pPMWuqy3wKUkp0Q9qWvBmGRFPsy+H6B1Q6IzoSwS3PricUidqea/tnB/CrIqeLT2dU8qtX3itoane2i6QhwUhZg5uJvozhi1T9ZGAWEMKscBgqpRwkpRwEDAXuBQ7G9nGbScAWKWWqwV+L6exNBBBCHCGEuCz2ryaEuBCYCbxZgjUl8dRbO3pceNuj93IwrNMRMZi3pJE9X3ay4/MDXPP4O0S6WXApqCnkOp3Kr4qsoRVFUfCr8L0zjuXN28/i7otOtNVE1qUp1H/BxKp4xYImYNroQXx/xrGoihKf4uFUx/xcQ3Nev7FEg9w/oOJXSzOQtp9P5YFLJ7s21CAv6U4hxJ+Bg1LKix0efxYok1JekPdCusIXM6WUr6c89ndgl5TyH1K2jwc2AtdKKR8VQvQHnscc7iqA94E7pZROYZFs3YcVFRUV5CLdOX3xX8izBNWjCPpy2MWnCBDFVWQUSj4aLX5VsGDOOBYv34KgKzTgtO+fbpzG+KoKNu1u42cvbGb1B5/nsS6BnoctsjxmnyKyXtiCfoVZNUfzvxs/wShCd8OJMp8Sr1zJE1ekO88E3sjw+OuxfdwmCNgpyHckPI6U8oCU8iwp5QApZYWU8iQng+wmzXtDaErPtn0ebvRGg6wK8za2PKCZky4c9osYPVO/rink7CUD8canS6eMIGrIeCLRzhMO65Jv3reaW559l7kPrM7LIAOOBlkVpDWl9PerfHPiMPyqyGkobWfYiAv2p77tfhW+UVtVcB25TxVxVUC3yNcoS2BchsdPKGItmQhh1hunUpbweN7EjLfjH5BTAfLwQUFTvcrjsKbMp/LDc47jG7VDTI2HEp3HJ2Dm8V/Br5ljkDQlfV6cHT84cwxBX+6tCVED7nh+E4+u/ijeDGVtt9/f1DYutMHFDl2me/ZRw+D5DR8T1qWt8H8qmex2QNO4etoxrJ53FjedNYaAln20lIWqCPyawuLlW1ybZA35G+WXgR8IIa4QouulCpMrgeti+7jNbswGkFSsbfZp0m6isryMRXW1Rc9i8+jbHAjr3POXrTzzTvYfaDFfFVVTuOuiWhacW0NEl/i09Hlxduc7d8IQx9rb3kpAU7huxrGU+cw7kIAm+OaJQ9PqwO3QFMF3po2izOdsZDtjpYCV5WVcPnUkl506glCOI4R0Q3KgU3d1kjXk/934IbALeBzYJYR4QwjxRmzbY5jG8YeurCyZ9UCNEOKIlO2nxv7dUIJz5kXd5GrWzJ/Ff5w/vqeX4tGDJM5uS8QKbVj6EqkDQHMhcRYcmF14Yd2wPWdqPsvATEbPm1287GUpcDJEUsLlU0ey8raZXDVtFCBY1vRJWh14Kn5V8D83TeO6M0fbhkYsmy6EoO7BNSyob+SMu17l0dUf2U5BCfozm0q3apQhf+2LFsxKh7sw25anxP72xbZNjO3jNs8BPuBaa4MQIoCpbbFaStmjnrJFZXkZtdUDMgq7FEqpvPBcYmkK8I0JQ0qzgMOEf//GeJbecDoylttJvMXP1PRR5lMo8ynccf54fnfNqay8bSZ1k6tp3hvK2JBh93V5em0z5UEt72GhxaKKzK/xitNGOL6WBXO64rX/9eZ2OqPJ3XzWheqSKcPj3nSZT2HxhebElkTxov5+Fb8qmHnc4HhIpDNq0BEx9a+dwi5Bn8KVp43KKBfqVo0yFKZ90QYsiP0VjRDix7H/tGLVlwshpgP7pZT3x875VyHEH4G7hRBVwAfAlZj6zVe5sQ63GD4omL1vNg9+eM5Yzh5/NC81fcJ9r77v3oFj5JJwMoCZ447mlc2tabfJJw2v4Kd1E3j4je38yaHY/nDnmmmjmFNbxYpNrfg0kVSlU6bZ6zAowO+unULQpzF8UDAtkdTUsj+jt6gqClGbeOue9jC5VFy5OUlHVUTGY/3u7Z1pTocizDK7CdUD2LS7jfU729IMd3+/yo/mjOfs8aaa21mxKS211RVJ75c1taV5b4igX+GC+/NTFzYkfHPSUJ5466O0x/oH1Lguh1vJvryNcgn4acr/XxP7dwdwf8L2K2L7XgEMxGzLniOldEW/WQjxus3m1HBJTpx/4lCec6mRZOunX/Lg6x/QWcAQwHzLxlRh3i7an0ly48wxaQIv65rb2P7Zl1w+dRTLN35S1GTgUnLBxCr+tN6pv8iZKaMGsG7n/oLn3wU0wYhB/Zhx92u2IjodUbOpIvWDMoBPv+hk1FHpP9HW9g4WL0/V3Up5vkPi+YzjBvOV8gD/9scNGQ2lU5XZGWOP4s1t+VVWhHWZVZ/YutYrgKqa39zfrNzOf725HV2a72OqJxvRDSaOqIiLFFldkgvnTkirGa4sL6OyvIx1O/aZx8/xiqMA1804lsFHBFg4dwK3PbchvlZNge9MO4bLpo50tfqix42ylLkN3JJSdmAObf230q6ocOobWuIKVm7xQmP+hsQiX0fHkM7PqQj62dNuL7B+49PrmTtxaE4eWHehCrMeOCIl/zrrOMYNPbIgo7z2o+IG6cz5alVcrN0JpxzdvCWNlGlqmqGxtIadFMwAEILTjh3I2zF9YoBLpgxnfFUF46sqqKkq55v3rXI8t9Mned3XRjOr5mgWLtsCQhLOodLCp4BPy22en0Fyt1xXmCE53BPWJYoimPvA6phwvoy/H/PrGzky6EvzmMG8w3CK+zut59HVH/Lwm9uZN7vGvAOJJUujBjy08gMumzoy5+PlQl7NI4cbQoj9FRUVFdmaR8D0Xmbc/Vqv9RSLQVPgx+eN5yfPb8po6Gd/9WiW//3TbltXJhLvEjQFzv3qEJ5v/CTTU0qCFUvN1InmVyCbnUhsUGht7+D0RX9xNKiJz3nkypPZ+XmIiSNMY5xIfUNLkueXDU2B/77udE4aOZDW9g4aW9q44fcNWSs/fAogMstl5krQZ3b+ZdPS6OdX0Q2D62aMZvaEIYTCBkG/Qt2Dawr+jfpVgU9T0sJGN501hlu/fnwhh3SlecTDAbtJCYcKl5w6nDtfzGyQAV7Z6KwZIHD/y5YpX5W41qhBjxhkMF9zttZgRVW49etj8SmCMp+CTyUtGZea3Rc53I0JATs/D3H2+Mo0gwxmrPV/bpqecwpEU5V4MquyvIyzxx3NojozidYvQ9mZgYh/9lbi0krMWUnxMp85KDXbZKmIbuQ0vur/b+/cw6Oo0vz/PVXV3WkwhEsGCSRGEDCESURUlOuAl1mIuiNhdnZGYQYdV10vM87jbxAYZ3/uOgZ0H9cRL+sVcdFx1pEwo3L5KYKGAIqAEMYQucklXIwIxAhJurvq/P44fTp1OVVdnXRIq+fzPHkgnepTp6q73nrrPe/7fU9HmJzmgtW7UfZYNW54/gNhLDmoElZh6QMFQIvA23+2am/a0uGADAhffFso6B3GqVaxkMo3nZc3HPRVCJGsTDadzxA/vSQfS7YcQmfU9inEPabqhkqYEbQ7jS0+vMM7Jg7GmPNyccWwsy0enRlz141VtQ0Iagqiunc4oDli4IFln+D+N/6OuWXDMHPsQMc2O440JT2D4YACCggXs/gi2tPv7cFCl96DukHBZ2oYNFGOffdVQxOLb80RdnzLa47g/jdrHWMEVQVR3WDntx251gnv2PZ5EALMm1aK2kONWPyBt36N22cZVNvaVKUDaZQznIBKUNQvG9sPfdVlc/BzCSRbVEyn6QxpBN2DWqeVKqvtMMo6BV7+5Sh8uPc4nnl/j2/hnKBK8Pjq3Xi2aq8ldlwxtcSxeOW364YZHj/lhs5smCu31OO3ryVP8Y8ZFL8rG+YpuPOnjQccr6mEdb02hwtCmpqYE198M8PnV7G8DhRMSzmoEhiGkZbvUEgloCAIaQqaozEYFLj/b58gahiYWzYMB46fdr25uJHOdDhAhi/SxsHjzenMhLOw8/MmXDqwV0Z+WATsonG7YFgniPSemTsnDcFLG/aldUwz7Qk5KgDCAQ33/PB8PDX9InQTlOqKGotGdIqIbqCpNWapDCsfmY+qWZOw6KZRqJo1CeOG5GJOpXvXDT9nmPe/A9gayJzKGl833Kipk7QIUeguS1Nw3zXFDqfWjwErK83DvGkl4EcV0Snc1hOvKc2z5CdfP6rAs2FqzKD4651j8OhPR1iU6VqiBuavrMNtE8/D/dcWI6gqCAeVpOdVVZB27QvpKaeJcFBBZ4jEcW/Q3vMtU/jRiP54Z8fniAgepa8a1he/+eFQXLugOm376x5UEdQU34tTZwoDwOubD6Cgdxil+TmOLuZBFXj8+pGoP34af1hW6zp/c8dmsyf5X29/6rmglkr/u77ZWaz4hChACt9at0d0UZPTmGFg/oo68Ns1j5Eny+flGUwK/C0M3j7pPPzbtcU4eLw5kc9991VDsXjDfuETi05Zx+yS/J4IqtbydH7uubf+4LIdSc+rbgAVK+rQIxxIm3RnJjpfXQIh5D37D1LIU26OGGe8UioTaGhqdS1iWBvvBpzOIINOKYIZKjLy6sZ6jH9oDap3HUPF1BJLQQQFwVfNUZSV5kH1WBA2x4437z+BhqYWNDS14JmqPR2en27QhJeaqohWS9RwLTU2V81xfQpCCFpjBlrjRtEwKCpvH+NpuBqaWhI9+E770J9QCFB3pAl9s7NwUWGvhLHvm52Fe354PmZcJk5Vq1heh+ZoDC0x6z7M537eijrfWteRmIF7l9R0mfaFxIV0xpQ6g2Slru1lvYdEowKCFduPphyfNWOfcfmFA1DQu1v7B+xkWmMsBFGUl20xvlGdNdasqW909OkzU37hAFTvOoYJD6/BzIUbMeHhNXh5w34E1fZJw2oKsWhmmA2XXUTL69sRUolnfq853PLkDRchS7PO1xxLdsMtg0kU9gFY3N9NCKihqQUvf7jfZU8Uv3hho8UTD6gkcX5e3rA/aZqfnahOffUQ9IMMX8ShlE60vxYXwffVqpp7C6nkfXYWmsIq8/h3TlNYd99kq+Tp5nRUxxMdLA232/PXNh3EmMF9OjRmZxNQFGw90Oh4PFYJwZ6GrxHxKA9csqUeS7YcYpoM8Ygv85JTv6EGVIK/3Tk2kdlgDhs0NLWgsE93vHHXOBw+2YLG5iiONbXisXd3CYs8iEKSOh483FJ7pNFh1PzEkgt6hx3nJqSRRNinYvkO4QKqKKxy8HgzgqqK1pgzI0o0BgEwbkiu51NJSFMwbeQA/GXTwXatO/hFGuU0Mm5IrqXip6v49RVDkZ2l4cHlO6ApBLpBse/YKYwbkosxg3PbbZQVADdcdg7+9OEB37oI6c6PiBmsgrArUOJl6MmOKWqw8l+7gTkV0fHoqp2I6QYCKoFKCFpsxkslisP+aoqCKSV5WPpxve9y76CmYH55iTA/mcdteWZH+YUDsCTexVrUXURTiG9tBz52qrFkAKjedQwR2/mYNjIfV8Y1LQaffRamP7/R8ndzWKWhqSURWy7oHXaEJziiKmuFAKtqG9C3R0hozCedn4vfTi5CcV4Opo8uFDaZrT9+2vP4/CIr+jxIpaIPYK3lZy7ciKYuzld2azcU0hRMv/QcvOBhlL0aU2ZpBAYlnp5eqhCwXNEz3KquXfzskgKMHdLH9aZgFqcB4PnURMBW7u1/Z5kDxOFpdg+piET1pB5aUCWYcVkhRg/OFZYZt6fyNKgSVM++PKlRrT3SiOueWG/5fpjbRCV7748eX+eI44Y0BWvvnZTQrbj+2Q2JODXAwip/umU09n95ynKjmT25yHNBVUQoLkJiUOe1Y/5sy0fmY9G6zxz51O1oCyUr+jobt7uzW0yso0wamit83a3dUGvMwOIP9rl+6Aq8W8G3pzNwMroFVU8Jykzi1Y8OenrpYwb1QeXtY9AjHMCcyu2eBoFCZJAVzCsvTVTJmTtgnGplBtnrTGUFFMy4rBCvbDyA3/x5KyY8vMbREaM9laea2lZNaF6ANFO5pR7XPbFOEH5oiyW7vXfRus/wjwKDDDCFOb7vgt5hENt3hSgE4aCSWCDk6W0Vy3cg5FFlKKI1ZqA1RmEYNNHVhWMXsy/J7+mQ8kyXprIMX6QZ+5MHAXDxub1SVtZKRjig4M4rhuIH5/fFH5btSKoFwIno7kpfBtgXwsu5SKeXzMbToRKCaFpHdaIqaLfSm1/e2dGA93cdg0acnm4yugVUPDD1+yjs0x0FvcOomjUJq2ob8ODyWkt2S1AlMADoOnV8TpRSvPzhAUs8eu7S7Rg3JDehmXHidCTlzzCms3jwonWfoWJ5HVSV7Z9XCfKsCdETFo8l20MmFVNLfK1xmDNG+LqNfZzmiOEQaFIV4qmsGIw/kdjDJQALbSgGRdkFeVj+d6ugPje8Bb3D0G3XXLqKSKRRTiMHjzcjHNAs4QsKpN0gA6ypZUHvMC4q7IXBZ5+FmQs/8m2YP9h7HFma4ohnArAk6fuR/mxPSbIZwwCiZyCEphssNur3HLVXTzgSMyDW0kvyPl3H3MrtrJw4bmxGnJODqO3ppFUwqXBAgW5QTL+0EK9tqrfcEBQQ1NQ34qvmNnlLkSHywqDAvGU7sJTrZcdtFH98L8nvCSo4r5rCmorW1DdiduV2REw3i1mvb4NC3ENlAFuonFdujUWbtZH54mVDU4sjT7o5aiSyjYSfJaWe3+2oQfHGtsOwP5tww+t2g0hHEYmMKXuQakz5TCrFBVSCh6aVAgDmVG6HQf13SPajSuYXTRAXzVQUJC8ZD6oKFIUtTB06eRp/XLUr7R62SoDLi/qiatcxBFUFEV13xDFVwsIGoBStOkVWQIFhMLlK8/crqBIYFAgFmKi9KB7K0t7835BSIagqWHjjxY4FOIDdsAOqAsMwUspWUAD8/tphuLq0v28jZ5bNtVc8BlWCm8YNxMLqfQhoBK1RHYQwkSTRTY6THdIwc+y5eG7tXqFWM1fKA5zC+j4RRqOkUfYgVaMMsC/H7CU1vrUPOF4LbO7vAQyanostnZ0m+HiqyspUU32U7whBVQEhVNjaJ6B4l1AHVeCp6RejND/HojPRwi9iwtKpAmrHZSi7h1TEdCYtWVrQE7/581bPBeKgSrDwxktw80ubPG/6moKkHqiILI1l6dAUjXf3EOv+cf8bf095n26ENAWEQChWb8eccQFAGPLJDmlYdNMoFPQO+5YbBdoW7gBYPHOAXedzKrdDIW3aICLBpyTIhT4vOlrRxykfmY+/3jk2pZ56CuAoyzUT0ggTd7EVf0R074U5EbyBp5nuQRX/dHH7SkTdjvP31xSj+t5JqCgvade47YGAadteWzrA0TgUSK5pcXVJf/TqFsSxr1stOhM6ZRVpFEBWgMAwqOfiZPeQipBGoHlsc6qVSUs+u3Yv+vfMEpapmwlpKsIBzVI5F1SJo4qU62+kurgc1Sl+d3UxHv5xqad2hON9MQPn9AkjjX0dEn3zvDpENzS14JG3P8X4h9qKbKp3HcOVxX0d4ZmWmJ4wqL26BT2LdzhBra3Yxlwx2NDUglU7Pses12vQGjPQHDUQ1Snuf7MWj7+7Ky3HL41yJ1B3pMnXIz2PeRkQhwAUwuJyRrwQJB2qaEFBO/qoruPK4n4OY+2HmAFcOrAXAipBVoBAU4B7rhqKstK8+Eo0FRpIv6RShEgBPPLOTry+pb5dce6lWw9j+vMf4JrHqh3nyAA7/y1RCp3Cschjnu8vxw7EnCnD4CdLO6Aw2cpbJ5znuR2PZfLKuUd/OgLzppXCfi+PGgbGD811nZ8bOgXmrajDuCG5WHvv5bjr8sEIqgTdAuwGw/WPzXFagOk63/zSJlx4Ti/LeAr8NeXlqMSZpSTKZuDGeNz81Xh89W60xqxiTse+bnUstpt/Dwed3387QZUA1MC+Y6csN4XKLfWY8PAa3PnKFqEz9Mg7O7Fo3We+j9kNudAXp6MVfRyuvpWMKcPPxrt1X0B04YY1BTdPGIRnq/am/dFfVBwQ0YFl2w4LG2364cPPTsT7+xGEAioee3cnFqzelTRmlwwF6Q2p+EF0flJBp20VeH4+Om5sp48uxDOCzzscVBCLGbhh1DmJ18w96XSDFaKY20YV5+VgbtkwRx5tUAV0StxvKPH0s4sKe2FgbncmpE8AUIKLz+2d0D/+29Z6/M8GJtXJ52tuOwUw8aM7Jw1x9HR0gwjCZ/ZsBhYyqHHtOs0rKe2L7eGAhoPHmxPnzV3TkMHDMAtW78YzVXsxr5xlivC0Oy8qltehrDSvQwt+0lNOM23qW968W8c6G4ugBLigoCe0TtCqcGPp1sMOrysVdMpCKadadcQM5lV2xCAHVeKrw0QyzuApNO1T8fV0YNZbYFoUJRYZyqkj+iMSZYtkL6zbh7HzV2PRus8sObkxgz1RPfrTEaiaNSkRg505dmBCgpJrX8yfdgH+88elrhooPP2Mp7i1xox4Bw/mhQKsxx03yF4EVRV9zgo5cnndUAhbEAyoJHH85myGtjm5f6d4JaU9EyNqGJZc5lS+l/zYa+obfeV3cyW+jiA95TQTDipo8aFwparE9a5bnNcD63cf82wh3xlkShJFUCX4wdDv4Z0d7u2l/HL3lUPw2Lu7OyXzwI3WmO7Lw1cIK83nmNO9wkGFNQU1jRPVKR5cvsOhkhdUVfTqFhQKxvMwknmRatyQXCzesB9PrdmdGN+cfrZ5/wlH3m9AUVBT34iK5Tt8nQNuIP2GUXjudEhjNxh7NkNNfSMUj8A1jwEX5+X4zmX2CzfGdmMvIqbTDucqS6OcRvjjlS/h8Kjh6pluOXASWw50rIuyF5mexhbRaVoM8mWDeuOuK4biSGML/rTxYBpm5g0PNVDERTKSEFRVi5iOOZPATe9YU4gjd9mraEHU3aNvdhZmjC7EBQU90dgcRY6t87OoMpX/Hoin8Lmeg4ACg1LcMn4Qcs8KoWJqCeZU1kAlzFlJnpLovMEkC1sAwLzytkwNv7nMdlTCQjj2LJKoYaA0P8fNukF+AAAgAElEQVSX4NgtEwZ1OFdZGuU04efxyoxBWRytE1rMeTK8fzZevHEU/vjOzjNiqLqCycP74ZoL+iEvpxtqjzSi8uNDnbq/gEpwTUkelm0/CiIwmm6YjSmvlgvEU9NmTy4S6h0bFJhbVoT5K+vaXbRgrq6L6DpLy8u3Lp3EbIYpplP075mVtA/jFcPOxju1n+PFdfvw7Nq9KL9wAAAWm1ZVgkCSakf7DcbvdZUTDlh+t9+M+mZnYfbkImH/P6At778oLxsLq/fhjW2HEVKt5zeZ4FhII5gxWqzhnArSKKcJring9/HIADrdIItau+/6/Gsc+7oVFeWlmD66EGt3HsPDK+t8L6gFVYKrS/LaqrvaOzeV5f1GdZp2MaKYYeCe12ri0pm6r87P7SGoEEwfXYjhA3Iwt3K7Z/kyl32s/PiQw5iaxW143cP8lXWYM2WYRVSHhxjKR+ZbwhIAE8Oyy3OKMAvJ8++qeUGrfGQ+auobHV9NCuDwyRZUTC3BvUu2QRShIwDerv0ckZiRMLzOG7+zsEVVFEslo/kY/FxXAZUkbirmpw37ueB6FeawYFhTcPdVQzF1JNOxLn9qPQKKAgKKG8eei+mjCxPjMDlQZ/ZGt6AKg9K0VfRJo5wmRC1xzjQaYWXSPPn+lvGD8MK6zyxdjyM6xXVPrMf8aewCzD0rhIdX1vkaX1WAuy4fgmH9e6A5qmPlJ5+3e64EwNI7xqLuSBPmVNZA9+jDliqr4qGPtounc+5+EYNi4bp9noU/QZXgth+cl7i4eQaD+bFaFKdVFYKS/J5YP+cKYdUY9wTNlWxR3XDtWs1xM3J8QWvckFxsiHeMcYO4FOBrqlhLwoy5calZA8PNkCa7rriHy88FD5Xo1MC88lJL8YlIr4ISYOrIAQDguFk9U7UH00cXJgx9OKg45hJUgV9fMQTjh+YmVcLzizTKacJcCx+NGWc8lSugEPz3jIvQv2eWpV27aLEword1x1i78xgUn9V8ugHfKU7J0BQFh0+yBqFFedn4xzT28TvTuEqdBhQ8cf3IhB4w4HysPni8WRinjcZowkhdOUzsfZm9Xs79b9bi4PHTuHXieb776XH4Yp6oY4cCoH/PLNz5py2ux+snj54oBEtvH+MQ3nfzMPl1JYopm2VBG5paTPFedi5nvb4tIcZkHsscurll/CAA4ptVa4xibuV2VO8+ZtGe5k87zdEYKAieWL0bj67a6asC0Q8yJS5OOir6ykfmo/L2MZ492DqLn48uRGl+DorzcnBRIUvkn+/hAVOD4toF1fFeZOmbh1cVm5nTUR13vLIFlVvq0RwxUj5nCpCxvfo4lMIRq7Uj6rYBsLgxAKHUJYctBjrP9wvr9mH8Q07ZTqDNMImq9rix1gSfhQFgySbneJoilqbtFlQT3aXNaX48Q4IvZvrpa1c+Mj9R0BLS2saaP6004Z3W1Dc6FuBiBhwtmnjxzcyx5wIgeHHdPkx4eA22158ULmCu2tFgkQSt/PgQKm8Xd8P2qkBMBekpp5nDJ1vYY1wnZrMRMIOkKW26DK9tqscrGw8kHgdX1TZ4lgJ3JIfYTpbG9BJ+fcVQDOgVTrpSnphD/JG58vYxMFJo4gkA3eNiMc9W7QWltN26C6kox6XKDZee4yvGaK9A42nEEx5e4xmWKOgdRtQljm0OR9jnwLMTHn1nJ/6yqR5BTUnEREvzc1zHfGn9PkeIieVJ24XpCRb87MJEuMUeshHJeCbzMBPNUEcXuoY6UuG5tXvj+snsWOevrMP0Sws9G0AATHXv8MmWRLm2qBt2R+PKme1qnEEopRPtPwC+TmWMyi31uOOVzcIeZ+0lqAK/HHuupeT1kZ9cgKV3jMGvrhgCVVUQM9ru1rNe34bxD63Bg8tqPfOcfZT/+4ILGYUCKp58j/Xje/KGi9At6K9ogJcZ/+vEwSntN2oY6NM9CIDCcDHIfnz2zsxffuXDA0Jv1QyXezWTFVBRsbwuob0RcdFW6Judhbllw1zH9hJd/+M7O/HqxoOIGRSnIzqmDO+H8pH5nmNqmuIonVaJteScLUaW4sphZ1vCBmbtCLsg/ewlNaje/YXrU4FZHN/euZpTmp/jmJt5AdCMSOg/oCgYPTg3qdQAf8LbXn9SWKSSDj1laZTTRKopcX6hIBg+IAfVsy/H4psvxdp7LwcAlD+1Ho+t2uVYWIkZzEuyyxfaSZccJYlnd5g7M/TvmeUpsGQmoreVGaeivTFleD/MX1mH1hiFSFstVe2FjqAqbH/23XFv1euRVhS+iOoGRM2rRdoKvHIvIHgqcjMStUcaHVkRS7ceRu2RxsSY91w11PE+wzAcBRw6ta5HEFBLQYwdkUGM6BTTn9+IG57/wNEthetNcNEht5tc3+wsPDStFCFNQbegipCmJBYA7R1PRHF1notsr6q8flSBI9TTGjMwb8UO3DDqHEs4JV3ZF9Iopwm3+F5H4a3pASRixdzTOO2jctCNkEo6XIKcpSmOuC73fGdPLmIymknGuOHScxLe3LzyEtcSYDvLth+B7nFnMeC+AJdudIPtT7Q7AiYn6WaYq3cdg24yEAGVYG7ZMOguc69YXucYa+bYgVg3py3m2j2oIqgygXmRkdh6oNHxmv31u64YgvuvLbZ8RygIpo0c4KlSF9GBxRvYQqGo/ZPXQmNLlKnDzamswaodn6N69xeYvcTqVXvd5MpH5mPpHWNw39XFWHrHGLbGIzDqPK7OjyOkkcSCH485L7ppFKpmTUJFeanwya81RvHqRweBeOqcucS9o0ijnCa84nsdRSUkcWG7lZuGAwpSUWskCsHdVw1BQGEXlsgWJjOPUd0QtsRZv+cLPLhsByg1PJPRVMIe8fkFAwB/vXOs0OtzvFchvhcoswJd1wOwOWrg39/8ROjl8acr88OOQoCy0jzXEIKbtgKPuc6ZUoSoThHQFMxfWYfKLfUO4zjiHPHio/31stI8JrQfJ6rTxELXoptG4a93jnXEwwHg2aq9WLTuM6GHy4s4vBaEW2MU//ryZkx/fqPjKcIrJFO5pR7lT63HvOU7UP7UeodOCDfqtUcaUdinOypvH+NY8ONG2xwiKc3PET75MV0QJr+aTqRRThPJ4nsd4VREx4PLazFm3ru4bfEmYcw6ZlD8w/B+vsYLagpmTy7Ck2v2IGpQtETFKXxJW0EpBHOmFFke96YM74dH3t7Fmre6GM1wQEFIIyDx6i7zBZN7VggTz/9e0mPQDepLszocVPBv13wf919bbJGePJO46QOLHuV56bVbCMHcs85OQ1ML5q+sQ0Q3EuEkvsZgNo7FeTm4flSB5b3Xjyqw5Nk2NLVgVW2DQxSLPwldVNgLxXk5QslRTSGJmLjdw63cUs+qEZN8Dm7pda0xHeGg84MXxaorlu9wLHZTg+Xpz1y4EVOfXIen39/j+A6KPPF/GT+IhUYEzVjT1TCVI7Mv0kiPcCDtHTw4ycSJojrF27VHfY01/dJzUJLfs90CLZygpqAkvyeqZk1CTX0jGpujuPf1bUnfFzUoSgfkOPQ9eJ5s1a4vHO/h7ZF45dfsyUWYtyJ5w1hKgSuL+6JvdhbKSvPwf17biqpO6JnoB/vqvOhR3mx07rpiCLKzNEv5tVvckhtRuxGKGazCkWcJzKmswbghuYmKzq0HGjHinByLQfZqrcTj1LygYnJJP4fkaFQ3ENAUSwYS/2z9yF96oSgE5U+tdxSdiPKMVUHJO8s6oq7Vl/bPyJwpQqmBK4r74f/VNljWciK6jhOnI4mFyI4ijXKa4HfqM100YiakafjeWSoONbZ6brf4g/348cX5Ha5A5AUOXKdW8akhHNOpUHCpJRrD+59+Ec9Ztt2ECPPKS/J7JtKheoQDuOe1bUKP3txSiF8o/7vxoNAgJ2sTZUYlzDC0p+GAfeHNXMxADSZ1Sghw3RPrEilwbkpvZryMqJ3WGMXiDftxzw/PR3FejqMKTVSQArBuKvymUL3rGOZUboeqsBuFvXx89uQiR4581DDQ2BzxVHrzA5/XrNe3WcqzZ08uEjZPnTqiP1Z8chQBhZXc23sc2jF/RqKS9Ddrjia6AGVpKpqjMRgU+M2ft/pO70uGNMpp4uDxZs+84DNBcyTm2eeNE9UpVmw/ioqpJfjtX7a1+0ZydQkLl3TU+0nMy2A3DBG6QTF/ZR2qZk2ySFBqLgZVpxSLZl6CcEBLPI4+9q64GnF4/xxsP9SYOA+ssen3UJTXA0FNwZNr9lhyanuEA/jVqx+nlPoY0sSr87yi8bon1gNo6y3ItTBmjh0oVHrjeBnRWDzmb79R/vd7ezDDpOlgRuRxdg+yPnxXFvcFAIydv9pyU/rL5nr87c6xliq9HuGAJRe5/MIBmFuZPDuJl2ETeDeMtT8BzF9ZhzsmDnZUnK745Cgq4xWE4aCC8qfWW/7upb3hVpKuUyZpcP+PhsePyUBUZ9edW254KkijnCa215/sVP3jcEBBVDdACEFAUYSZF6kY16ff34MXb7yEeTvttMrL/34ERXk9ztjNyP5oefB4MxQi7lFtGBQ3vvhRoiPHLeMHIaCKFb622qq+dApU7/kS7+88hrllw1A1a5JDr0K08BMOsnirnV+OPde17BkA0/oVFByZu1i4Ce0kM6LPvLfHURARMyhq6huF5duikIpOaSIEtGrH546nhKhOcfhki6Wc3K4NXf7UemGZNIU1fmwuww4HWSl+Y3MUcyu3C88dJ6AoCVF983VojoEDEGotu2lveGWKqETBsabWTikgkUY5DfAFls4kZlDcd3UxykrzXDvypmJaozrFzIUfCd/lV1G0NUbxx1U7HS2Ugiowb9oF2LTvOP73o4NpU4GzP/4X9A6DusTEDQoYOk14MLxFk1+4cb3/zVo0tcQwZnBb7q1dQ4E/Ppfk98T2+pOYv7IOBCy2+usrhuKuK4Z47sstc4dnWpjbP9kfkZMZ0dGDc5NWqZkRHVt782+5hy8Sze8WULHg+gvxVXPUsS8eUjHHc81tryK6DsPWs9JNVN/+nRFpLfO5up0LUYXq6aiOR97+1HGdpKOARBrlNOAlLyj249oU3fwS1dnje1lpHq4cdjamjRwg1ENOZaHRbZEs4KF6ZkfU046CoKk5iqUfH+qQQQ6oBMSmKGbXyP3xRQW+dKFVouBnowrw8of7Uy7weeSdnXi6ak8iplo+Mt/14k54ZMvrElWOA3qFPeOMPHPHrvWrG9TSxoh/v8yPyDzFzG0xsDQ/x9HUQFOcmhxmT9zt2No7nujGYYAmyrBF+xLFc81dSfiNipehz55c5Np1RCTy7/cmw8/F4g378cz7eyzXRURnGUAhzV16tD1IoxwnLkBkx5cgkddjTlBToFNquaurBAAh0Igz3ucFX8EGGrFki1i4PaApeOkXF+PAl83Y+Nkx/HXrEf87iJPMIKuKd0VgVKeoWL4DgQ4KBqkKScQDRYtcDU0tvgXsT0d1LP5gP2ZcVojTkRhe/chZGeYWfgDasl9mL6lBUV42ivNyhBe3OS2NhyP8xBm5roXduIraGJkfkc0pZtEY08iw3wD+deJgPPP+HmhKm6SlvbOHyBN38x4f/vEFDonMZON5GUvReRQ5OuauJOUj8/FVcxQVy+ugqiQhfzpz7EBPKdD2wHPALyjo6VhLCAc0PPrTEejVLZi2/UmjnAb6Zmeh/EKx59oisLosKyd1F7I5GsMdr2yB6tK9IaiyxaRxg78HDAb69gi1yygnw0/ShuqhrRvSCGJxcXt+FgIqQVG/bGw/9FViu/ILByD3rJAwB7ShqQVLtxxKISDBbjYvrNuHgMq8b/M57B5U8buri/Hl162e8qR2PWqz1m5zxMCJ0xFPI+qFKNNC1MbInJZmX+TjT1N24R9CgJvHD7SItgNijzTZTcTLk3Ybr2rWJEds3gu3UmhzZgS/+fFEHfviaLoRFZHw8ux07k8a5ThxASILhJCTAJIqVzc0tWBJEuEZP3QPqvjHC/KEXlx8Rq6tdJi27BjknhXCqh1MfL5/z6xO6cfnK94cNSwGUyHAjWPORc/uASxYtdsZYqEUdUebLC+9tukgXt9cD01REDN03PaDwZg+uhDVu44l7ZXmRTSeq2rGHIfNztLw4LIdiLrEXrge9VfNUcxfWZdIZ8sKKKDU2UkllTijqI2Rm5fp1uCU38TsxvHZtXsx3dauSOSR+rmJuIUAvMYTCQl5je/lXbvJGpgXR9NNOuPtXkijnAYOHm+G4nt5zJ2oriM3OwsThvRx5NNmaQoUhSBmW6IPqgoUha0q1x1pwqzXqy3tg35ycQGWbKmHYcDVyHQGBMTSz82gLN3NoOJYdlTQHysWf43r3C5YvRtPv78HBk0t7JMMlRDLxTVz7EDk9+7mmfZGwAyAuQiBe6zpjjO6eaZe3qRfY5vMI02VdI7n5ZEnWxztDKOcbE7pQhrlNFDQO4xYCnrABBB2+4jowOOrdwvfQ0EdzSwBwKAG5k1lDR+nPrnOYqyiOsVfNtfjxRtZvm5zNIZfvLAxpdS5kef0bFdnbVGDzXQIBHWGyJBOKfr2CFleK83PQSTmnuLYHDVcVe06I84o8kyTeW5+jGO6vb/OGM8ttu22OJoO+cz2zCldEJGgiIRBCDmZk5OTc/JkcqNkbn6ZjFQyL8JBBbpOMbdsGI6fimCBwGh3C7JCAUBstAIKwUM/Zv3KHnjzE2GKlNnPVwlQPnIAbhw3EMV5OYlOy6oK14Wwbwpui5SaQvDwj9t6ujU0tWB0xbvtKqzJCiiWIpfOxi2HORUxea+Go+mcU7qxdwFPV0umM4RwSUQaZQ9SMcoA+4I8uHwHqJGeJqAqYXrFoYCSaDv/wFu17TIUmgKsn3MFltccSXrzCGkEa++93LGA8/DKOry+2V+2g1+uKe2Hd2obPNvOpxOvkmqzMd28/wR+/sKHnmXLouozABllGM6UcexKvsHHKDTKMnyRRsyr5+/UHsVza/d2SEyeZ2nE4ulYD7xV65m0EVQJogYVbhMzgLW7jmHeiuRFLlypzP4Ff2NragZZJcz7dgtlhzSCf7t2OC4vOhtzl26HYbC2Tly3Ysrwfli69XBK++QoAECc+/aqBjfHXAt6h4UhGDPm6rPmaAwHvmx2iPt0NZ39qJ0JfNuOUUp3dgLvfdqAF9ftS4i8qwpBOK5Z7LcimQAOAXGdigtROBGdemoqr9nxuS+PVBR7PHi8GZrir8UTgPg8iMUoEsDSqYHnt5aPzMfsyUUACMJBlsEwe3IRpo8+F91D/vfJCaoE904p8qXLbKYlpieOm6c5uh6fShLVZ/u/PIWbX9qU0PFN1gJKIvFCespppHJLvbAkUzcoSgt7YtP+E75yfEMaM0yxdrjZbk/bmgK8Xfu5r32LFmaY5+hvPqoC3DZxsGPRkgKoKC/BuX26J/J6uViQPed0/so6VN4+Rri4mYyITnHiVCTl5rDmUF6ywhSFMEGk9uT5SiReSE85TSTr0ffhZyd8hzJaYwYiujgMYUcl3h9iSCMIqAS/GH0uQprT69QUkmiJ88ux5+LJG0YKe6z1zc7CvPJSi7C8m065lzzjsaZWbK8/ifKn1uPnL3yIcfNX45n39ggbWa7YftS1AMU8h9t+MNDx+sJ1nwmzI7oFVWQFFPzskgLH38MBLZHjKxKgN8NDPG5NONMpei75biE95TThpX/RXrICKlpjumtO7l2XD8YFBT1x1ytb0OyyUWuMIqASLP5gH6htXSGkESy9g0kuciGd1zbVu67U8xzNmriq2raDJ4UpfLpBsbD6M8frAPDHd3c6MjheWLfP0R4oout45v09STO/b580GBPP74vFHxywqIMxrV/rtiGNYMHPLkT98dOYt2KH4wbaHI0lwhfhoOIZ6okabSpmjsanaepqLPluIj3lNOGlfyGCwOppaoqz+7JOmTKcyCPVFGDG6EJW+pnEdEV1iojOxOXtMd3ivBwU9A5j/kpx+x47fbOzUJqfg17dgujdPSDcn0FhyVowdwN2S6kzDHbz6B5inuz0SwtBkgiihzSCGaMLWWjFsBtYA1eX9LO0qppXXorS/JxEF2w7fH+81xuJn9esgJL4fPhY5RcOQPlT6/GbP29NKJilu6ux5LuJ9JTThD1p/nQkZtF2sKMqBKoCzLy0EKMH56J/zyys2H4Uz1btRVBVENF13DJ+EMpK81BWmodn3tuDlzbsg6ayeLNZBGbayHxfSmkUwJTvn40Zowda0odSKbU1575GdB2EJJfxGHNeH2zcd8JTb9oAEFKYqM7lRX3xPxv2J61AnBb35GvqGzF1RH9Hefqy7UewMF44w49XVJrMydJUYcsiw6B4465xCR0Org/spmAmDbKkI0ijnEbswt5Tn1zv+ggci3eEeGXjART07oY7/7QlHpukuGxQb1TtOoYX1+3Ds2v3omJqCX5/7XDcOvE8ocShX6U0APjr1iMYUdArITEJ+C+NFS1qmWF51U6tjXW7v0zq9QJtMqArP0m+IAkwbYz//eiga952RKe46cVNmFtWlHjN64mGv2432iFNRXPE8NQHNiuYSSQdQYYv4hBC3rP/wKd0p5m+2aw9eXPEQFBNfnpVYu382xqjeGdHg7DDLh9b1HkiFR5ctsMSmuBeflZAQfegiqBKMHtykcXwb95/AjX1jZ770uILinYCquLQ2z3/7O5MtNwrhy9OSCOYcVmhY2EuZiTXjo7oBu5/sxY/X/ghJjy8BtW7jiWONRSPC2UFlETYoTQ/J+kNKt16ERKJGekpdxJ+Y8yizr92vFS7Uo1lA0yYiDfPBJjRLezTHXdMOg+Pv7sHAU3B/JV16BFmMWPeJDOme0evg6qK0YNz8fKHByxPCKciOj7ad8Ky7aefn0JQJbhh1DmO7e1QSvAP3z/btX+fH3joxC4jyVPzzE8fbtoNtUcaE92fz4RamOS7iTTKcToi3SlCFGMmhIAQiqje1m159uQizFuxw3Msu46sPYTxL+MH4dmqvQBl5b4hjSTtrvFs1V7MiMtgzqncDoW0hQ/4DWJOZQ10w6rqphA2d00Rt5/v3zMLt0xg8xFtYyaiU7yy8QDmTGHdj926Mc8tK0I4oCGkkqS5x4G4toWbiTfLSAIQlueKlMDmVtZY4vbXjypISR9YIvGLNMqdADecRXnZeOL6kQCY6tjymiOoWF6HcBDQdYo5U1gbeTehIc4dEwc7RMt5h2De2h2guPUH52FyST80Rwys3/NF3OslaG7VHUYqqLIuJvcuqXE0wuQQEEejUYO2FYDwNDrzfMqfWp+YT1lJHpb//ajnAl9AUVCS3xNVsyZh8Yb9eGrNbktX6d9fU4yZYweioakFxCatpylskZG/pCnA6POcsqdmeCrbI29/mlhUFaUAmkt3q3d/4VhI/dPGg5gez/zgOcnSMEvSgTTKaYYbTrPwOcA8YlHVWllpHqaPLsR/v7/H1TgueHcXsrO0RNoaX2DihsIuYl696xieXLMn0SLoN1cNxR9X7bTEX09HYmhsjrjuE4Cw8zMA5IQDuKiQLRZyrQ9RRsIb2w4nrcgzPwU8t3avZY7m/7tJQvK86Q27j2HxB/s9DXJQJRg/OBfXPbEuoabHwyZuVXiVW+ox6/VtwvEWVu/Dm9sOWdoiZYoQkeSbi1zoSyPm7AT+mN0SNdASNVCxfAdUW4GEOVb8u7JhruNGDYoH3qp1vN8O7+HH53AqoiOiUzyxZpcj+0GnwJ8/OJD0mOy7tDfJNC9s2hcBNUVh6RgelF84AH2zszw7SfBFyfKR+aiaNQmLbhqFytvHoLBPdwBsPq9sPOCptTxhSB8AwDs7GlzkTZ1VePzzdAt3L91Sj9YYxemojtYYxazXtwlzuyWSVJBGOY14ZUIEVAVRW5zX7CXOHDsQt/1gkGvPOZ0iacmxOaXLDAFxFKYAwEdJxOtjBjPKQZUgrCkIqMB9Vxf7XnCM6kbS7IrKjw+hoaklaScJTt/sLOz/8hTKn1qPmQs3YsLDa/Dyhv2eWSEKWJm7l9GO6DpOnI5YjKrX56nAmfkRM5CodpRI2os0ymnEKxMiZhiYcdk5loo6vmLf0NSCm1/6CE+/v9czu+HnlxVaKtSuH1Vg+d0tpaslZrRbrzikqfjh8H6IUYqgpjKhIIEKmjmtjs9nbtkwR6WdHfPTwlzB0wLvJMHT8mqPtD0J8JTBZ6r2OEqdzWQF1KRPGboB/ObPWzHh4TWJ4/P6PH2kXUsk7UKK3HuQqsg94N6BRFMIwgEVEV3HrRPOS3QVrtxSj9lLtnmmxAFs0euZn1+M/j2zLClcomyMVLqgtAcuBg84sxfs87FXABoUlji2vUsHbxSgKUz2c155CQAkxmiN6VAUYqm4yw5pmDn2XDy3dq8wg4OVebs3nXU7PvP8+RoBb0QrynBRFeDNu8ZllJ6yJKORIvdngpL8nugeUh0ZBzGDoqk1BqBtQY7HLN0MsrlFk06BO/+0BQDLo+UpXSKBb7c5pEJAIa5lzgFFweIN+/Hc2r2OVkP2+djTy6p3HfPM7+0RDkAhLOQCGImO0ZYqQt0ZBpoxuhAzRhfi4PFmR1ZIxVRm2EWyqm7HV1PfiF7drBk0WQEFN724CYBhGYd3MwmobLEzkzqPSL55SE/Zg/Z4yg1NLZjw8BqLJ2cnO6Rh0U2jAAAzF25MGGs7QZUI46DJesD5mYMZlVg7TwdV4J4fFmHB6l1Cwy7yPO1z8mrR4/Y30byDKkFAU1xvMAGV4KFpzqwH0T5W7fgcd/1pSyIf2w1NAVRFSeR987j4LeMH4cV1+yyfVzigIKpTSy73me7RJ/nGIvSUZUw5zdhjq1zP2Axf4POKWV41rK9Q/xhIrtfL5xD0aNoRDii4bcIgLP/1OPznP5VaYsHzp12AqSMHCNPZNIXg1gnnOUrI+ZwamlrwyNufYvxDaxILcfYYtKhcHBAvrAVUBREPI8rF5kXnwL6P0vwcYcjp1l0AABiQSURBVMze/OloClOLa405M2hEseuobiCoic+FRNIeZPiiE0jlkb38wgGWwoRLCnvh368bjtyzQhj/0Grh+BHdn84CIQoCiiHsS0cB3DR+IPpmZ6E4L8dRwVa5pV7Y+UQhFIW53dEq0BBev+cLLHhmdyLskSwH2O7Nim5SOqX4+ehCYQdugKXd2UvQ3Txxc66zSggiMR0UxOLlEgIENQVR3emZB1UVlw7qjVU7GiyvnxZUNkodDEl7kUY5zZgNAo/7isp2+bZ2hbfthxuRe1YIfbOzMP1SsTG6ZcIgz0fjti4oVgPHFxtFsVxzLJi/X5RBFtGB2bYqwIBKMGV4Pzzy9i7hfMyeo9eNqnxkvuMmVX7hANw68TxXfYxTER3rdx9LnGt71aM9vmv+LE6cjuA3f95qCUdoijN1se3YdazddczymnlK3UNqos29DF1I2os0ymnEyyCIFuS8dIyrdx3Dyx86BXhCGsGUkn7YvP+Eq+aCWyEGgbfmL7+hnDgdAfVIZbNXARJQvLX9iPv2ButscsPzHziyMPixz15Sg749Qo6bVOXHh3D3VUMxr5x5uIZuwK6T/8g7O5GdpaGsNM9Xvzz+WTQ0tTg88+aogakj+mPFJ0cd3bV5TFl0c+geVPG7smJcWdxXGmRJh5BGOU20p4GmmwRkOKgI+/2FNAXTRrbpS7i1bXJrZRQMKK6av+YbSmssljRFz4ymqKCgwpLtkKYkSszddJgBJk40c+FHwvzfg8ebEx7uqtoGPPDWJ47Fuorldcjv3c23WD/AjPPsyUWO9MEVnxxNKOaFgyShU1JWmodn1+4Vzl+nVBpkSVqQC31poj0NNEUFFxVTS4Qly90CKirKS1D58SHPtk28lZEiWNLihRhAm0ZyQ1MLGppaMKeyrSDDzSC7FWDo1IDIsf7l2HOx9t5JKMnv6UvzOWY4DXtLlN2kAHa+rizuKyxICcRbTqWqc8zTB80QCix4dzciuoHmCGtiO39lHQCWjmjvxxpQiQxZSNKG9JTTRHuFz0XxZtFjNVMypo5O0QoIauobceWwLIu3bqZbUIVB22Kd9jDL+MG5/gorKHuUt6eK8TxgvoAW1Q3MLRuGUYN6Y+mWQwhqxLPizouQSix9/Xjln927bY3qKM1PXedY2N9PcC74DXbckFyoimIRayKg6BEOJBoRSCQdQeYpe5BqnnKyRaZUEMl0LokL4NgJaQrmlZegsE93R95z96CK313dFuv0m8OswKlJnB3S8OhPR6BXt6BQHN68yPnHd3ZaFuwIWGeSLE0VVuW5kRVQUHn7mMS+AKYvcdvizZasCU0Bnp5xcUIsKRWd48ot9UkLS3ju8cHjzcLccvONTxaOSHwiK/o6G7vXC7CCBQApN9S09/srf2q9q9FojbEwRuXtYxwedswwLLFO0eKiEHM5YZyWWAz9e2a5lhHzBbTaI40O/WE+1KM/HYH+PbNQ/tR6x/tDmoKy7/fDsu1HEdAIdINaNJqbo6xRQEBRLAaZHSfwq1c/bpdhLB+Zjx7hAH716seO9Db7UwbgDJEAbWlxydYRJJJkyJhynHT36KvedQxj5r2Lm1/ahJtf2oSx81cLhXz8jCWKMYtojhiOmKdBgWpTGpff9lGiGHFUB370xLqkx7H1gFgpTSEEvboFUZyXg9mTixBUFXQPqQhpBL+6fDDmTCnCik+OJnSg75g42BJDjxksY+N0VBz0Ph3R0RI1MHtJDWqPOOdgjqPbKc3PgUHtC6sEC352IapmTbIY+X8ZPwghTUG3gLM6RxaOSDqKNMqdAFs4q7HksEZ1ijmVbYtyXgbCjh9DyhfEeMzTvF+zkRItLprV5oJJmpnaj0PEiHPEnrROKcJBBZVb6pk2Rdz4zpkyDNNHFyYyNLgO9OOrdyVVdxMR0Sl+9Pg6/NfbnybmWbmlHhMedlYZ8s8BgOO8zCsvxZXDzk54vXyMRev2AaD42agCRzPXiK7LwhFJh5Dhizjp7NHH8oQVJFqMxFEVkshBti+KzRw70HU8t9QtM8H4gtjB480IqtaUuIhOcd0T6zF/WrxAIx4a4dq//Xtm4fJhZyf+LwotiI7DrUIvHFQcehoA09iY+uQ6h1Lc/JUu6WyqgpZUcvNMRA2KBat345mqvYkegPZ0xa+ao5i3og6qwkIl88pLXPvuiVIeX/5wP+ZMGYY/LKtN3ID5k4mMK0vaizTKnUBB7zB06vRsdYMmcpDNi1zc2HoZ5uTKbyThoYm86ohuWOKd/MZgb1tVMbXEUoosamRq1jgWVei1xnQWE45ajbJbhgcPzThF8nVW92wz7t1DKqIxw9diYWvMwB+W1SJo0xFRCcEDy3ZYMi/uXVKDdbMvT1QHmhHF4ltjFAeOn7ZkY0R1KuPKkg4hwxedQN/sLEyzeUoqIZhXznKQk7U9EiFK3TIzt6wosdBWfuEA4TYqIVhV22ARireL7nCDUjVrEv7nl5fi/muLLTHqgMqOo3rXsUQ4YPxDq3HvkhpTnjP1rVAHAK0xHf17WsMqmsIUOu3HzCvn/nrnWN/jxwynPkVrTHeMHdWpa+eQgt5hRAR6GIs/2A/NJjgl48qSjiCNcpppaGrBqh2fY8kWa7mwpjI1M79tj+yYY8GhuBEIaQRBleD+a4sTXrZIT4NzKqLjD8tqcd0T61z3ww03AFxU2Aszxw7E+jlX4PlfXIznf3Ex1s2+HOOG5Fq6f7TGnEUfIZVAcwkHq4RlWvDjUBSSCJlUzZqER386AoqiQJTazCvnivNyhLHxkNtObSRpiJLAHHO+dcJ5jr9rCvFs8yWRpIoMX6QRnlusCLpcBFUVB48346LCXsLiB3O1nRv2NLnmiJH4lxcuJEt5S3iMLv3qTkV0PLi8Fv/x1icW4forh7U9im/efyJpWp1OKdPVFOzn9kmDMaWkH657Yj2ANq967tLtqJo1CTUHTwr7EQZVVrLNb16iwpu7rxqKxRv242mP7uCA2Cjbm8Lac8VnTy5CSLPG681aGX4LViQSL6RRThNu1XSc1pieKBfmXm3F8rpEPq7XhWyXojRnA9iLVcYNyfWV8hZUmVA9AYspm9sb8bi1W2xUlA3CheGDqrjtk5k+3YNojhgIaYql0o93/Hj6/d3C99007lxHRxF+0+D0zc7CPT88P2H0vSoJQyprORVQFejUwLzyUodSnnlhb/7KOtw5aTAeeWenZZwVnxzF87+4GAe+bMaIc3JkOyhJh5BGOU24eagBlSCq08QjOjckM8cORFlpXtLKM2547ZkabgJIVbMmoWJqCeZUboeqsPcoBILCE4K/3jkm4W2v2H4Uz1TttXioPJRhLj7hNwguMmS/IXCVuV+9+rFQkxhgxu35X1zseJrghl5TVEf8VlOAhdX7ENGTCz41NLWgOWLgrsudBtRyBhSCv5mqBc3jiD5PBQRBTXEsuFKD4qYXNyGkdbySUyKRRjlNiLzHoNpWyWZ+ROeGRCTnCVhTy9wyNbjIj0gRjUEByjxhgzrjrHPLihIe3dzKGkcFHuAMZXzVHLV497MnF6Ekv6fFoPXNzsJv/uysjDPDjRiJn52QRkApMHtyEUrzc4SZK4QQBDRiEUsSKcBZ1e50R7gBsFbpuXm1os/zdFTHI29/6gh9sMVSmvDKZfaFpCPIhb40IWoDdUF+L8cjfLKVeXORw3VPrINIm+SBt2oRDipJZT9PR3VEdIDGhYS6B1XHwqCoJNrMqVZWJffbv2zD/W/WIqIbidfmr6xzeJi1RxqxdOthz3PVqjMDxjM/WmMUmqpg/so6VO86hnnlpY6Mj9+VDXNmS5gW1PgCq1XtjjoMckABHrhuuKNKzw7/PO0LhxFbPz5VgaPYRmZfSDqC9JTTCF94WrxhP555fw8+iq/am/Fq5SQKSYjQKXD4ZItQEY2XZJvfHw60CQnZjahbSXRQUyyhDFF4WFRE4jaephKEPcSIzNoRVbMmYf2cKxLpaVw3pEc4kDjeiK7jlvGDAHgvsIZUAmqKnauqgt8t/buvEIOXJkbiHDjTqGX2haRDSKPcCTy3dq+wCzXg3cpJFMfUCOAmXuZX9jNqGOjfM8sigclxK4n2QzTmzBhxG2/RjZcgHNAS4kpucC/zosJelowPwHrTe7ZqL15ctw/PVO3xXFQEAeZMKULF8h0AxGEkL0SaGGZCmoaZY8/Fc2v3yuwLSVqQ4Ys0IxK754Q0ghmjC13fK4pjqiqBZtN/CKgkkbpl79os0rbgSmui7tLFeTm4flSBZfzrRxVgfrl3R26grWDFjNt44wZ/DxcV9rLkF3cXtNv242U+t3YvWmPuOdIAixtrCkt9+8+Vn6YcRuKYz6fbfGeMLkTVrElYdNOopGERiSQZUk/Zg1T1lAG46hVzzeNkF6wozQ1AIpuCazQkG8e8WFj+1HrLfLg2sD0WvPVAoyWlS1RG7VevQzSeaH7b608K09zc2Lz/hFDP2ExII6goL8XcSmfzWK9z4EV75yuReCCsdJJG2YP2GGXAalgjuo5bJ5yH6aMLUzYAbgLyqTwai4xYdkjDoptGCTUeUp1XOjDfQETpafZt7Tc9c440N5QiwX8gPWL0nXUeJN85pMj9mUIU600FUaqcW/pcMtrbpsrvvNKBWSApmQfKwwn2bYvysi2euSi2zvWRU204IJqDNMaSzkJ6yh6011PONNLZpqozEHm/ycILbgp15uPL9OOWfOeRnvJ3lY567p2NKOtEVBhihnurbpWN44bkZvxxSyQipFH+jpDJj9wdCbEkM+iZfNwSiQiZEifpckRpfH5zfdMZM5dIMgEZU/bg2xJT/qbQ3qwGGTuWfEORKXGpIo3yNweZpib5BiIX+iTfXmTsWPJtQcaUJRKJJIOQRlkikUgyCGmUJRKJJIOQRlkikUgyCLnQF4cQ8p7g5bPO9DwkEsl3G+kpSyQSSQYh85Q9IIQYAEhOjmwZL5FI0ktjY+MBSqmj64UMX3hjAFAaGxu/8rk9t97iRnVdT6bPD8j8Ocr5dYxMnx/QxXOUnnIaIYScBABKac+unouITJ8fkPlzlPPrGJk+P6Dr5yhjyhKJRJJBSKMskUgkGYQ0yhKJRJJBSKMskUgkGYQ0yhKJRJJBSKMskUgkGYQ0yhKJRJJByDxliUQiySCkpyyRSCQZhDTKEolEkkFIoyyRSCQZhDTKEolEkkFIo5wGCCEhQshDhJDDhJBmQsgHhJAr0jh+HiFkPiFkDSGkiRBCCSETXbb9R0LIFkJICyHkACHk/xJCHGqAhJCehJBnCSFfEEJOEUJWE0JGtGdMQsglhJAnCSG18bEOEEL+TAgZLBhrDCGkmhBymhBylBDyGCGkm2A73+c02ZiEkIsJIUsJIfvjYx0lhKwkhIzJhPm5vGdW/HPemglzJIRMjM9H9FPU1fMzbXcJIWQZIeQEIeRrQsg2QshM2zZn/BpJCUqp/OngD4BXAUQAPAzgFgDr47+PTtP4EwFQALsArIv/f6JguylgcqOrAPwLgAUAdACP27ZT4uN8BeDfANwB4BMAJwGcl+qYAF4HcCT+t5sB3AfgKIAmAMNM240A0AxgE4DbAPwBQAuAN9t7Tv2MCeCfAbwZn9cvAdwDYAuAGICrunp+grH7xT+brwFsbe94aT6H/Dv4KIDptp8eXT0/03c1AuBtAHcCuBXAIwB+39XXSErX+5k2YN+2HwCj4l/Wu02vZQHYDaAqTfvIBtAn/v/r4G6UPwGwGYBqeu0P8S/IENNrP4mPcZ3pte8BOAHgf1IdE8AYAEHb+4bEL5xFpteWA6gHcJbptZvjc7m8PefU75iCc9UN7MbxVqbND8AiAKsBvAenUe6SOaLNKF8nmnMGzC8HwOcAHksyvy65RlK63lN9g/xxfMgPg92dz7K9Pgfs7pmX5v0JjTKA4vjrt9he7x9/fbbptdcAHEI8T930+jNgnkEg1TFd5roZwIfx//cAEAVQYdsmCOZRP53qOU1lTJf5bQdQnUnzAzNWMTDv8D2YjHJXzhEmowzmJGiCuXfl/P4VQCuAnPjv2XB+vzPuGhH9yJhyx7kQQB2l9Gvb6xsBELCL60zNA2CPeAkopYfBvIwLbdtupvFvj4mNYF/mwabt/I5pgRBCAJwN4Fj8pRKwTjf2sSIAtgrm5+ecpjImCCHZhJBcQsj5hJAKAN8H8G4GzY8AeBzAS5RSRyw5E+YIYDGYUWomhLxNCCnJkPldCaAOQBkh5GB8jscJW4tRTfuEYKwuuUbckEa54+SBxVPt8Nf6n8F5mPdrn0t/27Z+5pzKmHZuADAAzOPIlPm9COALsIv3HgBPA6jIoPn9HMzzuk+wfVfPMQK2dvBrAD8C8O9gXn01IWRoBsxvMIACsNDPIgDTACwFcC9YXLmr5+cb2aOv44TBHpvstJj+fqbmAY+5dLNt62fOqYyZIL4a/ySAajDPys9Y5vOUrvnZz/2/gz1+5gOYASAEIBB/f5fOjxCSDWA+gPmUUtEFntJ46Z4jpXQ92EIc5w1CyJtgHuL/BbsJd+U5PAtAL7BwwUPx1yoJIWcBuJ0Q8gcfY52xa8QL6Sl3nGawi9tOlunvZ2oe8JhLs21bP3NOZUwAACGkH4BlYAsi/0QpNTJlfpTS7ZTSdyilLwL4BwAXgXlVmTC/+8C80f8SbGveb1fO0QKldBtYxgFPY+vK+fF/X7Vt8wrYjXdUF8/PN9Iod5wjaHuEMcNfO3wG52Her30uh23b+plzKmOCEJIDYAXYSvg/UEqPZtL8zFBKowD+BqCcEBLuyvkRQvIA3A32dHE2IeRcQsi5YBd1MP57r66cowcHAfRux1jpnh/f7nPbNvz3TD1/DqRR7jhbARTFH5PMXBr/d9sZnAcAXGx+kRDSH+xxfatt24viC0tmLgXLjd2d6piEkCywXOChAK6hlH5qG/vvYFkF9rGCYIs69vn5OaepjCkiDLaolN3F8zsbLJvgIQCfmX4uBTAs/v97u3iObgwCi9OnOla657c5/u8A23j58X+/QBdfI75JNV1D/jjSgC6FM98yBFboUd0J+/PKU94BFuMz50s+AJYvOdT02j/DmYOZCxZyeDnVMQGoYF5nFECZx9xXgHlW5nzTX8bncmV7zqmfMQF8TzCXHgD2ATiQAfPLiX+u9p+/gxnk6wAUZ+A5HBf/HizMgHN4Ufz3B03bEAArwYxoj668RlK6xtNpML6rP2AZBhEwT+cWsEqgCICxadzHffGfV+Jflhfiv99p2uYaWCuLHot/MZ6yjaUC2IC2aqXb4wagEcBg27ZJxwTwx/ic3oCz2sv8pR4JtvhhrsxqBrC8vefUz5hghRjL4+frZrAFvwPx4/pJV8/P4zN/D87ika48h28BmBsfa0H8PZ8DOKer5xff7qX4Z/ocWN7yW2Dfy9929TWS0rV+Jo3Xt/UHLPb3n2DxpRawXMYr07wP6vKzz7bddQA+js/jIJgBEiX69wLwPFge8SkAawCMdNm355hgxsPv/MbFL77m+AW9AED3jpzTZGMCuCk+xwYwb/4LsFDLD1IdqzPm5/GZvwebUe7Cc/grAB8C+DJ+Dg8BWAiTQe7qcwgWAnoA7IYbAUt9vDXV73NnXCOp/MjOIxKJRJJByIU+iUQiySCkUZZIJJIMQhpliUQiySCkUZZIJJIMQhpliUQiySCkUZZIJJIMQhpliUQiySCkUZZIOhlT09GZptfOjb92f9fNTJKJSKMskUgkGYQ0yhKJRJJBSKMskUgkGYQ0yhKJDUJIkBAyixCylRBymhDSSAjZRAi507ZdDiHkIULIbkJIKyHkC0LIq4SQQV01d8k3H9mjTyIxERdP/38AJgJ4G8DLYMpfJQDKATwR3y4HrGfdOWBqaZ+AdZq4HcCHhJCLKaX7z/T8Jd98pFGWSKzcDWaQ51FK55r/QAgxP1n+B1jXjcso61XHt1kEYDuYdOPMTp6r5FuINMoSiZUbwLpL/If9DzTeBDbeIugGAFUADhFCck2bnQLwAYAfdv5UJd9GpFGWSKwMAROWb/HY5nsA+oAZ3i9ctjFcXpdIPJFGWSJJHd5McxVYKyOJJG1IoyyRWNkJ1mU5RCltddnmCwAnwZpxrjpzU5N8F5ApcRKJlVfAerPdZ/8Dbzcfjy2/AmAUIeTHokEIIX07c5KSby/SU5ZIrDwG4FoA9xFCLgFLi2sBMBzA+QCujG/3OwBjAbxGCHkNbHEvAqAQQBmAzZDZF5J2II2yRGKCUhohhPwQwD0ArgdQAWaUdwF40bRdIyFkbHy7nwD4EYAYgHoA1WBdkCWSlJHdrCUSiSSDkDFliUQiySCkUZZIJJIMQhpliUQiySCkUZZIJJIMQhpliUQiySCkUZZIJJIMQhpliUQiySCkUZZIJJIMQhpliUQiySCkUZZIJJIM4v8DbN0L9S9Dy/0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "clusters = np.array(pd.read_csv('datasets/33k_pbmc/analysis/kmeans/10_clusters/clusters.csv')['Cluster'])\n",
    "counts_downsamplingExp = counts[clusters==2,:]\n",
    "\n",
    "#randomly pick half of the cells\n",
    "group1_size = int(np.ceil(counts_downsamplingExp.shape[0]/2))\n",
    "seed = 42\n",
    "np.random.seed(seed)\n",
    "cell_idx_shuffled = np.random.permutation(np.arange(counts_downsamplingExp.shape[0]))\n",
    "cell_idx_group1 = cell_idx_shuffled[:group1_size]\n",
    "cell_idx_group2 = cell_idx_shuffled[group1_size:]\n",
    "\n",
    "counts_group1 = counts_downsamplingExp[cell_idx_group1,:].astype(int)\n",
    "counts_group2 = counts_downsamplingExp[cell_idx_group2,:].astype(int)\n",
    "\n",
    "#downsample group 1 to 50% depth\n",
    "seed = 42\n",
    "np.random.seed(seed)\n",
    "counts_downsampled_group1 = np.random.binomial(counts_group1,0.5)\n",
    "\n",
    "#concatenate again\n",
    "counts_postdownsampling = np.concatenate((counts_downsampled_group1,counts_group2),axis=0)\n",
    "#remove zero rows and columns\n",
    "zero_genes_idx = np.sum(counts_postdownsampling,axis=0) == 0 \n",
    "zero_cells_idx = np.sum(counts_postdownsampling,axis=1) == 0\n",
    "print('removing', sum(zero_cells_idx),'cells and',sum(zero_genes_idx),'genes with zero counts')\n",
    "counts_postdownsampling = counts_postdownsampling[:,~zero_genes_idx]\n",
    "counts_postdownsampling = counts_postdownsampling[~zero_cells_idx,:]\n",
    "depths_postdownsampling = np.sum(counts_postdownsampling,axis=1)\n",
    "\n",
    "#plot depth per cell\n",
    "with sns.plotting_context('talk'):\n",
    "    plt.plot(depths_postdownsampling,'.')\n",
    "    plt.xlabel('cell')\n",
    "    plt.ylabel('depth')\n",
    "    plt.yscale('log')\n",
    "    sns.despine()\n",
    "\n",
    "#label downsampled vs nondownsampled cells\n",
    "downsampled_idx = np.arange(len(depths_postdownsampling)) < group1_size\n",
    "downsampling_labels=['Original counts', '50% downsampled']\n",
    "\n",
    "#save data as pickle for later use in plotting\n",
    "downsampling_data = dict(counts=counts_postdownsampling,\n",
    "                         depths=depths_postdownsampling,\n",
    "                         downsampled_idx=downsampled_idx,\n",
    "                         downsampling_labels=downsampling_labels)\n",
    "with open('pbmc_downsampling_experiment_data.pickle','wb') as f:\n",
    "    pickle.dump(downsampling_data,f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Apply GLM PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_glmpca(counts_postdownsampling,fam='nb',optimize_nb_theta=False,theta=100,dataset_label = 'pbmc_downsamplingExp_')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Apply other methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "normalized_postdownsampling = normalize_and_scale(counts_postdownsampling,scale=1)\n",
    "log1p_normalized_scaleMedian_postdownsampling = log_transform(counts_postdownsampling,scale_mode='median')\n",
    "log1p_normalized_scale10k_postdownsampling = log_transform(counts_postdownsampling,scale=10000,scale_mode='user_provided_scale')\n",
    "\n",
    "sqrt_normalized_postdownsampling = sqrt_lazy(counts_postdownsampling)\n",
    "pearson_residuals_postdownsampling = pearson_residuals(counts_postdownsampling,theta=100)\n",
    "with open('glmpca_results/pbmc_downsamplingExp_glmpca-py_nb_fixedTheta100_penalty1.pickle','rb') as f:\n",
    "    glmpca_postdownsampling = pickle.load(f)['factors']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compute tSNEs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get PCA of the pearson version for init\n",
    "_,PCAinit = PCA_sklearn(sqrt_normalized_postdownsampling,50,42)\n",
    "# compute tSNEs\n",
    "_ = kobak_tsne(log1p_normalized_scale10k_postdownsampling,name='pbmc_downsamplingExp_log1p_normalized_scale10k',init=PCAinit)\n",
    "_ = kobak_tsne(log1p_normalized_scaleMedian_postdownsampling,name='pbmc_downsamplingExp_log1p_normalized_scaleMedian',init=PCAinit)\n",
    "_ = kobak_tsne(normalized_postdownsampling,name='pbmc_downsamplingExp_normalized',init=PCAinit)\n",
    "_ = kobak_tsne(sqrt_normalized_postdownsampling,name='pbmc_downsamplingExp_sqrt_normalized',init=PCAinit)\n",
    "_ = kobak_tsne(pearson_residuals_postdownsampling,name='pbmc_downsamplingExp_residuals',init=PCAinit)\n",
    "_ = kobak_tsne(glmpca_postdownsampling,name='pbmc_downsamplingExp_glmpca',init=PCAinit,do_pca=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gene selection based on Analytical Pearson Residuals\n",
    "\n",
    "#### Select genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "residual_var = np.var(dataset['residuals_theta100'],axis=0)\n",
    "maxvar_idx = np.argsort(residual_var)[-1000:]\n",
    "counts_maxResidualVar = counts[:,maxvar_idx]\n",
    "genes_maxResidualVar = dataset['genes'][maxvar_idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compute sqrt(x), residuals, glmpca on reduced count matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_glmpca(counts_maxResidualVar,fam='nb',optimize_nb_theta=False,theta=100,dataset_label='pbmc_maxResidualVarSelectionExp_',penalty=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('glmpca_results/pbmc_maxResidualVarSelectionExp_glmpca-py_nb_fixedTheta100_penalty5.pickle','rb') as f:\n",
    "    glmpca_theta100_maxResSelection = pickle.load(f)['factors']\n",
    "residuals_theta100_maxResSelection = pearson_residuals(counts_maxResidualVar,100)\n",
    "sqrt_lazy_maxResSelection =sqrt_lazy(counts_maxResidualVar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compute tSNEs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "_,sqrt_lazy_maxResSelection_init = PCA_sklearn(sqrt_lazy_maxResSelection,50,42)\n",
    "\n",
    "_ = kobak_tsne(sqrt_lazy_maxResSelection,name='pbmc_maxResSelection_sqrt_lazy_initsqrtLazy',do_pca=True,init=sqrt_lazy_maxResSelection_init)\n",
    "_ = kobak_tsne(residuals_theta100_maxResSelection,name='pbmc_maxResSelection_residuals_theta100_initsqrtLazy',do_pca=True,init=sqrt_lazy_maxResSelection_init)\n",
    "_ = kobak_tsne(glmpca_theta100_maxResSelection,name='pbmc_maxResSelection_glmpca_nb_thetaFixed100_initsqrtLazy',do_pca=False,init=sqrt_lazy_maxResSelection_init)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext watermark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Last updated: 2021-05-21T18:34:47.603182+02:00\n",
      "\n",
      "Python implementation: CPython\n",
      "Python version       : 3.8.0\n",
      "IPython version      : 7.21.0\n",
      "\n",
      "Compiler    : GCC 8.3.0\n",
      "OS          : Linux\n",
      "Release     : 3.10.0-957.el7.x86_64\n",
      "Machine     : x86_64\n",
      "Processor   : x86_64\n",
      "CPU cores   : 40\n",
      "Architecture: 64bit\n",
      "\n"
     ]
    }
   ],
   "source": [
    "watermark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "numpy     : 1.20.1\n",
      "rpy2      : 3.4.2\n",
      "scipy     : 1.6.0\n",
      "matplotlib: 3.3.3\n",
      "seaborn   : 0.11.1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "watermark --iversions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}