05_Sampling_potential_energy_surfaces with runtime.ipynb

{
 "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 Nature to generate the dissociation profile of a molecule. We use the Born-Oppenheimer 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": 33,
   "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.utils import QuantumInstance, algorithm_globals\n",
    "from qiskit import Aer\n",
    "from qiskit.algorithms.minimum_eigen_solvers import NumPyMinimumEigensolver, VQE\n",
    "from qiskit.algorithms.optimizers import SLSQP\n",
    "from qiskit.circuit.library import ExcitationPreserving\n",
    "from qiskit import BasicAer\n",
    "from qiskit.algorithms import NumPyMinimumEigensolver, VQE\n",
    "from qiskit.algorithms.optimizers import SLSQP, SPSA\n",
    "\n",
    "# qiskit nature imports\n",
    "from qiskit_nature.problems.second_quantization import ElectronicStructureProblem\n",
    "from qiskit_nature.converters.second_quantization import QubitConverter\n",
    "from qiskit_nature.mappers.second_quantization import JordanWignerMapper, ParityMapper\n",
    "from qiskit_nature.algorithms import GroundStateEigensolver\n",
    "from qiskit_nature.drivers import UnitsType, Molecule\n",
    "from qiskit_nature.drivers.second_quantization import (\n",
    "    ElectronicStructureDriverType,\n",
    "    ElectronicStructureMoleculeDriver,\n",
    ")\n",
    "from qiskit_nature.algorithms.pes_samplers import BOPESSampler, Extrapolator\n",
    "\n",
    "import warnings\n",
    "\n",
    "warnings.simplefilter(\"ignore\", np.RankWarning)\n",
    "\n",
    "algorithm_globals.random_seed = 75"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we use the H2 molecule as a model system for testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "stretch1 = partial(Molecule.absolute_stretching, atom_pair=(1, 0))\n",
    "mol = Molecule(\n",
    "    geometry=[(\"H\", [0.0, 0.0, 0.0]), (\"H\", [0.0, 0.0, 0.3])],\n",
    "    degrees_of_freedom=[stretch1],\n",
    ")\n",
    "\n",
    "# pass molecule to PSYCF driver\n",
    "driver = ElectronicStructureMoleculeDriver(mol, driver_type=ElectronicStructureDriverType.PYSCF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "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": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.5])]\n",
      "[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.8999999999999999])]\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": "markdown",
   "metadata": {},
   "source": [
    "### VQEClient setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import IBMQ\n",
    "\n",
    "IBMQ.load_account()\n",
    "provider = IBMQ.get_provider(group=\"open\")  # replace by your runtime provider\n",
    "\n",
    "backend = provider.get_backend(\"ibmq_qasm_simulator\")  # select a backend that supports the runtime\n",
    "# backend = provider.get_backend(\"simulator_statevector\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.circuit.library import EfficientSU2\n",
    "\n",
    "ansatz = EfficientSU2(num_qubits=4, reps=1, entanglement=\"linear\", insert_barriers=True)\n",
    "\n",
    "np.random.seed(5)  # fix seed for reproducibility\n",
    "initial_point = np.random.random(ansatz.num_parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.runtime import VQEClient\n",
    "\n",
    "# from qiskit_nature.test.runtime import FakeRuntimeProvider\n",
    "\n",
    "# Pass an qiskit.algorithms.optimizer.Optimizer object (like the ``optimizer`` variable from before)\n",
    "# or SPSA/QN-SPSA as dictionary. Here we're using the dictionary version\n",
    "optimizer = {\n",
    "    #     \"name\": \"QN-SPSA\",  # leverage the Quantum Natural SPSA\n",
    "    \"name\": \"SPSA\",  # set to ordinary SPSA\n",
    "    \"maxiter\": 100,\n",
    "    \"resamplings\": {1: 100},  # 100 samples of the QFI for the first step, then 1 sample per step\n",
    "}\n",
    "\n",
    "vqe_c = VQEClient(\n",
    "    ansatz=ansatz,\n",
    "    optimizer=optimizer,\n",
    "    initial_point=initial_point,\n",
    "    backend=backend,\n",
    "    provider=provider,\n",
    "    shots=1024,\n",
    "    measurement_error_mitigation=True,\n",
    ")  # use a complete measurement fitter for error mitigation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "es_problem = ElectronicStructureProblem(driver)\n",
    "\n",
    "qubit_converter = QubitConverter(JordanWignerMapper())\n",
    "qubit_converter = QubitConverter(ParityMapper(), two_qubit_reduction=True)\n",
    "quantum_instance = QuantumInstance(backend=Aer.get_backend(\"aer_simulator\"))\n",
    "solver = VQE(quantum_instance=quantum_instance, optimizer=SPSA())\n",
    "\n",
    "# me_gsc = GroundStateEigensolver(qubit_converter, vqe_c)\n",
    "me_gsc = GroundStateEigensolver(qubit_converter, solver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance1 = partial(Molecule.absolute_distance, atom_pair=(1, 0))\n",
    "mol = Molecule(\n",
    "    geometry=[(\"H\", [0.0, 0.0, 0.0]), (\"H\", [0.0, 0.0, 0.3])],\n",
    "    degrees_of_freedom=[distance1],\n",
    ")\n",
    "\n",
    "# pass molecule to PSYCF driver\n",
    "driver = ElectronicStructureMoleculeDriver(mol, driver_type=ElectronicStructureDriverType.PYSCF)\n",
    "\n",
    "es_problem = ElectronicStructureProblem(driver)\n",
    "\n",
    "\n",
    "# Specify degree of freedom (points of interest)\n",
    "points = np.linspace(0.25, 2, 30)\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(gss=me_gsc, bootstrap=value, num_bootstrap=None, extrapolator=None)\n",
    "    # execute\n",
    "    res = bs.sample(es_problem, 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": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define numpy solver\n",
    "solver_numpy = NumPyMinimumEigensolver()\n",
    "me_gsc_numpy = GroundStateEigensolver(qubit_converter, solver_numpy)\n",
    "bs_classical = BOPESSampler(\n",
    "    gss=me_gsc_numpy, bootstrap=False, num_bootstrap=None, extrapolator=None\n",
    ")\n",
    "# execute\n",
    "res_np = bs_classical.sample(es_problem, 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": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Energy')"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABKeklEQVR4nO3dd5xU5dXA8d+Z2d7ZQu+I0gVBMGIBW7Bhb69E0VhQo6KviSYmxpaob4waexd7w4YJ0QhiQY0ISBcVEZSysI2F7Tszz/vHc2d3WGZnZ2bL7C7n+/nMZ+/ce+feZy7DnHnauWKMQSmllGqMK9YFUEop1b5poFBKKRWSBgqllFIhaaBQSikVkgYKpZRSIWmgUEopFZIGCtWhiMijIvKnGJz33yJyfpSv7SsiZSLibulytTQRuV1ECkUkv2G5ReQjEbko1mVUbU90HoVqL0RkA9AN8ABeYA3wHPC4McYXw6JFxHkfFxlj5sW6LJEQkb7At0A/Y8z2INs/Al4wxjzZ1mVTsaU1CtXenGiMSQf6AXcC1wNPxbZInYOIxDWxS1+gKFiQUHs3DRSqXTLGlBpj5gBnAeeLyAgAEZklIrc7y7ki8k8R2SEixSLyqYi4nG3Xi8hmEdklIt+KyJHO+kQRuU9EtjiP+0Qk0X9eETlJRJaJyE4R+UFEpjjr65pdRGSQiHwoIkVOM82LIpLlbHse+4X7rtNs8zsR6S8ixv9FLSI9RWSOU+Z1InJxwPlvFpHXROQ5p+yrRWRcY9fJOe5VIrLeKcvfAq7BdBH5TETuFZEi4GYRyXSOXSAiG0XkjyLiEpGjgA+Ank65ZzUsd5BzXygi34hIiYi8LyL9ovm3Vu2fBgrVrhljFgGbgEODbP5fZ1setsnqD4ARkf2A3wAHOrWTXwIbnNfcCBwEjAb2B8YDfwQQkfHYpq7fAlnAYQGvCyTAHUBPYCjQB7jZKe+vgJ+wNaM0Y8z/BXn9K065ewKnA38VkSMCtk919skC5gAPBjlGoFOAccABwEnAhQHbJgDrsdfnL8ADQCYwEDgcOA+4wGkmOxbY4pR7eqgTishJ2Ot9Kvb6fwq83EQ5VQelgUJ1BFuA7CDra4Ee2Db1WmPMp8Z2unmBRGCYiMQbYzYYY35wXnMucKsxZrsxpgC4BfiVs+3XwNPGmA+MMT5jzGZjzNqGJzXGrHP2qXaOcQ/2S7dJItIHmAhcb4ypMsYsA57EfmH7LTTGzDXGeIHnsQEtlLuMMcXGmJ+A+4BzArZtMcY8YIzxADXA2cDvjTG7jDEbgL8HvP9IzADuMMZ84xz7r8BorVV0ThooVEfQCygOsv5vwDrgP07Tyw1gv8iBmdhf+dtF5BUR6em8piewMeAYG511YGsGP9AEEenmHHOziOwEXgByw3wvPYFiY8yuBmXoFfA8P2C5Akhqon/h5wbH6tnItlwgnj3ff+C5w9UP+IfT7LcD++8jUR5LtXMaKFS7JiIHYr98Fjbc5vwq/l9jzEBsc821/r4IY8xLxphDsF9oBrjLedkWZ51fX2cd2C/VQWEU66/OMUcaYzKAadgvybqihXjtFiBbRNIblGFzGOdtTJ8Gx9oS8DywLIXYWljD9x/NuX8GLjXGZAU8ko0xn0dxLNXOaaBQ7ZKIZIjICdi2+heMMSuD7HOCiOwjIgKUYpucfCKyn4gc4XRSVwGVgH947cvAH0UkT0RygZuwNQKwo6suEJEjnQ7eXiIyJEjx0oEyoFREemH7NAJtw/YB7MEY8zPwOXCHiCSJyChsk9cLwfYP029FpIvTrHU18Goj5/YCrwF/EZF0p5no2ijP/SjwexEZDuB0kp8RXfFVe6eBQrU374rILuwv1hux7f8XNLLvYGAe9kv7C+BhY8wCbP/Endhf0PlAV+D3zmtuBxYDK4CVwFJnnb/j/ALgXmzg+Zjdf3373YLtOC4F/gW82WD7HdhgtENErgvy+nOA/thf/m8Bf27mnIt3gCXAMqc8oYYTXwmUYzu4FwIvAU9HekJjzFvYWtorTvPbKmxnuOqEdMKdUh2YiBhgsNMvo1Sr0BqFUkqpkDRQKKWUCkmbnpRSSoWkNQqllFIhNZUkrMPJzc01/fv3j3UxlFKqQ1myZEmhMSYv2LZOFyj69+/P4sWLY10MpZTqUERkY2PbtOlJKaVUSBoolFJKhaSBQimlVEidro9CKdU8tbW1bNq0iaqqqlgXRbWCpKQkevfuTXx8fNiv0UChlNrNpk2bSE9Pp3///th8i6qzMMZQVFTEpk2bGDBgQNiv06YnpdRuqqqqyMnJ0SDRCYkIOTk5EdcWNVAopfagQaLziubfVgNFNIrXw7r5sS6FUkq1CQ0U0Vh4H7x5caxLoZQKYvr06cyePbvZx1m2bBlz584Nuc+GDRt46aWXmn2ucDz66KM899xzbXKuhjRQRKNsO1SWgCZUVKrTam6g8Hg8LVqeGTNmcN5557XoMcOlgSIaFYVgfFBTFuuSKNXpbNiwgaFDh3LxxRczfPhwjjnmGCorKwH75X3QQQcxatQoTjnlFEpKSoIeY968eYwbN459992Xf/7zn4DtpL/gggsYOXIkY8aMYcGCBY2ur6mp4aabbuLVV19l9OjRvPrqq3z88ceMHj2a0aNHM2bMGHbt2sUNN9zAp59+yujRo7n33nuZNWsWU6dO5YgjjuDII4+krKyMI488kgMOOICRI0fyzjvv1L3HIUOGcO655zJ06FBOP/10KioqAJuG6He/+x0jR45k/PjxrFtn70l18803c/fddwMwadIkrr/+esaPH8++++7Lp59+CkBFRQVnnnkmw4YN45RTTmHChAktktJIh8dGo7zA/q3aCYnpsS2LUq3olndXs2bLzhY95rCeGfz5xOEh9/n+++95+eWXeeKJJzjzzDN54403mDZtGueddx4PPPAAhx9+ODfddBO33HIL99133x6v37BhA4sWLeKHH35g8uTJrFu3joceeggRYeXKlaxdu5ZjjjmG7777rtH1t956K4sXL+bBBx8E4MQTT+Shhx5i4sSJlJWVkZSUxJ133sndd99dF4xmzZrF0qVLWbFiBdnZ2Xg8Ht566y0yMjIoLCzkoIMOYurUqQB8++23PPXUU0ycOJELL7yQhx9+mOuus3fOzczMZOXKlTz33HPMnDmz7viBPB4PixYtYu7cudxyyy3MmzePhx9+mC5durBmzRpWrVrF6NGjm/EvVU9rFNEoL7J/q0pjWw6lOqkBAwbUfcmNHTuWDRs2UFpayo4dOzj88MMBOP/88/nkk0+Cvv7MM8/E5XIxePBgBg4cyNq1a1m4cCHTpk0DYMiQIfTr14/vvvuu0fUNTZw4kWuvvZb777+fHTt2EBcX/Hf20UcfTXZ2NmDnLfzhD39g1KhRHHXUUWzevJlt27YB0KdPHyZOnAjAtGnTWLhwYd0xzjnnnLq/X3zxRdDznHrqqbtdH4CFCxdy9tlnAzBixAhGjRoV9LWR0hpFpGqroGaXXa5u2V9aSrU3Tf3yby2JiYl1y263u67pKVwNh4C2xHDfG264geOPP565c+cyceJE3n///aD7paam1i2/+OKLFBQUsGTJEuLj4+nfv3/dHIZQZWxsOZD/Grnd7hbvD2lIaxSRqiisX9YahVJtJjMzky5dutS1xz///PN1tYuGXn/9dXw+Hz/88APr169nv/3249BDD+XFF18E4LvvvuOnn34KuT49PZ1du3bVHfOHH35g5MiRXH/99Rx44IGsXbt2j30aKi0tpWvXrsTHx7NgwQI2bqzP5P3TTz/V1RZeeuklDjnkkLptr776at3fX/ziF2Ffo4kTJ/Laa68BsGbNGlauXBn2a0PRGkWkygMDhdYolGpLzz77LDNmzKCiooKBAwfyzDPPBN2vb9++jB8/np07d/Loo4+SlJTE5ZdfzmWXXcbIkSOJi4tj1qxZJCYmNrp+8uTJ3HnnnYwePZrf//73LFy4kAULFuByuRg+fDjHHnssLpcLt9vN/vvvz/Tp0+nSpctu5Tj33HM58cQTGTlyJOPGjWPIkCF12/bbbz8eeughLrzwQoYNG8Zll11Wt62kpIRRo0aRmJjIyy+/HPb1ufzyyzn//PMZNmwYQ4YMYfjw4WRmZkZ4lffU6e6ZPW7cONOqNy5aNw9eOM0uH3c3jNf5FKpz+eabbxg6dGisi9GpbdiwgRNOOIFVq1btsc1/87Xc3NyIj+v1eqmtrSUpKYkffviBo446im+//ZaEhITd9gv2bywiS4wx44IdV2sUkfJ3ZIM2PSml2pWKigomT55MbW0txhgefvjhPYJENGIaKETkaeAEYLsxZkSQ7QL8AzgOqACmG2OWtm0pG/APjQXtzFZKRaV///5BaxNA3QimaKSnp7fKraBj3Zk9C5gSYvuxwGDncQnwSBuUKbSKQnDFQ2qe1iiUUnuFmAYKY8wnQHGIXU4CnjPWf4EsEenRNqVrRHkhpOZCUpZ2Ziul9gqxrlE0pRfwc8DzTc662CkvhJRcSMrQGoVSaq/QKTqzReQSbNMUffv2bd2TVRRCag6IS/solFJ7hfZeo9gM9Al43ttZtxtjzOPGmHHGmHF5eXmtW6LyQts/kag1CqVay4YNGxgxYo/xLRG777776pLtNWbWrFls2bKl2ecKx8EHH9wm52lp7T1QzAHOE+sgoNQYszWmJapresrUPgql2rnmBgqv19ui5fn8889b9HhtJaaBQkReBr4A9hORTSLyaxGZISIznF3mAuuBdcATwOUxKqrlqbZ5nlJztI9CqVbm8Xj2SMM9f/58xowZw8iRI7nwwguprq4GCLr+/vvvZ8uWLUyePJnJkyfj9XqZPn06I0aMYOTIkdx7773Mnj2bxYsXc+655zJ69GgqKyvp378/119/PQcccACvv/46TzzxBAceeCD7778/p512Wl3gmT59OjNmzNgjnfmsWbM46aSTmDRpEoMHD+aWW26pe09paWkAfPTRR0yaNInTTz+9Lt24f/Lz3LlzGTJkCGPHjuWqq67ihBNOaMvLHlRM+yiMMec0sd0AV7RRcZrmT9+R6jRveSrBUwNxzZ/QolS79O8bIL9l8gXV6T4Sjr2zyd0apuG+5557eOyxx5g/fz777rsv5513Ho888ggzZsxg+vTpe6yfOXMm99xzDwsWLCA3N5clS5awefPmuvkLO3bsICsriwcffJC7776bcePqJyXn5OSwdKmdslVUVMTFF9sMDH/84x956qmnuPLKK4Hg6cwBFi1axKpVq0hJSeHAAw/k+OOP3+34AF9//TWrV6+mZ8+eTJw4kc8++4xx48Zx6aWX8sknnzBgwIC6LLKx1t6bntoX/2S7FGd4LGiHtlKtpGEa7vnz5zNgwAD23XdfoD7N+Lfffht0fUMDBw5k/fr1XHnllbz33ntkZGQ0eu6zzjqrbnnVqlUceuihjBw5khdffJHVq1fXbQuWzhxsqvGcnBySk5M59dRTd0sh7jd+/Hh69+6Ny+Vi9OjRbNiwgbVr1zJw4EAGDBgA0G4CRacY9dRm/JljU3OhptwuV5Xa50p1RmH88m8tDdNrZ2VlUVRU1MjeTevSpQvLly/n/fff59FHH+W1117j6aefDrpvYKrw6dOn8/bbb7P//vsza9YsPvroo0bL6H8eTprzhqnUWztVeHNojSIS/jxPqXm2Mxu0n0KpVtIwDfe4cePYsGFDXfOOP834fvvtF3Q9sFsa8MLCQnw+H6eddhq33357XdNSU6nCd+3aRY8ePaitra1LR+4XLJ05wAcffEBxcTGVlZW8/fbbdTWjpuy3336sX7++Lo2HP914rGmNIhL+GkVKDiTZu1RpoFCqdTRMw33//fdz0EEHccYZZ+DxeDjwwAOZMWMGiYmJPPPMM3usB7jkkkuYMmUKPXv25L777uOCCy7A5/MBcMcddwD1ndLJyclB7yZ32223MWHCBPLy8pgwYcJuQSVYOnOwzUqnnXYamzZtYtq0aXv0TzQmOTmZhx9+mClTppCamsqBBx7YrGvYUjTNeCTm3QyfPwh/KoDta+CRg+HM52DYSa1zPqViQNOMh2f69OmccMIJnH766butnzVr1m732o5UWVkZaWlpGGO44oorGDx4MNdcc01LFLlOpGnGtekpEv48TyJ2wh1ojUIp1aKeeOIJRo8ezfDhwyktLeXSSy+NdZG06SkiFUV2xBNoH4VSe7lZs2YFXT99+nSmT58e9XGvueaaFq9BNJfWKCJRXmAn2wEkpNl8Tzo7WynVyWmgiIQ/zxOAywWJ6VqjUEp1ehooIhHY9AS2+Ukn3CmlOjkNFOHyVNug4G96AkjM1BqFUqrT00ARLn+ep4Y1Cu2jUEp1choowlXRICEgaAZZpdReQQNFuPwJAVMb1ig0UCjV0jZs2MDQoUO5+OKLGT58OMcccwyVlZVMmjQJ/4TawsJC+vfvD9ihqieffDJHH300/fv358EHH+See+5hzJgxHHTQQRQXFwMwadIkrr76akaPHs2IESNYtGgRPp+PwYMHU1Bg/4/7fD722WefuudK51GEz5/nKbDpKTEDqjVQqM7rrkV3sbZ4bYsec0j2EK4ff32T+33//fe8/PLLPPHEE5x55pm88cYbIfdftWoVX3/9NVVVVeyzzz7cddddfP3111xzzTU899xzzJw5E4CKigqWLVvGJ598woUXXsiqVauYNm0aL774IjNnzmTevHnsv//+tPrdMjsQrVGEKzBzrJ+/j8LJHaOUajkDBgxg9OjRAIwdO7YuUV5jJk+eTHp6Onl5eWRmZnLiiScCMHLkyN1e60/dfdhhh7Fz50527NjBhRdeyHPPPQfA008/zQUXXNDi76cj0xpFuMoLwBVfPyMbnGUDNWW2v0KpTiacX/6tpWEa7srKSuLi4uqS+lVVVTW6v8vlqnvucrl2S+EdLAV4nz596NatGx9++CGLFi3aI0vs3k5rFOEqL7RZYwM/ZEma70mpttS/f3+WLFkCwOzZs6M6hj9198KFC8nMzCQz0/74u+iii5g2bRpnnHEGbre7ZQrcSWigCFdF0e4jnqC+dqGT7pRqE9dddx2PPPIIY8aMobCwMKpjJCUlMWbMGGbMmMFTTz1Vt37q1KmUlZVps1MQmmY8XE8eDQkpcN479et+WADPnwwX/Bv6Hdzy51QqBjpzmvFJkybtcX9sv8WLF3PNNdfw6aefxqBkbUvTjLeW8oLdRzyBZpBVqpO48847Oe200+puZqR2p53Z4QrV9KSzs5XqEALvdx3ohhtu4IYbbmjbwnQgWqMIR7A8T6A1CtVpdbYmaVUvmn9bDRThCJbnCervcqeT7lQnkpSURFFRkQaLTsgYQ1FRUd29vcOlTU/hCJbnCSAuAeKStUahOpXevXuzadMmTWHRSSUlJdG7d++IXqOBIhzlQWZl+2kGWdXJxMfHM2DAgFgXQ7Uj2vQUjsaankAzyCqlOj0NFOGoa3rK2XObZpBVSnVyGijCUV4IrjhIytpzW2KGzsxWSnVqGijC4Z9s1yCZGKA1CqVUp6eBIhwVRcE7skE7s5VSnZ4GinCUF4YIFNqZrZTq3DRQhKOiMPiIJ7A1Cm811FYF366UUh2cBopwhKpR1M3O1uYnpVTnFNNAISJTRORbEVknIntk5BKR6SJSICLLnMdFbV7IujxPjdUosuxfbX5SSnVSMZuZLSJu4CHgaGAT8JWIzDHGrGmw66vGmN+0eQH9Kors30abnvx3udMahVKqc4pljWI8sM4Ys94YUwO8ApwUw/IEV+7kuwk16gmgakebFEcppdpaLANFL+DngOebnHUNnSYiK0Rktoj0aZuiBShvJCGgn94OVSnVybX3zux3gf7GmFHAB8CzwXYSkUtEZLGILG7xjJdNNT35O7O1j0Ip1UnFMlBsBgJrCL2ddXWMMUXGmGrn6ZPA2GAHMsY8bowZZ4wZl5fXyC//aNU1PQXJ8wR6lzulVKcXy0DxFTBYRAaISAJwNjAncAcR6RHwdCrwTRuWzwqV5wkgIRXErTUKpVSnFbNRT8YYj4j8BngfcANPG2NWi8itwGJjzBzgKhGZCniAYmB6mxfUP9kuWJ4nsOt1drZSqhOL6Y2LjDFzgbkN1t0UsPx74PdtXa7dhJps56cZZJVSnVh778yOvfJCSGmkf8JPM8gqpToxDRRNqShsfGisn2aQVUp1YhoomhJO05PWKJRSnZgGilD8eZ4am0Phl5SpfRRKqU5LA0Uo/sl24XRma41CKdVJaaAIpS59Rzg1il3g87V+mZRSqo1poAjFPyu7yaanDMBo85NSqlPSQBFKXdNTGKOeQJuflFKdkgaKUOqansKYRwFao1BKdUoaKEIpLwid58mvvWaQ3fET/Lwo1qVQSnVwGihCaSrPk197zSD78V3w6q9iXQqlVAengSKU8qKmRzxBwO1Q21mNYudWKN8OPm+sS6KU6sA0UIRSXtB0nieob5pqb30U5QVgfFBZEuuSKKU6MA0UoYST5wkgMd3+bW81Cv/w3vIWvuufUmqvooEilHCbntzxEJ/avgKFMRoolFItQgNFYzzVUF3a9GQ7v/Z286LKEvB57LIGCqVUM2igaEzdZLsw+iig/WWQDQwO5UWxK4dSqsPTQNGYusl2YfRRQPvLILtboNAahVIqehooGlPhBIpwm57aWwbZsu31yxoolFLNoIGiMeFmjvVrb3e58weH1DwNFEqpZtFA0ZiIA0U7q1GUF4C4IHe/+veilFJR0EDRmIrC8PI8+fk7s41p1WKFrWy7bTZL76Y1CqVUs4QVKETk7yIyvLUL0674Z2U3lefJLzEDfLXgqWrdcoWrvMA2O6XmaY1CKdUs4dYovgEeF5EvRWSGiGS2ZqHahfKi8Ec8Qfu7J0V5AaTl2aaz6lI7L0QppaIQVqAwxjxpjJkInAf0B1aIyEsiMrk1CxdTFYXh5Xnya28ZZMu2Q2rX+mBXoXMplFLRCbuPQkTcwBDnUQgsB64VkVdaqWyxVV4Qfkc2tM8aRWpe/fBe7adQSkUpLpydRORe4ERgPvBXY4z/bjh3ici3rVW4mIq26am6HQSKmnKorXCanpz3oIFCKRWlsAIFsAL4ozGmPMi28S1YnvbBUxNZnidoX3e580+2S+1aXyvSDm2lVJTCDRTLgf1k9xFApcBGY0w7+GZsYRVh3is7UHtqegqcbKc1CqVUM4UbKB4GDsDWLAQYAawGMkXkMmPMf1qpfLERaZ4nCLjLXTvozPYHhbQ8e68Md6IGCqVU1MLtzN4CjDHGjDPGjAXGAOuBo4H/a63CxUykeZ4A4lPsBL32UKMIbHoS0bkUSqlmCTdQ7GuMWe1/YoxZAwwxxqxvnWLFWKTpO8B+IbeXDLJ1TU+59X81UCilohRu09MaEXkE8A+FPctZlwjUtkrJYsn/pRrJPApoPxlkywts0IpLtM9Tc7XpSSkVtXBrFOcD64CZzmM9MB0bJDrfpLtI8zz5tZcMsv7Jdn7a9KSUaoYmaxTORLu5xpjJwN+D7FIW7clFZArwD8ANPGmMubPB9kTgOWAsUAScZYzZEO35wlbuzMp2RZgzsb1kkPVPtvPz1yiMCT93lVJKOZr8JjTGeAFfS+d3cgLQQ8CxwDDgHBEZ1mC3XwMlxph9gHuBu1qyDI0qL4ysI9uvvdwO1Z/nyS81DzyVdiKeUkpFKNw+ijJgpYh8ANR92xhjrmrGuccD6/wd4k4qkJOANQH7nATc7CzPBh4UETGmlXN5VxRG1pHtl9hOOrPLtsOAw+ufB86lSEyLTZmUUh1WuIHiTefRknoBPwc83wRMaGwfY4xHREqBHGyuqdZTXgg9R0f+uvZQo/DUQNWOBk1P/kBRCNkDYlIspVTHFVagMMY8KyLJQF9jTLvL7SQilwCXAPTt27f5B2xO01NNGXg94A43Brcw/xyQtAZ9FKAjn5RSUQn3xkUnAsuA95zno0VkTjPPvRnoE/C8t7Mu6D4iEgdkYju1d2OMedyZDDguLy+C2dTB+PM8RTIr288/OzuWzU+Bk+38/EGvQkc+KaUiF+6wnpuxfQo7AIwxy4CBzTz3V8BgERkgIgnA2UDD4DMHOzQX4HTgw9bvn3DiUCR5nvzqMsjGMFAE5nny0xqFUqoZwm0fqTXGlDZICuhrzomdPoffAO9jh8c+bYxZLSK3AouNMXOAp4DnRWQdUIwNJq3L/2UaTdNTe8ggG5jnyS8+GRLSdS6FUioq4QaK1SLyP4BbRAYDVwGfN/fkxpi5wNwG624KWK4CzmjueSJSEUVCQL/2kEE2WNMT6OxspVTUwm16uhIYDlQDLwM7sTO0O59yf9NTNJ3Z7SCDbHkBxCVDQuru61PzNFAopaIS7qinCuBG59G51TU9NaOPItZNT2l5e87ATs2DHRtjUyalVIcW7q1Q9wWuA/oHvsYYc0TrFCuGKgpB3JHneYL20ZndMM+TX2oubF7c9uVRSnV44fZRvA48CjwJeFuvOO1AuTMrO9I8T9B+OrMz++y5PjXXjujy+aJ7b0qpvVa4gcJjjHmkVUvSXlQURTfiCcDltqOLYt1H0euAPden5oHPY2dtp2S3ebGUUh1XuD8t3xWRy0Wkh4hk+x+tWrJYKS+Ibg6FXywzyPp8To0oWNNTQBoPpZSKQLg1Cv+kt98GrDM0f9Jd+xNtnie/pEz7qz0WKovBeIMP7Q2cdJe3b9uWSynVoYU76mnvySRXEWWeJ7/EjNh1ZgebbOcXmEFWKaUiELLpSUR+F7B8RoNtf22tQsWMp8Y2G0Uzh8IvlhlkG5tsBxoolFJRa6qPIjBlxu8bbJvSwmWJvYpmTLbzi+XtUIPlefJLzgZE+yiUUhFrKlBII8vBnnd8/vQdzWl6imVndl3TU5AahTvOjnbSDLJKqQg1FShMI8vBnnd8db/Im1mjqN5p70/d1sq2gyuu8cmCKZrvSSkVuaY6s/cXkZ3Y2kOys4zzPKlVSxYLdXmemnFPi8QMO1+htmLPfEuhrHoDfvoSjvu/6M9dXmCDQWMT6lLztOlJKRWxkDUKY4zbGJNhjEk3xsQ5y/7n8W1VyDZT1/TUnHkU/nxPEfZTLH0eFj9l744XLX+ep8ZoBlmlVBQ0l0Og8oLo8zz5JUWRxsMYyF9payKlP0V/7sbyPPlpBlmlVBQ0UARqTp4nv2gyyO7aWl+bKVof/bnLC0M3m6XmQWUJeGujP4dSaq+jgSJQc/I8+SVGkUF264r65eIfojuvMVC+vemmJ6gfBqyUUmHQQBGouXmeILoaRb4TKOKSoSjKQFG9CzxVTTc9gTY/KaUiooEiUHkz03dAdIFi63LIHgS5g6OvUYSabOdXl+9JRz4ppcIXblLAvUNFE2384YimMzt/BfQaa5e3Lo/uvKHyPPlpBlmlVBS0RuGorSnnE6kmPyG5eQeKSwJ3Qvh9FJUlsOMn6D7K1ipKNkbX2Rwqz5NfYAZZpZQKkwYKx6oNK7mie1eeKNrQvAOJ2El34dYo8lfavz1GQc4gmyZ8RxRDZMNpekrKsjO3NVAopSKggcIxcuCBxFf1559VmzDNTb8RSQZZ/4in7vvbGgVE16EdTvoREZ1LoZSKmAYKR1ycm1E5x1DBdr7Y9HXzDpaUEf7M7PwVkN7D9i3kOIEimg7tsu02Q6y7iQnzqbnaR6GUiogGigCXjTsV44vjsaWvN+9AkdYouo+yyyk5ttkq2hpFOB3xWqNQSkVIA0WA8f16klQ7iuUlC6htzuxlfwbZptRWQuF3tn8CbNNQ9kAojmJ2dnlB8PTiDaXkaqpxpVRENFAEEBGO6H0sXinn7W/nR3+gcDuzt62xndf+GgXY5qdom57CrlFooFBKhU8DRQOXTzgOnyeV51e9Gf1Bwr3LXb4zZ6JHQKDIHmRHPXlqIjtnU3me/FJzoaYMaioiO75Saq+lgaKB/jkZ5DKBHyu+ojTaO9UlZUJtedPzIbausLmhsvrVr8sZBMYHOzaGf77aKqguDT3Zzs8fTLT5SSkVJg0UQZw8eCqIh6eWvRPdAcK9J0X+Cug+0vZN+GUPtH8j6dCuGxobRh+F5ntSSkVIA0UQF4w7FF9NV+asmxPdARKdNB7VIWokXg9sW717sxPUz6WIpJ8inMl2fprGQykVIQ0UQWSlJDAw6TCKvN+yIZobCYWTGLBonc322r1BoEjJtq+PZORTXZ6ncGoUmsZDKRUZDRSNOG/EKQA8/FUUcyrCaXrypxZvWKMQsbWKSJqe6vI8hdmZDRoolFJh00DRiJNHjUCqBrFg878jT+kRTgbZrcvBnQi5++65LdIhspE0PSWkQnyKNj0ppcIWk0AhItki8oGIfO/87dLIfl4RWeY8ouwwiE6828WYLkdSxTa+2Lw0shcnhXGXu/wV0G1Y8JQb2YOgdBN4qsM7X3kBJKRBQkp4+2saD6VUBGJVo7gBmG+MGQzMd54HU2mMGe08prZd8axLx52C8cXxeKQpPRKbqFEYs3vqjoayB9ohsiUbwjtfuJPt/DSNh1IqArEKFCcBzzrLzwInx6gcIf1iQC8Sa0axrPijyFJ6JGYA0nigKP0Zqnbs2T/hlxNhFtlw8zz5aaBQSkUgVoGimzFmq7OcD3RrZL8kEVksIv8VkZPbpmj1RITJvWxKj3e+iyClh8sFiemNd2YHphYPxj+XItyRT+HmefLTpielVARaLVCIyDwRWRXkcVLgfsb2FDfWW9zPGDMO+B/gPhEZ1Mi5LnECyuKCgpb9pXzFhOPwedJ4buUbkb0wVAbZ/BUgLug2PPj2lGxI7hJ+h3a0TU/Nve+GUmqv0Gr3zDbGHNXYNhHZJiI9jDFbRaQHsL2RY2x2/q4XkY+AMcAe357GmMeBxwHGjRvXot9+A/IyyDHj+bHiE0qrSsn0d1Q3JVQG2a0rIGdw6M7ncIfI+rxQURRZoEjJBV+tDWTJWeG/Tim1V4pV09Mc4Hxn+Xxgj1wZItJFRBKd5VxgIrCmzUoY4CQnpcczy98O/0WhMsjmr2i8f8Iv3HTjFUWAibDpyZ/vqSj81yil9lqxChR3AkeLyPfAUc5zRGSciDzp7DMUWCwiy4EFwJ3GmJgEil+POxRfdVfeXvdu+C9qrOmpvAh2bm58xJNfjjNEtrYq9H6RTLbz00l3SqkItFrTUyjGmCLgyCDrFwMXOcufAyPbuGhBZaUmMCDpMDZ6ZrNhx0/0z+rb9IuSMmB7kLgWLLV4MNmDAAMlP0LXoY3vF8lkOz9NDKiUioDOzA7TuSNOwRjh4SWvhfeCxmoUdSOemqpRhDnyKZI8T34aKJRSEdBAEaZTRw1Hqgbx0aYwU3okZtjO7Ib75q+AjN52ZFMo2WHOpYim6Sklx/7VIbJKqTBooAhTYpyb0V2OoJLt/HdLGCk9kjLt7Oqast3X569sutkJ7GiklJymh8iWF4A7oT5tSDjiEiApS2sUSqmwaKCIwKXjTg4/pUewDLI15VD4fdPNTn7ZA5uuUfhnZQfe/CgcqbkaKJRSYdFAEYGJA3uTUDOKr4vCSOkRLIPsttWACa9GAbb5qak+ikgn2/ml5mnTk1IqLBooIiAiHN5zipPSY17onYNlkN3qjHgKt0aRM8gOpa2paHyfSPM8+WkaD6VUmDRQRMif0uP5VW+G3jExyF3u8lfY1ByZvcM7mT/nU6gsspHmefLTxIBKqTBpoIjQPl0zyfb9gvXlX7KqYHXjOwa7Hao/tXi4/Qk5Tdw/25hm1Cjy7Mxsnzfy1yql9ioaKKJw1dgZ+LypzJx/Ez7jC75Twz4Kb62dgBdu/wTU1yga69CuKgVvTfSBAgMVxZG/Vim1V9FAEYXTxwxmn7iz2Vb9HU98/XLwnRrevKjgW/ul3lhq8WCSMm0Cv8ZqFNFMtvPTNB5KqTBpoIiCiPDISRdjKgfyyIr7KaoMklwvPsneE9vfmZ3vzMiOpEYBtvmpqJGRT9Gk7/BL0UChlAqPBooo9cxK4ddDr8NDFTP/85fgOwWm8di6AuJTIGefyE6UPajxGkU0s7L9NI2HUipMGiia4epDDyHHczTLdnzAhz9+uecOSZn1E+7yV9gbFbnckZ0kZyDs2mon6zXUrKYnTTWulAqPBopmcLmEh467HlObxe8/uZlaX4NJeEnOPSl8Ppu6I9z5E4Hqbov6457bygsAqc/dFInkLvYue1qjUEo1QQNFM43omccvu8+ggk3c/NFju2/0Nz3t2GD7KiLtn4D65IDBmp/KttsgEWktBex9vVM0jYdSqmkaKFrAX395Nkk1I5iz8Rm+L/65foM/g2y4qcWDyQmRRTbayXZ+msZDKRUGDRQtIDHezV8O/xMGwxX/vrl+g79Gkb8CxA1dh0Vx8HRI7Rq8RhHtZDs/TQyolAqDBooWcsy+wxiZehpbPYt4aslcuzIpw3Zmb10BeUPskNloNDZENtqEgH4aKJRSYdBA0YIeOfEaXJ6u3L/8b+yoLLc1Ck8lbFkK3ZtxV9fGsshq05NSqg1ooGhBWckp/GbU7/C5C/nNv/5ubw4EdghqNB3ZftkDoCwfqgNuglRTYW+K1NwaRfVO8FRHfwylorV1Baz/ONalUGHQQNHCLj7wl3R3/4Jlu95k0c6AL/ZoOrL96pIDBtQqmjMr269u0p3WKlQbq9wBL5wKL57R9M25VMxpoGgFD0/5M0I8f/jxPerumN3cpifYvUO7OZPt/HR2toqVj+6wP1Dc8fCv/93z3vKqXdFA0QoG5/bihD4Xsi1hC++npkBWP3sP7GgFyyKrNQrVUeWvgkWPw7gL4aibYf0CWPVGrEulQtBA0Upum3wxmb7u/F92Ftu7DGnewRLTIK377k1Pzcnz5KcZZFVbMwbmXmf77474ow0WPcfA+3+wzVGqXdJA0Uri3HHcPv5aCt1uLqzayqYdpU2/KJTsgQ36KFogUGgGWdXWVrwGP31haxIp2TarwAn32s/gh7fHunSqERooWtGkYVO4MuMQNibt4OTZl1KwK0hiv3DlDGzQ9FRob7ca7dwMsJP53ImxDxQf3QUvnKbt1J1dVSn854/QayyM+VX9+p5j4MCL4asnYfPS2JVPNUoDRWsS4eJTH+XcQf9LdfxqTnjtYkoqKqI7VvYgW4vwZ6Mt217fdNSM8sV8LkVFMXx2H6ybBxsWxq4cqvV9dKf9UXLc3TbXWKAjboS0bvDPmXp73nZIA0UbuOGQ6Zza92oq4lZywquXsKuqKvKDNBwi29zJdn6puVARw0Dx1VNQWwEJ6fDfR2JXDtUijDGs3lKKx9vgFsHbVsOXj8HY6dDrgD1fmJQJU+6ArcttzUK1Kxoo2sgtky/iuB5XsNO1nONfvZTK2ggnuTUcItvcPE9+qXmxa3qqrYJFj8E+R8OES+HbucHTqasOwRjDnf9ey/H3L+Ssx//LppIK/wb413U2GBx5U+MHGH4KDDoS5t8GO7e2TaFVWDRQtKG7jpnBEbmXUMJSTnhlBtWemvBfnD3A/vXnfGpunie/SJueti6H1W81/7wAy1+2QWriVXDgRbZjc9HjLXNs1aaMMdzy7hoe+2Q9Rw3txrf5uzjuH5/y3qqtsPJ1+OlzOOrPtgO7MSJw/N323vLv/77tCq+apIGijf3j+Cv5ReaFbPctZuprl1PjDTNYJKRCeg/b9OSthcrilmt6Ki8IryO5aie8dBa8fkF96vRo+XzwxYPQYzT0PxQyethflF+/ANW7mnds1aZ8PsONb69i1ucbuOiQATxx3lj+ddUh9M9N5boXFrJrzg34ehwAY85r+mDZA+Gw39ofI+vmtX7hVVg0UMTAYyfNZEzqeWyp/ZLTZl+Jx+cJ74X++2f7b1/aIjWKXPBU2bxRTfnoDtiVb0dLvf+H5o1S+nYuFK2ztQkRu27CZTb31LKXoj+ualNen+F3b6zgpS9/4vJJg7jx+KGICP1yUpk942Ce6jeP1NpiZu46l3WFYQ7kmHgV5Ay2zVW1la37BjqL2ir7I3L72lY5fFyrHFWFJCI8e9p1nPmKh7VVL3Hmm1fx2qn3E+dq4p8jZyCsndsyk+38AtN4JKY3vt/W5fDlo3aCVNehdtLUt3NhyPHRnffz+yGrLww9qX5d77HQe7w9z4EX7zkyRoVn02LYuRmGTq0Pwq3A4/Xxv68v551lW7jmqH256sh9kIDzJRStZcL219m8z5ks3NCPDx5YyK0nDef0sb13228PcYlw/N/huanw6T12RFRrqKmAuKSIP2cer4/tu6rZWlrF1tJK8kur2Fayi8riLXh3biXZu4tUU06qr4xUU0aKr5w0Y5dTfeWkmnLSfGUgQkVcJuXxmZQnZFKZmElNUga1yRl4k9PxpWQgyamQlISnqghvRQG+ykK8VcV4q4sxNaV4a3biq92F11tFrUBGUg8uuPTLFr9UGihiRER45awbOOVFD9+Xv8a0d67hxZPuwx3qtqbZg+wIpaJ19nmLND0FpPHwpwppyOeFd2fa264eeRMkpMGiJ+A/f7Id0XEJkZ3zpy/h5y/h2P8Dd4OP4EEzYPaF8P1/YL8pEb+dvd7XL9h/K18t9Jtoh6J2i+KGWU2o9fq4+pWvmbsyn99N2Y/LJ+2z+w7GwNzfQlIGvU67k397Upj5yjJ+O3sFn60r5PZTRpKWGOLrZ+DhMOosWHgvjDwD8vZtfqF3bYONC2HDZ7DxMyhYC+LGpObiTc6hJjGXyoRsyuK6UOruQolkst2XziZvEkWV5dSWbcZXlU+cp4hU105SXLtIdpUR76qgl6uaCpeLKhEq44RSEbaJ2OcuNxXipsLlpkpcVItQJYnUio9aKcMrZcDm+nJWOo/iJt6PG0gGkhOBRMQIXTxZXND8K7WHmAQKETkDuBkYCow3xixuZL8pwD+wl+RJY8ydbVbINuB2CbPPuZETn/eweuebTJtzFY8ddwcZCRnBX+AfIvuz84uhpZqeIHSH9pJn7D01Tn2iPmfVL/8CL55uhzL+4vLIzvn5/ZDcBcZM23Pb0KmQ3hO+fEQDRSR8PvjwVvvFOnAyDD3BznR+9BD773P4DTYVTAuo9ni54sWvmffNNv54/FAuOjTID4yVs+2X8gn3Qko23YAXLprAwwvWce+871j28w4eOOcARvbObPxEx9wO370H/7oWzn837NqRMYbiijI2b15LyYZPqN66GNmxGl/1dspcLna4E/kxrisbcydQY2owUoWRanyejXh866lx+ahwCRUiVPlrG/FAl8bOmOI8IFHiSHInkhiXSHJcMolxKSTFJ5PkTiIrLolEdyJJ7iSSnGX/I8GdQIIrAZfXQFUVVFUiVRW4qspxVVciCV0wyV0huRuk9CAhIY0EVzyJcfEkuhNIcCcQ53aTHir4NkOsahSrgFOBxxrbQUTcwEPA0cAm4CsRmWOMWdM2RWwbCXEu3j73Jk54DlbueItJLx/LtQf8jnNHTN2zeu7/xf/Tf+3flm56CqZsO8y7FQYcZn/Z+e1zFAw6Aj6+E/Y/O/RolkCF62Dtv+Cw62wHfUPueBh/Mcy/BbataZVfw51OTQW8dSl8MwfGXgDH/c1ex2GnwPyb4fMHYOUbMOWvMOzkyJujin+ENe9A2TZqswZx79eG5RuTuW3qRH518IA996/aaWdg9xwDB5xft9rtEq48cjATBuZw9Stfc+ojn3HqmN4kJ+xeixaxX/ZeU8XgrucydsuTzH32Otak7EeFp4xKXxnV3jKqfeV4faUYsxNDGV6pwCM11LhqMQ3fYhZA4P+XCtz4SHClkODKItmdTEp8KqnxqfRISCUjPokubiFThDSfh2RXPCmpXUlO60FyWneSE9JIjksmJS6F5PhkkuNsMAjZItCBiYlh2gQR+Qi4LliNQkR+AdxsjPml8/z3AMaYO0Idc9y4cWbx4qAVlHatqtbLHfM/4I2N9yFJm+mdeAAP/fI2BnbpW79TTQX8tYe9/7Y7Hm7Mb34bdG0V/KUbHPEn++Xd0BsXw5q34bLPIXfw7tu2rYFHJ8L4S+DYu8I737tXw7KX4ZpVjTedVRTDPcNg1Jkw9f6I3k6H4/PZ4aDRpmLZlQ8vnwNbvra1vIMu3/Mz8fMi+6s8f6UN7sfdXV87bUzxj/bfffXbsHUZACYuCfEETBZNzLCfiZzB9m/uYMjdF5Y+ZydPXjQfeo/FGMOu2l2UVJVQUlVCcVUxm3cW8trXa9lYsh3jqgB3BbjKoW65AsQXrGQAxBtDhtdHhq/+ke78TTZxmMSeSNpAEnNHkNltGD2zcuiZnkVWUiapTkBwifaBBRKRJcaYccG2tec+il7AzwHPNwETYlSWVpcU7+aWKVO4qHgiV/7rQdZVzuakt0/hhL7TuXXSDOLd8ZCQAhm9bEdlaq+W6aiMT7L/4YM1Pa3/CFa+Bof9bs8gAfbX/tjptvnpwIuC7xOobLsNEqPPCd2/kpJtg8SKV+uTx3VGPy+Ct2bAzi2w/1l21FfXCDIN56+yw5UrS+Dsl2DIccH36zMeLv4IFj9lm6MePggmzoRDr4X45Pr9ggQHeh4AR9/GzoHHc/E7+Wza+AN3HZ7E2MwCigpWU1yyjqLNn1L047sUu9wUud0Uu12UDB5J8dK/UPJ5CTuqduAxwUf2JWUnkZmYSWZiJlmJXZ2/Wbv9zazYQeYHN5NpXGSk5JKe2pWk9B72M5TW1ab+8P9N7Wr70nQgRItqtRqFiMwDugfZdKMx5h1nn49ovEZxOjDFGHOR8/xXwARjzG+C7HsJcAlA3759x27cuLHF3kes/GvNGm7+7HaqElaSbPpyy8F/5th9x8OsE2DDp/Y/8CULWuZk/xhtE7Wd/lT9Ok81PHKw7ci+/Ivdv1AClRXAAwdAv4Phf14NfZ75t8Gnf4ffLIbcfULvu/0b+4V25E1w6P9G9HbaPU+1zXv02X2Q0Rv6HwKr37TDlAdOtrWCfY4K/WX33X9g9gU2yP/PK9Bj//DOvSvfDkJY+Rpk9aPi0BupKNhA0ndzSCxeRbHbzZrUQSxKGc4iVy82ebxUeIqplZ244neRl1lLla+UCk/woa7p7kSyXYlkZ/ajS3IeXZK62Edil/rlpC5kJ2aTlZRFclwjn6uGjGnVEVwqRjUKY8xRzTzEZqBPwPPe7DY0YLdzPQ48DrbpqZnnbReOHzaMY/Z9npvnv8KcTQ/z288v4omlU3gqs4/tU2uJEU9+wdJ4fPYPO7pq2huNBwmAtDz7RT7vz/DDAhg0Ofh+1WW25jHk+KaDBNghuAMnwaIn4eCrbFNbZ5C/0tYitq2ynfm/vAOSMmzH7ZJn7DV66QzI2QcmzID9z9mzE/rLx+C9G6DbCBucM3o2ejqvz/BzyS5W529mbeEm1pdsYUvFYOh2En29S6hadCMFcW7y0xLYlelv5qwClthHIiRIOl3isuiR2pX+XbqRk5xDdlI2OUk55CTn1P3NTsomwR3hCLhwaZCIqfbc9PQVMFhEBmADxNnA/8S2SG0rPs7NX355LpfumMIVc//Cd1XvcWx5An9NSWZSSl7LzZZMzYOSDfXPi36AT+6G4afaX7ZNmTADFj8N798IMz61qTga+voFqNphv/TDNeEyePks25E68vTwX2eMHb67fbVtYskO0uHa1rwe+PwfsOAOO+LrnFdgv2Prt6fm2D6iiVfb9/vfh+1clfm3wQG/sh38Gb1tgPjqCbz7HsfWI+/n56IqfvxhET+XbiW/fBvbK7dTUl1AaW0R5Z4iPLID3OWINPj9lOxisyuHHEkgK7kHY9L70DOtK30zu9ErvRt5KXnkpeS17pe/6jBi0pktIqcAD2CHIewAlhljfikiPbHDYI9z9jsOuA87PPZpY8xfmjp2R+3MDsdbaz7jni9+z464ErpWJ9Mj6ywuPuB0DhvUN/QEpqbMuQq+/Tf89nv7JfvCqfDzV/Cbr2xqjXCsfgtenw4n/sP2WwTyeuCBMTYFya//E365fD54cKxtc74ozHQO3lp7D+alzwJig9YB59u0EOG+l5ZWuA7engGbvrKjjo6/xwaGEGo8Xr5d/m9Klj2Fp+hLtrld/BCfSYWrkmXubmx0JyHxOxFXkBQw3jTiySI9LofspFy6p3ajX2Z39snpxZDcXnRP60aXxC6ddoSOik6opqeYjnpqDZ05UADUbF3BWy8dy1NdBrDVXYrxuUmoHsWhPY7nkvHHMLxHo4O9Gzf/Njv+/k+FtjNz9gV2MtyES8M/hjHwzLG2uerKpbY5xW/lbHjj13DWi3Z8fyS+fAz+/TtnBE3Qz3C9yhJ47Xz48WPbHDbu17ZPZOmz4Iqzo7MOuabtOsd9PtuU9MFN9bONR5yGD0NxVTHbKrbxc+kWVm77iR+KN7O5bCtFVQWUe4vwSgni2v2+DG5jSDOpxCX0JjMhl9ykrnRL6UrvjB70y+zBwC69GNilB4lxiW3z/lSnooGiM/H5YOHfYfQ0lpQX8fDil1lSNA+vVOCrzSTD8wtOHDiVaWPH0jcnJbxj/vdReO96uGoZPD0F0rvDxR8Gb0IKZfMSeOIIOORamykUbAB5/HCoKYcrvop8NEr1LjtUdvAxu3e2N1S8Hl480zahTb0fRge0Uhb/aDuPV7xq05QcfCUcdFnolCVR8hkfJVUlbNu+kvyPbmNb4Vry8waR32s022p3sqUsn+0V2/Ga2t1eZ4wbU5tJonQhMz6Pbind6JvZg/1y+zCqRz8GZPUiOylbh3SqVqOBopOr9lbz7vcf8NzK1/mx/GsQg6d8ID3dh3PCoGMY3DWbvtkp9MtOJTMlSKew/xf/wMl2SOzFHwa/uUw43rzEDq/8zVfQpR+s/9jm7AnWJBWu9/5g71sxc2XwjtuNn8Mr5wLG1lr6Twx+nO3f2OGha/9pm7P8tY5gcxiMsTWU4vV1D2/ReorKt7LNV2UfpoZtvlryxcs2fGxzGba7oLZBM6Bb4kiUbDw1GZSXp+GrzUS8WfTL7MH+Pfozvs8A9u/ZmwG5acS7NRCo2NBAsRfJL8/n+VWzeXvd2+z0bMP43PiqeuGt7IO3si/JZgD9MnvTLzuVPtkp9MtJYWT114yY76SAHn+JndkbrdLN8MBY21F7xjP2XthbV9gv+WgnlRX/CPePsV/sR/5p923LX4E5V9oEg//zWtMTycDWfObfBusX2Hkph1xDucvNtsJv2F76IwW7NrOtspDtppbtcW62u91si7NzBLwNgkC8EbJ9cWSbeLJ8CWSZRDJMEsm+VBZWHsSyoh4YbyrJ8fEc0C+LA/tnM35ANmP6dNljRrJSsaSBYi/kMz6WbFvC/I0fszR/Get2rqXWZ++qF2cycdX0o3xnL2or+jCwyvBBwp8oMFlclfM4fXp2Z99u6QzpnsF+3dPJTUuIrLN8wV/h47vghPvsPZCP+KPtTG6OV861NYdr19jhuj4fLPgLfHq3TS9y5nN2NFEAr89LSXUJRZVFFFQWUFhZSGFlIUWVReSXb2db0TpKd22kiFrKgzSJxXvdxHlS8XgyqarNxePJwngy8NVmYmozMZ4sxJdCnMuN2yXEuYQ4t+B2uUhwC8N6ZtQFhhG9MrW2oNo1DRSKWl8t35d8z4qCFfZRuIKNO+3ERBcuBnp8xMUNpMA7gB27kikrT8PUZuDzZJKVkMuQbtns1z2d/bqn0z8nlS6p8XRJSSAzOZ6k+Aa/jKvLbK2iLB/iU226jig7kI0x7KzyULJ6Hnn/OodPhv2W73LH0X/d/cTvWMr63AP4vsfBVPrKqPDspNxTWvcoqy3FR5A0EL5EfLXp+DzpGG8aXTwGryeLStOTzMTudE3Jo1dad3pmZtI1I4nuGUl0y0ikW0YSmSnxxLtcdYHB5dLx/apz0EChgiqpKmFl4UqWFyxndeFqtpRvYVv5tqCzbl0mFW9NBt7aDIwnFWPiwReHMfHESTzJ8YkkxyWRmpBEWkIyI2u/4dDCt1iffQhre5yAy+XD5fIiLh8iPhAPIj5EvBi8VHqqKK7cyY6qXeyq2UWFp4Jqbzm1phJcVeCq3nMugMMYwXiTMd5U+/CkBiynY7zpZMZn0zUljx5pXemVmUn3TBsAumcm0c35GzLttVKdnAYKFZGymjK2V2wnv8KO0NlWvo1tFdvYVr6NTbvyKakqocZXQ62vBo+vFh9h3qEvBOOLx/iSwJuImxQS3cmkxKWRnpBGVlI6OSkZ9KvYQN/175JJHHHjriG+3zGkxWeS5E7DGBden8HjM3h9PjxeQ3ycix6ZSeSlJRKnzT5KhdRRkwKqGElLSCMtIY2BWY3cyKgBr89Lja+GGm8N1d5qqr3V1HhrqPJUAW6Mz4XPuMHnxudz4TOC1+fG53Xh9Qm1HiEhzk3PLPvLPjGukU5eTzV81MNOWus5uqXerlKqCRooVLO5XW6SXcnhJ3iLVlyizSarlGpTWh9XSikVkgYKpZRSIWmgUEopFZIGCqWUUiFpoFBKKRWSBgqllFIhaaBQSikVkgYKpZRSIXW6FB4iUgBsbMYhcoHCFipOa+tIZYWOVd6OVFboWOXtSGWFjlXe5pS1nzEmL9iGThcomktEFjeW76S96UhlhY5V3o5UVuhY5e1IZYWOVd7WKqs2PSmllApJA4VSSqmQNFDs6fFYFyACHams0LHK25HKCh2rvB2prNCxytsqZdU+CqWUUiFpjUIppVRIGiiUUkqFtNcEChGZIiLfisg6EbkhyPZrRWSNiKwQkfki0i9gm1dEljmPOe2kvNNFpCCgXBcFbDtfRL53Hue3g7LeG1DO70RkR8C2Nr22IvK0iGwXkVWNbBcRud95LytE5ICAbW16XcMs77lOOVeKyOcisn/Atg3O+mUi0ur3Bw6jrJNEpDTg3/umgG0hP0MxKOtvA8q5yvmcZjvb2vS6OufsIyILnO+o1SJydZB9Wu+za4zp9A/ADfwADAQSgOXAsAb7TAZSnOXLgFcDtpW1w/JOBx4M8tpsYL3zt4uz3CWWZW2w/5XA0zG8tocBBwCrGtl+HPBvQICDgC9jcV0jKO/B/nIAx/rL6zzfAOS2o2s7Cfhncz9DbVHWBvueCHwYq+vqnLMHcICznA58F+Q7odU+u3tLjWI8sM4Ys94YUwO8ApwUuIMxZoExpsJ5+l+gdxuXMVCT5Q3hl8AHxphiY0wJ8AEwpZXKCZGX9Rzg5VYsT0jGmE+A4hC7nAQ8Z6z/Alki0oO2v65hldcY87lTHojx5zaMa9uY5nzeoxJhWWP6mQUwxmw1xix1lncB3wC9GuzWap/dvSVQ9AJ+Dni+iT0vcqBfYyOzX5KILBaR/4rIya1QvobCLe9pThVztoj0ifC1LSXs8znNeQOADwNWt/W1bUpj76etr2s0Gn5uDfAfEVkiIpfEqEwN/UJElovIv0VkuLOu3V5bEUnBfqm+EbA6ptdVRPoDY4AvG2xqtc9uXMSl7OREZBowDjg8YHU/Y8xmERkIfCgiK40xP8SmhHXeBV42xlSLyKXAs8ARMS5TU84GZhtjvAHr2uO17XBEZDI2UBwSsPoQ59p2BT4QkbXOL+lYWYr99y4TkeOAt4HBMSxPOE4EPjPGBNY+YnZdRSQNG7RmGmN2tsU5Ye+pUWwG+gQ87+2s242IHAXcCEw1xlT71xtjNjt/1wMfYaN5a2qyvMaYooAyPgmMDfe1LSyS851Ngyp8DK5tUxp7P219XcMmIqOwn4GTjDFF/vUB13Y78Ba2iSdmjDE7jTFlzvJcIF5EcmnH15bQn9k2va4iEo8NEi8aY94MskvrfXbbskMmVg9szWk9ttnD31k2vME+Y7AdaoMbrO8CJDrLucD3tH5HWzjl7RGwfArwX1PfcfWjU+4uznJ2LMvq7DcE2wkosby2zrn603iH6/Hs3iG4KBbXNYLy9gXWAQc3WJ8KpAcsfw5MiXFZu/v//bFfrj851zmsz1BbltXZnontx0htB9dVgOeA+0Ls02qf3b2i6ckY4xGR3wDvY0dYPG2MWS0itwKLjTFzgL8BacDrIgLwkzFmKjAUeExEfNga2J3GmDXtoLxXichUwIP9ME93XlssIrcBXzmHu9XsXm2ORVnB/jJ7xTifXEebX1sReRk7+iZXRDYBfwbinffyKDAXO3pkHVABXOBsa9PrGkF5bwJygIedz63H2Oyh3YC3nHVxwEvGmPdiXNbTgctExANUAmc7n4egn6EYlxXsD7D/GGPKA17a5tfVMRH4FbBSRJY56/6A/aHQ6p9dTeGhlFIqpL2lj0IppVSUNFAopZQKSQOFUkqpkDRQKKWUCkkDhVJKqZA0UKh2T0TKwthnppNuoaXOebKIDGvB433ejNeWOX97isjsEPtlicjl0Z5HqcZooFCdxUwgokAhIu4Qm08GWixQGGMOboFjbDHGnB5ilyxAA4VqcRooVIfh3M/gIycJ4loRedHJwX8V0BNYICILnH2PEZEvRGSpiLzu5Mjx30vgLhFZCpwhIheLyFdOoro3RCRFRA4GpgJ/c+45MEhERjuJC1eIyFsi0sU53kdi77exWES+EZEDReRNJ+//7QFlLwtYvl7s/QyWi8idQd7nAKfsKxsco784908QkeEissgp3woRGQzcCQxy1v1NRNLE3ltlqXOskwKO842IPCH23gb/EZFkZ9s+IjLPKdtSERnkrP+tc51WiMgtLfoPq9q/1p56rg99NPeBc88K7EzaUmyuGhfwBTZBGwTcIwCbDuQTnNQLwPXATQH7/S7g2DkBy7cDVzrLs4DTA7atAA53lm/FSaWAzU91l7N8NbAFe++ARGyWzpwG7+FYbNoH/71P9kilAMwBznOWrwh4bX+clBPAA8C5znICkEyDlBTYmcMZAddkHTa9Q3/sjP7RzrbXgGnO8pfAKc5yEraWdgzwuPNaF/BP4LBYfy700XaPvSKFh+pUFhljNgE4qQz6Awsb7HMQttnoMyfVQgI2qPi9GrA8wvnVnoVN4fJ+wxOKSCaQZYz52Fn1LPB6wC7+NCUrgdXGmK3O69Zjk7EVBex7FPCMce59YoKnUpgInOYsPw/cFWSfL4AbRaQ38KYx5nvnve5WdOCvInIY4MOmlu7mbPvRGLPMWV4C9BeRdKCXMeYtp2xVzvs4Bhssvnb2T8NmfY1lJlrVhjRQqI6mOmDZS/DPsGBv1HJOI8cIzN0zCzjZGLNcRKZjay3RlsnXoHy+RsoXjpC5dYwxL4nIl9hEcHPFpppf32C3c4E8YKwxplZENmBrCYFlBnsdk0OcToA7jDGPRVB+1YloH4XqLHZhbxEJ9k5vE0VkHwARSRWRfRt5XTqwVWwK53ODHc8YUwqUiMihzrZfAR8TnQ+AC/wjtMS5D3MDn2GTKNKgTHXE3r9jvTHmfuAdYBS7XwOw2U+3O0FiMtBvzyPVM/bOaZvEuYGUiCQ65XwfuDCgn6eX2HsxqL2EBgrVWTwOvCciC4wxBdhsui+LyApsM82QRl73J2y7/GfA2oD1rwC/FZGvnQ7d87Gd2yuA0dh+iogZm2l0DrDYaTq7LshuVwNXiMhKGr8T2ZnAKucYI7C3wCzCNretEpG/AS8C45zjnNfg/TXmV9jMxCuwfSndjTH/AV4CvnCONZvdA5Lq5DR7rFJKqZC0RqGUUiokDRRKKaVC0kChlFIqJA0USimlQtJAoZRSKiQNFEoppULSQKGUUiqk/wesXmgyFWLeaQAAAABJRU5ErkJggg==\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": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no bootstrapping\n",
      "Total evaluations for no bootstrapping:\n",
      "6000\n",
      "bootstrapping\n",
      "Total evaluations for bootstrapping:\n",
      "6000\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 != \"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 specifying an integer value to num_bootstrap - which is 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": 11,
   "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(gss=me_gsc, bootstrap=True, num_bootstrap=None, extrapolator=extrap)\n",
    "    # execute\n",
    "    res = bs_extr.sample(es_problem, points)\n",
    "\n",
    "    results_full[f\"{key}\"] = res.raw_results\n",
    "    results[f\"points_{key}\"] = res.points\n",
    "    results[f\"energies_{key}\"] = res.energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "poly\n",
      "diff_model\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Energy')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEWCAYAAACaBstRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd4VGX6//H3PSWNhCRA6IGAdAgECEUjTbALqAi6wkrEhrvrKvvTr7ruurquLq6u69orIoqo6FoXdQVBQVAM0nsxVClJSG9Tnt8fM4khBDJAJiflfl3XXMyc+pmTYe45zznnOWKMQSmllAqEzeoASiml6g8tGkoppQKmRUMppVTAtGgopZQKmBYNpZRSAdOioZRSKmBaNFS9IiIviMifLVjvZyIy9TTn7SAi+SJir+lcNU1ErhCRvf68/UVko4iM9I97QETetDiispjodRqqrhCRdKAV4AY8wCZgDvCSMcZrYbRT4n8fNxpjFlqd5VSJyE7gD8aYj6oY9wDQxRgzpdaDqTpD9zRUXTPWGBMFdARmAncDr1obqWEQEUcAk3UENgY7i6q/tGioOskYk2OM+Ri4GpgqIn0ARGS2iPzN/7yFiHwqItkikiUiS0XE5h93t4jsF5E8EdkqIqP9w0NF5EkROeB/PCkioWXrFZHxIrJGRHJFZKeIXOQfvkREbvQ/P0tEvhKRTBHJEJG5IhLjH/cG0AH4xN/E838ikiAipuxLW0TaisjH/sw7ROSmCut/QETeFZE5/uwbRST5RNvJv9zfi8guf5bHKmyDVBH5VkT+JSKZwAMiYhORP4nIbhE57F9PtH+75AN2YK1/jwMRSReRMSdY91ARWe7f/mvLmrFUw6ZFQ9VpxpiVwD5gWBWj/59/XBy+Zq0/AkZEugO/Awb591ouBNL989wHDAWSgH7AYOBPACIyGF9z2F1ADDC8wnwVCfB3oC3QE4gHHvDn/TWwB98eU6Qx5h9VzP+2P3db4CrgERE5r8L4cf5pYoCPgWeq3jrlrgCSgQHAeGBahXFDgF34ts/DQKr/MQroDEQCzxhjSowxkf55+hljzjrZCkWkHfBf4G9AM+BO4H0Riasmq6rntGio+uAAvi+mylxAG6CjMcZljFlqfAfpPEAo0EtEnMaYdGPMTv88k4G/GmMOG2OOAA8Cv/aPuwGYZYz50hjjNcbsN8ZsqbxSY8wO/zQl/mU8AYwI5I2ISDyQAtxtjCk2xqwBXgGuqzDZMmPMAmOMB3gDX3E7mUeNMVnGmD3Ak8CvKow7YIx52hjjNsYU+d//E8aYXcaYfOBe4JoAm64qmgIs8Of0GmO+BNKAS05xOaqe0aKh6oN2QFYVwx8DdgD/8zfP3AO+L3XgDny//g+LyNsi0tY/T1tgd4Vl7PYPA98ew06qISKt/MvcLyK5wJtAiwDfS1sgyxiTVylDuwqvD1Z4XgiEVfOlvrfSstqeYFzZ+iu/fwe+PZFT0RGY6G+ayhaRbOBcfEVcNWBaNFSdJiKD8H2hLqs8zhiTZ4z5f8aYzviadP5QduzCGPOWMeZcfF9uBnjUP9sB/7AyHfzDwPcFe9JmGb9H/MtMNMY0xferWypGO8m8B4BmIhJVKcP+ANZ7IvGVlnWgwuvKWap6/27g0Cmucy/whjEmpsKjiTFm5ikuR9UzWjRUnSQiTUXkMnxt+28aY9ZXMc1lItJFRATIwdcs5RWR7iJynv8AdzFQBJSdsjsP+JOIxIlIC+B+fHsK4DtL63oRGe0/YNxORHpUES8KyAdy/G37d1Uafwjf8YLjGGP2AsuBv4tImIj0xdcsdibXP9wlIrH+pq/bgXdOMu08YIaIdBKRSHwF8B1jjPsU1/kmMFZELhQRu/+9jBSR9qf3FlR9oUVD1TWfiEgevl+y9+E7XnD9CabtCizE9wW+AnjOGLMY3/GMmUAGvqaelvja7sF34DYNWAesB370Dys76H498C98Rehrjv1VXuZBfAedc/AdDP5PpfF/x1eYskXkzirm/xWQgO9X/wfAX87wmo6PgFXAGn+ek52iPAvfcZJvgJ/wFdXbTnWF/uI3Ht/JB0fw/b3uQr9TGjy9uE+pekxEDNDVfxxHqaDTXwVKKaUCpkVDKaVUwLR5SimlVMB0T0MppVTATvUq0DqvRYsWJiEhweoYSilVr6xatSrDGFNtNzANrmgkJCSQlpZmdQyllKpXRGR39VNp85RSSqlToEVDKaVUwLRoKKWUCliDO6ahlDozLpeLffv2UVxcbHUUFQRhYWG0b98ep9N5WvNr0VBKHWPfvn1ERUWRkJCAry9I1VAYY8jMzGTfvn106tTptJahzVNKqWMUFxfTvHlzLRgNkIjQvHnzM9qL1KKhlDqOFoyG60z/tlo0/HZtXMesm/7EJ889bXUUpZSqs7Ro+NnsNors55Gx7mD1Eyul6qzU1FTee++9M17OmjVrWLBgwUmnSU9P56233jrjdQXihRdeYM6cObWyrpPRouGX0KMPdlc+piTc6ihKqTrgTIuG232qN0M8uenTp3PdddfV6DJPhxaNChzuLIwnxuoYSjVq6enp9OzZk5tuuonevXtzwQUXUFRUBPi+yIcOHUrfvn254oorOHr0aJXLWLhwIcnJyXTr1o1PP/0U8B3gv/7660lMTKR///4sXrz4hMNLS0u5//77eeedd0hKSuKdd97h66+/JikpiaSkJPr3709eXh733HMPS5cuJSkpiX/961/Mnj2bcePGcd555zF69Gjy8/MZPXo0AwYMIDExkY8++qj8Pfbo0YPJkyfTs2dPrrrqKgoLCwFfV0j/93//R2JiIoMHD2bHDt/9tR544AEef/xxAEaOHMndd9/N4MGD6datG0uXLgWgsLCQSZMm0atXL6644gqGDBlS490q6Sm3FYjJxGtva3UMpeqMBz/ZyKYDuTW6zF5tm/KXsb1POs327duZN28eL7/8MpMmTeL9999nypQpXHfddTz99NOMGDGC+++/nwcffJAnn3zyuPnT09NZuXIlO3fuZNSoUezYsYNnn30WEWH9+vVs2bKFCy64gG3btp1w+F//+lfS0tJ45plnABg7dizPPvssKSkp5OfnExYWxsyZM3n88cfLC9Ps2bP58ccfWbduHc2aNcPtdvPBBx/QtGlTMjIyGDp0KOPGjQNg69atvPrqq6SkpDBt2jSee+457rzTd3fg6Oho1q9fz5w5c7jjjjvKl1+R2+1m5cqVLFiwgAcffJCFCxfy3HPPERsby6ZNm9iwYQNJSUln9Leqiu5pVCDOXFzOZhQVFFgdRalGrVOnTuVfeAMHDiQ9PZ2cnByys7MZMWIEAFOnTuWbb76pcv5JkyZhs9no2rUrnTt3ZsuWLSxbtowpU6YA0KNHDzp27Mi2bdtOOLyylJQU/vCHP/DUU0+RnZ2Nw1H1b+7zzz+fZs2aAb7rIv74xz/St29fxowZw/79+zl06BAA8fHxpKSkADBlyhSWLVtWvoxf/epX5f+uWLGiyvVceeWVx2wfgGXLlnHNNdcA0KdPH/r27VvlvGdC9zQqsDUpxZQ62bD8awadf4nVcZSyXHV7BMESGhpa/txut5c3TwWq8mmlNXEK8T333MOll17KggULSElJ4YsvvqhyuiZNmpQ/nzt3LkeOHGHVqlU4nU4SEhLKr5E4WcYTPa+obBvZ7fYaP35yMrqnUUFEK98fe8/aDRYnUUpVFh0dTWxsbHn7/RtvvFG+11HZ/Pnz8Xq97Ny5k127dtG9e3eGDRvG3LlzAdi2bRt79uw56fCoqCjy8vLKl7lz504SExO5++67GTRoEFu2bDlumspycnJo2bIlTqeTxYsXs3v3L72P79mzp3wv4q233uLcc88tH/fOO++U/3v22WcHvI1SUlJ49913Adi0aRPr168PeN5A6Z5GBW16duPIXsjfn211FKVUFV5//XWmT59OYWEhnTt35rXXXqtyug4dOjB48GByc3N54YUXCAsL4ze/+Q233noriYmJOBwOZs+eTWho6AmHjxo1ipkzZ5KUlMS9997LsmXLWLx4MTabjd69e3PxxRdjs9mw2+3069eP1NRUYmNjj8kxefJkxo4dS2JiIsnJyfTo0aN8XPfu3Xn22WeZNm0avXr14tZbby0fd/ToUfr27UtoaCjz5s0LePv85je/YerUqfTq1YsePXrQu3dvoqOjT3Ern1yDu0d4cnKyOd2zBXIyM3jzj2sI937NtJcfquFkStUPmzdvpmfPnlbHaNDS09O57LLL2LDh+FaNshvJtWjR4pSX6/F4cLlchIWFsXPnTsaMGcPWrVsJCQk5Zrqq/sYissoYk1zdOnRPo4Lo5i1wlmaDN8rqKEopdcoKCwsZNWoULpcLYwzPPffccQXjTGnRqMTuzcQQW/2ESil1mhISEqrcywDKz4Q6HVFRUUG/3bUeCD9OFh57c6tDKKVUnaRFoxIJKcAVEsORA3utjqKUUnWOFo1KnNG+EwPWf73Y4iRKKVX3aNGoJKq9r2nq8LZ0a4MopVQdpEWjks6D+gNQfETvj6yUVdLT0+nTp88ZL+fJJ58s7wjwRGbPns2BAwfOeF2BOOecc2plPcGkRaOSXoPPxeYpxVMYWv3ESqk67UyLhsfjqdE8y5cvr9HlWUGLRiUOpxOnKxM8Ta2OolSj5na7j+s6fNGiRfTv35/ExESmTZtGSUkJQJXDn3rqKQ4cOMCoUaMYNWoUHo+H1NRU+vTpQ2JiIv/617947733SEtLY/LkySQlJVFUVERCQgJ33303AwYMYP78+bz88ssMGjSIfv36MWHChPIilJqayvTp04/rgn327NmMHz+ekSNH0rVrVx588MHy9xQZGQnAkiVLGDlyJFdddVV5F+llF1ovWLCAHj16MHDgQH7/+99z2WWX1eZmr5Yl12mISDPgHSABSAcmGWOq7BhfRJoCm4APjTG/q5V83kyM6Gm3SvHZPXCwhvsvap0IF8+sdrLKXYc/8cQTvPjiiyxatIhu3bpx3XXX8fzzzzN9+nRSU1OPG37HHXfwxBNPsHjxYlq0aMGqVavYv39/+fUR2dnZxMTE8Mwzz/D444+TnPzLxdDNmzfnxx9/BCAzM5ObbroJgD/96U+8+uqr3HbbbUDVXbADrFy5kg0bNhAREcGgQYO49NJLj1k+wOrVq9m4cSNt27YlJSWFb7/9luTkZG655Ra++eYbOnXqVN7bbV1i1Z7GPcAiY0xXYJH/9Yk8BFTd/3Gw2HNwO5vjdrlqdbVKqV9U7jp80aJFdOrUiW7dugG/dI2+devWKodX1rlzZ3bt2sVtt93G559/TtOmJ25NuPrqq8ufb9iwgWHDhpGYmMjcuXPZuHFj+biqumAHX/fozZs3Jzw8nCuvvPKYbs/LDB48mPbt22Oz2UhKSiI9PZ0tW7bQuXNnOnXqBFAni4ZVV4SPB0b6n78OLAHurjyRiAwEWgGfA9X2iVJTbOFFeLzh7Fq/mm4DBtfWapWqewLYIwiWyl2Cx8TEkJmZedrLi42NZe3atXzxxRe88MILvPvuu8yaNavKaSt2b56amsqHH35Iv379mD17NkuWLDlhxrLXgXTNXrn799rs3vxMWLWn0coY87P/+UF8heEYImID/gncWd3CRORmEUkTkbQjR46ccbjQZk4Atq1cecbLUkqdnspdhycnJ5Oenl7eBFTWNXr37t2rHA4c03V5RkYGXq+XCRMm8Le//a28+am67s3z8vJo06YNLpervAv1MlV1wQ7w5ZdfkpWVRVFRER9++GH5HlN1unfvzq5du8q7EinrIr0uCdqehogsBFpXMeq+ii+MMUZEqupq9zfAAmPMvupuoGKMeQl4CXy93J5e4l80P6s9RzMg+6dDZ7oopdRpqtx1+FNPPcXQoUOZOHEibrebQYMGMX36dEJDQ3nttdeOGw5w8803c9FFF9G2bVuefPJJrr/+erxeLwB///vfgV8OaIeHh1d5l7yHHnqIIUOGEBcXx5AhQ44pMFV1wQ6+pqcJEyawb98+pkyZctzxjBMJDw/nueee46KLLqJJkyYMGjTojLZhMFjSNbqIbAVGGmN+FpE2wBJjTPdK08wFhgFeIBIIAZ4zxpzs+McZdY1eZu/2rXz8z/2Es4hpLzx8RstSqr7RrtEDk5qaymWXXcZVV111zPDZs2cfc2/xU5Wfn09kZCTGGH7729/StWtXZsyYURORy51J1+hWNU99DEz1P58KfFR5AmPMZGNMB2NMAr4mqjnVFYyaEt+1Ow5XHqY4ojZWp5RS5V5++WWSkpLo3bs3OTk53HLLLVZHOoZVB8JnAu+KyA3AbmASgIgkA9ONMTdalKuc3Z0JJsbqGEqpOmr27NlVDk9NTSU1NfW0lztjxowa37OoSZYUDWNMJjC6iuFpwHEFwxgzG5gd9GAVCFl47fG1uUqllKrz9IrwE3Hk4XI2o6igwOokSilVZ2jROAFHpAtjs7P2m0VWR1FKqTpDi8YJRLTy9RGzb/0mi5MopVTdoUXjBNr29p0BXPBzrsVJlFKq7tCicQJ9R4wG48WTZ9UJZkopVfdo0TiBqJgYQkqzMK5Iq6Mo1eikp6fTs2dPbrrpJnr37s0FF1xAUVERI0eOpOzi3YyMDBISEgDf6a+XX345559/PgkJCTzzzDM88cQT9O/fn6FDh5KVlQXAyJEjuf3220lKSqJPnz6sXLkSr9dL165dKeuCyOv10qVLF2qiS6KGSH9Gn4TNmwU0szqGUpZ5dOWjbMnaUqPL7NGsB3cPPq5/0uNs376defPm8fLLLzNp0iTef//9k06/YcMGVq9eTXFxMV26dOHRRx9l9erVzJgxgzlz5nDHHXcAUFhYyJo1a/jmm2+YNm0aGzZsYMqUKcydO5c77riDhQsX0q9fP+Li4mrk/TY0uqdxMnIUt13vq6GUFTp16kRSUhIAAwcOLO/E70RGjRpFVFQUcXFxREdHM3bsWAASExOPmbesu/Hhw4eTm5tLdnY206ZNY86cOQDMmjWL66+/vubfUAOhexonISH5uG3RHN67m5bxHa2Oo1StC2SPIFgqdx1eVFSEw+Eo73CwuLj4hNPbbLby1zab7Zhux6vqtjw+Pp5WrVrx1VdfsXLlyuN6s1W/0D2Nk3DG+D5c679ZbHESpRRAQkICq1atAuC99947rWWUdTe+bNkyoqOjiY6OBuDGG29kypQpTJw4EbvdXjOBGyAtGifRNN7XNHV4226LkyilAO68806ef/55+vfvT0ZGxmktIywsjP79+zN9+nReffXV8uHjxo0jPz9fm6aqYUnX6MFUE12jl9mwYhlfv15KhGMh1z/zSI0sU6m6riF3jT5y5Mjj7gdeJi0tjRkzZrB06VILktWu+tg1er3QI3kINk8JpjC0+omVUvXWzJkzmTBhQvmNmdSJ6YHwk3A4nThdmRhvtNVRlFI1oOL9vSu65557uOeeWrldT72nexrVEG8mxqan3SqlFGjRqJbYc3A5m+N2uayOopRSltOiUQ0JL8ZrD2PHuh+tjqKUUpbTolGNsBYhAGz/7juLkyillPW0aFSjRRffLV9z9mjnZUrVVRU7MlTBpUWjGr2HDQfAld2wrmdRSqnToUWjGm07dcXhysWURFgdRalGIz09nR49ejB58mR69uzJVVddRWFhIYsWLaJ///4kJiYybdo0SkpKjplv1qxZ5b3ZArz88svMmDGjtuM3aHqdRgDs7kyM0S7SVeNz8JFHKNlcs12jh/bsQes//rHa6bZu3cqrr75KSkoK06ZN44knnuDFF19k0aJFdOvWjeuuu47nn3/+mCIxadIkHn74YR577DGcTievvfYaL774Yo3mb+x0TyMAQhZeuxYNpWpTfHw8KSkpAEyZMoVFixbRqVMnunXrBsDUqVP55ptvjpknMjKS8847j08//ZQtW7bgcrlITEys9ewNme5pBECceZTam1GQm0OTpnp1uGo8AtkjCJbKXZjHxMSQmZlZ7Xw33ngjjzzyCD169NDOB4NA9zQCYI90g9hZt1S7SFeqtuzZs4cVK1YA8NZbb5GcnEx6ejo7duwA4I033mDEiBHHzTdkyBD27t3LW2+9VX7DJVVztGgEIKJ1FAD71m+2OIlSjUf37t159tln6dmzJ0ePHmXGjBm89tprTJw4kcTERGw2G9OnT69y3kmTJpGSkkJsbGwtp274tHkqAO0Te3I4HQoP5lkdRalGw+Fw8Oabbx4zbPTo0axevfq4aSt3RLhs2TI9aypIdE8jAH2HjQLjwZOvNVapuiw7O5tu3boRHh7O6NGjrY7TIOm3YACaNI0mpDQL442yOopSjUJCQgIbNmw45fliYmLYtm1bEBKpMrqnESCbJwuDnnarlGrctGgESGxZeBx6Xw2lVOOmRSNAElqI29mUg7t3Wh1FKaUso0UjQM4Y34VG679eYm0QpZSykBaNAEV3iAMgY+c+i5MopZR1tGgEqMugZACKj5RUM6VSqqY98MADPP7449x///0sXLgQgKVLl9K7d2+SkpIoKirirrvuonfv3tx1111BzzN79mx+97vfnfE09ZGechugrv0Hs8TzJaYkzOooSjVaf/3rX8ufz507l3vvvZcpU6YA8NJLL5GVlYXdbrcqXqOgRSNADqcTpysT49UOC1XjsfTdbWTsza/RZbaIj2TYpG7VTvfwww/z+uuv07JlS+Lj4xk4cCCpqalcdtllZGdn8+677/LFF1/w2WefkZeXR35+PgMHDuTee+/l6quvPm55qamphIeHs3r1ag4fPsysWbOYM2cOK1asYMiQIcyePRuAefPm8cgjj2CM4dJLL+XRRx8F4LXXXuPvf/87MTEx9OvXj9DQUACOHDnC9OnT2bNnDwBPPvlkee+8DZEWjVMg3kyMTU+7VSrYVq1axdtvv82aNWtwu90MGDCAgQMHlo+/8cYbWbZsGZdddhlXXXUV4OsWfc2aNSdd7tGjR1mxYgUff/wx48aN49tvv+WVV15h0KBBrFmzhpYtW3L33XezatUqYmNjueCCC/jwww8ZMmQIf/nLX1i1ahXR0dGMGjWK/v37A3D77bczY8YMzj33XPbs2cOFF17I5s0Nt586LRqnQBw5lNq743a5cDidVsdRKugC2SMIhqVLl3LFFVcQEeG7Y+a4ceNqZLljx45FREhMTKRVq1bl99ro3bs36enp7N69m5EjRxIX5zvxZfLkyeX37Kg4/Oqrry6/8nzhwoVs2rSpfB25ubnk59fs3lldYknREJFmwDtAApAOTDLGHK1iug7AK0A8YIBLjDHptRa0cp6IErzuULakfU+fs8+1KoZS6jSVNSnZbLby52Wv3W43ztP4Mej1evnuu+8IC2scxzutOnvqHmCRMaYrsMj/uipzgMeMMT2BwcDhWspXpbDmvg/ZzrQfrIyhVIM3fPhwPvzwQ4qKisjLy+OTTz6plfUOHjyYr7/+moyMDDweD/PmzWPEiBEMGTKEr7/+mszMTFwuF/Pnzy+f54ILLuDpp58uf11dE1l9Z1Xz1HhgpP/568AS4O6KE4hIL8BhjPkSwBhj+f5ey24dyToEuXurv3uYUur0DRgwgKuvvpp+/frRsmVLBg0aVCvrbdOmDTNnzmTUqFHlB8LHjx8P+E77Pfvss4mJiSEpKal8nqeeeorf/va39O3bF7fbzfDhw3nhhRdqJa8VxBhT+ysVyTbGxPifC3C07HWFaS4HbgRKgU7AQuAeY4yniuXdDNwM0KFDh4G7d+8OSu7De3cz/+GdhHsXMe2lh4OyDqWstnnzZnr27Gl1DBVEVf2NRWSVMSa5unmD1jwlIgtFZEMVj/EVpzO+qlVV5XIAw4A7gUFAZyC1qnUZY14yxiQbY5LLDlQFQ8v4jjhKczClkUFbh1JK1WVBa54yxow50TgROSQibYwxP4tIG6o+VrEPWGOM2eWf50NgKPBqUAIHyOHJBPQWkkrVVQ8//PAxxxwAJk6cyH333WdRoobFqmMaHwNTgZn+fz+qYpofgBgRiTPGHAHOA9JqL+KJZOG1JVgdQqmgMsbgazmuf+677z4tECdxpockrDp7aiZwvohsB8b4XyMiySLyCoD/2MWdwCIRWQ8I8LJFecuJM5/SkFjysrOtjqJUUISFhZGZmXnGXy6q7jHGkJmZeUanB1uyp2GMyQSOu4GvMSYN38HvstdfAn1rMVq17FFuKLKz7utFpIyfYHUcpWpc+/bt2bdvH0eOHLE6igqCsLAw2rdvf9rz6xXhp6hJm6bk74IDm7b6ThxWqoFxOp106tTJ6hiqjtKu0U9Rh6Q+ABQetPyyEaWUqnVaNE5RYsooxOvGk687aUqpxke/+U5ReJMmOF1ZGG9Tq6MopVSt0z2N02DzZGJoZnUMpZSqdVo0ToPYj+J2tsTtclkdRSmlapUWjdNgjy7E44hg+SfvWR1FKaVqlRaN09Au2XdjmvTlGyxOopRStUuLxmkYftU12N1FuDPCrY6ilFK1KqCiISL/EZFLRUSLDBASFoazdDdec/pXVSqlVH0UaBF4DrgW2C4iM0WkexAz1QsScpCSsHYc3L3T6ihKKVVrAioaxpiFxpjJwAB89/ReKCLLReR6ETn1m+o2ABHt7SB2vp3/vtVRlFKq1gTc3CQizfHdBOlGYDXwb3xF5MugJKvjksZeBEDujhyLkyilVO0J9JjGB8BSIAIYa4wZZ4x5xxhzG9Aob2PXY+AQQkoO4S1qYXUUpZSqNYF2I/KUMWZxVSMCuadsQ2Xz7sFt74bb5cLhbJStdEqpRibQohErIldWGpYDrDfGVHWr1kbB3jSHYlc0aQs/Y+jF46yOo5RSQRfoMY0bgFeAyf7Hy8DdwLci8usgZavz4hJ9p9zuWLLS4iRKKVU7Ai0aTqCnMWaCMWYC0AswwBB8xaNRGj7pWmyeUkqPaNOUUqpxCLRotDfGHKrw+jAQb4zJAhptr31RMTGElO7B69GL/JRSjUOgxzSWiMinwHz/6wn+YU2A7KAkqyfEvp8SxzkcPXKQ2LjWVsdRSqmgCnRP47fAa0CS/zEH+K0xpsAYMypY4eqD0DZgbE6+eXue1VGUUiroqt3TEBE7sNBfHPTy50q/Lwd/AAAgAElEQVR6nz+Mb9+GrM1HrI6ilFJBV+2ehjHGA3hFJLoW8tQ7SSPH4CzJwlugd/JTSjV8gR7TyAfWi8iXQEHZQGPM74OSqp6xe3fjsXW0OoZSSgVdoEXjP/6HqoKtSSbF3v6sXfoV/YadZ3UcpZQKmoCKhjHmdREJBzoYY7YGOVO9E9s9jsLNsPGLr7VoKKUatEA7LBwLrAE+979OEpGPgxmsPhl29STE66b4Z6uTKKVUcAV6yu0DwGD812QYY9YAnYOUqd5p3rodISV7Ma62VkdRSqmgCrRouIwxlW8c4a3pMPWZzbaf0tCOFOTq/TWUUg1XoEVjo4hcC9hFpKuIPA0sD2KuesfZshSvPYSv33nL6ihKKRU0gRaN24DeQAkwD8gF7ghWqPqoy3DfbUUOr9trcRKllAqeQO8RXmiMuc8YM8gYk+x/XhzscPXJoAsuxVGaiydXr4FUSjVcAZ1yKyLdgDuBhIrzGGP0/FI/h9OJw52O197B6ihKKRU0gV7cNx94Ad+NmDzBi1O/ScQRSunLth9X0m3AYKvjKKVUjQu0aLiNMc8HNUkDEH1WNEU74cdPPteioZRqkAI9EP6JiPxGRNqISLOyR1CT1UMpkyaA8VK4t9Hel0op1cAFuqcx1f/vXRWGGfQCv2O07ngWocXLMEZvxqSUapgC7XuqU7CDNBQ2215K7f0pLS4mJCzM6jhKKVWjTto8JSL/V+H5xErjHglWqPrMEVuIxxHO0vfftjqKUkrVuOqOaVxT4fm9lcZddLor9R8T+VJEtvv/jT3BdP8QkY0isllEnhIROd111pYOQ3sDsO+HbRYnUUqpmldd0ZATPK/q9am4B1hkjOkKLPK/PnbhIucAKUBfoA8wCBhxBuusFeeMm4DdXYAnO9LqKEopVeOqKxrmBM+ren0qxgOv+5+/Dlx+gnWHASFAKOAEDp3BOmtFSFgYztLdeIm3OopSStW46opGPxHJFZE8oK//ednrxDNYbytjTNndJw4CrSpPYIxZASwGfvY/vjDGbK5qYSJys4ikiUjakSNHziBWzZDQg5SEtmHPtk1WR1FKqRp10qJhjLEbY5oaY6KMMQ7/87LXzpPNKyILRWRDFY/xldZhqGKvRUS6AD2B9kA74DwRGXaCnC/5+8RKjouLq+YtB19kQhiIje/f/9DqKEopVaMCvU7jlBljxpxonIgcEpE2xpifRaQNcLiKya4AvjPG5Pvn+Qw4G1galMA1KPnycXz2dAZ5PxVYHUUppWpUoFeE17SP+eWCwanAR1VMswcYISIOEXHiOwheZfNUXdO5d19Cin/GFB/X6qaUUvWaVUVjJnC+iGwHxvhfIyLJIvKKf5r3gJ3AemAtsNYY84kVYU+HzezB7UzA7dIuRZRSDUfQmqdOxhiTCYyuYngacKP/uQe4pZaj1Rh7TB7FJZGs+OQDhl05yeo4SilVI6za02jw2iT5el75aflai5MopVTN0aIRJMMnXYvNU4wrU/ufUko1HJY0TzUG4U2aEFKyGyN6kZ9SquHQPY0gkpCfKQ1tx+G9u62OopRSNUKLRhCFt7NjbHaWzX/X6ihKKVUjtGgEUf+xF4Lxkr0p1+ooSilVI7RoBFGP5KGEFW3H4+mj12sopRoELRpB5mi+l9LQOBa88KzVUZRS6oxp0QiyEbdci3hdZPyoTVRKqfpPi0aQJfToQ2jJBly2vhTk5lgdRymlzogWjVoQFp+L29mUT598yuooSil1RrRo1IJLbr8Vu7uQgp2hVkdRSqkzokWjFsTGtcbpXktJSCIHd++0Oo5SSp02LRq1JKa3A689lC+ffc3qKEopddq0aNSSsbfdjrM0i9KDLa2OopRSp02LRi0JCQvDLuspCevJllXfWx1HKaVOixaNWtT+3PYYm53vZld1d1ullKr7tGjUogtSbyS0eD+e3E5WR1FKqdOiRaMWiQi2sC0Uh5/Fik8+sDqOUkqdMi0atazX+CEAbP4kzeIkSil16rRo1LKhF48jrGgHntJe2vOtUqre0aJhAXv0T5SGteF/s1+xOopSSp0SLRoWOCf1CsTr5uCKQ1ZHUUrVMo/X8MiCzdw8J40vNh7E5fFaHemUOKwO0Bh1GzCYpcX/xu3oS1FBAeFNmlgdSSlVC4pdHma8s4bPNhwkJsLJ/zYdokVkKBMGtGNicjxdWkZaHbFauqdhkdA2R3CFxPDp09rzrVKNQV6xi9TXVvLZhoP8+bJepN03hleuS6Z/hxheWfYTY574mqueX867aXspKHFbHfeExBhjdYYalZycbNLS6v6ZSYf37ub9v24kpHQtN7x+r9VxlFJBdCSvhNTXVrL1YB6PT+zH5f3bHTP+cF4x//lxP+/+sJddGQU0CbEztl9bJg2Kp398DCIS9IwissoYk1zddNo8ZZGW8R0JKX0bl7MvmQf307x1u+pnUkrVO3syC/n1rO85nFvCK1OTGdn9+P7nWkaFMX3EWdwyvDNpu4/yzg97+WjNAd7+YS+dWjSheZOQ8mkr/8yv+MO/W6soZk7oG6y3AmjzlKWiurnwOML5/KkXrI6ilDqBDftzWLz18GkdsN54IIcrn19OTpGLt24aUmXBqEhEGJTQjMcn9uOHP41h5pWJdGrRhFCnrfwR5rQR7rSXPyJCHDQJ9T3CnPbTfZsB0z0NC1162228ecdiivc3szqKUqqSwlI3//h8K6+vSMcYaBkVyqTkeK4eFE98s4hq51+xM5Ob56QRFebg7ZvPpkvLqFNaf2Sog2sGd+CawR0Cmr60qJCCjIOntI7ToUXDQk2aRuM06ygOHcKO9avpktjf6khKKWD5jgzu/s869mYVMfXsjpzTpQXv/rCX55bs4NklOxjWNY5rB8czumcrnPbjG2w+3/Azv397DR2aRTBn2mDaxoSf0vpLiwo5tG0dh3dsJ3v/IQqPFlKc58JVJLhLnHg9oXi94Rgi8Nqa4LE3wWMPJ6z4J254vXNNbYYqadGwWNygWPasc/LtK/Pp8m8tGkpZKa/YxSMLtjBv5R4Smkfw7i1nM7iTryXgwt6tOZBdxLtpe3n3h71Mf/NHWkSGMjG5PdcMiqdjc9+p8/NW7uG+D9bTLz6GWVMHEes/HuEuKeHg1rUc2LSRo3sPkZ9RSEkBuIudeNzhGNMEr0ThsUfidkSC2IB2/scvbFKMnQJsFGCjEKc5Sqgpxi4uIuJr4YC5nj1lLbfLxWs3v4fNk8MNc6ZbHUepRmvJ1sPc+5/1HMot5oZzO/GH87sTHlL1MQKP1/DNtiO8tXIPX205TLgrl/HRh2lbnEXB4TxiTAihJgKvOxKvicJji8HtbIqxHf873eYpxuHOw+bNw0Y+NnsRjpASnBEewmOcRDaPJLp1S5p16EjLLj2IiI0LyvvXs6fqCYfTid2xiaKQEaxa9BkDR19sdSSlGpWcQhcP/XcT763aR5eWkbx36zkM6BB7zDQZP20hPW0lGbv2kXe4mOI8O+7SSIaaGAbZY3E5W0Ce7yB3Exu4AK+7CAc52MghhF2ESxGh4W4imoUQ3bYZcZ0SaNMrieh6duakFo06oOvFiaz7n4117y3VoqFULfrfxoPc9+EGsvKKmNGtmH6unex+4Ss2Zwmu4kg8phkuRws8zkigvf8BggunHMVujhJittHEXkRYtMHRPIKEnl2I79ufZvFnWfregkWLRh0w7MpJbP1oFh7pRWlxMSFhYVZHUqpBSt+ygbVffkn2riO4chzgjuEWaYErJA7vyuZspLlvQuMlRLKwmwzCWU+Is4QmLRw06xBH+z59aN83BUdoqLVvxiJaNOqIkLid5OWP5u0/PcR1jz9sdRyl6i23y8WGb79m+4qV5O/Pw1sQAZ4WeBytcIXEAP0AEIcbpzcDu8mkidlDWHgJ0e2a0K5PN84aOpLwGD0VvipaNOqIa/72Z+bc8hbFpQP5OX0HbRK6WB1JqTpvz7ZN/Pjfz8jelYk3PwLjbY3L2QaPIxwYCoDdXoTDewi7ZxsO8ghr6aRd3+4kX3AJTZpGW/sG6iEtGnVESFgYLfpns39LAp89/DrTXn7I6khK1RmlxcV8t+AD9qZtouQImJLmeO1tKQ1tAfhOVbfbC3F69uN0ryY0rJDIdlF0GZpM4rkX4HA6rX0DDYieclvHzJr6OCUhfRg40c3gCy6zOo5Sta60uJjlH7/P3pWbcGWFYDxtKA2Jx2v3H0MwXkJKDmMzPyOhWYS3CeGslGQGjNLicCbq9Cm3IjIReADoCQw2xlT5LS8iFwH/BuzAK8aYmbUW0iJ9Jvdg1XxYP3ezFg3V4LldLlYu+Jhdy9dQesSG8bTB5YzH42gDtMFmLyXEs5dQ9w/YI4tp3r0tQ8aNJ67tGKujN1qW7GmISE/AC7wI3FlV0RARO7ANOB/YB/wA/MoYs+lky67vexoAs278M0WOUbTrsYbL7/iD1XGUqjF7tm3iu/kfkJ9ehClpjSukIx6H70pq8boIKdmPTQ5gjy2kTVJnUq6YqMcdakmd3tMwxmwGqusjfjCwwxizyz/t28B44KRFoyG46L7r+PSh9WSsjtFTcFW95Xa5WP7Jf0j/dh2uzAgM8ZSEtgU5G2xeQm0/E+Jeh71JAXF94km5ciLRzS+0OraqRl0+EN4O2Fvh9T5gSFUTisjNwM0AHToE1iNkXda2U1fCYmeTV6Cn4Kr6I/doNl/NeY2szRmYwpa4nQm4nXHAaOyOIpylPxHu3U5kQhhDJ15Bh27axFQfBa1oiMhCoHUVo+4zxnxUk+syxrwEvAS+5qmaXLZVrnn4l1NwD+7eSeuODfPqUlV/HT1ykCVvvMnRLdmY4raUhnbCa/ddAxFiO4TTs5mwyFzaJndhxMRrdY+5gQha0TDGnOnPiP1AfIXX7f3DGoWKp+Au+NtsPQVXWS7z4H4Wz3mD3O0FmNJ2lIYk4LUPABuE2vYT4k4jLM5F0rgL6T30V1bHVUFSl5unfgC6ikgnfMXiGuBaayPVrsvv+IP/FNwUVv7vUz2bqgHyeg1eY3BUcU8GqxUVFLBw9iscWXsYUxxPSVgnjG2w73gE+wn1fE9YKy/JV15KtwG/tjquqiVWnXJ7BfA0EAf8V0TWGGMuFJG2+E6tvcQY4xaR3wFf4DvldpYxZqMVea1UdgruBj0Ft8Hxeg03v5HG+v05PHPtAAYlWNtthdvlYvnH77FryXo8uXG4nF3xOBJ9exKyl1D3d0TEw6CrxtElUY9HNFZ6cV89oKfgNkyvLvuJhz7dREyEk/xiN/de0pNpKQnVnVVYozavXEba+59TerAJHltXXKG+wuUsycJuthPaqoCBEy6m5+CUWsukrFGnT7lVp0ZPwW14Nh7I4dHPtjCmZ0v+OSmJO+ev5aFPN/HjnqM8OqEvkaHB+a9ZkJvDl7NeIWtDNl7XWZSEdwDOw+4sxFm6nRDnj5w1Komzx16hV1erKumeRj0x5//dR17BaKIiF+kpuPVcUamHsc8sI7fIxed3DKdZkxC8XsOL3+zisS+20DkukhemDKBLy6gaWd+6ZYtZ8/6XuLKa43J2811MZ7yEFf+EhO4mrl8cY1JvJLxJkxpZn6qfAt3T0KJRT5QWFzPnlrfw2mMY95d+egpuPXbfB+uZ+/0e3rxhCOd2bXHMuOU7M/j9vNUUlXp49Kq+XNa37SkvvyA3hy9eeYnsjXl43WdREu47CdFZmo3ds5XQ1oUMnXIlXRL1nvTqF9o81cDoKbgNwxcbDzL3+z3cMrzzcQUD4JyzWvDpbcP4zdxV/O6t1fy4O5t7L+mBs5qzq3asTWP5mx/gPhxLqbMHHsdAcHgIc+0iwiyi9eAOnH/dNBzOK4P11lQjoXsa9UxZL7hJV5Ry9iXjrI6jTsHBnGIu+vc3tI8N5z+3phDiOHEhKHV7eWTBZmYvTye5YyzPTh5Aq6a/HMtyu1x88948dn+9BW9hR4rDO4PYcbhycbg3E9auiJTrJpHQo09tvDXVAGjzVAO18n+fsmq+Dacri3bXtuXiC86zOpIKgNdrmPLq96zek81/f38uneMiA5rv47UHuOf9dUSEOHjsko5kfD6fvO1uPHSnNDQOgNCivdicO2mR1JyLbrxFT5RQp0WbpxqowRdcxqEN/2Dvln7sn5vOP7MWcvvE8+rkxWHqFy8t3cXynZk8OiEx4IIBMCQ8kwfcy8lLb8rW9BI8jqHYHKWElGyjiXMNvccNY9D5U4OYXKljadGoh8b+4f/49O8Pszt9MDELdjDtcDGPpY45pvlC1R3r9mXz+BdbuSSxNZOS46udfsvi/7L24+/Iz2pHcdhZICNxOHIIca8lvF0p592YSqv4i2ohuVLH0+apeuzDPz/A/sPnEla0kVc7JfCPKSlVHlxV1ikocXPpU0spdXv57PbhREccf+2Dp9TFyndeY+e3ByguPouSsHYAhBb/TFjodhKGtmTotTfgCA2t7fiqEdHmqUbg8oceYP7/+zOHZRQ3pKcx9RUvvxvTm9+P7ordVntXFasTe+DjjezOKuTtm4YeUzCKsjNZ9uqrHNjkpsT0wBXSBUxnwsxOmocups/FA+hz0WQLkytVNS0a9dzEfz7EW7+7n6PhI/njoeU89KWdtN1ZPHl1f+Ki9JeplT5dd4D5q/Zx23ldGNK5OUd2bmb56++StSeSEkcPPI5kbLYSwkq3ENdsLUOmjKNtr1usjq3USWnzVAPgcbt56zcPk2sbRhPPEv7e6myahjt56pr+nH1Wc6vjNUr7jhZy8b+XMsrxEylH91OQ0Ybi0K4Ym++02DCzmVY9DOekptK05alfwKdUTdNTbhsZj9vNmzf/g/yQoUTbv2JW5zGkZxRw54XduXXEWbXaCV5j5na5WPTmbHYv+wmbuysl4R0BCCk+RFjINhIGN2fI5GmEhEdYnFSpY2nRaIRKiwqZe+tzFIYNIC5qMUsHXMknaw9w3dkdeWBsb2x6nCMo8rKz+fz558jb5sZDD0pDfScjhBX9RETUbnpe2JuksVdbnFKpk9MD4Y1QSHgEv3r6Rub97g2OmBGM2fIxbYZP4qVvdpFd6OLxif1OehWyCtyujev4ds67uA5GU+rsWeH6ia1Ehq6h/6Tz6Ztyg9UxlapxuqfRAOUePsD8uz6hJLQznTr/wO6zr2bmZ1sY0S2O56cMICJEfyucjuWffMDWL1bize1ASVgX//GJPBzuTYR3KGH0TVNp1aGT1TGVOi3aPNXIZaRv46O/fEtxWDyRrpU4xg3iz9+XkhQfw6zUQcREhFgdsc7LPXqU/730ArlbS/B4u1Ia1gaAkOKfsdu30axvNBfeeIt2Ka4aBC0aiqw9O/j0obfIs52Nw11IRMxK/hYylI4tmzJn2hBaR+sV5JWt+3Yxq9/7EndmC1wh3fE4whGvi9DiHdgi99NlTF+GXXmN1TGVqnFaNFS5H96dzbr/uikO70xY0U5WNfuZda3P5Y0bhtCpReP+lVxUUMCXs14mY10m3tKyO9n57j3h8G4mrH0pI2+YQttOXS1OqlRwadFQx3CXlPDxXx7myJFk3I5wnCXLeD++G/+cfhl92kVbHa/WuF0uvv1oPj99vQFPbhwuZ1c8jogKd7JLp82QsntP6O1OVeOhRUNV6efNq/nfY5+THzIEZ+lRskJXMHzGXZzTrZXV0YJm3beLWfPhQlyHI/HYu+IKaQaAsyQTu9lBSKsChlw9jm4DBlucVCnraNFQJ7X01WfZujSKkrD2hBVuotn5MYyfPJmMghIO5ZRwMLeYgzlF/n9LOJRbzMHcYg7lFBMd4eSW4Z2ZmBxPmNNu9Vs5zrrlX7NhwVcUH7Dh9XYu7wDQ7srH6d6OPTqDLqMHMPSSy3VvQik/LRqqWiX5efznT/8gJ38oXpuD0OKNHA3dzxdNz2K703cls90mxEWG0io6jNZNQ2ndNIwNB3JZtfsorZqGMn3EWfxqcAfLikdedjbfvD2Xwxv24c1rhtfWsfziOpunlJCSXUjEflr1b8t5U1L1TCelTkCLhgrY5m+/YsUry3FLL1whMf72/XTCInfT/bzu9B9/DXbHL9d2GGNYsTOTfy/azvc/ZdEiMpTpIzpz7ZAOQb0GxO1ysen7ZWz4YgnFB8C42lIa2hGv3Xf6sLM0G7s7HWmSQbPuLRh2zdU0b90uaHmUaki0aKhT5i4pIW3+HHZ+u4+iwk7lZxKFlGQQat9Gm96hnJt6PeGxv9yz47tdmTz91Xa+3ZFJ8yYh3DS8M78e2pEmoadfPApyc/hx4WfsW7uF4kOleIuigBa4HK3wOH13vROvm5CSfdhs+whpWUq30UMZMOpCbW5S6jRp0VBnbNd3i/nx/cXkHm5BSUg3vPYQbJ5iQkvTsdnysTuKcYS5CIuy4Q53sq1ESMsP53BUR649rx/Xnd2RqDDfl3hRQQEH0rdzZPducg4doiAzh5K8QtwFLjzFXrzFDiiNxWtrSWlIHMb2S9FxlObgcB9C7BlIeCGxXZtz7jVX06JNe6s2jVINjhYNVaNyft7HijlzOLRVcHta4bFF4XFElTcNVWZ3F2L3FGDEgccejtdezYWExkNIyRFs3sOI8yj2pi5iOseReN5oOvfuG4R3pJSqSDssVDUquk17Lrr7j8cM87jd5P68h5+3biBz917yDmdTdLSEkgIodTnxesMwuHFQCrgRpwdbKDjC7YREhhMWE0XTls2JbduOTr0GEBUTY82bU0oFTIuGOm12h4PY+M7Exne2OopSqpZoP9lKKaUCpkVDKaVUwLRoKKWUCpgWDaWUUgHToqGUUipgWjSUUkoFTIuGUkqpgGnRUEopFbAG142IiBwBdp/BIloAGTUUJ9jqU1aoX3nrU1aoX3nrU1aoX3nPJGtHY0xcdRM1uKJxpkQkLZD+V+qC+pQV6lfe+pQV6lfe+pQV6lfe2siqzVNKKaUCpkVDKaVUwLRoHO8lqwOcgvqUFepX3vqUFepX3vqUFepX3qBn1WMaSimlAqZ7GkoppQKmRUMppVTAGk3REJGLRGSriOwQkXuqGP8HEdkkIutEZJGIdKwwziMia/yPj+tI3lQROVIh140Vxk0Vke3+x9Q6kPVfFXJuE5HsCuNqdduKyCwROSwiG04wXkTkKf97WSciAyqMq9XtGmDeyf6c60VkuYj0qzAu3T98jYgE/R7IAWQdKSI5Ff7e91cYd9LPkAVZ76qQc4P/c9rMP65Wt6t/nfEistj/HbVRRG6vYpra+ewaYxr8A7ADO4HOQAiwFuhVaZpRQIT/+a3AOxXG5dfBvKnAM1XM2wzY5f831v881sqslaa/DZhl4bYdDgwANpxg/CXAZ4AAQ4Hvrdiup5D3nLIcwMVlef2v04EWdWjbjgQ+PdPPUG1krTTtWOArq7arf51tgAH+51HAtiq+E2rls9tY9jQGAzuMMbuMMaXA28D4ihMYYxYbYwr9L78D2tdyxoqqzXsSFwJfGmOyjDFHgS+Bi4KUE04966+AeUHMc1LGmG+ArJNMMh6YY3y+A2JEpA21v10DymuMWe7PAxZ/bgPYtidyJp/303KKWS39zAIYY342xvzof54HbAbaVZqsVj67jaVotAP2Vni9j+M3eEU34KvYZcJEJE1EvhORy4MRsJJA807w74a+JyLxpzhvTQl4ff4mv07AVxUG1/a2rc6J3k9tb9fTUflza4D/icgqEbnZokyVnS0ia0XkMxHp7R9WZ7etiETg+4J9v8JgS7eriCQA/YHvK42qlc+u43RnbKhEZAqQDIyoMLijMWa/iHQGvhKR9caYndYkLPcJMM8YUyIitwCvA+dZnKk61wDvGWM8FYbVxW1b74jIKHxF49wKg8/1b9uWwJcissX/C9sqP+L7e+eLyCXAh0BXC/MEYizwrTGm4l6JZdtVRCLxFbA7jDG5tbHOyhrLnsZ+IL7C6/b+YccQkTHAfcA4Y0xJ2XBjzH7/v7uAJfiqfDBVm9cYk1kh4yvAwEDnrWGnsr5rqLSbb8G2rc6J3k9tb9eAiUhffJ+B8caYzLLhFbbtYeADfM1AljHG5Bpj8v3PFwBOEWlBHd62nPwzW6vbVUSc+ArGXGPMf6qYpHY+u7V5MMeqB749ql34mkbKDrT1rjRNf3wH47pWGh4LhPqftwC2E/yDdIHkbVPh+RXAd+aXg14/+XPH+p83szKrf7oe+A4gipXb1r+uBE58sPZSjj2YuNKK7XoKeTsAO4BzKg1vAkRVeL4cuMjirK3L/v74vmj3+LdzQJ+h2szqHx+N77hHkzqwXQWYAzx5kmlq5bPbKJqnjDFuEfkd8AW+MzVmGWM2ishfgTRjzMfAY0AkMF9EAPYYY8YBPYEXRcSLb89spjFmUx3I+3sRGQe48X2wU/3zZonIQ8AP/sX91Ry7a21FVvD9Ynvb+D/FfrW+bUVkHr6zeFqIyD7gL4DT/15eABbgOwtlB1AIXO8fV6vb9RTy3g80B57zf27dxtfLaSvgA/8wB/CWMeZzi7NeBdwqIm6gCLjG/3mo8jNkcVbw/Rj7nzGmoMKstb5d/VKAXwPrRWSNf9gf8f1oqNXPrnYjopRSKmCN5ZiGUkqpGqBFQymlVMC0aCillAqYFg2llFIB06KhlFIqYFo0VJ0nIvkBTHOHv8uHmlrn5SLSqwaXt/wM5s33/9tWRN47yXQxIvKb012PUoHQoqEaijuAUyoaImI/yejLgRorGsaYc2pgGQeMMVedZJIYQIuGCiotGqre8N+PYYm/g8YtIjLXfw+B3wNtgcUistg/7QUiskJEfhSR+f4+e8ruhfCoiPwITBSRm0TkB38neu+LSISInAOMAx7z3zPhLBFJ8nequE5EPhCRWP/ylojvfiFpIrJZRAaJyH/89y34W4Xs+RWe3y2++zGsFZGZVbzPTv7s6ystI0H8938Qkf/f3tmFWFVFcfz398FpaExhCIMRHPADwUGEehCkkUCEniyKIB9y1NUAAANKSURBVCY1fRWZl6wHKUjElHnTp/RBJdS+KBoi1IIwGoYZSvHOiEEx9SAJisjQSyLO6mGtq2eO98qZSRmcWT/YsDh773XWOffCumfve/5rtaThiK8maQVwAFgWx/oktclrw1wIX5sLfq5IOiqvzXBOUmv0LZf0Q8R2QdKyOL477lNN0oeP9INNniwe9+vv2bL930bU3MDf4B3HtXPmAYO4eBwUahzgkiQ/EfIPwHvAB4Vx7xZ8txfsfcCusI8Drxf6asCGsPcScg64XtbBsHuBv/HaBy24mmh76RpexqUn6rVbHpBzAPqBrWHvLMztJGQvgMNAT9jzgVZKshj4G8vPFO7JH7jERCeuJLA2+j4H3gp7CHg17Kfwp7dNwJGYOw/4Fuie6e9Ftplpc0JGJJlVDJvZVYCQU+gEfi6NWYcvLQ2E3MN8PMHU+axgd8Wv+UW4jMzZ8gklLQQWmdn5OHQC+KIwpC6VMgJcNrNrMW8MF4q7WRi7EThmUbvFGss5rAdeC/sT4GCDMYPAHklLgK/M7Pe41kmhA/sldQMTuBz24uj708zqchS/Ap2SFgAdZvZ1xPZvXMcmPHFcjPFtuDrtTCrmJjNEJo3kSeN2wb5L4++w8KIzbzbxUdQSOg68YmaXJL2NP81MN6aJUnwTTeKrwkP1fczslKQhXKTuO7k8/lhpWA/wLPC8md2R9Bf+9FCMGfw+tj7kdAI+MrOPpxB/MkvJPY1ktvAPXgYTvILdeknLASQ9LWllk3kLgGty2emeRv7MbBy4JenF6NsCnGd6fA9sr//TS1F3usQALvBIKaZ7yOuPjJnZIeAbYA2T7wG4Suv1SBgvAUsf9HQf84pwVxXFsCS1RJxngR2FfaEOeS2JZA6SSSOZLRwBzkj60cxu4Kq/pyXV8KWcVU3mvY+v4w8AvxWOfwrslnQxNoO34RvjNWAtvq8xZcwVUfuBX2J57Z0Gw3qBnZJGaF5h7Q1gNHx04WU+b+JLcqOS+oCTwAvhZ2vp+pqxBVdQruF7L8+Z2TngFDAYvr5kcnJK5hCpcpskSZJUJp80kiRJkspk0kiSJEkqk0kjSZIkqUwmjSRJkqQymTSSJEmSymTSSJIkSSqTSSNJkiSpzH/jSrr/3FkMJQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fed5ff03160>"
      ]
     },
     "metadata": {},
     "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": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no bootstrapping\n",
      "Total evaluations for no bootstrapping:\n",
      "21482\n",
      "bootstrapping\n",
      "Total evaluations for bootstrapping:\n",
      "3908\n",
      "np\n",
      "Total evaluations for np:\n",
      "0\n",
      "poly\n",
      "Total evaluations for poly:\n",
      "1411\n",
      "diff_model\n",
      "Total evaluations for diff_model:\n",
      "3351\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 != \"np\":\n",
    "            sum += result_full[key].raw_result.cost_function_evals\n",
    "        else:\n",
    "            sum = 0\n",
    "    print(sum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td><code>qiskit-terra</code></td><td>0.19.0.dev0+cdea7de</td></tr><tr><td><code>qiskit-aer</code></td><td>0.8.1</td></tr><tr><td><code>qiskit-ignis</code></td><td>0.5.1</td></tr><tr><td><code>qiskit-aqua</code></td><td>0.10.0.dev0+649b495</td></tr><tr><td><code>qiskit-nature</code></td><td>0.2.0</td></tr><tr><td><code>qiskit-finance</code></td><td>0.3.0</td></tr><tr><td><code>qiskit-optimization</code></td><td>0.3.0</td></tr><tr><td><code>qiskit-machine-learning</code></td><td>0.3.0</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.6.2 |Continuum Analytics, Inc.| (default, Jul 20 2017, 13:51:32) \n",
       "[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]</td></tr><tr><td>OS</td><td>Linux</td></tr><tr><td>CPUs</td><td>1</td></tr><tr><td>Memory (Gb)</td><td>5.827335357666016</td></tr><tr><td colspan='2'>Fri Aug 13 14:31:56 2021 EDT</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, 2021.</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",
    "\n",
    "%qiskit_version_table\n",
    "%qiskit_copyright"
   ]
  },
  {
   "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.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}