PythonCH_Clifford_simulator.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import random\n",
    "import math\n",
    "import itertools\n",
    "from copy import deepcopy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def CliffordUpdate(psi,gate,qubits,n):\n",
    "    \"\"\"psi is a stabilizer state in CH form\n",
    "    gate is 'S','CZ','CX',or 'H'\n",
    "    qubits=[q] or qubits=[q1 q2] are the qubits the gate acts on\n",
    "    n is the total number of qubits\n",
    "    psi is a dict containing np.arrays/integers {n,s,r,Uc and p}\"\"\"\n",
    "    \n",
    "    if psi['n']!=n:\n",
    "        raise AssertionError(\"Inconsistent number of qubits, n \", n, psi['n'])\n",
    "    else:\n",
    "        if gate=='S':\n",
    "            psi['Uc']=LeftMultS(psi['Uc'],qubits,n)\n",
    "        elif gate=='CZ':\n",
    "            psi['Uc']=LeftMultCZ(psi['Uc'],qubits[0],qubits[1],n)\n",
    "        elif gate=='CX':\n",
    "            psi['Uc']=LeftMultCNOT(psi['Uc'],qubits[0],qubits[1],n)\n",
    "        elif gate=='X':\n",
    "            P=psi['Uc'][qubits[0],:].copy() #Pauli U_c^{\\dagger} X_j U_c\n",
    "            # Now conjugate by H(r)\n",
    "\n",
    "            numy=np.mod(np.count_nonzero(P[0:n]*P[n:2*n]*psi['r']),2)#Count number of times Y-->-Y when applying H(r)\n",
    "            P[2*n]=np.mod(P[2*n]+2*numy,4)#adjust phase accordingly\n",
    "\n",
    "            inds=np.nonzero(psi['r'])[0] #Apply H(r): swap X and Z at appropriate indices\n",
    "\n",
    "            for j in inds:\n",
    "                holder=P[j].copy()\n",
    "                P[j]=P[n+j].copy()\n",
    "                P[n+j]=holder\n",
    "\n",
    "            # Now compute P|s>\n",
    "            v1=P[0:n]\n",
    "            w1=P[n:2*n]\n",
    "            alpha=np.mod(P[2*n]+np.sum(v1*w1*(1+2*psi['s']))+2*np.sum((1-v1)*w1*psi['s']),4)\n",
    "            psi['p']=np.mod(psi['p']+2*alpha,8)\n",
    "            psi['s']=np.mod(psi['s']+P[0:n],2)\n",
    "\n",
    "\n",
    "\n",
    "        if gate=='H':\n",
    "\n",
    "        #Takes a stabilizer state psi in CH form applies the single-qubit H gate to the\n",
    "        #specified qubit\n",
    "\n",
    "\n",
    "        #Define Paulis\n",
    "        # P1=H(r) U_c ^{\\dagger} X_qubit U_c H(r)\n",
    "        # = [v1 | w1 | b1] (v1 is n bits, w1 is n bits, and b1 in the set 0,1,2,3)\n",
    "        # P2=H(r) U_c ^{\\dagger} Z_qubit U_c H(r)\n",
    "        # = [v2 | w2 | b2]\n",
    "        # Then the state after applying H is\n",
    "        #        |phi>= w^p U_c H(r) (P1|s>+P2|s>)/sqrt(2)\n",
    "            qubit=qubits[0]\n",
    "\n",
    "            P1=psi['Uc'][qubit,:].copy()\n",
    "\n",
    "            numy=np.mod(np.count_nonzero(P1[0:n]*P1[n:2*n]*psi['r']),2)\n",
    "            P1[2*n]=np.mod(P1[2*n]+2*numy,4)\n",
    "\n",
    "            P2=psi['Uc'][n+qubit,:].copy()\n",
    "            numy2=np.mod(np.count_nonzero(P2[0:n]*P2[n:2*n]*psi['r']),2)\n",
    "            P2[2*n]=np.mod(P2[2*n]+2*numy2,4)\n",
    "\n",
    "            inds=np.nonzero(psi['r'])[0] #Apply H(r): swap X and Z at appropriate indices\n",
    "            for j in inds:\n",
    "                holder=P1[j].copy()\n",
    "                P1[j]=P1[n+j].copy()\n",
    "                P1[n+j]=holder\n",
    "\n",
    "                holder=P2[j].copy()\n",
    "                P2[j]=P2[n+j].copy()\n",
    "                P2[n+j]=holder\n",
    "\n",
    "            b1=P1[2*n]\n",
    "            v1=P1[0:n]\n",
    "            w1=P1[n:2*n]\n",
    "\n",
    "            b2=P2[2*n]\n",
    "            v2=P2[0:n]\n",
    "            w2=P2[n:2*n]\n",
    "\n",
    "    #         print(type(b1))\n",
    "    #         print(type(v1))\n",
    "    #         print(type(w1))\n",
    "    #         print(type(psi['s']))\n",
    "            alpha=np.mod(b1+np.sum(v1*w1*(1+2*psi['s']))+2*np.sum((1-v1)*w1*psi['s']),4)\n",
    "            beta=np.mod(b2+np.sum(v2*w2*(1+2*psi['s']))+2*np.sum((1-v2)*w2*psi['s']),4)\n",
    "\n",
    "            s1=np.mod(psi['s']+v1,2)\n",
    "    #         print('s1=',s1)\n",
    "            s2=np.mod(psi['s']+v2,2)\n",
    "    #         print('s2=',s2)\n",
    "            a=np.mod(beta-alpha,4)\n",
    "    #         print('a=',a)\n",
    "\n",
    "            y=np.mod(s1+s2,2)\n",
    "    #         print(y)\n",
    "    #         print(\"=y\")\n",
    "\n",
    "            if np.count_nonzero(y)==0:\n",
    "                if a==1:\n",
    "                    psi['p']=np.mod(psi['p']+2*alpha+1,8)\n",
    "                else:\n",
    "                    psi['p']=np.mod(psi['p']+2*alpha-1,8)\n",
    "\n",
    "                psi['s']=s1\n",
    "\n",
    "            if np.count_nonzero(y)>0:\n",
    "                #Two cases\n",
    "                klist=np.nonzero(np.mod(psi['r']+1,2)*y)[0]#Indices where r_k=0 and y_k=1\n",
    "    #             print(klist)\n",
    "    #             print(\"=klist\")\n",
    "                mlist=np.nonzero(psi['r']*y)[0]#Indices where r_k=1 and y_k=1\n",
    "    #             print(mlist)\n",
    "    #             print(\"=mlist\")\n",
    "\n",
    "                if klist.size!=0:# if klist is nonempty\n",
    "                    k=klist[0]\n",
    "                    psi['p']=np.mod(psi['p']+2*alpha+2*s1[k]*a,8)\n",
    "                    psi['r'][k]=1\n",
    "                    if s1[k]==0:\n",
    "                        psi['s']=s1.copy()\n",
    "                    else:\n",
    "                        psi['s']=s2.copy()\n",
    "\n",
    "\n",
    "                    for j in klist[1::]: #All except k=klist[0]\n",
    "                        psi['Uc']=RightMultCNOT(psi['Uc'],k,j,n)\n",
    "\n",
    "                    for j in mlist:\n",
    "                        psi['Uc']=RightMultCZ(psi['Uc'],k,j,n)\n",
    "\n",
    "                    if s1[k]==0:\n",
    "    #                     print(\"a=\")\n",
    "    #                     print(a)\n",
    "                        for j in range(1,a+1):\n",
    "    #                         print(\"j=\")\n",
    "    #                         print(j)\n",
    "                            psi['Uc']=RightMultS(psi['Uc'],k,n)\n",
    "                    else:\n",
    "                        a2=np.mod(-a,4)\n",
    "    #                     print(\"a2=\")\n",
    "    #                     print(a2)\n",
    "                        for j in range(1,a2+1):\n",
    "    #                         print(\"j=\")\n",
    "    #                         print(j)\n",
    "                            psi['Uc']=RightMultS(psi['Uc'],k,n)\n",
    "    #                         print(\"psi['Uc']=\")\n",
    "    #                         print(psi['Uc'])\n",
    "\n",
    "\n",
    "                if klist.size==0:# if klist is empty\n",
    "                    k=mlist[0]\n",
    "\n",
    "                    # Update phase p\n",
    "                    if a==1:\n",
    "                        psi['p']=np.mod(psi['p']+2*alpha+1,8)\n",
    "                    if a==3:\n",
    "                        psi['p']=np.mod(psi['p']+2*alpha-1,8)\n",
    "                    if (a==0 or a==2):\n",
    "                        psi['p']=np.mod(psi['p']+2*alpha,8)\n",
    "\n",
    "\n",
    "                    B=np.nonzero(y)[0]#0-indexed in python\n",
    "\n",
    "                    # Update Uc\n",
    "                    if s1[k]==1:\n",
    "                        psi['Uc']=RightMultS(psi['Uc'],k,n)#Apply Z_k gate\n",
    "                        psi['Uc']=RightMultS(psi['Uc'],k,n)\n",
    "                    for j in B:\n",
    "    #                     print('B=',B)\n",
    "    #                     print('j=',j)\n",
    "    #                     print('k=',k)\n",
    "                        if j != k:\n",
    "                            psi['Uc']=RightMultCNOT(psi['Uc'],j,k,n)\n",
    "\n",
    "                    if (a==1 or a==3):\n",
    "                        psi['Uc']=RightMultS(psi['Uc'],k,n)\n",
    "\n",
    "\n",
    "                    # Update r\n",
    "    #                 print('r_in=',psi['r'])\n",
    "    #                 print('k=',k)\n",
    "                    if (a==0 or a==2):\n",
    "    #                     print('r_to_update=',psi['r'])\n",
    "    #                     print('k-th=',[k])\n",
    "    #                     print('update r')\n",
    "                        psi['r'][k]=np.mod(psi['r'][k]+1,2)\n",
    "    #                 print('r_out=',psi['r'])\n",
    "\n",
    "\n",
    "                    # Update s\n",
    "                    psi['s']=s1\n",
    "\n",
    "                    if (s1[k]==1 and (a==1 or a==2) or (s1[k]==0 and (a==0 or a==3))):\n",
    "                        pass\n",
    "                    else:\n",
    "                        psi['s'][k]=np.mod(psi['s'][k]+1,2)\n",
    "            \n",
    "            \n",
    "    return psi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def AmplitudeUpdates(Pauli, psi, j, n):\n",
    "    \"\"\"Compute amplitude \n",
    "    amp=w^p <0|U_c^{\\dagger} X_j U_c Pauli*H(r)|s>  where j is a bit index\n",
    "    Here eps=0,1\n",
    "    m=-n,...0\n",
    "    p =0, 1,2, ...,7 \n",
    "    w=exp(i*pi/4)\n",
    "\n",
    "    Note <x|psi>=<0|U_c^{\\dagger} X(x) U_c H(r)|s>w^p\n",
    "\n",
    "\n",
    "    X=[x zeros(1,n) 0];\n",
    "    P=PauliImage(X,psi.Uc,n);\"\"\"\n",
    "\n",
    "    P = MultiplyPauli(psi['Uc'][j,:], Pauli, n)\n",
    "    v1 = P[0:n]\n",
    "    eps = 1\n",
    "    if np.count_nonzero(np.mod(v1*(1-psi['r'])+psi['s']*(1-psi['r']),2)) > 0:    \n",
    "        #v1 and s should only differ where r is 1\n",
    "        eps = 0\n",
    "        m = 0\n",
    "        p = 0\n",
    "        amp = 0\n",
    "    else:\n",
    "        p = np.mod(psi['p']-2*np.count_nonzero(P[n:2*n]*P[0:n])+2*P[2*n], 8)     \n",
    "        #Number of Y paulis in P is nnz(P(n+1:2*n).*P(1:n)), each contributes -i\n",
    "        m = np.count_nonzero(psi['r'])\n",
    "\n",
    "        # Now compute sign(<v1|H(r)|s>)\n",
    "        p = np.mod(p+4*np.mod(np.sum(v1*psi['s']*psi['r']), 2), 8)    \n",
    "        #Look at places where r=v1=s=1. Each contributes a minus sign\n",
    "\n",
    "        amp = 2**(-m/2)*eps*np.exp(1j*p*np.pi/4)\n",
    "    return amp\n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Amplitude(x, psi, n):\n",
    "    \"\"\"Compute amplitude amp=<x|psi>=2^(-m/2)*eps* w^{p)  where x is an n-bit string and psi is a\n",
    "    stabilizer state in CH form\n",
    "\n",
    "    Here eps=0,1\n",
    "    m=-n,...0\n",
    "    p =0, 1,2, ...,7 \n",
    "    w=exp(i*pi/4)\n",
    "\n",
    "    Note <x|psi>=<0|U_c^{\\dagger} X(x) U_c H(r)|s>w^p\"\"\"\n",
    "\n",
    "    X = np.hstack((x,np.zeros(n,dtype=int),0))\n",
    "    P = PauliImage(X, psi['Uc'], n)\n",
    "    v1 = P[0:n]\n",
    "    eps = 1\n",
    "    if np.count_nonzero(np.mod(v1*(1-psi['r'])+psi['s']*(1-psi['r']),2)) > 0:    \n",
    "        #v1 and s should only differ where r is 1\n",
    "        eps = 0\n",
    "        m = 0\n",
    "        p = 0\n",
    "        amp = 0\n",
    "    else:\n",
    "        p = np.mod(psi['p']-2*np.count_nonzero(P[n:2*n]*P[0:n])+2*P[2*n], 8) \n",
    "        #Number of Y paulis in P is nnz(P(n+1:2*n).*P(1:n)), each contributes -i\n",
    "        m = np.count_nonzero(psi['r'])\n",
    "\n",
    "        # Now compute sign(<v1|H(r)|s>)\n",
    "        p = np.mod(p+4*np.mod(np.sum(v1*psi['s']*psi['r']), 2), 8)    \n",
    "        #Look at places where r=v1=s=1. Each contributes a minus sign\n",
    "\n",
    "        amp = 2**(-m/2)*eps*np.exp(1j*p*np.pi/4)\n",
    "    return amp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def CompBasisVector(z):\n",
    "    \"\"\"constructs a computational basis vector in C-H form\n",
    "    |psi>= w^p U_c H(r)|s>\n",
    "    where\n",
    "    s is a computational basis vector\n",
    "    H(r) is hadamards applied to qubits in r \n",
    "    U_c is a Clifford such that U_c|0>=|0> (stored by its tableau). \n",
    "    tableau format: first n columns are U^{\\dagger} X_j U plus phase bit, last\n",
    "    n columns are U^{\\dagger} Z_j U plus phase bit\n",
    "    p is an integer,  w=exp(i\\pi/4).\n",
    "    Here z vector of length n,  and basis is either 0 (z basis) or 1 (x basis)\"\"\"\n",
    "\n",
    "    if np.ndim(z)!=1:\n",
    "        z=np.squeeze(z) #Make z into a 1-d array by removing singletons\n",
    "\n",
    "    n = np.size(z)\n",
    "    psi={\n",
    "        'n':n,\n",
    "        's':np.asarray(z),\n",
    "        #'s':z,\n",
    "        'r':np.zeros(n,dtype=int),\n",
    "        'Uc':np.concatenate((np.eye(2*n,dtype=int), np.zeros([2*n,1],dtype=int)), axis=1),\n",
    "        'p':0,\n",
    "        'c':np.exp(0*1j*np.pi)\n",
    "    }\n",
    "    return psi\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def XBasisVector(z):\n",
    "    \"\"\"constructs an X-basis vector in C-H form\n",
    "    |psi>= w^p U_c H(r)|s>\n",
    "    where\n",
    "    s is a computational basis vector\n",
    "    H(r) is hadamards applied to qubits in r \n",
    "    U_c is a Clifford such that U_c|0>=|0> (stored by its tableau). \n",
    "    tableau format: first n columns are U^{\\dagger} X_j U plus phase bit, last\n",
    "    n columns are U^{\\dagger} Z_j U plus phase bit\n",
    "    p is an integer,  w=exp(i\\pi/4).\n",
    "    Here z vector of length n,  and basis is either 0 (z basis) or 1 (x basis)\"\"\"\n",
    "\n",
    "    if np.ndim(z)!=1:\n",
    "        z=np.squeeze(z) #Make z into a 1-d array by removing singletons\n",
    "\n",
    "    n = np.size(z)\n",
    "    psi={\n",
    "        'n':n,\n",
    "        's':np.asarray(z),\n",
    "        #'s':z,\n",
    "        'r':np.ones(n,dtype=int),\n",
    "        'Uc':np.concatenate((np.eye(2*n,dtype=int), np.zeros([2*n,1],dtype=int)), axis=1),\n",
    "        'p':0,\n",
    "        'c':np.exp(0*1j*np.pi)\n",
    "    }\n",
    "    return psi\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def RightMultCNOT(tableau, control, target, n):\n",
    "    \"\"\"Takes a tableau for a Clifford U and outputs the tableau for\n",
    "    U*CNOT_{control,target}\n",
    "     XI --> XX\n",
    "     YI --> YX\n",
    "     IZ --> ZZ\n",
    "     IY-->  ZY \n",
    "     XX --> XI\n",
    "     ZZ--> IZ\n",
    "     YY--> -XZ\n",
    "     XY--> YZ\n",
    "     YX--> YI\n",
    "     XZ--> -YY\n",
    "     YZ--> XY\n",
    "     ZY-->IY\n",
    "    tableau format: first n columns are U^{\\dagger} X_j U plus phase bit, \n",
    "    last n columns are U^{\\dagger} Z_j U plus phase bit\"\"\"\n",
    "\n",
    "    #Look at the columns of the tableau corresponding to qubits control/target\n",
    "\n",
    "    for j in range(2*n):\n",
    "        a = np.array([tableau[j, control], tableau[j, control + n], tableau[j, target], tableau[j, target + n]],dtype=int)\n",
    "        b = 1*a[0] + 2*a[1] + 4*a[2] + 8*a[3]\n",
    "        if b == 1:        # [1 0 0 0] % XI-->XX\n",
    "            tableau[j, target] = 1\n",
    "        elif b == 3:        #[1 1 0 0] % YI-->YX\n",
    "            tableau[j, target] = 1\n",
    "        elif b == 8:        # [0 0 0 1] %IZ-->ZZ\n",
    "            tableau[j, control + n] = 1\n",
    "        elif b == 12:        #[0 0 1 1] %IY-->ZY\n",
    "            tableau[j, control + n] = 1\n",
    "        elif b == 5:        # [1 0 1 0] % XX-->XI\n",
    "            tableau[j, target] = 0\n",
    "        elif b == 10:        #[0 1 0 1] %ZZ-->IZ\n",
    "            tableau[j, control + n] = 0\n",
    "        elif b == 15:        #[1 1 1 1] %YY-->-XZ\n",
    "            tableau[j, control + n] = 0\n",
    "            tableau[j, target] = 0\n",
    "            tableau[j, -1] = np.mod(tableau[j,-1] + 2, 4)\n",
    "        elif b == 13:        #[1 0 1 1] %XY-->YZ\n",
    "            tableau[j, control + n] = 1\n",
    "            tableau[j, target] = 0\n",
    "        elif b == 7:        # [1 1 1 0] %YX-->YI\n",
    "            tableau[j, target] = 0\n",
    "        elif b == 9:        # [1 0 0 1] %XZ-->-YY\n",
    "            tableau[j, control + n] = 1\n",
    "            tableau[j, target] = 1\n",
    "            tableau[j, -1] = np.mod(tableau[j, -1] + 2, 4)\n",
    "        elif b == 11:        #[1 1 0 1] %YZ-->XY\n",
    "            tableau[j, control + n] = 0\n",
    "            tableau[j, target] = 1\n",
    "        elif b == 14:        #[0 1 1 1] %ZY-->IY\n",
    "            tableau[j, control + n] = 0\n",
    "        else:\n",
    "            pass\n",
    "        \n",
    "    return tableau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def RightMultCZ(tableau, control, target, n):\n",
    "    \"\"\" Takes a tableau for a Clifford U and outputs the tableau for            \n",
    " U*CZ_{control,target}                                                   \n",
    " XI --> XZ                                                               \n",
    " YI --> YZ                                                               \n",
    " IY-->  ZY                                                               \n",
    " IX--> ZX                                                                \n",
    "                                                                         \n",
    " XX --> YY                                                               \n",
    " YY--> XX                                                                \n",
    "                                                                         \n",
    " XY--> -YX                                                               \n",
    " YX--> -XY                                                               \n",
    "                                                                         \n",
    " XZ--> XI                                                                \n",
    " ZX-->IX                                                                 \n",
    " YZ--> YI                                                                \n",
    " ZY-->IY                                                                 \n",
    "\"\"\"\n",
    "\n",
    "    #Look at the columns of the tableau corresponding to qubits control/target\n",
    "\n",
    "    for j in range(2*n):\n",
    "        a = np.array([tableau[j, control], tableau[j, control + n], tableau[j, target], tableau[j, target + n]],dtype=int)\n",
    "        b = 1*a[0] + 2*a[1] + 4*a[2] + 8*a[3]\n",
    "        if b == 1:        #  [1 0 0 0] % XI-->XZ\n",
    "            tableau[j,n + target] = 1\n",
    "        elif b == 3:        #[1 1 0 0] % YI-->YZ\n",
    "            tableau[j,n + target] = 1\n",
    "        elif b == 12:        # [0 0 1 1] %IY-->ZY\n",
    "            tableau[j, control + n] = 1\n",
    "        elif b == 4:        #[0 0 1 0] %IX-->ZX\n",
    "            tableau[j, control + n] = 1\n",
    "        elif b == 5:        # [1 0 1 0] % XX--> YY\n",
    "            tableau[j, control + n] = 1\n",
    "            tableau[j, target + n] = 1\n",
    "        elif b == 15:        #[1 1 1 1] %YY-->XX\n",
    "            tableau[j, control + n] = 0\n",
    "            tableau[j, target + n] = 0\n",
    "        elif b == 13:        #[1 0 1 1] %XY-->-YX\n",
    "            tableau[j, control + n] = 1\n",
    "            tableau[j, target + n] = 0\n",
    "            tableau[j, -1] = np.mod(tableau[j, -1] + 2, 4)\n",
    "        elif b == 7:        #[1 1 1 0] %YX-->-XY\n",
    "            tableau[j, control + n] = 0\n",
    "            tableau[j, target + n] = 1\n",
    "            tableau[j, -1] = np.mod(tableau[j, -1] + 2, 4)\n",
    "        elif b == 9:        #[1 0 0 1] %XZ-->XI\n",
    "            tableau[j, target + n] = 0\n",
    "        elif b == 6:        #[0 1 1 0] %ZX-->IX\n",
    "            tableau[j, control + n] = 0\n",
    "        elif b == 11:        #[1 1 0 1] %YZ-->YI\n",
    "            tableau[j, target + n] = 0\n",
    "        elif b == 14:        #[0 1 1 1] %ZY-->IY\n",
    "            tableau[j, control + n] = 0\n",
    "        else:\n",
    "            pass\n",
    "        \n",
    "    return tableau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def LeftMultCNOT(tableau, control, target, n):\n",
    "    \"\"\"Takes a tableau for a Clifford U and outputs the tableau for\n",
    "    CNOT_{control,target}*U\n",
    "\n",
    "    X_control --> X_control X_target\n",
    "    Z_target --> Z_control Z_target\"\"\"\n",
    "\n",
    "    newtableau = tableau.copy()\n",
    "\n",
    "    \n",
    "    Xprod = np.zeros(2 * n + 1,dtype=int)\n",
    "    Xprod[control] = 1\n",
    "    Xprod[target] = 1 #Pauli X_control*X_target\n",
    "\n",
    "    Zprod = np.zeros( 2 * n + 1,dtype=int)\n",
    "    Zprod[n + control] = 1\n",
    "    Zprod[n + target] = 1 #Pauli Z_control*Z_target\n",
    "\n",
    "    newtableau[control, :] = PauliImage(Xprod, tableau, n)\n",
    "    newtableau[n + target, :] = PauliImage(Zprod, tableau, n)\n",
    "    return newtableau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def LeftMultCZ(tableau, control, target, n):\n",
    "    \"\"\" Takes a tableau for a Clifford U and outputs the tableau for\n",
    "    CZ_{control,target}*U\n",
    "\n",
    "    X_control --> X_control Z_target\n",
    "    X_target --> Z_control X_target\"\"\"\n",
    "\n",
    "    newtableau = tableau.copy()\n",
    "\n",
    "    XZ = np.zeros(2 * n + 1,dtype=int)\n",
    "    XZ[control] = 1\n",
    "    XZ[target + n] = 1 #Pauli X_control*Z_target\n",
    "\n",
    "    ZX = np.zeros( 2 * n + 1,dtype=int)\n",
    "    ZX[n + control] = 1\n",
    "    ZX[target] = 1 #Pauli Z_control*X_target\n",
    "\n",
    "\n",
    "    newtableau[control, :] = PauliImage(XZ, tableau, n) #X_control --> U^{\\dagger} X_control Z_target U\n",
    "    newtableau[target, :] = PauliImage(ZX, tableau, n) #X_target -->U^{\\dagger} Z_control X_target U\n",
    "    return newtableau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def LeftMultS(tableau, qubit, n):\n",
    "    \"\"\" Takes a stabilizer tableau for a Clifford U and left multiplies the single-qubit S gate to the\n",
    "    specified qubit : U--> S_qubit*U\n",
    "\n",
    "    S^{\\dagger} X S=-Y\n",
    "    S^{\\dagger} Z S=Z \n",
    "\n",
    "    So must replace row qubit of tableau with image of -Y under U\"\"\"\n",
    "    newtableau = tableau.copy()\n",
    "    minusYj = np.zeros(2 * n + 1,dtype=int)\n",
    "    minusYj[2*n] = 2 #phase is -1\n",
    "    minusYj[qubit] = 1\n",
    "#     print(minusYj)\n",
    "#     print(type(_))\n",
    "#     print(qubit)\n",
    "    minusYj[qubit[0] + n] = 1\n",
    "    newtableau[qubit, :] = PauliImage(minusYj, tableau, n)# outputs -U^{\\dagger} Y_qubit U^{\\dagger}\n",
    "    return newtableau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def RightMultS(tableau, qubit, n):\n",
    "    \"\"\" Takes a stabilizer tableau for a Clifford U and left multiplies the single-qubit S gate to the\n",
    "    specified qubit : U--> U*S\n",
    "\n",
    "    Look at the columns of the tableau corresponding to qubit\"\"\"  \n",
    "\n",
    "\n",
    "\n",
    "    for j in range(2*n):\n",
    "        if tableau[j, qubit] == 1:        #X_qubit--> -Y_qubit; %Y_qubit--> X_qubit\n",
    "#             print(\"j=\")\n",
    "#             print(j)\n",
    "            tableau[j, qubit + n] = np.mod(tableau[j, qubit + n] + 1, 2)\n",
    "#             print(\"tableau[j, qubit]=\")\n",
    "#             print(tableau[j, qubit])\n",
    "            tableau[j, 2*n] = np.mod(tableau[j, 2*n] + 2 * tableau[j, qubit + n], 4)        #Flip phase if X-->-Y\n",
    "#             print(\"tableau[j, 2*n]=\")\n",
    "#             print(tableau[j, 2*n])\n",
    "            \n",
    "    return tableau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def PauliImage(P, tableau, n):\n",
    "    \"\"\"Takes as input a tableau for a Clifford U_c and a Pauli P and outputs\n",
    "    the pauli Uc^{\\dagger} P U_c\"\"\"\n",
    "\n",
    "\n",
    "    Pout = np.zeros(2 * n + 1,dtype=int)\n",
    "    Pout[-1] = P[-1].copy()\n",
    "    for j in range(n):\n",
    "        if P[j] == 1:        #Multiply by U_c X_j U_c\n",
    "            Pout = MultiplyPauli(Pout, tableau[j,:], n)\n",
    "\n",
    "        if P[n + j] == 1:        #Multiply by U_c Z_j U_c\n",
    "            Pout = MultiplyPauli(Pout, tableau[n + j,:], n)\n",
    "\n",
    "        if P[j] == 1 and P[n + j] == 1:        #Correct phase X_j*Z_j*(+i)=Y_j\n",
    "            Pout[-1] = np.mod(Pout[-1] + 1, 4)\n",
    "    return Pout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MultiplyPauli(P1, P2, n):\n",
    "    #Multiplies Paulis P3=P1*P2 as in Eq. 15.17 of Kitaev/Vyalyi/Shen book\n",
    "\n",
    "    P3 = np.mod(P1 + P2, 2)#Gets everything correct except phase bit\n",
    "    tau1 = 0\n",
    "    tau2 = 0\n",
    "    tau3 = 0\n",
    "    kap = 0\n",
    "\n",
    "    for j in range(n):\n",
    "        tau1 = tau1 + P1[j] * P1[n + j]\n",
    "        tau2 = tau2 + P2[j] * P2[n + j]\n",
    "        tau3 = tau3 + P3[j] * P3[n + j]\n",
    "        kap = np.mod(kap + P2[j] * P1[j + n], 2)\n",
    "\n",
    "    P3[-1] = np.mod(P1[-1] + P2[-1] + tau1 + tau2 - tau3 + 2 * kap, 4)\n",
    "    return P3\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bin_array(num, m):\n",
    "    \"\"\"Convert a positive integer num into an m-bit bit vector\"\"\"\n",
    "    return np.array(list(np.binary_repr(num).zfill(m))).astype(np.int8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Indicator(nparray, *length):\n",
    "    if not length:\n",
    "        n_values = np.max(np.asarray(nparray))+1\n",
    "        mat=np.eye(n_values,dtype=int)[nparray]\n",
    "    else:\n",
    "        n_values = length[0]+1\n",
    "        mat=np.eye(n_values,dtype=int)[nparray]\n",
    "    return mat  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def U_H(psi,i,n):\n",
    "    i+=1 #zero-indexing\n",
    "    h=np.asarray([[1,1],[1,-1]])*2**-0.5\n",
    "    H=np.kron(np.kron(np.eye(2**(i-1)),h),np.eye(2**(n-i)))\n",
    "    #print(H)\n",
    "    return H@psi\n",
    "\n",
    "def U_X(psi,i,n):\n",
    "    i+=1 #zero-indexing\n",
    "    x=np.asarray([[0,1],[1,0]])\n",
    "    X=np.kron(np.kron(np.eye(2**(i-1)),x),np.eye(2**(n-i)))\n",
    "    #print(X)\n",
    "    return X@psi\n",
    "\n",
    "def U_S(psi,i,n):\n",
    "    i+=1 #zero-indexing\n",
    "    s=np.asarray([[1,0],[0,1j]],dtype=complex)\n",
    "    S=np.kron(np.kron(np.eye(2**(i-1)),s),np.eye(2**(n-i)))\n",
    "    #print(S)\n",
    "    return S@psi\n",
    "\n",
    "\n",
    "def U_CX(psi,i,j,n):\n",
    "    i+=1 #zero-indexing\n",
    "    j+=1 #zero-indexing\n",
    "    # Makes local projector P_ij\n",
    "    swap=np.asarray([[1,0,0,0],[0,0,1,0],[0,1,0,0],[0,0,0,1]])\n",
    "    cnot=np.asarray([[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]])\n",
    "    if (j<i):\n",
    "        cnot=swap@cnot@swap\n",
    "    \n",
    "    i2=min([i,j])    \n",
    "    j2=max([i,j])\n",
    "    j=j2\n",
    "    i=i2\n",
    "    \n",
    "    CNOT=np.kron(np.kron(np.eye(2**(i-1)),cnot),np.eye(2**(n-i-1)))\n",
    "    \n",
    "    for k in range(i+1,j):\n",
    "        swk=np.kron(np.kron(np.eye(2**(k-1)),swap),np.eye(2**(n-k-1)))\n",
    "        CNOT=swk@CNOT@swk\n",
    "    \n",
    "    return CNOT@psi\n",
    "\n",
    "def U_CZ(psi,i,j,n):\n",
    "    i+=1 #zero-indexing\n",
    "    j+=1 #zero-indexing\n",
    "    # Makes local projector P_ij\n",
    "    swap=np.asarray([[1,0,0,0],[0,0,1,0],[0,1,0,0],[0,0,0,1]])\n",
    "    cz=np.asarray([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,-1]])\n",
    "    if (j<i):\n",
    "        cz=swap@cz@swap\n",
    "    \n",
    "    i2=min([i,j])    \n",
    "    j2=max([i,j])\n",
    "    j=j2\n",
    "    i=i2\n",
    "#     print('i,j=',i,j)\n",
    "    \n",
    "    CZ=np.kron(np.kron(np.eye(2**(i-1)),cz),np.eye(2**(n-i-1)))\n",
    "    \n",
    "#     print(list(range(i+1,j)))\n",
    "    for k in range(i+1,j):\n",
    "        swk=np.kron(np.kron(np.eye(2**(k-1)),swap),np.eye(2**(n-k-1)))\n",
    "        CZ=swk@CZ@swk\n",
    "#     print(CZ)\n",
    "    \n",
    "    return CZ@psi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def UnitaryUpdate(psi,gate,qubits,n):\n",
    "    \"\"\"psi is a stabilizer state in CH form\n",
    "    gate is 'S','CZ','CX',or 'H'\n",
    "    qubits=[q] or qubits=[q1 q2] are the qubits the gate acts on\n",
    "    n is the total number of qubits\n",
    "    psi is a dict containing np.arrays/integers {n,s,r,Uc and p}\"\"\"\n",
    "    if gate=='S':\n",
    "        psi=U_S(psi,qubits,n)\n",
    "    elif gate=='CZ':\n",
    "        psi=U_CZ(psi,qubits[0],qubits[1],n)\n",
    "    elif gate=='CX':\n",
    "        psi=U_CX(psi,qubits[0],qubits[1],n)\n",
    "    elif gate=='X':\n",
    "        psi=U_X(psi,qubits,n)\n",
    "    elif gate=='H':\n",
    "        psi=U_H(psi,qubits,n)        \n",
    "    \n",
    "    return psi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compUnitaryUpdate(psi,gate,qubits,n):\n",
    "    \"\"\"psi is a stabilizer state in CH form\n",
    "    gate is 'S','CZ','CX',or 'H'\n",
    "    qubits=[q] or qubits=[q1 q2] are the qubits the gate acts on\n",
    "    n is the total number of qubits\n",
    "    psi is a dict containing np.arrays/integers {n,s,r,Uc and p}\"\"\"\n",
    "    if gate=='S':\n",
    "        psi=U_S(psi,qubits[0],n)\n",
    "    elif gate=='CZ':\n",
    "        psi=U_CZ(psi,qubits[0],qubits[1],n)\n",
    "    elif gate=='CX':\n",
    "        psi=U_CX(psi,qubits[0],qubits[1],n)\n",
    "    elif gate=='X':\n",
    "        psi=U_X(psi,qubits[0],n)\n",
    "    elif gate=='H':\n",
    "        psi=U_H(psi,qubits[0],n)        \n",
    "    \n",
    "    return psi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fullvec(psi):\n",
    "    \"Write out psi in the computational basis using normal lexic ordering e.g. 000,001,...,111\"\n",
    "    \n",
    "    n=psi['n']\n",
    "    psivec=np.zeros((2**n,1),dtype=complex)\n",
    "    for j in range(2**n):\n",
    "        psivec[j]=Amplitude(bin_array(j, n),psi,n)\n",
    "    \n",
    "    if 'c' in psi:\n",
    "        psivec=psi['c']*psivec\n",
    "        \n",
    "    return psivec\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def normalize(v):\n",
    "    return v / np.linalg.norm(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "n=3\n",
    "#Make |100>. \n",
    "# Here the leftmost bit is the least significant\n",
    "# The qubits are labeled 0,1,2,3,...n-1 from left to right\n",
    "psi=CompBasisVector([1, 0, 0]); \n",
    "\n",
    "\n",
    "# Apply CNOT_{12} gate (1 is target, 2 is control)\n",
    "\n",
    "psi=CliffordUpdate(psi,'CX',[0, 1],n);\n",
    "# Apply CNOT_{23} gate \n",
    "psi=CliffordUpdate(psi,'CX',[1, 2],n);\n",
    "#State is now |111>\n",
    "\n",
    "#Apply Hadamards to qubits 2,3\n",
    " \n",
    "psi=CliffordUpdate(psi,'H',[1],n);   \n",
    "psi=CliffordUpdate(psi,'H',[2],n); \n",
    "\n",
    "# State is |1-->\n",
    "\n",
    "#Apply CZ_{12} and CZ_{13}\n",
    "psi=CliffordUpdate(psi,'CZ',[0, 1],n);\n",
    "psi=CliffordUpdate(psi,'CZ',[0, 2],n);\n",
    "\n",
    "# State is |1++>\n",
    "\n",
    "# Apply S gate to 1st qubit\n",
    "   \n",
    "psi=CliffordUpdate(psi,'S',[0],n);   \n",
    "# State is i*|1++>\n",
    "\n",
    "\n",
    "\n",
    "#Note that the ordering of bits is opposite to the tensor product\n",
    "#ordering. The state q is equal to \n",
    "v=(1/2)*1j*np.kron(np.kron(np.asarray([[0],[1]]),np.asarray([[1],[1]])),np.asarray([[1],[1]]))\n",
    "print(v)\n",
    "print(fullvec(psi))\n",
    "                         \n",
    "#Check q-v=0\n",
    "# norm(q-v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def update_and_compare(psi_in,phi_in,gtype,gqubits,n):\n",
    "    psi = CliffordUpdate(psi_in, gtype, gqubits, n)\n",
    "    phi = compUnitaryUpdate(phi_in, gtype, gqubits, n)\n",
    "    return psi, phi, np.linalg.norm(fullvec(psi) - phi)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def displayboth(psi,phi):\n",
    "    print(np.linalg.norm(fullvec(psi) - phi) )\n",
    "    print(np.squeeze(np.asarray([fullvec(psi),phi]),axis=2).T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def randomstabilizer(n):\n",
    "    psi=XBasisVector(np.zeros((n),dtype=int)); # Psi will be the state computed by the CH Clifford simulator\n",
    "    phi=np.ones((2**n,1),dtype=complex)*2**(-n/2)\n",
    "\n",
    "    m=100;   #Number of randomly chosen gates\n",
    "    r=np.random.randint(5, size=m);\n",
    "    qubits=[];\n",
    "    nrm=0;\n",
    "\n",
    "    for j in range(m):\n",
    "        if r[j] == 0:\n",
    "            q = np.random.randint(n)\n",
    "#             print('S',[q])\n",
    "            qubits = np.append(qubits, q)\n",
    "            psi, phi, nrm = update_and_compare(psi,phi,'S',[q],n)\n",
    "#             print('S',[q],nrm)\n",
    "\n",
    "        if r[j] == 1 or r[j] == 2:\n",
    "            collision = 0\n",
    "            while collision == 0:\n",
    "                q1 = np.random.randint(n)\n",
    "                q2 = np.random.randint(n)\n",
    "                if q1 != q2:\n",
    "                    collision = 1\n",
    "\n",
    "            qubits = np.append(qubits, [q1, q2])\n",
    "    #         print(qubits)\n",
    "            if r[j] == 1:\n",
    "#                 print('CZ', [q1, q2])\n",
    "                psi, phi, nrm = update_and_compare(psi,phi,'CZ', [q1, q2],n)\n",
    "#                 print('CZ', [q1, q2],nrm)\n",
    "\n",
    "            if r[j] == 2:\n",
    "#                 print('CX', [q1, q2])\n",
    "                psi, phi, nrm = update_and_compare(psi,phi,'CX', [q1, q2],n)\n",
    "#                 print('CX', [q1, q2],nrm)\n",
    "\n",
    "        if r[j] == 3:\n",
    "            q = np.random.randint(n)\n",
    "#             print('H',[q])\n",
    "            qubits = np.append(qubits, q)\n",
    "            psi, phi, nrm = update_and_compare(psi,phi,'H',[q],n)\n",
    "#             print('H',[q],nrm)\n",
    "\n",
    "        if r[j] == 4:\n",
    "            q = np.random.randint(n)\n",
    "#             print('X',[q])\n",
    "            qubits = np.append(qubits, q)\n",
    "            psi, phi, nrm = update_and_compare(psi,phi,'X',[q],n)\n",
    "#             print('X',[q],nrm)\n",
    "\n",
    "    q=fullvec(psi)\n",
    "\n",
    "    return psi, phi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psi, phi = randomstabilizer(4)\n",
    "displayboth(psi,phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def EquatorialA(n):\n",
    "    \"Randomly generate the A matrix corresponding to |phi_A> in Eq 61 \"\n",
    "    \n",
    "    A=np.zeros((n,n),dtype=np.int8)\n",
    "    offdiag=np.random.randint(2,size=int(n*(n-1)/2))\n",
    "    A[np.triu_indices(n, 1)]=offdiag\n",
    "    A=A+A.T\n",
    "    np.fill_diagonal(A,np.random.randint(4,size=n))\n",
    "    return A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ind2vec(ind, N=None):\n",
    "    ind = np.asarray(ind)\n",
    "    if N is None: \n",
    "        N = ind.max() + 1\n",
    "    return (np.arange(N) == ind[:,None]).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def binvec2dec(vec):\n",
    "    return np.dot((2**np.arange(vec.shape[0], dtype = np.uint64)[::-1]),vec).astype(np.int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Equatorialfullvec(A):\n",
    "    \"\"\"Write out psi in the computational basis using normal lexic ordering e.g. 000,001,...,111. \n",
    "    See Sec IV C of BBCCGH\"\"\"  \n",
    "    \n",
    "    n=A.shape[0]\n",
    "    psivec=np.zeros((2**n,1),dtype=complex)\n",
    "    for j in range(2**n):\n",
    "        x=bin_array(j,n)\n",
    "        psivec=psivec+1j**((x@A)@x.T)*ind2vec([binvec2dec(x)],2**n).T\n",
    "    psivec=2**(-n/2)*psivec\n",
    "        \n",
    "    return psivec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def InnerPsiA(psi,A):\n",
    "    \"\"\"Implementation of <phi|phi_A> as described in Lemma 3, Sec IV C of BBCCGH\"\"\"  \n",
    "    \n",
    "    n=psi['n']\n",
    "    r=psi['r']\n",
    "    s=psi['s']\n",
    "    gamma=psi['Uc'][0:n,-1]\n",
    "    F=psi['Uc'][0:n,0:n]\n",
    "    M=psi['Uc'][0:n,n:2*n]\n",
    "    G=psi['Uc'][n:2*n,n:2*n]\n",
    "    J=np.mod(M@F.T+np.diag(gamma),4)\n",
    "    mask=~np.eye(J.shape[0],dtype=bool)\n",
    "    J[mask]=np.mod(J[mask],2)\n",
    "    K=G.T@(A+J)@G\n",
    "    wtr=np.count_nonzero(r)\n",
    "    term1=2**(-(n+wtr)/2)\n",
    "    term2=1j**((s@K)@s.T)\n",
    "    term3=(-1)**np.dot(s,r)\n",
    "    phase=np.exp(-psi['p']*1j*np.pi/4)*psi['c'].conj()\n",
    "    locs=np.nonzero(r)[0]\n",
    "    sumtot=0\n",
    "    for j in range(2**(wtr)):\n",
    "        x=np.zeros(n,dtype=int)\n",
    "        x[locs]=bin_array(j, wtr)#n-bit strings x satisfying x_j <= r_j\n",
    "        sumtot=sumtot+1j**((x@K)@x.T+2*x@(s+s@K))\n",
    "    return phase*term1*term2*term3*sumtot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def checkInnerPsiA(psi,A):\n",
    "    aa=abs(InnerPsiA(psi,A))\n",
    "    bb= abs(np.vdot(fullvec(psi),Equatorialfullvec(A)))\n",
    "    return print([aa, bb, aa-bb])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def etaA(SUBSET,Amat):\n",
    "    n=Amat.shape[0]\n",
    "    olap=sum(list(map(lambda x: InnerPsiA(SUBSET[x],Amat), range(len(SUBSET)))))\n",
    "    return 2**n*abs(olap)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def eta_estimate(SUBSET,n,Lsamples):\n",
    "    return np.mean([etaA(SUBSET,EquatorialA(n)) for j in range(Lsamples)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def explicitetaA(SUBSET,Amat):\n",
    "    tot=0\n",
    "    for j in range(len(SUBSET)):\n",
    "        tot=tot+2**6*np.vdot(fullvec(SUBSET[j]),Equatorialfullvec(Amat))\n",
    "    return tot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def HiddenShiftCircuitBuilderCH(numtof,numqub,*hidstring):\n",
    "    \"\"\"Build the circuit for implimenting Hidden Shift for Bent functions\"\"\"\n",
    "    n=numqub\n",
    "    if not hidstring:\n",
    "        s=np.random.randint(2, size=n)\n",
    "    else:\n",
    "        s=np.asarray(hidstring)[0]\n",
    "    \n",
    "    q3=np.random.choice(int(n/2), 3, replace=False).copy()\n",
    "    Og=pd.DataFrame([{'H':[q3[2].copy()]},{'Toff':q3.copy()},{'H':[q3[2].copy()]}])\n",
    "    Og2=pd.DataFrame([{'H':[int(n/2)+q3[2].copy()]},{'Toff':int(n/2)+q3.copy()},{'H':[int(n/2)+q3[2].copy()]}])\n",
    "    for j in range(numtof-1):\n",
    "        for i in range(10):# The number of random Clifford gates Z or C-Z\n",
    "            r=np.random.randint(2)\n",
    "            if r==1:\n",
    "                #print('random Z')\n",
    "                q=np.random.randint(n/2)\n",
    "                Og=Og.append(pd.DataFrame([{'S':[q]},{'S':[q]}]),ignore_index=True) # Z gate applied to qubit q\n",
    "                Og2=Og2.append(pd.DataFrame([{'S':[int(n/2)+q]},{'S':[int(n/2)+q]}]),ignore_index=True)\n",
    "            elif r==0:\n",
    "                #print('random CZ')\n",
    "                q2=np.random.choice(int(60/2), 2, replace=False).copy()\n",
    "                Og=Og.append(pd.DataFrame([{'CZ':q2}]),ignore_index=True) # C-Z gate between qubits in q2\n",
    "                Og2=Og2.append(pd.DataFrame([{'CZ':int(n/2)+q2}]),ignore_index=True)\n",
    "        q3=np.random.choice(int(60/2), 3, replace=False).copy()\n",
    "        Og=Og.append(pd.DataFrame([{'H':[q3[2]]},{'Toff':q3},{'H':[q3[2]]}]),ignore_index=True)\n",
    "        Og2=Og2.append(pd.DataFrame([{'H':[int(n/2)+q3[2]]},{'Toff':int(n/2)+q3},{'H':[int(n/2)+q3[1]]}]),ignore_index=True)\n",
    "\n",
    "    Of=Og.copy()\n",
    "    Oftil=pd.DataFrame()\n",
    "    for j in range(int(n/2)):\n",
    "        Of=Of.append(pd.DataFrame([{'CZ':[j,j+int(n/2)]}]),ignore_index=True)\n",
    "        Oftil=Oftil.append(pd.DataFrame([{'CZ':[j,j+int(n/2)]}]),ignore_index=True)\n",
    "    Oftil=Oftil.append(Og2)\n",
    "\n",
    "    # This portion needs checking\n",
    "    Xs=pd.DataFrame()\n",
    "    FT=pd.DataFrame()\n",
    "    for i in range(n):\n",
    "        FT=FT.append(pd.DataFrame([{'H':[i]}]))\n",
    "        if s[i]==1:\n",
    "            Xs=Xs.append(pd.DataFrame(np.asarray([[[i],np.nan],[np.nan,[i]],[np.nan,[i]],[[i],np.nan]],dtype=object),columns=['H','S']),ignore_index=True)\n",
    "\n",
    "    U=FT.append([Of,Xs,FT,Oftil,FT],ignore_index=True)\n",
    "    return U, s "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def SamplesToTake(circ,delta):\n",
    "    numtof=0\n",
    "    numt=0\n",
    "    if 'CCZ' in circ.columns:\n",
    "        numtof=numtof+circ.count()['CCZ']\n",
    "    if 'Toff' in circ.columns:\n",
    "        numtof=numtof+circ.count()['Toff']\n",
    "    if 'T' in circ.columns:\n",
    "        numt=numt+circ.count()['T']\n",
    "    \n",
    "    print('numtof',numtof)\n",
    "    print('numt',numt)\n",
    "    c1T=np.cos(np.pi/8)**(-numt)\n",
    "    c1Toff=(3/4)**(-numtof)\n",
    "    ksamps=int(round((c1T*c1Toff/delta)**2))\n",
    "    #2**(-m/2)*eps*np.exp(1j*p*np.pi/4)\n",
    "    return ksamps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def CreateSubsetStates(circ,ksamps,psi_in):\n",
    "    cliffordgates={'H','S','CZ','CX'}\n",
    "\n",
    "\n",
    "    pathind=0\n",
    "\n",
    "    SUBSET=[]\n",
    "    for path in range(ksamps):\n",
    "        #print('path=',path)\n",
    "        #psi=XBasisVector([0,0,0])#create the |+++> state\n",
    "        psi=deepcopy(psi_in)\n",
    "        for ind, row in circ.iterrows():\n",
    "            tmp=row.dropna()\n",
    "            gtype=tmp.index[0]\n",
    "#             print(gtype)\n",
    "            gqubits=np.asarray(tmp.values[0])\n",
    "            if gtype in cliffordgates:\n",
    "                psi=CliffordUpdate(psi, gtype, gqubits, n)\n",
    "            elif gtype=='CCZ':\n",
    "                x=random.choices(list(range(0,8)), weights=np.ones(8)/6, k=1)[0]\n",
    "                #print('CCZ','x',x)\n",
    "                if x==0: # Identity\n",
    "                    pass\n",
    "                elif x==1: # CZ_{12}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,1]],n)\n",
    "                elif x==2: # CZ_{13}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,2]],n)\n",
    "                elif x==3: # CZ_{12,13}Z_{1}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,1]], n)\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[0]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[0]], n)\n",
    "                elif x==4: # CZ_{23}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[1,2]], n)\n",
    "                elif x==5: # CZ_{23,12}Z_{2}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[1,2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,1]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[1]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[1]], n)\n",
    "                elif x==6: # CZ_{13,23}Z_{3}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[1,2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[2]], n)\n",
    "                elif x==7: # -CZ_{12,13,23}Z_{1,2,3}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,1]], n)\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[1,2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[0]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[0]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[1]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[1]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[2]], n)\n",
    "                    psi['p']=np.mod(psi['p']+4,8) # This inserts the minus sign\n",
    "\n",
    "            elif gtype=='Toff':\n",
    "                psi=CliffordUpdate(psi, 'H', gqubits[[2]], n)\n",
    "                x=random.choices(list(range(0,8)), weights=np.ones(8)/6, k=1)[0]\n",
    "#                 x=feynmanpaths[pathind]\n",
    "#                 print('Toff','x',x)\n",
    "                pathind+=1\n",
    "                if x==0: # Identity\n",
    "                    pass\n",
    "                elif x==1: # CZ_{12}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,1]],n)\n",
    "                elif x==2: # CZ_{13}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,2]],n)\n",
    "                elif x==3: # CZ_{12,13}Z_{1}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,1]], n)\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[0]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[0]], n)\n",
    "                elif x==4: # CZ_{23}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[1,2]], n)\n",
    "                elif x==5: # CZ_{23,12}Z_{2}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[1,2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,1]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[1]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[1]], n)\n",
    "                elif x==6: # CZ_{13,23}Z_{3}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[1,2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[2]], n)\n",
    "                elif x==7: # -CZ_{12,13,23}Z_{1,2,3}\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,1]], n)\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[0,2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'CZ', gqubits[[1,2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[0]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[0]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[1]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[1]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[2]], n)\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits[[2]], n)\n",
    "                    psi['p']=np.mod(psi['p']+4,8) # This inserts the minus sign\n",
    "#                     psi['c']=(-1)*psi['c']\n",
    "                psi=CliffordUpdate(psi, 'H', gqubits[[2]], n)\n",
    "\n",
    "            elif gtype=='T':\n",
    "                x=random.choices(list(range(0,2)), weights=np.ones(2)/2, k=1)[0]\n",
    "                # T = c_0 Id + c_1 S where \n",
    "                # c_0 = 0.5*(1+1j)*(np.exp(1j*math.pi/4)-1j)\n",
    "                # c_1 = -0.5*(1+1j)*(np.exp(1j*math.pi/4)-1)\n",
    "                phi_0=np.angle(0.5*(1+1j)*(np.exp(1j*math.pi/4)-1j))\n",
    "                phi_1=np.angle(-0.5*(1+1j)*(np.exp(1j*math.pi/4)-1))\n",
    "                #print('T','x',x)\n",
    "                if x==0: # Identity\n",
    "                    pass\n",
    "                    psi['c']=np.exp(1j*phi_0)*psi['c']\n",
    "                elif x==1: # S\n",
    "                    psi=CliffordUpdate(psi, 'S', gqubits, n)\n",
    "                    psi['c']=np.exp(1j*phi_1)*psi['c']\n",
    "\n",
    "        SUBSET.append(psi)\n",
    "\n",
    "    return SUBSET"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "(circ,shift)=HiddenShiftCircuitBuilderCH(1,6,[0,0,0,1,1,1])\n",
    "print(shift)\n",
    "# (circ,shift)=HiddenShiftCircuitBuilderCH(1,6)\n",
    "# print(shift)\n",
    "\n",
    "SamplesToTake(circ,0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=shift.shape[0]\n",
    "def SumOverCliffords(circ):\n",
    "    print('n=',n)\n",
    "    ksamps=SamplesToTake(circ,0.1)\n",
    "    print(ksamps)\n",
    "    return CreateSubsetStates(circ,ksamps,CompBasisVector(np.asarray(np.zeros(n),dtype=int)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=3\n",
    "circ=pd.DataFrame([{'CCZ':[0,1,2]}])\n",
    "ksamps=SamplesToTake(circ,0.1)\n",
    "print('ksamps',ksamps)\n",
    "\n",
    "SUBSET=CreateSubsetStates(circ,ksamps,XBasisVector([0,0,0]))\n",
    "\n",
    "print(normalize(sum(list(map(fullvec,SUBSET)))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "circ=pd.DataFrame([{'Toff':[0,1,2]}])\n",
    "ksamps=SamplesToTake(circ,0.1)\n",
    "print('ksamps',ksamps)\n",
    "\n",
    "SUBSET=CreateSubsetStates(circ,ksamps,XBasisVector([0,0,0]))\n",
    "# SUBSET, phiSUBSET =DebugSubsetStates(circ,ksamps,CompBasisVector([1,1,1]))\n",
    "\n",
    "print(normalize(sum(list(map(fullvec,SUBSET)))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=3\n",
    "circ=pd.DataFrame([{'CCZ':[0,1,2]},{'Toff':[0,1,2]},{'T':[0]},{'T':[1]},{'T':[2]}])\n",
    "ksamps=SamplesToTake(circ,0.05)\n",
    "print('ksamps',ksamps)\n",
    "\n",
    "SUBSET=CreateSubsetStates(circ,ksamps,XBasisVector([0,0,0]))\n",
    "# SUBSET, phiSUBSET =DebugSubsetStates(circ,ksamps,CompBasisVector([1,1,1]))\n",
    "\n",
    "print(normalize(sum(list(map(fullvec,SUBSET)))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(circ,shift)=HiddenShiftCircuitBuilderCH(1,6,[0,0,0,1,1,1])\n",
    "print(shift)\n",
    "# (circ,shift)=HiddenShiftCircuitBuilderCH(1,6)\n",
    "# print(shift)\n",
    "\n",
    "SamplesToTake(circ,0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n=shift.shape[0]\n",
    "def SumOverCliffords(circ):\n",
    "    print('n=',n)\n",
    "    ksamps=SamplesToTake(circ,0.1)\n",
    "    print(ksamps)\n",
    "    return CreateSubsetStates(circ,ksamps,CompBasisVector(np.asarray(np.zeros(n),dtype=int)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "SUBSET=SumOverCliffords(circ)\n",
    "unnorm=sum(list(map(fullvec,SUBSET)))\n",
    "answer=normalize(unnorm)\n",
    "print(np.sum(np.round(answer)-fullvec(CompBasisVector(shift))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "A6= EquatorialA(6)\n",
    "print(A6)\n",
    "phiA=Equatorialfullvec(A6)\n",
    "print('<phi_A|phi_A>=',np.vdot(phiA,phiA))\n",
    "print('<phi_A|psi>=',np.vdot(unnorm,phiA))\n",
    "print('<psi|psi>=nrm_sq=',np.vdot(unnorm,unnorm))\n",
    "print('eta_estimate=',eta_estimate(SUBSET,6,100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}