import qiskit
qc = qiskit.QuantumCircuit(4, 2)
qc.h([0, 2]); qc.cx([0], [1]); qc.cx([2], [3]) # preparacao e distribuicao de emaranhamento
qc.barrier() # base de Bell -> base computacional
qc.cx([1], [2]); qc.h([1]) 
qc.barrier() # feedback baseado na info classica
qc.x([1, 2]); qc.cx([2], [3]); qc.cz([1], [3]); qc.x([1, 2]) 
qc.barrier() # para verificacao da funcionalidade do cq (tem que sair Psi-=11)
qc.cx([0], [3]); qc.h([0]); qc.measure([0, 3], [0, 1])
qc.draw(output = 'mpl')


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

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

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

#from qiskit.circuit.Instruction import c_if
from qiskit import *
q0 = QuantumRegister(1); q1 = QuantumRegister(1); q2 = QuantumRegister(1); q3 = QuantumRegister(1)
c0 = ClassicalRegister(1); c1 = ClassicalRegister(1); c2 = ClassicalRegister(1); c3 = ClassicalRegister(1)
qc = qiskit.QuantumCircuit(q0, q1, q2, q3, c0, c1, c2, c3)
qc.h(q0); qc.h(q2); qc.cx(q0, q1); qc.cx(q2, q3) # preparacao e distribuicao de emaranhamento
qc.barrier() # base de Bell -> base computacional
qc.cx(q1, q2); qc.h(q1); qc.measure(q1, c0); qc.measure(q2, c1)
qc.barrier() # feedback baseado na info classica
qc.x(q3).c_if(c1, 0); qc.z(q3).c_if(c0, 0)
qc.barrier() # para verificacao da funcionalidade do cq
qc.measure(q0, c2); qc.measure(q3, c3)
qc.draw(output = 'mpl')

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

# parece que não podemos executar no chip real a c_if
job = qiskit.execute(qc, backend = device, shots = nshots)
job_monitor(job)
plot_histogram(job.result().get_counts(qc))

import math

def E_Phi(p, q):
    return 2*math.sqrt(p*(1-p)*q*(1-q))/(p*q+(1-p)*(1-q))
def E_Psi(p, q):
    return 2*math.sqrt(p*(1-p)*q*(1-q))/(p*(1-q)+(1-p)*q)

import numpy as np
from matplotlib import pyplot as plt
import matplotlib
matplotlib.rcParams.update({'font.size':12})
plt.figure(figsize = (6,4), dpi = 100)

p = np.zeros(101)
Eq0 = np.zeros(101); Eq001 = np.zeros(101); Eq005 = np.zeros(101); Eq025 = np.zeros(101); Eq05 = np.zeros(101)
Eq0 = np.zeros(101); Eq099 = np.zeros(101); Eq095 = np.zeros(101); Eq075 = np.zeros(101)
for j in range(0, 101):
    p[j] = j*0.01
    if p[j] != 0 and p[j] != 1:
        q = 0.0; Eq0[j] = E_Phi(p[j], q); q = 0.01; Eq001[j] = E_Phi(p[j], q)
        q = 0.05; Eq005[j] = E_Phi(p[j], q); q = 0.25; Eq025[j] = E_Phi(p[j],q)
        q = 0.5; Eq05[j] = E_Phi(p[j],q); q = 0.75; Eq075[j] = E_Phi(p[j],q)
        q = 0.95; Eq095[j] = E_Phi(p[j],q); q = 0.99; Eq099[j] = E_Phi(p[j],q)
plt.plot(p, Eq0); plt.plot(p, Eq001); plt.plot(p, Eq005); plt.plot(p, Eq025); plt.plot(p, Eq05)
plt.plot(p, Eq075); plt.plot(p, Eq095); plt.plot(p, Eq099)
plt.xlabel('p'); plt.ylabel(r'$E(\rho_{AD})$')
plt.legend(('q=0','q=0.01','q=0.05','q=0.25','q=0.5','q=0.75','q=0.95','q=0.99'), 
           bbox_to_anchor=(1.05, 1.0), loc='upper left')
plt.title(r'$|\Phi_{\pm}\rangle$')
plt.show()

matplotlib.rcParams.update({'font.size':12})
plt.figure(figsize = (6,4), dpi = 100)
p = np.zeros(101)
Eq0 = np.zeros(101); Eq001 = np.zeros(101); Eq005 = np.zeros(101); Eq025 = np.zeros(101); Eq05 = np.zeros(101)
Eq0 = np.zeros(101); Eq099 = np.zeros(101); Eq095 = np.zeros(101); Eq075 = np.zeros(101)
for j in range(0, 101):
    p[j] = j*0.01
    if p[j] != 0 and p[j] != 1:
        q = 0.0; Eq0[j] = E_Psi(p[j], q); q = 0.01; Eq001[j] = E_Psi(p[j], q)
        q = 0.05; Eq005[j] = E_Psi(p[j], q); q = 0.25; Eq025[j] = E_Psi(p[j],q)
        q = 0.5; Eq05[j] = E_Psi(p[j],q); q = 0.75; Eq075[j] = E_Psi(p[j],q)
        q = 0.95; Eq095[j] = E_Psi(p[j],q); q = 0.99; Eq099[j] = E_Psi(p[j],q)
plt.plot(p, Eq0); plt.plot(p, Eq001); plt.plot(p, Eq005); plt.plot(p, Eq025); plt.plot(p, Eq05)
plt.plot(p, Eq075); plt.plot(p, Eq095); plt.plot(p, Eq099)
plt.xlabel('p'); plt.ylabel(r'$E(\rho_{AD})$')
plt.legend(('q=0','q=0.01','q=0.05','q=0.25','q=0.5','q=0.75','q=0.95','q=0.99'), 
           bbox_to_anchor=(1.05, 1.0), loc='upper left')
plt.title(r'$|\Psi_{\pm}\rangle$')
plt.show()

E_Phi(p,q), E_Psi(p,q)

p,q = symbols('p q'); factor(diff(E_Phi(p,q),p)), factor(diff(E_Psi(p,q),p))

E_Phi(1-q,q), E_Psi(q,q) # nos 2 casos, E=1

import numpy as np
from matplotlib import pyplot as plt

q = np.arange(0.001, 1.0, 0.001); PmaxE = 2*(1-q)*q; plt.plot(q, 1-PmaxE, label = r'$1-Pr(\max E)$'); 
pred = 1 + q*np.log2(q) + (1-q)*np.log2(1-q); plt.plot(q, pred, label = r'$P_{vn}(\rho_{A})$')
plt.legend(); plt.xlabel(r'$q$'); plt.show()

