%run init.ipynb

c00,c01,c10,c11 = symbols('c_{00} c_{01} c_{10} c_{11}')
psiAB = Matrix([[c00],[c01],[c10],[c11]]); psiAB

rhoAB = psiAB*conjugate(psiAB.T); rhoAB

rhoA = ptraceB(2,2,rhoAB); rhoB = ptraceA(2,2,rhoAB)
rhoA, rhoB

pauli(3)*pauli(1)*pauli(3), pauli(3)*pauli(2)*pauli(3), comm(pauli(3),pauli(1)), comm(pauli(3),pauli(2))

%run init.ipynb

def pTraceL_num(dl, dr, rhoLR):
    rhoR = np.zeros((dr, dr), dtype=complex)
    for j in range(0, dr):
        for k in range(j, dr):
            for l in range(0, dl):
                rhoR[j,k] += rhoLR[l*dr+j,l*dr+k]
            if j != k:
                rhoR[k,j] = np.conj(rhoR[j,k])
    return rhoR
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

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 += np.abs(rho[j,k])
    return 2*C

from numpy import linalg as LA
def Econc_(rho):
    s2 = np.array([[0,-1j],[1j,0]])
    R = rho@np.kron(s2,s2)@np.matrix.conjugate(rho)@np.kron(s2,s2)
    w, v = LA.eig(R)
    evm = max(abs(w[0]), abs(w[1]), abs(w[2]), abs(w[3]))
    Ec = 2*math.sqrt(abs(evm)) - math.sqrt(abs(w[0])) - math.sqrt(abs(w[1]))\
    - math.sqrt(abs(w[2])) - math.sqrt(abs(w[3]))
    if Ec < 0.0:
        Ec = 0.0
    return Ec

import scipy.linalg.lapack as lapak
def Econc(rho):
    s2 = np.array([[0,-1j],[1j,0]])
    R = rho@np.kron(s2,s2)@np.matrix.conjugate(rho)@np.kron(s2,s2)
    ev = lapak.zheevd(R)
    evm = max(abs(ev[0][0]), abs(ev[0][1]), abs(ev[0][2]), abs(ev[0][3]))
    Ec = 2*math.sqrt(abs(evm)) - math.sqrt(abs(ev[0][0])) - math.sqrt(abs(ev[0][1]))\
    - math.sqrt(abs(ev[0][2])) - math.sqrt(abs(ev[0][3]))
    if Ec < 0.0:
        Ec = 0.0
    return Ec

def concurrence_psi(rho):
    rho_r = pTraceL_num(2,2,rho)
    return math.sqrt(2*abs(1-np.matrix.trace(rho_r@rho_r)))

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

import qiskit
from qiskit import *
nshots = 8192
IBMQ.load_account()
#provider= qiskit.IBMQ.get_provider(hub='ibm-q-research-2',group='federal-uni-sant-1',project='main')
provider = qiskit.IBMQ.get_provider(hub='ibm-q', group='open', project='main')
device = provider.get_backend('ibmq_manila')
simulator = Aer.get_backend('qasm_simulator')
from qiskit.tools.monitor import job_monitor
from qiskit.ignis.verification.tomography import state_tomography_circuits, StateTomographyFitter
from qiskit.ignis.mitigation.measurement import complete_meas_cal, CompleteMeasFitter

q_exp = np.arange(0,1.1,0.1); N = len(q_exp)
th_xi = np.zeros(N); th_xi = 2*np.arccos(np.sqrt(1-q_exp))
th_eta = np.zeros(N); th_eta = 2*np.arccos(np.sqrt(q_exp))
E_sim_10 = np.zeros(N); C_sim_0 = np.zeros(N); P_sim_0 = np.zeros(N)
jobs_ids = []

for j in range(0,N):
    qr = QuantumRegister(5); qc = QuantumCircuit(qr)
    # prepara xi
    qc.u(th_xi[j],0,0,qr[1])
    qc.cx(qr[1],qr[2])
    qc.h([qr[1],qr[2]])
    # prepara eta
    qc.u(th_eta[j],0,0,qr[3])
    qc.cx(qr[3],qr[4])
    qc.h([qr[3],qr[4]])
    # tomografia
    qstc = state_tomography_circuits(qc, [qr[1],qr[2]])
    
    job_sim = qiskit.execute(qstc, backend = simulator, shots=nshots)
    qstf_sim = StateTomographyFitter(job_sim.result(), qstc)
    rho_10_sim = qstf_sim.fit(method='lstsq')
    E_sim_10[j] = concurrence_psi(rho_10_sim)
    rho_0_sim = pTraceL_num(2, 2, rho_10_sim)
    C_sim_0[j] = coh_l1(rho_0_sim)
    P_sim_0[j] = predict_jb(rho_0_sim)

qc.draw(output='mpl')

q_exp = np.arange(0,1.1,0.1); N = len(q_exp)
th_xi = np.zeros(N); th_xi = 2*np.arccos(np.sqrt(1-q_exp))
th_eta = np.zeros(N); th_eta = 2*np.arccos(np.sqrt(q_exp))
jobs_ids = []

for j in range(0,N):
    qr = QuantumRegister(5); qc = QuantumCircuit(qr)
    # prepara xi
    qc.u(th_xi[j],0,0,qr[1])
    qc.cx(qr[1],qr[2])
    qc.h([qr[1],qr[2]])
    # prepara eta
    qc.u(th_eta[j],0,0,qr[3])
    qc.cx(qr[3],qr[4])
    qc.h([qr[3],qr[4]])
    # tomografia
    qstc = state_tomography_circuits(qc, [qr[1],qr[2]])
    
    job_exp = qiskit.execute(qstc, backend = device, shots=nshots)
    jobs_ids.append(job_exp.job_id())
    print(job_exp.job_id())
    job_monitor(job_exp)

print(jobs_ids)

f = open("jobs_ids_ES_pre_bbm.txt", "w")
f.write(str(jobs_ids))
f.close()

f = open("jobs_ids_ES_pre_bbm.txt","r")
#jobs_ids_list = f.read()
list_ids = f.read().replace("'","").replace(" ","").replace("[","").replace("]","").split(",")
f.close()

list_ids[0]

print(list_ids)

print(jobs_ids_list)

len_ids = len('6368049443e1f0708afdaf73')
print(len_ids)

list_ids = []
for j in range(0,N):
    s = 0
    if j > 0:
        s = 2
    print(jobs_ids_list[(j*(len_ids+2+s)+2):(j+1)*(len_ids+2)+j*s])
    list_ids.append(jobs_ids_list[(j*(len_ids+2+s)+2):(j+1)*(len_ids+2)+j*s])

list_ids[0]

q_exp = np.arange(0,1.1,0.1); N = len(q_exp)
E_exp_10 = np.zeros(N); C_exp_0 = np.zeros(N); P_exp_0 = np.zeros(N)
for j in range(0,N):
    job = device.retrieve_job(list_ids[j])
    qstf_exp = StateTomographyFitter(job.result(), qstc)
    rho_10_exp = qstf_exp.fit(method='lstsq'); 
    E_exp_10[j] = Econc_(rho_10_exp)
    print(E_exp_10[j])
    rho_0_exp = pTraceL_num(2, 2, rho_10_exp)
    C_exp_0[j] = coh_l1(rho_0_exp)
    P_exp[j] = predict_jb(rho_0_exp)

import matplotlib
matplotlib.rcParams.update({'font.size':11}); plt.figure(figsize = (7,5), dpi = 100)
q = np.arange(0,1.01,0.01)
Ei2 = 4*q*(1-q)
Ci2 = (2*q-1)**2
Pi2 = np.zeros(len(q))
Pm = 1-2*q*(1-q)
plt.plot(q, Ei2, '-', label=r'$E_{conc}(\xi_{AC})^2$')
plt.plot(q, Ci2, '--', label=r'$C_{l_1}(\xi_C)^2$')
plt.plot(q, Pi2, '-.', label=r'$P(\xi_C)^2$')
plt.plot(q, Ei2+Ci2, ':', label=r'$E_{conc}(\xi_{AC})^{2}+C_{l_{1}}(\xi_{C})^{2}+P(\xi_C)^2$')
plt.plot(q_exp, E_sim_10**2, 'X', label=r'$E_{conc}(\xi_{AC})^2_{sim}$')
plt.plot(q_exp, C_sim_0**2, 'h', label=r'$C_{l_1}(\xi_C)^2_{sim}$')
plt.plot(q_exp, P_sim_0**2, '>', label=r'$P(\xi_C)^2_{sim}$')
plt.plot(q_exp, E_sim_10**2+C_sim_0**2+P_sim_0**2, '*', label=r'$E_{conc}(\xi_{AC})^{2}_{sim}+C_{l_{1}}(\xi_{C})^{2}_{sim}+P(\xi_C)^2_{sim}$')
plt.plot(q_exp, E_exp_10**2, 's', label=r'$E_{conc}(\xi_{AC})^2_{exp}$')
plt.plot(q_exp, C_exp_0**2, 'o', label=r'$C_{l_1}(\xi_C)^2_{exp}$')
plt.plot(q_exp, P_exp_0**2, 'v', label=r'$P(\xi_C)^2_{exp}$')
plt.plot(q_exp, E_exp_10**2+C_exp_0**2+P_exp_0**2, 'd', label=r'$E_{conc}(\xi_{AC})^{2}_{exp}+C_{l_{1}}(\xi_{C})^{2}_{exp}+P(\xi_C)^2_{exp}$')
plt.legend(bbox_to_anchor=(1.75, 1.0), loc='upper right')
plt.xlabel(r'$q$')  
plt.show()

# error mitigation
qr = QuantumRegister(5); qubit_list = [1,2]
meas_calibs, state_labels = complete_meas_cal(qubit_list = qubit_list, qr = qr)
job = qiskit.execute(meas_calibs, backend = device, shots = nshots)
print(job.job_id())
job_monitor(job)

job = device.retrieve_job('636d895d04e46a1feb6227a9')
meas_fitter = CompleteMeasFitter(job.result(), state_labels)

q_exp = np.arange(0,1.1,0.1); N = len(q_exp)
for j in range(0,N):
    job = device.retrieve_job(list_ids[j])
    mitigated_results = meas_fitter.filter.apply(job.result())
    qstf_exp = StateTomographyFitter(mitigated_results, qstc)
    rho_10_exp = qstf_exp.fit(method='lstsq'); 
    E_exp_10[j] = Econc_(rho_10_exp)
    print(E_exp_10[j])
    rho_0_exp = pTraceL_num(2, 2, rho_10_exp)
    C_exp_0[j] = coh_l1(rho_0_exp)
    P_exp[j] = predict_jb(rho_0_exp)

import matplotlib
matplotlib.rcParams.update({'font.size':11}); plt.figure(figsize = (7,5), dpi = 100)
q = np.arange(0,1.01,0.01)
Ei2 = 4*q*(1-q)
Ci2 = (2*q-1)**2
Pi2 = np.zeros(len(q))
Pm = 1-2*q*(1-q)
plt.plot(q, Ei2, '-', label=r'$E_{conc}(\xi_{AC})^2$')
plt.plot(q, Ci2, '--', label=r'$C_{l_1}(\xi_C)^2$')
plt.plot(q, Pi2, '-.', label=r'$P(\xi_C)^2$')
plt.plot(q, Ei2+Ci2, ':', label=r'$E_{conc}(\xi_{AC})^{2}+C_{l_{1}}(\xi_{C})^{2}+P(\xi_C)^2$')
plt.plot(q_exp, E_sim_10**2, 'X', label=r'$E_{conc}(\xi_{AC})^2_{sim}$')
plt.plot(q_exp, C_sim_0**2, 'h', label=r'$C_{l_1}(\xi_C)^2_{sim}$')
plt.plot(q_exp, P_sim_0**2, '>', label=r'$P(\xi_C)^2_{sim}$')
plt.plot(q_exp, E_sim_10**2+C_sim_0**2+P_sim_0**2, '*', label=r'$E_{conc}(\xi_{AC})^{2}_{sim}+C_{l_{1}}(\xi_{C})^{2}_{sim}+P(\xi_C)^2_{sim}$')
plt.plot(q_exp, E_exp_10**2, '^', label=r'$E_{conc}(\xi_{AC})^2_{exp}$')
plt.plot(q_exp, C_exp_0**2, 'o', label=r'$C_{l_1}(\xi_C)^2_{exp}$')
plt.plot(q_exp, P_exp_0**2, '.', label=r'$P(\xi_C)^2_{exp}$')
plt.plot(q_exp, E_exp_10**2+C_exp_0**2+P_exp_0**2, 'd', label=r'$E_{conc}(\xi_{AC})^{2}_{exp}+C_{l_{1}}(\xi_{C})^{2}_{exp}+P(\xi_C)^2_{exp}$')
plt.legend(bbox_to_anchor=(1.75, 1.0), loc='upper right')
plt.xlabel(r'$q$')
plt.savefig('fig_Eswap_pre_bbm.pdf')  
plt.show()

qc.draw(output='mpl')

def pTrace_3_num(rho_abc, da, db, dc):
    rho_ac = np.zeros(da*dc*da*dc, dtype=complex).reshape(da*dc,da*dc)
    for j in range(0,da):
        for l in range(0,dc):
            cj = j*dc+l
            ccj = j*db*dc+l
            for m in range(0,da):
                for o in range(0,dc):
                    ck = m*dc+o
                    cck = m*db*dc+o
                    for k in range(0,db):
                        rho_ac[cj,ck] = rho_ac[cj,ck] + rho_abc[ccj+k*dc,cck+k*dc]
    return rho_ac

import numpy as np
import math

# projetores para a pos-seleção
I2 = np.array([[1,0],[0,1]], dtype=complex)
P00 = np.kron(np.kron(I2,np.array([[1,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],dtype=complex)), I2) # PhiP
P01 = np.kron(np.kron(I2,np.array([[0,0,0,0],[0,1,0,0],[0,0,0,0],[0,0,0,0]],dtype=complex)), I2) # PsiP
P10 = np.kron(np.kron(I2,np.array([[0,0,0,0],[0,0,0,0],[0,0,1,0],[0,0,0,0]],dtype=complex)), I2) # PhiM
P11 = np.kron(np.kron(I2,np.array([[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,1]],dtype=complex)), I2) # PsiM
#P00+P01+P10+P11

q_exp = np.arange(0,1.1,0.1); N = len(q_exp)
q_exp[0] = 0.03; q_exp[N-1] = 1-q_exp[0]
th_xi = np.zeros(N); th_xi = 2*np.arccos(np.sqrt(1-q_exp))
th_eta = np.zeros(N); th_eta = 2*np.arccos(np.sqrt(q_exp))
E_phiP_04_sim = np.zeros(N); E_psiM_04_sim = np.zeros(N)
E_psiP_04_sim = np.zeros(N); E_phiM_04_sim = np.zeros(N)
C2_sim = np.zeros(N)

for j in range(0,N):
    qr = QuantumRegister(5); qc = QuantumCircuit(qr)
    # prepara xi
    qc.u(th_xi[j],0,0,qr[0])
    qc.cx(qr[0],qr[1])
    qc.h([qr[0],qr[1]])
    # prepara eta
    qc.u(th_eta[j],0,0,qr[3])
    qc.cx(qr[3],qr[4])
    qc.h([qr[3],qr[4]])
    # BBM
    qc.cx(qr[1],qr[3])
    qc.h(qr[1])
    # circuito para tomografia
    qstc = state_tomography_circuits(qc, [qr[4],qr[3],qr[1],qr[0]])
    
    # Simulacao
    job_sim = qiskit.execute(qstc, backend = simulator, shots=nshots)
    qstf_sim = StateTomographyFitter(job_sim.result(), qstc)
    rho_4310_sim = qstf_sim.fit(method='lstsq')
    # pos-selecao
    rho_PhiP = P00@rho_4310_sim@P00
    pr_phiP = np.real(np.trace(rho_PhiP))
    rho_PhiP = (1/pr_phiP)*rho_PhiP
    rho_PsiP = P01@rho_4310_sim@P01
    pr_psiP = np.real(np.trace(rho_PsiP))
    rho_PsiP = (1/pr_psiP)*rho_PsiP
    rho_PhiM = P10@rho_4310_sim@P10
    pr_phiM = np.real(np.trace(rho_PhiM))
    rho_PhiM = (1/pr_phiM)*rho_PhiM
    rho_PsiM = P11@rho_4310_sim@P11
    pr_psiM = np.real(np.trace(rho_PsiM))
    rho_PsiM = (1/pr_psiM)*rho_PsiM
    # estados reduzidos
    rho_PhiP_04 = pTrace_3_num(rho_PhiP,2,4,2)
    rho_PsiP_04 = pTrace_3_num(rho_PsiP,2,4,2)
    rho_PhiM_04 = pTrace_3_num(rho_PhiM,2,4,2)
    rho_PsiM_04 = pTrace_3_num(rho_PsiM,2,4,2)
    # emaranhamentos e probabilidades
    E_phiP_04_sim[j] = Econc_(rho_PhiP_04)
    E_phiM_04_sim[j] = Econc_(rho_PhiM_04)
    E_psiP_04_sim[j] = Econc_(rho_PsiP_04)
    E_psiM_04_sim[j] = Econc_(rho_PsiM_04)
    C2_sim[j] = 1 - 2*(pr_phiP+pr_psiP)

provider = qiskit.IBMQ.get_provider(hub='ibm-q', group='open', project='main')
device = provider.get_backend('ibmq_manila')

q_exp = np.arange(0,1.1,0.1); N = len(q_exp)
q_exp[0] = 0.03; q_exp[N-1] = 1-q_exp[0]
th_xi = np.zeros(N); th_xi = 2*np.arccos(np.sqrt(1-q_exp))
th_eta = np.zeros(N); th_eta = 2*np.arccos(np.sqrt(q_exp))
#jobs_ids2 = []

for j in range(N-1,N):
    qr = QuantumRegister(5); qc = QuantumCircuit(qr)
    # prepara xi
    qc.u(th_xi[j],0,0,qr[1])
    qc.cx(qr[1],qr[2])
    qc.h([qr[1],qr[2]])
    # prepara eta
    qc.u(th_eta[j],0,0,qr[3])
    qc.cx(qr[3],qr[4])
    qc.h([qr[3],qr[4]])
    # BBM
    qc.cx(qr[2],qr[3])
    qc.h(qr[2])
    # circuito para tomografia
    qstc = state_tomography_circuits(qc, [qr[4],qr[3],qr[2],qr[1]])
    
    # Experimento
    job_exp = qiskit.execute(qstc, backend = device, shots=nshots)
    #jobs_ids2.append(job_exp.job_id())
    print(job_exp.job_id())
    job_monitor(job_exp)

jobs_ids2.append('636e2c8910cc1925ad755b77')

qc.draw(output='mpl')

f = open("jobs_ids_ES_post_bbm.txt", "w")
f.write(str(jobs_ids2))
f.close()

f = open("jobs_ids_ES_post_bbm.txt","r")
list_ids2 = f.read().replace("'","").replace(" ","").replace("[","").replace("]","").split(",")
f.close()
print(list_ids2)

q_exp = np.arange(0,1.1,0.1); N = len(q_exp)
q_exp[0] = 0.03; q_exp[N-1] = 1-q_exp[0]
E_phiP_04_exp = np.zeros(N); E_psiM_04_exp = np.zeros(N)
E_psiP_04_exp = np.zeros(N); E_phiM_04_exp = np.zeros(N)
C2_exp = np.zeros(N)
for j in range(0,N):
    job = device.retrieve_job(list_ids2[j])
    qstf_exp = StateTomographyFitter(job.result(), qstc)
    rho_4310_exp = qstf_exp.fit(method='lstsq')
    # pos-selecao
    rho_PhiP = P00@rho_4310_exp@P00
    pr_phiP = np.real(np.trace(rho_PhiP))
    rho_PhiP = (1/pr_phiP)*rho_PhiP
    rho_PsiP = P01@rho_4310_exp@P01
    pr_psiP = np.real(np.trace(rho_PsiP))
    rho_PsiP = (1/pr_psiP)*rho_PsiP
    rho_PhiM = P10@rho_4310_exp@P10
    pr_phiM = np.real(np.trace(rho_PhiM))
    rho_PhiM = (1/pr_phiM)*rho_PhiM
    rho_PsiM = P11@rho_4310_exp@P11
    pr_psiM = np.real(np.trace(rho_PsiM))
    rho_PsiM = (1/pr_psiM)*rho_PsiM
    # estados reduzidos
    rho_PhiP_04 = pTrace_3_num(rho_PhiP,2,4,2)
    rho_PsiP_04 = pTrace_3_num(rho_PsiP,2,4,2)
    rho_PhiM_04 = pTrace_3_num(rho_PhiM,2,4,2)
    rho_PsiM_04 = pTrace_3_num(rho_PsiM,2,4,2)
    # emaranhamentos e probabilidades
    E_phiP_04_exp[j] = Econc_(rho_PhiP_04)
    E_phiM_04_exp[j] = Econc_(rho_PhiM_04)
    E_psiP_04_exp[j] = Econc_(rho_PsiP_04)
    E_psiM_04_exp[j] = Econc_(rho_PsiM_04)
    C2_exp[j] = 1 - 2*(pr_phiP+pr_psiP)
    print(C2_exp[j])

import matplotlib
from matplotlib import pyplot as plt
matplotlib.rcParams.update({'font.size':11}); plt.figure(figsize = (6,4), dpi = 100)
q = np.arange(0,1.01,0.01)
Ep = np.ones(len(q)); Em = (2*q*(1-q))/(1-2*q+2*q**2); C2 = 1-4*q*(1-q)
plt.plot(q, Ep, '-', label=r'$E(\hat{\phi}_+)=E(\hat{\psi}_+)$')
plt.plot(q, Em, '--', label=r'$E(\hat{\phi}_-)=E(\hat{\psi}_-)$')
plt.plot(q, C2, '-.', label=r'$1-2(Pr(\hat{\phi}_+)+Pr(\hat{\psi}_+))$')
plt.plot(q_exp, E_phiP_04_sim, 'o', label=r'$E(\hat{\phi}_+)_{sim}$')
plt.plot(q_exp, E_phiM_04_sim, '*', label=r'$E(\hat{\phi}_-)_{sim}$')
#plt.plot(q_exp, E_psiP_04_sim, 'X', label=r'$E(\hat{\psi}_+)_{sim}$')
#plt.plot(q_exp, E_psiM_04_sim, '.', label=r'$E(\hat{\psi}_-)_{sim}$')
plt.plot(q_exp, C2_sim, 'X', label=r'$1-2(Pr(\hat{\phi}_+)_{sim}+Pr(\hat{\psi}_+)_{sim})$')
plt.plot(q_exp, E_phiP_04_exp, '<', label=r'$E(\hat{\phi}_+)_{exp}$')
plt.plot(q_exp, E_phiM_04_exp, '>', label=r'$E(\hat{\phi}_-)_{exp}$')
#plt.plot(q_exp, E_psiP_04_exp, '^', label=r'$E(\hat{\psi}_+)_{exp}$')
#plt.plot(q_exp, E_psiM_04_exp, '.', label=r'$E(\hat{\psi}_-)_{exp}$')
plt.plot(q_exp, C2_exp, '^', label=r'$1-2(Pr(\hat{\phi}_+)_{exp}+Pr(\hat{\psi}_+)_{exp})$')
plt.xlabel(r'$q$')
plt.legend(bbox_to_anchor=(1.8, 1.0), loc='upper right')
plt.show()

# error mitigation
qr = QuantumRegister(5)#; qc = QuantumCircuit(qr)
qubit_list = [1,2,3,4]
meas_calibs, state_labels = complete_meas_cal(qubit_list = qubit_list, qr = qr)
job = qiskit.execute(meas_calibs, backend = device, shots = nshots)
print(job.job_id())
job_monitor(job)

job = device.retrieve_job('636e3407fb0a57272187586a')
meas_fitter = CompleteMeasFitter(job.result(), state_labels)

q_exp = np.arange(0,1.1,0.1); N = len(q_exp)
q_exp[0] = 0.03; q_exp[N-1] = 1-q_exp[0]
E_phiP_04_exp = np.zeros(N); E_psiM_04_exp = np.zeros(N)
E_psiP_04_exp = np.zeros(N); E_phiM_04_exp = np.zeros(N)
C2_exp = np.zeros(N)
for j in range(0,N):
    job = device.retrieve_job(list_ids2[j])
    mitigated_results = meas_fitter.filter.apply(job.result()) # error mitigation
    qstf_exp = StateTomographyFitter(mitigated_results, qstc) # error mitigation
    rho_4310_exp = qstf_exp.fit(method='lstsq')
    # pos-selecao
    rho_PhiP = P00@rho_4310_exp@P00
    pr_phiP = np.real(np.trace(rho_PhiP))
    rho_PhiP = (1/pr_phiP)*rho_PhiP
    rho_PsiP = P01@rho_4310_exp@P01
    pr_psiP = np.real(np.trace(rho_PsiP))
    rho_PsiP = (1/pr_psiP)*rho_PsiP
    rho_PhiM = P10@rho_4310_exp@P10
    pr_phiM = np.real(np.trace(rho_PhiM))
    rho_PhiM = (1/pr_phiM)*rho_PhiM
    rho_PsiM = P11@rho_4310_exp@P11
    pr_psiM = np.real(np.trace(rho_PsiM))
    rho_PsiM = (1/pr_psiM)*rho_PsiM
    # estados reduzidos
    rho_PhiP_04 = pTrace_3_num(rho_PhiP,2,4,2)
    rho_PsiP_04 = pTrace_3_num(rho_PsiP,2,4,2)
    rho_PhiM_04 = pTrace_3_num(rho_PhiM,2,4,2)
    rho_PsiM_04 = pTrace_3_num(rho_PsiM,2,4,2)
    # emaranhamentos e probabilidades
    E_phiP_04_exp[j] = Econc_(rho_PhiP_04)
    E_phiM_04_exp[j] = Econc_(rho_PhiM_04)
    E_psiP_04_exp[j] = Econc_(rho_PsiP_04)
    E_psiM_04_exp[j] = Econc_(rho_PsiM_04)
    C2_exp[j] = 1 - 2*(pr_phiP+pr_psiP)
    print(C2_exp[j])

import matplotlib
from matplotlib import pyplot as plt
matplotlib.rcParams.update({'font.size':11}); plt.figure(figsize = (6,4), dpi = 100)
q = np.arange(0,1.01,0.01)
Ep = np.ones(len(q)); Em = (2*q*(1-q))/(1-2*q+2*q**2); C2 = 1-4*q*(1-q)
plt.plot(q, Ep, '-', label=r'$E(\hat{\phi}_+)=E(\hat{\psi}_+)$')
plt.plot(q, Em, '--', label=r'$E(\hat{\phi}_-)=E(\hat{\psi}_-)$')
plt.plot(q, C2, '-.', label=r'$1-2(Pr(\hat{\phi}_+)+Pr(\hat{\psi}_+))$')
plt.plot(q_exp, E_phiP_04_sim, 'o', label=r'$E(\hat{\phi}_+)_{sim}$')
plt.plot(q_exp, E_phiM_04_sim, '*', label=r'$E(\hat{\phi}_-)_{sim}$')
#plt.plot(q_exp, E_psiP_04_sim, 'X', label=r'$E(\hat{\psi}_+)_{sim}$')
#plt.plot(q_exp, E_psiM_04_sim, '.', label=r'$E(\hat{\psi}_-)_{sim}$')
plt.plot(q_exp, C2_sim, 'X', label=r'$1-2(Pr(\hat{\phi}_+)_{sim}+Pr(\hat{\psi}_+)_{sim})$')
plt.plot(q_exp, E_phiP_04_exp, '<', label=r'$E(\hat{\phi}_+)_{exp}$')
plt.plot(q_exp, E_phiM_04_exp, '>', label=r'$E(\hat{\phi}_-)_{exp}$')
#plt.plot(q_exp, E_psiP_04_exp, '^', label=r'$E(\hat{\psi}_+)_{exp}$')
#plt.plot(q_exp, E_psiM_04_exp, '.', label=r'$E(\hat{\psi}_-)_{exp}$')
plt.plot(q_exp, C2_exp, '^', label=r'$1-2(Pr(\hat{\phi}_+)_{exp}+Pr(\hat{\psi}_+)_{exp})$')
plt.xlabel(r'$q$')
plt.legend(bbox_to_anchor=(1.8, 1.0), loc='upper right')
plt.savefig('fig_Eswap_post_bbm.pdf')
plt.show()

C = np.arange(0,1.1,0.1); 
E_PhiP = (1-C**2)/(1+C**2)
E_PhiM = np.ones(len(C))
plt.plot(C, E_PhiP, '*'); plt.plot(C, E_PhiM, '.')
plt.show()



from qiskit import *

qr = QuantumRegister(4); 



%run init.ipynb

p,q = symbols('p q')
ket_xi = Matrix([[sqrt(p)],[sqrt((1-p)/2)],[sqrt((1-p)/2)],[0]])
ket_eta = Matrix([[sqrt(q)],[sqrt((1-q)/2)],[sqrt((1-q)/2)],[0]])
ket_xi, ket_eta

rho_xi = ket_xi*ket_xi.T
rho_eta = ket_eta*ket_eta.T
rho_xi, rho_eta

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

rho_xi_A = ptraceB(2,2,rho_xi); rho_xi_C = ptraceA(2,2,rho_xi)
rho_xi_A, rho_xi_C

rho_eta_Cp = ptraceB(2,2,rho_eta); rho_eta_B = ptraceA(2,2,rho_eta)
rho_eta_Cp, rho_eta_B

c00,c01,c10,c11 = symbols('c_{00} c_{01} c_{10} c_{11}')
d00,d01,d10,d11 = symbols('d_{00} d_{01} d_{10} d_{11}')
def two_qb_basis():
    zz = Matrix([[1],[0],[0],[0]])
    zu = Matrix([[0],[1],[0],[0]])
    uz = Matrix([[0],[0],[1],[0]])
    uu = Matrix([[0],[0],[0],[1]])
    return zz,zu,uz,uu
zz,zu,uz,uu = two_qb_basis(); zz,zu,uz,uu

def psi_p(c00,c01,c10,c11,d00,d01,d10,d11):
    zz,zu,uz,uu = two_qb_basis()
    psi = (c00*d00+c01*d10)*zz + (c00*d01+c01*d11)*zu + (c10*d00+c11*d10)*uz + (c10*d01+c11*d11)*uu
    return psi/sqrt(2)

def psi_m(c00,c01,c10,c11,d00,d01,d10,d11):
    zz,zu,uz,uu = two_qb_basis()
    psi = (c00*d00-c01*d10)*zz + (c00*d01-c01*d11)*zu + (c10*d00-c11*d10)*uz + (c10*d01-c11*d11)*uu
    return psi/sqrt(2)

def phi_p(c00,c01,c10,c11,d00,d01,d10,d11):
    zz,zu,uz,uu = two_qb_basis()
    psi = (c00*d10+c01*d00)*zz + (c00*d11+c01*d01)*zu + (c10*d10+c11*d00)*uz + (c10*d11+c11*d01)*uu
    return psi/sqrt(2)

def phi_m(c00,c01,c10,c11,d00,d01,d10,d11):
    zz,zu,uz,uu = two_qb_basis()
    psi = (c00*d10-c01*d00)*zz + (c00*d11-c01*d01)*zu + (c10*d10-c11*d00)*uz + (c10*d11-c11*d01)*uu
    return psi/sqrt(2)

c00 = sqrt(p); c01 = sqrt((1-p)/2); c10 = sqrt((1-p)/2); c11 = 0
d00 = sqrt(q); d01 = sqrt((1-q)/2); d10 = sqrt((1-q)/2); d11 = 0
psip = psi_p(c00,c01,c10,c11,d00,d01,d10,d11); 
psim = psi_m(c00,c01,c10,c11,d00,d01,d10,d11)
phip = phi_p(c00,c01,c10,c11,d00,d01,d10,d11); 
phim = phi_m(c00,c01,c10,c11,d00,d01,d10,d11)
simplify(psip), simplify(psim), simplify(phip), simplify(phim)

psip_norm2 = psip.T*psip; simplify(psip_norm2)

psim_norm2 = psim.T*psim; simplify(psim_norm2)

phip_norm2 = phip.T*phip; simplify(phip_norm2)

phim_norm2 = phim.T*phim; simplify(phim_norm2)