cjnolet/hlca_lung_gpu_analysis-Copy1 (1).ipynb

## hlca_lung_gpu_analysis-Copy1 (1).ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RAPIDS & Scanpy Single-Cell RNA-seq Workflow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copyright (c) 2020, NVIDIA CORPORATION.\n",
    "\n",
    "Licensed under the Apache License, Version 2.0 (the \"License\") you may not use this file except in compliance with the License. You may obtain a copy of the License at\n",
    "\n",
    "    http://www.apache.org/licenses/LICENSE-2.0 \n",
    "\n",
    "Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an \"AS IS\" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates a single-cell RNA analysis workflow that begins with preprocessing a count matrix of size `(n_gene, n_cell)` and results in a visualization of the clustered cells for further analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For demonstration purposes, we use a dataset of ~70,000 human lung cells from Travaglini et al. 2020 (https://www.biorxiv.org/content/10.1101/742320v2) and label cells using the ACE2 and TMPRSS2 genes. See the README for instructions to download this dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import requirements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scanpy as sc\n",
    "import anndata\n",
    "\n",
    "import time\n",
    "import os\n",
    "\n",
    "import cudf\n",
    "import cupy as cp\n",
    "\n",
    "from cuml.decomposition import PCA\n",
    "from cuml.manifold import TSNE\n",
    "from cuml.cluster import KMeans\n",
    "\n",
    "import rapids_scanpy_funcs\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore', 'Expected ')\n",
    "warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the RAPIDS memory manager on the GPU to control how memory is allocated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rmm\n",
    "\n",
    "rmm.reinitialize(\n",
    "    managed_memory=True, # Allows oversubscription\n",
    "    pool_allocator=False, # default is False\n",
    "    devices=0, # GPU device IDs to register. By default registers only GPU 0.\n",
    ")\n",
    "\n",
    "cp.cuda.set_allocator(rmm.rmm_cupy_allocator)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the cell below, we provide the path to the `.h5ad` file containing the count matrix to analyze. Please see the README for instructions on how to download the dataset we use here.\n",
    "\n",
    "We recommend saving count matrices in the sparse .h5ad format as it is much faster to load than a dense CSV file. To run this notebook using your own dataset, please see the README for instructions to convert your own count matrix into this format. Then, replace the path in the cell below with the path to your generated `.h5ad` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_file = \"../data/krasnow_hlca_10x.sparse.h5ad\"\n",
    "\n",
    "if not os.path.exists(input_file):\n",
    "    print('Downloading import file...')\n",
    "    os.makedirs('../data', exist_ok=True)\n",
    "    wget.download('https://rapids-single-cell-examples.s3.us-east-2.amazonaws.com/krasnow_hlca_10x.sparse.h5ad',\n",
    "                  input_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# marker genes\n",
    "RIBO_GENE_PREFIX = \"RPS\" # Prefix for ribosomal genes to regress out\n",
    "markers = [\"ACE2\", \"TMPRSS2\", \"EPCAM\"] # Marker genes for visualization\n",
    "\n",
    "# filtering cells\n",
    "min_genes_per_cell = 200 # Filter out cells with fewer genes than this expressed \n",
    "max_genes_per_cell = 6000 # Filter out cells with more genes than this expressed \n",
    "\n",
    "# filtering genes\n",
    "n_top_genes = 5000 # Number of highly variable genes to retain\n",
    "\n",
    "# PCA\n",
    "n_components = 50 # Number of principal components to compute\n",
    "\n",
    "# t-SNE\n",
    "tsne_n_pcs = 20 # Number of principal components to use for t-SNE\n",
    "\n",
    "# k-means\n",
    "k = 35 # Number of clusters for k-means\n",
    "\n",
    "# KNN\n",
    "n_neighbors = 15 # Number of nearest neighbors for KNN graph\n",
    "knn_n_pcs = 50 # Number of principal components to use for finding nearest neighbors\n",
    "\n",
    "# UMAP\n",
    "umap_min_dist = 0.3 \n",
    "umap_spread = 1.0\n",
    "\n",
    "# Gene ranking\n",
    "ranking_n_top_genes = 50 # Number of differential genes to compute for each cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load and Prepare Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We load the sparse count matrix from an `h5ad` file using Scanpy. The sparse count matrix will then be placed on the GPU. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_load_start = time.time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 125 ms, sys: 228 ms, total: 354 ms\n",
      "Wall time: 353 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "adata = sc.read(input_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(65662, 26485)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adata.X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.diff(adata.X.indptr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1347, 1713, 1185, ...,  651, 1050, 2218], dtype=int32)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "52985"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(a[a<3000])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We maintain the index of unique genes in our dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 836 ms, sys: 586 ms, total: 1.42 s\n",
      "Wall time: 1.45 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "genes = cudf.Series(adata.var_names)\n",
    "sparse_gpu_array = cp.sparse.csr_matrix(adata.X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verify the shape of the resulting sparse matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(65662, 26485)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sparse_gpu_array.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the number of non-zero values in the matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "126510394"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sparse_gpu_array.nnz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total data load and format time: 1.8720729351043701\n"
     ]
    }
   ],
   "source": [
    "data_load_time = time.time()\n",
    "print(\"Total data load and format time: %s\" % (data_load_time-data_load_start))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "preprocess_start = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We filter the count matrix to remove cells with an extreme number of genes expressed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 535 ms, sys: 304 ms, total: 839 ms\n",
      "Wall time: 838 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "sparse_gpu_array = rapids_scanpy_funcs.filter_cells(sparse_gpu_array, min_genes=min_genes_per_cell, max_genes=max_genes_per_cell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some genes will now have zero expression in all cells. We filter out such genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.03 s, sys: 162 ms, total: 1.19 s\n",
      "Wall time: 1.19 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "sparse_gpu_array, genes = rapids_scanpy_funcs.filter_genes(sparse_gpu_array, genes, min_cells=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The size of our count matrix is now reduced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(65462, 22058)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sparse_gpu_array.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We normalize the count matrix so that the total counts in each cell sum to 1e4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 415 µs, sys: 599 µs, total: 1.01 ms\n",
      "Wall time: 747 µs\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "sparse_gpu_array = rapids_scanpy_funcs.normalize_total(sparse_gpu_array, target_sum=1e4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we log transform the count matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 42.7 ms, sys: 52.4 ms, total: 95.1 ms\n",
      "Wall time: 94.4 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "sparse_gpu_array = sparse_gpu_array.log1p()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select Most Variable Genes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now select the most variable genes in the dataset. However, we first save the 'raw' expression values of the ACE2 and TMPRSS2 genes to use for labeling cells afterward. We will also store the expression of an epithelial marker gene (EPCAM)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 208 ms, sys: 175 ms, total: 383 ms\n",
      "Wall time: 383 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "tmp_norm = sparse_gpu_array.tocsc()\n",
    "marker_genes_raw = {\n",
    "    (\"%s_raw\" % marker): tmp_norm[:, genes[genes == marker].index[0]].todense().ravel()\n",
    "    for marker in markers\n",
    "}\n",
    "\n",
    "del tmp_norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we convert the count matrix to an annData object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 178 ms, sys: 52.9 ms, total: 231 ms\n",
      "Wall time: 229 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "adata = anndata.AnnData(sparse_gpu_array.get())\n",
    "adata.var_names = genes.to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using scanpy, we filter the count matrix to retain only the 5000 most variable genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 977 ms, sys: 26.5 ms, total: 1 s\n",
      "Wall time: 1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "sc.pp.highly_variable_genes(adata, n_top_genes=n_top_genes, flavor=\"cell_ranger\")\n",
    "adata = adata[:, adata.var.highly_variable]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regress out confounding factors (number of counts, ribosomal gene expression)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now perform regression on the count matrix to correct for confounding factors -  for example purposes, we use the number of counts and the expression of ribosomal genes. Many workflows use the expression of mitochondrial genes (named starting with `MT-`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "ribo_genes = adata.var_names.str.startswith(RIBO_GENE_PREFIX)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now calculate the total counts and the percentage of ribosomal counts for each cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 881 ms, sys: 32.3 ms, total: 914 ms\n",
      "Wall time: 913 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "n_counts = adata.X.sum(axis=1)\n",
    "percent_ribo = (adata.X[:,ribo_genes].sum(axis=1) / n_counts).ravel()\n",
    "\n",
    "n_counts = cp.array(n_counts).ravel()\n",
    "percent_ribo = cp.array(percent_ribo).ravel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And perform regression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 653 ms, sys: 114 ms, total: 767 ms\n",
      "Wall time: 766 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "sparse_gpu_array = cp.sparse.csc_matrix(adata.X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 31.8 s, sys: 10.5 s, total: 42.3 s\n",
      "Wall time: 43.2 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "sparse_gpu_array = rapids_scanpy_funcs.regress_out(sparse_gpu_array, n_counts, percent_ribo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we scale the count matrix to obtain a z-score and apply a cutoff value of 10 standard deviations, obtaining the preprocessed count matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 119 ms, sys: 169 ms, total: 289 ms\n",
      "Wall time: 287 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "sparse_gpu_array = rapids_scanpy_funcs.scale(sparse_gpu_array, max_value=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total Preprocessing time: 48.95587778091431\n"
     ]
    }
   ],
   "source": [
    "preprocess_time = time.time()\n",
    "print(\"Total Preprocessing time: %s\" % (preprocess_time-preprocess_start))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster & Visualize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We store the preprocessed count matrix as an AnnData object, which is currently in host memory. We also add the expression levels of the marker genes as observations to the annData object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 199 ms, sys: 92.1 ms, total: 292 ms\n",
      "Wall time: 291 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "genes = adata.var_names\n",
    "adata = anndata.AnnData(sparse_gpu_array.get())\n",
    "adata.var_names = genes\n",
    "\n",
    "for name, data in marker_genes_raw.items():\n",
    "    adata.obs[name] = data.get()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reduce"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use PCA to reduce the dimensionality of the matrix to its top 50 principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.14 s, sys: 1.41 s, total: 2.55 s\n",
      "Wall time: 2.55 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "adata.obsm[\"X_pca\"] = PCA(n_components=n_components, output_type=\"numpy\").fit_transform(adata.X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### UMAP + Density clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also visualize the cells using the UMAP algorithm in Rapids. Before UMAP, we need to construct a k-nearest neighbors graph in which each cell is connected to its nearest neighbors. This can be done conveniently using rapids functionality already integrated into Scanpy.\n",
    "\n",
    "Note that Scanpy uses an approximation to the nearest neighbors on the CPU while the GPU version performs an exact search. While both methods are known to yield useful results, some differences in the resulting visualization and clusters can be observed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4.04 s, sys: 276 ms, total: 4.32 s\n",
      "Wall time: 4.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "sc.pp.neighbors(adata, n_neighbors=n_neighbors, n_pcs=knn_n_pcs, method='rapids')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The UMAP function from Rapids is also integrated into Scanpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: .obsp[\"connectivities\"] have not been computed using umap\n",
      "CPU times: user 246 ms, sys: 229 ms, total: 475 ms\n",
      "Wall time: 474 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "sc.tl.umap(adata, min_dist=umap_min_dist, spread=umap_spread, method='rapids')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we use the Louvain algorithm for graph-based clustering, once again using the `rapids` option in Scanpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 718 ms, sys: 465 ms, total: 1.18 s\n",
      "Wall time: 1.19 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "import pandas as pd\n",
    "from cuml.cluster import HDBSCAN\n",
    "hdbscan = HDBSCAN(min_samples=5, min_cluster_size=30)\n",
    "adata.obs['hdbscan_gpu'] = pd.Categorical(pd.Series(hdbscan.fit_predict(adata.obsm['X_pca'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 23.5 s, sys: 1.2 s, total: 24.7 s\n",
      "Wall time: 29.6 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "import pandas as pd\n",
    "from hdbscan import HDBSCAN as refHDBSCAN\n",
    "hdbscan = refHDBSCAN(min_samples=5, min_cluster_size=30, core_dist_n_jobs=-1)\n",
    "adata.obs['hdbscan_cpu'] = pd.Categorical(pd.Series(hdbscan.fit_predict(adata.obsm['X_pca'])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the cells using the UMAP visualization, and using the Louvain clusters as labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdoAAAEFCAYAAABEjCZVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAACU10lEQVR4nOy9dZwc2Xmv/5yCZpoeJg1oJI1YWkmrZSavGdaOITeOE4fhhh3n3ht0fuHkJnaSG0MSO7ZjhvUaFrwMWmnFDKNhnmYsOr8/ukUr7QpHM5qt5/ORputU1elzqrvrW+c973lfIaXExcXFxcXFZXZQ5roBLi4uLi4uCxlXaF1cXFxcXGYRV2hdXFxcXFxmEVdoXVxcXFxcZhFXaF1cXFxcXGYRV2hdXFxcXFxmEVdoXWYFIUS/EOKeCzxWCiF6qq//QwjxZ7PbOhcXF5erhyu0Li4uLi4us4grtC4uLi4uLrOIK7Qus8k6IcRuIURaCPEVIYQPQAjxO0KIMSHEqBDiI+c4r04I8ZgQIiuEeFoI0VE9Twgh/l4IMVmtc7cQYlV1n18I8bdCiIHqvueEEP7qvq8JIcar5c8IIVaeeKOqqfpTQohHqu+3RQix+HwdE0LcJ4Q4VK3zn6vt/Nnqvg8LIZ4XQvxTdf9BIcTdp517hlldCPFHQoj/usRr7OLiMs9xhdZlNnkv8ADQBawBPiyEeAD4beBeYAlwrnncDwJ/CtQBO4EvVsvvA24DlgIx4H3ATHXf3wAbgJuAOPC7gFPd94PqezUA20+r7wTvB/4YqAGOAp94vU4JIeqArwO/D9QCh6rvezqbgb5qH/4Q+KYQIv569bq4uCxMXKF1mU3+UUo5KqVMAA8D66iI779LKfdKKfPAH53jvEeklM9IKcvAHwA3CiHaARMIA72AkFIekFKOCSEU4CPAr0spR6SUtpTyher5SCk/J6XMVrf/CFgrhIie9n7flFK+LKW0qIjwuvP060Fgn5Tym9Vz/hEYf9Uxk8A/SClNKeVXqIjxm897xVxcXBYcrtC6zCani08BCAEtwNBp5QPnOO/kfillDkgALVLKHwOfBD4FTAgh/k0IEaEyavQBx15dkRBCFUL8hRDimBAiA/RXd9Wdp52vxxl9kJXMHMOvOmZEnpmxY6B6nouLyxsMV2hdrjZjQPtp24vOcczJ/UKIEBVT8CiAlPIfpZQbgJVUTMi/A0wDJeBcc6sfAN5OxUQdBTpPVH2ZfWg7rY3i9O0qrdXyEyw60QcgDwRO29d0GW1xcXGZ57hC63K1+SqVudoVQogAlfnLV/OgEOIWIYSHylztFinlkBBikxBisxBCpyJWJcCWUjrA54C/E0K0VEexNwohvFRMzWUqc7kB4M+vQB8eAVYLId4hhNCAX+ZssWwAfk0IoQshHgKWA9+v7tsJ/ER130bgPVegTS4uLvMUV2hdripSyh8A/wD8mIrj0Y/PcdiXqAhwgoqD0wer5RHg00CSiil2hooTFFQcrPYAW6vn/SWV7/fnq8eOAPuBl65AH6aBh4C/qrZhBbCNiqCfYAsVB6xpKs5V75FSnnDc+t9URt9JKk5YX7rcNrm4uMxfhJv43cXl8qg6Yw0DH5RSPimE+DDws1LKW+a2ZS4uLvMBd0Tr4nIJCCHuF0LEqubpj1OZ873s0bKLi8vCwxVaF5dzIIS4VQiRO9e/6iE3UvFyngbeCrxDSlmcswa7uLjMW1zTsYuLi4uLyyzijmhdXFxcXFxmEW02Kq2rq5OdnZ2zUbWLi4vLguWVV16ZllLWz3U7XK4ssyK0nZ2dbNu2bTaqdnFxcVmwCCHOFSnN5RrHNR27uLi4uLjMIq7Quri4uLi4zCKu0Lq4uLi4uMwirtC6uLi4uLjMIq7Quri4uLi4zCKu0Lq4uLi4uMwirtC6uLi4uLjMIq7Quri8wZCOxLHtuW6Gi8sbhlkJWOHi4jJ7ZEYn8EXCeEIBHMsCIRBCIB2H7MQ0+ckpjGyOzNAYetBPZmQCRVUJtzaSHZ0gN5NA1zVu+Z1fnuuuuLi8IXCF1sXlGmHvdx/lpcd30hkwAdDDIWqXdjN96ChSKHiCPorj0wBo0QihujiTBw6jaip22STdP3iyLnc86+Jy9XCF1sXlGmD3l7/N2M59aOLUT9bM5hjffQAsC6TEtKryqetY6QypdAYA27JxOG2eyKNzz5/87lVtv4vLGxlXaF1c5iEHBtIs74ie3C4mkqgC2vwWKAo4TiXVvGnijddQSiSxisWKmprmWfUpAKpg1fveRdOa3qvVDRcXF1yhdXGZl5SMU8bdlz/17yy65w52P7WdWqXMdP8wvmiQfXoHy6b3Uk4kcQAVwDl3ffGVS1n/wfcgFHE1mu/i4nIartC6uMxD1i+JM7ZjL/GeTnrffj+HfvwiM1Mp+suSHsUh1tZK/NAwOUfBrzgVkT0Hka52Nn30QwjFXWDg4jJXuELr4nKZGPkSo7t203rdWnSf95LqMIslJvcepHXTOka27kT3+2hau5KJ3ftID46SS2Ro9ks2blxMZtCL7vNw/R1rqO3pINjUSPL4AFMHjhBtbyVYX4vq0Qk3N2Lk8q7IurjMMa7QLkBKZZtkqsTeA0lMq1IWCqooqmRRS5hFbeG5beACYOZIH5Zh8vlv7GR94SiOA4e++xihxjp8sQhGJk/LdavJJzO0blqNJxDEGwqgqCrJvgHSI2NkhsaItDSi+rwMPreFQF0tTetW0bppHQADz7xEenAEo1jEnp4Ay2bwqQkUvw8Q9LzpbtIDQ4ipGUL1dTQsX3pWO73h0NW9MC4uLmfhCu0CwTBsdu+fIZk2CAU95PPGSZEFyOUrc357M2n2HUojZaU84INYxEtHR5SaqGcOWn5tUJhOkJucYnr/UUZHZhCpaRQBq/JFTAG6V2dCj+Nr76I82kdybJLMIxMIB4aee5mnk2FurcuTlR68VomgT0Pz6NR0t+OvidF5+03YzZ08uWWQjYvDDOQ12lsaeWX3KHX5FFin5mz1gB9pWXhDARpWLpu7i+Li4nJBuEK7ANi1b5qRsdLJ7WTKOOsYITgprif+AhRKUCiVGZ2cZNP6Wupr/bPd3GuGwnSC0e17CNbVsO9rD4NQiC/tpnt1N7bVzvBLO1FFZU3qgYRk86YoxZkpClMJVAQCiSWgVL+I9672kUlk2HDvrZilEs1rViGlg6qf+RP0psYZeXEb+WAdQ9ND1FszlDPZM45ZdNNG2m/YcBWvhIuLy+XgCu0CYCZROu8xp4vra9E/nCMS0vF639hfC7NQZODFbfQ/+Tydt9/E6M59aMEAoeZGSsksQ2OvYFoOIzMGXa0xjESKoOZw9NAYEVnA71XRVAelpYOaxnqali+mZUk7Rq6APx477Z3OnjttXNVL46rK8pv0cCsHv/0oZduhaMNUWWPlihaa165EUV/L/cnFxWW+8ca+oy4AxmaKRBsC2ONFTFOiCFBVMC3weASGcQEKW2VqqsxUQ4m25jfevF5mZIz+Z15i+sARykIn6FXR/F6GXtqG4tHRfT6yI+PULu1m5kiGdLaIEowiDQst6GcRRRS/jVD9eIMBWjevo23j+jNGrP74xZnmo20t+G+/m+wXP49HQFR3WPfBd+IJBa90911cXGYRV2ivcRpiXiIBnVFdZWKqiNejMDVTMR2rr7Fk0rYdHCT6OUZFu/elaKjz49EX/ojJyOUxCkW2f/ZLGOmKeTba1YYcnyZdBI9ZQAXsfJHwutUYR46ROHSMqYJDS30NjmkRbmkiMzpGTU8n2dEJet98D03rVl6xNq5Z3c4ry3pIDQ3zlt/8BbyuyLq4XHO4QnuNo6oKQb/CtuNJNvfEMco2QsBLexNoqmBR/bluzBLttVQYeGXXNDdsaECIhRvcwLEddn7+a2QGRyoFApCQPj4MgIfKg4ri8SAFTO7cc/LHEtcUSjNJEJCf8eBYDgMDMyxbu+SKiuwJNvz0+654nS4uLlcPV2gXCO+7s4NiycLrVTnal2FlR4hMuWJKdl5lPVbPM7+XyZgUSiZB/8L0Qi5lcjz35//3VIGAZwoN3OafPFmkenTaNqzFMgzMQpH81Ayl6URlp10NvyTBKZdp27iOG2/ZhC/iLptycXE5G1doFwiKIjh6PEtHa5DeJTEWd4bZdyhJLm+Rz1snteFCsB3Ysm2SO25uRVmAIfsmDx47+VqPRajt6eYXl3RhZDJofi+qx0Owvo5wcyO7vvxtBhpWED7+Qyyh4PH7adu4mrGdezBzBTb/6s+4a1VdXFxeF1dorxEM08GjK0xMFgmHdTy6gqKIM4RwZW+MmUSJkbEc6ayBpiq0NgVIpsqMT5Uv6v1KZRgcztK5KHKluzLnBGvC6KEgjmlhpjMUJqcI1tWw6OZNqLp+8jgjlyfW3spaM81kYx2l6STL3nE/TauXs/TBu+ewBy4uLtcSQl7Iuo+LZOPGjXLbtm1XvN43IsWSRf9ghuODhTPKa6I6nYtCNDcGMQybwdE8PZ0RjvRlaG3yMzFVJJEyaG32MzxSYHLm4oQWoLbGw+YNDVeqK/MOKeXJeehk3wADz72MLxpG9XhYdMv15ManiLa3ouiqu5zG5aoghHhFSrlxrtvhcmVxR7TzlMmZArv3JjDOzngGQDJtktyTZHyiQDpjkjcc2puDTOfKTB8q4REKK5fHKJUcVi6vwd+fZmyyiHF2LIvXZhYewuYTpzt71XR3EOtsR0p5UlRdk7CLi8uVwBXaeYjjSLbtSJx7HyClw8h0nnjYz9gk+H0Q8CrkCxYdzUGCfg3HkmiaQqlskEzZjE4WEfLMCFGKAprKa4r5TOo1dixQhKKw8GakXVxc5ho3rce8RNLTGeL0pCvBgODYTJZYSMWjKTTX+mms9eL1Qm2tnyMTORzHIR7xEg15qIl50VQFjy7YeyzJwHQeRYFQoFJpwAeqAo79Gk2oMn0BUadcXFzeGAgheoUQLwohykKI357r9lxphBCfE0JMCiH2nlYWF0I8JoQ4Uv1bc7H1uiPaeYiiKNTW+iiWLXI5k3zBxjAlq9silMo2tg26opHP29ywoY7thxJ87tFjbFgWpyyz3LSqMq9q2Q75vE1Pa4jl3WEGhvJk8w6qAsUyKAI0HazXmL6NhAQxN9GAi8s1hXrHlz4A/DmwCBgEPm4/9YEvXaHqE8CvAe+4QvVdMnvv+OWz+rnqqU9dbj//A/gk8PnTyj4GPCGl/AshxMeq2793MZW6I9p5iqbC2HiJbM4mW7RIZMoUihLDBMeprNXx+gRH+jLcc0ML//hzG2lvCHJdT5ypmRJlwyabMzl6PM3AcIGaiI+Av/JxOxLCQY1YRCMa1lFf41vQ0R5Be62dLi4u846qyH4a6KAShqUD+HS1/LKRUk5KKbcCczqvVBXZs/pZLb9kpJTPUHmYOJ23A/9Zff2fXMJDhnsXnYfMJEo8//I0JcNm/2Aar6YQ8GhY1VRpiqIgBDTV+TEMh4GRHEPJIsNjOV7aPsXW/VP0D+VIJMvYjsTvV+gfzFIqSRQVmhs9qKqgviGAz6ujnMOhtr5Wo73VDcDg4nKN8edA4FVlgWr5QuJq9rNRSjkGUP170UsxXKGdZziOZMv2aQA0VeD3aSiKgu2ApgmEAI9ecWiSCI4m8qTTZTavrsVxYMXSGB3NYQJ+FduujIBLJYdAQK9EiHLAtgU+v0Y6Y2BYzlmjVkWBpYvjc9B7FxeXy2TRRZZfq1xT/XSFdp5xenhhr1ehqzGIqkIoqNLaHGTFshiRSCWoQiym89NvXsyaFbU8t3uK7o4IXq9KXdhDMKAyNpEDIBxSKRYtdFUQCKqoimD18hp0VaFUtFCVM5fxKEolkbyLi8s1x+BFlp8XIcQvCyF2Vv+1XGo9V5gr3s/XYUII0QxQ/Tt5nuPPwnWGmmcIIXjwnraT28Wiydad04SDOl3tYSIRD+GghuOkmZoq0dZUWev58oFpUuky77unA01X+PHWcQJCQwjw+XV0VUH3OPR0hsjlbfJ5i5bmAKPjBaxXaWpNVKeu1nc1u+3i4nJl+DiVucvTzaqFavklIaX8FPCpy2zXleaK9/N1+C7wU8BfVP9+52IrcEe08xy/X+e2G5tZt7qWsllxggqFPPR0RehoPzWH+qlvHcbnVfnzL+6nZDk0hr0AtLcGaKoLUFfrxaMrDI0WmZwuMzicoybqpaGuIqjaafO0i7uiCzpzj4vLQqXqXfxRYACQ1b8fvVJex0KIJiHEMPCbwP8SQgwLIa56nNaqd/FZ/bxcr2MhxJeBF4Fl1b79DBWBvVcIcQS4t7p9cfW6IRivXcplG6+3opBSSh5+foQjwxnqwz7iXh0FQSSssnRxjJGxAooCNVEP2ZzJku4oHo9KNmvQN5hBUxUGhithHk8fUbu4uFw93BCMCxPXdHwNky9YjE8WqK/zE/BrbO6tZf3SGkp5i8HhHEJIDMOhULBoqPOhKDA+WSQS8tA/lKOl0U82b9HSFCIe8zI6XiAQcL8SLi4uLlcS13R8DROv8dLaHDwZNtCrKzTF/XS2h1AUwUiixKreGpob/diWxOPRkBJS6TKN9X4KJZuWpgDhoEbfYAafV+Pm65vmtE8uLi4uCw13+HKNo2kKmlZ5XsoXbECwe/8MNVEPhuGw/3Aav0+lsd7P0b40rU1BJmeKnJgyONKXQdMEqlDYuK5uDnvi4uLisjBxhXYB0docYMv2CbYfTdDbEcYwKsEq6hp8NDVUzMuimr/WMB0KBZNy2WZJ90WH7nRxcXFxuUBcoV1AWJaNZUmWNkcpFC10VSGfc5hKlrBNSSptEA1XIkGNjhdobQnS0uzGMnZxcXGZTdw52gVCMlXm6PEU6YxFTVSjscaLioIEEpMG5ZKFYzscH8wRDerYlsTvVZlOlkhlygyP5ue6Cy4uLi4LEldoFwiZrEHfQBGAfNFGOg52NfmAEDAyViSbM0FC/3COXMGkVLIZnyny0u4pQq638RXBLJf5zOcfIz85PNdNcXG5ogghHhBCHBJCHK1msVlwvEaavIeEEPuEEI4Q4pKWXrl312uYfMEinzcplW1mknkcx2Fpd4TJdAmjYOH3KhgmeHUoGeBYlfN8PpWjE3lWihiO4aBY4PWdI7OAy0UxNDROJpHgLeuCZCaGSE9P8+yRMg+9ZSOK6l5fl6vD//mLR89KH/cnH7vvcgM5qFSiQ90LDANbhRDflVLuv9z2XirxF288q5+JG1+cjTR5e4F3Af/vUit1R7TXMOmMgWU5DI/l0TWdpkY/9XV+MkmLUpkzQiu2NHpprK/ESE6mLKK6zsvbZ5icqiSj9XldIbhUpGOTHu3Hkz5OjUijKoJcWVDIZrit22Hm6G7MYm6um+nyBqAqsmelj6uWXw7XA0ellH1SSgP4byrp4+aEqsie1c9q+SVzrjR5UsoDUspDl1OvK7TXKOlMmYHhHIbpEK/Rqa31s7Q7ysBwjmyhjOM4OA6UDYuSAaMTZSamKikkQwGBrpz66O+7o80NuXgZpNMFxgYGsSyJrgBCIeaTRIM6mgDMEtPH9rBj6465bqrLwme20se1AkOnbQ9Xy+aKayodoCu01yC2LXn+5SlaGn1omoLjCCIhnS//eICWpiC+QCW1HkAyb7Dl8AxFwyKbNzg+kSNXkJwIvFlWHCanT62rdbk4vvRYH1PH9hAL6wit8rBiSwcbsJxT1zRfsGnwlMhNjcxRS13eIMxW+rhzPYnP5U3DTZPnMru8tG0MgEPHMkyliwT9GpYleei2RdTXelm6KELvkiiRsEZTTYDNS2vxezS8HgVVAcOyKZYro9s7NtQzkyhj267QXgxSSr742HHagmWOjBnki5JEujIJblsSadvsGSydPD7kVRACJvuPY5bLc9Vsl4XPbKWPGwbaT9tuA0Yvs87L4WqmybtsXKG9BmlvCyMEbFhbx7LOKB3tYWJRD/EaL0IIujsjdHeE6emMsKo3iqcyNYtH12iM+ZCOxO+tFP7ouXG2H53hwNEUhmlj2w7JVBnLcuawh/OfsZkiUsKm3jiburzYjmR4pjIpnshZCCGI+E/5GqqqQEGia4IvP7JrrprtsvD5OJV0cadzJdLHbQWWCCG6hBAe4CeopI+bK2arn7OCm73nGuSlPVMIW7J5XcMFn9M/mOGxbROsXxxjOl0m6FURUmHVihi5kgW2ZO+BFJoChgW9PWHKpoMiBIvaQvh9FdHIZA38PpVS2UZRIBjw0Deao6MxgKq+cZ7bHv7xAe66cQnl0QPY5TJIiWVbKDjsHjRY11lJP2iYNlv6yqxu1YmFqk88COLL1qPq3rnrgMu85Epk75kNr+Nq2x4E/gFQgc9JKT9xuXVeDrPhdVxNk3cHUAdMAH9IxTnqn4B6IAXslFLef1H1ukJ77WGYNpqqoCgX58BkmDbf+NEAi1uDpNImjfU+0hmDxZ0xWpr8JFMGT24bwy806ms9FCwbHMneo2meOzDJX/3yBtIpk3/93hEeXN9CMlvi+GSBuoiHn3loKbq2sIV297Ek9TEfO44kaPXnUXwRmhlFmiaKL4BTOtuzuGjYCCHQFFCqDyIK4KttJdw8L6eTXOYQN03ewsQV2jcg+w8lGJ0o4tgSy4b1q2qQCOrrfBQKFqGgjqoKnt8zRW6mzOGhDIWyxVi6xKqeGtZ2xDjcnybs96AIUBRlweewnUgUCfg0PJrATI3jWBZmIYMvVotdLuGLNZA+vvd163AcSJZsQrqCVxdEu1biCV71nNku8xhXaBcmC3sI4nJOViyLc8OGek44xX7uB8cYHs2RTpcoGzbFkoWUEiEgkSnTVOPnyHiWjzy4mLuua6KhKUDWsaqjaoVbb7xwE/a1xsTYFMf3bKc4sJNS/w7KiVFKM2M4ZonanjUEalvwRuuQ0gbl7PgvFmCetp55MmWTKVXmvzMjx65SL1xcXOYSNzLUG5RQ0MMDd7UxMV0kY5oMT+RJpQ1WL68hlTEYniygWBIHQSTq4d8+fhOKIjAMm77BLEvqQpQNuOPmJgL+hfk1cmyboUMHUISgvVYDxyGTyjFYqOe6pgD5iSE84Rjl1DSaPwBCUK5G4gJ4ZGeO7nqV5a3+SoECy1s8TKRMQCWRNnATE146jmGieHSk4yAUd8zgMn9ZmHdIlwumsc7PT7ypC4CnXxhjz4EkgaCKBbTUBig2SRIpixPxLCzboa8/hyJY8ObizOBhOur0UwW6B9O0WNVgY5WLBBrasYo5Qk2LKGVTYDsnPbwB3rwuBMDLx4ps6PQhhARFoTFWOagm4jpDXSz54TEyO/aB/SqveK8HLJu6OzajBQOu8LrMK9xvo8tJFndFqIl5UYUgmzF5dtcke48kyGUln/9GZcmcaUpWLa/hrlsb57i1s4tj25jl01cPCDyBKI5ZRug6vlgd+ckhCtNjfP6Rg1i5FINZgQMYr6rr+sV+jk2ZvHC4eEb515+fQjruMqoLwcoXGHv4cTLb9pwtsgBlA2yb6SdeYPx7j1Mcm7z6jXRxeQ1cZyiXM5BSsn3XDEKB5iYv//n9fpxCCFBZ0gG3bKyjufHVkc8WFtKxSfTtP8OL2LRsjk3aNLc2UKOXcUwDPRRBWiZ6tBanWKBczJI2PITIcK7I0U71vxODLcXrJ9q+BM0XvBrduuawTROnWKI8nSC7+yJDzQpBaMUSwks6Z6Vts4XrDLUwcUe0LmcghGAsU+D4SJZ9B9PcuqIeUAn5K4vnDh9L8cLWCbK5V4/bFgZSSlL9h3h2f4Yzo84Jels8RGUKx6r0/choEW+sHlXzYBtFvJFaWlobXlNkbRtk9RdnOJUlPq7Ino2Vy+OYFjMvvEJyx36yew+fsV9KyXTePqOsYEq2j5nYDtgAfh+5fYdJ7Nrnhhe9ApwrfdxC5DXS5P21EOKgEGK3EOJbQojYRdfrjmhdTseybA4PpSlkbUbHC+wfzrCkMYrfq7CsJ8LxwSxjMyWOjWXZvKyOuhoPybSBpoHugTUrGoiGPXPdjUvGLObIjg1gl4qgq1AuYdtwRpY7oRKob6FQMlGtfCXwhKKgqh7MUgGrnAfz7DCLjiMxbZuhGYeO1ih1HT14AqGr17l5irRsUBVK00lSL+8C06T2rhuZeeGVSgoqyz7rnKLp4NdPjRNsRzKRt0mXHJbXV75/u8YN6vyChKGw+abF+Brq8NbHr1q/LoUrMaJ98pl9ZwVyuPO2lZcbyOE2IAd8Xkq56nLqulL83f5PntXP31zxK1e8n0KI+4AfSyktIcRfAkgpf++i6nWF1uV0vv94JWH5DdfV8dmHjxLUVR64uYVEymBpd5RwSCeXt9h3KEE2Y+LxKeTzDo50UMSpG99NG+uJxa4tZ59icorcyDEUXwBV17BKZfK5AgHv2YYfNRRHERKrVEDRPAyMZmlpCKBIC0X34FTT4uWKNj5doL0qmIcnVk+0reeq9Gs+4pQNHOlgpbJY2TzlVAZjZPzMgzQNLOuc55u2ZDRj0xZVUKu2+EzRRtcU/PqZgVwsGzQVtPo4KIKaNcvRgvNz+uNyhbYqsp/mzMw2BeCjV0BsO4HvzQehrYrsOft5BcS2k9fopxDincB7pJQfvJg6XdOxyxm86e5WMlgMj+fZ3BtnSVuYYtEhnTY4cjzNkb4MqbRBMKghBbx8NEk8pvPq2+EL26bmpP2XQyE5Dkj0cBzHBmmWSebPfaO3cwn6+idJZsrYZomaqIYnGAbbwikX8cSbAIHQVMSrI3gpOuGmjlnvz9UkPzjK2A+fYuw7jzL2yBOM/eBJxh55nKlnt5LccwhpOziGgV2qjPQTh/rJH+nHLpWxHRszmYbqddo/ZYEAEfSfrN95ldOYlOAIyJThwHTlM4r41TNE1nAkRxMmWtUaYU0lsCZmmHpuK+WZJNJxKv8Wlmn5mkofdxnMVT8/AvzgYk9yl/e4nIEQguMjWaKKRteiAP0TBTweBQPJxGSZlNegpTHI2GSJhmY/P726lq17Z1jeE2FkuECpfOqGWCha18waW+k4OIXKKLQ8dSrtZmv8tc3gNQGBzw9YJkFdw8qnEKoOioqRGMcBdF2gUglcIQAVQai1C0XTX7Peaw0rXyCzvTKl5QBOyUSrhps0pmYoT87whe8e5KGVfvzLFxNd2s2zaT9vu7EFczpBYWgMaVkowQBOLk9nvU6mbDM0kGVlDHJlm5D3lO3ekRKPJuiIVb5bNf5zjxc8iqA7XrnO0wULXVGI+hQolkk8vw3h8SB8HlRNw9vcgB4L44lGULRr4zv7GlxT6eMug6veTyHEH1D5KX/xYs+9pr9R1zJSSgzDwes9l+vM3HB8MEsuV2ZVc5RoRCUS8uJPG8RiOk0lH2bRJhzx0lTvx+tViNdUTMN33tBEuWwjbYeRkQJtbQEWd8YQ50xhOT+xjMpI66z52Cr9KZVFEYvRpE1bbeVnEzrxECEUAnUtFGbGKtGhZGVO0TRhMm+QK0jaazVCXgUU9Rq6KheGFgzg7+kgcaAfVYDndDO5oqACD6304wDFA8coHjjGDV4PyS0TBNqaUfxeFOmA5QCCABK8KqtqfMh88QyRLRg2jlQJeSvJUF/z1yMAeSphal3gzFud7UjUUhlZKuMA5nTy1M5IiBlTJ1gXJSYMSoOj4NWpveE6PDXRy7xas84gcC5zybxMH3cZXNV+CiF+CngLcLe8BBOIK7RzhCMrI775IrTff3wIx5FIBD6vSveiGMGQRn3Bx6GhHKoNhbyNx2MRDGr4/Rq+09ru82pEI156e+bO2aQSNlIgpSRjZnjl2Dh3LV92QeeWUlN4Yg2UUqfWX6rBGnBM7FKRgFqmb8phKHFKaAEsCbbj4EwMVuZhdBVpWWRKEPFBe8wLMQfjhAVaKPiitVeu0/OE2Kpl/HgEPMMDbKiFaXRqMFH9PuxiCZWK9jlU5qu29uXY1GqQzuTQohEi3R2URscplUro0TB2sYynvgYrk8eaOSWCAU/1O6eAY1W8uDVgJOfQ1hBEeDSkI9EDfozxKVSgKFVGkmV6WgJQMrAARdNwLKtigxYC03QQCnhUBTI5rLxDKTlDyafgOA5KudKXa4CPc+65y3mZPu4yuGr9FEI8APwecLuU8tWp+S4IV2gvgVK+wPFXdoJ0wLLRAgEcAQ2L2og1NyHEa49ZEpkykxNFujvC1MwTZyHHkVi2BCkxbYdYg04kojM2XeDHOydpiHqp8Wh4dYV1q+dX0EDLsfizLf+KlBDMLKa3O8ArA0N0BNtYdRHzoLo/SG6sHz1Ug51LMp4yaSLJiSU+DWENv2aTN1UQGo5tYVgOugrCEdiKD82rgmMjpY2u2JQtgTcchmIWj66AdPDG5tf1u5K86/5lFEeimNk8nnQGLRrmhWMFbl4ZQQ8FsMtlPNEIQtd4q9+HXSyTfGUP0rbRIiFq2taR3L6Xpw6muHdTK2o4RJmKaRrbBrP6tOLxoAS8hAJBjJkZZNlkrKTQisRXV4tjmaB7GD8+RZ0fbNPC49GgVFmWZaGgWxaGCWXTQQqJV1UYStgsq1c4lrBYHNcwq4ExFFWh4f7bUH3z4/f6etx528ovPfnMPrjyXscn08cJIYaBP5RSfvYym3vJ/OaKX/nS3+3/JFx5r+Oz+gn8PuAFHqve21+SUv7CRdXreh1fOONH+0gODJ1zn2U72LYk3tZM+6rec4rtS9sHONg3zYffs2G2m3pRFEsWz7w4jmHYqKqKIRx0VeH5Q1P8/IM9ONJh14EUN15XT2P9/PPW/OrRbzGWTfHRle8n4KmMOvb3p1nReeFmPse2yE0M4Y/VkTp+AMuySeRMGqI6IEjkbVq7OigULcaHR2iOali2Q7FsEw6cGNmLSjQKCdMZA78uCPgULECaDh5doW7lDa/7IHatc3BnH2P9E6ztiBBc1IKi65SnZggt7sApV5yh9Gj4jHPyA8N4aqJooSDpXQcwWtpoaIwy/dxWtEgYRVMoj8+AR0eoCt66OP6WBuxiCQvByHAKZ3qaptYY5ZFx0HW0QAAhJeWZBDM5k3LZoSWkkAeCmgqqwmSqzL5JizsX+0BVoGyRtyWPHTN4x9KKqKoNtdTfeN1V+8zcgBULE3dEewEY+QLJsfHXFFkAPeTHKZZJjY6Tn5wk2tJMc+8SbNPkO08ehelRTMvhhpvXXb2GXwBTiSL7DiRprvcyPF4mmSvjSPjR7jE+9r4VTCfKvLRvkva6IHVx//krnAPe2/POs8ouRmQBFFUj0tKFbZRR/H60QgFVcDKLUSToQdomta3thMN+pBTkk5P41GylAqHib2jHNksYiXHiYRWcynhYFxp4JYruWdAiC/DkD/fxjuvr+OvHp/nAu1tZ2ebFW18xlSteD4r3bOcyX1M9QigIRSGyYgl22aA0MUVs7Qq0cCWgh9GYwrEt9EgYM51BWjaeeAxpOuQCkjU3N2Mm0/jq6yiPjOFf1IoaDLDz6T28mDHo1gt0tsQJNdZhZ3LYhSKL6nWa6lJ4aqJYpRLfOprkzR0qt7RpoKsctwJcv37lgv/MXGYfV2jPQ2psnLEjx06Zrc6FEPi8XjTNg5mt3HjT4xMYhTxmqYw3VWDrUInulV309syvlHLDowUUFNpaIwyPTxELenCE5CN3dzI4UaCjKUBXQyVnqqou/BuOWcrjFAuAQ23klGewPxoj2NCOY5RRfQHKyUlq2rpI9h8EywBp89WnR7l3VRAFgQewkJV5W2lRMqB58YXNF1/LvOfBpYiaGH+4UuVwTqJ4dBTP63tYq95TJtkTYqxHKoE8jJlUZamPR4eciZXJ4musx8oXKAyM4G9r5rqlcQpDo2jBAPt39NPQ2YDPtikcH2TaE+VXfmYxVqlMdvdBKJt44jXIYABfc0Ml84/Hg1Mo8sHFJWShSDCTRYuFWZHOo10b87Iu8xx3He3rYJsWY/sPnVtkVQUtUDWjSkkxlcHMZimi8uMhncGSj+/vLfDOTw6TyFu8e3Ocn3rryqvbgQtg/apabru5iXiNl+uvq0MIQblkUSpJtu2foW8gXz0uNrcNvVo4NnqkFtUXhGoAjrGkQTGdwMgmUT0+pGXir6+ETwy3dBFs6kD1h3nHei+xqI///c08g1acsulgVef5fH4Pmnd+WgSuJFORJr61M4uvsY41i2suuz414MOxHYrDY3jqa/A2VOa4jUQaPR6t+EkAgfYWhKLQef0KGjub8NbFCXZ3cM/1rdiFEukd+8GrowT8CFWhNDlNYWAEx6xkpiqNTaICfTNlsqEYvsYGYuuWX3b7XVzAnaN9XUb2HSIzXolWowcDmPkzHc6EpuLx+/FFw6SHx0BTUVWNaEsjUkrGRxI4ZglV87L65utQtPnhYfxaOI5kKlFk7/4kpiU5ESOgqcHHdWsWrhPP6VilApqv8gAlpaSUmmZ6chKfLBFu7nxNj2G7XET1+jELWfKWzlh/PzmthnprBL9Spr53A4p27YamnGucskF5Oom/tbEyeg34QVHQAv6zzNF2sYSVy1OeTiItG8WrU5qYQYjKiNnf0ogej5HeuQ9/ewuBtmasfIFtzx7iQE7lbRvqCXgE/tamq95Pd452YeIK7euQGh0n0lhP/ys76Vi/hkI6Q3JkDMc00X0+pCKI1Nfz/Sf2EfdJOmo9OJaJoqiouoZQVZqWLiZYc/lP9i7XFtK2OXx8EpmbpLa5jfrGhbekZzYpT87gbTjzmpXGJvE1n5p6sfIFpONQHJkgtHgRiq7jGAbZw/2EejpQvJU58czBY9iFIo5hogV9BBd3UByeqEbtUrAyORSvzsxYCm82iRoJ4YlFia29+iNaV2gXJu4c7eugeXTS45OE6+JMHDmGJxRAAL5YlEAkTDlfIDs5yYoGDR0LadvULGqjobOjslheOqj6wokA5HLhCFVlWU8z0DzXTbkm0SJnJ1uoZPUxUaq/KcXrwcoX8bc2Up5KAhK7UCS8rBurUGTmxVewyyaUyuDzEu3twSqVKY1M4qmJYBdL2GUDxedFCkGsLkzJyGOWTEp4iF3dLrssYFyhPQfFTBZFVcjMJEgPVxKeoyowXlnc7jPKpMcncMoGRcNB1xScWJze9ctRzwjfNr9NxS4Lh1/625f4fw/3YT/1gbluyhXh9DWrVi6P4vOix2NYmTye2hj/9ehx3r0qgJ0v4u9oRbVshKbhravBLpYoDI1gZgsojoMDDI3l6fYNoAaCWPkCqeOTZMoW5YLJ0iYfz49ZLPKZtMe8fGvcx9qAveBiFl4OQoh24PNAE5W4I/8mpfy/c9uqK48Q4nNUIkBNnpa950+Bt1Pp9yTwYSnl6EXV65qOzyQ1Ok45l6dcKJCfSQAgdB1FCCTgGAZC15G2DY5D2Ya61ibaVi58j1KX+YuUkly+TDi08LxkrXwR1ecBRaEwMIKvqf6kEEvHoTAwgrRthKIiVJXy1EzFA3k4x6L8xMl6bKhYmEyTvJT4dQ0su+K0pusEFcmWcYcNb95MvmSxYvHVN/dfCdPxgSeePit93PK7b7/cQA7NQLOUcrsQIgy8ArxDSrn/cuq9HF48/tJZ/byx64bZSJMXkVJmqq9/DVhxsQErXK/jKma5THZqGtMoUyrkMctlUFXC9XU09nTR2NOFZVQiy/ijYXAcEIKGNldkXeYeIcSCFFmoPNxmj/ZTGhnH39pEcWAEKSXliWnMZBpPTZTAolYUr85XHzsCHp3ixBTd3jKEKibokgOFskSe+A0LgerzomgaKUth73CRREly5/tvpyWksLTh2ryWVZH9NJU4wKL699PV8ktGSjkmpdxefZ0FDgCtl9ncS6Yqsmf1s1p+yUgpnwESryrLnLYZ5FQI7QvGFVrAKpc5+txLDO/ex/SxfgozSYxcHmybUqHAVN8ALz6392SqrsJMEj3gp2vTdbSscEXWxeVKY2ayJ197aqJ4ImHQNBRdw7YsZp5/hdJkAiOdpTAwgpFMo/q8tPV28C8/GKQoVFSfBz3gBY8HU9HYPmmRMyRKPELWgMlUGTUaon1FB4viHmprfHz/pVH0aPhkoIxrkFlPH1fN17oe2HKl6rwErmqaPCHEJ4QQQ8AHgf9zsee7c7Rw0rniXGSSWfwelXhARQkE6dm0Dt3jLtNwcblUStNJto+WaQ0rtC2qQ1XP8bz/qjGDr7kBpzoaLRwbAEdipjN4GmrRwiGePZpjWWuIBgqIlhbqF9chLRvHtnHKBnVNIW73J0gULRQp6Lh1HYXRCQKtjbxyJE28Nspn+xQ+9oH6q3AFZpVZTR8nhAgB3wD+56tGelebq5omT0r5B8AfCCF+H/gVKjGQL5h5J7R/9JeP4kh46/29bFp/ddwRFEVh+d23M3LoKKVUCtXvIxKvIZtKU+PzE29tQVEFmiuwLi6XjJSSvY9v5wuPH+e/9hr8+vU+ZGMzP/f+9fi9Kn7fqQfeV8dDBlA8HkoTUxXHRMcGAVYyg1MssbmhlmBDEHXc5lfva6XYN4AeDlaDWki0SBhvYz3sPoQWDvKNPTl+/Ow4k6lB/vAXNrH36Cj/6+du43svjvDWm872eL6GmLX0cUIInYrIflFK+c3Lre8ymat0gF8CHuFaF1qn+iT78I8O8vCPDvKBd6+lq6MWr2f2m9q6rOeM7XjbnE1BuLgsOF7YM8lffa2fZN7iJ1Z4iXkEibFxfuXj32NdV5Sfed9a0mqAR14YYkl7lAc2t5xxvmU5fH1HlnuDAcxEhp39Rcq1QcZHUzy0AYypBIFQkEBtFMVoQo+G0EJBrEIBM5HGshN4omE+s9fiF97ZTG/QZM+ETTwzhfe2NYwnSqzprmFwIs/IdIFNy2rRtGtudm1W0seJSsDnzwIHpJR/dzl1XSGuZpq8JVLKI9XNtwEHL7qO+eZ1vH3XMP/vKzupC2pUYnkLVEUQDOr8xi/ehj7Poyu5uLiczbefHeLX/mEr71gd5Pu70tzYIuiIabyr18POcYvvHDT46C0x/mqbxQNNNilfkFeGLe7b3ML/uK+L//uNg5iW5M8/vILEs1sw8mWOJy2eH7F408owfk2hNupB8fuIrVmOFvBj5Qoofi+5w8cJdrWDEKheD2MzRVLTWdqUPH/8g0neuqGBNDpL28PURbzUxXx8/od9+LwqD92x6KomFZjHXse3AM8Ce6gscwH4uJTy+5dT7+UwS17HJ9PkARNURq4PAsuo9HsA+AUp5chF1TvfhDaVLvCn//AUoDCTt4n5FSL+irj+7Ic2sajNjbLk8vo4xSmM6T0cKa5g9dKrH0bP5UyKZQvLcvj0w0eIywKDxybZuChAOZ2nI6bwlX1lVjeqxEM6e6cc6qIenhs0eOctbRwreRmczHP3+kbefFMbzz28BU9ymv6k5MY2jemipDagIJTK8jtVVdmTcLj9vnVokSBmIoXm9+Opq9w38v3D+Be1sONIija1yL/8aIDfuruR0JJOvvjEIB+6r2tOr5UbGWphMu9Mx7FogI1rW9i9b4JFcZ3aGh83bOxg0/pFKMrCzx4DcPBIEtt2qI372X4sya1rGwgH3AhTF0LmRx/COf4NAEQ+ypQ3jWx7Ew3vmOsppYWPdBymnt9OQniILe8hHvHyd/+9n+Yg7D0wwf4Jg1y2xP+5Pci3diX5jc0BdowbSAQp6SFT9pDxaviEw0+uVlgTt2jI2WyOqkTJMDCe4xu78/QoJtc16Tw9YHF3twfDdhjIChoWNRK1c6wUZXJ9AxybNlm1ogUiIaxsHi0cRPX7cIolVtWCGmvml9+m4w1XlvJ86L4upONgF0uofh9CuebMxi7zlHkntADdHTUcODh58vXmDeea816YTEzl6RvI40jJF37UT0nIs+aqXM5Gmjky33snmaGt7E12sKFximQ5SKMvzcu7jnLLeITQe55Bq1s3101dsKR27eeFXWPU+RW+9Ow4a9Yt4q2b6qlLjROILubXV9QxvKePY30TNAYF3zxY5N5uH0Gfw57xEtd162wKqay6ZQ12NsPvfWWQP3jvYr70zd2844ZGzJ272LQkhj8nKVtFVjVq/Mu2PBt6a1lWq/Jc2s+DS2LUeBUCbU14MwU8KmR2HyTQ0YYWDuJrrKM8naA4OkE4HCJeHz0ZE/n4WI6wUyakOpWy04RWSgmyEhnOzU/rcrHMO9PxCZLpIn//L8/yjgdXcN2ativUsvmNbTs89vQoU7kig1MFagIe3vdAJ2Vb0lL36iVjLiewJrcxsOWLaNmDJMYOENUL2LZFvb8IeCDcCaEmKCXwrfp5fCt/dq6bvCCZfPJFCjNp+hImlgO1zTWMaFGaawOs6Krh27vSvLPXi7+9mft+4bt8eBl0RlW2jNlcv6KOtrBK/1QRvbYGTyzEEr9JQfEwc3iA7aMG9y0JkLAUYo1xxg8OEvUIDCmY6VzGyoDBkL+OockCb76xhR9vGebuG9r58fZxVjdo/Pv3j/HLDy5CmhZCEQhNpzw+ieL18s2daT745qUcPzRCfXcr8YYo6X2HUeMxcvuO8NLxPMvCDjGfQsGCPRMmbVGV2oYQwSWLaei6cvGsXdPxwmTeCu0bDdt2eGXXDH6/yhOvTBDza9ge2HUkyf3XN/PgzW3uk/RpSMemtOvvscZewtDqeGL7FPfV/BjTNtEVm1emu2luaaCjZzX2zF6cchpSRyDURuCuf8PTcstcd2HBMTpd4MWvPEVLwOHwlM3dPV7S3jBaXS2HCjq1hSQrltbzyktH2TbhcFOsxH/sMvjzd7UwNp6lhKAzKDiedtgpa3jvygB2MkWueympRJ7euEKmaNE/kqXTTjBVhgZ/JTBFbN1K9GiYTN4k5BV8/fuHcPx+3nZrO6W9B9g5AxtaPGS9IV7eM05udJoHNzcTaG1gemiKfKiWTKZIaWiYXeMWG3pi7Nw7xsYVdQhNY3GDF2N8GiwbpMS0JcezEk1IYmtX0NbTTDR0+cv/XKFdmLhCO0/YdyjJN58ZoqslyLLWMIZh8eKBBK0xH1sOJQn6Vf7kl9ade3H/NYox/BSFH74frCKT1NF43c8Sue7XEdprJ0jPTxxkZ1+GZcd/mYnxEWJ+yVTRT4562kNphu1lxIrbaVp+O2pxFKw8Sv1GcocfJiCnTtYT/mgOVXU92K80xYlp/v4LO7i/yaIlpPLwoSJ3dHopyIrJNeIFECgeD8/TxHvv7sbOZHC8Pr63fYZ3bKojf7SfrPAS0yWhZd2kdx9ECwXITmeQPj8/2jHFTzywBGNqhi1H07QuqmPFhu6TD6JmNsdU1qKlJYZjWRjTSfLHh1BCQTSPzn//6ChCOrxlZQipKmBLHj2U54E1NeRmkui1tWiZNIl0mZmSpLY2yLgSoiGsooZDlMenmJgqsD+t8D82Rfn2zhQtDUFaexexbtXljW5doV2YuEI7D9h9KMGhY2laG32MjJeIhHR6uiL0D2TZfjhBS62fkF/nLfcsQlWv/VGtU06T+dqNDI5nOJxpYlXNMLrqEPNkAT+qP8zDdd9jfW8zOw+OUmvtJTDxfY7nmrghuoVHJu+krvA8t9bvpRS5DqNUoCFYopw8Ro0nj1TCCD0ARhJUH1gFTEDHYroUoEYvoPe8h8j9X5jrS7EgKY6MowaDlMYmGDg8ilrI49dVfH6NdHsPHfU+fE3157XQfOPpQd59+yIM0+a5PVMs9pYopPO0dzXw348cYuPyWhbX+9FCfixUykIjFvHiFIsowQCJnYcQjoFdKGHnizjAC2kfDVaO3jqdvJQEATwadiCEms8zjQ+/madxeRfPbR/lYH+alYuCqEKweyjPra0KU3nJjpTKDe0eHCHYcizPe5dpHJsx8Dc3cvvbr7/ka+cK7cJkXjpDvZGYThX4xhOD3Nhbx/VrGzkez/L8jkke3z6OqlQsVV1NYYI+lYXiBFk6+J+Q66fOq9LYmMKr2pQtsCzwakUoFukc+H2eedbk9s4Eh7Mt/Nj7e2SSx7jltlv4+d57mX7uD6nVF2NPbkFbfA/enndT2vFX2PkpEjNpIqqCLnKgBdiXaKArlkF3kpiiBv+6nyFw81/N9WVYsPhbK0uqfrB9khtu3URiMoVvYpBXpqG0c4Ce915/hsialsOxkSy9HdEz6nnLja1IKdGQrOupISwNipqNqgruXx4k3hLCTGV48UiKjjDMjKfo6Yjysa8M8Fs3hZnMmPTWKiArQd1Ny6HTU2JRVGPriEFTVCVhS8IBQa3fQUaDtAb8lIpeSuPTyLp63r6uk4G9A7R6LTqui6OZBkG/TXvUYsdIjt46hbd3CX54uEDBBHtqgtvebCMW2Hp/IYQPeAbwUtGNr0spLyo60rXAudLknbbvt4G/BuqllNMXU68rtHPMZ75xjA3dNUTDKoPDWV7YM8WOo0luXl5HyO9h6eIgPV0LZ+2wnRtheu+3EMRIlR1+NLaBe5p3Ml4Ic0PDQOUgvZ61jSqrG9J8Z08MW5Z4IPSP3PUnX8IaeRZVN4lrCRJ5jXwhTsv+z2Id/w5KpIPw3Z+Gh9+CXc4zUxTUReKI2l4M9RgBbw1L7vlPlNAb04vbtBxGp4scHcmyqbeWSHB2l4y9/a7FADTG/RyPhtkz1M/tSxvBdk4ec8Ki9mq7mrRtrNEx9NYmrHSWmngUu2ijR8LooSA1zXFK41Mofh8rGnXyAyN0BTXs6SS/cUOAmUyZPWMWvXEfEzkbRcCzAwZ3dHkBhZ4aDUs6hP2Chw8VeEuv5PkxyfVtWbJ4aaj101WcZMdWhbAwMG2TmGVweMbAQVAXEEwXbY4kJLV+hRvbPTTW+Gi477Y5F9nHP/aJswI53PMXf3BZgRyAMnCXlDJXDcX4nBDiB1LKly6z3ktmau+LZ/WzftWNl9vP/wA+SSX37kmq+Xjv5RJDPC6QMdK1yfBolt6WMIomUDXBdLLMirYobXUBNE1h49qaBSOydm6Eqc900/fpNQRTLxF0kpQdLx/ofJKwVmBVbBjLAdQwaaURxROhKGu4e2MH7+o6yAbP0xS2/X/sfuE7JL79DrTaVfiKh/HFFmEJL6h+PCt/kce/8x/MBG5CSkltUECkh9qAQV3ndUTe/kO0+rUo/ms+cPwFkdyxn89+dRefefgw+/vT/PzfbGF0usBNq+ouSWTLkzPVv9NIx3nN47765ACmdeb+ruYQv/+hVVy/tolC/zAAX3y0j9LIBLqmsPxVo1mhqmhBP4quIZEYM0msTB4zlcZMZ/HUxMj7w9i5POP7BxlNGuTyBpgWAWnTEBA8tNJHyZD8564iAU3wrhV+PB6F54cMDEcyk7dRdIUb23VeGixxQ4NEUxSaKGCOTfF3T0xxU6RAp8fg6IzJ53bmSJcl+yZMdoyZvH2JhzWNGktqVI4lHLBt9Dle618V2bPSx1XLLxlZIVfd1Kv/rvy84wVSFdmz+lktv2TOlSavyt8Dv8sl9tkV2jkiXzR5/pUpNFXBscEwHboXhQkGNX7mnYt515s6aai/ZlN1nYUsJzGtLI6UCAG2hJ7IJN8b2YRXNVGwKZgeMmXwGoNYE1tQnBxyeivYeZ4svYsnHn2U+qkv4o134UxtJ9RxK7LxJraGP0ZRb2X84JPcsraZ/ORBLBvw1iBndnLMXIV3+YdRAm+sKFFH9w/xL98+wv5Dk2RSef74I2u4cVU9fq+GMZPEMczXPPfpHeNMJIpnlGmRSrB9oeuMTuax7XOL7Xvv7EB/jRjBqs9LqLcbM53lg/d14287+zORUvLd7+/DylTu69J20GuijEo/ScfDmOUFRfDiC0cx8wUaQyq9dR68ekXkUoakL2EBoGuC374xQNBbaU9IhZvbPTTGvWRMsA2HLSMWlgOZkkQ3DSYyJl/YXeSW9orBT1ck6xo1Huj2siiisKpRY2WdTl/KJupVKDgKXztQ5mvJOIo+50bCWUsfJ4RQhRA7gUngMSnlGyJNnhDibcCIlHLXpdYx59+KNyovvjyBp+pBnCuaZEsmm9c3URv3LchlPNbQEwREiUDAomCqCCHQsXhP58tMW01EtBlsUyeqZbFtBSNv4y0m8bXcCizm3oDNUOtnUV98L587cjMPrZ6hb6aeJeXvsr6wH8PbQowdmC9/mUWqgdL1dnYcTdMTy3LffW9Fq79uri/BrGOXy1iFEk/vniLuhSeGHD56bxtdMUnHZB9iWpBfvIjdBR8rGj2EBDi2jZUroAV89B+foTmi8Z19OR7o1pH6mSbQb7w4hs+jsrY1wHRREq/x41cVpJR8+YkBPnBP53nbWJ6c4d++vp8RSyeVLvHHH72OpuYzR7Njg9M0+yX+rnaklBw7Nsli3cNEosQyn8BOjJM+luCWDh+JokMxa7Jt3KFoOrx7mYfmkIomJDvGyqxp1HFOztFKNBUSRZugaXN9S2U5zruX+7Bsp/q7qwxY3rPcQ9Cj4QCjGRuPptAU05nKWKRLDosigl2TkqW1gtqGEP/0sXU0dzVegU/xspm19HFSShtYJ4SIAd8SQqySUu693HovkauSJk8IEQD+ALjvcuo5r9AKISJUJn+Pvap8jZRy9+W8+RuZYtlB11QURRALeRE+yOVNamLeuW7aFccxi5Re/lNemOjkpoY+ArpdKQdG8nU0+6co25KIlqeIj6BHJ23pqE4OdXILKBraht/nmce/x5uWvonrEoMMBD6EUtjKEfN2ejb/LH/4cC1/v+ZvIXMcW4kwnlJ45nicmvf+Fq1vAJEFmPzB00zlbWocqAurvL1HRbFTBKRCOW0zmZc8v20Xt3V4sNZ08Z0fjlNO5uhu8OKRFv2GRrDVy63NrUzrEYaPJqmPefna48f40MYYu/ZO8ZM31zOU8XLLmoYz3vs9m2p55Hu7SVoK77urAz0S5olXxrl7QxPZvgH++6lhuut9TOoR3rk+wps/2cfmRkHy6CC7t5RZvbyRAzM2G+scQopCV0hSGpng+3vSpDJFyB6gu6MJr0cjOZiklMmj2g65nEVrSOXty3SG0iaKIpASvKrCmgYFAScfXNNlm7hfIagrWJYEvfIdtCxJzoKYpxL8KepV8KgCw3EAlcaIhl9XMC2Hsi2Jx3zUrOvhbfc2Up5O0LSkdT49HM96+jgpZUoI8RTwADBXQnu10uQtBrqAXdXPuA3YLoS4Xko5fqGVvK7QCiHeC/wDMFmdAP+wlHJrdfd/AG+MO9gV5oUtowgqJrJUzuAn3taFslBcis+FXQQtwE0NR08WHc/W0h6coc4zjZTg08B0BHkzTIkAtYtvQ6lbjbPrb0ELkkolWdlgMlMOs255A69M9rFEO8gz/g+yaekS3n37OMpMCaftDrLHniY5s5Xli+5l5erVc9jxq0zAR4wS3niMsckMzR0NlEfG2TFqcF2zzs6pMre2qISFxeOPH+Cvni/wq9cHyCVtXhg2mCwrtPoi7D98FKsuy9s7JMlRwebFETp7W/ktI4Ovzs/iljpyR/v50+9P8Hvv7EIr5jCTWbb05fiDD64ACaMHBjh2KMF6YxQlFOCD93ZTnpoh1NtG7lAfP/rYap7fPcnffPs4f/fxOzjy/F6Wew2ySVBk5UEsm0pxa9ALAYcX0j46PAq5IwOohQJBjwKGw6KoxhPHDTY1K7SGNaR0UBSFF4Yt3rTEi5TypAjWBSu3O48AX1VkFaBgSQZTNvEmHdsWBL1A0EcgHEbxeMhPpik31LP7wBh33dKGJxTC11x90IjPu4hts5Umrx4wqyLrB+4B/vJy6rxMrkqaPCnlHuDkU6UQoh/YeLFex6+7jrZqj3+TlHJMCHE9FU+sj0spvymE2CGlXH+u89x1tK/P9x8fPvm6tydC56Lwgk6YYCcOkP3qJsA+eXMrmiq6UvEGlRJUFVBD7DNXsyyaZTjp0OqfpGja+CkxaneyeMlK9jf/OYtL32VQvYHhyRzpTB470stP3NWBR7XJfOdeZCnN0UIXDVGFvuaPc8ctG+b2AswxUkoc0+T5x/ewqNaL3yryaJ9BTNhMpQrMmILPvZDhV+9qQNg27+z1kFuykrHBKXpkhsf6yixbXMfyRp2RY+MMeWq5bV0jA8cn2b17lKURh5cHiqxv1Yl54OlClIZSmjvuWEppchqnZBBdtZTyxDRaKEjuWD/hpYsxM1m+8/QAb1obx0okwDklfic4loNlHTWk8wbx+jBqIEDhSD9PHCuxslnjmWNl3rvyzAAnjuNgO5X88CMZi/bYmRGbLBs0FfJlB82r8MqYxZomL987VOLda4IMRNu5eXMHiqZhFUtse/Ygg2oNS5r9rF9x5kj+SnMl1tHOhtexEGIN8J+ASuUj+qqU8k8up87LZTa8js+VJk9K+dnT9vczC0K7R0q5+rTtZuB7VC74h6WU5xzRukL72uzYO8HYeMUJZeWyGB3toTlu0ewjbYP0l9ZAfuCc+1+YXEp3ZJKmJfeQ6nsMW6+hcfXbsA5+npyp8NToMlZ1RhGKQ4d3gB8NLUPWrOfWJSbj/jt59JVJ3qr+DdHmlYTz29g21Yrw1THY9Nu8886lbD+SYElbmPaGheNcdjU4MRq0LIeR6SINTo7csQFEYyOff/ggrTGdW5sEwrZI5y0CHgVVgKIIciWbkK86x6sqeFsbsdI5Qt2LKE4lkOUyvoZaclNJSpPTeKnMjkpe5aFZnTZVoiGcdA6lLo6TTFa86V4Dx6nUdJaVyKODqqB4dFRNR41HCTY3oNdEKY2M8/SgxQM3tXN4KMPS9siVvJQXjBuwYmFyPqF9AfjJ0+dnhRBh4NvALVLKc04oukL72pwYzd53RzPaAlvUfj7s3DBmYj+l7X8D4y8ANqBhSh8y0ITjrceX3cuM00CwcRXRlQ/x31/5Ks3BPP0zgq6aHJtDu0ARjPtvRyRe5snCu3Gygwhsnhnp4N2dr9C0eD0/yj7EUMpDb0cUw7Qpli3+v5+vPBeebk50uXiklEyOzOAc68NJpk6NQh0HEGdZZ0xbovs9oCoE21rQ62oo9A9hZyuJ2c2ZFCIYQGZzJ89xqv/Khk04Hkavi5PvH2akpNDZ4MdO57BteYZ3s207qKpSFVr4zuEy7+z1odZEwIHw8sUY0wlGUxYx3aF58/mnFVJZg7BWSZ3niccu99KdF1doFybnc4b6RSrPlCeRUmaFEA8A7521Vi1g7ripEb9fe0Pe6NVQG2qoDd+iczvw2blRpGMRKExgHf0KiubjbV2HKSeOsKIePCo8PH0/9zVux5t6nu9O3c9DTV+k0LScyYkJdk01sLxF4b9fHudDP1WiEFlDwKfiSPjkNw/x7v/9FL/yruXcHHoKvevtOPlRynojgYA70r0YhBA0ttXx/LM76PBINFXgnFxXK3GcymgyW7bxqgo/7DN421IBPg+KV8dKpAi0tWAXS/g7WkjvPohwJIYi0MIBjEQSpWjgb2pg15jDnYvrCSxqQQ8HMWcK+AIqubKJlSuhA3ZV4E/EATcdMIXgnesihJd0U+gfxt/ejK+xDlkqI0IhAvFTYwTDtPnO8yM8dMfZDqvP753iwRtaULyXnzDA5Y3LeWMdCyHeAfQAe6SUP7qQSt0RrcuVwi5MUTr8FVLDu3hxtJ194xp3+L+Kx8mzpLOZyZHDxGIxvtZ/A52hCZ4a7sQSYdbUDlHXcxctncvYfNNtAOw4kmAyWeL+61twzAKf+cxnidnHCScfo3Pxenre8gn0SPsc9/ja4bl//z6dkbPzs57TsU9AwRMg7FOIrFiKt6EWI5FCCwaqASgilCZnMNJZikOjUDbxd7SSiNTTXudDC1UehuxSmfLkDHapTGl0EjuTrS7NAbX6vkOOF9nQxGK9iAD0UBBfc/3JEemevhSKgOVtIYSmnvOhd3S6MCepKd0R7cLkfKbjfwZWAi8AdwMPSyn/9HyVukLrMps8+NuP8YsrX2K0FGbyyIvYpsH6ukFemF7D+oZhSpG1PLsnT9KMMCx7+cT742y65U08uXOCtd1RkjmLDcviSCODrQQ48s2fJzTxLcLROnY1/DHL1t1FfUSgBmbX8eVa5/BXf0jAU1lCg5QoQpzlzHQSrwfDFyBaE0QPhwl0tiFtGzOVwdtQe/Kw8uQ0+aFRrGIZPRggunZ5ZZ44m0colb/epnpyR44TWtKFMZOkMDxOaWSc/mmDzlgltrEIBpjOG3Su7+XQcJb1m3twiiWMdOX9ClNJjk8UWNLg46t7ivzUm7r54mPH+eC9XViWw5YDM9y8+upHEHOFdmFyvjUlt1GJb/n7VDyx3jHbDXJxOR9f/ZM7uPfW5SyO2zw31g04+NQy3bEEPcEhYtmXCIcCPLA8xy2deRJ6L//7jz/B4vA0mqay9WDFYfCHP/oRmqbR2dZA009u549fuoXki3/F4598F9v+cRPT/zo3DjHXCnLT9ejhEAqgVEeFClC2HI5VIzNRLXfKBuVCGSkl04k8dtnALpVORps6gV0yiG9YQ6SnA5TqaNmR6JEQWiiIXhPFmE4SXlpJi+eti1OzbgXxG9ZDLEp45TI8dXG8tTGaW+Ps6U9T192MUAQIgaJpqB4PTxwqUPSFeWZUEg5UwlF+8N4uikNjaJrCzavrOZ+1z8XlQjnfHK1RjQaClLIg3ogTiy7zjlBAh573cE8P3Pkehy//1+ford/FyqFH+K1nH2BnsoNfXP5jDgy3Elvawh/887NEg91cNxbhk598isHxAlv3T3Hbuopz1GdH3oc1lGRcW8++QgnFLlOjpdhxrJPRz+7iT35m7Rz3eH6yrCuO3bgeoWmky5I9r/SxMlDiyNEZbMqIgBehaORzJbzREOFigdLENC8MWLw56uWoHSaaGaLruqUouoZj26AIHMticqZIQ1szSElm7yECHa3osQiqz4viORWn+YdbRtnQHSbm0TAaW3luRsM+nuPmdWFUIYgLA9/+PRjelXjjMbRQxRwcKGZY192GMzODv/VUkgl/ezNSSsYPDvBXj8/w97/6xl4a5nJlOJ/puACciDIgqETJOFp9LaWUa851nms6dpkLEkee5Pmv/S/81gSPjGzEp1l8e+hGDG87+bLC//7war6/ZYT9fTO0BVP0rliNV1PoaApw3/WtPLblOG/d3IA8/Gn+6MtT3Nu6m2ishlLPr3LdxutZ0Rmb6y5eE0gpsQpFFF1jyxN76QiD7tUId7QwuOMoqt9DqJRDDwXR4lFCnW2oPh9WoURpbAKhKPQdT9DZU4+qCAJtzRQnpvA31mPni5jZHL6meoxECseykIaJlSuQ8EXpG8+zti2Ix8ij+nwomk5xfJLoqqUcHy/w1I5x7OkZckWD29p1Bg0P966Kk4s34pEW//DVA7ypU+W7fRbvvqeHjb215+/wFWQ+m46FECqwjUrc37fMdXtmg3OlyRNC/BHwUWCqetjHpZTfv6h6zyO05wpxdRIp5TkXRrpC6zKXWKbJY9/6Dx5+9jDf7l9DQ10NH36gm4+8uYeXD86w+1gKpETXBf/09UN0NodoiPvoaAiyeUU9m5Uv8eSOCUqTB0l5lpOfOsqAchNve88HectNrXPdvWuKUtlGs8powQBffqKfd66NVpJK5Arkjg8idB1ZLhPs7sBbFz854gQoTUxRGBrDE4tg5grooQBaMIAa8KNHwxjJNFo4xMuHkly/NIZVLDHz7FZi65bjqYtTHBpFCIHhCzDj6Ow6NMPNTYKyA488cZRDOcH1i3wkCpDJG/zy/W0UghHSowlsv4+6zhaaav2v07srz5UQ2r/33nJWIIffKD93uenjEEL8JrARiMwHoU0/+ndn9TN6329ebsCK24Ac8PlXCW1OSvk3l1zvpcxDCCFuBj4gpfzlc+13hdZlvvDk9gm6moN0NlfmArfsn2bzijosy+HX/u9WJlIlNi2L86lvHiYW8nJkJINlwy0dCVZ6tyEck5ZAkuU1I/zM1t/jTTcu4jO/ewN+r5uP41KQUiItG0XXKA5XQsU60qHQN1TJNRsNYWRzRFcuZXjbAY6NF4hSpmdZK0IRhHo6Se0+QKR3MWY6S2FwhERdG7GpIUQ4wvjBAcqWw65yhIdubWZPzsvze6dJZMr8wtuX8K9f2UXfSJ6ldTrLIibfOWLzrjs7kdk8NfEgj+xOUW+m6V7TxYfevOyqX5/LFdqqyJ4rNOFHL0dshRBtVAIVfQL4zbkW2qrInrOfV0BsO4HvXUmhveC7hRBiHfABKutnjwPfvNQ3dXG5Wtx5XSPfeHqQprgfn1dl84o6ADRN4Z9/azPHR7N89/lhjn/1nQxM5HjlUILxmSL/53M7GQ7eyUTCIhJQ8AX8lE2HLftm+O1/3s6nfuP6Oe7ZtYkQAlFNJXciRZ5dKKGoGorPQ3F0AiuZIXe0n0jYR+/0JGoogFkoYgyNkjtwlCfTAd7RkEKPhhC6TmToGGaxRD5T5pt7C7xrhY87ujx885EDPD2psudYio++dTF//V/7uLnOQTp+FoccGv0671uu8sCtHXzwT1/gP397McbQENfVKETWRF+vG/OZ10sfdzkC9A9U8rGGL6OOK8ls9fO1+BUhxP+gYjr/LSll8mJOPl9SgaXATwDvB2aAr1AZBd95iY11cbnqvPv2186c1dEU4hfevhRNU1jcGmFxa8XT+NceWs6BgTR/9vm9fOi+LhxHcnw0R3tDgH/65mE+8fk93LupmU29tW/I4CNXCuk45I4ex9/eSml0AscwQFPI9w2h+H2gqNi5PHY1apRlO0yOVHLpWsksQlXBNoks6cIzk+TnbozyR48lKYsRPnBjLT/pSzDeGmGZP8/BfI7nxgoEGuIcTzs0+hUWxcAcGOK3b4vwzWeG2H0wx6J1fjxjE0Tar8n8xVc8fZwQ4sSc5StCiDsutZ4rzFVJk1flX4A/pRIh9E+BvwU+cjEVnG9EexB4FnirlPIogBDiNy6+nS4u8xNFEXg95w6FubglxB9+ePXJuLeGaePRVVJ5k6a4j/9+YgC1Gm5ww7Kr6zSzUBCKQnTNcgA8NRGkZZM+cAQrmUbVVPTaGESClPtHUPwaU1MFblrfRDlfwihb1Ha3IHSNr2xP4SlkSc7kWVKrkciU+OyjwziOZHm9h2/tzqELKOoewuky7+tVOZQWrLdMiqPjJIbzdGg2NR06u8ZN8g2npWy5tpiN9HE3A28TQjwI+ICIEOK/pJQfuow6L5erlSYPKeXEiddCiE9Tifd/UZxPaN9NZUT7pBDih8B/86qQjC4uCxUpK/9O4KkmQv+pB7oBaIr7+Iv/2seOo0k+9T83cffG5rlo5oLCSKbxtzRipnOUimX2DJboDU7i0zUeO1rmniUB7GyShx+b5EeZCH94R5lkTRM//e5uvvTEAO01GdbWOEwdH+dHhwuYjsLWkRIfWOWnPqTw3eMOi/Q8IcvPA/U226Ykz7yQpNYHPlWwuklnc1Qj1HXNmo6vePq4ahyF3weojmh/e45FFq5SmjyoJNORUo5VN9/JJeTgvSBnKCFEkEqwivcDd1GZFP+WlPLRcx3vOkO5vJF4dOsoX3tykI//5Cqm02U29daSHTlKqGUxSAehvLGSR1wJjHSG3IFjSKA0lUAJeHkqF2O9PU6h5PDwoQIPvXcTS9rC5BQf2w/O8E/fOgT5PH/3Ex3883eO0uB1uKXTx0zeoiuukihKGoMqRQsMCUVLZXENYMNzwwZrWnwUOhaz8fquOev3fPY6rrbvDipCu1C9js9Kk1fdXkfFdNwP/Pxpwnth9V6s17EQIg48BLxPSnnXuY5xhdbljcTBgTRTqTK3rKmnVLbweTV+7i+eIzE+xo7+Ircu0fncJ96F6vXNdVOvKcxMlsSWnTi2BSUTiUSJhLAzeWwheHzazwffVlnKk5rO8rUXxrl/eYiGeIDnXh7gma3DfHiFTq5sEdAFW0ZtHryrC9Xr49DhSTqCkmRJ0rK0hUBrE1rg6i7lORfzeR2ty6VzvnW08dc7WUqZOFe5K7Qub1S+9tQgpuXw8sFpvIpgUfEVDk043NhuEWpfxlvedidCcZcGXSh2qczMS9uxU9nKNqB6PYR7F6NHQnhrawAwZpJo4RB2qUR/RpLKGixtj2Dt2YsxOU2qYBH2KJQchVBIJ7y6F9XnxRuPIc6VBGGOcIV2YXK+X/w0MAxUA5eeMT8rge7ZaJSLy7XKiVRrNWGdJa1h/uQ/CjQHh3hsRGV55gix9AGaIhaL3/RRFN98WSkxf1F9XhruuBHHtJjZshNmEghdwymV8HadyrTkqQqu4tFZcnqI6huvw7EsvAeOUp6eJmhJGu+6qeKt7OJylTif0P4TFfv088CXgeekG2nbxeW8vGlzJYLULWua2HNMZ9vRUTydDTzoG2Yqq6K+8AO673JTOl8oiq5Rf8ulDfQUTSO+uvcKt8jF5cJ5XaGVUv56NZHAHcBPAv8khHgU+Bcp5fGr0D4Xl2uaN93QQldzCFVRWdoW5P3fsNlYn2LZSIobaye5be01uojExcXlgjnvZFF1BPukEGIHlaU+fwocoeJa7fI6FMuWG6rvDU57Q5D2hiBL28Pc/As/4Hfe3cm//Gic3dMZbrzfOn8FLi4u1zzniwwVBN4OvA+opxJ28Top5dBVaNs1ybP7X2IqZ9A3UubHLyj87FuX8q7b2s9/osuCpqMpxNf+7A5a6gLEG+o5OJDi9z+9m//+oyh1Ma/7QObisoA53697ksro9ctU0uNJYJMQYhOAlNKNd3waUkos1cGRNu+94Tp++5212LYz181ymScIIWio8fKtZ4bwaQq5skmuaDGdKbN8UZRs0aQ+5uOLTx7khpX1LG5wo025zB+EEP1Alorzt7UQvaPPlSavWv6rwK9QcQx+REr5uxdV73mW9/x79eXpB53wPJZSynPGe3SX97i4vDbP75miNuLhr7+0j62Hk+z87IMkswZPvDKOlJJjM9MoQvCxh9yk4280rsTynjf/z+GzAjk88g9tVyJNXj+wUUo5fbl1XQlS/+o/q5+xXyjORpq8O4E/AN4spSwLIRqklJMXU+/5FpDtBfYA+6r/9gI/Bv7ktUTWxcUFvvnMEJ975OhZ5UeGs2zZP83e42k+87Ebed+dHSSzBsOijOUHzSuI+wKsWxbDcVxriMvFURXZT1OJAyyqfz9dLV8wVEX2rH5Wyy8ZKeUzwKvjQ/wi8BdSynL1mIsSWTi/0IaopEUKnfZ6I/ADIcRPXOybubi8UXjXbe185M09J7ellCQyZZLZMr/yrqWUDZOjw1merMvS+cRT/PTD2/n64/3827ePsqghwtY9Kb76/NlC7eJyHl4vfdzlIoFHhRCvCCF+7grUdznMZj9fzVLgViHEFiHE0yemTi+G8y3v+eNzlVcjRj1OJcmAi8s1w5d2HeWVvQke2NzMvqNptvZPs7Q1Qn2rl4m+Mjd1NvHE1lGWtUe5e2MTnU1XJqhErmhxYCDDzavrOTacRQiV23//cQordWrzgvRAkaW3dPDc2ARjM0XGEmUMYfPDLSMEfRq3rm28Iu1wWfDMZvq4m6WUo0KIBuAxIcTB6ghwLriaafI0oAa4AdgEfFUI0X0xMSUuydVRSpkQbhJOl2sMy7GZMPNE2hSOywyBHri1uxbTkBQMg7plGk/0DfOdl0dQXhzhF/9uK7/+7mVEwx5Cfp3/+VAvj7w4wptvbD2r7j19KRxHsran5uT28kURNK1iNAoHdG5eXQ/A4rYw//TwITKWSU1WpzHoY29imm++NEitz8u/f/8YPo+COgpD/SZP7pjk19+9hF9/aAWqOn/CBbrMS2YtfZyUcrT6d1II8S3gemCuhPaqpcmjEh3xm1VhfVkI4VBJOjB1oRVc0q9WCHEXcFEZ5l1c5pKxmSL/ue8wlnSwHQfbdpASVFUhENQJ+DU0j8Li3hCLFgVZuz7Gg3c3UTAsvvJEP//fF/bw11/ax7qqkAI8vu1UAo/V3bGTIgvg1RW+9Hj/OdtydDjD954d5oEljbyttpFO08u2f34TdlGy7UCCo8M5SmWHo30FpICVvRG+/PJxVv7Cw3z6u0dm7Rq5LAg+TiVd3Olcdvo4IURQCBE+8Rq4j0tIF3cFmZV+vgbfppK1DiHEUsBDJTzxBXM+r+M9nOlxDBAHRoH/IaU8eK7zXK9jl/lGwbL48tARSiUTn08/6WikVAPK2/YJ8RUIIfnxIzMUCja5lEn/QAHDdFjZGaNxkZd02qBjuZ/WRT7e3NPO2ob6C27HF5/p43BfjolUgX/4xQ1MpQ0iAZ1oyAPA//jEczz68jiNNT40HX73d3o5vCtP1jA4MJwiIn0IU+ETH11LR1PoilybgXyWA+kkRdugxR9kbU09PtVd1zsXzFevYyFEN/Ct6qYGfElK+YnLqfNymSWv43OlyfsC8DkqqfIMKmkCf3xR9Z5HaF89NJfAjJQy/3qVukLrMp8wHYenxocYLBaQQKFoEAx6Tu5/ZWuCUsnhxptrKZUsCnkLHMHEgMHuQynifi+hgIaRgaJtMjCSR9NhYLzI5hW1fOF3bkFRFIL+84vTWC5PMSPpbgnx9acGeddt7SjKqVmYisOUwcsHpvnnbx1mY28tH7ink5cPTvP++zvZnp7EyEHc6+H4kRLRkM71q+JEvV507cIMVIeG0xSCBobtcCSXpGxK0MBx4EQimybNS8YyeEtrN2GP5/UrdLliuNl7Fibnc4YauFoNcXGZLUzHYaxcwLFB0SEY9FAoWBw8kGXduhgtTX6SSYOwqqIEBOWS5OCBHJs21qD4QajQEQzz+AvjpIoG+49k8PtUPnRvF021fh55aZRE1uCX3rH0vG2Je3wcKVVSvr3njlN+G7mCyaGhLENTeRQhUISguzVMvmjxwwPDhHUPQhN85t8G6F4Z4B0bF/HYtjE8PthanKBGeAgXvdQ2e9EbHGK6F5+qsTZWx7Fcmqim84M9w+w7nCXUoNDY7Ec7MSA+qfMSx5FICeOUAfjqyDECCHqjcdbH3bjMLi6XwkUnfr8Q3BGty3xiulzkO6P9J7cNwwIEigKapmKkLBxF0F7jJ2Eb7N2ZwWupdK8LYhUdXnw+QXunnwP7sjiORJcK2/Ym6GmKUB/z4dgS6cCanhrWLK4hkzd5/z2dF9VGy3L4t4eP8EvvXIZlOSedqAD+7PN7aOzSGZ8uExE6O/en+MhHO5k+YvLjlye5aXUtn/jqXlauiLBuSQ3+JoVgUCHq1QnpHvoSWVRV4NFVTMtGKKBo4qTZHAccqqNZh3N6btSrHt7S3o3i+kDOKu6IdmHiTsS4LHiypnHy9UsvJAiGVDo7/ITDHnDAE6v8DJ7eOcWKpSGWL4+AF6QlmU4b3HhbLVu2zRCIKhw6mOf22+pYvCbI0QM5tJjDo98bJ1+QbOyt5cEbWgBI5ww09Uxz8u5jSdYsruFcHBrKcNOqyohRe5UJuLMpxFvWtfK/PreT93+kEwUFphQ2Latj/aoa6nU/92xs5rnhcfxxBWELhnN5LOnQn8gR82j4PDBjOui6CiqVIHonRFUBxQGkxEGCU5m7TiVNDhzIsHpNFELw3eE+3tG++Ep9LC4ubxjctQIuZ2CViuRnxpna+xJTe1+s/Nu/jVImgZQXHqnIsUzscnEWW3rhFK1KlhzHkay7LsLqNVE0j8CwbFCgkDYxTYfeJWFGpgxQIZMw+MHD42geyXS6QLxGZ9PmOD/90x1Y0kFoUNuko3sEb3lPM5Z0+Msv7uNrTw5g25LpdGWuFSrBKk7QP5ajWD47a8/KrhjrlpxbhD90XxcBn0pLPEBHIEQgpnBdVx37+zJ0RiIIXfAPDx/AUSVhXUdYgqBXpd7rZ1EsgK2BhYomFU5+girY0q68PiG41WYqisL0dJlQWOWpJ6YZHChUr6N5uR+Fi8sbEtd07AJAuVigMHYcq5A59wGqhj/eiGOZBGqbEaqGVcrjCcU4saRaSkk5l2JkcIzBjM5tG9pRPb6r2ItzkzXLfHW4D8eRzEyXqW841aaxsSKODTVxnbqQh6HJIu11Pmo8XnYcTJEumNgq1Nbo1Nf6ePrZKVpafQT8CqmURSCosntniu7uIF/90gixGo0lrTF+5cFl3Ht9C3uOpSgaFtcvrwPg8FCG1rrABTlOvZq+0Rwv7J3iQ/d1AVAoWfxoaIgmI8RkpsjKFRF2TyeYnCqzsiOCX2g82z9GYswmEFJp6fFiAY7toKhK1fNaoCDOeOS2bZvjfUVa2/z4/eoZTlI/07X8Uj8GlwvANR0vTNwRrQuJgUNkju06Q2QdKtbFE3gicUqWgqJ7MYo5cmP92JZBcXqM7Gg/R48OMb1vC9mBg0RkmlWxNLnpMZx5MAoqO5WeSClJJEws69TDpa4rtLb5CQQ0Dh7J0uYJIIUgZZisWBZl06paspMWXo+gYJrccFMNSxYHmZw06OkJEgxqlIuSsN/De+5YRFNdgANH0rz/T55nf3+KH2wZIR72nny/pe2RSxJZgO6W0EmRBQj4NFb761CaHFqXenjmpSmefWGa+1e2cHAqw/58gtbmIJ0rAgSiKl5UsgkDiYJtOwiUisezUhnt4wBOJctQW5ufocGKRWKwv0AhVxmFD+Wzl9R2F5c3Mq7QvsHJT49gZ18dQ7vyxVBPbAgFM5/j6a39JCanySRTDCVMbMOgmE1SSowRyg/iOA6mVTFOjk4XMBLjmMXXXQl2VTjRJikly3rDaFrVoceBujovKuBXVHSvgiegEJAaI6Ml0qbB8aksXp8gNWWjCYXhoRKKECxbHsbr06iJa7z3ve3EmzQiPi933l3Pm29rZsPSGMvaI/zuB1bS03ZmGMeDA2nGZi7frP5fjx6npy3MKy+naPAGSeglMhmTx0ZHMR0bS4GAEAR8KqoOxbSDP6CjqqC+6qcvpYODU5mvVRR8fpWlyypuyfV1HnTvKauFyxsTIURMCPF1IcRBIcQBIcSNc92mK40Q4nNCiEkhxN7Tyr4ihNhZ/dcvhNh5sfW6zlDzELNUIHV0Fyga8SVrUHXv+U+6BLJjA5jFApX1HWffQH+0O8/9a4Kk8haxYJlVTYJ/+cEYb14fosYPpalTI+AToQFHEiaNMQ0pAa8fbzg2K22/GMqOxdBgAVUVtLT6AbAsG9OQ+AMaNqALhd72ECMzRTwKNLd4GOwrUc5IbEuyrDdEtmzR0RJgdLzI+ES5MhL2q+hhOLA1T0x6qS36CMW99DflsYXky48eP2MUCtAY9+PzXN4zruNIlraH+dfvHGbzmjg/eH6MqOalVvHR4QkxWiiQG7YYCpfpVCKgK8RCKoPlAl4EjgpFy0ZTVHBAVdWT9Z6+rhfAH6y01Qu0BK5MkAyX2ePP/3nwrEAOH/+lRZedJg/4v8APpZTvEUJ4ODuo/1UltvPps/qZWnf75fbzP4BPAp8/USClfN+J10KIvwXSF1upO6KdZ0gpSQ9UA245Fom+/Uzte5l84qIzM70mqYED5CaGKc2MYRdSpybgTsMBbl5aEaVAVRQikQA9bRGaIipBv0YqbzGTs884rzGmUzQkrXGd2u7VV6zNl0NrIEz7ogChsMb+fZWHAyHg2adnKOUtgkCuZKJIga5LUpMWu3dl6W4Osajbx923N9HuC5FKmnh9Ck2tARYvDqEgcBwo5S1uuznOhvuixHs0FjUH6W2NsmN0+pyOTzVhD35v5Rn3vx49ftH92XU0ycHBDLURL7mCRdGy+fn7l/KuW9pZ1VXDjQ1NrI7V0t0cZG20FjNoIjTJRKmEV1dACEwbNNRzLOc5+4FLURV0VWFtrA7tHN8Vl/lDVWTPSh9XLb9khBAR4DbgswBSSkNKmbq81l46VZE9q5/V8kvmNdLkAVCN7/9e4MsXW6/7q5lnFGbGkGb5VIFZAmlTGD3G1MGtyGrowImZIgf7UxdVdymbZmrvFsxsiuToEGWzKgLOKbG0gHzJxjQddsxEQNGwqks/vH4v9y1XiAQEQQ/Eghp7hsoMJyxKpiSZs5EIogGVhw/qKKp6znZcbfLVeeJIRGfFyghQGcHdc38DvqDGgYE8QY9KCRt/WGfp8hC3bKzl+GQO0+NgOBbHihn27klj24AFmgdMS+ILqLy8NUEya7L/aJpkymB7cpoP3NrJo/tHeH7fBK8cnDnpgfxqXj3avRDW9tSwojOKz6PSWh/g+0+O8/yeSaYnDIRX8oWnj/LtZ4cIh3TGyTJSzJO2DXQvNHl8rIzUVm5Nr7JnOY5TWTp0unN59bUEpowyLvOe2Uof100liP6/CyF2CCE+U415PFdczTR5J7gVmJBSXnTA8QUptFJKPj/+HXJWgZSVQUrJvvy1kduzMD7AaaF6Kq9VHeH1gmXzr//xBGYxz5a9I7zjfz3DTb/0w/PW6dgWtmnw6FN72DdUmRvUFYktKzdX+7QE4xpQdFT+49ks3YEMOBYBT6UddrmARzmzdXcsDzCesjg2YTCRU7AVD9HFq/nIe+aP42TM6yOq6gSFipQS06zENd66pZIXo609QFbaFJE4gGE7TJpl6po9FMo2ecdGAZYuCxGtrkMt5R00vRK8IhTWkQLiDR6ylklHr5+XcpMU8ja33VnHP33rEI9uH7msPhwdzjKROHNe97k9Uzx0xyJ+5/0reHHfNIeHsvzw+THqF3mwQxZHxjIU8g7SUKgJ6JQtSJkG26amK798BwzTxrZPi/tcXVd7EgViKHT4Q9zV2HZZfXC5KsxW+jgNuA74FynleiAPfOwy67wcrmaavBO8n0sYzcICnKMt2wbNL9+OgsL/G/sq98Ru5L/GHyZJxWTYTjMxb4hP9/4JSwOdc9vYc6HqYJ/uqSvBtrBsiQq85/oQqWN7uLO3lf2ffytTqRIAmbxJ2K8iXmXas8pFbKOEHgixuTeOKKV4dE+ejV1eZnIOS5oqcWzTJaiJR3AKWWIBwQduieL1aKAIVJ8fu5ABy8HwhAj5PVjZZKVtQmXjkjC2r46SGqRj0YUH2L+a3NXUxg+OD/DsC1MoimDJ0jBr11WclE6fk1SAnO1UtMYj6A1FGDGKHDucp6XdT0HYlIsWNXVezLKNgUN7m5+gX8cRDkKAIyWFQpnJmTIv70gQ1XR2DiR4+u+m+Pm3LTkjy8+FUhP2oOuC4/kMXcHKqNysRpD64y/uYnKiRN60uO/ORo4eyfGrd6xkV3qKCbtAuWyTQyHqUck5NooNXgRlReLxvr7Vod0X5L7m2bx3uVxhZit93DAwLKXcUt3+OnMrtFczTR5CCA14F7DhUs5fcCPaL059DwAHh2OlAT41/mWSZPBmYsiiyhBj7Cke4Qs7X5rjlp4bb+xcQiURr9ouTQ1TTk8TsivTCU/tnCA7egwzX3mgMIt5chMDlJJTmPkMw9NlBsemcYTgvuvqiIe0kyILsHeoVHWIkWg4HB4p4cECRcExSqB6CDR18sxggEBNPZGu5fhqW6jtvQ5fvJFoNEDnPBVZqMQY3tBUz213NHDLbfU0NvnIZi2mJkpnHPfiiwlmpspk8iZaEZ7dP0VUqCxZGsTrrzoFeVUyaQNdV9E0hfpGH1JIJifK7HolzZ5daepiXiwc7r67gZ98Zyc33lZDtEvQ03Fua1uuYDKTfm3TbG3Ui+IRdAUjjBULPDM1ypK1AV46PsEN6+N0LAlgmw47+hL40NhTmKIx4qfB76Mp7kOakHZsbAOyRZPKIp9XYZ/yNvcJhQ+09bAqVnvpF91lLpiV9HFSynFgSAixrFp0N7D/cuq8TK5mmjyAe4CDUsrhSzl5wQntu+ruoVfrREFgYVeWLABlTw4pq3KlwBfFv/P4jku6ZrOGlJJy4kSO0zOlVdG8fP3lHAgFEFg2JEeOMTowxMzoMDe2lPDHm1F0D0Y+A45NsKEdf7yRFwZUokEPNdEoWihOoLapEilf9VRG0ELlllX1mKUTS3FUaiOVoA6WVaakRQjUt+CLxPjAPZ14I3EURcM2S+RzeZ48KvBfAwHne6Nx4tqph4v6Bj/1jZV+mmZlnnrzjXGmZsrMTJc5MpqjpcVHWtoYgHSgry8HiqAm6kFRKn5kKiAUweLmEDffVse6dTUM9hVYvjJEMeHw3w8Po9saPbVR/unJA7znL5+ibzR3RttyRYtU7tzzuCfYk5pmolTAI2D6oIVZjVB1d28Lq3ojHJ/MMzZVpKsnwEy2zP7MDFOFMkXLpsHnxYcg5NUIBDSkc2rZE0B/X570tEFEqCwJxXhLaxd+XafFP5fTcC4XS9W7+KPAAJWp9QHgo1fI6/hXgS8KIXZTSRk3m/Ohr0vVu/isfl6u13E1Td6LwDIhxLAQ4mequ36CSzQbwwKNDDVjprh398/Sb1TmxQJ4KGHgAB40DCyMnUv4xnt/jvtrb5mzdp6OY9skju9DlvJkTBVvIIhXlk86KgUaFyGEwCoVKCcnyRZNNF0lEPBhl8uUfbXEw150XwCEglXO4w3H0XwBBocmqA8KHMtE4qD7QxSnx8mUTLb3ldnQJpGan6BXoupezFIJzeslnytVRthGhmCshuMDk9RG/DQtXkIpNYW0bfy1TXN85S6eL/QfxDjte3/saJ7FPRVB2fJiAn9A5fDBLJmMxfRUmf/5W0vweCvPpCfCBHuoJKbEAo8paKzxMVoqUpqxUQKCbM4kFvVSKFkUiw6/eX3FA/tvv7eP99/eyUzGYH372aPFjGkQ1vST0bYA9qRnyJkG+7MpoLIExyg7bG5oYCCbZRlxHts1ilG0CXep7Nuf5brro7THA/QN5HjkmXGa23SuW1NLMKRwQs5Nw6ZYcoiEdOp0D03BIOtj9aiKQBUL7hn8msCNDLUwWZBCC/C3w//Op0e/zqR9Tk9tpAOt1iJuGn0fD94bpctcyt9+tp+fflMPd25ovsqthezkCKXJyvSCIzRQVBRsvDWNGJkZhKoRqK20q5xLYRtlsC2ccsV6oviC1HStxCpm8YRilNMzaIEwVjHPN54b46E7FyGEgmUUEUJBUVRSw0fBNJksSo5OwrImnWcP5DgwpfKW60IU8iVu/v/bu/PwOK/60OPf866zama0WJIl2ZL3LY6zr2RzFpoEQqApkBYoBErTS9gent57H2h72wIFWm4L96GktNCy3BKWhNxCQgIhq0niJCRxIjveZFuy9m00+8y7nfvHyLFjO4sjy7LG5/NPpHdmjs6ZWPrNWd7fb1W8mjCoXALdINGxBOn7mNEEv960k6vfsvIVQWG+GMzn+PnQfrpfyNDWHqap0X55fadU8hkeKpPLuviFgLMvSLFtZx7hS85YmaJg+IyNlgkJjfYFESrSJ6wbpH0XrQye8DFtnYGBEpGIQSbt0L4oyjs6OglpOmm3QsH36IzGX1FgfaJSJkDSaIUQQrA9m2ZfIUfBrVAo+bjWK39XpSdZEovTEArjFAN+8/QwiToTzYIXtmZwKgFGI+RLDu0dEZZ2xqvrJObBNqYmKpzb0cTCSIyl8cTsv/HKa1KBtjbVbKAF2FPazyXPv58iZfysjQg5aJZEQ7DQaOaM2Cq+vuyz7Crtwx9cwJWfeJBoyOAv/ng9H33bsiOqqMym4tQ4hf5DTo2bIdB1KBexG1owQhE0TceIxMkN9GDH6/EqJZz8FNKp7jMKw0LTdX7Z7fKe69bhYfLElv1cuLYBJ5cm2tROYWIIAQSOQzk9DKaNp0WZGB2ntTVJd88Ea1qt6hL1dAAVevW0rp1qZvu+NCtaTOILu9BNGyefQbfD6Ob8Kw7+u4kRns9WP4j5bsDouIOh83IuZBvY3VtkaLjC+eel2L+nSHNbCMPQ8GTAAtMiKz0mxxzqUhaVgkexEtCUtMiWPSIhozr1BRaFIyyJJXmse4QNK1M02mEsTSdhmvQX8uzKT5H3PMrpgLbmMBnXYcJzaNItfE1iawZ51yE3vcLheQFDg2U6GiKc3tpA0rLZl8+SDGw27x6ngIuLz67tBfpHSqxZFadrRYS6hIVJ9TauGBqr6xs4LdF44t985ahUoK1NNR1oD/hK33f40sC/AnBR9Ay+2PVJTou/fpHuE22s+4lXfUyPJbDiDfjlAl4xh26F0UydSmYSM5Yg8Hx0O4RXyiNCdTyy3eGtZ6X4zVN9rGqPsbC5DqFpvLhtgD0jZd6yDKp5KARoOr7nUXKqwaFQgbhN9WbRIACh8WRPBSfWxtUrJUhJrGURQtOpZNPYdcd+ivZksXVynGczYzhAetLB9XwWLKgm6tjTk0fTBJ1dUZ58fJLzz6s/+EIdXDdASjANQdK0yJRchA8rGxLsmswhQgEJYWAYGlnP510dS/jaA9u46awuknUWJd/lsdEh8oF3ROIIE3A5MPkUdMXiuIGk4HvEdZ1ndkyyuCNMVyJJIAP2dZc5bUmS2/7P0yQtk4/dsoLeYoaKGTA0UqI+ZTM0UCae0DEtjVSqmm3sxtZO6kPhE/JeK69PBdradEoEWqiWBPOkj62dvDOvSm6KbO9LR31sT1pnaWsYWc5XDzIdKHFm2FiROIHv0JOLMj40CpE4p6eKdA86dLXYWJpLNBTGsMIMT+Sot1x2DTssb7Z4ZXaC6ncDEy4dDdPri0In2tzBI90ZLt/QSCU/RTjVxE8en+QPr+oi8D00ff7fJRZIya8He/GACaeEC4yPVWhsevX0l6MjFRCwYIFNVNOpMy0c32dtop5K4DNSLmGKalGDih9QH4qQDGy66qu3Ff2ibx9Z8dpFFyJSIyQMGmMhbE0nYdoUfJf9xTydkQRZt4Su6XTZCR7dPsIvHx5kZ2+O9V1JrrixgdFSiXDEZNOD49Q3GiRSJmHToL7F5rLGVpaeBCkylYNUoK1Np8yJB13oJ3WQBbDjSZrWXYCVbD7isZd6MwyNTKfYlIekPfQqONlxvEKWxcY4Zy0yuGBJCF2DsAXP7c4hdBvpe/i+S3MqxEAGnu71GchWg+z9Wx229FVvLSmWA9pSJv6B+Ct9nGKO1U0+QxNZQokGrGiC91y6kMA/Mr3gfKUJwTVtnWxsaefipoU0GNYRQTYIXvmhdEGzTbnsMz5eIaqbFH2PgueyLTPJWKVCyfeY9FyGyyVSdohGK8Reb4qi7/KNTd0MjB5+d8KRfBGQitg02REGx0ossMOYmk6DFabJtmmLxFltN/DgkyM0iDBNjSE+c8sqrEbo3pZhx/YcMQzWrIyzZmUdDz84zqObxgkqAcmToIShopwK5v9UpAbVLVyMbF0EQkMGAaXMBO9IFfFKeaRbIfBdCEBKHxB4UsPUAFwk4JZyCA1Wtdg8uzvPyHiBPZOCK9ca+IGke8zi0R1pbn5LA3gO16ybTlogfcKWhqZxMNAC3bsnWN1uM14yMCNxZBDgOyV0OwxSQg3MaA8I6QZLYgmWxBIUPZe+Qp5nJodxqd4C09RkE00cPE1ULPgk6gxG3el9csCSAYYmiOomui5YE08y4VR48IUh6pIaH/nRE+AKrnxbI0FQXX4+UJThAB2I6gbZgkfGcEAITE2jv5Sn0Q5jCo2YYSI9uGtzL7dctZzNY8O8b+livvHjXVx4cYrlTXW8lMlQxGNgsExdwuKGdy4kZGt0b8liy33cvHKFyl+sKLOsdv5C1hCh6QfvotU0Yg3VGa4MAoSm8Vx3LwP7+jhvWYg9I0UWJg1y5YBYyMD1QafIzjHJ0vY6HMfHkB5RXVTTXmiCqCjwd3/Uzv3dBd66WkcKCAKfTBGKjk9bcwy9cjDd35olKex4kqeenKKzo3r4yfc9ZKmAXVd/5ABqRMQwWZVIsSoxvQfdBV4Q4Po+w+Ui+0s51pyRwNcg65aJmzbrk40kbJuC5xIEklLg0WSH6RnuZf3KJKvqUpzb3sLeQpaeQoYgkEc9tR3WNJZGE6xtayDtVMj7LlZM47m+CZYvS+JLyWNjQyRtiwvPbCDvuRT2w3OVMdaujvG5L3TzlksaueaaZjwgntQxYoKYY9IZj1JY6rPthRz5LpdkaHaqQym1YzpRxY8OubQE+Esp5T/NTY9mhxDiO8D1wKiUct30tQ3A7UCI6jnCP5NSPnUs7apAO48cSK94xrrFnLFuMeXcFA3lXZiaiyMEugZjJciOV1iywELHJ5kIM14y0XEYGslQ39LK2V0l/vWhES5aoiP9MJHWLh56oofGlMEj3XluaU6AYSKnk/Gb0QRWOMbvXyCwEw34Thk3O0WkaeFcvh1zwtA0DE2jy0zQ9Rq3w0SN6qw3Pn3seG1dA5auUWdabO4fxbU8TF1/1c2bvBuQ9xyynsPO/BQddpRASPbsLvKMOU5HS4SwozMwUsZLganpJDt1jOfAbYRwRKfiBOzpKbB0aZQlS2MEHtPZNQSJeotgGfzw0T3cevXq4/wuKXPp3gf6jygfd+2V7TNK5CCl3EE1SQVCCB0YAH42s57OzLf3vnTEOG/pWn3cy+QBXwH+Wkr5SyHEtdPfX3Ysjao1o3ksFE+y7MxzsJP1jGUDJnIe+0cqtDeYDE+54DnccHaCdR1hljfbtLbVMziS46Eem4aGOCN5gR6pQwqNhUnB0pTBTecm0KyDQRbAzad5qqdIyagGFqeQJdqyCCM0p+Uo55VF0ThIwfPpcaKuha3r1cwX8HJS/1fQoeB5DJWKBF7AY7tHuPeBIbrqYtxxTy9+IBnXyuRiZXqyGe7bsZ8nekawl8BAX4kPfrAT29KI+Sa6I3CKoBmQzVaQQpLQDJyCz6YdYyf2jVBm1XSQPaJ83PT142Uj0COl7D2ObR6T6SB7xDinr79pr1ImTwJ1018ngMFjbVfNaGtAvGURS8tFRsYLpKKCiCXQG6q3bAyMl0iEAQ1e2JUm2raSid07IZBMFSvkpybJjuUIWzqZUoXxrCRuVavaPLevzBmdIfxChg43zZbuRi7e0Mb+sQorUif3wbKTUUs4Qks4Ag2wKzeFpWvsKeZfsT9bynmEowbShb5Kgc5onNHRCjEsli4Ps3xRnMloicFMgXTGIzACpCUpedX9dQed05cnKZs+HfEIE16FdN6npd7CEQExzWKgUMDNQmtLmLOWNxJIiTYPk44oR/Va5eOORxpGmGE6wuPkRIzzgE8C9wsh/oHq5PTCY21AzWhrgGFHwK3QnDBY1lwNgCFdgO/RloBYLELUhJVNsCjYw9XL4Q/O0rh8bZRoQzMWPr9+bpLAC1jabPJET4Xn+1w6G01cr1pWLu8EpKIaW7bu46t376dQqp0Tx3NheTxJWzhOm/3KvxXjE2WG+kuUii51YZMXMxMUYhUc3WVLeYzvvrCLK5YtZPdUltyQx+r6BL4bMNBXwY7qTE1WmHAqlCYCdNdAZjXKBZ/NWybxZMADD40yPuSQShl4WkB/vsBzu4+ePU2Zl2a1fJwQwgLeDvzkeLQ3AyeyTN6twKeklB3Ap4BvH2sDKtDWCu3QxQmteq/t9DWvXKBQ8qqzFt8nbArSRVi8pBNT10jVmVx7dh1hWwPf44KlNvVhiS8lkwWPX3UXGc8ELF8YobW1iTNW1BMJnRxF3eezpfEEEcPEKB2cTdoRk2SDSSRlUgx8+iZKbN2eo1Dn0dIY4sKOBTiVgOV2PZevaiHjucTiFudvqKNSCDAMgWbAvrEC+azD/t4SD90/RlIL4TpwzTXNxCs2EcvEdyT79ud5Ymzw6MvXynz0amXijlf5uN8DnpVSjhyn9t6s2R7noT4A3DX99U+Ac4+1ARVoa4Q4kAJRM9myv4z0KgjLoOJBruizddDD9XyEgP1ph6EpifTK+E4Jz5O4VgpNSm5/aIpiJcA0BKmoxkjG47ozYqxbVg9I2tqbuenyrnmZ3/hkowtBnWmya3e1kk+x6JJK6YSjBz80pRIWq9dUt4fKgc+kKKGFJXp9QI+YYjxdIawZWJrJwliYRi/C5kfTdCwKUWn0WLe0jituaGCwlGP3jhweAVaHZLhSIhE2aGi0MJsEj3QPz8l7oBx3s10+7k0XPz/OTmSZvEHg0umvrwB2vcZzj0rt0dYAr1JCVooEgAxcTl8cRmg6mmEznisQS0Q5/8wke3v2s32owCUrw5g6OFNjZF0DFx0rpGFZMa49U0doPjuGy9gGdNQbbN5dYv1pLfzL/+vjMzcnSdgalewkVjylAu4MbUg1YV+i8Xh6lMAHPaRDAI7nY1mvPJVckgGuW+HZiVEcGVDxJFNZh7CvYTRZrG2ox4hq7Nib5xv/tIeJSYd33bQQ09DI5DwoaIjFsG1XluaWENGIyehQGcPWWbnEf/VO1hDp+xS276XUs5/K3n78sTQY0/UODYP6G64g0tmGHp2faSmvvbL9P+99oB+O86ljACFEBLgK+OhM25qpW7pW/+e3974Ex/nU8XSZvMuARiFEP/BXVMvxfW26+HsZ+JNjbvdUScFYyzynQnrns9Wvg+rfDXQTzbR4ZMsE6xZHqY8bbNtfojEqaYofXPb1MJmQMbpaYnjlEuNZn4goUqw4WIHHYzvLtCY0pkqSh3oEX//zKwAop8emE1YEmNG6o3VLeYOyToXfjvUz6FQL2AV+QCbrEAoZhMPVz8K+56MbOuWyRyg0/fk4gGYrzP2bhtiwPsmyBXW02TFKrktnOM6n73yaF57P0NoUxg0CpANOWXLzezooChfTFmzaNIEfSDrawnzh986p2eQVMggIKg79X/sB5F8/IxeWwcKP/xFm8sT+21YpGGuTmtHWAMM6mHDg5YJDvkvgu6xpt2lbvJDi2BCrO6LglnA8n3TexwtgQQIe7x6jOeIxOFHBSC1k+45xYqbEDhlccVqc53sr3PXMFKGwTcXxsS2d//Wj/Vy2LsHZzUWiTc2Ekyd/4feTVZ1lc01rF0+nR+nJpskD4bBxMKACBzKYHLhWLHhEwgYjXokNZyY5p6me3z45Rexck91TGf7tzh6e3Jymqz2MHRdsWJfioQdGGB/02LE7h2lCWfPxfMnlG5so5/yaDLLl/cOM3nEv0vXAccF7YzP3B4cFZ3/le6z92z9D02vvfVFOLBVoa0yh7GGaBtb0pLWlpZHADxiectGNgJgRIKXADTR0ARUXOhMO+XSWhTEbXYzRsjIGQuO5nRNkwhaJkE59Yz1/94Gl6NKhUDJZVh/Q1WjglfJkerNs37qLaNNiVqxon9s3YJ7SNI3zGlroitbxy6HeI1My6occPvOhUPDIZlxsU0cYgqeMCc69oImRqSKeA0sWxVh8ps2W57Ls3J1nz54ciYjNaLpA/9YSLUts8oWA9evrCHkat2yonaQVUkoqwxNMPfA4lZ79rwyuugDTrC7Jl10wDPDcg7U1NIFoqedt57fh9I9AIKtJPhRlBtTScY1I927Hy6Vf8zk5RyB9j7pw9S+HN53SsVDyGcrAyvYw+9Iay9ui7B/KIKRHVPfIFCXNy9eyoCmGZlgIIZCymjrw/ke3ES2PkoxqREyNrnMuUvu2x4GUkl/072HccwioTmglsCKSYKBUIDvp0ZyyyeDQ05Onodki6ptEogZvSS1kYLLIjvEpzl/RRE9fnj0jGTbtGuPOuwdZs6qOpiabH/7Pi2mssRJ5Yz/7NVue2cPysKwm7A5bkD24VKw1JpGOg5AC6iIYsShO71B1thuxIV8CQ6f1U+/HSsRxRiawmhtOWP/V0nFtUoG2RgS+x8T2Z19Z2ecwwgqTnirxm615Kr7G750Woq4uSrS+gXKpDMU0CJ1CEOKZ7WNcsNTE9TX6vCYuPrsTr1LBL+dJexFSZhkrGsct5qgEBrnB7ZhBtahqYulaQmrfdsaKrsu+Yo5F4Ri2YZD3HFKHVdyRUiKBsu8R1o3X/ZBz+907ed/VnUQjtZlwpPfL/wa5IsQiUC6DZUGxWvABy6jObuvrEGjsyQR0BEWspgRGPEL9VRcx+cDjEEjiZ67Bbm8m8+gzNL5j4wnrvwq0tUltPtQITTeIt3W9zrME4YjBuy5sZPXSRjb3CXonfSL1zYTCNpPUM2F3UiiWaIhpbBt0GK9YNDYkCFyX3MAuKvksA5MOZcdDAqFkE5ZwqUs1V09uEjC568UTMOLaFzFN1iTqiVkWpqYdEWQBhBBoQhAxzDe0kvCn71hRs0EWoPG6SzEWtZC4+Awi61ditTZRMm2Cpnqip68icc1F6HYIOZFGlMvsy7nYbQtw9g0w/P3/wh0aw+kdpLS7FyMRJ/W2yxn56f1zPSxlnlN7tDXECMdf83HdNLHdCvc8m+Oc0xayanGCiKUxNTaKLOWJmxpoee7dF9AeNcmXAza94PLfP5DCCBk0LN8AwDmHtRtOLSCcWkBkQQc7d+ymTmZnZ4CK8jqi65ZjLWyqlm9cvRSzIUn39x9hXWmM+usvJf/Ui4Q6mhGLWji9Pkmmu4fS77ZVX+wWOLC+tzvVxsN3Pc/iFS1sfmqCSxfuYP2FK+dsXMr8pma0NUS3bJLLN4BRnfkEQfByxp9SxccrZADBpStNEowRrkxAcZKeXfv49weGGRnLMDY6xnWn2WxJR/jYd0f50o/3ELLf2GkQwzRZs2417aefN0sjVJTXZ9YnqewfQVYcZCC59Op1NNx4BZU9/YRXdGEk4jjlMk/8aBP09h+lAYPTF8cJBkZYa5fxyw6xbJpc0T3yuacQIcSnhBBbhRDdQogfCiGOXGKZ54QQ3xFCjAohug+5droQ4gkhxItCiJ8LIY55X0zt0dYg33WY3LUFz3XQhEDTBJ4vMfTq0qIHuA6EwyZWJA66yXPbh1nbIgCBlWgg0bF8TsegKDM1/vOHKe7uRdd1vHwRPA+cgzm6dxZ0VkR9iIcxkwmEbVJ3ydmUd+wjtfF8JvpGSdZHKW7bwyNDkpG8x0c+OLsfIo/HHu3Q3b86onxc6zuunmkihzZgE7BGSlkSQvwYuFdK+R8zaXcmto0NHzHONU0tMx3nJUAe+N4h9WifBj4jpXxECPEhoEtK+RfH0q6a0dYg3bRILV2HoWtoWjW4GrrA96sHpSpln7AFhYJLNj2Jk51gZMqhgkZ80QoVZJWa0Pi2y9AkeKOT1QNRB4Js2OQ5ewEr6nVEIgaWRezsNURXL0WPhPAzeYQQNC1rY9gz2RJrw/cDzuh7kfSWHXM7qNcxHWSPKB83fX2mDCA8nSEpwpsoF3e8TAfZI8Y5ff1Ne5UyeSuBR6e//jXwrmNtVwXaGmXYYZrWXUBq+ZmISBI4eC+mMb0zX3QhnkyR7FrDH7zzEtrXnUeorn6Oeqwox1/7pz/w8l85e1UXsUvOovGGq7j6vReiGTrxizbQfut7KTy3Hc22mbz3MayFB5OvdCyIcsmGZt5+01nUV4pkf/XbORrJG/Za5ePeNCnlAPAPVGeOQ0BGSvmrmbQ5Q7MyzlfRTbViEcBNQMexNqAOQ9U4w7ZpXFJNRlCaHMWIxBGaQNMNmvTa+N/ffd2n4IwVrPv8rXPdFeUktPhvPg6AO5YGTTB5zyOkrr0UaZtk791E9t5NWKu6KO3uJSiWqPQNUlpQT3jZIoSh87PH9rPxjAWklyzhzFuum+PRvK5ZKR8nhEgBNwBdwBTwEyHEH0kpfzCTdmfgRJbJ+xDwdSHEXwL/BTjH2kBt/KVV3pBwfW2mSfQLDv7m7XPdDeUkZzalAGi86RrKewcQmvbyKWOnpw+vLo4WtqkMVCvAhToXIgyda85pJRIyOOvD189Rz49JH9Vl1KNdn4krgb1SyjEAIcRdVAugz1Wgna1xHkFKuR24GkAIsQI45k9baulYmfdOf/gbnPnrr811N5R5Qg+HiK5ZSvvHbqbpD69Ha0yC6xNMTOF7LtI0KG/fy+iP7wMgEjKQnkfgeq/d8MlhtsrH9QHnCyEionrD9kbgpRm2ORMnrEyeEGLB9H814HPA7cfahgq0iqKckoRhEFm9hI5Pvp8FH70JDB05PAmT1fvAKzt7+ekPnsDL5vGyBfypk//+8OnTxR8Beqlm7ewFPjLTU8dSys3AT4FngRepxo5vzay3b9706eIjxnkcTh3/EHgCWCmE6BdC3AK8VwixE9hO9QDYvx9zu+r2HkVRlGo6S6/iMPjFf3m5yIDjgxMOs/JTN2PEo7PeB5WCsTapGa2iKArVdJZmyGbx33ycltv+EABLh5hTYviuB/DfSB1bRTkKFWgVRVEOYzc3sOhvbwOzeh+6v6uX/I59c9spZd5SgVZRFOUohBAs/qvbaL6tmgNh6t5H5rhHynylbu9RFEV5DaHmRhZ//uNz3Q1lHlMzWkVRFEWZRSrQKoqiKMosUoFWURRFeV1CiE9Ml8jbKoT45Fz3Zz5Re7SKoig1pPdzXz+ifNziz398pokc1lFNEHEu1Vy/9wkh7pFS7pppf08FakarKIeo5PI88D++wPb71QlTZf6ZDrJHlI+bvj4Tq4EnpZRFKaUHPALcOMM2Txkq0CrKIex4DID+hzbx6Fe/Oce9UZRjNlvl47qBS4QQDUKICHAtb6Jc3KlKBVpFOcyVX/oskbZWMsOTbH1m51x3R1GOxayUj5NSvgR8mWrh8/uALcC8qLJwMlB7tIpyFBfe9iEm0wUe+Mq3ePaOu1m2opn29evoOP+sue6aoryWWSsfJ6X8NvBtACHEF4H+mbZ5qlAzWkV5FfWpKOtuvIY9WYPsrn623n0f5czJX8FFOaXNWvm4Q8rFLQLeCfxwpm2eKlSgVZTXsObcNVx/yWJ0HYIAfnvPE3heMNfdUpSjmj5dfET5uJmeOp52pxBiG/Bz4L9JKdPHoc1TgiqTpyhvwI57fsPW3zzJ9wYauP7mS/nwDavnuktKDVJl8mqT2qNVlDdg5XUbWX7NZax7eis/2DzGjx6M8O4rjrYVpiiK8kpq6VhR3iDN0FlxwXpGc74KsoqivGEq0CrKMbr9Ly4H4PHuMcoVf457oyjKyU4FWkV5g+7ZlOPux4ZwXJ/N28axTZ2nt4/PdbcURTnJqT1aRXmDvvDd3RhamMm0zou9fewfLfLXH1qPlJLiyCQ9T77EwP2b6LjqAta989K57q6iKCcJFWgV5Q3IFQI+9q7FXHtRgpF0gU3bXBAwPlXm55tK3Pmfv8Wqi7Hgh4+R/tHTLH3vZ8nYMW68/dOc9q7L8MoOoURsroehKMocULf3KMrreH5Hnk/+4yBlp4AUOdqbGpksDfDHV6+i+7m9MDhK4fEnad71PMOhFAUzxKqp/XSnOqlziyzPDaIBWkOKSFOCSNRm49c/Q8vZ6hYh5ZXU7T21Se3RKsqreGlvmes+2c9nvzlFyLSJRmxGM0XaWx1+//JFfP+B7UQWRNhnhtBvuo4/2XMXjddfxhJnkopusi7TS1duiIBq5oDhgs/Izn6GfreD/3vRR/jHyCU8/62fzfUwFUWZZWpGqyjTyhWfYsnhX++e4L9+OwFSIxqOIRCUnByDU30sa6snVyxy88bVbO0bobVJZ+3iJA89O055JE2mXGJV2GHynoc5J5QjrYUo9fRBsUjINqFUwg2qv3Pm9M99z2PfovXcNXM3cOWkoWa0tUnt0SqnrL6hErd+eZSy4+IHARWvQn0sxVQhQ9lxiUYqbDwvxPlr49z3VBrTbOftF7Txs00D7B8f5qsfW0/Yrv4KpfMOU88+yof/9t3c+r3dfOanf8kLPRmuHN3O0/98F8NTFSJJm5I0CbZ0Yx/Sjzve8ie0X3oWb//pF7HronPzZiiKMmvUjFY55Xz2G/08vwuCIMAPfISQtC8wufzsOJ2tFo9v7Wdhk8HIZJmHnx/hqnNa+c4vdrNjf5ZAwg0Xd3D6siTv1AdZffM1APzqH36MvXYld97yZezGBF0dSZK4/IROrpV9pPcMMOloROM2se3bERGbtGeQNiN0FcZe7ttFX7iVDX96I1bs8JKiyqlAzWhrkwq0yinnI5/vZ3AcvnRbPYubbXRdEA3P7LjCVM8AFcflqtvuJ9rfy7vow9JhzZ/dxH2/3k7bA/eQt6NEC3k0S0dWPGQ4RI+RwpQBUadIo5MDILZoAR/eeSdCiOMxXGUeUYG2NqlAqyizQAYBv7n3BTZ/826K+/pJ7tnJoy3rWe1P0TTWhxeASYAEYs0pyiOvLISSWt3F+576DrplHv0HKDVJBdrapPZoFWUWCE3jyus3cPrFq2hKhtj3u52U/v7nPLozx8apISqeBN8lb8fYVwyx7LDXp1/ay4v//gs2fPTGOem/oijHj7q9R1FmUVMyBEDnWSv49Ldu5d1/cAbp1sU807SKFxuWocuArtxQ9cn2wc+9AfDQx7/K0NNb56DXiqIcT2pGqygnSKwuwgc/8VbSN1/Aw9/8OcW+QbwXKuSSy7nqT6/nkb+/Az2fQQiNQt8IGBotZ6mkFooy36k9WuWk8v3797B+oUZraxLNCNE4PSOsVVJKdehJeZnao61NaulYOam0h0tMpgu80L2fT/7vTfh+MNddmlUqyCpK7VOBVjkp9Y6W+cSNnei6+ieqKMr8pvZolZPKZW9ZQ7FY5nxNIxy2X/8FiqIoJzkVaJWTihCCaDQ8191QFEU5btS6nKIoiqLMIhVoFUVRFGUWqUCrKIqiKLNIBVpFURRFmUUq0CqKoijKLFKBVlEURVFmkQq0iqIoijKLVKBVFEVRlFk0K0UFhBBjQO9xb1hRFKW2LZZSNs11J5Tja1YCraIoiqIoVWrpWFEURVFmkQq0iqIoijKLVKBVFEVRlFmkAq0yrwghOoUQ3Ydd+19CiM8IIf5DCFEUQsQPeexrQggphGg85NqN09dWHdZuSQjxvBBimxDidiGENv3YfUKIKSHEL07EGBVFqS0q0Cq1ZjdwA8B0oLwcGDjsOe8FNgHvOex6j5RyA7AeWAO8Y/r63wPvm53uKopS61SgVWrND4F3T399GfBbwDvwoBAiBlwE3MKRgRYAKaUHPA4sm/7+N0Bu1nqsKEpNU4FWqTW7gCYhRIrqzPWOwx5/B3CflHInMCmEOPPwBoQQEWAj8OIs91VRlFOACrTKfPNqN34fev0uqrPV84DHDnveocH3junvD1gqhHie6iz4HinlL2fcW0VRTnnGXHdAUY7RBJA67Fo9sPeQ7+8AngW+K6UMhBAACCEagCuAdUIICeiAFEL8+fTrDuzRKoqiHDdqRqvMK1LKPDAkhNgIIISoB95K9XDTgef0AZ8F/vmwl/8+8D0p5WIpZaeUsoNqgL74hHReUZRTkgq0ynz0fuBz08u8DwJ/LaXsOfQJUsp/Ofwa1WXinx127U7g5tf6YUKIx4CfABuFEP1CiGtm0nlFUU4tKtexoiiKoswiNaNVFEVRlFmkAq2iKIqizCIVaBVFURRlFqlAqyiKoiizSAVaRVEURZlFKtAqiqIoyixSgVZRFEVRZtH/B/unuJOal+L8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 597 ms, sys: 12.2 ms, total: 609 ms\n",
      "Wall time: 607 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "sc.pl.umap(adata, color=[\"hdbscan_gpu\"])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (cuml_2108_070921)",
   "language": "python",
   "name": "cuml_2108_070921"
  },
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}