`tstop` error with Piecewise Deterministic Markov Processes

Specific Domains Modelling & Simulations
, ,
1

I’m asking for help for an annoying, (ironically) non-deterministic error when I solve Jumpy-ODE’s (i.e. “Piecewise Deterministic Markov Processes”). It tends to start happening as I increase the number of processes and does not occur for similar smaller systems. or short times.

>>> p0 = ... # a jump problem with 128 independent jumps acting on top of an ODE problem of 64 dimensions
>>> out = solve(p0, Tsit5())

ERROR: Tried to add a tstop that is behind the current time. This is strictly forbidden

It’s a bit hard to provide reproducing code because of the size of the jump problem that I’ve defined. I can provide code where I define the constantrate jump problem. Not enough to run, just showing how I build the JumpSet.

function approx_piecewise_jump_problem(model; tspan  = (0.0,1.0), u0 = ones(model.nstates) )
    odeprob = theta_ode_prob(model; tspan, u0)
    jumps = []
    # selection increasing
    for i in 1:Int(model.nstates)
        # (u+ θ) * σ
        rate(u,p,t) = (2*model.Ne*fitness(model, t, i) 
            # + 2 + sum(mod(u[j],1)  - 2*model.Ne*mutation_matrix(model, t, i, j) for j in 1:model.nstates if j != i) 
            # mutation contribution
            ) * u[i] 
        function affect!(integrator)
            integrator.u[i] +=  1.0
            nothing
        end
        push!(jumps,ConstantRateJump(rate, affect!))
    end
    # drift decreasing
    for i in 1:Int(model.nstates)
        rate(u,p,t) = floor(u[i])*(sum(u) - 1 )
        function affect!(integrator)
            integrator.u[i] += -1.0
            nothing
        end
        push!(jumps,ConstantRateJump(rate, affect!))
    end
    jset = JumpSet(; constant_jumps=jumps)
    return JumpProblem(odeprob, Direct(), jset)
end

The underlying ODE is a well behaved quasi-linear problem. (You can tell by my 2*model.Ne that I am doing population genetics.)

Googling this error only yields some issues with thread safety in ensembles, but I’m only single threading for now. The stack trace is the usual unhelpful ODE type dumps, but hits the relatively standard code in the Gillespie thingy.

1. error(s::String) at error.jl
2. add_tstop!
3. register_next_jump_time! at ssajump.jl
4. AbstractSSAJumpAggregator at ssajump.jl
5. apply_discrete_callback! at callbacks.jl
6. handle_callbacks!
7. _loopfooter!
8. loopfooter
9. solve

Should I file an issue, or is a problem on my end? If it’s a DifferentialEquations.jl issue, can anyone think of a workaround in the meantime (alternate sovlers maybe?) while I come up with an MRE for what seems to be a subtle issue?

1 Like
2

Might be completely unrelated, but your affect! functions being defined every iteration inside that loop looks kind of weird. They are not a closure over anything, so you can probably define a plusone_affect! and a minusone_affect! function outside those loops but use them inside the loops.

Is your first line then something like p0 = JumpProblem(odeproblem, DIrect(), jump1, jump2, ..., jump128)?

Additionally, I’m not sure if your rate functions being defined in a way that allows backward timesteps (just speculating) could cause this, unless the jump counter function (the underlying Poisson process) in DifferentialEquations.jl explicitly prevents this whatever rate function it is provided with.

3

Hard to debug this without a minimum working example, but have you checked your rate functions are always returning non-negative values? If they ever return a negative value it could definitely cause problems.

1 Like
4

@tom-plaa they capture the variable index from the loop. There are different jump rates for each dimension.

@isaacsas thanks for the idea. Nothing should go negative… but the integrator might not stay in it’s proper domain, especially in high-D. I’ll try enforcing the domain more carefully.

5

I would add an explicit check to your rate functions that they aren’t returning a negative value. Then if you get this error and the check isn’t firing we know that isn’t the issue.

1 Like
6

you could also try implicit solvers and more precise ones, set the tolerances, etc. This is really important when simulating PDMP, this is kinda obvious:

https://rveltz.github.io/PiecewiseDeterministicMarkovProcesses.jl/dev/solver/

If the flow between jumps is analytical, you can really speed up the simulation but it rarely is.

1 Like
7

Hate to own up to this one, but yes, one of my rates was negative :clown_face: So moral of the story kids, if you see

ERROR: Tried to add a tstop that is behind the current time. This is strictly forbidden

check the sign of your rates. :sweat_smile:

thanks @isaacsas the error(“rate went negative”) check worked perfectly, and @rveltz the link you posted was also useful in other ways.

8

:+1: