From f46437a98d4cdfe836dbbb7e180f87eaa5421cfd Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 9 Jan 2018 07:32:32 -0800 Subject: [PATCH 01/10] use DifferentialEquations --- src/timeevolution_base.jl | 29 ++++++++++++++++++------- test/test_timeevolution_schroedinger.jl | 7 ++---- 2 files changed, 23 insertions(+), 13 deletions(-) diff --git a/src/timeevolution_base.jl b/src/timeevolution_base.jl index 6ec3e5d4..8590dae0 100644 --- a/src/timeevolution_base.jl +++ b/src/timeevolution_base.jl @@ -1,5 +1,7 @@ using ..ode_dopri +import OrdinaryDiffEq, DiffEqCallbacks + function recast! end """ @@ -13,21 +15,32 @@ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex12 df(t, state, dstate) recast!(dstate, dx) end - function fout_(t::Float64, x::Vector{Complex128}) + function fout_(t::Float64, x::Vector{Complex128},integrator) recast!(x, state) fout(t, state) end - ode(df_, tspan, x0, fout_; kwargs...) + + # TODO: Infer the output of `fout` instead of computing it + recast!(x0, state) + out = DiffEqCallbacks.SavedValues(Float64,typeof(fout(tspan[1], state))) + + # Build callback solve with DP5 + # TODO: Expose algorithm choice + cb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan) + sol = OrdinaryDiffEq.solve( + OrdinaryDiffEq.ODEProblem(df_, x0,(tspan[1],tspan[end])), + OrdinaryDiffEq.DP5(); + reltol = 1.0e-6, + abstol = 1.0e-8, + save_everystep = false, save_start = false, save_end = false, + callback=cb, kwargs...) + out.t,out.saveval end function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128}, state::T, dstate::T, ::Void; kwargs...) - tout = Float64[] - xout = T[] function fout(t::Float64, state::T) - push!(tout, t) - push!(xout, copy(state)) + copy(state) end integrate(tspan, df, x0, state, dstate, fout; kwargs...) - (tout, xout) -end \ No newline at end of file +end diff --git a/test/test_timeevolution_schroedinger.jl b/test/test_timeevolution_schroedinger.jl index b072b2ef..893bddc6 100644 --- a/test/test_timeevolution_schroedinger.jl +++ b/test/test_timeevolution_schroedinger.jl @@ -61,13 +61,10 @@ for (i, t) in enumerate(tout) @test tracedistance(rho_rot, full(R)*rho*dagger(full(R))) < 1e-5 end -t_fout = Float64[] -psi_fout = [] function fout(t, psi) - push!(t_fout, t) - push!(psi_fout, deepcopy(psi)) + deepcopy(psi) end -timeevolution.schroedinger(T, psi0, Hrot; fout=fout) +t_fout, psi_fout = timeevolution.schroedinger(T, psi0, Hrot; fout=fout) @test t_fout == tout && psi_fout == psi_rot_t end # testset From 402604bb9ac05d4e4cef1a0230aa0c915dbe3d10 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 9 Jan 2018 08:22:44 -0800 Subject: [PATCH 02/10] force in place --- src/timeevolution_base.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/timeevolution_base.jl b/src/timeevolution_base.jl index 8590dae0..bce499a1 100644 --- a/src/timeevolution_base.jl +++ b/src/timeevolution_base.jl @@ -28,7 +28,7 @@ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex12 # TODO: Expose algorithm choice cb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan) sol = OrdinaryDiffEq.solve( - OrdinaryDiffEq.ODEProblem(df_, x0,(tspan[1],tspan[end])), + OrdinaryDiffEq.ODEProblem{true}(df_, x0,(tspan[1],tspan[end])), OrdinaryDiffEq.DP5(); reltol = 1.0e-6, abstol = 1.0e-8, From e818c79ec6e19c5880f280205e8522d46f3e534d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 9 Jan 2018 11:56:01 -0800 Subject: [PATCH 03/10] steady state attempt --- src/steadystate.jl | 50 ++++++++++++++++++--------------------- src/timeevolution_base.jl | 11 +++++++-- 2 files changed, 32 insertions(+), 29 deletions(-) diff --git a/src/steadystate.jl b/src/steadystate.jl index a9316e5f..12cd47c2 100644 --- a/src/steadystate.jl +++ b/src/steadystate.jl @@ -2,7 +2,7 @@ module steadystate using ..states, ..operators, ..operators_dense, ..superoperators using ..timeevolution, ..metrics - +import OrdinaryDiffEq type ConvergenceReached <: Exception end @@ -36,35 +36,31 @@ function master(H::Operator, J::Vector; Jdagger::Vector=dagger.(J), fout::Union{Function,Void}=nothing, kwargs...) - t0 = 0. - rho0 = copy(rho0) - function fout_steady(t, rho) - if fout!=nothing - fout(t, rho) - end - dt = t - t0 - drho = metrics.tracedistance(rho0, rho) - t0 = t - rho0.data[:] = rho.data - if drho/dt < eps - throw(ConvergenceReached()) - end - end - try - timeevolution.master([0., Inf], rho0, H, J; rates=rates, Jdagger=Jdagger, - hmin=hmin, hmax=Inf, - display_initialvalue=false, - display_finalvalue=false, - display_intermediatesteps=true, - fout=fout_steady, kwargs...) - catch e - if !isa(e, ConvergenceReached) - rethrow(e) - end - end + cb = OrdinaryDiffEq.DiscreteCallback(SteadyStateCondtion(rho0,eps), + (integrator)->OrdinaryDiffEq.terminate!(integrator), + save_positions=(false,false)) + + timeevolution.master([0., Inf], rho0, H, J; rates=rates, Jdagger=Jdagger, + hmin=hmin, hmax=Inf, + display_initialvalue=false, + display_finalvalue=false, + display_intermediatesteps=true, + fout=fout, + callback = cb, kwargs...) return rho0 end +struct SteadyStateCondtion{T,T2} + rho0::T + eps::T2 +end +function (c::SteadyStateCondtion)(t,rho,integrator) + dt = integrator.dt + drho = metrics.tracedistance(c.rho0, rho) + c.rho0.data[:] = rho.data + drho/dt < c.eps +end + """ steadystate.eigenvector(L) steadystate.eigenvector(H, J) diff --git a/src/timeevolution_base.jl b/src/timeevolution_base.jl index bce499a1..e7950d1e 100644 --- a/src/timeevolution_base.jl +++ b/src/timeevolution_base.jl @@ -8,7 +8,8 @@ function recast! end df(t, state::T, dstate::T) """ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128}, - state::T, dstate::T, fout::Function; kwargs...) + state::T, dstate::T, fout::Function; callback = nothing, kwargs...) + function df_(t, x::Vector{Complex128}, dx::Vector{Complex128}) recast!(x, state) recast!(dx, dstate) @@ -27,13 +28,19 @@ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex12 # Build callback solve with DP5 # TODO: Expose algorithm choice cb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan) + + if callback == nothing + _cb = cb + else + _cb = OrdinaryDiffEq.CallbackSet(cb,callback) + end sol = OrdinaryDiffEq.solve( OrdinaryDiffEq.ODEProblem{true}(df_, x0,(tspan[1],tspan[end])), OrdinaryDiffEq.DP5(); reltol = 1.0e-6, abstol = 1.0e-8, save_everystep = false, save_start = false, save_end = false, - callback=cb, kwargs...) + callback=_cb, kwargs...) out.t,out.saveval end From 8d1c17fb3db17e2398aff0769fb124257dc03fcc Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 9 Jan 2018 12:32:05 -0800 Subject: [PATCH 04/10] steady state works --- src/steadystate.jl | 30 +++++---------------- src/timeevolution_base.jl | 56 +++++++++++++++++++++++++++++---------- 2 files changed, 49 insertions(+), 37 deletions(-) diff --git a/src/steadystate.jl b/src/steadystate.jl index 12cd47c2..542b202f 100644 --- a/src/steadystate.jl +++ b/src/steadystate.jl @@ -1,10 +1,7 @@ module steadystate using ..states, ..operators, ..operators_dense, ..superoperators -using ..timeevolution, ..metrics -import OrdinaryDiffEq - -type ConvergenceReached <: Exception end +using ..timeevolution """ @@ -31,34 +28,21 @@ Calculate steady state using long time master equation evolution. """ function master(H::Operator, J::Vector; rho0::DenseOperator=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))), - eps::Float64=1e-3, hmin=1e-7, + hmin=1e-7, rates::Union{Vector{Float64}, Matrix{Float64}, Void}=nothing, Jdagger::Vector=dagger.(J), fout::Union{Function,Void}=nothing, kwargs...) - cb = OrdinaryDiffEq.DiscreteCallback(SteadyStateCondtion(rho0,eps), - (integrator)->OrdinaryDiffEq.terminate!(integrator), - save_positions=(false,false)) - - timeevolution.master([0., Inf], rho0, H, J; rates=rates, Jdagger=Jdagger, + t,u = timeevolution.master([0., Inf], rho0, H, J; rates=rates, Jdagger=Jdagger, hmin=hmin, hmax=Inf, display_initialvalue=false, display_finalvalue=false, display_intermediatesteps=true, fout=fout, - callback = cb, kwargs...) - return rho0 -end - -struct SteadyStateCondtion{T,T2} - rho0::T - eps::T2 -end -function (c::SteadyStateCondtion)(t,rho,integrator) - dt = integrator.dt - drho = metrics.tracedistance(c.rho0, rho) - c.rho0.data[:] = rho.data - drho/dt < c.eps + steady_state = true, + eps = eps, kwargs...) + timeevolution.recast!(u[end],rho0) + rho0 end """ diff --git a/src/timeevolution_base.jl b/src/timeevolution_base.jl index e7950d1e..b82f035b 100644 --- a/src/timeevolution_base.jl +++ b/src/timeevolution_base.jl @@ -1,4 +1,4 @@ -using ..ode_dopri +using ..ode_dopri, ..metrics import OrdinaryDiffEq, DiffEqCallbacks @@ -8,7 +8,8 @@ function recast! end df(t, state::T, dstate::T) """ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128}, - state::T, dstate::T, fout::Function; callback = nothing, kwargs...) + state::T, dstate::T, fout::Function; + steady_state = false, eps = 1e-3, kwargs...) function df_(t, x::Vector{Complex128}, dx::Vector{Complex128}) recast!(x, state) @@ -24,24 +25,38 @@ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex12 # TODO: Infer the output of `fout` instead of computing it recast!(x0, state) out = DiffEqCallbacks.SavedValues(Float64,typeof(fout(tspan[1], state))) + scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan) # Build callback solve with DP5 # TODO: Expose algorithm choice - cb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan) - if callback == nothing - _cb = cb + prob = OrdinaryDiffEq.ODEProblem{true}(df_, x0,(tspan[1],tspan[end])) + + if steady_state + _cb = OrdinaryDiffEq.DiscreteCallback( + SteadyStateCondtion(copy(state),eps,state), + (integrator)->OrdinaryDiffEq.terminate!(integrator), + save_positions = (false,false)) + cb = OrdinaryDiffEq.CallbackSet(_cb,scb) + sol = OrdinaryDiffEq.solve( + prob, + OrdinaryDiffEq.DP5(); + reltol = 1.0e-6, + abstol = 1.0e-8, + save_everystep = false, save_start = false, + callback=cb, kwargs...) + # TODO: On v0.7 it's type-stable to return only sol.u[end]! + return sol.t,sol.u else - _cb = OrdinaryDiffEq.CallbackSet(cb,callback) + sol = OrdinaryDiffEq.solve( + prob, + OrdinaryDiffEq.DP5(); + reltol = 1.0e-6, + abstol = 1.0e-8, + save_everystep = false, save_start = false, save_end = false, + callback=scb, kwargs...) + return out.t,out.saveval end - sol = OrdinaryDiffEq.solve( - OrdinaryDiffEq.ODEProblem{true}(df_, x0,(tspan[1],tspan[end])), - OrdinaryDiffEq.DP5(); - reltol = 1.0e-6, - abstol = 1.0e-8, - save_everystep = false, save_start = false, save_end = false, - callback=_cb, kwargs...) - out.t,out.saveval end function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128}, @@ -51,3 +66,16 @@ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex12 end integrate(tspan, df, x0, state, dstate, fout; kwargs...) end + +struct SteadyStateCondtion{T,T2,T3} + rho0::T + eps::T2 + state::T3 +end +function (c::SteadyStateCondtion)(t,rho,integrator) + timeevolution.recast!(rho,c.state) + dt = integrator.dt + drho = metrics.tracedistance(c.rho0, c.state) + c.rho0.data[:] = c.state.data + drho/dt < c.eps +end From 8d306d350b2b5e529e038d90c550d5a62b6f99a7 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 9 Jan 2018 13:42:30 -0800 Subject: [PATCH 05/10] steady state passes --- REQUIRE | 2 ++ src/steadystate.jl | 2 -- src/timecorrelations.jl | 20 ++++++++---------- src/timeevolution_base.jl | 43 +++++++++++++++++++++------------------ test/test_steadystate.jl | 8 ++++---- 5 files changed, 37 insertions(+), 38 deletions(-) diff --git a/REQUIRE b/REQUIRE index 7d394851..73219631 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,3 +1,5 @@ julia 0.5 Compat 0.17 Roots +OrdinaryDiffEq 2.34.0 +DiffEqCallbacks 0.7.0 diff --git a/src/steadystate.jl b/src/steadystate.jl index 542b202f..b65a0280 100644 --- a/src/steadystate.jl +++ b/src/steadystate.jl @@ -41,8 +41,6 @@ function master(H::Operator, J::Vector; fout=fout, steady_state = true, eps = eps, kwargs...) - timeevolution.recast!(u[end],rho0) - rho0 end """ diff --git a/src/timecorrelations.jl b/src/timecorrelations.jl index 5208126b..0aab32d7 100644 --- a/src/timecorrelations.jl +++ b/src/timecorrelations.jl @@ -36,13 +36,12 @@ function correlation(tspan::Vector{Float64}, rho0::DenseOperator, H::Operator, J rates::Union{Vector{Float64}, Matrix{Float64}, Void}=nothing, Jdagger::Vector=dagger.(J), kwargs...) - exp_values = Complex128[] function fout(t, rho) - push!(exp_values, expect(op1, rho)) + expect(op1, rho) end - timeevolution.master(tspan, op2*rho0, H, J; rates=rates, Jdagger=Jdagger, + t,u = timeevolution.master(tspan, op2*rho0, H, J; rates=rates, Jdagger=Jdagger, fout=fout, kwargs...) - return exp_values + u end function correlation(rho0::DenseOperator, H::Operator, J::Vector, @@ -52,15 +51,12 @@ function correlation(rho0::DenseOperator, H::Operator, J::Vector, Jdagger::Vector=dagger.(J), kwargs...) op2rho0 = op2*rho0 - tout = Float64[0.] - exp_values = Complex128[expect(op1, op2rho0)] + exp1 = expect(op1, op2rho0) function fout(t, rho) - push!(tout, t) - push!(exp_values, expect(op1, rho)) + expect(op1, rho) end - steadystate.master(H, J; rho0=op2rho0, eps=eps, h0=h0, fout=fout, - rates=rates, Jdagger=Jdagger, kwargs...) - return tout, exp_values + t,u = steadystate.master(H, J; rho0=op2rho0, eps=eps, h0=h0, fout=fout, + rates=rates, Jdagger=Jdagger, save_everystep=true,kwargs...) end @@ -115,7 +111,7 @@ end function spectrum(H::Operator, J::Vector, op::Operator; rho0::DenseOperator=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))), eps::Float64=1e-4, h0=10., - rho_ss::DenseOperator=steadystate.master(H, J; eps=eps), + rho_ss::DenseOperator=steadystate.master(H, J; eps=eps)[end][end], kwargs...) tspan, exp_values = correlation(rho_ss, H, J, dagger(op), op, eps=eps, h0=h0, kwargs...) dtmin = minimum(diff(tspan)) diff --git a/src/timeevolution_base.jl b/src/timeevolution_base.jl index b82f035b..ea5ab1ac 100644 --- a/src/timeevolution_base.jl +++ b/src/timeevolution_base.jl @@ -9,7 +9,8 @@ df(t, state::T, dstate::T) """ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128}, state::T, dstate::T, fout::Function; - steady_state = false, eps = 1e-3, kwargs...) + steady_state = false, eps = 1e-3, save_everystep = false, + kwargs...) function df_(t, x::Vector{Complex128}, dx::Vector{Complex128}) recast!(x, state) @@ -25,7 +26,12 @@ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex12 # TODO: Infer the output of `fout` instead of computing it recast!(x0, state) out = DiffEqCallbacks.SavedValues(Float64,typeof(fout(tspan[1], state))) - scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan) + + scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan, + save_everystep=save_everystep, + save_start = false) + + # Build callback solve with DP5 # TODO: Expose algorithm choice @@ -33,30 +39,27 @@ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex12 prob = OrdinaryDiffEq.ODEProblem{true}(df_, x0,(tspan[1],tspan[end])) if steady_state + affect! = function (integrator) + !save_everystep && scb.affect!(integrator,true) + OrdinaryDiffEq.terminate!(integrator) + end _cb = OrdinaryDiffEq.DiscreteCallback( SteadyStateCondtion(copy(state),eps,state), - (integrator)->OrdinaryDiffEq.terminate!(integrator), + affect!; save_positions = (false,false)) cb = OrdinaryDiffEq.CallbackSet(_cb,scb) - sol = OrdinaryDiffEq.solve( - prob, - OrdinaryDiffEq.DP5(); - reltol = 1.0e-6, - abstol = 1.0e-8, - save_everystep = false, save_start = false, - callback=cb, kwargs...) - # TODO: On v0.7 it's type-stable to return only sol.u[end]! - return sol.t,sol.u else - sol = OrdinaryDiffEq.solve( - prob, - OrdinaryDiffEq.DP5(); - reltol = 1.0e-6, - abstol = 1.0e-8, - save_everystep = false, save_start = false, save_end = false, - callback=scb, kwargs...) - return out.t,out.saveval + cb = scb end + sol = OrdinaryDiffEq.solve( + prob, + OrdinaryDiffEq.DP5(); + reltol = 1.0e-6, + abstol = 1.0e-8, + save_everystep = false, save_start = false, + save_end = false, + callback=cb, kwargs...) + out.t,out.saveval end function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128}, diff --git a/test/test_steadystate.jl b/test/test_steadystate.jl index 3563a08b..c7ec3b4d 100644 --- a/test/test_steadystate.jl +++ b/test/test_steadystate.jl @@ -42,8 +42,8 @@ Jdense = map(full, J) tout, ρt = timeevolution.master([0,100], ρ₀, Hdense, Jdense; reltol=1e-7) -ρss = steadystate.master(Hdense, Jdense; eps=1e-4) -@test tracedistance(ρss, ρt[end]) < 1e-3 +tss, ρss = steadystate.master(Hdense, Jdense; eps=1e-4) +@test tracedistance(ρss[end], ρt[end]) < 1e-3 ρss = steadystate.eigenvector(Hdense, Jdense) @test tracedistance(ρss, ρt[end]) < 1e-6 @@ -57,8 +57,8 @@ Hp = η*(destroy(fockbasis) + create(fockbasis)) Jp = [sqrt(2κ)*destroy(fockbasis)] n_an = η^2/κ^2 -ρss = steadystate.master(Hp, Jp; rho0=ρ0_p, eps=1e-4) -nss = expect(create(fockbasis)*destroy(fockbasis), ρss) +tss,ρss = steadystate.master(Hp, Jp; rho0=ρ0_p, eps=1e-4) +nss = expect(create(fockbasis)*destroy(fockbasis), ρss[end]) @test n_an - real(nss) < 1e-3 ρss = steadystate.eigenvector(Hp, Jp) From 411a1fcba3299cf5116a4d2f73b49b406583d5ab Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 11 Jan 2018 06:21:16 -0800 Subject: [PATCH 06/10] change compatibilities --- .travis.yml | 1 - REQUIRE | 2 +- appveyor.yml | 2 -- 3 files changed, 1 insertion(+), 4 deletions(-) diff --git a/.travis.yml b/.travis.yml index 9fa9442c..87547890 100644 --- a/.travis.yml +++ b/.travis.yml @@ -3,7 +3,6 @@ os: - osx - linux julia: - - 0.5 - 0.6 - nightly matrix: diff --git a/REQUIRE b/REQUIRE index 73219631..18ecbe61 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,4 +1,4 @@ -julia 0.5 +julia 0.6 Compat 0.17 Roots OrdinaryDiffEq 2.34.0 diff --git a/appveyor.yml b/appveyor.yml index 20c6794b..3d66d031 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -1,7 +1,5 @@ environment: matrix: - - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.5/julia-0.5-latest-win32.exe" - - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.5/julia-0.5-latest-win64.exe" - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.6/julia-0.6-latest-win32.exe" - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.6/julia-0.6-latest-win64.exe" - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe" From ea937ce4a577da11de8159e6cc9e644349734a90 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 11 Jan 2018 06:25:19 -0800 Subject: [PATCH 07/10] update readme --- README.md | 8 -------- 1 file changed, 8 deletions(-) diff --git a/README.md b/README.md index 708f4202..de463189 100644 --- a/README.md +++ b/README.md @@ -8,8 +8,6 @@ More information, documentation and examples can be found on our website http:// ### Development status **Latest release:** - * [![Status of latest release on julia 0.4][pkg-0.4-img]][pkg-0.4-url] - * [![Status of latest release on julia 0.5][pkg-0.5-img]][pkg-0.5-url] * [![Status of latest release on julia 0.6][pkg-0.6-img]][pkg-0.6-url] **Master branch:** @@ -55,12 +53,6 @@ Additionally, you can contact us on Gitter: [![Chat on Gitter][gitter-img]][gitt [codecov-url]: https://codecov.io/gh/qojulia/QuantumOptics.jl [codecov-img]: https://codecov.io/gh/qojulia/QuantumOptics.jl/branch/master/graph/badge.svg -[pkg-0.4-url]: http://pkg.julialang.org/?pkg=QuantumOptics&ver=0.4 -[pkg-0.4-img]: http://pkg.julialang.org/badges/QuantumOptics_0.4.svg - -[pkg-0.5-url]: http://pkg.julialang.org/?pkg=QuantumOptics&ver=0.5 -[pkg-0.5-img]: http://pkg.julialang.org/badges/QuantumOptics_0.5.svg - [pkg-0.6-url]: http://pkg.julialang.org/?pkg=QuantumOptics&ver=0.6 [pkg-0.6-img]: http://pkg.julialang.org/badges/QuantumOptics_0.6.svg From 4f0f675f6bc266c8e68372db83450be02842004f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 11 Jan 2018 07:50:08 -0800 Subject: [PATCH 08/10] core inference and algorithm passing --- src/timeevolution_base.jl | 21 ++++++++++++--------- test/test_steadystate.jl | 2 +- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/src/timeevolution_base.jl b/src/timeevolution_base.jl index ea5ab1ac..e3a25eba 100644 --- a/src/timeevolution_base.jl +++ b/src/timeevolution_base.jl @@ -8,7 +8,8 @@ function recast! end df(t, state::T, dstate::T) """ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128}, - state::T, dstate::T, fout::Function; + state::T, dstate::T, fout::Function, + alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.DP5(); steady_state = false, eps = 1e-3, save_everystep = false, kwargs...) @@ -25,17 +26,15 @@ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex12 # TODO: Infer the output of `fout` instead of computing it recast!(x0, state) - out = DiffEqCallbacks.SavedValues(Float64,typeof(fout(tspan[1], state))) + + out_type = pure_inference(fout, Tuple{eltype(tspan),typeof(state)}) + + out = DiffEqCallbacks.SavedValues(Float64,out_type) scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan, save_everystep=save_everystep, save_start = false) - - - # Build callback solve with DP5 - # TODO: Expose algorithm choice - prob = OrdinaryDiffEq.ODEProblem{true}(df_, x0,(tspan[1],tspan[end])) if steady_state @@ -51,9 +50,11 @@ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex12 else cb = scb end + + # TODO: Expose algorithm choice sol = OrdinaryDiffEq.solve( prob, - OrdinaryDiffEq.DP5(); + alg; reltol = 1.0e-6, abstol = 1.0e-8, save_everystep = false, save_start = false, @@ -63,7 +64,7 @@ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex12 end function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128}, - state::T, dstate::T, ::Void; kwargs...) + state::T, dstate::T, ::Void, args...; kwargs...) function fout(t::Float64, state::T) copy(state) end @@ -82,3 +83,5 @@ function (c::SteadyStateCondtion)(t,rho,integrator) c.rho0.data[:] = c.state.data drho/dt < c.eps end + +Base.@pure pure_inference(fout,T) = Core.Inference.return_type(fout, T) diff --git a/test/test_steadystate.jl b/test/test_steadystate.jl index c7ec3b4d..242bffc9 100644 --- a/test/test_steadystate.jl +++ b/test/test_steadystate.jl @@ -75,6 +75,6 @@ nss = expect(create(fockbasis)*destroy(fockbasis), ρss) function fout_wrong(t, x) @assert x == t end -@test_throws AssertionError steadystate.master(Hdense, Jdense; fout=fout_wrong) +@test_throws MethodError steadystate.master(Hdense, Jdense; fout=fout_wrong) end # testset From 8c4c3e9dbfdb8c6ec8e251103b3767cfa9e3f928 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 11 Jan 2018 07:53:51 -0800 Subject: [PATCH 09/10] update REQUIRE for new DiffEqCallbacks --- REQUIRE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/REQUIRE b/REQUIRE index 18ecbe61..977ff55e 100644 --- a/REQUIRE +++ b/REQUIRE @@ -2,4 +2,4 @@ julia 0.6 Compat 0.17 Roots OrdinaryDiffEq 2.34.0 -DiffEqCallbacks 0.7.0 +DiffEqCallbacks 0.7.1 From b5c73c808fc7dea7904153d4958ba2354ca75691 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 11 Jan 2018 09:00:46 -0800 Subject: [PATCH 10/10] fix mcwf --- src/mcwf.jl | 45 ++++++++++++++++++++------------------- src/timeevolution_base.jl | 12 ++++++----- 2 files changed, 30 insertions(+), 27 deletions(-) diff --git a/src/mcwf.jl b/src/mcwf.jl index bea37f0c..dcd26065 100644 --- a/src/mcwf.jl +++ b/src/mcwf.jl @@ -4,7 +4,8 @@ export mcwf, mcwf_h, mcwf_nh, diagonaljumps using ...bases, ...states, ...operators, ...ode_dopri using ...operators_dense, ...operators_sparse - +using ..timeevolution +import OrdinaryDiffEq const DecayRates = Union{Vector{Float64}, Matrix{Float64}, Void} @@ -29,40 +30,40 @@ Integrate a single Monte Carlo wave function trajectory. and therefore must not be changed. * `kwargs`: Further arguments are passed on to the ode solver. """ -function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan, psi0::Ket, seed; - fout=nothing, - kwargs...) +function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan, + psi0::Ket, seed; fout=nothing, + display_beforeevent=false, display_afterevent=false, + kwargs...) tmp = copy(psi0) as_ket(x::Vector{Complex128}) = Ket(psi0.basis, x) as_vector(psi::Ket) = psi.data rng = MersenneTwister(convert(UInt, seed)) - jumpnorm = Float64[rand(rng)] - djumpnorm(t, x::Vector{Complex128}) = norm(as_ket(x))^2 - (1-jumpnorm[1]) - function dojump(t, x::Vector{Complex128}) + jumpnorm = Ref(rand(rng)) + djumpnorm(t, x::Vector{Complex128},integrator) = norm(as_ket(x))^2 - (1-jumpnorm[]) + function dojump(integrator) + x = integrator.u + t = integrator.t jumpfun(rng, t, as_ket(x), tmp) x .= tmp.data - jumpnorm[1] = rand(rng) - return ode_dopri.jump + jumpnorm[] = rand(rng) end + cb = OrdinaryDiffEq.ContinuousCallback(djumpnorm,dojump, + save_positions = (display_beforeevent,display_afterevent)) - tout = Float64[] - xout = Ket[] - function fout_(t, x::Vector{Complex128}) + function fout_(t, x::Ket) if fout==nothing - psi = copy(as_ket(x)) + psi = copy(x) psi /= norm(psi) - push!(tout, t) - push!(xout, psi) - return nothing + return psi else - return fout(t, as_ket(x)) + return fout(t, x) end end - dmcwf_(t, x::Vector{Complex128}, dx::Vector{Complex128}) = dmcwf(t, as_ket(x), as_ket(dx)) - ode_event(dmcwf_, float(tspan), as_vector(psi0), fout_, - djumpnorm, dojump; - kwargs...) - return fout==nothing ? (tout, xout) : nothing + + timeevolution.integrate(float(tspan), dmcwf, as_vector(psi0), + copy(psi0), copy(psi0), fout_; + callback = cb, + kwargs...) end """ diff --git a/src/timeevolution_base.jl b/src/timeevolution_base.jl index e3a25eba..1cb34041 100644 --- a/src/timeevolution_base.jl +++ b/src/timeevolution_base.jl @@ -8,10 +8,10 @@ function recast! end df(t, state::T, dstate::T) """ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128}, - state::T, dstate::T, fout::Function, - alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.DP5(); + state::T, dstate::T, fout::Function; + alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.DP5(), steady_state = false, eps = 1e-3, save_everystep = false, - kwargs...) + callback = nothing, kwargs...) function df_(t, x::Vector{Complex128}, dx::Vector{Complex128}) recast!(x, state) @@ -51,6 +51,8 @@ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex12 cb = scb end + full_cb = OrdinaryDiffEq.CallbackSet(callback,cb) + # TODO: Expose algorithm choice sol = OrdinaryDiffEq.solve( prob, @@ -59,12 +61,12 @@ function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex12 abstol = 1.0e-8, save_everystep = false, save_start = false, save_end = false, - callback=cb, kwargs...) + callback=full_cb, kwargs...) out.t,out.saveval end function integrate{T}(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128}, - state::T, dstate::T, ::Void, args...; kwargs...) + state::T, dstate::T, ::Void; kwargs...) function fout(t::Float64, state::T) copy(state) end