DeepLift + Mutagenesis Analysis Notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "266a8e62",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import random\n",
    "import os \n",
    "import torch\n",
    "import torch.nn as nn\n",
    "from torch import Tensor\n",
    "import torch.optim as optim\n",
    "from torch.utils.data import Dataset, DataLoader, SubsetRandomSampler, ConcatDataset\n",
    "import torch.nn.functional as F\n",
    "from torch.cuda.amp.grad_scaler import GradScaler\n",
    "from torch.cuda.amp import autocast\n",
    "\n",
    "from torchsummary import summary\n",
    "from sklearn.metrics import r2_score\n",
    "\n",
    "from ray import tune\n",
    "\n",
    "import json\n",
    "import itertools\n",
    "from itertools import groupby\n",
    "import gzip \n",
    "from io import BytesIO\n",
    "from time import time \n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# import pyBigWig\n",
    "from scipy.sparse import csc_matrix\n",
    "import math \n",
    "\n",
    "import torch.cuda.amp as amp\n",
    "\n",
    "import functools\n",
    "\n",
    "import captum \n",
    "from captum.attr import DeepLift, DeepLiftShap\n",
    "# import seaborn as sns\n",
    "from sklearn.metrics import confusion_matrix, precision_score, f1_score\n",
    "\n",
    "from copy import deepcopy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "32be34d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from modules import * "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "115972ab",
   "metadata": {},
   "source": [
    "To run a classification, we need to update the last block in the model (conv_block_5). Pytorch's cross entropy loss accepts logits as input, so we need the model's very last layer to be linear. \n",
    "I tried to update the conv_block_5 attribute using class inheritannce and attribute overriding, but it did not work for me. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5ba5d5c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "class upd_GELU(nn.Module):\n",
    "    def forward(self, input: Tensor) -> Tensor:\n",
    "        return torch.sigmoid(torch.Tensor([1.702]).cuda() * input) * input\n",
    "\n",
    "def ones_(tensor: Tensor) -> Tensor:\n",
    "    return torch.ones_like(tensor)\n",
    "\n",
    "def zeros_(tensor: Tensor) -> Tensor:\n",
    "    return torch.zeros_like(tensor)\n",
    "\n",
    "class BasenjiModel_Classif(nn.Module):\n",
    "    def __init__(self, num_targets, n_channel=4, max_len=128, \n",
    "                 conv1kc=64, conv1ks=15, conv1st=1, conv1pd=7, pool1ks=8, pool1st=1 , pdrop1=0.4, #conv_block_1 parameters\n",
    "                 conv2kc=64, conv2ks=5, conv2st=1, conv2pd=3, pool2ks=4 , pool2st=1, pdrop2=0.4, #conv_block_2 parameters\n",
    "                 conv3kc=round(64*1.125), conv3ks=5, conv3st=1, conv3pd=3, pool3ks=4 , pool3st=1, pdrop3=0.4, #conv_block_2 parameters\n",
    "                 convdc = 6, convdkc=32 , convdks=3, debug=False):                 \n",
    "        super(BasenjiModel_Classif, self).__init__()\n",
    "        \n",
    "        self.convdc = convdc\n",
    "        self.debug = debug\n",
    "        self.num_targets =  num_targets\n",
    "        ## CNN + dilated CNN\n",
    "        \n",
    "        self.conv_block_1 = nn.Sequential(\n",
    "            upd_GELU(),\n",
    "            nn.Conv1d(n_channel, conv1kc, kernel_size=conv1ks, stride=conv1st, padding=conv1pd, bias=False),\n",
    "            nn.BatchNorm1d(conv1kc, momentum=0.9, affine=True),\n",
    "            nn.MaxPool1d(kernel_size=pool1ks),\n",
    "            nn.Dropout(p=0.2))\n",
    "                \n",
    "        self.conv_block_2 = nn.Sequential(\n",
    "            upd_GELU(),\n",
    "            nn.Conv1d(conv1kc, conv2kc, kernel_size=conv2ks, stride=conv2st, padding=conv2pd, bias=False),\n",
    "            nn.BatchNorm1d(conv2kc, momentum=0.9, affine=True),\n",
    "            nn.MaxPool1d(kernel_size=pool2ks),\n",
    "            nn.Dropout(p=0.2))\n",
    "        \n",
    "        self.conv_block_3 = nn.Sequential(\n",
    "            upd_GELU(),\n",
    "            nn.Conv1d(conv2kc, round(conv2kc*1.125), kernel_size=conv3ks, stride=conv3st, padding=conv3pd, bias=False),\n",
    "            nn.BatchNorm1d(conv3kc, momentum=0.9, affine=True),\n",
    "            nn.MaxPool1d(kernel_size=pool3ks),\n",
    "            nn.Dropout(p=0.2))\n",
    "        \n",
    "\n",
    "        self.dilations = nn.ModuleList()\n",
    "        for i in range(convdc):\n",
    "            padding = 2**(i)\n",
    "            self.dilations.append(nn.Sequential(\n",
    "                upd_GELU(),\n",
    "                nn.Conv1d(conv3kc, 32, kernel_size=3, padding=padding, dilation=2**i, bias=False),\n",
    "                nn.BatchNorm1d(32, momentum=0.9, affine=True), \n",
    "                upd_GELU(),\n",
    "                nn.Conv1d(32, 72, kernel_size=1, padding=0, bias=False),\n",
    "                nn.BatchNorm1d(72, momentum=0.9, affine=True), \n",
    "                nn.Dropout(p=0.25)))\n",
    "            \n",
    "        self.conv_block_4 = nn.Sequential(\n",
    "            upd_GELU(),\n",
    "            nn.Conv1d(72, 64, kernel_size=1, stride=1, padding=0, bias=False),\n",
    "            nn.BatchNorm1d(64, momentum=0.9, affine=True), \n",
    "            nn.Dropout(p=0.1)) \n",
    "            \n",
    "        self.conv_block_5 = nn.Sequential(\n",
    "            upd_GELU(),\n",
    "            nn.Linear(64, self.num_targets, bias=True)) \n",
    "\n",
    "    \n",
    "        self.conv_block_1[1].weight.data = self.truncated_normal(self.conv_block_1[1].weight, 0.0, np.sqrt(2/60)) #4\n",
    "        self.conv_block_2[1].weight.data = self.truncated_normal(self.conv_block_2[1].weight, 0.0, np.sqrt(2/322)) # conv1kc\n",
    "        self.conv_block_3[1].weight.data = self.truncated_normal(self.conv_block_3[1].weight, 0.0, np.sqrt(2/322)) # conv1kc\n",
    "        self.conv_block_4[1].weight.data = self.truncated_normal(self.conv_block_4[1].weight, 0.0, np.sqrt(2/72)) # 72\n",
    "        self.conv_block_5[1].weight.data = self.truncated_normal(self.conv_block_5[1].weight, 0.0, np.sqrt(2/64)) # 64        \n",
    "        self.conv_block_1[2].weight.data = ones_(self.conv_block_1[2].weight)\n",
    "        self.conv_block_2[2].weight.data = ones_(self.conv_block_2[2].weight)\n",
    "        self.conv_block_3[2].weight.data = ones_(self.conv_block_3[2].weight)\n",
    "        self.conv_block_4[2].weight.data = ones_(self.conv_block_4[2].weight)\n",
    "\n",
    "        \n",
    "        for i in range(convdc):\n",
    "            self.dilations[i][1].weight.data = self.truncated_normal(self.dilations[i][1].weight, 0.0, np.sqrt(2/218)) # 72\n",
    "            self.dilations[i][-2].weight.data = self.truncated_normal(self.dilations[i][-2].weight, 0.0, np.sqrt(2/32)) # 32\n",
    "            self.dilations[i][2].weight.data = zeros_(self.dilations[i][2].weight)\n",
    "            self.dilations[i][-2].weight.data = ones_(self.dilations[i][-2].weight)\n",
    "\n",
    "    \n",
    "    def truncated_normal(self, t, mean, std):\n",
    "        torch.nn.init.normal_(t, mean, std)\n",
    "        while True:\n",
    "            cond = torch.logical_or(t < (mean - 2.28*std), t > (mean + 2.28*std))\n",
    "            if not torch.sum(cond):\n",
    "                break\n",
    "            t = torch.where(cond, torch.nn.init.normal_(torch.ones(t.shape), mean=mean, std=std), t)\n",
    "        return t\n",
    "\n",
    "\n",
    "    def forward(self, seq):\n",
    "        if self.debug: \n",
    "            print (seq.shape)\n",
    "        seq = self.conv_block_1(seq)\n",
    "        if self.debug: \n",
    "            print ('conv1', seq.shape)\n",
    "        seq = self.conv_block_2(seq)\n",
    "        if self.debug: \n",
    "            print ('conv2', seq.shape)\n",
    "        seq = self.conv_block_3(seq)\n",
    "        if self.debug: \n",
    "            print ('conv3', seq.shape)\n",
    "        for i in range(self.convdc):\n",
    "            if i == 0:\n",
    "                y = self.dilations[i](seq)\n",
    "            if i >= 1:\n",
    "                y = y.add(self.dilations[i](seq))\n",
    "            if self.debug: \n",
    "                print ('dil', i, self.dilations[i](seq).shape)\n",
    "        if self.debug:\n",
    "            print ('y', y.shape)\n",
    "        res = self.conv_block_4(y)\n",
    "        if self.debug: \n",
    "            print ('res', res.shape)\n",
    "        res_lin = res.transpose(1, 2)\n",
    "        if self.debug: \n",
    "            print ('res_lin', res_lin.shape)\n",
    "        res = self.conv_block_5(res_lin)\n",
    "        if self.debug: \n",
    "            print ('res', res.shape)\n",
    "        return res\n",
    "        \n",
    "    def compile(self, device='cpu'):\n",
    "        self.to(device)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7a24cb1",
   "metadata": {},
   "source": [
    "First step for the analysis is to initialize and load the model, trainer and get the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d0e055b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "param_vals = { \n",
    "\"mode\" : \"classification\",\n",
    "\"optimizer\" : \"Adam\", \n",
    "\"init_lr\": 0.001, \n",
    "\"optimizer_momentum\": 0.9, \n",
    "\"weight_decay\": 1e-3, \n",
    "\"loss\": \"bce\", \n",
    "\"num_targets\": [1],\n",
    "\"seq_len\": 128*128*8,\n",
    "\"target_window\": 1024*32,\n",
    "\"batch_size\": 1,\n",
    "\"cut\": 0.8,\n",
    "\"num_workers\": 0,\n",
    "\"num_epochs\": 1,\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3d29aad4",
   "metadata": {},
   "outputs": [],
   "source": [
    "memmap_data_contigs_dir = '/data/users/goodarzilab/darya/work/Datasets/hg38_memmaps'\n",
    "targets_memmap_data_dir_cl = '/data/users/goodarzilab/darya/work/Datasets/hg38_targets_memmaps_CL.ATAC'\n",
    "targets_memmap_data_dir_pdx = '/data/users/goodarzilab/darya/work/Datasets/hg38_targets_memmaps_PDX.ATAC'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "100bc9dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_files_dir = [targets_memmap_data_dir_cl, targets_memmap_data_dir_pdx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "265c8bd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = BasenjiModel_Classif(num_targets=1)\n",
    "model.compile('cuda')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c9e585d3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "init dsets\n"
     ]
    }
   ],
   "source": [
    "trainer = Trainer(param_vals, model, memmap_data_contigs_dir, target_files_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "47bd9f10",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# trainer = Trainer(param_vals, model, memmap_data_contigs_dir, targets_memmap_data_dir_cl, targets_memmap_data_dir_pdx, mode='classification')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "672e9870",
   "metadata": {},
   "outputs": [],
   "source": [
    "trainer.load_model(trainer.model, \"basenji_toy_upd_motif_32k_recfield_13epochs.pt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "945b02e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_loader, val_loader = trainer.make_loaders(augm=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "fc27f157",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "batch = next(iter(train_loader))\n",
    "X, y = trainer.get_input(batch)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c19af53",
   "metadata": {},
   "source": [
    "Then, for the baseline data, we shuffle out inputs. Since the trainer outputs one-hot incoded data, we reverse that step, randomly permute the values and one-hot encode them again. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "147da04d",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_ints = torch.argmax(X, dim=1).squeeze()#.cpu()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "a5d24f81",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = torch.randperm(X_ints.shape[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e09706b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# X_shuffled = torch.stack([X_ints[i][idx] for i in range(X_ints.shape[0])])\n",
    "# X_shuffled = torch.transpose(F.one_hot(X_shuffled), 2, 1)#.unsqueeze(0)\n",
    "\n",
    "X_shuffled = X_ints[idx]\n",
    "X_shuffled = torch.transpose(F.one_hot(X_shuffled), 1, 0).unsqueeze(0)\n",
    "# X_shuffled = F.one_hot(X_shuffled).unsqueeze(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22e5c8b3",
   "metadata": {},
   "source": [
    "Next, we use the DeepLift algorithm from the Captum package."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f55315c",
   "metadata": {},
   "source": [
    "### DeepLift Algorithm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "cc494bad",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/users/goodarzilab/darya/anaconda3/envs/basenji_pytorch/lib/python3.9/site-packages/captum/_utils/gradient.py:56: UserWarning: Input Tensor 0 did not already require gradients, required_grads has been set automatically.\n",
      "  warnings.warn(\n",
      "/data/users/goodarzilab/darya/anaconda3/envs/basenji_pytorch/lib/python3.9/site-packages/captum/attr/_core/deep_lift.py:320: UserWarning: Setting forward, backward hooks and attributes on non-linear\n",
      "               activations. The hooks and attributes will be removed\n",
      "            after the attribution is finished\n",
      "  warnings.warn(\n",
      "/data/users/goodarzilab/darya/anaconda3/envs/basenji_pytorch/lib/python3.9/site-packages/captum/attr/_core/deep_lift.py:452: UserWarning: An invalid module MaxPool1d(kernel_size=8, stride=8, padding=0, dilation=1, ceil_mode=False) is detected. Saved gradients will\n",
      "                be used as the gradients of the module's input tensor.\n",
      "                See MaxPool1d as an example.\n",
      "  warnings.warn(\n",
      "/data/users/goodarzilab/darya/anaconda3/envs/basenji_pytorch/lib/python3.9/site-packages/torch/nn/modules/module.py:1025: UserWarning: Using a non-full backward hook when the forward contains multiple autograd Nodes is deprecated and will be removed in future versions. This hook will be missing some grad_input. Please use register_full_backward_hook to get the documented behavior.\n",
      "  warnings.warn(\"Using a non-full backward hook when the forward contains multiple autograd Nodes \"\n",
      "/data/users/goodarzilab/darya/anaconda3/envs/basenji_pytorch/lib/python3.9/site-packages/captum/attr/_core/deep_lift.py:452: UserWarning: An invalid module MaxPool1d(kernel_size=4, stride=4, padding=0, dilation=1, ceil_mode=False) is detected. Saved gradients will\n",
      "                be used as the gradients of the module's input tensor.\n",
      "                See MaxPool1d as an example.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "dl = DeepLift(trainer.model)\n",
    "attribution = dl.attribute(X, baselines=X_shuffled.float(), target=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "526f4675",
   "metadata": {},
   "outputs": [],
   "source": [
    "nucs = np.arange(4)\n",
    "inp_vals = X.cpu().detach().numpy()\n",
    "attr = attribution.cpu().detach().numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11187339",
   "metadata": {},
   "source": [
    "Then, for each minibatch, we collect the attribute values. They are grouped for each nucleotide, so it's the same thing as reversing one-hot encoding. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "15ceb23b",
   "metadata": {},
   "outputs": [],
   "source": [
    "nucs_vals = np.zeros((inp_vals.shape[0], inp_vals.shape[-1]))\n",
    "attr_vals = np.zeros((inp_vals.shape[0], inp_vals.shape[-1]))\n",
    "\n",
    "\n",
    "for j in range(inp_vals.shape[0]):\n",
    "    for i in range(len(nucs)):\n",
    "        ids = np.where(inp_vals[j][i] == 1.)\n",
    "        nucs_vals[j][ids] = nucs[i]\n",
    "        for idx in ids: \n",
    "            attr_vals[j][idx] = attr[j][i][idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e763a81b",
   "metadata": {},
   "source": [
    "Finally, we need to cross-reference the deeplift values with the actual locations of the motifs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "26b03f47",
   "metadata": {},
   "outputs": [],
   "source": [
    "import re\n",
    "motif_str = '202131'\n",
    "\n",
    "def return_ids(test_str, motif_str):\n",
    "    ind_found = np.hstack([[m.start() + i] for m in re.finditer(motif_str, test_str) for i in range(len(motif_str))])\n",
    "    return ind_found\n",
    "\n",
    "def get_motif_loc(seq):\n",
    "    seq_subset_str = \"\".join([str(int(s)) for s in seq])\n",
    "    test_empty = np.zeros(len(seq_subset_str))\n",
    "    ids = return_ids(seq_subset_str, motif_str)\n",
    "    test_empty[ids] = 1\n",
    "    return ids, test_empty\n",
    "\n",
    "def get_true_targets(labels_again):\n",
    "    ids, test_empty = get_motif_loc(labels_again)\n",
    "    rand_ids = ids[0::6]\n",
    "    diff_ids = np.diff(rand_ids)\n",
    "    bins = np.arange(0, labels_again.shape[0], 128)\n",
    "    # bin_indices = np.digitize(rand_ids[np.where(diff_ids >= 180)[0]], bins) #rand_ids[diff_ids_range], bins)\n",
    "    bin_indices = np.digitize(rand_ids[np.where(diff_ids >= 160)], bins) #rand_ids[diff_ids_range], bins)\n",
    "    targets = torch.zeros(len(bins))\n",
    "    targets[bin_indices] = 1\n",
    "    return targets, rand_ids[np.where(diff_ids >= 160)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "fd647449",
   "metadata": {},
   "outputs": [],
   "source": [
    "targets = get_true_targets(X_ints)[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "4f70cfa0",
   "metadata": {},
   "outputs": [],
   "source": [
    "targets_empty = torch.zeros(X_ints.shape)\n",
    "targets_empty[targets] = 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "7a98cdca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD4CAYAAAD2FnFTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAx7klEQVR4nO2de5Bc1Xngfx8CiTcSZsAyD0s4snflVAoTFcab2OWNk1hiXVaSzUOqJBDsXUUVqErWW7WB9R/Opta1eaxDwoIlYyMeMSCwsY2CxUOIpwE9RuiBXiONpNHMSKPRSPN+T898+8e9Ld3pud19b/e9fW93f7+qqek+73P63Pudx3e+I6qKYRiGYeRyQdIFMAzDMNKJCQjDMAzDFxMQhmEYhi8mIAzDMAxfTEAYhmEYvlyYdAGi4JprrtEFCxYkXQzDMIyqYseOHWdUtSGff00IiAULFtDY2Jh0MQzDMKoKETleyN+WmAzDMAxfTEAYhmEYvpiAMAzDMHwxAWEYhmH4YgLCMAzD8CWQgBCRpSLSJCLNInKfj7+IyIOu/x4RudV1v1FE3hCRAyKyT0T+0hPnahHZJCKH3f/zPH73u2k1iciXo6ioYRiGEY6iAkJEZgEPA8uAxcBKEVmcE2wZsMj9WwWscd0zwH9X1X8P3A7c44l7H7BZVRcBm93vuP4rgE8DS4HvumUwDMMwKkiQGcRtQLOqHlXVcWA9sDwnzHLgSXXYAswVkfmq2qGqHwCo6gBwALjeE+cJ9/MTwO943Ner6piqHgOa3TLEx/g47NkTTVqqsHu3k2YhpqZg504nfBwcOgT9/efLtGsXZDLQ2QmtrfHkGaRMfX3J5A3T2yGNDAzA5s1w6lTSJfEnyuckCMPDsH9/fv8DB6ClBY4XVOUvnYEBaGoqHGbPHhgbiyf/FBBEQFwPtHm+t3P+JR84jIgsAD4DbHWdrlPVDgD3/7Uh8kNEVolIo4g0dnV1BahGATZuhJ/8BNraioctRmsr/PSn8PLLhcNt2QIvvAAffFB+nn48/TQ88ojzubkZfvYzeO01WLMG1q2LJ88wZUqCQ4ecdnj99eTKUIhHH4V33oG1a5MuiT8//7nznLS3Vya/9evhuedgcHCm3+goPPssPP44PPZYPPk/9hg880z+QdyJE057vPhiPPmngCACQnzcclusYBgRuRx4HvgrVe2PID9U9RFVXaKqSxoa8p4UD8bAgPO/2Kg/CNnRRDbNfAwNOf9HRsrPMx/ZBytbJr8HrdJk650Eo6PO/zS0gx+9vUmXoDDZGWkUz0kQsu0xOTnTb2oq/vy7uwv7Z9uh2LNexQQREO3AjZ7vNwAng4YRkYtwhMNTqvoTT5hOEZnvhpkPnA6Rn2EYhhEzQQTEdmCRiCwUkdk4G8gbcsJsAO50tZluB/pUtUNEBHgUOKCq/+QT5y73813ACx73FSIyR0QW4mx8bwtdM8MwDKMsihrrU9WMiNwLvALMAtap6j4RWe36rwU2AnfgbCgPA3e70X8N+FPgQxHZ5br9T1XdCPwd8JyIfB1oBf7ATW+fiDwH7MfRgrpHVX3mmIZhGEacBLLm6r7QN+a4rfV8VuAen3i/wH9PAVU9C3wpj9+3gW8HKZthGIYRD3aS2jAMw/DFBISXuM4kpIVar59hGJFiAiJqxHdFrfRwRnRYm5dHpduvUH6VKEuxPOqgP5mAMAzDMHwxAWEYhmH4YgLCMAzD8MUEhGEYhuGLCQjDMAzDFxMQcRFUpdRUTyuPtXl5VLr9kv69iuWfdPlixAQE1L66Wq3XzzCMWDABYRiGYfhiAsIwDKNEDp7qp7V7OOlixIYJCMMwjBLpGZrgQEexO9CqFxMQhmEYhi8mIAzDMAxfTEB4qWF1NaD262cYRqQEEhAislREmkSkWUTu8/EXEXnQ9d8jIrd6/NaJyGkR2ZsT51kR2eX+tWRvnBORBSIy4vFbSzVh1lzTi7V5edRb+5k11+I3yonILOBh4LeAdmC7iGxQ1f2eYMtw7o5eBHwWWOP+B3gceAh40puuqv6RJ4/vAH0e7yOqekvIuhiGUUskbe7bCDSDuA1oVtWjqjoOrAeW54RZDjypDluAuSIyH0BV3wa68yUuIgL8IfBMKRUwDMMw4iGIgLgeaPN8b3fdwobJx+eBTlU97HFbKCI7ReQtEfm8XyQRWSUijSLS2NXVFTCr2mI8M0XfyETSxTAMo0YJIiD85nK5u51BwuRjJdNnDx3ATar6GeAbwNMicuWMxFUfUdUlqrqkoaEhYFa1xY92tLHuF8eSLoZhGDVKEAHRDtzo+X4DcLKEMDMQkQuB3wOezbqp6piqnnU/7wCOAJ8MUM6643T/WNJFMAyjhgkiILYDi0RkoYjMBlYAG3LCbADudLWZbgf6VLUjQNq/CRxU1fasg4g0uBvjiMjNOBvfRwOkVZ2Y6mnlsTY3wlCkv0gN96eiWkyqmhGRe4FXgFnAOlXdJyKrXf+1wEbgDqAZGAbuzsYXkWeALwLXiEg78C1VfdT1XsHMzekvAH8rIhlgElitqnk3uSMhDo2INHUa0/gwoqTezH3XMUUFBICqbsQRAl63tZ7PCtyTJ+7KAun+mY/b88DzQcplGIZhxIedpDYMwzB8MQFhGIZh+GICwjAMw/DFBIRhGKEYHs/wQWsPapvHNU+gTeq6odY7fK3Xz6gIbzZ10XT5Bdx09aVcc/mcpItjxIjNIKLGrLmmF2vz8nDbbzwzBcBUrQ84zJqrCQjDSA118MIJhVlzTRwTEIZhGIYvJiAMwzAMX0xAGIZhGL6YgDAMwzB8MQGRNLWuCZJGrM2NMBTrLzXcn0xAQO1rRNR6/QzDiAUTEHFRw6MKo94xc9/1ggkIwzBKQnxvGjZqCRMQhmGUhFZ6JmFUnEACQkSWikiTiDSLyH0+/iIiD7r+e0TkVo/fOhE5LSJ7c+L8jYicEJFd7t8dHr/73bSaROTL5VTQMAzDKI2iAsK9H/phYBmwGFgpIotzgi3DuTt6EbAKWOPxexxYmif5B1T1Fvdvo5vfYpyrSD/txvtu9o5qwzAMo3IEmUHcBjSr6lFVHQfWA8tzwiwHnlSHLcBcEZkPoKpvA2HulF4OrFfVMVU9hnPP9W0h4pdOrW+G1Xr9DMOIlCAC4nqgzfO93XULG8aPe90lqXUiMi9MWiKySkQaRaSxq6srQFYVwqy5phdr8/Kot/Yza66BBIRfK+QORYOEyWUN8AngFqAD+E6YtFT1EVVdoqpLGhoaimRlGIaRHo6fHWLvib6ki1GUIBcGtQM3er7fAJwsIcw0VLUz+1lEvg+8WGpahlETiNgyoJcaNvf9kw9OAPDL11+VcEkKE2QGsR1YJCILRWQ2zgbyhpwwG4A7XW2m24E+Ve0olGh2j8Lld4GsltMGYIWIzBGRhTgb39sClNMwDMOIkKIzCFXNiMi9wCvALGCdqu4TkdWu/1pgI3AHzobyMHB3Nr6IPAN8EbhGRNqBb6nqo8A/iMgtOMtHLcCfu+ntE5HngP1ABrhHVScjqa1hGIYRmEB3UrsqqBtz3NZ6PitwT564K/O4/2mB/L4NfDtI2QzDSAY7SV372EnqpLE158pjbW6EoUh/qWUxaQICqn7Dqyi1Xj8jEczURu1jAiIubJRq1CqV7tv2LCWGCQjDMAzDFxMQhmEYhi8mIAzDMAxfTEAYhmEYvpiA8FLrm2G1Xj/DMCLFBETUmDXX9GJtXh457VfzB+XMmqsJCMMwDMMfExCGYRiGLyYgDCMt1MGSRShq2Nx3tWACwjCMkjBTG7WPCQjDMAzDFxMQSWOqp5XH2twIQ7H+UsP9KZCAEJGlItIkIs0icp+Pv4jIg67/HhG51eO3TkROi8jenDj/KCIH3fA/FZG5rvsCERkRkV3u31riptbXM2u9foZhxEJRASEis4CHgWXAYmCliCzOCbYM52rQRcAqYI3H73FgqU/Sm4BfVtVfAQ4B93v8jqjqLe7f6oB1qWmaTg3w0OuHyUxOJV0UwzDqhCAziNuAZlU9qqrjwHpgeU6Y5cCT6rAFmJu9c1pV3wa6cxNV1VdVNeN+3QLcUGolUknE0853DncxMakMT9jtq0ayiJn7rhuCCIjrgTbP93bXLWyYQnwNeMnzfaGI7BSRt0Tk834RRGSViDSKSGNXV1eIrIozMTnF1JR1SqM22Xeyj8aWGWO20NT8SWojkIDw6wW5b88gYfwTF/kmkAGecp06gJtU9TPAN4CnReTKGYmrPqKqS1R1SUNDQ5CsAvPQ6828vO9UpGkaRlp4dV8n7xw+k3QxjCogiIBoB270fL8BOFlCmBmIyF3AV4A/VnXmkao6pqpn3c87gCPAJwOUM1KaTg1UOkvDMIxUEURAbAcWichCEZkNrAA25ITZANzpajPdDvSpakehREVkKfDXwFdVddjj3uBujCMiN+NsfB8NXKNyqPW1zlqvn2EYkXJhsQCqmhGRe4FXgFnAOlXdJyKrXf+1wEbgDqAZGAbuzsYXkWeALwLXiEg78C1VfRR4CJgDbBJHDXOLq7H0BeBvRSQDTAKrVbX8BdNKYdZc04u1eXnUW/uZNdfiAgJAVTfiCAGv21rPZwXuyRN3ZR73X8rj/jzwfJBy1SO13yWD0TcyQUffCP/uozO2pwzDiIhAAsIw0saPGtsYGM3wqeuuQOpgJJdGzBZT7WOmNoyqZGA0UzxQtWGCbjrWHoljAsIwjESZMOsAqcUERJVgCkhGLdLRN8JDrzdztGsw6aIYPpiASBp781cea/NIiOIk9cneUQBau4eLhEyQIv2l4qZHKogJCKj9tc5ar59hGLFgAsIwDCNFDI9n6BoYS7oYgAmI+Ih42mmTACM1pNyaq6pWtbHNJ947zg+3HE+6GIAJiKqhhpc5jQqgqkyW8dJ88v0W3m2uDgN/P9t1gn/ZfDjpYpTMaIpM+puAMOqaU32j9AyNJ12M2PnwRB8Pbj7MwOhESfHPDo7TN1Ja3ELEMTNuOZPiDe8qwwSEUdc8s62Vx99rSboYsXPQtU7cO1zeS34w4gOKNjNONyYgvNR6b631+hmx8+GJvnOfzdRG7WMCwjAMw/DFBETUmLnv9GJtXhaaczCu5q8cNXPfJiBy2X+yP3DYXW29PPbusaLhVJWOvpFyinUOs1xq1BLWndONCYgcXglxF/UbB08H2vTbe6Kf9dvaOBLA3szZwTFaz5oWhlHfDI1l6C9R4yqNPLu9lYderz7V20ACQkSWikiTiDSLyH0+/iIiD7r+e0TkVo/fOhE5LSJ7c+JcLSKbROSw+3+ex+9+N60mEflyORVMA2eHnFORQYTJk+8f5/kP2uMukpFGYhxOV9tA/cn3j9N4vDfpYkTGyd5RJiarb1O/qIBw74d+GFgGLAZWisjinGDLcO6OXgSsAtZ4/B4HlvokfR+wWVUXAZvd77hprwA+7cb7bvaOaqO+mJoq73BX1IxlJtl69GxVn9JNG7mKdeo6pOmwWD0TZAZxG9CsqkdVdRxYDyzPCbMceFIdtgBzRWQ+gKq+DfjdKb0ceML9/ATwOx739ao6pqrHcO65vi1EnaqLgKqnQVQKa+2hemrrcR4sciK2JM3dEtV93zl0hveOnKW5wFLhgY5+nqiDcxXFmJic4oFNh/iwva94YJcjXYP882uHOTOYDjtE5yjWX2pYfTyIgLgeaPN8b3fdwobJ5TpV7QBw/18bJi0RWSUijSLS2NXVVbQStc7hzgHWvHkkss3wNHBmMF0nnDNTzsU2hWY1L+89RXcdnMwuxvC4M1jZ1uI3NjyPd1Wt+bQjeE/1jcZWLiMcQQSE3/Jl7hMSJExQAqWlqo+o6hJVXdLQ0FBiVtkcq22FdiZtPc7G9ul+n9FXDdTPMIzKE0RAtAM3er7fAJwsIUwundllKPf/6TLSqnlqXufcMOqA946cObfPUg0EERDbgUUislBEZuNsIG/ICbMBuNPVZrod6MsuHxVgA3CX+/ku4AWP+woRmSMiC3E2vrcFKGe68HSCKM4u1Nr+QjUT5PmuppdAMbYcPZtsAaq4LQ909LPv5Pl9mK1HuxmdqJ47uC8sFkBVMyJyL/AKMAtYp6r7RGS1678W2AjcgbOhPAzcnY0vIs8AXwSuEZF24Fuq+ijwd8BzIvJ1oBX4Aze9fSLyHLAfyAD3qGpVvx2jeFlk3HXvs4NjXD4n/89WvY+SkVbeP1KCgAjZEWu13768d+a5qmpa8S0qIABUdSOOEPC6rfV8VuCePHFX5nE/C3wpj9+3gW8HKVu9kd0g3Zr0qC4hhsYyPPL20YRyr6In2zAiwE5Se6miqex7JY3qqqd++TANodrFxG/6MAFRZVT/K94wZiJY304jJiAiwLsJFWSBsfn0IANj0V68YlpOAYho8bee7kGYtn9WTYvnUVCiNdeNH3bwQWtPDAWqPCYgIuDVfZ0F/XO1mP5t90le8tm8qjWq/fJ4o3LUUi9pOjXAW021cXjXBERS1NITkYeff9hR1ZfHG+VTy3OO/jKvb60GTECEoHtoPHHjccNjk4xn8utRp0n//nBncfPm1USYFZYU/QxGTKzf3lqRfAo973FjAiIgg2MZnnivhTebThcPnEOQl/boxCSvHyy8VAXw2oFOntlWmY5pVJhqWeOPSfrl1t5kLLR1D/PwG820dSdzR4wJiIBkTzKf6I3YGJ77sG1v6WZ3WzDLl6bqWSZlvuDqaXZQal1raiO/SCNIjHXN2lhr70nGCKcJiJRge7mGUR9Uy0QRTEA4xPyLVeQe6UJZVFOPNNJPiP5kXa8wYTX9Tg+McqhzIMYSTccERAE27K57I7JGSNL4QuwbnmA0wY3OStHWPcxgxOeL4mZ7Sw//svkwY5lg5uae2tLKz/cUs4MaHYFsMdUrR06f18IJ/dyHXLxNk/ZRNRFVq01NKe80n+FXPz4vrzHEIH0gjT/junePRZqeVLqSAfP78Y52rri4ul5pH55w9h1Hx6eYc2H6bla2GUSN0O5qOfTWgW52ubzVdJqf7Twxza2tZ5gPjvfw2v7immRBKOUVevzsEH01+vupKq8f7KQnZgWLgdHqmkH0jzi/9/ik/wwvaQsJJiBqhOz1nO1Ra1nVIEPjkxw7MzTNLTtInSpzdFzOElNbzwj7O/rLyj9qSm2N3GbsGhxjd1sfP/+wcssjSZPvHg2/l/5EHgGRNCYgvMQ0dU7N8lFaylEGcVchjXsIRnz0DI1z8FR+oXxmcIwDJ0sT2iXdo5EyTECEJPXv2NQXMJ1Yq9UnT7zfwksf5reLdrhzkE0Holl2LIekzpUEEhAislREmkSkWUTu8/EXEXnQ9d8jIrcWiysiz4rILvevRUR2ue4LRGTE47c2N79U4zME9VNz1Rwne0FVgDKnBxVRV04z2eqXMQhJWz8vWJUiv7fWQX8ouuUvIrOAh4HfAtqB7SKyQVX3e4Itw7k7ehHwWWAN8NlCcVX1jzx5fAfwHiM+oqq3lFWzmKiDPlHXBNkUtEmaUWmS2qwOMoO4DWhW1aOqOg6sB5bnhFkOPKkOW4C5IjI/SFxxhmV/CDxTZl0Mo2RSs08UM3VSzVRTTYPMIALieqDN873ddQsSJkjczwOdquq1C71QRHaKyFsi8nm/QonIKhFpFJHGrq7asL1eKvXycjMqT7l9S879l2yCvv5GOgkiIPx+w9xeky9MkLgrmT576ABuUtXPAN8AnhaRK2ckovqIqi5R1SUNDQ15Cx81sb2L6+gd/8CmQ2yK6LxBWintxWqvy2kUuLFtSw1oCIUhzZvU7cCNnu83ALk2KPKFKRhXRC4Efg94NuumqmOqetb9vAM4AnwyQDljY9/JPp58/3g8idfp6H/viWCWa+PA7yRw1qXQ9L8WXt+RLG9EuEaSbfepKSUT0CZR06kBth6rnIAoenK8hp/hIAJiO7BIRBaKyGxgBbAhJ8wG4E5Xm+l2oE9VOwLE/U3goKq2Zx1EpMHd3EZEbsbZ+D5aYv0iYcfx8PfLjmem6BsJcSq2jGeuhvtnKqn0aG5icoodx3uq6vrWsCXd095H06nKGaEzglFUQKhqBrgXeAU4ADynqvtEZLWIrHaDbcR5iTcD3wf+olBcT/IrmLk5/QVgj4jsBn4MrFbV7hLrF4wYdo3eOdzFul+EsIET57NfoV2xY2eGSrpQySjM+0fO8vahLpoqaMUzKrKqwWnYmJ2YnGLvib7E9+xS0BSBCWTZSlU34ggBr9taz2cF7gka1+P3Zz5uzwPPBylX1EQ5QovqatLM5FRR/fu0jCuz9o2++KlrEy5JeNI8CxtzLbFG0aeSrmeQ7Idissj6bvMZdrb2csnsWXyi4fJY8qg17CS1hz151sWTfKgefuMI33v7SHIFKJG4HvKgyzuTU8p7R85U/D5fP331l/d28MCmQxUtR6zE/EC8F2IDund4nJ2tPWQC2DIaGXdMaid5x3O1YQLCQ3+YPYMiRDWNnFJlbKJwhy51yjyemTp3lWrUvLw3v/mCSrD/ZD9bj3aH2Mx02jD7EvGj1GWSAx3lLQ0luSTim3MJ5SnreSiQX8vZYUYnQu73pZBmz9UCacIEREXJ39GTUGP7wS+OsubNeGYnYwmP0jJTU+7/YO16sncUgI6+0VD5HO4c4In3WhLfQH52eyvfeys9M82k1/nD0j00XrHf0G/JuBRFmEpgAsJDxS9CiYjApc6pX7GZSTkkvSkZ9pcMY+bbG/TV/Z10D40zMTW9LSvdk072jjJcYPYTJ6rKwGiwEby37R7YdIitRwvrn+xs7Sl5lvvzPR209wyf+16oTz7xXgvvNJ8pKZ9axgREAKppNFQ9JZ3JsTNDBU0v+1Hstwkqp3a29obKt1yefL+FXxyujRfS3hP9/OCdY3T2h5t9AUVf/qMTU7T35L/jpNDve6hzgB81ttMfUHidtLtUZmACImoKDp3l/AstwiH2SEz7COVQSu1+tvNEQdPLYfCVGyHbfGB0IjJttFzODo6zvaU07e3JqXAX3UdFPuul2VF6d8y3xZXKpn3TT+0HHu8V6y8RT5NPD4zO2GxPeiZuAsJDxu+EbcSzh0KjoVLxXlMZZX/q6BthaCwzbZpeL0xOKT945xiv7jsvsAqfsq7ck/zg5sOs3+6YOJua0siuSfUjqu4ft6n0NMycT/SOBLoZLl9LPLWllddScPeEl+q64TtGJqeUtu5huGq6+7vN0R7pz9eB0raKNToxyfpt5+0s/tcv3Mzlc9LXXX60o52Vt93k4+M0aKkvpuzM4UhXadol7xzu4sZ5l7LouitKil+M7HLOmcGxcxffB6F3JIJRfhkv++yAK8ll23zFzxbp6a2tXHzRBfxeiDRf3tvBgY4BFn9shtm4UIRVkogbm0G4ZKb8X9yHck6vltex47rSNPoccgVZED3zJDiVsgcqu9m9u62PF/dEc/9ylL1m84H4Trr7vXjzvYzDqP6WO/nIjV/sEe7sH+X42XCz5mx9zgyOhYoXmITkqQmIAKRscD8DPxXZ0YnJRF/qSa6dTkxOlX32oBCV7g9x7YMkSedAugR72hlP6Fk2AeGSbw059hddBFNtvyR2tvayu7237LSzVPpGq9yZWxjeauqia8AZyfmWusQ2z7bB0RKXnYpn4N/GWSN2LWeH/ApVeXLaL0hrllTMCo8y8p5Fitiaa6Fq5Uuq0lp2WUxAhERE2NnaE2ppo9V9sAdGM7GMPvOlOZ6JLrfRTDhNqXIFysYPgy/P7Dg+XRtoaDweMx9Zwi4/hKWte5iX986sf5AN0Kgp9QBnvlhZ94kIDlIOj2c8exrh4r6675Sz5xgxSV0NGhcmICD0SOXNpi6e2dYaOPyZQWdjsBQ98UjIqd/hzgEefqM5VBL5VBh7h8cZDGF3KeiBtDAP/NuH8p8niHsQGsdma3vPSKxLZLn0Do/T0VfZMwD7ToY77+LH7rY+RkId9nQ6g6LsO9nPj3e0Fwmfn6e3+j//SV3sExcmIEIS5QthPDN1blRYTrJhy3RmcHyawTKvKue5NHO+53vRPvZuC99/u7TrOs4MjtFyxmfZpATMAFvpPPZuyzSNtShJ13g67HOSf8M5scFehTEBUUFyp587jvfw1iHnPu0wqoq5lCuyohjNlcK/vn+cn+48wQu7TgSOk08W5tvILXfKPzFZWyPCOCnU1kkf+CqFAx0zn4vO/lFe8RlQRUXaelsgASEiS0WkSUSaReQ+H38RkQdd/z0icmuxuCLyNyJyQkR2uX93ePzud8M3iciXy61kmvEeciuVaWquZfawPe29PLDpEJPlvhjzvBD8nI92TZ9F/PNr+U1j/3RncGFSTQRRj6y19e1Kcf6ZKHKvSoAuv2HXSfYXGFAF+Y3SduapEEUFhHv958PAMmAxsFJEFucEW4ZzNegiYBWwJmDcB1T1FvdvoxtnMc5Nc58GlgLfzV5BmgS9wxOl/aBl9gLV4OYU/m137hXhpbPtmLPhm2u+o5Kd2i8vVY1EbbescxOFTlKX+e5O2jx6GHKNWhbqG2cGx6dtBsf1LOU7xwTQGuFmdLH9tiB7EFsivE9738m+WJe7gswgbgOaVfWoqo4D64HlOWGWA0+qwxZgrojMDxg3l+XAelUdU9VjONeY3haiTiUTxJprT4ARfxTaJntP9PMvmw+Xnc40InzLqypvu8tjfkQ91m083sP/ez3cxjrMfHEfK2PPo9AlSHEI0KZTA/QkbN8oino1Hu+OffYTxjJx0iP40/3lH6Z7cY8zKHx1X2feDfMoCCIgrge8O1jtrluQMMXi3usuSa0TkXkh8kNEVolIo4g0dnXlf1Elgd/aZZagnTOsVdOoKbZJ3TcyUVEb9gcLtCnEv8Z9emCUw52VvdRl44cdPPn+ca6fewkANzdcVtH885LT2JXYX9h8oLOg5l1aVm2iFoT5UqtUXwwiIPzKOOP9kSdMobhrgE8AtwAdwHdC5IeqPqKqS1R1SUNDg0+UZJhSGBrLc2ZAzqt5arH10BLz7x4a53QJU84DHf0FT+y+9OEpdrf1nvte9OxQ6BLEw7RWLuNN9tSWwqO0uF6SU6pceYljA+uiWTMf18ruS/jnFWZEXmq/2NPex3hmirNxmbLwofgzGn/bF2qv0xU4jR5EQLQDN3q+3wDkLnrnC5M3rqp2quqkqk4B3+f8MlKQ/CrGvEsvCtSps6qmH0Q1qi7jDfuUZ8r5wKZDNJ8urlP/8t5T5/YfwF919vWDwe34nIjBam0hol42qPQMLq6rX2MhwTWafKbt07J9PxbyQGkuuc9doaZ+bns8qsleggiI7cAiEVkoIrNxNpA35ITZANzpajPdDvSpakehuO4eRZbfBfZ60lohInNEZCHOxve2EutXNj3DE4Huqh7LTDE1pXQUGb0n8WwFVWMdHs+cs36a9DptWHJvdCuXoEbtvO1Uzn7BwGi8p7+N4kTR5XvL1ErMVasu1ZpwVBS136yqGRG5F3gFmAWsU9V9IrLa9V8LbATuwNlQHgbuLhTXTfofROQWnN+lBfhzN84+EXkO2A9kgHtUNfXDqzVvHmHBNZdyYU4v8564DjPKifZEZrCcgy6TnOgd4XCOraSXApjGaDkz5Jgz/vVg+YQhb9FjHlpm77ze3tJd9PpMozTSMjuoBN6Z5I7jPedsikEyJtIDGfh3VVA35rit9XxW4J6gcV33Py2Q37eBbwcpW5poOTPML+W45apVbtofTJ0xqRF8EFVSv6ntwVPnBUa+jlwpW/dhTH+USs/QOI+/13Lue1LG1KJmPDNF99A4H73q4tBxi1ymCJT2kov6UaiWC7BytQS3Hqv8AMROUleY3Cmk5On+SV0ckr34PqjNJD+e8Lw4oyJMabxqhH6buPnaPAzH/CyrlolGMFTuC7AcWoiX9nbwzLbWQHsi3UPx2nDy2wQu1aTKsMeAY1g7V4Lyb7tPTtujm+5fGfbmWFuI+5Y+MAEROYWWhirxg/rnGz5OOVcQBDkrkksYcxuGw+jEJK8fnH5F5bpfHCsrzeyMN8gAYeOHHazf1ubb5+OaAb+wqzR9le+9dbT4en6BQjefHuTd5vxGIStBEu+P9N0hmQQiVWkrJjBu3U73j3KiNZiWVaXXO3PNbUSF93c9WWHNqrh5Zd8pTvZGO9PMziDTpKQQ1aNpBh3DYzOIlNJwxZzI0zzQ0c9owBOnKXo/RMbLMRpZy0cUgjbfkkje8zbA1JQWPPnth9cuWDmlHhrLTDszA9VprC9t5DZhJQSezSAiJqqRV5QPVNCk4jp01Xx6kF+69vKS4/cOj3N2sLAKqYiwaX8ncy5M15hnf5ET4HHx+sHTfHiij3v+Y67aRH6GJ8rb3M/2nqe3tk5TFGg5Ux2bworTV43zpOtpqnHCvH6jsNeSSxjZFeUSw6sBNbfy8di7LUXDXHiBsPdE3wzzH0kPXIupvr5+sDOwUcZcCg0isuvthYzY5TLdKnDp2kaV0CLLl3e5RGn4Mmr6R8u3/BwWExAxEcTwX6XoGZ5gYnIqlEGzWllkGpmYZNP+zkgswRYi33S/mOHG3W19tEe8NzI4ljm3lxCGMZ86+Hfj8H2jHEH9WpkDjLRR8jWuCTyStsTkUsoD5UfP8DjXRpJSdExkppw14RDG3lIk38oiez5hfgl6/WljYnKKrUeLm4p+q+m8/nwxG1Jefua5a0Nx7sa+fM7MV8TwuL/QCyIEVGfq9xvppa4FxOjEJD/ccpzfHx7nwMl+WFh+muMFLtpJcqPOOX8R/K2fL+Tj75anRpkUSQm83GzPDI7xkctml5TWe0fOMuIZyATpTqUu90xOKj/e0c7cSy+a4Zdrbyi7HNpfwFxIdhVtcCxTUSvAaaRa9mSgzpeY2nuGGRjNcDjKjanAO8JJr47PZJdH8yTfZl0pZxxSQ8JtfqJ3hH99/zg7czR8gjIRUGslimpmz0F4bQtpnoSzh/MK3S9+PIaDhVEzYxBRrCFT+AxHTV0LiEpz/Gz1jByaToU7bZp2orVtVRpZNdLOEk/JZ8o5vRiSUnJq7R7Ou7kd1PihkS7qeonJqB+ivHayVJpdzaKDOcL3jaZoX55BhfvUlNKZ506BUpfkKijDjApQ1wIiLRuxuZZRjeg53DnIVQmXYTjPfkB3QteKvnfkLNtb/NVw0zDjCkMUBxKTsJaadup6iSnbHZJeSXxxT3FT2UZhSjEuGPT8QSk39PmRtiXret8szqWq99dioq4FhFEZJjLxj8zeP1Jc/TMoJ3unn0s4kNL9mHIFzsc/cml+zyobTEd9lsRwqGsBkR10VvZe3zyFqGHKMR0elD3tfcUDZSlSnOd3tJdXmFKp8BSjUHap6JUh2iPWU8ZFL2BPRWvFQiABISJLRaRJRJpF5D4ffxGRB13/PSJya7G4IvKPInLQDf9TEZnrui8QkRER2eX+rc3NLyrOrbNekLK5f9SkYG3jrRQejsr3WM/QdvQLU7vvBKD0+iW1jl/rv0dSFBUQIjILeBhYBiwGVorI4pxgy3Dujl4ErALWBIi7CfhlVf0V4BBwvye9I6p6i/u3utTKGenhgxpb7w5qSdP74sqa5r7i4mh0QxKd+aYMM7IXD0FmELcBzap6VFXHgfXA8pwwy4En1WELMFdE5heKq6qvqmpWrWMLcEME9QlF9uHtG05Gi8RIlnyj3Vxnv83cciZlN11dYO0/JZR6U5wN5GuLIALiesB7CXG76xYkTJC4AF8DXvJ8XygiO0XkLRH5vF+hRGSViDSKSGNXV2nLF9kXwWABu/pG/RHlnonf1Z1RqbV6r9GMmjeb0rckaFSeIALCdwk2YJiicUXkm0AGeMp16gBuUtXPAN8AnhaRK2ckovqIqi5R1SUNDQ1FqmAY6SGq+8bLNatd6MKhUrG9gNoiyGJoO3Cj5/sNQK7R9HxhZheKKyJ3AV8BvqTufF9Vx4Ax9/MOETkCfBJoDFDWUMR6GMielNRTzi90ZjD6+zrC4txRXHotipkiz5tvgb4dyzNlj1JiBJlBbAcWichCEZkNrAA25ITZANzpajPdDvSpakehuCKyFPhr4Kuqes4Ogog0uJvbiMjNOBvfR8uqZR5y+7nUek80oTWNQsblinG4M/lN0XKvnIxDuS2Oi66M5Cg6g1DVjIjcC7wCzALWqeo+EVnt+q8FNgJ3AM3AMHB3obhu0g8Bc4BNzkiILa7G0heAvxWRDDAJrFbVwtdypYpgT10+y5hGjOS0+VBEd4BUK8WucZ1BgD6bBptXUVHsGdU6eIQD6dup6kYcIeB1W+v5rMA9QeO67r6X5arq88DzQcplGEa6CHVg0Ug9dX2S+tor5yRdBCNB3jhoJqgNoxB1LSDsoJFhRIutnNYWdS0gDMMwjPyYgDAMIzJG6nzjv9YwAZEwhXTKjXiwNjfCIEW6Sy33JxMQUPMLp6ZSWx2o7YlNw/pt8tS1gKi2axUNwzAqSV0LCMMwDCM/JiAMwzAMX0xAGIZhGL6YgIiJmjf8Z9QvFe7atawllHbqWkDMuLilxjuiPWiGYYShrgVEueaS/QiummcqfJXGxGN51J8abpH61oEabl0LCMMwDCM/JiAMwzAMX+paQJg1V8MwjPwEEhAislREmkSkWUTu8/EXEXnQ9d8jIrcWiysiV4vIJhE57P6f5/G73w3fJCJfLreShmEYRniKCgj3fuiHgWXAYmCliCzOCbYM5+7oRcAqYE2AuPcBm1V1EbDZ/Y7rvwL4NLAU+G72jmrDMAyjcogWUX0Ukc8Bf6OqX3a/3w+gqv/HE+Z7wJuq+oz7vQn4IrAgX9xsGFXtEJH5bvxP5aYvIq+4abyfr4xLlizRxsbG0JU/uPswJ3/wQy7pPX/l9cjcq0On42XW+Bizh4eKphVlnoXSH5l7NReNDHHh2PTL5OPIM0yZkmD28BCzxscSLUMh4u4T5VLp8mXzm7j4EjIXXzLdU5VL+npiLU82/9Er56IXzBxLX5CZYM7gQGz5B6XnYx/nP9//tZLiisgOVV2Szz/IndTXA22e7+3AZwOEub5I3OtUtQPAFRLXetLa4pPWNERkFc5shZtuuilANWay6PqrOXLVPMYvvYyrTrbRfePNcEH5+xJXHz/CQMNHmbj0srxhRq6Yy9VtR+n92E1MXXRR2XnmcklvN2NXXMnIVfMYuWqeU6Zr5zNnsJ/Zw0OMXDWveCIxlikJsu3Qf93HZr5wUsDUrAu57Oxp9IILEmujQoxfchlXdbTRfdMnKqKlnZk9hytOdzBw3cd8/acJiBjaK/t7DM/7SN4wcwYHGLzmOsYvuzzy/IPyiU98NLa0gwgIv66QO+3IFyZI3FLyQ1UfAR4BZwZRJE1fZl3zEf7T//6rUqIahmHUPEE2qduBGz3fbwBOBgxTKG6nu7SE+z97g3yQ/AzDMIyYCSIgtgOLRGShiMzG2UDekBNmA3Cnq810O9DnLh8VirsBuMv9fBfwgsd9hYjMEZGFOBvf20qsn2EYhlEiRZeYVDUjIvcCrwCzgHWquk9EVrv+a4GNwB1AMzAM3F0orpv03wHPicjXgVbgD9w4+0TkOWA/kAHuUVW76NYwDKPCFNViqgZK1WIyDMOoZ4ppMdX1SWrDMAwjPyYgDMMwDF9MQBiGYRi+mIAwDMMwfKmJTWoR6QKOl5HENcCZiIqTJFaPdFEr9YDaqYvVYzofV9WGfJ41ISDKRUQaC+3kVwtWj3RRK/WA2qmL1SMctsRkGIZh+GICwjAMw/DFBITDI0kXICKsHumiVuoBtVMXq0cIbA/CMAzD8MVmEIZhGIYvJiAMwzAMX+paQIjIUhFpEpFmEbkv6fIAiMiNIvKGiBwQkX0i8peu+9UisklEDrv/53ni3O/WoUlEvuxx/1UR+dD1e1BExHWfIyLPuu5bRWRBjPWZJSI7ReTFaq2HiMwVkR+LyEH3d/lcldbjv7l9aq+IPCMiF1dLPURknYicFpG9HreKlF1E7nLzOCwi2SsKoqzHP7p9a4+I/FRE5qamHqpal3845sePADcDs4HdwOIUlGs+cKv7+QrgELAY+AfgPtf9PuDv3c+L3bLPARa6dZrl+m0DPodzS99LwDLX/S+Ate7nFcCzMdbnG8DTwIvu96qrB/AE8F/cz7OBudVWD5xre48Bl7jfnwP+rFrqAXwBuBXY63GLvezA1cBR9/889/O8iOvx28CF7ue/T1M9EnkJpuHPbdxXPN/vB+5Pulw+5XwB+C2gCZjvus0HmvzKjXP3xufcMAc97iuB73nDuJ8vxDmRKTGU/QZgM/AbnBcQVVUP4EqcF6vkuFdbPbL3w1/t5vGi+2KqmnoAC5j+Yo297N4wrt/3gJVR1iPH73eBp9JSj3peYso+MFnaXbfU4E4PPwNsBa5T55Y+3P/XusHy1eN693Ou+7Q4qpoB+oD8N7OXzj8D/wOY8rhVWz1uBrqAx9ylsh+IyGXVVg9VPQH8X5zLuTpwbn18tdrqkUMlyl7p98TXcGYE08qUk3fF6lHPAkJ83FKj8ysilwPPA3+lqv2Fgvq4aQH3QnEiQ0S+ApxW1R1Bo/i4JV4PnFHYrcAaVf0MMISznJGPVNbDXZ9fjrNU8THgMhH5k0JR8pQp6d8jCFGWvWJ1EpFv4tyi+VQZZYq0HvUsINqBGz3fbwBOJlSWaYjIRTjC4SlV/Ynr3Cki813/+cBp1z1fPdrdz7nu0+KIyIXAVUB3xNX4NeCrItICrAd+Q0R+WIX1aAfaVXWr+/3HOAKj2urxm8AxVe1S1QngJ8B/qMJ6eKlE2SvynnA3jb8C/LG6a0AF8q5YPepZQGwHFonIQhGZjbOhsyHhMuFqIzwKHFDVf/J4bQCymgd34exNZN1XuNoLC4FFwDZ3yj0gIre7ad6ZEyeb1u8Dr3s6ZSSo6v2qeoOqLsBp29dV9U+qsB6ngDYR+ZTr9CWc+9Krqh44S0u3i8ilbv5fAg5UYT28VKLsrwC/LSLz3FnYb7tukSEiS4G/Br6qqsM59Uu2HlFtIFXjH3AHjpbQEeCbSZfHLdOv40z99gC73L87cNYRNwOH3f9Xe+J8061DE642g+u+BNjr+j3E+ZPzFwM/AppxtCFujrlOX+T8JnXV1QO4BWh0f5Of4WiBVGM9/hdw0C3Dv+Jox1RFPYBncPZOJnBGw1+vVNlx9gWa3b+7Y6hHM87+wC73b21a6mGmNgzDMAxf6nmJyTAMwyiACQjDMAzDFxMQhmEYhi8mIAzDMAxfTEAYhmEYvpiAMAzDMHwxAWEYhmH48v8BrYzHMd8b60EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.arange(len(attr_vals[0])-10000), np.abs(attr_vals[0][10000:]), alpha=0.5)\n",
    "plt.plot(np.arange(len(attr_vals[0])-10000), 2*targets_empty[10000:]/(10**2), alpha=0.5, color='r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a26ed58",
   "metadata": {},
   "source": [
    "### Mutagenesis Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "b93e9df6",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_ints_copy = deepcopy(X_ints)#.cpu()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "d850687a",
   "metadata": {},
   "outputs": [],
   "source": [
    "ids, test_empty = get_motif_loc(X_ints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "8580e484",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_ids = int(ids.shape[0] / 6)\n",
    "ids_groups = torch.tensor(np.array(np.split(ids, num_ids))).cuda()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f79cec5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "preds_modified_lst, new_ids_lst = [], []\n",
    "for i in range(num_ids):\n",
    "    print (i)\n",
    "    repl_tensor = torch.tensor([np.random.choice(np.arange(4)) for i in range(len(motif_str))]).cuda()\n",
    "    X_ints_copy = X_ints_copy.index_copy_(0, ids_groups[i], repl_tensor)\n",
    "    new_ids = get_true_targets(X_ints_copy)[-1]\n",
    "    X_modified = torch.transpose(F.one_hot(X_ints_copy), 1, 0).unsqueeze(0).type(torch.FloatTensor).cuda()\n",
    "    with torch.no_grad():\n",
    "        preds_modified = torch.sigmoid(trainer.model(X_modified)).round().flatten()#.detach().numpy()\n",
    "    preds_modified_lst.append(preds_modified)\n",
    "    new_ids_lst.append(new_ids)\n",
    "preds_modified_lst = torch.stack(preds_modified_lst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "9d732b59",
   "metadata": {},
   "outputs": [],
   "source": [
    "true_vals_ids = ids[0::6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "f55fce12",
   "metadata": {},
   "outputs": [],
   "source": [
    "# torch.save(true_vals_ids, 'true_vals_ids.pt')\n",
    "# torch.save(preds_modified_lst, 'preds_modified_lst.pt')\n",
    "# torch.save(new_ids_lst, 'new_ids_lst.pt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "99bfe7ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff_ids = np.diff(true_vals_ids)\n",
    "bins = np.arange(0, 131072, 128)\n",
    "bin_indices = np.digitize(true_vals_ids, bins)\n",
    "bin_indices_preds = np.digitize(true_vals_ids[np.where(diff_ids >= 160)], bins) \n",
    "# preds_modified_lst = preds_modified_lst#.cpu().detach().numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "2b3309e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "upd_ids, upd_ids_preds = [], []\n",
    "for i in range(len(new_ids_lst)):\n",
    "    bin_indices = np.digitize(new_ids_lst[i], bins)\n",
    "    upd_ids.append(bin_indices)\n",
    "    upd_ids_preds.append(np.where(preds_modified_lst[i])[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "f9a954a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "targets_all = []\n",
    "for i in range(len(upd_ids_preds)):\n",
    "    targets_empty = torch.zeros(1024)\n",
    "    targets_empty[upd_ids_preds[i]] = 1.\n",
    "    targets_all.append(targets_empty)\n",
    "targets_all = torch.stack(targets_all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "ce00d428",
   "metadata": {},
   "outputs": [],
   "source": [
    "targets_all_mean = torch.mean(targets_all, dim=0)\n",
    "targets_true = torch.zeros(1024)\n",
    "targets_true[bin_indices_preds] = 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "feec497e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgdklEQVR4nO3dbYxc133f8e9/Z3aXz6IssbFMUqEMsEmIAK0VQpHtIgjiABGVIEzfSYAr12hBCLBqO20RyM2LtK+KFkbqqlVFC5IMq7alOLKS0A5jxXHsOLX1RFlPpCiKS1IiV1ySKz7s887jvy/m7nJ2dnbn7sydmTtnfh+A3Jl7z517zr3n/u+55z6ZuyMiIuEa6HYGRESkvRToRUQCp0AvIhI4BXoRkcAp0IuIBC7b7QzUc/PNN/uuXbu6nQ0RkZ7xyiuvfODu2+qNS2Wg37VrF0eOHOl2NkREeoaZvbfSOHXdiIgEToFeRCRwCvQiIoFToBcRCZwCvYhI4GIFejO7y8xOmNmImT1YZ/wvm9nzZpYzs/+4lmlFRKS9GgZ6M8sADwP7gD3AvWa2pybZFeDzwJebmFZERNooTov+DmDE3U+7ex54GthfncDdL7n7y0BhrdP2nMunYH6i9d+5eAyKudXTXDkDs1dan1e13DSMn1g+j/lJ+GAk2XmtplyGsdcrfzutHcu1FaUinHsJLh3vdk6WilNHW7VaPSjMV/Iw9jqUS8nOd+J9mLq48vhLb0N+Jtl5dlGcQL8dOFf1fTQaFkfsac3sgJkdMbMj4+PjMX++C974Nhz5Wmu/MXUR3joEJw6vnu71p+HFr7Y2r1pv/BkcfRaK+aXz+PnX4c0/T3Zeq3n/FXj7MIy92rl5LmjHcm3Fuz+BkR/Csb+Euavdzk3F5Filjr7z/fbOZ+y1Sj14v84Nkif+upKHtw/Dez9Ndr4/fxKOPFF/XH4Wjv0FvPlMsvPsojiB3uoMi/u2ktjTuvuj7r7X3fdu21b3Lt70KMy1Nn0pCrK5qdbzslbz1yp/vaYFlZvubD4Ks9HfFpdlCKpbjkm3XJvVqTq6sP4X6kO1+cnrn/N1xreLR+sgiSP3lIgT6EeBnVXfdwDnY/5+K9OKiEgC4gT6l4HdZnabmQ0B9wCHYv5+K9OKiEgCGj7UzN2LZvYA8ByQAZ5w92Nmdn80/qCZfRg4AmwBymb2RWCPu0/Wm7ZNZRERkTpiPb3S3Q8Dh2uGHaz6fIFKt0ysaUVEpHN0Z6yISOAU6EVEAqdALyISOAV6EZHAKdCLiAROgV5EJHAK9CIigVOgFxEJnAK9iEjgFOhlKY/7YFIR6RUK9CIigVOgFxEJnAK9iEjgFOhFRAKnQC9L6WSsSHAU6EVEAqdALyISOAX6tQiqWyOksojIahToRUQCp0AvIhI4BXqpoS4dkdAo0IuIBE6Bfi1COhkbUllEZFUK9CIigVOgFxEJnAK9LKUuHZHgKNCLiAROgX5NQmrthlQWEVmNAr2ISOAU6EVEAqdALzXUpSN9LsALEmIFejO7y8xOmNmImT1YZ7yZ2UPR+DfM7PaqcX9oZsfM7KiZPWVm65IsgIiIrK5hoDezDPAwsA/YA9xrZntqku0Ddkf/DgCPRNNuBz4P7HX3XwUywD2J5V5ERBqK06K/Axhx99PungeeBvbXpNkPPOkVLwBbzeyWaFwWWG9mWWADcD6hvHdeSId0IZVFRFYVJ9BvB85VfR+NhjVM4+7vA18GzgJjwIS7/229mZjZATM7YmZHxsfH4+ZfREQaiBPorc6w2uZg3TRmdiOV1v5twEeAjWb26XozcfdH3X2vu+/dtm1bjGxJW6ilL30vvG0gTqAfBXZWfd/B8u6XldL8NnDG3cfdvQA8C3yi+eyKiMhaxQn0LwO7zew2MxuicjL1UE2aQ8B90dU3d1Lpohmj0mVzp5ltMDMDPgUcTzD/IiLSQLZRAncvmtkDwHNUrpp5wt2Pmdn90fiDwGHgbmAEmAU+G4170cyeAX4OFIFXgUfbUZDOCOmQLqSyiMhqGgZ6AHc/TCWYVw87WPXZgc+tMO2fAH/SQh5FRKQFujNWaqilL30uwAsSFOhFRAKnQC8iEjgF+rVI8JBuNl+k3M1DxAAPT0WkvlgnYyVZE/MF3hqd4KbyNP+025kRkeCpRd8Fc/kSANdmC13OiYj0AwV6WUpdOtL3wtsGFOhFRAKnQL8mIe3pQyqLiKxGgV5E6ro6m6dYVoMgBLrqRkSW+WAmx8kLU2z2SX6125mRlqlFLzXUghPIFyv1YLZQ6nJOIlbvlRdtEuAFCQr0ImkSYJBJhJZLSxTo1yKkyhZSWURkVQr0IiKBU6AXEQmcAr0sFXqXTujlC1UnT8YGSIFepMvm8iVOXJyiVHZ01dMKtINuiQK9SJcdvzDFlZk8l2dy3c6KBEqBfk1CalWEVBYRWY0CvYhI4BTo+1ndfs/AW/rq6+1NujO2JQr0ImkSYJBJhJZLSxToRUQCp0C/FiG1KkIqi4RP19G3RIFeRNJPDZOWKNCnVScqdr15BL9BhV4+aV14dUSBXiRVwgsy0n0K9CIigQs20M8XSvyPH7zDm6MTCf6qWlsiXaGTsS2JFejN7C4zO2FmI2b2YJ3xZmYPRePfMLPbq8ZtNbNnzOxtMztuZh9PsgArmZwvAPDa6LVOzE5E2in4c0ft1TDQm1kGeBjYB+wB7jWzPTXJ9gG7o38HgEeqxv1P4Pvu/svAPwOOJ5DvhowebwF0pGLrzliRZQKsI3Fa9HcAI+5+2t3zwNPA/po0+4EnveIFYKuZ3WJmW4DfAB4HcPe8u19LLvsi0hfUddOSOIF+O3Cu6vtoNCxOmo8C48DXzOxVM3vMzDbWm4mZHTCzI2Z2ZHx8PHYBGgpw7ywBU32VNogT6OvtSmtr40ppssDtwCPu/jFgBljWxw/g7o+6+15337tt27YY2eqCxDbCFGzMCigifSNOoB8FdlZ93wGcj5lmFBh19xej4c9QCfwiIvGpYdKSOIH+ZWC3md1mZkPAPcChmjSHgPuiq2/uBCbcfczdLwDnzOyXonSfAt5KKvMiItJYtlECdy+a2QPAc0AGeMLdj5nZ/dH4g8Bh4G5gBJgFPlv1E/8O+Ga0kzhdM65tev/cjR6B0B6hly9QHd2gw6sjDQM9gLsfphLMq4cdrPrswOdWmPY1YG/zWRTpJ+EFGem+YO+MXZDsZpPQr6Wi1ZyGPEhqpaKOSlKCD/QiEgDteFqiQC8iErhgA33vn4vVIxDaQi3D3qSXg7ck2EC/IMB1JiFTha1Py6UlwQf6ROnOWOkbqh8hUaAXkfTr/RtjuirYQG+qGCIiQMCBvvfpztj2CL180rrw6kjwgd6DD1wSFtXXurQdtyT4QC8i0u8U6NckpFZFSGWRxKWtBa1zbi1RoBcRCZwCfVrpztj2SFtLVdInwDqiQC8iErjgA314+2YJWoCtSem+YAN9Wy6rTOo307AxpyEPklqmJlJQgg30IiJSoUCfWroztj1CL5+0Lrw6EmygD29ViYg0J9hAL9Kb1ESR5AUf6JPtiQjoefSpyIOkl+pHSIIN9MF3NYv0E23QLQk20IuISEXwgT6N7YBY1/h37REIgVPLsDfp5eAtCTbQez8GMel9AQYZ6b5gA31bJLQRWho25jTkQdJL9SMoCvQByhVLvHD6MuWyNlYRCTnQpzjGtTtrPx35gOdPXeadS1NtnpOI9IJwA33Pa353kC9Wpi01atHrEQgidYRXR4IP9Ho5uPQW1VdJXqxAb2Z3mdkJMxsxswfrjDczeyga/4aZ3V4zPmNmr5rZ95LKeCOp3lxSsfNJQx4kvVQ/QtIw0JtZBngY2AfsAe41sz01yfYBu6N/B4BHasZ/ATjecm5FRGTN4rTo7wBG3P20u+eBp4H9NWn2A096xQvAVjO7BcDMdgC/CzyWYL5XlSuWmJgrdGp2IiKplo2RZjtwrur7KPDrMdJsB8aArwB/BGxebSZmdoDK0QC33nprjGyt7KkXz3J1tscDvV4O3h6p6DaTVAuwjsRp0de797h2SdRNY2a/B1xy91cazcTdH3X3ve6+d9u2bTGytbK0B/mFetTBm7pFpI/FCfSjwM6q7zuA8zHTfBL4fTN7l0qXz2+Z2Teazm0MEykP8iKrSklrMiXZkITECfQvA7vN7DYzGwLuAQ7VpDkE3BddfXMnMOHuY+7+JXff4e67oun+3t0/nWQBao1NzrXvxxOr/Z3Zimy1YwZtybIa1Y+gNOyjd/eimT0APAdkgCfc/ZiZ3R+NPwgcBu4GRoBZ4LPty3I42r0p6cFuIgLxTsbi7oepBPPqYQerPjvwuQa/8WPgx2vOYYDiBWC9HLw9Qi+ftC68OhLcnbHBxymRDtDRYFiCC/TfP3phyXcFfuktqrCSvOACfXsltBGmYu+ThjyISCco0IvIcmoHBEWBPq10Z2x7pOJoSlItwDoSfKCfzhUplMrdzoaISNcEH+gBnj91udtZEIknwNakdF9fBPp8MaEWfUIbYae2Zd0ZK83S5ZVh6YtAb3p6mIj0MQX61NKdse0RevmkdeHVkT4J9D0Z6Vumw2+R9PlfPzzJd14Z7eg8+yPQdzsDNfTCckk71dH2KZads1dmOzrP/gj0ibXoe+sxxatLQx4ktRTog9Ifgb7bGRAR6aL+CPQpi/QLbaWUZUtEAtUfgb4XQ6oegdAe6pKIpa+XUoB1pD8CfQ/G+ST05A5ORBLXH4E+Zjp355lXRhm5NA3AuSuzPP3SWUplX0iQTIbS0GJIQx5kuZSsF+vvNn1w+iLQr6Vhe+7KLN99/TwAP3jrImMT80zNF9qUMZF0Ssn+RhLSF4E+bhdGbeVe6PJRpRcJSB9u0P0R6Huyq1qPQGiP0MsnrQuvjvRHoI+ZLrzVKyLSL4E+ZpO+9rbv5VMldTK2My9CWf1ZN9qtpVNa1ote1hOSPgn08dKlZRNbzQunL/PQD092Oxsi0kP6I9C3OL0D5bIzXyglkZ2WPH/q8vXLPUVk7YI/D7VcfwT6uC36Vdb/T06O8+Tz73UuyOrO2Pbow41c1ijAOtIXgT5um766T/vi5PyScScuTAFQSqASpOJVgiKrCDDW9bW+CPTNXF75rRfPLp7EbdezubUtyTKKsNIG/RHoY6ZruI0lthGmYGNWQJFVqX6EpD8CfW/eMbUqvQFIROKKFejN7C4zO2FmI2b2YJ3xZmYPRePfMLPbo+E7zexHZnbczI6Z2ReSLkAcSVx1k6R4v6c7Y6V7wq4FYZeunoaB3swywMPAPmAPcK+Z7alJtg/YHf07ADwSDS8C/8HdfwW4E/hcnWnbLu5qrY1xV2byi5/TdlCgeCwiccVp0d8BjLj7aXfPA08D+2vS7Aee9IoXgK1mdou7j7n7zwHcfQo4DmxPMP+JWulO0jQG1RRmSVpUqWdas5K8OIF+O3Cu6vsoy4N1wzRmtgv4GPBivZmY2QEzO2JmR8bHx2NkK75W+7MdX9wIE+kbT8WeIw15kNTSExCCEifQ1+u0qI0Sq6Yxs03Ad4AvuvtkvZm4+6Puvtfd927bti1GtpK3YvxNYUyMs8NZ/Vk3ItIv4gT6UWBn1fcdwPm4acxskEqQ/6a7P9t8VpsXu48+xvDZfIliOZnmzqrd/g0CeftCeOA7h1QcTUlXNaoDAdaROIH+ZWC3md1mZkPAPcChmjSHgPuiq2/uBCbcfcwq1zU+Dhx39z9NNOdtsFIr2f36ydi3L0xx7Hzdg5KOilMXdWesiABkGyVw96KZPQA8B2SAJ9z9mJndH40/CBwG7gZGgFngs9HknwT+FfCmmb0WDftP7n440VI0ECcouq/c0VE7Zjbf2sPNkmgvqFtG2kn1KywNAz1AFJgP1ww7WPXZgc/Vme7/0fpl7AloXGn/z49PrTx1NHliL0xOw/NyAjw8DUJq1kun8pGW8oYtVqDvB/niyv3uqorSCapn0i598QiEVhtJvfq4AR1+i9TTf9tFXwT6ViUd5xPpo2/Xj3R6p9bxnWj/beTN6NG2TULCK3xfBPrwVlu81np/b6wi6dOt3oG+CPSNNFr410cHdDJWUiotKzakN6lJXwT6RnWpUGp0c1KvVcZey6+kjqpQW3Rrv9YngX71pdvoTtfE++iTeFxO6z8h0p/68CiiLwJ9I43e992VatGwOymZ267iDWun0E/+xterV3cFp43roVtruC8CfeM3BHYiqCZr9RwtvOu2EzkRkbh0MraL4rbok7rFN4k7bHv/zthuzz+lur5eKrxjzylOR3lDF1ygzw4sD8eNH1YX96qbZHTq5G7vnUQWCZu6bhLSzIJsHMhTGDBTmCURWZ2uummjRi3bcodb9PG0fslnw3zrzlhJi47Wg4YbRkdy0UnBBfq69aXFq2rSuNpT0pUrgVL9ao9udacGF+ib0ekWfaeuo9e22lu86n8Jk7puElJvj9lqF/zibya2lpK46ma13+iB3YCajOnWqfWjetARwQX6ZjS8vFJ1UUR6WHCBvpnzi436zboS6BN4OXjjmzN0Z6ykRQfXi14O3tuaveusUYu+UR9+N6QwSyICjF6d5ZX3rtQd163tti9eJdi4xd64RW9pePVtlViXV3YgHyKy1J8fGQXg137xQ1yezrF+KMOGoUqo1VU3CVjpccOtHqnlS2Wmc8VUvRx89azE2Cl1/ZCg2/NPqa6vl0jH8pGS8rbJk8+/x5PPv9ftbIQV6J87dqGp6RrV6Z+8M97U7644v0R/bZX5hL0NhUfrK0hz+dLiZ11emYDT4zN1hzdats30wa/UB5ec1k/GNp6F7oyVlOiTO2P1rJsENHvStJmpfvLOB+SLnXrC33LxiqqgFtfkfIH5QqlxQpGYyo2u8uigoAL9ShqdbG12B1Eolfn2y+e4MDHf1PSt0JMpk/X4P57hiZ+e6XY2UkO167orM3muzebXPF29uKLn0Sfotps3ril9/IW/NN3YxDzvX5vjRycurWl+7X85uO6MbUau0L0jNEjbIxB0Z+yCr//sXb7203fXPF2xToteXTcJqt2TNvMQxzgWum4ydZ6BnwY9sA1JlbNXZsmXuruzkeQUUrQugwz0pao9aWbAGkb6ZrvScsVKn269l520LIE7Y2PMpG2/3Foe2jm7dO/9qq/Q6C9dujO2jRckFOtc7q2rbhJUvTDjtLab7e+ei07edaNFH6e7Kd0hLTnzhRLfe+M8M7liU9On6Z3AaclLSrLRU2ZyRSZmC4vf1aJvs5IvbdE3fPFIk+tjJlcJ9IOZtS3GRO6X0oa46PjYJCcvTvPSmeYuea3uS63dWRx6/Tyvnbu2pt9zd0YuTTV11UWKLtSQVZy7MstTL52lWBXMH6s5oV+otzLVok/OWq+iiduir70zduFMfHWgH7k0zdWZRmfo4+fvg+kck/OFxgnrKJacH7x1kZFL03Wy0O2Iktz8F3+pyQOr6kPsP3v53JJxpy5N86O313ay/eSlab77+hivnrvaRF7KuHsKWva9cWdsO5ZVvlTiwuQ8zspHWN8/eoELE/NMzl9vGNTGneqdwEJjoTrWdPLyy1iB3szuMrMTZjZiZg/WGW9m9lA0/g0zuz3utO1QvQCNxjFtpUcnNDIxVwnAC+cE3J3vvn6eb7yw9Jbn+UJpyUovRYcQ1qDLp1guc/LSNH/16vvLxsXZmb105gpH35/gu6+fT0Hg6ICYRSyWy0uWX6HqkG5hncLSDXUtpqOjgsm55V1J7s5sfuUupmKpzDdfPMs3Xjzb1LyT0q4rkCbnC7xw+nIi9fG1c9f4yt+d5FsvJbusjrx7lTMfzHD2yixf+buTi8NHLk3z9oVJ4Pr5ua//7N0Vy5Krus/mR29f4qUzV5bEmlwH78NpGOjNLAM8DOwD9gD3mtmemmT7gN3RvwPAI2uYNjGf+cQuYPnh7xujE5y/NrfidNPzzfXtTkXTzRdKTMwVGL1amUftZVWP/PgUf1EVrGfzS0/iXp3JL1aWhc9Hz09w9vIswJJWw8I0UzV5LpTKiy3/+WgjrQ5msyk7yTdfLDNyaZrjY5McfX9icXi+WI51BFMq+2J/6NJza86VaBlWf14cj/Pyu1d55+IUxVKZidnCikFtPsaGWCiVl+wcatXm4WenLvPVfzjN2cuzfP/ohWU3aRXdGZ/K8cFUruG8mzUxW+C5YxdWvOHvxIUp3ozWyXyh1FRfs7szNjG3pNHl7jz+j2d4/tRlxqdyXJvNE/enF+rFxGyBS5PzHPyHU4tHWpcmk11W16L1WXt/zHdfP8/fvHmBidnCkoD9/rW5usH+uaNLH8ny05EPuDxzPa9//eYY56/NrVp/kmKN9qxm9nHgP7v770TfvwTg7v+1Ks1XgR+7+1PR9xPAbwK7Gk1bz969e/3IkSNrLszszx7jhZELbBweZCa3fOFtHK7/sM5csdyw9bZ+MEMul2O4VOkGmRvcumr66nkt9PsuDCvlZrFCJYiX199ErlhiwIzBzAC5YqlyXqFUZLg4tTivhWkXfmsoM8BgdmDZPDYMZZnLl5Z1R60bzJAZMNblK90J+exG3DIMFyaXpMtnN1EeGFy1bEkozU9hxRzFgWEKmfXLyrfSulpwvbwZiiUnXyqTGRhgwCrBNxt1pxVLZQYzAwxFy8rKJUqzlWWwsA6Hs5nFFlr1vMvui1fBrJSf6uVu0QFaoVium5+hzMCy1nw2M8BQYRovRd19Q5uY82ysZdCshTwPZzNkM8uPKmdyRbLlHIOlSsOluv7FNV8oUyqXK3VuMANUdsgL5V9Y5oOlObZkixQHhilmNyz5jYW6CjBVzJDPLB1fbS35My8zXIh2ZEM3Ls/73BxDpcrjVBpt54u/SeNzgVDp5q2341yoPza0gTv+5edjzXNZHsxecfe99cbFWTrbgeqOy1Hg12Ok2R5z2oVMHqByNMCtt94aI1vLrd+6jZ07hti+dT0z+SL5YplcsczJi1PcuHFoxatjssDwYIZ8sUy+WGbAYHKuwPYb1zM1XyQzYNy0cbiy0V49x7WBrUzPw00bh7l50xCj1+aYzRUpufPhLetwll7i6ZkC6wczZAcrwSa7GbbOjzKxfifuzvRUjps2DWMG01M5tmwaZjpXJJvPUh7azA3rrt8A9qEtxsCALass69aVmS+WGFw/yBDGfKHE+qEMw9kBrs0WGI7mXfKbWD83RnHD9kpe5gcoDG7BLcvw/DjFDbd05MRNdvM25i+OUN56K1mM7HBmsRy5Yons+tV3Nhs3OJPzBQY3DjEIzEzluGHzMFBp4d2wZfnnBaVymcyGG9mQXcfMfJF1GwexQpm5fIlN67JLgt/cXLTusvWXyrr1lekGN1zPb5bKetxWJz+lmTxD2QE2DWe5MDHPh29YB/wC+UsnwTIM3bQTKzmZARiw9lzNtWVTlKcbhuuOv2FzZfzW3CizG29l/dDK5V/JJiqt8KGa6QrTOUoluOGGYWZn8mxcN8jQ/CjFjTuXBaMiN7NhpvLI33L2I5SLJQxjy4Ys45N5Nq3L4jhWVX/iGsgP45YhO7h52bjNm42Bifcob9mBF52BgQEuT+e4ceMQALO5Ih/aNIxROYdWLjs3bx5ePALYMJRhNl9i2+bKNp0rltk8nF2MC4OZAcruXJiYZ9O6LLlCebH+ZAbXr6kcccUJ9PVqW+2ua6U0caatDHR/FHgUKi36GPlaxvbsZ3fUMVS97/9oMz+2Bs3tlkQkrl/pdgZ6XJxAPwrsrPq+AzgfM81QjGlFRKSN4hyPvQzsNrPbzGwIuAc4VJPmEHBfdPXNncCEu4/FnFZERNqoYYve3Ytm9gDwHJABnnD3Y2Z2fzT+IHAYuBsYAWaBz642bVtKIiIidTW86qYbmr3qRkSkX6121U2Qd8aKiMh1CvQiIoFToBcRCZwCvYhI4FJ5MtbMxoH3Gias72bggwSzk2b9VFbor/L2U1mhv8rbrrL+ortvqzcilYG+FWZ2ZKUzz6Hpp7JCf5W3n8oK/VXebpRVXTciIoFToBcRCVyIgf7Rbmegg/qprNBf5e2nskJ/lbfjZQ2uj15ERJYKsUUvIiJVFOhFRAIXTKDvxkvI28nMdprZj8zsuJkdM7MvRMM/ZGY/MLOT0d8bq6b5UlT+E2b2O93LffPMLGNmr5rZ96LvQZbXzLaa2TNm9na0jj8ealkBzOwPo3p81MyeMrN1IZXXzJ4ws0tmdrRq2JrLZ2a/ZmZvRuMeMkvoNWMLL1Hu5X9UHoF8isrLpIaA14E93c5Xi2W6Bbg9+rwZeIfKC9b/O/BgNPxB4L9Fn/dE5R4GbouWR6bb5Wii3P8e+Bbwveh7kOUFvg782+jzELA14LJuB84A66Pv3wb+dUjlBX4DuB04WjVszeUDXgI+TuXtfH8D7Esif6G06O8ARtz9tLvngaeB/V3OU0vcfczdfx59ngKOU9lg9lMJEkR//yD6vB942t1z7n6GyrsB7uhopltkZjuA3wUeqxocXHnNbAuVwPA4gLvn3f0aAZa1ShZYb2ZZKm/6PE9A5XX3nwBXagavqXxmdguwxd2f90rUf7JqmpaEEuhXejl5EMxsF/Ax4EXgF7zy9i6iv/8kShbCMvgK8EdA9ZvPQyzvR4Fx4GtRN9VjZraRMMuKu78PfBk4C4xReQPd3xJoeaustXzbo8+1w1sWSqCP/RLyXmNmm4DvAF9098nVktYZ1jPLwMx+D7jk7q/EnaTOsF4pb5bKYf4j7v4xYIbKof1KermsRH3T+6l0U3wE2Ghmn15tkjrDeqa8MaxUvraVO5RAH+cF5j3HzAapBPlvuvuz0eCL0SEe0d9L0fBeXwafBH7fzN6l0vX2W2b2DcIs7ygw6u4vRt+foRL4QywrwG8DZ9x93N0LwLPAJwi3vAvWWr7R6HPt8JaFEuiDewl5dLb9ceC4u/9p1ahDwGeiz58B/qpq+D1mNmxmtwG7qZzY6Qnu/iV33+Huu6isv793908TYHnd/QJwzsx+KRr0KeAtAixr5Cxwp5ltiOr1p6iccwq1vAvWVL6oe2fKzO6MltN9VdO0pttnqxM86303lStTTgF/3O38JFCef0HlsO0N4LXo393ATcAPgZPR3w9VTfPHUflPkNDZ+i6V/Te5ftVNkOUF/jlwJFq/fwncGGpZo/z/F+Bt4Cjwf6lccRJMeYGnqJx/KFBpmf+bZsoH7I2W0SngfxM9vaDVf3oEgohI4ELpuhERkRUo0IuIBE6BXkQkcAr0IiKBU6AXEQmcAr2ISOAU6EVEAvf/AXZSC1idRtAVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.arange(1024), targets_all_mean, alpha=0.5)\n",
    "plt.plot(np.arange(1024), targets_true/10, alpha=0.5)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "basenji_pytorch",
   "language": "python",
   "name": "basenji_pytorch"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}