from sympy import * 
init_printing(use_unicode=True)

mu_u = Matrix([[-1],[+1],[+1]]); mu_u

A_R = Matrix([[1,1,-1,-1],[1,-1,1,-1],[1,-1,-1,1]]); A_R

p1,p2,p3,p4 = symbols('p_1 p_2 p_3 p_4')
p4 = 1-p1-p2-p3
cl = A_R[:,0]*p1 + A_R[:,1]*p2 + A_R[:,2]*p3 + A_R[:,3]*p4
mu_u - cl

s1 = Matrix([[0,1],[1,0]])
s2 = Matrix([[0,-1j],[1j,0]])
s1*s2 

import numpy as np

a = [] # gera uma lista com com todos os 512 estados ônticos
for j in range(-1,2,2):
    for k in range(-1,2,2):
        for l in range(-1,2,2):
            for m in range(-1,2,2):
                for n in range(-1,2,2):
                    for o in range(-1,2,2):
                        for p in range(-1,2,2):
                            for q in range(-1,2,2):
                                for r in range(-1,2,2):
                                    a.append(np.array([j,k,l,m,n,o,p,q,r]))

a[10], a[10][5], len(a)

mu = [] # gera, a partir dos estados ônticos, uma lista com os 512 vetores de medida, muitos dos quais são iguais
for j in range(0,2**9):
    r1 = a[j][0]*a[j][1]*a[j][2]
    r2 = a[j][3]*a[j][4]*a[j][5]
    r3 = a[j][6]*a[j][7]*a[j][8]
    c1 = a[j][0]*a[j][3]*a[j][6]
    c2 = a[j][1]*a[j][4]*a[j][7]
    c3 = a[j][2]*a[j][5]*a[j][8]
    mu.append([r1,r2,r3,c1,c2,c3])
mu[0], len(mu)

mu_ = [] # remove as repeticoes do vetor de medida
for j in range(0,2**6):
    if mu[j] not in mu_:
        mu_.append(mu[j])
mu_, len(mu_), mu_[1]

A_R = np.zeros((6,16))#; A_R
for j in range(0,16):
    A_R[:,j] = mu_[j]
print(A_R)
print(A_R.shape)

def rpv_zhsl(d): # vetor de probabilidades aleatório
  rn = np.zeros(d-1)
  for j in range(0,d-1):
    rn[j] = np.random.random()
  rpv = np.zeros(d)
  rpv[0] = 1.0 - rn[0]**(1.0/(d-1.0))
  norm = rpv[0]
  if d > 2:
    for j in range(1,d-1):        
      rpv[j] = (1.0 - rn[j]**(1.0/(d-j-1)))*(1.0-norm)
      norm = norm + rpv[j]
  rpv[d-1] = 1.0 - norm
  return rpv
rpv = rpv_zhsl(16)
print(rpv)

print(np.linalg.norm(rpv))

mu_Q = np.zeros((6,1)); mu_Q = np.array([1,1,1,1,1,-1]); print(mu_Q.shape)
p = np.zeros((16,1)); p =  rpv_zhsl(16); print(p.shape)
vec = mu_Q - A_R@p
print(vec)

def objective(x):
    mu_Q = np.array([1,1,1,1,1,-1])
    vec = mu_Q - A_R@x
    return np.linalg.norm(vec)

from scipy.optimize import minimize

# define o vinculo como sendo uma distribuicao de probabilidades
constraints = [{'type': 'eq', 'fun': lambda x: np.sum(x)-1},
               {'type':'ineq', 'fun': lambda x: x},
               {'type':'ineq', 'fun': lambda x: 1-x}]
# 'eq' quer dizer = 0 e 'ineq' quer dizer >=0

np.random.seed(130000)
x0 = rpv_zhsl(16)
print(x0)

result = minimize(objective, x0, constraints=constraints, method='trust-constr')
print('melhor distribuicao de probabilidades', result.x)

print(objective(result.x))

from qiskit import *
import numpy as np
import math
import qiskit
nshots = 8192
IBMQ.load_account()
provider = qiskit.IBMQ.get_provider(hub='ibm-q', group='open', project='main')
device = provider.get_backend('ibmq_manila')
simulator = Aer.get_backend('qasm_simulator')
from qiskit.tools.monitor import job_monitor
#from qiskit.ignis.verification.tomography import state_tomography_circuits, StateTomographyFitter
#from qiskit.ignis.mitigation.measurement import complete_meas_cal, CompleteMeasFitter
from qiskit.visualization import plot_histogram

qr = QuantumRegister(3)
cr = ClassicalRegister(1)
qc = QuantumCircuit(qr,cr)
qc.barrier() # mede XI
qc.h(0)
qc.cx(0,2) 
qc.barrier() # mede IX
qc.h(1)
qc.cx(1,2) 
qc.barrier() # mede XX
qc.cx(0,2)
qc.cx(1,2)
qc.barrier()
qc.measure(2,0)
qc.draw('mpl')

job_sim = execute(qc, backend=simulator, shots=nshots)
counts_sim = job_sim.result().get_counts(qc)
job_exp = execute(qc, backend=device, shots=nshots)
print(job_exp.job_id())
job_monitor(job_exp)
counts_exp = job_exp.result().get_counts(qc)

counts_exp

plot_histogram([counts_sim, counts_exp], legend=['sim', 'exp'])

qr = QuantumRegister(3)
cr = ClassicalRegister(1)
qc = QuantumCircuit(qr,cr)
qc.barrier() # mede YI
qc.sdg(0)
qc.h(0)
qc.cx(0,2) 
qc.barrier() # mede IY
qc.sdg(1)
qc.h(1)
qc.cx(1,2) 
qc.barrier() # mede YY
qc.cx(0,2)
qc.cx(1,2)
qc.barrier()
qc.measure(2,0)
qc.draw('mpl')

job_sim = execute(qc, backend=simulator, shots=nshots)
counts_sim = job_sim.result().get_counts(qc)
job_exp = execute(qc, backend=device, shots=nshots)
print(job_exp.job_id())
job_monitor(job_exp)
counts_exp = job_exp.result().get_counts(qc)
counts_exp

plot_histogram([counts_sim, counts_exp], legend=['sim', 'exp'])

qr = QuantumRegister(3)
cr = ClassicalRegister(1)
qc = QuantumCircuit(qr,cr)
qc.barrier()
qc.h(0)
qc.cx(0,2) # mede XI
qc.barrier()
qc.sdg(1)
qc.h(1)
qc.cx(1,2) # mede IX
qc.barrier()
qc.cx(0,2); qc.cx(1,2)
qc.barrier()
qc.measure(2,0)
qc.draw('mpl')

job_sim = execute(qc, backend=simulator, shots=nshots)
counts_sim = job_sim.result().get_counts(qc)
job_exp = execute(qc, backend=device, shots=nshots)
print(job_exp.job_id())
job_monitor(job_exp)
counts_exp = job_exp.result().get_counts(qc)
counts_exp

plot_histogram([counts_sim, counts_exp], legend=['sim', 'exp'])

qr = QuantumRegister(3)
cr = ClassicalRegister(1)
qc = QuantumCircuit(qr,cr)
qc.barrier()
qc.sdg(0)
qc.h(0)
qc.cx(0,2) # mede XI
qc.barrier()
qc.h(1)
qc.cx(1,2) # mede IX
qc.barrier()
qc.cx(0,2); qc.cx(1,2)
qc.barrier()
qc.measure(2,0)
qc.draw('mpl')

job_sim = execute(qc, backend=simulator, shots=nshots)
counts_sim = job_sim.result().get_counts(qc)
job_exp = execute(qc, backend=device, shots=nshots)
print(job_exp.job_id())
job_monitor(job_exp)
counts_exp = job_exp.result().get_counts(qc)
counts_exp

plot_histogram([counts_sim, counts_exp], legend=['sim', 'exp'])

I = Matrix([[1,0],[0,1]]); X = Matrix([[0,1],[1,0]]); 
Y = Matrix([[0,-1j],[1j,0]]); Z = Matrix([[1,0],[0,-1]])

I,X,Y,Z

XX = kronecker_product(X,X); XX.eigenvects()

YY = kronecker_product(Y,Y); YY.eigenvects()

ZZ = kronecker_product(Z,Z); ZZ.eigenvects()

qr = QuantumRegister(3)
cr = ClassicalRegister(1)
qc = QuantumCircuit(qr,cr)
qc.barrier() # mede ZZ
qc.cx(0,2)
qc.cx(1,2)
qc.barrier() # mede XX
qc.h([0,1])
qc.cx(0,2)
qc.cx(1,2)
qc.barrier() # mede YY
qc.h([0,1])
qc.sdg([0,1])
qc.h([0,1])
qc.cx(0,2)
qc.cx(1,2)
qc.barrier()
qc.measure(2,0)
qc.draw('mpl')

job_sim = execute(qc, backend=simulator, shots=nshots)
counts_sim = job_sim.result().get_counts(qc)
job_exp = execute(qc, backend=device, shots=nshots)
print(job_exp.job_id())
job_monitor(job_exp)
counts_exp = job_exp.result().get_counts(qc)
counts_exp

plot_histogram([counts_sim, counts_exp], legend=['sim', 'exp'])

XY = kronecker_product(X,Y); XY.eigenvects()

YX = kronecker_product(Y,X); YX.eigenvects()

qr = QuantumRegister(3)
cr = ClassicalRegister(1)
qc = QuantumCircuit(qr,cr)
qc.barrier() # mede ZZ
qc.cx(0,2)
qc.cx(1,2)
qc.barrier() # mede XX
qc.sdg(1)
qc.h([0,1])
qc.cx(0,2)
qc.cx(1,2)
qc.barrier() # mede YY
qc.h([0,1])
qc.sdg(0)
qc.s(1)
qc.h([0,1])
qc.cx(0,2)
qc.cx(1,2)
qc.barrier()
qc.measure(2,0)
qc.draw('mpl')

job_sim = execute(qc, backend=simulator, shots=nshots)
counts_sim = job_sim.result().get_counts(qc)
job_exp = execute(qc, backend=device, shots=nshots)
print(job_exp.job_id())
job_monitor(job_exp)
counts_exp = job_exp.result().get_counts(qc)
counts_exp

plot_histogram([counts_sim, counts_exp], legend=['sim', 'exp'])

cr1 = {'0': 7921, '1': 271}
cr2 = {'0': 7944, '1': 248}
cr3 = {'0': 7754, '1': 438}
cc1 = {'0': 7913, '1': 279}
cc2 = {'0': 7940, '1': 252}
cc3 = {'0': 610, '1': 7582}

r1a = (cr1['0']-cr1['1'])/(cr1['0']+cr1['1'])
r2a = (cr2['0']-cr2['1'])/(cr2['0']+cr2['1'])
r3a = (cr3['0']-cr3['1'])/(cr3['0']+cr3['1'])
c1a = (cc1['0']-cc1['1'])/(cc1['0']+cc1['1'])
c2a = (cc2['0']-cc2['1'])/(cc2['0']+cc2['1'])
c3a = (cc3['0']-cc3['1'])/(cc3['0']+cc3['1'])
print(r1a,r2a,r3a,c1a,c2a,c3a)

def objective(x):
    mu_Q = np.array([0.93,0.94,0.89,0.93,0.94,-0.85])
    vec = mu_Q - A_R@x
    return np.linalg.norm(vec)

constraints = [{'type': 'eq', 'fun': lambda x: np.sum(x)-1},
               {'type':'ineq', 'fun': lambda x: x},
               {'type':'ineq', 'fun': lambda x: 1-x}]

np.random.seed(130000)
x0 = rpv_zhsl(16)

result = minimize(objective, x0, constraints=constraints, method='trust-constr')
print('melhor distribuicao de probabilidades', result.x)

print(objective(result.x))

