import math; import numpy as np; phmax = 4*math.pi; npt = 100; dth = phmax/npt; ph = np.arange(0,phmax+dth,dth)
from matplotlib import pyplot as plt; PrD0 = 0.5*(1+np.cos(ph)); plt.plot(ph,PrD0, label=r'$Pr(D_0)$')
PrD0 = 0.5*(1-np.cos(ph)); plt.plot(ph,PrD0, label=r'$Pr(D_1)$'); plt.legend(); plt.xlabel(r'$\phi$'); plt.show()

from sympy import Matrix, sqrt
S = Matrix([[1,0],[0,1j]]); H = (1/sqrt(2))*Matrix([[1,1],[1,-1]])
S*H*S

Y = Matrix([[0,-1j],[1j,0]]); Z = Matrix([[1,0],[0,-1]]); Y*Z

from qiskit import QuantumCircuit

def qc_df(): # divisor de feixes 50:50
    qc = QuantumCircuit(1, name='DF')
    qc.s(0)
    qc.h(0)
    qc.s(0)
    return qc
qcdf = qc_df(); qcdf.draw('mpl')

def qc_espelhos():
    qc = QuantumCircuit(1, name='E')
    qc.z(0)
    qc.y(0)
    return qc
qce = qc_espelhos(); qce.draw('mpl')

def qc_imz(ph): # retorna o circuito quântico para o IMZ (sem as medidas)
    qc = QuantumCircuit(1, name='IMZ')
    qcdf = qc_df(); qc.append(qcdf, [0])
    qce = qc_espelhos(); qc.append(qce, [0])
    qc.p(ph, [0])
    qcdf = qc_df(); qc.append(qcdf, [0])
    return qc
qcimz = qc_imz(0); qcimz.draw('mpl')

# Simulação clássica
from qiskit import Aer, execute
simulator = Aer.get_backend('qasm_simulator')
nshots = 2**13
import numpy as np; import math
phmax = 2*math.pi
npe = 6
dth = phmax/npe
ph = np.arange(0,phmax+dth,dth)
PD0sim = np.zeros(len(ph))
for j in range(0,len(ph)):
    qc = QuantumCircuit(1,1)
    qcimz = qc_imz(ph[j])
    qc.append(qcimz,[0])
    qc.measure(0,0)
    job = execute(qc, backend=simulator, shots=nshots)
    counts = job.result().get_counts()
    if '0' in counts:
        PD0sim[j] = counts['0']/nshots
print(PD0sim)

npt = 100
dtht = phmax/npt
pht = np.arange(0,phmax+dtht,dtht)
PD0teo = 0.5*(1+np.cos(pht))
from matplotlib import pyplot as plt
plt.plot(pht,PD0teo, label=r'$Pr(D_0)_{teo}$')
plt.plot(ph,PD0sim, 'o', label=r'$Pr(D_0)_{sim}$')
plt.legend(); plt.xlabel(r'$\phi$')
plt.show()

import qiskit
qiskit.IBMQ.save_account('a221013b3d88dd070a495178c3f732ac683a494e5af4e641383f16d361e9fa7fa22af7f09cad3f88988a3e6e5621c7fe951d8ca6ff0c0371422b211ec0ae9709', 
                  overwrite = True)
qiskit.IBMQ.load_account()
#provider = qiskit.IBMQ.get_provider(hub='ibm-q', group='open', project='main')
provider = qiskit.IBMQ.get_provider(hub='ibm-q-research-2', group='federal-uni-sant-1', project='main')
device = provider.get_backend('ibm_nairobi')
nshots = 2**13

import numpy as np; import math
phmax = 2*math.pi
npe = 6
dth = phmax/npe
ph = np.arange(0,phmax+dth,dth)
job_ids = []
for j in range(0,len(ph)):
    qc = QuantumCircuit(1,1)
    qcimz = qc_imz(ph[j])
    qc.append(qcimz,[0])
    qc.measure(0,0)
    job = execute(qc, backend=device, shots=nshots)
    jobid = job.job_id()
    job_ids.append(jobid)
    print(jobid)
f = open("jobs_ids_IMZ.txt", "w")
f.write(str(job_ids))
f.close()

device = provider.get_backend('ibm_nairobi')
f = open("jobs_ids_IMZ.txt","r")
list_ids_IMZ = f.read().replace("'","").replace(" ","").replace("[","").replace("]","").split(",")
f.close()
PD0exp = np.zeros(len(ph))
for j in range(0,len(ph)):
    job = device.retrieve_job(list_ids_IMZ[j])
    counts = job.result().get_counts()
    if '0' in counts:
        PD0exp[j] = counts['0']/nshots
print(PD0exp)

npt = 100; dtht = phmax/npt; pht = np.arange(0,phmax+dtht,dtht)
PD0teo = 0.5*(1+np.cos(pht))
from matplotlib import pyplot as plt
plt.plot(pht,PD0teo, label=r'$Pr(D_0)_{teo}$')
plt.plot(ph,PD0sim, 'o', label=r'$Pr(D_0)_{sim}$')
plt.plot(ph,PD0exp, '*', label=r'$Pr(D_0)_{exp}$')
plt.legend()
plt.xlabel(r'$\phi$')
plt.show()

Vsim = (np.max(PD0sim)-np.min(PD0sim))/(np.max(PD0sim)+np.min(PD0sim)); print('Vsim = ',Vsim)
Vexp = (np.max(PD0exp)-np.min(PD0exp))/(np.max(PD0exp)+np.min(PD0exp)); print('Vexp = ',Vexp)

from sympy import Matrix, eye; P0 = Matrix([[1,0],[0,0]]); P1 = Matrix([[0,0],[0,1]]); X = Matrix([[0,1],[1,0]])
from sympy.physics.quantum import TensorProduct as tp; tp(eye(2),P0) + tp(X,P1)

def qc_bbo():
    from qiskit import QuantumCircuit
    qc = QuantumCircuit(2, name='BBO')
    qc.x(1)
    qc.h(0)
    qc.cx(0,1)
    return qc
qcbbo = qc_bbo(); qcbbo.draw('mpl')

def qc_dfp():
    from qiskit import QuantumCircuit
    qc = QuantumCircuit(2, name='DFP')
    qc.cz(0,1)
    qc.cy(0,1)
    return qc
qcdfp = qc_dfp()
qcdfp.draw('mpl')

def qc_pmo():
    from qiskit import QuantumCircuit
    qc = QuantumCircuit(2, name='PMO')
    qc.cx(1,0)
    return qc
qcpmo = qc_pmo()
qcpmo.draw('mpl')

def qc_imz_ic(ph): # não consideramos o caminho do fóton A
    from qiskit import QuantumCircuit
    qc = QuantumCircuit(3, name='IMZIC')
    qcbbo = qc_bbo(); qc.append(qcbbo, [0,1])
    qcdfp = qc_dfp(); qc.append(qcdfp, [1,2])
    qcpmo = qc_pmo(); qc.append(qcpmo, [1,2])
    qce = qc_espelhos(); qc.append(qce, [2])
    qc.p(ph, [2])
    qcdf = qc_df(); qc.append(qcdf, [2])
    return qc
qcimzic = qc_imz_ic(0); qcimzic.draw('mpl')

qcimzic.decompose().draw('mpl')

qcimzic.decompose().decompose().decompose().draw('mpl')

# Simulação clássica
from qiskit import Aer, execute
simulator = Aer.get_backend('qasm_simulator')
nshots = 2**13
import numpy as np; import math
phmax = 2*math.pi
npe = 5
dth = phmax/npe
ph = np.arange(0,phmax+dth,dth)
PD0sim = np.zeros(len(ph))
for j in range(0,len(ph)):
    qc = QuantumCircuit(3,2)
    qcimzic = qc_imz_ic(ph[j])
    qc.append(qcimzic,[0,1,2])
    qc.measure([0,2],[0,1])
    job = execute(qc, backend=simulator, shots=nshots)
    counts = job.result().get_counts()
    if '00' in counts:
        PD0sim[j] = counts['00']/nshots
    if '10' in counts:
        PD0sim[j] += counts['10']/nshots
print(PD0sim)

qc.draw('mpl')

import numpy as np; import math
phmax = 2*math.pi; npt = 100; dpht = phmax/npt
pht = np.arange(0,phmax+dtht,dtht)
PD0teo = 0.5*np.ones(len(pht))
from matplotlib import pyplot as plt
plt.plot(pht,PD0teo, label=r'$Pr(D_0)_{teo}$')
plt.plot(ph,PD0sim, 'o', label=r'$Pr(D_0)_{sim}$')
plt.ylim(0,1); plt.legend()
plt.xlabel(r'$\phi$')
plt.show()

import qiskit
qiskit.IBMQ.save_account('e20d3045c7c23a3be753606794714e9b47d4331ec46a41b2281e2d601bae27b2b3c1baa0b5cb54f3282a8b5248e552b14c9c556699010660fb94b11d00ff1073', 
                  overwrite = True)
qiskit.IBMQ.load_account()
#provider = qiskit.IBMQ.get_provider(hub='ibm-q', group='open', project='main')
provider = qiskit.IBMQ.get_provider(hub='ibm-q-research-2', group='federal-uni-sant-1', project='main')
device = provider.get_backend('ibm_nairobi')
nshots = 2**13

import numpy as np; import math
phmax=2*math.pi
npe=5
dth=phmax/npe
ph=np.arange(0,phmax+dth,dth)
job_ids=[]
for j in range(0,len(ph)):
    qc = QuantumCircuit(3,2)
    qcimzic = qc_imz_ic(ph[j])
    qc.append(qcimzic,[0,1,2])
    qc.measure([0,2],[0,1])
    job = execute(qc, backend=device, shots=nshots)
    jobid = job.job_id()
    job_ids.append(jobid)
f = open("jobs_ids_IMZIC.txt", "w")
f.write(str(job_ids))
f.close()

print(job_ids)

f = open("jobs_ids_IMZIC.txt","r")
list_ids_IMZIC = f.read().replace("'","").replace(" ","").replace("[","").replace("]","").split(",")
f.close()
phmax = 2*math.pi
npe=5
dth = phmax/npe
ph = np.arange(0,phmax+dth,dth)
PD0exp = np.zeros(len(ph))
device = provider.get_backend('ibm_nairobi')
for j in range(0,len(ph)):
    job = device.retrieve_job(list_ids_IMZIC[j])
    counts = job.result().get_counts()
    if '00' in counts:
        PD0exp[j] = counts['00']/nshots
    if '10' in counts:
        PD0exp[j] += counts['10']/nshots
print(PD0exp)

import numpy as np; import math
phmax = 2*math.pi; npt = 100; dpht = phmax/npt
pht = np.arange(0,phmax+dtht,dtht)
PD0teo = 0.5*np.ones(len(pht))
from matplotlib import pyplot as plt
plt.plot(pht,PD0teo, label=r'$Pr(D_0)_{teo}$')
plt.plot(ph,PD0sim, 'o', label=r'$Pr(D_0)_{sim}$')
plt.plot(ph,PD0exp, '*', label=r'$Pr(D_0)_{exp}$')
plt.ylim(0,1); plt.legend(); plt.xlabel(r'$\phi$'); plt.show()

Vsim = (np.max(PD0sim)-np.min(PD0sim))/(np.max(PD0sim)+np.min(PD0sim)); print('Vsim = ',Vsim)
Vexp = (np.max(PD0exp)-np.min(PD0exp))/(np.max(PD0exp)+np.min(PD0exp)); print('Vexp = ',Vexp)

import numpy as np; import math; phmax = 2*math.pi; npt = 100; dpht = phmax/npt; pht = np.arange(0,phmax+dpht,dpht)
PD0teoA0 = 0.5*(1+np.cos(pht)); PD0teoA1 = 0.5*(1-np.cos(pht))
from matplotlib import pyplot as plt
plt.plot(pht,PD0teoA0, label=r'$Pr(D_0|A=0)_{teo}$'); plt.plot(pht,PD0teoA1, label=r'$Pr(D_0|A=1)_{teo}$')
plt.xlim(0,2*math.pi); plt.xlim(-0.1,phmax+0.1); plt.ylim(-0.01,1.01); plt.legend(); plt.xlabel(r'$\phi$'); plt.show()

def qc_pqo():
    from qiskit import QuantumCircuit
    qc = QuantumCircuit(1, name='PQO')
    qc.h(0)
    qc.s(0)
    return qc
qcpqo = qc_pqo()
qcpqo.draw('mpl')

def qc_apagadorQ(ph):
    from qiskit import QuantumCircuit
    qc = QuantumCircuit(3)
    qcimzic = qc_imz_ic(ph)
    qc.append(qcimzic, [0,1,2])
    qcpqo = qc_pqo()
    qc.append(qcpqo, [0])
    return qc
qcaq = qc_apagadorQ(0)
qcaq.decompose().draw('mpl')

# Simulação clássica
from qiskit import Aer, execute; simulator = Aer.get_backend('qasm_simulator'); nshots = 2**13
import numpy as np; import math; phmax = 2*math.pi; npe = 6; dth = phmax/npe; ph = np.arange(0,phmax+dth,dth)
PA0sim = np.zeros(len(ph)); PA1sim = np.zeros(len(ph)); PD0A0sim = np.zeros(len(ph)); PD0A1sim = np.zeros(len(ph))
for j in range(0,len(ph)):
    qc = QuantumCircuit(3,2)
    qcaq = qc_apagadorQ(ph[j])
    qc.append(qcaq,[0,1,2])
    qc.measure([0,2],[0,1])
    job = execute(qc, backend=simulator, shots=nshots)
    counts = job.result().get_counts()
    if '00' in counts:
        PA0sim[j] = counts['00']/nshots
    if '01' in counts:
        PA0sim[j] += counts['01']/nshots
    if '10' in counts:
        PA1sim[j] = counts['10']/nshots
    if '11' in counts:
        PA1sim[j] += counts['11']/nshots
    if PA0sim[j] != 0 and '00' in counts:
        PD0A0sim[j] = (counts['00']/nshots)/PA0sim[j]
    if PA1sim[j] != 0 and '10' in counts:
        PD0A1sim[j] = (counts['10']/nshots)/PA1sim[j]
print(PD0A0sim, PD0A1sim)

import numpy as np; import math
phmax = 2*math.pi; npt = 100; dpht = phmax/npt; pht = np.arange(0,phmax+dpht,dpht)
PD0teoA0 = 0.5*(1+np.cos(pht)); PD0teoA1 = 0.5*(1-np.cos(pht))
from matplotlib import pyplot as plt
plt.plot(pht,PD0teoA0, label=r'$Pr(D_0|A=0)_{teo}$')
plt.plot(pht,PD0teoA1, label=r'$Pr(D_0|A=1)_{teo}$')
plt.plot(ph,PD0A0sim, 'o', label=r'$Pr(D_0|A=0)_{sim}$')
plt.plot(ph,PD0A1sim, 's', label=r'$Pr(D_0|A=1)_{sim}$')
plt.xlim(-0.1,2*math.pi+0.1); plt.ylim(-0.01,1.02); plt.legend(); plt.xlabel(r'$\phi$'); plt.show()

import qiskit
qiskit.IBMQ.save_account('dabd4674f7b3794f740813fec0c5478bfefdebb7d1dca9b7acd41e3f90ed2c8a8b80af26fb3236289fbfc70d0ea334829a836453ac46f96fee68abcb22495ef4', 
                  overwrite = True)
qiskit.IBMQ.load_account()
#provider = qiskit.IBMQ.get_provider(hub='ibm-q', group='open', project='main')
provider = qiskit.IBMQ.get_provider(hub='ibm-q-research-2', group='federal-uni-sant-1', project='main')
device = provider.get_backend('ibm_nairobi')
nshots = 2**13

import numpy as np; import math
phmax=2*math.pi
npe=6
dth=phmax/npe
ph=np.arange(0,phmax+dth,dth)
job_ids=[]
for j in range(0,len(ph)):
    qc = QuantumCircuit(3,2)
    qcaq = qc_apagadorQ(ph[j])
    qc.append(qcaq,[0,1,2])
    qc.measure([0,2],[0,1])
    job = execute(qc, backend=device, shots=nshots)
    jobid = job.job_id()
    job_ids.append(jobid)
    print(jobid)
f = open("jobs_ids_apagador.txt", "w")
f.write(str(job_ids))
f.close()

f = open("jobs_ids_apagador.txt","r")
list_ids_apagador = f.read().replace("'","").replace(" ","").replace("[","").replace("]","").split(",")
f.close()
phmax = 2*math.pi; npe=6; dth = phmax/npe; ph = np.arange(0,phmax+dth,dth)
PD0A0exp = np.zeros(len(ph)); PD0A1exp = np.zeros(len(ph)); PA0exp = np.zeros(len(ph)); PA1exp = np.zeros(len(ph))
device = provider.get_backend('ibm_nairobi')
for j in range(0,len(ph)):
    job = device.retrieve_job(list_ids_apagador[j])
    counts = job.result().get_counts()
    if '00' in counts:
        PA0exp[j] = counts['00']/nshots
    if '01' in counts:
        PA0exp[j] += counts['01']/nshots
    if '10' in counts:
        PA1exp[j] = counts['10']/nshots
    if '11' in counts:
        PA1exp[j] += counts['11']/nshots
    if PA0exp[j] != 0 and '00' in counts:
        PD0A0exp[j] = (counts['00']/nshots)/PA0exp[j]
    if PA1exp[j] != 0 and '10' in counts:
        PD0A1exp[j] = (counts['10']/nshots)/PA1exp[j]
print(PD0A0exp,PD0A1exp)

import numpy as np; import math
phmax = 2*math.pi; npt = 100; dpht = phmax/npt
pht = np.arange(0,phmax+dpht,dpht)
PD0teoA0 = 0.5*(1+np.cos(pht)); PD0teoA1 = 0.5*(1-np.cos(pht))
from matplotlib import pyplot as plt
plt.plot(pht,PD0teoA0, label=r'$Pr(D_0|A=0)_{teo}$'); plt.plot(pht,PD0teoA1, label=r'$Pr(D_0|A=1)_{teo}$')
plt.plot(ph,PD0A0sim, 'o', label=r'$Pr(D_0|A=0)_{sim}$'); plt.plot(ph,PD0A1sim, 's', label=r'$Pr(D_0|A=1)_{sim}$')
plt.plot(ph,PD0A0exp, '*', label=r'$Pr(D_0|A=0)_{exp}$'); plt.plot(ph,PD0A1exp, '+', label=r'$Pr(D_0|A=1)_{exp}$')
plt.xlim(-0.1,2*math.pi+0.1); plt.ylim(-0.02,1.02); plt.legend(bbox_to_anchor=(1.4, 1.0),loc='upper right')
plt.xlabel(r'$\phi$'); plt.show()

VsimA0 = (np.max(PD0A0sim)-np.min(PD0A0sim))/(np.max(PD0A0sim)+np.min(PD0A0sim))
print('VsimA0 = ',VsimA0)
VsimA1 = (np.max(PD0A1sim)-np.min(PD0A1sim))/(np.max(PD0A1sim)+np.min(PD0A1sim))
print('VsimA1 = ',VsimA1)
VexpA0 = (np.max(PD0A0exp)-np.min(PD0A0exp))/(np.max(PD0A0exp)+np.min(PD0A0exp))
print('VexpA0 = ',VexpA0)
VexpA1 = (np.max(PD0A1exp)-np.min(PD0A1exp))/(np.max(PD0A1exp)+np.min(PD0A1exp))
print('VexpA1 = ',VexpA1)

import qiskit
qiskit.IBMQ.save_account('e20d3045c7c23a3be753606794714e9b47d4331ec46a41b2281e2d601bae27b2b3c1baa0b5cb54f3282a8b5248e552b14c9c556699010660fb94b11d00ff1073', 
                  overwrite = True)
qiskit.IBMQ.load_account()
#provider = qiskit.IBMQ.get_provider(hub='ibm-q', group='open', project='main')
provider = qiskit.IBMQ.get_provider(hub='ibm-q-research-2', group='federal-uni-sant-1', project='main')
device = provider.get_backend('ibm_nairobi')
nshots = 2**13

from qiskit import QuantumRegister, execute
from qiskit.ignis.mitigation.measurement import complete_meas_cal
qr = QuantumRegister(3)
qubit_list = [0, 2]
meas_calibs, state_labels = complete_meas_cal(qubit_list = qubit_list, qr = qr)

job = execute(meas_calibs, backend=device, shots=nshots)
print(job.job_id())

from qiskit.ignis.mitigation.measurement import CompleteMeasFitter
job = device.retrieve_job('cmp2xzqyj8fg0080m0e0')
cal_results = job.result()
meas_fitter = CompleteMeasFitter(cal_results, state_labels)

f = open("jobs_ids_apagador.txt","r")
list_ids_apagador = f.read().replace("'","").replace(" ","").replace("[","").replace("]","").split(",")
f.close()
phmax = 2*math.pi; npe=6; dth = phmax/npe; ph = np.arange(0,phmax+dth,dth)
PD0A0exp_mit = np.zeros(len(ph)); PD0A1exp_mit = np.zeros(len(ph)); PA0exp_mit = np.zeros(len(ph)); PA1exp_mit = np.zeros(len(ph))
device = provider.get_backend('ibm_nairobi')
for j in range(0,len(ph)):
    job = device.retrieve_job(list_ids_apagador[j])
    mitigated_results = meas_fitter.filter.apply(job.result())
    mitigated_counts = mitigated_results.get_counts()
    if '00' in counts:
        PA0exp_mit[j] = mitigated_counts['00']/nshots
    if '01' in counts:
        PA0exp_mit[j] += mitigated_counts['01']/nshots
    if '10' in counts:
        PA1exp_mit[j] = mitigated_counts['10']/nshots
    if '11' in counts:
        PA1exp_mit[j] += mitigated_counts['11']/nshots
    if PA0exp_mit[j] != 0 and '00' in mitigated_counts:
        PD0A0exp_mit[j] = (mitigated_counts['00']/nshots)/PA0exp[j]
    if PA1exp[j] != 0 and '10' in mitigated_counts:
        PD0A1exp_mit[j] = (mitigated_counts['10']/nshots)/PA1exp[j]
print(PD0A0exp_mit, PD0A1exp_mit)

import matplotlib; matplotlib.rcParams.update({'font.size':11}); plt.figure(figsize = (7,5), dpi = 100)
import numpy as np; import math; phmax = 2*math.pi; npt = 100; dpht = phmax/npt; pht = np.arange(0,phmax+dpht,dpht)
PD0teoA0 = 0.5*(1+np.cos(pht)); PD0teoA1 = 0.5*(1-np.cos(pht)); from matplotlib import pyplot as plt
plt.plot(pht,PD0teoA0, label=r'$Pr(D_0|A=0)_{teo}$'); plt.plot(pht,PD0teoA1, label=r'$Pr(D_0|A=1)_{teo}$')
plt.plot(ph,PD0A0sim, 'o', label=r'$Pr(D_0|A=0)_{sim}$'); plt.plot(ph,PD0A1sim, 's', label=r'$Pr(D_0|A=1)_{sim}$')
plt.plot(ph,PD0A0exp, '*', label=r'$Pr(D_0|A=0)_{exp}$'); plt.plot(ph,PD0A1exp, 'v', label=r'$Pr(D_0|A=1)_{exp}$')
plt.plot(ph,PD0A0exp_mit, '^', label=r'$Pr(D_0|A=0)_{expmit}$'); plt.plot(ph,PD0A1exp_mit, '+', label=r'$Pr(D_0|A=1)_{expmit}$')
plt.xlim(-0.1,2*math.pi+0.1); plt.ylim(-0.03,1.03); plt.legend(bbox_to_anchor=(1.4, 1.0),loc='upper right')
plt.xlabel(r'$\phi$'); plt.show()

VsimA0 = (np.max(PD0A0sim)-np.min(PD0A0sim))/(np.max(PD0A0sim)+np.min(PD0A0sim)); print('VsimA0 = ',VsimA0)
VsimA1 = (np.max(PD0A1sim)-np.min(PD0A1sim))/(np.max(PD0A1sim)+np.min(PD0A1sim)); print('VsimA1 = ',VsimA1)
VexpA0 = (np.max(PD0A0exp)-np.min(PD0A0exp))/(np.max(PD0A0exp)+np.min(PD0A0exp)); print('VexpA0 = ',VexpA0)
VexpA1 = (np.max(PD0A1exp)-np.min(PD0A1exp))/(np.max(PD0A1exp)+np.min(PD0A1exp)); print('VexpA1 = ',VexpA1)
VexpA0_mit = (np.max(PD0A0exp_mit)-np.min(PD0A0exp_mit))/(np.max(PD0A0exp_mit)+np.min(PD0A0exp_mit))
print('VexpA0_mit = ',VexpA0_mit)
VexpA1_mit = (np.max(PD0A1exp_mit)-np.min(PD0A1exp_mit))/(np.max(PD0A1exp_mit)+np.min(PD0A1exp_mit))
print('VexpA1_mit = ',VexpA1_mit)

import math; import numpy as np
al_max = math.pi/4; npt = 100; dalt = al_max/npt
alt = np.arange(0,al_max+dalt,dalt)
Vteo = np.sin(2*alt)
from matplotlib import pyplot as plt
plt.plot(alt,Vteo, label=r'$V_{teo}(A=0)$')
plt.xlim(-0.01,al_max+0.01); plt.ylim(-0.01,1.01)
plt.xlabel(r'$\alpha$'); plt.legend()
plt.show()

from qiskit import QuantumCircuit

def qc_Psi(al): # prepara o estado parcialmente emaranhado
    qc = QuantumCircuit(2, name = r'$|\Psi\rangle$')
    qc.u(2*al,0,0, 0)
    qc.cx(0,1)
    qc.z(0)
    qc.x(1)
    return qc
qcpsi = qc_Psi(0.1); qcpsi.draw('mpl')

# simulação clássica
from qiskit import Aer, execute; simulator = Aer.get_backend('qasm_simulator'); nshots = 2**13
import math; import numpy as np
almax = math.pi/4; npa=2; dal = almax/npa; al = np.arange(0,almax+dal,dal)
phmax = 2*math.pi; npe=4; dph = phmax/npe; ph = np.arange(0,phmax+dph,dph); VA0sim = np.zeros(len(al))
for j in range(0,len(al)):
    PD0A0sim = np.zeros(len(ph)); PA0sim = np.zeros(len(ph))
    for k in range(0,len(ph)):
        qc = QuantumCircuit(3,2)
        qcpsi = qc_Psi(al[j]); qc.append(qcpsi, [0,1])
        qcdfp = qc_dfp(); qc.append(qcdfp, [1,2]); qcpmo = qc_pmo(); qc.append(qcpmo, [1,2])
        qce = qc_espelhos(); qc.append(qce, [2]); qc.p(ph[k], [2])
        qcdf = qc_df(); qc.append(qcdf, [2]); qcpqo = qc_pqo(); qc.append(qcpqo, [0])
        qc.measure([0,2],[0,1])
        job = execute(qc, backend=simulator, shots=nshots); counts = job.result().get_counts()
        if '00' in counts:
            PA0sim[k] = counts['00']/nshots
        if '01' in counts:
            PA0sim[k] += counts['01']/nshots
        if PA0sim[k] != 0:
            if '00' in counts:
                PD0A0sim[k] = (counts['00']/nshots)/PA0sim[k]
    VA0sim[j] = (np.max(PD0A0sim)-np.min(PD0A0sim))/(np.max(PD0A0sim)+np.min(PD0A0sim))
VA0sim

PD0A0sim

qc.draw('mpl')

import math; import numpy as np
al_max = math.pi/4; npt = 100; dalt = al_max/npt; alt = np.arange(0,al_max+dalt,dalt)
Vteo = np.sin(2*alt)
from matplotlib import pyplot as plt
plt.plot(alt,Vteo, label=r'$V_{teo}(A=0)$')
plt.plot(al,VA0sim, 'o', label=r'$V_{sim}(A=0)$')
plt.xlim(-0.01,al_max+0.01); plt.ylim(-0.01,1.02); plt.xlabel(r'$\alpha$'); plt.legend()
plt.show()

# Experimento (simulação quântica)
from qiskit import execute; import math; import numpy as np
almax = math.pi/4; npa = 2; dal = almax/npa; al = np.arange(0,almax+dal,dal)
phmax = 2*math.pi; npe=4; dph = phmax/npe; ph = np.arange(0,phmax+dph,dph)
job_ids=[]
for j in range(0,len(al)):
    for k in range(0,len(ph)):
        qc = QuantumCircuit(3,2); qcpsi = qc_Psi(al[j]); qc.append(qcpsi, [0,1])
        qcdfp = qc_dfp(); qc.append(qcdfp, [1,2]); qcpmo = qc_pmo(); qc.append(qcpmo, [1,2])
        qce = qc_espelhos(); qc.append(qce, [2]); qc.p(ph[k], [2])
        qcdf = qc_df(); qc.append(qcdf, [2]); qcpqo = qc_pqo(); qc.append(qcpqo, [0]); qc.measure([0,2],[0,1])
        job = execute(qc, backend=device, shots=nshots); jobid = job.job_id(); job_ids.append(jobid); print(jobid)
f = open("jobs_ids_apagador_parcial.txt", "w"); f.write(str(job_ids)); f.close()    

f = open("jobs_ids_apagador_parcial.txt","r")
list_ids_apagador_parcial = f.read().replace("'","").replace(" ","").replace("[","").replace("]","").split(",")
f.close()
import math; import numpy as np
almax = math.pi/4; npa=2; dal = almax/npa; al = np.arange(0,almax+dal,dal)
phmax = 2*math.pi; npe=4; dph = phmax/npe; ph = np.arange(0,phmax+dph,dph); VA0exp = np.zeros(len(al))
for j in range(0,len(al)):
    PD0A0exp = np.zeros(len(ph)); PA0exp = np.zeros(len(ph))
    for k in range(0,len(ph)):
        job = device.retrieve_job(list_ids_apagador_parcial[(npe+1)*j+k])
        counts = job.result().get_counts()
        if '00' in counts:
            PA0exp[k] = counts['00']/nshots
        if '01' in counts:
            PA0exp[k] += counts['01']/nshots
        if PA0exp[k] != 0:
            if '00' in counts:
                PD0A0exp[k] = (counts['00']/nshots)/PA0exp[k]
    VA0exp[j] = (np.max(PD0A0exp)-np.min(PD0A0exp))/(np.max(PD0A0exp)+np.min(PD0A0exp))
VA0exp

import math; import numpy as np
al_max = math.pi/4; npt = 100; dalt = al_max/npt; alt = np.arange(0,al_max+dalt,dalt); Vteo = np.sin(2*alt)
from matplotlib import pyplot as plt
plt.plot(alt,Vteo, label=r'$V_{teo}(A=0)$')
plt.plot(al,VA0sim, 'o', label=r'$V_{sim}(A=0)$')
plt.plot(al,VA0exp, '^', label=r'$V_{exp}(A=0)$')
plt.xlim(-0.01,al_max+0.01); plt.ylim(-0.01,1.02); plt.xlabel(r'$\alpha$'); plt.legend()
plt.show()

