Merge branch 'master' into compathelper/new_version/2022-11-22-01-09-…

…55-674-03677576554

Project.toml

@@ -1,12 +1,15 @@
 name = "QuantumOptics"
 uuid = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
-version = "v1.0.8"
+version = "1.0.15"
 
 [deps]
 Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
@@ -19,14 +22,17 @@ StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
 
 [compat]
-Arpack = "0.5.1 - 0.5.3, 0.5"
+Arpack = "0.5.1"
+DiffEqBase = "6.113"
 DiffEqCallbacks = "2"
 FFTW = "1"
+ForwardDiff = "0.10"
 IterativeSolvers = "0.9"
+KrylovKit = "0.6"
 LinearMaps = "3"
 OrdinaryDiffEq = "5, 6"
-QuantumOpticsBase = "0.3"
-RecursiveArrayTools = "2"
+QuantumOpticsBase = "0.3, 0.4"
+RecursiveArrayTools = "2, 3"
 Reexport = "0.2, 1.0"
 StochasticDiffEq = "6"
 WignerSymbols = "1, 2"

src/QuantumOptics.jl

@@ -19,6 +19,7 @@ module timeevolution
     export diagonaljumps, @skiptimechecks
 
     include("timeevolution_base.jl")
+    include("time_dependent_operators.jl")
     include("master.jl")
     include("schroedinger.jl")
     include("mcwf.jl")

src/bloch_redfield_master.jl

@@ -16,7 +16,8 @@ See QuTiP's documentation (http://qutip.org/docs/latest/guide/dynamics/dynamics-
                         This argument is only taken into account if use_secular=true.
 """
 function bloch_redfield_tensor(H::AbstractOperator, a_ops; J=SparseOpType[], use_secular=true, secular_cutoff=0.1)
-
+    _check_const(H)
+    _check_const.(J)
     # Use the energy eigenbasis
     H_evals, transf_mat = eigen(Array(H.data)) #Array call makes sure H is a dense array
     H_ekets = [Ket(H.basis_l, transf_mat[:, i]) for i in 1:length(H_evals)]

src/master.jl

@@ -10,6 +10,9 @@ function master_h(tspan, rho0::Operator, H::AbstractOperator, J;
                 Jdagger=dagger.(J),
                 fout=nothing,
                 kwargs...)
+    _check_const(H)
+    _check_const.(J)
+    _check_const.(Jdagger)
     check_master(rho0, H, J, Jdagger, rates)
     tmp = copy(rho0)
     dmaster_(t, rho, drho) = dmaster_h!(drho, H, J, Jdagger, rates, rho, tmp)
@@ -33,6 +36,10 @@ function master_nh(tspan, rho0::Operator, Hnh::AbstractOperator, J;
                 Jdagger=dagger.(J),
                 fout=nothing,
                 kwargs...)
+    _check_const(Hnh)
+    _check_const(Hnhdagger)
+    _check_const.(J)
+    _check_const.(Jdagger)
     check_master(rho0, Hnh, J, Jdagger, rates)
     tmp = copy(rho0)
     dmaster_(t, rho, drho) = dmaster_nh!(drho, Hnh, Hnhdagger, J, Jdagger, rates, rho, tmp)
@@ -76,6 +83,9 @@ function master(tspan, rho0::Operator, H::AbstractOperator, J;
                 Jdagger=dagger.(J),
                 fout=nothing,
                 kwargs...)
+    _check_const(H)
+    _check_const.(J)
+    _check_const.(Jdagger)
     isreducible = check_master(rho0, H, J, Jdagger, rates)
     if !isreducible
         tmp = copy(rho0)
@@ -145,6 +155,12 @@ H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k
 The given function can either be of the form `f(t, rho) -> (Hnh, Hnhdagger, J, Jdagger)`
 or `f(t, rho) -> (Hnh, Hnhdagger, J, Jdagger, rates)` For further information look
 at [`master_dynamic`](@ref).
+
+    timeevolution.master_dynamic(tspan, rho0, Hnh::AbstractTimeDependentOperator, J; <keyword arguments>)
+
+This version takes the non-hermitian Hamiltonian `Hnh` and jump operators `J` as time-dependent operators.
+The jump operators may be `<: AbstractTimeDependentOperator` or other types
+of operator.
 """
 function master_nh_dynamic(tspan, rho0::Operator, f;
                 rates=nothing,
@@ -155,6 +171,12 @@ function master_nh_dynamic(tspan, rho0::Operator, f;
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
 end
 
+function master_nh_dynamic(tspan, rho0::Operator, Hnh::AbstractTimeDependentOperator, J;
+    kwargs...)
+    f = master_nh_dynamic_function(Hnh, J)
+    master_nh_dynamic(tspan, rho0, f; kwargs...)
+end
+
 """
     timeevolution.master_dynamic(tspan, rho0, f; <keyword arguments>)
 
@@ -179,6 +201,12 @@ operators:
         permanent! It is still in use by the ode solver and therefore must not
         be changed.
 * `kwargs...`: Further arguments are passed on to the ode solver.
+
+    timeevolution.master_dynamic(tspan, rho0, H::AbstractTimeDependentOperator, J; <keyword arguments>)
+
+This version takes the Hamiltonian `H` and jump operators `J` as time-dependent operators.
+The jump operators may be `<: AbstractTimeDependentOperator` or other types
+of operator.
 """
 function master_dynamic(tspan, rho0::Operator, f;
                 rates=nothing,
@@ -189,6 +217,11 @@ function master_dynamic(tspan, rho0::Operator, f;
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
 end
 
+function master_dynamic(tspan, rho0::Operator, H::AbstractTimeDependentOperator, J;
+    kwargs...)
+    f = master_h_dynamic_function(H, J)
+    master_dynamic(tspan, rho0, f; kwargs...)
+end
 
 # Automatically convert Ket states to density operators
 for f ∈ [:master,:master_h,:master_nh,:master_dynamic,:master_nh_dynamic]

src/mcwf.jl

@@ -18,6 +18,9 @@ function mcwf_h(tspan, psi0::Ket, H::AbstractOperator, J;
         tmp=copy(psi0),
         display_beforeevent=false, display_afterevent=false,
         kwargs...)
+    _check_const(H)
+    _check_const.(J)
+    _check_const.(Jdagger)
     check_mcwf(psi0, H, J, Jdagger, rates)
     f(t, psi, dpsi) = dmcwf_h!(dpsi, H, J, Jdagger, rates, psi, tmp)
     j(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, rates)
@@ -42,6 +45,8 @@ function mcwf_nh(tspan, psi0::Ket, Hnh::AbstractOperator, J;
         seed=rand(UInt), fout=nothing,
         display_beforeevent=false, display_afterevent=false,
         kwargs...)
+    _check_const(Hnh)
+    _check_const.(J)
     check_mcwf(psi0, Hnh, J, J, nothing)
     f(t, psi, dpsi) = dschroedinger!(dpsi, Hnh, psi)
     j(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, nothing)
@@ -96,6 +101,9 @@ function mcwf(tspan, psi0::Ket, H::AbstractOperator, J;
         fout=nothing, Jdagger=dagger.(J),
         display_beforeevent=false, display_afterevent=false,
         kwargs...)
+    _check_const(H)
+    _check_const.(J)
+    _check_const.(Jdagger)
     isreducible = check_mcwf(psi0, H, J, Jdagger, rates)
     if !isreducible
         tmp = copy(psi0)
@@ -157,6 +165,12 @@ and therefore must not be changed.
 is returned. These correspond to the times at which a jump occured and the index
 of the jump operators with which the jump occured, respectively.
 * `kwargs...`: Further arguments are passed on to the ode solver.
+
+    mcwf_dynamic(tspan, psi0, H::AbstractTimeDependentOperator, J; <keyword arguments>)
+
+This version takes the Hamiltonian `H` and jump operators `J` as time-dependent operators.
+The jump operators may be `<: AbstractTimeDependentOperator` or other types
+of operator.
 """
 function mcwf_dynamic(tspan, psi0::Ket, f;
     seed=rand(UInt), rates=nothing,
@@ -172,8 +186,14 @@ function mcwf_dynamic(tspan, psi0::Ket, f;
         kwargs...)
 end
 
+function mcwf_dynamic(tspan, psi0::Ket, H::AbstractTimeDependentOperator, J; kwargs...)
+    f = mcfw_dynamic_function(H, J)
+    mcwf_dynamic(tspan, psi0, f; kwargs...)
+end
+
 """
     mcwf_nh_dynamic(tspan, rho0, f; <keyword arguments>)
+    mcwf_nh_dynamic(tspan, rho0, Hnh::AbstractTimeDependentOperator, J; <keyword arguments>)
 
 Calculate MCWF trajectory where the dynamic Hamiltonian is given in non-hermitian form.
 
@@ -192,6 +212,11 @@ function mcwf_nh_dynamic(tspan, psi0::Ket, f;
         kwargs...)
 end
 
+function mcwf_nh_dynamic(tspan, psi0::Ket, Hnh::AbstractTimeDependentOperator, J; kwargs...)
+    f = mcfw_nh_dynamic_function(Hnh, J)
+    mcwf_nh_dynamic(tspan, psi0, f; kwargs...)
+end
+
 """
     dmcwf_h_dynamic!(dpsi, f, rates, psi, dpsi_cache, t)
 

src/schroedinger.jl

@@ -5,7 +5,7 @@ Integrate Schroedinger equation to evolve states or compute propagators.
 
 # Arguments
 * `tspan`: Vector specifying the points of time for which output should be displayed.
-* `psi0`: Initial state vector (can be a bra or a ket) or initial propagator.
+* `psi0`: Initial state vector (can be a bra or a ket) or an Operator from some basis to the basis of the Hamiltonian (psi0.basis_l == basis(H)).
 * `H`: Arbitrary operator specifying the Hamiltonian.
 * `fout=nothing`: If given, this function `fout(t, psi)` is called every time
         an output should be displayed. ATTENTION: The state `psi` is neither
@@ -14,8 +14,10 @@ Integrate Schroedinger equation to evolve states or compute propagators.
 """
 function schroedinger(tspan, psi0::T, H::AbstractOperator{B,B};
                 fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B,T<:Union{AbstractOperator{B,B},StateVector{B}}}
+                kwargs...) where {B,Bo,T<:Union{AbstractOperator{B,Bo},StateVector{B}}}
+    _check_const(H)
     dschroedinger_(t, psi, dpsi) = dschroedinger!(dpsi, H, psi)
+    tspan, psi0 = _promote_time_and_state(psi0, H, tspan) # promote only if ForwardDiff.Dual
     x0 = psi0.data
     state = copy(psi0)
     dstate = copy(psi0)
@@ -30,23 +32,39 @@ Integrate time-dependent Schroedinger equation to evolve states or compute propa
 
 # Arguments
 * `tspan`: Vector specifying the points of time for which output should be displayed.
-* `psi0`: Initial state vector (can be a bra or a ket) or initial propagator.
+* `psi0`: Initial state vector (can be a bra or a ket) or an Operator from some basis to the basis of the Hamiltonian (psi0.basis_l == basis(H)).
 * `f`: Function `f(t, psi) -> H` returning the time and or state dependent Hamiltonian.
 * `fout=nothing`: If given, this function `fout(t, psi)` is called every time
         an output should be displayed. ATTENTION: The state `psi` is neither
         normalized nor permanent! It is still in use by the ode solver and
         therefore must not be changed.
+
+    timeevolution.schroedinger_dynamic(tspan, psi0, H::AbstractTimeDependentOperator; fout)
+
+Instead of a function `f`, this takes a time-dependent operator `H`.
 """
 function schroedinger_dynamic(tspan, psi0, f;
                 fout::Union{Function,Nothing}=nothing,
                 kwargs...)
     dschroedinger_(t, psi, dpsi) = dschroedinger_dynamic!(dpsi, f, psi, t)
+    tspan, psi0 = _promote_time_and_state(psi0, f, tspan) # promote only if ForwardDiff.Dual
     x0 = psi0.data
     state = copy(psi0)
     dstate = copy(psi0)
     integrate(tspan, dschroedinger_, x0, state, dstate, fout; kwargs...)
 end
 
+function schroedinger_dynamic(tspan, psi0::T, H::AbstractTimeDependentOperator;
+    kwargs...) where {B,Bp,T<:Union{AbstractOperator{B,Bp},StateVector{B}}}
+    promoted_tspan, psi0 = _promote_time_and_state(psi0, H, tspan)
+    if promoted_tspan !== tspan # promote H
+        promoted_H = TimeDependentSum(H.coefficients, H.static_op.operators; init_time=first(promoted_tspan))
+        return schroedinger_dynamic(promoted_tspan, psi0, schroedinger_dynamic_function(promoted_H); kwargs...)
+    else
+        return schroedinger_dynamic(promoted_tspan, psi0, schroedinger_dynamic_function(H); kwargs...)
+    end
+end
+
 """
     recast!(x,y)
 
@@ -102,3 +120,22 @@ function check_schroedinger(psi::Bra, H)
     check_multiplicable(psi, H)
     check_samebases(H)
 end
+
+
+function _promote_time_and_state(u0, H::AbstractOperator, tspan)
+    Ts = eltype(H)
+    Tt = real(Ts)
+    p = Vector{Tt}(undef,0)
+    u0data_promote = DiffEqBase.promote_u0(u0.data, p, tspan[1])
+    tspan_promote = DiffEqBase.promote_tspan(u0data_promote, p, tspan, nothing, Dict{Symbol, Any}())
+    if u0data_promote !== u0.data
+        u0_promote = rebuild(u0, u0data_promote)
+        return tspan_promote, u0_promote
+    end
+    return tspan_promote, u0
+end
+_promote_time_and_state(u0, f, tspan) = _promote_time_and_state(u0, f(first(tspan), u0), tspan)
+
+rebuild(op::Operator, new_data) = Operator(op.basis_l, op.basis_r, new_data)
+rebuild(state::Ket, new_data) = Ket(state.basis, new_data)
+rebuild(state::Bra, new_data) = Bra(state.basis, new_data)