ibmq-qsim-challenge_with_my_solution.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://en.wikipedia.org/wiki/Spin-1/2) 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\n",
    "# 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": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEnCAYAAAByjp6xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABHQ0lEQVR4nO3dd3xV9f348dc7mxAChEDYhBmGLBkyJSCKq+C2y7qq1fZbrbVDq7Wun6O2dtmhtdY666xbcRGFiCggyIawNwkjEMjO+/fHOdFrvEnuvbm5597k/Xw87uMmZ75vzs193/OZoqoYY4wxkRDndQDGGGNaD0s6xhhjIsaSjjHGmIixpGOMMSZiLOkYY4yJGEs6xhhjIsaSjjHGmIixpGOMMSZiLOl4SESyRURF5NFYPkdzEpFkEblLRDaKSKX7WkZ5HZdXRORW92+Q7XUsLY2IXOL+bS/xOpaWrNUlHZ8P4YYeeV7HGSgRyXVjvtXrWJrJz4Ebga3AvcBtwJ6mHjScf7douAYi8l0ReUhElohIuRvPec2wzzgReUNEDorIURH5RES+7We7M93jndnU19aaiMg6EakRkfb1rG/vrt8Q6djCJcHrADy0AXiqnnVbIhhHc9sJDAGKvQ4kRKcBJcAsVa30OpgodifQBygE9gK9wr2PiOQCc4EK4L8476lzgCdFJFtV7/LZ/F3gKHAW8FrgL6P1EpE0YCCwXlXr+38dCwiwNGKBhVlrTjrrVfVWr4Nobu4H9Vqv42iCbsB+SziN+j7Oe3qbe8f1m3DuIyIJwMOAAieq6mfu8luBj4HbROQ5Vd0AoKplIjIX+IaIxKlqTegvrdUYjZNQPm1gm3Huc8wmnVZXvBYMETnRLSL4az3rR7vrH/NZligivxCRlSJS6hZDvCUi0wI8Z73lynWLcdznee7q39QpIsx2t/FbpxNMnL7nFZGxIvKOiBwRkWIR+V8w9QuBnre27gLoC/QJpuhTRM4XkfkiUigiZSKy0y0SOjmIv1uSiFzjvtadIlIhIrtF5CkRGegbZ2PH8tn2XBHJc/9upSLymYhcGejfriGq+q6qbmvGfWYA/YGnahOOe4wS4HacL7CX1tnnJaALMCnQmIJ5/9ddFuh7U0TaisjvRWSXex2Wisi5jcQV0LWrE88UEXnX3edgAC//ePe5RSed1nynE4j5wDbgAhG5VlWr6qz/jvv8BICIxAH/A84AVgN/ATKAC4H3ReQ7qvrfMMaXB2QDFwMfuL/XOlTfTk2IcxzwC5wP2QdxvpmdBQwXkeNUtayhYIM8b+1r+Yn7/Ef3eUsj5/gR8ACwEXgGOAJ0ByYDs4B3COzvlgHcD3wIvIpTlJQDXACcKiJjVHVzgMdCRO4DfoZTN/UMcAw4GXhQRIao6nUNva4okOs+v+1n3Vvuc90vLK8BVcDZwILmCesLAb033ffgazivZynwGM7d9FM4742vCfHaTQZ+BbznxtM5gNfQKpIOqtqqHjgfEAqsB26t5zHBZ/u73e3PrHOcOJz6kt1AvLvsEnfbuUCCz7bDcN6oh4B0P7E86rOs9hiX+Ik91113a0PL6nm9/s4RaJy151DgwjrHf8xd/s0A/vZBndddtwXYEsT1Xepel1Q/6zoF8XdLBrr7WT4N54P04SCONctd/xKQ4rM80V2mwLgAX9+t7vbZAWxzXhB/twb3AZ5z14+pZ/0eYJ+f5e8BBUHEEer7P6D3JnCZz7WI81k+w+c4l/gsD+ra1Ynne4G+bnffFe5+z+F8ka37eMpdvzWY40bbozUXrw3EKcP295jgs92T7vN3+KpcnG/Q/1XVanfZxe7zDepzV6Sqq4BHgPY43768FmqcH6rqM3WWPeI+j6u7cRjPG6wKoLruQlXdH+gBVLVcVXf5Wf4Bzl3azCDi+RHOh8UP1OduUJ16qpvdXy8M4nheqG1NVV8Fd7HPNr5eAvqLyPDmCMpHoO/N77rPN6tPPZOqvo//O51Qr90SVX3Mz3K/RCQFGOz+eh7O503dx7fc9bF7l0PrLl57XVUbbc6pqitFZDkwW0TS1CnDhjpFa66RwBH1KfP2kYfzBh7ZhJjDJdQ4l/hZtsN97tCM5w3GM8A9wEoR+a973IWqeizYA4nIGJwim8k4dROJPqsrgjjUCTjFfFeLSN11tcccXHdFC/ES8GecIrYVzXieQN+bI4FiVV3pZ/sFOMVmvkK9dosbCtaPkTifx0+q6nf9bSAiN+CUvDRb0hGRJ4FNqvrr5jpHa046wXgCuA+neehjIpIMnAusVVXfN3s6sLmeY+zx2cZrocZ52M+y2juW+GY8bzB+CxwArsL5JnozUC4izwPXq+reQA4iIlNwioZqcIoDC3CaACtOEVCfIGLKwPlf+00D27QN4nheqL3D8Xc3U7v8a3dBqrpdRHYCE5srMFeg78321P8e9PfeCPXaBfQ+81Fbn7OsgW1Guc9fSzoisgJ4VVV/VWf5k0BbVT0rkCBUtW6Jjl8i8nPgCpwvYwqsBG5R1XmN7duai9eC8TTOh0/tBTkT5837ZJ3tDuNcBH+yfLZpSO0tv78vBPX9wwcrHHFG5XnV8U9VHeOe63ycyu/vAM8GcagbgSTgJFU9S1V/pqq/UaeZfWmQYR0G9qqqNPCYHuQxI622M+LAuitEJB3n+n2tw6KIDAB64DSyCERzv/+Lqb9SP8vPslCvnQYZV8hJR0TaAUPrLneNBhYFGUsgXgImqmoHIBP4K/CGiHRtbEdLOgFQ1Z04xTQnuX/U7+C8qeomnWVAuvgfpmWazzYNOeg+9/CzbrSfZbV1F4HcadRaRtPjDEVEz6uqhar6vKrOdo97onzZ07uxv1t/nP5BH/kuFJEsd52vxo71CZAVTNPyKFSbNE7xs+7UOtv4Ost9finA8wT7/g/WcqC9iBznZ90UP8side0aTDoikoqT8Peo6u46q8fifJbXTUapOC0uP/FZNkFE3haRve7jDyIS764bKc5QU23c39NE5C8isl1EStwuA7cCqOoGnzpSwfkfSMFpuNQgSzqBewLnQ+Vq4HTgI3WazPqqrTi8u/ZCAojIEJyOeMXAy42cZwluixu3crH2GAOBa/1sf8B9DqQHejjjDEWzn1f89zNqg/MtuYovE0Rjf7dtQIYbW+1xknCaYyfW2baxY/3Fff6XiHTwE192DCSk94BNwLd9vzSI04v+Fpy/7b/97HcWTjF0oB2Ug33/B6u2DvZOt/l07fFn8PX6HIjAtRORROA4YKeqFtWz2Qj8JBbXeOCQqm7ysw+4TbDFGZJoLk4n357ufqcBP3S3GwesUtXaO/k/4dzFjFTVNJzk9sX/p9sP6RBO/eazwIv4JLj6tOY6nUFS/1hZZap6T51lLwB/wyl2SeTrdzngfKieh9MPZZmIvAF0BL6J0wT3MlVtsPhIVXeJyNPAt4ElIvIWTjHR2Tj9Iep2YlsL7ML5Jy3HqTxV4C9a/1AaTY4zRJE478vuP8IinH4VKTj/WH2BB3wagjT2d3sA50MoX0SewflQnYlz7Zfz1QYPDR5LVd8Qkbtx3jsbxOmpvwOnmGcITmvJb9OE4ZdE5Pt8+U19lPv8I/ly7LOXVPWlUPdR1Sp3+7nAfPc9ehinnrMvTmuwrxSviUgXnLqc3wb6OkJ4/wfrUeAiYA7wqYi8g9NP50LgdZz3pm88zX7tcBJOEiHW5+A0dmgnInUTVhJOwj/sfmH6J07L0dpi5q3ue3sKTnIdx1cbQAwGtuP8D9WW+OysXamqC4AObvHeN4EkDWTkiUi2z46GB1/2W2nocaiefZ9111fg0+ejzjaJwC+BVUAZTt+TucC0BmJ5tM7yNjjfMva4x1iO88bOxU9/EPdNl4fzIVD7GrIbOUcwcfo9b0PHb+DvH/B53e23EFw/nauBV9z9ynDGFVuA80Ejgf7d3PUXAJ/h9CHag/OBleXuo8Ecy93mNJwPtiL3PbQTp0jqeiAzwNd3az3HfpSG39P+rl0o+4wH3nSv2zGcb9HfqSfWy93jnBDk/2jA7/9Q3ptAGk7H39049XOf4SSzS6i/j1BA166heBp4vbV/pzsb2OYf7jZn+1m3E6eFYM86j5eBf7vbTHf3L3avXe2jBHjQ3WYpTtPw2uMOwflyXYhzB3pVI69jjb/46j7E3dgYEwPkyzHS+qrqFm+jaZiIvIpTV9FT7YOmWYhID5w7rzNV9fU667YCd6vqP0TkAuC3qppdz3FScL4wTdSvtsitHcXhOzglFfW+78QZ+fpf+vVSoq+wOh1jTNiJSFuc4siXLeE0qxNwWvzl+y4UkV5Ab75subYE6CoiF4sz/mG8iPQRkVnu+lHucT539z9HnLElE3FKJ3rhFCHvcNf/WER6iKODiNSOWO5vmKSvsKRjjGkOJ+L0a3rB60BauBNwKv8P1Vk+FafocAWAqm7EKS6+Bqfhy36cYuie7vbjgM/1y9Hcx+E0DCjGaVQzFpihX44kMgWnaPUozpBiE4DTVbXRjqtWvGZMDIml4jVj/GnNrdeMiUV57vMhD2MwJmR2pwNkZmZqdnZ2SPsePXqUtm2jfQSTls2ugffsGnjPi2uwZMmSIlUNZNqGL9idDpCdnc3ixcGOz+fIy8sjNzc3vAGZoNg18J5dA+95cQ3cFnJBsYYExhhjIsaSjjHGmIixpGOMMSZioibpiEhPd0TThSJyTEQ00IH0RCRORG4UkS0iUiYiy0WkqWM0GWOMCbOoSTrAAJzOSweB+UHuewfOmFQP4IyR9DHwnIicHs4AjTHGNE00tV77UFWz4IvRb/3N2/E17ki2PwPuUdXfuYvniTN51D3AG80RrDHGmOBFzZ2OBjIktn+zcIbwfqLO8ieA4SLSt0mBGWOMCZuoSTpNMAwox5nD3tcq93loc5343/mbWbCzksIj5c11CmOMCbuS8irmrtrDb98KdG698Imm4rVQZeDMf1N3aIUDPuu/RkSuBK4EyMrKIi8vL6iTqioP55eys0R5eMW79G0fx4xeCUztWXdSSdPcSkpKgr5+JrzsGngvkGuwZn81r22qYO2BGqoV2iTAYNlNerJEJkhaRtIJiao+BDwEMHbsWA2lJ+/8acrjr73P4ba9eXPlHv618jDSoQe/On0IcXGRu4itnfWG955dA+81dg2e/XQ7v397BVnpKVw+tRvTc7owNrsjifGRLfBqCUnnIM6UqVLnbqf2DueAn33CIi5O6JMeT27uQH44fQB3vLaahxdsZsfBUv5w4SjaJMU316mNMSYgqsrv317PA/MKmDowk79953japXhXIhNSihORdiLyTRF5WkSeEJHzRCQt3MEFaBWQDPSvs7y2Lmd1JIKIjxNunT2MW84cytzVe7j00U+orrHBVI0x3rrjtTU8MK+Ab47rxSOXjPM04UAQScftvPlDEZmLM2f2g+7+KcC/gUIReVNEfiAi3ZonXL/eAipxplP19V1gpapujmAsXDalL/eeO4KPNx3gHx9sjOSpjTHmK+at3ccj+Zu5eGIf7j5neMSL0vxptHhNRDrjfLCPAvbgzDb3B+B9Va1wt0nGmZp2Ns4EU38VkcXALFUtDjQYETnP/XGM+3yaiBQChar6gbtNFfAfVb0cQFX3icj9wI0icgRYClwIzHDjibjzx/Rk/oYi/vDOeqYOzGREzw5ehGGMacX2l5Tz8+c/Z3DXdtx4+hBEoqOeOZA6nXjgTeAqVf3U3waqWg68DrwuIlfhTF06m+CL756r8/vf3OcPgFyfeOpWltwElADXAl2BdcAFqvpakOcPCxHhzjnHsXjLAX7y32W8ds0UUpNaQvWZMSYWqCq/fGEFh0sreeL740lJjJ765UaTgqruUdWb60s4frZXVV2oqjeq6sFgglFVqeeRW2ebS+rsV62qd6pqH1VNVtURqvp8MOcOt/apifz+gpFs3n+Uu95Y42UoxphW5plPt/Pumr384tQcBndN9zqcr/C+gK8Fm9Q/k0sn9eXJRdtYv/eI1+EYY1qBo+VV/HbuOsb3zeCyydE3IEswDQnaiMhPRGSeiOwVkQr3sddd9hMRSW3OYGPRj2cMoG1SAn94Z73XoRhjWoFHP9rCgaMV3Hja4KjsLxhQ0hGRXsDnwH2AAM8D97qP2mKs3wLLRaR3M8QZszq2TeLyKX15c+UeVu4MuE2FMcYErbi0kgc/2MjMIV0Y3buj1+H4FWjt9h+BUmCgqm7xt4E7981LOC3bbC4bH5dP7cujH23h/nfW88gl47wOxxjTQv1r/iYOl1Vx3cmDvA6lXoEWr80Ebqov4QC4625xtzU+0lMS+cG0fry/dh9LtgbVtsIYYwJypEL514LNnDG8G8O6t/c6nHoFmnSC6Vpv3fD9uGRSNplpSdz/zjqvQzHGtEBvbq6ktLKa604e6HUoDQo06bwL/L+G5qZxi9fuAN4JQ1wtTmpSAlee2I/8gv2s3nXY63CMMS3IsYoq5m2v5IwR3RnQpZ3X4TQo0KTzE6ANsF5E5ovI30XkLvfxdxH5EFjvbnNdM8Ua8y4c25uUxDge/3ir16EYY1qQlz7bRWmVU6IS7QJKOqq6AxiBMy10OXAWcL37OBtn7LOfA6PcbY0f7VMTmT2yOy99tpPDZZVeh2OMaQFUlccWbqF3uziO793B63AaFXA/HVUtVdU/qepMVe3m9vxPVtWuqnqSu+5YcwbbElw0IZvSympeXGK52RjTdEu3HWTtniPM6J0QNeOrNcRGJIiw4T3bM7JXBx7/eCtfn+zUGGOC8/jCrbRLTmBit9gY3zGsSUdEThSR98N5zJboogl92Fh4lIWb9nsdijEmhhWVlPPGij2cO6YnyQnRf5cD4b/T6QxMC/MxW5wzR3SjQ2oiT1iDAmNMEzy7eDsV1TV8d0Ifr0MJWED3Y0EMbdO5CbG0GimJ8Zw/pieP5G+h8Eg5ndslex2SMSbGqCpPf7KNif06MaBLGjsiMkdy0wV6p7MF2BzA42/17G/qOG9ML6prlDdW7PY6FGNMDFq67RDbD5Ry3pieXocSlEBrnkqBD/lycM/6jAWubFJErURO13YM7tqOl5ft5OIYaFtvjIkuLy/bSXJCHLOO6+p1KEEJNOksB6pV9V8NbSQih7CkE7A5o3pw71tr2bb/GL072awQxpjAVFbX8Prnu5k5NIu05NhotVYr0OK1JcCYALeNjSYUUeAbI7sB8MrynR5HYoyJJQsKith/tII5I7t7HUrQAk069wDfbGwjVX1BVa3vT4B6dkxlfHYGLy3bZX12jDEBe2XZLtq3SSQ3p4vXoQQt0GFwdqrqB80dTGs0e1R3CvaVsHq3DQJqjGncsYoq5q7aw+nDu5KUEHvf8WMv4hbmjOHdSIgTXlm2y+tQjDEx4N01+zhWUc2cUT28DiUkNiKBxzq2TWLaoM68snwXNTVWxGaMadgry3bSrX0K47MzvA4lJDYiQRSYPao7u4vLWLrNZhU1xtTvcFklH6wv5MwR3YiLi802WzYiQRSYPrgLifHC26v3MjZGv70YY5rfvLX7qKxWTo2xvjm+Am3gvYXApqGWALczPtJTEpnUP5O5q/Zw42mDY2J4cmNM5L29ei+ZacmM7tXR61BCZiMSRIlThmVx0/9WsmFfCYOyonu6WWNM5JVXVZO3dh+zR/WI2aI1sBEJosbJQ7O4+aWVzF25x5KOMeZrPirYz9GKak4ZluV1KE1iIxJEiS7tUhjdqwNvr97rdSjGmCj09uo9pCUnMKl/J69DaRIbkSCKzBrWlRU7i9l5qNTrUIwxUaS6Rnln9V5yczqTnBDvdThNYiMSRJFThjktUt5ZtcfjSIwx0eSzbQcpKqn44jMiltldSRTpm9mWQVlpzF1lRWzGmC/NXbWHxHhhek7s90qxpBNlThnalU+2HODg0QqvQzHGRAFV5e3Ve5nUP5N2KYleh9NkgXYO/V5TTqKqjwVwjl7AH4CTcRojvAv8RFW3BbBvb+AOYDpOB9XtwLPA3ap6tAmhR9zJQ7N4YF4BH24ojNmxlYwx4bOx8Chb9x/j+1P7eR1KWATaZPrRJpxDgQaTjoikAu8D5cDF7j53AvNEZERDiUNE2uIkqETg18A2YBxwGzAQuLAJsUfc8B7t6dQ2iXlr91nSMcaQt24fQIsoWoMAk04EWqRdAfQDclS1AEBEPgc2AD8A7m9g38k4yWWWqr7tLpsnIhnAz0QkVVWPNV/o4RUXJ0zL6cy8tfuorlHiY7gTmDGm6eat28egrDR6dmwZswtHS53ObODj2oQDoKqbgXxgTiP7JrnPdSekOYTz+mLuU3t6ThcOHqtk+Y5DXodijPFQSXkVn2w+wPQYnKytPtGSdIYBK/0sXwUMbWTfd3HuiO4VkaEikiYiM4BrgX/EWp0OwIkDOxMnkLd2n9ehGGM8lF9QRGW1xuQMofUJtE6nuWUA/sb1PwA0OLKdqpaJyBTgBZwkVeth4P/q209ErsQdsicrK4u8vLwgQ3aUlJSEvG9DBnSI45XFmzg+aXfYj93SNNc1MIGza9A8nlpZTko8HN36OXnbGy60iZVrEC1JJ2QikgI8A3QBLsJpSDAeuAWoAq72t5+qPgQ8BDB27FjNzc0N6fx5eXmEum9DVmkB981dx9AxE+jSLiXsx29JmusamMDZNQg/VeWGj95n+pBMZs5ofBSyWLkG0VK8dhD/dzT13QH5uhzIBU5X1SdU9UNV/R1wPXCViIwMa6QRUluG+8G6Qo8jMcZ4Ye2eI+w5XNai6nMgepLOKpx6nbqGAqsb2Xc4cFBVN9ZZ/on7PKSJsXliSLd2ZKUnk2dJx5hWaZ7bVHpaC2kqXSukpCMio0XkRREpEpEqETneXX6XiJwawiFfASaIyBe9n0QkG6c59CuN7LsH6CgiA+osP8F93hlCPJ4TEabndOHDDYVUVtd4HY4xJsLy1hYyrHs6Wektq3g96KTjVtovBAYDT9U5Rg1wVQhx/BNndtKXRWSOiMwGXsYZWeBBn3P3cZPcLT77PgocAd4QkYtFZLqI/Bz4Hc6UDPkhxBMVcnO6cKSsiqVbGythNMa0JMWllSzZdpDcFnaXA6Hd6dwDzMUpDvtpnXVLgeODPaDbrHkGsB54HHgS2AzMUNUSn00FiMcnblXdAkwAluGMYvAGTmfTh4CTVTVmbxMmDehEfJwwf0OR16EYYyJo4cYiqmuUaYNaVn0OhNZ67XjgHFVVEdE664pwxj4LmjvG2rmNbLMFP509VXU1cEEo541m6SmJjO7VgfkbCvnZrByvwzHGRMiHG4pomxTP6N4dvA4l7EK50ykD6huPoRtQHHo4pq6pAzvz+c5iG3XamFZCVflwfSET+2eSGB8tbb3CJ5RXtAD4iYj4Tl9Xe8dzOc7AnSZMpg7KRBXyN1oRmzGtwdb9x9hxsJQTB2V6HUqzCCXp/BqniG25+7MCF4vIPJy6ldvCF54Z0aM96SkJzF9vSceY1mD+BqebxNSBLa8RAYSQdFR1OXAisBe4CaeOpXa4mWmqui584ZmE+DgmD8hk/oZCVOtWoRljWpoPNxTRs2Mbsju1jFGl6wqpwFBVl6rqSUA7oCeQrqrTVfWzsEZnAOcbz67iMjYWxtzYpcaYIFRW17Bw436mDuyMSMwNkB+QJtVSqWqZqu6KpflqYtHUgU7Z7oINNjqBMS3Z8u2HKCmv4sSBLbM+B0LrHLqpvvHMROQ4EdnU9LCMr14ZqWR3SrX+Osa0cB9uKCJOYFJ/Szq+soHketalAH1CjsbUa+rAzizctJ+Kqpjt62qMacT8DYWM7NWB9qmJXofSbEItXquvRnsszoydJsymDszkWEU1S2xIHGNapOJjlSzffqjFtlqrFdCIBCJyHXCd+6sCr4pI3d6KbXCmIvhv+MIztSb2d4bEyS8oYmL/Tl6HY4wJs4WbiqhRmDKg5RatQeDD4GwC3nN/vhhYDNSt1S7HmYbg4fCEZny1S0lkZM/25G8s4mfYkDjGtDT5Bftb7NA3vgJKOqr6Ms6oz7XN+G5X1c3NGJfxY/KATP46r4DDZZWkp7TcMl9jWqP8giLG981okUPf+Aqlc+illnC8MXlAJjUKH2/c73Uoxpgw2nWolE1FR5ncwovWILRRpgFwm03n4LRY+wpVfawpQRn/RvfuQJvEeD7auJ9ThnX1OhxjTJjkFzjdIaa04P45tYJOOiLSAXgdZ5w1+HKqAd8WbZZ0mkFyQjzj+mawoMD66xjTkuQXFJGZlkROVjuvQ2l2oRQe3gV0whl/TYCzcSZgexKnwcH4sEVnvmbKgE4U7CthT3GZ16EYY8JAVcnfuJ9J/TNb7NA3vkJJOrNwEs/H7u87VDVPVb8HvAtcG67gzNfVlvl+ZFMdGNMibNhXQuGR8hbfVLpWKEmnG7BJVatxJnTzvR98ETgjHIEZ/4Z0TSejbZIVsRnTQixwh7ea3ArqcyC0pLMH6OD+vBWY6LNuQFMDMg2LixMm9u9EfkGRTXVgTAuQX1BEdqdUenRo43UoERHqzKG1jQgeB34jIg+KyF+B+4C54QrO+DdlQCZ7D5fbVAfGxLjK6hoWbT7QKppK1wqlyfRtQHf35/twGhVcCKQCrwA/Dk9opj61Zb/5BUUM6JLmcTTGmFB9vsOZyqA1JZ1QOoduVNX57s+Vqnq9qvZU1QxV/baqWs/FZtYrI5VeGW2+aNtvjIlN+QX7EYGJ/VrPeIo2n06MmtQvk4837ae6xup1jIlV+QVFDO2WTse2SV6HEjE2n06MmjSgE4fLqli5s9jrUIwxISitqOazbYdaVdEa2Hw6Mat2ZsF8669jTEz6dMsBKqprmNTKpioJKOmIyHUisk1EtvHlfDrb6jwKgb8CbzVnwMbRuV0yOVnt+KjAqtCMiUX5G4tIjBfG983wOpSIsvl0YtikAZ14atE2yiqrSUmM9zocY0wQPirYz+heHUlNCnnc5Zhk8+nEsMn9M/l3/haWbjv4RXGbMSb6HTpWwcpdxVx70kCvQ4m4UOp0Lge2+y4QkVkicr2IjApLVCYgJ/TLID5OrIjNmBjz8ab9qNLqGhFAaEnnaeCR2l9E5CrgTZyOootEZGaYYjONaJeSyAh3CmtjTOzIL9hPalI8I3t28DqUiAsl6UwA3vD5/ec49TjtcQb8vCkMcZkATe6fyec7ijlSVul1KMaYAOVvdKamTkpo2VNT+xPKK+4C7AQQkQFAX+ABVT0C/BsYHkogItJLRJ4XkWIROSwiL4pI7yD2HyIiz4lIkYiUisg6EWnx0yxMGtCJ6hpl0aYDXodijAnAnuIyNhUeZXIrrYcNJekcxhlvDSAXKFLVz93fq/EzfXVjRCQVeB8YjNM67iJgIDBPRNoGsP9YYBFOp9XvA6cDvwdafJOu43t3JDkhzorYjIkRtcNXTRrQuvrn1Aqlrd5HwA0iUgX8hK8WtQ0AdoRwzCuAfkCOqhYAiMjnwAbgB8D99e0oInE402O/p6pn+6yaF0IcMSclMZ7xfTOsMYExMSJ/YxEZbZMY0jXd61A8Ecqdzi9w7nRewbmrudVn3YXAwhCOORv4uDbhALhNsvOBOY3smwsMoYHE1NJN6p/Jur1H2HfEprA2JpqpKvkFRUzs34m4uJY/NbU/oYwyvUFVBwKdVXWAqm7xWX0tTlIK1jBgpZ/lq4Chjew7xX1OEZGPRaRSRPaJyJ9FpFXMilQ71cHCjXa3Y0w021h4lL2HW8/U1P6E3BXW3xQGqroixMNlAAf9LD8AdGxk39q5fZ4BHgBuwBkD7nagF3C2v51E5ErgSoCsrCzy8vKCDhqgpKQk5H3DpUaVtonw/PyVtD+0wdNYvBAN16C1s2sQmHe3Oq1ME4oKyMsL74D8sXINWsL4C7V3a0+o6i3uz3kiEg/cIyJDVHVN3Z1U9SHgIYCxY8dqbm5uSCfPy8sj1H3D6cSdS1i+/RDTpk2rHTWi1YiWa9Ca2TUIzJOPLaZXxmHOP31G2I8dK9cgWhqJH8T/HU19d0C+au+43qmz/G33eXQT4ooZkwdksqu4jC37j3kdijHGj6rqGj7etL9VF61B9CSdVTj1OnUNxRlEtLF9G1ITUkQxpnY4jQU2m6gxUWnFzmKOlFW1+nESoyXpvAJMEJF+tQtEJBuY7K5ryJs4I1zPqrP8VPd5cZhijGrZnVLp0aEN+Rss6RgTjb7on9PK5s+pK9D5dA6IyPHuz4+ISN8wx/FPYAvwsojMEZHZOKNabwce9Imjj4hUiUht3U1tg4a7gatE5C4RmSkiNwC3AP/xbYbdkokIkwd0YqFNYW1MVMov2M/Qbul0Sqtv4uXWIdA7nbZ8OUX1JUDncAahqkeBGcB64HHgSWAzMENVS3w2FZxRBurGfTtOU+0LcDqrXo0zAOkV4Ywz2k0ekElxaSWrdtkU1sZEk9KKapZsPcjkVjoKga9AW69tBa4QkdrEM1pE6h3uRlU/DDYQVd0GnNvINltwEk/d5YrTObTVdhCFL6ewXlBQxIhWOHqtMdGqdmrq1jiVQV2BJp17cIq5LsaZrvpv9Wwn7voWP+ZZNOrcLpnBXduxYEMRP8wd4HU4xhjXgoLWOTW1P4HOHPqIiLwJDMIZ0+wa4Gt9X4z3pg7M5D8fbaW0opo2SZb7jYkG8zcUMbZPRqubmtqfgP8Cqrob2C0i/wFet+mqo9PUgZ355/zNLNq8n9ycLl6HY0yrV3iknDW7D/PzWTlehxIVQhl77VJLONGrdmKo+dZ02pioUNtU+sSBYW1/FbNC6qcjIt1E5Hci8qmIbHSffysiXcMdoAlOSmI847MzWGBJx5io8OGGQjqmJjKse+ucyqCuoJOOiAwCluHU65QAn7jP1wLLRGRgOAM0wZs60JnqYO9hm+rAGC+pKgs2FDF5QGarncqgrlDudO7FmT10kKpOV9Vvqep0nEYGxe5646EpA51mmVbEZoy31u8tYd+Rcita8xFK0pkO/LrOPDqo6lacCd2mNz0s0xRDuqaTmZbEgg2FXodiTKs23/0frP0iaEJLOknAkXrWHXHXGw/FxQlTBmSyoKCIGhsSxxjPzN9QRP/ObeneoVXMJxmQUJLOMuDHIvKVfcWZxOWH7nrjsSkDO1NUUsGaPYe9DsWYVqmssppFm/cz1YrWviKUnkq3A68Ba0TkGWA30BU4HxgInBG+8EyopvrU6wzr3t7jaIxpfZZsPUhZZc0X/4vGEUo/nbeAM3GK0m4C/grcjNOC7UxVfbuB3U2EZKWnkJPVjg/XW72OMV74cH0hifHCCf1skE9fIfXTUdW3VHUs0A7oBbRT1fGqOjes0Zkmyc3pzKdbDlBSXuV1KMa0OnnrChmXnUFasg1946tJk7ip6jFV3amqNkdyFJqW05nKauUjm03UmIjadaiUdXuPkJtj9Tl1RcvMoaYZjO2TQdukeOatsyI2YyIpz/2fs/EPv86STguWlBDHlIGZfLBuH86UQ8aYSMhbt48eHdowsEua16FEHUs6LVxuThd2FZexYV9J4xsbY5qsoqqG/IIipuV0xulJYnxZ0mnhasuU89bt8zgSY1qHxVsOcLSimtxBVp/jTygDfn4kIhf5TF1toli39m3IyWr3RRmzMaZ55blNpSfZ1NR+hXKnUwH8B9glIveLyOAwx2TCzJpOGxM5eev2WVPpBoTSOTQXGIqTeL4HrBKRPBG5UEQSwxyfCYPcnC5UVusXk0kZY5rHzkOlrN9bwnRrtVavUDuHrlXVnwI9gEuAeOApYIeI3CMi/cIXommqsdkdSUtOsCI2Y5pZbd2p9c+pX1M7h5ar6uM4E7jNBzoDvwDWi8hzNpNodEiMj2PqwEzeX7vXmk4b04zeX+M0lR5gTaXrFXLSEZE2InKZiHwCfAp0wUk+3YGrgUnAk2GJ0jTZzCFZ7D1czoqdxV6HYkyLdKyiigUFRZw8NMuaSjcg6JouERkO/AD4DtAWeBn4parO89nsnyKyB3guLFGaJpsxuAtxAu+u3suInh28DseYFmf+hiLKq2o4ZWiW16FEtVDudJYDZwF/BPqo6vl1Ek6tAmBh6KGZcOrYNomx2Rm8vXqv16EY0yK9s3ov6SkJjOub4XUoUS2UpHMeTrK5TVV317eRqq5RVZu6OoqcMjSLtXuOsP2Ajc9qTDhV1yjvr93H9MFdSIy3PvcNCaXJ9IuqWt0cwZjmNXOIc9v/7hq72zEmnJZuO8iBoxWcbEVrjQplRIJNIjKynnXHicimpodlmkN2ZlsGdkmzpGNMmL27ei+J8cI0G/qmUaHcB2YD9Q2BkwL0CTka0+xmDs1i0aYDFJdWeh2KMS3GO6v3MqFfJ9qlWP/4xoRa+FhfZ4+xwKEQj2ki4OShWVTVqA0AakyYbCwsYVPRUStaC1BATaZF5DrgOvdXBV4VkYo6m7UBMoD/hi88E26jenYgMy2Zd1bvZc6oHl6HY0zMe8dtEVpbZ2oaFuidzibgPfchwGKf32sfL+AkpitCCUREeonI8yJSLCKHReRFEekdwnFuEBEVkQWhxNHSxcUJJw/NYt7afZRVWnsQY5rqzZV7GN6jPd07tPE6lJgQ0J2Oqr6M0wm0tqftHaoatgYDIpIKvA+UAxfj3E3dCcwTkRGqejTA4/QDbgas7KgBZ47oxtOfbCNvXSGnHmcjFRkTqu0HjrF8+yFuPM0G2w9UKE2mLw1nwnFdAfQDzlLVl9wkNxunUcIPgjjO33GG3lkT5vhalBP6ZtCpbRKvfb7L61CMiWlvrHC6Kp4+vJvHkcSOgJKOiFSLyHj35xr39/oeoUzaMhv4WFULaheo6mYgH5gTYIzfBo4Hbgzh/K1KQnwcpx7XlffW7KO0worYjAnV6yt2M7JXB3plpHodSswIdOy124EdPj+He6jiYbjFd3WsAs5vbGcR6Qj8AfiFqh6wwfYad8aIbjy5aBvz1u2zb2nGhGDb/mN8vqOYm04f4nUoMSXQOp3bfH6+tRniyAAO+ll+AOgYwP73AeuBRwM9oYhcCVwJkJWVRV5eXqC7fkVJSUnI+3qpRpX0JPj3e8tJ3b/O63CaJFavQUvSGq/B65ucBrwdj24hL2+bx9HEzjWI+flURWQqzgymx2sQk8Wo6kPAQwBjx47V3NzckM6fl5dHqPt6bU7xSp5bsp3xk6aQmhS7b4VYvgYtRWu8Br9bMZ/Rvdty3mmTvQ4FiJ1rEGidTmP1OE2t0zmI/zua+u6AfD0I/Atn1tIOItIBJ5nGu7/XN3pCq3fGiG6UVdbw/lpr7GdMMLYUHWXlzsOcYUXTQQumTqc5p5xchVOvU9dQYHUj+w5xH1f5WXcQp+/QH5sSXEs1LjuDzu2SeW35bs4c0d3rcIyJGa9bq7WQBVqnc2szx/EK8DsR6VfbHFtEsoHJwA2N7Otv+oQ/AvHAj3Hm9TF+xMcJZwzvxlOfbKP4WCXtU23cKGMao6r877OdjO3T0TqEhiBaJn74J7AFeFlE5ojIbJzWbNtxis8AEJE+IlIlIrfULlPVvLoPnPHfit3fd2Dqdd6YnlRU1fDK8p1eh2JMTFi2/RAF+0o4b0xPr0OJSVGRdNwRB2bgtEB7HKeD52ZghqqW+GwqOHcwURF3SzCsezqDu7bj+SWWm40JxHNLdpCSGMcZI6xoLRSBDvhZDUxU1U9EpIaG63dUVYNuCqWq24BzG9lmC07iaexYucGev7USEc4b05M7X1/D+r1HGJTVzuuQjIlaZZXVvLp8F6cd182mMQhRtHQONR46e3QP7nlzLc8v2cGvrKObMfWau2oPR8qqON+K1kIWLZ1DjYc6pSUzY3AXXly6k5/PyrE53o2px/NLdtCjQxsm9OvkdSgxK+RPFxFJF5FJInK++2zlMjHsvDE9KSop54N1hV6HYkxU2nWolAUFRZw7pidxcTbUVqhCSjpu67HtwHzgGWABTufMm8MYm4mg6YO70KltkjUoMKYeLy7dgSqcd7wVrTVF0ElHRG4DbsVJNicDw4GZwLPAbSJyaxjjMxGSGB/HOcf34N01e9lTXOZ1OMZElarqGp7+ZDsT+mXQu5ONKN0UodzpXAH8XlWvVNX3VXWV+3wFzkjPV4Y3RBMpF03IplqVJz7e6nUoxkSVd9fsZeehUi6ZlO11KDEvlKTTHphbz7q33PUmBvXulMpJg7N46pNtNpW1MT4eyd9Cjw5tOHmozbTbVKEknUXAuHrWjXPXmxh12ZRsDhyt4JVlNquoMQCrdhXzyeYDXDIpm3hrQNBkgY4yHVf7AK4BLhORn4tItoi0cZ9/AVwG/F9zBmya18R+nRjctR2P5G8miJkijGmx/p2/hdSkeC4Y18vrUFqEQO90qoBK97Ec6A/cA2wEStznu93ln4c/TBMpIsKlk7NZu+cICzft9zocYzxVeKScV5bt4tzje9K+jY1AEA7RMrWBiSJzRjkjFPw7fwuT+md6HY4xnnlq0TYqqmu4ZHK216G0GNEytYGJIimJ8Xx3Qh8emFdg47GZVutoeRWPLdxCbk5n+ndO8zqcFsPGOzF+XTa5L22TErj/7fVeh2KMJ/6dv5n9Ryu49qSBXofSogQ9GnQtERkJ5AApddep6mNNCcp4r2PbJC6f0pc/vbeBFTuKGd7TWsKb1qP4WCUPfriJmUO6MLp3R6/DaVGCTjoi0gF4HZhQu8h99q3zsaTTAlw+tS//WbiF+99Zx78vHe91OMZEzMMLNnGkrIqfnpzjdSgtTijFa3cBnYATcRLO2TgTsD0JbALs06mFSE9J5Acn9mfeukKWbD3gdTjGRMT+knIeWbCZM0Z0Y2j3dK/DaXFCSTqzcBLPx+7vO9xpob8HvAtcG67gjPcuntSHzLRk7pu7zvrtmFbhHx9spLSymutmDvI6lBYplKTTDdikqtVAGeDbtOlF4IxwBGaiQ2pSAv83vT8fbzrAmyv3eB2OMc2qYN8R/vPRVs45vicDuliLteYQStLZA3Rwf94KTPRZN6CpAZno890JfRjWPZ3fvLKK4mOVXodjTLOoqVFueGEFqcnx3HDaYK/DabFCSToL+LIRwePAb0TkQRH5K3Af9Q8GamJUQnwc9547ggNHK7j7zTVeh2NMs3jyk20s3nqQm04fQmZastfhtFihNJm+Deju/nwfTqOCC4FU4BXgx+EJzUST43q05/tT+vLgh5uYM6oHE/vbdL2m5dhTXMa9b65l8oBOnDfGJmlrTkHf6ajqRlWd7/5cqarXq2pPVc1Q1W+rqg3Y1UL9ZOYgemek8qv/reBYRZXX4RgTFqrKzS+tpKqmhrvOHo6IjSTdnEIekUBE0kVkkoicLyITRcTGSmnh2iTFc885w9my/yjXP7ucmhprzWZi39/yNvLumr387JQc+nRq63U4LV5ISUdEbgG2A/Nxpq3OB3aIyM1hjM1EoUkDMvnVaUN4c+Ue/vTeBq/DMaZJ3lq5h/vmrmPOqO5cPqWv1+G0CqGMSHAb8GvgYeC/wF4gC/gWcJuIJNgAoS3b96f2Zf3eI/zpvQ0MzErjzBHdG9+pGagqew6X8dm+KtZ+sJFDxyopLq3gWEU17dsk0iE1ic5pSYzo2YGh3dNJjLehBiNBVdlYWMLSbYfYU1zGgaMVFJdWkhgvX1yXXhmpHNc9nexObYnzaGK01bsO89NnlzGyVwfuPXeEFatFSCgNCa4Afq+qP/dZtgp4X0SKgSuBW8MQm4lSIsKdZx/H5iKnmC09JZETB3WOyLn3HSkjb10heev2sWjTAfYfrXDXrCUpPo72qYm0SYynuLSSw2WV1PZnbZMYz/F9OjB7ZHfOHNGdtskhDzto/KiqUd5etYcXlu5g0eYDHPJpWt8uOYH2qYlUVSvFpZWU+kyF3jYpnuP7dCQ3pwszBnehb2ZkircK9pXw/f98SnpKIv+8aAwpifEROa8JLem0p/5m0W8BV4cejokVyQnx/OOiMXz34UVc+uin/L+zjuOb43s3y7mKSsp5Y8VuXl62iyVbDwKQlZ5Mbk4XRvRsT+XejVx42omkJSd85dtqdY2y93AZS7cdZPGWg8zfUMgvX1jB7a+u5hsju3Plif3oZ0PWN8mBoxU8PH8TTy4spbh8CV3aJXPK0CzG9slgTHZHenVMJSnhq3eYZZXVbCo8yspdxazcWUx+QRF3vLaaO15bzYAuacwe2Z3ZI7uT3UwJaNGm/Vz5+BIS44X/XDaeLulfG7PYNKNQks4iYBzOkDd1jXPXm1YgMy2Z566ayI+e+owbXlzBtgPH+NkpOWEpLimrrObt1Xt5cekO5m8oorpGGdy1HdefPIgZQ7owtFv6FwkmL28L7VK+PqtjfJzQvUMbundow5kjuqOqLN12iGc+3cbLy3bx/JIdfGt8b645aSCd21m/jGCUVlTzSP5m/pG3kaMVVYzsHM+PTh1Nbk5nEhopxkxJjGdo93SGdk/ngrHOFNDb9h/j/bV7eWPlHu5/Zz33v7Oekb06cO7xPfjGiO50bJsUlrhf+mwnv3j+c3pltOHRS8fTKyM1LMc1gZNAxtMSEd930VDgf8BDwHN8WadzAU7R2xxVXR3+UJvP2LFjdfHixSHtm5eXR25ubngDijGV1TXc8vJKnv5kO0O7pXPTGUOYPCD4GUcrq2tYuHE/r32+izdX7OFIeRXd26dw1ugezBnVg5yu/htIhnINCo+U86f31vP0J9tJSYjjZ7NyuHhitmf1C7HkvTV7uel/K9lzuIyZQ7L45ak57FyzJGz/B7sOlfLa57t4celO1u45QmK8MD2nC2eO7M7MIV1ITQr+u/KGvUe46401zFtXyIR+GTz43bG0T21Z00978VkkIktUdWww+wR69ar46tQFAtzjPqiz/PMgjmtagMT4OO46ezgT+nXit2+t4zsPL2J6Tme+Ob43kwdkktZA/cm+w2Us3LSfBRuKeHfNXg4eqyQtOYFZw7py7vE9mNCvU7Mkgs7tkrnzrOFcNrkvd7y2mtteXc3bq/Zy3/kj6NnRvv36c6SskjteW82zi3cwuGs7/vyt0YzvmwHAzjAOVNG9QxuuPLE/V57Yn9W7DvPC0h28unwXb6/eS0piHNNzujBlYCaT+meS3Sm13gYAFVU1LN56gFeW7eLZxdtpm5zAjacN5tLJfb9W5GciJ9DkcDtfTTrGfIWIMGdUD2YN68p/PtrCX+cVMG9dIQlxwpg+HemVkUpacgJtk+M5cLSC7QdK2XbgGNsOHAMgPSWBaTldOHNEN6YN6hyxit1+ndN45JJxPLt4O7e/uppT/zifu84ZzuyR3rTIi1afbTvI/z31GbuLS/nR9P5ce9KgiHxwO8VwQ/nV6UP4dMsBXv98N++s3vvF4LNd2iWT3aktvTJSyUpPpqyyhqPlVew7UsYnmw9wtKKaxHjhexOzueakgWSEqZjOhC6gpBOJJtAi0gv4A3Ayzh3Tu8BPVHVbI/uNxWkxdyLQGyjC6T90s6pubtagzdekJMbzg2n9uXRyX5ZsPcgH6wvJLygiv6CIkvIqjpZXOU1mO7ZheM/2fPuE3kzun8nQ7unEe1S0JSJcOK43k/pnct0zy7jm6c9YubOYX8zKabR+ojV45tNt/PqlVXRJT+a5qyYxpk/kZ9KMjxMm9OvEhH6duH3OMDYXHSV/434+23qQ7QePsaCgkMIj5bRJjKdtcgIdUhM5a3QPcnO6MLF/pwbvtk1kRcWVEJFU4H2gHLgY567qTmCeiIxQ1aMN7P5NYBjwZ5ym2z1w+hEtFpFRqrq9WYM3fiUlxDGxf6evjdGmqlHbH6JXRipPXTGBO15bzUMfbmL1rsP85Vujw1aJHWsqqmq447XVPP7xVqYOzOTP34yOv4WI0K9zGv06p3HRhD5fLI/m95b5UkhJR0S6AdcD04AM4AAwD7hfVUOZdOUKoB+Qo6oF7jk+BzYAPwDub2Dfe1W1sE58+cBm97i3hBCPaSbR/qGQlBDHHWcdx/Ae7bn5pZWc8/ePeOyy1tfKqaS8iqufWML8DUX84MR+/DwG7vqi/b1lHEG/i0RkELAMuAYoAT5xn68FlonIwBDimA18XJtwANyisXxgTkM71k047rKtQCHOXY8xQbtgXC+euuIEDhyt4Oy/fcTKncVehxQx+46UceGDC/lo437uO28EN54+JOoTjokdobyT7gUOA4NUdbqqfktVpwODgGJ3fbCGASv9LF+F00Q7KCIyBOgC2OQvJmRjszN44eqJJCfEceGDC5m/4Wvfb1qczUVHOedvH7G56Cj/ungs57v9aIwJl4D66XxlB5FDwFWq+l8/674F/E1Vg6ppFJEKnKK5G+osvxO4QVUDLgYUkQTgPWAITnHdwXq2uxKnAQJZWVlj/vvfr72cgJSUlJCWZr3avdTc1+BgWQ33Lylnd0kNV49KZkxWVFSFht32IzXc92kpCvx0TAp92wfegtD+D7znxTWYPn16s/XT8ZUEHKln3RF3vZceACYBZ9SXcABU9SGcDq6MHTtWQ+1UZZ1DvReJazDjxEouefQT/ra8mN+fP4SzRresktvPth3kmkc+oW2bFJ74/gn0D3J4IPs/8F6sXINQiteWAT+uM0oB4tTi/dBdH6yDgL+7owx3XUBE5B6cu5fLVPXtEOIwxq/2qYk8fvkJjM/O4Lpnl/Hkoq1ehxQ2H2/az3cfXkSH1CSe/cHEoBOOMcEI5U7nduA1YI2IPAPsBroC5wMDgTNCOOYqnHqduoYCAQ2pIyI3Ab8Efqyqj4cQgzENSktO4N+XjuOHTy7lpv+tpKKqhksnx/YcLPM3FHLFY4vp2TGVJ79/Alk2+KVpZqFMV/0WTmI5AtwE/BW4GacF25kh3mG8AkwQkX61C0QkG5jsrmuQiFyD06/nJlV9IITzGxOQlMR4/vHdMcwalsVtr67mwQ82eh1SyOat3cfl/1lMdqe2PHPlBEs4JiKCSjoikiQi/wNK3cqjdkAvoJ2qjlfV+qY8aMw/gS3AyyIyR0RmAy/jzE76oM/5+4hIlTtzae2ybwJ/xJlW4X0RmeDzCLrlmzGNSUqI44FvH883Rnbn7jfX8pf3NhBsgxyvzV21hysfX0xOVjuevmICndJslG0TGUEVr6lqhYjMBP7k/n4MONbUIFT1qIjMwBkG53GcYXDewxkGp8RnUwHi+WqyPNVdfqr78PUBkNvU+IypKzE+jj9eOIrEeOH376znWGU1v5iVExMdFF9etpOfPrucET3b8+il42nfpmWNtmyiWyh1OvnABCAvnIG4Y6yd28g2W3ASjO+yS4BLwhmLMYGIjxN+d95I2iTG8/e8jZRWVHPLmUOjenqEZz7dxg0vruCEvhk8fPE4G5PMRFwo77jrgZdEpAR4CachwVfKFlS1pumhGRP94uKEO886jjaJ8Ty8YDNHy6u4+5zhUdmD/5EFm7n9tdVMG9SZf3x3DG2SbIpmE3mhJJ0V7vOf3EddGuJxjYlJIsJNZwyhbXICf3pvAwePVfKXb42Omg91VeW3c9fx97yNnDqsK3/61iiSE6IjNtP6hNpkOrZqTY1pZiLCdScPolNaEr95ZRUX/WsRD188lg6p3vaVrqyu4cYXV/D8kh18+4Te3DHnOM+mkDAGQkg6kZhbx5hY9b2J2WSmJfOT/y7jvH8s5F8Xj6VPp7aexFJcWsk1T3/GB+sLuW7mIK45aUBMNHQwLVvIBc8iki4ik0TkfPfZ/wT2xrQypw/vxmOXj6eopJzZD+SzYENRxGPYWFjC2X/NJ7+giHvOGc61MwdawjFRIaSk4/aT2Y4zQ+czwAJgh4jcHMbYjIlZE/p14pUfTaFregrfe2QRD8/fFLG+PO+v3ctZD+RTXFrJU1dM4Jvje0fkvMYEIpT5dG4DbsVJNicDw4GZwLPAbSJyaxjjMyZm9e6Uygs/nMTMIVnc+foaLn30U/YdLmu285VWVHPLyyu57NHF9MpI5eX/m8z4vhnNdj5jQhHKnc4VwO9V9UpVfV9VV7nPV+B07rwyvCEaE7vSkhN48KIx3D5nGB9v2s+sP37I65/vDvtdz7Lthzjjz/N5bOFWLp/Slxd/OImeHVvXbKcmNoSSdNoD9Q1385a73hjjEhG+NzGb1348lV4ZqfzoqaV8+5+LWLGj6bOR7jh4jJ8+s4yz/5ZPaWU1T33/BH595lBSEq1JtIlOoTSZXgSMA971s26cu94YU8eALmm8cPUknlq0jT+9t4FvPLCAb4zszsUT+zCmT8egKvoL9pXw1KJtPLFoKwJcNa0/V+f2Jz3FhrQx0S2UpHMN8D8RqQKeA/YCWcAFwGXAHN+5dmx0AmO+lBgfx8WTsjnn+B48+MEmHsnfzKvLd9Evsy3nHN+DsdkZDO2e/rXkUVZZzbo9R1i2/RAvLdvJZ9sOER8nnDWqB9efMojuHdp49IqMCU4oSedz9/ke9+FL+HLEArDRCYzxq11KIj+blcPVuf15Y8Vunluyg9+9vf6L9T06tCHVHdGgqkbZduAY1TVOPdCgrDRuOt2ZvbRzOxsd2sQWG5HAGA+1TU7g/LG9OH9sL4pKylm5s5hVuw6zfu8RKqudQgIR4cwR3RjWPZ1h3dvTs2Mb63NjYpaNSGBMlMhMSyY3pwu5OV28DsWYZhN9Q+EaY4xpsSzpGGOMiRhLOsYYYyLGko4xxpiIsaRjjDEmYizpGGOMiRhLOsYYYyLGko4xxpiIkUhNLBXNRKQQ2Bri7plA5KeGNL7sGnjProH3vLgGfVS1czA7WNJpIhFZrKpjvY6jNbNr4D27Bt6LlWtgxWvGGGMixpKOMcaYiLGk03QPeR2AsWsQBewaeC8mroHV6RhjjIkYu9MxxhgTMZZ0jDHGRIwlnRCISC8ReV5EikXksIi8KCK9vY6rtRCRniLyFxFZKCLHRERFJNvruFoTETlPRF4Qka0iUioi60TkbhFp53VsrYWIzBKR90Vkj4iUi8gOEXlWRIZ6HVtDrE4nSCKSCiwHyoGbcabuvhNIBUao6lEPw2sVRCQXeAZYAsQDpwB9VXWLd1G1LiLyMbANeBnYAYwGbgXWApNUtca76FoHEfkWcDywCCgEegM3AL2A4aoaaof3ZmVJJ0gici1wP5CjqgXusr7ABuAXqnq/l/G1BiISV/uhJiLfB/6JJZ2IEpHOqlpYZ9n3gP8AJ6nq+95E1rqJSA5O4v+Zqv7e63j8seK14M0GPq5NOACquhnIB+Z4FlUrYt+ivVc34bg+dZ97RDIW8xX73ecqT6NogCWd4A0DVvpZvgqI6rJUY5rZNPd5jadRtDIiEi8iSSIyEHgQ2AM87XFY9UrwOoAYlAEc9LP8ANAxwrEYExVEpAdwO/Cuqi72Op5WZhEwxv25AJihqvs8jKdBdqdjjGkSEUnDaVBQBVzqcTit0UXABODbwGHgnWhuzWlJJ3gH8X9HU98dkDEtloi0AV4F+gGzVHWHxyG1Oqq6RlUXqerTwElAGk4rtqhkxWvBW4VTr1PXUGB1hGMxxjMikgg8D4wFTlbVFR6H1Oqp6iERKQAGeB1LfexOJ3ivABNEpF/tAvdWdrK7zpgWT0TigCeBGcBZqvqxxyEZQESygMHARq9jqY/10wmSiLTF6RxaypedQ+8A2uF0Di3xMLxWQ0TOc388CbgK+CFOB7lCVf3As8BaCRH5O87f/f8Br9VZvcOK2ZqfiPwPWAp8jlOXMwi4DugKjFfV9R6GVy9LOiFwh7z5A3AyIMB7wE+sc2LkiEh9b9wPVDU3krG0RiKyBehTz+rbVPXWyEXTOonIL4ELgP5AErAdyAPujubPIks6xhhjIsbqdIwxxkSMJR1jjDERY0nHGGNMxFjSMcYYEzGWdIwxxkSMJR1jjDERY0nHmGYiImeJyE/rLMt1p9fO9SYqY7xl/XSMaSYi8igwU1V7+ixLxx2nT1UPexWbMV6xAT+NiSA30dg4ZabVsuI1Y5qBe5dzMdDDLU5TEdnir3hNRPJEZIGInCoiy0SkVEQ+E5ETRCRBRO4Skd0ickBEHnXH//M9V6qI3Csim0Wkwn2+yR2U05ioYnc6xjSPO4DOwDhgtrusHGhfz/YDgPtwBtAsAX6LM2r5Kzj/p5cAQ9xt9gG/ABCRBGAuTpHdHcAKnAm9fo0zx9P1YX1VxjSRJR1jmoGqbhSRQqDCd9j/BhoQdAImqeomd7s4nNk4+6rqTHebuSJyInA+btIBvgVMAaap6ofusvdEBOA3InJvNE9dbFofu/02Jjqsr004rrXu89w6260FeoqbVYBTga3AR25RXIJ79/M2kIhz12NM1LA7HWOiQ92pzisaWJ4AxANVQBecKQYq6zlup3AFaEw4WNIxJrbtBzbjzKviz5bIhWJM4yzpGNN8yoE2zXyOt4BzgRJVXdvYxsZ4zZKOMc1nNZAhIlcDi4GyZjjHk8ClOI0Hfo8zlXoSzmySs4GzVPVYM5zXmJBY0jGm+TyMU5F/F9ABp8L/knCeQFUrRWQWcANwJdAXOApsBF7ny7ohY6KCDYNjjDEmYqzJtDHGmIixpGOMMSZiLOkYY4yJGEs6xhhjIsaSjjHGmIixpGOMMSZiLOkYY4yJGEs6xhhjIub/A2c+SDDT+QNXAAAAAElFTkSuQmCC\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, ClassicalRegister, 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",
    "# Other\n",
    "from copy import deepcopy\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": [
    "MY TROTTERIZATION GOES HERE -- START"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9156980e-d3a7-4494-8ad4-7ec15ca92e7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameterize variables to be evaluated at t1=pi/6 and t2=pi/12 later\n",
    "t1 = Parameter('t1')\n",
    "t2 = Parameter('t2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9ba5b4cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build the subcircuit for XX(t1)YY(t1)ZZ(t1)\n",
    "def xxyyzz_30():\n",
    "    \n",
    "    # IMPORTANT! This subcircuit is derived for t1=pi/6 (30°)\n",
    "    qc = QuantumCircuit(2, name='xyz_30')\n",
    "    \n",
    "    # ZZ\n",
    "    qc.cx(0,1)\n",
    "    qc.rz(2*t1, 1)\n",
    "    \n",
    "    # XXYY\n",
    "    qc.crx(4*t1,1,0)\n",
    "    qc.cx(0,1)\n",
    "    \n",
    "    # This subcircuit represents \"the building block\" of a trotter step.\n",
    "    return qc.to_gate()\n",
    "\n",
    "# Build the subcircuit for XX(t2)YY(t2)ZZ(t2)\n",
    "def xxyyzz_15():\n",
    "    \n",
    "    # IMPORTANT! This subcircuit is derived for t2=pi/12 (15°)\n",
    "    qc = QuantumCircuit(2, name='xyz_15')\n",
    "    \n",
    "    # ZZ\n",
    "    qc.cx(0,1)\n",
    "    qc.rz(2*t2, 1)\n",
    "    \n",
    "    # XXYY\n",
    "    qc.crx(4*t2,1,0)\n",
    "    qc.cx(0,1)\n",
    "    \n",
    "    # This subcircuit represents \"the building block\" of a trotter step.\n",
    "    return qc.to_gate()"
   ]
  },
  {
   "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": 10,
   "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",
    "\n",
    "#\n",
    "# Make subcircuit for the angle = 30°\n",
    "#\n",
    "Trot_qc_30 = QuantumCircuit(Trot_qr, name='Trot_30')\n",
    "\n",
    "for i in range(0, num_qubits - 1):\n",
    "    Trot_qc_30.append(xxyyzz_30(), [Trot_qr[i], Trot_qr[i+1]])\n",
    "\n",
    "# Convert custom quantum circuit into a gate\n",
    "Trot_gate_30 = Trot_qc_30.to_instruction()\n",
    "\n",
    "#\n",
    "# Make subcircuit for the angle = 15°\n",
    "#\n",
    "Trot_qc_15 = QuantumCircuit(Trot_qr, name='Trot_15')\n",
    "\n",
    "for i in range(0, num_qubits - 1):\n",
    "    Trot_qc_15.append(xxyyzz_15(), [Trot_qr[i], Trot_qr[i+1]])\n",
    "\n",
    "# Convert custom quantum circuit into a gate\n",
    "Trot_gate_15 = Trot_qc_15.to_instruction()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0641237c-048e-41e7-b47c-88a364683fb7",
   "metadata": {},
   "source": [
    "MY TROTTERIZATION GOES HERE -- FINISH"
   ]
  },
  {
   "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": 11,
   "id": "c6135540-dcea-4313-8ce8-0eb8e3ee23e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The final time of the state evolution\n",
    "target_time = np.pi\n",
    "\n",
    "# Number of 30° and 15° trotter steps. Their total number must match the target_time.\n",
    "trotter_steps_30 = 4  \n",
    "trotter_steps_15 = 4\n",
    "trotter_steps = trotter_steps_30 + trotter_steps_15\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)\n",
    "# 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_30):\n",
    "    qc.append(Trot_gate_30, [qr[1], qr[3], qr[5]])\n",
    "    \n",
    "# Simulate time evolution under H_heis3 Hamiltonian\n",
    "for _ in range(trotter_steps_15):\n",
    "    qc.append(Trot_gate_15, [qr[1], qr[3], qr[5]])\n",
    "    \n",
    "# Apply the basic error correction for qubits #1 and #5.\n",
    "# We chose these qubits based on the physical connections between qubits in the \"ibmq_jakarta\" quantum computer.\n",
    "qc.cx(1,0)\n",
    "qc.cx(1,2)\n",
    "qc.cx(5,4)\n",
    "qc.cx(5,6)\n",
    "\n",
    "# Define the parameters t1 and t2 to specific values\n",
    "qc = qc.bind_parameters({t1: target_time/6})  # t1 = pi/6  (30°)\n",
    "qc = qc.bind_parameters({t2: target_time/12}) # t2 = pi/12 (15°)\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": "code",
   "execution_count": 12,
   "id": "87336572",
   "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",
       "q1_0: ─────────────────────────────────────────────────────────────────────────»\n",
       "           ┌───────────────┐┌───────────────┐┌───────────────┐┌───────────────┐»\n",
       "q1_1: ─────┤0              ├┤0              ├┤0              ├┤0              ├»\n",
       "           │               ││               ││               ││               │»\n",
       "q1_2: ─────┤               ├┤               ├┤               ├┤               ├»\n",
       "      ┌───┐│               ││               ││               ││               │»\n",
       "q1_3: ┤ X ├┤1 Trot_30(π/6) ├┤1 Trot_30(π/6) ├┤1 Trot_30(π/6) ├┤1 Trot_30(π/6) ├»\n",
       "      └───┘│               ││               ││               ││               │»\n",
       "q1_4: ─────┤               ├┤               ├┤               ├┤               ├»\n",
       "      ┌───┐│               ││               ││               ││               │»\n",
       "q1_5: ┤ X ├┤2              ├┤2              ├┤2              ├┤2              ├»\n",
       "      └───┘└───────────────┘└───────────────┘└───────────────┘└───────────────┘»\n",
       "q1_6: ─────────────────────────────────────────────────────────────────────────»\n",
       "                                                                               »\n",
       "c0: 3/═════════════════════════════════════════════════════════════════════════»\n",
       "                                                                               »\n",
       "c1: 4/═════════════════════════════════════════════════════════════════════════»\n",
       "                                                                               »\n",
       "«                                                                              »\n",
       "«q1_0: ────────────────────────────────────────────────────────────────────────»\n",
       "«      ┌────────────────┐┌────────────────┐┌────────────────┐┌────────────────┐»\n",
       "«q1_1: ┤0               ├┤0               ├┤0               ├┤0               ├»\n",
       "«      │                ││                ││                ││                │»\n",
       "«q1_2: ┤                ├┤                ├┤                ├┤                ├»\n",
       "«      │                ││                ││                ││                │»\n",
       "«q1_3: ┤1 Trot_15(π/12) ├┤1 Trot_15(π/12) ├┤1 Trot_15(π/12) ├┤1 Trot_15(π/12) ├»\n",
       "«      │                ││                ││                ││                │»\n",
       "«q1_4: ┤                ├┤                ├┤                ├┤                ├»\n",
       "«      │                ││                ││                ││                │»\n",
       "«q1_5: ┤2               ├┤2               ├┤2               ├┤2               ├»\n",
       "«      └────────────────┘└────────────────┘└────────────────┘└────────────────┘»\n",
       "«q1_6: ────────────────────────────────────────────────────────────────────────»\n",
       "«                                                                              »\n",
       "«c0: 3/════════════════════════════════════════════════════════════════════════»\n",
       "«                                                                              »\n",
       "«c1: 4/════════════════════════════════════════════════════════════════════════»\n",
       "«                                                                              »\n",
       "«      ┌───┐      ░          ┌─┐         \n",
       "«q1_0: ┤ X ├──────░──────────┤M├─────────\n",
       "«      └─┬─┘      ░ ┌─┐      └╥┘         \n",
       "«q1_1: ──■────■───░─┤M├───────╫──────────\n",
       "«           ┌─┴─┐ ░ └╥┘       ║ ┌─┐      \n",
       "«q1_2: ─────┤ X ├─░──╫────────╫─┤M├──────\n",
       "«           └───┘ ░  ║ ┌─┐    ║ └╥┘      \n",
       "«q1_3: ───────────░──╫─┤M├────╫──╫───────\n",
       "«      ┌───┐      ░  ║ └╥┘    ║  ║ ┌─┐   \n",
       "«q1_4: ┤ X ├──────░──╫──╫─────╫──╫─┤M├───\n",
       "«      └─┬─┘      ░  ║  ║ ┌─┐ ║  ║ └╥┘   \n",
       "«q1_5: ──■────■───░──╫──╫─┤M├─╫──╫──╫────\n",
       "«           ┌─┴─┐ ░  ║  ║ └╥┘ ║  ║  ║ ┌─┐\n",
       "«q1_6: ─────┤ X ├─░──╫──╫──╫──╫──╫──╫─┤M├\n",
       "«           └───┘ ░  ║  ║  ║  ║  ║  ║ └╥┘\n",
       "«c0: 3/══════════════╩══╩══╩══╬══╬══╬══╬═\n",
       "«                    0  1  2  ║  ║  ║  ║ \n",
       "«c1: 4/═══════════════════════╩══╩══╩══╩═\n",
       "«                             0  1  2  3 </pre>"
      ],
      "text/plain": [
       "                                                                               »\n",
       "q1_0: ─────────────────────────────────────────────────────────────────────────»\n",
       "           ┌───────────────┐┌───────────────┐┌───────────────┐┌───────────────┐»\n",
       "q1_1: ─────┤0              ├┤0              ├┤0              ├┤0              ├»\n",
       "           │               ││               ││               ││               │»\n",
       "q1_2: ─────┤               ├┤               ├┤               ├┤               ├»\n",
       "      ┌───┐│               ││               ││               ││               │»\n",
       "q1_3: ┤ X ├┤1 Trot_30(π/6) ├┤1 Trot_30(π/6) ├┤1 Trot_30(π/6) ├┤1 Trot_30(π/6) ├»\n",
       "      └───┘│               ││               ││               ││               │»\n",
       "q1_4: ─────┤               ├┤               ├┤               ├┤               ├»\n",
       "      ┌───┐│               ││               ││               ││               │»\n",
       "q1_5: ┤ X ├┤2              ├┤2              ├┤2              ├┤2              ├»\n",
       "      └───┘└───────────────┘└───────────────┘└───────────────┘└───────────────┘»\n",
       "q1_6: ─────────────────────────────────────────────────────────────────────────»\n",
       "                                                                               »\n",
       "c0: 3/═════════════════════════════════════════════════════════════════════════»\n",
       "                                                                               »\n",
       "c1: 4/═════════════════════════════════════════════════════════════════════════»\n",
       "                                                                               »\n",
       "«                                                                              »\n",
       "«q1_0: ────────────────────────────────────────────────────────────────────────»\n",
       "«      ┌────────────────┐┌────────────────┐┌────────────────┐┌────────────────┐»\n",
       "«q1_1: ┤0               ├┤0               ├┤0               ├┤0               ├»\n",
       "«      │                ││                ││                ││                │»\n",
       "«q1_2: ┤                ├┤                ├┤                ├┤                ├»\n",
       "«      │                ││                ││                ││                │»\n",
       "«q1_3: ┤1 Trot_15(π/12) ├┤1 Trot_15(π/12) ├┤1 Trot_15(π/12) ├┤1 Trot_15(π/12) ├»\n",
       "«      │                ││                ││                ││                │»\n",
       "«q1_4: ┤                ├┤                ├┤                ├┤                ├»\n",
       "«      │                ││                ││                ││                │»\n",
       "«q1_5: ┤2               ├┤2               ├┤2               ├┤2               ├»\n",
       "«      └────────────────┘└────────────────┘└────────────────┘└────────────────┘»\n",
       "«q1_6: ────────────────────────────────────────────────────────────────────────»\n",
       "«                                                                              »\n",
       "«c0: 3/════════════════════════════════════════════════════════════════════════»\n",
       "«                                                                              »\n",
       "«c1: 4/════════════════════════════════════════════════════════════════════════»\n",
       "«                                                                              »\n",
       "«      ┌───┐      ░          ┌─┐         \n",
       "«q1_0: ┤ X ├──────░──────────┤M├─────────\n",
       "«      └─┬─┘      ░ ┌─┐      └╥┘         \n",
       "«q1_1: ──■────■───░─┤M├───────╫──────────\n",
       "«           ┌─┴─┐ ░ └╥┘       ║ ┌─┐      \n",
       "«q1_2: ─────┤ X ├─░──╫────────╫─┤M├──────\n",
       "«           └───┘ ░  ║ ┌─┐    ║ └╥┘      \n",
       "«q1_3: ───────────░──╫─┤M├────╫──╫───────\n",
       "«      ┌───┐      ░  ║ └╥┘    ║  ║ ┌─┐   \n",
       "«q1_4: ┤ X ├──────░──╫──╫─────╫──╫─┤M├───\n",
       "«      └─┬─┘      ░  ║  ║ ┌─┐ ║  ║ └╥┘   \n",
       "«q1_5: ──■────■───░──╫──╫─┤M├─╫──╫──╫────\n",
       "«           ┌─┴─┐ ░  ║  ║ └╥┘ ║  ║  ║ ┌─┐\n",
       "«q1_6: ─────┤ X ├─░──╫──╫──╫──╫──╫──╫─┤M├\n",
       "«           └───┘ ░  ║  ║  ║  ║  ║  ║ └╥┘\n",
       "«c0: 3/══════════════╩══╩══╩══╬══╬══╬══╬═\n",
       "«                    0  1  2  ║  ║  ║  ║ \n",
       "«c1: 4/═══════════════════════╩══╩══╩══╩═\n",
       "«                             0  1  2  3 "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Add the classical register. The register is used for the error correction.\n",
    "qst_orig = deepcopy(st_qcs)\n",
    "cr = ClassicalRegister(4)\n",
    "for qst_ghz_circ in st_qcs:\n",
    "    qst_ghz_circ.add_register(cr)\n",
    "    qst_ghz_circ.measure([qr[0], qr[2], qr[4], qr[6]], [cr[0], cr[1], cr[2], cr[3]])\n",
    "    \n",
    "st_qcs[-1].draw()"
   ]
  },
  {
   "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": 13,
   "id": "79c3ef21-c143-418a-8b50-f0ab39c5ae35",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Job ID 623c9d67a2f72d0d3edab3a0\n",
      "Job ID 623c9d6b8293e98d2e1e5b0e\n",
      "Job ID 623c9d6de32b4234b5ecd736\n",
      "Job ID 623c9d6fa2f72d1b81dab3a1\n",
      "Job ID 623c9d71538eba81e0611255\n",
      "Job ID 623c9d74d97bff2cef69455d\n",
      "Job ID 623c9d7674de0ed0ee85acfb\n",
      "Job ID 623c9d7919e689048ac80bea\n"
     ]
    }
   ],
   "source": [
    "shots = 8192\n",
    "# reps = 8\n",
    "# backend = sim_noisy_jakarta\n",
    "reps = 8\n",
    "backend = jakarta\n",
    "\n",
    "jobs = []\n",
    "st_qcs_t = transpile(st_qcs, backend=backend, optimization_level=3)\n",
    "for _ in range(reps):\n",
    "    # execute\n",
    "    job = execute(st_qcs_t, 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": 14,
   "id": "1b142880-af5b-4913-8f50-9248e5d1646f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r",
      "Job Status: job has successfully run\n",
      "\r",
      "Job Status: job has successfully run\n",
      "\r",
      "Job Status: job has successfully run\n",
      "\r",
      "Job Status: job has successfully run\n",
      "\r",
      "Job Status: job has successfully run\n",
      "\r",
      "Job Status: job has successfully run\n",
      "\r",
      "Job Status: job has successfully run\n",
      "\r",
      "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": 15,
   "id": "1c813957",
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# POST-PROCESSING\n",
    "#\n",
    "# Remove corrupted results from the statistics\n",
    "def result_filter(raw_results, new_results):\n",
    "\n",
    "    for resultidx, _ in enumerate(raw_results.results):\n",
    "        old_counts = raw_results.get_counts(resultidx)\n",
    "        new_counts = {}\n",
    "\n",
    "        # Change the size of the classical register\n",
    "        new_results.results[resultidx].header.creg_sizes = [new_results.results[resultidx].header.creg_sizes[0]]\n",
    "        new_results.results[resultidx].header.clbit_labels = new_results.results[resultidx].header.clbit_labels[0:3]\n",
    "        new_results.results[resultidx].header.memory_slots = 3\n",
    "\n",
    "        # Filter the results\n",
    "        for reg_key in old_counts:\n",
    "            reg_bits = reg_key.split(' ')\n",
    "            \n",
    "            # DON'T FORGET, REVERSE ORDER OF BITS! FROM LEFT TO RIGHT!\n",
    "            # Our case:\n",
    "            #            4. Classical bit [0][0] should be the same as the 3. Tomography bit [1][0]\n",
    "            #            3. Classical bit [0][1] should be the same as the 3. Tomography bit [1][0]\n",
    "            #            2. Classical bit [0][2] should be the same as the 1. Tomography bit [1][2]\n",
    "            #            1. Classical bit [0][3] should be the same as the 1. Tomography bit [1][2]\n",
    "            if reg_bits[0][0]==reg_bits[1][0] and reg_bits[0][1]==reg_bits[1][0] and reg_bits[0][2]==reg_bits[1][2] and reg_bits[0][3]==reg_bits[1][2]:    \n",
    "                new_counts[reg_bits[1]]=old_counts[reg_key]\n",
    "\n",
    "        new_results.results[resultidx].data.counts = new_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f87c0c22-1f37-4f56-9ec5-7e3c5d7a8384",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "state tomography fidelity = 0.4367 ± 0.0073\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_qc):\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_qc)\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",
    "    \n",
    "    # Implement the post-processing of results\n",
    "    current_results = job.result()\n",
    "    filtered_results = deepcopy(current_results)    \n",
    "    result_filter(current_results, filtered_results)\n",
    "\n",
    "    # Apply the state tomography with the filtered results\n",
    "    fid = state_tomo(filtered_results, qst_orig)\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": 2,
   "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><th>System information</th></tr><tr><td>Python</td><td>3.8.8 (default, Apr 13 2021, 15:08:03) [MSC v.1916 64 bit (AMD64)]</td></tr><tr><td>OS</td><td>Windows</td></tr><tr><td>CPUs</td><td>4</td></tr><tr><td>Memory (Gb)</td><td>7.799663543701172</td></tr><tr><td colspan='2'>Fri Mar 25 16:50:55 2022 Střední Evropa (běžný čas)</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "%qiskit_version_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0b1e5c2",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.8"
  },
  "toc-autonumbering": false,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}