%run init.ipynb

from qiskit import *

nshots = 8192
qiskit.IBMQ.load_account()
provider = qiskit.IBMQ.get_provider(hub = 'ibm-q-research-2', group = 'federal-uni-sant-1', project = 'main')
simulator = Aer.get_backend('qasm_simulator')
device = provider.get_backend('ibmq_bogota')
from qiskit.tools.visualization import plot_histogram
from qiskit.ignis.verification.tomography import state_tomography_circuits, StateTomographyFitter
from qiskit.tools.monitor import job_monitor, backend_overview, backend_monitor

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

def qc_bomb():
    qc = QuantumCircuit(2, name = 'bomba')
    qc.cx(0,1)
    return qc
qc_bomb_ = qc_bomb(); qc_bomb_.draw(output = 'mpl')

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

qc = QuantumCircuit(2, 2)
qc_df_ = qc_df(); qc.append(qc_df_, [0])
qc_espelho_ = qc_espelho(); qc.append(qc_espelho_, [0])
qc_bomb_ = qc_bomb(); qc.append(qc_bomb_, [0,1])
qc.append(qc_df_, [0]); qc.barrier()
qc.measure([0,1],[0,1])
qc.draw(output = 'mpl')

job = qiskit.execute(qc, backend = simulator, shots = nshots)
plot_histogram(job.result().get_counts(qc))

job = qiskit.execute(qc, backend = device, shots = nshots)
print(job.job_id()); job_monitor(job)

plot_histogram(job.result().get_counts(qc))

T, R = symbols('T R'); Udfg = Matrix([[T,1j*R],[1j*R,T]]); Udfgd = Matrix([[T,-1j*R],[-1j*R,T]])
Udfg*Udfgd, Udfgd*Udfg # ok, é unitária

def qc_dfg(th):
    qc = QuantumCircuit(1, name = 'DFG')
    qc.rx(-2*th, 0)
    return qc
qc_dfg_ = qc_dfg(math.pi); qc_dfg_.draw(output = 'mpl')

def qc_imzg(th, ph):
    qc = QuantumCircuit(1, name = 'IMZG')
    qc_dfg_ = qc_dfg(th); qc.append(qc_dfg_, [0])
    qc_espelho_ = qc_espelho(); qc.append(qc_espelho_, [0])
    qc.p(ph, 0)
    qc_dfg_ = qc_dfg(th); qc.append(qc_dfg_, [0])
    return qc
qc_imzg_ = qc_imzg(math.pi, math.pi/2); qc_imzg_.draw(output = 'mpl')

phmax = 2*math.pi; dph = phmax/20; ph = np.arange(0, phmax+dph, dph)#; ph
th = math.pi/8
print('T = ', math.cos(th), ', R = ', math.sin(th))
d = ph.shape[0]; PD0teo = np.zeros(d); PD0sim = np.zeros(d);
for j in range(0, d):
    PD0teo[j] = 2*math.sin(th)**2*math.cos(th)**2*(1+math.cos(ph[j])) # theoretical
    qc = QuantumCircuit(1)
    qc_imzg_ = qc_imzg(th, ph[j]); qc.append(qc_imzg_, [0])
    qstc = state_tomography_circuits(qc, [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.plot(ph, PD0exp, '*', label = r'$P_{exp}(D_{0})$')
plt.legend(); plt.xlabel(r'$\phi$'); plt.show()

def V(pd):
    return (np.max(pd)-np.min(pd))/(np.max(pd)+np.min(pd))
def Vt(pd):
    return np.max(pd)*(np.max(pd)-np.min(pd))/(np.max(pd)+np.min(pd))
def pd_ph(th):
    phmax = 2*math.pi; dph = phmax/100; ph = np.arange(0, phmax+dph, dph)
    return 2*math.sin(th)**2*math.cos(th)**2*(1+np.cos(ph))

dth = 0.05; th = np.arange(dth, math.pi/2, dth); d = th.shape[0]
Vis = np.zeros(d); VisT= np.zeros(d); Vteo = np.zeros(d)
for j in range(0, d):
    Vteo[j] = 4*math.cos(th[j])**2*math.sin(th[j])**2
    pd = pd_ph(th[j]); Vis[j] = V(pd); VisT[j] = Vt(pd)
    
matplotlib.rcParams.update({'font.size':12}); plt.figure(figsize = (6,4), dpi = 100)
plt.plot(th, Vis, label = r'$V$'); plt.plot(th, VisT, label = r'$\tilde{V}$')
plt.plot(th, Vteo, '*', label = r'$V_{teo}$')
plt.xlabel(r'$\theta$'); plt.legend(); plt.show()

float(cos(pi/12)**2)

float((cos(pi/12)**2)**6), float((cos(pi/24)**2)**12), float((cos(pi/48)**2)**24), float((cos(pi/96)**2)**48)

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

T = cos(math.pi/12); T

th = +pi-pi/6; float((sin(th/2)**2)**6)

T = cos(th); R = sin(th); print('T = ', float(T), ', R = ', float(R)) # ??????????????

N = 6; th = math.pi-math.pi/N; qc = QuantumCircuit(2, N)
for j in range(0, N-1):
    qc.rx(th, 0)
    qc_espelho_ = qc_espelho(); qc.append(qc_espelho_, [0])
    qc.cx(0, 1); qc.measure(1, j); qc.reset(1)
    qc.barrier()
qc.rx(th, 0) # o último DFG é colocado só na última volta
qc.measure(0, N-1)
#qc.draw(output = 'mpl')
qc.decompose().draw(output = 'mpl')

matplotlib.rcParams.update({'font.size':12})
#plt.figure(figsize = (6,4), dpi = 100)
result = execute(qc, backend = simulator, shots = nshots).result()
plot_histogram(result.get_counts(qc))

job = execute(qc, backend = device, shots = nshots)
job_monitor(job)

plot_histogram(job.result().get_counts(qc))





import numpy as np
for j in range(0,20):
    print(np.random.randint(11))

