deuteron_vqe.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell is added by sphinx-gallery\n",
    "# It can be customized to whatever you like\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Accelerating VQEs with quantum natural gradient\n",
    "===============================================\n",
    "\n",
    ".. meta::\n",
    "    :property=\"og:description\": Accelerating variational quantum eigensolvers\n",
    "        using quantum natural gradients in PennyLane.\n",
    "    :property=\"og:image\": https://pennylane.ai/qml/_images/qng_example.png\n",
    "\n",
    ".. related::\n",
    "\n",
    "   tutorial_vqe Variational quantum eigensolver\n",
    "   tutorial_quantum_natural_gradient Quantum natural gradient\n",
    "\n",
    "*Authors: Maggie Li, Lana Bozanic, Sukin Sim (ssim@g.harvard.edu). Last updated: 8 Apr 2021.*\n",
    "\n",
    "This tutorial showcases how one can apply quantum natural gradients (QNG) [#stokes2019]_ [#yamamoto2019]_\n",
    "to accelerate the optimization step of the Variational Quantum Eigensolver (VQE) algorithm [#peruzzo2014]_.\n",
    "We will implement two small examples: estimating the ground state energy of (1) a single-qubit VQE\n",
    "problem, which we can visualize using the Bloch sphere, and (2) the hydrogen molecule.\n",
    "\n",
    "Before going through this tutorial, we recommend that readers refer to the\n",
    ":doc:`QNG tutorial </demos/tutorial_quantum_natural_gradient>` and\n",
    ":doc:`VQE tutorial </demos/tutorial_vqe>` for overviews\n",
    "of quantum natural gradient and the variational quantum eigensolver algorithm, respectively.\n",
    "Let's get started!\n",
    "\n",
    "\n",
    "(1) Single-qubit VQE example\n",
    "----------------------------\n",
    "\n",
    "The first step is to import the required libraries and packages:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from pennylane import numpy as np\n",
    "import pennylane as qml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this simple example, we consider the following single-qubit Hamiltonian: $\\sigma_x + \\sigma_z$.\n",
    "\n",
    "We define the device:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [],
   "source": [
    "dev = qml.device(\"default.qubit\", wires=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the variational ansatz, we use two single-qubit rotations, which the user may recognize\n",
    "from a previous :doc:`tutorial </demos/tutorial_qubit_rotation>` on qubit rotations.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [],
   "source": [
    "def circuit_4(params,wires):\n",
    "    t0 = params[0]\n",
    "    t1 = params[1]\n",
    "    t2 = params[2]\n",
    "    qml.PauliX(wires=0)\n",
    "    qml.RY(t2, wires=1)\n",
    "    qml.RY(t1, wires=2)\n",
    "    qml.RY(t0, wires=3)\n",
    "    qml.CNOT(wires=[3,0])\n",
    "    qml.CNOT(wires=[0,2])\n",
    "    qml.RY(-t1, wires=2)\n",
    "    qml.CNOT(wires=[0,2])\n",
    "    qml.CNOT(wires=[2,0])\n",
    "    #qml.CNOT(wires=[1,0])\n",
    "    qml.CNOT(wires=[0,1])\n",
    "    qml.RY(-t1, wires=2)\n",
    "    qml.CNOT(wires=[0,1])\n",
    "    qml.CNOT(wires=[1,0])\n",
    "def circuit_3(params,wires):\n",
    "    t0 = params[0]\n",
    "    t1 = params[1]\n",
    "    qml.PauliX(wires=0)\n",
    "    qml.RY(t1, wires=1)\n",
    "    qml.RY(t0, wires=2)\n",
    "    qml.CNOT(wires=[2,0])\n",
    "    qml.CNOT(wires=[0,1])\n",
    "    qml.RY(-t1, wires=1)\n",
    "    qml.CNOT(wires=[0,1])\n",
    "    qml.CNOT(wires=[1,0])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then define our cost function using the ``ExpvalCost`` class, which supports the computation of\n",
    "block-diagonal or diagonal approximations to the Fubini-Study metric tensor [#stokes2019]_. This tensor is a\n",
    "crucial component for optimizing with quantum natural gradients.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coeffs = [1, 1]\n",
    "obs = [qml.PauliX(0), qml.PauliZ(0)]\n",
    "coeffs = [15.531709 ,-2.1433, -2.1433, 0.21829, -6.125, -9.625, -3.91, -3.91]\n",
    "obs = [qml.Identity(0),qml.PauliX(0) @ qml.PauliX(1), qml.PauliY(0) @ qml.PauliY(1),qml.PauliZ(0),qml.PauliZ(1),qml.PauliZ(2),qml.PauliX(1) @ qml.PauliX(2),qml.PauliY(1) @ qml.PauliY(2)]\n",
    "\n",
    "H = qml.Hamiltonian(coeffs, obs)\n",
    "cost_fn = qml.ExpvalCost(circuit, H, dev)\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kron(i,j):\n",
    "    if(i==j):\n",
    "        return 1\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "def matrix_element(i,j):\n",
    "    \n",
    "    ele = 0.0\n",
    "    \n",
    "    ele = ((2.0*i + 1.5)*kron(i,j) - np.sqrt(i*(i+0.5))*kron(i,j+1) \\\n",
    "    - np.sqrt((i+1)*(i+1.5))*kron(i,j-1)) * 3.5\n",
    "    \n",
    "    ele += -5.68658111 * kron(i,0) * kron(i,j)\n",
    "    \n",
    "    \n",
    "    return ele\n",
    "\n",
    "def pauli_token_to_operator(token):\n",
    "    \"\"\"\n",
    "    DO NOT MODIFY anything in this function! It is used to judge your solution.\n",
    "\n",
    "    Helper function to turn strings into qml operators.\n",
    "\n",
    "    Args:\n",
    "        token (str): A Pauli operator input in string form.\n",
    "\n",
    "    Returns:\n",
    "        A qml.Operator instance of the Pauli.\n",
    "    \"\"\"\n",
    "    qubit_terms = []\n",
    "\n",
    "    for term in range(len(token)):\n",
    "        # Special case of identity\n",
    "        if token[term] == \"I\":\n",
    "            pass\n",
    "        else:\n",
    "            #pauli, qubit_idx = term, term\n",
    "            if token[term] == \"X\":\n",
    "                qubit_terms.append(qml.PauliX(int(term)))\n",
    "            elif token[term] == \"Y\":\n",
    "                qubit_terms.append(qml.PauliY(int(term)))\n",
    "            elif token[term] == \"Z\":\n",
    "                qubit_terms.append(qml.PauliZ(int(term)))\n",
    "            else:\n",
    "                print(\"Invalid input.\")\n",
    "    if(qubit_terms==[]):\n",
    "            qubit_terms.append(qml.Identity(0))\n",
    "    full_term = qubit_terms[0]\n",
    "    for term in qubit_terms[1:]:\n",
    "        full_term = full_term @ term\n",
    "\n",
    "    return full_term\n",
    "\n",
    "\n",
    "def parse_hamiltonian_input(input_data):\n",
    "    \"\"\"\n",
    "    DO NOT MODIFY anything in this function! It is used to judge your solution.\n",
    "\n",
    "    Turns the contents of the input file into a Hamiltonian.\n",
    "\n",
    "    Args:\n",
    "        filename(str): Name of the input file that contains the Hamiltonian.\n",
    "\n",
    "    Returns:\n",
    "        qml.Hamiltonian object of the Hamiltonian specified in the file.\n",
    "    \"\"\"\n",
    "    # Get the input\n",
    "    coeffs = []\n",
    "    pauli_terms = []\n",
    "    chunks = input_data.split(\"\\n\")\n",
    "    # Go through line by line and build up the Hamiltonian\n",
    "    for line in chunks:\n",
    "        #line = line.strip()\n",
    "        tokens = line.split(\" \")\n",
    "        # Parse coefficients\n",
    "        sign, value = tokens[0][0], tokens[1]\n",
    "\n",
    "        coeff = float(value)\n",
    "        if sign == \"-\":\n",
    "            coeff *= -1\n",
    "        coeffs.append(coeff)\n",
    "\n",
    "        # Parse Pauli component\n",
    "        pauli = tokens[3]\n",
    "        \n",
    "        pauli_terms.append(pauli_token_to_operator(pauli))\n",
    "\n",
    "    return qml.Hamiltonian(coeffs, pauli_terms)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  (-9.625) [Z0]\n",
      "+ (-6.125) [Z1]\n",
      "+ (0.21829055499999983) [Z2]\n",
      "+ (15.531709445) [I0]\n",
      "+ (-3.913118960624632) [Y0 Y1]\n",
      "+ (-3.913118960624632) [X0 X1]\n",
      "+ (-2.1433035249352805) [Y1 Y2]\n",
      "+ (-2.1433035249352805) [X1 X2]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/kyle/.local/lib/python3.10/site-packages/pennylane/vqe/vqe.py:153: UserWarning: ExpvalCost is deprecated, use qml.expval() instead. For optimizing Hamiltonian measurements with measuring commuting terms in parallel, use the grouping_type keyword in qml.Hamiltonian.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "from qiskit_nature.operators.second_quantization import FermionicOp\n",
    "from qiskit_nature.mappers.second_quantization import JordanWignerMapper, ParityMapper\n",
    "\n",
    "mapper = JordanWignerMapper()\n",
    "#mapper = ParityMapper()\n",
    "\n",
    "N = 3\n",
    "op=0\n",
    "for i in range(N):\n",
    "    for j in range(N):\n",
    "        op += matrix_element(i,j) * \\\n",
    "        FermionicOp([([(\"+\", i),(\"-\", j)], 1.0)])#,register_length=2)\n",
    "\n",
    "    #op = 30.2*FermionicOp([(\"+\",0),(\"-\",4)],register_length=2)\n",
    "#print(op)\n",
    "\n",
    "hamstr = \"+ \"+str(mapper.map(second_q_op=op))\n",
    "\n",
    "#qml.load(mapper.map(second_q_op=op),format='qiskit')\n",
    "\n",
    "hamiltonian = parse_hamiltonian_input(hamstr)\n",
    "\n",
    "print(hamiltonian)\n",
    "\n",
    "cost_fn = qml.ExpvalCost(circuit_3, hamiltonian, dev)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To analyze the performance of quantum natural gradient on VQE calculations,\n",
    "we set up and execute optimizations using the ``GradientDescentOptimizer`` (which does not\n",
    "utilize quantum gradients) and the ``QNGOptimizer`` that uses the block-diagonal approximation\n",
    "to the metric tensor.\n",
    "\n",
    "To perform a fair comparison, we fix the initial parameters for the two optimizers.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [],
   "source": [
    "#init_params = np.array([3.97507603, 3.00854038])\n",
    "init_params = np.array([2.5,4.5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will carry out each optimization over a maximum of 500 steps. As was done in the VQE\n",
    "tutorial, we aim to reach a convergence tolerance of around $10^{-6}$.\n",
    "We use a step size of 0.01.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_iterations = 500\n",
    "conv_tol = 1e-06\n",
    "step_size = 0.01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we carry out the VQE optimization using the standard gradient descent method.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration = 0,  Energy = 2.62754488 MeV,  Convergence parameter = 0.41528526 MeV\n",
      "Iteration = 20,  Energy = -1.45759287 MeV,  Convergence parameter = 0.07290178 MeV\n",
      "Iteration = 40,  Energy = -1.98314671 MeV,  Convergence parameter = 0.00676702 MeV\n",
      "Iteration = 60,  Energy = -2.03668718 MeV,  Convergence parameter = 0.00089606 MeV\n",
      "Iteration = 80,  Energy = -2.04435261 MeV,  Convergence parameter = 0.00013419 MeV\n",
      "Iteration = 100,  Energy = -2.04548435 MeV,  Convergence parameter = 0.00001928 MeV\n",
      "Iteration = 120,  Energy = -2.04564498 MeV,  Convergence parameter = 0.00000269 MeV\n",
      "\n",
      "Final value of the energy = -2.04566218 MeV\n",
      "Number of iterations =  131\n"
     ]
    }
   ],
   "source": [
    "opt = qml.GradientDescentOptimizer(stepsize=step_size)\n",
    "\n",
    "params = init_params\n",
    "\n",
    "gd_param_history = [params]\n",
    "gd_cost_history = []\n",
    "\n",
    "for n in range(max_iterations):\n",
    "\n",
    "    # Take step\n",
    "    params, prev_energy = opt.step_and_cost(cost_fn, params)\n",
    "    gd_param_history.append(params)\n",
    "    gd_cost_history.append(prev_energy)\n",
    "\n",
    "    energy = cost_fn(params)\n",
    "\n",
    "    # Calculate difference between new and old energies\n",
    "    conv = np.abs(energy - prev_energy)\n",
    "\n",
    "    if n % 20 == 0:\n",
    "        print(\n",
    "            \"Iteration = {:},  Energy = {:.8f} MeV,  Convergence parameter = {\"\n",
    "            \":.8f} MeV\".format(n, energy, conv)\n",
    "        )\n",
    "\n",
    "    if conv <= conv_tol:\n",
    "        break\n",
    "\n",
    "print()\n",
    "print(\"Final value of the energy = {:.8f} MeV\".format(energy))\n",
    "print(\"Number of iterations = \", n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then repeat the process for the optimizer employing quantum natural gradients:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration = 0,  Energy = 1.49924499 MeV,  Convergence parameter = 1.54358516 MeV\n",
      "Iteration = 20,  Energy = -2.04299636 MeV,  Convergence parameter = 0.00076858 MeV\n",
      "Iteration = 40,  Energy = -2.04565644 MeV,  Convergence parameter = 0.00000439 MeV\n",
      "\n",
      "Final value of the energy = -2.04566796 MeV\n",
      "Number of iterations =  46\n"
     ]
    }
   ],
   "source": [
    "opt = qml.QNGOptimizer(stepsize=step_size)\n",
    "\n",
    "params = init_params\n",
    "\n",
    "qngd_param_history = [params]\n",
    "qngd_cost_history = []\n",
    "\n",
    "for n in range(max_iterations):\n",
    "\n",
    "    # Take step\n",
    "    params, prev_energy = opt.step_and_cost(cost_fn, params)\n",
    "    qngd_param_history.append(params)\n",
    "    qngd_cost_history.append(prev_energy)\n",
    "\n",
    "    # Compute energy\n",
    "    energy = cost_fn(params)\n",
    "\n",
    "    # Calculate difference between new and old energies\n",
    "    conv = np.abs(energy - prev_energy)\n",
    "\n",
    "    if n % 20 == 0:\n",
    "        print(\n",
    "            \"Iteration = {:},  Energy = {:.8f} MeV,  Convergence parameter = {\"\n",
    "            \":.8f} MeV\".format(n, energy, conv)\n",
    "        )\n",
    "\n",
    "    if conv <= conv_tol:\n",
    "        break\n",
    "\n",
    "print()\n",
    "print(\"Final value of the energy = {:.8f} MeV\".format(energy))\n",
    "print(\"Number of iterations = \", n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualizing the results\n",
    "^^^^^^^^^^^^^^^^^^^^^^^\n",
    "\n",
    "For single-qubit examples, we can visualize the optimization process in several ways.\n",
    "\n",
    "For example, we can track the energy history:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAFYCAYAAAB+s6Q9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABHlklEQVR4nO3deZzN9eLH8df3rLMag0FClqwhpCxdUSqVW9269VOSirQvkpC0KYU21U1lKaSuoroqaREhawvFWMqaVBjrbGbO8v39cZwxY2bMdpY5zvvZYx4zZ/t+P/Nxmvf5LN/PxzBN00REREQqPUu4CyAiIiKlo9AWERGJEAptERGRCKHQFhERiRAKbRERkQih0BYREYkQtnAXoCR796YH9HjJyXEcOJAV0GNGG9VhxakOK051GBiqx4oLdB2mpCQW+1jUtbRtNmu4ixDxVIcVpzqsONVhYKgeKy6UdRh1oS0iIhKpFNoiIiIRQqEtIiISIRTaIiIiEUKhLSIiEiEU2iIiIhFCoS0iIhIhKv3iKiIi4fDHHzt55ZUXOXBgP2ByyimnMnjwMKpWrVrmY/3nP+Np1Kgx1apV56+//uSqq64p0+u/+24RHTt2wW63F/l4VlYW/fr1ZvbsT8tctopYuHA+559/YUjPGe3U0hYROY7H42HkyGH06XMjkyZNY9Kk6TRr1pxXXnmhwPO8Xm+ZjtupU5cyBzbAzJnv4nK5yvy6YHK5XLz//nvhLkbUUUtbROQ4q1atoEGDhrRt2z7vvuuvvxHTNAEYPfoJbDYbBw8eZOTIJ3jyyZFkZWWRk5PDAw88RMuWrfjyy895773ppKTUxOFw0qhRYz7//FO2bt3CPfcM4sMPP+Cbb77CNL2cf/6F/N//9WHy5DfIyspkx44d7Nq1k0GDHuLgwQOsX7+OIUPu4+WXX89rbWdmZjBy5DBycnI488x2eeX8+ec1TJz4GhaLhVq1ajN8+KOkpaUxatRIDMPA7Xbz6KOjqFWrNmPHPs3+/Xtxu708+ugokpOr8cILY9i16w9crlwGDryL9u07cPfdAznnnE788MMqDh8+zLhxLzFjxjS2bNnM88+PYciQ4aH9B4piIQ3t7Oxshg8fzr59+8jKyuLuu++mR48eITv/tAWrqJacwOXtWobsnCJSMU884eTTTwv/qbJYwOuNL9cxL7/czRNP5BT7+M6dO2jc+PTjzlewYzIpqSrDho1k587f6dnzMnr0uJiffvqBd9+dxtNPj2PixAlMmTKDhIQEBgy4scBr//77LxYvXshrr00C4M47B9C9ew8Mw2DPnt288MIrLF++lDlzPuLZZ59n8uQ3eP75Vwp0j3/55TwaN27CPfcMYsGC+Xz99RcAvPzy84wfP4EqVarw2msvs3DhfPbtS+Pssztyyy0D2bAhlX370li79meSkqry8ssv8t//zmbJkkXEx8dTrVp1hg0bycGDB7n//juZNu2/GIZBfHw8r776JhMmvMKiRQvp0+fGox8mFNihFNLQXrBgAa1atWLgwIHs2rWL/v37hzS0H1l1F6bXwuXtVoXsnCISedxuDx6PJ+/28OGDycjIYO/ePUybNhOAFi18H/6TkpL47rvFfPTRLHJycoiJieHQoUPExsbljX+3bn1mgeNv2rSR33/fwb333g5AVlYmf/31JwBt2rQFoFatWmRkFL9h0vbtW2nb9iwA2rXz9QgcPnyYnTt3MGLEEACOHDlCcnIy55zTmREjhpCRkUG3budz5pntWLhwfl4L/cILewIwfvzz/PjjKn7+eTUAOTk5ed3y/ufWqlWLQ4cOlak+JXBCGtq9evXK+/nvv/+mVq1aoTw9FruLI65csrIgLi6kpxaRcnriiZwiW8UpKYns3ZsZlHM2atSYOXM+yrs9ZsyLAFxzzeV549g2m6/VO2vWTFJSUnj88adJTV3H66+/gmmaGEb+I5qFztGpUxeGDRtZ4L4ffliF1Xps8wl/d3xRTJO8c+R/XvXqNfjPfyYWev7Uqf9l5cplvPrqi/TqdQUAXm/h49944y1cfPGlhe4vbbkkuMIyEe3aa69lyJAhPProoyE9r8NqA4uLbds0/05Eitehwzn8+eculi37Lu++det+ISsrC6u14N+Pw4cPceqp9QDfbGq3201SUhKZmZmkp6fjdrv55Zc1BV7TrFkLfvrpB44cOYJpmowf/zw5OUeKLY9hWAq0/AHq1z+NTZs2AvDjj98DUKVKFQzDYOvWzQDMnj2TzZt/Y/78L/n99+2cf/6FXH99PzZsWE+zZi1Ys+YnAJYuXcK0aVNo2bIVS5YsAuDAgf28+eZrZSqTBF9YJqLNmjWL1NRUBg8ezJw5cwqNFeWXnBwXsG3P4mIcHPa4SEuLJyUlIIeMWifa71VKR3VYccGsw7ffnsLo0aN59923cblcVK9enUmTJlK3bgoxMXaSkmJJSUmkT5/eDB06lOXLF3PDDTewbNlili9fyP3338f9999O3bp1adGiOfHxDiwWC3FxDlq3bsKAAf0ZNOgODMPgwgsvpG7dFOLjnSQkxJCSksiBA/E4HDZSUhLp3LkjgwbdwdSpU6lWrRoAffv25u6772bIkHvo0KEDhuGrjzFjnmXs2Gew2WzUrl2bAQNuokoVJ6NGjcLhcOB2u3n88cdp1KgRI0eupm/fvhiGwbhx40hJSWH9+jXce+9A3G439957LykpiTgcNpKT40lJSSQhIQaXy0nz5g0wTQ+jRz/K+PHjg/bvEClC9f+zYYawn2Pt2rVUr16dOnXqAL7u8unTp1O9evViX7N3b/FjOmV1zpSubD+8hRFGGoMG5QbsuNHG1y0ZuH+XaKQ6rDjVYWCoHisu0HV4og8AIe0nXr16NdOmTQMgLS2NzMxMkpOTQ3b+uBg7WF1s2aLucRERiTwhTa/rrruOtLQ0+vTpwx133MHjjz9+wq7xQIuP8Y1pK7RFRCQShXRM2+Fw8MILL5T8xGCd32oHi5et2zTzUUREIk9UNTltFt9nlP0HPRw4EObCiIiIlFFUhbbdcnQ1IauLrVuj6lcXEZGTQFQll816NLQ1ri0iIhEoqpJLLW0RKa1du/5g6NAHuPXWftx6az+ef/5ZsrOzA3b8NWt+OrrtZ+gtXDi/wse45prLycrKCkBpClq6dAmjRz/Bvn1pjBs3usyv37z5N37/fccJnzNgwI15y8aGynffLQrITm1RlVz2o2PaWBTaIlI8r9fLyJFD6d27D5MnT2fy5OnUqVOXZ58dFbBzzJ37SdhCe8aMaaV6Xlm3Hi1JWZYFqV69BkOHPlLmcyxatICdO38v8+uCLVDbq0bV1py2oy1tZ1yuusdFpFjff7+S+vVP46yzzs677/rr+3L99f9m//59vP76q3Tv3oNzz+3K0qVL+Pbbb3jkkSeYOXMGCxd+g9frpXPnc+nf/zY++2wO69b9woED+9mxYwd9+95ErVq1WLLkW7Zt28rTT49jwIAbmDv3GwBGjhzK1Vf/H6ZpMnv2+1itVjZt2sitt97Ot98uYOvWzQwbNpL27TvklW3Xrj949tlR1KlzKr/9tonmzc9g2LBHjm6d+Qw2mx2LxcJTT43hs8/msHnzr4wY8RDXXNObuXM/5tFHfS3aXr16MHfuN6XaerQoX301jxkzplK79ilUqZJE27btsVgsLF++lH370njiidF88MF/SU39hdxcF//617+5/PJ/sWXLZp5++jGqVEmiTp26APz115+MHDmMKVPeKXK70Z9/Xs3HH8/CMCxs27aVHj0u4rzzzmfOnI9YtGgBycnJBco5fvxzrFu3lvr1T8Pt9oVnWtpexo59mtzcXKxWK8OGjaR69RqMGvUoaWl7OXIkm1tuGUjXrt35739nsHjxQtxuN3feeS/t23fI217VajXo2vX8Mm+vWh5RFdr+7vG6p+WwZZ2lwIL7IlI5PbFsJJ9u+V+h+y0Wo8gNL0rj8sb/4okuTxf7+I4d2zn99KYF7jMMg0aNGrN9+7ZiX2ea8Morr+NwOPm///sXvXv3wWKxsGXLZt544y127vydJ598hLfffo/TT2/K4MFDqV27dpHHMgyDzZt/4733ZrNo0QImT36DGTNm8e233/Dll58XCG2LxcLGjet58slnSEqqyrXXXkF6ejoHDuzn7rsfoFWr1kyZ8iZffTWPPn368e6703jmmef46acfiv1dTrT16OjRzxXxu5u88cZ/mDJlBnFxcdx003V5Zdy7dw+vvz7l6HKwNZgwYQq5uTn07n0Vl1/+L6ZOncytt97Jued25YUXxuaFql9R241Wr16D9etTee+9D/F43PTufRW33DKQjh070717jwKBvW3bVlJT1zJp0jTS0vZy7bW+DVMmT36D3r1voEOHc1i27DumTZvClVf+m4MHDzBhwmQOHTrIsmXf8eefu1i4cD5vvPEWf/65i3feeZs6dU7N2141JSWRa675vzJvr1oeURXa/pZ23fo5bFllsGePQa1aumZbRApyuXKL3AzDNL1YrcX/2bTZbDzwwD1YLBYOHjzA4cOHAWjVqg1Wq5WaNWuRkZFR6nKcfvrp2O12qlWrTsOGjbDZbFSrVp3MzMK7m9WtW5/q1WsAHH1OBklJSbz++qvk5uaSlraXiy66pNTnPtHWo0U5dOgQcXHxeatc+rcYBWjevAWGYeBwODh8+BB33XUrVqs1b3hg+/ateSHbtm17VqxYmvfa4rYbrV69Bs2aNc8rj2kW35W/fftWmjc/A8MwSEmpSd26vg1eNmxYz44d25k6dTJer5fk5GTq1z+NjIx0nnrqUbp27c5FF13C0qWLad68BRaLhbp16/Hww4+xaNHCvO1VHQ5bubZXLY+oCm3/mPap9X3rjm/ZYqFWLe1SI1KZPdHl6SJbxcFcM/u00xoyd+6cAvd5vV62bdtG/fqnYeTrovOH+65dfzB79kzefvtd4uLiueGGa/KeU5ZtLfN/WMj/uvw/F7XVZ8HHfed5+eUXuOGGfnTu/A9mzJhKTk7BLU6N47oa85/7RFuPFuX47UjzH9t/rJ9++oE1a37k1VcnYrFYuPji846+Nv/zC/9uRW03+tNPPxT6nYtzfK9q/n+DUaOeJSWlZoHnT5o0nZ9/Xs2cOR+xcOF8LrjgoiJ7dfzbq+Z/L5Zle9XyiKqBXX9Lu05d3xtX49oiUpSzz+7Ili2bC2zNOXPmu5xxxhkkJycTH5/A4cOHANi0aQMA6enpJCdXIy4unnXr1rJnz54TTjyyWI5tbWm12sjNzcXtdrN58+aA/R7+bUNzcnJYtuw73G43cGwf7YSEBA4cXWlq164/iuwFKGrr0aIkJSWRnp7O4cOHyck5krft5/HHqlXrFGw2G4sWLcDrNXG5XEe3GfXVo3+bUb/ithstjmEYRW5j+uuvmzBNk927/+aPP3YCFNiK9Mcfv+frr79g06aNfPPN17Rv34FBgx5iw4b1NGnSjPXr1+LxeNi/fx8PP/xgQLZXLY8oa2n7QrtWnWMtbRGR4zmdTsaOfYnnn3+WSZMm4PV6admyVd5s5p49L2Ps2KdYs+Yn4uPj8Xq9NGnSlPj4BO6661ZatWrN1Vdfy/jxz3HhhT2LPEfbtu157LGHGT36Oa688moeeOBuGjRoSMOGDY+2Wis+4ebaa6/nkUce4tRT63L99X0ZP/55evS4mKZNm3H77bfw+utTcDgcPPLIQ9StW5/atU8pNGP80ksv5+mnH2PRogVcffX/sWTJt8yb91mhc1ksFm6+eQB33XUrDRo0pEWLM7BarQWC6uyzO/Lee+9w33130LnzP+jatRsvvjiOm24awDPPPMns2TOpU+fUQj0Cw4c/xrPPPoXNZqNmzZpcccXVrFv3S5G/85lntuPVV18kLi6ODh3OAeD005vQqFFjbr/9FurXP42mTZvj9XoZMOA2Ro9+gm+++QrDMBgx4nESE6swceIEPvnkI3Jzc7n11juoU+dULrroUu6661a8Xi+33XYXtWvXpnfvG7jnntuw26107twVp7PooQOAdu3ac889t/Hyy69TtWrV0v0DFiGkW3OWRyC7v8asfIoXf3yO6efPo1+3S7jkEhfTpxf/yUiKpq38Kk51WHGhqsPU1HWMH/8cb775dkg3OAqVQNbjwoXzOeuss6lSJYnBg++hf//baNWqTUCOXZmdtFtzhpu/ezw2wUVSkqlrtUWkRGec0Yo2bc6kf/++AVmU5GR25MgR7rvvTu64oz9169aLisAOtajsHveYLho39rJ2rQWPB0o5l0FEotS99w4OdxEiwqWX/pNLL/1nuItxUouqpqa/pe3yumnUyIvLZfD777pQW0REIkNUhbb/ki+Xx0WTJr7JFps3R1UViIhIBIuqxPLv8uX2ujj9dF9o//ZbVFWBiIhEsKhKLHte97iLpk0V2iIiElmiKrH8oe32umnY0IvVaiq0RUQkYkRVYtn8Y9peFw4HNGhg8ttvVir3leoiIiI+URXa+bvHAZo08XDggMG+fZpBLiIilV9UhbbNcmwiGpA3g1xd5CIiEgmiKq38l3zlKrRFRCQCRVVa5bW0Pb7Q1mVfIiISSaIqrQqPaSu0RUQkckRVWuW/5AsgKQlq1vRqVTQREYkIUZVWdmvBljZA06Zedu40yMoKV6lERERKJ7pC+7jZ4+Ab1zZNgy1boqoqREQkAkVVUtksRbe0QePaIiJS+UVVUtnzbc3ppxnkIiISKaIqqfzLmLqLaGlrMpqIiFR2UZVUx1/yBXDKKSZxcSa//hpVVSEiIhEoqpKqqIlohuG7XnvrVgseT7hKJiIiUrKoCm1bEWPa4AvtnByD33/XxiEiIlJ5RVVo24sY0watjCYiIpEhqlKqqEu+4NhktE2brCEvk4iISGlFVWgXdckXQPPmvsHsTZuiqjpERCTCRFVKFXXJF0CDBiZOp6nQFhGRSi2qUsowDGwWGy5PwdC2Wn2LrPz6qwWvN0yFExERKUFUhTb4usiPb2kDNGvmJTvbYMcOzSAXEZHKKfpC22ovNKYN0KKFfzJa1FWJiIhEiKhLqOJb2v7JaJpBLiIilVP0hbbVXuiSL/B1jwNs3Bh1VSIiIhEi6hLK19Iu3D1+2mkmsbGmQltERCotWzhO+uKLL7Jy5UpcLhcDBw7k0ksvDdm57VY72blHCt1vsfhWRtu0ybcGuVW95CIiUsmEvFn5/fffs2HDBt5//33eeustnn322ZCe324punscoHlz3xrk27drBrmIiFQ+IQ/tdu3aMX78eAASExNxuVx4Q3hxtN1a9EQ0yD+urWa2iIhUPiEPbZvNRnx8PACzZ8+mW7duWCyhK4avpV14TBu0nKmIiFRuYRnTBpg/fz4ffPABb7/99gmfl5wch80WuJavv6WdkpJY6LEuXXzft21zkpLiDNg5T0ZF1Z+Ujeqw4lSHgaF6rLhQ1WFYQnvJkiVMmDCBKVOmUKVKlRM+98CBrICe2z+mvXdveqHHYmMhLi6BX37xsndvYM97MklJSSyy/qT0VIcVpzoMDNVjxQW6Dk/0ASDk/cDp6emMGTOGiRMnkpycHOrTY7fa8ZpevGbhcXSLxTcZbfNmC66ih71FRETCJuQt7c8//5xDhw7xwAMP5N03duxY6tSpE5Lz2/Ptqe20Fu4Cb9bMy08/Wdm2zZK3z7aIiEhlEPLQ7t27N7179w71afPYrUdD25NbTGh7ADsbNyq0RUSkcom6adL5W9pFad5cy5mKiEjlFHXJlNfSLvayL19ob9gQdVUjIiKVXNQlk7+lXdwCK6ecYlK1qsmGDVpgRUREKpeoC22H1QEU3z1uGHDGGR62bTPIzAxlyURERE4s6kK7pJY2QMuWXkzT0Li2iIhUKlGXSiWNaYMvtAHWr1cXuYiIVB7RF9olzB4HX/c4wPr1UVc9IiJSiUVdKvlb2m5P8aHdtKkXi8UkNTXqqkdERCqxqEulYy3t4rvH4+KgUSMv69dbMc1QlUxEROTEoi+0rSVPRAPfuPbhwwa7dhmhKJaIiEiJoi+0SzGmDXDGGf7JaFFXRSIiUklFXSKVvqXtm4yWmqoZ5CIiUjlEX2iXYkwb8l/2FXVVJCIilVTUJdKx67RP3NKuW9ekShVToS0iIpVG1CVSaVZEA99ypi1betiyxUJ2dihKJiIicmLRF9qlbGmDr4vc6zX49deoqyYREamEoi6NjrW0TzymDcfGtbXIioiIVAZRl0ZlaWkfW85UM8hFRCT8oi+0SzmmDdCsmRfD0GQ0ERGpHKIujcrS0k5IgIYNTdat03KmIiISftEX2qW8TtuvdWsPBw8a7Nyp5UxFRCS8oi+0S7HLV36tW/smo61dq3FtEREJr6gLbYfVAUCuN7dUz2/TxjcZbe3aqKsqERGpZKIuiZxWJwC5ntKFtr+l/csvammLiEh4RV1ol7WlXb26yamnetXSFhGRsIu6JHLa/C3tnFK/pnVrD7t3W9i9W5PRREQkfKIutP0tbVcpJ6LBsS7ydeuirrpERKQSiboU8o9p55Shpe2fjKZxbRERCaeoC+28lnYpx7Qh/2S0qKsuERGpRKIuhfxj2jmlnD0OcMopJjVqeHWttoiIhFWpQjstLY3Vq1cD4HaXbiWxysrfPe4qQ2gbBrRq5eX33y0cPBikgomIiJSgxND++OOPueGGG3jssccAGD16NO+9917QCxYs/u7xnDJ0j8Oxce1169TaFhGR8CgxtN9//33mzJlD1apVARg+fDgffvhhsMsVNDaLDQOjTC1tgDZtNK4tIiLhVWICxcbGEhMTk3fb6XTidDqDWqhgMgwDh9VRpuu0AVq10gxyEREJL1tJT0hMTGTOnDnk5OSQmprKvHnzqF69eijKFjQOq5PcUmzNmV+DBiaJiaau1RYRkbApMYFGjRrFzz//TG5uLiNHjuTIkSM89dRToShb0Dgs9jK3tC0WX2v7t98sZGYGqWAiIiInUGJLu2rVqnmT0E4WDquz1BuG5NemjZfly22sXWulUydPEEomIiJSvBJDu1u3bhhGwTW3vV4vixcvDlqhgs1udZQrtNu18wX16tUWhbaIiIRciaGd//Iuj8fDL7/8wl9//RXUQgWb0+Igy5VR5tf5Q3vNGitQtjFxERGRiioxtE899dQCt+vXrx/x3eW+lnbZQ7dBA5PkZJPVqzWDXEREQq/E0F6+fHmB2/v372ft2rVBK1AoOMtxyRf4VkY780wP335r48ABSE4OQuFERESKUWJoT5gwIe9nwzBISkpi2LBhQS1UsNktDnLLuCKaX7t2vtBes8bK+edrXFtEREKnxNB+5513QlGOkHJanXhNL26vG5ulxCoooG1b38poq1crtEVEJLSKTaw+ffoUmjWe37vvvhuUAoWC3WoHINeTW+bQbt/ePxlNi6yIiEhoFZtYgwYNKvZFJwrzSOA4utNXrieHOHtcmV5bq5bJKad4NRlNRERCrtjm4jnnnJP3dcYZZ1C3bl3q1q1LzZo1eeGFFyp00l9//ZULL7yQGTNmVOg45eWw+Hb6KutSpn5t23rYvdvCX39F9ocXERGJLCX2DU+aNIk333yT3NxcYmJicLlcXHHFFeU+YVZWFk899RSdO3cu9zEqyr89Z3lmkAO0a+dl3jz46ScrvXpF9v7iIiISOUocmP3yyy9ZtmwZZ555JqtWreLZZ5+lYcOG5T6hw+Fg0qRJ1KxZs9zHqKhjLe3yzSBv21bj2iIiEnql2prT4XDg8fiC6pJLLmHJkiXlPqHNZiuw1Wc4HGtpVyy0Na4tIiKhVGL3ePXq1Zk1axannXYaDz30EM2aNePvv/8ORdkASE6Ow2YLbDhWTUgEIKGKnZSUxDK/PiUFmjSBX36xUaNGIhE+L69cylNvUpDqsOJUh4Gheqy4UNVhiaE9duxY0tLSuPTSS5k6dSp///03L774YijKBsCBA1kBPV5KSiLuow3sv9P2s9eWXq7jtG4dw0cf2Vm1KoNGjcwAlrDyS0lJZO/e8tWb+KgOK051GBiqx4oLdB2e6ANAid3jDz/8MJs3byYuLo577rmHkSNH0qxZs4AVLhz83eOucs4eh2Obh/z4o7rIRUQkNEoM7e7du/Pf//6XCy64gKeffrrC646vW7eOG2+8kY8//pjp06dz4403cvDgwQods6z8E9Fyyjl7HOCssxTaIiISWiV2j19xxRVcccUVpKen8/XXXzNhwgR27tzJZ599Vq4TtmrVKuxLo/oXV3GVcyIaQOvWXpxOk++/V2iLiEholOqaJdM0Wb9+PWvXrmXbtm00b9482OUKKofFt4xpTgVC2+mENm28rF9vITMzUCUTEREpXokt7ccee4xFixbRokULevXqxdChQ4mNjQ1F2YImbxlTb/m7xwE6dPDw/fdW1qyxcu652jxERESCq8TQbtasGYMGDaJatWqhKE9IVPQ6bb+zz/bw+uvw/fcKbRERCb4SQ/uGG24IRTlCKpChDfDDDxrXFhGR4IvKdTj9s8dd5VzG1K9WLZN69bz88IMFM7ou1RYRkTCIztA+OqZdkYlofmef7WH/fgtbt0bhsmgiIhJSJXaPb9q0iQ8//JBDhw5h5mtOjhs3LqgFCyaH1Td7vCKXfPl16ODho4/sfP+9lcaNteOXiIgET4mhPXjwYC677DJatGgRivKExLGWdsVmj4MvtME3rn3ddQptEREJnhJDu1q1atx9992hKEvI2C0VX8bU74wzvMTGapEVEREJvhLHtC+44AIWL15Mbm4uXq837yuSOa0VX8bUz273bdW5caOFdK25LyIiQVRiS/u1114jIyOjwH2GYbBhw4agFSrY7AGaPe7XoYOH5ctt/PSTlW7ddL22iIgER4mhvWrVKiyWk2uSudO/IloAJqIBdOjg63n4/nuFtoiIBE+JoZ2VlcXUqVNZt24dhmHQrl07+vXrR0xMTCjKFxT2o7PHAxXa/kVWVq7UuLaIiARPiU3oRx55hIyMDK6//np69+7N7t27GTlyZCjKFjTHWtoVH9MGqFHDpGlT3zrkbk0gFxGRICmxpb1//35efvnlvNvdu3enb9++QS1UsOUtYxqA2eN+nTp5mD7dytq1Ftq1i+yJeiIiUjmV2NLOzs4mKysr73ZmZiY5OYFpoYaLfxnTQLW0wRfaAMuXq4tcRESCo8SW9nXXXcell15KmzZtME2T1NRU7rvvvlCULWgcAZ6IBtC5sy+0V6ywctddgWvBi4iI+JUY2tdccw3nnnsuqampADz66KPUqlUr6AULJrvl6ES0AF3yBXDqqSb163tZscKG1wsn2YR7ERGpBIqNlkWLFgEwe/Zsli5dysGDBzl48CBLlixh9uzZIStgMBiGgcPiCGhLG3xd5AcPGmzapMQWEZHAK7alvWnTJrp168aPP/5Y5OPXXHNN0AoVCnZr4EO7c2cPH3xgZ/lyKy1aaDKaiIgEVrGhfdtttwHwj3/8g169ehV4bObMmcEtVQg4rY6ATkQD6NTJd73XihVW+vfXuLaIiARWsaG9YcMG1q1bx1tvvUV2dnbe/W63mwkTJnDdddeFpIDBYrc4AjqmDdCokUlKipfly62YJhjaYltERAKo2NB2OBzs27eP9PT0Al3khmEwZMiQkBQumJxWZ8C7xw3D10X+ySd2tm0zaNTILPlFIiIipVRsaDdu3JjGjRvTqVMnGjduTGJiIgBpaWnUqFEjZAUMFrvVTnZOdslPLCN/aK9caaVRIy2PJiIigVPiNOc1a9YwbNiwvNsPPPAA06dPD2qhQsFhcQa8exzyL7JS4tV0IiIiZVJiaM+dO5dXXnkl7/bbb7/NZ599FtRChYLD6sAV4O5xgBYtvCQlmSxbppXRREQksEoMbZvNhs12rNVonCSzqxxWBzkBnj0OvkVVunRx8/vvFnbsODnqSkREKocS+3C7detG7969ad++PaZpsnLlSi655JJQlC2oHBYHHtODx+vBaglsq/i88zzMm2dnyRIbp52mS79ERCQwSgztO+64g86dO7NmzRoAnnjiCc4888xglyvoju30lUusJTagx+7a1TeuvWSJlb59FdoiIhIYpVpv0+l00qxZM5o2bUpWVhbLly8PdrmCzh/awRjXbtLES61aXpYs8V2vLSIiEggltrTvvvtuNm/eTM2aNfPuMwyDzp07B7Vgwebf6SsnCKFtGPCPf3j48EM7GzdatKSpiIgERImh/eeff/LFF1+cNBPQ/IKxp3Z+553n5sMP7SxZonXIRUQkMErsHm/UqBG5uYFvjYZbjC0GgBzPkaAc/9i4tq7XFhGRwCgxUaxWK7169aJNmzYFLv0aN25cUAsWbP7QznYHJ7Tr1jVp2NDLsmVW3G6wKbtFRKSCSoySc845J+LHr4sSY/XNGM92ZwXtHF27upk+3cGaNRY6dFAXuYiIVEyJod2xY8dC93m9kR9A/pb2kSB1j4Pveu3p031d5B06nHxDDCIiElolhvZNN92EYRiYponH4yEtLY3GjRvzySefhKJ8QRNr87W0j7gDv2mI37nnHrte+4EHgnYaERGJEiWG9oIFCwrc/uuvv3jvvfeCVqBQibEGd0wboHp1k1atPHz/vZXsbIgN7BouIiISZUq1uEp+p5xyCrt37w5GWUIqJgQtbYBu3Tzk5BgsX64NREREpGJKbGmPHz++wDXa+/btY9u2bUEtVCiEYkwb4IIL3Lz2moMFC2xccIEnqOcSEZGTW6l2+bJarXlfTZo04bXXXgtF2YIqFGPaAB07eoiLM/nmG13zJSIiFVNskowaNYrHHnuM/fv389hjj4WyTCHhD+1gjmkDOBy+1dG++MLO9u0GDRpoMXIRESmfYkP7u+++48EHH2TlypVkZGQUejzyF1cJTUsb4PzzPXzxhZ0FC2z0769dv0REpHyKDe0pU6bw448/sn79+pN0cZXQjGmDb1wbYOFChbaIiJRfsaFdr1496tWrR/v27alfv37ATjh+/HhWrFhBbm4uTz75JK1btw7YscsilC3t004zOf10D0uWWMnJAacz6KcUEZGTUIkT0QIZ2CtWrGDdunXMnDmTMWPGMGbMmIAdu6xi/bPHgzym7dejh4esLIOVK3Xpl4iIlE+Zr9OuiJUrV9KjRw8AmjZtyp49e8jODn5LtyjH1h4PzfnPP9/XRa5Z5CIiUl4lhvbcuXML3Tdz5sxynWzv3r1Uq1Yt73a1atVIS0sr17EqKq97PARj2gBduniIjTVZuFAtbRERKZ9im33r168nNTWVt956q0Br2OVy8frrr3PdddeV+WR2u73AbdM0CyzcUpTk5DhstsAGXUpKInG5vs8rXouLlJTEgB6/ON27w7x5Vo4cSaRevZCcMmhCVWcnM9VhxakOA0P1WHGhqsNiQ9vpdLJv3z7S09P58ccf8+43DIMhQ4aU62QpKSns27cv7/b+/fupUaPGCV9z4EBgt85MSUlk7950vKZvp7LDWRns3Zse0HMUp2tXO/PmxTBz5hFuvjlyZ5H761DKT3VYcarDwFA9Vlyg6/BEHwCKDe3GjRvTuHFjOnXqRNu2bQFfyzgzM5OEhIRyFeS8887jpZdeok+fPqSmplKvXj1iYmLKdayKshgWHBYHRzyhG1O/+GI3I0bAl1/aIjq0RUQkPEoc0964cSMzZswgNzeXq666iq5duzJ58uRynaxVq1Y0b96cq666iscff5zhw4eX6ziBEmOLDfqKaPnVr2/SsqXv0q8i1qsRERE5oRJD++OPP+a6665j/vz5tGjRgpUrV/L111+X+4QPPfQQH3/8MbNnz6ZZs2blPk4gxNhiQnKddn6XXOImN9dg4ULNIhcRkbIpMbSdTic2m40lS5Zw2WWX4XA4cJ4kq4PE2GJDdp2236WX+i79mjdPoS0iImVTquu0n3jiCVasWEHHjh1Zs2YNLtfJMR4ba40J6Zg2QJs2XurU8TJ/vg23O6SnFhGRCFdiaD///PM0aNCASZMm4XA42LVrF0888UQIihZ84WhpGwb07Onm4EGtjiYiImVTYmjXrFmTli1bsnjxYqZOnUqdOnXCPhYdKDG2GLLd2ZhmaLfLvOQSXxP7iy/URS4iIqVXYmi/9NJLPPfcc6SlpbF7926eeuop3nzzzVCULehirDGYmOR6c0N63i5dPCQkmMybZyPEnxdERCSCldjUW7lyJe+//z4Wiy/fXS4Xffv25fbbbw964YItNt9OX05r6CbXOZ3Qo4ebOXPsbNhgoWVLb8jOLSIikatUE9H8gQ2+pUhLWno0UhwL7dCOa4O6yEVEpOxKTIxWrVpx++23849//AOApUuXhm0P7EDzbxqS5Q7sUqmlceGFbux2k08/tTF4cGi750VEJDKVGNojRoxg3rx5/PzzzwBcccUVXHrppUEvWCjEhHhP7fySkuD88z189ZWNzZsNTj9dg9siInJiJwztnTt3Uq9ePXr16kWvXr3Izs5m9+7dJ033uH9P7VCviuZ3xRUuvvrKxief2NXaFhGREhU7pr18+XKuv/560tOP7Vzy+++/M2DAANauXRuSwgVbrL+lHaI9tY93ySVuHA6TOXM0ri0iIiUrNrRfffVV3nrrLRITj20R1qxZM9544w1eeumlkBQu2Pxj2tlhamlXqQIXXOBmwwYrv/1WqjmBIiISxYpNCovFQtOmTQvd36RJk5NmGdNwjmn7XXGFbxb5J5+otS0iIidWbGhnZRU/ozrjJNlXMm9MO8Trj+fXs6cbp9NUaIuISImKDe1GjRrxwQcfFLp/ypQpNG/ePKiFCpXK0NJOTDzWRb5pk7rIRUSkeMU27x5++GHuvvtuPvroI1q1aoXX62XNmjXYbDYmTZoUyjIGTWyYx7T9rrzSzbx5dj75xMZDD2kWuYiIFK3Y0K5evTozZ85k+fLl/Pbbb3g8Hi666CI6duxYYIW0SBYTxhXR8rv4YjcxMSYff2xjyJBcTpIr6kREJMBKHEjt3LkznTt3DkVZQi7G6r/kK7wt7YQE3+Vf//ufnTVrLLRrp7XIRUSksJOjyVxOlaWlDfB//+ebkf/BB/Ywl0RERCqrqA7tvMVVwjymDdC9u4caNbx8/LGNXA1ri4hIEaI6tP2XfGVXgpa2zQb//reb/fstfPONLv8SEZHCoju0bZVjTNvvWBe5QltERAqL8tCuPGPaAK1aeWnRwrfz14ED4S6NiIhUNlEd2pVpTBvAMHytbZfL4H//04Q0EREpKMpDOw4I/+Iq+f37324sFlOzyEVEpJCoDm2bxUaMNYYMV3rJTw6R2rVNunXz8OOPWtZUREQKivpUSHAkkJFbuTZA6dvXNyHtnXfU2hYRkWMU2vZEMlyVK7QvucRNSoqX99+3k115eu5FRCTMFNqORNJzK0/3OIDdDjfc4OLQIUNbdoqISJ6oD+1ERyKZrgy8ZuVa77tvXxeGYTJtmiPcRRERkUoi6kM7wZ6AiUmWKzPcRSmgfn2T88/38MMPVtavj/p/JhERQaFNoiMRoNKNawP06+ebkDZ9uiakiYiIQpt4uy+0K9u4Nvj22a5d28usWXYyK1dHgIiIhEHUh3ZeS7sShrbN5puQlp5u8OGHam2LiES7qA/tBHsCAOmVaIGV/G66yYXdbjJxoh1v5ZorJyIiIabQdvhCu7ItsOJXu7bJv/7l5tdfrSxcaA13cUREJIyiPrQT7VUASM89HOaSFO+OO3IBeOMNXf4lIhLNoj6081ralXD2uF/r1l7OPdfNokU2Xf4lIhLFoj4BKvNEtPxuv93X2p44URPSRESiVdSHtv+Sr8rc0ga4+GIPDRt6mT3bzp49RriLIyIiYRD1oX1scZXK3dK2WOC223LJzTV46y21tkVEolHUh3beJV+VvHsc4LrrXFSv7mXyZAeHK++8ORERCZKoD+1jY9qVu3scID4e7rjDxeHDBlOmaCa5iEi0ifrQTvAvY1rJu8f9+vfPpWpVkzfftJNR+T9niIhIAEV9aNutdpxWJ5kR0D0OkJgIAwfmsn+/halTNbYtIhJNQh7aq1atonPnzixcuDDUpy5WoiMxIsa0/QYOzCUhwWTCBAdZWeEujYiIhEpIQ/v333/n7bffpn379qE8bYni7QmV/pKv/KpWhVtvzSUtzcKMGWpti4hEi5CGdkpKCv/5z39ITEwM5WlLlOioElEtbYDbb3cRF2fyyisObdspIhIlbKE8WWxsbJlfk5wch80W2I0yUlIKfmhIjksiNS2D6jXisRiRMcyfkgIPPghPPWUwY0YiI0eG+vyV64NXJFIdVpzqMDBUjxUXqjoMWmjPmjWLWbNmFbjv3nvvpWvXrmU6zoEDgR20TUlJZO/egq3qGCMOE5Mdf/5NgiNy3rw33wwTJsQzdqzBv/+dSY0aZkjOW1QdStmoDitOdRgYqseKC3QdnugDQNBC+9prr+Xaa68N1uEDyr/ASoYrI6JCOzERHnwwlxEjYhg/3sHTT+eEu0giIhJEkdEXHGQJDv/2nJH3abNfPxennebl7bft7NihNclFRE5mIQ3tb7/9lhtvvJElS5bw4osv0r9//1Cevlh5Le0IDG2HA0aMyMHlMhgzxhnu4oiISBCFdCJa9+7d6d69eyhPWSr+PbUjZVW04115pZvXXvPw4Yd2BgzIpUMHb7iLJCIiQaDucSJr/fGiWCzkjWcPHx6DxxPmAomISFAotIEqjiQADuceCnNJyq9TJw/XXuvil1+sTJ+uBVdERE5GCm2gWkx1APZl7wtzSSrmscdySEw0efZZJ2lpmpQmInKyUWgDNWJrAJCWvTfMJamYWrVMhg3L4eBBg9GjtXWniMjJRqENVD8a2vuOpIW5JBXXv7+LFi08vPuug1Wr9M8rInIy0V91jrW092VHfmjbbDBuXA6GYXL//bFkZ4e7RCIiEigKbSDBnojT6oz47nG/jh09DBzoYssWC2PH6tptEZGThUIbMAyD6jE1In4iWn4jRuTQoIGXN96w8/33+mcWETkZ6K/5UdVja5B2EnSP+8XFwSuvHME04f77Y9RNLiJyElBoH1UjtgZZ7kyyXIHdVSycOnXydZNv3mxl9Gh1k4uIRDqF9lEn0wzy/EaMyKFpUw8TJzr44ovA7ksuIiKhpdA+qvpJNIM8v7g4mDjxCDExJvfdF8sff2jRFRGRSKXQPiolNgWI/AVWitKypZfRo32LrtxxRwxud7hLJCIi5aHQPqp6jH9VtJOrpe3Xt6+Lf/3LxapVNp59VquliYhEIoX2Uce6x0+ey77yMwx44YUjNGzo5dVXnXz0UUh3ZRURkQBQaB91sqw/fiKJifDOO9kkJpoMGhTD6tX65xcRiST6q33UyTp7/HhNm3qZODGb3Fzo1y+Wv/7SxDQRkUih0D7qZFp/vCQ9enh4/PEcdu+20K9fLBkZ4S6RiIiUhkL7KP/649EQ2gB33OHihhty+flnKzffHEtOTrhLJCIiJVFoH+Vff/xknT1+PMOA557L4ZJLXCxebOOuu2LweMJdKhERORGFdj6+9cf3YppmuIsSEjYbvPnmEbp0cfPpp3aGDnUSJb+6iEhEUmjnc2rCqWS5s9h/ZH+4ixIysbEwfXo2rVt7eOcdB8OHO/F6w10qEREpikI7n9OSGgKw/fDWMJcktKpUgfffz6ZlSw9vv+1g8GCnuspFRCohhXY+DZMaAbD90LYwlyT0atQw+fjjLNq29fDeew7uvjsGlyvcpRIRkfwU2vk0qOJvaUdfaAMkJ8Ps2Vmcc46bjz6yc9NNuhxMRKQyUWjn0+Bo9/i2Q9HVPZ6fv6v8ggvczJ9v4/LL4/jzTy3AIiJSGSi086mXUB+rYY3K7vH84uNhxoxsbropl9RUKz17xvHzz3qriIiEm/4S52O32qmbWC9qu8fzs9lg3LgcnnzyCHv2GPzzn3G8+65dl4SJiISRQvs4Dao0ZE/WbjJdmeEuStgZBtx5p4t33skmJgYeeCCG++6LISsr3CUTEYlOCu3jNIjiGeTFufhiD/PnZ9K2rYf337dzzjmwdq3eOiIioaa/vMeJ9hnkxTntNJNPP82if/9cUlOhZ884XnjBocvCRERCSKF9nGi+VrskTieMGZPDF19ASorJ2LFOevWKU6tbRCRE9Nf2OP7Q3nLwtzCXpPLq2RMWL87k2mtdrFlj5aKL4nj4YSeHDoW7ZCIiJzeF9nGaJDclzhbPD7tXhbsolVpSErz22hE++CCLhg1Npkxx0LlzPNOn23G7w106EZGTk0L7ODaLjbNqn83G/Rs4EEUbh5RX9+4evv02k5Ejc8jKMhgyJIauXeP59FObLg8TEQkwhXYROtbuBMD3f68Mc0kig9MJ992Xy8qVmdx0Uy7btxsMGBBLjx5xzJlj0+YjIiIBotAuQsdTOgOw4q/lYS5JZKlVy+S553JYujSTq692sX69hYEDYzn33HjeecdOpi59FxGpEIV2Ec6qfTZWw8pKhXa5NGpk8sYbR1i2LJO+fXPZudPgwQdjOPPMBEaMcLJpk952IiLlob+eRUiwJ9C6RhvW7PmJbHd2uIsTsRo1MnnxxRx++CGTBx/MITbWZPJkB127xnPllbF88IFNu4iJiJSBQrsYHU/pjMvrYtmuJeEuSsQ75RSTYcNy+emnTKZMyea889wsX27jnntiadkygVtuiWHOHJuWRxURKYFCuxhXNbkGgHfWTwtzSU4edjtcfrmb2bOzWbEig6FDc6hf38vcuXYGDjwW4DNm2LUdqIhIEWzhLkBl1a7mWbSucSZfbv+cvzP/onb8KeEu0kmlUSOTIUNyefDBXDZutDBnjo3//c/O3Lm+L4AWLTxccIGHrl3dnH22h8TEMBdaRCTM1NIuhmEY3HRGfzymh3c3TA93cU5ahgEtWngZPjyXFSsyWbEig9Gjj3DBBW62bbPw2msOrrsujiZNEujRI45HHnHy6ac2/vjD0HXgIhJ1DNOs3H/69u5ND+jxUlISS33MDFcGbaY2w2l1sPi6VaTEpQS0LJGqLHVYEVlZsGKFleXLraxYYWX1aiu5uce6zatV89KqlZc2bTy0aeP73qCBiSUCPoqGqg5PZqrDwFA9Vlyg6zAlpfhuxZB2j3s8HkaOHMmOHTvIzc1l6NChnHPOOaEsQpkk2BMYds4IHl36MA9+ey/TLv0vhqGx1lCJi4MLLvB1kQMcOQJr1lhZudLKL79YWLvWyuLFNhYvtuV7jUmjRl6aNPHSuLGX00/3/dyokZf4+HD9JiIigRHS0P70009xOp289957bN68maFDh/LRRx+FsghlNrDNnXy1/Qu+2P45475/hqFnj1Bwh0lMDHTq5KFTp2NLrB0+DOvW+UL8l1+sbNxoYetWC+vWWQu9vkYNL/XqmZx6qpe6dU3q1fNy6qm+73Xreqla1dddLyJSWYU0tC+77DJ69uwJQHJyMpkRsESWxbDwao83uPzjnrzww1h2HN7OU+eOoXps9XAXTYAqVaBLFw9dungA3+beXi/8+afBb79Z2LLFkvd9504LqakWVq8uHOgATqdJSsqxr5o1vfl+NqlWzSQpySQ52aRqVZP4eIW8iIRWSEPb4XDk/Txt2jT++c9/hvL05VYn4VTm/XsB/T7vzexf3+eLbZ9zQ8t+XHX6vzkzpR1WS9EhIOFhsUDduiZ163o4//yCC597vbB3r8Effxj88YeFnTsNdu2y8McfFvbsMdizxyA11VJg7Lw4NpsvvH1fULWqL9QTE00SEnyhXvC7SUICxMeb1K8POTkGcXEmTqfvcjh9ABCRkgRtItqsWbOYNWtWgfvuvfdeunbtyrvvvss333zDG2+8USDIi+J2e7DZKkco5rhzmPjjRJ5e8jR7MvcAkORM4uxTz+aMlDNolNyI+kn1qZ1Qm5S4FJJjk6nirILNoivrIolpwqFD8PffsHu37+vvv2HfPjhwAPbvL/p7RbYktVh83f+xsb4v/88n+u5w+MLe/3X87bLcZ7OB1Xrsy2Ip/vaJHtMHD5HgCvns8VmzZjF37lzeeOMNYmJiSnx+OGePFyfHk8O3Oxcwb+tnrPhrGVsPbTnh82NtsUe/4nBanTitTmwWO3aLDZvFjs1iw2pYsRiWvC8Dw/fdMDAw8v4aGkf/A449lk9pxtuPf01RbBYbg84aQvNqLQo9ptmmhZkmZGbCwYMGGRkGGRmQmen7OTOTo9+P/ezxOEhLc5GdbXDkCEe/DHJyIDvb991/+8iRyElCwzALBfyx2ybG0bfyib78s/+LfvzYMex2K16v5wTP9R3LOPa/T5HHyf+8wr9PaX7nkl9T0nPK85ryH7fgn3yn005Ojiso5QulcJ3fYoFBg+ycccZJOHt8586dvPfee7z77rulCuzKyml10rPBpfRscCkAGbnp/HbgV35P38Ef6X+wN3sPB47s50DOATJy00nPTSfbncUR9xEyXBkcOLKfXK8Lt9eF2+vGY1bOvSu717ugyNCWwgwDEhJ83eBQ8ufglBQHe/ceKdWxvV7IzfUFe06OQXa2L8hdLl/r3vfdIDfXf9vIu993X8HnulyFb3s85H15vb4vj8c47nb+x489Vvi7Ucz9vu+m6bvG/vgv8H33n6+o55imkfc8AI/HUszzCh6vuOf4jyf2cBcgop1xhu8rFEIa2rNmzeLw4cPcfvvtefdNmTKlxC7yyi7BkUi7WmfRrtZZ5Xq9aZp4TA9e04vH9GCaJl68mKYX0zQxj4ZAgZ8xCy0uYpYiLErbsWK32Kgak1y2X0SCwt917vuc6//3q9TLK4SEr8cnMDvOeL2F7yv0/1cRVR6I5wTquKV5vKjj1qiRSFpa+gmfU/JxwvvhJ5yrjRgGNGuWQFpaaM4X0tAePHgwgwcPDuUpI4JhGNgMjXuLhEskLMgTLFWr+npcKia6P0SGsms+it+qIiIikUWhLSIiEiEU2iIiIhFCoS0iIhIhFNoiIiIRQqEtIiISIRTaIiIiEUKhLSIiEiEU2iIiIhFCoS0iIhIhFNoiIiIRIuRbc4qIiEj5qKUtIiISIRTaIiIiEUKhLSIiEiEU2iIiIhFCoS0iIhIhFNoiIiIRwhbuAoTS+PHjWbFiBbm5uTz55JO0bt063EWq9FauXMn9999PkyZNAGjatCl33XUXQ4cOJT09ndq1a/P888/jcDjCXNLK6ddff+Wuu+7i5ptvpm/fvuzbt6/Iuvv666+ZPHkyOTk59O3bl2uuuSbcRa80jq/D4cOHk5qaStWqVQEYMGAA3bt3Vx2ewIsvvsjKlStxuVwMHDiQc845R+/DMjq+DhctWhSe96EZJZYvX24OGDDANE3T3LRpk9mnT58wlygyrFixwrz33nsL3Dd06FBz7ty5pmma5pgxY8xZs2aFo2iVXmZmptm3b19z5MiR5jvvvGOaZtF1l56ebl544YXm4cOHzaysLLNnz55mRkZGOIteaRRVh8OGDTMXLFhQ4Hmqw+KtWrXKvPXWW03TNM0DBw6YXbt21fuwjIqqw3C9D6Ome3zlypX06NED8LUW9+zZQ3Z2dphLFZlWrVrFBRdcAECPHj347rvvwlyiysnhcDBp0iRq1qyZd19Rdbd27VpatWpFYmIisbGxtG/fnh9++CFcxa5UiqrDoqgOi9euXTvGjx8PQGJiIi6XixUrVuh9WAZF1aHX6y30vFDUYdSE9t69e6lWrVre7WrVqpGWlhbGEkWOzZs3c+utt3L99dezdOlSMjMziYmJAVSPJ2Kz2fLqya+oujv+vVm9enXV6VFF1SHAjBkz6Nu3L4MGDWL//v2qwxOw2WzEx8cDMHv2bLp160Z2drbeh2VQVB1aLJawvA+jZkzbbrcXuG2aJoZhhKk0kaNBgwbceeed9OrVi127dtGvXz/MfCvfqh7LJv/70F93em+WzZVXXklCQgKtW7dmypQpvPLKK3Ts2LHAc1SHhc2fP58PPviAt99+myVLluTdr/dh6eWvw9TU1LC8D6OmpZ2SksK+ffvybu/fv58aNWqEsUSRoVatWlx++eVYLBbq1atHjRo1yMrKyhtaSEtLK7HrUo6Jj48vVHfHvzdVpyfWuXPnvEmk3bt3Z/PmzarDEixZsoQJEyYwefJkqlSpovdhORxfh+F6H0ZNaJ933nl88803AKSmplKvXr0iu92koLlz5/Lqq68Cvg86+/bt45prrsmry6+//ppu3bqFs4gRpWvXroXqrk2bNmzatIn09HQyMzP5+eef6dChQ5hLWnndf//9bNy4EYAff/yRJk2aqA5PID09nTFjxjBx4kSSk5MBvQ/Lqqg6DNf7MKp2+XruuedYtmwZVquV0aNH06xZs3AXqdLLzMxk6NCh7Nu3D9M0ueuuu2jRogUPPvggWVlZNGzYkDFjxmCzRc1IS6mtW7eOsWPHsmvXLmw2G7Vq1eL5559nyJAhhepu3rx5vP7661gsFm699Vb++c9/hrv4lUJRdXjfffcxduxYYmNjiY+P55lnnqFatWqqw2K8//77vPrqqzRs2DDvvjFjxjB8+HC9D0upqDoM1/swqkJbREQkkkVN97iIiEikU2iLiIhECIW2iIhIhFBoi4iIRAiFtoiISIRQaIuEwP79+xk2bBhXXnkl1113HX369GH+/Pklvm737t0sX74cgI8++ohZs2aV+pyjR49m3bp1ZS5rRc5ZVnPmzAnasUVORrrkSyTITNPk+uuv57LLLqNfv34A/P333wwcOJBBgwblbWRTlE8++YQtW7bwwAMPhKq4ITunx+Phsssu48svvwzqeUROJloRQyTIVqxYgdfrzQtsgNq1azNkyBBef/11evToQZ8+fTjzzDNZu3Yt+/fv55577qF169aMHz8e0zSpWrUqGRkZuN1u7rvvPjp06MCdd97JggULcLlc3H777XzwwQds27aNp556ii5dunDjjTdy55134nA4ePnllwHYt28fHo+HL7/8kjVr1vDss8/icDjIycnh8ccfp0qVKkWe84EHHuDbb7/ltddew+l04nQ6GT16NLVr16Zbt27ccsstLFq0iF27dvH4449z7rnnFqiD6dOn87///Y/4+Hjsdjvjxo3jueeeY9euXfTv35+33nqLzz77jHfffTdvc4bRo0dTtWpVWrduzZ133snSpUvJyclh3LhxNGnShOeff54VK1bgdDqpVq0aY8eOJS4uLqT/tiIhF9CNPkWkkGnTppmjRo0qdP++ffvMtm3bmqZpmn379jWfeeYZ0zRNc+vWrWbnzp1Nj8djvvLKK+aLL75omqZZ4OemTZuay5Yty3vtww8/bJqmaX7wwQd5+5/37dvXXLp0ad75XC6Xee2115pffPGFaZqmOX/+fHP16tWmaZrmJ598kve6os6ZlZVldunSxdy1a1fe7zR06FDTNE3z/PPPN2fMmGGapml+/PHH5h133FHod23fvr25Z88e0zR9e7SvXbvW3Llzp9m1a1fTNE3z77//Ni+//HIzJyfHNE3TnDFjRl6dNW3a1Jw/f36B3+/gwYNm27ZtTbfbbZqmaX799dfmjh07SvqnEIl4ammLBJnb7cbj8RS63zRNLJZj00r8OwT5l0rcv3//CY/bvn17wLepS7t27fJ+Tk9PL/L5r776Kk2bNqVnz54AJCUlMX78eNxuN4cPHyYpKanYc23fvp2UlBTq1KkD+DbteP/99/Me79SpEwCnnHIKhw4dKvT6Xr16MWDAAHr27MnFF19MkyZN+OOPP/IeX7t2LXv37mXAgAEAuFyuAhst+Oumffv2vPXWWyQlJXHWWWdxww03cNFFF9GzZ0/q1q17wvoSORkotEWCrHnz5nz22WeF7l+3bh0tW7bMu23mm17i9XpL3NLParUW+bNZxDSV77//nq+//prZs2fn3Td06FBGjRrFP/7xD+bPn8+0adNK9wtR+ANH/rXnizr/qFGj2LFjB4sWLeLOO+9k8ODBtGnTpsBz2rRpw5tvvnnC8+avl8mTJ7Nx40YWL15Mnz59eOmllzjrrLNK/TuIRCLNHhcJso4dO2K1WpkyZUrefbt37+aFF17gtttuy7vPP2N727Zt2Gw2qlWrhmEYuN3uCp3/0KFDjBgxgnHjxhUY8z148CDNmzfH6/Xy+eefk5ubC1DkORs0aEBaWhp//fUX4Num8Mwzzyz1+f/zn/9Qv359+vXrx9VXX83q1auxWCx552ndujW//PILaWlpAHz11VcFJqj562b16tU0a9aMnTt3MnXqVJo3b85tt91G165dWbt2bTlrSCRyqKUtEmRWq5XJkyfz3HPPcdVVV+F0OrFYLNx///107do173lut5tbbrmFPXv28Oijj2IYBh06dGDw4ME4nc4SW97FmTlzJocOHWLs2LEAWCwWpk2bxh133MFNN91E3bp1ufnmmxk+fDjvvPNOkeeMjY1l9OjR3HfffTgcDhISEhg9enSpzp+UlER2djbXXnstsbGxADzzzDPUrFmTGjVqcM011zB9+nQeeeQRbr/9dmJiYnA4HIwZMybvGKmpqUydOpXs7Gyee+45ateuzaZNm+jduzd2u524uDiGDBlSrvoRiSS65EukEvDP9O7SpUu4i1LpNGvWjNTUVG3/KoK6x0VERCKGWtoiIiIRQi1tERGRCKHQFhERiRAKbRERkQih0BYREYkQCm0REZEIodAWERGJEP8PYoKzvE1QHC4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x396 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.style.use(\"seaborn\")\n",
    "plt.plot(gd_cost_history, \"b\", label=\"Gradient descent\")\n",
    "plt.plot(qngd_cost_history, \"g\", label=\"Quantum natural gradient descent\")\n",
    "\n",
    "plt.ylabel(\"Cost function value\")\n",
    "plt.xlabel(\"Optimization steps\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can visualize the optimization path in the parameter space using a contour plot.\n",
    "Energies at different grid points have been pre-computed, and they can be downloaded by\n",
    "clicking :download:`here<../demonstrations/vqe_qng/param_landscape.npy>`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_surface(theta0,theta1,cost_fn):\n",
    "    parameter_landscape = np.zeros((theta0.shape[0],theta1.shape[0]))\n",
    "    for i in range(theta0.shape[0]):\n",
    "        parameter_landscape[i,:] = cost_fn([theta0[i],theta1[:]])\n",
    "        #for j in range(theta1.shape[0]):\n",
    "        #    parameter_landscape[i,j] = cost_fn([theta0[i],theta1[j]])\n",
    "    return parameter_landscape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAF3CAYAAABE/OtwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA/1klEQVR4nO2dd3hUZfq/75kUEtIgARIRNAhSJCDFigVWxPKzsyrggquu7Cq76Bd2XcGGiA1sCDbEFVFXpKwrqyJIWVZAQExQCIQSpHdCSEhC6szvjzBhkswkU04/z+3ldYWZM+e8k3KfzzznfZ/jcLvdbgRBEARL49R7AIIgCIL6iOwFQRBsgMheEATBBojsBUEQbIDIXhAEwQaI7AVBEGxApN4D8MdtD28NaLsuvc5VeSThkdG5iar775Kar8p+08u3qLJfb+J2ZKl+DKGa4va9VD/GrujOquw353BzVfYLkL2lTLV9K0FO1u6gtp//bie/zxlW9t4YXei+EMk3jIheO8wservjy33BngA8GFb2ZhQ8qC95ENEL9kHNVA/Vf69GT/d1CdWNhpW92TCz5MHeZRvHgSPgcOA+q6XeQ1EcLVK9YA5E9mGiheRBRK8mzqzNgIOqm/vqPRRTIiUccyCyDxGRfOAYWfRgXdlLqhe8EdkHiRUkDyL6GsrKcWbnVn9dXgHRUfqORyG0Er2aqV7ter0HM9btQ0FkHyBaSR4kzWuJM3s7jvKKmq9dvS7QeUSCoA4i+wbQUvAgaV4PnFk5Xl9vtoTspXwTPHZI97KC1gcZnZtonuTVTvNaXYQ1k+jBU68//XXm5ga2DJz9xwu58O9T2bjnUK3H75k6m6dnL+bv/1xIaUWl39c/+tHXioxDbeTCrLmQZH8arVM8WCfJg3HTvHPVeqLen4OjoKjxbQ8eJeb2kT6fcyfFU/HHu3Fd0TOg47ZJTuS7Dbl0OycNgAP5hRQUlwIw6Xc3NPjaN++7OaBj+EJSfehYPd1rKvuvvvqKGTNm4Ha7efTRR+nXr5+Wh6+HHoIH9SUPUrLx4LqiJ2UZHYiaNpeIH9aHtI+qPj2p+NNdkJQQ8Gu6n3sWa3P31vx78YZc+nQ6l9LyCm586SP+NfoeZizPoqSsgp1Hj7P3WAFjbu/LFZ3Ope+z0/nfs8N58vPvaJEYR/aew5wqr2DIFd35cl0OZZWVfPjQb/n2523kHsrjrzdfSUlZOb99/TPmfdCL56ZMJ6V5Ejnbd1JSWspd/+9avlm2krLyct6eMIboqPAvQls11VtZ+JqVcYqLi5kxYwazZs3ivffeY8mSJVoduhaeEo1eSV6LNC+ir0NSAhV/f4Dyxx7AnRQf8MvcSfGUP/YAFX9/ICjRA0Q6HXRq3ZINu6tLOd/n7OLKTrVXPjodDg4VnOSdP9zG47ddzdw12fWej3I6+cdDA0lqGsPevAL+8dBAEmKakL23dokIwBUZXfO6yIgI3nrucZIS4tl/6AhvPfc48XFN2bx9Z1Dvw47oFQLVRrNkv3LlSvr27UuTJk1ITU3l+eef1+rQuv/wrJTkwUSSr0MwKT+UNF+X67p1YNGG7bRKiiOxaROaNqmfqHumtwYgtVkCJ0/VT5QXtGkFQIuEpnRu3bLm66LS8gaP3bl9OwBSmiVxfrtza74uPnUq5PdjJ6yY8DWT/cGDBzl+/DgPPvggJSUljBw5kssvv1yVY+ktd9BG8CCSD5qkBCqH3tyo7CuH3RKW6AEuO78tUxb+QOvmCfTPaO9zmwino8F9RDidPr92u8H7lYXndgfm1fw7MsLp82vc7sAGL1hO+JrJvry8nLy8PKZNm8bevXu57777WLZsGU5n+JUkI8jdgxUlDxYR/Wm8Z+D43SZzc9graqMiI+jcuiX//nEzMx7+LVsOHA1rf3WJj4nmREn1Rd+tO3Ypuu+GsGq93hcet1hB+prJvmXLlvTo0YOIiAjS09OJj48nPz+flJSUoPdlJLl7EMmbh4BknxW+7AEGdO9AfvEpEmKV/529tENbPvpfFn/7z4+c3+4cCe0qYgXpO9xubX5FDh8+zJNPPsn06dM5fvw4d955J0uXLvWb7MdML6352ohyB+0EDyJ5xSgrp8mwMTWrZj1TKoFaUzTd0VGUfTrRFO0TtJ5uqVWy16pdQigYVfovD4/x+5xmyT41NZUBAwZw7733UlxczFNPPdVgCceoggeRvJnxbo9Q9yKs98VbR3mFtE/wgZ1KOA3h7ScjiD8QX2qW7IPl0xXGGpaWggeRvFpETp9HxMrMBhdIeRZiVV3Vm8oH79R4hMFh1VTvwcjp3h9qyT8QoQ+9yv8Ff1lB2wgieYsRE03ZlCcanGnjmaIZ+dVyzYYlWAejViVE9j7QWvAgkteKymG3BrZhUgKVQ29RdzBhIq0RhGAQ2Z/GDoIHa0q+KLPhOfPxvQPrZyM0jNTrzY2tZa+H4EEkHyqNST2Q11lF/JLqhWCxnez1EjxIqSZYQpV7IPtUU/p7jp3gla9WcLzoFG7cnN08kSfu6EfzuNig9/Xa1yvpkJZCi4Sm7D9eyN2Xdwvq9SvWreeyHt2IivL9p15yqpSho57ii/deDXps4bB21TIuveKaRrfrkppvyou0RsQWsreT4MG8kldD7g0dK753T1xuF9+dzOHnkr30aNqW6xK64HSEvqq7yuXib598y+O3XU3v884GYMbyTCb9ZwUvDbmuZjuXy42zkVYJ3lxRp4laoHz+1SJ6Z3TxK/tAUbKEU1lRwYIvZwUke0E5LCt7PQUPIvlA0FLuvijMzOIPzX/i68KNNY/dnNiN2e0eDFn4q7ft4bzU5BrRA/z+6l64qZ5K/PTsxURGOMkvLuX5QQMY89kiTpVXUFpRyZjbrqbbOWl8nbWFj5ZnkZoUT3RUBB3SUpj/U05NO+OPt+WzZMaLuFxurulzMYNuvo7ps/5NSWkpu/cfZN/Bw4x+cCj5BSfZtG0Hf33hdaaM+3uN8ItLTvHEq29TXl7OhV061ozzl5xtvPfPfxHhdJLaMpmxIx7g2PETjH9zGmXOOCorKxkx+hlatExj+lsvcezIIRxOJyNGP0NSs2Q+fOcVDh/aT2VFOXcN/RNdu/fmuTEP073npWT/so6ik4U89syr/Gfex+zd/SsfvvMKD4x4LKTvsxA8lpK9HQUP5pK83oL3Zin7aoke4OvCjXx3MocbEruGtM9dR09wflrtFiDVCf5Mim8WF8u4O/uz++gJbu7ViRt6dGTdjn3MWJ7Fa8Nu5K2Fa5j16CASYpowZMrntfZ1IL+Q5WszeWfCWAD+9OQL/Obyi3A4HRw+dpzXnxrN6qwNfPndcl5+fCTTP/+C154cXSvZL/x+NR3ObcPI3w9m2Q/r+G7FGgAmf/gZb457jMT4ON6aOZtlP6wjL7+Ai7t35Zqhj7Nj22ZO5OexbctGEhKb8fCoZ/jh+8VkrllBbNM4kponM3zkWAoLTvDiUyN5eeonOBwOYpvG8dSLbzPro7f5cfVybh44lNxtmwMWvZRylMH0stdb8CApviGMJPe6bHDk+Xz8l5J9Icu+yuWiymud4qMffU1RaTmHC4qYO2oIAF1Pty1uFhfD/zbvZPbqjZRVVBIbHcWJklJim0TV1Pd7nG6B7OHnqgT27D/EX8ZNBKpr7gcOHwPgwi7nA9CqRTJFxSV+x7hr7356dq0uy/To2gmAwqJi9hw4xNhJUwE4VVpG86RELu2RwZhJUzlQNoWLL+9L5649WLtyGZ0zegDQ5+oBAMx8/w02/fITWzf9AkB5eRmVFdUrlTt1vRCA5BatKDpZGNw3VFAMU8reroIHc0jeyIL3prvbdxO+C5u2CXmfHdJSat2ExHOLwRtf+gjPYvWoiAgA/rnyZ1olxfPSPdezYc8h3lzwQ73Wxb4WuF/WqxtjH76/1mM/bdxcuwUy/legu+HMQbz2n9IsibefG1Nv+49fe46vso/zyQdv0m9A9ftxu1z1trv1rnu5st/19R6POP1+/b2fQJB0Hz6mueG45y5PeoteyztBeWP0m3kXZa6v+d8s9KcNN7jPqfXYzYnduC6hS8j7vKRDG/YfL+D7nDN3hPpl90GKyyrq9YIqKCmjbUoSUH3bwoqqKpo1jaGotJzCU2VUVFWxftfBWq/pdF46WdlbKC0rw+1288Y//klZmf8bmTgdTqrqiPmc1mk1LZF/2pgDQGJ8HE6nkx179gEwd8EScnftZfHKtaw50oRLr7iGmwf+jl+359Du/M7kZFf/nLN+XMmXsz+iQ8cLyFzzffX7OnGc2R+/53dMDocDV1WV3+cFdTB0stdb7B4kxfvGTGL3hRMHM139Wco+NjryuKTDJWHPxomKiOCtB25l0n++5/2l66iodJGSEMtb999CTJ0ZMbdd1IUnP/+OJdk7GNynO//d9CvfrN/Cw9ddygPv/ouzkxPpkJqMy+3G6XBQ0awVaS1TGHzL9Yx4+mUcDgdXX9KLJk2i/Y6nZ0Yn/vz0y0x59jGaJVa3iLix3xWMnTiVkc9O4sIuHWvS9tgR9/PSOzOIjIigZUpzbhvQlypXFS9Me53o6GgqKyu5/+G/cXabdLJ/XsdzYx7G4XAwYvQ4miWnsGlDJuMeG05VVRV33jPc75iaJbegsrKCKROf4pHHA79jnaT78DBsI7TMbcf1HoJI3gdmF3xjGHnRlV4LqYy0clZk7xtPMO7dMdnvNoZO9nohF1xrY3XBC+ZB0v0Zgq18iOxPIym+NiJ4YyGp/gx2FX64ZW3by14kfwYR/JmVtYJgBJS8bmlL2YvgzyCCF8yIVdO9mpNSbCV7kXw1InhzIR0ufWMV4Ws169AWshfJi+CF4DFivb4uZhS+XlPKTbOoKlg8i5/svgDKbAudjEC43699eQWMnPEV90yZzT1TZvP8F/+l5PRNzpVg/aatHC/Qp+3A2lXLwt7HI3+4g9JT/ts5BItnseW2X/7L/A/HUFRwlG8+fibo/Rzeu5W8Qzsb3OaDCb/lxLF9IY8xFNFnrl1R03oiHCwne70ED8aRvBlXs+qFy+1m8fa9vLbiZxZv34srzGUnLpeb0Z8sYOhVPfjskUF89sgg2qYkMW7OEoVGDN8sW0G+TrL/z7xPAtrO5aOdQjgEshyobbMiAOKTWnLTvc8FfYwtWd+Rd3hX0K/zhbfcw03yC76cRWVl+LK3TBlHL8GDcco1IvfgcLnd/H7OUhZu21Pz2A0dz2Hm3f1D3ufq7XtIb9mcSzu0rXns3qt7cusrn5J3soTJC1ZxbbcO9L2gHf/bvJMlG3OZMGgAH3+/nsUbcnG53VzV+VweGnAp/163mczj5eQXrGD3/oMMu+P/0apFCt//mMWve/fz4mN/4f7HnuXbj6qblz3xytvceWN/XG4X8xYsJSLCydYdu3lw8O0sX5PJjj37GPPQffTudqYdxP5DR3jh7Q85O60l23buoUv7dox5+D527N7LhOmvExkZhcPh4NExL/DfxV+xZ+d23nhxDNffcjfffT2P/xv7IgB/vOcG3v9sIe+9MYGIyEhOFp7goVHP8PYrz1Baeory8jJ+/6fRdOjou8HcyuWL+Grex7RomUZ8YhJdMnricDj5+acfOJGfx18ee46F8z9n+9ZsKioquPaGO/jN9beyZ1cu777xHPHxibQ6qw3NYio4cWwf8959lAef/hd7tmfy3y/ewOF0kpR8Frfc9wK7t/3ET8s/w+FwcvRALl0vvpFOPa8l83+zaZr1HXEJKZx9XveasS387Hn2/fozKantqDot3ZMnDvP1zGeorCjDGRHBqL8+RlLzFN557Vnyjx+jrPQUA4f8gYsuu5pv/v0Z61b/j6qqSgb/fgRdu/fmu2/msfr7JbjdLi694hpuvG0wcz99n9JTJRzYt5tDB/dx359GU1hwgtytm5j47GiefH4qkVFRIf9uml72dq/Hi+BDZ2nuvlqiB1i4bQ9Lc/dxx0WhXRTddTSfTme1qPWYw+GgQ1oyvx7xvyrc7XYz/U930CQygpsmfszQq3ridDjYsWc/0154kr0HDzFu8jRmvjqe89udw+gHh5LW0ncjN6fDSe6uvcya8iLL12YyfdYXfPbmC/x3TSaLvl9dS/ZOp5MtO3YyYfTDJCXG89uHHuNkcQn5BSf53QOPcH7nDOZ99gGrln/HLQOH8tW8Txn1xMts3uj/9z8hMYnhI8dy8MBerrzmRi6/6lo2b8w6/dqXfL732TPf5YXJM4iJbcrjfxnKBd16A3A87wjPTnqfysoKmiW3YNzEaVSUlzHqj3fxm+tv5d+fz+Cu3/2RXpdcyYx3X6GqooIOLQqJiaruvbNo1gsM/esMYuOSWDJ3Ept/Wkh8UksO/LqBEc9/i8tVxVtPXMfVt/yZ9hlX0qX39bVEf/RALvt//YU/PDmXkyeOMHXMtQD8vOhV7h70WzIuvJj161bx79kz6H/DHRQWnmDcxPc4WVjA+nUrOXLoAGtXLePZSe9z5NB+/jPvY1qltWbdD8t55uV3AXj273/kkiuuweF0knfsCI+Pf4Off1rN0oVfMvrJicz95/s8/uzrYYkeTCx7uyd5kXz4bDjku8XxxkN53BHiPssrq2q1OPbgdrtrdaWsS2REBA9N/5IIp5P8olMUlJQCkNGxPRERTlqlNNy2uC7tz21DVFQkKc2SaNf2bCIjq78uPnWq3rZtz0olpXl1Q7aU5kkUFZdQlJzBvI/eobKinON5R7nCRzdLf5x3/gUAJCQkkrl2BYu/+RflZWU0iYnxuf3JwgJiYpuSmFR9obXTBRee2VeHLjgcDqKioik6Wcj4xx8iIiKCwoLq0sj+vbvo0Kn600KXjJ78/NPqM+8rbjf5R3byzfSHKCmPpLz8FE0TUohPaknauV2JalLdRtrt9l9yOnogl64XnM8FaScgLZp5rc+mQ4tC5m/P4cC+3fz78xm4XC4Sk5pzVptzKCk6yTuvjeeiy67min7Xk7l2Je06dMbpdJLWui1/fORJ1q1ezsH9e3j+iT8DUHqqhKOHD5x+79UnmpSWrSgpLgr4ex4IppO9SF4krxTd03wn425+Hg+Edq2S+XLdplqPuVxudhzOJ71lc3CcaWDs6Ua5L6+AWat+Yfajg4mLieb2Vz+t2SYiwqttcSN16yqvTpLebYVrtxiu/zrv5z18/P4b3HLnMHpe1If/zPuE8vKyWs/XvaGiy3Xm2JGR1VpZ+J85JKe05C9/G0/u1k3M+uhtPyN34/D6vnjfrtGzr00bMsnJzuKpF9/B6XTyh9OlturvieP0XurTrHkKT7/0Tq3HNm/MYmtTV00tPcLhpktqPv+LKadts6JaNfaCpCIKvcbm/TN45PEJJKe0qrXv51//kC2bf2Hpwi9Zs2oZl13ZH7er/sgu7H05w0eOrfVY9i8/4VSgHbQ/THOB1s4XXuWCqzr079CGGzrWbnF8Q8dz6N+hTcjf68vPb8v2g3m1Whx/smI93c5JJTk+loSY6JrUvnn/UQAKT5WRHB9LXEw0G3Yf4vCJIiqqqihr6fu+sw6Hg6qq6hNFRIST8ooKKquq2LE7+Fki/iguKiTtrDaUl5exft2qmguEnhQcGxdP0ckCAA4f3OczhRYVFZJ6VvW9AdauXOr3ImN8QhLFRYUUFRVSXlZa0z651nhOFtKiZRqRkZGsW70cl9tFZUUFrducw84d1V7Y9MtPtfcbn4jT6WTv7h0ALPpqLnt25vp9zw6no9YJE6D12eeya8dW3G43eUcPc/hg9fe4faeuZK5ZUXPcH/73HTtzt7J6xRIu6NaL3/9xFL9uzyH9vI7kbs3GVVVFQf5xXnv+cdq178zmjZmUlZbidruZ+f4blJeV+h+Xw1lvXKFgimRv17q8yF1dnA4HM+/uz9LcfWw8lEe3tBT6d2iD0xH4jcDr0iQqkin338LzX/yXtxatweVy0+2cVJ75bfXNtW/u1Zln5y0lc+cB4ptE43K76dS6BfEx0dz3zjwuPPcsBvXpxsvzv+eaAQN8HqNn10488/q7vPjYX7j9ut/w6PhXade2NeltW+Nyu8Jq0ezhhlsG8caLY0k962xuuuMePn7/DfpcNYD08zryzN+G8+ykaURGRfHGi2NJa92WFq3S6s3Aubr/Tbz7+nh+/GE51938W35a+z3fL1tQ71hOp5M7Bt/P+L8/xNnnpNP+/AuIiIioJbiMnpfw9b//yfNP/oWeF/Xhosv6MmPaa9x+9/1Me/N5Fn01l1Zpret9Ahk+8gnen/IiERGRJLdoxTU33Mb2Ldl1hwBAp649+PSDN4lt2pSMCy8G4Jx2HWh7bnvGPTacs84+h/T2nXC5Xdw55A+8N/l5Vq+onmX10P89TVx8AnM+eY+lC7+ksrKCO383nFZprbmi3/U8+/hDuN0u7h72J1q0SuP/3TaECWNHgAMuuuxqopv4LnEBXNCtJxPGjuCJ56eSmNQsoJ+fLwzd4lgkL+hJuD1yNuw5xMT53/PJn++qVZoIFDs1P1u7ahldu/cmPiGJl575P357zx/o2Lmb5uMwO6ZscWxH0YvkrUX3c9LokX4Wg9/8nOH9L2ZA9w56D8mwlJWV8sKTI4lu0oT09p1E9Cpg2GSfl/2DpscTyQt10bP7pZ79cMzQJkGoT3r5FlIy+vh93rDJXktE9IJQjYjefARaBbG17EXygiCYmWDK3aaZeqk0InpBEMxMsNc1bZns9RK9SF4QBCUIZQKL7ZK9iF4wOnaacikET6gzFW0lez1ELytfzYv83ASjEc6UdNvIXi/RC4IgKEG4a49sIXsRvSAIdscWstcaEb1gNqReb2yU6Chg+dk4eve5MTs7l21qfCMv2l3j+05EgiCEhlKtYywteynfBE+wcm/o9SL+4NGzTYJgbSwte60xq+jDFXwg+xXxGxcp4dgDqdnbmJ3LNqkmel/HEgQhOJTs/qtJsl+7di2PPvoo559/PgAdO3bk6aefVvWYWpdwzJbq9ZDvzmWbJOELgk5oVsa55JJLmDJlilaHE/ygd8IW4QuCPkgZRwHMkur1Fr0Ho4xDkHq9ndBM9rm5uTz44IMMGTKEVatWaXVYwaCI8AVBWzQp46Snp/Pwww9z0003sX//fu69914WLVpEdHS0FocXELmaAZl2KaiJJsk+NTWVW265BafTSdu2bWnRogWHDx/W4tACxhW9UcdlF6SEY3yU/BlpIvtvvvmGqVOnAnD8+HHy8vJITU3V4tCCIAgCGpVx+vXrx4IFCxg8eDBut5tx48ZJCUcQBEFDNJF9XFwcb7/9thaHEkyGTMUUhIbZFd1ZkcVVlp16KRe7BME/Uq+3H9IbRwHie/c0zVx7QbAiOYeb13usS2q+DiNRByXSvaVlX9y+l7Q4proJmcx8EayCL7EHup2ZTwDhCt/SstcSSfeCoB6BCj6Y/ZhZ/KFg2Zq9B6ndVyMXQQUPZqrX5xxurpjo/e1brf2rQTg/O8vLHrQTfnzvnpocJ1SMKHwjjknQH60lbCbxhyp8W8heS0T4ghAeegvXDNIPRfi2kb2Uc84gwjce8vtpPMkabTx1CVb4tpE9SDnHGyMI3whjsBtGrdcbWapGln4wP09byR5E+N6IbAUjYFSR1sWodf1AhW872YMI35t213TVRfpyohGMKM5AMdrYAxG+LWUPIvy6aCl9Eb1gJFGGg5Gk35jwDSt7LWqLIvz6qC19Eb1gFDkqiZGk7w+H2+126z0IX2RuOw6gSLe3xtCqpYJZV9gq0WrBrJLX6kStRfAwwsVZowtRKfRandu7Y7Lf5wwvexDhG41A5W9WwXsjslcOu4jeG62l35DsTdEbx/NLqqb0PX9sakvfIw8zS98KEhcELfCc4IzQh8ewNXtfSB1fEMyJHVO9N0ao6ZtK9lAtfLWlL8I3FpVuF5O2ZzFwxbdM2p5Fpdul95CEINBbckZCT+mbombvD6nlW59Kt4sL587myNZTNY+16hTLL3cNItKhflaRmn14iOgbRunyTkM1e9Mle28k5Vuf13N/riV6gCNbT/F67s/6DMik6H1xVvCNlknfFBdoG0OpG/L6Q6s7XskNUOqz5tBhn49/sDwH1wE37Zolck7zBNLim5Ka0JSmUfV/pV24Wco+Njjy6O5OoT9tcOJQe+j1x+Fy8V1WDj//upce57Xlul5dcDqDy1sul4ulK9ewIWcb3bt0pP+VlwW9D62QVB84WlzItYTsQf0ZOzJbRx8uS0tlBQfrPV55rosp7g0k7owmMstJZZGbk0XlRDsjSI2PJTWhKanxsSQ3jWF50/3siC2EWCAWLo1JZVLU5SRERxMXHUlcdBTREU4cDnVOAC63i4UFm3ni9S/ZuuHMyevmS7ox6/EHWFK8lR9zNtC9WUf6p10GwNJDa9hwYhsZSR1wOBxsPLGdCxLOY/Jzn5C5enPNPm7odyUfv/mS4YQvog8NNe+kZeqavT/UruVreV9bu0u/oZp9saOCDeTxiyOPDRzjZ45xqKyEs0/G0fJkLM2KmlBYUs6q0kNwiur/S4ESOLssDkcFFJdXUlxeQZXLTdPoSGIiI4iOiCA60klFhIuq2AjiI2NIjmqK0+kkwukkwln7pOB2u8mrLOZkZSmxRJPgiMHlclFZ5aKiqoo9ZfkUl5RBUf3354hw4Ha7wU31/w5wOMDtACKojmPRQMzpFxyqv48rHumJw+HgrBMtueOS/lx75eX15K91GUdkrxzBSN/0i6pCxSrSF+G7eD33Z9YePsylqamM7tDD78XZIirYQQG5jgJyKeAbx25yHPX/WLq5k+nvbkMqTWlGE+IqI4kpj6BpVRRRlU6edf/I6qrDUAm4oE/T83g27WbcLqhyufB8CHC53Uw4tIA1JTvBATjhkrh0HjtrAFVOF+tKd/Pa0SXwE7DWx4D7AldR/VrPOcRdfUyqqD5+OdUnqTXABh/7aAqUnPln+05tWTHrY6Kjomse01L2Inp1CET6tpW9BzWlLynf2CxmL/dELK73+B9dF5BMDIcp4QRlnHCUU0AZJyjnGKcodFTUe01SRCxNHJG43G48/5W5Kyl2ldfbtkVEPMmRTSmsKuVQZSFsAz7zMcB7gI4Bvhl/+/BFLGxe/G9aJbUCRPZWoiHp2172ICnfrrhw83vnUhY69tQ8doP7HGa6+vu9SPua42dedtb/ef615bWMbPUbHDhwOMCJkylHljLpSP2TybNpN/N42vUsLNzEHb++V53UZwNbz2xzcc901t2yK/A5cT72QQvgmP+XbFk1nxaJLTSTvYheO3xJ3/TtEpRAiwu4Ws3YAZF+oDhxMNPVn6XsY6Mjj24BzMbp7k7x+fiVCR1IjUqs9djl8e3Bh+wvbNoGgOsSunBzYje+LtwIg4Bc6JSfyou9b+O6nl0YsvvD6udOc33aFTgcDhYeXOnrzdTsg0NAGtUln1n+3/8FN97OoRXf+99AMC3BzuCxTbKvi5R2BH/4+jRwc2I3Zrd7EGedawUut4tBOz+oJey627rcLr47mcMvJfu4sGkbrkvo4vO5Tuf1rTUbZ2PBdromtsfhcJBdkEvXxPa4cTM++122n9ztGSz8A9jv//3MevsVzr/2gfC/MY0gqV5fuqTmSxnHH1LaEfzhmZu/7ewm9QRdb9sGZB4MgS7gc7ldNVMzyyvLeS1nJrxE9cVcHzwx8o8M/MtLQY8nWET2+jP0Kv+fWA0r+09XVA9Li25xkvIFf2i5sjnU1dqllaW0ebl/dT3fB5Ls7UNDsjd8zV6LlWVq1vO1WowFUs+3KzGRMex6bBHpWdfD9trPXde3D/2vvIw9vl+qGGYTffaWspBfm9G5iYIj0Q7Dy96DSD9wRPr2I75JPPtmL+Wuj0azYdU2UqKTmDDkEW7se5XhVtfqQThyD3RfRj8JGL6M4w+1yztS2hHAHGWcQFB76qURk72Sgg8VrU8Api7j+EPtpC8pX5BOpIFhJNEbQfDe1B2PnunftLL3YHbpa5XyRfqC1TGa6H2hp/xNL3sPWkjf7CkfRPqC9TCD5P3hPXa1xW8Z2XtQU/pWKe2ASF+wBmYWfV3UFr9lL9OreQcYNe+QpdWdsTxIXVoIB73q9dlbyiwl+rp43p+S79Fyyb4uaid9SfmCoC1WlrwvlEr8lk32dVEr6Vst5UvS1w8tp+SaFbuJvi7hpH3LJ/u6qJX01arna53yQZK+YEzsLnpvQkn7tkn2dVGznq8Gxe17SdIXBKEegaZ9w8pei7O4lHYCQ4Qv6I2k+sZp7HukqexLS0u59tpr+eKLLwLaXqsr7maTvqR89bHTezU6Inpl0FT27777LklJSUG/zgrSVwO9Ur6IUBDMh2ay37FjBzt27KBfv34h70ONuae+kJTfOCJ9c6H2jXoE46OZ7CdNmsSYMWMU25/a0peUHxgifWWR6Ze1kRKOcmgi+y+//JKLLrqINm3aKL5vM0rfaikfRPqCYHQ0mWe/fPly9u3bx+LFizl06BDR0dGkpaXRp08fxY7hEb5azYRyDjdXZW6+2btp1iW+d0+Zny8ogqR6ZdFE9pMnT675eurUqZx99tmKit6b7C1lqgoflF2QZZWWC96YfVGWfEIRrIhh59mHgxalHSWx2rx8D1LaEQTjYNjbEo6ZXqrYvtTsE610aUetWRNGuPBnlqRvhBOUWidppUOFml0vzVDGycnaHdLruvQ6V+GRVPPy8Bi/z9miN46a9Xyla/lW6rFTF7OXdwQhVLkHuh+1TgJg0TKOP9RKCmrN2FEDPcs6HoyQnP1h5LEJ+pGTtVsx0QdyHDWOZyvZg7r1fBF+4Eg93xqodRtQo6CV5Bs7vhJjsJ3sPZhJ+Fabk++NSN83apXbZCVt4OgpeV+EK37byh7US/lS1gkeI0hf7+MLxsFooq9LKNK3tew9mCnlq4FRhA/GkL5gb4wuem+CSfsi+9OomfKVxA7CB+2lb6cTjNKlHKvX7c1CY8I3rOz1ujAiwjcWdk76RlgbYTfMlOqDxbCy96CH9M1Qx7fS7Q8DQU3p2/VkItgLw8veg9bSN0NZx6ptFhrCI32lBG1X0cusnPpYOdWDiWTvwSopX0nsJnwP4YjfLOUhs5Ry1Kjbq9nmxI6YTvagT8pXGhG+sniLvyGJm0Xygr0568Qu0gpkBW0NIvza2F343tSVv0i+PjIrx7j03vs9vfd8r+g+TS17EOHXRYRvPcxSylEDu5Zyeu1ZQa89KxTdp+llD9qWdUT4gpWQdH8GNTtOBkN0ZSkZB36k28EfiapUzjeWkL0HEf4ZRPjWQtK9fcg48CNNqspoclr6SmEp2YMI3xsRvhAIMg3TWHjX6pWs21tO9iDC90aEL2iNmadhGqGU412rV1L2hr0t4W0Pbw17H1r94NT4RVTyD0bN5Gbn8oIeqHmSNfotC7W6TaHaYbHPjoX8acVzNCs9HtZ+TsQkM+2qZ/ih/Q01j81/t5Pf7S2Z7D2YOeEriVrpHiThC9phlXT/Q/sbGDnoa1aed0PjG/th5XnV+/AWfWNYWvZgXuGbpZwDInwtUfOTlMzM0Y7C2GReuW4yEwdM5kRMcsCvOxGTzMQBk3nluskUxgb+OrCB7M2MCF8QamOVdO8hmJQfSpr3xhayN2u6BxG+UB+7p3urCb8wNplPL/m/Rrf75NJRQad5b2whey0xev1esAYifGsJv/fexmfdhDszxzayN3P7Ukn3gtkxc/1eC+EH0hoh3PYJtpG9lkg5R4SvBWZK92qg5cpaNYUfXWelrOcibN2Lt+G2T7CV7M2c7tVAhC80hJRzaqOW8D3tEaD2Rdi6F2/DbZ9gK9lriRnSvdqI8NVH7UVtIvzaqCH83nu+9zulsu4UzXDq9iJ7k2Gmcg6I8LXAbKuYrSB8JaV/Kqppo1MqPSm/NKppyMcR2auIWWbmqC18wdyoUb83u/BBuZT/6aWjA5pSWRibzKeXjgr5OCJ7EyLlHKEuZivngHWEb4TmaYEgslcZSffViPDVR4RfjR79780gfVvJ3ug/jGAwW7oXtEGEX01G5ya6Sd+o2Er2QsNIuhcCwSzCB0n53thG9kb85hsREb750WJ2jgi/cTzSN4p7wpZ9bm6uEuNQFb2/2WrV7aWUI/hDhF8bvco6Howg/bBl/+KLLyoxDtXQ+xtsRiTdWwMRfn30vnm5nmk/MtANe/XqRY8ePXC73TgcDgDcbjc5OTmqDS5c7CD6nMPNTdlkqrh9L9MtBjIjcTuyVD+5ppdvUTwgeH6n1fj06hG+3jPlvP2kRSuXgGWfnp7O5MmTSUxMrPX4/fffr/igwsUOklebXdGdTdEMS2gcswofqqWvVrkyo3MT3YXvQQvxB3TD8WPHjrF+/XrS09Np164dkZEBnyNCJtgbjhtd8Gp+fFQr2Wshe0n32qFF+UytEqDa16eMIn1/BHoCaOiG443KfsqUKcyZM4e0tDS2bNlCbGws/fv3Z9SoUaSmpgY34iBoSPZGF7sv1K4VivCFQDCz8EGkX5e6J4GGZN9oRJ8/fz5LliwhJiaG22+/nVdffZVvv/2WIUOGMHnyZLp37x7+iH1gRqELgtHRqqQD6khfzbIOGKeeHyjBeLLR2TgJCQlER0cD4HA46NChAyNHjmTatGm88MILoY/SRug9A8DoyOwcbYnbkWXamTpQLXy1JyXoPVVTDRqV/ZAhQ3jsscc4efJkrcfPP/98jh07FtBBTp06xaOPPsrQoUMZOHAgS5cuDW20guZIR0zrYmbhgza3OrSS9Bst4wwaNIikpCTuvPNO8vLyeOWVV4iNjSUzM5N27doFdJBly5aRkZHB8OHD2b9/Pw888AD9+/cPe/CCdZCpmPpghbIOqF/LN1t5xxcBzcaB6jn1WVlZbNy4kfz8fM477zyuv/56YmJigjpgZmYmb775Jh9//HGD242ZXhrUfo2KpjdVUDHpaDUNU4SvD1qV0tT+pKjVqnKjSv/l4f59HPAcSofDQe/evendu3fIA7nrrrs4duwY77//fsj7EARBeTwnWTOnfFD/Aq4HMyZ9TRuhzZ07l7feeovRo0fjcrm0PLQuWKXWpyVysVZftPpkpXYtX6tV5Z6avhn+1jWR/caNGzlw4AAAXbt2xeVykZ9vviX+dkYu1NoHLYVvFemD8S/maiL79evXM3PmTKB6NW5xcTHNm1u7Y6ORf+hGR9K9/mg1PROsK32jOSDgC7ThUF5eztixYzl48CDl5eX8+c9/5je/+U2DrzHzBVrd+mdr8AstF2rtiVUu4IJ+rcG1qO83dIFWE9mHglllr2vPbAvJHkT4RkPLT1xWlj6oJ35FZuMIjWO0j22CoCRazdgB9WftgHZz9H3h7QqtZvTY5raEamMX0Wt5oVZq98ZEy09catfzQfuafl28a/xqekSSfZjYRfKC4I2WKR+0Tfqgb4lHrdQvsg8DEb36SBsFY2NF6YO+JR5v6jomHPmL7INEBC8I9bG69EF/8YNv/wR6AjDsbJxPV1QPywjLkc0ieC3rjjIrR2gIra+3aHktyQjS98fQqxx+nzN8svcnWjVPAmaRu12QUo75sGrSB+Ol/UAxvOz9IUKujZ6zCQTBH1aWPhinth8IMvVSCAmte+XINExzo2X7BdBmyqY3numbRg5dpk32gmAlijLX+3w8vndPjUeiLnolfdA+7YOxEr/IXjANVqrd+5N7MNuZ+USgtfRB+xIPGEv8InsLYOSPjkJtApV8KPsyo/ztIn3Qv74vshdCZld0Z03romZHSdEHsn8zyd/7E5uVSzygX9oX2QumwqylHLVFH8gxzSJ/O6Z9UF/8InuTIyUc46OH6H3hPQ4ziF9P6YO+4gfl5S+yF0yHmdK9UURfFzOJX48SD+iX9j0onfpF9iZGUr2gBGYq99ipxOONEuIX2QthIRdp/WPUVN8YZkj9divxeBOq+EX2JsXuqd5MpRwzY3Tx27XE4yEY8RtW9nrPSRWEcDBrqm8Is4jfjmkfPM5M9vu8YWXvwUgr0IyC3VO9oD9GFr8e0gfjpH1/GLaffea24w0+b1fxG1X0etXtjVrKsWKyDwSjid+DXo30tBZ/744mTvb+UHtOqiCEil1FD8ZN/JL2TSz7utih3GPUVC8IvjCi+PWWPugnftOWcYLBCvI3uuj1nH5ptFKOnZN9IBhF/B6sVOJpqIxjC9l7Y0bxG130HqRuX43IPjBE+tUoKX1L1uxDxWy1frOIXqhGRB84Rivz6F3iUbu8YzvZ18Wo8hfJC3bCSOLXW/qgjvhtV8YJFj3kb1bRS91ekr2S6C19D3re/zhY6UsZJwy0TP5mlbwgqIFR0r5eSR+ULfFIsleAcE4AVhO83S/SSrJXFyOkfSMnfUn2KmM1YQuCUfGcTCXpB5/0JdkLimLnur2ken3QO+0bKek3lOydag9GEARBTYoy1+t6oo3bkaVb0AgmXInsBUUxQg8QwZ7YVfrp5VsCkr7IXhAES2EE6etBY8I3rOzlVndCsOhZOxWMh97C1/saUl0MK3sI/OOJYCyklCMYBbumfF8YWvYeRPiCIISDntI3ivBNIXsQ4QuCED56Cl9v6ZtG9iDCFwQhfOya8k0lexDhmwW96vZykVYIFLsJ33SyBxG+IAjKYCfhayr7119/nUGDBjFw4EC+/fbbsPYlwheMhLRKMC96lXW0Fr5msl+3bh05OTnMnj2bDz/8kJdeeinsfYrwBUFQCqsLXzPZ9+zZk8mTJwOQkJBARUUFLpdLq8MLOiDz7QWzYeVPaJq1OI6MjCQysvpw8+bNo2/fvjid4Z9r0su3iFSEWhS376X7NDezs3PZpoC3bXdNVxVHYn3idmRpMrFA8372S5YsYc6cOcyYMUPrQwuC0ADBCL6x15n5BFCUuV7ztslaCF/TC7QrVqzgnXfe4YMPPiAxMVHLQws6IZ+6jM/OZZtCFr2W+9QSK5ZzNEv2J0+e5OWXX2bmzJk0b67sfVyllCMIwaOFjD3HMHPStwqaJfsFCxZQUFDAqFGjGDZsGMOGDePAgQNaHV4QBC+0Tt1mTPpap3u1rzMZ9raEedk/BLW9JHtjo/U0Wa0v0JrpY7/e0jVTyte6dh9u3T4lo4/f50y5glYQhNDQW/RGGUOgmOkk3hgie0GwCUaSrJHGYhdE9oIlkYZoxkeEry2WkL3U642P/Iz0xahiNeq4rIglZC8IgiA0jMheEATBBphe9lIeMA/ys9IHKZUIYAHZC4IgqIWW8+zVnlSgeSM0JbFiUsw5fKaVRJfUfB1HIgiClTC17M2Ot9iDeV5OAoJgLSzZ4lgpzJjqG5N7KPsxm/h3RXeWO4wJNRi5dYLWrRLUxpSyN5PolRJ8Y/s3m/QFwcii1xKtFgCa7gKtWUSfc7i56qLX83iCeRCpBo/VUj2YTPZGF71HuHpKV4R/BmmZYFzkBFSNlr+jpinjGFn0RhOslHYEo2IGyWuV6rUOI6ZI9kYVvd4pvjGMOjaj/jytjBEka4QxNIZVRQ8mkL0RxWB0yXtjlnEK1kZEfwa9youGLuMYTfRmFaeUdQSoFq7WrRPMIHmwvujBwLI3iujNKnhf5BxuLsK3OVoJ3yySB3uIHkxQxtETK4negxXfkxAcaoq43TVdRfQ+0EL0jQVkwyZ7vbCDDI2Q8GUlrb54hKxEyjeT3L2xk+hBZF+DHSTvjRGEL+iPt6gDFb9Z5e6N3UQPInvbSV4Q/GEFiQeCVUQf7HVN28peJC/pXrAXVpE8hDaBxXYXaM00R14L5Hsh2AG7ix5slOxFav6RhC9YGauIPtzp6JaXvUjeuMiMHEFNrCJ5UGbdkWVlL5IPDkn3gpWwiuiVXFxqOdmL5ENHhC+YHZG8fywje5G8INgbEX3DmF72InllkXQvmA2RfGCYUvYieEEQtLx1oJqi16rpo6lkL5LXBkn3gtGRNB88ppC9SF4QBJA0Hw6Glr1IXj8k3QtGQiQfPoZtlyCiF8IlbkeW3kMQFEBErwyGTvaCINgXkXzgeMJx747+txHZh0j2lrKgX5PRuYkKI1EPKeUIeiCSD5xgKiAi+wYIRejB7s9sJwBBUBMRfWCEUuYW2Z9GabGHclwRv3mJ792Tosz1eg/DtIjkAyOca5m2lb1ecm8II4pfSjmCmojkA0OJCSu2kr0RBe+P7C1lhhG+ICiNlpIH9URvBsl70FT227ZtY8SIEdx3330MHTpU9eOZSe6+8IxfpC9YBatIHtQVvRpTzzWTfUlJCRMmTODyyy9X9ThmF7wvJOULVsAqojeb5D1oJvvo6GimT5/O9OnTFd+3FQVfFz1TvtTthXAQyTeOFotINZN9ZGQkkZHKHc4OgveFpHzBLFhF8qCe6LXsFGC6C7R2lbw3IvzGkVYJ+qG15MF8aV6PdjCmkL0IXhCMj0g+MPTq+2Vo2Yvk/SPpXjASIvrG0bu5o2ayz87OZuLEiezfv5/IyEgWLVrE1KlTadasme/tRfSNIsIX9EYk3zhaST57SxlcFeP3ec1kn5GRwSeffKLV4QRBc+zUMsFKkgd1RK+p5APA0GUcoXEk3QtaoofkQdK8L4KtfojsLYBZhZ9evkXvIQgBYjXJg3nTfKglbpG9YDlk2qVy6CV5MFeaN7LkPYjsBUHwiaT5wFBb9EpNVhHZWwSzlnIE4yFpPjDMInkPIntBUBAzz8jRU/JgHtGbTfIeRPaCYHNE8oGjpujVXltkO9nnZO0O6XVdep2r8EjsjVozceTibODoLXkwj+jNLHkPlpZ9qGJvbF8ifsHs6C16M12EVUv0WncJsJTslZR7oMcR8Qt1MXLdXm/Jg3nSPFhH9GBy2Wsl90DGYATpqzUjR25cYn6MIHkwj+itJHkPppO9EQTvi5ys3YYQvhmQer12GEXyIKLXu7mjaWRvVMl7o7fwZZ69cdC7lCOSDw01RK+35D0YXvZmkLw3egtfsDdGkjyYR/RWTfPeGFb2ZpO8NyJ8AbRN90aTPNhb9EaSvAfDyl4wBkpfnJV6vbIYUfIgoteaM+G4k99tRPaCoCJqpXujSh5E9FoSTAVEZK8SUsqpj1371yslfCMLHsyzUMoKkofgS90ie4tg5/n1ZijhhCN8o0seJM1rSajXM0X2gqAR3tJuSPxmkLs3InptCHfSishe8IlZLsyaFbMJ3R8iem1QYnaiU4FxCDpj58VUZijhWBURvTYoNQ1dkr2gOpLqrYeIXn2UXmskyV4lzDwTRy7MCg0holcfNRaViuxNjtFLOJLqrYWIXn3U6h4gslcBSfWCFVFzHr2SiOh9I7I3MXZO9VLC0RazLJhSAyuIHkT2iqNVqrfzIipBW8wkeqVTvVVEDyJ74TRmmlcvqV477Cx6qyGyVxAzp3pBqIvdRW+lVA8ie8Uws+gl1Qt1EdFbS/QgslcEmX1zBplqKZgdI954RAlE9mGipejtXr6RVK8Ndk/1WqH13fhE9mFgdtGbKdWL6LXB7qK3aqoHkX3IiOhrI+Ub82Mm0ZsdPe6xLbIPARG9tkiqVx+zrI71IKk+eKTrZRBofSHWLKKX8o3QEJLqjYEk+wAR0ftGyjfmx2zlG0n1oSHJvhH0mFYpoq9GUr36SPnGPkiybwARvX9E9OZHbdFL+cZYiOx90KXXuSL6BhDRC3qgVqrXo4Sjh1+kjOOFXith1VosZUbRC9ogqd5+SLJHvyQPIvq6SKpXH7PV6a2K1s7RLNlPnjyZNWvWUF5ezvjx4+nWrZtWh/aJEfrZiOhrI6K3BmaZgWMEuvQ6V7MFVprIfs2aNWRnZ/P555+zbds2xo8fzz//+U8tDl0LIwge1O1xI6IXGkLKN7UxwpRLrYSviezXrl1L//79AejYsSNHjhzh1KlTxMbGqn5sowge1G9kJqIXGsKs5RurpnpvtBC+JrI/evQonTufOeMnJydz7Ngx2rZtq8rxjCR4D2ZL8yCiF4LDbKneaHi8pZb0NZF9VFRUrX+73W4cDkeDr5n/bic1h2QxklXabx+V9utFhgbHEACIUXn/KSrtt3dHlXYMcJXa35VQUMd9mszGadmyJXl5eTX/Pn78OC1atNDi0IIgCAIayf7qq69m6dKlAGzatIm2bdsSE2PEM6ogCII10aSMk5GRQefOnbnjjjuIiIjghRde0OKwgiAIwmkcbrfbrfcgBEEQBHWRFbSCIAg2QGQvCIJgAwwn+8mTJzN48GAGDhzIxo0b9R6OYmzbto1rr72WTz/9VO+hKMbrr7/OoEGDGDhwIN9++63ewwmbU6dO8eijjzJ06FAGDhxYM6nACpSWlnLttdfyxRdf6D0URVi7di2XXXYZw4YNY9iwYUyYMEHvISnCV199xcCBA7njjjtYvny5ovs2VNdLo7RVUJqSkhImTJjA5ZdfrvdQFGPdunXk5OQwe/ZsTpw4wa233sqNN96o97DCYtmyZWRkZDB8+HD279/PAw88ULPy2+y8++67JCUl6T0MRbnkkkuYMmWK3sNQjOLiYmbMmMGsWbM4ceIEU6dOpV+/fort31Cy17OtgppER0czffp0pk+frvdQFKNnz55MnjwZgISEBCoqKnC5XDidhvuwGDA33XRTzdeHDh0iNTVVx9Eox44dO9ixY4ei4hCUZ+XKlfTt25cmTZqQmprK888/r+j+DfWXefToUZKTz6wG9bRVMDuRkZGWW1cQGRlJXFwcAPPmzaNv376mFr03d911F3/72994+umn9R6KIkyaNIkxY8boPQzFyc3N5cEHH2TIkCGsWrVK7+GEzcGDBzl+/DgPPvgg99xzD6tXr1Z0/4ZK9qG0VRD0ZcmSJcyZM4cZM2boPRTFmDt3Lps2bWL06NHMnz/f1CexL7/8kosuuog2bdroPRRFSU9P5+GHH+amm25i//793HvvvSxatIjo6Gi9hxYy5eXl5OXlMW3aNPbu3ct9993HsmXLFPv9M5Tspa2CuVixYgXvvPMO//jHP0hMTNR7OGGzceNGUlJSaN26NV27dsXlcpGfn09KilpdX9Rn+fLl7Nu3j8WLF3Po0CGio6NJS0ujTx9z9yRKTU3llltuAaBt27a0aNGCw4cPq9ZcUQtatmxJjx49iIiIID09nfj4eEV//wwVWaStgnk4efIkL7/8Mu+//z7Nm1ujBe369euZOXMmAMeOHaO4uNj0723y5MnMmzePOXPmcNdddzFixAjTix7gm2++YerUqUB1KMzLyzP9NZY+ffqwZs0a3G43eXl5iv/+GSrZW7WtQnZ2NhMnTmT//v1ERkayaNEipk6dSrNmzfQeWsgsWLCAgoICRo0aVfPYxIkTad26tY6jCo/BgwczduxY7rnnHsrLyxk3bpypSzhWpl+/fixYsIDBgwfjdrsZN26cqUs4UP1pZcCAAdx7770UFxfz1FNPKfr7J+0SBEEQbIDEFkEQBBsgshcEQbABIntBEAQbILIXBEGwASJ7QRAEGyCyFwRBsAEie0EQBBsgsheERpg9eza33XYb11xzDW+++abewxGEkBDZC0IDLFy4kLVr1zJv3jy+/vpr5s6dy5EjR/QeliAEjaHaJQiCkXC5XLzxxhvMmTOHqKgooqKiSE1N5ddff2XatGnk5ORw9dVX89BDD+k9VEFoFJG9IPghKyuLY8eOce+999Y8lpubS2RkJBEREXz22Wc88sgjHDt2TLqzCoZHyjiC4Ifs7Gzuvvtu5s+fz/z583n11VeJiYlh06ZNXHbZZQBcfPHFbNq0SeeRCkLjiOwFwQ/5+fm1bom5cOFC+vfvT0lJCfHx8QDExcVRWFio1xAFIWBE9oLgh3bt2pGZmQnA9u3b+de//sWoUaNITEykqKgIqL5JtBVu3CJYH5G9IPjh+uuvJzIykmuuuYYnnniCN954g7POOotu3bqxdu1aANatW0fXrl11HqkgNI70sxeEEBg/fjxbt27lyiuvZMSIEXoPRxAaRWQvCIJgA6SMIwiCYANE9oIgCDZAZC8IgmADRPaCIAg2QGQvCIJgA0T2giAINkBkLwiCYANE9oIgCDZAZC8IgmAD/j/Md8eotOqWwAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Discretize the parameter space\n",
    "theta0 = np.linspace(0.0, 2.0 * np.pi, 100)\n",
    "theta1 = np.linspace(0.0, 2.0 * np.pi, 100)\n",
    "#print(cost_fn([theta0[40],theta1[:]]))\n",
    "#parameter_landscape = gen_surface(theta0, theta1, cost_fn)\n",
    "# Load energy value at each point in parameter space\n",
    "#parameter_landscape = np.load(\"param_landscape.npy\")\n",
    "#parameter_landscape_old = np.load(\"param_landscape.npy\")\n",
    "#np.save(\"deut_pes.npy\", parameter_landscape)\n",
    "parameter_landscape = np.load(\"deut_pes.npy\")\n",
    "minloc=np.unravel_index(parameter_landscape.argmin(),parameter_landscape.shape)\n",
    "#parameter_landscape = np.zeros(parameter_landscape_old.shape)\n",
    "#print(parameter_landscape_old.shape)\n",
    "#for i in range(theta0.shape[0]):\n",
    "#    for j in range(theta1.shape[0]):\n",
    "#        parameter_landscape[i,j] = cost_fn([theta0[i],theta1[j]])\n",
    "    \n",
    "# Plot energy landscape\n",
    "fig, axes = plt.subplots(figsize=(6, 6))\n",
    "cmap = plt.cm.get_cmap(\"coolwarm\")\n",
    "contour_plot = plt.contourf(theta0, theta1, parameter_landscape.T, cmap=cmap)\n",
    "plt.xlabel(r\"$\\theta_0$\")\n",
    "plt.ylabel(r\"$\\theta_1$\")\n",
    "plt.plot(theta0[minloc[0]],theta1[minloc[1]],\"r*\",markersize=20,label=\"Minimum\")\n",
    "\n",
    "# Plot optimization path for gradient descent. Plot every 10th point.\n",
    "gd_color = \"g\"\n",
    "plt.plot(\n",
    "    np.array(gd_param_history)[::10, 0],\n",
    "    np.array(gd_param_history)[::10, 1],\n",
    "    \".\",markersize=12,\n",
    "    color=gd_color,\n",
    "    linewidth=2,\n",
    "    label=\"Gradient descent\",\n",
    ")\n",
    "plt.plot(\n",
    "    np.array(gd_param_history)[:, 0],\n",
    "    np.array(gd_param_history)[:, 1],\n",
    "    \"-\",\n",
    "    color=gd_color,\n",
    "    linewidth=1,\n",
    ")\n",
    "# Plot optimization path for quantum natural gradient descent. Plot every 10th point.\n",
    "qngd_color = \"k\"\n",
    "plt.plot(\n",
    "    np.array(qngd_param_history)[::10, 0],\n",
    "    np.array(qngd_param_history)[::10, 1],\n",
    "    \".\",markersize=12,\n",
    "    color=qngd_color,\n",
    "    linewidth=2,\n",
    "    label=\"Quantum natural gradient descent\",\n",
    ")\n",
    "plt.plot(\n",
    "    np.array(qngd_param_history)[:, 0],\n",
    "    np.array(qngd_param_history)[:, 1],\n",
    "    \"-\",\n",
    "    color=qngd_color,\n",
    "    linewidth=1,\n",
    ")\n",
    "\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the red regions indicate states with lower energies, and the blue regions indicate\n",
    "states with higher energies. We can see that the ``QNGOptimizer`` takes a more direct\n",
    "route to the minimum in larger strides compared to the path taken by the ``GradientDescentOptimizer``.\n",
    "\n",
    "Lastly, we can visualize the same optimization paths on the Bloch sphere using routines\n",
    "from `QuTiP <http://qutip.org/>`__. The result should look like the following:\n",
    "\n",
    ".. figure:: /demonstrations/vqe_qng/opt_paths_bloch.png\n",
    "    :width: 50%\n",
    "    :align: center\n",
    "\n",
    "where again the black markers and line indicate the path taken by the ``QNGOptimizer``,\n",
    "and the green markers and line indicate the path taken by the ``GradientDescentOptimizer``.\n",
    "Using this visualization method, we can clearly see how the path using the ``QNGOptimizer`` tightly\n",
    "\"hugs\" the curvature of the Bloch sphere and takes the shorter path.\n",
    "\n",
    "Now, we will move onto a more interesting example: estimating the ground state energy\n",
    "of molecular hydrogen.\n",
    "\n",
    "(2) Hydrogen VQE Example\n",
    "------------------------\n",
    "\n",
    "To construct our system Hamiltonian, we first read the molecular geometry from\n",
    "the external file :download:`h2.xyz </demonstrations/h2.xyz>` using the\n",
    ":func:`~.pennylane_qchem.qchem.read_structure` function (see more details in the\n",
    ":doc:`tutorial_quantum_chemistry` tutorial). The molecular Hamiltonian is then\n",
    "built using the :func:`~.pennylane_qchem.qchem.molecular_hamiltonian` function.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(15.531709) [I0 I1]\n",
      "+ (-2.1433) [X0 X1]\n",
      "+ (-2.1433) [Y0 Y1]\n",
      "+ (0.21829) [Z0]\n",
      "+ (-6.125) [Z1]\n",
      "+ (-9.625) [Z2]\n",
      "+ (-3.91) [X1 X2]\n",
      "+ (-3.91) [Y1 Y2]\n"
     ]
    }
   ],
   "source": [
    "#geo_file = \"h2.xyz\"\n",
    "\n",
    "#symbols, coordinates = qml.qchem.read_structure(geo_file)\n",
    "#hamiltonian, qubits = qml.qchem.molecular_hamiltonian(symbols, coordinates)\n",
    "\n",
    "#15.531709 - 2.1433 X0X1 - \"\n",
    "#        \"2.1433 Y0Y1 + .21829 Z0 - 6.125 Z1 - \"\n",
    "#        \"9.625 Z2 - 3.91 X1 X2 - 3.91 Y1 Y2\n",
    "        \n",
    "coeffs = [15.531709 ,-2.1433, -2.1433, .21829, -6.125, -9.625, -3.91, -3.91]\n",
    "obs = [qml.Identity(0)@qml.Identity(1),qml.PauliX(0) @ qml.PauliX(1), qml.PauliY(0) @ qml.PauliY(1),qml.PauliZ(0),qml.PauliZ(1),qml.PauliZ(2),qml.PauliX(1) @ qml.PauliX(2),qml.PauliY(1) @ qml.PauliY(2)]\n",
    "#coeffs = [-2.1433, -2.1433, .21829, -6.125, -9.625, -3.91, -3.91]\n",
    "#obs = [qml.PauliX(0) @ qml.PauliX(1), qml.PauliY(0) @ qml.PauliY(1),qml.PauliZ(0),qml.PauliZ(1),qml.PauliZ(2),qml.PauliX(1) @ qml.PauliX(2),qml.PauliY(1) @ qml.PauliY(2)]\n",
    "H = qml.Hamiltonian(coeffs, obs)\n",
    "hamiltonian=H\n",
    "print(H)\n",
    "qubits=3\n",
    "#print(\"Number of qubits = \", qubits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For our ansatz, we use the circuit from the\n",
    "`VQE tutorial <https://pennylane.ai/qml/demos/tutorial_vqe.html>`__\n",
    "but expand out the arbitrary single-qubit rotations to elementary\n",
    "gates (RZ-RY-RZ).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "dev = qml.device(\"default.qubit\", wires=qubits)\n",
    "hf_state = np.array([1, 1, 0, 0], requires_grad=False)\n",
    "'''\n",
    "def ansatz(params, wires=[0, 1, 2]):\n",
    "    #qml.BasisState(hf_state, wires=wires)\n",
    "    for i in wires:\n",
    "        qml.RZ(params[3 * i], wires=i)\n",
    "        qml.RY(params[3 * i + 1], wires=i)\n",
    "        #qml.RZ(params[3 * i + 2], wires=i)\n",
    "    qml.CNOT(wires=[2, 1])\n",
    "    qml.CNOT(wires=[2, 0])\n",
    "    #qml.CNOT(wires=[3, 1])\n",
    "'''\n",
    "def ansatz(params, wires=[0,1,2]):\n",
    "    #for i in wires:\n",
    "    qml.PauliX(0)\n",
    "    qml.RY(params[0],wires=1)\n",
    "    qml.RY(params[1],wires=2)\n",
    "    qml.CNOT(wires=[2,0])\n",
    "    qml.CNOT(wires=[0,1])\n",
    "    qml.RY(-params[0],wires=1)\n",
    "    qml.CNOT(wires=[0,1])\n",
    "    qml.CNOT(wires=[1,0])\n",
    "\n",
    "#print(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the qubit register has been initialized to $|1100\\rangle$, which encodes for\n",
    "the Hartree-Fock state of the hydrogen molecule described in the minimal basis.\n",
    "Again, we define the cost function using the ``ExpvalCost`` class.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "cost = qml.ExpvalCost(ansatz, hamiltonian, dev, diff_method=\"parameter-shift\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this problem, we can compute the exact value of the\n",
    "ground state energy via exact diagonalization. We provide the value below.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "exact_value = -2.22"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now set up our optimizations runs.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "init_params = np.random.uniform(low=0, high=2 * np.pi, size=3)#(qubits-1))\n",
    "max_iterations = 500\n",
    "step_size = 0.01\n",
    "conv_tol = 1e-09"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As was done with our previous VQE example, we run the standard gradient descent\n",
    "optimizer.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration = 0,  Energy = 7.08874905 Ha, C = 1.6357e+00\n",
      "Iteration = 20,  Energy = -1.72667384 Ha, C = 5.9942e-02\n",
      "Iteration = 40,  Energy = -2.03551209 Ha, C = 1.8440e-03\n",
      "Iteration = 60,  Energy = -2.04482410 Ha, C = 5.4870e-05\n",
      "Iteration = 80,  Energy = -2.04510127 Ha, C = 1.6346e-06\n",
      "Iteration = 100,  Energy = -2.04510953 Ha, C = 4.8713e-08\n",
      "Iteration = 120,  Energy = -2.04510978 Ha, C = 1.4519e-09\n",
      "\n",
      "Final convergence parameter = 0.00000000 Ha\n",
      "Number of iterations =  123\n",
      "Final value of the ground-state energy = -2.04510978 Ha\n",
      "Accuracy with respect to the FCI energy: 0.17489022 Ha (0.17489022 kcal/mol)\n",
      "\n",
      "Final circuit parameters = \n",
      " [3.50058904 6.02654419 3.78727399]\n"
     ]
    }
   ],
   "source": [
    "opt = qml.GradientDescentOptimizer(step_size)\n",
    "\n",
    "params = init_params\n",
    "\n",
    "gd_cost = []\n",
    "\n",
    "for n in range(max_iterations):\n",
    "    params, prev_energy = opt.step_and_cost(cost, params)\n",
    "    gd_cost.append(prev_energy)\n",
    "\n",
    "    energy = cost(params)\n",
    "    conv = np.abs(energy - prev_energy)\n",
    "\n",
    "    if n % 20 == 0:\n",
    "        print(\n",
    "            \"Iteration = {:},  Energy = {:.8f} Ha, C = {:.4e}\".format(n, energy, conv)\n",
    "        )\n",
    "\n",
    "    if conv <= conv_tol:\n",
    "        break\n",
    "\n",
    "\n",
    "print()\n",
    "print(\"Final convergence parameter = {:.8f} Ha\".format(conv))\n",
    "print(\"Number of iterations = \", n)\n",
    "print(\"Final value of the ground-state energy = {:.8f} Ha\".format(energy))\n",
    "print(\n",
    "    \"Accuracy with respect to the FCI energy: {:.8f} Ha ({:.8f} kcal/mol)\".format(\n",
    "        np.abs(energy - exact_value), np.abs(energy - exact_value) \n",
    "    )\n",
    ")\n",
    "print()\n",
    "print(\"Final circuit parameters = \\n\", params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we run the optimizer employing quantum natural gradients.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration = 0,  Energy = 3.07502059 Ha\n",
      "Iteration = 20,  Energy = -2.04510918 Ha\n",
      "\n",
      "Final convergence parameter = 0.00000000 Ha\n",
      "Number of iterations =  29\n",
      "Final value of the ground-state energy = -2.04510979 Ha\n",
      "Accuracy with respect to the FCI energy: 0.17489021 Ha (109.74413458 kcal/mol)\n",
      "\n",
      "Final circuit parameters = \n",
      " [3.50058187 6.02656594 3.78727399]\n"
     ]
    }
   ],
   "source": [
    "opt = qml.QNGOptimizer(step_size, lam=0.001, diag_approx=False)\n",
    "\n",
    "params = init_params\n",
    "prev_energy = cost(params)\n",
    "qngd_cost = []\n",
    "\n",
    "for n in range(max_iterations):\n",
    "    params, prev_energy = opt.step_and_cost(cost, params)\n",
    "    qngd_cost.append(prev_energy)\n",
    "\n",
    "    energy = cost(params)\n",
    "    conv = np.abs(energy - prev_energy)\n",
    "\n",
    "    if n % 20 == 0:\n",
    "        print(\n",
    "            \"Iteration = {:},  Energy = {:.8f} Ha\".format(n, energy)\n",
    "        )\n",
    "\n",
    "    if conv <= conv_tol:\n",
    "        break\n",
    "\n",
    "\n",
    "print(\"\\nFinal convergence parameter = {:.8f} Ha\".format(conv))\n",
    "print(\"Number of iterations = \", n)\n",
    "print(\"Final value of the ground-state energy = {:.8f} Ha\".format(energy))\n",
    "print(\n",
    "    \"Accuracy with respect to the FCI energy: {:.8f} Ha ({:.8f} kcal/mol)\".format(\n",
    "        np.abs(energy - exact_value), np.abs(energy - exact_value) * 627.503\n",
    "    )\n",
    ")\n",
    "print()\n",
    "print(\"Final circuit parameters = \\n\", params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualizing the results\n",
    "^^^^^^^^^^^^^^^^^^^^^^^\n",
    "\n",
    "To evaluate the performance of our two optimizers, we can compare: (a) the\n",
    "number of steps it takes to reach our ground state estimate and (b) the quality of our ground\n",
    "state estimate by comparing the final optimization energy to the exact value.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAFYCAYAAACoFn5YAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABGQUlEQVR4nO3dd3xUVf7G8c+UFAiQRqhJ6AoKUhWwAVJUEKWjFPdn3UWxoQgqKiAqIiqKii66FhBBKYIgIiBVxAKKgEjvNZ30NvP7Y5iQQEIoM3NnMs97X3mtU3LvN4chD+fcc88x2e12OyIiIuKTzEYXICIiIhdPQS4iIuLDFOQiIiI+TEEuIiLiwxTkIiIiPkxBLiIi4sOsRhdwMeLiUl16vPDw8iQlZbj0mL5I7eCgdnBQOzioHRzUDg5GtUNUVMUSX1OPHLBaLUaX4BXUDg5qBwe1g4PawUHt4OCN7aAgFxER8WEKchERER+mIBcREfFhCnIREREfpiAXERHxYQpyERERH6YgFxER8WE+uSCMiIgRDh06yDvvvElSUiJgp3r1mgwbNoKwsLALPta7706ibt16REREcvToEXr27HNB37927Spat76WgICAYl/PyMjg7rv7M3v2txdc26VYsWIZHTp08ug5/Z165CIi5yE/P59Ro0YwYMBgpk79jKlTP+fyyxvyzjtvFHmfzWa7oOO2aXPtBYc4wMyZX5Cbm3vB3+dOubm5zJo1w+gy/I565CIi5+HXX9dTu3YdmjVrUfDcXXcNxm63A/Dyy6OxWq0kJyczatRoxowZRUZGBtnZ2TzxxHCuuKIxS5Z8x4wZnxMVVYXAwCDq1q3Hd999y549uxk69HHmzPmK5ct/wG630aFDJ/r1G8BHH31ARkY6+/fv5/Dhgzz++HCSk5P4++8tPPXUo7z99pSCXnl6ehqjRo0gOzubpk2bF9S5adOf/Pe/72E2m6latRojRz5PfHw8Y8eOwmQykZeXx/PPj6Vq1Wq89to4jh07itls5vnnxxIeHsEbb4znxImjZGRk8sADD9GiRSsefvgBrrmmDb///isnT55kwoS3mD79M3bv3sXEieN56qmRnv0D8mN+H+Q7D+1g9Mf/44V7XsFs1gCFiC8YvW4U3+7+xqXH7F6vB6OvHVfi6wcP7qdevfpFnjvzd0ZoaBgjRozi4MED3HxzVzp27MLGjb/zxRefMW7cBP773/f5+OPpVKhQgfvuG1zke48dO8rq1St4772pAAwZch/t23fEZDJx4sRx3njjHX7++Sfmz5/Lq69O5KOPPmDixHeKDK0vWbKYevUaMHTo4/z44zKWLv0egLffnsikSe9TqVIl3nvvbVasWEZCQjxXX92ae+55gG3btpKQEM/mzZsIDQ1j1KgxLFu2hDVrVhESEkJERCQTJ77Gzp0HeeyxIXz22ZeYTCZCQkKYPPlD3n//HVatWsGAAYNP/QNDIe5Jfh/kz709gpWfLadmzVo8cMsQo8sRES+Vl5dPfn5+weORI4eRlpZGXNwJPvtsJgCNGl0BQGhoKGvXrmbu3K/Jzs4mODiYlJQUypUrX3A9vUmTpkWOv337Pxw4sJ9HHvk3ABkZ6Rw9egSAq65qBkDVqlVJSyt506h9+/bQrFlLAJo3d4wcnDx5koMH9/Pss08BkJWVRXh4ONdc05Znn32KtLQ02rXrQNOmzVmxYllBT75Tp5sBmDRpIhs2/Mq2bZvJyckjOzu7YEjf+d6qVauSkpJyoU0qLuL3QV47vDYAy7YuVZCL+IjR1447Z+/ZHerWrcf8+XMLHo8f/yYAffp0L7gubrU6esdffz2TqKgoXnxxHFu3bmHKlHew2+2YTIWPaD/rHG3aXMuIEaOKPPf7779isZzeqMM5lF8cu52CcxR+X2RkZd59979nvf/TT7/kl1/WMXnym3TrdjsANtvZxx88+B4GDux31s6T51uXuJffjyU3jrkKgN/3/YrNfmGTVETEf7RqdQ1Hjhxm3bq1Bc9t2fIXGRkZWCxFf5WePJlCzZoxgGMWd15eHqGhoaSnp5OamkpeXh5//fVnke+5/PJGbNz4O1lZWdjtdiZNmkh2dlaJ9ZhM5iIjBACxsbXYvv0fADZs+A2ASpUqYTKZ2LNnFwCzZ89k166dLFu2hAMH9tGhQyfuuututm37m8svb8Sff24E4Kef1vDZZx9zxRWNWbNmFQBJSYl8+OF7F1STuJ9X9MhPnDjByy+/zPXXX0/fvn09eu6I8EgAUlNOsvH477Sqdo1Hzy8ivsFqtTJx4tu8/fZEPv30I/LycgkPj2TChEkEBQUXee+tt3Zn3LgXWLXqR3r16seaNStZsuQ77r33QYYOfYDq1WtQp049bDZbwXX2atWq0b//QIYOfRCTCW64of1Zxy2sefMWDB36IG+/PaVguP6WW7rx7LNP8dhjQ7jqqmYFIwUjR77Aq6++hNVqpUqVKtx+ey/y8/N5883XCAwMJC8vj2HDRlC7dh1+//1Xhg59EIDnnx9LZGRlNm78jTvvvJOsrGzuvfffJdZUuXJl8vJyeeGFZxg79tWLb2y5ICa7G8dDduzYwUMPPcT//d//MWjQIAAmTZrE+vXrycnJYcyYMTRp0oT4+Hh27NjB4cOHzyvIzxzeuRQ//bSGnj27QTt4dNgwRrUd7bJj+5qoqIoubVtfpXZwUDs4qB0c1A4ORrVDVFTFEl9z29B6RkYGL730Em3bti14bv369WzZsoWZM2cyfvx4xo8fDzj+FVf4WosnhYaGAWDJtvD9vkWG1CAiInKx3Da0HhgYyNSpU5k6dWrBc7/88gsdO3YE4LLLLuPEiRNkZmZSrly5Czp2eHh5rFbXBH/9+o7rWNUt1dmRtJ0k81Eui7zMJcf2Ref6V58/UTs4qB0c1A4OagcHb2sHtwW51WrFai16+Li4OBo2bFjwOCIigvj4eA4dOsSXX35JamoqYWFhdO7c+ZzHTkrKcFmd+fmOWaYV80KBQ3yx4SuGNn/MZcf3JRo6c1A7OKgdHNQODmoHB28cWvfoZLcz1wR23I5hom3btkWG4D2pfPnyBAQEEJATgNlkZvHehX4b5CIi4ns8evtZVFQUCQkJBY8TExOpXLmyJ0s4i8lkIiIigrSUVK6p1obfj/3KiYwThtYkIiJyvjwa5DfeeCPLly8HYOvWrcTExBAcXPLtFZ4SERFBSkoyt9a5DTt2fti32OiSREREzovbgnzLli0MHjyYefPm8fnnnzN48GCio6Np2LAhPXv25MUXX2TkSO9Yjzc8PJykpCQ613IsSfjD/u8NrkhEvNHhw4d4+uknuP/+u7n//ruZOPFVMjMzXXb8P//ceGqLVM9bsWLZJR+jT5/uZGS4bg6T008/reHll0eTkBDPhAkvX/D379q1kwMH9p/zPffdN7hgSVxPWbt2lUt2sHPbNfLGjRszbdq0s54fPny4u0550SIiIrDZbFSxVKFeWH1WH1xJdn42QZYgo0sTES9hs9kYNepphg59gpYtrwZgxoxpvPrqWJctfrJo0QLuumsQ4eERLjnehZg+/bPz2ke88CI2ruCcK3U+IiMr8/TTz13wOVat+pGGDa8gNrbWBX+vO82c+QUtWlxd4p7y58srVnYzWkSE4y9NcnIynWrdzIeb3uPnIz/RPuYmgysTEW/x22+/EBtbqyDEAe66axB33dWbxMQEpkyZTPv2Hbnuuhv46ac1rFy5nOeeG83MmdNZsWI5NpuNtm2v4957H2Thwvls2fIXSUmJ7N+/n0GD/kXVqlVZs2Yle/fuYdy4Cdx330AWLXJcihw16ml69eqH3W5n9uxZWCwWtm//h/vv/zcrV/7Inj27GDFiFC1atCqo7fDhQ7z66lhq1KjJzp3badjwSkaMeO7UNqOvYLUGYDabeeml8SxcOJ9du3bw7LPD6dOnP3PnfsW4cRMA6NatI4sWLWfkyJHk5dnPuU1rcX74YTHTp39KtWrVqVQplGbNWmA2m/n5559ISIhn9OiX+eqrL9m69S9ycnLp0aM33bv3YPfuXYwb9wKVKoVSo0Y0AEePHmHUqBF8/PG0Yrdm3bTpD+bN+xqTyczevXvo2LEzN97Ygfnz57Jq1Y+Eh4cXqXPSpNfZsmUzsbG1yMtz9Izj4+N47bVx5OTkYLFYGDFiFJGRlRk79nni4+PIy8th8OB7ueGG9nz55XRWr15BXl4eQ4Y8QosWrS55K9qLoSCncJAn0Sm2Cx9ueo9l+5coyEW81OjRo/j2229ceszu3XswenTJG7Hs37+P+vWLrjFhMpmoW7ce+/btLfH77HZ4550pBAYG0a9fD/r3H4DZbGb37l188MH/OHjwAGPGPMcnn8ygfv3LGDbsaapVq1bssUwmE7t27WTGjNmsWvUjH330AdOnf83KlctZsuS7IkFuNpv555+/GTPmFUJDw+jb93ZSU1NJSkrk4YefoHHjJnz88Yf88MNiBgy4my+++IxXXnmdjRt/L/FnOdc2rS+//HoxP7udDz54l48/nk758uX517/uLKgxLu4EU6Z8TG5uLpGRlXn//Y/Jycmmf/+edO/eg08//Yj77x/CddfdwBtvvFYQtE7Fbc0aGVmZv//eyowZc8jPz6N//57cc88DtG7dlvbtOxYJ8b1797B162amTv2M+Pg4+vZ1bBrz0Ucf0L//QFq1uoZ169by2Wcfc8cdvUlOTuL99z/Cas3j22+/58iRw6xYsYwPPvgfR44cZtq0T6hRo+Ylb0V7MRTkOK6RAyQlJdH2yusICajA0v1LGHf9awZXJiLeIjc3p9gNQex2GxZLyb9KrVYrTzwxFLPZTHJyEidPngSgceOrsFgsVKlSlbS0tPOuo379+gQEBBAREUmdOnWxWq1ERESSnp5+1nujo2OJjHTcGeR4TxqhoaFMmTKZnJwc4uPj6Nz5lvM+97m2aS1OSkoK5cuHFPyOdW7HCtCwYSNMJhOBgYGcPJnCQw/dj8ViKZgjsG/fnoLgbdasBevX/1TwvSVtzRoZWZnLL29YUI/9HBth7du3h4YNr8RkMhEVVYXoaMfiYNu2/c3+/fv49NOPsNlshIeHExtbi7S0VF566Xm6dbuVzp1v4aefVtOwYSPMZjPR0TE888wLrFq14pK3or0YCnJO98hTUpIJtATSLroD3+39lt3JO6kX1sDg6kTkTKNHjztn79kdatWqw6JF84s8Z7PZ2Lt3L7GxtYpc53UG/uHDh5g9eyaffPIF5cuHMHBgn4L3XMgWoIX/AVH4+4oubX32Mc5c+tput/P2228wcODdtG17PdOnf0p2dnaR95x5vbrwuc+1TWtxzty6tfCxncfauPF3/vxzA5Mn/xez2UyXLjee+t7C7z/7Zytua9aNG38/7+W+C2/56qzVaezYV4mKqlLk/VOnfs6mTX+wZMm3zJ//LTfd1LnYLV8vdSvai+H325jC6SBPSkoCKJi9vmz/D4bVJCLe5eqrW7N7964i25jOnPkFV155JeHh4YSEVODkyRQAtm/fBkBqairh4RGULx/Cli2bOXHixDlnKZvNp7cBtVis5OTkkJeXx65du1z2czi3WM3OzmbdurXk5eUBp/chr1ChAikpjp/j8OFDxY4WFLdNa3FCQ0NJTU3l5MmTZGdnFWyReuaxqlatjtVqZdWqH7HZ7OTm5p7aktXRjs4tWZ1K2pq1JCaTqdgtX3fs2I7dbuf48WMcOnQQoMi2rRs2/MbSpd+zffs/LF++lBYtWjFq1Ci2bfubBg0u5++/N5Ofn09iYgLPPPOkS7aivRjqkXN6aD052RHknWp1AWDp/h/4d9OHDatLRLxHUFAQr732FhMnvsrUqe9js9m44orGBbOob765K6+99hJ//rmRkJAQbDYbDRpcRkhIBR566H4aN25Cr159mTTpdTp1urnYczRr1oIXXniGl19+nTvu6MUTTzxM7dp1qFOnzgXN7j6Xvn3v4rnnhlOzZjR33TWISZMm0rFjFy677HL+/e97mDLlYwICAnnuueFER8dSrVr1gu1QnYrbpnXx4oVnnctsNvN//3cfDz10P7Vr16FRoyuxWCxFwuvqq1szY8Y0Hn30P7Rtez033NCON9+cwL/+dR+vvDKG2bNnUqNGzbNGDorbmnXLlr+K/ZmbNm3O5MlvUr58eVq1cmxVXb9+A+rWrce//30PsbG1uOyyhthsNu6770Fefnk0y5f/gMlk4tlnX6RixUr897/vs2DBXOz2fO6//z/UqFGTzp1v5aGH7sdms/Hggw+5ZCvai+HWbUzdxdXr3O7Z8zdt2rTh4Ycf48UXXwKg09c3si1hK9vv3UeFQO9aIN9dtJayg9rBQe3gUFw7bN26hUmTXufDDz9x6a1Y3uxiPw8rViyjZcurqVQplGHDhnLvvQ/SuPFVbqjQM7xxrXX/+ASWovCsdadOsZ3JteWy6tBKg6oSEW915ZWNueqqptx77yCXLKRSlmVlZfHoo0P4z3/uJTo6xqdD3FtpaJ2is9adOtW6mTc3vM7y/T/QrW53o0oTES/1yCPDjC7BJ9x6623ceuttRpdRpqlHDgXXJgr3yJtXaUlkcCTLDvzg8hmGIiIirqIgx3GfZ6VKoSQnJxc8ZzFb6BDbiWPpR9masMW44kRERM5BQX5KWFhYkR45QMfYzgD8eGCpESWJiIiUSkF+SlhY+FlB3j6mIyZMup9cRES8loL8lLCw8IINAJwiy0XSomorfjv2CynZycYVJyIiUgIF+SmnJ7wlF3m+Y2xn8u35rNZtaCIi4oUU5KeEhRVd3c3JeZ1cw+siIuKNFOSnnF6mNbnI802rNKdyucr8eGCZbkMTERGvoyA/paQeudlkpkNMJ45nHGNLwmYjShMRESmRgvwU5zVy5164hXWs5RheX67hdRER8TIK8lOcPfKUlOSzXmsfcxNmk5nlup9cRES8jIL8lOLWW3eKCI6kRZVW/H7sV5Kzzn5dRETEKAryU0JDw4Czr5E7darVhXx7PqsOrfBgVSIiIuemID/l9Kz1koMcdBuaiIh4FwX5KadnrScX+3rjyldRpXxVlh9Yis1u82BlIiIiJVOQn1KuXDkCAwNL7JGbTWY6xnYmPjOOv+L+9GxxIiIiJVCQn2IymQgLCy92spuThtdFRMTbKMgLCQ8PL/b2M6d20R2wmq0sP6AgFxER76AgLyQ0NIzk5GRstuKvgVcKCqV1tbZsPL6B+Mx4D1cnIiJyNgV5IeHh4dhsNlJTT5b4no61umDHzooDyzxYmYiISPEU5IWUNnMdTl8n1/C6iIh4AwV5Iaf3JC95wtvl4Q2JrhDDjweWkWfL81BlIiIixVOQF+LskZ9r5rrJZKJjrS4kZyez4fjvnipNRESkWAryQs61cUphzuH1pfu+d3dJIiIi56QgL+RcG6cUdkPNdgRbglm6f4knyhIRESmRgryQ87lGDlA+oDw3RLdjW+JWDqYe8EBlIiIixVOQFxIeHgFAYmJiqe/tXOsWAH7Q8LqIiBhIQV5IREQkAElJ5xPkNwPww77Fbq1JRETkXBTkhURGOoI8MTGh1PfWrBjNlZFN+OnwGtJy09xdmoiISLEU5IWEhFQgICDgvIIc4Obat5Bjy2H1wZVurUtERKQkCvJCTCYTERGRJCScX5B3ru28Tq7hdRERMYaC/AwREZHnNdkNoHmVllQuF8XS/Uuw2YvfaEVERMSdFORniIyM5OTJFHJzc0t9r9lkpnOtm4nLPMGmE394oDoREZGiFORncM5cP99eufM2tCX7NbwuIiKepyA/Q0SE817y87tO3j6mA4HmQJbsVZCLiIjnKcjPcLpHfn5BXiGwItdH38jWhM0cOLnfnaWJiIicRUF+htP3kp/f0DrALbW7AbBk33duqUlERKQkCvIzXGiPHOCWOl0B+H6vglxERDxLQX6G0+utn3+QVwupTvMqLVh3ZC3JWefecEVERMSVFORnuJBlWgu7pXY38u35LDvwgzvKEhERKZaC/AzOofXzXd3N6ZY6juvkGl4XERFPUpCf4WKukQM0jGhErUq1WX5gKdn52e4oTURE5CwK8jOUL1+e4ODgCw5yk8nErXVuIz03jZ8Or3ZTdSIiIkUpyM/g3DjlQm4/c7r11PD6Yg2vi4iIhyjIi3EhO6AVdnW11kQER7Bk33faREVERDxCQV6MiIhI0tPTyM6+sGvdVrOVLrVv5Vj6Uf44scFN1YmIiJymIC9GZKTjXvKkpAsfXu9W93YAFu351qU1iYiIFEdBXoyLvQUNoF10B8pbQ1i0ZwF2u93VpYmIiBShIC/Gxazu5hRsDaZzrZvZm7KHbYl/u7o0ERGRIhTkxbjY1d2cutXtDsCiPQtcVpOIiEhxFOTFuJShdYBOtboQaA7UdXIREXE7BXkxLnZ1N6cKgRVpH3MTfydsYU/KbleWJiIiUoSCvBiXGuRwevb6d3sWuqQmERGR4ijIi3Gp18gButS+FYvJouvkIiLiVgryYjhnrV/sNXKAyHKRXFvjejYc/42jaUdcVZqIiEgRCvJilCtXjvLlQ0hKSrqk43TV7HUREXEzBXkJIiMjL2loHeC2urdjwsS3e+a7qCoREZGiFOQlCA+PuOQgrxpSjdbV27L+yDqOpx9zUWUiIiKnKchLEBERQUZGBhkZGZd0nO717sCOnUV7dU+5iIi4noK8BM5b0C5m45TCbqt7BwDf7vrmUksSERE5i4K8BK64BQ2geoUaXF2tNT8f/YkTGSdcUZqIiEgBBXkJLnWZ1sK617sDm93G4r1aHEZERFxLQV4CV6zu5uQcXl+w+5tLPpaIiEhhCvISuGpoHSC6Ygwtq7Zi3eE1xGfGX/LxREREnBTkJXDl0DpA93o9ybfna3hdRERcSkFeAlfNWne67dQmKvN3zXPJ8UREREBBXqKICMd6664YWgeIrVSLFlVa8tPh1RpeFxERl1GQl+D0ximu6ZED3FG/N/n2fL7bo8VhRETENRTkJQgKCqJChYokJLiu99y9nmP2+vzdGl4XERHXUJCfgys2TiksumIMrapew0+HVxOXEeey44qIiP9SkJ9DZGQkCQnx2O12lx3zjvo9sdltLNSOaCIi4gIK8nOIjKxMbm4uaWmpLjtm93o9AFig2esiIuICCvJzcN6CFh/vuuvkNSrU5JpqbVh3ZC3HM4677LgiIuKfSg3yxMRExo0bx+OPPw7A0qVLiYvzj+u7kZGVAdfdguZ0R/2e2LGzcLeG10VE5NKUGuSjRo2iTp06HD9+uvc4cuRItxblLZxB7sqZ6+AYXjdhYv6uuS49roiI+J9Sgzw9PZ2BAwditVoB6Ny5Mzk5OW4vzBs411t31TKtTtVCqtOmxrX8cvRnjqUfdemxRUTEv5Qa5Lm5ueTm5mIymQBHqGVmZrq9MG9wukfu2iAHuKN+L+zYNelNREQuSalBPnDgQPr06cOuXbv4z3/+w+233869997ridoM51ym1dVD6+DY2tRsMjNv1xyXH1tERPyHtbQ3dOvWjVatWrF582ZMJhNjx46lSpUqnqjNcO6a7AZQpXwVrqt5I2sOreRg6gFiKsa6/BwiIlL2ldoj3759O9OnT6dTp0507NiRN998k+3bt3uiNsNVruyeyW5OPer3ArQjmoiIXLxSg3zMmDG0a9eu4HGfPn0YM2aMW4vyFhUrViIgIMBtQd6tbnesZivfaHhdREQuUqlBbrFYaNWqVcHjVq1aYbFY3FqUtzCZTERERLplshtARHAk7aI78Ffcn+xJ2e2Wc4iISNlWapBXrFiRL774gp07d7J9+3Y+/vhjQkJCPFGbV4iMrOy2IAfH7HWA+Tt1T7mIiFy4UoN83Lhx7Ny5k+HDhzNixAgOHjzIK6+84onavEJkZCSpqSfddu981zq3EWgO1PC6iIhclFJnrUdERDB69GgPlOKdnIvCJCYmUK1adZcfv1JQKDfV6sz3exfxT+I2GkY0cvk5RESk7Cq1R75gwQJ69OhBhw4daN++fcGXv3DHxiln6lm/N4B65SIicsFK7ZG/++67vPLKK1Sv7vreqC9w573kTl1q30p5a3m+2TmHEVc/V7CKnoiISGlKDfLatWsXmbXub9y1cUphIQEhdKl9C9/smsvm+E1cFdXMbecSEZGypdQgb9GiBRMnTqR169YFG6cAtG3b1q2FeYvC18jdqUf9Pnyzay7zds5RkIuIyHkrNcjXrFmD2Wxm06ZNBc+ZTCY/CnJHj9yd18gBbortRMXASnyzaw7Ptx2D2VTq9AUREZHSg/yLL74AIC8vr0iP3F94YmgdINgaTNc6tzFr+wx+O/Yrrau3cev5RESkbCi127d582a6d+9Ot27dAHjnnXfYsGGD2wvzFs5Z64mJiW4/V88GfQD4Ztdst59LRETKhlKDfMKECbz11lsFO5716NGDCRMmuL0wb+HOrUzPdEPNdkQGRzJ/1zzybHluP5+IiPi+81prvX79+gWPY2NjCQgIcGtR3iQgIIDQ0DC3T3YDCLAEcFu9HsRnxrHuyFq3n09ERHzfec2oOnz4cMG9zatXr8Zut7u1KG8TGRnp9sluTs7FYebt1PC6iIiUrtQgHzFiBEOGDGHjxo20bNmSiRMn8vzzz3uiNq8RGVmZxMQEbDab28/VunpbqoVUZ+GeBWTnZ7v9fCIi4tvOa631BQsWkJqaislkokKFCp6oy6tERkaSn5/PyZMphIWFu/VcFrOFO+r34sNN77Hy4I/cXPtWt55PRER8W6k98ieffBJwbGfqjyEOnrsFzalXfcfs9Xk7v/bI+URExHeV2iOvV68ew4cPp3nz5gQGBhY836dPH7cW5k1OB3ki9eq5/3zNqrSgdqU6fL/3O9Jz0wkJ8J/930VE5MKU2iPPycnBarWyefNmNmzYUPDlT5z3knuqR24ymejVoA8ZeRn8sG+xR84pIiK+qdQe+auvvkpeXh5xcXF+vAOaZ4McoEeDPry54XXm7ZpTsFCMiIjImUrtka9evZrOnTtz7733AjBu3DiWLl3q9sK8iac2TimsYUQjGkVcyY/7l5KSneyx84qIiG8pNcinTJnCnDlzqFzZcZ146NChfPjhh24vzJt4auOUM/Vq0IccWw6L9nzr0fOKiIjvKDXIg4KCCpYpBQgLCyM4ONitRXkbZ5B7skcO0KOBY3GYuVocRkRESlBqkAcGBvL7778DkJKSwpdffklQUJDbC/Mmnp7s5lSrUm1aVr2atYdXcTzjuEfPLSIivqHUIH/hhReYOnUq//zzD126dGH16tWMHTvWE7V5jZCQEIKDgz0e5OAYXrfZbXy7a57Hzy0iIt6vxCCfO3cuAD///DMffvghv/32G7/88gtTpkyhZs2aHivQG5hMplPLtLp/K9Mz3V6/F2aTmTlaHEZERIpR4u1n7777Lrm5uXz66acFG6YU5k8LwoBjeH337l0eP2/V8lW5vmY7Vh9awb6UvdQOrePxGkRExHuV2CMfOXIkf/31F6mpqUUWgvHHBWHAcQtaRkY6mZmZHj937wZ9Afhm1xyPn1tERLxbiT3yvXv38vLLL/Pee+/x8MMPe7Imr+Sc8JaYmEDNmtEePXfXurcxfNXjzN35NY+3fMqj5xYREe9WYpDPmjWL9PR0Fi5cSF5e3lmvP/bYY24tzNs476M3IshDg8LoVOtmvtv7LX8nbOWKyCs9en4REfFeJQ6tv/HGGwX3i1sslrO+/I2zR+7pRWGcep1apnXuDk16ExGR00rskTdv3pzmzZvTunVrWrZs6cmavJJRi8I4da59CyEBFZi3azbPtXmx2AmIIiLif0oM8nHjxjFq1CjeeOONYkPjiy++cGth3saoRWGcylnL0a1ud77a/iW/HvuF1tXbGFKHiIh4lxKD3Hl72eOPP+6pWrxa4WvkRunVoC9fbf+SeTu/VpCLiAhwjmvkqamp/Pbbb5hMprO+7Ha7J2v0CqevkRsX5DdGt6dyuSjm75pLbn6uYXWIiIj3KLFH/tZbbwGQk5PDzp07qVevHjabjd27d9OsWTNat27tsSK9gdHXyAGsZis96vfio80fsvrQCjrW6mJYLSIi4h1K7JHPmDGDGTNmULduXZYtW8bcuXP55ptv+OGHH/xuiVZw7PpmMpkMu0bu1OvU4jCzd3xlaB0iIuIdSt005cCBA0RFRRU8rl69OocOHXJrUd7IarUSFhZmaI8coGXVq6lVqTaL9y4iPTfd0FpERMR4JQ6tO4WHh/Pkk0/SokULAP744w+/24/cKTKyMgkJxga5yWSid4O+vLnhdX7Yt5ieDfxrzXsRESmq1B75W2+9RevWrdm7dy979uyhRYsWvP32256ozetERESSlJSIzWYztI5eDfoBMEfD6yIifq/UHnlwcDD9+vXzRC1eLzKyMvn5+aSkJBMeHmFYHZdFXE6Tyk358eAyErMSiAiONKwWERExVqk9cjktMtK5KIyxw+vgmPSWZ8tjwa5vjC5FREQMpCC/AM5b0LwhyHs26I0JE3N2anhdRMSflRrkTzzxBOvWrfNELV7P6GVaC6tRoSbX17yRX47+zIGT+40uR0REDFJqkHfu3Jkvv/ySrl278t5773Hs2DFP1OWVIiIc18WNvgXNqfdljrkLc3dqRzQREX9VapB37dqVyZMn8/XXX1OrVi2GDRvGAw88wE8//eSJ+ryKc711b+iRA9xW93aCLEF8vX2mXy6bKyIi53mNPDMzk6VLlzJ79mxsNhsdOnRg2rRpfncb2umhde/okVcKCuXm2l3ZmbyDzfGbjC5HREQMUGqQP/PMM3Tp0oVNmzYxcuRIZs6cyYABA5gyZQpr1671RI1e4/RkN+/okcPp4XUt2Soi4p9KDfIWLVqwZMkSXnzxRRo2bFjwvMlkYuzYsW4tzts4e+Teco0coGNsZ8KDwpm782vybflGlyMiIh5W6oIwv//+Oxs3biz6TVYrtWvX5s4773RbYd4oJCSE4OBgr+qRB1oC6V6vJ5///T/WHF5F+5ibjC5JREQ8qNQeeUREBIcOHaJOnTrUrVuXAwcOEBISwuHDh3n66ac9UaPXMJlMREREkpiYaHQpRfS5vD8As3fMMrgSERHxtFKDfNeuXXz++ec8+OCDPPDAA0ybNo1Dhw7xwgsvkJyc7JIi/vrrL5599lmeeeYZr99ZLSIi0msmuzldU601sRVrsXD3Au2IJiLiZ0oN8qNHj5KSklLwOCMjg3379pGWlkZ6+rlDY8eOHXTq1Inp06cXPDdp0iTuvPNOevXqxebNmwH4+uuvGT16NA8//DBz58692J/FIyIjI0lPTyMrK8voUgqYTWb6XNaPjLx0Fu9daHQ5IiLiQaUG+Z133kmXLl3o1asXvXv3pmPHjtx+++0sW7aMvn37lvh9GRkZvPTSS7Rt27bgufXr17NlyxZmzpzJ+PHjGT9+PADZ2dkEBgYSFRVFfLz3XH8ujnO9dW+a8AbQ93LHfIWvt880uBIREfGkUie7DRw4kB49erB//35sNhuxsbGEhoaWeuDAwECmTp3K1KlTC5775Zdf6NixIwCXXXYZJ06cIDMzk+DgYLKzszl+/DhVq1a9hB/H/QrfglajRk2DqzmtXlgDWlZtxapDKziefoyqIdWMLklERDzgnEFut9sZOHAgM2bM4Morr7ywA1utWK1FDx8XF1fkFraIiAji4+O58847efHFFwF4/PHHSz12eHh5rFbLBdVTmqioiuf1vpiYGgDk52ee9/d4yr0t7+Hh735nyZEFPHntkxd1DG/7mYyidnBQOzioHRzUDg7e1g7nDHKTyUTjxo2ZNGkSzZs3JzAwsOC1wkPm5ysgIKDIY7vdjslk4oorrigYZj8fSUkZF3zuc4mKqkhcXOp5vTc42PEHuHv3AZo1O7/v8ZSO1boSYA7gk42fcXeDBy/4+y+kHcoytYOD2sFB7eCgdnAwqh3O9Y+HUofWt23bBsCGDRsKnjOZTBcV5FFRUUVmfCcmJhasX+4rvPUaOUBEcCQda3Xh+72L2Bq/hSsrNza6JBERcbNSg3zatGkA5OXlnTVUfqFuvPFG3nrrLQYMGMDWrVuJiYkhODj4ko7pad64TGth/S67i+/3LuLrHTO5svI4o8sRERE3K3XW+ubNm+nevTvdunUD4J133inSOy/Jli1bGDx4MPPmzePzzz9n8ODBREdH07BhQ3r27MmLL77IyJEjL/0n8LDTG6d416IwTp1r30xYUBhzdnylJVtFRPxAqV3sCRMm8NZbbzFmzBgAevTowfDhw5k169yriDVu3LigN1/Y8OHDL7JU73A6yL2zRx5kCeKO+r35bOvHrDq0gptiOxldkoiIuFGpPXKLxUL9+vULHsfGxp41ac2fREREAN55jdzpzoYDAJj1zxcGVyIiIu52XvuRHz58GJPJBMDq1aux2+1uLcqbWa1WwsLCvLZHDtCiSivqhzXgu70LSclONrocERFxo1KDfMSIEQwZMoSNGzfSsmVLJk6cyKhRozxRm9eKjKzsdeutF2Yymbiz4UCy87OZv2ue0eWIiIgblXqNvFGjRixYsIDU1FRMJhMVKlTwRF1eLSIikn379mKz2TCbz2tQw+P6XNafl9ePYeY/X3D3lfcYXY6IiLhJqUG+fft25syZQ0pKSpEh9QkTJri1MG8WGVmZ/Px8UlKSCQ+PMLqcYtWoUJN2MR1YefBHdifvpF5YA6NLEhERNyi1Ozls2DBCQ0Np06YNbdu2LfjyZ968KExh/S93Tnr70uBKRETEXUrtkUdERPDwww97ohaf4bwFLT4+gXr1vLene2ud26gYWImvtn/JiGuew2J27fr0IiJivFJ75DfddBOrV68mJycHm81W8OXPnKu7eXuPvHxAee6o15Mj6YdZc3iV0eWIiIgblNojf++990hLSwMcs6GdG50412D3R857yb35FjSnOxsOYvq2z/hy2zTax9xkdDkiIuJipQb5r7/+6rUzs40SFRUFQHx8nMGVlO7qatfQIOwyvtu7kKSsRMKDvXNynoiIXJwSE/p///uf4w2nQnzz5s0Frz377LNuLsu7RUVVAeDEieMGV1I6k8nEXY0Gk52fzdydXxtdjoiIuFiJQb5y5coij19//fWC/z548KDbCvIFziCPi/P+HjlAv8vvwmKyMGPbdKNLERERFysxyM9chrXwY+dyrf6qcmXH0Hpc3AmDKzk/VcpXoXPtW9gcv4nNcZuMLkdERFyoxCA/V1j781rrAAEBAURERPhMkAMMbHQ3AF9s+9zgSkRExJXOexZb4WD39x45OIbXfSnIO8Z2pkr5qszZ+TVZeVlGlyMiIi5S4qz1P/74g/bt2xc8TkhIoH379tjtdpKSkjxRm1eLiqrC9u3/kJ2dTVBQkNHllMpqttL/8gFM/uMtvtv7Lb0a9DW6JBERcYESg/z777/3ZB0+p0oVx4S3+Pg4ataMNria8zOg0SAm//EW0//+TEEuIlJGlBjkNWvW9GQdPuf0zPUTPhPk9cIacG2N61l7eDV7kndRN6y+0SWJiMgl0kovF6lwkPuSwVf8HwDTNelNRKRMUJBfJF+7l9ypW93biQiOYOY/08nJzzG6HBERuUQK8ovkXKbVF1Z3KyzYGky/ywcQnxnP93sXGV2OiIhcIgX5RapSpSrge0PrcHp4/fO/PzW0DhERuXQK8ovkq9fIARqEX0bbGtex+tAK9qbsMbocERG5BAryi3R6mVbfukbu5OyVf/G3Jr2JiPgyBflFCggIIDw83Cd75AC31b2D8KBwZvwzTZPeRER8mIL8EkRFVfG5yW5OwdZg7mw4iPjMOBbumW90OSIicpEU5JegSpWqJCcnk5Pjmz3afzW+F4D/bZ5qcCUiInKxFOSXwHkLWny8b14nrxtaj5tiO/HrsfVsid9sdDkiInIRFOSXwJdnrjvd0/gBAD7Z8pHBlYiIyMVQkF+CshDknWK7EFMxljk7ZpGclWx0OSIicoEU5JfAGeQnTvhukFvMFv515X1k5GXw+SbdiiYi4msU5JfAuZWpL/fIAQY0GkygOZD3f3sfu91udDkiInIBFOSXoCwMrQNULleZO+r3YnvCdlYe/NHockRE5AIoyC9BWQlygAeu+g8AU/+aYnAlIiJyIRTkl8DXl2ktrFmVFlwXcx3LDvzArqSdRpcjIiLnSUF+CQIDAwkLC/PZ1d3O9HibxwGYulm9chERX6Egv0RVqlQtE0PrAD0a9iC6Qgyz/plBclaS0eWIiMh5UJBfoqioKiQlJZGbm2t0KZfMarZyX5N/k5GXwfRtuhVNRMQXKMgvka8v03qmgY0GU94awsebPyTPlmd0OSIiUgoF+SUqSzPXAcKCw7mz4QAOpx1i0Z4FRpcjIiKlUJBfotOru5WNCW8AD141BBMm3vvjbS0QIyLi5RTkl6hKlapA2bgFzaluWH261u3On3F/sPbwaqPLERGRc1CQXyLnNfKyMrTu9EjzxwF4949JhtYhIiLnpiC/RGXtGrlTi6qtuK7GDaw4uJzN8X8ZXY6IiJRAQX6JymqQAzzS4nEA3lOvXETEaynIL1FUVBVMJhNHjx41uhSX6xDTiSsiG/PNrrnsP7nP6HJERKQYCvJLFBgYSJUqVTl8+JDRpbicyWRiaPPHsNltTPlzstHliIhIMRTkLhAdHcORI4fJz883uhSX61G/N7EVazFj2zSOpx8zuhwRETmDgtwFYmJiyMvL49ixsje8bjVbebTFMLLyszSDXUTECynIXSA6OhaAQ4fK3vA6wJ0NBxJdIYbPtv6P4xllZ+EbEZGyQEHuAtHRMQAcOnTA4ErcI9ASWNArf++Pt40uR0REClGQu0BMjDPIDxpcifvc1WgQNUJq8tnWjzmRUfZutRMR8VUKcheoWdMR5AcPlt0gD7IE8UiLJ8jMy9QMdhERL6Igd4HTPfKyObTuNLDR3VQLqc4nW6aqVy4i4iUU5C5QsWIlQkPDyuS95IUFW4N5vOVTZORlMGnD60aXIyIiKMhdJjo6hoMHD5b5bT8HNfoXtSrV5rOt/2Nfyl6jyxER8XsKcheJiYkhIyOdpKREo0txq0BLIM+0fp5cWy4TfnvF6HJERPyegtxFTt+CVnYnvDn1qN+bxpWvYs6Or9gav8XockRE/JqC3EX8Yea6k9lkZlSbF7Fj55VfxhhdjoiIX1OQu4i/zFx36hDTiWtrXM/S/Uv4+chPRpcjIuK3FOQucnpovWzPXHcymUw839bRGx+1diT5trK3YYyIiC9QkLvI6fXWy/7QulPLqlfT57L+bI7fxMx/vjC6HBERv6Qgd5GoqCiCg4P9KsgBnm8zhvLWEF7+ZQwns1OMLkdExO8oyF3EZDJRs2a031wjd6peoQaPt3yS+Mw43tQiMSIiHqcgd6GaNWNISEggPT3d6FI86j9NhxJbsRZT/5rC7uSdRpcjIuJXFOQu5Jy5XtaXaj1TsDWY0de+TK4tl+fWjijzq9uJiHgTBbkL+dOiMGfqVrc77aI78OOBZXyza47R5YiI+A0FuQv5c5CbTCZebzeJctZyPLd2BElZZXupWhERb6Egd6GYGP+7Ba2w2qF1eOrqZ4jPjGPMuueNLkdExC8oyF3I2SM/eNC/Zq4XNqTpUBpXvooZ/0xj7eHVRpcjIlLmKchdqHr1GpjNZr/tkQNYzVbeaPc2ZpOZp1Y+RkZuhtEliYiUaQpyFwoICKBatep+N2v9TM2rtuTBqx5iT8puxvw8yuhyRETKNAW5i8XExHLkyGGys7ONLsVQz7Z+gYYRjfhky0cs3/+D0eWIiJRZCnIXu/zyRthsNnbs2G50KYYKtgbzXqepBJgDeGzFwyRkJhhdkohImaQgd7HGjZsAsGXLXwZXYrwmla9ixDWjOJFxnKdWPaaFYkRE3EBB7mJNmlwFwObNmwyuxDs83OxR2lS/lkV7FvD5358YXY6ISJmjIHexRo2uxGw2s2XLZqNL8QoWs4X3O00lPCic59Y8zZ8nNhpdkohImaIgd7Hy5ctTv34DtmzZjM1mM7ocrxBdMYYpnT8m15bLfUvuJjFL18tFRFxFQe4GjRtfRVpaKvv37zO6FK9xU2wnhl/9DAdTDzBk6f3k2/KNLklEpExQkLtB48aO6+Sa8FbUsFZP0zG2MysOLue1X182uhwRkTJBQe4GzglvCvKizCYz73eaSq1KtZm0cSIztk0zuiQREZ+nIHcDZ49882YF+ZnCgyOYedscwoPCeWrVY6w8+KPRJYmI+DQFuRtERkZSo0ZNBXkJ6oU14LOuMzFj5r4ld/N3wlajSxIR8VkKcjdp0uQqjh8/xokTJ4wuxSu1qd6WyR0/IDXnJHct7M3elD1GlyQi4pMU5G5y5ZVa4a00PRv0Yex1r3A0/Qi953fnwMn9RpckIuJzFORu0qRJUwAtDFOK/zQdyqg2ozmUdpBeC7pzONW/d44TEblQCnI3Ob3mupZqLc2jLYbx9NXPcuDkPnrO76aeuYjIBVCQu0lsbC1CQ8PUIz9PT109kmGtnmbfyb10m9uZrfFbjC5JRMQnKMjdxGQy0bhxE3bv3kVaWprR5fiEkdeMYux1r3A84xi3f3MLPx1eY3RJIiJeT0HuRo0bN8Fut7N1q3qX5+s/TYfyQeePycrLpP+3PZm9Y5bRJYmIeDUFuRtdc01bAJYtW2JwJb6lV4O+zOg2m0BLEA8te4BRa0eQm59rdFkiIl5JQe5GnTp1ISSkAvPmzcFutxtdjk9pF9OBH/qs5PLwhvz3ryn0XtCd4xnHjS5LRMTrKMjdqFy5ctx6azcOHNjHhg2/GV2Oz6kf3oDFvZfTvV4P1h9dR4dZ17J47yKjyxIR8SoKcjfr1asPAPPmzTa4Et9UIbAiH3X5jLHXvUJqzkn+tfguHv1xCCezU4wuTUTEKyjI3axdu5uIiIhg/vx55OdrD+6LYTKZ+E/ToSztu5omlZsy858vaDerLQt3L9AlCxHxewpyNwsICOC223pw4sRx1q1ba3Q5Pq1hRCMW917OsFZPcyLjOPcuGUT/hT3ZlbTT6NJERAyjIPeAnj17Axped4VASyAjrxnFqv7raR9zEysP/ki7WW0YufpJjqUfNbo8ERGPU5B7QJs211KtWnUWLpxPTk6O0eWUCfXDGzDrtnl8cssX1KhQk/9tmco105vy/E/PcCJDO86JiP9QkHuAxWLhjjt6kZyczIoVy40up8wwmUx0q9uddXdt4M32k4ksV5kPN71Hi8+v4NEfh2iZVxHxCwpyD+nduy8A77zzpia9uViAJYBBV/yL9QP/4LUb3yS6Ygwz//mCDl9dS89vuvHV9i/JyM0wukwREbdQkHtI06bNueOOXvz22y9MnvyW0eWUSUGWIO5pfD/rBmxgetdZ3FCzHT8dWcPQ5f+m8acNGLbiEVYe/FGrxIlImWKy++D9O3FxqS49XlRURZcfszhJSYm0b38tcXEn+O67ZTRr1sLt57wQnmoHT9qTspuv/pnBrO1fcjjNsdd5WFAYN9fuSpfat3JDzRsJCw4v8j1lsR0uhtrBQe3goHZwMKodoqIqlviaghzP/sGsXr2SPn1up169+ixbtoaQkBCPnPd8lOW/qPm2fNYfXceiPQtYtOdbjqYfAcBsMtO8SktujG5H6+ptaVX1GupFR5fZdrgQZfnzcCHUDg5qBwcFuYv4cpADvPDCs3zwwbv06dOft956l6CgII+d+1z85S+qzW7jzxMb+fHAMlYe/JENx38j3+6Yt2DCRJOqTWgc3pSrqjTjqspNaRh5BRUCKhhctef5y+ehNGoHB7WDg4LcRXw9yLOysujatRNbtvxFo0ZX8M47U2jatLnHzl8Sf/2LejI7hd+O/cIvR9fzy7Gf+fPERjLzMou8J6ZiLJeHN6R++GXUCa1LndC61K5UhxoVahJoCTSocvfy18/DmdQODmoHBwW5i/h6kAOkpp5k7NgX+eyzj7FYLAwZ8giDBt1N3br1PVpHYfqL6hAeWY51OzbwV9yfbI7fxLbEbWxP3MaJYnZfM2GiWkh1oivGUC2kOtXKV6NqSHWqlK9C5XKVqVwuiojgSCKCIwgJqIDJZDLgJ7o4+jw4qB0c1A4OCnIXKQtB7rR69UqGDXuEAwf2A1C/fgM6d76FJk2uolat2tSqVYeoqCiPBID+ojqU1A5JWYnsTt7F3pQ97EnZzf6T+zicdohDqQc5kna4YHi+JAHmAMKCwqkUVIlKgZWoGBhKxcCKVAioQIXACoRYK1A+oDzlA8pTzlqectZylLOWI9gSTJA1mCBLMMGWIAItQQRaAhz/bw7Eag4g0BKA1RxAgDkAi8niks+LPg8OagcHtYODgtxFylKQA6SlpTF//lyWLFnM6tUryMgoes+z2WymUqVKVKoUSkhIBYKCAgkMDCIoKAiz2YzFYsFisWA2m0/9AjdhMp3+Aor8Yi/pl3xQkJWcHM/f437TTZ0YMGCwx89bkov5POTb8onPiud4+lGOpR8lPjOe+Mw44jLjSMiMJzkriaTsJJKzkziZfZKTOSlk52e76ScAq9mK1WTFYrZiNVuxmMyYTRYszi+zBbPJjMVU9P9NJjPmU1+BAVby82yYTWbA5Hjd+f/Oz9Wp/5lNZjA5HzmYTKcfnfn+s54r/PnkjM9nCa+d9b4i33KO187xfcV9f1CQlezsvFK/xxXOpzajeLIdvNn5tkPjyk14vOVTLjuvgrwURgd5YZmZmfz22y/s2rWT/fv3sW/fXuLj40hNPcnJkydJS0sjNzeH7OxsbDab0eW6xA03tGfOnAVGl1HAU5+HrLws0nLTSMtJJS03jfTcdDJy08nMyyQjL52svCyy8jLJzM8iOy+L7PxssvKzyMnPJjc/lxxbjuO/bXnk2XLJyc8hz55Pni2X3PxcbPb8U4/zyLflkW/PP/Vlw2bLL3hss9uw222n/tuOHTs2uw2bPR97ocenny8bnzsRd4qtVJvfBm5y2WiqgrwU3hTkFyIvL4/8/Hzy8/Ox2fKx2WyOX7yFvgAK/wmf6487MjKE+Pg0d5d9loiICCwWi8fPWxJf/Ty42vm0gzPo7fbTYQ8UPOf878LvLfwchZ5zvqfI8Ut4rfDzZ9V0rtfO47fdmd9fObIC8Qnu/3vh7b+KK1euYMjvB29zvu0QGhTq0omw5wpyq8vOIh5ntVqxWl33RxgVVRGTqZzLjidlX8HwuQkseM8/xlwpKqQiZAQbXYbhokIqYsrQ7wdvbAct0SoiIuLDFOQiIiI+TEEuIiLiwxTkIiIiPkxBLiIi4sMU5CIiIj5MQS4iIuLDFOQiIiI+TEEuIiLiwxTkIiIiPkxBLiIi4sN8ctMUERERcVCPXERExIcpyEVERHyYglxERMSHKchFRER8mIJcRETEhynIRUREfJjfB/mkSZO488476dWrF5s3bza6HI9688036d+/P7169WLx4sUkJCRw33330a9fPx599FFycnKMLtFjsrKy6NSpE3PnzvXbdvj222/p1asXPXv2ZOXKlX7ZDunp6Tz88MMMHjyYfv36sWrVKvbu3cugQYPo3bs3o0ePpqzfsbtjxw46derE9OnTAUr8HCxdupT+/fvTo0cPZs+ebWTJbnFmOxw/fpx7772XQYMGcffdd3P8+HHAS9rB7sd+/vln+3333We32+327du32wcMGGBwRZ7z66+/2u+//3673W63JyUl2W+44Qb7008/bV+0aJHdbrfbx48fb//666+NLNGj3nzzTXuvXr3sc+bM8ct2SEtLs/fs2dOelZVlP3bsmP25557zy3aYNm2a/fXXX7fb7Xb7sWPH7F26dLEPHDjQ/ueff9rtdrv9kUcesa9bt87IEt0qPT3dPmjQIPuoUaPs06ZNs9vt9mI/B6mpqfZOnTrZT548ac/IyLDffPPN9rS0NCNLd6ni2mHkyJEF7TB9+nT7+PHjvaYd/LpH/ssvv9CxY0cALrvsMk6cOEFmZqbBVXlG8+bNmTRpEgAVK1YkNzeX9evXc9NNNwHQsWNH1q5da2CFnrN79252795N+/btAfj111/9rh3Wrl1Lu3btCAoKomrVqowbN84v2yE8PJyEhAQAkpOTCQ8P58CBAzRt2hSAm266qUy3Q2BgIFOnTqVKlSoFzxX3Odi8eTONGzemYsWKlCtXjhYtWvD7778bVbbLFdcOo0aNokuXLoDjc5KWluY17eDXQR4XF0dERETB44iICOLj4w2syHOsVishISEAzJ49m3bt2pGZmUlwcDDgX20xYcIERo4cWfA4PT3d79rh6NGjJCYmcv/99zNgwAB+/vlnv2yHrl27cvToUW6++Wb+9a9/MXz4cEJDQwtej4yMLNPtYLVaC/7MnYr7HJz5u7OstUtx7RASEoLVaiU/P58ZM2bQrVs3r2kHq8fP6EUCAgKKPLbb7ZhMJoOqMcayZcv46quv+OSTT1izZk3B8/7SFt988w2tWrUiOjq64LnCnwt/aYecnBwSEhL48MMPOXjwIP/3f/+HxWIpeN1f2mH+/PnUrFmTTz/9lH/++YehQ4dSrly5gtf9pR0KK+7vg7/+7szPz+fpp5/mmmuuoU2bNixevLjI60a1g18HeVRUVMEwGkBiYiKVK1c2sCLPWrNmDe+//z4ff/wxlSpVIiQkhMzMTMqVK0d8fHyRYaWyauXKlRw6dIilS5dy7NgxAgMDCQoK8rt2iIqKolmzZlgsFmrXrk2FChUwm81+1w5//PEHN954IwANGzYkKyuLrKysgtf9pR0KK+73wpm/O+Pj42nTpo2BVXrGM888Q3R0NI8++ihwdoYY1Q5+PbR+4403snz5cgC2bt1KTEzMWcMpZVVqairjx4/nv//9L+Hh4QDccMMNBe2xdOlS2rVrZ2SJHjFp0iRmz57NV199Rd++fXnooYfo0KGD37XDtddey/r167Hb7SQkJJCenu6X7RAbG8uWLVsAxyzlkJAQGjduzB9//AH4TzsUVtzvhauuuort27eTmppKeno6mzZtolWrVgZX6l4LFizAbDbzxBNPFDznLe3g97ufvf7666xbtw6LxcLLL7/M5ZdfbnRJHjFr1iwmT55MnTp1Cp4bP348I0eOJCMjgzp16jB+/HisVv8ZtJk8eTI1a9bk+uuv58knn/S7dpg1axYLFy4suAWrSZMmftcO6enpjBw5kqSkJHJzc3n88ceJiorimWeeIT8/n2uuuabIfIqyZsuWLbz22mscPnwYq9VK1apVmThxIk899dRZn4PFixczZcoUzGYz999/P7fddpvR5btMce2QkJBAUFAQFSpUAKBevXqMHj3aK9rB74NcRETEl/n10LqIiIivU5CLiIj4MAW5iIiID1OQi4iI+DAFuYiIiA8r2/eSiMh5WbVqFR9++CFWq5WMjAxiYmIYM2YMu3btIioqipiYGKNLFJES6PYzET+Xk5PD9ddfz8KFCwtWLZswYQJRUVHs3r2brl27cu211xpcpYiURD1yET+Xk5NDRkYG6enpBc89/fTTLF26lPfee4+//vqrYGnKsWPHkp2dTU5ODkOGDKFdu3YMHz6ckJAQdu/eTWJiIn369OGee+4x8CcS8S8KchE/V6FCBR577DF69epF06ZNad26NTfffDOdO3fm888/Z8iQIbRt25YHH3yQBx98kKuvvpqEhAR69+7NDz/8gMViIT4+nmnTppGSkkKXLl3o0aNHwdK/IuJeCnIR4YEHHqBPnz6sW7eO9evX069fP4YNG1bkPRs3bmTSpEmYzY45skFBQcTFxQHQunVrAEJDQ4mNjeXAgQMKchEPUZCLCJmZmYSHh9OtWze6devGLbfcwvjx4wkLCyt4j8lkYvLkyUX2X3YqPNXGZrP5xZaWIt5Ct5+J+Lk1a9bQr18/0tLSCp47cuQItWvXxmQykZubC0DLli35/vvvAUhOTmbcuHEF7//5558BSElJ4dChQ9SuXdtzP4CIn9OsdRFh+vTpzJs3j/LlywOOfZafffZZ5syZw6xZs3jmmWdo1KgRL7zwAtnZ2WRmZvLggw9yyy23MHLkSEwmE8ePH+fEiRP079+fwYMHG/wTifgPBbmIXJKRI0fSsmVL+vbta3QpIn5JQ+siIiI+TD1yERERH6YeuYiIiA9TkIuIiPgwBbmIiIgPU5CLiIj4MAW5iIiID1OQi4iI+LD/B0NBc5HxqS0+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x396 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.style.use(\"seaborn\")\n",
    "plt.plot(np.array(gd_cost) - exact_value, \"g\", label=\"Gradient descent\")\n",
    "plt.plot(np.array(qngd_cost) - exact_value, \"k\", label=\"Quantum natural gradient descent\")\n",
    "plt.yscale(\"log\")\n",
    "plt.ylabel(\"Energy difference\")\n",
    "plt.xlabel(\"Step\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that by employing quantum natural gradients, it takes fewer steps\n",
    "to reach a ground state estimate and the optimized energy achieved by\n",
    "the optimizer is lower than that obtained using vanilla gradient descent.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Robustness in parameter initialization\n",
    "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n",
    "\n",
    "While results above show a more rapid convergence for quantum natural gradients,\n",
    "what if we were just lucky, i.e., we started at a \"good\" point in parameter space?\n",
    "How do we know this will be the case with high probability regardless of the\n",
    "parameter initialization?\n",
    "\n",
    "Using the same system Hamiltonian, ansatz, and device, we tested the robustness\n",
    "of the ``QNGOptimizer`` by running 10 independent trials with random parameter initializations.\n",
    "For this numerical test, our optimizer does not terminate based on energy improvement; we fix the number of\n",
    "iterations to 200.\n",
    "We show the result of this test below (after pre-computing), where we plot the mean and standard\n",
    "deviation of the energies over optimization steps for quantum natural gradient and standard gradient descent.\n",
    "\n",
    ".. figure:: ../demonstrations/vqe_qng/k_runs_.png\n",
    "    :align: center\n",
    "    :width: 60%\n",
    "    :target: javascript:void(0)\n",
    "\n",
    "We observe that quantum natural gradient on average converges faster for this system.\n",
    "\n",
    "<div class=\"alert alert-info\"><h4>Note</h4><p>While using QNG may help accelerate the VQE algorithm in terms of optimization steps,\n",
    "    each QNG step is more costly than its vanilla gradient descent counterpart due to\n",
    "    a greater number of calls to the quantum computer that are needed to compute the Fubini-Study metric tensor.</p></div>\n",
    "\n",
    "While further benchmark studies are needed to better understand the advantages\n",
    "of quantum natural gradient, preliminary studies such as this tutorial show the potentials\n",
    "of the method. 🎉\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "References\n",
    "--------------\n",
    "\n",
    ".. [#stokes2019]\n",
    "\n",
    "    Stokes, James, *et al.*, \"Quantum Natural Gradient\".\n",
    "    `arXiv preprint arXiv:1909.02108 (2019).\n",
    "    <https://arxiv.org/abs/1909.02108>`__\n",
    "\n",
    ".. [#yamamoto2019]\n",
    "\n",
    "    Yamamoto, Naoki, \"On the natural gradient for variational quantum eigensolver\".\n",
    "    `arXiv preprint arXiv:1909.05074 (2019).\n",
    "    <https://arxiv.org/abs/1909.05074>`__\n",
    "\n",
    ".. [#peruzzo2014]\n",
    "\n",
    "    Alberto Peruzzo, Jarrod McClean *et al.*, \"A variational eigenvalue solver on a photonic\n",
    "    quantum processor\". `Nature Communications 5, 4213 (2014).\n",
    "    <https://www.nature.com/articles/ncomms5213?origin=ppub>`__\n",
    "\n"
   ]
  }
 ],
 "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}