Basic framework for QuArrays and basis types

.travis.yml

@@ -0,0 +1,4 @@
+language: julia 
+  julia: 
+      - release 
+      - nightly 

REQUIRE

@@ -0,0 +1 @@
+julia 0.3-

src/QuBase.jl

@@ -1,51 +1,65 @@
 module QuBase
-
-	#############
-	# Constants #
-	#############
-		const lang = "\u27E8"
-		const rang = "\u27E9"
-		const otimes = "\u2297"
-		const vdots ="\u205E"
+
+    import Base: kron
+
+    ####################
+    # String Constants #
+    ####################
+        const lang = "\u27E8"
+        const rang = "\u27E9"
+        const otimes = "\u2297"
+        const vdots ="\u205E"
 
-	##################
-	# Abstract Types #
-	##################
-		abstract AbstractStructure
-		abstract AbstractQuantum{S<:AbstractStructure}
-		abstract AbstractState{S<:AbstractStructure} <: AbstractQuantum{S}
-		abstract AbstractOperator{S<:AbstractStructure} <: AbstractQuantum{S}
-		abstract AbstractBasis{S<:AbstractStructure} <: AbstractQuantum{S}
-		abstract QuantumScalar <: Number
+    ##################
+    # Abstract Types #
+    ##################
+        abstract AbstractStructure
+        abstract Orthonormal <: AbstractStructure
+        abstract AbstractQuantum{S<:AbstractStructure}
 
-	#############
-	# Functions #
-	#############
-		structure{S}(::AbstractQuantum{S}) = S
+        # Various constructor methods in this repo allow an argument 
+        # of type Type{BypassFlag} to be passed in in order to 
+        # circumvent value precalculation/checking. This is useful for
+        # conversion methods and the like, where you know the input 
+        # has already been vetted elsewhere. Don't use this unless
+        # you're sure of what you're doing, and don't export this.
+        abstract BypassFlag
 
-		for T=(:AbstractQuantum, :AbstractState, :AbstractOperator, :AbstractBasis)
-			@eval begin
-				structure{S}(::Type{($T){S}}) = S
-				structure(::Type{($T)}) = AbstractStructure
-			end
-		end
+    #############
+    # Functions #
+    #############
 
-	######################
-	# Include Statements #
-	######################
-		include("statelabel.jl")
-		include("diracstates.jl")
-		include("diracoperators.jl")
-		include("scalar.jl")
-		include("quarray.jl")
+        # This should be the only `structure` 
+        # method that needs to be defined for 
+        # type *instances*
+        structure{S}(::AbstractQuantum{S}) = S
+        structure(::DataType) = AbstractStructure
 
-	export AbstractStructure, 
-		AbstractQuantum,
-		AbstractState,
-		AbstractOperator,
-		AbstractBasis,
-		QuantumScalar,
-		structure
-
+        # ...and all relevant singleton types 
+        # should have it defined as well:
+        structure{S}(::Type{AbstractQuantum{S}}) = S
+
+        # an n-arity form of the tensor
+        # product, reduction is done via
+        # the binary definition of tensor()
+        # defined in the files included above. 
+        tensor(s...) = reduce(tensor, s)
+
+        # For the sake of convenience, kron() 
+        # is defined to be equivalent to 
+        # tensor() for quantum objects
+        kron(a::AbstractQuantum, b::AbstractQuantum) = tensor(a, b)
+
+    ######################
+    # Include Statements #
+    ######################
+        include("bases/bases.jl")
+        include("arrays/quarray.jl")
+
+    export AbstractStructure, 
+        AbstractQuantum,
+        Orthonormal,
+        structure,
+        tensor 
 end
-	
+    

src/arrays/ladderops.jl

@@ -0,0 +1,68 @@
+###########################
+# Ladder Operator Methods #
+########################### 
+    # n specifies which particle (a.k.a tensor product 
+    # factor) the operator acts on
+    raisematrix(lens, n=1) = laddermatrix(lens, RaiseOpFlag, n)
+    lowermatrix(lens, n=1) = laddermatrix(lens, LowerOpFlag, n)
+
+    raiseop(b::AbstractBasis, n=1) = QuArray(raisematrix(size(b), n), b, b)
+    raiseop(lens::Tuple, n=1) = raiseop(FiniteBasis(lens), n)
+
+    lowerop(b::AbstractBasis, n=1) = QuArray(lowermatrix(size(b), n), b, b)
+    lowerop(lens::Tuple, n=1) = lowerop(FiniteBasis(lens), n)
+
+    ##########################
+    # Helper Functions/Types #
+    ##########################
+    abstract LadderOpFlag
+    abstract RaiseOpFlag <: LadderOpFlag
+    abstract LowerOpFlag <: LadderOpFlag    
+
+    ladder_inds(n, ::Type{RaiseOpFlag}) = ([1:n-1], [2:n])
+    ladder_inds(n, ::Type{LowerOpFlag}) = ([2:n], [1:n-1])
+    laddercoeffs(n) = sqrt(linspace(1, n, n))
+
+    function fill_op_arr!(arr::AbstractMatrix, ladderflag)
+        if size(arr, 1) == size(arr, 2)
+            len = size(arr, 1)
+            inds = ladder_inds(len, ladderflag)
+            coeffs = laddercoeffs(len)
+            for i=1:len-1
+                arr[inds[1][i], inds[2][i]] = coeffs[i] 
+            end
+        else
+            error("Cannot generate ladder coefficients for non-square matrix")
+        end
+        return arr
+    end
+
+    # returns a coefficient matrix
+    # for a ladder operator for a
+    # single particle fock basis
+    gen_op_mat(len, ladderflag) = fill_op_arr!(spzeros(len, len), ladderflag)
+
+    # this could/should be further optimized,
+    # it uses the naive approach of taking the 
+    # kronecker product of identity matrices 
+    # and the relevant ladder operator matrix
+    function laddermatrix(lens, ladderflag, n=1)   
+        if n==1
+            arr = gen_op_mat(lens[1], ladderflag)
+        else
+            arr = speye(lens[1])
+            for i=2:n-1
+                arr = kron(speye(lens[i]), arr)
+            end
+            arr = kron(gen_op_mat(lens[n], ladderflag), arr)
+        end 
+        for i=n+1:length(lens)
+            arr = kron(speye(lens[i]), arr)
+        end
+        return arr
+    end
+
+export raisematrix,
+    lowermatrix,
+    raiseop,
+    lowerop

src/arrays/quarray.jl

@@ -0,0 +1,141 @@
+import Base: size,
+    length,
+    getindex,
+    similar,
+    in,
+    ctranspose,
+    transpose,
+    summary,
+    zeros,
+    eye,
+    #TODO: Implement the below operations
+    *,.*,
+    /,./,
+    +,.+,
+    -,.-,
+    kron
+
+abstract AbstractQuArray{B<:AbstractBasis,T,N} <: AbstractArray{T,N}
+
+###########
+# QuArray #                 
+###########
+    # Using an NTuple allows us to have a basis for each dimension, 
+    # and gives us a less ambiguous way to determine if a QuArray will act
+    # like a vector, matrix, or tensor using the dimension parameter N. 
+    type QuArray{B<:AbstractBasis,T,N,A} <: AbstractQuArray{B,T,N}
+        coeffs::A
+        bases::NTuple{N,B}
+        function QuArray(coeffs::AbstractArray{T}, bases::NTuple{N,B}) 
+            if checkbases(coeffs, bases) 
+                new(coeffs, bases)
+            else 
+                error("Coefficient array does not conform to input bases")
+            end
+        end
+    end
+
+    typealias QuVector{B<:AbstractBasis,T,A} QuArray{B,T,1,A}
+    typealias QuMatrix{B<:AbstractBasis,T,A} QuArray{B,T,2,A}
+
+    QuArray{T,N,B<:AbstractBasis}(coeffs::AbstractArray{T}, bases::NTuple{N,B}) = QuArray{B,T,N,typeof(coeffs)}(coeffs, bases)
+    QuArray(coeffs, bases::AbstractBasis...) = QuArray(coeffs, bases)
+    QuArray(coeffs) = QuArray(coeffs, basesfordims(size(coeffs)))
+
+    ############################
+    # Convenience Constructors #
+    ############################
+    statevec(i::Int, fb::FiniteBasis) = QuArray(single_coeff(i, length(fb)), fb)
+    statevec(i::Int, lens::Int...=i) = statevec(i, FiniteBasis(lens))
+
+    zeros(qa::QuArray) = QuArray(zeros(qa.coeffs), qa.bases)
+    eye(qa::QuArray) = QuArray(eye(qa.coeffs), qa.bases)
+
+    ######################
+    # Property Functions #
+    ######################
+    bases(qa::QuArray) = qa.bases
+    coeffs(qa::QuArray) = qa.coeffs
+    size(qa::QuArray, i...) = size(qa.coeffs, i...)
+
+    ########################
+    # Array-like functions #
+    ########################
+    similar{B,T}(qa::QuArray{B,T}, element_type=T) = QuArray(similar(qa.coeffs, T), qa.bases)
+    # Is there a way to properly define the below for
+    # any arbitrary basis? Obviously doesn't make sense
+    # for B<:AbstractInfiniteBasis, and I'm reluctant to
+    # enforce that every B<:AbstractFiniteBasis will have a 
+    # constructor B(::Int), which is how the below is constructing
+    # instances of FiniteBasis.
+    function similar{B<:FiniteBasis,T}(qa::QuArray{B,T}, element_type, dims)
+        return QuArray(similar(qa.coeffs, T, dims), basesfordims(dims, B))
+    end 
+    getindex(qa::QuArray, i::AbstractArray) = getindex(qa.coeffs, i)
+    getindex(qa::QuArray, i::Real) = getindex(qa.coeffs, i)
+    getindex(qa::QuArray, i) = getindex(qa.coeffs, i)
+    getindex(qa::QuArray, i...) = getindex(qa.coeffs, i...)
+
+    in(c, qa::QuArray) = in(c, qa.coeffs)
+
+    ctranspose(qa::QuVector) = QuArray(ctranspose(qa.coeffs), qa.bases)
+    ctranspose(qa::QuMatrix) = QuArray(ctranspose(qa.coeffs), reverse(qa.bases))
+    transpose(qa::QuVector) = QuArray(transpose(qa.coeffs), qa.bases)
+    transpose(qa::QuMatrix) = QuArray(transpose(qa.coeffs), reverse(qa.bases))
+
+    ######################
+    # Printing Functions #
+    ######################
+    function summary{B,T,N,A}(qa::QuArray{B,T,N,A})
+        return "$(sizenotation(size(qa))) QuArray\n" *
+               "...bases: $B,\n" * 
+               "...coeff: $A"             
+    end
+
+    ######################
+    # Include Statements #
+    ######################
+    include("ladderops.jl")
+
+    ####################
+    # Helper Functions #
+    ####################
+    sizenotation(tup::(Int,)) = "$(first(tup))-element"
+    sizenotation(tup::(Int...)) = reduce(*, map(s->"$(s)x", tup))[1:end-1] 
+
+    # checkbases() is overloaded for a single basis
+    # to handle the fact that row vectors are
+    # N=2
+    function checkbases(coeffs, bases::NTuple{1, AbstractBasis})
+        if ndims(coeffs) <= 2
+            if size(coeffs, 2) == 1
+                return checkcoeffs(coeffs, 1, first(bases))
+            elseif size(coeffs, 1) == 1
+                return checkcoeffs(coeffs, 2, first(bases))
+            end 
+        end
+        return false
+    end
+
+    function checkbases{N}(coeffs, bases::NTuple{N, AbstractBasis})
+        if ndims(coeffs) == length(bases)
+            return reduce(&, [checkcoeffs(coeffs, i, bases[i]) for i=1:N])
+        end
+        return false
+    end
+
+    # Assumes that every basis type passed in
+    # has a constructor B(::eltype(lens))
+    function basesfordims(lens::Tuple, B=ntuple(length(lens), x->FiniteBasis))
+        return ntuple(length(lens), n->B[n](lens[n]))
+    end
+    one_at_ind!(arr, i) = setindex!(arr, one(eltype(arr)), i)
+    single_coeff(i, lens...) = one_at_ind!(zeros(lens), i)
+
+export AbstractQuArray,
+    QuArray,
+    QuVector,
+    QuMatrix,
+    bases,
+    coeffs,
+    statevec

src/bases/bases.jl

@@ -0,0 +1,55 @@
+####################
+# Type Definitions #
+####################
+    abstract AbstractBasis{S<:AbstractStructure} <: AbstractQuantum{S}
+    abstract AbstractFiniteBasis{S<:AbstractStructure} <: AbstractBasis{S}
+    abstract AbstractInfiniteBasis{S<:AbstractStructure} <: AbstractBasis{S}
+
+#############
+# Functions #
+#############
+    # Note: All B<:AbstractBasis types should implement the following: 
+    #
+    # checkcoeffs(coeffs, dim, basis::B) -> checks whether a coefficient
+    #                                       array is valid for use with 
+    #                                       the given basis. This is used 
+    #                                       by QuArray to ensure that the 
+    #                                       coefficient array is not malformed 
+    #                                       with respect to the input bases. 
+    #                                       The second argument, `dim`, specifies 
+    #                                       the dimension of the coefficient 
+    #                                       array which corresponds to the provided 
+    #                                       basis.
+    #
+    # nfactors(basis::B) -> the number of factor bases for `basis`; this is `N`
+    #                       in `FiniteBasis{S,N}` 
+    #
+    # tensor(a::B, b::B) -> Take the tensor product of these two bases. This
+    #                       function should optimally return a basis of same 
+    #                       type as the input bases. We can and should also 
+    #                       implement promote_rules/conversion methods between
+    #                       bases.
+    #
+    # structure{S}(::Type{MyBasisType{S}}) -> returns S<:AbstractStructure for 
+    #                                         the provided basis type.
+
+    checkcoeffs(coeffs::AbstractArray, dim::Int, basis::AbstractBasis) = error("checkcoeffs(coeffs, dim, ::$B) must be defined!")
+
+    for basis=(:AbstractBasis, :AbstractFiniteBasis, :AbstractInfiniteBasis)
+        @eval begin
+            structure{S}(::Type{($basis){S}}) = S
+        end
+    end
+
+######################
+# Include Statements #
+######################
+    include("finitebasis.jl")
+
+export AbstractBasis,
+    AbstractFiniteBasis,
+    AbstractInfiniteBasis,
+    checkcoeffs,
+    structure
+
+