Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add derivative methods without caches #410

Closed
wants to merge 5 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 12 additions & 0 deletions src/master.jl
Original file line number Diff line number Diff line change
Expand Up @@ -280,6 +280,7 @@ end

"""
dmaster_h!(drho, H, J, Jdagger, rates, rho, drho_cache)
dmaster_h!(drho, H, J, rates, rho)

Update `drho` according to a master equation given in standard Lindblad form.
A cached copy `drho_cache` of `drho` is used as a temporary saving step.
Expand Down Expand Up @@ -332,8 +333,11 @@ function dmaster_h!(drho, H, J, Jdagger, rates::AbstractMatrix, rho, drho_cache)
return drho
end

dmaster_h!(drho, H, J, rates, rho) = dmaster_h!(drho, H, J, dagger.(J), rates, rho, copy(drho))

"""
dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates, rho, drho_cache)
dmaster_nh!(drho, Hnh, J, rates, rho)

Updates `drho` according to a master equation given in standard Lindblad form.
The part of the Liuovillian which can be written as a commutator should be
Expand Down Expand Up @@ -373,6 +377,8 @@ function dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates::AbstractMatrix, r
return drho
end

dmaster_nh!(drho, Hnh, J, rates, rho) = dmaster_nh!(drho, Hnh, dagger(Hnh), J, dagger.(J), rates, rho, copy(drho))

"""
dmaster_liouville!(drho,L,rho)

Expand All @@ -389,6 +395,7 @@ end

"""
dmaster_h_dynamic!(drho, f, rates, rho, drho_cache, t)
dmaster_h_dynamic!(drho, f, rates, rho, t)

Computes the Hamiltonian and jump operators as `H,J,Jdagger=f(t,rho)` and
update `drho` according to a master equation. Optionally, rates can also be
Expand All @@ -410,8 +417,11 @@ function dmaster_h_dynamic!(drho, f::F, rates, rho, drho_cache, t) where {F}
dmaster_h!(drho, H, J, Jdagger, rates_, rho, drho_cache)
end

dmaster_h_dynamic!(drho, f::F, rates, rho, t) where {F} = dmaster_h_dynamic!(drho, f, rates, rho, copy(drho), t)

"""
dmaster_nh_dynamic!(drho, f, rates, rho, drho_cache, t)
dmaster_nh_dynamic!(drho, f, rates, rho, t)

Computes the non-hermitian Hamiltonian and jump operators as
`Hnh,Hnh_dagger,J,Jdagger=f(t,rho)` and update `drho` according to a master
Expand All @@ -433,6 +443,8 @@ function dmaster_nh_dynamic!(drho, f::F, rates, rho, drho_cache, t) where {F}
dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates_, rho, drho_cache)
end

dmaster_nh_dynamic!(drho, f::F, rates, rho, t) where {F} = dmaster_nh_dynamic!(drho, f, rates, rho, copy(drho), t)


function check_master(rho0, H, J, Jdagger, rates)
isreducible = true # test if all operators are sparse or dense
Expand Down
6 changes: 6 additions & 0 deletions src/mcwf.jl
Original file line number Diff line number Diff line change
Expand Up @@ -251,6 +251,7 @@ end

"""
dmcwf_h_dynamic!(dpsi, f, rates, psi, dpsi_cache, t)
dmcwf_h_dynamic!(dpsi, f, rates, psi, t)

Compute the Hamiltonian and jump operators as `H,J,Jdagger=f(t,psi)` and
update `dpsi` according to a non-Hermitian Schrödinger equation.
Expand All @@ -270,6 +271,8 @@ function dmcwf_h_dynamic!(dpsi, f::F, rates, psi, dpsi_cache, t) where {F}
dmcwf_h!(dpsi, H, J, Jdagger, rates_, psi, dpsi_cache)
end

dmcwf_h_dynamic!(dpsi, f::F, rates, psi, t) where {F} = dmcwf_h_dynamic!(dpsi, f, rates, psi, copy(dpsi), t)

"""
dmcwf_nh_dynamic!(dpsi, f, psi, t)

Expand Down Expand Up @@ -532,6 +535,7 @@ end

"""
dmcwf_h!(dpsi, H, J, Jdagger, rates, psi, dpsi_cache)
dmcwf_h!(dpsi, H, J, rates, psi)

Update `dpsi` according to a non-hermitian Schrödinger equation. The
non-hermitian Hamiltonian is given in two parts - the hermitian part H and
Expand All @@ -557,6 +561,8 @@ function dmcwf_h!(dpsi, H, J, Jdagger, rates::AbstractVector, psi, dpsi_cache)
return dpsi
end

dmcwf_h!(dpsi, H, J, rates, psi) = dmcwf_h!(dpsi, H, J, dagger.(J), rates, psi, copy(dpsi))

"""
check_mcwf(psi0, H, J, Jdagger, rates)

Expand Down
12 changes: 11 additions & 1 deletion src/semiclassical.jl
Original file line number Diff line number Diff line change
Expand Up @@ -284,9 +284,10 @@ end

"""
dmaster_h_dynamic!(dstate, fquantum, fclassical!, rates, state, tmp, t)
dmaster_h_dynamic!(dstate, fquantum, fclassical!, rates, state, t)

Update the semiclassical state `dstate` according to a time-dependent,
semiclassical master eqaution.
semiclassical master equation.

See also: [`semiclassical.master_dynamic`](@ref)
"""
Expand All @@ -297,8 +298,13 @@ function dmaster_h_dynamic!(dstate, fquantum::F, fclassical!::G, rates, state, t
return dstate
end

function dmaster_h_dynamic!(dstate, fquantum::F, fclassical!::G, rates, state, t) where {F,G}
dmaster_h_dynamic!(dstate, fquantum, fclassical!, rates, state, copy(dstate.quantum), t)
end

"""
dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, tmp, t)
dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, t)

Update the semiclassical state `dpsi` according to a time-dependent, semiclassical
and non-Hermitian Schrödinger equation (MCWF).
Expand All @@ -312,6 +318,10 @@ function dmcwf_h_dynamic!(dpsi, fquantum::F, fclassical!::G, rates, psi, tmp, t)
return dpsi
end

function dmcwf_h_dynamic!(dpsi, fquantum::F, fclassical!::G, rates, psi, t) where {F,G}
dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, copy(dpsi.quantum), t)
end

function jump_dynamic(rng, t, psi, fquantum::F, fclassical!::G, fjump_classical!::H, psi_new, probs_tmp, rates) where {F,G,H}
result = fquantum(t, psi.quantum, psi.classical)
QO_CHECKS[] && @assert 3 <= length(result) <= 4
Expand Down
14 changes: 14 additions & 0 deletions test/test_semiclassical.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@ using Test
using QuantumOptics
using LinearAlgebra
using QuantumOpticsBase: IncompatibleBases
using OrdinaryDiffEq

@testset "semiclassical" begin

Expand Down Expand Up @@ -163,6 +164,19 @@ tout5, ut = semiclassical.mcwf_dynamic(T,ψ0,fquantum,fclassical,fjump_classical

@test ψt1[end].classical[2] == u0[2] + 1.0

# Test no cache derivative methods

t0, t1 = (0.0, 1.0)
rho0 = semiclassical.State(dm(spinup(ba)⊗fockstate(bf,0)),u0)

fmaster_h_dynamic!(dstate, state, p, t) = semiclassical.dmaster_h_dynamic!(dstate, fquantum_master, fclassical, nothing, state, t)
prob_master_h_dynamic! = ODEProblem(fmaster_h_dynamic!, rho0, (t0, t1))
@test_nowarn sol_master_h_dynamic! = solve(prob_master_h_dynamic!, DP5(); save_everystep=false)

fmcwf_h_dynamic!(dstate, state, p, t) = semiclassical.dmcwf_h_dynamic!(dstate, fquantum, fclassical, nothing, state, t)
prob_mcwf_h_dynamic! = ODEProblem(fmcwf_h_dynamic!, ψ0, (t0, t1))
@test_nowarn sol_mcwf_h_dynamic! = solve(prob_mcwf_h_dynamic!, DP5(); save_everystep=false)

# Test classical jump behavior
before_jump = findfirst(t -> !(t∈T), tout3)
after_jump = findlast(t-> !(t∈T), tout4)
Expand Down
26 changes: 26 additions & 0 deletions test/test_timeevolution_master.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
using Test
using QuantumOptics
using LinearAlgebra
using OrdinaryDiffEq

@testset "master" begin

Expand Down Expand Up @@ -38,6 +39,13 @@ J = [Ja, Jc]
Jlazy = [LazyTensor(basis, 1, sqrt(γ)*sm), Jc]

Hnh = H - 0.5im*sum([dagger(J[i])*J[i] for i=1:length(J)])
function Ht(t, psi)
return H*exp(-(5-t)^2), J, dagger.(J)
end
function Hnht(t, psi)
Hnhtime = H*exp(-(5-t)^2) - 0.5im * sum(i' * i for i in J)
return Hnhtime, dagger(Hnhtime), J, dagger.(J)
end

Hdense = dense(H)
Hlazy = LazySum(Ha, Hc, Hint)
Expand Down Expand Up @@ -109,6 +117,24 @@ tout, ρt = timeevolution.master_nh(T, ρ₀, Hnh_dense, J; reltol=1e-6)
tout, ρt = timeevolution.master_nh(T, Ψ₀, Hnh_dense, Jdense; reltol=1e-6)
@test tracedistance(ρt[end], ρ) < 1e-5

# Test no cache derivative methods

t0, t1 = (0.0, 1.0)
fmaster_h!(drho, rho, p, t) = timeevolution.dmaster_h!(drho, H, J, nothing, rho)
prob_master_h! = ODEProblem(fmaster_h!, ρ₀, (t0, t1))
@test_nowarn sol_master_h! = solve(prob_master_h!, DP5(); save_everystep=false)

fmaster_nh!(drho, rho, p, t) = timeevolution.dmaster_nh!(drho, Hnh, J, nothing, rho)
prob_master_nh! = ODEProblem(fmaster_nh!, ρ₀, (t0, t1))
@test_nowarn sol_master_nh! = solve(prob_master_nh!, DP5(); save_everystep=false)

fmaster_h_dynamic!(drho, rho, p, t) = timeevolution.dmaster_h_dynamic!(drho, Ht, nothing, rho, t)
prob_master_h_dynamic! = ODEProblem(fmaster_h_dynamic!, ρ₀, (t0, t1))
@test_nowarn sol_master_h_dynamic! = solve(prob_master_h_dynamic!, DP5(); save_everystep=false)

fmaster_nh_dynamic!(drho, rho, p, t) = timeevolution.dmaster_nh_dynamic!(drho, Hnht, nothing, rho, t)
prob_master_nh_dynamic! = ODEProblem(fmaster_nh_dynamic!, ρ₀, (t0, t1))
@test_nowarn sol_master_nh_dynamic! = solve(prob_master_nh_dynamic!, DP5(); save_everystep=false)

# Test explicit gamma vector
rates_vector = [γ, κ]
Expand Down
15 changes: 15 additions & 0 deletions test/test_timeevolution_mcwf.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
using Test
using QuantumOptics
using Random, LinearAlgebra
using OrdinaryDiffEq

@testset "mcwf" begin

Expand Down Expand Up @@ -115,6 +116,20 @@ tout, Ψt = timeevolution.mcwf_nh(T, Ψ₀, Hnh_dense, J; seed=UInt(1), reltol=1
tout, Ψt = timeevolution.mcwf_nh(T, Ψ₀, Hnh, J; seed=UInt(2), reltol=1e-6)
@test norm(Ψt[end]-Ψ) > 0.1

# Test no cache derivative methods

function Ht(t, psi)
H*exp(-(5-t)^2), J, dagger.(J)
end
t0, t1 = (0.0, 1.0)

fmcwf_h_dynamic!(dpsi, psi, p, t) = timeevolution.dmcwf_h_dynamic!(dpsi, Ht, nothing, psi, t)
prob_mcwf_h_dynamic! = ODEProblem(fmcwf_h_dynamic!, Ψ₀, (t0, t1))
@test_nowarn sol_mcwf_h_dynamic! = solve(prob_mcwf_h_dynamic!, DP5(); save_everystep=false)

fmcwf_h!(dpsi, psi, p, t) = timeevolution.dmcwf_h!(dpsi, H, J, nothing, psi)
prob_mcwf_h! = ODEProblem(fmcwf_h!, Ψ₀, (t0, t1))
@test_nowarn sol_mcwf_h! = solve(prob_mcwf_h!, DP5(); save_everystep=false)


# Test convergence to master solution
Expand Down