from sympy import Matrix, sqrt, init_printing; init_printing(use_unicode=True)
ket0 = Matrix([[1],[0]]); ket1 = Matrix([[0],[1]]); from sympy.physics.quantum import TensorProduct as tp
ket00 = tp(ket0,ket0); ket01 = tp(ket0,ket1); ket10 = tp(ket1,ket0); ket11 = tp(ket1,ket1)
Psim = (ket01-ket10)/sqrt(2); Psim.T

X = Matrix([[0,1],[1,0]]); Z = Matrix([[1,0],[0,-1]]); (tp(Z,Z)*Psim).T, (tp(X,X)*Psim).T

A1 = Z; A2 = X; B1 = (Z+X)/sqrt(2); B2 = (Z-X)/sqrt(2); O = tp(A1,B1)+tp(A1,B2)+tp(A2,B1)-tp(A2,B2)
O, sqrt(2)*(tp(Z,Z)+tp(X,X))

from qiskit import QuantumCircuit
def qc_Psim(): # função que cria o circuito quantico para preparar o estado de Bell |Psi->
    qc = QuantumCircuit(2, name=r'$|\Psi_-\rangle$')
    qc.x([0,1])
    qc.h(0)
    qc.cx(0,1)
    #qc.z(0)
    #qc.x(1)
    return qc
qc_Psim_ = qc_Psim()
qc_Psim_.draw('mpl')

import qiskit
from qiskit.quantum_info import Statevector
qc_Psim_ = qc_Psim()
state = Statevector(qc_Psim_)
qiskit.visualization.plot_state_city(state)

from sympy import Matrix
Psim = Matrix([[0],[1],[-1],[0]])/sqrt(2)
#Psim
from sympy.physics.quantum.dagger import Dagger
rho_Psim = Psim*Dagger(Psim)
rho_Psim

qc = QuantumCircuit(2,2)
qc_Psim_ = qc_Psim()
qc.append(qc_Psim_,[0,1])
qc.measure([0,1],[0,1])
qc.draw('mpl')

from qiskit import Aer, execute
simulator = Aer.get_backend('qasm_simulator')
nshots = 2**13
job = execute(qc, backend=simulator, shots=nshots)
counts = job.result().get_counts()
p00=0; p01=0; p10=0; p11=0
for j in (0,2):
    for k in (0,2):
        if '00' in counts:
            p00 = counts['00']/nshots
        if '01' in counts:
            p01 = counts['01']/nshots
        if '10' in counts:
            p10 = counts['10']/nshots
        if '11' in counts:
            p11 = counts['11']/nshots
ZZ_avg = (p00+p11)-(p01+p10)
ZZ_avg

qc = QuantumCircuit(2,2)
qc_Psim_ = qc_Psim()
qc.append(qc_Psim_,[0,1])
qc.h([0,1])
qc.measure([0,1],[0,1])
qc.draw('mpl')

from qiskit import Aer
simulator = Aer.get_backend('qasm_simulator')
nshots = 2**13
job = execute(qc, backend=simulator, shots=nshots)
counts = job.result().get_counts()
p00=0; p01=0; p10=0; p11=0
for j in (0,2):
    for k in (0,2):
        if '00' in counts:
            p00 = counts['00']/nshots
        if '01' in counts:
            p01 = counts['01']/nshots
        if '10' in counts:
            p10 = counts['10']/nshots
        if '11' in counts:
            p11 = counts['11']/nshots
XX_avg = (p00+p11)-(p01+p10)
XX_avg

O_avg = sqrt(2)*(ZZ_avg + XX_avg)
O_avg, float(O_avg)

import qiskit
qiskit.IBMQ.save_account('5a0140fc106f9f8ae43104d5751baf65c726195f8b45b1f42e161f0182b5fe79ca6dbb1405159301e1e44d40d680223ea572d9b8737205bfa6f90485f1f88123', 
                  overwrite = True)
qiskit.IBMQ.load_account()
provider = qiskit.IBMQ.get_provider(hub='ibm-q', group='open', project='main')

device = provider.get_backend('ibm_nairobi')
qc = QuantumCircuit(2,2); qc_Psim_ = qc_Psim(); qc.append(qc_Psim_,[0,1]); qc.measure([0,1],[0,1])
job = execute(qc, backend = device, shots = nshots); print(job.job_id())

device = provider.get_backend('ibm_nairobi'); job = device.retrieve_job('cmfbdw78vz2g008daj90')
counts = job.result().get_counts()
p00=0; p01=0; p10=0; p11=0
for j in (0,2):
    for k in (0,2):
        if '00' in counts:
            p00 = counts['00']/nshots
        if '01' in counts:
            p01 = counts['01']/nshots
        if '10' in counts:
            p10 = counts['10']/nshots
        if '11' in counts:
            p11 = counts['11']/nshots
ZZ_avg_exp = (p00+p11)-(p01+p10); ZZ_avg_exp

device = provider.get_backend('ibm_nairobi')
qc = QuantumCircuit(2,2); qc_Psim_ = qc_Psim(); qc.append(qc_Psim_,[0,1]); qc.h([0,1]); qc.measure([0,1],[0,1])
job = execute(qc, backend = device, shots = nshots); jobid = job.job_id()
print(jobid)

device = provider.get_backend('ibm_nairobi'); job = device.retrieve_job('cmfbekafysfg008t2y4g')
counts = job.result().get_counts()
p00=0; p01=0; p10=0; p11=0
for j in (0,2):
    for k in (0,2):
        if '00' in counts:
            p00 = counts['00']/nshots
        if '01' in counts:
            p01 = counts['01']/nshots
        if '10' in counts:
            p10 = counts['10']/nshots
        if '11' in counts:
            p11 = counts['11']/nshots
XX_avg_exp = (p00+p11)-(p01+p10); XX_avg_exp

O_avg_exp = sqrt(2.)*(ZZ_avg_exp + XX_avg_exp)
O_avg_exp

import numpy as np
import math
th_max = math.pi/2; npt = 100; dth = th_max/npt; th = np.arange(0,th_max+dth,dth)
Oteo = -math.sqrt(2)*(1+np.sin(th))
Odl = -2*np.ones(len(th))
from matplotlib import pyplot as plt
plt.plot(th, Oteo, label = r'$\langle O\rangle_{\Psi}^{teo}$')
plt.plot(th, Odl, label = 'LHV limit')
plt.xlabel(r'$\theta$')
plt.legend(); plt.show()

# circuito quântico para preparação do estado Psi
def qc_Psi(th):
    qc = QuantumCircuit(2)
    qc.u(th,0,0, 0)
    qc.cx(0,1)
    qc.z(0)
    qc.x(1)
    return qc
import math
th = math.pi/4
qc_Psi_ = qc_Psi(th)
qc_Psi_.draw('mpl')

# função para calcular as médias tipo <A\otimes B>, com A,B=+1,-1, dadas as probabilidades
def ABavg(result, nshots):
    avg = 0
    if '00' in result:
        avg += result['00']
    if '01' in result:
        avg -= result['01']
    if '10' in result:
        avg -= result['10']
    if '11' in result:
        avg += result['11']
    return avg/nshots # = (N(00)-N(01)-N(10)+N(11))/nshots

th_max = math.pi/2; npe = 5; dth = th_max/npe; th = np.arange(0,th_max+dth,dth)
Osim = np.zeros(len(th))
from qiskit import Aer, execute
simulator = Aer.get_backend('qasm_simulator')
nshots = 2**13
for j in range(0,len(th)):
    # mede ZZ 
    qc = QuantumCircuit(2,2)
    qc_Psi_ = qc_Psi(th[j]); qc.append(qc_Psi_, [0,1])
    qc.measure([0,1],[0,1])
    job = execute(qc, backend=simulator, shots=nshots)
    result = job.result().get_counts()
    ZZavg = ABavg(result, nshots)
    # mede XX 
    qc = QuantumCircuit(2,2)
    qc_Psi_ = qc_Psi(th[j]); qc.append(qc_Psi_, [0,1])
    qc.h([0,1])
    qc.measure([0,1],[0,1])
    job = execute(qc, backend=simulator, shots=nshots)
    result = job.result().get_counts()
    XXavg = ABavg(result, nshots)
    # média de O
    Osim[j] = math.sqrt(2)*(ZZavg + XXavg)
Osim

import numpy as np; import math
th_max = math.pi/2; dtht = th_max/npt; tht = np.arange(0,th_max+dtht,dtht)
Oteo = -math.sqrt(2)*(1+np.sin(tht)); Odl = -2*np.ones(len(tht))
dth = th_max/npe; th = np.arange(0,th_max+dth,dth)
from matplotlib import pyplot as plt
plt.plot(tht, Oteo, label = r'$\langle O\rangle_{\Psi}^{teo}$')
plt.plot(tht, Odl, label = 'LHV limit')
plt.plot(th, Osim, '*', label = r'$\langle O\rangle_{\Psi}^{sim}$')
plt.xlabel(r'$\theta$')
plt.legend(); plt.show()

import qiskit
qiskit.IBMQ.save_account('45b1db3a3c2ecae609e6c98184250919ccd31c9ea24100440abe0821fe6312c2b05a67239806c5b0020754eb3072f924547adb621e1e3931200aee13583a776f', 
                  overwrite = True)
qiskit.IBMQ.load_account()
provider = qiskit.IBMQ.get_provider(hub='ibm-q', group='open', project='main')

th_max = math.pi/2; dth = th_max/npe; th = np.arange(0,th_max+dth,dth)
device = provider.get_backend('ibm_nairobi')
nshots = 2**13
jobs_ids_ZZ = []; jobs_ids_XX = []
for j in range(0,len(th)):
    # mede ZZ 
    qc = QuantumCircuit(2,2); qc_Psi_ = qc_Psi(th[j]); qc.append(qc_Psi_, [0,1]); qc.measure([0,1],[0,1])
    job = execute(qc, backend=device, shots=nshots)
    print('ZZ',job.job_id()); jobs_ids_ZZ.append(job.job_id())
    # mede XX 
    qc = QuantumCircuit(2,2); qc_Psi_ = qc_Psi(th[j]); qc.append(qc_Psi_, [0,1]); qc.h([0,1]); qc.measure([0,1],[0,1])
    job = execute(qc, backend=device, shots=nshots)
    print('XX',job.job_id()); jobs_ids_XX.append(job.job_id())
f = open("jobs_ids_ZZ.txt", "w"); f.write(str(jobs_ids_ZZ)); f.close()
f = open("jobs_ids_XX.txt", "w"); f.write(str(jobs_ids_XX)); f.close()

f = open("jobs_ids_ZZ.txt","r")
list_ids_ZZ = f.read().replace("'","").replace(" ","").replace("[","").replace("]","").split(",")
f.close()
f = open("jobs_ids_XX.txt","r")
list_ids_XX = f.read().replace("'","").replace(" ","").replace("[","").replace("]","").split(",")
f.close()
th_max = math.pi/2; dth = th_max/npe; th = np.arange(0,th_max+dth,dth)
Oexp = np.zeros(len(th))
for j in range(0,len(th)):
    # média de ZZ
    job = device.retrieve_job(list_ids_ZZ[j])
    result = job.result().get_counts()
    ZZavg = ABavg(result, nshots)
    # média de XX
    job = device.retrieve_job(list_ids_XX[j])
    result = job.result().get_counts()
    XXavg = ABavg(result, nshots)
    # média de O
    Oexp[j] = math.sqrt(2)*(ZZavg + XXavg)
Oexp

import numpy as np; import math
th_max = math.pi/2; dtht = th_max/npt; tht = np.arange(0,th_max+dtht,dtht)
Oteo = -math.sqrt(2)*(1+np.sin(tht)); Odl = -2*np.ones(len(tht))
dth = th_max/npe; th = np.arange(0,th_max+dth,dth)
from matplotlib import pyplot as plt
plt.plot(tht, Oteo, label = r'$\langle O\rangle_{\Psi}^{teo}$')
plt.plot(tht, Odl, label = 'LHV limit')
plt.plot(th, Osim, '*', label = r'$\langle O\rangle_{\Psi}^{sim}$')
plt.plot(th, Oexp, '*', label = r'$\langle O\rangle_{\Psi}^{exp}$')
plt.xlabel(r'$\theta$'); plt.legend(); plt.show()

def f_Oavg(O,psi):
    from sympy.physics.quantum.dagger import Dagger
    return Dagger(psi)*O*psi

from sympy import Matrix, symbols, cos, sin, sqrt, simplify
from sympy.physics.quantum import TensorProduct as tp
Psim = Matrix([[0],[1],[-1],[0]])/sqrt(2)
X = Matrix([[0,1],[1,0]])
Z = Matrix([[1,0],[0,-1]])
th = symbols('theta')
A1B1 = tp(Z,cos(th)*Z+sin(th)*X)
A1B1avg = f_Oavg(A1B1,Psim)
A1B2 = tp(Z,(Z-X)/sqrt(2))
A1B2avg = f_Oavg(A1B2,Psim)
A2B1 = tp(X,cos(th)*Z+sin(th)*X)
A2B1avg = f_Oavg(A2B1,Psim)
A2B2 = tp(X,(Z-X)/sqrt(2))
A2B2avg = f_Oavg(A2B2,Psim)
Oavg = A1B1avg + A1B2avg + A2B1avg - A2B2avg
A1B1avg, A1B2avg, A2B1avg, A2B2avg, simplify(Oavg)

import numpy as np; import math
th_max = math.pi/4; npt = 100; dth = th_max/npt; tht = np.arange(-th_max,th_max+dth,dth)
Oteo = -math.sqrt(2)*(1+np.sin(tht+math.pi/4))
Odl = -2*np.ones(len(tht))
from matplotlib import pyplot as plt
plt.plot(tht, Oteo, label = r'$\langle O\rangle_{\Psi}^{teo}$')
plt.plot(tht, Odl, label = 'LHV limit')
plt.xlabel(r'$\theta$')
plt.legend(); plt.show()

th_max = math.pi/4; npe = 4; dth = 2*th_max/npe; th = np.arange(-th_max,th_max+dth,dth); Osim = np.zeros(len(th))
from qiskit import *; simulator = Aer.get_backend('qasm_simulator'); nshots = 2**13
for j in range(0,len(th)):
    # mede <A1B1> 
    qc = QuantumCircuit(2,2); qc_Psim_ = qc_Psim(); qc.append(qc_Psim_, [0,1]); qc.u(th[j],0,math.pi, 1); qc.measure([0,1],[0,1])
    job = execute(qc, backend=simulator, shots=nshots); result = job.result().get_counts(); A1B1avg = ABavg(result, nshots)
    # mede <A1B2>
    qc = QuantumCircuit(2,2); qc_Psim_ = qc_Psim(); qc.append(qc_Psim_, [0,1]); qc.u(-math.pi/4,0,math.pi, 1); qc.measure([0,1],[0,1])
    job = execute(qc, backend=simulator, shots=nshots); result = job.result().get_counts(); A1B2avg = ABavg(result, nshots)
    # mede <A2B1> 
    qc = QuantumCircuit(2,2); qc_Psim_ = qc_Psim(); qc.append(qc_Psim_, [0,1]); qc.h(0); qc.u(th[j],0,math.pi, 1); qc.measure([0,1],[0,1])
    job = execute(qc, backend=simulator, shots=nshots); result = job.result().get_counts(); A2B1avg = ABavg(result, nshots)
    # mede <A2B2>
    qc = QuantumCircuit(2,2); qc_Psim_ = qc_Psim(); qc.append(qc_Psim_, [0,1]); qc.h(0); qc.u(-math.pi/4,0,math.pi, 1); qc.measure([0,1],[0,1])
    job = execute(qc, backend=simulator, shots=nshots); result = job.result().get_counts(); A2B2avg = ABavg(result, nshots)
    # média de O
    Osim[j] = A1B1avg + A1B2avg + A2B1avg - A2B2avg
Osim

import numpy as np; import math
th_max = math.pi/4; dtht = 2*th_max/npt; tht = np.arange(-th_max,th_max+dtht,dtht)
Oteo = -math.sqrt(2)*(1+np.sin(tht+math.pi/4)); Odl = -2*np.ones(len(tht))
dth = 2*th_max/npe; th = np.arange(-th_max,th_max+dth,dth)
from matplotlib import pyplot as plt
plt.plot(tht, Oteo, label = r'$\langle O\rangle_{\Psi}^{teo}$')
plt.plot(tht, Odl, label = 'LHV limit')
plt.plot(th, Osim, '*', label = r'$\langle O\rangle_{\Psi}^{sim}$')
plt.xlabel(r'$\theta$')
plt.legend(); plt.show()

import qiskit
qiskit.IBMQ.save_account('45b1db3a3c2ecae609e6c98184250919ccd31c9ea24100440abe0821fe6312c2b05a67239806c5b0020754eb3072f924547adb621e1e3931200aee13583a776f', 
                  overwrite = True)
qiskit.IBMQ.load_account()
#provider = qiskit.IBMQ.get_provider(hub='ibm-q', group='open', project='main')
provider = qiskit.IBMQ.get_provider(hub='ibm-q-research-2', group='federal-uni-sant-1', project='main')

th_max = math.pi/4; npe = 4; dth = 2*th_max/npe; th = np.arange(-th_max,th_max+dth,dth)
from qiskit import *; device = provider.get_backend('ibm_nairobi'); nshots = 2**13
jobs_ids_A1B1 = []; jobs_ids_A1B2 = []; jobs_ids_A2B1 = []; jobs_ids_A2B2 = []
for j in range(0,len(th)):
    # mede <A1B1> 
    qc = QuantumCircuit(2,2); qc_Psim_ = qc_Psim(); qc.append(qc_Psim_, [0,1]); qc.u(th[j],0,math.pi, 1)
    qc.measure([0,1],[0,1])
    job = execute(qc, backend=device, shots=nshots); print('A1B1:',job.job_id()); jobs_ids_A1B1.append(job.job_id())
    # mede <A1B2>
    qc = QuantumCircuit(2,2); qc_Psim_ = qc_Psim(); qc.append(qc_Psim_, [0,1]); qc.u(-math.pi/4,0,math.pi, 1)
    qc.measure([0,1],[0,1])
    job = execute(qc, backend=device, shots=nshots); print('A1B2:',job.job_id()); jobs_ids_A1B2.append(job.job_id())
    # mede <A2B1> 
    qc = QuantumCircuit(2,2); qc_Psim_ = qc_Psim(); qc.append(qc_Psim_, [0,1]); qc.h(0); qc.u(th[j],0,math.pi, 1)
    qc.measure([0,1],[0,1])
    job = execute(qc, backend=device, shots=nshots); print('A2B1:',job.job_id()); jobs_ids_A2B1.append(job.job_id())
    # mede <A2B2>
    qc = QuantumCircuit(2,2); qc_Psim_ = qc_Psim(); qc.append(qc_Psim_, [0,1]); qc.h(0); qc.u(-math.pi/4,0,math.pi, 1)
    qc.measure([0,1],[0,1])
    job = execute(qc, backend=device, shots=nshots); print('A2B2:',job.job_id()); jobs_ids_A2B2.append(job.job_id())
f = open("jobs_ids_A1B1.txt", "w"); f.write(str(jobs_ids_A1B1)); f.close()
f = open("jobs_ids_A1B2.txt", "w"); f.write(str(jobs_ids_A1B2)); f.close()
f = open("jobs_ids_A2B1.txt", "w"); f.write(str(jobs_ids_A2B1)); f.close()
f = open("jobs_ids_A2B2.txt", "w"); f.write(str(jobs_ids_A2B2)); f.close()

device = provider.get_backend('ibm_nairobi')
f = open("jobs_ids_A1B1.txt","r")
list_ids_A1B1 = f.read().replace("'","").replace(" ","").replace("[","").replace("]","").split(","); f.close()
f = open("jobs_ids_A1B2.txt","r")
list_ids_A1B2 = f.read().replace("'","").replace(" ","").replace("[","").replace("]","").split(","); f.close()
f = open("jobs_ids_A2B1.txt","r")
list_ids_A2B1 = f.read().replace("'","").replace(" ","").replace("[","").replace("]","").split(","); f.close()
f = open("jobs_ids_A2B2.txt","r")
list_ids_A2B2 = f.read().replace("'","").replace(" ","").replace("[","").replace("]","").split(","); f.close()
th_max = math.pi/4; dth = 2*th_max/npe; th = np.arange(-th_max,th_max+dth,dth); Oexp = np.zeros(len(th))
for j in range(0,len(th)):
    # média de A1B1
    job = device.retrieve_job(list_ids_A1B1[j]); result = job.result().get_counts(); A1B1avg = ABavg(result, nshots)
    # média de A1B2
    job = device.retrieve_job(list_ids_A1B2[j]); result = job.result().get_counts(); A1B2avg = ABavg(result, nshots)
    # média de A2B1
    job = device.retrieve_job(list_ids_A2B1[j]); result = job.result().get_counts(); A2B1avg = ABavg(result, nshots)
    # média de A2B2
    job = device.retrieve_job(list_ids_A2B2[j]); result = job.result().get_counts(); A2B2avg = ABavg(result, nshots)
    # média de O
    Oexp[j] = A1B1avg + A1B2avg + A2B1avg - A2B2avg
Oexp

import numpy as np; import math
th_max = math.pi/4; dtht = 2*th_max/npt; tht = np.arange(-th_max,th_max+dtht,dtht)
Oteo = -math.sqrt(2)*(1+np.sin(tht+math.pi/4)); Odl = -2*np.ones(len(tht))
dth = 2*th_max/npe; th = np.arange(-th_max,th_max+dth,dth)
from matplotlib import pyplot as plt
plt.plot(tht, Oteo, label = r'$\langle O\rangle_{\Psi}^{teo}$')
plt.plot(tht, Odl, label = 'LHV limit')
plt.plot(th, Osim, '*', label = r'$\langle O\rangle_{\Psi}^{sim}$')
plt.plot(th, Oexp, '*', label = r'$\langle O\rangle_{\Psi}^{exp}$')
plt.xlabel(r'$\theta$'); plt.legend(); plt.show()

