#!pip install pylatexenc # para desenhar os circuitos quânticos

from qiskit import QuantumCircuit
qc = QuantumCircuit(1,1)
qc.barrier()
qc.measure(0,0)
qc.draw('mpl') # Para medidas de Z

qc = QuantumCircuit(1,1)
qc.barrier()
qc.h(0)
qc.measure(0,0)
qc.draw('mpl') # Para medidas de X

qc = QuantumCircuit(1,1)
qc.barrier()
qc.sdg(0)
qc.h(0)
qc.measure(0,0)
qc.draw('mpl') # Para medidas de Y

# Calculo da inversa da matrix de Hadamard
from sympy import Matrix
H2 = Matrix([[1,1],[1,-1]]); H2

H2.inv()

def tomo_1qb(avgs):
    import numpy as np
    s0 = np.array([[1,0],[0,1]]); s1 = np.array([[0,1],[1,0]])
    s2 = np.array([[0,-1j],[1j,0]]); s3 = np.array([[1,0],[0,-1]])
    return 0.5*(s0 + avgs[1]*s1 + avgs[2]*s2 + avgs[3]*s3)

import numpy as np
avgs = np.array([1,0.1,-0.2,0.5]); rho1qb = tomo_1qb(avgs); rho1qb

def qc_test_1qb():
    from qiskit import QuantumCircuit
    qc = QuantumCircuit(1, name='qc')
    qc.h(0)
    return qc

from qiskit import QuantumCircuit, execute
from qiskit import Aer
simulator = Aer.get_backend('qasm_simulator')
nshots = 2**13

qct = qc_test_1qb() # circuito quantico principal

# medida de X
qc = QuantumCircuit(1,1)
qc.append(qct,[0])
qc.h(0); qc.measure(0,0)
qcx = qc
job = execute(qc, backend = simulator, shots = nshots)
counts_x = job.result().get_counts()
print('counts_x=',counts_x)
if '0' in counts_x:
    Xavg = counts_x['0']
if '1' in counts_x:
    Xavg -= counts_x['1']
Xavg = Xavg/nshots
print(Xavg)

# medida de Y
qc = QuantumCircuit(1,1)
qc.append(qct,[0])
qc.sdg(0); qc.h(0); qc.measure(0,0)
qcy = qc
job = execute(qc, backend = simulator, shots = nshots)
counts_y = job.result().get_counts()
print('counts_y=',counts_y)
if '0' in counts_y:
    Yavg = counts_y['0']
if '1' in counts_y:
    Yavg -= counts_y['1']
Yavg = Yavg/nshots
print(Yavg)

# medida de Z
qc = QuantumCircuit(1,1)
qc.append(qct,[0])
qc.measure(0,0)
qcz = qc
job = execute(qc, backend = simulator, shots = nshots)
counts_z = job.result().get_counts()
print('counts_z=',counts_z)
if '0' in counts_z:
    Zavg = counts_z['0']
if '1' in counts_y:
    Zavg -= counts_z['1']
Zavg = Zavg/nshots
print(Zavg)

#qcx.draw('mpl') #ok
qcy.draw('mpl')
#qcz.draw('mpl')

avgs = [1, Xavg, Yavg, Zavg]
rho1qb = tomo_1qb(avgs)
rho1qb

from numpy import linalg as LA
evals, evecs = LA.eig(rho1qb)
evals

evecs

from sympy.physics.quantum import TensorProduct as tp
H4 = tp(H2,H2); H4

H4i = H4.inv(); H4i

def tomo_2qb(avgs):
    # avgs e uma lista com 16 valores medios na seguinte ordem [II,IX,IY,IZ,XI,XX,XY,XZ,YI,YX,YY,YZ,ZI,ZX,ZY,ZZ]
    import numpy as np
    s0 = np.array([[1,0],[0,1]]); s1 = np.array([[0,1],[1,0]])
    s2 = np.array([[0,-1j],[1j,0]]); s3 = np.array([[1,0],[0,-1]])
    paulis = [s0,s1,s2,s3] # lista com as matrizes de Pauli
    rho = np.zeros([4,4],dtype=complex)
    l = 0
    for j in range(0,3):
        for k in range(0,3):
            l += 1
            rho += avgs[l]*tp(paulis[j],paulis[k])
    return 0.25*rho       

avgs = [1,0.1,0,0,0.1,0.3,0,0,0.3,0,0.2,0,0,0.1,0,0.5] #[II,IX,IY,IZ,XI,XX,XY,XZ,YI,YX,YY,YZ,ZI,ZX,ZY,ZZ]
rho = tomo_2qb(avgs); rho

# Inversa da matriz de Hadamard 4x4 (sem o multiplo 1/4)
H4i = np.array([[1,1,1,1],[1,-1,1,-1],[1,1,-1,-1],[1,-1,-1,1]]); H4i

def qc_test_2qb():
    from qiskit import QuantumCircuit
    qc = QuantumCircuit(2, name='qc2qb')
    qc.h(0); qc.cx(0,1)
    return qc

def tomo_2qb(avgs):
    # avgs e uma lista com 64 valores medios na seguinte ordem [III,IIX,IIY,IIZ,IXI,IXX,IXY,IXZ,...]
    import numpy as np
    s0 = np.array([[1,0],[0,1]]); s1 = np.array([[0,1],[1,0]])
    s2 = np.array([[0,-1j],[1j,0]]); s3 = np.array([[1,0],[0,-1]])
    paulis = [s0,s1,s2,s3]
    rho = np.zeros([4,4],dtype=complex)
    m = 0
    for j in range(0,4):
        for k in range(0,4):
            for l in range(0,4):
                rho += avgs[m]*np.kron(np.kron(paulis[j],paulis[k]),paulis[l])
                m += 1
    return 0.25*rho

from qiskit import QuantumCircuit, execute
from qiskit import Aer
simulator = Aer.get_backend('qasm_simulator')
nshots = 2**13

qct = qc_test_2qb() # circuito quantico principal


# medida de XX
qc = QuantumCircuit(2,2); qc.append(qct,[0,1]); qc.h([0,1]); qc.measure([0,1],[0,1]); qcxx = qc
job = execute(qc, backend = simulator, shots = nshots); counts_xx = job.result().get_counts()
pr_xx = np.zeros(4); cb2_labels = ['00','10','01','11']
for j in range(0,4):
    if cb2_labels[j] in counts_xx:
        pr_xx[j] = counts_xx[cb2_labels[j]]/nshots
print('counts_xx=',counts_xx, pr_xx)   
XX_vec = np.matmul(H4i,pr_xx); print('XX_vec=',XX_vec) # =[II,IX,XI,XX]

qcxx.draw('mpl')

ss_vecs = [] # list with the averages of the observables
obs = ['x','y','z']
for j in range(0,3):
    lo = obs[j] # left observable
    for k in range(0,3):
        ro = obs[k] # right observable
        qc = QuantumCircuit(2,2)
        qc.append(qct,[0,1]) # the quantum circuit for preparing psi
        if lo == 'x' and ro = 'x':
            qc.h([0,1])
        elif lo == 'x' and ro = 'y':
            qc.sdg(1); qc.h([0,1])
        elif lo == 'y' and ro = 'x':
            qc.sdg(0); qc.h([0,1])
        elif lo == 'y' and ro = 'y':
            qc.sdg([0,1]); qc.h([0,1])        
        qc.measure([0,1],[0,1])
        job = execute(qc, backend = simulator, shots = nshots)
        counts = job.result().get_counts()
        pr = np.zeros(4); cb2_labels = ['00','10','01','11']
        for j in range(0,4):
            if cb2_labels[j] in counts_xx:
                pr_xx[j] = counts_xx[cb2_labels[j]]/nshots
        ss_vec = np.matmul(H4i,pr)
        ss_vecs.append(ss_vec)









# matriz de Hadamard
H8 = tp(H2,H4); H8

H8i = H8.inv(); H8i

from sympy import symbols
f000,f001,f010,f011,f100,f101,f110,f111 = symbols('f_{000},f_{001},f_{010},f_{011},f_{100},f_{101},f_{110},f_{111}')
f = Matrix([[f000],[f001],[f010],[f011],[f100],[f101],[f110],[f111]])
f.T

al = 8*H8i*f; al

import numpy as np
s0 = np.array([[1,0],[0,1]]); s1 = np.array([[0,1],[1,0]]); s2 = np.array([[0,-1j],[1j,0]]); s3 = np.array([[1,0],[0,-1]])
paulis = [s0,s1,s2,s3] # lista com as matrizes de Pauli
list_avgs = []
m = 0
for j in range(0,4):
        for k in range(0,4):
            for l in range(0,4):
                list_avgs.append((j,k,l))
                m += 1

 



















