EvaluationOfDiminishingReturns.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This notebook contains code for evalating the diminishing returns design against\n",
    "# prior designs (Figure 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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",
    "\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": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"./data.pkl\", 'rb') as fin:\n",
    "    mhc1_data, mhc2_data = pickle.load(fin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "779it [00:00, 929.25it/s] \n",
      "1200it [00:02, 550.79it/s]\n",
      "440it [00:00, 2962.11it/s]\n",
      "128it [00:00, 366.72it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([1195, 298])\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "537it [00:00, 1004.62it/s]\n",
      "920it [00:01, 636.28it/s]\n",
      "502it [00:00, 2767.39it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([4082, 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",
    "    credence_data, _, _, population = data\n",
    "    \n",
    "    credence_distribution = makeDiploidDistribution(population, credence_data)\n",
    "    credence_data = torch.tensor(1-credence_data[0], dtype = torch.float32), credence_data[2],\\\n",
    "        credence_distribution[0], torch.tensor(credence_distribution[1], dtype = torch.float32)\n",
    "    \n",
    "    print (credence_data[0].size())\n",
    "    return credence_data\n",
    "\n",
    "inputs1_credences = reformatData(mhc1_data)\n",
    "inputs2_credences = reformatData(mhc2_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Evaluation algorithm</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Utility functions\n",
    "\n",
    "# f(x) = min(n, x)\n",
    "def getThresholdUtility(n):\n",
    "    return torch.tensor(np.arange(0,n+1,1),dtype = torch.float32)\n",
    "\n",
    "def getMarginalImprovement(utility):\n",
    "    return torch.cat( (utility[1:] - utility[:-1], torch.zeros(1)) )\n",
    "\n",
    "def getMultiUtility(utilities):\n",
    "    maxlength = max([z.size(dim = 0) for z in utilities])\n",
    "    allUtilities = []\n",
    "    for utility in utilities:\n",
    "        extendedUtil = torch.ones(maxlength, dtype = utility.dtype) * utility[-1]\n",
    "        extendedUtil[:utility.size(dim = 0)] = utility\n",
    "        allUtilities.append(extendedUtil)\n",
    "    return torch.stack(allUtilities)\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\n",
    "\n",
    "def data_to_cuda(data):\n",
    "    candidates, seqs, columnIndex, columnWeights = data\n",
    "    return candidates.cuda(), seqs, columnIndex, columnWeights.cuda()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#candidates: [candidate, 1 - pMHC hit probability]\n",
    "#columnIndex: [diplotype, allele in diplotype]\n",
    "#columnWeights: [diplotype]\n",
    "#distributions: [dummy, diplotype, distribution]\n",
    "#marginalImprovement: [improvement (shifting from i to i+1, so last entry should be 0)]\n",
    "\n",
    "def evaluateCandidates(candidates, columnIndex, columnWeights, distributions, marginalImprovement):\n",
    "    probabilityOfHit = 1 - torch.prod(candidates[:, columnIndex], dim = 2).unsqueeze(2)\n",
    "    shiftedMass = distributions * probabilityOfHit\n",
    "    improvement = torch.sum( shiftedMass * marginalImprovement, dim = 2)\n",
    "    weightedImprovement = torch.sum( improvement * columnWeights, dim = 1)\n",
    "    return weightedImprovement\n",
    "\n",
    "def updateDistribution(newRow, columnIndex, distributions):\n",
    "    probabilityOfMiss = torch.prod(newRow[columnIndex], dim = 1).reshape(1, -1, 1)\n",
    "    shiftedMass = distributions * (1-probabilityOfMiss)\n",
    "    \n",
    "    convolution = distributions * probabilityOfMiss\n",
    "    convolution[:,:,1:] += shiftedMass[:,:,:-1]\n",
    "    convolution[:,:,-1] += shiftedMass[:,:,-1]\n",
    "    return convolution\n",
    "\n",
    "def evaluateDesign(candidates, seqs, columnIndex, columnWeights, design, utility):\n",
    "    design = dedupe_designs(design)\n",
    "    distributions = torch.zeros( (1, len(columnIndex), len(utility)) )\n",
    "    distributions[:, :, 0] = 1\n",
    "    \n",
    "    seqToIndex = {}\n",
    "    for i, seq in enumerate(seqs):\n",
    "        seqToIndex[seq] = i\n",
    "        \n",
    "    for seq in design:\n",
    "        if seq not in seqToIndex: continue\n",
    "        row = seqToIndex[seq]\n",
    "        distributions = updateDistribution(candidates[row], columnIndex, distributions)\n",
    "        \n",
    "    scores = torch.sum( distributions * utility.view(1,1,-1), dim = 2 ).reshape(-1)\n",
    "    return torch.sum(scores * columnWeights).item()\n",
    "\n",
    "def evaluateDesignTrajectory(candidates, seqs, columnIndex, columnWeights, design, utility):\n",
    "    design = dedupe_designs(design)\n",
    "    distributions = torch.zeros( (1, len(columnIndex), len(utility)) )\n",
    "    distributions[:, :, 0] = 1\n",
    "    \n",
    "    seqToIndex = {}\n",
    "    for i, seq in enumerate(seqs):\n",
    "        seqToIndex[seq] = i\n",
    "        \n",
    "    scores = []\n",
    "    for seq in design:\n",
    "        row = seqToIndex[seq]\n",
    "        distributions = updateDistribution(candidates[row], columnIndex, distributions)\n",
    "        currentScores = torch.sum( distributions * utility.view(1,1,-1), dim = 2 ).reshape(-1)\n",
    "        scores.append(torch.sum(currentScores * columnWeights).item())\n",
    "        \n",
    "    return scores\n",
    "\n",
    "def evaluateDesignTrajectoryMultiUtility_cpu(candidates, seqs, columnIndex, columnWeights, design, utility):\n",
    "    design = dedupe_designs(design)\n",
    "    distributions = torch.zeros( (1, len(columnIndex), utility.size(dim = 1) ) )\n",
    "    distributions[:, :, 0] = 1\n",
    "    \n",
    "    seqToIndex = {}\n",
    "    for i, seq in enumerate(seqs):\n",
    "        seqToIndex[seq] = i\n",
    "        \n",
    "    scores = []\n",
    "    for seq in tqdm(design, position = 0, leave = True):\n",
    "        row = seqToIndex[seq]\n",
    "        distributions = updateDistribution(candidates[row], columnIndex, distributions)\n",
    "        currentScores = torch.sum(\n",
    "            distributions.unsqueeze(0) * utility.view(utility.size(dim=0),1,1,utility.size(dim=1)), dim = 3 ).squeeze(1)\n",
    "        scores.append(torch.sum(currentScores * columnWeights.view(1,-1), dim = 1).numpy())\n",
    "        \n",
    "    return scores\n",
    "\n",
    "def evaluateDesignTrajectoryMultiUtility(candidates, seqs, columnIndex, columnWeights, design, utility):\n",
    "    design = dedupe_designs(design)\n",
    "    distributions = torch.zeros( (1, len(columnIndex), utility.size(dim = 1) ) ).cuda()\n",
    "    distributions[:, :, 0] = 1\n",
    "    \n",
    "    seqToIndex = {}\n",
    "    for i, seq in enumerate(seqs):\n",
    "        seqToIndex[seq] = i\n",
    "        \n",
    "    scores = []\n",
    "    for seq in design:\n",
    "        row = seqToIndex[seq]\n",
    "        distributions = updateDistribution(candidates[row], columnIndex, distributions)\n",
    "        currentScores = torch.sum(\n",
    "            distributions.unsqueeze(0) * utility.view(utility.size(dim=0),1,1,utility.size(dim=1)), dim = 3 ).squeeze(1)\n",
    "        scores.append(torch.sum(currentScores * columnWeights.view(1,-1), dim = 1).cpu().numpy())\n",
    "        \n",
    "    scores = np.array(scores).T\n",
    "    return scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Analysis</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "credence1 = []\n",
    "binarized1 = []\n",
    "credence2 = []\n",
    "binarized2 = []\n",
    "for threshold in range(1, 21):\n",
    "    with open(\"./outputs/mhc1_credences_threshold_{}.pkl\".format(threshold), 'rb') as fin:\n",
    "        credence1.append(pickle.load(fin))\n",
    "    with open(\"./outputs/mhc2_credences_threshold_{}.pkl\".format(threshold), 'rb') as fin:\n",
    "        credence2.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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "benchmark_designs_scores = {}\n",
    "for design in tqdm(benchmark_designs, position = 0, leave = True):\n",
    "    if \"MHC1\" in design:\n",
    "        data = inputs1_credences\n",
    "    elif \"MHC2\" in design:\n",
    "        data = inputs2_credences\n",
    "    else:\n",
    "        raise Exception(\"Unknown MHC Class\")\n",
    "    scores = [evaluateDesign(*data,benchmark_designs[design],getThresholdUtility(evaluate_threshold)) for\n",
    "                 evaluate_threshold in range(1,21)]\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_credences\n",
    "    elif \"MHC2\" in design:\n",
    "        data = inputs2_credences\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 = [evaluateDesign(*data,trajectory[designSize],getThresholdUtility(evaluate_threshold)) for\n",
    "                     evaluate_threshold in range(1,21)]\n",
    "        trajectoryScores[designSize] = scores\n",
    "    benchmark_trajectories_scores[design] = trajectoryScores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getAllTrajectories(designs, data):\n",
    "    trajectories = [[[] for i in range(20)] for j in range(20)]\n",
    "    with tqdm(position = 0, leave = True) as pbar:\n",
    "        for design_threshold, design in enumerate(designs, 1):\n",
    "            for evaluate_threshold in range(1,21):\n",
    "                trajectories[design_threshold-1][evaluate_threshold-1] =\\\n",
    "                    evaluateDesignTrajectory(*data,design, getThresholdUtility(evaluate_threshold))\n",
    "                pbar.update(1)\n",
    "    return trajectories\n",
    "\n",
    "mhc1_design_scores = getAllTrajectories(credence1, inputs1_credences)\n",
    "mhc2_design_scores = getAllTrajectories(credence2, inputs2_credences)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getSelectablePeptides(data):\n",
    "    _, _, peptides = data[0]\n",
    "    filtered = sorted(list( set(data[2]).intersection(set(peptides)) ))\n",
    "    print (len(filtered))\n",
    "    return filtered\n",
    "\n",
    "def calculateRandomTrajectories(candidates, size, rng, data, utility, iterations):\n",
    "    data = data_to_cuda(data)\n",
    "    utility = utility.cuda()\n",
    "    scores = []\n",
    "    for _ in tqdm(range(iterations), position = 0, leave = True):\n",
    "        randomTrajectory = [candidates[i] for i in rng.permutation(len(candidates))[:size]]\n",
    "        scores.append( evaluateDesignTrajectoryMultiUtility(*data, randomTrajectory, utility) )\n",
    "    return scores\n",
    "\n",
    "mhc1_candidates = getSelectablePeptides(mhc1_data)\n",
    "mhc2_candidates = getSelectablePeptides(mhc2_data)\n",
    "\n",
    "rng_mhc1 = np.random.RandomState(1054969295)\n",
    "random_scores1 = calculateRandomTrajectories(mhc1_candidates, 151, rng_mhc1,\n",
    "                            inputs1_credences, getMultiUtility([getThresholdUtility(z) for z in range(1,21)]), 1000)\n",
    "\n",
    "rng_mhc2 = np.random.RandomState(1800457161)\n",
    "random_scores2 = calculateRandomTrajectories(mhc2_candidates, 151, rng_mhc2,\n",
    "                            inputs2_credences, getMultiUtility([getThresholdUtility(z) for z in range(1,21)]), 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAC1CAYAAABifukJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2dd5hU1f2H3zNlZwudpSNFRUBkWZAiASkBsYIdo8SoRFGJBo09REUjsYsY/GlINJYQsQs2LLBEUEEEV0CKFKkqIAsLbJl2z++PM+fOnWGWbbOV8z7PPOzce+fOmbvs/cy3CyklBoPBYDAkG1dNL8BgMBgM9RMjMAaDwWCoEozAGAwGg6FKMAJjMBgMhirBCIzBYDAYqgQjMAaDwWCoEjw1vYDaRGZmpuzUqVNNL8NgMBjqFMuXL/9FStkifrsRGAedOnXi66+/rullGAwGQ51CCLE10XbjIjMYDAZDlWAExmAwGAxVghEYg8FgMFQJRmAMBoPBUCUYgTEYDAZDlWAExmAwGAxVghEYg8FgMFQJRmDqAA0aNABgy5YtnHTSSYftv/LKK+ncuTPZ2dn06dOHL7/8EoC7776brKwssrOzGTVqFD/++GO1rrsqcbvdZGdn06NHD3r16sUTTzyBZVkALFy4kMaNG9O7d2+6devGrbfeWsOrNRiOTkyhZZKZNWsWkydPZtu2bXTo0IGpU6cybty4Kn/fRx99lIsuuoiPP/6Ya6+9lpUrV3Lbbbfx17/+FYCnnnqK+++/n2effTbp733PPbBtW/LO16ED3H//kY9JS0sjNzcXgN27d3PZZZeRn5/PfffdB8Cpp57Ke++9R1FREb179+b8889n0KBByVtkPWTu+rn4w/6knc/n9jGm65gjHuN2u+nZsyehUIjOnTvz8ssv06RJk0q/95YtWzjnnHNYvXp1pc8FcNNNN3HBBRcwZMgQhg0bxmOPPUbfvn3p1KkTDRs2BCAcDnPBBRdw99134/P52LJlC927d6dr164EAgH69u3Lc889h9frZdWqVTz++OO88MILSVlfbcVYMElk1qxZTJgwga1btyKlZOvWrUyYMIFZs2ZV2xqGDBnCxo0bAWjUqJG9vaCgACFElbzntm3QqVPyHuUVq5YtWzJz5kxmzJhB/ITWtLQ0srOz2blzZ0U/3lFDMsWlrOfTXxRWr15Ns2bNePrpp5O6hmSQl5fHkiVLGDJkSML9OTk5rFq1iq+++orNmzczYcIEe99xxx1Hbm4uq1atYseOHbz22msA9OzZkx07drAtmd/MaiFGYMrJsGHD7G8dwWCQYcOG8Z///AeAu+66i8LCwpjjCwsLuemmmwD45ZdfGDZsGO+++y4AP//8c9LX9+6779KzZ0/7+eTJkznmmGOYNWsW95dmFtRhjj32WCzLYvfu3THb9+3bx4YNG0q8ORhqDwMHDrS/CBw6dIgRI0bQp08fevbsyZw5cwBsq+Caa66hR48ejBo1iqKiIgCWL19Or169GDhwYIxQFRcXc9VVV9GzZ0969+5NTk4OAC+88ALnnXceo0ePpnPnzsyYMYMnnniC3r17c8opp5CXlwfAG2+8wRlnnFHq+hs0aMCzzz7LO++8Y79W43a76d+/f8wXndGjRzN79uxKXLHajxGYJLJjx46E2/fu3Vvl733bbbeRnZ3NzJkzee655+ztU6dOZfv27YwbN44ZM2ZU+TpqEqf1smjRIrKysmjdujXnnHMOrVu3rsGVGUojHA4zf/58xoxRLrXU1FTefvttVqxYQU5ODrfccov9+92wYQN/+MMf+O6772jSpAlvvvkmAFdddRVPPfWUHYPUaLFZtWoVr7zyCldccQXFxcUArF69mv/+97989dVXTJ48mfT0dL755hsGDhzISy+9BMDnn3/OySefXKbP0ahRIzp37syGDRtithcXF7N06dIYoerbty+LFi0q76WqUxiBKScLFy7kyiuvBMDr9bJw4UJ++9vfAtChQ4eEr9HbMzMzWbhwIaNHjwZI6k3v0UcfJTc3l08++SRhIsBll11m/yHWRzZv3ozb7aZly5aAisGsXLmSVatW8cwzz9jxGkPtoqioiOzsbJo3b05eXh6nnXYaoL4s/PnPfyYrK4uRI0eyc+dOdu3aBWAntACcfPLJbNmyhfz8fPbv38/QoUMBuPzyy+33WLx4sf28W7dudOzYke+//x6A4cOH07BhQ1q0aEHjxo3tv82ePXuyZcsWAH766SdatDisUXCJOL/obNq0yf58HTp0ICsry97XsmXLepV4kwgjMElk6tSppKenx2xLT09n6tSpNbIe57eouXPn0q1btxpZR1WzZ88errvuOm644YbD4kwnnHACd911Fw8//HANrc5wJHQMZuvWrQQCAdvamDVrFnv27GH58uXk5ubSqlUr2+rw+Xz2691uN6FQCClliTHG+LicE+e5XC6X/dzlchEKhew16vcujYMHD7JlyxZOOOEEIBqD2bhxI0uWLGHu3Ln2scXFxaSlpZXpvHUVIzBJZNy4ccycOZOOHTsihKBjx47MnDkzqVlk69evp3379vbj9ddfL/HYO++8k5NOOomsrCw+/vhjpk+fnrR11DT6m2+PHj0YOXIko0aN4t5770147HXXXcdnn33GDz/8UM2rNJSVxo0b89RTT/HYY48RDAbJz8+nZcuWeL1ecnJy2Lo1YTd4myZNmtC4cWMWL14MEJNYM2TIEPv5999/z7Zt2+jatWuZ19a9e3c7ceZIHDp0iIkTJ3LeeefRtGnTmH1t2rThoYce4sEHH7S3ff/99wm9DfUJk6acZMaNG5f0tORDhw4Bal5NMBg8bP/FF1+c8HXV5RLr0AEi3oSkna80wuFwifuGDRvGsGHD7OdpaWkmi6wM+Ny+pKcpl4fevXvTq1cvZs+ezbhx4xg9ejR9+/YlOzu7TNb3v//9b8aPH096ejqnn366vX3ixIlcd9119OzZE4/HwwsvvBBjuZTG2WefzT/+8Q+uvvrqhPuHDx+OlBLLsjj//PO5++67Ex533nnnMWXKFBYtWsSpp55KTk4OZ599dpnXURcRRzIfjzb69u0rzcAxg8EQz+DBg3nvvfeSUqMD4Pf7GTp0KIsXL8bjqfvf84UQy6WUfeO3GxeZwWAwlMLjjz+e1JqVbdu28dBDD9ULcTkS9fvTGQwGQxIYMGBAUs/XpUsXunTpktRz1kaMBWMwGAyGKsEIjMFgMBiqBCMwBoPBYKgSjMDUASrarn/KlCm0a9eO7OxssrOz+eCDD6p13VXJjh07OPfcc+nSpQvHHXcckyZNIhAIkJubG/M5p0yZwmOPPVaDKzUYjl5MkD+JPPLII/Tr14/hw4fb23Jycli2bBm33357lb53onb9ADfffHPNz0OREpxV1vHPy306yQUXXMD111/PnDlzCIfDTJgwgcmTJ9OjRw++/vprzjrrrCQsXNXbuN3upJyrtvPtrm8JWaGknc/j8tCrVa8jHlOX2/XPmDGDgQMHcu2119rHvfPOO8ycOZMPPvigxM+2Z88eLr/8cubNm5eUtdVmjAWTRPr168fYsWPtbq05OTmMHTuWfv36VdsanO36awVTpsDNNytRAfXvzTer7RVkwYIFpKamctVVVwHqJjVt2jT+9a9/cfvtt/Pqq6+SnZ3Nq6++CsCaNWsYNmwYxx57LE899ZR9nv/85z/079+f7Oxsrr32Wrt4s0GDBtxzzz0MGDDgsMaJ9ZlkiktZz1eX2/Vfeumlh3VDnj17NpdeeilQ8mdr0aIFbdq04fPPP6+eD1CDGIEpJ0dq1z9gwADatm3Lueeeyz333MPFF19M27Zt2bdvH1Az7fpnzJhBVlYW48ePt9dRbUgJ+/fD9OlRkbn5ZvV8//6o6JST77777rDuto0aNaJTp0785S9/4ZJLLiE3N5dLLrkEgHXr1vHRRx/x1Vdfcd999xEMBlm7di2vvvoqn3/+Obm5ubjdbrudSEFBASeddBJLly5l8ODBlbsGhjJT19r1jxw5knXr1vHTTz8BajTHp59+ynnnnXfEzwaqqr8650TVFEZgkkzTpk057bTT+Otf/8r48eMP60lUVSRq13/99dezadMmcnNzadOmDbfccku1rMVGCJg2DSZNUqLicql/J01S2yvoJiupsWFJ288++2x8Ph+ZmZm0bNmSXbt2MX/+fJYvX06/fv3Izs5m/vz5bN68GVAW0YUXXlihtRkqRl1s1+92u7ngggvsIWJz5861uzMf6bPB0dGqH0wMptwsXLjQ/lm369ekp6dz7733MnbsWO6++26eeeYZXnvtNTsmo9v1a5Ldrv+iiy6K2daqVSv752uuuYZzzjknae9XZrTIOBttVkJcAHr06HFYn7UDBw6wffv2hPGSkrrvXnHFFTHNBzWpqalHTdylptFNS7ds2cLJJ598WLv+zz77DJfLVaF2/R9++CGg2vXfeOONQMnt+hs2bHhYu34dxzxSu/5LL72U2267jUmTJjF79mx+97vflfrZ4Oho1Q/GgkkqOuby2muvcf/99/Paa6/FxGSqG226A7z99ts107lVu8WcOGMyFWDEiBEUFhba3zDD4TC33HILV155Ja1ateLgwYNlOscbb7xhT8DMy8srtWOvIfnU9Xb9gwYN4qeffuLbb7/liy++iEkuKemzwdHRqh+MwCSVZcuWxVgsw4cP57XXXmPZsmVJe4/ytOu//fbb6dmzJ1lZWeTk5DBt2rSkraNMOGMukyaBZUXdZZUQGSEEb7/9Nq+//jpdunThhBNOIDU1lb/97W8MHz6cNWvWxAT5E3HiiSfywAMPMGrUKLKysjjttNNiBNlQvdTVdv1CCMaOHcsVV1zBWWedRWpqaqmfTa+jvrfqB+MiSyqJUpGHDx8ek7ZcESrarv/ll1+u1PtWGiGgSZPYmIsWuSZNKuUmO+aYY+xkCSc+n++Igu5MW73kkkvsRAAn+nofbXhcnqSnKZeHutqu/9JLL+XRRx/loYceKtNnu/zyy4+KVv1g2vXHYNr1VxFJroMxGKqbZLfrHzJkCHPmzKm2JKCqxrTrN9Qc8WJixMVQx0hmu/49e/bwpz/9qd6Iy5EwFowDY8EYDAZD+TEWjMFgMBiqFSMwBoPBYKgSjMAYDAaDoUowAlMHEEJw+eWX289DoRAtWrSIqcyfN28e/fv3p1u3bmRnZ3PJJZckdYa4wWAwlBdTB1MHyMjIYPXq1RQVFZGWlsYnn3xCu3bt7P2rV6/mxhtvZO7cuXTv3h1QfZG2bNlChw4dqnx9j33+GAXBgqSdL8Obwa2DjjxioC63ee/bty/5+fnceOONdkfdU045hRkzZtC0aVO7oWPXrl0JBAL07duX5557Dq/Xy6pVq3j88cfthqvJZNehXVjSStr5XMJFqwatjnhMXfw9Llq0iOuuuw6v18uXX35ZYkW+8/fdqVMnvv76axo1asTIkSNZsGABHk/9v/0aC6aOcOaZZ/L+++8D8Morr9gtwQEefvhh/vznP9viAjBmzJjD2otXFckUl7Kery63eQf4/e9/z7HHHsumTZvYtGkTxx9/PFdeeaW9/7jjjiM3N5dVq1axY8cOu6Fiz5492bFjR5VYp8kUl7Kery7+HmfNmsWtt95Kbm5uudu9pKSkMGLEiCN2mahP1H8JTSI33XQTubm5ST1ndnY2Tz75ZKnH/eY3v+H+++/nnHPOYeXKlYwfP97uxvrdd9/V/FCxGmTgwIF2Y8JDhw5x7rnnsm/fPoLBIA888ADnnnsuW7Zs4cwzz2Tw4MF88cUXtGvXjjlz5pCWlsby5cvtCnBne/7i4mKuv/56vv76azweD0888QTDhw/nhRde4J133iEcDrN69WpuueUWAoEAL7/8Mj6fjw8++IBmzZqV2OZ948aNLF++POYmc88993Dcccexfv36w3pt9e/fP6bV++jRo5k9e3aVD7GrburC7/Ff//oXr732Gh999BGffvop11xzDY899hjvvfceADfccAN9+/aN+bIQz3nnncddd93FuHHjqu5i1hKMBVNHyMrKYsuWLbzyyitHnNa4d+9esrOzOeGEE46KUcF1sc277pXm7Njsdrvp3bs3a9eujTm2uLiYpUuXxghVfWz1Xld+j1dffTVjxozh0UcfrfA8l5NOOimp/QlrM8aCKQdlsTSqkjFjxnDrrbeycOFC9u7da2/v0aMHK1asoFevXjRv3pzc3Fwee+yxet1Tqy63eT/SPBvNpk2byM7OZsOGDVx00UVkZWXZ++pTq/e6/HusKG63m5SUFA4ePHjY7Jj6hrFg6hDjx4/nnnvuiZlYCarJ5tSpU2O+/RYWFlb38qqVutzmvUePHnzzzTdYVjRGYVkWK1eupE+fPkA0BrNx40aWLFnC3Llz7WPrU6v3uvx7BPB4PDG/x5KOi8fv9yfsvFzfMAJTh2jfvj2TJk06bHvPnj2ZPn06v/vd7+jWrRuDBg1i7dq1XHbZZTWwyuqlLrZ5P/744+nduzcPPPCAve2BBx5gxIgRh2X9tWnThoceeihmMFp9bPVeF3+PAB07dmTNmjX4/X7y8/OZP39+qefbu3cvLVq0wOv1lnkNdRXjIqsDJHJ1DRs2jGHDhtnPzz777Bpr/53hzUh6mnJ5qItt3p9//nluvPFGjj/+ePLz8+nXr1/C8QOggsJTpkxh0aJFnHrqqVXW6t0lXElPUy4PdfH3eMwxxzB27FiysrLo0qULvXv3LvV8OTk5CeOoUkI4rMYmhcPqEQqp52V56HNIqZ4HAlBcHD1HKKS2BQKHn1dKGDIE2rQp82UpE6bZpQPT7NKQbMrS5n39+vWcddZZ/P3vfz9iAgco18rQoUNZvHjxUVFHUZ2EwxAMRm/Efr/6NxiEsWMvY8aMZ0lLa2RvD4XUMX6/Osb5moICKCxUj6KiqHBYFsyb9xH9+w+gUaMm9s3deRuOn2yht7lc6l/nfv1c73e5oq9zHiul2ud2H35+vW/cOIh4aMtNSc0ujcA4MAJjSDZLly4lLS0tJkhfGTZs2MDOnTtjrFeDukkGAupmXlQUvbEXF6ubvt528GD0cehQdH9hoXq9lEosnNZDOAx79+5DCDcZGY0AdUP2eNTD7Y7+qx/6Zq6PS02F9HTw+UKsWfMNgwf3Iz0d0tLUPp9P/ew8pxYM/TwlBbze2PPr4+NfV90YgSkDRmAMhtpBOKysgIMH1b8FBUoQ9u+HvDw4cEDty89X2505LVocgsHow7LUN3cp1U1Y36idAqGtgAYNoFEjaNhQiYIWAp9PbWvQILpdi0NKivpZn/doG3lUksAYG9tgMFQrRUVKKPLzlUjs2wd79sDu3Uo88vKUoGiXEETjE6FQbJxBSnWD93qj29PSlAhkZESFIjVVbW/SRG3Tr9HbtSWRlnb0iUNVYgTGYDAkDSmVNbFvnxKRffvg55/VY/du+OUX5ZZyCocWCYjGKbS1AdFjPR7IzFQikZkJzZtD69bQtGlUODIy1HGG2kG5fxVCiMFAf2C1lPLj5C/JYDDUZrSI/PKLsjx+/hm2b4+KSFFR1PoIh5WloAVEB5TD4digtd8PjRurLKbWraFFCyUc6enKGmneXO2vifiCoeKUKjBCiK+klP0jP18D/AF4G7hXCNFHSvlQFa/RYDDUACq4rURDi8iOHfDTT0pgXK5otpJ2MzVrpmIQOm6i020PHVL7unWD9u2jVkezZurf9HT1OkP9oiwWjLMaaAJwmpRyjxDiMWAJUGaBEUKcAUwH3MC/4sVJCNEReB5oAeQBv5VS7ojsexjQyf9/lVK+Gtn+a+AxIAVYDvxeShkSQjSNnOs4oBgYL6VMTu9ug6EeIaVyZ+3ercTjhx/UY+fOqJvKsqIi0rSpsjCEUAF0nZFVVBQ9Z7t20LMndOqkfm7TRlkihqOLsgiMK3KzdqGyzvYASCkLhBChsr6REMINPA2cBuwAlgkh5kop1zgOewx4SUr5YkQ4HgQuF0KcDfQBsgEf8D8hxIfAIeBFYISU8nshxP3AFcBzwJ+BXCnl+UKIbpH3HlHW9RoM9ZFQKOrW2rYNvv9eiYl2a4XDKtidkQFt20atilAomrW1b5861rLUsV26wIknquObNYOWLVVWlcFQFoFpjLIMBCCFEK2llD8LIRpEtpWV/sBGKeVmACHEbOBcwCkwJwI3R37OAd5xbP+flDIEhIQQ3wJnRI7xSym/jxz3CXAXSmBORAkUUsp1QohOQohWUspd5VizwVBnsSxllezcCZs2wdq1ys2lq76FUFZF06ZKFCBaB1JQoNxhOtju9cKxxyoXV9u2USumcWOTdWUomVIFRkrZqYRdFnB+Od6rHbDd8XwHMCDumG+BC1FutPOBhkKI5pHt9wohngDSgeEoYfoF8Aoh+kopvwYuAo5xnOsCYLEQoj/QEWgPGIEx1EsOHlSi8MMP8N13sHGjKh4ElVnVsKESB2e1d0GBcouBEgrLUtZLhw5wzjlw/PFKSDIyjJAYyk+FE/qklIXAD+V4SaL/nvFVnrcCM4QQVwKfATuBkJTyYyFEP+ALYA/wZWS7FEL8BpgmhPABHwPabfcQMF0IkQusAr5x7IsuSogJqNhStYwXNhiSgWVF3Vxr1qjH3r1REWjQQAmDTtm1LOUG27s3WpQopcrYGj0aunePpvqaWhBDsqjOjPEdRK0LUNZEzFALKeWPKKuDiAvuQillfmTfVGBqZN9/gQ2R7V8Cp0a2jwJOiGw/AFwV2S5QYniYIEopZwIzQVXyJ+WTGgxJJhhU1snmzbByJaxfH21t4vMpYejQIerSKipSBYtFRdFU4Fat4KSToGNHFXjv2FGJisFQVVSnwCwDugghOqMsk98AMf3khRCZQJ6U0kLFUp6PbHcDTaSUe4UQWUAWylpBCNFSSrk7YsHcQVSEmgCFUsoAcDXwWUR0DIZaTyikBGXjRlixQgXjdeW6rgvR3d4tS2Vx/fSTEh0hlPWSnQ1du0Lnzso1ZgoQDdVNWepgBgJLZCWblkVSh28APkKlKT8vpfwukvn1tZRyLjAMeFAIIVEusj9EXu4FFkUGCh1ApS9rd9dtQohzUFluz0gpF0S2dwdeEkKEUfGa31dm/QZDVRIOw48/qmB8bq5yeem2KA0bqjRfZ52IlCqbKz9fCUqnTnDyyXDCCSpu0rRpjX0Ug8Gm1GaXQohnURlg3wPzgHlSyp+rYW3Vjml2aagupFTpwps2KQtl1aqoyysjQ7mu4i2OUCja5FFK5eIaNUrVmzRqVDOfw2CASjS7lFJeFzlBN+BM4AUhRGNUivA84HMpZTjJ6zUY6h1FRSqGsmoVLFumYiRSqip2p8vLSUGBCsyDEpwuXZSgdO8ejbkYDLWVMntlpZTrgHWojK00VKrwxcATwGHKZTAc7Uip3F7r18NXX6k4Cqg04WbNShYIy1LWTVGROu6CCyAry8RRDHWPCv13lVIWAR9EHgaDIUJRUTSOsmxZ1J3VqJHqwXWkZo3FxaowUkro3RtOP11ZLKbBo6GuYr4PGQyVQErVVXjdOliyRFXLg7I0mjU7crA9HFaiUlioXGHp6TBmDAwapLLADIa6jhEYg6GcWJZqubJqFXz+uUoP1m1XSrNSDhxQ2V+gjmvVSvXxGjBAucFMDy9DfcIIjMFQBkIh1YJlxQr44gslFEIoC6VjxyMH250pxa1bw9VXq1TizEzTot5Qv6nIwLEhqKr8+1Et8mdIKT9L9sIMhpomGFSFjsuWwZdfKneW262EobQ6EymVCO3fH00p/v3vVQaYiakYjhYqYsFcimqZ/ydgP6pdvhEYQ70gFFKisnSpEhW/X8VTWrQo3X0lpapl2b9fxVU6dIDTTlPV9Mcea1KKDUcfFRGYHsBBKeVuACFEfnKXZDBUL5al3F9Ll6qYSmFh2UXl4EFVz6I7ETduDD16qALILl2MqBiObioiMHejWvVrPqrsIqpz0mVl12qoH0ipgvNffQULF6r4iNutRKUsGVwHD6rssRYt4Kqr4Jhj1EyVjIwqX7rBUGeoiMDcBBQLIU4F5ksp367MAmpg0qXhKCY/H775BubPjw7UysxUMZLSCAZVU8n8fNXK5Q9/UP2/TKDeYEhMuQUmMoI4BegJnCOEOFNKOaUSa6juSZeGo4xAQDWPXLhQtbqXUgnEkVqtFBYqIfH7o0F5n081lTzjDBg+XI0LNhgMJVPRSv4Ayu20XAgxvZJrqO5Jl4ajACnVqODPP1fCUlSkChlLqlNxZn2ByhIbPFgF51u1Uq6wRo1MTMVgKA9lFhghRAMp5aGIhbECNd9lK9C5kmuo7kmX8Z/LTLSsRxQXqzYtH32kAvdut4qNtGoVe5yUUSslGFTb2reHs85SQ7latzZikmyklGzfvp22bdviMU3VahQpJZa0CIaDhGQIy7JIcaeQ6k2uWV6m33Jk3srTQoiLgSuBXsDJwGhUfKQyVOuky3jMRMu6j7ZWFi1S1orfr7K5SiqA3L9fZX61bQtDhkC3bqrw0Ux3rBzBYJA1a9awefNmtm7dyq5du/jll184cOAAu3fvZuXKleTl5fHtt9+SlZVV08utt0gpCYQDWNJCSolEIqUkLMOErTAhK6T2IREIZOT7vCWtmhEY4NdAdylloRBiqJTyf8DSJK2hWiddGuoPgQB8+y3Mm6caTHo8ylqJTy22LGWpHIjMM23eHG6/XaUTGyulbEgp2bdvH0IImjZtSigU4ssvvyQ3N5eCggLWrVvHu+++S15env0aj8dDZmYmjRs3plmzZlx00UVkZ2fTKt6cNFSKsBUmEA7Yj5AVQkQcQ/IwZ1Aspe2vLGUVmKZSysLIz+OA/yVrATUw6dJQx9mzBxYvhk8+UW6uRo0SWyuhkLJshFCTHs88U40P7tzZ9PwqjUOHDvHhhx+yYMECvvjiC9atW0cgEACgZcuWBINB9ummakCTJk0YPXo0Z555Jl27dqVjx440a9YMYRS8wkipLA/twrJk9CGRhKyQbY04LREou3BY0iJkhQhbYdwi+emQpU60BBBCfIByY72Jco9dGLEm6hVmomXtxbJUx+KPP1ZWixDKWkmUyeX3R6vpzz1Xtb039SmHs3XrVjwkn+AAACAASURBVL799lt27NjB+vXrWbNmDfv370dKyXfffUdxcTGNGjXilFNOoVevXrRp04ZwOMz69euxLIszzzyTwYMH06RJE3w+nxGTSqBdWWErTHGomOJQMWErbLuxoHLWhnaRBcNBAlaAoBW0xQkJYRmmUUojurXoVqHzV3iiZYQfUHUlfVCpwb8IId4Cpkopf6jQigyGMlBQoIoh339fFTamp6uixvhMML9fzVIJh9UM+5NOUhbLscfWzLprM0VFRTzwwAM88sgjhELK4E9PT+fEE0+kZcuWAAwcOJALL7yQU089Fbcp9EkKlmVRHC5WN3srbN/ktbgAh1kiUD5h0RZJMBy0zx+WYcIyDBIC4QD+sJ/iUDGBcMCOy0gkzVObV1hgSqKsAjMH2B+pkr9DCHEiKl34P8CgpK7IYAB+/hlycmDBApXllZmpalCchMNqYqSUkJam2rMMGgTt2pnYSjybN29m2rRprFixgrVr17Jv3z6uuOIKJk6cSPv27WndujUu04UzqUgp7Ru8P+SnMFhY+mvKICY6AywkQ3bQPmgpQdEWiX4eCAfwh/wEwspqKcka0udLNmUSmEg6cHfH8zXAGiHE+0lfkeGoRUrYsAE+/FBV27vdKl040az6Q4dg1y4lKiNGqOPM/fFw1q1bx6OPPspLL72E2+22LZPLLruM4cOH1/Ty6gXa/RQIqZu4RBIKh454Qy/reS1pKbFwuLacbi29XQuJ3p/ofaWUBK0g/pDftmL8IT/BcJCwDNOmQZvkXBAHZU5Gl1KuTbAtN7nLMRyNhEKqwn7OHNiyRbnBOnQoWTC0K+xPf1KjhQ2Hs2zZMqZOncqcOXNITU1l4sSJ3HHHHbRt27aml1an0WKig+6BcICiUFHCQDuUIYvLYY1YlqWEKmJthGXYPsYfVkIQDAcpChURCKuEC/2ezvfRWWX+sJ9ASK1Pi59AIIRQIhWHfr9kYqqdDDVGcbGKr8yZo+IrTZooN1gi95ZlqXhMXp4K7t98s7JaDIr9+/fz1ltvsWbNGlasWEFOTg5Nmzblnnvu4YYbbqCFmcFcLpzuLf1zcaiYkKXcSAkthDKIiXZbOUVEIOyix4Cl0owD4YAdR0lkkUgp1fHhAMWhYoqCRfjDfjtw7xLq21m8kOiamOrCCIyh2jlwQBVFvv++yvTKzFSpw/FIqepX9u9X1kz79mq08FlnqfHEBsjPz+fOO+/kpZdeorCwEJ/Px/HHH8+DDz7IxIkTadSoUU0vsdajg+4hGbID5P6QP/aYMlom2sLRKcT2Iy5eErbCysKIiE1Jbi07YB8OUhQsojBYaAfnhRC2OMWTaFtNYATGUG3s2aO6GH/6qXJxtWqlrBEnUkYtFSmVq+yiiyA7++gWFcuy2LFjB+3atbOzun788UfOPPNM1qxZw+9+9zsmTpxI7969TbC+BOxYSThgx0icdSRQvliJzgQLhoP4LT+hcAgLdS6ngBSHim3rQruonJlj8bERf9BPcbiYYDiIlNJO/z7MGpGyQrEdTfxaUj3J795qBMZQ5ezcqartFy9W7q82bWID91Kq+Sp6vHCrVnDBBdC3r/r5aM0Iy8vLY/369SxevJh//vOfbNiwgYYNG5KdnU1aWhqrVq3i4MGDvP/++4waNaqml1tr0DdMS1p2XEMHtOPbo9ivKcUiCVrBqChFrBFd9BiyQvjDfgoCBbZ1AeDCFRMf0UWTgVDAPr4oVEQoHHG7lRAbqaxLyxl3cQkXXreXVE8qqZ5UfB4fPrcPj8tDs7RmlXqfRBiBMVQJUqqA/XvvwfLlSlDat4+dnVJcrDLBpFS1LWecoWbWH+2NJteuXcudd97J3Llz7W2DBw/m+uuvZ+PGjaxcuZIDBw7Qo0cPHnroIU4++eQaXG3NYUlL3aSDRXZswSkU2iqBssdK7FhIJEius8IsyyrRInGeNz6jTKcnO1u4xFsw9usrISTO8wJ43V5S3Cn43D58Hh8pHvWz2+WOeY1et0CQ6j7KLJiqmHRpqFqkhPXrVeB+7Vo1QyU+IywQUJlgQsDFF0O/fmWbIlnf+PHHH3nuuef44osvCIVCFBUV8fPPP/PDDz+QkZHBX/7yFwYMGECPHj3onChIdRSg4xZhGY5xa0HpcYYjCYnTMtGCooPuum5Ei0xJ7Vh01b0/5KcoVERxsDia3VWSW6sSQfZ4EfG4PEpEPD68bi8+j49UTyoelyf2NY5MM7dw2+JjP1wpCCFo4kt+t9daKzBVMelSSnmgej/F0YOUSlDeegu+/15V0zv7gxUUwN690aLIIUNg9Gg1d+VoYt++fbz++uu8/fbbfPLJJ4TDYXr16kV6ejqpqan079+fcePGHRWZX84Ov04BCVvKxeTs+AsVqyWxK9sj2VtaLJwFkM40Xp3V5Xy/sKWO0y1cdM1JWIZxCVeVWSNa1HQb/TRPmm2JeFweW8TirSiPy4PXFRURr9uLR3hwu9wIIXALN27hxuPy4BIu+/iGKQ0rvOaSqLUCQ9VMunytOhZ+NCGlEpQ331SWS8OGKiNMCFXf8vPPKqDftCmcd54K1pc09Ks+sGnTJn755Rc6duzITz/9xMcff0wgEGDo0KF8++233Hvvvezbt4/jjjuOW2+9lauvvprjjz++ppddZTibKeqMqJClguFSyhLdWE5KExZLWnYWmLNFio6RhK2wnf6rLQ5tBTnP77RIdOqvThMuSUgqk63lwgUCW0i0RZKRkkFGSgYp7pQSXVr62BR3Cl6XF7fLTYorRVkvAvt8Wjy8Lq8tKCX1jKuKXnK1WWCqYtKlIUlICRs3KmFZu1YJi7OGpahIBfdHjlTWSqdO9UtUCgoK2LBhAzt37uTAgQMcOnSIt956i3nz5h12rNOtMWLECB5++GH69OlT75tDFgWLyC/Ot58nEoryWCU6SK5FRFskEglSpfQ6q9oTZW45g/KBUMAWFH/YX3VC4jinx+XB5/aR5k0j1ZNKiicqEJpEQuJz+2x3mEdERcTtcuMRHltM9EPXwdQ0tVlgkj7pMuGbmImW5WbzZiUsq1er1GGnsPj9Kh3ZsuAPf4BTTql/Afvc3FxGjRrFnj17Yra3adOG+++/n+zsbLZu3Urjxo0ZOXIkPp+PRYsWkZGRwYgRI+q9sATCAQ76D6o02wqm0eriRmdhYkiGQKJSecP+mB5bJRUj6qC8TgZwxkiqSkgAO1MrzZtmC4TX7Y353ccIifDgcSsxSfWkquNdXhDqvNoK0e4up4ssWRxtFkyVTLqMx0y0LDs7d6oYy9dfq3YuWlh0J2MdXxk2DIYOVZlh9YlgMMiyZcs455xzaNCgAa+++iodOnSgSZMmpKWl0bZtW7yJGqcB5557bjWvtupwVrVblhWTwaXdUuVB14HorCu7TUqkaaM/7KcoWKRSehMM09KtVPyhyHER11bIClVJjCReSHxuFVxP9abaQfb4+IiOJWlXVoonhVR3qu0G0zERr1u5svRzt3An5cYff83sGIzLhUd4bEsq2dRmgamSSZeG8rNnD8ydq6rvU1JU8N7lUtlgP/2kMsXOOkvFVzp1ik1FrquEw2Heffdd/vnPf7Jp0yZ27drF/v37AejUqRMLFiyol5ldzl5bYRm2my2GwtExuxX9lq/dU3YrlIiQhK2wHSdJ5N6CaPaVnd3lyNwKWkHlEpJgYR32nhXB2WpFWxBaRHTqb0lCkuKKZnZ5XB5SXCrQ7hIue3+qJzXGvZVsEXEJly1OXpeXFE8KHqHe50hxmGRTawWmCiddGsrIwYOqQHLePGWpHHOMEg8pVfA+GFQFkb/+dd0e6HXgwAFycnJYvHgxy5YtIy8vj927d7Nr1y46dOjAgAEDaNWqFZmZmWRmZnLhhRfSuh40QtN1JDpzK37AFVS8C3CiFGBdT6IbROq6knghsaRlB9oD4QDFwWjdSUk9tioqJM6MLa9Lpfraj0jcIz4+om/gHqEsDp2xpa0Zp0XicXns527hrtTNPd4a0q1itGvNzg5z154YTJkmWh4tmImWimBQWSuvvabcX23bqkJJy1JNKQsK1Aji8eNVVX5dIhQK8c033/Dll1+ybt06Vq1axZIlSwiFQqSkpNCnTx9at25Nw4YNGTNmDOeddx4eT639HlYhtEtpf9H+SrUa0a4np1WihUBbJMFw0C5MdApYvHurOFhMYbDQLkosybVVUZw3Zy0kzkr2FE+KfVN2CizE1pukuJVry+Vy2ZlazpTgiri0EqVhOwXE7XLb8Rev25s0t1kyqexES8NRgK5leekl5fpq3VrFVIqL1XMpoVcvNYK4a9e6kxW2detW3njjDRYsWMCiRYs4ePAgAI0bN6Zbt27ceuutnHHGGQwYMIDURDOY6yj2zTsU7Wul3VwVwc6+iriptAtNWxv+kN9OB05YSxKJk9jHR0RHu46cX3aTYZEIBD6Pj3RvejRry51SYnxEWyA+j0+l/Lo9MZaI8+F1ectsJTgF1S3cuF1u24XlcrnsLDKXcNnrry8YgTEAyjJ59VVYulTVrHTqpNxgu3erTLFzzoFTT60bFffhcJjly5fz+eef88EHHzB//nyklHTt2pVx48YxdOhQTj31VNq2bVuv/pghEucIhygKqWA3lN/NpYPuzlhJKBwiKNW/xeFiOyNLFyZqgQhZocOKEnW34GSmADuFRFen66wtLSTatRUfm9AWjK4h0dXsOmai3VoeceS6kfj16PcAYgoYdQC9OmMftQUjMEc5fj988gm8846Ks+ixxFu3KjfYuHEq9lJbrZWioiIWLFjAypUrAWWtvP322+zevRuA4447jnvuuYcrrriiXgTl9bd8HavQBYzxbd/LKipaFOwBVTo+EsnmsgUmzirRRYm6IFFnbh1poFVlhcQZbE9PSbddXPqbv/2Z4lqi6FiKx+WxrRPtaqrMjV9fZ6/LS3pKeowlYlAYgTlKkRK++w5eeEFZL23bqgyx/ftVS5fBg+HKK1WGWG0jEAjw4Ycf8p///If333+foqIie19GRgZnn302559/PkOGDKmTExxDluq462xtksi1VVILlYQFjY7JibaghPwUhgptq0S70pzncbq2ikPFdsC9JIukor22nJ2HdbzB2fFXWyTOz6xFRNeY6Gp1Z2Bdb9PWSmXiIzrOomMjOkPL605+em99wQjMUci+fTB7Nnz5JTRvrppR6uB927ZqWmR2dvVYLbNmzWLy5Mls27aNDh06MHXqVMaNGxdzjJSSgoIC9u/fzzvvvMMjjzzC9u3badGiBePHj2fMmDH86le/wuv14vF47HkpdY2wFeaA/wDFoeIyHV+SkOgCRecsd123EgwHE2Zv6XRkfygScA8VUhwsPmK/rYpYJM4sMLdw2zESp5DoG7j+jHY/rsgxzqwtr/DicXvsmpHKZmppPC4PqZ5UO724NlXH1yWMwBxFhELw2Wcq1hIOKzHZvVtZLb17q3Tj7t2rr45l1qxZTJgwgcLCQkC5tyZMmADAhRdeyIsvvsj777/PkiVLYqrmBw0axNNPP80ZZ5xRYmFjbcauNYlYJpZlURgsJGgFS39x3Dl0fEQH12MC6ZGhVfENHJ2ZW84piZB4Jkl5hSS+Jb12baV502wx0Vlb8RaJtjZ0HCW+51ZZemqVdY2AbenorDJnSrKh8pg0ZQf1OU35hx+UO+yHH6BxY1XjkpYGo0apXmHNm1f/mpo2bWoXLzpJSUmhadOm7Nq1iy5dujBo0CC6d+9OkyZN6NGjB7/61a/qXLBUC8Ih/yGKQkWHpcKW1lreklZMIaKuXSkOq/TeRA0cnbNNAuFATKykpHqS8hAvJF63lzRPGmneNNsascUgTkh8Hp9dye5xe+xGjdoySWYle3wWl16nz+Orc/+PaismTfkopaBABfA//ljVsbjd6nHFFapPWFpa9a7HsixWrlxJVlZWQnEBFWPJzs7mzjvvZOjQoXXiJqDTdQNhVfuhb9z6xpZoLkhJ59Edh4NWMKaWRPfV0s0ena/R+4qCRYelACd6//IKi3aTCSHsmhDdM8uZJaU/mzOLSlsjXrfXDrA7s7aS2Q7FaTFp95ae4lgb60fqO0Zg6ilSQm4u/PvfKiOsYUPVJv/886FPn9iRxVXN5s2badmyJQ0aNODll1/myiuvPOKI3w4dOiTsSlxbkFLasQwd13Cm65b5HFZUSOz03nCxHeTXsRLnubWQ+UN+CoKqgaMOzFfWvZXIIkn1pJLuTbcFRferiu+t5bRYYlxajuaMVZWmG5N27C57fYqh6qmzAlOZaZf1nbw8mDULFiyIFkdeeKEaR1wdgfs9e/YghCAzM5OvvvqKAQMG8Nprr9G5c2cWLVqEz+dj/vz5jB07lrlz51JcHA1qp6en87e//a3qF1lGtGvLDpBH3FPlmeuurRJnZ2AdH9Gi4kwx1q9xWiXOupNkNHB0nsPn9pGRkmHHRlJcKTHptjHV7+5oi3lnXMTtihYkJusG77zGAmFbI7pBo8etAvFGUGovdVJgKjPtsvpXW32EwyqI/8ILatb98cerOpb+/ZMTuC8p46ugoIADBw7Qpk0b9u3bR+vWrfnrX//Krbfeys6dOznxxBO56qqrKCgosOMr7du359VXX2Xy5Mk89NBDWJaFy+XipptuOiyLrDoJW2G7/bt2TSWqLSkpg0sH3p2t5LXbTKcCB8KBw7K3dAqwrkGx55ggKtXAsaQW8olcW/oYZ91IijuFdE+63dxRV6InKzai0bEZfW7tYkvxVKz9iqF2UCcFhspNu6yXbN8OM2aoVi8tW8Idd8CvfpW8OpYjZXzdcccdjBgxghdffJGGDRty44038sEHH/D444+Tl5eHz+dDSsm///1vxowZA6gAf05ODjNnzuTTTz9l+PDh5OTkMHbsWEaOHMnw4cOTs/A4dGGiXVsSCaBbRGMfiepLElkLzsB7cUgF27XrrDhUHDOnxJJWzLjdomARxeFiwla4xOJEXRdSFhK1kE9LSSPVHU3tFULEZGvp4xIJSVUOrtJtWdK8aXbLlKqYb2KoeepkFpkQ4iLgDCnl1ZHnlwMDpJQ3OI75L7BUSjldCHEB8CaQKaXcG3cu58Cxk7du3VpdHyMp+P2qlf7bb4PHAxdfrKZIJjt4X1LGV2pqKm63m3A4TKNGjdi7dy/hsJqpftZZZzFx4kSysrJo2rQpDRo0iHntI488Qr9+/WLEJCcnh2XLlnH77bcnZd3OMb26FqQ8M96dc+NDMtoGRXcDdloqAhW/CFgBO05SGCy037eyDRydMRKICEmkjiRRjMTZhNH5iB+lm+y4Rby1p91b+n21K85Qfygpi6yuCszFwOlxAtNfSnmj45i2wAygM6qV/4VADz2QLBF1LU157Vr4+98hPx8GDYJLL1V9xJLNBx98wNlnn13i/j59+rBixQpOP/10+vbtS2ZmJgMGDGDgwIHJX0wJOBs5aivCH/KXuXWKLlAMyZBdV6KzsYLhYEyBonMYlm7yqGeT6FThktqllIX4PlvOrC09qMo5e133u9IxETu7yyEmVZGt5VyrnTXmjsZk3MJtXwtD/aa+pSlXatplXefAAXj6aVi3TgnK1KnRHmLJQErJ4sWLOemkk2jatCl+v9+2UuJp3rw5ixYtYtmyZQwcOJCUlKr9ZuqcnKhjFroNPCRun5LIvaWFRLuzdMsUHYB3Bt219aPdW7oHl07bTSQkpX1xi8/YchYY6hhJqifVdmvpwVXa3eVxe1QtR5yQJLNuRFse+nw6yK7jNyY2YiiNuiowFZ52WZexLPjf/1SbF5cLrr9eBfCTmRkWDod5++23ufjiiznllFPo0KEDS5cuTSguaWlpTJ8+nfT0dIYOHZq8RUTQ7im7uDCSDnzE1yBjXhsMBwhG3GQqVlJEUIZimjXqNGA7Bdgxflfvr2jmlgsXiEhrFMe4XF3op2MgiOgYW51y63P5SPWmRpszRjr8VkWQPcWdgsvlwoUrpmWLwVAZ6qTAVHLaZZ3k55/h4YdVW5cBA1ShZLLiLJZlcfrpp9tjgbdv3w5gt2jp1asXF110Efn5+Tz//PN2xtfNN99c6YwvnY5rzyvBsqcrliYmEJtGrOMtBaECioPFtHvin4j8A6y94wpCMoyUFtkPvYi/YRorrh1ji0llJyYmio04e2c555Bol5Y9CdETsVrcaXjc0eFVlQ2wJ+oU4BKu2Op6IWwXnBETQ1VQJwUGQEr5AfBB3LZ7HD+/AbxR3etKNoGACuK/954aS/yXv8Bxx5X+utIC6MuXL2f27Nm0atWK//u//+OHH34AYOTIkTz44IMA9O7dmxNPPNF+7dixYyud8SWlpCBQQEGwwL5pl7VAMb6eRLu3dIt5ZxBfWhYt9u6m26yPOOg/wII/jmbwtDfp8vrnLLlwAFv3bYFyFibGN2rULUecc9rtyvFI3YiuWvd5fHaWlj1C1+XGRcXnjWhhsy0aRzzEOVTLJVx2nYrBUJ3UySB/VVGbgvy6nf4zz8ChQ6pn2EUXlS3teNasWdxyyy3s2rULl8uFZVk0b96cgwcPcsstt9C6dWsef/xxtm3bBsCIESMYOHAgQ4YMYeTIkQlveBXJ+NIuJW2RBK0gRcGiI1om8RaJbuBYGCw8LGvLkhZIbLHR1ojddl5KTn96HgPeXGKff8mFA/joD2eo4TcJcIqIjjfoGe12hXrE+nCOsdUikuZOw+f12a6ssloiiQREC4Mexaufm9iHobZRr7LIqoraIjD79sGDD6rssKZN4Y83Stq2UzeURx55hH59+zL817+2j3fe6OPrVY5E8+bN+cc//sEFF1xQ7huWJS27nsSeNRIRhLAVjplfUlJasC4yDIaDBEIBisJFdv2Idl/ZNSuWFTOXpCgUqWy3wjFxjsOQknt/fZ/99L4F9yIi3/DjXVrOSnUtJB4RaX0S2abdWSmelGg8xJExVRLxzR61lRFfc2IGVhnqIvUti6xeEgrBnDlqwmRaGtx1F3R6YQo8uh+mTQMh6Ne3L2vPOIPlKSncXlhIy5YtKSwsZNasWUyfPp2HH374iOLSvHlz1q1bh8/nIyMjA1cZMwS0iBSHiikIFMS0gIcj15RY0rIr3HV1vG2ROGbGa5EKhAL2WF7t9jriyF2UJeNETxQ57enYnmbn/mMhX916Kb5Iyq8OsuvMLJ3mm+ZRo3d1Zpaez14SiWIezjkmLlfUdWasD8PRghGYWsKKFfDf/0JxMZx+OowZA26XVFH96dPVQdOm0eaRRxgeDPJkMIgEdu3aRUpKCh9//DEzZswo9X3y8vLIzMxMuM+2EiKdeJ3V7glboySwSBIF2/1hf0z7FOesdp3JFQqHbNFypu/Gn9+J03LR1eHaEvG5Uhg07Q1OfHMJG357FmsnTyDroX+T/eIcWmS05Of7byfNmx7t6OtoUeIkXjicn90pSkIIO56iLRNjiRiOdozA1DD79qnsML9f9Q4bP96ZHSaU5QJKZKZPpxvwJNEeOKDa25dFXABatmwZ8zxkhexmimXpCBwTIwmHKAoXURAssAsOnfNHYkQkWETAMdtdAFYJoqVdV1psIBoTsZsuRjKz9MPulRWJWfjcPhq36sjea34Lf7ubnm4vriencyg9k8aNm+Bu1D5mRogz5qHdV846EC00Oohu4iAGQ+kYgakhiopUx+MVKyA9HTIy/o8fd+6ge/f/2s0kL7v0Uor9fho1bswUx2tvLuGce3bvZt5HHzFhwoSYOfUan8/HkKFDOOg/aE9CjA+46xu8bo+iq9ULg4WqADEyK74oFO3w66ycDzpFJGKJDHlhAamHiu3gupQWo56eh79BGp9dpWJJWlR0Nbg90TBiIehsrRS36varg++6nsTOyorc/AHkffcTkJIGIjo/XUx7Eo/bSzMjFAZDlWMEppoJBlWh5NKl6vlFF8GwYZAz9BVWLV6M7oS2detWWj70EPuB+O4v04D9QBOiYtOqZUsyH3iA3zZtivjnP/nzn//Mtm3b7Ar89se058577+TCsRdy0H8wGheJiIOzt5Z2WxUHi2PbykcC8Lpy3iVcIDms2y9E3GcSBJB2yM+AN5ciECz442hGzZjHyW8uZdVvfk27hm1JiWRp6ap0/UjzqBG7KS4lNB4RcTtF9EALhMvlOswK8Qg3Llc0dRcpS8wcMxgMVYPJInOQlCwyx43skUceIT8vj1mzZ7N16480a/ZHTjjuLLypDbjttv6cemoRb74xCwEEbriB6/x+2/01DbgJ+AloA/DHP8KTT/Jh166cuWEDB4GGKHfZHV4vb3Y8hnM2bibQtw+/fDoXhFAWiFSBdd10sThcHDOsSlsd2jLRvba0SwoSZGdJaWdiAUjLQjoq0e10Xo8Xn8uHz51Cv8f+y/Evv2+f4uff/4a9UyeT6kmzA+i68E+LhJ6/7qzt0EOunO6zw5gyRcWuIokRSAk33wxNmqh9BoMhqZgssirErhH53//sG1vOwoWM+tvfaJufj5ds7uNKbszLYWLek+S1yuS7wNN06zqRd3bvBmAe8A1KVG6KnPcb1LQ0gEA4yKHCPIaeNgw2bKBICBpKqY4PBmHjZgB+7tGBb35cRmG4yO7oq4VDdxV2ZoAlirfo2IMOtOubutfl5dR/zye9IKAysbyp+Nwp9H30v4jGTdh5yzXRVvGetJiqcdffh8LL7aPvMW0abSLpv27htjO1Kh0Yl4cnRnDzzer5pEnGkjEYqpF6bcGUNvUynopaMDk5OZw7ZgyPWRYTCguZmZ7On6zm/FC8lxYUspt0WtGTQywlA9gNtAJeatKEyyMt8Fe53fRM0O9rucdN9xOOI33N9/a2vK4d2HJqFp6Dh8h6daG9fckFA/johpKLCONxzuLwuKPNFHWjxVRvKhnejOiAKpeX9lMep+k/XuLgdeM5+PD9NL7zPjKe+SdFE6+l+NG/4XIUF9r1HoDrT7cgnnoq+uaTJkUtjGSjLRYtMlX9R3C7UAAACvNJREFUfgbDUc5RZ8GUceplUti+fTv+QIBrAwEKacdNhTuZgKpF2Y2blhQiWWof35JI2UZEXPK7H0vPtZsPO+83wOZjW3KyQ1wAmq3fxvc9WhNIkELrFm7cbo9tceiMKq/HS6pbdepN96ST5k2zA+mp7lQ1S8SRsmuPp9VxjUisw/X080hfExo+9RQNn430D500ibRp00hLdPPWN/unnore5J03/6q46YtI9p1TYIy4GAzVTr21YIQQA4EpUsrTI8/vApBSPljSaypqwTRu3JgDB84CuqKcTvdH10FsDeCTRF1gELE6Jo7i3pF/TXjut3q051fFQVpv2mVv23vCMezu05Xusz897Pi91/yWvKl340tJjQmOe93emEwrZ6+qcmdRSRnbwtmyjnzzru6YiLFgDIZqpSQLRvnZ6+EDuAjlFtPPLwdmJDhuAvA18HWHDh1kRQAknCFhgJymbm/2Y1cpz9ddOkpu+O3ZMdv0w9/zRBno2yfhPqtVKylBhm+8UVrhsJR//KPaN2CAlJZVoc9RJixLykmTYtczaVLp7xm/v6rW6FyfXlf8c4PBkFRQXewPuw/XWxcZJPAfHdZQBKSUM4GZoCyYirxRq1at2LVrnp35pTPBdqPcYb+4XLS0LHYLQcs4i7HrKx+rH7KzITc35t+U3Fw18EVv02RnIxo1gksuwfXkk+pb+ZNPqn1Nm1bdt3SnZVBed1f89qpaoxDKMnJaLLpYtUkTY8EYDNVIfe5lUerUy2QxfPhwfD4f+4mKi8/nY09mJmRmkjl5MpZlkdmvn3pBu3bKrfTHP6rnrVvD6NHqprh8ufp3zBi1XwglLpMmqddMmqSe9+4de0PXIlOVabgl3bwnTapdN+8pUw6/NtOmmRRlg6Gaqc8xGA/wPTACNfVyGXCZlPK7kl5TmTqYWbNmMXnyZLZt3UqHjh2ZOnWqGsYVDoM70iRxyhTVG0ZbHVLCTTepm/N990VTaJ3/3ndf7avpiE/1Nam/BsNRzVHZrl8IcRbKqNBTL6ce6fhqaddfkZuzuaEbDIZazFGXpgyJp17WOBWJRVRX/MJgMBiSSL22YMqLEGIP2O3Aykom8EsVLCeZmDUmB7PGylPb1wdmjRWho5SyRfxGIzCVRAjxdSLTsDZh1pgczBorT21fH5g1JpP6nEVmMBgMhhrECIzBYDAYqgQjMJVnZk0voAyYNSYHs8bKU9vXB2aNScPEYAwGg8FQJRgLxmAwGAxVghGYSiCEOEMIsV4IsVEIcWctWM8xQogcIcRaIcR3QohJke3NhBCfCCE2RP6Nn8JcE2t1CyG+EUK8F3neWQixNLLGV4UQKTW8viZCiDeEEOsi13NgbbuOQoibI7/n1UKIV4QQqTV9HYUQzwshdgshVju2JbxuQvFU5O9npRCiTw2u8dHI73qlEOJtIUQTx767ImtcL4Q4vabW6Nh3qxBCCiEyI89r5DqWBSMwFcQxb+ZM4ETgUiHEiTW7KkLALVLK7sApwB8ia7oTmC+l7ALMjzyvaSYBax3PHwamRda4D/h9jawqynRgnpSyG9ALtdZacx2FEO2APwJ9pZQnobpV/Iaav44vAGfEbSvpup0JdIk8JgDP1OAaPwFOklJmoVpM3QUQ+fv5DdAj8pr/i/zt18QaEUIcg5pxtc2xuaauY6kYgak4/YGNUsrNUsoAMBs4tyYXJKX8SUq5IvLzQdRNsV1kXS9GDnsROK9mVqgQQrQHzgb+FXkugF8Db0QOqdE1CiEaAUOA5wCklAEp5X5q2XVEdeJIi/TdSwd+ooavo5TyMyAvbnNJ1+1c4KVIx/clQBMhRJuaWKOU8mMpZSjydAmqOa5e42wppV9K+QOwEfW3X+1rjDANuJ3YzvA1ch3LghGYitMO2O54viOyrVYghOgE9AaWAq2klD+BEiHUFIGa5EnUH4kVed4c2O/4A6/pa3kssAf4d8SN9y8hRAa16DpKKXcCj6G+yf4E5APLqV3XUVPSdautf0PjgQ8jP9eaNQohxgA7pZTfxu2qNWuMxwhMxSnTvJmaQAjRAHgTuElKeaCm1+NECHEOsFtKudy5OcGhNXktPUAf4BkpZW+ggNrhVrSJxDHOBToDbYEMlKsknlrxf7IEatvvHSHEZJSreZbelOCwal+jECIdmAzck2h3gm214vduBKbiVNu8mfIghPCixGWWlPKtyOZd2mSO/Lu7ptYHDALGCCG2oNyKv0ZZNE0irh6o+Wu5A9ghpVwaef4GSnBq03UcCfwgpdwjpQwCbwG/onZdR01J161W/Q0JIa4AzgHGyWj9Rm1Z43GoLxPfRv522gMrhBCtqT1rPAwjMBVnGdAlkrWTggoEzq3JBUViGc8Ba6WUTzh2zQWuiPx8BTCnutemkVLeJaVsL6XshLpmC6SU44Ac1JhrqPk1/gxsF0J0jWwaAayhFl1HlGvsFCFEeuT3rtdYa66jg5Ku21zgd5EsqFOAfO1Kq26EEGcAdwBjpJSFjl1zgd8IIXxCiM6oQPpX1b0+KeUqKWVLKWWnyN/ODqBP5P9qrbmOh5FojrJ5lO0BnIXKONkETK4F6xmMMo1XArmRx1moGMd8YEPk32Y1vdbIeocB70V+Phb1h7sReB3w1fDasoGvI9fyHaBpbbuOwH3AOmA18DLgq+nrCLyCigkFUTfB35d03VCunacjfz+rUBlxNbXGjag4hv67edZx/OTIGtcDZ9bUGuP2bwEya/I6luVhKvkNBoPBUCUYF5nBYDAYqgQjMAaDwWCoEozAGAwGg6FKMAJjMBgMhirBCIzBYDAYqgQjMAZDBRFChIUQuZGOxt8KIf4khKjQ35QQ4gNnB99kU9XnNxgSYdKUDYYKIoQ4JKVsEPm5JfBf4HMp5b01uzKDoXZgLBiDIQlIKXejWqXfEKmodkdmjCyLzOi4FlSrFCHEZxHLZ7UQ4tTI9i2O+R53R2aTfCLUnJdbI9sXCiEeFkJ8JYT4Xr/WSWnnF0JcF9mXK4T4QQiRE9k/SgjxpRBihRDi9Ug/O4OhUhiBMRiShJRyM+pvqiWqOjxfStkP6AdcE2k1chnwkZQyGzVnJtd5DiFEX+BCVCfsC4C+cW/jkVL2B24CEllKRzy/lPLZyL5+qArxJyLC9hdgpJSyD6qDwZ8qdhUMhiie0g8xGAzlQHe2HQVkCSF0X7DGqD5Wy4DnI01J35FS5sa9fjAwR0pZBCCEeDduv25guhzolOD9Szu/ZjqqD9y7kQ7XJwKfq7ZmpABflvpJDYZSMAJjMCQJIcSxQBjVLVgAN0opP0pw3BDUwLWXhRCPSilfcu4u5W38kX/DJPj7lVJ+Vsr5EUJcCXQEbnC85ydSyktLeW+DoVwYF5nBkASEEC2AZ4EZUmXOfARcH7EkEEKcIITIEEJ0RM3D+Seq83X8/PTFwGghRGokDnJ2OddxxPMLIU4GbgV+K6XUA9+WAIOEEMdHjkkXQpxQnvc1GBJhLBiDoeKkCSFyAS9qSNXLgB6T8C+UC2tFpJ3+HtSo4GHAbUKIIHAI+J3zhFLKZUKIucC3wFZUPCS/HGs64vlRVkszICfiDvtaSnl1xKp5RQjhixz3F1SncIOhwpg0ZYOhliGEaCClPBSZYvgZMEFKuaKm12UwlBdjwRgMtY+ZQogTgVTgRSMuhrqKsWAMBoPBUCWYIL/BYDAYqgQjMAaDwWCoEozAGAwGg6FKMAJjMBgMhirBCIzBYDAYqgQjMAaDwWCoEv4fiH/M4JZDzSgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAC1CAYAAACAovJwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2dd5gUVdaH39txEgxIkiBBBEEEBiWIqIAoioKYVlRWxYBhVxaMu4gBA4qIIi66fqy6ri4G1FVZV1QkLIKCBAcBAxIGHCSHgQk9ne73x+1bUz059PQE7vs89fR0VXXV7erp+vUJ9xwhpcRgMBgMhnjiqOkBGAwGg+HYw4iPwWAwGOKOER+DwWAwxB0jPgaDwWCIO0Z8DAaDwRB3jPgYDAaDIe64anoAdYGmTZvK9u3b1/QwDAaDoU6xZs2a/VLKZsVtM+JTDtq3b8/q1atrehgGg8FQpxBCbC9pm3G7GQwGgyHuGPExGAwGQ9wx4mMwGAyGuGPEx2AwGAxxx4iPwWAwGOKOER+DwWAwxB0jPgaDwWCIO0Z86jApKSkAZGRkcOqppxbZPmbMGDp06EBaWhqnnXYa33zzDQAPPfQQPXr0IC0tjaFDh/Lbb7/FddzVidPpJC0tjW7dutGzZ0+ee+45wuEwAEuWLCE1NZVevXrRpUsX7r333hoercFw7GImmcaJOXPmMGnSJHbs2EHbtm2ZMmUKo0ePrvbzPvPMM1x55ZV88cUX3HbbbXz//ffcd999PP744wC88MILPPbYY7z88ssxP/fDD8OOHbE7Xtu28Nhjpe+TmJhIeno6AHv37uXaa68lKyuLRx99FICzzz6bTz75hLy8PHr16sVll13GgAEDYjfIesq8n+eRH8qP2fG8Ti+XnHxJqfs4nU66d+9OMBikQ4cOvPnmmzRq1KjK587IyGD48OFs2LChyscCmDBhApdffjnnnHMOgwYNYvr06fTu3Zv27dvToEEDAEKhEJdffjkPPfQQXq+XjIwMunbtysknn4zf76d37968+uqruN1u1q9fz7PPPsvrr78ek/HVVozlEwfmzJnDrbfeyvbt25FSsn37dm699VbmzJkTtzGcc845bN68GYCGDRta63NychBCVMs5d+yA9u1jt1RUyJo3b87s2bOZNWsWhTv2JiYmkpaWxs6dOyv79o4pYik85T2e/iGxYcMGjjvuOF588cWYjiEWHDx4kBUrVnDOOecUu33x4sWsX7+eb7/9lq1bt3Lrrbda2zp27Eh6ejrr168nMzOTuXPnAtC9e3cyMzPZEctfbrUQIz4xYtCgQdYvlUAgwKBBg/jXv/4FwMSJE8nNzY3aPzc3lwkTJgCwf/9+Bg0axH/+8x8Adu/eHfPx/ec//6F79+7W80mTJnHCCScwZ84cHivLnKjDnHjiiYTDYfbu3Ru1/tChQ/zyyy8l3jQMtYv+/ftbPxSys7MZMmQIp512Gt27d+fjjz8GsKyJsWPH0q1bN4YOHUpeXh4Aa9asoWfPnvTv3z9KxHw+HzfeeCPdu3enV69eLF68GIDXX3+dSy+9lBEjRtChQwdmzZrFc889R69evTjjjDM4ePAgAO+//z4XXnhhmeNPSUnh5Zdf5qOPPrJeq3E6nfTt2zfqh9CIESN45513qnDFaj9GfOJAZmZmsesPHDhQ7ee+7777SEtLY/bs2bz66qvW+ilTpvDrr78yevRoZs2aVe3jqEnsVs9XX31Fjx49OP744xk+fDjHH398DY7MUB5CoRALFy7kkkuUmy4hIYEPP/yQtWvXsnjxYu655x7rM/7ll1/44x//yMaNG2nUqBEffPABADfeeCMvvPCCFffUaCFav349b7/9NjfccAM+nw+ADRs28NZbb/Htt98yadIkkpKS+O677+jfvz9vvPEGAMuXL+f0008v1/to2LAhHTp04Jdffola7/P5WLlyZZSI9e7dm6+++qqil6pOYcQnRixZsoQxY8YA4Ha7WbJkCb///e8BaNu2bbGv0eubNm3KkiVLGDFiBEBMb4jPPPMM6enpLFiwoNikhGuvvdb6gtZHtm7ditPppHnz5oCK+Xz//fesX7+ev/3tb1Z8yFD7yMvLIy0tjSZNmnDw4EHOP/98QP2YeOCBB+jRowfnnXceO3fuZM+ePQBWgg3A6aefTkZGBllZWRw+fJiBAwcCcN1111nnWLZsmfW8S5cutGvXjk2bNgEwePBgGjRoQLNmzUhNTbW+n927dycjIwOAXbt20axZsUWbi8X+Q2jLli3W+2vbti09evSwtjVv3rxeJQIVhxGfODBlyhSSkpKi1iUlJTFlypQaGY/9l9e8efPo0qVLjYyjutm3bx+33347d955Z5G4VufOnZk4cSJPP/10DY3OUBY65rN9+3b8fr9lpcyZM4d9+/axZs0a0tPTadGihWWteL1e6/VOp5NgMIiUssS4ZuFYoB37sRwOh/Xc4XAQDAatMepzl8XRo0fJyMigc+fOQEHMZ/PmzaxYsYJ58+ZZ+/p8PhITE8t13LqKEZ84MHr0aGbPnk27du0QQtCuXTtmz54d02y3n3/+mTZt2ljLe++9V+K+f/nLXzj11FPp0aMHX3zxBTNnzozZOGoa/Wu5W7dunHfeeQwdOpRHHnmk2H1vv/12li5dyrZt2+I8SkNFSE1N5YUXXmD69OkEAgGysrJo3rw5brebxYsXs317iVX7AWjUqBGpqaksW7YMICrR55xzzrGeb9q0iR07dnDyySeXe2xdu3a1EnlKIzs7mz/84Q9ceumlNG7cOGpby5YtmTp1Kk899ZS1btOmTcV6KuoTJtU6TowePTrmqdXZ2dmA6jcUCASKbP/d735X7Ovi5WZr2xYi3omYHa8sQqFQidsGDRrEoEGDrOeJiYkm262ceJ3emKdaV4RevXrRs2dP3nnnHUaPHs2IESPo3bs3aWlp5bLc//GPf3DTTTeRlJTEBRdcYK3/wx/+wO2330737t1xuVy8/vrrURZPWVx88cX83//9H7fcckux2wcPHoyUknA4zGWXXcZDDz1U7H6XXnopkydP5quvvuLss89m8eLFXHzxxeUeR11ElGZ2GhS9e/eWppmcwWAojrPOOotPPvkkJnOQAPLz8xk4cCDLli3D5arb9oEQYo2Usndx24zbzWAwGKrAs88+G9M5OTt27GDq1Kl1XnjKon6/O4PBYKhm+vXrF9PjderUiU6dOsX0mLURY/kYDAaDIe4Y8TEYDAZD3DHiYzAYDIa4Y8SnDlPZlgqTJ0+mdevWpKWlkZaWxqeffhrXcVc3mZmZjBw5kk6dOtGxY0fGjx+P3+8nPT096r1OnjyZ6dOn1+BIDYZjF5NwEAemTZtGnz59GDx4sLVu8eLFrFq1ivvvv79az11cSwWAu+66q+b72UgJ9pnnhZ9X6pCSyy+/nDvuuIOPP/6YUCjErbfeyqRJk+jWrRurV6/moosuquLAFaFQCKfTGZNj1XbW7VlHMByM2fFcDhc9W/QsdZ+63FJh1qxZ9O/fn9tuu83a76OPPmL27Nl8+umnJb63ffv2cd111/HZZ5/FZGy1GWP5xIE+ffpw1VVXWRVzFy9ezFVXXUWfPn3iNgZ7S4VaweTJcNddSnBAPd51l1pfBRYtWkRCQgI33ngjoG5gM2bM4JVXXuH+++/n3XffJS0tjXfffReAH374gUGDBnHiiSfywgsvWMf517/+Rd++fUlLS+O2226zJq+mpKTw8MMP069fvyJFKuszsRSe8h6vLrdUuOaaa4pUpX7nnXe45pprgJLfW7NmzWjZsiXLly+PzxuoQYz4xIjSWir069ePVq1aMXLkSB5++GF+97vf0apVKw4dOgTUTEuFWbNm0aNHD2666SZrHHFDSjh8GGbOLBCgu+5Szw8fLhCkSrBx48YiVYYbNmxI+/btefDBBxk1ahTp6emMGjUKgJ9++onPP/+cb7/9lkcffZRAIMCPP/7Iu+++y/Lly0lPT8fpdFolWHJycjj11FNZuXIlZ511VuWvgaFC1LWWCueddx4//fQTu3btAlQLlS+//JJLL7201PcGqtpBPHt91RRGfOJE48aNOf/883n88ce56aabitR3qi6Ka6lwxx13sGXLFtLT02nZsiX33HNPXMZiIQTMmAHjxyvBcTjU4/jxan0VXG8lFZEsaf3FF1+M1+uladOmNG/enD179rBw4ULWrFlDnz59SEtLY+HChWzduhVQltQVV1xR6fEZKk5dbKngdDq5/PLLrQZx8+bNs6pkl/be4NhopwAm5hMzlixZYv2tWypokpKSeOSRR7jqqqt46KGH+Nvf/sbcuXOtGJBuqaCJdUuFK6+8MmpdixYtrL/Hjh3L8OHDY3a+cqMFyF7UtIrCA9CtW7citeuOHDnCr7/+Wmx8pqQqyDfccENUoUdNQkLCMRPnqWl0kdiMjAxOP/30Ii0Vli5disPhqFRLhfnz5wOqpcK4ceOAklsqNGjQoEhLBR07La2lwjXXXMN9993H+PHjeeedd7j++uvLfG9wbLRTAGP5xAUd45k7dy6PPfYYc+fOjYoBxRvtCgD48MMPa6Z6rna12bHHgCrJkCFDyM3NtX6ZhkIh7rnnHsaMGUOLFi04evRouY7x/vvvW91PDx48WGblZEPsqestFQYMGMCuXbtYt24dX3/9dVSiS0nvDY6NdgpgxCcurFq1KsrSGTx4MHPnzmXVqlUxO0dFWircf//9dO/enR49erB48WJmzJgRs3GUC3uMZ/x4CIcLXHBVFCAhBB9++CHvvfcenTp1onPnziQkJPDkk08yePBgfvjhh6iEg+I45ZRTeOKJJxg6dCg9evTg/PPPjxJsQ3ypqy0VhBBcddVV3HDDDVx00UUkJCSU+d70OOp7OwUwbre4UFw69eDBg6NSrytDZVsqvPnmm1U6b5URAho1io7xaAFs1KjKrrcTTjjBSt6w4/V6SxV8e+rtqFGjrKQEO/qaH2u4HK6Yp1pXhLraUuGaa67hmWeeYerUqeV6b9ddd90x0U4BTEuFcmFaKlQT1TDPx2CIN7FuqXDOOefw8ccfxy0pqToxLRUMtZPCQmOEx1AHiWVLhX379nH33XfXC+EpC2P5lANj+RgMBkPFMZaPwWAwGGoVRnwMBoPBEHeM+BgMBoMh7hjxqcMIIbjuuuus58FgkGbNmkVVLPjss8/o27cvXbp0IS0tjVGjRsW037zBYDBUBjPPpw6TnJzMhg0byMvLIzExkQULFtC6dWtr+4YNGxg3bhzz5s2ja9eugKoxlZGRQdu2bat9fNOXTycnkBOz4yW7k7l3QOltIOpyGf7evXuTlZXFuHHjrKrGZ5xxBrNmzaJx48ZW4cyTTz4Zv99P7969efXVV3G73axfv55nn33WKm4ba/Zk7yEswzE7nkM4aJHSotR96uJn+dVXX3H77bfjdrv55ptvSqxUYP/M27dvz+rVq2nYsCHnnXceixYtwuWq/7dmY/nUcYYNG8Z///tfAN5++22rZDvA008/zQMPPGAJD8All1xSpPx7dRFL4Snv8epyGX6Am2++mRNPPJEtW7awZcsWTjrpJMaMGWNt79ixI+np6axfv57MzEyrcGX37t3JzMysNqs2lsJT3uPVxc9yzpw53HvvvaSnp1e4RI7H42HIkCGlVt+oT8RVXoUQFwIzASfwipRyaqHt7YDXgGbAQeD3UsrMyLanAT3t93Ep5buR9ecC0wEPsAa4WUoZFEI0jhyrI+ADbpJSboi8ZjwwFhDA36WUz1flfU2YMIH09PSqHKIIaWlpPP982cO6+uqreeyxxxg+fDjff/89N910k1URd+PGjTXfMK4G6d+/v1UAMjs7m5EjR3Lo0CECgQBPPPEEI0eOJCMjg2HDhnHWWWfx9ddf07p1az7++GMSExNZs2aNNSve3j7B5/Nxxx13sHr1alwuF8899xyDBw/m9ddf56OPPiIUCrFhwwbuuece/H4/b775Jl6vl08//ZTjjjuuxDL8mzdvZs2aNVE3n4cffpiOHTvy888/F6lb1rdv36hS/CNGjOCdd96p9gaFNUFd+CxfeeUV5s6dy+eff86XX37J2LFjmT59Op988gkAd955J7179476MVGYSy+9lIkTJzJ69Ojqu5i1hLhZPkIIJ/AiMAw4BbhGCHFKod2mA29IKXsAjwFPRV57MXAakAb0A+4TQjQUQjiAfwJXSylPBbYDN0SO9QCQHjnW9SjRQwhxKkp4+gI9geFCiE7V866rnx49epCRkcHbb79daofOAwcOkJaWRufOnY+J1tF1sQy/rjtnr5rtdDrp1asXP/74Y9S+Pp+PlStXRolYfS3FX1c+y1tuuYVLLrmEZ555ptL9eE499dSY1nyszcTT8ukLbJZSbgUQQrwDjAR+sO1zCqBLHS8GPrKt/5+UMggEhRDrgAsj++RLKTdF9lsATARejbzmKQAp5U9CiPZCiBZAV2CFlDI3Mo7/AZcB0yr7xspjoVQnl1xyCffeey9LlizhwIED1vpu3bqxdu1aevbsSZMmTUhPT2f69On1uj5ZXS7DX1ovIs2WLVtIS0vjl19+4corr6RHjx7WtvpWir8uf5aVxel04vF4OHr0aJHeP/WNeMZ8WgO/2p5nRtbZWQfoTl2XAQ2EEE0i64cJIZKEEE2BwcAJwH7ALYTQM2ivjKzXx7ocQAjRF2gHtAE2AOcIIZoIIZKAi2yvsRBC3CqEWC2EWL1v374qvO3q56abbuLhhx+O6lQKqqDplClTon415+bmxnt4caUul+Hv1q0b3333HeFwQTwkHA7z/fffc9pppwEFMZ/NmzezYsUK5s2bZ+1b30rx1+XPEsDlckV9liXtV5j8/PxiK2DXN+IpPsV9+oU/+XuBgUKI74CBwE4gKKX8AvgU+Bp4G/gmsl4CVwMzhBDfAkcBXXp3KtBYCJEOjAO+i7zmR+BplJX0GUqkipTrlVLOllL2llL2juUvm+qgTZs2jB8/vsj67t27M3PmTK6//nq6dOnCgAED+PHHH7n22mtrYJTxpS6W4T/ppJPo1asXTzzxhLXuiSeeYMiQIUWyE1u2bMnUqVOjGt7V11L8dfGzBGjXrh0//PAD+fn5ZGVlsXDhwjKPd+DAAZo1a4bb7S73GOoq8XS7ZRJtYbQBonwEUsrfKLBWUoArpJRZkW1TgCmRbW8Bv0TWfwOcHVk/FOgcWX8EuDGyXgDbIgtSyldRrjmEEE9GxlbnKM59NmjQIAYNGmQ9v/jii2usPHuyOznmqdYVoS6W4X/ttdcYN24cJ510EllZWfTp06fY9hCggtOTJ0/mq6++4uyzz67WUvwO4Yh5qnVFqIuf5QknnMBVV11Fjx496NSpE7169SrzeIsXLy41dlsRQiHw+yEQUH8Hg5CfH/3c7wefD/Lyohe/X+17+DC0bQu2JNqYEbfCokIIF7AJGIKyaFYB10opN9r2aQoclFKGhRBTgJCU8uFIskIjKeUBIUQP4C0gLZLV1lxKuVcI4UVZR1OklIuEEI2AXCmlXwgxFjhbSnl95Dz6NW2BL4D+UspDJY3dFBY1xJrylOH/+eefueiii/jrX/9a5g0pPz+fgQMHsmzZsmNijkhtIhYtFcJhdbP/3e+uY+LER2jX7iQCgQIR0IKQlwdHj0JWlhKGw4fhyBElJvYlHC7oyShEwaJx2LRfb3M4ovcTApxO6NoVbrutcu+rtMKica1qLYS4CHgelWr9mpRyihDiMWC1lHKeEOJKVJKABJYCf5RS5gshEoC1kcMcAW6XUqZHjvkMMBzlQvybTpsWQvQH3gBCqKSGm7XACCG+ApoAAeBuKWWp9rARH0OsWblyJYmJiVEJA1Xhl19+YefOnVFWryE+fPPNShyOJDp37k5eHuTkQG5ugdWhBSM7Gw4cgEOHlHjYrRBlpYTYvj2Djh07AkoMtEhoMdGCoBeXSy2NG6s+jA0bQmqqWhITweNRi9td8HfhxeWqvm4mtUZ86ipGfAyGY4NQSIlHdrYSkNxc5ZbKyVGiYV8OHFCuq0BAPQaD0VaHFg8h1A3e7VaL/jslRYlE48bqsUkT9ej1qn0SEyEpCRIS1HMtInWp7VVp4mPsc4PBUG+RUlkg2dlqyclRy+HDsG9fgRVy5IhacnLU/tpqCQQKXFHa0vB6lRAkJMDxx0Pz5gUC0qBBgWAkJhZdHKamjIURH4PBUCcJh5VYHDmixOTQIfW4fz/s3avE5eDBAjHRbi6/X1koxbmiUlKUgDRqpESlRQto1gySkwsEp2FDtZ9tLrChElRYfIQQZ6EmjG6IpEAbDAZDTPH7oy2SrCwlKDpmsnu3Wp+frxafTwmL3TLRbq7jjlOC0aCBcm01aaLWNWigRCQ5ueDRCEr8KFN8hBDfSin7Rv4eC/wR+BB4RAhxWuH6bAaDwVAaUqogvM7W0vGTPXuUqOzYocTE5yuwVqQssFS83gKXV7NmSkiaNlXWSmpqgWWSlKRcXXUpRnIsUR7Lxz7b6VbgfCnlPiHEdGAFajKnwWAwWOg5Ijo4v2cP7NwJu3apRYuLz6eExeVSgpKQoMTlpJOUoDRrphZtueglOdmISl2nPOLjiFSIdqCy4/YBSClzhBBFKgMYDIb6TyikXGFaXA4ehF9/hd9+U0KTlaWsG229uN0F2VvJydC+PbRsqSYwNm1aYLE0aqT2NdR/yiM+qahWBQKQQojjpZS7IxUIzG8Pg6GeEg4rEdm/Xy27dimX2K5dKpgfDhdkkOXnK5dYcrISmNRU6NwZTjgB2rRRLjIda6lAAQFDPaZM8ZFSti9hUxhV/LPK1JY+PwbDsYjPpyyXAweUqGzbBhkZSmT8fjXXJSdHiY0WmIQElTbctauyYtq0KcgOa9LECIyhbCqdah1pSbCtqgOw9fk5H1VjbZUQYp6U0t5qQff5+WdEVJ4CrivU58cL/E8IMR/IRvX5GSKl3BSponADqp6b7vNzmRCiS+TcQ6r6PgyG2oyUSkD27VNLZiZs3apcZYcPK4HRbjI9ATIpSYmIFpfWrZXlkpqq3GNNmyoxMhgqQ22Y51Mr+vxIKfdU1xs0GOJJIKAEZs8eJTK//KKEJitLlXnJy1Mxm+RktXi90LGjcpGdcIKyYFJTlbg0aWImRhqqh9ogPsX1+elXaB/d52cmRfv8PCKEeA5IQvX5+QFbnx8p5WqK7/OzrFCfHyM+hjqFlEpQdu1SmWRbtiiX2d69KgaTk6MsGY9HZYglJkKXLtCuHbRqpYSlcWPlKjNuMkO8Kc88n/6ozp/VVQSuvH1+ZgkhxqAKjlp9foQQfVB9fvZh6/MjhNB9fryoytX2Pj8zI31+1hPp81NkUELcikotL9JLxWCIN8FggTWzfTts2qTiMrpczNGjykJJSVEi07gx9OihXGY6s6xlS5NJZqg9lMfyuQF4UQixCdV87TMp5e4YjqHW9PkpdM7ZwGxQhUVj8k4NhnIQCqnJlpmZsHmzEprMTLX+yBGVYZaYqKwZt1vNiencWT22bKksmZQUMw/GULspT7bb7QCR4Pww4HUhRCoqrvIZsFxKGarCGFYBnYQQHVAWzdVAVKtNe58fVOzmtcj6wn1+eqCsHHvPHi/wZwoEyurzA9wCLI0IksEQd4LBggmYmzer+Myvvyp3WXa2it8kJhbUEuvRA3r2hBNPLIjJHAMdlw31kHLHfKSUPwE/oVxZiaj4yu+A54BiS2aX87hBIcSdwOcU9PnZaO/zAwwCnhJCWH1+Ii93A19F+rMfQaVgaxfafUIIe5+fRZH1XYE3hBBWn5/Kjt1gqAhSqvkyO3eqBICNG5ULLRxWIpOfX5AI0KABnHKKsma0RdOsmaoEYDDUB0w/n3Jg+vkYKoPPp2b8Z2bChg1KbPLylAhlZyuhadBAWTRuN3TooKyaHj1UQoBxmxnqOqafj8FQzYTDKstMZ51t2KBEB9RETVACo3/rnXIKnHmmEpxmzdS8GSM2hmMJIz4GQyUIBAoman7/Pfz8s3Kb6YoATqeyaoRQjyecoFKcTz4ZOnVS82gMhmMZIz4GQzkIh9V8mk2bYPVqJTa5ucp9BkpMXC712K+fsmx0TbPUVGPVGOoOUkr8IT+5gVyC4SACQdPkpjE/T2WayZ2DSo9+DFU3bZaUcmmsB2Yw1CRSqnk1mzdDerpaDhxQkzoTE5Wb7LjjoE8fVd9M1zZr2tRUBDDUDaSUBMIB8oP5hGSIcDhMWIYJyRBSSmSR6ZaxpTKWzzWoOmp3A4dRNdSM+BjqPDk5yo22fj18841KgT5yRNU4a9BAJQL06qXcZq1aqRYAxqIx1EbCMhy9hMNF1gXCgbiITElURny6AUellHsBhBBZsR2SwRAfwmGVjfbTT7BihXo8fFi5zxo1gjPOUGJz0knKfWaqAxhqC1JKguEggXCAYChIMBy0LJawDCORiELFYyojMtUpTpURn4dQ7RQ0n8doLAZDtXPkiMpGW7cOvv1WWTfZ2cqF1rmzitd06aISBIzYGGoKfdMPh8MEpRKXYDhIKBwiGA4SlmFLXEoSh/KIhnazhcIhQjJU5Dy+kA8pJa0bto7p+4PKic8EwCeEOBtYKKX8MMZjMhhihpSqVI12pWVkqDpoR46oMjRnnAF9+0L37kqADIZ4EgqHLAEIh8MEw0H8IT+BcACgVIEpj7hIKQkTjhaXUMRiCgcLhC2kzusP+5UlJQsErtZYPpE+OB6gOzBcCDFMSjk55iMjvk3mqmP8hpohFFKVA9atg6+/VlUF8vKUmy0pSbVuvuQSJTpJSTU9WsOxhJSSkAzhD/rJCeRY2WRQSEykBGG78Ueel3TMsAxHi4ltCcmQEpewXwlMyG9ts7vmdBzIH/JbYuQL+qrtWlQq1TpSF20NsEYIMTO2Q1LUQJM5Qx0mFFLtBL79FpYvVxM+s7JUmnNysspG69VLudU6dFDzcAyG6sIekwmE1KID/Toeo4WlsGWR8uR0HFlHODL1USU4UtLgLw8TatiAQ38Zr44bChQIS8RasmJAke2BcIBQOLrsZkiGCISUwARCAXxBH/nBfOt1AoEQonZluwkhUqSU2ZFmb2tRBUG3Ax2qaWzxbjJnqGOEw0pwVq1SggZ8Bw8AACAASURBVKPjNw0bqq6bV1+tJnW2bavWGQyxIizDlqjomIldXMqKyZQYpwmHkYcPk/zya+SH8tn56L20eGgqKa+8za6bR5GZ9Sv+cID8UH6UANnPFQ6HyQ/lKysn6McX8uEPKldeWIZxCDUXoDiBkUjiVXKtXOITKdD5ohDid8AYoCdwOjACZX1UB/FuMmeoA+gYzurVsGiRyk7z+VR16NatlSvttNNUOrSxbmqW/fv3c9999/Hoo4/WuZ5Y2pVlX7Q14wv6rBt+adZBadssN5nNOtJuseBfbqBDfhZtZr/JcbPfBGDT7y/i+7uugJw9lkD4Q35LhPKD+fiCPvwhf4HASAhH5YYVnLsiFM6aixXltXzOBbpKKXOFEAOllP8DVlbLiAqId5O56JObZnK1iqwsNdFz0SIVzxFCiYuUKlHg+uvVhE8zwTO+ZGZmMnXqVP75z3+Sn5+P1+vl4osvZsCAATzxxBMcOnSICy+8sNZ/h7TY5Ifyyc7PJhTpEmO/8RZnJZRGWKpAvz0OEwgpq0WLjrZMtAWlz7fnnlFc+Y8PrGN9Oe4i8o9kkh/Mt6wdbcEUJyYVERjtatOvEwg8Tg8JrgS8Li8Jrurp2VFe8WkspcyN/D0a+F+1jCaauDaZK4xpJlfz5OaqAp1Ll6qK0KASBBwO5XI76SQYMEDFclJSanas9REpJXv37qVhw4YEAgFWrFjB6tWrycjIYPv27WRkZLB161YArrnmGlq3bs3+/fv597//zbvvvkuvXr348ssv6d69e42+h8JuMXuGmf5bC0lha6ZUy0YH+qVKTdZxE52mHJQFWWR5wTwri80uaPYYTH4oH1/Ahz+Yz6CZ86LOddq0N/n8jxdGJR1URmC0q83pcOJxePC6vHhdXjxOD26nG6/Ti0M41D7CicvhqnHxaSGEWA98ADQTQjgijd2qk7g2mTPUDqRUVQYWLYKVK5U7TVszQqjEgqFD4dxzVTkbQ9XYv38/06dPZ9u2bVx55ZVceOGFeL1eli5dyp///GfWrl1b5DXNmjWjXbt29OzZkyuvvJKxY8fSvn17a/usWbNYt24daWlpeDyean8P2oVlFxIdDymc0VWWtVJsfCZShkYfzx/yW4KjBUyvt2eTgbrph2U4ykWmg/z+kL+oBSMlF7z4Gf0+WMmKK/rx+R8v5IIXP+OMD5SjqbAAFcaBAwSW683r8pLoSsTj9OBxeiyh0e9VIHA5XGqbUwmR2+m2hMflcOFxVM9nWF7x2YbKEjsNFczfL4T4NzBFSlmkBXUsqIEmc4YaJDcXvvsO5s9XnTy9XlVDLStLtRzo31+Vt2nf3jRUqwr5+fl88cUXrFy5kq1bt/LJJ5+QnZ1N06ZNmTt3btS+7dq14+mnnyYcDhMKhejbty/9+vWjYRnZG16vl759+1bn2wDUDfZo/lFyA7lVng8DBZaMziAL2AL72kKxx1aC0pYmLSVBaROgoCrMqYWotCyyKAtGCHwpCZbwIIR6BHwpCSCEdSwtMIWFQz93OpyWwNiFx+1w43Q4cTvc0X9HRMcpnDiEw3LFVRflaiYXcU8dllJ+G3l+CirAf5GUckC1jrAWYJrJVR+HD8PixfDZZ6olQYMGyrrJyVGTQG+4QVWINjXUKobf72fRokWsXLmSjRs3cvDgQUKhEOnp6Rw+fBiHw8EJJ5xA//79efDBB+nSpQuLFi1i9erVhMNhjj/+eH7/+9/j9Xpr9H3Y3WR2F5qOm1QGfRwdg/GHC9KUdQqyjs3YU5Ct8cgQ+cF88oJ55PnzLIEComInVcGBAJvAuIULj9triUxxAqMFRAuMy+my1mmh0daMy+GyLK7y0LJBy0q9j9KayZW7k6kQoquU8sdC69KklOmVGlUdwohP7Nm3D774QrnXfD6VPKA7eqalqeSBnj0hDl6bekNOTg4LFizg448/5uOPP+bQoUMIITjxxBNp0aIFTqeTE088kVGjRjFkyJC4uMTKg84i02VjCotLRdxmhY9rzXuxuc0s0YnEYnxBH8FwMKoemhYY+3wYX9AXJTJVFRi7NQRYMRdLXFzeIgJjucIibjQtNh6nB5dQ25yOiMtMuMptwRS2HAUi2m3nqtz/Skw6mRYWnsi6ei88htiyaxd8+il8+aWyepo0UUU8zzwTTj8dTjzRCE5l+Prrr7niiivYvXs3jRo1Yvjw4Vx11VWce+65JCcn1/TwgEhsJlK6RQuDP+gvMk+lMGXGaSLuMrvQ2GfxFxYZ+xwcHS+yLJlAHnnBPELh0rPJyv2jvZRAv84m8zq9uBwuSyS0G8zrVCKkLRiP06NiMdqiiVgyFbFgACVIEXHRx3I6lLvN6XBGZb9VJ8Z7bogLO3fCJ5+ozLW9e1VXz4svVhZOx44mjlNRAoEA2dnZbNmyhSVLljBp0iTatm3LggULGDhwIO4aqopauNqyPSBvtywqmrassYuMtkp0HbJgSBXCzAvkWZaT3ZKxJwT4AgVBfyjekqmIZWOfuClEQapygitBWSYu5Qqz9kPicqgbf4IrAa/TW+Aecxa4yrQ1U9EYjP06a8HyOpXgOR21YwKc+cobqpU9e+DjjwssnVat4NZb4aKLVNkbQ9ns2LGDw4cPk5WVxXvvvcc777zDvn37ovY577zzePfddzkujtVRpZTkBnKjLAqoejFMiBYZXSrGXkZGZ4/peS9a2PQkUH/IbyUH6GZpsbRkwjJsub8S3YkFIuP0WNaDfr96vwRnQhFXmcehLBudXVZeK0anhAtElNXiEA5r0RZUPKyYymDEx1AtHDgA//2vyl47dEi1KbjlFujd21SPLkwoFGLp0qV88sknbNu2jX379tG+fXvatGnD/PnzWbdunbWvx+Nh5MiRdO/enZSUFNq1a0fXrl3p0qVL3G4y2rrJ8mVZN/6o7RWsCWa3lrRo6CQAPffFF/JFZY2FwiGr8GVeIM+ydnSAvsyssjLQxwCirAY9L0bPh7G/Z6dw4na6SXQlqrhNxL2m3WPamnEKZ6U+K+0es2I8NlddXcSIjyGm5ObC55/Dv/+tyuC0aQN33QWDB5tYDkBGRgb//e9/+e2339i9ezfbt29n48aN7N69m4SEBDp06ECzZs1YsmQJmZmZ9OvXj+eee4527drh8Xg488wz42rdQEHg3hfwRbmzKioy1ox/e1xGFgT/fUGflZ5sFd+0tXrW23UhzFiJjC5D43a6LVeZXgpbI/aMMi0A9nkybmdBRplTVNy9Zb+uDuHAKZwkuBJI8iRVOLZT2zHiY4gJoZDql/P665CZqebm3HqrmhCamlrTo4s/R48e5bffVEGOQ4cOsXr1aubPn8/8+fORUuJ0OmnevDlt27bl3HPPZeTIkQwfPpwkW48Hv98f94w0nYasqzCHwiHygnlFbvJlCU8oHMIf9kelMuuUaS0kvqAPf9hvJQDoyZi5gVzLmollTMaeVZbgSiDBnWAF9b1Or2VFFHaV6fkzluXhiBaeioiC3R2ng/5WdlqhNOi6bNWUByM+hiohJWzaBC+/rCpMN2sGEyaoSaHHmuh89913vPzyy8yfP59ff/21yPY2bdrw4IMPMmbMGNq3b4+jjEJ01S08Wmj0jb60MjOloTPGtIVizZuJzJXJD+Zb1oomFA5Zloy2dnQaMxSTkFCJmIxTOPG4PCS6EklyJ5HoTizWktGuNJ1RZlk3kXiMS7iseEpFsMdltJjp7LK67jKLBXVWfKqj0ZyhYhw4AG+9Bf/7n0oeuOMOGDRIVSc4FsjPz2fhwoV8+umnfPnll/z8888kJiYyYsQIbr/9dtq1a4fD4SA5OZlevXrRpk2buN1wojpk2mqaFV5fksiUJDwlCY0/WJDOrGMv+le+3W2WnZ9tVYV2CEel58oUFhmvyxsV+Pe6vFZKsUSqiZpaaCLWToIzwZpbo4WhKq4yLTJ2S0ZbS8e60BRHnRSf6mg0J6U8Et93UXfx+2HhQiU8+fkwbBhccYWar1NfOXLkCD/88AO7du3ixx9/ZNWqVSxatIgjR46QnJzMwIEDGTduHKNHj6ZRHC+ElDJqxr6Ondhv/lBK/5gyyv5r68UX8lnCFQgFik1p1jXV9FyZXH9uqUkA5RUeu8tMZ5dpodEiYz++UzjxOpUYeZzK8tEpxtqqcYnyWx52V5l9jowWGbtLzlB+6qT4UD2N5qILWxmK5eef4cUXlYutY0cYN0491ld+/fVXXnrpJV588UWOHj1qre/UqRNXXnklV1xxBUOGDIlbGRoppbIg/NmEwqEoNxlEi0lFEgKs5mi6xEw4YM3VyQ/lF1RcDvuLJAJoodGxoarEZ+xC43V5SXYnK6FxJxSZJ6MtGW3BaFeZvcSMx+mpVEzGLihaZCqbpWYonroqPtXRaM5QCkePwpw5yuJJTITbb1fJBDU0l7HaCAQCfPbZZ8yZM4evvvqK3377DSEEo0aN4tprr6V169Z06NCBxo0bx21MVun+cJBsfzaBUKBK6c06oUBXANCuMrt7zO46k1LiD6uJmVpofEEfUHwiQLnqRRaaL5PgSoiyZjwOj6oMELE2dNkZ7SIrXMesonGU4rLKnA6nVXmgvmWW1UbqqvjEvNFckROYZnKA6puzdKkSnvx8GDkSRoyoX8kER48e5fnnn+fLL79k7dq1VpXnCy64gD59+jBs2DA6dy627VO1oLtUFhaayqQ3awvFHqOx2isHfVafGY1uqJYfzCcvkGelNkMVhEYI69upM82S3EkkuFU6s46zCCEsa0XHbuxFMa1Z/xWwQuwxHz2R0+10G2umFlBXxadaGs0Vev0x30xuzx6YNk21NejZE0aPrh8TRKWUrFq1im3btrFjxw6effZZ9uzZQ9++fRkzZgwXXHABF1xwQbWXqNEZYYWLapbUork8Nc7ssZ/8UMEM/9xgwfwYe2pzMBwkL5BHjj+HnEAO/pC/xBbM5RGa4uIzOtNMu8DsCQB6Po3ONtPiowP3dldbaRR2O+oJn/a4kKF2UVfFp1oazRkUoRDMnauy2JKSYPJkNVm0PvDZZ5/x6KOPsmLFCmvdGWecwbx586q1B422ZnSZmPxgfqWyzQofU2ee6Xkz2mrRYqOFTM+h0fvlBVRmmpQqS8suNJWNzyS5k0jyJJHoSrTKuujYjL0ts72njNdR8Uwzu9BogdM9aPRirJnaT50Un2psNHfMs2MHzJypXGxXXAFDhhR0Eq3L5OXlMX78eP7+97/Trl07XnrpJQYOHEijRo1o2bJlTG5W9mZk4XA4qjGZdl1VNlZjd5/ZS9D4g35yg2pSpm7RLKUkP5RvzaHJC+QVtAIoJDTlHYMWGiEEXqeXZI9KBCgcn3E5VCmaZHdygTWjKwJEkgHK4+oqPBmz8ETM8h7HUHspdz+fY5ljoZ9PMAizZ8NPPylr54EHoIyGlbWSzMxM5s+fz7Jly0hISKBRo0bs3r2br7/+ms2bN/OXv/yFxx57rMouNW3J2Nshl6ewZoWOHy4ojpkTyLFSnO3VmAEraSA3kEuuP5fcYK4VY6nMGPT8Gz1/JtmdTII7gURXIi6ny7LWrESBSC0znQygRUa7zIoTiMLXSLvhrOwyUTBfxghM3SUm/XwM9ZeMDPjrX1W1ghtvhF69anpE5WPhwoW88MILeDweHA6HFccBaNGiBVJKDh8+TIsWLejYsSN//etfufDCCyt8Hj1/RE/WDIaCHPUfLbauGFT8hh8Kh5RVYytFkx9SAf9sf7aVWQYFCQFaaLTFU9kJm/p1boebBHckEcBW10y/H51xpoUmxZMSVcusPLEZXfHZXvvMuMmOXYz4HMOEQvDGG7BmjZog+uCDkJBQ06Mqm5ycHCZNmsTMmTNp1aoVqamp+Hw+TjvtNMaNG8fQoUM55ZRTKnxD0wF7bdHo+maF59JA5a0a7ZqzWzXamtHxmkA4YMVSguEguYFccvw5VpzG3pxMUx7hsVs0ulhlklvFaHS3TIdwRNU0czlcalKnKzGqW2Z5rq1DOCw3XH0sjGmoGkZ8jlEyM2HGDAgE4LrroF/hWVK1kGAwyKuvvsrkyZPZvXs348aNY+rUqVHFOCuCnlhpiY2OmVSiGkBp59DC4g/7rYZnenKmDvqDEpq8YB6+gM/qkxMKh4qkOJc368w+h0ZnnBW2aHQFAN0qWWecWWJTjNAU13JZ95RxOVxWCwIz499QGkZ8jjFCIVV5et06Fdt55JHaH9uRUrJo0SLGjx/Pxo0bGTBgAB988AFnnnlmma8Nh8P4Qj4rpTksw5b1YZXtr2RVgOLGaU91ttc60xlmWuDsBT1z/Dml1jsrS2x0MoBDOEh0J5LsSY5ynRVObdbzZ7TAWLXNSqnQrAVHJxkYl5mhqhjxOYbYswemT1cCdN11qoV1bcXn8zFlyhTeeOMN9u7di8/no0OHDnzwwQdcdtll5brhBUIBDuQeKFVQqio2VsMzm/tMWy06KUBKadVC0+6zkuI0ZbnP9Gt08UxtzSS4VPmZ4uIzViKAbdJmeYTGOlZk1r+97YDBUFWM+BwDSKm6in7xRcG8ndrYwlpKydatW1m8eDEPP/wwu3btAqBhw4Zcc801vPTSSyQUE5TSFoe9YnMwpFxYsUSXpbE6aEYsm9xAblSBTX/Ib1UHsFcIACoUp7HPo0lwJZDsSbbiNE5REKOx4iruJBp4GpDoTrQyx1yi5N4wdqvPKZxWuX+nw4kDh1UJwGCoDoz41HNycuCJJ8Dng4suggsugNr443Xjxo2MGzeOxYsXF9l25MgR3nrrLfLy8nj77beBgliKnsGvXWhQNWtGY5/AqbPL8oJqvowWHu2+0/XOcvw5pRbXLI3C82jslQHcDjcIJRB6sqYVx3ElFvSeqUAigHa/6aoDxqIxxBsjPvWYX35RFag9HnjoodpZGmfbtm1Mnz6d2bNn06BBA6ZNm8bTTz/NgQMHovbLz89n4aKFHPYdxh/0F1uCpiqioysA6BYCOYEcy3rxhQpaFGirJtufTY4/p8R2zqXFaexZczrrzJ5RZrUFcHmtZAA9adNeHaAs7HXN7PXSjDVjqA3UO/GpSpO5+kI4DG+/rcrjnHACTJwIrlr0SWdkZPD++++zcOFCFixYgMPh4Oabb2bs2LGsXr26iPBo9u/bT16gwJVWlXRn7ULzBX3kBJXQWOnOoXzLqtHzaXICOQUlcQpZNeWN02jXWII7UsHZmQCCIq0Bkj3JlhC5HeWfya9FTc+lMWJjqM3UoltS1alKk7n4j7Z6yM6Gxx+H3Fy49loYOLD2uNlyc3N56qmnmDZtGn6/n5NOOomzzz6b6dOnc/rppzNx4kSmTp1Kk6ZNOLC/qAC1atOqUufVYpMfVBZNTiAnKlajG7LZs9Jy/KrIZkUqORfXJkC7zpLcSQihmpB5HB7L0knyJOFxeAoSAiI10cpD4RI0DbwNrPMYDLWdeiU+VK3JXJ1n0yZ45RUlNo8/Hr/OotOmTaNPnz4MHjzYWrd48WJWrVrF/fffz9GjR5k9ezbPPPMMe/bsYciQIbzyyissW76M635/Hf948x+069KOYZcPY036Gk7vezrPTX0Ov7+ghIzH42HiIxPLNR57YkBeII+jgaNWL5r8UL5VaFNbNfYqAoVrn5XmPis8l0bXO0t0JVpxFB2/SXAlkORS1QMqGqPR44KCKgFO4cThcFgp1CZuY6hr1DfxqXSTOSll8b6eOsK8efDNN9C2LfzhD/EtBtqnTx9GjhxJUlISe/fuJTU1ldzcXIYMGcL06dM5cOAA4XCY/v37s2fPHlxuF42Pb8yZ55/JLXfcwmWjLiMQDtCpSyf+OfefLF+6HI/HQ8PUhhzYf4AmTZvgy/PR4vgWRc5t76ipa6D5Qj58AVW9OSiDIMEfUgU4s/OzyQ3kEgwHi59TU4Irz54QkOhKjBYahyMqGcDj9JDkTiLZnWzNoymv68tuzegsNoejIKPNCIyhvlDfxKfSTeaKHKiONJMLhwuy2S65BPr3j/8YfvvtN/x+v9Vm+vDhwwAsX76cI0eOcNxxx/Hf//6Xfv368fLfX6b3mb3JDeSSmJjIY08/VuR46WvT+cc7/2DAOQOsdcuXLue7Nd/RZ0AfKykgL1jQ7MzeFC0UDllpztn+7AqnOjtwWPtZrQIi7jNdUTnBmWBZMUmugn41Hqen1BhNYVeZQziiimlqIdP9bAyG+kq9qmothOgPTJZSXhB5PhFASvlUCfunAD9JKUvtVlNbq1pnZ8PUSDrFH/4ArSoXEqkSo0aN4r333ivWPdWuXTvef/99unXrBi7I8mVVKElAl7/REzS1wOiYTViGQWLFanL8OVEZaOVJddZxGimlJTQ6RqMFQE/YTHSrrDMtOqUlA9gz8QTCcs3puTOmHYDhWOBYqmpd6SZzdY3Nm+HVV8HpjG9B0Jdeeom///3vPPTQQyxbtoy5c+eWuO+OHTvo2K0jWcEsCJadnRYKh6yCm9qq0a40f8ivytKEQ6p/jb9oxeeyrBp7irPb6bbERNc1czlcVuWAJHcSDT0NrXk2JSUCaJHRVoy2jMpqKWAwHOvUK/GpYpO5OsO8efD118rSGTeuerPZPv/8cx544AEWLFhAMBhk2rRpbN++nSuuuAKAq666ikWLFrF///4ir23StEmUONixF/W0qgXYLJxgOKjqskWy0rT7rKT6Z4WxZ57pbpfJ7mTLPeZ1eklwJ5DgVIkCKe4Uq4qz3Sopru+MUzitdGZdssZUbDYYKka9Eh8AKeWnwKeF1j1s+/t94P14jysWBIPw6KPqccQIGDCg7NcUZtq0aWRlZTFnzhx27NhB27ZtGT16NKmpqdx///388MMPjB07lmHDhrFz507mz5/P9u3badOmDT6fD6fTyc233Mz1N11P+47tSUpOYuyYsXw679Mi2Wk6ZqOrBfjDfktM8oJ5loWjWxdYac6RCZ4lJQUUfq5v/GEZtgpm6rkyui6ZXpLdySoRwFV0Do22YpzCaVkvWoycDqclaAaDoerUq5hPdVEbYj5798Lzz6u///QnOP74yh1n0qRJPPnkk0XWDx8+nBtuuIEZM2bw9ddfA6qmWteuXUlNTaVT504kN0jm4ssvpvPJnYu40P4999889ehT/Jb5G63atOL+h+5n2OXDyAnkWO4xHbMJh6MbouUGckusFFAYu0WjJ20meZKUq8ylHnVlgBR3CknuJKuNs0M4isRihBBWpprH6cHj8hgrxmCIEaXFfIz4lIOaFp8lS1RhUK9XtUCobAfoOXPmcP311xMOF3VZORwOwuEwrVq1onv37oy9dSxDhw1FCmkVzSwOexuBvEAeWflZVkKAP+RXS9BvxXDsFZ3LEpqowpqFumymuFOU2yzyvIGngVXZ2T5/pnDqsi5TY2IxBkP1cywlHNQr/H7lZvP5YNAg5WqLQsrogE/h5xHC4TDPP/88Dz30ULHCo/d58uknue2PtxGQyg2WE8wpUq8sJEPWvJqj/qOWVaPFxp4sYHedldZ5s3BlAGtipjtBWTCeFKuIpn7UCQIu4SriNtPZZVYLAWf5S9QYDIb4YMSnlrJhg2pxLQTcfTe0bl1oh8mT4fBh1Y5UCKY9/TQDPviAVZs2cfeRI1GxnHXr1vHWW2+Ver7jmhzHjXfcSH64YE6MvTFatj+bo/6jltDoBmg6RmNPBihJaOxpzQgKysxEGpSluFOsagBep1e1EHAn4XV6o7LNCicBuB1uK5HAxGYMhrqBEZ9aRl4ePPUU5OeragV33AFOhyRq/mw4rIRn5kz1fMYMBnzwAQNWrWIValbt9u3befLJJ2nevHmJhTrtXDvmWiv4r+M0uUFlvRzNP2oJjY7blNQ2QD+3u8zcTrdV50w3QNPxGW3N6CQAl1Bza+xxF23RuB1uy5rRc3BMfMZgqJuYmE85iFfM57PPYNkyZe3ccQe0ailZPnQoZGUxes8eKzttTosW0LAhZ3brhnjhBev1z1NQtE6TkJDAmLFjmPvWXA4eOMgjQCPbfg6H4JNTO5OfnEDe//3ZsmLyAnlWrKek+EyUJYOqO5bg9JLoScLj9CjLxZlIgkdZM7o9QII7ISrTLCoJIJIAoGf563kzFamDZjAYagcm5lNbicRoNm6Ed9+FcEhy5gDBRReh3GqHDsHhwwxYvZoJkZecsWMHZ+zYwUwhGLhwYVRdoMLCA6oPzp8e/hOtOrdiyn1P0Mjnt451v9fFx2knMmzlz6y4oh8LDvyCpPjJoFpkdGaYdpVpd1myJ5kus97Bm+1j/xMPkOROxu1wcdykxyG1EUcfuEcdx15SRijLJSobzQiMwXBMYMSnhlh+/vkE9uVy6bYryTpyhEapgv90/ATxbSoM+4KcnTtJfuUV3ELwHViCAfAdcEhKvnO51KSfCDMoKkCNWzTm882f40nzMOKeS5j89y9g7xEmABPyg7ByEyuu6McXdw6LCEzYaj5mT1tOciepwH+kMrNOcfY4I64y4aBB4EOSXv0nxyU0Jv+ZqSTe/wDuv71K+E/jSEpujjDiYjAYIhyz4lNW07lYETWpc/t22rZrx4jhd3DT0q308m9lMiHu4kfmZx3hjLWwITmZvv36smrVatYCfYtxi/ZCle9uHgxy4OS2vP/GX2g65ikm/KgKetsFqNuQbmT5sgA49bxTOfW8UzkigXMnW/us+/MY2rrVBMyG3oaqs6Yr0RIanQDgcrqsApg63hJlrcx6GdyJeGfOxPviy2rd+PE4IkkRBoPBoDkmxaecTediQlZWFk8++WQk1tKCu7bfwKwXczmJpvRiKxNYGWXVBHJyyM3P4dXUBvTKOlricZsDu046ntkvjYHcvbx/VmfC4RC+zAOQm0+j5o3oNbQXKQ1TOD7leBLdiTTwNKCBJ4VuT74Sdazf/2MNedOexOPyWj1iHMJhiUy5rRUhVPadToQAKxvPYDAY7ByT4kP5ms5VmczMTF544e/AHTRiLRNYCXzCXWxF9eeUgwAACd1JREFUcrjI/t+hrJoN3/8IwK6Ox9Nyy+4Sj//e6/fTMkFVWe42oRsN7kvhHG9DvvMkkeRKsrLLdNDeiQPPvX/G+dpc5J/+BDNmIO6+m+SZM0n2JFddKKSEuwo5/u66ywiQwWAowrEqPuVpOldlhBBkZ3cAmnMXKQBMYK1l6Wix0fyv0POWW3azpUUTnMEQ7Q8UFas/vLsFnpuBw2adlGmlHNcExo9HaEGYMUOtb9QoNsIzcyaMH6+Oq5+DESCDwRDFsSo+ZTadi0UzudatW9Oy9a/s2vkooGIxdhdbr0L7Tyj0nLQ0Op59Nvz1r5CWBunpqqjb7t2Qlobrry+Cw1WxG/vkydGVELQAVVUYhFACpoUnlsJmMBjqHceq+GQCJ9ietwF+s+8gpZwNzAY1z6eyJ7rxhhutQp4zitn+WnIyt+TksN7tplsgQHZKCilZWaqswcyZqqDbn/6kbuDnnKNu6HffDampMHBg5W7shfePlTBUl7AZDIZ6x7EqPmU2nYsVqampPDBxIh1nzeKmo0d5rUEDttx5J3/8179o9euv3HTzzdz0/PPqxj1vHimXXAIOR7TV8Mgj6gaub+z6hl5CLbcapbqEzWAw1CuO2QoHQoiLUEUBdNO5KSXtG5MKB4VqsSElTJgAjRurbaDK5jhs5WJqo7gYDAZDOTEtFapIzMrrlLMKtcFgMNQHShMfU5UxnhiXlMFgMADG8ikXQoh9wPZKvLQpsD/Gw4kltX18YMYYK8wYY4MZY8VoJ6VsVtwGIz7ViBBidUkmZ22gto8PzBhjhRljbDBjjB3G7WYwGAyGuGPEx2AwGAxxx4hP9TK7pgdQBrV9fGDGGCvMGGODGWOMMDEfg8FgMMQdY/kYDAaDIe4Y8akGhBAXCiF+FkJsFkL8pabHAyCEOEEIsVgI8aMQYqMQYnxk/XFCiAVCiF8ij41reJxOIcR3QohPIs87CCFWRsb3rhDCU5Pji4ypkRDifSHET5Hr2b82XUchxF2Rz3iDEOJtIURCbbiOQojXhBB7hRAbbOuKvW5C8ULkO/S9EOK0GhrfM5HP+XshxIdCiEa2bRMj4/tZCHFBdY+vpDHatt0rhJBCiKaR53G/hhXBiE+MsTWqGwacAlwjhDilZkcFQBC4R0rZFTgD+GNkXH8BFkopOwELI89rkvHAj7bnTwMzIuM7BNxcI6OKZibwmZSyC9ATNd5acR2FEK2BPwG9pZSnospHXU3tuI6vAxcWWlfSdRsGdIostwJ/q6HxLQBOlVL2ADYBEwEi352rgW6R17wU+e7XxBgRQpyAao65w7a6Jq5huTHiE3usRnVSSj+gG9XVKFLKXVLKtZG/j6JumK1RY/tnZLd/ApfWzAhBCNEGuBh4JfJcAOcC70d2qdHxAQghGgLnAK8CSCn9UsrD1KLriCoYnCiEcAFJwC5qwXWUUi4FDhZaXdJ1Gwm8IRUrgEZCiJbxHp+U8gspZTDydAWqAr4e3ztSynwp5TZgM+q7X62UcA1BFc2/n+jWMHG/hhXBiE/sKa5RXesaGkuxCCHao9oJrQRaSCl3gRIoVIfumuJ51BcoHHneBDhs+/LXhmt5IrAP+EfEPfiKECKZWnIdpZQ7gemoX8C7gCxgDbXvOmpKum618Xt0EzA/8netGZ8Q4hJgp5RyXaFNtWaMxWHEJ/aU2aiuJhFCpAAfABOklEdqejwaIcRwYK+Uco19dTG71vS1dAGnAX+TUvYCcqh5V6VFJGYyEugAtAKSUe6XwtT0dSyLWvXZCyEmoVzXc/SqYnaL+/iEEEnAJODh4jYXs67WfO5GfGJPmY3qagohhBslPHOklP+OrN6jTfHI494aGt4A4BIhRAbKVXkuyhJqFHEfQe24lplAppRyZeT5+ygxqi3X8Txgm5Ryn5QyAPwbOJPadx01JV23WvM9EkLcAAwHRsuCuSm1ZXwdUT801kW+O22AtUKI46k9YywWIz6xx2pUF8kouhqYV8Nj0vGTV4EfpZTP2TbNA26I/H0D8HG8xwYgpZwopWwjpWyPumaLpJSjgcXAlTU9Po2UcjfwqxDi5MiqIcAP1JLriHK3nSGESIp85np8teo62ijpus0Dro9kbJ0BZGn3XDwRQlwI/Bm4REqZa9s0D7haCOEVqillJ+DbeI9PSrleStlcStk+8t3JBE6L/J/WimtYIlJKs8R4AS5CZcZsASbV9HgiYzoLZXJ/D6RHlotQcZWFwC+Rx+NqwVgHAZ9E/j4R9aXeDLwHeGvB+NKA1ZFr+RHQuDZdR+BR4CdgA/Am4K0N1xF4GxWHCqBukjeXdN1QLqMXI9+h9ajsvZoY32ZU3ER/Z1627T8pMr6fgWE1dQ0Lbc8AmtbUNazIYiocGAwGgyHuGLebwWAwGOKOER+DwWAwxB0jPgaDwWCIO0Z8DAaDwRB3jPgYDAaDIe4Y8TEYYowQIiSESI9Ull4nhLhbCFGp75oQ4lN7JeVYU93HNxhKwqRaGwwxRgiRLaVMifzdHHgLWC6lfKRmR2Yw1B6M5WMwVCNSyr2ocvZ3RmaaOyM9YlZFeqzcBqq0jBBiacRi2iCEODuyPsPWn+WhSG+ZBUL16bk3sn6JEOJpIcS3QohN+rV2yjq+EOL2yLZ0IcQ2IcTiyPahQohvhBBrhRDvRWoDGgxVxoiPwVDNSCm3or5rzVGz5rOklH2APsDYSHmWa4HPpZRpqB5B6fZjCCF6A1egqpFfDvQudBqXlLIvMAEozsIq9fhSypcj2/qgZs4/FxG9B4HzpJSnoao63F25q2AwROMqexeDwRADdIXhoUAPIYSus5aKqgu2CngtUvz1IylleqHXnwV8LKXMAxBC/KfQdl0odg3Qvpjzl3V8zUxUXb3/RCqNnwIsV2Xi8ADflPlODYZyYMTHYKhmhBAnAiFUxWYBjJNSfl7Mfuegmum9KYR4Rkr5hn1zGafJjzyGKOZ7LaVcWsbxEUKMAdoBd9rOuUBKeU0Z5zYYKoxxuxkM1YgQohnwMjBLquyez4E7IhYIQojOQohkIUQ7VD+jv6Oqj59W6FDLgBFCiIRI3OXiCo6j1OMLIU4H7gV+L6XUzfxWAAOEECdF9kkSQnSuyHkNhpIwlo/BEHsShRDpgBvVgOxNQLexeAXlFlsbaXmwD9U6ehBwnxAiAGQD19sPKKVcJYSYB6wDtqPiL1kVGFOpx0dZO8cBiyMuttVSylsi1tDbQghvZL8HURXbDYYqYVKtDYY6ghAiRUqZHeleuRS4VUq5tqbHZTBUBmP5GAx1h9lCiFOABOCfRngMdRlj+RgMBoMh7piEA4PBYDDEHSM+BoPBYIg7RnwMBoPBEHeM+BgMBoMh7hjxMRgMBkPcMeJjMBgMhrjz/yVnD7hEplTmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This block of code was used for generating Figure 3. viewDesignSizesUpTo can be adjusted to zoom\n",
    "# in on lower x values, and evaluationThreshold can be adjusted to change which evaluation function\n",
    "# is used to plot the y position. For additional analysis, benchmark_designs_scores, benchmark_trajectories_scores\n",
    "# mhc1_design_scores, mhc2_design_scores, random_scores1, and random_scores2 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 a utility score array.\n",
    "# The ith entry (zero indexed) of the array contains the score evaluated by $\\mathcal{F}_{U_{i+1}}$.\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 a utility score array.\n",
    "\n",
    "# mhc1_design_scores, and mhc2_design_scores are multidimensional lists. The entry in [i][j][k] denotes\n",
    "# the score evaluated by the objective $\\mathcal{F}_{U_{j+1}}$ of a design of size k that was optimized\n",
    "# via the diminishing returns framework with the objective by $\\mathcal{F}_{U_{i+1}}$.\n",
    "\n",
    "# random_scores1 and random_scores2 are multidimensional lists. The entry in [i][j][k] denotes\n",
    "# the score evaluated by the objective $\\mathcal{F}_{U_{j+1}}$ of a design of size k. i indexes the random\n",
    "# samples.\n",
    "\n",
    "\n",
    "\n",
    "# threshold is off by 1 (set threshold to 2 if we want min(n,3))\n",
    "def sortOut(benchmark_designs_scores,\n",
    "            benchmark_designs, benchmark_trajectories_scores, mhccode, design_scores, traj, threshold, xmax):\n",
    "    fc = lambda x : -np.log( 1-(x/(threshold+1)) )/np.log(10)\n",
    "    n = xmax\n",
    "    benchmark_designs_points = {}\n",
    "    for k in benchmark_designs_scores:\n",
    "        if mhccode in k:\n",
    "            benchmark_designs_points[k.split(mhccode + '_')[1]] = (\n",
    "                len(benchmark_designs[k]), fc(benchmark_designs_scores[k][threshold]))\n",
    "    benchmark_trajectories_curves = {}\n",
    "    for k in benchmark_trajectories_scores:\n",
    "        if mhccode in k:\n",
    "            xs = []\n",
    "            ys = []\n",
    "            for i in range(1,n+1):\n",
    "                if i in benchmark_trajectories_scores[k]:\n",
    "                    xs.append(i)\n",
    "                    ys.append(benchmark_trajectories_scores[k][i][threshold])\n",
    "            benchmark_trajectories_curves[k.split('_' + mhccode)[0]] = (np.array(xs), fc(np.array(ys)))\n",
    "    \n",
    "    diminishingReturns = [ (np.arange(1,n+1,1), fc(np.array(design_scores[z][threshold][:n]))) for z in range(20)]\n",
    "    random = np.sort( np.array( [tr[threshold][:n] for tr in traj] ), axis = 0 )\n",
    "    random = fc(random)\n",
    "    q1 = random.shape[0]//4\n",
    "    q2 = random.shape[0]//10\n",
    "    q3 = random.shape[0]//20\n",
    "    \n",
    "    _, ax = plt.subplots(1,1, figsize = (6,2.5))\n",
    "    ax.plot(*benchmark_trajectories_curves[\"ILP_n3\"], linestyle = ':', marker = 'o', label = \"ILP3\", color = \"black\")\n",
    "    ax.plot(*benchmark_trajectories_curves[\"ILP_n5\"], linestyle = ':', marker = 'x', label = \"ILP5\", color = \"black\")\n",
    "    ax.plot(*benchmark_trajectories_curves[\"MarginalGreedy\"], label = \"MG\", color = \"black\")\n",
    "            \n",
    "    dr_distr = np.array([z[1] for z in diminishingReturns]).T\n",
    "    ax.fill_between(np.arange(1,n+1,1), np.max(dr_distr, axis = 1), np.min(dr_distr, axis = 1),\n",
    "                    color = \"blue\", label = \"DR\", alpha = 0.5, linewidth = 1, zorder = 10)\n",
    "    \n",
    "    pts = np.array([benchmark_designs_points[k] for k in benchmark_designs_points]).T\n",
    "    ax.scatter(*pts, marker = 'x', color = \"red\", zorder = 10, label = \"Other\")\n",
    "    \n",
    "    col = np.array([1,0.5,1])\n",
    "    ax.fill_between(np.arange(1,n+1,1), random[q1], random[-q1],\n",
    "                    zorder = -5, color = 1-col*0.5, label = \"Random(IQR)\")\n",
    "    ax.fill_between(np.arange(1,n+1,1), random[q2], random[-q2],\n",
    "                    zorder = -6, color = 1-col*(1/3), label = \"Random(IDR)\")\n",
    "    ax.fill_between(np.arange(1,n+1,1), random[q3], random[-q3],\n",
    "                    zorder = -7, color = 1-col*(0.2), label = \"Random(IVR)\")\n",
    "    ax.fill_between(np.arange(1,n+1,1), random[1], random[-1],\n",
    "                    zorder = -8, color = 1-col*(0.1), label = \"Random(full)\")\n",
    "    \n",
    "    ticks = np.arange(0,ax.get_ylim()[1],1)\n",
    "    ax.set_yticks(ticks)\n",
    "    ax.set_yticklabels([\"0\"] + [\"0.\" + ('9'*int(x)) for x in ticks[1:]])\n",
    "    ax.legend(ncol = 3, framealpha = 0, loc = \"upper left\")\n",
    "    ax.set_xlabel(\"Design size\")\n",
    "    ax.set_ylabel(\"$\\\\mathcal{F}_{U5}$ / 5\")\n",
    "    ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1]*1.5)\n",
    "    \n",
    "    plt.show()\n",
    "    \n",
    "# Which utility to evaluate against (1-20)\n",
    "evaluationThreshold = 5\n",
    "\n",
    "# How far to extend the x axis\n",
    "viewDesignSizesUpTo = 150\n",
    "\n",
    "sortOut(benchmark_designs_scores,\n",
    "        benchmark_designs,\n",
    "        benchmark_trajectories_scores,\n",
    "        \"MHC1\",\n",
    "        mhc1_design_scores,\n",
    "        random_scores1,\n",
    "        evaluationThreshold-1,\n",
    "        viewDesignSizesUpTo)\n",
    "sortOut(benchmark_designs_scores,\n",
    "        benchmark_designs,\n",
    "        benchmark_trajectories_scores,\n",
    "        \"MHC2\",\n",
    "        mhc2_design_scores,\n",
    "        random_scores2,\n",
    "        evaluationThreshold-1,\n",
    "        viewDesignSizesUpTo)"
   ]
  },
  {
   "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
}