%run init.ipynb

import qiskit
from qiskit import *
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

def state_prep_psip():
    qr = QuantumRegister(2); qc = QuantumCircuit(qr, name = 'sprep')
    qc.h(qr[0]); qc.cx(qr[0],qr[1]); qc.x(qr[1]) # prepara psi+_01
    return qc
state_prep_psip_ = state_prep_psip(); state_prep_psip_.draw(output = 'mpl')

def qc_dfp():
    qr = QuantumRegister(2); qc = QuantumCircuit(qr, name='DFP')
    qc.cz(qr[0],qr[1]); qc.cy(qr[0], qr[1])
    return qc
qc_dfp_ = qc_dfp(); qc_dfp_.draw(output = 'mpl')

def qc_df():
    qr = QuantumRegister(1); qc = QuantumCircuit(qr, name='DF')
    qc.s(qr[0]); qc.h(qr[0]); qc.s(qr[0])
    return qc
qc_df_ = qc_df(); qc_df_.draw(output = 'mpl')

def qc_pmo():
    qr = QuantumRegister(2); qc = QuantumCircuit(qr, name='PMO')
    qc.cx(qr[1],qr[0])
    return qc
qc_pmo_ = qc_pmo(); qc_pmo_.draw(output = 'mpl')

def qc_pqo():
    qr = QuantumRegister(1); qc = QuantumCircuit(qr, name='PQO')
    qc.h(qr[0]); qc.s(qr[0])
    return qc
qc_pqo_ = qc_pqo(); qc_pqo_.draw(output = 'mpl')

def qc_espelho():
    qr = QuantumRegister(1); qc = QuantumCircuit(qr, name='espelho')
    qc.z(qr[0]); qc.y(qr[0])
    return qc
qc_espelho_ = qc_espelho(); qc_espelho_.draw(output = 'mpl')

import math
def qc_fase(ph):
    qr = QuantumRegister(1); qc = QuantumCircuit(qr, name='fase')
    qc.p(ph, qr[0])
    return qc
qc_fase_ = qc_fase(math.pi); qc_fase_.draw(output = 'mpl')

def qc_apagador(ph):
    qr = QuantumRegister(3); qc = QuantumCircuit(qr, name = 'apagador')
    state_prep_psip_ = state_prep_psip(); qc.append(state_prep_psip_, [qr[0],qr[1]]) # sprep
    qc.barrier()
    qc_dfp_ = qc_dfp(); qc.append(qc_dfp_, [qr[1],qr[2]]) # DFP
    qc_pmo_ = qc_pmo(); qc.append(qc_pmo_, [qr[1],qr[2]]) # PMO
    qc_pqo_ = qc_pqo(); qc.append(qc_pqo_, [qr[0]]) # PQO
    qc_espelho_ = qc_espelho(); qc.append(qc_espelho_, [qr[2]]) # espelho
    qc_fase_ = qc_fase(ph); qc.append(qc_fase_, [qr[2]]) # fase
    qc_df_ = qc_df(); qc.append(qc_df_, [qr[2]]) # DF
    return qc
qc_apagador_ = qc_apagador(math.pi/3)
qc_apagador_.draw(output = 'mpl')
qc_apagador_.decompose().draw(output = 'mpl')

import numpy as np
from numpy import random
phmax = 2*math.pi; dph = phmax/20; ph = np.arange(0, phmax+dph, dph)#; ph
d = ph.shape[0]; PD0teo = np.zeros(d); PD0sim = np.zeros(d); PD0exp = np.zeros(d)
ri = random.randint(0,d)
for j in range(0, d):
    PD0teo[j] = 1/2 # theoretical
    qr = QuantumRegister(3); qc = QuantumCircuit(qr)
    qc_apagador_ = qc_apagador(ph[j]); qc.append(qc_apagador_, [qr[0],qr[1],qr[2]]) # apagador q
    qstc = state_tomography_circuits(qc, qr[2]) # 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')
    PD0sim[j] = rho[0,0].real
    #if j == ri:
    #    print(rho)
    # 1º fiz as simulações, só depois adicionei o código para os experimentos
    '''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')
    PD0exp[j] = rho[0,0].real'''

import matplotlib 
from matplotlib import pyplot as plt
matplotlib.rcParams.update({'font.size':12}); plt.figure(figsize = (6,4), dpi = 100)
plt.plot(ph, PD0teo, label = r'$P_{teo}(D_{0})$')
plt.plot(ph, PD0sim, 'o', label = r'$P_{sim}(D_{0})$')
#plt.plot(ph, PD0exp, '*', label = r'$P_{exp}(D_{0})$')
plt.legend(); plt.xlabel(r'$\phi$'); plt.ylim(-0.02,1.02); 
plt.title('IMZ com info de caminho'); plt.show()
# Nesse caso, não foram feitas medidas projetivas para obtenção de informação de caminho 

qr = QuantumRegister(1); cr = ClassicalRegister(1); qc = QuantumCircuit(qr,cr)
qc.h(qr[0]); qc.measure(qr[0],cr[0])
job = qiskit.execute(qc, Aer.get_backend('qasm_simulator'), shots=10)
res = job.result().get_counts()
if '0' in res:
    print('z = 0, ', res, ', ', res['0'])
elif '1' in res:
    print('z = 1, ', res, ', ', res['1'])
qc.draw(output = 'mpl')

for j in range(0,1):
    qr = QuantumRegister(2); cr = ClassicalRegister(2); qc = QuantumCircuit(qr,cr)
    qc.h(qr[0]);
    qc.measure(qr[0],cr[0])
    qc.h(qr[1]); qc.measure(qr[1],cr[1])
    job = qiskit.execute(qc, Aer.get_backend('qasm_simulator'), shots=10)
    res = job.result().get_counts()
    print(res)
    if '00' in res:
        P00 = 1; print('P00',P00)
    elif '01' in res:
        P01 = 1; print('P01',P01)
qc.draw(output = 'mpl')

# O apagador quântico mesmo, com as medidas
def apagadorM(ph):
    qr = QuantumRegister(3); cr = ClassicalRegister(2); qc = QuantumCircuit(qr, cr)
    qc_apagador_ = qc_apagador(ph); qc.append(qc_apagador_, [qr[0],qr[1],qr[2]])
    qc.measure(qr[0],cr[0]) # medida da polarização de A
    qc.measure(qr[2],cr[1]) # medida do caminho de B
    return qc
apagadorM_ = apagadorM(math.pi/8); apagadorM_.draw(output = 'mpl')

qc = apagadorM(math.pi/8)
job = qiskit.execute(qc, Aer.get_backend('qasm_simulator'), shots=1)
res = job.result().get_counts(); print(res)

nshots = 1000
phmax = 2*math.pi; dph = phmax/20; ph = np.arange(0, phmax+dph, dph); d = ph.shape[0]
PD0teo0 = np.zeros(d); PD0sim0 = np.zeros(d); PD0exp0 = np.zeros(d)
PD0teo1 = np.zeros(d); PD0sim1 = np.zeros(d); PD0exp1 = np.zeros(d)
for j in range(0, d):
    PD0teo0[j] = (1 + math.cos(ph[j]))/2; PD0teo1[j] = (1 - math.cos(ph[j]))/2
    qc = apagadorM(ph[j])
    N0 = 0; N1 = 0; N0e = 0; N1e = 0
    for k in range(0, nshots):  # não usar esse loop
        job = qiskit.execute(qc, Aer.get_backend('qasm_simulator'), shots=1)
        res = job.result().get_counts()#; print(res)
        if '00' in res or '10' in res:
            N0 += 1
        elif '01' in res or '11' in res:
            N1 += 1
        if '00' in res:
            PD0sim0[j] += 1
        elif '01' in res:
            PD0sim1[j] += 1
        '''job = qiskit.execute(qc, Aer.get_backend(device), shots=1)
        res = job.result().get_counts()
        if '00' in res or '10' in res:
            N0e += 1
        elif '01' in res or '11' in res:
            N1e += 1
        if '00' in res:
            PD0exp0[j] += 1
        elif '01' in res:
            PD0exp1[j] += 1'''
PD0sim0 = PD0sim0/N0; PD0sim1 = PD0sim1/N1
#PD0exp0 = PD0exp0/N0e; PD0exp1 = PD0exp1/N1e

# sem o loop em k

matplotlib.rcParams.update({'font.size':12}); plt.figure(figsize = (5,3.3), dpi = 100)
plt.plot(ph, PD0teo0, '-', label = r'$P_{teo}(0,0)$')
plt.plot(ph, PD0teo1, '-.', label = r'$P_{teo}(1,0)$')
plt.plot(ph, PD0sim0, 'o', label = r'$2P_{sim}(0,0)$')
plt.plot(ph, PD0sim1, '*', label = r'$2P_{sim}(1,0)$')
#plt.plot(ph, PD0exp0, '+', label = r'$2P_{exp}(0,0)$')
#plt.plot(ph, PD0exp1, 'd', label = r'$2P_{exp}(1,0)$')
plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')
plt.xlabel(r'$\phi$'); plt.ylim(-0.05,1.05); 
plt.title('IMZ com apagador de info de caminho'); plt.show()

def qc_MBB():
    qr = QuantumRegister(2); qc = QuantumCircuit(qr)
    qc.h(qr[0]); qc.cx(qr[0],qr[1])
    return qc
qc_MBB_ = qc_MBB(); qc_MBB_.draw(output = 'mpl')

def qc_apagadorE(ph):
    qr = QuantumRegister(3); qc = QuantumCircuit(qr, name = 'apagador')
    state_prep_psip_ = state_prep_psip(); qc.append(state_prep_psip_, [qr[0],qr[1]]) # sprep
    qc_dfp_ = qc_dfp(); qc.append(qc_dfp_, [qr[1],qr[2]]) # DFP
    qc_pmo_ = qc_pmo(); qc.append(qc_pmo_, [qr[1],qr[2]]) # PMO
    qc_pqo_ = qc_pqo(); qc.append(qc_pqo_, [qr[0]]) # PQO
    qc_espelho_ = qc_espelho(); qc.append(qc_espelho_, [qr[2]]) # espelho
    qc_fase_ = qc_fase(ph); qc.append(qc_fase_, [qr[2]]) # fase
    qc_df_ = qc_df(); qc.append(qc_df_, [qr[2]]) # DF
    return qc
qc_apagador_ = qc_apagador(math.pi/3); qc_apagador_.draw(output = 'mpl')