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.

Sign up for GitHub

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

Generalized CTranspose + conversion methods and promotion rules #33

Merged
merged 9 commits into from
May 20, 2015
Merged

Generalized CTranspose + conversion methods and promotion rules #33

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
dd09f6e
simplified similar_type usage by auto extracting the type from a valu…
jrevels May 14, 2015
730b968
generalized CTranspose and related conversion/promotion methods
jrevels May 14, 2015
a1a6e89
added space for clarity
jrevels May 14, 2015
1fcfc0f
generalized LabeledQuArray, bases descriptions now show properly for …
jrevels May 14, 2015
949482c
rm'd basis restriction for tensor funcs, added nary tensor test
jrevels May 19, 2015
3f20292
resolve ambiguity warnings by specifying basis supertype
jrevels May 19, 2015
08d3d24
rm'd parameterization of factor number from basis types, fixed tensor…
jrevels May 20, 2015
8e65edc
added unit tests for basic QuArray and CTranspose functions
jrevels May 20, 2015
114e25e
added type promotion tests
jrevels May 20, 2015
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
41 changes: 17 additions & 24 deletions src/arrays/arraymath.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,13 +27,13 @@ end
# operator * ket -> ket
function *{B<:OrthonormalBasis}(op::AbstractQuMatrix{B}, ket::AbstractQuVector{B})
vc = rawcoeffs(op)*rawcoeffs(ket)
QAT = similar_type(typeof(ket))
QAT = similar_type(ket)
return QAT(vc, bases(op,1))
end

function *{B<:OrthonormalBasis}(op::DualMatrix{B}, ket::AbstractQuVector{B})
vc = Ac_mul_B(rawcoeffs(op), rawcoeffs(ket))
QAT = similar_type(typeof(ket))
QAT = similar_type(ket)
return QAT(vc, bases(op,1))
end

Expand All @@ -50,29 +50,27 @@ end

function *{B<:OrthonormalBasis}(dm::DualMatrix{B}, qm::AbstractQuMatrix{B})
mc = Ac_mul_B(rawcoeffs(dm), rawcoeffs(qm))
QAT = similar_type(typeof(qm))
QAT = similar_type(qm)
return QAT(mc, bases(dm,1), bases(qm,2))
end

function *{B<:OrthonormalBasis}(qm::AbstractQuMatrix{B}, dm::DualMatrix{B})
mc = A_mul_Bc(rawcoeffs(dm), rawcoeffs(qm))
QAT = similar_type(typeof(qm))
QAT = similar_type(qm)
return QAT(mc, bases(qm,1), bases(dm,2))
end

function *{B<:OrthonormalBasis}(qm1::AbstractQuMatrix{B}, qm2::AbstractQuMatrix{B})
mc = rawcoeffs(qm1)*rawcoeffs(qm2)
QAT = similar_type(promote_type(typeof(qm1), typeof(qm2)))
QAT = similar_type(qm1, qm2)
return QAT(mc, bases(qm1,1), bases(qm2,2))
end



# addition and subtraction
function +{B<:AbstractBasis}(qarr1::AbstractQuArray{B}, qarr2::AbstractQuArray{B})
if bases(qarr1) == bases(qarr2)
sc = coeffs(qarr1) + coeffs(qarr2)
QAT = similar_type(promote_type(typeof(qarr1), typeof(qarr2)))
QAT = similar_type(qarr1, qarr2)
return QAT(sc, bases(qarr1))
else
error("Bases not compatible")
Expand All @@ -82,7 +80,7 @@ end
function -{B<:AbstractBasis}(qarr1::AbstractQuArray{B}, qarr2::AbstractQuArray{B})
if bases(qarr1) == bases(qarr2)
sc = coeffs(qarr1) - coeffs(qarr2)
QAT = similar_type(promote_type(typeof(qarr1), typeof(qarr2)))
QAT = similar_type(qarr1, qarr2)
return QAT(sc, bases(qarr1))
else
error("Bases not compatible")
Expand All @@ -92,15 +90,14 @@ end
+(ct1::CTranspose, ct2::CTranspose) = (ct1.qarr+ct2.qarr)'
-(ct1::CTranspose, ct2::CTranspose) = (ct1.qarr-ct2.qarr)'


# scaling
Base.scale!(num::Number, qarr::AbstractQuArray) = (scale!(num, rawcoeffs(qarr)); return qarr)
Base.scale!(num::Number, ct::CTranspose) = CTranspose(scale!(num', ct.qarr))
Base.scale!(qarr::AbstractQuArray, num::Number) = scale!(num, qarr)

function Base.scale(num::Number, qarr::AbstractQuArray)
fc = scale(num, rawcoeffs(qarr))
QAT = QuBase.similar_type(typeof(qarr))
QAT = similar_type(qarr)
return QAT(fc, bases(qarr))
end
Base.scale(num::Number, ct::CTranspose) = CTranspose(scale(num', ct.qarr))
Expand All @@ -110,7 +107,6 @@ Base.scale(qarr::AbstractQuArray, num::Number) = scale(num, qarr)
*(qarr::AbstractQuArray, num::Number) = scale(qarr, num)
/(qarr::AbstractQuArray, num::Number) = scale(1/num, qarr)


# matrix operations returning a scalar
# normalization
Base.norm(qarr::AbstractQuArray) = vecnorm(rawcoeffs(qarr))
Expand All @@ -122,48 +118,45 @@ end

normalize(qarr::AbstractQuArray) = normalize!(copy(qarr))


# matrix operations returning an array
# sparse to dense
function Base.full(qarr::AbstractQuMatrix)
fc = full(rawcoeffs(qarr))
QAT = similar_type(typeof(qarr))
QAT = similar_type(qarr)
return QAT(fc, bases(qarr))
end
#Base.full(ct::CTranspose) = full(ct.qarr)'

# exponential of dense matrix
function Base.expm(qarr::AbstractQuMatrix)
fc = expm(full(rawcoeffs(qarr)))
QAT = similar_type(typeof(qarr))
QAT = similar_type(qarr)
return QAT(fc, bases(qarr))
end
#Base.expm(ct::CTranspose) = expm(ct.qarr)'


##################
# Tensor Product #
##################

# General tensor product definitions for orthonormal bases
function tensor{B<:OrthonormalBasis,T1,T2,N}(qarr1::AbstractQuArray{B,T1,N}, qarr2::AbstractQuArray{B,T2,N})
function tensor{B,T1,T2,N}(qarr1::AbstractQuArray{B,T1,N}, qarr2::AbstractQuArray{B,T2,N})
tc = kron(coeffs(qarr1), coeffs(qarr2))
QAT = similar_type(promote_type(typeof(qarr1), typeof(qarr2)))
QAT = similar_type(qarr1, qarr2)
return QAT(tc, map(tensor, bases(qarr1), bases(qarr2)))
end

function tensor{B<:OrthonormalBasis,T1,T2,N}(qarr1::CTranspose{B,T1,N}, qarr2::CTranspose{B,T2,N})
return QuArray(kron(rawcoeffs(qarr1), rawcoeffs(qarr2)),
map(tensor, rawbases(qarr1), rawbases(qarr2)))'
end
# defined to resolve ambiguity warnings
tensor{B}(ct1::DualVector{B}, ct2::DualVector{B}) = tensor(ct1.qarr, ct2.qarr)'
tensor{B}(ct1::CTranspose{B}, ct2::CTranspose{B}) = tensor(ct1.qarr, ct2.qarr)'

function tensor{B<:OrthonormalBasis}(ket::QuVector{B}, bra::DualVector{B})
function tensor{B}(ket::AbstractQuVector{B}, bra::DualVector{B})
return QuArray(kron(coeffs(ket), coeffs(bra)),
Copy link
Contributor

Choose a reason for hiding this comment

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

I guess this one you forgot?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Good catch. I fixed it, and in the process of doing so uncovered an issue.

I found out that n-ary calls of tensor weren't dispatching properly because of the restriction that all arguments share a basis type. Since we parameterize the number of factors in our basis types, computing things like the following wouldn't work:

a = QuArray(rand(Complex{Float64},4), FiniteBasis(2,2))
b = QuArray(rand(Complex{Float64},1), FiniteBasis(1,1))
c = QuArray(rand(Complex{Float64},3))

ab = tensor(a,b)
tensor(ab, c) # fails, ab has a different basis type than c since the number of factors are different

After removing the basis type restriction, I added a test similar to the above. It's currently failing because type promotion for bases of different factor length is ambiguous:

julia> promote_type(typeof(ab), typeof(c))
ERROR: type: QuArray: in B, expected B<:AbstractBasis{S<:AbstractStructure}, got Type{FiniteBasis{S,N}}
 in promote_type at promotion.jl:110

I don't know what the right decision is here...perhaps we should consider abandon parameterizing factor length for basis types? Will any functions in the future require factor length as a parameter for type stability? Otherwise, we can just pull it from instances using the nfactors function (requiring it to be defined for all B<:AbstractBasis), and treat factor length as a value.

Copy link
Contributor

Choose a reason for hiding this comment

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

It feels like we had this discussion already? Anyways, at least at the moment we don't have anything that needs the factor length as a parameter, so we can just treat it as a value for now and switch back if necessary :-) I think when we introduced this I had in mind that doing partial traces/transpose could profit from it, but IMO this is to vague to stop this PR.

bases(ket,1),
bases(bra,1))
end

tensor{B<:OrthonormalBasis}(bra::DualVector{B}, ket::QuVector{B}) = tensor(ket, bra)
tensor{B}(bra::DualVector{B}, ket::AbstractQuVector{B}) = tensor(ket, bra)

###############
# Commutators #
Expand Down
1 change: 1 addition & 0 deletions src/arrays/constructors.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
QAT = similar_type(typeof(qa))
return QAT(fc, bases(qa))
end

function Base.eye(qa::AbstractQuArray)
fc = eye(coeffs(qa))
QAT = similar_type(typeof(qa))
Expand Down
62 changes: 39 additions & 23 deletions src/arrays/quarray.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,8 @@
typealias AbstractQuVector{B<:AbstractBasis,T} AbstractQuArray{B,T,1}
typealias AbstractQuMatrix{B<:AbstractBasis,T} AbstractQuArray{B,T,2}

similar_type{Q<:AbstractQuArray}(::Q) = similar_type(Q)
similar_type{A<:AbstractQuArray,B<:AbstractQuArray}(::A,::B) = similar_type(promote_type(A,B))

###########
# QuArray #
Expand Down Expand Up @@ -61,7 +63,6 @@
# and a tuple of bases (see QuArray(coeffs, bases) above)
similar_type{Q<:QuArray}(::Type{Q}) = QuArray


######################
# Accessor functions #
######################
Expand All @@ -74,30 +75,26 @@

Base.copy(qa::QuArray) = QuArray(copy(qa.coeffs), copy(qa.bases))


##############
# CTranspose #
##############
immutable CTranspose{B<:AbstractBasis,T,N,A} <: AbstractQuArray{B,T,N}
qarr::QuArray{B,T,N,A}
CTranspose(qarr::QuVector{B,T,A}) = new(qarr)
CTranspose(qarr::QuMatrix{B,T,A}) = new(qarr)
CTranspose(qarr::QuArray) = error("Conjugate transpose is unsupported for QuArrays of dimension $N")
immutable CTranspose{B<:AbstractBasis,T,N,Q} <: AbstractQuArray{B,T,N}
qarr::Q
CTranspose(qarr::AbstractQuVector{B,T}) = new(qarr)
CTranspose(qarr::AbstractQuMatrix{B,T}) = new(qarr)
CTranspose(qarr::AbstractQuArray) = error("Conjugate transpose is unsupported for AbstractQuArrays of dimension $N")
end

CTranspose{B<:AbstractBasis,T,N,A}(qa::QuArray{B,T,N,A}) = CTranspose{B,T,N,A}(qa)
CTranspose{B<:AbstractBasis,T,N}(qa::AbstractQuArray{B,T,N}) = CTranspose{B,T,N,typeof(qa)}(qa)
CTranspose{B<:AbstractBasis,T,N}(coeffs::AbstractArray{T,N}, bases::NTuple{N,B}) = CTranspose(QuArray(coeffs, bases))

typealias DualVector{B,T,A} CTranspose{B,T,1,A}
typealias DualMatrix{B,T,A} CTranspose{B,T,2,A}
typealias DualVector{B,T,Q} CTranspose{B,T,1,Q}
typealias DualMatrix{B,T,Q} CTranspose{B,T,2,Q}

similar_type{QC<:CTranspose}(::Type{QC}) = CTranspose


######################
# Accessor functions #
######################
coefftype{B,T,N,A}(::CTranspose{B,T,N,A}) = A
coefftype(ct::CTranspose) = coefftype(ct.qarr)

rawcoeffs(ct::CTranspose) = rawcoeffs(ct.qarr)
coeffs(ct::CTranspose) = rawcoeffs(ct)'
Expand All @@ -109,6 +106,10 @@

bases(ct::CTranspose, i) = rawbases(ct, revind(ndims(ct), i))

qarr_type{B,T,N,Q}(::Type{CTranspose{B,T,N,Q}}) = Q
qarr_type(ct::CTranspose) = qarr_type(typeof(ct))

similar_type{QC<:CTranspose}(::Type{QC}) = CTranspose

########################
# Array-like functions #
Expand All @@ -130,36 +131,51 @@
Base.ctranspose(qarr::QuArray) = CTranspose(qarr)
Base.ctranspose(ct::CTranspose) = ct.qarr

eager_ctranspose(qarr::AbstractQuArray) = similar_type(qarr)(coeffs(qarr)', reverse(bases(qarr)))

################
# LabelQuArray #
################
typealias LabelQuArray{B<:LabelBasis,T,N,A} QuArray{B,T,N,A}
typealias LabeledQuArray{B<:LabelBasis,T,N} AbstractQuArray{B,T,N}
typealias TupleArray{T<:Tuple,N} Array{T,N}

Base.getindex(larr::LabelQuArray, tups::Union(Tuple,TupleArray)...) = getindex(larr, map(getindex, bases(larr), tups)...)
Base.setindex!(larr::LabelQuArray, x, tups::Union(Tuple,TupleArray)...) = setindex!(larr, x, map(getindex, bases(larr), tups)...)
# redudant definitions added to resolve ambiguity warnings
Base.getindex{B<:LabelBasis}(larr::DualVector{B}, tups::Union(Tuple,TupleArray)) = getindex(larr, bases(larr,1)[tups])
Base.getindex{B<:LabelBasis}(larr::CTranspose{B}, tups::Union(Tuple,TupleArray)...) = getindex(larr, map(getindex, bases(larr), tups)...)
Base.getindex(larr::LabeledQuArray, tups::Union(Tuple,TupleArray)...) = getindex(larr, map(getindex, bases(larr), tups)...)

# redudant definitions added to resolve ambiguity warnings
Base.setindex!{B<:LabelBasis}(larr::DualVector{B}, x, tups::Union(Tuple,TupleArray)) = setindex(larr, x, bases(larr,1)[tups])
Base.setindex!{B<:LabelBasis}(larr::CTranspose{B}, x, tups::Union(Tuple,TupleArray)...) = setindex!(larr, x, map(getindex, bases(larr), tups)...)
Base.setindex!(larr::LabeledQuArray, x, tups::Union(Tuple,TupleArray)...) = setindex!(larr, x, map(getindex, bases(larr), tups)...)

###################
# Promotion rules #
###################
########################
# Conversion/Promotion #
########################
Base.promote_rule{B1,B2,T1,T2,N,A1,A2}(::Type{QuArray{B1,T1,N,A1}}, ::Type{QuArray{B2,T2,N,A2}}) = QuArray{promote_type(B1,B2),promote_type(T1,T2),N,promote_type(A1,A2)}
Base.promote_rule{C<:CTranspose,B,T,N,A}(::Type{C}, ::Type{QuArray{B,T,N,A}}) = promote_type(qarr_type(C), QuArray{B,T,N,A})
Base.promote_rule{DV<:DualVector, QV<:QuVector}(::Type{DV}, ::Type{QV}) = typejoin(DV,QV)

Base.promote_rule{B1,B2,T1,T2,N,A1,A2}(::Type{CTranspose{B1,T1,N,A1}}, ::Type{QuArray{B2,T2,N,A2}}) = QuArray{promote_type(B1,B2),promote_type(T1,T2),N,promote_type(A1,A2)}
Base.convert{B,T,N,A}(::Type{QuArray{B,T,N,A}}, qa::QuArray{B}) = QuArray(convert(A, coeffs(qa)), bases(qa))
Base.convert{B,T,N,A}(::Type{QuArray{B,T,N,A}}, qa::QuArray) = QuArray(convert(A, coeffs(qa)), map(i -> convert(B,i), bases(qa)))

Base.convert{C<:CTranspose,B,T,N,Q<:AbstractQuArray}(::Type{C}, ct::CTranspose{B,T,N,Q}) = CTranspose(convert(qarr_type(C),ct.qarr))
Base.convert{B,T,N,Q<:AbstractQuArray}(::Type{Q}, ct::CTranspose{B,T,N,Q}) = eager_ctranspose(ct.qarr)
Base.convert{QV<:AbstractQuVector, DV<:DualVector}(::Type{QV}, ::DV) = error("Cannot convert Bra to Ket. Use ctranspose if you would like to take the conjugate transpose.")
Base.convert{DV<:DualVector,B,T,QV<:AbstractQuVector}(::Type{DV}, ct::DualVector{B,T,QV}) = CTranspose(convert(qarr_type(DV),ct.qarr))

######################
# Printing Functions #
######################
Base.summary{B}(qarr::AbstractQuArray{B}) = "$(sizenotation(size(qarr))) $(typerepr(qarr)) in $B"
Base.summary(larr::LabelQuArray) = "$(sizenotation(size(larr))) $(typerepr(larr)) in labeled bases $(bases(larr))"
Base.summary(larr::LabeledQuArray) = "$(sizenotation(size(larr))) $(typerepr(larr)) in labeled bases $(bases(larr))"

function Base.show(io::IO, qarr::AbstractQuArray)
println(io, summary(qarr)*":")
println(io, "...coefficients: $(coefftype(qarr))")
print(io, repr(coeffs(qarr)))
end


####################
# Helper Functions #
####################
Expand Down