examples/measurement_reduction/classical_shadows.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classical Shadows\n",
    "\n",
    "The output of a quantum computer is a histogram of measurements, corresponding to the different outcomes observed, usually expressed as bitstrings. The cost and duration of a quantum experiment is roughly linear with the number of shots used to build such histograms, which also correlates with the accuracy of the results. This measurement overhead can be prohibitive for state tomography and computation of observables for quantum chemistry applications.\n",
    "\n",
    "The emergent method *classical shadows* ([Nat. Phys. 2020, 16, 1050–1057](https://arxiv.org/abs/2002.08953)) has been developed to mitigate the measurement overhead by offloading quantum tasks to the pre- and post-processing steps. This prediction protocol exhibits logarithmic scaling with the number of terms to measure, in order to compute the expectation value of a term within a given accuracy.\n",
    "\n",
    "In short, gates representing a unitary operation are appended to the state prepared by a quantum circuit. Upon measurement, the output bitstring (ex: 100101...) is stored in a conventional register. This is repeated $N$ times, where $N$ defines the shadow size and is equivalent to the number of shots run on the quantum hardware. The classical shadow is built from the bitstrings stored in the classical register, and can be used to calculate observables. Different flavours of classical shadow exist, as they each come with a specific protocol for choosing a set of unitaries. The figure below, from [Nat. Phys. 16, 1050–1057 (2020)](https://arxiv.org/abs/2002.08953), illustrates the randomized shadow protocol.\n",
    "\n",
    "![Classical Shadows overview](../img/classical_shadow_overview.png \"Classical Shadow\")\n",
    "\n",
    "The randomized, derandomized and adaptive shadow protocols using the single-qubit Pauli basis as a set of unitaries  are currently available in Tangelo. This introduction will shed light on how to leverage their use in your own investigation, and highlight the main differences between these approaches. At the end of this notebook, a comparison is made between energies predicted with those techniques and the one computed by using a Hamiltonian partitioning approach relying on qubitwise commutativity.\n",
    "\n",
    "This notebook assumes that you already have installed Tangelo in your Python environment, or have updated your Python path so that the imports can be resolved. If not, executing the cell below installs the minimal requirements for this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Installation of tangelo if not already installed.\n",
    "try:\n",
    "    import tangelo\n",
    "except ModuleNotFoundError:\n",
    "    !pip install git+https://github.com/goodchemistryco/Tangelo.git@develop --quiet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## System and Reference Values\n",
    "To demonstrate the classical shadow capabilities, we first define a molecular system composed of two hydrogen atoms in the 3-21G basis. After using the scBK qubit mapping, the Hamiltonian spans 6 qubits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of qubits to describe the system: 6\n"
     ]
    }
   ],
   "source": [
    "from tangelo.molecule_library import mol_H2_321g\n",
    "from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping\n",
    "from tangelo.toolboxes.operators import count_qubits\n",
    "\n",
    "# Encoding options.\n",
    "encoding = \"scBK\"\n",
    "up_then_down = True\n",
    "\n",
    "# Converting the fermionic Hamiltonian according to the encoding options.\n",
    "qubit_ham = fermion_to_qubit_mapping(mol_H2_321g.fermionic_hamiltonian, \n",
    "                                     encoding, \n",
    "                                     n_spinorbitals=mol_H2_321g.n_active_sos, \n",
    "                                     n_electrons=mol_H2_321g.n_active_electrons, \n",
    "                                     up_then_down=up_then_down)\n",
    "\n",
    "print(f\"Number of qubits to describe the system: {count_qubits(qubit_ham)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, the Qubit Coupled-Cluster (QCC) ansatz has been chosen ([J. Chem. Theory Comput. 2018, 14(12), 6317–6326](https://pubs.acs.org/doi/10.1021/acs.jctc.8b00932)). While the circuit depth rapidly increases with the size of the basis (number of spin orbitals), the QCC ansatz admits a low-depth quantum circuit compared to the widely used Unitary Coupled-Cluster Single and Double (UCCSD) ansatz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of gates in the circuit: 79\n"
     ]
    }
   ],
   "source": [
    "from tangelo.toolboxes.ansatz_generator.qcc import QCC\n",
    "\n",
    "ansatz = QCC(mol_H2_321g, mapping=encoding, up_then_down=up_then_down)\n",
    "ansatz.build_circuit()\n",
    "\n",
    "print(f\"Number of gates in the circuit: {ansatz.circuit.size}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Variational Quantum Eigensolver (VQE) has been introduced in another notebook as a hybrid quantum–classical algorithm for simulating quantum systems. In the next cell, we focus on VQE within the context of optimizing the QCC variational parameters to minimize the 3-21G dihydrogen energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy from statevector simulation: -1.1473 hartree\n"
     ]
    }
   ],
   "source": [
    "from tangelo.algorithms.variational import VQESolver\n",
    "\n",
    "vqe = VQESolver({\"qubit_hamiltonian\": qubit_ham, \"ansatz\": ansatz.circuit})\n",
    "vqe.build()\n",
    "energy_ref = vqe.simulate()\n",
    "print(f\"Energy from statevector simulation: {energy_ref:.4f} hartree\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can get the optimized circuit with the `vqe.optimal_circuit` attribute, and we also know that with the provided Hamiltonian, the ideal energy is -1.1475 hartree."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classical Shadows\n",
    "In this section, we will revisit the computation of the 3-21G dihydrogen energy using the QCC ansatz and the classical shadow protocols. Three flavours of classical shadows are currently available in Tangelo. Each flavour has a distinct class: `RandomizedClassicalShadow`, `DerandomizedClassicalShadow` and `AdaptiveClassicalShadow`. Common methods `build` and `simulate` are defined to respectively construct a set of unitaries and simulate results with a given backend. As a side note, a user who wants to compare noisy simulation can do so by providing the appropriate `noise_model` to the `Simulator` object defined in the next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tangelo.linq import get_backend, Circuit\n",
    "\n",
    "# Definition of the shots budget (shadow size) and backend.\n",
    "n_shots_budget = 1000\n",
    "backend = get_backend(\"cirq\", noise_model=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Randomized Single-Qubit Pauli Classical Shadows\n",
    "Randomized classical shadows are the orginal implementation of this algorithm. Random bases are constructed from the tensor product of random single-qubit Pauli gates ($X$, $Y$ or $Z$). The post-processing consists of inverting the samples with the quantum channel $\\mathcal{M}^{-1} = \\bigotimes_{i=1}^n \\mathcal{M}_i^{-1}$. Observables can then be predicted from the median of means. More details can be found in [Nat. Phys. 16, 1050–1057 (2020)](https://arxiv.org/abs/2002.08953)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy from randomized single-Pauli classical shadow: -1.0640 hartree\n",
      "Error vs statevector simulation: 0.083 hartree\n"
     ]
    }
   ],
   "source": [
    "from tangelo.toolboxes.measurements import RandomizedClassicalShadow\n",
    "\n",
    "random_cs = RandomizedClassicalShadow(vqe.optimal_circuit)\n",
    "random_cs.build(n_shots_budget)\n",
    "random_cs.simulate(backend)\n",
    "\n",
    "energy_random = random_cs.get_observable(qubit_ham)\n",
    "print(f\"Energy from randomized single-Pauli classical shadow: {energy_random:.4f} hartree\")\n",
    "print(f\"Error vs statevector simulation: {abs(energy_ref-energy_random):.3f} hartree\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Derandomized Single-Qubit Pauli Classical Shadows\n",
    "The derandomized procedure aims to fix a drawback of the randomized procedure. As the number of qubits increases, the probability of matching random tensor products of single-Pauli bases with terms in a qubit Hamiltonian decreases greatly. Typically, bigger shadows are needed to compute an observable within a reasonable accuracy. The derandomized procedure takes as input a qubit operator to steer the basis choices, thus removing the randomness. This results in bases hitting the relevant operator systematically, improving the accuracy for a given shadow size versus a randomized classical shadow. More details on this algorithm can be found in [Phys. Rev. Lett. 127, 030503](https://arxiv.org/abs/2103.07510)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy from derandomized single-Pauli classical shadow: -1.1636 hartree\n",
      "Error vs statevector simulation: 0.016 hartree\n"
     ]
    }
   ],
   "source": [
    "from tangelo.toolboxes.measurements import DerandomizedClassicalShadow\n",
    "\n",
    "derandom_cs = DerandomizedClassicalShadow(vqe.optimal_circuit)\n",
    "derandom_cs.build(n_shots_budget, qubit_ham)\n",
    "derandom_cs.simulate(backend)\n",
    "\n",
    "energy_derandom = derandom_cs.get_observable(qubit_ham)\n",
    "print(f\"Energy from derandomized single-Pauli classical shadow: {energy_derandom:.4f} hartree\")\n",
    "print(f\"Error vs statevector simulation: {abs(energy_ref-energy_derandom):.3f} hartree\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adaptive Single-Qubit Pauli Classical Shadows\n",
    "Another procedure for choosing single-qubit Pauli bases is an adaptive one. It uses an iterative process to draw a Pauli gate for each qubit, one by one. The probability distributions are built so the next single-qubit Pauli gates drawn will depend on the measurement bases already assigned to previous qubits. This shifts the single-qubit Pauli basis choice towards terms hitting a given operator, while keeping a fair amount of randomness to extent the state characterization coverage. For more details on this algorithm, we refer to [arXiv:2105.12207](https://arxiv.org/abs/2105.12207)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy from adaptive single-Pauli classical shadow: -1.0768 hartree\n",
      "Error vs statevector simulation: 0.071 hartree\n"
     ]
    }
   ],
   "source": [
    "from tangelo.toolboxes.measurements import AdaptiveClassicalShadow\n",
    "\n",
    "adaptive_cs = AdaptiveClassicalShadow(vqe.optimal_circuit)\n",
    "adaptive_cs.build(n_shots_budget, qubit_ham)\n",
    "adaptive_cs.simulate(backend)\n",
    "\n",
    "energy_adaptive = adaptive_cs.get_observable(qubit_ham)\n",
    "print(f\"Energy from adaptive single-Pauli classical shadow: {energy_adaptive:.4f} hartree\")\n",
    "print(f\"Error vs statevector simulation: {abs(energy_ref-energy_adaptive):.3f} hartree\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Energy with Hamiltonian Partitioning\n",
    "For a given measurement budget, we compare the accuracy of total energy computation obtained with the Hamiltonian partitioning approach relying on qubitwise commutativity, to the energy predicted by classical shadow protocols. The `group_qwc` method defines the minimal number of bases to measure by identifying commutating terms. Then, the number of shots can be evenly distributed or according to a heuristic that considers the term coefficients.\n",
    "\n",
    "The reference value computed in the next cell only contains statistical errors and no quantum noise, \n",
    "even though it is always present in hardware. This shot-based simulation culminates with an energy containing an  error $\\epsilon$, inherent to the finite numbers of measurements. This error $\\epsilon$ will be compared to the ones coming from classical shadow protocols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy with qubitwise commutativity and equally distributed shots: -1.1619 hartree\n",
      "Error vs statevector simulation: 0.015 hartree\n"
     ]
    }
   ],
   "source": [
    "from tangelo.toolboxes.operators import QubitOperator\n",
    "from tangelo.toolboxes.measurements import group_qwc, exp_value_from_measurement_bases\n",
    "from tangelo.linq import Circuit\n",
    "from tangelo.linq.helpers import measurement_basis_gates\n",
    "\n",
    "# Performing Hamiltonian partitioning.\n",
    "qubitwise_measurements = group_qwc(qubit_ham)\n",
    "\n",
    "# Distributing the available number of shots evenly across identified bases.\n",
    "n_shots_per_basis = round(n_shots_budget / len(qubitwise_measurements))\n",
    "\n",
    "# Shot-based backend.\n",
    "backend = get_backend(\"cirq\", n_shots=n_shots_per_basis, noise_model=None)\n",
    "\n",
    "# Computing the histograms for the grouped bases.\n",
    "qubitwise_results = dict()\n",
    "for basis in qubitwise_measurements.keys():\n",
    "    basis_circuit = Circuit(measurement_basis_gates(basis))\n",
    "    qubitwise_results[basis], _ = backend.simulate(vqe.optimal_circuit + basis_circuit)\n",
    "\n",
    "# Computing the energy from the histograms.\n",
    "energy_shots = exp_value_from_measurement_bases(qubitwise_measurements, qubitwise_results)\n",
    "\n",
    "print(f\"Energy with qubitwise commutativity and equally distributed shots: {energy_shots:.4f} hartree\")\n",
    "print(f\"Error vs statevector simulation: {abs(energy_ref-energy_shots):.3f} hartree\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Averages of 10 runs\n",
    "Due to the nature of the classical shadow processes, there is a random component during prediction of observables. This is why your mileage may vary if you try to rerun all code cells in this notebook. To better portray the situation, we preran 10 simulations of 10 000 shots each for all the above methods, and compared the results to the ideal energy obtained with VQE run with an exact statevector simulator. The data is shown in the boxplot below and the ideal energy is represented by a dashed line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '(Noiseless Energy - Statevector Energy) / Hartree')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnMAAAFsCAYAAACn2D/dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAx9ElEQVR4nO3deZhkVX3/8feHbZAddUR2UVBcEjGMqIkLyho31KDggmhQxLhETRRFo4gQITEuiVFBQHFBRI2C/FQEFdyVGQUVN3aGVWSTff3+/ri3oabp7ro9U9VdPfN+PU89fe+5S32r6tbtb517zrmpKiRJkjQ3rTTbAUiSJGnpmcxJkiTNYSZzkiRJc5jJnCRJ0hxmMidJkjSHmcxJkiTNYSOTzCXZNckfkpyb5O0TLH9qkl8kuTPJ7uOW3ZXkzPZx4sxFLUmSNLsyCuPMJVkZ+COwE3AJcAbw4qr6bc86DwHWAf4VOLGqvtyz7MaqWmtGg5YkSRoBq8x2AK3tgHOr6nyAJMcBuwH3JHNVdWG77O7ZCFCSJGkUjcpl1o2BxT3zl7RlXa2eZGGSnyZ53kAjkyRJGmGjUjO3rDavqkuTPBT4bpJfV9V541dKsi+wL8Caa6657dZbbz3TcUqSJE3bokWL/lxV8ydaNirJ3KXApj3zm7RlnVTVpe3f85OcBjwOuE8yV1VHAEcALFiwoBYuXLgMIUuSJM2MJBdNtmxULrOeAWyVZIskqwF7Ap16pSZZP8m8dvqBwN/R09ZOkiRpeTYSyVxV3Qm8HjgZ+B1wfFWdneSgJM8FSPL4JJcALwQOT3J2u/kjgYVJzgK+Bxza2wtWkiRpeTYSQ5PMBi+zSpKkuSLJoqpaMNGykaiZkyRJ0tIxmZMkSZrDTOYkSZLmMJM5SZKkOcxkTpIkaQ4zmZMkSZrDTOYkSZLmMJM5SZKkOcxkTpIkaQ4zmZMkSZrDTOYkSZLmsFW6rphka5qb3D+4ql7Xzq9WVb8aWnSSJEmaUqeauSQvBL4PbAzs1RavBXxwSHFJkiSpg66XWQ8Cdqqq/YC72rKzgMcOJSpJkiR10jWZexAwdjm1ev7WxKtLkiRpJnRN5hZx7+XVMXsCPx9sOJIkSZqOrh0g3gh8O8k+wJpJTgYeDuw8tMgkSZLUV6dkrqp+3/ZefTZwErAYOKmqbhxmcJIkSZpa56FJqurmJD8BLq6qnw4xJkmSJHXUdWiSzZL8CPg9cGpbtnuSI4cZnCRJkqbWtQPE4cD/A9YG7mjLTgF2GkZQkiRJ6qbrZdbtgGdV1d1JCqCqrk+y7vBCkyRJUj9da+auBLbsLUjyKODigUckSZKkzromcx8ATkrySmCVJC8GvggcNrTIJEmS1FfXoUmOTnI18BqaYUleDvxbVX1tiLFJkiSpj77JXJKVge8Au1TVCcMPSZIkSV31vcxaVXcBW3RZV5IkSTOra4L2XuDjSTZPsnKSlcYewwxOkiRJU+s6NMnY4MB79ZQFKGDlgUYkSZKkzromc1sMNQpJkiQtla6XSV9YVReNfwD/MMzgJEmSNLWuydy7Jyl/16ACkSRJ0vRNeZk1yTPayZWTPJ2mndyYhwI3DCswSZIk9devzdxR7d/VgaN7ygu4AnjDMIKSJElSN1Mmc1W1BUCSz1fVS2cmJEmSJHXVt81ceweIFySZNwPxSJIkaRq63gHij8ADhh+OJEmSpqPrOHOfB05K8hHgEpo2cwBU1XeHEZgkSZL665rMvbb9e+C48qLp1SpJkqRZ0CmZG+sIIUmSpNHStWZO0nIi2+4z2yHcRy06qv9KkqQJdboDRJJ1knwwyaIkFyW5eOwxqECS7JrkD0nOTfL2CZY/NckvktyZZPdxy/ZOck772HtQMUmSJI26rjVzHwM2AQ4CPge8DHgr8JVBBNEOf/K/wE40HSzOSHJiVf22Z7WLgVcA/zpu2/sD7wEW0LThW9Rue+0gYpOWN4OoBRur3bNGTZJmX9d7s+4M/ENVnQDc1f7dA9hrQHFsB5xbVedX1e3AccBuvStU1YVV9Svg7nHb7gKcUlXXtAncKcCuA4pLkiRppHVN5lYCrm+nb0yyLnA5sOWA4tgYWNwzf0lbNtBtk+ybZGGShVddddVSBSpJkjRKuiZzZwFPa6d/QHPZ9eM0gwnPGVV1RFUtqKoF8+fPn+1wJEmSllnXZO7VwIXt9D8DtwDrAS8fUByXApv2zG/Slg17W0mSpDmt6zhz5/dM/wl41YDjOAPYKskWNInYnsBLOm57MvDvSdZv53cG3jHg+CRJkkbSlMlckmf028EgbudVVXcmeT1NYrYycHRVnZ3kIGBhVZ2Y5PHAV4H1geckeW9VPbqqrknyPpqEEOCgqrpmWWOSJEmaC/rVzI0fd2BTluxsMLDbeVXVN4BvjCt7d8/0GTSXUCfa9mjg6EHEIUmSNJdMmcyNv41Xkmu9tZckSdLo6NoBYkwNJQpJkiQtlekmc5IkSRohJnOSJElzWL/erD9gyUurayf5fu86VfXUYQQmSZKk/vr1Zj1y3Lx31ZYkSRoh/XqzHjNTgUiSJGn6bDMnSZI0h3W6nZek0XDrbXew+rxVZzsMatFotbgYlfdFkmaDyZw0h6w+b1Wy7T6zHcbIGbXkUpJm0pSXWZMcnORvk2SmApIkSVJ3/drM3QgcBlye5PNJXprkATMQlyRJkjqYMpmrqkOr6inA1sCJwC7Ab5L8JMm7kvzNTAQpSZKkiXXqzVpV11XVF6vq5cBGwD/TtLc7PMklSfYYZpCSJEma2LQ7QFRVAT9vHwcmeRCw7qADkyRJUn/L3Ju1qv4E/GkAsUiSJGmaHDRYkiRpDjOZkyRJmsNM5iRJkuawvm3m2uFHngU8FlgPuA44C/hmVS0cZnCSJEma2qQ1c0l2SbIQ+AKwKfAj4Lj276bA55MsTLLrjEQqSZKk+5iqZu7VwGur6ozJVkjyeGB/4FuDDkySJEn9TZrMVdXu/TZuE72+60mSJGk4OnWASPKhJNsMORZJkiRNU9ferCsDJyf5TZL9k2wyzKAkSZLUTdd7s76R5p6sbwe2AX6X5NQkL0+y1hDjkyRJ0hQ6jzNXVXdV1UlV9WLgicB84NPAFUmOTLLxkGKUJEnSJDonc0nWSbJPku8B3wd+BjwFeCRwI/DN4YQoSZKkyfQdNBggyZeBXWiSuE8AX6uq23qWvwW4figRSpIkaVKdkjngp8Drq+qKiRZW1d1JNhhcWJIkSeqiUzJXVR/osM7Nyx6OJEmSpqPrZdbFQE2w6DbgEuD/gI9X1Z0DjE2SJEl9dL3M+t/Ay9q/i4HNgNcBXwKuAf6F5n6tbxtCjJIkSZpE12TuFcBOVXXZWEGSbwLfrqpHtz1cT8VkTpIkaUZ1HZpkQ5rhR3rdRDOQMMAfgfUGFJMkSZI66loz93XghCSH0LSR2wR4R1sO8CTgwoFHp06y7T6zHcJ91KKjZjsESZJWCF1r5l5DM0jw4cAv279nAPu1y88HnjXw6CRJkjSlvjVzSVamSd72raq3T7TOZOPPaWYMqhZsrIbPWjVJkuaOvjVzVXUXsDNw9/DDkSRJ0nR0vcz6IeC9SVYbZjCSJEmanq7J3BuAtwJ/SbI4ycVjj0EFkmTXJH9Icm6S+1zOTTIvyRfb5T9L8pC2/CFJbklyZvv4xKBikiRJGnVde7O+bJhBtO3y/hfYiaa37BlJTqyq3/astg9wbVVtmWRP4DBgj3bZeVW1zTBjlCRJGkVd7816+pDj2A44t6rOB0hyHLAb0JvM7QYc2E5/Gfhokgw5LkmSpJHW6TJre4nzkCTnJ7m+Lds5yesHFMfGNLcJG3NJWzbhOu09YK8HHtAu2yLJL5OcnuQpU7yOfZMsTLLwqquuGlDokiRJs2c6HSAeA7wUqLbsbOC1wwhqmi4HNquqxwFvAY5Nss5EK1bVEVW1oKoWzJ8/f0aDlCRJGoaubeaeD2xZVTcluRugqi5NMr72bGldCmzaM79JWzbROpckWQVYF7i6qgq4rY1pUZLzgIcDCwcUmyRJ0sjqWjN3O+MSvyTzgasHFMcZwFZJtmiHP9kTOHHcOicCe7fTuwPfrapKMr/tQEGShwJb0dyRQpIkabnXNZn7EnBMki0AkmwIfBQ4bhBBtG3gXg+cDPwOOL6qzk5yUJLntqsdBTwgybk0l1PHhi95KvCrJGfSdIzYr6quGURckiRJo67rZdYDaIYC+TWwBnAO8EngoEEFUlXfAL4xruzdPdO3Ai+cYLuvAF8ZVBySJElzSdehSW4H3gy8ub28+ue2rZokSZJmUdeaOZKsCzwCWKudB6CqvjuUyCRJktRXp2QuySto7tBwI3Bzz6ICHjr4sCRJktRF15q5Q4Ddq+qbwwxmRXXrbXew+rxVZzsMatFRsx3CPUblPZEkadR1TeZWAb49zEBWZKvPW5Vsu89shzFSRimxlCRplHUdmuQw4F1Juq4vSZKkGdC1Zu7NwIOBtyVZYqDgqtps4FFJkiSpk67J3MuGGoUkSZKWStdx5k4fdiCSJEmavinbwCU5Ydz8e8fNnzGMoCRJktRNvw4NTx83/4Zx81sPMBZJkiRN03R7p2bcvLf0kiRJmkXTTeZM3iRJkkZIvw4QqyZ5JffWyM1L8o/T2F6SJElD1C8Z+xnw8p75nwN7jVsuSZKkWTJlMldV289QHJI68lZnkqRe3p5LkiRpDrPNmzTHZNt9ZjuEkWNtpaQVmTVzkiRJc1jfZC7JSkmekWS1mQhIkiRJ3fVN5qrqbuCEqrp9BuKRJEnSNHS9zPr9JE8caiSSJEmatq4dIC4CvpnkBGAxPXeCqKp3DyMwSZI0d4xi56wVpXNU12TufsDX2ulNesq9vZckSdIs6pTMVdUrhx2IJEmauwZRCzZWu7ei1KgNSudx5pJsBbwY2Bi4FPhCVZ0zrMAkSZLUX6cOEEmeAywCtgauAR4BLEzy3CHGJkmSpD661sz9O7BbVX1vrCDJ9sBHgRMHH5YkSZK66JrMbQL8YFzZD1myM4SWge0DJEnS0ug6ztyZwL+MK3tLWy5JkqRZ0rVm7p+AE5P8M804c5sCNwPPGVZgK5pRHJ9nNllTKUlSN12HJvldkkcCTwQ2Ai4DflZVdwwzOEmSJE2tUzKX5ISq2o2mnVxv+f9V1QuGEpkkSZL66nqZ9emTlG8/oDgkSSNmFJt/2ARDuq8pk7kkB7WTq/VMj3kozT1bJUmSNEv61cxt2v5dqWcamnuyLgYOHEJMkqQRMKhaMG/RJA3XlMnc2D1Zk/y4qj45MyFJkiSpq67jzN2W5K97C5I8NsleQ4hJkiRJHXVN5t5Hc1m112Lg4MGGI0mSpOno2pt1HeAv48quB9YbaDQrqFtvu8O2JOPcetsdrD5v1dkOQ5rTRuV7NErnt1F5T6RB6prM/Rb4B+D4nrLnA78bVCBJdgU+AqwMHFlVh45bPg/4DLAtcDWwR1Vd2C57B7APcBfwxqo6eVBxzYRRObGMUiPlUXlPpLls9XmrjuTwIrNpFM5v0qB1Teb2B76RZA/gPGBLYAfgmYMIIsnKwP8COwGXAGckObGqftuz2j7AtVW1ZZI9gcOAPZI8CtgTeDTN3SlOTfLwqrprELFJkiSNsk5t5qrqh8BfAWcAawI/Bx5TVT8aUBzbAedW1flVdTtwHLDbuHV2A45pp78M7JAkbflxVXVbVV0AnNvuT5IkabnXtQMEVXUR8B/AwVV1aFWN7xCxLDZmyQ4Wl7RlE65TVXfStNl7QMdt72PRokVLzG+77bYkmfCx7777LrHdZOslWWK/++6776Trbbvttks8/1T7POKII+5Z74gjjphy3WV5TZPFMpdf0/L2Oamb2f6cRuXY0+RG6XNaHo+9pX1Nfk6Tv6apdL0363rAx4DdgTuANZM8F9iuqt7VZR+jIMm+wL59V5QkaQVywYUX85DNN51wWVV12se+++67RNIyla4/TrfddtvOz3/EEUcskYhNpes+d3/RHrzmNa/ptO5sSpcXlOQ44FrgIOC3VbV+kvnAj6tqq2UOInkScGBV7dLOvwOgqt7fs87J7To/SbIKcAUwH3h777q96031nAsWLKiFCxcua+jLlVHqAKHJ2aD9vjxmJ+fxsiSPlcl5rNzXKB0vSRZV1YKJlnW9zLoDTS/Ry2lu5UVVXQU8aDAhcgawVZItkqxG06HhxHHrnAjs3U7vDny3mkz0RGDPJPOSbAFsRdOmT5IkabnXtTfr9cADgcvHCpJs1ju/LKrqziSvB06mGZrk6Ko6O8lBwMKqOhE4CvhsknOBa2gSPtr1jqcZPuVO4HX2ZJWkxijVLEgajq7J3JHAV5K8E1ipvSz678AnBhVIVX0D+Ma4snf3TN8KvHCSbQ8BDhlULJIkSXNF12TuMOAWmrHgVgWOBg6nGeRXs2zQ7RwGsT9rA6TRYDuoJXlu0vKoUzLXtk37CCZvkiRJI6Xr0CRnAZ8DvlBVlww3JE2XvzQlSVpxde3NeiDweOD3SU5P8pok9x9eWJIkSeqi6+28vlpVLwI2pGkv93xgcZLxw4dIkiRpBnXtAAFAVd2Q5FjgOmA14JnDCEqSJEnddG0zF+AZwEtoauUuAo7l3kF8JUkj5tbb7rBN7Ti33nYHq89bdbbDkAaqa83cZcCNwHHA31XV74YXkiRpEEYlaRmlWwWOynsiDVLXZG63qvIWWZIkSSOm6zhzP0+yNc0dGB5cVa9r51erql8NNUJJkiRNqlNv1iQvBH4AbAzs1RavBXxwSHFJkiSpg67jzB0E7FhV+wFjN7E/C3jsUKKSJElSJ12TuQcBY5dTq+dvTby6JEmSZkLXZG4R915eHbMnYKcISZKkWdS1N+sbgW8n2QdYM8nJwMOBnYcWmSRJmjGjMHSMlk7X3qy/b3uvPhs4CVgMnFRVNw4zOEmSJE2t6x0g/ruq3ggcP678w1X1pmEEJkmSZs7Y4M6611yprex6mfUVNJdax9sLeNOggpEkjY5B/3MfxP7myj9XaSZNmcwl+cex9XqmxzwU+PNQopIkSVIn/WrmxnqwrsaSvVkLuBLYexhBSZJmn7Vg0twwZTJXVU8HSHJwVb1rZkKSJEkz6dbb7jB5n8Ctt93B6vNWne0w+uram/WeRC5JgPQsu3sIcUmSpBkyKgnLWLvKUUksR+V96afrvVk3SvLVJFcDdwJ39DwkSZI0S7reAeJw4HZgB+BG4G+AE4H9hhSXJEmSOug6NMnfAptV1U1JqqrOau8G8WPgk8MLT5IkSVPpWjN3F83lVYDrkswHbgI2HkpUkiRJ6qRrzdzPgGcCXwVOBr4I3AIsHFJckiZgj7OJzZUeZ5I0DF2Tub24txbvTcC/AGsDHxpCTJImMSoJiz3OJGl0dL3MulNVXQNQVbdU1cFVtT/w5OGFJkmSpH66JnOT/fw+YlCBSJIkafr63Zv1oe3kSkm2oGewYJp7s946rMAkSZLUX782c+fS3Ic1wHnjll0BHDiEmCRJktRRv3uzrgSQ5PSqetrMhCRJkqSuOrWZ603kkqyX5HFJ1hheWJIkSepiymQuyduSvKBnfldgMbAIWJzkiUOOT5IkSVPoVzP3j8Bveub/u32sDXwQ+PchxSVJkqQO+iVzG1bVHwGSbAlsDry/qm4CPgD89ZDjkyRJ0hT6JXM3J1mnnX4y8KuqurGdv5vud5CQJEnSEPRL5r4BHJHkucC/Al/pWfZYmvZzkiRJmiX9krm3ADcDhwA/Ycl7se4KHDekuCRJktRBv3HmrqfpBDHRsoOHEpEkSZI663pv1qFJcv8kpyQ5p/27/iTr7d2uc06SvXvKT0vyhyRnto8HzVz0kiRJs2vWkzng7cB3qmor4Dvt/BKS3B94D/AEYDvgPeOSvpdW1Tbt408zEbQkSdIoGIVkbjfgmHb6GOB5E6yzC3BKVV1TVdcCp9C02ZMkSVqhjUIyt0FVXd5OXwFsMME6G7Nkz9lL2rIxn2ovsf5bkkz2REn2TbIwycKrrrpqmQOXJEmabdNO5pL8eim2OTXJbyZ47Na7XlUVUNPc/Uur6q+Ap7SPvSZbsaqOqKoFVbVg/vz5030ZkiRJI2dpBv3dfLobVNWOky1LcmWSDavq8iQbAhO1ebsU2L5nfhPgtHbfl7Z/b0hyLE2bus9MN0ZJkqS5aGkus056GXMpnQiM9U7dGzhhgnVOBnZOsn7b8WFn4OQkqyR5IECSVYFns+S9ZCVJkpZrS1Mz9/cDjuFQ4Pgk+wAXAS8CSLIA2K+qXlVV1yR5H3BGu81BbdmaNEndqsDKwKnAJwccnyRJ6iPb7jNy+6pFRw1kP6Nu2slcVf1wkAFU1dXADhOULwRe1TN/NHD0uHVuArYdZDySJElzydLUzEmSJC1hRakFG0WjMDSJJEmSlpLJnCRJ0hzWKZlL8tUkz2s7GkiSJGlEdK2Z+wHwbuCKJB9P8rdDjEmSJEkddUrmquqDVfU3wFOB64AvJDknybuTPGyYAUqSJGly02ozV1VnV9U7gJcBNwPvAX7R3q7rscMIUJIkSZPrnMwleUSS9yU5DzgC+CLwEGAD4BvA14YRoCRJkibXaZy5JAtpErcvAi+pqp+NW+WDSd4w4NgkSZLUR9dBgw8FTqyq2ydboaq2GExIkiRJ6qprMvd/AEmWuCxbVXcPPCJJkiR11rXN3J3AHeMfSW5LckGS/0qy1rCClCRJ0sS6JnNvAL4L7Aw8EtgF+A7wNuC1wN8CHx5CfJIkSZpC18usbwH+pqqub+f/2HaKWFRVD0vya2DRUCKUJEnSpLrWzK0DrDGubA1g3Xb6CuB+gwpKkiRJ3XStmfsMcEqSjwCLgU2AfwaOaZfvDPxh8OFJkiRpKl2TubcC5wB7AhsBlwP/C3yyXf494LRBBydJkqSp9U3mkqxM09lhl6r6xETrVNWtgw5MkiRJ/fVtM1dVdwFbABl+OJIkSZqOrh0g3gt8IsnmSVZOstLYY5jBSZIkaWpd28wd2f7dq6csQAErDzQiSZIkddY1mfO+q5IkSSOoUzJXVRfBPfdm3aCqLh9qVJIkSeqkU5u3JOslORa4FTi3LXtukoOHGZwkSZKm1rUDwyeA64HNgdvbsp8AewwjKEmSJHXTtc3cDsBGVXVHkgKoqquSPGh4oUmSJKmfrjVz1wMP7C1IshnNnSAkSZI0S7omc0cCX0nydGClJE+iuS/rhHeEkCRJ0szoepn1MOAWmvuxrgocDRwOfGRIcUmSJKmDrkOTFE3iZvImSZI0QrrWzJHkEcBjgbV6y6vq6EEHJUmSpG46JXNJDgDeDZwF3NyzqGguuUqSJGkWdK2ZexOwXVX9aoixSJIkaZq69ma9Bfj9MAORJEnS9HVN5v4N+J8kGyZZqfcxzOAkSZI0ta6XWT/d/n1VT1lo2sytPMiAJEmS1F3XZG6LoUYhSZKkpdJ1nLmLhh2IJEmSpm/KNm9Jfjlu/lPj5v80jKAkSZLUTb8ODFuOm99t3Pz9BhiLJEmSpqlfMlfj5tNn+bQluX+SU5Kc0/5df5L1vpXkuiQnjSvfIsnPkpyb5ItJVlvWmCRJkuaK6Q4tsszJ2wTeDnynqrYCvtPOT+Q/gb0mKD8M+FBVbQlcC+wzhBglSZJGUr8OEPOSHNQzf79x84OoBdsN2L6dPgY4Ddh//EpV9Z0k2/eWJQnwDOAlPdsfCHx8AHFJkiSNvH7J3LHApj3zx42b/8IAYtigqi5vp68ANpjGtg8ArquqO9v5S4CNBxCTJEnSnDBlMldVrxzEkyQ5FXjwBIveOe75KskwLuWOxbEvsC/AZpttNqynkSRJmjFdBw1eJlW142TLklyZZMOqujzJhsB0hju5GlgvySpt7dwmwKVTxHEEcATAggULhpY0SpIkzZRRuLfqicDe7fTewAldN6yqAr4H7L4020uSJM11o5DMHQrslOQcYMd2niQLkhw5tlKSHwBfAnZIckmSXdpF+wNvSXIuTRu6o2Y0ekmSpFk0I5dZp1JVVwM7TFC+EHhVz/xTJtn+fGC7oQUoSZI0wjrVzCV5epIt2ukNkxyT5FNJJurUIEmSpBnS9TLrx4C72un/AlYF7qbtTCBJkqTZ0fUy68ZVdXGSVYBdgM2B24HLhhaZJEmS+uqazP0lyQbAY4DfVtWN7T1QVx1eaJIkSeqnazL3P8AZNLfvelNb9nfA74cQkyRJkjrqlMxV1WFJvgrcVVXntcWX0tPbVJIkSTOv89AkVfXHsekkTwfurqrThxKVJEmSOuk6NMnpSf6und4fOA44NskBwwxOkiRJU+s6NMljgJ+2068Gng48EdhvGEFJkiSpm66XWVcCKsnDgFTVbwGSrD+0yCRJktRX12Tuh8BHgQ2BrwK0id2fhxSXJEmSOuh6mfUVwHXAr4AD27KtgY8MPCJJkiR11nVokquBA8aV/b+hRCRJkqTOuvZmnZfkkCTnJ7m+Lds5yeuHG54kSZKm0vUy64doerS+FKi27GzgtcMISpIkSd107QDxfGDLqropyd0AVXVpko2HF5okSZL66VozdzvjEr8k84GrBx6RJEmSOuuazH0JOCbJFgBJNqQZquS4YQUmSZKk/romcwcAFwC/BtYDzgEuA947nLAkSZLURdehSW4H3gy8ub28+ueqqj6bSZIkacgmTeaSPHSK7dZOAkBVnT/ooCRJktTNVDVz59IMQ5Ip1ilg5YFGJEmSpM4mTeaqqmt7OkmSJM2SaSVsSTZN8sRhBSNJkqTp6Xo7r82S/Aj4PXBqW7Z7kiOHGZwkSZKm1rVm7nDg/wFrA3e0ZacAOw0jKEmSJHXT9XZe2wHPqqq7kxRAVV2fZN3hhSZJkqR+utbMXQls2VuQ5FHAxQOPSJIkSZ11TeY+AJyU5JXAKkleDHwROGxokUmSJKmvrneAODrJ1cBrgMXA3sC/VdXXhhibJEmS+ujaZo6qOgE4YYixSJIkaZq6Dk3y4iSPbKcfnuT0JN9LsvVww5MkSdJUuraZOxi4pp3+L+AM4HTgY8MISpIkSd10vcw6v6quTLI68GRgd5rx5v48tMgkSZLUV9dk7qokWwJ/BZxRVbclWQPI8EKTJElSP12TufcBi4C7gD3ash2Bs4YRlCRJkrrpOjTJp5Mc307f3Bb/FNhzWIFJkiSpv0mTuSSpqmqnVwJu7ZkG28tJkiTNuqlq5q4H1mmn7wRq3PK0ZSsPIS5JkiR1MFUy9+ie6S2GHYgkSZKmb9JkrqoW90xfBPdcYt0AuLKq7h5+eJIkSZpK1ztArJPkMzTt5i4FbklyTJJ1lzWAJPdPckqSc9q/60+y3reSXJfkpHHln05yQZIz28c2yxqTJEnSXNH1DhD/DawJPAa4H814c2u05cvq7cB3qmor4Dvt/ET+E9hrkmVvrapt2seZA4hJkiRpTug6ztyuwEN7hiX5Y5JXAucNIIbdgO3b6WOA04D9x69UVd9Jsv34ckmSpBVZ15q5W4H548oeCNw2gBg2qKrL2+kraNrkTdchSX6V5ENJ5k22UpJ9kyxMsvCqq65aqmAlSZJGSdeauSOBU5J8ELgI2Bx4M3BEl42TnAo8eIJF7+ydqapKMn4IlH7eQZMErtbGsz9w0EQrVtUR7TosWLBgus8jSZI0cromc4cAlwEvATZqp/8DOLrLxlW142TLklyZZMOqujzJhsCfOsY0tu+xWr3bknwK+NfpbC9JkjSXdb2dV9Ekbp2St2k6EdgbOLT9e8J0Nu5JBAM8D/jNwCOUJEkaUVMmc0le3m8HVfWZZYzhUOD4JPvQXMJ9UfvcC4D9qupV7fwPgK2BtZJcAuxTVScDn08yn+aOFGcC+y1jPJIkSXNGv5q5V09SXsAjgfsDy5TMVdXVwA4TlC8EXtUz/5RJtn/Gsjy/JEnSXDZlMjdRApXkr4H3tbOTjQknSZKkGdB1aBKSbJXkCzTjwC2iGXfuP4cVmCRJkvrrm8wl2SzJ0cBC4GJgy6o6qKpuGHp0kiRJmtKUyVySjwK/Bm4Atqqq/avqmhmJTJIkSX316wDxT8BNwAuA5zejfyypqjYbQlySJEnqoF8y9/QZiUKSJElLpV9v1tNnKhBJkiRN36TJXJI3AodX1W1TrDMPeE1V/fcwgpM0eNl2n5HbVy06aiD7kaQV0VQ1cw8Gzk3yDeB04A80HSHWBh4ObA/8Pcs4aLAkSZKW3qTJXFUdkOSDwCuAfYC/AtYDrgV+BXwDOKC9g4OkOcJaMElavvRrM/dn4APtQ5IkSSOm8x0gJEmSNHpM5iRJkuYwkzlJkqQ5zGROkiRpDjOZkyRJmsNM5iRJkuYwkzlJkqQ5zGROkiRpDjOZkyRJmsNM5iRJkuYwkzlJkqQ5zGROkiRpDjOZkyRJmsNSVbMdw6xIchVw0WzHMYIeCPx5toPQnOCxounweFFXHisT27yq5k+0YIVN5jSxJAurasFsx6HR57Gi6fB4UVceK9PnZVZJkqQ5zGROkiRpDjOZ03hHzHYAmjM8VjQdHi/qymNlmmwzJ0mSNIdZMydJkjSHmcx1lOTCJDtOsuwpSf4wgOcYyH4ESV6a5Nsd1z07yfbDjWjuSHJgks8NYb+bJbkxycoD3u+k303NrmF9NkkOSHLkgPe5fZJLBrnPFVmSTyc5eEj77nx+X1GsMMlcklck+XWSm5NckeRjSdYdxL6r6gdV9Yie51qqE9j4/ayo2s/qh9NY/yFJKskqY2VV9fmq2rnL9lX16Ko6rd3XUBKZQWiPq1vahOiK9mS51mzH1VVVXVxVa1XVXbMdy/Ki55i4Icl1SX6cZL8ky/W5var+vapeNdtxrKiSnJbk2iTzZuC5lun8vqJYrr/wY5L8C3AY8FZgXeCJwEOAbydZdRZDk6brOVW1FrAN8DjgHbMbjkbAc6pqbWBz4FBgf+Co6e6k95+lNJkkDwGeAhTw3NmNRmOW+2QuyTrAe4E3VNW3quqOqroQeBHwUOAl7XpLVAlPUuX++CS/bX+RfCrJ6uPXTfJZYDPg620NytuSHNMmlCTZuP2V8bp2/mFJrkmy0vjnTLJ/kkvbX91/SLJDW75SkrcnOS/J1UmOT3L/Kd6D3ZKcmeQv7Ta7tuUbJTmxff5zk7y6Z5sDk3wpyefa5/91kocneUeSPyVZnGTnnvVPS3JwWzNwY5KvJ3lAks+3z3tGexKY8JdWu/2rkjwS+ATwpHY/17XLn5Xkl+2+Fic5sOclfr/9e127zZN6a/eSfDzJB8a9JyckeUs7fWGSHdv35QBgj3Y/ZyV5YZJF47Z9S5ITJnu/Z0JVXQGcTJPU0XM83NAeo88fW3fsvUjygfbYvSDJ3/cs3yLJ6e22p9CMvk7P8uemuRR9Xfs5PbJn2YVJ3prkV0luSnJUkg2SfLPd36lJ1m/Xvedzbz+jG3setya5sF1vyuM7yV5JLmqXvXMob/AcVFXXV9WJwB7A3kkek2Re+7lfnOTKJJ9Icj+497yV5jxzBfCpJOsnOSnJVe2xclKSTcaeo/3835fkR+3n++0kD+xZPuln08by4SSXtY8Pp63Z6YnlbWnOL5cneV6SZyb5Y5pz1AE9+7qnBj3JR8cdS3eOnR/SnOO+0r6eC5K8sWcf90tz3r82yW+Bxw/hY1kevRz4KfBpYO+xwiSPS/KL9rj4IrB6z7Iux9X7k/w8zTn+hJ7v/LKe3yc9BpYrVbVcP4BdgTuBVSZYdgzw+Xb608DBPcu2By7pmb8Q+A2wKXB/4Edj60+y7o498/8IfL2dfglwHvDFnmUnjN8P8AhgMbBRO/8Q4GHt9D/TfJk2AeYBhwNfmOT1bwdcD+xEk7xvDGzdLvs+8DGaL902wFXAM9plBwK3ArsAqwCfAS4A3gmsCrwauKDneU4DzgUeRlP7+Vvgj8COPdt/que1VO9n0m7/qnb6FcAPx72O7YG/al/DXwNXAs+bYn/37AN4avtejvXeXh+4pee9vefzal/353r2Mw+4BnhkT9kvgX+YhWO5N85NgF8DH2nnXwhs1L4/ewA3ARv2vBd3tJ/ZysBrgct63o+fAB9sX+tTgRvG3gPg4e2+dmo/97e1n/NqPTH9FNiA5tj6E/ALmlrD1YHvAu+Z7HNqy1cFTgfe3+/4Bh4F3NjGOa+N+056vm8r0oNx55qe8ovbz/lDwIk056y1ga/3vM/bt+/dYe17eT/gAcA/AGu0638J+Nq47+l57XFxv3b+0C6fDXBQ+7k+CJgP/Bh437hY3s2955ergGPbOB5N853dYqLvaU9827TbPY7mu7Co3edqND/ezwd2adc9FPhB+95sSnN+v2RZP5Pl/dF+//8J2JbmvLJB+/5eBLy5/fx2b5eN/Y/sclxdCjwGWBP4Cveegx7CUp7f+x0Dy9Nj1gOYgQPvZcAVkyw7FPh2O/1p+idz+/XMPxM4b4p1e5O5hwHXtgfWJ4DXcG/SdgzwlvH7Abak+ce4I7DquLh/B+zQM79h+8WZKGE9HPjQBOWbAncBa/eUvR/4dDt9IHBKz7Ln0JyoV27n126/YOu186cB7+xZ/7+Ab47b/sx2eqIv52lMkcxNEP+Hx15Xhy97aP65PbWdfzXw3Yk+Lyb4JwF8HDiknX50+1nOm4Vj+cL2M7ihfb3fGXv/J1j3TGC3nvfi3J5la7TbP5imFvlOYM2e5cdy74n034Dje5atRHPS3b4nppf2LP8K8PGe+TfQnrQn+px63t+TgJX6Hd80J+XjepatCdyOydz48p/S/PC6ifZHYFv+JNofYTTnm9uB1afY/zbAtT3zpwHv6pn/J+Bb7fSUnw1NEvjMnuW7ABf2xHIL9z2/PKFn/UXc+wPuQO77PZ3fvh97tvNPAC4et847uPdH5fnArj3L9sVkrt/x9uT2u/jAdv73NAncU+n5gdgu+zE9/1M7HFeH9sw/qj12Vp7ovEHH83u/Y2B5eiz3l1lpbtb7wEzcHmRDpncz38U90xfRZP59VdV5NCfVbWjaGpwEXJbkEcDTaGolxm9zLvAmmpPWn5Icl2Ts+TYHvprmstd1NP/87qL5hTTepjQn0fE2Aq6pqhvGvaaNe+av7Jm+Bfhz3dt4/Zb271pTrD9+fqkb6yd5QpLvtVXl1wP7Me5y4GSq+QYfB7y4LXoJ8PlpPP0xwEuSBNiLJrm5bRrbD9LzqmkftT2wNe17kOTlaS6ljx0Tj2HJ9+eKsYmqurmdXIvmOLi2qm7qWfeinumNeuer6m6a78FUx0nnzz3Ja9rX8pJ23zD18b0RPd/DNu6rJ9v/CmxjmuR3DWBRz3v5LZqkZ8xVVXXr2EySNZIc3l4q/QtN7f16WbIH8hU90zdz7+fb77NZ4ljivufQqyc4v3Q6ltK0ff4ycGxVHdcWbw5sNPba29d/APeeJ5eId1xsmtjeNBUgY/83j23LNgIubc+1Y+55PzseV+M/i1XpcI7vc37vdwwsN1aEZO4nwG3AC3oL0/QC/HuaXwTQJFtr9Kzy4An2tWnP9GY0v0QmUhOUnU5T9bxaVV3azu9NUyV85oQ7qTq2qp5Mc0AWzeUQaA76v6+q9Xoeq7f7HW8xTc3geJcB90+y9rjXNNE+Bm0scZjs/Z7o/TuW5nLRplW1Lk0NZ6ZYf7wvALsn2Zzm19pXJlnvPvuqqp/S/Ep8Cs2J4rMdnm+oqup0mtrkD7Sv6ZPA64EHVNV6NJeMMukO7nU5sH6SNXvKNuuZvozm+AOgTWg3ZQDHSZKnAO+jqUH8S8+iqY7vy+n5HiZZg+YSjlpJHk+TzH2NJgF6dM/7uG41HWjGjD/e/4WmiccTqmodmhoX6H4sTfXZLHEsMfU5dLr+B/gL8K6essU0tZC9x9HaVfXMieJlyeNe46Rpa/ki4GlpetNfQVMr91ia93Lj9vwwpvf97HJcjf8s7qCpbFmW83u/Y2C5sdwnc1V1PU0HiP9JsmuSVdM0xD+e5kAZy+DPBJ6Z5P5JHkxTKzbe65Js0jbMfCfwxUme9kqaa/O9Tqf5ZzvWmPO0dv6HNcFQDUkekeQZbQPhW2lOymM1F58ADmkPXJLMT7LbJLEcBbwyyQ5pGpZvnGTrqlpMUw3+/iSrJ/lrYB9g6MNyVNVVNMnAy5KsnOQfWTLhvBLYJMlqPWVr09Qk3ppkO9qOK62raN6b8e9573P+kubzPhI4uaqum2TVK4GH5L5DO3wG+ChwR1V1HjZlyD5M05ZtPZoT3lUASV5JUzPXV1VdBCwE3ptktSRPprkkPuZ44Fnt8bMqzUn5NppjZ6kl2bTd98ur6o/jFk91fH8ZeHaSJ7fHx0GsAOexLpKsk+TZNLUUn6uqs2iS/A8leVC7zsZJdpliN2vTnGuua89z75lGCP0+my8A72o/zwfSXJZd5vNNW7v7NJrL/Xf3LPo5cEOaDh73a881j2mTXWiOv3ekaZy/CU2TAE3ueTQ15I+iucq0DfBImnaHz6NprvHG9n/sC2jaa4/pcly9LMmj2h8BBwFfbv83Lsv5vd8xsNxYIU6CVfUfNFWrH6Bpb3QBTa3Qjj2Xlz4LnEXT5uLbTJyoHdsuO5/m0uVkAyK+n+akdV2Sf23LTqc5oMeSuR+2MXx/gu2haUB8KM0BegVNo+GxYSg+QlNL9e0kN9C0j3nCJK/958AraRpCX9/GMfbr+MU07REuA75K01D91EniGbRX0wwVczVNO7Te5OC7wNnAFUnGqvP/CTiofb3vpjkRA/dcNjwE+FH7nj9xkuc8lqYN4rFTxPWl9u/VSX7RU/5ZmgRpZMaga5Piz9C8H/9FUwt9JU1HkR9NY1cvoTl+rqE5yX6m5zn+QNPu9H9ojsXn0AyFcfsyhr8DzaWOL+feXohnt8smPb6r6mzgdTSf4eU07RdX9IFev96+T4tpfmR+kOY7D80wJecCP20vb51KU0MymQ/TdGz4M837/q2uQXT4bA6m+eHwK5rOO79g8nPodLyY5h/9ZT3H0gFtIvBsmqTjAu79Z79uu917aS7nXUBzXp/1GvcRtzdNW7OLq+qKsQfNj9wX01z9egXNeWQP4P96tv0w/Y+rz9JcbbiCpvPUG2HZzu8djoHlxgp5b9a25uIg4O+q6uLZjkejr73E8Cfgb6rqnNmOR5KWF0lOo6lNHuhdPVYkK+QgkVX1qSR3An9L0wtG6ue1wBkmcpKkUbNCJnMAVWWVujpJM5htaNqFSJI0UlbIy6ySJEnLixWiA4QkSdLyymROkiRpDjOZkyRJmsNM5iRJkuYwkzlJkqQ5zGROkiRpDvv/Ns6H55RkHOEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "data = np.array([\n",
    "    [ 0.00508,  0.02309, -0.03742, -0.02013],\n",
    "    [ 0.01601, -0.16360, -0.02321, -0.02002],\n",
    "    [ 0.04252,  0.03412,  0.00546, -0.05518],\n",
    "    [ 0.01322,  0.10274,  0.00458,  0.02129],\n",
    "    [ 0.05665, -0.01826, -0.00687, -0.01218],\n",
    "    [ 0.05709, -0.08834, -0.01307,  0.00555],\n",
    "    [-0.00185,  0.07029,  0.00910, -0.02183],\n",
    "    [ 0.03663,  0.09825, -0.00555,  0.00459],\n",
    "    [-0.00593,  0.09175, -0.00332, -0.04029],\n",
    "    [ 0.00574, -0.03616,  0.00986, -0.02638]\n",
    "])\n",
    "\n",
    "_, ax = plt.subplots(figsize=(10,6))\n",
    "VP = ax.boxplot(data, positions=[2, 4, 6, 8], \n",
    "                widths=1, patch_artist=True,\n",
    "                showmeans=False, showfliers=False,\n",
    "                medianprops={\"color\": \"white\", \"linewidth\": 2},\n",
    "                boxprops={\"facecolor\": \"#003262\", \"edgecolor\": \"white\", \"linewidth\": 1},\n",
    "                whiskerprops={\"color\": \"#003262\", \"linewidth\": 2},\n",
    "                capprops={\"color\": \"#003262\", \"linewidth\": 2})\n",
    "\n",
    "ax.set(ylim=(-.18, .15), xlim=(1, 9))\n",
    "ax.hlines(y=0., xmin=0, xmax=10, linewidth=2, color=\"k\", linestyles=\"dashed\")\n",
    "\n",
    "ax.set_xticklabels(labels=[\"Qubitwise commutativity\", \"Randomized\", \"Derandomized\", \"Adaptive\"],fontsize=12)\n",
    "ax.set_ylabel(r\"(Noiseless Energy - Statevector Energy) / Hartree\", fontsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this use case, the plot shows that the derandomized and adaptive single-qubit Pauli classical shadows seem to be more accurate than the other methods. However, it is not clear whether randomized classical shadows are more appealing than qubitwise commutativity. From those results, the scaling advantage is not apparent at first sight, at least for this system (H2 3-21G).\n",
    "\n",
    "Classical shadows may show benefits for larger problems, requiring many more qubit terms to be evaluated. Several studies have taken a look at the noise implication of predicting observables with shadows ([arXiv:2011.09636](https://arxiv.org/abs/2011.09636) and [arXiv:2011.11580](https://arxiv.org/abs/2011.11580)). The authors demonstrate noise resilience inherent to classical shadows, especially when comparing to simpler sampling methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Closing words\n",
    "We demonstrated the use of several classical shadow protocols in this notebook. As their theoretical scaling is logarithmic with the number of linear functions in a qubit operator, they show potential for helping mitigate the measurement problem. Using Tangelo, we showed how straightforward it is to incorporate classical shadows in your own research projects. \n",
    "\n",
    "Many avenues have yet to be explored. For example, active research areas include implementation of classical shadows on NISQ devices, and some works report an inherent noise resilience ([arXiv:2011.11580](https://arxiv.org/abs/2011.11580)). What will you do with classical shadows and Tangelo ?\n",
    "\n",
    "## References\n",
    "1. Chen, S., Yu, W., Zeng, P. & Flammia, S. T. Robust shadow estimation. arXiv:2011.09636 [quant-ph] (2020).\n",
    "2. Koh, D. E. & Grewal, S. Classical Shadows with Noise. arXiv:2011.11580 [math-ph, physics:quant-ph] (2020).\n",
    "3. Huang, H.-Y., Kueng, R. & Preskill, J. Efficient estimation of Pauli observables by derandomization. arXiv:2103.07510 [quant-ph] (2021).\n",
    "4. Hadfield, C. Adaptive Pauli Shadows for Energy Estimation. arXiv:2105.12207 [quant-ph] (2021).\n",
    "5. Huang, H. Y., Kueng, R. & Preskill, J. Predicting many properties of a quantum system from very few measurements. Nature Physics 16, 1050–1057 (2020)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.0"
  },
  "vscode": {
   "interpreter": {
    "hash": "95050af2697fca56ed7491a4fb0b04c1282c0de0a7e0a7cacd318a8297b0b1d8"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}