Quantum State Tomography.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantum State Tomography\n",
    "\n",
    "Below is a simple code capable of performing quantum state tomography on a random quantum state. First, the quantum state is generated by sending the $|0\\rangle$ state through Z, Y, and Z rotations parametrized by random $\\theta_1, \\theta_2,$ and $ \\theta_3$, respectively.  Then, we perform measurements in the Pauli matrices' eigenvector basis so that we can compute the Stokes parameters. Remember that the Stokes parameters are defined as $S_0 = 1$ and $S_i = P_{|\\psi_i\\rangle} - P_{|\\psi_i^\\bot\\rangle}$ for $i = 1, 2, 3\\,$, where $|\\psi_i\\rangle$ is the eigenvector of eigenvalue 1 of the $i^{th}$ Pauli matrix and $P_{|\\psi_i\\rangle}$ the probability of measuring said state. We can finally compute the density matrix of the final state as $\\rho = \\frac{1}{2}\\sum_{i = 0}^3 S_i \\sigma_i\\,$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Real density matrix:\n",
      "[[0.96239578+0.j         0.07626519-0.17428071j]\n",
      " [0.07626519+0.17428071j 0.03760422+0.j        ]]\n",
      "q0: ─RZ─RY─RZ─\n",
      "\n",
      "Measurement of S_1:\n",
      "q0: ─RZ─RY─RZ─H─M─\n",
      "Counter({'0': 575906, '1': 424094})\n",
      "\n",
      "Measurement of S_2:\n",
      "q0: ─RZ─RY─RZ─S─H─M─\n",
      "Counter({'1': 673505, '0': 326495})\n",
      "\n",
      "Measurement of S_3:\n",
      "q0: ─RZ─RY─RZ─M─\n",
      "Counter({'0': 962833, '1': 37167})\n",
      "\n",
      "Tomographic approximation of the density matrix\n",
      "[[0.962833 (0.075906+0.173505j)]\n",
      " [(0.075906-0.173505j) 0.037167000000000006]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from qibo import Circuit, gates\n",
    "\n",
    "# Definition of parameters\n",
    "rng = np.random.default_rng()\n",
    "theta = rng.uniform(0, np.pi, 3)\n",
    "S = np.ones(4)\n",
    "Pauli = np.array([np.eye(2),[[0, 1], [1, 0]],[[0, -1j], [1j, 0]],[[1, 0], [0, -1]]], dtype = np.matrix)\n",
    "\n",
    "# We first run the circuit without performing any measurements so that we can compare the real density matrix with our approximation\n",
    "c0 = Circuit(1, density_matrix = True)\n",
    "c0.add(gates.RZ(0, theta[0]))\n",
    "c0.add(gates.RY(0, theta[1]))\n",
    "c0.add(gates.RZ(0, theta[2]))\n",
    "print(\"Real density matrix:\")\n",
    "print(c0().state())\n",
    "print(c0.draw())\n",
    "\n",
    "# We perform measurements in each of the Pauli matrices eigenvector basis\n",
    "for i in range(1, 4):\n",
    "    c = Circuit(1)\n",
    "    c.add(gates.RZ(0, theta[0]))\n",
    "    c.add(gates.RY(0, theta[1]))\n",
    "    c.add(gates.RZ(0, theta[2]))\n",
    "\n",
    "    # Change of basis\n",
    "    if(i == 1): c.add(gates.H(0))\n",
    "    elif(i == 2):\n",
    "        c.add(gates.S(0).dagger)\n",
    "        c.add(gates.H(0))\n",
    "\n",
    "    # Measurement\n",
    "    c.add(gates.M(0))\n",
    "    result = c(nshots = 1000000)\n",
    "\n",
    "    freq = result.frequencies()\n",
    "    S[i] = (freq['0'] - freq['1'])/(freq['0']+freq['1'])\n",
    "    print()\n",
    "    print(\"Measurement of S_\", i, \":\",  sep = '')\n",
    "    print(c.draw())\n",
    "    print(freq)\n",
    "\n",
    "# Compute of density matrix\n",
    "rho = 0.5*np.eye(2)\n",
    "for i in range(1, 4):\n",
    "    rho = rho + 0.5*S[i]*Pauli[i]\n",
    "\n",
    "print()\n",
    "print(\"Tomographic approximation of the density matrix\")\n",
    "print(rho)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For n-qubit states, tomography becomes much more complicated. Instead of 3 Stokes parameters, we need $4^n - 1$. First, a random state is generated by sending the $|0\\rangle$ state through a contraption parametrized by random $\\theta_i$. Then, we loop over all possible Stokes parameters; to compute $S_{i_1, \\dots, i_n}$, we will make projective measurements in the basis $\\{|\\psi_{i_1, j_1}\\rangle\\otimes\\cdots\\otimes|\\psi_{i_n, j_n}\\rangle\\}_{j_1, \\dots, j_n}$, where $\\{|\\psi_{i, j}\\rangle\\}_j$ are the eigenvectors of the $i^{th}$ Pauli matrix. Remember that $ S_{i_1, \\dots, i_n} = \\bigotimes_{k = 1}^n \\left(P_{|\\phi_{i_k}\\rangle} \\pm P_{|\\phi_{i_k}^\\bot \\rangle }\\right) = \\sum_{k_1, \\dots, k_n} \\pm P_{|\\phi_{i_{k_1}}\\rangle\\cdots|\\phi_{i_{k_n}}\\rangle}\\,$. On the tensor product, the + sign is used for 0 indeces and the - sign for 1, 2 or 3 indeces. Thus, on the summatory, the sign will be -1 to the power of the number of $\\bot$ terms that don't correspond to a 0 subindex. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Real density matrix:\n",
      "[[ 7.51925440e-02+6.93889390e-18j -3.28357040e-02-1.01944908e-01j\n",
      "  -1.52163401e-01-3.59224253e-02j  6.11042352e-02+2.47133755e-02j\n",
      "  -8.36030970e-03-1.92875672e-02j -1.23475702e-01-8.84679652e-02j\n",
      "   1.28935023e-03-3.54445498e-03j -4.91483267e-02+5.77456443e-02j]\n",
      " [-3.28357040e-02+1.01944908e-01j  1.52554323e-01+3.81639165e-17j\n",
      "   1.15151055e-01-1.90613923e-01j -6.01895230e-02+7.20521509e-02j\n",
      "   2.98006397e-02-2.91212586e-03j  1.73863916e-01-1.28773553e-01j\n",
      "   4.24247400e-03+3.29590344e-03j -5.68281673e-02-9.18514009e-02j]\n",
      " [-1.52163401e-01+3.59224253e-02j  1.15151055e-01+1.90613923e-01j\n",
      "   3.25087038e-01-2.08166817e-17j -1.35460141e-01-2.08193375e-02j\n",
      "   2.61327686e-02+3.50372403e-02j  2.92136234e-01+1.20039027e-01j\n",
      "  -9.15868695e-04+7.78870992e-03j  7.18716601e-02-1.40337063e-01j]\n",
      " [ 6.11042352e-02-2.47133755e-02j -6.01895230e-02-7.20521509e-02j\n",
      "  -1.35460141e-01+2.08193375e-02j  5.77780489e-02+7.80625564e-18j\n",
      "  -1.31331003e-02-1.29260232e-02j -1.29417490e-01-3.13098327e-02j\n",
      "  -1.17175012e-04-3.30412290e-03j -2.09606037e-02+6.30797181e-02j]\n",
      " [-8.36030970e-03+1.92875672e-02j  2.98006397e-02+2.91212586e-03j\n",
      "   2.61327686e-02-3.50372403e-02j -1.31331003e-02+1.29260232e-02j\n",
      "   5.87697934e-03+2.52619106e-17j  3.64215225e-02-2.18362915e-02j\n",
      "   7.65827878e-04+7.24821473e-04j -9.34770560e-03-1.90274600e-02j]\n",
      " [-1.23475702e-01+8.84679652e-02j  1.73863916e-01+1.28773553e-01j\n",
      "   2.92136234e-01-1.20039027e-01j -1.29417490e-01+3.13098327e-02j\n",
      "   3.64215225e-02+2.18362915e-02j  3.06849968e-01+3.12250226e-17j\n",
      "   2.05296014e-03+7.33743303e-03j  1.27670167e-02-1.52651258e-01j]\n",
      " [ 1.28935023e-03+3.54445498e-03j  4.24247400e-03-3.29590344e-03j\n",
      "  -9.15868695e-04-7.78870992e-03j -1.17175012e-04+3.30412290e-03j\n",
      "   7.65827878e-04-7.24821473e-04j  2.05296014e-03-7.33743303e-03j\n",
      "   1.89188773e-04+4.55364912e-18j -3.56479816e-03-1.32658993e-03j]\n",
      " [-4.91483267e-02-5.77456443e-02j -5.68281673e-02+9.18514009e-02j\n",
      "   7.18716601e-02+1.40337063e-01j -2.09606037e-02-6.30797181e-02j\n",
      "  -9.34770560e-03+1.90274600e-02j  1.27670167e-02+1.52651258e-01j\n",
      "  -3.56479816e-03+1.32658993e-03j  7.64719099e-02-7.80625564e-18j]]\n",
      "q0:     ─RY─RZ─o───────────────X─RY─RZ─o───────────────X─RY─RZ─o────────────── ...\n",
      "q1:     ───────X─RY─RZ─o───────|───────X─RY─RZ─o───────|───────X─RY─RZ─o────── ...\n",
      "q2:     ───────────────X─RY─RZ─o───────────────X─RY─RZ─o───────────────X─RY─RZ ...\n",
      "\n",
      "q0: ... ─X─RY─RZ─o───────────────X─RY─RZ─o───────────────X─RY─RZ─o──────────── ...\n",
      "q1: ... ─|───────X─RY─RZ─o───────|───────X─RY─RZ─o───────|───────X─RY─RZ─o──── ...\n",
      "q2: ... ─o───────────────X─RY─RZ─o───────────────X─RY─RZ─o───────────────X─RY─ ...\n",
      "\n",
      "q0: ... ───X─RY─RZ─o───────────────X─RY─RZ─o───────────────X─RY─RZ─o────────── ...\n",
      "q1: ... ───|───────X─RY─RZ─o───────|───────X─RY─RZ─o───────|───────X─RY─RZ─o── ...\n",
      "q2: ... RZ─o───────────────X─RY─RZ─o───────────────X─RY─RZ─o───────────────X─R ...\n",
      "\n",
      "q0: ... ─────X─RY─RZ─o───────────────X─\n",
      "q1: ... ─────|───────X─RY─RZ─o───────|─\n",
      "q2: ... Y─RZ─o───────────────X─RY─RZ─o─\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Qibo 0.2.9|ERROR|2024-07-25 16:28:55]: Symbolic Hamiltonian should be a ``sympy`` expression but is <class 'float'>.\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "Symbolic Hamiltonian should be a ``sympy`` expression but is <class 'float'>.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[3], line 80\u001b[0m\n\u001b[0;32m     77\u001b[0m     rho \u001b[38;5;241m=\u001b[39m rho \u001b[38;5;241m+\u001b[39m (\u001b[38;5;241m1\u001b[39m\u001b[38;5;241m/\u001b[39m(\u001b[38;5;241m2\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mn))\u001b[38;5;241m*\u001b[39mS\n\u001b[0;32m     79\u001b[0m \u001b[38;5;66;03m#Print results\u001b[39;00m\n\u001b[1;32m---> 80\u001b[0m rho \u001b[38;5;241m=\u001b[39m \u001b[43mSymbolicHamiltonian\u001b[49m\u001b[43m(\u001b[49m\u001b[43mrho\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241m.\u001b[39mmatrix\n\u001b[0;32m     81\u001b[0m \u001b[38;5;28mprint\u001b[39m(rho)\n\u001b[0;32m     82\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mFidelity = \u001b[39m\u001b[38;5;124m'\u001b[39m, (rho\u001b[38;5;129m@c0\u001b[39m()\u001b[38;5;241m.\u001b[39mstate())\u001b[38;5;241m.\u001b[39mtrace())\n",
      "File \u001b[1;32mc:\\Users\\alvar\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\qibo\\hamiltonians\\hamiltonians.py:326\u001b[0m, in \u001b[0;36mSymbolicHamiltonian.__init__\u001b[1;34m(self, form, nqubits, symbol_map, backend)\u001b[0m\n\u001b[0;32m    323\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mbackend \u001b[38;5;241m=\u001b[39m _check_backend(backend)\n\u001b[0;32m    325\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m form \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m--> 326\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mform\u001b[49m \u001b[38;5;241m=\u001b[39m form\n\u001b[0;32m    327\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m nqubits \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m    328\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mnqubits \u001b[38;5;241m=\u001b[39m nqubits\n",
      "File \u001b[1;32mc:\\Users\\alvar\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\qibo\\hamiltonians\\hamiltonians.py:357\u001b[0m, in \u001b[0;36mSymbolicHamiltonian.form\u001b[1;34m(self, form)\u001b[0m\n\u001b[0;32m    353\u001b[0m \u001b[38;5;129m@form\u001b[39m\u001b[38;5;241m.\u001b[39msetter\n\u001b[0;32m    354\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mform\u001b[39m(\u001b[38;5;28mself\u001b[39m, form):\n\u001b[0;32m    355\u001b[0m     \u001b[38;5;66;03m# Check that given form is a ``sympy`` expression\u001b[39;00m\n\u001b[0;32m    356\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(form, sympy\u001b[38;5;241m.\u001b[39mExpr):\n\u001b[1;32m--> 357\u001b[0m         \u001b[43mraise_error\u001b[49m\u001b[43m(\u001b[49m\n\u001b[0;32m    358\u001b[0m \u001b[43m            \u001b[49m\u001b[38;5;167;43;01mTypeError\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[0;32m    359\u001b[0m \u001b[43m            \u001b[49m\u001b[38;5;124;43mf\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mSymbolic Hamiltonian should be a ``sympy`` expression but is \u001b[39;49m\u001b[38;5;132;43;01m{\u001b[39;49;00m\u001b[38;5;28;43mtype\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mform\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;132;43;01m}\u001b[39;49;00m\u001b[38;5;124;43m.\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[0;32m    360\u001b[0m \u001b[43m        \u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    361\u001b[0m     \u001b[38;5;66;03m# Calculate number of qubits in the system described by the given\u001b[39;00m\n\u001b[0;32m    362\u001b[0m     \u001b[38;5;66;03m# Hamiltonian formula\u001b[39;00m\n\u001b[0;32m    363\u001b[0m     nqubits \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0\u001b[39m\n",
      "File \u001b[1;32mc:\\Users\\alvar\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\qibo\\config.py:46\u001b[0m, in \u001b[0;36mraise_error\u001b[1;34m(exception, message)\u001b[0m\n\u001b[0;32m     39\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Raise exception with logging error.\u001b[39;00m\n\u001b[0;32m     40\u001b[0m \n\u001b[0;32m     41\u001b[0m \u001b[38;5;124;03mArgs:\u001b[39;00m\n\u001b[0;32m     42\u001b[0m \u001b[38;5;124;03m    exception (Exception): python exception.\u001b[39;00m\n\u001b[0;32m     43\u001b[0m \u001b[38;5;124;03m    message (str): the error message.\u001b[39;00m\n\u001b[0;32m     44\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m     45\u001b[0m log\u001b[38;5;241m.\u001b[39merror(message)\n\u001b[1;32m---> 46\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m exception(message)\n",
      "\u001b[1;31mTypeError\u001b[0m: Symbolic Hamiltonian should be a ``sympy`` expression but is <class 'float'>."
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import qibo\n",
    "from qibo import Circuit, gates\n",
    "from qibo.hamiltonians import SymbolicHamiltonian\n",
    "from qibo.symbols import I, X, Y, Z\n",
    "\n",
    "#Returns the number such that its binary representation in n bits is the inverse of l's binary representation \n",
    "def invert_binary(l, n):\n",
    "    s = 0\n",
    "    for i in range(n):\n",
    "        s += (l%2)*(2**(n - 1 - i))\n",
    "        l //= 2\n",
    "    return s\n",
    "\n",
    "n = 3 #number of qubits\n",
    "m = 10 #number of rng layers\n",
    "shots = 1000000 #number of shots\n",
    "rng = np.random.default_rng()\n",
    "theta = rng.uniform(0, np.pi, 2*n*m)\n",
    "rho = 0\n",
    "c0 = Circuit(n, density_matrix = True)\n",
    "Pauli = np.array([I, X, Y, Z], dtype = qibo.symbols.Symbol)\n",
    "\n",
    "\n",
    "#Generates the state and prints its density matrix\n",
    "for i in range(m): \n",
    "    for j in range(n):\n",
    "        c0.add(gates.RY(j, theta[2*n*i+2*j]))\n",
    "        c0.add(gates.RZ(j, theta[2*n*i+2*j+1]))\n",
    "        c0.add(gates.CNOT(j, (j+1)%n))\n",
    "print(\"Real density matrix:\")\n",
    "print(c0().state())\n",
    "print(c0.draw())\n",
    "\n",
    "#Loop over all possible n-quatrit quaternary codes. \n",
    "for k in range(0, 4**n):\n",
    "    c = Circuit(n)\n",
    "\n",
    "    #Generates state\n",
    "    for i in range(m): \n",
    "        for j in range(n):\n",
    "            c.add(gates.RY(j, theta[2*n*i+2*j]))\n",
    "            c.add(gates.RZ(j, theta[2*n*i+2*j+1]))\n",
    "            c.add(gates.CNOT(j, (j+1)%n))\n",
    "\n",
    "    #Change of basis\n",
    "    for l in range(n):\n",
    "        if(((k//(4**l))%4) == 1) :\n",
    "            c.add(gates.H(l))\n",
    "        elif(((k//(4**l))%4) == 2) :\n",
    "            c.add(gates.S(l).dagger())\n",
    "            c.add(gates.H(l))\n",
    "    \n",
    "    #Measurement\n",
    "    c.add(gates.M(*range(n)))\n",
    "    result = c(nshots = shots)\n",
    "    freq = result.frequencies(binary = False)\n",
    "\n",
    "    #Computes the k-th Stokes parameter\n",
    "    S = 0\n",
    "    for l in range(2**n):\n",
    "        sign = 1\n",
    "        for j in range(n):\n",
    "            if((l//(2**j))%2 == 1 and (k//(4**j))%4 != 0):\n",
    "                sign *= -1\n",
    "\n",
    "        S += sign * freq[invert_binary(l, n)]\n",
    "    S /= shots\n",
    "\n",
    "    #Adds the k-th Pauli string term to rho\n",
    "    aux = k\n",
    "    P = Pauli[aux%4](0)\n",
    "    aux //= 4\n",
    "    for i in range(1, n):\n",
    "        P = P*Pauli[aux%4](i)\n",
    "        aux //= 4\n",
    "    rho = rho + (1/(2**n))*S\n",
    "\n",
    "#Print results\n",
    "rho = SymbolicHamiltonian(rho).matrix\n",
    "print(rho)\n",
    "print('Fidelity = ', (rho@c0().state()).trace())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code is much faster than the previous one. Instead of measuring in a new basis for each of the $4^n - 1$ Stokes parameters, it only uses $3^n$ settings, corresponding to all the possible Stoke parameters subindeces that don't include 0s. The Stokes parameters whose subindeces include 0s can be computed from the same set of probabilities as those that don't, but share the same non-zero terms. Thus, after each set of probabilities is obtained, the code loops over those subindeces that can be obtained by substituting some of the quatrits by 0. An array that keeps track of the n-quatrit codes already visited helps avoid duplicates. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Qibo 0.2.9|INFO|2024-08-02 11:24:13]: Using numpy backend on /CPU:0\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Real density matrix:\n",
      "[[ 8.17384857e-01-5.55111512e-17j -2.48602824e-03-5.67311974e-03j\n",
      "  -2.48169470e-03+1.55231663e-03j  5.61812989e-02-3.82182654e-01j]\n",
      " [-2.48602824e-03+5.67311974e-03j  4.69358145e-05-3.38813179e-21j\n",
      "  -3.22603851e-06-2.19456649e-05j  2.48169470e-03+1.55231663e-03j]\n",
      " [-2.48169470e-03-1.55231663e-03j -3.22603851e-06+2.19456649e-05j\n",
      "   1.04828165e-05-8.47032947e-22j -8.96387198e-04+1.05366462e-03j]\n",
      " [ 5.61812989e-02+3.82182654e-01j  2.48169470e-03-1.55231663e-03j\n",
      "  -8.96387198e-04-1.05366462e-03j  1.82557725e-01+1.38777878e-17j]]\n",
      "q0: ─RY─RZ─o───────X─\n",
      "q1: ───────X─RY─RZ─o─\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Qibo 0.2.9|WARNING|2024-08-02 11:24:14]: Calculating the dense form of a symbolic Hamiltonian. This operation is memory inefficient.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 8.17727e-01+0.00000e+00j -2.28050e-03-5.77250e-03j\n",
      "  -2.29400e-03+2.05650e-03j  5.57995e-02-3.82582e-01j]\n",
      " [-2.28050e-03+5.77250e-03j  1.99000e-04+0.00000e+00j\n",
      "   6.65000e-05+1.13000e-04j  2.84200e-03+1.30250e-03j]\n",
      " [-2.29400e-03-2.05650e-03j  6.65000e-05-1.13000e-04j\n",
      "  -1.42000e-04+0.00000e+00j -1.03950e-03+1.29850e-03j]\n",
      " [ 5.57995e-02+3.82582e-01j  2.84200e-03-1.30250e-03j\n",
      "  -1.03950e-03-1.29850e-03j  1.82216e-01+0.00000e+00j]]\n",
      "Fidelity =  (1.0004821481600734-4.301572119336239e-17j)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import qibo\n",
    "from qibo import Circuit, gates\n",
    "from qibo.hamiltonians import SymbolicHamiltonian\n",
    "from qibo.symbols import I, X, Y, Z\n",
    "\n",
    "#Returns the number such that its binary representation in n bits is the inverse of l's binary representation \n",
    "def invert_binary(l, n):\n",
    "    s = 0\n",
    "    for i in range(n):\n",
    "        s += (l%2)*(2**(n - 1 - i))\n",
    "        l //= 2\n",
    "    return s\n",
    "\n",
    "#Interpreting n as a code in tertiary and m as a code in binary, returns a code in quaternary\n",
    "#such that the i-th character is 0 if m's i-th character is, and n's i-th character + 1 otherwise\n",
    "def and_3_2(n, m):\n",
    "    if(m == 0): return 0\n",
    "    return 4*and_3_2(n//3, m//2) + (n%3 + 1)*(m%2)\n",
    "\n",
    "#Returns the number of ones in the binary representation of n\n",
    "def number_of_ones(n):\n",
    "    if(n < 2): return n\n",
    "    return number_of_ones(n//2) + n%2\n",
    "\n",
    "n = 2 #number of qubits\n",
    "m = 1 #number of rng layers\n",
    "shots = 1000000 #number of shots\n",
    "rng = np.random.default_rng()\n",
    "theta = rng.uniform(0, np.pi, 2*n*m)\n",
    "rho = 0\n",
    "c0 = Circuit(n, density_matrix = True)\n",
    "Pauli = np.array([I, X, Y, Z], dtype = qibo.symbols.Symbol)\n",
    "visited = np.zeros(4**n, dtype = bool)\n",
    "\n",
    "#Generates the state and prints its density matrix\n",
    "for i in range(m): \n",
    "    for j in range(n):\n",
    "        c0.add(gates.RY(j, theta[2*n*i+2*j]))\n",
    "        c0.add(gates.RZ(j, theta[2*n*i+2*j+1]))\n",
    "        c0.add(gates.CNOT(j, (j+1)%n))\n",
    "print(\"Real density matrix:\")\n",
    "print(c0().state())\n",
    "print(c0.draw())\n",
    "\n",
    "#Loop over all the posible n-trit tertiary codes.\n",
    "for k in range(0, 3**n):\n",
    "    c = Circuit(n)\n",
    "    \n",
    "    #Generates state\n",
    "    for i in range(m): \n",
    "        for j in range(n):\n",
    "            c.add(gates.RY(j, theta[2*n*i+2*j]))\n",
    "            c.add(gates.RZ(j, theta[2*n*i+2*j+1]))\n",
    "            c.add(gates.CNOT(j, (j+1)%n))\n",
    "\n",
    "    #Change of basis\n",
    "    for l in range(n):\n",
    "        if(((k//(3**l))%3) == 0) :\n",
    "            c.add(gates.H(l))\n",
    "        elif(((k//(3**l))%3) == 1) :\n",
    "            c.add(gates.S(l).dagger())\n",
    "            c.add(gates.H(l))\n",
    "    \n",
    "    #Measurementt \n",
    "    c.add(gates.M(*range(n)))\n",
    "    result = c(nshots = shots)\n",
    "    freq = result.frequencies(binary = False)\n",
    "\n",
    "    #Computes as many of the 4^n - 1 Stokes parameters as possible and sums the corresponding\n",
    "    #term to rho, but only if the parameter hasn't already been calculated\n",
    "    for i in range(2**n):\n",
    "        s = and_3_2(k, i)\n",
    "        if(not visited[s]):\n",
    "            visited[s] = True\n",
    "            S = 0\n",
    "            for l in range(2**n):\n",
    "                sign = (-1)**(number_of_ones(i&l)%2)\n",
    "                S += sign * freq[invert_binary(l, n)]\n",
    "            S /= shots\n",
    "            \n",
    "            P = Pauli[s%4](0)\n",
    "            s //= 4\n",
    "            for j in range(1, n):\n",
    "                P = P*Pauli[s%4](j)\n",
    "                s //= 4\n",
    "            rho = rho + (1/(2**n))*S*P\n",
    "\n",
    "\n",
    "#Print results\n",
    "rho = SymbolicHamiltonian(rho).matrix\n",
    "print(rho)\n",
    "print('Fidelity = ', (rho@c0().state()).trace())"
   ]
  }
 ],
 "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.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}