qc-mentorship-program/Task2/task.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from scipy import linalg as la\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from qiskit import Aer, execute, QuantumCircuit, QuantumRegister, ClassicalRegister\n",
    "from qiskit.circuit import Parameter\n",
    "from qiskit.aqua.operators import MatrixOp\n",
    "\n",
    "from qiskit.aqua.algorithms import VQE\n",
    "from qiskit.aqua.components.optimizers import SPSA\n",
    "from qiskit.aqua import QuantumInstance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Matrix and compute its eigenvalues and eigenvectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues:  [ 1. -1.  1.  1.]\n",
      "Eigenvectors:  [[ 0.          0.          1.          0.        ]\n",
      " [-0.70710678  0.70710678  0.          0.        ]\n",
      " [ 0.70710678  0.70710678  0.          0.        ]\n",
      " [ 0.          0.          0.          1.        ]]\n"
     ]
    }
   ],
   "source": [
    "M = np.array([[1,  0,  0, 0], \n",
    "              [0,  0, -1, 0], \n",
    "              [0, -1,  0, 0], \n",
    "              [0,  0,  0, 1]])\n",
    "eigvals, eigvecs = la.eig(M)\n",
    "# The eigenvalues of A are\n",
    "print(\"Eigenvalues: \", eigvals.real)\n",
    "# The eigenvectors of A are\n",
    "print(\"Eigenvectors: \", eigvecs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NOTE: Our lowest eigenvalue is -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Process\n",
    "### 1. Decompose M into sum of Pauli terms\n",
    "\n",
    "Decomposition involves only terms consisting of the same matrices, i.e.: II, XX, YY, ZZ, and the coefficients are from the set [-1/2, -1, 0, 1, 1/2].\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SummedOp(\n",
      "[0.5 * II,\n",
      "-0.5 * XX,\n",
      "-0.5 * YY,\n",
      "0.5 * ZZ])\n"
     ]
    }
   ],
   "source": [
    "H = MatrixOp(M).to_pauli_op()\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Choose ansatz\n",
    "\n",
    "The ansatz you can use is: (RX I) CX (HI) |00>, where angle in RX is your variational parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ansatz():\n",
    "    \"\"\"\n",
    "    Creates ansatz and returns a quantum circuit.\n",
    "    :return: quantum circuit\n",
    "    \"\"\"\n",
    "    phi = Parameter('φ')\n",
    "    qc = QuantumCircuit(2)\n",
    "    qc.h([0])\n",
    "    qc.cx([0], [1])\n",
    "    qc.rx(phi, [0])\n",
    "    \n",
    "    return qc, phi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR0AAAB7CAYAAABJsFSIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAALpElEQVR4nO3df0zU9x3H8ecdCtT6owUUBISKgKtXDwWdlrYyN6u49YerorGdS62bVG02tUv6w9CsoeMPZzuzZa2zbWa2WtJAaNfWH8uaIWqwK4RWy7pK5YeAIiqKeoqocPsDpZwCYr37fEFej+QSeN/x+b7vcrzu8/18v3dnc7vdbkREDLFb3YCI9C8KHRExSqEjIkYpdETEKIWOiBil0BERoxQ6ImKUQkdEjFLoiIhRCh0RMUqhIyJGKXRExCiFjogYpdAREaMUOiJilEJHRIxS6IiIUQodETFKoSMiRil0RMQohY6IGKXQERGjFDoiYpRCR0SMUuiIiFEKHRExSqEjIkYNsLoBMWPlZmu2u/4Ja7ZrpU/WWbPdGb+xZrs3SjMdETFKoSMiRil0RMQohY6IGKWFZBGLLMgMx3XuJHa7H3a7H+Oi7mVV2puMuGOU1a35lGY6IhY4fuoQJ07X8frKEj76nYu/P19B49mjvLXlOatb8zmFjogF9tcUMShwKFEjvgfA0NuDiY1IpNF11OLOfE+hI2KB/TVFjI2cjM1mo6XlEkVfb6dg73v8cOLjVrfmc1rTkVvOpRZougD+AyGglz7D99cUUVq1mzkZd9B0wUXgwEE8m/Y2KQnzrW7N5zTT6aC1tZV169YRFxdHYGAgCQkJFBQUMHbsWJYuXWp1e0a99Uw4pflvedTcbjdv/GIoB4ret6ir7p1wQc5n8GIOZOTBc+/Bxnw4UG91Z9f6praYFx/P5oPMRt5dU0PIsEiqj35tdVtGKHQ6WLJkCZmZmaSnp7Nt2zbmz5/PwoULqaioICkpyer2jHGdOMTZxjqGR0/wqJ86WsGF82cIjZlkUWddq2uEddug8Bu40PJt/X918OdP4LMK63q72qHjBzjTdJK4yEQAgoeOZF7Ks3z86QZaW1uBtoD/1Z/u5fDxco6fOsyq1x/gyIkqC7v2nl46+TQvOzubTZs2sWPHDlJSUgCYPn06JSUl5OXlkZiYaHGH5tRXFGGz+xEceY9H/Xj1XgYNC2VIcO86pOt2w193te1SuTu5DiD7UxgzAoIHG2/vGvtrihgyKIjQO6Pba8mOR1mfu5TSqt04Y6Zhs9l4anYWGz5azbnm06ycu5GwoLusa9qLNNO5LCsri9TU1PbAuSI2NpaBAwfidDoBqKqqIiUlhfj4eMaPH8+uXbusaNen6iuKuDMsngH+gR71Y9V7GTG6981yDtTD0dPXBk5HbnfbLKg3KKspIjZ8okdt2O0hjLsrmV37cttrMSOdlFbuYtbkp4gOvdt0mz6jmQ5QW1tLaWkpq1atuua66upqHA4HAQEBAKSnp7NgwQKWL19OYWEhaWlpVFZW4u/v75VebDabV8a52q/f6e5f0lN9RRGN9Qf4y9MhHvWLzS4mPfzCDW3XV/eno6lzX+b7czK63Zbb7WbzR0U8kjjF5/386/fdP9ZPP/Jap/U/LP/2Bayp2cUr7yxg+aN/ZHvR2zyYtOi62zXxWHfH7e7Zc0yhQ1voAISFhXnUm5qaKCgoYPbs2QAcP36c3bt38+GHHwKQnJxMeHg4+fn5zJo1y2zTPlRfWcyUx37L3ff/3KO++YXxhPbCmY7N7nf929hs2P36xtP9wqVmMt+Zz+M/WsOE2OnsLc9n95fvc//4n1rdmldo9woICWl7RS8rK/Oor127lrq6uvZF5OrqakJDQ9tnPQCjR4/m4MGDXuvF7Xb75NJTjUcO0Hz2JNHOWQwJjmy/tFw8T/O5RkJjJveK+9PxsuHVNdd9lbcBj85INNLPzfIfEEDWkq1MiJ0OwLPz3+5R4Ji4b964330j+n0sJiYGp9NJVlYWQUFBREREkJuby9atWwH61ZGr+ooiBgQMYnhUgke97ptCBgePYtCwERZ11jVHBAwNhDPnu17XcQP3xZnsSrqimQ5gt9vJycnB4XCwbNkyFi9eTEhICCtWrMDPz699ETkqKor6+nqam5vb/7ayspLo6Oiuhu5z6iuKCB09+ZpdkboDe3rlrhWAnx1+dh/Y7W0zmo6u/D7zHogMMt2ZdMbm9sZ88Ba1aNEi9u7dy759+9prM2fOZM6cOe0LyfPmzaOqqsprC8m+0h8+rrTyGGzZ63kyYPBgeNABU8aAqXVWfVxp97R71Y3i4mKmTp3qUduwYQNPPvkk69evx9/fn+zs7F4fOP3F6OHwzAxocEHmP9pqax4Bu7UHdeQqCp0uuFwuysrKWL58uUc9JiaGnTt3WtSV9ETHEwAVOL2PQqcLgwcPpqWl5fo3FPGC02cbePlvcxl3VzLF+7fz/MLNt9QJgR1pIVmkF9jyn43MmryYrw9+ysq5G8nd+arVLfmMQkekFzh6spr4yEm0uluJDh1Hw6lDVrfkMwodkV4gPCSW8sNfAFB1pJSRwWMs7sh3FDoivcCPp/yS/C+yKT/8BW98uJK0lD5y/Ps7UOiI9AK3Bw7llac+Zkz4BF5btvOW+RiLzih0RMQoHTIX6UVeXbbD6hZ8TqHTT5h8O0J/11fejmAV7V6JiFEKHRExSqEjIkYpdETEKIWOiBil0BERoxQ6ImKUQkdEjFLoiIhRCh0RMUqhIyJGKXRExCiFjogYpdAREaMUOiJilEJHRIxS6IiIUQodETHK5na73VY3IXKzzl+EL2vgYAPUNUL50bb6pNEwKgjuifT8jnOxjkJH+rSzzbD9S/isHJovdX07GzAuAn6SAOF3GmtPOqHQkT7rq0OQ/SmcOd/zv/Gzwczx8OA9YLf5rjfpmkJH+qTPKiB7D3zXJ+/kGFg4VcFjBS0kS59TduTmAgegqAK2fOG1luQGKHSkTzl/Ed7tQeCsf+L63/X176+g8pjXWpMeUuh00Nrayrp164iLiyMwMJCEhAQKCgoYO3YsS5cutbo9oS0oGs95Zyw3kFfsnbGk5/QNnx0sWbKEvLw8MjIySEpKorCwkIULF3Ls2DFWr15tdXv93qUWKDzg3TFrTkB1A0QFe3dc6ZpC57Ls7Gw2bdrEjh07SElJAWD69OmUlJSQl5dHYmKixR3KN/XguoEjVT1VXKnQMUm7V5dlZWWRmpraHjhXxMbGMnDgQJxOJwAvvfQS8fHx2O12cnNzrWi136pu8M24NSd8M650TqED1NbWUlpaSlpa2jXXVVdX43A4CAgIACA1NZXt27czbdo00232e0dO+WbcukbfjCud0+4VbaEDEBYW5lFvamqioKCA2bNnt9eSk5N92ovNphNHuvLQqg8Yk/SoR+16R6i6un7l5m9/dp1rxmYLvMnupKen/GmmA4SEhABQVlbmUV+7di11dXUkJSVZ0ZZcpeWiDxZ0fDiudE4zHSAmJgan00lWVhZBQUFERESQm5vL1q1bAYyGjk4Q79q2ffDPLz1rHWcsHV2Z4XR1fUfxUcP0uBukmQ5gt9vJycnB4XCwbNkyFi9eTEhICCtWrMDPz699EVmsFRXUt8aVzmmmc1l8fDz5+fketUWLFjFu3Dhuu+02i7qSjuLCYJA/nLvg3XEnRnt3POmeZjrdKC4uvmbXKiMjg8jISPbs2UN6ejqRkZGUl5db1GH/4j8Apozx7pgj74DRw707pnRPodMFl8tFWVnZNScFZmZmUltbS3NzMw0NDdTW1jJmjJf/E6RLMxwwxIsHmh5LAh0wNEsfbSF9zn9r4c2Cmx/ngbEwd9LNjyM3RjMd6XMckTcfFs5RMEfvbLGEZjrSZ5VUQU4RNN3AwrIN+MHd8NAE8NNLriUUOtKnnWqCjz+HkoPQ0tr9bWOGw8MTtXBsNYWO3BJc5+Hzg23fBnG4EZovts1kRgxt+zYI56i2I1ViPYWOiBilvVoRMUqhIyJGKXRExCiFjogYpdAREaMUOiJilEJHRIxS6IiIUQodETFKoSMiRil0RMQohY6IGKXQERGjFDoiYpRCR0SMUuiIiFEKHRExSqEjIkYpdETEKIWOiBil0BERoxQ6ImLU/wHsHmQ2pZdQ9wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 358.792x144.48 with 1 Axes>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qc, param = ansatz()\n",
    "qc.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. VQE implementation\n",
    "#### 3.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "backend = Aer.get_backend('qasm_simulator')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def measure(circuit):\n",
    "    \"\"\"\n",
    "    Takes the quantum circuit and performs measurements and returns the expectation value\n",
    "    :param circuit: Quantum circuit\n",
    "    :return: expectation value\n",
    "    \"\"\"\n",
    "    result = execute(circuit, backend=backend, shots=1024).result()\n",
    "    counts = result.get_counts()\n",
    "    \n",
    "    expectation = 0.0\n",
    "    for key in ['00', '11']:\n",
    "        expectation += counts.get(key, 0.0)\n",
    "    for key in ['01', '10']:\n",
    "        expectation -= counts.get(key, 0.0)\n",
    "    return expectation / 1024"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# XX expectation\n",
    "def XX(phi):\n",
    "    \"\"\"\n",
    "    Gives the expectation value of XX sub hamiltonian from measurement\n",
    "    on parametric state\n",
    "    :param phi: angle in radians\n",
    "    :return: expectation value of XX \n",
    "    \"\"\"\n",
    "    circ = qc.bind_parameters({param: phi[0]})\n",
    "    # circ.h(qr[0])\n",
    "    # circ.h(qr[1])\n",
    "    circ.ry(-np.pi/2, 0)\n",
    "    circ.ry(-np.pi/2, 1)\n",
    "    circ.measure_all()\n",
    "    \n",
    "    expectation = measure(circ)\n",
    "    return expectation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# YY expectation\n",
    "def YY(phi):\n",
    "    \"\"\"\n",
    "    Gives the expectation value of YY sub hamiltonian from measurement\n",
    "    on parametric state\n",
    "    :param phi: angle in radians\n",
    "    :return: expectation value of YY \n",
    "    \"\"\"\n",
    "    circ = qc.bind_parameters({param: phi[0]})\n",
    "#     circ.sdg(qr[0])\n",
    "#     circ.h(qr[0])\n",
    "#     circ.sdg(qr[1])\n",
    "#     circ.h(qr[1])\n",
    "    circ.rx(np.pi/2, 0)\n",
    "    circ.rx(np.pi/2, 1)\n",
    "    circ.measure_all()\n",
    "    \n",
    "    expectation = measure(circ)\n",
    "    return expectation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ZZ expectation\n",
    "def ZZ(phi):\n",
    "    \"\"\"\n",
    "    Gives the expectation value of ZZ sub hamiltonian from measurement\n",
    "    on parametric state\n",
    "    :param phi: angle in radians\n",
    "    :return: expectation value of ZZ \n",
    "    \"\"\"\n",
    "    circ = qc.bind_parameters({param: phi[0]})\n",
    "    circ.measure_all()\n",
    "    \n",
    "    expectation = measure(circ)\n",
    "    return expectation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost(phi):\n",
    "    return 0.5 + (-0.5 * XX(phi)) + (- 0.5 * YY(phi)) + (0.5 * ZZ(phi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lowest eigenvalue:  -0.5009765625\n"
     ]
    }
   ],
   "source": [
    "phi = np.random.uniform(0, 2 * np.pi)\n",
    "\n",
    "print('Lowest eigenvalue: ', cost(phi=[phi]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenvalues = []\n",
    "angles = np.linspace(0,2*np.pi,250)\n",
    "for phi in angles:\n",
    "    eigenvalues.append(cost(phi=[phi]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnAAAAGDCAYAAACr/S2JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOzdd3hUZd7G8e8vhRoivfdOkKYBLChFUFGaSFGwICKKWFdd264dy+q6K7q6ILqigoJgQaqgFEUJht4EAoK00EE6JHnePzLwRqQEyOTJTO7Pdc2VOWUyd8aSO8855znmnENEREREQkeE7wAiIiIicmZU4ERERERCjAqciIiISIhRgRMREREJMSpwIiIiIiFGBU5EREQkxKjAiUjQmVlPM/vGd45TMbNpZtbHd47sZma9zOwH3zlE5MyowIlIljCzNWZ2wMz2Zni8BeCcG+acu9J3RhGRcBHlO4CIhJX2zrkpvkOIiIQ7jcCJSNAdf5jOzK40s+VmttvM3jaz6RkPX5pZbzNbZmY7zWySmVXKsM2Z2V1mttLMdpnZfyxd3sDy+Rn2LREYFSxpZkXMbKyZbQ1837FmVv4keZ8xs48zLFcOvG9UYPk8M3vPzDaZ2QYze8HMIk/wfcoG3r9ohnWNzGybmUWbWfXAz747sG7EKT7Dz8wsObDvDDOrm2HbB4HPYZyZ7TGzBDOrltnP+7j3qW1mk81sR+A13TJsu8bMlgbeY4OZPXyyvCISXCpwIpKtzKw4MAp4HCgGLAcuybC9I/AE0BkoAXwPfHLct2kHNAbqA92Aq5xzh4DPgRsz7NcNmO6c20L6/+/+B1QCKgIHgLfO8sf4AEgBqgONgCuBPxUi59xG4Cfg+gyrewCjnHNHgOeBb4AiQHngzVO85wSgBlASmAsMO277DcCzge+VBAyA03/eGZlZQWAyMDzwPjcAb5tZXGCX94A7nXOFgPOB706RV0SCSAVORLLSl4FRsKOPO06wzzXAEufc5865FGAgkJxh+13AS865ZYHtLwINM47CAS8753Y5534DpgINA+uHk146juoRWIdzbrtzbrRzbr9zbg/pBaf5mf6AZlYq8DM84JzbFyiH/zrufTMaTqBUmpkF9hse2HaE9EJZ1jl30Dl30osJnHPvO+f2BIrqM0ADMzsvwy5fOOdmBz6zYfz/Z3K6zzujdsAa59z/nHMpzrl5wGiga4a8cWYW65zb6Zybe7K8IhJcKnAikpU6OecKZ3i8e4J9ygLrji445xywPsP2SsAbR0sgsAMwoFyGfTIWkP1ATOD5VKCAmTU1s8qkl5gvAMysgJkNMrO1ZvY7MAMofKJDn6dRCYgGNmXIOIj0EasTGQ1cbGZlgMuBNNJHFQH+GvjZZpvZEjPrfaJvYGaRZvayma0KZF8T2FQ8w24n+0xO93kf/7M1zVjCgZ5A6cD260kvhGsDh2EvPsn3EZEg00UMIpLdNpF+uBA4NiqV8Vy0dcAA59zxhwhPyzmXamYjSR/x2gyMDYy2ATwE1AKaOueSzawhMI/0AnW8fUCBDMulMzxfBxwCigdGtE6XaaelT6HSHagDfBooUTjnkoE7AMysGTDFzGY455KO+zY9gI5Aa9LL23nAzpNkP97pPu+M1pF+yLnNSX6Wn4GOZhYN3AOMBCpkIoOIZDGNwIlIdhsH1DOzToGLAvrzx4L0X+DxoyfpBy4Y6HqC73Myw0kvSz35/0OVAIVIP+9tV+CigqdP8T3mA5ebWcXAYcrHj25wzm0i/by1f5pZrJlFmFk1MzvV4djhwC1Al4yZzKxrhgspdgKO9BG64xUivTRuJ71YvniK9zre6T7vjMYCNc3s5sBFFtFm1tjM6phZHkufz++8wPl7v58kq4hkAxU4EclKX9sf54H74vgdnHPbSD+n6h+kF5I4IJH0goJz7gvgFeDTwOHCxUDbzAZwziWQPoJWlvQT/4/6N5Af2AbMAiae4ntMBkYAC4E5pBebjG4B8gBLSS9eo4Ayp4g1hvQLEJKdcwsyrG8MJJjZ3sA+9zvnVp/g9R8Ca4ENgfecdYr3Ov5nOeXnfdy+e0i/IOMGYCPph2VfAfIGdrkZWBP453IX6SVZRDywwEi+iIgXZhZB+jlZPZ1zU33nCXf6vEXCg0bgRCTbmdlVZlbYzPKSPmWIcQajSnJm9HmLhB8VOBHx4WJgFemHM9uTfvXqAb+Rwpo+b5Ewo0OoIiIiIiFGI3AiIiIiIUYFTkRERCTE5KqJfIsXL+4qV67sO4aIiIjIac2ZM2ebc67EibblqgJXuXJlEhMTfccQEREROS0zW3uybTqEKiIiIhJiVOBEREREQowKnIiIiEiIUYETERERCTEqcCIiIiIhRgVOREREJMSowImIiIiEGBU4ERERkRCjAiciIiISYrwWODN738y2mNnik2w3MxtoZklmttDMLsiw7VYzWxl43Jp9qUVERET88j0C9wFw9Sm2twVqBB59gXcAzKwo8DTQFGgCPG1mRYKaVERERCSH8HovVOfcDDOrfIpdOgIfOuccMMvMCptZGaAFMNk5twPAzCaTXgQ/CW7iU1uyZQmb9m6iaP6i5I3My86DO4mwCOLLxpMnMo/PaCIiInIWUtJSWL1zNdv3b2f3od2kVxK4uMLFFM5X2FuunH4z+3LAugzL6wPrTrb+T8ysL+mjd1SsWDE4KQMGzRnEm7Pf/NP6QnkK0aRcE/JG5SUmTwwtK7ekbfW2VCpcKah5RERE5MzsOriLKaunkLA+gYQNCczZNIf9R/b/ab+EPgk0KdfEQ8J0Ob3AnTPn3GBgMEB8fLwL5ns9cskjdI3ryo4DOziYcpCi+Yuy5/AeJiVNYl7yPNwhx+Itixm5ZCQAtYvX5prq19D9/O40LtsYMwtmPBEREQGSdiQxfNFw5myaA0BURBRF8xVl+4HtjF85nkOph8gbmZdGZRrRp1EfLihzAaViSnFe3vOIsPSzz+JKxPn8EXJ8gdsAVMiwXD6wbgPph1Ezrp+WbalOosJ5FahwXoU/re9cp/Ox5845VmxfwYSkCUxImsBbP7/F67Nep3rR6rSr0Y52NdvRqkorlTkREZEs9OvOXxk8ZzATkiawYPMCDCOuRBzRkdEcTj3MzgM7iYqI4q74u+hWt1uOP/3Jjh7L9RYg/Ry4sc6580+w7VrgHuAa0i9YGOicaxK4iGEOcPSq1LnAhUfPiTuZ+Ph4l5iYmIXpz92ug7v4fNnnjFwykmlrpnEo9RDd63ZnSIchxOSJ8R1PREQk5K3asYpm/2vGtv3buLTCpbSr2Y4bzr+B8rHlfUc7JTOb45yLP9E2ryNwZvYJ6SNpxc1sPelXlkYDOOf+C4wnvbwlAfuB2wLbdpjZ88DPgW/13OnKW05VOF9hejfqTe9Gvdl/ZD//nvVv/j717yzYvIDral9HneJ16Fa3G3mj8vqOKiIiEnLW7V5H649acyT1CPPvnE/dknV9R8oS3kfgslNOHIE7kW9Xf8v9E+9n+fblpKSl0KpKK77s/iWF8hbyHU1ERCQkLN6ymH/P+jfDFg0jOiKa7279jviyJxzMyrFONQLnex44OYErql7B4rsXc+DJA/yv4/+YvmY6rT5sxfb9231HExERybGcc/y84We6ftaVeu/UY/ii4dxS/xZ+vuPnkCtvp6MCl4NFRUTRq2EvvrzhSxZtXkTbYW3Zc2gPCesTaPxuY4YvGu47ooiIiHfOOQbPGUyd/9ShyZAmTEqaxN8u+xvrHlzHoPaDqFW8lu+IWU6HUEPE2BVjuW7EddQuXpvl25aT5tIwM8b3GE+bam18xxMREfHiSOoR7p1wL4PmDOKi8hfRu2Fvutbt6nWS3aySYy9ikMxrV7MdH3b6kJ6f9+Tq6lfz9rVv0+GTDlw/8np61utJZEQk/eL7hc3JmSIiIifyzLRn2LZ/G1dVu4p1v6/jvXnvMXfTXB699FFevOLFY/O0hTuNwIWY5L3JlCxYkgiLYP3v67luxHWs3bWWXQd30aRcE76/7XvNISciImFp3qZ5XDD4AiIsgjSXBkDdEnV5vNnj9Kzf03O6rKcRuDBSOqb0seflY8vz8x3pM6n8N/G/9BvXjwlJE7imxjW+4omIiATNKzNfITZvLCvvXcnSrUspkq8I9UvVz5UDF7ljnDEX6N2oN1WLVOVv3/3t2F8lIiIi4SJpRxKfLf2MfvH9KFmwJC0qt6BB6Qa5sryBClzYyBOZh2eaP8O85Hnc+fWdjFo6il0Hd/mOJSIikiVenfkq0RHR3N/0ft9RcgQVuDDSo14PusR14cOFH9L1s66Ue70cd429i8VbFvuOJiIictaWbV3G/+b/j9sa3kaZQmV8x8kRdBFDGDqcepjEjYm8N/c9hi8ezsGUg7Ss3JLXrnyNC8pccPpvICIikkM452j1YSsWJC9g+T3LKVGwhO9I2eZUFzGowIW5bfu3MWTuEP41618UyVeEJXcvITIi0ncsERGRP0lNS6XXV72IiY6hbY221C5em+lrptN3bF8GtRtE3wv7+o6YrVTgAnJjgTvqsyWf0W1UN0Z0GUG3ut18xxEREfmTKaun0OajNkRHRHMk7cix9U3LNeXH23/MNXO8HaVpRITOdTpTu3htXpjxAl3iuuS6/whERCTnG75oOLF5Y1n34DrmJ89n/e/r2XVwF9fWuFa/t46jApdLREZE8kSzJ7jly1sYmDCQyypeRp0SdSgQXcB3NBEREQ6mHGT0stF0rtOZ2LyxXF7pct+RcjTV2Vzkxno3Uqd4HR6c9CDx78YT95845myc4zuWiIgI41eO5/dDv9Pj/B6+o4QEFbhcJCoiill9ZvH9bd/zyfWf4HBc+v6lfLr4U9/RREQkF1n/+3rK/rMsbT5qw/vz3mfJliV8uOBDShUsRcsqLX3HCwk6hJrLxOaNpVnFZgC0rtqajp92pO/XfWlTtQ3FChTznE5ERHKD9+a+R/LeZPJH5+f2MbcfW39fk/uIilA1yQx9SrlY8QLFebf9u5z/9vm8MvMV/tHmH74jiYhImEtNS+X9+e/TplobJvacyILNC1i6dSnrdq/jlga3+I4XMlTgcrm4EnHcVP8m3pr9Fg9c9ABlC5X1HUlERMLYlNVT+G33b7zW5jXMjIalG9KwdEPfsUKOzoETnmnxDEfSjjBgxgDfUUREJMwNmTeE4gWK06FWB99RQpoKnFC1SFV6N+zNkHlD2PD7Bt9xREQkDDnnGLtiLF/98hW31L+FvFF5fUcKaSpwAsBjzR4jNS2VV398FYAhc4fw3PTnPKcSEZFwsGjzIq78+Eraf9KeqkWq8sBFD/iOFPJ0DpwAUKVIFW5ucDOD5gyicuHKPDjpQSIsgn7x/XLVjYNFRCTrHEw5yIMTH2Tw3MGcl/c8Bl49kLvi7yI6Mtp3tJCnETg55vFmj3Mo5RAPTnqQuBJxpLk0xq4Y6zuWiIiEqP/M/g//nfNf+jfuT9J9Sdzb9F6VtyyiAifH1CxWkzsvvJP4svH8cNsPVDyvIl8u/9J3LBERCUEHjhzg1R9fpXXV1gxsO5Ci+Yv6jhRWdAhV/uDta98GwMzoVKsTg+cOZt/hfRTMU9BzMhERCSXvzXuPzfs2M+KyEb6jhCWNwMkfmBlmBkCn2p04mHKQSasm8cH8D3hsymOkuTTPCUVEJKc7nHqYV2a+wmUVL6N55ea+44QljcDJSV1W6TKK5i/K7WNuZ9fBXQCUiSnD/Rfd7zmZiIjkVNv3b6fH5z1Y//t63u/wvu84YcvrCJyZXW1my80sycweO8H2f5nZ/MBjhZntyrAtNcO2MdmbPHeIioiiS50u7Du8j1fbvEq7mu14dMqjLN6y2Hc0ERHJgX7b/RsXDr6QaWumMbjdYNpUa+M7Utgy55yfNzaLBFYAbYD1wM/Ajc65pSfZ/16gkXOud2B5r3Mu5kzeMz4+3iUmJp5b8FzmwJED7D60m9Ixpdmybwv13qlHbN5YXmj5Ap3rdNbVRCIicsyjkx/l9VmvM7P3TJqUa+I7TsgzsznOufgTbfM5AtcESHLOrXbOHQY+BTqeYv8bgU+yJZkckz86P6VjSgNQsmBJRnQZgXOOG0bfQIP/NuBgykHPCUVEJCdITUtl2KJhXF39apW3bOCzwJUD1mVYXh9Y9ydmVgmoAnyXYXU+M0s0s1lm1il4MSWjFpVbsOLeFQxpP4Rl25YxcslI35FERMSTpB1JPD31aQ6mHGT62uls2LOBm+vf7DtWrhAqFzHcAIxyzqVmWFfJObfBzKoC35nZIufcquNfaGZ9gb4AFStWzJ60YS7CIujdqDev/vgqb//8Nrc0uMV3JBERyWZHUo/Q7bNuzEuex6a9mziSdoTYvLG0r9ned7RcwecI3AagQobl8oF1J3IDxx0+dc5tCHxdDUwDGp3ohc65wc65eOdcfIkSuiVUVjEz7m58NwkbEpizcY7vOCIiks1e+uEl5iXPo3XV1rw7910+WvARXep0IX90ft/RcgWfBe5noIaZVTGzPKSXtD9dTWpmtYEiwE8Z1hUxs7yB58WBS4ETXvwgwXNLg1soEF2AdxLf8R1FRESy0YLkBTw/43l61OvBxJ4TaV21NakulZvq3+Q7Wq7h7RCqcy7FzO4BJgGRwPvOuSVm9hyQ6Jw7WuZuAD51f7xctg4wyMzSSC+hL5/s6lUJnsL5CtOzXk8+mP8BeSPzcl/T+6hVvJbvWCIiEkTOOe6beB9F8hVh4NUDiYyIZFTXUcxYO4MWlVv4jpdreJtGxAdNI5L1tu7bysOTH+bTxZ8SaZGsvHcl5WJPeC2KiIiEgS+WfUHnkZ1559p3uCv+Lt9xwlpOnUZEwkCJgiUY2mkoiXckciDlAJ8v+9x3JBERCZLDqYd5ZPIjxJWIo88FfXzHydVU4CRL1CtVj7ol6jJq2SjfUUREJEgGzxnMqp2r+OeV/yQqIlQmsghPKnCSZbrEdeH7td+TvDfZdxQREcliaS6NNxLe4KLyF3F19at9x8n1VOAky3SJ64LD8cWyL3xHERGRLDYxaSJJO5K4r8l9vqMIKnCSheqWqEutYrUYvWy07ygiIpLF3pz9JqVjSnN93PW+owgqcJKFzIwucV2YumYqwxYOIzdd4SwiEs5WbF/BxKSJ9IvvR57IPL7jCCpwksXubXIvjcs25qYvbqLdJ+3Yf2S/70giInIOnHM88e0TREdE0/fCvr7jSIAKnGSpUjGlmNl7Jq9f+TrjV47n5R9e9h1JRETOwX8T/8voZaMZ0GoApWNK+44jASpwkuUiIyJ58OIHufH8G/nHzH/w685ffUcSEZGzMGfjHB6c9CBtq7floUse8h1HMlCBk6D5R5t/EBkRyT0T7mHI3CH8/bu/s+/wPt+xREQkE0YuGUnzD5ofm7A9wlQZchLNwidBUz62PE9e9iRPfvck41eOP7buzvg7PScTEZFTeW76czw97WkuqXAJn3X9jBIFS/iOJMdRnZag+uulf2Vcj3GsvHcldYrXYfji4b4jiYjIKfzzx3/y9LSnubXBrUy9dSplC5X1HUlOQAVOgioqIopralxD9aLV6VGvBzPWzmDd7nW+Y4mIyAl8sugTHp78MF3juvJeh/c0ZUgOpgIn2ebG828E4NPFn+KcY/XO1ZorTkQkB3n1x1dpWLohH3f+mMiISN9x5BRU4CTbVCtajablmjJ0wVBuHH0j1QZW010bRERyiLW71jIveR496/XUyFsIUIGTbNWjXg+WbF3C6GWjickTw/BFOidORCQn+Gr5VwB0rNXRcxLJDF2FKtnqlga3sGzrMno36s2HCz5kyLwh7D28l5g8Mb6jiYjkal/+8iVxJeKoUayG7yiSCRqBk2xVOF9h3mn3Do3LNaZLXBcOphxkwsoJvmOJiORq2/dvZ8baGXSq1cl3FMkkFTjxplnFZpQsWJJRy0b5jiIikquNXTGWVJdKp9oqcKFCBU68iYyI5Lra1zFuxTjd9F5ExJP1v6/nxR9epHxseS4se6HvOJJJKnDiVZe4Luw7so8+Y/qwcc9G33FERHKVFdtXcOn7l7JpzyaGdR6m22WFEP2TEq9aVWnFE82eYPSy0dR8syY/rvvRdyQRkVzBOcdNn9/E/iP7md5rOpdXutx3JDkDKnDiVYRFMOCKASy9eykFogvw2o+v+Y4kIpIrJGxI4OeNP/Nsi2dpVKaR7zhyhlTgJEeoVrQavRr24usVX5O8N9l3HBGRsPfm7DeJzRvLLQ1u8R1FzoIKnOQYtze6nZS0FIbOH+o7iohIWEvem8xnSz7jtoa3aR7OEKUCJzlGreK1uKziZQyZN0T3SBURCZLdB3fzzLRnOJJ2hP6N+/uOI2dJBU5ylD4X9CFpRxIz1s7wHUVEJOw8PfVpSr1WikFzBtGzXk/ddSGEqcBJjtIlrguxeWMZMm+I7ygiImElcWMiz814jrY12pLQJ4GPrvvIdyQ5B14LnJldbWbLzSzJzB47wfZeZrbVzOYHHn0ybLvVzFYGHrdmb3IJlgLRBehZryejlo5i54GdvuOIiIQF5xx/mfQXShQowdBOQ2lSrglm5juWnANvN7M3s0jgP0AbYD3ws5mNcc4tPW7XEc65e457bVHgaSAecMCcwGv1Gz8M9LmgD+8kvsPwRcPp30TnZ4iInK1hC4eRtCOJyIhIvv/te/577X+JzRvrO5ZkAW8FDmgCJDnnVgOY2adAR+D4AnciVwGTnXM7Aq+dDFwNfBKkrJKNLihzAY1KN+Ldue9yd+O79VeiiMhZcM7Rf3x/dh/aDUDdEnW5/YLbPaeSrOLzEGo5YF2G5fWBdce73swWmtkoM6twhq+VENXngj4s2LyA73/73ncUEZGQ9OuuX9l9aDf/vPKfTLt1GpNumkRUhM9xG8lKOf0ihq+Bys65+sBk4IwnCDOzvmaWaGaJW7duzfKAEhw96vWgeIHitBzakps+v4lHJz9Kt8+6MfXXqb6jiYiEhHmb5gFwWcXLaF65OeViNc4RTnxW8Q1AhQzL5QPrjnHObc+wOAT4R4bXtjjutdNO9CbOucHAYID4+HhNLhYiCucrzOJ+i3ntx9d4O/FtUtJSiLRIVmxfwbw75+mwqojIacxLnkekRVKvVD3fUSQIfI7A/QzUMLMqZpYHuAEYk3EHMyuTYbEDsCzwfBJwpZkVMbMiwJWBdRJGSsWU4tUrX2XHX3ew/4n9DGw7kAWbFzB97XTf0UREcry5m+YSVyKOfFH5fEeRIPBW4JxzKcA9pBevZcBI59wSM3vOzDoEdrvPzJaY2QLgPqBX4LU7gOdJL4E/A88dvaBBwk/eqLxERkTSs15Pihcozr9m/ct3JBGRHG9e8jzdpD6MeT2b0Tk3Hhh/3LqnMjx/HHj8JK99H3g/qAElR8kfnZ+7LryLAd8PIGlHEtWLVvcdSUQkR9q0ZxPJe5NpVFoFLlzl9IsYRP7g7sZ3ExURxTs/v+M7iohIjjUvOf0ChgvKXOA5iQSLCpyElDKFytC6amu+Wv6VbngvInKcZVuXkbw3+dgVqA1LN/ScSIJFE8JIyGlXsx39x/dn+fbl1C5e23ccEZEc4cCRAzQd0hQzo1TBUlQrUk13XQhjGoGTkNOuZjsAvl7+teckIiJ+DUwYyLKt6RM0TFszjT2H91C2UFlW7ljJhWUv9JxOgkkFTkJOxfMq0qBUA75eoQInIrlX0o4k7p94P49/m36t37iV48gflZ+5fecyrsc4Xr7iZc8JJZhU4CQktavZjpnrZrLjgGaPEZHcadyKcQCMXTGWTXs2MW7lOK6oegX5o/NzTY1rqFKkiueEEkwqcBKS2tdsT5pL46mpT/HU1Kf4ad1PviOJiGSrcSvHUbJgSVJdKo9OeZQ1u9ZwbY1rfceSbKKLGCQkNS7XmPKx5fnPz/8BYPLqyfx0u0qciOQOew/vZfra6dzb5F4SNyby0cKPALimxjWek0l20QichKQIi2B2n9msvHclT13+FAnrE0jem+w7lohItpiyegqHUw9zbY1r6XNBHwDql6pPxfMqek4m2UUFTkJWmUJlqF60OtfHXY/DMXbFWN+RRESyxbgV44jNG0uzis24vs71lC1Ulhvq3uA7lmQjHUKVkFevZD0qF67MmOVjjv0lKiISrg6nHmZ80niurHYl0ZHRREdGs/q+1URHRvuOJtlII3AS8syMDjU7MHn1ZPYd3uc7johI0DjnuGf8PWzcs5FeDXodW583Ki8Rpl/puYn+aUtY6Fi7IwdTDjJl9RTfUUREguat2W/x7tx3ebzZ41xbU1ec5mYqcBIWLqt4GeflPY/Ry0b7jiIiEhSfL/ucByY9QMdaHXmh1Qu+44hnKnASFqIjo+lRrwcjloxg456NvuOIiGSpyasmc+PoG2larikfd/5Yh0tFBU7Cx8OXPExKWgr/nvVv31FERLLM9v3buW7EddQuXptxPcYRkyfGdyTJAVTgJGxULVKV7nW7807iO+w8sNN3HBGRLPHT+p/Yd2QfA68eSJH8RXzHkRxCBU7CyqOXPsrew3t5I+EN31FERLJEwvoEIi2S+LLxvqNIDqICJ2GlQekGdInrwrPTn+Wt2W/5jiMics5mb5zN+SXPp2Cegr6jSA6iiXwl7Hx03UekpKVw74R7OZx6mL9c/BffkUREzkqaS2P2htl0i+vmO4rkMBqBk7CTLyofn3X9jLbV2/LCjBc4nHrYdyQRkbOycvtKdh3cRZNyTXxHkRxGBU7CUlREFP0b92fnwZ1MXjXZdxwRkbMye8NsAJqWb+o5ieQ0KnASttpUa0ORfEX4dMmnvqOIiJyVhA0JxOSJoU7xOr6jSA6jAidhK09kHq6vcz1f/vIlB44c8B1HROSMzd4wm/iy8URGRPqOIjmMCpyEtRvOv4G9h/cyfuV431FERM7I6p2rmZ88n6bldPhU/kwFTsJa88rNKVmwJO/OfZc0l+Y7jojIae06uIvWH7am2sBqpLpU2tds7zuS5EAqcBLWoiKiePjih5m0ahJ3jb1LJU5Ecrz3573Pt79+y7MtnmXVfau4tOKlviNJDqR54CTsPXzJw+w+tJsB3w/AOceg9oN0I2gRyZGccwyeM5iLy1/MU82f8h1HcjCvBbxVNeAAACAASURBVM7MrgbeACKBIc65l4/b/hegD5ACbAV6O+fWBralAosCu/7mnOuQbcElpJgZz7d8HsN44fsXcDgGtx+sEiciOc70tdNZvn05H3T8wHcUyeG8FTgziwT+A7QB1gM/m9kY59zSDLvNA+Kdc/vNrB/wD6B7YNsB51zDbA0tIcvMeK7lc0RYBM/NSP86qN0gzMx3NBGRYwbNGUThfIXpVld3XpBT8zkC1wRIcs6tBjCzT4GOwLEC55ybmmH/WcBN2ZpQwoqZ8WzLZzmSdoSXfniJBqUa0L9Jf9+xREQA2Lx3M6OXjqZffD/yR+f3HUdyOJ/HkMoB6zIsrw+sO5nbgQkZlvOZWaKZzTKzTsEIKOHphVYv0K5mOx6Y9AA/rvvRdxwREQAe+uYhAP1hKZkSEicBmdlNQDzwaobVlZxz8UAP4N9mVu0kr+0bKHqJW7duzYa0ktNFWAQfXfcRlc6rxE2f36QrU0XEu4lJExm2aBiPN3ucmsVq+o4jIcBngdsAVMiwXD6w7g/MrDXwJNDBOXfo6Hrn3IbA19XANKDRid7EOTfYORfvnIsvUaJE1qWXkFY4X2GebfEsv+76lYT1Cb7jiEgutu/wPu4aexe1i9fmicue8B1HQoTPAvczUMPMqphZHuAGYEzGHcysETCI9PK2JcP6ImaWN/C8OHApGc6dE8mMdjXbkScyD6OXjfYdRURysQ/mf8Da3WsZ1G4QeaPy+o4jIcJbgXPOpQD3AJOAZcBI59wSM3vOzI5OCfIqEAN8ZmbzzexowasDJJrZAmAq8PJxV6+KnNZ5+c6jTdU2jFo6Cuec7zgikks459i2f9ux54PmDOLCMhdyeaXLPSeTUOJ1Hjjn3Hhg/HHrnsrwvPVJXvcjUC+46SQ36BLXhXErxzFn0xziy8b7jiMiucCY5WPoPLIzY24YQ9H8RVm0ZRGD2w32HUtCjO7EILlah1odiIqIYvTS0SpwIpItJiZNJM2lceuXt9KkXBMK5SnEjfVu9B1LQkxIXIUqEixF8xfliipX8NnSz3QYVUSyxbS102hQqgEHUw4yIWkCPev1JCZPjO9YEmJU4CTX6163O6t2rmLW+lm+o4hImNu8dzO/bPuFG8+/kcHtB1MsfzHN+yZnRQVOcr0ucV0oEF2AoQuG+o4iImFuxtoZADSv3Jwe9Xqw+eHNnF/yfM+pJBSdcYEzswgziw1GGBEfCuUtROc6nRmxZAQHUw76jiMiYWz62ukUjC7IhWUuBCAyItJzIglVmSpwZjbczGLNrCCwGFhqZo8EN5pI9rm1wa3sOriLMcvHnH5nEZGzNH3tdC6teCnRkdG+o0iIy+wIXJxz7negE+n3I60C3By0VCLZrGXllpQrVE6HUUUkaLbt38biLYtpXqm57ygSBjJb4KLNLJr0AjfGOXcE0CV7EjYiIyK5qf5NTEqaxPb9233HEZEwNG3NNAAVOMkSmS1wg4A1QEFghplVAn4PVigRHzrX6UyqS2VC0gTfUUQkDH208CNKx5SmSbkmvqNIGMhUgXPODXTOlXPOXePSrQVaBjmbSLaKLxtPmZgyx86Dm7ByAs0/aM6BIwc8JxORULdxz0bGrRhHrwa9dP6bZInMXsRQyszeM7MJgeU44NagJhPJZhEWQfua7ZmYNJGDKQd5ePLDzFg7g69XfO07moiEuA/mf0CqS6V3o96+o0iYyOwh1A9Iv+l82cDyCuCBYAQS8alDrQ7sObyH+yfcz9KtS4m0SD5a+JHvWCISwtJcGu/Ne48WlVtQo1gN33EkTGS2wBV3zo0E0gCccylAatBSiXjSqkorCkQXYPDcwdQoWoMHLnqAiUkT2bpvq+9oIhKipq2Zxuqdq+nTqI/vKBJGMlvg9plZMQJXnprZRcDuoKUS8SR/dH6urHYlAE9c9gS9GvYiJS2FEUtGeE4mIqFq5JKRFIwuSOc6nX1HkTASlcn9/gKMAaqZ2UygBNAlaKlEPLq/6f3ki8pHz3o9iY6MpkGpBny88GPuaXKP72giEmLSXBpfLf+KtjXakj86v+84EkYyVeCcc3PNrDlQCzBgeWAuOJGw06JyC1pUbnFsuUe9Hjw65VE27tlI2UJlT/5CEZHjzN4wm+S9yXSq1cl3FAkzmSpwZnbLcasuMDOccx8GIZNIjnJx+YsBWJC8QAVORM7Il798SVREFNfUuMZ3FAkzmT0HrnGGx2XAM0CHIGUSyVHOL3k+AAs3L/ScRERyspS0FHYd3PWHdV/+8iUtK7ekSP4inlJJuMrsRL73ZnjcAVwAxAQ3mkjOUCR/ESrEVmDRlkW+o4hIDvbaj69R6rVSfLHsCwB+2fYLy7cvp1NtHT6VrJfZixiOt4/0G9qL5Ar1S9XXCJyInNLMdTM5nHqYLp914baGtzF+5XgiLZIOtXTASrJeZs+B+5r/v3l9BBAHjAxWKJGcpl7JekxaNYnDqYfJE5nHdxwRyYEWb1lM+5rtOZx6mPfmvcfV1a/mr5f8lfKx5X1HkzCU2RG41zI8TwHWOufWByGPSI5Uv1R9UtJSWL5tOfVK1fMdR0RymD2H9rBm1xr6NOrDY80eY8u+LZQpVMZ3LAljmZ1GZHqwg4jkZPVL1QfSL2RQgROR4y3duhRIv+gpMiJS5U2C7pQXMZjZHjP7/QSPPWb2e3aFFPGtZrGaREdEs3DzQiYmTaTmmzX5ZdsvvmOJSA6xeMti4P+vWhcJtlMWOOdcIedc7AkehZxzsdkVUsS36Mho4krEMXPdTPqM6cPKHSt5ZPIjvmOJSA6xeMti8kflp0oRXd8n2SOz88ABYGYlzazi0UewQonkRPVK1WPmupls3LORnvV6MnbFWKasnuI7lojkAIu3LqZuybpE2Bn9WhU5a5n6N83MOpjZSuBXYDqwBpgQxFwiOU79kunnwd3f9H6GdBhC5cKVeeibh0hNS/WcTER8W7xlsQ6fSrbK7J8KzwMXASucc1WAK4BZQUslkgN1rduVexrfw/OtnidfVD5euuIlFm5eyNgVY31HExGPtu/fTvLeZM4voQIn2SezBe6Ic247EGFmEc65qUD8ub65mV1tZsvNLMnMHjvB9rxmNiKwPcHMKmfY9nhg/XIzu+pcs4icTuXClXnzmjeJyZN+E5IucV0oHVOa/83/n+dkIuLTkq1LAF3AINkrswVul5nFADOAYWb2Bul3YzhrZhYJ/AdoS/rEwDeaWdxxu90O7HTOVQf+BbwSeG0ccANQF7gaeDvw/USyTVREFDfXv5mxK8ayee9m33FExJOjd2lRgZPslNkC1xHYDzwITARWAe3P8b2bAEnOudXOucPAp4H3Of59hwaejwKuMDMLrP/UOXfIOfcrkBT4fiLZ6raGt5HqUvl44ce+o4hINtu2fxv9x/XnoW8eouJ5FSlbqKzvSJKLZLbA3QmUcc6lOOeGOucGBg6pnotywLoMy+sD6064j3MuBdgNFMvka0WCrk6JOlxU/iLen/8+zrnTv0BEwsb9E+9n8NzB9GrQixm9ZpA+viCSPTJb4AoB35jZ92Z2j5mVCmaorGRmfc0s0cwSt27d6juOhKHeDXuzdOtS5m6a6zuKiGSTw6mHGbtiLL0a9GJQ+0FUKlzJdyTJZTJV4Jxzzzrn6gL9gTLAdDM71wmwNgAVMiyXD6w74T5mFgWcB2zP5GuPZh/snIt3zsWXKFHiHCOL/Fmn2p0AmJg00XMSEcku09dM5/dDv9OhVgffUSSXOtMZB7cAyaSXqJLn+N4/AzXMrIqZ5SH9ooQxx+0zBrg18LwL8J1LP041BrghcJVqFaAGMPsc84iclRIFS9CodCO+Wf0NAFv3baX1h61Zvm2552QiEixjlo8hf1R+rqh6he8okktldiLfu81sGvAt6eeg3eGcq38ubxw4p+0eYBKwDBjpnFtiZs+Z2dE/ad4DiplZEvAX4LHAa5cAI4GlpF9U0d85p9lUxZsrq13JT+t+Yu/hvfxv/v/49tdvGblkpO9YIhIEzjnGrBjDldWupEB0Ad9xJJeKyuR+FYAHnHPzs/LNnXPjgfHHrXsqw/ODQNeTvHYAMCAr84icrTZV2/DKzFeYtmbasXnhZvw2w3MqEQmGhZsX8tvu33i6+dO+o0gulqkC55x73Mwizaxsxtc4534LWjKREHJpxUvJF5WPF79/kV+2/ULpmNL8uO5HjqQeIToy2nc8EclCXy3/CsO4tsa1vqNILpbZQ6j3AJuBycC4wEP3DxIJyBeVj+aVmvPT+p8oEF2AF1u9yP4j+3VlqkiYOXDkAO8kvkOrKq0oFRMyEzJIGMrsRQwPALWcc3Wdc/UCj3M6B04k3LSp2gaArnFduabGNQDMWKvDqCLhZPCcwSTvTebvl//ddxTJ5TJb4NaRPomuiJxEp9qdKB1Tmv6N+1MqphS1itXSeXAiYeTAkQO8PPNlWlZuSfPKzX3HkVwusxcxrAammdk44NDRlc6514OSSiQEVStajU0PbTq2fHmlyxm5ZCSpaalERuhWvSKh7PdDv/P37/5O8t5kRnQZ4TuOSKZH4H4j/fy3PKTfleHoQ0RO4vJKl7P70G4WbVnkO4qInIO3f36bcq+XY+Dsgdx4/o1cXuly35FEMn0V6rMAZlbAObc/uJFEwkOLyi0AmLByAg1LN/QbRkTOyoLkBdw/8X6aV2rOS1e8RONyjX1HEgEyfxXqxWa2FPglsNzAzN4OajKREFc+tjxNyzVl5FJN6CsSilLTUunzdR+K5i/KyK4jVd4kR8nsIdR/A1eRfgstnHMLAI0hi5xG97rdmZ88nxXbV/iOIiJnaGDCQBI3JjLw6oEUzV/UdxyRP8j0vVCdc+uOW6VbV4mcRte66TcS0W21RELLrzt/5W9T/0a7mu3oVreb7zgif5LpaUTM7BLAmVm0mT1M+v1LReQUyseWp1nFZoxYoqvWREKFc45+4/oRYRG8fc3bmJnvSCJ/ktkCdxfQHygHbAAaBpZF5DS6xXVj8ZbFLN261HcUEcmEYYuGMWnVJF664iUqnFfBdxyRE8pUgXPObXPO9XTOlXLOlXTO3eSc2x7scCLhoEtcFwxjxGKNwomEgldmvsKFZS6kX3w/31FETipT04iY2cATrN4NJDrnvsraSCLhpUyhMjSv3JyRS0fyTItndDhGJAfbc2gPS7Ys4enmT2sCbsnRMnsINR/ph01XBh71gfLA7Wb27yBlEwkb3et255dtv2hSX5EcLnFjIg5H0/JNfUcROaXMFrj6QEvn3JvOuTeB1kBt4DrgymCFEwkX19e5ngiL0GFUkRwuYUMCAE3KNfGcROTUMlvgigAxGZYLAkWdc6lkuDeqiJxYiYIlaFWlFSOXjsQ55zuOiJzE7A2zqV60uuZ9kxwvswXuH8B8M/ufmX0AzANeNbOCwJRghRMJJ93rdidpRxLT1kwDIGF9Aq0/bM2Y5WP8BhORYxI2JNC0nA6fSs6X2Xuhvmdm44GjY8pPOOc2Bp4/EpRkImGmc53OPDblMVp92IrGZRsfO9dm35F9dKjVwXc8kVxv/e/r2bhnow6fSkg45QicmdUOfL0AKAOsCzxKB9aJSCYVzV+Upf2X8kLLF9h/ZD93N76bv132N2atn0XSjiTf8URyvdkbZgNoBE5CwulG4B4C7gD+eYJtDmiV5YlEwljJgiV58vInefLyJ4H0v/gHfD+Ajxd+zDMtnvEbTiSXS1ifQHRENA1KN/AdReS0TlngnHN3BL62zJ44IrlL+djytKzSko8XfszTzZ/WHHEinqS5NL5b8x0NSzckX1Q+33FETut0h1D/muF51+O2vRisUCK5yc31b2bVzlXHpi8Qkez32JTHSNyYyO2NbvcdRSRTTncV6g0Znj9+3LarsziLSK7UuU5n8kfl542EN3xHEcmV3p3zLq/++Cp3x99N3wv7+o4jkimnK3B2kucnWhaRsxCbN5aHL3mYTxd/yg+//eA7jkiuMmX1FO4efzdXV7+aN9q+odMYJGScrsC5kzw/0bKInKXHmj1GhdgK3DvhXlLTUn3HEQlLS7cuZeOejX9Y7jKyC7WL12ZElxFERWRqZi2RHOF0/7Y2MLPfSR9tyx94TmBZZ3mKZJEC0QV47crX6D6qO+/Pe587LrzDdySRsHLgyAEav9uYw6mHua72daSkpTB59WQKRhdk7I1jic0b6zuiyBk55Qiccy7SORfrnCvknIsKPD+6HJ1dIUVyg65xXWlYuiHvzXvPdxSRsJO4MZH9R/ZzVbWrmLx6MnM2zaFnvZ58d+t3VCpcyXc8kTPmZbzYzIoCI4DKwBqgm3Nu53H7NATeAWKBVGCAc25EYNsHQHNgd2D3Xs65+dmRXSRYzIzudbvz+LeP89vu36h4XkXfkUTCxtHzS4d2GnrsPqc6301CWWbvhZrVHgO+dc7VAL4NLB9vP3CLc64u6Ve8/tvMCmfY/ohzrmHgofImYeH6OtcD8Pmyzz0nEQkvM9fNpHbx2hQrUAwzU3mTkOerwHUEhgaeDwU6Hb+Dc26Fc25l4PlGYAtQItsSinhQo1gN6peqz6ilo3xHEQkbaS6NH9f9yKUVLvUdRSTL+CpwpZxzmwLPk4FSp9rZzJoAeYBVGVYPMLOFZvYvM8sbpJwi2a5LnS7MXDeTDb9v8B1FJCz8su0Xdh7cqQInYSVoBc7MppjZ4hM8OmbczznnOMWUJGZWBvgIuM05lxZY/ThQG2gMFAUePcXr+5pZopklbt269Vx/LJGg6xLXBdBhVJFzkZqWSrfPujFm+Rhm/jYTgEsrqsBJ+LD0/pTNb2q2HGjhnNsUKGjTnHO1TrBfLDANeNE5d8JjSmbWAnjYOdfudO8bHx/vEhMTzym7SHaIHxzPut/XseCuBZSOKe07jkjIWbF9BbXeqkXeyLzUK1WPtbvWsvnhzTr3TUKKmc1xzsWfaJuvQ6hjgFsDz28Fvjp+BzPLA3wBfHh8eQuUPiz9v8ROwOKgphXJZkM7DWXPoT3c9PlNpKalkpKWgo8/tkRC1dKtS4H0ORYTNyZySYVLVN4krPgqcC8DbcxsJdA6sIyZxZvZkMA+3YDLgV5mNj/waBjYNszMFgGLgOLAC9kbXyS46pasy5tt3+TbX7+lwr8qkO+FfHQf1d13LJGQcbTATb11KlUKV+G62td5TiSStbzMA+ec2w5ccYL1iUCfwPOPgY9P8vpWQQ0okgP0btSbrfu3Mj95PgdSDvDZ0s/44bcfaFaxme9oIjne0q1LqRBbgQalG7DqvlUafZOwoxu/ieRQZsZjzdKnSNx/ZD813qzBXyf/lZm9Z+qXkchpLNu2jLgScYAm7JXw5OsQqoicgQLRBXi2xbP8tP4nvlr+p1NGRSSDNJfGsq3/X+BEwpEKnEiI6NWwF7WL1+apqU/pggaRU1i7ay0HUg5Qp3gd31FEgkYFTiREREVE8dilj7FoyyK+WfWN7zgiOdbRCxg0AifhTAVOJITcWO9GyhYqy2s/veY7ikiOdbTA1SmhETgJXypwIiEkT2Qe7m96P1NWT2F+8nzfcURypKXbllI6pjRF8xf1HUUkaFTgREJM3wv7EpMnhkcmP8Lug7t9xxHxJu3Y3RX/aOnWpTp8KmFPBU4kxBTOV5h/tP4HU3+dSsNBDZm1fpbvSCLZ7tedvxLzYgzfr/3+D+udc+lXoBZXgZPwpgInEoL6Ne7H97el/+JqN7wduw7u8pxIJHv98NsPHEg5wLBFw/6wfvn25ew5vId6pep5SiaSPVTgRELUxRUuZnS30ew4sIMXv3/RdxyRbHX0HNCvV3z9h0OpX/2SPk9i2+ptveQSyS4qcCIh7IIyF3Brw1t5I+ENVu9c7TuOSLaZv3k+hrFxz0bmbJxzbP2Xy78kvmw8Fc6r4DGdSPCpwImEuAGtBhAVEcVD3zykCX4lV3DOMT95Pp3rdCbCIhizfAwAm/ZsYtb6WXSq1clzQpHgU4ETCXFlC5XlmebP8OUvX/LW7Ld8xxEJuvW/r2fHgR1cUeUKmlVsxpgV6QXuaJHrVFsFTsKfCpxIGHjokodoX7M9f/nmL8z8babvOCJBdfT8t4alG9KhZgcWbl7ItDXT+PyXz6letLqmEJFcQQVOJAxEWAQfXvchlc6rRO8xvX3HEQmq+cnp57/VK1WP6+Oup1CeQrQc2pJvVn1Dp1qdMDPfEUWCTgVOJEwUzleYBy56gBXbV/Drzl99xxEJmvmb51OjWA1i8sRQuXBlku5LYminofSL70f/Jv19xxPJFlG+A4hI1mlZuSUAU9dMpUqRKp7TiATH/OT5xJeNP7ZcsmBJbmlwC7c0uMVjKpHspRE4kTASVyKOkgVL8t2v3/mOIhIUuw/uZvXO1TQs1dB3FBGvVOBEwoiZ0apKK6aumaopRSQs/fDbDwA0KtPIcxIRv1TgRMJMy8ot2bhnIyu2r+CLZV/wwMQHVOYkbAyaM4hSBUvRqkor31FEvNI5cCJh5ugvttd/ep2hC4ZyKPUQbau35arqV3lOJnJu1uxaw9gVY3nysifJE5nHdxwRrzQCJxJmqhWpRvnY8gyeO5iSBUtStlBZXp75su9YIuds8JzBmBl9L+zrO4qIdypwImHGzGhTtQ15IvMwqtsoHr74Yaatmcas9bN8RxM5a4dSDjFk7hA61Oqg+5yKoAInEpZebfMqc/rOoUm5Jtxx4R0UyVeEl3/QKJyErk8Xf8rW/VvpF9/PdxSRHEEFTiQMFStQjPNLng9ATJ4Y+jfuz1fLv+K33b95TiZy5tJcGq/MfIX6perTpmob33FEcgQVOJFcoFfDXgB8sugTv0FEzsLXy79m2bZlPHbpY7pNlkiACpxILlCtaDUuLn8xHy38SFOKSEhxzvHSDy9RpXAVutbt6juOSI6hAieSS9xc/2aWbF3Cgs0LfEcRybSJSRNJ2JDAI5c8QlSEZr4SOcpLgTOzomY22cxWBr4WOcl+qWY2P/AYk2F9FTNLMLMkMxthZpoQSOQ0utXtRnRENB8v/Nh3FJGTcs6xbvc6IP22WXeOvZPaxWtzW6PbPCcTyVl8jcA9BnzrnKsBfBtYPpEDzrmGgUeHDOtfAf7lnKsO7ARuD25ckdBXrEAxrqlxDcMWDeNgykHfcURO6KOFH1Hx3xW56fObuHv83WzYs4EPOn5Avqh8vqOJ5Ci+ClxHYGjg+VCgU2ZfaOlnsLYCRp3N60Vys3ub3Evy3mQGzBjgO4rICX21/Cti8sQwcslIhi8aziOXPELT8k19xxLJcXwVuFLOuU2B58lAqZPsl8/MEs1slpkdLWnFgF3OuZTA8nqgXBCzioSNK6pewc31b+blmS+zaPMi33FE/iAlLYVvV3/LDXVvYN6d8xjQagDPtHjGdyyRHCloZ4Sa2RSg9Ak2PZlxwTnnzOxkl8VVcs5tMLOqwHdmtgjYfYY5+gJ9ASpWrHgmLxUJS69f9ToTkiZwx9d38OPtPxJhupZJcoaE9QnsPrSbq6pfRd2Sdalbsq7vSCI5VtD+z+2ca+2cO/8Ej6+AzWZWBiDwdctJvseGwNfVwDSgEbAdKGxmR8tneWDDKXIMds7FO+fiS5QokWU/n0ioKl6gOC+0fIGEDQnM2TjHdxyRY75Z9Q0RFsEVVa7wHUUkx/P1p/cY4NbA81uBr47fwcyKmFnewPPiwKXAUpc+idVUoMupXi8iJ9clrgsRFsGY5WNOv7NINpm0ahJNyjWhSP4TTkwgIhn4KnAvA23MbCXQOrCMmcWb2ZDAPnWARDNbQHphe9k5tzSw7VHgL2aWRPo5ce9la3qREFesQDGaVWzGV8v1t4/kDDsO7ODnjT9zVbWrfEcRCQleZkV0zm0H/jRG7pxLBPoEnv8I1DvJ61cDTYKZUSTcdazVkYe+eYhfd/5KlSJVfMeRXGz1ztU8P+N50lyaCpxIJunsZZFcqn3N9gB8veJrz0kkN/tg/gdUH1idjxZ8xO2NbqdJOf1tLpIZui+JSC5Vo1gN6hSvw7BFw1iyZQnT1k6jX3w/7m58N3kidXMTCb6t+7by4KQHubTipXx6/aeUi9WMUCKZpRE4kVysY62OzN4wmw8WfEBs3lgenPQgDf/bkO37t/uOJrnA498+zt7DexncbrDKm8gZUoETycUevPhBXr7iZVbdt4rZfWbzebfPWbZtGUPmDjn9i0XOQcL6BN6b9x4PNH2AOiXq+I4jEnIsfVaO3CE+Pt4lJib6jiGSo7Uc2pI1u9aQdG8SkRGRvuNIGEpNS6XpkKZs2ruJX/r/QqG8hXxHEsmRzGyOcy7+RNs0Aicif3B3/N2s2bWGiUkTfUeRMPXu3HeZs2kO/7zynypvImdJBU5E/qBT7U6UjinN24lv+44iYeTjhR9z0ZCLGDBjAE98+wQtK7eke93uvmOJhCwVOJH/a+/Oo6uq7r+Pv78ZIIBMGkRkkIAIOD3MFkFFJpEUUaQqDj+sNCiDiNRStbW0jxP6aEF/ikKhVVRQilJFQEADCgrKJLMKCgoICKhVDGPyff64R5sigUCGk3Pzea11V+7ZOfeezz2LFb5373P2lv+SnJhM32Z9mbFuBhu+2RB2HIkTf1v6Nz7c9iF/nPNHvt//Pf976f9iZmHHEoksFXAi8jN9mvXBcSasnBB2FIkDWQeyWLBpAYPOG/TTDTNaqF6kYFTAicjP1Klch7Z12vLCyhcoTTc6SdF494t3OZBzgIvrXky9qvVoWqNp2JFEIk8FnIgc1rVnX8vanWtZsX1F2FEk4uZsnENSQhJt67QNO4pI3FABJyKH9auzfkVSQpKGUaXAMjdk0vLUlrrjVKQQqYATkcNKLZ9K5/qdmbhqIjmeE3Yciajv9n3H4i8X0z6tfdhRROKKCjgRydO1Z1/Lpu82MWfDnLCjSETN+3we2Z7NxXUvDjuKSFxRAScieerRuAfVyldjxMIRYUeRiPrnmn9SJrEM59c+P+wo2yIpIQAAFl9JREFUInFFBZyI5KlccjkGthrItHXTWLNjTdhxJGLGLh3Ls8ufZWDLgZRLLhd2HJG4ogJORI6of8v+pCSl8NcFfw07ikTI2xvfpv+0/nSu35mHOj0UdhyRuKMCTkSOKLV8Kr9u8mueW/EcW7/fGnYciYCZ62fSdUJX6p9Yn5d6vkRSQlLYkUTijgo4ETmqIa2HkJ2TzQPzHgg7ipRwr370Kt0mdqPBiQ2Y23suVVKqhB1JJC6pgBORozr9xNPJaJbB00ueZv3X68OOIyXUngN7uGXaLZxT/Rzm3jiX6idUDzuSSNxSASci+TKs3TDKJJbhD5l/CDuKlFCjl4xm2+5tjLhkhHreRIqYCjgRyZdTTjiF37b+LZNWT+KsUWfRbHQzlm1dFnYsKSGyDmQxfP5w2qe158LTLgw7jkjcUwEnIvk2tM1QBrYcSOPUxqzduZa/Lf1b2JGkhBi9eDTbf9jOny/6c9hRREoFc/ewMxSbFi1a+OLFi8OOIRIXek7qybub3mXLkC0kmL4Llmb7s/eT9lgaDU9qSGbvzLDjiMQNM1vi7i0O9zv91RWR49KjcQ+27d7Gws0Lw44iIZu4ciJffv8lv2/z+7CjiJQaKuBE5LikN0gnOSGZV9a+EnYUCZG788iCRzj75LPpXL9z2HFESg0VcCJyXCqnVKZjvY68svYVStOlGBLz6defkrkhkwkrJ7Dqq1Xc0foOzCzsWCKlhqbHFpHj1qNxDzKmZrDoy0W0qtkq7DhSTMYuHcuA6QPYn70fgFMrnkqvc3qFnEqkdAmlB87MTjSz2Wa2LvhZ9TD7XGxmH+Z67DWzy4PfPWNmG3L9rknxfwoR6d6wO5XLVqbD+A48+cGTZOdkhx1JilCO59B/Wn8ypmZw0WkXMev6WYz+5WimXD2FMollwo4nUqqEcheqmT0MfO3uw83sTqCqu+d59auZnQisB2q5e5aZPQO87u6Tj+W4ugtVpPBt+GYDt0y7hVmfziKtShr9W/ZnYKuBpCSlhB1NCpG7M/iNwTz+wePc0foOhnccTmJCYtixROJaSbwLtTvwbPD8WeDyo+zfE5jh7llFmkpEjlla1TTeuO4NXr7qZWpVqsXvZv+OYXOGhR1LCkmO57D4y8X0n9afxz94nNt/cTsPd3pYxZtIyMLqgfvW3asEzw345sftPPbPBP7q7q8H288ArYF9wFvAne6+72jHVQ+cSNHr8VIP5n8xn81DNmtYLQ50GN+BzA2xud36tejHk12f1M0KIsUklB44M3vTzFYd5tE9934eqyDzrCLNrAZwDjAzV/NdQCOgJXAicKTh175mttjMFu/YsaMgH0lE8iGjWQY7snYw9eOpYUeRAvrsm8/I3JDJoFaD2H7Hdkalj1LxJlJCFFkB5+4d3f3swzxeBbYHhdmPBdpXR3irq4Ap7n4g13tv9Zh9wD+APG9/c/cx7t7C3VtUq1atcD6ciOSpc/3O1KpUi7HLxoYdRQpo9qezAejXsh8nVzg55DQikltY18C9BvQOnvcGXj3Cvr2AibkbchV/Ruz6uVVFkFFEjkNiQiI3NbmJmetnsuGbDWHHkQKY/dlsaleqTcOTGoYdRUQOEVYBNxzoZGbrgI7BNmbWwsx++tpuZnWB2sDbh7z+BTNbCawEUoH7iiGziOTTTU1vAqDe4/Wo8EAFmo9pTr/X+/HRzo9CTib5lZ2TTeaGTDrV66RhU5ESKJSJfN19F9DhMO2Lgd/k2t4I1DzMfu2LMp+IFMxpVU5jaq+pLN++nJ1ZO1mxfQXjV4xnzsY5LL9lOWWTyoYdUY5iydYlfLP3GzrV7xR2FBE5DK3EICJFIv2MdNLPSP9pe8a6GXSd0JVHFzzK3Rfczdodazmp/Em6tqqEmvXpLAA6pP3su7aIlABaC1VEisWlDS7lysZXcu8793LN5Gs4c9SZ9JzUM+xYkodZn86i6SlNqVZBN3+JlEQq4ESk2IzsMpKkhCT+9dG/aFe3HfO+mMfiLzU3Y0ni7vzhrT8w74t59GjcI+w4IpIHFXAiUmxqVarFooxFrLt1Ha9e8yoVy1RkxMIRYceSXAbNGMQD8x8go1kGd7W9K+w4IpIHFXAiUqwapTaiduXaVCpbiT5N+zBp9SS2fLcl7FgCZG7I5IlFT3Dbebcx+pejtVyWSAmmAk5EQjPovEHkeA7D5g4jOyc77Dilyp4De8i9lGKO5zB09lBqV6rN8I7DNXWISAmnAk5EQpNWNY0BLQcwbtk4OozvwKZ/bwo7Uqmw+qvVpP6/VNo9245lW5cBMHnNZJZsXcK9F99LSlJKyAlF5GhCWcw+LFrMXqTkcXfGLx/PwBkDKZNYholXTqRz/c5hx4pb2TnZtP1HWz7e+TGJCYnsytpF9ROqs3v/btKqpLHs5mUaOhUpIUJZzF5EJD/MjN5NerO071JOrXgqXZ7vwkPzH6I0fbksTk8tfoqFmxfy+KWPs+7Wddzf/n66ndGNrg26MqbbGBVvIhGhHjgRKTF+2P8Dv5n6G15c9SKDzxvMo5c8SoLpe2Zh2fLdFho92Yg2tdsw47oZus5NpIQ7Ug+cVmIQkRKjQpkKvNDjBapXqM7I90eye/9uxnQbo0KjkAx9cygHsg8wKn2UzqlIxKmAE5ESJcESGHHJCMonl+fB+Q/S5JQmDGg1IOxYkTf/i/lMWDmBP17wR+pVrRd2HBEpII1NiEiJY2bc1/4+0hukM3jmYN7b9F7YkSIt60AWt864ldqVanPXBZqcVyQeqIATkRIpwRJ47ornqFO5DukT0pm0elLYkSLpzc/e5JynzuHDbR8ysstIyieXDzuSiBQCFXAiUmJVLVeVN294k4YnNeTqyVczYNoA3Z16DB557xE6PdeJREtkTu85WttUJI6ogBOREi2tahrzfj2PwecNZtTiUby46sWwI5V4u7J2MXT2UH43+3dcddZVLL9lOe3qtgs7logUIt3EICIlXnJiMo90foSFWxYyYPoA2tVtR42KNcKOVeJMWTuFO9+6k092fQLAzc1v5smuT2puN5E4pB44EYmExIREnun+DHsO7qHds+04f9z59P5Xbw5kHwg7Wuj2HtxL36l96TGpB+WSyjG8w3DeufEdnkp/SsWbSJxSD5yIREbD1IaM7TaWx95/jKSEJMYvH0/9qvX500V/CjtaaNydPq/1YcLKCdzZ5k7+cvFfKJNYJuxYIlLEtBKDiETWda9cx6TVk1iUsYgmpzQJO04oHpz3IHdn3s0D7R/QFCEiceZIKzGogBORyNqVtYuzRp3FwZyDVCxbkZMrnMxr17xG9ROqhx2tSG3bvY2RC0eycPNC3v78ba4951qev+J5ra4gEme0mL2IxKWTyp/E5Ksmc+FpF9K2TltWbl/J1ZOv5mDOwbCjFZn1X6/n/HHn8+iCR9lzcA9Dzx/K2G5jVbyJlDK6Bk5EIq1tnba0rdMWgEvqX8INU25gwLQBXH/u9dSrWo+alWqGnLDwLPlyCekT0jmYc5D3bnqPljVbhh1JREKiAk5E4sb1517P+5vf54lFTzBm6RiSE5J5+8a3aV27ddjRjsveg3v5YMsH1KxYk7kb5zJg+gCqn1CdOdfNoXG1xmHHE5EQ6Ro4EYkr7s7Kr1aybfc2MqZmkJKUwoc3f0i55HJhRzsm2TnZXPL8Jby14a2f2jrW68jEKyeSWj41xGQiUlyOdA2ceuBEJK6YGedWP5dzq5/LuMvG0em5TgyZOYSuDbqyL3sflze6nKSEkv+n75459/DWhre4v/391DihBmbGDefeoHndRARQD5yIxLl+r/fj6SVP/7R9Y5MbGXfZOBLsOO/hcocPPoBHHoHp02HPHihXDtLT4Y47oGVLOI4bCj7/9nMenP8gCzYvoFLZSsz/Yj4ZzTIY023M8eUUkcjTNCIBFXAipc/eg3uZvm46NSvW5PVPXue+efdxS/Nb6FCvA3sP7qXL6V3yPyR54AD8z//Aa6/B3r2Qk/Of3yUkxAq5bt1g/HhITs53xofmP8Q9c+7BzOiQ1oGsA1nUqVyHMd3GkJKUcoyfWETihQq4gAo4kdLN3Rkycwgj3x/5U1tKUgpXNLqC/dn72bp7K9efcz0ZzTN+PszqDtdeGyvesrLyPki5ctC9O0yYcNieOHdnzJIxfP7vz7ntvNuYuGoit8+8nZ5n9mTEJSOoValWYX1cEYm4ElfAmdmvgD8DjYFW7n7YqsrMugCPAYnAWHcfHrSnAS8CJwFLgBvcff/RjqsCTkTcnRXbV5BgCezL3sfYpWOZvGYyqeVTKZtUlhXbV3BWtbPo07QPXU7vQqPURrE51t5/Hzp0gB9+OPpBKlSAzExo1eq/mr/f9z03v34zE1dNBGLF496De+nRuAcv9XwpEtfmiUjxKYkFXGMgBxgN3HG4As7MEoFPgE7AZmAR0Mvd15jZJOAVd3/RzJ4Glrv7U0c7rgo4ETkSd2fKR1MYNncYq75aBUBq+VTOq3kej47bzBlzV2K5h03zkJNgrLnoTIb2qcPOrJ04zs6snXzx7y8AuO/i++jRuAfD3x3OgewDjLtsHGWTyhbpZxOR6ClxBdxPBzebS94FXGvgz+5+SbD94yJ/w4EdwCnufvDQ/Y5EBZyI5NfGbzcy+9PZLNi8gAWbF7D49o+ocCD/r9+dDK0fO5talWphGFVSqtDgxAZ0Ob1LZOelE5HiFdVpRGoCm3JtbwbOIzZs+q27H8zVnudU62bWF+gLUKdOnaJJKiJxp26VumQ0zyCjeQYAfmsCkP8vvBUOGiv7rSyidCJS2hXZWqhm9qaZrTrMo3tRHfNw3H2Mu7dw9xbVqlUrzkOLSByxcsc2EbCVL19ESUREirAHzt07FvAttgC1c23XCtp2AVXMLCnohfuxXUSk6KSnw8sv//fUIXlJSIjtLyJSRIqsB64QLAIamFmamZUBrgFe89hFe3OAnsF+vYFXQ8ooIqXFb38bmyIkP1JSYvuLiBSRUAo4M7vCzDYDrYFpZjYzaD/VzKYDBL1rA4GZwFpgkruvDt7i98AQM1tP7Jq4ccX9GUSklGnVKjZJ79GKuHLl4LLLYisyiIgUEU3kKyKSXz+uxDB1amwJrUNXYkhJiRVvx7gSg4jI4RzpLtSSPIQqIlKyJCfHVljIzIQrr4xN2JuQEPvZsyfMnQsTJ6p4E5EiV5KnERERKXnMYsOpkyaFnURESjH1wImIiIhEjAo4ERERkYhRASciIiISMSrgRERERCJGBZyIiIhIxKiAExEREYkYFXAiIiIiEVOqVmIwsx3A50V8mFRgZxEfI97pHBaczmHB6PwVnM5hwekcFlzUz+Fp7l7tcL8oVQVccTCzxXkteyH5o3NYcDqHBaPzV3A6hwWnc1hw8XwONYQqIiIiEjEq4EREREQiRgVc4RsTdoA4oHNYcDqHBaPzV3A6hwWnc1hwcXsOdQ2ciIiISMSoB05EREQkYlTAFSIz62JmH5vZejO7M+w8UWNmfzezr8xsVdhZosjMapvZHDNbY2arzey2sDNFjZmlmNkHZrY8OId/CTtTFJlZopktM7PXw84SRWa20cxWmtmHZrY47DxRZGZVzGyymX1kZmvNrHXYmQqbhlALiZklAp8AnYDNwCKgl7uvCTVYhJjZhcBuYLy7nx12nqgxsxpADXdfamYVgSXA5fo3mH9mZkAFd99tZsnAfOA2d18YcrRIMbMhQAugkrv/Muw8UWNmG4EW7h7l+ctCZWbPAvPcfayZlQHKu/u3YecqTOqBKzytgPXu/pm77wdeBLqHnClS3P0d4Ouwc0SVu29196XB8++BtUDNcFNFi8fsDjaTg4e+5R4DM6sFpANjw84ipZOZVQYuBMYBuPv+eCveQAVcYaoJbMq1vRn95ykhMbO6QFPg/XCTRE8w/Pch8BUw2911Do/NSGAokBN2kAhzYJaZLTGzvmGHiaA0YAfwj2Aof6yZVQg7VGFTAScSZ8zsBOBlYLC7fxd2nqhx92x3bwLUAlqZmYbz88nMfgl85e5Lws4ScW3dvRlwKTAguLxE8i8JaAY85e5NgR+AuLsuXQVc4dkC1M61XStoEyk2wXVbLwMvuPsrYeeJsmDIZQ7QJewsEdIGuCy4hutFoL2ZPR9upOhx9y3Bz6+AKcQu0ZH82wxsztV7PplYQRdXVMAVnkVAAzNLCy6YvAZ4LeRMUooEF+CPA9a6+1/DzhNFZlbNzKoEz8sRuynpo3BTRYe73+Xutdy9LrG/gZnufn3IsSLFzCoENyERDPt1BnRn/jFw923AJjNrGDR1AOLuZq6ksAPEC3c/aGYDgZlAIvB3d18dcqxIMbOJQDsg1cw2A8PcfVy4qSKlDXADsDK4hgvgbnefHmKmqKkBPBvcVZ4ATHJ3TYUhxak6MCX2fYwkYIK7vxFupEi6FXgh6FD5DPh1yHkKnaYREREREYkYDaGKiIiIRIwKOBEREZGIUQEnIiIiEjEq4EREREQiRgWciIiISMSogBORyDOzy83MzaxRAd7jRjN7Ip/H+tPxHid4j7lm1iJ4Pv3HueeO430GmtlNBckiItGkAk5E4kEvYH7ws6gNBUYd2mhmxzWvprt3LcBC238nNt+ViJQyKuBEJNKCtV/bAn2Izf6PmbULerkmm9lHZvZCsFIFZtY1aFtiZo+b2c8m6g1WZHjZzBYFjzZB+xnAPnffGWw/Y2ZPm9n7wMNm1srMFgQLaL/340zwZlbOzF40s7VmNgUol+tYG80sNXj+ryDX6tyLmJvZbjO738yWm9lCM6sO4O5ZwEYz01JLIqWMCjgRibruwBvu/gmwy8yaB+1NgcHAmUA9oI2ZpQCjgUvdvTlQLY/3fAwY4e4tgSuBsUF7G2DpIfvWAs539yHElt26IFhA+0/AA8E+/YAsd28MDAOac3g3BblaAIPM7KSgvQKw0N3/D/AOkJHrNYuBC/J4PxGJU1pKS0SirhexggtiC6j3Al4HPnD3zQDB0mJ1gd3AZ+6+Idh/ItCXn+sInBl02gFUCnr6agA7Dtn3n+6eHTyvTGwprgaAA8lB+4XA4wDuvsLMVuTxWQaZ2RXB89pAA2AXsD/4TABLiK3R+qOvgOO+9k9EokkFnIhElpmdCLQHzjEzJ7YOsQPTgH25ds3m2P7eJQC/cPe9hxxvD7EiLbcfcj2/F5jj7leYWV1gbn4PaGbtiBWOrd09y8zmAinBrw/4f9Y9PPSzpAB78nscEYkPGkIVkSjrCTzn7qe5e113rw1sIO8hxY+BekFxBXB1HvvNItfNAWbWJHi6Fjj9CHkqA1uC5zfman8HuDZ4r7OBc/N47TdB8dYI+MURjpPbGcCqfO4rInFCBZyIRFkvYMohbS+Tx92o7r4H6A+8YWZLgO+Bfx9m10FACzNbYWZrgFuC9neAppZrbPUQDwMPmtky/ruX7CngBDNbC/xfYsOgh3oDSAr2GQ4szOMYh2oDzM7nviISJ+w/vfIiIvHPzE5w991BEfYksM7dRxzD6x8Dprr7m0UWMv9ZmgJD3P2GsLOISPFSD5yIlDYZwU0Nq4kNW44+xtc/AJQv9FTHJxW4J+wQIlL81AMnIiIiEjHqgRMRERGJGBVwIiIiIhGjAk5EREQkYlTAiYiIiESMCjgRERGRiFEBJyIiIhIx/x96T+oEEr7o0gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,6))\n",
    "plt.plot(angles, eigenvalues, 'g')\n",
    "plt.plot(angles[eigenvalues.index(-1.0)], min(eigenvalues), \"ro\", markersize=12)\n",
    "plt.title(\"Eigenvalue vs angles\")\n",
    "plt.xlabel(\"Angle(radian)\")\n",
    "plt.ylabel(\"Eigenvalues\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Using `scipy` minimize function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lowest eigenvalue:  -1.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "     fun: -1.0\n",
       "   maxcv: 0.0\n",
       " message: 'Optimization terminated successfully.'\n",
       "    nfev: 17\n",
       "  status: 1\n",
       " success: True\n",
       "       x: array(9.40225711)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# random initial choice\n",
    "phi = np.random.uniform(0, 2 * np.pi)\n",
    "\n",
    "res = minimize(cost, phi, method='COBYLA')\n",
    "\n",
    "print('Lowest eigenvalue: ', res.fun)\n",
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. VQE in Qiskit\n",
    "\n",
    "We use an already built VQE instance in Qisit Aqua"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lowest eigenvalue:  (-1+0j)\n"
     ]
    }
   ],
   "source": [
    "optimizer = SPSA()\n",
    "\n",
    "var_form, _ = ansatz()\n",
    "vqe = VQE(operator=H, var_form=var_form, optimizer=optimizer)\n",
    "vqe_results = vqe.run(QuantumInstance(backend=Aer.get_backend('qasm_simulator')))\n",
    "\n",
    "print(\"Lowest eigenvalue: \", vqe_results.eigenvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PyCharm (Qiskit)",
   "language": "python",
   "name": "pycharm-f8797ed6"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}