ibmq-qsim-challenge.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "852ab4ac-e501-4694-bf67-90cd4122140e",
   "metadata": {},
   "source": [
    "# 1 Introduction\n",
    "\n",
    "### 1-1 Quantum Simulation\n",
    "One of the leading uses for quantum computers will be to simulate quantum systems such as molecules or engineered materials. Actually executing a quantum simulation on a current quantum computer, however, can be difficult and error prone. Your objective for this Open Science Prize is to significantly improve the state-of-the-art for quantum simulation (details in Section 2).\n",
    "\n",
    "The quantum system you will simulate is a quantum spin-1/2 model, where each qubit represents a quantum [spin-1/2](https://en.wikipedia.org/wiki/Spin-1/2) particle in a 1D chain. Quantum spin models have some amazing properties and uses. Computationally speaking, [certain optimization problems can be mapped](https://arxiv.org/pdf/1712.05771.pdf) to spin models and hence minimized on a quantum computer. Physically speaking, spin models show a variety of quantum behavior such as large entangled states, quantum phases of matter, quantum many-body affects, and many other open questions in physics research.\n",
    "\n",
    "There are many different spin models, but we will focus on the [$XXX$ Heisenberg spin model](https://en.wikipedia.org/wiki/Quantum_Heisenberg_model#XXX_model). Though it is not as well known as the [transverse-field Ising model](https://en.wikipedia.org/wiki/Transverse-field_Ising_model), the $XXX$ Heisenberg model is of practical and theoretical interest. It can also be difficult to simulate, which is where you come in.\n",
    "\n",
    "The notebook is organized as follows:\n",
    "1) Section 1 (labeled 1-X) provides you with brief introductions to the $XXX$ Heisenberg model, time evolving a Hamiltonian, and decomposing a time evolution operator into gates for a quantum computer to execute.\n",
    "2) Section 2 (labeled 2-X) outlines the specific requirements and boundaries associated with completing and competing in the Open Science Prize. Continuing from Section 1, the notebook showns an example solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d1548b05-e691-40a5-8b35-c0eaa87a2a72",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams.update({'font.size': 16})  # enlarge matplotlib fonts\n",
    "\n",
    "# Import qubit states Zero (|0>) and One (|1>), and Pauli operators (X, Y, Z)\n",
    "from qiskit.opflow import Zero, One, I, X, Y, Z\n",
    "\n",
    "# Suppress warnings\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c764af20-017f-4209-9483-4f89f10cb562",
   "metadata": {},
   "source": [
    "### 1-2 The $XXX$ Heisenberg Spin Model\n",
    "\n",
    "To simulate a quantum system, we first need to define how the quantum particles that make up the system interact with each other. In quantum physics, this is done by defining a [Hamiltonian operator](https://en.wikipedia.org/wiki/Hamiltonian_(quantum_mechanics)). For this Open Science Prize, you will work with the $XXX$ Heisenberg Hamiltonian $H_{\\text{Heis}}$ as defined below\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "H_{\\text{Heis}} = \\sum_{\\langle ij \\rangle}^{N} J \\left(\\sigma_x^{(i)}\\sigma_x^{(j)} + \\sigma_y^{(i)}\\sigma_y^{(j)} + \\sigma_z^{(i)}\\sigma_z^{(j)}\\right).\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "$N$ is the number of spin-1/2 particles in model. The operators $\\sigma_x$, $\\sigma_y$, and $\\sigma_z$ are [Pauli operators](https://en.wikipedia.org/wiki/Pauli_matrices) where the $i$ and $j$ superscripts label which qubit they act on. For example, $\\sigma_x^{(1)}$ would be the $\\sigma_x$ operator acting on only qubit 1 (which is the 2nd qubit since indexing starts at 0). This version of the general Heisenberg spin model is called $XXX$ because the same $J$ value multiplies each pair of Pauli operators. The sum notation $\\langle ij \\rangle$ means the sum is over nearest neighbors (only qubits next to each other interact), and $J$ is the interaction strength, which we will set $J=1$.\n",
    "\n",
    "You will work with the explicit case of $N=3$ with the 3 spins arranged in a line. Written out fully, the Hamiltonian is\n",
    "$$\n",
    "H_{\\text{Heis3}} = \\sigma_x^{(0)}\\sigma_x^{(1)} + \\sigma_x^{(1)}\\sigma_x^{(2)} + \\sigma_y^{(0)}\\sigma_y^{(1)} + \\sigma_y^{(1)}\\sigma_y^{(2)} + \\sigma_z^{(0)}\\sigma_z^{(1)} + \\sigma_z^{(1)}\\sigma_z^{(2)}.\n",
    "$$\n",
    "Now that we have a Hamiltonian ($H_{\\text{Heis3}}$), we can use it to determine how the quantum system of 3 spin-1/2 particles changes in time.\n",
    "\n",
    "*See the supplementary material for further details on $H_{\\text{Heis}}$, $H_{\\text{Heis3}}$ and its matrix repsentation.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e6d371fd-f5c4-446a-83f4-ac7593fdcee1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Returns the matrix representation of the XXX Heisenberg model for 3 spin-1/2 particles in a line\n",
    "def H_heis3():\n",
    "    # Interactions (I is the identity matrix; X, Y, and Z are Pauli matricies; ^ is a tensor product)\n",
    "    XXs = (I^X^X) + (X^X^I)\n",
    "    YYs = (I^Y^Y) + (Y^Y^I)\n",
    "    ZZs = (I^Z^Z) + (Z^Z^I)\n",
    "    \n",
    "    # Sum interactions\n",
    "    H = XXs + YYs + ZZs\n",
    "    \n",
    "    # Return Hamiltonian\n",
    "    return H"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35323a1f-9aac-493a-87a8-116c3b7980ca",
   "metadata": {},
   "source": [
    "### 1-3 Time Evolution\n",
    "\n",
    "Knowing the Hamiltonian, we can determine how quantum states of that system evolve in time by solving the [Schrödinger equation](https://en.wikipedia.org/wiki/Schrödinger_equation)\n",
    "$$\n",
    "i\\hbar \\dfrac{d}{dt}|\\psi(t)\\rangle = H |\\psi(t)\\rangle\n",
    "$$\n",
    "\n",
    "For simplicity, let's set $\\hbar = 1$. We know that the Hamiltonian $H_{\\text{heis3}}$ does not change in time, so the solution to the Schrödinger equation is an exponential of the Hamiltonian operator\n",
    "$$\n",
    "\\begin{align}\n",
    "U_{\\text{Heis3}}(t) &= e^{-it H_\\text{Heis3}} = \\exp\\left(-it H_\\text{Heis3}\\right) \\\\\n",
    "U_{\\text{Heis3}}(t) &= \\exp\\left[-it \\sum_{\\langle ij \\rangle}^{N=3} \\left(\\sigma_x^{(i)}\\sigma_x^{(j)} + \\sigma_y^{(i)}\\sigma_y^{(j)} + \\sigma_z^{(i)}\\sigma_z^{(j)}\\right) \\right] \\\\\n",
    "U_{\\text{Heis3}}(t) &= \\exp\\left[-it \\left(\\sigma_x^{(0)}\\sigma_x^{(1)} + \\sigma_x^{(1)}\\sigma_x^{(2)} + \\sigma_y^{(0)}\\sigma_y^{(1)} + \\sigma_y^{(1)}\\sigma_y^{(2)} + \\sigma_z^{(0)}\\sigma_z^{(1)} + \\sigma_z^{(1)}\\sigma_z^{(2)}\\right) \\right]\n",
    "\\end{align}\n",
    "$$.\n",
    "\n",
    "Now that we have the time evolution operator $U_{\\text{Heis3}}(t)$, we can simulate changes in a state of the system ($|\\psi(t)\\rangle$) over time $|\\psi(t)\\rangle = U_{\\text{Heis3}}(t)|\\psi(t=0)\\rangle$. \n",
    "\n",
    "*See the supplementary material for more information on exponentiating an operator or matrix.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ccdedbc8-f335-468c-9db5-b594aa2e50bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Returns the matrix representation of U_heis3(t) for a given time t assuming an XXX Heisenberg Hamiltonian for 3 spins-1/2 particles in a line\n",
    "def U_heis3(t):\n",
    "    # Compute XXX Hamiltonian for 3 spins in a line\n",
    "    H = H_heis3()\n",
    "    \n",
    "    # Return the exponential of -i multipled by time t multipled by the 3 spin XXX Heisenberg Hamilonian \n",
    "    return (t * H).exp_i()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1668a459-1c19-4769-9263-ff0d8de330af",
   "metadata": {},
   "source": [
    "### 1-4 Classical Simulation of $H_{\\text{Heis3}}$\n",
    "\n",
    "Now that we have the operator $U_{\\text{Heis3}}(t)$ which describes how any state of the quantum system changes in time, we can simulate the $XXX$ Heisenberg spin model! Just for fun, let's simulate the system using a classical computer first.\n",
    "\n",
    "(Fun fact: Notice the size of $U_{\\text{Heis3}}(t)$ in matrix form. It's representated by an $8$x$8$ matrix. This is because there are $2^3=8$ states in the $N=3$ system. If the simulation were of $50$ particles ($N=50$), $U_{\\text{Heis}}(t)$ would be approximately $10^{15}$x$10^{15}$! Well beyond the capacity of today's computers.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0813fede-9945-4477-9304-9f9891ce50a7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEnCAYAAACdTPWPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABHsElEQVR4nO3deXhU5fXA8e/JAiEsIWGJ7AECCMgmoKwSFNRWC1q3auvaarW2aqu1trZuta2trbW21qrtr7YudVdoXQGJLAqyyCJ72AkQCIQlJGQ9vz/ujY5hksxMJnNnMufzPPNMcteTuZM58y73fUVVMcYYY5pagtcBGGOMiQ+WcIwxxkSEJRxjjDERYQnHGGNMRFjCMcYYExGWcIwxxkSEJRxjjDERYQnHGGNMRFjC8ZiI5IiIish9sXyOpiQiLUXk1yKyWUQq3L9luNdxeUFEnhGRqL9bW0QudK/TOV7HEigRucaN+RqvY2mu4jLhiEiW+8aq75HrdZyBivWEEoAfAz8FtgO/Be4H9jb2oOF83aLhGojIt0TkSRFZJiJlbjwXhGt7d5+xIvKeiBwWkWIRmS8i5/nZ9D2gFJjeqD8qzojIRhGpEpHWdazPdK/T8kjHFg5JXgfgsQ3Ai3Ws2xbBOJraJ8BAoNDrQEL0FaAYOEdVK7wOJoo9CPQC9gMFQI9wbi8iU4B3gBLgP+7zpcD/RORaVX2mZltVLRGR2cA0EblZbQytBolIWyAbWKuqx+rY7DT3+dPIRBVe8Z5w1qvqfV4H0dRUtQRY73UcjdAFOGDJpkHfATaq6g63pHVvuLYXkWTgaaASmKCqq93lv8b58PuTiMxU1YM+u80AvgaMBJaG9ifFlRGAAEvq2aYm4cRkCScuq9SCISKT3CLsX+pYP85d/6TPsmQRuUNEVotIqYgUici7IjIxwHPWWZdce537QTHXXX2vb7Wgzz5+q3uCidP3GCIyWkRmu1UqB0XkeRHpFMjfFsx53XMp0BvoFUx1p4hcKiILRKTQPcd2EXlTRCYE8rqJSAsRuUVEZolIvoiUi8geEXlBRPrVjrO+Y/lsd5GI5LrVUaUi8qmI3BDo69YQVZ2tqjuaaPuzgCzguZpk4x6jEPg90A64pNY+/wWqCaJaLZj3vrss6PeliLQWkT+IyG73OiwXkYsCiK3B61crnjNE5AMROSIiWwP48091n+tLOKPd55hMOPFewgnEPGAHcKmI3KaqlbXWX+E+Pw8gIgK8CkwD1gJ/BtKBy4C5InK5qr4SxvhycT4IrgY+dH9vUCPiHA3cCcwG/gaMx3kN+ojIuIaqToI8b83fcpv7/Kj7vK2Bc9wM/AXYjFNlWgx0BSYBOcACGn7dMoBHcK7/f4HDwACcKqRzRWSkqtZ8iDR0LETkYeAOnHaol3Cqo6YCT4rIQFX9YX1/UxSY5D7P8rPuPZ9tPv/ipar7RGQRTsL5RdOGF9j7UkQSgP/hvA+WA//GeW+8gP+/DXe/YK/feOBnwPvAE0DbAP6GQBNONbAqgONFH1WNuwfOh4PiVDPdV8djjM/2D7nbf7XWcZKAfThvQnGXXe1u+x6Q5LPtQOAYcAho67M8x93+Pp9l17jLrvET+wnr/B2j1j7+zhFqnApc5LM8AZjjLh8bwGsf1HnddduAbUFc3+XALiC11nIBMgJ53YCWQFc/yyfhVCv9vaHX2GfdOe66N4EUn+XJ7jIFRgf4tz0DaADb3ece94IAj1vv9jhfEhQ41c+6ZJwPwU/8rLvT3a93gHGE+t4P6H0JXOdzLRJ8lp/lc5xrap034OtXK54rAn3Puvt+5u73MvCcn8fL7vq1wRw3mh7xXqU2AKfe2t9jjM92z7nPV/BlZwOdgBfUfcfgfKAC/Fh9SkOqug6nDjwNuCB8f0LIQo3zQ1V9zWf7apxviQCjmvC8wSrHSQyfU8fBOrb/ElUtU9XdfpZ/iFMymxJELDfjfFDcoKrHfY5VAfzc/fWyII7nhXbu85HaK9y/oxTn2tU2w31u6t5qgb4vv+U+/8zdpmb7OdRdwgnl+i1V1RcCDV5EWgEnu79eAnzTz6OmyrLJqtNE5BUReaipjh/vVWozVPWChjZS1c9EZBVwgYikqtMID86bAL5ISADDgSJV9VfkzQVudbd5NsSYw2U4ocXpr3dMvvvcvgnPG4yXcEqln4nISzjVXB9r3T1//BKRkTjf0McDnXG+0dYoD+JQp+N8UH/PqVH8kppjnlx7RZSpCTyo3maqukFECnBKho+GOygfgb4vhwEHVXWtn+0X4FST1RbK9Qu2k8QwIBF4RlWv9beBiNwP3IOfhCMia4DXVPWeWstnAsWqWvvLcl1uxKkurJeI3IVT4uyJ82VjIc6XyA317RfvCScYz+PcAzId+I84/eSnAytUdY3Pdu2ATXUcY6/PNl4LNc7DfpbVlCQSm/C8wfgdUATchPMN9OfAcRF5EfiRqhY1dAC3c8EcnKqi94A8nCo/xflH6xVEPBk4/2v31rON3/suokjNdT+hFCNOD7ZW+HlviEgvIBNY3KTRBf6+TKPu919BHctDuX776tnWn5r2m5X1bDPcff5ScnVLRwPwX/I5FfhToEGo6oEAN90MfN99bodTJfuuiPTxqe05QbxXqQXjBZwPm5pSzQU4b7Lna213BOcfzJ9Mn23qU1PU9/cBHq5kFY44o/K8btXZU6o6AjgJp7rjQ5xE8UyAh/kp0AI4U1UvUNU7VPVedbrRlwYZ0hEgX1WlnsfkII8ZaTUf0tl+1vXFKQHl+Vl3gfs8w886f5r6vX8Yp7TqT13vy1CuX1AlQRqRcIAhOK/XlxKO20Ovm+/2ItJfRGa6Pfn2iMiv3I48iMhwcW46TXV/T3bXbxWR4yKyQ0TuAFDVV9Tp5bhVVVfilLyycDpg1MkSToBUdRdOj6WzRaQDTuKpxrkBztcKIF1ETvFzmEk+29TnkPvczc+6EX6WVbnPgZQwaqyg8XGGIqLnVdUCVX0Z+CrOh+ZXRKSmZF/f69YX596fj30Xikimu662+o71CdBNRBq6ETOafeg++6tyqhm+Zp6fddOBTW4bXSAOuc+BvveDtRLIEJFBftZNqGOfSFy/ehOOiKTjVF9tUdVDtVaPwHmv1u7iXnPMT91jnAJ8jFN1PRznM+xa4Ep3u6HAZp8mg/uBM93tTga+DZxwHd0S1jU4XzjqHQHEEk5wnsOps70Z5x9vrqrm19qmpqHyNyLy+YePiPQHbsD5htXQt73lON+QLhORlj7HOI0vSli+ahrC/f2T1iUccYaiyc8rImf7HtuVCrTBaXup+RZd3+u2A+eDaaDPcVvgdLdO9rN9fcf6s/v8DxHxVyXVW0Sy/OwXTebg9Bb8pu+XBRHpCNyOUwr4Ujd690NyIk5PrkAF+94PVk1766/cLtI1xz8L/8kUmvj6uVWSg4Gd9XRqGe4++2urGgF0kFrDcwHvusesqSb7G04vykdUNU9VP8Cpoclx1w/ly92tp+C0C32kqttUdZaqvuUT9/kiUoxT1XwecLaqVlGPeG/DOVnqHvvqkKo+WmvZqzgfOD/Hee1qV6eB84F6MXA+8KmIvMMX95m0Ar6pqvVWGanqLhF5zT3OEhGZhfNBNh14ixN7cW0A9gDfEJFjuI2lqlpfb5NGxxmiSJz3ZaBYRBbgdFlPxfmH6AI86NM7qb7X7S84H0AL3Y4HlTj/gMk430KH1TpnncdS1bdF5Dc41XR5IvIeTrftTjjdwcfg9IDc1pg/WkS+wxff0oe7z7fIF+Oj/V1VF4SyvapWiHOT49s4r4nv0DbdgGv91P+fh/N/EvAXiBDe+8F6Bucb/QU+x++K8/57y425dkxNff2G4FTfrqhnm+Hus792mhHAP4HHai1/CDgOIM7NyuOBEe7fUqMFX1QzD8XpOFHjf8BDIjIO53/qTd9eejg3Ow/Hqba+E3hZRCaoalmdf0Ugfaeb24Mv7sOp77Gtjn1fd9eXAu3q2CYZ+Amwxr3gh3Aanif52TYHP/dv4LQPPYHT+FiK0+h6LnXcpwCMA+bj3OSo+NyrUc85Gh1nQ+sa+/q422+r63rUsf1NODdrbnePvw+nuucyP9vW97pdivONsgSnquAZnHr+XN/tAjmWu/4rOB9qhTglrXycqqrbgY4B/m3P+Du377p6HrXfM0Ft7/M3vo9TojmG0zvp/DriecV97RMC+dtCee+H8r7EKek+gvMFodS9xhf5O36w16++eOr5e7/t7vPLhq47cG6t5Ynu+/Obfvb5DLjX/fkiYCtOG1ztR4a7zV5q3YeFkwzvxSnxrwNa1BFfC/f98PX6/taamxWNMTFARJ4BrlbVE/rnRhO3Omw/8IqqftvreJorty1qDTBYfbp6u+0qR3ESwEwRmQb8C+ikJ46WUtPBYB/QV1W31HOek1T1hN58bnXzIeAbqjqzrnjjvUrNGNM0zsQZzqUp2gHNF0bglNJq3/9S03Otps3nY5yehE+KyB9w2jEH44wA8XucKuJinFIQIvITYCfO/UQtgbuBRTXJRkR+h1PbsxunSu0u4AD+O458zhKOMaYpfIjTZrbf60CauRHAKj2xsb6m59pOAFXdLyJfw7mX8BOcKsE1OPesgZOgVusXVV4pOL3UeuDc0/YeX24/6w28BnTEucbzcW4hOFRfsFalZkwMiZUqNWP8sRKOMbHlTZrX5IAmjlgJB+jYsaNmZWWFtO+xY8do3TraRyVp3uwaeM+ugbe8eP2XLVtWqKoBz4MFVsIBICsri6VLQ5uQMDc3l5ycnPAGZIJi18B7dg285cXrLyLbg93HRhowxhgTEZZwjDHGRIQlHGOMMRERNQlHRLqLyJ9F5GMRKXEHoMsKcN8UEXnYHW671D3GGU0csjHGmCBETcLBGdPnUpybjOYHue8/gOtx5mQ4H2eMpPdEZHg4AzTGGBO6aOqlNk9VM+HzUWzPDmQnERmGM1Lrdar6T3fZhzh30T4ATGuacI0xxgQjako4+sWQ8cGaBlTgzGNfc6xK4EXgHN85NYwxxngnahJOIwwGtuoXs9TVWIMzZLa/KXHD4k+zN/HJnkoOl1Y01SmMMSbsio6VM2NFPk/kbo7oeaOpSi1UGTjtPrUd9Fl/AncyqRsAMjMzyc3NDeqkZVXK3z8s4WgF/G3V+/Rrn8BXeiczvHNzeEljS3FxcdDXz4SXXQNvBfr6L9lbyXvbKth8qBoFMlKE/tU7SEyIzNB8zeHTUXAmJvK3vE6q+hTwFMCoUaM0lLt0z8qp5p8z53IotRtvr97Lnz49xt1f7cO3J/RGxMZWjBS7y917dg281dDrr6o8NiePx1dspF/nNvzgrC5MHtCJod3bRyzZQPNIOAeBnn6Wp/usbxJJiQn0S08kJ+dkfnBmP3740goefGsdOw6WcM/5g0hKbA41lsaYWFZeWc3P3ljNq8t28fUR3XjooqG0SPLmsymks4pIOxH5hoj8R0SeE5GLRaRNuIML0Bqgt4ik1lo+CGfOh7xIBJGSnMjjV5zKDWf04d8fb+fO11ZF4rTGGFMnVeUH/1nOq8t2cduUfvzh0mGeJRsIIuG4N2Z+T0Tew5mK9El3/xTgn0ChiLwjIt8VkS5NE65fM4Fk4BKfWJOAy4D3VbUsUoEkJAg/++pAbjmrH68vz2fGivxIndoYY07w/OIdvLemgLu/OpDbpvT3vKq/wYQjIp1EZBmwHWea0S04M791UtXLVPVinFnfLgJ2APcCu0RkkYikBROMW1K6GBjpLvqKu2ySu76XiFSKyD01+6jqCpwu0Y+KyHdE5CycLtG93Vgi7taz+jGqVzo/f/MzdhXV7jxnjDFNb/P+Yh58ay1n9O/Etyf09jocILASTiLwDjBGVbup6k2q+q6qltdsoKplqvqWqn5XVbsCE4C5AR7f1yvu40b397+6v9/v/i5uPLWPey1OKetB4C2caVHPVdXlQZ4/LBIThD9eNhxV+NHLK6mqtjmHjDGRU15ZzW0vrqBVciIPXzyUhAh2DKhPg50GVHUv8PNgDqqqHwMfBxtMQ9Pmquo2/PQ+U9VS4EfuIyr0yEjlvmmDueOVlTw9fws3TurrdUjGmDjx2JxNrM4/zN++dSqZ7VK8Dudz1o2qCV10ajfOHpTJY3M2caA4Yk1Jxpg4ln+olKfmbeHCEd0495RINqc3LJhOA61E5DYRmSsiBSJS7j4K3GW3+ekpFtdEhDvPPZnjFVX87cPI3tFrjIlPf/lgEwB3nDPA40hOFFDCEZEewCrgYZwqrVeB3wK/c3/G/XmliPi7JyZuZXduw4UjuvPvj7dTcOS41+EYY5qxbYXHeHnpLq44vSfd2rfyOpwTBHrj56NAKdDPbUc5gTt3zZvAH3F6rBnXrWf1Y8aKfB6fm8cD00/xOhxjTDP12JxNJCcK38uJzjbjQKvUpgB315Vs4PMG/XvcbY2Pnh1SuXR0D/7zyQ7rJm2MaRJ5+47yxop8rh6bReco6ijgK9CEE0y/XusD7McPzsxGRPjznIgMfGCMiTN/nL2J1OREvhvFPWIDTTizgV+JSJ13D7lVar8EZoUhrmanS1orvjG6B298mk+h9VgzxoRRYWk176zew1Xjssho3cLrcOoUaMK5DWgFbBSR+SLyhIj8WkR+5f48D9jobvPDJoo15l01NovyqmpeWrLT61CMMc1I7s5KAK4c08vjSOoXUMJR1V3AUOAOoAxnaJvb3d8vxJlx88fAcHdb40d25zaM69uBFxbvsNEHjDFhUVZZxYe7KpgyMJOuUdgzzVfA9+Goaqmq/klVp6hqF1Vt6T5OUtWz3HXWIt6AK8f0Iv9QKXPX7/M6FGNMM/DuZ3s5Wg5Xjo3u0g3YSAMRN2VQJpntWvLsou1eh2KMaQae/Xg7manC+L4dvQ6lQWFNOCJyhoh8EM5jNjfJiQlcflpPPty4n+0HjnkdjjEmhq3dfYSl24s4s2dy1AzQWZ9wl3A6AZPCfMxm5/LTepKYIDy/eIfXoRhjYthzi7eTkpzAhG6xMXlzQFEGMVxNp0bEEjcy26Vw9qBMXl66kzvOHuDpDHzGmNhUUl7Jm5/m87WhXWmdXOR1OAEJNC1uI7AbOiXA7eLexSO7885ne5m/aT9nDcz0OhxjTIyZtbaAkvIqLh7ZndIdzSvhlALz+GKgzrqMAm5oVERx4oz+nUhPTWbGit2WcIwxQZuxYjdd01IYnZXBvBipnQ804awEqlT1H/VtJCKHsIQTkOTEBL46pAuvL8/nWFklrVvGRh2sMcZ7B4+VM2/jfr49sXdMdBaoEWjjwTJgZIDbxs5f77Hpw7tRWlHFrLUFXodijIkhb63eQ2W1Mn1YN69DCUqgCech4BsNbaSqr6mqtYAHaFSvdLq1b8WMFfleh2KMiSEzV+TTP7MNA7u09TqUoAQ6tE2+qn7Y1MHEm4QE4WvDujJvU6FNQW2MCciuohKWbCti+vBuiMRWhZKVRjx2wYiuVFUrb6/e43UoxpgYMHPlbgCmDevqcSTBs5EGPHbySe0YkNmWGSt2ex2KMSYGzFyxm5G90umRkep1KEGzkQaiwLThXVm6vYg9h0u9DsUYE8Xy9hWzfu/RmCzdQIAJR0R6BvLARhoIyTmDTwKw3mrGmHq9t2YvAGcPjs1792ykgSiQ3bkNfTu15v01BVw1NsvrcIwxUer9tQUM655Gl7TonvemLjbSQJQ4e/BJPD1vC4dLKkhLTfY6HGNMlNl7+Dgrdx7ix+cM8DqUkNlIA1HinMEn8UTuZj7YUMCFI7p7HY4xJsrMWutUp50To9VpYCMNRI2h3dLIbNeS99dYO44x5kTvry2gT6fWZHeOrZs9fdlIA1EiIUE4e9BJ5G7Yz/GKKq/DMcZEkcMlFXy8+QBnDzrJ61AaxUYaiCJnD86ktKKKBZsKvQ7FGBNF5m7YR2W1xmzvtBpWGokiY/p0oG1K0uddH40xBpzu0J3btmR49/Zeh9IolnCiSHJiAmed3JnZ6wqorKr2OhxjTBQ4XlHFhxv3M3VQZkxNReBPoFNMX9WYk6jqvwM4Rw/gj8BUnI4Hs4HbVLXBqYXcm05/CUwGOgK7gJeB36jqsUaEHnFTB53Emyt2s3LXIUb2yvA6HGOMxxZtOUBJeRVTBsV2dRoE3i36mUacQ4F6E46IpAIfAGXA1e4+DwJzRWRofUlDRFrjJKdk4BfADmA0cD/QD7isEbFH3IR+HUlMEOau328JxxhD7ob9pCQnMLZPB69DabSAEk4Eep5dD/QBBqhqHoCIrAI2Ad8FHqln3/E4ieUcVX3fXTZXRDKAO0QkVVVLmi708EprlczIXunM3bCPO2L4Bi9jTOOpKh+s38e4vh1JSU70OpxGi5Y2nGnAoppkA6CqW4GFwPQG9m3hPh+ptfwQzt8Xc5Wekwd0Zs3uIxQcOe51KMYYD20tPMaOgyVMHtA8hqmMloQzGPjMz/I1wKAG9p2NUxL6rYgMEpE2InImcCvwt1hrwwGYfLLz5vpww36PIzHGeGmu+xmQM6Czx5GER6BtOE0tAyjys/wgkF7fjqp6XEQmAK/hJKgafwe+X9d+InID7jA8mZmZ5ObmBhmyo7i4OOR966KqZKQILy9YQ+djm8N67OaoKa6BCY5dg6bx+pJSurYWNq/6hPo+CWLl9Y+WhAP+R5lusDpMRFKAl4DOwJU4nQZOA+4BKoGb/J5M9SngKYBRo0ZpTk5OSEHn5uYS6r71OadoNf9duZvxE88gOTFaCqLRqamugQmcXYPwO1ZWyaZZs7h6XBY5OfVX9MTK6x8tn2RFOKWc2tLxX/Lx9W0gB/iqqj6nqvNU9ffA7cCNIjIsrJFGyOQBnSguq2Tptob+fGNMc/TR5gOUV1UzuZlUp0H0JJw1OO04tQ0C1jaw7xCgSFVrlzg/cZ8HNjI2T4zP7khyopC7YZ/XoRhjPDB3wz5at0hkVFbzuT0ipIQjIiNE5HURKRSRShE51V3+axE5N4RDzgTGiEgfn3Nk4XR5ntnAvnuBdBHJrrX8dPc5P4R4PNe6ZRKn9+7AXEs4xsQdVSV3/T4m9OtIi6RoKRc0XtB/idtA/zFwMvBCrWNUAzeGEMfTOLOKzhCR6SIyDZgB7ASe9Dl3LzfB3eOz7zPAUeBtEblaRCaLyI+B3+NMq7AwhHiiQs6ATmwsKCb/UKnXoRhjImhjQTG7Dx9vNr3TaoSSOh8C3sOpAvtRrXXLgVODPaDbdflMYCPwLPA8sBU4U1WLfTYVIBGfuFV1GzAGWIEzOsHbODeSPgVMVdWYHZRsUn+ne/T8jdY92ph4Mn+T8z9f8xnQXITSS+1U4OuqqiJSu2dZIRDSK+SOmXZRA9tsw0/PNVVdC1waynmjWXbnNpzULoX5mwr5xmk9vQ7HGBMh8zYV0rdTa7q2b+V1KGEVSgnnOJBax7ouwOHQwzG+RISJ/TqyIK+Qqmp/vcaNMc3N8YoqFm85wMR+zat0A6ElnAXAbSLiO7BPzafht3EG4TRhMrF/Jw6XVrA63/K4MfFg6bYiyiqrOaN/R69DCbtQEs4vcKrVVro/K3C1iMzFaUu5P3zhmQnZHRGxdhxj4sX8TftJThRO7x37o0PXFnTCUdWVwBlAAXA3TptKzRAyk1R1Q/jCMxmtW3BK1zTm27TTxsSFeZsKGdkrndYto2kgmPAIqYO3qi5X1bOAtkB3oJ2qTlbVT8ManQFgYr+OLN9RRHFZpdehGGOa0P6jZazbc6RZtt9AI0caUNXjqro7luabiUUT+nWkslpZtPmA16EYY5rQwjynJuMMSzgOEdlS1/hkInKKiGxpfFjG18he6bRKTvy8b74xpnmat2k/6anJDO7azutQmkQoJZwsoGUd61KAXiFHY/xqmZTImD4Z1o5jTDOmqszfVMiEfp1ISIi5eSMDEmqVWl03hYzCmWnThNnEfp3YUniMnQet9tKY5mhDwVH2Hy1jYr/m1x26RkAJR0R+KCI7RGQHTrL5b83vPo/9wOPAu00ZcLyqeRPW1PEaY5qXBW4NxoTs5ptwAu13twWY4/58NbAUqN2gUIYzlcDfwxOa8ZXduQ2d27Zk4eYDNsyNMc3QwrxC+jTD4Wx8BZRwVHUGzujNiAjAA6q6tQnjMrWICOOzOzJv436qq7XZ1vEaE4/KK6tZvPUgF53a3etQmlQoN35ea8nGG+OzO3LgWDnr9x71OhRjTBit2HmIkvIqxjfj6jQIbbRoANyu0QNweqZ9iar+uzFBGf/GZztDXXy0uZBBzbTbpDHxaGFeIQkCY/s0v+FsfAWdcESkPfAWzrhp8MV0Ab491yzhNIEuaa3o06k1C/IK+c7EPg3vYIyJCQvzChnSLY201GSvQ2lSoXSL/jXQAWc8NQEuxJk87XmczgWnhS06c4IJ2R1ZvOUg5ZUxO6+cMcZHcVklK3YeavbVaRBawjkHJ+kscn/fpaq5qnoVMBu4NVzBmRONz+5IaUUVK3Ye8joUY0wYfLL1AJXV2qy7Q9cIJeF0AbaoahXOZGxtfda9DpwXjsCMf2P6dCBBYIHdj2NMs7Bg0wFaJiVwaq90r0NpcqEknL1Ae/fn7cBYn3XZjQ3I1C+tVTJDure3G0CNaSYW5hUyOiuDlOTEhjeOcaHO+FmTZJ4F7hWRJ0XkceBh4L1wBWf8m5DdgRU7D3H0eIXXoRhjGmHf0eNsKDgaF+03EFrCuZ8vhq95GGc4m/OAy4GZwA/CE5qpy/jsjlRVK59sPeh1KMaYRvjYnXKk5paH5i6UGz83q+p89+cKVb1dVburaoaqXqGqNmlLEzu1ZzotkxJYmGcvtTGxbGFeIWmtkhncNc3rUCLC5sOJQSnJiYzKSuejzdaOY0ysUlUW5h1gTJ8MEuNkqCqbDydGjevbkfV7j1JYXOZ1KMaYEOw4WEL+odK4ab8Bmw8nZtW8ST+yaaeNiUk1VeLj+sZPwgloaBsR+SHwQ/fXmvlwymtt1grIAF4MX3imLkO6pdE2JYmP8gqZNqyr1+EYY4K0cHMhme1a0rdTa69DiRibDydGJSYIY/p0YKG14xgTc6qrlY83HyCnf6eaKV/igs2HE8PG9+3ArLUF7DxYQo+MVK/DMcYEaP3eoxw8Vs64OGq/gdDacL4N7PRdICLniMjtIjI8LFGZgNS049ioA8bElpoepvFy/02NUBLOf4D/q/lFRG4E3sG5CXSxiEwJU2ymAb7TThtjYsfCvEL6dGxNl7TmO520P6EknDHA2z6//xin3SYNZ/DOu8MQlwmAiDCubwc+3lyIal0dB40x0aSiqppPth5kXJyVbiC0hNMZyAcQkWygN/AXVT0K/BMYEkogItJDRF4VkcMickREXheRnkHsP1BEXhGRQhEpFZENItLsp0oYl92RwuJyNhTYtNPGxIKVOw9xrLyK8XHUHbpGKAnnCM4EbAA5QKGqrnJ/r8LPlNMNEZFU4APgZJxecFcC/YC5ItJgn0ERGQUsxrkh9TvAV4E/AM1++NUv2nGsWs2YWLAw7wAizlQj8SboKaaBj4C7RKQSuI0vV69lA7tCOOb1QB9ggKrmAYjIKmAT8F3gkbp2FJEE4F/AHFW90GfV3BDiiDnd2reid8fWfJRXyLcn9PY6HGNMAxZuLuSUrmmkt27hdSgRF0oJ506cGzxn4pRm7vNZdxnwcQjHnAYsqkk2AG6364XA9Ab2zQEGUU9Sau7G9e3Aoi0HqKiyaaeNiWYl5ZV8uqMoLttvILTRojepan+gk6pmq+o2n9W34iSkYA0GPvOzfA1OMqnPBPc5RUQWiUiFiOwTkcdEJC66gEzI7six8ipW7TrkdSjGmHp8svUgFVXxMZ20P6FUqQHgbxoCVV0d4uEygCI/yw8CDc27WjOuy0vAX4C7cMZ0ewDoAVzobycRuQG4ASAzM5Pc3NyggwYoLi4Oed9wqSpXBHhu1lKOZsdfMT0arkG8s2sQmBfXl5GUAKU7PiM3P3wjDMTK6x9ywmkC/vr1BnJFakppz6nqPe7PuSKSCDwkIoNUde0JJ1N9CngKYNSoUZqTkxNCyJCbm0uo+4bTkxsWkF+ZSE7O2IY3bmai5RrEM7sGgfntyvmMzkrm7LPGhPW4sfL6hzpadLgV4ZRyakvHf8nHV01Ja1at5e+7z8NDDyt2jM/uyKc7izhWVul1KMYYPw4Ul7FuzxEm9IvP6jSInoSzBqcdp7ZBOAOCNrQvnFhCqikdxUVL+vjsDlRUKZ9ss2mnjYlGNVOJjOsbnx0GIHoSzkxgjIj0qVkgIlnAeHddfd7BGan63FrLz3Gfl4Ypxqg2OiuDFkkJLNxk46oZE40W5hXSNiWJId3iYzppfwJKOCJyUEROdX/+PxEJ9w0fTwPbgBkiMl1EpuGMTr0TeNInjl4iUikiNW01NZ0XfgPcKCK/FpEpInIXcA/wL9+u1s1ZSnIio3ql27hqxkSphZsLGdunA0mJ0fI9P/IC/ctb88W00tcAncIZhKoeA84ENgLPAs8DW4EzVbXYZ1PBGT2gdtwP4HTHvhTnRtSbcAYTvT6ccUa78dkdWbfniE07bUyU2XGghJ0H42s6aX8C7aW2HbheRGqSzggRqXMIG1WdF2wgqroDuKiBbbbhp+eaOiNXPkIc3/wJTsJ5+L0NfLT5gM0CakwUWZBXMx2BJZxAPIRTtXU1TuP8X+vYTtz1zX4Ms2g0pFsa7VKSWLBpvyUcY6LIgrz9nNQuJa6mk/Yn0Bk//09E3gH644xRdguwrikDM8FLTBDGZ3dkwSZnuoJ4mrrWmGhVVa0szDvA2YMy4/5/MuAbP1V1D7BHRP4FvGVTTEenif068c5ne9m8/xjZndt4HY4xcW91/mEOl1bE9f03NUIZS+1aSzbRa6L7pp6/ab/HkRhjABa4/4vxOn6ar5D654lIFxH5vYgsEZHNIvKJiPxORE4Kd4AmOD0yUsnqkMoCux/HmKgwb1Mhp3RrR4c2LRveuJkLOuGISH9gBU47TjHwCXAMZ6ToFSLSL5wBmuBN7NeJj7ccoLwyLgZZMCZqFZdVsnx7EROyw3onScwKpYTzW5xZP/ur6mRVvVxVJ+N0KDjsrjcemtCvIyXlVSzf0dAwdMaYprR4ywEqq5UzrP0GCC3hTAZ+UWseHFR1O85kbJMbH5ZpjLF9O5CYIFatZozH5m8qJCU5gZFZDc2yEh9CSTgtgKN1rDvqrjceapeSzIge7a3jgDEem79pP6f37kDLJLs1EUJLOCuAH4jIl/YVp4P599z1xmMT+nVkVf5hio6Vex2KMXEp/1Apm/cf+7znqAkt4TwATAHWicgDInKTiNyPM03AVOD+cAZoQjOxXydUnQEDjTGRV9MdemI/6zBQI5T7cN4FzsepPrsbeBz4OU6PtfNV9f16djcRMqy7M8zNvI1WrWaMF+ZtLKRz25b0z7QbsGuENMW0m3TeFZFU3Fk5VbUkrJGZRklKTGBiv058uHG/DXNjTIRVVlUzb9N+vnLKSfa/56NREzOoaomq5luyiU6TBnSi4EgZ6/bU1cfDGNMUlu84xNHjleQM6Ox1KFElfmcCigM5/Z2647kb9nkciTHxJXfDvs8H0zVfsITTjHVul8Lgru34cIO14xgTSbkb9jOyVzpprZK9DiWqWMJp5nIGdGLZjiIOl1Z4HYoxcaHgyHHW7jlCzgDrnVabJZxmLmdAZ3c+DusebUwk1NQo5PS39pvaQhm8c6GIXOkz3bSJYiN6tKddShK51o5jTETkbtxHZruWDOzS1utQok4oJZwK4F/AbhF5RERODnNMJoxqukfnbnC6Rxtjmk5lVTXzNxUyqX8n6w7tRyg3fuYAA3GSzlXAGhHJFZHLRMRayKJQzoBO7Dtaxto9R7wOxZhmraY79GTrDu1XSG04qrpBVX8EdAOuARKBF4BdIvKQiPQJX4imsSa5jZe51lvNmCY1d8M+khKE8TZ+ml+NvfGzTFWfxZl8bT7QCbgT2Cgir9gMoNGhc9sUhnRL44P11o5jTFP6YN0+RvZKp12KVfb4E3LCEZFWInKdiHwCLMFJNrcCXYGbgHHA82GJ0jTalIGZLN9RxP6jZV6HYkyztONACRsKjjJ1UKbXoUStUHqpDRGRvwC7gb8B24EpqjpYVf+sqntV9WngRmB8eMM1oZo6KBNVmGulHGOaxKx1BQCWcOoRSglnJXAB8CjQS1UvUdW5frbLAz4OPTQTTgO7tKVb+1a8v7bA61CMaZZmrd1L/8w29OrQ2utQolYoCecSnERzv6ruqWsjVV2nqjbddJQQEaYOymRB3n5Ky6u8DseYZuVQSTlLthVZ6aYBoXSLfk1V7RMrBk0ZmMnximoW2KgDxoTV3A37qKpWpg6yflL1CaUNZ4uIDKtj3SkisqXxYZmmcHqfDNqmJDHbqtWMCavZa/fRuW1LhnZL8zqUqBZKlVoWUNewNilAr5CjMU0qOTGBnAGdmbO+gKpqG3XAmHAoq6wid8M+zhqYSUKCjS5Qn1C7Rdf1aTUKOBTiMU0ETB2USWFxOSt2HvI6FGOahUVbDnKsvIqpg2x0gYYENMW0iPwQ+KH7qwL/FZHyWpu1AjKAF8MXngm3Sf07kZQgzFpbwMhe6V6HY0zMm7V2L6ktEhnX10YXaEigJZwtwBz3IcBSn99rHq/hJKXrQwlERHqIyKsiclhEjojI6yLSM4Tj/FREVEQWhBJHc5fWKpmxfTvw7md7bDBPYxqpqlp5b00Bk/p3IiU50etwol5AJRxVnQHMAGpGQP2lqoatc4CIpAIfAGXA1TilqAeBuSIyVFWPBXicPsDdgN3dWI/zh3bhJ6+tZs3uI5xijZzGhGzJtoPsP1rG+UO7eh1KTAilW/S14Uw2ruuBPsAFqvqmm+Cm4XRA+G4Qx3kCZziddWGOr1k5e9BJJCUI/1tV521UxpgAvLVqD62SE5l8ss3uGYiAEo6IVInIae7P1e7vdT0qQ4hjGrBIVfNqFqjqVmAhMD3AGK8ATgV+GsL540p66xaMz+7IW6t3W7WaMSGqqlbe+WwPZw7sTGqLgCqL4l6gr9IDwC6fn8P9KTUYt8quljU4IxvUS0TSgT8Cd6rqQZv4qGHnDe3Cna+uYnX+YYZ2b+91OMbEnMVbD1BYXM75Q7p4HUrMCLQN536fn+9rgjgygCI/yw8CgXSlehjYCDwT6AlF5AbgBoDMzExyc3MD3fVLiouLQ97XS60rlESBJ976hMsGtPA6nEaJ1WvQnMTjNfjXmjJaJkLCvvXk5m7wNJZYef2jqRzor9TUYFFFRCbizDx6qgZRP6SqTwFPAYwaNUpzcnIC3fVLcnNzCXVfr7266xNWFRTz10mTYno63Fi+Bs1FvF2Dyqpqbp8/h7NP6co5Z43wOpyYef0DbcNpqN2msW04RTilnNrS8V/y8fUk8A+c2Ubbi0h7nESa6P5e16gIce+8oV3JP1TKyl2HvQ7FmJiyeOtBDhwr5zyrTgtKMG04Tdm6vAanHae2QcDaBvYd6D5u9LOuCOfeoEcbE1xzNXVQJsmJwv9W7mZ4j/Zeh2NMzPjfqt20bpFIzgDrnRaMQNtw7mviOGYCvxeRPjVdrkUkC2cCt7sa2NffFAiPAonAD3Dm5TF+pLVKZlL/zsxcuZu7vnIySYmNmnHcmLhwvKKKt1btYeqgTLvZM0jR8gnzNLANmCEi00VkGk6vtZ04VWYAiEgvEakUkXtqlqlqbu0Hznhuh93fd2HqdPHIbuw7Wsb8TTZlgTGBmLW2gCPHK7l4ZA+vQ4k5UZFw3JEEzsTpafYszs2bW4EzVbXYZ1PBKblERdzNwZknZ5KemsyryywvGxOIV5btomtaCuP6dvA6lJgT6OCdVcBYVf1ERKqpvz1HVTXo3m+qugO4qIFtthFAzzVVzQn2/PGqRVIC04d344XFOzhUUk771NjuIm1MU9p7+DgLNu3n5snZNhVBCKLlxk/joUtGdeeZj7Yxc+Vurhqb5XU4xkSt15bvolrh4pHdvQ4lJkXLjZ/GQ4O7pjGoSzteWbrLEo4xdVBVXl22i9N6Z9CrQ2uvw4lJIbeFiEg7ERknIpe4z+3CGZiJrItHdmd1/mHW7z3idSjGRKXlO4rYWnjMSjeNEFLCcXuJ7QTmAy8BC4CdIvLzMMZmIuiCEd1IThReXWqdB4zx55Wlu0htkWg3ezZC0AlHRO4H7sNJNFOBIcAU4GXgfhG5L4zxmQjJaN2CKQMzeW35Lo5XVHkdjjFR5XBpBTNX7uarQ7rQumU0jQgWW0Ip4VwP/EFVb1DVD1R1jft8Pc6IzTeEN0QTKVeNzaKopII3P833OhRjosorS3dSUl7FNeOyvA4lpoWScNKA9+pY96673sSgMX0yOPmktvxz4TabJ8cYV1W18sxH2xidlW4z5DZSKAlnMTC6jnWj3fUmBokI103ozYaCo3y8+YDX4RgTFWavK2BXUSnXje/tdSgxL9DRohNqHsAtwHUi8mMRyRKRVu7zncB1wPebMmDTtKYN60qH1i34v4XbvA7FmKjwz4Vb6da+FVMHZXodSswLtPWrki/f7CnAQ+6DWstXBXFcE2VSkhO54vSe/GVuHtsPHLP7DUxcW7P7MIu2HOSnNrhtWETL9AQminxrTC+eyN3MMx9t496v+Zs1wpj48MzCbbRKTuQbo3t6HUqzEC3TE5goktkuha8N68pLS3by/cnZdGhjc9iZ+JN/qJQZK3Zz6ejupKUmex1Os2BlROPXzZOzOV5RxRO5m70OxRhP/HnOJgBuysn2OJLmI+S2FhEZBgwAUmqvU9V/NyYo473szm24cER3nl20ne9M7MNJaSdcZmOarW2Fx3hl2S6uHNOLbu1beR1OsxF0whGR9sBbwJiaRe6zbxuPJZxm4LYp/ZixIp/H5+bxywtO8TocYyLmT3M2kZwofG9yX69DaVZCqVL7NdABOAMn2VyIM3na88AW4LSwRWc81SMjlctG9+DFJTvYebDE63CMiYiNBUd5c0U+V4/LonNbK9mHUygJ5xycpLPI/X2XO5XzVcBs4NZwBWe894Mz+yEi/MmtzzamufvjrI20bpHEjWdY6SbcQkk4XYAtqloFHAfa+qx7HTgvHIGZ6HBSWgpXjenFa8t38emOIq/DMaZJLcwr5J3P9vKdib1Jb22z34ZbKAlnL9De/Xk7MNZnnXXnaIZundKPzLYp/PT11ZRXVnsdjjFNorS8ip+9sZqsDqncOMlKN00hlISzgC+SzLPAvSLypIg8DjxM3QN7mhjVNiWZX15wCuv3HuWpedZN2jRPj87ZyPYDJfz660NISU70OpxmKZRu0fcDXd2fH8bpQHAZkArMBH4QntBMNJk6KJPzhnThsQ/y+MqQLvTt1MbrkIwJm8/yD/P3+Vu5bFQPxvXt6HU4zVbQJRxV3ayq892fK1T1dlXtrqoZqnqFqtoww83UvdMGkZKUwE9fW01VtY10ZJqH8spq7np9FempLfjZVwd6HU6zFvJIAyLSTkTGicglIjJWRNqFMzATfTq3TeGerw3mk20HeeiddV6HY0yjqSr3zvyMz/KP8OAFg20ImyYW0kgDInIPcDvQhi9u/DwqIg+r6oPhCs5En4tHdmf1rkM8PX8r/TLbcumoHl6HZEzInvloG//5ZCc3T+7Luad08TqcZi+UkQbuB34B/B14ESgAMoHLgftFJMkG+2zefnH+IDbvP8bdb6ymd8fWjM7K8CQOVWXnwVKW7K1kbW4eh0sqOFRSQVllFe1TW9A+NZnMdimM6Nme/p3bkpAgDR/UNFpVtbJ+7xFW7DzEviNlHCop53BpBSnJiaSlJtO+VQt6d0xlcNc0uqe3QsSb6zJv435++b+1nD0ok9unDvAkhngTSgnneuAPqvpjn2VrgA9E5DBwA3BfGGIzUSopMYHHrziVC/66kO8+u4x/X3daxKbe3VVUwtz1+5i7YT9Lth3k6PFKd80GWiYl0D41mRZJCRwqqfBZB21TkjgtK4MLT+3G1EGZtEyyXkjhVF6lvPHpLt78dDfLthdRXPbFa98uJYm01GTKKqo5VFrxpa71aa2SGZ2VwZknd2byyZ3okhaZccuWbT/IzS8sp39mW/542XD7MhIhoSScNOru+vwucFPo4ZhYkZaazP9dM5pv/X0xlz75MX++fARnDWyaGRF3FZXw35V7mLlyN+v2HAGgZ0Yq5w/tytDuaZTt2cRlX8mhVYsvJ5GKqmryi0pZtr2IpduLmLdxP99/4VPSU5P5+qndueGMPmS2s6FLGmNXUQlPzdvCK0tKKK1cSc+MVC4Y0ZVRvTIY2SudLmkpJ0xcdqyskk37ivks/zCrdx1mQV4hs9cVADCkWxrTh3fl/KFdm2zA2LdW7eGHL6+ga1oKf796FK1b2nyRkSKqwfU2EpEPgFmq+hs/634KTFXVM8MUX0SMGjVKly5dGtK+ubm55OTkhDegGLLvyHGu+9cS1u4+wn3TBnPlmF5hqSI5cryCt1ft4fXl+Xyy7SAAI3q257whXZh8cmf6dGz9+XkCvQbV1crCzYW8uGQn7322l6RE4fqJfbjhjD60TbHG4mAcLqng8dw8nvloGyiM7Czccv5oTu+dEXRpQVXZtK+YOev28fbqPazOP4wIjO3Tga+f2p2vnHJSWJKCqvLUvC385p31jOqVzlNXjSKjmYwm4MXnkIgsU9VRwewT0FUUEd+vKLcAb4hIJfAKX7ThXApcB0wPJgAT2zq3S+GlG8Zyy38+5Z4Za3h79R7u/uoghnQPvorteEUVuRv2879Vu5m1toCyymr6dGrNHWf3Z9qwbvTskNqoWBMShIn9OjGxXyd2HCjh4fc38OcP8nhh8Q7unTaYacO6NnyQOKeqvLJ0F796ex1Hjldw0and+dHU/mxcsZixfTuEdEwRoX9mW/pntuWmnL5s2V/MjBW7eXNFPne8spJfvPkZ5wzO5PyhXZnYv2NI1aHLthfxq7fWsnzHIc4f2oXfXzLMbu70QEAlHBGp5svTDwj+p5wWoFpVY6qMaiWcxquqVp5fvJ1HZ2/i4LFypg/vytdP7c7pvTPq/cfeebCEjzYXsiDvAHPX76O4rJKM1i04b0gXLhrZnWHd0xosMTXmGqzadYh7Zqxhxc5DnDe0C7+cfkqz+dYbbvuOHuenr61mzvp9nN47g/umDWZgF+duiKb4P1BVlm0v4rXlu3h79V4Ol1bQNiWJKQMzGZ/dkfHZHept8yktr2LRlgO8smwnb6/eS6e2Lbnj7P5cMrJHs2uzaVYlHOAB/CcYYwBITBCuGpvFBSO68UTuZv65cCszVuwmJTmB0VkZdElLoXXLJFolJ7L/aBk7i0rYfqCEPYePA9CxTUvOG9KF84d1YWyfDifU+zeVod3b8+qNY3ly3hYenb2RxVsO8tg3hjMu2+429/XB+gJuf3klJeVV3HP+IK4Zl9XkH9oiwqisDEZlZXD/tFNYuLmQ/63cw9wN+3jj03wAurVvRa8OqfRIT6Vj2xaUlFdRfLyS3YdLWbKtiPLKalq3SOS2Kf24fmIfa6/xWECvfiS6OYtID+CPwFScktJs4DZV3dHAfqNwesadAfQECoH5wM9VdWuTBm1O0C4lmZ+cezK3nNmPRVsP8OGG/SzacoCNBUc5VlZFSXklHdq0pEd6K07rncHwHu0Zn92Rfp3beNY9NikxgZsnZzN5QGduefFTvvWPxfzsqwP59oTensUULaqrlcfn5vHI7I0MPKkdj10+guzOkR/WqEVSApMHdGbygM5UVysbCo6yMK+QlbsOs/NgCXPWF3DwWDmtWyTRumUSGa1bcOWYXuQM6MTorPpL2SZyoiLdi0gq8AFQBlyNU5p6EJgrIkNV9Vg9u38DGAw8htM9uxvOfUJLRWS4qu5s0uCNX61aJH7+AeFLVaP2Q3xQ13a8efN4bn95BQ++tY7P8g/zm68PPaH3W7woLqvk9pdX8N6aAi4c0Y3fRMmglgkJwsAu7T6vzqsRze8t4wh1pIEuOCMNTAIygANALvCIqu4N4ZDXA32AAaqa555jFbAJ+C7wSD37/lZV99eKbyGw1T3uPSHEY5pItH8gtGmZxBPfHMlfc/P4w6yNbD9Ywj+uHh137Tr7jhznmn8uYUPBUX5x/iCuG58V9dcu2uMzIYylJiL9gRU4vdWKgU+AYzgzfa4QkX4hxDENWFSTbADc6rCFNNDrrXaycZdtB/bjlHaMCUpCgvD9M/vxxDdHsmb3ES5+4qO4mmJ78/5ivv7ER2w7cIx/XD3KqhZN2ITSMvtb4AjQX1Unq+rlqjoZ6A8cdtcHazDwmZ/la4BBwR5MRAYCnQEbYdKE7NxTTuL575zOgWPlXPjXj1iz+7DXITW55TuKuOiJjzheUcVLN4wlp1aVqDGNEcqNn4eAG1X1RT/rLgf+qqrpQR6zHKc67q5ayx8E7gqmm7WIJAFzgIE4VXR+50UWkRtwOhuQmZk58sUXT/hzAlJcXEybNjY3jJea+hrsLq7m90uPU1qp3D4yhex079sxmsLaA1X8aflx2rcUbh+VQufUwL+P2v+Bt7x4/SdPntxk3aJ9tQCO1rHuqLs+FHXd1xOsvwDjgPPqSjYAqvoU8BQ49+GE2ofd7sPxXiSuwaQJpXzr74t55NPjPH3VMMY3s27Ts9cW8Ojs5fTu2JZnv3MandsGN6yM/R94K1Ze/1Cq1FYAP6g1+gDiVPJ+z10frCKczge1pbvrAiIiv8EptVynqu+HEIcxfnVr34qXvjuGHumpXPvMEmavLfA6pLD578rd3PjcMgae1JaXvjsm6GRjTKBCSTgPAFOAdSLygIjc5E5ZsAbnHpr7QzjmGpx2nNoGAWsDOYCI3A3cBdyqqs+GEIMx9ercNoWXvjuGgSe15cbnlvHO6j1eh9Rory7bxa0vfsqpvdJ57jun0z41vnrjmcgKZYrpd4HzcKrP7gYeB36O02Pt/BBLFjOBMSLSp2aBiGQB49119RKRW3Du27lbVf8cwvmNCUj71BY8+53TGdajPd//z6fMWJHvdUghe/GTHfz41ZWM7duBZ64dbQOYmiYXVMIRkRYi8gZQ6jYWtQV6AG1V9TRVrWvagoY8DWwDZojIdBGZBswAdgJP+py/l4hUujOO1iz7BvAoztQIH4jIGJ9H0D3cjGlIu5Rk/n3daYzqlc4PX1rBq8t2eR1S0P798Tbuen01k/p34h9Xjya1RVTcA26auaASjqqW41SnJbi/l6hqvqo26iYFdySBM4GNwLPA8zg3bp6pqsU+mwqQWCvuc93l5wIf13r8tTFxGVOX1i2TeOba0xif3ZE7XlnJs4u2ex1SwJ78cDP3zFjD1EGZPHnlyKgYPcDEh1C+1iwExuCMLBA27phpFzWwzTZq9VxT1WuAa8IZizGBaNUikaevGsX3X1jOL978jOPlVVx/Rp+Gd/SIqvLH2Zt4bM4mvjasK49cOozkCA2SagyE1mngduDbIvJ9EekuIokikuD7CHeQxkSrlOREnvjWSM4f2oVfvb2OR2ZtJNh72yKhulr51VvreGzOJi4Z2Z1HLxtuycZEXCglnNXu85/cR20a4nGNiUnJiQn86RsjaJWcyGNzNlFYXMYvp59CYpTMuVJRVc2dr67ijU/zuWZcFvecP6jZzQdjYkMoicHmxjGmlsQE4XcXD6Vj25Y8kbuZg8XlPPqN4Z63jxwrq+Sm55czb+N+fnzOAL6X09fGRTOeCTrhRGJuHGNikYjwk3NPplObljzwv7Vc+Y/FPPGtkXRs09KTePYePs71/17Kmt2H+e1FQ7hsdE9P4jCmRsiVuCLSTkTGicgl7nO7hvcypvm7bkJv/nLFCFbnH2b6XxbyWX7kB/1ctr2Ir/1lAVv2F/P0VaMs2ZioEFLCce+D2Ykzs+ZLwAJgp4j8PIyxGROzzh/alVdvHEe1Khf/7aOI3SCqqry8ZCeXP7WI1BaJvHHzeM4amBmRcxvTkFDmw7kfuA8n0UwFhuDcm/MycL+I3BfG+IyJWad0S2Pm9ydwStc0bn1xBT96eQVHjlc02fmKjpXz/Rc+5c7XVjG6dzozbh5P/8y2TXY+Y4IVSqeB64E/qOqPfZatwbnL/zDO4Jn3hSE2Y2Jep7Yt+c8NY/jzB3k8PjePxVsO8vAlQxnXN7yjTc/dsI+fvLqKopJy7jx3AN89o2/U9JIzpkYoVWppQF1D2LzrrjfGuJITE/jR1P68cuNYkhOFK55ezE3PLWNr4bFGH3tjwVGue2YJ1/5zCe1Tk3nz5vF8Lyfbko2JSqGUcBYDo4HZftaNdtcbY2o5tWc6b986kb/P38qTH25m1toCLhvdg2+e3otBXQPvc6OqrM4/zLMfb+e15bto3TKJu75yMteMy/K8G7Yx9Qkl4dwCvCEilcArQAGQCVwKXAdM9x1tQFWrwxGoMc1BaoskbjmrH5ef1pM/zdnIy0t28fziHQzu2o4LR3RjRM90BnZpe8JgmsfKKlm35whLtxfxxvJ8NhQcpWVSAteM680PzswmvbVNK2CiXygJZ5X7/JD78CV8MRIB2KgDxvjVqW1LHrxgCLdPHcDMlbt5ZdlOHnxrHQAJAj0yUmnhDj1TVlnNzqISakbMGd6jPb+68BTOH9qVtFY2pYCJHTbSgDEeSm/dgqvHZXH1uCx2Hyplze4jrNl9mLx9xVS7GSYxIYGLR3bnlG7tGNw1jcx2NiOniU020oAxUaJr+1Z0bd+KqYPsvhnTPNlwscYYYyLCEo4xxpiIsIRjjDEmIizhGGOMiQhLOMYYYyLCEo4xxpiIsIRjjDEmIizhGGOMiQhRtUEDRGQ/sD3E3TsChWEMxwTProH37Bp4y4vXv5eqdgpmB0s4jSQiS1V1lNdxxDO7Bt6za+CtWHn9rUrNGGNMRFjCMcYYExGWcBrvKa8DMHYNooBdA2/FxOtvbTjGGGMiwko4xhhjIsISjjHGmIiwhBMCEekhIq+KyGEROSIir4tIT6/jihci0l1E/iwiH4tIiYioiGR5HVc8EZGLReQ1EdkuIqUiskFEfiMibb2OLV6IyDki8oGI7BWRMhHZJSIvi8ggr2Ori7XhBElEUoGVQBnwc5zpth8EUoGhqnrMw/DigojkAC8By4BE4Gygt6pu8y6q+CIii4AdwAxgFzACuA9YD4xT1WrvoosPInI5cCqwGNgP9ATuAnoAQ1Q11JvZm4wlnCCJyK3AI8AAVc1zl/UGNgF3quojXsYXD0QkoeYDTUS+AzyNJZyIEpFOqrq/1rKrgH8BZ6nqB95EFt9EZABO0r9DVf/gdTy1WZVa8KYBi2qSDYCqbgUWAtM9iyqO2Ldn79VONq4l7nO3SMZivuSA+1zhaRR1sIQTvMHAZ36WrwGitu7UmAiY5D6v8zSKOCMiiSLSQkT6AU8Ce4EXPQ7LrySvA4hBGUCRn+UHgfQIx2JMVBCRbsADwGxVXep1PHFmMTDS/TkPOFNV93kYT52shBMafw1fEvEojIkCItIGp/NAJXCtx+HEoyuBMcAVwBFgVrT22rSEE7winFJOben4L/kY02yJSAowE+gDnKOquzwOKe6o6jpVXayq/wHOAtrg9FaLOlalFrw1OO04tQ0C1kY4FmM8IyLJwGvAacAUVV3tcUhxT1UPiUgekO11LP5YCSd4M4ExItKnZoFbfB3vrjOm2RORBOB5nG/U01V1kcchGUBEMoGTgc1ex+KP3YcTJBFpjXPjZylf3Pj5S6Atzo2fxR6GFzdE5GL3x7OAG4Hv4dz8tl9VP/QssDghIk/gvO6/Av5Xa/Uuq1preiLyBrAcWIXTdtMf+CFwEnCaqm70MDy/LOGEwB3G5o/AVJzOAnOA2+zGw8gRkbreuB+qak4kY4lHIrIN6FXH6vtV9b7IRROfROQnwKVAX6AFsBPIBX4TrZ9FlnCMMcZEhLXhGGOMiQhLOMYYYyLCEo4xxpiIsIRjjDEmIizhGGOMiQhLOMYYYyLCEo4xTURELhCRH9ValuNOiZ3jTVTGeMfuwzGmiYjIMzhjjHX3WdYOd9w9VT3iVWzGeMEG7zQmgtwkY+OOmbhkVWrGNAG3dHM10M2tQlMR2eavSk1EckVkgYicKyIrRKRURD4VkdNFJElEfi0ie0TkoIg8447n53uuVBH5rYhsFZFy9/lud4BNY6KGlXCMaRq/BDoBo4Fp7rIyIK2O7bOBh3EGwywGfocz+vhMnP/Ta4CB7jb7gDsBRCQJeA+nmu6XwGqcybh+gTNv0+1h/auMaQRLOMY0AVXdLCL7gXLfofvr6SzQARinqlvc7RJwZtHsrapT3G3eE5EzgEtwEw5wOTABmKSq89xlc0QE4F4R+W20Tjds4o8VuY2JDhtrko1rvfv8Xq3t1gPdxc0owLnAduAjt/otyS31vA8k45R2jIkKVsIxJjrUnp68vJ7lSUAiUAl0xpkmoKKO43YIV4DGNJYlHGNi2wFgK868KP5si1woxtTPEo4xTacMaNXE53gXuAgoVtX1DW1sjJcs4RjTdNYCGSJyE7AUON4E53geuBano8AfcKY/b4EzC+Q04AJVLWmC8xoTNEs4xjSdv+M02v8aaI/TuH9NOE+gqhUicg5wF3AD0Bs4BmwG3uKLtiBjPGdD2xhjjIkI6xZtjDEmIizhGGOMiQhLOMYYYyLCEo4xxpiIsIRjjDEmIizhGGOMiQhLOMYYYyLCEo4xxpiI+H+p210dpKPeZQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define array of time points\n",
    "ts = np.linspace(0, np.pi, 100)\n",
    "\n",
    "# Define initial state |110>\n",
    "initial_state = One^One^Zero\n",
    "\n",
    "# Compute probability of remaining in |110> state over the array of time points\n",
    " # ~initial_state gives the bra of the initial state (<110|)\n",
    " # @ is short hand for matrix multiplication\n",
    " # U_heis3(t) is the unitary time evolution at time t\n",
    " # t needs to be wrapped with float(t) to avoid a bug\n",
    " # (...).eval() returns the inner product <110|U_heis3(t)|110>\n",
    " #  np.abs(...)**2 is the modulus squared of the innner product which is the expectation value, or probability, of remaining in |110>\n",
    "probs_110 = [np.abs((~initial_state @ U_heis3(float(t)) @ initial_state).eval())**2 for t in ts]\n",
    "\n",
    "# Plot evolution of |110>\n",
    "plt.plot(ts, probs_110)\n",
    "plt.xlabel('time')\n",
    "plt.ylabel(r'probability of state $|110\\rangle$')\n",
    "plt.title(r'Evolution of state $|110\\rangle$ under $H_{Heis3}$')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04db5341-8ba0-433b-b671-d929aba1fbe0",
   "metadata": {},
   "source": [
    "### 1-5 Decomposition of $U_{\\text{Heis3}}(t)$ Into Quantum Gates\n",
    "\n",
    "To execute $U_{\\text{Heis3}}(t)$ on a [circuit model](https://qiskit.org/documentation/apidoc/circuit.html) quantum computer, we must decompose $U_{\\text{Heis3}}(t)$ into a product of single and two-qubit gates that are native to that quantum computer. There are a [variety of techniques](https://en.wikipedia.org/wiki/Hamiltonian_simulation) to do this. For this Open Science Prize, you will only use [Trotterization](https://en.wikipedia.org/wiki/Hamiltonian_simulation#Product_Formulas) (which is sometimes refered to as Suzuki-Trotterization decomposition, Trotter formulas, product formulas, and other names).\n",
    "\n",
    "Below, we show an example Trotterization as outlined in \\[1-2\\].\n",
    "\n",
    "Since the [Pauli operators do not commute](https://en.wikipedia.org/wiki/Pauli_matrices#Commutation_relations) with each other, the exponential $U_{\\text{Heis3}}(t)$ cannot be split into a product of simpler exponentials. However, we can approximate $U_{\\text{Heis3}}(t)$ as a product of simpler exponentials through Trotterization. Consider a subsystem of 2 spin-1/2 particles within the larger 3 spin system. The Hamiltonian on spins $i$ and $j$ ($i,j \\in \\{0,1,2\\}$) would be $H^{(i,j)}_{\\text{Heis2}} = \\sigma_x^{(i)}\\sigma_x^{(j)} + \\sigma_y^{(i)}\\sigma_y^{(j)} + \\sigma_z^{(i)}\\sigma_z^{(j)}$. Rewritting $U_{\\text{Heis3}}(t)$ in terms of the two possible subsystems within the total $N=3$ system you will simulate,\n",
    "\n",
    "$$\n",
    "U_{\\text{Heis3}}(t) = \\exp\\left[-i t \\left(H^{(0,1)}_{\\text{Heis2}} + H^{(1,2)}_{\\text{Heis2}} \\right)\\right].\n",
    "$$\n",
    "\n",
    "$H^{(0,1)}_{\\text{Heis2}}$ and $H^{(1,2)}_{\\text{Heis2}}$ do not commute, so $U_{\\text{Heis3}}(t) \\neq \\exp\\left(-i t H^{(0,1)}_{\\text{Heis2}}\\right) \\exp\\left(-i t H^{(1,2)}_{\\text{Heis2}} \\right)$. But, this product decomposition can be approximated with Trotterization which says $U_{\\text{Heis3}}(t)$ is approximately a short evolution of $H^{(0,1)}_{\\text{Heis2}}$ (time = $t/n$) and followed by a short evolution of $H^{(1,2)}_{\\text{Heis2}}$ (time = $t/n$) repeated $n$ times\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "U_{\\text{Heis3}}(t) &= \\exp\\left[-i t \\left(H^{(0,1)}_{\\text{Heis2}} + H^{(1,2)}_{\\text{Heis2}} \\right)\\right] \\\\\n",
    "U_{\\text{Heis3}}(t) &\\approx \\left[\\exp\\left(\\dfrac{-it}{n}H^{(0,1)}_{\\text{Heis2}}\\right) \\exp\\left(\\dfrac{-it}{n}H^{(1,2)}_{\\text{Heis2}} \\right)\\right]^n.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "$n$ is the number of Trotter steps, and as $n$ increases, the approximation becomes more accurate. (Note that how a untiary is split up into subsystems for Trotterization is not necessarily unique.) The decomposition goes further. Within each 2 spin subsystems, the Pauli operator pairs ($\\sigma_x^{(i)}\\sigma_x^{(j)}$, $\\sigma_y^{(i)}\\sigma_y^{(j)}$, and $\\sigma_z^{(i)}\\sigma_z^{(j)}$) commute. This means we can decompose the exponential of a subsystem Hamiltonian ($H^{(i,j)}_{\\text{Heis2}}$) into a product of even simpler exponentials getting us closer to a gate implementation of $U_{\\text{Heis3}}(t)$\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "U_{\\text{Heis3}}(t) &\\approx \\left[\\exp\\left(\\dfrac{-it}{n}H^{(0,1)}_{\\text{Heis2}}\\right) \\exp\\left(\\dfrac{-it}{n}H^{(1,2)}_{\\text{Heis2}} \\right)\\right]^n \\\\\n",
    "U_{\\text{Heis3}}(t) &\\approx \\left[\\exp\\left(\\dfrac{-it}{n}\\left(\\sigma_x^{(0)}\\sigma_x^{(1)} + \\sigma_y^{(0)}\\sigma_y^{(1)} + \\sigma_z^{(0)}\\sigma_z^{(1)}\\right)\\right) \\exp\\left(\\dfrac{-it}{n}\\left(\\sigma_x^{(1)}\\sigma_x^{(2)} + \\sigma_y^{(1)}\\sigma_y^{(2)} + \\sigma_z^{(1)}\\sigma_z^{(2)}\\right)\\right)\\right]^{n} \\\\\n",
    "U_{\\text{Heis3}}(t) &\\approx \\left[\\exp\\left(\\dfrac{-it}{n}\\sigma_x^{(0)}\\sigma_x^{(1)}\\right) \\exp\\left(\\dfrac{-it}{n}\\sigma_y^{(0)}\\sigma_y^{(1)}\\right) \\exp\\left(\\dfrac{-it}{n}\\sigma_z^{(0)}\\sigma_z^{(1)}\\right) \\exp\\left(\\dfrac{-it}{n}\\sigma_x^{(1)}\\sigma_x^{(2)}\\right) \\exp\\left(\\dfrac{-it}{n}\\sigma_y^{(1)}\\sigma_y^{(2)}\\right) \\exp\\left(\\dfrac{-it}{n}\\sigma_z^{(1)}\\sigma_z^{(2)}\\right)\\right]^{n}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "For simplicity, and to use more common notation, let's rename the products with $XX(2t) = \\exp\\left(-it \\sigma_x\\sigma_x\\right)$, $YY(2t) = \\exp\\left(-it \\sigma_y\\sigma_y\\right)$, and $ZZ(2t) = \\exp\\left(-it \\sigma_z\\sigma_z\\right)$ and rewrite the Trotterized $U_{\\text{Heis3}}(t)$\n",
    "$$\n",
    "U_{\\text{Heis3}}(t) \\approx \\left[XX\\left(\\frac{2t}{n}\\right)^{(0,1)} YY\\left(\\frac{2t}{n}\\right)^{(0,1)} ZZ\\left(\\frac{2t}{n}\\right)^{(0,1)} XX\\left(\\frac{2t}{n}\\right)^{(1,2)} YY\\left(\\frac{2t}{n}\\right)^{(1,2)} ZZ\\left(\\frac{2t}{n}\\right)^{(1,2)} \\right]^{n}\n",
    "$$\n",
    "And that's it! We have now approximately decomposed $U_{\\text{Heis3}}(t)$ into two-qubit gates $XX(t)$, $YY(t)$, and $ZZ(t)$. These gates are not native to superconducting qubits, but in Section 2, they will be further decomposed into native single and two-qubit gates with *further details found in the supplementary material.*\n",
    "\n",
    "\n",
    "\\[1\\] Y. Salathe, et al., *Digital Quantum Simulation of Spin Models with Circuit Quantum Electrodynamics*, [Phys. Rev. X **5**, 021027 (2015)](https://link.aps.org/doi/10.1103/PhysRevX.5.021027)\n",
    "\n",
    "\\[2\\] F. Tacchino, et al., *Quantum Computers as Universal Quantum Simulators: State-of-the-Art and Perspectives*, [Adv. Quantum Technol. *3* 3 (2020)](https://doi.org/10.1002/qute.201900052) \\[[free arxiv version](https://arxiv.org/abs/1907.03505)\\]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "509e3761-6196-4b8f-aa5d-434bf1070eeb",
   "metadata": {},
   "source": [
    "# 2 The Open Science Prize\n",
    "### 2-1 Contest Details\n",
    "The Open Science Prize requirements and judging criteria can be found in the terms and conditions when registering on the [website](https://ibmquantumawards.bemyapp.com/#/event). Below is an outline of those requirements.\n",
    "\n",
    "You are to time evolve the state $|110\\rangle$ to time $t=\\pi$ under the $XXX$ Heisenberg model Hamiltonian $H_{\\text{heis3}}$ as defined in Section 1. This will be done on the ibmq_jakarta quantum computer. From left to right, the qubits in $|110\\rangle$ correspond to ibmq_jakarta's qubits 5, 3, and 1. The evolution must be done using Trotterization with at least 4 steps ($n\\geq4$). You will be judged on the clarity, originiality, and creativity of your solution, as well as the fidelity of the evolved state at time $t=\\pi$ using [state tomography](https://qiskit.org/documentation/tutorials/noise/8_tomography.html). State tomography is a way of determining the exact quantum state--even phase information--at the end of a quantum circuit by running the circuit multiple times and measuring in different bases. The closer your final quantum state to the expected state, the higher the fidelity. Ibmq_jakarta's qubits 0, 2, 4, and 6 will not be evaluated by the state tomography calculation, so you may use them as you wish. You do not need to follow the same decomposition given in the example above, but your solution must be a Trotterization of $U_{\\text{heis3}}(t)$ with at least 4 steps ($n\\geq4$). You may find it useful to implement techniques that reduce errors and improve qubit coherence."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c794b8f1-6de7-497b-a000-0b37f9d094b3",
   "metadata": {},
   "source": [
    "### 2-2 Import Qiskit\n",
    "\n",
    "Feel free to import packages as needed. However, only free packages and those obtained through ```pip install``` or ```conda install``` are allowed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8b3e598e-1156-4263-abc6-890e554cf4fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing standard Qiskit modules\n",
    "from qiskit import QuantumCircuit, QuantumRegister, IBMQ, execute, transpile\n",
    "from qiskit.providers.aer import QasmSimulator\n",
    "from qiskit.tools.monitor import job_monitor\n",
    "from qiskit.circuit import Parameter\n",
    "\n",
    "# Import state tomography modules\n",
    "from qiskit.ignis.verification.tomography import state_tomography_circuits, StateTomographyFitter\n",
    "from qiskit.quantum_info import state_fidelity\n",
    "\n",
    "# suppress warnings\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "900d2ab8-b568-4d28-87d8-473962fdc998",
   "metadata": {},
   "source": [
    "### 2-3 Quantum Devices\n",
    "Connect to IBM provider and connect to a real or simulated backend. Final submissions must be run on a real backend, but simulated devices are faster for debugging and testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7345d1e9-4803-4bd8-a578-0af5294cb22a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load IBMQ Account data\n",
    "\n",
    "# IBMQ.save_account(TOKEN)  # replace TOKEN with your API token string (https://quantum-computing.ibm.com/lab/docs/iql/manage/account/ibmq)\n",
    "provider = IBMQ.load_account()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "78e0af5d-6da9-4a08-912a-5a172535e03b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get backend for experiment\n",
    "provider = IBMQ.get_provider(hub='ibm-q-community', group='ibmquantumawards', project='open-science-22')\n",
    "jakarta = provider.get_backend('ibmq_jakarta')\n",
    "# properties = jakarta.properties()\n",
    "\n",
    "# Simulated backend based on ibmq_jakarta's device noise profile\n",
    "sim_noisy_jakarta = QasmSimulator.from_backend(provider.get_backend('ibmq_jakarta'))\n",
    "\n",
    "# Noiseless simulated backend\n",
    "sim = QasmSimulator()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4e8a24c-5846-4de9-95ea-3033653c4037",
   "metadata": {},
   "source": [
    "### 2-4 Decomposition of $U_{\\text{Heis3}}(t)$ into Quantum Gates (Example)\n",
    "\n",
    "The following circuit code is written based on the example given in Section 1. This is where you write your solution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "207beaa7-58e9-40c1-a35c-f359ddb7aa8b",
   "metadata": {},
   "source": [
    "YOUR TROTTERIZATION GOES HERE -- START (beginning of example)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9156980e-d3a7-4494-8ad4-7ec15ca92e7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameterize variable t to be evaluated at t=pi later\n",
    "t = Parameter('t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ec328a9f-809c-49d6-8400-fa09579c1d8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build a subcircuit for XX(t) two-qubit gate\n",
    "XX_qr = QuantumRegister(2)\n",
    "XX_qc = QuantumCircuit(XX_qr, name='XX')\n",
    "\n",
    "XX_qc.ry(np.pi/2,[0,1])\n",
    "XX_qc.cnot(0,1)\n",
    "XX_qc.rz(2 * t, 1)\n",
    "XX_qc.cnot(0,1)\n",
    "XX_qc.ry(-np.pi/2,[0,1])\n",
    "\n",
    "# Convert custom quantum circuit into a gate\n",
    "XX = XX_qc.to_instruction()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "52977709-1b68-4ebf-98de-dc5c88265acb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build a subcircuit for YY(t) two-qubit gate\n",
    "YY_qr = QuantumRegister(2)\n",
    "YY_qc = QuantumCircuit(YY_qr, name='YY')\n",
    "\n",
    "YY_qc.rx(np.pi/2,[0,1])\n",
    "YY_qc.cnot(0,1)\n",
    "YY_qc.rz(2 * t, 1)\n",
    "YY_qc.cnot(0,1)\n",
    "YY_qc.rx(-np.pi/2,[0,1])\n",
    "\n",
    "# Convert custom quantum circuit into a gate\n",
    "YY = YY_qc.to_instruction()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4fd739df-0dc4-48ef-823e-230fb0c8fd64",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build a subcircuit for ZZ(t) two-qubit gate\n",
    "ZZ_qr = QuantumRegister(2)\n",
    "ZZ_qc = QuantumCircuit(ZZ_qr, name='ZZ')\n",
    "\n",
    "ZZ_qc.cnot(0,1)\n",
    "ZZ_qc.rz(2 * t, 1)\n",
    "ZZ_qc.cnot(0,1)\n",
    "\n",
    "# Convert custom quantum circuit into a gate\n",
    "ZZ = ZZ_qc.to_instruction()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4335d6d3-4a8e-4140-aad0-97b9b189905a",
   "metadata": {},
   "source": [
    "Combine subcircuits into single gate representing one ($n=1$) trotter step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "8ef03f79-2c47-49da-93d5-0c0ae3ea5872",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine subcircuits into a single multiqubit gate representing a single trotter step\n",
    "num_qubits = 3\n",
    "\n",
    "Trot_qr = QuantumRegister(num_qubits)\n",
    "Trot_qc = QuantumCircuit(Trot_qr, name='Trot')\n",
    "\n",
    "for i in range(0, num_qubits - 1):\n",
    "    Trot_qc.append(ZZ, [Trot_qr[i], Trot_qr[i+1]])\n",
    "    Trot_qc.append(YY, [Trot_qr[i], Trot_qr[i+1]])\n",
    "    Trot_qc.append(XX, [Trot_qr[i], Trot_qr[i+1]])\n",
    "\n",
    "# Convert custom quantum circuit into a gate\n",
    "Trot_gate = Trot_qc.to_instruction()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0641237c-048e-41e7-b47c-88a364683fb7",
   "metadata": {},
   "source": [
    "YOUR TROTTERIZATION GOES HERE -- FINISH (end of example)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9fb21ba-a2ea-4595-ba76-a4cef87ea666",
   "metadata": {},
   "source": [
    "### 2-5 Trotterized Time Evolution\n",
    "\n",
    "Time evolve the state $|110\\rangle$ from $\\theta=0$ to $\\theta=\\pi$ under $H_{\\text{heis3}}$. Although jakarta has 7 qubits, only the qubits 1, 3, and 5 will be used for the state $|110\\rangle$. From left to right, the qubits in $|110\\rangle$ corresponds to the qubits 5, 3, and 1 on jakarta. State tomography will evaluate how well the quantum simulation matches to the expected state at $\\theta= \\pi$. (Note qubits 0, 2, 4, and 6 are not part of the state tomography calculation.)\n",
    "\n",
    "You choose how many trotter steps (```trotter_steps```) to take in the simulation, but it must be at least 4. The code below is written to be generic enough for you may add to it if necessary for your solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c6135540-dcea-4313-8ce8-0eb8e3ee23e7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">                                                                    ░          \n",
       "q4_0: ──────────────────────────────────────────────────────────────░──────────\n",
       "           ┌────────────┐┌────────────┐┌────────────┐┌────────────┐ ░ ┌─┐      \n",
       "q4_1: ─────┤0           ├┤0           ├┤0           ├┤0           ├─░─┤M├──────\n",
       "           │            ││            ││            ││            │ ░ └╥┘      \n",
       "q4_2: ─────┤            ├┤            ├┤            ├┤            ├─░──╫───────\n",
       "      ┌───┐│            ││            ││            ││            │ ░  ║ ┌─┐   \n",
       "q4_3: ┤ X ├┤1 Trot(π/4) ├┤1 Trot(π/4) ├┤1 Trot(π/4) ├┤1 Trot(π/4) ├─░──╫─┤M├───\n",
       "      └───┘│            ││            ││            ││            │ ░  ║ └╥┘   \n",
       "q4_4: ─────┤            ├┤            ├┤            ├┤            ├─░──╫──╫────\n",
       "      ┌───┐│            ││            ││            ││            │ ░  ║  ║ ┌─┐\n",
       "q4_5: ┤ X ├┤2           ├┤2           ├┤2           ├┤2           ├─░──╫──╫─┤M├\n",
       "      └───┘└────────────┘└────────────┘└────────────┘└────────────┘ ░  ║  ║ └╥┘\n",
       "q4_6: ──────────────────────────────────────────────────────────────░──╫──╫──╫─\n",
       "                                                                    ░  ║  ║  ║ \n",
       "c0: 3/═════════════════════════════════════════════════════════════════╩══╩══╩═\n",
       "                                                                       0  1  2 </pre>"
      ],
      "text/plain": [
       "                                                                    ░          \n",
       "q4_0: ──────────────────────────────────────────────────────────────░──────────\n",
       "           ┌────────────┐┌────────────┐┌────────────┐┌────────────┐ ░ ┌─┐      \n",
       "q4_1: ─────┤0           ├┤0           ├┤0           ├┤0           ├─░─┤M├──────\n",
       "           │            ││            ││            ││            │ ░ └╥┘      \n",
       "q4_2: ─────┤            ├┤            ├┤            ├┤            ├─░──╫───────\n",
       "      ┌───┐│            ││            ││            ││            │ ░  ║ ┌─┐   \n",
       "q4_3: ┤ X ├┤1 Trot(π/4) ├┤1 Trot(π/4) ├┤1 Trot(π/4) ├┤1 Trot(π/4) ├─░──╫─┤M├───\n",
       "      └───┘│            ││            ││            ││            │ ░  ║ └╥┘   \n",
       "q4_4: ─────┤            ├┤            ├┤            ├┤            ├─░──╫──╫────\n",
       "      ┌───┐│            ││            ││            ││            │ ░  ║  ║ ┌─┐\n",
       "q4_5: ┤ X ├┤2           ├┤2           ├┤2           ├┤2           ├─░──╫──╫─┤M├\n",
       "      └───┘└────────────┘└────────────┘└────────────┘└────────────┘ ░  ║  ║ └╥┘\n",
       "q4_6: ──────────────────────────────────────────────────────────────░──╫──╫──╫─\n",
       "                                                                    ░  ║  ║  ║ \n",
       "c0: 3/═════════════════════════════════════════════════════════════════╩══╩══╩═\n",
       "                                                                       0  1  2 "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The final time of the state evolution\n",
    "target_time = np.pi\n",
    "\n",
    "# Number of trotter steps\n",
    "trotter_steps = 4  ### CAN BE >= 4\n",
    "\n",
    "# Initialize quantum circuit for 3 qubits\n",
    "qr = QuantumRegister(7)\n",
    "qc = QuantumCircuit(qr)\n",
    "\n",
    "# Prepare initial state (remember we are only evolving 3 of the 7 qubits on jakarta qubits (q_5, q_3, q_1) corresponding to the state |110>)\n",
    "qc.x([3,5])  # DO NOT MODIFY (|q_5,q_3,q_1> = |110>)\n",
    "\n",
    "# Simulate time evolution under H_heis3 Hamiltonian\n",
    "for _ in range(trotter_steps):\n",
    "    qc.append(Trot_gate, [qr[1], qr[3], qr[5]])\n",
    "\n",
    "# Evaluate simulation at target_time (t=pi) meaning each trotter step evolves pi/trotter_steps in time\n",
    "qc = qc.bind_parameters({t: target_time/trotter_steps})\n",
    "\n",
    "# Generate state tomography circuits to evaluate fidelity of simulation\n",
    "st_qcs = state_tomography_circuits(qc, [qr[1], qr[3], qr[5]])\n",
    "\n",
    "# Display circuit for confirmation\n",
    "# st_qcs[-1].decompose().draw()  # view decomposition of trotter gates\n",
    "st_qcs[-1].draw()  # only view trotter gates"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1742151c-1261-4b42-80a4-da2c0ad90ac5",
   "metadata": {},
   "source": [
    "### 2-6 Execute\n",
    "For your final submission, you will need to execute your solution on a real backend with 8 repetitions. For faster debugging, considering using a simulated backend and/or 1 repetition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "79c3ef21-c143-418a-8b50-f0ab39c5ae35",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Job ID 69578ef6-8d0d-4413-88c8-46179605d4d6\n",
      "Job ID 0d5068ed-6957-4a5e-9956-e5c924377a26\n",
      "Job ID 568ff236-6b75-4ad5-b6a3-c67a733b1f30\n",
      "Job ID 8f795bf9-9092-423c-bce0-8bbafd0b588e\n",
      "Job ID bd3b6064-dca9-4611-9817-1050c1ff6c41\n",
      "Job ID eb0ffe57-fabd-46e6-a7bd-468d5a5645e7\n",
      "Job ID 2a44d28d-2fee-4b40-b576-be625eee43fc\n",
      "Job ID df446ee1-4eb5-48ec-a246-903d8ca63c6d\n"
     ]
    }
   ],
   "source": [
    "shots = 8192\n",
    "reps = 8\n",
    "backend = sim_noisy_jakarta\n",
    "# reps = 8\n",
    "# backend = jakarta\n",
    "\n",
    "jobs = []\n",
    "for _ in range(reps):\n",
    "    # execute\n",
    "    job = execute(st_qcs, backend, shots=shots)\n",
    "    print('Job ID', job.job_id())\n",
    "    jobs.append(job)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be7acc76-86f0-4b1f-94e4-e4e1014f3956",
   "metadata": {},
   "source": [
    "We can monitor the status of the jobs using Qiskit's job monitoring tools."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "1b142880-af5b-4913-8f50-9248e5d1646f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n"
     ]
    }
   ],
   "source": [
    "for job in jobs:\n",
    "    job_monitor(job)\n",
    "    try:\n",
    "        if job.error_message() is not None:\n",
    "            print(job.error_message())\n",
    "    except:\n",
    "        pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4186e06-db3d-4d0d-aef4-51906f41df61",
   "metadata": {},
   "source": [
    "### 2-7 Results Analysis\n",
    "Extract the results for the completed jobs and compute the state tomography fidelity for each repetition. You may choose to include other post-processing analyses here as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f87c0c22-1f37-4f56-9ec5-7e3c5d7a8384",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "state tomography fidelity = 0.0378 ± 0.0013\n"
     ]
    }
   ],
   "source": [
    "# Compute the state tomography based on the st_qcs quantum circuits and the results from those ciricuits\n",
    "def state_tomo(result, st_qcs):\n",
    "    # The expected final state; necessary to determine state tomography fidelity\n",
    "    target_state = (One^One^Zero).to_matrix()  # DO NOT MODIFY (|q_5,q_3,q_1> = |110>)\n",
    "    # Fit state tomography results\n",
    "    tomo_fitter = StateTomographyFitter(result, st_qcs)\n",
    "    rho_fit = tomo_fitter.fit(method='lstsq')\n",
    "    # Compute fidelity\n",
    "    fid = state_fidelity(rho_fit, target_state)\n",
    "    return fid\n",
    "\n",
    "# Compute tomography fidelities for each repetition\n",
    "fids = []\n",
    "for job in jobs:\n",
    "    fid = state_tomo(job.result(), st_qcs)\n",
    "    fids.append(fid)\n",
    "    \n",
    "print('state tomography fidelity = {:.4f} \\u00B1 {:.4f}'.format(np.mean(fids), np.std(fids)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a7a0994c-f8bf-4613-8b1e-c78eb2864b8e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td><code>qiskit-terra</code></td><td>0.18.3</td></tr><tr><td><code>qiskit-aer</code></td><td>0.9.1</td></tr><tr><td><code>qiskit-ignis</code></td><td>0.6.0</td></tr><tr><td><code>qiskit-ibmq-provider</code></td><td>0.18.1</td></tr><tr><td><code>qiskit-aqua</code></td><td>0.9.5</td></tr><tr><td><code>qiskit</code></td><td>0.32.1</td></tr><tr><td><code>qiskit-nature</code></td><td>0.1.3</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.8.8 (default, Apr 13 2021, 12:59:45) \n",
       "[Clang 10.0.0 ]</td></tr><tr><td>OS</td><td>Darwin</td></tr><tr><td>CPUs</td><td>8</td></tr><tr><td>Memory (Gb)</td><td>64.0</td></tr><tr><td colspan='2'>Mon Nov 29 09:12:47 2021 EST</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "%qiskit_version_table"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "toc-autonumbering": false,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}