import numpy as np
zero_ket = np.array([[1], [0]])
print("|0> ket:\n", zero_ket)
print("<0| bra:\n", zero_ket.T.conj())

zero_ket.T.conj() @ zero_ket

one_ket = np.array([[0], [1]])
zero_ket.T.conj() @ one_ket

zero_ket @ zero_ket.T.conj()

ψ = np.array([[1], [0]])/np.sqrt(2)
Π_0 = zero_ket @ zero_ket.T.conj()
ψ.T.conj() @ Π_0 @ ψ

from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister

from qiskit_aer import AerSimulator
from qiskit.visualization import plot_histogram
from qiskit.primitives import BackendEstimator, BackendSampler

backend = AerSimulator()
sampler = BackendSampler(backend)
q = QuantumRegister(1)
c = ClassicalRegister(1)
circuit = QuantumCircuit(q, c)
circuit.h(q[0])
circuit.measure(q, c)

job = sampler.run(circuit, shots=100)
result = job.result()
counts = result.quasi_dists[0]
plot_histogram(counts)


ψ = np.array([[np.sqrt(2)/2], [np.sqrt(2)/2]])
Π_0 = zero_ket @ zero_ket.T.conj()
probability_0 = ψ.T.conj() @ Π_0 @ ψ
Π_0 @ ψ/np.sqrt(probability_0)


c = ClassicalRegister(2)
circuit = QuantumCircuit(q, c)
circuit.h(q[0])
circuit.measure(q[0], c[0])
circuit.measure(q[0], c[1])

job = sampler.run(circuit, shots=100)
result = job.result()
counts = result.quasi_dists[0]
plot_histogram(counts)



q = QuantumRegister(2)
c = ClassicalRegister(2)
circuit = QuantumCircuit(q, c)
circuit.h(q[0])
circuit.measure(q, c)

job = sampler.run(circuit, shots=100)
result = job.result()
counts = result.quasi_dists[0]
plot_histogram(counts)


q = QuantumRegister(2)
c = ClassicalRegister(2)
circuit = QuantumCircuit(q, c)
circuit.h(q[0])
circuit.cx(q[0], q[1])
circuit.measure(q, c)

job = sampler.run(circuit, shots=100)
result = job.result()
counts = result.quasi_dists[0]
plot_histogram(counts)



ψ = np.array([[1], [1]])/np.sqrt(2)
ρ = ψ @ ψ.T.conj()
Π_0 = zero_ket @ zero_ket.T.conj()
np.trace(Π_0 @ ρ)

probability_0 = np.trace(Π_0 @ ρ)
Π_0 @ ρ @ Π_0/probability_0

zero_ket = np.array([[1], [0]])
one_ket = np.array([[0], [1]])
ψ = (zero_ket + one_ket)/np.sqrt(2)
print("Density matrix of the equal superposition")
print(ψ @ ψ.T.conj())
print("Density matrix of the equally mixed state of |0><0| and |1><1|")
print((zero_ket @ zero_ket.T.conj()+one_ket @ one_ket.T.conj())/2)