search

SciML / JumpProcesses.jl

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)
https://docs.sciml.ai/JumpProcesses/stable/
Other
136 stars 35 forks source link

Type conversion error when adding a `VariableRateJump` to an `EnsembleProblem` and specifying `saveat` #322

Closed benmcdonough20 closed 11 months ago

benmcdonough20 commented 1 year ago

This problem seems to regard several parts of the DifferentialEquations.jl ecosystem, so I hope this is the right place to post it!

An error is thrown when an EnsembleProblem is created using a JumpProblem which includes aVariableRateJump and the saveat keyword is specified. Here is a minimum working example (using the base problem from the documentation here: 

using DifferentialEquations

function f!(du, u, p, t)
    du[4] = u[2] * u[3] /100000 - u[1]*u[4]/100000
    nothing
end
u0 = [999.0, 10.0, 0.0, 100.0]

p = (.1/1000, .01)
tspan = (0.0, 250)
save_times = range(tspan...; step = 1)

variable_rate(u, p, t) = 1e-2*u[4]
function variable_affect!(integrator)
    integrator.u[2] += 1
    nothing
end

jump = VariableRateJump(variable_rate, variable_affect!)

#breaks if saveat is specified and interpolation is required
prob = ODEProblem(f!, u0, tspan, p, saveat = save_times)
jump_prob = JumpProblem(prob, Direct(), jump)
ensemble_prob = EnsembleProblem(jump_prob)
sol = solve(ensemble_prob, Tsit5(); trajectories = 2)

Here is the stack trace: (pastebin)

The issue appears to be that JumpProcesses keeps up with the variable rate by appending the rates to the existing solution to create an ExtendedJumpArray. The problem only appears when saveat is specified because the integrator tried to interpolate to time in between the timesteps, and the interpolate method returns the ExtendedJumpArray and attempts to save it alongside the raw solution in an output array. I was able to implement a hack-ish fix by manually type checking and replacing the copyat_or_push! method here:

https://github.com/SciML/OrdinaryDiffEq.jl/blob/f56981814afbc058a95471d2f1b8ac7abdb8ff41/src/integrators/integrator_utils.jl#LL69C9-L75C84

I am new to Julia, so please let me know if I can provide more info!

isaacsas commented 1 year ago

I'm going to move this to JumpProcesses since it appears to be another issue with the ExtendedJumpArray interface needing custom dispatches (as we also see in https://github.com/SciML/JumpProcesses.jl/issues/321).

If you think you have a fix via adding a custom dispatch for ExtendedJumpArrays, a PR that adds it to JumpProcesses in

https://github.com/SciML/JumpProcesses.jl/blob/master/src/extended_jump_array.jl

with a test would be very much appreciated!

isaacsas commented 11 months ago

Closed by https://github.com/SciML/JumpProcesses.jl/pull/340