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

cos(pi/8), sin(pi/8)

U = (1/sqrt(2))*Matrix([[1,1],[1,-1]]); V1 = Matrix([[1,0],[0,1]]); V2 = V1
D1 = Matrix([[cos(pi/8),0],[0,cos(pi/8)]]); D2 = Matrix([[sin(pi/8),0],[0,sin(pi/8)]])
U, U*U.T, D1, D2

M1 = V1*D1*U; M2 = V2*D2*U; M1, M2

M1*M1, M2*M2, M1*M1 + M2*M2 # não são projetores, mas são operadores de medida

2*float((1/4)*(1+1/sqrt(2))), 2*float((1/4)*(1-1/sqrt(2)))

qr = QuantumRegister(2); cr = ClassicalRegister(1); qc = QuantumCircuit(qr, cr)
qc.h(qr[0]); qc.x(qr[0]); qc.cry(math.pi/4, 0, 1); qc.x(qr[0]); qc.cry(math.pi/4, 0, 1)
# V1=V2=I
qc.measure(1,0)
qc.draw(output = 'mpl')

job_sim = execute(qc, backend = simulator, shots = nshots)
job_exp = execute(qc, backend = device, shots = nshots)
job_monitor(job_exp)

plot_histogram([job_sim.result().get_counts(), job_exp.result().get_counts()], 
                legend = ['sim', 'exp'])

def U(xi, et, ga):
    return Matrix([[(cos(xi)+1j*sin(xi))*cos(ga), (cos(et)+1j*sin(et))*sin(ga)],
                    [-(cos(et)-1j*sin(et))*sin(ga), (cos(xi)-1j*sin(xi))*cos(ga)]])
xi, et, ga = symbols('xi eta gamma'); U(xi, et, ga)

def D1(th, ph):
    return Matrix([[cos(th/2),0],[0,cos(ph/2)]])
def D2(th, ph):
    return Matrix([[sin(th/2),0],[0,sin(ph/2)]])
th,ph = symbols('theta phi'); D1(th, ph), D2(th, ph)

xi1,et1,ga1,xi2,et2,ga2,xi3,et3,ga3 = symbols('xi_1 eta_1 gamma_1 xi_2 eta_2 gamma_2 xi_3 eta_3 gamma_3')
simplify(U(xi1,et1,ga1)*D1(th,ph)*U(xi2,et2,ga2) - (1/2)*sqrt(1+1/sqrt(2))*Matrix([[1,1],[1,-1]]))

simplify(U(xi3,et3,ga3)*D2(th,ph)*U(xi2,et2,ga2) - (1/2)*sqrt(1+1/sqrt(2))*Matrix([[1,1],[1,-1]]))

cos(pi/3), sin(pi/3), cos(2*pi/3), sin(2*pi/3)

M1 = sqrt(2/3)*Matrix([[1,0],[0,0]])
M2 = (1/(4*sqrt(6)))*Matrix([[1,sqrt(3)],[sqrt(3),3]])
M3 = (1/(4*sqrt(6)))*Matrix([[1,-sqrt(3)],[-sqrt(3),3]])
M1,M2,M3

M1.T*M1 + M2.T*M2 + M3.T*M3

M1 = sqrt(2/3.)*Matrix([[1,0],[0,0]]); M2 = sqrt(2/3.)*(1/4)*Matrix([[1,sqrt(3.)],[sqrt(3.),3]])
M3 = sqrt(2/3.)*(1/4)*Matrix([[1,-sqrt(3.)],[-sqrt(3.),3]])
M1, M2, M3

M1.T*M1 + M2.T*M2 + M3.T*M3

2/3, 1/6, 2/3+2*(1/6)

#th1 = acos(sqrt(2/3)); ph1 = pi; th2 = pi/2; ph2 = pi/2
D11 = Matrix([[sqrt(2/3),0],[0,0]])
D21 = Matrix([[sqrt(1/3),0],[0,1]]) 
D12 = Matrix([[1/sqrt(2),0],[0,1/sqrt(2)]])
D22 = Matrix([[1/sqrt(2),0],[0,1/sqrt(2)]])
U = Matrix([[1,0],[0,1]])
V11 = Matrix([[1,0],[0,1]])
V21 = (1/sqrt(2))*Matrix([[1,1],[-1,1]])
V12 = (1/2)*Matrix([[1,-sqrt(3)],[sqrt(3),1]])
V22 = -(1/2)*Matrix([[sqrt(3),-1],[1,sqrt(3)]])

M1 = V11*D11*U
np.array(M1).astype(np.float64), np.array(Matrix([[sqrt(2/3),0],[0,0]])).astype(np.float64)

M2 = V12*D12*V21*D21*U
np.array(M2).astype(np.float64), np.array((1/4)*sqrt(2/3)*Matrix([[1,sqrt(3)],[sqrt(3),3]])).astype(np.float64)
# não é o resultado que precisamos

M3 = V22*D22*V21*D21*U
np.array(M3).astype(np.float64), np.array((1/4)*sqrt(2/3)*Matrix([[1,-sqrt(3)],[-sqrt(3),3]])).astype(np.float64)
# não é o resultado que precisamos

np.array(M1.T*M1 + M2.T*M2 + M3.T*M3).astype(np.float64) # esta ok a relacao de completeza

cos(pi/3), sin(pi/3), cos(pi/6)

def qc_ry(th):
    qr = QuantumRegister(1); qc = QuantumCircuit(qr, name = 'RY')
    qc.ry(th, 0)
    return qc
def qc_u(th,ph,lb):
    qr = QuantumRegister(1); qc = QuantumCircuit(qr, name = 'U')
    qc.u(th,ph,lb, 0)
    return qc

qr = QuantumRegister(3); cr = ClassicalRegister(2); qc = QuantumCircuit(qr, cr)
# U = I
qc.x(qr[0]); qc.cry(math.acos(math.sqrt(2/3)), 0, 1); qc.x(qr[0]); qc.cry(math.pi/2, 0, 1)
# V11 = I
qc.cu(math.pi/2,math.pi,math.pi,0, 1,0)
qc.barrier()
qc_ry_ = qc_ry(0); ccry = qc_ry_.to_gate().control(2) # cria a ctrl-ctrl-RY
qc.x(0); qc.append(ccry, [0,1,2]); qc.x(0)
qc_ry_ = qc_ry(math.pi/2); ccry = qc_ry_.to_gate().control(2)
qc.append(ccry, [0,1,2]) # os primeiros sao o controle, os ultimos sao o target
qc.x(2)
qc_u_ = qc_u(2*math.pi/3,0,0); ccu = qc_u_.to_gate().control(2)
qc.append(ccu, [1,2,0])
qc.x(2)
qc_u_ = qc_u(4*math.pi/3,0,0); ccu = qc_u_.to_gate().control(2)
qc.append(ccu, [1,2,0])
qc.barrier()
qc.measure(1,0); qc.measure(2,1)
qc.draw(output = 'mpl')

qc.decompose().draw(output = 'mpl')

job_sim = execute(qc, backend = simulator, shots = nshots)
job_exp = execute(qc, backend = device, shots = nshots)
job_monitor(job_exp)

plot_histogram([job_sim.result().get_counts(), job_exp.result().get_counts()], 
                legend = ['sim', 'exp'])

# sequencia: 000 001 010 011 100 101 110 111 = 0 1 2 3 4 5 6 7
Phi_ACD = [math.sqrt(2/3), (1/4)*math.sqrt(2/3), (1/4)*math.sqrt(2/3), 0,
           0, (1/4)*math.sqrt(2), -(1/4)*math.sqrt(2), 0]
qr = QuantumRegister(3); cr = ClassicalRegister(2); qc = QuantumCircuit(qr,cr)
qc.initialize(Phi_ACD, [qr[2],qr[1],qr[0]])
qc.draw(output='mpl')

backend = BasicAer.get_backend('statevector_simulator')
job = backend.run(transpile(qc, backend))
qc_state = job.result().get_statevector(qc)
qc_state

sqrt(2/3), 0.25*sqrt(2/3), 0.25*sqrt(2.)

state_fidelity(Phi_ACD,qc_state)

qc.measure([1,2],[0,1])
qc.draw(output='mpl')

job_sim = execute(qc, backend = simulator, shots = nshots)
job_exp = execute(qc, backend = device, shots = nshots)
job_monitor(job_exp)

# tava dando errado o resultado. Me parece que na inicializacao a ordem dos qubits foi trocada
plot_histogram([job_sim.result().get_counts(), job_exp.result().get_counts()], 
                legend = ['sim', 'exp'])

# sequencia: 000 001 010 011 100 101 110 111 = 0 1 2 3 4 5 6 7
Phi_ACD = [1/math.sqrt(2), 1/(2*math.sqrt(2)), 0, 1/(2*math.sqrt(2)), 
           0, 1/(2*math.sqrt(2)), 0, -1/(2*math.sqrt(2))]
qr = QuantumRegister(3); cr = ClassicalRegister(2); qc = QuantumCircuit(qr,cr)
qc.initialize(Phi_ACD, [qr[2],qr[1],qr[0]])
qc.measure([1,2],[0,1])
qc.draw(output='mpl')

qc.decompose().decompose().decompose().decompose().decompose().draw(output='mpl')

job_sim = execute(qc, backend = simulator, shots = nshots)
job_exp = execute(qc, backend = device, shots = nshots)
job_monitor(job_exp)

plot_histogram([job_sim.result().get_counts(), job_exp.result().get_counts()], 
                legend = ['sim', 'exp'])

qr = QuantumRegister(4); cr = ClassicalRegister(3); qc = QuantumCircuit(qr, cr)
# U = I
qc.x(qr[0]); qc.cry(math.acos(math.sqrt(2/3)), 0, 1); qc.x(qr[0]); qc.cry(math.pi/2, 0, 1)
# V11 = I
qc.cu(math.pi/2,math.pi,math.pi,0, 1,0)
qc.barrier()
qc_ry_ = qc_ry(0); ccry = qc_ry_.to_gate().control(2) # cria a ctrl-ctrl-RY
qc.x(0); qc.append(ccry, [0,1,2]); qc.x(0)
qc_ry_ = qc_ry(math.pi/2); ccry = qc_ry_.to_gate().control(2)
qc.append(ccry, [0,1,2]) # os primeiros sao o controle, os ultimos sao o target
qc.x(2)
qc_u_ = qc_u(2*math.pi/3,0,0); ccu = qc_u_.to_gate().control(2)
qc.append(ccu, [1,2,0])
qc.x(2)
qc_u_ = qc_u(4*math.pi/3,0,0); ccu = qc_u_.to_gate().control(2)
qc.append(ccu, [1,2,0])
qc.barrier()
qc_ry_ = qc_ry(math.pi); cccry = qc_ry_.to_gate().control(3) # cria a ctrl-ctrl-ctrl-RY
qc.x(0); qc.append(cccry, [0,1,2,3]); qc.x(0)
qc_ry_ = qc_ry(math.pi/2); cccry = qc_ry_.to_gate().control(3)
qc.append(cccry, [0,1,2,3])
qc_u_ = qc_u(math.pi,0,0); cccu = qc_u_.to_gate().control(3)
qc.x(3); qc.append(cccu, [3,2,1,0]); qc.x(3)
qc.append(cccu, [3,2,1,0])
qc_u_ = qc_u(math.pi/3,0,0); cccu = qc_u_.to_gate().control(3)
qc.barrier()
qc.measure(1,0); qc.measure(2,1); qc.measure(3,2)
qc.draw(output = 'mpl')

qc.decompose().draw(output = 'mpl')



%run init.ipynb

# 2 element 1-qubit povm (Yordanov)
M1 = (1/2)*Matrix([[1,0],[0,sqrt(3)]]); M2 = M1
M1, M2

M1.T*M1 + M2.T*M2

# 3 element 1-qubit povm (Yordanov)
M1 = sqrt(2/3)*Matrix([[1,0],[0,0]])
M2 = (1/(4*sqrt(6)))*Matrix([[1,sqrt(3)],[sqrt(3),1]])
M3 = (1/(4*sqrt(6)))*Matrix([[1,-sqrt(3)],[-sqrt(3),1]])
M1,M2,M3

M1.T*M1 + M2.T*M2 + M3.T*M3

# 2-element 1-qubit POVM (Douglas)
M1 = (1/(2*sqrt(2)))*Matrix([[1,0],[sqrt(3),2]])
M2 = (1/(2*sqrt(2)))*Matrix([[1,0],[-sqrt(3),2]])
M1,M2

M1.T*M1 + M2.T*M2

