qosf-simulator-task.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Task: Quantum Circuit Simulator\n",
    "\n",
    "The goal here is to implement simple quantum circuit simulator.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "Before we start coding:\n",
    "\n",
    "\n",
    "### Qubit\n",
    "\n",
    "Qubit is the basic unit of quantum information. It is a two-state (or two-level) quantum-mechanical system, and can be represented by a linear superposition of its two orthonormal basis states (or basis vectors). The vector representation of a single qubit is: ${\\vert a\\rangle =v_{0} \\vert 0\\rangle +v_{1} \\vert 1\\rangle \\rightarrow {\\begin{bmatrix}v_{0}\\\\v_{1}\\end{bmatrix}}}$,\n",
    "Here, ${\\displaystyle v_{0}}v_{0}$ and ${\\displaystyle v_{1}}v_{1}$ are the complex probability amplitudes of the qubit. These values determine the probability of measuring a 0 or a 1, when measuring the state of the qubit.\n",
    "\n",
    "Code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Qubit in |0> state (100% probability of measuring 0)\n",
    "q0 = [1, 0]\n",
    "\n",
    "# Qubit in |1> state (100% probability of measuring 1)\n",
    "q1 = [0, 1] \n",
    "\n",
    "# Qubit |+> state (superposition: 50% probability of measuring 0 and 50% probability of measuring 1)\n",
    "q2 = [0.7071067811865475, 0.7071067811865475]\n",
    "\n",
    "# Qubit |-> state (superposition: 50% probability of measuring 0 and 50% probability of measuring 1) with phase pi\n",
    "q3 = [0.7071067811865475, -0.7071067811865475]\n",
    "\n",
    "# Qubit |i> state (superposition: 50% probability of measuring 0 and 50% probability of measuring 1) with phase pi/2\n",
    "q3 = [0.7071067811865475, 0+0.7071067811865475j]\n",
    "\n",
    "# Qubit |-i> state (superposition: 50% probability of measuring 0 and 50% probability of measuring 1) with phase -pi/2\n",
    "q4 = [0.7071067811865475, 0-0.7071067811865475j]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that vector contains probability amplitudes - not probabilities. Probability amplitude is complex number and can be negative. Probability is calculated as absolute value squared:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.5 0.5]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "q4 = np.array([0.7071067811865475+0j, 0-0.7071067811865475j])\n",
    "p4 = np.abs(q4)**2\n",
    "print(p4)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### State vector\n",
    "\n",
    "The combined state of multiple qubits is the tensor product of their states (vectors). The tensor product is denoted by the symbol ${\\displaystyle \\otimes }$.\n",
    "\n",
    "The vector representation of two qubits is:\n",
    "\n",
    "${\\displaystyle \\vert ab\\rangle =\\vert a\\rangle \\otimes \\vert b\\rangle =v_{00}\\vert 00\\rangle +v_{01}\\vert 01\\rangle +v_{10}\\vert 10\\rangle +v_{11}\\vert 11\\rangle \\rightarrow {\\begin{bmatrix}v_{00}\\\\v_{01}\\\\v_{10}\\\\v_{11}\\end{bmatrix}}}$\n",
    "\n",
    "Example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 0 0]\n"
     ]
    }
   ],
   "source": [
    "# Qubit in |0> state (100% probability of measuring 0)\n",
    "q0 = [1, 0]\n",
    "\n",
    "# Qubit in |1> state (100% probability of measuring 1)\n",
    "q1 = [0, 1] \n",
    "\n",
    "combined_state = np.kron(q0, q1)\n",
    "\n",
    "print(combined_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, what this vector tells us?\n",
    "\n",
    "It will be more clear if we write vector elements in a column with index expressed in binary format:\n",
    "\n",
    "```\n",
    "Index (dec)  Index (bin)  Amplitude  Probability\n",
    "================================================\n",
    "0            00           0          0 (  0%)\n",
    "1            01           1          1 (100%)\n",
    "2            10           0          0 (  0%)\n",
    "3            11           0          0 (  0%)\n",
    "```\n",
    "\n",
    "- First element (binary: 00) is probability of measuring 0 on both qubits.\n",
    "- Second element (binary: 01) is probability of measuring 0 on first qubit and 1 on second qubit.\n",
    "- Third element (binary: 10) is probability of measuring 1 on first qubit and 0 on second qubit.\n",
    "- Fourth element (binary: 11) is probability of measuring 1 on both qubits.\n",
    "\n",
    "#### Endianness\n",
    "\n",
    "It is important to say that different quantum programming frameworks use different orientation of bitstrings (endianness). In previous example, left bit belongs to first qubit and right bit belongs to second qubit. This enconding is called \"big endian\".\n",
    "\n",
    "But, in some frameworks (like Qiskit), encoding is opposite: rightmost bit belongs to first qubit and leftmost bit belongs to last qubit. This is called \"little endian\".\n",
    "\n",
    "So, vector from our example in Qiskit's \"little endian\" encoding will look like this:\n",
    "\n",
    "```\n",
    "Index (dec)  Index (bin)  Amplitude  Probability\n",
    "================================================\n",
    "0            00           0          0 (  0%)\n",
    "1            01           0          0 (  0%)\n",
    "2            10           1          1 (100%)\n",
    "3            11           0          0 (  0%)\n",
    "```\n",
    "\n",
    "\"Little endian\" encoding:\n",
    "\n",
    "- First element (binary: 00) is probability of measuring 0 on both qubits.\n",
    "- Second element (binary: 01) is probability of measuring 0 on second qubit and 1 on first qubit.\n",
    "- Third element (binary: 10) is probability of measuring 1 on second qubit and 0 on first qubit.\n",
    "- Fourth element (binary: 11) is probability of measuring 1 on both qubits.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quantum gates\n",
    "\n",
    "Quantum gates are basic units of quantum processing. Gates are represented as unitary matrices. The action of the gate on a specific quantum state is found by multiplying the vector ${\\displaystyle \\vert \\psi _{1}\\rangle }$  which represents the state, by the matrix ${\\displaystyle U}$ representing the gate. The result is a new quantum state ${\\displaystyle \\vert \\psi _{2}\\rangle }$ \n",
    "\n",
    "${\\displaystyle U\\vert \\psi _{1}\\rangle =\\vert \\psi _{2}\\rangle }$\n",
    "\n",
    "Quantum gates (usually) act on small number of qubits. We have single-qubit and multi-qubit gates. n-qubit gate is represented as $2^n\\times2^n$ unitary matrix.\n",
    "\n",
    "Examples:\n",
    "\n",
    "#### Single qubit gates\n",
    "\n",
    "X (aka NOT) gate:\n",
    "\n",
    "$X = \\begin{bmatrix} 0 & 1 \\\\ 1 & 0 \\end{bmatrix}$\n",
    "\n",
    "Hadamard gate:\n",
    "\n",
    "$H = \\tfrac{1}{\\sqrt{2}}\\begin{bmatrix} 1 & 1 \\\\ 1 & -1 \\end{bmatrix}$\n",
    "\n",
    "General single qubit rotation gate:\n",
    "\n",
    "$U_3(\\theta, \\phi, \\lambda) = \\begin{bmatrix} \\cos(\\theta/2) & -e^{i\\lambda}\\sin(\\theta/2) \\\\\n",
    "            e^{i\\phi}\\sin(\\theta/2) & e^{i\\lambda+i\\phi}\\cos(\\theta/2)\n",
    "     \\end{bmatrix}$\n",
    "\n",
    "#### Two-qubit gates:\n",
    "\n",
    "Controlled-X (aka CNOT) gate:\n",
    "\n",
    "${CNOT} = \\begin{bmatrix} 1 & 0 & 0 & 0 \\\\\n",
    "                         0 & 1 & 0 & 0 \\\\\n",
    "                         0 & 0 & 0 & 1 \\\\\n",
    "                         0 & 0 & 1 & 0 \\\\\n",
    "         \\end{bmatrix}$\n",
    "\n",
    "SWAP gate:\n",
    "\n",
    "${SWAP} = \\begin{bmatrix} 1 & 0 & 0 & 0 \\\\\n",
    "                         0 & 0 & 1 & 0 \\\\\n",
    "                         0 & 1 & 0 & 0 \\\\\n",
    "                         0 & 0 & 0 & 1 \\\\\n",
    "         \\end{bmatrix}$\n",
    "\n",
    "#### Examples\n",
    "\n",
    "Let's see how single-qubit gate modifies state of the qubit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial state:\t [1 0]\n",
      "Final state:\t [0 1]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Let's start with qubit in state |0> (100% probability of measuring 0)   \n",
    "\n",
    "q0 = np.array([1, 0])\n",
    "\n",
    "print(\"Initial state:\\t\", q0)\n",
    "\n",
    "# Define X (NOT) gate:\n",
    "\n",
    "X = np.array([\n",
    "[0, 1],\n",
    "[1, 0]\n",
    "])\n",
    "\n",
    "# Now apply X gate to a qubit (matrix-vector dot product):\n",
    "\n",
    "q0 = np.dot(X, q0)\n",
    "\n",
    "print(\"Final state:\\t\", q0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After applying X gate, qubit flips from state $|0\\rangle$ to state $|1\\rangle$.\n",
    "\n",
    "Now, let's see how Hadamard gate works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial state:\t [1 0]\n",
      "Final state:\t [0.70710678 0.70710678]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Let's start with qubit in state |0> (100% probability of measuring 0)\n",
    "    \n",
    "q0 = np.array([1, 0])\n",
    "\n",
    "print(\"Initial state:\\t\", q0)\n",
    "\n",
    "# Define H (Hadamard) gate:\n",
    "\n",
    "H = np.array([\n",
    "[1/np.sqrt(2), 1/np.sqrt(2)],\n",
    "[1/np.sqrt(2), -1/np.sqrt(2)]\n",
    "])\n",
    "\n",
    "# Now apply H gate to a qubit (matrix-vector dot product):\n",
    "\n",
    "q0 = np.dot(H, q0)\n",
    "\n",
    "print(\"Final state:\\t\", q0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After applying Hadamard gate on qubit in state $|0\\rangle$ it evolves to state $|+\\rangle$ which is equal superposition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix operator\n",
    "\n",
    "Quantum program eveloves quantum state by multiplying state vector with each gate's unitary matrix (dot product). Note that dimension of the state vector and dimension of the unitary matrix describing a gate usually don't match. For example: 3-qubit quantum circuit's state vector has $2^n=2^3=8$ elements, but single-qubit gate has $2^n\\times2^n=2^1\\times2^1=2\\times2$ elements. In order to perform matrix-vector multiplication, we need to \"resize\" gate's matrix to the dimension of the state vector. Let's call that matrix a **matrix operator**.\n",
    "\n",
    "Note that size of the matrix operator is $2^n\\times2^n$ where $n$ is total number of qubits in the circuit, so storing it into memory and calculating it requires a lot of memory and cpu power for bigger circuits. Optimizing this code is most interesting and challenging part, but for our purpose it is enough if you make it work smoothly with 8 qubits (the more - the better).\n",
    "\n",
    "#### Matrix operator for single-qubit gates\n",
    "\n",
    "Matrix operator for single-qubit gate can be calculated by performing tensor product of gate's unitary matrix and $2\\times2$ identity matrices in correct order.\n",
    "\n",
    "Example for single-qubit gate $U$ in 3-qubit circuit:\n",
    "\n",
    "- gate on qubit 0: ${O=U \\otimes I \\otimes I}$\n",
    "- gate on qubit 1: ${O=I \\otimes U \\otimes I}$\n",
    "- gate on qubit 2: ${O=I \\otimes I \\otimes U}$\n",
    "\n",
    "Example matrix operator for X gate acting on third qubit in 3-qubit circuit can be calculated like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial state: [1, 0, 0, 0, 0, 0, 0, 0]\n",
      "\n",
      "Operator:\n",
      "\n",
      " [[0. 1. 0. 0. 0. 0. 0. 0.]\n",
      " [1. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 0. 0.]\n",
      " [0. 0. 1. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 0. 0.]\n",
      " [0. 0. 0. 0. 1. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 1.]\n",
      " [0. 0. 0. 0. 0. 0. 1. 0.]] \n",
      "\n",
      "Final state: [0. 1. 0. 0. 0. 0. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Let's define state vector of the 3-qubit circuit in \"ground state\" (all qubits in state |0>)\n",
    "\n",
    "psi = [1, 0, 0, 0, 0, 0, 0, 0]\n",
    "print(\"Initial state:\", psi)\n",
    "\n",
    "\n",
    "# Define X (NOT) gate:\n",
    "\n",
    "X = np.array([\n",
    "[0, 1],\n",
    "[1, 0]\n",
    "])\n",
    "\n",
    "# Define 2x2 identity\n",
    "\n",
    "I = np.identity(2)\n",
    "\n",
    "# Calculate operator for X gate acting on third qubit in 3-qubit circuit\n",
    "\n",
    "O = np.kron(np.kron(I, I), X)\n",
    "\n",
    "print(\"\\nOperator:\\n\\n\", O, \"\\n\")\n",
    "\n",
    "\n",
    "# And finally, apply operator\n",
    "\n",
    "psi = np.dot(psi, O)\n",
    "print(\"Final state:\", psi)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are dealing with \"big endian\" encoding, so this result is correct: third qubit is flipped to state $|1\\rangle$  and other qubits are not changed.\n",
    "\n",
    "**Note**: if we want vector in \"little endian\" encoding (like Qiskit), then order in which we perform tensor product to calculate operator is opposite. Instead ${O=I \\otimes I \\otimes U}$ we would do ${O=U \\otimes I \\otimes I}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Matrix operator for multi-qubit gates\n",
    "\n",
    "If we want to apply two qubit gate on subsequent qubits ( 0-1, 1-2, 2-3 etc.) then we can use the same technique like we do with single qubit gates:\n",
    "\n",
    "For 3-qubit circuit, CNOT gate:\n",
    "\n",
    "- acting on first and second qubit, operator is ${O=CNOT \\otimes I}$\n",
    "\n",
    "- acting on second and third qubit, operator is ${O=I \\otimes CNOT}$\n",
    "\n",
    "But, multi-qubit gates can be applied to qubits which are not consequent, so this is not that trivial.\n",
    "\n",
    "The main feature of a controlled-$U$ operation, for any unitary $U$, is that it (coherently) performs an operation on some qubits depending on the value of some single qubit. The way that we can write this explicitly algebraically (with the control on the first qubit) is:\n",
    "\n",
    "$\\mathit{CU} \\;=\\; \\vert{0}\\rangle\\!\\langle{0}\\vert \\!\\otimes\\! \\mathbf 1 \\,+\\, \\vert{1}\\rangle\\!\\langle{1}\\vert \\!\\otimes\\! U$\n",
    "\n",
    "where ${\\mathbf 1}$ is an identity matrix of the same dimension as $U$. Here, ${\\ket{0}\\!\\bra{0}}$ and ${\\ket{1}\\!\\bra{1}}$ are projectors onto the states ${\\ket{0}}$ and ${\\ket{1}}$ of the control qubit &mdash; but we are not using them here as elements of a measurement, but to describe the effect on the other qubits depending on one or the other subspace of the state-space of the first qubit.\n",
    "\n",
    "We can use this to derive the matrix for the gate ${\\mathit{CX}_{1,3}}$ which performs an $X$ operation on qubit 3, coherently conditioned on the state of qubit 1, by thinking of this as a controlled-${(\\mathbf 1_2 \\!\\otimes\\! X)}$ operation on qubits 2 and 3:\n",
    "\n",
    "$\\begin{aligned}\n",
    "\\mathit{CX}_{1,3} \\;&=\\;\n",
    "\\vert{0}\\rangle\\!\\langle{0}\\vert \\otimes \\mathbf 1_4 \\,+\\, \\vert{1}\\rangle\\!\\langle{1}\\vert \\otimes (\\mathbf 1_2 \\otimes X)\n",
    "\\\\[1ex]&=\\;\n",
    "\\begin{bmatrix}\n",
    "  \\mathbf 1_4 & \\mathbf 0_4  \\\\\n",
    "  \\mathbf 0_4 & (\\mathbf 1_2 \\!\\otimes\\! X) \n",
    "\\end{bmatrix}\n",
    "\\;=\\;\n",
    "\\begin{bmatrix}\n",
    "  \\mathbf 1_2 & \\mathbf 0_2 & \\mathbf 0_2 & \\mathbf 0_2 \\\\\n",
    "  \\mathbf 0_2 & \\mathbf 1_2 & \\mathbf 0_2 & \\mathbf 0_2 \\\\\n",
    "  \\mathbf 0_2 & \\mathbf 0_2 & X & \\mathbf 0_2 \\\\\n",
    "  \\mathbf 0_2 & \\mathbf 0_2 & \\mathbf 0_2 & X\n",
    "\\end{bmatrix},\n",
    "\\end{aligned}$\n",
    "\n",
    "where the latter two are block matrix representations to save on space (and sanity).\n",
    "\n",
    "Better still: we can recognise that &mdash; on some mathematical level where we allow ourselves to realise that the order of the tensor factors doesn't have to be in some fixed order &mdash; the control and the target of the operation can be on any two tensor factors, and that we can fill in the description of the operator on all of the other qubits with $\\mathbf 1_2$. This would allow us to jump straight to the representation\n",
    "\n",
    "$\\begin{aligned}\n",
    "\\mathit{CX}_{1,3} \\;&=&\\;\n",
    "\\underbrace{\\vert{0}\\rangle\\!\\langle{0}\\vert}_{\\text{control}} \\otimes \\underbrace{\\;\\mathbf 1_2\\;}_{\\!\\!\\!\\!\\text{uninvolved}\\!\\!\\!\\!} \\otimes \\underbrace{\\;\\mathbf 1_2\\;}_{\\!\\!\\!\\!\\text{target}\\!\\!\\!\\!}\n",
    "&+\\,\n",
    "\\underbrace{\\vert{1}\\rangle\\!\\langle{1}\\vert}_{\\text{control}} \\otimes \\underbrace{\\;\\mathbf 1_2\\;}_{\\!\\!\\!\\!\\text{uninvolved}\\!\\!\\!\\!} \\otimes \\underbrace{\\; X\\;}_{\\!\\!\\!\\!\\text{target}\\!\\!\\!\\!}\n",
    "\\\\[1ex]&=&\\;\n",
    "\\begin{bmatrix}\n",
    "  \\mathbf 1_2 & \\mathbf 0_2 & \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} \\\\\n",
    "  \\mathbf 0_2 & \\mathbf 1_2 & \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} \\\\\n",
    "  \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} \\\\\n",
    "  \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2}\n",
    "\\end{bmatrix}\n",
    "\\,&+\\,\n",
    "\\begin{bmatrix}\n",
    "  \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} \\\\\n",
    "  \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} \\\\\n",
    "  \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} & X & \\mathbf 0_2 \\\\\n",
    "  \\phantom{\\mathbf 0_2} & \\phantom{\\mathbf 0_2} & {\\mathbf 0_2} & X\n",
    "\\end{bmatrix}\n",
    "\\end{aligned}$\n",
    "\n",
    "and also allows us to immediately see what to do if the roles of control and target are reversed:\n",
    "\n",
    "$\\begin{aligned}\n",
    "\\mathit{CX}_{3,1} \\;&=&\\;\n",
    "\\underbrace{\\;\\mathbf 1_2\\;}_{\\!\\!\\!\\!\\text{target}\\!\\!\\!\\!} \\otimes \\underbrace{\\;\\mathbf 1_2\\;}_{\\!\\!\\!\\!\\text{uninvolved}\\!\\!\\!\\!} \\otimes \\underbrace{\\vert{0}\\rangle\\!\\langle{0}\\vert}_{\\text{control}}\n",
    "\\,&+\\,\n",
    "\\underbrace{\\;X\\;}_{\\!\\!\\!\\!\\text{target}\\!\\!\\!\\!} \\otimes \\underbrace{\\;\\mathbf 1_2\\;}_{\\!\\!\\!\\!\\text{uninvolved}\\!\\!\\!\\!} \\otimes \\underbrace{\\vert{1}\\rangle\\!\\langle{1}\\vert}_{\\text{control}}\n",
    "\\\\[1ex]&=&\\;\n",
    "{\\scriptstyle\\begin{bmatrix}\n",
    " \\!\\vert{0}\\rangle\\!\\langle{0}\\vert\\!\\! & & & \\\\\n",
    " & \\!\\!\\vert{0}\\rangle\\!\\langle{0}\\vert\\!\\! & & \\\\\n",
    "& & \\!\\!\\vert{0}\\rangle\\!\\langle{0}\\vert\\!\\! & \\\\\n",
    "& & & \\!\\!\\vert{0}\\rangle\\!\\langle{0}\\vert\n",
    "\\end{bmatrix}}\n",
    "\\,&+\\,\n",
    "{\\scriptstyle\\begin{bmatrix}\n",
    " & & \\!\\!\\vert{1}\\rangle\\!\\langle{1}\\vert\\!\\! & \\\\\n",
    " & & & \\!\\!\\vert{1}\\rangle\\!\\langle{1}\\vert \\\\\n",
    "\\!\\vert{1}\\rangle\\!\\langle{1}\\vert\\!\\! & & &  \\\\\n",
    "& \\!\\!\\vert{1}\\rangle\\!\\langle{1}\\vert & &\n",
    "\\end{bmatrix}}\n",
    "\\\\[1ex]&=&\\;\n",
    "\\left[{\\scriptstyle\\begin{matrix}\n",
    "1 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 1 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 1 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 1\n",
    "\\end{matrix}}\\right.\\,\\,&\\,\\,\\left.{\\scriptstyle\\begin{matrix}\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 1 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 1 \\\\\n",
    "1 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 1 & 0 \\\\\n",
    "0 & 0 & 0 & 0\n",
    "\\end{matrix}}\\right].\n",
    "\\end{aligned}$\n",
    "\n",
    "But best of all: if you can write down these operators algebraically, you can take the first steps towards dispensing with the giant matrices entirely, instead  reasoning about these operators algebraically using expressions such as $\\mathit{CX}_{1,3} =\n",
    "\\vert{0}\\rangle\\!\\langle{0}\\vert \\! \\otimes\\!\\mathbf 1_2\\! \\otimes\\! \\mathbf 1_2 +\n",
    " \\vert{1}\\rangle\\!\\langle{1}\\vert \\! \\otimes\\!  \\mathbf 1_2 \\! \\otimes\\!  X$\n",
    "and\n",
    "$\\mathit{CX}_{3,1} =\n",
    "\\mathbf 1_2 \\! \\otimes\\! \\mathbf 1_2 \\! \\otimes \\! \\vert{0}\\rangle\\!\\langle{0}\\vert +\n",
    "X \\! \\otimes\\! \\mathbf 1_2 \\! \\otimes \\! \\vert{1}\\rangle\\!\\langle{1}\\vert$.\n",
    "\n",
    "\n",
    "For Example, let's calculate operator for Controlled-X (CNOT) on first qubit as control and third qubit as target in 3 qubit quantum circuit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CNOT(0, 2) for 3-qubit circuit, operator is:\n",
      "\n",
      "[[1. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 1. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 1. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 0. 0.]\n",
      " [0. 0. 0. 0. 1. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 1.]\n",
      " [0. 0. 0. 0. 0. 0. 1. 0.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Define X gate (CNOT is controlled-X):\n",
    "\n",
    "X = np.array([\n",
    "[0, 1],\n",
    "[1, 0]\n",
    "])\n",
    "\n",
    "# Define 2x2 Identity\n",
    "\n",
    "I = np.identity(2)\n",
    "\n",
    "\n",
    "# Define projection operator |0><0|\n",
    "\n",
    "P0x0 = np.array([\n",
    "[1, 0],\n",
    "[0, 0]\n",
    "])\n",
    "\n",
    "# Define projection operator |1><1|\n",
    "\n",
    "P1x1 = np.array([\n",
    "[0, 0],\n",
    "[0, 1]\n",
    "])\n",
    "\n",
    "# And now calculate our operator:\n",
    "\n",
    "O = np.kron(np.kron(P0x0, I), I) + np.kron(np.kron(P1x1, I), X)\n",
    "\n",
    "print(\"CNOT(0, 2) for 3-qubit circuit, operator is:\\n\")\n",
    "print(O)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to implement simulator, it is best if you have function which returns operator for any unitary targeting any qubit(s) for any circuit size, something like:\n",
    "\n",
    "```\n",
    "get_operator(total_qubits, gate_unitary, target_qubits)\n",
    "```\n",
    "\n",
    "But this is not trivial so **it is enough if you can implement it for any 1-qubit gates and CNOT only.**\n",
    "\n",
    "If you are still enthusiastic and you wish to implement universal operator function then please refer to:\n",
    "\n",
    "- [qosf-simulator-task-additional-info.pdf](https://github.com/quantastica/qosf-mentorship/blob/master/qosf-simulator-task-additional-info.pdf)\n",
    "\n",
    "\n",
    "- Book *Nielsen, Michael A.; Chuang, Isaac (2000). Quantum Computation and Quantum Information, 10th Anniversary Edition, Section 8.2.3, Operator-Sum Representation*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Measurement\n",
    "\n",
    "State vector of the real quantum computer cannot be directly observed. All we can read out of qubit is a single classical bit. So best we can get as output from quantum computer is bitstring of size $n$ where $n$ is number of qubits. Reading the state from a qubit is called \"measurement\". When qubit is in superposition, measurment puts qubit in one of two classical states. If we read 1 from qubit, it will \"collapse\" to state |1> and will stay there - superposition is \"destroyed\", and any subsequent measurement will return the same result.\n",
    "\n",
    "Measurement is non-unitary operation on the state vector. But for simplicity, and because we can access state vector of our simulator, it is easier if we do it with a trick:\n",
    "\n",
    "We can simulate measurement by choosing element from the state vector with weighted random function. Elements with larger probability amplitude will be returned more often, and elements with smaller probability amplitude will be returned less often. Elements with zero probability will never be returned.\n",
    "\n",
    "For example, this state vector:\n",
    "\n",
    "```\n",
    "Index (dec)  Index (bin)  Amplitude   Probability\n",
    "=================================================\n",
    "0            00           0.70710678  0.5 (50%)\n",
    "1            01           0           0   ( 0%)\n",
    "2            10           0           0   ( 0%)\n",
    "3            11           0.70710678  0.5 (50%)\n",
    "```\n",
    "\n",
    "Our random function should return elements 00 and 11 equaly often and never return 01 and 10. If we execute it 1000 times (1000 shots) we should get something like this:\n",
    "\n",
    "```\n",
    "{\n",
    "  \"00\": 494,\n",
    "  \"11\": 506\n",
    "}\n",
    "```\n",
    "*(this is random, so it usually is not exactly 500/500 and that is completelly fine)*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Requirements\n",
    "\n",
    "It is expected that simulator can perform following:\n",
    "\n",
    "- initialize state\n",
    "\n",
    "- read program, and for each gate:\n",
    "    - calculate matrix operator\n",
    "    - apply operator (modify state)\n",
    "    \n",
    "- perform multi-shot measurement of all qubits using weighted random technique\n",
    "\n",
    "It is up to you how you will organize code, but this is our suggestion:\n",
    "\n",
    "### Input format (quantum program)\n",
    "\n",
    "It is enough if simulator takes program in following format:\n",
    "\n",
    "```\n",
    "[\n",
    "  { \"unitary\": [[0.70710678, 0.70710678], [0.70710678, -0.70710678]], \"target\": [0] }, \n",
    "  { \"unitary\": [ [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0] ], \"target\": [0, 1] }\n",
    "  ...\n",
    "]\n",
    "```\n",
    "\n",
    "Or, you can define unitaries in the code, and accept program with gate names instead:\n",
    "\n",
    "```\n",
    "[\n",
    "  { \"gate\": \"h\", \"target\": [0] }, \n",
    "  { \"gate\": \"cx\", \"target\": [0, 1] }\n",
    "  ...\n",
    "]\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulator program\n",
    "\n",
    "Not engraved into stone, but you can do something like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_ground_state(num_qubits):\n",
    "    # return vector of size 2**num_qubits with all zeroes except first element which is 1\n",
    "    return\n",
    "\n",
    "def get_operator(total_qubits, gate_unitary, target_qubits):\n",
    "    # return unitary operator of size 2**n x 2**n for given gate and target qubits\n",
    "    return\n",
    "\n",
    "def run_program(initial_state, program):\n",
    "    # read program, and for each gate:\n",
    "    #   - calculate matrix operator\n",
    "    #   - multiply state with operator\n",
    "    # return final state\n",
    "    return\n",
    "\n",
    "def measure_all(state_vector):\n",
    "    # choose element from state_vector using weighted random and return it's index\n",
    "    return\n",
    "\n",
    "def get_counts(state_vector, num_shots):\n",
    "    # simply execute measure_all in a loop num_shots times and\n",
    "    # return object with statistics in following form:\n",
    "    #   {\n",
    "    #      element_index: number_of_ocurrences,\n",
    "    #      element_index: number_of_ocurrences,\n",
    "    #      element_index: number_of_ocurrences,\n",
    "    #      ...\n",
    "    #   }\n",
    "    # (only for elements which occoured - returned from measure_all)\n",
    "    return\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example usage\n",
    "\n",
    "If your code is organized as we suggested, then usage will look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "None\n"
     ]
    }
   ],
   "source": [
    "# Define program:\n",
    "\n",
    "my_circuit = [\n",
    "{ \"gate\": \"h\", \"target\": [0] }, \n",
    "{ \"gate\": \"cx\", \"target\": [0, 1] }\n",
    "]\n",
    "\n",
    "\n",
    "# Create \"quantum computer\" with 2 qubits (this is actually just a vector :) )\n",
    "\n",
    "my_qpu = get_ground_state(2)\n",
    "\n",
    "\n",
    "# Run circuit\n",
    "\n",
    "final_state = run_program(my_qpu, my_circuit)\n",
    "\n",
    "\n",
    "# Read results\n",
    "\n",
    "counts = get_counts(final_state, 1000)\n",
    "\n",
    "print(counts)\n",
    "\n",
    "# Should print something like:\n",
    "# {\n",
    "#   \"00\": 502,\n",
    "#   \"11\": 498\n",
    "# }\n",
    "\n",
    "# Voila!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus requirements\n",
    "\n",
    "If you have implemented simulator as described above: congratulations!\n",
    "\n",
    "Now, if you wish you can continue improving it, first and useful thing to do would be to allow parametric gates:\n",
    "\n",
    "\n",
    "### Parametric gates\n",
    "\n",
    "For example, following gate:\n",
    "\n",
    "```\n",
    "[\n",
    "  [\"cos(theta/2)\", \"-exp(i * lambda) * sin(theta / 2)\"],\n",
    "  [\"exp(i * phi) * sin(theta / 2)\", \"exp(i * lambda + i * phi) * cos(theta / 2)\"]\n",
    "]\n",
    "```\n",
    "\n",
    "Contains strings with expressions, and expressions can contain variables (usually angles in radians).\n",
    "\n",
    "When your program gets gate like this, it should parse and evaluate expressions (with variables) and make a \"normal\" unitary matrix with values, which is then applied to the state vector.\n",
    "\n",
    "Example program with parametric gates:\n",
    "\n",
    "```\n",
    "[\n",
    "  { \"unitary\": [[\"cos(theta/2)\", \"-exp(i * lambda) * sin(theta / 2)\"], [\"exp(i * phi) * sin(theta / 2)\", \"exp(i * lambda + i * phi) * cos(theta / 2)\"]], \"params\": { \"theta\": 3.1415, \"phi\": 1.15708, \"lambda\": -3.1415 }, \"target\": [0] }\n",
    "  ...\n",
    "]\n",
    "```\n",
    "\n",
    "Or, if you have defined unitaries somewhere in the program, then:\n",
    "\n",
    "```\n",
    "[\n",
    "  { \"gate\": \"u3\", \"params\": { \"theta\": 3.1415, \"phi\": 1.5708, \"lambda\": -3.1415 }, \"target\": [0] }\n",
    "  ...\n",
    "]\n",
    "```\n",
    "\n",
    "Which your program translates to:\n",
    "\n",
    "```\n",
    "[\n",
    "  [ 0+0j,  1+0j],\n",
    "  [ 0+1j,  0+0j]\n",
    "]\n",
    "```\n",
    "\n",
    "\n",
    "### Allow running variational quantum algorithms\n",
    "\n",
    "With support for parametric gates, all you need to do is to allow global params - and your simulator will be able to run variational quantum algorithms!\n",
    "\n",
    "In that case, parametrized gates in your program will contain strings instead parameter values:\n",
    "\n",
    "```\n",
    "[\n",
    "  { \"gate\": \"u3\", \"params\": { \"theta\": \"global_1\", \"phi\": \"global_2\", \"lambda\": -3.1415 }, \"target\": [0] }\n",
    "  ...\n",
    "]\n",
    "```\n",
    "\n",
    "Notice `global_1` and `global_2` instead angle values, which you pass to `run_program` method:\n",
    "\n",
    "```\n",
    "final_state = run_program(my_qpu, my_circuit, { \"global_1\": 3.1415, \"global_2\": 1.5708 })\n",
    "```\n",
    "\n",
    "And that way you can use it in variational algorithms:\n",
    "\n",
    "```\n",
    "mu_qpu = [...]\n",
    "my_circuit = [...]\n",
    "\n",
    "def objective_function(params):\n",
    "    final_state = run_program(my_qpu, my_circuit, { \"global_1\": params[0], \"global_2\": params[1] })\n",
    "\n",
    "    counts = get_counts(final_state, 1000)\n",
    "    \n",
    "    # ...calculate cost here...\n",
    "    \n",
    "    return cost\n",
    "\n",
    "# initial values\n",
    "params = np.array([3.1415, 1.5708])\n",
    "\n",
    "# minimize\n",
    "minimum = minimize(objective_function, params, method=\"Powell\", tol=1e-6)\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional info/help\n",
    "\n",
    "Any questions? Ping us on Slack!\n",
    "\n",
    "To be blamed: Petar Korponaić,\n",
    "Quantastica\n",
    "\n",
    "May the force be with you!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}