#import sys
#import time
import numpy as np
import matplotlib.pyplot as plt
from qiskit_ibm_runtime import QiskitRuntimeService, Sampler
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator



#service = QiskitRuntimeService("ibm_quantum",token=tk)
#device_name="ibm_kyiv"
#backend= service.backend(device_name)
qc=[]
for i in range(4):
    qc.append(QuantumCircuit(2,2))
    qc[-1].sx(0)
    qc[-1].sx(1)
    qc[-1].s(0)
    qc[-1].s(1)
    qc[-1].cz(0,1)
    qc[-1].sx(0)
    qc[-1].sx(1)
    #qc[-1].save_statevector()
    qc[-1].rz((4*(i%2)-3)*np.pi/8,0)
    qc[-1].rz((4*(i//2)-3)*np.pi/8,1)
    qc[-1].sx(0)
    qc[-1].sx(1)
    qc[-1].measure([0,1],[0,1])
#fig=qc.draw("mpl")
#fig.savefig(".pdf")
backend = AerSimulator()
sampler = Sampler(backend)
n=1000
job=sampler.run(qc, shots=n)
print(job.job_id())
result=job.result()
#v=np.asarray(result)
#for i in range(4):
#    print(i,v[i])
b=[]
for r in result:
    a=r.data.c.get_counts()
    print(a)
    b.append((a["00"]+a["11"]-a["01"]-a["10"])/n)
    
print(*b)
print("Bell = " ,b[0]+b[1]+b[2]-b[3])
print("error = ",np.sqrt((4-b[0]**2-b[1]**2-b[2]**2-b[3]**2)/n))