Loosen type constraints

Project.toml

@@ -25,7 +25,7 @@ FFTW = "1"
 IterativeSolvers = "0.9"
 LinearMaps = "3"
 OrdinaryDiffEq = "5"
-QuantumOpticsBase = "0.2.7"
+QuantumOpticsBase = "0.3"
 RecursiveArrayTools = "2"
 Reexport = "0.2, 1.0"
 StochasticDiffEq = "6"

src/bloch_redfield_master.jl

@@ -15,7 +15,7 @@ See QuTiP's documentation (http://qutip.org/docs/latest/guide/dynamics/dynamics-
 * `secular_cutoff=0.1`: Cutoff to allow a degree of partial secularization. Terms are discarded if they are greater than (dw\\_min * secular cutoff) where dw\\_min is the smallest (non-zero) difference between any two eigenenergies of H.
                         This argument is only taken into account if use_secular=true.
 """
-function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secular=true, secular_cutoff=0.1)
+function bloch_redfield_tensor(H::AbstractOperator, a_ops; J=SparseOpType[], use_secular=true, secular_cutoff=0.1)
 
     # Use the energy eigenbasis
     H_evals, transf_mat = eigen(Array(H.data)) #Array call makes sure H is a dense array
@@ -45,7 +45,7 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
     end
 
     #Transform interaction operators to Hamiltonian eigenbasis
-    A = Array{ComplexF64}(undef, N, N, K)
+    A = Array{eltype(transf_mat)}(undef, N, N, K)
     for k in 1:K
         A[:, :, k] = to_Heb(a_ops[k][1], transf_mat).data
     end
@@ -54,8 +54,8 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
     W = H_evals .- transpose(H_evals)
 
     #Array for spectral functions evaluated at transition frequencies
-    Jw = Array{Float64}(undef, N, N, K)
-    #Loop over all a_ops and calculate each spectral density at all transition frequencies
+    Jw = Array{eltype(transf_mat)}(undef, N, N, K)
+    # Jw = zeros(Complex{Float64}, N, N, K)
     for k in 1:K
         Jw[:, :, k] .= a_ops[k][2].(W)
     end
@@ -106,7 +106,6 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
 end #Function
 
 
-
 """
     timeevolution.master_bloch_redfield(tspan, rho0, R, H; <keyword arguments>)
 
@@ -127,9 +126,9 @@ Time-evolution according to a Bloch-Redfield master equation.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
 function master_bloch_redfield(tspan,
-        rho0::T, L::SuperOperator{Tuple{B,B},Tuple{B,B}},
+        rho0::Operator{B,B}, L::SuperOperator{Tuple{B,B},Tuple{B,B}},
         H::AbstractOperator{B,B}; fout::Union{Function,Nothing}=nothing,
-        kwargs...) where {B<:Basis,T<:Operator{B,B}}
+        kwargs...) where {B}
 
     #Prep basis transf
     evals, transf_mat = eigen(dense(H).data)
@@ -142,28 +141,28 @@ function master_bloch_redfield(tspan,
     L_ = isa(L, SparseSuperOpType) ? SparseOperator(basis_comp, L.data) : DenseOperator(basis_comp, L.data)
 
     # Derivative function
-    dmaster_br_(t, rho::T2, drho::T2) where T2<:Ket = dmaster_br(drho, rho, L_)
+    dmaster_br_(t, rho, drho) = dmaster_br(drho, rho, L_)
 
     return integrate_br(tspan, dmaster_br_, rho0_eb, transf_op, inv_transf_op, fout; kwargs...)
 end
 master_bloch_redfield(tspan, psi::Ket, args...; kwargs...) = master_bloch_redfield(tspan, dm(psi), args...; kwargs...)
 
 # Derivative ∂ₜρ = Lρ
-function dmaster_br(drho::T, rho::T, L::DataOperator{B,B}) where {B<:Basis,T<:Ket{B}}
+function dmaster_br(drho, rho, L)
     QuantumOpticsBase.mul!(drho,L,rho)
 end
 
 # Integrate if there is no fout specified
-function integrate_br(tspan, dmaster_br::Function, rho::T,
-                transf_op::T2, inv_transf_op::T2, ::Nothing;
-                kwargs...) where {T<:Ket,T2<:Operator}
+function integrate_br(tspan, dmaster_br, rho,
+                transf_op, inv_transf_op, ::Nothing;
+                kwargs...)
     # Pre-allocate for in-place back-transformation from eigenbasis
     rho_out = copy(transf_op)
     tmp = copy(transf_op)
     tmp2 = copy(transf_op)
 
     # Define fout
-    function fout(t, rho::T)
+    function fout(t, rho)
         tmp.data[:] = rho.data
         QuantumOpticsBase.mul!(tmp2,transf_op,tmp)
         QuantumOpticsBase.mul!(rho_out,tmp2,inv_transf_op)
@@ -174,9 +173,9 @@ function integrate_br(tspan, dmaster_br::Function, rho::T,
 end
 
 # Integrate with given fout
-function integrate_br(tspan, dmaster_br::Function, rho::T,
-                transf_op::T2, inv_transf_op::T2, fout::Function;
-                kwargs...) where {T<:Ket,T2<:Operator}
+function integrate_br(tspan, dmaster_br, rho,
+                transf_op, inv_transf_op, fout::Function;
+                kwargs...)
     # Pre-allocate for in-place back-transformation from eigenbasis
     rho_out = copy(transf_op)
     tmp = copy(transf_op)
@@ -185,7 +184,7 @@ function integrate_br(tspan, dmaster_br::Function, rho::T,
     tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
 
     # Perform back-transfomration before calling fout
-    function fout_(t, rho::T)
+    function fout_(t, rho)
         tmp.data[:] = rho.data
         QuantumOpticsBase.mul!(tmp2,transf_op,tmp)
         QuantumOpticsBase.mul!(rho_out,tmp2,inv_transf_op)

src/master.jl

@@ -1,20 +1,18 @@
-const DecayRates = Union{Vector, Matrix, Nothing}
-
 """
     timeevolution.master_h(tspan, rho0, H, J; <keyword arguments>)
 
 Integrate the master equation with dmaster_h as derivative function.
 
 Further information can be found at [`master`](@ref).
 """
-function master_h(tspan, rho0::T, H::AbstractOperator{B,B}, J::Vector;
-                rates::DecayRates=nothing,
-                Jdagger::Vector=dagger.(J),
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Operator{B,B}}
+function master_h(tspan, rho0::Operator, H::AbstractOperator, J;
+                rates=nothing,
+                Jdagger=dagger.(J),
+                fout=nothing,
+                kwargs...)
     check_master(rho0, H, J, Jdagger, rates)
     tmp = copy(rho0)
-    dmaster_(t, rho::T, drho::T) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
+    dmaster_(t, rho, drho) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
 end
 
@@ -29,15 +27,15 @@ H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k
 ```
 Further information can be found at [`master`](@ref).
 """
-function master_nh(tspan, rho0::T, Hnh::AbstractOperator{B,B}, J::Vector;
-                rates::DecayRates=nothing,
+function master_nh(tspan, rho0::Operator, Hnh::AbstractOperator, J;
+                rates=nothing,
                 Hnhdagger::AbstractOperator=dagger(Hnh),
-                Jdagger::Vector=dagger.(J),
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Operator{B,B}}
+                Jdagger=dagger.(J),
+                fout=nothing,
+                kwargs...)
     check_master(rho0, Hnh, J, Jdagger, rates)
     tmp = copy(rho0)
-    dmaster_(t, rho::T, drho::T) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
+    dmaster_(t, rho, drho) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
 end
 
@@ -73,21 +71,21 @@ non-hermitian Hamiltonian and then calls master_nh which is slightly faster.
         be changed.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
-function master(tspan, rho0::T, H::AbstractOperator{B,B}, J::Vector;
-                rates::DecayRates=nothing,
-                Jdagger::Vector=dagger.(J),
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Operator{B,B}}
+function master(tspan, rho0::Operator, H::AbstractOperator, J;
+                rates=nothing,
+                Jdagger=dagger.(J),
+                fout=nothing,
+                kwargs...)
     isreducible = check_master(rho0, H, J, Jdagger, rates)
     if !isreducible
         tmp = copy(rho0)
-        dmaster_h_(t, rho::T, drho::T) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
+        dmaster_h_(t, rho, drho) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
         return integrate_master(tspan, dmaster_h_, rho0, fout; kwargs...)
     else
         Hnh = nh_hamiltonian(H,J,Jdagger,rates)
         Hnhdagger = dagger(Hnh)
         tmp = copy(rho0)
-        dmaster_nh_(t, rho::T, drho::T) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
+        dmaster_nh_(t, rho, drho) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
         return integrate_master(tspan, dmaster_nh_, rho0, fout; kwargs...)
     end
 end
@@ -132,12 +130,12 @@ The given function can either be of the form `f(t, rho) -> (Hnh, Hnhdagger, J, J
 or `f(t, rho) -> (Hnh, Hnhdagger, J, Jdagger, rates)` For further information look
 at [`master_dynamic`](@ref).
 """
-function master_nh_dynamic(tspan, rho0::T, f::Function;
-                rates::DecayRates=nothing,
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Operator{B,B}}
+function master_nh_dynamic(tspan, rho0::Operator, f;
+                rates=nothing,
+                fout=nothing,
+                kwargs...)
     tmp = copy(rho0)
-    dmaster_(t, rho::T, drho::T) = dmaster_nh_dynamic(t, rho, f, rates, drho, tmp)
+    dmaster_(t, rho, drho) = dmaster_nh_dynamic(t, rho, f, rates, drho, tmp)
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
 end
 
@@ -166,22 +164,20 @@ operators:
         be changed.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
-function master_dynamic(tspan, rho0::T, f::Function;
-                rates::DecayRates=nothing,
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Operator{B,B}}
+function master_dynamic(tspan, rho0::Operator, f;
+                rates=nothing,
+                fout=nothing,
+                kwargs...)
     tmp = copy(rho0)
-    dmaster_(t, rho::T, drho::T) = dmaster_h_dynamic(t, rho, f, rates, drho, tmp)
+    dmaster_(t, rho, drho) = dmaster_h_dynamic(t, rho, f, rates, drho, tmp)
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
 end
 
 
 # Automatically convert Ket states to density operators
-master(tspan, psi0::Ket{B}, args...; kwargs...) where B<:Basis = master(tspan, dm(psi0), args...; kwargs...)
-master_h(tspan, psi0::Ket{B}, H::AbstractOperator{B,B}, J::Vector; kwargs...) where B<:Basis = master_h(tspan, dm(psi0), H, J; kwargs...)
-master_nh(tspan, psi0::Ket{B}, Hnh::AbstractOperator{B,B}, J::Vector; kwargs...) where B<:Basis = master_nh(tspan, dm(psi0), Hnh, J; kwargs...)
-master_dynamic(tspan, psi0::Ket{B}, f::Function; kwargs...) where B<:Basis = master_dynamic(tspan, dm(psi0), f; kwargs...)
-master_nh_dynamic(tspan, psi0::Ket{B}, f::Function; kwargs...) where B<:Basis = master_nh_dynamic(tspan, dm(psi0), f; kwargs...)
+for f ∈ [:master,:master_h,:master_nh,:master_dynamic,:master_nh_dynamic]
+    @eval $f(tspan,psi0::Ket,args...;kwargs...) = $f(tspan,dm(psi0),args...;kwargs...)
+end
 
 # Non-hermitian Hamiltonian
 function nh_hamiltonian(H,J,Jdagger,::Nothing)
@@ -207,18 +203,16 @@ function nh_hamiltonian(H,J,Jdagger,rates::AbstractMatrix)
 end
 
 # Recasting needed for the ODE solver is just providing the underlying data
-function recast!(x::T, rho::Operator{B,B,T}) where {B<:Basis,T}
+function recast!(x::T, rho::Operator{B,B,T}) where {B,T}
     rho.data = x
 end
-recast!(rho::Operator{B,B,T}, x::T) where {B<:Basis,T} = nothing
+recast!(rho::Operator{B,B,T}, x::T) where {B,T} = nothing
 
-function integrate_master(tspan, df::Function, rho0,
-                        fout::Union{Nothing, Function}; kwargs...)
-    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
+function integrate_master(tspan, df, rho0, fout; kwargs...)
     x0 = rho0.data
     state = deepcopy(rho0)
     dstate = deepcopy(rho0)
-    integrate(tspan_, df, x0, state, dstate, fout; kwargs...)
+    integrate(tspan, df, x0, state, dstate, fout; kwargs...)
 end
 
 
@@ -231,9 +225,7 @@ end
 # the type of the given decay rate object which can either be nothing, a vector
 # or a matrix.
 
-function dmaster_h(rho::T, H::AbstractOperator{B,B},
-                    rates::Nothing, J::Vector, Jdagger::Vector,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_h(rho, H, rates::Nothing, J, Jdagger, drho, tmp)
     QuantumOpticsBase.mul!(drho,H,rho,-eltype(rho)(im),zero(eltype(rho)))
     QuantumOpticsBase.mul!(drho,rho,H,eltype(rho)(im),one(eltype(rho)))
     for i=1:length(J)
@@ -248,9 +240,7 @@ function dmaster_h(rho::T, H::AbstractOperator{B,B},
     return drho
 end
 
-function dmaster_h(rho::T, H::AbstractOperator{B,B},
-                    rates::Vector, J::Vector, Jdagger::Vector,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_h(rho, H, rates::AbstractVector, J, Jdagger, drho, tmp)
     QuantumOpticsBase.mul!(drho,H,rho,-eltype(rho)(im),zero(eltype(rho)))
     QuantumOpticsBase.mul!(drho,rho,H,eltype(rho)(im),one(eltype(rho)))
     for i=1:length(J)
@@ -265,9 +255,7 @@ function dmaster_h(rho::T, H::AbstractOperator{B,B},
     return drho
 end
 
-function dmaster_h(rho::T, H::AbstractOperator{B,B},
-                    rates::Matrix, J::Vector, Jdagger::Vector,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_h(rho, H, rates::AbstractMatrix, J, Jdagger, drho, tmp)
     QuantumOpticsBase.mul!(drho,H,rho,-eltype(rho)(im),zero(eltype(rho)))
     QuantumOpticsBase.mul!(drho,rho,H,eltype(rho)(im),one(eltype(rho)))
     for j=1:length(J), i=1:length(J)
@@ -282,9 +270,7 @@ function dmaster_h(rho::T, H::AbstractOperator{B,B},
     return drho
 end
 
-function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOperator{B,B},
-                    rates::Nothing, J::Vector, Jdagger::Vector,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_nh(rho, Hnh, Hnh_dagger, rates::Nothing, J, Jdagger, drho, tmp)
     QuantumOpticsBase.mul!(drho,Hnh,rho,-eltype(rho)(im),zero(eltype(rho)))
     QuantumOpticsBase.mul!(drho,rho,Hnh_dagger,eltype(rho)(im),one(eltype(rho)))
     for i=1:length(J)
@@ -294,9 +280,7 @@ function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOper
     return drho
 end
 
-function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOperator{B,B},
-                    rates::Vector, J::Vector, Jdagger::Vector,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_nh(rho, Hnh, Hnh_dagger, rates::AbstractVector, J, Jdagger, drho, tmp)
     QuantumOpticsBase.mul!(drho,Hnh,rho,-eltype(rho)(im),zero(eltype(rho)))
     QuantumOpticsBase.mul!(drho,rho,Hnh_dagger,eltype(rho)(im),one(eltype(rho)))
     for i=1:length(J)
@@ -306,9 +290,7 @@ function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOper
     return drho
 end
 
-function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOperator{B,B},
-                    rates::Matrix, J::Vector, Jdagger::Vector,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_nh(rho, Hnh, Hnh_dagger, rates::AbstractMatrix, J, Jdagger, drho, tmp)
     QuantumOpticsBase.mul!(drho,Hnh,rho,-eltype(rho)(im),zero(eltype(rho)))
     QuantumOpticsBase.mul!(drho,rho,Hnh_dagger,eltype(rho)(im),one(eltype(rho)))
     for j=1:length(J), i=1:length(J)
@@ -323,9 +305,7 @@ function dmaster_liouville(rho,drho,L)
     return drho
 end
 
-function dmaster_h_dynamic(t, rho::T, f::Function,
-                    rates::DecayRates,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_h_dynamic(t, rho, f, rates, drho, tmp)
     result = f(t, rho)
     QO_CHECKS[] && @assert 3 <= length(result) <= 4
     if length(result) == 3
@@ -338,9 +318,7 @@ function dmaster_h_dynamic(t, rho::T, f::Function,
     dmaster_h(rho, H, rates_, J, Jdagger, drho, tmp)
 end
 
-function dmaster_nh_dynamic(t, rho::T, f::Function,
-                    rates::DecayRates,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_nh_dynamic(t, rho, f, rates, drho, tmp)
     result = f(t, rho)
     QO_CHECKS[] && @assert 4 <= length(result) <= 5
     if length(result) == 4
@@ -354,36 +332,35 @@ function dmaster_nh_dynamic(t, rho::T, f::Function,
 end
 
 
-function check_master(rho0::Operator{B,B}, H::AbstractOperator{B,B}, J::Vector, Jdagger::Vector, rates::DecayRates) where B<:Basis
-    # TODO: clean up type checks by dispatch; make type of J known
+function check_master(rho0, H, J, Jdagger, rates)
     isreducible = true # test if all operators are sparse or dense
     if !(isa(H, DenseOpType) || isa(H, SparseOpType))
         isreducible = false
     end
     for j=J
-        @assert isa(j, AbstractOperator{B,B})
+        @assert isa(j, AbstractOperator)
         if !(isa(j, DenseOpType) || isa(j, SparseOpType))
             isreducible = false
         end
         check_samebases(rho0, j)
     end
     for j=Jdagger
-        @assert isa(j, AbstractOperator{B,B})
+        @assert isa(j, AbstractOperator)
         if !(isa(j, DenseOpType) || isa(j, SparseOpType))
             isreducible = false
         end
         check_samebases(rho0, j)
     end
     @assert length(J)==length(Jdagger)
-    if isa(rates,Matrix)
+    if isa(rates,AbstractMatrix)
         @assert size(rates, 1) == size(rates, 2) == length(J)
-    elseif isa(rates,Vector)
+    elseif isa(rates,AbstractVector)
         @assert length(rates) == length(J)
     end
     isreducible
 end
 
 get_type(rho, H, rates::Nothing, J, Jdagger) = promote_type(eltype(rho),eltype(H),eltype.(J)...,eltype.(Jdagger)...)
-get_type(rho, H, rates::Union{Vector,Matrix}, J, Jdagger) = promote_type(eltype(rho),eltype(H),eltype(rates),eltype.(J)...,eltype.(Jdagger)...)
+get_type(rho, H, rates::AbstractVecOrMat, J, Jdagger) = promote_type(eltype(rho),eltype(H),eltype(rates),eltype.(J)...,eltype.(Jdagger)...)
 get_type(rho, Hnh, Hnhdagger, rates::Nothing, J, Jdagger) = promote_type(eltype(rho),eltype(Hnh),eltype(Hnhdagger), eltype.(J)...,eltype.(Jdagger)...)
-get_type(rho, Hnh, Hnhdagger, rates::Union{Vector,Matrix}, J, Jdagger) = promote_type(eltype(rho),eltype(Hnh),eltype(Hnhdagger),eltype(rates),eltype.(J)...,eltype.(Jdagger)...)
+get_type(rho, Hnh, Hnhdagger, rates::AbstractVecOrMat, J, Jdagger) = promote_type(eltype(rho),eltype(Hnh),eltype(Hnhdagger),eltype(rates),eltype.(J)...,eltype.(Jdagger)...)