%run init.ipynb

# tem algum problema com a funcao produto tensorial do sympy (implementar eu mesmo ...)
k000 = Matrix([1,0,0,0,0,0,0,0]); k001 = Matrix([0,1,0,0,0,0,0,0])
k010 = Matrix([0,0,1,0,0,0,0,0]); k011 = Matrix([0,0,0,1,0,0,0,0])
k100 = Matrix([0,0,0,0,1,0,0,0]); k101 = Matrix([0,0,0,0,0,1,0,0])
k110 = Matrix([0,0,0,0,0,0,1,0]); k111 = Matrix([0,0,0,0,0,0,0,1])
#k000,k001,k010,k011,k100,k101,k110,k111, k001*k001.T

p = symbols('p')
#p = 0
Psi0 = sqrt((4-3*p)/4)*k000 + sqrt(p/4)*(k101+k010+k111)
Psi1 = sqrt((4-3*p)/4)*k100 + sqrt(p/4)*(k001-k110-k011)
#Psi0.T, Psi1.T

r00,r01,r10,r11 = symbols('r_{00} r_{01} r_{10} r_{11}')
rhoA = Matrix([[r00,r01],[r10,r11]]); rhoA, rhoA[0,0]
#rhoA = Matrix([[2/3,1/3],[1/3,1/3]]); #rhoA, rhoA[0,0]

def rhoABt_s(rhoA,p):
    Psi0 = sqrt((4-3*p)/4)*k000 + sqrt(p/4)*(k101+k010+k111)
    Psi1 = sqrt((4-3*p)/4)*k100 + sqrt(p/4)*(k001-k110-k011)
    return rhoA[0,0]*Psi0*Psi0.T + rhoA[0,1]*Psi0*Psi1.T + rhoA[1,0]*Psi1*Psi0.T + rhoA[1,1]*Psi1*Psi1.T
rhoABt_ = rhoABt_s(rhoA,p); rhoABt_
# não foi possivel diagonalizar com sympy

def ptraceB(da, db, rho):
    rhoA = zeros(da,da)
    for j in range(0, da):
        for k in range(0, da):
            for l in range(0, db):
                rhoA[j,k] += rho[j*db+l,k*db+l]
    return rhoA

rhoAt = ptraceB(2, 4, rhoABt_); simplify(rhoAt) # ok!

rhoA = Matrix([[2/3,1/3],[1/3,1/3]])
p = np.arange(0,1.1,0.1); N = len(p)
Cl1 = np.zeros(N); Pjb = np.zeros(N)
for j in range(0,N):
    rhoABt_ = rhoABt_s(rhoA,p[j])
    rhoA_ = ptraceB(2, 4, rhoABt_)
    Cl1[j] = coh_l1_s(rhoA_)
    Pjb[j] = predict_jb_s(rhoA_)

# calculo feito a partir de rhoAB_til
import matplotlib
matplotlib.rcParams.update({'font.size':12})
plt.figure(figsize = (6,4), dpi = 100)
plt.plot(p,Cl1,label=r'$C_{l_1}$')
plt.plot(p,Pjb,label=r'$P_{jb}$')
plt.xlabel(r'$p$')
plt.legend()
plt.show()

evv = rhoABt.eigenvects(); evv

def tp_vv_s(psi, csi): # tensor product, of vectors, symbolic
    M = psi.shape[0]; N = csi.shape[0]
    eta = zeros(M*N,1)
    for j in range(0,M):
        for k in range(0,N):
            eta[j*N+k] = psi[j]*csi[k]
    return eta

def cb(d,j): # estados da base computacional
    v = zeros(d,1)
    v[j] = 1
    return v
cb(2,0)

def PhiABCt_s(rhoA,p):
    rhoABt = rhoABt_s(rhoA,p)
    eig = rhoABt.eigenvects()
    d = rhoABt.shape[0]; Phi = zeros(d*d,1)
    ne = 0; j = 0; l = -1
    while ne < d:
        mult = eig[j][1]; ne += mult
        for k in range(0,mult):
            l += 1
            Phi += sqrt(abs(eig[j][0]))*tp_vv_s(eig[j][2][k],cb(d,l))
        j += 1
    for j in range(0,d*d):
        if im(Phi[j]) < 10**-5:
            Phi[j] = re(Phi[j])
    return Phi

def coh_l1_s(rho):
    d = rho.shape[0]; C = 0
    for j in range(0,d-1):
        for k in range(j+1,d):
            C += abs(rho[j,k])
    return 2*C

def predict_jb_s(rho):
    return abs(rho[0,0]-rho[1,1])

def proj_s(psi): # simbolic projector
    d = psi.shape[0]
    proj = zeros(d,d)
    for j in range(0,d):
        for k in range(0,d):
            proj[j,k] = psi[j]*conjugate(psi[k])
    return proj

rhoA = Matrix([[2/3,1/3],[1/3,1/3]])
p = np.arange(0,1.1,0.1); N = len(p)
Cl1 = np.zeros(N); Pjb = np.zeros(N)
for j in range(0,N):
    Phi = PhiABCt_s(rhoA,p[j]); PPhi = proj_s(Phi)#; print(PPhi)
    rhoA_ = ptraceB(2, 2**5, PPhi)#; print(rhoA_[0,1])
    Cl1[j] = coh_l1_s(rhoA_)
    Pjb[j] = predict_jb_s(rhoA_)

# calculo feito a partir de PhiABC
import matplotlib
matplotlib.rcParams.update({'font.size':12})
plt.figure(figsize = (6,4), dpi = 100)
plt.plot(p,Cl1,label=r'$C_{l_1}$')
plt.plot(p,Pjb,label=r'$P_{jb}$')
plt.xlabel(r'$p$')
plt.legend()
plt.show()

def rho_AB_til(rhoA,p):
    rhoAbc = np.zeros((2**3,2**3), dtype=complex)#; print(rhoAbc)
    ket0 = np.array([[1],[0]]); ket1 = np.array([[0],[1]])#; print(ket0,ket1)
    ket00 = np.kron(ket0,ket0); ket01 = np.kron(ket0,ket1); ket10 = np.kron(ket1,ket0)
    ket11 = np.kron(ket1,ket1); #print(ket00,'',ket01,'',ket10,'',ket11)
    ket000 = np.kron(ket0,ket00); ket100 = np.kron(ket1,ket00)
    ket001 = np.kron(ket0,ket01); ket101 = np.kron(ket1,ket01)
    ket010 = np.kron(ket0,ket10); ket110 = np.kron(ket1,ket10)
    ket011 = np.kron(ket0,ket11); ket111 = np.kron(ket1,ket11)
    Psi0 = math.sqrt((4-3*p)/4)*ket000 + math.sqrt(p/4)*(ket101+ket010+ket111)
    Psi1 = math.sqrt((4-3*p)/4)*ket100 + math.sqrt(p/4)*(ket001-ket110-ket011)
    rhoAbc = rhoA[0,0]*Psi0@Psi0.T + rhoA[0,1]*Psi0@Psi1.T\
    + rhoA[1,0]*Psi1@Psi0.T + rhoA[1,1]*Psi1@Psi1.T
    return rhoAbc

def pTraceR_num(dl, dr, rhoLR):
    rhoL = np.zeros((dl, dl), dtype=complex)
    for j in range(0, dl):
        for k in range(j, dl):
            for l in range(0, dr):
                rhoL[j,k] += rhoLR[j*dr+l,k*dr+l]
        if j != k:
            rhoL[k,j] = np.conj(rhoL[j,k])
    return rhoL

rhoA = np.array([[2/3,1/3],[1/3,1/3]]); print(rhoA) # estado inicial
p = 0.
rhoAbc = rho_AB_til(rhoA,p)#; print(rhoAbc)
rhoA_ = pTraceR_num(2, 4, rhoAbc); print(rhoA_)

def coh_l1(rho):
    d = rho.shape[0]; 
    #d = rho.dims()[0]
    C = 0
    for j in range(0,d-1):
        for k in range(j+1,d):
            C += np.abs(rho[j,k])
    return 2*C

def predict_jb(rho):
    return abs(rho[0,0]-rho[1,1])

p = np.arange(0,1.1,0.1); #print(p)
N = len(p)#; print(N)
Cl1 = np.zeros(N); Pjb = np.zeros(N)
for j in range(0,N):
    rhoAbc = rho_AB_til(rhoA,p[j])
    rhoA_ = pTraceR_num(2, 4, rhoAbc)
    Cl1[j] = coh_l1(rhoA_)
    Pjb[j] = predict_jb(rhoA_)

# calculo feito a partir de rhoAB_til
import matplotlib
matplotlib.rcParams.update({'font.size':12})
plt.figure(figsize = (6,4), dpi = 100)
plt.plot(p,Cl1,label=r'$C_{l_1}$')
plt.plot(p,Pjb,label=r'$P_{jb}$')
plt.xlabel(r'$p$')
plt.legend()
plt.show()

from numpy import linalg
w, v = linalg.eig(rhoAbc)
print(v[0][:])
print(np.shape(v[:][0]))

ket0 = np.array([[1],[0]]); ket1 = np.array([[0],[1]])
ket00 = np.kron(ket0,ket0); ket01 = np.kron(ket0,ket1)
ket10 = np.kron(ket1,ket0); ket11 = np.kron(ket1,ket1)
ket000 = np.kron(ket0,ket00); ket100 = np.kron(ket1,ket00)
ket001 = np.kron(ket0,ket01); ket101 = np.kron(ket1,ket01)
ket010 = np.kron(ket0,ket10); ket110 = np.kron(ket1,ket10)
ket011 = np.kron(ket0,ket11); ket111 = np.kron(ket1,ket11)

p = np.arange(0,1.1,0.1)
N = len(p)
Cl1 = np.zeros(N); Pjb = np.zeros(N)
for j in range(0,N):
    rhoAbc = rho_AB_til(rhoA,p[j])
    w, v = linalg.eig(rhoAbc); w = np.abs(w)
    PhiAbcdef = math.sqrt(w[0])*np.kron(v.T[0],ket000) + math.sqrt(w[1])*np.kron(v.T[1],ket001)\
                + math.sqrt(w[2])*np.kron(v.T[2],ket010) + math.sqrt(w[3])*np.kron(v.T[3],ket011)\
                + math.sqrt(w[4])*np.kron(v.T[4],ket100) + math.sqrt(w[5])*np.kron(v.T[5],ket101)\
                + math.sqrt(w[6])*np.kron(v.T[6],ket110) + math.sqrt(w[7])*np.kron(v.T[7],ket111)
    rhoAbcdef = np.outer(PhiAbcdef,np.conj(PhiAbcdef))#; print(np.shape(rhoAbcdef))
    rhoA_ = pTraceR_num(2, 2**5, rhoAbcdef); print(rhoA_)#; print(np.shape(rhoA_))
    Cl1[j] = coh_l1(rhoA_)
    Pjb[j] = predict_jb(rhoA_)

# calculo feito a partir da purificacao rhoAB_til
import matplotlib
matplotlib.rcParams.update({'font.size':12})
plt.figure(figsize = (6,4), dpi = 100)
plt.plot(p,Cl1,label=r'$C_{l_1}$')
plt.plot(p,Pjb,label=r'$P_{jb}$')
plt.xlabel(r'$p$')
plt.legend()
plt.show()

rho_A = Matrix([[2/3,1/3],[1/3,1/3]])
rho_A.eigenvects()

rho_A*Matrix([[0.85],[0.52]])/0.87, rho_A*Matrix([[-0.5257],[0.85]])/0.127

w, v = linalg.eig(rhoA) # os autovetores são as colunas de v
print(w, v, v.T[0], v.T[1], np.shape(v.T[1])) # nao pode usar 

import qiskit
from qiskit import *
nshots = 8192
IBMQ.load_account()
provider = qiskit.IBMQ.get_provider(hub='ibm-q', group='open', project='main')
device = provider.get_backend('ibmq_quito')
simulator = Aer.get_backend('qasm_simulator')
from qiskit.tools.monitor import job_monitor
from qiskit_experiments.library import StateTomography

r00 = 2/3; r01 = 1.33/3; r10 = 1.33/3; r11 = 1/3 # initial state
r = math.sqrt((r00-r11)**2 + abs(2*r01)**2) # raio de Bloch
th = math.acos((r00-r11)/r)
ph = math.acos(re(2*r01)/(r*sin(th))) # angulos de Bloch
r0 = (1+r)/2.; r1 = (1-r)/2. # autovetores
print(r, th, ph, r0, r1)

pt = np.arange(0,1.01,0.01) # for the theoretical results
Ct = (1-pt)*(2*1.33/3)
Pt = (1-pt)*(1/3)

p = np.arange(0,1.1,0.1)
d = len(p)
Csim = np.zeros(d); Psim = np.zeros(d)

for j in range(0,d):
    pI = (4-3*p[j])/4; pX = p[j]/4; pZ = p[j]/4; pY = p[j]/4 # depolarizing
    # sequencia: 0000 0001 0010 0011 0100 0101 0110 0111 1000 1001 1010 1011 1100 1101 1110 1111  
    #           = 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
    Phi_ABCD = [math.sqrt(r0*pI)*math.cos(th/2), 
                math.sqrt(r0*pX)*(math.cos(ph)+1j*math.sin(ph))*math.sin(th/2), 
                math.sqrt(r0*pZ)*math.cos(th/2), 
                -1j*math.sqrt(r0*pY)*(math.cos(ph)+1j*math.sin(ph))*math.sin(th/2), 
                math.sqrt(r1*pI)*math.sin(th/2), 
                -math.sqrt(r1*pX)*(math.cos(ph)+1j*math.sin(ph))*math.cos(th/2), 
                math.sqrt(r1*pZ)*math.sin(th/2), 
                1j*math.sqrt(r1*pY)*(math.cos(ph)+1j*math.sin(ph))*math.cos(th/2), 
                math.sqrt(r0*pI)*(math.cos(ph)+1j*math.sin(ph))*math.sin(th/2), 
                math.sqrt(r0*pX)*math.cos(th/2), 
                -math.sqrt(r0*pZ)*(math.cos(ph)+1j*math.sin(ph))*math.sin(th/2), 
                1j*math.sqrt(r0*pY)*math.cos(th/2), 
                -math.sqrt(r1*pI)*(math.cos(ph)+1j*math.sin(ph))*math.cos(th/2), 
                math.sqrt(r1*pX)*math.sin(th/2), 
                math.sqrt(r1*pZ)*(math.cos(ph)+1j*math.sin(ph))*math.cos(th/2), 
                1j*math.sqrt(r1*pY)*math.sin(th/2)]
    qr = QuantumRegister(4) 
    qc = QuantumCircuit(qr)
    qc.initialize(Phi_ABCD, [qr[3],qr[2],qr[1],qr[0]])
    job_sim = StateTomography(qc, measurement_qubits = [0])
    data = job_sim.run(simulator, shots=nshots).block_for_results()
    rho_sim = data.analysis_results(0).value
    rho = rho_sim.to_operator().data
    Csim[j] = coh_l1(rho)
    Psim[j] = predict_jb(rho)

p = np.arange(0,1.1,0.1); d = len(p)
Cexp = np.zeros(d); Pexp = np.zeros(d)

for j in range(0,d):
    pI = (4-3*p[j])/4; pX = p[j]/4; pZ = p[j]/4; pY = p[j]/4
    Phi_ABCD = [math.sqrt(r0*pI)*math.cos(th/2), 
                math.sqrt(r0*pX)*(math.cos(ph)+1j*math.sin(ph))*math.sin(th/2), 
                math.sqrt(r0*pZ)*math.cos(th/2), 
                -1j*math.sqrt(r0*pY)*(math.cos(ph)+1j*math.sin(ph))*math.sin(th/2), 
                math.sqrt(r1*pI)*math.sin(th/2), 
                -math.sqrt(r1*pX)*(math.cos(ph)+1j*math.sin(ph))*math.cos(th/2), 
                math.sqrt(r1*pZ)*math.sin(th/2), 
                1j*math.sqrt(r1*pY)*(math.cos(ph)+1j*math.sin(ph))*math.cos(th/2), 
                math.sqrt(r0*pI)*(math.cos(ph)+1j*math.sin(ph))*math.sin(th/2), 
                math.sqrt(r0*pX)*math.cos(th/2), 
                -math.sqrt(r0*pZ)*(math.cos(ph)+1j*math.sin(ph))*math.sin(th/2), 
                1j*math.sqrt(r0*pY)*math.cos(th/2), 
                -math.sqrt(r1*pI)*(math.cos(ph)+1j*math.sin(ph))*math.cos(th/2), 
                math.sqrt(r1*pX)*math.sin(th/2), 
                math.sqrt(r1*pZ)*(math.cos(ph)+1j*math.sin(ph))*math.cos(th/2), 
                1j*math.sqrt(r1*pY)*math.sin(th/2)]
    qr = QuantumRegister(4); qc = QuantumCircuit(qr)
    qc.initialize(Phi_ABCD, [qr[3],qr[2],qr[1],qr[0]])
    qcst = StateTomography(qc, measurement_qubits = [0])
    data = qcst.run(device)
    print(data.experiment_id)
    rho = data.block_for_results().analysis_results(0).value
    rhoM = rho.to_operator().data
    Cexp[j] = coh_l1(rhoM)
    Pexp[j] = predict_jb(rhoM)
    print(Cexp,Pexp)

import matplotlib
matplotlib.rcParams.update({'font.size':12})
plt.figure(figsize = (6,4), dpi = 100)
plt.plot(pt,Ct,label=r'$C_{l_1}^{the}$')
plt.plot(pt,Pt,label=r'$P_{jb}^{the}$')
plt.plot(p,Csim,'*',label=r'$C_{l_1}^{sim}$')
plt.plot(p,Psim,'o',label=r'$P_{jb}^{sim}$')
plt.plot(p,Cexp,'^',label=r'$C_{l_1}^{exp}$')
plt.plot(p,Pexp,'+',label=r'$P_{jb}^{exp}$')
plt.xlabel(r'$p$')
plt.legend()
plt.show()

# para o chip quito, sem mitigacao de erro

# para o chip manila sem mitigacao de erro