Molecular-Meta-VQE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Non-Linear Molecular Meta-VQE <img src=\"img/ML_logo.png\" width=\"20%\" align=\"right\">\n",
    "_______________________________________________________\n",
    "\n",
    "_Authors:_ Alba Cervera-Lierta, Jakob S. Kottmann, Alán Aspuru-Guzik.\n",
    "\n",
    "This notebook presents a demo code illustrate a non-linear molecular meta-VQE as introduced in the article \"The Meta-Variational Quantum Eigensolver (Meta-VQE): Learning energy profiles of parameterized Hamiltonians for quantum simulation\" (http://arxiv.org/abs/2009.13545).\n",
    "\n",
    "We will explicitly illustrate how to construct a non-linear meta-VQE for molecular systems.  \n",
    "In this Notebook we will use the simplest system available, the Hydrogen molecule in a minimal representation.  \n",
    "The molecular coordinates, or the basis-set can however be replaced (as long as there is just a single meta-parameter $R$ in the coordinates), but note that the training will take substantially longer, depending on the system/basis you chose, and a jupyter notebook is probably not the best environment anymore.\n",
    "\n",
    "In order to run the notebook you need to have `tequila` installed.  \n",
    "Just follow the intructions from the [AAG github repository](https://github.com/aspuru-guzik-group/tequila) and let us know if you experience any troubles. Feel also free to check out the `tequila` [tutorials](https://github.com/aspuru-guzik-group/tequila/tree/master/tutorials). For this notebook we recommend [basic usage](https://github.com/aspuru-guzik-group/tequila/blob/master/tutorials/BasicUsage.ipynb) and [basic chemistry](https://github.com/aspuru-guzik-group/tequila/blob/master/tutorials/Chemistry.ipynb).\n",
    "\n",
    "We start by defining our non-linear encoder. As in the paper we will use a set of floating Gaussians to approximate our unknown non-linear function. Each parameter of the VQE is then encoded, depending on the meta variable $R$ (in this example the bond distance of H$_2$; in the paper the intermolecular distance of H$_4$), as\n",
    "$$\n",
    "\\theta\\left(R\\right) = \\sum_n \\alpha_n \\exp\\left( -\\beta_n (\\gamma_n - R)^2\\right) + \\delta.\n",
    "$$\n",
    "\n",
    "In the following cell, we define this encoder as abstract class can produce `tequila` objectives which will be used as parameters for quantum circuits later.\n",
    "\n",
    "The `DummyEncoder` class just gives back the original variable, we will use it to test the meta-VQE later.  \n",
    "In the same way the `evaluate` function will later be used to initialize the $\\theta$ parameters to test the meta-VQE after training.\n",
    "\n",
    "If you have any questions left, feel free to contact us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# tested with tq git revision: ab252adc4090465ff47f182de32d69b84d7fc13f\n",
    "import tequila as tq\n",
    "import numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class DummyEncoder:\n",
    "    \"\"\"\n",
    "    No encoding\n",
    "    \"\"\"\n",
    "    def __call__(self, variable, *args, **kwargs):\n",
    "        return variable\n",
    "        \n",
    "class GaussianEncoder:\n",
    "    \"\"\"\n",
    "    Represent an unknown linear function by n floating Gaussians    \n",
    "    \"\"\"\n",
    "    def __init__(self, n):\n",
    "        self._n = n\n",
    "\n",
    "    def __call__(self, variable, R):\n",
    "        # this allows us to just pass names as parameters\n",
    "        variable = tq.assign_variable(variable)\n",
    "        \n",
    "        # here we create the non-linear encoded parameter\n",
    "        # as described in the paper\n",
    "        # as abstract tequila objective\n",
    "        mapped = tq.Variable((\"delta\", variable))\n",
    "        for n in range(self._n):\n",
    "            a = tq.Variable((\"alpha\", variable, n))\n",
    "            b = tq.Variable((\"beta\", variable, n))\n",
    "            c = tq.Variable((\"gamma\", variable, n))\n",
    "            mapped += a * ((-(b * (c - R) ** 2)).apply(tq.numpy.exp))\n",
    "        return mapped\n",
    "    \n",
    "    def evaluate(self, variable, R, trained_variables):\n",
    "        # evaluates a specficic variale at point R\n",
    "        # according to alpha*exp(-beta*(gamma - variable)**2) + delta\n",
    "        # needs the trained variables which are (\"alpha\", variable, n)\n",
    "        # and similar for the other meta parameters\n",
    "        result = 0.0\n",
    "        for n in range(self._n):\n",
    "            a = trained_variables[(\"alpha\", variable, n)]\n",
    "            b = trained_variables[(\"beta\", variable, n)]\n",
    "            c = trained_variables[(\"gamma\", variable, n)]\n",
    "            d = trained_variables[(\"delta\", variable)]\n",
    "            result += a*numpy.exp(-b*(c-R)**2) + d\n",
    "        return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The VQE Ansatz\n",
    "\n",
    "Now its time to define the ansatz for our molecular VQE.  \n",
    "As in the paper, we will use the [UpCCGSD](https://pubs.acs.org/doi/10.1021/acs.jctc.8b01004) ansatz.\n",
    "\n",
    "The quantum circuit for this ansatz has the form \n",
    "$$\n",
    "U(\\boldsymbol{\\theta}) = \\left(\\prod_k e^{-i\\frac{\\theta_k}{2} G_k}\\right) U_\\text{ref} \n",
    "$$\n",
    "where the generators $G_k$ are qubit encoded products of fermionic creation and anihialtion operators; In this case, single excitations and pair-restricted double excitations of the form\n",
    "$$\n",
    "G_{pq} = i(a_p^\\dagger a_q - a^\\dagger_q a_p)\n",
    "$$\n",
    "and\n",
    "$$\n",
    "G_{pqrs} = i(a_p^\\dagger a_q a_r^\\dagger a_s - h.c.).\n",
    "$$\n",
    "The pair-restriction in the doubles requires the spin orbitals $p,q$ and $r,s$ to have the same spatial part. This is enforced in the code block below by iterating over spatial indices and generating $p,q$ and $r,s$ from them (`tequila` will enumerate spin-up with even and spin-down with off numbers).\n",
    "\n",
    "The qubit encoded generators are sums of tensor-products of Pauli matrices. Their explicit form depends on the qubit encoding chosen below (the default here is `jordan_wigner`). If you want to know the explicit form of some generators, just construct them as illustrated below, and then print them out. Note that the qubit encoding is set over the `transformation` key when the molecule object is created (see cells further below).\n",
    "\n",
    "The $U_\\text{ref}$ part of the circuit creates the Hartree-Fock reference. The explicit form also depends on the chosen qubit encoding and the `tequila` molecule object handles the correct initialization. In our specfic H$_2$\\STO-3G example below in `jordan_wigner` representation, the $U_\\text{ref}$ circuit consists of two $X$ gates on qubits 0 and 1 (the spin-up and spin-down encoded occupied HF orbital).\n",
    "\n",
    "Note that for the Hydrogen system we can neglect the singles for point-group symmetry reasons, and for the same reasons will only end up with one pair-ecitation for the doubles anyway.  \n",
    "The code-block below will however generate the ansatz independent of the molecular that is passed to it.  \n",
    "Note that the encoder is also passed to the function, so in principle you can also use other encoders as the ones defined above; They just need to have a call operator defined which takes the name of the variable to encode and the meta paramerer $R$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def encoded_upccgsd(molecule, encoder, R, layers=1, include_singles=True):\n",
    "    U = molecule.prepare_reference()\n",
    "    for layer in range(layers):\n",
    "        for i in range(molecule.n_orbitals):\n",
    "            for a in range(i+1,molecule.n_orbitals):\n",
    "                # First the pair-restricted doubles\n",
    "                # create the encoded angle objective\n",
    "                angle = encoder(\"D_{}_{}_{}\".format(layer, i, a), R)\n",
    "                G = molecule.make_excitation_generator(indices=[(2*i,2*a),(2*i+1,2*a+1)])\n",
    "                U += tq.gates.Trotterized(generators=[G], angles=[angle], steps=1)\n",
    "                if not include_singles:\n",
    "                    continue\n",
    "                # now the same for the singles; we will use the same parameter\n",
    "                # for spin-up and spin-down singles to enforce singlet to enforce singlet symmetry\n",
    "                angle = encoder(\"S_{}_{}_{}\".format(layer, i, a), R)\n",
    "                G1 = molecule.make_excitation_generator(indices=[(2*i,2*a)])\n",
    "                G2 = molecule.make_excitation_generator(indices=[(2*i+1,2*a+1)])\n",
    "                U += tq.gates.Trotterized(generators=[G1,G2], angles=[angle,angle], steps=1)\n",
    "    return U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Construction and Training of the meta-VQE  \n",
    "\n",
    "Now we will construct our meta-VQE objective by summing up expectationvalues with meta parameters defined by a set of training points.  \n",
    "\n",
    "After the `tequila` objective is created we minimize it with the `BFGS` optimizer of the `scipy` package.  \n",
    "If you want to try another optimizer, see the tequila [tutorials](https://github.com/aspuru-guzik-group/tequila/blob/master/tutorials/Optimizer_Tutorial.ipynb) or type `tq.show_available_optimizers()` to see which optimizers are available on your system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initial values are:\n",
      " {('alpha', D_0_0_1, 0): 0.0, ('delta', D_0_0_1): 0.0, ('gamma', D_0_0_1, 0): 1.0, ('beta', D_0_0_1, 0): 1.0}\n",
      "Optimizer: <class 'tequila.optimizers.optimizer_scipy.OptimizerSciPy'> \n",
      "backend         : qulacs\n",
      "device          : None\n",
      "samples         : None\n",
      "save_history    : True\n",
      "noise           : None\n",
      "\n",
      "Method          : BFGS\n",
      "Objective       : 5 expectationvalues\n",
      "grad instr      : 2-point\n",
      "gradient        : scipy numerical 2-point\n",
      "hessian         : scipy numerical None\n",
      "\n",
      "active variables : 4\n",
      "\n",
      "Optimization terminated successfully.\n",
      "         Current function value: -1.010502\n",
      "         Iterations: 13\n",
      "         Function evaluations: 110\n",
      "         Gradient evaluations: 22\n",
      "Finished Training of Meta VQE\n",
      "optimized meta-variables are:\n",
      " ('alpha', D_0_0_1, 0) : -1.3507766463040127\n",
      "('delta', D_0_0_1) : 1.5022768577876158\n",
      "('gamma', D_0_0_1, 0) : 0.6887929726422215\n",
      "('beta', D_0_0_1, 0) : 0.8901106480850393\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# define the encoder\n",
    "encoder = GaussianEncoder(n=1)\n",
    "\n",
    "# define details on the VQE ansatz\n",
    "layers = 1\n",
    "include_singles = False\n",
    "transformation = \"jordan_wigner\"\n",
    "geometry = \"H 0.0 0.0 0.0\\nH 0.0 0.0 {R}\" # gan also give path to an xyz file\n",
    "basis_set = \"sto-3g\"\n",
    "\n",
    "# define the training points\n",
    "training_points = [0.5, 0.7, 1.5, 3.0, 5.0]\n",
    "\n",
    "# define the test points\n",
    "start = 0.3\n",
    "end = 7.0\n",
    "steps = 25\n",
    "test_points = [start + step/steps*(end-start) for step in range(steps)]\n",
    "\n",
    "# Details on the used optimizer:\n",
    "optimizer_arguments = {}\n",
    "# see tq.show_available_optimizers() for other options\n",
    "optimizer_arguments[\"method\"] = \"bfgs\"\n",
    "# increase to see more output in the optimization\n",
    "optimizer_arguments[\"print_level\"] = 1\n",
    "optimizer_arguments[\"gradient\"] = \"2-point\"\n",
    "# eps for scipy < 1.5, and finite_diff_rel_step for scipy > 1.5\n",
    "optimizer_arguments[\"method_options\"] = {\"finite_diff_rel_step\":1.e-5, \"eps\":1.e-5, \"gtol\":1.e-5}\n",
    "\n",
    "\n",
    "#construct the tequila objective by accumulating expectation values of training points\n",
    "objective = 0.0\n",
    "for R in training_points:\n",
    "    molecule = tq.chemistry.Molecule(geometry=geometry.format(R=R), basis_set=basis_set, transformation=transformation)\n",
    "    H = molecule.make_hamiltonian()\n",
    "    U = encoded_upccgsd(molecule, encoder, R, layers, include_singles)\n",
    "    \n",
    "    E = tq.ExpectationValue(H=H, U=U)\n",
    "    objective += E\n",
    "\n",
    "objective = 1.0/len(training_points)*objective\n",
    "    \n",
    "# initialize in a way that the encoded angles start with zero\n",
    "# initialize beta and gamma as 1.0, since they are width and shift of the Gaussians\n",
    "initial_values = {k:(1.0 if \"beta\" in k.name[0] or \"gamma\" in k.name[0] else 0.0) for k in objective.extract_variables()}\n",
    "\n",
    "# if we have more than one floating Gaussian, we should not initialize the shifts (gamma) to the same values\n",
    "# otherwise the optimizer can get confused\n",
    "if encoder._n > 1:\n",
    "    initial_values = {k:(numpy.random.uniform(training_points[0], training_points[-1], 1)[0] if \"gamma\" in k.name[0] else initial_values[k]) for k in objective.extract_variables()}\n",
    "\n",
    "print(\"initial values are:\\n\", initial_values)\n",
    "    \n",
    "# minimize the objective formed from the training points\n",
    "trained = tq.minimize(objective=objective, initial_values=initial_values, **optimizer_arguments)\n",
    "\n",
    "print(\"Finished Training of Meta VQE\")\n",
    "print(\"optimized meta-variables are:\\n\", trained.variables)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing the performance of the meta-VQE\n",
    "\n",
    "Here we gonna test the performance of the meta-VQE by computing the expectation values of different points with the parameters of the quantum circuit initialized by the trained meta-VQE. For this we will need the `evaluate` function of the encoder. Furthermore we will compute the `fci` energy of each test point, and optimize the individual points with a standard VQE initialized with the meta-VQE values (opt-meta-VQE).  \n",
    "See the paper for further details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "processing test point +0.3000\n",
      "processing test point +0.5680\n",
      "processing test point +0.8360\n",
      "processing test point +1.1040\n",
      "processing test point +1.3720\n",
      "processing test point +1.6400\n",
      "processing test point +1.9080\n",
      "processing test point +2.1760\n",
      "processing test point +2.4440\n",
      "processing test point +2.7120\n",
      "processing test point +2.9800\n",
      "processing test point +3.2480\n",
      "processing test point +3.5160\n",
      "processing test point +3.7840\n",
      "processing test point +4.0520\n",
      "processing test point +4.3200\n",
      "processing test point +4.5880\n",
      "processing test point +4.8560\n",
      "processing test point +5.1240\n",
      "processing test point +5.3920\n",
      "processing test point +5.6600\n",
      "processing test point +5.9280\n",
      "processing test point +6.1960\n",
      "processing test point +6.4640\n",
      "processing test point +6.7320\n"
     ]
    }
   ],
   "source": [
    "fci = {}\n",
    "meta = {}\n",
    "meta_angles = {}\n",
    "meta_opt = {}\n",
    "vqe = {}\n",
    "# we are using this to be able to use the same codeblocks\n",
    "dummy = DummyEncoder()\n",
    "for R in test_points:\n",
    "    print(\"processing test point {:+2.4f}\".format(R))\n",
    "    molecule = tq.chemistry.Molecule(geometry=geometry.format(R=R), basis_set=basis_set, transformation=transformation)\n",
    "    H = molecule.make_hamiltonian()\n",
    "    U = encoded_upccgsd(molecule, dummy, R, layers, include_singles)\n",
    "    E = tq.ExpectationValue(H=H, U=U)\n",
    "    \n",
    "    # evaluate regular VQE with meta VQE parameters\n",
    "    values = {k:encoder.evaluate(variable=k, R=R, trained_variables=trained.variables) for k in E.extract_variables()}\n",
    "    energy = tq.simulate(E, variables=values)\n",
    "    meta[R] = energy\n",
    "    meta_angles[R] = values\n",
    "    \n",
    "    # optimize from meta VQE starting point\n",
    "    result = tq.minimize(objective=E, initial_values=values, silent=True, **optimizer_arguments)\n",
    "    meta_opt[R] = result.energy\n",
    "    \n",
    "    # try fci (you will need psi4 for this)\n",
    "    try:\n",
    "        fci[R] = molecule.compute_energy(\"fci\")\n",
    "    except:\n",
    "        if molecule.n_orbitals < 5:\n",
    "            fci[R] = numpy.linalg.eigvalsh(H.to_matrix())[0]\n",
    "        else:\n",
    "            print(\"skipping fci for point \", R)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The collected data of the test above is plotted here.  \n",
    "Since we only have one parameter in our simple guinny pig system here, we will also plot the quantum circuit's parameter with respect to the meta parameter $R$.  \n",
    "If you have changed the molecule or the basis_set above, you might need to adapt the plot for the parameter in order to see the parameter you are interested in.\n",
    "\n",
    "For our guinny-pig system and the default settings in this notebook you see that the meta-VQE predicts the whole PES accurately with exception of the point at very close distance. The reason can be seen in the plot of the parameter. For this specific system and this specfic VQE ansatz, the ideal parameter will have a sigmoidal dependence on the meta parameter $R$ and the single floating Gaussian is not able to perfectly adapt to it. As explained in the paper: Chosing floating Gaussians was a pragmatic choice since they are fairly good approximators for generalized functions.\n",
    "\n",
    "\n",
    "If you are interested you can try to change the encoder above, to a sigmoidal function, then you will be able to fit this test system almost perfectly. Another way to improve is to increase the number of floating Gaussians. \n",
    "Feel free to play around with it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU9b3/8ddnJpNkEpJMWGQVocUiiwga8VJbF8CCW0GtIte2aGu1tt5frUt/Uqu2d1Fard1vvdaKWKvYKihULIr1d12qYhAUwSJoERNWwWxkmUzm+/tjJgiSkIRJciZz3s/HI4/MTL7nfD+J+M7Jd875HHPOISIi/hDwugAREek+Cn0RER9R6IuI+IhCX0TERxT6IiI+kuV1AYfSt29fN2zYMK/LEBHpMVatWvWhc65fa19P69AfNmwYpaWlXpchItJjmNn7h/q6lndERHxEoS8i4iMKfRERH0nrNX0RyQyNjY2UlZVRX1/vdSkZIzc3lyFDhhAKhTq0nUJfRLpcWVkZBQUFDBs2DDPzupwezznH7t27KSsrY/jw4R3aNuNC//HV5dyxfANbK+oYFAlzw7SRzJww2OuyRHytvr5egd+JzIw+ffqwa9euDm+bUaH/+Opy5i5aS11jEwDlFXXMXbQWQMEv4jEFfuc63J9nRr2Re8fyDfsCv1ldYxN3LN/gUUUiIuklo0J/a0Vdh14XEWmPNWvWsGzZsg5ts2DBAmbPnn3Aax9++CH9+vWjoaGBaDTKNddcw4gRIxgxYgTnnHMOW7Zs2Tc2GAwyfvz4fR/z5s3rlO8lo5Z3BkXClLcQ8IMiYQ+qEZHDlW7vza1Zs4bS0lLOOuusdm9z3nnncd1111FbW0teXh4Ajz76KOeeey45OTlcf/31VFdXs2HDBoLBIPPnz2fGjBmsWrWKQCBAOBxmzZo1nf69ZNSR/g3TRhIOBQ94LRwKcsO0kR5VJCId1fzeXHlFHY6P35t7fHV5SvvdvHkzxxxzDJdeeimf+cxnuOSSS1ixYgUnn3wyRx99NCtXrmTv3r187WtfY+LEiUyYMIEnnniCaDTKLbfcwiOPPML48eN55JFHWLlyJZMmTWLChAl89rOfZcOGg5eQCwsLOfXUU1m6dOm+1xYuXMjs2bOpra1l/vz5/OxnPyMYTGTWZZddRq9evVixYkVK32dbMupIv/lI4HuPvUk0FmdwGhwhiMiBfrR0Heu3VrX69dVbKog2xQ94ra6xie89+iYPr9zS4jajBxVy67lj2px706ZN/PnPf+a+++7jxBNP5KGHHuLFF19kyZIl3HbbbYwePZrJkydz3333UVFRwcSJE5k6dSr//u//TmlpKb/+9a8BqKqq4oUXXiArK4sVK1bw/e9/n8cee+yg+WbPns0f//hHZs2axdatW3nnnXeYPHky69atY+jQoRQWFh4wvqSkhPXr1/OFL3yBuro6xo8fv+9rc+fOZdasWW1+j23JqNCHRPAvX7edTTtreObaU70uR0Q66JOB39brHTF8+HCOPfZYAMaMGcOUKVMwM4499lg2b95MWVkZS5Ys4c477wQSp5ruv87erLKykjlz5rBx40bMjMbGxhbnO/vss/nWt75FVVUVf/rTn7jgggv2Hdm3pauWdzIu9AGKwiEq61r+jyAi3mrriPzkeX9r8b25wZEwj1w5KaW5c3Jy9j0OBAL7ngcCAWKxGMFgkMcee4yRIw9cEn711VcPeH7zzTdz+umns3jxYjZv3sxpp50GJJZoVq9ezaBBg1i2bBnhcJjp06ezePFiFi5cyF133QXApz/9abZs2UJ1dTUFBQX79rtq1SouuOCClL7HtmTUmn6zonCICoW+SI/k5Xtz06ZN41e/+hXOOQBWr14NQEFBAdXV1fvGVVZWMnhwYtn4/vvv3/f6/PnzDzrTZ/bs2dx1113s2LGDSZMSv7Ty8/OZM2cO1157LU1NidPMH3jgAXJzczn55JO79HvMyNAvDIeIxuLUf+KcfRFJfzMnDOb2849lcCSMkTjCv/38Y7vlvbmbb76ZxsZGxo0bx5gxY7j55psBOP3001m/fv2+N3K/973vMXfuXCZMmEAsFjvkPs844wy2bt3KrFmzDrig6vbbbyccDjNy5EgGDx7MXXfdxRNPPLFvTPOafvPHjTfe2CnfozX/RktHJSUl7nBuovLHV9/npsVv8er3p9C/MLcLKhORjnj77bcZNWqU12Wkre3bt3PmmWdy1VVXccUVV7R7u5Z+rma2yjlX0to2GbumD1BZ16jQF5G0N2DAgH1LSV0tpeUdM+ttZs+Y2cbk5+JWxg01s6fN7G0zW29mw1KZty3NoV9Rq3V9EZH9pbqmfyPwrHPuaODZ5POWPADc4ZwbBUwEdqY47yHtf6QvIiIfSzX0ZwALko8XADM/OcDMRgNZzrlnAJxzNc652hTnPaRIOBtQ6IuIfFKqa/r9nXPbko+3A/1bGPMZoMLMFgHDgRXAjc65rjm1xjkGLZrJvwWPJFxWA3Xb4HPXdMlUIiI9TZuhb2YrgAEtfOmm/Z8455yZtXQqUBbweWACsAV4BLgU+H0r810BXAEwdOjQtspraQcEd67lutCrsPpRmLO07W1ERHyizeUd59xU59zYFj6eAHaY2UCA5OeW1urLgDXOufecczHgceD4Q8x3j3OuxDlX0q9fv8P6piynoO1BIuIrn2xVvHnzZgBWrlzJKaecwsiRI5kwYQKXX345tbW13H///Vx99dXeFt0FUl3eWQLMAeYlPz/RwpjXgIiZ9XPO7QImAx0/+b4jLMh6O5q/DfkWV5e/DsNP6dLpRCT9tdTLZseOHVx44YUsXLhw39Wyjz766AFX32aaVN/InQecYWYbganJ55hZiZndC5Bcu78eeNbM1gIG/C7FeQ/tiFFkBYzSwFit54tIq37zm98wZ86cfYEP8KUvfYn+/Vt6ezIzpHSk75zbDUxp4fVS4PL9nj8DjEtlrg7J600R63X2jkg6eupG2L627XGVZZDTK/G4oQaKhrQ+dsCxcOah7yy1f6vi4cOHs3jxYt566y3mzJnT3sozQkZekUu4NwWuWqEv0pPl9IIdbyUe9x+b8u66qlVxT5OhoV9MXryGvbX1XlciIp/UxhH5Pv98Hhacm3g8/fYueW9uzJgxrFq1ihkzZnT6vtNVRnbZJK83APG6CtK5oZyIHEL564lTrucsTTzuAldffTULFiw4oF/+okWL2LFjR5fMlw4y9kgfoMBVUxttIj8nM79NkYy2/0kYXXQGXv/+/Vm4cCHXX389O3fuJBAIcMoppzB9+vQumS8dZGYahhNH+hFqqKxrVOiLCDU1NS2+PmnSJF544YWDXr/00ku59NJLu7iq7pehyzuJI/2I1ejNXBGR/WRm6CeXd4qpUXtlEZH9ZGjoJ5d3TKdtiqQLnVTRuQ7355mZoZ9bhLMgEdtLlUJfxHO5ubns3r1bwd9JnHPs3r2b3NyO3xkwM9/hNMPlRihu1JG+SDoYMmQIZWVl7Nq1y+tSMkZubi5DhhziKuVWZGboA5bXm+KaGrbVRb0uRcT3QqEQw4cP97oMIVOXd0iEft9grY70RUT2k7GhT7iY3oEaKutiXlciIpI2Mjj0exOhhopaLe+IiDTL3NDPS3Ta1Nk7IiIfy9zQD0fIdfXU1tZ6XYmISNrI4NBPXKBF3Ufe1iEikkYyOPQTrRiCDR8Rj+uCEBERyOTQT/bUL3I11ER1Bo+ICGRy6O/rv1NDpZquiYgAGR36aq8sIvJJmRv6yeWdYtR/R0SkWeaGfiiPeDCHiO1V6IuIJGVu6Cc7bUZ0pC8isk/mhj7JTpta0xcR2ccXoa9bJoqIJKQU+mbW28yeMbONyc/FLYw53czW7PdRb2YzU5m33fWFi+kT0Jq+iEizVI/0bwSedc4dDTybfH4A59xzzrnxzrnxwGSgFng6xXnbJ683EatR0zURkaRUQ38GsCD5eAHQ1hH8l4CnnHPd0wUt3JtCV01FbUO3TCciku5SDf3+zrltycfbgf5tjL8YePhQA8zsCjMrNbPSlO+nGS4mRIyGuurU9iMikiHavEeuma0ABrTwpZv2f+Kcc2bWamczMxsIHAssP9R8zrl7gHsASkpKUuuUlrxAi9o9Ke1GRCRTtBn6zrmprX3NzHaY2UDn3LZkqO88xK4uAhY757pvgT3ZiiFQX9FtU4qIpLNUl3eWAHOSj+cATxxi7GzaWNrpdMmma6FoBU1qrywiknLozwPOMLONwNTkc8ysxMzubR5kZsOAI4H/TXG+jkku70TYS3W9zuAREWlzeedQnHO7gSktvF4KXL7f883A4FTmOizJ5Z1iS7RiiORld3sJIiLpJKOvyN3XUx+1YhARgUwP/axsmrLyiagVg4gIkOmhD8RzI2q6JiKSlPGhT15vLe+IiCRlfOgH8nvrlokiIkkZH/rBvN70VuiLiAA+CH3CyRup6I1cEREfhH5ebwrYS5U6bYqI+CD0w8UEidNYq/47IiI+CP3EBVpxddoUEfFD6Dd32vzI40JERLyX+aGfbLqW1aDlHRGRzA/95PJOTmMljU1xj4sREfGWD0K/udOmbpAuIuKD0I/gMPXfERHBD6EfCBLLLqRI/XdERHwQ+kBTTqLTZoVCX0R8zhehT7g3xVRrTV9EfM8XoR/I702R7dXyjoj4ni9CPys/caSvpmsi4ne+CP1Afh+t6YuI4JPQJ9ybAqujurbO60pERDzlk9BPXKAV26umayLib/4I/WT/HafQFxGf80fohyOJz3XqtCki/uaT0G/utKnQFxF/Szn0zay3mT1jZhuTn4tbGfcTM1tnZm+b2S/NzFKdu92SyzuhqNori4i/dcaR/o3As865o4Fnk88PYGafBU4GxgFjgROBUzth7vZJvpGb11RNQ6yp26YVEUk3nRH6M4AFyccLgJktjHFALpAN5AAhYEcnzN0+OYXELYtiq9ZVuSLia50R+v2dc9uSj7cD/T85wDn3MvAcsC35sdw593ZLOzOzK8ys1MxKd+3a1QnlAWZEs4uIsFf9d0TE17LaM8jMVgADWvjSTfs/cc45M3MtbD8CGAUMSb70jJl93jn3wifHOufuAe4BKCkpOWhfh6spJ0KkVkf6IuJv7Qp959zU1r5mZjvMbKBzbpuZDQR2tjDsPOAV51xNcpungEnAQaHfVVy4mOKPaqhQ/x0R8bHOWN5ZAsxJPp4DPNHCmC3AqWaWZWYhEm/itri801UsrzcRddoUEZ/rjNCfB5xhZhuBqcnnmFmJmd2bHPMo8C6wFngDeMM5t7QT5m63rPzeRPRGroj4XLuWdw7FObcbmNLC66XA5cnHTcCVqc6VilCvvhTrloki4nP+uCKXxI1UwhalpqbG61JERDzjm9BvvkCrqWa3x4WIiHjHR6GfaMUQr1WnTRHxLx+FfrIlkDptioiP+Sf0k03XAvUKfRHxL/+Efri506ZCX0T8y0ehn1jeyY5WeVyIiIh3/BP62XnEAjkUuCrqG9VeWUT8yT+hD0RDiU6b6r8jIn7lq9CP5UTUU19EfM1XoR8PFxMxtWIQEf/yVehbuDcRaqiojXpdioiIJ3wV+sH83hTrSF9EfMxXoZ/Vqy8RaqjUkb6I+FTKrZV7kuyCvgSsibq9lV6XIiLiCV8d6QfyEhdoNVar06aI+JOvQr+5/05TrUJfRPzJX6Hf3GmzVv13RMSffBb66rQpIv7mr9BPLu9kKfRFxKf8Ffq5EQCyG3X2joj4k79CPyubhmAeubFKnHNeVyMi0u38FfpAQyhCETXURtVeWUT8x3ehH8suSvTfUSsGEfEh34V+PLc40X9HPfVFxId8F/ourzdFqOmaiPhTSqFvZr3N7Bkz25j8XNzKuB+b2VvJj1mpzJmqYF5zp001XRMR/0n1SP9G4Fnn3NHAs8nnBzCzs4HjgfHAScD1ZlaY4ryHLatXH4rYS1VtvVcliIh4JtXQnwEsSD5eAMxsYcxo4HnnXMw5txd4E5ie4ryHLdFp01FXrQu0RMR/Ug39/s65bcnH24H+LYx5A5huZnlm1hc4HTiytR2a2RVmVmpmpbt27UqxvIPlFPQBIFr9YafvW0Qk3bXZT9/MVgADWvjSTfs/cc45Mzvoiifn3NNmdiLwd2AX8DLQ6knyzrl7gHsASkpKOv0KKstLhH5TjTptioj/tBn6zrmprX3NzHaY2UDn3DYzGwjsbGUf/wX8V3Kbh4B3DrPe1CX778TVaVNEfCjV5Z0lwJzk4znAE58cYGZBM+uTfDwOGAc8neK8hy/ZXlmdNkXEj1K9XeI84E9m9nXgfeAiADMrAb7pnLscCAEvmBlAFfBl51wsxXkPXzL0sxoU+iLiPymFvnNuNzClhddLgcuTj+tJnMGTHnKLiGOEouq0KSL+47srcgkEqQ8WkKv2yiLiQ/4LfaA+VEReUyXxuNori4i/pLqm3yPFsiMU7d1LTTRGYW7I63LEI4+vLueO5RvYWlHHoEiYG6aNZOaEwZ5vo7pUV0e26Sh/hn5uhGL7gMraRoV+N0jH/ykfX13O3EVrqWtMXDJSXlHH3EVrATzdRnWpro5sczgsne8gVVJS4kpLSzt9v+Xzv0r8n3+n8spVjB1c1On7l4998h8yQDgU5Pbzj233P/79tzl33EAaY400RuuJRaPEGuuJNUb53/XlLHjhHVxTI0HiGHFyg8YFxw9g3KBC4i6Oa4rjXAzXFOd/nt9EbX2UswIv874bADg+Zdt5KTSJC08YQhwHzuGcw7k4zjmWrtlKXTTGGYFS3neJi8+Psh28lDWRM0YdgYPkNuASz3huw07qo3GmBFaxxR0BwFDbyUvBEzl5RB8s+f+fS368/O6H1DfGmRxY/YnxJ3DS8D4H3PGt+eFrm/fQEGthm8AJnHBU7+SePx4P8PqWj1rdZsLQyEHj13xQ0er4446MtPjfcc0HFURb2ObFQ2zzxqG2GXLwNm+UtT5+XAvjAd7shm0OOb6VzHmzvHLfNu+6QbzmjmGcvcdfCi7ipRsnt7hNS8xslXOupNWv+zH0tz3yHXqtf4Q3v/IWJ4/o2+n7z3QdOaI+ed7fKK+oI5cGBtoeCtlLgdXRP7uB6SPCWEMV1lBNMFpNqLGaUFMN8doKelHLYHbhDLJoIocYYISIETj4wm+RjDQ7ehOvxMfwz3lnt3ubtkLfl8s7Wb36UGB1VNbsBRT6HdHyn6BvsmdvA6P7hago+wf1OzbBnvfIq36fn9ZtYVjOdgZYC9dFbEp8imPUEqY2kE9dIJ8d5LDDFVNBPp8LrAPgyaaJbHH9OX54fwiGsKxsCIQgKxsLZvPwqu00ksVw28a1occA+HHjLN52R3HD9NFYIIAFggSSn3+8fCMf1TYy2jbzn9n3AzA3+jXK80bx4y8dRyBgGJYYb2AW4OqHV7OrppFxtomfZt8DwHXRK9maP4b7v3YiiWtRjICBmREw419/9yo7qusZb5u4K/tuAL4bvYptvcaw8IpJiR9A4hoWAGb9z8v7xv88+7cAXBO9im29xvLIlZNa/G/SHduoLm/qAhgUCbc4/nD5MvSzCxJBX1+1GzjK22J6mJ/89R9Y417uzJpPlCz620eMsK3kPNN4ULBXBiJssv68GD+WmAtwcdb/A+DW6Fd4N/94FnxrKsFwEYHsAnoFAvRKbvevyb8OJgXW8blgIvQfbDqDLYUlXPX1lv/MvWbjx9s0W+NGsKWwhDGnHLzNuaGxzF20lhPjG5gdTbSROiHrn5x09tkMPKblv1ouOusI5i5ay+T46gO2+fxZU8kZ2PI2s87MZ+6itZwRX3XANqeeORn6HrzN7DNzmbtoLf3jpQeMP+3M06FPy3N0xzaqy7u6Zk27pMXxh8uXoZ+bDP1otZqutUd9YxMvr/snW1cu4oe1f+WUnDfJtY/vPLYhPoTX4sdw4vEl5A/8DJEhx5DV99MU5Rbywepybl60lvHxN7mY/wfA5uAwvnTWNILFLf/jv2HaSOYuWsu4+Hvt/sff0W2al6PuWJ69b5lq1rRLDvmGWXdso7pUV0e2ORy+XNN3m/6GPXgeD466my/Pmt3p++9JWlufr43G+PvajWx/bTFDtj3DZ3mDbGtiu+vNU00n8oHrxy2hB4HEuuOWwpJW32x6fHU55U/O48XaI+mbn8OVR1cy9qJbDquuzt5GJNPojdyWbF0D95zKA0fdxlcv+3bn77+HaOksmf7BKr5a/BbHVT/PSawjZE18FBpAzafPZsC/XMSTewYzd/E6vhp/nDfdp4DEEfWI836ggBVJA3ojtyXJpmtW5++ma3cs30BdYxN9qeTM4KucGVjJSYG3CdY4ducOYceIbzDgXy6ieMjxFCffbJw5DLBAl/8JKiJdw5+hn+ypH/R5p82tFbXMDv6NH2YtIMdilMX78nJ8NLfFvsyyW6864KyS/c2cMFghL9JD+bL3Dtm9iBEk1FDhdSWeeaF0DQ9kz+P20O95Jz4EgCGBD/lN00wqi45pNfBFpGfz55G+GbXBInJ82Gmzui7Ksj/8lDPLf0FWwHFr09d5p6k/DwdvAyA7GOCGaSM9rlJEuoo/Q59Ep81wvb9C/7W162hc/G/Miq9iS9HxDPjK75lQlsMRT85jdu1N9M3P4YbRlYzV0o1IxvJt6Eezi+i1t4amuCMYyOyljLqGGE899HOmbP4pORbjg3/5IUO/8B0IBJjZD5jwK/x7DpOIv/g29JtyIkTsXarrG4nkZXtdTpdZ+48NVP35as5vWsmWXuM44iv3ceSAo70uS0Q84tvQj+f2ptjeoKI2c0J//4uTBhbmcGnR61y48xfkWZR/Hv99hp9zPQSCXpcpIh7ybehbXjERathQ19j24B5g/wut+lDJD+rmc1Z0Je9kj2TQpfczfHD63KZYRLzj29DPyu9DrjVSXV0FtNxDuydpvtDqtqzfcXbwFcI08Lem47g152ZeUOCLSJI/z9MHQgV9AKir+tDjSjrH1oo6Rtn7nBd8kSKrI9vi/K7pHMoqo16XJiJpxLeh39xpsyFDQn90rxruy76DveQe8Hpn9+IWkZ7Nt8s74Ug/AGI1PT/0d+/ZzZ2x2ymglkeaTmNF/ASga3pxi0jP5tvQb76RSry2Z/ffqW+I8s+7L2aC28zSMT/jvveOUiM0EWlVSqFvZhcCPwRGAROdcy32QTaz6cAvgCBwr3NuXirzdopkp01X23NvpOKc49W7v8mp0ZWsnXArM2dexkyvixKRtJbqmv5bwPnA860NMLMg8BvgTGA0MNvMvD+dJBn6wfqee6T/v3/4D0796DHWDLmEY2de63U5ItIDpHSk75x7G0jeELpVE4FNzrn3kmMXAjOA9anMnbJQmAZyyGromf13XnnqQT7/7l2sLfw8x132S6/LEZEeojvO3hkMfLDf87Lka57bGywkuwd22ly/6nnGvXItm7NHMPKqh7Ggb9+aEZEOajMtzGwFMKCFL93knHuiswsysyuAKwCGDh3a2bs/QF1WEeHGntVTv2zzO/Rb+lWqAoX0/cZisvMKvC5JRHqQNkPfOTc1xTnKgSP3ez4k+Vpr890D3AOJe+SmOPchRbOLyK+v7sopOlXlR3toeOBCiqin8uInKTriyLY3EhHZT3cs77wGHG1mw80sG7gYWNIN87Yplh2h0FXT2BT3upQ2RaNR3rv7Io5q2kL51LsZMvIEr0sSkR4opdA3s/PMrAyYBDxpZsuTrw8ys2UAzrkYcDWwHHgb+JNzbl1qZXeOeLiYiNVQleZN11w8zqq7v8GEhtd447hbOOZzOjFTRA5PqmfvLAYWt/D6VuCs/Z4vA5alMleXCCc6bb5fG6VPrxyvqzlIc6vk6dWPcnPocZYXzWLa+d/1uiwR6cF823sHIJjfhyyLU121x+tSDtLcKnls1fPclPUQy5om8t09M3h8datvh4iItMnXoZ/VK9lpszL9+u/csXwDQ2Lv8+vQL3nXDWRh0+l8Jb6UO5Zv8Lo0EenBfB36OYWJ/jvRNOy0ua1iL/NCv6OOEEcHtvJA9o95032KrRV1XpcmIj2Yr0M/XJQM/er0C/2v5z7HCYGNPBg744DX1SpZRFLh60s58yNHANC0N73W9Mve38R33EO86MZSST6zozcBapUsIqnzdeiH8hNr+q4ufULfOcfWh/8PfYmx45Tb+UtpTK2SRaTT+Dr0mzttBurSp9Pmi0vv5/P1L7F65DVcMPUULkj1emgRkf34ek2fYBZ7ySPYkB79d3bs3MlnVv2IzVmf4rgLb/K6HBHJQP4OfaAmWEhO1PvQd86x7g/X0Y8Kss/7NYFQttcliUgG8n3o12UVkRvzvr3yS88t47Sqpaw7cjaDxpzsdTkikqF8H/oNoULCTd522txTVcOA5/8vHwb7MuqSH3tai4hkNt+HfmN2hIJ4lac1vPKHWxjBBzRMu4OscKGntYhIZvN96Mdzi4lQTUOsyZP5X1n5ClN2LmBDn6kcedJ5ntQgIv7h+9B34WKKrJbKvfXdPnd1XZTsp64lajkM+8qvun1+EfEf34d+IHmBVvVH3d+KYcVDP+V4t449J99MTmRQt88vIv7j+9DPSoZ+beXObp339fUbmLzll7zfazxHTbmyW+cWEf/y9xW5QHZBIvTrq3Z325z1jU18tOg6whYl+5K7IeD7370i0k18nzbhon4ANHZjp80lj97PlNgLbBt3NeGBo7ptXhER3x/p5ydDP1bTtU3Xmm99WFHxEctzbuOD0FCO+qJaLYhI9/L9kX6v4kR75Xht14V+860Pyyvq+G7WnxliH/K9hq/z+NpdXTaniEhLfB/6wXARTc4I1Hdd6N+xfAN1jU0ca+9xWfCv/CE2lZcbj9atD0Wk2/k+9AkEqLZeBOu7runa1oo6QsS4N3QnleTxXHw8VwaX6taHItLtFPpATaCQUBd22hxQlMt3sx6lf6CC3raX+7Lv5E33Kd36UES6nUIfqM0q7NJOm9PyN/HN4FJWNI3f91p2MMAN00Z22ZwiIi3x/dk7AA1ZRYTru+birKdKN3D57p+wKzSI99woZtedTd/8HG4YXclY3fpQRLqZQh9ozInQt3ZTp+/3gz21xJZex0Dbg2Rkk9kAAAeSSURBVPvqcq4YOpErOn0WEZH2S2l5x8wuNLN1ZhY3s5JDjLvPzHaa2VupzNdVYjkRCl01zrnO22dTnD/d/3POtReoPuk6soZO7LR9i4gcrlTX9N8Czgeeb2Pc/cD0FOfqEo+vLufV7Y58a+C025fz+OryTtnvfcte5PLKX7Gn+DgiX7ixU/YpIpKqlELfOfe2c67Nk82dc88DXXvJ62FovmhqR2MeALVVHzJ30dqUg/+Vd3cx9rUbyQ04en95PgS1iiYi6SHtzt4xsyvMrNTMSnft6torVpsvmjrB3gHg9MAavhp/PKWLpipqo7z60H/y2cA63JnzoM+nO6tcEZGUtRn6ZrbCzN5q4WNGVxTknLvHOVfinCvp169fV0yxT/PFUS/FxwDwk+zf8aErOOyLppxz/PKhxVwVe5DKo6aRe+KcTqtVRKQztLnu4Jyb2h2FeGFQJEx5RR3lfPzL5b9C8wkRoi46nXB2sEP7W/jyRi7c8h/EciMUXfRbMOvskkVEUpJ2yzvd6YZpIwmHgoyz95gdvYlvRr/DThdhnv2GZ39yMa+/u63d+3pnRzV1f/0howIfkHvBbyF5cxYRkXSS6imb55lZGTAJeNLMlidfH2Rmy/Yb9zDwMjDSzMrM7OupzNtZZk4YzO3nH8tfCi7ilfgY1haexuoZz1I25irOiT1N7oJp3LP46TZvml7f2MR9D9zP1wJPUjv+awRGfqGbvgMRkY6xzjw3vbOVlJS40tJST+auXbcMt+hK4rFGfpH/Hc7/8rcZPaiwxbHzHvs7c968hMKiCPlXvwTZed1crYhIgpmtcs61et2Ur5d3DiVvzFnk/9vfiff5DD+oncdrv/0G/71iPbGm+AHjVqzbztg1P+IIqyT/4vkKfBFJazqB/FAiR1L0rRXUP/UD5qz6H9Y8v5Fvrv8BJx53HA+8/D5bK+qYGXyRn4VeIXb6LTBofNv7FBHxkI7025KVTe65P4GLHmBM9g7u3H01ryx/mPKKOgbbLn6UNZ/X3DH8peBCrysVEWmTjvTba/QMQv3HsuPXFzA/+w7+O/ZFTgi8gwHfjV6Fe3oTM48f6nWVIiKHpCP9jujzaWbU3cpDscl8K2sJJwX+wQOxMzg78IrugiUiPYJCv4P6RIr4fuxyft54HgDfDi3RXbBEpMdQ6HdQ8wVdr7rR+17TXbBEpKfQmn4HzUze7ar8ySeZXXuT7oIlIj2KQv8wzJwwGCb8im97XYiISAdpeUdExEcU+iIiPqLQFxHxEYW+iIiPKPRFRHxEoS8i4iNp3U/fzHYBe4EPva7lMPVFtXulJ9ev2r3Tk+tvrv0o51yrNxhP69AHMLPSQ90QIJ2pdu/05PpVu3d6cv3trV3LOyIiPqLQFxHxkZ4Q+vd4XUAKVLt3enL9qt07Pbn+dtWe9mv6IiLSeXrCkb6IiHQShb6IiI+kbeib2XQz22Bmm8zsRq/r6Qgzu8/MdprZW17X0lFmdqSZPWdm681snZl9x+ua2svMcs1spZm9kaz9R17X1FFmFjSz1Wb2F69r6Sgz22xma81sjZmVel1PR5hZxMweNbN/mNnbZjbJ65raw8xGJn/ezR9VZnbNIbdJxzV9MwsC7wBnAGXAa8Bs59x6TwtrJzM7BagBHnDOjfW6no4ws4HAQOfc62ZWAKwCZvaEn72ZGZDvnKsxsxDwIvAd59wrHpfWbmZ2LVACFDrnzvG6no4ws81AiXOux13cZGYLgBecc/eaWTaQ55yr8LqujkjmZjlwknPu/dbGpeuR/kRgk3PuPedcFFgIzPC4pnZzzj0P7PG6jsPhnNvmnHs9+bgaeBvoEbcFcwk1yaeh5Ef6HdW0wsyGAGcD93pdi5+YWRFwCvB7AOdctKcFftIU4N1DBT6kb+gPBj7Y73kZPSR4MomZDQMmAK96W0n7JZdH1gA7gWeccz2mduDnwPeAuNeFHCYHPG1mq8zsCq+L6YDhwC5gfnJp7V4zy/e6qMNwMfBwW4PSNfTFY2bWC3gMuMY5V+V1Pe3lnGtyzo0HhgATzaxHLK+Z2TnATufcKq9rScHnnHPHA2cC304uc/YEWcDxwG+dcxNI9Pvqae8jZgNfBP7c1th0Df1y4Mj9ng9JvibdILke/hjwR+fcIq/rORzJP8+fA6Z7XUs7nQx8MbkuvhCYbGYPeltSxzjnypOfdwKLSSzT9gRlQNl+fxU+SuKXQE9yJvC6c25HWwPTNfRfA442s+HJ32AXA0s8rskXkm+G/h542zl3l9f1dISZ9TOzSPJxmMSJAP/wtqr2cc7Ndc4Ncc4NI/Hv/W/OuS97XFa7mVl+8o1/kksjXwB6xNlrzrntwAdmNjL50hQg7U9c+ITZtGNpBxJ/1qQd51zMzK4GlgNB4D7n3DqPy2o3M3sYOA3oa2ZlwK3Oud97W1W7nQx8BVibXBsH+L5zbpmHNbXXQGBB8iyGAPAn51yPO/Wxh+oPLE4cM5AFPOSc+6u3JXXIvwF/TB5kvgdc5nE97Zb8JXsGcGW7xqfjKZsiItI10nV5R0REuoBCX0TERxT6IiI+otAXEfERhb6IiI8o9EVEfEShLyLiI/8fgKBUbIQt13EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de3hU5bn38e89k8mBZJIgUEDAgoIIqJyiqBQUwYJWBbVWU23R1qptsVqLVov7bWu3u3bTTX19Sw9WWsGqqAiord14rlQ8gYCAgIgihjNoyPk087x/zIRySMhhkqxk1u9zXXMxh7XWcyfALyvPWute5pxDRET8IeB1ASIi0nYU+iIiPqLQFxHxEYW+iIiPKPRFRHwkxesCjqZr166ub9++XpchItJhrFixYq9zrlt9n7fr0O/bty/Lly/3ugwRkQ7DzD452uea3hER8RGFvoiIjyj0RUR8pF3P6YuIN6qrqykoKKCiosLrUqQe6enp9O7dm1Ao1KT1Egp9MzsGeBzoC2wBvuac+7yO5Y4DHgT6AA64wDm3JZGxRaT1FBQUEA6H6du3L2bmdTlyGOcc+/bto6CggH79+jVp3USnd+4AXnLODQBeir+uyzxgpnNuEHA6sDvBceu1eOU2Rt/7Mv3u+Duj732ZxSu3tdZQIkmroqKCLl26KPDbKTOjS5cuzfpNLNHQnwzMjT+fC0w5fAEzGwykOOdeAHDOlTjnyhIct06LV27jzoVr2FZYjgO2FZZz58I1Cn6RZlDgt2/N/ftJdE6/u3NuR/z5TqB7HcucCBSa2UKgH/AicIdzLpLg2EeYuWQj5dU1zAvdy+vRIbznTuDU6EfMXJLKlOG9Wno4EZEOp8E9fTN70czW1vGYfPByLtaYv67m/CnAGGA6cBpwPHDNUca73syWm9nyPXv2NOVrYXthOWCMDGzkztB8Hku9h/fc8fH3RSSZ3XfffZSVtcokAosXL+b9999v0jrXXnstf/zjH4/Yzvnnnw/EjptMnjyZAQMGcPzxxzNt2jQqKysBePXVV8nJyWHYsGEHHi+++GKLfC0Nhr5zboJz7uQ6Hk8Du8ysJ0D8z7rm6guAVc65j5xzNcBiYMRRxnvAOZfnnMvr1q3eK4nrdGxuBgBlpNf5voi0jvZwLK29hX5+fj7z588/5L358+eTn5+Pc45LL72UKVOmsGnTJjZt2kR5eTm33377gWXHjBnDqlWrDjwmTJjQIl9LonP6zwBT48+nAk/Xscw7QK6Z1Sb4uUDTvnuNdNvEgWSEgkRcgBWR/uRXzWBkysfcNnFgawwnIrTusbRZs2Zx8sknc/LJJ3PfffexZcsWTjrpJK666ioGDRrEV7/6VcrKyrj//vvZvn0748aNY9y4cUds56GHHmLKlCmcd9559O3bl9/+9rfMmjWL4cOHc8YZZ/DZZ58BsHnzZiZNmsTIkSMZM2YMGzZsYNmyZTzzzDPcdtttDBs2jM2bN/OnP/2J0047jaFDh3LZZZfV+cNm/PjxbNiwgR07YjPgpaWlvPjii0yZMoWXX36Z9PR0rr32WgCCwSC/+c1vmDdvHiUlJQl/344m0Tn9e4EnzOzbwCfA1wDMLA+40Tl3nXMuYmbTgZcsduRhBfCnBMetU+28/ceLepFi1WzNzuOKiVdpPl8kAT9/dh3vby+q9/OVWwupikQPea+8OsLtC97jsbe31rnO4GOz+elFQ4467ooVK/jLX/7CW2+9hXOOUaNGcfbZZ7Nx40bmzJnD6NGj+da3vsXvfvc7pk+fzqxZs3jllVfo2rVrndtbu3YtK1eupKKigv79+/OrX/2KlStX8sMf/pB58+Zxyy23cP311/OHP/yBAQMG8NZbb/G9732Pl19+mYsvvpgLL7yQr371qwDk5ubyne98B4C77rqLOXPmcNNNNx0yXjAY5LLLLuOJJ57g5ptv5tlnn+Wcc84hOzubdevWMXLkyEOWz87Opm/fvnz44YcALF26lGHDhh34/KmnnuKEE0446vesMRIKfefcPmB8He8vB6476PULwKmJjNVYU4b34t0Xu3BM+ce8fse5bTGkiK8dHvgNvd9Y//rXv7jkkkvIzMwE4NJLL2Xp0qX06dOH0aNHA3D11Vdz//33M3369Aa3N27cOMLhMOFwmJycHC666CIATjnlFN577z1KSkpYtmwZl19++YF1aufYD7d27VruuusuCgsLKSkpYeLEiXUul5+fz/Tp07n55puZP38+3/jGNxr99Y8ZM4a//e1vjV6+sZLyitya1GyySou9LkMkKTS0Rz763pfZVsfJEr1yM3j8hjNbvJ7DT1Ws69TFRYsW8fOf/xyABx98EIC0tLQDnwcCgQOvA4EANTU1RKNRcnNzWbVqVYM1XHPNNSxevJihQ4fy0EMP8eqrrxKJRA7svV988cXcfffdnHXWWezYsYPVq1ezbNmyA3P8gwcPZsGCBYdss6ioiJ07dzJw4EDeeuutxn47miwpe+9E0nMJu1KvyxDxhdpjaQfLCAUTPpY2ZswYFi9eTFlZGaWlpSxatIgxY8awdetW3njjDQAeffRRvvSlLwEQDocpLo7t7F1yySUHDoDm5eU1arzs7Gz69evHk08+CcSuel29evUR2wYoLi6mZ8+eVFdX88gjjwCx6ZzaMe+++24g9gPpiiuuYOrUqZx//vmkp8dOMhk/fjxlZWXMmzcPgEgkwo9+9COmTZtGRkbrnniSlKFPei5pVk1FWeseEBGR2JTqLy89hV65GRixPfxfXnpKwsfSRowYwTXXXMPpp5/OqFGjuO666+jcuTMDBw5k9uzZDBo0iM8//5zvfve7AFx//fVMmjSpzgO5jfXII48wZ84chg4dypAhQ3j66di5KVdeeSUzZ85k+PDhbN68mV/84heMGjWK0aNHc9JJJx11m/n5+axevZr8/PwD75kZixYtYsGCBQwYMIAuXboQCASYMWPGgWVq5/RrH4f/ZtBcFju9vn3Ky8tzzbmJyltP/ppR637B7u+s4gu9mtaXQkRg/fr1DBo0yOsyjrBlyxYuvPBC1q5d63UpLWrZsmXk5+ezaNEiRoyo94z2I9T192RmK5xz9f56k5Rz+imZXQAo3b8XFPoi0s6dddZZfPLJUW941WKScnonNesYAMqL9nlciYi0pL59+ybdXn5bS8rQT8+OhX5l0V6PKxERaV+SMvQzc2IX/1aXHtHaX0TE15Iz9HNjV+RFyxT6IiIHS8rQD2d3JuoMV67QFxE5WFKGfiAYpMgyCVQUel2KiDRTMBg85Dz1LVu2APD2228zduxYBg4cyPDhw7nuuusoKyvjoYceYtq0ad4W3QEk5SmbACWWRbCq/iZRItK+ZWRkHNESYdeuXVx++eXMnz+fM8+MtXhYsGDBIVfLytElbeiXB8OkVu33ugwRaUGzZ89m6tSpBwIfOND5UhoneUM/JZu0Gu3piyTsH3fAzjUNL7e/ANKyYs8rSyCnd/3L9jgFzr/3qJsrLy8/0Fq4X79+LFq0iLVr1zJ16tSjridHl7ShXx3KpnPljoYXFJGWkZYFu+IXTnU/OeHN1TW9I4lL2tCvSc0hy6nhmkjCGtgjP+Dj12BurEc9k34J/ca2eClDhgxhxYoVTJ48ueGFpU5JefYOQDQ9l7ArwUUTu5GDiDTStndh6rOxx7Z3W2WIadOmMXfu3EP6zS9cuJBdu3a1ynjJKGn39C0jlxSLUlKyn6zszl6XI5L8vnTLv5+3wl4+QPfu3Zk/fz7Tp09n9+7dBAIBxo4dy6RJk1plvGSUtKEfyIgFffHnuxX6Ih1QfTcIP/PMM1m6dOkR719zzTVcc801rVxVx5e00zspWbH2ymX71WlTRKRW0oZ+Wm175WKFvohIraQN/fTs2J5+lUJfpFna8131pPl/P0kb+rWdNmvUXlmkydLT09m3b5+Cv51yzrFv374DN1pviqQ9kBtWe2WRZuvduzcFBQXs2bPH61KkHunp6fTufZSrnuuRtKHfKTObahdUe2WRZgiFQvTrp/tLJ6Oknd6xQIAiyyJQqaZrIiK1kjb0AUoDWaQo9EVEDkjq0C8PhEmtVuiLiNRK6tCvCGWTXqObK4iI1Eoo9M3sGDN7wcw2xf88ot+BmY0zs1UHPSrMbEoi4zZWdSibTlGFvohIrUT39O8AXnLODQBeir8+hHPuFefcMOfcMOBcoAx4PsFxGyWSmkNY7ZVFRA5INPQnA3Pjz+cCDe3BfxX4h3OuLMFxGyWa0ZksV0Y0EmmL4URE2r1EQ7+7c6729lQ7ge4NLH8l8NjRFjCz681suZktT/TCEMvIJWCO4v2fJbQdEZFk0eDFWWb2ItCjjo9mHPzCOefMrN5rts2sJ3AKsORo4znnHgAeAMjLy0voGvBgp9ghhpLC3eQc0y2RTYmIJIUGQ985N6G+z8xsl5n1dM7tiIf67qNs6mvAIudcdTPqbJaQ2iuLiBwi0emdZ4DaW9NPBZ4+yrL5NDC109LSwrH2yhXqtCkiAiQe+vcC55nZJmBC/DVmlmdmD9YuZGZ9gT7APxMcr0ky1F5ZROQQCTVcc87tA8bX8f5y4LqDXm8BeiUyVnNk5cTm8dVeWUQkJqmvyM3qHG+vrE6bIiJAkod+ekYmFS6ElRd6XYqISLuQ1KEPUGRhApUKfRER8EHolwWySKkq8roMEZF2IelDvzwYJk3tlUVEAB+EfkUoh/SIOm2KiIAPQr8mlE1mRJ02RUTAB6EfScsh7LSnLyICPgh9l55LplVQXVXpdSkiIp5L+tAPxDttFhfu9bgSERHvJX3o/7u9cmK9+UVEkkHSh34oK9Zps6xITddERJI+9NPinTYri3X3LBGRpA/9TvHQr1Z7ZRGR5A/9zJxYp81ImTptiogkfehnd4711I+q06aISPKHfig1jVKXjqmnvohI8oc+QLFlEaxU0zUREV+EfmkwTEq12iuLiPgi9CuCYdLVXllExB+hXxXKJkPtlUVE/BH61ak5dIqqvbKIiC9CP5qWS7baK4uI+CP0Sc8l3aqpKC/1uhIREU/5IvStttPm52qvLCL+5ovQT8mMh/5+tVcWEX/zReinxtsrl6u9soj4nC9CPz071nRN7ZVFxO8SDn0zO8bMXjCzTfE/O9ez3H+b2TozW29m95uZJTp2Y3XKibdXLtGevoj4W0vs6d8BvOScGwC8FH99CDM7CxgNnAqcDJwGnN0CYzdKltori4gALRP6k4G58edzgSl1LOOAdCAVSANCwK4WGLtRsnK6EHWGU3tlEfG5lgj97s65HfHnO4Huhy/gnHsDeAXYEX8scc6tr2tjZna9mS03s+V79rTM2TbBlBRKrBMBtVcWEZ9LacxCZvYi0KOOj2Yc/MI558zM1bF+f2AQ0Dv+1gtmNsY5t/TwZZ1zDwAPAOTl5R2xreYqsSyCVWq6JiL+1qjQd85NqO8zM9tlZj2dczvMrCewu47FLgHedM6VxNf5B3AmcETot5bSQJiQ2iuLiM+1xPTOM8DU+POpwNN1LLMVONvMUswsROwgbp3TO62lIiVMukJfRHyuJUL/XuA8M9sETIi/xszyzOzB+DILgM3AGmA1sNo592wLjN1o1aFsMqJquiYi/tao6Z2jcc7tA8bX8f5y4Lr48whwQ6JjJaI6LZfMErVXFhF/88UVuQDRtByyXQkuGvW6FBERz/gm9C09l5BFKCvVvL6I+JdvQj+QGWu6Vlyo9soi4l++Cf3a9sqlCn0R8THfhH5te+UKtVcWER/zTehn1LZXLlF7ZRHxL9+Efm175Rq1VxYRH/NN6GfldgMgqqZrIuJj/gn9cC41LqD2yiLia74JfQsEKLYsAhUKfRHxL9+EPsTaK6eovbKI+JivQr8sqPbKIuJvvgr9ypRsMmoU+iLiX74K/apQNhkRtVcWEf/yVehH0nLIcmqvLCL+5avQj6bnku1KiUYiXpciIuIJX4W+ZeQSMEdxkS7QEhF/8lXoBzrFmq6p06aI+JWvQj9U2155/x6PKxER8YavQj8tHGu6VlGkTpsi4k++Cv2MnFh75Sq1VxYRn/JV6GfGQ7+mVO2VRcSffBX64dxY6EfLdPaOiPiTr0I/PSOTShcCddoUEZ/yVejH2itnElToi4hP+Sr0AUoDYVKq1HRNRPzJd6FfFgyTqk6bIuJTvgt9tVcWET/zXehXp2bTSe2VRcSnEgp9MzvGzF4ws03xPzvXs9yvzGxt/HFFImMmKpKWq/bKIuJbie7p3wG85JwbALwUf30IM/sKMAIYBowCpptZdoLjNptLzyVs5dRUV3lVgoiIZxIN/cnA3PjzucCUOpYZDLzmnKtxzpUC7wGTEhy32SwjF4DiQl2VKyL+k2jod3fO7Yg/3wl0r2OZ1cAkM+tkZl2BcUCf+jZoZteb2XIzW75nT8t3w0zJjLVXLtmv9soi4j8pDS1gZi8CPer4aMbBL5xzzszc4Qs55543s9OAZcAe4A2g3ltXOeceAB4AyMvLO2J7iQplxQ47lBWqvbKI+E+Doe+cm1DfZ2a2y8x6Oud2mFlPYHc927gHuCe+zqPAB82sN2FpWbH2ypXF6rQpIv6T6PTOM8DU+POpwNOHL2BmQTPrEn9+KnAq8HyC4zZbRk43AKpKFfoi4j8N7uk34F7gCTP7NvAJ8DUAM8sDbnTOXQeEgKVmBlAEXO2cq0lw3GbLzInt6UcU+iLiQwmFvnNuHzC+jveXA9fFn1cQO4OnXcjuHNvTV3tlEfEj312Rm5qWTplLwyr2e12KiEib813oAxRbFsFKtVcWEf/xZeiXBcKkVGlPX0T8J9EDuR1SWUqYtBo1XUs2i1duY+aSjWwvLOfY3AxumziQKcN7tdjybbWO6lJdTVmnqXwZ+lUp2eRWFHhdhhxFc/6D3blwDeXVsev+thWWc+fCNQB1rtfU5dtqHdWlupqyTnOYcy1+0WuLycvLc8uXL2/x7b59Xz59C9/kCz/7uMW3LUdqboBXV1eSQRUhaginRPj+2OMY9cVsaqoqiFRXEqmu/bOSh5Z+QFVFBRMD71DguhIwx3G2mxXBYXx5cDdwUVw0Ai4KLsqyTbupqq7hjMD77HKdMRw97HPWBE9iRJ9cwGHO4XDgHOBYW1BIdSTKcNvEbmI9nL5AIe8HBjD42FgPQTvs/9P6HUVURaIMs83sJie+zn7WB/pzUo/wEV/7hp3FVEUcQxu5fFuto7ratq6trhsLo2M51T7ib+Gv8fod59a5Tl3MbIVzLq/ez/0Y+m/+/kZO3bmQTj+v8wJiaUBTQnzxym3MWLiKY2p20dv2kk0ZnQNljDo2hZ6pFVjlfgJVRYSqikirKSY9UkJ6pJhsSsm0yjb+yo4UdUYs8i3+AHCkWhSAKhfEEeDg/0UOO+J5gCjpVg1AuQsRJXDQUgevy4HlMxqxfFuto7q8qSu/agZvRofw8b1fqWetIzUU+r6c3nEZuXSySiorykhL7+R1OR3K0X4FHX98Bjs2r2H/p+uo2f0BaYWbGVT8Me8GdpKWVn3ohnbHArXEOlFiWZQFsqgMZlGY0YdN+4PsJ5NcSrg85TUAHqi5gI3R48gf3Z9gSiqBUBqBlDSCoXQCoVR+8b+b2V3qONk+5jepvwfg+1U38UnmUB769igCgSCBQBALphAIBJg8exnbiqo5PbCBh1PvBeDrVT/hk/AI/vXj8Vggdo7DwWc6jL73ZbYVlnNmYB2Ppd4DwNTqO9ianVfvnlhd63yr+vZ612nq8m21jurypi6AY3Mz6ly+uXwZ+oGMWNO14sK9pPU4zuNqOpaZSzZSU13JL1P+QjUpHGt7GWDbSF9cRdgKqf3FtcYF2BHowQbXg1ejpwJwQ8rfAbipahr/jA5j1T2Xkh0McvjNFaYd9I//cmKh/0p0OFuz8/ifr9T9H+by4EncuXAN50ZXkl8V6wU4MuVjzrvgLLr2OPK3kB+cP5w7F65hcPSTQ5b/2qSrDwT+4W6bOJA7F67h1OhHh6xzxcSr6v1+NXWdthhDdSVPXc3hy9APZsXaK5cW7qWrQr9RopEI69/8X75f8icuSHubXCs98NnGaC/ejpxE7xOHkd5jIMccN5ie/QbTJy2dKw8K8BuIhf5ecgjndiEQDNY5VnP+8ddOL81cknpg2umKiVfVO+3U1OXbah3Vpbqask5z+HJO/71Xn+LUV7/FhvOf5KRRX27x7ScLF42yec0y9i77K8fvWsIX+IxSl8aS6Glsivbix6HHgdi8Y32/ttZOB30zupj33PFALMD7X3JXi5/uJiKa069Tenas6VpFiZqu1RWuI8OfUfDaw/Qq+Bv9o9s4zgVZlzmKT4dcxpYuY7nr7x/xTRY3ai+8uXsvU4b3UsiLtAJfhn6n7K4AVPs89A8+KNuNz/ly8T/ot+h1+gQ+opcz1qedwtsnXsfAcVcxvEvspmgjgWBaZpN/bVWAi7QPvgz9rNxY6Pu9vfLMJRvJrP6Mxan30N8KCBrsiHbm13yDb3znhwzpfUKd6ynERTouX/beCcdD35X7u+la/6I3+Ufaj+lnOwjGTx6+teZ7zK44n+71BL6IdGy+DP1gSgpFdMIq/Bn6lRVlvPn7G5ib+iv2uRx+Uv3tQz5v6fOCRaT98OX0DkCJZZHiw/bKWz9YRdXj3+KMyGaez7yI24ou50r7R6ueFywi7YdvQ78sECZUXeR1GW3GRaMsX/z/GLL6HiotlVWjf8+Xz/s6ZSu3MXNJuFXPCxaR9sO3oV/ho/bK+z/fy4dzvs1pJa+yLm0o3abOZVivfoAOyor4jW9DvyqUTbjsI6/LaHUb3n6BnOe+y1C3jzeO/z6nX3U3wRTf/rWL+J5v//fXpOaQWVLidRkt6uALrXrlhLg98zku2DeX3YFubL7wKc7Ma3x7VhFJTr4N/UhaZ7JdCS4arbfBVkdy8IVWPdnH/5TPZlTlBv6ZPo4R3/szx+Yc43WJItIO+Db0ycgl1WooKyumU1aO19UkbOaSjZRXR/ivlAe5OLiMIBGerTmDe9N/wOsKfBGJ6/i7uM0UyIjd+ai4cK/HlbSM7YXlfNF2Mjn4OllWQYZV82h0PNsLy70uTUTaEd+GfijeXrlsf3KE/vHZjgdCs4gcdl8eXWglIgfz7fRO6oHQ3+dxJYlz0Sh3M5v+to3HI+fwbPQsQBdaiciRfBv6te2Vq5Kg0+abD9/F6KplLO7+PWYXnacLrUSkXr4N/QPtlTt4p83VrzzJqI9+x/KcCUy+8R6mJMGZSCLSehJKCDO73MzWmVnUzOq9U4uZTTKzjWb2oZndkciYLSWrczcAoh049As+XEu/f97Mxyl9GXLDQ0lx6qmItK5EU2ItcCnE715dBzMLArOB84HBQL6ZDU5w3IRlhXOJOOuw7ZVLiwupeTSfKEbG1fPJyAw3vJKI+F5Coe+cW++c29jAYqcDHzrnPnLOVQHzgcmJjNsSAsEgxZZJoHK/16U0mYtG2fjHb9An8imfnvs7ju13ktcliUgH0RbzAb2ATw96XRB/r05mdr2ZLTez5Xv27GnVwkos3CHbK7/58F2MKHmNdwbcwiljPf/5KSIdSIMHcs3sRaBHHR/NcM493dIFOeceAB4AyMvLcy29/YOVBTtee+X3XlkQO3CbPZ5RX/8/XpcjIh1Mg6HvnJuQ4BjbgD4Hve4df89zFSlh0ms6TugXfLiWvv/8QezA7Y1zdeBWRJqsLVLjHWCAmfUzs1TgSuCZNhi3QdWhHDpFOkZP/dLiQqof/boO3IpIQhI9ZfMSMysAzgT+bmZL4u8fa2bPATjnaoBpwBJgPfCEc25dYmW3jJrUbDJd+2+vXHvg9rjIVj4997c6cCsizZbQxVnOuUXAojre3w5ccNDr54DnEhmrNUTTc8l2JUQjEQLBoNflHKG2P/7FxY/z49BrLO52A1PGXuJ1WSLSgfl6UtgyOhM0R0lx+zuDp7Y/fv+iN7kt5XGejZzBnbvOYfHKdnE4REQ6KF+HfrBTZwBK2mF75ZlLNpJbvZs/hGbxqevGgsgYvhl9hplLGrosQkSkfv4O/cz22155e2E5M0J/JYDji4HdzE2dyXvuePXHF5GE+Dr008KxTpsVRe0v9Cd22siFwbd4JnLWIe+rP76IJMK3XTYBMtppe+Xqqkpuc3/m02g3trge5FfNANQfX0QS5+vQ75QTC/2a0s89ruRQKxbM5Az3KY8P+CWLt52k/vgi0mJ8HfrZ8fbKkbL2E/r7dhUw+IPZvJc+kq9ddSNX6KpbEWlBvk6UjE5hqlwQ2lF75c3zbyfDVZJzySy1WRCRFufrVLFAgGLLItBOOm1+8O4/yfvsOVb0vJIvDhzmdTkikoR8HfoAJYEwoSrvQz8aicBz0/nMchiS/59elyMiScr3oV8eCJPaDtorL39mNifWfMDHw24nnHOM1+WISJLyfehXhsKk13jbabOocB8nrP41G1IGMfKiGz2tRUSSm+9DvzqUQ6eot6H//mM/obMrIuWiX7fLxm8ikjx8H/qRtByyPGyv/Mn6FYzc+STvdLmI/kO/5FkdIuIPvg/9aHou2ZQRqalp87FdNErx4lsps3ROzP/vNh9fRPzH96FvGbFOm8UedNpc+fw8Tq5cxYZBP6Bzt55tPr6I+I/vQ9+r9srlpcUc++Z/8lGgLyMvvbVNxxYR//J96Iey4u2Vi/a16bir5v+MHuyhfMIvSQmltunYIuJfvg/9dA/aK2//eAMjts5lRfhchpx1QcMriIi0EN+H/oH2yqVt115514IfESFA7yt+3WZjioiAz7tsAmTmdAUg0srtlWtvct6v6G3+mvovnu56HZN7n9CqY4qIHM73e/pZnWOhHy1rvT392puc7y4s5mcpc/k42p27dp+tm5yLSJvzfeinZ2RS7lKxVmyvPHPJRsqrI/whNIv+ge08ETmHr0ef003ORaTN+T70gXh75f2ttv3theUcy17OCqwD4Mehx3WTcxHxhEIfKA2ECVW1Xugfm5PGf4XmHPm+bnIuIm3M9wdyAcqDYVJrWq+98i3dVnBOwWqej4zkL5FJgG5yLiLeUOgDlaFsciq2t8q29+7cypcL/i/vBQbx87Q72b6/Sjc5FxHPKPSBmtQcMk5BYREAAAemSURBVMta56Dq1oe/xxBXRe6Vf+D1E3ULRBHxVkJz+mZ2uZmtM7OomeUdZbk/m9luM1ubyHitJZKWQ7gV2iu/+4+/MKJ0Ke+e8F2OU+CLSDuQ6IHctcClwGsNLPcQMCnBsVqNS8+lk1VSVVnRYtss3LuTL771UzYF+3Na/n+02HZFRBKRUOg759Y75xqcF3HOvQa0XZ+DJli8chuvfFINwJRf/73FLpja9PBNZLsSAlNmq6GaiLQbvj5ls/ZK2Z1V6QBUluzjzoVrEg7+1S8/wWn7n2dFn2s44ZQzWqJUEZEW0eCBXDN7EehRx0cznHNPt3RBZnY9cD3Acccd19KbP0TtlbKjUtYDcGfKo6yInsjMJanNPrOmeP9n9HjtDrYE+jD86v9syXJFRBLW4J6+c26Cc+7kOh4tHvjx8R5wzuU55/K6devWGkMcUHtF7N+isb3xCcGVfCXwJuWFu5q9zffn3UJX9xmVF9xPWnqnFqlTRKSl+Hp6599XxNqB904MFPB02n/w8bq3mry9da//nVH7nuadHlcyMO/cFqpSRKTlJHrK5iVmVgCcCfzdzJbE3z/WzJ47aLnHgDeAgWZWYGbfTmTclnLbxIFkhIKcah+RXzWD/KoZPB6dQKrV0P2Ji3h3ycON3lZ5aTE5L95KgfVk6DdntmLVIiLNl9DFWc65RcCiOt7fDlxw0Ov8RMZpLbXz9jOXpMaaouVmEJ54FYHu1RT8+XJGvDGNN7atYdTUXxIIBo+6rdXzbuMMt5N1X36M3pnhtihfRKTJzDnndQ31ysvLc8uXL/dk7IryUtb84VpO27+EdzPHMvDGh8kM59a57IblL3His5fxTtfJjLppbhtXKiLyb2a2wjlX78Wyvp7TP5r0jEzybp7PmwNuZWjJUnb95my2bznykoTKijIynruZ3daFwd/8jQeViog0nkL/KCwQ4Iyrfsq6cXPoGt1DxkMTWLfsuUOWefevM/hi9FN2nX0v4ZxjPKpURKRx1HCtEU495zI+7dUf91g+Jy65mqfW3MKsz79Ezv4NPJ06l39mnMvZ4y73ukwRkQZpT7+R+gwYSu4PlrIyNJzLdvwP3y35Lb8K/ZFCsri99Ou6362IdAgK/SbIzu3CrcEf8/uai7g65SVOCWzh4ZoJTIm+pPvdikiHoNBvom37q/lVTT6zqi8D4NbQU7rfrYh0GAr9Jqq9ivcdd1Kd74uItGcK/Saq6yrekSkfc9vEgV6XJiLSIJ2900R1XcWr+92KSEeh0G+GKcN7KeRFpEPS9I6IiI8o9EVEfEShLyLiIwp9EREfUeiLiPhIu+6nb2Z7gFJgr9e1NFNXVLtXOnL9qt07Hbn+2tq/6Jyr9wbj7Tr0Acxs+dFuCNCeqXbvdOT6Vbt3OnL9ja1d0zsiIj6i0BcR8ZGOEPoPeF1AAlS7dzpy/ardOx25/kbV3u7n9EVEpOV0hD19ERFpIQp9EREfabehb2aTzGyjmX1oZnd4XU9TmNmfzWy3ma31upamMrM+ZvaKmb1vZuvM7Gava2osM0s3s7fNbHW89p97XVNTmVnQzFaa2d+8rqWpzGyLma0xs1VmttzreprCzHLNbIGZbTCz9WZ2ptc1NYaZDYx/v2sfRWZ2y1HXaY9z+mYWBD4AzgMKgHeAfOfc+54W1khmNhYoAeY55072up6mMLOeQE/n3LtmFgZWAFM6wvfezAzIdM6VmFkI+Bdws3PuTY9LazQzuxXIA7Kdcxd6XU9TmNkWIM851+EubjKzucBS59yDZpYKdHLOFXpdV1PEc3MbMMo590l9y7XXPf3TgQ+dcx8556qA+cBkj2tqNOfca8BnXtfRHM65Hc65d+PPi4H1QIe4eYCLKYm/DMUf7W+vph5m1hv4CvCg17X4iZnlAGOBOQDOuaqOFvhx44HNRwt8aL+h3wv49KDXBXSQ4EkmZtYXGA685W0ljRefHlkF7AZecM51mNqB+4DbgajXhTSTA543sxVmdr3XxTRBP2AP8Jf41NqDZpbpdVHNcCXwWEMLtdfQF4+ZWRbwFHCLc67I63oayzkXcc4NA3oDp5tZh5heM7MLgd3OuRVe15KALznnRgDnA9+PT3N2BCnACOD3zrnhxPp9dbTjiKnAxcCTDS3bXkN/G9DnoNe94+9JG4jPhz8FPOKcW+h1Pc0R//X8FWCS17U00mjg4vi8+HzgXDP7q7clNY1zblv8z93AImLTtB1BAVBw0G+FC4j9EOhIzgfedc7tamjB9hr67wADzKxf/CfYlcAzHtfkC/GDoXOA9c65WV7X0xRm1s3McuPPM4idCLDB26oaxzl3p3Out3OuL7F/7y875672uKxGM7PM+IF/4lMjXwY6xNlrzrmdwKdmNjD+1nig3Z+4cJh8GjG1A+30xujOuRozmwYsAYLAn51z6zwuq9HM7DHgHKCrmRUAP3XOzfG2qkYbDXwDWBOfGwf4iXPuOQ9raqyewNz4WQwB4AnnXIc79bGD6g4siu0zkAI86pz7X29LapKbgEfiO5kfAdd6XE+jxX/Ingfc0Kjl2+MpmyIi0jra6/SOiIi0AoW+iIiPKPRFRHxEoS8i4iMKfRERH1Hoi4j4iEJfRMRH/j+NLcy29N/i6QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de3iU9Z338fc3k5AECERIOCUgQSAcRMGNrC62oliDh4L12e3iPmuvbWvtQfvY2qXC012t3YMH2l7W1T2w2tN2Vx7bIrCrEhW0tp6hIBAggKCSiSEBDQdJSDL5Pn8k0EBOEzLJPTP5vK6LS+ae3z3zEZ0Pv/zmPpi7IyIiiS8l6AAiIhIbKnQRkSShQhcRSRIqdBGRJKFCFxFJEqlBvXFOTo6PHz8+qLcXEUlIGzduPOjuue09F1ihjx8/ng0bNgT19iIiCcnM3uvoOS25iIgkCRW6iEiSUKGLiCSJwNbQ29PQ0EB5eTl1dXVBR0kYGRkZ5Ofnk5aWFnQUEQlYXBV6eXk5WVlZjB8/HjMLOk7cc3cOHTpEeXk5BQUFQccRkYDFVaHX1dWpzLvBzBg+fDjV1dVBR0kaqzaFWVZSRkVNLWOyM1lcXMgNs/JiNr6v9lGu5MjVXXFV6IDKvJv059Wxs/lQL125ldqGCADhmlqWrtwK0O5+3R3fV/soV3LkOhsW1OVzi4qK/Mzj0Hfs2MHUqVMDyZPI+sOfW0/LGSAzLYV7F0znqmmjONEYoa6h6Q//bIjwtf/8PYc+rm/zWtmZadx1zRSa3Glqcpocmtx56IVdHK5tbDN+SEYqX778PNydkx8vB9zhsd/t5Whd232yMlL5wpyCU2Nb+8kr+9rfJz2Vz88Z32b7T155l6Mnoh/fV/soV/v75GVn8sqSK9vdpz1mttHdi9p9ToWe+JL9z639cg7xj585nyumjODAkRMcOFLHgSN1VB1t/v2TG/ZT19AUYOq+0d4PaJ19pDv6ga4v9lGuDsYD++6/ruMXbPP6HRd63C25dEdvrEeFQiFmzJhBQ0MDqampfO5zn+Ob3/wmKSkdH+F533338fjjjxMKhXj44YcpLi7ucOzatWu54447iEQi3HLLLSxZsqTDsY888ggPPfQQ77zzDtXV1eTk5PTo3y1RPViy87QyB6htiPDNJ99ud/zQzLROy/yeT08jIy1EemoKGWkhMtJSSE8NcceKTRw81naGPnJIOqtum0PIDDMjlGKkGFzzo9/yweG2R2SNGZrBS4uvAJo/3Cc/32bGJx5cT0VN2306m6XNuX894ZraqPfp7vi+2ke52t9nTHZmu+PPRsIeh35y1hauqcX5w3rUqk3hHr1uZmYmmzdvprS0lOeff55nn32We++9t8Px27dvZ8WKFZSWlrJ27Vq+9rWvEYlE2h0biUS47bbbePbZZ9m+fTtPPPEE27dv7/C158yZwwsvvMC5557bo3+neLRqU5g596+nYMnTzLl//Wn/3dyddw9+zC9ef4+v/mJjuwV40t9eP41H/mIWv/zKpby8+Ap2/t183r7navI6+JDkZWfy+TkF3DR7HDdelM+1M0Zz5ZSRzJmYw99cN43MtNBp4zPTQiy9Ziqjh2YyYkgGuVnpDBs0gOyBA7hr/pR2x397/hQGpKYwIDWFtFAKqS2/QinGt4vb32dxcWGH/46Liwu7tU93x/fVPsrV/X26K25n6Pf+dynbK450+Pym92uoj5w+C6ttiPDtX23hiTffb3efaWOGcM+np0edYcSIESxfvpyLL76Y7373u+1+Abl69WoWLVpEeno6BQUFTJw4kTfffJNLL720zdg333yTiRMnMmHCBAAWLVrE6tWrmTZtWrvvP2vWrKizJpL2vhxasnILG9//iPqGJn635+CpmcyYoRkMHBDieH3bvyTzsjP54mXtH665uLiw3WWazj48J3+6i/anvu6O76t9lCs5cp2NuC30rpxZ5l1tP1sTJkwgEolQVVXFyJEj2zwfDoe55JJLTj3Oz88nHG7/p4RwOMzYsWNPG/vGG2/ENG8iWFZS1mYJpa6hif947T2GZKTyJ+fl8JW553HZxBzGDx/I6s0VvV7Orffrzgesu+P7ah/lSo5c3RW3hd7VTLqzNaz/9+W2s2OJHxXt/HeD5rXmTXdfTSjl9J+E+qqcRRJd3BZ6V87mR+qzsXfvXkKhECNGjGj3+by8PPbv33/qcXl5OXl57ZdId8Ymo6N1Dfzohd1tDss7aUx2ZpsyP0nlLNK1hP1S9IZZedx34wzysjMxmmfm9904I6Yf+urqar7yla9w++23d3gCz4IFC1ixYgUnTpxg37597N69m9mzZ7c79uKLL2b37t3s27eP+vp6VqxYwYIFC2KWN165O2vermDeD37D46/s49IJw8hIO/1/vd74y1ikv0nYGTr0zqyttraWmTNnnjps8eabb+bOO+/scPz06dP57Gc/y7Rp00hNTeXRRx8lFAq1OzY1NZVHHnmE4uJiIpEIX/jCF5g+veOlpYcffpgHH3yQyspKLrjgAq699loee+yxHv879qU9VUe5e3Upr75ziPPzhrD8c0XMHJvd66dAi/RHOrEoCcTjn9vx+kYeXreHx3+3t3n2PX8KfzF7XIdLKiISnaQ9sUjiw+mz7QyunjaSktIDVByu40//KJ8l10whZ3B60DFFkp4KPUolJSXcddddp20rKCjgqaeeajP20KFDzJs3r832devWMXz48DbbP/OZz7Bv377Ttj3wwAOdnnEaL9oeU17HT159j9FD0vnVVy6laPywgBOK9B9dFrqZ/Ri4Hqhy9/M7GXcx8BqwyN1/dbaB3D0uryBYXFwcdcEOHz6czZs3R/3a7f2lEK2glsxOau+YcgBLMZW5SB+L5iiXnwLzOxtgZiHgAeC5noTJyMjg0KFDgZdUojh5g4uMjIzAMnR0TPkHnZyuLyK9o8sZuru/bGbjuxj2deDXwMU9CZOfn095eblu2NANJ29BF5Rhgwa0e8nZWF5wSESi0+M1dDPLAz4DXEEXhW5mtwK3AowbN67N82lpabqVWgLZU3WMYycaME6/hreOKRcJRixOLHoIuMvdu7yIirsvd/cidy/Kzc2NwVtLUD76uJ4v/uwtsjLS+Nvrp/bqCV4iEp1YHOVSBKxo+SIzB7jWzBrdfVUMXlviUH1jE1/+xUY+OFzHilsv4aJx5/CFyyYEHUuk3+txobv7qTUSM/sp8D8q8+Tl7nznqa28ue9DfrRoJheNOyfoSCLSIprDFp8A5gI5ZlYO3AOkAbj7v/ZqOok7//byXn65sZw75k1i4Uwtq4jEk2iOcrkp2hdz97/qURqJayWllTywdifXXzCab1w1Keg4InKGhL3aovStbeHDfGPFZi7Iz+b7f3ZhXJ78JdLfqdClSweO1HHLzzZwzsA0/v1zf0RGWvtXkxSRYOlaLtKp2voIX/r5Bo7UNfDrr/4JI7KCOytVRDqnQpcONTU53/rlZraGD/PvNxcxdfSQoCOJSCdU6HKa1pfCHZyeytETjXzn2qlcNa3tDbJFJL6o0OWUMy+Fe/REI6EUI2fwgICTiUg09KWonNLepXAjTc73n9sVUCIR6Q4VupzS0aVwO9ouIvFFhS6ndHTJW10KVyQxqNDllL++ejJnni+kS+GKJA4VupwyJDMNd8gemKZL4YokIB3lIkDzVRQfemE344YNZN23LictpL/rRRKNPrUCwPqdVWwNH+b2KyeqzEUSlD65ctrs/DNaXhFJWCp00excJEno09vPaXYukjxU6P2cZuciyUOf4H5Ms3OR5KJC78c0OxdJLvoU91OanYskHxV6P/VimWbnIsmmy0+ymf3YzKrMbFsHz/9vM9tiZlvN7FUzuzD2MSWWNDsXSU7RTM1+Cszv5Pl9wOXuPgP4O2B5DHJJL3qxrIot5ZqdiySbLq/l4u4vm9n4Tp5/tdXD14H8nseS3qLZuUjyivX07IvAsx09aWa3mtkGM9tQXV0d47eWaGh2LpK8YvaJNrMraC70uzoa4+7L3b3I3Ytyc3Nj9dYSJc3ORZJbTArdzC4AHgMWuvuhWLymxJ5m5yLJrcefajMbB6wEbnZ33U04Tml2LpL8uvxS1MyeAOYCOWZWDtwDpAG4+78CdwPDgX+25vuXNbp7UW8FlrNzcnb+4J9eoNm5SJKK5iiXm7p4/hbglpglkpjT7Fykf9BULcmt2hSm6O9fYEv5YY7UNfD0lg+CjiQivUT3FE1iqzaFWbpyC7UNTQDUHG9g6cqtALrxs0gS0gw9iS0rKTtV5ifVNkRYVlIWUCIR6U0q9CRWUVPbre0ikthU6Els9NCMdrePyc7s4yQi0hdU6Ems+PxRbbZlpoVYXFwYQBoR6W0q9CRWdeQEg9NDjMnOwIC87Ezuu3GGvhAVSVI6yiVJHa1r4IUdB1h08VjuXXh+0HFEpA9ohp6knis9wInGJhbM1GxcpL9QoSep1W9XkH9OJheNyw46ioj0ERV6Ejp47ASv7DnIggvH0HJ9HRHpB1ToSeiZrR8QaXIWarlFpF9RoSeh1ZsrmDIqi8JRWUFHEZE+pEJPMvs/PM7G9z5iwcwxQUcRkT6mQk8y/72lAoBPX6BCF+lvVOhJZs3mCv7o3HMYO2xg0FFEpI+p0JPIzsoj7Kw8ykItt4j0Syr0JLJmcwWhFOPaGaODjiIiAVChJwl3Z83bFcyZmEPO4PSg44hIAFToSeL379dQ/lEtCy/UcotIf6VCTxJrNodJT03h6ukjg44iIgHpstDN7MdmVmVm2zp43szsYTPbY2ZbzOyi2MeUzjRGmnh66wfMmzqCrIy0oOOISECimaH/FJjfyfPXAJNaft0K/EvPY0l3vPrOIQ4eq2fBhTrVX6Q/67LQ3f1l4MNOhiwEfu7NXgeyzUyHWfSh1ZsryMpIZW5hbtBRRCRAsVhDzwP2t3pc3rKtDTO71cw2mNmG6urqGLy11DVEKCmt5JrzR5GRFgo6jogEqE+/FHX35e5e5O5FubmaTcbC+p1VHDvRqOUWEYlJoYeBsa0e57dskz6wenOYnMHpXHre8KCjiEjAYlHoa4DPtRztcglw2N0/iMHrShcO1zbwYlk1n75wNKEU3chCpL/r8ibRZvYEMBfIMbNy4B4gDcDd/xV4BrgW2AMcBz7fW2HldCWlldQ3NulGFiICRFHo7n5TF887cFvMEknU1myu4NzhA7kwf2jQUUQkDuhM0QRVdbSOV9/RfUNF5A9U6Anq6S0f0OToUrkicooKPUGt3lzBtNFDmDhC9w0VkWYq9AT03qGP2by/RvcNFZHTqNAT0JrNLfcN1aVyRaQVFXqCcXdWv13B7PHDyMvODDqOiMQRFXoCWbUpzB//4zr2VB2j7MBRVm3SCbki8gddHocu8WHVpjBLV26ltiECNJ8lunTlVgBumKUTi0REM/SEsayk7FSZn1TbEGFZSVlAiUQk3qjQE0RFTW23totI/6NCTxBjOvgCtKPtItL/qNATxOLiQs48wz8zLcTi4sJgAolI3FGhJ4jLJuXgDlkZqRiQl53JfTfO0BeiInKKjnJJEC+VNd+y74kvXcL5ebq6ooi0pRl6gli/8wAjh6QzfcyQoKOISJxSoSeA+sYmXt51kCunjNClckWkQyr0BPDWux9y7EQjV04ZGXQUEYljKvQEsG5HFQNSU5gzUTeCFpGOqdDjnLuzbucB/uS84QwcoO+wRaRjKvQ4t/fgx7x36DjzpowIOoqIxLmoCt3M5ptZmZntMbMl7Tw/zsxeNLNNZrbFzK6NfdT+ad2OAwBcoUIXkS50WehmFgIeBa4BpgE3mdm0M4b9DfCku88CFgH/HOug/dW6HVVMGZVF/jkDg44iInEumhn6bGCPu+9193pgBbDwjDEOnDxAeihQEbuI/dfh4w1seO8j5k3V7FxEuhZNoecB+1s9Lm/Z1tp3gb80s3LgGeDr7b2Qmd1qZhvMbEN1dfVZxO1ffrO7mkiT63BFEYlKrL4UvQn4qbvnA9cC/2FmbV7b3Ze7e5G7F+Xm5sborZPX+h0HGDZoADPHZgcdRUQSQDSFHgbGtnqc37KttS8CTwK4+2tABpATi4D9VWOkiZd2VTO3MJdQis4OFZGuRVPobwGTzKzAzAbQ/KXnmjPGvA/MAzCzqTQXutZUemDT/hpqjjcwT8stIhKlLgvd3RuB24ESYAfNR7OUmtn3zGxBy7BvAV8ys7eBJ4C/cnfvrdD9wbodVaSmGJ+YrB90RCQ6UZ166O7P0PxlZ+ttd7f6/XZgTmyj9W/rdx5gdsEwhmSkBR1FRBKEzhSNQ/s/PM6uA8e4UicTiUg3qNDj0PqdVQDMm6r1cxGJngo9Dq3bWcWEnEEU5AwKOoqIJBAVepz5+EQjr79zSMstItJtKvQ487s9B6mPNHGlTvcXkW5SoceZdTsOkJWRysXjhwUdRUQSjAo9jjQ1Oet3VnP55FzSQvpPIyLdo9aII1vDhzl47ISurigiZ0WFHkfW7awixeDyySp0Eek+FXocWb/zABeNO4dhgwYEHUVEEpAKPU5UHq5jW/iIjm4RkbOmQo8TL5a1nB2qqyuKyFlSoceJdTuqyMvOZPLIwUFHEZEEpUKPA3UNEV7Zc5B5U0dgpptZiMjZUaHHgdf2HqK2IaLT/UWkR1TocWD9jioy00JcMmF40FFEJIGp0APm7qzfWcVlk3LISAsFHUdEEpgKPWBlB44SrqllnpZbRKSHVOgBW7ej+XBFrZ+LSE9FdU9Rib1Vm8IsKykjXFNLWsh49Z1D3DArL+hYIpLAVOgBWLUpzNKVW6ltiADQEHGWrtwKoFIXkbMW1ZKLmc03szIz22NmSzoY81kz225mpWb2X7GNmVyWlZSdKvOTahsiLCspCyiRiCSDLmfoZhYCHgU+BZQDb5nZGnff3mrMJGApMMfdPzIzLQh3oqKmtlvbRUSiEc0MfTawx933uns9sAJYeMaYLwGPuvtHAO5eFduYyWVMdma3touIRCOaQs8D9rd6XN6yrbXJwGQze8XMXjez+e29kJndamYbzGxDdXX12SVOAouLCxkQOv0U/8y0EIuLCwNKJCLJIFaHLaYCk4C5wE3Av5tZ9pmD3H25uxe5e1Fubm6M3jrx3DArjwvyh2IGBuRlZ3LfjTP0haiI9Eg0R7mEgbGtHue3bGutHHjD3RuAfWa2i+aCfysmKZNMQ6SJXQeOceOsfH7w2QuDjiMiSSKaGfpbwCQzKzCzAcAiYM0ZY1bRPDvHzHJoXoLZG8OcSeX1vYc4UtfI/PNHBR1FRJJIl4Xu7o3A7UAJsAN40t1Lzex7ZragZVgJcMjMtgMvAovd/VBvhU50a7dVMnBAiE9Mygk6iogkkahOLHL3Z4Bnzth2d6vfO3Bnyy/pRKTJKSk9wBWFI3QxLhGJKV3LpY9tev8jDh47QbGWW0QkxlTofWzttkoGhFK4orD/HuUjIr1Dhd6H3J21pZVcNimHrIy0oOOISJJRofeh0oojlH9Uy/zpWm4RkdhTofehktJKUgyumjYy6CgikoRU6H1o7bZK/rhgOMMGDQg6iogkIRV6H9lTdYzdVcd0MpGI9BoVeh8pKa0E4OrpWm4Rkd6hQu8jJaWVzBybzeihukSuiPQOFXofCNfUsqX8sJZbRKRXqdD7QMm25uWWYh2uKCK9SIXeB9aWVjJlVBYFOYOCjiIiSUyF3suqj57grXc/1OxcRHqdCr2XvbDjAO5o/VxEep0KvZet3VbJucMHMmVUVtBRRCTJqdB70eHaBl595yDzp4/CzLreQUSkB1TovejFnVU0RFzXPheRPqFC70Vrt1Uyckg6M/Ozg44iIv2ACr2X1NZHeGlXFcXTR5GSouUWEel9KvRe8ptd1dQ1NOna5yLSZ6IqdDObb2ZlZrbHzJZ0Mu5/mZmbWVHsIiamktJKsgemMbtgWNBRRKSf6LLQzSwEPApcA0wDbjKzae2MywLuAN6IdchEU9/YxAs7DvCpqSNJDemHIBHpG9G0zWxgj7vvdfd6YAWwsJ1xfwc8ANTFMF9Cem3vIY7WNepkIhHpU9EUeh6wv9Xj8pZtp5jZRcBYd386htkS1tptlQwaEGLOxJygo4hIP9Lj9QAzSwF+CHwrirG3mtkGM9tQXV3d07eOS5Em5/ntlVwxZQQZaaGg44hIPxJNoYeBsa0e57dsOykLOB94yczeBS4B1rT3xai7L3f3Incvys3NPfvUcWzjex9x8Fi9lltEpM9FU+hvAZPMrMDMBgCLgDUnn3T3w+6e4+7j3X088DqwwN039EriOLd2WyUDUlOYWzgi6Cgi0s90Weju3gjcDpQAO4An3b3UzL5nZgt6O2AicXdKSiv55KQcBqenBh1HRPqZqFrH3Z8Bnjlj290djJ3b81iJaVv4COGaWr5x1aSgo4hIP6RpZAys2hRmWUkZ4ZpaoPk4dBGRvqZC76FVm8IsXbmV2obIqW1///QOBqWncsOsvE72FBGJLZ3G2EPLSspOK3OA2oYIy0rKAkokIv1VQs3QTy5tVNTUMiY7k8XFhYHPgitallmi3S4i0lsSZoZ+cmkjXFOLA+GaWpau3MqqTeEu9+1NY7Izu7VdRKS3JEyhx+vSxuLiQkJnXO48My3E4uLCYAKJSL+VMIUer0sbk0dmEXHISk/FgLzsTO67cUbgS0Ei0v8kzBr6mOzMU4cFnrk9SD98fhdDMlL57V1XMjQzLdAsItK/JcwMfXFxIZlnXOwqNcUCXdp4e38NL+w4wJc+MUFlLiKBS5gZ+skljJNHuWSkhahvjAR6R6AfPL+Lcwam8fnLCgLLICJyUsLM0KG51F9ZciX77r+O5+/8JKFQSmBfir717oe8vKuar849T9dtEZG4kFCF3lr+OQO55bICntoU5u39NX363u7O90vKyM1K5+ZLxvfpe4uIdCRhCx3gq3PPI2fwAP7h6R24e5+976vvHOKNfR9y29zzyBygm1iISHxI6ELPykjjm5+azJvvfkhJaWWfvKe78/3nyhg9NINFs8f1yXuKiEQjoQsd4M+LxjJpxGDue3Znn1zl8KWyaja9X8PXr5ykW8yJSFxJ+EJPDaXwneum8t6h4/z8tXd79b3cnR88X8a4YQP5s6L8Xn0vEZHuSvhCB5hbOIJPTMrhn9bvoeZ4fa+9T0npAbaFj/B/5k0iLZQUf3QikkSSppW+c91UjtY18KN1u3vl9SNNzg+fL2NC7iBumDmmV95DRKQnkqbQp4wawp9fPJb/eO099lYfi/nr/8+WCnYdOMY3rppMqmbnIhKHkqqZvvmpyaSnpnD/sztj+rqNkSZ+9MJuCkdmcf2M0TF9bRGRWEmqQh+RlcFX557Hc9sP8PreQzF73ac2hdl78GPuvHoyKSnW9Q4iIgGIqtDNbL6ZlZnZHjNb0s7zd5rZdjPbYmbrzOzc2EeNzi2fmMCYoRn8/dPbaWrq+clG9Y1N/GjdbmbkDeXqaSNjkFBEpHd0WehmFgIeBa4BpgE3mdm0M4ZtAorc/QLgV8CDsQ4arYy0EIvnF7ItfIRVm3t+N6NfbtxP+Ue13Hn1ZMw0OxeR+BXNDH02sMfd97p7PbACWNh6gLu/6O7HWx6+DgR6kPbCC/O4IH8oD64to7Y+0vUOHahriPBP6/Zw0bhs5k7OjWFCEZHYi6bQ84D9rR6Xt2zryBeBZ9t7wsxuNbMNZrahuro6+pTdlJJi/M1106g8Usdjv9171q/zX2+8T+WROv766kLNzkUk7sX0uq9m9pdAEXB5e8+7+3JgOUBRUVGvXk1rdsEw5k8fxb/85h3+/OKxjBiSEfW+qzaFeWDtTj44XMeA1BSqjp7oxaQiIrERTaGHgbGtHue3bDuNmV0FfAe43N3jogGXXDOF57ZXMvf7L1FbH2FMdiaLiws7vd/nqk1hlq7cQm1D83Vh6hubWLpyK4DuEyoicS2aQn8LmGRmBTQX+SLgL1oPMLNZwL8B8929KuYpz9Lm/TWYGcdb1tHDNbUsWbmFisO1nD9mKFVHT1B1tI6qIyeobvn979+vIXLG0TG1DRGWlZSp0EUkrnVZ6O7eaGa3AyVACPixu5ea2feADe6+BlgGDAZ+2bLW/L67L+jF3FFZVlLWppzrGpp4cO3pdzkanJ7KiKx0crPS24w/qaKdG1SLiMSTqNbQ3f0Z4Jkztt3d6vdXxThXTHRWwk9++VJGZKUzYkg6Awf84Y9hzv3rCbez35jszF7JKCISK0l1puiZOirhvOxMZhcMY3zOoNPKHGBxcSGZZ1znPDMtxOLiwl7LKSISC0ld6GdTzjfMyuO+G2eQl52J0Vz+9904Q+vnIhL3kvp29SdLeFlJGRU1tVEd5XJyPxW4iCSapC50UDmLSP+R1EsuIiL9iQpdRCRJqNBFRJKECl1EJEmo0EVEkoS59+pFDzt+Y7Nq4GPgYCABYiOHxM2v7MFJ5PzKHozW2c9193Zv0BBYoQOY2QZ3LwosQA8lcn5lD04i51f2YESbXUsuIiJJQoUuIpIkgi705QG/f08lcn5lD04i51f2YESVPdA1dBERiZ2gZ+giIhIjKnQRkSQRWKGb2XwzKzOzPWa2JKgc3WVmPzazKjPbFnSW7jKzsWb2opltN7NSM7sj6EzdYWYZZvammb3dkv/eoDN1l5mFzGyTmf1P0Fm6w8zeNbOtZrbZzDYEnae7zCzbzH5lZjvNbIeZXRp0pmiYWWHLn/nJX0fM7Bsdjg9iDd3MQsAu4FNAOc03or7J3bf3eZhuMrNPAseAn7v7+UHn6Q4zGw2Mdvffm1kWsBG4IRH+3AGs+Ya1g9z9mJmlAb8D7nD31wOOFjUzuxMoAoa4+/VB54mWmb0LFLl7Qp6YY2Y/A37r7o+Z2QBgoLvXBJ2rO1p6Mwz8sbu/196YoGbos4E97r7X3euBFcDCgLJ0i7u/DHwYdI6z4e4fuPvvW35/FNgBJMzF4r3ZsZaHaS2/EuZbfTPLB64DHgs6S39iZkOBTwKPA7h7faKVeYt5wDsdlTkEV+h5wP5Wj8tJoGJJBmY2HpgFvBFsku5pWbLYDFQBz7t7IuV/CPg20BR0kLPgwHNmttHMbg06TDcVANXAT1qWu4Ribd8AAAG8SURBVB4zs0FBhzoLi4AnOhugL0X7ITMbDPwa+Ia7Hwk6T3e4e8TdZwL5wGwzS4hlLzO7Hqhy941BZzlLl7n7RcA1wG0tS4+JIhW4CPgXd59F8zWkEuZ7O4CWZaIFwC87GxdUoYeBsa0e57dsk17Wsvb8a+A/3X1l0HnOVsuPzC8C84POEqU5wIKWtegVwJVm9otgI0XP3cMt/6wCnqJ52TRRlAPlrX6a+xXNBZ9IrgF+7+4HOhsUVKG/BUwys4KWv3kWAWsCytJvtHyp+Diww91/GHSe7jKzXDPLbvl9Js1fqu8MNlV03H2pu+e7+3ia/39f7+5/GXCsqJjZoJYv0WlZqrgaSJijvNy9EthvZoUtm+YBCXEgQCs30cVyCwR0k2h3bzSz24ESIAT82N1Lg8jSXWb2BDAXyDGzcuAed3882FRRmwPcDGxtWYcG+L/u/kyAmbpjNPCzlm/7U4An3T2hDv9LUCOBp5rnA6QC/+Xua4ON1G1fB/6zZQK5F/h8wHmi1vKX6KeAL3c5Vqf+i4gkB30pKiKSJFToIiJJQoUuIpIkVOgiIklChS4ikiRU6CIiSUKFLiKSJP4/sGhmQrKuHnwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(list(meta.keys()), list(meta.values()), label=\"meta-VQE\", marker=\"o\")\n",
    "plt.plot(list(fci.keys()), list(fci.values()), label=\"FCI\", marker=\"x\", markersize=3)\n",
    "plt.legend()\n",
    "plt.show()\n",
    "plt.figure()\n",
    "\n",
    "plt.plot(list(meta_opt.keys()), list(meta_opt.values()), label=\"opt-meta-VQE\", marker=\"o\")\n",
    "plt.plot(list(fci.keys()), list(fci.values()), label=\"FCI\", marker=\"x\", markersize=3)\n",
    "plt.legend()\n",
    "plt.show()\n",
    "plt.figure()\n",
    "\n",
    "# plot specfic angles\n",
    "key = tq.assign_variable(\"D_{}_{}_{}\".format(0, 0, 1))\n",
    "angles = [angles[key] for angles in meta_angles.values()]\n",
    "plt.plot(list(meta_angles.keys()), angles, label=key.name, marker=\"o\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The example from the Paper:\n",
    "In the paper we used the rectangular H$_4$ system with STO-3G as basis_set.  \n",
    "Below we give the structural data for this system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "geomstring = \"H 0.0 0.0 0.0\\nH 0.0 0.0 1.23\\nH {R} 0.0 0.0\\nH {R} 0.0 1.23\"\n",
    "include_singles=True\n",
    "training_points = [0.5, 1.0, 1.5, 2.0, 2.5]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tequila-3.6",
   "language": "python",
   "name": "tequila-3.6"
  },
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}