import numpy as np

from qiskit import QuantumCircuit, BasicAer, execute, Aer
from qiskit.quantum_info import Operator, DensityMatrix, Pauli
from qiskit.extensions import YGate
from qiskit.visualization import plot_histogram

Cnot = Operator([[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]])
Cnot

#returns Numpy array
Cnot.data

#total input and output dimension
in_dim, out_dim = Cnot.dim
in_dim, out_dim

ghz = QuantumCircuit(3)
ghz.h(0)
ghz.cx([0,0],[1,2])
ghz.draw('mpl')

U_ghz = Operator(ghz) 
np.around(U_ghz.data,3)

#np.around specifies how many decimals is desired in the oputput 

#Aer has also a unitary_simulator backend.
back_uni = Aer.get_backend('unitary_simulator')
job = execute(ghz, back_uni)
result = job.result()

#getting the unitary with 3 decimals
U_qc = result.get_unitary(decimals=3)

U_qc 

op1 = [[1+0.j, 0.5+0.j], [0.5+0.j, 1+0.j]]
op2 = [[0.5+0.j, 1+0.j], [0.5+0.j, 1+0.j]]

matrix = DensityMatrix(op1)

#doing the tensor product between op1 and op2
matrix.tensor(op2)

# Pauli('Y') generates the Pauli Y-matrix object 
# On the other hand, Operator generates the operator of Pauli('Y') 

Operator( Pauli('Y')) == Operator(YGate())

Operator(YGate()) == np.exp(1j * 0.4) * Operator(YGate())

ZZ = Operator(Pauli('ZZ'))
ZZ

# Add to a circuit
qc = QuantumCircuit(2, 2)
qc.append(ZZ, [0, 1])
qc.measure([0,1], [0,1])

qc.draw('mpl')

backend = BasicAer.get_backend('qasm_simulator')
result = backend.run(qc).result()

counts = result.get_counts()
plot_histogram(counts)

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

op1 = Operator(Pauli('Y'))
op2 = Operator(Pauli('Z'))
op1.compose(op2, front=True)

op1 = Operator(Pauli('X'))
op2 = Operator(Pauli('Y'))
op1.tensor(op2)

