Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add qexpand functionality and cleanup scalings #59

Merged
merged 6 commits into from
Jul 2, 2024
Merged
Show file tree
Hide file tree
Changes from 4 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion .typos.toml
Original file line number Diff line number Diff line change
@@ -1,2 +1,3 @@
[default.extend-words]
ket = "ket"
ket = "ket"
BA = "BA"
1 change: 1 addition & 0 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ function main()
pages = [
"QuantumSymbolics.jl" => "index.md",
"Qubit Basis Choice" => "qubit_basis.md",
"Express Functionality" => "express.md",
"API" => "API.md",
]
)
Expand Down
48 changes: 18 additions & 30 deletions src/QSymbolicsBase/QSymbolicsBase.jl
Original file line number Diff line number Diff line change
@@ -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,
Expand All @@ -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₋ᵢ,
Expand All @@ -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
##
Expand Down Expand Up @@ -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")

Expand Down
32 changes: 15 additions & 17 deletions src/QSymbolicsBase/basic_ops_homogeneous.jl
Original file line number Diff line number Diff line change
Expand Up @@ -81,8 +81,8 @@
_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 = flattenop(+,[c*obj for (c,obj) in d])
length(d)==1 ? first(terms) : SAdd{S}(d,Set(terms),terms)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks a bit suspicious. There is no one-to-one correspondence between the dictionary, set, and list forms anymore. I think what you want is actually countmap_flatten as seen here

isempty(nonzero_terms) ? f : SAdd{T}(countmap_flatten(nonzero_terms, SScaled{T}))

Given that there is already flattening done in the + operation, I think we can just not do any flattening here. Just take the dictionary and do not worry about processing it. Any normal use of SAdd, usually through the use of + will already take care of flattening. This is a fairly general statement: fancy processing should be done in calls to + and *, not in constructors.

Also, something I did not notice previously, the variable names are misleading here. It should be (obj, coef) in d not (coef, obj) in d.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For this particular case, unlike the other constructors, we probably still want to keep the inner constructor that makes sure to create the set and vector.

end
isexpr(::SAdd) = true
iscall(::SAdd) = true
Expand All @@ -99,6 +99,11 @@
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

Check warning on line 105 in src/QSymbolicsBase/basic_ops_homogeneous.jl

View check run for this annotation

Codecov / codecov/patch

src/QSymbolicsBase/basic_ops_homogeneous.jl#L103-L105

Added lines #L103 - L105 were not covered by tests
end
const SAddKet = SAdd{AbstractKet}
function Base.show(io::IO, x::SAddKet)
ordered_terms = sort([repr(i) for i in arguments(x)])
Expand All @@ -109,11 +114,6 @@
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

Expand All @@ -128,7 +128,7 @@
terms
function SMulOperator(terms)
coeff, cleanterms = prefactorscalings(terms)
coeff*new(cleanterms)
coeff*new(flattenop(*,cleanterms))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not sure this is the best order of operations. You might have (c*A) * (c*((c*A)*(c*A))) -- for this particular example, flattening followed by prefactor extraction would work better. But I am not sure how general my statement is.

A second request: I know this state is kinda inherited already, but now that we are doing more and more fancy processing on construction, I want us to be a bit more disciplined about what goes in the constructor and what goes in the * operator. May we just remove the custom constructor and move all of this to *? Any reason that is a bad idea?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't see any reason why this would be a bad idea. My original thought process is that using inner constructor methods would be more idiomatic, but perhaps you're right in that it is too fancy for our purposes. It would also make the source code more clear to future developers if we do the processing in the calls to *, +, etc.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks! Please proceed with that plan. For more context: one thing that is pestering the back of my mind is how fancy constructors plagued the development of SymPy. At some point, for various reasons related to more advanced simplification rules, it became useful to be able to construct the most raw unsimplified versions of expressions, but all the __init__ methods were already fancy, so a lot of work had to be put in using another set of "pre-init" methods.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, that's interesting. Just out of curiosity, would there be any performance benefits in simplifying the constructors and moving the processing to conversion functions? I can't find any discussion about this in Julia discourse or on GitHub issue pages.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I do not think there would be anything (besides potentially avoiding situation where we happen to write repetitive canonicalization steps). Nothing inherent to the language though (constructors and other functions are treated the same way by the compiler, which is also why we are permitted to return object of type A from the constructor of type B (strictly speaking we are not really constructing type B anymore))

However, this separation of concerns we are discussing above might in some situation make the functions more "type stable" which permits the compiler to make a bunch of optimizations (mostly removing runtime type checks and allocations).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, this is good to know.

end
end
isexpr(::SMulOperator) = true
Expand All @@ -155,22 +155,22 @@
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)
coeff * new{S}(flattenop(⊗,cleanterms))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same question about order of operations and moving this to

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]

Check warning on line 173 in src/QSymbolicsBase/basic_ops_homogeneous.jl

View check run for this annotation

Codecov / codecov/patch

src/QSymbolicsBase/basic_ops_homogeneous.jl#L173

Added line #L173 was not covered by tests
function ⊗(xs::Symbolic{T}...) where {T<:QObj}
zero_ind = findfirst(x->iszero(x), xs)
isnothing(zero_ind) ? STensor{T}(collect(xs)) : SZero{T}()
Expand All @@ -182,9 +182,9 @@
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)),"⊗")*")")

Check warning on line 187 in src/QSymbolicsBase/basic_ops_homogeneous.jl

View check run for this annotation

Codecov / codecov/patch

src/QSymbolicsBase/basic_ops_homogeneous.jl#L187

Added line #L187 was not covered by tests

"""Symbolic commutator of two operators

Expand All @@ -202,7 +202,7 @@
op1
op2
function SCommutator(o1, o2)
coeff, cleanterms = prefactorscalings([o1 o2], scalar=true)
coeff, cleanterms = prefactorscalings([o1 o2])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same question about dropping the custom constructor and moving this to commutator

cleanterms[1] === cleanterms[2] ? SZeroOperator() : coeff*new(cleanterms...)
end
end
Expand All @@ -218,7 +218,6 @@
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

Expand All @@ -233,7 +232,7 @@
op1
op2
function SAnticommutator(o1, o2)
coeff, cleanterms = prefactorscalings([o1 o2], scalar=true)
coeff, cleanterms = prefactorscalings([o1 o2])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same question about dropping the custom constructor and moving this to anticommutator

coeff*new(cleanterms...)
end
end
Expand All @@ -249,4 +248,3 @@
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
2 changes: 2 additions & 0 deletions src/QSymbolicsBase/basic_ops_inhomogeneous.jl
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,8 @@
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)

Check warning on line 89 in src/QSymbolicsBase/basic_ops_inhomogeneous.jl

View check run for this annotation

Codecov / codecov/patch

src/QSymbolicsBase/basic_ops_inhomogeneous.jl#L89

Added line #L89 was not covered by tests
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"""
Expand Down
117 changes: 117 additions & 0 deletions src/QSymbolicsBase/linalg.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,117 @@
##
# Linear algebra operations on quantum objects.
##

"""Projector for a given ket

```jldoctest
julia> SProjector(X1⊗X2)
𝐏[|X₁⟩|X₂⟩]

julia> express(SProjector(X2))
Operator(dim=2x2)
basis: Spin(1/2)
0.5+0.0im -0.5-0.0im
-0.5+0.0im 0.5+0.0im
```"""
@withmetadata struct SProjector <: Symbolic{AbstractOperator}
ket::Symbolic{AbstractKet} # TODO parameterize
end
isexpr(::SProjector) = true
iscall(::SProjector) = true

Check warning on line 21 in src/QSymbolicsBase/linalg.jl

View check run for this annotation

Codecov / codecov/patch

src/QSymbolicsBase/linalg.jl#L21

Added line #L21 was not covered by tests
arguments(x::SProjector) = [x.ket]
operation(x::SProjector) = projector
head(x::SProjector) = :projector
children(x::SProjector) = [:projector,x.ket]
projector(x::Symbolic{AbstractKet}) = SProjector(x)

Check warning on line 26 in src/QSymbolicsBase/linalg.jl

View check run for this annotation

Codecov / codecov/patch

src/QSymbolicsBase/linalg.jl#L25-L26

Added lines #L25 - L26 were not covered by tests
projector(x::SZeroKet) = SZeroOperator()
basis(x::SProjector) = basis(x.ket)
function Base.show(io::IO, x::SProjector)
print(io,"𝐏[")
print(io,x.ket)
print(io,"]")
end

"""Dagger, i.e., adjoint of quantum objects (kets, bras, operators)

```jldoctest
julia> @ket a; @op A;

julia> dagger(2*im*A*a)
(0 - 2im)|a⟩†A†

julia> @op B;

julia> dagger(A*B)
B†A†

julia> ℋ = SHermitianOperator(:ℋ); U = SUnitaryOperator(:U);

julia> dagger(ℋ)

julia> dagger(U)
U⁻¹
```
"""
@withmetadata struct SDagger{T<:QObj} <: Symbolic{T}
obj
end
isexpr(::SDagger) = true
iscall(::SDagger) = true

Check warning on line 61 in src/QSymbolicsBase/linalg.jl

View check run for this annotation

Codecov / codecov/patch

src/QSymbolicsBase/linalg.jl#L61

Added line #L61 was not covered by tests
arguments(x::SDagger) = [x.obj]
operation(x::SDagger) = dagger
head(x::SDagger) = :dagger
children(x::SDagger) = [:dagger, x.obj]

Check warning on line 65 in src/QSymbolicsBase/linalg.jl

View check run for this annotation

Codecov / codecov/patch

src/QSymbolicsBase/linalg.jl#L65

Added line #L65 was not covered by tests
dagger(x::Symbolic{AbstractBra}) = SDagger{AbstractKet}(x)
dagger(x::Symbolic{AbstractKet}) = SDagger{AbstractBra}(x)
dagger(x::Symbolic{AbstractOperator}) = SDagger{AbstractOperator}(x)
dagger(x::SScaledKet) = SScaledBra(conj(x.coeff), dagger(x.obj))
dagger(x::SAdd) = (+)((dagger(i) for i in arguments(x))...)
dagger(x::SScaledBra) = SScaledKet(conj(x.coeff), dagger(x.obj))
dagger(x::SZeroOperator) = x
dagger(x::SHermitianOperator) = x
dagger(x::SHermitianUnitaryOperator) = x

Check warning on line 74 in src/QSymbolicsBase/linalg.jl

View check run for this annotation

Codecov / codecov/patch

src/QSymbolicsBase/linalg.jl#L74

Added line #L74 was not covered by tests
dagger(x::SUnitaryOperator) = inv(x)
dagger(x::STensorBra) = STensorKet(collect(dagger(i) for i in x.terms))
dagger(x::STensorKet) = STensorBra(collect(dagger(i) for i in x.terms))
dagger(x::STensorOperator) = STensorOperator(collect(dagger(i) for i in x.terms))
dagger(x::SScaledOperator) = SScaledOperator(conj(x.coeff), dagger(x.obj))
dagger(x::SApplyKet) = dagger(x.ket)*dagger(x.op)
dagger(x::SApplyBra) = dagger(x.op)*dagger(x.bra)
dagger(x::SMulOperator) = SMulOperator(collect(dagger(i) for i in reverse(x.terms)))
dagger(x::SBraKet) = SBraKet(dagger(x.ket), dagger(x.bra))
dagger(x::SOuterKetBra) = SOuterKetBra(dagger(x.bra), dagger(x.ket))
dagger(x::SDagger) = x.obj
basis(x::SDagger) = basis(x.obj)

Check warning on line 86 in src/QSymbolicsBase/linalg.jl

View check run for this annotation

Codecov / codecov/patch

src/QSymbolicsBase/linalg.jl#L86

Added line #L86 was not covered by tests
function Base.show(io::IO, x::SDagger)
print(io,x.obj)
print(io,"†")
end

"""Inverse Operator

```jldoctest
julia> @op A;

julia> inv(A)
A⁻¹

julia> inv(A)*A
𝕀
```
"""
@withmetadata struct SInvOperator <: Symbolic{AbstractOperator}
op::Symbolic{AbstractOperator}
end
isexpr(::SInvOperator) = true
iscall(::SInvOperator) = true

Check warning on line 108 in src/QSymbolicsBase/linalg.jl

View check run for this annotation

Codecov / codecov/patch

src/QSymbolicsBase/linalg.jl#L108

Added line #L108 was not covered by tests
arguments(x::SInvOperator) = [x.op]
operation(x::SInvOperator) = inv
head(x::SInvOperator) = :inv
children(x::SInvOperator) = [:inv, x.op]
basis(x::SInvOperator) = basis(x.op)

Check warning on line 113 in src/QSymbolicsBase/linalg.jl

View check run for this annotation

Codecov / codecov/patch

src/QSymbolicsBase/linalg.jl#L111-L113

Added lines #L111 - L113 were not covered by tests
Base.show(io::IO, x::SInvOperator) = print(io, "$(x.op)⁻¹")
Base.:(*)(invop::SInvOperator, op::SOperator) = isequal(invop.op, op) ? IdentityOp(basis(op)) : SMulOperator(invop, op)
Base.:(*)(op::SOperator, invop::SInvOperator) = isequal(op, invop.op) ? IdentityOp(basis(op)) : SMulOperator(op, invop)

Check warning on line 116 in src/QSymbolicsBase/linalg.jl

View check run for this annotation

Codecov / codecov/patch

src/QSymbolicsBase/linalg.jl#L116

Added line #L116 was not covered by tests
inv(x::Symbolic{AbstractOperator}) = SInvOperator(x)
2 changes: 1 addition & 1 deletion src/QSymbolicsBase/literal_objects.jl
Original file line number Diff line number Diff line change
Expand Up @@ -92,4 +92,4 @@ symbollabel(x::SZero) = "𝟎"
basis(x::SZero) = nothing

Base.show(io::IO, x::SZero) = print(io, symbollabel(x))
Base.iszero(x::SymQObj) = isa(x, SZero)
Base.iszero(x::Union{SymQObj, Symbolic{Number}, Symbolic{Complex}}) = isa(x, SZero)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is piracy. Why is this needed?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My bad about the piracy. It was intended for when we process * for SScaled objects, although it's kinda of a hackish way to deal with this issue. Basically, I was getting errors before when I would multiply an inner product (which is a subtype of Symbolic{Complex}) or a real number plus an inner product (subtype of Symbolic{Number}), as iszero was not defined for SBraKet. I'll fix this.

Loading
Loading