Change internals to Diffeq

.travis.yml

@@ -3,7 +3,6 @@ os:
     - osx
     - linux
 julia:
-    - 0.5
     - 0.6
     - nightly
 matrix:

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
 

REQUIRE

@@ -1,3 +1,5 @@
-julia 0.5
+julia 0.6
 Compat 0.17
 Roots
+OrdinaryDiffEq 2.34.0
+DiffEqCallbacks 0.7.1

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"

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
 
 """

src/steadystate.jl

@@ -1,10 +1,7 @@
 module steadystate
 
 using ..states, ..operators, ..operators_dense, ..superoperators
-using ..timeevolution, ..metrics
-
-
-type ConvergenceReached <: Exception end
+using ..timeevolution
 
 
 """
@@ -31,38 +28,19 @@ 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...)
-    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
-    return rho0
+    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,
+                        steady_state = true,
+                        eps = eps, kwargs...)
 end
 
 """

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))

src/timeevolution_base.jl

@@ -1,33 +1,89 @@
-using ..ode_dopri
+using ..ode_dopri, ..metrics
+
+import OrdinaryDiffEq, DiffEqCallbacks
 
 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;
+            alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.DP5(),
+            steady_state = false, eps = 1e-3, save_everystep = false,
+            callback = nothing, kwargs...)
+
     function df_(t, x::Vector{Complex128}, dx::Vector{Complex128})
         recast!(x, state)
         recast!(dx, dstate)
         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_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)
+
+    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),
+                                affect!;
+                                save_positions = (false,false))
+        cb = OrdinaryDiffEq.CallbackSet(_cb,scb)
+    else
+        cb = scb
+    end
+
+    full_cb = OrdinaryDiffEq.CallbackSet(callback,cb)
+
+    # TODO: Expose algorithm choice
+    sol = OrdinaryDiffEq.solve(
+                prob,
+                alg;
+                reltol = 1.0e-6,
+                abstol = 1.0e-8,
+                save_everystep = false, save_start = false,
+                save_end = false,
+                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; 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
+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
+
+Base.@pure pure_inference(fout,T) = Core.Inference.return_type(fout, T)

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)
@@ -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

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