adding Commutators, Anticommutator, Daggers (#40)

.gitignore

@@ -1 +1,2 @@
 Manifest.toml
+.vscode

src/QSymbolicsBase/QSymbolicsBase.jl

@@ -6,7 +6,7 @@ using TermInterface
 import TermInterface: isexpr, head, iscall, children, operation, arguments, metadata
 
 using LinearAlgebra
-import LinearAlgebra: eigvecs
+import LinearAlgebra: eigvecs, ishermitian, inv
 
 import QuantumInterface:
     apply!,
@@ -17,24 +17,30 @@ import QuantumInterface:
     AbstractKet, AbstractOperator, AbstractSuperOperator, AbstractBra
 
 export SymQObj,QObj,
-       AbstractRepresentation, AbstractUse,
-       QuantumOpticsRepr, QuantumMCRepr, CliffordRepr,
-       UseAsState, UseAsObservable, UseAsOperation,
+       AbstractRepresentation,AbstractUse,
+       QuantumOpticsRepr,QuantumMCRepr,CliffordRepr,
+       UseAsState,UseAsObservable,UseAsOperation,
        apply!,
        express,
        tensor,⊗,
-       dagger,projector,
-       X,Y,Z,σˣ,σʸ,σᶻ,Pm,Pp,σ₋,σ₊,
+       dagger,projector,commutator,anticommutator,expand,
+       I,X,Y,Z,σˣ,σʸ,σᶻ,Pm,Pp,σ₋,σ₊,
        H,CNOT,CPHASE,XCX,XCY,XCZ,YCX,YCY,YCZ,ZCX,ZCY,ZCZ,
        X1,X2,Y1,Y2,Z1,Z2,X₁,X₂,Y₁,Y₂,Z₁,Z₂,L0,L1,Lp,Lm,Lpi,Lmi,L₀,L₁,L₊,L₋,L₊ᵢ,L₋ᵢ,
        vac,F₀,F0,F₁,F1,
-       N,n̂,Create,âꜛ,Destroy,â,
-       SProjector,MixedState,IdentityOp,
-       STensorKet,STensorOperator,SScaledKet,SScaledOperator,SAddKet,SAddOperator,SScaledBra,SAddBra,STensorBra,SDagger,
-       HGate, XGate, YGate, ZGate, CPHASEGate, CNOTGate,
-       XBasisState, YBasisState, ZBasisState,
-       NumberOp, CreateOp, DestroyOp,
-       XCXGate, XCYGate, XCZGate, YCXGate, YCYGate, YCZGate, ZCXGate, ZCYGate, ZCZGate
+       N,n̂,Create,âꜛ,Destroy,â,SpinBasis,FockBasis,
+       SBra,SKet,SOperator,
+       SAdd,SAddBra,SAddKet,SAddOperator,
+       SScaled,SScaledBra,SScaledOperator,SScaledKet,
+       STensorBra,STensorKet,STensorOperator,
+       SProjector,MixedState,IdentityOp,SInvOperator,SHermitianOperator,SUnitaryOperator,SHermitianUnitaryOperator,
+       SApplyKet,SApplyBra,SMulOperator,SSuperOpApply,SCommutator,SAnticommutator,SDagger,SBraKet,SOuterKetBra,
+       HGate,XGate,YGate,ZGate,CPHASEGate,CNOTGate,
+       XBasisState,YBasisState,ZBasisState,
+       NumberOp,CreateOp,DestroyOp,
+       XCXGate,XCYGate,XCZGate,YCXGate,YCYGate,YCZGate,ZCXGate,ZCYGate,ZCZGate,
+       qsimplify,qsimplify_pauli,qsimplify_flatten,qsimplify_commutator,qsimplify_anticommutator,
+       isunitary
 
 function countmap(samples) # A simpler version of StatsBase.countmap, because StatsBase is slow to import
     counts = Dict{Any,Any}()
@@ -127,17 +133,29 @@ newwithmetadata(x) = x
 # Basic Types
 ##
 
-const QObj = Union{AbstractKet,AbstractOperator,AbstractSuperOperator,AbstractBra}
+const QObj = Union{AbstractBra,AbstractKet,AbstractOperator,AbstractSuperOperator}
 const SymQObj = Symbolic{<:QObj} # TODO Should we use Sym or Symbolic... Sym has a lot of predefined goodies, including metadata support
 Base.:(-)(x::SymQObj) = (-1)*x
 Base.:(-)(x::SymQObj,y::SymQObj) = x + (-y)
 
-function Base.isequal(x::X,y::Y) where {X<:SymQObj, Y<:SymQObj}
+function _in(x::SymQObj, y::SymQObj)
+    for i in arguments(y)
+        if isequal(x, i)
+            return true
+        end
+    end
+    false
+end
+function Base.isequal(x::X,y::Y) where {X<:Union{SymQObj, Symbolic{Complex}}, Y<:Union{SymQObj, Symbolic{Complex}}}
     if X==Y
         if isexpr(x)
             if operation(x)==operation(y)
                 ax,ay = arguments(x),arguments(y)
-                (length(ax) == length(ay)) && all(zip(ax,ay)) do xy isequal(xy...) end
+                if (operation(x) === +) && (length(ax) == length(ay))
+                    all(x -> _in(x, y), ax)
+                else
+                    all(zip(ax,ay)) do xy isequal(xy...) end
+                end
             else
                 false
             end
@@ -151,12 +169,13 @@ end
 
 # TODO check that this does not cause incredibly bad runtime performance
 # use a macro to provide specializations if that is indeed the case
-propsequal(x,y) = all(n->getproperty(x,n)==getproperty(y,n), propertynames(x))
+propsequal(x,y) = all(n->isequal(getproperty(x,n),getproperty(y,n)), propertynames(x))
 
 ##
 # Most symbolic objects defined here
 ##
 
+include("literal_objects.jl")
 include("basic_ops_homogeneous.jl")
 include("basic_ops_inhomogeneous.jl")
 include("predefined.jl")

src/QSymbolicsBase/basic_ops_homogeneous.jl

@@ -1,24 +1,23 @@
-"""This file defines the symbolic operations for quantum objects (kets, operators, and bras) that are homogeneous in their arguments."""
+##
+# This file defines the symbolic operations for quantum objects (kets, operators, and bras) that are homogeneous in their arguments.
+##
 
-struct SKet <: Symbolic{AbstractKet}
-    name::Symbol
-    basis::Basis
-end
-struct SOperator <: Symbolic{AbstractKet}
-    name::Symbol
-    basis::Basis
-end
-const SymQ = Union{SKet, SOperator}
-isexpr(::SymQ) = false
-metadata(::SymQ) = nothing
-basis(x::SymQ) = x.basis
+"""Scaling of a quantum object (ket, operator, or bra) by a number
+
+```jldoctest
+julia> k = SKet(:k, SpinBasis(1//2))
+|k⟩
 
-symbollabel(x::SymQ) = x.name
-Base.show(io::IO, x::SKet) = print(io, "|$(symbollabel(x))⟩")
-Base.show(io::IO, x::SOperator) = print(io, "$(symbollabel(x))")
-Base.show(io::IO, x::SymQObj) = print(io, symbollabel(x)) # fallback that probably is not great
+julia> 2*k
+2|k⟩
 
-"""Scaling of a quantum object (ket, operator, or bra) by a number."""
+julia> A = SOperator(:A, SpinBasis(1//2))
+A 
+
+julia> 2*A
+2A
+````
+"""
 @withmetadata struct SScaled{T<:QObj} <: Symbolic{T}
     coeff
     obj
@@ -54,13 +53,21 @@ end
 const SScaledBra = SScaled{AbstractBra}
 function Base.show(io::IO, x::SScaledBra)
     if x.coeff isa Number
-        print(io, "$(x.obj)$(x.coeff)")
+        print(io, "$(x.coeff)$(x.obj)")
     else
-        print(io, "$(x.obj)($(x.coeff))")
+        print(io, "($(x.coeff))$(x.obj)")
     end
 end
 
-"""Addition of quantum objects (kets, operators, or bras)."""
+"""Addition of quantum objects (kets, operators, or bras)
+
+```jldoctest
+julia> k₁ = SKet(:k₁, SpinBasis(1//2)); k₂ = SKet(:k₂, SpinBasis(1//2));
+
+julia> k₁ + k₂
+(|k₁⟩+|k₂⟩)
+```
+"""
 @withmetadata struct SAdd{T<:QObj} <: Symbolic{T}
     dict
     SAdd{S}(d) where S = length(d)==1 ? SScaled{S}(reverse(first(d))...) : new{S}(d)
@@ -76,13 +83,61 @@ Base.:(+)(xs::Vararg{Symbolic{<:QObj},0}) = 0 # to avoid undefined type paramete
 basis(x::SAdd) = basis(first(x.dict).first)
 
 const SAddKet = SAdd{AbstractKet}
-Base.show(io::IO, x::SAddKet) = print(io, "("*join(map(string, arguments(x)),"+")::String*")") # type assert to help inference
+function Base.show(io::IO, x::SAddKet)
+    ordered_terms = sort([repr(i) for i in arguments(x)])
+    print(io, "("*join(ordered_terms,"+")::String*")") # type assert to help inference
+end
 const SAddOperator = SAdd{AbstractOperator}
-Base.show(io::IO, x::SAddOperator) = print(io, "("*join(map(string, arguments(x)),"+")::String*")") # type assert to help inference
+function Base.show(io::IO, x::SAddOperator) 
+    ordered_terms = sort([repr(i) for i in arguments(x)])
+    print(io, "("*join(ordered_terms,"+")::String*")") # type assert to help inference
+end
 const SAddBra = SAdd{AbstractBra}
-Base.show(io::IO, x::SAddBra) = print(io, "("*join(map(string, arguments(x)),"+")::String*")") # type assert to help inference
+function Base.show(io::IO, x::SAddBra)
+    ordered_terms = sort([repr(i) for i in arguments(x)])
+    print(io, "("*join(ordered_terms,"+")::String*")") # type assert to help inference
+end
+
+"""Symbolic application of operator on operator
+
+```jldoctest
+julia> A = SOperator(:A, SpinBasis(1//2)); B = SOperator(:B, SpinBasis(1//2));
 
-"""Tensor product of quantum objects (kets, operators, or bras)."""
+julia> A*B 
+AB
+```
+"""
+@withmetadata struct SMulOperator <: Symbolic{AbstractOperator}
+    terms
+    function SMulOperator(terms)
+        coeff, cleanterms = prefactorscalings(terms)
+        coeff*new(cleanterms)
+    end
+end
+isexpr(::SMulOperator) = true
+iscall(::SMulOperator) = true
+arguments(x::SMulOperator) = x.terms
+operation(x::SMulOperator) = *
+head(x::SMulOperator) = :*
+children(x::SMulOperator) = [:*;x.terms]
+Base.:(*)(xs::Symbolic{AbstractOperator}...) = SMulOperator(collect(xs))
+Base.show(io::IO, x::SMulOperator) = print(io, join(map(string, arguments(x)),""))
+basis(x::SMulOperator) = basis(x.terms)
+
+"""Tensor product of quantum objects (kets, operators, or bras)
+
+```jldoctest
+julia> k₁ = SKet(:k₁, SpinBasis(1//2)); k₂ = SKet(:k₂, SpinBasis(1//2));
+
+julia> k₁ ⊗ k₂
+|k₁⟩|k₂⟩
+
+julia> A = SOperator(:A, SpinBasis(1//2)); B = SOperator(:B, SpinBasis(1//2));
+
+julia> A ⊗ B 
+A⊗B
+```
+"""
 @withmetadata struct STensor{T<:QObj} <: Symbolic{T}
     terms
     function STensor{S}(terms) where S
@@ -107,3 +162,62 @@ const STensorSuperOperator = STensor{AbstractSuperOperator}
 Base.show(io::IO, x::STensorSuperOperator) = print(io, join(map(string, arguments(x)),"⊗"))
 const STensorBra = STensor{AbstractBra}
 Base.show(io::IO, x::STensorBra) = print(io, join(map(string, arguments(x)),""))
+
+"""Symbolic commutator of two operators
+
+```jldoctest
+julia> A = SOperator(:A, SpinBasis(1//2)); B = SOperator(:B, SpinBasis(1//2));
+
+julia> commutator(A, B)
+[A,B]
+
+julia> commutator(A, A)
+0
+```
+"""
+@withmetadata struct SCommutator <: Symbolic{AbstractOperator}
+    op1
+    op2
+    function SCommutator(o1, o2) 
+        coeff, cleanterms = prefactorscalings([o1 o2], scalar=true)
+        cleanterms[1] === cleanterms[2] ? 0 : coeff*new(cleanterms...)
+    end
+end
+isexpr(::SCommutator) = true
+iscall(::SCommutator) = true
+arguments(x::SCommutator) = [x.op1, x.op2]
+operation(x::SCommutator) = commutator
+head(x::SCommutator) = :commutator
+children(x::SCommutator) = [:commutator, x.op1, x.op2]
+commutator(o1::Symbolic{AbstractOperator}, o2::Symbolic{AbstractOperator}) = SCommutator(o1, o2)
+Base.show(io::IO, x::SCommutator) = print(io, "[$(x.op1),$(x.op2)]")
+basis(x::SCommutator) = basis(x.op1)
+expand(x::SCommutator) = x == 0 ? x : x.op1*x.op2 - x.op2*x.op1
+
+"""Symbolic anticommutator of two operators
+
+```jldoctest
+julia> A = SOperator(:A, SpinBasis(1//2)); B = SOperator(:B, SpinBasis(1//2));
+
+julia> anticommutator(A, B)
+{A,B}
+```
+"""
+@withmetadata struct SAnticommutator <: Symbolic{AbstractOperator}
+    op1
+    op2
+    function SAnticommutator(o1, o2) 
+        coeff, cleanterms = prefactorscalings([o1 o2], scalar=true)
+        coeff*new(cleanterms...)
+    end
+end
+isexpr(::SAnticommutator) = true
+iscall(::SAnticommutator) = true
+arguments(x::SAnticommutator) = [x.op1, x.op2]
+operation(x::SAnticommutator) = anticommutator
+head(x::SAnticommutator) = :anticommutator
+children(x::SAnticommutator) = [:anticommutator, x.op1, x.op2]
+anticommutator(o1::Symbolic{AbstractOperator}, o2::Symbolic{AbstractOperator}) = SAnticommutator(o1, o2)
+Base.show(io::IO, x::SAnticommutator) = print(io, "{$(x.op1),$(x.op2)}")
+basis(x::SAnticommutator) = basis(x.op1)
+expand(x::SAnticommutator) = x == 0 ? x : x.op1*x.op2 + x.op2*x.op1