{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sampling the potential energy surface "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction \n",
    "\n",
    "This interactive notebook demonstrates how to utilize the Potential Energy Surface (PES) samplers algorithm of qiskit chemistry to generate the dissociation profile of a molecule. We use the Born-Oppenhemier Potential Energy Surface (BOPES)and demonstrate how to exploit bootstrapping and extrapolation to reduce the total number of function evaluations in computing the PES using the Variational Quantum Eigensolver (VQE). \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import common packages\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from functools import partial\n",
    "\n",
    "# qiskit\n",
    "from qiskit.aqua import QuantumInstance\n",
    "from qiskit import BasicAer\n",
    "from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE, IQPE\n",
    "from qiskit.aqua.components.optimizers import SLSQP\n",
    "from qiskit.circuit.library import ExcitationPreserving\n",
    "from qiskit import BasicAer\n",
    "from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE\n",
    "from qiskit.aqua.components.optimizers import SLSQP\n",
    "\n",
    "# chemistry components\n",
    "from qiskit.chemistry.components.initial_states import HartreeFock\n",
    "from qiskit.chemistry.components.variational_forms import UCCSD\n",
    "from qiskit.chemistry.drivers import PySCFDriver, UnitsType\n",
    "from qiskit.chemistry.core import Hamiltonian, TransformationType, QubitMappingType\n",
    "from qiskit.aqua.algorithms import VQAlgorithm, VQE, MinimumEigensolver\n",
    "from qiskit.chemistry.transformations import FermionicTransformation\n",
    "from qiskit.chemistry.drivers import PySCFDriver\n",
    "from qiskit.chemistry.algorithms.ground_state_solvers import GroundStateEigensolver\n",
    "from qiskit.chemistry.algorithms.pes_samplers.bopes_sampler import BOPESSampler\n",
    "from qiskit.chemistry.drivers.molecule import Molecule\n",
    "from qiskit.chemistry.algorithms.pes_samplers.extrapolator import *\n",
    "\n",
    "import warnings\n",
    "warnings.simplefilter('ignore', np.RankWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we use the H2 molecule as a model sysem for testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "ft = FermionicTransformation()\n",
    "driver = PySCFDriver()\n",
    "solver = VQE(quantum_instance=\n",
    "             QuantumInstance(backend=BasicAer.get_backend('statevector_simulator')))\n",
    "me_gsc = GroundStateEigensolver(ft, solver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "stretch1 = partial(Molecule.absolute_stretching, atom_pair=(1, 0))\n",
    "mol = Molecule(geometry=[('H', [0., 0., 0.]),\n",
    "                        ('H', [0., 0., 0.3])],\n",
    "                       degrees_of_freedom=[stretch1],\n",
    "                       )\n",
    "\n",
    "# pass molecule to PSYCF driver\n",
    "driver = PySCFDriver(molecule=mol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.3])]\n"
     ]
    }
   ],
   "source": [
    "print(mol.geometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make a perturbation to the molecule along the absolute_stretching dof"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.8999999999999999])]\n",
      "[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.5])]\n"
     ]
    }
   ],
   "source": [
    "mol.perturbations = [0.2]\n",
    "print(mol.geometry)\n",
    "\n",
    "mol.perturbations = [0.6]\n",
    "print(mol.geometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate bond dissociation profile using BOPES Sampler\n",
    "\n",
    "Here, we pass the molecular information and the VQE to a built-in type called the BOPES Sampler. The BOPES Sampler allows the computation of the potential energy surface for a specified set of degrees of freedom/points of interest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First we compare no bootstrapping vs bootstrapping \n",
    "\n",
    "Bootstrapping the BOPES sampler involves utilizing the optimal variational parameters for a given degree of freedom, say r (ex. interatomic distance) as the initial point for VQE at a later degree of freedom, say r + $\\epsilon$. By default, if boostrapping is set to True, all previous optimal parameters are used as initial points for the next runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "stretch1 = partial(Molecule.absolute_stretching, atom_pair=(1, 0))\n",
    "mol = Molecule(geometry=[('H', [0., 0., 0.]),\n",
    "                        ('H', [0., 0., 0.3])],\n",
    "                       degrees_of_freedom=[stretch1],\n",
    "                       )\n",
    "\n",
    "# pass molecule to PSYCF driver\n",
    "driver = PySCFDriver(molecule=mol)\n",
    "\n",
    "# Specify degree of freedom (points of interest)\n",
    "points = np.linspace(0.1, 2, 20)\n",
    "results_full = {} # full dictionary of results for each condition\n",
    "results = {} # dictionary of (point,energy) results for each condition\n",
    "conditions = {False: 'no bootstrapping', True: 'bootstrapping'}\n",
    "\n",
    "\n",
    "for value, bootstrap in conditions.items():\n",
    "    # define instance to sampler\n",
    "    bs = BOPESSampler(\n",
    "        gss=me_gsc\n",
    "        ,bootstrap=value\n",
    "        ,num_bootstrap=None\n",
    "        ,extrapolator=None)\n",
    "    # execute\n",
    "    res = bs.sample(driver,points)\n",
    "    results_full[f'{bootstrap}'] =  res.raw_results\n",
    "    results[f'points_{bootstrap}'] = res.points\n",
    "    results[f'energies_{bootstrap}'] = res.energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare to classical eigensolver\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define numpy solver\n",
    "solver_numpy = NumPyMinimumEigensolver()\n",
    "me_gsc_numpy = GroundStateEigensolver(ft, solver_numpy)\n",
    "bs_classical = BOPESSampler(\n",
    "               gss=me_gsc_numpy\n",
    "               ,bootstrap=False\n",
    "               ,num_bootstrap=None\n",
    "               ,extrapolator=None)\n",
    "# execute\n",
    "res_np = bs_classical.sample(driver, points)\n",
    "results_full['np'] =  res_np.raw_results\n",
    "results['points_np'] = res_np.points\n",
    "results['energies_np'] = res_np.energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Energy')"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEWCAYAAABIVsEJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABA9klEQVR4nO3dd3gU5drH8e+92RRKhBRa6L0jSEARULCCiiio6AGleRR7OSjYu4J6FLvySrEAgiiCCiooiBTFgAFChxBaQEJC6Cm7+7x/7MCJIYFN29kk9+e69mJ26m8nS+48U54RYwxKKaVUQTnsDqCUUqp00gKilFKqULSAKKWUKhQtIEoppQpFC4hSSqlC0QKilFKqULSAqDJDRD4Ukads2O48ERlcyGXrichREQkq7lzFTUReFJEDIrIvd24RWSQit9udUfmX6H0gqjQQkSSgBuAC3MB64FNgvDHGY2O0ArE+x+3GmAV2ZykIEakHbALqG2P25zF9EfC5MeZjf2dT9tEWiCpN+hhjwoH6wBhgFDDB3khlg4g4zzJLPSA1r+Khyi8tIKrUMcYcMsbMAQYAg0WkDYCITBaRF63haBH5TkTSRSRNRH4TEYc1bZSI7BGRIyKySUQutcaHisg4EUm2XuNEJPTkdkWkr4jEi8hhEdkmIr2s8acO34hIYxH5RURSrcM9U0SkqjXtM7y/iL+1Dv88KiINRMSc/AUuIjEiMsfKvFVE/p1j+8+KyAwR+dTKvk5EYvPbT9Z67xeRRCvLazn2wRARWSoib4pIKvCsiFSx1p0iIjtE5EkRcYjIZcB8IMbKPTl37jy2PUxENojIQRH5UUTqF+ZnrQKbFhBVahljVgC7ge55TP6PNa0a3kNfjwNGRJoD9wKdrNbMlUCStcwTwAVAe+BcoDPwJICIdMZ7yOwRoCpwUY7lchLgFSAGaAnUBZ618t4K7MTbkqpsjHk1j+W/sHLHADcAL4vIJTmmX2vNUxWYA7ybxzpyuh6IBc4D+gLDckw7H0jEu39eAt4BqgCNgIuB24Ch1uG23kCylXvImTYoIn3x7u9+ePf/b8C0s+RUpZAWEFXaJQOReYzPBmrhPWafbYz5zXhP+LmBUKCViAQbY5KMMdusZQYCzxtj9htjUoDngFutacOBicaY+cYYjzFmjzFmY+6NGmO2WvNkWut4A+8v47MSkbpAV2CUMSbDGBMPfIz3F/lJS4wxc40xbuAzvIXuTMYaY9KMMTuBccAtOaYlG2PeMca4gCzgZuAxY8wRY0wS8N8cn78gRgCvGGM2WOt+GWivrZCyRwuIKu1qA2l5jH8N2Ar8ZB3CGQ3eX/DAg3hbBftF5AsRibGWiQF25FjHDmsceFsS2zgLEalhrXOPiBwGPgeiffwsMUCaMeZIrgy1c7zfl2P4OBB2lvMXu3KtKyafadFAMKd//pzb9lV94C3r8GE63p+PFHJdKoBpAVGlloh0wvtLaUnuadZf0f8xxjTCe9jn4ZPnOowxU40x3fD+ojPAWGuxZGvcSfWsceD9ZdvYh1gvW+tsa4w5BxiE95fnqWhnWDYZiBSR8FwZ9viw3fzUzbWu5Bzvc2Y5gLfVlvvzF2bbu4A7jTFVc7wqGGOWFWJdKoBpAVGljoicIyLX4D0X8LkxZm0e81wjIk1ERIBDeA9deUSkuYhcYp0czwBOACcvA54GPCki1UQkGngabwsCvFd7DRWRS60Ty7VFpEUe8cKBo8AhEamN95xJTn/jPcdwGmPMLmAZ8IqIhIlIO7yHzj7Pa34fPSIiEdbhsQeA6fls2w3MAF4SkXDrcNPDhdz2h8BjItIawDo5f2Ph4qtApgVElSbfisgRvH/hPoH3/MLQfOZtCizA+8t8OfC+MWYh3vMfY/D+xb0PqA48Zi3zIhAHrAHWAquscSdP2A8F3sRbkH7ln3+tn/Qc3hPWh4Dvga9zTX8Fb5FKF5GReSx/C9AAb0thFvBMEe8ZmQ2sBOKtPGe67Pk+4BjeE+tLgKnAxIJu0BgzC2+r7gvrMF4C3pPwqozRGwmVKqNExABNrfM+ShU7bYEopZQqFC0gSimlCkUPYSmllCoUbYEopZQqlLN1oFamREdHmwYNGtgdQymlSpWVK1ceMMZUyz2+XBWQBg0aEBcXZ3cMpZQqVURkR17j9RCWUkqpQtECopRSqlC0gCillCoULSBKKaUKRQuIUkqpQtECopRSqlC0gCillCoULSA+eOfLhxg79Xa7YyilVEDRAuKDlWlL+fXY73bHUEqpgKIFxAcRQVVJcYLH7bY7ilJKBQwtID6IDqtFhkNI3LPB7ihKKRUwtID4oFaVxgBsSFphcxKllAocWkB80LBmawB2HdhocxKllAocWkB80KZJFwD2HUmyN4hSSgWQctWde2FVi4ihittDqmu/3VGUUipgaAvER9HuIA5y2O4YSikVMLSA+CjSU5FUR5bdMZRSKmBoAfFRpDOCFCe4XNl2R1FKqYCgBcRH0RViyBZh047VdkdRSqmAoAXER7WrNgVg8059prpSSoEWEJ81rn0uADsP6N3oSikFWkB81rrx+QDsP7rL5iRKKRUY9D4QH1WpHEmUy0OqO8XuKEopFRC0BVIA0W4nBzlqdwyllAoIWkAKIIJKei+IUkpZtIAUQKQzigNOISPzuN1RlFLKdlpACqB6xTq4RVi/7U+7oyillO20gBRA7UjvvSBb9vxlcxKllLKfFpACaFL7PAD2pG22OYlSStlPL+MtgFaNO+H407D/2G67oyilVL48bjdbdyewacef7Ni/gb+P7uDePm9RI6p2sW5HC0gBVAyrRLTbkOZOtTuKUqo8M4bd+xJZt/13dvydwN4jSaRm/s1BzyFSHZnsd0KmQ/43v0DPbcupEXVDscbQAlJAUe5g0vReEKWUH2zduZa4jfNJ3L+alIy9HHQfJM2RQUqQh6NB/zwDUTnIQzUcVPdUpIU7gqjgWsRUbUyjWufSpsn5RFWpUez5tIAUUCThbAxKszuGUqoMSU3fy+9rf2Tjnj/YfXQre90H2BOcRXqOIhEaZKhuhChPKA2oQrSzJjWrNKRBjTa0bHA+dWs29HtuLSAFFBkcTaojnSPH0gmvVNXuOEqpUiQj8zhx639h7fbf2JG+gb3Ze9kbdIJ9TjDiPeQUFmSo43HQxhVF7eD6NKreno7NL6Vp3bY4goJs/gT/pAWkgKpXqgsntpKw7Xe6tOtldxylVIByubJZvGo2K7bMJfHYJvbJYfY4DVnWuQmHwxATBHXclekUFEPDyDa0a3Qx57W4CKcz2Ob0vtECUkB1o1vAroVs27NaC4hS6pQjx9L5ecUXrNr5C9szEtkWfIIj1iGoqCAPtV1h9HRXp27lZrSs24Xz21xBlcqRNqcuGi0gBdS07nmwC5IPbrU7ilLKRskpSfz85zQS9i1lu2s324Jdp1oXdYIMHdzRNAtvT/c2/WnfrGvAHX4qDrYUEBGJBKYDDYAk4CZjzME85hsLXG29fcEYM90aPxm4GDhkTRtijIkv0dCWFvU74FxqSDmxxx+bU0oFiE2JK/klfjqbUleSyH52BBs8IgQ5DI1wcLG7Nq2izueSDgNoVLe13XH9wq4WyGjgZ2PMGBEZbb0flXMGEbkaOA9oD4QCi0RknjHmsDXLI8aYmX7MDEBISCjVXZDm1iuxlCrLMjKOMXvxRyzf8T3r5W/2BlsnuZ2GxtkhXG3q0y7mIi7vfAtRVWvanNYedhWQvkAPa/gTYBG5CgjQClhsjHEBLhFZA/QCZvgpY74iPSEcFO2RV6myJv3IAb5e9A5/7v2Ztc6DHApyEOo0tMquQDdHc2IbXkGP2P5UDKtkd9SAYFcBqWGM2WsN7wPyusNlNfCMiPwXqAj0BNbnmP6SiDwN/AyMNsZk5rUhEbkDuAOgXr16xRI+UsJZE3SgWNallLLX7v1JzPr1Lf5KW0ZCyFFOOByEOz20cVUlNqon/XvcVyI34ZUFJVZARGQBkFe77omcb4wxRkRM7pmMMT+JSCdgGZACLAfc1uTH8BaeEGA83tbL83nlMMaMt+YhNjb2tO0URlRIddJJI+VgMtUiYopjlUopP9q0/S++Wf4eq4+sYkNoFi4RopwezndVp0u9q7ju4ru1leGDEisgxpjL8psmIn+LSC1jzF4RqQXsz2cdLwEvWctMBTZb40+2XjJFZBIwsljDn0X1yvXh6EbWJ/7JxR37+nPTSqlCWrF2AfP++pg1GRvYEuLGiBATZLjUXZeLmt7AVRfeVmruvwgUdh3CmgMMBsZY/87OPYOIBAFVjTGpItIOaAf8ZE07WXwEuA5I8FdwgHrVWsLRH0ncu4aL0QKiVKDatXcLk+Y/w4rMtewI8Y5rJNDHNOPyNrdxUYc+ZfLyWn+xq4CMAWaIyHBgB3ATgIjEAiOMMbcDwcBv3hrBYWCQdUIdYIqIVAMEiAdG+DN8i/qdYDskp2/z52aVUj6at/QzZq8bT1zwQTIdQjPjYICjHVd1uoPzWnS3O16ZYUsBMcakApfmMT4OuN0azsB7JVZey19SogHPokmd1oR6DKnZe88+s1LKLw4dTWPS3GdYnP4rW0INFYI9dHFV5/pz7+GSzsXbjbny0jvRC8ERFOS9F8Rz2r2PSik/i1u3iC+Wj+V3x04OBTmo4zD8yxnL4CueI6ZafbvjlWlaQAop0hNGmuOE3TGUKpeysjKZOv9VFuyezZrQDBxO6JBZmV51b+HGS+7T8xp+ogWkkCId55AY9LfdMZQqV7bv2cgnC55laXYC+4KFSKeHq00TBl30FK0bx9odr9zRAlJIUaE1OGJS2L0/iTrVG9gdR6kyy+N28/3ST/hu40TigtPJcggtPU76hV/J4Kue1vs1bKQFpJBqntMQDiWwfttyLSBKlZDvl0zms/VvsS7URcVgD11dNejX4T56xF5vdzSFFpBCq1+9NRz6lqT9688+s1KqQJbEz2XCH88RF3acqk4PAxzn8e9rxlAjqrbd0VQOWkAKqXXDzrAF9h1OtDuKUmXGum1xvPfzf1gWkkpoiKGPpxkP9H1XC0eA0gJSSLWrN6KSx0Nq1j67oyhV6u3au4W35t7PIscu3CHQI7sm9145jib12tgdTZ2BFpBCcgQFUc3lIM1z6OwzK6XydPBQCm9+cw8L3Os4GiRckFWFEReN0bvFSwktIEUQ6QnjgCPD7hhKlTrHM47x7qyHmHt8KalOB+2zKzC0w2i9Y7yU0QJSBJFBEawPSsbjduuNS0r5wON289GcJ5md+i17goXmbif3NbyT/pfcbXc0VQhaQIogKrQmGZ69JCVvLDfPQFaqsKb++F9m7PiEbaGGekZ4KPIGhlz1pP7xVYppASmCmKqNIe0vNuz4UwuIUvn4deVsPlz5HAmh2VQP8jC8wiXcfcvrhISE2h1NFZEWkCJoWLMNpM1kh94LotRpMjKP8/IXg/nObKCC03CTowP3D3ibKpUj7Y6miokWkCJo1egCWA9/H9lhdxSlAspvq+bwRtyTbA01dMqqxOO9JugluWWQFpAiqBFVmypuD6lu7VRRKfC2Ol75YijfedZRwWm4O/wa7ho81u5YqoRoASmiaHcQBz1H7I6hlO2WxM/ljRWj2RJqiM2qxGO9/o9m9dvZHUuVIC0gRRRpKpAcdMzuGErZJisrk1e+GMwcdwJhTsOI8Ku4a9AYvbqqHNACUkSRjkjincdwubJxOoPtjqOUXy1bPY///jGazaEeOmZV5PFe42lWv73dsZSfaAEpouiKMWRn72bLzjW0bNTR7jhK+UVWViZjvhjKHPcaQpyGOypdyT2DXtNWRznjsDtAaRdTtTEAm3bE2ZxEKf9YvuYH/jW5M1+atbTKqsCknp9z3w1vaPEoh7QFUkSNYtpDyjR2HthgdxSlSlRWViZjpw9ntiueEKfh3xUv595B/9XCUY5pASmiNo07w2rYf3SX3VGUKjEr1i7gteUj2Rjqpn1WGI9f/pEeslVaQIqqang0kS4Pqe79dkdRqkS8Mf0eph7/lWCnYXiFS7l/0Jva6lCAFpBiEe12cpCjdsdQqlhlZB5n9KfX8nPI37TJDubJS/+P1o1j7Y6lAogWkGIQSSWSHIftjqFUsdmRvJlR3w1gXaiLy7Jr8crgbwgLrWh3LBVg9CqsYhDhjCTF6f2LTanSblHcLP49tx+bQ7IZXqEHb97+kxYPlSctIMWgesU6uEVYv32l3VGUKpKP5zzNqLVPkiWG5xo9zIM3vWN3JBXAtIAUgzqRzQDYuusvm5MoVTguVzZPTLqetw7OIiY7iA96fkqfi4bZHUsFOD0HUgwa1z4X9sLutE12R1GqwFIOJvPojOuJCztOl8xzeO1f3+ozO5RPtAVSDFo3vgCHMaQc03tBVOmyauNvDJt5JatCj3GjtOXD4Yu1eCifaQukGFQMq0S025DqTrU7ilI++3LBu7y18wPcQfBIjdsY1HuU3ZFUKaMFpJhEuYM5iHbrrkqHMVOHMz3rD2p6hKc7v06Xdr3sjqRKIS0gxSSSymwKOmh3DKXO6NjxI4ya0odfQ1JpnxnG2H5fEVOtvt2xVCml50CKSWRwNVKDhCPH0u2OolSetu5MYMjn3fk1JJXernpMGLpUi4cqEi0gxaR6pToYEdYlrrA7ilKn+Wn5NO6cP4DtwS5GhPfi1eHfExISancsVcrZUkBE5EYRWSciHhHJt3MdEeklIptEZKuIjM4xvqGI/GGNny4iIf5Jnr86US0A2LYn3t4gSuXy8ZxneHzjSxjgpRaPc0+/1+yOpMoIu1ogCUA/YHF+M4hIEPAe0BtoBdwiIq2syWOBN40xTYCDwPCSjXt2zep5u7bec3CrzUmU+p/3vn6Ed9O+ol62gw8vm8aVXf5ldyRVhthSQIwxG4wxZ7vrrjOw1RiTaIzJAr4A+oqIAJcAM635PgGuK7GwPmpRvwNOY0g5vsfuKEoB3m7YPzo8j2ZZTj7oP49m9dvZHUmVMYF8DqQ2kPPOvN3WuCgg3RjjyjU+TyJyh4jEiUhcSkpKiYUNCQmlmgvS3Gkltg2lfDVmyjAmZSymTWYIHw2YT42ofP+LKFVoJXYZr4gsAGrmMekJY8zsktpubsaY8cB4gNjYWFOS24ryhOi9IMp2z382kC89a+iQGca7g37knEp6Z7kqGSVWQIwxlxVxFXuAujne17HGpQJVRcRptUJOjrddBOEkOA/YHUOVY09O6s9sx2Y6ZVTinVvnU6liuN2RVBkWyIew/gSaWldchQA3A3OMMQZYCNxgzTcY8FuL5kyiQqpzMMhBavo+u6OocsbjdvPoxGuY7dhMl8xzeH/IL1o8VImz6zLe60VkN9AF+F5EfrTGx4jIXACrdXEv8COwAZhhjFlnrWIU8LCIbMV7TmSCvz9DXmqEe2/KWq/3gig/8rjd/GdSb+YF7eCizCjeH7ZIHwCl/MKWrkyMMbOAWXmMTwauyvF+LjA3j/kS8V6lFVDqRreAoz+ybe9qunOt3XFUOeByZfPgpMv5NSSVS7Nq8MbwH3EEBdkdS5UTgXwIq9Rp2cB7T+Te9G02J1HlQUbmce6Z2JNfQ1Lp5arHG8O0eCj/8qmAiMh/RaR1SYcp7ZrUaUuox3DgxF67o6gy7njGMe6ZfAnLQg9xracpY4fM0eKh/M7XFsgGYLzVfcgIEalSkqFKK0dQENXckOZOtzuKKsOOHEtnxCc9WBF2jBukDS8N/VqLh7KFTwXEGPOxMaYrcBvQAFgjIlNFpGdJhiuNotyhpDmO2x1DlVEHD6Vw55RL+Cssg385O/LMbdPsjqTKMZ/PgVh9U7WwXgeA1XivhPqihLKVSpGOKqQEeeyOocqgv1P3MGL6FSSEZDEktBuPDZxsdyRVzvl6DuRNYBPeK6ReNsZ0NMaMNcb0ATqUZMDSJiq0OkeCHCSn7LA7iipDdu9P4u6vrmJTSDb/rnwF/7n5A7sjKeVzC2QNcK4x5k5jTO6bHALuclo71QxvCMC6bb/bnESVFTuSN3PP7D4khri5u+q13HfDG3ZHUgrw/T6Q1UBzb0e4pxwCdhhjDhV7qlKsfo02cPg7duxfd/aZlTqLv1P38MD3N7DbaXiw2s0MvvpJuyMpdYqvBeR94Dy8LREB2gDrgCoicpcx5qcSylfqtGzQCbbA3kOJdkdRpdyRY+k8+FUfkkI83B91oxYPFXB8PYSVDHQwxsQaYzriPe+RCFwOvFpS4UqjujUaU9Hj4UCG9oelCs/lyubBKb1JCM3mtgoXM6zPM3ZHUuo0vhaQZjn6ocIYsx5oYXUponJwBAVR3eXgoB7ZU4Xk7duqFytCj9Kf1jw84D27IymVJ18PYa0XkQ/wPhUQYIA1LhTILpFkpViEJ4w0R4bdMVQp9dSnN/JLyH4uz47h6aFT7I6jVL58bYEMBrYCD1qvRGAI3uKhNxPmEhlUlf1ODx632+4oqpQZO/V25ji20CXzHF4d8p3eYa4C2llbINYNhHONMT2B/+Yxy9FiT1XKRYfW4oRnHzv2baFh7RZ2x1GlxEezHmNK1u+cmxnGuME/4XQG2x1JqTM6awvEGOMGPNr/le9qVvHeC7Jhuz4XRPnmi5/e5KND39I0K4h3bp5HxbBKdkdS6qx8PQdyFFgrIvPhfw/9NsbcXyKpSrlGNdvCwa/1XhDlkx+XT+XNPR9T0yWM6zOTiCrV7I6klE98LSBfWy/lg1aNzocN8PcR7c5Endnva3/ixQ0vUcnAqz0nULdWU7sjKeUznwqIMeYTEakA1DPGbCrhTKVezei6nOP2cMD9t91RVADbtP0vnvzjIYwDnu/4Gm2anG93JKUKxNfOFPsA8cAP1vv2IjKnBHOVetVcQRz0HLY7hgpQu/cnMXLBbRwOgseaPUK39ledfSGlAoyvl/E+i7fTxHQAY0w80KhEEpURkVQgLSjT7hgqAKUfOcBD31zH7mDDA7UGc3W3IXZHUqpQfC0g2Xl0mqgPvTiDCEcE+4O8XVIodVJG5nEemNqbTSEubg/vxcBej9odSalC87WArBORfwFBItJURN4BlpVgrlKvWoUYshzClp1r7I6iAoTH7eahT65kVVgGA5yx3NP/dbsjKVUkvhaQ+4DWQCYwDTiM9450lY+YiCYAbN650uYkKlCMmnwtS0LTudrdkCcGTbY7jlJF5utVWMeBJ6yX8kGjWu0gZRo7UzbYHUUFgBc+G8QPzp1clBnFy8Nn2R1HqWLhUwERkWbASKBBzmWMMZeUTKzSr02TC2AN/H10p91RlM3e/WokMzyric2oyJtDf9T+rVSZ4euNhF8CHwIfA9pDoA+qhkcT6fKQ5k6xO4qy0ZzFE5h0+AdaZDsZN3AeISGhdkdSqtj4WkBcxpgPSjRJGRTtdpLGEbtjKJus2byM17e8QZQRXr36C6pUjrQ7klLFyteT6N+KyN0iUktEIk++SjRZGRBhKpLq0Mt4y6OUg8k88eudZDngqfNe1l6ZVZnkawtksPXvIznGGfRmwjOKDI5ipeMIGZnHCQutaHcc5ScuVzaPzLiOnaGGkdVvpft519odSakS4VMLxBjTMI+XFo+zqFahNi4RNiatsjuK8qPHP7mOlWEnuDHoPG69arTdcZQqMWcsICLyaI7hG3NNe7mkQpUVtSObAbB5pxaQ8uKtLx9gnnW57uP/mmR3HKVK1NlaIDfnGH4s17RexZylzGlapz0Au9M22xtE+cXsReP55NjPtMp08tpt3+vluqrMO1sBkXyG83qvcmnVsDNiDCnHd9kdRZWw+E1L+G/iW0S74NVrpusTBVW5cLaT6Caf4bzeq1wqVQwn2m1IdR+wO4oqQX+n7uHJxXeRHQQvx75K/ZhmdkdSyi/OVkDOFZHDeFsbFaxhrPdhJZqsjIh2B3Pwf08BVmWMy5XNIzOvY3eo4ZGag/W5HqpcOeMhLGNMkDHmHGNMuDHGaQ2ffB9c2I2KyI0isk5EPCISe4b5eonIJhHZKiKjc4yfLCLbRSTeerUvbJaSFkFlUoP0XpCy6rFP+vJXWAY3OmO1a3ZV7vh6I2FxSwD6AYvzm0FEgoD3gN5AK+AWEWmVY5ZHjDHtrVd8SYYtisjgaA4ECUeOpdsdRRWzcTPu4wfnLi7OiuaxWybYHUcpv7OlgBhjNvjwbPXOwFZjTKIxJgv4Auhb8umKV41KdTEibEiMszuKKkbfLPyIT48vpHWmk1dv/U6vuFLlkl0tEF/UBnJevrTbGnfSSyKyRkTeFJF8e6gTkTtEJE5E4lJS/N+xYe3I5gBsTf7L79tWJWPVxt/47/a3qeaC1/rM0CuuVLlVYgVERBaISEIer+JoRTwGtAA6AZHAqPxmNMaMN8bEGmNiq1WrVgybLpjm9c4DYHfaFr9vWxW/v1P38PRvd+MGnun0GnVrNbU7klK28bUvrAIzxlxWxFXsAermeF/HGocxZq81LlNEJuF9VklAatGgI85lhgPH99gdRRVRVlYmI2f2ZXeo4dFaQ7jw3N52R1LKVoF8COtPoKmINBSRELx3xc8BEJFa1r8CXIf3pHxACgkJpZoL0txpdkdRRfT4p9cRH5bJzcGd+deVj5x9AaXKOFsKiIhcLyK7gS7A9yLyozU+RkTmAhhjXMC9wI/ABmCGMWadtYopIrIWWAtEAy/6+zMURJQnmAOi94KUZm9Mv4cfg3fTM6saowdOtDuOUgGhxA5hnYkxZhZw2oOhjTHJwFU53s8F5uYxX6l6lG59Z13mORLZtXeLHjMvhWYt/JDPTvxKm6wQXh38nd1xlAoYgXwIq8zo2uQ6PCLMXvaR3VFUAa3Z8jtvbn+H6i549dqZ+lwXpXLQAuIHvS+8lapuD6tTl9kdRRXAoaNpPLPwDjIFnuo4hro19RE4SuWkBcQPnM5gWrsiWO88RFZWpt1xlA88bjejpvZla6jh31E30K3DNXZHUirgaAHxk3OrdedwkIO5SyfbHUX5YOwXw1kams41nsbcfu1zdsdRKiBpAfGTvt1G4DSGpYlz7I6izuLLBe8yPTuODhlhvHDrl3bHUSpg2XIVVnkUU60+zbOC2SA77Y6izmDN5mW8vfMDanmEV/t/jdNZ6E6nlSrztAXiR63CWrIjxPtLSgWe9CMHeHrRCLIEnoh9lZrRdc++kFLlmBYQP7qk3SAA5sZp19+BxuN2M2radWwLNdwZfZM+GEopH2gB8aML215JrWzDuiOr7Y6icnll2jCWhR6ij6cJw/o8Y3ccpUoFLSB+5AgKopWpyYaQDA4e8n/X8ipvM+a/zZeulZyXEcbzt86wO45SpYYWED/rVPcKMh3C7CUf2h1FAfGblvD2ro+o5YKxetJcqQLRAuJn13a/gwoeDyuTf7E7SrmXfuQAzyy+C5fAU51e05PmShWQFhA/C69UlVZZlVjn2I/H7bY7Trnlcbt5dFpftgcb7qx2sz7bQ6lC0AJig7ZVOpLidPDrqm/sjlJuvTx1CMtDD9PHNGPoNU/ZHUepUkkLiA2uOf8OABaun25zkvLpi5/eZKb7LzpmVOC5W/VnoFRhaQGxQfOGHWiUBeszN9kdpdxZtfE33t39MTHZMPaGWXrSXKki0AJik1bORmwJcbMjebPdUcqNg4dSePa3u3EDT53/OjWiatsdSalSTQuITbo1vR6PCHOWfmB3lHLB43bz6PS+JAUb7qxxM13a9bI7klKlnhYQm1x5wUAi3B5Wp/1ud5Ry4cWpt/F76BGuNc0ZcrWeNFeqOGgBsYnTGUwrVyTrgw+TkXnc7jhl2hc/vcFX7tV0zKjAs7d+YXccpcoMLSA26lCtG0eCHHyvD5kqMX+u+5m3d0+gtp40V6rYaQGxUR/rIVO/b//O7ihl0r4Du3hm2QMAPHfhW3rSXKlipgXERqceMuXZZXeUMsflyubRr/uRHAz31xlOp9aX2h1JqTJHC4jNWlVoxY4Qb6d+qvg8/un1/BWawYDgWG6+4mG74yhVJmkBsdll1kOmfoibaHOSsuPdr0YyL2gHF2VFMupmfXiXUiVFC4jNLmhzBTHZhoSj+pCp4vDDsilMPvwDLTKDeHXgdziCguyOpFSZpQXEZo6gIFpSi40hmfqQqSLavCOeVze8TBWP4ZVen1OpYrjdkZQq07SABIDOdXuR6RC++U0fMlVYR46l89iPgzniEEa3fpom9drYHUmpMs9pdwAFfbvfwVvTJ7Jy788MRe+SLiiP280jU/qwOdTD3ef04fILBtgdqczKzs5m9+7dZGRk2B1FlYCwsDDq1KlDcLBv90tpAQkAlSqG0zKrEuudKXjcbj1uX0AvTx3C0tB0rvE04q7rX7E7Tpm2e/duwsPDadCgASJidxxVjIwxpKamsnv3bho2bOjTMnoIK0C0qxpLitPBopWz7I5Sqkz54VW+tJ7t8cKtM+2OU+ZlZGQQFRWlxaMMEhGioqIK1LrUAhIgrrlAHzJVUL+v/Yn3kj+hbrbw2k3faDclfqLFo+wq6M9WD2EFiGb129M4U9iAPh/EF8kpO3juj4dxOOD5ru9SLSLG7khKlTvaAgkgLYMbsTnEzfY9G+2OEtCysjJ59Ov+7HPC/XXv5LxWF9sdSZUSQ4YMYebMoh/qjI+PZ+7cuWecJykpialTpxZ5W7748MMP+fTTT/2yrZy0gASQ7k37YUT4dpleznsmT3x2PavDMrk55AJuuvx+u+OocqioBcTlchVrnhEjRnDbbbcV6zp9oYewAsgVF9zCmE1jWXPwD7ujBKx3vnyIH5y76JEVzajBH9sdp1x77tt1rE8+XKzrbBVzDs/0aZ3v9KSkJHr37k23bt1YtmwZtWvXZvbs2VSoUIH4+HhGjBjB8ePHady4MRMnTiQiIuK0dSxYsIAxY8Zw+PBh3njjDa655hoyMjK46667iIuLw+l08sYbb9CzZ888x3ft2pWnn36aEydOsGTJEh577DFq1qzJAw94e34WERYvXszo0aPZsGED7du3Z/DgwURERPD1119z9OhR3G4333//PX379uXgwYNkZ2fz4osv0rdvX5KSkujVqxcdO3Zk1apVtG7dmk8//ZSKFSvSoEEDbrrpJubNm0eFChWYOnUqTZo04dlnn6Vy5cqMHDmSHj16cP7557Nw4ULS09OZMGEC3bt35/jx4wwZMoSEhASaN29OcnIy7733HrGxsYX+ednSAhGRG0VknYh4RCTf9CIyUUT2i0hCrvGRIjJfRLZY/57+LSmFnM5gWrkjWe/Uh0zlZe6ST5l8dD6tMp2MvVW7wC+vtmzZwj333MO6deuoWrUqX331FQC33XYbY8eOZc2aNbRt25bnnnsuz+WTkpJYsWIF33//PSNGjCAjI4P33nsPEWHt2rVMmzaNwYMH5zve4/Hw/PPPM2DAAOLj4xkwYACvv/467733HvHx8fz2229UqFCBMWPG0L17d+Lj43nooYcAWLVqFTNnzuTXX38lLCyMWbNmsWrVKhYuXMh//vMfjDEAbNq0ibvvvpsNGzZwzjnn8P7775/KX6VKFdauXcu9997Lgw8+mOdndLlcrFixgnHjxp3aD++//z4RERGsX7+eF154gZUrVxb5Z2FXCyQB6Ad8dJb5JgPvArkP7o0GfjbGjBGR0db7UcUd0g4dqndn6aFv+W7JJG649B674wSMTdv/4tVNY4kwMObqaVQMq2R3pHLvTC2FktSwYUPat28PQMeOHUlKSuLQoUOkp6dz8cXe82GDBw/mxhtvzHP5m266CYfDQdOmTWnUqBEbN25kyZIl3HfffQC0aNGC+vXrs3nz5nzH59a1a1cefvhhBg4cSL9+/ahTp06e27788suJjIwEvPddPP744yxevBiHw8GePXv4+++/Aahbty5du3YFYNCgQbz99tuMHDkSgFtuueXUvycLU279+vX7x/4BWLJkyalWUps2bWjXrl2eyxaELS0QY8wGY8wmH+ZbDKTlMakv8Ik1/AlwXfGls1ffbnd5HzKVpH9hn3ToaBqPzR/McQc8ce7zNKzdwu5IykahoaGnhoOCggp8PiH3parFcVny6NGj+fjjjzlx4gRdu3Zl48a8L4SpVOl/f/hMmTKFlJQUVq5cSXx8PDVq1Dh1D8aZMuY3nNPJfVSY/VMQpfUkeg1jzF5reB9QI78ZReQOEYkTkbiUlMDvrLBmdF1aZAWzXh8yBXivuHp46tVsCTXcHtmPnp362x1JBaAqVaoQERHBb7/9BsBnn312qjWS25dffonH42Hbtm0kJibSvHlzunfvzpQpUwDYvHkzO3fuPOP48PBwjhw5cmqd27Zto23btowaNYpOnTqxcePG0+bJ7dChQ1SvXp3g4GAWLlzIjh07Tk3buXMny5cvB2Dq1Kl069bt1LTp06ef+rdLly4+76OuXbsyY8YMANavX8/atWt9XjY/JXYIS0QWADXzmPSEMWZ2cW3HGGNExJxh+nhgPEBsbGy+8wWSVhVaM8OzmlXrfy3Xl6h63G7+80kvVoQe5QZpzR19X7A7kgpgn3zyyamT6I0aNWLSpEl5zlevXj06d+7M4cOH+fDDDwkLC+Puu+/mrrvuom3btjidTiZPnkxoaGi+43v27MmYMWNo3749jz32GEuWLGHhwoU4HA5at25N7969cTgcBAUFce655zJkyJDTTugPHDiQPn360LZtW2JjY2nR4n8t6+bNm/Pee+8xbNgwWrVqxV133XVq2sGDB2nXrh2hoaFMmzbN5/1z9913M3jwYFq1akWLFi1o3bo1VapUKeBezsUYY9sLWATEnmWeBkBCrnGbgFrWcC1gky/b69ixoykNlq/50bSZ3Ma8/Plgu6PY6tEJ15g2k9uYkR/3tjuKsqxfv97uCGXe9u3bTevWrfOcVr9+fZOSklKo9bpcLnPixAljjDFbt241DRo0MJmZmafNl9fPGIgzefxOLa2X8c4BBgNjrH+LrUUTCC5oewUxfxgSMorexCytXvhsEHODkrgoK5Kxw761O45Spd7x48fp2bMn2dnZGGN4//33CQkJKdI6bSkgInI98A5QDfheROKNMVeKSAzwsTHmKmu+aUAPIFpEdgPPGGMm4C0cM0RkOLADuMmOz1GSWhHD4pBkUtP3EVU1ryOBZdcb0+9lhmc1nTIq8ebQn7R3YlWuNGjQgISEhDynnbyiqjDCw8OJi4sr9PJ5sesqrFnGmDrGmFBjTA1jzJXW+OSTxcN6f4sxppYxJtiaf4I1PtUYc6kxpqkx5jJjTF5XapVq59frRZZDmP3b2a50Lls+nvM0n5xYRNuMYN4a9AMhIaFnX0gpZYvSehVWmden27+p6PGwcu9Cu6P4zfT543g/7WuaZAXx1k3fEV6pqt2RlFJnoAUkQFWqGE6rrMqsD/I+ZKqsm7f0M97Y/X/EZAvj+szU3nWVKgW0gASwtlVjOeB08POfZftBSUv++o6XNo3hHLfw2qWTqVurqd2RlFI+0AISwK654N8A/Lpxhs1JSs6azct4etUoHAZeOn8cLRt1tDuSCmBJSUm0adOmyOsZN24cx4+fub+5yZMnk5ycXORt+eLCCy/0y3aKmxaQANasfnuaZAobsrbaHaVEbN2ZwKOL7yBD4Nm2z9G57WV2R1LlRFELiLuYDysvW7asWNfnL6X1PpByo2VwY76TLSTuWkejuvZ0XlcSklN28J8fb+GAE55scD+XdL7B7kiqoOaNhn3FfK9SzbbQe8wZZ3G5XAwcOPAfXZ0vX76ckSNH4nK56NSpEx988AGhoaH8/PPPp43/6KOPSE5OpmfPnkRHR7NgwQKGDx9OXFwcIsKwYcOoW7cucXFxDBw4kAoVKrB8+XJatmzJgAEDmD9/Po8++ihHjhxh/PjxZGVl0aRJEz777DMqVqzIkCFDCAsLIy4u7h9dxk+ePJlZs2Zx6NAh9uzZw6BBg3jmmWcAqFy5MkePHmXRokU8++yzREdHk5CQQMeOHfn8888REebOncvDDz9MpUqV6Nq1K4mJiXz3nb195mkLJMB1b2Y9ZOr38XZHKTbpRw7w4Ky+7Aw2PFBjENf1vNPuSKoUyd3V+RtvvMGQIUOYPn06a9euxeVy8cEHH5CRkZHn+Pvvv5+YmBgWLlzIwoULiY+PZ8+ePSQkJLB27VqGDh3KDTfcQGxsLFOmTCE+Pp4KFSoAEBUVxapVq7j55pvp168ff/75J6tXr6Zly5ZMmDDhVMa8uowHWLFiBV999RVr1qzhyy+/zPO+jL/++otx48axfv16EhMTWbp0KRkZGdx5553MmzePlStXEij9+mkLJMBdfv7NRG4cQ3za73ZHKRbHM45x/7TebAxxcec5vbn1qtF2R1KFdZaWQknJ3dX5Cy+8QMOGDWnWrBng7cr9vffeo2fPnnmOz/0MjUaNGpGYmMh9993H1VdfzRVXXJHvtgcMGHBqOCEhgSeffJL09HSOHj3KlVdeeWpaXl3Gg7c796ioKMDb5fqSJUtOe6BT586dT3UH3759e5KSkqhcuTKNGjWiYcOGgLcr9/Hj7f+jUlsgAc7pDKajqUNc2HHGzbjP7jhF4nJl8+CnV/BXaAYDQy7gnn6v2R1JlUK5uzCvWrVqkdYXERHB6tWr6dGjBx9++CG33357vvPm7I59yJAhvPvuu6xdu5ZnnnnmVCsjr4wn3/vSlXxRu6v3Jy0gpcBzt0ynTWYwk48vZPL3pbNHWo/bzchJvVkeepi+phmj/qWPo1WFk7ur89jYWJKSkti61Xuxycmu3Js3b57neOAfXa0fOHAAj8dD//79efHFF1m1atVp8+TlyJEj1KpVi+zs7FNdvp+UV5fxAPPnzyctLY0TJ07wzTffnGpJnU3z5s1JTEw81ZXJyS7d7aYFpBQIr1SVcf2/pX628O7+6cxZPOHsCwWYpz69gZ9D/uby7Biev7XsXpasSt7Jrs5btmzJwYMHeeihh5g0aRI33ngjbdu2xeFwMGLECMLCwvIcD3DHHXfQq1cvevbsyZ49e+jRowft27dn0KBBvPLKK4C3hTFixAjat2/PiRMnTsvxwgsvcP7559O1a9d/dMUO/+syvnfv3qe6jAfv4an+/fvTrl07+vfv7/PzyCtUqMD7779/6lnp4eHhRe+KvRiIMaXiERnFIjY21hR3Z2L+tHnHGu6d/y+OOQxjO4ylW4dr7I50Vi5XNs9+fjOzZTNdM6vy/vBF2jliKbZhwwZatmxpd4yANmTIEK655hpuuOGfVxZOnjyZuLg43n333UKt9+jRo1SuXBljDPfccw9NmzbN95G2RZHXz1hEVhpjTqt22gIpRZrVb8dLXd4iCHh61SjWbQvsYrgjeTNDJ3Rhtmzm/IzKjBv8oxYPpQrp//7v/2jfvj2tW7fm0KFD3Hmn/VcvagukFPpx+VSe3vgS1VwOPugzm7o1G9kd6TRzFk/grS1vkBYk3BQcy6ibJ2jxKAO0BVL2aQukjLuyy794qPZw9gQbHp7Tj/QjB+yOdIrH7eaFzwbxTOKbiIGXm47isYGTtXgoVQZpASmlbr7iYf59zlVsCnHx4NSryMg8c7cM/rB7fxLDPu7CDM9q2mVWYNI139K76612x1JKlRAtIKXY3f1eZYAzlpVhJ/jPJ1fZ2u37vKWfMXTONawOPc5NjnOZdPvvAXloTSlVfLSAlHJPDJpMb1c9Foem8sSn/fy+fY/bzStThvLElrG4xfBco4d46tbP9ZCVUuWAFpAyYMyQOXTLrMp3jkTGTs3/Ltritu/ALm7/uCtTXXG0zgxjUu9ZXHvRcL9tXyllLy0gZYAjKIg3B/9I+4xQpmT9zkffPF7i25z/+3QGf9ObVaFH6U9rJg1fTv2YZiW+XaVU4NDOFMuIsNCKvH3zPG6ffhkfpc8hckEMN152b7Fvx+N28/qMu5iesYxwMTxT/z6u7zmi2LejAt/YFWPZmLaxWNfZIrIFozqPynd6UlISvXv3plu3bixbtozatWsze/Zsevfuzeuvv05sbCwHDhw41b3J5MmT+eabbzh27Bhbtmxh5MiRZGVl8dlnnxEaGsrcuXOJjIykR48enHvuufz666+4XC4mTpxIbGwszZs3Z9myZVSrVg2Px0OzZs1Yvnw51apVK9bPXVppC6QMiahSjXF9ZlIzG/676wMW/PFlsa4/5WAyd07ozmdZy2mRFcrEK2Zo8VB+t2XLFu655x7WrVtH1apV+eqrr844f0JCAl9//TV//vknTzzxBBUrVuSvv/6iS5cufPrpp6fmO378OPHx8bz//vsMGzYMh8PBoEGDTvVztWDBAs4991wtHjloC6SMqVurKWN6fMwDi4fzQsKzRFapyXktuhd5vYviZjE2/imSQ6Cvac6zw7/A6QwuhsSqtDpTS6EkNWzYkPbt2wPQsWPHUx0M5qdnz56Eh4ef6j+qT58+ALRt25Y1a9acmu+WW24B4KKLLuLw4cOkp6czbNgw+vbty4MPPsjEiRMZOnRoiXym0koLSBnUrukFPHv4ZR5b/TiPL7mL9ytN9/lphkeOpbMxaSXbk9exN30bB44lk5aVwh/BB6gkhidq38lNl99fwp9Aqfzl7u78xIkTOJ1OPB4PwD+6Vc89v8PhOPXe4XD8o6v0vLpar1u3LjVq1OCXX35hxYoVp/W6W95pASmjLu7Yl0cP7+PFpHcYOe8WPrrxB8IrVmXzzjVs37OGPWlbSTm6i4OZKaR7DpMuGaQ53RwMynVUU6CK00PrrIo80WsCzeq3s+cDKXUGDRo0YOXKlXTu3JmZM2cWah3Tp0+nZ8+eLFmyhCpVqpzq7fb2229n0KBB3HrrrQTp5en/oAWkDLuu552kfbuPt1K/5NpZl3PcIXhy/ZVV0ekhyiVEeEKJcUUQIVFEVYwhJqIR9Wu2okX9jkRU0WO+KrCNHDmSm266ifHjx3P11VcXah1hYWF06NCB7OxsJk6ceGr8tddey9ChQ/XwVR60M8Vy4OM5TxO3byFVnRFEVqhJzSoNqFe9Jc3qn0dMtfp2x1OlSFntTLFHjx6nruLKLS4ujoceeojffvvNhmT+V5DOFLUFUg7cfu3z+O/2QqXKjjFjxvDBBx/ouY98aAFRSpV7ixYtynP86NGjGT16tH/DlCJ6H4hSqkDK02Hv8qagP1stIEopn4WFhZGamqpFpAwyxpCamnrq+e2+0ENYSimf1alTh927d5OSkmJ3FFUCwsLCqFOnjs/zawFRSvksODiYhg0b2h1DBQg9hKWUUqpQtIAopZQqFC0gSimlCqVc3YkuIinADrtz5CMaOGB3iDPQfEWj+YpG8xVNUfPVN8ac1qdRuSoggUxE4vLqKiBQaL6i0XxFo/mKpqTy6SEspZRShaIFRCmlVKFoAQkc4+0OcBaar2g0X9FovqIpkXx6DkQppVShaAtEKaVUoWgBUUopVShaQPxARHqJyCYR2Soipz1cQEQeFpH1IrJGRH4Wkfo5prlFJN56zbEp3xARScmR4/Yc0waLyBbrNdimfG/myLZZRNJzTCvR/SciE0Vkv4gk5DNdRORtK/saETkvxzR/7Luz5Rto5VorIstE5Nwc05Ks8fEiUiKP8vQhXw8ROZTjZ/h0jmln/F74Kd8jObIlWN+3SGuaP/ZfXRFZaP3+WCciD+QxT8l9B40x+irBFxAEbAMaASHAaqBVrnl6AhWt4buA6TmmHQ2AfEOAd/NYNhJItP6NsIYj/J0v1/z3ARP9uP8uAs4DEvKZfhUwDxDgAuAPf+07H/NdeHK7QO+T+az3SUC0zfuvB/BdUb8XJZUv17x9gF/8vP9qAedZw+HA5jz+/5bYd1BbICWvM7DVGJNojMkCvgD65pzBGLPQGHPcevs74Ht/yn7IdwZXAvONMWnGmIPAfKCXzfluAaYVc4Z8GWMWA2lnmKUv8Knx+h2oKiK18M++O2s+Y8wya/vg/++eL/svP0X53vqsgPn8+t0DMMbsNcassoaPABuA2rlmK7HvoBaQklcb2JXj/W5O/wHnNBzvXwsnhYlInIj8LiLX2Zivv9X8nSkidQu4rD/yYR36awj8kmN0Se+/s8kvvz/2XUHl/u4Z4CcRWSkid9iUCaCLiKwWkXki0toaF1D7T0Qq4v3l+1WO0X7dfyLSAOgA/JFrUol9B/V5IAFERAYBscDFOUbXN8bsEZFGwC8istYYs83P0b4FphljMkXkTuAT4BI/Z/DFzcBMY4w7x7hA2H8BT0R64i0g3XKM7mbtu+rAfBHZaP1F7k+r8P4Mj4rIVcA3QFM/Z/BFH2CpMSZna8Vv+09EKuMtXg8aYw6XxDbyoi2QkrcHqJvjfR1r3D+IyGXAE8C1xpjMk+ONMXusfxOBRXj/wvBrPmNMao5MHwMdfV3WH/lyuJlchxD8sP/OJr/8/th3PhGRdnh/rn2NMaknx+fYd/uBWXgPG/mVMeawMeaoNTwXCBaRaAJo/1nO9N0r0f0nIsF4i8cUY8zXecxSct/BkjzBoy8D3lZeIt5DKydP9rXONU8HvCcEm+YaHwGEWsPRwBaK+UShj/lq5Ri+HvjdGo4Etls5I6zhSH/ns+Zrgfekpfhz/1nrbkD+J4Gv5p8nMFf4a9/5mK8esBW4MNf4SkB4juFlQC8b8tU8+TPF+wt4p7UvffpelHQ+a3oVvOdJKvl7/1n74lNg3BnmKbHvoB7CKmHGGJeI3Av8iPfKkYnGmHUi8jwQZ4yZA7wGVAa+FBGAncaYa4GWwEci4sHbWhxjjFlvQ777ReRawIX3P8oQa9k0EXkB+NNa3fPmn014f+UD71+AXxjrf4alxPefiEzDe6VQtIjsBp4Bgq3sHwJz8V4FsxU4Dgy1ppX4vvMx39NAFPC+9d1zGW+vrTWAWdY4JzDVGPODDfluAO4SERdwArjZ+hnn+b2wIR94/6j6yRhzLMeiftl/QFfgVmCtiMRb4x7H+4dBiX8HtSsTpZRShaLnQJRSShWKFhCllFKFogVEKaVUoWgBUUopVShaQJRSShWKFhBVqonIUR/medDqaqK4tnmdiLQqxvUtK8KyR61/Y0Rk5hnmqyoidxd2O0rlRQuIKg8eBApUQEQk6AyTrwOKrYAYYy4shnUkG2NuOMMsVQEtIKpYaQFRZYL13IhFVmePG0VkivUchPuBGGChiCy05r1CRJaLyCoR+dLqR+jk8xvGisgq4EYR+beI/Gl15PeViFQUkQuBa4HXrOc8NBaR9lZnjWtEZJaIRFjrWyTeZ5XEicgGEekkIl9bz154MUf2ozmGR4n3GRKrRWRMHp+zoZV9ba51NBDrmRUi0lpEVlj51ohIU2AM0Nga95qIVBbvs2dWWevqm2M9G0Tk/8T7fImfRKSCNa2JiCywsq0SkcbW+Ees/bRGRJ4r1h+sCmwlceu/vvTlrxfW8z7w3i18CG9/Pg5gOd7O7CDHcxnwdmmyGKvbCWAU8HSO+R7Nse6oHMMvAvdZw5OBG3JMWwNcbA0/j9WtBN6+t8Zaww8AyXif3xCKt+fTqFyfoTfeLi9OPhvmtG4lgDnAbdbwPTmWbYDV3QbwDjDQGg4BKpCrOw68d0efk2OfbMXb1UUDvD0OtLemzQAGWcN/ANdbw2F4W3VXAOOtZR3Ad8BFdn8v9OWfl3ZlosqSFcaY3QBWtw4NgCW55rkA7+GnpVY3EyF4i81J03MMt7H+yq+Kt6uZH3NvUESqAFWNMb9aoz4Bvswxy8muVtYC64wxe63lEvF2ZJeaY97LgEnGejaMybtbia5Af2v4M2BsHvMsB54QkTrA18aYLdZn/Ud04GURuQjw4O3Gu4Y1bbsxJt4aXgk0EJFwoLYxZpaVLcP6HFfgLSJ/WfNXxttbrr977VU20AKiypLMHMNu8v5+C96H6NySzzpy9mc0GbjOGLNaRIbgbeUUNpMnVz5PPvl8ccb+h4wxU0XkD7yd6M0Vbxf8iblmGwhUAzoaY7JFJAlvqyJnZvDuxwpn2JwArxhjPipAflVG6DkQVR4cwfu4T/A+da+riDQBEJFKItIsn+XCgb3i7S57YF7rM8YcAg6KSHdr2q3ArxTOfGDoySvGxHq2di5L8XYcSa5Mp4j32SeJxpi3gdlAO/65D8Dbg+x+q3j0BOqfKZjxPu1ut1gP5RKRUCvnj8CwHOeRaov3+ReqHNACosqD8cAPIrLQGJOCtzfhaSKyBu/hnhb5LPcU3uP+S4GNOcZ/ATwiIn9ZJ5IH4z2pvgZoj/c8SIEZb2+tc4A46xDcyDxmewC4R0TWkv/T424CEqx1tMH7ONNUvIftEkTkNWAKEGut57Zcny8/t+LtmXkN3nM1NY0xPwFTgeXWumbyz0KlyjDtjVcppVShaAtEKaVUoWgBUUopVShaQJRSShWKFhCllFKFogVEKaVUoWgBUUopVShaQJRSShXK/wPLqpE5CedvfwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "for value, bootstrap in conditions.items():\n",
    "    plt.plot(results[f'points_{bootstrap}'], results[f'energies_{bootstrap}'], label = f'{bootstrap}')\n",
    "plt.plot(results['points_np'], results['energies_np'], label = 'numpy')\n",
    "plt.legend()\n",
    "plt.title('Dissociation profile')\n",
    "plt.xlabel('Interatomic distance')\n",
    "plt.ylabel('Energy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare number of evaluations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no bootstrapping\n",
      "Total evaluations for no bootstrapping:\n",
      "1714\n",
      "bootstrapping\n",
      "Total evaluations for bootstrapping:\n",
      "1266\n",
      "np\n",
      "Total evaluations for np:\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "for condition, result_full in results_full.items():\n",
    "    print(condition)\n",
    "    print('Total evaluations for ' + condition + ':')\n",
    "    sum = 0\n",
    "    for key in result_full:\n",
    "        if condition is not 'np':\n",
    "            sum += result_full[key]['raw_result']['cost_function_evals']\n",
    "        else:\n",
    "            sum = 0\n",
    "    print(sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extrapolation \n",
    "\n",
    "Here, an extrapolator added that will try to fit each (param,point) set to some specified function and suggest an initial parameter set for the next point (degree of freedom). \n",
    "\n",
    "- Extrapolator is based on an external extrapolator that sets the 'window', that is, the number of previous points to use for extrapolation, while the internal extrapolator proceeds with the actual extrapolation.\n",
    "- In practice, the user sets the window by specifies an integer value to num_bootstrap - which also the number of previous points to use for bootstrapping. Additionally, the external extrapolator defines the space within how to extrapolate - here PCA, clustering and the standard window approach. \n",
    "\n",
    "In practice, if no extrapolator is defined and bootstrapping is True, then all previous points will be bootstrapped. If an extrapolator list is defined and no points are specified for bootstrapping, then the extrapolation will be done based on all previous points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Window Extrapolator: An extrapolator which wraps another extrapolator, limiting the internal extrapolator's ground truth parameter set to a fixed window size\n",
    "2. PCA Extrapolator: A wrapper extrapolator which reduces the points' dimensionality with PCA, performs extrapolation in the transformed pca space, and untransforms the results before returning.\n",
    "3. Sieve Extrapolator: A wrapper extrapolator which performs an extrapolation, then clusters the extrapolated parameter values into two large and small clusters, and sets the small clusters' parameters to zero.\n",
    "4. Polynomial Extrapolator: An extrapolator based on fitting each parameter to a polynomial function of a user-specified degree.\n",
    "5. Differential Extrapolator: An extrapolator based on treating each param set as a point in space, and performing regression to predict the param set for the next point. A user-specified degree also adds derivatives to the values in the point vectors which serve as features in the training data for the linear regression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Here we test two different extrapolation techniques"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define different extrapolators\n",
    "degree = 1\n",
    "extrap_poly = Extrapolator.factory(\"poly\", degree = degree)\n",
    "extrap_diff = Extrapolator.factory(\"diff_model\", degree = degree)\n",
    "extrapolators = {'poly': extrap_poly, 'diff_model': extrap_diff}\n",
    "\n",
    "for key in extrapolators:\n",
    "    extrap_internal = extrapolators[key]\n",
    "    extrap = Extrapolator.factory(\"window\", extrapolator = extrap_internal)\n",
    "    # define extrapolator\n",
    "    # BOPES sampler\n",
    "    bs_extr = BOPESSampler(\n",
    "        gss=me_gsc\n",
    "        ,bootstrap=True\n",
    "        ,num_bootstrap=None\n",
    "        ,extrapolator=extrap)\n",
    "    # execute\n",
    "    res = bs_extr.sample(driver, points)\n",
    "    \n",
    "    results_full[f'{key}']= res.raw_results\n",
    "    results[f'points_{key}'] = res.points\n",
    "    results[f'energies_{key}'] = res.energies\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "poly\n",
      "diff_model\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Energy')"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEWCAYAAABIVsEJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABFjklEQVR4nO3dd3hUZfr/8fc9ySQzoYXeMaB0AkGaCNJElLXgqqAuKJF1FXtZXAVU7KLrqj8VRVYFQUHFiopfBQWpigFD7xg6EkoSUiaZyTy/P+aQDSGB9DNJ7td1zcXMOWfO+eRkyD3PKc8jxhiUUkqponLYHUAppVTFpAVEKaVUsWgBUUopVSxaQJRSShWLFhCllFLFogVEKaVUsWgBUZWGiEwVkcds2O53IjK6mO9tISKpIhJS2rlKm4g8IyJHRORQ3twislhEbrU7oypfoveBqIpARBKAhoAPyAY2ATOBacYYv43RisT6OW41xiy0O0tRiEgLYCtwjjHmcD7zFwMfGGPeKe9syj7aAlEVyZXGmBrAOcBk4GHgXXsjVQ4iEnqWRVoAR/MrHqrq0gKiKhxjTLIxZh5wPTBaRDoBiMgMEXnGel5PRL4RkSQROSYiS0XEYc17WET2i8gJEdkqIhdb08NF5FUROWA9XhWR8JPbFZFhIhIvIikislNELrOm5xy+EZFzReQnETlqHe75UEQirXmzCPwh/to6/PMvEYkSEXPyD7iINBGReVbmHSLyj1zbf0JEPhGRmVb2jSLSvaD9ZK33XhHZZWX5d659ECsiy0XkFRE5CjwhIrWsdSeKyG4ReVREHCIyGFgANLFyz8ibO59tjxGRzSJyXES+F5FzivO7VsFNC4iqsIwxq4B9wEX5zP6nNa8+gUNfEwAjIm2Bu4EeVmvmUiDBes9E4AIgBugC9AQeBRCRngQOmT0ERAL9cr0vNwGeB5oA7YHmwBNW3puAPQRaUtWNMS/m8/6PrNxNgOuA50RkUK75V1nLRALzgDfyWUdufwW6A+cDw4Axueb1AnYR2D/PAq8DtYBWQH/gZuAW63DbUOCAlTv2TBsUkWEE9vc1BPb/UmDOWXKqCkgLiKroDgB18pnuBRoTOGbvNcYsNYETftlAONBBRJzGmARjzE7rPSOBp4wxh40xicCTwE3WvL8D7xljFhhj/MaY/caYLXk3aozZYS2Taa3jZQJ/jM9KRJoDfYCHjTEeY0w88A6BP+QnLTPGzDfGZAOzCBS6M3nBGHPMGLMHeBW4Mde8A8aY140xPiALuAEYb4w5YYxJAP6T6+cvirHA88aYzda6nwNitBVS+WgBURVdU+BYPtP/DewAfrAO4TwCgT/wwP0EWgWHReQjEWlivacJsDvXOnZb0yDQktjJWYhIQ2ud+0UkBfgAqFfIn6UJcMwYcyJPhqa5Xh/K9TwdcJ3l/MXePOtqUsC8eoCT03/+3NsurHOA/2cdPkwi8PuRYq5LBTEtIKrCEpEeBP4oLcs7z/oW/U9jTCsCh30ePHmuwxgz2xjTl8AfOgO8YL3tgDXtpBbWNAj8sT23ELGes9YZbYypCYwi8MczJ9oZ3nsAqCMiNfJk2F+I7RakeZ51Hcj1OneWIwRabXl//uJsey9wuzEmMtfDbYxZUYx1qSCmBURVOCJSU0SuIHAu4ANjzPp8lrlCRM4TEQGSCRy68otIWxEZZJ0c9wAZwMnLgOcAj4pIfRGpBzxOoAUBgau9bhGRi60Ty01FpF0+8WoAqUCyiDQlcM4ktz8JnGM4jTFmL7ACeF5EXCLSmcChsw/yW76QHhKR2tbhsfuAjwvYdjbwCfCsiNSwDjc9WMxtTwXGi0hHAOvk/PDixVfBTAuIqki+FpETBL7hTiRwfuGWApZtDSwk8Md8JfCmMWYRgfMfkwl84z4ENADGW+95BogD1gHrgTXWtJMn7G8BXiFQkH7m1G/rJz1J4IR1MvAt8Hme+c8TKFJJIjIun/ffCEQRaCl8AUwq4T0jXwGrgXgrz5kue74HSCNwYn0ZMBt4r6gbNMZ8QaBV95F1GG8DgZPwqpLRGwmVqqRExACtrfM+SpU6bYEopZQqFi0gSimlikUPYSmllCoWbYEopZQqlrN1oFap1KtXz0RFRdkdQymlKpTVq1cfMcbUzzu9ShWQqKgo4uLi7I6hlFIViojszm+6HsJSSilVLFpAlFJKFYsWEKWUUsWiBUQppVSxaAFRSilVLFpAlFJKFYsWEKWUUsWiBaQQZo2byIx7xp99QaWUqkKq1I2ExZV1uD7GcdpNmEopVaVpC6QQJDQZr7MOPq/X7ihKKRU0tIAUgrgz8YeEsyXuV7ujKKVU0NACUgiueuEA7NJ+tJRSKocWkEKo26opACn7jticRCmlgocWkELo1K8/AFnJOviWUkqdpAWkEJq0bE2o9wTGE2F3FKWUChp6GW8hhfiOgYm0O4ZSSgUNbYEUkphj+EPq2B1DKaWChhaQQhJnCl5nHbI8HrujKKVUUNACUkiOCC/G4WTDiiV2R1FKqaCgBaSQ3PXdAOxeu87mJEopFRy0gBRSw3atADix/7jNSZRSKjhoASmk6H4DAPAli71BlFIqSGgBKaS6jZrizErCZFWzO4pSSgUFLSBFEJJ9FExtu2MopVRQ0AJSJEl6L4hSSlm0gBSBOE+QFVabtJRku6MopZTttIAUQUh1H0gI637+ye4oSillOy0gRRDRsDoA+zdtsTmJUkrZTwtIETTp0BaAtEMnbE6ilFL20954i6Bz/4uJnx9H9gndbUqp4OXzetm0ajl/rPmdlD1H8SbDJQ/dRPPWbUt1O/qXsAhqREbizEoCf3W7oyilqjDj97Nr/Vq2/LqSpISDZB334/e4wV8Lv9TF56yDPyQM6BJ4gwM2LVuiBcRuIf6jGPReEKVU2TLGsPGXZWxevIzUfSfwZ7gguxZGauMLrUt2qBtoZz0gJCSNUP9RHP4/Cc/ehoRn4qrnov55LejUrz+NWgwq9YxaQIpI5Dg+Rxu7YyilKpGDu3cQ9818jm3/k+xkF8ZfH19oE3zO6kAvABwhmTj9RxD/ccKyExBnBmF1nNRp2Zg2F/Tm3E6lXyDORgtIUTnT8IZEcjzxELXrN7I7jVKqAklLSeLXb+dxIH47WUcF461LtqMx3vB6QCegE47QTMKyDhDq20hY+AmqN69Gm7696XjBYEKdTrt/hFNoASmi0Jp+SIN1P/9E/+v+ZnccpVSQyvR4WDxnFgfX/EF2SiRGGpEV1hDjaAY0g5Bswn2HCcneS2j2WlwNQ2nerSO9hl5FmMtld/xC0QJSRNWbRnJiGxzassvuKEqpIHI88RA/f/Ahx7Yew6Q3whvWkuzQc4FzcYYkEeLbj8u3k9A6Xuq3b07PK66ibsNL7I5dIlpAiqh5dEcOboOMw2l2R1FK2Shhy3p++eRL0vd4Md5mZIafg3F0BSDccZAw31pCI9NoPaAnPS69MugOP5UGWwqIiNQBPgaigARghDHmtJGaROQF4HLr5dPGmI+t6TOA/sDJTqlijTHxZRra0unC/vw2dzn+tLDy2JxSKkisXbKQ9d/+TObhMIy/BZmupiB9EEc2YezF5VuJq4mh61VDaN9zpN1xy4VdLZBHgB+NMZNF5BHr9cO5FxCRy4HzgRggHFgsIt8ZY1KsRR4yxnxajpkBcFerhtN7DOOvWd6bVkqVo7TkJOa/9SYntvrIpi1Z4fWBgYGT3JkJuP0/U6NVBBf9bQSNzqnYh6KKy64CMgwYYD1/H1hMngICdACWGGN8gE9E1gGXAZ+UU8YCObKPYUTvBVGqskk8sJeFU98hY48bb0gHfM4LcIRmEZa5FbdjLfU7NWbA326iRuRQu6MGBbsKSENjzEHr+SGgYT7LrAUmich/gAhgILAp1/xnReRx4EfgEWNMZn4bEpHbgNsAWrRoUSrhxXEcX0inUlmXUspeuzauY/n0j8k6XJes8Pb4Q/oTEppOmHcT1RulcvHtsTRqcZndMYNSmRUQEVkI5HejxMTcL4wxRkRM3oWMMT+ISA9gBZAIrASyrdnjCRSeMGAagdbLU/nlMMZMs5ahe/fup22nOCQ8HZ/U4MAf22nSsnVprFIpVY7WLv2J3z9ZQHZyMzJdbTCOi3GGJhHuXU31KGHoXXdRI/IKu2MGvTIrIMaYwQXNE5E/RaSxMeagiDQGDhewjmeBZ633zAa2WdNPtl4yRWQ6MK5Uw59FaC2BFNi0fLkWEKUqiKWff8KOBfH401vhcUWBXEJYyGFc2cuoGx3J0H/cQZjrGrtjVih2HcKaB4wGJlv/fpV3AREJASKNMUdFpDPQGfjBmney+AhwNbChvIID1Gxel5SNkLh9d3luVilVRDvWxrF06udkZ7QPXDXFYMLZg9u/mOYXtWLgjTdVystry4tdBWQy8ImI/B3YDYwAEJHuwFhjzK2AE1gaqBGkAKOsE+oAH4pIfUCAeGBseYZv2e189m3MxnM039MuSimbfTP1Df5ceYKssK74QwbjMruIcCyk3eU96X15rN3xKg1bCogx5ihwcT7T44BbreceAldi5ff+8u81LJcOPS5k+fQf8WeF2xlDKZXLkYP7mf/vKWQdPY9MdwccYZmE++JoemFDLr3lH3bHq5T0TvRiCHU6CfUeBaP3gihlt+VffcbWr9aRJd3Idg4mTA4SEbqAgfePJKqdXm5blrSAFJPDfwwjdeyOoVSVlJGWxlf/fom0XXXwuNuDsy8uzzrqtcni6gfHEeqsGneC200LSDGJIxlvaEu7YyhVpWxZ/Ssr/zsPX2YXssIvwhmahNu/mJgb+3H+oH/aHa/K0QJSTOJOJ9tEsGvjOlp17Gx3HKUqLZ/Xy/ypU0iMyyQzPAbjuBhX9jZqRa5n2CP/pEakXnprFy0gxRRWO4T0Y7BlxXItIEqVka+nvMbhVWF43J0JcWYQ7l1Fi0FNueSmcr3wUhVAC0gxRZ7TiKRjcDzhkN1RlKp0fvpoFn/MP4InoguhoSeIcCxk8MOjad768rO/WZUbLSDF1PaCXiT8nkLmca/dUZSqNFb/9ANrZ8aREd4TR3hd3NmLuOSRUTRvPczuaCofWkCKKapDF0J8P2Cy3XZHUarC27E2jiX/72syQ3pjwrrhzlpOr7EX07HX03ZHU2egBaSYQp1OQn1HwV/L7ihKVViH9+7m22emkZndm2znRbg8q+kwoi29L59kdzRVCFpASkD8xzCOenbHUKrCOZGUxGePv0hWane8YRfj8m6gVV83Q2LzDgukgpkWkBKQ0GSyQtri83q1QzalCsGXlcVHjz9Nxp8dyAofjCt7J/Xa7OSaBx+yO5oqBi0gJSCuTPz+cLb/vor2PfvYHUepoDb3uWdJ3tqUTPdAws0B6jRcwvBHJ+qXrwpMC0gJuOqFkX4YdvwWpwVEqQIsnDWd3Qsy8ET0xuk4Ro1qCxn+n4m4q1WzO5oqIS0gJVD33GYcOwzJexLtjqJU0ElLSeaTh14kw/TD4cwkwrGQqybfTd2G19kdTZUSLSAl0P7Ci9i+cj9ZyaUyUq5SlcaPH87gjx/8ZLovxuX5nV539KLTBc/ZHUuVMi0gJdC8dVtCvVsx2RF2R1EqKKSlJPPxv17E4+9HSGgmtSIXMWqq3stRWWkBKaEQ31Ewei+IUos+msWu77LwuC/G5VlLjzu60bm3Fo/KTAtICQnH8Dua2B1DKdtkpKXx0bhnycjujyM0k5q1FnHj64/r1VVVgBaQEpKQFLyhncjyeAhzueyOo1S5WvzJB+z81oPHPRhX5jp63N6Fzhdqq6Oq0AJSQo5qXozXyaZflhEzYLDdcZQqFydbHR5fP8TppUaNH/nb609oq6OK0QJSQq76LtIOQMLva7WAqCrh509ns+PrtECrI2s93W+NpstFz9odS9lAC0gJNWjTkqMHIGX/MbujKFWmMtLS+OihZ/F4+yHOmtSorq2Oqk4LSAlF9xvA5sXb8aWI3VGUKjNLP/+YbV8m44kYjCtrA93GdCCmv7Y6qjotICVUv0lzQrPiMNnaLYOqnGY+MIG01H5IWDVqVNNWh/ofLSClIDT7KFDb7hhKlaq05CQ+vv81MsIH48rcStdbzuP8gdrqUP+jBaRUHMfvOMfuEEqVmm1rVrH01dV4IvriylrKDW/cT7WaesOsOpUWkFIgzhNkhdQmLSVZ/5OpCm/BzHdJWFQdr6sVNaot5Oap2oeVyp/D7gCVQUg1H0gI65YusjuKUiUy+9FJ7FjaGCNOmkVv5ub/aPFQBdMCUgoiGlUHYP/GLTYnUap4sjwe3vvHoxw/0h9n1iF63lyTq+653+5YKsjpIaxS0LhDaw7vhrSDKXZHUarIDvyxne8mfY0nYhCujN+4+oVR1G3U1O5YqgLQFkgp6NL/EjB+sk+E2B1FqSJZ+e2XzH9yJR53Z9yykNHvPKjFQxWatkBKQY3ISJxZSRh/dbujKFVon734Ike2tsWE1qFe81+4/lE936GKRgtIKQnJPgqi94Ko4GeM4f17J5CeOQhndiJtrhD6X/eo3bFUBaQFpLTIcbJDWtudQqkzSj56hM/++TYZrktweTZw8YRBRLXrZHcsVUFpASklEpaG11GL44mHqF2/kd1xlDrNhpXL+HXqJjzu3ri9i7nxrYdxV9MueFTx6Un0UhJaIxvEwXq9F0QFofn/fYsV0w6RFdaCmrV+Ysy7T2nxUCVmSwERkeEislFE/CLS/QzLXSYiW0Vkh4g8kmt6SxH51Zr+sYiElU/yglVrEgnAoc077Q2iVB6zH53E7lUtQYQWPXZx0wvP2B1JVRJ2tUA2ANcASwpaQERCgCnAUKADcKOIdLBmvwC8Yow5DzgO/L1s455diy6B48jph9NtTqLU/3zw8KMcT+yLM2s/F9xaj8tvu9PuSKoSsaWAGGM2G2O2nmWxnsAOY8wuY0wW8BEwTEQEGAR8ai33PnB1mYUtpE4X9kf8PrJTtZtrFRxmPjCB5KQBuDy7uPTRvnS+sL/dkVQlE8znQJoCe3O93mdNqwskGWN8eabnS0RuE5E4EYlLTEwss7DuatVweo+Dr0aZbUOpwppx9wROZAzG5dnKlc9cRvPWbe2OpCqhMrsKS0QWAvldjjTRGPNVWW03L2PMNGAaQPfu3U1ZbsuRfRSj94Iom02/YwLpZjCu9PX89aXh1GnQxO5IqpIqswJijBlcwlXsB5rnet3MmnYUiBSRUKsVcnK67cSRhC+kw9kXVKoMGGOYfttjZIQMxpX+O9e9MppadevZHUtVYsF8COs3oLV1xVUYcAMwzxhjgEXAddZyo4Fya9GcUVgaPmdNDu3WK7FU+fJ5vUz/xyQyQgIdIt7wxq1aPFSZs+sy3r+KyD6gN/CtiHxvTW8iIvMBrNbF3cD3wGbgE2PMRmsVDwMPisgOAudE3i3vnyE/zkgBYOOypTYnUVWJz+tl5m3PkRE6AJdnJTe9fZ8ObKbKhS13ohtjvgC+yGf6AeAvuV7PB+bns9wuAldpBZUaTeuSkgKJ23fbHUVVEVkeDx+M/Q8ZrotwZy7j5v+OJ9SpVwKq8hHMh7AqnFbdYwDIOOKxN4iqEtJSkpl1+/8jw9Ubl/dnbp6mxUOVr0IVEBH5j4h0LOswFV2Hnn1xZGdh0sPtjqIquRNJSXx89zt43D1wZ//E6KmPavFQ5a6wLZDNwDSr+5CxIqIHWPMR6nQS6j2G8enuUWXneOIhPrn3fTIiuhIhCxnz32e0eChbFKqAGGPeMcb0AW4GooB1IjJbRAaWZbiKyOHXe0FU2Tm8dzef//MzPBHRRIQu4Ja3dBAoZZ9CnwOx+qZqZz2OAGsJXAn1URllq5DEkYzPWdfuGKoS2rt9K19P/A6Puy3VXQu55Y3n7Y6kqrjCngN5BdhK4Aqp54wx3YwxLxhjrgS6lmXAikbc6WSHViNhywa7o6hKZNfGdfzwzFI87nOpUXMRo1/VloeyX2FbIOuALsaY240xq/LMC7rLae0UFhkCwOZly2xOoiqLbWtW8dOLa8gMb0Fk3SXc/O9n7Y6kFFD4+0DWAm0DHeHmSAZ2G2OSSz1VBRYZ1Yik43A84aDdUVQlsHf7Vpb8v41khTehTtOV3DDpabsjKZWjsAXkTeB8Ai0RAToBG4FaInKHMeaHMspX4bTu0YOE31PJPOa1O4qq4I4nHuKHp38k092aOg2Xc8OkJ+yOpNQpCnsI6wDQ1RjT3RjTjcB5j13AJcCLZRWuImoV3ZUQXwYmw2V3FFWBZXk8fD7uAzwR7agesZgbn3rC7khKnaawBaRNrn6oMMZsAtpZXYqoXEKdTkJ9RzHZei+IKh6f18usO/6Nx30+EfzI6Ff0nIcKToU9hLVJRN4iMCogwPXWtHBAj9XkIf7jGEcdu2OoCmrmnU/iCR+EK2spN739hN1xlCpQYVsgo4EdwP3WYxcQS6B46M2EeUhIMl5nPXxera2qaGbcMz6nS/ab3nxY7zBXQe2sLRDrBsL5xpiBwH/yWSS11FNVcOL24PeHs2PdGtp162V3HFVBfDj+MdKyLsaVsZERr91OmEvPo6ngdtYWiDEmG/Br/1eFF1438K1x+695b5lRKn9zJz9P8tG+hHsSuPLZq6gRGWl3JKXOqrDnQFKB9SKyAEg7OdEYc2+ZpKrg6rZqxvFESN592O4oqgL4dtqbHN3RGafvCAMejKFB83PsjqRUoRS2gHxuPVQhtO/Thx2/HsKb7Lc7igpySz6fw75fm+IwGZx/Ux3O69Ld7khKFVqhCogx5n0RcQMtjDFbyzhThdeiTQdCvDswfrfdUVQQW7v0J7Z+HQIO4bwhGXS7+Dq7IylVJIXtTPFKIB74P+t1jIjMK8NcFV6o7xgmO9LuGCpI7dq4jt/e3YcvtAZNzv+DQTfcZHckpYqssJfxPkGg08QkAGNMPNCqTBJVEsJR/HoviMpH4oG9LHrxF7LCGlE3ajVX3nWf3ZGUKpbCFhBvPp0m6gH+M5CQE3jD6pLl0fHR1f+kpSQz75HP8LhaUSNyCSMmPmp3JKWKrbAFZKOI/A0IEZHWIvI6sKIMc1V4jmpZGIeTTb9ot+4qwOf18tE9b+KJ6EyE8yduevEZuyMpVSKFLSD3AB2BTGAOkELgjnRVAFf9wE1gCfHrbE6igsX7Y5/B4+6F27eYW97QAaFUxVfYq7DSgYnWQxVCg9ZRHD0AJ/YdtTuKCgLT75yAxzkYl2clN//3MbvjKFUqClVARKQNMA6Iyv0eY8ygsolV8UX3H8jmn7fjTZGzL6wqtVn/epR0/2Bc6Wv521v3a/9WqtIo7I2Ec4GpwDtAdtnFqTzqN2lOaFYcxl/N7ijKRl+9/gonki7ClbmDq18cgbuafh5U5VHYAuIzxrxVpkkqoZDsY2Ai7Y6hbLLqh2849HtLnCaJPvd2oW6jpnZHUqpUFfYk+tcicqeINBaROicfZZqsEhCO4Q+pa3cMZYMDf2xn3YfHMA4nrQZnaK/MqlIqbAtktPXvQ7mmGfRmwjMS5wm8jtqkpSRTraZ2ZlxVZHk8fDfpWzLdnajXdCUXj9ST5qpyKlQLxBjTMp+HFo+zCKnmwzhCWL/8Z7ujqHL0wV3PBe71CPmJ6x/X4qEqrzMWEBH5V67nw/PM0wvZzyKiYXUA9q3fZHMSVV5mjptIhnMALs9Kbvp/T9odR6kydbYWyA25no/PM++yUs5S6TTu2BqAtIMpNidR5eHLV18mNaUfrvRtjHj1Lr1cV1V6ZysgUsDz/F6rPKL7DgTjJzs1xO4oqoz98t08/lzfCqc3iYvuP19HFFRVwtkKiCngeX6vVR616tbD6U3CZFW3O4oqQ3u3b2XDx8kYCeXcS7Noc35PuyMpVS7OdhVWFxFJIdDacFvPsV67yjRZJRHiOwpS2+4YqoxkeTz88PT/kenuQP0WvzLoBu1dV1UdZ2yBGGNCjDE1jTE1jDGh1vOTr4t9gFdEhovIRhHxi0iBY3iKyGUislVEdojII7mmzxCRP0Qk3nrEFDdLmZMksvVekErrg7uexxMRTYRzkXbNrqqcwt5IWNo2ANcASwpaQERCgCnAUKADcKOIdMi1yEPGmBjrEV+WYUtCwlLxOmtxPPGQ3VFUKZv54AQynP1xZ67gplf0iitV9dhSQIwxmwsxtnpPYIcxZpcxJgv4CBhW9ulKV2iNbBAHG5YWWCtVBfTlK/8hNbU/rvStDH/lbr3iSlVJdrVACqMpsDfX633WtJOeFZF1IvKKiIQXtBIRuU1E4kQkLjExsayyFqha40gADm7eXu7bVmVj5bdf8ueG83B6j9H/nz30iitVZZVZARGRhSKyIZ9HabQixgPtgB5AHeDhghY0xkwzxnQ3xnSvX79+KWy6aJp3CRx1Sz+cVu7bVqVv7/atbJx7AiMhtB7q57wuBZ7CU6rSK2xfWEVmjBlcwlXsB5rnet3MmoYx5qA1LVNEphMYqyQoRfcZSNynK/CnltmuVuUkIy2NH57+nkx3e+q3WMWAETq+mqragvmv2m9AaxFpSaBw3AD8DUBEGhtjDoqIAFcTOCkflNzVquH0HgN/TbujqBKac8+LeCL6U825kBETq2ZPPl6vl3379uHxeOyOosqAy+WiWbNmOAt5Ts+WAiIifwVeB+oD34pIvDHmUhFpArxjjPmLMcYnIncD3wMhwHvGmI3WKj4UkfoE7keJB8aW/09ReI7sY/hFL+WtyN5/YAIZYYNxZS4ndmrVLB4A+/bto0aNGkRFRRH4/qYqC2MMR48eZd++fbRs2bJQ77GlgBhjvgC+yGf6AeAvuV7PB+bns1yFGkpXwg7iCenHjrVxesy8AvrilZdISxuAy7OFG9641+44tvJ4PFo8KikRoW7duhTlYqNgvgqr0qjXpQ6Ig5Uffml3FFVEvy2Yz+EN5+H0HqX/Q711XBfQ4lGJFfV3qwWkHFz2j7GEek+Q9af+8alIjh7az9oPDmHEyblDvJwX3dXuSEoFFS0g5SDM5SLUtxlvaFsy0vRy3orA5/Xy1SOzyHRHUbPxKgbdeLPdkVQpiI2N5dNPPy3xeuLj45k//7Sj66dISEhg9uzZJd5WYUydOpWZM2eWy7Zy0wJSTsKbpJHtrM7306baHUUVwqwHHifD1RN39iL+9rR2U6JOVdIC4vP5SjXP2LFjufnm8v+SE8yX8VYqfUdfz/yXD3J0gw4uFew++/eLpHsH4cpYz6hpeq9HQZ78eiObDpTu57lDk5pMurJjgfMTEhIYOnQoffv2ZcWKFTRt2pSvvvoKt9tNfHw8Y8eOJT09nXPPPZf33nuP2rVP7wl74cKFTJ48mZSUFF5++WWuuOIKPB4Pd9xxB3FxcYSGhvLyyy8zcODAfKf36dOHxx9/nIyMDJYtW8b48eNp1KgR9913HxA4j7BkyRIeeeQRNm/eTExMDKNHj6Z27dp8/vnnpKamkp2dzbfffsuwYcM4fvw4Xq+XZ555hmHDhpGQkMBll11Gt27dWLNmDR07dmTmzJlEREQQFRXFiBEj+O6773C73cyePZvzzjuPJ554gurVqzNu3DgGDBhAr169WLRoEUlJSbz77rtcdNFFpKenExsby4YNG2jbti0HDhxgypQpdO9e/At7tAVSTqLadSI8cycm+zy7o6gzWPXDNxzZ0pawzEQuefQSwlw6akGw2b59O3fddRcbN24kMjKSzz77DICbb76ZF154gXXr1hEdHc2TT+bfckxISGDVqlV8++23jB07Fo/Hw5QpUxAR1q9fz5w5cxg9enSB0/1+P0899RTXX3898fHxXH/99bz00ktMmTKF+Ph4li5ditvtZvLkyVx00UXEx8fzwAMPALBmzRo+/fRTfv75Z1wuF1988QVr1qxh0aJF/POf/8SYwDBLW7du5c4772Tz5s3UrFmTN998Myd/rVq1WL9+PXfffTf3339/vj+jz+dj1apVvPrqqzn74c0336R27dps2rSJp59+mtWrV5f4d6EtkHLkiNhDuhnMqh++oeeQK+yOo/JIPLCXdR8ewTjrc+5l2bRo0+Hsb6rCztRSKEstW7YkJiYGgG7dupGQkEBycjJJSUn0798fgNGjRzN8+PB83z9ixAgcDgetW7emVatWbNmyhWXLlnHPPfcA0K5dO8455xy2bdtW4PS8+vTpw4MPPsjIkSO55ppraNasWb7bvuSSS6hTpw4QuO9iwoQJLFmyBIfDwf79+/nzzz8BaN68OX369AFg1KhRvPbaa4wbF+hw48Ybb8z592Rhyuuaa645Zf8ALFu2LKeV1KlTJzp37pzve4tCWyDl6Jz+7QHY9PUKm5OovHxeL/PGzybT3YJaTX5j0A032R1JFSA8/H99p4aEhBT5fELeS1VL47LkRx55hHfeeYeMjAz69OnDli1b8l2uWrVqOc8//PBDEhMTWb16NfHx8TRs2DDnDv8zZSzoeW4n91Fx9k9RaAEpR/2uvYGwzESykxvZHUXlMev+SXjcPXBnL+LGp56wO44qolq1alG7dm2WLl0KwKxZs3JaI3nNnTsXv9/Pzp072bVrF23btuWiiy7iww8/BGDbtm3s2bPnjNNr1KjBiRMncta5c+dOoqOjefjhh+nRowdbtmw5bZm8kpOTadCgAU6nk0WLFrF79+6ceXv27GHlypUAzJ49m759++bM+/jjj3P+7d27d6H3UZ8+ffjkk08A2LRpE+vXry/0ewuih7DKUajTSQhbyQzvyeG9u2nQ/By7Iyng0xcnk+4biCtjnZ40r8Def//9nJPorVq1Yvr06fku16JFC3r27ElKSgpTp07F5XJx5513cscddxAdHU1oaCgzZswgPDy8wOkDBw5k8uTJxMTEMH78eJYtW8aiRYtwOBx07NiRoUOH4nA4CAkJoUuXLsTGxp52Qn/kyJFceeWVREdH0717d9q1a5czr23btkyZMoUxY8bQoUMH7rjjjpx5x48fp3PnzoSHhzNnzpxC758777yT0aNH06FDB9q1a0fHjh2pVatk96bJyZM2VUH37t1NXFycrRk+efYZEvdeSN1GS7nhiUm2ZlHwy3fzWPupwZGdxqUTYvS8x1ls3ryZ9u3b2x2jUktISOCKK65gw4bT+4iNiooiLi6OevXqFXm92dnZeL1eXC4XO3fuZPDgwWzdupWwsLBTlsvvdywiq40xp12upS2QcnbJbbfy0YTfSU8IsTtKlZd4YC8bPjqGcdan9V+MFg9VqaWnpzNw4EC8Xi/GGN58883TikdRaQEpZ7XrNyIscyu+0Hb4vF4dCtUmPq+Xr8fPIdN1PnUaLWfACG0NquAQFRWVb+sDyLmiqjhq1KhBaR+B0ZPoNgitcxhvWB1+mv2+3VGqrFn3TSLD3R23fzE3PqnFQ6ni0AJig5hrA4M17lv5h81Jqqa5k58nPXsgrvR1jHr9UbvjKFVhaQGxQZeLBhGesQeTEWV3lCpn5bdfcmx7R8IyDzPksaF6p7lSJaAFxCaOsJ14XC3ZtmaV3VGqjMN7d7Pxk2SMOGhzhdC8dVu7IylVoWkBsUn9mHogDn6d/ZXdUaoEn9fL149+QqarKZHNVtP/ur/ZHUkVQ0JCAp06dSrxel599VXS09PPuMyMGTM4cOBAibdVGBdeeGG5bKe0aQGxyaV/v51QbwpZiZF2R6kSZt43CY+7G26zWO+/USUuINnZ2aWaZ8WKitm9kV7Ga5Mwlwtn9haynB1JS0nWoVLL0NznnyMjexAuz1pGva0nzUvNd4/AoZJ3h3GKRtEwdPIZF/H5fIwcOfKUrs5XrlzJuHHj8Pl89OjRg7feeovw8HB+/PHH06a//fbbHDhwgIEDB1KvXj0WLlzI3//+d+Li4hARxowZQ/PmzYmLi2PkyJG43W5WrlxJ+/btuf7661mwYAH/+te/OHHiBNOmTSMrK4vzzjuPWbNmERERQWxsLC6Xi7i4uFO6jJ8xYwZffPEFycnJ7N+/n1GjRjFpUuDLTPXq1UlNTWXx4sU88cQT1KtXjw0bNtCtWzc++OADRIT58+fz4IMPUq1aNfr06cOuXbv45ptvSnf/F5G2QGwU3jiN7NBqfP/ONLujVFrLvprL0Z3RhGf+yZBH/6InzSuBvF2dv/zyy8TGxvLxxx+zfv16fD4fb731Fh6PJ9/p9957L02aNGHRokUsWrSI+Ph49u/fz4YNG1i/fj233HIL1113Hd27d+fDDz8kPj4et9sNQN26dVmzZg033HAD11xzDb/99htr166lffv2vPvuuzkZ8+syHmDVqlV89tlnrFu3jrlz5+Z7X8bvv//Oq6++yqZNm9i1axfLly/H4/Fw++23891337F69WoSExPLZ2efhbZAbNTHGmQqaWOq3VEqpT3bNrHliywIddHur2F60ry0naWlUFbydnX+9NNP07JlS9q0aQMEunKfMmUKAwcOzHd63jE0WrVqxa5du7jnnnu4/PLLGTJkSIHbvv7663Oeb9iwgUcffZSkpCRSU1O59NJLc+bl12U8BLpzr1u3LhDocn3ZsmWnDejUs2fPnO7gY2JiSEhIoHr16rRq1YqWLVsCga7cp02z/4untkBsFNWuE+GeHfh1kKlSl+XxsOCZBWSGN6DuuevpOyz/sSFUxZO3C/PIyMgSra927dqsXbuWAQMGMHXqVG699dYCl83dHXtsbCxvvPEG69evZ9KkSTmtjPwynnxdmK7kS9pdfXnSAmIzqb6XTFdTfvlunt1RKpUP7n4OT0Q0Ec6fGD5+gt1xVCnK29V59+7dSUhIYMeOHcD/unJv27ZtvtOBU7paP3LkCH6/n2uvvZZnnnmGNWvWnLZMfk6cOEHjxo3xer05Xb6flF+X8QALFizg2LFjZGRk8OWXX+a0pM6mbdu27Nq1K6crk5NduttNC4jNWg0IjOq2Zf4vNiepPGb961EyQgfg9qzkpleesjuOKmUnuzpv3749x48f54EHHmD69OkMHz6c6OhoHA4HY8eOxeVy5Tsd4LbbbuOyyy5j4MCB7N+/nwEDBhATE8OoUaN4/vnngUALY+zYscTExJCRkXFajqeffppevXrRp0+fU7pih/91GT906NCcLuMhcHjq2muvpXPnzlx77bWFHo/c7Xbz5ptv5oyVXqNGjRJ3xV4atDt3m/m8Xqbf9hmO7ET+PvMeu+NUeN++PYU9cecSlrWH616+hlp1i97ttSqYdud+drGxsVxxxRVcd911p0yfMWMGcXFxvPHGG8Vab2pqKtWrV8cYw1133UXr1q0LHNK2JIrSnbu2QGwWGGRqG1nhrTm8d/fZ36AKtG75Ivb/2ogQXyq9xnbQ4qEqlf/+97/ExMTQsWNHkpOTuf322+2OpC2QYDD3uec4vOcC6jRcqj3DFtPxxEN8Pu4bssKaEtVrN0NvHWt3pEpJWyCVn7ZAKpght9+GI9tD+m69qro4fF4vXzw0A4+7FTXqrdDioVQ50QISBGrVrUdY1layHYFBplTRzLpvEhmunrizFzHq+aftjqNUlaEFJEiE1knEG1abHz+YYXeUCuWTZ58hPXsQroy1jHp9ot1xlKpStIAEia7XXQLA/pUJ9gapQJZ8Podjf8QQnnmQoU9cod2UKFXOtIAEic59BgYGmcqMsjtKhZCwZQNb5xnE+Ol4bXWatGxtdySlqhwtIEHk5CBTW1b/aneUoJaRlsaPzy3CG1afum030fvKv9odSakqSS/7CSL1uzZgzzoHq+bMo123XnbHCVpz7n0RT0R/qjkXcN2/nrc7TpX1wqoX2HJsS6mus12ddjzc8+EC5yckJDB06FD69u3LihUraNq0KV999RVDhw7lpZdeonv37hw5ciSne5MZM2bw5ZdfkpaWxvbt2xk3bhxZWVnMmjWL8PBw5s+fT506dRgwYABdunTh559/xufz8d5779G9e3fatm3LihUrqF+/Pn6/nzZt2rBy5Urq169fqj93RaUtkCBy6Zh/EJqVgi+xtt1RgtascRPJcPbH7VlB7OtaPKqi7du3c9ddd7Fx40YiIyP57LPPzrj8hg0b+Pzzz/ntt9+YOHEiERER/P777/Tu3ZuZM2fmLJeenk58fDxvvvkmY8aMweFwMGrUqJx+rhYuXEiXLl20eORiSwtERIYDTwDtgZ7GmHzv7hOR94ArgMPGmE65ptcBPgaigARghDHmeNmmLnthLheh2ZvJckbrIFP5+PrN1zmR0g+XZxvDX7vb7jhV3plaCmWpZcuWxMTEANCtW7ecDgYLMnDgQGrUqJHTf9SVV14JQHR0NOvWrctZ7sYbbwSgX79+pKSkkJSUxJgxYxg2bBj3338/7733HrfcckuZ/EwVlV0tkA3ANcCSsyw3A7gsn+mPAD8aY1oDP1qvKwVX03SyQyP4v2lv2x0lqKxd+hMH45ri9CbT576u1ChhF96q4sqvu/PQ0FD8fj/AKd2q513e4XDkvHY4HKd0lZ5fV+vNmzenYcOG/PTTT6xatYqhQ4eW+s9TkdlSQIwxm40xWwux3BLgWD6zhgHvW8/fB64uvXT2uij2RsTvI3lzmt1RgsbRQ/v57Z0/yA5xc86AFD0/pE4TFRXF6tWrAfj000+LtY6TXaQvW7aMWrVq5fR2e+uttzJq1CiGDx9OSEhI6QSuJCrqOZCGxpiD1vNDQMOCFhSR20QkTkTigmUYyDNp0aYD4Z6d+LPb2B0lKGSkpfHlw3PIdLekZoNfGDK64MF+VNU1btw43nrrLbp27cqRI0eKtQ6Xy0XXrl0ZO3bsKcPTXnXVVaSmpurhq3yUWWeKIrIQaJTPrInGmK+sZRYD4wo6B2ItEwV8k+ccSJIxJjLX6+PGmLOeeQ7WzhTzmn7nBNL9gzn/8uQqfYmqz+tl5j9eIMN1IRGykFvees7uSFVeZe1MccCAATlXceUVFxfHAw88wNKlS21IVv6CojNFY8xgY0ynfB5flcLq/xSRxgDWv4dLYZ1B49zB0QBs+X6VzUnsNfOOp8lwXYjb+7MWD2WLyZMnc+211+YMMqVOVVEPYc0DRlvPRwOlUZSCRr9rbiQs80/8KU3tjmKb6XdOyBlV8Oapj9odR1Vyixcvzrf18cgjj7B792769u1rQ6rgZ0sBEZG/isg+oDfwrYh8b01vIiLzcy03B1gJtBWRfSLyd2vWZOASEdkODLZeVyohso1MV2sO7d5pd5Ry9/4DgUN4rvTfuXHK/YQ6nXZHUkrlw66rsL4wxjQzxoQbYxoaYy61ph8wxvwl13I3GmMaG2Oc1vLvWtOPGmMuNsa0tg6V5XelVoVWo7UT43Cy6N1ZdkcpV7MfnURq+iBc6Zu55j8jcVerZnckpVQBKuohrEpvyG23EeLLIH1P1fn2/enk50k6fCHhnt0MfXIItevndw2GUipYaAEJUrXq1sOZta3KDDL1zVuvk7izC2FZhxnwYBftXVepCkALSBALrRsYZGrB+++efeEK7Kc5M9kfdw6hvlS6xzbmvC6nn8xUqjgGDBhARbh0v6LSAhLETg4ydXDVXpuTlJ1V33/NjgXVwPhpdyXEDBhsdySlVCFpd+5BrHOfgaya9h5+WtkdpUxsWLmMtXNS8YfWoGWfg1x0zT/sjqSK4NBzz5G5uXS7cw9v345GEyYUOD8hIYHLLruMbt26sWbNGjp27MjMmTNZuXIl48aNw+fz0aNHD956661T+sB67733WLduHa+++ioA//3vf9m0aROvvPJKqeavarQFEuQcYX+Q6TqHzauW2x2lVCVs2cAvU3fgc9amUcetXHqLFg9VOFu3buXOO+9k8+bN1KxZk5dffpnY2Fg+/vhj1q9fj8/n46233jrlPSNGjODrr7/Ga51PnD59OmPGjLEjfqWiLZAg16BbQ3bHO/jtk/m079nH7jilIvHAXn58bjlZrlbUbfYLVz/wmN2RVDGcqaVQlpo3b06fPoH/C6NGjeLpp5+mZcuWtGkT6D9u9OjRTJkyhfvvvz/nPdWrV2fQoEF88803tG/fHq/XS3R0tB3xKxVtgQS5IbG3EpqVjPdwXbujlIoTSUnMe+RzPO5zqVH7Z65/XIuHKpq83a5HFrJr/1tvvZUZM2Ywffp07RixlGgBCXJhLhehrMMTEcPMB+35xldasjwe5t43DU9ENNXCfuSmF56xO5KqgPbs2cPKlSsBmD17ds7wtTt27ABg1qxZ9O/f/7T39erVi7179zJ79uycwaNUyWgBqQCueeF2XOlbOJE2kI+eeNLuOMXi83r54I4XyXB3x+3/UYejVcXWtm1bpkyZQvv27Tl+/DgPPPAA06dPZ/jw4URHR+NwOBg7dmy+7x0xYgR9+vShdm0dNro06DmQCqB2/UYMeexivn9mJcf39+Kr119h2D0P2B2rSGbe+SQZ4YNwZS3l5refsDuOqsBCQ0P54IMPTpl28cUX8/vvv5+27OLFi095vWzZMh54oGL93wlm2gKpIJq3bkvP21sS6kvmYPy5/DRnpt2RCiXL4+G92yaSETIIt2cVo9+eoJ0jqnKXlJREmzZtcLvdXHzxxXbHqTS0BVKBdL6wPymJc9k4z8+OBRFENvyB8wcNsTtWgbatWcXSV3/BE3ExrozVXP/6HVo8VIlERUWxYcOGIr8vMjKSbdu2lUGiqk1bIBVM32HDadbrANkhEayZ+Sc71p/ebA8GX73+Cotf30Omqx0RoQsY/c79VKtZy+5YSqlSpAWkArr8tjup1yqerLCG/PzSbyQeCJ6uTnxeL9PvnMD+9Z0APy3O38YtbzyvLQ+lKiEtIBXU8PETqFl7CR5XK+Y98hlpKcl2R2LXxrW8//c3SPcPJtyzkUEPtuGKsXfbHUspVUa0gFRgo154hgjnT3giOvPRvW/Y2u37N1Pf4MeXtpLp6kSEYyGj372L86JjbMujlCp7WkAquFveeA63dzEeV29m3vlUuW/f5/Uy4+4J7FnTGiMOmkZv4JY3n9NDVkpVAVpAKoGbpz6GK+NXMkIGMuOe8eW23T1bN/H+318jzTeYcM9WBtxzToW7P0VVbE888QQvvfQSjz/+OAsXLgRg6dKldOzYkZiYGDIyMnjooYfo2LEjDz30UJnnmTFjBnfffebDtoVZpqLQy3grgVCnkxtev5OP7n6fNPfFfDjhMUY+93SZbvO7d6ayb3kkWe7ORPAjI6c9RpjLVabbVMFl6SfbOLI3tVTXWa95dS4a0abI73vqqf+1vj/88EPGjx/PqFGjAJg2bRrHjh0jJCSk1HKqAG2BVBLVatbiymeH4fIkkHykL5/9+8Uy2Y7P62XGfRP4Y1VLjITRpN1abpn6rBYPVW6effZZ2rRpQ9++fdm6dSsAsbGxfPrpp7zzzjt88sknPPbYY4wcOZKrrrqK1NRUunXrxscff5zv+mJjY7njjju44IILaNWqFYsXL2bMmDG0b9+e2NjYnOXmzJlDdHQ0nTp14uGHH86ZPn36dNq0aUPPnj1Zvvx/wy4kJiZy7bXX0qNHD3r06HHKvMpCWyCVSIPm59D/wRgWvbKVxK0d+b/3pnHZmNtKbf0H/tjOd098icc9GJdnExfe3ZX2Pa8ptfWriqU4LYWSWr16NR999BHx8fH4fD7OP/98unXrljP/1ltvZdmyZVxxxRVcd911QKAr9/j4+DOu9/jx46xcuZJ58+Zx1VVXsXz5ct555x169OhBfHw8DRo04OGHH2b16tXUrl2bIUOG8OWXX9KrVy8mTZrE6tWrqVWrFgMHDqRr164A3HfffTzwwAP07duXPXv2cOmll7J58+Yy2zd20AJSyZzXpTvHRx7m9zkZ7F5Wn5UNv6T35VeXeL0LZr5Lwk9uslxdcft/ZJQeslI2WLp0KX/961+JiIgA4KqrriqV9V555ZWICNHR0TRs2DBnrJCOHTuSkJDA7t27GTBgAPXr1wdg5MiRLFmyBOCU6ddff33OHe8LFy5k06ZNOdtISUkhNbV0D/nZTQtIJdTjkr+Q/OcMdiyux4a5J4isv7zQg1EdTzzExhVLObR1J+mHT+BLAZMVQZazGw5HBg3PjeO6fz1bxj+BUuXr5PC3DofjlKFwHQ4HPp8PZzGuKvT7/fzyyy+4KvEXLS0gldTgUbGkJv6Hg5s6snzKVmrVb0Ctug3Y9Oty9m/cQurBY3hT/BiPC5NdEyORZIdG4nPWBOpaD8ABoSGphGVuocfYrnS+8BE7fyxVxfXr14/Y2FjGjx+Pz+fj66+/5vbbby/z7fbs2ZN7772XI0eOULt2bebMmcM999xDz549ue+++zh69Cg1a9Zk7ty5dOnSBYAhQ4bw+uuv51z9FR8fT0xMTJlnLU9aQCqxqx/4J3Mef4Jjf/Zl3rNbyQ7ZDeICYnKWCQnJINQkIf7jOLMP4HSkE1LdT0SD6jRoHUX7Xn1o0Pwc234GpXI7//zzuf766+nSpQsNGjSgR48e5bLdxo0bM3nyZAYOHIgxhssvv5xhw4YBgUuJe/fuTWRk5CkF4rXXXuOuu+6ic+fO+Hw++vXrx9SpU8slb3kRY4zdGcpN9+7dTVxcnN0xyt3sRyeRsdeFODNwVPPhquuiblRTWvfsRVS7TnbHUxXI5s2bad++vd0xVBnK73csIquNMd3zLqstkCrgb89UzFEMlVLBTQuIUqpKePbZZ5k7d+4p04YPH87EiRNtSlTxaQFRShWJMQYRsTtGkU2cOFGLxVkU9ZSG3omulCo0l8vF0aNHi/yHRgU/YwxHjx4t0mXH2gJRShVas2bN2LdvH4mJiXZHUWXA5XLRrFmzQi+vBUQpVWhOp5OWLVvaHUMFCT2EpZRSqli0gCillCoWLSBKKaWKpUrdiS4iicBuu3MUoB5wxO4QZ6D5SkbzlYzmK5mS5jvHGFM/78QqVUCCmYjE5ddVQLDQfCWj+UpG85VMWeXTQ1hKKaWKRQuIUkqpYtECEjym2R3gLDRfyWi+ktF8JVMm+fQciFJKqWLRFohSSqli0QKilFKqWLSAlAMRuUxEtorIDhE5bVBxEXlQRDaJyDoR+VFEzsk1L1tE4q3HPJvyxYpIYq4ct+aaN1pEtluP0TbleyVXtm0ikpRrXpnuPxF5T0QOi8iGAuaLiLxmZV8nIufnmlce++5s+UZaudaLyAoR6ZJrXoI1PV5EymQoz0LkGyAiybl+h4/nmnfGz0U55XsoV7YN1uetjjWvPPZfcxFZZP392Cgi9+WzTNl9Bo0x+ijDBxAC7ARaAWHAWqBDnmUGAhHW8zuAj3PNSw2CfLHAG/m8tw6wy/q3tvW8dnnny7P8PcB75bj/+gHnAxsKmP8X4DtAgAuAX8tr3xUy34UntwsMPZnPep0A1LN5/w0Avinp56Ks8uVZ9krgp3Lef42B863nNYBt+fz/LbPPoLZAyl5PYIcxZpcxJgv4CBiWewFjzCJjTLr18heg8P0pl0O+M7gUWGCMOWaMOQ4sAC6zOd+NwJxSzlAgY8wS4NgZFhkGzDQBvwCRItKY8tl3Z81njFlhbR/K/7NXmP1XkJJ8bgutiPnK9bMHYIw5aIxZYz0/AWwGmuZZrMw+g1pAyl5TYG+u1/s4/Rec298JfFs4ySUicSLyi4hcbWO+a63m76ci0ryI7y2PfFiH/loCP+WaXNb772wKyl8e+66o8n72DPCDiKwWkdtsygTQW0TWish3ItLRmhZU+09EIgj88f0s1+Ry3X8iEgV0BX7NM6vMPoM6HkgQEZFRQHegf67J5xhj9otIK+AnEVlvjNlZztG+BuYYYzJF5HbgfWBQOWcojBuAT40x2bmmBcP+C3oiMpBAAemba3Jfa981ABaIyBbrG3l5WkPgd5gqIn8BvgRal3OGwrgSWG6Myd1aKbf9JyLVCRSv+40xKWWxjfxoC6Ts7Qea53rdzJp2ChEZDEwErjLGZJ6cbozZb/27C1hM4BtGueYzxhzNlekdoFth31se+XK5gTyHEMph/51NQfnLY98Vioh0JvB7HWaMOXpyeq59dxj4gsBho3JljEkxxqRaz+cDThGpRxDtP8uZPntluv9ExEmgeHxojPk8n0XK7jNYlid49GEg0MrbReDQysmTfR3zLNOVwAnB1nmm1wbCref1gO2U8onCQuZrnOv5X4FfrOd1gD+snLWt53XKO5+1XDsCJy2lPPefte4oCj4JfDmnnsBcVV77rpD5WgA7gAvzTK8G1Mj1fAVwmQ35Gp38nRL4A7zH2peF+lyUdT5rfi0C50mqlff+s/bFTODVMyxTZp9BPYRVxowxPhG5G/iewJUj7xljNorIU0CcMWYe8G+gOjBXRAD2GGOuAtoDb4uIn0BrcbIxZpMN+e4VkasAH4H/KLHWe4+JyNPAb9bqnjKnNuHLKx8EvgF+ZKz/GZYy338iMofAlUL1RGQfMAlwWtmnAvMJXAWzA0gHbrHmlfm+K2S+x4G6wJvWZ89nAr22NgS+sKaFArONMf9nQ77rgDtExAdkADdYv+N8Pxc25IPAl6ofjDFpud5aLvsP6APcBKwXkXhr2gQCXwzK/DOoXZkopZQqFj0HopRSqli0gCillCoWLSBKKaWKRQuIUkqpYtECopRSqli0gKgKTURSC7HM/VZXE6W1zatFpEMprm9FCd6bav3bREQ+PcNykSJyZ3G3o1R+tICoquB+oEgFRERCzjD7aqDUCogx5sJSWMcBY8x1Z1gkEtACokqVFhBVKVjjRiy2OnvcIiIfWuMg3As0ARaJyCJr2SEislJE1ojIXKsfoZPjN7wgImuA4SLyDxH5zerI7zMRiRCRC4GrgH9b4zycKyIxVmeN60TkCxGpba1vsQTGKokTkc0i0kNEPrfGXngmV/bUXM8flsAYEmtFZHI+P2dLK/v6POuIEmvMChHpKCKrrHzrRKQ1MBk415r2bxGpLoGxZ9ZY6xqWaz2bReS/Ehhf4gcRcVvzzhORhVa2NSJyrjX9IWs/rRORJ0v1F6uCW1nc+q8PfZTXA2u8DwJ3CycT6M/HAawk0Jkd5BqXgUCXJkuwup0AHgYez7Xcv3Ktu26u588A91jPZwDX5Zq3DuhvPX8Kq1sJAn1vvWA9vw84QGD8hnACPZ/WzfMzDCXQ5cXJsWFO61YCmAfcbD2/K9d7o7C62wBeB0Zaz8MAN3m64yBwd3TNXPtkB4GuLqII9DgQY837BBhlPf8V+Kv13EWgVTcEmGa91wF8A/Sz+3Ohj/J5aFcmqjJZZYzZB2B16xAFLMuzzAUEDj8tt7qZCCNQbE76ONfzTta3/EgCXc18n3eDIlILiDTG/GxNeh+Ym2uRk12trAc2GmMOWu/bRaAju6O5lh0MTDfW2DAm/24l+gDXWs9nAS/ks8xKYKKINAM+N8Zst37WU6IDz4lIP8BPoBvvhta8P4wx8dbz1UCUiNQAmhpjvrCyeayfYwiBIvK7tXx1Ar3llnevvcoGWkBUZZKZ63k2+X++hcAgOjcWsI7c/RnNAK42xqwVkVgCrZziZvLnyecvIF9hnLH/IWPMbBH5lUAnevMl0AX/rjyLjQTqA92MMV4RSSDQqsidGQL70X2GzQnwvDHm7SLkV5WEngNRVcEJAsN9QmDUvT4ich6AiFQTkTYFvK8GcFAC3WWPzG99xphk4LiIXGTNuwn4meJZANxy8ooxscbWzmM5gY4jyZMphwTGPtlljHkN+ArozKn7AAI9yB62isdA4JwzBTOB0e72iTUol4iEWzm/B8bkOo/UVALjX6gqQAuIqgqmAf8nIouMMYkEehOeIyLrCBzuaVfA+x4jcNx/ObAl1/SPgIdE5HfrRPJoAifV1wExBM6DFJkJ9NY6D4izDsGNy2ex+4C7RGQ9BY8eNwLYYK2jE4HhTI8SOGy3QUT+DXwIdLfWc3Oen68gNxHomXkdgXM1jYwxPwCzgZXWuj7l1EKlKjHtjVcppVSxaAtEKaVUsWgBUUopVSxaQJRSShWLFhCllFLFogVEKaVUsWgBUUopVSxaQJRSShXL/weOjduqvPcoeAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "for value, bootstrap in conditions.items():\n",
    "    plt.plot(results[f'points_{bootstrap}'], results[f'energies_{bootstrap}'], label = f'{bootstrap}')\n",
    "plt.plot(results['points_np'], results['energies_np'], label = 'numpy')\n",
    "for condition in extrapolators.keys():\n",
    "    print(condition)\n",
    "    plt.plot(results[f'points_{condition}'], results[f'energies_{condition}'], label = condition)\n",
    "plt.legend()\n",
    "plt.title('Dissociation profile')\n",
    "plt.xlabel('Interatomic distance')\n",
    "plt.ylabel('Energy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare number of evaluations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no bootstrapping\n",
      "Total evaluations for no bootstrapping:\n",
      "1714\n",
      "bootstrapping\n",
      "Total evaluations for bootstrapping:\n",
      "1266\n",
      "np\n",
      "Total evaluations for np:\n",
      "0\n",
      "poly\n",
      "Total evaluations for poly:\n",
      "715\n",
      "diff_model\n",
      "Total evaluations for diff_model:\n",
      "928\n"
     ]
    }
   ],
   "source": [
    "for condition, result_full in results_full.items():\n",
    "    print(condition)\n",
    "    print('Total evaluations for ' + condition + ':')\n",
    "    sum = 0\n",
    "    for key in results_full[condition].keys():\n",
    "        if condition is not 'np':\n",
    "                sum += results_full[condition][key]['raw_result']['cost_function_evals']\n",
    "        else:\n",
    "                sum = 0\n",
    "    print(sum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td>Qiskit</td><td>0.23.0</td></tr><tr><td>Terra</td><td>0.16.0</td></tr><tr><td>Aer</td><td>0.7.0</td></tr><tr><td>Ignis</td><td>0.5.0</td></tr><tr><td>Aqua</td><td>0.8.0</td></tr><tr><td>IBM Q Provider</td><td>0.11.0</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.6.12 |Anaconda, Inc.| (default, Sep  8 2020, 17:50:39) \n",
       "[GCC Clang 10.0.0 ]</td></tr><tr><td>OS</td><td>Darwin</td></tr><tr><td>CPUs</td><td>4</td></tr><tr><td>Memory (Gb)</td><td>16.0</td></tr><tr><td colspan='2'>Mon Oct 19 16:12:19 2020 CEST</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><h3>This code is a part of Qiskit</h3><p>&copy; Copyright IBM 2017, 2020.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "%qiskit_version_table\n",
    "%qiskit_copyright"
   ]
  }
 ],
 "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.6.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}