#!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

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

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

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

H2.inv()

from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import SamplerV2 as Sampler
import numpy as np

def tomo_rho_1qb_sim(nqb, qbt, qct, nshots, backend): 
    # qct = quantum circuit for tomography
    # nqb = number of qubits
    # qbt = qubit to be tomographed    
    sampler = Sampler(backend=backend)
    avgs = []
    for j in range(0,3):
        qc_list = []
        qc = QuantumCircuit(nqb,1)
        qc.append(qct,list(range(nqb)))
        if j == 0: # medida de X
            qc.h(qbt)
        elif j == 1: # medida de Y
            qc.sdg(qbt)
            qc.h(qbt)
        qc.measure(qbt,0)
        qc_list.append(qc.decompose())
        job = sampler.run(qc_list, shots=nshots)
        counts = job.result()[0].data.c.get_counts()
        if '0' in counts:
            avg = counts['0']
        if '1' in counts:
            avg -= counts['1']
        avg = avg/nshots
        avgs.append(avg)        
    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]]) 
    rho_1qb = 0.5*(s0 + avgs[0]*s1 + avgs[1]*s2 + avgs[2]*s3)
    return rho_1qb

nqb = 2; qbt = 0
qct = QuantumCircuit(nqb)
qct.h(0)
qct.cx(0,1)
nshots = 2**13
backend = AerSimulator()
rho = tomo_rho_1qb(nqb, qbt, qct, nshots, backend)
print(rho)

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

H4i = H4.inv(); H4i

def hadarmard_inv():
    return (1/4)*np.array([[1,1,1,1],[1,-1,1,-1],[1,1,-1,-1],[1,-1,-1,1]])
                           
#hadarmard_inv()                        

def tomo_rho_2qb(nqb, qbt, qct, nshots, backend): 
    # qct = quantum circuit for tomography
    # nqb = number of qubits
    # qbt = list with the 2 qubits to be tomographed    
    sampler = Sampler(backend=backend)
    avgs = []
    for j in range(1,4):
        for k in range(1,4):
            qc_list = []
            qc = QuantumCircuit(nqb,2)
            qc.append(qct,list(range(nqb)))
            if j == 1: # medida de X_0
                qc.h(qbt[0])
            elif j == 2: # medida de Y_0
                qc.sdg(qbt[0])
                qc.h(qbt[0])
            if k == 1: # medida de X_1
                qc.h(qbt[1])
            elif k == 2: # medida de Y_1
                qc.sdg(qbt[1])
                qc.h(qbt[1])
            qc.measure(qbt,[0,1]) # medida de ZZ
            qc_list.append(qc.decompose())
            job = sampler.run(qc_list, shots=nshots)
            counts = job.result()[0].data.c.get_counts()
            labels = ['00','01','10','11']
            probs = np.zeros(4)
            for l in range(0,4):
                if labels[l] in counts:
                    probs[l] = counts[labels[l]]/nshots
            print(probs)            
    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]]) 
    rho_2qb = 0.5*(s0 + avgs[0]*s1 + avgs[1]*s2 + avgs[2]*s3)
    return rho_2qb

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 

nqb = 2; qbt = [0,1]
qct = QuantumCircuit(nqb)
qct.h(0)
qct.cx(0,1)
nshots = 2**13
backend = AerSimulator()
rho = tomo_rho_2qb(nqb, qbt, qct, nshots, backend)
print(rho)

print(np.zeros(4))



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

 



















