nTimesCoverage.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This notebook contains code that analyzes our designs using nTimesCoverage (Figure 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pickle\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "import random\n",
    "import gzip\n",
    "import math\n",
    "\n",
    "from copy import deepcopy\n",
    "from matplotlib import pyplot as plt\n",
    "import Levenshtein as lev\n",
    "\n",
    "import torch\n",
    "from torch.nn import parallel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Import Data</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"./data.pkl\", 'rb') as fin:\n",
    "    mhc1_data, mhc2_data = pickle.load(fin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "779it [00:00, 1128.42it/s]\n",
      "1200it [00:01, 639.29it/s]\n",
      "440it [00:00, 4539.60it/s]\n",
      "333it [00:00, 2423.57it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([1100, 234])\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "537it [00:00, 2199.85it/s]\n",
      "920it [00:01, 890.77it/s] \n",
      "502it [00:00, 1232.79it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([4195, 281])\n"
     ]
    }
   ],
   "source": [
    "# Script for processing data\n",
    "\n",
    "def makeProductDistribution(population, dataset):\n",
    "    _,alleles,_ = dataset\n",
    "    intmap = {}\n",
    "    for i,allele in enumerate(alleles):\n",
    "        intmap[allele] = i\n",
    "    default = len(alleles)\n",
    "    population = [(frozenset([intmap.get(x[1], default) for x in genotype]),p) for genotype, p in population]\n",
    "    diploidDistribution = {}\n",
    "    for i,(haplotype1,p1) in tqdm(enumerate(population), position = 0, leave = True):\n",
    "        if haplotype1 not in diploidDistribution:\n",
    "            diploidDistribution[haplotype1] = 0\n",
    "        diploidDistribution[haplotype1] += p1 ** 2\n",
    "        for (haplotype2,p2) in population[:i]:\n",
    "            diplotype = haplotype1.union(haplotype2)\n",
    "            if diplotype not in diploidDistribution:\n",
    "                diploidDistribution[diplotype] = 0\n",
    "            diploidDistribution[diplotype] += 2 * p1 * p2\n",
    "    return diploidDistribution\n",
    "    \n",
    "def makeDiploidDistribution(populations, dataset):\n",
    "    default = len(dataset[1])\n",
    "    distributions = [makeProductDistribution(population, dataset) for population in populations]\n",
    "    for distribution in distributions[1:]:\n",
    "        for genotype in distribution:\n",
    "            if genotype not in distributions[0]:\n",
    "                distributions[0][genotype] = distribution[genotype]\n",
    "            else:\n",
    "                distributions[0][genotype] += distribution[genotype]\n",
    "    distribution = distributions[0]\n",
    "    indexes = []\n",
    "    weights = []\n",
    "    for genotype in distribution:\n",
    "        weights.append( distribution[genotype]/3 )\n",
    "        genotype = tuple(genotype)\n",
    "        if len(genotype) != 6:\n",
    "            genotype = genotype + ( (default,) * (6-len(genotype)) )\n",
    "        indexes.append(genotype)\n",
    "    indexes = np.array(indexes, dtype = np.int16)\n",
    "    weights = np.array(weights)\n",
    "    return indexes, weights\n",
    "    \n",
    "def reformatData(data):\n",
    "    _, binary_data, _, population = data\n",
    "    \n",
    "    binary_distribution = makeDiploidDistribution(population, binary_data)\n",
    "    binary_data = torch.tensor(binary_data[0], dtype = torch.float32).cuda(), binary_data[2],\\\n",
    "        binary_distribution[0], torch.tensor(binary_distribution[1], dtype = torch.float32).cuda()\n",
    "    \n",
    "    print (binary_data[0].size())\n",
    "    return binary_data\n",
    "\n",
    "inputs1_binarized = reformatData(mhc1_data)\n",
    "inputs2_binarized = reformatData(mhc2_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Implementation of nTimesCoverage</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Utility functions\n",
    "\n",
    "def dedupe_designs(seqs):\n",
    "    seen = set()\n",
    "    seqs2 = []\n",
    "    for seq in seqs:\n",
    "        if seq not in seen: seqs2.append(seq)\n",
    "        seen.add(seq)\n",
    "    return seqs2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#candidates: [candidate, 1 - pMHC hit probability]\n",
    "#columnIndex: [diplotype, allele in diplotype]\n",
    "#columnWeights: [diplotype]\n",
    "#distributions: [dummy, diplotype, distribution]\n",
    "\n",
    "def nTimesCoverage(candidates, seqs, columnIndex, columnWeights, design, thresholds):\n",
    "    design = set(design)\n",
    "    designIndex = [i for i,seq in enumerate(seqs) if seq in design]\n",
    "    \n",
    "    candidates_in_design = candidates[designIndex, :]\n",
    "    numberOfHits = torch.sum( torch.sum(candidates_in_design[:, columnIndex], dim = 2).unsqueeze(2), dim = 0 )\n",
    "    passed = (numberOfHits >= thresholds.view(1,-1)).float()\n",
    "    return torch.sum( passed * columnWeights.view(-1,1), dim = 0 ).cpu().numpy()\n",
    "\n",
    "def nTimesCoverage_trajectory(candidates, seqs, columnIndex, columnWeights, design, thresholds):\n",
    "    design = dedupe_designs(design)\n",
    "    seqToIndex = {}\n",
    "    for i, seq in enumerate(seqs):\n",
    "        seqToIndex[seq] = i\n",
    "        \n",
    "    designIndex = [seqToIndex[seq] for seq in design]\n",
    "    candidates_in_design = candidates[designIndex, :]\n",
    "    numberOfHits = torch.cumsum( torch.sum(candidates_in_design[:, columnIndex], dim = 2).unsqueeze(2), dim = 0 )\n",
    "    stack = torch.stack([torch.sum((row >= thresholds.view(1,-1)).float() * columnWeights.view(-1,1), dim = 0)\n",
    "         for row in numberOfHits]).cpu().numpy()\n",
    "    return stack"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Analysis</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "credence1 = []\n",
    "binarized1 = []\n",
    "credence2 = []\n",
    "binarized2 = []\n",
    "for threshold in range(1, 21):\n",
    "    with open(\"./outputs/mhc1_binarized_threshold_{}.pkl\".format(threshold), 'rb') as fin:\n",
    "        binarized1.append(pickle.load(fin))\n",
    "    with open(\"./outputs/mhc2_binarized_threshold_{}.pkl\".format(threshold), 'rb') as fin:\n",
    "        binarized2.append(pickle.load(fin))\n",
    "        \n",
    "with open(\"./benchmarks.pkl\", 'rb') as fin:\n",
    "    benchmark_trajectories, benchmark_designs = pickle.load(fin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold_20 = torch.tensor(np.arange(1,21,1)).cuda()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 29/29 [00:00<00:00, 73.76it/s]\n",
      "100%|██████████| 8/8 [00:00<00:00, 50.67it/s]\n",
      "100%|██████████| 8/8 [00:00<00:00, 86.06it/s]\n",
      "100%|██████████| 101/101 [00:04<00:00, 22.65it/s]\n",
      "100%|██████████| 7/7 [00:00<00:00, 88.22it/s]\n",
      "100%|██████████| 7/7 [00:00<00:00, 54.30it/s]\n",
      "100%|██████████| 96/96 [00:01<00:00, 66.75it/s] \n",
      "20it [00:03,  4.86it/s]\n",
      "20it [00:01, 13.31it/s]\n"
     ]
    }
   ],
   "source": [
    "# Calculate nTimesCoverage over all benchmarks\n",
    "\n",
    "benchmark_designs_scores = {}\n",
    "for design in tqdm(benchmark_designs, position = 0, leave = True):\n",
    "    if \"MHC1\" in design:\n",
    "        data = inputs1_binarized\n",
    "    elif \"MHC2\" in design:\n",
    "        data = inputs2_binarized\n",
    "    else:\n",
    "        raise Exception(\"Unknown MHC Class\")\n",
    "    scores = nTimesCoverage(*data, benchmark_designs[design], threshold_20)\n",
    "    benchmark_designs_scores[design] = scores\n",
    "    \n",
    "benchmark_trajectories_scores = {}\n",
    "for design in benchmark_trajectories:\n",
    "    if \"MHC1\" in design:\n",
    "        data = inputs1_binarized\n",
    "    elif \"MHC2\" in design:\n",
    "        data = inputs2_binarized\n",
    "    else:\n",
    "        raise Exception(\"Unknown MHC Class\")\n",
    "        \n",
    "    trajectory = benchmark_trajectories[design]\n",
    "    trajectoryScores = {}\n",
    "    for designSize in tqdm(trajectory, position = 0, leave = True):\n",
    "        scores = nTimesCoverage(*data, trajectory[designSize], threshold_20)\n",
    "        trajectoryScores[designSize] = scores\n",
    "    benchmark_trajectories_scores[design] = trajectoryScores\n",
    "    \n",
    "# Calculate nTimesCoverage over our designs\n",
    "\n",
    "def getAllTrajectories(designs, data):\n",
    "    trajectories = [[] for j in range(20)]\n",
    "    with tqdm(position = 0, leave = True) as pbar:\n",
    "        for design_threshold, design in enumerate(designs, 1):\n",
    "            trajectories[design_threshold-1] = nTimesCoverage_trajectory(*data,design, threshold_20)\n",
    "            pbar.update(1)\n",
    "    return trajectories\n",
    "\n",
    "mhc1_design_scores = getAllTrajectories(binarized1, inputs1_binarized)\n",
    "mhc2_design_scores = getAllTrajectories(binarized2, inputs2_binarized)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAACpCAYAAADDeP+kAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2de/wc0/nH3x8RhIgQJQSJn7qVChGUal1/qlp1r9JW0FZbLVq9uLVFb2iL0lbboMS1LeLaVqWtVGldksiN6E/rLiLELQgSnt8f56xMxl5md2d3dvf7vF+vee3MOc885+zsfPbMnDPzHJkZjuM4juP0JksVXQHHcRzHcVqHN/SO4ziO08N4Q+84juM4PYw39I7jOI7Tw3hD7ziO4zg9jDf0juM4jtPDeEPvOI7TB5B0oqQLcvJlkt6dh688kHSopNuLrken4g29kwuSHpG0azXBSZoo6TVJL0t6VtJ4SWu0u66OUzRRL29IWjWVPjU2oiPyLtPMfmhmn83bL4Ck0ZJukvS8pBck3S/pB5JWbkV5Tn14Q++0my+b2UBgA2AwcHbB9XGcongYOKi0Iem9wIBGHElaOq9KNVD2dsBE4A5gIzMbDOwOLAJGVtinsPr2RbyhdwrBzJ4DrgE2LboujlMQlwKHJLbHAJeUNiR9RNK9kl6S9LikUxJ5I+Kd/2ckPQb8LaYfIulRSfMkfbvU0xbzTpF0WWr/MZIeiz1sJyX8by3pX/Hu/ClJP5e0TIXv8SPgIjM7zcyeBjCzx8zsZDObGP0dKukOSWdLeg44JaYfLmlW7An4s6ThiTpsJGmCpOck/VvSxxN5QyTdEI/N3cB6ibxfSDozWUFJN0r6Ss1fpEfxht4phNhluR9wb9F1cZyCuBMYJGljSf2AA4HLEvmvEC4EBgMfAb4oae+Ujx2AjYEPSXoPcB7wSWANYCVgWI06bA9sCOwCfEfSxjH9TeCrwKrAtjH/yPTOklaI+ddk+L7bAA8BqwE/iN/lRGBf4F3AP4ArE34nAFdE+4OA8yRtEn39Angtfs/D41JiHHCQpKWir1Vj/a/MUMeexBt6p92cK+kFYBrwFHBswfVxnCIp3dX/L/AA8GQpw8wmmtkMM3vLzKYTGqodUvufYmavmNkCYH/gRjO73czeAL4D1JrM5FQzW2Bm0wiaHBnLnmxmd5rZIjN7BPh1mbIBVia0I3NKCZJ+FHsCXpH0rYTtbDP7WfS5APg8cJqZzTKzRcAPgc3jXf1HgUfM7KJoP4VwMbF/vCjaD/hO/O4zCY176bjdDbxIaNwBPgFMLPU29EW8oXfazdFmNtjMhpnZJ83smaIr5DgFcilwMHAoiW57AEnbSLpV0jOSXgS+QLjDTvJ4Yn3N5LaZvQrMq1H+nMT6q8DAWPYG8eG6OZJeIjTC6bIBngfeItxZl8r9ZhynvxZIjsU/ntp3OHBOvCh4AXgOEKEXYjiwTSkv5n8SGEq4+1865e/RlO9xwKfi+qcIx7nP4g294zhOQZjZo4SH8vYAxqeyrwBuANY2s5WAXxEawiVcJNafAtYqbUgaAAxpsGq/JPQwrG9mgwhd7OmyMbNXgLsI3e+1SPcuPA58Pl74l5YBZvbPmPf3VN5AM/si8AzhQb+1E77WSfm+DNhL0kjC0MZ1GerXs3hD77QCSVouuRRdIcfpYD4D7BwbzSQrAs+Z2WuStibc+VfjamBPSdvFB+dOpUzjnJEVgZeAlyVtBHyxiu03gcMlHS9pNQBJawHr1ijjV8AJpXF3SStJOiDm3QRsIOnTkvrHZStJG5vZm4SLolMkLR+fTRiTdGxmTwD3EO7kr4lDBX0Wb+idVrAdsCC5+Os0jlMeM/uvmU0qk3Uk8F1J8wnj7b+v4ec+4Cjgt4S7+/nAXOD1Bqr1dcKFxXzgfOB3Vcq9HdgZ+CDwf7Gb/WbCK3c/q7LftcAZwG/j8MBM4MMxbz6wG2F8fTZhiOEMYNm4+5cJwwxzgIuBi8oUMQ54L3282x5AZrWe1XAcx3G6DUkDgRcI3e8PF12fdiPpg4Qu/BFm9lbR9SkSv6N3HMfpESTtGbuzVwB+AswAHim2Vu1HUn/gGOCCvt7Igzf0juM4vcRehK7u2cD6wCesj3XbxlgALxDeBPhpwdXpCLzr3nEcx3F6GL+jdxzHcZwexht6x3Ecx+lheu6Vp8GDB9vw4cNr2k2fPr0NtXGcjuZZM3tX0ZWoRhY9L1y4kFmzZrWpRo7TsVTUc8c39JJ+Q4h7PNfMas50Nnz4cKZNm5bFbw61c5yuJh02tOW0Qs+zZ89m2LBac7c4Ts9TUc/d0HV/MWFuY8dxup+LcT07TlvJ/Y4+Tg04kjDBwgLgvmZmDTKz2ySNyKd2juMUievZcdpPbg29pPWA44BdgQcJEw8sR4hX/CphmsNxHrzAcRzHcdpHnnf03wfOI8xGtMTL+XGig4OBT5OYNzgvJB0BHAH4WJ3jdDmuZ8fJl9zG6M3sIOAOYNsyeXPN7KdmlnsjH/2PNbPRZjZ6yJBGZ2V0HKcTcD07Tr7k+jBe7JY/M0+fjuM4juM0Tiueur9F0n7K6f01SVcC/wI2lPSEpM/k4ddxnPbjenac9tOK9+iPBVYA3pS0ABBgZjaoEWdxSCAzCxcuZPbs2TXtJk+eDMDQoUOr2vkYoePkRyv0PGfOHCZPnlxTy3PmzAFqa37OnDlsueWW9VTTcTqa3Bt6M1sxb5+O4ziO4zRGK96jF/BJYF0z+56ktYE1zOzuBv2tDVwCDAXeAsaa2TmV7Pv378+aa66Z2X89to7jNE69Wob69Jy3neP0ChUbekmjqu1oZlMqZJ1HEPHOwPeAl4FfAFs1WMdFwNfMbIqkFYHJkiaY2f0N+nMcpxhcy45TANXu6EtPzy8HjAamEcbbNwPuAravsN82ZjZK0r0AZva8pGUaraCZPQU8FdfnS5oFDAP8z8FxugjXsuMUQ8Wn7s1sJzPbiRAof1R8r3VLYAvgP1V8LpTUDzAASe8i3OE3TQyduQXhQsNxnC7Ftew47SPL63UbmdmM0oaZzQQ2r2J/LnAtsJqkHwC3A6c1VUtA0kDgGuArZvZSKu8ISZMkTZo3b16zRTmO00KqaTnmu54dJ0eyPIw3S9IFwGWEu/RPARUnfzazyyVNBnYhdPXvbWZNTRYtqT/hj+FyMxtfpsyxwFiAkSNHWjrfcZzOoJaWwfXsOHmTpaE/DPgicEzcvg34ZSVjSZea2aeBB8qk1U18iv9CYJaZndWID8dxise17DjFULOhN7PXgLPjkoVNkhtxvL6Z6BPvJ0yGM0PS1Jh2opn9sQmfjuO0H9ey4xRAtdfrZhAfqCuHmW2Wsj8BOBEYIOklQrc9wBvEbrhGMLPbE75qkjUyXilKVi2yRtBrRdStvKN9ZbFrhc9esWtV2XlHYct63rQ7+lu9WobskfGyUI9dPdrzKHpOp1Ptjv6j9Tgys9OA0ySdZmYnNFctx3Ecx3HyoGJDb2aPltYlrc7igDd3m9ncKj5PkvQpcoqMVy+tioxXRNStoUOHFlq/osrudLtW+cyLes6bTqfIyHi9cgwdp+brdZI+DtwNHAB8HLhL0v5VdvkFYU76g+N2KTJeQ0haTtLdkqZJuk/SqY36chynOFzLjlMMWZ66PwnYqnQXHwPg/AW4uoJ9rpHxgNeBnc3s5fhqzu2S/mRmdzbh03Gc9uNadpwCyNLQL5Xqqp9H9Z6AXCPjmZkRegUA+sfF3611nC7Dtew4xZAlMt7Nkv4s6VBJhwJ/AKq9DlMuMt4Pm6mkpH7xdZy5wAQzuyuV75G0HKcLqKXlaON6dpwcqdnQm9k3CK/HbQaMJEwteVwV+8uBbxLC3j5FiIx3VTOVNLM3zWxzYC1ga0mbpvLHxlj8o4cMGdJMUY7jtJBaWo42rmfHyZFM89Gb2TWEsJVZeRr4R/Q/QNKoKtPaZsbMXpA0EdgdmNmsP8dxisG17Djto2ZDL2lf4AxgNUKwCxGG2wZVsP8ecCjwXxaPvxlhfvq6iWP8C+MfwwBg11gfx3G6CNey4xRDljv6HwF71jExzceB9czsjcartQRrAOPiA35LAb83s5sqGecdGa8ou6y2RZXbF+2KLjsreZ43OVOXlqG4yHh52jlO0WRp6J+uc/a5mcBgwsM2TWNm0wnzVudKK8KOOk4jtCLUcSfSKi07+ZDneViPbZHhovsK1WLd7xtXJ0n6HXAd4T1YACpNMUl4CO9eSTNT9h9rvrq1qTcyXl6UTlKPjNf7dkWW3deiunlkvPbRiv+bemz72vFuJ9Xu6PdMrL8K7JbYNqBSQz+OMO42gyben08Tu/smAU+aWV1x+B3H6Rxcy47TXqrFuj+sQZ/Pmtm5De5bjWOAWUDZhwAdx+kaXMuO00ayxLr/kaRBkvpL+qukZ+OkNZWYLOk0SdtKGlVamqmkpLWAjwAXNOPHcZxicS07TvvJ8jDebmb2TUn7AE8QJre5Fbisgn3pYZv3JdIafr0u8lNCEJ4Vy2VKOgI4AmDYsGFNFOM4ToupqmVwPTtO3mRp6PvHzz2AK83sOUkVjc1spzwqVkLSR4G5ZjZZ0o4VyhxLiN7HyJEjPXa243QgWbQMrmfHyZssDf2Nkh4AFgBHxqAXr1XbQdJHgE2A5UppZvbdBuv4fuBjkvaI/gZJuszMqg0fOI7TebiWHacAssS6P54wv/xoM1tIeAJ/r0r2kn4FHAgcRYiidwAwvNEKmtkJZraWmY0APgH8zf8YHKf7cC07TjFkCYG7PPAlYB3CuNmawIZApYhW25nZZpKmm9mpks6k8qt4uZM1Ml7eeGS83rcrsuy+GtXNI+O1j7wjKvrx7hyydN1fBEwGtovbTwBXUbmhL3XrvyppTcL89es2U8kSZjYRmFjNZtasWR3/AE89EagcpxfJomUoTs/1RLp88skna9q0IqpbX4momKQV37nddq0qu9r5laWhX8/MDpR0EICZLVC1p/HCmP5g4MfAFMIT9+dnKKerqTcyXqfbFVl2p9sVWbZHGWsPrYoSlydFRs7si1rphmNTiSwN/RtxpikDkLQeidC2SSQtBfzVzF4ArpF0E7Ccmb3YTCUlPQLMB94EFpnZ6Gb8OY5TDK5lx2k/WRr6k4GbgbUlXU54cvbQcoZm9lYck982br9OhYuCBtjJzJ7NyZfjOMXhWnacNlK1oY9d9A8A+xIC4Ag4poZIb5G0HzDezPwdWMdxHMcpkKoNvZmZpOvMbEvgDxl9HgusACyS9Brh4sDMrJm41ka4gDDg1zGgxtskI2k5jtPRVNUyuJ4dJ2+ydN3fKWkrM7sni0Mzqxjasgneb2azJa0GTJD0gJndlijz7Uha8Q/EcZzOpKqWwfXsOHmTpaHfCfi8pEeBV1h8h75ZpR0krQysz5KR8W6rZF8LM5sdP+dKuhbYGmjYn+M4xeBadpz2k6Wh/3A9DiV9ljAN5VrAVMLY/r9ocFIbSSsAS5nZ/Li+G9BoOF3HcQrCtew4xVCzoTezRyX1A1bPYk9o5LcC7jSznSRtBJzaRB1XB66Nr+4vDVxhZjc34a8ldHp0rm6IaNXpdkWW3SNRxrpCy1BclLh6KCJypmulebsiys4SAvcowit2TwNvxWQDKnXdv2Zmr0lC0rJm9oCkDTPVtgxm9hAwstH9O5GsEbLqic7lOJ1OL2oZ8p9K1yNnOnmT9Q59QzObl9HnEzEy3nWEh22eB9offL7N1BsZL6vPXora1Ct2RZbd6dHaeoV6tFdk2Z1uV2TZnW7XKp/lyNLQPw5kjmxnZvvE1VMk3QqsRAi40zDxwuECYFNCb8LhZvavZnw6jtN+XMuO034qNvSSjo2rDwETJf2BRJQ7Mzuryr6jgO0JQr7DzN5osp7nADeb2f6SlgGWb9Kf4zjF4Fp2nDZTbT76FePyGDABWCaRNrDSTpK+A4wDhgCrAhdJ+lajFZQ0CPggcCGAmb0RY+k7jtNFuJYdpxgq3tGb2akAkg4ws6uSeZIOqOLzIGALM3st2p5OmMXu+w3W8X+AZwgXDCMJU+YeY2avJOrjkbQcp/OpqWVwPTtO3lS7oy9xQsa0Eo+QCJQDLAv8t446pVkaGAX80sy2IATtOT5pYGZjzWy0z4TlOB1NTS2D69lx8qbaGP2HgT2AYZLOTWQNAhZV8fk6cJ+kCYQx+v8Fbi/5MLOj66zjE8ATZnZX3L6aMn8OjuN0PK5lxymAak/dzwYmAR8jdLGVmA98tcp+18alxMRGKwdgZnMkPS5pQzP7N7ALcH8zPh3HaT+uZccphmpj9NOAaZKuMLOFWR2a2bhcarYkRwGXx6d0HwIOa0EZTVFU9KReiNrULXZFlt0jkfGgC7QM+UfGK6ps10rn2RVRdpYQuJkb+VZhZlOBPjde5xH0nF6jr2rZcYokS8AcJwOtiIxXT9ndELWpV+yKLNsj47UHj4yXj12RZXe6Xat8lqPmU/eSNm2qhCaRtKGkqYnlJUlfKbJOjuPUj2vZcYohyx39r+J42sWE2aaqBriQtAHwDWB40r+ZNTRNbXxoZ/Poux/wJEs+7Oc4ThfgWnacYsgyRr+9pPWBw4FJku4GLjKzCRV2uQr4FXA+8GZuNQ3sAvzXzB7N2a/jOO3Ftew4bSLTGL2ZPRjD2E4CzgW2UJhU+kQzG58yX2Rmv8y5niU+AVzZIt+O47QP17LjtIksY/SbSTobmAXsDOxpZhvH9bPL7HKjpCMlrSFpldLSbEXj8MHHCD0G6bwjJE2SNKnZchzHaS3VtBzzXc+OkyNZ7uh/TuiGP9HMFpQSzWx2hclqxsTPbyTSjBDnuhk+DEwxs6fTGWY2FhgLIMmaLMdxnNZSUcvgenacvMnS0I83s0uTCZKOMbNz0ukAZrZubrVbkoPwrj7H6QVcy47TRrI09IcAP02lHUqYV/odSNq3TPKLwAwzm1tX7Rb7XJ4QM//zjezvOE5n4Fp2nPZTbVKbg4CDgXUl3ZDIWhGYV8XnZ4BtgVvj9o7AncAGkr5brhegFmb2KmF++45l2LBhhZWdd7jOTg8h2ZdCVxZtlzfdoGUoVs9Zy84SEdO10nl2RZRd7Y7+n8BTwKrAmYn0+cD0Kvu9BWxcGn+TtDrwS2Ab4Dag7obecRzHcZzGqDapzaPAo4S783oYkXrIZi6wgZk9J6mhuPmSvgp8lvBQ3wzgMDN7rRFfvYiHwG2vXZFld3sIXNdyfnRDqNwiy+50u1b5LEfF1+sk3R4/58dQlaVlvqSXqvj8h6SbJI2RNAa4HrhN0gpA1ah6FeoxDDgaGG1mmwL9CO/gOo7TRbiWHacYqt3Rbx8/V6zT55eA/YD3AwIuAa4xMwN2aqKeA2KPwPLA7Ab9OI5TLK5lx2kz1R7Gqxrkxsyeq5BuwNVxaRoze1LST4DHgAXALWZ2Sx6+HcdpH65lxymGag/jTSaMo6lM3jsC4Ei6PcbFnx/z384itP+DGqmgpJWBvYB1CV3/V0n6lJldlrA5AjiiEf+O47SHLFqOdq5nx8mRal33dQW+aaKrvxa7Ag+b2TMAksYD2wFv/zl4JC3H6Qpqahlcz46TN9W67jcyswckjSqXb2ZTKuy3HvCEmb0uaUdgM+CSWtPbVuEx4H0x0MYCwqxXHgPbcboP17LjFEC1rvtjCd1nZ5bJM8KkNuW4Bhgt6d3AhcANwBXAHo1U0MzuknQ1MAVYBNxLvNp3HKd7cC07TjFU67ovjZF9OP2eq6Tlqvh8y8wWSdoH+KmZ/UzSvc1U0sxOBk5uxkcv45Hx2mNXZNk9FBnPtZwTeeretdI+uyLKrjlNLSFCXpa0Egtj+NwxwE0xrX+GchzHcRzHyZlqY/RDgWGEd163YPHT94MI779W4jDgC8APzOxhSeuSetimXiQdA3wu1uF8M0tPstOn8ch47bUrsuweiIznWs4Jj4zX3Xat8lmOamP0HyLMUrcWcFYi/SXgxEo7mdn9ko4D1onbDwOnN1pBSZsS/hi2Bt4Abpb0BzN7sFGfjuO0H9ey4xRDxa57MxtnZjsBh5rZTollLzMbX2k/SXsCU4Gb4/bmqdnv6mVj4E4ze9XMFgF/B/Zpwp/jOMXgWnacAqjWdX9sat2AZ4Hb4116JU4hXLFPBDCzqbH7vlFmAj+QNITwSs4e+Cs5jtONuJYdpwCqdd2XC3wzAjhJ0ilm9tsK+y0ysxelJQLqNRz0wsxmSToDmAC8DEwjvJrzNh5Jy3E6nyxaBtez4+RNtdfrTi2XHmPg/wWo1NDPlHQw0E/S+oTZqqo9pV8TM7uQ8E4+kn4IPJHK90hajtMF1NJytHE9O06OZHm9bgniZDbl4t+XOArYBHgduJLw8N5XGqpdRNJq8XMdYN/o13GcLsO17Djtp1rXfVkk7Qw8XynfzF4FTopLXlwTx/UWAl8ys4rlO47T0biWHafNVHsYbwbvHFtfhTB/9CEV9hkDHANsGJNmAeea2SXNVNLMPtDM/r2OR8Zrj12RZfdQZDzXck54ZLzutCui7Gp39B9NbRswz8xeKWcs6RBCF/2xhFjWAkYBP5ZEs4294ziO4zj1U+1hvEfr9HUksI+ZPZJI+5uk/QgP7nlD3yI8Ml577Yosu9sj4zn54ZHxutuuVT7LUffDeFUYlGrkAYhpg2rtLOk3kuZKmplIW0XSBEkPxs+Vc6yv4zgtwLXsOJ1Fng39ggbzSlwM7J5KOx74q5mtD/w1bjuO09lcjGvZcTqGup+6r8LGkqaXSRfwP7V2NrPbJI1IJe8F7BjXxxGi7R3XcA0dx2k5rmXH6Sxybehz9FVidTN7CsDMniq9g5vGI2k5TseTScvgenacvMmzoX/MzKpGsZKkWjaN4JG0HKd3cD07Tr7kOUZ/q6SjYsSrt5G0jKSdJY0DxtTp82lJa0Q/awBzc6qr4zjtxbXsOAWRZ0O/O/AmcKWk2ZLul/QQ8CBwEHC2mV1cp88bWHxxMAa4Pq/KOo7TVlzLjlMQuXXdm9lrwHnAeZL6A6sCC8zshSz7S7qS8LDOqpKeAE4GTgd+L+kzwGPAAXnVt5fwyHjtsSuy7G6KjOdabg8eGa877YooWy0YMi8USc8A6WA/qwLPZnSR1bZX7Iosu9Ptiiy7HXbDzexdGfYtjCb17OdN59kVWXan2zXrs7KezaznF2BS3ra9YtcNdfRj097v3OlLNxzDTq+jH5vOs2uVTzPLdYzecRzHcZwOwxt6x3Ecx+lh+kpDP7YFtr1iV2TZnW5XZNlFfudOpxuOYafX0Y9N59m1ymfvPYznOI7jOM5i+sodveM4juP0SbyhdxzHcZwexht6x3Ecx+lh+nRDL2kjSbtIGphKT8+ljaStJW0V198j6VhJe9Twf0nGemwf/e2WSt9G0qC4PkDSqZJulHSGpJUSdkdLWjtDOctIOkTSrnH7YEk/l/SlGM0wabuepK9LOkfSmZK+kCzTcTqNrHpuRMvRtqaem9VyzHc9O7nS5x7Gk3SYmV0k6WjgS8AsYHPgGDO7PtpMMbNRiX1OBj5MCBk8AdiGMJ/2rsCfzewHkm5IFwXsBPwNwMw+lvB3t5ltHdc/F+txLbAbcKOZnR7z7gNGmtkiSWOBV4GrgV1i+r7R7kXgFeC/wJXAVWb2TJnvfnn8DssDLwADgfHRn8xsTLQ7GtgT+DuwBzAVeB7YBzjSzCZmOtgdjKTVzCy3iVUkDTGzeXn5a4b4B34CsDdQipQ1lxBf/nTLGJa60ylpOa5n0nMWLUe7THrOW8vR1vVcB72sZchJz/VE1+mFhTCdLsAMYGBcHwFMIvw5ANyb2mcG0I8gqJeAQTF9ADA9rk8BLiPE+N4hfj4V13dI+bs3sX4P8K64vgIwI5E3K7E+JeVjatIfoXdmN+BC4BngZsLkISsm7Ep1XRp4GugXt1XKS37fuL48MDGur1Pm2KxEiGP+ADAvLrNi2uCMv8mfEuuDgNOAS4GDU3bnJdaHAr8EfgEMAU6J9f49sEZqv1VSyxDgEWBlYJWE3e6p73UhMB24gjCfeinvdGDVuD4aeAj4DyFU6w4JuynAt4D1anz/0cCt8fxZm9AAvRjPjS0SdgOB7wL3xfxngDuBQ1P+/gwcBwxNHa/jgAlFazBvLSfO2Zp6JoOWE79dTT2Ts5ZLPilAz+SgZWuxnulwLSf8dJSeCxdrK5b4g5ZbZgCvR5v7U/sMjGI6q5zwyq3H7anxcyngq/FH3TymPVShftPiiTmEVCjDVFlXAYfF9YuA0XF9A+Ce5EmY8tEf+BjhbuCZRPpMYJlY9vySMIDlWPKPaAawbFxfGZic9NHISQiMqrBsCTyVsLsmim9vwoxn1yTqMiVhdzNwFHB8/G2PI/xxHQVcn6rjW8DDqWVh/Hyo3HEELgC+DwyPv+t1yeOTWL8V2Crxu0xK5D0M/IQwicvd0c+aZc6Huwl3mQcBjwP7x/RdgH8l7K4HDgXWAo4Fvg2sD4wDfpiw+3cVbVTM68SFDFquR89k0HI9eiZnLRepZ+poUChIz3S4ljtVz4ULuRUL4ep28/jDJpcRwOxo8zeigBP7LQ1cAryZSr8LWD6uL5VIX4l3inItgqh/TuKOI2XzCOHK8eH4OTSmD2TJP5uVgIsJXXh3xRP6IUIX3MiE3b1VjsWAxPpX4/6PAkcDfwXOJ/wRnJywOyaKbSzh6r70B/Uu4LasJ1oyjzCF8d+imNLLgoRd+iLrJOAOwh/plHLfOX2cy/j4OuGP5L2JtIfL1HdKFR/J3+UBYOm4fmfKbkYFfx8gzO44J37nIzJ+l2TetFTePaVzEnggkX4L8E2WvHNZnfDn+Zd267GZhQxajnaZ9EwdWo7pVfVMzlpO/+ZlymuZnqmjQaEgPdPhWs7wXQrRc+FCbsVC6KbZvkLeFfFzLRJXrimb96e2l61gt2ryhEvlfYTEVVnGei8PrFsmfUVgJOFqefUy+RvUUcaaxCtRYDCwP7B1GbtNYt5GNfxlOgkJdx/rV/DxeGJ9Fok/4Jg2htC99WgibVpi/fsp+xllyij9YZ8Vj2e5u7MnCFfWXyP8gSqRl+wKPSp+750JXYw/BT4InApcmrAr13D0A3YHLkqk/YvQTXsA4U9775i+A0veVfyzdF4Txlz/nLbWkacAAAhdSURBVMhLXlStDJxB+BN7HnguHtczSHRvdsOSRcuJ37emnhvRcsyvS8+NajnaFKLnrFqO6YXpuZO1HNM7Ts+FC9mX7l5SJ+FzqZNw5YTd/sCGFXzsnVj/EbBrGZvdgQcT298ljsmm7N4NXF2lvnsSxsHmlMk7ObWUxluHApekbHcEfkcYT50B/BE4gnh3EG1+m/EYjiR0m/4J2Ag4h/Bw1X3Adgm7zQjdgi8AtxMbBMKd2dEpnxsRHjAbmErfPUudfOl7S1YtR9vC9dyJWo62Hafnwk8uX3p3IXYRdpod4cGrTTu5js3YEbpx/w1cR+ha3iuR9447E198qbVkPQ/rsc3Drpu03KjPPPTc516vc9qHpMfMbJ1OteuGOjZiJ2kGsK2ZvSxpBOE1rkvN7BxJ95rZFrX8OU6Svq6pPOwa9ZmHnpfOUjnHqYSk6ZWyCON7hdp1Qx1b8J37mdnLAGb2iKQdgaslDY+2jvMOXFPN27XIZ9N69obeaZbVgQ8RHhJJIsLDJkXbdUMd87abI2lzM5sKEO8EPgr8BngvjlMe11Tzdq3w2bSevaF3muUmwgMiU9MZkiZ2gF031DFvu0OARcl8M1sEHCLp1+l9HSfimmrerhU+m9azj9E7juM4Tg/Tpye1cRzHcZxexxt6x3Ecx+lhvKF3HMdxnB7GG3rHcRzH6WG8oa8TSYMlHZlKS79e4dRA0gBJf5fUL26/KWlqYjm+Qb8vx8+Jkj6UyvuKpPOar33ZckdImlkhr6HzI89zTdIykm6T5G/aRFzL+eBazuSzUC17Q18/g4ElfjAz266gurQMBVp5fhwOjDezN+P2AjPbPLGc3qT/K4FPpNI+EdPbShPnR27nmpm9QZjd7MAG69KLuJbzwbVcm2K1nDVGbzcvhPmfjwLencH2WsLcxf8gTEO4ayr/t8ACYCrw45j2cvwcQZgQ4gLC7E6XEyYiuAN4kMSsUsCnCBMaTAV+TZgJaQXgD4Q5rmcCB5ap3yGEKSenseTsSsfGfWYCX4lpZwBHJmxOAb5WpfwRhEksziNM8DCcEF95MmFChuTUqt+O33UCQXBfr/bdynyPfwIjEtsvl7GpVv9K9Sr9FkOAZ1g89/UIwlzSqrQ/sFU8tsvF3+I+YhztVL3KHevSbz8u+riaxdOhvlzr2JT7Xal+rpU9NtWOPWGyjT8WrUfXsmsZ13Ly+7Vcy4ULtx0LYQakYwjzGE8nzCb0IWC5MrYPlk50YF/eOQXhCGBmKi3557CIEK1oqXjy/YYQ6Wgv4LpotzFwI9A/bp8XT479gPMTfldKlbMJYXKDVeP2KvFzS8KsSysQ5sG+D9giLn9P7H8/sE6V8kcAbwHvS+xTKmNAFMMQYHQ8+QYQpolMHrOyvlPfYxlSM04R5reemlgOrFT/SvUqI8Q/ECeAAI4nCqzG/t8HfgL8AjihzPlR6ViPAIw4JWr83b+erFOV417pdx1B5XOt3LHZodqxJzQAzxStR9eyaxnXctJXy7XcJ8brzOw/wDmSfg7sQryCAj5OuOoHQNLywErA2TFpacIUgvXwsJnNiP7uA/5qZhYnJhgRbXYhnGT3SIJwgs4FrgB+IukM4CYz+0fK986EKRufjd/ruZi+PXCtmb0Syx0PfMDMzpW0mqQ1CVMfPm9mj0n6coXybyPMEX1nosyjJe0T19cG1gfeB1xvZgtieTcm7Ct9tySr8s7jusDMNk8fzHL1r1KveandS11+18fPw2t8r3mE6TLvAV4jzBqVpuyxBm4gzMN9R7S7LO7/k8S+lY7NSpT/XStiZvemjw2hUap47M3sTUlvSFrRzObXKqMTcS27lnEt163lPtHQSzoI2Jtw8O4CTgNuMbP0CboJMNkWjzVtRrhCrIfXE+tvJbbfYvHxFjDOzE4oU9ctgT2A0yTdYmbfTWYTrjTfsVuV+lxNmDt6KKH7qGL5CjMjvZLY3pHQXbmtmb0awzIuV6O8it8twYLoJwvvqH+VeqW5DjhL0ihggJlNybD/KoSr+/4x7ZWUz2rfPf3bpLcrHfejy9hmIX1sshz7ZQl/fF2Ja9m1jGu5RHYtZ7317+aFcMU/ijimU8XuMOC0xPYNwFYpmyGEK+VK3X0zE+kXA/un84D3ELrIVovbqxDG0NYkdkES/syuS5WzCfB/LO6aKnULjSJ0Yy5P6IaaCWyR2Oefcb81apSfrv9ewI1xfaN4Uu1IGP+aQhDPQEJX1der+S5zrB8n0d1KmXG9KvUvW69yfoDfE7oPT6n1vRK/+cHAScDPy9Sn7LFmcXffttHufBaPQb5c7dhU+V0rnmvljk2tYx/9zSpaj65l1zKu5SW+X/rY1Dr21KnlwoXbjgX4C4sfuEgue6XsziLx0AzwEOHqMe3virh/uQd4av45xO0D40k7nTD+9z7CWOP0mH4PMLpM2WNi2dOAixPp73ioJJE3A7g1lVau/HQdlwX+FG2uAiYmRHQK4U/hFsKDSp+r5rvM97iQxMNRvHNc7/RK9a9Rr/Sfwz4E0W5Ua3/CGNv4aNOPcMe4c5m6V3qA537gV9HvNSx+gGd+rWNT5Xcte65VOTYVjz3hjuHMovXoWnYtu5bbq2Wf1MZpCEkDLUyXuDxhPPAIi91pGfffAjjWzD7dskp2AJKGAFPMbHgH1GU84aGkfxddF6dzcC1no5u13CfG6J2WMFbSewhdfuPq+WOAtx9AuVVSP1s8jtpTxIdrJrLkQzyFIGkZQvexN/JOGtdyDbpdy35H7ziO4zg9jEfGcxzHcZwexht6x3Ecx+lhvKF3HMdxnB7GG3rHcRzH6WG8oXccx3GcHsYbesdxHMfpYbyhdxzHcZwexht6x3Ecx+lh/h/haxDtso6lGwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAACpCAYAAADDeP+kAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2dedgcVZWH3x8hQCCEJQhhD4PsSCAJIMjIOogosougEsAZHHEAZRzZZgTcAEdAGEVlEcLqCITVEQlKRFSWfCGQQHAY2Q1hCVsCYQmc+ePeNpWiu7q6u7qru7/zPk89XXXr1Lm3qutXt+rWrXNlZjiO4ziO058sUXYBHMdxHMdpH17RO47jOE4f4xW94ziO4/QxXtE7juM4Th/jFb3jOI7j9DFe0TuO4zhOH+MVveM4ziBA0kmSLirIl0n6YBG+ikDSYZLuKrsc3YpX9E4hSHpC0m5ZgpM0RdKbkuZLelHSJEmrd7qsjlM2US9vS1ollT49VqKji87TzL5rZv9YtF8ASeMl3SLpZUmvSHpY0nckrdSO/JzG8Ire6TT/YmbDgQ2BFYFzSi6P45TF48DBlQVJHwKGNeNI0pJFFaqJvLcHpgB/ADY2sxWBPYCFwJga25RW3sGIV/ROKZjZS8B1wOZll8VxSuJy4NDE8gTgssqCpE9Iul/Sa5KelnRqYt3o+OT/BUlPAb+N6YdKelLSXEn/UWlpi+tOlXRFavsJkp6KLWwnJ/xvI+lP8en8WUk/lLRUjf34HnCJmZ1uZs8BmNlTZnaKmU2J/g6T9AdJ50h6CTg1ph8haVZsCfi1pHUTZdhY0mRJL0n6s6RPJ9aNlHRTPDb3Ausn1v1I0lnJAkq6WdJX6v4jfYpX9E4pxCbL/YH7yy6L45TE3cAISZtIGgIcBFyRWP864UZgReATwJck7ZPysSOwCfAxSZsC5wOfBVYHVgDWrFOGHYCNgF2Bb0jaJKa/C3wVWAXYLq4/Kr2xpOXi+uty7O+2wGPAqsB34r6cBOwHfAD4PXB1wu9k4KpofzBwvqTNoq8fAW/G/TwiThUmAgdLWiL6WiWW/+ocZexLvKJ3Os15kl4BHgCeBY4ruTyOUyaVp/p/AB4B/lpZYWZTzGyGmb1nZg8SKqodU9ufamavm9kC4ADgZjO7y8zeBr4B1BvM5DQzW2BmDxA0OSbmPWBmd5vZQjN7AvhplbwBViLUI3MqCZK+F1sCXpf07wnb2Wb2X9HnAuCLwOlmNsvMFgLfBbaMT/WfBJ4ws0ui/TTCzcQB8aZof+Abcd9nEir3ynG7F3iVULkDfAaYUmltGIx4Re90mmPMbEUzW9PMPmtmL5RdIMcpkcuBQ4DDSDTbA0jaVtIdkl6Q9Crwz4Qn7CRPJ+bXSC6b2RvA3Dr5z0nMvwEMj3lvGDvXzZH0GqESTucN8DLwHuHJupLv1+N7+uuB5Lv4p1PbrgucG28KXgFeAkRohVgX2LayLq7/LDCK8PS/ZMrfkynfE4HPxfnPEY7zoMUresdxnJIwsycJnfL2BCalVl8F3ASsbWYrAD8hVISLuUjMPwusVVmQNAwY2WTRfkxoYdjAzEYQmtjTeWNmrwP3EJrf65FuXXga+GK88a9Mw8zsj3Hd71LrhpvZl4AXCB391k74Wifl+wpgb0ljCK82bshRvr7FK3qnHUjSMsmp7AI5ThfzBWCXWGkmWR54yczelLQN4ck/i2uBvSRtHzvOnUaVyjknywOvAfMlbQx8KcP268ARkk6QtCqApLWA9erk8RPgxMp7d0krSDowrrsF2FDS5yUNjdPWkjYxs3cJN0WnSlo29k2YkHRsZs8A9xGe5K+LrwoGLV7RO+1ge2BBcvLPaRynOmb2FzObWmXVUcA3Jc0jvG//RR0/DwFHAz8nPN3PA54H3mqiWF8j3FjMAy4E/jsj37uAXYCPAv8bm9lvJXxy918Z210PnAn8PL4emAl8PK6bB+xOeL8+m/CK4Uxg6bj5vxBeM8wBLgUuqZLFROBDDPJmewCZ1eur4TiO4/QakoYDrxCa3x8vuzydRtJHCU34o83svbLLUyb+RO84jtMnSNorNmcvB3wfmAE8UW6pOo+kocCxwEWDvZIHr+gdx3H6ib0JTd2zgQ2Az9gga7aNsQBeIXwJ8IOSi9MVeNO94ziO4/Qx/kTvOI7jOH2MV/SO4ziO08f03SdPkvxdhOPk40Uz+0DZhcjC9ew4uamp565/opf0M0nPS5pZdlkcp89Ihw1tO65nx2kbNfXc9RU9IRjCHmUXwnGcQrgU17PjdJTCm+7j0IBjCAMsLAAeamXUIDO7U9LoYkrnOE6ZuJ4dp/MUVtFLWh84HtgNeJQw8MAyhHjFbxCGOZzowQscx3Ecp3MU+UT/beB8wmhEi3WgiQMdHAJ8nsS4wUUh6UjgyKL9Oo7TeVzPjlMshQbMic32H47DDBbpdzRwi5ltnsPWe+k6Tj4GzGx8pzN1PTtOW6ip50I748Vm+bOK9Ok4juM4TvO0o9f9bZL2l9TsOMiLIelq4E/ARpKekfSFIvw6jtN5XM+O03kKj3Ufx05eDniX0OtegJnZiEIzqsGmm25qt99+e127OXPmADBq1KhMuzXXXLOQciUZGBiom++cOXMYN25cx/Ot5A31j81gs+uFMjZiN27cuFKa7hshj54b3OfCyuYUR95rInS3pvLYtSvvLD0X/nmdmS1ftE/HcRzHcZqjHd/RC/gssJ6ZfUvS2sDqZnZvk/7WBi4DRgHvAReY2bm17IcOHcoaa6yR238jtkUxatSonsg3r+1gsysz73bsS6doVMvQmJ67cZ+dfDRybeoXu3b5rEbNil7S2KwNzWxajVXnE0S8C/AtYD7wI2DrJsu4EPhXM5smaXlgQNJkM3u4SX+O45SDa9lxSiDrib7Se34ZYDzwAOF9+xbAPcAONbbb1szGSrofwMxelrRUswU0s2eBZ+P8PEmzgDUBvzg4Tg/hWnaccqjZ697MdjaznQmB8sea2XgzGwdsBfxfhs93JA0BDEDSBwhP+C0Tv7/dinCj4ThOj+JadpzOkefzuo3NbEZlwcxmAltm2J8HXA+sKuk7wF3A6S2VEpA0HLgO+IqZvZZad6SkqZKmzp07t9WsHMdpI1lajutdz45TIHk6482SdBFwBeEp/XPArFrGZnalpAFgV0JT/z5mVtM+D5KGEi4MV5rZpCp5XgBcADBmzBiPpOU4XUo9LYPr2XGKJk9FfzjwJeDYuHwn8ONaxpIuN7PPA49USWuY2Iv/YmCWmZ3djA/HccrHtew45VC3ojezN4Fz4pSHzZIL8X19K1EqPkIYDGeGpOkx7SQz+58WfDqO03lcy45TAlmf180gdqirhpltkbI/ETgJGCbpNUKzPcDbxGa4ZjCzuxK+6vLOO+8we/bsunaViENFMTAwAOSLYJQn70bKlzeqVF6fbtd9eZe5z0XRqJYhn54b2eeyolI62RR5Tex2uzLyznqi/2SuHCJmdjpwuqTTzezERrZ1HMdxHKc91KzozezJyryk1VgU8OZeM3s+w+fJkj5HQZHxGqWsyHiVp4Qyoni1I6pUI7aDza7MvAdblLh2RMbrl2PTT3hkvOJ8VqPu53WSPg3cCxwIfBq4R9IBGZv8CNgOOCQuVyLjNYWkZSTdK+kBSQ9JOq1ZX47jlIdr2XHKIU+v+5OBrStP8TEAzu3AtTXsC42MB7wF7GJm8+OnOXdJ+pWZ3d2CT8dxOo9r2XFKIE9Fv0SqqX4u2S0BhUbGszCO7vy4ODRO/m2t4/QYrmXHKYc8kfFulfRrSYdJOgz4JZD1OUy1yHjfbaWQkobEz3GeByab2T2p9R5Jy3F6gHpajjauZ8cpkLoVvZn9G+HzuC2AMYShJY/PsL8S+Doh7O2zhMh417RSSDN718y2BNYCtpG0eWr9BTEW//iRI0e2kpXjOG2knpajjevZcQok13j0ZnYdIWxlXp4Dfh/9D5M0NmNY29yY2SuSpgB7ADNb9ec4Tjm4lh2nc9St6CXtB5wJrEoIdiHC67YRNey/BRwG/IVF79+MMD59w8R3/O/EC8MwYLdYHsdxegjXsuOUQ54n+u8BezUwMM2ngfXN7O3mi7UYqwMTYwe/JYBfmNkttYzLiozX7VGWujlqU6/YlZl3P0TGo0EtQ/GR8Yq0c4rDI+O11y5PRf9cg6PPzQRWJHS2aRkze5AwbnUuZs2axZprrllE1qWTN1yn4/QCjWq5TBrRXtHhcosO01uGv0Z9Ou0lK9b9fnF2qqT/Bm4gfAcLQK0hJgmd8O6XNDNl/6nWi9u9tCsyXj9FbeoXuzLzHmzR38qMjFfWMWwkSlw/+GvUZ7/YtctnNbKe6PdKzL8B7J5YNqBWRT+R8N5tBi18P58mNvdNBf5qZg3F4Xccp3twLTtOZ8mKdX94kz5fNLPzmtw2i2OBWUDVToCO4/QMrmXH6SB5Yt1/T9IISUMl/UbSi3HQmloMSDpd0naSxlamVgopaS3gE8BFrfhxHKdcXMuO03nydMbb3cy+Lmlf4BnC4DZ3AFfUsK90tvlwIq3pz+siPyAE4Vm+2kpJRwJHtuDfcZzOkKllWFzP/dKx1nHKJE9FPzT+7glcbWYvSappbGY7F1GwCpI+CTxvZgOSdqqR5wWE6H1I8tjZjtOF5NEyLK7nMWPGuJ4dp0XyVPQ3S3oEWAAcFYNevJm1gaRPAJsBy1TSzOybTZbxI8CnJO0Z/Y2QdIWZZb0+cByn+3AtO04J5Il1fwJhfPnxZvYOoQf+3rXsJf0EOAg4mhBF70Bg3WYLaGYnmtlaZjYa+AzwW78wOE7v4Vp2nHLIEwJ3WeDLwDqE92ZrABsBtSJabW9mW0h60MxOk3QWtT/F6xu6PRpTN0dt6hW7MvMerNHfujkyXruOYRmRLsvUQF7bbv+fu/k6kqfp/hJgANg+Lj8DXEPtir7SrP+GpDUI49evlyOfupjZFGBKEb4cxymPvFouK9LlwMBAx/NslKIjZ7YjEqdH9+wO8lT065vZQZIOBjCzBcrqjRfe6a8I/CcwjdDj/sLWi9rdeGS8wWNXZt7dHtWtX2hH9Ld25N3tdmXm3e127fJZjTwV/dtxpCkDkLQ+idC2SSQtAfzGzF4BrpN0C7CMmb3aSiElPQHMA94FFprZ+Fb8OY5TDq5lx+k8eSr6U4BbgbUlXUnoOXtYNUMzey++k98uLr9FjZuCJtjZzF4syJfjOOXhWnacDpJZ0ccm+keA/QgBcAQcW0ekt0naH5hkZv4NrOM4juOUSGZFb2Ym6QYzGwf8MqfP44DlgIWS3iTcHJiZtRLX2gg3EAb8NAbU+BseGc9xeoZMLYPr2XGKJk/T/d2Stjaz+/I4NLOaoS1b4CNmNlvSqsBkSY+Y2Z2JPD0ynuP0BplaBtez4xRNnop+Z+CLkp4EXmfRE/oWtTaQtBKwAYtHxruzln09zGx2/H1e0vXANkDT/hzHKQfXsuN0njwV/ccbcSjpHwnDUK4FTCe82/8TTQ5qI2k5YAkzmxfndweaDafrOE5JuJYdpxzqVvRm9qSkIcBqeewJlfzWwN1mtrOkjYHTWijjasD18dP9JYGrzOzWFvy1hW6PxtTNUZt6xa7MvHslqlsdekLLUHw0ubLydq10n10ZeecJgXs04RO754D3YrIBtZru3zSzNyUhaWkze0TSRrlKWwUzewwY0+z2juN0B81oeZNNNuH222/PtKlc5PJGYMtjN27cuAZKmQ+PEueURd4n9I3MbG5On8/EyHg3EDrbvAxkB6vuAzwy3uCxKzPvwRYZb+jQoX1zbPop4t1g1EovHJta5KnonwZyR7Yzs33j7KmS7gBWIATcaZp443ARsDmhNeEIM/tTKz4dx+k8rmXH6Tw1K3pJx8XZx4Apkn5JIsqdmZ2dse1YYAeCkP9gZm+3WM5zgVvN7ABJSwHLtujPcZxycC07TofJGo9++Tg9BUwGlkqkDa+1kaRvABOBkcAqwCWS/r3ZAkoaAXwUuBjAzN6OsfQdx+khXMuOUw41n+jN7DQASQea2TXJdZIOzPB5MLCVmb0Zbc8gjGL37SbL+HfAC4QbhjGEIXOPNbPXE+XxSFqO0/3U1TIsrucyhqh1nH4j64m+wok50yo8QSJQDrA08JcGypRmSWAs8GMz24oQtOeEpIGZXWBm430kLMfpaupqGRbX88iRIztdRsfpO7Le0X8c2BNYU9J5iVUjgIUZPt8CHpI0mfCO/h+Auyo+zOyYBsv4DPCMmd0Tl6+lysXBcZyux7XsOCWQ1et+NjAV+BShia3CPOCrGdtdH6cKU5otHICZzZH0tKSNzOzPwK7Aw634dByn87iWHaccst7RPwA8IOkqM3snr0Mzm1hIyRbnaODK2Ev3MeDwNuTREt0ejambozb1il2ZefdJZDxoUMvvvPMOs2dnh+HolWPTDxHvXCut25WRd54QuLkr+XZhZtMBf//uOD1Oo1qeNWuWd8hznBbJEzDHyYFHxhs8dmXm3e3R35zaeGQ8t+uEz2rU7XUvafOWcmgRSRtJmp6YXpP0lTLL5DhO47iWHacc8jzR/yS+T7uUMNpUZoALSRsC/wasm/RvZk0NUxs77WwZfQ8B/srinf0cx+kBXMuOUw553tHvIGkD4AhgqqR7gUvMbHKNTa4BfgJcCLxbWEkDuwJ/MbMnC/brOE5ncS07TofI9Y7ezB6NYWynAucBWykMKn2SmU1KmS80sx8XXM4KnwGubpNvx3E6h2vZcTpEnnf0W0g6B5gF7ALsZWabxPlzqmxys6SjJK0uaeXK1GpB4+uDTxFaDNLrjpQ0VdLUVvNxHKe9ZGk5rnc9O06B5Hmi/yGhGf4kM1tQSTSz2TUGq5kQf/8tkWaEONet8HFgmpk9l15hZhcAFwBIshbzcRynvdTUMrieHado8lT0k8zs8mSCpGPN7Nx0OoCZrVdY6RbnYLypz3H6Adey43SQPBX9ocAPUmmHEcaVfh+S9quS/Coww8yeb6h0i3wuS4iZ/8VmtnccpztwLTtO58ka1OZg4BBgPUk3JVYtD8zN8PkFYDvgjri8E3A3sKGkb1ZrBaiHmb1BGN++a2kketfAwEBdm8EUnrFX7MrMu1fCvNajF7TcLjwErtuVlXfWE/0fgWeBVYCzEunzgAcztnsP2KTy/k3SasCPgW2BO4GGK3rHcRzHcZoja1CbJ4EnCU/njTA61cnmeWBDM3tJUlNx8yV9FfhHQqe+GcDhZvZmM766gV4IhVlm3t1uV2bevR4Ct9+03Ai9oHvXSufs2uWzGjU/r5N0V/ydF0NVVqZ5kl7L8Pl7SbdImiBpAnAjcKek5YDMqHo1yrEmcAww3sw2B4YQvsF1HKeHcC07TjlkPdHvEH+Xb9Dnl4H9gY8AAi4DrjMzA3ZuoZzDYovAskD2uJWO43QrrmXH6TBZnfEyg9yY2Us10g24Nk4tY2Z/lfR94ClgAXCbmd1WhG/HcTqHa9lxyiGrM94A4T2aqqx7XwAcSXfFuPjz4vq/rSLU/yOaKaCklYC9gfUITf/XSPqcmV2RsDkSOLIZ/47jdIY8Wo52rmfHKZCspvuGAt+00NRfj92Ax83sBQBJk4Dtgb9dHDySluP0BHW1DK5nxymarKb7jc3sEUljq603s2k1tlsfeMbM3pK0E7AFcFm94W0zeAr4cAy0sYAw6pXHwHac3sO17DglkNV0fxyh+eysKuuMMKhNNa4Dxkv6IHAxcBNwFbBnMwU0s3skXQtMAxYC9xPv9h3H6R1cy45TDllN95V3ZB9Pf+cqaZkMn++Z2UJJ+wI/MLP/knR/K4U0s1OAU1rx0U10c4SsMvPudrsy8+6jyHh9peVG6Gbdu1Y6Z1dG3nWHqSVEyMuTVuGdGD53AnBLTBuaIx/HcRzHcQom6x39KGBNwjevW7Go9/0IwvevtTgc+GfgO2b2uKT1SHW2aRRJxwL/FMtwoZmlB9npKXohQlaZeXe7XZl590FkvL7SciP0gu5dK52za5fPamS9o/8YYZS6tYCzE+mvASfV2sjMHpZ0PLBOXH4cOKPZAkranHBh2AZ4G7hV0i/N7NFmfTqO03lcy45TDjWb7s1sopntDBxmZjsnpr3NbFKt7STtBUwHbo3LW6ZGv2uUTYC7zewNM1sI/A7YtwV/juOUg2vZcUogq+n+uNS8AS8Cd8Wn9FqcSrhjnwJgZtNj832zzAS+I2kk4ZOcPfFPchynF3EtO04JZDXdVwt8Mxo4WdKpZvbzGtstNLNXpcUC6jUd9MLMZkk6E5gMzAceIHya8zc8kpbjdD95tAyuZ8cpmqzP606rlh5j4N8O1KroZ0o6BBgiaQPCaFVZvfTrYmYXE77JR9J3gWdS6z2SluP0APW0HG1cz45TIHk+r1uMOJhNtfj3FY4GNgPeAq4mdN77SlOli0haNf6uA+wX/TqO02O4lh2n82Q13VdF0i7Ay7XWm9kbwMlxKorr4nu9d4Avm1nN/B3H6Wpcy47TYbI6483g/e/WVyaMH31ojW0mAMcCG8WkWcB5ZnZZK4U0s79vZftuo5sjZJWZd7fblZl3H0XG6ystN0I369610jm7MvLOeqL/ZGrZgLlm9no1Y0mHEprojyPEshYwFvhPSbRa2TuO49RiYGCAUaNGZdrMmTOHcePGdahEjtM9ZHXGe7JBX0cB+5rZE4m030ran9Bxzyv6SC9EyCoz7263KzPvXo+M1y4a0VRZ9ILuXSuds2uXz2o03BkvgxGpSh6AmDai3saSfibpeUkzE2krS5os6dH4u1KB5XUcpw24lh2nuyiyol/Q5LoKlwJ7pNJOAH5jZhsAv4nLjuN0N5fiWnacrqHhXvcZbCLpwSrpAv6u3sZmdqek0ankvYGd4vxEQrS945suoeM4bce17DjdRaEVfYG+KqxmZs8CmNmzlW9w03gkLcfpenJpGVzPjlM0RVb0T5lZZhQrSapn0wweSctx+gfXs+MUS5Hv6O+QdHSMePU3JC0laRdJE4EJDfp8TtLq0c/qwPMFldVxnM7iWnackiiyot8DeBe4WtJsSQ9Legx4FDgYOMfMLm3Q500sujmYANxYVGEdx+kormXHKYnCmu7N7E3gfOB8SUOBVYAFZvZKnu0lXU3orLOKpGeAU4AzgF9I+gLwFHBgUeUtk26OkFVm3t1uV2bevRQZrwwtF6mpdtHNunetdM6ujLzVhlfmpSLpBSAd7GcV4MWcLvLa9otdmXl3u12ZeXfCbl0z+0CObUujRT37edN9dmXm3e12rfqsrWcz6/sJmFq0bb/Y9UIZ/dh0dp+7feqFY9jtZfRj03127fJpZoW+o3ccx3Ecp8vwit5xHMdx+pjBUtFf0AbbfrErM+9utysz7zL3udvphWPY7WX0Y9N9du3y2X+d8RzHcRzHWcRgeaJ3HMdxnEGJV/SO4ziO08d4Re84juM4fcygruglbSxpV0nDU+npsbSRtI2kreP8ppKOk7RnHf+X5SzHDtHf7qn0bSWNiPPDJJ0m6WZJZ0paIWF3jKS1c+SzlKRDJe0Wlw+R9ENJX47RDJO260v6mqRzJZ0l6Z+TeTpOt5FXz81oOdrW1XOrWo7rXc9OoQy6zniSDjezSyQdA3wZmAVsCRxrZjdGm2lmNjaxzSnAxwkhgycD2xLG094N+LWZfUfSTemsgJ2B3wKY2acS/u41s23i/D/FclwP7A7cbGZnxHUPAWPMbKGkC4A3gGuBXWP6ftHuVeB14C/A1cA1ZvZClX2/Mu7DssArwHBgUvQnM5sQ7Y4B9gJ+B+wJTAdeBvYFjjKzKbkOdhcjaVUzK2xgFUkjzWxuUf5aIV7ATwT2ASqRsp4nxJc/w3KGpe52KlqO87n0nEfL0S6XnovWcrR1PTdAP2sZCtJzI9F1+mEiDKcLMAMYHudHA1MJFweA+1PbzACGEAT1GjAipg8DHozz04ArCDG+d4y/z8b5HVP+7k/M3wd8IM4vB8xIrJuVmJ+W8jE96Y/QOrM7cDHwAnArYfCQ5RN2lbIuCTwHDInLqqxL7m+cXxaYEufXqXJsViDEMX8EmBunWTFtxZz/ya8S8yOA04HLgUNSducn5kcBPwZ+BIwETo3l/gWwemq7lVPTSOAJYCVg5YTdHqn9uhh4ELiKMJ56Zd0ZwCpxfjzwGPB/hFCtOybspgH/DqxfZ//HA3fE82dtQgX0ajw3tkrYDQe+CTwU178A3A0clvL3a+B4YFTqeB0PTC5bg0VrOXHO1tUzObSc+O/q6pmCtVzxSQl6pgAtW5v1TJdrOeGnq/RculjbMcU/tNo0A3gr2jyc2mZ4FNPZ1YRXbT4uT4+/SwBfjX/qljHtsRrleyCemCNJhTJM5XUNcHicvwQYH+c3BO5LnoQpH0OBTxGeBl5IpM8Elop5z6sIA1iGxS9EM4Cl4/xKwEDSRzMnITC2xjQOeDZhd10U3z6EEc+uS5RlWsLuVuBo4IT43x5PuHAdDdyYKuN7wOOp6Z34+1i14whcBHwbWDf+rzckj09i/g5g68T/MjWx7nHg+4RBXO6Nftaocj7cS3jKPBh4Gjggpu8K/ClhdyNwGLAWcBzwH8AGwETguwm7P2doo+a6bpzIoeVG9EwOLTeiZwrWcpl6poEKhZL0TJdruVv1XLqQ2zER7m63jH9schoNzI42vyUKOLHdksBlwLup9HuAZeP8Eon0FXi/KNciiPqHJJ44UjZPEO4cH4+/o2L6cBa/2KwAXEpowrsnntCPEZrgxiTs7s84FsMS81+N2z8JHAP8BriQcCE4JWF3bBTbBYS7+8oF6gPAnXlPtOQ6whDGv41iSk8LEnbpm6yTgT8QLqTTqu1z+jhX8fE1woXkQ4m0x6uUd1qGj+T/8giwZJy/O2U3o4a/vyeM7jgn7vOROfclue6B1Lr7Kuck8Egi/Tbg6yz+5LIa4eJ5e6f12MpEDi1Hu1x6pgEtx/RMPVOwltP/eZX82qZnGqhQKEnPdLmWc+xLKXouXcjtmAjNNDvUWHdV/F2LxJ1ryuYjqeWla9itkjzhUus+QeKuLGe5lwXWq5K+PDCGcLe8WpX1GzaQxxrEO1FgReAAYJsqdpvFdRvX8ZfrJCQ8fWxQwwdn0W4AAAjASURBVMfTiflZJC7AMW0CoXnryUTaA4n5b6fsZ1TJo3LBPjsez2pPZ88Q7qz/lXABVWJdsin06LjfuxCaGH8AfBQ4Dbg8YVet4hgC7AFckkj7E6GZ9kDCRXufmL4jiz9V/LFyXhPeuf46sS55U7UScCbhIvYy8FI8rmeSaN7shSmPlhP/b109N6PluL4hPTer5WhTip7zajmml6bnbtZyTO86PZcuZJ96e0qdhC+lTsKVEnYHABvV8LFPYv57wG5VbPYAHk0sf5P4TjZl90Hg2ozy7kV4DzanyrpTUlPlfeso4LKU7U7AfxPep84A/gc4kvh0EG1+nvMYjiE0m/4K2Bg4l9C56iFg+4TdFoRmwVeAu4gVAuHJ7JiUz40JHcyGp9L3yFMmnwbflFfL0bZ0PXejlqNt1+m59JPLp/6diE2E3WZH6Hi1eTeXsRU7QjPun4EbCE3LeyfWve/JxCef6k15z8NGbIuw6yUtN+uzCD0Pus/rnM4h6SkzW6db7XqhjM3YSZoBbGdm8yWNJnzGdbmZnSvpfjPbqp4/x0ky2DVVhF2zPovQ85J5Cuc4tZD0YK1VhPd7pdr1QhnbsM9DzGw+gJk9IWkn4FpJ60Zbx3kfrqnW7drks2U9e0XvtMpqwMcInUSSiNDZpGy7Xihj0XZzJG1pZtMB4pPAJ4GfAR/Ccarjmmrdrh0+W9azV/ROq9xC6CAyPb1C0pQusOuFMhZtdyiwMLnezBYCh0r6aXpbx4m4plq3a4fPlvXs7+gdx3Ecp48Z1IPaOI7jOE6/4xW94ziO4/QxXtE7juM4Th/jFb3jOI7j9DFe0TeIpBUlHZVKS39e4dRB0jBJv5M0JC6/K2l6YjqhSb/z4+8USR9LrfuKpPNbL33VfEdLmlljXVPnR5HnmqSlJN0pyb+0ibiWi8G1nMtnqVr2ir5xVgQW+8PMbPuSytI2FGjn+XEEMMnM3o3LC8xsy8R0Rov+rwY+k0r7TEzvKC2cH4Wda2b2NmF0s4OaLEs/4louBtdyfcrVct4Yvb08EcZ/Phr4YA7b6wljF/+eMAzhbqn1PwcWANOB/4xp8+PvaMKAEBcRRne6kjAQwR+AR0mMKgV8jjCgwXTgp4SRkJYDfkkY43omcFCV8h1KGHLyARYfXem4uM1M4Csx7UzgqITNqcC/ZuQ/mjCIxfmEAR7WJcRXHiAMyJAcWvU/4r5OJgjua1n7VmU//giMTizPr2KTVf5a5ar8FyOBF1g09vVowljSqrU9sHU8tsvE/+IhYhztVLmqHevKfz8x+riWRcOhzq93bKr9r2Sfa1WPTdaxJwy28T9l69G17FrGtZzcv7ZruXThdmIijIB0LGEc4wcJowl9DFimiu2jlRMd2I/3D0E4GpiZSkteHBYSohUtEU++nxEiHe0N3BDtNgFuBobG5fPjybE/cGHC7wqpfDYjDG6wSlxeOf6OI4y6tBxhHOyHgK3i9LvE9g8D62TkPxp4D/hwYptKHsOiGEYC4+PJN4wwTGTymFX1ndqPpUiNOEUY33p6YjqoVvlrlauKEH9JHAACOIEosDrbfxv4PvAj4MQq50etYz0aMOKQqPF//1qyTBnHvdb/Opra51q1Y7Nj1rEnVAAvlK1H17JrGddy0lfbtTwo3teZ2f8B50r6IbAr8Q4K+DThrh8AScsCKwDnxKQlCUMINsLjZjYj+nsI+I2ZWRyYYHS02ZVwkt0nCcIJ+jxwFfB9SWcCt5jZ71O+dyEM2fhi3K+XYvoOwPVm9nrMdxLw92Z2nqRVJa1BGPrwZTN7StK/1Mj/TsIY0Xcn8jxG0r5xfm1gA+DDwI1mtiDmd3PCvta+JVmF9x/XBWa2ZfpgVit/RrnmpjavNPndGH+PqLNfcwnDZd4HvEkYNSpN1WMN3EQYh/sP0e6KuP33E9vWOjYrUP1/rYmZ3Z8+NoRKqeaxN7N3Jb0taXkzm1cvj27EtexaxrXcsJYHRUUv6WBgH8LBuwc4HbjNzNIn6GbAgC1617QF4Q6xEd5KzL+XWH6PRcdbwEQzO7FKWccBewKnS7rNzL6ZXE2403zfZhnluZYwdvQoQvNRzfwVRkZ6PbG8E6G5cjszeyOGZVymTn419y3BgugnD+8rf0a50twAnC1pLDDMzKbl2H5lwt390Jj2espn1r6n/5v0cq3jfkwV2zykj02eY7804cLXk7iWXcu4livk13LeR/9engh3/GOJ73Qy7A4HTk8s3wRsnbIZSbhTrtXcNzORfilwQHodsCmhiWzVuLwy4R3aGsQmSMLF7IZUPpsB/8uipqlKs9BYQjPmsoRmqJnAVolt/hi3W71O/uny7w3cHOc3jifVToT3X9MI4hlOaKr6WpbvKsf6aRLNrVR5r5dR/qrlquYH+AWh+fDUevuV+M8PAU4GflilPFWPNYua+7aLdhey6B3k/Kxjk/G/1jzXqh2besc++ptVth5dy65lXMuL7V/62NQ79jSo5dKF24kJuJ1FHS6S094pu7NJdJoBHiPcPab9XRW3r9aBp+7FIS4fFE/aBwnv/z5MeNf4YEy/DxhfJe8JMe8HgEsT6e/rVJJYNwO4I5VWLf90GZcGfhVtrgGmJER0KuGicBuho9I/Zfmush8Xk+gcxfvf651Rq/x1ypW+OOxLEO3G9bYnvGObFG2GEJ4Yd6lS9lodeB4GfhL9XseiDjzz6h2bjP+16rmWcWxqHnvCE8NZZevRtexadi13Vss+qI3TFJKGWxgucVnC+8AjLTan5dx+K+A4M/t82wrZBUgaCUwzs3W7oCyTCJ2S/lx2WZzuwbWcj17W8qB4R++0hQskbUpo8pvYyIUB/tYB5Q5JQ2zRe9S+InaumcLinXhKQdJShOZjr+SdNK7lOvS6lv2J3nEcx3H6GI+M5ziO4zh9jFf0juM4jtPHeEXvOI7jOH2MV/SO4ziO08d4Re84juM4fYxX9I7jOI7Tx3hF7ziO4zh9jFf0juM4jtPH/D96bblhWfzxKAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This block of code was used for generating Figure 4. For additional analysis, benchmark_designs_scores,\n",
    "# benchmark_trajectories_scores, mhc1_design_scores, and mhc2_design_scores contain all necessary data\n",
    "# from the experiment\n",
    "\n",
    "# benchmark_designs_scores is a dictionary that maps from a previous design (indicated by a string\n",
    "# denoting the source of the design and the MHC class it is for) to an nTimesCoverage array.\n",
    "# The ith entry (zero indexed) of the array contains the (i+1)-times coverage score.\n",
    "\n",
    "# benchmark_trajectories_scores is a dictionary of dictionaries. The first index indicates a design\n",
    "# strategy from Liu et al., while the second index indicates the size of the design. The value mapped\n",
    "# to by the keys is a an nTimesCoverage array.\n",
    "\n",
    "# mhc1_design_scores, and mhc2_design_scores are multidimensional lists. The entry in [i][j][k] denotes\n",
    "# the (k+1)-times coverage score of a design of size j+1 that was optimized via the diminishing returns framework\n",
    "# with the threshold function U(x) = min(i+1,x).\n",
    "\n",
    "def compare(optivax, diminishingReturns):\n",
    "    dr19 = []\n",
    "    for design in diminishingReturns:\n",
    "        if len(design) < 19:\n",
    "            dr19.append(design[-1])\n",
    "        else:\n",
    "            dr19.append(design[18].reshape(-1))\n",
    "            \n",
    "    _, axs = plt.subplots(1,2, figsize = (8,3))\n",
    "    for ax, optivax_design in zip(axs, optivax):\n",
    "        matr = np.array([(optivax_design < design) for design in dr19], dtype = np.float32)\n",
    "        matr = matr[:10, :]\n",
    "        ax.imshow(matr, cmap = \"gray\")\n",
    "        ax.set_yticks([])\n",
    "        ax.set_xticks(np.arange(0,20,1))\n",
    "        ax.set_xticklabels(np.arange(1,21,1), rotation = 90)\n",
    "        ax.set_xlabel(\"$\\geq n$ times coverage (EvalVax objective)\")\n",
    "        \n",
    "        ax.set_xticks(np.arange(-0.5, 20, 1), minor=True)\n",
    "        ax.set_yticks(np.arange(-0.5, 10, 1), minor=True)\n",
    "        ax.grid(which='minor', color='black', linestyle='-', linewidth=2, alpha = 0.1)\n",
    "        ax.tick_params(axis='x', which='minor', bottom=False)\n",
    "        ax.tick_params(axis='y', which='minor', left=False)\n",
    "        \n",
    "        ax.set_yticks(np.arange(0,10,1))\n",
    "        ax.set_yticklabels(np.arange(1,11,1), rotation = 0)\n",
    "    axs[0].set_ylabel(\"Utility threshold\\n(Design parameter)\")\n",
    "    axs[0].set_title(\"ILP\")\n",
    "    axs[1].set_title(\"MarginalGreedy\")\n",
    "\n",
    "optivax1 = benchmark_trajectories_scores[\"ILP_n3_MHC1\"][19], benchmark_trajectories_scores[\"MarginalGreedy_MHC1\"][19]\n",
    "optivax2 = benchmark_trajectories_scores[\"ILP_n3_MHC2\"][19], benchmark_trajectories_scores[\"MarginalGreedy_MHC2\"][19]\n",
    "compare(optivax1, mhc1_design_scores)\n",
    "compare(optivax2, mhc2_design_scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}