diff --git a/src/arrays/arraymath.jl b/src/arrays/arraymath.jl index 4d6d48c..9865b01 100644 --- a/src/arrays/arraymath.jl +++ b/src/arrays/arraymath.jl @@ -1,4 +1,4 @@ -import Base: *, +, - +import Base: *, +, -, / ################## # Multiplication # @@ -11,12 +11,12 @@ import Base: *, +, - typealias OrthonormalBasis{S<:Orthonormal} AbstractBasis{S} # bra * ket -> scalar -function *{B<:OrthonormalBasis}(bra::DualVector{B}, ket::QuVector{B}) +function *{B<:OrthonormalBasis}(bra::DualVector{B}, ket::AbstractQuVector{B}) return dot(rawcoeffs(bra), rawcoeffs(ket)) end # bra * operator -> bra -function *{B<:OrthonormalBasis}(bra::DualVector{B}, op::QuMatrix{B}) +function *{B<:OrthonormalBasis}(bra::DualVector{B}, op::AbstractQuMatrix{B}) return QuArray(Ac_mul_B(rawcoeffs(op), rawcoeffs(bra)), bases(op,2))' end @@ -25,12 +25,16 @@ function *{B<:OrthonormalBasis}(bra::DualVector{B}, op::DualMatrix{B}) end # operator * ket -> ket -function *{B<:OrthonormalBasis}(op::QuMatrix{B}, ket::QuVector{B}) - return QuArray(rawcoeffs(op)*rawcoeffs(ket), bases(op,1)) +function *{B<:OrthonormalBasis}(op::AbstractQuMatrix{B}, ket::AbstractQuVector{B}) + vc = rawcoeffs(op)*rawcoeffs(ket) + QAT = similar_type(typeof(ket)) + return QAT(vc, bases(op,1)) end -function *{B<:OrthonormalBasis}(op::DualMatrix{B}, ket::QuVector{B}) - return QuArray(Ac_mul_B(rawcoeffs(op), rawcoeffs(ket)), bases(op,1)) +function *{B<:OrthonormalBasis}(op::DualMatrix{B}, ket::AbstractQuVector{B}) + vc = Ac_mul_B(rawcoeffs(op), rawcoeffs(ket)) + QAT = similar_type(typeof(ket)) + return QAT(vc, bases(op,1)) end # ket * bra -> operator @@ -41,70 +45,101 @@ function *{B<:OrthonormalBasis}(ket::QuVector{B}, bra::DualVector{B}) end # operator * operator -> operator -function *{B<:OrthonormalBasis}(qm::QuMatrix{B}, dm::DualMatrix{B}) - return QuArray(A_mul_Bc(rawcoeffs(qm), rawcoeffs(dm)), - bases(qm,1), - bases(dm,2)) +*(dm1::DualMatrix, dm2::DualMatrix) = (dm2.qarr*dm1.qarr)' +*{B<:OrthonormalBasis}(dm1::DualMatrix{B}, dm2::DualMatrix{B}) = (dm2.qarr*dm1.qarr)' + +function *{B<:OrthonormalBasis}(dm::DualMatrix{B}, qm::AbstractQuMatrix{B}) + mc = Ac_mul_B(rawcoeffs(dm), rawcoeffs(qm)) + QAT = similar_type(typeof(qm)) + return QAT(mc, bases(dm,1), bases(qm,2)) end -function *{B<:OrthonormalBasis}(dm::DualMatrix{B}, qm::QuMatrix{B}) - return QuArray(Ac_mul_B(rawcoeffs(dm), rawcoeffs(qm)), - bases(dm,1), - bases(qm,2)) +function *{B<:OrthonormalBasis}(qm::AbstractQuMatrix{B}, dm::DualMatrix{B}) + mc = A_mul_Bc(rawcoeffs(dm), rawcoeffs(qm)) + QAT = similar_type(typeof(qm)) + return QAT(mc, bases(qm,1), bases(dm,2)) end -function *{B<:OrthonormalBasis}(qm1::QuMatrix{B}, qm2::QuMatrix{B}) - return QuArray(rawcoeffs(qm1)*rawcoeffs(qm2), - bases(qm1,1), - bases(qm2,2)) +function *{B<:OrthonormalBasis}(qm1::AbstractQuMatrix{B}, qm2::AbstractQuMatrix{B}) + mc = rawcoeffs(qm1)*rawcoeffs(qm2) + QAT = similar_type(promote_type(typeof(qm1), typeof(qm2))) + return QAT(mc, bases(qm1,1), bases(qm2,2)) end -*(dm1::DualMatrix, dm2::DualMatrix) = (dm2.qarr*dm1.qarr)' -function +(qarr1::AbstractQuArray, qarr2::AbstractQuArray) + +# addition and subtraction +function +{B<:AbstractBasis}(qarr1::AbstractQuArray{B}, qarr2::AbstractQuArray{B}) if bases(qarr1) == bases(qarr2) - return QuArray(coeffs(qarr1)+coeffs(qarr2), bases(qarr1)) + sc = coeffs(qarr1) + coeffs(qarr2) + QAT = similar_type(promote_type(typeof(qarr1), typeof(qarr2))) + return QAT(sc, bases(qarr1)) else error("Bases not compatible") end end -function -(qarr1::AbstractQuArray, qarr2::AbstractQuArray) +function -{B<:AbstractBasis}(qarr1::AbstractQuArray{B}, qarr2::AbstractQuArray{B}) if bases(qarr1) == bases(qarr2) - return QuArray(coeffs(qarr1)-coeffs(qarr2), bases(qarr1)) + sc = coeffs(qarr1) - coeffs(qarr2) + QAT = similar_type(promote_type(typeof(qarr1), typeof(qarr2))) + return QAT(sc, bases(qarr1)) else error("Bases not compatible") end end ++(ct1::CTranspose, ct2::CTranspose) = (ct1.qarr+ct2.qarr)' +-(ct1::CTranspose, ct2::CTranspose) = (ct1.qarr-ct2.qarr)' + + # scaling -Base.scale!(num::Number, qarr::QuArray) = (scale!(num, rawcoeffs(qarr)); return qarr) +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::Union(QuArray,CTranspose), num::Number) = scale!(num, qarr) +Base.scale!(qarr::AbstractQuArray, num::Number) = scale!(num, qarr) -Base.scale(num::Number, qarr::QuArray) = QuArray(scale(num, rawcoeffs(qarr)), rawbases(qarr)) +function Base.scale(num::Number, qarr::AbstractQuArray) + fc = scale(num, rawcoeffs(qarr)) + QAT = QuBase.similar_type(typeof(qarr)) + return QAT(fc, bases(qarr)) +end Base.scale(num::Number, ct::CTranspose) = CTranspose(scale(num', ct.qarr)) -Base.scale(qarr::Union(QuArray,CTranspose), num::Number) = scale(num, qarr) +Base.scale(qarr::AbstractQuArray, num::Number) = scale(num, qarr) -*(num::Number, qarr::Union(QuArray,CTranspose)) = scale(num, qarr) -*(qarr::Union(QuArray,CTranspose), num::Number) = scale(qarr, num) -/(qarr::Union(QuArray,CTranspose), num::Number) = scale(1/num, qarr) - -# sparse to dense -Base.full(qarr::AbstractQuMatrix) = QuArray(full(coeffs(qarr)),bases(qarr)) +*(num::Number, qarr::AbstractQuArray) = scale(num, qarr) +*(qarr::AbstractQuArray, num::Number) = scale(qarr, num) +/(qarr::AbstractQuArray, num::Number) = scale(1/num, qarr) -# exponential of dense matrix -Base.expm(qarr::AbstractQuMatrix) = QuArray(expm(full(coeffs(qarr))),bases(qarr)) +# matrix operations returning a scalar # normalization -Base.norm(qarr::Union(QuArray,CTranspose)) = vecnorm(rawcoeffs(qarr)) +Base.norm(qarr::AbstractQuArray) = vecnorm(rawcoeffs(qarr)) -function normalize!(qarr::Union(QuArray,CTranspose)) +function normalize!(qarr::AbstractQuArray) scale!(1/norm(qarr), rawcoeffs(qarr)) return qarr end -normalize(qarr::Union(QuArray,CTranspose)) = normalize!(copy(qarr)) +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)) + 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)) + return QAT(fc, bases(qarr)) +end +#Base.expm(ct::CTranspose) = expm(ct.qarr)' + ################## # Tensor Product # @@ -112,8 +147,9 @@ normalize(qarr::Union(QuArray,CTranspose)) = normalize!(copy(qarr)) # General tensor product definitions for orthonormal bases function tensor{B<:OrthonormalBasis,T1,T2,N}(qarr1::AbstractQuArray{B,T1,N}, qarr2::AbstractQuArray{B,T2,N}) - return QuArray(kron(coeffs(qarr1), coeffs(qarr2)), - map(tensor, bases(qarr1), bases(qarr2))) + tc = kron(coeffs(qarr1), coeffs(qarr2)) + QAT = similar_type(promote_type(typeof(qarr1), typeof(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}) diff --git a/src/arrays/constructors.jl b/src/arrays/constructors.jl index bcbebdf..ad93181 100644 --- a/src/arrays/constructors.jl +++ b/src/arrays/constructors.jl @@ -7,8 +7,16 @@ ########################### # Array-like Constructors # ########################### - Base.zeros(qa::AbstractQuArray) = QuArray(zeros(coeffs(qa)), bases(qa)) - Base.eye(qa::AbstractQuArray) = QuArray(eye(coeffs(qa)), bases(qa)) + function Base.zeros(qa::AbstractQuArray) + fc = zeros(coeffs(qa)) + 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)) + return QAT(fc, bases(qa)) + end #################### # Helper Functions # @@ -21,7 +29,7 @@ function coherentstatevec_inf(n::Int,alpha::Number) s = zeros(typeof(float(alpha)),n) s[1] = one(alpha) - for i in [2:n] + for i in 2:n s[i] = alpha/sqrt(i-1)*s[i-1] end z = QuArray(s) diff --git a/src/arrays/ladderops.jl b/src/arrays/ladderops.jl index e1587ee..594fd6c 100644 --- a/src/arrays/ladderops.jl +++ b/src/arrays/ladderops.jl @@ -16,11 +16,11 @@ lowermatrix(lens, particle) = eye_sandwich(lens, particle, lower_single(lens[particle])) laddercoeffs(n) = sqrt(linspace(1, n, n)) - lower_single(n) = sparse([1:n-1], [2:n], laddercoeffs(n-1), n, n) - raise_single(n) = sparse([2:n], [1:n-1], laddercoeffs(n-1), n, n) + lower_single(n) = sparse([1:n-1;], [2:n;], laddercoeffs(n-1), n, n) + raise_single(n) = sparse([2:n;], [1:n-1;], laddercoeffs(n-1), n, n) - before_eye(lens, pivot) = speye(prod(lens[[1:pivot-1]])) - after_eye(lens, pivot) = speye(prod(lens[[pivot+1:length(lens)]])) + before_eye(lens, pivot) = speye(prod(lens[1:pivot-1])) + after_eye(lens, pivot) = speye(prod(lens[pivot+1:length(lens)])) function eye_sandwich(lens, pivot, op) return kron(before_eye(lens, pivot), op, diff --git a/src/arrays/quarray.jl b/src/arrays/quarray.jl index 471959b..3c57ca1 100644 --- a/src/arrays/quarray.jl +++ b/src/arrays/quarray.jl @@ -1,11 +1,42 @@ -########### -# QuArray # -########### +################### +# AbstractQuArray # +################### abstract AbstractQuArray{B<:AbstractBasis,T,N} + # all subtypes of AbstractQuArray have to implement the following functions + # - coefftype + # - rawcoeffs + # - rawbases + # - copy + # - similar_type + + coeffs(qarr::AbstractQuArray) = rawcoeffs(qarr) + + # works generally as long as the single index form is defined + bases(qarr::AbstractQuArray, i) = rawbases(qarr, i) + bases(qarr::AbstractQuArray) = ntuple(ndims(qarr), i->bases(qarr, i)) + + ######################## + # Array-like functions # + ######################## + Base.size(qarr::AbstractQuArray, i...) = size(rawcoeffs(qarr), i...) + Base.ndims(qarr::AbstractQuArray) = ndims(rawcoeffs(qarr)) + Base.length(qarr::AbstractQuArray) = length(rawcoeffs(qarr)) + + Base.getindex(qarr::AbstractQuArray, i...) = getindex(rawcoeffs(qarr), i...) + Base.setindex!(qarr::AbstractQuArray, i...) = setindex!(rawcoeffs(qarr), i...) + + Base.in(c, qarr::AbstractQuArray) = in(c, rawcoeffs(qarr)) + + Base.(:(==))(a::AbstractQuArray, b::AbstractQuArray) = coeffs(a)==coeffs(b) && bases(a)==bases(b) + typealias AbstractQuVector{B<:AbstractBasis,T} AbstractQuArray{B,T,1} typealias AbstractQuMatrix{B<:AbstractBasis,T} AbstractQuArray{B,T,2} + +########### +# QuArray # +########### type QuArray{B<:AbstractBasis,T,N,A} <: AbstractQuArray{B,T,N} coeffs::A bases::NTuple{N,B} @@ -25,54 +56,44 @@ typealias QuVector{B<:AbstractBasis,T,A} QuArray{B,T,1,A} typealias QuMatrix{B<:AbstractBasis,T,A} QuArray{B,T,2,A} + # returns an appropriate outer constructor for the given (sub)type; + # the constructor should construct an instance from a coefficient container + # and a tuple of bases (see QuArray(coeffs, bases) above) + similar_type{Q<:QuArray}(::Type{Q}) = QuArray + + ###################### # Accessor functions # ###################### coefftype{B,T,N,A}(::QuArray{B,T,N,A}) = A rawcoeffs(qarr::QuArray) = qarr.coeffs - coeffs(qarr::QuArray) = rawcoeffs(qarr) + rawbases(qarr::QuArray) = qarr.bases rawbases(qarr::QuArray, i) = qarr.bases[i] - bases(qarr::QuArray, i) = rawbases(qarr, i) - rawbases(qarr::AbstractQuArray) = qarr.bases - - bases(qarr::QuArray, i) = rawbases(qarr, i) - # works generally as long as the single index form is defined - bases(qarr::AbstractQuArray) = ntuple(ndims(qarr), i->bases(qarr, i)) - - ######################## - # Array-like functions # - ######################## Base.copy(qa::QuArray) = QuArray(copy(qa.coeffs), copy(qa.bases)) - Base.size(qarr::QuArray, i...) = size(rawcoeffs(qarr), i...) - Base.ndims(qarr::QuArray) = ndims(rawcoeffs(qarr)) - Base.length(qarr::QuArray) = length(rawcoeffs(qarr)) - - Base.getindex(qarr::QuArray, i...) = getindex(rawcoeffs(qarr), i...) - Base.setindex!(qarr::QuArray, i...) = setindex!(rawcoeffs(qarr), i...) - - Base.in(c, qarr::QuArray) = in(c, rawcoeffs(qarr)) - - Base.(:(==))(a::AbstractQuArray, b::AbstractQuArray) = coeffs(a)==coeffs(b) && bases(a)==bases(b) ############## # CTranspose # ############## - immutable CTranspose{B,T,N,A} <: AbstractQuArray{B,T,N} + 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") end - CTranspose{B,T,N,A}(qa::QuArray{B,T,N,A}) = CTranspose{B,T,N,A}(qa) + CTranspose{B<:AbstractBasis,T,N,A}(qa::QuArray{B,T,N,A}) = CTranspose{B,T,N,A}(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} + similar_type{QC<:CTranspose}(::Type{QC}) = CTranspose + + ###################### # Accessor functions # ###################### @@ -88,6 +109,7 @@ bases(ct::CTranspose, i) = rawbases(ct, revind(ndims(ct), i)) + ######################## # Array-like functions # ######################## @@ -108,6 +130,7 @@ Base.ctranspose(qarr::QuArray) = CTranspose(qarr) Base.ctranspose(ct::CTranspose) = ct.qarr + ################ # LabelQuArray # ################ @@ -117,6 +140,13 @@ 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)...) + +################### +# Promotion rules # +################### + + 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)} + ###################### # Printing Functions # ###################### @@ -129,9 +159,11 @@ print(io, repr(coeffs(qarr))) end + #################### # Helper Functions # #################### + typerepr(qa::AbstractQuArray) = "$(typeof(qa).name)" typerepr(::QuVector) = "QuVector" typerepr(::QuMatrix) = "QuMatrix" typerepr(::DualVector) = "DualVector" @@ -148,7 +180,9 @@ return false end -export QuArray, +export QuArray, QuVector, QuMatrix, CTranspose, DualVector, DualMatrix, rawcoeffs, - coeffs - + coeffs, + coefftype, + rawbases, + bases diff --git a/src/arrays/spinops.jl b/src/arrays/spinops.jl index e5ce429..298ac50 100644 --- a/src/arrays/spinops.jl +++ b/src/arrays/spinops.jl @@ -12,9 +12,9 @@ spin(s::HalfSpin) = s spin(j::Integer) = HalfSpin(2*j) function spin(j) if isinteger(j) - return spin(int(j)) # j = 0.0, 1.0, 2.0... + return spin(@compat(Int(j))) # j = 0.0, 1.0, 2.0... elseif isinteger(2*j) - return HalfSpin(int(2*j)) # j = .5, 1.5, 2.5... + return HalfSpin(@compat(Int(2*j))) # j = .5, 1.5, 2.5... else error("Spin must be a multiple of 1/2.") end diff --git a/src/bases/labelbasis.jl b/src/bases/labelbasis.jl index c5ce295..30722ac 100644 --- a/src/bases/labelbasis.jl +++ b/src/bases/labelbasis.jl @@ -146,10 +146,10 @@ ###################### # Accessor Functions # ###################### - ind_value(n, range, denom, modulus) = range[(div(n, denom) % modulus)+1] + ind_value(n, range, denom, modulus) = range[@compat(Int(div(n, denom) % modulus))+1] tuple_at_ind(b::LabelBasis, i) = ntuple(nfactors(b), x->ind_value(i-1, ranges(b,x), b.denoms[x], size(b,x))) pos_in_range(r::Range, i) = i in r ? (i-first(r))/step(r) : throw(BoundsError()) - getpos(b::LabelBasis, label) = int(sum(map(*, map(pos_in_range, ranges(b), label), b.denoms))+1) + getpos(b::LabelBasis, label) = @compat Int(sum(map(*, map(pos_in_range, ranges(b), label), b.denoms))+1) Base.in(label, b::LabelBasis) = reduce(&, map(in, label, ranges(b)))