ising-ham-trotter-error-study.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import trace, dot, pi, exp\n",
    "\n",
    "from scipy.stats import sem\n",
    "from scipy.linalg import norm, expm, kron\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "from functools import reduce\n",
    "from itertools import product\n",
    "\n",
    "import random\n",
    "\n",
    "from qiskit import Aer, QuantumCircuit, QuantumRegister, ClassicalRegister, execute\n",
    "from qiskit.quantum_info import Operator\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "s0 = np.array([1, 0]) # |0>\n",
    "s1 = np.array([0, 1]) # |1>\n",
    "P_0 = np.array([[1,0],[0,0]]) # |0><0|\n",
    "P_1 = np.array([[0,0],[0,1]]) # |1><1|\n",
    "X_np = np.array([[0.,1.],\n",
    "                 [1.,0.]]) # X Pauli matrix\n",
    "Y_np = np.array([[0.,-1.j],\n",
    "                 [1.j, 0.]]) # Y Pauli matrix\n",
    "Z_np = np.array([[1., 0.],\n",
    "                 [0.,-1.]]) # Z Pauli matrix\n",
    "I_np = np.array([[1.,0.],\n",
    "                 [0.,1.]]) # 2x2 identity matrix\n",
    "H_np = (X_np + Z_np)/np.sqrt(2) # Hadamard gate\n",
    "S_np = np.array([[1.,0. ],\n",
    "                 [0.,1.j]]) # phase gate\n",
    "T_np = np.array([[1.,0. ],\n",
    "                 [0.,exp(1j*pi/4)]]) # phase gate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ising_model_dict(qubits_number, J, hx):\n",
    "    '''\n",
    "    returns a Hamiltonian represented as a dictionary\n",
    "    example: {'ZZ': -1, 'XI': 0.5, 'IX': 0.5}\n",
    "    '''\n",
    "\n",
    "    ham = {}\n",
    "\n",
    "    line = 'Z' + 'Z' + 'I' * (qubits_number - 2)\n",
    "    for i in range(qubits_number):\n",
    "        term = line[-i:] + line[:-i]\n",
    "        ham[term] = J\n",
    "        \n",
    "    line = 'X' + 'I' * (qubits_number - 1)\n",
    "    if hx != 0:\n",
    "        for i in range(qubits_number):\n",
    "            term = line[-i:] + line[:-i]\n",
    "            ham[term] = hx\n",
    "\n",
    "    return ham"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def string_to_matrix(operator_string):\n",
    "    \n",
    "    qubits_number = len(operator_string)\n",
    "\n",
    "    if operator_string[0] == 'I':\n",
    "        operator = I_np\n",
    "    elif operator_string[0] == 'X':\n",
    "        operator = X_np\n",
    "    elif operator_string[0] == 'Y':\n",
    "        operator = Y_np\n",
    "    elif operator_string[0] == 'Z':\n",
    "        operator = Z_np\n",
    "\n",
    "    for k in range(1, qubits_number):\n",
    "        if operator_string[k] == 'I':\n",
    "            operator = kron(operator, I_np)\n",
    "        elif operator_string[k] == 'X':\n",
    "            operator = kron(operator, X_np)\n",
    "        elif operator_string[k] == 'Y':\n",
    "            operator = kron(operator, Y_np)\n",
    "        elif operator_string[k] == 'Z':\n",
    "            operator = kron(operator, Z_np)\n",
    "            \n",
    "    return operator\n",
    "\n",
    "\n",
    "def hamiltonian_matrix_to_dict(n_qubits, hermitian_matrix):\n",
    "    \n",
    "    d = 2**n_qubits\n",
    "        \n",
    "    paulies = ['I', 'X', 'Y', 'Z']\n",
    "    bases = list(product(*[paulies for i in range(n_qubits)]))\n",
    "    \n",
    "    hermitian_dict = {}\n",
    "    for basis in bases:\n",
    "        observable = string_to_matrix(basis)\n",
    "        result = trace(hermitian_matrix.dot(observable)).real\n",
    "        \n",
    "        if result != 0:\n",
    "            label = ''.join(basis) # pythonic hacks for concatinating chars in a list\n",
    "            hermitian_dict[label] = result / d\n",
    "    \n",
    "    return hermitian_dict\n",
    "\n",
    "\n",
    "def hamiltonian_dict_to_matrix(hamiltonian_dictionary):\n",
    "    '''\n",
    "    transforms a Hamiltonian dictionary into the NumPy matrix\n",
    "    '''\n",
    "\n",
    "    hamiltonian_m = 0\n",
    "\n",
    "    for operator, multiplier in hamiltonian_dictionary.items():\n",
    "\n",
    "        hamiltonian_term = [] \n",
    "        for op in operator:\n",
    "            if op == 'I':\n",
    "                hamiltonian_term.append(I_np)\n",
    "            elif op == 'X':\n",
    "                hamiltonian_term.append(X_np)\n",
    "            elif op == 'Y':\n",
    "                hamiltonian_term.append(Y_np)\n",
    "            elif op == 'Z':\n",
    "                hamiltonian_term.append(Z_np)\n",
    "            \n",
    "        hamiltonian_m = hamiltonian_m + multiplier*reduce(kron, hamiltonian_term)\n",
    "\n",
    "    return hamiltonian_m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hea(qreg, creg, n_layers, pars):\n",
    "  # generates a circuit for the hardware-efficient ansatz\n",
    "    \n",
    "    circuit = QuantumCircuit(qreg, creg)\n",
    "    shift = 0\n",
    "    for l in range(n_layers):\n",
    "        for q in range(n_qubits):\n",
    "            circuit.rx(pars[3*q+shift], q)\n",
    "            circuit.rz(pars[3*q+shift+1], q)\n",
    "            circuit.rx(pars[3*q+shift+2], q)\n",
    "        shift += 3*n_qubits\n",
    "        if n_qubits > 1:\n",
    "            for q in range(n_qubits - 1):\n",
    "                circuit.cry(pars[q+shift], q, q+1)\n",
    "            circuit.cry(pars[q+shift+1], n_qubits-1, 0)\n",
    "        shift += n_qubits\n",
    "        circuit.barrier()\n",
    "        \n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Constructors of the exact and approximated unitary evolutions of the Ising Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ising_unitary_trotter(qreg, creg, t, k, J, h):\n",
    "  # returns the k-step Trotter decomposition of the Ising Hamiltonian with the coupling strength J, transverse field h and simulation time t\n",
    "    \n",
    "    n_qubits = qreg.size\n",
    "    \n",
    "    T = 2*t / k # multiplied by two because of the definition of rotation gates in qiskit: R_A(x) = e^(-i*A*x/2)\n",
    "    \n",
    "    circuit = QuantumCircuit(qreg, creg)\n",
    "    for i in range(k):\n",
    "        for j in range(n_qubits - 1):\n",
    "            circuit.rzz(J * T, j, j + 1)\n",
    "        if n_qubits > 2:\n",
    "            circuit.rzz(J * T, 0, n_qubits - 1)\n",
    "        for j in range(n_qubits):\n",
    "            circuit.rx(h * T, j)\n",
    "\n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ising_unitary_exact(qreg, creg, t, J, h):\n",
    "  # returns the exact unitary evolution of the Ising Hamiltonian with the coupling strength J, transverse field h and simulation time t\n",
    "    \n",
    "    n_qubits = qreg.size\n",
    "    \n",
    "    circuit = QuantumCircuit(qreg, creg)\n",
    "    matrix = expm(-1j * t * hamiltonian_dict_to_matrix(ising_model_dict(n_qubits, J, h)))\n",
    "    circuit.append(Operator(matrix), [i for i in range(n_qubits)])\n",
    "    \n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Test: Hilbert-Schmitd distance between the approximate and the exact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hs_distance(U, V):\n",
    "  # Hilbert-Schmidt distance\n",
    "    d = len(U)\n",
    "    return 1 - abs(trace(dot(U, V.conj()))) / d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmx0lEQVR4nO3deZxU1Zn/8c/TC0sD3exrb6Cg0SCoLUJVjDEYt2iMWzRBRxknxKgxTswkGjKTGAeXX0yGZFwZI2okLjEmgmJciCbKooCigEZFFLrBBVyapVl6eX5/VHVTvVVXI9W3lu+b133V3ercp+vV9FPnnHvONXdHRESyV07QAYiISLCUCEREspwSgYhIllMiEBHJckoEIiJZLi/oADpr4MCBXl5eHnQYIiJpZfny5ZvdfVBbx9IuEZSXl7Ns2bKgwxARSStmtq69Y2oaEhHJckoEIiJZTolARCTLKRGIiGQ5JQIRkSyXFYlgzso5lM8sJ+fqHMpnljNn5ZygQxIRSRlpd/toZ81ZOYdp86ZRU1sDwLrqdUybNw2AKWOnBBmaiEhKyPgawfQF05uSQKOa2hqmL5geUEQiIqkl4xPB+ur1ndovIpJtMj4RlBaVdmq/iEi2yfhEMGPyDAryC5rtK8gvYMbkGQFFJCKSWjI+EUwZO4VZp8yirKgMgBzL4faTb1dHsYhIVFITgZmdYGZvmNkaM7uyjeNFZjbPzF4xs9VmNjUZcUwZO4V3L3+X333tdzR4AxXDK5JxGRGRtJS0RGBmucDNwInAQcA3zeygFqddArzm7uOALwG/MrNuyYopVBICYOH6hcm6hIhI2klmjWACsMbd17r7buB+4NQW5zjQx8wM6A18DNQlK6ADBhxA/579WVS5KFmXEBFJO8lMBCOAypjtqui+WDcBnwM2AiuB77t7Q7ICMjNCJSEWVqpGICLSKJmJwNrY5y22jwdWAMOB8cBNZlbYqiCzaWa2zMyWbdq06TMFFS4J88ZHb7C5ZvNnKkdEJFMkMxFUASUx28VEvvnHmgo87BFrgHeAA1sW5O6z3L3C3SsGDWrzSWsJa+wnWFK15DOVIyKSKZKZCJYCo81sZLQD+Bxgbotz1gOTAcxsCHAAsDaJMXHE8CPIy8lTh7GISFTSJp1z9zozuxR4AsgF7nT31WZ2UfT4bcA1wF1mtpJIU9KP3T2pbTY983ty2LDDWFSlDmMREUjy7KPuPh+Y32LfbTHrG4HjkhlDW8IlYW5ddiu763fTLTdpd6uKiKSFjB9Z3JZQSYiddTtZ8f6KoEMREQlc1iYC0MAyERHI0kQwvM9wyvuWq59ARIQsTQQQqRUsXL8Q95ZDG0REskvWJoJwSZj3tr3Huup1QYciIhKorE0E6icQEYnI2kQwdvBYenfrrQnoRCTrZW0iyM3JZWLxRHUYi0jWy9pEAJF+glc/eJWtu7YGHYqISGCyOhGESkI0eAMvbHgh6FBERAKT1YlgYvFEDFOHsYhktaxOBIXdCxk7ZKz6CUQkq2V1IoBIP8HiysXUN9QHHYqISCCyPhGESkJs3b2V1ZtWBx2KiEgglAg0sExEslzWJ4KRfUcytPdQ9ROISNbK+kRgZk0T0ImIZKOsTwQQ6TB+59N3eG/re0GHIiLS5ZQI2NNPoHmHRCQbKREAhw07jO653ZUIRCQrKREA3XK7ccSII9RhLCJZSYkgKlwSZvnG5eyo3RF0KCIiXUqJICpUEqK2oZbl7y0POhQRkS6lRBA1qXgSoIFlIpJ9lAiiBvUaxJgBY9RPICJZR4kgRqgkxKLKRbh70KGIiHSZhBKBmfU0swOSHUzQwiVhNtds5q2P3wo6FBGRLtNhIjCzU4AVwF+j2+PNbG6S4wqEJqATkWyUSI3g58AE4FMAd18BlCcroCAdOPBA+vXop4FlIpJVEkkEde5enfRIUkCO5TCpZBILK1UjEJHskUgiWGVm3wJyzWy0mf0vkLFfmcMlYV7f/Dof7/g46FBERLpEIonge8DBwC7gD0A1cHkSYwpUYz/BkqolAUciItI1OkwE7l7j7tPd/Yjo8lN339kVwQVhwogJ5FquOoxFJGskctfQU2bWN2a7n5k9kdSoAlSQX8Chww7VwDIRyRqJNA0NdPdPGzfc/RNgcNIiSgGh4hAvVL1AbX1t0KGIiCRdIomgwcxKGzfMrAzI6KG34dIwO+p28MoHrwQdiohI0uUlcM504Hkz+3t0+4vAtOSFFLzYgWUVwysCjkZEJLkS6Sz+K3AY8ADwIHC4u2dsHwFAcWExpUWl6icQkayQ6KRz3YGPidw6epCZfTF5IaWGUEmIhesXagI6Ecl4idw1dAOwkEgT0X9Elx8mUriZnWBmb5jZGjO7sp1zvmRmK8xsdUzzU+DCJWE2bN1A5ZbKoEMREUmqRPoIvg4c4O67OlOwmeUCNwNfAaqApWY2191fizmnL3ALcIK7rzezlLkbKbafoHRsaQdni4ikr0SahtYC+XtR9gRgjbuvdffdwP3AqS3O+RbwsLuvB3D3D/fiOklxyJBD6JXfSxPQiUjGS6RGUAOsMLMFRKaZAMDdL+vgfSOA2HaVKuDIFueMAfLN7FmgD/Abd7+nZUFmNo3onUqlpV3z7TwvJ48ji4/UBHQikvESSQRzo0tnWRv7Wva85gGHA5OBnsBiM1vi7m82e5P7LGAWQEVFRZf13oZLwsx4bgbbdm+jd7feXXVZEZEu1WEicPe797LsKqAkZrsY2NjGOZvdfTuw3cz+AYwD3iQFhEpCNHgDL254kS+P/HLQ4YiIJEUidw2NNrOHzOw1M1vbuCRQ9lJgtJmNNLNuwDm0rlk8AhxlZnlmVkCk6ej1zv4QyTKxeCKGaQI6EcloiTQNzQZ+BvwPcAwwlbabfZpx9zozuxR4AsgF7nT31WZ2UfT4be7+upn9FXgVaADucPdVe/ej7Ht9e/Tl4MEHa2CZiGS0RBJBT3dfYGbm7uuAn5vZc0SSQ1zuPh+Y32LfbS22fwn8shMxd6lQcYgHVj9AgzeQY4mOvxMRSR+J/GXbaWY5wFtmdqmZnUaGzz4aK1wapnpXNa9teq3jk0VE0lAiieByoAC4jMgdPucC/5LEmFJK7MAyEZFMlEgiKHf3be5e5e5T3f0MIGuG2u7Xbz8G9xqsfgIRyViJJIKrEtyXkcysaQI6EZFM1G5nsZmdCJwEjDCz38YcKgTqkh1YKgmXhPnLP//CB9s+YEjvIUGHIyKyT8WrEWwElgE7geUxy1zg+OSHljoa+wk075CIZKJ2awTu/grwipn9wd1rIfLgeqAk+tzirHH4sMPpltuNRZWLOO1zpwUdjojIPpVIH8FTZlZoZv2BV4DZZvbrJMeVUrrndadieIU6jEUkIyWSCIrcfQtwOjDb3Q8Hjk1uWKknVBxi2cZl7KzbGXQoIiL7VCKJIM/MhgHfAB5NcjwpK1waZnf9bl5676WgQxER2acSSQS/IDJf0Bp3X2pmo4C3khtW6plUPAnQwDIRyTyJTEP9R+CPMdtrgTOSGVQqGtJ7CPv331/9BCKSceKNI/iRu/8/M/tfWj9QJpEnlGWcUEmIx996HHfHrMMJWEVE0kK8pqHG5wIso/k4gsYl64RLwmyq2cTbn7wddCgiIvtMvHEE86Kve/uEsowTOwHd/v33DzgaEZF9I17T0DzaaBJq5O5fS0pEKeygQQdR1L2IRZWLOH/8+UGHIyKyT8TrLL4x+no6MBS4N7r9TeDdJMaUsnIsh0klk1hYqTuHRCRztNtH4O5/d/e/A4e6+9nuPi+6fAv4QteFmFpCxSFWb1rNpzs/DToUEZF9IpFxBIOiYwcAMLORwKDkhZTawqVhABZXLg44EhGRfSORRPDvwLNm9qyZPQs8A3w/qVGlsAkjJpBruZqJVEQyRiIDyv5qZqOBA6O7/unuu5IbVurq3a0344aO08AyEckYidQIcPdd7v5KdMnaJNAoVBzihaoXqGvIqufziEiGSigRSHPh0jDba7fz6gevBh2KiMhnpkSwF2IHlomIpLsOE4GZLUhkXzYpLSqluLBY/QQikhHijSzuARQAA6OPqGycZa0QGN4FsaW0UElINQIRyQjxagTfITK53IHAS+yZbO4R4Obkh5bawiVhKrdUUlldGXQoIiKfSbyRxb9x95HAD919ZMwyzt1v6sIYU1JjP4HGE4hIuovXNHR6dHVDzHoTd384aVGlgXFDxlGQX8CiykWc/fmzgw5HRGSvxRtQdkr0dTAQAv4W3T4GeBbI6kSQn5vPhBETNAGdiKS9eE1DU919KpGpqA9y9zPc/Qzg4C6LLsWFikOseH8F23dvDzoUEZG9lsg4gnJ3fy9m+wNgTJLiSSvh0jD1Xs/SjUuDDkVEZK8lkgieNbMnzOwCMzsfeIzIxHNZb2LxREADy0QkvSUy6dyl0c7io6K7Zrn7n5MbVnro37M/Bw06SAPLRCStdZgIoOkOoazuHG5PqDjEQ68/RIM3kGOasUNE0k+7f7nM7Pno61Yz2xKzbDWzLV0XYmoLl4b5dOen/HPzP4MORURkr8S7a+gL0dc+7l4Ys/Rx98KuCzG1aQI6EUl3CbVlmFk/MzvEzA5rXJIdWLoY3X80AwsGqp9ARNJWh30EZnYNcAGwFmiI7nbgy8kLK32YmSagE5G0lkiN4BvAfu5+tLsfE10SSgJmdoKZvWFma8zsyjjnHWFm9WZ2ZqKBp5JQcYi3Pn6LTds3BR2KiEinJZIIVgF9O1uwmeUSmaX0ROAg4JtmdlA7590APNHZa6SKcGkY0AR0IpKeEkkE1wEvRweVzW1cEnjfBGCNu691993A/cCpbZz3PeBPwIcJR51iDh92OPk5+UoEIpKWEhlHcDeRb+wr2dNHkIgRQOxk/VXAkbEnmNkI4DQi/Q1HtFeQmU0DpgGUlpZ2IoSu0TO/J4cPP1wT0IlIWkokEWx299/uRdnWxj5vsT0T+LG715u1dXr0Te6zgFkAFRUVLctICaHiEDcvvZlddbvontc96HBERBKWSNPQcjO7zswmdfL20SqgJGa7GNjY4pwK4H4zexc4E7jFzL6eQNkpJ1waZlf9Ll5+/+WgQxER6ZREagSHRl8nxuxL5PbRpcBoMxsJbADOAb4Ve0L0CWgAmNldwKPu/pcEYko5sQPLGiejExFJB4lMOnfM3hTs7nVmdimRu4FygTvdfbWZXRQ9ftvelJuqhvYeyqh+o1hUtYgruCLocEREEpbIgLLuwBlAeez57v6Ljt7r7vOB+S32tZkA3P2CjspLdaGSEE+9/RTuTrw+DxGRVJJIH8EjRG77rAO2xyzSQrgkzAfbP+CdT98JOhQRkYQl0kdQ7O4nJD2SDBDbTzCq36iAoxERSUwiNYJFZjY26ZFkgIMHHUxh90INLBORtNJujcDMVhK5OygPmGpma4FdRMYHuLsf0jUhpo/cnFwmFk/UwDIRSSvxmoZO7rIoMkioOMTVf7+a6p3VFPUoCjocEZEOxXswzTp3XwcMAz6O2f4YGNpVAaabcGkYx1lStSToUEREEpJIH8GtwLaY7e3RfdKGI0ccSY7lqJ9ARNJGIonA3L1pfh93byDBh95noz7d+3DIkEP0xDIRSRuJJIK1ZnaZmeVHl+8TeVqZtCNUHGJJ1RLqGuqCDkVEpEOJJIKLgBCR+YI2EJlKeloyg0p34dIw23ZvY9WHq4IORUSkQ4nMNfQhkQnjJEGxA8vGDx0fbDAiIh1ot0ZgZt82s9HRdTOzO82s2sxeTXAa6qxVVlTGsN7D1E8gImkhXtPQ94F3o+vfBMYBo4AfAL9JbljpzcwIl4ZZuF4Dy0Qk9cVLBHXuXhtdPxm4x90/cvengV7JDy29hYpDrKtex4YtG4IORUQkrniJoMHMhplZD2Ay8HTMsZ7JDSv9hUvDABpPICIpL14i+C9gGZHmobnuvhrAzI5Gt492aPzQ8fTI66FEICIpr927htz9UTMrA/q4+ycxh5YBZyc9sjTXLbcbE0ZM0AR0IpLy4o4jcPe62CRgZrPcfbu7b4v3PokIFYd4+f2XqamtCToUEZF2JTKgLFZFUqLIUOHSMHUNdSzdsDToUERE2tXZRPBhUqLIUJOKJwHqMBaR1NZhIjCzsxrXGx9ZGbtP2jegYAAHDjxQA8tEJKUlUiO4KsF90oZQcYhFlYto8IagQxERaVO8R1WeCJwEjDCz38YcKgQ0rWaCQiUh7lxxJ29+9CYHDjww6HBERFqJVyPYSORW0Z3A8phlLnB88kPLDI0DyzTdhIikqnjjCF4xs1XAce5+dxfGlFHGDBhD/579WVS5iAsPuzDocEREWuloHEE9MMDMunVRPBknx3IIlYQ0sExEUlYij5xcByw0s7lEnlcMgLv/OmlRZZhQcYhH33yUzTWbGVgwMOhwRESaSeSuoY3Ao9Fz+8QskqDGfoLFlYsDjkREpLVEnlB2NYCZ9XL37R2dL61VDK8gLyePRZWLOOWAU4IOR0SkmUQGlE0ys9eA16Pb48zslqRHlkEK8gs4bNhh6icQkZSUSNPQTCK3i34EkbuJgC8mMaaMFCoOsXTjUnbX7w46FBGRZhKaa8jdK1vsqk9CLBktXBpmZ91OVry/IuhQRESaSSQRVJpZCHAz62ZmPyTaTCSJC5WEAA0sE5HUk0giuAi4BBgBbADGR7elE4b3GU5ZUZkmoBORlJPIXUObgSldEEvGC5eGeeadZ3B3zCzocEREgMTuGhplZvPMbJOZfWhmj5jZqK4ILtOEikO8t+091lWvCzoUEZEmiTQN/QF4EBgGDAf+CNyXzKAyVfWuagBG/mYk5TPLmbNyTsARiYgklgjM3X8ffX5xnbvfC3iyA8s0c1bOYcZzM5q211WvY9q8aUoGIhK4dhOBmfU3s/7AM2Z2pZmVm1mZmf0IeCyRws3sBDN7w8zWmNmVbRyfYmavRpdFZjZu73+U1DZ9wfRWD7Gvqa1h+oLpAUUkIhIRr7N4OZFv/o29mt+JOebANfEKNrNc4GbgK0AVsNTM5rr7azGnvQMc7e6fRB+EMws4snM/QnpYX72+U/tFRLpKvOcRjPyMZU8A1rj7WgAzux84FWhKBO4eey/lEqD4M14zZZUWlbbZSVxcmLE/soikiXiPqjw93hvd/eEOyh4BxI5IriL+t/0LgcfbiWUaMA2gtLS0g8umphmTZzBt3rRWzUPdcruxafsmBvUaFFBkIpLt4jUNxZsm04GOEkFbN8q32clsZscQSQRfaPNi7rOINBtRUVGRlh3VU8ZGhmJMXzCd9dXrKS0q5esHfp3bl9/OxN9N5PEpjzNmwJiAoxSRbGTuyfm7amaTgJ+7+/HR7asA3P26FucdAvwZONHd3+yo3IqKCl+2bFkSIg7GkqolfO2+r1Hv9fzl7L9wVNlRQYckIhnIzJa7e0Vbx+LdNXRu9PUHbS0JXHcpMNrMRkYfdXkOkQffx16jlEjN4rxEkkAmmlg8kcUXLmZgwUCO/f2xPLDqgaBDEpEsE28cQa/oa592lrjcvQ64FHiCyCR1D7r7ajO7yMwuip72X8AA4BYzW2FmmfNVvxP2678fiy9czJEjjuScP53DDc/fQLJqaiIiLSWtaShZMq1pKNbOup1MfWQq96+6n2mHTePmr95MXk4ij5UWEYkvXtNQvLuGfhuvUHe/7LMGJs31yOvBnNPnMKrvKK59/lrWb1nPg2c+SJ/uekS0iCRPvKah5THL11psL09+aNkpx3KYMXkG/3fK//HU209x1OyjqNpSFXRYIpLBEmoaMrOX3f3QLoinQ5ncNNTSE2ue4Kw/nkVh90LmT5nPIUMOCTokEUlTe3XXUAvp1ZGQIY7f/3iem/ocAF+48ws8seaJgCMSkUyUaCKQgIwbOo4l/7aEkf1G8tU/fJU7Xroj6JBEJMPEG0ew1cy2mNkW4JDG9cb9XRhj1isuLOa5qc9x7Khj+fa8bzN9wXQavCHosEQkQ7SbCNy9j7sXRpe8mPU+7l7YlUEKFHYvZN435/Htw77Ntc9fy7kPn8uuul1BhyUiGUA3qaeR/Nx8bj/5dkb1G8VVC66iaksVfz77zwwoGBB0aCKSxtRHkGbMjCu/cCX3nXEfL2x4gdCdIdZ+sjbosEQkjSkRpKlzPn8OT5/3NJtrNjPxjoksqVoSdEgikqaUCNLYUWVHsehfF9Gnex+OufsYHn69o5nBRURaUyJIcwcMPIAlFy5h/NDxnPngmfzP4v/RhHUi0ilKBBlgUK9B/O1f/sbpnzudHzz5Ay57/DLqG+qDDktE0oQSQYbomd+TB896kCsmXcFNS2/itAdOY/vu7UGHJSJpQIkgg+RYDjcedyM3nXgTj731GEffdTTvb3s/6LBEJMUpEWSgSyZcwiPnPMLrm19n4h0TWf3h6qBDEpEUpkSQoU4eczL/uOAf7KrfRfjOMH97529BhyQiKUqJIIMdPvxwlly4hBGFIzjh3hO455V7gg5JRFKQEkGGK+tbxsJ/XchRZUdx/l/O5+pnr9btpSLSjBJBFujboy+PT3mc88edz8///nOOnn00ZTPLyLk6h/KZ5cxZOSfoEEUkQJp0Lkt0y+3G7FNns3X31mYjkNdVr2PavGkATBk7JajwRCRAqhFkETNj+cbWj5uuqa3hx0/9OICIRCQVKBFkmfXV69vcv2HrBsb87xgunX8pc9+Yy5ZdevaQSLZQ01CWKS0qZV31ulb7+/Xox5gBY5i9YjY3L72ZvJw8JhZP5LhRx3HcfsdRMbyC3JzcACIWkWSzdLuDpKKiwpctWxZ0GGlrzso5TJs3jZramqZ9BfkFzDplFlPGTmFX3S4WVy3mybef5Mm3n+Sl917Ccfr26MvkkZM5br/j+MqorzCy38gAfwoR6SwzW+7uFW0eUyLIPnNWzmH6gumsr15PaVEpMybPaLejeHPNZhasXRBJDGufpGpLFQD799+/qbZwzMhjKOyup5eKpDIlAtkn3J03Pnqjqbbw7LvPsr12O7mWG2lG2m9PM1JejlodRVKJEoEkxe763SyuXNxUW1i+cXlTM9KXR365qcagZiSR4CkRSJf4qOYjFryzoKnGULmlEoD9+u3XVFs4pvwYinoUAZ1rohKRz0aJQLpcYzPSU28/xZNrn+SZd55pakY6svhIhvYaymNvPcau+l1N74nttBaRfUuJQALX2Iz01NqnePLtJ1m6cWmb5/Xv2Z/7zriPsqIySotK6Znfs4sjFclMSgSScnKuzsHp+HdvcK/BlBaVUlZUFln6ljUlibK+ZfTr0Q8z64KIRdJbvESgWzskEO0NbBvRZwR/OOMPrPt0Heuq17Hu03Ws37KeVR+u4rG3HmNn3c5m5/fu1rspQZQWljYlisbXob2HdjgQTn0Vku2UCCQQMybPaHNg2w1fuYEvln0Rylq/x93ZVLMpkhyq1zclinXVkWVJ1RI+3vFxs/fk5+RTXFi8J0HE1CbKispYWLmQS+Zf0hSHJuGTbKREIIFo/CPbmW/iZsbgXoMZ3GswR4w4os1ztu7a2ixJNK1Xr+PptU+zcevGDpukampruPixi3l/6/v07dGXfj370bdH38h6j8h6YffCfT7lhmomEhT1EUhWqa2vpWpLVVOiuOCRC/a6rKLuRXsSRBvJotV6zDm98ns169voaOqPrqJklLnURyASlZ+bz8h+I5sGuf3s2Z+12VdRWlTKyu+u5JMdn/Dpzk+blk927tn+ZMcnfLprz/rbH7/ddM623dvixpGXk9csWaz8cGWr/o+a2hq+N/97bNm5hYL8gqalZ37PPet5PZvt757bfa87z1smoyCbyZSQupYSgWS19voqrp18LYXdCynsXkhZWx0WHahrqKN6Z3XbyaONpNIyCTT6ZOcnXDz/4oSva1hCCaMgr3Viuf7565t9DhBJRlc8cQX799uf7nnd6ZbbrdnSPXfPvrycvH1yB5cSUtfHoaYhyXqp8J+9fGZ5mzWTksISXvz2i+yo3UFNbQ01tTXsqNuzXlNb0+xYy+Pxzm081l4S6izDWieKFskjNnG0Op4T2b7nlXvYuntrq/L79ejHNcdcQ15OHnk5eeTn5jet5+XkkZ/TYjvO8Y7ea2Yp1Vy3L+IIbByBmZ0A/AbIBe5w9+tbHLfo8ZOAGuACd38pXplKBJKJgvyj0+AN7KjdwYE3H9g0u2yswb0GM/vU2eyu383u+t3sqtvVtN60r775vqZzGtp+T7vn1+/mox0fJfXnTUSu5VLv9e0eKy4sJjcnl1zL3SeveTl5kfU2js9eMbvNxFhWVMa7l7+b8M8USB+BmeUCNwNfAaqApWY2191fizntRGB0dDkSuDX6KpJV9uYuqn0lx3Lo1a0X1x97fZvJ6NfH/5qTRp+U9DgatVc7Ki4s5qVpL1HXUEdtQy11DXWR9fo96y2PfZbjM56b0WZ89V7P0eVHU99QT73XJ/S623e3e7yuoS7ue9tKAtD+0wb3RjL7CCYAa9x9LYCZ3Q+cCsQmglOBezxSLVliZn3NbJi7v5eMgL70pS+12veNb3yDiy++mJqaGk46qfUv+wUXXMAFF1zA5s2bOfPMM1sd/+53v8vZZ59NZWUl5513XqvjV1xxBaeccgpvvPEG3/nOd1od/+lPf8qxxx7LihUruPzyy1sdv/baawmFQixatIif/OQnrY7PnDmT8ePH8/TTT/Pf//3frY7ffvvtHHDAAcybN49f/epXrY7//ve/p6SkhAceeIBbb7211fGHHnqIgQMHctddd3HXXXe1Oj5//nwKCgq45ZZbePDBB1sdf/bZZwG48cYbefTRR5sd69mzJ48//jgA11xzDQsWLGh2fMCAAfzpT38C4KqrrmLx4sXNjhcXF3PvvfcCcPnll7NixYpmx8eMGcOsWbMAmDZtGm+++Waz4+PHj2fmzJkAnHvuuVRVNf82PGnSJK677joAzjjjDD76qPk31cmTJ/Of//mfAJx44ons2LGj2fGTTz6ZH/7wh0Div3vl0X8AtX1rYSxd+rtXur2Udz55h111uxhy6hB+9W+/YuTWkW3Gn6zfvZ7be5LzUQ4NZzZAL+BlyHk1h94DenPWI2c1nZfs372yo8oiCelpoHLPse553al/q77Tv3v55Dcd78zvXq9xvaipjibnqXvOKS0qbfUz761kPrN4BM0+Pqqi+zp7DmY2zcyWmdmyTZs27fNARSRiSK8hTCyeyNHlR/Pw2Q8H0jE6pNcQxgwYQ3FhMYYxoGAAYwaMYUivIV0ax4zJMyjIL2i2Lycnp8unVT940MHk5DT/U12QX8CMyW3XWPZG0voIzOws4Hh3/7fo9nnABHf/Xsw5jwHXufvz0e0FwI/cfXl75aqPQES6SircSLCv4ghqHEEVUBKzXQxs3ItzREQCMWXslJQYv5DsOJLZNLQUGG1mI82sG3AOMLfFOXOBf7GIiUB1svoHRESkbUmrEbh7nZldCjxB5PbRO919tZldFD1+GzCfyK2ja4jcPjq1vfJERCQ5kjqy2N3nE/ljH7vvtph1By5JZgwiIhJfMpuGREQkDSgRiIhkOSUCEZEsp0QgIpLl0m72UTPbBLSeiCS9DAQ2Bx1ECtHn0Zw+jz30WTT3WT6PMncf1NaBtEsEmcDMlrU3wi8b6fNoTp/HHvosmkvW56GmIRGRLKdEICKS5ZQIgjEr6ABSjD6P5vR57KHPormkfB7qIxARyXKqEYiIZDklAhGRLKdE0IXMrMTMnjGz181stZl9P+iYgmZmuWb2spk92vHZmS36qNaHzOyf0d+RSUHHFCQz+/fo/5NVZnafmfUIOqauZGZ3mtmHZrYqZl9/M3vKzN6KvvbbF9dSIuhadcAV7v45YCJwiZkdFHBMQfs+8HrQQaSI3wB/dfcDgXFk8ediZiOAy4AKd/88kanszwk2qi53F3BCi31XAgvcfTSwILr9mSkRdCF3f8/dX4qubyXyH73VM5qzhZkVA18F7gg6lqCZWSHwReB3AO6+290/DTSo4OUBPc0sDyggy55e6O7/AD5usftU4O7o+t3A1/fFtZQIAmJm5cChwAsBhxKkmcCPgIaA40gFo4BNwOxoU9kdZtYr6KCC4u4bgBuB9cB7RJ5e+GSwUaWEIY1PcYy+Dt4XhSoRBMDMegN/Ai539y1BxxMEMzsZ+NDdlwcdS4rIAw4DbnX3Q4Ht7KNqfzqKtn2fCowEhgO9zOzcYKPKXEoEXczM8okkgTnu/nDQ8QQoDHzNzN4F7ge+bGb3BhtSoKqAKndvrCE+RCQxZKtjgXfcfZO71wIPA6GAY0oFH5jZMIDo64f7olAlgi5kZkakDfh1d/910PEEyd2vcvdidy8n0gn4N3fP2m987v4+UGlmB0R3TQZeCzCkoK0HJppZQfT/zWSyuPM8xlzg/Oj6+cAj+6LQpD6zWFoJA+cBK81sRXTfT6LPdhb5HjDHzLoBa4GpAccTGHd/wcweAl4icrfdy2TZdBNmdh/wJWCgmVUBPwOuBx40swuJJMuz9sm1NMWEiEh2U9OQiEiWUyIQEclySgQiIllOiUBEJMspEYiIZDklAhGRLKdEICKS5ZQIJC2Z2QAzWxFd3jezDTHb3RIso6+ZXdze9j6MNSnliuwrGlAmac/Mfg5sc/cb2zhmRH7PW81wGp0B9tHofPetthO8drvlt3edrmBmue5e31XXk/SmGoFkHDMrjz7h6xYiUxSUmNkPok+6WmVml0dPvR7YL1qL+GUb25jZuWb2YnTf7dEnqrUqP+bavczsMTN7JXqtsztZ7j/N7G4zezX6tLKCOOW2/Ln/aGa/NrNngKuS9flK5tFcQ5KpDgCmuvvFZnY4kXl7jgQMeMHM/k5kmufPu/t4aPrmHrv9OeBsIOzutdE//FOAf8SW3+K6JwAb3f2r0TKKiDxzojPlXujuC83sTuBiIvPyt1VuS2OJTGh4zN5/bJKNVCOQTLXO3ZdE178A/Nndt7v7NiJTGh+VQBmTgcOBpdFJAicTeYBMy/JjrQSONbMbzOwod6/uZLmV7r4wun5vNPYOy40+z7c/8IsEfi6RZlQjkEy1PWbd9rIMA+5292bNLNGaw/a23uDub0ZrICcB15nZk8A9nSi3Zaedt1euu8f+0T8YeMHd6zrx84kAqhFIdvgH8PXo3Pa9gNOA54CtQJ+Y81puLwDONLPBAGbW38zK4l3IzIYDNe5+L5EmncM6WW6pmU2Krn8TeD5OubHGAq/G/xhE2qYagWQ8d3/JzO4CXozuusPdXwYws4Vmtgp43N3/o43tnwJPmlkOUAtcArwf53JjgV+aWUP0/O+6+0edKPd14Hwzux14C7i1vXLbuO6LiOwF3T4qkiKCuM1UBNQ0JCKS9VQjEBHJcqoRiIhkOSUCEZEsp0QgIpLllAhERLKcEoGISJZTIhARyXJKBCIiWe7/Ax04U9PQuItoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_qubits = 6\n",
    "qreg = QuantumRegister(n_qubits)\n",
    "creg = ClassicalRegister(n_qubits)\n",
    "\n",
    "t = 1\n",
    "J = 1\n",
    "h = 1\n",
    "ks = [i for i in range(1, 11)]\n",
    "\n",
    "hs_distances = []\n",
    "for k in ks:\n",
    "    circuit_approx = ising_unitary_trotter(qreg, creg, t, k, J, h)\n",
    "    circuit_exact = ising_unitary_exact(qreg, creg, t, J, h)\n",
    "    U = execute(circuit_approx, backend=Aer.get_backend('unitary_simulator')).result().get_unitary()\n",
    "    V = execute(circuit_exact, backend=Aer.get_backend('unitary_simulator')).result().get_unitary()\n",
    "    hs_distances.append(hs_distance(U, V))\n",
    "    \n",
    "plt.plot(ks, hs_distances, color='green', marker='o')\n",
    "plt.plot(ks, [0]*len(ks), color='black', linestyle='--')\n",
    "plt.xlabel(r'Trotter steps $r$')\n",
    "plt.ylabel('Hilbert-Schmidt distance')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $F_H$ and $F_Q$ cost functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def uni_cost(n_qubits, n_layers, x0, shots, t, k, J, h, seed_simulator=None):\n",
    "  # calculates the apporximate and exact values of the cost function F_Q\n",
    "    \n",
    "    qreg = QuantumRegister(n_qubits)\n",
    "    creg = ClassicalRegister(n_qubits)\n",
    "        \n",
    "    ansatz = hea(qreg, creg, n_layers, x0)\n",
    "    ansatz_dg = ansatz.inverse()\n",
    "    \n",
    "    u_circ_approx = ising_unitary_trotter(qreg, creg, t, k, J, h)\n",
    "    u_circ_exact = ising_unitary_exact(qreg, creg, t, J, h) \n",
    "    \n",
    "    measurer = QuantumCircuit(qreg, creg)\n",
    "    for i in range(n_qubits):\n",
    "        measurer.measure(i, i)\n",
    "    \n",
    "    circuit_approx = ansatz + u_circ_approx + ansatz_dg + measurer\n",
    "    circuit_exact = ansatz + u_circ_exact + ansatz_dg\n",
    "    \n",
    "    job = execute(circuit_approx, shots=shots, seed_simulator=seed_simulator, backend=Aer.get_backend('qasm_simulator'))   \n",
    "    result = job.result()\n",
    "    counts = result.get_counts()\n",
    "    \n",
    "    zeros_str = '0' * n_qubits\n",
    "    if zeros_str in counts.keys():\n",
    "        p0_approx = counts[zeros_str] / shots\n",
    "    else:\n",
    "        p0_approx = 0\n",
    "        \n",
    "    job = execute(circuit_exact, backend=Aer.get_backend('unitary_simulator'))   \n",
    "    result = job.result()\n",
    "    u = result.get_unitary()\n",
    "    \n",
    "    zerostate = np.array([1] + [0]*(2**n_qubits - 1))\n",
    "    p0_exact = abs( (zerostate.conj().T).dot(u).dot(zerostate) )**2\n",
    "    \n",
    "    gate_count = ansatz.size() + u_circ_approx.size() + ansatz_dg.size() # excluding measurements\n",
    "    \n",
    "    return p0_approx, p0_exact, gate_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def var_cost(n_qubits, n_layers, x0, hamiltonian_dict, hamiltonian_dict_sq, shots, t, J=1, h=1, seed_simulator=None):\n",
    "  # calculates the apporximate and exact values of the cost function F_H\n",
    "    \n",
    "    qreg = QuantumRegister(n_qubits)\n",
    "    creg = ClassicalRegister(n_qubits)\n",
    "    \n",
    "    ansatz = hea(qreg, creg, n_layers, x0) \n",
    "    \n",
    "    gate_count_1 = 0\n",
    "    \n",
    "    \n",
    "  # calculate <H>^2 \n",
    "\n",
    "    circuits = []\n",
    "    multipliers = []\n",
    "    for operator, multiplier in hamiltonian_dict.items():\n",
    "\n",
    "        circuit = QuantumCircuit(qreg, creg)\n",
    "        circuit = circuit + ansatz\n",
    "\n",
    "        gate_count_1 += circuit.size()\n",
    "        \n",
    "        for j in range(n_qubits):\n",
    "            if operator[j] == 'X':\n",
    "                circuit.h(j)\n",
    "                gate_count_1 += 1\n",
    "            if operator[j] == 'Y':\n",
    "                circuit.s(j)\n",
    "                circuit.h(j)\n",
    "            if operator[j] != 'I':\n",
    "                circuit.measure(j, j)\n",
    "                \n",
    "        circuits.append(circuit)\n",
    "        multipliers.append(multiplier)\n",
    "\n",
    "    job = execute(circuits,\n",
    "                  backend=Aer.get_backend('qasm_simulator'),\n",
    "                  shots=shots,\n",
    "                  seed_simulator=seed_simulator)    \n",
    "    result = job.result()\n",
    "    counts = result.get_counts()\n",
    "    \n",
    "    L = len(multipliers)\n",
    "    h1 = 0\n",
    "    for l in range(L):\n",
    "        h1_local = 0\n",
    "        for key in counts[l].keys():\n",
    "            h1_local += counts[l][key] * (-1)**key.count('1') / shots\n",
    "        h1 += h1_local * multipliers[l]\n",
    "        \n",
    "\n",
    "  # calculate <H^2>        \n",
    "\n",
    "    gate_count_2 = 0\n",
    "    circuits = []\n",
    "    multipliers = []\n",
    "    for operator, multiplier in hamiltonian_dict_sq.items():\n",
    "\n",
    "        circuit = QuantumCircuit(qreg, creg)\n",
    "        circuit = circuit + ansatz\n",
    "        gate_count_2 += circuit.size()\n",
    "\n",
    "        for j in range(n_qubits):\n",
    "            if operator[j] == 'X':\n",
    "                circuit.h(j)\n",
    "                gate_count_2 += 1\n",
    "            if operator[j] == 'Y':\n",
    "                circuit.s(j)\n",
    "                circuit.h(j)\n",
    "            if operator[j] != 'I':\n",
    "                circuit.measure(j, j)\n",
    "                \n",
    "        circuits.append(circuit)\n",
    "        multipliers.append(multiplier)\n",
    "        \n",
    "    job = execute(circuits,\n",
    "                  backend=Aer.get_backend('qasm_simulator'),\n",
    "                  shots=shots,\n",
    "                  seed_simulator=seed_simulator)    \n",
    "    result = job.result()\n",
    "    counts = result.get_counts()\n",
    "    \n",
    "    L = len(multipliers)\n",
    "    h2 = 0\n",
    "    for l in range(L):\n",
    "        h2_local = 0\n",
    "        for key in counts[l].keys():\n",
    "            h2_local += counts[l][key] * (-1)**key.count('1') / shots\n",
    "        h2 += h2_local * multipliers[l]\n",
    "        \n",
    "        \n",
    "    # calculate approximate (<H^2> - <H>^2)\n",
    "    fun_value_approx = h2 - h1**2\n",
    "       \n",
    "    # calculate exact (<H^2> - <H>^2)\n",
    "    circuit = QuantumCircuit(qreg, creg)\n",
    "    circuit = circuit + ansatz\n",
    "    job = execute(circuit,\n",
    "                  backend=Aer.get_backend('statevector_simulator'))    \n",
    "    statevector = job.result().get_statevector()\n",
    "    fun_value_exact = ( statevector.conj().T.dot(ham_sq).dot(statevector) - (statevector.conj().T.dot(ham).dot(statevector))**2 ).real\n",
    "        \n",
    "    gate_count = gate_count_1 + gate_count_2    \n",
    "        \n",
    "    return fun_value_approx, fun_value_exact, gate_count"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare $F_H$ and $F_Q$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_qubits = 5\n",
    "n_layers = 4\n",
    "\n",
    "# Hamiltonian and unitary evolution parameters\n",
    "J = 1\n",
    "h = 1\n",
    "t = 1\n",
    "\n",
    "shots = 1000 # number of measurements per circuit\n",
    "ks = [i for i in range(1, 11)] # list of Trotter steps number\n",
    "\n",
    "n_instances = 10 # number of instances to average over\n",
    "seeds = random.sample(range(int(1e10)), n_instances) # simulator seeds for each instance\n",
    "n_pars = 4 * n_qubits * n_layers # number of parameters in HEA\n",
    "x0s = [[random.uniform(-4*pi, 4*pi) for j in range(n_pars)] for i in range(n_instances)] # initial ansatz parameters for each instance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate the maximum possible value of $F_H$ by maximization with a cheaty ansatz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The anzatz is the vector with the components of the form $r_j e^{\\imath \\phi_j}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def constrvec(parameters):\n",
    "    state_vector = []\n",
    "    for i in range(int(len(parameters)/2)):\n",
    "        state_vector.append(parameters[2*i]*np.exp(1j*parameters[2*i+1]))\n",
    "    return np.array(state_vector) / np.linalg.norm(state_vector)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function returns $-F_H$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def maxvar(parameters, ham):\n",
    "    statevector = constrvec(parameters)\n",
    "    t1 = statevector.conj().T.dot(ham).dot(ham).dot(statevector)\n",
    "    t2 = (statevector.conj().T.dot(ham).dot(statevector))**2\n",
    "    return -(t1 - t2).real"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F_H maximal: 39.86345818809669\n"
     ]
    }
   ],
   "source": [
    "hermitian_dict = ising_model_dict(n_qubits, J, h)\n",
    "ham = hamiltonian_dict_to_matrix(hermitian_dict)\n",
    "ham_sq = ham.dot(ham)\n",
    "hermitian_dict_sq = hamiltonian_matrix_to_dict(n_qubits, ham_sq)\n",
    "\n",
    "x0 = [random.uniform(0, 2*pi) for i in range(2**(n_qubits + 1))]\n",
    "res = minimize(maxvar, x0, args=(ham), method='BFGS')\n",
    "\n",
    "# this quantity is the maximal possible value of F_H and is used for normalization of plots\n",
    "normalizer = abs(res.fun)\n",
    "print('F_H maximal:', normalizer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate averaged approximate $F_H$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\kardashin\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:19: DeprecationWarning: The QuantumCircuit.__add__() method is being deprecated.Use the compose() method which is more flexible w.r.t circuit register compatibility.\n",
      "C:\\Users\\kardashin\\anaconda3\\lib\\site-packages\\qiskit\\circuit\\quantumcircuit.py:933: DeprecationWarning: The QuantumCircuit.combine() method is being deprecated. Use the compose() method which is more flexible w.r.t circuit register compatibility.\n",
      "  return self.combine(rhs)\n"
     ]
    }
   ],
   "source": [
    "var_difs = []\n",
    "for i in range(n_instances):\n",
    "    var_val_approx, var_val_exact, var_gate_count = var_cost(n_qubits, n_layers, x0s[i], hermitian_dict, hermitian_dict_sq, shots, t=t, J=J, h=h, seed_simulator=seeds[i])\n",
    "    var_difs.append(abs(var_val_approx - var_val_exact) / normalizer)\n",
    "var_std = sem(var_difs)\n",
    "var_dis = np.var(var_difs)\n",
    "var_dif = np.mean(var_difs) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate averaged approximate $F_Q$ for each Trotter step $r$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k: 10\r"
     ]
    }
   ],
   "source": [
    "uni_difs_k = []\n",
    "uni_diss_k = []\n",
    "uni_stds_k = []\n",
    "uni_gate_counts = []\n",
    "for k in ks:\n",
    "    print('k:', k, end='\\r')\n",
    "    uni_difs = []\n",
    "    for i in range(n_instances):\n",
    "        uni_val_approx, uni_val_exact, uni_gate_count = uni_cost(n_qubits, n_layers, x0s[i], shots, t=t, k=k, J=J, h=h, seed_simulator=seeds[i])\n",
    "        uni_difs.append(abs(uni_val_approx - uni_val_exact))\n",
    "    uni_difs_k.append(np.mean(uni_difs))\n",
    "    uni_diss_k.append(np.var(uni_difs))\n",
    "    uni_stds_k.append(sem(uni_difs))\n",
    "    uni_gate_counts.append(uni_gate_count)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgYAAAETCAYAAABNzM4VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8Z0lEQVR4nO3deXxU9b3/8dcnk4SQRHZZFSgIuFFRcK0KLhSrvdVr93qtXq/lV7eqVFvrtVZrK/W2orTVa7W1WL0Ub28XuymoNe6iRLGAgKIsCgKykwWSTD6/P85kHIYZspjMmZO8n4/HPJL5nu+Zec8QMp98v99zjrk7IiIiIgAFYQcQERGR/KHCQERERJJUGIiIiEiSCgMRERFJUmEgIiIiSYVhB8gH/fr18+HDh7d6v+rqasrKyto/UA5ENXtUc0N0s0c1N0Q3e1RzQ3SzRzU3tD17ZWXlJnfff68N7t7lb+PHj/e2eOqpp9q0Xz6Iavao5naPbvao5naPbvao5naPbvao5nZve3ZggWf4TNRUgoiIiCSpMBAREZEkFQYiIiKSpMJAREREklQYiIiISJIOVxQRkU5lx44dbNy4kfr6+hbv07NnT5YuXdqBqTpOpuxFRUX079+fHj16tPrxVBiIiEinsWPHDjZs2MCQIUPo3r07Ztai/Xbu3Ml+++3Xwek6Rnp2d6e2tpa1a9cCtLo40FSCiIh0Ghs3bmTIkCGUlpa2uCjobMyM0tJShgwZwsaNG1u9v0YM2uCOx99k5pNvBXce+1uL97vytFFcPXl0B6USEZH6+nq6d+8edoy80L1791ZNpzRRYdAGV08ezdWTR1NRUcGkSZP22Db8ur+x6kdnhRNMRES67EhBura+D5pKEBERkSSNGIiISJe2x/RwK3TW6WEVBiIi0qVdPXk0Fx83KONRCV1xelhTCSIiIpKkwkBERCTNrvo4jywMzgPwyMK17KqP5+R5Tz31VMxsr9u7776bk+cHTSWIiIjs4fV3t3H+/fOJxx2A6/+wiO8+spgHLzqWIw7s1aHP/dprr/HDH/6Qiy66aI/2gQMHdujzplJhICIikrCrPs75989nR21Dsq26LhgtOP/++bx8/emUFMU65Lnffvtttm3bxsSJE3NaCKTTVIKIiEjC3CXrkyMF6eJxZ+6S9R323JWVlcRiMY488sgOe46WCG3EwMwuBa4FBgFLgKvc/dl99B8L/Bw4BtgC/AK4xd09pU8xcANwPjAY2AD8xN1/2lGvQ0RE8t/w61p+ltpsquviXDlnIVfOWdii/q09mqGyspJ4PE7//v2TbcOGDWPJkiXJ+w899BB33nkn9fX1lJWV8dOf/pQxY8a06nmaE0phYGZfBGYClwLPJb4+amaHuvuaDP17AI8DzwBHA2OAWUA1cHtK198CBwJTgbeAAYDOjSki0sU19yHddCGiRxau5fo/LEpOH6QqK45x67ljOXvckA7JWFlZyec+9zmmT5+ebEs9vfP999/PQw89xLx58+jTpw8vvfQSn/3sZ3nllVfa9QJQYU0lTANmuft97r7U3a8A3gcuydL/PKAUuMDdF7v774HbgGmWOOejmX0SOB04090fd/dV7j7f3Ss6/NWIiEinMOWwgcRimU8lHIsZUw7ruLn/1157jRNPPJGDDjooeRsyJChCduzYwU033cRvf/tb+vTpA8Bxxx3H4MGDWbRoUbvmyHlhkBjuHw/MS9s0Dzghy27HA8+6e21K21yC6YLhifvnAK8QFAvvmdlbZvZTMytvr+wiItK5lRTFePCiY+nRvZCy4mCRYVlxjB7dC3nwomM7bOHhypUr2bJlC0cddVTG7Q8//DAnn3wyAwYM2KO9uLiY2trajPu0VRhTCf2AGMH8f6oNBH/xZzIQeC9D/6ZtK4ERwInAbuCzQC/gZwTFw+fSH9DMphJMOTBgwAAqKipa9yqAqqqqjPu15bFyLVv2fBfV3BDd7FHNDdHNHtXcEH72nj17snPnzlbvF4/Hk/uN6BXjySuO5cnlm7nukWV891OjOP3gfnQrLGjTY7fEc889h5kxcuTIjM+xcOFCRo8evce2+vp6lixZwogRI7Lm2rVrV6v/PcI8XDF92adlaGuuf2p7QeL7r7j7dgAzuxyYa2YD3H2PQsTd7wXuBZgwYYKnXyWxJTJdXZHH/rZ3Wx7KmD0Copobops9qrkhutmjmhvCz7506dI2zbc3rTFosh/wpeN7cd0jy/jS8SPbMWFmb7zxBqNGjWLw4MEZt/fr1y/IlZJx9uzZTJgwgQMPPDDray4pKWn1UQ5hrDHYBMQJ/tJP1Z+9RxGarM/Sn5R93gfWNhUFCUsTX4e2LaqIiEjHmz59OsuXL8+6/VOf+hS//e1v2bp1KxAsVLz55pv52c9+1u5Zcj5i4O51ZlYJTAZ+l7JpMvD7LLu9CNxmZiXuviul/zpgVeL+88Dnzazc3asSbU2XvVrdXvlFRKRzae7qitkOdczl1RWPO+44rr76aiZNmkRtbS2lpaX88Y9/ZNSoUe0+vRHWVMIM4EEze5ngA/3rBGsB7gEws+nAMe5+WqL/bOB7wCwz+wHBB/51wM0p5zGYDXwX+LWZ3USwxmAm8H/uvjEXL0pERKJnX1dXzCdTp05l6tSpvPPOO3z605+mvLxj1taHUhi4+8Nm1pfgZESDgMUEhxk2/WU/CBiZ0n+7mU0G7gIWAFsJzl8wI6VPlZmdTrDg8JVEnz8RFBAiIiKdwogRI3jjjTc67PFDW3zo7ncDd2fZdmGGtkXAyc085nLgk+2RT0REpCvStRJEREQkSYWBiIiIJKkwEBERkSQVBiIiIpKkwkBERESSVBiIiIhIUpjXShAREQndBz/7OZvuuqvV+/W77DL2v+LyDkgULhUGIiLSpe1/xeWUXHhBxjMfLj34EA5ZtjTDXp2XphJERETyxKmnnoqZ7XV79913c5ZBIwYiIiJ54rXXXuOHP/whF1100R7tAwemX2C446gwEBERSfHmiScR37QpeX/pwYckv4/168fo557tkOd9++232bZtGxMnTsxpIZBOUwkiIiIpUouC1mz7qCorK4nFYhx55JEd9hwtocJAREQkD1RWVhKPx+nfvz/l5eWUl5dz2GGHJbfPnDmT73znO8n77k7v3r3ZtWtXu+bQVIKIiHR6qdMBuXqs1h7NUFlZyec+9zmmT5+ebOvevXvy+8WLF3Paaacl769YsYIBAwZQUlLSqudpjgoDERHp9Jr7kN65c2fycMXmPvg76vDF1157jRtvvJGDDjoo4/ZFixYxbdq05P2FCxdy+OGHt3sOFQYiIiIhW7lyJVu2bOGoo47KuN3dWbp0KWeccQZmBsC2bdu4+uqr2z2LCgMREZGQVVZWYmaMGzcu4/ZVq1bxsY99jIULFybbzjnnnA4ZMdDiQxERkRSxfv3atO2jqKysZNSoURnPvgjB+oLUhYgAb7zxBmPHjm33LBoxEBERSZF6noJcnRJ5+vTpeyw6TLdo0SIOPfTQ5P3du3ezbt06Ro4cSU1NTbtm0YiBiIhInlu8ePEehcGyZcsYPXo0sVis3Z8rtBEDM7sUuBYYBCwBrnL3rKeTMrOxwM+BY4AtwC+AW9zdE9snAU9l2PUQd1/WruFFRKTTaO7qitmOUsjl1RVnz569x/0jjjiCV199tUOeK5TCwMy+CMwELgWeS3x91MwOdfc1Gfr3AB4HngGOBsYAs4Bq4Pa07ocRFA5NPmjv/Jnsqo8zd8l6AB5ZuJYphw2kpKj9KzkREWlf+7q6YlcU1lTCNGCWu9/n7kvd/QrgfeCSLP3PA0qBC9x9sbv/HrgNmGZNx218aKO7r0+5xTvsVSS8/u42jrn1Ca7/wyIArv/DIo659Qlef3dbRz+1iIhIu7LESHzuntCsGKgBvuzuv0tpvws43N0nZtjnN0Bfdz8rpe1o4GVghLuvTJlKWA10A94AfuDumaYXMLOpwFSAAQMGjJ8zZ06rX0tVVRXF3cu46qkaahr23l5aCHeeUkpxLL12CV9VVRXl5eVhx2i1qOaG6GaPam6Ibvao5obws/fs2TPrCYL2JR6Pd8h8fS7sK/uKFSvYvn17xm2nnHJKpbtPSG8PYyqhHxADNqS1bwBOz7LPQOC9DP2btq3kwxGHV4Bi4HzgSTOb5O7PpD+gu98L3AswYcIEnzRpUqtfSEVFBdt7jcIKFgF7D0xYQYzavqP55LghrX7sjlZRUUFbXnPYopobops9qrkhutmjmhvCz7506dI2TQmknvkwavaVvaSkpNUXZQrzcMX0oQrL0NZc/2S7uy8Hlqdsf9HMhgPXEKxN6BCrN9dQU5d5tqKmLs6aze17GImIiOybu7P3LHPX09YZgTDWGGwi+PM6/WLT/dl7FKHJ+iz92cc+APOBUa0N2BrD+pZSWpx5CKe0OMbQvqUd+fQiIpKiqKiI2trasGPkhdraWoqKilq9X84LA3evAyqByWmbJgMvZNntReAkMytJ678OWLWPpxtHMMXQYaYcNpBYljUEsZgx5bD0ekZERDpK//79Wbt2LTU1NW3+iznq3J2amhrWrl1L//79m98hTVhTCTOAB83sZeB54OvAYOAeADObDhzj7k3Xl5wNfA+YZWY/AEYD1wE3p5zH4CqCImEJwRqDfwPOAT7bkS+kpCjGgxcdy/n3zyced6rr4hjBaMGDFx2rQxZFRHKoR48eAKxbt476+voW77dr1652v3xxrmTKXlRUxIABA5LvR2uEUhi4+8Nm1he4geAER4uBM919daLLIGBkSv/tZjYZuAtYAGwlOH/BjJSHLQZ+AgwBagkKhLPc/e8d/HI44sBevHz96cxdsp4r5yzky8cOZdn7O/j4AT07+qlFRCRNjx49Wv2BWFFR0epFevmivbOHdkpkd7/b3Ye7ezd3H5965IC7X+juw9P6L3L3k929xN0HuXtytCCx/b/c/SB37+7ufdz9pFwUBU1KimKcnTj64JazD2fnrgaeWr4xV08vIiLSLnSthA4QKzCumTKGH899k8bGrjnHJSIi0aTCoIN88tABdCss4C//XBd2FBERkRZTYdBBzIxvnTGGGY+/SX28Mew4IiIiLaLCoAOdMLIfQ/uU8vAr74YdRUREpEVUGHSwa6eM4Wf/eIvaLGdHFBERyScqDDrYxw/oxfhhvZn1wqqwo4iIiDRLhUEOTJs8hvuefYftNS0/2YaIiEgYVBjkwEH9y5l8yAB+8czbYUcRERHZJxUGOXLl6aOY/fIaNu7cFXYUERGRrFQY5MjgXt357FEH8PN/rAg7ioiISFYqDHLo0kkj+cvr61izuSbsKCIiIhmpMMihvuXduOCE4dzxxJthRxEREclIhUGOXXzSCJ596wOWrd8RdhQREZG9qDDIsfJuhVwy6SB+Mnd52FFERET2osIgBOcdO5Sl7++kcvWWsKOIiIjsQYVBCEqKYlx52ihue2w57ross4iI5A8VBiE596ghbK7azdNvfhB2FBERkSQVBiEpjBVwzSfH8OO5y2ls1KiBiIjkBxUGITrj8IHECoy/L34/7CgiIiJAiIWBmV1qZivNbJeZVZrZSc30H2tmT5tZrZmtNbMbzcyy9D3RzBrMbHHHpG8fZsa3phzM7fPepD7eGHYcERGRcAoDM/siMBO4FTgSeAF41MyGZunfA3gc2AAcDXwDuBaYlqFvb+A3wJMdEr6dnTiqH4N6lvB/le+FHUVERCS0EYNpwCx3v8/dl7r7FcD7wCVZ+p8HlAIXuPtid/89cBswLcOowa+AB4AXOyh7u7t2yhh++uRb7KqPhx1FRES6uJwXBmZWDIwH5qVtmgeckGW344Fn3b02pW0uMBgYnvLYlwIDgR+0V95cOHJobz5+QE9+8+KqsKOIiEgXZ7k+jt7MBgNrgYnu/kxK+43Aee4+JsM+84D33P2ilLahwGrgBHd/0czGAk8Ax7n7SjO7Cficux+eJcdUYCrAgAEDxs+ZM6fVr6Wqqory8vI92i58rJpZZ5S1+rHWVjXyo5drue2kUkqLMi6daFeZskdBVHNDdLNHNTdEN3tUc0N0s0c1N7Q9+ymnnFLp7hPS2wvbJVXbpFcklqGtuf4AbmbdgDnANe6+skVP7n4vcC/AhAkTfNKkSS3ZbQ8VFRXstd9jf9u7rYVeq32dZV7CtEl71UbtLmP2CIhqbohu9qjmhuhmj2puiG72qOaG9s8exhqDTUCcYMg/VX+CxYWZrM/Sn8Q+g4BDgV8njkZoAG4EDkvc/2S7JO9gV50+it+8tJoPdu4OO4qIiHRROS8M3L0OqAQmp22aTHB0QiYvAieZWUla/3XAKoKpibHAuJTbPcCKxPfZHjevHNC7lHPGDeGup1aEHUVERLqosI5KmAFcaGYXm9khZjaTYCHhPQBmNt3MUg83nA3UALPM7HAzOxe4DpjhgfrE0QrJG7AR2J24X5Xbl9d2l596EH9auJZ3t9SEHUVERLqgUAoDd38YuAq4AVgInAic6e6rE10GASNT+m8nGCEYDCwA7gJuJygwOpV+5d346nHDuPOJt8KOIiIiXVBoZz5097vdfbi7d3P38alHKLj7he4+PK3/Inc/2d1L3H2Qu9/s+zikwt1vynZEQr67+OQRVCzfyJsbdoYdRUREuhhdKyEP9Sgp4usTR/KTucvDjiIiIl2MCoM8df7xw1i0djuvrdkadhQREelCVBjkqZKiGN84bRT/9dhycn0SKhER6bpUGOSxz48/gA07dvHcik1hRxERkS5ChUEeK4wVMO2To/nxXI0aiIhIbqgwyHNnHj6IeKPz2OL1YUcREZEuQIVBnisoML51xsH8eN5yGuKNYccREZFOToVBBJw8qh/7l3fjD6+uDTuKiIh0cioMIsAsGDW484k32VUfDzuOiIh0YioMImL8sN4cOrgnD720uvnOIiIibVQYdoAouuPxN5n5ZOJaBo/9ba/tw6/buw3gytNGcfXk0W1+3munjOG8X77El44ZSnk3/dOJiEj706dLG1w9eTRXTx5NRUUFkyZNytnzjhm4HyeP2p9fPvsOV53e9gJDREQkG00lRMzVk0cz64VVbK7aHXYUERHphFQYRMyBfUr5zBGDubvi7bCjiIhIJ7TPwsDMfmxm30y5f6yZPWRmj5nZL8xsipnFOj6mpLr81IP4/avvsXZbbdhRRESkk2luxOCLwCoAMxsDPAF8Hhie+Poo8Hpim+RI//1KOO/Yocx84s2wo4iISCfTXGHQH3g/8f3lwNvAAe5+sLv3AU4GdgIvmdmIjosp6aaePJInlm5kxcaqsKOIiEgn0lxhsBPolvj+LOBWd/+gaaO7PwecCLwB3NwhCSWjnt2LmHryCGY8vjzsKCIi0ok0VxgsAj5tZscBwwgKgD24exy4A5jS/vFkXy44fjiVq7fyz/e2hR1FREQ6ieYKgxnAVcDzwHrgFjMrztAvDuzXvtGkOd2LY1xx6ih+PFejBiIi0j72WRi4+18JpgpuSnxdA7xmZv9hZsPNrMTMxgHfI8Nowr6Y2aVmttLMdplZpZmd1Ez/sWb2tJnVmtlaM7vRzCxl+0Qze8HMNif6LDOza1qTKYq+ePSBrNlSwwsrNoUdRUREOoFmz3zo7i8CLwKY2VXAD4G70/atBz7T0ic1sy8CM4FLgecSXx81s0PdfU2G/j2Ax4FngKOBMcAsoBq4PdGtCvgpwfRHDfAJ4BdmVuPud7c0W9QUxQqYNnk0t81dzp9G9iWlVhIREWm1Vp3gyAPXExyu+DXgFuBbwGHuPq8VDzUNmOXu97n7Une/guDoh0uy9D8PKAUucPfF7v574DZgWtOogbtXuvscd1/i7ivd/SFgLrDPkYjO4F8+Ppi6hkbmvbEh7CgiIhJxbTrzobu/7+73u/tN7n67u7f4NHyJNQrjgfRCYh5wQpbdjgeedffUM/rMBQYTFCmZnufIxOM93dJsUVVQYHxryhh+Mnc58UYPO46IiESYuef2g8TMBgNrgYnu/kxK+43Aee6+18mSzGwe8J67X5TSNhRYDZyQmO5oan8P2J9gquNmd/9+lhxTgakAAwYMGD9nzpxWv5aqqirKy8tbvV9HcHemv7yLkw8o5MQhRc32z6fsrRHV3BDd7FHNDdHNHtXcEN3sUc0Nbc9+yimnVLr7hL02uHtObwR/5TtwUlr794BlWfaZB/wqrW1Y4nGOS2v/GDCWYKpjC3B+c5nGjx/vbfHUU0+1ab+O8vLKzX7C9Cd9V31Ds33zLXtLRTW3e3SzRzW3e3SzRzW3e3SzRzW3e9uzAws8w2diGBdR2kRweOPAtPb+QLZJ8vVZ+pO+jwfrCxa5+30Eh1ve9JHSRsjRw/swZuB+zJ6/1/pNERGRFsl5YeDudUAlMDlt02TghSy7vQicZGYlaf3XkbiWQxYFfHjmxi7hmk+O4a6n3qZ6d0PYUUREJILCuuzyDOBCM7vYzA4xs5kEUwz3AJjZdDN7MqX/bIJDEGeZ2eFmdi5wHTAjMRyCmV1hZp82s1GJ238A1wAP5fKFhe3QwT34xEF9uf+5lWFHERGRCGr2PAYdwd0fNrO+wA3AIGAxcKa7r050GQSMTOm/3cwmA3cBC4CtBOcvmJHysDGCQxiHAw0EF3y6jkSx0ZVMmzyac+56nn87bhi9yzKdqFJERCSzUAoDAA9OOpTxxEPufmGGtkUEV3PM9nh3Ane2T7poG9a3jDPHDuK/n36b6888JOw4IiISIWFNJUgH+8Zpo/jfBe/y/vba5juLiIgkqDDopAb0KOFLRw/lp0++FXYUERGJkNCmEqTjXTJxJKfcXsHXTqpixP7l3PH4m8xMLRQe+1uLHufK00Zx9eTRHZRSRETyiQqDTqxnaRH/ceLHmPH4m/z8K0dx9eTRyQ/4iooKJk2alOw7/Lq/sepHZ4WUVERE8oWmEjq5f//EcF5euYXFa7eHHUVERCJAhUEnV1pcyOWnHsSP5y4PO4qIiESACoMu4EtHD+WdTVW89M7msKOIiEieU2HQBRQXFjBt8mj+67FlTReaEhERyUiFQRfxmSOGUL07zpNLN4YdRURE8pgKgy4iVmBcO2UMP567nHijRg1ERCQzFQZdyGmH9KesW4w/v7427CgiIpKnVBh0IWbGt884mBmPv0mDRg1ERCQDFQZdzLEj+jKiXzlPv9cQdhQREclDKgy6oGunjOEvb9dTU6fiQERE9qTCoAs6fEhPRvcu4NfPrwo7ioiI5BkVBl3UuaOK+dVzK9leUx92FBERySMqDLqogWUFTDlsAP/99NthRxERkTyiwqAL+8Zpo/jty6v5zYurAHhk4Vp21cfDDSUiIqFSYdCFbdyxm9r6Rr7/lzcAuP4Pizjm1id4/d1t4QYTEZHQqDDoourizvn3z6euoTF5ToPqujg7ahs4//75GjkQEemiQisMzOxSM1tpZrvMrNLMTmqm/1gze9rMas1srZndaGaWsv1cM5tnZh+Y2U4zm29mn+n4VxJNlRvixOOZT3IUjztzl6zPcSIREckHoRQGZvZFYCZwK3Ak8ALwqJkNzdK/B/A4sAE4GvgGcC0wLaXbROAfwFmJx/w78MfmCo6uamNNIzV1mUcFaurirNlck+NEIiKSDwpDet5pwCx3vy9x/wozOwO4BPhOhv7nAaXABe5eCyw2s0OAaWY2wwNXpu1zs5mdBZwDPNshryLC+pcWUFocozpDcVBgUFSoWSYRka7I3HN7znwzKwZqgC+7++9S2u8CDnf3iRn2+Q3Q193PSmk7GngZGOHuK7M811Lgf9z9Bxm2TQWmAgwYMGD8nDlzWv1aqqqqKC8vb/V++WDrjir+82WjJsPJD4sKoFsMJgwo5JxRRfTqlj9FQpTf86hmj2puiG72qOaG6GaPam5oe/ZTTjml0t0npLeHMWLQD4gRTAuk2gCcnmWfgcB7Gfo3bdurMDCzy4ADgAczPaC73wvcCzBhwgSfNGlSC6LvqaKigrbslw8qKir47f8bx/n3zyced6rr4pQVx4jFjAcvOpbhfcv4+VNvcfP897joEx/j4pNG0L04FnbsyL/nUcwe1dwQ3exRzQ3RzR7V3ND+2cP8UzB9qMIytDXXP1M7ZvZZ4MfAee6+us0JO7kjDuzFy9efzq3njgXg1nPH8vL1p3PEgb3oWVrEf551KI9cdiLLNuzk1Nsr+L/K92jUVRlFRDq1MAqDTUCc4C/9VP3ZexShyfos/UnfJ1EUPAh81d3//NGidn4lRTHOHjcEgLPHDaGkaM9RgaF9S7nrK0fx868cxez5q/n0z57jhRWbwogqIiI5kPPCwN3rgEpgctqmyQRHJ2TyInCSmZWk9V8HrGpqMLMvAA8BF7r7/7VXZoHxw3rz+0tO4PJTD+K6PyziP2a9woqNO8OOJSIi7SysqYQZwIVmdrGZHWJmM4HBwD0AZjbdzJ5M6T+bYMHiLDM73MzOBa4DZnhi9aSZfQn4n0T7M2Y2MHHrk8PX1amZGWeOHcTj007m+JF9+cIvXuKGPy1iU9XusKOJiEg7CaUwcPeHgauAG4CFwInAmSnrAQYBI1P6bycYIRgMLADuAm4nKDCafJ1gMeWdwPsptz902AvporoVxrj4pBH845sTKY7FmDzjae56aoXOligi0gmEdR4D3P1u4O4s2y7M0LYIOHkfjzepvbJJy/QqLebGfzmUrx4/jNseW8Zptz/NNVNGc/YRQygosOYfQERE8k7+HKAukTW8Xxn//W/jmfmlcTzwwmrOvut5Xnpnc9ixRESkDVQYSLuZMLwPf7z0BKaePIJrfvc6Fz+wgLc/qAo7loiItIIKA2lXZsa/HDGYJ6ZN5Ojhvfn8PS/yvUcWs1kLFEVEIiG0NQbSuZUUxfh/E0fy+QkH8tMn32LyHc8w9eQRXHjC8L3OldASdzz+JjOffCu489jfWrzflaeN4urJo1v9fCIiXZUKA+lQfcqKuekzhyUXKD744mq+dcYY/uXjg1u1QPHqyaO5evLojKf+HH7d31j1o7My7ygiIq2iqQTJiRH7l/OL8ycw4wtH8KvnVvKvdz/PK6u2hB1LRETSqDCQnDp2RF/+dOknuOjEj3HVnIV8/cFKVm6qDjuWiIgkqDCQnCsoMM4eN4QnvzmRIw7sxWf/+wVu/ssStlbXhR1NRKTLU2EgoSkpinHJpJE8fvXJxBud02Y8zX3PvMPuBp1BUUQkLFp82IXssbIf9lrdP/y6zKv9O3plf9/ybnz/7MP56vHD+dGjy/jNS6v49hkHc9bYQZjpDIoiIrmkwqALaVrZD2Rc3R+2g/qX88sLJvDC25u49e9L+dVzK7nhrEMYP0zXwRIRyRVNJUjeOWFkP/582Ymcf9wwrpj9Gpf9z6us3qwFiiIiuaDCQPJSQYFx7lEH8OQ3J3Ho4B6cc9fz3PLXN6iq82SfXfVxHlm4FoBHFq7V1R1FRNqBphIkr3UvjnHZKQfxhQkHcucTb/Kdl2tY3/0dxh3Yi4seeIV4PCgUrv/DIr77yGIevOhYjjiwV7ihRUQiTCMGEgn779eNH/7rWK47pjvPvrWJL/ziRXbUNlBdF4wSVNfF2VHbwPn3z9fIgYjIR6DCQCJlSHkB5x41hOJY5h/deNyZu2R9jlOJiHQemkqQyFm9uYbdDY0Zt1XXxfnTa2sZ0a+cQwbtR2GWAiJXmjtENBtd/ElEwqLCQCJnWN9SSotjyWmEVN0KC9hV38g3f7eQtVtr+fgBvZgwvDdHDevNUUN707N7UU6z7usQUV38SUTykQqDNB/87OdsuuuuFvUdACxNfN/vssvY/4rLOyxXc1qTGz7MHnZuaP173rugEDvjRigu3Wt7t6ICfv3vR1NSFGN7TT2vvruVV1dv5d6n3+Gf721jSO/ujB/Wh/HDejNhWG+G9S1t80mUutJ7rp/zj07vee7pPW8jd+/yt/Hjx3tLvDHm4D3uP/XUUy3aL2zpud2jm70p98I1W33sTY/5od991Id9+69+6Hcf9bE3PeYL12zN+lj1DXH/57vb/P7n3vHL/qfSj7v1CR9/yzz/2gOv+D0VK3zBqs1eW9fQIblTszcZ9u2/tstztTf9nOee3vPc03vuDizwDJ+JoY0YmNmlwLXAIGAJcJW7P7uP/mOBnwPHAFuAXwC3JF4cZjYIuB04ChgFPOjuF3bka5DwHHFgL16+/nTmLlnPlXMWcuu5Y5ly2EBKimJZ9ymMFTD2gJ6MPaAn//6JjwGwblstC1YHowrf+/MS3t5YzaGDezB+WO/krV95t1y9LBGR0IVSGJjZF4GZwKXAc4mvj5rZoe6+JkP/HsDjwDPA0cAYYBZQTVAMAHQDNgE/Aqa2V9Y3TzyJ+KZNyftLDz4k+X2/Hj3g5fnt9VTtal+5Y/36wQ9uCSNWi7T0PS8pinH2uCFcOWchZ48b0qbnGtyrO5/p1Z3PHDEYgOrdDbz+7jYqV29l9vw1XPO71+lbVsxRw3ozITEFMap/OQUFe08/dIX3PN/oPc89vee5l+v3PKwRg2nALHe/L3H/CjM7A7gE+E6G/ucBpcAF7l4LLDazQ4BpZjYjMSqyCvgGgJl9rr2Cpv5jpIvt2NFeT9Pu9pV7X9vyQZjveVm3Qk44qB8nHNQPgMZG562NVSxYvYUFq7fwi2feZmt1HUcODdYojB/Wm3FDe1FaXNip3/NsF9jKJJdHVHTm9zxf6T3PvVy/5zkvDMysGBgP/CRt0zzghCy7HQ88mygKmswFbgGGAyvbkGMqiZGFAQMGUFFRkbHfgGYeJ9t+YWsud1VVVWSzZ8rd0a9lCDCkP3ymv7F9dzErtu1gyVvbeOSVFazZ2cigsgLuaeYxmt7zurhTuSE4omL67CcYPyBGcSzcq0g2957POqNsj/sXPla9V9uH1lFRsa5dcjWnq/2c5wO957mX6/fcElP0OWNmg4G1wER3fyal/UbgPHcfk2GfecB77n5RSttQYDVwgru/mNb/r8Cmlq4xmDBhgi9YsCDjttQhGxERkXyz4Z7/btPVcs2s0t0npLeHebhiekViGdqa65+pPacOWba0+U4haK6gaesPUi40lz39Pc+X8wE0l/uqf7+dVTst48mZenQv5OXrT9/n4smO1Fnf8870c54v9J7nXq7/QA2jMNgExIGBae39gQ1Z9lmfpT/72Eckr4zoFeOdHZmv41CzO843/3chJ47an37l3ehbXky/suBrWTedbqStLnysWmebFGmlnP/Gcfc6M6sEJgO/S9k0Gfh9lt1eBG4zsxJ335XSfx2wqqOyirSnkpjREM88wNXQ6Hyws47X1mxlU1Udm6t2s6mqjk1Vuykwo295MX3Lu9GvrDhZOPQt70a/8pT7Zd3oXVrUoaeB3lUfT16L4pGFa5s9RDRss84o09kmRVoprD9FZgAPmtnLwPPA14HBEKzfMrPpwDHuflqi/2zge8AsM/sBMBq4DrjZUxZJmNm4xLc9gMbE/Tp3f6OtQWP9+mVd9Rnv0aOtD9vh9pU71q9fjtO0Tmd9z/uXFmQ9lXNZcYzzjhu612GX7k5NXZxNiUJhc9VuNlcHX9/dUsNra7ayuaqOzdW72VxVx/baenp2L0oWCn0ThUO/RCHRtyz4un+imCgtjmFmLXrPX393G+ffPz+vLnWtn/Pc03uee7l+z0MpDNz9YTPrC9xAcIKjxcCZ7r460WUQMDKl/3YzmwzcBSwAthKcv2BG2kO/lnb/XwgWKA5va9bRz314zqWlBx+yxxxUvq5ghX3nBlgXsex7XIwow6Fz2Q6ny+XwcHPv+fgnn+K3b2WeSojFjCmHpc+WgZlR1q2Qsm6FDOub7SiADzXEG9laU58sFFILinfXbGNz9W4+SI5I7Aagb1k3+n35R/QrK6ZveTHMfoAx064Iiomybqxc9k96bqnh/F/NZ8euhuRzNRU4598/P7T1ES39Oc/HkQ79bsk9vectE9rkpbvfDdydZduFGdoWASc385jhHvMlHabpYkTpFyKKkuKY8eBFxyb/6q6ui1NWHCOWaG+PD6rCWAH779eN/fdr2dkaa+oa2FxVxwdVQSGxuWo3Sxt2897WGha+u43NVbtZs2EXty14hqrd2ddHXPzAAg4b0oOy4kJKi2OUFhdS1i3xtThGabfga/fiWNCnW4ziWEGbr1PRGvk40pHNHgVwC9dGgNZHSPvK+eGK+Sj1cMXWXryiSdQuutEk7NwQ3extzV1+6eW8evI5XDlnITO/NC6Uv15bm332mNN58OBPgmVev3Dqwf05engfauoaqN4dp6augZq6+B73q+vi1OwOvtbWxWl0p7Q4Rlm3wg8LhsT90sT97sWxPQqM+LNPE5/3KN3idXRvqKN7w25KGnbTPV5HScNuShrqiKUcqLS7oJB/O+NGqjJccCvXR4J81J/zsNZHRPX/J0Q3e65yZztcUYUB+z6Pwb5E+a/XqGaPam6I7mWXKyoq2N5rFNf/YVHW9RG3nju21aelrmtopLYuTnVdQ0oBsXcR8WExkbifWmzsUYQE3xfFCpLFRU1NLdvqoDHDr7niWAHnHTuUKYcPpE9ZMb1Li+lVWkRRBy7ebKlMP+f5/POyx0hHK+TTSEdn+t3SUvl4HgMRiYgphw3ku48szrgt2/qI5hQXFlBcWEDP0qKPGi/J3dlV3xgUG7vj3PTws/xjTUPGvnXxRp5d8QFL1u1gS00dW6vr2FZbT1lxLCgUyorpU1pMn7LiPe73LiumT1kRvRPbepQUZbx2RnvJx/UR6Zqm+iB6BbCmb/amwkBEmlVSFOvw9RHtwczonljLQDmM6hVj/nrPOtJxxamj9hjpaGx0duyqZ0t1HVtr6thSXc/W6jo2V9expbqOFRur2FpdlywktlTXUV0Xp1f3oj0KifTiIX1bWeJokOZEaX1EuigUNJB5/VJT9jCn+5qz1yhNO56vQ4WBiLRIWy51HbbxA2KtOhKkoMDoVVpMr9LiFj9HfbyRbTX1iUKibo/CYd22XSxeu50tNfXJQmJrTR0Ncad3WRF9yrrtWUAkvr7/fgONSzdyxZxXqU5Z9JkPR4K0RJQLmqhkzzRK014FjQoDEWmx9rjUdS7l4kiQolYeCQLBX6QfFhL1e4xAvP1BFcs2NPCPx5ZSk2GkA2BHbQPH/PAJ9ispoihmFMUKErfg+8K0tsJYAcWxAgoLjKLCAooKmvoVUJzYnrp/02MU7/FYiX0KCiguNAoL9n7O4lgBW2obueJX89mZZ4e2tsSu+jjn3z+fHbXRy96eBY0KA5EO1NxwXz6cf6Gzy8eRjpKiGIN6dmdQz+4Zt1dUbOKf8cHc8fibGbcb8OVjhvJvxw2jodGpjzcmbnt+35DW3hB36uKNifam7526hmBdRn2D09CYeZ9M+zc9V0OjU9/QSH2js6NmN1mObGVHbQMTfvA45d2KkgVHYUFTYWJp3+9Z2BQVGIV7fJ8oWgoK9nqspmKnsODDoiV4vOz7r6tq5KGXVlPfkHlBfjzuPLb4fc458oC2/JN3qLq4t2tBo8JApAPta1GW5E7URjoAhvUtzXqmzNLiGIcO7sGBffY+DDNsV/9yHn9aUZ/x6nYGnHfcMC44fjgNcae+MSg2moqLpoKjoTG1CGkqVhLbE/1S92+IN1LV0LBHwbLH9pT9mwqpYPuHj7V95y5q/C1q6zNXNdV1ca56+HWu+d0/iRUERUwsUWA03W8qQPbYnvyaaI9lbw/aCvbcN1t7cnsBj7+xm1312QuauUvWt+rnXoWBiGTU3GptjXZ0rI44EiQX9nXq79LiGIcO6sHgXplHSsLUksNyf3DO4Xz6iMHEG52GRifeVMSk3U9uT34Nio9M7fXxtH5N2zO018cbqa1P3R4USqt3BqM+mdTUxVmzuaZV74UKAxHJKKpnm+ws0zdRORIkXWsXfOaT5oqxT40dlJjeyHGwZkyv3cBDyxqyFmND+7ZuZEmFgYh0Kp1p+iYf10c0JxcLPjuKirFA+Kf4EhGRrJrWRwCcPW5I3n44pWoqaG49dywAt547lpevPz2vDvfLJorZm4qxHt2DU4dDMPXRo3thmwoajRiIiOSJqK7raG765so5CzPuF3bubKK4WLU9R5dUGIiI5ImoruuI8vSNirG9qTAQEZEuS8XY3rTGQERERJJUGIiIiEiSCgMRERFJUmEgIiIiSaEVBmZ2qZmtNLNdZlZpZic103+smT1tZrVmttbMbrS0C5qb2cTEY+0ys3fM7Osd+ypEREQ6l1AKAzP7IjATuBU4EngBeNTMhmbp3wN4HNgAHA18A7gWmJbS52PA3xOPdSQwHfiZmX22416JiIhI5xLWiME0YJa73+fuS939CuB94JIs/c8DSoEL3H2xu/8euA2YljJq8HVgnbtfkXjM+4AHgGs69qWIiIh0Hjk/j4GZFQPjgZ+kbZoHnJBlt+OBZ929NqVtLnALMBxYmegzL22/ucAFZlbk7vXZMi1fvrxNx4Bu27aNXr16tXq/fBDV7FHNDdHNHtXcEN3sUc0N0c0e1dzQ/tnDOMFRPyBGMC2QagNwepZ9BgLvZejftG1l4usTGfoUJp7z/dQNZjYVmApQVFTEtm3bWvwCmsTj8Tbtlw+imj2quSG62aOaG6KbPaq5IbrZo5ob2j97mGc+9LT7lqGtuf7p7S3pEzS43wvcCzBhwgRfsGDBPsNmErUzZaWKavao5oboZo9qbohu9qjmhuhmj2puaHv2tPX7SWGsMdgExAn+wk/Vn71HEZqsz9KflH2y9WkANrcpqYiISBeT88LA3euASmBy2qbJBEcUZPIicJKZlaT1XwesSumTPhUxGViwr/UFIiIi8qGwjkqYAVxoZheb2SFmNhMYDNwDYGbTzezJlP6zgRpglpkdbmbnAtcBM9y9aZrgHuAAM7sz8ZgXAxey9yJHERERySKUNQbu/rCZ9QVuAAYBi4Ez3X11ossgYGRK/+1mNhm4C1gAbAVuJygwmvqsNLMzgTsIDntcB3wjcWijiIiItEBoiw/d/W7g7izbLszQtgg4uZnHfBo4qj3yiYiIdEW6VoKIiIgkqTAQERGRJBUGIiIikmQfLurvuszsA2B1sx331o/gvAxRFNXsUc0N0c0e1dwQ3exRzQ3RzR7V3ND27MPcff/0RhUGH4GZLXD3CWHnaIuoZo9qbohu9qjmhuhmj2puiG72qOaG9s+uqQQRERFJUmEgIiIiSSoMPpp7ww7wEUQ1e1RzQ3SzRzU3RDd7VHNDdLNHNTe0c3atMRAREZEkjRiIiIhIkgoDERERSVJhICIiIkkqDFrJzE42sz+b2VozczO7MOxMLWFm3zGzV8xsh5l9YGZ/MbPDw87VEmZ2mZn9M5F9h5m9aGZnhZ2rtczs+sTPzM/DztIcM7spkTX1tj7sXC1hZoPM7IHEz/kuM3vDzCaGnas5ZrYqw3vuZva3sLPti5nFzOwWM1uZeL9XmtkPzCy0i/S1hpntZ2Z3mtlqM6s1sxfM7Oiwc6Vq7nPHAjeZ2brEa6gws8Pa+nwqDFqvnOAy0VcCtSFnaY1JBFezPAE4FWgAnjCzPmGGaqH3gG8TXDlzAvAP4E9m9vFQU7WCmR0HfA34Z9hZWmE5wSXQm25jw43TPDPrBTwPGHAWcAhwBbAxxFgtdTR7vt9HAQ78b5ihWuDbwGXAN4CDCX43XgZ8J8xQrfBLYApwAcHP+DyC341DQk21p+Y+d74FfJPgZ/1ogp/3x81sv7Y8mY5K+AjMrAq43N1nhZ2ltcysHNgOnOPufwk7T2uZ2RbgO+7+i7CzNMfMegKvEhQGNwKL3f3ycFPtm5ndBHzO3SMxqtTEzG4FJrr7J8LO8lGZ2X8C1wKD3b0m7DzZmNlfgc3ufkFK2wNAX3f/dHjJmmdm3YGdwGfd/ZGU9krgUXe/IbRwWaR/7piZAeuAn7v7DxNt3QmKg2va8jtSIwZd134E//5bww7SGolhyy8RVNAvhJ2nhe4F/s/d/xF2kFYakRi6XGlmc8xsRNiBWuAcYL6ZPWxmG81soZldnvjlGRmJvP8BPJTPRUHCc8ApZnYwgJkdSjAq+fdQU7VMIRADdqW11wIn5j5Om3wMGEgw0gGAu9cCzxCMELdaJOaApEPMBBYCL4aco0XMbCxB1hKgCvhXd18UbqrmmdnXgIOA88PO0krzgQuBZUB/4AbgBTM7zN03hxmsGSOAS4E7gB8B44CfJbbl/dqOFJMJfuH/MuwgLXAbwR8ab5hZnOBz5Yfufne4sZrn7jvN7EXgBjNbDKwHvgwcD6wINVzLDUx83ZDWvgFo03SICoMuyMxmEFTDJ7p7POw8LbSc4Jd8L+CzwANmNsndF4cZal/MbAxwK3CSu9eFnac13P3R1Ptm9hLwDsE87IxQQrVMAbDA3Zvmt18zs1EEc95RKgy+Brzi7gvDDtICXwS+CnwFWELw/3Smma1091+FGayFzgfuJ1jLFCeY9vstwRqPKElfF2AZ2lpEUwldjJndQVARn+ru74Sdp6Xcvc7dV7h70y/9hcDVIcdqzvEEl0NdbGYNZtYATAQuTdzvFm68lnP3KoJf+qPCztKM94E30tqWAkNDyNImZtYfOBu4L+wsLfRj4CfuPsfdF7n7gwTFYyQWH7r72+4+kWB68kB3PwYoAlaGm6zFmo4WGpjW3p+9RxFaRIVBF2JmMwmq+lPdfVnYeT6iAiDfP1j/RLDKeVzKbQEwJ/F9ZEYRzKyEYMX5+2FnacbzwJi0ttHA6hCytNW/A7sJfk6ioJTgL+1UcSL2+eLu1e7+vpn1JjhK4ZHm9skTKwmKg8lNDYn/ryfRxnVYmkpopcRq/oMSdwuAoWY2Dtji7mtCC9YMM7uLYMjsHGCrmTVVl1WJvwbzlpn9CPgb8C7BXOZXCA6/zOtzGbj7NmBbapuZVRP8rOTtFAiAmf0E+AuwhuAvj+8CZcADYeZqgTsI1kL8J/AwcCTBYXTXh5qqhRKLDi8G5rj7zrDztNBfgOvMbCXBqNKRwDTgN6GmaiEzm0Lwu3wZwe/2HxNMXf46zFypmvvcMbM7gf80s2XAmwRrgqqA2W16QnfXrRU3gg8kz3CbFXa2ZnJnyuzATWFna0H2WQR/8e0mOATnCWBK2Lna+FoqCA4rCj1LMznnEBwCVQesBX4PHBp2rhZmPwt4nWCl+ZsEhYGFnauF2U9J/L88Juwsrci8H3Bn4v9oLcFalFuBkrCztTD/F4C3E79f3idYi9Iz7FxpGff5uUOwnuCmRP5dwNPA4W19Pp3HQERERJIiNQckIiIiHUuFgYiIiCSpMBAREZEkFQYiIiKSpMJAREREklQYiIiISJIKAxEREUlSYSAiIiJJKgxE8pSZeQtuq9rpuc4xs2ktbQ9LvuUR6Yx05kORPGVmx6U1/ZHgVL83pbTtdvfX2uG5ZgGnu/sBLWkPS77lEemMdBElkTzl7i+l3jez3cCm9PZMzKybu+/usHAfQT5n6wiJCyMVunt92FlEWkJTCSIRZ2Y3JaYVDjezuWZWBfxvyvYzzOxFM6s1s+1m9iczG5OyfRZwATAkdYoiW3vKfkeY2Z/NbGvisZ83s5Naky3DaxltZn80s41mtsvM1pjZ78yssJ3zjDWzp8ysxszeN7Pvm1lBS3K04N/jdTP7tZldYmZLgXrgjOb2E8kXGjEQ6TweAX4F3AY0QlAUEFyy+h/AF4Fy4PvAc2Y2zt3XArcA+wNHA59JPNZuYEeWdszsKOBZ4DXga0AN8HXgCTM7wd0rm8uWxV8JLlV9CbAJGAKcSfBHTLacbcnzJ+B+YDowheCy0o18OE2zrxxZmVkxcEgi58jE41YBr+xrP5F8osJApPP4qbvPTGv7AcFlcD/l7g0AZvYiweWIvwlMc/e3zewDoC7D9EXGdoJr1q8BTnX3ukTfucBigg/Dc1qQbQ9m1g8YBZzt7n9O2dR0TfmsOduQ5z53/1Hi+3lm1gP4ZuK69oXN5NiXw4Eigvf3NHePt2AfkbyiqQSRzuOPqXfMrAw4Cni4qSgAcPeVwPPAxLY8iZl1T+z7O6AxMcxfSHBN+CeAk5vLlsVmgiLmR2b2NTMb1YF50qcz5hCMphze1hwJRya+XqOiQKJKhYFI5/F+2v3eBB+O6e0A64E+bXyePkCM4C/x+rTb5UDv1Pn6LNn24sEhUpOBBQRD/G+a2TtmdkkH5NmQ5f6Qj5ADgsJgnbsvaEFfkbykqQSRziP92OOtibaBGfoOJPjLuC22EczH3wX8JmMQ9/R1BC06Ltrd3wG+mljJfwTBB/vdZrbK3R9txzwDCEYFUu8DrP0IOSAoDLSeQCJNIwYinZS7VwOVwOfNLNbUbmbDgBOAp1O67wa6Z3iYvdoTj/sswQfmq+6+IP3WDtnd3RcCTSczOryd83wh7f6XCBYJLm5hjr0kRiU+TvCei0SWRgxEOrfvEhyV8Fczu5tgHv1mYDtwe0q/N4A+ieHyBcAud1+0j/ZpwDPAXDP7FcFUQT+CNQ0xd7+utUHN7OPATOBhYAXB9MCFQAPBURX7ytnaPF9LfJC/QnBUwsXATe6+rYU5MhlF8P6+2trXLpJPVBiIdGLu/piZnQV8j2DBXR1QAXzL3deldP0lcBxwK9ALWA0Mz9bu7q+a2dGJx/0p0BP4gOBD8Z42xl1PcGTBNOAAYBewCPh0yuGG7ZXnbOBnBIXTdoKjN25pRY5MmhYeqjCQSNMpkUWkyzCzmwiKh6LUIzVE5ENaYyAiIiJJKgxEREQkSVMJIiIikqQRAxEREUlSYSAiIiJJKgxEREQkSYWBiIiIJKkwEBERkSQVBiIiIpKkwkBERESS/j9iT1ACrwFpOQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgoAAAETCAYAAABX8C6bAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1WUlEQVR4nO3deZwdZZ3+/c+VJpB9gZAEgmwSMQKCJI6EYUkiCA7jyCPqKAoqCIwKvxlkfiqoYxwXGHEw4OgzgCAaHgw6jrIpYUsrstOsgQACYQ1JgJClk87W+T5/VHVSOTnVXafT3aerc71fr/NKn6q7qq5zaPr+VtVdVYoIzMzMzKrpV+8AZmZm1nu5UDAzM7NcLhTMzMwslwsFMzMzy+VCwczMzHJtV+8AvdGoUaNizz33rHm5lStXMnjw4K4P1M3KmhvKm72suaG82cuaG8qbvay5obzZO5u7qanpjYjYudo8FwpV7Lnnnjz44IM1L9fY2MiUKVO6PlA3K2tuKG/2suaG8mYva24ob/ay5obyZu9sbkkv5s3zqQczMzPL5ULBzMzMcrlQMDMzs1wuFMzMzCyXCwUzMzPL5aseavDMYYfT+sYbVec1jBoF3/1ODycqrr3so4YNg/vv6+FExfg773n+znuev/Oe5++8OB9RqEHef5iO5vUG7eVrWL68B5PUxt95z/N33vP8nfc8f+fFuVAwMzOzXIqIemfodSZNmhTVbrg0750T6pDGzMysuAlPzat5GUlNETGp2jyPUehCi/77/+21d/LqqMjpzC9WT+got7/zrufvvOf5O+95ffk772o+9WBmZma5XCiYmZlZLhcKNWgYNapT83qD9vK1DhvWg0lq4++85/k773n+znuev/PiPEahBu/4y53tzl/Q2NgzQTqhveyNJc0N/s67g7/znufvvOf5Oy/ORxTMzMwslwsFMzMzy9WjhYKkL0l6TNLy9HWPpOMy86+SFBWveyvWsYOkH0t6Q9JKSddL2q2izUhJMyUtS18zJY3ooY9pZmbWZ/T0EYVXgK8CBwOTgDuA30t6d6bNbcAumdffVaxjBnAC8EngcGAYcKOkhkyba9JtfBA4Nv15Zhd/FjMzsz6vRwczRsR1FZO+LukLwGTgsXTamohYWG15ScOBU4HPRcSt6bSTgBeBo4DZkiaQFAeHRcTdaZszgDsl7RsRT3f15zIzM+ur6jZGQVKDpE8AQ4C7M7MOk7RY0jOSLpc0OjNvItAfuKVtQkS8DMwDDk0nTQaaK9Z5F7Ay08bMzMwK6PFnPUg6ALgHGEDSoX8qIm5K530CWAXMB/YEvgs0ABMjYo2kE4FfAv0jE1zSHcBfI+IMSecBn4+IvSu2+zxweUScn5PrdOB0gDFjxkycNWtWzZ+tubmZIUOG1LxcvZU1N5Q3e1lzQ3mzlzU3lDd7WXNDebN3NvfUqVNzn/VARPToC9ge2IdkjML5wBvA/jltdwXWAR9J358IrCctcDLt5gD/nf58HvBclXXNB75WJOPEiROjM+bMmdOp5eqtrLkjypu9rLkjypu9rLkjypu9rLkjypu9s7mBByOnT+zxUw8RsTYino2IByPiXOAR4OyctgtIBkCOTyctJDnCUHlbqtHAokyb0ZLUNjP9eedMGzMzMyugN9xHoR+wQ7UZkkYB44DX0klNJEcYjs602Q2YwKYxCfeQjHuYnFnVZGAwm49bMDMzsw706FUPki4AbgJeBoaSnEqYAhwnaQgwHfgtSWGwJ8mpicXA7wAiYpmkK4ALJS0G3gQuIrli4ra0zTxJNwOXSjoNEHApcGP4igczM7Oa9PSzHsYCV6f/LiPp4D8YEbMlDQQOAE4GRpAUC3OAj0fEisw6ziYZp3AtMBC4HTg5IlozbT4FXMKmqyOuB87sps9kZmbWZ/X0fRQ+2868FuCYAutYDZyVvvLaLAE+3YmIZmZmltEbxiiYmZlZL+VCwczMzHK5UDAzM7NcLhTMzMwslwsFMzMzy+VCwczMzHK5UDAzM7NcLhTMzMwslwsFMzMzy+VCwczMzHK5UDAzM7NcLhTMzMwslwsFMzMzy+VCwczMzHK5UDAzM7NcLhTMzMwsV48WCpK+JOkxScvT1z2SjsvMl6TpkhZIapHUKGm/inXsIOnHkt6QtFLS9ZJ2q2gzUtJMScvS10xJI3roY5qZmfUZPX1E4RXgq8DBwCTgDuD3kt6dzv8KcA5wFvBeYDFwq6ShmXXMAE4APgkcDgwDbpTUkGlzTbqNDwLHpj/P7J6PZGZm1ndt15Mbi4jrKiZ9XdIXgMmSHgf+BbggIn4LIOkzJMXCicClkoYDpwKfi4hb0zYnAS8CRwGzJU0gKQ4Oi4i70zZnAHdK2jcinu7uz2lmZtZX1G2MgqQGSZ8AhgB3A3sBY4Fb2tpERAvwZ+DQdNJEoH9Fm5eBeZk2k4HmdJ1t7gJWZtqYmZlZAYqInt2gdABwDzCApEP/VETcJOlQkg59j4h4KdP+SmBcRBwj6UTgl0D/yASXdAfw14g4Q9J5wOcjYu+K7T4PXB4R5+fkOh04HWDMmDETZ82aVfNna25uZsiQITUvV29lzQ3lzV7W3FDe7GXNDeXNXtbcUN7snc09derUpoiYVG1ej556SD0NHASMIBlr8AtJUzLzKysXVZlWqbJNtfbtriciLgMuA5g0aVJMmTIlr2muxsZGOrNcvZU1N5Q3e1lzQ3mzlzU3lDd7WXNDebN3R+4eP/UQEWsj4tmIeDAizgUeAc4GFqZNxlYsMhpYlP68EGgARnXQZrQktc1Mf94508bMzMwK6A33UegH7ADMJ+nkj26bIWkAyZUNbeMNmoB1FW12AyZk2txDMu5hcmYbk4HBbD5uwczMzDrQo6ceJF0A3AS8DAwluZphCnBcRISkGSRXQjwFPAN8g2QcwzUAEbFM0hXAhZIWA28CFwGPAbelbeZJupnkKonTSE45XArc6CsezMzMatPTYxTGAlen/y4j6eA/GBGz0/k/AAYCPwFGAvcBH4iIFZl1nA2sB65N294OnBwRrZk2nwIuYdPVEdcDZ3bHBzIzM+vLevo+Cp/tYH4A09NXXpvVJDdkOqudNkuAT3cmo5mZmW3SG8YomJmZWS/lQsHMzMxyuVAwMzOzXC4UzMzMLJcLBTMzM8vlQsHMzMxyuVAwMzOzXC4UzMzMLJcLBTMzM8vlQsHMzMxyFSoUJJ0saaeceTtKOrlrY5mZmVlvUPSIws+Bt+fM2yudb2ZmZn1M0UJB7cwbTPI0RzMzM+tjcp8eKekg4ODMpA9J2r+i2UDgE8Bfuz6amZmZ1Vt7j5n+MPCt9OcAvp7T7k3g1K4MZWZmZr1De6ceZpCMP9ib5NTDR9L32deuwOiIuL7IxiSdK+kBScslvS7phsqjFJKukhQVr3sr2uwg6ceS3pC0UtL1knaraDNS0kxJy9LXTEkjiuQ0MzOzRO4RhYhYBiwDkLQX8FpErN3K7U0Bfgo8QFJ8/Dtwm6R3RcSSTLvbgJMy7yu3O4PkiMcnSY5oXATcKGliRLSmba4Bdgc+SHJE5GfATOBDW/kZzMzMthntnXrYKCJebPtZ0mhgQJU2LxVYzzHZ95JOIilG/ha4ITNrTUQsrLYOScNJTnV8LiJuzaznReAoYLakCcCxwGERcXfa5gzgTkn7RsTTHWU1MzOz4vdRGCbp55JWAa8B86u8OmNomuGtiumHSVos6RlJl6fFSZuJQH/glrYJEfEyMA84NJ00GWgG7s4sdxewMtPGzMzMOqCI6LiRNBM4AbgCeBxYU9kmIn5R88alXwPjgUltpwwkfQJYRVJ87Al8F2gAJkbEGkknAr8E+kcmvKQ7gL9GxBmSzgM+HxF7V2zveeDyiDi/SpbTgdMBxowZM3HWrFm1fhyam5sZMmRIzcvVW1lzQ3mzlzU3lDd7WXNDebOXNTeUN3tnc0+dOrUpIiZVnRkRHb6AxcCXirQt+iIZV7AA2LuDdrsC64CPpO9PJLlvgyrazQH+O/35POC5KuuaD3yto2wTJ06MzpgzZ06nlqu3suaOKG/2suaOKG/2suaOKG/2suaOKG/2zuYGHoycPrGWGy512Xl9ST8iGYg4LSKeb69tRCwAXiE58gCwkOQIw6iKpqOBRZk2oyVtvFFU+vPOmTZmZmbWgaKFwiy66GoBSReTHBWYFhFPFWg/ChhHMjYCoInkCMPRmTa7ARPYNCbhHmAIyViFNpNJ7iKZHbdgZmZm7Sh01QPJwMEZkoYCfwCWVDaIiDs6Womkn5Bc9ng88Jaksems5oholjQEmA78lqQw2BM4n+TUx+/S7SyTdAVwoaTFbLo88jGSyyqJiHmSbgYulXQayRGRS4Ebw1c8mJmZFVa0ULgu/Xcv4LOZ6UHSCQfJ6YCOfDH99/aK6d8mKRBagQOAk4ERJMXCHODjEbEi0/5sknEK15LcRvp24OTYdA8FgE8Bl7Dp6ojrgTMLZDQzM7NU0UJhaldsLCLae7gUEdECHNNem7TdauCs9JXXZgnw6VozmpmZ2SZFb7j0p+4OYmZmZr1P0SMKZmZmpbN8+XIWL17MunXralpu+PDhzJs3r5tSdZ9qufv378/o0aMZNmxYp9ZZqFBIb2bUnoiI93cqgZmZWTdYvnw5ixYtYty4cQwcOJDMFfMdWrFiBUOHDu3GdN2jMndE0NLSwquvvgrQqWKh6OWR/UgGLWZfo0ie0fCO9L2ZmVmvsXjxYsaNG8egQYNqKhL6EkkMGjSIcePGsXjx4k6to+gYhSk5Ad4O/B74fqe2bmZm1k3WrVvHwIED6x2jVxg4cGDNp1/aFD2iUFVEPAdcAFy4NesxMzPrDtvqkYRKW/M9bFWhkHqd5PSDmZmZ9TFbVShI2hH4MvBc18QxMzOz3qToVQ/zSe6+mLU9MCb9+YSuDGVmZma9Q9H7KPyJLQuF1cCLwG/SsQpmZmbWxxS96uGz3ZzDzMysV1q9rpXZTyzkxTdXscdOgzhmv7EM6F/k8UZbZ9q0acyZM2eL6S+99BJve9vbun37bWq+M2P6hMeRwJKIWNn1kczMzHqHR19eyklX3kdra7BqbSuDtm/gm9fNZeYp7+PAt43o1m0//PDDfO973+OUU07ZbPrYsWNzlugehQczSjpG0oPAUuAFYJmk+yUd3U3ZzMzM6mb1ulZOuvI+lresZ+XaVgJYubaV5S3rOenK+1i9rrXDdXTWc889x9KlSznyyCMZO3bsZq+eVqhQkHQMcBMwBPgOyeOivwsMBf7gYsHMzPqa2U8spLW1cnheorU1mP3Ewm7bdlNTEw0NDbznPe/ptm0UVfTUw3TgFuDvI2JD20RJ/w7cCHwbuLXL05mZmXWxPb9201avY+XaVv551iP886xHCrV/4YLjalp/U1MTra2tjB49euO0PfbYgyeeeGLj+6uvvpoZM2awbt06Bg8ezCWXXMK+++5b03aKKFooHAh8LFskAETEBkk/BX7d5cnMzMy6QZFOe8WKFdzx3HLO+9/HWbl2y1MMg7dv4PsfOYAPHzSuOyLS1NTERz/6Uc4///yN07K3o77yyiu5+uqrueWWW9hxxx259957OeGEE3jggQe6/GFWRccorAHyHjk1NJ3fIUnnSnpA0nJJr0u6QdL+FW0kabqkBZJaJDVK2q+izQ6SfizpDUkrJV0vabeKNiMlzZS0LH3NlDSi4Oc1M7Nt3DH7jaWhofqtjxsaxDH7dd94gYcffpjDDjuMffbZZ+Nr3LikKFm+fDnTp0/nV7/6FTvuuCMAhxxyCLvuuiuPP/54l2cpWig0At+RtFd2oqTdSU5LbHn9RnVTgJ8ChwLTgPXAbekdHtt8BTgHOAt4L7AYuFVStkSaQXKTp08Ch5MUMTdKyl6vcg1wMPBB4Nj055kFc5qZ2TZuQP8GZp7yPoYN3I7B2zcgkiMJwwZux8xT3tdtl0jOnz+fJUuWcPDBB1edf+2113LEEUcwZsyYzaZvv/32tLS0dHmeoqcevgrcBTwt6V7gNWAscAjJVRBfLbKSiDgm+17SScAyksdV36DkqRX/AlwQEb9N23yGpFg4EbhU0nDgVOBzEXFrZj0vAkcBsyVNICkODouIu9M2ZwB3Sto3Ip4u+LnNzGwbduDbRnD/eUcx+4mFvPTmKnbvgfsoNDU1IYmDDjqo6vwnn3ySAw44YLNp69at48knn2T8+PFdnqfoDZeekfRukj39w0n2zpcAFwM/iojXOrn9oSRHNd5K3+9FUoDcktl2i6Q/kxyFuBSYCPSvaPOypHlpm9nAZKAZuDuzrbuAlWkbFwpmZlbIgP4N3TYWoZqmpibGjx+fO9Zg2LBhrFmz+Rn/3/zmN0ycOJFdd921y/MoovqlHz1B0q+B8cCkiGiVdChJh75HRLyUaXclMC4ijpF0IvBLoH9kwku6A/hrRJwh6Tzg8xGxd8X2ngcuj4jzqSDpdOB0gDFjxkycNWtWzZ+nubmZIUOG1LxcvZU1N5Q3e1lzQ3mzlzU3lDd7vXMPHz6cffbZp1PLtra20tDQ/Xdf7Iz777+fL37xi9x6662MHDmShx9+mFNPPZVf//rX7LXXXrm5n332WZYtW1Z13tSpU5siYlK1eUUfCvUOYJeI+FOVeUcAr0XEX4usK7PcRcBhJKcHKoeUVlYvqjJti1VWtKnWPnc9EXEZcBnApEmTYsqUKR1sbkuNjY10Zrl6K2tuKG/2suaG8mYva24ob/Z65543b16nrwBYsWJFl1890FXe//73c8455/ChD32IlpYWBg0axHXXXce73vWudnMPGDCgU/dlKDqYcQbwoZx5fw/8qJaNSvoRyUDEaRHxfGZW290rKoeSjgYWZdo0AKM6aDM6HfPQtk0BO2famJmZldLpp5/Oo48+ys0338zatWu79chN0UJhEvDnnHl/Jrk6oRBJF5MMTJwWEU9VzJ5P0skfnWk/gGRcRNt4gyZgXUWb3YAJmTb3kNxFcnJm3ZOBwWw+bsHMzKy09t57b5588kl23333bttG0asehpI8VrqadcDwIiuR9BPgJOB44C1JbUcOmiOiOSJC0gzg65KeAp4BvkEyMPEagIhYJukK4EJJi4E3gYuAx4Db0jbzJN1McpXEaSSnHC4FbvQVD2ZmZsUVLRSeB95P5kqDjGkkD4kq4ovpv7dXTP82yf0YAH4ADAR+QvKUyvuAD0TEikz7s0nuwXBt2vZ24OSKsQ6fAi7JZL4eOLNgTjMzM6N4ofBLkhsuvQT8LCLWSNoB+DzJfQ+mF1lJRFS/xdXmbSJdX+46I2I1yQ2ZzmqnzRLg00VymZmZWXVFC4UfkoxD+DFwsaQlwI4kYxx+C/xH98QzMzOzeip6w6VW4KOSppEMItwJeAO4JSIauy+emZmZ1VPRIwoARMQdwB3dlMXMzMx6maKXR5qZmdk2yIWCmZmZ5XKhYGZmZrlcKJiZmfVC06ZNQ9IWr5dffrlHc9Q0mNHMzMx6xsMPP8z3vvc9TjnllM2mjx1b+Tik7lVToSDp3cARJJdHXhoRCyXtAyyquHOimZlZ6T1z2OG0vvFG1XkNo0bxjr/c2S3bfe6551i6dClHHnlkjxcGlYo+ZnoH4GrgI2x6VPMNJA9w+gHJMxm+1k0ZzczM6iKvSOho3tZqamqioaGhU4+F7mpFxyh8DziK5IFOY0iKhTZ/BI7p4lxmZmbbrKamJlpbWxk9ejRDhgxhyJAh7LfffhvnX3zxxZx77rkb30cEI0eOZPXqvOc3dl7RUw+fBL4REddIaqiYNx/Ys0tTmZmZdZN575zQ4+ua8NS8mtbb1NTERz/6Uc4///yN0wYOHLjx57lz5/L+979/4/tnn32WMWPGMGDAgJq2U0TRQmEnIO9T9gN26Jo4ZmZm3atIp71ixQqGDh3aYSFQawFQ1MMPP8y//du/sc8++1Sd//jjj/PlL3954/tHHnmE/fffv1uyFC0U5gOTqX775r8Bnu6yRGZmZtuw+fPns2TJEg4++OCq8yOCefPmceyxxyIlIwGWLl3K2Wef3S15annM9HmSXgD+N50WkqYCZ1PwMdNmZmbWvqamJiRx0EEHVZ3/wgsvsNdee/HII49snHb88cd32xGFooMZfwDcBMwElqTT/gLcBtwcET/uhmxmZmZ11TBqVKfmbY2mpibGjx/P0KFDq86fO3fuZgMbAZ588kkOOOCAbslTy2OmPyHpJyRXOIwG3iQpEv5UywYlHQH8KzAR2BX4XERclZl/FfCZisXui4hDMm12AH5IMshyIHA78MWIeCXTZiRwCfAP6aTrgbMiYmktec3MbNvVXfdJaM/555+/2SDGSo8//jjvete7Nr5fs2YNCxYs4O1vfzurVq3q8jxF76OwO/BaRNwJ3Fkxbztg14h4qeA2hwBzSU5n/DKnzW0kl2K2WVsxfwbwYZJC4U3gIuBGSRPTogbgGmB34IMk9334GckRkQ8VzGlmZtbrzJ07l4997GMb3z/11FO84x3voKGh8qLErlHrYMb7q8w7MJ1eKGFE/AH4A2w8elDNmohYWG2GpOHAqSRHIm5Np50EvEhyr4fZkiYAxwKHRcTdaZszgDsl7RsRHnxpZmaldM0112z2/sADD+Shhx7qtu0pIjpuJG0ADomILQoFSYcAd0ZE/5o3LjUDZ1Y59XA8yVGEpcCfgK9HxOJ0/jSSUw2jI+L1zHJPAP8TEd+SdApwMTAs0g+oZGjoCpLTDz+vkuV04HSAMWPGTJw1a1atH4fm5maGDBlS83L1VtbcUN7sZc0N5c1e1txQ3uz1zj18+PDcyws70tra2m176N2pvdzPPvssy5Ytqzpv6tSpTRExqdq83CMKkkYAO2YmjZO0d0WzgSTjCaru/XfSzSRXVrTdyOm7wB3paYU1wFigFai8d+aidB7pv69HpgqKiJC0ONNmMxFxGXAZwKRJk2LKlCk1B29sbKQzy9VbWXNDebOXNTeUN3tZc0N5s9c797x583IHBHak7T4KZdNe7gEDBnTqltDtnXr4Z+BbJOf3A/ifnHZK23WJiMjuyj8uqYnktMJxbLo0My9H9vBItUMllW3MzMysHe0VCr8HXiDpXK8k2bN/rqLNGuDJiHisO8IBRMQCSa8A49NJC0nGQ4wCXs80HQ38OdNmtCRVnHrYmeTIg5mZmRWQWyhExKPAowCSArgxIt7sqWBtJI0CxgGvpZOagHXA0SRXNiBpN2ACcHfa5h6SqysmZ6ZNBgZn3puZWR8XERvvXrgtKzIeMU/R+yj8otNbqCBpCNA2uqQfsLukg0hu5LSE5C6PvyUpDPYEzgcWA79LsyyTdAVwYTrmoO3yyMdILqskIuZJuhm4VNJpJEdFLiUpdnzFg5nZNqB///60tLQwaNCgekepu5aWFvr3r/maA6D45ZFI2p/kssR9gcrHU0VEvH/LpaqaBMzJvP92+voF8AXgAOBkYARJsTAH+HhErMgsczawHriWTTdcOjlzDwWAT5HccOmW9P31wJkFM5qZWcmNHj2aV199lXHjxjFw4MBt8shCRNDS0sKrr77KmDFjOrWOojdceh/JZYovkIwVeAwYSXJDo1eAZ4tuMCIaSfbw8xxTYB2rgbPSV16bJcCni+YyM7O+ZdiwYQAsWLCAdevW1bTs6tWru+WRzd2tWu7+/fszZsyYjd9HrYoeUfg+yRUHJ5GMDzg1Ih5K72kwk2Sgo5mZWa8ybNiwTnWQjY2NnbqUsN66I3fRh0K9G7iaTZcWNgBExB0kRUL+TanNzMystIoWCv2BlRGxgWTA4S6ZeU8D3fNsSzMzM6urooXCcySXKEIyPuEUSf0k9QM+R9femdHMzMx6iaJjFG4AppDct+D7wE3AcpJbKQ8B/k93hDMzM7P6KnofhemZn29LHwT1UZJLE2+OiFvyljUzM7PyKnwfhayIeBh4uIuzmJmZWS9Tc6GQjkvIitiae0OamZlZr5U7mFHSWEk3STo5M62B5D4K2ddSSZ273ZOZmZn1au0dUfgicDDJWIQsAZcDC9Kf/xH4J5LbMJuZmVkf0l6hcCxweUS0VEwP4NKIeAhA0uskz2ZwoWBmZtbHtHcfhX2p/kjmyuc0PJO2NTMzsz6mvSMKA4Dm7ISIaJW0C/BGZvJqtnyapJmZmfUB7R1RWAzsXTkxIhZVPM55L+D1rg5mZmZm9ddeofAXkqdFduRk4K6uiWNmZma9SXuFwiXANEk/lLTFKQpJ20m6iOTWzhcX3aCkIyRdL+lVSSHpsxXzJWm6pAWSWiQ1Stqvos0Okn4s6Q1JK9P17VbRZqSkmZKWpa+ZkkYUzWlmZmbtFAoRcQ/wFeBs4JW0o/1e+poJvELyjIdz07ZFDQHmAv8MVF5RQbrNc4CzgPeSnAK5VdLQTJsZwAnAJ4HDgWHAjel9HtpcQ3J55wdJruA4GJhZQ04zM7NtXrt3ZoyI/5T0EPBVko65bdDiauDPwA8i4o5aNhgRfwD+ACDpquw8SQL+BbggIn6bTvsMSbFwInCppOHAqcDnIuLWtM1JwIvAUcBsSRNIioPDIuLutM0ZwJ2S9o2Ip2vJbGZmtq3q8DHTETEnIo4FhgJjgV2AoRFxbK1FQgF7pdvY+JCp9D4OfwYOTSdNBPpXtHkZmJdpM5nkio3s5Z13ASszbczMzKwDhZ/1kF7psLgbs0BSJAAsqpi+CBiXadPK5pdotrUZm2nzevYZFBERkhZn2mxG0unA6QBjxoyhsbGx5vDNzc2dWq7eypobypu9rLmhvNnLmhvKm72suaG82bsjd6eeHtkDKh8ypSrTKlW2qdY+dz0RcRlwGcCkSZNiypQphYJmNTY20pnl6q2suaG82cuaG8qbvay5obzZy5obypu9O3J3eOqhhy1M/63c6x/NpqMMC4EGYFQHbUanYx6AjeMfdmbLoxVmZmaWo7cVCvNJOvmj2yZIGkByZUPbeIMmkqdWZtvsBkzItLmH5OqKyZl1TwYGU/221GZmZlZFj596kDQE2Cd92w/YXdJBwJKIeEnSDODrkp4ieY7EN0gGJl4DEBHLJF0BXJiOOXgTuAh4DLgtbTNP0s0kV0mcRnLK4VLgRl/xYGZmVlw9xihMAuZk3n87ff0C+CzwA2Ag8BNgJHAf8IGIWJFZ5mxgPXBt2vZ24OSKW0t/iuSmUW1XR1wPnNnFn8XMzKxP6/FCISIa2fIJlNn5AUxPX3ltVpPckOmsdtosAT7dyZhmZmZG7xujYGZmZr2ICwUzMzPL5ULBzMzMcrlQMDMzs1wuFMzMzCyXCwUzMzPL5ULBzMzMcrlQMDMzs1wuFMzMzCyXCwUzMzPL5ULBzMzMcrlQMDMzs1wuFMzMzCyXCwUzMzPL5ULBzMzMcrlQMDMzs1y9rlCQNF1SVLwWZuYrbbNAUoukRkn7VaxjB0k/lvSGpJWSrpe0W89/GjMzs3LrdYVC6mlgl8zrgMy8rwDnAGcB7wUWA7dKGpppMwM4AfgkcDgwDLhRUkO3JzczM+tDtqt3gBzrI2Jh5URJAv4FuCAifptO+wxJsXAicKmk4cCpwOci4ta0zUnAi8BRwOwe+QRmZmZ9gCKi3hk2I2k6yVGDt4C1wH3AeRHxvKS9geeAv4mIBzLL3AS8ERGfkTQNuB0YHRGvZ9o8AfxPRHwrZ7unA6cDjBkzZuKsWbNqzt7c3MyQIUNqXq7eypobypu9rLmhvNnLmhvKm72suaG82Tube+rUqU0RManqzIjoVS/gg8DHgXeTHAFoBBYCOwGHAgHsXrHMlcDs9OcTgfWkRVCmzR3ApUUyTJw4MTpjzpw5nVqu3sqaO6K82cuaO6K82cuaO6K82cuaO6K82TubG3gwcvrEXnfqISL+mH0v6V7geeAzwL1tzSoWU5VplYq0MTMzs4zeOphxo4hoBp4AxpMcWQAYW9FsNLAo/Xkh0ACMaqeNmZmZFdDrCwVJA4B3Aq8B80kKgaMr5h8O3J1OagLWVbTZDZiQaWNmZmYF9LpTD5J+CNwAvERyFOCbwGDgFxERkmYAX5f0FPAM8A2gGbgGICKWSboCuFDSYuBN4CLgMeC2Hv44ZmZmpdbrCgVgN+BXJKcOXicZl3BIRLyYzv8BMBD4CTCS5KqID0TEisw6ziYZ0Hht2vZ24OSIaO2RT2BmZtZH9LpCISI+0cH8AKanr7w2q0luyHRWV2YzMzPb1vT6MQpmZmZWP73uiIKZmZltvdXrWpn9xEJefHMVe+w0iGP2G8uA/rU/ycCFgpmZdbuu6rR6WllzP/ryUk668j5aW4NVa1sZtH0D37xuLjNPeR8Hvm1ETetyoWBm25yy/vEva+6u7LR6Ullzr17XyklX3sfylvUbp61cm4zlP+nK+7j/vKNq+r1xoWBWZ2X941/W3GX941/W3F3dafWUta3R6dwRwYaA1g3BhghaNwStEcQGaE3fb5ye+Tn5d8vlNmxsR870dLkI5i5Yz603PsnqdRuqZmttDWY/sZAPHzSu8HfhQsE2U9Y//mXNXdY//mXNvS12Wnlisw4n04FlO6F0WqSdXuX0DVU6vo3LBzz6+nrumv00a3I6rTXrNvCt65/ggHHDM50lFduo3C5bdLJtHWW17JvaslmHvKGdz9C6IVj01ipWrKn+3S1vWc/E797Kdv36Vd3WhgAJGiT69RMNEg39RD+x8f1m0/tVa9s2TzSIjdPyprdNe/P19SzjLdaur/6dr1rbyktvrqrpd8WFQjdwp9WzyprbnVb7unKvbEMEf32rlYfmPJv7B3TNug1c8MenOHiPkZvtqWX31jZUbGvTtE3zN8Tm06vtAW7qzNi0zGbrZbO2Ly5qoXl19UfVrGhZz2H/cQcDt29IOr2qHfqmgqAtXwT0q+hoNnZWG6dt6sD6ZTq7zZbZYvlN85e+tZ5mFrMm7ztfv4EnXl222ba0xXo3ZeivfhvzbepY2Szfps/AltM262QzHXeVz/Cfv7+XO15aXzW3gJMO2YMzjnh7RR42bkNSp37vt1ZjYyPLRryd8/738Y3/X2YN2r6B3XcaVNM6XSh0MXdaPaunOq08teyVtc3fEMGrzRv4+V3zczuttes3cPHtz/C+vXbatLdUdU+JnL2ntk6OwntlG5fJdopV9sqeeaWFlTmdVvPq9Rz9oz8xbED/Lfawqm6rIl937pU1L19LC6/kHpJds34Ddz37Bq83r8n8sWeLziXbgVVO2367fpvl23LPsKIDq+iAt2ybzP/pHx5k0arqnRbAhw7clc8eumehzju7ve7uzJJOa3xupzV4+wZOO2Lvmg6D95TxIxq4b2HkdrYTdhnGyMHb1yFZx47ZbyzfvG5u1XkNDeKY/Sofl9Q+FwpdqN6dVjV5e2UbMn+831q9gVn3v8S69dX/+K9dv4HL/vw8h40fVdO5ssq9pw1RfK9ss72zdvbKHnu+hZU5hwdXrl7PP/zXX9hp8A5brHdjZ1m5rbaOsoO9srZpUNteWdv8lpbVrGV+bqe1et0GbnrsNea+unzjMmpbV7VtZfaiKqcV3SvbfE8rf6/s50vf4KUV1W9yGgGHj9+ZT753942dU+7eXNu6K76b7torK9JpnTltn17ZaTXu1MBDr+d3Wge9bQR77DS4Dsk61tWdVk+ZOKaBX/21+u95b84NMKB/AzNPed8WO60NDWLmKe+ruR9yodCFmha10tpavbNduaaVf7zsHsYOG1DzubL298oy66pxr6xtz2PdunW08gwt66r/T7F63Qb+v3tf5PanFm/sUKTKP/xbnivbYu+r6t5e9b2yLZbL2St78ZWFtG6onntDwAHjhvORg3fL7Gll9toq9r5qOaSa7QA7o0indc4H9u2VndZDj27H42+S22m9b68dOWC34XVI1jF3Wj2vqzutnrJ9mq9sudsc+LYR3H/eUcx+YiEvvbmK3X0fhd5h8aoNrKryxxOSDv5tIwfx9+/epWrHl92LqnrYsvJwZgfLZQ+ZdrRXVqTTOu+4Cb2y03r2mad5aml+p3XEO3bmb/epfOJ47+BOq+e506qPruy0elJZc7cZ0L+hS/5uu1DoQqMH9WPQ9g25ne3R7xrDsfvvUodkHXOn1fPcadVHWf/4lzV3m67qtHpaWXN3JRcKXcidVs9zp1UfZc3dpqx//Mua28rNhUIXcqdVH2XN3aasf/zLmtvMauNCoYu506qPsuY2M+vt+nyhIOmLwP8FdgGeAP4lIu7szm260zIzs76iX70DdCdJ/whcDHwfeA9wN/BHSbvXNZiZmVlJ9OlCAfgycFVEXB4R8yLiLOA14At1zmVmZlYKffbUg6TtgYnADytm3QIc2t6yTz/9NFOmTKl5m0uXLmXEiBE1L1dvZc0N5c1e1txQ3uxlzQ3lzV7W3FDe7N2Ru88WCsAooAFYVDF9EXBUZWNJpwOnA/Tv35+lS5fWvMHW1tZOLVdvZc0N5c1e1txQ3uxlzQ3lzV7W3FDe7N2Ruy8XCm0q76msKtOIiMuAywAmTZoUDz74YM0bamxs7NSRiHora24ob/ay5obyZi9rbihv9rLmhvJm72zu9u7g25fHKLwBtAKVdzkazZZHGczMzKyKPlsoRMRaoAk4umLW0SRXP5iZmVkH+vqph4uAmZLuB+4C/gnYFfjvuqYyMzMriT5dKETEtZJ2Ar5BcsOlucDfRcSL9U1mZmZWDn26UACIiJ8CP613DjMzszLqs2MUzMzMbOu5UDAzM7NcLhTMzMwslyK2uPfQNk/S60BnBjyOIrl/Q9mUNTeUN3tZc0N5s5c1N5Q3e1lzQ3mzdzb3HhGxc7UZLhS6kKQHI2JSvXPUqqy5obzZy5obypu9rLmhvNnLmhvKm707cvvUg5mZmeVyoWBmZma5XCh0rcvqHaCTypobypu9rLmhvNnLmhvKm72suaG82bs8t8comJmZWS4fUTAzM7NcLhTMzMwslwsFMzMzy+VCYStJOkLS9ZJelRSSPlvvTEVIOlfSA5KWS3pd0g2S9q93ro5I+pKkx9LcyyXdI+m4eueqlaTz0t+X/6p3lo5Imp5mzb4W1jtXUZJ2kfSL9Pd8taQnJR1Z71ztkfRCle88JN1U72wdkdQg6TuS5qff93xJ35XU6x9CKGmopBmSXpTUIuluSe+td65KHfU7SkyXtCD9HI2S9uvs9lwobL0hJI+v/megpc5ZajGF5KmahwLTgPXAbZJ2rGeoAl4BvgocDEwC7gB+L+nddU1VA0mHAKcBj9U7Sw2eJnlUe9vrgPrGKUbSCOAuQMBxwATgLGBxHWMV8V42/74PBgL4dT1DFfRV4EvA/wHeSfK38UvAufUMVdDPgGOAz5D8jt9C8ndxXF1TbamjfucrwDkkv+vvJfl9v1XS0M5szFc9dCFJzcCZEXFVvbPUStIQYBlwfETcUO88tZC0BDg3Ii6td5aOSBoOPERSKPwbMDcizqxvqvZJmg58NCJ6/RGnSpK+DxwZEX9b7yxbQ9LXgf8L7BoRq+qdpz2SbgTejIjPZKb9AtgpIv6+fsnaJ2kgsAI4ISKuy0xvAv4YEd+oW7h2VPY7kgQsAP4rIr6XThtIUiz8a2f+TvqIgrUZSvL78Fa9gxSVHuL8BEl1fXe98xR0GfA/EXFHvYPUaO/0MOd8SbMk7V3vQAUdD9wn6VpJiyU9IunM9I9pKaRZTwWu7u1FQuovwFRJ7wSQ9C6So5Z/qGuqjm0HNACrK6a3AIf1fJxO2wsYS3I0BICIaAH+THIEuWa9/pyR9ZiLgUeAe+qco0OSDiDJOQBoBv6fiHi8vqk6Juk0YB/gpHpnqdF9wGeBp4DRwDeAuyXtFxFv1jNYAXsDXwR+BFwAHAT8OJ3X68eHpI4m+eP/s3oHKeg/SHY8npTUStLPfC8iflrfWO2LiBWS7gG+IWkusBD4JDAZeLau4WozNv13UcX0RUCnTqG4UDAkXURSMR8WEa31zlPA0yR/8EcAJwC/kDQlIubWM1R7JO0LfB84PCLW1jtPLSLij9n3ku4Fnic5j3tRXUIV1w94MCLazo8/LGk8yTnzshQKpwEPRMQj9Q5S0D8CJwMnAk+Q/L96saT5EXFFPYMVcBJwJclYqFaS04S/IhkjUjaV4wpUZVohPvWwjZP0I5KqeVpEPF/vPEVExNqIeDYi2jqAR4Cz6xyrI5NJHv86V9J6SeuBI4Evpu93qG+84iKimaQDGF/vLAW8BjxZMW0esHsdstRM0mjgw8Dl9c5SgwuBH0bErIh4PCJmkhSUvX4wY0Q8FxFHkpzOfFtE/A3QH5hf32Q1absiaWzF9NFseZShEBcK2zBJF5NU/dMi4ql659kK/YDe3tH+nmQU9UGZ14PArPTn0hxlkDSAZDT7a/XOUsBdwL4V094BvFiHLJ3xOWANye9JWQwi2RvPaqVE/U1ErIyI1ySNJLkK4rqOlulF5pMUC0e3TUj/nz2cTo7l8qmHrZReLbBP+rYfsLukg4AlEfFS3YJ1QNJPSA6zHQ+8Jamt+mxO9xh7JUkXADcBL5OcBz2R5FLPXn0vhYhYCizNTpO0kuT3pNeeMgGQ9EPgBuAlkr2SbwKDgV/UM1dBPyIZT/F14FrgPSSX7Z1X11QFpIMYPw/MiogV9c5TgxuAr0maT3Lk6T3Al4Ff1jVVAZKOIfk7/hTJ3/ULSU51/ryeuSp11O9ImgF8XdJTwDMk44qagWs6tcGI8GsrXiSdVFR5XVXvbB3krpY5gOn1ztZB7qtI9gbXkFzucxtwTL1zdfKzNJJcwlT3LB3knEVyudVa4FXgt8C76p2rhvzHAY+SjGZ/hqRQUL1zFcg9Nf1/8m/qnaXG3EOBGen/py0k41m+Dwyod7YC2T8OPJf+fXmNZBzL8HrnqpKz3X6HZDzC9PQzrAb+BOzf2e35PgpmZmaWqzTnjMzMzKznuVAwMzOzXC4UzMzMLJcLBTMzM8vlQsHMzMxyuVAwMzOzXC4UzMzMLJcLBTMzM8vlQsGsBCRFgdcLXbSt4yV9uej0eultecz6Kt+Z0awEJB1SMel3JLclnp6ZtiYiHu6CbV0FHBURuxWZXi+9LY9ZX+WHQpmVQETcm30vaQ3wRuX0aiTtEBFrui3cVujN2bpD+qCn7SJiXb2zmBXlUw9mfYik6elpiP0lzZbUDPw6M/9YSfdIapG0TNLvJe2bmX8V8BlgXPaURt70zHIHSrpe0lvpuu+SdHgt2ap8lndI+p2kxZJWS3pJ0m8kbdfFeQ6QNEfSKkmvSfp3Sf2K5Cjw3+NRST+X9AVJ84B1wLEdLWfWm/iIglnfdB1wBfAfwAZIigSSR3TfAfwjMAT4d+Avkg6KiFeB7wA7A+8F/iFd1xpgec50JB0M3Ak8DJwGrAL+CbhN0qER0dRRthw3kjya+wvAG8A44O9IdnDycnYmz++BK4HzgWNIHqO9gU2nddrLkUvS9sCENOfb0/U2Aw+0t5xZb+NCwaxvuiQiLq6Y9l2SR/5+MCLWA0i6h+TRy+cAX46I5yS9Dqytcrqj6nTgQuAlYFpErE3bzgbmknSOxxfIthlJo4DxwIcj4vrMrGvSf3NzdiLP5RFxQfrzLZKGAedImkHyN7K9HO3ZH+hP8v2+PyJaCyxj1uv41INZ3/S77BtJg4GDgWvbigSAiJgP3AUc2ZmNSBqYLvsbYEN6WmA7QMBtwBEdZcvxJklRc4Gk0ySN78Y8lac/ZpEcbdm/szlS70n//VcXCVZmLhTM+qbXKt6PJOksK6cDLAR27OR2dgQaSPbU11W8zgRGZs/352TbQiSXYx0NPEhySuAZSc9L+kI35FmU837cVuSApFBYEBEPFmhr1mv51INZ31R53fNb6bSxVdqOJdlz7oylJOfzfwL8smqQiMpxCIWuyY6I54GT0ysFDiTp6H8q6YWI+GMX5hlDctQg+x7g1a3IAUmh4PEIVno+omC2DYiIlUAT8DFJDW3TJe0BHAr8KdN8DTCwymq2mJ6u906SDvShiHiw8tUF2SMiHgHabq60fxfn+XjF+0+QDDqcWzDHFtKjFu8m+c7NSs1HFMy2Hd8kuerhRkk/JTkP/21gGfCfmXZPAjumh9cfBFZHxOPtTP8y8GdgtqQrSE4tjCIZE9EQEV+rNaikdwMXA9cCz5KcTvgssJ7kqo32ctaa57S0Y3+A5KqHzwPTI2JpwRzVjCf5fh+q9bOb9TYuFMy2ERFxs6TjgG+RDOBbCzQCX4mIBZmmPwMOAb4PjABeBPbMmx4RD0l6b7reS4DhwOskneR/dzLuQpIrF74M7AasBh4H/j5zeWNX5fkw8GOSQmoZydUh36khRzVtAxldKFjp+RbOZrZNkjSdpJjon70SxMw25zEKZmZmlsuFgpmZmeXyqQczMzPL5SMKZmZmlsuFgpmZmeVyoWBmZma5XCiYmZlZLhcKZmZmlsuFgpmZmeVyoWBmZma5/n8Npu1tmTigDQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dif = (ks[1] - ks[0]) / 10\n",
    "l = len(ks)\n",
    "\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.errorbar(np.array(ks) + dif, uni_difs_k, yerr=uni_stds_k, color='tab:blue', label=r'$F_Q$', fmt='o-', capsize=7, linewidth=1, markersize=7)\n",
    "plt.errorbar(np.array(ks) - dif, [var_dif]*l, yerr=[var_std]*l, color='tab:red', label=r'$F_H$', fmt='s-', capsize=7, linewidth=1, markersize=7)\n",
    "plt.hlines(y=0, xmin=ks[0]-0.3, xmax=ks[-1]+0.3, colors='black', linestyles='solid')\n",
    "plt.legend(fontsize=14)\n",
    "plt.grid()\n",
    "plt.xlim(ks[0]-0.3, ks[-1]+0.3)\n",
    "plt.xticks(ks, fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "plt.xlabel(r'Trotter steps $r$', fontsize=16)\n",
    "plt.ylabel(r'$\\delta$', fontsize=16)\n",
    "# plt.savefig(str(n_qubits) + 'q-' + str(n_instances) + 'n' + '-delta_k.png', bbox_inches='tight')\n",
    "# plt.savefig(str(n_qubits) + 'q-' + str(n_instances) + 'n' + '-delta_k.pdf', bbox_inches='tight')\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot(np.array(ks) + dif, uni_gate_counts, color='tab:blue', label=r'$F_Q$', marker='o', linewidth=1, markersize=7)\n",
    "plt.plot(np.array(ks) - dif, [var_gate_count]*l, color='tab:red', label=r'$F_H$', marker='s', linewidth=1, markersize=7)\n",
    "plt.hlines(y=0, xmin=ks[0]-0.3, xmax=ks[-1]+0.3, colors='black', linestyles='solid')\n",
    "plt.legend(fontsize=14)\n",
    "plt.grid()\n",
    "plt.xlim(ks[0]-0.3, ks[-1]+0.3)\n",
    "plt.xticks(ks, fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "plt.xlabel(r'Trotter steps $r$', fontsize=16)\n",
    "plt.ylabel(r'Gate count', fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Number of terms in the Ising Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N^2_terms expected: 36\n",
      "N^2_terms actual: 36\n",
      "\n",
      "N_x expected: 35\n",
      "N_x actual: 35\n"
     ]
    }
   ],
   "source": [
    "n_qubits = 5\n",
    "ham = hamiltonian_dict_to_matrix(ising_model_dict(n_qubits, 1, 1))\n",
    "\n",
    "print('N^2_terms expected:', 1 - 3*n_qubits + 2*n_qubits**2)\n",
    "print('N^2_terms actual:', len(hamiltonian_matrix_to_dict(n_qubits, ham.dot(ham))))\n",
    "print()\n",
    "\n",
    "n_X = 0\n",
    "for key in hamiltonian_matrix_to_dict(n_qubits, ham.dot(ham)).keys():\n",
    "    for char in key:\n",
    "        if char == 'X':\n",
    "            n_X += 1\n",
    "print('N_x expected:', n_qubits * (2*n_qubits - 3))\n",
    "print('N_x actual:', n_X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}