%run init.ipynb

Udfp = (1/sqrt(2))*Matrix([[1,1j,0,0],[1j,1,0,0],[0,0,1,1j],[0,0,1j,1]]); Udfp

Udfp*conjugate(Udfp), conjugate(Udfp)*Udfp # ok!

Udf = (1/sqrt(2))*Matrix([[1,1j],[1j,1]]); Udf

X = Matrix([[0,1],[1,0]]); C_Udf = tp(proj(cb(2,0)),id(2)) + tp(proj(cb(2,1)),Udf)
C_Udf, tp(X,id(2))*C_Udf*tp(X,id(2)), tp(X,id(2))*C_Udf*tp(X,id(2))*C_Udf, C_Udf*tp(X,id(2))*C_Udf*tp(X,id(2))

C_Udf*tp(X,id(2))*C_Udf*tp(X,id(2))

def P(al):
    return Matrix([[1,0],[0,cos(al)+1j*sin(al)]])
S = Matrix([[1,0],[0,1j]]); Sd = Matrix([[1,0],[0,-1j]]); H = (1/sqrt(2))*Matrix([[1,1],[1,-1]])
X = Matrix([[0,1],[1,0]]); Y = Matrix([[0,-1j],[1j,0]]); Z = Matrix([[1,0],[0,-1]])
def U(th,ph,lb):
    return Matrix([[cos(th/2),-(cos(lb)+1j*sin(lb))*sin(th/2)],[(cos(ph)+1j*sin(ph))*sin(th/2),cos(th/2)]])
th, ph, lb = symbols('theta phi lambda'); U(pi,-pi/2,pi/2)*S*X*H*X*S

A = U(pi,-pi/2,pi/2)*S; B = H; C = S; A*B*C  # = Udf^dagger != I

a,b,c,d, x, y, z, w = symbols('a,b,c,d, x, y, z, w')
A = Matrix([[a,b],[c,d]]); B = Matrix([[x,y],[z,w]]); A, B, A*B

a = -1/sqrt(2); b = 1j/sqrt(2); c = b; d = a # para trocar o sinal dos elementos da diagonal principal
solve([a*x+b*z+a, a*y+b*w-b, c*x+d*z-c, c*y+d*w+d], [x,y,z,w])

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(th, ph):
    qr = QuantumRegister(2); qc = QuantumCircuit(qr, name = 'prep')
    qc.u(th, ph, 0, qr[0])
    qc.u(th, ph, 0, qr[1])
    return qc

def qc_dfp():
    qr = QuantumRegister(2); qc = QuantumCircuit(qr, name='DFP')
    '''qc.x(qr[1])
    qc.cu(math.pi/2, math.pi/2, -math.pi/2, 0.0, qr[1], qr[0]) # Udf controlada
    qc.x(qr[1])
    qc.cu(math.pi/2, math.pi/2, -math.pi/2, 0.0, qr[1], qr[0])'''
    qc.cu(math.pi, math.pi/2, -math.pi/2, 0.0, qr[0], qr[1])
    return qc
qc_dfp_ = qc_dfp(); qc_dfp_.draw()

###### Para ver a matriz unitário implementada por um certo circuito quântico ##########
backend = Aer.get_backend('unitary_simulator'); qc_dfp_ = qc_dfp(); job = execute(qc_dfp_, backend)
result = job.result(); print(result.get_unitary(qc_dfp_, decimals=3))
# agora tá certo, mas ao contrário em relação aos qubits...

qr = QuantumRegister(2); qc = QuantumCircuit(qr, name='DFP')
#qc.cx(qr[0],qr[1]) # Sai a matrix para C21
qc.cx(qr[1],qr[0]) # Sai a matrix para C12 
job = execute(qc, backend); result = job.result(); print(result.get_unitary(qc, decimals=3))

qr = QuantumRegister(2); qc = QuantumCircuit(qr, name='DFP')
qc.h(qr[0])
#qc.cx(qr[0], qr[1]) # sai o estado certo, |00>+|11>
qc.cx(qr[1], qr[0]) # sai |+0>
qstc = state_tomography_circuits(qc, [qr[0],qr[1]])
job = qiskit.execute(qstc, Aer.get_backend('qasm_simulator'), shots=nshots)
qstf = StateTomographyFitter(job.result(), qstc); rho = qstf.fit(method='lstsq')
rho

# Verificação da Udf controlada
qr = QuantumRegister(2); qc = QuantumCircuit(qr, name='DFP')
qc.x(qr[1])
qc.cu(math.pi/2, math.pi/2, -math.pi/2, 0.0, qr[1], qr[0]) # a sequência é target depois control
qc.x(qr[1])
job = execute(qc, backend); result = job.result(); print(result.get_unitary(qc, decimals=3))

qc.draw()

# Verificação da Udf controlada
qr = QuantumRegister(2); qc = QuantumCircuit(qr, name='DFP')
qc.cx(qr[1], qr[0]) # a sequência é target depois control
job = execute(qc, backend); result = job.result(); print(result.get_unitary(qc, decimals=3))

def qc_mzi_dfp(ph): # quantum circuit for the Mach-Zehnder interferometer com DFP
    qr = QuantumRegister(2); qc = QuantumCircuit(qr)
    state_prep_ = state_prep(0, 0); qc.append(state_prep_, [qr[0],qr[1]]) # estado prep
    qc_dfp_ = qc_dfp(); qc.append(qc_dfp_, [qr[0],qr[1]]) # DFP
    qc.z(qr[1]); qc.y(qr[1]) # espelho
    qc.p(ph, qr[1]); # fase
    qc.s(qr[1]); qc.h(qr[1]); qc.s(qr[1]) # DF
    return qc
qcmzi = qc_mzi_dfp(math.pi/8); qcmzi.draw()

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)
for j in range(0, d):
    PD0teo[j] = 1/2 # theoretical
    qr = QuantumRegister(2); qc = QuantumCircuit(qr)
    state_prep_ = state_prep(0, 0); qc.append(state_prep_, [qr[0],qr[1]]) # state preparation (estado |00>)
    qc_mzi_dfp_ = qc_mzi_dfp(ph[j]); qc.append(qc_mzi_dfp_, [qr[0],qr[1]]) # apply sim mzip 
    qstc = state_tomography_circuits(qc, qr[1]) # 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
    #print(rho[0,0].real)
    # 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'''

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.show()

def state_prep_psip():
    qr = QuantumRegister(3); qc = QuantumCircuit(qr, name = 'prep')
    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()

def qc_dfp():
    qr = QuantumRegister(2); qc = QuantumCircuit(qr, name='DFP')
    qc.cu(math.pi, math.pi/2, -math.pi/2, 0.0, qr[0], qr[1])
    return qc
qc_dfp_ = qc_dfp(); qc_dfp_.draw()

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()

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()

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()

def qc_apagador(ph):
    qr = QuantumRegister(3); qc = QuantumCircuit(qr)
    state_prep_psip_ = state_prep_psip(); qc.append(state_prep_psip_, [qr[0],qr[1],qr[2]])
    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]])
    qc.z(qr[2]); qc.y(qr[2]) # espelho
    qc.p(ph, 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()

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)
for j in range(0, d):
    PD0teo[j] = (2+math.cos(ph[j])+math.sin(ph[j]))/4 # 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
    #print(rho[0,0].real)
    # 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'''

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.show()

H = (1/sqrt(2))*Matrix([[1,1],[1,-1]]); S = Matrix([[1,0],[0,1j]]); X = Matrix([[0,1],[1,0]])
def U(th, ph, lb):
    return Matrix([[cos(th/2),-(cos(lb)+1j*sin(lb))*sin(th/2)],
                   [(cos(ph)+1j*sin(ph))*sin(th/2),(cos(ph+lb)+1j*sin(ph+lb))*cos(th/2)]])
def dfp():
    return tp(proj(cb(2,0)), id(2)) + tp(proj(cb(2,1)), U(math.pi, math.pi/2, -math.pi/2))
def pmo():
    return tp(id(2), proj(cb(2,0))) + tp(X, proj(cb(2,1)))

th, ph, lb = symbols('theta phi lambda');# U(th, ph, lb)

psi0 = tp(tp(cb(2,0),cb(2,0)),cb(2,0))#; psi0
dfp_ = dfp(); psi1 = tp(id(2), dfp_)*psi0; psi1
pmo_ = pmo(); psi2 = tp(id(2), pmo_)*psi1; psi2

