sciml / jumpprocesses.jl Goto Github PK

Build and simulate jump equations like Gillespie simulations and jump diffusions with constant and state-dependent rates and mix with differential equations and scientific machine learning (SciML)

Home Page: https://docs.sciml.ai/JumpProcesses/stable/

License: Other

Julia 100.00%

gillespie ssa jump-diffusion stochastic differential-equations hybrid-differential-equation ode sde stochastic-jump-equations neural-ode

jumpprocesses.jl's People

Contributors

Stargazers

Watchers

jumpprocesses.jl's Issues

Tau-Leaping and beyond

We should add more discrete aggeregators! Explicit (Euler) Tau-leaping is the obvious next step, with implicit Euler leaping and Trapezoid as the next steps (from one will come the others).

To implement these, one just needs to implement a new DiscreteAggregator. A new type needs to be declared:

https://github.com/JuliaDiffEq/DiffEqJump.jl/blob/master/src/aggregators/aggregators.jl

And the logic for building the callback needs to be implemented. For example, here's the direct:

https://github.com/JuliaDiffEq/DiffEqJump.jl/blob/master/src/aggregators/direct.jl

The DiscreteCallback for tau-leaping should work similarly. It should store a tstop for the next timepoint where it will be called. For tau-leaping, this is just the $$\Delta t$$. This could then keep the

@inline function (p::DirectJumpAggregation)(t,u,integrator) # condition
  p.next_jump==t
end

line. But the affect! function will have to change. But that change would be to add some Poissons to the values dependent on the rates. The constructor will be very similar to the one for direct.

Without adapting the time to the next jump, a first implementation should be very short.

Simulating a large number of independent jumps

We want to simulate a large number of independent jumps in a SDE problem. Code is below. When N = 300, it's reasonable, however for N = 1000 (for some cases we might need N = 10000) it's stalling. As discussed in SciML/DiffEqDocs.jl#258 (comment)
SSAStepper() is no use here, but even if wrongly someone use it, with N = 1000 it's stalling.

This type of model we are using is for continuous-time models of an economy with heterogeneous firms.

Here is a sample code:

using DifferentialEquations ,BenchmarkTools,  StochasticDiffEq, DiffEqJump, Plots

struct AffectIndex{F1, F2}
  affect_index!::F1
  index::F2
end
function (p::AffectIndex)(integrator)
    return p.affect_index!(integrator, p.index)
end

####
function μ_SDE(du,u,p,t)
  du .= p.μ
end

function σ_SDE(du,u,p,t)
  du .= p.σ
end

function test_SRIW1()
    p = (μ = 0.01, σ = 0.1, N = 500) # if all constant
    T = 10.0  # maximum time length
    x_iv = rand(p.N)  # just draws from the inital condition

    prob = SDEProblem(μ_SDE, σ_SDE, x_iv ,(0.0, T), p)
    rate(u,p,t) = 0.2
    affect_index!(integrator, index) = (integrator.u[index] = max(integrator.u[index], integrator.u[rand(1:integrator.p.N)]) )
    jumps = [ConstantRateJump(rate,AffectIndex(affect_index!, i)) for i in 1:p.N]
    jump_prob = JumpProblem(prob,Direct(),JumpSet((), jumps, nothing, nothing))
    @btime solve($jump_prob, $SRIW1());
end

List of SSAs to consider implementing

A nice review that summarizes many SSA methods is in [1] Ch. 3.

Some SSAs they we might consider implementing (there are others that could be added that they don't examine):

Many of these share common features, like dependency graphs or jump rate data structures, so we should be able to reuse a number of subcomponents in the methods.

[1] reports large network benchmarks for which RSSA/RSSA-CR and the PDM/PDM-CR give the best performance. (The RSSA methods are by the authors of [1].) [2] compares open source simulators and generally finds that the simulator based on PDM/PDM-CR gives the best performance as the number of molecules and network complexity increases, with BioNetgen coming in second (the latter uses SDM -- but it's not clear exactly why they get such good performance).

[1] Marchetti et al. Simulation Algorithms for Computational Systems Biology, see ch 3.
[2] Gupta and Mendes, An Overview of Network-Based and -Free Approaches for Stochastic Simulation of Biochemical Systems, Computation, 6, 9, 2018.

Write a Poisson rng

Distributions.jl is a pretty large dependency. We are only using it for one line of code:

https://github.com/JuliaDiffEq/DiffEqJump.jl/blob/480d6ab5106f3332ac676a4f8dce210848661880/src/simple_regular_solve.jl#L38

For Poisson random numbers. Other tests show that not using the Rmath method would actually be faster anyways:

JuliaStats/Distributions.jl#253

This blog suggests that combining two algorithms could be fastest:

https://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values/

Those two are already created as samplers in Distributions.jl

https://github.com/JuliaStats/Distributions.jl/blob/master/src/samplers/poisson.jl

This is all the code that we need to write a fast bindep-free Poisson RNG

Update `executerx!` and `evalrxrate` to support different type than stoichiometry representations

Right now MassActionJump doesn't work in hybrid models involving ODEs as it requires the species population vector to have the same type as the species index and stoichiometry coefficient. executrerx! and evalrxrate need to be updated to allow a different type for the species population.

Package compatibility caps

Ref: https://discourse.julialang.org/t/package-compatibility-caps/15301

Solving JumpProblem with EM() produces strange discontinuities after jump

The following kinds of output is produced when solving JumpProblems with the EM() solver.

This is clearly wrong as the solution jumps again after the jump of sets u to 0. I did check that this problem went away when the recommended SRIW1 solver was used instead (output is below).

The minimal working example that generated the plots is the following

using DifferentialEquations

rate(u,p,t) = 5
affect!(integrator) = (integrator.u[1] = 0.0)
jump = VariableRateJump(rate,affect!)

function f(du,u,p,t)
  du[1] = 0.
end
function g(du,u,p,t)
  du[1] = 0.05
end
prob = SDEProblem(f,g,[0.0],(0.0,1.0))

jump_prob = JumpProblem(prob,Direct(),jump)
sol_EM = solve(jump_prob, EM(), dt=0.5e-3, adaptive=false, seed=1);
sol_SRI = solve(jump_prob, SRIW1(), dt=0.5e-3, adaptive=false);

using Plots
plot(sol_EM, linewidth=1)
plot(sol_SRI, linewidth=1)

variable rate tests failing on latest release and master

Looks like the variable rate tests are failing on the current release and master:

Test Summary:       | Pass  Total
Constant Rate Tests |    2      2
  8.976325 seconds (26.17 M allocations: 1.428 GiB, 5.33% gc time)
Variable Rate Tests: Error During Test at /Users/isaacsas/.julia/packages/DiffEqJump/5ka8T/test/runtests.jl:5
  Got exception outside of a @test
  LoadError: MethodError: no method matching __init(::SDEProblem{ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(DiffEqJump, Symbol("#jump_f#51")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing},JumpSet{Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Tuple{},Nothing,Nothing}},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},Nothing,Nothing}, ::SRIW1, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}; callback=CallbackSet{Tuple{ContinuousCallback{getfield(DiffEqJump, Symbol("##62#64")),getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},ContinuousCallback{getfield(DiffEqJump, Symbol("##66#68")),getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}}((ContinuousCallback{getfield(DiffEqJump, Symbol("##62#64")),getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}(getfield(DiffEqJump, Symbol("##62#64"))(Core.Box(1)), getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}}(Core.Box(1), VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}(getfield(Main, Symbol("##11#12"))(), getfield(Main, Symbol("##13#14"))(), nothing, true, 1000, (true, true), 1.0e-12, 0)), getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}}(Core.Box(1), VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}(getfield(Main, Symbol("##11#12"))(), getfield(Main, Symbol("##13#14"))(), nothing, true, 1000, (true, true), 1.0e-12, 0)), DiffEqBase.INITIALIZE_DEFAULT, nothing, true, 1000, Bool[true, true], 1.0e-12, 0), ContinuousCallback{getfield(DiffEqJump, Symbol("##66#68")),getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}(getfield(DiffEqJump, Symbol("##66#68"))(Core.Box(2)), getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}}(Core.Box(2), VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}(getfield(Main, Symbol("##11#12"))(), getfield(Main, Symbol("##13#14"))(), nothing, true, 1000, (true, true), 1.0e-12, 0)), getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}}(Core.Box(2), VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}(getfield(Main, Symbol("##11#12"))(), getfield(Main, Symbol("##13#14"))(), nothing, true, 1000, (true, true), 1.0e-12, 0)), DiffEqBase.INITIALIZE_DEFAULT, nothing, true, 1000, Bool[true, true], 1.0e-12, 0)), ()))
  Closest candidates are:
    __init(!Matched::DiffEqBase.AbstractJumpProblem{P,J} where J, ::DiffEqBase.DEAlgorithm, ::Any, ::Any, ::Any, ::Type{Val{recompile_flag}}; callback, seed, kwargs...) where {P, recompile_flag} at /Users/isaacsas/.julia/packages/DiffEqJump/5ka8T/src/solve.jl:17
    __init(!Matched::DiffEqBase.AbstractJumpProblem{P,J} where J, ::DiffEqBase.DEAlgorithm, ::Any, ::Any, ::Any) where P at /Users/isaacsas/.julia/packages/DiffEqJump/5ka8T/src/solve.jl:17 got unsupported keyword argument "callback"
    __init(!Matched::DiffEqBase.AbstractODEProblem, !Matched::algType<:OrdinaryDiffEqAlgorithm, ::Any, ::Any, ::Any, ::Type{Val{recompile_flag}}; saveat, tstops, d_discontinuities, save_idxs, save_everystep, save_timeseries, save_on, save_start, save_end, callback, dense, calck, dt, adaptive, gamma, abstol, reltol, qmax, qmin, qsteady_min, qsteady_max, qoldinit, fullnormalize, failfactor, beta2, beta1, maxiters, dtmax, dtmin, internalnorm, internalopnorm, isoutofdomain, unstable_check, verbose, force_dtmin, timeseries_errors, dense_errors, advance_to_tstop, stop_at_next_tstop, initialize_save, progress, progress_steps, progress_name, progress_message, userdata, allow_extrapolation, initialize_integrator, alias_u0, kwargs...) where {algType<:OrdinaryDiffEqAlgorithm, recompile_flag} at /Users/isaacsas/.julia/packages/OrdinaryDiffEq/gGw8L/src/solve.jl:62
    ...
  Stacktrace:
   [1] #init#426(::Base.Iterators.Pairs{Symbol,CallbackSet{Tuple{ContinuousCallback{getfield(DiffEqJump, Symbol("##62#64")),getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},ContinuousCallback{getfield(DiffEqJump, Symbol("##66#68")),getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},Tuple{Symbol},NamedTuple{(:callback,),Tuple{CallbackSet{Tuple{ContinuousCallback{getfield(DiffEqJump, Symbol("##62#64")),getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},ContinuousCallback{getfield(DiffEqJump, Symbol("##66#68")),getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}}}}}, ::Function, ::SDEProblem{ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(DiffEqJump, Symbol("#jump_f#51")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing},JumpSet{Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Tuple{},Nothing,Nothing}},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},Nothing,Nothing}, ::SRIW1, ::Vararg{Any,N} where N) at /Users/isaacsas/.julia/packages/DiffEqBase/xzHPw/src/solve.jl:20
   [2] (::getfield(DiffEqBase, Symbol("#kw##init")))(::NamedTuple{(:callback,),Tuple{CallbackSet{Tuple{ContinuousCallback{getfield(DiffEqJump, Symbol("##62#64")),getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},ContinuousCallback{getfield(DiffEqJump, Symbol("##66#68")),getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}}}}, ::typeof(init), ::SDEProblem{ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(DiffEqJump, Symbol("#jump_f#51")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing},JumpSet{Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Tuple{},Nothing,Nothing}},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},Nothing,Nothing}, ::SRIW1, ::Array{Any,1}, ::Vararg{Any,N} where N) at ./none:0
   [3] #__init#71(::Nothing, ::UInt64, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::JumpProblem{SDEProblem{ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(DiffEqJump, Symbol("#jump_f#51")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing},JumpSet{Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Tuple{},Nothing,Nothing}},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},Nothing,Nothing},Direct,CallbackSet{Tuple{ContinuousCallback{getfield(DiffEqJump, Symbol("##62#64")),getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},ContinuousCallback{getfield(DiffEqJump, Symbol("##66#68")),getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},Nothing,Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Nothing,Nothing}, ::SRIW1, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/isaacsas/.julia/packages/DiffEqJump/5ka8T/src/solve.jl:20
   [4] __init(::JumpProblem{SDEProblem{ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(DiffEqJump, Symbol("#jump_f#51")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing},JumpSet{Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Tuple{},Nothing,Nothing}},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},Nothing,Nothing},Direct,CallbackSet{Tuple{ContinuousCallback{getfield(DiffEqJump, Symbol("##62#64")),getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},ContinuousCallback{getfield(DiffEqJump, Symbol("##66#68")),getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},Nothing,Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Nothing,Nothing}, ::SRIW1, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/isaacsas/.julia/packages/DiffEqJump/5ka8T/src/solve.jl:17
   [5] #init#426(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::JumpProblem{SDEProblem{ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(DiffEqJump, Symbol("#jump_f#51")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing},JumpSet{Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Tuple{},Nothing,Nothing}},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},Nothing,Nothing},Direct,CallbackSet{Tuple{ContinuousCallback{getfield(DiffEqJump, Symbol("##62#64")),getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},ContinuousCallback{getfield(DiffEqJump, Symbol("##66#68")),getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},Nothing,Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Nothing,Nothing}, ::SRIW1, ::Vararg{Any,N} where N) at /Users/isaacsas/.julia/packages/DiffEqBase/xzHPw/src/solve.jl:20
   [6] init(::JumpProblem{SDEProblem{ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(DiffEqJump, Symbol("#jump_f#51")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing},JumpSet{Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Tuple{},Nothing,Nothing}},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},Nothing,Nothing},Direct,CallbackSet{Tuple{ContinuousCallback{getfield(DiffEqJump, Symbol("##62#64")),getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},ContinuousCallback{getfield(DiffEqJump, Symbol("##66#68")),getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},Nothing,Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Nothing,Nothing}, ::SRIW1, ::Array{Any,1}, ::Vararg{Any,N} where N) at /Users/isaacsas/.julia/packages/DiffEqBase/xzHPw/src/solve.jl:8
   [7] #__solve#70(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::JumpProblem{SDEProblem{ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(DiffEqJump, Symbol("#jump_f#51")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing},JumpSet{Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Tuple{},Nothing,Nothing}},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},Nothing,Nothing},Direct,CallbackSet{Tuple{ContinuousCallback{getfield(DiffEqJump, Symbol("##62#64")),getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},ContinuousCallback{getfield(DiffEqJump, Symbol("##66#68")),getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},Nothing,Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Nothing,Nothing}, ::SRIW1, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/isaacsas/.julia/packages/DiffEqJump/5ka8T/src/solve.jl:6
   [8] __solve(::JumpProblem{SDEProblem{ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(DiffEqJump, Symbol("#jump_f#51")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing},JumpSet{Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Tuple{},Nothing,Nothing}},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},Nothing,Nothing},Direct,CallbackSet{Tuple{ContinuousCallback{getfield(DiffEqJump, Symbol("##62#64")),getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},ContinuousCallback{getfield(DiffEqJump, Symbol("##66#68")),getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},Nothing,Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Nothing,Nothing}, ::SRIW1, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/isaacsas/.julia/packages/DiffEqJump/5ka8T/src/solve.jl:6 (repeats 5 times)
   [9] #solve#427(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::JumpProblem{SDEProblem{ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(DiffEqJump, Symbol("#jump_f#51")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing},JumpSet{Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Tuple{},Nothing,Nothing}},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},Nothing,Nothing},Direct,CallbackSet{Tuple{ContinuousCallback{getfield(DiffEqJump, Symbol("##62#64")),getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},ContinuousCallback{getfield(DiffEqJump, Symbol("##66#68")),getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},Nothing,Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Nothing,Nothing}, ::SRIW1) at /Users/isaacsas/.julia/packages/DiffEqBase/xzHPw/src/solve.jl:39
   [10] solve(::JumpProblem{SDEProblem{ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(DiffEqJump, Symbol("#jump_f#51")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing},JumpSet{Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Tuple{},Nothing,Nothing}},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(DiffEqJump, Symbol("#jump_g#52")){SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,getfield(Main, Symbol("##15#16")),getfield(Main, Symbol("##21#22")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},getfield(Main, Symbol("##21#22")),Nothing,Nothing}},Nothing,Nothing},Direct,CallbackSet{Tuple{ContinuousCallback{getfield(DiffEqJump, Symbol("##62#64")),getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##63#65")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},ContinuousCallback{getfield(DiffEqJump, Symbol("##66#68")),getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},getfield(DiffEqJump, Symbol("##67#69")){VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},Nothing,Tuple{VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64},VariableRateJump{getfield(Main, Symbol("##11#12")),getfield(Main, Symbol("##13#14")),Nothing,Float64,Int64}},Nothing,Nothing}, ::SRIW1) at /Users/isaacsas/.julia/packages/DiffEqBase/xzHPw/src/solve.jl:27
   [11] top-level scope at none:0
   [12] include at ./boot.jl:326 [inlined]
   [13] include_relative(::Module, ::String) at ./loading.jl:1038
   [14] include(::Module, ::String) at ./sysimg.jl:29
   [15] include(::String) at ./client.jl:403
   [16] top-level scope at /Users/isaacsas/.julia/packages/DiffEqJump/5ka8T/test/runtests.jl:5
   [17] top-level scope at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/Test/src/Test.jl:1083
   [18] top-level scope at /Users/isaacsas/.julia/packages/DiffEqJump/5ka8T/test/runtests.jl:5
   [19] top-level scope at util.jl:156
   [20] top-level scope at /Users/isaacsas/.julia/packages/DiffEqJump/5ka8T/test/runtests.jl:5
   [21] top-level scope at util.jl:156
   [22] include at ./boot.jl:326 [inlined]
   [23] include_relative(::Module, ::String) at ./loading.jl:1038
   [24] include(::Module, ::String) at ./sysimg.jl:29
   [25] include(::String) at ./client.jl:403
   [26] top-level scope at none:0
   [27] eval(::Module, ::Any) at ./boot.jl:328
   [28] exec_options(::Base.JLOptions) at ./client.jl:243
   [29] _start() at ./client.jl:436
  in expression starting at /Users/isaacsas/.julia/packages/DiffEqJump/5ka8T/test/variable_rate.jl:46
Test Summary:       | Pass  Error  Total
Variable Rate Tests |    5      1      6
ERROR: LoadError: Some tests did not pass: 5 passed, 0 failed, 1 errored, 0 broken.
in expression starting at /Users/isaacsas/.julia/packages/DiffEqJump/5ka8T/test/runtests.jl:3
ERROR: Package DiffEqJump errored during testing

Register

@JuliaRegistrator register()

Allow changing reaction rate in MassActionJumps

Right now rate constants in MassActionJumps are set at construction, and can't be modified if the user passes a different parameter vector.

This could be handled by having MassActionJump's take an optional extra parameter vector: rate_idxs_in_params or such. Then when the JumpProblem is constructed the corresponding parameter from prob.p could be used instead.

Such functionality would probably be useful to allow changing parameters down the line without having to create a new JumpProblem from scratch. (I'm not sure if remake would just work or need a dispatch here.)

Variable rate tests failing

See https://travis-ci.org/github/SciML/DiffEqJump.jl/builds/687850081

Seems to be using SRW1 for an SDEProblem as the line in question is:
https://github.com/SciML/DiffEqJump.jl/blob/16c483c477bd122e2161f9de834ed466fd9f5a26/test/variable_rate.jl#L46

Coupling two JumpProblems

A number of variance reduction techniques involve constructing a coupling between two jump processes. I was thinking this could be implemented with a function called construct_coupling that takes two jump problems and produces a 2*length(u) dimensional jump problem representing the coupled process. Does this seem like a reasonable approach?

Keyword argument passing problems in solving a JumpProblem

The following minimal working example shows there is some issues with passing a keyword argument into the solve function.

The example is basically comprised of code in the documentation, except that EM() which requires "dt" was used instead of SRIW1.

using DifferentialEquations

function f(du,u,p,t)
  du[1] = u[1]
end

function g(du,u,p,t)
  du[1] = 0.1
end

prob = SDEProblem(f,g,[0.2],(0.0,10.0))

rate(u,p,t) = 2
affect!(integrator) = (integrator.u[1] = integrator.u[1]/2)
jump = ConstantRateJump(rate,affect!)

jump_prob = JumpProblem(prob,Direct(),jump)

sol_SRIW1 = solve(jump_prob, SRIW1()); #This works
sol_EM = solve(jump_prob, EM(), dt=1e-4); #Does not work

The error message and stack trace is given by

Fixed timestep methods require a choice of dt or choosing the tstops

Stacktrace:
 [1] #init#82(::Float64, ::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Void, ::Bool, ::Bool, ::Void, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Rational{Int64}, ::Float64, ::Float64, ::Rational{Int64}, ::Rational{Int64}, ::Int64, ::Rational{Int64}, ::Bool, ::Rational{Int64}, ::Rational{Int64}, ::Rational{Int64}, ::Float64, ::Float64, ::Float64, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK, ::DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::Bool, ::DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE, ::String, ::Void, ::DiffEqBase.CallbackSet{Tuple{},Tuple{DiffEqBase.DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star}}}}, ::Bool, ::Bool, ::Bool, ::Bool, ::UInt64, ::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void}, ::StochasticDiffEq.EM{true}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/onoderat/.julia/v0.6/StochasticDiffEq/src/solve.jl:81
 [2] (::DiffEqBase.#kw##init)(::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void}, ::StochasticDiffEq.EM{true}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0
 [3] #init#56(::Void, ::Int64, ::Tuple{Bool,Bool}, ::Array{Any,1}, ::Function, ::DiffEqJump.JumpProblem{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.Direct,DiffEqBase.CallbackSet{Tuple{},Tuple{DiffEqBase.DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},Tuple{},Void,Void}, ::StochasticDiffEq.EM{true}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/onoderat/.julia/v0.6/DiffEqJump/src/solve.jl:19
 [4] init(::DiffEqJump.JumpProblem{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.Direct,DiffEqBase.CallbackSet{Tuple{},Tuple{DiffEqBase.DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},Tuple{},Void,Void}, ::StochasticDiffEq.EM{true}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/onoderat/.julia/v0.6/DiffEqJump/src/solve.jl:19
 [5] #solve#55(::Array{Any,1}, ::Function, ::DiffEqJump.JumpProblem{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.Direct,DiffEqBase.CallbackSet{Tuple{},Tuple{DiffEqBase.DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},Tuple{},Void,Void}, ::StochasticDiffEq.EM{true}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/onoderat/.julia/v0.6/DiffEqJump/src/solve.jl:6
 [6] solve(::DiffEqJump.JumpProblem{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.Direct,DiffEqBase.CallbackSet{Tuple{},Tuple{DiffEqBase.DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate},Tuple{#affect!},RandomNumbers.Xorshifts.Xoroshiro128Star},Tuple{},Void,Void}, ::StochasticDiffEq.EM{true}) at /home/onoderat/.julia/v0.6/DiffEqJump/src/solve.jl:6

Large SSA support

We need to work out how we want to do large SSAs. We can build hybrid algorithms over the RegularJumps, but @isaacsas mentioned that large SSAs want to be able to do re-ordering so that's not the right data structure. But our current setup with tuples is optimized for less than 16 jumps, so if we want to use ConstantRateJump for this then we need to somehow use FunctionWrappers and allow putting them in an array. However, there might be more information needed, in which case it might be useful to write SSAs on yet another jump type.

JumpArray Typing

Instead of Float64, generalize via triangular dispatch.

Problem Library

We need a bunch of Gillespie and Reaction problems in DiffEqProblemLibrary.jl in order to run tests and benchmarks with.

Feedback before putting in the mass action update

I see a few ways I can add in a mass action direct method. My current code

Adds a MassActionJump type and updates JumpProblem and all the JumpSet stuff to handle this fourth type. This new jump is passed down as a parameter to aggregate and is stored in the JumpAggregation as a new field.
I then have a new DirectMA aggregator algorithm which supports both constant jumps (stored as function wrappers internally) and mass action jumps. It works with a mix of both, or either on their own. (So it would also be useful to people who won't have mass action jumps, but will have > ~10 constant jumps.) It conforms to the new SSA interface @alanderos91 put through yesterday.

I could just put this through as is. Alternatively, I could revise the existing direct method SSA jump routines @alanderos91 updated yesterday to support mass action jumps and function wrappers (the latter through dispatch on the time to next jump function as discussed in the other thread). However, that would mean for now that users don't really have any way to use the function wrappers since the current JumpSet processing code will always result in a tuple of constant jumps being passed into aggregate.

My instinct is to just submit the DirectMA aggregator algorithm as is, and if down the line someone wants to update the JumpSet code to allow other ways to pass in constant jumps (general iterable containers?) then we can look into merging Direct and DirectMA. That would be a good update anyways since the current tuple-based code may have issues with larger numbers of constant jumps... More generally, I will probably just stick with all future SSAs that support both mass action and constant jumps using function wrappers internally to store constant jumps. This will be a bit slower for small numbers of constant jumps, but faster once there are ~10 or more...

Any thoughts?

Hybrid systems

Hi,

This is not really an issue (so you can change the tag maybe), but I have a package PDMP.jl that does similar things for jump processes without diffusion. Many exact algorithms are implemented like rejection algorithms but also a trick, called CHV, which I discovered and described in Arxiv.

CHV requires to have good stiff solvers and there were not available in JuliaDiffEq at the time of writing the package.

Best regards

Romain

bitcast: target type not a leaf primitive type when using interpret for JumpProblem

I am encountering an issue with using JumpProblems when I interpret a real vector as a complex vector. This is useful for quantum mechanical simulations as the equations are written with complex numbers and gives me access to better solvers which are often only implemented in the real case. This is achieved with the following function,

function ascomplex(x::AbstractVector{T}) where {T<:Real}
    return reinterpret(Complex{T}, x)
end

using interpret does not cause any issue when an SDEProblem is solved.

using DifferentialEquations

function f(du,u,p,t)
     du_com = ascomplex(du)
     du_com .= im
end

function g(du,u,p,t)
    du .= 0.1

end

sde_prob = SDEProblem(f,g,[0.0, 0.0],(0.0,1.0))
sde_sol = solve(sde_prob, EM(), dt=1e-3);

There is an issue when I then use this with a JumpProblem,

rate(u,p,t) = 5
affect!(integrator) = (integrator.u[1] = 0.0)
jump = VariableRateJump(rate,affect!)

jump_prob = JumpProblem(sde_prob,Direct(),jump)
sol_jump = solve(jump_prob, EM(), dt=1e-3, adaptive=false, seed=1);

The error message and stack trace is given by,

bitcast: target type not a leaf primitive type

Stacktrace:
 [1] f at ./In[75]:2 [inlined]
 [2] (::DiffEqJump.#jump_f#37{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.JumpSet{Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Tuple{},Void,Void}})(::Array{Float64,1}, ::DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}}, ::Void, ::Float64) at /Users/tatsuhiroonodera/.julia/v0.6/DiffEqJump/src/problem.jl:67
 [3] perform_step!(::StochasticDiffEq.SDEIntegrator{StochasticDiffEq.EM{true},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Float64,Float64,Void,Float64,Float64,Float64,DiffEqNoiseProcess.NoiseProcess{Float64,2,Float64,Array{Float64,1},Void,Void,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_DIST,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_BRIDGE,true,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Void}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Void}},DiffEqNoiseProcess.RSWM{:RSwM1,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},DiffEqBase.RODESolution{Float64,2,Array{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},1},Void,Void,Array{Float64,1},DiffEqNoiseProcess.NoiseProcess{Float64,2,Float64,Array{Float64,1},Void,Void,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_DIST,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_BRIDGE,true,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Void}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Void}},DiffEqNoiseProcess.RSWM{:RSwM1,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},DiffEqBase.SDEProblem{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Float64,true,Void,Void,Void,DiffEqJump.#jump_f#37{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.JumpSet{Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Tuple{},Void,Void}},#g,Void,UniformScaling{Int64},Void},StochasticDiffEq.EM{true},StochasticDiffEq.LinearInterpolationData{Array{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},1},Array{Float64,1}}},StochasticDiffEq.EMCache{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Array{Float64,1},Array{Float64,1},Array{Float64,1},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}}},Void,DiffEqJump.#jump_f#37{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.JumpSet{Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Tuple{},Void,Void}},#g,StochasticDiffEq.SDEOptions{Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{DiffEqBase.ContinuousCallback{DiffEqJump.##51#53,DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqBase.#INITIALIZE_DEFAULT,Float64,Int64,Void}},Tuple{}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Void,Void,Float64,Float64,Float64,Float64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void}, ::StochasticDiffEq.EMCache{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Array{Float64,1},Array{Float64,1},Array{Float64,1},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}}}, ::Function) at /Users/tatsuhiroonodera/.julia/v0.6/StochasticDiffEq/src/perform_step/low_order.jl:25
 [4] solve!(::StochasticDiffEq.SDEIntegrator{StochasticDiffEq.EM{true},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Float64,Float64,Void,Float64,Float64,Float64,DiffEqNoiseProcess.NoiseProcess{Float64,2,Float64,Array{Float64,1},Void,Void,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_DIST,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_BRIDGE,true,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Void}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Void}},DiffEqNoiseProcess.RSWM{:RSwM1,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},DiffEqBase.RODESolution{Float64,2,Array{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},1},Void,Void,Array{Float64,1},DiffEqNoiseProcess.NoiseProcess{Float64,2,Float64,Array{Float64,1},Void,Void,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_DIST,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_BRIDGE,true,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Void}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Void}},DiffEqNoiseProcess.RSWM{:RSwM1,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},DiffEqBase.SDEProblem{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Float64,true,Void,Void,Void,DiffEqJump.#jump_f#37{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.JumpSet{Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Tuple{},Void,Void}},#g,Void,UniformScaling{Int64},Void},StochasticDiffEq.EM{true},StochasticDiffEq.LinearInterpolationData{Array{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},1},Array{Float64,1}}},StochasticDiffEq.EMCache{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Array{Float64,1},Array{Float64,1},Array{Float64,1},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}}},Void,DiffEqJump.#jump_f#37{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.JumpSet{Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Tuple{},Void,Void}},#g,StochasticDiffEq.SDEOptions{Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{DiffEqBase.ContinuousCallback{DiffEqJump.##51#53,DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqBase.#INITIALIZE_DEFAULT,Float64,Int64,Void}},Tuple{}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Void,Void,Float64,Float64,Float64,Float64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void}) at /Users/tatsuhiroonodera/.julia/v0.6/StochasticDiffEq/src/solve.jl:384
 [5] #solve#55(::Array{Any,1}, ::Function, ::DiffEqJump.JumpProblem{DiffEqBase.SDEProblem{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Float64,true,Void,Void,Void,DiffEqJump.#jump_f#37{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.JumpSet{Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Tuple{},Void,Void}},#g,Void,UniformScaling{Int64},Void},DiffEqJump.Direct,DiffEqBase.CallbackSet{Tuple{DiffEqBase.ContinuousCallback{DiffEqJump.##51#53,DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqBase.#INITIALIZE_DEFAULT,Float64,Int64,Void}},Tuple{}},Void,Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Void,Void}, ::StochasticDiffEq.EM{true}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/tatsuhiroonodera/.julia/v0.6/DiffEqJump/src/solve.jl:7
 [6] (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqJump.JumpProblem{DiffEqBase.SDEProblem{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Float64,true,Void,Void,Void,DiffEqJump.#jump_f#37{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.JumpSet{Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Tuple{},Void,Void}},#g,Void,UniformScaling{Int64},Void},DiffEqJump.Direct,DiffEqBase.CallbackSet{Tuple{DiffEqBase.ContinuousCallback{DiffEqJump.##51#53,DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqBase.#INITIALIZE_DEFAULT,Float64,Int64,Void}},Tuple{}},Void,Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Void,Void}, ::StochasticDiffEq.EM{true}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0 (repeats 2 times)

I understand this may be an esoteric application that DifferentialEquation was not designed for. If you could give me some pointers on roughly how to debug this issue - I can also try to solve this problem myself.

Re-enable FunctionWrappers in v0.7 after update

DirectFW(), FRMFW(), SortingDirect(), NRM() are all disabled until FunctionWrappers.jl updates for v0.7.

Fast broadcast for ExtendedJumpArray and CoupledArray does not work with AbstractArray

6197312

Could not find out how to get the views right to split them from the generated function.

Jump DDEs error

Hi,

I tried to run a Jump DDE and, given the supreme coolness of DifferentialEquations.jl, I rather expected this to just work. However, it did not.

M(not)WE:

using DifferentialEquations

function dde(du, u, h, p, t)
    du[1] = h(p, t-p[1], idxs=1)
end
h(p, t; kwargs...) = 0
p = [ 5.]
tspan = (0., 15)
u0 = [1.]

prob = DDEProblem(dde, u0, h, tspan, p)

### Add some Gillespie action
prod_rate(u,p,t) = (1 + 9 * u[1] / (5 + u[1]))/10
prod_affect!(i) = (i.u[1] += 1)
prod_jump = ConstantRateJump(prod_rate, prod_affect!)

deg_rate(u,p,t) = u[1]/10
deg_affect!(i) = (i.u[1] -= 1)
deg_jump = ConstantRateJump(deg_rate, deg_affect!)

jump_prob = JumpProblem(prob, Direct(), prod_jump, deg_jump)

sol = solve(jump_prob, MethodOfSteps(Tsit5()))

MethodError: no method matching __init(::DDEProblem{Array{Float64,1},Tuple{Float64,Float64},Array{Any,1},Array{Any,1},true,Array{Float64,1},DDEFunction{true,typeof(dde),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(h),Nothing}, ::MethodOfSteps{Tsit5,Nothing,Nothing,Nothing,false}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}; callback=CallbackSet{Tuple{},Tuple{DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}}}}((), (DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}}(DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}(0, Inf, 15.0, [6.93198e-310, 6.93186e-310], 0.0, MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}}(Float64[], Array{Pair{Int64,Float64},1}[], Array{Pair{Int64,Float64},1}[]), (prod_rate, deg_rate), (prod_affect!, deg_affect!), (true, true), RandomNumbers.Xorshifts.Xoroshiro128Star(0x388bd4d0a1d82370, 0x4a55a747b07e9e37)), DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}(0, Inf, 15.0, [6.93198e-310, 6.93186e-310], 0.0, MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}}(Float64[], Array{Pair{Int64,Float64},1}[], Array{Pair{Int64,Float64},1}[]), (prod_rate, deg_rate), (prod_affect!, deg_affect!), (true, true), RandomNumbers.Xorshifts.Xoroshiro128Star(0x388bd4d0a1d82370, 0x4a55a747b07e9e37)), DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}(0, Inf, 15.0, [6.93198e-310, 6.93186e-310], 0.0, MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}}(Float64[], Array{Pair{Int64,Float64},1}[], Array{Pair{Int64,Float64},1}[]), (prod_rate, deg_rate), (prod_affect!, deg_affect!), (true, true), RandomNumbers.Xorshifts.Xoroshiro128Star(0x388bd4d0a1d82370, 0x4a55a747b07e9e37)), Bool[true, true]),)))
Closest candidates are:
  __init(::DiffEqBase.AbstractDDEProblem{uType,tupType,lType,iip}, ::algType<:DelayDiffEq.AbstractMethodOfStepsAlgorithm, ::Any, ::Any, ::Any; d_discontinuities, dtmax, dt, saveat, tstops, save_idxs, save_everystep, save_start, save_end, save_on, dense, minimal_solution, discontinuity_interp_points, discontinuity_abstol, discontinuity_reltol, initial_order, initialize_integrator, initialize_save, callback, alias_u0, kwargs...) where {uType, tupType, lType, iip, algType<:AbstractMethodOfStepsAlgorithm} at /home/Niklas/.julia/packages/DelayDiffEq/NcMLM/src/solve.jl:22
  __init(!Matched::DiffEqBase.AbstractJumpProblem{P,J} where J, ::DiffEqBase.DEAlgorithm, ::Any, ::Any, ::Any, ::Type{Val{recompile_flag}}; callback, seed, kwargs...) where {P, recompile_flag} at /home/Niklas/.julia/packages/DiffEqJump/cDHIH/src/solve.jl:17
  __init(!Matched::DiffEqBase.AbstractODEProblem, !Matched::algType<:OrdinaryDiffEqAlgorithm, ::Any, ::Any, ::Any, ::Type{Val{recompile_flag}}; saveat, tstops, d_discontinuities, save_idxs, save_everystep, save_timeseries, save_on, save_start, save_end, callback, dense, calck, dt, adaptive, gamma, abstol, reltol, qmax, qmin, qsteady_min, qsteady_max, qoldinit, fullnormalize, failfactor, beta2, beta1, maxiters, dtmax, dtmin, internalnorm, internalopnorm, isoutofdomain, unstable_check, verbose, force_dtmin, timeseries_errors, dense_errors, advance_to_tstop, stop_at_next_tstop, initialize_save, progress, progress_steps, progress_name, progress_message, userdata, allow_extrapolation, initialize_integrator, alias_u0, kwargs...) where {algType<:OrdinaryDiffEqAlgorithm, recompile_flag} at /home/Niklas/.julia/packages/OrdinaryDiffEq/TUKTm/src/solve.jl:62
  ...

Stacktrace:
 [1] #init#455(::Base.Iterators.Pairs{Symbol,CallbackSet{Tuple{},Tuple{DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},Tuple{Symbol},NamedTuple{(:callback,),Tuple{CallbackSet{Tuple{},Tuple{DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}}}}}}}, ::Function, ::DDEProblem{Array{Float64,1},Tuple{Float64,Float64},Array{Any,1},Array{Any,1},true,Array{Float64,1},DDEFunction{true,typeof(dde),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(h),Nothing}, ::MethodOfSteps{Tsit5,Nothing,Nothing,Nothing,false}, ::Vararg{Any,N} where N) at /home/Niklas/.julia/packages/DiffEqBase/6ewdP/src/solve.jl:20
 [2] (::getfield(DiffEqBase, Symbol("#kw##init")))(::NamedTuple{(:callback,),Tuple{CallbackSet{Tuple{},Tuple{DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}}}}}}, ::typeof(init), ::DDEProblem{Array{Float64,1},Tuple{Float64,Float64},Array{Any,1},Array{Any,1},true,Array{Float64,1},DDEFunction{true,typeof(dde),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(h),Nothing}, ::MethodOfSteps{Tsit5,Nothing,Nothing,Nothing,false}, ::Array{Any,1}, ::Vararg{Any,N} where N) at ./none:0
 [3] #__init#75(::Nothing, ::UInt64, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::JumpProblem{DDEProblem{Array{Float64,1},Tuple{Float64,Float64},Array{Any,1},Array{Any,1},true,Array{Float64,1},DDEFunction{true,typeof(dde),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(h),Nothing},Direct,CallbackSet{Tuple{},Tuple{DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},Tuple{},Nothing,Nothing}, ::MethodOfSteps{Tsit5,Nothing,Nothing,Nothing,false}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/Niklas/.julia/packages/DiffEqJump/cDHIH/src/solve.jl:20
 [4] __init(::JumpProblem{DDEProblem{Array{Float64,1},Tuple{Float64,Float64},Array{Any,1},Array{Any,1},true,Array{Float64,1},DDEFunction{true,typeof(dde),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(h),Nothing},Direct,CallbackSet{Tuple{},Tuple{DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},Tuple{},Nothing,Nothing}, ::MethodOfSteps{Tsit5,Nothing,Nothing,Nothing,false}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/Niklas/.julia/packages/DiffEqJump/cDHIH/src/solve.jl:17
 [5] #init#455(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::JumpProblem{DDEProblem{Array{Float64,1},Tuple{Float64,Float64},Array{Any,1},Array{Any,1},true,Array{Float64,1},DDEFunction{true,typeof(dde),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(h),Nothing},Direct,CallbackSet{Tuple{},Tuple{DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},Tuple{},Nothing,Nothing}, ::MethodOfSteps{Tsit5,Nothing,Nothing,Nothing,false}, ::Vararg{Any,N} where N) at /home/Niklas/.julia/packages/DiffEqBase/6ewdP/src/solve.jl:20
 [6] init(::JumpProblem{DDEProblem{Array{Float64,1},Tuple{Float64,Float64},Array{Any,1},Array{Any,1},true,Array{Float64,1},DDEFunction{true,typeof(dde),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(h),Nothing},Direct,CallbackSet{Tuple{},Tuple{DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},Tuple{},Nothing,Nothing}, ::MethodOfSteps{Tsit5,Nothing,Nothing,Nothing,false}, ::Array{Any,1}, ::Vararg{Any,N} where N) at /home/Niklas/.julia/packages/DiffEqBase/6ewdP/src/solve.jl:8
 [7] #__solve#74(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::JumpProblem{DDEProblem{Array{Float64,1},Tuple{Float64,Float64},Array{Any,1},Array{Any,1},true,Array{Float64,1},DDEFunction{true,typeof(dde),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(h),Nothing},Direct,CallbackSet{Tuple{},Tuple{DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},Tuple{},Nothing,Nothing}, ::MethodOfSteps{Tsit5,Nothing,Nothing,Nothing,false}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/Niklas/.julia/packages/DiffEqJump/cDHIH/src/solve.jl:6
 [8] __solve(::JumpProblem{DDEProblem{Array{Float64,1},Tuple{Float64,Float64},Array{Any,1},Array{Any,1},true,Array{Float64,1},DDEFunction{true,typeof(dde),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(h),Nothing},Direct,CallbackSet{Tuple{},Tuple{DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},Tuple{},Nothing,Nothing}, ::MethodOfSteps{Tsit5,Nothing,Nothing,Nothing,false}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/Niklas/.julia/packages/DiffEqJump/cDHIH/src/solve.jl:6 (repeats 5 times)
 [9] #solve#456(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::JumpProblem{DDEProblem{Array{Float64,1},Tuple{Float64,Float64},Array{Any,1},Array{Any,1},true,Array{Float64,1},DDEFunction{true,typeof(dde),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(h),Nothing},Direct,CallbackSet{Tuple{},Tuple{DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},Tuple{},Nothing,Nothing}, ::MethodOfSteps{Tsit5,Nothing,Nothing,Nothing,false}) at /home/Niklas/.julia/packages/DiffEqBase/6ewdP/src/solve.jl:39
 [10] solve(::JumpProblem{DDEProblem{Array{Float64,1},Tuple{Float64,Float64},Array{Any,1},Array{Any,1},true,Array{Float64,1},DDEFunction{true,typeof(dde),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(h),Nothing},Direct,CallbackSet{Tuple{},Tuple{DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},DiffEqJump.DirectJumpAggregation{Float64,MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{typeof(prod_rate),typeof(deg_rate)},Tuple{typeof(prod_affect!),typeof(deg_affect!)},RandomNumbers.Xorshifts.Xoroshiro128Star},Tuple{},Nothing,Nothing}, ::MethodOfSteps{Tsit5,Nothing,Nothing,Nothing,false}) at /home/Niklas/.julia/packages/DiffEqBase/6ewdP/src/solve.jl:27
 [11] top-level scope at In[59]:23

I can achieve what I want using callbacks, so I'm not desperate for a fix, but I thought I should report the problem.

Jump Diffusion Monte Carlo simulations: the first jump occurs in the same time

discussion here
rate(u,p,t) = 2 affect!(integrator) = (integrator.u[1] = integrator.u[1]/2) jump = VariableRateJump(rate,affect!) jump_prob = JumpProblem(prob,Direct(),jump) monte_prob = MonteCarloProblem(jump_prob) sol = solve(monte_prob,SRIW1(),num_monte=3,parallel_type=:none,dt=0.001,adaptive=false)

Support for marked point processes

Following the discussion at https://discourse.julialang.org/t/differentialequations-jl-marked-point-processes-or-storing-user-data/19453, I suggest adding support for marked point processes. Specifically, for simulating the process alongside an SDE affected by it, and retrieving times and marks from the solution.

As a suggested interface (slightly different than the one I suggested in the discourse thread above), the constructors for VariableRateJump and ConstantRateJump can have an additional optional (keyword?) argument mark specifying a mark sampling function. The affect! function would then take an additional argument for the mark.

total_rate(u,p,t) = ...
sample_mark(u,p,t) = ...
affect!(du,u,p,t,y) = ...
mpp = VariableRateJump( total_rate, affect!, mark=sample_mark )
prob = JumpProblem( ode_or_sde, mpp )
sol = solve(prob)
times, marks = get_point_events(sol)

which would mean something like

function _affect!(du,u,p,t)
    y = sample_mark(u,p,t)
    store_time_and_mark_somewhere!(t,y,somewhere)
    affect!(du,u,p,t,y)
end
jump = VariableRateJump( total_rate, _affect! )
prob = JumpProblem( ode_or_sde, jump )

Mathematically, this is enough to support equations of the form
du_t=<drift and diffusion terms> +∫_y h(u_t,p,t,y)N(dt,dy)
where N is a MPP with an intensity kernel of the form λₜ(dy)=λ(u_t,p,t,dy) (total_rate is ∫λₜ(dy), sample_mark samples y from λₜ(dy)/∫λₜ(dy))

Jump PDEs

Any plans for Jump PDEs? Here are some typical problems: arXiv:1602.07398. I haven't looked too closely at the PDE code, but callbacks are not yet implemented for PDEs, right?

extend_problem for AbstractSDEProblem does not extend diffusion (g) function

The extrend_problem function extends drift (f) function to work with the ExtendedJumpArray state but it does not extend the diffusion (g) function.

MWE is below

using DifferentialEquations

rate(u,p,t) = 10
affect!(integrator) = (integrator.u .= 0.0)
jump = VariableRateJump(rate,affect!)

function f(du,u,p,t)
    du .= 1.0
end

function g(du,u,p,t)
    du .= 2*u
end

u0 = zeros(4)
sde_prob = SDEProblem(f,g,u0,(0.0,1.0))

jump_prob = JumpProblem(sde_prob,Direct(),jump)
sol_jump = solve(jump_prob, EM(), dt=1e-3, adaptive=false, seed=1);

This gives the following error

MethodError: objects of type Int64 are not callable

Stacktrace:
 [1] broadcast_c at /Users/tatsuhiroonodera/.julia/v0.6/DiffEqJump/src/extended_jump_array.jl:63 [inlined]
 [2] broadcast at ./broadcast.jl:455 [inlined]
 [3] * at ./arraymath.jl:45 [inlined]
 [4] g at ./In[9]:6 [inlined]
 [5] perform_step!(::StochasticDiffEq.SDEIntegrator{StochasticDiffEq.EM{true},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Float64,Float64,Void,Float64,Float64,Float64,DiffEqNoiseProcess.NoiseProcess{Float64,2,Float64,Array{Float64,1},Void,Void,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_DIST,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_BRIDGE,true,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Void}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Void}},DiffEqNoiseProcess.RSWM{:RSwM1,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},DiffEqBase.RODESolution{Float64,2,Array{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},1},Void,Void,Array{Float64,1},DiffEqNoiseProcess.NoiseProcess{Float64,2,Float64,Array{Float64,1},Void,Void,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_DIST,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_BRIDGE,true,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Void}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Void}},DiffEqNoiseProcess.RSWM{:RSwM1,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},DiffEqBase.SDEProblem{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Float64,true,Void,Void,Void,DiffEqJump.#jump_f#37{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.JumpSet{Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Tuple{},Void,Void}},#g,Void,UniformScaling{Int64},Void},StochasticDiffEq.EM{true},StochasticDiffEq.LinearInterpolationData{Array{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},1},Array{Float64,1}}},StochasticDiffEq.EMCache{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Array{Float64,1},Array{Float64,1},Array{Float64,1},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}}},Void,DiffEqJump.#jump_f#37{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.JumpSet{Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Tuple{},Void,Void}},#g,StochasticDiffEq.SDEOptions{Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{DiffEqBase.ContinuousCallback{DiffEqJump.##51#53,DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqBase.#INITIALIZE_DEFAULT,Float64,Int64,Void}},Tuple{}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Void,Void,Float64,Float64,Float64,Float64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void}, ::StochasticDiffEq.EMCache{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Array{Float64,1},Array{Float64,1},Array{Float64,1},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}}}, ::Function) at /Users/tatsuhiroonodera/.julia/v0.6/StochasticDiffEq/src/perform_step/low_order.jl:35
 [6] solve!(::StochasticDiffEq.SDEIntegrator{StochasticDiffEq.EM{true},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Float64,Float64,Void,Float64,Float64,Float64,DiffEqNoiseProcess.NoiseProcess{Float64,2,Float64,Array{Float64,1},Void,Void,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_DIST,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_BRIDGE,true,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Void}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Void}},DiffEqNoiseProcess.RSWM{:RSwM1,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},DiffEqBase.RODESolution{Float64,2,Array{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},1},Void,Void,Array{Float64,1},DiffEqNoiseProcess.NoiseProcess{Float64,2,Float64,Array{Float64,1},Void,Void,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_DIST,DiffEqNoiseProcess.#INPLACE_WHITE_NOISE_BRIDGE,true,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Void}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Void}},DiffEqNoiseProcess.RSWM{:RSwM1,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},DiffEqBase.SDEProblem{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Float64,true,Void,Void,Void,DiffEqJump.#jump_f#37{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.JumpSet{Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Tuple{},Void,Void}},#g,Void,UniformScaling{Int64},Void},StochasticDiffEq.EM{true},StochasticDiffEq.LinearInterpolationData{Array{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},1},Array{Float64,1}}},StochasticDiffEq.EMCache{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Array{Float64,1},Array{Float64,1},Array{Float64,1},DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}}},Void,DiffEqJump.#jump_f#37{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.JumpSet{Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Tuple{},Void,Void}},#g,StochasticDiffEq.SDEOptions{Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{DiffEqBase.ContinuousCallback{DiffEqJump.##51#53,DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqBase.#INITIALIZE_DEFAULT,Float64,Int64,Void}},Tuple{}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Void,Void,Float64,Float64,Float64,Float64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void}) at /Users/tatsuhiroonodera/.julia/v0.6/StochasticDiffEq/src/solve.jl:384
 [7] #solve#55(::Array{Any,1}, ::Function, ::DiffEqJump.JumpProblem{DiffEqBase.SDEProblem{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Float64,true,Void,Void,Void,DiffEqJump.#jump_f#37{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.JumpSet{Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Tuple{},Void,Void}},#g,Void,UniformScaling{Int64},Void},DiffEqJump.Direct,DiffEqBase.CallbackSet{Tuple{DiffEqBase.ContinuousCallback{DiffEqJump.##51#53,DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqBase.#INITIALIZE_DEFAULT,Float64,Int64,Void}},Tuple{}},Void,Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Void,Void}, ::StochasticDiffEq.EM{true}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/tatsuhiroonodera/.julia/v0.6/DiffEqJump/src/solve.jl:7
 [8] (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqJump.JumpProblem{DiffEqBase.SDEProblem{DiffEqJump.ExtendedJumpArray{Array{Float64,1},Array{Float64,1}},Float64,true,Void,Void,Void,DiffEqJump.#jump_f#37{DiffEqBase.SDEProblem{Array{Float64,1},Float64,true,Void,Void,Void,#f,#g,Void,UniformScaling{Int64},Void},DiffEqJump.JumpSet{Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Tuple{},Void,Void}},#g,Void,UniformScaling{Int64},Void},DiffEqJump.Direct,DiffEqBase.CallbackSet{Tuple{DiffEqBase.ContinuousCallback{DiffEqJump.##51#53,DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqJump.##52#54{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},DiffEqBase.#INITIALIZE_DEFAULT,Float64,Int64,Void}},Tuple{}},Void,Tuple{DiffEqJump.VariableRateJump{#rate,#affect!,Void,Float64,Int64}},Void,Void}, ::StochasticDiffEq.EM{true}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0 (repeats 2 times)

Broken interpolation when solving with `save_idxs`

Hi,
the interpolation of a jump simulation is broken if save_idxs is specified in solve. The error occurs when one tries to evaluate the simulation object at the first timepoint. It seems to be fine for the rest.

function f(du,u,p,t)
  du[1] = u[1]
end

prob = ODEProblem(f,[0.2],(0.0,10.0))
rate(u,p,t) = 2
affect!(integrator) = (integrator.u[1] = integrator.u[1]/2)
jump = ConstantRateJump(rate,affect!)
jump_prob = JumpProblem(prob,Direct(),jump)
sol = solve(jump_prob; save_idxs=1)

sol(0.) # broken
plot(sol) # broken

ERROR: MethodError: Cannot `convert` an object of type Array{Float64,1} to an object of type Float64
Closest candidates are:
  convert(::Type{T}, ::Ratios.SimpleRatio{S}) where {T<:AbstractFloat, S} at /home/Niklas/.julia/packages/Ratios/iJ67w/src/Ratios.jl:16
  convert(::Type{T}, ::T) where T<:Number at number.jl:6
  convert(::Type{T}, ::Number) where T<:Number at number.jl:7
  ...
Stacktrace:
    Function                  Module               Signature
    ────────                  ──────               ─────────
[1] push!                     Base                 (::Array{Float64,1}, ::Array{Float64,1})
at: ~/programs/julia-1.3.0/bin/../share/julia/base/array.jl:868
[2] copyat_or_push!           RecursiveArrayTools  (::Array{Float64,1}, ::Int64, ::Array{Float64,1}, ::Type{Val{true}})
at: ~/.julia/packages/RecursiveArrayTools/euiv5/src/utils.jl:68
[3] copyat_or_push!           RecursiveArrayTools  (::Array{Float64,1}, ::Int64, ::Array{Float64,1})
at: ~/.julia/packages/RecursiveArrayTools/euiv5/src/utils.jl:58
[4] addsteps!                 OrdinaryDiffEq       (::Array{Float64,1}, ::Float64, ::Float64, ::Float64, ::Float64, ::ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}, ::DiffEqBase.NullParameters, ::OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}, ::Bool, ::Bool, ::Bool)
at: ~/.julia/packages/OrdinaryDiffEq/8Pn99/src/dense/low_order_rk_addsteps.jl:139
[5] addsteps!                 OrdinaryDiffEq       (::Array{Float64,1}, ::Float64, ::Float64, ::Float64, ::Float64, ::Function, ::DiffEqBase.NullParameters, ::OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}})
at: ~/.julia/packages/OrdinaryDiffEq/8Pn99/src/dense/low_order_rk_addsteps.jl:124
[6] ode_interpolation         OrdinaryDiffEq       (::Float64, ::Function, ::Nothing, ::Type, ::DiffEqBase.NullParameters, ::Symbol)
at: ~/.julia/packages/OrdinaryDiffEq/8Pn99/src/dense/generic_dense.jl:281
[7] #_#325 [i]
at: ~/.julia/packages/OrdinaryDiffEq/8Pn99/src/interp_func.jl:74
[8] ODECompositeSolution [i]
at: ~/.julia/packages/OrdinaryDiffEq/8Pn99/src/composite_solution.jl:16
[9] top-level scope
at: REPL[36]:1

Cheers!

VariableRateJumps cause a broadcast error with certain types of solvers

Minimum reproduction:

using DifferentialEquations

function rate(u,p,t)
    return t * u[1]
end 
function prod!(integrator)
    integrator.u[1] += 1
end

constant_jump = ConstantRateJump(rate, prod!)
variable_jump = VariableRateJump(rate, prod!)
p = [1.0]
u0 = [0.0]
tspan = (0.0, 2.0)
prob = DiscreteProblem(u0, tspan, p)
constant_prob = JumpProblem(prob, Direct(), constant_jump)
sol1 = solve(constant_prob, FunctionMap())
println(sol1) #goes fine
prob2 = DiscreteProblem(u0, tspan, p)
variable_prob = JumpProblem(prob2, Direct(), constant_jump, variable_jump)
sol2 = solve(variable_prob, Rodas5()) #doesn't go fine
println(sol2)

Stack trace:

LoadError: MethodError: no method matching __broadcast!(::getfield(Base.Broadcast, Symbol("##1#2")){Base.Broadcast.Broadcasted{Base.Broadcast.ArrayStyle{ExtendedJumpArray},Tuple{Base.One
To{Int64}},typeof(muladd),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.ArrayStyle{ExtendedJumpArray},Nothing,typeof(DiffEqBase.ODE_DEFAULT_NORM),Tuple{ExtendedJumpArray{Array{Float64,1},Array{Float64,1}}}},Float64,Float64}},getfield(Ba
se.Broadcast, Symbol("##7#8")){Base.Broadcast.Broadcasted{Base.Broadcast.ArrayStyle{ExtendedJumpArray},Nothing,typeof(DiffEqBase.ODE_DEFAULT_NORM),Tuple{ExtendedJumpArray{Array{Float64,1},Array{Float64,1}}}},getfield(Base.Broadcast, Symbo
l("##9#10")){getfield(Base.Broadcast, Symbol("##11#12"))},getfield(Base.Broadcast, Symbol("##13#14")){getfield(Base.Broadcast, Symbol("##15#16"))},getfield(Base.Broadcast, Symbol("##5#6")){getfield(Base.Broadcast, Symbol("##5#6")){getfiel
d(Base.Broadcast, Symbol("##5#6")){getfield(Base.Broadcast, Symbol("##3#4"))}}}}}, ::Array{Float64,1}, ::ExtendedJumpArray{Array{Float64,1},Array{Float64,1}}, ::Float64, ::Float64)

Package versions:

  [a134a8b2] BlackBoxOptim v0.4.0                                                                                      │
  [34f1f09b] ClusterManagers v0.3.2                                                                                    │
  [2b5f629d] DiffEqBase v4.27.0+ #master (https://github.com/JuliaDiffEq/DiffEqBase.jl.git)                            │
  [0c46a032] DifferentialEquations v5.2.1                                                                              │
  [aaf54ef3] DistributedArrays v0.5.1                                                                                  │
  [28b8d3ca] GR v0.34.1                                                                                                │
  [033835bb] JLD2 v0.1.2                                                                                               │
  [1dea7af3] OrdinaryDiffEq v4.13.0                                                                                    │
  [91a5bcdd] Plots v0.19.3                                                                                             │
  [295af30f] Revise v0.7.12

Julia version 0.7.0.

Error happens on Rodas5() but not on Tsit5().

Accessing an ArrayPartition through x doesn't work with VariableRateJump

Example:

function f(du,u,p,t)
  #du[1] = -u[1]
  du.x[1][1] = -u.x[1][1]  # ERROR here
end

u0 = ArrayPartition([1.0],[1.0])
ode = ODEProblem(f,u0,(0.0,1.0))

rate(u,p,t) = 10*u[1]
function affect!(integrator)
  #integrator.u[2] += randn()
  integrator.u.x[2] .+= randn() # ERROR here
end
jump = VariableRateJump(rate,affect!)
jump_prob = JumpProblem(ode,Direct(),jump)
sol = solve(jump_prob,Tsit5())

Results in:
type SubArray has no field x
inside the function f, where apparently u is of type SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true} rather than an ArrayPartition.
Replacing the offending line with the commented-out line in f, a similar exception appears in affect!:
type ExtendedJumpArray has no field x.
In the offending line is similarly replaced in affect!, it works as expected.

Proposal: Abstraction for joint distribution sampling

Disclaimer: I cannot speak to whether this approach is valid/generalizes to more general jumps. The discussion that follows is predicated on pure Markov jump processes (and therefore applies to ConstantRateJump).

Gillespie's paper identifies a so-called reaction probability density that serves as the basis of his simulation algorithm (Equation 18). The main work of any Gillespie algorithm is therefore sampling a random pair (τ, μ) where

τ is the time to the next reaction, and
μ is the index of the next reaction (assuming each reaction is assigned a label 1, 2, 3, ..., m.

I am probably not offering anything new here, but I want to highlight the three main computational steps in Gillespie-like algorithms:

Sample from the joint probability distribution of (τ, μ).
Update the Markov chain state.
Compute the new joint probability distribution of (τ, μ).

The probability distribution changes because the jump rates change. Perhaps this is all obvious, but Step 3 is obscured in the literature and in implementations as a consequence of how the algorithms are described in articles. All the Gillespie-like algorithms can be summarized by how they implement Step 3, and so I believe what is needed is a flexible abstraction around this task.

The aggregator system with the new SSAStepper seem like they partially address this problem. I am proposing something of the following form (apologies for the bad names):

Have the aggregator track both the next jump index and the next jump time. The aggregator already "encodes" the probability distribution with cur_rates, sum_rate, and rates. All we need is the index.

mutable struct DirectJumpAggregation{T,F1,F2} <: AbstractJumpAggregator
  next_jump::Int    # the index of the next jump
  next_jump_time::T # was next_jump! sorry!
  end_time::T
  cur_rates::Vector{T}
  sum_rate::T
  rates::F1
  affects!::F2
  save_positions::Tuple{Bool,Bool}
end

Provide a function that samples from the distribution. Assume the aggregator already determined the next random pair:

# not really sampling...
sample_distribution(p::DirectJumpAggregation) = p.next_jump_time, p.next_jump

Provide a function that determines the next random pair. This is what distinguishes all the SSAs. New aggregators should implement. For the existing aggregator:

# maybe called `update_pair`?
function update_distribution!(p, integrator)
  rng_val = rand()

  p.sum_rate, ttnj = time_to_next_jump(integrator.u,integrator.p,integrator.t,p.rates,p.cur_rates)
  p.next_jump = searchsortedfirst(p.cur_rates,rng_val)

  p.next_jump_time = integrator.t + ttnj
  if p.next_jump_time < p.end_time
      add_tstop!(integrator,p.next_jump_time)
  end
end

Rewrite the current affect into a general procedure:

function (p::DirectJumpAggregation)(integrator) # affect!
  # 1. obtain the next random pair
  τ, μ = sample_distribution(p)

  # 2. update simulation by applying the affects for each jump
  @inbounds p.affects![i](integrator)

  # 3. update the jump rates and determine the next random pair
  update_distribution!(p, integrator)

  nothing
end

The benefits of this "abstraction" are:

New algorithms are new aggregators that implement an update_distribution! method. That's it.
The sample_distribution function does not care whether the algorithm couples or decouples the random pair. This is all handled by the update procedure, where it actually matters.
It may be possible to provide some default algorithms for some of the subtasks in "updating the distirbution". For example, most papers (and implementations) determine the next jump index using a linear search rather than bisection. This requires that cur_rates is a vector of the jump rates rather than a cumulative sum, but that's an implementation detail specific to a particular SSA. The point is, new SSAs can be built up by combining the smaller algorithms for specific tasks.
If (3) is possible, then we get new algorithms for free!

Relevant: #16

Retrodictive SSA

Keeping it from #8

https://www.sciencedirect.com/science/article/pii/S0021999110000501

Where to pass in additional network info?

The last thing we would need to support other SSAs is a mechanism for passing in additional network info (for example, a jump dependency graph).

For pure MassActionJump systems this actually doesn't need to be passed in; the jump type stores everything that is needed to build a dependency graph. So this could just be handled with the aggregate or initialize routines for a specific SSA. That said, if we allow it to be passed in it should be possible for a number of the new SSAs to support both MassActionJumps and ConstantJumps. While I don't myself need non-mass action systems, I like the idea of supporting both jump types within as many of the new SSAs as possible.

I see a few ways this could be handled; having the new SSA aggregator algorithms take this info as parameters when constructed, or extending JumpProblem to forward named arguments to aggregate and letting the specialized aggregate define parameter fields for the needed network info. But maybe there is a better way to handle this?

Then we could set everything up so that for a pure mass-action system everything is calculated for users, but for a hybrid system they are required to pass the necessary additional SSA info.

Initial time points with jumps are double saved

I am solving a Jump Diffusion SDE with the following code:

using DifferentialEquations
Nsim=10000;
Nsteps=100;
#Starting Point
r=0.02
sigma=0.2
T=1.0;
d=0.01;
u0=0.0;
#Drift part
f(u,p,t) = (r-d-sigma^2/2.0)
#Diffusion part
g(u,p,t) = sigma
#Time step
Dt =T/Nsteps;
#Time Window
tspan = (0.0,T)
#Definition of the SDE
prob = SDEProblem(f,g,u0,tspan)
#Definition of the Jump
rate(u,p,t) = 2
affect!(integrator) = (integrator.u = integrator.u/2)
jump = ConstantRateJump(rate,affect!)
#Definition of the Jump Problem
jump_prob = JumpProblem(prob,Direct(),jump)
#Definition of the MonteCarlo Problem
monte_prob = MonteCarloProblem(jump_prob)
sol = solve(monte_prob,SOSRI(),num_monte=Nsim,parallel_type=:threads,dt=Dt,adaptive=false)

When I look into the elements of the paths of the solution I get:

lengthT=length(sol[1].t) # =109 != Nsteps+1 and random
lengthUniqueT=length(unique(sol[1].t)) # =103 != Nsteps+1 and random

The problem is not only the repetitions of the elements (for example the initial time is repeated twice) but also the time discretization step that is not uniform.

The problem shows up also if I don't put the flag on the parallel type parallel_type=:threads, but I don't know which is the default argument for such flag.

Improve typing on tstops kwarg

Currently it's set to Float64[]. It should be tType[], but getting tType requires triangular dispatch, so this will have to be done in v0.6.

Create A Type for the Generate Rate Functions for Inlining

Anonymous functions don't inline, and this lack of inlining plus the fact that the generated rate equations are usually extremely simple is likely the reason why the jump equations are not at optimal speeds (about 2x-4x from Gillespie.jl). However, if a callable immutable type is created, its call can be inlined at the time of compilation, which could drastically speed up the rate and affect! computations, eliminating this difference. Then the only remaining difference would be the cost due to using a binary heap for tstop, which should be minimal, and the resulting code should compile to something close to an optimal "just a loop" Gillespie simulation.

Identical Ensemble solutions

Hi,
using the Ensemble interface with Jump equations yields multiple identical solutions when they should (most likely) not be.

using DifferentialEquations
using DiffEqBiological
using Plots

rn = @reaction_network ModelTest begin
    (10., 1.), 0 ↔ X
end

u0 = [0]
prob = DiscreteProblem(u0, (0., 100.), [])
jump_prob = JumpProblem(prob, Direct(), rn)
sol = solve(EnsembleProblem(jump_prob); trajectories=3)

plot(sol)

The output of this is three overlapping lines.

It works fine with SDE problems, so I guessed that the issue should be filed here.

Cheers!

Fictitious jump methods for time-dependent rates

There are thinning approaches that allow SSAs to "exactly" handle time-dependent rates by adding a fictitious jump to the system that makes the overall jump rate constant. Perhaps we could add support for one of these methods?

When I have some time I'll poke around for references, but please feel free to post any you might know of.

some SSAs giving segmentation faults

So far I've found RSSA, DirectCR now seem to crash on an inlined getindex call. MWE:

rn = @reaction_network begin
    β, S + I --> 2I
    γ, I --> R
end β γ
u0 = [999,1,0]; p = (0.1/1000,0.01); tspan = (0.,250.)
dprob2 = DiscreteProblem(rn, u0, tspan, p)
jprob2 = JumpProblem(dprob2, ssa, rn)

Nsims = 30000
function getmean(jprob,Nsims)
  m = 0.0
  for i = 1:Nsims
    sol = solve(jprob, SSAStepper())
    m += sol[end,end]
  end
  m/Nsims
end
m2 = getmean(jprob2,Nsims)

Stochastic splitting methods for regular leaping

http://www.math.kit.edu/iwrmm/seite/preprints/media/preprint%20nr.%2011-1.pdf

investigate dropping FunctionWrappers

In redoing benchmarks with stock Julia 1.1 I'm finding that on a Cartesian grid random walk simulation (32^3 mesh, so ~6*32^3 jumps) Direct consistently beats DirectFW. So it seems like Julia can now handle very deep recursion efficiently (how the jumps are evaluated to ensure type stability).

I'm not sure how to exploit this in methods using dependency graphs, where only a few jumps would be evaluated per time-step and we don't know them at compile time, but it would be great to be able to ditch function wrappers for ConstantRateJumps in the non-Direct SSAs.

SIR test gives funny answers

Hi @ChrisRackauckas

Your version of my SIR model in test/ gives funny answers: the epidemic never seems to take off, and it doesn't give plots similar to that in the introduction, even over multiple simulations.

Add support for passing multiple `MassActionJumps` to a `JumpSet`

As a convenience for users we should handle merging multiple MassActionJumps into one.

Performance

As came up over here: #9 (comment), the performance could use improvements. The basic design is simply using callbacks on ODE solvers, and for pure Gillespie models, using the Order 0 Discrete ODE solver. This works, gives interpolations, event handling, etc., and the biggest thing is that this makes it naturally work with ODEs/SDEs/PDEs/etc. So this is a good long-term solution and I don't want to do an alternative design. Instead, I want to find out how to make this work better.

The good thing is, we have @sdwfrost 's Gillespie.jl as a benchmark target. For very large equations we should be faster because we have a sparse representation of our jump stoichiometry vs their dense matrix. So let's not focus there. Instead we should focus on the performance of the small test cases vs Gillespie.jl. The test case from that comment is a good one.

There are two areas which likely need to be optimized. First of all, we are definitely hitting splatting penalties:

JuliaLang/julia#13359

I don't think we should work around what is hopefully a bug that should be patched up before 1.0. The other thing is we should check the generated code in the OrdinaryDiffEq solvers (and equivalently, StochsticDiffEq solvers which handle this almost the same) for the event handling. Essentially, a benchmark of jumps is just a benchmark of the differential equation solvers' event handling. There is nothing theoretically that should make this have a large performance gap from a tight loop here other than the fact that tstops, the structure for the next timepoint to hit, is a Binary Heap. This has to be a Binary Heap in order to allow event handling to push timepoints into tstops in an arbitrary order. However, I wouldn't expect it to be that bad...

So the aim is: don't lose any generality, but get within 2x of Gillespie.jl, likely closer to 1.5x. We may need to use generated functions for the rate equations to do so, and the splatting recursions may need to be checked. This should be doable.

Move abstract types to DiffEqBase

Zeroth order mass action jumps

Do we handle this in MassActionJump? @Msturroc

Tests against BioSimulator.jl

using DiffEqBase, DiffEqBiological, DiffEqJump, OrdinaryDiffEq,BenchmarkTools
rs = @reaction_network dtype begin
    k1, DNA --> mRNA + DNA
    k2, mRNA --> mRNA + P
    k3, mRNA --> 0
    k4, P --> 0
    k5, DNA + P --> DNAR
    k6, DNAR --> DNA + P
end k1 k2 k3 k4 k5 k6
tf        = 20000.0
num_save  = 25000
dtsave    = tf / num_save
params    =  (.5, (20*log(2.)/120.), (log(2.)/120.), (log(2.)/600.), .025, 1.)
prob      = DiscreteProblem([1,0,0,0], (0.0, tf), params)
# DiffEqsJumps Direct method, saving at fixed positions
jump_prob = JumpProblem(prob, Direct(), rs, save_positions=(false,false))
println("DiffEqJump, Direct, saving at ", num_save, " times:")
@btime solve($jump_prob, FunctionMap(), saveat=$dtsave)

# DiffEqsJumps Direct method, saving everything
jump_prob = JumpProblem(prob, Direct(), rs)
println("DiffEqJump, Direct, saving at all jump times:")
@btime solve($jump_prob, FunctionMap())

# DiffEqsJumps Direct method, saving at fixed positions
jump_prob1b = JumpProblem(prob, Direct(), rs, save_positions=(false,false))
println("DiffEqJump, Direct-SSAStepper, saving at ", num_save, " times:")
@btime solve($jump_prob1b, SSAStepper(), saveat=$dtsave)

# DiffEqsJumps Direct method, saving everything
jump_prob1b = JumpProblem(prob, Direct(), rs)
println("DiffEqJump, Direct-SSAStepper, saving at all jump times:")
@btime solve($jump_prob1b, SSAStepper())

@profile for i in 1:200 solve(jump_prob1b, SSAStepper(), saveat=dtsave) end
Profile.clear()
Juno.profiler()

tf        = 20000.0
num_save  = 25000
dtsave    = tf / num_save
using BioSimulator
model = Network("Gene Expression with Neg Feedback")
model <= Species("DNA", 1)
model <= Species("mRNA", 0)
model <= Species("P", 0)
model <= Species("DNAR", 0)
model <= Reaction("transcription", .5, "DNA --> mRNA + DNA")
model <= Reaction("translation", 20*log(2.)/120., "mRNA --> mRNA + P")
model <= Reaction("mRNA degradation", (log(2.)/120.), "mRNA --> 0")
model <= Reaction("protein degradation", (log(2.)/600.), "P --> 0")
model <= Reaction("DNA repression", .025, "DNA + P --> DNAR")
model <= Reaction("DNA freedom", 1., "DNAR --> DNA + P")
result  = simulate(model, algorithm=BioSimulator.SSA, time=tf, epochs=num_save, trials=1)
gui( plot(result.t_index, result[1,:,1], reuse=false) )
println("BioSimulator, Direct, saving at ", num_save, " times:")
@btime simulate($model, BioSimulator.SSA, time=$tf, epochs=$num_save, trials=1);

remake for JumpProblems

Broken as described here.

Automated depedency graph generation

We should hook into the new sparsity detection tools and/or extend them to generate jump dependency graphs from affect! functions. This would be great for users, and should allow the non-Direct SSAs to be used as broadly as Direct (for ConstantRateJumps at least).

Add discrete callbacks and tstops to SSAStepper

Recently requested by a collaborator who wants to be able to simulate cell division in SSA models at predetermined times.

Support for graphs of SSAs

Since this is being talked about in MT here, we should figure out how to make this work efficiently for jumps, preferably with a general API that allows us to compose different SSAs.

I think the primary challenge is that we need to support the ability for jumps at one location to also change species at neighboring locations on the graph, with data structures, rate functions and dependencies at that and other neighboring sites getting updated as appropriate. For spatial bio problems spatial jumps are usually representable as mass action reactions, i.e. A_i -> A_j or A_i + B_j -> C_k. In many cases though, for example Cartesian grids, it can be more efficient to do one jump for a diffusive hop away from a site, and then sample the particular placement site within the affect. These will often be the most executed jumps, sometimes as much as 99% of the events, so I would think performance will strongly depend on optimizations like we did for MassActionJumps.

Note, if the spatially-distributed jumps were all just thrown into one SSA there is nothing to do, but the idea is there are hierarchical SSAs where the top level stores next event times (or total propensities) for each site, and then each site has an associated SSA for determining what happens there. So we should aim to be able to support this type of composition.

sciml / jumpprocesses.jl Goto Github PK

jumpprocesses.jl's People

Contributors

Stargazers

Watchers

Forkers

jumpprocesses.jl's Issues

Recommend Projects

Recommend Topics

Recommend Org