classical_shadows_numpy.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cTG5BdpXkr9t"
      },
      "source": [
        "#Classical Shadows using Numpy"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Introduction"
      ],
      "metadata": {
        "id": "3dKbn3pAhR4p"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "This tutorial explains how to compute the [Classical Shadows](https://arxiv.org/abs/2002.08953) of mixed quantum states and then use these shadows to estimate the expectation of observables, all implemented from scratch in Numpy."
      ],
      "metadata": {
        "id": "auxo1fjchQEP"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let $\\Phi: \\ \\mathcal{L}(\\mathbb{C} ^ {2 ^ n}) \\rightarrow \\mathcal{L}(\\mathbb{C} ^ {2 ^ n})$ be a quantum channel acting on $n$-qubit operators and let $f$ be a probability distribution defined over the set of all unitaries acting on $\\mathbb{C} ^ {2 ^ n}$. Define $T_{(f, \\Phi)}: \\  \\mathcal{L}(\\mathbb{C} ^ {2 ^ n}) \\rightarrow \\mathcal{L}(\\mathbb{C} ^ {2 ^ n})$ as $T_{(f, \\Phi)} (\\rho) =  \\int \\limits_{U} U  ^ {\\dagger} \\Phi (U \\rho U ^ {\\dagger}) U  \\ \\text{d}f(U) $. This process is called twirling the channel $\\Phi$ using the distriubtion $f$. Let $\\Phi$ be the measurement (in the computational basis) channel and let $f$ be a distribution such that $T_{(f, \\Phi)}$ is invertible. Consider a randomized protocol defined as follows: (i) sample an $n$-qubit unitary $U$ according to $f$, (ii) apply $U$ on $\\rho$, (iii) measure all qubits of the resulting state, giving the classical $n$-bit string $b$ as output, (iv) compute $\\hat{\\rho} = T_{(f, \\Phi)} ^ {-1} (U ^ {\\dagger} |b\\rangle \\langle b| U)$ classically."
      ],
      "metadata": {
        "id": "Pgm3txUchmFU"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "    \n",
        "$\\hat{\\rho}$ is called a *classical shadow* of $\\rho$. It is a random variable that reproduces the original state in expectation. For any observable $O$, the variance of the random variable $\\text{tr}(\\hat{\\rho}O)$ is upper bounded by a quantity called the *shadow norm* of $O$, denoted as $||O||_{\\text{sh}}$, which is independent of the state $\\rho$ but dependent on $f$, and hence we have the following theorem\n"
      ],
      "metadata": {
        "id": "P8GzsJZUiS1M"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "  **Theorem**: Let $\\rho \\in \\mathbb{C} ^ {2 ^ n}$ be a quantum state, $O_1, O_2, \\dots, O_m \\in \\mathbb{C} ^ {2 ^ n}$ be Hermitian matrices and $\\epsilon, \\delta \\in \\left[ 0, 1\\right]$ be accuracy parameters. By using $ T = \\mathcal{O} (\\frac{\\log{(M/ \\delta)}}{\\epsilon ^ 2} \\max \\limits_i ||O_i - \\frac{\\text{tr}(O_i)}{2 ^ n} \\mathbf{I}||_{\\text{sh}})$ classical shadows $\\left\\{ \\hat{\\rho}_1, \\hat{\\rho}_2, \\dots, \\hat{\\rho}_T \\right\\}$, we have $\\left| \\overline{o}_i - \\text{tr} (\\rho O_i) \\right| \\leq \\epsilon \\ \\forall \\ i$ with probability at least $1 - \\delta$. Here, $\\overline{o}_i$ is a median of means estimator of samples $\\text{tr}\\left(\\hat{\\rho}_1 O_i\\right), \\text{tr}\\left( \\hat{\\rho}_2 O_i\\right), \\dots, \\text{tr}\\left( \\hat{\\rho}_T O_i\\right)$."
      ],
      "metadata": {
        "id": "Hy9FNeWZi2dD"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "  This is remarkable because the number of shadows required to estimate the expectation of the observables depends only on (i) the shadow norm of the observable, (ii) the precision to which we want to estimate the expectation, and (iii) the error probability, and not on the state itself. Also, to estimate the expectation of $M$ observables, we only require $\\mathcal{O}(\\log M)$ classical shadows loaded. To make this procedure efficient, we should have $f$ defined over unitaries that enables computing $T_{(f, \\Phi)} ^ {-1} (U ^ {\\dagger} | b \\rangle \\langle b |U) $ classically feasible, with classical cost polynomial in $n$. Two such distributions $f$ have been proposed, each resulting in a different range of values for the shadow norm"
      ],
      "metadata": {
        "id": "NsUEU28ejDYF"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "\n",
        "*   *Clifford ensemble*: Here, $f$ is the discrete uniform distribution over all $n$-qubit Clifford unitaries. In this case, $ T_{(f, \\Phi)} ^ {-1}(\\rho) = \\left(2 ^ n + 1\\right) \\rho - \\mathbf I$. The stabilizer formalism can be used to work with classical shdows in this scenario. In this case, $|| O||_{\\text{sh}} = \\mathcal{O}(\\text{tr}(O ^ 2))$.\n",
        "*   *Pauli ensemble*: Here, $f$ is the discrete uniform distribution over the set $\\left\\{ \\mathbf I, H, HS ^ {\\dagger}\\right\\} ^ {\\otimes n}$. In this case, $ T_{(f, \\Phi)} ^ {-1}(\\rho) = \\bigotimes \\limits_{i = 1} ^ n \\left( 3 \\rho - \\mathbf I\\right)$ and $|| O||_{\\text{sh}} = \\mathcal{O}(4 ^ k)$, where $k$ is the total number of qubits which is non-trivially acted upon by $O$\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "    "
      ],
      "metadata": {
        "id": "hwroLJogjRy_"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "  For example, if the observable is a quantum state (fidelity with a pure state), one can choose to use the Clifford ensemble, and if the observable is a local observable (or can be written as a linear combination of a few local observables), then we can use the Pauli ensemble. Many relevant observables in nature can be written as linear combinations of a few observables that act non-trivially on only a few qubits. This makes classical shadows a very powerful ``storage facility'' for states as well. In many contexts, this is why classical shadows are considered an alternative to conventional quantum state tomography, which requires resources exponential in the number of qubits. The key fact here is that to estimate many properties of quantum states, up to acceptable precision and quality, one does not need to know the full classical description of the states. Similar complexity guarantees are extended to other important problems involving quantum states such as finding entanglement witnesses, direct fidelity estimation, estimating the output of non-linear functions, etc. "
      ],
      "metadata": {
        "id": "HBFerbH_lLqX"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "In this tutorial we shall load classical shadows of quantum states using the Clifford ensemble as well as the Pauli ensemble. Using the Pauli ensemble, we will estimate the expectation of local Pauli observables and with the Clifford ensemble, we will estimate the fidelity with a target pure stabilizer state."
      ],
      "metadata": {
        "id": "0li79vDak-GV"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Importing packages"
      ],
      "metadata": {
        "id": "Rbp-VOdQgrYn"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "RorGxLgvkuV7"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import math\n",
        "import random"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XyqId2UZqbn3"
      },
      "source": [
        "# Required functions"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Matrix initialisations and basic functions"
      ],
      "metadata": {
        "id": "i6J4DmiG1kvk"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Gates"
      ],
      "metadata": {
        "id": "gw267oVzmJSH"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "H_gate = (0.5 ** (0.5)) * np.matrix([[1,1], \n",
        "               [1,-1]])\n",
        "\n",
        "X_gate = np.matrix([[0 + 0j, 1 + 0j], \n",
        "               [1 + 0j, 0+ 0j]])\n",
        "\n",
        "Y_gate = np.matrix([[0 + 0j, -1j], \n",
        "               [1j, 0 + 0j]])\n",
        "\n",
        "Z_gate = np.matrix([[1 + 0j, 0 + 0j], \n",
        "               [0 + 0j, -1 + 0j]])\n",
        "\n",
        "Id = np.matrix([[1 + 0j, 0 + 0j], \n",
        "                [0 + 0j, 1 + 0j]])\n",
        "\n",
        "\n",
        "Y_basis = (1/(2 ** 0.5)) * np.matrix([[1 + 0j, -1j], \n",
        "                                      [1 + 0j, 1j]]) # This is the gate HS ^ \\dagger\n",
        "\n",
        "Y_basis_dag = np.conj(Y_basis).T\n",
        "\n",
        "\n",
        "CNOT = np.matrix([[1 + 0j, 0 + 0j, 0 + 0j, 0 + 0j], \n",
        "                  [0 + 0j, 1 + 0j, 0 + 0j, 0 + 0j],\n",
        "                  [0 + 0j, 0 + 0j, 0 + 0j, 1 + 0j],\n",
        "                  [0 + 0j, 0 + 0j, 1 + 0j, 0 + 0j]])\n",
        "\n",
        "S_gate = np.matrix([[1 + 0j, 0 + 0j], \n",
        "                    [0 + 0j, 0 + 1j]])\n",
        "\n",
        "S_dag_gate = np.matrix([[1 + 0j, 0 + 0j], \n",
        "                    [0 + 0j, 0 - 1j]])\n",
        "\n",
        "clifford_gate_dict = {'H': H_gate, 'S': S_gate, 'X':X_gate, 'Y': Y_gate, 'Z':Z_gate, \n",
        "                      'CX': CNOT, 'Y_basis': Y_basis, 'Y_basis_dag': Y_basis_dag, 'I': Id, 'S_dag': S_dag_gate}"
      ],
      "metadata": {
        "id": "L22M3hONmLgP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Tables"
      ],
      "metadata": {
        "id": "Yc6QZa4kmdgv"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "For each of the Clifford gates, we shall store their actions in tables so that this can be used to update the stabilizer tables of the stabilizer states we will be working with. More information regarding stabilizer tables and updating them by action of Clifford gates, as well as measuring them, can be seen in [Nielsen and Chuang](http://mmrc.amss.cas.cn/tlb/201702/W020170224608149940643.pdf) Stabilizer codes section. We will be following the techniques explained there."
      ],
      "metadata": {
        "id": "44Ca2S83mfza"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "As an example, consider the Hadamard gate. This gate maps $\\mathbf{I}$ to itself, $X$ to $Z$, $Y$ to $-Y$ and $Z$ to $X$. By labelling $\\mathbf{I}$ as index $0$, $X$ as index $1$, $Y$ as index $2$ and $Z$ as index $3$, we store the action of $H$ gate on Paulis using the martix \"H_table\" given in the next cell. The first row is $[0,1]$, meaning that $\\mathbf{I}$ is mapped to $\\mathbf{I}$ and a sign of $+1$ is multiplied. Similarly, all the actions discussed above can be seen in the table. For all the one qubit Clifford gates, we have stored their actions in similar tables.\n",
        "\n",
        "When it comes to the CNOT gate, we need to store its action on all possible pairs of Paulis. Hence we store its actions in a dictionary with $2$-tuples as keys and $3$-tuples as values. This dictionary can also be read in a similar way to the tables of the one-qubit gates. For example: $(1, 3): (2, 2, -1)$ means CNOT maps $X \\otimes Z$ to $-Y \\otimes Y$."
      ],
      "metadata": {
        "id": "ZRcmhQOOm-sU"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We also have a multiplcation table \"pauli_multi_table\" which can be read as follows: pauli_multi_table$[i][j]$ = $(ij, \\text{induced phase})$. Here, $i$ and $j$ refers to the indices that we have given to the Paulis, and the \"induced phase\" is the phase that could come as a result of multiplication between Paulis. For example, $XZ$ gives us $Y$ with a phase of $-i$. "
      ],
      "metadata": {
        "id": "Xcwd6pBOowgF"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "An $n$-qubit pure state stabilized by a generator set $\\langle P_1, P_2, \\dots, P_n \\rangle$ is stored using a $n \\times (n + 1)$ matrix, with its $i ^ {th}$ row denoting the Pauli string $P_i$, stored as a row of numbers each belonging in $\\left\\{ 0, 1, 2, 3\\right\\}$ (similar indexing as previously discussed). The last column stores the sign of the Pauli, that is, one of $\\left\\{ 1, -1, i, -i\\right\\}$. "
      ],
      "metadata": {
        "id": "kOJlHZZg4zUZ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "H_table = np.array([[0, 1], \n",
        "                    [3, 1], \n",
        "                    [2, -1],\n",
        "                    [1, 1]], dtype = 'complex128')\n",
        "\n",
        "S_table = np.array([[0, 1], \n",
        "                    [2, 1], \n",
        "                    [1, -1],\n",
        "                    [3, 1]], dtype = 'complex128')\n",
        "\n",
        "S_dag_table = np.array([[0, 1], \n",
        "                        [2, -1], \n",
        "                        [1, 1],\n",
        "                        [3, 1]], dtype = 'complex128')\n",
        "\n",
        "X_table = np.array([[0, 1], \n",
        "                    [1, 1], \n",
        "                    [2, -1],\n",
        "                    [3, -1]], dtype = 'complex128')\n",
        "\n",
        "Y_table = np.array([[0, 1], \n",
        "                    [1, -1], \n",
        "                    [2, 1],\n",
        "                    [3, -1]], dtype = 'complex128')\n",
        "\n",
        "Z_table = np.array([[0, 1], \n",
        "                    [1, -1], \n",
        "                    [2, -1],\n",
        "                    [3, 1]], dtype = 'complex128')\n",
        "\n",
        "\n",
        "pauli_multi_table = [[(0 + 0j, 1 + 0j), (1 + 0j, 1 + 0j), (2 + 0j, 1 + 0j), (3 + 0j, 1 + 0j)],\n",
        "                    [(1 + 0j, 1 + 0j), (0 + 0j, 1 + 0j), (3 + 0j, 0 + 1j), (2 + 0j, 0 - 1j)],\n",
        "                    [(2 + 0j, 1 + 0j), (3 + 0j, 0 - 1j), (0 + 0j, 1 + 0j), (1 + 0j, 0 + 1j)],\n",
        "                    [(3 + 0j, 1 + 0j), (2 + 0j, 0 + 1j), (1 + 0j, 0 - 1j), (0 + 0j, 1 + 0j)]]\n",
        "\n",
        "cnot_table = {(0 + 0j, 0 + 0j): (0 + 0j, 0 + 0j, 1 + 0j),\n",
        "              (0 + 0j, 1 + 0j): (0 + 0j, 1 + 0j, 1 + 0j),\n",
        "              (0 + 0j, 2 + 0j): (3 + 0j, 2 + 0j, 1 + 0j),\n",
        "              (0 + 0j, 3 + 0j): (3 + 0j, 3 + 0j, 1 + 0j),\n",
        "              (1 + 0j, 0 + 0j): (1 + 0j, 1 + 0j, 1 + 0j),\n",
        "              (1 + 0j, 1 + 0j): (1 + 0j, 0 + 0j, 1 + 0j),\n",
        "              (1 + 0j, 2 + 0j): (2 + 0j ,3 + 0j, 1 + 0j),\n",
        "              (1 + 0j, 3 + 0j): (2 + 0j ,2 + 0j, -1 + 0j),\n",
        "              (2 + 0j, 0 + 0j): (2 + 0j, 1 + 0j, 1 + 0j),\n",
        "              (2 + 0j, 1 + 0j): (2 + 0j, 0 + 0j, 1 + 0j),\n",
        "              (2 + 0j, 2 + 0j): (1 + 0j, 3 + 0j, -1 + 0j),\n",
        "              (2 + 0j, 3 + 0j): (1 + 0j, 2 + 0j, 1 + 0j),\n",
        "              (3 + 0j, 0 + 0j): (3 + 0j, 0 + 0j, 1 + 0j),\n",
        "              (3 + 0j, 1 + 0j): (3 + 0j, 1 + 0j, 1 + 0j),\n",
        "              (3 + 0j, 2 + 0j): (0 + 0j, 2 + 0j, 1 + 0j),\n",
        "              (3 + 0j, 3 + 0j): (0 + 0j, 3 + 0j, 1 + 0j)}\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "clifford_table_dict = {'X': X_table, 'Y': Y_table, 'Z': Z_table, 'S':S_table, 'S_dag':S_dag_table, 'H':H_table}\n",
        "\n",
        "\n",
        "# This function is used to convert an integer to a list of binary values.\n",
        "# 1. x: integer value\n",
        "# 2. BitNo: how many binary bits for x's expansion is required\n",
        "def binary(x, BitNo):\n",
        "  format(x, 'b').zfill(BitNo)\n",
        "  Binlist = [int(y) for y in list(format(x, 'b').zfill(BitNo))]\n",
        "  return Binlist\n"
      ],
      "metadata": {
        "id": "Jphb2gqa1m_u"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Tensor operations"
      ],
      "metadata": {
        "id": "t1tCcflg1xhQ"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We can use tensor operations to model the operation of quantum gates on $n$-qubit quantum states, treated as mode $n$-tensors. The details of this can be found in the matricsation and mode-$n$ multiplication section in this [paper](https://www.kolda.net/publication/TensorReview.pdf)."
      ],
      "metadata": {
        "id": "hFLm_TsR4RD3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# This function applies a gate on specific qubits of the state and returns the resulting state using matricisation of the tensor. \n",
        "\n",
        "# 1. state: Numpy array (vector)\n",
        "# 2. qubits_to_act: list of integers (0 is for the first qubit)\n",
        "# 3. gate: numpy matrix\n",
        "\n",
        "\n",
        "def across_qubits_gate(state, qubits_to_act, gate):\n",
        "  NQubits = int(np.log2(len(state)))\n",
        "\n",
        "  rest_of_qubits = []\n",
        "  for i in range(NQubits):\n",
        "    if i not in qubits_to_act:\n",
        "      rest_of_qubits.append(i)\n",
        "\n",
        "  state_tensor = np.reshape(state, tuple([2] * NQubits))\n",
        "  matricised_tensor = np.reshape(np.ndarray.transpose(state_tensor, qubits_to_act + rest_of_qubits), (2 ** len(qubits_to_act), 2 ** len(rest_of_qubits)))\n",
        "  answer = np.matmul(gate, matricised_tensor)\n",
        "  answer = np.asarray(answer)\n",
        "  p = qubits_to_act + rest_of_qubits\n",
        "\n",
        "  s = np.zeros(len(p), dtype = int)\n",
        "  \n",
        "  for i in list(range(len(p))):\n",
        "      s[p[i]] = i\n",
        "\n",
        "  Z = np.ndarray.transpose(np.reshape(answer, tuple([2] * NQubits)), s)\n",
        "  Z = np.reshape(Z, int(2 ** NQubits))\n",
        "\n",
        "  return(Z)\n"
      ],
      "metadata": {
        "id": "q16lyxiM10i3"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Pauli Ensemble"
      ],
      "metadata": {
        "id": "FDhIGIhbp190"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "In this section, we shall use the Pauli ensemble to generate Classical shadows and use it to estimate expectation with $1$ and $2$ local Pauli strings. The procedure is to simply measure each qubit in a uniformly randomly sampled Pauli basis, and then apply (classically) the function $T(X) = 3X-\\mathbf I$ on each qubit separately. So, that means that each shadow will be a fully separable matrix and hence, for an $n$-qubit state, a single shadow can be stored using a $4n$ dimensional vector.\n",
        "\n",
        "Let $\\rho \\in \\mathcal{L}(\\mathbb{C} ^ {2^n})$ be an $n$-qubit quantum state and let $\\left\\{ \\rho_i = \\bigotimes \\limits_{j = 1} ^ n \\rho_i ^ {(j)} \\Big| \\ i = 1, 2, \\dots, T\\right\\}$ be its $T$ classical shadows loaded using the Pauli ensemble. This set can be stored using a $T \\times 4n$ dimensional block matrix as follows:\n",
        "\n",
        "\\begin{align*}\n",
        "  \\begin{bmatrix}\n",
        "    \\text{vec}(\\rho_1 ^ {(1)}) ^ T & \\text{vec}(\\rho_1 ^ {(2)}) ^ T & \\dots & \\text{vec}(\\rho_1 ^ {(n)}) ^ T \\\\\n",
        "    \\text{vec}(\\rho_2 ^ {(1)}) ^ T & \\text{vec}(\\rho_2 ^ {(2)}) ^ T & \\dots & \\text{vec}(\\rho_2 ^ {(n)}) ^ T \\\\\n",
        "    \\vdots \\\\\n",
        "    \\text{vec}(\\rho_T ^ {(1)}) ^ T & \\text{vec}(\\rho_T ^ {(2)}) ^ T & \\dots & \\text{vec}(\\rho_T ^ {(n)}) ^ T\n",
        "  \\end{bmatrix}\n",
        "\\end{align*}\n",
        "\n",
        "To estimate the expectation of a Pauli string, with a single non-identity term $P$ featuring in the $i^{\\text{th}}$ qubit, we simply multiply the $i^{th}$ block column of this matrix with the vector $\\overline{\\text{vec}(P)}$. This results in a $T \\times 1$ column vector and a median of means estimator of the entries of the vector gives us the estimate of the expectaion value. The reason that we choose this way is that for any two matrices $X$ and $Y$, to compute the value $\\text{tr}(X^{\\dagger} Y)$, it is inefficient to compute the matrix $X ^ {\\dagger} Y$. $\\text{vec}(X) ^ {\\dagger} \\text{vec}(Y)$ gives the same answer with much less floating point computations. Even though we are computing the conjugate of this value, since the expectation is a real value, it will not change anything.\n",
        "\n",
        "If the Pauli string has multiple non-identity terms, simply compute the $T \\times 1$ column vector associated with each term, and then compute an element-wise product of the column vectors. Then, a median of means estimator results in the estimate of the expectation value."
      ],
      "metadata": {
        "id": "rRNKJ_mkp5Io"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# This function computes T Pauli Classical Shadows of a mixed state and returns it in the matrix form presented above.\n",
        "# 1. mixture: A list of 2-tuples of the form [(p_1, \\ket{\\psi_1}), (p_2, \\ket{\\psi_2}), ... ,(p_t, \\ket{\\psi_t})]. Here, p_i are probabilities and \\ket{\\psi_i} are pure states in the mixture.\n",
        "# 2. T: Number of Classical shadows to be computed\n",
        "\n",
        "def pcs(mixture, T):\n",
        "  num_qubits = int(np.log2(len(mixture[0][1])))\n",
        "  shadow_list = []\n",
        "  gate_set_list = ['H', 'Y_basis']\n",
        "  \n",
        "  for j in range(T):\n",
        "    temp_mixture = mixture.copy()\n",
        "    final_prob_array = np.zeros(2 ** num_qubits)\n",
        "    gate_list = []\n",
        "    for i in range(num_qubits):\n",
        "      temp = np.random.choice(3, 1, p=[1/3, 1/3, 1/3])[0]\n",
        "      if temp != 2:\n",
        "        gate = gate_set_list[temp]\n",
        "        gate_list.append((gate, i))\n",
        "        for k in range(len(temp_mixture)):\n",
        "          temp_state = across_qubits_gate(temp_mixture[k][1], [i], clifford_gate_dict[gate])\n",
        "          temp_mixture[k] = (temp_mixture[k][0], temp_state)\n",
        "      else:\n",
        "        gate_list.append(('I', i))\n",
        "    for pair in temp_mixture:\n",
        "      prob_array = np.array([np.abs(pair[1][k]) ** 2 for k in range(len(pair[1]))])\n",
        "      final_prob_array += pair[0] * prob_array\n",
        "    outcome = np.random.choice(2 ** num_qubits, 1, p = list(final_prob_array))[0]\n",
        "    outcome_bin = binary(outcome, num_qubits)\n",
        "    final_outcome = np.array([])\n",
        "    for i in range(len(gate_list)):\n",
        "      if gate_list[i][0] == 'Y_basis':\n",
        "        temp1 = clifford_gate_dict['Y_basis_dag']\n",
        "      else:\n",
        "        temp1 = clifford_gate_dict[gate_list[i][0]]\n",
        "      temp = temp1[:, outcome_bin[i]]\n",
        "      final_outcome = np.concatenate([final_outcome, np.array((3 * np.outer(temp, np.conj(temp))) - Id).flatten()], axis = 0)\n",
        "\n",
        "    shadow_list.append(final_outcome)\n",
        "  shadow_list = np.vstack(shadow_list)\n",
        "  return shadow_list\n",
        "\n",
        "# This function returns the expectaion value of a Pauli observable using the method explained above.\n",
        "# 1. shadow_matrix: Shadows as a matrix. Output of the previous function.\n",
        "# 2. observable_list: A list of single qubit Paulis and the indices on which they are acting. For example [(0, 'X'), (3, 'Z')] denotes the Pauli string X I I Z.\n",
        "# 3. N: Parameter in median of means estimation. The means will be taken of sublists of size N.\n",
        "def pcs_expectation(shadow_matrix, observable_list, N):\n",
        "  temp_mat = np.zeros((shadow_matrix.shape[0], len(observable_list)), dtype = 'complex128')\n",
        "  k = 0\n",
        "  for pair in observable_list:\n",
        "    temp = shadow_matrix[:, pair[0] * 4: (pair[0] + 1) * 4]\n",
        "    flattened_observable = np.conj(np.array(clifford_gate_dict[pair[1]].flatten())[0])\n",
        "    temp_mat[:, k] = temp @ flattened_observable\n",
        "    k += 1\n",
        "\n",
        "  exp_mat = np.prod(temp_mat, axis = 1)\n",
        "  chunks = [np.mean(exp_mat[x:x+N]) for x in range(0, len(exp_mat), N)]\n",
        "  return np.median(chunks)\n",
        "\n"
      ],
      "metadata": {
        "id": "c8ZXkGm3vj66"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now, let's prepare a mixture of $| + \\rangle ^ {\\otimes 4}$ with $0.8$ probability and a Haar random state with $0.2$ probability."
      ],
      "metadata": {
        "id": "65uXUwgfxoXU"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "n_qubits = 4 # Number of qubits the state should be prepared on \n",
        "\n",
        "global_mixture = [] # The mixture is stored in this list\n",
        "\n",
        "state = np.ones(2 ** n_qubits) # Plus state\n",
        "state = state/np.linalg.norm(state)\n",
        "\n",
        "global_mixture.append((0.8, state))\n",
        "\n",
        "state = np.random.uniform(-1, 1, 2 ** n_qubits) + (1j * np.random.uniform(-1, 1, 2 ** n_qubits)) # Haar random state\n",
        "state = state/np.linalg.norm(state)\n",
        "\n",
        "global_mixture.append((0.2, state))\n",
        "\n"
      ],
      "metadata": {
        "id": "V9ovH1ZzyBRp"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let the local Pauli string be $X I I I$. So the acual expectation value of our mixture with this Pauli string is"
      ],
      "metadata": {
        "id": "Jgamny9lykua"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "exp_val = 0\n",
        "for pair in global_mixture:\n",
        "  temp = across_qubits_gate(pair[1], [0], X_gate)\n",
        "  exp_val += pair[0] * np.dot(np.conj(pair[1]), temp)\n",
        "\n",
        "print(\"Actual expectation value: {}\".format(exp_val))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "b1_36Dgnygan",
        "outputId": "d4c3910b-4bd1-4b69-9f8f-193886c76731"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Actual expectation value: (0.8792299857815542+5.551115123125783e-18j)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now, let's estimate this value using $4000$ Classical shadows."
      ],
      "metadata": {
        "id": "zjf5Kux8y65m"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "shadows = pcs(global_mixture, 4000)"
      ],
      "metadata": {
        "id": "3dfGVYK8y20Q"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "pcs_expectation(shadows, [(0, 'X')], N = 50, K = 80)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "9aGbyTgWzCB1",
        "outputId": "aeaf327d-ac72-4fc7-c9fd-1073bca2643d"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(0.8400000000000003+0j)"
            ]
          },
          "metadata": {},
          "execution_count": 71
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "As we can see, the estimated value is close to the true value."
      ],
      "metadata": {
        "id": "pA8NS6Nv0XtC"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now we shall choose a different observable, the Pauli string $X I Y I$. Then, the actual expectation value is given as "
      ],
      "metadata": {
        "id": "s76teZE-EU2I"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "exp_val = 0\n",
        "for pair in global_mixture:\n",
        "  temp = across_qubits_gate(pair[1], [0], X_gate)\n",
        "  temp = across_qubits_gate(temp, [2], Y_gate)\n",
        "  exp_val += pair[0] * np.dot(np.conj(pair[1]), temp)\n",
        "\n",
        "print(\"Actual expectation value: {}\".format(exp_val))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "shBlc1RPEjEk",
        "outputId": "e3c568e0-282b-47ec-e4c5-e3d9658bdfa6"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Actual expectation value: (0.013211411954024485-3.4694469519536144e-19j)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Using Classical shadows, we get the estimate"
      ],
      "metadata": {
        "id": "L8_wFKgzEtMN"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "pcs_expectation(shadows, [(0, 'X'), (2, 'Y')], N = 50, K = 80)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "kTk20AxQFWFh",
        "outputId": "46bc59fc-9f64-4548-9ba3-7823b3626a4d"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0j"
            ]
          },
          "metadata": {},
          "execution_count": 73
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "This is also pretty close to the actual answer."
      ],
      "metadata": {
        "id": "KSUV-ScrFcdr"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Clifford Ensemble"
      ],
      "metadata": {
        "id": "9yU5tXeM0oVj"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "In this section, we shall use the Classical shadows generated using the Clifford ensemble to estimate the fidelity between a mixed state $\\rho$ and a pure stabilizer state $| \\psi \\rangle$, that is, the value $\\langle \\psi | \\rho | \\psi \\rangle $. Let $\\left\\{ \\rho_i = (2 ^ n + 1) U_i ^ {\\dagger} | b_i \\rangle \\langle b_i | U_i - \\mathbf{I}\\ | \\ i = 1,2, \\dots, T \\right\\}$ be shadows of $\\rho$ generated using the Clifford ensemble. Then, each $\\rho_i$ is stored as a stabilizer table. The fidelity between $|\\psi\\rangle$ is then estimated using the median of means estimator of the values $\\left\\{ (2 ^ n + 1)\\left| \\langle \\psi | U_i ^ {\\dagger} | b_i \\rangle \\right| ^ 2 - 1 \\ \\Big| \\ i = 1, 2, \\dots, T \\right\\}$. So, we should be able to compute these fidelity values efficiently.\n",
        "\n",
        "On top of $| \\psi \\rangle$ being a stabilizer state, we also assume that we know the Clifford circuit that can be used to prepare $| \\psi \\rangle$. This will allow us to compute the fidelity values efficiently. To compute the pure state fidelity between two stabilizer states $|\\alpha \\rangle = U |0 \\rangle$ and $|\\beta \\rangle = V |0 \\rangle$, first we compute the stabilizer table of the state given as $V ^ {\\dagger} U | 0 \\rangle$. Then, by using property of stabilizer formalism, we simply compute the probability that all qubits will be measured to be in $0$. This gives us the required fidelity. \n",
        "\n",
        "Such a measurement will require us to check whether a given Pauli string is an element of the subgroup generated by the Paulis in the stabilizer table. This is implemented using the Row-Reduced Echelon Form technique discussed in this [work](https://studentutsedu-my.sharepoint.com/personal/afrad_m_basheer_student_uts_edu_au/Documents/Microsoft%20Teams%20Chat%20Files/Dimitris-Kakofengitis-Dissertation.pdf). To generate a circuit uniformly sampled from the Clifford group, we shall use the technique given in this [work](https://arxiv.org/abs/2008.06011)."
      ],
      "metadata": {
        "id": "QRtCqTP50tq6"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# This function generates a circuit which is uniformly sampled from the Clifford group.\n",
        "# 1. n_qubits: Number of qubits the circuit should act on\n",
        "\n",
        "def generate_random_clifford_circuit(n_qubits):\n",
        "\n",
        "  gate_list = []\n",
        "\n",
        "  for i in range(n_qubits):\n",
        "    m = n_qubits - i\n",
        "    while 1 == 1:\n",
        "      \n",
        "      a_x = np.random.choice(2, m, p = [0.5, 0.5])\n",
        "      b_x = np.random.choice(2, m, p = [0.5, 0.5])\n",
        "\n",
        "      a_z = np.random.choice(2, m, p = [0.5, 0.5])\n",
        "      b_z = np.random.choice(2, m, p = [0.5, 0.5])\n",
        "      count = 0\n",
        "      for j in range(m):\n",
        "        if (a_x[j], a_z[j]) != (b_x[j], b_z[j]):\n",
        "          if (a_x[j], a_z[j]) != (0, 0) and (b_x[j], b_z[j]) != (0, 0):\n",
        "            count += 1\n",
        "      if (count % 2) == 1:\n",
        "        break\n",
        "    x = np.vstack([a_x, b_x]) \n",
        "    z = np.vstack([a_z, b_z])\n",
        "    \n",
        "    for j in range(m):\n",
        "      if z[0, j] == 1:\n",
        "        if x[0, j] == 0:\n",
        "          gate_list.append(('H', [j + i]))\n",
        "          temp = x[:, j].copy()\n",
        "          x[:, j] = z[:, j]\n",
        "          z[:, j] = temp\n",
        "        else:\n",
        "          gate_list.append(('S', [j + i]))\n",
        "          z[:, j] += x[:, j]\n",
        "          z = z % 2\n",
        "\n",
        "    j_list = []\n",
        "    for j in range(m):\n",
        "      if x[0, j] == 1:\n",
        "        j_list.append(j)\n",
        "\n",
        "    flag = 0\n",
        "    while flag == 0:\n",
        "      m_1 = len(j_list)\n",
        "      new_j_list = []\n",
        "      for j in range(0, m_1, 2):\n",
        "        if j != m_1 - 1:\n",
        "          gate_list.append(('CX', [j_list[j] + i, j_list[j+1] + i]))\n",
        "          z[:, j_list[j]] += z[:, j_list[j + 1]]\n",
        "          z = z % 2\n",
        "\n",
        "          x[:, j_list[j + 1]] += x[:, j_list[j]]\n",
        "          x = x % 2\n",
        "        new_j_list.append(j_list[j])\n",
        "\n",
        "      j_list = new_j_list\n",
        "      if len(j_list) == 1:\n",
        "        flag += 1\n",
        "\n",
        "    if j_list[0] != 0:\n",
        "      gate_list.append(('CX', [0 + i, j_list[0] + i]))\n",
        "      z[:, 0] += z[:, j_list[0]]\n",
        "      z = z % 2\n",
        "\n",
        "      x[:, j_list[0]] += x[:, 0]\n",
        "      x = x % 2\n",
        "\n",
        "      gate_list.append(('CX', [j_list[0] + i, 0 + i]))\n",
        "      z[:, j_list[0]] += z[:, 0]\n",
        "      z = z % 2\n",
        "\n",
        "      x[:, 0] += x[:, j_list[0]]\n",
        "      x = x % 2\n",
        "\n",
        "      gate_list.append(('CX', [0 + i, j_list[0] + i]))\n",
        "      z[:, 0] += z[:, j_list[0]]\n",
        "      z = z % 2\n",
        "\n",
        "      x[:, j_list[0]] += x[:, 0]\n",
        "      x = x % 2\n",
        "    flag1 = 0\n",
        "    for j in range(1, m):\n",
        "      if z[1, j] == 1:\n",
        "        flag1 += 1\n",
        "\n",
        "    flag1 += ((1 + z[1, 0]) % 2)\n",
        "    flag1 += sum(x[1,:])\n",
        "\n",
        "    if flag1 != 0:\n",
        "\n",
        "      gate_list.append(('H', [0 + i]))\n",
        "      temp = x[:, 0].copy()\n",
        "      x[:, 0] = z[:, 0]\n",
        "      z[:, 0] = temp\n",
        "\n",
        "      for j in range(m):\n",
        "        if z[1, j] == 1:\n",
        "          if x[1, j] == 0:\n",
        "            gate_list.append(('H', [j + i]))\n",
        "            #x[:, j], z[:, j] = z[:, j], x[:, j]\n",
        "            temp = x[:, j].copy()\n",
        "            x[:, j] = z[:, j]\n",
        "            z[:, j] = temp\n",
        "          else:\n",
        "            gate_list.append(('S', [j + i]))\n",
        "            z[:, j] += x[:, j]\n",
        "            z = z % 2\n",
        "\n",
        "      j_list = []\n",
        "      for j in range(m):\n",
        "        if x[1, j] == 1:\n",
        "          j_list.append(j)\n",
        "      flag = 0\n",
        "      m_1 = 0\n",
        "\n",
        "      while flag == 0:\n",
        "        m_1 = len(j_list)\n",
        "        new_j_list = []\n",
        "        for j in range(0, m_1, 2):\n",
        "          if j != m_1 - 1:\n",
        "            gate_list.append(('CX', [j_list[j] + i, j_list[j+1] + i]))\n",
        "            z[:, j_list[j]] += z[:, j_list[j + 1]]\n",
        "            z = z % 2\n",
        "\n",
        "            x[:, j_list[j + 1]] += x[:, j_list[j]]\n",
        "            x = x % 2\n",
        "          new_j_list.append(j_list[j])\n",
        "\n",
        "        j_list = new_j_list\n",
        "\n",
        "        if len(j_list) == 1:\n",
        "          flag += 1\n",
        "\n",
        "      gate_list.append(('H', [0 + i]))\n",
        "      temp = x[:, 0].copy()\n",
        "      x[:, 0] = z[:, 0]\n",
        "      z[:, 0] = temp\n",
        "    a_s = np.random.choice(2, 1, p = [0.5, 0.5])[0]\n",
        "    b_s = np.random.choice(2, 1, p = [0.5, 0.5])[0]\n",
        "\n",
        "    if a_s == 0:\n",
        "      if b_s == 1:\n",
        "        gate_list.append(('X', [0 + i]))\n",
        "\n",
        "    else:\n",
        "      if b_s == 0:\n",
        "        gate_list.append(('Z', [0 + i]))\n",
        "      else:\n",
        "        gate_list.append(('Y', [0 + i]))\n",
        "\n",
        "  return gate_list\n",
        "\n",
        "# This function converts a binary string (standard basis vector) into its corresponding stabilizer table.\n",
        "# 1. bin_list: List of binary values\n",
        "def convert_bit_string(bin_list):\n",
        "  n_qubits = len(bin_list)\n",
        "  answer = np.zeros((n_qubits, n_qubits + 1), dtype = 'complex128')\n",
        "  for i in range(n_qubits):\n",
        "    if bin_list[i] == 0:\n",
        "      answer[i][i] = 3 + 0j\n",
        "      answer[i][n_qubits] = 1 + 0j\n",
        "    else:\n",
        "      answer[i][i] = 3 + 0j\n",
        "      answer[i][n_qubits] = -1 + 0j\n",
        "  return answer\n",
        "\n",
        "# This function applies a Clifford gate on a stabilizer state given as a stabilizer table, and returns the updated table.\n",
        "# 1. stab: Stabilizer table of the state. n x n+1 numpy matrix\n",
        "# 2. gate: The Clifford gate to act given as a letter. For example 'H' for Hadamard\n",
        "# 3. qubits_to_act: The qubit on which the gate should be applied. Note that the first qubit is 0 and the qubits should always be given in a list. \n",
        "def apply_gate_stabilizer(stab, gate, qubits_to_act):\n",
        "  n_qubits = stab.shape[0]\n",
        "  if len(qubits_to_act) == 1:\n",
        "    for i in range(n_qubits):\n",
        "      stab[i][qubits_to_act[0]] = clifford_table_dict[gate][int(np.real(stab[i][qubits_to_act[0]]))][0]\n",
        "      \n",
        "      stab[i][n_qubits] = stab[i][n_qubits] * clifford_table_dict[gate][int(np.real(stab[i][qubits_to_act[0]]))][1]\n",
        "  else:\n",
        "    for i in range(n_qubits):\n",
        "      temp = cnot_table[(stab[i][qubits_to_act[0]], stab[i][qubits_to_act[1]])]\n",
        "      stab[i][qubits_to_act[0]] = temp[0]\n",
        "      stab[i][qubits_to_act[1]] = temp[1]\n",
        "      stab[i][n_qubits] = stab[i][n_qubits] * temp[2]\n",
        "\n",
        "  return stab\n",
        "\n",
        "# This function checks if a given Pauli string commutes with element in the generator set (rows of the stabilizer table). \n",
        "# It returns (does it commute with everyone, updated stabilizer table with only one generator anti commuting with the pauli string, the index of the generator that it anti-commutes with)\n",
        "# 1. stab: Stabilizer table of th state. n x n+1 numpy matrix\n",
        "# 2. pauli_string: The Pauli string given as a list of previously mentioned indices for Paulis\n",
        "def commutator_check(stab, pauli_string):\n",
        "  n_qubits = len(pauli_string)\n",
        "  for i in range(n_qubits):\n",
        "    sign = 1\n",
        "    for j in range(n_qubits):\n",
        "      if stab[i][j] != pauli_string[j] and (stab[i][j] * pauli_string[j]) != 0:\n",
        "        sign *= -1\n",
        "\n",
        "    if sign != 1:\n",
        "      \n",
        "      for k in range(i + 1, n_qubits):\n",
        "        sign = 1\n",
        "        for l in range(n_qubits):\n",
        "          if stab[k][l] != pauli_string[l] and (stab[k][l] * pauli_string[l]) != 0:\n",
        "            sign *= -1\n",
        "\n",
        "        if sign != 1:\n",
        "          stab = row_addition(stab, [i, k])\n",
        "      return 0, i, stab\n",
        "\n",
        "  return 1, 1000, stab\n",
        "\n",
        "# This function can be used to multiply a specific generator on another generator. It returns the updated stabilizer table.\n",
        "# 1. stab: Stabilizer table of the state. n x n+1 numpy matrix\n",
        "# 2. rows: If P_0 is row 0 and P_2 is row 2, to make row 2 P_1P_2, give this argument [0, 2]\n",
        "def row_addition(stab, rows):\n",
        "  n_qubits = stab.shape[1] - 1\n",
        "  for i in range(n_qubits):\n",
        "    stab[rows[1]][i] = pauli_multi_table[int(np.real(stab[rows[0]][i]))][int(np.real(stab[rows[1]][i]))][0]\n",
        "    stab[rows[1]][n_qubits] *= pauli_multi_table[int(np.real(stab[rows[0]][i]))][int(np.real(stab[rows[1]][i]))][1]\n",
        "\n",
        "  stab[rows[1]][n_qubits] *= stab[rows[0]][n_qubits]\n",
        "  return stab\n",
        "\n",
        "# The rref subroutine. The returns the rref form of the stabilizer table.\n",
        "# 1. stab: Stabilizer table of the state. n x n+1 numpy matrix\n",
        "def rref(stab):\n",
        "  n_qubits = stab.shape[1] - 1\n",
        "  number_of_generators = stab.shape[0]\n",
        "  current_row = 0\n",
        "  for column in range(n_qubits):\n",
        "    uni = np.unique(stab[current_row:, column], return_counts=True)[0]\n",
        "    uni_non_zero = [x for x in uni if x != 0]\n",
        "    if len(uni_non_zero) == 1:\n",
        "      ks = []\n",
        "      for row in range(current_row, number_of_generators):\n",
        "        if stab[row][column] != 0:\n",
        "          ks.append(row)\n",
        "\n",
        "      stab[[ks[0], current_row]] = stab[[current_row, ks[0]]]\n",
        "      for row in range(current_row + 1, number_of_generators):\n",
        "        if stab[row][column] != 0:\n",
        "          stab = row_addition(stab, [current_row, row])\n",
        "      current_row += 1\n",
        "      \n",
        "    elif len(uni_non_zero) >= 2:\n",
        "      ks = []\n",
        "      for row in range(current_row, number_of_generators):\n",
        "        \n",
        "        if stab[row][column] != 0:\n",
        "          if len(ks) < 1:\n",
        "            ks.append((stab[row][column], row))\n",
        "\n",
        "          elif len(ks) < 2:\n",
        "            if ks[0][0] != stab[row][column]:\n",
        "              ks.append((stab[row][column], row))\n",
        "\n",
        "      stab[[ks[0][1], current_row]] = stab[[current_row, ks[0][1]]]\n",
        "      stab[[ks[1][1], current_row + 1]] = stab[[current_row + 1, ks[1][1]]]\n",
        "\n",
        "      if current_row + 2 <= number_of_generators-1:\n",
        "        for row in range(current_row + 2, number_of_generators):\n",
        "          if stab[row][column] == ks[0][0]:\n",
        "            stab = row_addition(stab, [current_row, row])\n",
        "          elif stab[row][column] == ks[1][0]:\n",
        "            stab = row_addition(stab, [current_row + 1, row])\n",
        "          else:\n",
        "            if stab[row][column] != 0:\n",
        "              stab = row_addition(stab, [current_row, row])\n",
        "              stab = row_addition(stab, [current_row + 1, row])\n",
        "\n",
        "      current_row += 2\n",
        "  return stab\n",
        "\n",
        "\n",
        "# This function computes T classical shadows of a state using the Clifford ensemble. It returns a list of T stabilizer tables.\n",
        "# 1. mixture: A list of 2-tuples of the form [(p_1, \\ket{\\psi_1}), (p_2, \\ket{\\psi_2}), ... ,(p_t, \\ket{\\psi_t})]. Here, p_i are probabilities and \\ket{\\psi_i} are pure states in the mixture.\n",
        "# 2. T: number of shadows to load\n",
        "def ccs(mixture, T):\n",
        "  num_qubits = int(np.log2(len(mixture[0][1])))\n",
        "  dim = 2 ** num_qubits\n",
        "  stab_list = []\n",
        "  for j in range(T):\n",
        "    outer_prob_array = np.zeros(dim)\n",
        "    randomly_generated_clifford_circuit = generate_random_clifford_circuit(num_qubits)\n",
        "    m = len(randomly_generated_clifford_circuit)\n",
        "    for pair in mixture:\n",
        "      Ustate = np.copy(pair[1])\n",
        "      for i in range(m):\n",
        "        Ustate = across_qubits_gate(Ustate, randomly_generated_clifford_circuit[i][1], clifford_gate_dict[randomly_generated_clifford_circuit[i][0]])\n",
        "      outer_prob_array += pair[0] * np.array([np.abs(Ustate[k]) ** 2 for k in range(dim)])\n",
        "    outcome = np.random.choice(dim, 1, p = list(outer_prob_array))[0]\n",
        "    outcome_bin = binary(outcome, num_qubits)\n",
        "    stab = convert_bit_string(outcome_bin)\n",
        "    randomly_generated_clifford_circuit.reverse()\n",
        "    for i in range(m):\n",
        "      if randomly_generated_clifford_circuit[i][0] == 'S':\n",
        "        stab = apply_gate_stabilizer(stab, 'S_dag', randomly_generated_clifford_circuit[i][1])\n",
        "      else:\n",
        "        stab = apply_gate_stabilizer(stab, randomly_generated_clifford_circuit[i][0], randomly_generated_clifford_circuit[i][1])\n",
        "    stab_list.append(stab)\n",
        "  return stab_list\n",
        "\n",
        "# This function returns the probability that the given stabilizer state will result in all zeros after measrement of all qubits.\n",
        "# 1. stab: Stabilizer table of the state. n x n+1 numpy matrix\n",
        "def measure_clifford_prob_0(stab):\n",
        "  prob_list = []\n",
        "  for qubit in range(stab.shape[0]):\n",
        "    pauli_string = list(3 * np.eye(stab.shape[0], dtype = 'complex128')[qubit])\n",
        "    temp, index, stab = commutator_check(stab, pauli_string)\n",
        "    \n",
        "    if temp == 1:\n",
        "      \n",
        "      updated_pauli_string = np.concatenate([pauli_string, np.array([1 + 0j])], axis = 0)\n",
        "      stab1 = np.concatenate([stab, np.array([updated_pauli_string])], axis = 0)\n",
        "      \n",
        "      stab1 = rref(stab1)\n",
        "      flag = 0\n",
        "      for i in range(stab1.shape[0]):\n",
        "        if np.linalg.norm(stab1[i,:-1]) == 0 and stab1[i,-1] == 1:\n",
        "          flag = 1\n",
        "          \n",
        "      if flag == 1:\n",
        "        prob_list.append(1)\n",
        "      else:\n",
        "        prob_list.append(0)\n",
        "\n",
        "    else:\n",
        "      prob_list.append(0.5)\n",
        "      stab[index] = np.concatenate([np.array(pauli_string), np.array([1 + 0j])], axis = 0)\n",
        "\n",
        "  return np.prod(np.array(prob_list))\n",
        "\n",
        "# This function computes the fidelity between two pure stabilizer states, one given as a stabilizer table, and the other given as a circuit that prepares the state.\n",
        "# 1. stab: Stabilizer table of the state. n x n+1 numpy matrix\n",
        "# 2. circuit: The circuit that prepares the other state. For example: [('H', [0]), ('CNOT', [1,5]), ('S', [3])]\n",
        "def compute_fidelity_with_circuit(stab, circuit):\n",
        "  for i in range(len(circuit)-1, -1, -1):\n",
        "    if circuit[i][0] == 'S':\n",
        "      stab = apply_gate_stabilizer(stab, 'S_dag', circuit[i][1])\n",
        "    else:\n",
        "      stab = apply_gate_stabilizer(stab, circuit[i][0], circuit[i][1])\n",
        "\n",
        "  return measure_clifford_prob_0(stab)\n",
        "\n",
        "# This function estimates the fidelity of a state given as its shadows with another stabilizer state whose preparation circuit is given\n",
        "# 1. shadows: A list of stabilizer tables. Output of the ccs function.\n",
        "# 2. circuit: The circuit that prepares the other state. For example: [('H', [0]), ('CNOT', [1,5]), ('S', [3])]\n",
        "# 3. N: Parameter in median of means estimation. The means will be taken of sublists of size N.\n",
        "def compute_fidelity_with_shadows(shadows, circuit, N):\n",
        "  exp_list = []\n",
        "  num_qubits = shadows[0].shape[0]\n",
        "  for stab in shadows:\n",
        "    temp = compute_fidelity_with_circuit(stab, circuit)\n",
        "    exp_list.append((((2 ** num_qubits) + 1) * temp) - 1)\n",
        "  chunks = [np.mean(exp_list[x:x+N]) for x in range(0, len(exp_list), N)]\n",
        "  return np.median(chunks)"
      ],
      "metadata": {
        "id": "GS0vvFOq2JKj"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now, let's prepare a mixture of a randomly generated pure stabilizer state $| \\psi \\rangle$ with probability $0.8$ and a Haar random state with probability $0.2$. We shall try to estimate the fidelity between this mixture and $| \\psi \\rangle$."
      ],
      "metadata": {
        "id": "w3_7KgPd-ff3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "n_qubtis = 4\n",
        "cliff_state = np.array([1] + [0] * ((2 ** n_qubtis)-1)) # initial state ket 0\n",
        "circ = generate_random_clifford_circuit(4) # randomly generated clifford circuit\n",
        "for pair in circ:\n",
        "  cliff_state = across_qubits_gate(cliff_state, pair[1], clifford_gate_dict[pair[0]])\n",
        "\n",
        "haar_random_state = np.random.uniform(-1, 1, 2 ** n_qubits)\n",
        "haar_random_state = haar_random_state/np.linalg.norm(haar_random_state)\n",
        "\n",
        "global_mixture = [(0.8, cliff_state), (0.2, haar_random_state)]"
      ],
      "metadata": {
        "id": "HCnkA8ie-13X"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "The actual fidelity value is given as"
      ],
      "metadata": {
        "id": "NU8zcQl0_b9H"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "fid_val = 0\n",
        "for pair in global_mixture:\n",
        "  fid_val += pair[0] * np.abs(np.dot(np.conj(pair[1]), cliff_state)) ** 2\n",
        "\n",
        "print(\"Actual Fidelity: {}\".format(fid_val))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "8TUplLt9_l_I",
        "outputId": "1153cb14-2c98-4f6d-b7a6-4b17e3d1166e"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Actual Fidelity: 0.8020757052299634\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now, we shall use Classical shadows to estimate this quantity."
      ],
      "metadata": {
        "id": "T2je_C9X_6uP"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "shadows = ccs(global_mixture, 4000)"
      ],
      "metadata": {
        "id": "GVxYJatkABGE"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "compute_fidelity_with_shadows(shadows, circ, 50, 80)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "EpGTf_v7AJ1Z",
        "outputId": "016d528a-795d-4f37-805b-61c18a7de1ef"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.785"
            ]
          },
          "metadata": {},
          "execution_count": 58
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "As you can see, we have received a good estimate of the fidelity value."
      ],
      "metadata": {
        "id": "38DQv5shBZ4z"
      }
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "name": "classical_shadows_numpy.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}