%run init.ipynb

S = Matrix([[1,0],[0,1j]]); Sd = Matrix([[1,0],[0,-1j]]); H = (1/sqrt(2))*Matrix([[1,1],[1,-1]])#; S, H
S*H*S

from qiskit import *

def qc_mzi_sim(ph): # quantum circuit for the Mach-Zehnder interferometer
    qr = QuantumRegister(1); qc = QuantumCircuit(qr)
    qc.u(math.pi/2, math.pi/2, -math.pi/2, qr[0]); # DF
    #qc.s(qr[0]); qc.h(qr[0]); qc.s(qr[0]) # DF
    qc.z(qr[0]); qc.y(qr[0]) # espelho
    qc.p(ph, qr[0]); # fase
    qc.u(math.pi/2, math.pi/2, -math.pi/2, qr[0]) # DF
    #qc.s(qr[0]); qc.h(qr[0]); qc.s(qr[0]) # DF
    return qc
qcmzi = qc_mzi_sim(math.pi/8); qcmzi.draw()

def state_prep(th, ph):
    qr = QuantumRegister(1); qc = QuantumCircuit(qr)
    qc.u(th, ph, 0, qr[0])
    return qc

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

# for error mitigation
from qiskit.ignis.mitigation.measurement import complete_meas_cal, CompleteMeasFitter
qr = qiskit.QuantumRegister(1)
qubit_list = [0] # os qubits para os quais aplicaremos calibracao de medidas
meas_calibs, state_labels = complete_meas_cal(qubit_list = qubit_list, qr = qr)
job = execute(meas_calibs, backend = device, shots = nshots); job_monitor(job)
cal_results = job.result()
meas_fitter = CompleteMeasFitter(cal_results, state_labels)

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+math.cos(ph[j]))/2 # theoretical
    qr = QuantumRegister(1); qc = QuantumCircuit(qr)
    qc_sp = state_prep(0, 0); qc.append(qc_sp, [qr[0]]) # state preparation (estado |0>)
    qc_mzi = qc_mzi_sim(ph[j]); qc.append(qc_mzi, [qr[0]]) # apply sim mzi 
    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')
    PD0sim[j] = 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.show()

# gráfico para estado inicial |+>
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):
    c0 = 1/math.sqrt(2); c1 = c0; 
    PD0teo[j] = (abs(c0)**2*(1+math.cos(ph[j])) + abs(c1)**2*(1-math.cos(ph[j]))\
                 - 2*((c0*np.conjugate(c1)).real)*math.sin(ph[j]))/2 
    qr = QuantumRegister(1); qc = QuantumCircuit(qr)
    qc_sp = state_prep(math.pi/2, 0); qc.append(qc_sp, [qr[0]]) # state preparation (estado |+>)
    qc_mzi = qc_mzi_sim(ph[j]); qc.append(qc_mzi, [qr[0]]) # 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')
    PD0sim[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.legend(); plt.xlabel(r'$\phi$'); plt.show()

phmax = 4*math.pi; dph = phmax/50; ph = np.arange(0, phmax+dph, dph); d = ph.shape[0]
w = 1; pw1 = (1+(2*w-1)*np.cos(ph))/2; w = 0.25; pw025 = (1+(2*w-1)*np.cos(ph))/2
w = 0.5; pw05 = (1+(2*w-1)*np.cos(ph))/2; w = 0.75; pw075 = (1+(2*w-1)*np.cos(ph))/2
w = 0; pw0 = (1+(2*w-1)*np.cos(ph))/2
matplotlib.rcParams.update({'font.size':12}); plt.figure(figsize = (6,4), dpi = 100)
plt.plot(ph, pw1, label = r'$w=1$'); plt.plot(ph, pw075, label = r'$w=0.75$')
plt.plot(ph, pw05, label = r'$w=0.5$'); plt.plot(ph, pw025, label = r'$w=0.25$')
plt.plot(ph, pw0, label = r'$w=0$')
plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left'); plt.xlabel(r'$\phi$'); plt.show()

def shannon_num(pv):
    d = pv.shape[0]; SE = 0.0; j = -1
    while (j < d-1):
        j = j + 1
        if pv[j] > 10**-15 and pv[j] < (1.0-10**-15):
            SE -= pv[j]*math.log(pv[j], 2)
    return SE
import scipy.linalg.lapack as lapak
def von_neumann_num(rho):
    d = rho.shape[0]; b = lapak.zheevd(rho)
    return shannon_num(b[0])

def V(pd):
    return (np.max(pd)-np.min(pd))/(np.max(pd)+np.min(pd))
def pd_ph(w):
    phmax = 4*math.pi; dph = phmax/100; ph = np.arange(0, phmax+dph, dph)
    return (1+(2*w-1)*np.cos(ph))/2
def rhow(w):
    return np.array([[w,0],[0,1-w]])
w = np.arange(0, 1.05, 0.05); d = w.shape[0]; Vis = np.zeros(d); Svn = np.zeros(d)
for j in range(0, d):
    pd = pd_ph(w[j]); Vis[j] = V(pd); rho = rhow(w[j]); Svn[j] = von_neumann_num(rho)
matplotlib.rcParams.update({'font.size':12}); plt.figure(figsize = (6,4), dpi = 100)
plt.plot(w, Vis, label = 'V'); plt.plot(w, Svn, label = r'$S_{vn}$') 
plt.xlabel(r'$w$'); plt.legend(); plt.show()

