Projected Variational Quantum Dynamics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Projected Variational Quantum Dynamics\n",
    "\n",
    "The projected Variational Quantum Dynamics (p-VQD) algorithm is a quantum algorithm for real time evolution. It's a variational algorithm that projects the state at time $t + \\Delta_t$, as calculated with Trotterization, onto a parameterized quantum circuit. \n",
    "\n",
    "For a quantum state $|\\phi(\\theta)\\rangle = U(\\theta)|0\\rangle$ constructed by a parameterized quantum circuit $U(\\theta)$ and a Hamiltonian $H$, the update rule can be written as\n",
    "\n",
    "$$\n",
    "\\theta_{n+1} = \\theta_n + \\arg\\min_{\\delta\\theta} 1 - |\\langle\\phi(\\theta_n + \\delta\\theta)|e^{-i\\Delta_t H}|\\phi(\\theta_n)\\rangle|^2,\n",
    "$$\n",
    "\n",
    "where $e^{-i\\Delta_t H}$ is calculated with a Trotter expansion (using e.g. the `PauliEvolutionGate` in Qiskit!).\n",
    "\n",
    "The following tutorial explores the p-VQD algorithm in Qiskit, which is implemented available as `qiskit.algorithms.time_evolvers.PVQD`. For details on the algorithm, see the original paper: [Barison et al. Quantum 5, 512 (2021)](https://quantum-journal.org/papers/q-2021-07-28-512/#)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The example we're looking at is the time evolution of the $|00\\rangle$ state under the Hamiltonian\n",
    "\n",
    "$$\n",
    "    H = 0.1 Z_1 Z_2 + X_1 + X_2,\n",
    "$$\n",
    "\n",
    "which is an Ising Hamiltonian on two neighboring spins, up to a time $T=1$ and we want to keep track of the total magnetization $M = Z_1 + Z_2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.quantum_info import SparsePauliOp\n",
    "\n",
    "final_time = 1\n",
    "hamiltonian = SparsePauliOp.from_sparse_list([\n",
    "    (\"ZZ\", [0, 1], 0.1), (\"X\", [0], 1), (\"X\", [1], 1),\n",
    "], num_qubits=2)\n",
    "observable = SparsePauliOp([\"ZI\", \"IZ\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to choose the parameterized ansatz we project the update steps onto. We have different choices here, but for real time evolution an ansatz that contains building blocks of the evolved Hamiltonian usually performs very well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAB7CAYAAAB9yZWUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAASZklEQVR4nO3de1RVdd7H8fc5clNBE1FpvCEjN7l4IU00Q7OU7MmecdSy0UmZFcqYuSZpuRw1x4cpS1ktbRopa/JSWhNoT+RtplROJmY6olwUQQWJHvAGKCAK55z9/EGiW25Hg7M38n2tddaS396b/T3yPZ+z9z5n721QFEVBCCF+ZtS6ACGEvkgoCCFUJBSEECoSCkIIFQkFIYSKhIIQQkVCQQihIqEghFCRUBBCqEgoCCFUJBSEECoSCkIIFQkFIYSKhIIQQkVCQQihIqEghFCRUBBCqEgoCCFUJBSEECoSCkIIFQkFIYSKhIIQQkVCQQihIqEghFCRUBBCqEgoCCFUHLQuQO8OLV1PcWauJut2D+zHw7GzNFm3sJ1WPdJS/SGh0ITizFzOHzyhdRlCx+63HpHdByGEioSCEEJFQkEIoSKhIIRQkQONzSBi63K6hfpiNZtRLFbK8y9wfM02zm0/qHVpQidaU49IKDST46sTSVu9FUM7IwGRTxK+dj5fZORSllekdWlCJ1pLj8juQzNTLFayP/kGo6MD7oFeWpcjdEjvPSKh0MyMjg74vTAegKtnCzWuRuiR3ntEdh+aScjLkwiaMxEHVxeUagsHXllLyclzALh5efLo2vnsnLgExWwh+OVJAKS/s03Lku1KURSqrlSgmC04u7thMLa996PW0iO6/stYrVbi4uLw8fHBxcWFgQMHYjKZ8PPzIyoqSuvyVNLe2cYW/xf4LDCSgr2peD4SXDutLK+I/N2HCZ77DG5envSJGErG2i81rNZ+FEUh59M9JI2N4dOAmXwW/AcSHppD2pqtmK/d0Lo8u2otPaLrUIiMjCQ2NpbZs2eza9cupk6dyrRp0zh79iyhoaFal1evqisVHFgQT6/HBtN7/NDa8Yy1X9J7/FAe/ft8fli6HsVs0bBK+1AUhZQF8Rx4JZ6SU/m149eKSjj65qfsnryM6opKDSvUht57RLehsGXLFjZu3EhSUhIxMTGMGTOGxYsXExYWhtlsrg2F8+fPM27cODp06MDAgQNJTU3VuHKoKi0nc912hix6HgwGABSzhfOHTuLo6sLF/2RrXKF9nEkwkfPp3pofrMqtCUrNvy+lnuboii0aVKY9PfeIbkNhxYoVREREEB4erhrv378/jo6OBAfXbHpFR0fj7+/P5cuXmTt3LpMnT8Zi0f5d+OQHO+jQ/QH6T6mp/wHfXnR/yI+ilEz6P/eYxtXZx4kPd4DR0Og8OZ/upbq87W0tgH57xKAoitL0bPZVUFBA7969+cc//kFkZKRq2rRp08jKyiI1NZWysjI8PDz46aef8PDwAMDLy4vNmzczcuTIRtdhMDTerDctdH8Uf6du9/ZEbhPxxf/ww2vrKcst4skv/8q/pi7nxuWrjS6TVXWRt4q//cXr1oKzoR3v9fhvm+Z987KJU9WXWragFqRVj9xtf9j6UtfllkJBQQEAnp6eqvHKykpMJlPtrkNOTg5du3atDQSA4OBgTpzQ12ms/jMjuHTsNMXpuVSXV3L87QSG/eUFrctqUUZsC10Ao40BfT/TU4/oMhRuvsizs9X7VStXrqSwsJAhQ4YAUFFRQadOnVTzdOrUifLy8ibXoSiKTY87d1/uRdaG3RxZvqn253M7vmf/vL81uVx4eLjNdertUWGpwrVP99r95YYYHNrxn59yNK/3lzy06pG77Q9b6fJ7Ct7e3oSEhPDGG2/g7u5Oz549SUxMZOfOnQC1WwodO3akrKxMtezVq1dxdXW1e81CzWAwEDDrSQ4v39jwTEYD/Z4Oo71HZ/sVJpqkyy0Fo9FIQkICgYGBREdHM2vWLDw8PJg7dy4ODg6EhIQA4OPjw6VLl7h8+XLtshkZGQwYMECr0sVt/GdF0COsgb+F0UBHz6489Nrv7VuUaJIuQwHA19eXffv2UVFRQX5+PrGxsaSnpxMQEED79u0BcHNz46mnniI2Npbr16/z4YcfYjAYGD58uMbVC4B2zo48sXkxQXOfwdGtQ+240dGBX08axVM7V9DB013DCkV9dLn70JAjR47UecHHx8czffp0unTpgo+PD1u3bqVdu3YaVSju5NDemYeWzGDQgql84v07AKYe+wAXdzeNKxMN0e2Wwp3Ky8vJzs6uPch4U48ePfj666+prKwkLS2NwYMH26Ue117d8BwZhGuvbox692XVNM+wQCYfjsdr4ggAvCeNYkLS64zdtAhH15qtnBFxc5iU0vTBxvuFQ3vn2n+3lUC4mx6J+OJ/iNi2nHGfL8Ola83Bc616pNWEgqurKxaLhXnz5mldCgCuvbvz4CNBDU4/k2giLykFg0M7/H4/jl2/WcqZRBO+M54AICXmPSovltqpWqEFW3sE4F9TlrN70jLOJCTz66mjAe16pFXtPuiJ7/TH6T7Mn26hfhgMMPqDBbj26c7emW+p5uvs/SAlJ8+hWKwU7k8jbNUcjSoW9mZrjwC15zk4tHem9NSP9i5VpdVsKehN9iffcCbRRMora3Hq7Epy1NuceH87fZ9SH/Nw6tyR6rKar/FWXb2Gc+eOWpQrNGBrjwB07OnBhK9ex39WBCVZ+fX8NvuRUGgGpTkFoChcKyzGqZP6RV91pQJHt5rjCI5uHai6ek2LEoXGGusRgIqfLrHz6cUci/snQdETNajwFgmFe2Q1mzG0+/m/77Zvi915TsWVs4U84N8Hg9HIr0YFt5kzJIXtPWJwuPVpWVVZJZbKKrvU1xA5pnCPSrJ+ZMiff4dbnx5YGznvXTFbyN78DU/+byw3rpTz7R/X2LFKoSVbe6RDjy6M+ts8FKuCtaqa7+b/3Y5V1iWhcI+qy66x+zevqcaKDmZSdDCTbkN86P14KCUn88lLSuFs4recTVSfzTYibg7o7vxU0Zzupkd2T1pWZ3mtekRCoQVcPJpD0hOvNjpPSsx7dqpG6JGee0RCoQnugf3a5LqF7bT6O7XUenV5kRVxf9rw4GQAZhYmalyJaIx8+iCEUJFQELqxbt06Ro8ezejRowkPD8fJyYk1a9bUGauoqFAtFxcXR2pqKnl5eezdW3Oh2OTkZJYsWaKab8OGDfj5+WEymTh06BAjRoxg1KhR/OlPf6qd5/Tp0wwaNIglS5Zw/PhxVq1a1fJPXGckFIRuREVFkZycTHJyMlOmTGHhwoXMnz+/zljHjre+/GO1Wjlw4ACDBw9WhUJDXn31VcLDw+nbty979+5l//79XLhwgfT0dKDmwsCrV68GYODAgaSkpGC1WlvsOeuRhILQndzcXDZv3szSpUsbHQM4fvw4vXr1Amq2ND7++GPGjh0LQHp6Ok8//TQjR46sc4k+T09PXFxcAHBwcGjwdHsfHx+OHTvWXE+tVZBQELqiKAqzZ8/m3XffxcnJqcGxm3JycvDy8gJqtjRmzJjBnj17aqd/9dVXTJgwQTV2u7S0NC5dutTg1bq8vb3JyspqhmfWeshHkkJX4uPjGTp0qOoOYPWN3dTYh2dBQTWnLffs2ZPS0tI604uLi3nppZf4/PPPG/wdiqLYfDuA+4WEgtCNvLw8Pv74Y0wmU6Njt/P19eW7774DwNHRUXUjoNtfzHeGh9lsZvr06axatarOrQRul5uby3PPPXdPz6e1klAQurFy5UouXrzIuHHjase8vb3rjG3atIk+ffoANQcDf/yx5voDQUFBLFq0iGeffZbo6OhG15WQkMDhw4dZuHAhUHNHsrCwsDrznTp1ikGDBv3Sp9a6KELYyXrP3yrrPX/b7L931apVytGjR5ucLyEhQQkNDVWSk5MbnCcnJ0d5+OGHlbi4OOXYsWPKm2++2ZyltgryjUZhN/KNxtZBPn0QQqhIKAghVCQUhBAq8ulDEw4tXU9xZq4m63YP7MfDsbM0WbewnVY90lL9IaHQhOLMXM4f1Net7YW+3G89IrsPQggVCQUhhIqEghBCRY4pNIOIrcvpFuqL1WxGsVgpz7/A8TXbOLf9oNalCZ1oTT0iodBMjq9OJG31VgztjAREPkn42vl8kZFLWV6R1qUJnWgtPSK7D81MsVjJ/uQbjI4OuAd6aV2O0CG994iEQjMzOjrg98J4AK6eLdS4GqFHeu8R2X1oJiEvTyJozkQcXF1Qqi0ceGUtJSfPAeDm5cmja+ezc+ISFLOF4JcnAZD+zjYtSxZ21lp6RNdbClarlbi4OHx8fHBxcWHgwIGYTCb8/PyIiorSujyVtHe2scX/BT4LjKRgbyqejwTXTivLKyJ/92GC5z6Dm5cnfSKGkrH2Sw2rtS/LjWrObL1127zMddu5UVreyBL3p9bSI7oOhcjISGJjY5k9eza7du1i6tSpTJs2jbNnz9Z7aS49qLpSwYEF8fR6bDC9xw+tHc9Y+yW9xw/l0b/P54el61EaueHo/eTC4Sw+D53N/pfeqR07vGwD/xz0Iqc/T9auMA3pvUd0Gwpbtmxh48aNJCUlERMTw5gxY1i8eDFhYWGYzebaUFi2bBkDBgzAaDSSmKiP8/SrSsvJXLedIYueh58vCaaYLZw/dBJHV5c2czv6K6d/4t/PxVJVUlZnmrXKzHfz3yX/X4c1qEx7eu4R3YbCihUriIiIIDw8XDXev39/HB0dCQ6u2fTy8fFhzZo1DBs2TIsyG3Tygx106P4A/afU1P+Aby+6P+RHUUom/Z97TOPq7CMjPgnz9SoUaz3X8VEUMBhIfeuzRi++ej/Ta4/o8kBjQUEBGRkZqjv33JSfn09gYCDOzs4ATJ8+HYDXX3/drjXebvdv695GvLq8kk8H3DqDbfhbUXz/5w8pyy3iyS//yo9fH+HG5av2LNOurGZLzXGE+gLhJkWh5OQ5SrPy6RLQ137FaaA19YgutxQKCgoA6lxlt7KyEpPJ1CzHEwwGg02Phq4ifDf8Z0Zw6dhpitNzqS6v5PjbCQz7ywtNLmcymWyuU28PN6f2WG9U2/T/88igYZrX+0seWvXI3faHrXS5peDh4QFAdnY2EyZMqB1fuXIlhYWFDBkyRKvS7knWht2qn8/t+J5zO77XqBr7uKGYMSsWHAz133npduXWG3aoSN/01CO63FLw9vYmJCSEN954g02bNrFnzx6io6P56KOPAJplS0FRFJsedx7TsKfw8HCb69Tbw6xY8ZsyBoyNvEMZDHT26UleVYnm9f6Sh1Y9crf9YStdhoLRaCQhIYHAwECio6OZNWsWHh4ezJ07FwcHB0JCQrQuUdgg6I/PYHR0aDgYFIVBMc+2uTsw6Z0udx+g5s4/+/btU43NmDGDgIAA2rdvXztWXV2NxWLBarVSXV3N9evXcXZ2lkbTgS4BfXl80yL2/WEV1eWVGIw170GK1YrBaGRY7Cz6TRyhcZXiTroNhfocOXKE4cOHq8ZefPFFNm7cCMD+/fuBmlt93bzpqNDWrx4NYcrR9zm79VsKv8vAarbQNbgfPs+PpeODXbUuT9RDl7sP9SkvLyc7O7vOQcYNGzbU2XeyRyC49uqG58ggXHt1Y9S7L6umeYYFMvlwPF4/vws+tn4hz2dt5MFRt77WOiJuDpNS/tbideqBk1sH/GdGMObDGMZuWMigBVPbRCDcTY8AuAf3Y2ZhIoZ2NS9LrXqk1YSCq6srFouFefPmaV0KAK69u/PgI0ENTj+TaCIvKQWAgwvf58QHO1TTU2Leo/JiaUuWKDR2Nz0CP38smXam9meteqRV7T7oie/0x+k+zJ9uoX4YDDD6gwW49unO3plv1Zm38kKp/QsUmrubHnnAtxfX/u8ynbwavgO2vbSaLQW9yf7kG84kmkh5ZS1OnV1JjnqbE+9vp+9Tw5teWLQJd9MjA6L+i5Prd2lQZV0SCs2gNKcAFIVrhcU4deqodTlChxrrEbd+nlSVXeNGcd0Tx7QgoXCPrGZz7QEhbvtiiHwUKm6ytUe6BPTFY1B/ntiymC4BfQh7S9trhcgxhXtUkvUjQ/78O9z69MDaxHnvw2Ij6f1EKL3HP8SpTf8m+5Nv7FSl0JKtPZK/8xD5Ow8BNVd9Prhwnb1KrJeEwj2qLrvG7t+8phorOphJ0cFMug3xoffjoZSczCcvKYUfln7ED0s/Us07Im4OtM0zhtuMu+mRm24/m1KrHpFQaAEXj+aQ9MSrjc6TEvOenaoReqTnHpFQaIJ7YL82uW5hO63+Ti21XoPSVi97I4Sol3z6IIRQkVAQQqhIKAghVCQUhBAqEgpCCBUJBSGEioSCEEJFQkEIoSKhIIRQkVAQQqhIKAghVCQUhBAqEgpCCBUJBSGEioSCEEJFQkEIoSKhIIRQkVAQQqj8PwpwFbCvA7WLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 327.252x144.48 with 1 Axes>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiskit.circuit import QuantumCircuit, ParameterVector\n",
    "\n",
    "theta = ParameterVector(\"th\", 5)\n",
    "ansatz = QuantumCircuit(2)\n",
    "ansatz.rx(theta[0], 0)\n",
    "ansatz.rx(theta[1], 1)\n",
    "ansatz.rzz(theta[2], 0, 1)\n",
    "ansatz.rx(theta[3], 0)\n",
    "ansatz.rx(theta[4], 1)\n",
    "\n",
    "# you can try different circuits, like:\n",
    "# from qiskit.circuit.library import EfficientSU2\n",
    "# ansatz = EfficientSU2(2, reps=1)\n",
    "\n",
    "ansatz.draw(\"mpl\", style=\"iqx\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this ansatz, the $|00\\rangle$ state is prepared if all parameters are 0. Hence we'll set the initial parameters to $\\theta_0 = 0$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "initial_parameters = np.zeros(ansatz.num_parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before running we need to select the backend and how we want to calculate the expectation values. Here, we'll use exact statevector simulations (which is still very fast as we investigate a 2 qubit system)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.primitives import Sampler, Estimator\n",
    "from qiskit.algorithms.state_fidelities import ComputeUncompute\n",
    "\n",
    "# the fidelity is used to evaluate the objective: the overlap of the variational form and the trotter step\n",
    "sampler = Sampler()\n",
    "fidelity = ComputeUncompute(sampler)\n",
    "\n",
    "# the estimator is used to evaluate the observables\n",
    "estimator = Estimator()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since p-VQD performs a classical optimization in each timestep to determine the best parameters for the projection, we also have to specify the classical optimizer. As a first example we're using BFGS, which typically works well in statevector simulations, but later we can switch to gradient descent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.algorithms.optimizers import L_BFGS_B\n",
    "\n",
    "bfgs = L_BFGS_B()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can define p-VQD and execute it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.algorithms.time_evolvers.pvqd import PVQD\n",
    "\n",
    "pvqd = PVQD(\n",
    "    fidelity,\n",
    "    ansatz, \n",
    "    initial_parameters, \n",
    "    estimator=estimator,\n",
    "    num_timesteps=100, \n",
    "    optimizer=bfgs\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The p-VQD implementation follows Qiskit's time evolution interface, thus we pack all information of the evolution problem into an input class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.algorithms.time_evolvers.time_evolution_problem import TimeEvolutionProblem\n",
    "\n",
    "problem = TimeEvolutionProblem(hamiltonian, time=final_time, aux_operators=[hamiltonian, observable])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then run the algorithm!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = pvqd.evolve(problem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can have a look at the results, which are stored in a `PVQDResult` object. This class has the fields\n",
    "\n",
    "* `evolved_state`: The quantum circuit with the parameters at the final evolution time.\n",
    "* `times`: The timesteps of the time integration. At these times we have the parameter values and evaluated the observables.\n",
    "* `parameters`: The parameter values at each timestep.\n",
    "* `observables`: The observable values at each timestep.\n",
    "* `fidelities`: The fidelity of projecting the Trotter timestep onto the variational form at each timestep.\n",
    "* `estimated_error`: The estimated error as product of all fidelities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The energy should be constant in a real time evolution. However, we are projecting the time-evolved state onto a variational form which might violate this rule. Ideally the energy is still more or less constant. In this evolution here we observe shifts of ~5% of the energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Energy over time')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEYCAYAAAB2qXBEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAw7UlEQVR4nO3dd5iU1dnH8e89s3W2so2OoCJFKcIiKAoqiqgoGDWKXaO8Ro1RkxjTTGI0idEUjSYGS9TYWxS7KBZEEXYRkCKKWOh1ey/n/WNGnays7DJtZ/f3ua65duap99ndmXvOOc9zjjnnEBER2VOeWAcgIiLxTYlERERCokQiIiIhUSIREZGQKJGIiEhIlEhERCQkSiQigpndYWa/inUcEp9M95FIvDOzz4DuQFPQ4nudc5fFJqKOzczOAy50zh0a61ikc0iIdQAiYXKCc+7VSJ7AzBKcc42RPEe4xWPMEn/UtCWdmpmdZ2Zvm9nNZlZiZp+a2bFB67PM7G4z22RmG8zsejPzBu0738z+amY7gd+YWa6ZPWtm5Wa2KLD924HtbzezP7c4/7NmdkUrsR0SOEZZ4OchgeWnm1lRi22vNLPZgefJgfJ8YWZbAs1SqYF1h5vZejP7qZltBv7d4jhDgDuAg82s0sxKA8vvNbPrWxzjajPbGvjdTDez48zsIzPbaWY/Dzqmx8yuMbNPzGyHmT1mZjnt/2tJvFIika5gLLAayAP+BNxtZhZYdx/QCOwLHAhMBi5sse9aoAC4AbgdqAJ6AOcGHgQda4aZeQDMLA+YBDzcMqDAB+3zwK1ALvAX4HkzywVmA4PMbGDQLmcADwWe3wjsB4wMxN0buDZo2x5ADrAXMDP4vM65VcDFwLvOuXTnXPY3f11fHSMl6Nh3AmcBo4HDgGvNbO/AtpcD04GJQC+gJPB7kq7COaeHHnH9AD4DKoHSoMdFgXXnAWuCtvUBDv8HZXegDkgNWj8DeD1o3y+C1nmBBmBQ0LLrgbeDXq8Cjg48vwx4oZWYzwYWtlj2LnBe4PkDwLWB5wOBikDshj+R7RO038HAp4HnhwP1QMq3/L7OC445sOxe4PqgY9QA3sDrjMDvbGzQ9sXA9KAyTwpa1zPwe0qI9f+GHtF5qI9EOovprvU+ks1fPnHOVQcqI+n4v7UnApu+rqDgAdYF7Rv8PB9/v2Jr68FfKzkLmBP4eUsrMfUCPm+x7HP8NQDw1z7+DFyHvzbydCD2AvwJpTgoZsOf5L60zTlX28p522qHc+7LixdqAj+3BK2vwf87BH/N579m1hy0vgl/ot4QYhwSB5RIpCtbh79Gkuda75AOvqxxG/5msD7AR4FlfVts/wCw3MxGAEOAp1s57kb8H8DB+gEvBZ6/AuSZ2Uj8taQrA8u34/8Q398519qH9O4uxQz3pZrrgAucc/PDfFyJE+ojkS7LObcJ/wf2n80sM9BpvI+ZTWxl+ybgKfyd7j4zGwyc02Kb9cAi4D/Ak865mm8eCYAXgP3M7AwzSzCz04ChwHOB4zQCTwA34a85zQksb8bfX/HXQO0EM+ttZse0o+hbgD5mltSOfb7NHcANZrZXIJ58M5sWpmNLHFAikc7i2cBVSF8+/tvG/c4BkoCV+DuJn8Dfxt+ay4As/M1l/8HfkV7XYpv7gGGB9bvknNsBTAV+BOwArgamOue2B232EHAU8HiLGtNPgTXAAjMrB14FBn1rKf/XXGAFsNnMtu9u4za4Bf8FAq+YWQWwAP9FCtJF6IZEkRCY2Y1AD+fcuUHLJuBv4uofqEGIdGqqkYi0g5kNNrPh5ncQ8D3gv0HrE4EfAncpiUhXoUQi0j4Z+PtJqoDH8F9Z9Qx8dbNfKf6msb/FJjyR6FPTloiIhEQ1EhERCUmXvI8kLy/P9e/fP9ZhiIjEleLi4u3OufyWy7tkIunfvz9FRUW731BERL5iZi1HYwDUtCUiIiFSIhERkZAokYiISEiUSEREJCRKJCIiEhIlEhERCYkSiYiIhKRL3kciEi8amxwbtzWycXsjVTXNVNc2U13raGxyNAWGhOydn8DevRPp2z2RBK99+wFFIkCJRKQDaWxyLPmolvlLa1i8upYN2xppbuMYwokJMGJgCuMOSGXcsFR65entLdGh/zSRDmLZmlp+dcd2KqqbSU40Rg1OYcKBPvoWJNArP5HMNA++FCM12UNiguH1QlMTrNvSwNoNDaz+op6FK2q47fESbnu8hDFDUzhzSibD902JddGkk+uSo/8WFhY6DZEiHYlzjh/cvIXtpU384LRujB6cQkrSnnVhbtjawBuLq3nq9QpKKpoZtm8yl5yczaC9ksMctXQ1ZlbsnCtsuVyd7SIdwLI1daz8tJ4ZkzMZP9y3x0kEoHdBImdOyeLB3/XislO7sXFbI5fetIW7nymlvqHrfXGUyFMiEekAHn6lnOx0D1MOTgvbMVOSPHzniAz+/aueTB6bxoMvl3PxHzfzyfr6sJ1DBJRIRGJuzbp6Fq6o5eQjM0gOoSbSmnSfh6vPzuUPl+ZTUd3MD27ewhvFVWE/j3RdSiQiMfbwnHJ8Kca0CRkRPc/Y/VP51zU92LdvEtfdvYNZT5fS1KymLgmdEolIDG3Y1sCbxdWceFg66b7Ivx1zsrz8+YcFnHBYOo+8Us4f7t2hZCIh0+W/IjH08MvleL1w8pGZUTtnYoJx5YwceuQmcOfTpaQkGVedkYPHo5sZZc8okYjEyIZtDby0oIrpE9LJzfJG/fwzJmdSW9fMf14sJyXJuPTUbpgpmUj7KZGIxMgDL5aT4DVmHJMVsxjOm5pFTZ3jibkVdMv0cuaU2MUi8SuqfSRmNsXMVpvZGjO7ZhfrB5vZu2ZWZ2Y/bsu+ZnaTmX1oZsvM7L9mlh2FooiEZN2WBua8V8W0GNVGvmRmfP/kbI4a4+OeZ8soWlUTs1gkfkUtkZiZF7gdOBYYCswws6EtNtsJXA7c3I595wAHOOeGAx8BP4tYIUTC5D8vlJGUaJx2dPT6RlpjZlx5Rg579Ujk+nt2sGVnY6xDkjgTzRrJQcAa59xa51w98AgwLXgD59xW59wioKGt+zrnXnHOffmfvwDoE8lCiITq800NvFZUzfSJ6eRkxq42Eiw12cNvZubR2OS47q7tNDTqSi5pu2gmkt7AuqDX6wPLwrnvBcCLexSdSJTcPdt/pdR3j4p9bSRYv+6JXH12Lqs+q+fu2aWxDkfiSDQTya4uB2nr157d7mtmvwAagQd3eQCzmWZWZGZF27Zta+NpRcJr6Ue1vL20hhmTM8nO6Bi1kWATDvQx9dB0nphbwcfrNJSKtE00E8l6oG/Q6z7AxnDsa2bnAlOBM10rwxk752Y55wqdc4X5+fntClwkHJqbHf94soSCbl5OnRTZu9hDcdH0bLLSPfzloZ26WVHaJJqJZBEw0MwGmFkScDowO9R9zWwK8FPgROdcdQTiFgmLVxdW8fG6Bi6clh2RMbXCJcPn4dJTurH683qeebMy1uFIHIjaf3OgQ/wy4GVgFfCYc26FmV1sZhcDmFkPM1sPXAX80szWm1lma/sGDn0bkAHMMbMlZnZHtMok0la19c3cNbuMQXslcWShL9bh7NYRo32MGZrCPc+Wsq1EV3HJt9PEViJRcM+zpTzwYjm3XFXAsDiZsXDj9kYu+N0mxo9I5VcX5MU6HOkANLGVSIys+qyOh14uZ/LYtLhJIgC98hI4dVIGrxdV89EX6niX1imRiERQbX0zf7h3B3lZXi77brdYh9Nupx2dSWaah7ueKY11KNKBKZGIRNCs/5ayfmsjV5+TS3pq/L3d0lM9nHFMJkWralm8ujbW4UgHFX//2SJxomhVDU+/WcnJR2YwalD8NGm1NH1iBgXdvNz5dCldsU9Vdk+JRCQC6hscf324hD4FCVx4YnyPqJuUaJw7NYvVn9czb4kGdZRvUiIRiYBH55SzaXsjPzw9p0PfM9JWk8emsVePBO5/vky1EvmG+P8PF+lgNu9o5MGXy5lwYCqjB8dvk1Ywr8c/UvHajQ0UrVJfifwvJRKRMLv9iRI8BpecHH9XaX2bSWPSyM3y8uic8liHIh2MEolIGC1cUcP8pTWcNSWTgpzONQFpYoJx8hEZLF5dp/tK5H8okYiESU1tM397ZCd9ChI4ZVLHGiI+XKYemk5qsvH4a6qVyNeUSETCZNYzpWzZ2cRPzsohKXFXMx/Ev3Sfh6mHpvN6cTWbd2gMLvFTIhEJg6Uf1fLMm5WcNDE9roZB2RMnH5GBAU++XhHrUKSDUCIRCVFtfTM3PbiTnnkJfG9adqzDibiCnAQmjPLxyoIq6ht0KbAokYiE7I4nS9m4rZGfnJVDanLXeEsdd0g6FdXNzFuiKYBEiUQkJM/Oq2D2vEq+e1QGI/fr3E1awQ7cL5keuV5eercq1qFIB6BEIrKH3l9dy62PlnDQ/ilcND071uFElcdjTDk4neIPa9XpLkokIntiw9YGfnvXdvoUJPDLC/LwejrnVVrf5phxaZjBS+9qOt6uTolEpJ1q65u5dtZ2AK7/fn5cDg8fDt1zEigcksJL71bR1KxO966sa74DRPaQc45bHinhs00N/OL8XHrnJ8Y6pJg69pB0tpY0sfhDjb/VlXWuMRxEIuzFd6p4eUEV5xyXyZihqbEOJ+YOGZZKZpqHG/69g4IcL75kD2mpHjLTPKT7PPQtSGDiKB9Z6d5YhyoRpEQi0kafrK/n1sdKGD04hbOPi+85RsIlKdG4ckYO85ZUU13bTHWtY8vORj5Z30xFdTM1dY7bHi/h4GGpHD8+nTFDUzDrev1JnZ0SiUgblFU2ce2s7WT4PPz8/Nwu2bnemomjfEwc5dvlujXr6nn5vSpeW1jFvCU1DNs3mZnTs9l/7+QoRymRZF1xkprCwkJXVFQU6zAkTjQ0Oq7++1ZWflrHX6/sztAB+hBsr8Ymx/PzK7n/hTJKypsZPyKVC07IYkCvpFiHJu1gZsXOucKWy1UjEfkWzjn+/lgJSz+u42fn5iqJ7KEErzFtQgaTD0rjibkVPPpqORcuq2HSGB/nHp/V5S9aiHe6akvkWzz1egXPvV3JjMmZHD02LdbhxL3UFA9nH5fFg9f14rSjMpj3fg0XXr+Zt97XUCvxTIlEpBWPvVrO7U+UMn54Kt87UZ3r4ZSV7mXmSd24/7c92adPIr+5czsPv1Ku+eDjlBKJSAvOOe5+ppQ7nipl4igf116Yh0ed6xGRn53An39YwBGFPu58upQ/P7hTNzfGIfWRiAT5sk/k6TcrOX58GlfMyNEVWhGWnOThl+fn0js/gQde9M+8eNUZOUrecUSJJAwWrazhgRfL6ZWfwF49Ezlg72QO2EedsvHo/hfKefrNSk6dlMHF38nWPQ9RYmZccEI2BvznxXJSkj1ceop+//FCiSQMHp1Tzicb6tm4vZGXF/iH1f75ebkcdZA6Z+PJs/MquO/5Mo4Zl6YkEiPnTc2ius7x5NwKfCn+5CIdnxJJiHaWNbHkozrOnJLJ+SdkU17VxK9nbeemB3aQn+1lRBeaoyKevb2kmlseKWHs/in86MwcJZEYMTMuOTmb2rpmHnixnN75CRwzLj3WYcluqLM9RG8srqbZwRGF/tpHZpqX6/4vn565CVw7aztfbG6IcYSyO2+9X811d29nv72SuPbCPBK8SiKxZGZccXoOBw5K5i8P7eSjL+pjHZLshhJJiOYWVbF370T69/z6hqoMn4c/XFpAghd+dvtWauqaYxihfJuXF1Ry3V3b2a9fEjdeVtBlpsrt6Lxe41cX5NEt08u1s7ZRWtEU65DkW+hdE4LNOxpZ+Wk9RxZ+c5yhnnn+CY827WjixXc0HWlH9PSbFdx4/04OHJTCTT8oIMOnt0NHkp3h5bqZ+ZSUN/G7e7brsuAOTO+cEMwt8ieIIwt33al+4KAUhu2TzGOvldPYpDdBR/LU6xXc+mgJ44encsP380lN0VuhI9qvXxJXzMjh/dV1PP5qRazDkVbo3ROC14urGTogiR65rV+zMGNyJlt3NjG3SENAdBRPvV7BbY+XcNjIVH59UR5JieoT6cimjEvjsJGp/Pu5Uj7dqP6SjkiJZA99vqmBT9Y3tFob+dLYA1IY0CuRR14pp1lV85h7Ym75V0nkV99Tx3o8MDOumJGDL8XDH+/bodp9BxTVRGJmU8xstZmtMbNrdrF+sJm9a2Z1ZvbjtuxrZjlmNsfMPg787BbpcjQ1O+5/oQwzWp2HISg+Tj86k882NbBgeU2kQ5NWNDQ6/vbwTv7xRKmSSBzqluHlyhk5fLyugQdfKo91ONJC1BKJmXmB24FjgaHADDMb2mKzncDlwM3t2Pca4DXn3EDgtcDriGlscvz+3zt4vbia752QRW7W7qcQPaLQR/ccLw+9XE59g75NRdvOsiZ+dMtWZs+r5PSjM3SJb5yacKCPSWN8PPBimZq4Opho1kgOAtY459Y65+qBR4BpwRs457Y65xYBLW+++LZ9pwH3BZ7fB0yPUPzUNzh+e+d2Xi+u5uLvZHPGlLaNCJvgNWZMzmTlp/Wc/NP1/OHe7Sz4oEZXoURYc7PjpXcruegPm1izrp5fXZDLzJO6aeysOHbZqd3wpXj4+2MlGim4A4nmne29gXVBr9cDY8Owb3fn3CYA59wmMysINdDW/PXhncxfVsPlp3Vj+sSMdu17wmHp9MxL4PXiat5eUs2chdX0zk/glEkZTBmXRnKSuqvCxTnH0o/ruOOpUj76op4h/ZO46owc9umj2fjiXVa6lwtOyOKWR0t46/2a3TYtS3REM5Hs6mtgW79ShLKv/wBmM4GZAP369WvPrl+ZMTmTAwelMHkPJjgyM8YMTWXM0FSunJHD/GU1PDqnnFseKeG+58o49/gsph6ajldNLntsw7YGXl1YzWuLqli/tZH8bC+/OD+XI0b7NJJsJzL1sHSem1/JP58qYewBKaToS1jMRTORrAf6Br3uA2wMw75bzKxnoDbSE9i6qwM452YBs8A/Z3t7Av9Svx6J9OsR+pSgiQnG4aN8TDwwlWVr6rjvuTJuebSE2fMqueSUbowerPG52qO6tpl/P1fGf1+vwAEjBiZz+tGZHDnGpw+ZTsjrMS47tRtX/nUrj7xSznlTs2MdUpcXzUSyCBhoZgOADcDpwBlh2Hc2cC7wx8DPZ8IZdCSZGSMGpvDnK5J5e2kNdzxVyk9u3cox49K49NRupKfqQ/DbOOd454Ma/v5oCdtKmzjh0HTOnJJJfjeNRdrZjRiYwpGFPh5+pZwpB6d/671cEnkWzQ4rMzsO+BvgBe5xzt1gZhcDOOfuMLMeQBGQCTQDlcBQ51z5rvYNHDMXeAzoB3wBnOqc2/ltcRQWFrqioqIIlDA09Q2OB14q46GXy8nL8nL1ObmMGqTaya6sWVfPv/5bSvGHtQzolchVZ+Sw/96aA6Yr2bqzkbN/s5GjD0rjx2flxjqcLsHMip1zhd9Y3hWvfOioieRLqz6r4w/37mD91kbOnJLJeVOzdKVRwMbtjdz/fBlzFlaR4fNw9rGZTJuYoct5u6jbHi/h6TcruPfanvQpCL3ZWb5da4lEbScd0JD+ycz6eQ+OPTiNB18q5xf/2EZFddceQXjj9kZu+s8OzvnNRl4vruK7kzJ44Le9OPnITCWRLuyMyZkkJRj3P18W61C6NDUsdlApSR5+fFYOQwYkc+ujO7n4j5u54fv5/zNcfVewvbSR/7xQzvPvVOL1wPQJ6Zw+OZO8bP3rCuRkeZl+eAaPzilnxjH1DOilS7xjQTWSDszMmHpoOn+7qjt19c1cfvNmlq2pjXVYUbGtpJE7ny7l7F9v4oV3Kjnh0HQevK4Xl303R0lE/sdpR2WQmmzc+5xqJbGid2QcGDogmdt+0oNrbtvKT27dys/Py+uUN2LV1DUzb0kNc96rYvHqWpyDSWN8nDc1i975XasmJm2Xle7llCMzuP+Fctasq2ffvqqVRJs62+NIeVUTv7xjOyvW1nHF6d044bD23V3fETnn+OiLel54p4q5i6qoqnX0zPVy9Ng0jh6bpgQibVJZ3cxpv9jAIcNT+cX5ebEOp9NqrbNdNZI4kpnm5aYf5PPbu7bz14dLaGyCkw6Pz2SycXsjcxdV8dqiKj7f3EhyojFxlI9jD0lj+L7JmKkDXdou3edh6qHpPPl6Bd87sVH3lUSZfttxJjnJw29n5nPd3dv5+2MlNDY5Tp2UGeuwdmvdlgYWrqhh1Wf1rPqsnk3bGwEYtm8yV5yewZGFaaRrqlsJwclHZPDU6xU8+XoFl54S8dkkJIgSSRxKTDB+fWEeN/x7B/98shSgwyUT5xyfb25k/tJq3lxczZr1/gGd87O9DBmQxLQJ6Uw40KdvjhI2BTkJHFHo4/n5lZxzXBYZ+mISNXoXx6kEr/HL8/138/7zyVJ8KR6OH58e05hq65t5f3UdC5bX8N6KGrbubAJg6IAkLjklm8NG+uieo385iZzTjsrk1YXVPPtWRZuneZDQ6V0dx7xe4+fn5VJT18xfHtqJL9k4YjdT/4bb5h2NvLeihveW17B4dR31DY6UZGP0oBTOmpLKQUNTKFDykCjZp08ShUNSeOqNCk6ZlElSovraokHv8DiXmGD85qI8fnrbNn5/7w5Wf1HPwcNS2X/v5Ijc8V1W2cTSj+tY9nEtxavr+HyTv8mqZ14CU8enMW5YKsP3TdEbWGLmtKMz+cmtW3mjuIrJ42JbS+8qlEg6gZQkD7//fj5/vG8HT86t4LFXK0hL9Y8sPGyfZIbvm8w+fZL26MO9qdmx+vN6FiyvYcHyGtas8yeOpERj2D7JHHdIGmMPSKVvQYKutJIOYdSgZPp1T2D2vEolkijZbSIxs/2dcyuiEYzsubRUD7+7OJ+qmmaKP6xl4Yoaln5cxzvLagDweKBvQQIDeiXRMy+B3CwvedleUpMNr9fweqC23lFW0URpZTMbtjbyyYZ61m5ooLbe4THYf+9kLjghi5H7pTBoryQSE5Q4pOMxM044LJ3bnyjVDYpRstsbEs1ssXNuVOD5hc65u4LW+Zxz1RGOMezi9YbEPbGjrInln9SxZn09n21sYO3GBraVNNLY9O37Zfg87NM7kb37JDGkfxJjhqaQmeaNTtAiIaqobua7P9vA5LFpXHlGTqzD6TRCuSEx+GvnJcBdQa/nAaNDjE0iKDfLy8RRvv8ZUqW52VFe1cyOsiZq6hxNzY6mJkhOMrLTPWRleElLMTVVSdzK8Hk4otDHnEVVzDwpmzRNEhdRbUkkwVWWlp8s+uvEIY/HyM7wkp2hGoZ0Xiccls5L71bx6sIqpk2MzxEg4kVbEkEPMzvPzA7km4mk6w3UJSJxYfBeSQzsm8iz8yrpimMKRlNbEslvgEL809z2MbMVZvakmd0AaHQ0EemQzIwTJ2SwdmMDK9bWxzqcTm23icQ5N8s5d5lzbqJzLg84BrgbKAfeinSAIiJ76sjRPlKSjZcXVMY6lE6t3X0czrn1zrkXnHM3OufOikRQIiLhkJriYcJIH28UV1NX37Wnq44kdZaLSKc2eWwaVbXuq3uqJPyUSESkUxu5XzIF3by8/F5VrEPptNqcSMzsMjPTIP8iElc8HuOog9IoWlnLzrLd3Ikre6Q9NZIewCIze8zMppjuVhORODF5bBrNDl4rUq0kEtqcSJxzvwQG4r9i6zzgYzP7vZntE6HYRETCol+PRAb3T+LlBUokkdCuPhLnv6tnc+DRCHQDnjCzP0UgNhGRsJk8No21GxpYs073lIRbe/pILjezYuBPwHxgmHPu+/jH2jo5QvGJiITFEaN9eD0wtzjuxpnt8NpTI8kDvuOcO8Y597hzrgHAOdcMTI1IdCIiYZKV7mXU4BTeLK7SkClh1p4+kmudc5+3sm5V+EISEYmMw0f52LSjiY++UPNWOLV5hkQzu2oXi8uAYufckrBFJCISIeNHpPKXh+DNxdUM2is51uF0Gu1p2ioELgZ6Bx4zgcOBO83s6vCHJiISXplpXkYPSeGNxdVq3gqj9iSSXGCUc+5Hzrkf4U8s+cAE/JcDi4h0eIeP8rF5RxOrP1fzVri0J5H0A4J/8w3AXs65GqAurFGJiETI+BE+ErzwxmJdvRUu7UkkDwELzOzXZvZr/JcAP2xmacDKiEQnIhJmGT4Powen8Kaat8KmTYkkMBzKvcBFQCn+TvaLnXPXOeeqnHNnRixCEZEwO3y0jy07m/jwMzVvhUObrtpyzjkze9o5NxoojnBMIiIRNX64D69nJ/OWVDNkgK7eClV7mrYWmNmYiEUiIhIl6T4PI/dLYb7mKAmL9iSSI/Ank0/MbJmZfWBmyyIVmIhIJB0yPJV1Wxr5YktDrEOJe+1JJMcCewNHAifgHxblhPacLDD8/GozW2Nm1+xivZnZrYH1y8xsVNC6H5rZcjNbYWZXBC0faWYLzGyJmRWZ2UHtiUlEuqZDhqcC8M5S1UpC1Z5E8gVwGHBuYKgUB3Rv685m5gVux5+QhgIzzGxoi82OxT9U/UD8Nzz+M7DvAfg7+g8CRgBTzWxgYJ8/Ab91zo0Erg28FhH5Vt1zEti3byLzl+ky4FC1J5H8AzgYmBF4XYE/MbTVQcAa59xa51w98AgwrcU204D7nd8CINvMegJDgAXOuWrnXCPwJnBSYB8HZAaeZwEb2xGTiHRh44f7WPlpPTvLNXNiKNqTSMY65y4FagGccyVAUjv27w2sC3q9PrCsLdssByaYWa6Z+YDjgL6Bba4AbjKzdcDNwM92dXIzmxlo+iratm1bO8IWkc5q/PBUnIMFH6h5KxTtSSQNgeYpB2Bm+UBzO/bf1dS8Le8G2uU2gdGFbwTmAC8BS/FPrAXwfeBK51xf4Er8Mzju6iCznHOFzrnC/Pz8doQtIp3VPn0S6Z7j1dVbIWpPIrkV+C/Q3cxuAN4Gft+O/dfzdS0CoA/fbIZqdRvn3N3OuVHOuQnATuDjwDbnAk8Fnj+OvwlNRGS3zIxDhqdS/GEtNXXt+V4swdozH8mDwNX4k8dGYLpz7vF2nGsRMNDMBphZEnA6MLvFNrOBcwJXb40DypxzmwDMrCDwsx/wHeDhwD4bgYmB50fydYIREdmt8SN81Dc4ilfVxjqUuNWe+UiSgVH4O7QTgFPNDOfcdW3Z3znXaGaXAS8DXuAe59wKM7s4sP4O4AX8/R9rgGrg/KBDPGlmufgHi7w00EcD/qu5bjGzBPz9NzPbWiYRkeH7JuNLMRaurOXQkb5YhxOX2pxIgGcITGTFHo7265x7AX+yCF52R9BzB1zayr6HtbL8bfzzxouItFuC1zhwvxSKV6mfZE+1J5H0cc5NiVgkIiIxMnqIf7iUDdsa6J2fGOtw4k57OtvfMbNhEYtERCRGRg9OAaBI/SR7pD2J5FBgcWCIE421JSKdRp+CBApyvOpw30Ptado6NmJRiIjEkJlRODiFN9+vpqnJ4fXu6pY2aU3UxtoSEenIRg9OoarGsfoLTXbVXtEca0tEpMMaNTgFM9S8tQeiOdaWiEiHlZXuZWDfJIo+VCJpr2iOtSUi0qGNHpzCyrV1VNfqo6099mSsrYI9HGtLRKRDGz0khaZmeODFMt5YXM2ilTV8vqmB+oaW48tKsDZfteWce9DMioFJ+EfpnR4YlVdEpFM4YO9kumV4eGROBf5uYD+PQffcBAb1S2LEfsmMHJhCvx4JmOnqLmjf5b845z4EPoxQLCIiMZWUaDz0u16UVjZTXdtMZXUzW3Y2sX5rA19sbmT52jreWOyfUbFXfgKTx6Zx9EFp9Mxr10dpp9O1Sy8i0kJykofuOV+3+gcP5+GcY+P2Rt5fXcfcoirufa6Me58rY9wBKfzfSd3Yq2fXHF7F/OMkdi2FhYWuqKgo1mGISJzbvKORlxdU8cRr5dTUO6aOT+fcqVl0y/DGOrSIMLNi51xhy+Xt6WwXEZEgPXITOPf4LB64rhcnHpbO8/MrufD6TRR3sUuIlUhEREKUle7l8tNymPXzHmSme7n671u597lSmpq7RouPEomISJgM6JXEP67uzuSxadz/Qjk/u31bl7gnRYlERCSMUpM9/PScXK46I4fFq2u56m9bKaloinVYEaVEIiISAVMPTef6/8vn800NXH7zFjZub4x1SBGjRCIiEiHjhqVy8w8LqKhu5sq/bGFbSedMJkokIiIRtP/eydx8eQFVtc387PZtVNZ0vj4TJRIRkQjbt28Sv7koj883N3Dtv7Z1urG7lEhERKKgcEgqPzk7lyUf1XHzgzvoTDeDa4gUEZEomTw2jS07G/n3s2UM6Z/MSYdnxDqksFCNREQkis48JpNxB6TwzydLWPVpXazDCQslEhGRKPJ4jGvOzSUv28tv79pOWWX832OiRCIiEmWZaV5+fWEeJRVN/OG++O8vUSIREYmBQXsl838ndWPhilpeea8q1uGERIlERCRGpk9MZ/+9k/jHE6XsLI/fJi4lEhGRGPF4jB+fmUttfTO3P14S63D2mBKJiEgM7dUzkbOmZPF6cTXvLKuOdTh7RIlERCTGTp+cyYBeidzySAk1dfE3hIoSiYhIjCUmGFec3o1tpU089mpFrMNpNyUSEZEOYNi+KRw+yscjr5TH3SjBSiQiIh3EzJOyaXaOO58pjXUo7aJEIiLSQfTITeDUSZm8urA6roZPUSIREelAzjgmk26ZHm5/oiRu7nhXIhER6UB8KR7On5rNyk/rWbSyNtbhtElUE4mZTTGz1Wa2xsyu2cV6M7NbA+uXmdmooHU/NLPlZrbCzK5osd8PAsddYWZ/ikJRREQi5phxaRTkeLn/hbK4qJVELZGYmRe4HTgWGArMMLOhLTY7FhgYeMwE/hnY9wDgIuAgYAQw1cwGBtYdAUwDhjvn9gdujnxpREQiJzHBOGNyJis/rWfx6o7fVxLNGslBwBrn3FrnXD3wCP4EEGwacL/zWwBkm1lPYAiwwDlX7ZxrBN4ETgrs833gj865OgDn3NZoFEZEJJKmHJxOXra/VtLRRTOR9AbWBb1eH1jWlm2WAxPMLNfMfMBxQN/ANvsBh5nZe2b2ppmN2dXJzWymmRWZWdG2bdvCUBwRkchJSjROPzqTD9bUsfSjjt1XEs1EYrtY1rLxb5fbOOdWATcCc4CXgKXAl3fsJADdgHHAT4DHzOwbx3HOzXLOFTrnCvPz8/ewCCIi0XP8+DRyMj3c18FrJdFMJOv5uhYB0AfY2NZtnHN3O+dGOecmADuBj4P2eSrQHLYQaAbyIhC/iEhUJSd5+O5RmSz5qI7Vn3fcvpJoJpJFwEAzG2BmScDpwOwW28wGzglcvTUOKHPObQIws4LAz37Ad4CHA/s8DRwZWLcfkARsj3BZRESi4vjx6fhSjCfmdtwxuKKWSAKd5JcBLwOrgMeccyvM7GIzuziw2QvAWmANcCdwSdAhnjSzlcCzwKXOuS8H778H2NvMluPvwD/XxcP1ciIibZCW6uG4Q9J5o7i6w47BZV3xM7ewsNAVFRXFOgwRkTbZvKORs67dyHePzmTm9OyYxWFmxc65wpbLdWe7iEgH1yM3gUNHpvL825Udcr4SJRIRkThwypGZVFQ388qCqliH8g1KJCIicWD/vZMY3D+JJ+ZW0NzcsboklEhEROKAmXHyERls2NZI8Ycd6wZFJRIRkThx2EgfWekenp9fGetQ/ocSiYhInEhKNI4Zl8b8pTXsLGuKdThfUSIREYkjx49Pp6kZXnq349RKlEhEROJI3+6JjNwvmefnV3aYTnclEhGRODP10HQ27WjqMJ3uSiQiInHm0BH+Tvfn3u4YzVtKJCIicearTvdlNewsj32nuxKJiEgcmnJwOs3N8EZxdaxDUSIREYlH/Xsmsk+fRF5bFPshU5RIRETi1KQxaaz6rJ4NWxtiGocSiYhInJpU6MMM5hbFtnlLiUREJE7ld0tgxL7JvLqoiljOLaVEIiISx44ck8a6LY18vC52zVtKJCIicWziKB+JCcS0012JREQkjmX4PBy0fypzi6ppitGQKUokIiJxbtKYNHaUNfHBmrqYnF+JREQkzo3dP4XkROOt92Nz9ZYSiYhInEtN9jBmaArzltTEZERgJRIRkU5gwoE+dpQ1seqz+qifW4lERKQTGDcslQQvMWneUiIREekE0lM9jB6cwrwl1VG/OVGJRESkk5hwoI/NO5qifnOiEomISCdxyPBUPB6YF+XmLSUSEZFOIivdy4iBybwV5eYtJRIRkU5kwkgf67Y08tmm6DVvKZGIiHQih4xIBWDB8tqonVOJRESkE8nPTmDfvoks+KAmaudUIhER6WQOPiCVFWvrKKtsisr5lEhERDqZg4el0uxg4YroNG8pkYiIdDL79UuiW6aHBcuj07ylRCIi0sl4PMa4/VNZuLKGxqbIXwasRCIi0gkdPCyVqhoXlTlKlEhERDqh0YNTSEwgKs1bUU0kZjbFzFab2Rozu2YX683Mbg2sX2Zmo4LW/dDMlpvZCjO7Yhf7/tjMnJnlRbgYIiIdXmqKh5H7pfDuBzXU1Daz/JM6nnmzgp1l4b+SKyHsR2yFmXmB24GjgfXAIjOb7ZxbGbTZscDAwGMs8E9grJkdAFwEHATUAy+Z2fPOuY8Dx+4bOO4X0SqPiEhHN+6AVP7+WAlTf7SeL0dMycnycthIX1jPE7VEgj8JrHHOrQUws0eAaUBwIpkG3O/8g8QsMLNsM+sJDAEWOOeqA/u+CZwE/Cmw31+Bq4FnolISEZE4MGmMj7Ub6inISWDfPkns2yeRvGxv2M8TzUTSG1gX9Ho9/lrH7rbpDSwHbjCzXKAGOA4oAjCzE4ENzrmlZtbqyc1sJjAToF+/fiEVREQkHmSmefnRmbkRP080E8muPuVbXpe2y22cc6vM7EZgDlAJLAUazcwH/AKYvLuTO+dmAbMACgsLoz+psYhIJxXNzvb1QN+g132AjW3dxjl3t3NulHNuArAT+BjYBxgALDWzzwLbLzazHhEpgYiIfEM0E8kiYKCZDTCzJOB0YHaLbWYD5wSu3hoHlDnnNgGYWUHgZz/gO8DDzrkPnHMFzrn+zrn++BPRKOfc5iiVSUSky4ta05ZzrtHMLgNeBrzAPc65FWZ2cWD9HcAL+Ps/1gDVwPlBh3gy0EfSAFzqnCuJVuwiItI6i/Yk8R1BYWGhKyoqinUYIiJxxcyKnXOFLZfrznYREQmJEomIiIREiURERELSJftIzGwb8Pke7p4HbA9jOPFAZe4aVOauIZQy7+Wcy2+5sEsmklCYWdGuOps6M5W5a1CZu4ZIlFlNWyIiEhIlEhERCYkSSfvNinUAMaAydw0qc9cQ9jKrj0REREKiGomIiIREiUREREKiRNKKUOaXj1dtKPOZgbIuM7N3zGxELOIMp92VOWi7MWbWZGanRDO+cGtLec3scDNbYmYrArORxrU2/F9nmdmzZrY0UObzd3WceGJm95jZVjNb3sr68H5+Oef0aPHAPzrxJ8DeQBL+ibSGttjmOOBF/JNxjQPei3XcUSjzIUC3wPNju0KZg7abi3906lNiHXeE/8bZ+Ke/7hd4XRDruKNQ5p8DNwae5+Of7ygp1rGHWO4JwChgeSvrw/r5pRrJrn01v7xzrh74cn75YF/NL++cWwB8Ob98vNptmZ1z77ivh+9fgH8isXjWlr8zwA+AJ4Gt0QwuAtpS3jOAp5xzXwA457pCmR2QYf65utPxJ5LG6IYZXs65t/CXozVh/fxSItm11uaOb+828aS95fke/m808Wy3ZTaz3sBJwB1RjCtS2vI33g/oZmZvmFmxmZ0Ttegioy1lvg0Ygn821g+AHzrnmqMTXsyE9fMrmnO2x5M9nl8+ArFES5vLY2ZH4E8kh0Y0oshrS5n/BvzUOdfk/8Ia19pS3gRgNDAJSAXeNbMFzrmPIh1chLSlzMcAS4Aj8U/fPcfM5jnnyiMcWyyF9fNLiWTXQppfPk61qTxmNhy4CzjWObcjSrFFSlvKXAg8EkgiecBxZtbonHs6KhGGV1v/r7c756qAKjN7CxgBxGsiaUuZzwf+6PydB2vM7FNgMLAwOiHGRFg/v9S0tWshzS8fp3ZbZjPrBzwFnB3H31CD7bbMzrkBzrn+zrn+wBPAJXGaRKBt/9fPAIeZWYKZ+YCxwKooxxlObSnzF/hrYJhZd2AQsDaqUUZfWD+/VCPZBRf6/PJxp41lvhbIBf4R+Ibe6OJ45NQ2lrnTaEt5nXOrzOwlYBnQDNzlnNvlJaTxoI1/498B95rZB/ibfH7qnIvroeXN7GHgcCDPzNYDvwYSITKfXxoiRUREQqKmLRERCYkSiYiIhESJREREQqJEIiIiIVEiERGRkCiRiIhISJRIREQkJEokImFmZtlmdknQ63cidJ4+ZnZaJI4t0h5KJCLhlw18lUicc4dE6DyT8M85IRJTSiQi4fdHYJ/ALIM3mVklgJn1N7MPzewuM1tuZg+a2VFmNt/MPjazg748gJmdZWYLA8f4l5l5g09gZocCfwFOCWwzIKolFAmiIVJEwszM+gPPOecOCLyudM6lB5avAQ4EVuAfUHAp/iH5TwTOd85NN7MhwJ+A7zjnGszsH8AC59z9Lc7zEvDjeB4LSzoHDdooEl2fOuc+ADCzFcBrzjkXGDCwf2CbSfjnBFkUGBwzlV3PzjgIWB3xiEV2Q4lEJLrqgp43B71u5uv3owH3Oed+1tpBzCwX/9DfDRGJUqQd1EciEn4VQEYI+7+Gv++jAMDMcsxsrxbbDCC+J1KTTkSJRCTMAjNHzg90qN+0B/uvBH4JvGJmy4A5QM8Wm32If66J5WYWqavCRNpEne0iIhIS1UhERCQkSiQiIhISJRIREQmJEomIiIREiUREREKiRCIiIiFRIhERkZD8P+2gSXLgS5h1AAAAAElFTkSuQmCC\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",
    "energies = np.real(result.observables)[:, 0]\n",
    "\n",
    "plt.plot(result.times, energies, color=\"royalblue\")\n",
    "plt.xlabel(\"time $t$\")\n",
    "plt.ylabel(\"energy $E$\")\n",
    "plt.title(\"Energy over time\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we also kept track of the total magnetization of the system, we can plot that quantity too. However let's first compute exact reference values to verify our algorithm results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy as sc\n",
    "\n",
    "def exact(final_time, timestep, hamiltonian, initial_state):\n",
    "    \"\"\"Get the exact values for energy and the observable.\"\"\"\n",
    "    O = observable.to_matrix()\n",
    "    H = hamiltonian.to_matrix()\n",
    "    \n",
    "    energ, magn = [], []  # list of energies and magnetizations evaluated at timesteps timestep\n",
    "    times = []  # list of timepoints at which energy/obs are evaluated\n",
    "    time = 0\n",
    "    while time <= final_time:\n",
    "        # get exact state at time t\n",
    "        exact_state = initial_state.evolve(sc.linalg.expm(-1j * time * H))\n",
    "        # store observables and time\n",
    "        times.append(time)\n",
    "        energ.append(exact_state.expectation_value(H))\n",
    "        magn.append(exact_state.expectation_value(observable))\n",
    "\n",
    "        # next timestep\n",
    "        time += timestep\n",
    "\n",
    "    return times, energ, magn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.quantum_info import Statevector\n",
    "\n",
    "initial_state = Statevector(ansatz.bind_parameters(initial_parameters))\n",
    "exact_times, _, exact_magnetizations = exact(final_time, 0.01, hamiltonian, initial_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jul/opt/miniconda3/envs/py10/lib/python3.10/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return np.asarray(x, float)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Magnetization over time')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEYCAYAAABLOxEiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAw4klEQVR4nO3dd3gU9drG8e+T0HvvVUXEDkawYsOKHEXxWOAIWBBBRcReEAs2LIiiiGJBFGwgKHAQRUQEVECanSpIlxJ6y/P+seM5++YEyJLdzCa5P9c1V3bKztyTwD47v5n5jbk7IiIi2ZUSdgAREclbVDhERCQmKhwiIhITFQ4REYmJCoeIiMREhUNERGKiwiESxczGmln7BKz3RzM7Pd7rTSZm1tbMPgs7hySe6T4OSRQzWwzUAGq4+9qo6bOAY4D67r44lHCRHL2AQ9y9XZzX+yawzN3vj+d6k4mZ1QMWAYXdfXfIcSSX6YhDEm0RcOXfI2Z2FFA8vDgSKzMrFHYGSS4qHJJobwNXR423BwZHL2BmLc3sBzNLN7OlwZFA9PyrzWyJmf1lZg+Y2WIzaxHM62Vm75vZYDPbFDQJpUW9t4aZfWRma8xskZndEkw/D7gXuNzMNpvZ7GD6RDO7Lng9O5j39+B/NzeZ2QdmttLMNprZJDM7IpjeCWgL3Bm855NgenTmombW18yWB0NfMysazDvdzJaZWQ8zW21mK8ys495+ucH+jTKzdWY238yuj5q+zcwqRC3b2MzWmlnhYPwaM/vZzNab2Tgzqxu1rJtZVzP7Hfg9i01PCn5uCPbzRDPrYGaTM62ji5n9HvxtHjGzg81savC3ft/MikQtf6GZzTKzDWY2xcyO3tt+S8jcXYOGhAzAYqAF8CvQCEgFlgJ1AQfqBcudDhxF5IvM0cAq4OJg3uHAZuAUoAjwNLALaBHM7wVsBy4I1v84MC2YlwLMAHoG7z0IWAicG/XeIZkyTwSuy2JfOgG/AGWC8WuA0kBRoC8wK2rZN4FHs/pdBK8fBqYBVYDKwBTgkajfxe5gmcLBfm0Fyu/ld/wV8BJQDDgWWAOcFcybAFwftWwfYEDw+mJgfvB3KQTcD0yJWtaB8UAFoHgW260XLFMoaloHYHKmdYwCygBHADuAL4K/Q1ngJ6B9sGwTYDXQLPg7tg9+Z0XD/nes4X8HHXFIbvj7qONsIh++f0bPdPeJ7j7X3TPcfQ4wFDgtmN0G+MTdJ7v7TiJFIPOJucnuPsbd9wTbOiaYfjxQ2d0fdved7r4QeBW4IpbwZnYK8CjwD3dPDzK/7u6b3H0HkQJ0jJmVzeYq2wIPu/tqd18DPAT8K2r+rmD+LncfQ6RwNswiV20iBfUud9/u7rOA16LW9S5BM6GZGZH9fjeYdwPwuLv/7JFzFI8Bx0YfdQTz17n7tmzuV1aedPd0d/8RmAd85u4L3X0jMBZoHCx3PfCKu3/r7nvc/S0iheaEHGxbEkSFQ3LD28BVRL6RDs4808yamdmXQXPSRqAzUCmYXYPIUQoA7r4V+CvTKlZGvd4KFAva5esCNYKmjw1mtoFI81TV7AYPPpzfJ/LN+LdgWqqZPWFmC8wsncg3Y6Iy708NYEnU+JJg2t/+8v9/wnkrUGov61nn7psyratm8PpD4EQzqwE0J1Jwvw7m1QWej/q9rAMs6r0Q9XvPgVVRr7dlMf73ftUFemT6W9Xm//9eJEnopJcknLsvMbNFRJpdrs1ikXeBF4Hz3X27mfXlvx/CK4j6tm1mxYGK2dz0UmCRuzfYW7R9vTnY1sdAX3cfGzXrKuAiIs1wi4k0u6wn8sG73/UCy4l8UP4YjNcJpsVqOVDBzEpHFY86BEd07r7BIpfH/pNIk9RQd/8721Kgt7u/s4/172s/4n055t95esd5vZIAOuKQ3HItcKa7b8liXmki35y3m1lTIh/Mf/sQaGVmJwUnUh/ivx/Q+/MdkG5md5lZ8eBI4UgzOz6YvwqoZ2Z7+3/wOvCLuz+VRd4dRI58ShBp5om2ikg7/t4MBe43s8pmVolI89uQbO7Tf7j7UiLnRx43s2LByeRrgehi8C6RZsJL+W8zFcAA4J6ok/plzeyyGDa/Bshg3/sZi1eBzsHRp5lZSYtcNFE6TuuXOFLhkFzh7gvcffpeZncBHjazTUQ+RN+Pet+PwM3AMCJHH5uInETdkY1t7gFaETlpvAhYS+QcwN/nIj4Ifv5lZjOzWMUVQOtMV1adSqS5bQmRb/Y/ETnRHW0QcHjQ5PJxFut9FJgOzAHmAjODaQfiSiInqpcDI4AH3X181PxRQANglbvP/nuiu48AngSGBc1t84Dzs7vRoMmwN/BNsJ85OhcR/Nu4nsiR53oiJ+475GSdkji6AVDyFDMrBWwAGrj7opDjiBRIOuKQpGdmrcyshJmVJHI57lz+e0JaRHKZCofkBRcRaYpZTqTZ5QrXobJIaNRUJSIiMdERh4iIxKRA3MdRqVIlr1evXtgxRETyjBkzZqx198pZzSsQhaNevXpMn763K0FFRCQzM1uyt3lqqhIRkZiocIiISExUOEREJCYqHCIiEhMVDhERiUnSFA4zqx08k+Fnizz+s1sWy5iZ9QsekTnHzJqEkVVEpCBLpstxdwM93H1m0JXyDDMb7+4/RS1zPpEuJxoQecTky8FPERHJJUlTONx9BZFus3H3TWb2M5GnkUUXjouAwUE/RdPMrJyZVQ/eG3c/9uzLW1Mm0Oa4E2lQpQYrtm1i5OzvadXsZOrUqs36PTuZtWwxxx+XRqU6NdlTohgp5UpTtFTJRMQREUkKSVM4oplZPSLPIv4206ya/P/HWS4Lpv1P4TCzTkAngDp16hxQjt8GvkufVd/T8PsFVCxSgdm7NnLvxrnU+mwmxYuU45udG2ifPo+hZY/i+MJl+Wrneq5N/5GPapxAWs16zC68k1eX/0yvCy6jfqOGrCtZhGXs5LgzT6NknRpEHgMtIpK3JF3hCJ638BFwq7unZ56dxVuy7KXR3QcCAwHS0tIOqCfH1iu/Y09GBmRkYBlO7W3bWbNuHSVSC1M4A0qvXMUXc+bQqGotSu1xjv31V+6dNIHDGhxL0W27WT/vB35auYwtIyew7v0JvLd9Jfdtns/E8mnULl2eL0vvZvim5Tzf7gaqHtOILdUrUKrRIZSqVOFA4oqI5Iqk6h3XzAoDnwLj3P3ZLOa/Akx096HB+K/A6ftrqkpLS/OwuxzZk76ZpT/MZcakbzijSm0y/ljF0C8/44VZXzOq9DEUznCe37KEl7ct48fj2lDmyAbMLm1sqlyGi67vSOHqWXYZIyKSEGY2w93TspyXLIXDIu02bxF59vSte1mmJXATcAGRk+L93L3p/tadDIVjX3znLnYtWc7kj0fz7TeT6VjxEHb+OJ+bZn3O9F3pfF3heFKrVeL9UtvZU60it97eg6KNG5FSoljY0UUkn8orheMU4GsiT3fLCCbfC9QBcPcBQXF5ETgP2Ap03MdzrP8j2QvH3mxcsYqFX02l/sZd7Jj1M+3fe5XNmzfzVtkjoVAqjxVeS/1GDenevTvFmh1NavkyYUcWkXwiTxSORMqrhSMr21etZc/s39j+/Twu7/cY9bZmcE+xumBGF1/Cxc3PoP3NXSl2wtGkFCsadlwRyaNUOPJR4cgsY/sOdsz8mdUTpnD5c4/Sak9prihcha1FU7ktYym3te1Aixs7UvigWmFHFZE8RIUjHxeOzDK2bGPb1Nn88P5I2r/dn15F6nJSkXIsq1WeD0tso9tdd3JwyzOx1NSwo4pIElPhKECFI7Odi/5k2+dTeXfQG3T7dixflE+jdrVqLGpSn01HH8wFt9xAITVpiUgmKhwFuHBE27RyNUydw9YxX3PrsEF8umUl3x/UggotT2dT82Ooe+GZFCqqIiIiKhwqHFnYsn4DP7z9IQ1+WcGWsZO5fNk0ChcpzKhb7qfUZedQtHEj3dkuUoCpcKhw7FPG9h0Mffgptk+bzRm/rWf39h1cs2s+17e+jKsevZ/CtaqGHVFEctm+CkfSdTkiuS+lWFHaPvYAAHs2bmLh4BHs7nU36cM/548xs9l18tF8f0hF/tnzDoqWLhVyWhEJm444ZK92LV7OpvfGMqj/S9y5dCYj6pzMKR2vpFTblhRrUC/seCKSQGqqUuHIkd07d/JF/0EcPXsZW8d9w9MbF/Bz6RQ+fuV1yl54GlZIB64i+Y2aqiRHChUpwrndbwRg96q/OKhLd7ZM+oa/ru/FhuqVmZpWhzPu7Eq1wxqEnFREckPSPDpW8oZCVStyx0eDeWvlb1R7+3G21a9Ghzf6cWfamay++TF2zJsfdkQRSTAVDjkglppKyfNOodHIl/hu9Dhua9eBzaO+ZOKpl3N61brMevN9CkIzqEhBpMIhOXbMBWdz7IBHqDt7OBsvP4tF69eys8ezLG3engWvDGXn1m1hRxSROFLhkLhJLVeafw54ioXp62j0ci/MjOtv7kpapRqsH/A+GVu3hx1RROJAV1VJwrg7w3s/w7KPxnHhsu1YhTKMbVKDq/o8RPlaNcKOJyL7sK+rqnTEIQljZlx6/+10+2E8NT7pz4J6lbhp2Gu80Phs1vV5gz0bN4UdUUQOgI44JFdNefdDaoz5jj3jpvJVoW38flQt7h/yGiWrVAo7mohE0RGHJI2TrmpDvSFPUWvC68ysVJj3JnzG8hPbsv7ZwezZtCXseCKSDSocEoqiRzXgxblT+e7ryZQ64VhWPzaQZlVrM+jabmRs3xF2PBHZBxUOCVWlE5tQ/Z0nKfZub4qWKsnuD8bzR7OrWPfWSHbvUAERSUYqHJIU6p5zGpNX/UG7TwdTqHol+na9jSMqVOWP90frRkKRJKPCIUnDzCjZPI2aYwdw1G3X0aRkJXZ2eZzlrW5i+YQpYccTkYAKhyQdM+Oyh+/h3T9/ofLTt7P29wUc2eI0HjzhHHYtWR52PJECT4VDkpYVLkTZDhdT/6vBdDjtHJou3sgfJ7Vj8X3PsXHFqrDjiRRYKhyS9MpUq8KzX47m3JmfUOris3j42ac5rE49lr/2AZ6REXY8kQJHhUPyjEI1qlC1/320H/A01x18LNvu6cef597AwjETwo4mUqAkVeEws9fNbLWZzdvL/NPNbKOZzQqGnrmdUcJ3WvsrefjnKVR55UF++2MxjVq24LkzW7N75dqwo4kUCElVOIA3gfP2s8zX7n5sMDycC5kkCZkZpS9pwXETh3DTaedx0k9r+OPEtix6ZhC7tqkXXpFESqrC4e6TgHVh55C8o2z1qjwzcQyNvxlKsaZHcf29d3JitTps+eaHsKOJ5FtJVTiy6UQzm21mY83siL0tZGadzGy6mU1fs2ZNbuaTEBQ5uDbVh/Wh4+23clnZ2qy8+BZW3/I46X/8GXY0kXwnrxWOmUBddz8GeAH4eG8LuvtAd09z97TKlSvnVj4JkZnRtvcD3PHzZMrd3JbPhgyj/kEHMan387r7XCSO8lThcPd0d98cvB4DFDYz9cct/09KyeJU7NmZI996iuZV61LxufdYfvEtbPtlYdjRRPKFPFU4zKyamVnwuimR/H+Fm0qS1bGtzmX40l+o3fdetv74Oycd3ZgnL26H79gZdjSRPC2pCoeZDQWmAg3NbJmZXWtmnc2sc7BIG2Cemc0G+gFXuNogZB8sJYUy7S6k4mevUKd2bUpOmMHSs65l+/dZXvEtItmgJwBKgbJl/FTW3vE07y+cx4KG1Xju3x9TvGL5sGOJJB09AVAkUPLsE6k9+W3+PKI2M2fNYuXZndj6lb5UiMRChUMKnJRSJeg79XPGT5hAapHC/HbJzdx8/OlsXL4y7GgieYIKhxRYZU89jlpfvsGss49iwPRJTDjtSrZ++V3YsUSSngqHFGgpxYvSYdhAfvr0M46rXJMV/+zBsEuuJX3l6rCjiSQtFQ4RoEHLFtSaMIidHVrSfsQbdDv6ZLZNnhl2LJGkpMIhEkgpVpRGfe5mzCtvcGuto1jeuhu/dX+cres3hh1NJKmocIhkclan9hw1+V1KX9OaG/r34aQ6B7Pth5/DjiWSNFQ4RLKQUqIYVZ68jVsf7cW/ytVl+QU3sv6Zt8jYtSvsaCKhy1bhMLPyZna1mY0wsx/N7FMzu97MqiY6oEiYLrrzFm6dO5FS/ziDUQ/34cSqdVn6nbpsl4Jtv4XDzIYDw4EqwF3ufgTQBSgGvG1mExOaUCRkqeVKU/WVBynU6VJ2bd3Kpja3k/7uaPW4KwXWfrscMbNy7r7hQOcnA3U5IvGyc+lK1tz8GJsmz2Bg7VTu+/gdKtatHXYskbjLUZcj+ysKyV40ROKpSO1q1PjoOX6/4nRemDWZj07/J9umzAo7lkiuivnkuJmdbWavmtmxwXinuKcSSWKWmkqrF3sz95NxnFe5Hstbd2Ny94fYtV3POpeC4UCuquoC3AG0M7MzgWPjmkgkjzjswrOpNWEQOy48mfOff4TOjZqya9mqsGOJJNyBFI417r7B3W8HzgGOj3MmkTwjpVQJDh/0GM/dfDvtd5Vj2ekdSP9kYtixRBLqQArH6L9fuPvdwOD4xRHJm657/glO/HoohevX4vrLrqBTWnP2bFPTleRP2bkcd62Ztfx73N1HRs939xcSEUwkrylcvybVP3mRCsceTqGfFrH8gi7sXLA07FgicZedI45dwMtmdl3mGcGjXkUkkFqsKC9Pm8BTI4ax+89VfHHqP3mzx/1hxxKJq+wUjhVAc+B2M+uVaV7DuCcSyQdKnXsytb98nTdS1nPnc0+x6JZHydi+I+xYInGRrXMc7r4YOBk4x8xeMzP1cSWyH4VqVmXIbz/wyU33kTF0HMvO68wf02aEHUskx7JTAAzA3f8CziTS9cgoMyvx9zwRyVqREsVp1u9Bqr3zJIN/+p4jTjqBGS/rehLJ27JTOP7To5u7bwcuBv4EJgJlEpJKJJ8pec5JXDHyDdrXO4IyDwxk7QMv4Lt2hx1L5IBkp8uRazKNZ7j7DUQuy62fqGAi+c0hJx7PC798T7nrLmXxS+/Qum4jlsycE3YskZgd8LkKd38IFQ6RmFiRwlR+ojururXhq1VL+OHSrnpEreQ5OTrJ7e5LzOyueIURKSha9rydBdN/oEn1uiy/tDvjezxERkZG2LFEsuVAOjl8P2r4APif+ztEZP8qND6CWp8NZMlJh3Hes714pNm5ZGzaEnYskf06kCOOdHf/ZzBcBnwerzBm9rqZrTazeXuZb2bWz8zmm9kcM2sSr22LhCGlVAlO/+hl+rbtROsl21l2Tid2/roo7Fgi+xRzlyNA70yL3BfHPG8C5+1j/vlAg2DoBLwcx22LhCIlJYWbh7xCgxH92LlhIxc0bso79zwcdiyRvYq5yxF3/8/XITMb6u7r4hXG3ScB+1rfRcBgj5gGlDOz6vHavkiYip/cmNIfPM36IimsePEd/npkAL5nT9ixRP5HXutypCYQ3WvcsmDa/zCzTmY23cymr1mzJlfCieRU1SMP49vli2l/w/Vs6PcOn57TjjULF4cdS+T/yWtdjmR1p3qWD01394HunubuaZUrV05wLJH4KVSiOJWfvp1ST9zCNROH0+64U9jx04KwY4n8R17rcmQZUDtqvBawPJcziOSKqtdexgevvs79lRvx5/k3smnUl2FHEgHyXpcjo4Crg6urTgA2uvuKXM4gkmtOv6YtJ04aSpEjDqbr5W3pcUZLMnTeQ0KWVF2OBM/3mAo0NLNlZnatmXU2s87BImOAhcB84FUizz8XydcKVatE9Y+eI+OQWqR/O5tV7e/T/R4SKnPP8hTBfxcwS3H3LG9pNbO6wNK9zU8WaWlpPn369LBjiORIRkYGGwcNZ90DL7KsRhkq9bmdI85qHnYsyafMbIa7p2U1r1A23j/ezNYCHwOj3T09OL9xHtAaOBJoHK+wIpK1lJQUyl/fhqKNDuKKc1qw9vwvmD1mPGVanBR2NClgstNUdRbwEFAPGG1m04AvgKOB59xdRUMkF5U4pQnvfDaG5444jTVt72HDKx+onyvJVfttqvqfN5gVC06S5xlqqpL8KGPzVlZ3fZRBwz/gp7rleH3aRIqWKhl2LMkn9tVUFfP9GHmtaIjkVymlSlD1jUfZ3Oxw/vxtASsvv4M9a9eHHUsKAD07XCQPs5QUHv98JB+/O5Q9c37lpxbX8MOof4cdS/I5FQ6RfKBcm3OpMfIFHl4ygzNat+LP4SoekjgqHCL5RLEmh/P0F5/S54jT2N75MTa8PIxYz2GKZEe2T46bWVHgUiJXV/3nMl53T/r+n3VyXAqSjK3bWd31Ub4cMYqxtYvz6tQJOmkuMYvXyfGRRLo13w1siRpEJImklChG1UEP8+uJhzLl57ksuPw29qxPDzuW5COxHHHMc/cjE5wnIXTEIQXVqiEj2XzX86TUrMzup7vRoLluFpTsidcRxxQzOypOmUQkF1RtdxE1hvel78IfSDvjNH4foZPmknOxFI5TgBlm9mvwvO+5ZjYnUcFEJD6KNzuariMG06XWkaR0fZL0oWPCjiR5XCxNVXWzmu7uS+KaKAHUVCUCezZuYtU1DzD/y2/4d5MaPPLvj0gtlJ3u6qQgiktTVVAgygGtgqFcXigaIhKRWrY01Yc9zbjDK/HihNHMbNuDjG07wo4leVC2C4eZdQPeIfIEwCrAEDO7OVHBRCT+rHAhHpk0mm96PUuFCTNZfkk30pf8GXYsyWNiOcdxLdDM3Xu6e0/gBOD6xMQSkURJSUnhqJ63UPX1R3jv20kceWhDfp/4TdixJA+JpXAYEP3Myj3k/jPHRSROSrU6nWZ9e9K4WDkyrn2Ybd/8sP83iRBb4XgD+NbMeplZL2AaMCghqUQkV5xy9RV8OOc7ilerzPxLuzHszl5hR5I8IJaT488C1wDrgPVAR3fvm6BcIpJLCtetQc3RL/NW+Z207fMQ39/1hPq4kn2K6Vo8d58BzEhQFhEJSWq50vT+/gtOaHszFV4fzZptTuVn7sAK63Jd+V/7PeIws8nBz01mlh41bDIzdYAjkk8ULVWSSz8eRPnbO/D14Pc4p86hrFuqK67kf2XnmeOnBD9Lu3uZqKG0u5dJfEQRyS1mRoW7rmXjv85h0ZpVLPlnd3YvXx12LEkysdzH8WR2polI3nfNC08xY+x4yq7cyNJzb+CXcV+GHUmSSCxXVZ2dxbTz4xVERJJL2bNPouYn/Xll7e+knX82894ZEXYkSRL7PfNlZjcCXYCDMnVqWBqYkqhgIhK+okcewo2fvoO36UixHn3ZVKg4pS8/L+xYErL9dnJoZmWB8sDjwN1Rsza5+7oEZosbdXIokjN70jezqsN9LJ44lQmnHsK9I4aQkqInT+dnOerk0N03uvtid78SSAeqAnWBI82seZyDnhd02z7fzO7OYv7pZrbRzGYFQ894bl9EspZaphTVhz3Nx4eU4olP3ueHTvfhu3eHHUtCku2LtM3sOqAbUAuYRaSvqqnAmfEIYmapQH8i51KWAd+b2Sh3/ynTol+7+4Xx2KaIZJ8VKcwT0z6nTY9HKTdkPCu33Uellx+gcJlSYUeTXBbLsWY34HhgibufATQG1sQxS1NgvrsvdPedwDAizzgXkSSRkpJC0+d6Uump2/jg01E0rX0QK3+ZH3YsyWWxFI7t7r4dwMyKuvsvQMM4ZqkJLI0aXxZMy+xEM5ttZmPN7Ii9rczMOpnZdDObvmZNPOubiJTt2Joa3TtQZmcG6e3uZdci3ShYkMRSOJaZWTngY2C8mY0ElscxS1Y97WY+cz8TqOvuxwAvBFmy5O4D3T3N3dMqV64cv5QiAkCbXnfx76++pPCmrSw473q+e+/jsCNJLomlk8PW7r7B3XsBDxDpGTeeTUnLgNpR47XIVJjcPd3dNwevxwCFzaxSHDOISAyKNz2KmqNf4qn18znjyjYs/EDPMy8IYrlz/AszuwDA3b9y91FETmbHy/dAAzOrb2ZFgCuAUZkyVDMzC143JZL/rzhmEJEYFTmkDo9MGMkTjU7Gb36K9KEqHvldLE1V9YG7zOzBqGlZXuN7INx9N3ATMA74GXjf3X80s85m1jlYrA0wz8xmA/2AK1z9P4uEruaRjeg6bQzFT2nMxC738eCF/yQjIyPsWJIg+70B8D8Lms0kcuVTPyJNSu2AL929SeLixYduABTJHb5zFzcefxqj5s1kYtf7afDcPVhqatix5ADk6AbA6PW4+2537wJ8BEwGqsQjoIjkD1akMC/N/Jrx3R8kdeg4Vl7bk+3pm8KOJXEWS+EY8PcLd38T6AB8Fuc8IpLHpaSmcsTT91Dx4a489cEQmtdryMY/V4QdS+IolquqXsk0PsPdr4l/JBHJD8rdeAWNb7yag7bDxivvYvfKtWFHkjjREwBFJGGufvoRXh87it1/rGR2iw7MHvt52JEkDvQEQBFJqBKnpVFz5AvctWQ65/2jFRu+mRl2JMkhPQFQRBKu6DENeXnMx/Sp34x1V97JlvFTw44kOaAnAIpIrmhw6glcMXk4hRvUZchl1/PiDbeGHUkOUHbOcdxoZnOBhmY2J2pYBMxNfEQRyS8KValAzY/7MbLkdl57fRBr+g5G9/DmPdl5Hse7wFjy8BMARSR5pJQuyYhfZrG468Ok936V3WvWU+WhrqQWyvbjgSRk2X4CINAWOBVo7+5LgFJBf1EiIjEpWroUh775BGWuvYTOz/SmTaPG7Nm+I+xYkk2xnOPoD5wIXBmMbyK+nRyKSAFiKSlUevxWjjvnTA5avolVV99LxpZtYceSbIilcDRz967AdgB3Xw8USUgqESkQzIyeo4bx4IAX2PbVdCad254Vv/wedizZj1gaFXcFzwV3ADOrDKj7SxHJsTL/aoWXK8nZl7SiStOvmPTjbIrUrhZ2LNmLWI44+gEjgKpm1ptIJ4ePJSSViBQ4ZVudyaAX+9OzTAOWt+zCzl8XhR1J9iKWvqreAe4kUiyWAxe7+weJCiYiBc8ZN3TgrHFvwp49PH7KhXw+8K2wI0kWYrlzvCjQBCgLVAQuM7OeiQomIgVT0SMOoeLw5xiyeRn9u9+tu8yTUCxNVSOJPGN8N7AlahARiavSDQ9i8g8z6NP0XFZefQ/rho4OO5JEieXkeC13Py9hSUREolQ97BAqj+rPwqvu4Jyrr+KiD1vxwEdDwo4lxHbEMcXMjkpYEhGRTFJKl6TWO09Sr04dyo+fwV+PDFAXJUkglsJxCjDDzH4N+qqaa2ZzEhVMRASgWJnSfPDbLC7vdC0b+r3DxLa3snObbhQMUyxNVeoJV0RCYampVOrTg7VFjVaP30XbWVN5ecbXpBQvGna0AinbhSPon0pEJBRmRqPePXhmxRKOHvktKy7vQbUhT5BaplTY0QqcbBcOM7sti8kbgRnuPituiURE9uGG1/ux6aPxrOr6KHcf05ybRgym7rFHhh2rQInlHEca0BmoGQydgNOBV83szvhHExHJWulLz2Z7n1sYsGQOAy5sy67Fy8OOVKDEUjgqAk3cvYe79yBSSCoDzYEOCcgmIrJXR/7rUqaPHMP1xWvyZ8sb2Tbnt7AjFRixFI46wM6o8V1AXXffBqgjfRHJdQ1bnUOt0S+x2ndx3PFpjH3h1bAjFQixFI53gWlm9qCZPQh8Aww1s5LAT/EIY2bnBZf7zjezu7OYb2bWL5g/x8yaxGO7IpJ3FTm0HlUHP07xIkXZ/uhAtoz7JuxI+V4snRw+AlwPbCByUryzuz/s7lvcvW1OgwRdtvcnctnv4cCVZnZ4psXOBxoEQyfg5ZxuV0Tyvrppx/Dtot84/pjGrGx/H98/PSDsSPlaLEccAAuBqcBMoISZNY9jlqbAfHdf6O47gWFE+saKdhEw2COmAeXMrHocM4hIHlWoUnlqDO/L7MOqcMIdNzKo481hR8q3Yukd9zpgEjAOeCj42SuOWWoCS6PGlwXTYl0GADPrZGbTzWz6mjVr4hhTRJJVSqkSXPjpG9yVdibNPpnBXw+/rC5KEiCWI45uwPHAEnc/A2gMxPMT2bKYlvkvnp1lIhPdB7p7mrunVa5cOcfhRCRvKFKqJL2nfUbVa9qwqt8Qep5yPju2bg07Vr4SS+HY7u7bIfJsDnf/BWgYxyzLgNpR47WIPDAq1mVEpICz1FQqPdmdmS2b0HvKOIb/oyMZ23TxZ7zEUjiWmVk54GNgvJmNJL4f2t8DDcysvpkVAa4ARmVaZhRwdXB11QnARndfEccMIpJPmBnt3nqBKQ89R9M5K1hxeQ92b9wUdqx8IZarqlq7+wZ37wU8ALzG/568PmDuvhu4ici5k5+B9939RzPrbGadg8XGEDlBPx94FegSr+2LSP50Qs9bqfJKT2ZN/ZYT6hzMopnq1DunYumrKg24D6gbvM+A3sDR8Qrj7mOIFIfoaQOiXjvQNV7bE5GCoXTrFhT5YwHpd/RgRYd7qTXyFQrXz/K6GsmGWJqq3gHeAC4FWgEXBj9FRJLe2d1u4IfJU6i+3Vl2wY38Nm5i2JHyrFgKxxp3H+Xui9x9yd9DwpKJiMRZyaZHU/PT/gxJ/4Njz2/BzCEfhR0pT4rlQU4PmtlrwBdE9U3l7sPjnkpEJEGKHFqPjp++zV//aEuZu/qzpUJlSl4Qz3uZ879Yjjg6AscC5xFpovq7uUpEJE+p1fgonpn9NcWOasBP7e9mwA3dw46Up8RyxHGMux+VsCQiIrkotUJZanzUlz7Hn8GAgc9zUrW6HNWrG2ZZ3Wcs0WI54piWRaeDIiJ5VkrJ4jw9/Us+vaQTJV/6iL8eeAHPyAg7VtKLpXCcAswKuj2fY2ZzzUwXRItInla4RHFafPASZTtdxscvvMJFDY5mW7puFNyXWJqqzktYChGREFlKChUfvZkNv85g8cjh/NH+HhoM6UNKyeJhR0tKsdw5viSrIZHhRERyi5lx+4dv8eXAwRSaMpclrW9h5W8Lwo6VlGJ9HoeISL5WsWNrqr3xCA9MGUuzo49l3S8qHpmpcIiIZFLyguZ06vsEbUvWJP3yO9n5y6KwIyUVFQ4RkSw0v64dD076FHbvZtK57fn3S6+FHSlpqHCIiOxF0SMOoeaYl3kifQHX3HIT6z75MuxISUGFQ0RkHwrXrcF7079h8PEtWX/9Q6S/OzrsSKFT4RAR2Y+KB9fjzPFvU/zUJjzd6RZuO/sfZBTgGwVVOEREsiGlVAmqDXmCJXXK8+OkKay59/kCe5d5LDcAiogUaClFi/DG3G9Zef/zbH5tOOkrVlKj7z2ULF8u7Gi5SoVDRCQGKampVH+sO+trVKVlj67YmHf4esEvpJYpFXa0XKOmKhGRGJkZFW5uS9ebbuLyXWVYfvEt7F71V9ixco0Kh4jIAWr37KNc9+Hr7FqwlA+bt2Hm2PFhR8oVKhwiIjlQssUJVP3wWR5YPJ0bLr2CbTN/CjtSwqlwiIjkUMnjj2Ls55/x4sEnsaL1rWz94tuwIyWUCoeISBwccuqJHPf5mxSqX5Mb/nEJT3fsEnakhFHhEBGJk0JVK1J1+HOsLVuc398bxfq+b+PuYceKOxUOEZE4KlqhHJ8smEevf13Hut4DmdulFzu2bQs7VlypcIiIxFmh4sWo9nJPit3QhosGPsXlhx9HxrYdYceKm6QoHGZWwczGm9nvwc/ye1lucfCs81lmNj23c4qIZJelpFDz0W7c0e4a2m4owopLb2XPuo1hx4qLpCgcwN3AF+7eAPgiGN+bM9z9WHdPy51oIiIHrstb/Tn/zefYPvtX+je7gLkTJoUdKceSpXBcBLwVvH4LuDi8KCIi8VXqojMoO7g3jy2awQOtr2LHnN/CjpQjyVI4qrr7CoDgZ5W9LOfAZ2Y2w8w67WuFZtbJzKab2fQ1a9bEOa6ISGwqnXUik8aN57E6x/HnP25i8+dTw450wHKtcJjZ52Y2L4vhohhWc7K7NwHOB7qaWfO9LejuA909zd3TKleunOP8IiI5dehZp3Lo+EFQpzoXtGzJU+07hx3pgORa4XD3Fu5+ZBbDSGCVmVUHCH6u3ss6lgc/VwMjgKa5lV9EJB4KVatE1Q+foXSlivjwCazr80aeu9cjWZqqRgHtg9ftgZGZFzCzkmZW+u/XwDnAvFxLKCISJyWrVOKThT/S4er2rH/qdb5q152t6ZvCjpVtyVI4ngDONrPfgbODccyshpmNCZapCkw2s9nAd8Bod/93KGlFRHIopWgRqrx4L6ldL+OSof3peMyJZGzeGnasbLG8doh0INLS0nz6dN32ISLJacgt91Br6AQOPvIIqg99ikLVKoUdCTObsbfbHpLliENEpMBq1+9xmg7rx65Fy7jj2OZMGvZR2JH2SYVDRCQJlDirGSXfeZwR6//g7RtuZ+tXydtKosIhIpIkqp58HN/NmM59R5/CiituZ+lr7yflFVcqHCIiSaTKkYdRe/TL7G52BGd07kD3M1omXfFQ4RARSTKpZUpxyHvPcf6xTTn+h6WsvvERfMfOsGP9hwqHiEgSSi1ahH7ff8kFD93J5o/G88YpF/PHvJ/DjgWocIiIJC0zo/yt/6JY39vpMWM83U+/gJ0LloYdS4VDRCTZ1Wx7EeOHDKNnuUP58/zObJ0yK9Q8KhwiInlA2lWXcvhng6B8aS4662z6XNs1tCwqHCIieUThg2pRecTzpJYrw9b3xrHuyUGhXHGlwiEikoeUqlGVMYt/5sZrrmX902/y5WVdWL8qyw7FE0aFQ0Qkj0kpWoTKz99N8Ts7ctWI17ni6GbsXr0u97afa1sSEZG4MTNq3HENgx56nDtSqvHnuZ3Y8dOCXNm2CoeISB7W8v7bOPXfb+K793DHSecwovczCd+mCoeISB5X9JiGVB3dn5XFU6n13IdseHlYQk+aq3CIiOQDxevUYNj82dRudRZ/9ezPmlufxHfuSsi2VDhERPKJlJLFqTroYcrf1p5N745meZvuCXmqYKG4r1FEREJjKSlUuOc6Cjesx7avpmMli8d9GyocIiL5UOlLWlD6khYJWbeaqkREJCYqHCIiEhMVDhERiYkKh4iIxESFQ0REYqLCISIiMVHhEBGRmKhwiIhITCyMp0flNjNbAyw5wLdXAtbGMU5eoH3O/wra/oL2OVZ13b1yVjMKROHICTOb7u5pYefITdrn/K+g7S9on+NJTVUiIhITFQ4REYmJCsf+DQw7QAi0z/lfQdtf0D7Hjc5xiIhITHTEISIiMVHhEBGRmKhwAGZ2npn9ambzzezuLOabmfUL5s8xsyZh5IynbOxz22Bf55jZFDM7Joyc8bS/fY5a7ngz22NmbXIzXyJkZ5/N7HQzm2VmP5rZV7mdMd6y8W+7rJl9Ymazg33uGEbOeDGz181stZnN28v8+H9+uXuBHoBUYAFwEFAEmA0cnmmZC4CxgAEnAN+GnTsX9vkkoHzw+vyCsM9Ry00AxgBtws6dC3/ncsBPQJ1gvErYuXNhn+8FngxeVwbWAUXCzp6DfW4ONAHm7WV+3D+/dMQBTYH57r7Q3XcCw4CLMi1zETDYI6YB5cysem4HjaP97rO7T3H39cHoNKBWLmeMt+z8nQFuBj4CVudmuATJzj5fBQx39z8A3D2v73d29tmB0mZmQCkihWN37saMH3efRGQf9ibun18qHFATWBo1viyYFusyeUms+3MtkW8sedl+99nMagKtgQG5mCuRsvN3PhQob2YTzWyGmV2da+kSIzv7/CLQCFgOzAW6uXtG7sQLRdw/vwrlKE7+YFlMy3yNcnaWyUuyvT9mdgaRwnFKQhMlXnb2uS9wl7vviXwZzfOys8+FgOOAs4DiwFQzm+buvyU6XIJkZ5/PBWYBZwIHA+PN7Gt3T09wtrDE/fNLhSNSfWtHjdci8k0k1mXykmztj5kdDbwGnO/uf+VStkTJzj6nAcOColEJuMDMdrv7x7mSMP6y+297rbtvAbaY2STgGCCvFo7s7HNH4AmPnACYb2aLgMOA73InYq6L++eXmqrge6CBmdU3syLAFcCoTMuMAq4Ork44Adjo7ityO2gc7XefzawOMBz4Vx7+9hltv/vs7vXdvZ671wM+BLrk4aIB2fu3PRI41cwKmVkJoBnwcy7njKfs7PMfRI6wMLOqQENgYa6mzF1x//wq8Ecc7r7bzG4CxhG5IuN1d//RzDoH8wcQucLmAmA+sJXIN5Y8K5v73BOoCLwUfAPf7Xm4Z9Fs7nO+kp19dvefzezfwBwgA3jN3bO8rDMvyObf+RHgTTObS6QZ5y53z7PdrZvZUOB0oJKZLQMeBApD4j6/1OWIiIjERE1VIiISExUOERGJiQqHiIjERIVDRERiosIhIiIxUeEQEZGYqHCIiEhMVDhEcsjMyplZl6jxKQnaTi0zuzwR6xaJhQqHSM6VA/5TONz9pARt5ywiz10QCZUKh0jOPQEcHDxFr4+ZbQYws3pm9ouZvWZm88zsHTNrYWbfmNnvZtb07xWYWTsz+y5Yxytmlhq9ATM7BXgWaBMsUz9X91AkirocEckhM6sHfOruRwbjm929VDB9PtAY+JFIB3yziXRT/w+go7tfbGaNgKeAS9x9l5m9BExz98GZtvNv4Pa83JeU5A8FvpNDkQRb5O5zAczsR+ALd/egg716wTJnEXkmxvdBh5LFyfoJhA2BXxOeWGQ/VDhEEmtH1OuMqPEM/vv/z4C33P2eva3EzCoS6Q57V0JSisRA5zhEcm4TUDoH7/+CyLmLKgBmVsHM6mZapj55++Fhko+ocIjkUPB0xG+CE+B9DuD9PwH3A5+Z2RxgPFA902K/EHnewjwzS9RVWyLZopPjIiISEx1xiIhITFQ4REQkJiocIiISExUOERGJiQqHiIjERIVDRERiosIhIiIx+T8QGLOA7/LhXgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "magnetizations = np.real(result.observables)[:, 1]\n",
    "\n",
    "plt.plot(result.times, magnetizations, color=\"crimson\")\n",
    "plt.plot(exact_times, exact_magnetizations, \":\", color=\"k\")\n",
    "plt.xlabel(\"time $t$\")\n",
    "plt.ylabel(r\"magnetization $\\langle Z_1 Z_2 \\rangle$\")\n",
    "plt.title(\"Magnetization over time\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks pretty good!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient-based optimizations\n",
    "\n",
    "The `PVQD` class also implements parameter-shift gradients for the loss function and we can use a gradient descent optimization routine\n",
    "\n",
    "$$\n",
    "    \\theta_{k+1} = \\theta_{k} - \\eta_k \\nabla\\ell(\\theta_k).\n",
    "$$\n",
    "\n",
    "Here we're using a learning rate of \n",
    "\n",
    "$$\n",
    "    \\eta_k = 0.1 k^{-0.602}\n",
    "$$\n",
    "and 100 optimization steps in each timestep."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from qiskit.algorithms.optimizers import GradientDescent\n",
    "\n",
    "maxiter = 100\n",
    "learning_rate = 0.1 * np.arange(1, maxiter + 1) ** (-0.602)\n",
    "gd = GradientDescent(maxiter, lambda: iter(learning_rate))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "pvqd.optimizer = gd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Warning:** The following cell will run a few minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_gd = pvqd.evolve(problem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Energy over time')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEYCAYAAABPzsEfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA86UlEQVR4nO3dd3zTdf7A8dc7adp0b4ZAAZEtUDY4UQRRObfnOkXv1ENFTz3Peec6uXOd6yfquReuQxE8OQeKoIiyRPaeZXXRPdIkn98f31BrLZCWJmma9/PxyIPku/L+tiXvfLYYY1BKKaUOxRbqAJRSSoUHTRhKKaX8oglDKaWUXzRhKKWU8osmDKWUUn7RhKGUUsovmjCUiiAi8ryI/C3UcajwJDoOQ4ULEdkKtAU8dTa/ZoyZFJqIWjYRuQK4yhhzXKhjUa1DVKgDUKqRfmOMmR3INxCRKGOMO5Dv0dzCMWYVfrRKSrUKInKFiHwrIo+JyD4R2SIip9XZnywiL4vIbhHZKSIPioi9zrnzReQJESkE7hORdBH5WERKRGSR7/hvfcdPEZF/1Xv/j0XkpgPEdozvGsW+f4/xbb9IRBbXO/ZmEZnpex7ju5/tIrLXV50U69s3SkRyROR2EdkDvFrvOr2B54GRIlImIkW+7a+JyIP1rnGbiOT6fjZni8jpIrJeRApF5K4617SJyB0isklECkTkfRFJa/xvS4UrTRiqNRkOrAMygEeAl0VEfPteB9zAUcBAYCxwVb1zNwNtgMnAFKAcaAdM8D2oc62LRcQGICIZwGjgnfoB+T5QPwGeBtKBx4FPRCQdmAn0FJHudU65BHjb9/xhoAeQ7Yu7A3BPnWPbAWlAZ+Cauu9rjFkDTAQWGGMSjDEpv/5x1V7DWefaLwK/AwYDxwP3iMiRvmNvBM4GTgSOAPb5fk4qUhhj9KGPsHgAW4EyoKjO42rfviuAjXWOjQMM1gdiW6AaiK2z/2JgTp1zt9fZZwdqgJ51tj0IfFvn9RpgjO/5JGDWAWK+DFhYb9sC4Arf87eAe3zPuwOlvtgFK2F1q3PeSGCL7/kowAU4D/LzuqJuzL5trwEP1rlGJWD3vU70/cyG1zl+CXB2nXseXWdfe9/PKSrUfxv6CM5D2zBUuDnbHLgNY8/+J8aYCl/hIgHrW7gD2P1zgQMbsKPOuXWfZ2K17x1oP1iljN8BX/j+feoAMR0BbKu3bRvWN3qwShP/Ah7AKl185Iu9DVbiWFInZsFKZvvlGWOqDvC+/iowxuzvRFDp+3dvnf2VWD9DsEoy00XEW2e/Bysh7zzMOFQY0IShIsEOrBJGhjlww3Dd7oJ5WNVXHYH1vm2d6h3/FrBSRAYAvYGPDnDdXVgftHVlAZ/6nn8OZIhINlap52bf9nysD+u+xpgDfRgfqotjc3eB3AH83hgzv5mvq8KEtmGoVs8Ysxvrg/lfIpLka7ztJiInHuB4D/AhVuN3nIj0Ai6vd0wOsAh4E/jAGFP56ysBMAvoISKXiEiUiFwI9AH+67uOG5gGPIpVEvrCt92L1Z7whK+0gYh0EJFTG3Hre4GOIhLdiHMO5nlgsoh09sWTKSJnNdO1VRjQhKHCzce+Xj/7H9P9PO9yIBpYjdVYOw2rDv5AJgHJWNVcb2I1aFfXO+Z1oJ9vf4OMMQXAeODPQAFwGzDeGJNf57C3gVOA/9QrAd0ObAS+F5ESYDbQ86B3+UtfAauAPSKSf6iD/fAUVkP95yJSCnyP1VlARQgduKeUH0TkYaCdMWZCnW0nYFVNdfGVCJRq1bSEoVQDRKSXiPQXyzDgD8D0OvsdwJ+AlzRZqEihCUOphiVitWOUA+9j9WSaAbWD4oqwqrSeDE14SgWfVkkppZTyi5YwlFJK+aVVj8PIyMgwXbp0CXUYSikVNpYsWZJvjMlsaF+rThhdunRh8eLFhz5QKaUUACJSf2aCWlolpZRSyi+aMJRSSvlFE4ZSSim/tOo2jIbU1NSQk5NDVdXhTvKpWiqn00nHjh1xOByhDkWpViXiEkZOTg6JiYl06dKFOtNGq1bCGENBQQE5OTl07do11OEo1apEXJVUVVUV6enpmixaKREhPT1dS5BKBUDEJQxAk0Urp79fpQIjIhOGUkq1Vq6awE33pAkjBOx2O9nZ2QwYMIBBgwbx3XffAbB161ZiY2PJzs6ufbhcLgA+/fRThg0bRq9evcjOzubCCy9k+/btAHz//fcMHz6c7OxsevfuzX333ReqW1NKhVCVy8tVk3cz7auSgFw/qI3eIjIOaxEWO9a00A/V298LeBUYBNxtjHnMt90JzANifDFPM8bcG8zYm1NsbCzLli0D4LPPPuPOO+9k7ty5AHTr1q12334rV67khhtuYObMmfTu3RuAmTNnsnXrVrKyspgwYQLvv/8+AwYMwOPxsG7dumDejlKqhTAGhvZx0r1jcy2y+EtBSxgiYgemAGOAHGCRiMw0xqyuc1ghcCNwdr3Tq4GTjTFlvnUIvhWR/xljvg9C6AFVUlJCamrqQY95+OGHueuuu2qTBcCZZ55Z+zw3N5f27a3F4+x2O3369AlMsEqpFi02xsYNv00L2PWDWcIYBmw0xmwGEJF3gbOwlswEwBiTC+SKyBl1TzTWHOxlvpcO3+OwK+qe+c8+NuW4Dvcyv9CtYzSTLjh4AqisrCQ7O5uqqip2797NV199Vbtv06ZNZGdnA3DssccyZcoUVq1axa233nrA691888307NmTUaNGMW7cOCZMmIDT6WyW+1FKtXxuj2HyqwWcf3IifY+MCdj7BLMNowOwo87rHN82v4iIXUSWAbnAF8aYH5o3vODZXyW1du1aPv30Uy6//HL2r0uyv0pq2bJlTJky5VfnFhQUkJ2dTY8ePXjssccAuOeee1i8eDFjx47l7bffZty4cUG9H6VUaL33RQlzl1ZQVOoJ6PsEs4TRUF9Hv0sJxhgPkC0iKcB0ETnaGLPyV28icg1wDUBWVtZBr3mokkAwjBw5kvz8fPLy8g54TN++fVm6dCkDBgwgPT2dZcuW8dhjj1FWVlZ7TLdu3bj22mu5+uqryczMpKCggPT09GDcglIqhLbvreGNWcWcOCiOYwfEBfS9glnCyAE61XndEdjV2IsYY4qAr4EGv0YbY14wxgwxxgzJzGxwSvcWZe3atXg8noN+uN92221MnjyZNWvW1G6rqKioff7JJ5/UllA2bNiA3W4nJSUlYDErpVoGr9fwr7cKcUbbuOG3gf8CHMwSxiKgu4h0BXYCFwGX+HOiiGQCNcaYIhGJBU4BHg5YpAG2vw0DrKksXn/9dex2+wGP79evH0899RSXX345paWlpKenk5WVxf333w/Am2++yc0330xcXBxRUVFMnTr1oNdTSrUOXy+tYMWmam67LI20pMD/nw/qmt4icjrwJFa32leMMZNFZCKAMeZ5EWkHLAaSAC9WQ3cfoAvwuu88G/C+MeaBQ73fkCFDTP0FlNasWfOL3kaqddLfs4oEHo9h7o8VnDQ4rtlmOBCRJcaYIQ3tC+o4DGPMLGBWvW3P13m+B6uqqr7lwMDARqeUUuHBGEN5lSEh1sbJQ+KD9r460lsppcLMnCUVXH7vLrbvqQnq+2rCUEqpMFJc5uGZ9/fRNj2KDm2Cu0KFJgyllAojz07bR2mFl1svTcNuC+7MzJowlFIqTCxcVckXCyu4+NQkugVovqiD0YShlFJh4vuVlWS1jeJ345JD8v6aMEJg7969XHLJJRx55JEMHjyYkSNHMn369MO65n333feLqUJmz57dpOssW7aMWbNmHfpAYNSoUdTvthxMH330EatXrz70gUq1Ejf8NpWnb21LtCM0i4RpwggyYwxnn302J5xwAps3b2bJkiW8++675OTk/OpYt9vdpPd44IEHOOWUU5p0bmMSRqhpwlCRYuMOFzm5NYgISfGhG5SrCSPIvvrqK6Kjo5k4cWLtts6dO3PDDTcA8Nprr3HBBRfwm9/8hrFjx1JWVsbo0aMZNGgQ/fr1Y8aMGbXnTZ48mZ49e3LKKaf8Yg2MK664gmnTpgGwZMkSTjzxRAYPHsypp57K7t27Aat0cPvttzNs2DB69OjBN998g8vl4p577uG9994jOzub99577xexV1ZWctFFF9G/f38uvPBCKisra/d9/vnnjBw5kkGDBnHBBRfUznN1xx130KdPH/r371874+7evXs555xzGDBgAAMGDKhdQOqtt95i2LBhZGdn88c//hGPx5pILSEhgbvvvpsBAwYwYsQI9u7dy3fffcfMmTP5y1/+QnZ2Nps2bWqeX5BSLYyrxjD51Xzufi4Przd4A60bEtw+WS3QzU/s/dW2UYPiOOvERKpcXu6c8utJAU8dEc+4kQkUl3m478X8X+x74ua2B32/VatWMWjQoIMes2DBApYvX05aWhput5vp06eTlJREfn4+I0aM4Mwzz2Tp0qW8++67/Pjjj7jdbgYNGsTgwYN/cZ2amhpuuOEGZsyYQWZmJu+99x533303r7zyCmCVYBYuXMisWbO4//77mT17Ng888ACLFy/mmWee+VVczz33HHFxcSxfvpzly5fX3kd+fj4PPvggs2fPJj4+nocffpjHH3+cSZMmMX36dNauXYuIUFRUBMCNN97IiSeeyPTp0/F4PJSVlbFmzRree+895s+fj8Ph4LrrrmPq1KlcfvnllJeXM2LECCZPnsxtt93Giy++yF//+lfOPPNMxo8fz/nnn3/Qn6dS4Wzqp8Vs2+PmoeszsQW5V1R9EZ8wQu3666/n22+/JTo6mkWLFgEwZswY0tKsRVCMMdx1113MmzcPm83Gzp072bt3L9988w3nnHMOcXHW7JR1F1Tab926daxcuZIxY8YA4PF4ahdaAjj33HMBGDx4MFu3bj1krPPmzePGG28EoH///vTv3x+wlohdvXo1xx57LAAul4uRI0eSlJSE0+nkqquu4owzzmD8+PGAVcp64403AGvBp+TkZN58802WLFnC0KFDAas006ZNGwCio6Nrzx08eDBffPGFXz9bpcLdphwXb39WwinD4hjWNzbU4WjCOFiJwBltO+j+5AT7IUsU9fXt25cPPvig9vWUKVPIz89nyJCfp26Jj/95qP/UqVPJy8tjyZIlOBwOunTpQlVVFcAh544xxtC3b18WLFjQ4P6YGGuhFbvd7nd7SUPvaYxhzJgxvPPOO7/at3DhQr788kveffddnnnmmV8sFlX/GhMmTOCf//znr/Y5HI7a921MrEqFM4/H8NhbhSTG2bj+/NAvxQDahhF0J598MlVVVTz33HO12+pOVV5fcXExbdq0weFwMGfOHLZt2wbACSecwPTp06msrKS0tJSPP/74V+f27NmTvLy82oRRU1PDqlWrDhpfYmIipaWlDe474YQTmDp1KmCtM758+XIARowYwfz589m4cWPt/axfv56ysjKKi4s5/fTTefLJJ2vXKh89enTt/Xs8HkpKShg9ejTTpk0jNzcXgMLCwtp7bUqsSoU7jxf6dovhxgtTSU5oGbNPa8IIMhHho48+Yu7cuXTt2pVhw4YxYcIEHn644dnaL730UhYvXsyQIUOYOnUqvXr1AmDQoEFceOGFZGdnc95553H88cf/6tzo6GimTZvG7bffzoABA8jOzq5tYD6Qk046idWrVzfY6H3ttddSVlZG//79eeSRRxg2bBgAmZmZvPbaa1x88cX079+fESNGsHbtWkpLSxk/fjz9+/fnxBNP5IknngDgqaeeYs6cOfTr14/BgwezatUq+vTpw4MPPsjYsWPp378/Y8aMqW2gP5CLLrqIRx99lIEDB2qjt2p1oh3CpAtSGTU4eJMLHkpQpzcPNp3ePHLp71mFK6/X8Pg7hYwbkcDR3QK3PveBHGx6cy1hKKVUC/LJ/DJmzS9nW5BnovWHJgyllGohcgvd/Ht6EQN7xnD6MS2nKmq/iEwYrbkaTunvV4UnYwxPvFOI1wt/vjS92VbQa04RlzCcTicFBQX6odJKGWMoKCjA6XSGOhSlGmXBikp+WFXF789M5oiMljnioWVGFUAdO3YkJyeHvLxfj+BWrYPT6aRjx4ZW+lWq5Rp+dCx3XJ7G6GEtrypqv4hLGA6Hg65du4Y6DKWUqlVe6SU+1sbYEQmhDuWgIq5KSimlWpJ5P1bwu3t3sWWXK9ShHJImDKWUCpGScg9PvVdIZoqdTm0doQ7nkDRhKKVUiDz3QRElZV5uuyydKHvL6xVVnyYMpZQKgR9WVfLZ9+VcNDaJozoFf33uptCEoZRSIfDDyko6t3dw2WmhWZ+7KSKul5RSSrUEN/w2lbJKE7L1uZtCSxhKKRVEa7ZUs2OvtT53Ylx4fQSHV7RKKRXGKqu8/P2VfB54OT8sZ5vQhKGUUkHywowi9hZ6uPG3qS1yrqhD0YShlFJBsGx9FTPmlnHuSYn0Oyo85zrThKGUUgFWWe3l0TcLOCIzij+cGT69ourTXlJKKRVgIjCyfxzHZ8fijA7f7+maMJRSKsCc0TYmXZAa6jAOW/imOqWUauEqq73c8Uwua7ZUhzqUZqEJQymlAuTFj4pYuLoKlzv8utA2RBOGUkoFwLL1VXw0t4xzRyUwoHt49oqqTxOGUko1s8qqOr2izkoJdTjNRhOGUko1sxnflLGn0MNtl6URG9N6Pma1l5RSSjWzC0Yn0qtzNP3DdIDegQQ19YnIOBFZJyIbReSOBvb3EpEFIlItIrfW2d5JROaIyBoRWSUifwpm3Eop5Y/ySi/7Sj3YbUJ2j9aVLCCICUNE7MAU4DSgD3CxiPSpd1ghcCPwWL3tbuDPxpjewAjg+gbOVUqpkHrug31cPXk3FVXeUIcSEMEsYQwDNhpjNhtjXMC7wFl1DzDG5BpjFgE19bbvNsYs9T0vBdYAHYITtlJKHdr3KyqZ9V05p45IIM7Zetot6grmXXUAdtR5nUMTPvRFpAswEPjhAPuvEZHFIrI4Ly+vKXEqpVSjFJd5eGxqAUce4WDCGeE7V9ShBDNhNDSXb6NGs4hIAvABcJMxpqShY4wxLxhjhhhjhmRmZjYhTKWUapz/e38fxWVe7piQHlYr6DVWMBNGDtCpzuuOwC5/TxYRB1aymGqM+bCZY1NKqSapcRvcHsPlZyRzVKfoUIcTUMHsVrsI6C4iXYGdwEXAJf6cKNZKIy8Da4wxjwcuRKWUahxHlHDvVRmE4QJ6jRa0EoYxxg1MAj7DarR+3xizSkQmishEABFpJyI5wC3AX0UkR0SSgGOBy4CTRWSZ73F6sGJXSqn6jDG8+FFR7frcNlvrrYraL6gD94wxs4BZ9bY9X+f5Hqyqqvq+peE2EKWUColZ35XzzuclpCXZ6NTWEepwgqJ19v1SSqkA2pXv5tlp+xjYM4ZzRiWGOpyg0YShlFKN4PEaHn6jAJvAbZelR0RV1H6aMJRSqhE++baMFRurueG3qbRNi6zp+CLrbpVS6jCNHRFPlF0YMzw+1KEEnSYMpZTyg6vG4PEaYmNsnH5sQqjDCQmtklJKKT+8PLOIP/5zT6udWNAfmjCUUuoQflxXxbSvShncy9lqJxb0R+TeuVJK+aG0wstDrxfQITOKa85JCXU4IaUJQymlDsAYwxNvF1JY4uGvv89oVcutNkVk371SSh1EVbUhr8jNleOT6ZHVuicW9If2klJKqQOIddp48ua2OjGRj5YwlFKqHo/HmliwqNSD3S7YI2g098FowlBKqXre+rSEdz4v4cd1VaEOpUXRhKGUUnWs3FTNm7OKOWVYHCcNibzR3AejCUMppXzKK73847V82qbZ+dOFaaEOp8XRhKGUUj4vzSwid5+Hu67MID5WPx7r015SSinlM+H0ZAZ0d9L3yJhQh9IiaQpVSkW8whIPHo8hJdHOqEFxoQ6nxdKEoZSKaDVuw51Tcvnbv/NCHUqLpwlDKRXRXppRxIYdNZx+TGROWd4Yh0wYItI3GIEopVSwfb+ikv98WcpZJyZwXLZWRR2KPyWMN/c/EZGr6u4QEf0JK6XCUl6Rm4ffLODIDg6uPTc11OGEBX8SRt0x8dfV2/dNM8ailFJBU15pyEy1c88fMoh26NQf/vCnW62p87z+T1XbQJRSYalLewf/vqMdIpos/OXPB347EblCRAby64RhGjpBKaVaqmXrq3ji7UKqXV5NFo3kTwnjPmAIcCXQUURWAWt9j4zAhaaUUs2rsMTDg6/kEx9rwxu5S3M32SEThjHmhbqvRaQj0B/oB8wLUFxKKdWsvF7DP18roKzS8MgNGcRG8NrcTdXoqUGMMTlADjCr+cNRSqnAmPpZCUvWVnHLJWkc2UFXz2sKTbFKqVavpNzDf2aXMHpoHGccq1OWN5VOPqiUavWS4u1Mua0dGcl2beg+DH6XMERkkojo6BalVNjweAzzfqzAGEOntg5ttzhMjfnptQMWicj7IjJONE0rpVq41/5bzH0v5vPThupQh9Iq+J0wjDF/BboDLwNXABtE5B8i0i1AsSmlVJN9t7yCqZ+VcPqx8WT3cIY6nFahUeUzY4wB9vgebiAVmCYijwQgNqWUapJd+W4eer2A7p0c3PhbXWq1ufjd6C0iNwITgHzgJeAvxpgaEbEBG4DbAhOiUkr5z+s1PPBSPgD3XZ2p80Q1o8b0ksoAzjXGbKu70RjjFZHxzRuWUko1jc0mXDE+GRFon6EdQZuT3z9NY8w9B9m3pnnCUUqppiss8ZCWZGfE0bGhDqVVakyV1C0NbC4GlhhjljVbREop1QRrtlRzy5O53HlFOicM1KV6AqExjd5DgIlAB9/jGmAU8KKI+NV+4euOu05ENorIHQ3s7yUiC0SkWkRurbfvFRHJFZGVjYhZKRUB9pV6uPfFfFKTbGT3iAl1OK1WYxJGOjDIGPNnY8yfsRJIJnACVjfbgxIROzAFOA3oA1wsIn3qHVYI3Ag81sAlXgPGNSJepVQEcHsMf38pn5JyL/dfk0lSvD3UIbVajUkYWYCrzusaoLMxphLwZ1TMMGCjMWazMcYFvAucVfcAY0yuMWaR79rU2zcPK6EopVSt5z8sYtmGam65JI3unXRSwUBqTBeCt4HvRWSG7/VvgHdEJB5Y7cf5HYAddV7nAMMb8f5KKfULxhgyUuxcMDqRscN1UsFA8yth+KYBeQ1rSvPjsFbem2iMWew75FJ/LtPAtmZfsU9ErsFqXyErK6u5L6+UaiG8XoPNJlw0JinUoUQMv6qkfCO8PzLGLDHGPGWMebJOsvBXDtCpzuuOwK5GXuOQjDEvGGOGGGOGZGZmNvfllVItQGGxh2v+sYel66pCHUpEaUwbxvciMvQw3msR0F1EuopINHARMPMwrqeUikCuGsO9L+axM89NUpzOPhtMjWnDOAmYKCJbgXKsKiZjjOnvz8nGGLeITAI+A+zAK8aYVSIy0bf/eRFpBywGkgCviNwE9DHGlIjIO1jdeDNEJAe41xjzciPiV0qFOWMMT79fyKrNLu75QzpHaSN3UDUmYZx2uG9mjJlFvaVdjTHP13m+B6uqqqFzLz7c91dKhbeP5pYxa345l45LYtRgbeQOtsaU57YDxwMTfPNJGaBtQKJSSql6jDGs3+7imP6xXDk+OdThRKTGlDCeBbzAycADQCnwAXA47RpKKeUXEeG2y9KocVsTDKrga0wJY7gx5nqgCsAYsw/QCkSlVECVlHu469lccnJrEBGdrjyEGpMwanzTexgAEcnEKnEopVRAuD2G+1/KZ8naKopK9eMm1BqTMJ4GpgNtRWQy8C3wj4BEpZSKeMYYnnynkB/XVfPnS9I4uptOKhhqjVkPY6qILAFG+zadretgKKUC5b3Zpcz6zuoRNXZEQqjDUTRuPYwYYBCQ7DvvAhHBGPNAoIJTSkUmt8fw7bIKRg2K0x5RLUhjeknNwLdgEv7NTquUUk0SZRf+9ac2gPaIakkakzA6GmN0PQqlVMDsynfz6sdF3HRRGvGxOu1HS9OYhPGdiPQzxqwIWDRKqYhVUu7hzim5FJV6KSzxaMJogRqTMI4DrhSRzVhVUo2aS0oppQ7EVWP427/z2VPg5pEb2tCprSPUIakGBHUuKaWUqs/rNfzz9QJWbKzm7ivTGdDdGeqQ1AHoXFJKqZDKL/KwclM1E89NYfRQnVCwJdO5pJRSIdUmLYqX/9qORF3bosXTuaSUUiHxxQ/lPDttH16vISnejrUStGrJdC4ppVTQ/bCqkkfeLGBjjgu3J9TRKH81ZS6pNjqXlFKqqdZsqeb+F/Pp2sHB3/+YqbPPhpGmziUl6FxSSqlG2rq7hjufzSMt2c5D17fRsRZhpjGN3hhj1gJrAxSLUqqV253vJs4pPHJDG9KS7KEORzVSoxKGUko1hcdjsNuFkf1iGdzrCK2GClNaHlRKBVRJuYfrHtnD7IXlAJoswpgmDKVUwJRXerljSh5bd9eQqlVQYU8ThlIqIKpcXu5+Lo8N213ce1UGg3vplB/hThOGUqrZuT2Ge1/IZ+Wmau66Ip1j+seFOiTVDLTRWynV7Ow26NU5mlGD4zhpiM4P1VpowlBKNRu3x7C30E2HTAdX/iYl1OGoZqZVUkqpZuH2GB58JZ9Jj+6lpFzn+2iNNGEopQ7b/mQx78dKLjk1iaR47RHVGmmVlFIh5PUaqmsMNW6Dq8bg9oDHa/B6wRjfA7AJxEQLMQ4hNkaIdkiLmd21brK49rwULhidFOqQVIBowlAqSIwxLFtfzfSvS9m0s4ayCi/llV68pvHXirJDQqyNxHgb6Ul20pLtZKbYaZ8RRfuMKDq2cdA2zY7NFvik8v4XJcz7sZLrzk/h/JM1WbRmmjCUCrAql5c5iyv48OtSNuXUkJxgY3AvJ4nxNhLjbMTG2IhxWKUGux3sAjabIIL1ALwGql1WaaSy2kt5paGs0ktJuZfCYg9rt7n4ZpmbGvfP7xvjELLaRXFkh2h6do6md5dojuwQjSOqeZPIeScn0rGtgxMGatfZ1k4ThlIBsmWXi0++LeOzH8oprzR0PcLBrZemMXpoHDHRzd986PUaCoo97MpzsyPXzfY9NWzdXcPC1ZV89r01LUeMQ+h3VAyDejkZ2tvJkR0cTaraqnZ5eeXjYi47LZmEOJsmiwihCUOpZrSv1MNXi8r5/IdyNuyowREFJwyM4zfHJ9CvW0xA2x1sNiEzNYrM1CgG9Ph5uzGG3H0e1mx1sXJjFUvWVfPC9CJemA7t0+2cMDCOEwbF0atztF/xVVZ5+eu/81i2vpp+3WI4LluTRaQQY5pQgRomhgwZYhYvXhzqMFQrt2WXi++WV7JgRSVrtrowBnpkRTN2eDwnD4kjJbHl9RjKL3KzaHUVc3+sYOnaKtweyGobxbhjEhg7PP6AU4+XVni5c0oua7e5uP2ydMYM10F5rY2ILDHGDGlwnyYMpRovJ7eGOUsqmLO4gq27awDo2Tmakf1iOT47lq5HhM9y96UVXr75sYL/LShj1WYXdhucNDiOC05Jonunn+9jX6mH2/4vl+17avjr7zM4XksWrdLBEoZWSSnlB4/HsH67i+9WVDL/p8raJNHvqBj+dGEqx2XHkZ7c8koS/kiMs3H6sQmcfmwC2/fU8PG3ZcyaX8bsRRUM7BnDhDOS6X+Uk5oag9tteHBiJkP7xIY6bBUCWsJQqgFer2FjTg1L11axbEMVKzdVU1FlsNmg/1ExHNs/luOz42iT1jq/c5VVePnv/DI++KqUgmIPw/rEcPXZqXRp78BubxnjP1RgaJWUUodgjGFHrpuf1lexbH01S9dVUVzmBaBzuyj6d3cyoHsMg3s5SU4Iz5JEUyxbX8mdU/LxGlPbZddmsyYXtImvG7BNiI8VkuLtJCfYSEuy0ybVTmZqFBkp1viQjBQ7SfG2FjPYUB2YVkk1k4oqL6UV3to+87ExLWe0rWqc/CI3q7e42LDdxYYc6999pVaCSE+2M6yPk8G9Yxncyxm2VU2H64dVldz/Uj5pyXbuuzqDZeurKKvw4vWCx1ilMK/XGum9f0xIcZmXLbtqKCz2/GpAYmKcjS5HODjyCAfdO0XTu2s0We0c2IMwuFA1j6AmDBEZBzwF2IGXjDEP1dvfC3gVGATcbYx5zN9zA8ntMUz/upTXPymmourn/wXJCTaO7hZD/6NiGHF0LJ3aOoIVkmqEsgovG3NcbNjhYv12F6s2V7OnwJocz2aDLu0cDO0Ty9HdYsjuEUOHzKiI/yLwvwVl/GtqIUce4eCh69uQlmz/RQP4oXg81piQ/GIPefs85BW52bHHzeZdLr5YWM6MeWUAxDmFPl1jGNTTycCeMRzVKVoTSAsWtCopEbED64ExQA6wCLjYGLO6zjFtgM7A2cC+/QnDn3MbcrhVUsVlHlZsrOaVj4vZuruGYX2cHD8wjhq3ocpl2La7huUbq9mdb5XVh/d1cv7oJAb1DGx/e9UwY6wPqW173Gze6WLdNuuxM+/n4c8ZKXb6dI3m6G4x9O0aQ7eO0brGdD35RW4uu3c3R3eL4b6rM4iPbd5Bhl6vISfPzdot1aze4uKnjdVs83UiSE6wMeLoWEb2i2VobyexTp0fNdhaSpXUMGCjMWazL6h3gbOA2g99Y0wukCsiZzT23OZS4zbc9WwuKzZV47L+holzCmOHx3HJqUlktfv1t6zcQjf/W1DOzHml/OXpXOKdQmK8jYQ4GykJdtKS7KQl2UhPttMmLYq2aVG0S7frjJ6NYIyhsMTL7nw3efvclJRbVSD7Sq1vsHv3udmT76as8ucvQJkpdnp0jubUEfH0yIrmqE7RBxxfoKwPcptNyEiJ4omb29CtY/NPIwLWAMOstg6y2joYO8LaVljs4cf1VfywspLvllsj02McwjH9Yzl5aBzD+sQGJBbVOMFMGB2AHXVe5wDDm/tcEbkGuAYgKyur0UE6ooSyCi/RDiHeKdhsUFlt+PyHCo4bEEdWu2h+2lDF/72/j75HxjCsj5OBPZ1MOCOZi8cmMWdJORu2uyip8FJa7qW43Mv2vTXsK/H8Yp4fsOp0j8iMolObKLLaOejc3np0yIyKqGK512soKvWSV+Qmv8iqxsgvsh55+9zk7vOQu8+Dq+bXpeH4WKFtahSZqXb6dImp/Rl2ae/Q5NAI5ZVeHng5nxMGxnHGsQn06hIT1PdPS7Yzemg8o4fG4/EYVmyqZu7SCr5eWsGcJRUkxds4dUQ8ZxyXQJZW/YZMMBNGQ5+A/taH+X2uMeYF4AWwqqT8vP4vPHdH+19tq6jyEuXrTmgTq2F09sJyPv6mDLsNBvZ08pffpXHqiAROHdFgXBSXecnd52FvoZvd+W525bnZmedm+aZqZi+qqD3WGS10OcLBUR2j6ZFlPboe4Wjx37A8XkNZhZeKKkOVy0tltaHc1xhqNYh6KC6zGkaLy63nJeVeiko9uOutt2MT60OkTaqdbh2jGdnPmom1XXoUbdPsJMfbSYy3tfifSTjILXRz13N5bN1dw6hBoR+MZ7cL2T2cZPdwcv0FqSxeU8WnC8r4cE4p//mylIE9Y7jg5CSG9XUGZTZe9bNgJowcoFOd1x2BXUE4t1nE1alL7XeUk4cnOalxG1ZtrmbhqkqWrquunQJi0epKMlLsvxjtKyKkJNpJSbTTI+vX1VqVVVZJZPOuGjbvrGFTjouvl5Tz32+txkFHFL6eJTF07xRNcoKNhFgb8bE2YmOEOKf1b1QT+8h7vYaKalM75XZFlfWBX+H70C+v9FK2/98KL+VVvmOrft5WtzroQJLibSTF20hOsNEuPYoeWTZSEn3dMH3dLzNSokhJtEVUKStU1myt5m/P51HtMvzzupY3IC/KLow4OpYRR8dSWOzh0wVlzJhXxl3P5dG5XRS/PSWJMcPjm/x3rxonmI3eUVgN16OBnVgN15cYY1Y1cOx9QFmdRm+/z60rFOMwjDFc+cButu91M6S3kwvHNL0R3BjD7gIP67dVs3abizVbXazf5qK6gaqZ/ey2nxfasdsFh936xrb/7QXweK3SgNsDrhpTO232oYhAvFNIiLMSVbzT92+skBhnq33Ex9pwxthwRgsJsVaCSIy3kRRn00FfLUhekZvL791NapKNyddmhs10Jm6P4eslFbz/ZQkbd9TQPt3O705L1sTRTFrMwD0ROR14Eqtr7CvGmMkiMhHAGPO8iLQDFgNJgBcoA/oYY0oaOvdQ7xeqgXvFZR4++baM6XPLKCj2kN09honnpTZYsmgst8ewK8/9i2/7ldXWGgkV1QaX78O/usbg8Vgrue2v7tm/epvdZn1zs9sgOlpwOoSYaPlFEohzWqWWOKfUbouNEa0CaGX+t6CMkUfHtsgJEg/FGMMPK6t4/ZNi1m130SEziqvPTuH47FjtpXgYWkzCCLZQj/R21Rg+mV/GW58Wc9tl6Qzv27KK+yryVFZ5eWxqIWefmEC/o5yhDqdZGGNYsKKSl2ZY3d/7dYth4nkp9A5yw31roQkjxKpc1uhwEeHNWcW4PYaLT03CGYBFdJQ6kF35bu553mrcvuniNMYflxDqkJqVx2P434JyXv24iKIyL2ccm8BVZyVr9/VGainjMCLW/sRgjGFXvpvPvi9n9qIK/vK7NLJ7tI5veaplW7S6kgdfKcAYwz+vb3mN283BbhfGH5fAyUPieP2TYj6YU8q3yyqYeF4qY4bFaTVVM9ASRgj8uK6Kf71dyK48N2efmMDVZ6XoiFYVMCs2VnHTE7l0be/g/j9m0CEzMsYxbMpx8cQ7haze4mJ4Xye3XJpGZop+Rz4UrZJqgapcXl6aUczH35Ty7G3t6NYxPHqoqPBhjEFE8HoNH8wpZfxxCcTGRNYXE6/X8NHcMl78qIioKJh0QSpjh8draeMgNGG0YHn73GSmWt96Fq2uZHAvHYykDt+GHS6efKeQe6/KaLVrdjTGztwaHnmrkBUbqzllWBw3X5SmpfoDOFjC0J9YiO1PFqs2V3P7M3ncMSWPolLPIc5SqmHGGP77bRmTHt1DXpGHIt+aHpGuQxsHT9zUhivHJ/Plogque2QPW3a58HoNrflLc3PTEkYLYYzhk/nl/N/7haQmWusPBHs+HxXeyiu9PP5OIXMWVzC4l5O7r0wPy/EVgbZ0bRWTX82vXf8ErGnuox1CdJS11k1cjBDrtPkmErWT7JudICPFmkC0TaqdDpmOVjnTsVZJhZH1213c+0IehSUebro4jdNGtq6ujypwnp22jw/nlHLlb6yJMLVq88Dyi9zMXliBy23weg0eD7jcBpfbmvmgstpQWWVNh1NaYc1/Vlrhpe7HpQi0S7PTqZ21KFS3jtEc1TGaTm2jwvpnrwkjzBSXeZj8agGnDI1j7AhNGOrAvF5DSbmXlEQ75ZXWandHd9OSaSDsXxRq/wSiOblutu+tYfueGrbtrqmdUSE+Vuh7ZAxHHxlDdg8nvbtEh9WUOJowwtD+Hi4AS9ZW0SMrmsQ4bXJSPyss9vDQGwXsK/Xw7G3tdObeEKpxG7bvqWHDDhert7hYuamarb5FoeKd1uy7w4+O5dgBsaS28GpCHbgXhvYni7IKL/e9kEdqkp1/XpdJhzaR0YdeHdw3yyr419RCqlyG689PIaplfwa1eo4ooVvHaLp1jGbcSGtbcZmHZeurWbK2isVrKpm/vJIn34Gju8Vw0uA4Rg2OIzkhvH5xWsIIAys2VvG3f+cD8ODETK1yiGCV1V6eeX8f/1tQTvdODu66IoPO7fVLREtnjGHzzhq+WVbBvB8r2bq7hig7DOsby2kj4xnRL7bFTOevVVKtwM7cGu58No+9hW7unJDOqMHxoQ5JhUCN23D9o3sY3ieWy89I1mqoMLUpx8UXC8v5clEFBcUe2qTZOfO4BM44LiHkpQ5NGK1EcZmHe/6dT68u0Vx7Xmqow1FBUuXy8s5nJVxwShIJsTZq3EYTRSvh8RjmL6/ko7mlLFtfTYxDOO2YeH57ShLt0kPTYqAJoxVx1Rii7GCzCbmFbjJS7GHdhU8d3E8bqnjsrUJ25rm5+8p0Rg/VkmVrtWWXi/dnlzJ7YTkGOGVoPJednhT0ub80YbRCJeUe/vDgHvofFcPtl6e3ygFEkayiysuLM4qYMbeM9hlR3HppGgN76szGkSBvn5v/fFnKzG/K8HgM40bGc9lpyUGb4kUTRitkjOG92aW8ML2IgT1jeOCaTOJjtdtta/HwGwV8/kM5552UyJW/SY64SQMVFBR7mPppMf/9tgy7TbhobBIXjUkkJsDr6GjCaMU+/6GcR98soGsHBw9d34a0pPDqpqd+VljiweMxZKZGsafATUGxh75Hao+4SLenwM2/pxcxd2kFbdLsXHdeakCXodXJB1uxscPjmXxtJjl73bwysyjU4agm8HoNM+eVMuH+XTz13j4A2qVHabJQgPW3cO9VGTx+UxsSYm3c92I+f/t3PrmF7qDHoiWMVmLDDhcdMqOIc9p+MUpctWz7pyFfs9XFwJ4x/OmiNLLa6rgK1TCPxzDtq1Je+28xNhtcc3YKvzk+oVk7vmgJIwJ07xRNnNNGZZWX2/4vj6XrqkIdkjqEuUsruPahPewpcHPXFek8dmMbTRbqoOx24cIxSbzyt/b0PTKGp97bx+3P5JFXFJzShpYwWpmCYg9/eTqXnXk13HlFBqMGxYU6JFWH12soLPGQkRJFWYWXtz8r5uJTk3WeMNVoxhg+/qaM5z8sIsoO54xKRASqXNY4nT+cmdKk62qjd4QprfBy93N5rNpczQ2/TeXsExNDHZICVm6q5pn/7KPGbXjhznZhNYOparlycmt4+I0CVm12AeCMFjJT7bx+7xFNup5OPhhhEuNsPHpDJn9/pYCn39uH22M4/+SkUIcVsfYWunnxoyK+WlxBerKdP56Tgk0LFKqZdGzj4Ok/t8VVY4h2SEDbLzVhtFIx0TbuvzqDl2YUMbJfbKjDiVhrtlRz85O5AFw6LolLxibpWtKq2YkIMdGBL7FqwmjF7Hbhj+dac04ZY/jPl6WccWyCDvALMLfHsGNvDV2PiKZHVjRnn5jAOaMSaRukkbpKBYp+ckSIDTtqePGjIv70+F7y9gW//3YkMMYw78cKfv/33dz6VC6V1V7sdmHiuamaLFSroAkjQvTIiuaf12eyp8DNdY/sZf12V6hDalWWrq3i+kf3ct+L+dht8OdL03AGoYpAqWDShBFBhvSO5f9ubUuUHW56fC/zf6oIdUitwpqt1dz6dC6FxR7+clkaL93dnmP6x+ngSdXqaDk5wnQ9Ipopt7Xj/hfzdUK7w7Bhh4tNOS7GjUygV+do/vr7dI4bEKezBqtWTRNGBEpLsvPkLW1qvwHP+7GCob2d2nvHDxt2uHjjk2LmL68kPdnOyUPiiXYIJw/RdSpU66cJI0LtTxZ7Ctz8/eV8sto5uP+aDDq20akpGrJ9bw0vTC/iu+WVJMQKV4xP5txRiVqiUBFFv1JGuHbpUTw0qQ0FxR4mPrSHr5eUhzqkFsMYQ5XLaz33WiO1rxifzNsPduDy05NJ0Ok8VITRqUEUYI1G/vvL+aze4uK8kxO5/vzIXTPc6zUsWFHJO5+XkJFi576rMwFqR9Iq1Zrp1CDqkNqmRfHkLW15aUZRxI4ZcNUYvlhYzrQvS9i2x027dDtjhsXXThevyUJFusj8ZFANivINMtvv6yXl7Mxzc9GYpIiYKO/dL0p47b/FHNXRwd1XpjNqUFxE3LdS/tKEoQ7ox3XVfPxtGfN+rOCWS9Lo2bl1rQC3fruLD+eUcuKgOEb2i2X8sQkc3S2GgT1idAyFUg3QhKEO6OZL0hjYy8kz7xdy3SN7OeuEBH5/ZgoJYTwXVWW1l6+XVDDzmzLWbXMRGyP07hINQFqynbRkXRNdqQMJ6v98ERknIutEZKOI3NHAfhGRp337l4vIoDr7/iQiK0VklYjcFMy4I9moQXG8du8RnHVCAjPmlbFyY3WoQzostzyZy6NvFVLlMky6IJX3/tGBs3S9EKX8ErQShojYgSnAGCAHWCQiM40xq+scdhrQ3fcYDjwHDBeRo4GrgWGAC/hURD4xxmwIVvyRLCHWxo0XpnH2iYl0amv9yXw4p5Q4pzB6aDyOqJZZfZNb6OarJRUsWFHJoze0IdohXHZaEglxNvp102onpRormFVSw4CNxpjNACLyLnAWUDdhnAW8Yay+vt+LSIqItAd6A98bYyp8584FzgEeCWL8ES+rnTWoz+s1zFlSzqrNLl79uJjzRydyxrEJxLWAkeLFZR7mLKlg7tIKlm+sxhjo3SWavCI3HTIdHNNfl6xVqqmCmTA6ADvqvM7BKkUc6pgOwEpgsoikA5XA6YAOsAgRm014+s9tWbi6inc/L+G5D4p49b/F3HRhKmNHJAQ1Fq/XsDGnBmeMkNXWwe58N0+/t4/O7aK44oxkTh4aR4dMHb2uVHMIZsJoqPxff9Rgg8cYY9aIyMPAF0AZ8BPQ4KIOInINcA1AVlZW06NVByUiDO8by/C+sazZWs0n88vocoTVeLxyUzWzF5Yzsn8s2d1jiIlu3pLH5p0ulm+sZsXGapZtqGJfiZczj0/gpovT6Nk5mtfuaV9bGlJKNZ9gJowcoFOd1x2BXf4eY4x5GXgZQET+4Tv2V4wxLwAvgDXSuzkCVwfXu0sMvbv83OV2254aPv+hnJnflGG3wVEdo+ndNZqrz04hNsZGZbWXGIdgszXchmCMobLakF/sYVeem115bqIdwvjjrNLLXc/lkVvoISPFzqCeTob2iWVIbydgJTJNFkoFRjATxiKgu4h0BXYCFwGX1DtmJjDJ174xHCg2xuwGEJE2xphcEckCzgVGBi901RhnHJvAmGHx/Li+ihUbq1m9pZrvllcy6QJrUOBT7+7ji4XlxDmFxDgb0VFCcqKdp25pC8AdU/JYtLrqF9c8soOjNmHcflk67dKjaJdu14ZrpYIoaAnDGOMWkUnAZ4AdeMUYs0pEJvr2Pw/Mwmqf2AhUAFfWucQHvjaMGuB6Y8y+YMWuGi/a8XOVFVA7vQbAcQNiaZtmp6zSS3mlobrGkBT/c7XVmGHxDOzpJC3RxhGZDo7IiCI16ef9A3s6g3szSilAJx9USilVx8EmHwx9P0illFJhQROGUkopv2jCUEop5RdNGEoppfyiCUMppZRfNGEopZTyiyYMpZRSftGEoZRSyi+teuCeiOQB25p4egaQ34zhhAO959Yv0u4X9J4bq7MxJrOhHa06YRwOEVl8oNGOrZXec+sXafcLes/NSauklFJK+UUThlJKKb9owjiwF0IdQAjoPbd+kXa/oPfcbLQNQymllF+0hKGUUsovmjCUUkr5JaIThoiME5F1IrJRRO5oYL+IyNO+/ctFZFAo4mxOftzzpb57XS4i34nIgFDE2ZwOdc91jhsqIh4ROT+Y8QWCP/csIqNEZJmIrBKRucGOsbn58bedLCIfi8hPvnu+sqHrhAsReUVEckVk5QH2N//nlzEmIh9Yy8RuAo4EooGfgD71jjkd+B8gwAjgh1DHHYR7PgZI9T0/LRLuuc5xX2EtE3x+qOMOwu85BVgNZPletwl13EG457uAh33PM4FCIDrUsR/GPZ8ADAJWHmB/s39+RXIJYxiw0Riz2RjjAt4Fzqp3zFnAG8byPZAiIu2DHWgzOuQ9G2O+Mz+vl/490DHIMTY3f37PADcAHwC5wQwuQPy550uAD40x2wGMMeF+3/7cswESxVpcPgErYbiDG2bzMcbMw7qHA2n2z69IThgdgB11Xuf4tjX2mHDS2Pv5A9Y3lHB2yHsWkQ7AOcDzQYwrkPz5PfcAUkXkaxFZIiKXBy26wPDnnp8BegO7gBXAn4wx3uCEFxLN/vkVdVjhhDdpYFv9Psb+HBNO/L4fETkJK2EcF9CIAs+fe34SuN0Y47G+fIY9f+45ChgMjAZigQUi8r0xZn2ggwsQf+75VGAZcDLQDfhCRL4xxpQEOLZQafbPr0hOGDlApzqvO2J982jsMeHEr/sRkf7AS8BpxpiCIMUWKP7c8xDgXV+yyABOFxG3MeajoETY/Pz92843xpQD5SIyDxgAhGvC8OeerwQeMlYF/0YR2QL0AhYGJ8Sga/bPr0iukloEdBeRriISDVwEzKx3zEzgcl9vgxFAsTFmd7ADbUaHvGcRyQI+BC4L42+bdR3yno0xXY0xXYwxXYBpwHVhnCzAv7/tGcDxIhIlInHAcGBNkONsTv7c83asEhUi0hboCWwOapTB1eyfXxFbwjDGuEVkEvAZVg+LV4wxq0Rkom//81g9Zk4HNgIVWN9Qwpaf93wPkA486/vG7TZhPNOnn/fcqvhzz8aYNSLyKbAc8AIvGWMa7J4ZDvz8Pf8deE1EVmBV19xujAnbac9F5B1gFJAhIjnAvYADAvf5pVODKKWU8kskV0kppZRqBE0YSiml/KIJQymllF80YSillPKLJgyllFJ+0YShlFLKL5owlFJK+UUThlJNJCIpInJdndffBeh9OorIhYG4tlKNoQlDqaZLAWoThjHmmAC9z2isdQ+UCilNGEo13UNAN9+qdY+KSBmAiHQRkbUi8pKIrBSRqSJyiojMF5ENIjJs/wVE5HcistB3jX+LiL3uG4jIccDjwPm+Y7oG9Q6VqkOnBlGqiUSkC/BfY8zRvtdlxpgE3/aNwEBgFdbEeD9hTRd/JnClMeZsEekNPAKca4ypEZFnge+NMW/Ue59PgVvDea4n1TpE7OSDSgXYFmPMCgARWQV8aYwxvonvuviOGY21JsUi30SPsTS84l9PYF3AI1bqEDRhKBUY1XWee+u89vLz/zsBXjfG3Hmgi4hIOta01DUBiVKpRtA2DKWarhRIPIzzv8Rqm2gDICJpItK53jFdCe9Fu1QroglDqSbyrUY439ew/WgTzl8N/BX4XESWA18A7esdthZrvYOVIhKoXlhK+UUbvZVSSvlFSxhKKaX8oglDKaWUXzRhKKWU8osmDKWUUn7RhKGUUsovmjCUUkr5RROGUkopv/w//SONyqyBBAUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "energies_gd = np.real(result_gd.observables)[:, 0]\n",
    "\n",
    "plt.plot(result.times, energies, \"-\", color=\"royalblue\", label=\"BFGS\")\n",
    "plt.plot(result_gd.times, energies_gd, \"--\", color=\"royalblue\", label=\"Gradient descent\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.xlabel(\"time $t$\")\n",
    "plt.ylabel(\"energy $E$\")\n",
    "plt.title(\"Energy over time\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can observe here, that the energy does vary quite a lot! But as we mentioned before, p-VQD does not preserve the energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Magnetization over time')"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEYCAYAAABLOxEiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9iUlEQVR4nO3deZyN9fvH8dc1Y8YYhpFd1iS7sUxCthJRJJWQQiUpaVV8W4hQ4peSEllKyVKyVBQqRZGlhqyRdaxjMMYyzIzr98c5NGlwDufMfWbmej4e5zHnXs593vcM5zr3/bnvz0dUFWOMMcZTQU4HMMYYk7lY4TDGGOMVKxzGGGO8YoXDGGOMV6xwGGOM8YoVDmOMMV6xwmFMGiIyT0S6+GG760Skia+3G0hEpJOIzHc6h/E/sfs4jL+IyHagOFBcVQ+mmR8DRAFlVXW7I+FcOV4FrlXV+3283Y+AWFV92ZfbDSQiUgbYBoSoaorDcUwGsyMO42/bgI5nJ0SkGpDLuTjGWyKSw+kMJrBY4TD+9gnQOc10F2BS2hVE5HYR+UNEjorILveRQNrlnUVkh4jEi8grIrJdRG5xL3tVRKaLyCQRSXSfEopO89riIjJDROJEZJuIPOme3wJ4EWgvIsdEZLV7/iIR6eZ+vtq97OxDz55uEpHPRWSfiCSIyM8iUsU9vzvQCXjB/Zqv3PPTZs4pIm+LyB73420Ryele1kREYkXkORE5ICJ7ReTBC/1y3fs3R0QOicgWEXkkzfyTInJVmnVrishBEQlxTz8kIhtE5LCIfCcipdOsqyLSU0Q2A5vTeeuf3T+PuPeznoh0FZEl523jcRHZ7P7bvCYi5URkqftvPV1EQtOs30pEYkTkiIj8KiLVL7TfxmGqag97+OUBbAduATYBlYBgYBdQGlCgjHu9JkA1XF9kqgP7gTvdyyoDx4AGQCgwHEgGbnEvfxVIAm5zb/91YJl7WRCwCujnfu01wFbg1jSv/fS8zIuAbunsS3dgI5DXPf0QEAHkBN4GYtKs+xEwKL3fhfv5QGAZUBgoBPwKvJbmd5HiXifEvV8ngPwX+B3/BLwPhAE1gDigqXvZD8AjadYdBnzgfn4nsMX9d8kBvAz8mmZdBRYAVwG50nnfMu51cqSZ1xVYct425gB5gSrAKeB7998hH7Ae6OJetxZwALjB/Xfs4v6d5XT637E9/vuwIw6TEc4edTTD9eG7O+1CVV2kqn+q6hlVXQNMARq7F98DfKWqS1T1NK4icH7D3BJVnauqqe73inLPvx4opKoDVfW0qm4FPgQ6eBNeRBoAg4A7VPWoO/MEVU1U1VO4ClCUiOTzcJOdgIGqekBV44ABwANplie7lyer6lxchbNCOrlK4iqofVQ1SVVjgHFptvUZ7tOEIiK49vsz97JHgddVdYO62iiGADXSHnW4lx9S1ZMe7ld6hqrqUVVdB6wF5qvqVlVNAOYBNd3rPQKMUdXfVDVVVT/GVWjqXsF7Gz+xwmEywifAfbi+kU46f6GI3CAiP7pPJyUAPYCC7sXFcR2lAKCqJ4D48zaxL83zE0CY+7x8aaC4+9THERE5guv0VBFPg7s/nKfj+mb8l3tesIi8ISJ/i8hRXN+MSZP5UooDO9JM73DPOyte/93gfALIc4HtHFLVxPO2dbX7+RdAPREpDjTCVXAXu5eVBt5J83s5BEia10Ka3/sV2J/m+cl0ps/uV2ngufP+ViX59+/FBAhr9DJ+p6o7RGQbrtMuD6ezymfAKKClqiaJyNv88yG8lzTftkUkF1DAw7feBWxT1fIXinaxF7vfaxbwtqrOS7PoPqANrtNw23GddjmM64P3ktsF9uD6oFznni7lnuetPcBVIhKRpniUwn1Ep6pHxHV57L24TklNUdWz2XYBg1V18kW2f7H98PXlmGfzDPbxdo0f2BGHySgPAzer6vF0lkXg+uacJCJ1cH0wn/UF0FpE6rsbUgfwzwf0pSwHjopIHxHJ5T5SqCoi17uX7wfKiMiF/h9MADaq6pvp5D2F68gnHNdpnrT24zqPfyFTgJdFpJCIFMR1+u1TD/fpHFXdhat95HURCXM3Jj8MpC0Gn+E6TXg3/5ymAvgA+F+aRv18ItLOi7ePA85w8f30xodAD/fRp4hIbnFdNBHho+0bH7LCYTKEqv6tqisvsPhxYKCIJOL6EJ2e5nXrgF7AVFxHH4m4GlFPefCeqUBrXI3G24CDuNoAzrZFfO7+GS8iv6eziQ5A2/OurGqI63TbDlzf7NfjauhOazxQ2X3KZVY62x0ErATWAH8Cv7vnXY6OuBqq9wAzgf6quiDN8jlAeWC/qq4+O1NVZwJDganu021rgZaevqn7lOFg4Bf3fl5RW4T738YjuI48D+NquO96Jds0/mM3AJpMRUTyAEeA8qq6zeE4xmRLdsRhAp6ItBaRcBHJjety3D/5p0HaGJPBrHCYzKANrlMxe3CddumgdqhsjGPsVJUxxhiv2BGHMcYYr2SL+zgKFiyoZcqUcTqGMcZkGqtWrTqoqoXSW5YtCkeZMmVYufJCV4IaY4w5n4jsuNAyO1VljDHGK1Y4jDHGeMUKhzHGGK9kizYOY8zlSU5OJjY2lqSkJKejGD8JCwujRIkShISEePwaKxzGmAuKjY0lIiKCMmXK4BrSw2Qlqkp8fDyxsbGULVvW49cFzKkqESnpHpNhg7iG/3wqnXVEREa6h8hcIyK1nMhqTHaRlJREgQIFrGhkUSJCgQIFvD6iDKQjjhTgOVX93d2V8ioRWaCq69Os0xJXlxPlcQ0xOdr90xjjJ1Y0srbL+fsGTOFQ1b24us1GVRNFZAOu0cjSFo42wCR3P0XLRCRSRIq5X+tza58dwvhvZnJXyQqUiyjAvpQkvt67hTbNbqVsVDUOJRxh5fLlRFerTsFyZTlTMB85ShYjV/HCSFDAHMwZY4xPBUzhSEtEyuAai/i38xZdzb+Hs4x1z/tP4RCR7kB3gFKlSl1Wjr8+ncnbcSupuu8ExUOvYlNyAv0T/qTs5gPkCY3kl9NH6HJ0LVPyVeP6kHz8dPowDx9dxxeFr+f60uWIOZ3AmK0xvHhtHcqUKsXhQhHsjQzjhi73ElHlOvsmZ4wHgoODqVatGqpKcHAwo0aNon79+mzfvp1KlSpRocI/w7EvX76c0NBQvv32W/r168fRo0cJCwujQoUKDBs2jFKlSrFs2TKeeuopTp06xalTp2jfvj2vvvqqczuYCQVc4XCPtzADeFpVj56/OJ2XpNtLo6qOBcYCREdHX1ZPjncdWEHqmTNoaiqSkkqJxOPs3x9H7vBchAaHkHvbdr5dtISKufOT+9gpKv39Ny+siaRCnUbkDArh0LJf2HgygdQde0jaFc+cpH28dGwLi8bPoWREfr4PPsaMA1sZVv0mClUuT1LF0uS76QYia1a2IxZj3HLlykVMTAwA3333Hf/73//46aefAChXrty5ZWetXbuWXr16MWfOHCpVqgTAnDlz2L59O6VKlaJLly5Mnz6dqKgoUlNT2bRpU0buTpYQUIVDREJwFY3JqvplOqvE4hrA/qwSXN5YzR4LCgqCoCAICSFnrjAKF/5nuOuiZYpT9Kb656YLAvXSvLYr/wxhlnI4gfaLllJu/g9UrVID3XuQ47NmsPXUMUL+2snxzbG889kORp+MZW2NtuSLqsiqowc4rCm0eeRBcje+nqDcufy5q8YEvKNHj5I/f/6LrjN06FBefPHFc0UD4I477jj3/MCBAxQrVgxwHc1UrlzZP2GzsIApHOI6bzMe2KCqb11gtTnAEyIyFVejeIK/2jd8LUf+fJRr24JybVucm/f04Cd5GtDTyZxat4VWH0+l0NJfyVe1CqfXbWFMzEJWJh+l3pItAEzReFILRvL8uFGE1apMUHiYMztjsqWDL43k1NrNPt1mzqrlKTj4yYuuc/LkSWrUqEFSUhJ79+7lhx9+OLfs77//pkaNGgDceOONvPfee6xbt47evXtfcHvPPPMMFSpUoEmTJrRo0YIuXboQFmb/l7wRMIUDuBF4APhTRGLc814ESgGo6gfAXOA2XOMRnwAezPiYviehIYTVrETzmgNonmb+pD372DJnPhEb95L0+zoWr9zE8a17ub/tU5AjmEEnt1MyIpKebdqR5+5mhDesheQIpD+pMVcu7amqpUuX0rlzZ9auXQukf6oqrfj4eJo2bcqJEyfo3r07vXv3pl+/fnTq1In58+fz2WefMWXKFBYtWuT/HclCAuZTRlWXkH4bRtp1FOiZMYmcF1m8KNE9Op+b/g5I2n+Q1NV/kbRiLTuH9iM4/iDHpn3LsWnf0uPoBlpXrE63d4YSVrc6QWE5nQtvspxLHRlkhHr16nHw4EHi4uIuuE6VKlX4/fffiYqKokCBAsTExDB8+HCOHTt2bp1y5crx2GOP8cgjj1CoUCHi4+MpUKDABbdp/s1aYDOZsCIFyd28PgVe6s73R/cx9tDfFPl4MNqyHgk5IHFbLHvbPcvacs1pmqsQX9ZpQ8L4Lzlz7ITT0Y25Yhs3biQ1NfWiH/IvvPACgwcPZsOGDefmnTjxz7//b775hrMjn27evJng4GAiIyP9ljkrCpgjDnN5gsJykue2Rlx7WyOWT3qTM8dPcnLpalZNnMq26Ss5vmkrB/uOYOXzQ5geksjTL71IxZ73I8HBTkc3xiNn2zjA1UXGxx9/TPBF/v1Wq1aNd955h86dO5OYmEiBAgUoVaoUAwYMAOCTTz7hmWeeITw8nBw5cjB58uSLbs/8V7YYczw6Olqz60BOpzbvJHHiTKZOmcKz25fzff5oShYtyuaS+TmUU2g95CXCq5Z3OqYJUBs2bPjX1Ukma0rv7ywiq1Q1Or317Ygji8tZvhQ5hzzFE0Oeosu+A7B0DSfmLmbCZx/y9ckD/PbrZnLlz0diw+qUf6E7uSpe43RkY0yAszaObCSiaGEi2t5CkQ8HMC52PV8PHE5EzUqcOXqchz8Zzc3VaxHXdwRJv68n9aR1o22MSZ8Vjmwqd/5IGr/8NCUXjKPsjvk82eMxHoqqR+KnX7OzeXcaRRbjg+sacmTCTPTUaafjGmMCiBUOQ1BYTh549w0eWTGP0utmoX27kBwSTPK+g8T3eYvVpW5mTL3WHFtjXTMYY6xwmPME54ugQu9HWH7sII/+sYA89zRnfvIheiz7mgUtu3Dw1fc4vnyNXd5rTDZmhcNcUGi5EhQZ/QrP7oth3uAR3HDLzSSM+ZzeTW6ncYES7Gz3LKfW+bYLCmNM4LPCYS4pR2goLV58muIfD6F0zAxK14+mTGhukhetILbJQ4wv14Bt46c7HdNkUfv37+e+++7jmmuuoXbt2tSrV4+ZM2de0TZfffVVhg8fDkC/fv1YuHDhZW0nJiaGuXPnerRukyZNcPK2gFmzZrF+/fpLr+gBKxzGKzmKFKDvwplMOrKdwu+/zPFShXl866/07/UMB3oN4dTaLXYay/iMqnLnnXfSqFEjtm7dyqpVq5g6dSqxsbH/WTclJeWy3mPgwIHccsstl/VabwqH06xwGMdJcDAR7W6l+qoZ/DZ7Ls90uJ9jc35kUcP2NC5QgiVtHuX09t1OxzSZ3A8//EBoaCg9evQ4N6906dL06tULgI8++oh27drRunVrmjdvzrFjx2jatCm1atWiWrVqzJ49+9zrBg8eTIUKFbjlllv+NQZH165d+eKLLwBYtWoVjRs3pnbt2tx6663s3evqfLtJkyb06dOHOnXqcN1117F48WJOnz5Nv379mDZtGjVq1GDatGn/yn7y5Ek6dOhA9erVad++PSdPnjy3bP78+dSrV49atWrRrl27c/1o9e3bl8qVK1O9evVzPfzu37+ftm3bEhUVRVRUFL/++isAn376KXXq1KFGjRo8+uijpKamApAnTx5eeukloqKiqFu3Lvv37+fXX39lzpw5PP/889SoUYO///77iv4udgOguWI17mgBd7Qg9Ugiyx/vw87pGwhdsppd13fgaNWyXPtGbyJuqO50TOMDu9v0+s+8PG1uJt9DbTlzIom9HZ//z/KIDi3J2/E2UuOPsO+hV/617OrZ7170/datW0etWrUuus7SpUtZs2YNV111FSkpKcycOZO8efNy8OBB6tatyx133MHvv//O1KlT+eOPP0hJSaFWrVrUrl37X9tJTk6mV69ezJ49m0KFCjFt2jReeuklJkyYALiOaJYvX87cuXMZMGAACxcuZODAgaxcuZJRo0b9J9fo0aMJDw9nzZo1rFmz5tx+HDx4kEGDBrFw4UJy587N0KFDeeutt3jiiSeYOXMmGzduREQ4cuQIAE8++SSNGzdm5syZpKamcuzYMTZs2MC0adP45ZdfCAkJ4fHHH2fy5Ml07tyZ48ePU7duXQYPHswLL7zAhx9+yMsvv8wdd9xBq1atuOeeey76+/SEFQ7jM8GREXT47H3aTRhB4qTZHP6/STz3y9fENZjJT2+NIfLhu2wMEXNFevbsyZIlSwgNDWXFihUANGvWjKuuugpwndp68cUX+fnnnwkKCmL37t3s37+fxYsX07ZtW8LDw4F/D+x01qZNm1i7di3NmjUDIDU19dyATwB33XUXALVr12b79u2XzPrzzz/z5JOuHoWrV69O9equL0/Lli1j/fr13HjjjQCcPn2aevXqkTdvXsLCwujWrRu33347rVq1AlxHXZMmTQJcA0/ly5ePTz75hFWrVnH99dcDrqObwoULAxAaGnrutbVr12bBggUe/W69YYXD+FxwWE4iu99Lvkfa8Vi/19n9+TwODxzN4Xc/ZXbQUdr370PxDq1szPVM6GJHCEHhYRddHlwg8pJHGOerUqUKM2bMODf93nvvcfDgQaKj/+lCKXfu3OeeT548mbi4OFatWkVISAhlypQhKcnVC8Kl/r2pKlWqVGHp0qXpLs+Z0zVMQXBwsMftKem9p6rSrFkzpkyZ8p9ly5cv5/vvv2fq1KmMGjXqX4NWnb+NLl268Prrr/9nWUhIyLn39SarN6yNw/iNiNDutRd5euNiin/1Hn8Xz8tzm35h3CPPsKNaWxJnLCA7dLJpLt/NN99MUlISo0ePPjcvbRfp50tISKBw4cKEhITw448/smPHDgAaNWrEzJkzOXnyJImJiXz11Vf/eW2FChWIi4s7VziSk5NZt27dRfNFRESQmJiY7rJGjRoxefJkwDUO+po1awCoW7cuv/zyC1u2bDm3P3/99RfHjh0jISGB2267jbfffvvcAFVNmzY9t/+pqakcPXqUpk2b8sUXX3DgwAEADh06dG5fLyert6xwmAyRq251mi+aypJJU7mvdn1S98cz/aGnea5IJRJWb7j0Bky2JCLMmjWLn376ibJly1KnTh26dOnC0KFD012/U6dOrFy5kujoaCZPnkzFihUBqFWrFu3bt6dGjRrcfffdNGzY8D+vDQ0N5YsvvqBPnz5ERUVRo0aNcw3RF3LTTTexfv36dBvHH3vsMY4dO0b16tV58803qVOnDgCFChXio48+omPHjlSvXp26deuyceNGEhMTadWqFdWrV6dx48aMGDECgHfeeYcff/yRatWqUbt2bdatW0flypUZNGgQzZs3p3r16jRr1uxcQ/6FdOjQgWHDhlGzZs0rbhy3btWNI5JiNvLMnffyXewWvi3VgEJP3E94q0aEli9tp7ACiHWrnj142626HXEYR4TVqMjo7WtYvngJEfVqcmDIWG6sEsWocvU48ZMVeWMCmRUO46iC9WpRbPJQck56jdDwXIQePMree55ha4P7Of77xc8vG2OcYYXDBIQyt93EL0f20mnGOIKvLsxHMUupWacOW9+aYA3oxgQYKxwmYIgIeZrVp/QfX1D18c7UylOQM0MmsKf1E2z7v/GkHvHNFSHGmCtjhcMEHBHh3jf7MyXubwoN703cpi3Ufr4HL5eqQfyQD9HTyU5HNCZbs8JhApaE5CBf1zu5ZvEndL6+ITfmyMeREZP4s9Jt7PnsKzuFZYxDrHCYgJe3aGHe/u0HWqz9lrCGtXl7z3pqPXAPsUPHoWfOOB3P+FlwcDA1atQ493jjjTd8tu3M1LttILEuR0ymkaN4Ya7+8m26jI+m2KARnP6/Sez+fjkJ11eg4vPdCc6f1+mIxg9y5cp17i5qX4uJiWHlypXcdtttftl+VhVQRxwiMkFEDojI2gssbyIiCSIS4370y+iMxnmNH76fQVtXUnhMfzZt3Uqt119g2DU3cHjkZNTdtbTJ2hISEqhQocK57tE7duzIhx9+CLju2I6OjqZKlSr079//3GtWrFhB/fr1iYqKok6dOiQkJFy0W3RzEaoaMA+gEVALWHuB5U2Ar73dbu3atdVkTUf27NOnr79JlxW4QbcUbKB/VGmlR5escjpWlrF+/fp/TTdu3FgnTpyoqqqnT5/Wxo0b6yeffKKqqsePH9fGjRvr1KlTVVX1yJEj2rhxY50xY4aqqsbFxWnjxo11zpw5qqq6d+9ejzIEBQVpVFTUucfZ7c+fP1/r1q2rU6ZM0VtvvfXc+vHx8aqqmpKSoo0bN9bVq1frqVOntGzZsrp8+XJVVU1ISNDk5GSdOHGi9uzZ83J+NVnK+X9nVVVgpV7gMzWgTlWp6s8iUsbpHCbzyFesCCOW/8CpLTvZ/3A/ei/7isM33cSib+YR0ay+0/GMD1zoVFWzZs34/PPP6dmzJ6tXrz43f/r06YwdO5aUlBT27t3L+vXrERGKFSt2rhvyvHnttOaVCKjC4aF6IrIa2AP0VtV0by8Wke5Ad4BSpUplYDzjhJzXlqLkool0faE/sZ/O5sB9fTjRoSVapwpFOrVGggLqrGymtWjRonPPQ0JC/jUdHh7+r+l8+fL9a7pgwYL/mi5atOgVZTlz5gwbNmwgV65cHDp0iBIlSrBt2zaGDx/OihUryJ8/P127diUpKQlVtT7QfCiz/W/6HSitqlHAu8CsC62oqmNVNVpVowsVKpRR+YyDRIT7hw3khS3LiOzViQWfTqVyl3uYG9WSUxu3OR3P+NiIESOoVKkSU6ZM4aGHHiI5OZmjR4+SO3du8uXLx/79+5k3bx4AFStWZM+ePecGf0pMTCQlJcWnXY1nJ5mqcKjqUVU95n4+FwgRkYIOxzIBJih3Lgr060GVj4fSKP/VlNmbSGyjLux79k3OJJ1yOp7x0smTJ/91OW7fvn3566+/GDduHP/3f/9Hw4YNadSoEYMGDSIqKoqaNWtSpUoVHnrooXOj7IWGhjJt2jR69epFVFQUzZo1Iykp6aLdopsLC7hu1d1tHF+ratV0lhUF9quqikgd4AtcRyAX3QnrVj370jNnSBjzOfsHvk+Hg3/Qtti1vLJ5GZIz1OlomYJ1q549ZOpu1UVkCrAUqCAisSLysIj0EJEe7lXuAda62zhGAh0uVTRM9iZBQUQ+1p5Cv0yiZJEiXHXkJLuaPszJX2NITbBTFMZcjoBqHFfVjpdYPgoYlUFxTBZy1TWlmRO7meMLlnLw+eGMuvVeNkgSb374AQXb281fxngjoI44jPG33M3qUXLJJ+wqV5i1JxM43HMIu9v0IiXusNPRApYd1Gdtl/P3tcJhsp2gPOGMjFnC/PkLyFEgkrglK3nimhrsnT7P6WgBJywsjPj4eCseWZSqEh8fT1hYmFev8+hUlYjkB1oDbYHrgG3AbGCOqu73MqsxASHypjrkjZnBR/c8zPi5v3FrnyHcWqgA4TfVcTpawChRogSxsbHExcU5HcX4SVhYGCVKlPDqNZe8qkpEvgTyA9/gKhR/iUgpoA2uYhKqqk0uK3EGsauqzKVs/nohOYd8RPLmHfxSNIxmvR+nWJe2TscyxjEXu6rKk8IRqapHLnd5ILDCYTxxJukUG/83jJpv96NjWFEGt7iHouNfI/iqfE5HMybDXdHluJcqCoFeNIzxVFBYTiqPeJm5o8fzZPEqJC35g5hqd3Dwy/lORzMmoHjdOC4izUTkQxGp4Z7u7vNUxjioaY8HqfHnV4S3uYnecX/StP3dHP9phdOxjAkYl3Mfx+PAg8DLInIVUMOniYwJAEHhYRQbN5Cny+Vn88iP2dfhefI/15XwFjcSVrW80/GMcdTlXI4bp6pHVLU30By43seZjAkYbf73DM9sWEyeO27iq9eG06hmNGufGIAmpzgdzRjHXE7h+ObsE1XtC0zyXRxjAk9wZARFxvQnqOsdJAvI1O/YWe8+Tv+9y+loxjjikoVDRA6KyO1np1V1dtrlqvquP4IZE2geeGcIv/29iYiK5Ti5fTfPRTVgx5ipTscyJsN5csSRDIwWkW7nL3B3SmhMthFasiglf/qITa1v4IPju/h24HBO/hrjdCxjMpQnhWMvrrHAe4vIq+ctq+DzRMYEOAkOps3Et1gz51taXl2ePW2f4vsWXUi0K69MNuFRG4eqbgduBJqLyDgRsT6uTLZXqXVzSvwwnqSW9bhz/mR6tWjDwVdGoampTkczxq88KQACoKrxwM1AYWCOiISfXWZMdhWUJ5wqH73BW48+zYN5SpHwwTR2NOpKyl7r28lkXZ4Ujj/OPlHVJOBOYDewCMjrl1TGZDKPjB5O/eUzyFGqGH1WLaDbdbVIPmhdtZusyZMuRx46b/qMqj6K67Lcsv4KZkxmE1L2akr88gn5K15LeCrsa/ccp//ehZ4543Q0Y3zqstsqVHUAVjiM+ZfgsJx8ELOEN2dOJWX3fhbWv4f3rmtI8q59TkczxmeuqJFbVXf4KogxWUmeW2+k5I8TmKhxDNy6gvV17uXYV4ucjmWMT1xWJ4f+CGJMVpPj6iJM3r6O2V2fJk+Ksu/Bl/nzkZesuxKT6V3OEcdQn6cwJosKDc9FvQlvUmTiIKakxFF//FCWPT3A6VjGXBG7H8OYDJCnVWM6zvuUzsUrUHDajxx85V1S9sc7HcuYy+JR4RCRiSIyQUQmAqXczyeIyAQ/5zMmy7i23vWM2hpDZLe72f7+ZNqWrsSfzw2xq65MpuPpeBwfpXneAPjY91GMyfokNIRCbzzDshxJLBnSj83jp5N/3Q6KTR1OcGSE0/GM8cglxxz/zwtEflfVWn7K4xc25rgJRPG/r+Xo/S9xJu4Qa8JSafHVeHLVqOR0LGOAKxxzPB2nrzCPMQYoUKsqpZdPYWu1UtwTu5QhrTtyJvG407GMuSSvC4eq1vVHEAB3u8kBEVl7geUiIiNFZIuIrBGRTHXkY8z5gvKE03ThJ4xo15V2SXmIbd6dE4uWo6eTnY5mzAUF2lVVHwEtLrK8JVDe/egOjM6ATMb4VVBQEE9On0j5We9y6vARWt/akrFVbyZl30GnoxmTroAqHKr6M3DoIqu0ASapyzIgUkSKZUw6Y/wr1401yfvF/3E4REjce4CdN3Tk5K9/XPqFxmSwgCocHrgaSDvQc6x73n+ISHcRWSkiK+PirItrkzkUqVqR5Qd2cf+97dETSXxz+4NsfWMs3l7EYow/XXHhEJE+vgji6dulMy/d/1GqOlZVo1U1ulChQn6OZYzv5AjPRdHxAwnv/yiPJaynx4CXOPXnZqdjGXOOp/dxnCMi09NOAjXIuG5IYoGSaaZLAHsy6L2NyVDFnrifz3NB0Osfsaf1ExQY8iS5b76BHMXsi5Bx1uUccRxV1Xvdj3bAQl+Huog5QGf31VV1gQRV3ZuB729Mhmry8P3UXzKV0Crl6NWtO09VvIETy1Y7Hctkc5csHCJyUERuTzNr8HmrvOSrMCIyBVgKVBCRWBF5WER6iEgP9ypzga3AFuBD4HFfvbcxgSpH0YIUmzGC1JJFOJ18mr139CJhwkynY5ls7JJ3jovIXiAZGKiq485bNkVVO/oxn0/YneMmKzhz5gyH3/2MI4PGsD31JOEtG1L342FIDq/POBtzSVd65/heoBHQW0RePW9ZhSvMZozxUFBQEAWeup+iX7xF3+NbuG/aGBLm/Oh0LJMNedTGoarbgRuB5iIyTkQy22W8xmQZuRtfz+Tvv+Wtyg2Jf2wQh9+fyqktO52OZbIRTwqAAKhqPHAzUBiYIyLhpH95rDHGz65rUJc2v84kd4sbeb/PKzxQ/QYOz/7B6Vgmm/CkcJy7dVVVk4A7gd3AIiCvX1IZYy4pKE84RSYOIrHmtRxIOcWBh/txaNgEu1nQ+N0lC4eqPnTe9BlVfRT4Bijrr2DGmEuToCDe+HkeMz/6mBzBQex840N+avMIeso6sTb+48nluOmuo6oDgLLW3mGM8/Lf15ri37zPoFM7ufPrj9n54TSnI5kszJMP/QUiMk1EOopIXgARCReRu4BBwCq/JjTGeCRXdFX+7+d5vFmpASkDP+TI6Kkk77N+2ozveXKqqikwACgDfCMiy4DvgerACFWt6deExhiPla5VnW4r5pH79kbM+98QOl8TxeGv7JJd41se3TmkquuB9cDrIhLmbiQ3xgSgoPAwiowfyMYda1jx8xa2d30R+vUk8slOiNiFkObKXc4IgFY0jAlwEhTEgB/msPiDSeQNDiHutQ9Ydf+zaEqK09FMFmAN28ZkYcW63UOx2SMZdXo3Tae8y/phY52OZLIAKxzGZHHh9WrwxDdTeKx4ZcLe/ZyjU+ba5brmiljhMCYbuK5hXQav+4Vc9aKI6dmf54pW4fhva5yOZTIpjwuHiOQUkftE5EUR6Xf24c9wxhjfCc4XQbGpw/n2mgjGJWxj1e2PkDhjgdOxTCbkzRHHbKANkAIcT/MwxmQSEpKDQcsXsrjPEK6WUA70GMiuge9ZNyXGK9505F9CVVv4LYkxJkMEBQUR9foLJF5fkzFdezFywHN8nTOI6n0eczqaySS8OeL4VUSq+S2JMSZDRdzVjBve7U+NXPkJHTOTk7/8cekXGYN3haMBsEpENonIGhH5U0Ssdc2YTKxhl47MWLeK8GKF2XLXk4yrfBMpu/c7HcsEOG9OVbX0WwpjjGNCShfn6m9G807Nxry94Q/K3nA3N86dQFj165yOZgKUx0ccqroDiARaux+R7nnGmEwuODKCIX8u5rOb21PmlLD71u4cn/+r07FMgPLmctyngMm4RgAsDHwqIr38FcwYk7Fy5snNvQunkK/HvcScPELrVq2Infq107FMAPKmjeNh4AZV7aeq/YC6wCP+iWWMcYKIUPC1XhxpfzM7zyRxYMTHpOw54HQsE2C8KRwCpKaZTsXGHDcmS3p4zFusmreQyANH2dmsOyseeRE9nex0LBMgvCkcE4HfRORVEXkVWAaM90sqY4zj8jWrz9VfvceY/Ru4afybLGnUkdSjx5yOZQKAN43jbwEPAYeAw8CDqvq2n3IZYwJAzqrX8vh303m8SAWKbdnHroZd7NSV8a6TQ1VdpaojVfUdVbW7hYzJBkrXjuL1Tb+Rs1p59uzaRf9qjUhat8XpWMZBlywcIrLE/TNRRI6meSSKyFFfhhGRFu4bDLeISN90ljcRkQQRiXE/rJNFYzJAcN48lPhuLDNLhPLO4S38OXiUDQqVjXky5ngD988IVc2b5hGhqnl9FUREgoH3cN1oWBnoKCKV01l1sarWcD8G+ur9jTEXJ6EhDI35me+f7MdVC1axt/OLJC5c6nQs4wBv7uMY6sm8K1AH2KKqW1X1NDAVV2+8xpgAERQUxA1v96Pgm8/yxddf0bBlc/56ZYT1rpvNeNPG0Sydeb7shuRqYFea6Vj3vPPVE5HVIjJPRKpcaGMi0l1EVorIyri4OB/GNMbke7AtxZ56gHxBIaSM/oIDvYagqamXfqHJEjxp43hMRP4EKrg7Nzz72Ab86cMs6d0Tcv7XmN+B0qoaBbwLzLrQxlR1rKpGq2p0oUKFfJfSGANAu9deZN5PP5IzVxjxU+eysHlnziSdcjqWyQCeHHF8hqtvqjn8009Va6C2qnbyYZZYoGSa6RLAnrQrqOpRVT3mfj4XCBGRgj7MYIzxQnjdKEr+OJFhKbu584epbB75sdORTAbwpHE8QVW3q2pH4ChQBCgNVBWRRj7MsgIoLyJlRSQU6ICrWJ0jIkVFRNzP67jzx/swgzHGS6HXlmLQknm8fl19crz1GUenzOXMiSSnYxk/8qZxvBvwM/AdMMD981VfBVHVFOAJ93Y3ANNVdZ2I9BCRHu7V7gHWishqYCTQQa1VzhjHXV21Ek+s/JZcDWqy6PGX6FuyOknr/3Y6lvET8fRz193OcT2wTFVriEhFYICqtvdnQF+Ijo7WlStXOh3DmCxPTyfTo1o9vtqylq8KX0/lGSPJVb+m07HMZRCRVaoand4yb66qSlLVJPcGc6rqRqCCLwIaY7IGCQ1h9PrfmN+tN5EpsPvOp4j//DunYxkf86ZwxIpIJK4rmRaIyGzOa7w2xpig4GCqjhlE/j4P8d7xHdx6373E/bTc6VjGh7zp5LCtqh5R1VeBV3D1jGs36Blj0nVV7wep0e0+rg0J5+Qro0jZd9DpSMZHvGkc/15EbgNQ1Z9UdQ6uLkKMMSZdXd55nQnzviJl5z5iGt/Pzx2ftBsFswBvTlWVBfqISP8089JtODHGmLPCG0dz9ex36btrFfdOH8PWO3rajYKZnDeF4wjQFCgiIl+JSD7/RDLGZDU5oyow+ruveLNYDVi+jthmj5CakOh0LHOZvBo6VlVTVPVxYAawBCjsn1jGmKymfMO63Pf7PIJLFWPe6pUMr9qY1PgjTscyl8GbwvHB2Seq+hHQFZjv4zzGmCwsR+GrKLVoIrNyHmfKvs3EfzzLetbNhHJ4uqKqjjlvehWuoWSNMcZjQRG5mbV1Hdt7DiTx9fEk79pH/k6tCI+u6nQ046GAGgHQGJM95IzIw3UfvUHeh+/i8dFv0b7BzSR+u8TpWMZDATMCoDEme5GgIAq+/jS1GjfkuuBwDnR+kaOfzXU6lvFAII0AaIzJZkSE/t/N4JURw0GVZT1f4a+Bo5yOZS4hkEYANMZkU5Hd76HAB6/QI3ED9w9+maTNO5yOZC4ikEYANMZkY5F3N2fcyPd4tWBl9t39DKfWbkZTUpyOZdJxyW7V3Tf65QdeB/qmWZSoqof8mM1nrFt1YzKPU+u2sKfds3y0ewNR1apz58LPCArL6XSsbOeKulU/OwIg0AloCHRR1R1AHvcofMYY4zM5q1xLwZlv81nSPj5ZuZjY5t1JPXrM6VgmDW/aON4D6gEd3dOJWCeHxhg/iKhwDUtW/8EbFRuQvGEr2xp1JmW/jRIdKLwpHDeoak8gCUBVDwOhfklljMn2ilS8lvKLPyGlYhk6rPuefjWb2F3mAcKbwpEsIsGAAohIIeCMX1IZYwyuu8xLf/M+pYsWo/ixFA4NGmPFIwB4UzhGAjNx9Y47GFcnh0P8ksoYY9zC8kbwxfZ1tO/+MEdGTmbBje1I+GGZ07GyNW/6qposIqtwda0OcKeqbvBPLGOM+YcEB1Nw2HPEBaVy1/CXuLfVMt6ZMomINjc7HS1b8rhwiEhOoBaQz/26diKCqg70VzhjjDlLRKj8Zh+G795G1Ld/cKBbf1LjjxD50F1OR8t2vDlVNRvXGOMpwPE0D2OMyTA9Jn9A5ZH9UJSBvZ5hbd9hTkfKdjw+4gBKqGoLvyUxxhgP5et0OztPH2fsw/cROn4Cr/XoREiZ4k7Hyja8OeL4VUSq+S2JMcZ4odqD97Jy1jd0jyzL7tsfI+HLhdZFSQbxpnA0AFaJyCZ3X1V/isgafwUzxphLqXjHrZT45n32p5yi3r138vmNd3Hm5CmnY2V53hSOlkB5oDnQGmjl/ukzItLCXZi2iEjfdJaLiIx0L18jIrV8+f7GmMwn9LoyFJ08lFyhoeTYtJPYZt2sixI/87hwqOqO9B6+CuK+ufA9XAWqMtBRRCqft9rZ4lUe6A6M9tX7G2Myr9LRUfy2829ql69I8qbt/FinrXVR4kfeXI77bDqzE4BVqhrjgyx1gC2qutX9flNxXcW1Ps06bYBJ6rp1dJmIRIpIMVXd64P3N8ZkYjkK5qfkjxP4slE72v8xl1GtH+Cx5TaioD94c6oqGugBXO1+dAeaAB+KyAs+yHI1sCvNdKx7nrfrACAi3UVkpYisjIuL80E8Y0ygC8oTzh0/T+eFag1ouPUI8QNHcyY11elYWY43haMAUEtVn1PV53AVkkJAI6CrD7JIOvPO75TGk3VcM1XHqmq0qkYXKlToisMZYzKH0Dy5GfLHIoo8dA/7R37K/4pU4ciPvzkdK0vxpnCUAk6nmU4GSqvqScAXlzHEAiXTTJcA9lzGOsaYbE6Cgyk49BlW3VSZYfGbmNm2G4mzf3Q6VpbhTeH4DFe7Qn8R6Q/8AkwRkdz8ux3icq0AyotIWREJBToAc85bZw7Q2X11VV0gwdo3jDHpEREe+PxDlvR5nUYhkRzo1o/D4790OlaW4M1VVa8BjwBHcDWK91DVgap6XFU7XWkQVU0BngC+AzYA01V1nYj0EJEe7tXmAluBLcCHwONX+r7GmKyt/ht9KTiyLxtTT9C0R1fW/d+HTkfK9LzpcgRcH9rBQBgQLiKNVPVnX4VR1bm4ikPaeR+kea5AT1+9nzEme8jX8XZC9uzgWN/nOTJxFsl33UZI2XSvqzEe8PiIQ0S6AT/jOiIY4P75qn9iGWOMbzV/7nH++GUpxU/DrpY9+KPXAOui5DJ508bxFHA9sENVbwJqAnadqzEm08hdpzpXf/0enx7axo2jBrLo5vuti5LL4E3hSFLVJHCNzaGqG4EK/olljDH+EXpdGR76bgrdCl5LifWx1kXJZfCmcMSKSCQwC1ggIrOxS2GNMZlQiZrVGLFpOaHXlGD/hs28XbUJKQcOOR0r0/Dmqqq2qnpEVV8FXgHG4eoCxBhjMp3gq/JR8seJfBKRxIu7fufP197Fdf2NuRRv+qqKBl4CSrtfJ8BgoLp/ohljjH8F5c7F/61fSqtOT5N36kIO5o0kskd7QkoWdTpaQPPmVNVkYCJwN37qVt0YYzJaSHgums0YTb7u7Zg9aixtrqvOofm/OB0roHlTOOJUdY6qbvNHt+rGGOMUCQqiwKBeHGpQjdiUk+zs9DyJM793OlbA8qZw9BeRcSLSUUTuOvvwWzJjjMlAIsIL30zl+7fHEEEwux/px5bh452OFZC8uXP8QaAiEAKccc9TwDp/McZkGYV63kdY8cL0uK8zi/s8wYqoihRpdqPTsQKKN4UjSlWr+S2JMcYEiIi2t/DoiNcp9b/BnOz9FqenFSe0YlmnYwUMb05VLUtnKFdjjMmSGvXowqtL5kJKCj817cT0Ju3Q08lOxwoI3hSOBkCMiGwSkTUi8qeIrPFXMGOMcVrOKtdy9dzRvJHwN08v/oq/m3fjzLETTsdynDenqlr4LYUxxgSokNLFmb76N1a16o6s28quxl0p8d0YggvmdzqaY7y5c3xHeg9/hjPGmEBQoFwZblkxi9Cq1zJ+0yp6Vbie5MMJTsdyjDenqowxJtsKyhPO1d+OYVvh3Px97Ajxb4xHz5y59AuzICscxhjjoaCcoXy8+Q8mPvsixyfMZNPdT3LwuyVOx8pwVjiMMcYLQcHBFBvyDPn79eDBbz7l9tatOPr5d07HylDeDh1rjDHZnohwVa9OPLYhhsOffkXc44NI3RdP/l73OR0tQ9gRhzHGXKbO7w/joc/GQlAQc19+nZ8efiFbdM1uhcMYY65ARKtGFJn1Dv1PbKX3J2M5ufxPpyP5nRUOY4y5Qnnq1WDe9wt4t3wD9t3bm2Ozf8zSY5lb4TDGGB8o36g+0T98THCZ4jzWoRODKtUn9fBRp2P5hRUOY4zxkRxFClB05tvE5c7B7gP72Vm/E8m79jkdy+escBhjjA/lvCqSr3f9xcst7ubMwSOsqNuOxD/WOR3Lp6xwGGOMj+XIFUaJL0YQcs8tdNq3kvsa3kLqiSSnY/lMQBQOEblKRBaIyGb3z3R7DxOR7e5eeWNEZGVG5zTGGE9JUBClRvfnubs60jlHYfbd8wwp8UecjuUTAVE4gL7A96paHvjePX0hN6lqDVWNzphoxhhz+Z74fAItP36bpJiNjKrcmGXPDnI60hULlDvH2wBN3M8/BhYBfZwKY4wxvpSnzU2czJmDN26/hRveH8nYJKXwqJeQoED57u6dQEldRFX3Arh/Fr7AegrMF5FVItL9YhsUke4islJEVsbFxfk4rjHGeKdQi4b8NG8+gwpX5djn84m980n01GmnY12WDCscIrJQRNam82jjxWZuVNVaQEugp4g0utCKqjpWVaNVNbpQoUJXnN8YY65UheaNqfjbNM4Uzk+nbyczsEpDpyNdlgwrHKp6i6pWTecxG9gvIsUA3D8PXGAbe9w/DwAzgToZld8YY3whpFghiv84gdx58xK2/zCHhk3MdP1bBcqpqjlAF/fzLsDs81cQkdwiEnH2OdAcWJthCY0xxkdyFy7IN7s307VzFw6/OYHvm3Tk8LIYp2N5LFAKxxtAMxHZDDRzTyMixUVkrnudIsASEVkNLAe+UdVvHUlrjDFXKChnKIVHvUhQ97u4d/EMHrmpBcfnLXY6lkcksx0iXY7o6GhdudJu+zDGBKZJ3Z6h1JdLKBkcRoHXehHZ416nIyEiqy5020OgHHEYY0y21XncCKI/HgZBQfR/7nnmPfCU05EuygqHMcYEgDytm5B7+pvMOhXHjFkzOfFT4J4lscJhjDEBomjjG1j+xypeimrE3g692dZ/JCkJiU7H+g8rHMYYE0CKVKtEqbkfkFy7Is0G96XXNTVJjg2srtmtcBhjTIAJzpuH8jNG0qJ8VRqk5GJXvU4krQycuw+scBhjTAAKzhnKu+t/o3nPbmjSaT66pT2bPpzqdCzACocxxgQsEaHQkKcIGdCD/x39i77PPMvpv3c5HcsKhzHGBLpSj3fiu/Gf0L9IVXa37EHiNz+jqamO5bHCYYwxmUCdBztQZcEENF9u7rqrLa9VupEzx044ksUKhzHGZBIh15Sg8Ox3CcqdC91zkJ03dCR5T7p9wvqVFQ5jjMlE8hQvwry9W3mkRWtSDxxifs3W7PtpWYZmsMJhjDGZTFDOUIp9/hYhXe/g4YMxPHBbG5L3x2fc+2fYOxljjPEZEaHUsOcZ+8IrvBBRlj0tHiUpZmOGjO1hhcMYYzKxO15/mUbffcyZ5BSeb3ArkxvdjSan+PU9rXAYY0wmlzOqAkW+HsXeHEqF9XvYeeP9pBxK8Nv7WeEwxpgsILxMCabv3kTRWtVI2babndHtObVhq1/eywqHMcZkEUG5c1HiuzFEtG+BJh5n/6MD/HKvRw6fb9EYY4xjJCiIwqNeIqxBLZKWrkZy5/L5e1jhMMaYLChvh5bk7dDSL9u2U1XGGGO8YoXDGGOMV6xwGGOM8YoVDmOMMV6xwmGMMcYrVjiMMcZ4xQqHMcYYr1jhMMYY4xXJiC54nSYiccCOy3x5QeCgD+NkBrbPWV9221+wffZWaVUtlN6CbFE4roSIrFTVaKdzZCTb56wvu+0v2D77kp2qMsYY4xUrHMYYY7xihePSxjodwAG2z1lfdttfsH32GWvjMMYY4xU74jDGGOMVKxzGGGO8YoUDEJEWIrJJRLaISN90louIjHQvXyMitZzI6Use7HMn976uEZFfRSTKiZy+dKl9TrPe9SKSKiL3ZGQ+f/Bkn0WkiYjEiMg6EfkpozP6mgf/tvOJyFcistq9zw86kdNXRGSCiBwQkbUXWO77zy9VzdYPIBj4G7gGCAVWA5XPW+c2YB4gQF3gN6dzZ8A+1wfyu5+3zA77nGa9H4C5wD1O586Av3MksB4o5Z4u7HTuDNjnF4Gh7ueFgENAqNPZr2CfGwG1gLUXWO7zzy874oA6wBZV3aqqp4GpQJvz1mkDTFKXZUCkiBTL6KA+dMl9VtVfVfWwe3IZUCKDM/qaJ39ngF7ADOBARobzE0/2+T7gS1XdCaCqmX2/PdlnBSJERIA8uApHSsbG9B1V/RnXPlyIzz+/rHDA1cCuNNOx7nnerpOZeLs/D+P6xpKZXXKfReRqoC3wQQbm8idP/s7XAflFZJGIrBKRzhmWzj882edRQCVgD/An8JSqnsmYeI7w+edXjiuKkzVIOvPOv0bZk3UyE4/3R0RuwlU4Gvg1kf95ss9vA31UNdX1ZTTT82SfcwC1gaZALmCpiCxT1b/8Hc5PPNnnW4EY4GagHLBARBar6lE/Z3OKzz+/rHC4qm/JNNMlcH0T8XadzMSj/RGR6sA4oKWqxmdQNn/xZJ+jganuolEQuE1EUlR1VoYk9D1P/20fVNXjwHER+RmIAjJr4fBknx8E3lBXA8AWEdkGVASWZ0zEDOfzzy87VQUrgPIiUlZEQoEOwJzz1pkDdHZfnVAXSFDVvRkd1Icuuc8iUgr4EnggE3/7TOuS+6yqZVW1jKqWAb4AHs/ERQM8+7c9G2goIjlEJBy4AdiQwTl9yZN93onrCAsRKQJUALZmaMqM5fPPr2x/xKGqKSLyBPAdrisyJqjqOhHp4V7+Aa4rbG4DtgAncH1jybQ83Od+QAHgffc38BTNxD2LerjPWYon+6yqG0TkW2ANcAYYp6rpXtaZGXj4d34N+EhE/sR1GqePqmba7tZFZArQBCgoIrFAfyAE/Pf5ZV2OGGOM8YqdqjLGGOMVKxzGGGO8YoXDGGOMV6xwGGOM8YoVDmOMMV6xwmGMMcYrVjiMMcZ4xQqHMVdIRCJF5PE007/66X1KiEh7f2zbGG9Y4TDmykUC5wqHqtb30/s0xTXugjGOssJhzJV7AyjnHkVvmIgcAxCRMiKyUUTGichaEZksIreIyC8isllE6pzdgIjcLyLL3dsYIyLBad9ARBoAbwH3uNcpm6F7aEwa1uWIMVdIRMoAX6tqVff0MVXN456/BagJrMPVAd9qXN3U3wE8qKp3ikgl4E3gLlVNFpH3gWWqOum89/kW6J2Z+5IyWUO27+TQGD/bpqp/AojIOuB7VVV3B3tl3Os0xTUmxgp3h5K5SH8EwgrAJr8nNuYSrHAY41+n0jw/k2b6DP/8/xPgY1X934U2IiIFcHWHneyXlMZ4wdo4jLlyiUDEFbz+e1xtF4UBROQqESl93jplydyDh5ksxAqHMVfIPTriL+4G8GGX8fr1wMvAfBFZAywAip232kZc4y2sFRF/XbVljEescdwYY4xX7IjDGGOMV6xwGGOM8YoVDmOMMV6xwmGMMcYrVjiMMcZ4xQqHMcYYr1jhMMYY45X/B7JEm99STz85AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "magnetizations_gd = np.real(result_gd.observables)[:, 1]\n",
    "\n",
    "plt.plot(result.times, magnetizations, \"-\", color=\"crimson\", label=\"BFGS\")\n",
    "plt.plot(result_gd.times, magnetizations_gd, \"--\", color=\"crimson\", label=\"Gradient descent\")\n",
    "plt.plot(exact_times, exact_magnetizations, \":\", color=\"k\", label=\"Exact\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.xlabel(\"time $t$\")\n",
    "plt.ylabel(r\"magnetization $\\langle Z_1 + Z_2 \\rangle$\")\n",
    "plt.title(\"Magnetization over time\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The magnetization, however, is computed very precisely."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "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.23.0.dev0+1b4fed3</td></tr><tr><td><code>qiskit-aer</code></td><td>0.11.1</td></tr><tr><td><code>qiskit-ignis</code></td><td>0.7.1</td></tr><tr><td><code>qiskit-ibmq-provider</code></td><td>0.19.2</td></tr><tr><td><code>qiskit-nature</code></td><td>0.5.0</td></tr><tr><td><code>qiskit-optimization</code></td><td>0.5.0</td></tr><tr><td><code>qiskit-machine-learning</code></td><td>0.6.0</td></tr><tr><th>System information</th></tr><tr><td>Python version</td><td>3.10.4</td></tr><tr><td>Python compiler</td><td>Clang 12.0.0 </td></tr><tr><td>Python build</td><td>main, Mar 31 2022 03:38:35</td></tr><tr><td>OS</td><td>Darwin</td></tr><tr><td>CPUs</td><td>4</td></tr><tr><td>Memory (Gb)</td><td>32.0</td></tr><tr><td colspan='2'>Fri Dec 09 16:49:42 2022 CET</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, 2022.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "%qiskit_version_table\n",
    "%qiskit_copyright"
   ]
  },
  {
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}