39 - vqe.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "805ba3c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "from functools import reduce\n",
    "from itertools import product\n",
    "import numpy as np\n",
    "import random\n",
    "import math\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt\n",
    "from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, Aer, execute\n",
    "\n",
    "random.seed()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "3f0ef083",
   "metadata": {},
   "outputs": [],
   "source": [
    "PAULIS = {\n",
    "    \"I\" : np.eye(2, dtype = complex),\n",
    "    \"X\" : np.array([[0, 1], [1, 0]], dtype = complex),\n",
    "    \"Y\" : np.array([[0, -1j], [1j, 0]], dtype = complex),\n",
    "    \"Z\" : np.array([[1, 0], [0, -1]], dtype = complex)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "85b28488",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pauli_decomposition(H):\n",
    "\n",
    "    n = int(math.log2(len(H)))\n",
    "    dims = 2 ** n\n",
    "\n",
    "    components = defaultdict(int)\n",
    "\n",
    "    keys = [\"\".join(k) for k in product(PAULIS.keys(), repeat = n)]\n",
    "\n",
    "    for i, val in enumerate(product(PAULIS.values(), repeat = n)):\n",
    "\n",
    "        basis_mat = reduce(np.kron, val)\n",
    "\n",
    "        coeff = np.dot(H.reshape(-1), basis_mat.reshape(-1)) / dims\n",
    "        coeff = np.real_if_close(coeff).item()\n",
    "\n",
    "        if not np.allclose(coeff, 0):\n",
    "            components[keys[i]] = coeff\n",
    "            \n",
    "    return components\n",
    "      \n",
    "def ansatz(qc, qr, thetas, rounds = 3):\n",
    "    for r in range(rounds):\n",
    "        for i in range(len(qr)):\n",
    "            offset = r * len(qr) * 2 + i * 2\n",
    "            \n",
    "            qc.ry(thetas[offset], qr[i])\n",
    "            qc.rz(thetas[offset + 1], qr[i])\n",
    "            \n",
    "            if i > 0 and r < rounds - 1:\n",
    "                qc.cx(qr[i - 1], qr[i])\n",
    "    return qc\n",
    "\n",
    "def matrix_decomposition_circuits(H, add_ansatz):\n",
    "    \n",
    "    n = int(np.log2(len(H)))\n",
    "    dims = 2 ** n\n",
    "    \n",
    "    if H.shape != (dims, dims):\n",
    "        raise ValueError(\"The imput must be a 2^n x 2^n dimensional matrix.\")\n",
    "        \n",
    "    components = pauli_decomposition(H)\n",
    "    \n",
    "    circuits = {}\n",
    "    for (paulistring, weight) in components.items():\n",
    "        paulis = list(enumerate(reversed(paulistring)))\n",
    "        \n",
    "        qr = QuantumRegister(n, 'qr')\n",
    "        cr = ClassicalRegister(1, 'cr')\n",
    "        qc = QuantumCircuit(qr, cr)\n",
    "        \n",
    "        add_ansatz(qc, qr)\n",
    "        \n",
    "        if all(x[1] == \"I\" for x in paulis):\n",
    "            circuits[paulistring] = (qc, weight)\n",
    "            continue\n",
    "            \n",
    "        for idx, letter in paulis:\n",
    "            if letter == \"X\":\n",
    "                qc.h(qr[idx])\n",
    "            elif letter == \"Y\":\n",
    "                qc.h(qr[idx])\n",
    "                qc.sdg(qr[idx])\n",
    "                \n",
    "        for idx, (first, second) in enumerate(zip(paulis, paulis[1:])):\n",
    "            if second[1] == \"I\":\n",
    "                qc.swap(first[0], second[0])\n",
    "            elif first[1] != \"I\":\n",
    "                paulis[idx + 1] = (paulis[idx + 1][0], paulis[idx][1])\n",
    "                qc.cx(first[0], second[0])\n",
    "                \n",
    "        qc.measure(qr[-1], cr[0])\n",
    "        circuits[paulistring] = (qc, weight)\n",
    "        \n",
    "    return circuits\n",
    "\n",
    "def run_circuit(circuit, backend, shots = 10000):\n",
    "    result = execute(circuit, backend, shots = shots).result()\n",
    "\n",
    "    counts = result.get_counts()\n",
    "\n",
    "    return (\n",
    "        (counts['0'] if '0' in counts.keys() else 0) - (counts['1'] if '1' in counts.keys() else 0)\n",
    "    )/shots\n",
    "\n",
    "def expectation(H, rounds, thetas):\n",
    "    backend = Aer.get_backend('statevector_simulator')\n",
    "    circuits = matrix_decomposition_circuits(\n",
    "        H,\n",
    "        lambda qc, qr: ansatz(qc, qr, thetas, rounds)\n",
    "    )\n",
    "    \n",
    "    return sum(map(lambda key: run_circuit(circuits[key][0], backend) * circuits[key][1], circuits.keys()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "c99b8662",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an optimizer, SPSA\n",
    "class SPSA:\n",
    "    \n",
    "    def __init__(self, a, c, A, alpha, gamma, loss_function):\n",
    "        \n",
    "        self.a = a\n",
    "        self.c = c\n",
    "        self.A = A\n",
    "        self.alpha = alpha\n",
    "        self.gamma = gamma\n",
    "        self.loss_function = loss_function\n",
    "        \n",
    "        self.t = 0\n",
    "        \n",
    "    def step(self, current_estimate):\n",
    "        \n",
    "        a_t = self.a / (self.t + 1 + self.A) ** self.alpha\n",
    "        c_t = self.c / (self.t + 1) ** self.gamma\n",
    "        \n",
    "        delta = np.random.randint(0, 2, np.array(current_estimate).shape) * 2 - 1\n",
    "        \n",
    "        loss_plus = self.loss_function(current_estimate + delta * c_t)\n",
    "        loss_minus = self.loss_function(current_estimate - delta * c_t)\n",
    "        \n",
    "        g_t = (loss_plus - loss_minus) / (2.0 * delta * c_t)\n",
    "        \n",
    "        current_estimate = current_estimate - a_t * g_t\n",
    "        \n",
    "        self.t += 1\n",
    "        \n",
    "        return current_estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "5d12c914",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Variational Quantum Eigensolver\n",
    "def vqe(H, rounds, max_iter, thetas = None, save_steps = None):\n",
    "    \n",
    "    cnt_qubits = int(np.log2(len(H)))\n",
    "    \n",
    "    if thetas is None:\n",
    "        thetas = list(map(\n",
    "            lambda x: random.random(),\n",
    "            [0] * ((1 + rounds) * cnt_qubits * 2)\n",
    "        ))\n",
    "        \n",
    "    optimizer = SPSA(\n",
    "        a = 0.9,\n",
    "        c = 1.0,\n",
    "        A = max_iter / 10,\n",
    "        alpha = 0.6,\n",
    "        gamma = 0.1,\n",
    "        loss_function = lambda t_thetas: expectation(H, rounds, t_thetas)\n",
    "    )\n",
    "    \n",
    "    for i in range(max_iter):\n",
    "\n",
    "        thetas = optimizer.step(thetas)\n",
    "        \n",
    "        if save_steps is not None and i % save_steps == 0:\n",
    "            yield(i, thetas, expectation(H, rounds, thetas))\n",
    "    \n",
    "    return (thetas, expectation(H, rounds, thetas))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "ff6adeec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Step 0. Current expectation value:  2.915\n",
      "Step 3. Current expectation value:  2.349\n",
      "Step 6. Current expectation value:  2.356\n",
      "Step 9. Current expectation value:  1.686\n",
      "Step 12. Current expectation value:  1.457\n",
      "Step 15. Current expectation value:  1.483\n",
      "Step 18. Current expectation value:  1.539\n",
      "Step 21. Current expectation value:  1.396\n",
      "Step 24. Current expectation value:  1.432\n",
      "Step 27. Current expectation value:  1.432\n",
      "Step 30. Current expectation value:  1.411\n",
      "Step 33. Current expectation value:  1.411\n",
      "Step 36. Current expectation value:  1.425\n",
      "Step 39. Current expectation value:  1.421\n",
      "Step 42. Current expectation value:  1.430\n",
      "Step 45. Current expectation value:  1.389\n",
      "Step 48. Current expectation value:  1.408\n",
      "Step 51. Current expectation value:  1.404\n",
      "Step 54. Current expectation value:  1.403\n",
      "Step 57. Current expectation value:  1.395\n",
      "Step 60. Current expectation value:  1.392\n",
      "Step 63. Current expectation value:  1.390\n",
      "Step 66. Current expectation value:  1.394\n",
      "Step 69. Current expectation value:  1.399\n",
      "Step 72. Current expectation value:  1.388\n",
      "Step 75. Current expectation value:  1.381\n",
      "Step 78. Current expectation value:  1.389\n",
      "Step 81. Current expectation value:  1.391\n",
      "Step 84. Current expectation value:  1.380\n",
      "Step 87. Current expectation value:  1.376\n",
      "Step 90. Current expectation value:  1.375\n",
      "Step 93. Current expectation value:  1.383\n",
      "Step 96. Current expectation value:  1.382\n",
      "Step 99. Current expectation value:  1.380\n",
      "The solution is: [\n",
      "\t0.388,\t\t1.034,\t\t1.762,\n",
      "\t1.02,\t\t0.23,\t\t1.609,\n",
      "\t0.123,\t\t0.598,\t]\n",
      "The lowest expectation value is:  1.380\n"
     ]
    }
   ],
   "source": [
    "H = np.array([[3, 1], [1, 2]], dtype = complex)\n",
    "\n",
    "result = None\n",
    "step_size = 3\n",
    "\n",
    "for step in vqe(H, rounds = 3, max_iter = 100, save_steps = step_size):\n",
    "    print('Step {step}. Current expectation value: {ev: .3f}'.format(step = step[0], ev = step[2]))\n",
    "    result = step\n",
    "    \n",
    "print('The solution is: [\\n'\n",
    "     + \"\".join([\"\\t\" + str(round(t, 3))\n",
    "    + \",\\n\" if i % 3 == 2 else \"\\t\" + str(round(t, 3))\n",
    "    + \",\\t\" for i, t in enumerate(result[1])]) + \"]\")\n",
    "                \n",
    "print('The lowest expectation value is: {ev: .3f}'.format(ev = result[2]))            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "68082ed2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAGtCAYAAAAcWBLLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABG1klEQVR4nO3dd5hkVZ3/8c+3qjrnnCbnxMQmDmAPaUBRQAVBzIFddVVcdV3DimL47eKKsKuICBiRYRUEBAFJQ5A8zMDknKcn9HT3dI51fn9U9dgM09M901V9K7xfz1PPdN+6de+3+nTVfPrUueeYc04AAABAsvF5XQAAAADgBYIwAAAAkhJBGAAAAEmJIAwAAICkRBAGAABAUgp4deLi4mI3btw4T87d2tqqrKwsT86NkUVbJw/aOnnQ1smDtk4e0W7rZcuW1TnnSo7c7lkQHjdunF577TVPzr106VLV1NR4cm6MLNo6edDWyYO2Th60dfKIdlub2fajbWdoBAAAAJISQRgAAABJiSAMAACApEQQBgAAQFIiCAMAACApEYQBAACQlAjCAAAASEoEYQAAACQlgjAAAACSEkEYAAAASYkgDAAAgKREEAYAAEBSIggDAAAgKRGEAQAAkJQCXhcwkoJBp80HWtTS5bwuBQAAAB5Lqh7hLXWtOv8nz2rZ/h6vSwEAAIDHkioITyzJUkFmijY1BL0uBQAAAB5LqiBsZlowtlAbGnq9LgUAAAAeS6ogLEknjyvQvjanupZOr0sBAACAh5IuCFePK5AkvbatweNKAAAA4KWkC8KzqvIU8EnLttd7XQoAAAA8lHRBOC3g14Q8n16lRxgAACCpJV0QlqTJ+X6t3nNI7V1cNAcAAJCskjMIF/jU3ev0xq5Gr0sBAACAR5IyCE/K90uSlm1neAQAAECySsognJ1qmlyarVe3ccEcAABAskrKICxJ1eMKtWx7g4JB53UpAAAA8EDyBuGxBWru6NGG/c1elwIAAAAPJG0QPnlcoSQW1gAAAEhWSRuERxdmqCQnTa8xThgAACApJW0QNjOdPK5ArzFzBAAAQFJK2iAsSQvGFmpXQ7v2HurwuhQAAACMsKQOwiePK5Akvbad4REAAADJJqmD8PSKXGWk+LlgDgAAIAkldRBO8fs0b0w+PcIAAABJKKmDsBSaT3jNnia1dPZ4XQoAAABG0KBB2MxGm9nTZrbWzFab2RcH2K/GzFaE93km8qVGR/W4QgWdtGJHo9elAAAAYAQNpUe4R9KXnXPTJZ0m6XNmNqP/DmaWL+kWSe9xzs2UdHmkC42WeWPy5TPpVeYTBgAASCqDBmHnXK1z7vXw182S1kqqOmK3D0q6zzm3I7zf/kgXGi056SmaVp6rZcwnDAAAkFTMOTf0nc3GSXpW0iznXFO/7TdJSpE0U1KOpJudc789yuOvkXSNJJWVlS1YsmTJcGo/YS0tLcrOzj78/e/WdOr53T265dxM+X3mSU2IjiPbGomLtk4etHXyoK2TR7TbetGiRcucc9VHbg8M9QBmli3pXknX9g/B/Y6zQNK5kjIkvWhmLznnNvTfyTl3m6TbJKm6utrV1NQc15OIlKVLl6r/uZsK9ujJu5erbOp8zarK86QmRMeRbY3ERVsnD9o6edDWycOrth7SrBFmlqJQCL7LOXffUXbZJelR51yrc65OoV7jOZErM7qqx4YW1nhlK+OEAQAAksVQZo0wSXdIWuucu3GA3R6QdJaZBcwsU9KpCo0ljguV+RkaX5ylp9fHzdBmAAAADNNQhkYslPRhSSvNbEV42zckjZEk59ytzrm1ZvaopDclBSXd7pxbFYV6o2bxzHLd/twWHWrrVl5mitflAAAAIMoGDcLOueclDXoFmXPuR5J+FImivHDRrHLd+sxmPb52n96/YJTX5QAAACDKkn5luT6zR+WpKj9Dj6ys9boUAAAAjACCcJiZafHMcj23sU7NHd1elwMAAIAoIwj3c9FJ5erqDeqpdVw0BwAAkOgIwv0sGFOg0pw0Pbpqr9elAAAAIMoIwv34fKHhEUvXH1B7V6/X5QAAACCKCMJHuGhWudq7e/XMBoZHAAAAJDKC8BFOGV+ogswU/XUlwyMAAAASGUH4CAG/TxfMKNdT6/ars4fhEQAAAImKIHwUF55UrpbOHj2/sc7rUgAAABAlBOGjWDixWDnpAT3C7BEAAAAJiyB8FKkBn86fXqbH1+xTd2/Q63IAAAAQBQThAVw4q1yH2rv10paDXpcCAACAKCAID+DsKSXKTPUzewQAAECCIggPID3Fr0XTSvX4mr3qDTqvywEAAECEEYSP4aJZ5apr6dKr2+q9LgUAAAARRhA+hkVTS5UW8OlRZo8AAABIOAThY8hKC2jemHyt2NnodSkAAACIMILwIKaU5WjT/hY5xzhhAACAREIQHsSk0my1dPZoX1On16UAAAAgggjCg5hUmi1J2ri/2eNKAAAAEEkE4UFMLs2RJG3c1+JxJQAAAIgkgvAgirNTlZeRok0HCMIAAACJhCA8CDPT5NJsbaJHGAAAIKEQhIdgclk2Y4QBAAASDEF4CCaWZKuhrVsHW5g5AgAAIFEQhIdgcln4grn9DI8AAABIFAThIZh8eAo1gjAAAECiIAgPQUVeurJS/dpMEAYAAEgYBOEhMDNNKuWCOQAAgERCEB6iSaU5LKoBAACQQAjCQzSpNFv7mzt1qL3b61IAAAAQAQThIeq7YG4T44QBAAASAkF4iCaX9QVhxgkDAAAkAoLwEI0qyFRqwEePMAAAQIIgCA+R32eaWJLNXMIAAAAJgiB8HCaXZjNzBAAAQIIgCB+HyaXZ2t3YrrauHq9LAQAAwDARhI/DpPDMEZv3t3pcCQAAAIaLIHwc+maOYIU5AACA+EcQPg5ji7IU8BkzRwAAACQAgvBxSPH7NK44i5kjAAAAEgBB+DhNLs2mRxgAACABEISP0+TSbG0/2KrOnl6vSwEAAMAwEISP08TSbAWdtLWOmSMAAADi2aBB2MxGm9nTZrbWzFab2RePse/JZtZrZu+PbJmxY3JpjiSxsAYAAECcCwxhnx5JX3bOvW5mOZKWmdnjzrk1/XcyM7+k/5L0WBTqjBkTSrLkMzFOGAAAIM4N2iPsnKt1zr0e/rpZ0lpJVUfZ9fOS7pW0P6IVxpj0FL9GF2YShAEAAOKcOeeGvrPZOEnPSprlnGvqt71K0h8knSPpDkkPOef+dJTHXyPpGkkqKytbsGTJkmEVf6JaWlqUnZ19wo+/aVmHDrQH9YMzMyNYFaJhuG2N+EFbJw/aOnnQ1skj2m29aNGiZc656iO3D2VohCTJzLIV6vG9tn8IDrtJ0tecc71mNuAxnHO3SbpNkqqrq11NTc1QTx9RS5cu1XDO/VL7Ot3x/BadedbZCvi53jCWDbetET9o6+RBWycP2jp5eNXWQwrCZpaiUAi+yzl331F2qZa0JByCiyW908x6nHP3R6rQWDKpNFvdvU7b69s0sYS/VAEAAOLRoEHYQun2DklrnXM3Hm0f59z4fvv/WqGhEfdHqMaYM7k0FH437mshCAMAAMSpofQIL5T0YUkrzWxFeNs3JI2RJOfcrdEpLXZNDAfhzQe4YA4AACBeDRqEnXPPSxp44O/b9//YcAqKB9lpAVXmpWvjvmavSwEAAMAJ4kqvEzSxNFub6BEGAACIWwThE1SRl679TZ1elwEAAIATRBA+QYVZaWpo69LxzMMMAACA2EEQPkFFWanq7nVq6ujxuhQAAACcAILwCSrMSpUkNbR2eVwJAAAATgRB+AQVZoeC8EGCMAAAQFwiCJ+gonCPcD1BGAAAIC4RhE9Q4eEgzMwRAAAA8YggfIL6gjBDIwAAAOITQfgEZaYGlJ7iU30LQRgAACAeEYSHoSgrjTHCAAAAcYogPAyFWakMjQAAAIhTBOFhKMxKVUMbQRgAACAeEYSHoSgrVQcZIwwAABCXCMLDUJiVyhhhAACAOEUQHobC7FS1d/eqvavX61IAAABwnAjCw1CY2TeXMItqAAAAxBuC8DAUsswyAABA3CIID0NRNqvLAQAAxCuC8DAUZqVJkhoIwgAAAHGHIDwMDI0AAACIXwThYchNDyjFbwyNAAAAiEME4WEwMxVkpqqeRTUAAADiDkF4mAqzUukRBgAAiEME4WEKrS7HPMIAAADxhiA8TCyzDAAAEJ8IwsNUxNAIAACAuEQQHqbCrDQ1d/SouzfodSkAAAA4DgThYSoMry7HohoAAADxhSA8TEVZLLMMAAAQjwjCw8TqcgAAAPGJIDxM9AgDAADEJ4LwMBX09Qi3MJcwAABAPCEID1NBZqrMGBoBAAAQbwjCw+T3mfIzUlTfRhAGAACIJwThCGB1OQAAgPhDEI6Aoqw0HWwhCAMAAMQTgnAE0CMMAAAQfwjCEVCYTRAGAACINwThCCjKSlVDW5eCQed1KQAAABgignAEFGSmKuikxvZur0sBAADAEBGEI6Aou2+ZZRbVAAAAiBcE4Qgo7FtdrpUeYQAAgHgxaBA2s9Fm9rSZrTWz1Wb2xaPsc7WZvRm+vWBmc6JTbmz6RxCmRxgAACBeBIawT4+kLzvnXjezHEnLzOxx59yafvtslfQO51yDmV0k6TZJp0ah3phUlJUmSTrIzBEAAABxY9Ag7JyrlVQb/rrZzNZKqpK0pt8+L/R7yEuSRkW4zphWkJUiSapnUQ0AAIC4Yc4NfcovMxsn6VlJs5xzTQPs8xVJ05xznzrKfddIukaSysrKFixZsuREah62lpYWZWdnR/SYn3miVWdWBXT19LSIHhfDE422RmyirZMHbZ08aOvkEe22XrRo0TLnXPWR24cyNEKSZGbZku6VdO0xQvAiSZ+UdObR7nfO3abQsAlVV1e7mpqaoZ4+opYuXapIn7v01aeVkZ+vmpp5ET0uhicabY3YRFsnD9o6edDWycOrth5SEDazFIVC8F3OufsG2Ge2pNslXeScOxi5EuMDyywDAADEl6HMGmGS7pC01jl34wD7jJF0n6QPO+c2RLbE+FCYSRAGAACIJ0PpEV4o6cOSVprZivC2b0gaI0nOuVslfVtSkaRbQrlZPUcbh5HICrNStXrPUUeMAAAAIAYNZdaI5yXZIPt8StLbLo5LJoXZoR5h55zCfwwAAAAghrGyXIQUZaWqqzeols4er0sBAADAEBCEI6QwvKgG44QBAADiA0E4QorCyyyzuhwAAEB8IAhHSGE4CLO6HAAAQHwgCEfI4SBMjzAAAEBcIAhHyOEg3EYQBgAAiAcE4QjJTPUrLeCjRxgAACBOEIQjxMxUlJWqg4wRBgAAiAsE4QgKLarR6XUZAAAAGAKCcAQVZqUxNAIAACBOEIQjqCgrlXmEAQAA4gRBOIIKs1LpEQYAAIgTBOEIKsxKVVtXrzq6e70uBQAAAIMgCEcQi2oAAADED4JwBBGEAQAA4gdBOIKKwkGYC+YAAABiH0E4gv7RI8xcwgAAALGOIBxBRVlpksTqcgAAAHGAIBxBuRkBBXzGGGEAAIA4QBCOIDNTAXMJAwAAxAWCcIQVEYQBAADiAkE4wgoyCcIAAADxgCAcYYXZBGEAAIB4QBCOsKKsVB1o7lR7F8ssAwAAxDKCcISdPqFIzZ09etf/Pqc3dzV6XQ4AAAAGQBCOsItOqtBdnzpVbZ29eu8tL+inT21Ub9B5XRYAAACOQBCOgoWTivXYtWfropMq9N9/26ArfvGidhxs87osAAAA9EMQjpK8zBT971XzdPOVc7VhX7MuuvlZ/WnZLq/LAgAAQBhBOMoumVulR689W7Oq8vSVP76hdXubvC4JAAAAIgiPiKr8DP30g/NlJj26aq/X5QAAAEAE4RFTkpOm6rEFemz1Pq9LAQAAgAjCI2rxzHKtrW3SznounAMAAPAaQXgEXTCjXJL02GqGRwAAAHiNIDyCxhRlalp5jv7G8AgAAADPEYRH2OKZ5Xp1e73qWjq9LgUAACCpEYRH2OKZ5XJOemINvcIAAABeIgiPsOkVORpdmME4YQAAAI8RhEeYmWnxjHL9fdNBNXd0e10OAABA0iIIe2DxrHJ19Qa1dP0Br0sBAABIWgRhD8wfU6Di7FSGRwAAAHiIIOwBv8903vQyLV1/QJ09vV6XAwAAkJQIwh5ZPLNcLZ09emHzQa9LAQAASEoEYY+cMalI2WkB/Y3hEQAAAJ4gCHskLeBXzdQSPb5mn3qDzutyAAAAks6gQdjMRpvZ02a21sxWm9kXj7KPmdn/mNkmM3vTzOZHp9zEsnhmuepauvT6jgavSwEAAEg6Q+kR7pH0ZefcdEmnSfqcmc04Yp+LJE0O366R9POIVpmgaqaWKNXv02OrGB4BAAAw0gYNws65Wufc6+GvmyWtlVR1xG6XSPqtC3lJUr6ZVUS82gSTk56ihZOK9NiavXKO4REAAAAjyY4ngJnZOEnPSprlnGvqt/0hSf/pnHs+/P2Tkr7mnHvtiMdfo1CPscrKyhYsWbJk2E/gRLS0tCg7O9uTcx/pmZ3d+tXqLn1vYYZG5zBkO9Jiqa0RXbR18qCtkwdtnTyi3daLFi1a5pyrPnJ7YKgHMLNsSfdKurZ/CO67+ygPeVvCds7dJuk2SaqurnY1NTVDPX1ELV26VF6d+0jle5v0q9XPqWDsNNXMrvS6nIQTS22N6KKtkwdtnTxo6+ThVVsPqQvSzFIUCsF3OefuO8ouuySN7vf9KEl7hl9e4ivJTpMk1TV3elwJAABAchnKrBEm6Q5Ja51zNw6w24OSPhKePeI0SYecc7URrDNhFWSmyu8zHWghCAMAAIykoQyNWCjpw5JWmtmK8LZvSBojSc65WyX9VdI7JW2S1Cbp4xGvNEH5fKbCrFTVNXd5XQoAAEBSGTQIhy+AO9oY4P77OEmfi1RRyaY4O0119AgDAACMKKYpiAElOQRhAACAkUYQjgHF2ak6wMVyAAAAI4ogHANKstNU19LFohoAAAAjiCAcA4qz09TVG1RTR4/XpQAAACQNgnAMKMkJzyXMOGEAAIARQxCOAcXhRTUYJwwAADByCMIxoDgnVRI9wgAAACOJIBwDillmGQAAYMQRhGNA3zLLdS2sLgcAADBSCMIxwB9eZpkxwgAAACOHIBwjWGYZAABgZBGEY0RxdipBGAAAYAQRhGNESU4aY4QBAABGEEE4RpRkp+lAcyfLLAMAAIwQgnCMYJllAACAkUUQjhEsqgEAADCyCMIxoiQ7XRKLagAAAIwUgnCM6OsRPkCPMAAAwIggCMcIllkGAAAYWQThGFGQmSqfiSnUAAAARghBOEb4faYiVpcDAAAYMQThGFIcnksYAAAA0UcQjiEsswwAADByCMIxpCSbZZYBAABGCkE4hpTkpOlAC8ssAwAAjASCcAwpzk5TVw/LLAMAAIwEgnAMYZllAACAkUMQjiEsqgEAADByCMIxpCQnHIS5YA4AACDqCMIxpK9H+EBzh8eVAAAAJD6CcAxhmWUAAICRQxCOIX6fqTCLZZYBAABGAkE4xpTkEIQBAABGAkE4xhRnp+oAs0YAAABEHUE4xrDMMgAAwMggCMeYYpZZBgAAGBEE4RhTEl5mubmTZZYBAACiiSAcY/qWWWacMAAAQHQRhGMMyywDAACMDIJwjDkchLlgDgAAIKoIwjGmJKcvCNMjDAAAEE0E4RjTt8wyY4QBAACiiyAcY1hmGQAAYGQQhGNQcXYqQRgAACDKBg3CZnanme03s1UD3J9nZn8xszfMbLWZfTzyZSaXkpw0HeBiOQAAgKgaSo/wryVdeIz7PydpjXNujqQaST82s9Thl5a8SrLTmD4NAAAgygYNws65ZyXVH2sXSTlmZpKyw/uyLNowsMwyAABA9NlQwpaZjZP0kHNu1lHuy5H0oKRpknIkfcA59/AAx7lG0jWSVFZWtmDJkiUnXvkwtLS0KDs725NzD8UjW7t1z/ou3XJupjJTzOty4lqstzUih7ZOHrR18qCtk0e023rRokXLnHPVR24PRODYiyWtkHSOpImSHjez55xzTUfu6Jy7TdJtklRdXe1qamoicPrjt3TpUnl17qFoyNule9a/oWlzT9aEEt4AhiPW2xqRQ1snD9o6edDWycOrto7ErBEfl3SfC9kkaatCvcM4QX2ryzGXMAAAQPREIgjvkHSuJJlZmaSpkrZE4LhJi2WWAQAAom/QoRFmdrdCs0EUm9kuSddJSpEk59ytkr4n6ddmtlKSSfqac64uahUnAZZZBgAAiL5Bg7Bz7qpB7t8j6YKIVYTDyywThAEAAKKHleViUN8yy4wRBgAAiB6CcIximWUAAIDoIgjHKJZZBgAAiC6CcIximWUAAIDoIgjHKJZZBgAAiC6CcIwqzk5VV09QzZ09XpcCAACQkAjCMerwXMIMjwAAAIgKgnCMYnU5AACA6CIIx6i+IMxcwgAAANEx6Mpy8EZ5brr8PtOPH1+vXuf0rpMq5PeZ12UBAAAkDHqEY1RBVqp+fvV8+c30hbuX64KfPKP7l+9Wb5BZJAAAACKBIBzDLphZrseuPVu3XD1fKX6frr1nhc6/8Rndu2yXenqDXpcHAAAQ1wjCMc7nM73zpAr99Qtn6dYPLVB6il9f/uMbWvTjpbrr5e3q6O71ukQAAIC4RBCOEz6f6cJZ5Xr4C2fqlx+pVmFWmr7551U6+4andftzW9TKfMMAAADHhSAcZ8xM588o0/2fPUN/+NSpmlyWre8/vFYL/+sp3fTEBh1q7/a6RAAAgLjArBFxysx0xqRinTGpWMt3NOiWpZt10xMbtWl/i376wflelwcAABDz6BFOAPPGFOiXH6nWRbPKtaa2yetyAAAA4gJBOIGMKsjQnsZ2OccUawAAAIMhCCeQqvwMdXQHdbCVZZkBAAAGQxBOIFUFmZKk3Q3tHlcCAAAQ+wjCCaQqP0OStLuRIAwAADAYgnACqSoIBeFdDW0eVwIAABD7CMIJJC8jRTlpAYZGAAAADAFBOMFUFWQwNAIAAGAICMIJZlRBhnbRIwwAADAognCCqcqnRxgAAGAoCMIJpqogQ80dPWrq6Pa6FAAAgJhGEE4wVfnMJQwAADAUBOEE848p1AjCAAAAx0IQTjCHF9VgLmEAAIBjIggnmOLsVKUFfFwwBwAAMAiCcIIxM2aOAAAAGAKCcAKqKsjgYjkAAIBBEIQT0ChWlwMAABgUQTgBVeVnqK6lSx3dvV6XAgAAELMIwgmobwo1eoUBAAAGRhBOQH2LajCXMAAAwMAIwgnocI8wQRgAAGBABOEEVJaTJr/PtLuRRTUAAAAGQhBOQAG/T+W56fQIAwAAHANBOEExhRoAAMCxEYQTFItqAAAAHBtBOEGNys/Q3qYOdfcGvS4FAAAgJhGEE1RVQYaCTtp7qMPrUgAAAGLSoEHYzO40s/1mtuoY+9SY2QozW21mz0S2RJwI5hIGAAA4tqH0CP9a0oUD3Wlm+ZJukfQe59xMSZdHpDIMC6vLAQAAHNugQdg596yk+mPs8kFJ9znndoT33x+h2jAMFXnpklhUAwAAYCDmnBt8J7Nxkh5yzs06yn03SUqRNFNSjqSbnXO/HeA410i6RpLKysoWLFmy5IQLH46WlhZlZ2d7cu6R9MWn2zS72K9PnpTmdSmeSZa2Bm2dTGjr5EFbJ49ot/WiRYuWOeeqj9weiMCxA5IWSDpXUoakF83sJefchiN3dM7dJuk2SaqurnY1NTUROP3xW7p0qbw690iasObvCqb6VVNzmteleCZZ2hq0dTKhrZMHbZ08vGrrSMwasUvSo865VudcnaRnJc2JwHExTFX5zCUMAAAwkEgE4QcknWVmATPLlHSqpLUROC6GqaogQ3saOxQMDj78BQAAINkMOjTCzO6WVCOp2Mx2SbpOoTHBcs7d6pxba2aPSnpTUlDS7c65Aadaw8gZlZ+hrt6gDrR0qiw33etyAAAAYsqgQdg5d9UQ9vmRpB9FpCJETN8Uarsa2gnCAAAAR2BluQTWt6gGcwkDAAC8HUE4gR1eVIML5gAAAN6GIJzAstMCystI0e7GNq9LAQAAiDkE4QTHFGoAAABHRxBOcKMKMhgjDAAAcBQE4QRXVZChXQ3tGspS2gAAAMmEIJzgqvIz1NbVq8a2bq9LAQAAiCkE4QQ3qm/mCIZHAAAAvMWgC2ogvvXNJbyroV2zqvKG9Jj2rl49urpWew91qqGtSwdbulTf2qn61i71BJ2+d+kszR9TEM2yAQAAoo4gnOCqjrNHeO+hDn36t69p5e5DkqS0gE9FWakqzE5VYVaa1u9t0uf/sFyPXHuWctNTolY3AABAtBGEE1xBZooyUvxDmkLtzV2N+tRvXlNrZ49+fvV8nT2lRJmpfpnZ4X1e39Ggy299Ud++f5VuunJeNEsHAACIKsYIJzgzU1VBxqCLavzljT26/NYXleL36d7PnqGLTqpQVlrgLSFYkuaPKdAXz52s+1fs0QMrdkezdAAAgKgiCCeBUQUZWrW7SQ+/WavaQ2/tGQ4GnX7y+AZ9/u7lOqkqTw/8y0JNK8895vE+WzNRC8YW6Ft/XqWd9axaBwAA4hNDI5LAudNK9eLmg/rcH16XJFXkpWv+mALNG5Ov5Tsa9fDKWr1v/ij98L2zlBbwD3q8gN+nmz4wVxfd/Jz+9f9WaMk1p8vvs0EfBwAAEEsIwkngw6eP0wdOHqO1tU16fUeDlm1vOByAzaSvXzRN15w94W3DII5ldGGmvnfpTH3pnjf086Wb9C/nTI7iMwAAAIg8gnCSSA34NGd0vuaMztfHF46XJO1r6lBnd1BjijJP6JiXzq3SU+sO6KYnNurMySWaOzr/Lfd3dPdqa12rRhdmKjuNXzUAABBbSCdJrCw3fViPNzN9/9JZen17g65dslzXnjdFG/Y1a+P+Fm3c16wd9W0KOun0CUX6w6dPPa4eZwAAgGjjYjkMS15Gim68Yo6217fp2ntW6LZnt2hrXatmVObq8+dM1qfOHK8XtxzUva8zwwQAAIgt9Ahj2E6dUKSHP3+W/D7T+OIspQb+8fdVMOi0fGejfvDwGp0zrVSFWakeVgoAAPAP9AgjImZU5mpqec5bQrAk+XymH152kpo7evTDv671qDoAAIC3Iwgj6qaW5+iasyfoT8t26cXNB70uBwAAQBJBGCPkC+dO1pjCTH3zzyvV0d3rdTkAAAAEYYyM9BS/vn/pLG2pa9XPl272uhwAAACCMEbO2VNKdMncSv186WZt2t/idTkAACDJEYQxor71rhlKT/Hpm39eKeec1+UAAIAkRhDGiCrJSdPX3zldL2+t1x+X7fK6HAAAkMQIwhhxH6gerflj8vWTxzeopzfodTkAACBJEYQx4nw+0zVnT1TtoQ49tW6/1+UkBeccs3UAAHAEgjA8cd70UpXnput3L233upSkcMNj6zX3+r/pV3/fqmCQsdkAAEgEYXgk4Pfpg6eO0XMb67S1rtXrchLa+r3Nuu3ZLcpOS9F3/7JGV/3yJe042OZ1WQAAeI4gDM9cefJoBXymu+gVjhrnnL79wCrlpAf0ty+drRveP1tr9jTpwpuf1e9f2s7MHQCApEYQhmdKc9O1eFa5/rhsl9q74m/86t5DHfrLG3vU2tnjdSkDevCNPXp5a72+uniqCrNSdUX1aD36pbO1YGyBvnX/Kn34jle0u7Hd6zIBAPBEwOsCkNw+fNpYPfxmrf7y5h5dUT3a63KOy/ceXqOH36xVdlpAl8yt1FWnjNGsqjyvyzqsuaNbP3h4rWaPytOVJ485vL0qP0O//cQp+sMrO/SDh9fqnP9eqjmj8jV3TL7mjg7dKvLSZWYeVg8AQPQRhOGpU8cXanJptu56aXtcBeGWzh49sWafzptepryMFP1p2S7d9fIOzR6Vpw+eMkbvnlOprDRvX143P7FRB1o69cuPVMvve2uoNTNdfepYnT25RHf+fauW72jUr/++TV3h6exKc9J08vhCXXfxDJXmpntRPgAAUUcQhqfMTB86bayue3C13tjZqDmj870uaUj+tnqvOnuC+kzNBC0YW6hvXzxD96/YrT+8vEP/ft9K3fDYej3wuYUaXZjpSX3r9zbrVy9s05Unjz7mz3R0Yaaue/dMSVJnT6/W1jZrxY4GrdjZqEdW7VVnd1C//MgCeocBAAmJMcLw3GXzq5SZ6tfv4+iiuQdW7NGoggzNH1MgScrLTNFHzxinR689S3/859PV0d2rbz+wypOL0Zxz+o/wBXJfXTxtyI9LC/g1d3S+PrZwvG66cp6+fMEUPbF2nx5ZtTeK1QJAfOjpDWrNnqYhz8ne1RPUq9vqtbWulQuTYxg9wvBcbnqKLp1XpXuX7dI33zVd+Zmpb9untbNH2w62amal92Nw61o69fymOv3T2RPe1lNqZjp5XKG+fMFUfe+hNXp4Za0unl05ovU9+MYevbK1Xj+4bJYKs97+sxyqTywcrwdW7NF1D67WwonFystMiWCVABA/DrV361/+8Lqe21injBS/zpxcrPOml2rRtFKV5vxj+FhjW5eWrj+gx9fu0zPrD6glfDF1cXaqFowtUPXYQi0YV6BZlXlK8Zsa27q151C79jR2qPZQu3Y3tmtCcZYuXzBaPh+fxI0EgjBiwodOHas/vLxDf1q2S586a8Lh7c45Pbpqr65/aI1qD3Xo2xfP0CfOHO9hpdLDb9aqN+h0ydyqAff56Olj9eflu/Tdv6zRWZNLlJcxMiGyuaNb3z/KBXInIuD36T/fO1uX/Ox5/eej6/T/3ntShKoEgPix5UCLPvXb17Szvk1fuWCK9jV16sm1+/T4mn2SpDmj8nTahCK9satRr25rUG/QqSQnTe+eU6F3TClVQ1uXXt1Wr2XbG/TY6tBj0gI++czUfkTvcsBn6gk6/WnZLv3n+2ZrYkn2iD/fZEMQRkyYUZmr6rEF+v1L2/WJhePl85m21bXqugdX65kNBzS9IlfTynN0/UNrFHTuLWF5pD2wYremledoannOgPv0hcj3/PR53fDoOv3gsmOHyK11rfKZVJmfoRT/wCOWnHOqb+3Sjvo21R7q0P6mDu1v7tS+pk7tb+7Q9oNtqmvp1O1HuUDuRJw0Kk+fPHO8fvncVl06t1KnTiga9jEBIF48v7FOn71rmQJ+n37/yVMPvwdef8lMra1t1pNr9+mJdft123NbNLUsR595x0SdN6NMs6vy3tKje9UpoY6J/c0dWratQcu2NyjopMr8dFXlZ6gyP0MV+ekqzkrTn5fv1vUPrdFFNz+nL503RZ8+a7wCx/h/YbiCQafapg7trG/T9IrcEeu4iRUEYcSMD58+Vl9cskJPrduvVXsO6Zalm5Xq9+nbF8/QR04fKyfp2iUr9P2H16o36PRP75g44jXuONim13c06msXDj72dlZVnj6+cLzueH6r3ju/SgvGFr5tn2DQ6cbHN+inT2+SJPlMqsjLUFVBhkYXZGpUQYbaunq0/WCbdtS3aWd9m1q73t6DUJKTptKcNE0py9HnFk2M6EWHXzp/ih5ZtVdf//NK/fULZyk9xR+xYwM4cbsa2rRxf4tqppRwQWsU/PbFbfruX9ZoUkm2bv9o9VsufjYzzajM1YzKXH3+3Mnq6gkqNTB4WC3NSddFJ1XoopMqBtznfQtG6awpxbrugdX6r0fX6eGVe3TD++ZoRmXukGt3zulAc6eaO3vU2tmj1s5etXX1qLWrV80d3dpZ366tdS3aVtembQdb1dkTmjGoKCtVX3/ndL1vflXS/E4RhBEzLpxVrqKsVH36d6/JOeni2RX6j4tnqKzf9F03XzlXZtL/e2Sdep3TZ2smjWiNf3lzjyTp3XMGfhPr71/Pn6JHVtbqG/et0kNfOPMtvb0tnT360j0r9Piafbp8wSidMr5QOxvatau+TTsb2vT3TXXa19yhVL9PYwozNaYwU6dNKDr8dWV+hkpz01SYmRrVsWSZqQH98LKT9JE7X9EtT2/Sv14wNWrnAjA0K3Y26hO/flX1rV1695xK/eCyWcpNH15PXkd3rzYfaNHGfS3asK9ZG/a1aNvBVp01uVhfXTxVmanJERm6eoK6/qHV+v1LO3Te9FLddOU8ZQ8yHeZQQvDxKM1J188/tEB/XVmrbz+wSu/56fO6vHq0JpdmqyIvXRX5GarMS1dxdprMpB31bVq5+5BW7j6kVbsPadXuJh1q7x7w+Cl+09iiLI0rytLZU4o1vjhbRdmp+sUzm/WVP76he17doesvmaXpFUMP3/EqOX6rERfSAn59btEk3bd8l7524TSdNbnkbfsE/D7d9IG58vtMNzy6Xr29Tp8/d/KI1Oec0/3Ld+vkcQUaVTC0adGy0gK6/pJZ+tRvX9Ntz27R5xaFgvuOg2361G9f1eYDrfrOu2foo2eMO+pf3109QQV85vlFE2dPKdFl86r082c26+I5lZpSNvCwkMF09vTq+Y11enTVXtW3dumU8YVaOKlYMypyPX+eQDx4at0+fe6u5SrOSdUHTp6o257dohU7G/S/V83X3AE+DXLO6YXNB/WXN/aoqaNb7V296ugOqr27Vx3dvWrp7NGexnYFw5MbBHymccVZKs9N16/+vk1Prt2v/3rfbJ0+MXGHR3V09+qPr+3Urc9s0e7Gdn2mZqK+csHUiAwzO1HvPKlCp08o0vcfXqv7Xt91uOe2T8BnSgv4Dn9SmOI3TS3P0TtPKte08lzlZ6YoKzWgzDS/slIDykrzKystoNKc9KM+r/Onl+lPy3bp/z2yVhf/7/P66Onj9KXzJytnmH9kxTKCMGLKJ84cP+jFcAG/TzdeMVd+M/348Q3q7g3qk2dOiPqsBmtrm7Vxf4u+d+ms43rceTPKdOHMcv3Pkxt18ewK7W5s12fvel3OSb/5+Ck6c3LxgI+NdC/DcHzrXdO1dP1+/fu9b+pP/3zGcYXW1s4ePbPhgB5ZtVdPr9uvls4e5aQFVJKTpifX7Zck5Wem6PQJRTpjUrFqppR4Ngcz4k9Pb1D3Ld+tk8cVanxxltflRNWSV3bom/ev0oyKXN35sZNVkpOm86aX6gt3r9D7f/6CvrJ4qq45a8Lh12dHd68eXLFHd/59q9btbVZOekBluenKSPErI8Wv3IwUleWmKSPFrzHzqjS5LEdTynI0vjjr8PvPy1sO6t/ufVNX/fIlfeT0sfrahdM8XzAokpo7unXXyzt0+3NbVdfSqXlj8vWDy2apZmqp16VJkgqyUvXjK+bovy+frYa2btUealdtY4dqmzpU29iutq5eTS3P0UlVeZpSljOs/zd8PtMVJ4/W+TPK9KO/rdevXtiqh97co+9fOksXzCyP4LOKHYnzm4yk4veZfnT5HPl8pv95apP+56lNh8fITi7L1pSyHE2vyI3o3I0PvLFbAZ/pXccY2zWQ77xnpp6/sU4f/9Wr2l7fpvHFWbr9I9UaF0f/aRdlp+k/Lp6hf/2/N3T+T55R9dhCzRuTr3ljCjSpNPtw70Jv0GlrXYtW72nSmj1NWr2nSa9uq1dnT1CFWam6eHaFLpxVrjMmFis14NPeQx16YXOd/r7poF7YXKdHVu1VwGe6/pJZ+uCpw5v5AolvX1OHPn/3cr2ytV7pKT59/aLp+vBpY0/o04VDbd2665XtWrOnSVX5GRpVkKFR4bH6VQUZng4NcM7ppic26uYnN+odU0p0y9XzD4fRBWML9dcvnKV/v+9N/ecj6/T3TXX6xjun65GVtbrr5R062NqlaeU5uuH9s/WeOZXHPc7/1AlFevSLZ+tHj4WC0dPrQ73DZ0wsVk9vUPuaO7WnsV27G0LTf5XmpOm980d52pM6FI1tXbrz+a369Qvb1NTRo7MmF+uzNfN02oTCmBwfa2YqzEpVYVZq1KcSLchK1Q8vO0lXVI/WN/+8Utf8bpmuOXuC/m3x1KheuOcFGywomNmdki6WtN85N2BXmJmdLOklSR9wzv1psBNXV1e711577TjLjYylS5eqpqbGk3MjsoJBp+c31WlNbZM27GvWxn0t2rS/5fCUNO8cn6Jb/umCiJznzP96StPCvTAn4rcvbtO3H1itc6aV6uYr58blR03OOf3upe16et1+Ld/ZqMa20Bi07LSATqrKU3t3r9btbVJHd+jju1S/T1PKs1U9tlAXzipX9diCY76JOue0ta5V1z+0RkvXH9BHTh+r/7h4xjFn0ujD6zp59LX18xvr9MUly9XW1atvvHOanli7X89sOKAzJxXrhvfPVmV+xpCOt7O+TXf+favueXWn2rp6VZWfoQPNnYeXHO9TlJV6OBT3BeRRBRkqykrT/uZO7Wpo066GUCDc1dimuuYuvfOkCn120UQVZ6ed8PPt6Q3qm39epXte26nLF4zSD9970lFfE8453f3KTn33L6vV2ROUmXTutFJ9YuF4nT6xKCLh7pWt9fq3P72hbQfbVJmXrn3NneoNvj1HzKzM1fWXzNKCsQXDOl80XtfOOd2/Yre+99Ba1bd2afHMMn22ZlLcrGw60jp7evX9h9bqdy9t1ynjC/XTq+aptN+1O0dqaO3S3qYOTSvPOa7fuWi/h5vZMudc9du2DyEIny2pRdJvBwrCZuaX9LikDkl3EoThpWDQaXdju3761Cbd89pO/eQDc3TZvFHDOuYrW+t1xS9e1M1Xzj3m/MHH4pzTqt1NmlGZG/M9JUPhnNO2g21avqNBy3c06o1djcpI8WtmZZ5mVOZqZmWuJpZkn9DHdL1BpxseXadfPLtFp08o0i1Xz1fBIIuD8LqOLuecag91aNP+Fm050KI5o0OfBnjhqaef1hs9VfqfpzZqUkm2brl6viaX5cg5pz+8skPff2itAn7T9ZfM1KVzB776fdXuQ/rFs1v015W1MknvmVOpT589QdMrchUMOtW1dIYuYA0H3MNBt7Fduxra1XXEeE1JSk/xHQ7JqX6fnli7T+kpfn3yzPH69NkThnxBW+2hdr28pV4vb63XC5vrtP1gm75wziR96fwpg4aLDfua9fS6/bpgZnlUhoq0d/XqlqWbtLO+TVUFGarKzwz/m67K/Aw9tW6/vv/QWu1t6tAV1aP0tQunqegE/xCI9Ot6+8FWfev+VXpuY53mjg4NgYiFhZriwf3Ld+vr961UVlpAP/vgvLdMp+mc00tb6nX3Kzv06Kq96uoNalRBhi6eXal3z6nQjIrcQX9vYzYIhx88TtJDxwjC10rqlnRyeD+CMDzX3RvUu3/8mLY0Sf/3T6cPeBHJUHzjzyv159d367VvnZdQY+Ni3b3Ldunr961UWV6abv/Iycecuzmar+u+98lY+rj0UFu3Ont737Kq1fFwzinopKBz6ul1OtTercb2LjW0dutQe5ca2rpV39qlrXWt2ri/RZv3txxeJavP4pll+uriaZpUOnKT/h9o7tRHb31Kaw4G9d75Vfr+pbPeNmRhW12rvvzHN7Rse4MunFmud82uCM+33XH4tvdQh7YdbFNOWkAfPHWMPrZwnCryhtaDLIX+4K5r7dSuhnbVNXeqLDddowoyVJiV+pbfk80HWnTj4xv08Ju1ystI0WdqJuqjp49TRmpoeEJrZ492N/b1IrfrjZ2NennrQe2sb5ck5aQHdMq4Ql06r0rvnjOyq1QOR2tnj/7nqY2647mtykz166uLp+qDp4497k6ASL2uu3uDuv25rbrpiQ1K8fv0bxdO1dUnUE+yW7+3WZ/5/TJtr2/T1y6cqvfNH6X7Xt+tu1/ZoS11rcpJD+h980dpWnmOHl29V89trFNv0GlCSZbePbtS755TOeD7RdwGYTOrkvQHSedIukPHCMJmdo2kaySprKxswZIlS47nOURMS0uLsrNZrSUZ7G1o0X+/6VNPULru9HQVpB9/72RP0OmLT7fppGK//nnOiYUOnLhNjb363+Wd6uxx+qc5aZpXevQ/RCL5uu4NOu1sDmpjQ1AbG3u1oSGorl6n6vKATqsIaFphaFWokeSc064WpzcO9OjNA73a2BCUk1SeaZpW6Ne0Ir+mFfqUn+Y7vP++NqfNjb3a3BjU5kNB7W4JhsPv0M+bn2aqzDZVZPlUle1TZbZPxRmmv+/u0SNbu9XZK501KqDLJqWc0OtrqA60BfX3PT16akeP2nuC+vCMNJ1VFRjwj5Ogc3pka7fu29it3vDzTfGFnk9Buik/zTQhz6+zRwWUmRL9ttze1Kt7N3Trzbpe5aaGzn+wI6jWI2a4yk6Rphb6NaUg1J6jc0b+dy2S9rQE9bs1nVpbH1RhuikrxeQ3KeCT/Cb5fZLfTCZJJpnCNwvdnx/o0fjCNFVl+1SR7VOa/x8/i+6g056WoHY0BbWzOahd4d/vrBTrd5PS/aZndvVoZ3NQC8r8+tD01Kj+ria69h6nO1Z26rV9vTJJTtKkfJ9qRgd0cnngLW3U3OW0bF+PXqrt0fr60HvWD87MUFX223/+0c5mixYtiloQ/qOkHzvnXjKzX4seYcSQpUuXqnzafL33lhc0uSxH91xz2nFfKPLk2n365G9e0x0frda508uiVCmOZe+hDl3zu9f05q5DumhWuf71/CmafMQUbsN9XXf3BnXPqzv1yKpaLd/RqLbwdESVeelaMK5QAZ/pb6v3qrWrV2W5aXr37EpdMrdKs6oG/8hvqLp6gjrU3h2+damxrVsNbd16fUeDlq7brz2HOiRJMypydc60UuVmBPTSlnq9srX+cG/tpPA8oyt3H3rLGO45o/M0vTxX6Sl++SzUu+33mXwWulI8LyNFBZmpys9IUX5mqvIzQ9/39VwezcGWTv306U36/Uvb5TPTxxaO04dOHauKvPSIXFDT0tmjv66s1b3LdunlrfUykxZOLNbislZ9+N3nDOkYtYfa1dTeo/LcdOVmDBycR8orW+t127Ob1Rt0bxtWUJWfqdKctISbRtA5p4ferNWjq/eqqyeont6geoIu9HXQqac3FJBc+BMK50LhqqunV9vrWtUTjilm0uiCTI0vztK+ptAwnZ7wX3XpKT5NLctRWsB/+NONxrbuw9ONleem67uXzNTiBJ35YKT1jUfffrBV750/6pif1vXZ19Shpev364rq0Ud9HcZzj/BWhf6Ak6RiSW2SrnHO3X+sYxKEMRL62vrRVXv1z79fpvfOq9KPr5jzlhdhMOi0fGejnli7TzsOtqmpo1vNHT2H/z3U1q2sNL9e+eZ5Q7poC9HR0d2rny/drDue36rWrh5dNrdK1543RWOKQtOsnejr2jmnR1ft1Q2PrdfWulZNKcvWaROKtGBsgarHFaqq3wVX7V29enLdPj2wYo+Wrt+v7l6nirx0leamKz8jRXkZKcrPDP2bl5Gi3PQU5WYElJMe+jonPaDMVL9qD3Vo28FWba0L3baF/23q6DlqjZmpfp05qVjnTCtVzdRSlee99ZOJnt6gVu9p0otbDurFzQe1v7lTs6vyjjqrRzTsrG/TjY9v0P0rdsu50KwuFXnpb5l1ITstoKaOHjW1d6spHPabOrrV0R1UeopPGakBZab4lZnqV0aqX80dPXp8zT61d/dqfHGW3je/SpfNH6Wq/Azew5PIk089rbGzTtbG8AIfG/c3a8uBVpXmpmlGRa6mV4RWdxtXlHXU3/GO7l4dau9WfmaK0gKsihnLvArCwx7s6Jw7POlrvx7h+4d7XCCSLgz3It74+AZNq8jRR88Ypxc3H9Tf1uzT42v26UBzpwI+05iizHBoCagqP+NwiFk4qZgQ7LH0FL++dP4UffSMcfrFM5v16xe26cE39uiKk0fr8+dMUm/44qbGti7Vt4bGtx5q71J5XoZmVOSqJOftF+u8uq1eP/zrWi3f0agpZdm682PVWjS1dMBew4xUvy6eXamLZ1fqUFu3HllVqxc2H1Rje7ca27q0/WCrGsNBbyjDD8ykyrwMTSjJ0nvmVqosJ115mX2BOvVwoK7MTz/mf+IBv09zRudrzuh8/bMHS4+PLszUTz4wV5+tmahl2xvecmHZ8xtDKyQ6F3q+OWkB5fb7Q6E4O6D27l4dautSbVev2rp61d4d+sj10nmVev+CUZo/psDznlx4w+8zTSrN1qTSbF100vE/Pj3Fz7LwOKZBg7CZ3S2pRlKxme2SdJ2kFElyzt0a1eqACPr8OZO0fm+z/t8j6/S/T25Sc2ePMlP9qplaosUzy1UztVR5GfE3pVmyKcxK1dffOV2fOHO8fvb0Jt39yg4teWVHKHj+7YkBH1eSk6aZlbmaUZGrKWU5enhlrR5fs09luWm64X2z9b4FxzfvaV5miq48ZYyuPOXtcx0Hg07NnT1q7uhWU3v4347Qv62dPSrLTdf44iyNLsxMqP+kJ5flvG3IihSafqmjK6ic9EDCfewPIL4NGoSdc1cN9WDOuY8NqxogisxMP7p8tpycctNTdMHMMp0xsTihgkgyKctN1/WXzNKnz5qg/3ttp3Zs3675M6eoICtVhZmpKsgK9TjuamjXmtomrd5zSGv2NOn5jXXqCTrlpAX01cVT9YmF4485DvZE9I25zctIkbyZYSympAX8fCwNICYxDxSSSmZqQLdcvcDrMhBBowsz9eULpmrp0lrVnDHuqPefPvEf81129vRq8/5WVeanKz/z2HMTAwASG0EYQFJJC/g1ozLX6zIAADGAq38AAACQlAjCAAAASEoEYQAAACQlgjAAAACSEkEYAAAASYkgDAAAgKREEAYAAEBSIggDAAAgKRGEAQAAkJQIwgAAAEhKBGEAAAAkJYIwAAAAkhJBGAAAAEmJIAwAAICkRBAGAABAUjLnnDcnNjsgabsnJ5eKJdV5dG6MLNo6edDWyYO2Th60dfKIdluPdc6VHLnRsyDsJTN7zTlX7XUdiD7aOnnQ1smDtk4etHXy8KqtGRoBAACApEQQBgAAQFJK1iB8m9cFYMTQ1smDtk4etHXyoK2ThydtnZRjhAEAAIBk7REGAABAkiMIAwAAICklVRA2swvNbL2ZbTKzf/e6HkSOmY02s6fNbK2ZrTazL4a3F5rZ42a2Mfxvgde1IjLMzG9my83sofD3tHUCMrN8M/uTma0Lv75Pp60Tk5l9Kfz+vcrM7jazdNo6MZjZnWa238xW9ds2YNua2dfDWW29mS2OZm1JE4TNzC/pZ5IukjRD0lVmNsPbqhBBPZK+7JybLuk0SZ8Lt++/S3rSOTdZ0pPh75EYvihpbb/vaevEdLOkR51z0yTNUajNaesEY2ZVkr4gqdo5N0uSX9KVoq0Txa8lXXjEtqO2bfj/7islzQw/5pZwhouKpAnCkk6RtMk5t8U51yVpiaRLPK4JEeKcq3XOvR7+ulmh/yyrFGrj34R3+42kSz0pEBFlZqMkvUvS7f0209YJxsxyJZ0t6Q5Jcs51OecaRVsnqoCkDDMLSMqUtEe0dUJwzj0rqf6IzQO17SWSljjnOp1zWyVtUijDRUUyBeEqSTv7fb8rvA0JxszGSZon6WVJZc65WikUliWVelgaIucmSf8mKdhvG22deCZIOiDpV+FhMLebWZZo64TjnNst6b8l7ZBUK+mQc+5voq0T2UBtO6J5LZmCsB1lG3PHJRgzy5Z0r6RrnXNNXteDyDOziyXtd84t87oWRF1A0nxJP3fOzZPUKj4aT0jh8aGXSBovqVJSlpl9yNuq4JERzWvJFIR3SRrd7/tRCn3sggRhZikKheC7nHP3hTfvM7OK8P0VkvZ7VR8iZqGk95jZNoWGOJ1jZr8XbZ2Idkna5Zx7Ofz9nxQKxrR14jlP0lbn3AHnXLek+ySdIdo6kQ3UtiOa15IpCL8qabKZjTezVIUGYj/ocU2IEDMzhcYRrnXO3djvrgclfTT89UclPTDStSGynHNfd86Ncs6NU+h1/JRz7kOirROOc26vpJ1mNjW86VxJa0RbJ6Idkk4zs8zw+/m5Cl3rQVsnroHa9kFJV5pZmpmNlzRZ0ivRKiKpVpYzs3cqNLbQL+lO59wPvK0IkWJmZ0p6TtJK/WPc6DcUGif8f5LGKPRGe7lz7sgB+4hTZlYj6SvOuYvNrEi0dcIxs7kKXRSZKmmLpI8r1IlDWycYM/uupA8oNAvQckmfkpQt2jrumdndkmokFUvaJ+k6SfdrgLY1s29K+oRCvwvXOuceiVptyRSEAQAAgD7JNDQCAAAAOIwgDAAAgKREEAYAAEBSIggDAAAgKRGEAQAAkJQIwgAQR8xsm5l9xes6BmNm48zMmVm117UAwEAIwgDikpn9Ohy0jry95HVt/YXrfOgEHvcdM1t1lLtOlnTL8CsDAAS8LgAAhuEJSR8+YluXF4WMFOfcAa9rAIBEQY8wgHjW6Zzbe8Stb2Wid5hZd3j1OYW3/bOZNZnZhPD3S83sVjO72cwawrcfmZmv32NSzey/zGyXmbWa2atmtrh/EWY2zcweNLNDZtZiZi+a2Ulm9h2Flg59V78e65rwY/7TzNabWXt4uMMNZpYevu9jCq28NLPf4z4Wvu8tQyPMbIyZ/dnMmsO3+8xsVL/7v2Nmq8zsSjPbHN7nfjMrHuiHGq7/x0dsyw3Xeln4+w+FfxbNZrbfzP5oZlXHOGZN+HkU99v2tuETZjbDzB7ud9y7zax8oOMCwHAQhAEkJOfcM5J+JOl3ZlZoZtMk/VjS551zW/rterVC74WnS/onSddIurbf/b+S9A5JH5R0kqTfSPqLmc2RJDOrlPS8JCfpfEnzJf1MoaXc/1uhJUSfkFQRvr0QPm6rQkuITpf0WUlXSvpm+L57wrWu7/e4e458jmZmCi1TWibpHEmLJFVKuj98X59xCi1de5mkCyTNk3SsJeZ/L+nK/n8QSHqfpHZJD4e/T1UorM+RdLFCS6fefYxjDsrMKiQ9K2mVpFMknafQErsPHlELAEQEQyMAxLMLzazliG0/c859Lfz1dQqF09sVCoMPOed+c8T+tZK+4ELrza8zsymS/lXSjWY2UdJVksY553aE9/+pmZ2nUGj+rKTPKRRqL3fO9Q3L2NB3cDNrV7jnuv9JnXPf6/ftNjP7oaSvSPoP51x7+Hn1HPm4I5ynUBCd6JzbFj7fByVtknSuQgFcCr3Xf8w5dyi8z22SPn6M4y6R9BOFgvWT4W1XS/pj33N0zt3Zb/8tZvYZSWvNbJRzbtcxjn0sn5H0Rr/2k5l9RFK9pGpJr5zgcQHgqAjCAOLZswr14PbX2PeFc647HAxXS9qvUK/pkV4Kh+A+L0r6npnlKtS7a5LWvLWDVWmSngp/PU/S8/1C8JCY2fsV6nmepFCvpz98Ox7TJe3pC8GS5JzbYmZ7JM3QP4Lw9r4QHLZHUulAB3XOHTSzxxQKv0+Ge2oXSfpuv/rnK/SHxlxJhQr9nCRpjKQTDcILJJ19lD9uJGmiCMIAIowgDCCetTnnNg2yz2kKDX3Il1SifkF5CHwKDXk4WVL3Efe1h/81HSczO02hXtfvSvpSuKb3KDSU4rgOFa7vaPpvP7J2p8GHxv1e0m1m9lmFesV3KjQERGaWJekx/eNixf0KDY14TqEhE0cT7Fdzn5Qj9vEpNPTiaNPD7RukXgA4bgRhAAnLzMZJ+qlCwxculHSXmZ3hnOvpt9upZmb9eoVPU6iXtcnMlisU3Mqdc08PcJrXJX3IzFIH6BXu0tt7ehdK2t1/eISZjR3C4460RlKVmY3rNzRigkLjhNcM8tjBPCDpNoXG/14t6a5+P6NpCgXfbzjntobP+95Bjtc320VFv6/nHrHP65KuUKgH+8jwDgARx8UHAOJZmpmVH3ErkSQz8yvUq/mMc+4Xkj4laZRCH+f3VynpJjObGh6u8FWFxsfKObdB0l2Sfm1m7zezCWZWbWZf6Rf8blFoaMP/mdnJZjbJzK4ys7nh+7dJmhU+frGZpSg0hrjKzK4OH/MzCvW69rdN0lgzmx9+XNpRnv8Tkt5QKOAvCM++cJdCgfKpo+w/ZM65Dkn3SfqWQkNEft/v7h2SOiX9S7j+d0n63tuP8habFOpV/o6ZTTGzC8LH7u9nkvIk3WNmp4aPfZ6Z3WZmOcN5PgBwNARhAPHsPIUudut/Wx6+7xsKjb/9pBQa96rQVGb/bmZn9jvGXQr1vL4s6ZeS7lA4CId9XKGZI26QtE7SQ5LOlrQ9fNzd4e9TJT0dPv/nJfX1Ov9S0lpJrynUE7rQOfcXhWa0uEnSmwpd0PftI57bvZL+qtDFagf09qCscA/tpeH7l4bPv1fSpUeMez5Rv1PoYrzXnXNr+533gEI/y0sV6nm+TqELDAcU7uG9UtIEhcL7dxVqo/777FGotzwo6VGFxnb/TKHQ3RmB5wMAb2GRea8EgPhjZkslrXLO/YvXtQAARh49wgAAAEhKBGEAAAAkJYZGAAAAICnRIwwAAICkRBAGAABAUiIIAwAAICkRhAEAAJCUCMIAAABISv8fhH6O3+2KsAwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "steps = []\n",
    "evs = []\n",
    "\n",
    "for step in vqe(H, rounds = 3, max_iter = 100, save_steps = 1):\n",
    "    steps = steps + [step[2]]\n",
    "    evs = evs + [step[0]]\n",
    "    \n",
    "plt.figure(figsize = (12, 7))\n",
    "plt.plot(evs, steps)\n",
    "\n",
    "plt.xlabel(\"Number of executions\", fontsize = 14)\n",
    "plt.xlabel(\"Expectation value\", fontsize = 14)\n",
    "plt.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "3a1c1fa3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues:  [3.61803399+0.j 1.38196601+0.j]\n"
     ]
    }
   ],
   "source": [
    "eigenvalues, eigenvectors = scipy.linalg.eig(H)\n",
    "\n",
    "print(\"Eigenvalues: \", eigenvalues)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4af0310",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}