get_concrete_problem does not preserve ExtendedJumpArray when it is initialized as an integer

[gzagatti]

JumpProcess.jl does some tricks to handle VariableJumps. However, I might have found a small bug in the current implementation.

Let's say we have the following model taken from the tutorial:

using JumpProcesses, OrdinaryDiffEq

p = (λ = 2.0, μ = 1.5)

deathrate(u,p,t) = p.μ * u[1]
deathaffect!(integrator) = (integrator.u[1] -= 1; integrator.u[2] += 1)
deathvrj = VariableRateJump(deathrate, deathaffect!)

rate1(u,p,t) = p.λ * (sin(pi*t/2) + 1)
affect1!(integrator) = (integrator.u[1] += 1)
vrj = VariableRateJump(rate1, affect1!)

function f!(du, u, p, t)
    du .= 0
    nothing
end
u₀ = [0., 0.]
oprob = ODEProblem(f!, u₀, tspan, p)
jprob = JumpProblem(oprob, Direct(), vrj, deathvrj)

sol = solve(jprob, Tsit5())

This works perfectly. However, if we modify u0 from a Float64 to an Int array u0 = [0, 0]. Then solve raises an exception which basically complains that there is no method matching (::JumpProcess)(::Vector{Float64}, ::Vector{Float64}, ...) since it is looking for a method matching (::JumpProcess)(::ExtendedJumpArray, ::ExtendedJumpArray, ...). This makes total sense since JumpProcess.jl remakes oprob to use an ExtendedJumpArray.

So what's going on? After some digging around the source code, it turns out that the culprit is located in this line of DiffEqBase.jl/src/solve.jl.

isadapt && eltype(u0) <: Integer && (u0 = float.(u0))

When dealing with jump processes, it is common to be interested in modelling counts so having u₀ as an array of integers is a reasonable choice.

It's probably the case that you only want to convert the extended part of the array to float. However, I'm not sure if that's possible.

In any case, you should preserve u0 as an ExtendedJumpArray that it works even after conversion.

[ChrisRackauckas]

I don't think you can, because the state has to be continuous in order to have a continuous integration for the jump time. ConstantRateJump is able to do this because by the constant rate you don't need an integration process in order to choose the jump times, but for variable rate you do, and for type-stability in the definition, if it's all contained in a single state array then it needs to be floating point valued. In theory, we could maybe make a split state form, though that would increase a lot of the type information and improve a lot of extra constraints on its usage (such as not being able to use a continuous rate dependent on the state values) and thus it wasn't done because it would be a lot of work for little gain.

[gzagatti]

I agree with you. In this case, I would just fix this issue by maintaining the ExtendedJumpArray type when converting to Float. Would that make sense?

[ChrisRackauckas]

Oh it doesn't? Yeah that's an issue. That means that the broadcast similar needs an overload: i.e. broadcast in general is returning an Array, so that should get fixed.

[gzagatti]

It doesn’t. That’s why I opened the issue. But I think I overcomplicated my initial explanation. Apologies for that.

[ChrisRackauckas]

https://docs.julialang.org/en/v1/manual/interfaces/#man-interfaces-broadcasting

Specifically this section:

https://docs.julialang.org/en/v1/manual/interfaces/#Selecting-an-appropriate-output-array

That must be missing from the definition of the broadcast in JumpProcesses. I'm transferring it over to JumpProcesses.jl as this is pretty contained to there.

[isaacsas]

@gzagatti seems closed by https://github.com/SciML/JumpProcesses.jl/pull/340.

[gzagatti]

Awesome! Thanks @meson800!