Created
October 31, 2021 17:33
-
-
Save tianyu-lu/4e445cef3262055dca5422bc9cc60cff to your computer and use it in GitHub Desktop.
2-Discriminative-Models.ipynb
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "nbformat": 4, | |
| "nbformat_minor": 0, | |
| "metadata": { | |
| "colab": { | |
| "name": "2-Discriminative-Models.ipynb", | |
| "provenance": [], | |
| "collapsed_sections": [], | |
| "toc_visible": true, | |
| "authorship_tag": "ABX9TyPyus4koHrqTFyr0l7Pv0XV", | |
| "include_colab_link": true | |
| }, | |
| "kernelspec": { | |
| "name": "python3", | |
| "display_name": "Python 3" | |
| }, | |
| "accelerator": "GPU", | |
| "widgets": { | |
| "application/vnd.jupyter.widget-state+json": { | |
| "6ad08c7883434b4ba5f8b212e2551ad1": { | |
| "model_module": "@jupyter-widgets/controls", | |
| "model_name": "HBoxModel", | |
| "model_module_version": "1.5.0", | |
| "state": { | |
| "_view_name": "HBoxView", | |
| "_dom_classes": [], | |
| "_model_name": "HBoxModel", | |
| "_view_module": "@jupyter-widgets/controls", | |
| "_model_module_version": "1.5.0", | |
| "_view_count": null, | |
| "_view_module_version": "1.5.0", | |
| "box_style": "", | |
| "layout": "IPY_MODEL_2c61c420ee32463bbfb2d6e8de59b5a8", | |
| "_model_module": "@jupyter-widgets/controls", | |
| "children": [ | |
| "IPY_MODEL_636706f165014540b755cd9b6c69a70a", | |
| "IPY_MODEL_28cb297422c044edbf288bd5c2d16d80" | |
| ] | |
| } | |
| }, | |
| "2c61c420ee32463bbfb2d6e8de59b5a8": { | |
| "model_module": "@jupyter-widgets/base", | |
| "model_name": "LayoutModel", | |
| "model_module_version": "1.2.0", | |
| "state": { | |
| "_view_name": "LayoutView", | |
| "grid_template_rows": null, | |
| "right": null, | |
| "justify_content": null, | |
| "_view_module": "@jupyter-widgets/base", | |
| "overflow": null, | |
| "_model_module_version": "1.2.0", | |
| "_view_count": null, | |
| "flex_flow": null, | |
| "width": null, | |
| "min_width": null, | |
| "border": null, | |
| "align_items": null, | |
| "bottom": null, | |
| "_model_module": "@jupyter-widgets/base", | |
| "top": null, | |
| "grid_column": null, | |
| "overflow_y": null, | |
| "overflow_x": null, | |
| "grid_auto_flow": null, | |
| "grid_area": null, | |
| "grid_template_columns": null, | |
| "flex": null, | |
| "_model_name": "LayoutModel", | |
| "justify_items": null, | |
| "grid_row": null, | |
| "max_height": null, | |
| "align_content": null, | |
| "visibility": null, | |
| "align_self": null, | |
| "height": null, | |
| "min_height": null, | |
| "padding": null, | |
| "grid_auto_rows": null, | |
| "grid_gap": null, | |
| "max_width": null, | |
| "order": null, | |
| "_view_module_version": "1.2.0", | |
| "grid_template_areas": null, | |
| "object_position": null, | |
| "object_fit": null, | |
| "grid_auto_columns": null, | |
| "margin": null, | |
| "display": null, | |
| "left": null | |
| } | |
| }, | |
| "636706f165014540b755cd9b6c69a70a": { | |
| "model_module": "@jupyter-widgets/controls", | |
| "model_name": "FloatProgressModel", | |
| "model_module_version": "1.5.0", | |
| "state": { | |
| "_view_name": "ProgressView", | |
| "style": "IPY_MODEL_2bc74267a9e84d07bfa690f768435221", | |
| "_dom_classes": [], | |
| "description": "100%", | |
| "_model_name": "FloatProgressModel", | |
| "bar_style": "success", | |
| "max": 200, | |
| "_view_module": "@jupyter-widgets/controls", | |
| "_model_module_version": "1.5.0", | |
| "value": 200, | |
| "_view_count": null, | |
| "_view_module_version": "1.5.0", | |
| "orientation": "horizontal", | |
| "min": 0, | |
| "description_tooltip": null, | |
| "_model_module": "@jupyter-widgets/controls", | |
| "layout": "IPY_MODEL_ef6cd16bfebb4a9185d1bd824b870bbf" | |
| } | |
| }, | |
| "28cb297422c044edbf288bd5c2d16d80": { | |
| "model_module": "@jupyter-widgets/controls", | |
| "model_name": "HTMLModel", | |
| "model_module_version": "1.5.0", | |
| "state": { | |
| "_view_name": "HTMLView", | |
| "style": "IPY_MODEL_604bce716ef249f1b1e536759c0b5976", | |
| "_dom_classes": [], | |
| "description": "", | |
| "_model_name": "HTMLModel", | |
| "placeholder": "", | |
| "_view_module": "@jupyter-widgets/controls", | |
| "_model_module_version": "1.5.0", | |
| "value": " 200/200 [00:09<00:00, 21.72it/s]", | |
| "_view_count": null, | |
| "_view_module_version": "1.5.0", | |
| "description_tooltip": null, | |
| "_model_module": "@jupyter-widgets/controls", | |
| "layout": "IPY_MODEL_4063eabcefb041d1b6aaa7c4dc984e3d" | |
| } | |
| }, | |
| "2bc74267a9e84d07bfa690f768435221": { | |
| "model_module": "@jupyter-widgets/controls", | |
| "model_name": "ProgressStyleModel", | |
| "model_module_version": "1.5.0", | |
| "state": { | |
| "_view_name": "StyleView", | |
| "_model_name": "ProgressStyleModel", | |
| "description_width": "initial", | |
| "_view_module": "@jupyter-widgets/base", | |
| "_model_module_version": "1.5.0", | |
| "_view_count": null, | |
| "_view_module_version": "1.2.0", | |
| "bar_color": null, | |
| "_model_module": "@jupyter-widgets/controls" | |
| } | |
| }, | |
| "ef6cd16bfebb4a9185d1bd824b870bbf": { | |
| "model_module": "@jupyter-widgets/base", | |
| "model_name": "LayoutModel", | |
| "model_module_version": "1.2.0", | |
| "state": { | |
| "_view_name": "LayoutView", | |
| "grid_template_rows": null, | |
| "right": null, | |
| "justify_content": null, | |
| "_view_module": "@jupyter-widgets/base", | |
| "overflow": null, | |
| "_model_module_version": "1.2.0", | |
| "_view_count": null, | |
| "flex_flow": null, | |
| "width": null, | |
| "min_width": null, | |
| "border": null, | |
| "align_items": null, | |
| "bottom": null, | |
| "_model_module": "@jupyter-widgets/base", | |
| "top": null, | |
| "grid_column": null, | |
| "overflow_y": null, | |
| "overflow_x": null, | |
| "grid_auto_flow": null, | |
| "grid_area": null, | |
| "grid_template_columns": null, | |
| "flex": null, | |
| "_model_name": "LayoutModel", | |
| "justify_items": null, | |
| "grid_row": null, | |
| "max_height": null, | |
| "align_content": null, | |
| "visibility": null, | |
| "align_self": null, | |
| "height": null, | |
| "min_height": null, | |
| "padding": null, | |
| "grid_auto_rows": null, | |
| "grid_gap": null, | |
| "max_width": null, | |
| "order": null, | |
| "_view_module_version": "1.2.0", | |
| "grid_template_areas": null, | |
| "object_position": null, | |
| "object_fit": null, | |
| "grid_auto_columns": null, | |
| "margin": null, | |
| "display": null, | |
| "left": null | |
| } | |
| }, | |
| "604bce716ef249f1b1e536759c0b5976": { | |
| "model_module": "@jupyter-widgets/controls", | |
| "model_name": "DescriptionStyleModel", | |
| "model_module_version": "1.5.0", | |
| "state": { | |
| "_view_name": "StyleView", | |
| "_model_name": "DescriptionStyleModel", | |
| "description_width": "", | |
| "_view_module": "@jupyter-widgets/base", | |
| "_model_module_version": "1.5.0", | |
| "_view_count": null, | |
| "_view_module_version": "1.2.0", | |
| "_model_module": "@jupyter-widgets/controls" | |
| } | |
| }, | |
| "4063eabcefb041d1b6aaa7c4dc984e3d": { | |
| "model_module": "@jupyter-widgets/base", | |
| "model_name": "LayoutModel", | |
| "model_module_version": "1.2.0", | |
| "state": { | |
| "_view_name": "LayoutView", | |
| "grid_template_rows": null, | |
| "right": null, | |
| "justify_content": null, | |
| "_view_module": "@jupyter-widgets/base", | |
| "overflow": null, | |
| "_model_module_version": "1.2.0", | |
| "_view_count": null, | |
| "flex_flow": null, | |
| "width": null, | |
| "min_width": null, | |
| "border": null, | |
| "align_items": null, | |
| "bottom": null, | |
| "_model_module": "@jupyter-widgets/base", | |
| "top": null, | |
| "grid_column": null, | |
| "overflow_y": null, | |
| "overflow_x": null, | |
| "grid_auto_flow": null, | |
| "grid_area": null, | |
| "grid_template_columns": null, | |
| "flex": null, | |
| "_model_name": "LayoutModel", | |
| "justify_items": null, | |
| "grid_row": null, | |
| "max_height": null, | |
| "align_content": null, | |
| "visibility": null, | |
| "align_self": null, | |
| "height": null, | |
| "min_height": null, | |
| "padding": null, | |
| "grid_auto_rows": null, | |
| "grid_gap": null, | |
| "max_width": null, | |
| "order": null, | |
| "_view_module_version": "1.2.0", | |
| "grid_template_areas": null, | |
| "object_position": null, | |
| "object_fit": null, | |
| "grid_auto_columns": null, | |
| "margin": null, | |
| "display": null, | |
| "left": null | |
| } | |
| } | |
| } | |
| } | |
| }, | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "view-in-github", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "<a href=\"https://colab.research.google.com/gist/tianyu-lu/4e445cef3262055dca5422bc9cc60cff/2-discriminative-models.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "AspycAxD1_QE" | |
| }, | |
| "source": [ | |
| "## 2. Discriminative Sequence Models" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "vmKK-dmN0bo5", | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "outputId": "7355ad64-020a-49bb-cefd-a76d74318207" | |
| }, | |
| "source": [ | |
| "!git clone https://github.com/igemto-drylab/CSBERG-ML.git\n", | |
| "%cd CSBERG-ML\n", | |
| "from util import *\n", | |
| "import tqdm.notebook as tq" | |
| ], | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "text": [ | |
| "Cloning into 'CSBERG-ML'...\n", | |
| "remote: Enumerating objects: 176, done.\u001b[K\n", | |
| "remote: Counting objects: 100% (176/176), done.\u001b[K\n", | |
| "remote: Compressing objects: 100% (169/169), done.\u001b[K\n", | |
| "remote: Total 176 (delta 83), reused 23 (delta 4), pack-reused 0\u001b[K\n", | |
| "Receiving objects: 100% (176/176), 17.11 MiB | 21.08 MiB/s, done.\n", | |
| "Resolving deltas: 100% (83/83), done.\n", | |
| "/content/CSBERG-ML\n" | |
| ], | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "W4eP8HJyiYY2" | |
| }, | |
| "source": [ | |
| "Learning outcome: able to implement simple neural network architectures and train them on data. Students should be able to know how to write code for a simple neural network in PyTorch and how it’s trained." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "03M13bDgRyNB" | |
| }, | |
| "source": [ | |
| "A discriminative model is a parameterization for $p(y|\\mathbf{x})$. In our case, $\\mathbf{x}$ is a protein sequence and $y$ is a real number. More precisely, we can say that $\\mathbf{x} \\in \\mathbb{R}^{B \\times L \\times 20}$, where $B$ is the batch size $L$ is the sequence length, and 20 is the number of canonical amino acids. $y \\in \\mathbb{R}$ can be any quantitative property of a protein, e.g. an enzyme's catalytic activity, thermostability, $\\Delta \\Delta G$ when folding, etc. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "HUqhDKqkTGzT" | |
| }, | |
| "source": [ | |
| "What does it mean to \"train\" a discriminative model? In most cases, this means you want to find parameters $\\boldsymbol{\\theta}$ such that you satisfy \n", | |
| "\n", | |
| "$$\\min_{\\boldsymbol{\\theta}} \\mathbb{E}_{y}[\\sum_{i = 1}^{N} \\big((y_i - p(y|\\mathbf{x}_i))^2 \\big)] \\;\\; \\text{given } \\mathcal{D} = \\{ (\\mathbf{x}_1, y_1), \\cdots , (\\mathbf{x}_N, y_N) \\}$$\n", | |
| "\n", | |
| "In words, this means you are given a dataset of $N$ sequences $\\mathbf{x}_1 \\cdots \\mathbf{x}_N$, each having a corresponding property values $y_1 \\cdots y_N$, you feed these sequences into a neural network $p(y|\\mathbf{x})$, compute its difference with the actual property values and minimize this squared difference. Some terms to describe this function are the loss function, objective function, error function." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "-CE1wU0o1cNg" | |
| }, | |
| "source": [ | |
| "Before diving into any code, here's a very useful playground to deepen your intuition on neural networks and what they do: https://playground.tensorflow.org/ . Play around with the nonlinearities and see how well they do on the circle dataset." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "ht884nFGqs-1" | |
| }, | |
| "source": [ | |
| "### Feed-forward Neural Network" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "pPSHDtZYdBW-" | |
| }, | |
| "source": [ | |
| "Random numbers on a computer are not truly random. A random number generator is deterministic given a seed. They are useful when we want to fix the randomness in our model weight initialization, train/test splits, etc. by setting seeds. However, it's probably not good if you have to rely on picking 'good seeds' for your model to train." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "XKBXVv97dBHA" | |
| }, | |
| "source": [ | |
| "seed_val = 105716977\n", | |
| "random.seed(seed_val)\n", | |
| "np.random.seed(seed_val)\n", | |
| "torch.manual_seed(seed_val)\n", | |
| "torch.cuda.manual_seed_all(seed_val)" | |
| ], | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "Vo49wDP2iZ9G" | |
| }, | |
| "source": [ | |
| "class SeqfuncData(Dataset):\n", | |
| "\n", | |
| " def __init__(self, X, y, scale_X=True):\n", | |
| " if not torch.is_tensor(X):\n", | |
| " self.X = torch.from_numpy(X)\n", | |
| " else:\n", | |
| " self.X = X\n", | |
| " if not torch.is_tensor(y):\n", | |
| " self.y = torch.from_numpy(y)\n", | |
| " else:\n", | |
| " self.y = y\n", | |
| " self.X = self.X.type(torch.FloatTensor)\n", | |
| " self.y = self.y.type(torch.FloatTensor)\n", | |
| "\n", | |
| " def __len__(self):\n", | |
| " return len(self.X)\n", | |
| "\n", | |
| " def __getitem__(self, idx):\n", | |
| " return self.X[idx], self.y[idx]\n", | |
| "\n", | |
| "\n", | |
| "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")" | |
| ], | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "nTA7hHj1k1g4" | |
| }, | |
| "source": [ | |
| "Here is the code for training a simple feed-forward neural network as a protein sequence to function model. The operations of the neural network can be written as, for a single input protein sequence, \n", | |
| "\n", | |
| "$$\\hat{y} = \\tanh(\\mathbf{W}_2(\\text{ELU}(\\mathbf{W}_1\\mathbf{x} + \\mathbf{b}_1)) + \\mathbf{b}_2\n", | |
| "$$\n", | |
| "\n", | |
| "where $\\mathbf{W}_1 \\in \\mathbb{R}^{6258 \\times 20}$, $\\mathbf{W}_2 \\in \\mathbb{R}^{20 \\times 1}$ are the weight matrices and $\\mathbf{b}_1$, $\\mathbf{b}_2$ are the bias vectors.\n", | |
| "\n", | |
| "See [ELU](https://pytorch.org/docs/stable/generated/torch.nn.ELU.html) and [Tanh](https://pytorch.org/docs/stable/generated/torch.nn.Tanh.html#torch.nn.Tanh) for their definitions." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "ePJ7zx2H2GDc" | |
| }, | |
| "source": [ | |
| "def train_oracle(X, y, save_name, data_dim=298*21, hid_dim=20, batch_size=10, epochs=200,\n", | |
| " learning_rate=1e-4, show_preds=False):\n", | |
| " \"\"\"Raw values of y are fitted, they should be normalized prior to this\n", | |
| " function if you prefer.\n", | |
| " \"\"\"\n", | |
| " petase_dataset = SeqfuncData(X, y)\n", | |
| " train, test = train_test_split(list(range(X.shape[0])), test_size=.2)\n", | |
| "\n", | |
| " petase_train = DataLoader(petase_dataset, batch_size=batch_size,\n", | |
| " sampler=SubsetRandomSampler(train))\n", | |
| " petase_test = DataLoader(petase_dataset, batch_size=batch_size,\n", | |
| " sampler=SubsetRandomSampler(test))\n", | |
| " \n", | |
| " train_corrs, test_corrs = [], []\n", | |
| "\n", | |
| " for i in range(1):\n", | |
| " random.seed(i)\n", | |
| " model = torch.nn.Sequential(\n", | |
| " torch.nn.Linear(data_dim, hid_dim),\n", | |
| " torch.nn.ELU(),\n", | |
| " torch.nn.Linear(hid_dim, 1),\n", | |
| " torch.nn.Tanh(),\n", | |
| " )\n", | |
| " model.to(device)\n", | |
| " loss_fn = nn.MSELoss(reduction='sum')\n", | |
| "\n", | |
| " optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)\n", | |
| " scheduler = torch.optim.lr_scheduler.CyclicLR(optimizer, 1e-4, 1e-3, cycle_momentum=False)\n", | |
| " train_loss, test_loss = [], []\n", | |
| " best_test_loss = 1e3\n", | |
| " img_idx = 0\n", | |
| " for epoch in tq.tqdm(range(epochs)):\n", | |
| " total_train_loss, total_test_loss = 0, 0\n", | |
| " # for x, y in ps_train:\n", | |
| " for x, y in petase_train: \n", | |
| " model.train()\n", | |
| " optimizer.zero_grad()\n", | |
| "\n", | |
| " y_hat = model(x.reshape(-1,data_dim).to(device))\n", | |
| " loss = loss_fn(y.reshape(-1,1).to(device), y_hat)\n", | |
| " total_train_loss += loss.item()\n", | |
| "\n", | |
| " loss.backward()\n", | |
| " optimizer.step()\n", | |
| " scheduler.step()\n", | |
| " train_loss.append(total_train_loss / len(petase_train))\n", | |
| "\n", | |
| " for x, y in petase_test:\n", | |
| " model.eval()\n", | |
| " y_hat = model(x.reshape(-1,data_dim).to(device))\n", | |
| "\n", | |
| " loss = loss_fn(y.reshape(-1,1).to(device), y_hat)\n", | |
| " total_test_loss += loss.item()\n", | |
| " test_loss.append(total_test_loss / len(petase_test))\n", | |
| "\n", | |
| " # if (epoch+1) % 10 == 0:\n", | |
| " # print(\"Training loss: {}\".format(train_loss[-1]))\n", | |
| " # print(\"Validation loss: {}\".format(test_loss[-1]))\n", | |
| "\n", | |
| " if test_loss[-1] < best_test_loss:\n", | |
| " plt.clf()\n", | |
| " torch.save(model.state_dict(), \"{}_model.pt\".format(save_name))\n", | |
| " \n", | |
| " model.eval()\n", | |
| " train_preds = []\n", | |
| " train_reals = []\n", | |
| " test_preds, test_reals = [], []\n", | |
| " for x, y in petase_train: \n", | |
| " preds = model(x.reshape(-1,data_dim).to(device))\n", | |
| " preds = preds.detach().cpu().numpy()\n", | |
| " train_preds.extend(list(preds))\n", | |
| " reals = y\n", | |
| " train_reals.extend(list(reals))\n", | |
| "\n", | |
| " for x, y in petase_test:\n", | |
| " preds = model(x.reshape(-1,data_dim).to(device))\n", | |
| " preds = preds.detach().cpu().numpy()\n", | |
| " test_preds.extend(list(preds))\n", | |
| " reals = y\n", | |
| " test_reals.extend(list(reals))\n", | |
| " train_reals = [a.item() for a in train_reals]\n", | |
| " train_preds = [a.item() for a in train_preds]\n", | |
| " test_reals = [a.item() for a in test_reals]\n", | |
| " test_preds = [a.item() for a in test_preds]\n", | |
| "\n", | |
| " if show_preds:\n", | |
| " plt.scatter(train_reals, train_preds,label=\"Train\")\n", | |
| " plt.scatter(test_reals, test_preds,label=\"Test\")\n", | |
| " plt.xlabel(\"Actual log activity\")\n", | |
| " plt.ylabel(\"Predicted log activity\")\n", | |
| " plt.legend()\n", | |
| " # plt.savefig(\"scatter.png\", dpi=300)\n", | |
| " plt.show()\n", | |
| " # plt.savefig(\"images/{}.png\".format(img_idx), dpi=100)\n", | |
| " # img_idx += 1\n", | |
| " train_corrs.append(stats.spearmanr(train_reals, train_preds)[0])\n", | |
| " test_corrs.append(stats.spearmanr(test_reals, test_preds)[0])\n", | |
| " plt.clf()\n", | |
| " plt.plot(train_loss,label=\"Train\")\n", | |
| " plt.plot(test_loss,label=\"Test\")\n", | |
| " plt.legend()\n", | |
| " plt.xlabel(\"Epoch\")\n", | |
| " plt.ylabel(\"Mean Squared Error\")\n", | |
| " plt.title(\"{}\".format(save_name))\n", | |
| " plt.savefig(\"{}.png\".format(save_name), dpi=300)\n", | |
| " plt.show()\n", | |
| " print(\"\\nMinimum test loss:\", min(test_loss))\n", | |
| " print(\"\\nBest train correlation:\", max(train_corrs))\n", | |
| " print(\"\\nBest test correlation\", max(test_corrs))" | |
| ], | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "YyW8jygQoIMA" | |
| }, | |
| "source": [ | |
| "Let's get the required data from our Github repo. Note that the final activation function of our model is $\\tanh()$ which squishes outputs between $-1$ and $1$. Thus we must normalize our values within this range." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "8N7F23892H87", | |
| "colab": { | |
| "base_uri": "https://localhost:8080/", | |
| "height": 451, | |
| "referenced_widgets": [ | |
| "6ad08c7883434b4ba5f8b212e2551ad1", | |
| "2c61c420ee32463bbfb2d6e8de59b5a8", | |
| "636706f165014540b755cd9b6c69a70a", | |
| "28cb297422c044edbf288bd5c2d16d80", | |
| "2bc74267a9e84d07bfa690f768435221", | |
| "ef6cd16bfebb4a9185d1bd824b870bbf", | |
| "604bce716ef249f1b1e536759c0b5976", | |
| "4063eabcefb041d1b6aaa7c4dc984e3d" | |
| ] | |
| }, | |
| "outputId": "9ec8e8e0-9b6d-4e61-f3d3-3221f429458b" | |
| }, | |
| "source": [ | |
| "import time\n", | |
| "\n", | |
| "df = pd.read_csv(\"petase_activity.csv\")\n", | |
| "X_activity = get_X(list(df['sequence'])[:212])\n", | |
| "X_activity = X_activity.reshape(-1, 298, 21)\n", | |
| "y_activity = np.array(list(df['relative_activity']))[:212]\n", | |
| "y_activity = np.log(y_activity + 0.001)\n", | |
| "y_activity = ((y_activity - np.min(y_activity)) / (np.max(y_activity) - np.min(y_activity))) - 0.5 # don't let it saturate the sigmoid\n", | |
| "df = pd.read_csv(\"petase_stability.csv\")\n", | |
| "X_stability = get_X(list(df['sequence'])[:159])\n", | |
| "X_stability = X_stability.reshape(-1, 298, 21)\n", | |
| "y_stability = np.array(list(df['stability']))[:159]\n", | |
| "y_stability = np.log(y_stability + 0.001)\n", | |
| "y_stability = ((y_stability - np.min(y_stability)) / (np.max(y_stability) - np.min(y_stability))) - 0.5\n", | |
| "\n", | |
| "# train oracles\n", | |
| "train_oracle(X_activity, y_activity, \"catalytic\", show_preds=False)\n", | |
| "# time.sleep(1)\n", | |
| "# train_oracle(X_stability, y_stability, \"stability\", show_preds=False)" | |
| ], | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "display_data", | |
| "data": { | |
| "application/vnd.jupyter.widget-view+json": { | |
| "model_id": "6ad08c7883434b4ba5f8b212e2551ad1", | |
| "version_minor": 0, | |
| "version_major": 2 | |
| }, | |
| "text/plain": [ | |
| "HBox(children=(FloatProgress(value=0.0, max=200.0), HTML(value='')))" | |
| ] | |
| }, | |
| "metadata": { | |
| "tags": [] | |
| } | |
| }, | |
| { | |
| "output_type": "stream", | |
| "text": [ | |
| "\n" | |
| ], | |
| "name": "stdout" | |
| }, | |
| { | |
| "output_type": "display_data", | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd3hc1Zn/P++Myqh3WbbkXjFgGzC9hBpICIQUFgKpm01CNglks6QHQgokJJvyC0sKZEkhCRAgEIcOCaYX29gYbGPjKsu2itXrqMz5/fHeO3MlzUgjWaPiOZ/nmWdm7ty599wp7/e85ZwjxhgsFovFkrz4JroBFovFYplYrBBYLBZLkmOFwGKxWJIcKwQWi8WS5FghsFgsliTHCoHFYrEkOVYILJYxRER+LyLfH+V7rxSRJ8a6TRbLcFghsFiGQER2i8i5CTjuHBExIpLibjPG/NkY886xPpfFMhxWCCwWiyXJsUJgSRpEZKaI/E1E6kSkXkT+V0Tmi8i/nOcHReTPIpLv7H8nMAv4h4i0ichXnO33iki1iDSLyLMicmSM870pIhd5nqc65zgGeNbZ3OQc+2QR+biIPO/Z/0gReVJEGkSkRkS+kajPxpLcWCGwJAUi4gceAvYAc4By4G5AgB8AM4AjgJnADQDGmI8AlcBFxphsY8yPnMM9CiwESoHXgD/HOO0fgQ97nr8bOGCMWQ+c4WzLd4790oD25gBPAY85bVsA/HMUl26xDEvK8LtYLIcFJ6AG9cvGmF5nm9v73u7c14nIT4FvD3UgY8wd7mMRuQFoFJE8Y0zzgF3/BFwnIrnGmBbgI8Cdcbb3PUC1MeYnzvMu4JU432uxjAjrEViShZnAHo8IACAi00TkbhHZJyItqPEujnUQEfGLyA9FZIez/27npUHvMcbsB14APuCEm95FbO8hWnt3xLmvxXJIWCGwJAt7gVneKh2HmwADHG2MyUVDOeJ5feD0vFcA7wXOBfLQMBMD3uPlD84xLwVeMsbsi3HcaO2dN8w+FsuYYIXAkiy8ChwAfigiWSISEJFTgRygDWgWkXLgywPeV0N/g5wDBIF6IBMVkqF4EDgWuAbNGbjUASFiG/uHgOki8kURSReRHBE5cbiLtFhGgxUCS1JgjOkDLkKTrpVAFXAZ8B3UUDcDDwN/G/DWHwDfEpEmEbkWNeZ7gH3AZuDlYc7bCdwPzPUe2xjTAdwIvOAc+6QB72sFznPaXA28DZw14gu3WOJA7MI0FktiEZHrgUXGmA8Pu7PFMgHYqiGLJYGISCHwSbRiyGKZlNjQkMWSIETkU2jS91FjzLPD7W+xTBQ2NGSxWCxJjvUILBaLJcmZcjmC4uJiM2fOnIluhsVisUwp1q1bd9AYUxLttSknBHPmzGHt2rUT3QyLxWKZUojInliv2dCQxWKxJDlWCCwWiyXJSagQiMgFIrJVRLaLyNeivP4zEdng3LaJSFMi22OxWCyWwSQsR+DM/34rOky+ClgjIquMMZvdfYwx/+XZ/wvAMYlqj8ViSV56enqoqqqiq6tropuScAKBABUVFaSmpsb9nkQmi08AthtjdgKIyN3orI2bY+z/IYaZB95isVhGQ1VVFTk5OcyZMweRWBPFTn2MMdTX11NVVcXcuXPjfl8iQ0Pl6KhKlypn2yBEZDY6Kde/Yrz+aRFZKyJr6+rqxryhFovl8Karq4uioqLDWgQARISioqIRez6TJVl8OXCfM0PkIIwxtxljVhpjVpaURC2DtVgsliE53EXAZTTXmUgh2IeusuRS4WyLxuXAXQlsC2t2N/Cjx94iFLJTalgsFouXRArBGmChiMwVkTTU2K8auJOILAEKgJcGvjaWvL63iV+u3kFrsHf4nS0Wi2UMqa+vZ8WKFaxYsYKysjLKy8vDz7u7u4d879q1a7n66qsT2r6EJYuNMb0i8nngccAP3GGM2SQi3wXWGmNcUbgcuNskePa7/Mw0AJo6usnLiD+bbrFYLIdKUVERGzZsAOCGG24gOzuba6+9Nvx6b28vKSnRzfHKlStZuXJlQtuX0CkmjDGPAI8M2Hb9gOc3JLINLoVZavwb2ruZXZQ1Hqe0WCyWmHz84x8nEAiwfv16Tj31VC6//HKuueYaurq6yMjI4He/+x2LFy9m9erV/M///A8PPfQQN9xwA5WVlezcuZPKykq++MUvjom3MOXmGhotEY+gZ4JbYrFYJpLv/GMTm/e3jOkxl87I5dsXHTni91VVVfHiiy/i9/tpaWnhueeeIyUlhaeeeopvfOMb3H///YPe89Zbb/H000/T2trK4sWL+exnPzuiMQPRSBohKHCEoLFj6HicxWKxjBeXXnopfr8fgObmZj72sY/x9ttvIyL09ETvtF544YWkp6eTnp5OaWkpNTU1VFRUHFI7kkgIVDEbrUdgsSQ1o+m5J4qsrEiY+rrrruOss87igQceYPfu3Zx55plR35Oenh5+7Pf76e099AKYyTKOIOHkBlLxiSaLLRaLZbLR3NxMebmOuf39738/rudOGiHw+YS8jFQbGrJYLJOSr3zlK3z961/nmGOOGZNe/kiYcmsWr1y50ox2YZqzf7KaI8pyufXKY8e4VRaLZTKzZcsWjjjiiIluxrgR7XpFZJ0xJmodatJ4BKAJY+sRWCwWS3+STAhSbbLYYrFYBpBUQpCfmWaTxRaLxTKApBIC9QisEFgsFouXpBKC/Mw0unpCdHZHne3aYrFYkpKkEoLCLDu62GKxWAaSNCOLwTu6uJsZ+RkT3BqLxZIs1NfXc8455wBQXV2N3+/HXWTr1VdfJS0tbcj3r169mrS0NE455ZSEtC+phMBOPGexWCaC4aahHo7Vq1eTnZ2dMCFIqtCQnXjOYrFMFtatW8c73vEOjjvuOM4//3wOHDgAwC9+8QuWLl3KsmXLuPzyy9m9eze//vWv+dnPfsaKFSt47rnnxrwtSeUR2InnLBYLj34Nqt8Y22OWHQ3v+mHcuxtj+MIXvsDf//53SkpKuOeee/jmN7/JHXfcwQ9/+EN27dpFeno6TU1N5Ofnc9VVV43YixgJSSUE4dBQu/UILBbLxBEMBnnzzTc577zzAOjr62P69OkALFu2jCuvvJJLLrmESy65ZFzak1RCkJbiIyc9hYNtwYluisVimShG0HNPFMYYjjzySF56afBS7Q8//DDPPvss//jHP7jxxht5440x9l6ikFQ5AoCF07LZcqB1opthsViSmPT0dOrq6sJC0NPTw6ZNmwiFQuzdu5ezzjqLm2++mebmZtra2sjJyaG1NXF2K+mEYFlFPm/ub6YvNLVmXbVYLIcPPp+P++67j69+9assX76cFStW8OKLL9LX18eHP/xhjj76aI455hiuvvpq8vPzueiii3jggQdssnisWFaRx+9f3M2OujYWTcuZ6OZYLJYk44Ybbgg/fvbZZwe9/vzzzw/atmjRIjZu3JiwNiWhR5AHwOt7mya4JRaLxTI5SKgQiMgFIrJVRLaLyNdi7PNvIrJZRDaJyF8S2R6AecXZZKX5eWNfc6JPZbFYLFOChIWGRMQP3AqcB1QBa0RklTFms2efhcDXgVONMY0iUpqo9rj4fMJR5Xm8XmWFwGJJJowxiMhENyPhjGbVyUR6BCcA240xO40x3cDdwHsH7PMp4FZjTCOAMaY2ge0Js3xmPlv2t9DdGxqP01kslgkmEAhQX18/KiM5lTDGUF9fTyAQGNH7EpksLgf2ep5XAScO2GcRgIi8APiBG4wxjw08kIh8Gvg0wKxZsw65YQtKsunuC1HT0sXMwsxDPp7FYpncVFRUUFVVRV1d3UQ3JeEEAgEqKipG9J6JrhpKARYCZwIVwLMicrQxpl8m1xhzG3Ab6OL1h3rS3Ay97JYuO9WExZIMpKamMnfu3IluxqQlkaGhfcBMz/MKZ5uXKmCVMabHGLML2IYKQ0LJDeicQy2dvYk+lcVisUx6EikEa4CFIjJXRNKAy4FVA/Z5EPUGEJFiNFS0M4FtAiDHEYJW6xFYLBZL4oTAGNMLfB54HNgC/NUYs0lEvisiFzu7PQ7Ui8hm4Gngy8aY+oQ06I374I53QajPExqyHoHFYrEkNEdgjHkEeGTAtus9jw3wJeeWWDoaoPJF6GggJ6CDyqxHYLFYLMk0sjirWO/b68gJOB6BzRFYLBZLMgmBrg9Kex2pfh8ZqX7rEVgsFgtJKgSgJaS2fNRisViSUggOAlo51GqTxRaLxZJEQpBRAOKDDhWC3ID1CCwWiwWSSQh8Psgs9oSGrEdgsVgskExCABoe8oSGWjqtR2CxWCxJJgQejyCQYj0Ci8ViIYmFICeQSktXz2E/La3FYrEMR5IJQSQ0lJuRQk+fIWjXJLBYLElOkglBMQRboKcrPPGczRNYLJZkJ8mEwBlL0HGQ3ICdeM5isVggWYWgvS6yJoEdS2CxWJKcJBWC+vBU1LZyyGKxJDtJJgTeGUhtjsBisVgg6YRgcGjIegQWiyXZSS4hSMuGlAC010bWJLA5AovFkuQMKQQi4hORU8arMQlHBLKnQWs1mWl+/D6xaxJYLJakZ0ghMMaEgFvHqS3jQ14FNO9DRMgJpNhVyiwWS9ITT2jonyLyARGRhLdmPMgth5YqAPIzUmno6J7gBlksFsvEEo8QfAa4F+gWkRYRaRWRlgS3K3HklUPLAQiFmJGfwf6mzolukcVisUwowwqBMSbHGOMzxqQaY3Kd57nxHFxELhCRrSKyXUS+FuX1j4tInYhscG7/MZqLGBG55RDqgfZaKgoy2NdohcBisSQ3KfHsJCIXA2c4T1cbYx6K4z1+NL9wHlAFrBGRVcaYzQN2vccY8/kRtPnQyKvQ++Z9lOfnUtsapKunj0Cqf9yaYLFYLJOJYT0CEfkhcA2w2bldIyI/iOPYJwDbjTE7jTHdwN3Aew+lsWNCbrnet1RRXpABwIHmrglskMVisUws8eQI3g2cZ4y5wxhzB3ABcGEc7ysH9nqeVznbBvIBEdkoIveJyMxoBxKRT4vIWhFZW1dXF8eph6CfR6BCYMNDFoslmYl3QFm+53HeGJ7/H8AcY8wy4EngD9F2MsbcZoxZaYxZWVJScmhnzCiAlAxo2UeF4xHsa+o4tGNaLBbLFCaeHMFNwHoReRoQNFcwKPEbhX2At4df4WwLY4yp9zz9LfCjOI57aIhA7gxorqIsL4BPrEdgsViSmyGFQER8QAg4CTje2fxVY0x1HMdeAywUkbmoAFwOXDHg+NONMQecpxcDW0bQ9tGTVw4t+0j1+yjLDVBlS0gtFksSM6QQGGNCIvIVY8xfgVUjObAxpldEPg88DviBO4wxm0Tku8BaY8wq4GqnIqkXaAA+PpqLGDG5FbBzNQDlBRlUWY/AYrEkMfGEhp4SkWuBe4B2d6MxpmG4NxpjHgEeGbDtes/jrwNfj7u1Y0VeObRVQ18v5fkZrNndOO5NsFgslslCPEJwmXP/Oc82A8wb++aME7nlYELQVk15QQb/2HiA3r4QKf7kmozVYrFYIL4cwdeMMfeMU3vGh5wyvW+roTy/hL6QoaY1GC4ntVgslmQintlHvzxObRk/skv1vq2WmYVq/Pc22BJSi8WSnMQTC3lKRK4VkZkiUujeEt6yRJLteASt1cwqzASg0gqBxWJJUpIzR+AuWdlWy4z8DHwCVVYILBZLkjKsEBhj5o5HQ8aVlDTIKIS2GlL9PmbkZ1iPwGKxJC0xQ0Mi8hXP40sHvHZTIhs1LuSUQVsNALMKM60QWCyWpGWoHMHlnscDa/0vSEBbxpfs0rAQzCzIpLLBDiqzWCzJyVBCIDEeR3s+9cj2eARFmRxsC9LRPfL1ix9cv4/aVjuNtcVimboMJQQmxuNoz6ce2aXQWgPGMNOpHBrpVBPNnT188Z4N3LeuKhEttFgslnFhKCFY7q5RDCxzHrvPjx6n9iWO7GnQF4Su5kgJaf3I8gRNzsL3zZ09Y968w46eLvjteVC1dqJbYrFYBhBTCIwxfs8axSnOY/d56ng2MiGERxfXjnosQUtnr3NvhWBY2mqg6lXYv36iW2KxWAaQvJPrhEcXV1OQmUp2esqIhcD1BFxBsAxBb9C5t/kUi2WykcRCME3v22oREWaOooQ0LARd1iMYlj4rBBbLZMUKgVM5NK8ki511bSM6RMQjsEIwLGGPoHti22GxWAaRvEIQyAN/OrTqYmvzi7OobOgg2NsX9yFcIbDJ4jiwoSGLZdISc4oJpzooZpmoMSY3IS0aL0Q0YdyqK2XOL80mZGBPfQeLpuXEdYhIaMjmCIbFFQBXECwWy6QhphAYY3IAROR7wAHgTnQg2ZXA9HFpXaIpmA2NuwGYX5INwM66tpELQWcPxhhEpv44u4ThCkCfFQKLZbIRT2joYmPML40xrcaYFmPMr4D3Jrph40LBXGjYBcDc4iwAdtS1D/WOfri5gd6QoaM7/pBSUhJOFlshsFgmG/EIQbuIXCkifhHxiciVeNYuntIUzoWOg9DVQlZ6CtPzAuyojT9h7M0N2MqhYTgccwT711thsxwWxCMEVwD/BtQ4t0udbVOfQmdJhUb1CuaXZLMjRuXQhr1N9IX6p0z6CUGixxJseQja6xN7jkTSe5h5BO31cPvZ8ObfJrolFsshM6wQGGN2G2Pea4wpNsaUGGMuMcbsjufgInKBiGwVke0i8rUh9vuAiBgRWTmCth86Bc5SCw2uEGSxs64dY/ob/B11bVxy6ws8ubmm3/bmzh6Ks9OABHsEwTa450rY8OfEnSPRHG5C0NkIJgQdU1icLRaHYYVARBaJyD9F5E3n+TIR+VYc7/MDtwLvApYCHxKRpVH2ywGuAV4ZaeMPmUJHCByPYF5JNq3BXupa+xur3Qc1ElbV2H/AWXNnDxUFOj1Fc0cChaDbicR1j2ycw6TicMsR9DjfSa+dvtwy9YknNHQ7uh5BD4AxZiP91yqIxQnAdmPMTmNMN3A30ZPM3wNuBsY/eJyeA5nF0LATgIXTtHJo04GWfrvtb9I/e11bxIiFQoaWrp7wPEUJ9Qh6OvrfT0XC5aOHSY7AFeceKwSWqU88QpBpjHl1wLZ4AuLlwF7P8ypnWxgRORaYaYx5eKgDicinRWStiKytq6uL49QjoHBeODS0YmY+KT7hlZ0N/XapcoXA4ym0BnsxBmYWZgAJHl3sGs+pbHQOt9BQtyvOh4mwWZKaeITgoIjMxxlcJiIfRMcVHBIi4gN+Cvz3cPsaY24zxqw0xqwsKSk51FP3p3BueCxBZloKyyryeGVX/7jv/ib9s3uFwDX84dBQIpPFYY/gMBCCqTCOINgGT90w9OfthoamspdmsTjEIwSfA34DLBGRfcAXgavieN8+YKbneYWzzSUHOApYLSK7gZOAVROSMG6uChuqk+YV8UZVM+3BiGHfP8AjeOzN6vC2wqw0stL8CQ4NOQZpKhudqVQ+uutZeP5nsOu52Pt0HwbibLE4DCkETsL3P40x5wIlwBJjzGnGmD1xHHsNsFBE5opIGppXWOW+aIxpdiqR5hhj5gAvo4PXxnflksJ5gAmHh06cV0RvyLBuT2N4l33OymUH24JU1ndw1Z/W8bOntgGQG0glNyM1saGhnsMgNDSVksVdzXrfsCP2PjZZnBi6mmHn6oluRdIxpBAYY/qA05zH7caY1ngPbIzpBT4PPA5sAf5qjNkkIt8VkYsPoc1jy7Qj9b76DQBWzi7A75NweKinL0RNaxepfqG+vZutNfoRvLyzATAs2PQLjkyrscni4ZhKOYKuJr2vjyIEtVsgFLIeQaJY/ye4830QjNvUWMaAeEJD60VklYh8RETe797iObgx5hFjzCJjzHxjzI3OtuuNMaui7HvmuHsDACWLdRbS6tcByErXPMGLO1QIqpu7MAaOmJ6LMbB2TySRXEArJa/9nHPMy4mdgTQcGophdIyZ/EnLKSUEMTyC+h3wy5Pg7ScOj7zNZKSjQcdnWCEYV+IRggBQD5wNXOTc3pPIRo0r/lSYthQObAxvOm1BMRurmmnu7GGfkwtYXpEPwJpdDfh9Orlcvmh4IN/fldiRxb3DCMG2x+DH8yMGbDLizRGYmJPaji2hPrjrCtj9wsje1xnDIzj4tt637rflo4nCHSvTfXjMYjNViGdk8Sei3P59PBo3bpQtgwOvhw3U6QtL6AsZXtpRH04KL5+pQvDGvmaOmpHLrMJMCn36Y82VzolNFtdv1z9Qa0301ycD4SSxgb5xmpepswm2Pgy7nx/Z+1xBbd7b34Npdqqhu1qsR5Aogq4QTOHBk1OQeEYWB0TkcyLySxG5w72NR+PGjenLNC7s/NGPmZVPVpqf596uiwhBRR4APX2GucVZvP/Yco4s0BlHc6UzwaGhYYxOlzMAbjK7032elcnGq4Q02NL/Pl5cITChcGkxAE17IsfrHoe8zZ6XYMNfEnf8yUi38xu2HsG4Ek9o6E6gDDgfeAYtA53EFmcUTF+h9054KNXv4+T5RTz39kH21HdQlJUWHi8AMKc4i2vOWch3ztPxcdnSSWtXL109CZqKOlw1FMPouAIQnIShoYZd6ml5y0bHK0/g9ipHKpBdTZCq05L3Cw81VTqvt3iqhhKYm1nzW/jn9xJ3/MlIeDqVKVwYMQWJRwgWGGOuA9qNMX8ALgROTGyzxpnSpSA+DQ85nL6wBGncyb3r9rKsIo+MND856bqOz9ziLEQEcapLctAf7cC5iOKhq6ePGx/eTOtQoaXhPALX0HWNsOebaJoq4ZZjYfs/+xv/8RpLEBbIkQpBM8xwOgcNUYQg2Do+HkF329TpGTfuCS/7ekjY0NCEEI8QuBaqSUSOAvKA0sQ1aQJIy4TixXBgQ3jTJf7neSb9S/zspA7+94pjASjJSQcii9jQqWMNMo0agz31IzcKr+1p5PbndvHSjiFmsXQFoK8b+qIkpV1PYLKFhpr3aXilpWqAEIxXaGiUQtDZBPmzIaNQ8y8uTU6OIOjNESRQ1IJtGioZr+T6ofC3T8GjXxm8/bGvw18ui/84Nlk8IcQjBLeJSAFwHTogbDPwo4S2aiIoPw6q1uqfrv0geauvA+B9FW1kOZ5AcU46x8jbLH30g7ooiSME6X2jF4KGDo2dD7nucb+wShSvIGzwJplH0OmU2gZbNS+QEtDnk10IupohI7/fPFR0t+siRqCel2uo+oJanZQIultVSKfCaOy22ujFCtVvwM5ndOxFPARtjmAiiKdq6LfGmEZjzDPGmHnGmFJjzK/Ho3HjSsVKNVyNu3SemWAriD/SCwTO9b3GPenfI2XfGtjxdFgIfN0tZKenUNkwciFobHeEYKhkszf8EC085IaEJltoyPl8CLaq8Q9owj1uw7b+TyMv/fQSHEUSva9XDXAgD7KKI4PLPL8Dgs3Dfycumx6E2rfiP7+X4BTqHXe3RS9f7mrSzkvT7viP4723jAsxF693EZHro203xnx37JszgVQcr/e7noU374cVV+hQ9+aIAfho39/oy50JwXpo2R+uN5dgKzMLMthTP/I/bKOzjsGQ5adeQxPNKIy255toOhyPoKtFjX9mMbTVxOcRGAOPfQMWnA1zTo3vfBvvhdpNcO4N+nw0npK7byBPpymv26rP3d9BwRy9Hq8X0NsF6dnRj7fqC3DU++Gi/xd/G1y8ye6s4pG/fzwJtoEvijnpdMSh9q3IioDDHQemhvgdRsS1ZrHn1ocuNDMngW2aGEqP0EqR536ivb2j3g95M3VCOoCuZgI168lY/gHIq4DWA5Eer+ljQYGfPaPwCBocj2DI8lOvEETrfYZ7vpOsaqifR9ANgVx9Ho9H0F6n1zMSL2fTA/DanZHnoxEC1wMI5KsQuMdwS0enHeXkCNoh1akki5UwDoV039EO9JsqRjHUp73+qB6BKwSbhz9OX0+ktHgqT6cyBYknNPQTz+1G4EwgDmmfYvj8UH6sVoZkFsHs0yB/ZiQksOs5MH0w/2zInQEt+yKGDliYF6KqoXPQusbD0ejmCIYamTysEEzSqqFwjsDxCNIdIfCOKYjFwW2R98ZLW42e041HBz096ngTru6oYtcjCAtBJfjToGhBpGoos0hfixUacmviRyMEoVCkRHWyC4HrufR0qOC7hPoinZPaLcMfx+vR2tDQuBKPRzCQTHQsweGHGx5a8h7wp6hH0Lpfeyo7n1aPoeJ4yJnuhIYaw8ZtTk4f3X0hqltGlthzPYJhQ0OuER3YU+rrjWybdMnixsi96RtZjsCdzmEk4tZWq8lVt1fvGhYTir+H6RrtDMcj6AtqKKtpr3qCGfkQ6tXtwwmBe35XXIaiZT/87wmRcQteQ9g9yUJ+Awl62ur9DXoFMB4h6HfNk1z8DjPiGVn8hohsdG6bgK3AzxPftAlgthOLPvqDep9X4ZQ/7tfk8JzTICUNcsvV6HQ2apkhMCtTY8YjzRM0OTmCIUNDvZ2QUaCPBxod7x9vsuYI3AXew6GhOHIEbtlmvOJmjHoE3vOO5rNxjVcgD9Id4Qq2acVQVmlEkCEStx92fEccHsHuF+Dg1shYFq8hnOxG0ds+r+i51509TT284aYWCVohmCji8QjeQ2SyuXcCM4wx/5vQVk0UC86Bz74Ec8/Q5/nOujp7XtCBRfPP0ue50wGjvVxnnxkZ+iPfGyNPcLAtyE2PbKGnr38ZXVtbK19P+TO9HUMYi55OT+9zwPG9Bm7ShYYco9DulF2mj0AIwh5BnGGVruZIfNlbtuoStxAMCA2BCkpno4qx69VA7O9k4DnjuYbaTf3bPpV6x16PxXut7mc56yQI9YTXBo99nCl0zYcZ8QhBq+fWCeSKSKF7S2jrxhsRnYnUJW+W3r/o6N7id+l9rmfp5XzdpyglSKpf2HnQ+QFvuCu8xgFtdTy1YSe3PbuTzfv7G+t5nW/wmZSHWdqxJna7eoaIR7u93ozCyRsacuvvAzpxX1yhoXpHCHo64pukrt2zlrXrgfQTgjg/m7BHkD9ACJoj4SKXTMcjiHU94bLe5uFzFG7opMOTYA+3fRzi5aEQPH1T/zLZePG2r8vjEbgdgVkn6/1w4SH3mgN5NkcwzsQjBK8BdcA24G3n8TrnNv7rB4wneY7Br92kA84K5ujznBmgbvQAACAASURBVOmRfZzQkL+nlQWlObx1oFWTZP+4Gl52hlv84SIWbbwZUM/Apaunj7xeNVq53bVqLNb+bvAfv6cTMh3NjeUR5JUfemgoFNIxFPHEc+PB7d0axwuKNzTUG9QpC8KhmTiuq80zmMkrBFkl8R8D1HiJH9KyPELQGvEIvKGh8HcSKzTkCEGoZ/hZSmucqpqoHsE4GMWGHfDMzfDmfZFtob74RMjbe4/mEUxfrvfNw4iMe53Z0/rPNbTlIdi3bvh2WEZNPELwJHCRs6xkERoqesIYM9cYc/hVD3lJzYgYkiM9a/Hkzog8djwCgq0cMT2HLQdatLS0rxuaKzWZe3Ab2S3qFnuFoKmjh1LRP0tRXx29VevgoS/C63dFjh9yRpbG8gjcXmduRfy951g07tK1ejf/ffTHcOnp0vb40yLb0uMsH23YpWG3cp3aI67QSiwhcL23keQIMvLVO3SFoKPBGWSWHxEziD9HMNw1dLXob8U9F4x/vNydR8k72+qa3+pcUcN5M16h6icEzuP82Vpo0byPIQl6hcCd1C+o01c8/N/DXoJl9MQjBCcZYx5xnxhjHgVOSVyTJhl5Tp7gyPdFtmUURKZLcPMIXS0snZ5LbWuQ5v1OWKO5SquOTB853Wqo6lojQtDQ3s000VDAdKknWOtUjOx9JXIu12jG6n16PQLv89EQDmXVjv4YLm5YyP38IDLoarjyUTcsVLFS7+MJ63jb7ApBd9vohMDNA7j37liSQR7BcEIQo4JmIHWekccT5RG4vXWvENS9pQI7nBD3E7wooaGMfP19tlQNfRzX+OeURR7veUE7FPvXw8Htsd9rOSTiEYL9IvItEZnj3L4J7E90wyYN885UbyDPkxcQiXgFWSU6sCjYwpIyNRI1e9zRqFUa4gCK+g4ihPoJQWNHNyWORzBD6umuc4Xg1ci5XCOTnqsjNweFhpw/qWvw2ut0YFWsXlxvt1am7Hxm8AR2rhB44+2jxRUC12MCSMlQAR3OI3CNkRtSiJYENwa2PREZM9BWA75U7U12OGMJgq2R72k4IQi2aSivcXdEAFyPwO0tZwzwCOJNFsPQxrTGSRQXzvd4BM57fSnjIwRNUYTATfJ3DDEhIgwdGvKl6P8jryIiqDGP41xz9jS9ZmPg7Scdr1Lgjb/GcyWj5437oPLlxJ5jkhKPEHwIKAEecG6lzrbk4Nxvw6W/G7w9xzEwbi/RCQ0BtFU7PZe+7nBsM016KaKVg22R3nBDe3c4NDRD6jENu/WFpj1aYfG3T2tPCDRMlZoVEYaGXfDCLyJG0u15r/kt3PvxiHEZyD+ugd+cAX+8GN56qP9rNW/q/ZgIgWPQCmZHtqWk6224HEHzPkjLiVxTNI9g76vwl0th26P6vK0WsktVmDsanMFYxqnwinEML1tWwWNfhX1rI0ltVwjc3nJGgbYLXao0LASxhC1eIajdAmnZKnwDPYKs0sSEhl69Hf7pmSXGvcamvZEOgitKnoGTUXHbmlE4ODQUcMJsueXxhYZ8qer9mj79nbz9BMx9B8w9Hd64N7EzsT7+TXjp8CyIHI54RhY3GGOuMcYcg65b/EVjTMNw7zvsyZ2hvdvUDGcEagtF2emU5qRHDDpA5UvhhzN89f08gqaObqahf7Jimklp2KYGAeDB/4SN98Cmv+nz1Ew9l9v7fOU38OR1arx9KZDt5DJ2Pav3zXu1V7xzdf8/T90WXZoTBpfzVY+lEETzCNLBnz68R9C8Vz0wt/cdzYi6bQ+Hs2pUCDIKtAfrGuGMQv2ehiutbdwDCJz0OVjprMSaEtDP1p1eIpAPPl9EINKz9XpieQRdcYSGjNHwR+lSFRbX+LrGP7skMVVDmx6A9X+OPHc9AtMXCeF0xOkRBNv0d5tR0P86O5vUiwL1CNprh+4EdLdpkt79D9S8qeNJFr4Tjv43/c63PR7/NY6EUJ+2r32Yaz1MiSkEInK9iCxxHqeLyL+A7UCNiJw7Xg2ctBz1fjj+P/SxZyqCI6bnEmivivQWPUJwXH5Hv2RxQ5t6BL3p+fjEkF3/pv7o/WmR97kVPK7ouB6B+/quZ/X8buzajTc3V8GOf8If3wtVntLU5n0w4xg1kF5XvaNBDYAvBdrGQAhcg5Y/0CMI9J+GIBot+7QH6VYNRTPibg/W9WLaarX3nFnUXwjSc/pPFRGLpkoV9wtugqUX6zYR/VzDoaGCyDHBI85D5AgynNxOV4zRxW89rNdw3Me1J9zVHKnWScvW8yfCI2irgbbqyHoKzXsj4UU3PBR3aKjVEYL8AQPKmiJhNvfYLUNElYNt+tm6czhtdVKTC86BZZdB8SJ49MuJWb2so16r28aiEzQFGcojuAwdRQzwMWffUuAdwE3xHFxELhCRrSKyXUS+FuX1q5yRyxtE5HkRWRrtOJOSxe+C82/Ux4HcsLFaMj2Hkt4DdEw/SV/raqY5oH+CY/Pb+3kEna0NZEg3vWXHAOAzPVCyWA01qOFxZ79MzdA/SE+HGrVqXVaT9jo1Ft6BTqDG1B2U5Y7S7e3WXk9u+eCYrWtQZ56oeYdDXTMg7BF4hMDvhoaG8wiqtH2uRxAtrOMaZzcE5oaGMos0vBIWgtz4hcCb2HZJz+k/7YR7TNDea2rm0FVDbm4pmkfg1u4XLVBDl1EIGN3XNa5pWYnJEbhrBzRXaSioZT/MOV23NezStrlhqo562P4U/Pbc6CLe3a7tDOQN9gjcMJv7ObQMER7yXjPAvtf08y2YqyP6L/yJfk8v/mL01x2L8Kj0g2N/7CnAUELQbUw4pnA+cJcxps8Ys4X4pq/2A7eis5UuBT4UxdD/xRhztDFmBbrYzU9HfAWTAY+huXx5MaXSxJ178gj69Qe93b+ALtKo8DXQGoysbdzXokv7+dzqGNAf/YlXwWlf0gnu3InHvL3PqjXae3FLM9Nz+1ezgP6x3ZCGk7Cm1emN5ZX3n1nVmEhIyR09fag9o84GNfzZnsXs4skR9Ab13HkV4E/VBHM0I+oKQcMu/ezb6zTJmFmkIuT2wNOz4xcCbxjLxfu5ukYtkKuxbH8qpAaG9giySpzQVJRr2POCjlE54ys6t5VbGdbR4PSOs9UwjrUQdLdHErNNe8KVbcw6Ua+rcbd+fu74j456/X1UrYle+eO2daAQuKW4EBHZofIE3msGXTGwaL6G40BH/M85HbY9NupLj4krBJ2Nh1aCPUUZSgiCInKUiJQAZwFPeF7LjPEeLycA240xO40x3cDdwHu9OxhjvF29LGAKrMkXhfQ87Wnf/ynmbv8DALt6i9nVo3/sNc05NPqLKepT4+p6Bb42FYLUWV4hmKNhp3O/3b83neoJDVW+rGssL73EOX9O/xGvGQX6h3NdfFcQXLc8d0bEI+jphNvPhmd/rMnK0iN1n0MWAndKBo8hTQmoEPQNIQRujzHPmdcwkBvdI2je6xgMo/P0mD5HCArVgLki54bNhhKCvl49b1QhcD7XtBw11qDHS3P+AsN5BK63Fk0IXE/NXW/BDSN1NjjxctcjiBEaeuwbOivuSPGOuWjaE8kPFMzRz6BxdyQsBCpMrgGPNvK4u00/n0FC4PEIwqGhISqHvNcM+hsqXtR/n7JlurbBWK8KF638OIkYSgiuAe4D3gJ+ZozZBSAi7wbWx3HscsD7q6lytvVDRD4nIjtQj+DqaAcSkU+LyFoRWVtXNwljeOk5+uN546/wr+8BcMPHLmTBwiMAqE8pI5g5nbwe/bHVtQXp7g3RXKcfjxQtoBnnx184N3JcdyQzOB6BExqqfAnKjtZJ8ECNZWpAPYSUgPacWvZFPAH33v0zu6GhYLMmk/e/Bmd9Ez7xWGQA3aHmCdpq1Sh7e9ThHMEQQuAacNdwBPIG5whCjqGff7Y+f/EWvS+cF8nNuCIYT47A7REPJQRufsBtk9trTQlEXz4UHCGIYiBdmio1J+OOVHfP4XoEaa5HEEUImvfBy7fC498YeSWNd0nJpspIviVvlv7mGnf3N4Yd9RGBjlYC6iZ5A/mOJ2H01unJEaRlOh2UIYQg7BFkRbYVLey/T+kR+nl7y1zHgtbqyOMkzBPEFAJjzCvGmCXGmCJjzPc82x8xxoxZ+agx5lZjzHzgq8C3YuxzmzFmpTFmZUlJyVideuwoXaK90Q/do6EdIFAyj5RCNSxfuexcZs5ZSEan/tgOtgZ5ZVd9eHoJcqZR5yshKAHIKqG5s4cDzZ39Sy9TM/TWfhD2rtH5W9xcQriKJVerT/Jn9Q8NuWEU98/segQAW/6h98d9XP+s7mjZQ/kz9PXAnhd1QJjPr2Wv4AkNDcgR7H4+0gNtHuARpOcONqJtNVqaO+d0Pfae56HiBBUGN7wSFgI3RzBE1ZD7+UQTAtejyfDkYE68Cs5zSi+HShZ3tQztETRV6nX6/Po80xGCTmckc3q23rrbBq/5u8+Z3aV6Y7+ChLhwPQLxaRvcXn5ehQpBw65IrDw9r79HEM2Qe0NDfd36/Xa3q7i6oSHQ0e9DhYZcz8IVWYDiAULgzgUWqzx6tHg9gvYR5gna6oYX46a9k292YA+jWY8gXvYB3uxbhbMtFncDlySwPYnjuI/Dl96CxRfARx/UZQmzS8PGLLVoDv78ClI7avARoq4tyD+31DLd34RJ0xj2gZSZ7E2dCyJ8/6HNfOT/Xh1QceMki1uqtEd09KXaO0rJiIQU5pwGR1ykf7i+YGSyupZ9muRr2ad/7HRPjf7WR3R/N5bv3rcfwujiva+o4V34Tn3uCpU/Svno/g3w+/do0hQ8HoEzTiNaaCi8bORs/QwQeNfNGksOewSOCKbFkSMYSgiieQQzj49MVR4rNBTq0/xOINfpKccSAs9fJGNAjsAbJhlYolq1Rj3AQD68/KvY1xYNVwhKl+rn1LBDK65SA1CyRD1Ft5S4eKF2Ctz8UrT5gtyQjtv772ruv9KbS1750MniYOtgj2CgEJQ433c8K56NhLaayFKbIxGCXc/BTxZHvNJoGKOJ9n99P75jjnXYKw4SKQRrgIUiMldE0oDLgVXeHUTE+y1fiE5qNzVxE1oFc1QYAJZcBMd+VKtC8soR08dPU3/J0tdvYtHGmzkzsAPJngbAn4qv4cacbwKwsaqZ3QfbCeWU6wRoEPEIAKav0Enw/Knw4fvhVI2ofS/zq9ybcWn/uZDmnA4Y/QO37I9Ub7g97s5GmLEisr9bCTPSXhFoz37XszoIyJeqA4FAjan4Ncaekh6pPDFGQxsY2P6kPm+p0qkb3GtNzx0cGvIa7lOvURFw5yVyp32ofsOpUkqL9Mh7YlQreXvEA3GFwGvQvMRKFnvLV2N5BM17+4t9IE8/p84G7VF7jeLA8FDVWv0dHPcxHRjoluuCjhxvOeC5vkrY+mik19pWo+eZsUJr87c+CvOc76pUw5nsfl7vixfpPiFnkFnU0FB7pHwU9Fq9K725FM7XhXcGjmgH9SLdnEqaJwVZtKD/fmmZGj5NhEfg5iPi9YY7G+GBz6jn88LPY+dy2uu0VHd/HBH1+h1w43StmBpHEiYExphe4PPA48AW4K/GmE0i8l0RcQq1+byIbBKRDcCX0DLVw4fiBXDxLWqw550Fc07nZP9bLNz/Dz7Q8whzurdpVQSQW1DK600BuntD7KhrozdkONgZihhurxAc/x9a4w6aaHQM2APr9/H01tr+02G4ays07dE/sSsS2dPUWEMkxOSSVdLfVY4nBn3wbfjTB+HO98GGv8DskyNhlfScyNxM3ikmtqzSypmK49VY1W/X0IG3/dE8AlcI8mZqzf+Jn4m8llcB59+kQrnwPOczOl3/rG8PGIwUCukfr6lS4/Qp6YOvK+wRxBKCGB6B2+ZYQtAb1MkJvV6IiDMgzpssds7vrRzq61GjUnG8rqZnQpGqr72vwu3nwJ8v1Z7lfZ+Enx8Nd12u40pAcwTZpRrG7GzQ3vuyy/Q1VwiqXtXz587Qzw6cieMGCEEopG1Lz454TW01kQF/XnGdsUK92YNbGcTBbXoed+1wUE/V6x24lC5NgEdQrd6HLyV+IVh9s17ru36keZS1d0Tfz112tXbL8P+lva+oN7/nxfjbPgbEJQQicoqIXCEiH3Vv8bzPyScsMsbMd9Y7xhhzvTFmlfP4GmPMkcaYFcaYs4wxYyzzk4jCufDxh/hEwR85qut2Plr2d9o+/Qq8/zYAls/Mp6G9m+e319HrrHt8oLlLPQx/msaRS5bon/eoDww6fF/I0NTRTWtXb//1EsJCUKkegfuazxcxuNGEwP0zrP8z3HJc/wTjvtciveiOBtjxL50hMjWgyb32ukhYCBwhcEpd3fLRni544jqtUnrfb/S1t590xhB4J6qL4hE073Wm9shmECJw8ufgk4/D5X+OfAZZpTpFgZfX/qCza755f/SwkHt+6B8a8pKaoUZ/YAzfO47BFQKvEXAN6sDzZhZqfL6no39oyCsENW+qmFashBnHqljsXK3f0T0f0TbVvAF3X6HTSp/wGf0N7XxG399WrR0B1xvJLNaOCmiOKKtUY/2ZRZFQG0DFcdpu73W45c1p2ZER61VrI6GrsqMj+053PM9oPWN3hHjZ0Y7nGNCOVDSmHalCM9zU3rHY+2p/Dwq045MzXT+LeMcSbH1Ef+cnfka93+d/Ft1j8q6/Pdx8S67AjbXQDUM8S1XeCfwPcBpwvHNbOeSbLDH5+eUrePBzp3LPZ08le8aSsIFZXqE9zr+uifxQDjR3aiWM2ys99iNw9fr+rrNDS2cPIQNtwV79I/tSILtM3+9L0R67O5jMxTW4A4Ugu1SNebAVnrxeY8hPXq+v1WyGO86H370LDmyEX56kXsD+DXDRL+DD96nHstxTT+D1CLKKtR13X6FeygU3qVdUtFAHCtW9pYPqXAJ52ov01nYffDv64K9Y+Pwqntue6D/ydesjakT7goNj0eG2O0IQKzQ06xTtUb/2h8g2YwaHhkI9Gs5w479uIj9/wHVkFEZENlZoyJ2UsOJ4NZpzToNdz8DTN2oP/xOPqNHd9ph6Qxf8UPd1vYa2GkcIHBE6+oOR0liIeAVZxZHkO8DMk/S78FYUue1Kz9b9ixdpr7ZqjZYje72sogUqGPs3DP4cD2zUfJcbCsqfBeUxzMy0I9ULunkOPHytbqvfoccYjo4G/e164/XdHWqk3bmq4gmLNuzS79AV0Hf/WEOed18xeOSzd9bU4Qy8uy7FWIe+hmHYgWGo0V/qGVxmOQQWTcuJun1xWQ5pKT6e2hLpeR9o7tLBRq7bDpGQ0AAaOjTu3tbVq739nBnq1vv86p67c7R4wy7TjtQ/hvfPDvqH3vuKLlTScRAWvxs23q2CseHPapzaauC2MzU08qF79FiuUbvwJ/2PV7IkMtL41Gs0Efn243rceWfq9gXnwiu/gkXvgtOvjbzXNcRP36QCkTMddj+n5a4j4ehL9fgb/gIn/6d6JLue0xzOSVfFNvTRksX9jvtBFYEnvw2LLlBD+cdLImXA6blqdAF+faoatw/fFztBnVcBmx/Ux275KETmGwqFNARRckQk7DLvHTr5XuMeWPkJ7VW/+8c6idrFv9Dfw5zT9fvsbFTPYcYxaqhXXAknfbZ/G0qXqrB4PYKUQP8FZtzqMrddbjtnnaTrWfQGI3M2ufh86jUc2KB5gs5GPY6IVj9NWxqpoPrU09FDdaC98PO+C289ot/n+TfCg59VAf3S5sh/pG6bdoTc7+S0L6oXEeqN5KTuviLyHWRPg6yioUND//yu9vBdAXAHYJYshg/+H/zlMvUMzvb8Pg9uU++raY8KwaLzYx/fnVKmzhkr4X4eCSYeIXgTKAMODLejZfSkpfg4ckYu6yubmF+Sxd7GTqqbuyBvbn/jHYPGdkcIgk4i7uT/jPyJy1dGVp5ye3ugf6ZoawMUzNFe34u3aAz6/bfDb8/R2TkBLvuTutJPfAsu/T0sHGbqqbO/FflzBvLgQ3fD5gciyWSAM66FsqNg2eX9e6dunuF5Z9B5ep6275QvDH3OgZQfq8bwiW+qUc+ZpgZiwTlqLGIxXI5ABN7zczXyt5+lxqtpj2eiulz1RjKLNAfy5PXwh4t1Kg/xR2axdXnn9/Wz3/m0zsvvhr/c0NDmB9VIfPCOyGc670y99/nhtP/SxzNPgP94MnLcuWfAMz9U8es4qN5iagAu+eXga3J/I5kejyC3PCL0r98NL/w/zX+5I5TDQnAyvPZHfVwRpUc/Y4UK2f3/roKRVQLvvVVDQ941P6KF/VxSM7RDUbwY7roMNq9yvCSjgjJ9uXodt52pyfScGVpu+8zNkXEJTZWw/k+R+YzAEYKS2GMUutt1ssfuNj1fbnn/ZPai87UzsO73cMaXI+HQg9vUIwv1Rnr80ehs1Oqs4sWaR2nYFTs8NsbEIwTFwGYReRUIjwQyxlwc+y2W0bC8Ip/1lU0sKculN2TY36xJ1VDIcPPjb/GBYytiehQN7R6PAPr38t5/uya0fL7+PVu3rn8gp1yjIY+WfWpk0jLhquc1JNPVpL0+gGM+EvmxD8VAL8bnG5znyCqGYz48+L3zz9Fe+4or4c2/wau3wQdujyTO40UErrhHe4APXgXTjtIYtjsoLxalR2gPe2D4zEvxAvj3x+GBq3SVt0t+Dau+oOEgNz+y6J3AOzUEddeH1GDlz+oveqDTZn/kAa38KVsWGYn7yLXwzI/UiBcvjowqB/W4ihaoEYpW+QRqlFMCsP5ODat4p/4YdM1OrX6WxyNwpyUBeMVZgnX6ioixdw33zBM95zx+8LFnHKP5jc1/1/Bh1Vr93Lqa+ucT4mHOafodPnk94UkJtj2hQvDiLbpt3e/1ujOLVGxqt6go7npWK9bSslVAu5o9oaEYI4vfelhFIKtUPeIVHx782z7+P9Q727JKvcWeThWdFVfoOdwe/4HXddW1nk71Ks75dkQkjv6ghvlqN08qIbgh0Y2wKCtmaq9z0bQc6tuDVDdrMuzV3Q385pmdpKf4+dJ50YWg0Q0NdfcSChl8Ps8P1OfTP3W8+FO06seLz68D57zEIwKHSs407XmCCtBZ34jdOx+OtCwNY937MY2fzzszelWKl+xSFcHhmLECPvOsGrTsUv2jr71jcMhp4XkarnnwszqSNxoikZLe3Ao49Ys68jXYqp/52df1DxmIwOdeJbxOQjRS0mHJhZoYh/7rbg+kdIkz2duciBDkVjjrMWQ7r83WufvdgXXu51g4Tw2lSPQ8jiuo886E9/5SBwP+4SLd5iab4yU9Wz2QXc9o6CWzSL/X5ZfpNNsnfFqT6PXb4cpHNTwUbNbkecsBXQnv2I9qp+eJb0Wut7tVDfTAzsbrd6l4v+82+P2FOm5oIPPP1uO8/CudmLJhF2C0A9Ddrm3talHxa6/TDsmLt6gH4xZ2HPl+DYUe2KCiOf+ckf1/R8GwQmCMeSahLbCEOWFuIRmpfk6cV8ju+nZe3aWVDate18E8tS2xZ+1saNdkqjHQ0dNHdno8Gj8FGa0IuKQGNLT17I8jIZWxIiUt0tN+5/fh+E9GTeyz4grtleeUDX9Mnw/O+04c+8URS77k1+pZHdigIbFYpOfAF9Zp79iXoqGVsqPVuL//NmcltXr4/btV0AL54RH1iKg3akLR81nFCzWsNe8svba5Z2heaPtTkVHDI2HBuWpcl7xHw46rfwD3fkLPfeo1cMrVWuxQfqyGidb9Xr/3BeeqEBzzUR0guPxyfY87xcqd79eR1z6/GvO+bvUCzviydpKufXtwbg30mk77EvzjarhlpYY7QZPo6XlaEPHzo9Q7uOKvGk56/W5YdbXmvgJ5WjxROBeec3JtWSVw+n/rffmxQ4cyR0k8s4ieBNwCHAGkAX6g3RiTO+QbLSNmRn4Gb37nfPw+4ZltddS0dBHs7eORNzQ9UzOEELgeAUBrV8/hKwRjgT9VPYuEniMldiUSRA+DJZqUNBWAoUTAxTso8erXIjPdLrkwsv3Yj2rP46xv9jeKp39p6GMPDAu+71ea2B3OO4vGERdpqGr5ZWq4V9+koZj3/CwSJnNzG+dcDyd/Xj2JU74AJYsioS1XtCqO19BbX7d6Ud1tmm9JDWhi3E2AD9VDP+5javj/9X3Y85LmHooWqJheeT88/CU44uJI0nj55fr63VdoWEsEZp+qYnHWN2Hd7+AxZxb/9/wsIUIgwxUDichadFTwvWgF0UeBRcaYr495a+Jg5cqVZu3atRNx6nHlzpd2c93fN3HzB47mq/e/QU56CjMLM3nkmtOj7n/tva9z3zqNJz/5X2ewMEYuwWI5rKl9S0MzqYGJbolijHpHXo/NtbkDPaaeTq0USs+OjL5PSdNKsZZ9OrYkqyS6JxIHIrLOGBO1JjeuAWXGmO2A31mP4HdAlOCYZSwpy9P45I8e20phVhoXHFU2tEfQ7vEI3MohiyXZKF0yeUQA1NgPDNuJRA+bpWZEku4paZEcnM+nXk3J4lGLwHDEEz/ocOYK2iAiP0LLSBM5R5EFmJ6nP+bGjm5+94kT2FDZRH17N929IdJSBn/8DR3dZKT66ezpi1QOWSwWSxzEY9A/4uz3eaAdnVF08BwHljFlTnEW5fkZfPe9R/GORSVMy9UyzzrPmsdeGtu7mVWoick26xFYLJYREE/V0B4RyQCmG2PiKF+wjAXZ6Sm88LWzw8+n5aqHUNPSRXn+4Br6hvZuTphbxNaaVusRWCyWERHPXEMXARuAx5znK0Rk1dDvsow1pY5HEK2EtKcvREtXb9gjsDkCi8UyEuIJDd2Arj/cBGCM2QDMHeoNlrEn4hEEaXRyBaClovVtmiieWaiewlAewYPr9/F/z+9KcGstFstUIp5kcY8xpln6Z7ntBHTjTGFmGik+YWddG+/48dPkZ6Zx2sJi7ltXxbJyXfyjJCedzDQ/bcGemMe58+U91LUG+eRpVsstFosSj0ewSUSuAPwislBEbgHGDrRzzwAAIABJREFUd9UECz6fUJqTzgPr99HS1UtfyHDXq5XMKsxk7R6d2bMwM43s9JQhk8WVDR3UtnZhJ5O1WCwu8XgEXwC+iU44dxe64tj3hnyHJSGU5gbYv7eJstwAq798Jp09ffT1Gc740dO0BnspyEojO5Cii9NEobO7j7pWrTpq6eolLyN1PJtvsVgmKcN6BMaYDmPMN40xxxtjVjqPY49ssiQMt4T0wmXTSfX7yA2kUpCVxmfPmk+qXyjLDZAzhEewtzGyYEZdq/0KLRaLEtMjGK4yyE5DPf64CeP3LOs/c+Rn3zGf9x1TPqxHUFkfEYKaliALSu00FBaLZejQ0MnAXjQc9ApDznFrGQ8uOKqM7t5QeLpqFxFhujMlRXZ6CgdbO6K9ncqGyPZa6xFYLBaHoYSgDDgP+BBwBfAwcNdhvcD8JOeU+cWcMr94yH2y01NjhoYqGzpIS/HR3RuipiX6CGWLxZJ8xMwROBPMPWaM+RhwErAdWC0in4/34CJygYhsFZHtIvK1KK9/SUQ2i8hGEfmniMwe1VVYwuQEUmjtil4+urehg/kl2WSnp1BrhcBisTgMWTUkIunAhahXMAf4BfBAPAcWET9wK+pVVAFrRGSVMca7aOd6YKUxpkNEPgv8CLhs8NEs8eKWjxpjGDD2g8qGDuaVZBHs6aPGhoYsFotDTI9ARP4IvAQcC3zHqRr6njFmX5zHPgHYbozZaYzpBu4G3uvdwRjztDHGDVy/DMRYcNUSL9mBFEIGOnv6+m03xlDZ0MGswkxKctKpsx6BxWJxGKp89MPAQuAa4EURaXFurSLSEsexy9Fks0uVsy0WnwQejfaCiHxaRNaKyNq6uro4Tp28uCuT/fa5XayvbAxvr20NEuwNMasoi2m5AesRWCyWMEPlCHzGmBznluu55Yz1MpUi8mF09bMfx2jLbc4YhpUlJSVjeerDjsIsXczip09u4wePvhXevqOuDYDZhZmU5qRT2xK0o4stFguQ2AVm9qFrF7hUONv6ISLnoiOXLzbG2HjFIXLOEaX88spjeddRZWyraQ0b+41VzQAcXZ7HtNwAnT19dpZSi8UCJFYI1gALRWSus8LZ5UC/QWoicgzwG1QEahPYlqQhPcXPu4+ezknzimjq6KHWmVJiQ2UTs4syKchK80xpPVh373x5D/uaOse1zRaLZWJJmBAYY3rRVc0eB7YAfzXGbBKR74qIOyr5x0A2cK+IbLDrHIwdi8t01PBb1a0AvF7VFB6IVpqjI5QHrm1Q1xrkugff5I8v7h6/hloslgknnknnRo0x5hHgkQHbrvc8PjeR509mFk9TIdha3cKSshwONHexvEKFYEa+CkFlQweneN7jjjzetD+eWgCLxXK4YBehP0wpyEqjNCedrdVtrK9sAmDFLBWCmQWZ5ARS2Livud97qhpdIWi2iWSLJYmwQnAYs7gsh601Lbxe1USqX1g6XYu9fD5heUU+G6ua+u1f1ai5gcaOHg402/JSiyVZsEJwGLOkLIdtNW2s2rCfpTPyCKT6w68tq8jjrQOtdHkGnu31TEr35gBvwWKxHL5YITiMWTQth+7eEB3dvdxw0dJ+ry2ryKc3ZNh8IJIP2NvYweJpOfjE5gkmkubOHt7/yxfYVtM60U2xJAlWCA5jzls6jStPnMX9nz2FY2YV9HvNrSDauDcSHtrb0MnCadnMK8lm037rEUwUG6uaeK2yibW7G4ff2WIZAxJaNWSZWPIz07jxfUdHfa0sL0BpTnp4oFlfyLC/qZMLl00nxSe8vLNhPJtq8bCzrh2A+jY7vtIyPliPIIlZPjOf15z5iKpbuugNGWYWZLK4LJfqlq6Y6xpYEos7HUh9e/cEt8SSLFghSGLOWFjM7voOtte2hhPFMwszmF2UCcCe+vaJbF7SEvYIrBBYxgkrBEnMO48sA+DRN6ojQlCQyaxCFQLvGsdTgVd21nP3q5UT3YxDZqfrEdjQkGWcsEKQxEzLDXDc7AIe21TNjrp2RGB6fiDiETSMvRBsrW5l1ev7x/y4AL96ZgffevBNDk5hA9rR3ct+ZwxHfZv1CCzjgxWCJOeCI8vYtL+F257dwUlzi0hP8ZMTSKUwK409Y+QR/Gr1Dj7yf69gjOH7D2/mK/e9npCRy28daKU3ZPj7hsQIzXjghoXyM1Opb5+6gmaZWlghSHIuOKqMNL+PMxeX8tuPrQxvn1WYSWXD2OQIXtlVz3NvH+SpLbW8sP0gXT0hWrrGNhHd2N5NtTOJ3v3rqsb02ONBU0c3V925jic2VQNw/JxCGtq7CYXsVB+WxGOFIMmZWZjJC187m99+dCVZ6ZFq4tlFmYfkEfzllUpu+efbAFQ7oY6v3Pc6rl0bOPPpoeLOsnrOklI2H2hhy4HRDYjbfbA93N7xZPXWOh7bVM0v/rUdEThudgEhA02dPePeFkvyYYXAQklOOj5f/4XuZxdmsr+pk+7e0KiOee+6vfzFSdy6ayI0dvSQ6tfz1Izxmslbq9XwX3v+YkTgyc01gCaQW7riN6ZfuGs91977+pi2LR7W7mkgkOojLcVHeX4GM/IzAJswtowPVggsUZlVlEXIMOpFairrO6hp6aI92EtDezfnHlEKwEXLZwBQO8ZrJr9V3UpBZipLynJYPC2HV3c1UNXYwWW3vcyfX46/kmhPfTtrdjcQ7O0bfucxZO3uRo6fU8gdHzueb190JMXOkqMHbcLYMg5YIbBExa0c2j2KsQRtwV7q27sJGV0QB+CdS8v4y6dO5FsX6pxHY+0RbKluZUlZLiLCSfOKWLenkcc3qVdQGWf1U3uwl5auXoK9ofCI6/GgpauHrTWtHDe7gNMWFnPe0mkUZesqcg12LIFlHLBCYInK7ChjCXbUtfGDR7fQ2zd0uMj7HncthNLcdE6ZX0xhVhrZ6Slj6hGEQoZt1a3hVdlOmFtIZ08ftz27A4D9cXo11Z68xSs768esfcOxobIJY2Dl7MLwtqJs9QimYuVQe7CXjm47Kn0qYYXAEpWSnHSm5abz9w37MMbQ0d3LZ+5cx2+e2cm6PUNPhuatNnL3LcsLhLeV5qZHXS95tOxt7KCzp48lHiGAiNcRb3jLTRL7x3mupbV7GvFJZOEggILMNESmVmjIGMOtT2/nxJv+yTV3bwhvb+7o4byfPmOnNp/EWCGwREVE+O93Lua1yibufHkPX7lvIzvq2vAJPL/94JDv9VYbuXMZTcvxCEFOOjVjWDW0z1lQxx0RXZydzvySLEDXZNjf1BnXuAV3MZ7TFxazbk8jPcN4PmNBXWuQv2/Yx9IZuWR7qrb8PqEgM21KJYtf3tnAjx/fSk9fiO21beHtb9e28nZtG+sr7WyqkxUrBJaYfPDYCo4uz+P6v2/ioY0H+K9zF3HMrAKefbu/EPSFDE9urglXGFU2dJCfmUpJTjpNHT2kpfjIz0wN7z8tNxCuJBoLXAM+3am0ATh9YQmFWWm8d0U5Hd19NMdRhlndrIJyyYpyOnv6eH1v0zDvODRaunr40O0vU9sS5LoLlw56vSgrbUrlCFxxP2V+Ub/y4Drnux7L79wytlghsMTE5xP+59LlfPK0uTz+xTO4+pyFnLagmI1VTTR1RAzU/a9V8ak/ruW6B9/EGENlQwezCzMpdwzztNx0RCLlqdNyA9S0dI3Z6GI3tl+WG/E6vnrBEh65+nTmOEnveMJDB5q7KMhM5azFpfgEntlWNybti8Vz2w6yvbaNWz50DCfOKxr0elF22pSaZsL9TSyalkN7dx/tzuy1rgDUjVIIvnTPBm5/dufYNNISlYQKgYhcICJbRWS7iHwtyutniMhrItIrIh9MZFsso2NxWQ7XvWdpOBF7xqJijIEXtkeSqY+/WY3fJ9yzdi//9/wu9tR3MLMwk/ICFQKvgQYNDQV7x2508f6mTvIzU8lIiyzFmZHmpywvEK7Hd8NHQ1Hd3EVZXgZ5makcO6uA1VsTKwRuRdbJ8weLAEBRdjp1Uyg05A5+m1+aDQwWgNF6BE9tqeExZ8S1JTEkTAhExA/cCrwLWAp8SEQG+r+VwMeBvySqHZaxZXlFPjmBFP66di+hkKEt2Mtzbx/koyfP5oIjy/jBo29R1djB7KL/396Zh8VZnQv89zIDA4Rh35dAFkjMvhsTo8ao0aiJu3XptVVv21xbbb3aq7W3tt7WR622vVqrrVr3/apVU61ZTJNoVrITEhJCgIQ1LMMQYFhmzv1jvpkMhDUBBsv5PQ8P35zZXs73cd7vXU8oqcYiHN9RERiPPe6DVqfrjKyD8joHSREhnT7nUUa9yRwqq3OQbAS1F46PZ29JXb/XO/hSWNVAvNXSrqLbl3EJVgqrG9pZX2fC53vL+HSAGv4B2BpbCQ82kxTR/vx65vB0LAJHqxO7o42D5fUopcgpqfPu16DpPwbSIpgD5CulCpRSLcA7wDLfFyilCpVSe4CBj8pp+gWzKYCfXJTFuoPHeeKLPNbsr6DF6WLJ5CSeuH4KyZHBuBSkR4/o1iIAd1ZPq9PFgsfX8uczMP3L6hzexacjMSOCCDIHeDt67i+z852Xt7L49+v5aGf7nkTldoc3u+n8rDgA1h9sHw/ZXlSDs5/6/xRVN5IRM6LL588ZE4NS9FsG03PrDvOrT/cNWP8iW2MLkaFBxBuJAadaBH1Xqp731je3UVbn4D/e3MGvPs3tJ4k1HgZSEaQAR30eHzPG+oyIfE9EskUk+/jxgTXXNT3z3fkZ3DQnjefXHebe93YTG2ZhxsgowoMD+eNNM0iPCWVGelS7GIEvCR6LoN7B/jI75XYHb20p7rNV4FnQyuqaulQEIkJKZIg3RvB/24/xdX4VFfUOPthe4n2do9VJTUOL93MmJocTZ7Xwz7xK72v2ldZx7XObeGNzEQCbDlfjaD39CuTC6gZv4V5nTE2NJCTQxOZ+qmmosDuoOtHCrmMDEwS3NbUSGRroVfReRWC4t6pO9L2Jnq/y+GfecYprGjlSpS2C/uYbESxWSv1FKTVLKTUrLi7O3+IMe0SE/1k2iT/cOI2F4+L44cIxmIxeRVPTIll3/0LGxocxOs7tK07vcNebGB5MgLgL1Dx1BsU1jezqJEunrK6JX3ycc8qCu+uojYkPf8GWgmpqG1u7VAQAyZHB3hhBTkkdE5MjuHJKMjuLa713956Ml0TDxSQizMmIbldhvO2I+878/7YfY3NBNTe9sNmrFPpKQ3MblfXNZMR2bREEmQOYlRHFpsNnrgicLuW9u15t9GHqb2yNrUSEBBIZGkiQKcC7iFfamxFxy1DTRzeXb72JZ65LbY5BSe0dTgykIigB0nwepxpjmn8BzKYArpqewou3zeY780d1+ppRsSP4/J4FXHxWQrvxkCATszOiWZ1bSXZRLbFhFoLMAfxtZwm7j9qoazyZ6vnapiJe21TEmv2V7T7jyS/yaGp18tom9+KQ2EWMACA5IoRSWxMulyK31M6klHBmpkfR0OIkz+ha6k1B9VEo4xKtFNc0erNfthtV0ntL6vj533KA03fbeGotunMNAcwdHUNeRf0Z1xNUn2j2dn5dvX9gFEFdUyuRoUGICHFWC8ftzThdiuqGFkYbCq+vhYQeqyIk0ESu0VHW6VK9rhbX9I6BVATbgEwRGSUiQcC3gE8G8Ps0Q5CzksJP6WwKsHhiInkV9fzzQCXnjIlh0fh4Xt1UxLJnv+bWl7bgaHWilOLve8qA9ovXloJqvsqvwhwg3i6jyd1YBKPiRlBZ30x2US31zW1MSo5gZnoUcLLg7YCxyKRFnXTVeCqVD1a4lcWOolrmjo7GFCDkV57AajGzrbDmtHzunv2gM2K7dg3ByYyijWdoFXiqrOeOjuZgxYkB2Y+6trGFKKNeJM5qobK+mZqGFpwuxcTkCIA+Z0FV1jswB4j3fEWEuD+/vzZN0rgZMEWglGoDfgh8AewH3lNK7RORR0RkKYCIzBaRY8D1wJ9FZN9AyaMZWlw8wW0lNLQ4mZUexfILxrB0ajJ3L8okp7SOBz/cy96SOoprGokICWRtXiVtThetThePfrafOKuFO84dRYvhIkjsRhFcYnzXk1/kATApJYLUqBBiwyzsMFxTn+WUkxkfxsgYX0UQDri31yyvc1Bia+KSCYlcOD6eqNBA7r90HHVNrRyq7LvP+oixEHd0m3VkSkoEcVYLK/acWbaPx/V1zfRUAHb3c1M9l0u5LQJjoY63Wqisd3jdUROT3XPZ0z4Ue47ZuP2Vbfzq030cPn6CSnszsWEWr1JeanSv7W0jQU3v6DxvrZ9QSn0GfNZh7Bc+x9twu4w0w4y06FAmJoezr9TOzPQoJqVE8PRN0wEIMglPrjzI+oPHMQcID1w2ngc/3Mv2olrW5h1n97E6/nTLDGJGBHmzjbpKHwUYG28lKyGMrYU1BJqEzIQwRISZ6ZFsL66l0u5gW2ENd1+Y2e59qVEhhAaZOFBeT7ixwM1Mj+KaGSnUe2sg9rH1SLW3zqIj9763i/My47hqevs8iaKqRmLDLO3aSnSG2RTA0qnJvLapkNqGFqKM9tR9pcLw1882+jAV97NFUO9oQymICHXLFx9uYWthjdcCmGAogp4sghV7ylibV4lJhOoTLdQ1tRIfbmF8kvv9V05N5t3so1oR9DPfiGCx5l+TG2alkRYd4r3b83DXwrH8/PKzqG1s4dzMWK6cmkyQKYC73trB8+sOc9OckSyZnMT0kVGEBJpOKSbrjCWTkwB31avF7H7trPRoiqobeehvOSgFl09JaveegAAhK8HKgXI72YW1BAcGMCE5nMjQINKiQ0mNCiExPJgtRzqPE9Q0tPDhjhLe3nrqfghHqhsY1YNbyMPV01NodSpW7C3r1es7o8II2KZFhZAQbqGwn10rtiZ3EPikRRCMrbGVY7Xu70mPHoHVYu4xRlBw/ARZ8VYuOiuB3cdsVNY3E2+1cOXUJP50ywxmZ0QxMjp0QFxbwxmtCDR+47Z5GWz46YWYTe0vQxHhzgWj+eyeBfz2uqmEWczcsWAUWQlW7lmUycNXuusSg8wBLMiMZYyRndQdlxuKYJLhqwa4cU4aczKiWZVbQWZ8GFkJp97Vj0+0sq/EzvvbjzJ/TCyBPrKKCHNGRbNiTxlZD33Ox7va50J4ehXtPGprt9FNS5uLvcfqmGDc5fbExORwxiVY+XDHqXsx/yOnnCue2dBjGmul3UFsmAWzKYD06BHtWoX3BzYjwO/pKeVJIc0tdcde4qwWdwC5B4vg8PEGxsSPYNrISIqqGymsaiDOGozFbGLJ5CREhPToUIprdLC4PxlQ15BGcyZ4fPTg7h3UGU/eMJU2Z8/B2swEK/cvHuctFAMIDw7k9Tvn8MyafKb7tID2ZVyilXe2HSXIHMBDl591yvP3XpzF2Pgw3ss+yqsbC1k27aQLyJMO61n4Z2W43TI7i2tpanUyf2xsj3KDW+FcOzOFRz87wMGKetblHeflr4/w+Y/P4+k1h8gts7PxcBUXjk/o8jMq7A5vPcfImFA2HOq8Hqe2oYWXNxayYncpoRYTF46L5ycXZ7XrFdXp+4y0UK8iML4ru7AWq8VMSJDJm0nUFS1tLoprGrl8chJTU93no6nV6VUqHtKiQ9lcUI1SCqXoNBlB0ze0RaD5RhMeHEh0L/3mdy0cy6SUiHZjFrOJ+xaPY9FZnS+intf/cOFYb12ELxmxI7h7USa3zk1nR7GNI1UnXRa7jtq82UxbjtSwOreC3FI7X+dXESB02miuK66dkUqQKYAXNxTwx7X5lNY5uOvNHd6UylW5ld2+v8Le7G0Fnh4dSoW9maaW9laEy6X4/uvbeebLQyRHhmAxm3j6y3y2duH68sXT3TXSiBFkJVgJs5jJq6gnzljI48OD223+05HimgacLsWY+BFMTo3Ao3viOxQkpseE0tDiZN5jX/KDN7b3KJumZ7Qi0Gi6YVZ6FG/deTZ3LRzb7euumpaCCHy00+0eUkqx+5iN87LiyIwP460txfz769nc8eo21hyoZHJqpDcVsjfEhFlYPCmR97KPUdfUyvSRkXyVX4XVYub8rDjW7K/oNo21st7h7fHkyYzyBFxPNLeRV17PKxsL2VpYw+PXTuGNO8/mzTvPJio0kJe+OtKjfF7XkPE3pUaFsunBC3n82sk8vHQiAJnxYRytbaTe0XlL8PxKtxIdExdGmMVMptG8Lt7aPiMsM97twqttbGH1/gqq/NiYr9Lu4OxHV3dZpKeUGvT9r08HrQg0mm4QEeaNjfVWTndFYkQw546N5YPt7oW6qLoRW2MrU9MimTMqmhJbEymRIZTVOdhXaufcsb23BjzcPGck4N4450+3zCA0yMSNs9NYNi2Zyvpm9ho7gFXaHe22imx1uqg60eJ1DXmK2DwB1/ve283iP6znkRW5LMiM5fqZ7kS+4EATt5ydzqr9FRRXN6KU4q9fHel0nwaPIvBVbtbgQG6cPdLrjpuSGoFSkFNi7/Tv8zST81heHvdQR9fQ/LExfLD8HN7//jxcClbuG5gCud7w8sZCKuzNfLK71LtDm++mPO9sO8qc36zp1X4Y/kQrAo2mn1h+wRgq6x18+6UtPPNlPgDT0iJZMjmJUbEjeP2Os72L7PwxvYsP+DJ3dDT/eXEWD185kaSIEP553wX89NLxXDg+HlOA8MGOY1TYHVz0u3Xc/so2b+8mTy6/p8eTp79RUXUjJbYmVuaWc/nkJP7z4iyevH5qu3jAt89Jxxwg3Pf+bn6/+hCPrMjlN3/ff4pstY0tWC3mUwL/vkwxFvY9Rq8ju6OVH72909tm5PDxEySGB3tTas8ZE0OQKYDUqPapwe7U32gmpYS7q9dzTj+b6kyod7R62158lV9FdlEtv/0ij9+vOuh9zXvZR6lrah3wvS3OFB0s1mj6iXljYnn+1pksf2MH+0rtnJ8VR1aCFVOAsPa+CwD47ysnMH1kFHP7EB/wICL8aNHJWgePqyfIHMSNs9N4bVMRW4/UYHe0sbmghg2HqliQGevtVeSxCCJDgwgPNlNU08DbW4pRwINLxpMadWo6a0J4ME/dMI17393F1sIaokID2VpYQ6mtybvXA7hjBBGh3bu6okcEkRoVwp6SOpRS3PfeblbmVtDU0saLt82m4HgDo+NOFthdPT2Fc8bEEBNm6fTzRITLJiXy5/UF1DS09CpWVGJrIq/czoyRUd54xunQ6nTx1MqD1Dva+O78DF7+upBfr3B3RV2ZW071iWaa21zsNNqSrNlf4S2GG4poRaDR9COLzkpg888WYQqQTmMA4cGB3Hz2yH7/3l9cMYE9x2zklNi5f/E43tpSzC8/2Ycl0MT+MjtJEcFeVwu4K5o3HKqi3tHGovHxnSoBD0unJhMREsg/csr59tx0ljy9gb/vKeO6manUO9pIjAjG1thCVC8W1impEew5ZuPFDUdYmVtBekwo6w9WUVnv4HDlCZZNP7lYiki3hYIAS6cl89y6wzzwwR6eu3UmrU4XwYEna0pcLoWjzUlokJkX1hfwm8/c1szsjCje/ve5mE0BZBfW8MbmIn5z9WTK7Q6e/CKPR6+e7C3ee+zzA3ydX8XM9CjuWzwOl1Jc/9wm8irquW5mKnctHMvLXxey+1gd09Ii2XXU5o0Veb5r7YFKWp2udunHQwmtCDSafqa3WUz9SXCgiZdum83K3ApunjOSpIhg7n1vN5NTInj06slcMyOl3QJ5flYcL35VQHCgie+dN6bHzz8/K87r65+cEsELGwr43aqDNLU6EYEAEc7phZUzJTWSz/aW89sv8rh4QgL/ccEYrv7TRm5/ZRv1zW3ewr/eMj4xnIevmMAvP81l/mNfUm53cM7oGG6bl0FqVAg/+2gvpbYmXrv9bP53zSEWZMYyf2wsj31+gKdWHeSm2SP5/uvbqW5oYWJyBLuO2vg8p5yUyBB+fsUESm1NvLChgJTIEF7ZWEic1YI12J0N9cebp3PFFLfi8lTJ3794HE+tzOOFDQUEmgKYkBTOHeeO4gdv7OCOV7PJLbXz5PVTuGBcPC1tLoLMQ0MxSH/tGztYzJo1S2VnZ/tbDI1myNNbd0lf+etXR3hkRS6LxsdzycQESmwOSm1NLJ6Y6O0h1RUb86u4+cUtRIUGsvIn5xMbFsS5j6+lxNbEZZMSee7Wmacl0x9WH2RHsY1xCWF8vKvU27U0PNiM06VwKkVzm4svfnweWQlW7n13Fx8ad+1hFjPpMaGU2pqwNbVitZhxtLpYe/8FvLaxkBc2FLD+pwv5rw/2UFjVSHhIIAECf797Qbs5+Xh3KR8tn8fWwhp+8XEOZTYHDywZz1XTUpj+P6twuhSJRgptUkQwlfZmXrhtVrvaFnDHHh75NJd1B48jAvddMo7rZ6VxpojIdqXUrE6f04pAo9H0BadLkVdez1lJ1h4LzTrS0NzGDX/exN2LMlk8MRFwNwN86asjrPzJeaRF967tRnc0tznZWWzjQJmdiyYksL2olnve2cXSqcneflatThercyvYWljDxWcloIBbXtxCcGAA739/Htc+v5H06FDK6xycNy6OZ2+ewed7y1j+5g4AHlk2kX87J6PXMm06XE3UiEBSo0L59YpcahpaOFLVQFmdg8UTE9lWWMN352cwMjqUR1bkcrSmkSunJnO0ppEdxTaumZGCxRzAsmkppxVfAq0INBrNEKalzUVdU6u38Gwg+Dq/iimpEViDOw9oK6VY/sYOJiaH86NFmazYU8pf1hewv8zOB8vnMSU1klani/mPfYmtqZVtP7uox+B4T5Tamrjq2a+pa2olK8HqTf/NiAnlieumMmdUNC1tLh78cC//yCkj1GLmgUvHc+3M0+vTqRWBRqPRnAZOl2pXQ7Iqt4KahmZunN0/Af/ahhYCRAgPMfPRzhKa21xcNzN1QILK3SkCHSzWaDSaLuhYSNhTDKSv+LYVv2aG/zryD42QtUaj0Wj8hlbeW56DAAAGoElEQVQEGo1GM8zRikCj0WiGOVoRaDQazTBHKwKNRqMZ5mhFoNFoNMMcrQg0Go1mmKMVgUaj0QxzvnGVxSJyHCg6zbfHAlX9KE5/MlRl03L1DS1X3xmqsv2ryZWulIrr7IlvnCI4E0Qku6sSa38zVGXTcvUNLVffGaqyDSe5tGtIo9FohjlaEWg0Gs0wZ7gpgr/4W4BuGKqyabn6hpar7wxV2YaNXMMqRqDRaDSaUxluFoFGo9FoOqAVgUaj0Qxzho0iEJFLRSRPRPJF5AE/ypEmImtFJFdE9onIPcb4L0WkRER2GT9L/CBboYjsNb4/2xiLFpFVInLI+B01yDKN85mTXSJiF5Ef+2u+ROSvIlIpIjk+Y53Okbh52rjm9ojIjEGW67cicsD47o9EJNIYzxCRJp+5e36Q5ery3InIg8Z85YnI4oGSqxvZ3vWRq1BEdhnjgzJn3awPA3uNKaX+5X8AE3AYGA0EAbuBCX6SJQmYYRxbgYPABOCXwH1+nqdCILbD2BPAA8bxA8Djfj6P5UC6v+YLOA+YAeT0NEfAEuBzQIC5wJZBlusSwGwcP+4jV4bv6/wwX52eO+P/YDdgAUYZ/7OmwZStw/NPAb8YzDnrZn0Y0GtsuFgEc4B8pVSBUqoFeAdY5g9BlFJlSqkdxnE9sB9I8YcsvWQZ8Kpx/CpwlR9lWQQcVkqdbmX5GaOUWg/UdBjuao6WAa8pN5uBSBFJGiy5lFIrlVJtxsPNwKDvhdjFfHXFMuAdpVSzUuoIkI/7f3fQZRMRAW4A3h6o7+9Cpq7WhwG9xoaLIkgBjvo8PsYQWHxFJAOYDmwxhn5omHd/HWwXjIECVorIdhH5njGWoJQqM47Lgf7dtLVvfIv2/5j+ni8PXc3RULrubsd95+hhlIjsFJF1IrLAD/J0du6G0nwtACqUUod8xgZ1zjqsDwN6jQ0XRTDkEJEw4APgx0opO/AcMAaYBpThNksHm3OVUjOAy4C7ROQ83yeV2xb1S76xiAQBS4H3jaGhMF+n4M856goReQhoA940hsqAkUqp6cC9wFsiEj6IIg3Jc9eBm2h/0zGoc9bJ+uBlIK6x4aIISoA0n8epxphfEJFA3Cf5TaXUhwBKqQqllFMp5QJeYABN4q5QSpUYvyuBjwwZKjympvG7crDlMrgM2KGUqjBk9Pt8+dDVHPn9uhOR7wBXALcYCwiG66XaON6O2xefNVgydXPu/D5fACJiBq4B3vWMDeacdbY+MMDX2HBRBNuATBEZZdxZfgv4xB+CGL7Hl4D9Sqnf+Yz7+vWuBnI6vneA5RohIlbPMe5AYw7uebrNeNltwMeDKZcP7e7Q/D1fHehqjj4B/s3I7JgL1PmY9wOOiFwK/BRYqpRq9BmPExGTcTwayAQKBlGurs7dJ8C3RMQiIqMMubYOllw+XAQcUEod8wwM1px1tT4w0NfYQEfBh8oP7uj6Qdya/CE/ynEubrNuD7DL+FkCvA7sNcY/AZIGWa7RuDM2dgP7PHMExABrgEPAaiDaD3M2AqgGInzG/DJfuJVRGdCK2x97R1dzhDuT41njmtsLzBpkufJx+48919nzxmuvNc7xLmAHcOUgy9XluQMeMuYrD7hssM+lMf4K8IMOrx2UOetmfRjQa0y3mNBoNJphznBxDWk0Go2mC7Qi0Gg0mmGOVgQajUYzzNGKQKPRaIY5WhFoNBrNMEcrAo2mAyLilPYdT/utW63RxdKfNQ8azSmY/S2ARjMEaVJKTfO3EBrNYKEtAo2mlxj96Z8Q954NW0VkrDGeISJfGk3U1ojISGM8Qdz7AOw2fuYZH2USkReMfvMrRSTEb3+URoNWBBpNZ4R0cA3d6PNcnVJqMvBH4A/G2DPAq0qpKbgbuz1tjD8NrFNKTcXd936fMZ4JPKuUmgjYcFetajR+Q1cWazQdEJETSqmwTsYLgQuVUgVGY7BypVSMiFThbpPQaoyXKaViReQ4kKqUavb5jAxglVIq03j8X0CgUurXA/+XaTSdoy0CjaZvqC6O+0Kzz7ETHavT+BmtCDSavnGjz+9NxvFG3B1tAW4BNhjHa4DlACJiEpGIwRJSo+kL+k5EozmVEDE2LTf4h1LKk0IaJSJ7cN/V32SM/Qh4WUTuB44D3zXG7wH+IiJ34L7zX46726VGM6TQMQKNppcYMYJZSqkqf8ui0fQn2jWk0Wg0wxxtEWg0Gs0wR1sEGo1GM8zRikCj0WiGOVoRaDQazTBHKwKNRqMZ5mhFoNFoNMOc/wdwldekhFU2PAAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "tags": [], | |
| "needs_background": "light" | |
| } | |
| }, | |
| { | |
| "output_type": "stream", | |
| "text": [ | |
| "\n", | |
| "Minimum test loss: 0.25792052820324896\n", | |
| "\n", | |
| "Best train correlation: 0.9643175626613927\n", | |
| "\n", | |
| "Best test correlation 0.921231229609948\n" | |
| ], | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "Mi3XxCeW0tZb" | |
| }, | |
| "source": [ | |
| "**Activity**\n", | |
| "\n", | |
| "Do hyperparameter tuning on this model. Compete with other groups. The group with the best test loss wins." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "Q5Sp9lZCrQDZ" | |
| }, | |
| "source": [ | |
| "Exercise: \n", | |
| "\n", | |
| "1. Modify the example of a random forest regressor [here](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html) and use it to train a discriminative sequence to function model fo PETase catalytic activity and stability.\n", | |
| "2. How many trainable parameters are there in the feed-forward neural network? In the random forest?\n", | |
| "3. Based on the performance of the two models and your answer to 2, which is a better model? What does better mean?\n", | |
| "4. Interpret the trained random forest with `feature_importances_` and `decision_path`\n", | |
| "5. What are some issues with how we split our data into train/test?" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "e4EQDb13siOJ" | |
| }, | |
| "source": [ | |
| "### Rough Mount Fuji Model" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "jdOIU_Srq-7i" | |
| }, | |
| "source": [ | |
| "How can we overcome the problem of having limited labeled data? We can do a systematic evaluation of the power of our discriminative models with a benchmark sequence-to-function landscape, otherwise called a fitness landscape. One such landscape that was shown to match natural fitness landscapes is the [Rough Mount Fuji (RMF) model](https://iopscience.iop.org/article/10.1088/1742-5468/2013/01/P01005). It's a model with tunable epistasis. Much like how optimization algorithms are benchmarked on a variety of test functions, we will do the same on the RMF. This process removes the uncertainty around what the ground truth for unseen sequences are: we know them since we defined the RMF function. It also enables us to explore aspects of model generalization that would be difficult to do otherwise." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "nS9b-z1hzQJT" | |
| }, | |
| "source": [ | |
| "An RMF works by fixing a reference sequence $\\boldsymbol{\\sigma}^*$ which is assigned some arbitrary fitness value, say $\\xi$. Then for every sequence $\\boldsymbol{\\sigma} \\neq \\boldsymbol{\\sigma}^*$, we compute its fitness as $f(\\boldsymbol{\\sigma}) = \\xi - c\\cdot d(\\boldsymbol{\\sigma}, \\boldsymbol{\\sigma}^*) + \\mathcal{N}(0, a)$." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "IKwryJ9A09JD" | |
| }, | |
| "source": [ | |
| "where $d(\\cdot, \\cdot)$ is the Hamming distance and $c, a$ are the parameters of the RMF. Intuitively, the fitness of a sequence drops linearly by distance. The epistatic component is contributed by the Gaussian. Thus, $c$ tunes how smooth/linear the landscape is and $a$ tunes how rough/epistatic the landscape is. A perfectly smooth landscape has $c/a \\gg 1$. A completely epistatic landscape, otherwise known as a House of Cards (HoC) model, has $c/a \\ll 1$." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "pFIIcM6kv6ka" | |
| }, | |
| "source": [ | |
| "import numpy as np\n", | |
| "import random\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "\n", | |
| "random.seed(105716977)\n", | |
| "\n", | |
| "optim = np.zeros((263, 20)) \n", | |
| "for i in range(optim.shape[0]):\n", | |
| " rand_aa = random.randint(0, 19)\n", | |
| " optim[i][rand_aa] = 1\n", | |
| "\n", | |
| "def rmf(optim, curr, c=1.0, a=0.01, best_val=100.):\n", | |
| " \"\"\"calculates the fitness of curr according to the Rough Mount Fuji model\"\"\"\n", | |
| " dist = float(np.sum(np.abs(curr - optim))/2) # should be the Hamming distance\n", | |
| " noise = np.random.normal(loc=0., scale=a)\n", | |
| " return best_val - c*dist + noise\n" | |
| ], | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "4nZMpL1Vxr91" | |
| }, | |
| "source": [ | |
| "To quantify the ruggedness of a fitness landscape, one of the measures we can use is the roughness to slope ratio ($r/s$). It describes how well a linear model can explain the landscape, so that more rugged landscapes have a higher $r/s$. A perfectly linear landscape has $r/s = 0$." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "ENE-ALuDsj3h" | |
| }, | |
| "source": [ | |
| "# roughness to slope ratio (r/s)\n", | |
| "from sklearn.linear_model import RidgeCV\n", | |
| "def get_rs_ratio(x, y):\n", | |
| " \"\"\" Expects x to be of shape (batch_size, seq_len, 20)\n", | |
| " \"\"\"\n", | |
| " x = x.reshape(x.shape[0], -1)\n", | |
| " clf = RidgeCV(alphas=[1e-3, 1e-2, 1e-1, 1], cv=5).fit(x, y)\n", | |
| " s = clf.coef_ # 237*20 coefficients\n", | |
| " # print(s.shape)\n", | |
| " l = len(s)\n", | |
| " s = np.mean(np.abs(s))\n", | |
| " r = clf.predict(x)\n", | |
| " r = np.sqrt(np.mean((r - y)**2))\n", | |
| " print(\"r/s: \", r/s)" | |
| ], | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "e1iP6J7G8gKy" | |
| }, | |
| "source": [ | |
| "This defines a fitness landscape for all possible sequences! We can use this to construct a simulated dataset, and test *any* algorithm. However, using it as a benchmark in CbAS was difficult in practice. (Why?)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "l2Q4CAA68J4k" | |
| }, | |
| "source": [ | |
| "# plot c/a versus r/s \n", | |
| "# smooth is c/a >> 1, rough is c/a << 1\n", | |
| "a = [0., 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.01]\n", | |
| "cs = [1., 10., 5., 3.6, 2.3, 1.0, 0.5, 0.]\n", | |
| "rs = []\n", | |
| "for i in range(len(a)):\n", | |
| " # a = c/ratio\n", | |
| " X_rmf = make_data(optim)\n", | |
| " y_rmf = []\n", | |
| " for x in X_rmf:\n", | |
| " score = rmf(optim, x, c=cs[i], a=a[i])\n", | |
| " y_rmf.append(score)\n", | |
| " y_rmf = np.array(y_rmf)\n", | |
| " X_rmf = np.argmax(X_rmf, axis=2) # convert to integers to save space\n", | |
| " np.save(\"X_rmf_{}.npy\".format(1), X_rmf)\n", | |
| " np.save(\"y_rmf_{}.npy\".format(1), y_rmf)\n", | |
| " rs.append(get_rs_ratio(X_rmf, y_rmf))\n", | |
| "\n", | |
| "plt.plot(range(1, 9), rs)\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "_O817tP48ql_" | |
| }, | |
| "source": [ | |
| "Exercise:\n", | |
| "\n", | |
| "1. What is the relationship between the RMF model and Kauffman's NK landscape?" | |
| ] | |
| } | |
| ] | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment