finalize the upgrade

Author: ChrisRackauckas

Project.toml

@@ -23,14 +23,14 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 DataStructures = "0.17"
 DelayDiffEq = "5"
 DiffEqBase = "6.29"
-DiffEqNoiseProcess = "3.3.1"
+DiffEqNoiseProcess = "5"
 RandomNumbers = "1"
 RecursiveArrayTools = "1,2"
 Reexport = "0.2"
 StaticArrays = "0.11, 0.12"
 StochasticDiffEq = "6.19"
 UnPack = "0.1, 1.0"
-julia = "1.3"
+julia = "1.4"
 
 [extras]
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"

src/integrators/interface.jl

@@ -298,4 +298,24 @@ end
     end
   end
   return saved, savedexactly
-end
+end
+
+@inline function DiffEqNoiseProcess.setup_next_step!(integrator::SDDEIntegrator)
+  !isnothing(integrator.W) && DiffEqNoiseProcess.setup_next_step!(integrator.W,integrator.u,integrator.p)
+  !isnothing(integrator.P) && DiffEqNoiseProcess.setup_next_step!(integrator.P,integrator.u,integrator.p)
+end
+
+@inline function DiffEqNoiseProcess.reject_step!(integrator::SDDEIntegrator,dtnew = integrator.dtnew)
+  !isnothing(integrator.W) && reject_step!(integrator.W,dtnew,integrator.u,integrator.p)
+  !isnothing(integrator.P) && reject_step!(integrator.P,dtnew,integrator.u,integrator.p)
+end
+
+@inline function DiffEqNoiseProcess.accept_step!(integrator::SDDEIntegrator,setup)
+  !isnothing(integrator.W) && accept_step!(integrator.W,integrator.dt,integrator.u,integrator.p,setup)
+  !isnothing(integrator.P) && accept_step!(integrator.P,integrator.dt,integrator.u,integrator.p,setup)
+end
+
+@inline function DiffEqNoiseProcess.save_noise!(integrator::SDDEIntegrator)
+  !isnothing(integrator.W) && DiffEqNoiseProcess.save_noise!(integrator.W)
+  !isnothing(integrator.P) && DiffEqNoiseProcess.save_noise!(integrator.P)
+end

src/integrators/type.jl

@@ -13,9 +13,10 @@ mutable struct HistorySDEIntegrator{algType,IIP,uType,tType,SolType,CacheType} <
   end
 
 mutable struct
-    SDDEIntegrator{algType,IIP,uType,uEltype,tType,P,eigenType,tTypeNoUnits,uEltypeNoUnits,randType,rateType,solType,cacheType,F,G,OType,noiseType,EventErrorType,CallbackCacheType,H,IType} <: AbstractSDDEIntegrator{algType,IIP,uType,tType}
+    SDDEIntegrator{algType,IIP,uType,uEltype,tType,P,eigenType,tTypeNoUnits,uEltypeNoUnits,randType,randType2,rateType,solType,cacheType,F,G,F6,OType,noiseType,EventErrorType,CallbackCacheType,H,IType} <: AbstractSDDEIntegrator{algType,IIP,uType,tType}
     f::F
     g::G
+    c::F6
     noise::noiseType
     uprev::uType
     tprev::tType
@@ -53,6 +54,7 @@ mutable struct
     callback_cache::CallbackCacheType
     sqdt::tType
     W::randType
+    P::randType2
     opts::OType
     iter::Int
     success_iter::Int
@@ -64,11 +66,11 @@ mutable struct
     history::H
     integrator::IType # history integrator
   end
-  
+
   function (integrator::SDDEIntegrator)(t, deriv::Type=Val{0}; idxs=nothing)
       StochasticDiffEq.current_interpolant(t, integrator, idxs, deriv)
   end
-  
+
   function (integrator::SDDEIntegrator)(val::AbstractArray, t::Union{Number,AbstractArray},deriv::Type=Val{0}; idxs=nothing)
       StochasticDiffEq.current_interpolant!(val, t, integrator, idxs, deriv)
-  end
+  end

src/solve.jl

@@ -76,9 +76,9 @@ function DiffEqBase.__init(prob::AbstractSDDEProblem,# TODO DiffEqBasee.Abstract
         error("Bridge function must be given for adaptivity. Either declare this function in noise process or set adaptive=false")
     end
 
-    if !StochasticDiffEq.alg_compatible(prob, getalg(alg))
-        error("The algorithm is not compatible with the chosen noise type. Please see the documentation on the solver methods")
-    end
+    #if !StochasticDiffEq.alg_compatible(prob, getalg(alg))
+    #    error("The algorithm is not compatible with the chosen noise type. Please see the documentation on the solver methods")
+    #end
 
     if haskey(kwargs, :initial_order)
         @warn "initial_order has been deprecated. Please specify order_discontinuity_t0 in the DDEProblem instead."
@@ -171,6 +171,41 @@ function DiffEqBase.__init(prob::AbstractSDDEProblem,# TODO DiffEqBasee.Abstract
         end
     end
 
+    #= TODO: Jump handling
+    if typeof(_prob) <: JumpProblem && _prob.regular_jump !== nothing
+
+      if !isnothing(_prob.regular_jump.mark_dist) == nothing # https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/250
+        error("Mark distributions are currently not supported in SimpleTauLeaping")
+      end
+
+      jump_prototype = zeros(_prob.regular_jump.numjumps)
+      c = _prob.regular_jump.c
+
+      if isinplace(_prob.regular_jump)
+        rate_constants = zeros(_prob.regular_jump.numjumps)
+        _prob.regular_jump.rate(rate_constants,u./u,prob.p,tspan[1])
+        P = CompoundPoissonProcess!(_prob.regular_jump.rate,t,jump_prototype,
+                                    computerates = !alg_control_rate(alg) || !adaptive,
+                                    save_everystep=save_noise,
+                                    rng = Xorshifts.Xoroshiro128Plus(_seed))
+        alg_control_rate(alg) && adaptive && P.cache.rate(P.cache.currate,u,p,tspan[1])
+      else
+        rate_constants = _prob.regular_jump.rate(u./u,prob.p,tspan[1])
+        P = CompoundPoissonProcess(_prob.regular_jump.rate,t,jump_prototype,
+                                   save_everystep=save_noise,
+                                   computerates = !alg_control_rate(alg) || !adaptive,
+                                   rng = Xorshifts.Xoroshiro128Plus(_seed))
+        alg_control_rate(alg) && adaptive && (P.cache.currate = P.cache.rate(u,p,tspan[1]))
+      end
+
+    else
+    =#
+      jump_prototype = nothing
+      c = nothing
+      P = nothing
+      rate_constants = nothing
+    #end
+
     # tstops_internal, saveat_internal, d_discontinuities_internal =
     #   tstop_saveat_disc_handling(tstops, saveat, d_discontinuities, tspan) # TODO add delays to discontinuities
     # retrieve time stops, time points at which solutions is saved, and discontinuities
@@ -284,14 +319,14 @@ function DiffEqBase.__init(prob::AbstractSDDEProblem,# TODO DiffEqBasee.Abstract
 
     # create a history function
     history = build_history_function(prob, alg, reltol_internal,
-                                    rate_prototype, noise_rate_prototype, W, _seed, dense;
+                                    rate_prototype, noise_rate_prototype, jump_prototype, W, _seed, dense;
                                     dt = dt, adaptive = adaptive,
                                     internalnorm = internalnorm)
     f_with_history, g_with_history = wrap_functions_and_history(f, g, history)
 
     sde_integrator = history.integrator;
 
-    cache = StochasticDiffEq.alg_cache(getalg(alg), prob, u, W.dW, W.dZ, p, rate_prototype, noise_rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, uprev, f_with_history, t, dt, Val{isinplace(prob)})
+    cache = StochasticDiffEq.alg_cache(getalg(alg), prob, u, W.dW, W.dZ, p, rate_prototype, noise_rate_prototype, jump_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, uprev, f_with_history, t, dt, Val{isinplace(prob)})
 
     # id = StochasticDiffEq.LinearInterpolationData(timeseries,ts)
     id = StochasticDiffEq.LinearInterpolationData(sde_integrator.sol.u, sde_integrator.sol.t)
@@ -364,12 +399,19 @@ function DiffEqBase.__init(prob::AbstractSDDEProblem,# TODO DiffEqBasee.Abstract
     success_iter = 0
     q = tTypeNoUnits(1)
 
-    integrator = SDDEIntegrator{typeof(getalg(alg)),isinplace(prob),uType,uBottomEltype,tType,typeof(p),typeof(eigen_est),QT,uEltypeNoUnits,typeof(W),rateType,typeof(sol),typeof(cache),FType,GType,typeof(opts),typeof(noise),typeof(last_event_error),typeof(callback_cache),typeof(history),typeof(sde_integrator)}(f_with_history, g_with_history, noise, uprev, tprev,
+    integrator = SDDEIntegrator{typeof(getalg(alg)),isinplace(prob),uType,
+                                uBottomEltype,tType,typeof(p),typeof(eigen_est),
+                                QT,uEltypeNoUnits,typeof(W),typeof(P),
+                                rateType,typeof(sol),typeof(cache),FType,GType,typeof(c),
+                                typeof(opts),typeof(noise),typeof(last_event_error),
+                                typeof(callback_cache),typeof(history),
+                                typeof(sde_integrator)}(f_with_history,
+                                g_with_history, c, noise, uprev, tprev,
                       order_discontinuity_t0, tracked_discontinuities,
                       t, u, p, tType(dt), tType(dt), tType(dt), dtcache, tspan[2], tdir,
                       just_hit_tstop, isout, event_last_time, vector_event_last_time, last_event_error, accept_step,
                       last_stepfail, force_stepfail, dtchangeable, u_modified, saveiter, getalg(alg), sol,
-                      cache, callback_cache, tType(dt), W,
+                      cache, callback_cache, tType(dt), W, P,
                       opts, iter, success_iter, eigen_est, EEst, q, QT(qoldinit), q11, history, sde_integrator)
 
     if initialize_integrator
@@ -387,7 +429,7 @@ function DiffEqBase.__init(prob::AbstractSDDEProblem,# TODO DiffEqBasee.Abstract
     integrator.sqdt = integrator.tdir * sqrt(abs(integrator.dt))
 
     integrator.W.dt = integrator.dt
-    DiffEqNoiseProcess.setup_next_step!(integrator.W)
+    DiffEqNoiseProcess.setup_next_step!(integrator)
 
     integrator
 end

src/utils.jl

@@ -107,8 +107,8 @@ depends on the parameter settings of the numerical solver.
 """
 function Base.sizehint!(sol::RODESolution, alg, tspan, tstops, saveat;
                         save_everystep = isempty(saveat),
-                        adaptive = StochasticDiffEq.isadaptive(getalg(alg)), 
-                        internalnorm = DiffEqBase.ODE_DEFAULT_NORM, 
+                        adaptive = StochasticDiffEq.isadaptive(getalg(alg)),
+                        internalnorm = DiffEqBase.ODE_DEFAULT_NORM,
                         dt = zero(eltype(tspan)))
   # obtain integration time
   t0 = first(tspan)
@@ -134,7 +134,7 @@ function Base.sizehint!(sol::RODESolution, alg, tspan, tstops, saveat;
   nothing
 end
 
-function build_history_function(prob, alg, reltol, rate_prototype, noise_rate_prototype, W, _seed, dense;
+function build_history_function(prob, alg, reltol, rate_prototype, noise_rate_prototype, jump_prototype, W, _seed, dense;
   dt = zero(eltype(prob.tspan)),
   adaptive = StochasticDiffEq.isadaptive(getalg(alg)),
   calck = false,
@@ -172,7 +172,7 @@ function build_history_function(prob, alg, reltol, rate_prototype, noise_rate_pr
                                       save_start = true)
 
   # # obtain cache (we alias uprev2 and uprev)
-  sde_cache = StochasticDiffEq.alg_cache(getalg(alg),prob, sde_u, W.dW, W.dZ, p, rate_prototype,noise_rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, sde_uprev, f, t0, dt, Val{isinplace(prob)})
+  sde_cache = StochasticDiffEq.alg_cache(getalg(alg),prob, sde_u, W.dW, W.dZ, p, rate_prototype, noise_rate_prototype, jump_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, sde_uprev, f, t0, dt, Val{isinplace(prob)})
 
   # build dense interpolation of history
   id = StochasticDiffEq.LinearInterpolationData(sde_timeseries,sde_ts)
@@ -253,4 +253,4 @@ function DiffEqBase.unwrap_cache(integrator::SDDEIntegrator, is_stiff)
   else
     return cache.caches[integrator.cache.current]
   end
-end
+end