Merge pull request #29 from acroy/abstractqa

Consistently use AbstractQuArray and define its interface in terms of functions rawcoeffs, rawbases, coefftype, similar_type and copy. Fix deprecation warnings for 0.4-dev.

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,79 +45,111 @@ 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 #
 ##################
 
 # 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})

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)

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,