diff --git a/molpipeline/estimators/nearest_neighbor.py b/molpipeline/estimators/nearest_neighbor.py index 456d90ef..82540fd3 100644 --- a/molpipeline/estimators/nearest_neighbor.py +++ b/molpipeline/estimators/nearest_neighbor.py @@ -2,19 +2,22 @@ from __future__ import annotations +import multiprocessing from typing import Any, Callable, Literal, Sequence, Union +from joblib import Parallel, delayed + try: from typing import Self except ImportError: from typing_extensions import Self - import numpy as np import numpy.typing as npt from scipy import sparse from sklearn.neighbors import NearestNeighbors +from molpipeline.utils.kernel import tanimoto_similarity_sparse from molpipeline.utils.value_checks import get_length __all__ = ["NamedNearestNeighbors"] @@ -202,3 +205,206 @@ def fit_predict( """ self.fit(X, y) return self.predict(X, return_distance=return_distance, n_neighbors=n_neighbors) + + +class NearestNeighborsRetrieverTanimoto: # pylint: disable=too-few-public-methods + """k-nearest neighbors between data sets using Tanimoto similarity. + + This class uses the Tanimoto similarity to find the k-nearest neighbors of a query set in a target set. + The full similarity matrix is computed and reduced to the k-nearest neighbors. A dot-product based + algorithm is used, which is faster than using the RDKit native Tanimoto function. + + For handling larger datasets, the computation can be batched to reduce memory usage. In addition, + the batches can be processed in parallel using joblib. + """ + + def __init__( + self, + target_fingerprints: sparse.csr_matrix, + k: int | None = None, + batch_size: int = 1000, + n_jobs: int = 1, + ): + """Initialize NearestNeighborsRetrieverTanimoto. + + Parameters + ---------- + target_fingerprints: sparse.csr_matrix + Fingerprints of target molecules. Must be a binary sparse matrix. + """ + self.target_fingerprints = target_fingerprints + if k is None: + self.k = self.target_fingerprints.shape[0] + else: + self.k = k + self.batch_size = batch_size + if n_jobs == -1: + self.n_jobs = multiprocessing.cpu_count() + else: + self.n_jobs = n_jobs + if self.k == 1: + self.knn_reduce_function = self._reduce_k_equals_1 + elif self.k < self.target_fingerprints.shape[0]: + self.knn_reduce_function = self._reduce_k_greater_1_less_n + else: + self.knn_reduce_function = self._reduct_to_indices_k_equals_n + + @staticmethod + def _reduce_k_equals_1( + similarity_matrix: npt.NDArray[np.float64], + ) -> tuple[npt.NDArray[np.int64], npt.NDArray[np.float64]]: + """Reduce similarity matrix to k=1 nearest neighbors. + + Uses argmax to find the index of the nearest neighbor in the target fingerprints. + This function has therefore O(n) time complexity. + + Parameters + ---------- + similarity_matrix: npt.NDArray[np.float64] + Similarity matrix of Tanimoto scores between query and target fingerprints. + + Returns + ------- + npt.NDArray[np.int64] + Indices of the query's nearest neighbors in the target fingerprints. + """ + topk_indices = np.argmax(similarity_matrix, axis=1) + topk_similarities = np.take_along_axis( + similarity_matrix, topk_indices.reshape(-1, 1), axis=1 + ).squeeze() + return topk_indices, topk_similarities + + def _reduce_k_greater_1_less_n( + self, + similarity_matrix: npt.NDArray[np.float64], + ) -> tuple[npt.NDArray[np.int64], npt.NDArray[np.float64]]: + """Reduce similarity matrix to k>1 and k tuple[npt.NDArray[np.int64], npt.NDArray[np.float64]]: + """Reduce similarity matrix to k=n nearest neighbors. + + Parameters + ---------- + similarity_matrix: npt.NDArray[np.float64] + Similarity matrix of Tanimoto scores between query and target fingerprints. + + Returns + ------- + tuple[npt.NDArray[np.int64], npt.NDArray[np.float64]] + Indices of the query's k-nearest neighbors in the target fingerprints and + the corresponding similarities. + """ + indices = np.fliplr(similarity_matrix.argsort(axis=1, kind="stable")) + similarities = np.take_along_axis(similarity_matrix, indices, axis=1) + return indices, similarities + + def _process_batch( + self, query_batch: sparse.csr_matrix + ) -> tuple[npt.NDArray[np.int64], npt.NDArray[np.float64]]: + """Process a batch of query fingerprints. + + Parameters + ---------- + query_batch: sparse.csr_matrix + Batch of query fingerprints. + + Returns + ------- + tuple + Indices of the k-nearest neighbors in the target fingerprints and the corresponding similarities. + """ + # compute full similarity matrix for the query batch + similarity_mat_chunk = tanimoto_similarity_sparse( + query_batch, self.target_fingerprints + ) + + # reduce the similarity matrix to the k nearest neighbors + return self.knn_reduce_function(similarity_mat_chunk) + + def predict( + self, query_fingerprints: sparse.csr_matrix + ) -> tuple[npt.NDArray[np.int64], npt.NDArray[np.float64]]: + """Predict the k-nearest neighbors of the query fingerprints. + + Parameters + ---------- + query_fingerprints: sparse.csr_matrix + Query fingerprints. + + Returns + ------- + tuple[npt.NDArray[np.int64], npt.NDArray[np.float64]] + Indices of the k-nearest neighbors in the target fingerprints and the corresponding similarities. + """ + if query_fingerprints.shape[1] != self.target_fingerprints.shape[1]: + raise ValueError( + "The number of features in the query fingerprints does not match the number of features in the target fingerprints." + ) + if self.n_jobs > 1: + # parallel execution + with Parallel(n_jobs=self.n_jobs) as parallel: + # the parallelization is not optimal: the self.target_fingerprints (and query_fingerprints) are copied to each child process worker + # -> joblib does some behind the scenes mmapping but copying the full matrices is probably a memory bottleneck. + # If Python removes the GIL this here would be a good use case for threading with zero copies. + res = parallel( + delayed(self._process_batch)( + query_fingerprints[i : i + self.batch_size, :] + ) + for i in range(0, query_fingerprints.shape[0], self.batch_size) + ) + result_indices_tmp, result_similarities_tmp = zip(*res) + result_indices = np.concatenate(result_indices_tmp) + result_similarities = np.concatenate(result_similarities_tmp) + else: + # single process execution + result_shape = ( + (query_fingerprints.shape[0], self.k) + if self.k > 1 + else (query_fingerprints.shape[0],) + ) + result_indices = np.full(result_shape, -1, dtype=np.int64) + result_similarities = np.full(result_shape, np.nan, dtype=np.float64) + for i in range(0, query_fingerprints.shape[0], self.batch_size): + query_batch = query_fingerprints[i : i + self.batch_size, :] + ( + result_indices[i : i + self.batch_size], + result_similarities[i : i + self.batch_size], + ) = self._process_batch(query_batch) + + return result_indices, result_similarities diff --git a/molpipeline/utils/kernel.py b/molpipeline/utils/kernel.py index e949ebb9..b25ca221 100644 --- a/molpipeline/utils/kernel.py +++ b/molpipeline/utils/kernel.py @@ -25,8 +25,12 @@ def tanimoto_similarity_sparse( Matrix of similarity values between instances of A (rows/first dim) , and instances of B (columns/second dim). """ intersection = matrix_a.dot(matrix_b.transpose()).toarray() - norm_1 = np.array(matrix_a.multiply(matrix_a).sum(axis=1)) - norm_2 = np.array(matrix_b.multiply(matrix_b).sum(axis=1)) + norm_1 = np.array(matrix_a.sum(axis=1)) + if matrix_a is matrix_b: + # avoid calculating the same norm twice + norm_2 = norm_1 + else: + norm_2 = np.array(matrix_b.sum(axis=1)) union = norm_1 + norm_2.T - intersection # avoid division by zero https://stackoverflow.com/a/37977222 return np.divide( diff --git a/notebooks/advanced_04_dataset_similarity.ipynb b/notebooks/advanced_04_dataset_similarity.ipynb new file mode 100644 index 00000000..7ac91930 --- /dev/null +++ b/notebooks/advanced_04_dataset_similarity.ipynb @@ -0,0 +1,750 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e6afe352-c6de-4469-9508-a24cbd914d57", + "metadata": {}, + "source": [ + "# Analyzing similarity of molecular dataset\n", + "\n", + "This notebook illustrates how the `NearestNeighborsRetrieverTanimoto` can be used for analyzing the Tanimoto similarities of two datasets. \n", + "\n", + "Such analysis can be useful for many applications. For example, for analyzing how similar new molecules are to the training set to assess the applicability domain when making predictions. Alternatively the similarity of training and test set can be evaluated to understand how well the model generalizes. \n", + "\n", + "The notebook has the following sections:\n", + "\n", + "**How to compute dataset similarities?**\n", + "\n", + "**How to analyze similarities between train and test set?**\n", + "\n", + "**Comparison to native RDKit Tanimoto computation**" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "4b74dcad-0865-41a5-b380-fe6fdea89506", + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "\n", + "import matplotlib.pyplot as plt\n", + "from rdkit import DataStructs\n", + "from sklearn.model_selection import train_test_split\n", + "import seaborn as sns\n", + "\n", + "from molpipeline import Pipeline, ErrorFilter\n", + "from molpipeline.any2mol import AutoToMol\n", + "from molpipeline.mol2any import MolToMorganFP\n", + "\n", + "from molpipeline.utils.kernel import tanimoto_similarity_sparse\n", + "from molpipeline.estimators.nearest_neighbor import NearestNeighborsRetrieverTanimoto" + ] + }, + { + "cell_type": "markdown", + "id": "17a66d14-8a18-4363-9016-199feeefbff6", + "metadata": {}, + "source": [ + "For this notebook we use 20k molecules from ChEMBL35 as a dataset." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "dd318e60-226e-4fcc-a245-cbada6ac36bb", + "metadata": {}, + "outputs": [], + "source": [ + "df = pd.read_csv(\"example_data/chembl_35_20k.smi.gz\", index_col=\"index\")" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "d1013e21-f57f-4917-a764-ebf24c586b51", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
smileschembl_id
index
0Cc1cc(-c2csc(N=C(N)N)n2)cn1CCHEMBL153534
1CC[C@H](C)[C@H](NC(=O)[C@H](CC(C)C)NC(=O)[C@@H...CHEMBL440060
2CCCC[C@@H]1NC(=O)[C@@H](NC(=O)[C@H](CC(C)C)NC(...CHEMBL440245
3CC(C)C[C@@H]1NC(=O)CNC(=O)[C@H](c2ccc(O)cc2)NC...CHEMBL440249
4Brc1cccc(Nc2ncnc3ccncc23)c1NCCN1CCOCC1CHEMBL405398
.........
19995NS(=O)(=O)c1ccc(NC(=O)c2ccccc2)cc1CHEMBL23559
19996Cn1cncc1C(O)(C#Cc1ccc(C#N)cc1-c1cc(Cl)cc(Cl)c1...CHEMBL23578
19997CC(C)(C)C(=O)Nc1nnc(S(N)(=O)=O)s1CHEMBL23579
19998COC(=O)NCC(=O)N[C@@H](CC(C)C)C(=O)NC(Cc1ccccc1...CHEMBL23580
19999CC(C)c1nc2sc3c(c2c(-c2ccc(F)cc2)c1/C=C/[C@@H](...CHEMBL23581
\n", + "

20000 rows × 2 columns

\n", + "
" + ], + "text/plain": [ + " smiles chembl_id\n", + "index \n", + "0 Cc1cc(-c2csc(N=C(N)N)n2)cn1C CHEMBL153534\n", + "1 CC[C@H](C)[C@H](NC(=O)[C@H](CC(C)C)NC(=O)[C@@H... CHEMBL440060\n", + "2 CCCC[C@@H]1NC(=O)[C@@H](NC(=O)[C@H](CC(C)C)NC(... CHEMBL440245\n", + "3 CC(C)C[C@@H]1NC(=O)CNC(=O)[C@H](c2ccc(O)cc2)NC... CHEMBL440249\n", + "4 Brc1cccc(Nc2ncnc3ccncc23)c1NCCN1CCOCC1 CHEMBL405398\n", + "... ... ...\n", + "19995 NS(=O)(=O)c1ccc(NC(=O)c2ccccc2)cc1 CHEMBL23559\n", + "19996 Cn1cncc1C(O)(C#Cc1ccc(C#N)cc1-c1cc(Cl)cc(Cl)c1... CHEMBL23578\n", + "19997 CC(C)(C)C(=O)Nc1nnc(S(N)(=O)=O)s1 CHEMBL23579\n", + "19998 COC(=O)NCC(=O)N[C@@H](CC(C)C)C(=O)NC(Cc1ccccc1... CHEMBL23580\n", + "19999 CC(C)c1nc2sc3c(c2c(-c2ccc(F)cc2)c1/C=C/[C@@H](... CHEMBL23581\n", + "\n", + "[20000 rows x 2 columns]" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "df" + ] + }, + { + "cell_type": "markdown", + "id": "824a47b9-2ecc-4505-bf73-df06a70795e0", + "metadata": {}, + "source": [ + "## How to compute dataset similarities? " + ] + }, + { + "cell_type": "markdown", + "id": "8125630a-a8d0-4247-8929-bde009b6f550", + "metadata": {}, + "source": [ + "To start the comparison we need the fingerprints as sparse matrices." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "491e9e8b-9210-4d03-8b9b-42d8831389d8", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 1.82 s, sys: 657 ms, total: 2.47 s\n", + "Wall time: 16.9 s\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%%time\n", + "error_filter = ErrorFilter()\n", + "\n", + "fingerprint_pipeline = Pipeline(\n", + " [\n", + " (\"auto2mol\", AutoToMol()),\n", + " (\"error_filter\", error_filter),\n", + " (\"morgan2_2048\", MolToMorganFP(n_bits=2048, radius=2, return_as=\"sparse\")),\n", + " ],\n", + " n_jobs=-1,\n", + ")\n", + "\n", + "fp_matrix = fingerprint_pipeline.transform(df[\"smiles\"])\n", + "fp_matrix" + ] + }, + { + "cell_type": "markdown", + "id": "ae7865fe-67d7-441f-bd3d-1c1d7f6a4e84", + "metadata": {}, + "source": [ + "The resulting fingerprint matrix has the shape (19999, 2048) showing that 1 molecule could not be processed." + ] + }, + { + "cell_type": "markdown", + "id": "6425965c-9263-41ca-87e2-729693bf4a2e", + "metadata": {}, + "source": [ + "To make a data set comparison we need to define the target and the query data set. The `NearestNeighborsRetrieverTanimoto` will retrieve the k most similar molecules in the target data sets for every query fingerprint. In this example we use the same matrix as target and query data set and compute their 3-nearest neighbors using `k=3`. " + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "18455107-0f22-479e-9359-8a81ce66934f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 116 ms, sys: 14.2 ms, total: 130 ms\n", + "Wall time: 6.46 s\n" + ] + } + ], + "source": [ + "%%time\n", + "target_fps = fp_matrix\n", + "query_fps = target_fps\n", + "\n", + "retriever = NearestNeighborsRetrieverTanimoto(target_fps, k=3, n_jobs=-1)\n", + "indices, similarities = retriever.predict(query_fps)" + ] + }, + { + "cell_type": "markdown", + "id": "ef7b6df2-8741-44e1-aaba-94f4534db708", + "metadata": {}, + "source": [ + "The output of the retriever are a list of `indices` corresponding to the hits in the target dataset and a list of the hits' Tanimoto similarities" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "f87925c5-ec5c-4a7b-8ca3-9c550cd6cd2a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 0, 13038, 13544],\n", + " [ 2744, 1, 111],\n", + " [ 2, 2984, 24],\n", + " ...,\n", + " [19996, 3854, 10457],\n", + " [19997, 11806, 1485],\n", + " [19998, 19881, 19690]])" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "indices" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "a87fcc03-8941-4198-bc11-655b2c540ee0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(19999, 3)" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "indices.shape" + ] + }, + { + "cell_type": "markdown", + "id": "5fe38ce9-f778-4c68-83ee-a48e8535490b", + "metadata": {}, + "source": [ + "The `indices` array contains one row for each query fingerprint and three columns for the 3-nearest neighbors. The hits of each query are sorted from left to right in descending order. The `similarities` array has the same shape as `indices` but contains the Tanimoto scores. " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "b3977853-ccc3-4647-b9c5-90acb7c5fd95", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1. , 0.60784314, 0.30909091],\n", + " [1. , 1. , 1. ],\n", + " [1. , 0.97986577, 0.87179487],\n", + " ...,\n", + " [1. , 0.68571429, 0.575 ],\n", + " [1. , 0.41860465, 0.41836735],\n", + " [1. , 0.96969697, 0.63291139]])" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "similarities" + ] + }, + { + "cell_type": "markdown", + "id": "3cd1973b-4563-4dd9-8a66-1f93ef0ba4c0", + "metadata": {}, + "source": [ + "Since we used the same dataset for the query and the target dataset, we always find a molecule with a similarity of 1.0 because the query itself is contained in the target dataset. However, sometimes there are multiple hits with the same Tanimoto score of 1.0." + ] + }, + { + "cell_type": "markdown", + "id": "d8f41e49-a307-476a-a168-4a01e3276581", + "metadata": {}, + "source": [ + "## How to analyze similarities between train and test set?\n", + "\n", + "The nearest neighbors can be used for analyzing the similarity between training and test set which can be an essential tool to better understand the generalization capabilities of machine learning models. In addition, this information can be used to select an appropriate data splitting strategy." + ] + }, + { + "cell_type": "markdown", + "id": "2cd5fc30-f176-4a1e-8235-3d77bfc0ed69", + "metadata": {}, + "source": [ + "First we make a train/test split with our ChEMBL data and a dummy y vector because we don't use the labels in this example. " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "48c61bce-af90-48b9-9222-b0d984be6fa5", + "metadata": {}, + "outputs": [], + "source": [ + "# let's use dummy values for y\n", + "y = np.zeros(fp_matrix.shape[0], dtype=np.int64)\n", + "\n", + "X_train, X_test, y_train, y_test = train_test_split(\n", + " fp_matrix, y, test_size=0.33, random_state=42\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "a6706966-fce9-4aca-b175-ce65827fdbb9", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "X_train" + ] + }, + { + "cell_type": "markdown", + "id": "e18d7248-18c5-46bd-bfcf-d15856a0871b", + "metadata": {}, + "source": [ + "We use the `NearestNeighborsRetrieverTanimoto` to get the 1-nearest neighbors of the test compounds in the training set" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "3d3a9ab0-891a-4990-a5f6-2ce616fd2bc6", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.7826087 , 0.74736842, 0.85074627, ..., 0.7752809 , 0.60869565,\n", + " 0.81538462])" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "retriever = NearestNeighborsRetrieverTanimoto(X_train, k=1, n_jobs=-1)\n", + "indices, similarities = retriever.predict(X_test)\n", + "similarities" + ] + }, + { + "cell_type": "markdown", + "id": "4f123918-6595-4d95-8396-fb3a8eac3185", + "metadata": {}, + "source": [ + "Let's look at the mean similarities of the most similar compounds in the training set" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "11b2789c-9940-4958-889c-ecf107d26516", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.6575015090173907" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.mean(similarities)" + ] + }, + { + "cell_type": "markdown", + "id": "ab4e809c-ff50-4871-b4b0-1d290f71c9f7", + "metadata": {}, + "source": [ + "We can also plot the distribution of similarities to get a better impression how similar the train and test set are to each other." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "401ee1a7-04c7-4ac6-be9f-0eada88a1021", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0.5, 1.0, '1-nearest neighbor Tanimoto similarities to training data')" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGzCAYAAADJ3dZzAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAATTJJREFUeJzt3XlcVNX/P/AXIAw7iAIDyiK4oiiGG5o7iXsq5pIa7qaoKW65C5qUWfrJUDNLrDT7aC5p5L5lYqlJ4pIfUZRUQE0H3Fjn/P7wx/06MuwDM1xfz8djHjX3nrnzvmfujC/uPfdeIyGEABEREZFMGeu7ACIiIqLyxLBDREREssawQ0RERLLGsENERESyxrBDREREssawQ0RERLLGsENERESyxrBDREREssawQ0RERLLGsEOvBCMjIyxatKjUr504cWKR7RYtWgQjIyPcv3+/VO9TGeWt86vG09MTw4cP1+kyX95Go6OjYWRkhBs3buj0fTp06IAOHToUq+3w4cPh6emp0/eXo7J8VkePHoWRkRGOHj2q87pKS47fa4adMnr8+DEWLlyIrl27wsHBAUZGRoiOjtZ3WQYlJiam1EGDiifvB7M4D7ngdqUbd+7cwaJFixAXF6fvUgAAq1evLpff0KVLl2Lnzp06Xy5pKq/Pr6yq6LuAyu7+/fuIiIiAu7s7mjRpYlDp3FDExMQgKipKr/8wPXv2DFWqyHdzb9CgAb799luNabNnz4a1tTXmzp1bbu87b948vP/+++W2/MLoc7u6cuUKjI11+7diRW2j+/fv13h+584dhIeHw9PTE35+fhrzvvzyS6jV6nKv6UWrV69G9erVdb7nbOnSpejfvz/69Omj0+UCwLBhwzBo0CAoFIoSv7Zdu3Z49uwZzMzMdF6XPpTX51dW8v31ryAuLi5ITk6GUqnEmTNn0Lx5c32XVCY5OTlQq9Wy+eLlMTc313cJOiGEQEZGBiwsLDSmOzs7Y+jQoRrTPvzwQ1SvXj3fdF2qUqWKrENkQUrzj1pRynsbffr0KSwtLUv03TY1NS3HigzXkydPYGVlVez2JiYmMDExKdV7GRsby+b3yZDxMFYZKRQKKJXKMi0jb0zIzp070ahRIygUCjRs2BB79+7N1/b27dsYOXIknJ2dpXZff/21RpusrCwsWLAA/v7+sLOzg5WVFdq2bYsjR45otLtx4waMjIywfPlyrFy5Et7e3lAoFLh06RIA4O+//0b//v3h4OAAc3NzNGvWDD/99JPGMrKzsxEeHo46derA3Nwc1apVw+uvv44DBw4AeH7MPyoqSlrP4hxK8fT0RM+ePXHixAm0aNEC5ubm8PLywjfffJOvrUqlwpQpU+Dm5gaFQoHatWvjo48+yvfXqLYxO0ePHkWzZs1gbm4Ob29vfPHFF4Ueqy7O5wM839s3YMAA2Nraolq1anjvvfeQkZGh0SYnJweLFy+W+tzT0xNz5sxBZmam1r7Yt28fmjVrBgsLC3zxxReF9l9BSrNdrFu3TqqxefPmOH36tEZbbf2Vtz1v3boVPj4+sLCwQEBAAOLj4wEAX3zxBWrXrg1zc3N06NBB6ziHrVu3wt/fHxYWFlJgu337tjS/qO3qyZMnmDZtmrRd1KtXD8uXL4cQosh+unr1KoKDg6FUKmFubo6aNWti0KBBSEtLk9q8PGYnb8zGiRMnMHnyZDg6OsLe3h7jxo1DVlYWVCoV3nnnHVStWhVVq1bFzJkz89VSnHFlu3btQo8ePeDq6gqFQgFvb28sXrwYubm5Gu06dOiARo0a4ezZs2jXrh0sLS0xZ84caV7emJ2jR49Kf6CNGDFC6se8wxDaxuyo1WqsXLkSDRs2hLm5OZydnTFu3Dg8fPhQo92ZM2cQFBSE6tWrw8LCArVq1cLIkSMLXT9PT09cvHgRx44dk2p5cXzR9evX8dZbb8HBwQGWlpZo1aoVfv7550KXCTzv2ydPnmDjxo3ScvM+v7xt+NKlS3j77bdRtWpVvP766wCA8+fPY/jw4fDy8oK5uTmUSiVGjhyJf//9V2P52sbsFPd3TNuYnbzP79KlS+jYsSMsLS1Ro0YNLFu2LN+63bx5E71794aVlRWcnJwwdepU7Nu3r9jjgE6cOIHmzZtr/A5qs2HDBnTq1AlOTk5QKBTw8fHBmjVrNNoU9vk9ePAA06dPh6+vL6ytrWFra4tu3brhr7/+KrJGXXj1/iQzUCdOnMD27dsxYcIE2NjY4LPPPkNwcDCSkpJQrVo1AEBqaipatWol/WPi6OiIX375BaNGjUJ6ejqmTJkCAEhPT8f69esxePBgjBkzBo8ePcJXX32FoKAg/PHHH/l2VW/YsAEZGRkYO3YsFAoFHBwccPHiRbRp0wY1atTA+++/DysrK/z3v/9Fnz598OOPP6Jv374Anv9QREZGYvTo0WjRogXS09Nx5swZ/Pnnn3jjjTcwbtw43LlzBwcOHMh3mKUwCQkJ6N+/P0aNGoWQkBB8/fXXGD58OPz9/dGwYUMAz/9Sbd++PW7fvo1x48bB3d0dJ0+exOzZs5GcnIyVK1cWuPxz586ha9eucHFxQXh4OHJzcxEREQFHR8dSfz55BgwYAE9PT0RGRuLUqVP47LPP8PDhQ40fudGjR2Pjxo3o378/pk2bht9//x2RkZG4fPkyduzYobG8K1euYPDgwRg3bhzGjBmDevXqFbsfX1TS7WLz5s149OgRxo0bByMjIyxbtgz9+vXD9evXi/yL/9dff8VPP/2E0NBQAEBkZCR69uyJmTNnYvXq1ZgwYQIePnyIZcuWYeTIkTh8+LD02ujoaIwYMQLNmzdHZGQkUlNT8Z///Ae//fYbzp07J4WIgrYrIQR69+6NI0eOYNSoUfDz88O+ffswY8YM3L59GytWrCiw7qysLAQFBSEzMxOTJk2CUqnE7du3sWfPHqhUKtjZ2RW63nmvCQ8Px6lTp7Bu3TrY29vj5MmTcHd3x9KlSxETE4OPP/4YjRo1wjvvvFPo8l4WHR0Na2trhIWFwdraGocPH8aCBQuQnp6Ojz/+WKPtv//+i27dumHQoEEYOnQonJ2d8y2vQYMGiIiIwIIFCzB27Fi0bdsWANC6desCaxg3bpz0GU2ePBmJiYn4/PPPce7cOfz2228wNTXF3bt30aVLFzg6OuL999+Hvb09bty4ge3btxe6fitXrsSkSZM0Dr/m1Z2amorWrVvj6dOnmDx5MqpVq4aNGzeid+/e2LZtm/SbpM23334r/UaNHTsWAODt7a3R5q233kKdOnWwdOlSKYgeOHAA169fx4gRI6BUKnHx4kWsW7cOFy9exKlTp4r8w604v2MFefjwIbp27Yp+/fphwIAB2LZtG2bNmgVfX19069YNwPNQ36lTJyQnJ+O9996DUqnE5s2b8/0BU5D4+Hjpc1q0aBFycnKwcOFCrdvKmjVr0LBhQ/Tu3RtVqlTB7t27MWHCBKjVaul7Xtjnd/36dezcuRNvvfUWatWqhdTUVHzxxRdo3749Ll26BFdX12LVXGqCdOb06dMCgNiwYUOJXgdAmJmZiYSEBGnaX3/9JQCIVatWSdNGjRolXFxcxP379zVeP2jQIGFnZyeePn0qhBAiJydHZGZmarR5+PChcHZ2FiNHjpSmJSYmCgDC1tZW3L17V6N9586dha+vr8jIyJCmqdVq0bp1a1GnTh1pWpMmTUSPHj0KXb/Q0FBRkk3Nw8NDABDHjx+Xpt29e1coFAoxbdo0adrixYuFlZWV+N///qfx+vfff1+YmJiIpKQkaRoAsXDhQul5r169hKWlpbh9+7Y07erVq6JKlSr5ai3u57Nw4UIBQPTu3Vvj9RMmTBAAxF9//SWEECIuLk4AEKNHj9ZoN336dAFAHD58OF9f7N27t+AOK0DDhg1F+/btpecl3S6qVasmHjx4IE3ftWuXACB2796db51fBEAoFAqRmJgoTfviiy8EAKFUKkV6ero0ffbs2QKA1DYrK0s4OTmJRo0aiWfPnknt9uzZIwCIBQsWSNMK2q527twpAIglS5ZoTO/fv78wMjLS+Bxfdu7cOQFAbN26tcA2Qjz/XEJCQqTnGzZsEABEUFCQUKvV0vSAgABhZGQk3n33XWlaTk6OqFmzpsZnI0T+bTRvmS/2Y953/EXjxo0TlpaWGt/V9u3bCwBi7dq1+dq3b99e470L+90KCQkRHh4e0vNff/1VABCbNm3SaLd3716N6Tt27BAAxOnTp/Mtsygvb7d5pkyZIgCIX3/9VZr26NEjUatWLeHp6Slyc3MLXa6VlZXGZ5YnbxsePHhwvnna+vv777/P9/uk7bMq7u/YkSNHBABx5MgRaVre5/fNN99I0zIzM4VSqRTBwcHStE8++UQAEDt37pSmPXv2TNSvXz/fMrXp06ePMDc3Fzdv3pSmXbp0SZiYmOT7bmnri6CgIOHl5aUxraDPLyMjI99nlJiYKBQKhYiIiCi0Tl3gYSwDERgYqPGXRuPGjWFra4vr168DeP7X6o8//ohevXpBCIH79+9Lj6CgIKSlpeHPP/8E8Pz4cd5xebVajQcPHiAnJwfNmjWT2rwoODhYY4/GgwcPcPjwYQwYMACPHj2S3ufff/9FUFAQrl69Kh1SsLe3x8WLF3H16lWd9oePj4/0VyYAODo6ol69elJ/AM8PdbRt2xZVq1bV6I/AwEDk5ubi+PHjWpedm5uLgwcPok+fPhp/TdSuXVv6i+llRX0+L8r7KyfPpEmTADwfUPvif8PCwjTaTZs2DQDy7ZavVasWgoKCtNZVEiXdLgYOHIiqVatKz/M+D23r/LLOnTtrHP5o2bIlgOfbmo2NTb7pecs8c+YM7t69iwkTJmiMY+jRowfq169frEMWMTExMDExweTJkzWmT5s2DUII/PLLLwW+Nm/Pzb59+/D06dMi3+tlo0aN0vhrv2XLlhBCYNSoUdI0ExMTNGvWrFj9+LIXx2rlfTfbtm2Lp0+f4u+//9Zoq1AoMGLEiBK/R2G2bt0KOzs7vPHGGxrfOX9/f1hbW0t7FOzt7QEAe/bsQXZ2tk7eOyYmBi1atJAOMQGAtbU1xo4dixs3bkiH30vr3XffzTftxf7OyMjA/fv30apVKwDQ+p15WXF+xwpibW2tMd7OzMwMLVq00Hjt3r17UaNGDfTu3VuaZm5ujjFjxhS5/NzcXOzbtw99+vSBu7u7NL1BgwZaf29e7Iu0tDTcv38f7du3x/Xr1zUO8RZEoVBIg/pzc3Px77//wtraGvXq1StWX5YVw04FSUtLQ0pKivR48OCBxvwXN7Y8VatWlY6D37t3DyqVCuvWrYOjo6PGI+8H7e7du9JrN27ciMaNG0vjaBwdHfHzzz9r3Shr1aql8TwhIQFCCMyfPz/fey1cuFDjvSIiIqBSqVC3bl34+vpixowZOH/+fBl6qnj9ATwfW7F37958NQYGBubrjxfdvXsXz549Q+3atfPN0zatuPXkqVOnjsZzb29vGBsbS8fzb968CWNj43zvpVQqYW9vj5s3b2pMf/nzKYuSbBcvr3Ne8NG2zkW9Ni9EuLm5aZ2et8y8ddd2qK5+/fr5+kabmzdvwtXVVSNUAc9/xF98D21q1aqFsLAwrF+/HtWrV0dQUBCioqKK9WMOlGy9i9OPL7t48SL69u0LOzs72NrawtHRUfoH8eUaa9SoofMTDa5evYq0tDQ4OTnl+949fvxY+s61b98ewcHBCA8PR/Xq1fHmm29iw4YN+caklcTNmze1bhfF+VyLQ9v37MGDB3jvvffg7OwMCwsLODo6Su2Ks02U5HfjZTVr1sx3mOzl1968eRPe3t752hX0O/aie/fu4dmzZ/l+rwDt37/ffvsNgYGBsLKygr29PRwdHaVxYMXpC7VajRUrVqBOnTpQKBSoXr06HB0dcf78+WJ/v8qCY3YqyHvvvYeNGzdKz9u3b68xeKygkfzi/x87zhtwO3ToUISEhGht27hxYwDAd999h+HDh6NPnz6YMWMGnJycYGJigsjISFy7di3f614+syfvvaZPn17gHoW8L1O7du1w7do17Nq1C/v378f69euxYsUKrF27FqNHj9b62uIoqj/y6nzjjTcwc+ZMrW3r1q1b6vcvTT0FKei4fnGvefPy51NaJd0uyrLOBb22LMusKJ988gmGDx8ubdOTJ0+Wxl/VrFmz0NeWZL1Lus4qlQrt27eHra0tIiIi4O3tDXNzc/z555+YNWtWvkH5utpuXqRWq+Hk5IRNmzZpnZ+3h9jIyAjbtm3DqVOnsHv3buzbtw8jR47EJ598glOnTsHa2lrntZWVtv4aMGAATp48iRkzZsDPzw/W1tZQq9Xo2rVrsU7JL4/vkD6+K9euXUPnzp1Rv359fPrpp3Bzc4OZmRliYmKwYsWKYvXF0qVLMX/+fIwcORKLFy+Gg4MDjI2NMWXKlAq5vAHDTgWZOXOmxi7JFw8PFIejoyNsbGyQm5sr7bkoyLZt2+Dl5YXt27dr/IOat1emKF5eXgCen3Za1HsBgIODA0aMGIERI0bg8ePHaNeuHRYtWiSFnfK6kJ23tzceP35crBpf5OTkBHNzcyQkJOSbp21aSV29elXjr8SEhASo1WrpsI6HhwfUajWuXr0q/VUKPB+AqVKp4OHhUeYatCnrdlER8tb9ypUr6NSpk8a8K1euaPRNQduVh4cHDh48iEePHmns3ck7zFOc/vX19YWvry/mzZuHkydPok2bNli7di2WLFlS4nXSlaNHj+Lff//F9u3b0a5dO2l6YmJimZZbku+nt7c3Dh48iDZt2hQrTLVq1QqtWrXCBx98gM2bN2PIkCHYsmVLoX8IFfa5XrlyJd/04n6uJf0devjwIQ4dOoTw8HAsWLBAmq7rQ/Zl4eHhgUuXLkEIobF+xfkdc3R0hIWFhdb1ebmfd+/ejczMTPz0008ae6u0DYQuqJ+3bduGjh074quvvtKYrlKpUL169SLrLSsexqogPj4+CAwMlB7+/v4ler2JiQmCg4Px448/4sKFC/nm37t3T6MtoPkXwO+//47Y2NhivZeTkxM6dOiAL774AsnJyYW+18unYFpbW6N27doau6vzrlehUqmK9f7FNWDAAMTGxmLfvn355qlUKuTk5Gh9nYmJCQIDA7Fz507cuXNHmp6QkFDoeI7iyjslOs+qVasAQBoP1L17dwDId7bYp59+CuD5+JTyUNbtoiI0a9YMTk5OWLt2rcY29Msvv+Dy5csafVPQdtW9e3fk5ubi888/15i+YsUKGBkZFTguC3h+xtrL242vry+MjY3LdAhGF7R9fllZWVi9enWZlluS7+eAAQOQm5uLxYsX55uXk5MjLePhw4f59kDkne1XVD9aWVlpraV79+74448/NLbXJ0+eYN26dfD09ISPj0+pllsQbf0N5P/e6lNQUBBu376tcUmQjIwMfPnll0W+1sTEBEFBQdi5cyeSkpKk6ZcvX873m6qtL9LS0rBhw4Z8yy2on01MTPL15datWzUuKVGeuGdHBz7//HOoVCrpH87du3fj1q1bAJ4PTi3qdNXi+vDDD3HkyBG0bNkSY8aMgY+PDx48eIA///wTBw8elMYB9ezZE9u3b0ffvn3Ro0cPJCYmYu3atfDx8cHjx4+L9V5RUVF4/fXX4evrizFjxsDLywupqamIjY3FrVu3pGsj+Pj4oEOHDvD394eDgwPOnDmDbdu2adxLKi/YTZ48GUFBQTAxMcGgQYPK3B8zZszATz/9hJ49e0qncz558gTx8fHYtm0bbty4UeBfDIsWLcL+/fvRpk0bjB8/XvrHsVGjRmW+bH5iYiJ69+6Nrl27IjY2Ft999x3efvttNGnSBADQpEkThISEYN26ddKhiT/++AMbN25Enz590LFjxzK9f0F0sV2UN1NTU3z00UcYMWIE2rdvj8GDB0unnnt6emLq1KlS24K2q169eqFjx46YO3cubty4gSZNmmD//v3YtWsXpkyZku+U4xcdPnwYEydOxFtvvYW6desiJycH3377rfTHhj61bt0aVatWRUhICCZPngwjIyN8++23ZT6s4e3tDXt7e6xduxY2NjawsrJCy5YttY5had++PcaNG4fIyEjExcWhS5cuMDU1xdWrV7F161b85z//Qf/+/bFx40asXr0affv2hbe3Nx49eoQvv/wStra2UtgviL+/P9asWYMlS5agdu3acHJyQqdOnfD+++/j+++/R7du3TB58mQ4ODhg48aNSExMxI8//ljkFa39/f1x8OBBfPrpp3B1dUWtWrWkAfLa2Nraol27dli2bBmys7NRo0YN7N+/v8x70nRp3Lhx+PzzzzF48GC89957cHFxwaZNm6TB/UXtzQoPD8fevXvRtm1bTJgwATk5OVi1ahUaNmyoMfayS5cuMDMzQ69evTBu3Dg8fvwYX375JZycnPL9QVzQ59ezZ09ERERgxIgRaN26NeLj47Fp0ybpSEK5K/fzvV4BeacYanu8eCpiQQCI0NBQrct9+VTJ1NRUERoaKtzc3ISpqalQKpWic+fOYt26dVIbtVotli5dKjw8PIRCoRBNmzYVe/bsyXcaad4pxh9//LHWuq5duybeeecdoVQqhampqahRo4bo2bOn2LZtm9RmyZIlokWLFsLe3l5YWFiI+vXriw8++EBkZWVJbXJycsSkSZOEo6OjMDIyKvI0dA8PD62ns798yqwQz089nT17tqhdu7YwMzMT1atXF61btxbLly/XqAEvndYrhBCHDh0STZs2FWZmZsLb21usX79eTJs2TZibm2u0K+7nk3cK66VLl0T//v2FjY2NqFq1qpg4caLGadRCCJGdnS3Cw8NFrVq1hKmpqXBzcxOzZ8/WOH24sL4ojpdPAdXFdvFyPxZ06vnL/VXQMvNOu335VO8ffvhBNG3aVCgUCuHg4CCGDBkibt26pdGmsO3q0aNHYurUqcLV1VWYmpqKOnXqiI8//ljjtHBtrl+/LkaOHCm8vb2Fubm5cHBwEB07dhQHDx7UaFfQqecvn2qd1z/37t3TmB4SEiKsrKw0pr3ct9pOZ/7tt99Eq1athIWFhXB1dRUzZ84U+/bt03rqcsOGDbWuo7bv0a5du4SPj4906YW809Bf3jbyrFu3Tvj7+wsLCwthY2MjfH19xcyZM8WdO3eEEEL8+eefYvDgwcLd3V0oFArh5OQkevbsKc6cOaO1phelpKSIHj16CBsbGwFAo9Zr166J/v37C3t7e2Fubi5atGgh9uzZU+QyhRDi77//Fu3atRMWFhYCgPT5FfQZCSHErVu3RN++fYW9vb2ws7MTb731lrhz506xPqvi/o4VdOq5ts9P2+dx/fp10aNHD2FhYSEcHR3FtGnTxI8//igAiFOnThXZL8eOHRP+/v7CzMxMeHl5ibVr12r9Xv/000+icePGwtzcXHh6eoqPPvpIfP311/nWu6DPLyMjQ0ybNk24uLgICwsL0aZNGxEbG6t1eywPRkIY0MhAIj3r06dPuZxKT0RUUVauXImpU6fi1q1bqFGjhr7LMQgcs0OvrGfPnmk8v3r1KmJiYjQuT09EZMhe/h3LyMjAF198gTp16jDovIBjduiV5eXlJd335ubNm1izZg3MzMwKPJWdiMjQ9OvXD+7u7vDz80NaWhq+++47/P333wVeHuBVxbBDr6yuXbvi+++/R0pKChQKBQICArB06VKtF9kiIjJEQUFBWL9+PTZt2oTc3Fz4+Phgy5YtGDhwoL5LMygcs0NERESyxjE7REREJGsMO0RERCRrHLOD5/d7uXPnDmxsbMrt1gZERESkW0IIPHr0CK6uroVeWJJhB8CdO3fy3ZWYiIiIKod//vmn0Bv1MuwA0s0C//nnH9ja2uq5GiIiIiqO9PR0uLm5adz0VxuGHfzf/UNsbW0ZdoiIiCqZooagcIAyERERyRrDDhEREckaww4RERHJGsfsFFNubi6ys7P1XQa9IkxNTWFiYqLvMoiIZIFhpxgeP36MW7dugXfWoIpiZGSEmjVrwtraWt+lEBFVegw7RcjNzcWtW7dgaWkJR0dHXnSQyp0QAvfu3cOtW7dQp04d7uEhIiojhp0iZGdnQwgBR0dHWFhY6LscekU4Ojrixo0byM7OZtghIiojDlAuJu7RoYrE7Y2ISHcYdoiIiEjWeBirlJKSknD//v0Ke7/q1avD3d29wt6PiIhILhh2SiEpKQn1GzTAs6dPK+w9LSwt8ffly7IJPEePHkXHjh3x8OFD2Nvbl3o5np6emDJlCqZMmQLg+eGfHTt2oE+fPmWqr0OHDvDz88PKlSsLbBMdHY0pU6ZApVKV6b2IiKh8MeyUwv379/Hs6VMMmfUxnN29y/39UpOuYdNHM3D//v0ShZ3jx4/j448/xtmzZ5GcnKyTEKArrVu3RnJyMuzs7Mq0nNOnT8PKykpHVf2f7du3w9TUVHr+cqgCgIEDB6J79+46f28iItIthp0ycHb3Rs06DfVdRoGePHmCJk2aYOTIkejXr5++y9FgZmYGpVJZ5uU4OjrqoJr/k5WVBTMzMzg4OBTZ1sLCgmfoERFVAhygLGPdunXDkiVL0LdvX63zPT09sXTpUowcORI2NjZwd3fHunXrpPk3btyAkZERtm/fjo4dO8LS0hJNmjRBbGxssd7/5s2b6NWrF6pWrQorKys0bNgQMTExAJ4fxjIyMpIOAUVHR8Pe3h579uxBvXr1YGlpif79++Pp06fYuHEjPD09UbVqVUyePBm5ubka61DYoaZZs2ahbt26sLS0hJeXF+bPn69xJexFixbBz88P69evR61atWBubg7g+WGsvL04HTp0wM2bNzF16lQYGRlJZ0rl1fyiXbt24bXXXoO5uTm8vLwQHh6OnJwcAM+vn7No0SK4u7tDoVDA1dUVkydPLlZfEhFR6XHPzivuk08+weLFizFnzhxs27YN48ePR/v27VGvXj2pzdy5c7F8+XLUqVMHc+fOxeDBg5GQkIAqVQrffEJDQ5GVlYXjx4/DysoKly5dKvSKwE+fPsVnn32GLVu24NGjR+jXrx/69u0Le3t7xMTE4Pr16wgODkabNm0wcODAYq2fjY0NoqOj4erqivj4eIwZMwY2NjaYOXOm1CYhIQE//vgjtm/frvWaNtu3b0eTJk0wduxYjBkzpsD3+vXXX/HOO+/gs88+Q9u2bXHt2jWMHTsWALBw4UL8+OOPWLFiBbZs2YKGDRsiJSUFf/31V7HWg6g8+fo1RUpycqFtlC4uiI87V0EVEekWw84rrnv37pgwYQKA53tBVqxYgSNHjmiEnenTp6NHjx4AgPDwcDRs2BAJCQmoX79+octOSkpCcHAwfH19AQBeXl6Fts/OzsaaNWvg7f18HFT//v3x7bffIjU1FdbW1vDx8UHHjh1x5MiRYoedefPmSf/v6emJ6dOnY8uWLRphJysrC998802Bh8QcHBxgYmICGxubQg+9hYeH4/3330dISIi0vosXL8bMmTOxcOFCJCUlQalUIjAwEKampnB3d0eLFi2KtR5E5SklORlzvjteaJulQ9tVUDVEusfDWK+4xo0bS/9vZGQEpVKJu3fvFtjGxcUFAPK10Wby5MlYsmQJ2rRpg4ULF+L8+fOFtre0tJSCDgA4OzvD09NTY2+Qs7Nzsd47zw8//IA2bdpAqVTC2toa8+bNQ1JSkkYbDw8PnYz9+euvvxAREQFra2vpMWbMGCQnJ+Pp06d466238OzZM3h5eWHMmDHYsWOHdIiLiIjKD8POK+7FM46A54FHrVYX2CZvvMrLbbQZPXo0rl+/jmHDhiE+Ph7NmjXDqlWrSlRLceorSGxsLIYMGYLu3btjz549OHfuHObOnYusrCyNdro6m+vx48cIDw9HXFyc9IiPj8fVq1dhbm4ONzc3XLlyBatXr4aFhQUmTJiAdu3aaYwhIiIi3eNhLCpXbm5uePfdd/Huu+9i9uzZ+PLLLzFp0qQKee+TJ0/Cw8MDc+fOlabdvHmzVMsyMzPTGBitzWuvvYYrV66gdu3aBbaxsLBAr1690KtXL4SGhqJ+/fqIj4/Ha6+9Vqq6iIioaAw7ZZCadM2g3+fx48dISEiQnicmJiIuLg4ODg4VcnHCKVOmoFu3bqhbty4ePnyII0eOoEGDBuX+vnnq1KmDpKQkbNmyBc2bN8fPP/+MHTt2lGpZnp6eOH78OAYNGgSFQoHq1avna7NgwQL07NkT7u7u6N+/P4yNjfHXX3/hwoULWLJkCaKjo5Gbm4uWLVvC0tIS3333HSwsLODh4VHWVSUiokIw7JRC9erVYWFpiU0fzaiw97SwtNT6D2xhzpw5g44dO0rPw8LCAAAhISGIjo7WZXla5ebmIjQ0FLdu3YKtrS26du2KFStWlPv75unduzemTp2KiRMnIjMzEz169MD8+fOxaNGiEi8rIiIC48aNg7e3NzIzMyGEyNcmKCgIe/bsQUREBD766COYmpqifv36GD16NADA3t4eH374IcLCwpCbmwtfX1/s3r0b1apVK+uqEhFRIYyEtl/tV0x6ejrs7OyQlpYGW1tbjXkZGRlITEzUuAYLwHtjUfkqaLsjKg+OzspinY11LzWlgioiKp7C/v1+EffslJK7uzvDBxERUSXAs7Go1Lp166ZxmvWLj6VLl+q7PCIiIgDcs0NlsH79ejx79kzrvOLcW4qIiKgiMOxQqdWoUUPfJRBRJcNbU5A+MOwUE8dxU0Xi9kZyxVtTkD5wzE4R8m4M+fJVd4nKU972pu3GpEREVDLcs1OEKlWqwNLSEvfu3YOpqSmMjZkPqXyp1Wrcu3cPlpaWRd5ZnoiIisZf0iIYGRnBxcUFiYmJpb7VAFFJGRsbw93dXboXGRFp4tgfKgmGnWIwMzNDnTp1eCiLKoyZmRn3IhIVgmN/qCQYdorJ2NiYV7IlIiKqhPinIxEREckaww4RERHJGsMOERERyRrDDhEREckaww4RERHJGsMOERERyRrDDhEREckaww4RERHJml7Dzpo1a9C4cWPY2trC1tYWAQEB+OWXX6T5GRkZCA0NRbVq1WBtbY3g4GCkpqZqLCMpKQk9evSApaUlnJycMGPGDOTk5FT0qhAREZGB0mvYqVmzJj788EOcPXsWZ86cQadOnfDmm2/i4sWLAICpU6di9+7d2Lp1K44dO4Y7d+6gX79+0utzc3PRo0cPZGVl4eTJk9i4cSOio6OxYMECfa0SERERGRi93i6iV69eGs8/+OADrFmzBqdOnULNmjXx1VdfYfPmzejUqRMAYMOGDWjQoAFOnTqFVq1aYf/+/bh06RIOHjwIZ2dn+Pn5YfHixZg1axYWLVoEMzMzfawWERERGRCDGbOTm5uLLVu24MmTJwgICMDZs2eRnZ2NwMBAqU39+vXh7u6O2NhYAEBsbCx8fX3h7OwstQkKCkJ6erq0d0ibzMxMpKenazyIiEheVKo0ODorC334+jXVd5lUAfR+I9D4+HgEBAQgIyMD1tbW2LFjB3x8fBAXFwczMzPY29trtHd2dkZKSgoAICUlRSPo5M3Pm1eQyMhIhIeH63ZFiIjIoKjVat4ZnQAYwJ6devXqIS4uDr///jvGjx+PkJAQXLp0qVzfc/bs2UhLS5Me//zzT7m+HxEREemP3vfsmJmZoXbt2gAAf39/nD59Gv/5z38wcOBAZGVlQaVSaezdSU1NhVKpBAAolUr88ccfGsvLO1srr402CoUCCoVCx2tCRES6kHf4qag2RMWl97DzMrVajczMTPj7+8PU1BSHDh1CcHAwAODKlStISkpCQEAAACAgIAAffPAB7t69CycnJwDAgQMHYGtrCx8fH72tAxERlV5xDj9N7+5bQdWQHOg17MyePRvdunWDu7s7Hj16hM2bN+Po0aPYt28f7OzsMGrUKISFhcHBwQG2traYNGkSAgIC0KpVKwBAly5d4OPjg2HDhmHZsmVISUnBvHnzEBoayj03REREBEDPYefu3bt45513kJycDDs7OzRu3Bj79u3DG2+8AQBYsWIFjI2NERwcjMzMTAQFBWH16tXS601MTLBnzx6MHz8eAQEBsLKyQkhICCIiIvS1SkREslScQ0tKFxfEx52roIqIik+vYeerr74qdL65uTmioqIQFRVVYBsPDw/ExMToujQiInoBz2yiykzvZ2MRERERlSeGHSIiIpI1hh0iIiKSNYYdIiIikjWDu84OERFVTrwYIBkqhh0iItIJXgyQDBUPYxEREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrFXRdwFERET6olKlwdFZWWgbpYsL4uPOVVBFVB4YdoiI6JWlVqsx57vjhbZZOrRdBVVD5YWHsYiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1qrouwAiIiodX7+mSElOLrSN0sUF8XHnKqgiIsPEsENEVEmlJCdjznfHC22zdGi7CqqGyHDxMBYRERHJGsMOERERyZpew05kZCSaN28OGxsbODk5oU+fPrhy5YpGmw4dOsDIyEjj8e6772q0SUpKQo8ePWBpaQknJyfMmDEDOTk5FbkqREQGSaVKg6OzstCHSpWm7zKJypVex+wcO3YMoaGhaN68OXJycjBnzhx06dIFly5dgpWVldRuzJgxiIiIkJ5bWlpK/5+bm4sePXpAqVTi5MmTSE5OxjvvvANTU1MsXbq0QteHiMjQqNXqIsf1TO/uW0HVEOmHXsPO3r17NZ5HR0fDyckJZ8+eRbt2/zeoztLSEkqlUusy9u/fj0uXLuHgwYNwdnaGn58fFi9ejFmzZmHRokUwMzMr13UgIiJ6lVWGswIN6mystLTnu1IdHBw0pm/atAnfffcdlEolevXqhfnz50t7d2JjY+Hr6wtnZ2epfVBQEMaPH4+LFy+iadOm+d4nMzMTmZmZ0vP09PTyWB0iIiLZqwxnBRpM2FGr1ZgyZQratGmDRo0aSdPffvtteHh4wNXVFefPn8esWbNw5coVbN++HQCQkpKiEXQASM9TUlK0vldkZCTCw8PLaU2IiIjIkBhM2AkNDcWFCxdw4sQJjeljx46V/t/X1xcuLi7o3Lkzrl27Bm9v71K91+zZsxEWFiY9T09Ph5ubW+kKJyIiIoNmEKeeT5w4EXv27MGRI0dQs2bNQtu2bNkSAJCQkAAAUCqVSE1N1WiT97ygcT4KhQK2trYaDyIiIpInvYYdIQQmTpyIHTt24PDhw6hVq1aRr4mLiwMAuLi4AAACAgIQHx+Pu3fvSm0OHDgAW1tb+Pj4lEvdREREVHno9TBWaGgoNm/ejF27dsHGxkYaY2NnZwcLCwtcu3YNmzdvRvfu3VGtWjWcP38eU6dORbt27dC4cWMAQJcuXeDj44Nhw4Zh2bJlSElJwbx58xAaGgqFQqHP1SMiIiIDoNc9O2vWrEFaWho6dOgAFxcX6fHDDz8AAMzMzHDw4EF06dIF9evXx7Rp0xAcHIzdu3dLyzAxMcGePXtgYmKCgIAADB06FO+8847GdXmIiIjo1aXXPTtCiELnu7m54dixY0Uux8PDAzExMboqi4iIiGTEIAYoExEREZUXhh0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjWGHSIiIpK1KvougIiIyJCpVGlwdFYW2kbp4oL4uHMVVBGVFMMOERFRIdRqNeZ8d7zQNkuHtqugaqg0eBiLiIiIZI1hh4iIiGSNYYeIiIhkjWGHiIiIZI1hh4iIiGSNYYeIiIhkjWGHiIiIZI3X2SEiMkC+fk2RkpxcaBuVKq2CqiGq3Bh2iIgMUEpycpEXspve3beCqiGq3HgYi4iIiGRNr2EnMjISzZs3h42NDZycnNCnTx9cuXJFo01GRgZCQ0NRrVo1WFtbIzg4GKmpqRptkpKS0KNHD1haWsLJyQkzZsxATk5ORa4KERERGSi9hp1jx44hNDQUp06dwoEDB5CdnY0uXbrgyZMnUpupU6di9+7d2Lp1K44dO4Y7d+6gX79+0vzc3Fz06NEDWVlZOHnyJDZu3Ijo6GgsWLBAH6tEREREBkavY3b27t2r8Tw6OhpOTk44e/Ys2rVrh7S0NHz11VfYvHkzOnXqBADYsGEDGjRogFOnTqFVq1bYv38/Ll26hIMHD8LZ2Rl+fn5YvHgxZs2ahUWLFsHMzEwfq0ZEREQGwqDG7KSlPT+zwMHBAQBw9uxZZGdnIzAwUGpTv359uLu7IzY2FgAQGxsLX19fODs7S22CgoKQnp6Oixcvan2fzMxMpKenazyIiIhIngwm7KjVakyZMgVt2rRBo0aNAAApKSkwMzODvb29RltnZ2ekpKRIbV4MOnnz8+ZpExkZCTs7O+nh5uam47UhIiIiQ2EwYSc0NBQXLlzAli1byv29Zs+ejbS0NOnxzz//lPt7EhERkX4YxHV2Jk6ciD179uD48eOoWbOmNF2pVCIrKwsqlUpj705qaiqUSqXU5o8//tBYXt7ZWnltXqZQKKBQKHS8FkRERGSI9LpnRwiBiRMnYseOHTh8+DBq1aqlMd/f3x+mpqY4dOiQNO3KlStISkpCQEAAACAgIADx8fG4e/eu1ObAgQOwtbWFj49PxawIERERGSy97tkJDQ3F5s2bsWvXLtjY2EhjbOzs7GBhYQE7OzuMGjUKYWFhcHBwgK2tLSZNmoSAgAC0atUKANClSxf4+Phg2LBhWLZsGVJSUjBv3jyEhoZy7w0RERHpN+ysWbMGANChQweN6Rs2bMDw4cMBACtWrICxsTGCg4ORmZmJoKAgrF69WmprYmKCPXv2YPz48QgICICVlRVCQkIQERFRUatBREREBkyvYUcIUWQbc3NzREVFISoqqsA2Hh4eiImJ0WVpREREJBMGczYWERERUXlg2CEiIiJZY9ghIiIiWWPYISIiIlkziIsKEhG9Snz9miIlObnQNipVWgVVQyR/pQo7Xl5eOH36NKpVq6YxXaVS4bXXXsP169d1UhwRkRylJCdjznfHC20zvbtvBVVDJH+lOox148YN5Obm5puemZmJ27dvl7koIiIiIl0p0Z6dn376Sfr/ffv2wc7OTnqem5uLQ4cOwdPTU2fFEREREZVVicJOnz59AABGRkYICQnRmGdqagpPT0988sknOiuOiIiIqKxKFHbUajUAoFatWjh9+jSqV69eLkURERER6UqpBignJibqug4iIiKiclHqU88PHTqEQ4cO4e7du9Ienzxff/11mQsjIiIi0oVShZ3w8HBERESgWbNmcHFxgZGRka7rIiIiItKJUoWdtWvXIjo6GsOGDdN1PUREREQ6Vaqwk5WVhdatW+u6FiIiokpJpUqDo7Oy0DZKFxfEx52roIroRaUKO6NHj8bmzZsxf/58XddDRERU6ajV6iKvir10aLsKqoZeVqqwk5GRgXXr1uHgwYNo3LgxTE1NNeZ/+umnOimOiIiIqKxKFXbOnz8PPz8/AMCFCxc05nGwMhERERmSUoWdI0eO6LoOIiIionJRqhuBEhEREVUWpdqz07Fjx0IPVx0+fLjUBRERERHpUqnCTt54nTzZ2dmIi4vDhQsX8t0glIiIiEifShV2VqxYoXX6okWL8Pjx4zIVRERERKRLOh2zM3ToUN4Xi4iIiAyKTsNObGwszM3NdblIIiIiojIp1WGsfv36aTwXQiA5ORlnzpzhVZWJiIjIoJQq7NjZ2Wk8NzY2Rr169RAREYEuXbropDAiIiIiXShV2NmwYYOu6yAiIiIqF6UKO3nOnj2Ly5cvAwAaNmyIpk2b6qQoIiIiIl0pVdi5e/cuBg0ahKNHj8Le3h4AoFKp0LFjR2zZsgWOjo66rJGIiIio1Ep1NtakSZPw6NEjXLx4EQ8ePMCDBw9w4cIFpKenY/LkybqukYiIiKjUSrVnZ+/evTh48CAaNGggTfPx8UFUVBQHKBMREZFBKdWeHbVaDVNT03zTTU1NoVary1wUERERka6UKux06tQJ7733Hu7cuSNNu337NqZOnYrOnTvrrDgiIiKisipV2Pn888+Rnp4OT09PeHt7w9vbG7Vq1UJ6ejpWrVql6xqJiIiISq1UY3bc3Nzw559/4uDBg/j7778BAA0aNEBgYKBOiyMiIiIqqxKFncOHD2PixIk4deoUbG1t8cYbb+CNN94AAKSlpaFhw4ZYu3Yt2rZtWy7FEhHpk69fU6QkJxfaRunigvi4cxVUEREVR4nCzsqVKzFmzBjY2trmm2dnZ4dx48bh008/ZdghIllKSU7GnO+OF9pmZs8mcHRWFtpGpUrTZVlEVIQShZ2//voLH330UYHzu3TpguXLl5e5KCKiykqtVhcZiKZ3962gaogIKOEA5dTUVK2nnOepUqUK7t27V+aiiIiIiHSlRGGnRo0auHDhQoHzz58/DxcXlzIXRURERKQrJQo73bt3x/z585GRkZFv3rNnz7Bw4UL07NlTZ8URERERlVWJws68efPw4MED1K1bF8uWLcOuXbuwa9cufPTRR6hXrx4ePHiAuXPnFnt5x48fR69eveDq6gojIyPs3LlTY/7w4cNhZGSk8ejatatGmwcPHmDIkCGwtbWFvb09Ro0ahcePH5dktYiIiEjGSjRA2dnZGSdPnsT48eMxe/ZsCCEAAEZGRggKCkJUVBScnZ2LvbwnT56gSZMmGDlyJPr166e1TdeuXbFhwwbpuUKh0Jg/ZMgQJCcn48CBA8jOzsaIESMwduxYbN68uSSrRkRERDJV4osKenh4ICYmBg8fPkRCQgKEEKhTpw6qVq1a4jfv1q0bunXrVmgbhUIBpVL7aZyXL1/G3r17cfr0aTRr1gwAsGrVKnTv3h3Lly+Hq6triWsiIiIieSnV7SIAoGrVqmjevDlatGhRqqBTXEePHoWTkxPq1auH8ePH499//5XmxcbGwt7eXgo6ABAYGAhjY2P8/vvvBS4zMzMT6enpGg8iIiKSp1KHnYrQtWtXfPPNNzh06BA++ugjHDt2DN26dUNubi4AICUlBU5OThqvqVKlChwcHJCSklLgciMjI2FnZyc93NzcynU9iIiISH9KdW+sijJo0CDp/319fdG4cWN4e3vj6NGjZbq7+uzZsxEWFiY9T09PZ+AhIiKSKYMOOy/z8vJC9erVkZCQgM6dO0OpVOLu3bsabXJycvDgwYMCx/kAz8cBvTzQmYiIqDypVGlF3kqE91YrH5Uq7Ny6dQv//vuvdOHCgIAAqFQqnD17Fv7+/gCe36xUrVajZcuW+iyViIhIQ3FuJbJ0aLsKqubVotew8/jxYyQkJEjPExMTERcXBwcHBzg4OCA8PBzBwcFQKpW4du0aZs6cidq1ayMoKAgA0KBBA3Tt2hVjxozB2rVrkZ2djYkTJ2LQoEE8E4uISqQ4dzTnDTyJKie9hp0zZ86gY8eO0vO8cTQhISFYs2YNzp8/j40bN0KlUsHV1RVdunTB4sWLNQ5Bbdq0CRMnTkTnzp1hbGyM4OBgfPbZZxW+LkRUuRXnjua8gSdR5aTXsNOhQwfpwoTa7Nu3r8hlODg48AKCREREVCCDPvWciIiIqKwYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1vR6I1AiosL4+jVFSnJyoW2ULi6IjztXQRURUWXEsENEBislORlzvjteaJulQ9tVUDVEVFnxMBYRERHJGsMOERERyRrDDhEREckax+wQkewVZ6CzSpVWQdUQUUVj2CEi2SvOQOfp3X0rqBoiqmg8jEVERESyxj07RAasOIdfAF5rhoioMAw7RAasOIdfAF5rhoioMDyMRURERLLGPTtEesIzhIiIKgbDDpGeyPkMId7TiogMCcMOEZVIcfdILdvzV6FtOM6IiCoKww4RlYic90gRkTxxgDIRERHJGsMOERERyRrDDhEREckaww4RERHJGsMOERERyRrPxiJ6RfDaN0T0qmLYIXpFFOeUcV77hojkiIexiIiISNYYdoiIiEjWeBiLqISKM/bl8ZMnsLayKrSNId7kU6VKg6Ozssg2RESVCcMOUQkV93YJc3acLbKNoVGr1bwVBBHJDg9jERERkawx7BAREZGsMewQERGRrHHMDhHpha4GQ3NQNREVhWGHiPRCV4OhOaiaiIqi18NYx48fR69eveDq6gojIyPs3LlTY74QAgsWLICLiwssLCwQGBiIq1evarR58OABhgwZAltbW9jb22PUqFF4/PhxBa4FERERGTK9hp0nT56gSZMmiIqK0jp/2bJl+Oyzz7B27Vr8/vvvsLKyQlBQEDIyMqQ2Q4YMwcWLF3HgwAHs2bMHx48fx9ixYytqFYiIiMjA6fUwVrdu3dCtWzet84QQWLlyJebNm4c333wTAPDNN9/A2dkZO3fuxKBBg3D58mXs3bsXp0+fRrNmzQAAq1atQvfu3bF8+XK4urpW2LoQ6RPHrRARFcxgx+wkJiYiJSUFgYGB0jQ7Ozu0bNkSsbGxGDRoEGJjY2Fvby8FHQAIDAyEsbExfv/9d/Tt21frsjMzM5GZmSk9T09PL78VIaoAHLdCRFQwgz31PCUlBQDg7OysMd3Z2Vmal5KSAicnJ435VapUgYODg9RGm8jISNjZ2UkPNzc3HVdPREREhsJgw055mj17NtLS0qTHP//8o++SiIiIqJwYbNhRKp+PP0hNTdWYnpqaKs1TKpW4e/euxvycnBw8ePBAaqONQqGAra2txoOIiIjkyWDDTq1ataBUKnHo0CFpWnp6On7//XcEBAQAAAICAqBSqXD27P/dcPHw4cNQq9Vo2bJlhddMREREhkevA5QfP36MhIQE6XliYiLi4uLg4OAAd3d3TJkyBUuWLEGdOnVQq1YtzJ8/H66urujTpw8AoEGDBujatSvGjBmDtWvXIjs7GxMnTsSgQYN4JhYREREB0HPYOXPmDDp27Cg9DwsLAwCEhIQgOjoaM2fOxJMnTzB27FioVCq8/vrr2Lt3L8zNzaXXbNq0CRMnTkTnzp1hbGyM4OBgfPbZZxW+LkRERGSY9Bp2OnToACFEgfONjIwQERGBiIiIAts4ODhg8+bN5VEeERERyYDBjtkhIiIi0gWGHSIiIpI1hh0iIiKSNYYdIiIikjWDvTcWERHRq6Y4N/VVurggPu5cBVUkDww7REREBqI4N/VdOrRdBVUjHzyMRURERLLGsENERESyxrBDREREssawQ0RERLLGsENERESyxrBDREREssawQ0RERLLGsENERESyxrBDREREssawQ0RERLLGsENERESyxrBDREREssawQ0RERLLGsENERESyxrBDREREssawQ0RERLJWRd8FEBkSX7+mSElOLrSNSpVWQdUQEZEuMOwQvSAlORlzvjteaJvp3X0rqBoiItIFHsYiIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWePtIoiIiCoRlSoNjs7KQts8fvIE1lZWhbZRurggPu6cLkszWAw7RERElYharS7WPfzm7DhbaJulQ9vpsiyDxsNYREREJGsMO0RERCRrDDtEREQkaww7REREJGsMO0RERCRrDDtEREQkaww7REREJGsGHXYWLVoEIyMjjUf9+vWl+RkZGQgNDUW1atVgbW2N4OBgpKam6rFiIiIiMjQGHXYAoGHDhkhOTpYeJ06ckOZNnToVu3fvxtatW3Hs2DHcuXMH/fr102O1REREZGgM/grKVapUgVKZ/7LYaWlp+Oqrr7B582Z06tQJALBhwwY0aNAAp06dQqtWrSq6VCIiIjJABh92rl69CldXV5ibmyMgIACRkZFwd3fH2bNnkZ2djcDAQKlt/fr14e7ujtjY2ELDTmZmJjIzM6Xn6enp5boOZBh8/ZoiJTm50DYqVVoFVUNEZPjk8rtp0GGnZcuWiI6ORr169ZCcnIzw8HC0bdsWFy5cQEpKCszMzGBvb6/xGmdnZ6SkpBS63MjISISHh5dj5WSIUpKTi3U/GSIiek4uv5sGHXa6desm/X/jxo3RsmVLeHh44L///S8sLCxKvdzZs2cjLCxMep6eng43N7cy1UpERESGyeAHKL/I3t4edevWRUJCApRKJbKysqBSqTTapKamah3j8yKFQgFbW1uNBxEREclTpQo7jx8/xrVr1+Di4gJ/f3+Ympri0KFD0vwrV64gKSkJAQEBeqySiIiIDIlBH8aaPn06evXqBQ8PD9y5cwcLFy6EiYkJBg8eDDs7O4waNQphYWFwcHCAra0tJk2ahICAAJ6JRURERBKDDju3bt3C4MGD8e+//8LR0RGvv/46Tp06BUdHRwDAihUrYGxsjODgYGRmZiIoKAirV6/Wc9VERERkSAw67GzZsqXQ+ebm5oiKikJUVFQFVURERESVTaUas0NERERUUgw7REREJGsGfRiLqLjkcpVPIiLSPYYdkgW5XOWTiIh0j4exiIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1njqORER0StIpUqDo7OyyDZywLBDRET0ClKr1a/M9cl4GIuIiIhkjWGHiIiIZI1hh4iIiGSNYYeIiIhkjWGHiIiIZI1hh4iIiGSNYYeIiIhkjWGHiIiIZI1hh4iIiGSNV1Amg+fr1xQpycmFtpHLJc2JiEj3GHbI4KUkJ78ylzQnIiLd42EsIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjXeLoL0ive9IiKi8sawQ3rF+14REVF542EsIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjXZhJ2oqCh4enrC3NwcLVu2xB9//KHvkl55vn5N4eisLPTBCwYSEVF5k8VFBX/44QeEhYVh7dq1aNmyJVauXImgoCBcuXIFTk5O+i7vlcULBhIRkSGQxZ6dTz/9FGPGjMGIESPg4+ODtWvXwtLSEl9//bW+SyMiIiI9q/R7drKysnD27FnMnj1bmmZsbIzAwEDExsZqfU1mZiYyMzOl52lpzw+lpKen67y+Vm1eR2pKSqFtnJVKnPrtRKV6r+JQq9XIePK40DZCCLYpYxtDrIlt2IZt2OZFarW6XP6NzVumEKLwhqKSu337tgAgTp48qTF9xowZokWLFlpfs3DhQgGADz744IMPPviQweOff/4pNCtU+j07pTF79myEhYVJz9VqNW7evAk/Pz/8888/sLW11WN1r5b09HS4ubmx3/WAfa8f7Hf9Yd/rR3n2uxACjx49gqura6HtKn3YqV69OkxMTJCamqoxPTU1FUqlUutrFAoFFAqFxjRj4+fDl2xtbfkl0AP2u/6w7/WD/a4/7Hv9KK9+t7OzK7JNpR+gbGZmBn9/fxw6dEiaplarcejQIQQEBOixMiIiIjIElX7PDgCEhYUhJCQEzZo1Q4sWLbBy5Uo8efIEI0aM0HdpREREpGeyCDsDBw7EvXv3sGDBAqSkpMDPzw979+6Fs7NzsZehUCiwcOHCfIe3qHyx3/WHfa8f7Hf9Yd/rhyH0u5EQRZ2vRURERFR5VfoxO0RERESFYdghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWXulwk5UVBQ8PT1hbm6Oli1b4o8//ii0/datW1G/fn2Ym5vD19cXMTExFVSpvJSk37/88ku0bdsWVatWRdWqVREYGFjk50QFK+k2n2fLli0wMjJCnz59yrdAmSppv6tUKoSGhsLFxQUKhQJ169bl700plbTvV65ciXr16sHCwgJubm6YOnUqMjIyKqhaeTh+/Dh69eoFV1dXGBkZYefOnUW+5ujRo3jttdegUChQu3ZtREdHl2+Rurkdp+HbsmWLMDMzE19//bW4ePGiGDNmjLC3txepqala2//222/CxMRELFu2TFy6dEnMmzdPmJqaivj4+AquvHIrab+//fbbIioqSpw7d05cvnxZDB8+XNjZ2Ylbt25VcOWVX0n7Pk9iYqKoUaOGaNu2rXjzzTcrplgZKWm/Z2ZmimbNmonu3buLEydOiMTERHH06FERFxdXwZVXfiXt+02bNgmFQiE2bdokEhMTxb59+4SLi4uYOnVqBVdeucXExIi5c+eK7du3CwBix44dhba/fv26sLS0FGFhYeLSpUti1apVwsTEROzdu7fcanxlwk6LFi1EaGio9Dw3N1e4urqKyMhIre0HDBggevTooTGtZcuWYty4ceVap9yUtN9flpOTI2xsbMTGjRvLq0TZKk3f5+TkiNatW4v169eLkJAQhp1SKGm/r1mzRnh5eYmsrKyKKlG2Str3oaGholOnThrTwsLCRJs2bcq1TjkrTtiZOXOmaNiwoca0gQMHiqCgoHKr65U4jJWVlYWzZ88iMDBQmmZsbIzAwEDExsZqfU1sbKxGewAICgoqsD3lV5p+f9nTp0+RnZ0NBweH8ipTlkrb9xEREXBycsKoUaMqokzZKU2///TTTwgICEBoaCicnZ3RqFEjLF26FLm5uRVVtiyUpu9bt26Ns2fPSoe6rl+/jpiYGHTv3r1Can5V6ePfV1ncLqIo9+/fR25ubr7bRzg7O+Pvv//W+pqUlBSt7VNSUsqtTrkpTb+/bNasWXB1dc33xaDClabvT5w4ga+++gpxcXEVUKE8labfr1+/jsOHD2PIkCGIiYlBQkICJkyYgOzsbCxcuLAiypaF0vT922+/jfv37+P111+HEAI5OTl49913MWfOnIoo+ZVV0L+v6enpePbsGSwsLHT+nq/Enh2qnD788ENs2bIFO3bsgLm5ub7LkbVHjx5h2LBh+PLLL1G9enV9l/NKUavVcHJywrp16+Dv74+BAwdi7ty5WLt2rb5Lk72jR49i6dKlWL16Nf78809s374dP//8MxYvXqzv0kjHXok9O9WrV4eJiQlSU1M1pqempkKpVGp9jVKpLFF7yq80/Z5n+fLl+PDDD3Hw4EE0bty4PMuUpZL2/bVr13Djxg306tVLmqZWqwEAVapUwZUrV+Dt7V2+RctAabZ5FxcXmJqawsTERJrWoEEDpKSkICsrC2ZmZuVas1yUpu/nz5+PYcOGYfTo0QAAX19fPHnyBGPHjsXcuXNhbMz9AeWhoH9fbW1ty2WvDvCK7NkxMzODv78/Dh06JE1Tq9U4dOgQAgICtL4mICBAoz0AHDhwoMD2lF9p+h0Ali1bhsWLF2Pv3r1o1qxZRZQqOyXt+/r16yM+Ph5xcXHSo3fv3ujYsSPi4uLg5uZWkeVXWqXZ5tu0aYOEhAQpXALA//73P7i4uDDolEBp+v7p06f5Ak1e6BS8R3a50cu/r+U29NnAbNmyRSgUChEdHS0uXbokxo4dK+zt7UVKSooQQohhw4aJ999/X2r/22+/iSpVqojly5eLy5cvi4ULF/LU81Ioab9/+OGHwszMTGzbtk0kJydLj0ePHulrFSqtkvb9y3g2VumUtN+TkpKEjY2NmDhxorhy5YrYs2ePcHJyEkuWLNHXKlRaJe37hQsXChsbG/H999+L69evi/379wtvb28xYMAAfa1CpfTo0SNx7tw5ce7cOQFAfPrpp+LcuXPi5s2bQggh3n//fTFs2DCpfd6p5zNmzBCXL18WUVFRPPVcl1atWiXc3d2FmZmZaNGihTh16pQ0r3379iIkJESj/X//+19Rt25dYWZmJho2bCh+/vnnCq5YHkrS7x4eHgJAvsfChQsrvnAZKOk2/yKGndIrab+fPHlStGzZUigUCuHl5SU++OADkZOTU8FVy0NJ+j47O1ssWrRIeHt7C3Nzc+Hm5iYmTJggHj58WPGFV2JHjhzR+rud19chISGiffv2+V7j5+cnzMzMhJeXl9iwYUO51mgkBPfVERERkXy9EmN2iIiI6NXFsENERESyxrBDREREssawQ0RERLLGsENERESyxrBDREREssawQ0RERLLGsENERESyxrBDREREssawQ0RERLLGsENERESy9v8AkkltvKuJz34AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "sns.histplot(pd.DataFrame({\"1nn_similarities\": similarities}), bins=50)\n", + "plt.title(f\"1-nearest neighbor Tanimoto similarities to training data\")" + ] + }, + { + "cell_type": "markdown", + "id": "ba609b94-1a37-4691-98cb-48533a9fa5c5", + "metadata": {}, + "source": [ + "As the histogram shows, the similarity between the test and training set is relatively high with most compounds having a similarity >0.6 and even ~250 molecules with a Tanimoto score of 1. However, this is just a hypothetical example. If a real-world dataset would have such high similarities we would probably use cluster or time splits to reduce the similarity and data leakage. " + ] + }, + { + "cell_type": "markdown", + "id": "4439964e-d9e7-4e5c-af77-10e16db7d842", + "metadata": {}, + "source": [ + "## Comparison to native RDKit Tanimoto computation" + ] + }, + { + "cell_type": "markdown", + "id": "61516a5d-a871-4b30-a5eb-3fe57282a9d1", + "metadata": {}, + "source": [ + "`NearestNeighborsRetrieverTanimoto` performs an exhaustive comparison to find the k-nearest neighbors. To do this, the full similarity matrix must be computed. MolPipeline's algorithm for finding these Tanimoto similarity scores differs from the approach in RDKit. In MolPipeline, we use an implementation based on sparse matrices that exploits the sparse matrix dot product algorithm from scipy. The central function is `tanimoto_similarity_sparse` which computes the full similarity matrix." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "d54ef845-9b6b-474e-85a8-4368b69387cf", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 19.3 s, sys: 9.39 s, total: 28.7 s\n", + "Wall time: 28.7 s\n" + ] + }, + { + "data": { + "text/plain": [ + "(19999, 19999)" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%%time\n", + "sim_matrix = tanimoto_similarity_sparse(fp_matrix, fp_matrix)\n", + "sim_matrix.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "ab1a9fdf-ad68-4c53-9d66-5b26c1c4694b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1. , 0.08783784, 0.05747126, ..., 0.09836066, 0.06730769,\n", + " 0.08602151],\n", + " [0.08783784, 1. , 0.45212766, ..., 0.05405405, 0.19760479,\n", + " 0.10465116],\n", + " [0.05747126, 0.45212766, 1. , ..., 0.04678363, 0.17368421,\n", + " 0.09230769],\n", + " ...,\n", + " [0.09836066, 0.05405405, 0.04678363, ..., 1. , 0.07070707,\n", + " 0.10344828],\n", + " [0.06730769, 0.19760479, 0.17368421, ..., 0.07070707, 1. ,\n", + " 0.12903226],\n", + " [0.08602151, 0.10465116, 0.09230769, ..., 0.10344828, 0.12903226,\n", + " 1. ]])" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sim_matrix" + ] + }, + { + "cell_type": "markdown", + "id": "7f588705-1d0e-48d8-b5df-8f987298425a", + "metadata": {}, + "source": [ + "To get the full similarity matrix with RDKit using `BulkTanimotoSimilarity`, we have to have the fingerprints as a different datastructure, for example as `ExplicitBitVect`. " + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "c53e1b2f-3001-49ba-bc82-152098facd0d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 1.27 s, sys: 143 ms, total: 1.42 s\n", + "Wall time: 1.73 s\n" + ] + }, + { + "data": { + "text/plain": [ + "[,\n", + " ,\n", + " ,\n", + " ]" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%%time\n", + "error_filter = ErrorFilter()\n", + "\n", + "fingerprint_pipeline2 = Pipeline(\n", + " [\n", + " (\"auto2mol\", AutoToMol()),\n", + " (\"error_filter\", error_filter),\n", + " (\n", + " \"morgan2_2048\",\n", + " MolToMorganFP(n_bits=2048, radius=2, return_as=\"explicit_bit_vect\"),\n", + " ),\n", + " ],\n", + " n_jobs=-1,\n", + ")\n", + "\n", + "fp_matrix_explicit = fingerprint_pipeline2.transform(df[\"smiles\"])\n", + "fp_matrix_explicit[:4]" + ] + }, + { + "cell_type": "markdown", + "id": "83e28e13-b375-48d1-8160-2a9b08edba73", + "metadata": {}, + "source": [ + "Now, let's compute the full similarity matrix using RDKit's `BulkTanimotoSimilarity`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "797df24b-d431-47d5-860c-7eb7bd933564", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "sim_mat_rdkit = np.full((len(fp_matrix_explicit), len(fp_matrix_explicit)), np.nan)\n", + "for i, query_fp in enumerate(fp_matrix_explicit):\n", + " sim_mat_rdkit[i, :] = DataStructs.BulkTanimotoSimilarity(\n", + " query_fp, fp_matrix_explicit\n", + " )\n", + "sim_mat_rdkit.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "92b5b3fa-7365-439a-a490-877156f039a8", + "metadata": {}, + "outputs": [], + "source": [ + "assert np.allclose(sim_matrix, sim_mat_rdkit)" + ] + }, + { + "cell_type": "markdown", + "id": "442e7536-2e67-4f80-938c-12b958ac2d89", + "metadata": {}, + "source": [ + "Based on this simple comparison MolPipeline's similarity matrix computation is about ~2-3 times faster than RDKit's. However, of course there are many other things to consider that are not touched in this notebook. For example, `tanimoto_similarity_sparse` uses more memory since it needs intermediate matrices while `BulkTanimotoSimilarity` uses almost no memory. In addition, for both approaches different strategies for parallelization come to mind (one is implemented in `NearestNeighborsRetrieverTanimoto`), which can be beneficial in different scenarios. Lastly, while the here discussed functions are useful for easy analysis in Python, there are highly optimized tools for similarity search, like [Artor](https://www.nextmovesoftware.com/arthor.html) which should probably be used when search speed is essential. " + ] + } + ], + "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.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tests/test_estimators/test_nearest_neighbors.py b/tests/test_estimators/test_nearest_neighbors.py index 88bc14fa..cda24c8f 100644 --- a/tests/test_estimators/test_nearest_neighbors.py +++ b/tests/test_estimators/test_nearest_neighbors.py @@ -3,11 +3,13 @@ from unittest import TestCase import numpy as np +from scipy.sparse import csr_matrix from sklearn.base import clone from molpipeline import ErrorFilter, FilterReinserter, Pipeline, PostPredictionWrapper from molpipeline.any2mol import SmilesToMol from molpipeline.estimators import NamedNearestNeighbors, TanimotoToTraining +from molpipeline.estimators.nearest_neighbor import NearestNeighborsRetrieverTanimoto from molpipeline.mol2any import MolToMorganFP from molpipeline.utils.kernel import tanimoto_distance_sparse @@ -34,6 +36,15 @@ ] ) +FOUR_NN_SIMILARITIES = np.array( + [ + [1.0, 3 / 14, 0.0, 0.0], + [1.0, 3 / 14, 0.038461538461538464, 0.0], + [1.0, 4 / 9, 0.0, 0.0], + [1.0, 4 / 9, 0.038461538461538464, 0.0], + ] +) + class TestNamedNearestNeighbors(TestCase): """Test the NamedNearestNeighbors class if correct names are returned.""" @@ -209,3 +220,95 @@ def test_fit_and_predict_invalid_with_distance(self) -> None: self.assertTrue( 1 - np.allclose(distances[1:, :].astype(np.float64), TWO_NN_SIMILARITIES) ) + + +class TestNearestNeighborsRetrieverTanimoto(TestCase): + """Test nearest neighbors retriever with tanimoto.""" + + example_fingerprints: csr_matrix + + @classmethod + def setUpClass(cls) -> None: + """Set up the tests.""" + morgan_pipeline = Pipeline( + [ + ("mol", SmilesToMol()), + ("fingerprint", MolToMorganFP()), + ] + ) + cls.example_fingerprints = morgan_pipeline.transform(TEST_SMILES) + + def test_k_equals_1(self) -> None: + """Test the k=1 retrieval.""" + target_fps = self.example_fingerprints + query_fps = self.example_fingerprints + + retriever = NearestNeighborsRetrieverTanimoto(target_fps, k=1) + indices, similarities = retriever.predict(query_fps) + self.assertTrue(np.array_equal(indices, np.array([0, 1, 2, 3]))) + self.assertTrue(np.allclose(similarities, np.array([1, 1, 1, 1]))) + + # test parallel + retriever = NearestNeighborsRetrieverTanimoto( + target_fps, k=1, n_jobs=2, batch_size=2 + ) + indices, similarities = retriever.predict(query_fps) + self.assertTrue(np.array_equal(indices, np.array([0, 1, 2, 3]))) + self.assertTrue(np.allclose(similarities, np.array([1, 1, 1, 1]))) + + def test_k_greater_1_less_n(self) -> None: + """Test k>1 and k None: + """Test k=n retrieval.""" + target_fps = self.example_fingerprints + query_fps = self.example_fingerprints + + retriever = NearestNeighborsRetrieverTanimoto(target_fps, k=target_fps.shape[0]) + indices, similarities = retriever.predict(query_fps) + self.assertTrue( + np.array_equal( + indices, + np.array([[0, 1, 3, 2], [1, 0, 3, 2], [2, 3, 1, 0], [3, 2, 1, 0]]), + ) + ) + self.assertTrue(np.allclose(similarities, FOUR_NN_SIMILARITIES)) + + # test parallel + retriever = NearestNeighborsRetrieverTanimoto( + target_fps, k=target_fps.shape[0], n_jobs=2, batch_size=2 + ) + indices, similarities = retriever.predict(query_fps) + self.assertTrue( + np.array_equal( + indices, + np.array([[0, 1, 3, 2], [1, 0, 3, 2], [2, 3, 1, 0], [3, 2, 1, 0]]), + ) + ) + self.assertTrue(np.allclose(similarities, FOUR_NN_SIMILARITIES)) + + # [ + # [1.0, 3 / 14, 0.0, 0.0], + # [1.0, 3 / 14, 0.038461538461538464, 0.0], + # [1.0, 4 / 9, 0.0, 0.0], + # [1.0, 4 / 9, 0.038461538461538464, 0.0], + # ]