from qiskit import *
%matplotlib inline
import numpy as np

# Podemos encontrar o operador de densidade de estados

# Para isso vamos primeiramente usar o statevector_simulator para encontrar o vetor de 
# estado desejado

qc = QuantumCircuit(1)
qc.ry(3*np.pi/2,0)
qc.draw(output='mpl')

backend = Aer.get_backend('statevector_simulator')
job = execute(qc,backend)

result = job.result()
statevector = result.get_statevector(qc, decimals = 3)

# Aqui está o nosso vetor de estado
print(statevector)

# Agora para encontrar o seu operador de densidade
from qiskit.quantum_info import *
op_densidade = DensityMatrix(statevector)

print(op_densidade.data)

# É possível também visualizar a matriz de densidade por meio de um histograma

from qiskit.visualization import plot_state_city

plot_state_city(statevector)

# Também é possivel visualizar as matrizes de Pauli que compõem o seu operador de densidade

from qiskit.tools.visualization import plot_state_paulivec

plot_state_paulivec(statevector)

# Podemos ver o traço do nosso operador
op_densidade.trace()

estado = [np.sqrt(2),np.sqrt(2)]
print(np.inner(np.inner(estado, op_densidade.data), estado))

rho_2 = np.matmul(op_densidade.data,op_densidade.data)
trace = rho_2.trace()
print(trace)

def bloch_vector(op):
    pauli_x = [[0,1],[1,0]]
    pauli_y = [[0,-1j],[1j,0]]
    pauli_z = [[1,0],[0,-1]]
    
    x = np.matmul(op,pauli_x).trace()
    y = np.matmul(op,pauli_y).trace()
    z = np.matmul(op,pauli_z).trace()
    
    return [np.absolute(x),np.absolute(y),np.absolute(z)]

from qiskit.tools.visualization import plot_bloch_vector

plot_bloch_vector(bloch_vector(op_densidade.data))

# Exemplo de estado misto:


estado_misto = [[0.5, 0.25],[0.25,0.5]]
plot_bloch_vector(bloch_vector(estado_misto))