P1 = matrix(QQ,[[1,0,0,1],[0,0,0,0],[0,0,0,0],[1,0,0,1]])/2
P2 = matrix(QQ,[[1,0,1,0],[0,0,0,0],[1,0,1,0],[0,0,0,0]])/2
var('rho00 rho01 rho10 rho11')
rho_gen = matrix(SR,[[rho00,rho01],[rho10,rho11]])
ket0bra0 = matrix(QQ,[[1,0],[0,0]])
rho = ket0bra0

print html('Definitions: $P_1={}$, $P_2={}$, $\rho_\mathit{{generic}}={}$, $\lvert0\rangle\langle0\rvert={}$, $\rho={}$'.format(latex(P1),latex(P2),latex(rho_gen),latex(ket0bra0),latex(rho)))

P1_is_proj = (P1==P1.conjugate_transpose()) and (P1==P1*P1)
print html("$P_1$ is a projector: {}".format(P1_is_proj))
P2_is_proj = (P2==P2.conjugate_transpose()) and (P2==P2*P2)
print html("$P_2$ is a projector: {}".format(P2_is_proj))

pr_P1 = (P1*ket0bra0.tensor_product(rho_gen)).trace()
pr_P2 = (P2*ket0bra0.tensor_product(rho_gen)).trace()
print html('$\operatorname{{tr}}\Bigl(P_1\bigl(\lvert0\rangle\langle0\rvert\otimes\rho_\mathit{{generic}}\bigr)Bigr) = {}$'.format(latex(pr_P1)))
print html('$\operatorname{{tr}}\Bigl(P_2\bigl(\lvert0\rangle\langle0\rvert\otimes\rho_\mathit{{generic}}\bigr)Bigr) = {}$'.format(latex(pr_P2)))
print html("Measurement success is independent of whether we use $P_1$ or $P_2$: {}".format(bool(pr_P1==pr_P2)))

def traceout_A(M):
    s = 0
    for i in (0,1):
        ket_i = matrix(QQ,2,1); ket_i[i,0] = 1
        id2 = matrix.identity(2)
        s += ket_i.conjugate_transpose().tensor_product(id2) * M * ket_i.tensor_product(id2)
    return s

pms1 = traceout_A(P1*ket0bra0.tensor_product(rho)*P1)
pms2 = traceout_A(P2*ket0bra0.tensor_product(rho)*P2)

print html("Post-measurement state $\operatorname{{tr}}_A \Bigl(P_1\bigl(\lvert0\rangle\langle0\rvert\otimes\rho\bigr)P_1Bigr)$ is ${}$".format(latex(pms1)))
print html("Post-measurement state $\operatorname{{tr}}_A \Bigl(P_2\bigl(\lvert0\rangle\langle0\rvert\otimes\rho\bigr)P_2Bigr)$ is ${}$".format(latex(pms2)))

print html("The post measurement states are different: {}".format(pms1 <> pms2))