%run init.ipynb

lbda = symbols('lambda', real=True); E0 = Matrix([[1,0],[0,sqrt(1-lbda)]]); E0

# coefficients of the expansion of E0 in the Pauli basis
c0 = trace((id(2)/sqrt(2))*E0); c1 = trace((pauli(1)/sqrt(2))*E0)
c2 = trace((pauli(2)/sqrt(2))*E0); c3 = trace((pauli(3)/sqrt(2))*E0)
simplify(c0), simplify(c1), simplify(c2), simplify(c3)

E1 = Matrix([[0,0],[0,sqrt(lbda)]]); E1

# coefficients for the expansion of E1 in the Pauli basis
c0 = trace((id(2)/sqrt(2))*E1); c1 = trace((pauli(1)/sqrt(2))*E1)
c2 = trace((pauli(2)/sqrt(2))*E1); c3 = trace((pauli(3)/sqrt(2))*E1)
simplify(c0), simplify(c1), simplify(c2), simplify(c3)

V00,V10,V20,V30,W00,W01,W02,W03,W10,W11,W12,W13=symbols('V00 V10 V20 V30 W00 W01 W02 W03 W10 W11 W12 W13',real=True)
#nonlinsolve([V00**2+V10**2+V20**2+V30**2-1, W00**2+W01**2+W02**2+W03**2-1,
#             W10**2+W11**2+W12**2+W13**2-1,W10*W00+W11*W01+W12*W02+W13*W03,
#             W00*V00-(1+sqrt(1-lbda))/2, W01*V10, W02*V20, W03*V30-(1-sqrt(1-lbda))/2,
#             W10*V00-sqrt(lbda)/2, W11*V10, W12*V20, W13*V30+sqrt(lbda)/2],
#            [V00,V10,V20,V30,W00,W01,W02,W03,W10,W11,W12,W13])
# took too long for solving (did not return any result)

V = Matrix([[sqrt((1+sqrt(1-lbda))/2),0,0,sqrt((1-sqrt(1-lbda))/2)],[0,1,0,0],[0,0,1,0],
            [sqrt((1-sqrt(1-lbda))/2),0,0,-sqrt((1+sqrt(1-lbda))/2)]])
V

V*V.T, V.T*V # verification of unitarity

W = Matrix([[0,1,0,0],[sqrt((1-sqrt(1-lbda))/2),0,0,sqrt((1+sqrt(1-lbda))/2)],[0,0,1,0],
            [sqrt((1+sqrt(1-lbda))/2),0,0,-sqrt((1-sqrt(1-lbda))/2)]])
W

W*W.T, W.T*W # verification of unitarity

r00,r01,r10,r11 = symbols('r00,r01,r10,r11',real=True)
rhoS0 = Matrix([[r00,r01],[r10,r11]])
rhoS0 # initial state

rhoS = E0*rhoS0*E0 + E1*rhoS0*E1; simplify(rhoS) # evolved state using the Kraus' operators

# Outside these functions, initialize: rhos = zeros(ds,ds), s=A,B
def ptraceA(da, db, rho):
    rhoB = zeros(db,db)
    for j in range(0, db):
        for k in range(0, db):
            for l in range(0, da):
                rhoB[j,k] += rho[l*db+j,l*db+k]
    return rhoB

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

# tests the circuit for the PDC
rhoSA0 = tp(rhoS0,proj(cb(4,0)))#; rhoSA0 
rhoSA1 = tp(id(2),V)*rhoSA0*tp(id(2),V.T)#; rhoSA1
Uc = tp(id(2),proj(cb(4,0))) + tp(id(2),proj(cb(4,1))) + tp(id(2),proj(cb(4,2))) + tp(pauli(3),proj(cb(4,3)))
#Uc*Uc.T, Uc.T*Uc
rhoSA2 = Uc*rhoSA1*Uc.T; #rhoSA2
rhoSA3 = tp(id(2),W)*rhoSA2*tp(id(2),W.T) # Não precisava pois só mudará a base de A
rhoS = ptraceB(2, 4, rhoSA3); simplify(rhoS) # ok!

V = Matrix([[sqrt((1+sqrt(1-lbda))/2),sqrt((1-sqrt(1-lbda))/2)],
            [sqrt((1-sqrt(1-lbda))/2),-sqrt((1+sqrt(1-lbda))/2)]])
W = V; V

# tests the (http://arxiv.org/abs/1704.05593) circuit
rhoSA0 = tp(rhoS0,proj(cb(2,0))); rhoSA0 
rhoSA1 = tp(id(2),V)*rhoSA0*tp(id(2),V.T); rhoSA1
Uc = tp(id(2),proj(cb(2,0))) + tp(pauli(3),proj(cb(2,1))); Uc, Uc*Uc.T, Uc.T*Uc
rhoSA2 = Uc*rhoSA1*Uc.T; rhoSA2
rhoSA3 = tp(id(2),W)*rhoSA2*tp(id(2),W.T)
rhoS = ptraceB(2, 2, rhoSA3); simplify(rhoS) # also works!

lb = np.arange(0,1.05,0.05); th = 2*np.arccos(np.sqrt((1+np.sqrt(1-lb))/2))
th_ = 2*np.arcsin(np.sqrt((1-np.sqrt(1-lb))/2))
plt.plot(lb,th); plt.plot(lb, th_, '*') 
plt.xlabel(r'$\lambda$'); plt.ylabel(r'$\theta$'); plt.show()

# phase damping channel, exact
def phase_damping(rho, lb):
    E0 = Matrix([[1,0],[0,sqrt(1-lb)]]); E1 = Matrix([[0,0],[0,sqrt(lb)]])
    return E0*rho*E0 + E1*rho*E1

lb = symbols('lambda')
psi = (cb(2,0)+cb(2,1))/sqrt(2); rho = proj(psi); rhopd = phase_damping(rho, lb)
rho, rhopd

def coh_l1(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

coh = coh_l1(rhopd); coh

# simulation
from qiskit import *

def state_prep():
    qr = QuantumRegister(1); qc = QuantumCircuit(qr, name = 'state_prep')
    qc.h(qr[0])
    return qc
qc_sp = state_prep(); qc_sp.draw()

def qc_sim_pd(lb):
    qr = QuantumRegister(2); qc = QuantumCircuit(qr, name = 'pd_sim')
    th = 2*math.acos(math.sqrt((1+math.sqrt(1-lb))/2))
    qc.u(th, 0, math.pi, qr[1]); qc.cz(qr[1],qr[0]); qc.u(th, 0, math.pi, qr[1])
    return qc
lb = 0.5; qc_pd = qc_sim_pd(lb); qc_pd.draw()

qr = QuantumRegister(2); qc = QuantumCircuit(qr)
qc_sp = state_prep(); qc.append(qc_sp, [qr[0]]); qc.barrier()
lb = 0; qc_pd = qc_sim_pd(lb); qc.append(qc_pd, [qr[0],qr[1]]); qc.barrier()
qc.draw()

from qiskit.ignis.verification.tomography import state_tomography_circuits, StateTomographyFitter
nshots = 8192
qiskit.IBMQ.load_account()
provider = qiskit.IBMQ.get_provider(hub='ibm-q', group='open', project='main')
device = provider.get_backend('ibmq_belem')
from qiskit.tools.monitor import job_monitor

lb = np.arange(0,1.05,0.05); d = lb.shape[0]; Cteo = np.zeros(d); Csim = np.zeros(d); Cexp = np.zeros(d)
for j in range(0, d):
    Cteo[j] = np.abs(math.sqrt(1-lb[j])) # theoretical
    qr = QuantumRegister(2); qc = QuantumCircuit(qr)
    qc_sp = state_prep(); qc.append(qc_sp, [qr[0]]) # state preparation
    qc_pd = qc_sim_pd(lb[j]); qc.append(qc_pd, [qr[0],qr[1]]) # apply sim phase damping 
    qstc = state_tomography_circuits(qc, qr[0]) # circuit for state tomography
    job = qiskit.execute(qstc, Aer.get_backend('qasm_simulator'), shots=nshots)
    qstf = StateTomographyFitter(job.result(), qstc); rho = qstf.fit(method='lstsq')
    Csim[j] = coh_l1(rho)
    # 1st I did the simulation, only afterwards I added the code for the experiments
    job = qiskit.execute(qstc, backend = device, shots = nshots)
    print(job.job_id()); job_monitor(job)
    qstf = StateTomographyFitter(job.result(), qstc); rho = qstf.fit(method='lstsq')
    Cexp[j] = coh_l1(rho)

matplotlib.rcParams.update({'font.size':12}); plt.figure(figsize = (6,4), dpi = 100)
plt.plot(lb, Cteo, label = r'$C_{teo}$'); plt.plot(lb, Csim, 'o', label = r'$C_{sim}$')
plt.plot(lb, Cexp, '*', label = r'$C_{exp}$')
plt.legend(); plt.xlabel(r'$\lambda$'); plt.show()

2960/8

750/15

