Support time-dependent operators in timeevolution (#366)

Support time-dependent operators in timeevolution by providing `x_dynamic` methods that take those operators directly, as an alternative to the function-based interface.

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/master.jl

@@ -145,6 +145,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 +161,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 +191,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 +207,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

@@ -157,6 +157,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 +178,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 +204,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

@@ -37,6 +37,10 @@ Integrate time-dependent Schroedinger equation to evolve states or compute propa
  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,
@@ -49,6 +53,11 @@ function schroedinger_dynamic(tspan, psi0, f;
  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}}}
+ schroedinger_dynamic(tspan, psi0, schroedinger_dynamic_function(H); kwargs...)
+end
+
 """
  recast!(x,y)
 

src/time_dependent_operators.jl

@@ -0,0 +1,100 @@
+# Convert storage of heterogeneous stuff to tuples for maximal compilation
+# and to avoid runtime dispatch.
+_tuplify(o::TimeDependentSum) = TimeDependentSum(Tuple, o)
+_tuplify(o::LazySum) = LazySum(eltype(o.factors), o.factors, (o.operators...,))
+_tuplify(o::AbstractOperator) = o
+
+"""
+ schroedinger_dynamic_function(H::AbstractTimeDependentOperator)
+
+Creates a function `f(t, state) -> H(t)`. The `state` argument is ignored.
+
+This is the function expected by [`timeevolution.schroedinger_dynamic()`](@ref).
+"""
+function schroedinger_dynamic_function(H::AbstractTimeDependentOperator)
+ _getfunc(op) = (@inline _tdop_schroedinger_wrapper(t, _) = (op)(t))
+ Htup = _tuplify(H)
+ return _getfunc(Htup)
+end
+
+_tdopdagger(o) = dagger(o)
+function _tdopdagger(o::TimeDependentSum)
+ # This is a kind-of-hacky, more efficient TimeDependentSum dagger operation
+ # that requires that the original operator sticks around and is always
+ # updated first (though this is checked).
+ # Copies and conjugates the coefficients from the original op.
+ o_ls = QuantumOpticsBase.static_operator(o)
+ facs = o_ls.factors
+ c1 = (t)->(@assert current_time(o) == t; conj(facs[1]))
+ crest = (((_)->conj(facs[i])) for i in 2:length(facs))
+ odag = TimeDependentSum((c1, crest...), dagger(o_ls), current_time(o))
+ return odag
+end
+
+"""
+ master_h_dynamic_function(H::AbstractTimeDependentOperator, Js)
+
+Returns a function `f(t, state) -> H(t), Js, dagger.(Js)`.
+The `state` argument is ignored.
+
+This is the function expected by [`timeevolution.master_h_dynamic()`](@ref),
+where `H` is represents the Hamiltonian and `Js` are the (time independent) jump
+operators.
+"""
+function master_h_dynamic_function(H::AbstractTimeDependentOperator, Js)
+ Htup = _tuplify(H)
+ Js_tup = ((_tuplify(J) for J in Js)...,)
+
+ Jdags_tup = _tdopdagger.(Js_tup)
+ function _getfunc(Hop, Jops, Jdops)
+ return (@inline _tdop_master_wrapper_1(t, _) = ((Hop)(t), set_time!.(Jops, t), set_time!.(Jdops, t)))
+ end
+ return _getfunc(Htup, Js_tup, Jdags_tup)
+end
+
+"""
+ master_nh_dynamic_function(Hnh::AbstractTimeDependentOperator, Js)
+
+Returns a function `f(t, state) -> Hnh(t), Hnh(t)', Js, dagger.(Js)`.
+The `state` argument is currently ignored.
+
+This is the function expected by [`timeevolution.master_nh_dynamic()`](@ref),
+where `Hnh` is represents the non-Hermitian Hamiltonian and `Js` are the
+(time independent) jump operators.
+"""
+function master_nh_dynamic_function(Hnh::AbstractTimeDependentOperator, Js)
+ Hnhtup = _tuplify(Hnh)
+ Js_tup = ((_tuplify(J) for J in Js)...,)
+
+ Jdags_tup = _tdopdagger.(Js_tup)
+ Htdagup = _tdopdagger(Hnhtup)
+
+ function _getfunc(Hop, Hdop, Jops, Jdops)
+ return (@inline _tdop_master_wrapper_2(t, _) = ((Hop)(t), (Hdop)(t), set_time!.(Jops, t), set_time!.(Jdops, t)))
+ end
+ return _getfunc(Hnhtup, Htdagup, Js_tup, Jdags_tup)
+end
+
+"""
+ mcfw_dynamic_function(H, Js)
+
+Returns a function `f(t, state) -> H(t), Js, dagger.(Js)`.
+The `state` argument is currently ignored.
+
+This is the function expected by [`timeevolution.mcwf_dynamic()`](@ref),
+where `H` is represents the Hamiltonian and `Js` are the (time independent) jump
+operators.
+"""
+mcfw_dynamic_function(H, Js) = master_h_dynamic_function(H, Js)
+
+"""
+ mcfw_nh_dynamic_function(Hnh, Js)
+
+Returns a function `f(t, state) -> Hnh(t), Js, dagger.(Js)`.
+The `state` argument is currently ignored.
+
+This is the function expected by [`timeevolution.mcwf_dynamic()`](@ref),
+where `Hnh` is represents the non-Hermitian Hamiltonian and `Js` are the (time
+independent) jump operators.
+"""
+mcfw_nh_dynamic_function(Hnh, Js) = master_h_dynamic_function(Hnh, Js)

test/runtests.jl

@@ -6,6 +6,7 @@ names = [
  "test_timeevolution_master.jl",
  "test_timeevolution_mcwf.jl",
  "test_timeevolution_bloch_redfield.jl",
+ "test_timeevolution_tdops.jl",
 
  "test_timeevolution_twolevel.jl",
  "test_timeevolution_pumpedcavity.jl",

test/test_timeevolution_tdops.jl

@@ -0,0 +1,86 @@
+using Test
+using QuantumOptics
+
+@testset "time-dependent operators" begin
+
+b = FockBasis(7)
+
+a = destroy(b)
+
+H0 = number(b)
+Hd = (a + a')
+H = TimeDependentSum(1.0=>H0, cos=>Hd)
+
+ts = [0.0, 0.4]
+ts_half = 0.5 * ts
+
+_h(t, H0, Hd) = H0 + cos(t)*Hd
+_getf = (H0, Hd) -> (t,_) -> _h(t, H0, Hd)
+fman = _getf(H0, Hd)
+
+psi0 = basisstate(b, 1)
+ts_out, psis = timeevolution.schroedinger_dynamic(ts, psi0, H)
+# check this is not trivial
+@test !(psis[1].data ≈ psis[end].data)
+
+ts_out2, psis2 = timeevolution.schroedinger_dynamic(ts, psi0, fman)
+@test psis[end].data ≈ psis2[end].data
+
+_getf = (H0, Hd, a) -> (t,_) -> (_h(t, H0, Hd), (), ())
+fman = _getf(H0, Hd, a)
+
+ts_out, rhos = timeevolution.master_dynamic(ts, psi0, H, [])
+ts_out2, rhos2 = timeevolution.master_dynamic(ts, psi0, fman)
+@test rhos[end].data ≈ rhos2[end].data
+
+Js = [TimeDependentSum(cos=>a), 0.01 * a', 0.01 * LazySum(a' * a)]
+Jdags = dagger.(Js)
+
+_js(t, a) = (cos(t)*a, 0.01*a', 0.01*a'*a)
+
+_getf = (H0, Hd, a) -> (t,_) -> (_h(t, H0, Hd), _js(t, a), dagger.(_js(t, a)))
+fman = _getf(H0, Hd, a)
+
+ts_out, psis = timeevolution.mcwf_dynamic(ts, psi0, H, Js; seed=0)
+ts_out2, psis2 = timeevolution.mcwf_dynamic(ts, psi0, fman; seed=0)
+@test psis[end].data ≈ psis2[end].data
+
+ts_out, rhos = timeevolution.master_dynamic(ts, psi0, H, Js)
+ts_out2, rhos2 = timeevolution.master_dynamic(ts, psi0, fman)
+@test rhos[end].data ≈ rhos2[end].data
+
+Hnh = H - 0.5im * sum(J' * J for J in Js)
+
+_getf = (H0, Hd, a) -> (t,_) -> (
+ _h(t, H0, Hd) - 0.5im * sum(dagger.(_js(t, a)) .* _js(t, a)),
+ _js(t, a),
+ dagger.(_js(t, a)))
+
+fman = _getf(H0, Hd, a)
+
+ts_out, psis = timeevolution.mcwf_nh_dynamic(ts, psi0, Hnh, Js; seed=0)
+ts_out2, psis2 = timeevolution.mcwf_nh_dynamic(ts, psi0, fman; seed=0)
+@test psis[end].data ≈ psis2[end].data
+
+_getf = (H0, Hd, a) -> (t,_) -> (
+ _h(t, H0, Hd) - 0.5im * sum(dagger.(_js(t, a)) .* _js(t, a)),
+ _h(t, H0, Hd) + 0.5im * sum(dagger.(_js(t, a)) .* _js(t, a)),
+ _js(t, a),
+ dagger.(_js(t, a)))
+
+fman = _getf(H0, Hd, a)
+
+ts_out, rhos = timeevolution.master_nh_dynamic(ts, psi0, Hnh, Js)
+ts_out2, rhos2 = timeevolution.master_nh_dynamic(ts, psi0, fman)
+@test rhos[end].data ≈ rhos2[end].data
+
+# for sparse operators, we should not be allocating at each timestep
+allocs1 = @allocated timeevolution.schroedinger_dynamic(ts, psi0, H)
+allocs2 = @allocated timeevolution.schroedinger_dynamic(ts_half, psi0, H)
+@test allocs1 == allocs2
+
+allocs1 = @allocated timeevolution.master_nh_dynamic(ts, psi0, Hnh, Js)
+allocs2 = @allocated timeevolution.master_nh_dynamic(ts_half, psi0, Hnh, Js)
+@test allocs1 == allocs2
+
+end