A brief tutorial how to simulate quantum teleportation with QuTip quantum computing library in Python. More detailed explanation available here: http://mareknarozniak.com/2020/03/22/simulating-quantum-teleportation/

teleport.py

import numpy as np

import itertools

from qutip import basis, tensor, rand_ket, snot, cnot, rx, rz, qeye

# create random state
psi = rand_ket(2)

# initial state
psi0 = tensor([psi, basis(2, 0), basis(2, 0)])

# unitary time-evolution for quantum teleportation
psi1 = snot(N=3, target=1)*psi0
psi2 = cnot(N=3, control=1, target=2)*psi1
psi3 = cnot(N=3, control=0, target=1)*psi2
psi4 = snot(N=3, target=0)*psi3

# all the possible outcomes of measurements on first two qubits
confs = list(itertools.product([0, 1], repeat=2))

# projection operators
Ps = []
for m0, m1 in confs:
    P = tensor([
        basis(2, m0).proj(),
        basis(2, m1).proj(),
        qeye(2)])
    Ps.append(P)

# simulate measurement by performing the projection
psis_proj = []
for P in Ps:
    psi_proj = (P*psi4).unit()
    psis_proj.append(psi_proj)

# classical correction operators
X = rx(np.pi, N=3, target=2)
Z = rz(np.pi, N=3, target=2)

# perform classical correction
psis_corr = [
    psis_proj[0],
    X*psis_proj[1],
    Z*psis_proj[2],
    Z*X*psis_proj[3]
]

# produce reference states for fidelity calculation
psis_ref = []
for m0, m1 in confs:
    psi_ref = tensor([basis(2, m0), basis(2, m1), psi])
    psis_ref.append(psi_ref)

print()
print('classically corrected outcomes fidelities')
print('{0:2} {1:2} {2:8}'.format('m0', 'm1', 'fidelity'))
for conf, psi_corr, psi_ref in zip(confs, psis_corr, psis_ref):
    fidelity = np.round(np.abs(psi_corr.overlap(psi_ref))**2., 3)
    m0, m1 = conf
    print('{0:2} {1:2} {2:8}'.format(m0, m1, fidelity))

print()
print('fidelities without classical correction')
print('{0:2} {1:2} {2:8}'.format('m0', 'm1', 'fidelity'))
for conf, psi_proj, psi_ref in zip(confs, psis_proj, psis_ref):
    fidelity = np.round(np.abs(psi_proj.overlap(psi_ref))**2., 3)
    m0, m1 = conf
    print('{0:2} {1:2} {2:8}'.format(m0, m1, fidelity))