Add qexpand functionality and cleanup scalings (#59)

* utils cleanup

* add qexpand and tests

* Update tests

* create express doc page

* convert inner constructors to conversion functions

* update changelog.md and project.toml

.typos.toml

@@ -1,2 +1,3 @@
 [default.extend-words]
-ket = "ket"
+ket = "ket"
+BA = "BA"

CHANGELOG.md

@@ -1,5 +1,11 @@
 # News
 
+## v0.3.3 - 2024-07-02
+
+- Introduced `qexpand` function that manually expands expressions containing quantum objects.
+- Organized automatic scaling and flattening procedures.
+- Added `express.md` to docs.
+
 ## v0.3.2 - 2024-06-28
 
 - Added documentation for `express`.

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumSymbolics"
 uuid = "efa7fd63-0460-4890-beb7-be1bbdfbaeae"
 authors = ["QuantumSymbolics.jl contributors"]
-version = "0.3.2"
+version = "0.3.3"
 
 [deps]
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"

docs/make.jl

@@ -26,6 +26,7 @@ function main()
     pages = [
         "QuantumSymbolics.jl" => "index.md",
         "Qubit Basis Choice" => "qubit_basis.md",
+        "Express Functionality" => "express.md",
         "API" => "API.md",
     ]
     )

src/QSymbolicsBase/QSymbolicsBase.jl

@@ -1,20 +1,20 @@
 using Symbolics
 import Symbolics: simplify
 using SymbolicUtils
-import SymbolicUtils: Symbolic, _isone, flatten_term, isnotflat, Chain, Fixpoint, Prewalk
+import SymbolicUtils: Symbolic,_isone,flatten_term,isnotflat,Chain,Fixpoint,Prewalk
 using TermInterface
-import TermInterface: isexpr, head, iscall, children, operation, arguments, metadata, maketerm
+import TermInterface: isexpr,head,iscall,children,operation,arguments,metadata,maketerm
 
 using LinearAlgebra
-import LinearAlgebra: eigvecs, ishermitian, inv
+import LinearAlgebra: eigvecs,ishermitian,inv
 
 import QuantumInterface:
     apply!,
     tensor, ⊗,
-    basis, Basis, SpinBasis, FockBasis,
+    basis,Basis,SpinBasis,FockBasis,
     nqubits,
-    projector, dagger,
-    AbstractKet, AbstractOperator, AbstractSuperOperator, AbstractBra
+    projector,dagger,
+    AbstractBra,AbstractKet,AbstractOperator,AbstractSuperOperator
 
 export SymQObj,QObj,
        AbstractRepresentation,AbstractUse,
@@ -23,7 +23,7 @@ export SymQObj,QObj,
        apply!,
        express,
        tensor,⊗,
-       dagger,projector,commutator,anticommutator,expand,
+       dagger,projector,commutator,anticommutator,
        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₋ᵢ,
@@ -35,36 +35,16 @@ export SymQObj,QObj,
        SScaled,SScaledBra,SScaledOperator,SScaledKet,
        STensorBra,STensorKet,STensorOperator,
        SZeroBra,SZeroKet,SZeroOperator,
-       SProjector,MixedState,IdentityOp,SInvOperator,SHermitianOperator,SUnitaryOperator,SHermitianUnitaryOperator,
+       SProjector,MixedState,IdentityOp,SInvOperator,
        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,
+       qsimplify,qsimplify_pauli,qsimplify_commutator,qsimplify_anticommutator,
+       qexpand,
        isunitary
 
-function countmap(samples) # A simpler version of StatsBase.countmap, because StatsBase is slow to import
-    counts = Dict{Any,Any}()
-    for s in samples
-        counts[s] = get(counts, s, 0)+1
-    end
-    counts
-end
-
-function countmap_flatten(samples, flattenhead)
-    counts = Dict{Any,Any}()
-    for s in samples
-        if isexpr(s) && s isa flattenhead # TODO Could you use the TermInterface `operation` here instead of `flattenhead`?
-            coef, term = arguments(s)
-            counts[term] = get(counts, term, 0)+coef
-        else
-            counts[s] = get(counts, s, 0)+1
-        end
-    end
-    counts
-end
-
 ##
 # Metadata cache helpers
 ##
@@ -167,13 +147,21 @@ Base.isequal(::Symbolic{Complex}, ::SymQObj) = false
 # use a macro to provide specializations if that is indeed the case
 propsequal(x,y) = all(n->isequal(getproperty(x,n),getproperty(y,n)), propertynames(x))
 
+
+##
+# Utilities
+##
+
+include("utils.jl")
+
 ##
 # Most symbolic objects defined here
 ##
 
 include("literal_objects.jl")
 include("basic_ops_homogeneous.jl")
 include("basic_ops_inhomogeneous.jl")
+include("linalg.jl")
 include("predefined.jl")
 include("predefined_CPTP.jl")
 

src/QSymbolicsBase/basic_ops_homogeneous.jl

@@ -22,7 +22,6 @@ julia> 2*A
 @withmetadata struct SScaled{T<:QObj} <: Symbolic{T}
     coeff
     obj
-    SScaled{S}(c,k) where S = _isone(c) ? k : new{S}(c,k)
 end
 isexpr(::SScaled) = true
 iscall(::SScaled) = true
@@ -31,10 +30,14 @@ operation(x::SScaled) = *
 head(x::SScaled) = :*
 children(x::SScaled) = [:*,x.coeff,x.obj]
 function Base.:(*)(c, x::Symbolic{T}) where {T<:QObj} 
-    if iszero(c) || iszero(x)
+    if (isa(c, Number) && iszero(c)) || iszero(x)
         SZero{T}()
+    elseif _isone(c)
+        x
+    elseif isa(x, SScaled)
+        SScaled{T}(c*x.coeff, x.obj)
     else 
-        x isa SScaled ? SScaled{T}(c*x.coeff, x.obj) : SScaled{T}(c, x) 
+        SScaled{T}(c, x) 
     end
 end
 Base.:(*)(x::Symbolic{T}, c) where {T<:QObj} = c*x
@@ -81,8 +84,8 @@ julia> k₁ + k₂
     _arguments_precomputed
 end
 function SAdd{S}(d) where S 
-    xs = [c*obj for (c,obj) in d]
-    length(d)==1 ? first(xs) : SAdd{S}(d,Set(xs),xs)
+    terms = [c*obj for (obj,c) in d]
+    length(d)==1 ? first(terms) : SAdd{S}(d,Set(terms),terms)
 end
 isexpr(::SAdd) = true
 iscall(::SAdd) = true
@@ -99,6 +102,11 @@ end
 Base.:(+)(xs::Vararg{Symbolic{<:QObj},0}) = 0 # to avoid undefined type parameters issue in the above method
 basis(x::SAdd) = basis(first(x.dict).first)
 
+const SAddBra = SAdd{AbstractBra}
+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
 const SAddKet = SAdd{AbstractKet}
 function Base.show(io::IO, x::SAddKet)
     ordered_terms = sort([repr(i) for i in arguments(x)])
@@ -109,11 +117,6 @@ 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}
-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
 
@@ -126,10 +129,6 @@ 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
@@ -139,7 +138,13 @@ head(x::SMulOperator) = :*
 children(x::SMulOperator) = [:*;x.terms]
 function Base.:(*)(xs::Symbolic{AbstractOperator}...) 
     zero_ind = findfirst(x->iszero(x), xs)
-    isnothing(zero_ind) ? SMulOperator(collect(xs)) : SZeroOperator()
+    if isnothing(zero_ind)
+        terms = flattenop(*, collect(xs))
+        coeff, cleanterms = prefactorscalings(terms)
+        coeff * SMulOperator(cleanterms)
+    else
+        SZeroOperator()
+    end
 end
 Base.show(io::IO, x::SMulOperator) = print(io, join(map(string, arguments(x)),""))
 basis(x::SMulOperator) = basis(x.terms)
@@ -155,25 +160,27 @@ julia> k₁ ⊗ k₂
 julia> @op A; @op B;
 
 julia> A ⊗ B 
-A⊗B
+(A⊗B)
 ```
 """
 @withmetadata struct STensor{T<:QObj} <: Symbolic{T}
     terms
-    function STensor{S}(terms) where S
-        coeff, cleanterms = prefactorscalings(terms)
-        coeff * new{S}(cleanterms)
-    end
 end
 isexpr(::STensor) = true
 iscall(::STensor) = true
 arguments(x::STensor) = x.terms
 operation(x::STensor) = ⊗
 head(x::STensor) = :⊗
-children(x::STensor) = pushfirst!(x.terms,:⊗)
+children(x::STensor) = [:⊗; x.terms]
 function ⊗(xs::Symbolic{T}...) where {T<:QObj}
     zero_ind = findfirst(x->iszero(x), xs)
-    isnothing(zero_ind) ? STensor{T}(collect(xs)) : SZero{T}()
+    if isnothing(zero_ind)
+        terms = flattenop(⊗, collect(xs))
+        coeff, cleanterms = prefactorscalings(terms)
+        coeff * STensor{T}(cleanterms)
+    else
+        SZero{T}()
+    end
 end
 basis(x::STensor) = tensor(basis.(x.terms)...)
 
@@ -182,9 +189,9 @@ Base.show(io::IO, x::STensorBra) = print(io, join(map(string, arguments(x)),""))
 const STensorKet = STensor{AbstractKet}
 Base.show(io::IO, x::STensorKet) = print(io, join(map(string, arguments(x)),""))
 const STensorOperator = STensor{AbstractOperator}
-Base.show(io::IO, x::STensorOperator) = print(io, join(map(string, arguments(x)),"⊗"))
+Base.show(io::IO, x::STensorOperator) = print(io, "("*join(map(string, arguments(x)),"⊗")*")")
 const STensorSuperOperator = STensor{AbstractSuperOperator}
-Base.show(io::IO, x::STensorSuperOperator) = print(io, join(map(string, arguments(x)),"⊗"))
+Base.show(io::IO, x::STensorSuperOperator) = print(io, "("*join(map(string, arguments(x)),"⊗")*")")
 
 """Symbolic commutator of two operators
 
@@ -201,24 +208,22 @@ julia> commutator(A, A)
 @withmetadata struct SCommutator <: Symbolic{AbstractOperator}
     op1
     op2
-    function SCommutator(o1, o2) 
-        coeff, cleanterms = prefactorscalings([o1 o2], scalar=true)
-        cleanterms[1] === cleanterms[2] ? SZeroOperator() : 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)
+function commutator(o1::Symbolic{AbstractOperator}, o2::Symbolic{AbstractOperator})
+    coeff, cleanterms = prefactorscalings([o1 o2])
+    cleanterms[1] === cleanterms[2] ? SZeroOperator() : coeff * SCommutator(cleanterms...)   
+end
 commutator(o1::SZeroOperator, o2::Symbolic{AbstractOperator}) = SZeroOperator()
 commutator(o1::Symbolic{AbstractOperator}, o2::SZeroOperator) = SZeroOperator()
 commutator(o1::SZeroOperator, o2::SZeroOperator) = SZeroOperator()
 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
 
@@ -232,21 +237,19 @@ julia> anticommutator(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)
+function anticommutator(o1::Symbolic{AbstractOperator}, o2::Symbolic{AbstractOperator})
+    coeff, cleanterms = prefactorscalings([o1 o2])
+    coeff * SAnticommutator(cleanterms...)
+end
 anticommutator(o1::SZeroOperator, o2::Symbolic{AbstractOperator}) = SZeroOperator()
 anticommutator(o1::Symbolic{AbstractOperator}, o2::SZeroOperator) = SZeroOperator()
 anticommutator(o1::SZeroOperator, o2::SZeroOperator) = SZeroOperator()
 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

src/QSymbolicsBase/basic_ops_inhomogeneous.jl

@@ -86,6 +86,8 @@ Base.:(*)(b::SZeroBra, k::Symbolic{AbstractKet}) = 0
 Base.:(*)(b::Symbolic{AbstractBra}, k::SZeroKet) = 0
 Base.:(*)(b::SZeroBra, k::SZeroKet) = 0
 Base.show(io::IO, x::SBraKet) = begin print(io,x.bra); print(io,x.ket) end
+Base.hash(x::SBraKet, h::UInt) = hash((head(x), arguments(x)), h)
+maketerm(::Type{<:SBraKet}, f, a, t, m) = f(a...)
 Base.isequal(x::SBraKet, y::SBraKet) = isequal(x.bra, y.bra) && isequal(x.ket, y.ket)
 
 """Symbolic application of a superoperator on an operator"""