Start operator type overhaul

src/QuantumOptics.jl

@@ -8,12 +8,12 @@ export bases, Basis, GenericBasis, CompositeBasis, basis,
                 dagger, normalize, normalize!,
         operators, AbstractOperator, expect, variance, identityoperator, ptrace, embed, dense, tr,
             sparse,
-        operators_dense, DenseOperator, projector, dm,
-        operators_sparse, SparseOperator, diagonaloperator,
+        operators_dense, projector, dm, mul!,
+        operators_sparse, diagonaloperator,
         operators_lazysum, LazySum,
         operators_lazyproduct, LazyProduct,
         operators_lazytensor, LazyTensor,
-        superoperators, SuperOperator, DenseSuperOperator, SparseSuperOperator,
+        superoperators, AbstractSuperOperator, SuperOperator,
                 spre, spost, liouvillian,
         fock, FockBasis, number, destroy, create,
                 fockstate, coherentstate, displace,

src/manybody.jl

@@ -20,14 +20,14 @@ The basis has to know the associated one-body basis `b` and which occupation sta
 should be included. The occupations_hash is used to speed up checking if two
 many-body bases are equal.
 """
-mutable struct ManyBodyBasis <: Basis
+mutable struct ManyBodyBasis{B<:Basis} <: Basis
     shape::Vector{Int}
-    onebodybasis::Basis
+    onebodybasis::B
     occupations::Vector{Vector{Int}}
     occupations_hash::UInt
 
     function ManyBodyBasis(onebodybasis::Basis, occupations::Vector{Vector{Int}})
-        new([length(occupations)], onebodybasis, occupations, hash(hash.(occupations)))
+        new{B<:Basis}([length(occupations)], onebodybasis, occupations, hash(hash.(occupations)))
     end
 end
 
@@ -91,8 +91,8 @@ end
 
 Creation operator for the i-th mode of the many-body basis `b`.
 """
-function create(b::ManyBodyBasis, index::Int)
-    result = SparseOperator(b)
+function create(b::B, index::Int) where B<:ManyBodyBasis
+    result = Operator{B,B,SparseMatrixCSC{ComplexF64,Int}}(b)
     # <{m}_i| at |{m}_j>
     for i=1:length(b)
         occ_i = b.occupations[i]
@@ -113,8 +113,8 @@ end
 
 Annihilation operator for the i-th mode of the many-body basis `b`.
 """
-function destroy(b::ManyBodyBasis, index::Int)
-    result = SparseOperator(b)
+function destroy(b::B, index::Int) where B<:ManyBodyBasis
+    result = Operator{B,B,SparseMatrixCSC{ComplexF64,Int}}(b)
     # <{m}_j| a |{m}_i>
     for i=1:length(b)
         occ_i = b.occupations[i]
@@ -136,7 +136,7 @@ end
 Particle number operator for the i-th mode of the many-body basis `b`.
 """
 function number(b::ManyBodyBasis, index::Int)
-    result = SparseOperator(b)
+    result = sparse(Operator(b))
     for i=1:length(b)
         result.data[i, i] = b.occupations[i][index]
     end
@@ -149,7 +149,7 @@ end
 Total particle number operator.
 """
 function number(b::ManyBodyBasis)
-    result = SparseOperator(b)
+    result = sparse(Operator(b))
     for i=1:length(b)
         result.data[i, i] = sum(b.occupations[i])
     end
@@ -185,7 +185,7 @@ end
 Operator ``|\\mathrm{to}⟩⟨\\mathrm{from}|`` transferring particles between modes.
 """
 function transition(b::ManyBodyBasis, to::Int, from::Int)
-    result = SparseOperator(b)
+    result = sparse(Operator(b))
     # <{m}_j| at_to a_from |{m}_i>
     for i=1:length(b)
         occ_i = b.occupations[i]
@@ -228,22 +228,22 @@ where ``X`` is the N-particle operator, ``x`` is the one-body operator and
 ``|u⟩`` are the one-body states associated to the
 different modes of the N-particle basis.
 """
-function manybodyoperator(basis::ManyBodyBasis, op::T)::T where T<:AbstractOperator
-    @assert op.basis_l == op.basis_r
-    if op.basis_l == basis.onebodybasis
-        result =  manybodyoperator_1(basis, op)
-    elseif op.basis_l == basis.onebodybasis ⊗ basis.onebodybasis
-        result = manybodyoperator_2(basis, op)
-    else
-        throw(ArgumentError("The basis of the given operator has to either be equal to b or b ⊗ b where b is the 1st quantization basis associated to the nparticle basis."))
-    end
-    result
-end
-
-function manybodyoperator_1(basis::ManyBodyBasis, op::DenseOperator)
+# function manybodyoperator(basis::ManyBodyBasis, op::T) where {T<:AbstractOperator{B,B},B<:Basis}
+#     # @assert op.basis_l == op.basis_r
+#     if op.basis_l == basis.onebodybasis
+#         result =  manybodyoperator_1(basis, op)
+#     elseif op.basis_l == basis.onebodybasis ⊗ basis.onebodybasis
+#         result = manybodyoperator_2(basis, op)
+#     else
+#         throw(ArgumentError("The basis of the given operator has to either be equal to b or b ⊗ b where b is the 1st quantization basis associated to the nparticle basis."))
+#     end
+#     result
+# end
+
+function manybodyoperator(basis::ManyBodyBasis{B}, op::Operator{B,B,T}) where {B<:Basis,T<:Matrix{ComplexF64}}
     N = length(basis)
     S = length(basis.onebodybasis)
-    result = DenseOperator(basis)
+    result = Operator(basis)
     @inbounds for n=1:N, m=1:N
         for j=1:S, i=1:S
             C = coefficient(basis.occupations[m], basis.occupations[n], [i], [j])
@@ -255,10 +255,10 @@ function manybodyoperator_1(basis::ManyBodyBasis, op::DenseOperator)
     return result
 end
 
-function manybodyoperator_1(basis::ManyBodyBasis, op::SparseOperator)
+function manybodyoperator(basis::ManyBodyBasis{B}, op::Operator{B,B,T}) where {B<:Basis,T<:SparseMatrixCSC{ComplexF64,Int}}
     N = length(basis)
     S = length(basis.onebodybasis)
-    result = SparseOperator(basis)
+    result = sparse(Operator(basis))
     M = op.data
     @inbounds for colindex = 1:M.n
         for i=M.colptr[colindex]:M.colptr[colindex+1]-1
@@ -275,12 +275,13 @@ function manybodyoperator_1(basis::ManyBodyBasis, op::SparseOperator)
     return result
 end
 
-function manybodyoperator_2(basis::ManyBodyBasis, op::DenseOperator)
+function manybodyoperator(basis::ManyBodyBasis{BC}, op::Operator{B,B,T}) where {B<:Basis,BC<:CompositeBasis{Vector{Basis}},T<:Matrix{ComplexF64}}
+    @assert basis.onebodybasis == op.basis_l^2
     N = length(basis)
     S = length(basis.onebodybasis)
     @assert S^2 == length(op.basis_l)
     @assert S^2 == length(op.basis_r)
-    result = DenseOperator(basis)
+    result = Operator(basis)
     op_data = reshape(op.data, S, S, S, S)
     occupations = basis.occupations
     @inbounds for m=1:N, n=1:N
@@ -292,10 +293,11 @@ function manybodyoperator_2(basis::ManyBodyBasis, op::DenseOperator)
     return result
 end
 
-function manybodyoperator_2(basis::ManyBodyBasis, op::SparseOperator)
+function manybodyoperator(basis::ManyBodyBasis{BC}, op::Operator{B,B,T}) where {BC<:CompositeBasis{Vector{Basis}},B<:Basis,T<:SparseMatrixCSC{ComplexF64,Int}}
+    @assert basis.onebodybasis == op.basis_l^2
     N = length(basis)
     S = length(basis.onebodybasis)
-    result = SparseOperator(basis)
+    result = sparse(Operator(basis))
     occupations = basis.occupations
     rows = rowvals(op.data)
     values = nonzeros(op.data)
@@ -321,9 +323,7 @@ end
 
 Expectation value of the one-body operator `op` in respect to the many-body `state`.
 """
-function onebodyexpect(op::AbstractOperator, state::Ket)
-    @assert isa(state.basis, ManyBodyBasis)
-    @assert op.basis_l == op.basis_r
+function onebodyexpect(op::AbstractOperator{B,B}, state::Ket{MB}) where {B<:Basis,MB<:ManyBodyBasis{B}}
     if state.basis.onebodybasis == op.basis_l
         result = onebodyexpect_1(op, state)
     # Not yet implemented:
@@ -335,10 +335,7 @@ function onebodyexpect(op::AbstractOperator, state::Ket)
     result
 end
 
-function onebodyexpect(op::AbstractOperator, state::AbstractOperator)
-    @assert op.basis_l == op.basis_r
-    @assert state.basis_l == state.basis_r
-    @assert isa(state.basis_l, ManyBodyBasis)
+function onebodyexpect(op::AbstractOperator{B,B}, state::AbstractOperator{BM,BM}) where {B<:Basis,BM<:ManyBodyBasis{B}}
     if state.basis_l.onebodybasis == op.basis_l
         result = onebodyexpect_1(op, state)
     # Not yet implemented
@@ -351,7 +348,7 @@ function onebodyexpect(op::AbstractOperator, state::AbstractOperator)
 end
 onebodyexpect(op::AbstractOperator, states::Vector) = [onebodyexpect(op, state) for state=states]
 
-function onebodyexpect_1(op::DenseOperator, state::Ket)
+function onebodyexpect_1(op::Operator{B,B,T}, state::Ket{BM}) where {B<:Basis,T<:Matrix{ComplexF64},BM<:ManyBodyBasis{B}}
     N = length(state.basis)
     S = length(state.basis.onebodybasis)
     result = complex(0.)
@@ -368,7 +365,7 @@ function onebodyexpect_1(op::DenseOperator, state::Ket)
     result
 end
 
-function onebodyexpect_1(op::DenseOperator, state::DenseOperator)
+function onebodyexpect_1(op::Operator{B,B,T}, state::Operator{BM,BM,T}) where {B<:Basis,T<:Matrix{ComplexF64},BM<:ManyBodyBasis{B}}
     N = length(state.basis_l)
     S = length(state.basis_l.onebodybasis)
     result = complex(0.)
@@ -385,7 +382,7 @@ function onebodyexpect_1(op::DenseOperator, state::DenseOperator)
     result
 end
 
-function onebodyexpect_1(op::SparseOperator, state::Ket)
+function onebodyexpect_1(op::Operator{B,B,T}, state::Ket{BM}) where {B<:Basis,T<:SparseMatrixCSC{ComplexF64,Int},BM<:ManyBodyBasis{B}}
     N = length(state.basis)
     S = length(state.basis.onebodybasis)
     result = complex(0.)
@@ -406,7 +403,7 @@ function onebodyexpect_1(op::SparseOperator, state::Ket)
     result
 end
 
-function onebodyexpect_1(op::SparseOperator, state::DenseOperator)
+function onebodyexpect_1(op::Operator{B,B,T}, state::Operator{BM,BM,T2}) where {B<:Basis,T<:SparseMatrixCSC{ComplexF64,Int},BM<:ManyBodyBasis{B},T2<:Matrix{ComplexF64}}
     N = length(state.basis_l)
     S = length(state.basis_l.onebodybasis)
     result = complex(0.)

src/operators.jl

@@ -24,7 +24,7 @@ For fast time evolution also at least the function
 implemented. Many other generic multiplication functions can be defined in
 terms of this function and are provided automatically.
 """
-abstract type AbstractOperator end
+abstract type AbstractOperator{BL<:Basis,BR<:Basis} end
 
 
 # Common error messages