from sympy import *
init_printing()

from sympy.physics.quantum import *
from sympy.physics.quantum.spin import *

from sympy.physics.quantum.pauli import *
from sympy.physics.quantum.boson import *
from sympy.physics.quantum.density import *
from sympy.physics.quantum.operatorordering import *

s_ops = sx, sy, sz, sp, sm = SigmaX(), SigmaY(), SigmaZ(), SigmaPlus(), SigmaMinus()

#s_ops = sx, sy, sz, sp, sm = Jx, Jy, Jz, Jplus, Jminus  # does not work well, many expressions unevaluated

sx1, sy1, sz1 = SigmaX(1), SigmaY(1), SigmaZ(1)
sx2, sy2, sz2 = SigmaX(2), SigmaY(2), SigmaZ(2)

qsimplify_pauli(sx * sy)

qsimplify_pauli(sx1 * sy1 * sx2 * sy2)



from IPython.display import HTML

def print_table(data):
    t_table = "<table>\n%s\n</table>"
    t_row = "<tr>%s</tr>"
    t_col = "<td>%s</td>"
    table_code = t_table % "".join([t_row % "".join([t_col % ("$%s$" % latex(col))
                                                     for col in row])
                                    for row in data])
    return HTML(table_code)

data = [s_ops, [represent(s) for s in s_ops]]

print_table(data)

data = [[r"{\rm Expression}", r"{\rm Evaluated\;value}", r"{\rm Numerical\;test}"]] + \
       [[Commutator(s1, s2),
         Commutator(s1, s2).doit(),
         represent(Commutator(s1, s2)) == represent(Commutator(s1, s2).doit())]
        for s1 in s_ops for s2 in s_ops if s1 != s2]

print_table(data)

def elim_zero_matrix(e):
    return 0 if e == Matrix([[0, 0], [0, 0]]) else e

I_op = IdentityOperator(2)

data = [[r"{\rm Expression}", r"{\rm Evaluated\;value}", r"{\rm Numerical\;test}"]] + \
       [[AntiCommutator(s1, s2),
         AntiCommutator(s1, s2).doit(),
         elim_zero_matrix(represent(AntiCommutator(s1, s2))) == represent(I_op * AntiCommutator(s1, s2).doit())]
        for s1 in s_ops for s2 in s_ops if s1 != s2]

print_table(data)

Dagger(sx), Dagger(sy), Dagger(sz)

Dagger(sm), Dagger(sp)

sm, sp

s = [sx, sy, sz, sm, sp]

data = [[s2] + [(qsimplify_pauli(s1 * s2), represent((qsimplify_pauli(s1 * s2) * I_op).expand()) == represent(s1) * represent(s2))
         for s1 in s] for s2 in s]

print_table(data)

represent(sm) * represent(sm)

[sx ** n for n in range(10)]

[sy ** n for n in range(10)]

[sz ** n for n in range(10)]

[sm ** n for n in range(10)]

[sp ** n for n in range(10)]

up = SigmaZKet(0)
down = SigmaZKet(1)

ket = (up + down)/sqrt(2)

ket, qapply(sx * ket), qapply(sy * ket), qapply(sz * ket)

ket, qapply(sm * ket), qapply(sp * ket)

represent(qapply(sx * ket)) == represent(sx) * represent(ket)

represent(qapply(sy * ket)) == represent(sy) * represent(ket)

represent(qapply(sz * ket)) == represent(sz) * represent(ket)

represent(qapply(sm * ket)) == represent(sm) * represent(ket)

represent(qapply(sp * ket)) == represent(sp) * represent(ket)

from sympy_quantum_extra import *

rc = recursive_commutator(sy, sz, 5).doit().expand()

rc

A, B = Operator("A"), Operator("B")

Commutator(A, B).doit().doit()
# Add.doit() - > Mul.doit() -> ...

Commutator(sx, sy).doit()
# sx * sy - sy * sx

A = Matrix([[1, 2], [3, 4]])
B = Matrix([[1, 2], [3, 4]])
A * B, MatMul(A, B)

1 * 3

IdentityOperator(2).is_commutative

qsimplify(qsimplify(rc))

represent(sx * sx)

represent(sx) ** 2 == Integer(1) * represent(IdentityOperator(2))



theta = symbols("theta")

U = exp((I * theta / 2)* sy)

unitary_transformation(U, sx)

unitary_transformation(U, sz)

unitary_transformation(U, sm)

unitary_transformation(U, sp)

%reload_ext version_information

%version_information sympy