esm_2_clustering.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clustering Proteins with ESM-2 and Persistent Homology"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook will focus on analyzing protein sequences using a transformer model called ESM-2 from Facebook AI. We'll use the attention matrices from the transformer model to compute Jensen-Shannon distance matrices. Then, we'll compute persistent homology diagrams using the Gudhi library. Finally, we'll cluster these diagrams using various clustering algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's start by importing necessary libraries\n",
    "from transformers import EsmModel, AutoTokenizer\n",
    "import gudhi as gd\n",
    "import numpy as np\n",
    "import torch\n",
    "from scipy.spatial.distance import pdist, jensenshannon, squareform\n",
    "from sklearn.cluster import KMeans, AgglomerativeClustering, DBSCAN\n",
    "from sklearn.metrics import silhouette_score\n",
    "from gudhi.hera import wasserstein_distance"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's define some helper functions for our analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper function to get the attention matrix of a specific layer and head for a given input sequence\n",
    "def get_attention_matrix(tokenizer, model, layer, head, input_sequence):\n",
    "    model.config.output_attentions = True\n",
    "    encoded_input = tokenizer([input_sequence], return_tensors='pt', padding=True, truncation=True)\n",
    "    with torch.no_grad():\n",
    "        model_output = model(**encoded_input)\n",
    "    attention_weights = model_output.attentions\n",
    "    specific_attention_weights = attention_weights[layer][0, head]\n",
    "    return specific_attention_weights\n",
    "\n",
    "# Helper function to compute the Jensen-Shannon distance matrix\n",
    "def compute_js_distance_matrix_scipy(attention_matrix):\n",
    "    js_distances = pdist(attention_matrix, metric=jensenshannon)\n",
    "    js_distance_matrix = squareform(js_distances)\n",
    "    return js_distance_matrix\n",
    "\n",
    "# Helper function to compute the persistent homology of a given distance matrix\n",
    "def compute_persistent_homology(distance_matrix, max_dimension=3):\n",
    "    max_edge_length = np.max(distance_matrix)\n",
    "    rips_complex = gd.RipsComplex(distance_matrix=distance_matrix, max_edge_length=max_edge_length)\n",
    "    st = rips_complex.create_simplex_tree(max_dimension=max_dimension)\n",
    "    persistence = st.persistence()\n",
    "    return st, persistence\n",
    "\n",
    "# Helper function to compute the Wasserstein distances between all pairs of persistence diagrams\n",
    "def compute_wasserstein_distances(persistence_diagrams, dimension):\n",
    "    n_diagrams = len(persistence_diagrams)\n",
    "    distances = np.zeros((n_diagrams, n_diagrams))\n",
    "    filtered_diagrams = [[point for point in diagram if point[0] == dimension] for diagram in persistence_diagrams]\n",
    "    for i in range(n_diagrams):\n",
    "        for j in range(i+1, n_diagrams):\n",
    "            X = np.array([p[1] for p in filtered_diagrams[i]])\n",
    "            Y = np.array([p[1] for p in filtered_diagrams[j]])\n",
    "            distance = wasserstein_distance(X, Y)\n",
    "            distances[i][j] = distance\n",
    "            distances[j][i] = distance\n",
    "    return distances\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our helper functions ready, let's now load our transformer model and tokenizer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Some weights of the model checkpoint at facebook/esm2_t6_8M_UR50D were not used when initializing EsmModel: ['lm_head.dense.weight', 'lm_head.dense.bias', 'lm_head.bias', 'lm_head.layer_norm.bias', 'lm_head.layer_norm.weight']\n",
      "- This IS expected if you are initializing EsmModel from the checkpoint of a model trained on another task or with another architecture (e.g. initializing a BertForSequenceClassification model from a BertForPreTraining model).\n",
      "- This IS NOT expected if you are initializing EsmModel from the checkpoint of a model that you expect to be exactly identical (initializing a BertForSequenceClassification model from a BertForSequenceClassification model).\n",
      "Some weights of EsmModel were not initialized from the model checkpoint at facebook/esm2_t6_8M_UR50D and are newly initialized: ['esm.pooler.dense.weight', 'esm.pooler.dense.bias']\n",
      "You should probably TRAIN this model on a down-stream task to be able to use it for predictions and inference.\n"
     ]
    }
   ],
   "source": [
    "# Load the tokenizer and model\n",
    "tokenizer = AutoTokenizer.from_pretrained(\"facebook/esm2_t6_8M_UR50D\")\n",
    "model = EsmModel.from_pretrained(\"facebook/esm2_t6_8M_UR50D\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute the attention matrix, Jensen-Shannon distance matrix, and persistent homology for a sample protein sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define layer and head to be used\n",
    "layer = 3\n",
    "head = 2\n",
    "\n",
    "# Define a sample protein sequence\n",
    "input_sequence = \"MGWGRKRR\"\n",
    "\n",
    "# Compute the attention matrix\n",
    "attention_matrix = get_attention_matrix(tokenizer, model, layer, head, input_sequence)\n",
    "\n",
    "# Compute the Jensen-Shannon distance matrix\n",
    "distance_matrix = compute_js_distance_matrix_scipy(attention_matrix)\n",
    "\n",
    "# Compute the persistent homology\n",
    "st, persistence = compute_persistent_homology(distance_matrix, max_dimension=4)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've computed the persistent homology for our sample sequence, let's compute it for a list of protein sequences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a list of protein sequences\n",
    "sequences = [\n",
    "    # Protein family 1: Hemoglobin\n",
    "    \"VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFF\",\n",
    "    \"VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFES\",\n",
    "    \"VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESF\",\n",
    "\n",
    "    # Protein family 2: Insulin\n",
    "    \"MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGER\",\n",
    "    \"MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERG\",\n",
    "    \"MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGF\",\n",
    "\n",
    "    # Protein family 3: Albumin\n",
    "    \"MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGER\",\n",
    "    \"MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERG\",\n",
    "    \"MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGF\",\n",
    "\n",
    "    # Protein family 4: Myoglobin\n",
    "    \"GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDR\",\n",
    "    \"GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDRF\",\n",
    "    \"GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDRFG\",\n",
    "\n",
    "    # Protein family 5: Lysozyme\n",
    "    \"KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNR\",\n",
    "    \"KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNRN\",\n",
    "    \"KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNRNT\"\n",
    "]\n",
    "\n",
    "\n",
    "# Initialize list to store persistent diagrams\n",
    "persistent_diagrams = []\n",
    "\n",
    "# Compute persistent homology for each sequence\n",
    "for sequence in sequences:\n",
    "    attention_matrix = get_attention_matrix(tokenizer, model, layer, head, sequence)\n",
    "    distance_matrix = compute_js_distance_matrix_scipy(attention_matrix.numpy())\n",
    "    _, persistence_diagram = compute_persistent_homology(distance_matrix)\n",
    "    \n",
    "    # Store the persistent diagram\n",
    "    persistent_diagrams.append(persistence_diagram)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With persistent diagrams in hand, we can now compute the Wasserstein distances between all pairs of persistence diagrams, which provide a measure of similarity between the diagrams."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the Wasserstein distances between all pairs of persistence diagrams\n",
    "wasserstein_distances = compute_wasserstein_distances(persistent_diagrams, 0)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, compute the persistent homology of the Wasserstein distance matrix. We already have a helper function compute_persistent_homology for this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the persistent homology of the Wasserstein distance matrix\n",
    "st_2, persistence_2 = compute_persistent_homology(wasserstein_distances)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, plot the persistence diagram using the Gudhi library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: title={'center': 'Persistence diagram'}, xlabel='Birth', ylabel='Death'>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAHMCAYAAAAwHmdPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABXBUlEQVR4nO3deXwTdf4/8Fd6crYhKLqs7C7TddkV/QFpq+AqiqTeCkhCBZselCbIDUJD/eIqqBQqHt/1bFl1VTzahLC7311dbdZj12PXtgFddVehgxyCHE3TlKNXMr8/2kZK09L0mmT6ej4efWgyk5l3pqV59T2f+YxKkiQJRERERIQIuQsgIiIiChUMRkREREQtGIyIiIiIWjAYEREREbVgMCIiIiJqwWBERERE1ILBiIiIiKgFgxERERFRCwYjIiIiohYMRkQhymAwIDExESqVCiqVCgaDoc1XSkoKzGYz3G63bDW63W4kJCSgoKBAthrCicPhQGJiIkaMGAGDwdBmGY8lUWiIkrsAIgrMarUCAFQqFbRarf/xmcxmM0aMGIHS0lLodLr+LhEulwuiKKKsrKzb27BYLNi8eXMvVhW6dDodKioqkJiY2G5ZbxxLIuo5BiOiMFZYWAiHwwGDwYC9e/dCrVb36/4FQUBPb7coimIvVRM+NBpNu+d641gSUc/xVBpRmNPpdHC73XA4HHKXEjSbzSbrqUAiorMxGBGRLERRRE5OjtxlDCgDsTtHFCyeSiMKc62dojPHGLndblgsFiQkJKCqqgqiKCIvLw9ardb/GovFAlEUYTKZkJqaCofDgdLSUpjNZuj1eoiiiMLCQiQkJPi7OpWVlTCbzdBqtRBFEWazGeXl5RAEARUVFf79n+u1NpsNxcXFAIDy8nL/QGRBENqMNzrX+3A6ncjJyYEoitDpdNi6dSuKioqgVqtRWlrabntnKioqQkVFhf/0Y0pKSrtxWufaf1cUFBSgsrISCQkJUKvVEASh3TqdHcszj+fIkSNRVVUFAB2+L1EUYbFYkJyc3Ga/Go0GxcXFSE1NhVqtPuf3v6v7Pdf3QKPRoLCwEG63G0VFRQCAsrKyTr83RLKSiCikAZC0Wm3AZZs3b5YASKWlpf7nKisrJbVaLVVUVHT6nCRJkiAIkslkkjZv3ux/rNfrJUmSAu5Tp9O124ZOp2u3bjCv1el0Ad9bMO+jdTut76OVWq1u95wkSZJer5dMJlOb53Jzc7t9HDui1Wql3NzcNs+VlpZKarXaf5zPfh+Bjl2gWgOtV11dLanV6jbvQ6fT+V9fUVEhVVdX+5d19v0PZr+t++noe5Cbm9vl7w2R3BiMiEIcAP+Hy5lfrR/ulZWVbdbX6XQBP3T1en27EKLT6SS1Wu3/sGz9b0VFhSQIQrttWK3WdqHAZDK1+bAM5rWdBaNg3oder5cABDwWZ3+QFxYWSgDaBARJah8Kgtl/ILm5uZJarQ64TKvVBtz22cdSkpqD1Nnht7q6WgIgWa3WNuu2BuUzdfR+Janj73+w+5Wkzr8HXf3eEIUCjjEiCgOtpx3O/LJarSgsLGxzaqZ1EHZycnK7baSkpKC8vDzgtltPJ7X+VxAEiKKIlJSUNoO69Xr9OU8j9eS13X0fGo0m4GmqQKetLBYL9Hp9uyv4tFotUlNTu7X/QIqKijBnzpyAywJdldYRQRCg0+navJfW2s8eM1RZWdnhMehofFGg73+w+wU6/x509DwH3lMo4hgjIgVp/cCurKz0j+c4U6AxHYE+pNVqNaxWK3JycpCSkgKgOThs3rz5nPMl9eS1vf0+zuZ2u+F2uwMGpjPnierO/s8kiqJ/wsaeEgQBpaWlAJrrF0XRH0xax/20SkxMbFdva/joKJR2dNyC2e+5thVMECSSG4MRkYK0fgClpKT4B9CeS0dzH+n1euj1ev+gXJvNhpSUlC5NJtnd17rdbv9A4d56H2dyuVwAgJEjR3a6Xnf235dsNhsKCwv9Xa1AHS8AmDNnDiwWCxwOh/84FxcXIzc3t8Ntd3bcurrfc22rv+fXIuoJBiMiBWntCvT0suwzr3Rr/dq8eTMMBgOsVmun4aYnr229Kqu33sfZWjtFlZWVna7X0/13dT9dUVRUBIvFgoqKioCdrjOp1WrMmTMHVqvV3+1JTk7uNBj1xn6JlIRjjIgUJjc3F4WFhQGXmc3mLm+no1uQtHZdeuO1Z3cS3G63v1vTW+/jbK2drEDOPA3V0/3n5uZ2uB+Xy9Xl8TUWiwVz5sxpF07OfH3r/dUcDgdSUlJQWFjoH4vWnVAU7H6JlITBiCgMBDNIdfPmzVCr1e0+tGw2W7sbl3a27aKionYdk9LSUv8A5c620dXXnj2QWRRFf1gK5n10FDQCPb9161YAzR/8Z3I6nW3GwgR7HM/W+nqbzdbm+dZj01HAPLtejUbT7ljabDb/jOdnEgQB+fn5cDgccDqdbcYFdaSj738w+wWC+x60Pk8UilSSxJvzEIWi1gn/nE4ngOZTU62DmLui9YO/dTxN6+uB5hCQn5/v/9DW6/VtTrm0frCq1eo2H2qCILSZ/K91PIvb7YZer0deXh5cLtc5X3t2nU6nEykpKdBqte1OtXX2PgLVYDabIQhCm+d1Ol2biQtbt9s6QLr1qqlAp/k6239XWCwW/2urqqqQmpoKi8XiP224detW/4SLZx9LrVYLt9uNnJwcuN1u/2D21vfROmFmXl5em4kqA3Wq9Hq9f1/n+v4D6PJ+XS5Xl74HwXxviOTEYEREpBBmsxkJCQntAk55ebk/4J09qzYRtcVgRESkAK23eeko+LjdbowYMQLV1dW8SoyoExxjREQ0AASaxJGI2mMwIiJSAJ1OB41G0+GVYgaDgTdtJeoCnkojIlIQm82G0tLSNrNuV1ZWwmAwdHnmcaKBjMGIiIiIqAVPpRERERG14C1BguDz+XDo0CEMHz4cKpVK7nKIiIioCyRJQm1tLUaPHo2IiM57QgxGQTh06BDGjBkjdxlERETUDQcOHMBFF13U6ToMRkEYPnw4gOYDGxcXJ3M1RERE1BUejwdjxozxf453hsEoCK2nz+Li4hiMiIiIwkxXhsFw8DURERFRCwYjIiIiohYMRkREREQtGIyIiIiIWjAYEREREbVgMCIiIiJqwWBERERE1ILzGCmM1NiIhooK1H/4EbyHDkOlViN2ymTETpmMiGHD5C6PiIgopIV1x6igoEDuEkKKVF+P2meehWfLo6j/+GN4jxxB4+ef48Qzz8KzaTO8x4/LXSIREVFIC+tgVFVV1elyURRRVFQEp9PZbpnT6URRURFEUeyr8vrd6bfeQv377yNi1AWI+vnFiPzxjxElCIj86U/R8O9/4+TLr0CSJLnLJCIiCllhHYw6Y7PZ4HA4oNFokJOTA4PB4F9msVjgcrkwZ84ciKIYMDiFG6muDnWOd6EaPKTdKTNVTAwiL7gQDTt3wnvggEwVEhERhT5FjzEymUwAAL1ej8TERJjNZqjVapjNZgiCAADQ6XSw2WzQarXtXl9fX4/6+nr/Y4/HA6D5JrJduRFdf2o6cACe/fugGqFBREudZ5JUKjQdOQLXvz7FIBnqIyIikkttbW2X1w2bYOR2u5Gfn9/mOYfDAYvF4n88cuRI5Obmwul0QqfTtVn3b3/7G8aOHQuTyeQPRa3UanXAfebn52P9+vXtnt+yZQtiYmK6+U76hq+2Fg1ffQXExkAVFR1gDQm+Gg9itm1DpKO03+sjIiKSS0NDQ5fXVUlhPOjEYrFg8+bN7Z5vHTd0dgBKSUmBw+FAZWVlm2U2mw16vb7ddgJ1jMaMGYMvvvgi5DpGUmMjPJsL0HToMKIuuqjdcq/LBZXPi7i8tYg8/3wZKiQiIpJHbW0tLr30UtTU1CAuLq7TdcOmYxQMQRBgs9nahB+LxQKr1YqioiKkpKSgoqLC3ylyuVwBtxMbG4vY2Nh2z48ZM+acB1YOp/V6nCgsgqqxEREaDVQqFQDAd+IEvPV1GHzLLRiemChzlURERP3LE2CISUcUGYwAQKvVoqioCG63G1VVVf7xRbm5uRAEAdOnT4cgCBAEIWDXKRwN0k2H99hR1L35V3i/+RqIjgGaGoGYGAyaOhVD75ond4lEREQhTZGn0vqKx+NBfHx8l1pxcpEkCU2Vlaj/17/gO3oMEXFxiEnUIvrSS6GKUmwOJiIi6lAwn99hHYzcbneHA6f7QjgEIyIiImormM/vsJ7HqD9DERERESlfWAcjIiIiot7EYERERETUgsGIiIiIqAWDEREREVELBiMiIiKiFgxGRERERC0YjIiIiIhaMBgRERGRYnm9Xjz++ONdXp/3iCAiIiJF2rNnD9LT0/HJJ590+TXsGBEREZGiSJKEZ599FhMmTMDBgwfx3HPPdfm17BgRERGRYnz33XeYP38+3nnnHcyZMwf33HMPfD5fl1/PYERERERhT5IkvPbaa1iyZAliYmLw7LPP4qqrrgIAnDhxosvb4ak0IiIiCmvHjx+HwWBAWloapkyZArvd7g9FwWLHiIiIiMLWn//8Z2RnZ6OhoQGPPPIIbrzxxh5tjx0jIiIiCjsejwfZ2dm47bbb8Mtf/hJ2u73HoQhgx4iIiIjCzAcffID09HRUVVXhgQcewB133AGVStUr22bHiIiIiMLC6dOnsXLlSkybNg2jRo2CzWbD7Nmzey0UAewYKZIkSfDu2wdfdTVUQ4YgKiEBqih+q4mIKHyVl5fDaDRCFEXcc889MBqNiIjo/f4OPy0VpnFPJU4VF6Pxq68gna6DKiYakT/7GYbMmonYK66QuzwiIqKgNDY24qGHHsLDDz+McePGoaSkBAkJCX22PwYjBWkS96L20UfhPXQYERdeiIgLfwSpvh5Nu/eg9qmngSYvYn99pdxlEhERdcl//vMfpKWl4bPPPkNOTg5MJhOio6P7dJ8cY6Qgp/70JzR9dwiRF1+MiPh4qCIjEdFyKk2qq8dJqxVSfb3cZRIREXXK5/Ph8ccfx6RJk+B2u7Ft2zYsXry4z0MRwGCkGN5jx9C4cxciR42CKsA518jRo+E9cBCNX34lQ3VERERd8+2332LatGlYtWoVDAYDiouLcemll/bb/nkqTSGk2hOQ6usRERcXcLkqNhbweeGr9fRzZUREROcmSRJefPFFrFixAsOGDcPzzz+Pyy+/vN/rYMdIIVRxw6GKjYV0+nTA5VJ9PRAZ2WFwIiIiksv333+P2267DdnZ2Zg+fTrsdrssoQhgx0gxIs87D9GTJqLuvfegUqvbnU7zHjqEqDFjEH3JJTJVSERE1J7NZsPChQsBAL/97W8xbdo0Wethx0hBhsycgajRo+HdvRu+mhpIXi98J0+iqbISqkGxGGIwNJ9SIyIikll1dTXuuusuGAwGTJo0CXa7XfZQBLBjpChRP/sZhq9ejVMlJWj88iv4jh+DKjoG0b+4GIPvmIXY5GS5SyQiIsI777yDrKws1NbWYuPGjbj11lt7dfbqnmAwUpjoBAFxay3w7t/fPPP10KGIEgSoIiPlLo2IiAa4kydPYvXq1XjuuecwZcoUbNiwARdeeKHcZbXBYKRAKpUKUT/9KfDTn8pdChEREQDg448/Rnp6Or777jv8z//8D1JTU0OmS3QmjjEiIiKiPlNfX4+1a9fi6quvxrBhw2Cz2XDnnXeGZCgC2DEiIiKiPvLZZ5/BaDTiv//9L5YuXYqsrCxEhvjQDnaMiIiIqFc1NTVh48aNSE5ORn19PV5//XUsWLAg5EMRwI4RERER9aLdu3cjPT0dn376KbKysrBo0SLExMTIXVaXMRgRERFRj0mShGeffRZr1qzByJEj8dJLL2HixIlylxU0BiMiIiLqkYMHDyIrKwsOhwOpqalYtWoVhgwZIndZ3cJgRERERN0iSRJeffVVLFmyBLGxsXjuuefw61//Wu6yeoSDr4mIiChox44dg16vh9FoxFVXXQW73R72oQhgx4iIiIiC9Kc//Qk5OTloaGjAo48+iuuvv17uknoNO0ZERETUJR6PB1lZWZgxYwZ+9atfYceOHYoKRQA7RkRERNQF7733HjIzM1FVVYUNGzZg5syZITt7dU+wY0REREQdOn36NFasWIHrrrsOo0aNwvbt2zFr1ixFhiIgBDtGoijCZrNBEASIogiTyQS1Wt3h+g6HA6IoQhAEAIBOpzvndoLdBxER0UBUVlYGo9GIvXv3Ys2aNUhLS0NEhLJ7KiEXjAwGAyoqKgA0B5icnBxYrdaA6zocDlitVhQWFkIURaSkpKCysvKc2wlmH0RERANNY2MjHnzwQWzcuBHjxo1DSUkJEhIS5C6rX4RUMBJFsc1jQRDgcDg6XN9sNvsDjiAIKC0tPed2gt0HERHRQPLVV18hLS0Nn3/+OcxmMxYsWIDo6Gi5y+o3IdUPczgc0Gg0bZ7TaDRwOp3t1hVFES6XC2q1Gk6nE2632386rbPtBLOP+vp6eDyeNl9ERERK5PV68eijj0Kr1aKmpgavvvoq7r777gEVioAQC0Zutzvg8y6Xq91zTqcTGo3GP1aoqKgINpvtnNsJZh/5+fmIj4/3f40ZM6ZL74OIiCic7N27F9OmTcOaNWswZ84cFBcXY/z48XKXJYuQOpXWkUBhxuVyQRRF6HQ6qNVqmEwmjBgxApIkBbWdzpbl5eVh1apV/scej4fhiIiIFEOSJDz//PNYuXIl4uLi8PzzzyM5OVnusmQVUsFIrVa369y0ni47myAIUKvV/mWt/3U6nZ1uJ5h9xMbGIjY2ttvvh4iIKFQdPnwYCxYswJtvvok77rgDa9aswbBhw+QuS3YhdSqt9VL7syUlJbV7rnU8UbDbCWYfRERESmS1WnHppZfi008/xZNPPon169czFLUIqWB0dtgRRRFJSUltukGtV5UJgoCkpCT/KbDWuYy0Wm2n2znXPoiIiJTK5XJh7ty5mDNnDhITE2G323HttdfKXVZICalTaUBzirVYLEhOTkZZWVmb+YXy8/ORnJyM3NzcNusmJiaioqLCf7n+ubbT2TIiIiIlevvtt5GVlYWTJ09i06ZNuPnmmxU7e3VPqKTORitTGx6PB/Hx8aipqUFcXJzc5RAREZ3TiRMnsHr1ahQWFuLKK6/Ehg0bcMEFF8hdVr86ceIEpkyZ0qXP75DrGBEREVHv+Oijj2A0GvH9999j3bp1mDNnDrtE5xBSY4yIiIio5+rr65Gbm4urr74a8fHxsNlsSE1NZSjqAnaMiIiIFGTXrl1IS0vDN998g+XLlyMzMxORkZFylxU22DEiIiJSgKamJjz88MO4/PLL0djYiNdffx3Z2dkMRUFix4iIiCjMffPNN0hPT0dZWRmys7MH5D3OeguDERERUZjy+Xx45plnkJubi/PPPx8vvfQSJk6cKHdZYY3BiIiIKAwdOHAAWVlZ+Nvf/oY777wTK1euxJAhQ+QuK+wxGBEREYURSZKwbds2LF26FIMGDfLPT0S9g4OviYiIwsSxY8dwxx13ID09HVdffTW2b9/OUNTL2DEiIiIKA3/84x+Rk5ODxsZGPPbYY0hJSZG7JEVix4iIiCiE1dTUIDMzEzNnzsT48eOxY8cOhqI+xI4RERFRiHr33XeRmZkJl8uFBx98EDNmzODs1X2MHSMiIqIQc+rUKSxbtgzTp0/HhRdeCLvdjpkzZzIU9QN2jIiIiELIp59+CqPRiH379sFisWDevHmIiGAfo7/wSBMREYWAhoYG3HfffbjyyisRHR2NkpISpKWlMRT1M3aMiIiIZPbll18iLS0NX3zxBcxmMxYsWMBbesiEMZSIiEgmXq8XW7ZsQWJiImpra7Ft2zbe50xm7BgRERHJQBRFZGRk4KOPPkJ6ejqWLl2K2NhYucsa8BiMiIiI+pEkSdi6dStWrVoFtVqNF154AUlJSXKXRS0YjIiIiPrJ4cOHMX/+fPz1r3/F7NmzsWbNGgwdOlTusugMDEZERET9oLi4GIsWLUJkZCSefvppTJ06Ve6SKAAOviYiIupDLpcLd955J+68804kJSXBbrczFIUwdoyIiIj6yFtvvYX58+fj9OnT2Lx5M2666SbOXh3i2DEiIiLqZSdOnIDJZMLNN9+MhIQEbN++HTfffDNDURhgx4iIiKgXffjhh0hPT8eRI0dw3333wWAwMBCFEXaMiIiIekFdXR3WrFmDqVOnQq1Ww2azYc6cOQxFYYYdIyIioh7auXMn0tLSsGfPHqxYsQIZGRmIjIyUuyzqBnaMiIiIuqmpqQkPPvggLr/8cni9Xrz++uuYP38+Q1EYY8eIiIioG77++msYjUZUVFQgOzub9zhTCHaMiIiIguDz+fDkk09i0qRJOHr0KF5++WUsW7aMoUgh2DEiIiLqov379yMzMxPvvfce5s6di5UrV2Lw4MFyl0W9iMGIiIjoHCRJ8neGhgwZgqKiIkyZMkXusqgP8FQaERFRJ44ePYpZs2YhMzMT1157LbZv385QpGDsGBEREXVgx44dMJlM8Hq9eOKJJzB9+nS5S6I+FnLBSBRF2Gw2CIIAURRhMpmgVqsDrut0OgEAWq0WoijC7XZDq9WeczvB7IOIiAYet9uNZcuW4ZVXXsG0adPwm9/8Buedd57cZVE/CLlgZDAYUFFRAaA5wOTk5MBqtQZct7CwEEVFRQAAnU7XZr3OthPMPoiIaGBxOBzIzMxETU0NHnroIdx+++2cvXoACalgJIpim8eCIMDhcHS4fmJiIqqrqwGgTcens+0Euw8iIhoYTp06BYvFgqeeegpXXHEFXnzxRfzoRz+SuyzqZyE1+NrhcECj0bR5TqPR+E+ZBaJWq9udButsO93ZBxERKdu//vUvTJw4EVu3bsXatWtRVFTEUDRAhVTHyO12B3ze5XJ1uL7NZgMAlJWVwWw2QxCETrcTzD7q6+tRX1/vf+zxeDounoiIwk5DQwPWr1+PTZs24ZJLLkFJSQkEQZC7LJJRSAWjjnQUZs4cNC0IAlJSUlBZWRn0djpalp+fj/Xr1wdRKRERhYsvvvgCaWlp+PLLL7Fo0SJkZ2cjKiosPhapD4XUqTS1Wt2uc+NyuTq8YuzM8UKtV5iJotjpdoLZR15eHmpqavxfBw4c6N4bIyKikOH1elFQUIDExEScPHkSr776KsxmM0MRAQixYKTT6QI+n5SU1O45p9MZcD4JjUbT6XaC2UdsbCzi4uLafBERUfiqrKzENddcg7Vr12LevHl44403cMkll8hdFoWQkIrHZ5/XFUURSUlJ/m6O0+mEWq2GIAgQBAGbN2/2r+twOKDX6wMOxj5zO50tIyIiZZIkCUVFRbjnnnswYsQIvPjii0hMTJS7LApBIRWMAMBqtcJisSA5ORllZWVt5hfKz89HcnIycnNzoVarkZSUhIKCAqjValRWVrZZt7PtdLaMiIiU5dChQ5g/fz7efvtt6PV6rF69GkOHDpW7LApRKkmSJLmLCBcejwfx8fGoqanhaTUiojDwxhtvYNGiRYiKisIDDzyAqVOnyl0SyeDEiROYMmVKlz6/Q2qMERERUW+oqqrCnDlzMHfuXFxxxRWw2+0MRdQlIXcqjYiIqCfefPNNzJ8/H3V1dSgoKMBNN90kd0kURtgxIiIiRaitrUVOTg5uueUWXHzxxbDb7QxFFDR2jIiIKOz9/e9/R0ZGBo4ePYrf/OY30Ov1vPErdQs7RkREFLbq6upwzz334Nprr4VGo4HNZoPBYGAoom5jx4iIiMKS0+lEWloaKisrsWrVKhiNRkRGRspdFoU5doyIiCisNDU1YcOGDbjiiisANF+Sn5mZyVBEvYIdIyIiChv//e9/YTQasXPnTmRnZ2PhwoWIjo6WuyxSEHaMiIgo5Pl8Pvzv//4vJk2ahOPHj+Pll1/G0qVLGYqo17FjREREIW3fvn3IzMzE+++/j7vuugvLly/H4MGD5S6LFIrBiIiIQpIkSXjppZewbNkyDBs2DFu3bsXkyZPlLosUjqfSiIgo5Bw5cgQzZsxAVlYWpk2bhu3btzMUUb9gx4iIiEKK3W6H2WyGz+fDE088genTp8tdEg0g7BgREVFIcLvdMBqNmD17NiZMmAC73c5QRP2OHSMiIpJdaWkpsrKyUFNTg4cffhi33XYbZ68mWbBjREREsjl58iQWL16M66+/HhdddBHsdjtuv/12hiKSDTtGREQki08++QRGoxHfffcd8vLycOeddyIign+vk7z4E0hERP2qoaEB9957L6666ioMHToUVqsV8+bNYyiikMCOERER9ZvPP/8cRqMRX331FRYvXoz58+cjKoofRRQ6GM+JiKjPeb1ebNq0CcnJyTh9+jRee+01mEwmhiIKOfyJJCKiPrVnzx5kZGTgk08+QWZmJpYsWYKYmBi5yyIKiMGIiIj6hCRJKCwsxD333IORI0fi97//PbRardxlEXWKwYiIiHrdd999h/nz5+Odd96BwWDA6tWrMWTIELnLIjonBiMiIuo1kiTh9ddfx+LFixEdHY1nnnkGV199tdxlEXUZB18TEVGvOH78OObMmYO77roLU6ZMwY4dOxiKKOywY0RERD325z//GdnZ2WhoaMAjjzyCG2+8Ue6SiLqFHSMiIuo2j8eD7Oxs3HbbbRg3bhzsdjtDEYU1doyIiKhbPvjgA2RkZOD48eO4//77MXv2bN7jjMIeO0ZERBSUuro6rFq1CtOmTcPIkSNhs9mg1+sZikgR2DEiIqIuKy8vh9FohCiKuOeee5CWlobIyEi5yyLqNewYERHROTU2NuKBBx7A5MmTERERgeLiYmRkZDAUkeKwY0RERJ36z3/+A6PRiF27diEnJwcmkwnR0dFyl0XUJ9gxIiKigHw+H5544glotVq4XC5s27bNP3EjkVKxY0RERO18++23yMzMxAcffIC0tDQsX74cgwYNkrssoj7HYERERH6SJOHFF1/EihUrMGzYMPzud7/DFVdcIXdZRP0m5IKRKIqw2WwQBAGiKMJkMkGtVp/zdRaLBXl5ef51O9tOd/dBRKRk33//PXJycvDnP/8ZM2fORG5uLoYPHy53WUT9SiVJkiR3EWdKTExERUUFgOYAY7FYYLVaO32N0+lEYmIiqqur/QGns+10Zx9A8wyv8fHxqKmpQVxcXHffIhFRyNm+fTtMJhNUKhXuv/9+TJs2Te6SiHrNiRMnMGXKlC59fvf54Gu73d7ldUVRbPNYEAQ4HI4uvU4QhC5tp7v7ICJSourqatx1113Q6/XQarWw2+0MRTSg9XkwKiws7PK6DocDGo2mzXMajQZOp7PD17TOuNrV7QSzj/r6eng8njZfRERK8c477+DSSy/F//3f/2Hjxo147LHH2v1+JBpoejzGaMuWLSguLg64zO12t+vQdMbtdgd83uVydbh+oLFBnW0nmH3k5+dj/fr1AdcnIgpXJ0+exJo1a/Dss89i8uTJePDBB3HhhRfKXRZRSOhRMFq7di2KioqQlJTU5lRWq6qqqg5DTTA6CjMlJSUwmUw93k5Hy/Ly8rBq1Sr/Y4/HgzFjxnR5f0REoebjjz9Geno6vvvuO9x7771ITU1FRASntCNq1aNgJIriOYPPnDlzurw9tVrdbnsulytgV8jhcHS47c62E8w+YmNjERsb2+X6iYhCVX19Pe6//3488sgjuOyyy2C1WvGzn/1M7rKIQk6P/kxITk4+5zqbN2/u8vZ0Ol3A55OSkgI+X1JSgqKiIhQVFUEUReTn58PpdHa6nWD3QUQU7j7//HMkJyfjsccew9KlS/H73/+eoYioA30+j9HevXsxduzYLq179uk4URSRlJTk7+Y4nU6o1WoIgtAu4JjNZpjN5oCn9M7cztmdobP3QUSkFF6vFwUFBbj//vsxduxYvP766xg3bpzcZRGFtB51jEwmE7Zs2YJvv/22w3WCuSoNAKxWKywWC2w2GwoLC9vML5Sfnw+bzdZmfbfbjYKCAgDN3anWq8s6205ny4iIlGD37t246qqrsG7dOhiNRoYioi7q8gSPN9xwQ8DxRJIkYefOnVCr1dBoNG06L61XpXm93l4rWE6c4JGIQp0kSXj22WexZs0ajBw5Eg8//DAmTZokd1lEsgpmgscun0orKyvr8Oqzjsbn9NZVaUREdG4HDx7E/PnzUVpaitTUVKxatQpDhgyRuyyisNLlYCQIAt55552gdxDMVWlERBQ8SZLw2muvYfHixYiNjcVzzz2HX//613KXRRSWujzGaOvWrd3aQTBXpRERUXCOHTsGvV6PtLQ0/PrXv4bdbmcoIuqBLneMgjlH/e2330Kj0SAuLq7LV6QREVFw/vSnP2HBggVobGzEli1bcMMNN8hdElHY69FVaXl5ee2eq6mpQWVlJUpLS7Flyxa8++67PdkFERGdxePxYP78+ZgxYwYuueQS7Nixg6GIqJf0aB6jysrKds/Fx8dj+vTp/sdbtmzBdddd15PdEBFRi/fffx8ZGRmoqqrChg0bMHPmTKhUKrnLIlKMHnWMuvKPsbS0tCe7ICIiAKdPn8aKFSswbdo0jBo1CjabDbNmzWIoIuplQXWM1q5dC1EUUVNTAwAoLy/vtH1bXl4e1E1eiYiovfLycqSlpWHv3r1Ys2YN0tLSeONXoj4SVDDatGkTAMBms8FkMkGlUqGj+SHVajU2bdqEnJycnldJRDQANTY24qGHHsLDDz+McePGoaSkBAkJCXKXRaRo3RpjpNfrodVqsXbtWpSUlPR2TUREA95XX30Fo9GIzz77DCaTCTk5OYiOjpa7LCLF63YvVhAEpKam9mYtREQDns/nw2OPPQatVgu3241XX30VixYtYigi6ic9uipt9uzZHS575JFHoFKpYDKZeF8xIqIu+Pbbb5GRkYF//OMfSEtLw7JlyzBo0CC5yyIaUHoUjDqzZs0aAM1zHeXn5/fVboiIwp4kSXjhhRewYsUKxMXF4fnnn0dycrLcZRENSL0SjOx2O0RRRFVVVZvn3W43RFHsjV0QESnS999/j+zsbLz55puYNWsWcnNzMWzYMLnLIhqwehyMkpKS4HQ6/Y/VajWA5lCUkpKC5557rqe7ICJSJKvVioULFyIiIgJPPvkkrr32WrlLIhrwejQRxtq1a6HT6VBdXQ2fzwer1QqXywWXywWfz+e/pJ+IiH5QXV2NefPmYc6cOUhMTITdbmcoIgoRPZ4hbNOmTYiPjwfQfKXarl27/Mtmz54Nm83W010QESnG22+/jfHjx+Mvf/kL8vPz8eijj2LEiBFyl0VELXoUjM4777w2jwVBQHFxcY8KIiJSohMnTmDhwoW48cYb8bOf/Qzbt2/Hrbfeyq46UYjpUTA6fvw4AODdd9/Ft99+i/j4eFRUVGDfvn3+dXivNCIa6D766CNMmDABL7/8MtatW4fCwkJceOGFcpdFRAH0KBiZzWYsXLgQOp0OZrMZAGAymaDVanH33Xfjhhtu8A/GJiIaaOrr62GxWDB16lQMHz4cNpsNqamp7BIRhbAeXZU2duxYPPfcczAYDBAEAUDz7UJcLhfWrl2LkSNHoqKiolcKJSIKJ7t27YLRaMTXX3+NpUuXIisrC5GRkXKXRUTnoJI6ugsstePxeBAfH4+amhrO5k1EATU1NWHz5s1Yv349xo4di40bN2LcuHFyl0U0oJ04cQJTpkzp0ud3n818TUQ00OzevRtGoxFlZWWYP38+7r77bsTExMhdFhEFoceX6+/atQvXX389Ro4cid/97nf+5xcuXIh33323p5snIgp5Pp8PTz/9NCZMmIDDhw/jpZdewvLlyxmKiMJQj4LRzp07cd1110GtVmPTpk1tlj333HOorq5uM68REZHSHDhwANdffz2WLFmC22+/HVarFRMnTpS7LCLqph6dStu0aRMqKiowduxYAGjTMQKaJ3jcsmULf0kQkeJIkoRt27Zh6dKlGDRoEAoLC3HllVfKXRYR9VCPOkZjx471hyIiooHi2LFjmD17NtLT03H11Vdj+/btDEVECtGjjtHZM18HusCtqqqqJ7sgIgopf/zjH5GTk4PGxkY89thjSElJkbskIupFPeoY7dmzB5999pn/8dmTlm3ZsqUnmyciChk1NTXIzMzEzJkzMX78eOzYsYOhiEiBejzGSBAEpKSkIDk5GZWVldBoNBBFEYWFhVCr1SgrK+utWomIZPHee+8hIyMDLpcLGzZswMyZMzl7NZFC9ahjpFarUV5eDpfLhdzcXBQWFkKv1yM3NxezZ89mKCKisHb69GksX74c1113HS688EJs374ds2bNYigiUrAeT/AoCAJKS0tRU1OD8vJyaDQaTJo0qTdqIyKSzaeffgqj0Yh9+/YhNzcXd911FyIiejz1GxGFuF6b+To+Ph7Tp0/vrc0REcmisbERGzZsQH5+Pn75y1+ipKTEfy9IIlK+HgejXbt2weVyQaPRQBAE3kOMiMLWl19+CaPRiH//+98wm81YsGABoqOj5S6LiPpRt4KRx+NBTk4ObDZbu2UpKSkoLCzET3/6024VJIoibDYbBEGAKIowmUxQq9UB13U4HAAAt9uNsrIypKamQqvVnnM7weyDiJTP6/Xi8ccfx7p16/DjH/8Y27Ztw/jx4+Uui4hkEHQw2rJlCywWCyRJgiAIUKvV0Gg0cLlccLvdeOeddyAIAiwWCzZu3Bh0QQaDARUVFQCaA0xOTg6sVmuH6/7tb3+DTqeDy+WCwWBAZWXlObcTzD6ISNlEUURGRgY++ugjGI1G/0zWRDQwBRWMtm7dio0bN2LTpk0wmUyIj49vt05NTQ2Ki4uxdu1ajBw5Evfcc0+Xty+KYpvHgiD4u0KBWK1Wf4cIQJuOUEfbCXYfRKRMkiThd7/7HVauXAm1Wo3nn38eycnJcpdFRDLr8iUWO3fuxObNm1FRUYE1a9YEDEVA8yBsk8kEURTx8MMPt5kA8lwcDgc0Gk2b5zQaDZxOZ8D1dTqd//+tVivMZvM5txPsPohIeQ4fPoxbbrkFJpMJN954I7Zv385QREQAgugYbdq0CYWFhV2+N5parUZJSQk2btyI4uLiLr3G7XYHfN7lcnX4GqfTieLiYqSkpMBkMp1zO8Hso76+HvX19f7HHo+nwzqIKDyUlJTg7rvvRkREBJ5++mlMnTpV7pKIKIR0uWNUWVkZ9OX4Op0O1dXVQRd1to7CDABotVrk5eWhsrIy4GDwrm4n0LL8/HzEx8f7v8aMGdPFioko1LhcLtx5551ITU1FUlIS7HY7QxERtdPlYNTdeTyCudpLrVa369y4XK5zbkOtVsNgMMBgMMDtdne6nWD2kZeXh5qaGv/XgQMHuvxeiCh0vPXWWxg/fjzeeustbN68GVu2bMGIESPkLouIQlCXg1F3p8A/ezxPZ84cM3SmpKSkds85HI42v9hag5soip1uJ5h9xMbGIi4urs0XEYWPEydOwGw24+abb0ZCQgLsdjtuvvlm3tKDiDrU5TFGkiR1awfB/AI6uysliiKSkpL83Ryn0wm1Wg1BEKDRaNqEnNZlZ16lFmg7Z3eGzt4HESnDhx9+iPT0dHz//fe47777YDAYGIiI6Jy6HIy6+wuls4HTgVitVlgsFiQnJ6OsrKzN/EL5+flITk5Gbm4utFotUlNTUVRUBAAoLS31z010ru10toyIwltdXR3uu+8+PProo5gwYQKeeuop/OQnP5G7LCIKEyqpi62giIgIJCUlBd05cjqd8Hq93Sou1Hg8HsTHx6Ompoan1YhC0K5du5CWlobdu3dj8eLFyMjIQGRkpNxlEZHMTpw4gSlTpnTp8zuoCR7Ly8uDLoatayLqa01NTdi0aRPWr1+PhIQEvP766/jFL34hd1lEFIa6HIy0Wm23glGgQc1ERL3l66+/Rnp6OsrLy5GdnY27776bN34lom7r88v1u/s6IqLO+Hw+PPnkk5g0aRKOHDmCl156CcuWLWMoIqIe6fIYI+IYI6JQceDAAWRmZuLdd9/F3LlzsWLFCgwZMkTusogoRPXZGCMiIjlJkoSXX34Zy5Ytw5AhQ1BUVIQpU6bIXRYRKUiXT6UREcnp6NGjmDVrFjIzM3HNNddg+/btDEVE1OvYMSKikPeHP/wBOTk58Hq9ePzxxzucwZ6IqKfYMSKikFVTU4OMjAzMmjULl112Gex2O0MREfUpdoyIKCT97W9/Q2ZmJqqrq/Hggw9ixowZnBeNiPocO0ZEFFJOnTqFpUuXQqfTYfTo0bDb7Zg5cyZDERH1C3aMiChk/Otf/4LRaMT+/fuxdu1azJ07FxER/PuNiPoPf+MQkewaGhqwbt06XHnllYiNjUVJSQnuuusuhiIi6nfsGBGRrL744gsYjUZ88cUXuPvuu7FgwQJERfFXExHJg3+OEZEsvF4vCgoKkJiYiNraWrz66qtYuHAhQxERyYq/gYio34miiPT0dHz88cfIyMjAkiVLEBsbK3dZREQMRkTUfyRJwtatW7Fq1SqMGDECL7zwApKSkuQui4jIj8GIiPrFoUOHkJ2djb/+9a+YPXs21qxZg6FDh8pdFhFRGwxGRNTn3njjDSxatAiRkZF4+umnMXXqVLlLIiIKiIOviajPVFVVITU1FXPnzsUVV1yBHTt2MBQRUUhjx4iI+sSbb76J+fPno66uDgUFBbjpppvkLomI6JzYMSKiXlVbW4ucnBzccsst+PnPfw673c5QRERhgx0jIuo1//jHP5Ceno6jR4/iN7/5DfR6Pe9xRkRhhR0jIuqxuro6rF69Gtdccw00Gg1sNhsMBgNDERGFHXaMiKhHnE4njEYj9uzZg5UrVyI9PR2RkZFyl0VE1C3sGBFRtzQ1NWHDhg244oor4PP58MYbbyArK4uhiIjCGjtGRBS0r7/+GkajERUVFViwYAEWLlyI6OhoucsiIuoxdoyIqMt8Ph9++9vfYuLEiTh27BheeeUVLF26lKGIiBSDHSMi6pL9+/cjIyMD77//PubNm4cVK1Zg8ODBcpdFRNSrGIyIqFOSJOGll17C8uXLMWTIEGzduhWTJ0+Wuywioj7BU2lE1KGjR49i5syZyMrKwrXXXgu73c5QRESKxo4REQVkt9thMpkgSRKeeOIJTJ8+Xe6SiIj6HDtGRNSG2+2G0WjE7NmzMWHCBNjtdoYiIhow2DEiIj+Hw4HMzEzU1NTg4Ycfxm233cbZq4loQGHHiIhw6tQpLFmyBCkpKbjoootgt9tx++23MxQR0YDDjhHRAPfPf/4TRqMRBw8exNq1azF37lxERPBvJiIamEIuGImiCJvNBkEQIIoiTCYT1Gp1wHWdTiccDgcAoKysDFu3bvWv29l2gtkHkVI1NDTggQcewObNmzF+/HiUlJRg7NixcpdFRCSrkAtGBoMBFRUVAJoDTE5ODqxWa8B1HQ4HcnNzAQAFBQWYPn26/7WdbSeYfRAp0b///W+kpaXhq6++wuLFizF//nxERYXcrwMion4XUv1yURTbPBYEwd8ROpvT6UR+fr7/sV6vh9PphCiKnW4nmH0QKY3X68WmTZuQlJSE06dP47XXXoPJZGIoIiJqEVLByOFwQKPRtHlOo9HA6XS2W1er1WLr1q3+x263279+Z9sJZh/19fXweDxtvojCVWVlJaZOnYp7770Xd911F15//XX86le/krssIqKQElLBqDXcnM3lcgV8Xq/X+/+/uLgYOp0OarW60+0Es4/8/HzEx8f7v8aMGdNp/aHEe+QIGv/7XzTt3w9JkuQuh2QkSRKee+45TJgwAQcOHMCLL76IVatWITY2Vu7SiIhCTlj0zzsKM2cut9ls/nFD3dlOoGV5eXlYtWqV/7HH4wn5cNR04ABObbej0bkTvrrTUEVHI+riX2DIrBmImTBB7vKon3333XfIzs7G22+/DYPBgNWrV2PIkCFyl0VEFLJCKhip1ep2nRuXy3XOK8YsFgtKS0v963W2nWD2ERsbG1Z/VTcdPAhPwSNo2rcfkaNGIVKjgVRXh8bPdqF2714MW7IYsclJcpdJ/UCSJLzxxhtYtGgRoqOj8cwzz+Dqq6+WuywiopAXUqfSdDpdwOeTkjr+MC8oKIDFYoEgCHC73XC73Z1upzv7CBen//hHNO3bh6hf/AIRGg1UMTGIiItD5M8vhu/kSZwqLobU2Ch3mdTHjh8/jjlz5mDevHmYMmUKduzYwVBERNRFIRWMBEFo81gURSQlJfm7Oa1XnbWy2WzQarX+UFRSUgK1Wt3pds61j3DldbnQUFaBiPNHQXXW5HwqlQqRo0ejad8+NH71lUwVUn/4y1/+gksvvRQOhwOPPPIICgoKEB8fL3dZRERhI6ROpQGA1WqFxWJBcnIyysrK2swvlJ+fj+TkZOTm5kIURRgMhjavVavVMJlM59xOZ8vCleSugVRXh4jzzw+4XDVoENDUBN85xmtRePJ4PFi5ciVeeOEFXH311Vi/fj3O7+BngYiIOqaSeMlSl3k8HsTHx6OmpgZxcXFyl9OG98gRuNdYgNhYRIwY0W651NAA78EDiM/LQ0yiVoYKqa988MEHyMjIwPHjx7F69WrMnj2b9zgjIjrDiRMnMGXKlC59fofUqTTqvohRoxA94f/Bd/RIwMvzvYcPI3L0aESPv0SG6qgv1NXVYdWqVZg2bRpGjhwJq9UKvV7PUERE1AMhdyqNukelUmHw7beh8b9fw7t7NyJGj4Zq6FCgoQHew4ehiorEkDtmNZ9So7BXUVEBo9GIPXv24J577kFaWhoiIyPlLouIKOyxY6Qg0RdfjLiVKxB92aWQXC54d++G99B3iPrxaAxbaEbsNdfIXSL1UGNjI9avX4/JkycDAEpKSpCRkcFQRETUS9gxUpjoS36F+PUPoGn3bviqq6EaMgTR48ZBFUbzMVFg//nPf2A0GrFr1y4sWLAAZrMZ0dHRcpdFRKQoDEYKpIqIQPS4cXKXQb3E5/Pht7/9LfLy8vCjH/0Ir7zyCi677DK5yyIiUiQGI6IQtm/fPmRkZOCDDz5AWloali1bhsGDB8tdFhGRYjEYEYUgSZLw+9//HsuXL8ewYcPwu9/9DldccYXcZRERKR4HXxOFmCNHjuD222/H/Pnzcd1112H79u0MRURE/YQdI6IQsn37dpjNZgDA//7v/+K6666TuSIiooGFHSOFkurr4XW54Dt1Su5SqAvcbjfS0tKg1+sxceJE2O12hiIiIhmwY6QwPrcbp995B/Uf/B1S7QkgNgYxl1+OwTfegKgxY+QujwJ45513kJWVhdraWmzcuBG33norZ68mIpIJO0YK4nVVw/PIFpzc9hp87hpg8GBIjU2o+78/w5O/CY179shdIp3h5MmTWLRoEW644Qb85Cc/gd1ux2233cZQREQkI3aMFOT0n/6Ihn//G1EJP4fqjIn/pPPOg3fPHpx8+RXEP3A/VBHMw3L75JNPYDQa8d133+Hee+9FamoqIvh9ISKSHX8TK4TP40H9hx8hYoSmTSgCmid8jBg9Gk3ffIOmr7+WqUICgPr6euTl5eGqq67CsGHDYLVaMXfuXIYiIqIQwY6RQviOV0E6cRIR558fcHnEsGFoOnQI3mPHEP2rX/VzdQQAn3/+OdLS0vDf//4XS5YsQVZWFqKi+E+QiCiU8M9UpYiNBaKiIDU2BlwsNTUBESreM00GXq8X+fn5SEpKQl1dHV577TXk5OQwFBERhSD+ZlaIyNE/QvTFP0f9rs8QERfXbrnv6FFEnH8+oi+5RIbqBq49e/YgPT0d//znP5GVlYXFixcjJiZG7rKIiKgD7BgphEqlwuBbb0XEsKFo2rfP3zmSfD54jxyBdPoUBt14IyKGD5e50oFBkiQ888wzmDBhAg4ePIiXXnoJK1euZCgiIgpxDEYKEpOoxTCTCZEjNfDu+xaNu7+Bd88eqCIiMOTOOzHk9tvkLnFAOHjwIG644QYsXrwYt956K2w2GyZNmiR3WURE1AU8laYwg66+CjGTJqJh5y74qqsRMWQIoidOQOR558ldmuJJkoTXXnsNS5YsQUxMDJ599llcddVVcpdFRERBYDBSoIhhwxD9i4vhq66GauhQRIwcKXdJinf8+HEsXLgQ27dvx80334x7770X8fHxcpdFRERBYjBSmKYDB3Bqux2NO3fCd/o0VNExiB73CwyecTtiJkyQuzxF+r//+z8sWLAA9fX12LJlC2644Qa5SyIiom7iGCMFaTp4EJ6CR1D37ntAdAwifzQaquHD0bBzJ2offwL15RVyl6goHo8H8+fPx+23345f/vKX2LFjB0MREVGYY8dIQU7/8Y9o2rcPURf/AqrISACAKjYWquHD4d27F6feeAMxE/5fu5mxKXjvv/8+MjIyUFVVhfXr12PWrFm8xxkRkQKwY6QQXpcLDWUViDh/lD8UtVKpVIgcPRpN+/ah8auvZKpQGU6fPo0VK1Zg2rRpOP/882Gz2XDHHXcwFBERKQQ7RgohuWsg1dV1eEsQ1aBBQFMTfG53/xamIOXl5UhLS8PevXuxevVqGI1G3uOMiEhh+FtdIVRDhwAxMZDq6gIulxoagMhIRAzjBI/BamxsxP3334/JkycjMjISJSUlyMjIYCgiIlIgdowUImLUKET/v8vQ8I9/APHxwOnTQH09EBXVPMbo8GFEjh6N6PG8JUgwvvrqKxiNRnz22WcwmUzIyclBNMdoEREpFoORQqhUKgyZcTsanE40fPA+0ORtXQCoVIj80YUYcoe5+ZQanZPP58MTTzyBe++9Fz/+8Y+xbds2XHrppXKXRUREfYzBSEFUMTFQRUQCjU2QfD4AEgBV88BgCcDgwTJXGB6+/fZbZGRk4O9//zvS0tKwfPlyDGKgJCIaEBiMFOT0H/8En8eDWN10SLW1P5xKi4+Hd99+nC6xInbSJKh4I9OAJEnCCy+8gBUrVmD48OF4/vnncfnll8tdFhER9SOOHlUIb1UVGpxORJx/PlQRkVANj4MqXt38FRnVcrn+fl6u34Hvv/8et912GxYsWICUlBTY7XaGIiKiAYgdI4WQajyQ6uqhih+Epr174T10CFJjI1QREYi44AJEXnQR4G2Cr6ZG7lJDjs1mg9lshkqlwpNPPolrr71W7pKIiEgm7BgphGr4MCBChcbPP0fTN99Aqq+HKiICkteLpm+/RUNFBaSGBl6uf4bq6mrMmzcPBoMBiYmJsNvtDEVERANcyHWMRFGEzWaDIAgQRREmkwlqtbrD9Z1OJ3JyclBR0fY+YJ1tJ9h9hIPI888HoqLgrapCxKhRiGi9JQgAafBg+I4ehS82FlG/+qW8hYaIt99+G/Pnz0dtbS3y8/Nxyy23cPZqIiIKvWBkMBj8IUcUReTk5MBqtQZctzXcOJ3OoLYTzD7Cha+2FlJTEzB4MHzHjsHX2Aj4fM2X60dHA7GxUA0aBO/+/Yj45cANRydPnsTq1avx3HPPYcqUKdiwYQMuvPBCucsiIqIQEVLBSBTFNo8FQYDD4ehwfb1eH/R2gt1HuPAdOw6psQkqnw/SyZOAJP2wsKGhOSAB8B49iugBGow++ugjpKen4/Dhw/if//kfpKamsktERERthNQYI4fDAY1G0+Y5jUYTsCPU3e301j5CTmwMvPv3Q6qpASIi2n5FRgKnTqHp22+bu0cDTH19PSwWC6ZOnYrhw4fDarXizjvvZCgiIqJ2Qqpj5O7gBqcul6vXthPMPurr61FfX+9/7PF4gqqjP6mGDIPP5Wo+fdYaiFSq5s5R61ddHaTTge+lplSfffYZ0tLS8PXXX2Pp0qXIyspCZMv4KyIiorOFVMeoIx2Fmd7cTqBl+fn5iI+P93+NGTOmV+roCw1ffgE0NjY/kCT/qbPWW4IAAHw+NH7+uTwF9rOmpiZs3LgRycnJaGhowOuvv44FCxYwFBERUadCKhip1ep2nRuXyxX0FWOdbSeYfeTl5aGmpsb/deDAgaDq6E++6urmQBQV1dwt8np/+JIkoGW2a++xYzJX2vd2796Nq666Cvfddx/S09Px+uuvY9y4cXKXRUREYSCkgpFOpwv4fFJSUq9tJ5h9xMbGIi4urs1XqIoe+7PmUCRJzcHoTBER/uejf/5zWerrDz6fD08//TQmTpyIQ4cO4aWXXsKKFSsQw1ugEBFRF4VUMBIEoc1jURSRlJTk7+Y4nc52V5W1OvNUWGfbOdc+wlXURRch4rzzgKam5lNqZw6+9nqbr0wbNAgxycGFzHBx8OBB3HDDDViyZAluu+022Gw2TJw4Ue6yiIgozITU4GsAsFqtsFgsSE5ORllZWZv5hfLz85GcnIzc3FwAzVeflZaWtlnWegl/Z9vpbFm4ilCrEXXxz9Fw+HDzAGyvt+3g64gIRJ53HqISEuQutVdJkoRt27Zh6dKlGDRoEAoLC3HllVfKXRYREYUplSSdOeENdcbj8SA+Ph41NTUhd1rN63KhevFSNO7ZA5/HA9TX/3CFWkwMVMOHI2rUKKg3bUTMhAlyl9srjh07BrPZjB07duCWW25BXl4e4uPj5S6LiIhCzIkTJzBlypQufX6HXMeIukdy1wCShMgfj4avxg14mwAJ/nAUed5IIDYGvl66wk9uf/rTn7BgwQI0Njbi0UcfxfXXXy93SUREpAAhNcaIuk81dAi8tbVo/PcXQH0DEBMLDBoExMYCPh+avtkN39FjUA0dJnepPeLxeJCVlYUZM2Zg/Pjx2LFjB0MRERH1GnaMFEJ1/vnwVR1vPoUWG/vDZfpA88zXjY3wHj+OyJ/+RN5Ce+C9995DRkYGXC4XNmzYgJkzZ3L2aiIi6lXsGCmEd88e+GpPND+oq2u+Mq2pqTkg1dc3hySvF/Xvvidvod1w+vRpLF++HNdddx0uuOACbN++HbNmzWIoIiKiXseOkUJ4D3/fHIhaw0Jrt+jMsfVNTWg6eqT/i+uBsrIyGI1G7N27F7m5ubjrrrsQcfY8TURERL2EnzAKoRo+rLlLJEnNp86iotp+AYDXC5XPJ2+hXdTY2Ij77rsPU6ZMQVRUFKxWK4xGI0MRERH1KXaMFEI1eMgPkzkGPMUkAaoIqEae1++1BevLL7+E0WjE559/DrPZjAULFiA6OlrusoiIaADgn98K4atxQxUb2xyOmpqaL9NvndzR6wWgau4cnTwpd6kd8nq92LJlCxITE1FTU4NXX30Vd999N0MRERH1G3aMFCLygguhiouDpFI13/6j9bQaAERFAoMGQRURiaifjJG30A7s3bsXGRkZ+PDDD2E0Gv0zWRMREfUnBiOFiPrxaESPG4eGnTuh0owAoGqe5DEiEoiMgHTyFCJHj0bMhIlyl9qGJEl4/vnnsXLlSsTFxeH5559HcnKy3GUREdEAxVNpCqGKicGQeXMROer85tNoTS0zXzc1QWryIiJuOAbffhsiR/9I7lL9Dh8+jFtvvRU5OTm4/vrrsX37doYiIiKSFTtGCjL4+hT4vv8ep998E95jxwCooJIkqOLiEPvrKzEsM0PuEv1KSkqwcOFCREZG4qmnnsI111wjd0lEREQMRkqiiorC0KxMxF5xBer/+U94v/8eEfHxiElOQsykSVDFxMhdIlwuFxYtWoTi4mKkpKTgvvvuw4gRI+Qui4iICACDkeKoVCpEj78E0eMvkbuUdv76178iKysLp06dwqZNm3DzzTdz9moiIgopHGNEfe7EiRMwm8246aabkJCQALvdjltuuYWhiIiIQg47RtSnPvzwQ6Snp+P777/HfffdB4PBwEBEREQhix0j6hP19fXIzc3F1KlTER8fD5vNhjlz5jAUERFRSGPHiHrdrl27kJaWht27d2P58uXIzMxEZGSk3GURERGdEztG1Guamprw0EMPITk5GY2NjXjttdeQnZ3NUERERGGDHSPqFd988w2MRiPKy8uRnZ3Ne5wREVFYYseIesTn8+Gpp57CxIkTceTIEbz00ktYtmwZQxEREYUldoyo2w4cOIDMzEy8++67uPPOO7Fy5UoMGTJE7rKIiIi6jcGIgiZJEl555RUsW7YMgwYNQmFhIa688kq5yyIiIuoxnkqjoBw7dgx33HEHMjIyMHXqVNjtdoYiIiJSDHaMqMv+8Ic/ICcnB16vF48//jh0Op3cJREREfUqdozonGpqapCRkYFZs2bh0ksvhd1uZygiIiJFYseIOvXuu+8iIyMD1dXVePDBBzFjxgzOXk1ERIrFjhEFdOrUKSxbtgzTp0/H6NGjYbfbMXPmTIYiIiJSNHaMqJ1PP/0UaWlp2L9/PywWC+bNm4eICGZoIiJSPn7akV9DQwPWrVuHK6+8ErGxsSgpKUFaWhpDERERDRjsGBEA4IsvvoDRaMQXX3yBhQsXYsGCBYiK4o8HERENLGwFDHBerxcFBQVITExEbW0ttm3bhoULFzIUERHRgMRPvwFMFEWkp6fj448/Rnp6OpYuXYrY2Fi5yyIiIpINg9EAJEkStm7dilWrVkGtVuOFF15AUlKS3GURERHJjsFogDl06BCys7Px17/+FbNnz8aaNWswdOhQucsiIiIKCSEXjERRhM1mgyAIEEURJpMJarU66HW7u0zJiouLcffddyMyMhJPP/00pk6dKndJREREISXkgpHBYEBFRQWA5gCTk5MDq9Ua9LrdXaZEVVVVWLRoEUpKSnDDDTdg3bp1AyIIEhERBSukgpEoim0eC4IAh8MR9LrdXaZEb731FubPn4/Tp0+joKAAN910k9wlERGRjCRJgiRJcpfR61QqVa/cnSGkLtd3OBzQaDRtntNoNHA6nUGt291lZ6uvr4fH42nzFS5qa2uRk5ODm2++GQkJCbDb7QxFREQDnCRJaGpq8ocjJX15vV7/e+uJkOoYud3ugM+7XK6g1u3usrPl5+dj/fr1AdcPZf/4xz+Qnp6Oo0eP4r777oPBYOA9zoiIBrjW8DB48GBccMEFivpckCQJdXV1OHbsGLxeb4/m4gupYNSRjsJMsOsGuywvLw+rVq3yP/Z4PBgzZkyXa+lvdXV1WLduHR577DFMnDgRzzzzTEjXS0RE/UulUmHkyJEYNGiQ3KX0utb3dOTIEUiS1O3gF1LBSK1Wt+vcuFyugAOFO1u3u8vOFhsbGzYTHu7cuRNpaWnYvXs3Vq5cifT0dERGRspdFhERhYjWU0xKvrPBoEGDoFKpehSMQmqMkU6nC/h8oMkHO1u3u8vCUVNTEx588EFcfvnl8Pl8eOONN5CVlcVQREREASnpFNrZeuO9hVRsFAShzWNRFJGUlOTv5jidTqjVagiC0Om6Z3d/uros3Hz99dcwGo2oqKhAdnY27r77bkRHR8tdFhERhRHvoUOQXNX9tj+VZgQiR4/ut/0FK6SCEQBYrVZYLBYkJyejrKyszfxC+fn5SE5ORm5u7jnX7e6ycODz+fDUU0/BYrHgggsuwMsvv4wJEybIXRYREYUZ76FDqEm5Hqhv6L+dxsYgvvSdLoejvXv3YseOHRg7diz27t2L+fPn92kzQyUpcTKDPuLxeBAfH4+amhrExcXJUsP+/fuRmZmJ9957D/PmzcOKFSswePBgWWohIqLw4fP5IEkSfvrTn/rHzzZ98SU8M2f1ey1xf9iBqEvHd2ndK6+8Eh9//DGA5pC0bt06vPrqqwHXra+vx759+6BSqRAR8cNooRMnTmDKlCld+vwOuY4RBSZJEl5++WUsW7YMQ4YMwdatWzF58mS5yyIiIuoze/fubfN47NixeO+99/p0nyE1+JoCO3r0KGbOnInMzExce+212L59O0MREREp3rvvvttuUuYRI0Zg586dfbZPdoxC3I4dO2AymeDz+fDEE09g+vTpcpdERETUL2pqagI+X13dd4PF2TEKUW63G+np6bjjjjtw2WWXwW63MxQRERGh48DUG9gxCkEOhwOZmZmoqanBQw89hNtvv13R804QEREFEh8f325S5urqasTHx/fZPtkxCiGnTp3CkiVLkJKSgosuugh2ux0zZsxgKCIiogHpuuuuC/i8Vqvts32yYxQi/vnPfyI9PR379+/H2rVrMXfu3DaXGhIREQ00Y8eObfN479690Gq1fTqPEYORzBoaGrB+/Xps2rQJ48ePh9VqbfeDQERENFC9+uqrWLduHRITE1FRUYFt27b16f4YjGT073//G0ajEV9++SUWLVqE7OxsRd/cj4iIQotKMwKIjen3ma9VmhFdXn3s2LF46KGHAACzZvX9ZJT8FJaB1+vFli1b8Jvf/AY/+clP8Nprr+FXv/qV3GUREdEAEzl6NOJL3+G90s7AYNTPKisrkZ6ejk8++QSZmZlYvHixf2p2IiKi/hY5ejQQwkGlvzEY9RNJklBYWIjVq1djxIgRePHFF5GYmCh3WURERHQGBqN+8N133yE7Oxtvv/029Ho9Vq9ejaFDh8pdFhEREZ2FwaiPvfHGG1i0aBGioqLw9NNPY+rUqXKXREREA5gkSXKX0Gd6471xopw+UlVVhTlz5mDu3LmYPHkyduzYwVBERESyaZ0suKmpSeZK+k5dXR0kSerRxMjsGPWBv/zlL8jOzkZdXR0eeeQR3HjjjXKXREREBEmSUFVVhaioKEXdVUGSJNTV1eHYsWMAwGAUKmpra7Fy5Uo8//zzuOqqq7B+/XqMGjVK7rKIiIigUqkQGRmJ06dPY//+/XKX0+taT6NFRkb2aDsMRr3k73//O9LT03Hs2DHcf//9mD17tqLSOBERhT+VSuWfSFhpY40iIiJ65XOXwaiH6urqcO+99+KJJ57ApEmT8Oyzz2LMmDFyl0VERBRQa3jgH++BMRj1QEVFBYxGI/bs2YNVq1bBaDT2uIVHRERE8uFVad3Q2NiI9evXY/LkyQCA4uJiZGZmMhQRERGFOXaMuiElJQWff/45FixYALPZjOjoaLlLIiIiol7AYBSE1oFqx44dQ2FhIcaPH4/6+nrU19fLXBkRERF15OTJkwC6NuBcJSltWHofOnjwIAdWExERhakDBw7goosu6nQdBqMg+Hw+HDp0CMOHDw+L0fwejwdjxozBgQMHEBcXJ3c5suFxaMbj8AMei2Y8Ds14HH6g1GMhSRJqa2sxevRoRER0Pryap9KCEBERcc6kGYri4uIU9QPeXTwOzXgcfsBj0YzHoRmPww+UeCzi4+O7tB6vSiMiIiJqwWBERERE1ILBSMFiY2Nx//33IzY2Vu5SZMXj0IzH4Qc8Fs14HJrxOPyAx4KDr4mIiIj82DEiIiIiasFgRERERNSCwYiIiIioBecxCnOiKMJms0EQBIiiCJPJBLVa3eN1w02w783pdCInJwcVFRX9V2Q/COY4OJ1OOBwOAEBZWRm2bt2qmJ8HILhj0Xoc3G43ysrKkJqaCq1W24/V9p3u/ru3WCzIy8sbsD8TTqcTAKDVaiGKItxu94D9mXA4HBBFEYIgAAB0Ol0/VSoTicKaVqv1/39lZaWk1+t7Zd1wE8x7s1qtUkVFhaTEH/9gjsPmzZvb/P+Zr1WCYI6FWq2WKioqJEmSpMLCQkkQhD6vr790599967+P6urqPqys/wVzLEwmkwRAAiDpdDpFHYtgjkNpaalkMpn86yrp30ZHlPfJMIBUVla2+zBTq9U9XjfcdPe9KS0YBXMcKioq2iyrrKyUAEiVlZV9WmN/CfZnorS01P//hYWFigmJ3f23YbVaJUEQFBUGgj0WhYWFUnV1taKOgSQFfxzO/jlQyu+IznCMURhzOBzQaDRtntNoNP4WcHfXDTdKfm/BCOY4aLVabN261f/Y7Xb711eCYH8mzjw1YLVaYTab+7S+/tKdfxs2mw16vb6vS+t33TkWarVaUacSgeCOgyiKcLlcUKvVcDqdcLvd/tNpSsYxRmGs9cPsbC6Xq0frhhslv7dgBHsczvzwKy4uhk6nU8yHQHd+JpxOJ4qLi5GSkgKTydRHlfWvYI+D2+1WzM/A2bpzLGw2G4DmMXhms1kRoSCY4+B0OqHRaGCz2aDT6VBUVARBEBQZnM/EYKRAHf3g93TdcKPk9xaMcx2H1g8ApQ1ED6SzY6HVaiEIAiwWi2K7Jq06Og4lJSWKCYVd1dGxOHNAsiAISElJQWVlZf8V1s8CHQeXywVRFP1/NJlMJowYMQKSwueF5qm0MKZWq9ul/Na2Z0/WDTdKfm/B6O5xsFgsKC0tVdTx6u6xUKvVMBgMMBgMigjWwRwHh8OBOXPm9FNl/S/YnwlRFP3/33r11pnPhatgjoMgCG1OJ7b+V+nDFBiMwlhHl0wmJSX1aN1wo+T3FozuHIeCggJYLBYIggC3262IMAAEdywcDgdGjBjhf9x6ukQJH4LB/kyUlJSgqKgIRUVFEEUR+fn5ivkQDOZYOJ1OTJ8+vd3zShiDF8xxUMKpw+5gMApjZ//QiqKIpKSkNqm+9Zf7udYNZ8Ech7MpJQgAwR8Hm83mP33kdrtRUlKiiJ8HILhjodFo2nxYOJ1OqNVqRcxZE8xx0Ol0MJlM/i8AMJvNijgOQPC/Lzdv3uxf1+FwQK/XK+LfR7DHISkpyf97snUuI6X8THSEN5ENc6IoorCwEMnJySgrK2szIZvBYEBycjJyc3PPuW64C+Y4OBwOlJaWoqCgALm5uUhOTlbMeJKuHgdRFJGQkNDmtWq1GtXV1TJU3TeC+Zmw2Wz+0wulpaXYvHmzYv5aDuY4AM1/LBQVFcFiscBkMikqHAVzLFonQFWr1aisrGwTlMJdMMfB7XbDYrEgMTERFRUV/g6zkjEYEREREbXgqTQiIiKiFgxGRERERC0YjIiIiIhaMBgRERERtWAwIiIiImrBYERERETUgsGIiIiIqAWDERGFDYPBgMTERKhUKqhUKv99zVq/UlJSYDAYAt7Gwu12IyEhAQUFBb1Wj9PpRGJiIkaMGIGUlJRe2y4RySdK7gKIiLrKarUCAFQqFXQ6nf/xmQoKCpCYmIjCwsI2d4pvvVN4WVlZt/ZtsVjazX6s1WpRUVGBlJSUdjfmJKLwxI4RESlKbm4utFotzGZzuzukS5IUMEx1RWc3lVX6LRKIBhIGIyJSnNabwjocjl7Zns1mU9QNh4moYwxGRKQ4Z94dvDe2lZOT0+PtEFF44BgjIlIUURRhs9mg0+n8nSNRFGE2m1FeXg5BEFBRUQGguaNksVggiiJMJhNSU1PhcDhQWloKs9kMACguLgYAlJeXw2AwAGgOXIHutt56Z3oAKCsr63A9IgpdKkmSJLmLICIKRuvg69LSUv9zbrfbH3T0en3AQNI6SLo1GLVKSEiATqdDQkICcnNzkZCQAK1W6x+P1HrF2Zn7O1Nr6EpNTUVubq7/+REjRiAvL6/Nc0QU2tgxIqKwJIpiu0vvKysrodPpkJycHPA1giAEvHpMEASUlJRg7969AICKigqo1eqg6nE6ne0GdiclJaG4uJjBiCiMMBgRUVgSBCFg4HC73UhMTER+fn67ztC5ttcahoINRa2vOXtMU0dBjIhCFwdfE5GiqNVqFBYWwul0+scJdYVGo+nRfnv6eiIKDQxGRKQ4rYOuWwdCd0UwXaJAl+53p8tERKGHwYiIKEjl5eWdTvhIROGLwYiIFMdmswEA9Hp9r2zv7G6Q2+3mqTMihWIwIiJFcTqdyMnJgSAI2Lp1a7vlHc1g3dnM1ikpKSgvL/c/FkWxTVhyuVwBX9/R80QUujiPERGFDYPBAFEU4XQ6ATR3hFo7N603iQXQbj4hURRhsVjgcDjgdruh1+uRl5cHAMjPz2/TYUpOTg54tZvFYoHT6URKSgq0Wi10Ol3A7ZrNZgiC0OZ5nU4Hs9ncax0sIuo7DEZERERELXgqjYiIiKgFgxERERFRCwYjIiIiohYMRkREREQtGIyIiIiIWjAYEREREbVgMCIiIiJqwWBERERE1ILBiIiIiKgFgxERERFRCwYjIiIiohYMRkREREQtGIyIiIiIWvx/5tqA+80HWZUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the persistence diagram\n",
    "gd.plot_persistence_diagram(persistence_2)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can use various clustering algorithms to group similar protein sequences together based on their persistence diagrams."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Silhouette Coefficient for KMeans: 0.781\n",
      "Silhouette Coefficient for Agglomerative Clustering: 0.781\n",
      "Silhouette Coefficient for DBSCAN: 0.781\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/amelieschreiber/vscode_projects/vision_persistent_homology/.conda/lib/python3.11/site-packages/sklearn/cluster/_kmeans.py:870: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning\n",
      "  warnings.warn(\n",
      "/Users/amelieschreiber/vscode_projects/vision_persistent_homology/.conda/lib/python3.11/site-packages/sklearn/cluster/_agglomerative.py:983: FutureWarning: Attribute `affinity` was deprecated in version 1.2 and will be removed in 1.4. Use `metric` instead\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "# Perform KMeans clustering\n",
    "kmeans = KMeans(n_clusters=3, random_state=0).fit(wasserstein_distances)\n",
    "labels_kmeans = kmeans.labels_\n",
    "print(\"Silhouette Coefficient for KMeans: %0.3f\" % silhouette_score(wasserstein_distances, labels_kmeans))\n",
    "\n",
    "# Perform Agglomerative Clustering\n",
    "agg_clustering = AgglomerativeClustering(n_clusters=3, affinity='precomputed', linkage='average').fit(wasserstein_distances)\n",
    "labels_agg = agg_clustering.labels_\n",
    "print(\"Silhouette Coefficient for Agglomerative Clustering: %0.3f\" % silhouette_score(wasserstein_distances, labels_agg))\n",
    "\n",
    "# Perform DBSCAN clustering\n",
    "dbscan = DBSCAN(metric=\"precomputed\", eps=0.170, min_samples=1).fit(wasserstein_distances)\n",
    "labels_dbscan = dbscan.labels_\n",
    "if len(set(labels_dbscan)) > 1:  # More than 1 cluster\n",
    "    print(\"Silhouette Coefficient for DBSCAN: %0.3f\" % silhouette_score(wasserstein_distances, labels_dbscan))\n",
    "else:\n",
    "    print(\"Cannot compute Silhouette Coefficient for DBSCAN as there is only one cluster.\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Clusters for KMeans:\n",
      "Cluster 0: ['VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFF', 'VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFES', 'VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESF', 'GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDR', 'GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDRF', 'GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDRFG', 'KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNR', 'KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNRN', 'KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNRNT']\n",
      "Cluster 1: ['MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGER', 'MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERG', 'MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGF']\n",
      "Cluster 2: ['MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGER', 'MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERG', 'MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGF']\n",
      "\n",
      "Clusters for Agglomerative Clustering:\n",
      "Cluster 0: ['VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFF', 'VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFES', 'VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESF', 'GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDR', 'GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDRF', 'GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDRFG', 'KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNR', 'KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNRN', 'KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNRNT']\n",
      "Cluster 1: ['MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGER', 'MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERG', 'MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGF']\n",
      "Cluster 2: ['MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGER', 'MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERG', 'MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGF']\n",
      "\n",
      "Clusters for DBSCAN:\n",
      "Cluster 0: ['VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFF', 'VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFES', 'VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESF', 'GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDR', 'GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDRF', 'GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDRFG', 'KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNR', 'KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNRN', 'KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNRNT']\n",
      "Cluster 1: ['MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGER', 'MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERG', 'MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGF']\n",
      "Cluster 2: ['MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGER', 'MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERG', 'MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGF']\n"
     ]
    }
   ],
   "source": [
    "# Create a list of protein sequences that corresponds to the order of the Wasserstein distance matrix\n",
    "sequences = [\n",
    "    # Protein family 1: Hemoglobin\n",
    "    \"VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFF\",\n",
    "    \"VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFES\",\n",
    "    \"VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESF\",\n",
    "\n",
    "    # Protein family 2: Insulin\n",
    "    \"MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGER\",\n",
    "    \"MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERG\",\n",
    "    \"MALWMRFLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGF\",\n",
    "\n",
    "    # Protein family 3: Albumin\n",
    "    \"MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGER\",\n",
    "    \"MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERG\",\n",
    "    \"MKWVTVRPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGF\",\n",
    "\n",
    "    # Protein family 4: Myoglobin\n",
    "    \"GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDR\",\n",
    "    \"GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDRF\",\n",
    "    \"GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDRFG\",\n",
    "\n",
    "    # Protein family 5: Lysozyme\n",
    "    \"KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNR\",\n",
    "    \"KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNRN\",\n",
    "    \"KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNRNT\"\n",
    "]\n",
    "\n",
    "\n",
    "\n",
    "# Print the clusters for KMeans\n",
    "print(\"Clusters for KMeans:\")\n",
    "kmeans_clusters = {i: [] for i in set(labels_kmeans)}\n",
    "for sequence, label in zip(sequences, labels_kmeans):\n",
    "    kmeans_clusters[label].append(sequence)\n",
    "for label, cluster in kmeans_clusters.items():\n",
    "    print(f\"Cluster {label}: {cluster}\")\n",
    "\n",
    "# Print the clusters for Agglomerative Clustering\n",
    "print(\"\\nClusters for Agglomerative Clustering:\")\n",
    "agg_clusters = {i: [] for i in set(labels_agg)}\n",
    "for sequence, label in zip(sequences, labels_agg):\n",
    "    agg_clusters[label].append(sequence)\n",
    "for label, cluster in agg_clusters.items():\n",
    "    print(f\"Cluster {label}: {cluster}\")\n",
    "\n",
    "# Print the clusters for DBSCAN\n",
    "print(\"\\nClusters for DBSCAN:\")\n",
    "dbscan_clusters = {i: [] for i in set(labels_dbscan)}\n",
    "for sequence, label in zip(sequences, labels_dbscan):\n",
    "    dbscan_clusters[label].append(sequence)\n",
    "for label, cluster in dbscan_clusters.items():\n",
    "    print(f\"Cluster {label}: {cluster}\")\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compare the clusters given by different clustering algorithms, we can use a measure called the Adjusted Rand Index (ARI). The ARI is a function that measures the similarity of the two assignments, ignoring permutations and with chance normalization. The maximum ARI score is 1, indicating that the two sets of clusters are identical. Next, we can calculate the ARI for each pair of clustering results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adjusted Rand Index for KMeans and Agglomerative Clustering: 1.0\n",
      "Adjusted Rand Index for KMeans and DBSCAN: 1.0\n",
      "Adjusted Rand Index for Agglomerative Clustering and DBSCAN: 1.0\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import adjusted_rand_score\n",
    "\n",
    "# Calculate ARIs\n",
    "ari_kmeans_agg = adjusted_rand_score(labels_kmeans, labels_agg)\n",
    "ari_kmeans_dbscan = adjusted_rand_score(labels_kmeans, labels_dbscan)\n",
    "ari_agg_dbscan = adjusted_rand_score(labels_agg, labels_dbscan)\n",
    "\n",
    "# Print ARIs\n",
    "print(f\"Adjusted Rand Index for KMeans and Agglomerative Clustering: {ari_kmeans_agg}\")\n",
    "print(f\"Adjusted Rand Index for KMeans and DBSCAN: {ari_kmeans_dbscan}\")\n",
    "print(f\"Adjusted Rand Index for Agglomerative Clustering and DBSCAN: {ari_agg_dbscan}\")\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ARI scores will give you a measure of the similarity between the clusters found by the different algorithms. Please note that while the ARI can give a general sense of the similarity between different sets of clusters, it doesn't provide a detailed comparison of the clusters. For example, it won't tell you which specific clusters in one algorithm correspond to the clusters in another algorithm. To do that, you'd need to look more closely at the specific elements in each cluster and see where they appear in the other set of clusters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}