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} = \\arg\\min_{\\delta\\theta} 1 - |\\langle\\phi(\\theta + \\delta\\theta)|e^{-i\\Delta_t H}|\\phi(\\theta)\\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 `PVQD` in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.algorithms import PVQD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the necessary imports for this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sc\n",
    "\n",
    "from qiskit import BasicAer\n",
    "from qiskit.algorithms.evolvers import EvolutionProblem\n",
    "from qiskit.algorithms.optimizers import SPSA, L_BFGS_B, GradientDescent, COBYLA\n",
    "from qiskit.circuit import ParameterVector, QuantumCircuit\n",
    "from qiskit.opflow import X, Z, I, MatrixExpectation, Gradient\n",
    "from qiskit.quantum_info import Statevector\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "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": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_time = 1\n",
    "hamiltonian = 0.1 * (Z ^ Z) + (I ^ X) + (X ^ I)\n",
    "observable = Z ^ Z"
   ]
  },
  {
   "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": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAB7CAYAAAB9yZWUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAASPUlEQVR4nO3de1xUdf7H8dcMF1HRikAxRYXkoiC6kOYld3TRRNtsa9Oy1fLyWF0ys4utD9dLj+SRFbGPh+6vn5ZdZO1CBdYvMnX3kcp4gTXKC94I4uIshYoKJkjAzJzfH+jokdtgOOcAn+fjMY+HfM8Zzmfkw5vzPefMGYOiKApCCHGZUesChBD6IqEghFCRUBBCqEgoCCFUJBSEECoSCkIIFQkFIYSKhIIQQkVCQQihIqEghFCRUBBCqEgoCCFUJBSEECoSCkIIFQkFIYSKhIIQQkVCQQihIqEghFCRUBBCqEgoCCFUJBSEECoSCkIIFQkFIYSKhIIQQkVCQQihIqEghFBx17oAvdu/YiPnjxVqsm2f8EDujp+tybaF87TqkZvVHxIKzTh/rJDTmce1LkPoWHvrEZk+CCFUJBSEECoSCkIIFQkFIYSKHGhsBbGbX8IvOgS71Ypis1NhOcPhtZ9xckum1qUJnWhLPSKh0EoOr0kle81mDG5GBs6ZhGndIj4/WsjFolNalyZ0oq30iEwfWplis5P7wdcYPdzxCe+vdTlCh/TeIxIKrczo4U7oExMB+LmgRONqhB7pvUdk+tBKIp9+iIi/TMHd2wul1sa+59ZRduIkAN36+2N661m2/n4Z9lor4XFT8PDuzKHXP9G4atdRFIWaC5UoVhudfLphMHa8v0dtpUd0/ZOx2+0kJiYSHByMl5cXQ4YMwWw2Exoayrx587QuTyX7H5/xUdgTfBw+h+KdB/EfHeFYdrHoFCe/2s/ghQ/iHdCDwD+MJnvtZg2rdR1FUchL3kFazGKSB87i48FzSbnrL2Sv3Yz1UrXW5blUW+kRXYfC3LlziY+PZ/78+Wzbto1p06Yxffp0CgoKiI6O1rq8BtVcqGTf8+vpExNFwMRhjvGj674gYEI0pvXP8M3KJOw1Vg2rdA1FUch4fj37nltP2fcWx/ilU2UceDWZ7Q+/SG1llYYVakPvPaLbUEhOTiYpKYm0tDQWL17MuHHjWLZsGSNHjsRqtRIVFaV1iY2qKa/g2IYtRC19DAwGABSrjVP/OY7nLV05s/+ExhW6Rn6KmbzknXVf2JWrC5S6f589+AMHXvlIg8q0p+ce0W0orF69mtjYWEwmk2p8wIABeHh4EBkZCUBRUREmk4mQkBAGDx7Mnj17tCi3nhNvf0WXHrcyYGpd/beG9KHnsDB+2nOEkBnjNa7ONY6/8xUYDU2uk5e8k9qKjre3APrtEYOiKErzq7lWcXExAQEBvPvuu8yZM0e1bPr06eTk5HDw4EEAJk6cyAMPPMCTTz5JRkYGU6dOpbCwEE9Pzya3YTA03axXLPH5LWGefjf2Qq5ujEmfr+KblRu5kP8T9335Mv+atopfzl5o8mk5NaW8dn73r9u2RjoZ3Hiz5x+cWvfVc2a+rz17cwu6ibTqkZb2h7O/6rrcUyguLgbA399fNV5VVYXZbHZMHc6ePcvevXuZO3cuAKNGjeKOO+5g165dri24GWFPTORcdj7nsguwVv7CwYSPGb6qfd8nwYhzoQtgdDKg2zM99YguQ8HX1xeA3Nxc1XhCQgIlJSWOg4wWi4WePXvSqVMnxzqBgYGcPHmy2W0oiuLU4/rpy43ISdrONyuTHF9btmex+8k1zT7PZDI5XafeHpW2Grz79nDMlxtjcHfjux/zNK/31zy06pGW9oezdHmdQlBQEJGRkaxevRofHx969+5NamoqW7duBdDtmQdxlcFgYODsSWS99M/GVzIaCLx/JJ19b3FdYaJZutxTMBqNpKSkEB4eTlxcHLNnz8bX15cFCxbg5ubmOMjYt29fTp8+TXX11fPdhYWF9OvXT6vSxTXCZsfSc+SghhcaDXT1v527Vj7u2qJEs3QZCgAhISHs2rWLyspKLBYL8fHxHDlyhEGDBtG5c2egbpoxevRo3n33XQAyMjL48ccfGTdunJali8vcOnkw4cNlRCx4AI9uXRzjRg937nxoDPdtfYUu/j4aVigaosvpQ2O+/fZbRowYoRp78803mTVrFmvWrMHT05Pk5ORmzzwI13Hv3Im7ls9k6PPT+CDoTwBMO/Q2Xj7dNK5MNEa3ewrXq6ioIDc3t95FS0FBQezevZvc3FyOHj3aKgd9nOHdxw//0RF49/FjzBtPq5b5jwzn4az19J8yqq7Gh8YwOe1lYjYtxcO7bi8nZtNSJn0R75Ja9cC989WDwR0lEFrSI7GfryL2s5e499MX8bq9O6Bdj7SZUPD29sZms7Fw4UKtSwHAO6AHve6JaHR5fqqZorQMDO5uhD5+L9seXEF+qpmQmRMA2PH4K64qVWjE2R4B+NfUl9j+0Ivkp6Rz57SxgHY90qamD3oSMmM8PYaH4RcdipunO2PfWYx3gB87Z72mWq97UC/KTpxEsdkp2ZPNqMQ4jmlUs3AtZ3sE6i5xhro9qvLv/+vqUlXazJ6C3uR+8DX5qWYynluHexcv0v/8d46/tYV+96mPeXh270rtxbrLeGt+voRn9y4NfTvRDjnbIwBde/sy+cuXCZsdS1mOpYHv5joSCq2gPK8YFIVLJefx7N5Vtaz250o8utUdR/Do1oWany9pUaLQWFM9AlD541m23r+MQ4mfEBE3RYMKr5JQuEF2qxWD2+X/vmuuFrv+PRUXCkq4NawvBqORO8YMpvQ79VWaov1ytkcM7m6Of9dcrMJWVeOS+hojxxRuUFnOf4n625/o1rcn9svzwYYoVhu5H37NpP+Lp/pCBbufXOvCKoWWnO2RLj1vY8z/LESxK9hratm76H9dWGV9Ego3qPbiJbY/uFI1dirzGKcyj+EXFUzA+GjKTlgoSsugIHU3Banqd7PFbFpK1elyF1YsXK0lPbL9oRfrPV+rHpFQuAlKD+SRNuGFJteRU5Idm557REKhGT7hgR1y28J5Wv2cbtZ2dXmTFdE+JfV6GIBZJakaVyKaImcfhBAqEgpCNzZs2MDYsWMZO3YsJpMJT09P1q5dW2+ssrJS9bzExEQOHjxIUVERO3fW3Sg2PT2d5cuXq9ZLSkoiNDQUs9nM/v37GTVqFPfccw/PPvusY50ffviBoUOHsnz5cg4fPkxCQsLNf+E6I6EgdGPevHmkp6eTnp7O1KlTWbJkCYsWLao31rXr1Yt/7HY7+/bt4ze/+Y0qFBrzwgsvYDKZ6NevHzt37mTv3r2cOXOGI0eOAHU3Bl6zZg0AQ4YMITMzs0V3LWoPJBSE7hQWFvLhhx+yYsWKJscADh8+zIABA4C6PY3333+fmJgYAL777jvuv/9+Ro8eTUVFhep5/v7+eHl5AeDh4YGbmxsNCQ4OdtwkuKOQUBC6oigK8+fP54033nDcF6OhsSvy8vLo378/ULenMXPmTHbs2AGAp6cnX375JZMnT3aMXS87O5vS0lIGDWr4DlFBQUHk5OS00qtrG+SUpNCV9evXM2zYMNV9OBsac0ZERN3blnv37k15eXm95efPn+epp57i008//VU1tzcSCkI3ioqKeP/99zGbzU2OXSs4OJisrCygbhpgs129nPja9xhcf1zAarUyY8YMEhMT632UwLUKCgp49NFHb+j1tFUyfRC6kZCQQGlpKffee6/jjMOqVavqjVksV99aPGTIEMdHAURERLBv3z4eeeSRZreVkpJCVlYWf/3rXxk7diyZmZkNrpebm8vQoUNb5fW1GYoQLrLR/4/KRv8/tvr3ff3115UDBw40u15KSooSHR2tpKenN7pOXl6ecvfddyuJiYnKoUOHlFdffbU1S20T5IpG4TJyRWPbINMHIYSKhIIQQkVCQQihIqckm7F/xUbOHyvUZNs+4YHcHd++P526PdCqR25Wf0goNOP8sUJOZx7XugyhY+2tR2T6IIRQkVAQQqhIKAghVOSYQiuI3fwSftEh2K1WFJudCssZDq/9jJNbGr50VnQ8balHJBRayeE1qWSv2YzBzcjAOZMwrVvE50cLuVh0SuvShE60lR6R6UMrU2x2cj/4GqOHOz7h/bUuR+iQ3ntEQqGVGT3cCX1iIgA/F5RoXI3QI733iEwfWknk0w8R8ZcpuHt7odTa2PfcOspOnASgW39/TG89y9bfL8NeayU8bgoe3p059PonGlctXKmt9Iiu9xTsdjuJiYkEBwfj5eXFkCFDMJvNhIaGMm/ePK3LU8n+x2d8FPYEH4fPoXjnQfxHRziWXSw6xcmv9jN44YN4B/Qg8A+jyV67WcNqXctWXUv+5qsfm3dswxaqyyuaeEb71FZ6RNehMHfuXOLj45k/fz7btm1j2rRpTJ8+nYKCghbfmstVai5Usu/59fSJiSJg4jDH+NF1XxAwIRrT+mf4ZmUS9hqrhlW6zpmsHD6Nns+ep/7hGMt6MYlPhv6ZHz5N164wDem9R3QbCsnJySQlJZGWlsbixYsZN24cy5YtY+TIkVitVqKiorQusVE15RUc27CFqKWPweVbgilWG6f+cxzPW7pyZv8JjSt0jQs//Mi/H42npuxivWX2Git7F72B5V9ZGlSmPT33iG5DYfXq1cTGxmIymVTjAwYMwMPDg8jISABWrlxJSEgIRqOR1FT93LzjxNtf0aXHrQyYWlf/rSF96DksjJ/2HCFkxniNq3ONo+vTsP5Sg2Jv4D4+igIGAwdf+7jDfa7CFXrtEV0eaCwuLubo0aOqT+65wmKxEB4eTqdOnQCIjY1l1qxZzJkzx9VlOmz/Y/2PEa+tqCJ50OV3sBkMjEyYz3/+9g4X8n/ivi9fxrI9i1/OXnBxpa5jt9rqjiM0FAhXKAplJ05SnmPhtoH9XFecBtpSj+hyT6G4uBig3l12q6qqMJvNqqnDqFGjCAoKavE2DAaDU4/G7iLcEmFPTORcdj7nsguwVv7CwYSPGb6q+be8ms1mp+vU26ObZ2fs1bVO/f/cM3S45vX+modWPdLS/nCWLvcUfH19gbo76U6ePNkxnpCQQElJiW4PMjYmJ2m76mvL9iws29v3XLpasWJVbLgbGv7kpWtV2KtdUJG+6alHdLmnEBQURGRkJKtXr2bTpk3s2LGDuLg43nvvPYBWCQVFUZx6XH9Mw5VMJpPTdertYVXshE4dB8Ym/kIZDNwS3JuimjLN6/01D616pKX94SxdhoLRaCQlJYXw8HDi4uKYPXs2vr6+LFiwADc3N8dBRqFvEU8+gNHDvfFgUBSGLn6kRbu24ubT5fQBICQkhF27dqnGZs6cyaBBg+jcubNGVYmWuG1gP8ZvWsquua9TW1GFwVj3N0ix2zEYjQyPn03glFEaVymup9tQaMi3337LiBEjVGMrVqxg48aNlJaWcuTIEZ555hnMZjN33nmnRlWKa93x20imHniLgs27Kdl7FLvVxu2DAwl+LIauvW7XujzRAF1OHxpSUVFBbm5uvYuW4uPjKS4uprq6mnPnzlFcXOySQPDu44f/6Ai8+/gx5o2nVcv8R4bzcNZ6+l/+K/i7jUt4LOef9Boz2LFOzKalTPoi/qbXqQee3boQNiuWce8sJiZpCUOfn9YhAqElPQLgMziQWSWpGNzqfi216pE2Ewre3t7YbDYWLlyodSkAeAf0oNc9EY0uz081U5SWAUDmkrc4/vZXquU7Hn/lptYntNeSHoG605Jns/MdX2vVI21q+qAnITPG02N4GH7Robh5ujP2ncV4B/ixc9Zr9datOlPu+gKF5lrSI7eG9OFSyXm6B/bSoFK1NrOnoDe5H3xNfqqZjOfW4d7Fi/Q//53jb22h330jmn+y6BBa0iOD5v2eExu3aVBlfRIKraA8rxgUhUsl5/Hs3lXrcoQONdUj3QL9qbl4ierz9d84pgUJhRtkt1odB4S45sIQOecurnC2R24b2A/foQOY8NEybhvYl5GvaXuvEDmmcIPKcv5L1N/+RLe+PbFbbU2uOzx+DgETogmYeBffb/o3uR987aIqhZac7RHL1v1Ytu4H6u76nLlkg6tKbJCEwg2qvXiJ7Q+uVI2dyjzGqcxj+EUFEzA+mrITForSMvhmxXt8s+I91boxm5ZSdbrchRULV2tJj1xx7bspteoRCYWboPRAHmkTXmhyHTkl2bHpuUckFJrhEx7YIbctnKfVz+lmbdegdNTb3gghGiRnH4QQKhIKQggVCQUhhIqEghBCRUJBCKEioSCEUJFQEEKoSCgIIVQkFIQQKhIKQggVCQUhhIqEghBCRUJBCKEioSCEUJFQEEKoSCgIIVQkFIQQKhIKQgiV/weiJvkS1j0V/gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 327.252x144.48 with 1 Axes>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "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": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "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": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "backend = BasicAer.get_backend(\"statevector_simulator\")\n",
    "expectation = MatrixExpectation()"
   ]
  },
  {
   "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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "bfgs = L_BFGS_B()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can define p-VQD and execute it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "pvqd = PVQD(\n",
    "    ansatz, \n",
    "    initial_parameters, \n",
    "    timestep=0.01, \n",
    "    optimizer=bfgs, \n",
    "    quantum_instance=backend,\n",
    "    expectation=expectation\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": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "problem = EvolutionProblem(hamiltonian, time=final_time, aux_operators=[hamiltonian, observable])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then run the algorithm!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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 ~10% of the energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Energy over time')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEYCAYAAAB2qXBEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3NElEQVR4nO3dd3yb9bX48c+RbHnEI3biDOJsskMSyIBMApRAmGUV6CDhQin00gu0lHLb3kIp9NfFKC0UaEtpSwOlAwg0EAJlZUH2dvZyEhInjuMlWZJ1fn9IBtVxEjtalnTer5dekR494zyO5aPvFlXFGGOMOVmORAdgjDEmuVkiMcYYExFLJMYYYyJiicQYY0xELJEYY4yJiCUSY4wxEbFEYoxBRJ4Skf9LdBwmOYmNIzHJTkR2AF2BxrDNz6nq7YmJqH0TkZnAzao6KdGxmNSQkegAjImSS1X17VheQEQyVNUfy2tEWzLGbJKPVW2ZlCYiM0Vkvoj8QkQOi8h2EZke9n6hiPxeRPaJyB4ReVBEnGHHLhCRR0XkEHC/iHQSkddEpFpEloT2nx/a/wkRebjZ9WeLyF3HiG1C6BxHQv9OCG2/VkSWNtv3LhGZHXqeFbqfXSKyP1QtlRN6b6qIlIvId0TkE+APzc4zBHgKGC8itSJSFdr+nIg82Owc94jIgdDP5vMicpGIbBKRShH5btg5HSJyr4hsFZFDIvKSiBSfxH+XSVKWSEw6OBPYCHQGfgb8XkQk9N5zgB84FTgdmAbc3OzYbQSrzh4CngDqgG7AjNCjyR+B60XEASAinYHPAbOaBxT6Q/sv4HGgE/AI8C8R6QS8BgwSkQFhh3wx7Dw/AQYCo0Jx9wB+ELZvN6AY6A3cEn5dVd0A3AosUtU8Ve141E/rs3Nkh537t8CXgdHAZOD/RKRvaN9vAJ8HzgZOAQ6Hfk4mXaiqPeyR1A9gB1ALVIU9vhp6byawJWzfXEAJ/qHsCjQAOWHvXw+8G3bsrrD3nIAPGBS27UFgftjrDcD5oee3A3OOEfNXgI+bbVsEzAw9fx74Qej5AKAmFLsQTGT9w44bD2wPPZ8KeIHs4/y8ZobHHNr2HPBg2DncgDP0Oj/0MzszbP9lwOfD7vm8sPe6h35OGYn+3bBHfB7WRmJSxef12G0knzQ9UdX6UGEkj+C39kxg32cFFBzA7rBjw5+XEGxXPNb7ECyVfBmYF/r3l8eI6RRgZ7NtOwmWACBY+ngYeIBgaeSVUOxdCCaUZWExC8Ek16RCVT3HuG5rHVLVps4L7tC/+8PedxP8GUKw5POyiATC3m8kmKj3RBiHSQKWSEw6202wRNJZj90gHd6tsYJgNVgpsCm0rWez/Z8H1orISGAI8MoxzruX4B/gcL2AN0PP5wElIjKKYCmpqZ3lIME/4sNU9Vh/pE/UFTPaXTV3A/+lqguifF6TJKyNxKQtVd0HvAU8LCIFoUbj/iJy9jH2bwT+SbDRPVdEBgM3NNunHFgC/Bn4h6q6jz4TAHOAgSLyRRHJEJFrgaHA66Hz+IC/AT8nWHKaF9oeINhe8WiodIKI9BCRC9pw6/uBUhFxteGY43kKeEhEeofiKRGRy6N0bpMELJGYVPFaqBdS0+PlVh53A+AC1hNsJP47wTr+Y7kdKCRYXfZn4AWCpZpwfwROC73fIlU9BFwCfAs4BNwDXKKqB8N2m0Wwsf5vzUpM3wG2AItFpBp4Gxh03Lv8T/8G1gGfiMjBE+3cCr8EZgNviUgNsJhgJwWTJmxAojEREJGfAt1UdUbYtikEq7h6q33ATBqwEokxbSAig0VkhASNA24CXg57PxO4A/idJRGTLiyRGNM2+QTbSeqAvxLsWfUqfDrYr4pg1dhjiQnPmPizqi1jjDERsRKJMcaYiKTlOJLOnTtrnz59Eh2GMcYklWXLlh1U1ZLm29MykfTp04elS5eeeEdjjDGfEpHmszEAVrVljDEmQpZIjDHGRMQSiTHGmIhYIjHGGBMRSyTGGGMiYonEGGNMRCyRGGOMiUhajiMxJp0cOtJI+QEfeyv8HKkNcMbgbAb0zCRshUVjImKJxJgUoqpUVgfYtMvLkvVulm7wUH7g6MUfe5RkcM6YXK4+N5+CDs4WzmRM61kiMSYFbNrl5VcvVbJjn486d3Ai1qxMYdTALC6dnEef7pn06JJJTpawcLWb95bVM+vNal7/sJZbryri/HG5VkIxJ80SiTEpYO7iWjbt8nLRhDx6dcukT/dMhvXLwpV5dHK4eGIeF0/MY2u5l0dfqOQnfzzEm4tq+d6NnelUaKUT03bW2G5MCli7tYFh/bK447pirpiaz+mDsltMIuH6l7p4/Ftduev6Isp2ernzkf0cqDy6GsyYE7FEYkySc3sCbN3jY3j/rDYf63AIl07O5+ff6EJVTSN3PrqffQctmZi2sURiTJLbsMNLIADD+7U9kTQZ1i+LX9zRhTq3cucj+9m13xfFCE2qs0RiTJJbu60BERgaQSIBGNQ7i0fu7ILPr/zPL/azYUdDlCI0qc4SiTFJbu3WBvp2zyQvJ/KPc/9SF4/f3ZUO2cK3fnmAJevdUYjQpDpLJMYkscaAsn57A8NOon3kWEq7ZPL43d3oUZLBd5+sYHmZJ2rnNqnJEokxSWzHXh/1Ho2ofaQlnQqdPHpXV7p3zuDhWZV4vIGont+kFkskxiSxtVuD7Rgn02PrRPJyHHzzi8XsO+jn+Teqo35+kzoskRiTxNZua6BToZNunWIzkHDUwGwuHN+Bv86rZtseb0yuYZKfJRJjklhwIKIrptObfO2KjuTlOnhkViWBgMbsOiZ5WSIxJklVVPnZX9kYk2qtcIV5Tm67siPrt3t5bX5tTK9lklNcE4mIXCgiG0Vki4jc28L7U0RkuYj4ReTqZu/NEJHNoceMFo6dLSJrYxm/Me3Juqb2kSg3tLfk/DM7MHpwNr99pYqKwzby3fynuCUSEXECTwDTgaHA9SIytNluu4CZwKxmxxYD9wFnAuOA+0SkKOz9KwH7qmTSyqrNDWRnCaf2dMX8WiLCXV8sJhCAx148jKpVcZnPxLNEMg7YoqrbVNULvAhcHr6Dqu5Q1dVA876GFwDzVLVSVQ8D84ALAUQkD/gm8GCsb8CY9mRZmYeRp2aR4YzP9O+ndM7gxksLWbQmOA29MU3imUh6ALvDXpeHtkV67I+Ah4Hj/maLyC0islREllZUVLTyssa0Twcq/ZQf8HPG4Oy4XvfKc/IZ1NvFr146zJHaxrhe27RfSd3YLiKjgP6q+vKJ9lXVZ1R1jKqOKSkpiX1wxsTQ8o3B0eZnDIpvInE6hLu/VEx1fYAX3rKxJSYonolkD9Az7HVpaFskx44HxojIDmA+MFBE3os4UmPaueVlHoryHfQ9JTPu1+5f6mLKqFzmLKjF3WAj3k18E8kSYICI9BURF3AdMLuVx84FpolIUaiRfRowV1V/o6qnqGofYBKwSVWnxiB2Y9oNVWX5Rg+nD8rG4UjM8rhXnpNPrVt566O6hFzftC9xSySq6gduJ5gUNgAvqeo6EXlARC4DEJGxIlIOXAM8LSLrQsdWEmwLWRJ6PBDaZkza2bHPR2V1IO7VWuGG9XMxsJeLl9+tsUGKJr5rtqvqHGBOs20/CHu+hGC1VUvHPgs8e5xz7wCGRyVQY9qxptl4493QHk5EuPKcfH7yx0MsK/MwdmhOwmIxiZfUje3GpKPlGxvoUZJBt05x/R54lKln5FJU4OCf79YkNA6TeJZIjEki/kZl1WZPQqu1mrgyhUsn5fHROg+7bWnetGaJxJgkUrbDS71HE1qtFe6yyflkOOGV961Uks4skRiTRJaVeRCBUQNjP79WaxQXOpl6Ri5zF9dR77GuwOnKEokxSWTRGjdD+rgozIvN+iMn44qp+dR7lLmLrStwurJEYkySqDjsZ9MuLxNGtK8eUkP6ZjGot4tX3reuwOnKEokxSWLhajcAE0fkJjiSo10xNZ/d+/0sC3VNNunFEokxSWLhGjc9SjLo1S2x3X5bMvWMXIryHbz8njW6pyNLJMYkgTp3gBUbPUwYkRPTZXVPlitTuDjUFXjvQVv4Kt1YIjEmCSxZ78bfCBNHtq/2kXCXTs7DIfC6LcebdiyRGJMEFqx2U9DBwbA4LKt7sko6ZnDG4Gw+WFFvKyimGUskxrRz/kblo7Vuxp+WgzNBs/221uRRueyt8LNtj410TyeWSIxp51ZvaaDWre2u229Lgm048OFKW4o3nVgiMaadW7zGTWYGjBnSPqZFOZ7iAien9c/iw5XuRIdi4sgSiTHt3PIyD6f1zyInKzk+rpNH5bB9r88mckwjyfGbaUyaOlzTyLa9Pk5vB7P9ttbkUcEBk1a9lT4skRjTjq3aFFrEKokSSZfiDAb1dvHhCqveSheWSIxpx5ZvbCA3WxjYy5XoUNpkyum5bNzlZX+lDU5MB5ZIjGnHVmz0MOLULJzO9t3tt7lJoYGTVr2VHiyRGNNO7a/0s6fC324WsWqLnl0z6dM9k4WrrHorHVgiMaadWhlqHzl9YPIlEoCJI3JYvbWB6rrGRIdiYswSiTHt1PIyD4V5DvqekpnoUE7KhJE5BALw0VqbWj7VWSIxph1SVVZsauD0gdk42vm0KMcyqJeL4gIHC9dY9Vaqs0RiTDtUfsDPwapGTh/UfidpPBGHQxh/Wg4fr3Pj9dkkjqnMEokx7dDyjaH2kSQaP9KSiSNycTcoqzZb9VYqs0RiTDu0dIOHLkVOepS0v9UQ2+L0QVlku+TTZYJNarJEYkw74/UpyzZ4OGt4+1wNsS2yXA7GDMlm4Wq3rVGSwiyRGNPOrNzkweNVzjqt/U8b3xoTRuRQUdXI5t02iWOqskRiTDuzaI2bbJck1fxax3PW8BwcAgtX2yj3VGWJxJh2RFVZtMbN6MHZuDKTu1qrScd8J4P7uFhWZg3uqcoSiTHtyLY9Pg4cbmR8ilRrNRk1MJuyHV7cDYFEh2JiIK6JREQuFJGNIrJFRO5t4f0pIrJcRPwicnWz92aIyObQY0bY9jdFZJWIrBORp0TEGY97MSYWmgbvnTU8tRLJyAFZNAZg3baGRIdiYiBuiST0B/4JYDowFLheRIY2220XMBOY1ezYYuA+4ExgHHCfiBSF3v6Cqo4EhgMlwDWxugdjYm3xGjeDersoLkyt70PD+2XhdMDKTZZIUlE8SyTjgC2quk1VvcCLwOXhO6jqDlVdDTQv/14AzFPVSlU9DMwDLgwdUx3aJwNwAdbH0CSlyupGynZ6mZBi1VoAOdkOBvV2fToRpUkt8UwkPYDdYa/LQ9siPlZE5gIHgBrg7y2dQERuEZGlIrK0oqKiLXEbExcfrXWjSsq1jzQZNTCbjTu9uD3WTpJqUqKxXVUvALoDWcC5x9jnGVUdo6pjSkpK4hqfMa2xeK2bko5O+pcm52y/J9LUTrLW2klSTjwTyR6gZ9jr0tC2qByrqh7gVZpVlxmTDPyNyvIyD2OHZSf9aPZjaWonWbXZEkmqiWciWQIMEJG+IuICrgNmt/LYucA0ESkKNbJPA+aKSJ6IdAcQkQzgYqAsBrEbE1NlO7zUeZSxQ1KzWgusnSSVxS2RqKofuJ1gUtgAvKSq60TkARG5DEBExopIOcGeV0+LyLrQsZXAjwgmoyXAA6FtHYDZIrIaWEmwneSpeN2TMdGyZL0bh5CUy+q2hbWTpKa4Ti2qqnOAOc22/SDs+RKC1VYtHfss8GyzbfuBsdGP1Jj4WrLBw+A+LvJzU6LZ8phGDshi1txgO8nYoalb+ko3qf1ba0wSOFLbyMadXsYMSe3SCFg7SaqyRGJMgi0v86BKWnxDb2onWb3FEkkqsURiTIIt3eAhL0cY3NuV6FDiYli/LDbubLDld1OIJRJjEkhVWbLBwxmDs3E6U7Pbb3PD+2fh88OmXd5Eh2KixBKJMQm0Y5+Pg1WNaVGt1WRYvywA1m616q1UYYnEmARauiE4piIdGtqbFBcE16K3Ee6pwxKJMQn08ToPvbpm0LU4rj3xE254/yzWbWuwddxThCUSYxKkuq6RFZs8TByZm+hQ4m5YvyyO1AbYfcCf6FBMFFgiMSZBFq1xEwjAlNPTp32kyfD+wXaSddZOkhIskRiTIB+scNOl2MnAXunR7Tdcr64Z5Oc6rME9RVgiMSYB6j0Blm5wM3lUbsrO9ns8DocwrJ/LGtxThCUSYxLg43VufH6YPCr9qrWaDO+fxe79fo7UNiY6FBMhSyTGJMAHK9wU5Ts+HVORjpraSaxUkvwskRgTZ16fsnidm0kjc3E60q9aq8mgXi4ynLDW5t1KepZIjImzpRvceBqUSWlcrQWQ5XIwoKeL9dttqpRkd8JEIiLD4hGIMeniw5Vu8nKEUQPTZzT7sQzu42Lzbi+NjTYwMZm1pkTy56YnInJz+Bsikn4jqYyJgLshwAcr6pk8KpfMjPSt1moyuHcWHq+y8xNfokMxEWjNvAzhv+1fB34X9vpDYHRUIzImhX2woh53gzJ9Ql6iQ2kXBvUJjqHZuNNLvx4nP56mtj7Apl1eDh1ppLK6Ea9P6dYpgx5dMujVNZO8FF95MtFak0jCy5zNv0LZ/44xbfDGwjpKu2QwrF/6DUJsSWlJBh2yhbIdXqZPaNuxNfUB3ltWx4cr3azc5MF/jF7ETgeMPy2HSyblMXpIdlp3cIiV1iSSbiIyE1jF0YnEKjaNaaU9B3ys3tLAzZcXpuUgxJY4HMLA3i7Kdra+51a9J8A/3q3hpberqXMrPUoyuOqcfMYMzaGkyEmnAieZGcK+g372VAR/5m8trmP+KjfdO2dw53VFaTVtfzy0JpHcT7D66kagVETWAxuAMqBz7EIzJrW8ubgOh8C0MzskOpR2ZXDvLF56uxqvT3FlHjvBqipvLKzjd69WUVUbYOKIHL5yUSEDema2mJh7d8+kd/dMJozI5abLOrJgtZvnXqviO7+u4NJJeXztyo7kZlulSjScMJGo6jPhr0WkFDgNGAF8EKO4jEkpjQFl7uI6xg7NpnPH9Joy/kQG9XbRGICt5V6G9G15gOaR2kZ+8XwlC1a7GXFqFg9d0fGY+7YkM0OYekYu44dn8+xrR/j7v2tYusHN//vvLvTqlhmtW0lbbU7Hqlquqm+o6k9V9cuxCMqYVLO8zMPBqkYuHG+N7M01rVVftrPl8STLN3q46cF9fLzezW1XdeSRO7u0KYmEy3I5uO2qIh67qwser3LXo/vZvtfGsUTKynXGxMGbi+oo6OBg/GlWN99cSZGTogIHG1tIJG8uquU7vzpAfq6DJ+/pxjXnFeCIQmP5aadm8+hdXXE4hLsePWDrx0fIEokxMdbgDbBojZuzz8g9bhtAuhIRBvVysTGswV1V+dOcI/zsz5WMGpjNr77djf6l0e3p1qtbJo99swvZWcLdv9zPjn02luVktTqRiMg3RKQolsEYk4qWb2zA41UmjbTSyLEM6ZPFrv1+6twBVJXHXjzMc68f4YKzOvDjr5eQlxOb77w9SjJ57K6uZGYI9z1TQZ07EJPrpLq2/O90BZaIyEsicqFY/0VjWmXBqno6ZNuUKMczqLcLVdi0y8uv/3aY1z6s5bppBdzzleKYzwDQrVMG/3dTZ/ZU+Pn584dsHfmT0OpEoqrfBwYAvwdmAptF5Mci0j9GsRmT9BoDyqI1bsYNz7EpUY5jUKjB/ZFZlbz8Xi3XnJfPV+M43mbUwGy+enlHPljh5m/v1MTlmqmkTeVFDabqT0IPP1AE/F1EfhaD2IxJehu2ezlcE2DSCKvWOp7CPCfdOznZU+Hn8il53Hplx7gP2vzC5/KZcnoOz7xSxYbtNrV9W7SljeQOEVkG/AxYAJymqrcRHKx4VYziMyapzV9VT4YTxg2zRHIiV51bwHXTCvjGF4oSMvJfRPj2lztRXODkkRcqbUbiNmhLiaQYuFJVL1DVv6mqD0BVA8AlMYnOmCSmqixY5eb0Qdl0iFFjcSq58px8bvl8x6h07z1ZHXIc3H5NEVvLffzzPaviaq22tJHcp6o7j/HehtacI9RIv1FEtojIvS28P0VElouIX0SubvbeDBHZHHrMCG3LFZF/iUiZiKwTkZ+09n6MibWdn/jZU+FnolVrJZXJo3I4a3g2f3j9CPsr/YkOJym0pWrrmy08bhKRUa083gk8AUwHhgLXi8jQZrvtItiQP6vZscXAfcCZwDjgvrCuyL9Q1cHA6cBEEZne2nsyJpYWrKoHYIIlkqQiIvzPtcVoAH790uFEh5MU2lLeHgPcCvQIPb4GXAj8VkTuacXx44AtqrpNVb3Ai8Dl4Tuo6g5VXQ0078x9ATBPVStV9TAwD7hQVetV9d3QsV5gOVDahnsyJmbmr3QzqLfL5tZKQt06ZTDj4kIWrHbz8Tp3osNp99qSSEqBM1T1W6r6LYKN7F2AKQRLESfSA9gd9ro8tK01TnisiHQELgXeaekEInKLiCwVkaUVFRWtvKwxJ2fXfh8bd3k5Z7QtIpqsrjo3n26dnPzhtSM2tuQE2pJIugDhfeJ8QFdVdTfbHncikgG8ADyuqtta2kdVn1HVMao6pqSkJL4BmrQz76PglPHnjbUp45NVZobwlemFbNzlZdEaK5UcT1sSyV+Aj0TkPhG5j2AX4Fki0gFY34rj9wA9w16Xhra1xomOfQbYrKqPtfJ8xsRMIKDM+7iO0UOy6VToTHQ4JgLnn9mBU0oyeO51K5UcT6sSSWg6lOeAW4Cq0ONWVX1AVetU9UutOM0SYICI9BURF3AdMLuVcc4FpolIUaiRfVpoGyLyIFAI3NnKcxkTU2u2NHCgstEWsEoBGU7hhukFbCn3MX+VlUqOpVWJJDSifY6qLlXVX4YeS9tyIVX1A7cTTAAbgJdUdZ2IPCAilwGIyFgRKQeuAZ4WkXWhYyuBHxFMRkuAB1S1MrTI1vcI9gJbLiIrReTmtsRlTLS99VEdOVnCRJukMSWcN7YDpV0y+OO/jhAIWKmkJW3pTrJcRMaq6pKTvZiqzgHmNNv2g7DnSzhGrytVfRZ4ttm2co5eR96YhPF4A7y/op4pp+eS7bJBiKnA6RRmXFzIQ384xILVbiaPsg4UzbXlN/1MYLGIbBWR1SKyRkRWxyowY5LRwtVu6j3K+VatlVKmjs6la7GTl9+10e4taUuJ5IKYRWFMinjrozq6FDkZNeDkloI17ZPTIVw+JZ9nXqli2x4v/XpEd5GtZNeWEskuYDIwIzRVihJco8QYA+yp8LFkvYdpZ3VI6HxRJjYumtiBrEzhZZuD6yhtSSRPAuOB60OvawhOeWKMAV5+rxaHwGWT8xIdiomBgg5OzhuXy9sf11Nd15jocNqVNrWRqOp/Ax6A0FQlVr4zBqhzB3hzUS3njM61KVFS2BVn59PgU95YWJfoUNqVtiQSX2jiRQUQkRKOnhPLmLT0xqJa6j3KVefmJzoUE0P9S12MHJDFqx/U0GhdgT/VlkTyOPAy0FVEHgLmAz+OSVTGJJHGgPLyuzUM75/FoN7WyJ7qrpiazyeHGlm81gYoNml1GVxV/xJaIfG80KbPt3YdEmNS2aLVbvYdauSWK4pOvLNJehNH5FBU4OCtxXVMHGFjSqBt65FkAWcQnI6kE3CNiPzg+EcZk/r+8W4NXYudTLKR7GnB6RTOG9OBRWvc1uge0paqrVcJrh/iB+rCHsakrZ37fKza3MBlU/JxOq3Lb7o4/8wO+BvhvWX1iQ6lXWhL95JSVb0wZpEYk4TmLKzF6YALzrKR7Onk1NJM+p6SyVsf1XHZFOtg0ZYSyUIROS1mkRiTZLw+5a2P6pg4MofiApsuPp2ICOeP68D67V7KD/gSHU7CtSWRTCI4ceNGm2vLGFi4up4jtQEummADENPReeNyEYG3P7Ya/rZUbU2PWRTGJKF/LaijS7GT0UOyEx2KSYCSjhmcMSibeR/VMePiQoLLNqUnm2vLmJOw76CfZWUepo/Pw2nzaqWt88/swL5DjazZmtDVxhPO5toy5iTMWRicV2v6eGtkT2eTR+bgypS0771lc20Z00Y+v/LmojrGDs2mS7HNq5XOcrIdjB2SzYJV7rRe093m2jKmjd7+uI5DRxq5Yqp1+zQwaVQOFVWNbNzpTXQoCXMyc211sbm2TLpqDCgvvFXNqT0zGTvUGtkNjD8tB4cD5q9K37m3Wp1IVPUvwD3A/wP2EZxr62+xCsyY9ujDlW7KD/j54gXp3UvHfKagQ3BFzA9Xpm87SZsqeFW1DCiLUSzGtGuqyqy5RyjtksHkUTavlvnMpFG5PP7Xw+zc56N398xEhxN3banaMiatLVnvYctuH9dPK7Auv+Y/TAxN2Dl/VXqWSiyRGNMKwdJINSUdnXxunHX5Nf+ppGMGg/u4mL8yPdtJLJEY0wqL1rhZvaWBa88vIDPDSiPmaJNH5bJxl5f9lf5EhxJ3lkiMOQGvT3ni71X07pbBZVNsXi3Tsqb1aBakYe8tSyTGnMDf3qlm30E//31NERm25og5hp5dMyntksGS9ZZIjDFhKg77+cub1UwamcOYIdZTyxzfmCHZrNrcgNeXXqPcLZEYcxxPv1xFQOHrV9t67ObExg7NwePVtJvE0RKJMcewfa+Xfy+t5wufy6dbJ5tTy5zYqAFZZGaQdtVblkiMOYZ/vltDVqZw1Tk2p5ZpnZxsB8P7Z7FkvSfRocSVJRJjWnCktpF5H9dz/pkdKMyzZXRN640dksP2vT4qqtKnG3BcE4mIXBhaqneLiNzbwvtTRGS5iPhF5Opm780Qkc2hx4yw7Q+JyG4RqY3HPZj08K/5tXh9yhVTrbuvaZumyTyXbkifUkncEkloCvonCC7ZOxS4XkSGNtttFzATmNXs2GLgPuBMYBxwn4g0tX6+FtpmTFT4G5VXPqhl9OBs+p5iS+6YtunXI5NOhU6WrLNEEgvjgC2quk1VvcCLwOXhO6jqDlVdzdHrnFwAzFPVytCCWvOAC0PHLFbVfbEP36SLD1bUc7CqkavOtbYR03Yiwpgh2Swr89AYSI9uwPFMJD2A3WGvy0PbYn0sACJyi4gsFZGlFRUVbTnUpJl//LuG0i4ZjLP1RsxJGjs0m5r6AJvSZLGrtGlsV9VnVHWMqo4pKSmJ2XWqahpZvtHDpl1e9lT48PnT4xtJqli12cOGHV6uPCcfh83wa07S6MHZiMCSNGkniWfn+D1Az7DXpaFtrT12arNj34tKVFH28F8qWbD6sz7k/UszeeLb3XBl2h+lZPD8G9UUFTiYPt5m+DUnrzDPSf8emazc5OGGiwoTHU7MxbNEsgQYICJ9RcQFXAfMbuWxc4FpIlIUamSfFtrW7mzf5+OMQVn86NbO3Hx5IVvLffz9nepEh2VaYe3WBpaVebj2cwVkudKmsG5iZMSAbNZv96bFdClx+7Soqh+4nWAC2AC8pKrrROQBEbkMQETGikg5cA3wtIisCx1bCfyIYDJaAjwQ2oaI/Cx0TK6IlIvI/fG6p+Z8fuWTg36G9sti4ohcvnhBIZNH5fDnN6rTcmrpZPP8G0cozHNw6WTr8msiN2pAFl6fsnFn6k+XEtevXao6R1UHqmp/VX0otO0Hqjo79HyJqpaqagdV7aSqw8KOfVZVTw09/hC2/Z7QMY7Qv/fH857C7TvkJ6BQ2uWzpTab5mh68u+HExWWaYWyHQ18vN7DNeflk5NlpRETudNOzQJg5WZLJKYN9hwIljpKu3zW9NS1OIOvTC/gw5VuPl6XXvPvJJM/v1FNQQcHnz/buvya6CjMc9LvlExWWyIxbVF+wAdAj5L/7MNw9XkFlHbJ4Kl/VqGa+vWlyWbjzgYWrXFz1Tn55GbbR8JEz4gBWazb1oC/MbU/9/apiaLyA37ycx1Hzc3kyhSuOjefHft87PrE2kram9+9GmwbudImZzRRNmpgNh6vsjHFx5NYIomiPQd8/1GtFW7CaaFlOFfXxzMkcwLLyjwsK/PwpQsL6JBjHwcTXSNC7SSrNqX2eBL75ERReYWfHsdIJCVFGQzs5WLhamsnaS9Uld++UkWXYieXTbbSiIm+jvlOenfPZNWW1G4nsUQSJQ3eABWHG/+jx1ZzE0fksGGHl8rqxjhGZo7lgxVuNu3yMvPiQhswamJm1IAs1m5N7XYSSyRRsvegH9WjG9rDTRiRgyosWmOlkkTzNyq/n11F7+6ZnH+mjWI3sTNyQBbuBmXzrtRtJ7FEEiXloa6/Pbseu0TSr0cmXYudLFhl7SSJ9tK8asoP+Pnq5YU4bU4tE0MjTg1O/rkqhbsBWyKJkqYxJMcrkYgIE0fksHxjA+6G5jPlm3jZvtfLH+ccYcrpOUwYkZvocEyKKy50Utolg3XbLJGYEyiv8FGU7zhhz58JI3Px+pRlaTIraHvT2Kj8/PlKcrIc3HFdcaLDMWliaN8s1m9vSNlxZJZIomTPAT89jtPQ3mTEqVnk5ch/zBBs4udv79RQtsPLHdcWUZRva7Gb+Bja18XhmgCfHErNjjaWSKKk/ID/mGNIwmU4hbOG5zB/ZT279/viEJlpsueAjz+8XsWkkTlMHW1VWiZ+hvYNjidZvz01q7cskUSB2xPg0JFGSo/TPhJuxiWFZGQI//tkBVU1qfkNpT36/ewjOJ3CHdcVI2IN7CZ++p6SSXaWpGw7iSWSKNhTEWpob0XVFkCPkkwevLWEg1WN/N/TFWmxXkGile1o4L3l9XzhvHw6FVqVlokvp1MY3NvFhu2p2QXYEkkUlFccPevviQzrl8V3Z3Zi/XYvP/3zoZRthGsPVJVnXqmiY56Da84rSHQ4Jk0N7ZvFlnIvDd7U67FpiSQKykNtHae0smqryZTTc7nxkkLeXVrPR2utF1esLFnvYeWmBr48vdDm0zIJM7Svi8YAbErBgYn2qYqCPRV+Ond0ntSCSNdNK6Bn1wx+88/DKT2FQqIEAsHSSPfOGbbyoUmopgb3dSlYvWWJJAq27/XRs2vbSiNNMpzCrVcWsXu/n9kf1EY5MvP6/Fq27fFx02WFZGZYA7tJnI75Tk4pyWB9Cja4WyKJUL0nwJbdXob3yzrpc5w1PJvRg7P547+OUF1nvbii5ZNDfp5+uYrRg7M5x7r7mnZgaF9XSg5MtEQSoXXbGggojBiQfdLnEBFuu6ojde4Af5pTHcXo0lcgoPz8+UOIwN1ftu6+pn0Y1jeLyuoA+ytT6wujJZIIrd7cgNMR/KYRiX49XFw0MY9X36/hYJWtohip1+fXsmJjA7deWUTX4pOrdjQm2oaE2kk2pNjAREskEVq9pYGBvVwn1dDe3LXn5xPQ4B9Bc/L2VPg+rdK6eKJNEW/aj/49Msl2Sco1uFsiiUCDN0DZzoZPl9OMVI+STM4cls1r82vx+VOrDjVeauoDfO/JCjIzhG99yaq0TPvidAqn9nSlXBdgSyQRKNvhxeePrH2kuc+fnc/h6gAfrLA1S9rK36jc/9sK9h7088NbOtOtk1VpmfZnYM9Mtuz20hhInS+LlkgisHpLAyIwvH90SiQAY4ZkU9olg5ffq4naOdOBqvLYi5Ws2NjA3V8qZmQUk7sx0TSwlwuPV9m9P3XaQi2RRGD1lgb69cgkPzd6P0aHQ7h8Sh7rt3tTrvgbK40B5cl/VDFnQR1fnl7AtLNs4KFpvwb0CnbMSaXPtyWSk+RvVNZti177SLgLxueRnSW88r6VSk6kwRvggd8d5B//ruHKc/KZeXFhokMy5rh6dQs2uFsiMWze5cXj1U/XY46mvBwH08Z14J0lddTWp94Eb9FSUeXnW788wPxVbr5+dUduv6YIh62/bto5p0PoX5rJZkskZtWWYD/w02JQIgG4YHwHfH6Yv8oa3Zvz+ZUX51Uz84f72FLu476bO3P1uTarr0keA3q62FyeOg3ulkhO0tL1bnp2zaC4IDZrWwzu7aJ7JyfvLbNEEm7nPh+3/Hgfz7xcxcgBWfz++92YcrpNf2KSy8BeLjwNSnmKNLhbIjkJb31Ux/KNDVxwZuwGu4kIZ4/uwLIyD0dqU2s6hZPl9SkP/P4g1XUBHrqthB9/vQs9Slq3mJgx7cnAFGtwj2siEZELRWSjiGwRkXtbeH+KiCwXEb+IXN3svRkisjn0mBG2fbSIrAmd83GJ8Qi08gM+fvliJSNOzeLaabGtTpl6Ri6NAZi/yh3T6ySL516vYvteH9/+SifGn5aT6HCMOWm9u2WSlSls2m2JpE1ExAk8AUwHhgLXi8jQZrvtAmYCs5odWwzcB5wJjAPuE5Gi0Nu/Ab4KDAg9LozRLeDzKw8+e4jMDOF7N3bCGeOG3QE9M+lRkmHVW8DarQ389e0aLp7YgbOGWxIxyc3pDDa4p0qJJJ5Df8cBW1R1G4CIvAhcDqxv2kFVd4Tea95V6QJgnqpWht6fB1woIu8BBaq6OLT9T8DngTdicQO/e7WKTbu8PHBLZ0qKYv+jExGmjs7lhbnVHK5ppCg/Pdcad3sC/ORPh+hW7OS2q4pOfIAxSWBALxdvLa4jEND/6G1Y7wmwZmsD5ft9VNUEOFzbSMc8J+eMzqVfj8x2Oe1PPBNJD2B32OtygiWMkz22R+hR3sL2o4jILcAtAL169WrlZT/jb1Q27fJy+ZQ8Jo2KX+PuOaNz+cub1Xy4op7LpuTH7brtharyy78eZt9BP4/e1YXcbGvWM6lhYC8Xr75fS/kBP907B2ez+GBFPWU7vQRCX6UdDuiY56CqNsCsudX06prB9Al5XDE1H1dm+0koaTMZkao+AzwDMGbMmDb3uctwCr+4owuNcW737ntKJr26Bqu30jGRvPJ+LW99VMeMiwtjMmbHmEQZ2DPY4P7K+zUsL/Owa7+fwX1cXH9+AacPyqZ/aXDWDIdDqKpp5IMV9byztJ6nX67i5fdruPmyjpw7JrddjJ2KZyLZA/QMe10a2tbaY6c2O/a90PbSkzxnmzkdgjPOX4hFhHPGdOBPc47w6AuVfPGCgrRZX2PNFg9P/v0w40/L4SvTbZyISS19umfiyhReeb+W7p0z+PHXS47Z/tcx38llU/K5bEo+KzZ6eOqfh/nxc4d4fX4tD95WQl5OYkvq8bz6EmCAiPQVERdwHTC7lcfOBaaJSFGokX0aMFdV9wHVInJWqLfWDcCrsQg+kb5wXj6XTMrjjYW1fOW+vTz6QmXKL351sMrPD393kO6dM/jfmZ3axbcuY6LJ6RRmXlzIjZcW8uz3u7W6E8npg7L5zXe6cfeXilm3rYF7fnUg4TNgSDzXDhaRi4DHACfwrKo+JCIPAEtVdbaIjAVeBooAD/CJqg4LHftfwHdDp3pIVf8Q2j4GeA7IIdjI/g09wU2NGTNGly5dGu3bi7n9lX5emFvNnIW1OB3C1efmc+20goR/G4k2n1+569H9bNvr44lvd6XvKZGtPmlMqlqwqp4f/u4g/Xu4+Nn/dInqBLItEZFlqjrmqO2ptgh9ayRrImmy96CfP7xWxTtL6ino4GDmJYVcOikPpzM1vrX/8q+VvPp+LT+4uTNTz7BR68Ycz8LVwWTSr4eLx+7qQpYrdsnkWIkktb7KpolTOmfwvRs789S93ejXI5PH/3qYm3/8CUvWJ//AxbcW1/Lq+7Vcc16+JRFjWmHCiFx+cFNnNu3y8vCsShJROLBEksQG9nLx8B1d+NHXOuPzK9/5dQXf+00Few74Eh3aSdmy28sjLxxm1IAsbvl8x0SHY0zSmDgyl5mXFPL2x/X84934Lz9hiSTJiQgTR+by7Pe7c8sVHVm5ycN/PbiP375ShduTPFPQH65p5PtPV1DQwcH3b+qcMtV0xsTLly4oYPKoHJ76ZxXLyzxxvbYlkhThyhSuO7+AP91/CueM7sALb1Vzww/38c6SuoQUddvC51fuf+YgVTUBfvS1zjGbUdmYVOZwCN+5oRM9u2byo2cPUhPHnlyWSFJMp0In987oxK/u7kqnQicP/eEQdz5yoN3O6RMcuV7Jmq0N3POVYgb1js36Lsakg9xsB9+/sRM1dQGee70qbte1RJKihvXL4sl7unL3l4rZvd/HbT/9hF88f4jK6vY1Jf1Lb9cwZ0EdX7qggHPHxG5afmPSRf9SF5dOzuPVD2rZtic+XyAtkaQwh0O4aGIef/rhKVx9bj5zF9dxw/17mTX3CA3exLefzHrzCE+/XMXZZ+Ry46W21rox0TLzkkLychz8+qXDcanatkSSBvJyHNx2VRHP/l93Rg3M5nevHuGG+/fx5qLahCz1qar8/tUqfjf7CJ8bl8v3b7SR68ZEU2Gek/+6tJCVmxt4f0XshwVYIkkjPbtm8uCtJTx6ZxeKC5387M+V3PzgPv69NDiVdTx4vAEenlXJX+ZWc/HEDtx7QyfroWVMDFw8KY/+pZk89Y/DeH2x/XxbIklDIwdm8+Q9Xbnv5s6ICA8+e4ibHvqEeR/V4fPH7hdua7mX237yCXMW1PHFCwr45heLrSRiTIw4HcKtVxZx4HAjbyysjem1bIqUNBcIKO8vr+dPb1Szc5+PToVOrpiax0UT8ugYpYW03J4A/3i3hj+/cYSCDk6+c0MxY4bYKofGxJqqcucjB9h3yM/zPzwl4jVMbK6tMJZIjhYIKEs3ePjbOzUsK/PgdMCZw3OYdmYHzhyWfVLz99R7Arz6QS1/nVdNdV2AyaNyuOv64qglKGPMiS0r8/Dtxw9wx7VFXH52ZGsaHSuRpMfCFuaEHA5h3LAcxg3LYfteL28truPtJfUsXO0mMwNO65/F6CE5DOnjolfXTIoKHEct+enzK58c8rNqcwMLV9ezrMyDzw/jhmUz46JChvS1MSLGxNsZg7IY3j+LWXOrmT4hLyYrK1qJxBxTY0BZuamBj9e5WbrBw/a9n83hlZcjdMx34nQKmU6org9QcbiRpl+n7p2cjB+Ry3ljcxnSxxKIMYn0aankuiIuj2ClVSuRmDZzOoTRg7MZPTi4xG3lkUa27fWy6xM/u/f7qHUH8PkVnx/6niJ075xB984ZDOzlok/3zKNKLMaYxPi0VPJmNdPHR79UYonEtFpxoZPiwhzGDEl0JMaYthARZlxcyJuLanE3BHBlRred0hKJMcakgfDahWizcSTGGGMiYonEGGNMRCyRGGOMiYglEmOMMRGxRGKMMSYilkiMMcZExBKJMcaYiFgiMcYYE5G0nGtLRCqAnSd5eGfgYBTDSQZ2z+kh3e453e4XIr/n3qpa0nxjWiaSSIjI0pYmLUtlds/pId3uOd3uF2J3z1a1ZYwxJiKWSIwxxkTEEknbPZPoABLA7jk9pNs9p9v9Qozu2dpIjDHGRMRKJMYYYyJiicQYY0xELJEcg4hcKCIbRWSLiNzbwvtZIvLX0PsfiUifBIQZNa2432+KyHoRWS0i74hI70TEGU0nuuew/a4SERWRpO8q2pp7FpEvhP6v14nIrHjHGG2t+N3uJSLvisiK0O/3RYmIM1pE5FkROSAia4/xvojI46Gfx2oROSPii6qqPZo9ACewFegHuIBVwNBm+3wdeCr0/Drgr4mOO8b3ew6QG3p+WzLfb2vvObRfPvABsBgYk+i44/D/PABYARSFXndJdNxxuOdngNtCz4cCOxIdd4T3PAU4A1h7jPcvAt4ABDgL+CjSa1qJpGXjgC2quk1VvcCLwOXN9rkc+GPo+d+B80RE4hhjNJ3wflX1XVWtD71cDJTGOcZoa83/McCPgJ8CnngGFyOtueevAk+o6mEAVT0Q5xijrTX3rEBB6HkhsDeO8UWdqn4AVB5nl8uBP2nQYqCjiHSP5JqWSFrWA9gd9ro8tK3FfVTVDxwBOsUluuhrzf2Gu4ngN5pkdsJ7DhX5e6rqv+IZWAy15v95IDBQRBaIyGIRuTBu0cVGa+75fuDLIlIOzAG+EZ/QEqatn/cTyogoHJN2ROTLwBjg7ETHEksi4gAeAWYmOJR4yyBYvTWVYKnzAxE5TVWrEhlUjF0PPKeqD4vIeODPIjJcVQOJDixZWImkZXuAnmGvS0PbWtxHRDIIFokPxSW66GvN/SIinwO+B1ymqg1xii1WTnTP+cBw4D0R2UGwLnl2kje4t+b/uRyYrao+Vd0ObCKYWJJVa+75JuAlAFVdBGQTnNwwVbXq894WlkhatgQYICJ9RcRFsDF9drN9ZgMzQs+vBv6toZasJHTC+xWR04GnCSaRZK83hxPcs6oeUdXOqtpHVfsQbBe6TFWXJibcqGjN7/UrBEsjiEhnglVd2+IYY7S15p53AecBiMgQgomkIq5Rxtds4IZQ762zgCOqui+SE1rVVgtU1S8itwNzCfb6eFZV14nIA8BSVZ0N/J5gEXgLwYat6xIXcWRaeb8/B/KAv4X6FOxS1csSFnSEWnnPKaWV9zwXmCYi64FG4Nuqmqwl7dbe87eA34rIXQQb3mcm8ZdCROQFgl8GOofafe4DMgFU9SmC7UAXAVuAeuDGiK+ZxD8vY4wx7YBVbRljjImIJRJjjDERsURijDEmIpZIjDHGRMQSiTHGmIhYIjHGGBMRSyTGGGMiYonEmCgTkY4i8vWw1wtjdJ1SEbk2Fuc2pi0skRgTfR0JrlcDgKpOiNF1ziO47oQxCWWJxJjo+wnQX0RWisjPRaQWQET6iEiZiDwnIptE5C8i8rnQlO2bRWRc0wlE5Msi8nHoHE+LiDP8AiIyieDsxFeH9ukX1zs0JowlEmOi715gq6qOUtVvN3vvVOBhYHDo8UVgEnA38F34dOLAa4GJqjqK4JxXXwo/iarOJzgh4eWh6yTzxIomydmkjcbE13ZVXQMgIuuAd1RVRWQN0Ce0z3nAaGBJaILMHKClGZcHAWUxj9iYE7BEYkx8ha/jEgh7HeCzz6MAf1TV/z3WSUJTvB8Jrc5pTEJZ1ZYx0VdDcGGsk/UOwbaPLgAiUiwivZvt04ckX1vcpA5LJMZEWWj9jgUislZEfn4Sx68Hvg+8JSKrgXlA92a7lRFcb2KtiMSqV5gxrWLrkRhjjImIlUiMMcZExBKJMcaYiFgiMcYYExFLJMYYYyJiicQYY0xELJEYY4yJiCUSY4wxEfn/8vCxE6H1zAAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "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": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exact(final_time, timestep, hamiltonian, initial_state):\n",
    "    \"\"\"Get the exact values for energy and the observable.\"\"\"\n",
    "    observable = (Z ^ hamiltonian.num_qubits).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": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "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": 14,
   "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": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEYCAYAAABY7FHWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyyUlEQVR4nO3dd3gU1ffH8ffZFHovoQsi+hNFEQNYUEEUKUpREJEOAoJgA0XpTSyIYkFBLCgqiCCKglgoIipKUKRaqNIJPZQQkpzfHzvoJt8Au5DdySbn9Tz7sDszO/uZJORk7p25V1QVY4wx5hSP2wGMMcZkLVYYjDHGpGGFwRhjTBpWGIwxxqRhhcEYY0waVhiMMcakYYXB5Cgi8qWIdAzCfteISN3M3m9WIiJtReRrt3OY4BO7j8GcKxHZDJQByqjqXp/lvwHVgUqqutmVcN4cw4CLVLVdJu93MrBNVQdl5n6zEhGpCGwColQ12eU4JsTsjMGcr01Am1MvRKQakNe9OCZQIhLpdgaTtVhhMOdrCtDB53VH4D3fDUSkiYj8JiKHRWSr85e87/oOIrJFRPaJyGAR2SwitzjrhonIdBF5T0QSnCabWJ/3lhGRmSISLyKbRORBZ3lDYADQWkSOiMjvzvJFInKf8/x3Z92ph55qDhKRj0Vkl4gcEpHFInKZs7w70BZ43HnP585y38y5RGSciOxwHuNEJJezrq6IbBORviKyR0R2ikjn031xneObLSL7RWS9iHTzWX5cRIr6bHuViOwVkSjndRcRWSciB0TkKxG5wGdbFZEHRORv4O8MPnqx8+9B5zivFZFOIrIk3T56icjfzvdmpIhUFpEfne/1dBGJ9tn+dhFZISIHnW2uON1xG5epqj3scU4PYDNwC/AncCkQAWwDLgAUqOhsVxeohvcPkSuA3UBzZ11V4AhQB4gGngdOArc464cBiUBjZ/9PA0uddR5gOTDEee+FwEbgNp/3vp8u8yLgvgyOpTvwB1DQed0FKADkAsYBK3y2nQyMyuhr4TwfASwFSgIlgB+BkT5fi2RnmyjnuI4BRU7zNV4MvAbkxts8Fw/c7KxbAHTz2XYMMMF53gxY73xfIoFBwI8+2yrwDVAUyJPB51Z0ton0WdYJWJJuH58BBYHLgBPAfOf7UAhYC3R0tr0K2APUdr6PHZ2vWS63f47t8b8PO2MwmeHUWcOtwDpgu+9KVV2kqqtUNVVVVwJTgZuc1S2Bz1V1iaom4f0ln77ja4mqzlXVFOezrnSW1wRKqOoIVU1S1Y3AJOCeQMKLSB1gFNBUVQ87md9W1QRVPYG3wFwpIoX83GVbYISq7lHVeGA40N5n/Uln/UlVnYu3MF6SQa7ywPVAf1VNVNUVwJv8d4b2IU4znogI3uP+0Fl3P/C0qq5Tbx/BaKC671mDs36/qh7387gy8pyqHlbVNcBq4GtV3aiqh4Av8RYE8Bbeiar6s6qmqOq7eAvJNefx2SZIrDCYzDAFuBfvX5TvpV8pIrVFZKHT3HMI7y+t4s7qMsDWU9uq6jFgX7pd7PJ5fgzI7bSLXwCUcZomDorIQbzNRzH+Bnd++U7H+5ftX86yCBF5RkQ2iMhhvH/Z4pP5bMoAW3xeb3GWnbJP03boHgPyn2Y/+1U1Id2+yjrPZwLXikhp4EYgFfjeWXcB8JLP12U/ID7vBZ+v+3nY7fP8eAavTx3XBUDfdN+r8qT9upgswjqdzHlT1S0isglvs0jXDDb5EHgVaKSqiSIyjv9+ye7E569lEckDFPPzo7cCm1S1yuminenNzmd9CoxT1S99Vt2LtynmFrxFoRBwAO8v1rPuF9iB9xfhGud1BWdZoHYARUWkgE9xqIBzRqaqB8R7+WhrvE1G01T1VLatwFOq+sEZ9n+m48jsyxVP5Xkqk/drgsDOGExm6Yq37ftoBusK4P3LN1FEauH9xXvKDOAOEbnO6agcxn+/gM/mFyBBRPqLSB7nL/3LRaSms343UFFETvdz/jbwh6o+l0HeE3jPXPLibYbxtRtvO/rpTAUGiUgJESmOt3nsfT+P6V+quhVv/8TTIpLb6aztmm5fH+JtWmrJf81IABOAJ306zQuJSKsAPj4e7xnImY4zEJOA+52zRxGRfOK9KKFAJu3fZCIrDCZTqOoGVY07zepewAgRScD7S3K6z/vWAH2AaXjPHo7g7aQ84cdnpgC34+2U3QTsxdsGf6ov4GPn330i8msGu7gHaJHuyqQb8DaHbcH7l/lavB3Jvt4CqjpNIp9msN9RQBywElgF/OosOxdt8HYE7wBmAUNV9Vuf9bOBKsAuVf391EJVnQU8C0xzmsNWA438/VCnSe8p4AfnOM+rL8D52eiG98zxAN6O8U7ns08TPHaDm8lSRCQ/cBCooqqbXI5jTI5kZwzGdSJyh4jkFZF8eC9XXcV/Hb7GmBCzwmCygmZ4m0p24G0WuUftVNYY11hTkjHGmDTsjMEYY0waYX8fQ/HixbVixYpuxzDGmLCyfPnyvapaIqN1YV8YKlasSFzc6a6SNMYYkxER2XK6ddaUZIwxJg0rDMYYY9KwwmCMMSYNKwzGGGPSsMJgjDEmjZAVBhF525nKcPVp1ouIvOxMX7hSRGqEKpsxxpj/hPKMYTLQ8AzrG+EdDqEK3tmeXg9BJmOMMemE7D4GVV0sIhXPsEkz4D1njJylIlJYREqr6s5g5Dm+dCVTxr1CoqbS4cZbkDy5WbZ7K5GFCnDtTTcQVaE0ngL5gvHRxhiTpWWlG9zKknaqwW3Osv8pDCLSHe9ZBRUqVDinD0uMW80Hn8zkkJ6kyYK/ABh0aBUpqkwrfAUAvRL/plDBQrzasjPRVSuzIa9w8W11KVDS3xkejTEm/GSlwuA3VX0DeAMgNjb2nEYBLNL7Xhb2uoeUxBNEpCqpxxJ5e8XvHN+zn5j8RTj5z04un/oO+Y6c4OjXP5Dw4Rya7V/GFVEFeePG5uSuXY1VJfJQq11L8hUtnJmHZ4wxrspKhWE73snBTynnLAsaj8eDJ28e7/P8eanaoF6a9WMfbAuAqpKyZz+vvPkuuf/Zg+w6yva3Z3DrjsV0HzSKQY1akq/JjZy4rholL6oUzMjGGBN0WakwzAZ6i8g0oDZwKFj9C4ESESJjitFi4KP/Lit2KIEZr71F6X/2cjJuPYsfGkqLgyt467pm3PlIL/I1uQlPnlwupjbGmHMTsvkYRGQqUBcojncy9aFAFICqThARwTsfbEPgGND5DHMI/ys2NlbdHkRPVVn/9XeMf+oZOh7ORcHt+/ku8jiLY6J5ZsJ4Sl1zlav5jDEmPRFZrqqxGa4L94l6skJh8KWpqST+uIIx/Qcy+ZfFzC1UnYI3xpLQ4kYuatOMiMisdJJmjMmpzlQY7M7nTCYeD3nq1GDID1+ybsdWYob0ImnDVpp0bs8dMReSMH0empzsdkxjjDktO2MIgZSkJN57cjhRX//MtbtOkFohhu+vu5AOY0cRGR3tdjxjTA5kZwwui4iOpvPYp2j7+9eUmvI0X5/YR9dXxzC9RmOOzl1MuBdnY0z2YoUhhMTjIV/DOvRcOZ/PRz7P9flLsqvjQN6q2Zifps50O54xxgDWlOQqTU7m0PtfENujPaU0kk969afowB5ElijidjRjTDZnTUlZlERGUrhTc35Z/wev9nqUhI/msab23bzW7SFSrIPaGOMSKwxZQNHyZak+bjDlv5vMJ3mO0+fNl/mufnuS/tzkdjRjTA5khSELib64IkNXLebbEWOptPMwW2/uypcPD+VkYqLb0YwxOYgVhizG4/FQb/CjlF8yhb3XXErTl0bS//LrSfp7i9vRjDE5hBWGLCqyZFFqfvwKb/UbSHstyrb6Xdnxxkekpqa6Hc0Yk81ZYcjCPB4PHcaMpNoPU4muVY3Wve/nrsrVSNl30O1oxphszApDGIgsVZzS08bQuFlTah5IYWu9Lhz/eaXbsYwx2ZQVhjARERnJkx+/S9/Fs/HkzsX0hu0Y3ryNXdZqjMl0VhjCTK4rLqbst5NYEBPFB1/MZkv7J0g5fMTtWMaYbMQKQxiKKJifd9Yu46vnXkUXxrH51m6s+/Y7t2MZY7IJKwxhyuPxUOnRzpSZ+SLjNiynVoP6/PXeDLdjGWOyASsMYS7P9VfR98uPeKxKTSL6vsSBl9630VqNMefFCkM2ULHmVQxasYj8zW/mt+HjuO/K60g8nOB2LGNMmLLCkE148uSi5BvDWFHvMmatjmNFiwdIOXDY7VjGmDBkhSEbEREenfkuKyZPp8Qf29jeuCfxK9e5HcsYE2asMGRDFTq0oMyMF5m6cSWX1qjOmtlfuR3JGBNGrDBkU3muvZJGk8dxW5FyRPUew7FFy9yOZIwJE1YYsrErGtXn3dW/kOfCcvx9z6PMHDTa7UjGmDBghSGbi4wpRpnPXmFi/qO0eWoQa8a84XYkY0wWZ4UhB4goVIAxvyzg/fqtyf3cFA6Mm+J2JGNMFmaFIYfIV7Qwrb6cQv67buWrYWN4/OY7bG4HY0yGrDDkIBIVScnxA/mucmE+WTyfLY+NsbukjTH/I9LtACa0JCKC8csWsfHxMaS++wXxngiKPf0wEZH2o2CM8bIzhhzI4/FQeczjFOrdhsdffYHONa63eR2MMf+ywpBDiQhFB99PsdrVifx7G3sfHYNan4MxBmtKytE8Hg9jFnzBgefe5uDYdzmYfIKLXh5kzUrG5HAhPWMQkYYi8qeIrBeRJzJYX0FEForIbyKyUkQahzJfTuTxeCjavyue+1vS6I3n6HPtLdYhbUwOF7LCICIRwHigEVAVaCMiVdNtNgiYrqpXAfcAr4UqX04mIlwwvA931alL3b8PsG/gy1YcjMnBQnnGUAtYr6obVTUJmAY0S7eNAgWd54WAHSHMl6N5PB6eWziHm/vcx6FJM/ihzxC3IxljXBLKwlAW2OrzepuzzNcwoJ2IbAPmAn0y2pGIdBeROBGJi4+PD0bWHElEKDaiN381qM5N40cxvk03tyMZY1yQ1a5KagNMVtVyQGNgioj8T0ZVfUNVY1U1tkSJEiEPmZ2JCA3efZFB1zbghm/WcHDix25HMsaEWCgLw3agvM/rcs4yX12B6QCq+hOQGygeknTmXxGRkQxbPIeSt9dj58BxzBnyrNuRjDEhFMrCsAyoIiKVRCQab+fy7HTb/APUBxCRS/EWBmsrcoFERhIzcShvlkyl2cgnWPnWVLcjGWNCJGSFQVWTgd7AV8A6vFcfrRGRESLS1NmsL9BNRH4HpgKd1C6PcY3kimbogs947cpbyD/0TY4v+dXtSMaYEJBw/70bGxurcXFxbsfI1lL2H2J7sz78sXE9hV94nFqtm7sdyRhznkRkuarGZrQuq3U+mywoomghSk8fyyMH1tG5Y0eSNm5zO5IxJoisMBi/RJUuwdRPZvBa6avZ2bofyXv2ux3JGBMkVhiM365oWJ/aH79Cyu69vFq3BYd37XE7kjEmCPwqDCJSREQ6iMgsEVkjIl+ISDcRKRnsgCZryR17GQeGdqXfusU826AletKG6zYmuzlrYRCRT4BPgJJAf1W9DOiF91LS90VkUVATmizn6q5tmDvoGTruSCW+3/M2rpIx2Yw/4yt3UdWDvgtU9R/gFeAVESkchFwmi2swsj/78xVj85i3+GTPRu6f+obbkYwxmeSsZwzpi0Kg6032VaR/V6ZUysWD095k9SuT3Y5jjMkkAXc+i8itIjJJRKo7r7tneioTFkSEZ76by+x6rcnz9Lt2A5wx2cS5XJXUBXgM7yioNwPVMzWRCSu58uXl1llvEHVhOWbefT8rv5zvdiRjzHk6l8KQoKoHVbUf0AComcmZTJiJKFSAIu+M5LH4VfRt08HucTAmzJ1LYZhz6omqPgG8l3lxTLgqUKUin3/4Ec/nr8KujgNITTzhdiRjzDny53LVvSLS5NRrVf3Md72qvhKMYCb8xLZqSuUJwzm6bDWvNWpDamqq25GMMefAnzOGk8DrInJf+hUiYmMxmzTy334TixtfQZ9Fs5jV43G34xhjzoE/9zHsBO4E5olIOVUd5rPukqCkMmGt6+RXKZqQRPVPfyah8bcUaHGL25GMMQHwq49BVTcDdYAGIvJmRtNtGnOKx+Phzo9eJ/e1VxLXazA/T53pdiRjTAD8+QUvAKq6F+/saiWB2SKS99Q6Y9KT6Chi3hrBg4f+oFPnLiRt3+12JGOMn/wpDL+deqKqx4HmeOdqXgQUDEoqky1ElijKu9M+ZHyxauzuMIDUY4luRzLG+MGfITG6pHudqqo98F62WilYwUz2UKNpI2q/8yxJq/5mVqvudqWSMWHgnPsKVHU4VhiMH/I1uJ7fW15Hy7lTeLvjA27HMcacxXl1IqvqFhHpn1lhTPbV/NXRvFyvBTd8uYqj85a4HccYcwbnMojedJ/Hx8D/3N9gTHoej4cH5kwl31WXsqH7EP745ju3IxljTsOf+xjSO6yq/xYDEXk9E/OYbMyTJxel3htN64uqsqdpU1Zv3Uyu4kXcjmWMSSfgITGAp9JtMjBzI5nsLLJ0CUa8+DxD8lZiX8+RaEqK25GMMekEPCSGqm46tUJEpqqqDaVpAlK/e0eajRvJ8UXLWP3YM27HMcak409h2AncCPQTkWHp1tmQGOacFGx/B2tvvZLYFwcztf9wt+MYY3zYkBjGNQ0mPUeXC6tzyfsLOLF6vdtxjDEOGxLDuCZXvry88tO3FC1WjB0dnuTw1h1uRzLGYENiGJdFlixKzDsj6bXuO5rXvI6UpCS3IxmT49mQGMZ1ea6+jMZt7+GGIxEcGD3J7TjG5HhnvY9BRDyq+j8D3KjqcBGZfLr1xgSiz5svE1+kAofGTyPqiksodKfN4WCMW/xpSvpGRD4SkTYiUhBARPKKyJ3AKOBXfz9MRBqKyJ8isl5EnjjNNneLyFoRWSMiH/q7bxP+io96kD+qlOSq1s1Y+eW3bscxJsfypympPjAcqAjMEZGlwHzgCuBFVa3uzweJSAQwHmgEVAXaiEjVdNtUAZ4ErlfVy4CH/T0QE/4kVzRVXx5Msajc7On/IimHj7gdyZgcyd/LVdeq6tOqegNQT1WvVdVhqur32QJQC1ivqhtVNQmYBjRLt003YLyqHnA+d08A+zfZQKXY6nz37XwqxR9lT++nUBum25iQC/h+BOfKpHNRFtjq83qbs8zXxcDFIvKDiCwVkYbn+FkmjOW9rjpFhvRk1IwPGNmirdtxjMlxstqNapFAFaAu0AaYJCKF028kIt1FJE5E4uLj40Ob0IRE4R6t+KdUAf74+juOLvzF7TjG5CihLAzbgfI+r8s5y3xtA2ar6klnTKa/8BaKNFT1DVWNVdXYEiVKBC2wcY/H42HGrz/xbM1b2XP/CE5u3eV2JGNyjFAWhmVAFRGpJCLRwD3A7HTbfIr3bAERKY63aWljCDOaLCRX4YLEvDOK7YlH6HJtPRIPJ7gdyZgcwe/CICK5ROReERkgIkNOPfx9v6omA72Br4B1wHRVXSMiI0SkqbPZV8A+EVkLLAQeU9V9/h+OyW6iK5dna4db+WzH3/z4wGC34xiTI4iq+rehyDzgELAc+HcQfVUdG5xo/omNjdW4uDg3I5gQ+HvAWDyTPqXES09Q8N4mZ3+DMeaMRGS5qsZmtC6QGdzKqapdJWRccdHIh9n551be7d2fa6NSiG3V9OxvMsack0D6GH4UkWpBS2LMGUhEBPnG9mVUwgZG39+HlIPW32BMsARSGOoAy50hLVaKyCoRWRmsYMakV6hieb6dNoOnoiuy54FRdvObMUESSFNSo6ClMMZP1e5qwqEDiWx9Yizf93iMuye52sVlTLbk9xmDqm4BCgN3OI/CzjJjQqpg1zuZUM5Duzdf5M/pX7gdx5hsJ5DLVR8CPsA7g1tJ4H0R6ROsYMacjojw1JczmXzZzUQPmUDyDhtSy5jMFEgfQ1egtqoOUdUhwDV4B70zJuQKlY6h5ay30MQkvr/nQZKOHnM7kjHZRiCFQfC5f8F5bnM+G9dEV7mA4wM60ej7jxnS6C634xiTbQRSGN4BfhaRYSIyDFgKvBWUVMb46bLu9zLytru4c81+jsya73YcY7KFQDqfXwC6APudR2dVHRekXMb4rd/sDyhTuwa7H36GXct+dzuOMWHP7yExsiobEsMAJO/YQ/uqsaw4mUDc5vXkK1HM7UjGZGlnGhLjrGcMIrLE+TdBRA77PBJE5HBmhzXmXESWKUnrAf1oLkVIGPQq4f4HjzFuOusNbqpax/m3QPDjGHPumj/+IDdFFGL/M2+yP/YyinVr6XYkY8JSIPcxPOvPMmPcVPiR9myOrUxszw789MEMt+MYE5YCuSrp1gyW2TAZJksRj4fLXxlKTJ78HBgxgZT9h9yOZEzY8aePoaeIrAIucQbPO/XYBKwKfkRjAhNzUSUWfb+YS46ksrvnSBtsz5gA+XPG8CHesZFm8984SXcAV6tq2yBmM+ac5a7+fxQd1YfnvviY4Xfc43YcY8LKWQuDqh5S1c2q2gY4DMQAFwCXi8iNwQ5ozLkq1LEZm8sUZNX87zi64Ge34xgTNgLpfL4PWIx3Xubhzr/DghPLmPPn8XiYtuwHXqx5G3t6jiR5+263IxkTFgLpfH4IqAlsUdV6wFXAwWCEMiaz5C5ckFKTR7HrWAL3X3eLDbZnjB8CKQyJqpoIICK5VPUP4JLgxDIm80RXrsDfbeoy9Z+1fN97sNtxjMnyApnBbZuIFAY+Bb4RkQOATdRjwkKH54YTezKaXO/PI2HWtxRocYvbkYzJsvwuDKrawnk6TEQWAoWAL4OSypgguPS5/uz4azsfdu/H9XnGcnnD+m5HMiZLCqTzeb6INAZQ1e9UdTYwPmjJjMlkEhVJnnGPMeDAOoa270bqEetvMCYjgfQxVAL6i8hQn2UZjsxnTFZVtEol5r09hZFSlj0PP2uD7RmTgUAKw0GgPhAjIp+LSKHgRDImuGp2uJvSg+5n/6ffMqfvMLfjGJPlBDS1p6omq2ovYCawBCgZnFjGBFfhPvfycsmT3PXiSNbP/sbtOMZkKYEUhgmnnqjqZKAT8HUm5zEmJESEkXNnMLHKDUQOGE/y7n1uRzImywhkas+J6V4vV9UumR/JmNAoUbECbWa9RerhIyy992GST5xwO5IxWYLN4GZytFyXXcTR/u1psOBDhjZu5XYcY7IEfwbR+3cGN1Ut6PMooKoFgx/RmOCq9kBHnqh3O7ev2M2R2QvdjmOM60I6g5uINBSRP0VkvYg8cYbt7hIRFRG7HNaExOB5Myhf6yp29xnNnriVbscxxlUhm8FNRCLw3hDXCKgKtBGRqhlsVwDvgH02TrIJGYmOotRbIxhxeD033XQTR3bHux3JGNeEcga3WsB6Vd2oqknANKBZBtuNBJ4FEgPYtzHnLbJMSe4a1JfGUpjD/cfZzW8mxwrlDG5lga0+r7c5y/4lIjWA8qo650w7EpHuIhInInHx8faXnck8Tfv2ZvBzT3N8zmL2vjTF7TjGuMLvGdyAtsANQEdV3QLkF5FamRVERDzAC0BfPzK9oaqxqhpbokSJzIpgDACFerZmx01XUKvf/XzzyiS34xgTcoH0MYwHrgXaOK8TCGwQve1AeZ/X5ZxlpxQALgcWichm4BpgtnVAm1ATES5/aRBl8xfi6PPvcnLrLrcjGRNSgRSG2qr6AE7bv6oeAKIDeP8yoIqIVBKRaOAevM1TOPs7pKrFVbWiqlYElgJNVTUugM8wJlMUKVuaBct/5sqI/OzqNJCUY9blZXKOQArDSefKIgUQkRJAqr9vVtVkoDfeuaLXAdNVdY2IjBCRpgHkMCYkoitXoOSEwby6dAH31a5rndEmxwikMLwMzMI7uupTeAfRGx3Ih6nqXFW9WFUrq+pTzrIhztwO6beta2cLxm35GlxPYq1L2fPXBvZNnO52HGNCIpAZ3D4QkeV4h94GaK6q64ITy5isY+zXn7Gn40AODX+dfFdcQp7rqrsdyZigCuTO51xADbxTehYDWonIkGAFMyariIyKImbCEA6VKcrdjZqw9Ve7M9pkb4E0JX2G94a0ZOCoz8OYbC+iYH50aDd+PBLPoq6Pk3rcRmI12ZffTUlAOVVtGLQkxmRxVzdtxNoZX3Ck+0ji+z1PyVcHICJuxzIm0wVSGH4UkWqqGsgwGMZkKzEtGhC1fjuzRj1PYsI2HnjvdbcjGZPpAmlKqgMsd0ZHXSkiq0TEGltNjlP40Q58VDCJSVM/IGGBjfVosp9Azhj8HknVmOzMExHBR0sXs+fuvuztMZzcX71B1IXl3I5lTKYJZGrPLRk9ghnOmKyqSNnSXPjB8yShDKx3hw3TbbIVv88YROTRDBYfApar6opMS2RMmIiqWIZ/HrqLsX26ceFd99Hju0+QiAi3Yxlz3gLpY4gF7sc7VHZZoAfQEJgkIo8HIZsxWd6tve9j6ciXuPXP/ewfbSOxmuwhoMtVgRqqegRARIYCc4AbgeXAc5kfz5isL3ZgH/YeTCLuhUkkHN9Ly9GD3I5kzHkJ5IyhJOB7V89JIEZVj6dbbkyOIiIUf/phno7eS9/nniJh6Qq3IxlzXgIpDB8AP4vIUOds4QfgQxHJB6wNSjpjwoRERfLBwq/5sGp99nUZysltu92OZMw5C+SqpJFAd+Cg87hfVUeo6tEAp/g0JlsqW/USas54ldTjiUxs0Iqj8fvcjmTMOQnkjAFgA/AT8BuQV0RuzPxIxoSv6IsrsmdAex5c9z1jGrZCU1LcjmRMwAIZXfU+YDHeiXaGO/8OC04sY8LXNd3aMbvfMO7dksS+4TZkhgk/gZwxPATUBLaoaj3gKrxNSsaYdG4fM5Si3Vryz/gPmNnXRqc34SWQwpCoqongnZtBVf8ALglOLGPCX7FRfXihyDE6vjiarZ987XYcY/wWSGHYJiKFgU+Bb0TkM8CGxDDmNCQigpcWzuX96g052XcsJ9asdzuSMX4J5KqkFqp6UFWHAYOBN/FO3GOMOY1i5cty+xdv4ymQj0/u6MyOlWvcjmTMWQXS+RwrIrNE5FfgFWA03juejTFnEFmmJLlfH0j3zUt55LbmpCbYxIcmawv0Brd3gLuAO4DbnX+NMWdR6vqrmfnCqzyZEsOuLoPRk8luRzLmtAIpDPGqOltVN9mw28YErsHDPbjwxSdJWPgLk25vR2pqqtuRjMlQIIVhqIi8KSJtROTOU4+gJTMmGyrYtglf33wxPb7+iDk9+rsdx5gMBTK6amfg/4Ao4NSfOgp8ktmhjMnOek2dRMlkoeqnSzl03SwKdW7hdiRj0gikMNRUVbtvwZjzFBERwd3TJ7Kr0yDi+o0mccdmmgx8xO1YxvwrkKakH0WkatCSGJODSGQkMW8MY7js4v6hAzi08Ge3Ixnzr0AKwzXAChH5U0RWisgqEVkZrGDGZHeevLmZumQBU6rfxr4uQzjx+59uRzIGCKwwNASqAA2wy1WNyRRl/68KdT6fREThArzUsDUbvrczB+O+QO583pLRI5jhjMkJIkuXwDPhSUbuWcPou9qRvN0m+THuCnQ+BmNMEJSvfTXfz/iMJ/JVYkfLR0mOP+B2JJODhbQwiEhDp49ivYg8kcH6R0VkrdOHMV9ELghlPmPcVL1FY8pNfZ4DW7fTs8YNHN6xy+1IJocKWWEQkQhgPNAIqAq0yeAqp9+AWFW9ApgBPBeqfMZkBXmuuYLNfVrw3ra1zG1+H6lHjrkdyeRAoTxjqAWsV9WNqpoETCPd6KyqulBVT/1PWAqUC2E+Y7KEZv0fYtW706n5zxF2tnmc1KPH3Y5kcphQFoaywFaf19ucZafTFfgyoxUi0l1E4kQkLj4+PhMjGpM1XNyhJTGvD+a7Jd9zd9UanDh8xO1IJgfJkp3PItIOiAXGZLReVd9Q1VhVjS1RokRowxkTIvlb1Gdnizqs2L6Fv9s/TmriCbcjmRwilIVhO1De53U5Z1kaInILMBBoqqr2P8HkaH3fHs8Pk94nz4+r2NlhAElHbC4HE3yhLAzLgCoiUklEooF7gNm+G4jIVcBEvEVhTwizGZNlxXS+k+IvPM6Az6fSqurVnDxqHdI5napycPxUUvYfCsr+Q1YYVDUZ6A18BawDpqvqGhEZISJNnc3GAPmBj0VkhYjMPs3ujMlRCrW/gyrNG1I8PoE9nQaRetxOpnMqTU1lb/8X2DfsNRKmzwvKZ4iqBmXHoRIbG6txcXFuxzAmJA5/MIf4R57l8NUX83/vjyFPsSJuRzIhpCkpxD86hoQP51C4T1uKDu6BiJzTvkRkuarGZrQukGG3jTEuK9i2CUkepUH7u7ms2iI+/fN3PAXyuR3LhIAmJ7Onz2iOzPiGIv06UeTxLudcFM4mS16VZIw5veJtbqffA31od6IAO1o8RMpeGz4juzt28BA9r7ye7dO/pOiAbhTt3zVoRQGsMBgTlnq+9CxNpo4n6c9NTK7Tgp0r17gdyQRJ6pFjzG/ejclr41jd6jqKPNIh6J9phcGYMJWvwXVETRzEQ38toe/Nd5C04R+3I5lMlhS/nx0tH6HqH3v4/bV36TBhbEg+1wqDMWGsfON6zHtvKgOKXcL2Jr1I/HWt25FMJtmwNI4rLriQhXG/UOrtEVzSs13IPtsKgzFhrk67Vlw8dyIp+XJz9w038824iW5HMufpxNoNHO0ylMLqoczIB8nX+MaQfr4VBmOygejKFcg9eQRrSeSXIS9w+IM5bkcy52jJpCn8c3svCkbmYtHyZdTv2TnkGawwGJNNlKtWlZVbNtL5tibEP/wMqx57htSUFLdjmQAsf34i9bp35M3UPZSd+zq5q1Z2JYcVBmOykbzFi1L6w+dIaH4DN44dzJO1brG7pMOApqayb/QkCj/7Ps9efQtPLJlDVLkY1/JYYTAmm5GoSKpNGEmvxs25bVMCO1o8SPKuvW7HMqdxeHc8915yFb8/P5EC7W7nkZ/mUqRcGVczWWEwJhuKiIjgqS8+ptaUF0hat4mB1W/i1xk29FhWc3LrLlY17ck3G9bx1x21KfHC40iU+wNSuJ/AGBM0+ZvcyKEPRvFOvZs42LkPz6dEU6B1Q7djGeD3qbMoOOJtyiQls/rTeZRqerPbkf5lZwzGZHNl69Tit99XMKDeHezp/RSrew21obtdpKrMeWQINe69k3kpByn71cQsVRTACoMxOULZqpdQ4ZNx5O5xF3e+8Tx3XVSN5O273Y6V46QmHGV31yFcNGUBj9eoS9vvZhB9UQW3Y/0PKwzG5BASGUnZUQ/zSO8+tNIibK3XhaNf/eB2rBzjp6kzqV+uCrs/X0DM8AcYHbeAgqXdu/LoTKwwGJPD9Bz3DK1/mEVkuRjGtexEvzoNST523O1Y2ZaqcnDCdDb1HsWGxEOkvNyPwg+0CeroqOfLCoMxOVB05fKU+3ICGy6O4edffmZ7o56cWLfR7VjZzpbfVvFWnebsG/wK9Ro35M8tm6nRpoXbsc7KCoMxOZTkiubt5Uv49JOZaPx+Vt/ckdfa9yQ1OdntaNnCkc8X0femRvRbOpfooT0o9d5o8pQq4XYsv1hhMCaHK3b7zZRf/C4zS0fw0PsT+eHWTjaE93nYsfZPVrd/jN1dBjO4Rj2WzPuG8r3bZemmo/SsMBhjiChehJHL5vP1qBcpt3Uf227qzLyHh5ByIsntaGFDVdk7bQ6xV1bnsRnvUuSJrlzxzTtcfmtdt6MFzAqDMQbw3i1db+BDlF8yhW01q9DkpZEMqXo9ib+scjtalrf151/Z2bofh/o8w6Cq1/PMx+9TtG+nLHEX87mwwmCMSSOyVHFu+OR1JvYdyN25S7G9SS9+6/Q4e9dvdjtalpN6LJFPuj5K5Wti+f777yn21EP0/PUrrry9gdvRzouoqtsZzktsbKzGxcW5HcOYbCn1yDEOvPgezUYNYKcm8ctzr1Pk/rvx5M7ldjRXJZ88yaa3phP12kyO7NjFK2Vh8JRJlLviMrej+U1ElqtqbIbrrDAYY87ml8/msGrsm9Rdt5eIsiVZ26w2DQc/iicyPJtKzpWqcnzRMpq1asmuQweYU/ceSozoTZ5rr3Q7WsCsMBhjMsXxJb8yq88A2q74iperXE/nZ4eTr1k9xJP9W6Xnv/YWVebEcfKX1cwvkEJks7p0en4knogIt6OdEysMxphMczIpifeeGM5NSzchf2/lp5hcpNxSi3bPDCUimzUxaWoqx775iTlDx9D65894pXws7YcNoGDbJkiuaLfjnRcrDMaYTKcpKRyZvZDW3bvy14F4vr70Nor0aEW+1g2JLlHU7XjnJeXIMV576HE8S36n4X7wlC3JnCtj6PLS0+QtXMjteJnCCoMxJmiSk5PZMH0OBWYsIGHxchof+o2u19bj0eefJlfNy8Pqxq6tS34hz5wfSfhoHq22/kTRQoWYOeFN8je7OWwvPT0dKwzGmJDY9dOv9O3Rk5t3n6Ruaj6Oli/BxzFC92EDKV+7htvxMpS8M54js+Yz7NlneHvrapaWup5iTeuTfGc9St9yPZ5s2n9ypsKQvUqgMcZVpa6twQcrfyb1yDGOfLaQhS+/xqh5c6m5dCN6eTX2176UxKuqUPPuZkRERbmSUVVJWreR7ydNof87r/OipwJlI3JT/8KKlLq2BuWeHUHBC8q5ki2rCGlhEJGGwEtABPCmqj6Tbn0u4D3gamAf0FpVN4cyozHm/Hny56Vg2ybc17YJNy/7jSLL/+bYVz/w+uuv8/axbfw2qD7Frr+ajZWK4bn4Amq3bIYnT3A6rvVkMklrN7Bl0Y88+cqLtDyZn1rHIpCU45zUVE52up3yvTpT+aIK3BGUBOEnZE1JIhIB/AXcCmwDlgFtVHWtzza9gCtU9X4RuQdooaqtz7Rfa0oyJnzs+HM9S96bRt3j0Rz//ld6rPyWtclHWFjiGqL/ryLTUvYRUbIo3dq1J7JMSRLyRFL4wguIzp/vjPvVk8mk7NnHsoWLyXX4GJWSPBz9YyMtpr1GYwrRISqGRE2hweEVPHljYzp07ULeujWJLB0eo50GQ1ZpSqoFrFfVjU6oaUAzYK3PNs2AYc7zGcCrIiIa7h0hxhgAylxyEXc/Nejf1+OX/cqWxUspfCKKE7//yRfz5qHHk2i22Ds3xF0HV1BIonin5FVI/jz02PUb5XPlZ0iF6qBKi7XfUtWTl+F5KgHQdN/P3BxdlNEFqhBRugQx+QtR/OpaxHTuQq4al7K1XEy27TPITKEsDGWBrT6vtwG1T7eNqiaLyCGgGLDXdyMR6Q50B6hQIevNl2qM8U/lmjWoXPO/TunveJ4Th4/A9j2k7Iynz8wZeI4cp1Dly0g9coyKX+2nbL6C5L7yMhChXtQBLigZQ5EbbyGyRFHe2/0P5f/vYirVq4Mnf16+cPHYwlkom5JaAg1V9T7ndXugtqr29tlmtbPNNuf1BmebvRntE6wpyRhjzsWZmpJCeU61HSjv87qcsyzDbUQkEiiEtxPaGGNMiISyMCwDqohIJRGJBu4BZqfbZjbQ0XneElhg/QvGGBNaIetjcPoMegNf4b1c9W1VXSMiI4A4VZ0NvAVMEZH1wH68xcMYY0wIhfQ+BlWdC8xNt2yIz/NEoFUoMxljjEnLrtsyxhiThhUGY4wxaVhhMMYYk4YVBmOMMWmE/bDbIhIPbDnHtxcn3V3VOYAdc85gx5wznM8xX6CqGQ4WFfaF4XyISNzp7vzLruyYcwY75pwhWMdsTUnGGGPSsMJgjDEmjZxeGN5wO4AL7JhzBjvmnCEox5yj+xiMMcb8r5x+xmCMMSYdKwzGGGPSyBGFQUQaisifIrJeRJ7IYH0uEfnIWf+ziFR0IWam8uOYHxWRtSKyUkTmi8gFbuTMTGc7Zp/t7hIRFZGwv7TRn2MWkbud7/UaEfkw1Bkzmx8/2xVEZKGI/Ob8fDd2I2dmEZG3RWSPM5FZRutFRF52vh4rRaRGRtsFRFWz9QPvEN8bgAuBaOB3oGq6bXoBE5zn9wAfuZ07BMdcD8jrPO+ZE47Z2a4AsBhYCsS6nTsE3+cqwG9AEed1Sbdzh+CY3wB6Os+rApvdzn2ex3wjUANYfZr1jYEvAQGuAX4+38/MCWcMtYD1qrpRVZOAaUCzdNs0A951ns8A6ouIhDBjZjvrMavqQlU95rxcindGvXDmz/cZYCTwLJAYynBB4s8xdwPGq+oBAFXdE+KMmc2fY1agoPO8ELAjhPkynaouxjs/zek0A95Tr6VAYREpfT6fmRMKQ1lgq8/rbc6yDLdR1WTgEFAsJOmCw59j9tUV718c4eysx+ycYpdX1TmhDBZE/nyfLwYuFpEfRGSpiDQMWbrg8OeYhwHtRGQb3vlf+oQmmmsC/f9+ViGdqMdkPSLSDogFbnI7SzCJiAd4AejkcpRQi8TbnFQX71nhYhGppqoH3QwVZG2Ayao6VkSuxTsr5OWqmup2sHCRE84YtgPlfV6Xc5ZluI2IROI9/dwXknTB4c8xIyK3AAOBpqp6IkTZguVsx1wAuBxYJCKb8bbFzg7zDmh/vs/bgNmqelJVNwF/4S0U4cqfY+4KTAdQ1Z+A3HgHm8uu/Pr/HoicUBiWAVVEpJKIROPtXJ6dbpvZQEfneUtggTq9OmHqrMcsIlcBE/EWhXBvd4azHLOqHlLV4qpaUVUr4u1Xaaqqce7EzRT+/Gx/ivdsAREpjrdpaWMIM2Y2f475H6A+gIhcircwxIc0ZWjNBjo4VyddAxxS1Z3ns8Ns35Skqski0hv4Cu8VDW+r6hoRGQHEqeps4C28p5vr8Xby3ONe4vPn5zGPAfIDHzv97P+oalPXQp8nP485W/HzmL8CGojIWiAFeExVw/Zs2M9j7gtMEpFH8HZEdwrnP/REZCre4l7c6TcZCkQBqOoEvP0ojYH1wDGg83l/Zhh/vYwxxgRBTmhKMsYYEwArDMYYY9KwwmCMMSYNKwzGGGPSsMJgjDEmDSsMxhhj0rDCYIwxJg0rDMacgYgUFpFePq9/DNLnlBOR1sHYtzGBssJgzJkVxjtfBwCqel2QPqc+3jH3jXGdFQZjzuwZoLKIrBCRMSJyBEBEKorIHyIyWUT+EpEPROQWZ3jrv0Wk1qkdiEg7EfnF2cdEEYnw/QARqYN35NeWzjYXhvQIjUnHCoMxZ/YEsEFVq6vqY+nWXQSMBf7PedwL1AH6AQPg30HcWgPXq2p1vOMVtfXdiaouwTs4XDPnc8J5kDuTDWT7QfSMCaJNqroKQETWAPNVVUVkFVDR2aY+cDWwzBmsMA+Q0Wi2lwB/BD2xMX6wwmDMufOdwyLV53Uq//3fEuBdVX3ydDtxhsM+5MweaIzrrCnJmDNLwDvJz7maj7fvoCSAiBQVkQvSbVORMJ+X2GQvVhiMOQNn7oIfRGS1iIw5h/evBQYBX4vISuAbIP1E7X/gHWt/tYgE66onY/xm8zEYY4xJw84YjDHGpGGFwRhjTBpWGIwxxqRhhcEYY0waVhiMMcakYYXBGGNMGlYYjDHGpPH/MbO9VBnTbc0AAAAASUVORK5CYII=\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": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "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": 16,
   "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": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_gd = pvqd.evolve(problem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Energy over time')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEYCAYAAAB2qXBEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABBrUlEQVR4nO3deXxU1dnA8d8zk30PENawKrtAhLAVwQVRVEBrtWK1gm+rde372lcral8X1LZWa2srtVXrrkWlBWlFUdyXUgmIyC4gQpAlCWRfZ+Z5/7iXOMQACUlmMsnz/Xzmk7nnLvPcZDLPnHPuPUdUFWOMMeZYecIdgDHGmMhmicQYY0yTWCIxxhjTJJZIjDHGNIklEmOMMU1iicQYY0yTWCIxxiAifxaR/wt3HCYyid1HYiKdiGwHugD+oOKnVPW68ETUuonIbODHqnpSuGMxbUNUuAMwpplMV9VlLfkCIhKlqr6WfI3mFokxm8hjTVumTROR2SLyoYg8ICIHRORLETkraH2qiPxVRHaLyC4RuUdEvEH7fiQivxORAuBOEekoIv8UkWIRWeFu/6G7/TwR+W2d118sIjccJrbvuMcocn9+xy2/SERy6mx7g4gsdp/HuuezQ0T2us1S8e66U0QkV0RuFpE9wJN1jjMY+DMwXkRKRaTQLX9KRO6pc4yfi8g+93dznoicLSKbRWS/iNwadEyPiMwRka0iUiAiL4lIh2P4c5kIZYnEtAdjgU1AJ+A3wF9FRNx1TwE+4HjgROAM4Md19t2G03R2LzAPKAO6ArPcx0FPAxeLiAdARDoBpwMv1A3I/aB9FfgD0BF4EHhVRDoC/wQGikj/oF1+EHScXwMDgCw37h7A7UHbdgU6AL2BK4NfV1U3AFcB/1bVJFVN+9Zv65tjxAUd+zHgUmAUMBH4PxHp6257PXAecDLQHTjg/p5Me6Gq9rBHRD+A7UApUBj0uMJdNxvYErRtAqA4H5RdgCogPmj9xcA7QfvuCFrnBWqAgUFl9wAfBi1vAKa4z68Dlhwm5h8Cn9Qp+zcw233+HHC7+7w/UOLGLjiJ7Lig/cYDX7rPTwGqgbgj/L5mB8fslj0F3BN0jArA6y4nu7+zsUHbrwTOCzrnyUHrurm/p6hwvzfsEZqH9ZGYtuI8PXwfyZ6DT1S13K2MJOF8a48Gdn9TQcED7AzaN/h5Bk6/4uHWg1MruRR40/350GFi6g58VafsK5waADi1j98Cc3FqI4vc2DvjJJSVQTELTpI7KE9VKw/zug1VoKoHL16ocH/uDVpfgfM7BKfms1BEAkHr/TiJelcT4zARwBKJac924tRIOunhO6SDL2vMw2kGywQ2u2U962z/HLBWREYAg4FFhznu1zgfwMF6Aa+7z98EMkQkC6eWdLCfJR/nQ3yoqh7uQ/pol2I296WaO4H/UtWPmvm4JkJYH4lpt1R1N/AG8FsRSXE7jY8TkZMPs70f+AdOp3uCiAwCLquzTS6wAngW+LuqVnz7SAAsAQaIyA9EJEpELgKGAP9yj1MDvAzcj1NzetMtD+D0V/zOrZ0gIj1E5MxGnPpeIFNEYhqxz5H8GbhXRHq78WSIyLnNdGwTASyRmLbin+5VSAcfCxu432VADLAep5N4AU4b/+FcB6TiNJc9C/wNp1YT7GlgmLu+XqpaAEwD/hcoAH4OTFPV/KDNXsDprH+5To3pZmALsFxEioFlwMAjnuWh3gbWAXtEJP9oGzfAQ8Bi4A0RKQGW41ykYNoJuyHRmCYQkfuArqo6K6hsEk4TV2+1fzDTDliNxJhGEJFBIjJcHGOAHwELg9ZHA/8NPG5JxLQXlkiMaZxknH6SMuBFnCurXoHam/0KcZrGfh+e8IwJPWvaMsYY0yRWIzHGGNMk7fI+kk6dOmmfPn3CHYYxxkSUlStX5qtqRt3ydplI+vTpQ05OztE3NMYYU0tE6o7GAFjTljHGmCayRGKMMaZJLJEYY4xpknbZR1KfmpoacnNzqaxs6qCpprWKi4sjMzOT6OjocIdiTJtiicSVm5tLcnIyffr0IWh4btNGqCoFBQXk5ubSt2/fo+9gjGmwkDZtichUEdkkIltEZE496yeJyCoR8YnIBfWsT3GnAH04qGyUiHzuHvMPcoxZoLKyko4dO1oSaaNEhI4dO1qN05gWELJE4s6DPQ84C2e47ItFZEidzXbgzN72ralJXXcD79cpewS4AmcWuf7A1CbEeKy7mghgf19jWkYoayRjcKY83aaq1cB84JA5C1R1u6quAQJ1dxaRUTgzrr0RVNYNSFHV5e4Aec/gzB1tjDEmSEl5AH+gZYbECmUi6cGh05Lm8s20okckIh6cwfFurOeYuQ05pohcKSI5IpKTl5fX4KBDyev1kpWVxYgRIxg5ciQff/wxANu3byc+Pp6srKzaR3V1NQCvv/46Y8aMYdCgQWRlZXHRRRexY8cOAJYvX87YsWPJyspi8ODB3HnnneE6NWNMGJWUB/jZ7/fy2+f3t8jxI6Wz/RpgiarmHmvzhKo+CjwKkJ2d3SpHqoyPj2f16tUALF26lFtuuYX33nsPgOOOO6523UFr167l+uuvZ/HixQwePBiAxYsXs337dnr16sWsWbN46aWXGDFiBH6/n02bNoXydIwxrUB5ZYA5D+9jx54arjwvrUVeI5SJZBeHzm+d6ZY1xHhgoohcAyQBMSJSijMzW+YxHrNVKy4uJj09/Yjb3Hfffdx66621SQRgxowZtc/37dtHt27OZH9er5chQ+p2SRlj2ro/vnSATTuqufPHnRg9JL5FXiOUiWQF0F9E+uJ82M8EftCQHVX1koPPRWQ2kK2qc9zlYhEZB/wHZ9rUPzY10IdfPsDW3OqmHuYQx2XGcN2FR04MFRUVZGVlUVlZye7du3n77bdr123dupWsrCwAJkyYwLx581i3bh033li3te8bN9xwAwMHDuSUU05h6tSpzJo1i7i4uGY5H2NMZPivGamcNCKeCSMSWuw1QtZH4s45fR2wFNgAvKSq60RkrojMABCR0SKSC1wI/EVE1jXg0NcAj+PMYb0VeK1FTiAEDjZtbdy4kddff53LLruMg/PFHGzaWr16NfPmzfvWvgUFBWRlZTFgwAAeeOABAG6//XZycnI444wzeOGFF5g69ZgvaDPGRBB/QPnXh6X4A0pGWlSLJhEIcR+Jqi4BltQpuz3o+QoObaqq7xhPAU8FLecAJzRnnEerOYTC+PHjyc/P50gXBgwdOpRVq1YxYsQIOnbsyOrVq3nggQcoLS2t3ea4447j6quv5oorriAjI4OCggI6duwYilMwxoSBqvK7F/az5OMy0lM8TBjeskkEbKytVmvjxo34/f4jfuj//Oc/595772XDhg21ZeXl5bXPX3311doazRdffIHX6yUtLa3FYjbGhJeq8ud/FLLk4zIuPSslJEkEIueqrXbhYB8JOG+Ip59+Gq/Xe9jthw0bxkMPPcRll11GcXExnTp1olevXtx1110APPvss9xwww0kJCQQFRXF888/f8TjGWMi29+WFvPyWyWce3ISl09LDdnrWiJpRfx+f73lffr0Ye3atfWuO+ecczjnnHPqXTd//vxmi80Y07rlFfp47vViJo9O4PoL00M6koMlEmOMaQMy0qJ4+KYu9OoajccT2uGArI/EGGMi2Ir1FSx+vwSAfj1iiPKGfkw5SyTGGBOhNnxZxR2P5vPPD0up8YVvwA5LJMYYE4G2767hlj/l0SHVy33XdiY6KnyjW1siMcaYCLN3v4+b/7iPKC/85vrOdEgN79WY1tlujDERZuWGSsqrAvz+hi507xT+j3GrkbQie/fu5Qc/+AH9+vVj1KhRjB8/noULFzbpmHfeeechQ6YsW7bsmI6zevVqlixZcvQNgVNOOYWcnJxjep3msGjRItavXx+21zempZ09IYln7+zOcZkx4Q4FsETSaqgq5513HpMmTWLbtm2sXLmS+fPnk5ub+61tfT7fMb3G3LlzOf30049p38YkknCzRGLaIp9fufuJfD77wpkuOi259dxcbImklXj77beJiYnhqquuqi3r3bs3119/PQBPPfUUM2bM4LTTTmPy5MmUlpYyefJkRo4cybBhw3jllVdq97v33nsZMGAAJ5100iFzkMyePZsFCxYAsHLlSk4++WRGjRrFmWeeye7duwGnNnHzzTczZswYBgwYwAcffEB1dTW33347L774IllZWbz44ouHxF5RUcHMmTMZPHgw3/3ud6moqKhd98YbbzB+/HhGjhzJhRdeWDsO2Jw5cxgyZAjDhw+vHcF47969fPe732XEiBGMGDGidmKv5557jjFjxpCVlcVPfvKT2hs3k5KSuO222xgxYgTjxo1j7969fPzxxyxevJibbrqJrKwstm7d2jx/IGPCKBBQfvNsAe/klLNr37F9kWxJ4W9ca6Vu+N3eb5WdMjKBc09OprI6wC3zvj2Y4pnjEpk6PomiUj93PpZ/yLrf3dDliK+3bt06Ro4cecRtVq1axZo1a+jQoQM+n4+FCxeSkpJCfn4+48aNY8aMGaxatYr58+ezevVqfD4fI0eOZNSoUYccp6amhuuvv55XXnmFjIwMXnzxRW677TaeeOIJwKnxfPLJJyxZsoS77rqLZcuWMXfuXHJycnj44Ye/FdcjjzxCQkICGzZsYM2aNbXnkZ+fzz333MOyZctITEzkvvvu48EHH+Taa69l4cKFbNy4ERGhsLAQgJ/+9KecfPLJLFy4EL/fT2lpKRs2bODFF1/ko48+Ijo6mmuuuYbnn3+eyy67jLKyMsaNG8e9997Lz3/+cx577DF+8YtfMGPGDKZNm8YFF1xwxN+nMZHi0UWFLPuknMunp3L2hKRwh/MtlkhaqWuvvZYPP/yQmJgYVqxYAcCUKVPo0KED4DSF3Xrrrbz//vt4PB527drF3r17+eCDD/jud79LQoIzWFvwRFcHbdq0ibVr1zJlyhTAGZrl4ARYAOeffz4Ao0aNYvv27UeN9f333+enP/0pAMOHD2f48OGAM9Xv+vXrmTBhAgDV1dWMHz+e1NRU4uLi+NGPfsS0adOYNm0a4NTKnnnmGcCZiCs1NZVnn32WlStXMnr0aMCp/XTu3BmAmJiY2n1HjRrFm2++2aDfrTGR5KVlxby0zBk/69KpKeEOp16WSA7jSDWIuBjPEdenJnmPWgOpa+jQofz973+vXZ43bx75+flkZ2fXliUmJtY+f/7558nLy2PlypVER0fTp08fKisrG/RaqsrQoUP597//Xe/62NhYwPkwP9b+mIOvM2XKFP72t799a90nn3zCW2+9xYIFC3j44YcPmcSr7jFmzZrFr371q2+ti46Orh1PqKmxGtMaqSqbd1Qz6cR4rgvx+FmNEdI+EhGZKiKbRGSLiMypZ/0kEVklIj4RuSCovLdbvlpE1onIVUHr3nWPudp9dA7V+TSn0047jcrKSh555JHasuAh4esqKiqic+fOREdH88477/DVV18BMGnSJBYtWkRFRQUlJSX885///Na+AwcOJC8vrzaR1NTUsG7dkecQS05OpqSkpN51kyZN4oUXXgCceeTXrFkDwLhx4/joo4/YsmULAGVlZWzevJnS0lKKioo4++yz+d3vfsdnn30GwOTJk2vP3+/3U1RUxOTJk1mwYAH79u0DYP/+/bXneiyxGhMpVBUR4dbZHbnt8k54Qzx+VmOELJGIiBeYB5wFDAEuFpG6k4jvAGYDL9Qp3w2MV9UsYCwwR0S6B62/RFWz3Me+loi/pYkIixYt4r333qNv376MGTOGWbNmcd9999W7/SWXXEJOTg7Dhg3jmWeeYdCgQQCMHDmSiy66iBEjRnDWWWfVNgkFi4mJYcGCBdx8882MGDGCrKys2o7twzn11FNZv359vZ3tV199NaWlpQwePJjbb7+9tk8mIyODp556iosvvpjhw4czfvx4Nm7cSElJCdOmTWP48OGcdNJJPPjggwA89NBDvPPOOwwbNoxRo0axfv16hgwZwj333MMZZ5zB8OHDmTJlSu2FAYczc+ZM7r//fk488UTrbDcRadNXVVxz31727vfh8UhY71pvCDk48VGLv5DIeOBOVT3TXb4FQFW/1WYhIk8B/1LVBfWs6wh8CoxT1a9F5F3gRnemxAbJzs7Wuvc5bNiwgcGDBzf8hExEsr+zae1y99Xw0wf2EhsjPHxTVzqG+a71YCKyUlWz65aHsmmrB7AzaDnXLWsQEekpImvcY9ynql8HrX7Sbdb6PzlMI6KIXCkiOSKSc6Tpa40xJlz2F/m5+Y/7CCjcd33nVpVEjiRi7iNR1Z2qOhw4HpglIgd7sy9R1WHARPfxw8Ps/6iqZqtqdkZGRmiCNsaYBiqvDDDnT/s4UBLgV9dk0KtLdLhDarBQJpJdQM+g5Uy3rFHcmshanKSBqu5yf5bg9K2MOdYAQ9XMZ8LD/r6mNauqUaK8wh1XdGJw39hwh9MooUwkK4D+ItJXRGKAmcDihuwoIpkiEu8+TwdOAjaJSJSIdHLLo4FpOEmm0eLi4igoKLAPmzZKVSkoKCAuLi7coRhziEBA8fuV9GQvD9/YhbFD48MdUqOF7D4SVfWJyHXAUsALPKGq60RkLpCjqotFZDSwEEgHpovIXao6FBgM/FZEFBDgAVX9XEQSgaVuEvECy4DHjiW+zMxMcnNzsf6TtisuLo7MzMxwh2HMIR5dVMiOPTXcdWVGq78663BCekOiqi4BltQpuz3o+QqcJq+6+70JDK+nvAwYVbf8WERHR9O3b9/mOJQxxjTIy285d61/95QkoiKjX71eEdPZbowxbclbK8p45O+FTDoxnmsuaL13rTeEJRJjjAmxlRsrue+ZAkb0j+XW2a37rvWGsLG2jDEmxJLihRH947jjx52IiY7sJAKWSIwxJmTKKwMkxHkY2DuW+38akcMC1suatowxJgQOlPj5ya/28PzrReEOpdlZIjHGmBZWURng1j/lkV/oJ2tA27uXyRKJMca0IJ9fuevxfL7YWc3tP+7E0H6Rddd6Q1giMcaYFvTAc/v5ZH0lP7u4A+OHRd5d6w1hne3GGNOCsgfH0atLVKuca725WCIxxpgWkF/oo1NaFKePSTz6xhHOmraMMaaZLfukjEvv2M3nWyrDHUpIWCIxxphmtGJ9Bfc9U8CQPjEM7N32OtbrY4nEGGOaycbtVdzxWD59ukcz96qMNnHXekNYIjHGmGaQV+jjlj/lkZ7k4dfXdiYpvv18vFpnuzHGNIOOKV7OOzmZyaMTImau9eZiicQYY5qgpDxAaXmAbp2imHVOarjDCYuQ1r1EZKqIbBKRLSIyp571k0RklYj4ROSCoPLebvlqEVknIlcFrRslIp+7x/yDRPKg/saYiFJZHeC2R/L42e/3Ul3TfqfpDlkiEREvMA84CxgCXCwiQ+pstgOYDbxQp3w3MF5Vs4CxwBwR6e6uewS4AujvPqa2RPzGGBPM51fmPp7Pum1VXHV+ervpWK9PKGskY4AtqrpNVauB+cC5wRuo6nZVXQME6pRXq2qVuxiLG7eIdANSVHW5qirwDHBey56GMaa9CwSU+58tYPnaSv5nZjonj0wId0hhFcpE0gPYGbSc65Y1iIj0FJE17jHuU9Wv3f1zG3JMEblSRHJEJCcvL6/RwRtjzEEL3y3hzU/KuXx6KtMnJoc7nLCLmM52Vd0JDHebtBaJyIJG7v8o8ChAdnZ2+23MNMY02VnfSSI6Spg+se2On9UYoayR7AJ6Bi1numWN4tZE1gIT3f0zm3pMY4xpiH9/XkGFO8vhjEnJ2LU9jlAmkhVAfxHpKyIxwExgcUN2FJFMEYl3n6cDJwGbVHU3UCwi49yrtS4DXmmZ8I0x7dm7K8v4xZ/zePa1tjfDYVOFLJGoqg+4DlgKbABeUtV1IjJXRGYAiMhoEckFLgT+IiLr3N0HA/8Rkc+A94AHVPVzd901wOPAFmAr8FqozskY0z6sWF/BL58q4IR+sVzWTu8VORJxLnZqX7KzszUnJyfcYRhjIsDarVX8/I/7yOwcxYP/04WkhPYz9EldIrJSVbPrlrff34gxxhyFP6D89vkCOqV5+fV1ndt1EjmSiLlqyxhjQs3rEe65OgOvR+iQ0r7Gz2oMS6/GGFPH3v0+nn2tCFWlR0Y0XTvad+4jsURijDFB9hf7uekP+3hpWTF79/vDHU5EsERijDGu4jI/P//jPvIL/fzqms5WE2kgSyTGGAOUVwaYMy+PnXtrmPuTTpxwXPuYJrc5WCIxxhhg845qtu+u4fYfdyJ7cHy4w4koVm8zxrRrqoqIkDUgjhfmdict2a7OaiyrkRhj2i2fX7nr8XzeWlEGYEnkGFkiMca0S36/8ssnC3j/0wqKywJH38EcliUSY0y7Ewgov3m2gHdXlXPV+Wl89xSbU6QpLJEYY9oVVeXBv+2vnZjq+6enhDukiGeJxBjT7qQnebn0rBR+eJaN5Nsc7KotY0y7oKoUFPnplBbFf82wBNKcrEZijGnzVJVH/l7IFb/cQ94BHyJisxs2o5AmEhGZKiKbRGSLiMypZ/0kEVklIj4RuSCoPEtE/i0i60RkjYhcFLTuKRH5UkRWu4+sEJ2OMSYCqCp//kchC94uYXJ2Ap3S7BLf5haypi0R8QLzgClALrBCRBar6vqgzXYAs4Eb6+xeDlymql+ISHdgpYgsVdVCd/1NqrqgRU/AGBNxVJW/LCzk5bdKOO/kJK69MN1qIi0glH0kY4AtqroNQETmA+cCtYlEVbe76w65qFtVNwc9/1pE9gEZQGGLR22MiVivfVzGS8tKOPfkJK7/viWRlhLKRNID2Bm0nAuMbexBRGQMEIMzP/tB94rI7cBbwBxVrapnvyuBKwF69erV2Jc1xkSgyaMTqPYp505KsiTSgiKqs11EugHPAper6sFayy3AIGA00AG4ub59VfVRVc1W1eyMjIyQxGuMCT1V5e9vF1NaHiA2xsN5JydbEmlhoUwku4CeQcuZblmDiEgK8Cpwm6ouP1iuqrvVUQU8idOEZoxphwIB5Y8vHWDegkKWLi8NdzjtRigTyQqgv4j0FZEYYCawuCE7utsvBJ6p26nu1lIQ5yvHecDa5gzaGBMZAgHld3/bz6L3SrlwcjLnn2rDnoRKyBKJqvqA64ClwAbgJVVdJyJzRWQGgIiMFpFc4ELgLyKyzt39+8AkYHY9l/k+LyKfA58DnYB7QnVOxpjWwR9Q7n9uP69+VMYlU1O46vw0a84KIVHVI28gMlRV1x1xowiTnZ2tOTk54Q7DGNNM9hf7ueY3ezjnO0n88Gy7a72liMhKVc2uW96QGsmzQQf5cZ2DJjRDbMYYc0yqa5RAQOmQ4uXxW7tZEgmThiSS4PrhNXXWfdCMsRhjTINVVgf4xZ/zeGj+AQCSEiLqItQ2pSG/+eC2r7qNjvaXM8aEXGlFgJsfzmPlxkoG9Y0JdzjtXkNuSOwqIrOBz/h2IjlyB4sxxjSzwhI/Nz+8j227avjF5R05NTsx3CG1ew1JJHcCo4DLgUwRWY9z1dVGnKukjDEmJAIB5eaH9/HVHh93X5XBuBPiwx2SoQGJRFUfDV4WkUxgGDAceL+F4jLGmG/xeIQfzUgjNkYY0T8u3OEYV6PH2lLVXJxxsl5r/nCMMebbNu+o5suvqzlzXBJjhlotpLWxGRKNMa3ap5sq+b+/5JGa5OWUkQnExtg1Pq2NJRJjTKv17qpyfvVUPt0zovnN9RmWRFqpBv9VROR6EUlvyWCMMeagRe+VcPdf8xnYO5aHftaZjDT73ttaNSa9d8GZ1fAld8pcG8jGGNNiyisCjB8Wz/3XZ5CSaNPjtmZHHWvrkI2d5HEGzqXA2cBLwF9VdesRd2xlbKwtY1qnGp+Su6+Gvt1jUFUCCl6PfWdtLZoy1lYtdbLOHvfhA9KBBSLym2aJ0hjTbpVWBLhl3j7+58F9FJX6ERFLIhGiwY2OIvLfwGVAPvA4cJOq1oiIB/gC+HnLhGiMaev27fdx65/y+GpPDTdd2oHUJGvKiiSN6b3qAJyvql8FF6pqQESmNW9Yxpj2YvOOam57JI/KqgC/vq4zowbZjYaRpsFNW6p6R90kErRuQ0OO4XbSbxKRLSIyp571k0RklYj4ROSCoPIsEfm3iKwTkTUiclHQur4i8h/3mC+6sykaYyLEK++V4PXCH27sYkkkQjWmaetn9RQXAStVdXUD9vcC84ApOHfGrxCRxaq6PmizHcBs4MY6u5cDl6nqFyLSHVgpIktVtRC4D/idqs4XkT8DPwIeaeh5GWNCT1UprVCSEzz898wOlFYE6JBizVmRqjGd7dnAVUAP9/ETYCrwmIg0pH9kDLBFVbepajUwHzg3eANV3a6qa4BAnfLNqvqF+/xrYB+Q4V5FdhpwcB73p3HmbTfGtFI1PuW3L+zn+vv3UF4ZICZaLIlEuMYkkkxgpKr+r6r+L86IwJ1x51JvwP49gJ1By7luWaOIyBggBtgKdAQK3fngj/mYxpjQKCr1c/Mf97HkozImnphAXIxdldUWNKazvTNQFbRcA3RR1QoRqTrMPs1KRLrhTP07y+3kb8y+VwJXAvTq1atlAjTGHNaXX1fzi0fyyC/yM2dWR84Ya/OItBWNSSTPA/8RkVfc5enACyKSCKw//G61dgE9g5Yz3bIGEZEU4FXgNlVd7hYXAGkiEuXWSg57THc4/EfBuSGxoa9rjGkeD798gGof/P6GLgzuGxvucEwzalAicfsinsIZOn6CW3yVqh68PfySBhxmBdBfRPrifNjPBH7QwNePARYCz6jqwf4QVFVF5B3gApw+l1nAK/UfxRgTaoGAUlWjxMd6uGVWRxRszKw2qMFDpIjI56o6rEkvJnI28HvACzyhqveKyFwgR1UXi8honISRDlQCe1R1qIhcCjwJrAs63GxVXS0i/XCSSAfgU+BSVT1iU5sNkWJMyyurCPCrpwvw+5V7r87AY3epR7zDDZHSmK8Gq0RktKquONYgVHUJsKRO2e1Bz1fgNE/V3e854LnDHHMbzhVhxphWYvvuGu54NI9deT6u+V46NsRr29aYRDIWuFREtgNlgOC0Lg1vicCMMZHpnZwy7n9+P/ExwgM/7UzWALvJsK1rTCI5s8WiMMa0CRWVAf7090L6dY/mjis6WX9IO9GYv/IOnE71fqo6V0R6AV2BeodNMca0H/mFPtKTvcTHeXjwfzrTtWMU0VHWntVeNOaGxD8B44GL3eUSnCFPjDHt2H/WVfDje/fwzJIiAHp2ibYk0s40qo9EVUeKyKcAqnrABkg0pv3y+ZW/vlLIi8tK6NcjmtPH2A2G7VVjEkmNO/CiAohIBnXGxDLGtA97Cnzc/dd8NmyvZsbEJK7+XhqxMY2aJ8+0IY1JJH/Auceji4jci3MT4C9aJCpjTKtWXBZgT4GPO37ciZNHJoQ7HBNmDU4kqvq8iKwEJrtF5zV0HhJjTOQrrwzwwepyzhyXxIBeMbxwd3erhRigcfORxAIjgVR3vwtFBFWd21LBGWNah7Vbq/jV0wXsLfAxuE8svbpGWxIxtRrTtPUK7kRWHDoKsDGmjarxKU+/WsT8N4rp3MHL727oTK+u0eEOy7QyjUkkmao6tcUiMca0KqrKzX/cx+ovqjj7O4lcc0E6CXFWCzHf1phE8rGIDFPVz1ssGmNM2Pn9iscDIsL0iUl877RkJoywDnVzeI1JJCcBl4vINpymLRtry5g2Ztuuau5/dj9nfSeRGZOSOTXb7g0xR9eYRHJWi0VhjAmrGp/y/OtFPP96MckJHtKSbQ5103A21pYx7dzG7VX85tn9bN9dw+ljErj2gnRSkyyRmIZrTCL5E86d7KcBc3HG2vo7MLoF4jLGhEhpRYDyygC/vDqDccPiwx2OiUCNuQRjrKpeizNzIap6AGjUWFsiMlVENonIFhGZU8/6SSKySkR8InJBnXWvi0ihiPyrTvlTIvKliKx2H1mNicmY9kZVeWtFGfPfLAYge3A8z9zZ3ZKIOWYhG2vL3XceMAXIBVaIyGJVXR+02Q5gNnBjPYe4H0gAflLPupuC53I3xtTvq901/OGl/Xy6qYohfWO48LRkvF4hJtpG6zXH7ljG2up8jGNtjQG2uFPjIiLzgXOB2kSiqtvddd9KUKr6loic0ojXM8a4yioCPPtaEf94p4S4GOF/ZqZzzklJeG0eddMMjnWsLaHxY231AHYGLefiTN/bHO4VkduBt4A5qvqtO+9F5ErgSoBevXo108saExkKivwsfLeE08ckcsV5aaTbVVmmGTVqHkxV3QhsbKFYjtUtwB6c/ppHgZtxLgY4hKo+6q4nOztbQxmgMeHw2ReVrFhfyY/PTaNX12ien9udTjb1rWkBoXxX7QJ6Bi1numVNoqq73adVIvIk9fevGNNu5O6r4bFFhXywuoKMNC8XTk4mNclrScS0mFC+s1YA/UWkL04CmQn8oKkHFZFuqrpbRAQ4D1jb1GMaE4mKy/w8sbiIf31USky0cPn0VC6cnEycjdJrWljIEomq+kTkOmAp4AWeUNV1IjIXyFHVxSIyGqdDPx2YLiJ3qepQABH5ABgEJIlILvAjVV0KPO9eQSbAauCqUJ2TMa2BquJ8j4L3Py1n+klJ/PDsVDqkWD+ICQ1RbX/dBdnZ2ZqTkxPuMIxpkorKAAvfLSFnYyUP/LQzHo9QURkg3kboNS1ERFaqanbdcms0NSbCVFQGWPR+KS8tK6aoNMC4E+Ioq1SSE8SSiAkLSyTGRJBtu6r534f2UVQaYMyQOC47J5UhfWPDHZZp5yyRGNPKHSjxk7u3hmHHx9GrazTjTohnxsQkBlsCMa2EJRJjWqncfTW8vKyE15eXkprk5W93dyfKK9x8Wcdwh2bMISyRGNPKfLGzmmdeLeLjzyuI8sKZ45K4cLIzJpYxrZElEmNagcrqADU+SE7wUFjiZ+22Ki6ZmsJ5JyfbZbym1bNEYkwY7cqr4Z8flPLax2VMHZ/I1d9LJ3twHPPv6U6s3UhoIoQlEmPC4MPV5Sx6r4RVm6rweuCkrAQmZiUAICLExlgzlokclkiMCZGde2vo2SUagA9Wl5Ob5+Py6amcNT7RxsEyEc3evca0oMISP2/nlPPmf8rYtKOax27tynGZMVz//Q7Ex4nNB2LaBEskxrSAPQU+/vDiflasr8QfgON7RnPNBWlkpDsd50kJ1v9h2g5LJMY0g+oaZcX6CjweYfyweFITPeTu8/G905I5Y2wi/XrEhDtEY1qMJRJjjlF1jbJiQwXvryrn4zUVlFUqWQNiGT8snvg4D0/f0a12VF5j2jJLJMY0QnWNEhPtJIe7n8jno88qSE7wMPHEBE4dlcCJA+Nqt7UkYtoLSyTGHMWBEj//XlPBh5+V8+mmKp6b252Oqc7MgzMmJnHiwDii7K5z046FNJGIyFTgIZyJrR5X1V/XWT8J+D0wHJipqguC1r0OjAM+VNVpQeV9gflAR2Al8ENVrW7hUzHtwJad1fzxpQOs3VaFKnTu4OWcCYkEAs4cPsOPjzvKEYxpH0KWSETEC8wDpgC5wAoRWayq64M22wHMpv551+8HEoCf1Cm/D/idqs4XkT8DPwIeaebwTRvnDygbvqzmo8/KGdIvlolZCaQmeaioCvDDs1I4aUQCx2VGW3OVMfUIZY1kDLBFVbcBiMh84FygNpGo6nZ3XaDuzqr6loicElzmztN+Gt/M/f40cCeWSEwDqCoffVbBx59XsPzzCgpLA3g9cEmMMDErgYz0KB69tVu4wzSm1QtlIukB7AxazgXGNvGYHYFCVfUFHbNHfRuKyJXAlQC9evVq4suaSHWgxM/W3GqyB8cjIjzzWhF78n2MGRrPhOHxjB4aT1K83eNhTGO0m852VX0UeBScOdvDHI4JoT0FPj5YXc6Hn1WwbmsV0VHCovt7EBvjYe6VGXRK81pnuTFNEMpEsgvoGbSc6ZY1RQGQJiJRbq2kOY5p2gBVRUT4xzslPPzyAQD6dY/mUre/4+AlvF07tpvvUsa0mFD+F60A+rtXWe0CZvJN38YxUVUVkXeAC3Cu3JoFvNLUQE1k2lPg491V5by7spzLp6cydmg8owbFceV5aUzMiqdH5+hwh2hMmxSyRKKqPhG5DliKc/nvE6q6TkTmAjmqulhERgMLgXRguojcpapDAUTkA2AQkCQiucCPVHUpcDMwX0TuAT4F/hqqczLhV+NTli4v481Pyvh8SxUAA3vFcLChqne3aHp3swRiTEsS1fbXXZCdna05OTnhDsMcI79f2ZXno1fXaPwBZeZtX5MUL0wencip2Qn0yLDEYUxLEJGVqppdt9waiE3E2LWvhlc/LmPp8lI8Isy/pzter/CXOV1JT/HYPR7GhIklEtPqrdtWxZP/LGTVpio8Hhg7NJ6zxidysP2qQ6rNaW5MOFkiMa3SgRI/qtAhxUtVtdOUdfn0VKaOTyTDZhM0plWx/0jTqnz5dTV/f7uENz8pY8akZK69IJ0TB8by3NzuNpugMa2UJRLTKqzZUskLS4v5ZF0lsdHC1PFJTD8pCXCGY7f7BY1pvSyRmLA5eNMgwJKPyvhiRzX/NT2V6ROTSE2yfg9jIoUlEhNygYDy4WcVPPdaETf9sCP9e8Zw1flp3HBxOrExNs6VMZHGEokJGVXl4zUVPPmvIrbtqiGzcxSl5c5Az2nJVgMxJlJZIjEhoarc+Id9fLqpih4ZUdwyqyOnjU6wDnRj2gBLJKZFffl1NX26ORNCTRiewOmjEzljbCJe6z03ps2wRGJaxO58H4+/Usg7K8u5+yedmDAigfNPTQ53WMaYFmCJxDSr0vIAz71exMJ3S/CIcOlZKWQNsLnNjWnLLJGYZqOq/Oz3e9m6q4YzxyVy+fRUuwvdmHbA/stNk332RSVD+sYSHSVccV4aacle+veMCXdYxpgQsURijtne/T4e+fsB3v+0gv++KJ1zT05m9JD4cIdljAkxSySm0aprlBeXFfPC68UAXD7NGUzRGNM+hfQ2YhGZKiKbRGSLiMypZ/0kEVklIj4RuaDOulki8oX7mBVU/q57zNXuo3MozqU9++VT+Tz5zyLGnhDHk7d344dnp9od6ca0YyGrkYiIF5gHTAFygRUislhV1wdttgOYDdxYZ98OwB1ANqDASnffA+4ml6iqTXnYgvbu95EQ5yE5wcPFZ6Qw7aQksgdbM5YxJrQ1kjHAFlXdpqrVwHzg3OANVHW7qq4BAnX2PRN4U1X3u8njTWBqKIJu72p8yt/eKObyubt58p+FAAzsHWtJxBhTK5R9JD2AnUHLucDYJuzbI2j5SRHxA38H7tF6JqIXkSuBKwF69erViLDbr9WbK3lo/n6+2uNjwoh4vn96SrhDMsa0Qm2hs/0SVd0lIsk4ieSHwDN1N1LVR4FHAbKzs7+VaMyh/vFOCQ+/fIBuHb388uoMxg2zGkgkU1X8AQgEIMoLHhvjzDSjUCaSXUDPoOVMt6yh+55SZ993AVR1l/uzRERewGlC+1YiMUfnDyhlFQFSEr18Z3g8RaV+Lj4zhTjrSG/VanzKvz4sZenyMiqrAvgC4PMr1dVKVY1SXeMkkWBxMUJcjJCc6CEl0UNqkpeMNC+dO0TROd1L53QvGelRdErzEmXjopmjCGUiWQH0F5G+OIlhJvCDBu67FPiliKS7y2cAt4hIFJCmqvkiEg1MA5Y1c9ztwsbtVfx+/gFSkzz8+toMunaM4vLpaeEOyxyBP6C8vaKcp/5VyO4CPwN7x9C3RwxRXojyCjHRQmy08zPKC16v4BHn8u3KaqWySikuD1Bc5mdPgY+1W6soLjs044hAj4woBvSKoX9P53FcZrRNPGYOEbJEoqo+EbkOJyl4gSdUdZ2IzAVyVHWxiIwGFgLpwHQRuUtVh6rqfhG5GycZAcx1yxKBpW4S8eIkkcdCdU5tQXGZn7++UsS/PiolPdnD909PP/pOJqyqa5Sly0t5aVkJu/J8HN8zml/P7MDoIXG1M04eq/LKAPsO+Mk74CPvgJ+9+31s+7qGtduqeDunvHa7zuleTjguljFD4xk9JI50m0+mXZN6+qXbvOzsbM3JsauF122r4rZH8igtD3DeKcnMnpZKUrw1Y7VWewp8LPm4lFc/KuVAcYCBvWKYeUYKE7PiQ9LnUVTqZ0tuDVt2VrN5ZzWrN1VyoMSpwQzsFcOYoXGMPSGegb1jbJ6ZNkpEVqpqdt3yttDZbhqpukaJiRZ6d4tm+PGxzDonleMybWys1khVWbG+kn+8U8KKDZUAjBkSx/dPTyFrQGyTayCNkZrkZdQgL6MGOaM5BwLKltwa/rOugk/WVfD868U8+1oxndK8nDk2kanjE+nROTpk8ZnwsRpJO1JY4ueviwvZsL2aP8/pap2orZg/oHzwaTkvvFHMlp01dEz1cs6ERKaOT6Jrx9b5/a+4zM+K9ZUs+6SMFesrCSgMPz6W6ROTmJiVQEy0vd8indVI2jGfX1n0XglPv1pEZZXy3VOS8fnVEkkrVFDkd5qvPixl3wE/mZ2juOnSDpw+JpHoqNb990pJ9DJ5dCKTRyeSX+jjjf+UseTjMu59soCUxAOc/Z1ELpicQocU609pa6xG0sbtKfBx88P72LnXx+ghcVzzvXR6d7PmhtZmw/YqFrxdwvuryvEHYNSgOKZPTGLCiPiI7m8IBJTVm6tY/EEJH66uIDpKOOekJGZOSaaTzVUTcaxG0s6UVQRIjPeQkealV9dorjo/nXEnNP2qHtN8qqoDfPhZBa+8X8rarVUkxgnnn5rM9IlJZLaRvgWPRxg5KI6Rg+LYubeGF5YWs+i9Ev71YSnnn5LEzDNSSEm0GkqksxpJG7M738eT/ypk5cZKnr6ju12F1coEAsr6L6tZ9kkZb+WUUVahdOsUxfmnJHHWd5JIiGv7f6+v830882oRb35SRkKc8IMzU/neqcnH3Ifi9zs3XlbVKIGA8zsOKHg9zr0zUV4hPlasKbcZWI2kjdtf7OeF14tY/EEpHo/zzdb+bVoHVWXNlireySnnw8/K2V8cICZamHRiPFPHJ5HVP7ZdDVnSvVMUc2Z15PunJ/P4K4U8tqiQJR+Vcu0F6YwbFk8goOwv9rN3v5/8Qj95hT72F/k5UBLgQImfotIApeUBSisClFUE8Pkb9rrRUZAY597Jn+wlPdlDp7QounTw0qVDFD27RNEjI9ouCjgGViNpA/YU+Lh87m5q/MpZ4xO57BybK7012F/k541PyljyUSm5+3zExQhjhsYxMSuBcSfEk2i1RQBWrK9g3ssH2LHXR+cOXg4U+6nxHbpNdBSkJ3tJS/aSmuQhKcFDcryHhHgP8TFCbIxzB79HnOY0j0BAnQtNanzOnfwVlQHKKpWiUj+FblLKO+Cnsvqbz0CPQNdOURyfGc2AnjH07xXD4L6xVrN3Ha5GYokkQu0v8rN2WxWTTkwA4IWlRUw6MaHNtK1Hqv1Ffj78rJx3V5Wz5osqAgonHBfLORMSOXlkgo1bdhg1PuWV90vY8GU1ndO9dO0YRZcOUWSke+mU5iUl0dMi/XuqSkl5gD0FfnburWHn3hq2765hS24NX+c52cwj0K+Hc8/VyEFxnDgwjvjY9vl3tEQSJJITya68Gha8VcKSj0vxeoSXftXDvi2FkT+grN1axSfrKlmxvoItuTUA9OwSxamjEjh1VKJdJRehSssDbN5Rzedbq1izpZL126qpqlGio2D48XGcNCKeiScmtKvLmS2RBInERLJzbw2PLSrkozUVeD1wxrhELp6SYncOh0FxmZ/Vm6tYvraCf39eQVFpAK/HqXmMHhzHuGHx9O0ebVfItTHVNcrnW6v4ZJ3zd8/d58MjMLx/LFPGJHLKyATi2/jFEpZIgkRKIqmoDFBcHqBLhyh25dVw/f17OeekJM47OZmOqe3nW1A4qSpf5/vYuL2aTV9V89kXlWzJrUEVEuOF8SfEM2FEAtmD46zPox1RVbbvruG9VeW8nVPu9IHFCqeOTGD6xCQG9YkNd4gtwhJJkNacSFSdbz1vLC/j3VXlZA2I456rMgCnHbm1393cFhSW+Fm5sZL/rKtg5YZvBiaMiRYG94lh5ECnnXxQnxi7pNSgqqzbVs1rH5fyzqpyKquUof1iuOC0FE4aEY+3Db1HLJEEaalEUlUd4MEX9vNWTjkCeDyQEOehf88YBvaKYajb9HG4N9Yr75Xw0lsl7M53vt2cfKLz7WZI37b57aa1yC/0sXabU9v4bHMV23c7/RypSR5GD45jeP84BvWOoU/3aEsc5ojKKgK8/u9S/vFuKbvzfXTPiOKSqSlMGZPYJt47lkiCtEQi2V/k5//+ksfGr6o5Z0ISqYkeAuoMvb15RzXbvq4hEICOqV7O+k4ip45KYH+RnxUbKrl8WiqxMR6efa2IzzZXcsbYRCZmtf321nAoLPHzxc5qtu6qYWtuNeu2VbGnwLkRIS5WGHZcLCOOj2XU4Dj694xpV/d3mObjDygfr6ngudeK+GJnDd06RXHZ2SmcPiYxooe8aRWJRESmAg/hTEL1uKr+us76ScDvgeHATFVdELRuFvALd/EeVX3aLR8FPAXEA0uA/9ajnFRzJ5L1X1Yx9/F8issC3DK7IxOzEr61TVV1gLdzynjxzRJ27D30Ivnhx8cyeojTVDK4T2y7uLs5FCoqA7UJY/2XVaz7srr2kk6AjDQvg/vGMLRfLCccF0v/ntZUZZqXqrJ8bSVPv1rE5h3VHJfpDFd0cCj+SBP2RCIiXmAzMAXIxZnt8GJVXR+0TR8gBbgRWHwwkYhIByAHyAYUWAmMUtUDIvIJ8FPgPziJ5A+q+tqRYmmuRLJ5RzVP/auQ5Wsr6ZTm5d6rMzg+M5qi0gBf7alhx54atu2qYcKIeLIHx7Ppqyquvm8v3TOi6JTqXBvv8yu5+3zk7nOvWffA0H6xjB0ax7gT7OqfI6nxOTeX5Rc6j30H/HydV0NunvP73J3v4+DbOz3Fwwn9YhnSN5YBvWLo18OmizWhEwgo764q5/FXCtlT4GfcCXFc//0OdOsUWTcOt4YhUsYAW1R1mxvQfOBcoDaRqOp2d12gzr5nAm+q6n53/ZvAVBF5F0hR1eVu+TPAecARE8mx+nRTJa+8V8LeAz4KCv3kFwWIjoJTRsbzv5d0RIFpP8ulouqb5JwQJ/TqGk32YDi+ZwwLftWDDvVccVVSHmDj9irWbHEuL3z8lSIef6WI3t2imTLGuR+hcwdvRFeLj0VFZYCd+3xuUq5m264aduytoag0cMjv+aC4WCEzI4r+PWM4c2wi/TKjOa5HDF06eC0hm7DxeITTshM5aUQC/3i3hGeWFPFfd+/m0rNSOGdCEqUVAYrLApRXBqiuUap94PMpgYDiV+emyLgYD7ExQkKckJzgITnBQ0qit1UM6RLKRNID2Bm0nAuMbcK+PdxHbj3l3yIiVwJXAvTq1auBL3uol98qZvnaykPKorxC5w5RJMZ7UFWmnZRERrqX3l2j6dU1mow0b207u9cj9SYRgOQED6OHxDN6SDw/mpFGQZGfjz4r560V5bVJBZzElBjnvImS3DdTcuLBN1XQz0Rv7XJqooe4WGm1H6SBgLLvgJ/cfTXs3Otz7zB2kkde4TcDKUVHQa+u0QzpG0t6inN+qYkeOqZ5yUiLolOaM35Saz1PY2KihZlTUjhtVALzFhzgr4uL+OvioiYds0OKhy4doujaKYqenaPo2cX57OkbwotDIqte1QSq+ijwKDhNW8dyjP+Z2YGKKiUlyUN0lBAbLYdcjisiXP299GaJt2OqlxmTkpkxKZmv8318ss658a204psB60rKAuzO97F5R4CS8sAhYwbVFeWF5EQPKUEJqPZnvIfEg484ITHeQ3ysh/hYIS5oHKMY93yPVivy+ZWqaudRUeXGXKGUlPnZXxygsNTP/iJnUL59B3zs3e+nuubbtbisAbH06hpd+4+R2TnK+jBMm9C5QxR3XZnByo2VbNtVTWqih5QkL4lxQmyMh+go50uqx+N8AQ0oVFY5/+Pllc5nQEl5gAMlAfYW+Ni738fGL6t4d2V5bXNudBQcnxnD4D4xDDs+juH9Y0lPbpnm3FAmkl1Az6DlTLesofueUmffd93yzGM8ZqN17hCevNu9UxTnnZx81O2qa5xxg4rL/BSXOYmm2H2UuG+84jI/JeUB8ov8bN9dQ2m5M5BdY7rKRJzEJCKIu6zqXKniD3DUY3kE0pKdb1H9esQwfpiXzM5OosjsHEXHVGuGMu3DqEFxzdrxXl2j7MqrYfvXNWz8qpqNX1Xz6kdl/OPdUgB6d43ijisy6NPMw/aE8pNxBdBfRPrifNjPBH7QwH2XAr8UkYNf988AblHV/SJSLCLjcDrbLwP+2MxxR4yYaKFjqrfRd70HAkp5pfNNp6wyQFmFUlnt9EFUVKnTZuvO91DjU3x+xecHDSjBOSPKK3g9EBXl1mSihbhYD0nxUlv7SU/2kpzoaXd9PcaEQky00Ld7DH27x3BqdiLgXJSyeUc1a7ZUseaLSjqnN3+tJGSJRFV9InIdTlLwAk+o6joRmQvkqOpiERkNLATSgekicpeqDnUTxt04yQhg7sGOd+Aavrn89zVaqKO9LfN4hKQE58PeGNO2REcJQ/vFMrRfLBefkdIir2E3JBpjjGmQw13+a19BjTHGNIklEmOMMU1iicQYY0yTWCIxxhjTJJZIjDHGNIklEmOMMU1iicQYY0yTtMv7SEQkD/jqGHfvBOQ3YziRwM65fbBzbvuaer69VTWjbmG7TCRNISI59d2Q05bZObcPds5tX0udrzVtGWOMaRJLJMYYY5rEEknjPRruAMLAzrl9sHNu+1rkfK2PxBhjTJNYjcQYY0yTWCIxxhjTJJZIDkNEporIJhHZIiJz6lkfKyIvuuv/IyJ9whBms2rAOf9MRNaLyBoReUtEeocjzuZ0tHMO2u57IqIiEtGXijbkfEXk++7feZ2IvBDqGJtbA97XvUTkHRH51H1vnx2OOJuTiDwhIvtEZO1h1ouI/MH9nawRkZFNekFVtUedB84MjluBfkAM8BkwpM421wB/dp/PBF4Md9whOOdTgQT3+dXt4Zzd7ZKB94HlQHa4427hv3F/4FMg3V3uHO64Q3DOjwJXu8+HANvDHXcznPckYCSw9jDrz8aZTVaAccB/mvJ6ViOp3xhgi6puU9VqYD5wbp1tzgWedp8vACaLSCRPRH7Uc1bVd1S13F1cDmSGOMbm1pC/M8DdwH1AZSiDawENOd8rgHmqegBAVfeFOMbm1pBzVuDgHLSpwNchjK9FqOr7wP4jbHIu8Iw6lgNpItLtWF/PEkn9egA7g5Zz3bJ6t1FVH1AEdAxJdC2jIecc7Ec432gi2VHP2a3y91TVV0MZWAtpyN94ADBARD4SkeUiMjVk0bWMhpzzncClIpILLAGuD01oYdXY//cjimpyOKbdEZFLgWzg5HDH0pJExAM8CMwOcyihFIXTvHUKTo3zfREZpqqF4QyqhV0MPKWqvxWR8cCzInKCqgbCHViksBpJ/XYBPYOWM92yercRkSicKnFBSKJrGQ05Z0TkdOA2YIaqVoUotpZytHNOBk4A3hWR7ThtyYsjuMO9IX/jXGCxqtao6pfAZpzEEqkacs4/Al4CUNV/A3E4gxu2ZQ36f28oSyT1WwH0F5G+IhKD05m+uM42i4FZ7vMLgLfV7cWKUEc9ZxE5EfgLThKJ9LZzOMo5q2qRqnZS1T6q2genX2iGquaEJ9wma8j7ehFObQQR6YTT1LUthDE2t4ac8w5gMoCIDMZJJHkhjTL0FgOXuVdvjQOKVHX3sR7Mmrbqoao+EbkOWIpz1ccTqrpOROYCOaq6GPgrThV4C06n1szwRdx0DTzn+4Ek4GX3uoIdqjojbEE3UQPPuc1o4PkuBc4QkfWAH7hJVSO2pt3Ac/5f4DERuQGn4312hH8pRET+hvOFoJPb93MHEA2gqn/G6Qs6G9gClAOXN+n1Ivz3ZYwxJsysacsYY0yTWCIxxhjTJJZIjDHGNIklEmOMMU1iicQYY0yTWCIxxhjTJJZIjDHGNIklEmOamYikicg1Qcsft9DrZIrIRS1xbGMawxKJMc0vDWe+GgBU9Tst9DqTceacMCasLJEY0/x+DRwnIqtF5H4RKQUQkT4islFEnhKRzSLyvIic7g7Z/oWIjDl4ABG5VEQ+cY/xFxHxBr+AiJyEMzLxBe42/UJ6hsYEsURiTPObA2xV1SxVvanOuuOB3wKD3McPgJOAG4FboXbgwIuACaqahTPm1SXBB1HVD3EGJDzXfZ1IHljRRDgbtNGY0PpSVT8HEJF1wFuqqiLyOdDH3WYyMApY4Q6OGQ/UN9ryQGBji0dszFFYIjEmtILncAkELQf45v9RgKdV9ZbDHcQd4r3InZ3TmLCypi1jml8JzqRYx+otnL6PzgAi0kFEetfZpg9tYG5x0zZYIjGmmbnzd3wkImtF5P5j2H898AvgDRFZA7wJdKuz2UacuSbWikhLXRVmTIPYfCTGGGOaxGokxhhjmsQSiTHGmCaxRGKMMaZJLJEYY4xpEkskxhhjmsQSiTHGmCaxRGKMMaZJ/h8AAc5lmlnp8AAAAABJRU5ErkJggg==\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": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Magnetization over time')"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEYCAYAAABY7FHWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA95ElEQVR4nO3deZyN9fvH8dc1m7GPZUjGnj1GtiihRJKUb8nSQpFC2hSShPhWPxJFm5Jkj5RKlCJJliFkJ0sGMQbDMMMs1++Pc/jOTIM5zDn3LNfz8TgP517Ofd73nDHXue/PfX8+oqoYY4wx5/k5HcAYY0zWYoXBGGNMKlYYjDHGpGKFwRhjTCpWGIwxxqRihcEYY0wqVhhMriIi34tIVy9sd7OINM/s7WYlIvKgiPzgdA7jfWL3MZgrJSJ7gWuBa1X1aIr5fwB1gAqquteRcK4cQ4HrVPWhTN7uZCBSVQdn5nazEhEpD+wBAlU10eE4xsfsiMFcrT1A5/MTIlILyOdcHOMpEQlwOoPJWqwwmKv1OfBIiumuwJSUK4jIXSLyh4icFJH97m/yKZc/IiL7RCRaRF4Rkb0icrt72VARmS0iU0TklPuUTf0Ur71WROaKSJSI7BGRp93zWwODgI4iEisiG9zzl4pID/fzDe5l5x96/nSQiHwhIv+ISIyILBORmu75PYEHgf7u13zjnp8ycx4RGSsiB92PsSKSx72suYhEikg/ETkiIodE5NGL/XDd+zdfRI6JyC4ReTzF/DgRKZpi3RtE5KiIBLqnHxORrSJyXEQWiUi5FOuqiPQRkZ3AznTeepn73xPu/WwsIt1EZHmabfQWkZ3uz+Y1EakkIivcn/VsEQlKsX5bEVkvIifc69S+2H4bh6mqPexxRQ9gL3A7sB2oDvgDkUA5QIHy7vWaA7VwfRGpDRwG7nUvqwHEAk2AIGA0kADc7l4+FIgH2ri3/zqw0r3MD1gLDHG/tiKwG7gjxWunpsm8FOiRzr70BLYBhdzTjwEFgTzAWGB9inUnAyPS+1m4nw8HVgIlgFBgBfBaip9FonudQPd+nQGKXORnvAx4DwjGdXouCrjNvexn4PEU644CPnA/vwfY5f5cAoDBwIoU6yrwI1AUyJvO+5Z3rxOQYl43YHmabXwNFAJqAmeBn9yfQ2FgC9DVve4NwBHgRvfn2NX9M8vj9O+xPf79sCMGkxnOHzW0BLYCB1IuVNWlqvqnqiar6kZgBtDMvfh+4BtVXa6q53D9kU/b8LVcVReoapL7vcLd8xsAoao6XFXPqepuYCLQyZPwItIEGAG0U9WT7syTVPWUqp7FVWDCRaRwBjf5IDBcVY+oahQwDHg4xfIE9/IEVV2AqzBWTSdXGeBmYICqxqvqeuBj/neENh33aTwREVz7Pd297EngdVXdqq42gv8CdVIeNbiXH1PVuAzuV3r+T1VPqupmYBPwg6ruVtUY4HtcBQFchfdDVV2lqkmq+hmuQtLoKt7beIkVBpMZPge64PpGOSXtQhG5UUSWuE/3xOD6o1XcvfhaYP/5dVX1DBCdZhP/pHh+Bgh2nxcvB1zrPjVxQkRO4Dp9VDKjwd1/fGfj+ma7wz3PX0TeEJG/ROQkrm+2pMh8OdcC+1JM73PPOy9aUzfongEKXGQ7x1T1VJptlXY/nws0FpFSQFMgGfjVvawcMC7Fz+UYICleCyl+7lfhcIrncelMn9+vckC/NJ9VGVL/XEwWYY1O5qqp6j4R2YPrtEj3dFaZDowH7lTVeBEZy//+yB4ixbdlEckLFMvgW+8H9qhq5YtFu9SL3e/1FTBWVb9PsagLrlMxt+MqCoWB47j+sF52u8BBXH8IN7uny7rneeogUFRECqYoDmVxH5Gp6nFxXT7aEdcpo5mqej7bfmCkqk67xPYvtR+Zfbni+TwjM3m7xgvsiMFklu64zn2fTmdZQVzffONFpCGuP7znzQHuFpGb3A2VQ/nfH+DLWQ2cEpEBIpLX/U3/ehFp4F5+GCgvIhf7PZ8EbFPV/0sn71lcRy75cJ2GSekwrvPoFzMDGCwioSJSHNfpsakZ3KcLVHU/rvaJ10Uk2N1Y2z3NtqbjOrV0P/87jQTwAfBSikbzwiLSwYO3j8J1BHKp/fTEROBJ99GjiEh+cV2UUDCTtm8ykRUGkylU9S9VjbjI4t7AcBE5heuP5OwUr9sM9AVm4jp6iMXVSHk2A++ZBLTF1Si7BziK6xz8+baAL9z/RovIunQ20Qlon+bKpFtwnQ7bh+ub+RZcDckpfQLUcJ8S+Sqd7Y4AIoCNwJ/AOve8K9EZV0PwQWAe8KqqLk6xfD5QGfhHVTecn6mq84A3gZnu02GbgDsz+qbuU3ojgd/c+3lVbQHu343HcR05HsfVMN7tarZpvMducDNZiogUAE4AlVV1j8NxjMmV7IjBOE5E7haRfCKSH9flqn/yvwZfY4yPWWEwWcE9uE6VHMR1WqST2qGsMY6xU0nGGGNSsSMGY4wxqWT7+xiKFy+u5cuXdzqGMcZkK2vXrj2qqqHpLcv2haF8+fJERFzsKkljjDHpEZF9F1tmp5KMMcakYoXBGGNMKlYYjDHGpJLt2xiMMVcuISGByMhI4uPjnY5ivCQ4OJiwsDACAwMz/BorDMbkYpGRkRQsWJDy5cvjGtLB5CSqSnR0NJGRkVSoUCHDr/PZqSQRmeQeynDTRZaLiLzjHr5wo4jU9VU2Y3Kr+Ph4ihUrZkUhhxIRihUr5vERoS/bGCYDrS+x/E5c3SFUxjXa0/s+yGRMrmdFIWe7ks/XZ6eSVHWZiJS/xCr3AFPcfeSsFJEQESmlqoe8kSdu5UY+ffEV4s6cpkuVuki+PKw9e4LAa0Jp9mgXAsuWwq9gfm+8tTHGZGlZqY2hNKmHGox0z/tXYRCRnriOKihbtuwVvVl8xCZmrfqVGE3g3r9dXf+/FvMnSarM/OxHAHqf2kZB/0DeqtaEwMrl2FuuCDW6PUBI9euu6D2NMan5+/tTq1YtVBV/f3/Gjx/PTTfdxN69e6levTpVq/5vKOzVq1cTFBTEwoULGTJkCCdPniQ4OJiqVasyatQoypYty8qVK3nmmWc4e/YsZ8+epWPHjgwdOtS5HcymslJhyDBV/Qj4CKB+/fpX1Atgkae68HOvjiTFnsHvbAJJx0/y0e+riYuKpmTlqiT8fYhq496iwMl4EiMPk/j3IR44tobaY9/mo2btCb6xFmtPHKHx049RqHLGG3WMMf+TN29e1q9fD8CiRYt46aWX+OWXXwCoVKnShWXnbdq0ib59+zJ//nyqV68OwPz589m7dy9ly5ala9euzJ49m/DwcJKSkti+fbsvdyfHyEqF4QCuwcHPC3PP8xp/f3/8C7tGFgwoUZTaVcunWj7u6QcBV8t+wp4DjBn1LvkOn0CSgoj85AvuOvQrPT+ZzAvX1CT45hvQrm0pffvN3oxsTI518uRJihQpcsl13nzzTQYNGnShKAC0a9fuwvMjR45QqlQpwPX/u0aNGt4Jm8NlpcIwH3hKRGYCNwIx3mpf8JSIEFQxjAfef/PCvKLHY5j56v9Rcv1udPcR1n67kPZT32Bitebc/9pL5L+rGX558ziY2hjPHH35Hc5u2pmp28xzfWWKj3z6osvj4uKoU6cO8fHxHDp0iJ9//vnCsr/++os6deoAcPPNNzNhwgQ2b97MCy+8cNHtPffcc1StWpXmzZvTunVrunbtSnBwcKbtT27hy8tVZwC/A1VFJFJEuovIkyLypHuVBcBuXGPBTsQ1TnCWlbdIYe5/ZyS3LJtBxf2LqTzhVXpVrU+4f36O9HqNT8s24qHQiux993M0OdnpuMZkSedPJW3bto2FCxfyyCOPcH6MmPOnktavX8+ECRP+9dro6Gjq1KlDlSpVGD16NABDhgwhIiKCVq1aMX36dFq3vtSFkOZifHlVUufLLFegj4/iZCoRoVqXe3m3y71ocjLxK9bz9+NPsWLnduKHfcCeUZ8R06oB149+icCQQk7HNSZdl/pm7wuNGzfm6NGjREVFXXSdmjVrsm7dOsLDwylWrBjr169n9OjRxMbGXlinUqVK9OrVi8cff5zQ0FCio6MpVqyYL3Yhx7C+kjKZ+PmRt0ldhm1dwZZ9ewj5TyuSExPpMHkc7Updx6nZC9HERKdjGpPlbNu2jaSkpEv+Ee/fvz8jR45k69atF+adOXPmwvPvvvvuwhHHzp078ff3JyQkxGuZc6qs1MaQ4wSXLsk1H71K0rmXGPjYMwSt2sSRPiM59PrHLC6YSPe5n5An1L7JmNzrfBsDuC7y+Oyzz/D397/o+rVq1WLcuHE88sgjnDx5kuLFi1O2bFmGDRsGwOeff85zzz1Hvnz5CAgIYNq0aZfcnklfth/zuX79+ppdBurR5GTO/LCCSX368/Tu35kSUotWXR6g+P/1wz9/XqfjmVxo69atqa7wMTlTep+ziKxV1frprW+nknxI/PzI37oJfXYuZ17fQdxUsASxsxfxUZmGLHl+GNm9SBtjcgYrDA7w8/Pj3ndGUnH3Qgr2eoC3Tu7mlfFjOfLsGyRGHXc6njEml7PC4CAJCKDE8L6s+msr7zzah9jZi9jU4H7ebnAHCcdjnI5njMmlrDBkAcXKlaXuhyMo88tk5nKcFyJ+YHGVlpyYOMfpaMaYXMgKQxYSVKU8w3ZHsLDXS1Txz0v0oHF8eUMb4g8cdjqaMSYXscKQxfj5+dHyvf9Sdu0XHCpbjI7rFzIg/BbO7dzndDRjTC5hhSGLCixdgpvWfMlHjz5F14JliGzRnf2jPiHxZOzlX2xMNnL48GG6dOlCxYoVqVevHo0bN2bevHlXtc2hQ4em6iZj8eLFV7Sd9evXs2DBggyt27x5c5y8dP6rr75iy5YtmbItKwxZmJ+fH49Oeofav88isMH1dHn5Bf5T6jrOLF3jdDRjMoWqcu+999K0aVN2797N2rVrmTlzJpGRkf9aN/EKewwYPnw4t99++xW91pPC4DQrDLlMwDXFuXbWaO64pSmN/ApyqMPzRL08zu57MNnezz//TFBQEE8++eSFeeXKlaNv374ATJ48mXbt2nHbbbfRokULYmNjadGiBXXr1qVWrVp8/fXXF143cuRIqlSpQpMmTVKNw9CtWzfmzHFdyLF27VqaNWtGvXr1uOOOOzh0yNWBc/PmzRkwYAANGzakSpUq/Prrr5w7d44hQ4Ywa9Ys6tSpw6xZs1Jlj4uLo1OnTlSvXp327dsTFxd3YdkPP/xA48aNqVu3Lh06dLjQl9PAgQOpUaMGtWvXvtBL7OHDh2nfvj3h4eGEh4ezYsUKAKZOnUrDhg2pU6cOTzzxBElJSQAUKFCAl19+mfDwcBo1asThw4dZsWIF8+fP58UXX6ROnTr89ddfV/W5WJcY2YR/QACDf/qauOXrOPTgABaMn8j26VMYtvYnAouGOB3P5BAH7un7r3kF7rmNwo+1J/lMPIc6v/iv5QU73Umhzm1Iij7BP4+9kmpZ6a/fveT7bd68mbp1615ynXXr1rFx40aKFi1KYmIi8+bNo1ChQhw9epRGjRrRrl071q1bx8yZM1m/fj2JiYnUrVuXevXqpdpOQkICffv25euvvyY0NJRZs2bx8ssvM2nSJMB1RLJ69WoWLFjAsGHDWLx4McOHDyciIoLx48f/K9f7779Pvnz52Lp1Kxs3brywH0ePHmXEiBEsXryY/Pnz8+abbzJmzBj69OnDvHnz2LZtGyLCiRMnAHj66adp1qwZ8+bNIykpidjYWLZu3cqsWbP47bffCAwMpHfv3kybNo1HHnmE06dP06hRI0aOHEn//v2ZOHEigwcPpl27drRt25b777//kj/PjLDCkM3kbVKXshu/ZHH1+qzav4fuTwyj/Cev4V+ogNPRjLlqffr0Yfny5QQFBbFmjeuUacuWLSlatCjgOvU0aNAgli1bhp+fHwcOHODw4cP8+uuvtG/fnnz58gGpB+85b/v27WzatImWLVsCkJSUdGFQH4D//Oc/ANSrV4+9e/deNuuyZct4+mlXj7S1a9emdu3aAKxcuZItW7Zw882uQbvOnTtH48aNKVy4MMHBwXTv3p22bdvStm1bwHXUNGXKFMA1uFDhwoX5/PPPWbt2LQ0aNABcRyclSpQAICgo6MJr69Wrx48//pihn60nrDBkQwGFC/JZ5Fb2vjUJfWsqe1p0J7ZDc+r07+V0NJPNXeobvl++4Esu9y8WctkjhLRq1qzJ3LlzL0xPmDCBo0ePUr/+/7rwyZ8//4Xn06ZNIyoqirVr1xIYGEj58uWJj4/P0HupKjVr1uT3339Pd3mePK6Btfz9/a+4PeP8+7Rs2ZIZM2b8a9nq1av56aefmDNnDuPHj081MFHabXTt2pXXX3/9X8sCAwMRkUzJejHWxpBN+fn5UfHFHlw7923G7lhDk4F9WdehL+o+D2lMdnDbbbcRHx/P+++/f2Feym6004qJiaFEiRIEBgayZMkS9u1zXcbdtGlTvvrqK+Li4jh16hTffPPNv15btWpVoqKiLhSGhIQENm/efMl8BQsW5NSpU+kua9q0KdOnTwdcY1Fv3LgRgEaNGvHbb7+xa9cuAE6fPs2OHTuIjY0lJiaGNm3a8Pbbb7NhwwYAWrRocWH/k5KSiImJoUWLFsyZM4cjR44AcOzYsQv7eiVZPWWFIZvLe/MN9Fs4i36h1Si8dD37b+lK0qnTTscyJkNEhK+++opffvmFChUq0LBhQ7p27cqbb76Z7voPPvggERER1KpViylTplCtWjUA6tatS8eOHQkPD+fOO++8cAompaCgIObMmcOAAQMIDw+nTp06Fxp6L+bWW29ly5Yt6TY+9+rVi9jYWKpXr86QIUMutGmEhoYyefJkOnfuTO3atWncuDHbtm3j1KlTtG3bltq1a9OkSRPGjBkDwLhx41iyZAm1atWiXr16bNmyhRo1ajBixAhatWpF7dq1admy5YWG8ovp1KkTo0aN4oYbbrjqxmfrdjuHSI47y4G2vdn5xwY+S4pi/IrFFKxZ2elYJouzbrdzB+t2O5fyy5uHsMUfs7Z2ab45dZA/e75C0vGTTscyxmRDVhhyEBHhxaVf88cnMyi59wiRdz7JgbkLnY5ljMlmrDDkQOUevZ9r57zNzB3rCe/QjtXPDnc6kjEmG7HCkEPlbRxO60mjaZWvBIWn/sDhXq/ZndLGmAyxwpCDhd/bhs+2ryOoSCH+mb2AqU3utctZjTGXZYUhhwssXYIya2byvv8xuq/4ho0vvOF0JGNMFmeFIRcICCnEW9tW8fnN91Jg6kKOj/3cTiuZLMHf3586depceLzxRuZ9cclOPaNmNdYlRi6Rv2gIDyyZzZG+/2XR0FGsfPN1Rm9bRUDhgk5HM7lY3rx5Wb9+vVe2vX79eiIiImjTpo1Xtp+T2RFDLiKBAZSY8DJLQgP4NmovWxs+QKLd62CymJiYGKpWrXqh6+zOnTszceJEwHW3cf369alZsyavvvrqhdesWbOGm266ifDwcBo2bEhMTMwlu8w2l6Gq2fpRr149NZ5JSkrSLR2e0V3Fm+hfldvo2cNHnY5kHLJly5ZU082aNdNPP/1UVVXPnTunzZo1088//1xVVU+fPq3NmjXTmTNnqqrqiRMntFmzZjp37lxVVY2KitJmzZrp/PnzVVX10KFDl31/Pz8/DQ8Pv/A4v+0ffvhBGzVqpDNmzNA77rjjwvrR0dGqqpqYmKjNmjXTDRs26NmzZ7VChQq6evVqVVWNiYnRhIQE/fTTT7VPnz5X+qPJUdJ+zqqqQIRe5O+qnUrKhfz8/Kg2623+eWwwL838lLOVb2DK35sILFzI6Wgml7nYqaSWLVvyxRdf0KdPnwudzQHMnj2bjz76iMTERA4dOsSWLVsQEUqVKnWhf6RChez3+GpZYcilRISSn7xGyPpfObntL6IHj6fkuIGIn51dzM2WLl164XlgYGCq6Xz58qWaLly4cKrp4sWLp5q+5pprrjhHcnIyW7duJV++fBw/fpywsDD27NnD6NGjWbNmDUWKFKFbt24Z7nLbeMb+CuRifn5+jFm7lFGvDuP0zO/Z1nUA545EOx3LGN5++22qV6/O9OnTefTRR0lISODkyZPkz5+fwoULc/jwYb7//nvA1Z32oUOHLgzsc+rUKRITEzO1G+rcxqeFQURai8h2EdklIgPTWV5WRJaIyB8islFE7HICL/Pz86PogO7I4+25a8YE+lRrYA3Sxmfi4uJSXa46cOBAtm/fzscff8xbb73FLbfcQtOmTRkxYgTh4eHccMMNVKtWjS5dulwYIS0oKIhZs2bRt29fwsPDadmyJfHx8ZfsMttcms+63RYRf2AH0BKIBNYAnVV1S4p1PgL+UNX3RaQGsEBVy19qu9btduZITk7mxXrNuHHPSRqWKkeZldPwt0tZczzrdjt3yMrdbjcEdqnqblU9B8wE7kmzjgLnW44KAwd9mC9X8/PzY/S6ZdzSrg1JR4+zuEF7kk/HOR3LGOMAXxaG0sD+FNOR7nkpDQUeEpFIYAHQN70NiUhPEYkQkYioqChvZM2VRIRrPvsvW2qF0WbnT4y7ua3TkYwxDshqjc+dgcmqGga0AT4XkX9lVNWPVLW+qtYPDQ31ecicTERo88PnDLq+CS0i4zjx4RdORzJe5qvTycYZV/L5+rIwHADKpJgOc89LqTswG0BVfweCgeI+SWcu8A8IYPgfSyjR9lYOvTyWL+98GE1OdjqW8YLg4GCio6OtOORQqkp0dDTBwcEevc6X9zGsASqLSAVcBaET0CXNOn8DLYDJIlIdV2Gwc0UOkIAASn74Km9W/oZ3F05lWbs8NP5mIiLidDSTicLCwoiMjMROyeZcwcHBhIWFefQanxUGVU0UkaeARYA/MElVN4vIcFy3Zs8H+gETReQ5XA3R3dS+yjhG8gQxdNWPVLnxbkqu2s6RPiMp+d5gp2OZTBQYGEiFChWcjmGyGJ9druotdrmq9yVGn+DvGzuzM/oIeTq2ovkHrzsdyRhzlbLK5aommwooFkLYkk95/vQOen/yLvGbdzkdyRjjRVYYTIYElbmG6XO+4L3S9TncbTCJR445HckY4yVWGEyGhbdrTaM5E0j8J4p3bmhB1NJVTkcyxnhBhgqDiBQRkUdEZJ6IbBaRb0XkcREp4e2AJmsJrl+TY/260P/geka160LCnrRXHBtjsrvLFgYR+RL4EigBDFDVmkBvXJeSThWRpV5NaLKcBs9255teA+kRVIrIlj2s0z1jcpjLXpUkIiGqeuJKl3ubXZXknMO9X+PArAUsK5jEM9uXI3mCnI5kjMmgq7oq6XJ/9J0sCsZZJSYM5rOiCfTfF8GGwWOcjmOMySQeNz6LSEsRmSgiddzTPTM9lckWRIT/W7+Mr5vcR8Fpi4hbvs7pSMaYTHAlVyU9BryIqxfU24A6mZrIZCt58uej1beTCKwYxpz2PVg18A2nIxljrtKVFIZTqnpCVV8AWgENMjmTyWb8Cxck5JNhvHhsMy+/PYpT85c4HckYcxWupDB8d/6Jqg4EpmReHJNdFapakfmfTmF0waoc6TmU+E07nY5kjLlCGblc9aiI3HV+WlW/TrlcVd/1RjCT/TR8pCMVRvcnMTGJCc3bkxBld0cbkx1l5IghAXhfRHqkXSAiMzI/ksnOCne7l58bleOF6M3M7dzb6TjGmCuQkcJwCGgKvCAiQ9Msq5rpiUy298T8z/mi9UM02HCAU/MWOx3HGOOhDLUxqOpeoAnQSkQ+Tm+4TWPO8/Pz476vJxHcOJw1Tw5mae+XnY5kjPFARv7AC4CqHsU1uloJYL6I5Du/zJi0JCiQkp8M5+njW+j10VhOTP/u8i8yxmQJGSkMf5x/oqpxwL24huZcChTySiqTIwSEFmXytKmML1yD6OfeJP6PrU5HMsZkQEa6xHgszXSyqj6B67JVGxPQXFL9Du1oMGEYJCtzWz9EwtHjTkcyxlzGFbcVqOowrDCYDCjUqQ3rbqvBQ0fW8dGdnZ2OY4y5jKtqRFbVfSIyILPCmJzrvhnvM7Zha27fE8vphcudjmOMuYQr6URvdorHF8C/7m8wJi0/Pz/6Lv2K/DdUZ1ePV1g/crzTkYwxFxFwBa85qaoXioGIvJ+JeUwO5pc3D9dM+S8PVKjKkVdXsbZBHQq3auJ0LGNMGh53iQGMTLOKXaRuMiygVCjD33idV/NX4mi3V0jYf8jpSMaYNDzuEkNV95xfICIzVNU6xDEeuf3ZJ7hr8POQkEhEi64kx591OpIxJgXrEsM4ouiLj/JnrdI037mYya27OB3HGJOCdYlhHHPnN5N49NpqhO+M4uymXU7HMca4WZcYxjF58udj/B/LKVasOAceGsjRJSudjmSMwbrEMA4LKFGUkp++Rp9tv3DfnXdxdk+k05GMyfWsSwzjuLz1atK63d00Dwjh4F29rTHaGIdl5HLVdNc53yWGtTeYzPDM7Ek80eZekqOOc6Dzi07HMSZXy8gf9R9FZJaIdBaRQgAikk9E/gOMANZl9M1EpLWIbBeRXSIy8CLrPCAiW0Rks4hMz+i2TfZX6vPX2VzYn2bffMzK/v91Oo4xuVZGTiW1AIYB5YHvRGQl8BNQG3hbVetk5I1ExB+YANwJ1AA6i0iNNOtUBl4CblbVmsCzGd0Rk/1JniCun/4WRf2DiP3mF5JOxjodyZhcKaOXq25R1ddV9RbgVlVtrKpDVTXDRwtAQ2CXqu5W1XPATOCeNOs8DkxQ1ePu9z3iwfZNDlCxYT2WLVlKhRNnOdxnBEmnTjsdyZhcx+P2AfeVSVeiNLA/xXSke15KVYAqIvKbiKwUkdZX+F4mG8t3Ux2KDOnFyLnTGVL9JjQx0elIxuQqWa3hOACoDDQHOgMTRSQk7Uoi0lNEIkQkIioqyrcJjU+EPNGBfYUC2BN9hH8eH+p0HGNyFV8WhgNAmRTTYe55KUUC81U1wd0n0w5chSIVVf1IVeurav3Q0FCvBTbO8fPzY+629YwsU5cz3/5CzKR5TkcyJtfwZWFYA1QWkQoiEgR0AuanWecrXEcLiEhxXKeWdvswo8lC8oQU4tr54znIOXr16cPJdZucjmRMrpDhwiAieUSki4gMEpEh5x8Zfb2qJgJPAYuArcBsVd0sIsNFpJ17tUVAtIhsAZYAL6pqdMZ3x+Q0eapV4O8Hb+fbs0dY+fJop+MYkyuIqmZsRZGFQAywFkg6P19V3/JOtIypX7++RkREOBnB+MCOgaPw/2Q+oeMGUqjLXZd/gTHmkkRkrarWT2+ZJyO4hamqXSVkHFF55PMc3BHJ5N4vEr7sN5p9YDfAGeMtnrQxrBCRWl5LYswliL8/+d/qx8hTf/H25I85/ePvTkcyJsfypDA0Ada6u7TYKCJ/ishGbwUzJq2QCmX54bPpvF6gMocffZmEf+xSZWO8wZNTSXd6LYUxGRT+0H1E7z3KobGfMa3Jf+i6bQkS4MmvsTHmcjJ8xKCq+4AQ4G73I8Q9zxifKvpyT94rdIaef/3Gpjc+cDqOMTmOJ5erPgNMwzWCWwlgqoj09VYwYy5GRHj9t0VMrtqMAlMXknjQutQyJjN50sbQHbhRVYeo6hCgEa5O74zxucKlSnL/N5PR+HP80qY7p//c4XQkY3IMTwqDkOL+BfdzG/PZOCaocjlO9+vCXRsWMPjWu0k+E+90JGNyBE8Kw6fAKhEZKiJDgZXAJ15JZUwG1erbjaENW9JZQzjU6QUyesOmMebiPGl8HgM8BhxzPx5V1bFeymVMhg349RtKhJUmbsV6dg8f73QcY7K9DHeJkVVZlxgGIGH/IR6uWo/152JYuWoVIfWudzqSMVnapbrEuOwRg4gsd/97SkROpnicEpGTmR3WmCsRWKYUHfs/S/s8JYj7cI6dUjLmKlz2ziBVbeL+t6D34xhz5doP7U/zgiU59sbHRFUqQ+iLjyJ+WW0sKmOyPk/uY3gzI/OMcVLIcw+zu2YYNw56ip+7Pe90HGOyJU++TrVMZ551k2GyFPHzo9bEkZQMDObMN78Q+90vTkcyJtvJSBtDLxH5E6jq7jzv/GMP8Kf3IxrjmWsqV+TnJUuoEVCAw48PJeHAYacjGZOtZOSIYTquvpHm879+ku4G6qnqg17MZswVy9e4DoX6deXdmN280uB2NCHR6UjGZBuXLQyqGqOqe1W1M3ASKAmUA64XkabeDmjMlSrW/zF2hwSx89gRTv+8yuk4xmQbnjQ+9wCW4RqXeZj736HeiWXM1fPz82P25rW83eAOop5+nYTIf5yOZEy24Enj8zNAA2Cfqt4K3ACc8EYoYzJLcEghrpk8gkOxJ+lZvSGxG7Y5HcmYLM+TwhCvqvEAIpJHVbcBVb0Ty5jME1SpLDvuuZHZsQdYcncPkk/HOR3JmCzNk8IQKSIhwFfAjyLyNWAD9Zhsodt7o/j9gaeoESccfKCf3RltzCV40olee1U9oapDgVdw9ax6j7eCGZPZak19i4DSJfn216WsfmGk03GMybI8aXz+SUTaAKjqL6o6H5jgtWTGZDIJDKDArDcZdHoXb3wwgaSTsU5HMiZL8uRUUgVggIi8mmJeuj3zGZNVFataie/f/4SR+SoR9fwoO6VkTDo8KQwngBZASRH5RkQKeyeSMd51Y8+HKfVKL459tZi5zTugSUmXf5ExuYhHQ3uqaqKq9gbmAsuBEt6JZYx3hfTtwth8MTy47Es2dB/kdBxjshRPCsMH55+o6mSgG/BDJucxxidEhJHLFvBeiToU/G4FJ2ctdDqSMVmGJ1clfZhmeq2qPpb5kYzxjdDyZXnw+6kgwpo+rxK3Y4/TkYzJEmwEN5OrBdetQWzf+7k7OoJXb2vndBxjsoSMdKJ3YQQ3VS2U4lFQVQt5P6Ix3hX+ytP0b3Ab7ePyETt/idNxjHGcT0dwE5HWIrJdRHaJyMBLrHefiKiI2OWwxideXb6AMg1v4J+nRrJ/8lyn4xjjKJ+N4CYi/rhuiLsTqAF0FpEa6axXEFeHfdZPsvEZCQrkmk+GM/z4dlr2eISjPyx3OpIxjvHlCG4NgV2qultVzwEzSb9LjdeAN4F4D7ZtzFULuLYE9w14lrvyhHLisSEkHI52OpIxjvDlCG6lgf0ppiPd8y4QkbpAGVX97lIbEpGeIhIhIhFRUVEeRDDm0u4Z3I+Xnu+HnE3g7zZPoIk28pvJfTI8ghvwIHAL0FVV9wEFRKRhZgURET9gDNAvA5k+UtX6qlo/NDQ0syIYA0Cx4U+xv+q13L7+O75u38PpOMb4nCdtDBOAxkBn9/QpPOtE7wBQJsV0mHveeQWB64GlIrIXaATMtwZo42siQvissVwbXADW7SBhv438ZnIXTwrDjaraB/e5f1U9DgR58Po1QGURqSAiQUAnXKencG8vRlWLq2p5VS0PrATaqWqEB+9hTKYoUroUSzauJTxPIf7p9jLnDtopS5N7eFIYEtxXFimAiIQCyRl9saomAk/hGit6KzBbVTeLyHARsTuLTJYTVKksJT54hfErf6Z71bokHotxOpIxPuFJYXgHmIerd9WRuDrR+68nb6aqC1S1iqpWUtWR7nlD3GM7pF23uR0tGKflb3UzZ6qX49i5ePa16WU9sZpcwZO+kqYB/XEVg4PAvar6hbeCGZNVvP37j7x/Yxv4az9H+tjIbybn8+TO5zxAXaAwUAzoICJDvBXMmKwiIDCQ0t9O4Fj+ALp/+i7b3/7E6UjGeJUnp5K+xnVDWiJwOsXDmBzPv1AB5K3nWZEQw8r3p5Acd9bpSMZ4TYAH64apamuvJTEmi6t/391smfMNp58YQdQLowl99yX8/Dz5bmVM9uDJb/UKEanltSTGZAPX/OcOivR/jG+nTmdMlZvszmiTI3lSGJoAa929o24UkT9FZKO3ghmTVYU8/wgzAk8xdd9mDjz2itNxjMl0npxKynBPqsbkZH7+/nyxYRUHm3fj7PfLOT5uKkWeecjpWMZkGk8uV92X3sOb4YzJqoqULkXl7z/mnL8wYtBgji781elIxmSaDB8xiMjz6cyOAdaq6vpMS2RMNhF0XRn+fr4DYwc8TcVnXqb3liWIv7/TsYy5ap60MdQHnsTVVXZp4AmgNTBRRPp7IZsxWd4d/fvy+5C3aH1MOfbfiU7HMSZTeFIYwoC6qtpPVfsB9YASQFOgmxeyGZMtNBj6HIW63kPEmIlMaXwPmpzhLsSMyZI8KQwlgJR39SQAJVU1Ls18Y3IVEaH468/y3+SDDF79Awefed3pSMZcFU+uSpoGrBKRr93TdwPTRSQ/sCXTkxmTjUhgANNWLGFP68eJn7mQE+HVCOlxn9OxjLkinlyV9BrQEzjhfjypqsNV9bSHQ3wakyOF1axGw4Wfon7C5OdfJnrxCqcjGXNFPL2f/y/gd+APIJ+INM38SMZkX3lqXkfUwId5/tR2xnazbrpN9uRJ76o9gGW4BtoZ5v53qHdiGZN9NX7ucb7uPZCH4wsSPex9p+MY4zFPjhieARoA+1T1VuAGXKeUjDFp3D3hdYo+fj9/T5jG1Mb3WJ9KJlvxpDDEq2o8uMZmUNVtQFXvxDIm+ys2oi+j80Tz5Krv2NLpeVTV6UjGZIgnhSFSREKAr4Af3VcnWZcYxlyE+Pvz7qolfFa2EcG//EH00PecjmRMhnhyVVJ7VT2hqkOBV4CPcQ3cY4y5iGJlStNu+RwIDmLh2A/YNf4zpyMZc1mS0cNbEakPvAyUw3X/gwCqqrW9F+/y6tevrxEREU5GMOayDi35nSotbqF1vhLMPLgD/0IFnI5kcjkRWauq9dNb5smppGnAp8B9uG5ua+v+1xhzGaVubcwXI0fzct4KHO4+BE2wxmiTdXlSGKJUdb6q7rFut43xXOuXnqXi2EGcWrKa92o049w/R52OZEy6PCkMr4rIxyLSWUT+c/7htWTG5ECFHryL7+tdy1O7VjCv8b0kx55xOpIx/+JJX0mPAtWAQOB895EKfJnZoYzJyfoumEGJW5Oov/UIkXf0pMzSyUigJ/8VjfEuT34bG6iq3bdgzFXy9/en88+ziLy9Bzs2bmJls/vpsPxLxM/THmqM8Q5PfhNXiEgNryUxJheRgABKL/iAoQn76bdqIcfmLHI6kjEXeFIYGgHrRWS7iGwUkT9FZKO3ghmT0/nlC2bG6l+ZUqslMQPGcnbDdqcjGQN4VhhaA5WBVtjlqsZkirDqVbjl+0n4FS7A2Fvb8+crbzsdyZiMtzHYpanGeEdAqVD8JgxgxC03sW30aMZWKE/hx9o7HcvkYtbaZUwWUPbmhvwyfQ4DC1Xi6IAxnJq72OlIJhfzaWEQkdbuNopdIjIwneXPi8gWdxvGTyJSzpf5jHFS3U73EjZ9FKc0iWce7saR+T85HcnkUj4rDCLiD0wA7gRqAJ3TucrpD6C+u/+lOcD/+SqfMVlB/hY3sqdPO6bFHWRx/xF2A5xxhC+PGBoCu1R1t6qeA2aSpndWVV2iquf/J6wEwnyYz5gsof2wgWycNJMbTyiHOvcnMeqY05FMLuPLwlAa2J9iOtI972K6A9+nt0BEeopIhIhEREVFZWJEY7KGqo8+QMn3X+GXX5fRoVwNTq2yK8ON72TJxmcReQioD4xKb7mqfqSq9VW1fmhoqG/DGeMjBdq34ECrevx57iTb7ulF/B9bnY5kcglfFoYDQJkU02HueamIyO24xn1op6pnfZTNmCzpxdmf8st/36Vokh+RbXoRu+ZPpyOZXMCXhWENUFlEKohIENAJmJ9yBRG5AfgQV1E44sNsxmRZpft3J+SlHrx2YgedmrYgbvf+y7/I5GiqyokJM0g6FuOV7fusMKhqIvAUsAjYCsxW1c0iMlxE2rlXGwUUAL4QkfUiMv8imzMmVyn2fFcq3d6UkgQS3X8MyXF2MJ1bJScmEnl7D6KHvsep2Qu98h4ZHtozq7KhPU1ucnLad0Q99yYnapSl0pCnCLmtkdORjA8lnzvH/lu6krg7kuCb6nDtV+8gIle0rUsN7WmdwBuTjRR68C7OSjKtHulItTY/8MWCb8h3241OxzI+kBwXz/6bHyZx/z8E31SHUvPGXXFRuJwseVWSMebiQrvczfPdetA1uBSHOr3IqXnWfUZOF3voCD3DahK9L5K8tzag9Nfv4ufF8TusMBiTDfX5+B1aT3RdzT2969P8Ne5ThxMZb0mOPcPPD/Rm6rG9rAsvzbWzx3j9Pa0wGJNNFerQioD3XqJf7A4GDxrM2V3WAXJOE7dtNwfaP03NnUdZ/+6n9Ph5jk/e1wqDMdlYufvvZMHHnzHo2locbNuHMyv+ILtfUGJcts9dQJ1atVkSsYprJg2n2lOP+Oy9rTAYk801fexBqi78iMR8eeh02x183aIzmpjodCxzFU59uZgzT44gRPwp8WRn8rdp6tP3t8JgTA4QVKkswZOHs1nPsHnNWiJb9LCeWbOpn54fxsGeQynkH8ji+d/R+vVBPs9ghcGYHKJM7Zr8eeBvHq59I+e2/MXKOvdwbv8hp2MZD0S88R6t3x7GxwmHKP39hxRs3cSRHFYYjMlB8hUvSpmlkzneuAZ37V7KS3Vvtbuks4Hkc+c4OuJDirw1g9evb8bAiMUE31DNsTxWGIzJYSQwgBvmvceTTe/grqQCHGz/NAmHrHv6rOr4zj10LHEdG8dMpOBDbem37keKVa7gaCYrDMbkQP7+/ry+9Fsafj6Gs1t283L1xix7+Dk0OdnpaCaFM79EsKFJF5ac/IfNtcIIHdMfCXS+QworDMbkYAXuakqeKa8xOXY/0+Z+QeTtPUiKOeV0LAOsGvAGhzo8RxkNZN2oD+j7y1de6+LCU1YYjMnhwpo3Yt2GDfSveTPn/txJRK12Nq6Dg1SVb57oz03/9xLfJ52g1LxxlO33mNOxUrHCYEwuEFazGhV/n0bAPc14+OBqOjVtQYKN6+BziQeOcLj7EKrMXcGL1Rrx8NoF5Lv5Bqdj/YsVBmNyCQkIoNzHI3jmscd5oGAYka2fJHbhcpLPxDsdLVf45ZX/4/aK1Tg8/ydKDuvD61tWULh8mcu/0AFWGIzJZfp8NJZOv39NQFhJxnV4lGfDruf06o1Ox8qxks+e40D7Zzjw7lR2J8VxdlA3Qvp0zjLtCemxwmBMLhRUqQxh33/AzrDC/BEbzcG7ehP18jtoUpLT0XKUXV8t4oOyDYlfvo6bK1Zl666dNHy+p9OxLssKgzG5lOQJYvLWCL787HP88wTx9wfTeatiI+ulNZPEfrOUAQ89yqCozfBAS8qunU3+8mFOx8oQKwzG5HKhndtS9o85fFEkiZf+XsvKzs9y7q+/nY6Vbe1dsIQ/WnTl8GOvMKR+C379aj6VJgzJ0qeO0rLCYIwhILQoI7evZNErb1D2eDyRzR5lfuuHiVu/zelo2Uby2XPs6/YSje9uzaDl3xHSryu1f/yUWu3udDqax6wwGGMA193Stw3vT5nlnxNZtwLtF01l2M13cviJYdZT62Xs+mA6e6vcReJ3yxlUogYjJ31IsYE9ssRdzFfCCoMxJpWAa4pzy9cT+aD3C3QsWo7YLxezokor9r03zQYBSiP5TDxzOvWiRq+HWBVzhIIPt+OpyAjqPnif09GuimT3D7p+/foaERHhdAxjcqTk2DNEDRjDfRPf4lDyWVYOeYti/R/DLziP09Ecde7YCXYMHkOB5ZuIPfgP40LiGDLjU8o0qON0tAwTkbWqWj/dZVYYjDGXs+rLb9k4cjy3/X0a/9IlWHdNXtqMGkJwrSpOR/OppDPxHHt1PJ3eeYMjyWf5tvF/KPHm8+RtHO50NI9ZYTDGZIq45euY12sAD25azNsFq3L/zc0JHdOfPNUrOh3Nq/TsOb7r8QJVFm/EPzGJH5Nj8L/zZrp/Nh4/f3+n412RSxUGa2MwxmRY3iZ16bD2Oyb27sed11TkbMRmZjRqx/vX38rZ3ZFOx8t0yUlJnF70G180uY+7p77LwvgoCvXswBMHI3h86vvZtihcjh0xGGOuiCYlcXLad3Ts8wS74mJYVLUVRXs9QFCD68kbXjVbt0PErdvC+O59ybP/CHdJCH6lS/DNdYXoPvFt8hcJcTpeprBTScYYr0lMTGTXrG8oNHcpp5atpc2JdTxUIIzeHR8k5JmHyFOjktMRMyQ59gw7x0wizxc/kfTPUTqe2EiRoGDmfPgJhTq2zraXnl6MnUoyxnhNQEAA1R5sz7VfjqPQ/HE0LH8dZSQPsV8uZtMtD9K/RA12jJnkdMx0JZ+J59yuvznx3kyeva4eDYY+y5lDRwi6vjJfzpjFolOHKPxQ2xxXFC7HjhiMMZkuOfYMJ6d9x6y3xvHkzt/4OqQOda6vRVSNMhw/cYJGT3YlX+NwR043ndtzgNgvFrF86mxe3bqCcQWrUto/mM3li/BnkUCe+mhslu0OOzNd6ojBp2VQRFoD4wB/4GNVfSPN8jzAFKAeEA10VNW9vsxojLl6fgXyEfJEB554ogMt1/xBSMQO4n5YwUefT2FS3AHW/vQH+f0D2VUiH341KtJ8zKsEli6R6TmS48+S9E80yTGn2Pvzcvq/MpgO/sVpHBRC/qQ4EgP8iG9zE2WG9qPSdWVpl+kJsiefHTGIiD+wA2gJRAJrgM6quiXFOr2B2qr6pIh0AtqrasdLbdeOGIzJPg5u38WvE6dwy75Y4tdtoXdkBFsSY1lS/EaCqldg6sEd6Ok4uta7mYCypThdsghFa1cl5J4WACQdPU5y3FnXxlTRuHg0MYmgKuVJOhLNr29MwP+vA5Q/FkfcoSN0ObiatsGhPBJ8LfGaRKsTf/BitRt5pPtjFOzYmsBrM78YZRdZ5YihIbBLVXe7Q80E7gG2pFjnHmCo+/kcYLyIiGb3813GGACurXodHUcPvzA9Yc069v70GyHJwZzdsJ2FaxagCQl0/L0A/L6B+06sp7AE8mmJG5ACeen592rCNJDBBVz3Tdx/YgM1AvIzvMB1ADwQvYrbgory34KVCQzOwzUFClOkRnVKvvg8eepWZ39YSfz8rGn1cnxZGEoDKQeZjQRuvNg6qpooIjFAMeBoypVEpCfQE6Bs2bLeymuM8bJKDepSqUHdC9O/MJqzJ2PhwBESI/+hz/QZ+B07SeF6jUiOPUO5eQe41j+YPBVqAELzvQmULR5KkXvvIyC0KFMO7qVMzWpUuLUJfgXy8Z1zu5at+fJU0v1Aa1Xt4Z5+GLhRVZ9Ksc4m9zqR7um/3OscTW+bYKeSjDHmSmSVy1UPACmb+sPc89JdR0QCgMK4GqGNMcb4iC8LwxqgsohUEJEgoBMwP80684Gu7uf3Az9b+4IxxviWz9oY3G0GTwGLcF2uOklVN4vIcCBCVecDnwCfi8gu4Biu4mGMMcaHfHofg6ouABakmTckxfN4oIMvMxljjEnNrtsyxhiTihUGY4wxqVhhMMYYk4oVBmOMMalk+95VRSQK2HeFLy9OmruqcwHb59zB9jl3uJp9LqeqoektyPaF4WqISMTF7vzLqWyfcwfb59zBW/tsp5KMMcakYoXBGGNMKrm9MHzkdAAH2D7nDrbPuYNX9jlXtzEYY4z5t9x+xGCMMSYNKwzGGGNSyRWFQURai8h2EdklIgPTWZ5HRGa5l68SkfIOxMxUGdjn50Vki4hsFJGfRKScEzkz0+X2OcV694mIiki2v7QxI/ssIg+4P+vNIjLd1xkzWwZ+t8uKyBIR+cP9+93GiZyZRUQmicgR90Bm6S0XEXnH/fPYKCJ101vPI6qaox+4uvj+C6gIBAEbgBpp1ukNfOB+3gmY5XRuH+zzrUA+9/NeuWGf3esVBJYBK4H6Tuf2wedcGfgDKOKeLuF0bh/s80dAL/fzGsBep3Nf5T43BeoCmy6yvA3wPSBAI2DV1b5nbjhiaAjsUtXdqnoOmAnck2ade4DP3M/nAC1ERHyYMbNddp9VdYmqnnFPrsQ1ol52lpHPGeA14E0g3pfhvCQj+/w4MEFVjwOo6hEfZ8xsGdlnBQq5nxcGDvowX6ZT1WW4xqe5mHuAKeqyEggRkVJX8565oTCUBvanmI50z0t3HVVNBGKAYj5J5x0Z2eeUuuP6xpGdXXaf3YfYZVQ1p4wRn5HPuQpQRUR+E5GVItLaZ+m8IyP7PBR4SEQicY3/0tc30Rzj6f/3y/LpQD0m6xGRh4D6QDOns3iTiPgBY4BuDkfxtQBcp5Oa4zoqXCYitVT1hJOhvKwzMFlV3xKRxrhGhbxeVZOdDpZd5IYjhgNAmRTTYe556a4jIgG4Dj+jfZLOOzKyz4jI7cDLQDtVPeujbN5yuX0uCFwPLBWRvbjOxc7P5g3QGfmcI4H5qpqgqnuAHbgKRXaVkX3uDswGUNXfgWBcnc3lVBn6/+6J3FAY1gCVRaSCiAThalyen2ad+UBX9/P7gZ/V3aqTTV12n0XkBuBDXEUhu593hsvss6rGqGpxVS2vquVxtau0U9UIZ+Jmioz8bn+F62gBESmO69TSbh9mzGwZ2ee/gRYAIlIdV2GI8mlK35oPPOK+OqkREKOqh65mgzn+VJKqJorIU8AiXFc0TFLVzSIyHIhQ1fnAJ7gON3fhauTp5Fziq5fBfR4FFAC+cLez/62q7RwLfZUyuM85Sgb3eRHQSkS2AEnAi6qabY+GM7jP/YCJIvIcrobobtn5i56IzMBV3Iu7201eBQIBVPUDXO0obYBdwBng0at+z2z88zLGGOMFueFUkjHGGA9YYTDGGJOKFQZjjDGpWGEwxhiTihUGY4wxqVhhMMYYk4oVBmOMMalYYTDmEkQkRER6p5he4aX3CRORjt7YtjGessJgzKWF4BqvAwBVvclL79MCV5/7xjjOCoMxl/YGUElE1ovIKBGJBRCR8iKyTUQmi8gOEZkmIre7u7feKSINz29ARB4SkdXubXwoIv4p30BEmuDq+fV+9zoVfbqHxqRhhcGYSxsI/KWqdVT1xTTLrgPeAqq5H12AJsALwCC40IlbR+BmVa2Dq7+iB1NuRFWX4+oc7h73+2TnTu5MDpDjO9Ezxov2qOqfACKyGfhJVVVE/gTKu9dpAdQD1rg7K8wLpNebbVVgm9cTG5MBVhiMuXIpx7BITjGdzP/+bwnwmaq+dLGNuLvDjnGPHmiM4+xUkjGXdgrXID9X6idcbQclAESkqIiUS7NOebL5uMQmZ7HCYMwluMcu+E1ENonIqCt4/RZgMPCDiGwEfgTSDtS+DVdf+5tExFtXPRmTYTYegzHGmFTsiMEYY0wqVhiMMcakYoXBGGNMKlYYjDHGpGKFwRhjTCpWGIwxxqRihcEYY0wq/w+iMqiVqCp7ZAAAAABJRU5ErkJggg==\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": 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
}