{

 "cells": [

  {

   "cell_type": "code",

   "execution_count": 43,

   "metadata": {

    "id": "in3VVDRqlMUE",

    "tags": []

   },

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "Time in seconds since beginning of run: 1693283058.6479008\n",

      "Tue Aug 29 04:24:18 2023\n"

     ]

    }

   ],

   "source": [

    "# This cell is added by sphinx-gallery\n",

    "# It can be customized to whatever you like\n",

    "%matplotlib inline\n",

    "# !pip install pennylane\n",

    "# !pip install neural_tangents\n",

    "# !pip install networkx\n",

    "import time\n",

    "seconds = time.time()\n",

    "print(\"Time in seconds since beginning of run:\", seconds)\n",

    "local_time = time.ctime(seconds)\n",

    "print(local_time)"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "cfjAkbwNlMUF"

   },

   "source": [

    "Machine learning for quantum many-body problems\n",

    "===============================================\n",

    "\n",

    "::: {.meta}\n",

    ":property=\\\"og:description\\\": Machine learning for many-body problems\n",

    ":property=\\\"og:image\\\":\n",

    "<https://pennylane.ai/qml/_images/ml_classical_shadow.png>\n",

    ":::\n",

    "\n",

    "::: {.related}\n",

    "tutorial\\_classical\\_shadows Classical Shadows\n",

    "tutorial\\_kernel\\_based\\_training Kernel-based training with\n",

    "scikit-learn tutorial\\_kernels\\_module Training and evaluating quantum\n",

    "kernels\n",

    ":::\n",

    "\n",

    "*Author: Utkarsh Azad --- Posted: 02 May 2022. Last Updated: 09 May\n",

    "2022*\n",

    "\n",

    "Storing and processing a complete description of an $n$-qubit quantum\n",

    "mechanical system is challenging because the amount of memory required\n",

    "generally scales exponentially with the number of qubits. The quantum\n",

    "community has recently addressed this challenge by using the\n",

    "`classical shadow <tutorial_classical_shadows>`{.interpreted-text\n",

    "role=\"doc\"} formalism, which allows us to build more concise classical\n",

    "descriptions of quantum states using randomized single-qubit\n",

    "measurements. It was argued in Ref. that combining classical shadows\n",

    "with classical machine learning enables using learning models that\n",

    "efficiently predict properties of the quantum systems, such as the\n",

    "expectation value of a Hamiltonian, correlation functions, and\n",

    "entanglement entropies.\n",

    "\n",

    "![Combining machine learning and classical\n",

    "shadows](/demonstrations/ml_classical_shadows/class_shadow_ml.png){.align-center\n",

    "width=\"80.0%\"}\n",

    "\n",

    "In this demo, we describe one of the ideas presented in Ref. for using\n",

    "classical shadow formalism and machine learning to predict the\n",

    "ground-state properties of the 2D antiferromagnetic Heisenberg model. We\n",

    "begin by learning how to build the Heisenberg model, calculate its\n",

    "ground-state properties, and compute its classical shadow. Finally, we\n",

    "demonstrate how to use\n",

    "`kernel-based learning models <tutorial_kernels_module>`{.interpreted-text\n",

    "role=\"doc\"} to predict ground-state properties from the learned\n",

    "classical shadows. So let\\'s get started!\n",

    "\n",

    "Building the 2D Heisenberg Model\n",

    "--------------------------------\n",

    "\n",

    "We define a two-dimensional antiferromagnetic [Heisenberg\n",

    "model](https://en.wikipedia.org/wiki/Quantum_Heisenberg_model) as a\n",

    "square lattice, where a spin-1/2 particle occupies each site. The\n",

    "antiferromagnetic nature and the overall physics of this model depend on\n",

    "the couplings $J_{ij}$ present between the spins, as reflected in the\n",

    "Hamiltonian associated with the model:\n",

    "\n",

    "$$H = \\sum_{i < j} J_{ij}(X_i X_j + Y_i Y_j + Z_i Z_j) .$$\n",

    "\n",

    "Here, we consider the family of Hamiltonians where all the couplings\n",

    "$J_{ij}$ are sampled uniformly from \\[0, 2\\]. We build a coupling matrix\n",

    "$J$ by providing the number of rows $N_r$ and columns $N_c$ present in\n",

    "the square lattice. The dimensions of this matrix are $N_s \\times N_s$,\n",

    "where $N_s = N_r \\times N_c$ is the total number of spin particles\n",

    "present in the model.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 44,

   "metadata": {

    "id": "LWfXJdIBlMUG"

   },

   "outputs": [],

   "source": [

    "import itertools as it\n",

    "import pennylane.numpy as np\n",

    "import numpy as anp\n",

    "\n",

    "def build_coupling_mats(num_mats, num_rows, num_cols):\n",

    "    num_spins = num_rows * num_cols\n",

    "    coupling_mats = np.zeros((num_mats, num_spins, num_spins))\n",

    "    coup_terms = anp.random.RandomState(24).uniform(0, 2,\n",

    "                        size=(num_mats, 2 * num_rows * num_cols - num_rows - num_cols))\n",

    "    # populate edges to build the grid lattice\n",

    "    edges = [(si, sj) for (si, sj) in it.combinations(range(num_spins), 2)\n",

    "                        if sj % num_cols and sj - si == 1 or sj - si == num_cols]\n",

    "    for itr in range(num_mats):\n",

    "        for ((i, j), term) in zip(edges, coup_terms[itr]):\n",

    "            coupling_mats[itr][i][j] = coupling_mats[itr][j][i] = term\n",

    "    return coupling_mats"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "ZV-J9mFJlMUG"

   },

   "source": [

    "For this demo, we study a model with four spins arranged on the nodes of\n",

    "a square lattice. We require four qubits for simulating this model; one\n",

    "qubit for each spin. We start by building a coupling matrix `J_mat`\n",

    "using our previously defined function.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 45,

   "metadata": {

    "id": "GSqBR9p2lMUG"

   },

   "outputs": [],

   "source": [

    "Nr, Nc = 2, 2\n",

    "num_qubits = Nr * Nc  # Ns\n",

    "J_mat = build_coupling_mats(1, Nr, Nc)[0]"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "LMJScKXHlMUG"

   },

   "source": [

    "We can now visualize the model instance by representing the coupling\n",

    "matrix as a `networkx` graph:\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 46,

   "metadata": {

    "colab": {

     "base_uri": "https://localhost:8080/",

     "height": 254

    },

    "id": "PnCxl2OUlMUH",

    "outputId": "8383df52-8539-489f-edfc-3064c3cf6eb0"

   },

   "outputs": [

    {

     "data": {

      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAGjCAYAAACBlXr0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAiVElEQVR4nO3deXxU9b3/8fdkmYQkEEIgCcgmKJsLCFYRtFbFRFup9KfWixT4qbUqVq/VWi7eytZSRavVKqIX3C+KdW2pWyySgooiILggBMMSARMgBLIM2WbO/YNFQxKSGZL5nElez8ejj0dzzplzvkzT7yvnzJkZj+M4jgAAMBRlPQAAAIgRAMAcMQIAmCNGAABzxAgAYI4YAQDMESMAgDliBAAwR4wAAOaIEQDAHDECAJgjRgAAc8QIAGCOGAEAzBEjAIA5YgQAMEeMAADmiBEAwBwxAgCYI0YAAHPECABgjhgBAMwRIwCAOWIEADBHjAAA5mKsB9BUjuOowh9QlT+ggCNFeSRvdJTaxURbDw0AXCGS50lXx6i8qkb5Jfu1p6JaxRVVqvI7dbbxRnuUEu9Vp/hY9ezQToleV/+TAKBZtZZ50uM4Tt2RG3IcR4Xllcor9qnQVymPpKYM8NB26Qlx6puSoPTEOHk8npYdLAAYaI3zpKti5Kv2a3XBXu30VTX5yT3SocelJXg1NKOjEmLdf3oKAE3VWudJ18Ro6z6f1hSWKOA4IT25R/JIivJ4NCS9g3olJzTDHgHAVmueJ81j5DiOvioq0/qishY7xoDUJA1MTXLN6SgABKMtzJPmt3a39BMsSeuLyvRVCx8DAFpKW5gnTWO0dZ+vxZ/gQ9YXlWnrPl9YjgUAzaWtzJNmMfJV+7WmsCSsx1xTWCJftT+sxwSAULWledIkRo7jaHXBXgXC/HJV4OBxXXLPBgA0qK3NkyYxKiyv1E5fVbPcDRIMR9JOX5UKyyvDfGQACE5bmydNYpRX7FOw92tUV1XquT//Ub885zSNHdxH//Xzn2jtB/8O+tieg8cHADcLZZ7cX16uhX+9T3/45VWaeOYgXTagm9579cWgj20xT4Y9RuVVNSr0VQZd+4f/61Ytevp/dM7on+nqO2cqKipKs64fr69WfRzUfhxJhb5KlVfVNLptaV6e8l9+WfvWrQtytADwnX3r1in/5ZdVmpfXpO1DnSdLi/fopUf/om2bNqpX/0HBD/SgYObJ5hL2GOWX7A+69hs/+1QfvPl3jfvNFE383VRlXvkLTX/mJXXp1l3P3ffHoMfgOTiOhuxbt06f3HijcjIztXbyZC0dPVr5L70U9HEQeZYuXarRo0erW7du8ng8ev311xt9zJw5czRw4EC1a9dO/fv317PPPltr/bx583TOOecoJSVFKSkpGjVqlFasWNFC/wK4Tf7f/qalo0dr7eTJysnM1Cc33tjoH7ihzJOSlJKWpvnL1ujx9z7RhDvuCm3ABzU2Tza3sMdoT0V10LVf/s4/FRUdrQuv/MXhZd64eF1w2VhtWLNKu7/dHtT+nIPjONKhCC0dPVoF2dm11n39+ONBjhqRqLy8XIMHD9acOXOatP3cuXM1ZcoUTZ8+XV9++aVmzJihm266SYsWLTq8TU5OjsaOHaslS5Zo+fLl6tGjhzIzM7V9e3C/t4hMR84dBdnZWjp69FGjFMo8KUmx3jildEkL4ZF1NTRPtpSwfnSr4zgqrqgK+nGbv/pC3Xr3UUJS+1rLTzh1yMH1X6pz1+OC2ueu3Xu0cvY0SVJVcbH2fvaZ/OXlDW5fvnmz3h05MriBw5X6XHut+l5zTb3rLr74Yl188cVN3tdzzz2n66+/XldeeeWBfffpo08++USzZ8/W6NGjJUkLFiyo9Zj58+frlVde0eLFizVhwoR695v35JPa9MQTTR4H3KuioKDe5QXZ2SrIzlZMYqKSTz1V3pQUSQcisOuXv5Hi24VxlPUrbq0xOvA9G8H3vnjXTqV0Sa+z/NCy4p2FQe8z0C5RO5Z/LE9xUZMf09AvFSJLTVnzvYGwsrJS8fHxtZa1a9dOK1asUHV1tWJjY+s8xufzqbq6Wp06dTrqGPl9axtqystVtHz54Z+dlFQFfm0fIkmq8ge0v8Yflu9DCutluip/ILTHVVQoxuutszw2Lu7A+sqK0AbUoUNojwMOysrK0vz587Vq1So5jqOVK1dq/vz5qq6u1u7du+t9zOTJk9WtWzeNGjUqzKNFREhOth5BLaHO28EK65lRIMQb5r3x8aqpqnt5r7rywH3w3rj4OuuaJLZu4IBg3HXXXSooKNDw4cPlOI7S09M1ceJE3XvvvYqKqvu33j333KOFCxcqJyenzhkVIEmKqXs2bSnUeTtYYT0zigrxw2BTuqSpeFfdS3GHlqWk1b2E1yTVwb9+BXxfu3bt9OSTT8rn82nLli3Kz89X79691b59e3Xp0qXWtn/+8591zz33KDs7W6eeeqrRiOF6NeF7naYpQp23gxXWMyNvdGjt6z3gJH3x8YfylZXWuolh49pPJUnHDzwppP2eu/AFxR98pn3btilv3jwVvPuuHH/dz2WKSUrSSG7vbhW8qanNvs/Y2Fh1795dkrRw4UJdcskltc6M7r33Xs2aNUvvvPOOTj/99Eb312vcOHW96KJmHyfC7/0rrpC/vtcpo6KUkZmpE667TgkHf3ckqSLgaHFx+N7f05hQ5+1ghTVG8dFR8kZ7gr6J4aysS/SPJx/Tuy/+ry699kZJBz6R4b3XXtSJg4cGfSeddOAJTk777i/XuM6ddfqcOfJt26aNjz6qb155RU7Nd78Qaeeeqw79+gV9HESWsrIyff3114d/3rx5s9asWaNOnTqpZ8+emjJlirZv3374vUS5ublasWKFzjzzTBUXF+uBBx7QF198oWeeeebwPmbPnq2pU6fq+eefV+/evVVw8MaEpKQkJSUl1TuO+NRUxbdANBF+6eeeqx1vvHH4Z09MjHpcdplOnDSpVoQO8TqOvCWFId3s1dy80VFhuXlBCnOMPB6PUuK9QX/mUb/BQ3XWRaO14C93a9+e3croebxyXv+bdm3/RpP+eH9IY0mJr/+6bEL37hr8pz/pxEmTtOmJJ7Rn9Woln3SSBk2ZEtJxEFlWrlyp88477/DPt912myRp4sSJevrpp/Xtt98qPz//8Hq/36/7779fGzZsUGxsrM477zx9+OGH6t279+Ft5s6dq6qqKl1++eW1jjVt2jRNnz69Rf89sHfqrFmKSUrSvi+/VKehQ9Xn2mvrjdAhoc6Th7z5v0/KV1qiPQfvMl655F3tKfxWknTxL65RYvum37jV0DzZEsL+Ta9f7S7V+qKyoN/QVVVZoRceuldLF72q8n371Kv/QP3HLb/Taef8KOgxeHTwWw07t290WwAIt1DnSUm64fwztGvHtnrXzf3Xx0rr3qNJ+wn3PBn2GJVX1eidzbvCech6ZR3fRYnesJ4YAkCTtMV5MuwfB5TojVF6QlxIn7vUHDyS0hPiCBEA12qL86TJV0j0TUkI+3d0HOIcPD4AuFlbmydNYpSeGKe0BG/Yq++RlJbgVXpiXJiPDADBaWvzpEmMPB6PhmZ0VJQnvE9z1MHjesJ8XAAIVlubJ01iJEkJsdEakh7ez4Ybkt5BCbHhuWceAI5VW5onzWIkSb2SEzQgtf43/TW3AalJ6pXMa0UAIktbmSdNYyRJA1OTWvyJHpCapIFh+h8TAJpbW5gnw/4+o4Zs3efTmsISBRynWe4g8ejAtc8h6R04IwLQKrTmedI1MZIkX7Vfqwv2aqevSh4ppCf70OPSErwamtGR14gAtCqtdZ50VYykA19NXlheqbxinwp9lU1+sg9tl54Qp74pCUpPjOOuOQCtUmucJ10Xo+8rr6pRfsl+7amoVlFJuWqi674b2BsdpZT4WHWKj1XPDu34ZAUAbcr358k95RWqruedSZEwT7o6Rt9XsHixVkyecuCrwmO9SjnlFJ1+z5/C9vHmAOB25Vu2aPHlPz88T6q6SqkDB2rknIeth9Yo9+XxKDzFRVJxkSQpuksnQgQAR/j+PClJ0T0b/roKNzG/tbupSjZsqPXzvs8/NxoJAESO6vq+ZdaFIiZG2157rdbPNWVl8m3fbjQaAHCf8m11v8eo5KuvDEYSvIiI0b5161S+aVOd5RsffdRgNADgTpufeqrOskBFhfatW2cwmuBERIxyH67/xbdvXn6ZsyMAkOTbtk27li2rd11Dc6ibuD5G+9atU0F2dr3rnJoazo4AQAeuFDl+f73rCrKzXX925PoYNVZ0zo4AtHW+bdv0zSuvHHUbt58duTpGjuNo1wcfHH2bmhoVffxxmEYEAO6TN2+enJqao25TkJ2t0tzcMI0oeK6OkRxH7bp1a3Sz+PT0MAwGANxp14cfNmm73R991MIjCZ2rY+SJitKwhx7ScWPGKC4trfa66GilDh+uU2fNUpeRI41GCAD2Env0aNJ2CT17tvBIQhcxHwe0/oEHtHHOnMM/x3bsqItWrTIcEQC4Q/mWLVrxq1+pLC+v3vWemBj1vuoqnXTXXfJEufMcJKI+DggAUFdi79760TvvaPs//6lPb7211rqYDh104bJlikly9xeMujORAICgeDwexXXqVHd5dLTrQyQRIwCACxAjAIA5YgQAMEeMAADmiBEAwBwxAgCYI0YAAHPECABgjhgBAMwRIwCAOWIEADBHjAAA5ogRAMAcMQIAmCNGAABzxAgAYI4YAQDMESMAgDliBAAwR4wAAOaIEQDAHDECAJgjRgAAc8QIAGCOGAEAzBEjAIA5YgQAMEeMAADmiBEAwBwxAgCYI0YAAHPECABgjhgBAMwRIwCAOWIEADBHjAAA5ogRAMAcMQIAmCNGAABzxAgAYI4YAQDMESMAgDliBAAwR4wAAOaIEQDAHDECAJgjRgAAc8QIAGCOGAEAzBEjAIA5YgQAMEeMAADmiBEAwBwxAgCYI0YAAHPECABgjhgBAMwRIwCAOWIEADBHjAAA5ogRAMAcMQIAmCNGAABzxAgAYI4YAQDMESMAgDliBAAwR4wAAOaIEQDAHDECAJgjRgAAc8QIAGCOGAEAzBEjAIA5YgQAMEeMAADmiBEAwBwxAgCYI0YAAHPECABgjhgBAMwRIwCAOWIEADBHjAAA5ogRAMAcMQIAmCNGAABzxAgAYI4YAQDMESMAgDliBAAwR4wAAOaIEQDAHDECAJgjRgAAc8QIAGCOGAEAzBEjAIA5YgQAMEeMAADmiBEAwBwxAgCYI0YAAHPECABgjhgBAMwRIwCAOWIEADBHjAAA5ogRAMAcMQIAmCNGAABzxAgAYI4YAQDMESMAgDliBAAwR4wAAOaIEQDAnMdxHMd6EE3hBAJyAoHvFng8ioqOthsQALiM4/crUF1de6HHo+i4OJsBBSFiYgQAaL24TAcAMEeMAADmiBEAwBwxAgCYI0YAAHPECABgjhgBAMxFZIxycnK0f/9+62EAgCusW7dOkyZN0mmnnaauXbuqa9euOu200zRp0iStW7fOenhNEpFvevV6vVq7dq0GDhxoPRQAMPXWW29pzJgxGjp0qLKyspSeni5JKiws1LvvvqtVq1bp73//u7KysoxHenSujtHQoUPrXb5mzRoNGDBA8fHxkqTVq1eHc1gA4BqDBw/WpZdeqpkzZ9a7fvr06Xr11Vf12WefhXlkwYmxHsDRfP755xo1apSGDx9+eJnjOFq7dq3OO+88paWlGY4OAOzl5uZq3LhxDa4fO3asZs+eHcYRhcbVMcrJydHEiRN1xhlnaNq0aYqKOvAS16xZs3TTTTdp0KBBxiMEAFu9e/fWG2+8of79+9e7/o033lCvXr3CPKrguTpGI0eO1KpVq3TDDTdoxIgRWrBggfr27Ws9LABwjZkzZ+qqq65STk6ORo0aVes1o8WLF+vtt9/W888/bzzKxrk6RpKUnJysF154QU899ZTOPvtszZgxQx6Px3pYAOAKV1xxhY477jj99a9/1f3336+CggJJUkZGhs466yzl5OTorLPOMh5l41x9A8ORNm7cqHHjxmnlypX64osvuEwHAK1ERMVIkgKBgEpLS9WhQwfOkACglYi4N71GRUUpOTmZEAFAE9x555265pprrIfRqIiL0fdNnDhR559/vvUwAMC1tm/fri1btlgPo1Guv4HhaLp163b4dm8AQF3PPPOM9RCaJOJeMwIAtD4RfVrxzTffRMS1UABoSfv379f7779f74eiVlRU6NlnnzUYVXAi+sxo7dq1Gjp0qPx+v/VQAMBEbm6uMjMzlZ+fL4/Ho7PPPlsLFy5U165dJR1482u3bt1cP0+6+jWjf/zjH0ddv2nTpjCNBADcafLkyTr55JO1cuVK7d27V7feeqtGjhypnJwc9ezZ03p4TebqM6OoqCh5PB4dbYgej8f1xQeAlpKenq5//etfOuWUUyQd+DDpSZMm6c0339SSJUuUmJgYEWdGrn7NqGvXrnr11VcVCATq/Q9fHQGgrdu/f79iYr67yOXxeDR37lyNHj1a5557rnJzcw1H13SujtGwYcO0atWqBtc3dtYEAK3dgAEDtHLlyjrLH3nkEV166aX66U9/ajCq4Lk6RnfccYdGjBjR4PoTTjhBS5YsCeOIAMBdfvazn+mFF16od90jjzyisWPHRsQf7a5+zQgA0Da4+swIANA2uPrW7u9zAgE5gcB3CzweRUVH2w0IAFzG8fsVqK6uvdDjUXRcnM2AghAxMdrw4IPaOGfO4Z9jO3bURUe5uQEA2prdH32kjyZMqLUsNiVFF9Vzg4PbcJkOAGCOGAEAzBEjAIA5YgQAMEeMAADmiBEAwBwxAgCYI0YAAHPECABgjhgBAMwRIwCAOWIEADBHjAAA5ogRAMAcMQIAmCNGAABzxAgAYI4YAQDMESMAgDliBAAwR4wAAOaIEQDAHDECAJgjRgAAc8QIAGCOGAEAzBEjAIA5YgQAMEeMAADmiBEAwBwxAgCYI0YAAHPECABgjhgBAMwRIwCAOWIEADBHjAAA5ogRAMAcMQIAmCNGAABzxAgAYI4YAQDMESMAgDliBAAwR4wAAOaIEQDAHDECAJgjRgAAc8QIAGCOGAEAzBEjAIA5YgQAMEeMAADmiBEAwBwxAgCYI0YAAHPECABgjhgBAMwRIwCAOWIEADBHjAAA5ogRAMAcMQIAmCNGAABzxAgAYI4YAQDMESMAgDliBAAwR4wAAOaIEQDAHDECAJgjRgAAc8QIAGCOGAEAzBEjAIA5YgQAMEeMAADmiBEAwBwxAgCYI0YAAHPECABgjhgBAMwRIwCAOWIEADBHjAAA5ogRAMAcMQIAmCNGAABzxAgAYI4YAQDMESMAgDliBAAwR4wAAOaIEQDAHDECAJgjRgAAc8QIAGCOGAEAzBEjAIA5YgQAMEeMAADmiBEAwBwxAgCYI0YAAHPECABgjhgBAMwRIwCAOWIEADBHjAAA5ogRAMAcMQIAmCNGAABzxAgAYI4YAQDMESMAgDliBAAwR4wAAOaIEQC0Ao7jqHLPnror/H7VlJWFf0BBIkYAEOHKNm9WTlaWPr311jrrqktK9PawYfp8+nQ5gUD4B9dEro9RyYYNWn377cp/6aVay2tKS/XhVVdp68KFRiMDAHf48g9/UFleXoPrnZoabXnuOe1cujSMowpOjPUAjsYJBLTqlltU9vXXddf5/Sr6+GMVffyxEnr0UJeRIw1GCAD2yr/5pknb+fLzW3gkoXP3mZHHo/3fftvoZhWFhWEYDAC4U5cRI5q0Xefhw1t4JKFzdYw8Hk+jZzyemBilnnlmmEYEAO7T97rr5Ik5+oWujMxMte/XL0wjCp6rYyRJ/W6++ajre1x+uRKOOy5MowEA90no3l09LrvsqNs0Npdac32MkgcNUkZmZr3rPDExOnHSpDCPCADc58RJk6So+qf0jMxMJQ8aFOYRBcf1MZIaLjpnRQBwQEL37kr74Q/rXef2syIpQmKUPGiQEvv0qbOcsyIA+M7xV19dZ1l0fLzrz4qkCImRJHUfM6bWzzFJSZwVAcD3JHbvXmdZhwgIkRRBMeowYECtn5NPOcVoJAAQOWISE62H0CQRE6MjRTVyGyMAIHJEzIzuSHJSUqXkZCkmVtVd0rW/xq92MdHWQwMAV3Cc2vOkaqrlT2xvPawm8TiO41gPoiHlVTXKL9mvPRXV2lO2X9Weuidy3miPUuK96hQfq54d2inRGzF9BYBjVmue9FWqup4ZPRLmSdfFyHEcFZZXKq/Yp0JfpTw6cFbUmEPbpSfEqW9KgtIT4+TxeFp2sABgoDXOk66Kka/ar9UFe7XTV9XkJ/dIhx6XluDV0IyOSojlMh6A1qO1zpOuidHWfT6tKSxRwHFCenKP5JEU5fFoSHoH9UpOaIY9AoCt1jxPmsfIcRx9VVSm9UUt902EA1KTNDA1yTWnowAQjLYwT5rf2t3ST7AkrS8q01ctfAwAaCltYZ40jdHWfb4Wf4IPWV9Upq37fGE5FgA0l7YyT5rFyFft15rCkrAec01hiXzV/rAeEwBC1ZbmSZMYOY6j1QV7FQjzy1WBg8d1yT0bANCgtjZPmsSosLxSO31VzXI3SDAcSTt9VSosrwzzkQEgOG1tnjR5G25esS+o++O//nyNlrz2N32x4kPt2v6N2ndM0YmDh+mq//yduh3fN6hjew4ePyMpPthhA0DYBDtP5m/coL89cr/yvvxMe3fvVFx8O3U/oZ8uveZG/eD8+r+gtCEW82TYY1ReVaNCX3DFfW3eHK3/9BONyLpEvfoP1N7du/TWgqd0x2VZunvhP9Wz34DGd3KQI6nQV6nyqhpXfiQGAIQyT+7asU37y8t03pgrlJKWocqK/foo+w3dM+n/6/oZ9yrzyl80eV8W82TY32f01e5SrS8qC+rUc/3qT9T35MGK9XoPL9uxZZNu++kFOivrJ/rP+x4JagweHbynvnNkfIAgwm/OnDm67777VFBQoMGDB+vhhx/WGWec0ejjFi5cqLFjx+rSSy/V66+/Xu82N9xwgx5//HH95S9/0a233tq8A0erEMo8WR+/36/fXZalqspKPfzWsqAeG+55MuyvGe2pqA76CR4w9Ae1QiRJ3Xr3UY8T+mlb3sagx+AcHEd9fNu26bOpU7V8wgSV5uYGvW9EvhdffFG33Xabpk2bptWrV2vw4MHKysrSzp07j/q4LVu26Le//a3OOeecBrd57bXX9NFHH6lbt27NPWxEgNLcXC2fMEGfTZ0q37ZtDW4XyjxZn+joaKVmdJOvNPg78o42T7aEsMbIcRwVV1Q12772Fu1W+5ROIT2++Ign2bdtm9beeafeu+ACbV2wQLs/+ECbn3uuOYaKCPPAAw/ouuuu09VXX61BgwbpscceU0JCgp588skGH+P3+zVu3DjNmDFDffr0qXeb7du36+abb9aCBQsUGxvbUsOHi21+9lnt/uADbV2wQO9dcIHW3nlnnSgd6zxZ4fOppLhIBflbtOjp/9Gny5bolOFnh7SvI+fJlhTWF00q/AFV+ZvnquDSRa9qT+G3+o9bfhvS46v8AZXu3Sf/tm+UN3++vn3zTTn+2vfWV+7eLd+OHc0xXLhIVFyc4lNT611XVVWlVatWacqUKd9tHxWlUaNGafny5Q3uc+bMmUpLS9O1116rZcvqXg4JBAIaP3687rjjDp100kmNjrGiqEiBSu76bG0qi4oO/3enpkb5L76o/JdfVrcf/1h9r71Wib17H/M8+czsGcp+8cAf0lFRUTrzwh/rl3fNCmlfVf5A2L43LqwxqvIHmmU/2zZt1PyZd6r/kGH60Zifh7yfJf/vMnm2bm5wfUF2tgqys0PeP9wp48IL9YPHHqt33e7du+X3+5Wenl5reXp6utavX1/vY95//3098cQTWrNmTYPHnD17tmJiYnTLLbc0aYyf//d/q+Ddd5u0LSKc368dixZpx6JFkiSndx/p0WdD3t1PJv5Sw7N+ouKdhfrwrUUKBPyqqQ79DKfKHwhLjMJ6mS7QDCdFxbt26k/XT1BC+/b67UPzFB19DE9SrLfxbYCjKC0t1fjx4zVv3jx17ty53m1WrVqlhx56SE8//TQf1ovGxRzbJdzufU7U4BE/1I/GXKE7H39WFeXluvvGiSG/ibU55u2mCOuZUdQx/v+wvLREs341TuUlJfrjgtfUKT3jmPbXadgwlXy7Tf7y8nrXRyckyNsptNek4F6xDVyik6TOnTsrOjpahYWFtZYXFhYqI6Pu71teXp62bNmi0aNHH14WCBy4AhATE6MNGzZo2bJl2rlzp3r27Hl4G7/fr9tvv10PPvigtmzZUu8Y23XvHuw/DS5XtWeP/L76P/stOjFRnUeMUKB7LzXniwPDsy7R49N+px2b83RcnxOCfvyxzttNFdYYeaNDPxGrqqzQ3TdO1I4tmzTtyRfV44R+xzyeYTOnq92f/qB969Yp9+GH61ySy8jM1ND77z/m4yByeL1eDRs2TIsXL9aYMWMkHYjL4sWL9etf/7rO9gMGDNDnn39ea9nvf/97lZaW6qGHHlKPHj00fvx4jRo1qtY2WVlZGj9+vK6++up6xzFkVmjX+OFuq2+/XduPuOU/IzNT/W6+WcmDBkmS9tf4tSPv6HduBqOqskKS5CsrDenxxzJvByOsMYqPjpI32hP0i3N+v18P/OYG5a5ZpclznlL/004/5rF4o6MOXwdNHjRIP5g7t06Ukhq4Kwqt22233aaJEyfq9NNP1xlnnKEHH3xQ5eXlh8MxYcIEHXfccbr77rsVHx+vk08+udbjO3bsKEmHl6empir1iLOx2NhYZWRkqH///i3/D4JrfH9OOTJCh4Q6T+4r2q3k1NqXimuqq/Xv11+SNz5e3fsG/wf89+fJlhbWGHk8HqXEe4P+zKNnZs/QJ+9l6/TzLlTZvr369z9eqbX+3J9eFvRYUuLrXpc9FKXS3FyVb92q9PPPD3q/iHxXXnmldu3apalTp6qgoEBDhgzR22+/ffimhvz8fEVFmX8VGCLQCddfrw79+imxVy+171d/HEKdJx+b9jvtLyvToNPPVKf0DO3dvUtLF72q7Zu+1sTJ09QuMTHo8dY3T7aUiPgEhqnjL9OXnzR8W+0r64O7wsonMABws1DmyfffeF2LX3lB+bnrVbq3WO0Sk9TnpFP0419cox+cnxX0GMI9T4Y9RuVVNXpn865wHrJeWcd34bPpALhSW5wnw36tIdEbo/SEOFnd4OqRlJ4QR4gAuFZbnCdNLnz3TUkI+3d0HOIcPD4AuFlbmydNYpSeGKe0BG/Yq++RlJbgVXpiXJiPDADBaWvzpEmMPB6PhmZ0VFSY340edfC4vAsegNu1tXnS7P7UhNhoDUnvENZjDknvoITY8NwzDwDHqi3Nk6ZvluiVnKABqUlhOdaA1CT1Sua1IgCRpa3Mk+bv3BuYmtTiT/SA1CQNDNP/mADQ3NrCPBn29xk1ZOs+n9YUlijgOM1yB4lHB659DknvwBkRgFahNc+TromRJPmq/VpdsFc7fVXySCE92Ycel5bg1dCMjrxGBKBVaa3zpKtiJB34yt3C8krlFftU6Kts8pN9aLv0hDj1TUlQemIcd80BaJVa4zzpuhh9X3lVjfJL9mtPRbWKK6rq/RRbb3SUUuJj1Sk+Vj07tOOTFQC0Ka1lnnR1jI60v8avKn9AAefAFz6F8+PNASASROo8GVExAgC0Tua3dgMAQIwAAOaIEQDAHDECAJgjRgAAc8QIAGCOGAEAzBEjAIA5YgQAMEeMAADmiBEAwBwxAgCYI0YAAHPECABgjhgBAMwRIwCAOWIEADBHjAAA5ogRAMAcMQIAmCNGAABzxAgAYI4YAQDMESMAgDliBAAw938cSHKPAvaDvQAAAABJRU5ErkJggg==\n",

      "text/plain": [

       "<Figure size 400x400 with 1 Axes>"

      ]

     },

     "metadata": {},

     "output_type": "display_data"

    }

   ],

   "source": [

    "import matplotlib.pyplot as plt\n",

    "import networkx as nx\n",

    "\n",

    "G = nx.from_numpy_matrix(np.matrix(J_mat), create_using=nx.DiGraph)\n",

    "pos = {i: (i % Nc, -(i // Nc)) for i in G.nodes()}\n",

    "edge_labels = {(x, y): np.round(J_mat[x, y], 2) for x, y in G.edges()}\n",

    "weights = [x + 1.5 for x in list(nx.get_edge_attributes(G, \"weight\").values())]\n",

    "\n",

    "plt.figure(figsize=(4, 4))\n",

    "nx.draw(\n",

    "    G, pos, node_color=\"lightblue\", with_labels=True,\n",

    "    node_size=600, width=weights, edge_color=\"firebrick\",\n",

    ")\n",

    "nx.draw_networkx_edge_labels(G, pos=pos, edge_labels=edge_labels)\n",

    "plt.show()"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "t0E6JAIIlMUH"

   },

   "source": [

    "We then use the same coupling matrix `J_mat` to obtain the Hamiltonian\n",

    "$H$ for the model we have instantiated above.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 47,

   "metadata": {

    "colab": {

     "base_uri": "https://localhost:8080/",

     "height": 0

    },

    "id": "A_XDopwwlMUH",

    "outputId": "adc70e0d-9aa5-4011-b669-1ef80c5451cd"

   },

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "Hamiltonian =\n",

      "  (0.44013459956570355) [X2 X3]\n",

      "+ (0.44013459956570355) [Y2 Y3]\n",

      "+ (0.44013459956570355) [Z2 Z3]\n",

      "+ (1.399024099899152) [X0 X2]\n",

      "+ (1.399024099899152) [Y0 Y2]\n",

      "+ (1.399024099899152) [Z0 Z2]\n",

      "+ (1.920034606671837) [X0 X1]\n",

      "+ (1.920034606671837) [Y0 Y1]\n",

      "+ (1.920034606671837) [Z0 Z1]\n",

      "+ (1.9997345852477584) [X1 X3]\n",

      "+ (1.9997345852477584) [Y1 Y3]\n",

      "+ (1.9997345852477584) [Z1 Z3]\n"

     ]

    }

   ],

   "source": [

    "import pennylane as qml\n",

    "\n",

    "def Hamiltonian(J_mat):\n",

    "    coeffs, ops = [], []\n",

    "    ns = J_mat.shape[0]\n",

    "    for i, j in it.combinations(range(ns), r=2):\n",

    "        coeff = J_mat[i, j]\n",

    "        if coeff:\n",

    "            for op in [qml.PauliX, qml.PauliY, qml.PauliZ]:\n",

    "                coeffs.append(coeff)\n",

    "                ops.append(op(i) @ op(j))\n",

    "    H = qml.Hamiltonian(coeffs, ops)\n",

    "    return H\n",

    "\n",

    "print(f\"Hamiltonian =\\n{Hamiltonian(J_mat)}\")"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "JD6aZtNilMUH"

   },

   "source": [

    "For the Heisenberg model, a property of interest is usually the two-body\n",

    "correlation function $C_{ij}$, which for a pair of spins $i$ and $j$ is\n",

    "defined as the following operator:\n",

    "\n",

    "$$\\hat{C}_{ij} = \\frac{1}{3} (X_i X_j + Y_iY_j + Z_iZ_j).$$\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 48,

   "metadata": {

    "id": "UVBY4jUglMUH"

   },

   "outputs": [],

   "source": [

    "def corr_function(i, j):\n",

    "    ops = []\n",

    "    for op in [qml.PauliX, qml.PauliY, qml.PauliZ]:\n",

    "        if i != j:\n",

    "            ops.append(op(i) @ op(j))\n",

    "        else:\n",

    "            ops.append(qml.Identity(i))\n",

    "    return ops"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "O0-ZNiNXlMUH"

   },

   "source": [

    "The expectation value of each such operator $\\hat{C}_{ij}$ with respect\n",

    "to the ground state $|\\psi_{0}\\rangle$ of the model can be used to build\n",

    "the correlation matrix $C$:\n",

    "\n",

    "$${C}_{ij} = \\langle \\hat{C}_{ij} \\rangle = \\frac{1}{3} \\langle \\psi_{0} | X_i X_j + Y_iY_j + Z_iZ_j | \\psi_{0} \\rangle .$$\n"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "VhgAsFyJlMUH"

   },

   "source": [

    "Hence, to build $C$ for the model, we need to calculate its ground state\n",

    "$|\\psi_{0}\\rangle$. We do this by diagonalizing the Hamiltonian for the\n",

    "model. Then, we obtain the eigenvector corresponding to the smallest\n",

    "eigenvalue.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 49,

   "metadata": {

    "id": "W-r82wjjlMUH"

   },

   "outputs": [],

   "source": [

    "import scipy as sp\n",

    "\n",

    "ham = Hamiltonian(J_mat)\n",

    "eigvals, eigvecs = sp.sparse.linalg.eigs(ham.sparse_matrix())\n",

    "psi0 = eigvecs[:, np.argmin(eigvals)]"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "Bft7PBtVlMUH"

   },

   "source": [

    "We then build a circuit that initializes the qubits into the ground\n",

    "state and measures the expectation value of the provided set of\n",

    "observables.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 50,

   "metadata": {

    "id": "aMmNdfzAlMUH"

   },

   "outputs": [],

   "source": [

    "dev_exact = qml.device(\"lightning.qubit\", wires=num_qubits) # for exact simulation\n",

    "\n",

    "def circuit(psi, observables):\n",

    "    psi = psi / np.linalg.norm(psi) # normalize the state\n",

    "    qml.QubitStateVector(psi, wires=range(num_qubits))\n",

    "    return [qml.expval(o) for o in observables]\n",

    "\n",

    "circuit_exact = qml.QNode(circuit, dev_exact)"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "7_OTMOuMlMUH"

   },

   "source": [

    "Finally, we execute this circuit to obtain the exact correlation matrix\n",

    "$C$. We compute the correlation operators $\\hat{C}_{ij}$ and their\n",

    "expectation values with respect to the ground state $|\\psi_0\\rangle$.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 51,

   "metadata": {

    "id": "M2SqFe0IlMUI"

   },

   "outputs": [],

   "source": [

    "coups = list(it.product(range(num_qubits), repeat=2))\n",

    "corrs = [corr_function(i, j) for i, j in coups]\n",

    "\n",

    "def build_exact_corrmat(coups, corrs, circuit, psi):\n",

    "    corr_mat_exact = np.zeros((num_qubits, num_qubits))\n",

    "    for idx, (i, j) in enumerate(coups):\n",

    "        corr = corrs[idx]\n",

    "        if i == j:\n",

    "            corr_mat_exact[i][j] = 1.0\n",

    "        else:\n",

    "            corr_mat_exact[i][j] = (\n",

    "                np.sum(np.array([circuit(psi, observables=[o]) for o in corr]).T) / 3\n",

    "            )\n",

    "            corr_mat_exact[j][i] = corr_mat_exact[i][j]\n",

    "    return corr_mat_exact\n",

    "\n",

    "expval_exact = build_exact_corrmat(coups, corrs, circuit_exact, psi0)"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "wRhl6GSvlMUI"

   },

   "source": [

    "Once built, we can visualize the correlation matrix:\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 52,

   "metadata": {

    "colab": {

     "base_uri": "https://localhost:8080/",

     "height": 337

    },

    "id": "elPohZB_lMUI",

    "outputId": "de61464c-fc27-4900-bd66-880b847547f7"

   },

   "outputs": [

    {

     "data": {

      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAFACAYAAACm+Ov/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1tElEQVR4nO3dfVxUZcI//s8ZkIGAQTEKTIQ07C7FzIwsH4AsMLRV0xZtRc0kde/MFtpSfxb0oOTeqLVupS4toK2p+9J0t76KpoD5mPm0bnu7ISYiPiSKICBPM9fvD5u5HWeAc4aBc2A+79frvGrOueacawZnPnNd17nOkYQQAkRERLfRqV0BIiLSJgYEERHZxYAgIiK7GBBERGQXA4KIiOxiQBARkV0MCCIisosBQUREdjEgiIjILgYEtbqoqChIktTuj9GRZWVlQZIkZGVlqV0V0hBNB8SZM2cgSVKTS2hoqNrVtCs0NLRFdSspKcG8efMwYMAAdO7cGR4eHggKCsLIkSORlZWFuro651W2HUhNTYUkScjLy1O7KrKYA0uSJHz11VeNlnvssccs5Vry2vLy8iBJElJTUx3eB9Ht3NWugBy9evXCpEmT7G7r3Llz21amDXzxxRd46aWXcOPGDTzyyCOYNGkS/Pz8cPHiRezatQsvvvgi1qxZg507d6pdVc1YvXo1qqur1a6GDXd3d/zlL3/BqFGjbLb98MMP+O677+Du7o6GhgYVavd/xo4di0GDBiEoKEjVepC2tIuAuO+++1zml9G2bdswadIkdO7cGVu2bMHTTz9ttV0Igc2bNyMjI0OlGmpTjx491K6CXc888wy++uorXL58GQEBAVbbPvvsM+h0OsTGxuLrr79WqYY3+fn5wc/PT9U6kAYJDfvpp58EABEbGyurfFpamgAgZsyY0ei2mTNnWtaVlJSIt99+Wzz22GMiICBAeHh4iJCQEDFr1ixx6dIlu8eora0VS5cuFQMHDhQ+Pj7C29tbPPDAA+J3v/uduHr1qqXO9paUlJQm69/Q0CB69uwpAIhvvvmmybI1NTVWj+vr68WSJUtEv379hKenpzAYDCIqKkr8/e9/t3luZmamACAyMzPF3//+d/HEE08IHx8fERISIoQQYsqUKQKAKCwsFOnp6eKBBx4QHh4eYsqUKZZ9XLp0Sbz22muiV69ewsPDQ3Tt2lU899xz4sSJEzbHi4yMFLf/U7t27Zr44IMPxLBhw0RQUJDo1KmTCAoKEgkJCeLUqVN2n3/7Yq5vY8doyfuSk5MjHn/8ceHl5SX8/f3F5MmTRWlpqc1zGmOuz5dffikAiCVLllhtr6urEwEBAWLEiBFixowZAoDIzc21KvPZZ5+JX/3qVyIkJETo9XrRpUsXERMTI3bt2mVVLiUlpdF/cz/99JMQovm/6a2v3Uzp54k6ng4VEEajUTz55JOWD6bZwYMHRadOncSDDz4oqqurLeu/+OIL4e3tLX71q1+JV199VSQnJ1ue37NnT3Ht2jWr/VdXV4vBgwcLACIsLEzMnj1bvP7662L06NHijjvuEEePHhVlZWUiJSVF+Pn5CT8/P5GSkmJZbv8CuN2OHTsEAPHEE0/Ifo+EEMJkMonRo0cLAKJ3794iOTlZzJw5U3Tp0kUAEEuXLrUqb/4yiIuLE+7u7mLMmDHijTfesHzYzV8mcXFxwt/fXyQkJIg33nhDpKenCyGEOHXqlOjevbsAIGJiYkRycrJISEgQd9xxh/D29hYHDhywOp69L+/9+/cLDw8PERsbK37729+K3//+9+LZZ58Vbm5uwt/fX5w5c8aqvuZ9TJkyxfJ+Llu2rMljOPq+jB07Vnh4eIhx48aJ5ORk8eijjwoAYvDgwbL/Jub6XLhwQTz44IOib9++Vts3btwoAIgNGzY0GhCenp7iscceEy+99JKYO3euSEhIEL6+vkKn04nNmzdbyuXm5lr+ZpGRkVb/5srKyoQQzf9N7QWE0s8TdTztIiB69epl9Y/+1mXr1q1Wzzl37pzo2rWr8Pf3F+fOnRMVFRWiV69eQq/Xi+PHj1uVvXTpkrh+/brNcbOzswUA8f7771utT05OFgBEQkKCaGhosNp27do1q32FhIRY/cKVIzU1VQAQCxYsUPQ8c30jIyNFbW2tZX1RUZG48847hbu7uygsLLSsN38Z6HQ6sWPHDpv9mb9MunfvLoqKimy2P/HEE8LNzU1s27bNav1//vMf4evrK8LDw63WN9aCuHLlis2+d+3aJXQ6nZg+fbrVevOv5MZC1t4xHH1f3N3dxZ49eyzrGxoaRFRUlAAg9u/fb/f4jdXnwoULIj09XQAQ3333nWV7XFyc6Nq1q6itrW00IE6fPm2z3/Pnz4tu3bqJsLAwq/W5ublNtlKb+5vaCwghlH2eqONpFwHR1DJnzhyb523evFkAEFFRUWLSpEkCgPjoo49kH9dkMlm6Iszq6+uFr6+v8PPzE1evXm12H44ExMyZMwUAsWLFCkXPM//KO3jwoM22hQsXCgDi3Xfftay79ZeyPeYvE3vv2ZEjRwQAMW3aNLvPTUpKEgCsupoa6/5pTHh4uAgNDbVa50hAOPq+TJ482aa8edsf//hHWa/h1oC4dOmS6NSpk6WFVlJSItzc3Cz/dhsLiMbMnj1bALBqZckNiMY+B40FhBAt/zxR+9UuBqljY2Oxbds22eVHjx6NmTNnYsWKFQCAuLg4vPrqq3bLbtq0CStXrsSRI0dQVlYGo9Fo2Xb+/HnL/588eRLXr1/HU089hS5dujj4SlrH0aNHcccddyAiIsJmW3R0NADg2LFjNtvslW9u+4EDBwAAly5dsnviwMmTJy3/7du3b5P7z8vLw4cffoiDBw+itLTU6kweDw+PJp8rh6PvyyOPPGKzrnv37gCAa9euKa7HXXfdhZEjR2LdunVYtmwZsrOzYTQaMW3atCafd/r0aaSlpWHXrl0oKSlBbW2t1fbz588jJCREUV2a+5vbo+TzRB1LuwgIR4wdO9byD/qVV16xW2bJkiV4/fXXERAQgJiYGHTv3h1eXl4AgA8//NDqA1leXg4AuOeee1qtzoGBgQBuzoFQoqKiAsHBwXa3mU9brKiosNl29913N7lfe9uvXr0KAPj666+bPPOmqqqqyX3/7W9/Q3x8PHx8fBAbG4vQ0FDccccdlslaRUVFTT5fDkffF4PBYLPO3f3mR+XWHxBKTJs2DZs3b8bGjRuRmZmJRx55BP369Wu0/KlTpxAREYGKigpER0fj2WefhcFggE6nQ15eHvLz820CQ47m/uaNkfN5oo6nQwbEtWvXkJiYCG9vbxiNRsyePRtHjx6Fr6+vpUxDQwPee+89BAUF4dixY7jrrrss24QQ+MMf/mC1T/N8C6Vf3koMHjwYALBz5068++67sp9nMBjw888/29128eJFS5nbNTfz2N52836WL1/eoi+K1NRUeHp64vDhwwgLC7Patm7dOof3eytH35fWEBcXh6CgILz55psoKSnBJ5980mT5ZcuWoaysDGvWrLGZAzRz5kzk5+c7VA9HZpvL+TxRx6TpmdSOevnll3H27Fl89NFH+J//+R8UFhbiv//7v63KlJaWory8HI8//rhVOADA999/jxs3blitu//++2EwGHDo0CGUlZU1Wwc3NzfFvzajo6PRs2dP7Nu3D7m5uU2WvfXX48MPP4zq6mp89913NuXMs3P79++vqC6NeeyxxwAA+/fvb9F+CgsL8cADD9iEw4ULF3D69Gmb8m5ubgCU/YJvy/elOW5ubpg8eTJKSkrg6emJiRMnNlm+sLAQwM3unVsJIbB37167+wccb+E0Rc7niTqmDhcQn332Gf72t7/h+eefx0svvYRXXnkFo0aNwpo1a7B27VpLubvuugteXl44cuSI1QzcsrIyzJ4922a/7u7umDFjBsrLyzFnzhybD2J5eTkqKystj/39/VFaWoqamhrZdXdzc8PHH38MnU6HX//619i1a5fdcv/4xz8wfvx4y+MpU6YAAObNm4f6+nrL+uLiYixduhTu7u74zW9+I7seTYmIiMBjjz2GL774AuvXr7fZbjKZZP26DQkJwalTp3Dp0iXLupqaGsyaNcvqNZj5+/sDuPma5GrL90WOpKQkfPnll8jJyWn2CgDmsYU9e/ZYrf/ggw/wr3/9y6a8I++PHHI/T9QxtYsuplOnTjU5k3ru3Lnw9PTEjz/+iDlz5iA4OBirVq2ybP/LX/6Cfv36YdasWXj88cdx7733QqfT4be//S2WLFmChx56CM8++ywqKiqwdetWhISEoFu3bjbHeffdd3HgwAGsWbMGBw4cwDPPPAO9Xo/Tp09j27Zt2LNnj+UX6ZNPPonvv/8ezzzzDIYOHQoPDw8MGzYMw4YNa/K1jhgxAmvWrMH06dMxfPhwDBw4EI8//jh8fX1x6dIl5OXlobCwEE899ZTlOQkJCdi0aRO2bNmCfv36YdSoUaiqqsL69etx9epVLFmyBD179lT2pjfhiy++QHR0NCZMmIAPP/wQAwYMgJeXF86ePYv9+/fj8uXLzQbj7NmzMXv2bDz88MMYP348GhoasGPHDggh8NBDD+H48eNW5aOjoyFJEubPn48ffvgBfn5+6Ny5c5PdXG39vjTnrrvuwpgxY2SVnTlzJjIzMzFu3Dj8+te/RteuXXHgwAEcOXIEI0eOtBn/+a//+i9069YN69atg16vR/fu3SFJEmbPnu3wDGklnyfqoFQ+i6pJck5zBSDKyspEbW2tGDBggNDpdCI/P99mX9u3bxeSJIlBgwaJ+vp6IcTN2awLFy4UYWFhQq/Xix49eojk5GRx/fr1Rk9TrampEenp6aJ///7Cy8tL+Pj4iAcffFAkJydbJiUJIcT169dFYmKiCAoKEm5ubrJmUt/q3Llz4s033xQPP/ywMBgMwt3dXdx9991ixIgRIjMzU9TV1VmVr6+vF+np6SI8PFzo9Xrh6+srIiMjxZYtW2z23dQpjUL83ymR5lm49ly9elUsWLBA9O3b1/I+hIWFiRdeeEFs2rTJqmxjk9hWrFgh+vTpIzw9PUVgYKB46aWXxM8//9zoabFZWVmW1wcFM6md8b40dxrp7W49zbU5jZ3mmpubKwYPHix8fX1F586dRVxcnDh8+HCjp/weOHBAREZGCl9f30ZnUjf2N739tTvyeaKORxJCiDZLIyIiajc63BgEERE5BwOCiIjsYkAQkcv6/PPPMWPGDAwcOBB6vd7hu+qZTCYsX74c4eHh8PLyQkBAACZOnGj3lO32hAFBRC5rwYIFWLVqFYqKilp0s6QZM2bg1VdfhRACr776KkaMGIFNmzbh0UcfRUFBgRNr3LYYEETksjIyMnDmzBlcvnwZM2fOdGgfubm5yMjIwLBhw3DkyBEsXrwYa9aswebNm3H16tV2fWmSdjEPgoioNdw6n8hRf/7znwEA7733ntVFJp955hlERUVh+/btOHv2rGbvetgUVQPCZDLh/Pnz8PX1degaMUSkXUIIXL9+Hd26dYNO13E7K/Ly8uDt7W25ltqtYmNjLRdXTEhIUKF2LaNqQJw/f77Rq20SUcdQXFxsuVy6UjU1Nairq5NdXghh82NTr9dDr9c7dPzmVFVV4cKFC+jbt6/leli3Ml9rrL2OQ6gaEOarQbo9+GtIbp3UrEqb+/c4b7WroIpv/5indhVU4/X//qF2FdpUdVUlpg9/xOGrvtbU1MDL1x9ouNF84V/4+PhYXRMNAFJSUpq8VE9LmG8D0NjlTMxXCzaXa29UDQhz0ktunSC5tfwGMe2Jr6drvV6zO3S2v7JchZePa14e29Hu47q6OqDhBjr1nQjI+QFprEflv75AcXGx1WXcW6v14Ao4SE1EmiZ18pT1A1L88uPDYDC02X0+zC2HxloI5htSOXrBRLUxIIhI0ySdGyQ5LU/R9q1Tb29vBAUF4aeffoLRaLQZhzCPPdx+35P2ouOeWkBEHYIkuVlCoslFUqf7MjIyElVVVXZv5JSTkwMAzV7mX6sYEESkaZKbDpKbm4yldb/OSktLcfLkSZSWllqtf/nllwEAb731ltUZV1u3bkVeXh5iYmIsN4Bqb9jFRESappPZxSQcOAEiIyPDcte+EydOWNaZb0k7ZMgQTJ8+HQDwpz/9Ce+8847NWVHR0dGYPn06MjIyMGDAAIwcORIXLlzA+vXr4e/vj+XLlyuul1YwIIhI02SPQTgQEHv27EF2drbVur1791p1F5kDoikrV65EeHg4Vq1ahY8++gg+Pj4YO3YsFi5ciF69eimul1aoesOgiooK+Pn5wT38Ny53muvZCa45DyIvfafaVVCNV55rvfbqyut4YdD9KC8vd+isIvP3g+/QJEjuzZ+qKhpqcf3bpQ4fj2yxBUFEmibpdJDkXKqjA1/OQy0MCCLStNbsYqKmMSCISNNutiDkBARbEM7GgCAiTTPPg2i+IFsQzsaAICJt+2WeQ3OEiQHhbAwIItI0uWMQsloZpAgDgog0zc3dAzr35k+Dl4SpDWrjWhgQRKRpcgepZZ0KS4owIIhI09jFpB4GBBFpGgNCPQwIItI0BoR6GBBEpGly50GodT+IjowBQUSaJsmcByGnDCnDgCAiTeNZTOphQBCRpnEMQj0ORe6hQ4cQFxeHzp07w9vbG4MGDcKGDRucXTciInn3o5Z7xVdSRHELIjc3F7GxsfD09MSECRPg6+uLjRs3Ij4+HsXFxUhOTm6NehKRi9LpJOh0koyCMsqQIooCoqGhAYmJidDpdNi9ezf69+8PAHj77bcRERGB+fPnY/z48e32Bt1EpD2SToIk48tfThlSRlEX065du1BYWIgXXnjBEg4A4Ofnh/nz56Ours7m/q5ERC0hSZLshZxLUQsiLy8PABATE2OzLTY2FgCQn5/f8loREf1CktnFJNiCcDpFAVFQUAAACAsLs9kWGBgIHx8fSxl7amtrUVtba3lcUVGh5PBE5IIkSWYXE1sQTqeoi6m8vBzAzS4lewwGg6WMPWlpafDz87MswcHBSg5PRC7IPAYhZyHnatOZJfPmzUN5ebllKS4ubsvDE1E7pJMk2Qs5l6IuJnPLobFWQkVFBbp06dLo8/V6PfR6vZJDEpGL41lM6lHUgjCPPdgbZ7h48SIqKyvtjk8QETmKXUzqURQQkZGRAIDt27fbbMvJybEqQ0TkDOaJcnIWci5FATF8+HD07NkTa9euxbFjxyzry8vLsWjRInh4eGDy5MnOriMRuTBJJ38h51I0BuHu7o6MjAzExsZi2LBhVpfaKCoqQnp6OkJDQ1upqkTkiuROguNprs6n+FpM0dHR2LNnD1JSUrB+/XrU19cjPDwcixcvRnx8fGvUkYhcmE4HmRPl2qAyLsahy31HRERg69atzq4LEZENnsWkHt4Pgog0jTOp1cOAICJNkzsJTjAgnI4BQUSaJrnroHNvfoDBJKMMKcOAICJNkzvHgfMgnI8BQUSaxtNc1cOAICJNkzsJjhPlnI8BQUSaxi4m9TAgiEjTOA9CPQwIItI0jkGohwFBRJrGLib1MCCISNM4k1o9DAgi0jQ3nQQ3WRfrY0A4GwOCiDRNJzMgTAwIp2NAEJGmyW1BMCCcjwFBRJrGgFAPA4KINI0BoR5OTiciTXPXAe46Scbi2P4PHTqEuLg4dO7cGd7e3hg0aBA2bNgg+/lZWVmWuRr2lry8PMcqpgFsQRCRprVmCyI3NxexsbHw9PTEhAkT4Ovri40bNyI+Ph7FxcVITk6Wva/Ro0ejf//+NutDQ0MV10srGBBEpGlyz2IyKgyIhoYGJCYmQqfTYffu3ZYv97fffhsRERGYP38+xo8fj5CQEFn7GzNmDKZOnaqoDlrHLiYi0jQ3SQc3nYxF4eVcd+3ahcLCQrzwwgtWv/z9/Pwwf/581NXVITs728mvpn1hC4KINE1uF5OcMrcyjw3ExMTYbIuNjQUA5Ofny97f0aNHceXKFTQ0NCA0NBRPPfUUunbtqqhOWqOJgPj3OG/4enqoXY021WNdldpVUMXf8naqXQXV3IgarnYV2tQNk9Ep+2mtgCgoKAAAhIWF2WwLDAyEj4+PpYwcf/zjH60ee3l5ISUlBW+++aaiemkJu5iISNPMASFnAYCKigqrpba21u5+y8vLAdzsUrLHYDBYyjTl3nvvxfLly/Hjjz+iuroa586dw+rVq+Hv74+5c+di+fLlDr5y9TEgiEjT3CRJ9gIAwcHB8PPzsyxpaWmtWr/IyEi88sorCAsLg5eXF+655x4kJCQgJycHnp6eSE1NRUNDQ6vWobVooouJiKgxcs9iMl/uu7i4GAaDwbJer9fbLW9uOTTWSqioqECXLl2UVteiT58+GDJkCL755hv87//+L8LDwx3el1rYgiAiTVPaxWQwGKyWxgLCPPZgb5zh4sWLqKystDs+ocSdd94JAKiqap9jjgwIItI0ebOoby5KREZGAgC2b99usy0nJ8eqjCOMRiO+//57AJA9l0JrGBBEpGke7jrZixLDhw9Hz549sXbtWhw7dsyyvry8HIsWLYKHhwcmT55sWX/hwgWcPHnSpkvq8OHDNvs2Go2YO3cuTp06hejoaAQFBSl70RrBMQgi0jQ3SeZprgrvKOfu7o6MjAzExsZi2LBhVpfaKCoqQnp6utVlMubNm4fs7GxkZmZazZgeOHAg+vXrh379+uGee+7B1atXkZ+fjx9//BHdu3dHRkaGonppCQOCiDRN6SC1EtHR0dizZw9SUlKwfv161NfXIzw8HIsXL0Z8fLysfSQnJ+PAgQPYsWMHrl69Cg8PD9x3331YsGABkpKSWjTQrTYGBBFpWmtNlDOLiIjA1q1bmy2XlZWFrKwsm/Xp6ekOHbc9YEAQkaa1dkBQ4xgQRKRpbjp5X/5uPOXG6RgQRKRpbEGohwFBRJrGgFAPA4KINK01z2KipjEgiEjTbr0QX3PlyLkYEESkaTpJgk7Gl7+cMqQMA4KINM0NgJuM7363Vq+J62FAEJGm6XSSrPEFjkE4HwOCiDSNYxDqYUAQkaZxDEI9DAgi0jSdJG8Mgj1MzseAICJN4xiEehgQRKRp7GJSDwOCiDTNTWYXk5wypAwDgog0jS0I9Si+QO7nn3+OGTNmYODAgdDr9ZAkye5NNIiInMF8sT45CzmX4hbEggULUFRUhDvvvBNBQUEoKipqjXoREQFgC0JNilsQGRkZOHPmDC5fvoyZM2e2Rp2IiCzMYxByFnIuxS2Ip556qjXqQURklySzBSGxBeF0HKQmIk3jDYPU06YBUVtbi9raWsvjioqKtjw8EbVDOsibJc1bUjtfm76naWlp8PPzsyzBwcFteXgiaoc66XSyF3KuNn1H582bh/LycstSXFzclocnonbITSd/Iedq0y4mvV4PvV7flockonZOJ8k7hZVDEM7HQWoi0jSdzPtBcB6E8zEgiEjTOFFOPQwIItI0ueMLHINwPsUBkZGRgT179gAATpw4YVmXl5cHABgyZAimT5/uvBoSkUtjC0I9igNiz549yM7Otlq3d+9e7N271/KYAUFEziJJNxc55ci5FAdEVlYWr95KRG1GBwk6yGhByChDynAMgog0jS0I9TAgiEjTbs6DkFeOnIsBQUSaxhaEehgQRKRpHINQDwOCiLRNZguC+eB8DAgi0jSOQaiHAUFEmiZBXuOA+eB8DAgi0jTOpFYPA4KINE2CzLOYWr0mrocBQUSapoO8O5vxWn3Ox4AgIk2TJAmSjCaEnDKkDAOCiDSNZzGphwFBRJrGmdTqYUAQkaZxDEI9DAgi0jSOQaiHAUFEmsYxCPWwVUZEmifJWBx16NAhxMXFoXPnzvD29sagQYOwYcMGRfuora3Fu+++i7CwMHh6eqJbt254+eWX8fPPP7egZupjC4KINM1NJ8FNRvNATpnb5ebmIjY2Fp6enpgwYQJ8fX2xceNGxMfHo7i4GMnJyc3uw2QyYfTo0cjJycGgQYMwbtw4FBQUICMjAzt37sSBAwcQEBCguG5awBYEEWmauYtJzqJEQ0MDEhMTodPpsHv3bqxatQpLlizB8ePH0bt3b8yfPx9FRUXN7ic7Oxs5OTmYOHEi9u3bhw8++AAbN27EJ598gtOnT2PBggUOvnL1MSCISNPkdC850s20a9cuFBYW4oUXXkD//v0t6/38/DB//nzU1dUhOzu72f38+c9/BgCkpaVZDZTPmDEDPXv2xF//+lfcuHFDYe20gQFBRJpmvlifnEWJvLw8AEBMTIzNttjYWABAfn5+k/uoqanBwYMHcf/99yMkJMRqmyRJePrpp1FVVYXvv/9eUd20QhNjEN/+MQ936NzUrkab+lveTrWroIrnJ/1/aldBNWdfH652FdrU9Zo6IPXHFu9H6US5iooKq/V6vR56vd6mfEFBAQAgLCzMZltgYCB8fHwsZRpTWFgIk8lkdx+37rugoABDhw5t9jVoDVsQRKRpkhCyFwAIDg6Gn5+fZUlLS7O73/LycgA3u5TsMRgMljKNkbOPW8u1N5poQRARNUqYbi5yygEoLi62fDEDsNt6IHnYgiAiTZOESfYC3PzVfuvSWECYf/U39uu+oqKi0ZaBkn3cWs4ZDh8+jNmzZ6NPnz7o0qUL9Ho9QkJCMHHiROTk5DjtOAADgoi0ztyCkLMocOv4wO0uXryIysrKRscWzHr27AmdTtfoWEVT4xxK1dfXY+bMmXj00UexYsUK9OjRA1OmTMGrr76KPn364Msvv8SIESPwpz/9qcXHMmMXExFpmxA3FznlFIiMjERaWhq2b9+OCRMmWG0z/xKPjIxsch9eXl6IiIjAgQMHUFRUZHUmkxACO3bsgLe3NwYOHKiobrczGo0YN24c/vGPf+DJJ5/E6tWrcc8991iVOXfuHGbPno3Q0NAWHetWbEEQkba1Ugti+PDh6NmzJ9auXYtjx45Z1peXl2PRokXw8PDA5MmTLesvXLiAkydP2nQnvfzyywCAefPmQdwSUitXrsTp06fxm9/8Bl5eXg688P+TmppqCYecnBybcACA7t27Y9OmTXj66adbdKxbsQVBRJp28wyl5r/8JYUtCHd3d2RkZCA2NhbDhg2zutRGUVER0tPTrX6Nz5s3D9nZ2cjMzMTUqVMt66dMmYL169fjiy++wE8//YTIyEicOnUKmzZtwr333ov3339fUb1uV1hYiLS0NPj5+eGvf/0r3N0b/9qWJMmpg/JsQRCRtrVSCwIAoqOjsWfPHgwePBjr16/Hp59+irvvvhvr1q2TdR0mANDpdNiyZQtSU1Nx+fJlLFu2DHv37sVLL72E/fv3t/g6TEuXLoXRaMSsWbMQGBjYon0pJQmhMHadyHyWwGf+vV1uopwXJ8q5nLMTvNWuQpu6XlOH3qmrUF5ebnXaqVzm74fLZ36EweAro/x1BIT2dvh4WhUcHIxz587hxIkT6Nu3r+znJSQkwGAw4OOPPwYAzJo1C9XV1bIuH2LGLiYi0jaF8yA6ksuXL+PcuXPw9fVFnz59FD33o48+gqenp+WxeVxFCQYEEWmbMAEm1wyI0tJSAEBAQIDiO+b5+/tbPe7SpYvi43MMgog0TelEuY7E2/tmt+SlS5egZDTgP//5DyRJQllZGQCgpKQEkiTh3Llzio7PgCAibWvFQWqtCw4Oxt13342qqiocOnSoybKmW1pZx48fR3BwsKXVcPz4cfj7+6N79+6Kjs+AICJtM0+Uk7N0MJIk4ZVXXgFw8/4S9loAN27cwKeffmp1UcJ//vOfeOihhyyPjx8/jn79+ik+PscgiEjbXHiQGgDmzp2LI0eO4Msvv8T999+PESNGoFevXjAajfjxxx/x7bffory8HGvXrrU85/ZAcDQg2IIgIk0zT5Rrful4LQjg5oS+TZs2Yd26dRg2bBi+/fZbLFu2DGvWrMG5c+cQHx+Pr7/+GuPGjbM8hy0IInINLt6CMIuPj0d8fHyz5a5du4azZ89aAqGmpgYFBQVWgSEXA4KItI0Bocg///lPeHl5Wa4g+69//QsAFM+jABgQRKR1DAhFjh8/jj59+sDNzc3y+L777nPogoEcgyAiTXPleRCOmD17ttUpsYcOHcKjjz7q0L4YEESkbUaj/IUsampqcPjwYWzatAmxsbEO7YMBQUTa5sIT5Vrik08+wfDhwzFu3DibGyLJxTEIItI0ud1H7GKylpSUhKSkpBbtgwFBRNrGQWrVMCCISNuEkBkQHXOinJoUjUGUlJTgww8/RExMDHr06AEPDw8EBgZi3LhxOHjwYGvVkYhcmTACJhmL4CC1sylqQSxfvhyLFy9Gr169EBMTg4CAABQUFGDz5s3YvHkz1q5dK2umHxGRXMJkgpBxPwg5ZUgZRQERERGBvLw8REZGWq3/9ttvMXz4cMyaNQtjxoxx6k2zicjFmVsIcsqRUynqYnruuedswgEAhg4diujoaJSVleHEiRNOqxwRkazuJbkhQoo4bZC6U6dON3foznFvInIeYTRCyJgEJ6cMKeOUb/OzZ8/im2++QVBQEMLDwxstV1tbi9raWsvjiooKZxyeiDoyk8x7UnMMwulaPJO6vr4eCQkJqK2txeLFiy0XiLInLS0Nfn5+liU4OLilhyeijs5kktnFxIBwthYFhMlkwtSpU7F7924kJiYiISGhyfLz5s1DeXm5ZSkuLm7J4YnIBQiTUfZCzuVwF5PJZMK0adOwdu1aTJo0CStWrGj2OXq9nmc4EZEyQmYXE2dSO51DAWEymfDiiy9i9erVmDhxIrKysqDT8bp/ROR8clsHbEE4n+KAuDUc4uPjsWbNmibHHYiIWoTzIFSjKCDM3UqrV6/G888/j88//5zhQESti2cxqUZRQLz77rvIzs6Gj48Pevfujffff9+mzJgxY9C/f39n1Y+IXBznQahHUUCcOXMGAFBZWYmFCxfaLRMaGsqAICLnMZ/mKqccOZWigMjKykJWVlYrVYWIyA6OQaiG18UgIk3j1VzVw4AgIm1jC0I1DAgi0jYhMyB4wyCnY0AQkaaxi0k9DAgi0jaexaQaBgQRaRvHIFTDgCAiTRP19RD1zV+xQdTXt0FtXAsDgoi0jS0I1TAgiEjTeDVX9TAgiEjTeBaTehgQRKRpwiQgjHICQrRBbVwLA4KINE0YTfICQkYZUoYBQUSaxi4m9TAgiEjT2IJQDwOCiDSNAaEendoVICJqijAaYZKxqHFHuYqKCiQlJSEkJAR6vR6hoaH4/e9/j8rKSkX7kSSp0WXq1KmtU3kZ2IIgIk0TQuYYhGjbFkRVVRUiIyNx7NgxxMTEYOLEiTh69CjS09ORn5+P3bt3w9PTU/b+QkJC7IaBmnfoZEAQkaZptYvpD3/4A44dO4Y333wTH3zwgWX93LlzsXjxYixbtgzz5s2Tvb/Q0FCkpqa2Qk0dxy4mItI0c0DIWdqsTkIgIyMDPj4+eOutt6y2vfXWW/Dx8UFGRkab1ae1sAVBRJomTELmaa5tN1GuoKAA58+fR2xsLLy9va22eXt7Y/DgwcjJyUFxcTGCg4Nl7fPatWtYtWoVSktL4e/vj8GDByM8PLw1qi8bA4KINM1kNMEko3Ugp4yzFBQUAADCwsLsbg8LC0NOTg4KCgpkB8Tx48cxY8YMq3UjRoxAdnY27rrrrpZV2EHsYiIiTVPaxVRRUWG11NbWOr1O5eXlAAA/Pz+72w0Gg1W55iQnJ2Pfvn0oLS1FRUUF9u3bh2eeeQbbtm3DqFGjYFThDC1AIy0Ir//3D3j5+KpdjTZ1I2q42lVQxdnXXfN1A0CPdVVqV6FNCWOdk/ajbJD69l/sKSkpjQ7+JicnKwqQOXPmNNpqaIn09HSrx48//ji++uorPPnkk8jPz8eWLVvw3HPPOf24zdFEQBARNUbpaa7FxcWWX/AAoNfrG33OypUrUVUlP7jHjx+PsLAwS8uhsRZCRUUFgMZbGHLodDokJiYiPz8fe/fuZUAQEd1OaQvCYDBYBURTlE5oMzO3IsxjEbdrboxCrjvvvBMAFIWYMzEgiEjTtDgPIiwsDN26dcPevXtRVVVldSZTVVUV9u7di3vvvVf2AHVjDh48CODmHAk1cJCaiDTNZDLJXtqKJEmYPn06Kisr8d5771lte++991BZWYnExESr9dXV1Th58iTOnj1rtf7EiROot3M/7X379mHx4sXo1KkTnn/+eee/CBnYgiAiTdNiCwIA3njjDWzZsgWLFy/G0aNHMWDAABw5cgTbt2/Ho48+itdee82q/HfffYfo6GhERkYiLy/Psn7JkiX4+uuvMWTIEAQHB6NTp0744YcfsH37dkiShI8//hi9evVq09dmxoAgIk27GRAy7kndxgHh7e2N/Px8pKamYuPGjcjNzUVQUBCSk5ORkpICLy8vWfsZPXo0rl27huPHj2PHjh2oq6tDYGAgJkyYgNdeew0RERGt/Eoax4AgIk3T8g2D/Pz8sGzZMixbtqzZslFRURDCdrb32LFjMXbs2NaoXosxIIhI04RJZhcT7yjndAwIItI2uRfi4w2DnI4BQUSaZqxvgFGSZJUj52JAEJGmafUsJlfAgCAiTRNGAWFs/lLecsqQMgwIItI0k0nm5b45SO10DAgi0rSbNwyS0YJowxsGuQoGBBFpmskImHTNf/mb1LllQofGgCAiTRNGE4SOg9RqYEAQkaYJo4CQ0YLgILXzMSCISNNMRiGzi4kB4WwMCCLSNHYxqYcBQUSaZhICJhlnKJnsXAiPWoYBQUTaZhQQkowvf3YxOR0Dgog0zWQ0wSTJmCjHLianY0AQkaYJmS0InsXkfAwIItI0BoR6dEoK19TUICkpCcOGDUO3bt3g6emJwMBADB48GJmZmXZvvE1E1BImo0n2Qs6lKCAqKyvx6aefQpIkjBw5EklJSRg7dixKSkowbdo0jBo1ihfMIiKnEkJYrsfU5MKzmJxOUReTv78/ysvL4eHhYbW+oaEBTz/9NLZv346tW7di5MiRTq0kEbkuk1HABE6UU4OiFoROp7MJBwBwd3e33HT71KlTzqkZERHM94MwyVgYEM7mlEFqk8mEbdu2AQD69u3rjF0SEQH4JSBktCAYEM7nUEDU1dVh0aJFEELgypUr2LlzJ06ePIkXX3wRw4cPb/R5tbW1qK2ttTyuqKhw5PBE5ELYxaQehwPinXfesTyWJAmvv/460tLSmnxeWlqa1fOIiJojTCYISZJVjpxL0RiEmY+PD4QQMBqNKC4uxscff4yMjAxERUU12SqYN28eysvLLUtxcbHDFSci12AyCtkLOZdDAWF5sk6H7t27Y9asWVi1ahX27t2LhQsXNlper9fDYDBYLURETREm8ctAdTMLbznqdE6bSR0TEwMAyMvLc9YuiYgAowlCNN/FBHYxOZ3TAuL8+fMAgE6dOjlrl0REMNaZYNQ1HxBGBoTTKepi+ve//43q6mqb9dXV1UhKSgIAxMXFOadmRET45X4QMhdyLkUtiA0bNmDp0qUYMmQIQkNDYTAYUFJSgq1bt+LKlSsYOnQofve737VWXYnIBRmFgFHGl7+cMqSMooAYNWoUzp8/j3379mH//v2orKyEn58f+vXrhwkTJmDatGlwd+cFYonIeYxC3r2AeBKT8yn6Nh84cCAGDhzYWnUhIrLBFoR6+HOfiDSNLQj1MCCISNNMMlsQHKR2PgYEEWmaETJbEK1eE9fDgCAiTTMKAaOMi/VxDML5GBBEpGlGIa91wDEI52NAEJGmMSDUw4AgIk1jF5N6GBBEpGkmmS0IXszV+RgQRKRpbEGohwFBRJrGMQj1MCCISNNuBoScFkQbVMbFMCCISNPYglAPA4KINI1jEOphQBCRpgkAcu4Vx3hwPgYEEWkaWxDqUXTLUSKitma+3LecpS0dO3YM8+fPR2xsLAICAiBJEqKiohze36FDhxAXF4fOnTvD29sbgwYNwoYNG5xXYQewBUFEmqbVFsTmzZuRlpYGDw8P9O7dG6WlpQ7vKzc3F7GxsfD09MSECRPg6+uLjRs3Ij4+HsXFxUhOTnZizeVjC4KINE2rLYjnn38ehw8fRmVlJXbs2OHwfhoaGpCYmAidTofdu3dj1apVWLJkCY4fP47evXtj/vz5KCoqcmLN5WNAEJGmmW85KmdpS3369MGAAQPQqVOnFu1n165dKCwsxAsvvID+/ftb1vv5+WH+/Pmoq6tDdnZ2C2vrGHYxEZGmdfRrMeXl5QEAYmJibLbFxsYCAPLz89uyShYMCCLStBvCKKt1UPfLybAVFRVW6/V6PfR6favUzRkKCgoAAGFhYTbbAgMD4ePjYynT1lQNCPHLH726qlLNaqjihsk1b5B4vaZO7SqoRhhd67ULY/3N/zrY9ePh4YHAwECsvlgi+zk+Pj4IDg62WpeSkoLU1FSH6tAWysvLAdzsUrLHYDBYyrQ1VQPi+vXrAIDpwx9RsxrUllJ/VLsG1MauX7/e6JdfUzw9PfHTTz+hrk5+sAohIEmS1bqmWg/Jycmora2Vvf85c+bY/aXfUakaEN26dUNxcTF8fX1t/qitraKiAsHBwSguLobBYGjTY6uJr5uvu60IIXD9+nV069bN4X14enrC09PTibWytnLlSlRVVckuP378eKcHhDk8G2slVFRUoEuXLk49plyqBoROp0P37t3VrAIMBoNLfWGY8XW7FrVetyMth7ZUWal+97Y5cAoKCvDII9a9KRcvXkRlZSUiIiLUqBpPcyUiUlNkZCQAYPv27TbbcnJyrMq0NQYEEVEbqK+vx8mTJ1FYWGi1fvjw4ejZsyfWrl2LY8eOWdaXl5dj0aJF8PDwwOTJk9u4tje57Gmuer0eKSkpmj79rTXwdfN1k3OcPHkSH3zwAQDgxo0blnVTp061lMnKyrL8f0lJCR544AGEhITgzJkzlvXu7u7IyMhAbGwshg0bZnWpjaKiIqSnpyM0NLQNXpEtSTh6DhoRkQvLy8tDdHR0k2Vu/Xo9c+YM7r33XpuAMPvuu++QkpKCffv2ob6+HuHh4UhKSkJ8fLyzqy4bA4KIiOziGAQREdnFgCAiIrsYEEREZJfLBYQW79rU2j7//HPMmDEDAwcOhF6vhyRJVmdXdEQlJSX48MMPERMTgx49eliu6zNu3DgcPHhQ7eq1mpqaGiQlJWHYsGHo1q0bPD09ERgYiMGDByMzMxP19fVqV5HaEZcapG7srk3mU8nUumtTawsNDUVRURHuvPNOeHt7o6ioCJmZmVan43U0c+fOxeLFi9GrVy9ERUUhICAABQUF2Lx5M4QQWLt2rapnh7SW0tJSBAcHIyIiAr1790ZAQADKysqwdetWFBUVISYmBlu3boVO53K/DckRwkXU19eLXr16Cb1eL44ePWpZf+3aNdG7d2/h4eEhzpw5o14FW9GOHTssry0tLU0AEJmZmepWqpVt3LhR5OXl2azfvXu36NSpk+jSpYuoqalRoWaty2g0itraWpv19fX1IioqSgAQX331lQo1o/bIZX5GaPmuTa3tqaeeQkhIiNrVaFPPPfec3csTDB06FNHR0SgrK8OJEydUqFnr0ul08PDwsFnv7u6OsWPHAgBOnTrV1tWidsplAkLLd22itmW+RaS7u+tcSMBkMmHbtm0AgL59+6pcG2ovXOYTouW7NlHbOXv2LL755hsEBQUhPDxc7eq0mrq6OixatAhCCFy5cgU7d+7EyZMn8eKLL2L48OFqV4/aCZcJCC3ftYnaRn19PRISElBbW4vFixfDzc1N7Sq1mrq6OrzzzjuWx5Ik4fXXX0daWpqKtaL2xmW6mMi1mUwmTJ06Fbt370ZiYiISEhLUrlKr8vHxgRACRqMRxcXF+Pjjj5GRkYGoqCibezYTNcZlAkLOXZu0fnMTcozJZMK0adOwdu1aTJo0CStWrFC7Sm3GfFOuWbNmYdWqVdi7dy8WLlyodrWonXCZgLj1rk23M9+1yZXuNesqTCYTXnzxRWRnZ2PixInIyspy2TkA5hM0zCdsEDXHZT4pWr5rE7UOczisXr0a8fHxWLNmTYced2jO+fPnAfzfWVxEzXGZgNDyXZvI+czdSqtXr8bzzz+Pzz//3CXC4d///jeqq6tt1ldXVyMpKQkAEBcX19bVonaKl9pwgUttZGRkYM+ePQCAEydO4MiRIxg8eDDuu+8+AMCQIUMwffp0NavodKmpqXjnnXfg4+ODOXPm2J3zMGbMGKtJkx1Bamoqli5diiFDhiA0NBQGgwElJSXYunUrrly5gqFDhyInJwdeXl5qV5XaA3Uncre9gwcPihEjRgiDwSC8vLxERESEWLdundrValVTpkwRABpdpkyZonYVna6514wOermRQ4cOicTERNGnTx/RuXNn4e7uLrp27Sqio6PFypUrRX19vdpVpHbEpVoQREQkn8uMQRARkTIMCCIisosBQUREdjEgiIjILgYEERHZxYAgIiK7GBBERGQXA4KIiOxiQBARkV0MCCIisosBQUREdjEgiIjILgYEERHZ9f8Dp69DnhC1TLQAAAAASUVORK5CYII=\n",

      "text/plain": [

       "<Figure size 400x400 with 2 Axes>"

      ]

     },

     "metadata": {},

     "output_type": "display_data"

    }

   ],

   "source": [

    "fig, ax = plt.subplots(1, 1, figsize=(4, 4))\n",

    "im = ax.imshow(expval_exact, cmap=plt.get_cmap(\"RdBu\"), vmin=-1, vmax=1)\n",

    "ax.xaxis.set_ticks(range(num_qubits))\n",

    "ax.yaxis.set_ticks(range(num_qubits))\n",

    "ax.xaxis.set_tick_params(labelsize=14)\n",

    "ax.yaxis.set_tick_params(labelsize=14)\n",

    "ax.set_title(\"Exact Correlation Matrix\", fontsize=14)\n",

    "\n",

    "bar = fig.colorbar(im, pad=0.05, shrink=0.80    )\n",

    "bar.set_label(r\"$C_{ij}$\", fontsize=14, rotation=0)\n",

    "bar.ax.tick_params(labelsize=14)\n",

    "plt.show()"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "X7AxsPAnlMUI"

   },

   "source": [

    "Constructing Classical Shadows\n",

    "==============================\n"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "NTPT4XUFlMUI"

   },

   "source": [

    "Now that we have built the Heisenberg model, the next step is to\n",

    "construct a\n",

    "`classical shadow <tutorial_classical_shadows>`{.interpreted-text\n",

    "role=\"doc\"} representation for its ground state. To construct an\n",

    "approximate classical representation of an $n$-qubit quantum state\n",

    "$\\rho$, we perform randomized single-qubit measurements on $T$-copies of\n",

    "$\\rho$. Each measurement is chosen randomly among the Pauli bases $X$,\n",

    "$Y$, or $Z$ to yield random $n$ pure product states $|s_i\\rangle$ for\n",

    "each copy:\n",

    "\n",

    "$$|s_{i}^{(t)}\\rangle \\in \\{|0\\rangle, |1\\rangle, |+\\rangle, |-\\rangle, |i+\\rangle, |i-\\rangle\\}.$$\n",

    "\n",

    "$$S_T(\\rho) = \\big\\{|s_{i}^{(t)}\\rangle: i\\in\\{1,\\ldots, n\\},\\ t\\in\\{1,\\ldots, T\\} \\big\\}.$$\n",

    "\n",

    "Each of the $|s_i^{(t)}\\rangle$ provides us with a snapshot of the state\n",

    "$\\rho$, and the $nT$ measurements yield the complete set $S_{T}$, which\n",

    "requires just $3nT$ bits to be stored in classical memory. This is\n",

    "discussed in further detail in our previous demo about\n",

    "`classical shadows <tutorial_classical_shadows>`{.interpreted-text\n",

    "role=\"doc\"}.\n"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "voRlU_8blMUI"

   },

   "source": [

    "![](/demonstrations/ml_classical_shadows/class_shadow_prep.png){.align-center\n",

    "width=\"100.0%\"}\n"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "nJd4CxtClMUI"

   },

   "source": [

    "To prepare a classical shadow for the ground state of the Heisenberg\n",

    "model, we simply reuse the circuit template used above and reconstruct a\n",

    "`QNode` utilizing a device that performs single-shot measurements.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 53,

   "metadata": {

    "id": "puR1Cx8WlMUI"

   },

   "outputs": [],

   "source": [

    "dev_oshot = qml.device(\"lightning.qubit\", wires=num_qubits, shots=1)\n",

    "circuit_oshot = qml.QNode(circuit, dev_oshot)"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "729nluzslMUI"

   },

   "source": [

    "Now, we define a function to build the classical shadow for the quantum\n",

    "state prepared by a given $n$-qubit circuit using $T$-copies of\n",

    "randomized Pauli basis measurements\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 54,

   "metadata": {

    "colab": {

     "base_uri": "https://localhost:8080/",

     "height": 0

    },

    "id": "2lKMEDwRlMUI",

    "outputId": "73e9bbda-afdd-4cee-fce0-2c611ca063e7"

   },

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "First five measurement outcomes =\n",

      " [[-1.  1.  1. -1.]\n",

      " [ 1. -1. -1.  1.]\n",

      " [ 1. -1. -1.  1.]\n",

      " [ 1.  1.  1. -1.]\n",

      " [-1. -1. -1.  1.]]\n",

      "First five measurement bases =\n",

      " [[2 2 2 1]\n",

      " [0 1 0 0]\n",

      " [2 2 2 2]\n",

      " [1 2 2 2]\n",

      " [1 2 2 1]]\n"

     ]

    }

   ],

   "source": [

    "def gen_class_shadow(circ_template, circuit_params, num_shadows, num_qubits):\n",

    "    # prepare the complete set of available Pauli operators\n",

    "    unitary_ops = [qml.PauliX, qml.PauliY, qml.PauliZ]\n",

    "    # sample random Pauli measurements uniformly\n",

    "    unitary_ensmb = np.random.randint(0, 3, size=(num_shadows, num_qubits), dtype=int)\n",

    "\n",

    "    outcomes = np.zeros((num_shadows, num_qubits))\n",

    "    for ns in range(num_shadows):\n",

    "        # for each snapshot, extract the Pauli basis measurement to be performed\n",

    "        meas_obs = [unitary_ops[unitary_ensmb[ns, i]](i) for i in range(num_qubits)]\n",

    "        # perform single shot randomized Pauli measuremnt for each qubit\n",

    "        outcomes[ns, :] = circ_template(circuit_params, observables=meas_obs)\n",

    "\n",

    "    return outcomes, unitary_ensmb\n",

    "\n",

    "\n",

    "outcomes, basis = gen_class_shadow(circuit_oshot, psi0, 100, num_qubits)\n",

    "print(\"First five measurement outcomes =\\n\", outcomes[:5])\n",

    "print(\"First five measurement bases =\\n\", basis[:5])"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "aM9kSgoElMUI"

   },

   "source": [

    "Furthermore, $S_{T}$ can be used to construct an approximation of the\n",

    "underlying $n$-qubit state $\\rho$ by averaging over $\\sigma_t$:\n",

    "\n",

    "$$\\sigma_T(\\rho) = \\frac{1}{T} \\sum_{1}^{T} \\big(3|s_{1}^{(t)}\\rangle\\langle s_1^{(t)}| - \\mathbb{I}\\big)\\otimes \\ldots \\otimes \\big(3|s_{n}^{(t)}\\rangle\\langle s_n^{(t)}| - \\mathbb{I}\\big).$$\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 55,

   "metadata": {

    "id": "uhlXzt7IlMUI"

   },

   "outputs": [],

   "source": [

    "def snapshot_state(meas_list, obs_list):\n",

    "    # undo the rotations done for performing Pauli measurements in the specific basis\n",

    "    rotations = [\n",

    "        qml.matrix(qml.Hadamard(wires=0)), # X-basis\n",

    "        qml.matrix(qml.Hadamard(wires=0)) @ qml.matrix(qml.adjoint(qml.S(wires=0))), # Y-basis\n",

    "        qml.matrix(qml.Identity(wires=0)), # Z-basis\n",

    "    ]\n",

    "\n",

    "    # reconstruct snapshot from local Pauli measurements\n",

    "    rho_snapshot = [1]\n",

    "    for meas_out, basis in zip(meas_list, obs_list):\n",

    "        # preparing state |s_i><s_i| using the post measurement outcome:\n",

    "        # |0><0| for 1 and |1><1| for -1\n",

    "        state = np.array([[1, 0], [0, 0]]) if meas_out == 1 else np.array([[0, 0], [0, 1]])\n",

    "        local_rho = 3 * (rotations[basis].conj().T @ state @ rotations[basis]) - np.eye(2)\n",

    "        rho_snapshot = np.kron(rho_snapshot, local_rho)\n",

    "\n",

    "    return rho_snapshot\n",

    "\n",

    "def shadow_state_reconst(shadow):\n",

    "    num_snapshots, num_qubits = shadow[0].shape\n",

    "    meas_lists, obs_lists = shadow\n",

    "\n",

    "    # Reconstruct the quantum state from its classical shadow\n",

    "    shadow_rho = np.zeros((2 ** num_qubits, 2 ** num_qubits), dtype=complex)\n",

    "    for i in range(num_snapshots):\n",

    "        shadow_rho += snapshot_state(meas_lists[i], obs_lists[i])\n",

    "\n",

    "    return shadow_rho / num_snapshots"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "wB7CZQJXlMUI"

   },

   "source": [

    "To see how well the reconstruction works for different values of $T$, we\n",

    "look at the\n",

    "[fidelity](https://en.wikipedia.org/wiki/Fidelity_of_quantum_states) of\n",

    "the actual quantum state with respect to the reconstructed quantum state\n",

    "from the classical shadow with $T$ copies. On average, as the number of\n",

    "copies $T$ is increased, the reconstruction becomes more effective with\n",

    "average higher fidelity values (orange) and lower variance (blue).\n",

    "Eventually, in the limit $T\\rightarrow\\infty$, the reconstruction will\n",

    "be exact.\n",

    "\n",

    "![Fidelity of the reconstructed ground state with different shadow sizes\n",

    "$T$](/demonstrations/ml_classical_shadows/fidel_snapshot.png){.align-center\n",

    "width=\"80.0%\"}\n"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "KuzuhucBlMUI"

   },

   "source": [

    "The reconstructed quantum state $\\sigma_T$ can also be used to evaluate\n",

    "expectation values $\\text{Tr}(O\\sigma_T)$ for some localized observable\n",

    "$O = \\bigotimes_{i}^{n} P_i$, where $P_i \\in \\{I, X, Y, Z\\}$. However,\n",

    "as shown above, $\\sigma_T$ would be only an approximation of $\\rho$ for\n",

    "finite values of $T$. Therefore, to estimate $\\langle O \\rangle$\n",

    "robustly, we use the median of means estimation. For this purpose, we\n",

    "split up the $T$ shadows into $K$ equally-sized groups and evaluate the\n",

    "median of the mean value of $\\langle O \\rangle$ for each of these\n",

    "groups.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 56,

   "metadata": {

    "id": "29rEs_S1lMUI"

   },

   "outputs": [],

   "source": [

    "def estimate_shadow_obs(shadow, observable, k=10):\n",

    "    shadow_size = shadow[0].shape[0]\n",

    "\n",

    "    # convert Pennylane observables to indices\n",

    "    map_name_to_int = {\"PauliX\": 0, \"PauliY\": 1, \"PauliZ\": 2}\n",

    "    if isinstance(observable, (qml.PauliX, qml.PauliY, qml.PauliZ)):\n",

    "        target_obs = np.array([map_name_to_int[observable.name]])\n",

    "        target_locs = np.array([observable.wires[0]])\n",

    "    else:\n",

    "        target_obs = np.array([map_name_to_int[o.name] for o in observable.obs])\n",

    "        target_locs = np.array([o.wires[0] for o in observable.obs])\n",

    "\n",

    "    # perform median of means to return the result\n",

    "    means = []\n",

    "    meas_list, obs_lists = shadow\n",

    "    for i in range(0, shadow_size, shadow_size // k):\n",

    "        meas_list_k, obs_lists_k = (\n",

    "            meas_list[i : i + shadow_size // k],\n",

    "            obs_lists[i : i + shadow_size // k],\n",

    "        )\n",

    "        indices = np.all(obs_lists_k[:, target_locs] == target_obs, axis=1)\n",

    "        if sum(indices):\n",

    "            means.append(\n",

    "                np.sum(np.prod(meas_list_k[indices][:, target_locs], axis=1)) / sum(indices)\n",

    "            )\n",

    "        else:\n",

    "            means.append(0)\n",

    "\n",

    "    return np.median(means)"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "MJrMeks1lMUI"

   },

   "source": [

    "Now we estimate the correlation matrix $C^{\\prime}$ from the classical\n",

    "shadow approximation of the ground state.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 57,

   "metadata": {

    "id": "rjSaQJwYlMUI"

   },

   "outputs": [],

   "source": [

    "coups = list(it.product(range(num_qubits), repeat=2))\n",

    "corrs = [corr_function(i, j) for i, j in coups]\n",

    "qbobs = [qob for qobs in corrs for qob in qobs]\n",

    "\n",

    "def build_estim_corrmat(coups, corrs, num_obs, shadow):\n",

    "    k = int(2 * np.log(2 * num_obs)) # group size\n",

    "    corr_mat_estim = np.zeros((num_qubits, num_qubits))\n",

    "    for idx, (i, j) in enumerate(coups):\n",

    "        corr = corrs[idx]\n",

    "        if i == j:\n",

    "            corr_mat_estim[i][j] = 1.0\n",

    "        else:\n",

    "            corr_mat_estim[i][j] = (\n",

    "                np.sum(np.array([estimate_shadow_obs(shadow, o, k=k+1) for o in corr])) / 3\n",

    "            )\n",

    "            corr_mat_estim[j][i] = corr_mat_estim[i][j]\n",

    "    return corr_mat_estim\n",

    "\n",

    "shadow = gen_class_shadow(circuit_oshot, psi0, 1000, num_qubits)\n",

    "expval_estmt = build_estim_corrmat(coups, corrs, len(qbobs), shadow)"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "AtpcKs2dlMUJ"

   },

   "source": [

    "This time, let us visualize the deviation observed between the exact\n",

    "correlation matrix ($C$) and the estimated correlation matrix\n",

    "($C^{\\prime}$) to assess the effectiveness of classical shadow\n",

    "formalism.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 58,

   "metadata": {

    "colab": {

     "base_uri": "https://localhost:8080/",

     "height": 371

    },

    "id": "MAviSXyklMUJ",

    "outputId": "9187eaa8-7340-49b6-91ac-871d95d1c191"

   },

   "outputs": [

    {

     "data": {

      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAFiCAYAAADC2W5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABBrUlEQVR4nO3de3xMZ/4H8M+ZXCaRy0TcEkSCRttVl6paXSGCJor+4roulYYuiq26tbv40bAtoUX1YlsaK6k2Lb+l9LIkSi5EqSJht1IREpGUiktGornNPL8//GZ+ppPLmZjMjJPP+/U6r1fznGfO+c5I55vncp5HEkIIEBERNTKVvQMgIqKmgQmHiIhsggmHiIhsggmHiIhsggmHiIhsggmHiIhsggmHiIhsggmHiIhsggmHiIhsggnnASNJEgYOHGjvMKwqNTUVkiRh+fLl9g5Flvj4eEiShPj4eHuHIsvy5cshSRJSU1PtHQo1cYpMOHl5eZAkqc4jKCjI3mE2KQ9SojT8/kyZMsXeocjyoMVLTZezvQNoTJ07d8bkyZNrPOfj42PbYKzk7NmzaNasmb3DsKo+ffrg7NmzaNmypb1DkWXUqFHo27cv/P397R0K0QNF0QnnoYceemC6aeR65JFH7B2C1TVr1uyBel8ajQYajcbeYRA9cBTZpdYQhi6fwsJCPP/88/Dz84NKpUJqaqrJGMORI0cQHh4OHx8fSJJkfH1ZWRliYmLwyCOPwM3NDb6+vhg+fDgyMjLM7nVvn3p8fDx69eqFZs2ayepyqqlrasqUKZAkCRcvXsS7776LRx55BGq1GoGBgVixYgX0er1Fn8Xp06cxYcIE+Pv7w9XVFYGBgZgzZw6uX79uVjclJQXPPPMM2rZtC7VajTZt2qB///7YvHkzgP8fnwGAtLQ0k25NwxhIbWM4QUFBCAoKQklJCWbNmgV/f394eHhgwIABOHnyJACgqKgIkydPRuvWreHu7o7w8HDk5OSYxfnFF19g4sSJeOihh9CsWTNoNBr0798fO3fuNKkXHx+Pjh07AgASEhJM4jWMgdQ2hmP4t7l69Sqio6PRsmVLuLu7o2/fvrWOn5w+fRrDhg2Dl5cXNBoNhg0bhn//+9/Gf9O8vLza/plkx3uvxMRE9OzZE+7u7vD398fcuXPx66+/1njt9PR0PPvss2jZsiXUajWCg4OxdOlS3Llzp86YiGqj6BaOpa5fv46nnnoKvr6+mDBhAsrLy+Ht7Q2tVgsAOHLkCFatWoWwsDDMmDEDly5dAgCUl5dj0KBB+P7779GrVy/MmzcPV69exfbt25GUlITPPvsM48aNM7vfW2+9hZSUFERGRiI8PBxOTk73Ff+rr76KtLQ0jBgxAhEREdi9ezeWL1+OyspKrFy5UtY1vvzyS/zxj3+ESqVCZGQkAgIC8OOPP+L9999HUlISjh07hubNmwMAvvnmGzz77LPw8fFBZGQk/P39ce3aNWRlZWHbtm2YMWMGgoKCEBMTgxUrViAwMNBknKFnz571xlNZWYmnn34a5eXlGD9+PK5evYodO3ZgyJAhOHLkCCIiIuDv74/Jkyfj/Pnz+OqrrzB8+HCcPXvW5PNcvHgxXF1dERISYozzyy+/xNixY/Huu+9izpw5xpjmzp2Ld955Bz169MDIkSON15Az7nfr1i2EhIRAo9EgKioKv/zyC7Zv346IiAicOHECjz32mLFuVlYW+vfvj7KyMowePRrBwcH44YcfEBISgh49etR7L0vjff/997Fv3z5ERkZi0KBB2LdvH959910UFxfj008/Nan7wQcf4M9//jN8fHzw7LPPonXr1vjhhx+wcuVKpKSkICUlBa6urrJiJDISCnTx4kUBQHTu3FnExMTUeOzdu9fkNQAEADF16lRRXV1tci4lJcV4/h//+IfZ/VasWCEAiOeee07o9Xpj+cmTJ4Wrq6vw8fERWq3WWB4TEyMACA8PD3H69GmL3hsAERoaalIWHR0tAIiOHTuKoqIiY/m1a9eEj4+P8PLyEhUVFfVeu7i4WHh7e4t27dqJvLw8k3OfffaZACBeeuklY9no0aMFAJGZmVnjteqL28Dw+cbExJiUBwYGCgBi3Lhxoqqqyli+Zs0aAUD4+PiI+fPnm3zms2bNEgDEzp07Ta6Vm5trdt/bt2+Lbt26CY1GI8rKyozlht+f6OjoGuPdunWrACC2bt1q9h4BiNmzZwudTmcsj4uLEwDEiy++aFI/JCREABCffvqpSfmyZcuM17p48WKNMdyrvngNv28ajUZkZ2cby+/cuSO6dOkiVCqVKCwsNJb/5z//Ec7OzqJHjx5m/46xsbECgFi7dm29cRH9lqITTl3H3LlzTV4DQLi6uopr166ZXc/whdirV68a79epUyfh4uIiCgoKzM5Nnz5dABAff/yxsczwBTB//nyL31tdCaemZGg4JyexrV+/3izWe/Xq1Uu0bNnS+LMh4fz0008NitugvoSTn59vUn7p0iUBQHh6epokCiGESE9PFwDEa6+9Vm9MQgixbt06AUCkpqYay+4n4Xh4eIjbt2+blFdVVQlnZ2eT35+8vDwBQPTo0cPs+qWlpaJ58+ZWTzg1fSaGc19++aWx7OWXXxYARHp6ull9nU4nWrVqJZ544ol64yL6LUV3qUVERGDfvn2y63fs2LHOmVJPPvmkWZlWq8WFCxfw6KOPon379mbnw8LC8NFHHyEzMxNRUVEm5/r06SM7NjmeeOIJszJDTLdu3ar39UePHgUAHDt2DLm5uWbny8vLUVxcjOLiYrRs2RITJkzArl270LdvX0yaNAmDBw9G//79rTrbrHnz5ujQoYNJmWF2WHBwsNmMPcO5oqIik/JffvkFq1evxt69e5Gfn282bvHb+g3VpUsXeHp6mpQ5OzujTZs2Jv8GWVlZAIB+/fqZXcPDwwM9e/ZESkqKVWIykPv7Yfg9SEpKwoEDB8xe4+LiguzsbKvGRk2DohOOpdq0aWPxecP4Tm2vNXwBGupZcj9LeXt7m5U5O9/9J9bpdPW+/saNGwCAjRs31lmvrKwMLVu2xLhx47B7926sX78eH374ITZu3AhJkhAWFoZ169bJGqOpT13vqa5zVVVVxrIbN27gySefxKVLl9CvXz8MGTIEPj4+cHJyQmZmJvbs2YOKior7jrW2mAxx3ftvYPh9aN26dY31rf27UVtsNf1+GH4P5I77EcnFhHOPe2edyT1v+J/46tWrNb7mypUrJvUsuZ+tGWI8c+aMyeB2XSIjIxEZGYnbt28jIyMDu3btwpYtWzB06FBkZ2c7xPNOW7ZswaVLl/D6669j6dKlJudWr16NPXv22Dwmw2f9yy+/1Hi+tt8nWzDEptVq4eXlZbc4SHk4Lfo+eXt7o1OnTjh//jwKCwvNzhumplrjr/3G9vvf/x4A8N1331n8Wi8vLwwdOhSbN2/GlClTcPXqVRw7dsx4XqVSyWplNQZD92BkZKTZuUOHDpmVGWa3NWa8hlloR44cMTt3584dY5ebHNaO1/B7YOhaI7IWJhwriI6ORlVVFRYvXgwhhLH89OnTiI+Ph0ajMZmu6qimTp0KLy8v/Pd//zf+85//mJ2/c+eOyZdQenp6jV9yhr/a3dzcjGW+vr64fPlyI0Rdv8DAQADA4cOHTcoTExPxr3/9y6x+8+bNIUkSCgoKGjWmfv36ITMzE9u3bzc599Zbbxm7teSwdryzZ8+Gs7Mz5syZY5z6f69bt27h1KlTVrkXNS2K7lI7f/58nSsNLFq0yORLsaH+8pe/4JtvvsG2bdtw9uxZDB482Pj8RXV1NT766KMHomuiVatWxmeGevTogaFDh+KRRx5BRUUF8vLykJaWhj/84Q/GiRgvv/wyioqKEBISgqCgIEiShMOHD+P7779H3759ERISYrz2oEGDsGPHDowcORKPP/44nJyc8F//9V/o3r17o7+vqKgorFmzBnPmzEFKSgoCAwORlZWFAwcOYPTo0di1a5dJfU9PTzz55JNIT09HVFQUgoODoVKpEBUVZUxe1vDee+9hwIABeO6557Bz50489NBDOHnyJI4ePYoBAwYgPT0dKlX9fxNaO97HHnsMf//73zFr1iw8/PDDGDZsGDp37ozbt2/jwoULSEtLw5QpU/Dhhx825G1TE6bohJObm4sVK1bUen7evHlWSThubm44ePAg1qxZg+3bt+Ptt99Gs2bNEBoaiiVLlph88Tq64cOH49SpU3jrrbfw7bffYv/+/fDw8ED79u0xdepUk7XpFi9ejF27duHEiRNISkqCi4sLgoKCsGbNGsyePdvkwct33nkHAHDw4EF89dVX0Ov1aN++vU0STvv27ZGWloa//OUv+Pbbb1FdXY1evXohOTkZBQUFZgkHALZt24b58+fj66+/RklJCYQQCAkJsWrCefzxx3Ho0CEsWrQIe/fuhSRJCAkJweHDh7F48WIAtU9CaOx4p0+fjp49e2L9+vVIT0/HV199BY1Ggw4dOmD+/PmIjo5u0HWpaZPEvX1ARGR3Op0OnTt3xq+//mrXyQNE1sYxHCI7qa6uRnFxsVn56tWrkZ+f/0CM+xFZgi0cIju5desW2rRpg6effhpdunRBVVUVjh07huPHj8Pf3x8nTpzgFgikKEw4RHZSWVmJefPm4eDBgygqKkJ5eTn8/f3xzDPPYNmyZWjXrp29QySyKiYcIiKyCY7hEBGRTTDhEBGRTTDhNGGGHTUf9Hs0BbXtMkr0IGHCofsidyvkpq62bbSJmhJFrzRA9lfTfipkuVGjRqFv376cJk0PNCYcalSdO3e2dwiKoNFooNFo7B0G0X1hl1oDpaenY+TIkWjTpg3UajUCAgIwevRosxWJy8rKEBMTg0ceeQRubm7w9fXF8OHDkZGRYXbN5cuXQ5IkpKamIj4+Hr169UKzZs0wcOBAAMDAgQMhSRLKy8uxdOlSdO7cGS4uLibdNBcvXsS0adPQoUMHqNVq+Pv7Y8qUKcjPz5f1voqKihATE4O+ffuidevWUKvVCAoKwuzZs832bgkKCkJCQgKAu7ulSpIESZKM8Rrq1DSG09DPJTExET179oS7uzv8/f0xd+5cs90762KIr7CwEJMmTULLli3h5eWF4cOH48KFCwCAs2fPYuTIkfD19YWXlxfGjh1b4xIz//jHPxAZGYmgoCDje4iIiDDbqXP58uUICwsDAKxYscL4Od3bFWnomrxw4QLWrVuH3/3ud1Cr1ZgyZQqAmsdwVq9eDUmSMHPmTLPYDOdmzZol+7Mhamxs4TTAO++8g/nz58Pd3R2jRo1Chw4dUFhYiMOHD+Of//yncbHO8vJyDBo0CN9//z169eqFefPm4erVq9i+fTuSkpKMKzP/1ltvvYWUlBRERkYiPDzcZBFMABgzZgyysrIwdOhQ+Pj4oGPHjgDubg0dERGBsrIyjBgxAsHBwcjLy8Onn36KvXv34rvvvkOnTp3qfG/p6elYt24dBg8ejN///vdwcXHBqVOn8MEHHyApKQknT540/qU9b948xMfHIysrC3PnzjVutlbfJIGGfi7vv/8+9u3bh8jISAwaNAj79u3Du+++i+LiYnz66ad13vNeN2/eREhICPz8/BAdHY1z587h66+/RnZ2Nvbs2YP+/fvjiSeewAsvvIATJ05g586duHHjBg4ePGhynT//+c/o0aMHhgwZglatWqGwsBC7d+/GkCFDsGvXLuP+OwMHDkReXh4SEhIQGhpqkpB/u0HdnDlzcPToUQwfPhzPPvtsrTuCAndXKd+/fz82bdqEoUOHGpfC+f777/Haa6/hd7/7HdavXy/7cyFqdIIskpmZKVQqlWjbtq24ePGiyTm9Xi8KCwuNP69YsUIAEM8995zQ6/XG8pMnTwpXV1fh4+MjtFqtsTwmJkYAEB4eHuL06dNm9w4NDRUARM+ePcX169dNzlVWVoqgoCDh5eUlTp48aXLu0KFDwsnJSYwYMcKkPDAwUAQGBpqUXb16Vdy+fdvs3gkJCQKAeOONN0zKo6OjBQCzz6KuezT0c9FoNCI7O9tYfufOHdGlSxehUqlMPve6ABAAxPz5803KZ82aJQAIHx8fsWHDBmO5Xq8Xw4YNEwDEiRMnTF5z4cIFs+sXFRWJtm3biuDgYJPylJQUAUDExMTUGJfhc2zfvr3Iz883O79161YBQGzdutWk/PLly6JFixbC19dXXL58WWi1WtG5c2ehVqtFVlZWXR8Fkc2xS81CmzZtgl6vxxtvvGH2l7wkSWjbtq3x54SEBLi4uBi7Nwwef/xxREdH49atW9i9e7fZPWbMmIFu3brVGsOKFSvg6+trUvb1118jLy8Pr776Kh5//HGTcyEhIYiMjMS//vUvaLXaOt9f69at4enpaVYeFRUFb29vfPvtt3W+Xo6Gfi5z587Fww8/bPzZ3d0dEydOhF6vx4kTJ2Tf39PTE2+88YZJ2cSJEwEALVq0wMsvv2wslyQJEyZMAACzXTgNLct7+fv7Y8yYMcjJyZHdjXmvV199FR06dJBdv127dtiyZQtu3LiByZMnY/bs2cjNzcWbb75pk60fiCzBLjULff/99wCA8PDwOutptVpcuHABjz76KNq3b292PiwsDB999BEyMzMRFRVlcq5Pnz51Xrum84adOH/66acap95euXIFer0e586dQ+/eveu8/q5du7Bp0yacPHkSN2/eNNnVs6ioqM7X1ud+PpcnnnjCrL7hGrdu3ZIdQ3BwMJo1a2ZSZpj91b17d5MkeO+53773CxcuIDY2FgcPHkRhYSEqKipMzhcVFVm8H019//Y1iYyMxMyZM40bog0bNswkaRI5CiYcC5WUlECSpHqnpxpaEm3atKnxvOH1NbU4antNXecNWxLXN5ZRVlZW5/l169bhlVdeQatWrRAeHo727dvD3d0dALBhwwazL1VL3c/nUtNmZM7Od3+Fa9rqujZ1Xaeuc1VVVcay8+fPo0+fPtBqtQgLC8Ozzz4Lb29vqFQqpKamIi0trUGfVX3/9rUZNWqUMeG89NJLDboGUWNjwrGQj48PhBD4+eef61zN1/DFVdsGWleuXDGpd6/f/oUt57zhOl999RVGjBhR5+trU11djddffx3+/v7IzMw0GbAWQuDNN99s0HVrirMhn4sjefvtt3Hz5k1s27bNZBdUAJg5cybS0tIadN36/u1rcuvWLUyfPh0eHh7Q6XSYM2cOTp069UBsa05NC8dwLGTo8khOTq6znre3Nzp16oTz58+jsLDQ7HxqaioAoGfPnlaJ6/e//z0A4LvvvmvwNYqLi1FSUoKnnnrKbHbUDz/8UOP0Y8MMOrktDFt/Lo0lNzcXAIwz0QyEEDVO7bb0c7LEjBkzcOnSJbzzzjt46623kJubiz//+c9Wvw/R/WLCsdDMmTPh5OSEpUuXmg0KCyFM+vmjo6NRVVWFxYsXQ9yzC8Tp06cRHx8PjUZjtV0dIyMj0aFDB+Me9L9VVVVl9ozQb7Vu3Rru7u44efIk7ty5Yyy/efMm5syZU+NrDJMXCgoKZMdqy8+lsRjGZn77ma5evRr//ve/zeo35HOSY8uWLfif//kfjBs3Dn/605/w0ksvYcSIEdi2bRsSExOtei+i+8UuNQt169YNGzZswMsvv4yuXbti5MiRCAwMxJUrV5Ceno7hw4djw4YNAO4+J/HNN99g27ZtOHv2LAYPHoxffvkF27dvR3V1NT766COrdXuo1Wr885//xDPPPIPQ0FAMGjQI3bp1gyRJyM/Px6FDh9CiRQtkZ2fXeg2VSoXZs2dj3bp16NGjB5599llotVrs3bsXgYGBJjPwDAYNGoS1a9dixowZGDNmDDw8PBAYGGg24H8vW34ujWXmzJnYunUrxowZgz/+8Y9o0aIFjh49ipMnT2L48OH45ptvTOo/8sgjaNu2LT7//HOo1Wq0b98ekiRhzpw5DV5B4Ny5c5g7dy4CAgKwefNmY/k//vEPdO/eHbNmzcJTTz1V42w6Iruw55zsB1lKSooYMWKE8PX1Fa6urqJ9+/ZizJgxIiMjw6ReaWmpWLZsmejSpYvxGZNnnnlGHDp0yOyahudNUlJSaryn4Tmculy+fFnMnTtXBAcHC7VaLby9vcWjjz4qpk2bJg4cOGBSt6ZnZCorK8XKlSuNr+/QoYNYuHChuH37do31hRDizTffFMHBwcLFxUUAEKGhoXXew5qfS23Pp9Tmt/EZXLx4UQAQ0dHRZudqe4YmJSVF9OvXT3h5eQkfHx8xbNgwceLEiVrjPXr0qAgNDRVeXl7G54EMzy/V9zzTb99nRUWF6NWrl1CpVCItLc2sfnJyspAkSfTt21dUVVXV86kQ2QZ3/CQiIpvgGA4REdkEEw4REdkEEw4REdkEEw4REdkEEw4REdkEEw4REdmEXR/81Ov1KCoqgpeXV4PWkCIi5RBC4Pbt22jbti1UKv4trER2TThFRUUICAiwZwhE5GAKCgpq3LqCHnx2TTiG5UtycnIcfimTxqTjo7fQVlh/UcsHkYdL0/3L/vbt2+j6SBeH/y745JNPcOjQIZw4cQJnzpxBZWUltm7diilTplh0Hb1ej40bN2Lz5s04f/48PD09MWTIEKxcubLereAfVHZNOIZuNC8vL4dfjr4xMeEAggkHAODZhBOOgaN3rxsW7m3ZsiX8/f0btLMrALz44ouIi4tD165d8fLLL6OoqAg7duxAcnIyjh49iuDgYCtHbn/87SYiskBcXBzy8vJw7do1zJw5s0HXSElJQVxcHAYMGICTJ09izZo12LZtG3bv3o0bN24odhM9rhZNRGSBIUOG3Pc1PvroIwDA66+/DldXV2P5M888g4EDByI5ORmXLl1Chw4d7vtejoQtHCIiG0tNTYWHhwf69etndi4iIgIAGrxrrCNjC4eIFKW8vByVlZWy6wshzMaN1Go11Gq1tUMDAJSVleHnn3/GY489ZtwJ9l6GsZucnJxGub89MeEQkWKUl5fD3csXqDbfDr02np6eKC0tNSmLiYnB8uXLrRzdXSUlJQBQ68Z7hglUhnpKwoRDRIpRWVkJVP8Kl8cmAk4u9b9AV4XSf3+GgoICk5myjdW6aeqYcIhIcSQXN0hOrvXWE6q7XVre3t42ezTD0LKprQWj1WpN6ikJEw4RKY6kcoKkMh8fMSNk1LEyDw8P+Pv74+LFi9DpdGbjOIaxGz6HQ0T0AJAkJ2PSqfOQbJ9wACA0NBRlZWXIyMgwO5eUlAQAGDBggK3DanRMOESkOJKTCpKTk4yjcb8Ci4uLkZ2djeLiYpPyGTNmAACWLVtmMqNu7969SE1NRXh4OAIDAxs1NntglxoRKY5KZpeakNPt9htxcXE4fPgwAODMmTPGstTUVABASEgIpk2bBgB4//33sWLFCrNZb2FhYZg2bRri4uLQq1cvDB8+HD///DO2b98OX19fvPfeexbH9SBgwiEixZE9htOAhHP48GEkJCSYlGVkZJh0jxkSTl02bdqEbt26YfPmzXjnnXfg6emJUaNGYeXKlejcubPFcT0IJCGE3ZaO1Gq10Gg0uHLlChfvbOJKuHgngKa9eKdWq0WHdv4oKSlp8PeB4TvFq/8CSM71T20W1RW4fWj9fd2T5GMLh4gUR1KpIMnZxI0bvdkUEw4RKY7K2QWSs4zncKC3QTRkwIRDRIpjmBZdf0X7TItuqphwiEh5/m/ac32EngnHlphwiEhx5M5Sk9UKIqthwiEixWHCcUxMOESkOCqVE1SN9BwONRwTDhEpzt1p0XJaOJwWbUtMOESkOOxSc0xMOESkOEw4jqlB7cnjx49j2LBh8PHxgYeHB/r27YsdO3ZYOzYiogaRtTWB3PXWyGosbuGkpKQgIiICbm5umDBhAry8vLBz506MHz8eBQUFWLhwYWPESUQkm9wHP+21H05TZVHCqa6uxvTp06FSqZCeno6ePXsCAF577TX06dMHS5YswdixYxW5jwMRPTgkmQ9+yqlD1mNRl9rBgweRm5uLSZMmGZMNcHfv7SVLlqCystJs2W4iIlszzFKr/+AsNVuyqIVj2GAoPDzc7FxERAQAIC0t7f6jIiK6D5w04JgsSjg5OTkAgODgYLNzfn5+8PT0NNapSUVFBSoqKow/a7VaS25PRCQLE45jsqg9WVJSAuBuF1pNvL29jXVqEhsbC41GYzwCAgIsuT0RkSwqlST7INuxaQfm4sWLUVJSYjwKCgpseXsiaiIklST7INuxqEvN0LKprRWj1WrRvHnzWl+vVquhVte/7SsR0f2QJAmSVH8ykVOHrMeiFo5h7KamcZorV66gtLS0xvEdIiJbkmR2p7GFY1sWJZzQ0FAAQHJystm5pKQkkzpERPYiSTK71NjCsSmLEs7gwYPRqVMnJCYmIjMz01heUlKCVatWwdXVFc8//7y1YyQisgjHcByTRWM4zs7OiIuLQ0REBAYMGGCytE1+fj7Wrl2LoKCgRgqViEgelSRBJaP1ItjCsSmL11ILCwvD4cOHERMTg+3bt6OqqgrdunXDmjVrMH78+MaIkYjIInJbL2zh2FaDtifo06cP9u7da+1YiIisggnHMXE/HCJSHLkPdQomHJtiwiEixVE5SVA5yUg4MuqQ9TDhEJHiqFSQ2cKxQTBkxIRDRIrDMRzHxIRDRIpjePBTTj2yHSYcIlIcPofjmJhwiEh55K4iwC41m2LCISLF4RiOY2LCISLFkfscDjdgsy0mHCJSHO6H45iYcIhIcSTV3UNOPbIdJhwiUhx2qTkmJhwiUhxOGnBMTDhEpDgcw3FMTDhEpDjsUnNMTDhEpDhc2sYxMeEQkeI4qSQ4cT8ch8OEQ0SKo5KZcPRMODbFhENEiiO3hcOEY1tMOESkOEw4jokJh4gUhwnHMXFhByJSHGcV4KySZBwNu/7x48cxbNgw+Pj4wMPDA3379sWOHTtkvz4+Pt74rFBNR2pqasMCc3Bs4RCR4jRmCyclJQURERFwc3PDhAkT4OXlhZ07d2L8+PEoKCjAwoULZV8rMjISPXv2NCsPCgqyOK4HgUMkHJ24ezRVTmzVo/Wdy/YOwSFUNA+0dwh2IydByCV3lprOwntWV1dj+vTpUKlUSE9PNyaL1157DX369MGSJUswduxYBAbK+3ccOXIkpkyZYlEMDzJ2qRGR4jhJKjipZBwWLhd98OBB5ObmYtKkSSYtE41GgyVLlqCyshIJCQlWfjfK4RAtHCIia5LbpWZpq8owthIeHm52LiIiAgCQlpYm+3qnTp3C9evXUV1djaCgIAwZMgQtWrSwKKYHCRMOESmOpQlHq9WalKvVaqjVarP6OTk5AIDg4GCzc35+fvD09DTWkePdd981+dnd3R0xMTH461//KvsaDxJ2qRGR4rg6qWQfABAQEACNRmM8YmNja7xuSUkJgLtdaDXx9vY21qlLx44d8d577+HcuXO4c+cOLl++jI8//hi+vr5YtGgR3nvvvQa+c8fGFg4RKY7cSQOG1aILCgrg7e1tLK+pdWNNoaGhCA0NNf7crl07REVFoVevXujduzeWL1+OWbNmwdlZWV/RbOEQkeIYutTkHMDdlsm9R20Jx9Cyqa0Vo9Vqa239yNG1a1eEhITgxo0bOHv2bIOv46iYcIhIceQ99Hn3sIRh7KamcZorV66gtLS0xvEdS7Rs2RIAUFZWdl/XcURMOESkOJa2cOQydIMlJyebnUtKSjKp0xA6nQ4//PADAMh+ludBwoRDRIrTWAln8ODB6NSpExITE5GZmWksLykpwapVq+Dq6ornn3/eWP7zzz8jOzvbrAvuxIkTZtfW6XRYtGgRzp8/j7CwMPj7+1v2ph8AyhqRIiIC4CTJnBZt4Y6fzs7OiIuLQ0REBAYMGGCytE1+fj7Wrl1rsizN4sWLkZCQgK1bt5qsKNC7d290794d3bt3R7t27XDjxg2kpaXh3LlzaN++PeLi4iyK60HBhENEimPpLDVLhIWF4fDhw4iJicH27dtRVVWFbt26Yc2aNRg/frysayxcuBBHjx7F/v37cePGDbi6uuKhhx7C0qVLsWDBAjRv3tziuB4EkhDCbquYGWZ0FP58xWRKYlPDtdQA5xv59g7BITTltdS0Wi3a+fuhpKSkwd8Hhu+UpXtOwM3Ds9765WWleCPyifu6J8nHFg4RKU5jLW1D94cJh4gUx0klL5k4cdqUTTHhEJHisIXjmJhwiEhxmHAcExMOESlOY85So4ZjwiEixXGSJFnP2Fj6HA7dHyYcIlIclSRBJSOZyKlD1sOEQ0SK4wR5z7c5NXokdC8mHCJSHJVKkjU+wzEc22LCISLF4RiOY2LCISLF4RiOY2LCISLFUUnyxnDYo2ZbTDhEpDgcw3FMTDhEpDjsUnNMFi9d98knn+DFF19E7969oVarIUkS4uPjGyE0IqKGcZLkH2Q7Frdwli5divz8fLRs2RL+/v7Iz+c+JkTkWNjCcUwWt3Di4uKQl5eHa9euYebMmY0RExHRfTEs3innINuxuIUzZMiQxoiDiMhqXFQSXGQkEzl1yHo4aYCIFEeS2aUmsUvNpphwiEhxuB+OY7JpwqmoqEBFRYXxZ61Wa8vbE1EToYK8hzq5w7Rt2fTzjo2NhUajMR4BAQG2vD0RNRGGtdTkHGQ7Nk04ixcvRklJifEoKCiw5e2JqIkwTIuWc5Dt2LRLTa1WQ61W2/KWRNQEOanuHnLqke1w0gARKY5KkvdQJ+cM2BYTDhEpjkrm+Ay71GzL4oQTFxeHw4cPAwDOnDljLEtNTQUAhISEYNq0adaLkIjIQlzaxjFZnHAOHz6MhIQEk7KMjAxkZGQYf2bCISJ74hiOY7I44cTHx3N1aCJyaGzhOCaO4RCR4kjS3UNOPbIdJhwiUhwVJKggo4Ujow5ZDxMOESkOWziOiQmHiBTn7nM48uqR7TDhEJHisIXjmJhwiEhxOIbjmJhwiEh5ZLZwmG9siwmHiBSHYziOiQmHiBRHgrzGC/ONbTHhEJHicKUBx8SEQ0SKI0HmLLVGj4TuxYRDRIqjgrztjLl2p20x4RCR4kiSBElGE0dOHbIeJhwiUhxuT+CYmHCISHHYpeaYmHCISHHYpeaYmHCISHH44KdjYsIhIkViLnE87MIkIsUxtHDkHA1x/PhxDBs2DD4+PvDw8EDfvn2xY8cOi65RUVGBv/3tbwgODoabmxvatm2LGTNm4JdffmlYUA8AtnCISHEacwwnJSUFERERcHNzw4QJE+Dl5YWdO3di/PjxKCgowMKFC+u9hl6vR2RkJJKSktC3b1+MGTMGOTk5iIuLw4EDB3D06FG0atXK4tgcHVs4RKQ4jdXCqa6uxvTp06FSqZCeno7Nmzdj3bp1yMrKQpcuXbBkyRLk5+fXe52EhAQkJSVh4sSJOHLkCFavXo2dO3fi73//Oy5cuIClS5c28J07NiYcIlIcyYLDEgcPHkRubi4mTZqEnj17Gss1Gg2WLFmCyspKJCQk1Hudjz76CAAQGxtr0sp68cUX0alTJ3z66af49ddfLYzO8THhEJHiGBbvlHNYIjU1FQAQHh5udi4iIgIAkJaWVuc1ysvLcezYMTz88MMIDAw0OSdJEp5++mmUlZXhhx9+sCi2B4FDjOFoK3QQFTp7h2E3re9ctncIdlftG1h/pSZApxP2DsFudHrrvXdLt5jWarUm5Wq1Gmq12qx+Tk4OACA4ONjsnJ+fHzw9PY11apObmwu9Xl/jNe69dk5ODvr371/ve3iQsIVDRIojCSH7AICAgABoNBrjERsbW+N1S0pKANztQquJt7e3sU5t5Fzj3npK4hAtHCIiqxL6u4ecegAKCgqMX/QAamzd0P1jC4eIFEcSetkHcLdVce9RW8IxtEpqa31otdpaWy6WXOPeenK98MILkCQJLVq0QEVFhazXnDhxAnPmzEHXrl3RvHlzqNVqBAYGYuLEiUhKSrLo/nIw4RCR8hhaOHIOC9w7vvJbV65cQWlpaa1jMwadOnWCSqWqdaynrnGi2ty+fRs7duyARqPBjRs3sHv37jrrV1VVYebMmXjyySfx4YcfokOHDoiOjsbLL7+Mrl274osvvsDQoUPx/vvvy45BDiYcIlIeIeQfFggNDQUAJCcnm50ztAgMdWrj7u6OPn364KeffjJ7ZkcIgf3798PDwwO9e/eWHdf27dtRVlaGdevWwcPDA1u2bKm1rk6nw5gxY7Bp0yaEhYUhLy8Pe/fuxYYNG/DWW2/hX//6F86fP4+RI0ciKChIdgxyMOEQkfI0Ugtn8ODB6NSpExITE5GZmWksLykpwapVq+Dq6ornn3/eWP7zzz8jOzvbrPtsxowZAIDFixdD3JP0Nm3ahAsXLuC5556Du7u77Li2bNkCb29vTJo0CSNHjsSBAwdqfQB1+fLl+OqrrzBo0CAkJSWhXbt2ZnXat2+PXbt24emnn5YdgxxMOESkOHdnoMkZw7GshePs7Iy4uDjo9XoMGDAAM2bMwMKFC9GjRw+cO3cOq1atMmkVLF68GI8++ii++OILk+tER0cjIiICn332Gf7whz9g0aJFGDt2LGbPno2OHTvijTfekB3Tjz/+iKNHj2LMmDFwd3dHVFQU9Ho9tm7dalY3NzcXsbGx0Gg0+PTTT+HsXPu8MUmSrD55ggmHiJSnkVo4ABAWFobDhw+jX79+2L59Oz744AO0adMGn3/+uax11ABApVJhz549WL58Oa5du4a3334bGRkZ+NOf/oTvvvvOonXUDN1nUVFRAIAhQ4bAz88PW7duhV5v+v7Wr18PnU6HWbNmwc/PT/Y9rEUSwsIUb0WGGR1n8wrhdc+UxKaGD37ywU+Dyib84KdWq0WHdv4oKSkxmaJs6TU0Gg2u5Z2Dt7eXjPq30Sqoy33d056qqqrQrl07qNVqXLp0ybhMzoIFC/D2228jKSnJZFWEgIAAXL58GWfOnMFjjz0m+z5RUVHw9vbGxo0bAQCzZs3CnTt3ZC3jcy8+h0NEymPhczgPqj179uDatWv461//arImW1RUFN5++21s2bLFmHCuXbuGy5cvw8vLC127drXoPu+88w7c3NyMPxvGqyzFhENEyiP0gF75CcfQnTZ58mST8scffxy/+93vsGfPHty4cQO+vr4oLi4GALRq1cribRl8fX1Nfm7evHmD4uUYDhEpjqUPfj6ICgoKkJycjJ49e9bYPRYVFYWKigp88sknAAAPDw8AwNWrV2HJSMpPP/0ESZJw8+ZNAEBhYSEkScLly5YPBTDhEJHyNOKkAUcRHx8PvV5vnCzwW8899xwkSTK2ggICAtCmTRuUlZXh+PHjdV773skGWVlZCAgIMLZqsrKy4Ovri/bt21scMxMOESlPIz346SiEENi6dSucnJwwceLEGusEBAQgNDQUp0+fxg8//ABJkvDSSy8BuLvvTk0tlF9//RUffPCByeKlp0+fRo8ePYw/Z2VloXv37g2Km2M4RKQ8Cp80cPDgQVy8eBE+Pj5YtmxZrfUM67Jt2bIFvXv3xqJFi3Dy5El88cUXePjhhzF06FB07twZOp0O586dw6FDh1BSUoLExETjNX6bYJhwiIjuIel1kPTVsuo9iAzdZLdu3apzGRuDzz77DOvXr4e7uzt27dqF7du3Iz4+HocOHcKXX36J5s2bo127dhg/fjwiIyMxZMgQ42tPnz5t0m2XlZWFV155pUFxM+EQkfIovIWTmJho0gqx1Pjx4zF+/Ph66926dQuXLl0ytmjKy8uRk5Nj0sVmCSYcIlIehSccWzl9+jTc3d2NK1f/+9//BgCLn+MxYMIhIsWRO+X5QZ4WbQtZWVno2rUrnJycjD8/9NBDFi0sei/OUiMi5dHr5R9Uqzlz5phMoT5+/DiefPLJBl+PCYeIlEfh06Jtrby8HCdOnMCuXbsQERHR4Osw4RCR8jSBBz9t6e9//zsGDx6MMWPGYMKECQ2+DsdwiEhxOIZjXQsWLMCCBQvu+zpMOESkPJyl5pCYcIhIeYSQmXA4hmNLTDhEpDxCB8hZRUA8mCsNPKiYcIhIcYReDyFjyrOcOmQ9Fs1SKywsxIYNGxAeHo4OHTrA1dUVfn5+GDNmDI4dO9ZYMRIRWUavk3+QzVjUwnnvvfewZs0adO7cGeHh4WjVqhVycnKwe/du7N69G4mJibLW5yEialRykwkTjk1ZlHD69OmD1NRUhIaGmpQfOnQIgwcPxqxZszBy5Eio1WqrBklEZAmh00Ho6k8mcuqQ9VjUpTZ69GizZAMA/fv3R1hYGG7evIkzZ85YLTgiogbh0jYOyWqTBlxcXO5e0JnzEIjIzvR6mV1qTDi2ZJXscOnSJXz77bfw9/dHt27daq1XUVGBiooK48+G3eiIiKxJ6HUQMhKOnDpkPfe9llpVVRWioqJQUVGBNWvWGJexrklsbCw0Go3xCAgIuN/bExGZEzK707jSgE3dV8LR6/WYMmUK0tPTMX36dJNtSGuyePFilJSUGI+CgoL7uT0RUY0MLRw5B9lOg7vU9Ho9XnjhBSQmJmLy5Mn48MMP632NWq3mDDYianycFu2QGpRw9Ho9pk6dio8//hgTJ05EfHw8VCrudEBEDkLuDDROGrApixPOvclm/Pjx2LZtW53jNkREtsbncByTRQnH0I328ccfY9y4cfjkk0+YbIjI8XBatEOyKOH87W9/Q0JCAjw9PdGlSxe88cYbZnVGjhyJnj17Wis+IiLLcQzHIVmUcPLy8gAApaWlWLlyZY11goKCmHCIyK64WrRjsijhxMfHIz4+vpFCISKyEl0lUC2ju19X2fixkBHXoSEixWELxzEx4RCR8nDSgENiwiEi5eGkAYfEhENEisPncBwTEw4RKQ9XGnBITDhEpDzsUnNITDhEpDjcD8cxMeEQkeJwWrRjYsIhIsURegGhk5NwhA2iIQMmHCJSHKHTy0s4MuqQ9TDhEJHisEvNMTHhEJHisIXjmJhwiEhxmHAcExMOESmO0Omg50oDDkdl7wCIiKxNCL1xHKfOQ9i+haPVarFgwQIEBgZCrVYjKCgIr776KkpLSy26jiRJtR5TpkxpnODvE1s4RKQ4jtqlVlZWhtDQUGRmZiI8PBwTJ07EqVOnsHbtWqSlpSE9PR1ubm6yrxcYGFhjcnHUTTCZcIhIcRw14bz55pvIzMzEX//6V6xevdpYvmjRIqxZswZvv/02Fi9eLPt6QUFBWL58eSNE2jjYpUZEiiP0Ql6Xmg0f/BRCIC4uDp6enli2bJnJuWXLlsHT0xNxcXE2i8ce2MIhIsXR6/TQy2i9yKljLTk5OSgqKkJERAQ8PDxMznl4eKBfv35ISkpCQUEBAgICZF3z1q1b2Lx5M4qLi+Hr64t+/fqhW7dujRG+VTDhEJHiWNqlptVqTcrVajXUarVVY8rJyQEABAcH13g+ODgYSUlJyMnJkZ1wsrKy8OKLL5qUDR06FAkJCWjduvX9BdwI2KVGRIpjSDhyDgAICAiARqMxHrGxsVaPqaSkBACg0WhqPO/t7W1Srz4LFy7EkSNHUFxcDK1WiyNHjuCZZ57Bvn37MGLECOgccMo3WzhEpDiGadFy6gFAQUGB8QsfQJ2tm4ULF6KiokJ2LHPnzq21VXM/1q5da/LzU089ha+//hqDBg1CWloa9uzZg9GjR1v9vveDCYeIFMfSLjVvb2+ThFOXTZs2oaysTHYsY8eORXBwsLFlU1sLxtCtV1sLSA6VSoXp06cjLS0NGRkZTDg18XBRwdOl6fbuVTQPtHcIdqfTcZl4AHB1kuwdgt1Y87035rRoSx/QNDC0cgxjOb9V3xiPXC1btgQAi5KirThEwiEisiZdVTV0qvr/iNVVVdsgmruCg4PRtm1bZGRkoKyszGSmWllZGTIyMtCxY0fZEwZqc+zYMQB3n9FxNE23WUFEinW3haOTcdhuWrQkSZg2bRpKS0vx+uuvm5x7/fXXUVpaiunTp5uU37lzB9nZ2bh06ZJJ+ZkzZ1BVVWV2jyNHjmDNmjVwcXHBuHHjrP8m7hNbOESkOI66H85f/vIX7NmzB2vWrMGpU6fQq1cvnDx5EsnJyXjyyScxb948k/rff/89wsLCEBoaitTUVGP5unXr8M033yAkJAQBAQFwcXHBf/7zHyQnJ0OSJGzcuBGdO3e26XuTgwmHiBRH6GWO4dg44Xh4eCAtLQ3Lly/Hzp07kZKSAn9/fyxcuBAxMTFwd3eXdZ3IyEjcunULWVlZ2L9/PyorK+Hn54cJEyZg3rx56NOnTyO/k4aRhBB2G63VarXQaDS4VPiz7BkiSuSkaroDxQY67i0PoGlPGtBqtfDz80NJSUmDvw8M3yknX/gveLq61Fu/tLIKvf7x5X3dk+RjC4eIFMcRl7YhJhwiUiBHHcNp6phwiEhxHHV7gqaOCYeIFEfoBISMh4nl1CHrYcIhIsXR62WO4bBLzaaYcIhIce5uwCajhcPZkTbFhENEiqPXAXpV/clE73gr+CsaEw4RKY7Q6SFUnDTgaJhwiEhxhE5AyGjhcNKAbTHhEJHi6HVCZpcaE44tMeEQkeKwS80xMeEQkeLohYBexgw0vf2WkmySmHCISHl0AkKSkUzYpWZTTDhEpDh6nR56iYt3OhomHCJSHCGzhcNZarbFhENEisOE45iYcIhIcdil5piYcIhIcYSQuZYaZ6nZFBMOESmOXiegBx/8dDRMOESkOEInICDnwU8mHFtiwiEixbmbcDhpwNGoLH1BeXk5FixYgAEDBqBt27Zwc3ODn58f+vXrh61bt6Kqqqox4iQikk1XpZN9kO1YnHBKS0vxwQcfQJIkDB8+HAsWLMCoUaNQWFiIF154ASNGjOAuekRkV3qdkH2Q7Vjcpebr64uSkhK4urqalFdXV+Ppp59GcnIy9u7di+HDh1stSCIiSwi9zC417vhpUxa3cFQqlVmyAQBnZ2eMGjUKAHD+/Pn7j4yIqKF0+rsrRtdzgM/h2JTVJg3o9Xrs27cPAPDYY49Z67JERBbT64SslaDlrChN1tPghFNZWYlVq1ZBCIHr16/jwIEDyM7OxtSpUzF48OAaX1NRUYGKigrjz1qttqG3JyKqldAJWQ91skvNtu4r4axYscL4syRJeOWVVxAbG1vra2JjY01eQ0TUGPRCZguHKw3YlMVjOAaenp4QQkCn06GgoAAbN25EXFwcBg4cWGvLZfHixSgpKTEeBQUFDQ6ciKg2OiFkH2Q7DU44xguoVGjfvj1mzZqFzZs3IyMjAytXrqyxrlqthre3t8lBRGRtOiH/INux6koD4eHhAIDU1FRrXpaIyCJyWy9s4diWVRNOUVERAMDFxcWalyUisojc1gtbOLZlcZfajz/+iDt37piV37lzBwsWLAAADBs27P4jIyJqIL3M8RtOGrAti1s4O3bswPr16xESEoKgoCB4e3ujsLAQe/fuxfXr19G/f3/Mnz+/MWIlIpJFB5ktnEaPhO5lccIZMWIEioqKcOTIEXz33XcoLS2FRqNB9+7dMWHCBLzwwgtwduYi1ERkPzohoJOxtA3HcGzL4szQu3dv9O7duzFiISKyCp2Q13rhGI5tsSlCRIrDhOOYmHCISHHYpeaYmHCISHH0Mls4XErNtphwiEhx2MJxTEw4RKQ4HMNxTEw4RKQ4dxOOnBaODYIhIyYcIlIctnAcExMOESkOx3AcExMOESmOAKCXWY9shwmHiBSHLRzHdN8bsBERORpH3YAtMzMTS5YsQUREBFq1agVJkjBw4MAGX+/48eMYNmwYfHx84OHhgb59+2LHjh3WC9jK2MIhIsWp1AtIUv3ZpNLGLZzdu3cjNjYWrq6u6NKlC4qLixt8rZSUFERERMDNzQ0TJkyAl5cXdu7cifHjx6OgoAALFy60YuTWwRYOESmOnL1w5O4Kak3jxo3DiRMnUFpaiv379zf4OtXV1Zg+fTpUKhXS09OxefNmrFu3DllZWejSpQuWLFmC/Px8K0ZuHUw4RKQ4epndabZe2qZr167o1avXfe+KfPDgQeTm5mLSpEno2bOnsVyj0WDJkiWorKxEQkLCfUZrfexSIyLFUfqkgdTUVABAeHi42bmIiAgAQFpami1DkoUJh4gU51foZU0IqPy/ydNardakXK1WQ61WN0ZoVpGTkwMACA4ONjvn5+cHT09PYx1HwoRDRIrh6uoKPz8/fHqlUPZrPD09ERAQYFIWExOD5cuXWzk66ykpKQFwtwutJt7e3sY6joQJh4gUw83NDRcvXkRlZaXs1wghIEmSSVldrZuFCxeioqJC9vXnzp1bY0ukKWLCISJFcXNzg5ubW6Ndf9OmTSgrK5Ndf+zYsVZPOIaWTW2tGK1Wi+bNm1v1ntbAhENEZIHS0lJ7h2BMYDk5OXjiiSdMzl25cgWlpaXo06ePPUKrE6dFExE9YEJDQwEAycnJZueSkpJM6jgSJhwiIgdVVVWF7Oxs5ObmmpQPHjwYnTp1QmJiIjIzM43lJSUlWLVqFVxdXfH888/bONr62bVLTfzfHPjbt2/bMwy7c1JJ9VdSOB03lwcAuDo13d8Fw/eAeECfjZEjOzsbq1evBgD8+uuvxrIpU6YY68THxxv/u7CwEI8++igCAwORl5dnLHd2dkZcXBwiIiIwYMAAk6Vt8vPzsXbtWgQFBdngHVlI2FFBQYHA3RXCefDgwUMAEAUFBfb8WmpUKSkp9b7/e128eFEAEIGBgTVe79ixY2Lo0KHC29tbuLu7iz59+ojPP//cBu+kYSQh7PfnhF6vR1FREby8vMymJdqKVqtFQEAACgoK4O3tbZcY7I2fAT8DwP6fgRACt2/fRtu2baFSsbdfiezapaZSqdC+fXt7hmDk7e3dZL9oDPgZ8DMA7PsZ1PYgIykD/4wgIiKbYMIhIiKbaPIJR61WIyYmxqEX6mts/Az4GQD8DKjx2XXSABERNR1NvoVDRES2wYRDREQ2wYRDREQ2wYRDREQ20WQTzvHjxzFs2DD4+PjAw8MDffv2xY4dO+wdls188sknePHFF9G7d2+o1WpIkmSyhpPSFRYWYsOGDQgPD0eHDh2MO0WOGTMGx44ds3d4NlFeXo4FCxZgwIABaNu2Ldzc3ODn54d+/fph69atqKqqsneIpDBNcpZaSkoKIiIi4ObmVuOidwsXLrR3iI0uKCgI+fn5aNmyJTw8PJCfn4+tW7eaLCKoZIsWLcKaNWvQuXNnDBw4EK1atUJOTg52794NIQQSExMxfvx4e4fZqIqLixEQEIA+ffqgS5cuaNWqFW7evIm9e/ciPz8f4eHh2Lt3L5eZIeux3zJu9lFVVSU6d+4s1Gq1OHXqlLH81q1bokuXLsLV1VXk5eXZL0Ab2b9/v/F9xsbGCgBi69at9g3Khnbu3ClSU1PNytPT04WLi4to3ry5KC8vt0NktqPT6URFRYVZeVVVlRg4cKAAIL7++ms7REZK1eT+dDl48CByc3MxadIk9OzZ01iu0WiwZMkSVFZWIiEhwX4B2siQIUMQGBho7zDsZvTo0TVuUNW/f3+EhYXh5s2bOHPmjB0isx2VSgVXV1ezcmdnZ4waNQoAcP78eVuHRQrW5BJOamoqACA8PNzsXEREBAAgLS3NliGRg3FxcQFw94u3KdLr9di3bx8A4LHHHrNzNKQkTe7/qJycHAD/vyf4vfz8/ODp6WmsQ03PpUuX8O2338Lf3x/dunWzdzg2UVlZiVWrVkEIgevXr+PAgQPIzs7G1KlTMXjwYHuHRwrS5BJOSUkJgNqXQff29jbWoaalqqoKUVFRqKiowJo1a+Dk5GTvkGyisrISK1asMP4sSRJeeeUVxMbG2jEqUqIm16VGVBO9Xo8pU6YgPT0d06dPR1RUlL1DshlPT08IIaDT6VBQUICNGzciLi4OAwcOhFartXd4pCBNLuEYWja1tWK0Wi03gWpi9Ho9XnjhBSQmJmLy5Mn48MMP7R2SXRg2RJw1axY2b96MjIwMrFy50t5hkYI0uYRjGLupaZzmypUrKC0trXF8h5RJr9dj6tSpSEhIwMSJExEfH8/nTvD/k2oMk2yIrKHJ/Z9lmAqbnJxsdi4pKcmkDimbIdl8/PHHGD9+PLZt29Zkxm3qU1RUBOD/Z+wRWUOTSziDBw9Gp06dkJiYiMzMTGN5SUkJVq1aBVdXVzz//PP2C5BswtCN9vHHH2PcuHH45JNPmlyy+fHHH3Hnzh2z8jt37mDBggUAgGHDhtk6LFIwLm3TRJe2iYuLw+HDhwEAZ86cwcmTJ9GvXz889NBDAICQkBBMmzbNniE2quXLl2PFihXw9PTE3Llza3zmZuTIkSYPByvN8uXLsX79eoSEhCAoKAje3t4oLCzE3r17cf36dfTv3x9JSUlwd3e3d6ikFPZd6MB+jh07JoYOHSq8vb2Fu7u76NOnj/j888/tHZbNREdHCwC1HtHR0fYOsVHV9/7RBJb6OX78uJg+fbro2rWr8PHxEc7OzqJFixYiLCxMbNq0SVRVVdk7RFKYJtnCISIi22tyYzhERGQfTDhERGQTTDhERGQTTDhERGQTTDhERGQTTDhERGQTTDhERGQTTDhERGQTTDhERGQTTDhERGQTTDhERGQTTDhERGQTTDhERGQT/wtiTvJLVTVnAQAAAABJRU5ErkJggg==\n",

      "text/plain": [

       "<Figure size 420x400 with 2 Axes>"

      ]

     },

     "metadata": {},

     "output_type": "display_data"

    }

   ],

   "source": [

    "fig, ax = plt.subplots(1, 1, figsize=(4.2, 4))\n",

    "im = ax.imshow(expval_exact-expval_estmt, cmap=plt.get_cmap(\"RdBu\"), vmin=-1, vmax=1)\n",

    "ax.xaxis.set_ticks(range(num_qubits))\n",

    "ax.yaxis.set_ticks(range(num_qubits))\n",

    "ax.xaxis.set_tick_params(labelsize=14)\n",

    "ax.yaxis.set_tick_params(labelsize=14)\n",

    "ax.set_title(\"Error in estimating the\\ncorrelation matrix\", fontsize=14)\n",

    "\n",

    "bar = fig.colorbar(im, pad=0.05, shrink=0.80)\n",

    "bar.set_label(r\"$\\Delta C_{ij}$\", fontsize=14, rotation=0)\n",

    "bar.ax.tick_params(labelsize=14)\n",

    "plt.show()"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "CvaBmbD-lMUJ"

   },

   "source": [

    "Training Classical Machine Learning Models\n",

    "==========================================\n"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "P9L3yo0tlMUJ"

   },

   "source": [

    "There are multiple ways in which we can combine classical shadows and\n",

    "machine learning. This could include training a model to learn the\n",

    "classical representation of quantum systems based on some system\n",

    "parameter, estimating a property from such learned classical\n",

    "representations, or a combination of both. In our case, we consider the\n",

    "problem of using\n",

    "`kernel-based models <tutorial_kernel_based_training>`{.interpreted-text\n",

    "role=\"doc\"} to learn the ground-state representation of the Heisenberg\n",

    "model Hamiltonian $H(x_l)$ from the coupling vector $x_l$, where\n",

    "$x_l = [J_{i,j} \\text{ for } i < j]$. The goal is to predict the\n",

    "correlation functions $C_{ij}$:\n",

    "\n",

    "$$\\big\\{x_l \\rightarrow \\sigma_T(\\rho(x_l)) \\rightarrow \\text{Tr}(\\hat{C}_{ij} \\sigma_T(\\rho(x_l))) \\big\\}_{l=1}^{N}.$$\n",

    "\n",

    "Here, we consider the following kernel-based machine learning model:\n",

    "\n",

    "$$\\hat{\\sigma}_{N} (x) = \\sum_{l=1}^{N} \\kappa(x, x_l)\\sigma_T (x_l) = \\sum_{l=1}^{N} \\left(\\sum_{l^{\\prime}=1}^{N} k(x, x_{l^{\\prime}})(K+\\lambda I)^{-1}_{l, l^{\\prime}} \\sigma_T(x_l) \\right),$$\n",

    "\n",

    "where $\\lambda > 0$ is a regularization parameter in cases when $K$ is\n",

    "not invertible, $\\sigma_T(x_l)$ denotes the classical representation of\n",

    "the ground state $\\rho(x_l)$ of the Heisenberg model constructed using\n",

    "$T$ randomized Pauli measurements, and $K_{ij}=k(x_i, x_j)$ is the\n",

    "kernel matrix with $k(x, x^{\\prime})$ as the kernel function.\n",

    "\n",

    "Similarly, estimating an expectation value on the predicted ground state\n",

    "$\\sigma_T(x_l)$ using the trained model can then be done by evaluating:\n",

    "\n",

    "$$\\text{Tr}(\\hat{O} \\hat{\\sigma}_{N} (x)) = \\sum_{l=1}^{N} \\kappa(x, x_l)\\text{Tr}(O\\sigma_T (x_l)).$$\n",

    "\n",

    "We train the classical kernel-based models using $N = 70$ randomly\n",

    "chosen values of the coupling matrices $J$.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 59,

   "metadata": {

    "id": "awX0qjXKlMUJ"

   },

   "outputs": [],

   "source": [

    "# imports for ML methods and techniques\n",

    "from sklearn.model_selection import train_test_split, cross_val_score\n",

    "from sklearn import svm\n",

    "from sklearn.kernel_ridge import KernelRidge"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "bQWsJyqGlMUJ"

   },

   "source": [

    "First, to build the dataset, we use the function `build_dataset` that\n",

    "takes as input the size of the dataset (`num_points`), the topology of\n",

    "the lattice (`Nr` and `Nc`), and the number of randomized Pauli\n",

    "measurements ($T$) for the construction of classical shadows. The\n",

    "`X_data` is the set of coupling vectors that are defined as a stripped\n",

    "version of the coupling matrix $J$, where only non-duplicate and\n",

    "non-zero $J_{ij}$ are considered. The `y_exact` and `y_clean` are the\n",

    "set of correlation vectors, i.e., the flattened correlation matrix $C$,\n",

    "computed with respect to the ground-state obtained from exact\n",

    "diagonalization and classical shadow representation (with $T=500$),\n",

    "respectively.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 60,

   "metadata": {

    "colab": {

     "base_uri": "https://localhost:8080/",

     "height": 0

    },

    "id": "oNjK4cvrlMUJ",

    "outputId": "64d55c2a-9d38-4349-dbd6-04ab57db94fc"

   },

   "outputs": [

    {

     "data": {

      "text/plain": [

       "((100, 4), (100, 16), (100, 16))"

      ]

     },

     "execution_count": 60,

     "metadata": {},

     "output_type": "execute_result"

    }

   ],

   "source": [

    "def build_dataset(num_points, Nr, Nc, T=500):\n",

    "\n",

    "    num_qubits = Nr * Nc\n",

    "    X, y_exact, y_estim = [], [], []\n",

    "    coupling_mats = build_coupling_mats(num_points, Nr, Nc)\n",

    "\n",

    "    for coupling_mat in coupling_mats:\n",

    "        ham = Hamiltonian(coupling_mat)\n",

    "        eigvals, eigvecs = sp.sparse.linalg.eigs(ham.sparse_matrix())\n",

    "        psi = eigvecs[:, np.argmin(eigvals)]\n",

    "        shadow = gen_class_shadow(circuit_oshot, psi, T, num_qubits)\n",

    "\n",

    "        coups = list(it.product(range(num_qubits), repeat=2))\n",

    "        corrs = [corr_function(i, j) for i, j in coups]\n",

    "        qbobs = [x for sublist in corrs for x in sublist]\n",

    "\n",

    "        expval_exact = build_exact_corrmat(coups, corrs, circuit_exact, psi)\n",

    "        expval_estim = build_estim_corrmat(coups, corrs, len(qbobs), shadow)\n",

    "\n",

    "        coupling_vec = []\n",

    "        for coup in coupling_mat.reshape(1, -1)[0]:\n",

    "            if coup and coup not in coupling_vec:\n",

    "                coupling_vec.append(coup)\n",

    "        coupling_vec = np.array(coupling_vec) / np.linalg.norm(coupling_vec)\n",

    "\n",

    "        X.append(coupling_vec)\n",

    "        y_exact.append(expval_exact.reshape(1, -1)[0])\n",

    "        y_estim.append(expval_estim.reshape(1, -1)[0])\n",

    "\n",

    "    return np.array(X), np.array(y_exact), np.array(y_estim)\n",

    "\n",

    "X, y_exact, y_estim = build_dataset(100, Nr, Nc, 500)\n",

    "X_data, y_data = X, y_estim\n",

    "X_data.shape, y_data.shape, y_exact.shape"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "V18wEzs_lMUK"

   },

   "source": [

    "Now that our dataset is ready, we can shift our focus to the ML models.\n",

    "Here, we use two different Kernel functions: (i) Gaussian Kernel and\n",

    "(ii) Neural Tangent Kernel. For both of them, we consider the\n",

    "regularization parameter $\\lambda$ from the following set of values:\n",

    "\n",

    "$$\\lambda = \\left\\{ 0.0025, 0.0125, 0.025, 0.05, 0.125, 0.25, 0.5, 1.0, 5.0, 10.0 \\right\\}.$$\n",

    "\n",

    "Next, we define the kernel functions $k(x, x^{\\prime})$ for each of the\n",

    "mentioned kernels:\n"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "Seuq-zVflMUK"

   },

   "source": [

    "$$k(x, x^{\\prime}) = e^{-\\gamma||x - x^{\\prime}||^{2}_{2}}. \\tag{Gaussian Kernel}$$\n",

    "\n",

    "For the Gaussian kernel, the hyperparameter\n",

    "$\\gamma = N^{2}/\\sum_{i=1}^{N} \\sum_{j=1}^{N} ||x_i-x_j||^{2}_{2} > 0$\n",

    "is chosen to be the inverse of the average Euclidean distance $x_i$ and\n",

    "$x_j$. The kernel is implemented using the radial-basis function (rbf)\n",

    "kernel in the `sklearn` library.\n"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "38VnvOMGlMUK"

   },

   "source": [

    "$$k(x, x^{\\prime}) = k^{\\text{NTK}}(x, x^{\\prime}). \\tag{Neural Tangent Kernel}$$\n",

    "\n",

    "The neural tangent kernel $k^{\\text{NTK}}$ used here is equivalent to an\n",

    "infinite-width feed-forward neural network with four hidden layers and\n",

    "that uses the rectified linear unit (ReLU) as the activation function.\n",

    "This is implemented using the `neural_tangents` library.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 61,

   "metadata": {

    "colab": {

     "base_uri": "https://localhost:8080/",

     "height": 0

    },

    "id": "Ekn9atIZlMUK",

    "outputId": "5b0e44dc-a7aa-40a6-f945-2b56e0c57df4"

   },

   "outputs": [],

   "source": [

    "from neural_tangents import stax\n",

    "init_fn, apply_fn, kernel_fn = stax.serial(\n",

    "    stax.Dense(32),\n",

    "    stax.Relu(),\n",

    "    stax.Dense(32),\n",

    "    stax.Relu(),\n",

    "    stax.Dense(32),\n",

    "    stax.Relu(),\n",

    "    stax.Dense(32),\n",

    "    stax.Relu(),\n",

    "    stax.Dense(1),\n",

    ")\n",

    "kernel_NN = kernel_fn(X_data, X_data, \"ntk\")\n",

    "\n",

    "for i in range(len(kernel_NN)):\n",

    "    for j in range(len(kernel_NN)):\n",

    "        kernel_NN.at[i, j].set((kernel_NN[i][i] * kernel_NN[j][j]) ** 0.5)"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "GJbMkW7klMUK"

   },

   "source": [

    "For the above two defined kernel methods, we obtain the best learning\n",

    "model by performing hyperparameter tuning using cross-validation for the\n",

    "prediction task of each $C_{ij}$. For this purpose, we implement the\n",

    "function `fit_predict_data`, which takes input as the correlation\n",

    "function index `cij`, kernel matrix `kernel`, and internal kernel\n",

    "mapping `opt` required by the kernel-based regression models from the\n",

    "`sklearn` library.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 62,

   "metadata": {

    "id": "7qbJKRGglMUK"

   },

   "outputs": [],

   "source": [

    "from sklearn.metrics import mean_squared_error\n",

    "\n",

    "def fit_predict_data(cij, kernel, opt=\"linear\"):\n",

    "\n",

    "    # training data (estimated from measurement data)\n",

    "    y = np.array([y_estim[i][cij] for i in range(len(X_data))])\n",

    "    X_train, X_test, y_train, y_test = train_test_split(\n",

    "        kernel, y, test_size=0.3, random_state=24\n",

    "    )\n",

    "\n",

    "    # testing data (exact expectation values)\n",

    "    y_clean = np.array([y_exact[i][cij] for i in range(len(X_data))])\n",

    "    _, _, _, y_test_clean = train_test_split(kernel, y_clean, test_size=0.3, random_state=24)\n",

    "\n",

    "    # hyperparameter tuning with cross validation\n",

    "    models = [\n",

    "        # Epsilon-Support Vector Regression\n",

    "        (lambda Cx: svm.SVR(kernel=opt, C=Cx, epsilon=0.1)),\n",

    "        # Kernel-Ridge based Regression\n",

    "        (lambda Cx: KernelRidge(kernel=opt, alpha=1 / (2 * Cx))),\n",

    "    ]\n",

    "\n",

    "    # Regularization parameter\n",

    "    hyperparams = [0.0025, 0.0125, 0.025, 0.05, 0.125, 0.25, 0.5, 1.0, 5.0, 10.0]\n",

    "    best_pred, best_cv_score, best_test_score = None, np.inf, np.inf\n",

    "    for model in models:\n",

    "        for hyperparam in hyperparams:\n",

    "            cv_score = -np.mean(\n",

    "                cross_val_score(\n",

    "                    model(hyperparam), X_train, y_train, cv=5,\n",

    "                    scoring=\"neg_root_mean_squared_error\",\n",

    "                )\n",

    "            )\n",

    "            if best_cv_score > cv_score:\n",

    "                best_model = model(hyperparam).fit(X_train, y_train)\n",

    "                best_pred = best_model.predict(X_test)\n",

    "                best_cv_score = cv_score\n",

    "                best_test_score = mean_squared_error(\n",

    "                    best_model.predict(X_test).ravel(), y_test_clean.ravel(), squared=False\n",

    "                )\n",

    "\n",

    "    return (\n",

    "        best_pred, y_test_clean, np.round(best_cv_score, 5), np.round(best_test_score, 5)\n",

    "    )"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {

    "id": "IGPcW5sglMUK"

   },

   "source": [

    "We perform the fitting and prediction for each $C_{ij}$ and print the\n",

    "output in a tabular format.\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 63,

   "metadata": {

    "colab": {

     "base_uri": "https://localhost:8080/",

     "height": 0

    },

    "id": "vxZOUxk9lMUK",

    "outputId": "07ed6444-ca91-48c2-fe8c-64b9bba549c8"

   },

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "              Correlation                    Gaussian kernel              Neural Tangent kernel\n",

      "               \t C_00 \t|                           [-0.  0.]                          [-0.  0.]\n",

      "               \t C_01 \t|                   [0.08991 0.07814]                  [0.1165  0.08399]\n",

      "               \t C_02 \t|                   [0.10099 0.06564]                  [0.10664 0.08351]\n",