Continuous callbacks checking on u works on standard ODE version but does not work in the Jump Process

Specific Domains Modelling & Simulations
,
1

I am trying to solve a system both deterministically and stochastically.

The solutions match for my system in both cases without callback. But the callback does not seem to be working.

function condition1(u,t,integrator)
   u[1]-100
end

function affect1!(integrator)
    integrator.u[1] =1000
end
function condition2(u,t,integrator)
   u[2]-100
end

function affect2!(integrator)
    integrator.u[2] =1000
end
cb1=ContinuousCallback(conditions,affect1!)
cb2=ContinuousCallback(conditions,affect2!)
cbs=CallbackSet(cb1,cb2)
sol=solve(prob;callback=cbs)

This works for standard ODE problem. But when I try to do the same for the stochastic version It just does not get activated.
What am I doing wrong?

jprob=JumpProblem(prob,umap,tspan,pmap)
sol=solve(jpron,callback=cbs)

I saw in the docs a way to use time related callback.. But I want to change this u value to make sure I have enough of the reactant in the system forever. I feel like I am missing something fundamental here about JumpProcesses.

Thanks in advance.

2

Can you post a full MWE that can be run when cut and pasted?

3

Sure, Here you go. The following example programmatically generate a vector of species using same reactants A and B.

using Catalyst,OrdinaryDiffEq,Plots, JumpProcesses

function generate_creation(n)
    # Creates a vector `X` with n species.
    @species X(t)[1:n]
	@species A(t),B(t)
	@parameters k
#	@parameters α,K,m,k
    rxs = [
		[Reaction(k, [A,B], [X[i]]) for i in 1:n];
    ]
    # Assembly and return a complete `ReactionSystem` 
    @named lp_n = ReactionSystem(rxs, t)
    return complete(lp_n)
end

function condition1(u, t, integrator) # Event when condition(u,t,integrator) == 0
    u[1]-100
end
function condition2(u, t, integrator) # Event when condition(u,t,integrator) == 0
    u[2]-100
end

function affect1!(integrator)
    integrator.u[1] = 1000
	reset_aggregated_jumps!(integrator)
end
function affect2!(integrator)
    integrator.u[2] = 1000
	reset_aggregated_jumps!(integrator)
end

begin
	cb1 = ContinuousCallback(condition1, affect1!)
	cb2 = ContinuousCallback(condition2, affect2!)
	cbs = CallbackSet(cb1,cb2)
end

N=3
params  = (:k =>0.1)
X0 = fill(0.0, N)
X0[1] = n * 1.0
prob1=generate_creation(N)
u_0 = [prob1.X => X0,prob1.A => 1000, prob1.B => 1000]

odeprob1=ODEProblem(prob1, u_0, (0.0, 0.5), params)
sol1=solve(odeprob1)
plot(sol1)

Final plot commands produced expected result. Where X1,X2 and X3 increases untill they saturate at same value as A and B goes to zero.

image
image600×400 14.5 KB

Adding Callback and solving gives us:

sol2=solve(odeprob1;callback=cbs)
plot(sol2)

image
image600×400 25.7 KB

So callbacks are working fine that lets the system generate X1,X2 and X3 well.

Now on to the Jump Process stochastic Equivalent:

begin
	jumpinputs = JumpInputs(prob1, u_0, (0,0.05), params)
	jumpprob = JumpProblem(jumpinputs)
	# now, let's solve and plot the jump process:
	soljump = solve(jumpprob,callback=cbs)
	plot(soljump)
end

image
image600×400 18.3 KB

As you can see, it shows the right behaviour but the callbacks are not working. My instinct tells me that the condition u==0 is not working out at any time. But since, this is a stochastic simulation, the values are jumping by one. It should work.
P.S. I have adapted these examples from the SciMl documentation. Except the callback for jumpprocess. The documentation seems to be lacking there.

4

I have instead of using continuous callback, used discrete callback and managed to make it work. But why it is not working still seems to be confusing to me.

5

@ChrisRackauckas I tried a pure JumpProcesses version of this, and strangely while the condition in the callback shows as true for a bunch of timesteps (as it should) the affect seems to never get called. Should that be happening?

Here is my MWE

using JumpProcesses, OrdinaryDiffEqTsit5

# jump
rate(u, p, t) = 1.0
affect(integ) = (integ.u[1] += 1; nothing)
vrj = VariableRateJump(rate, affect)

# continuous event
cond(u, t, integ) = (@show (u[1] - 10 == 0); return  u[1] - 10)
affect1(integ) = (@show "here"; integ.u[1] = 1; reset_aggregated_jumps!(integ); nothing)
cb = ContinuousCallback(cond, affect1)

oprob = ODEProblem((du,u,p,t) -> (du .= 0), [1.0], (0.0, 50.0))
jprob = JumpProblem(oprob, Direct(), vrj)
sol = solve(jprob, Tsit5(); callback = cb)

u[1] - 10 == 0 ends up showing as true for ~10 timesteps, but “here” is never printed.

6

Hello,

Since the jump process is discrete a discrete callback should suffice. We do not have a continuous crossing over zero to detect in the ContinuousCallback.
This worked for me.

using JumpProcesses, OrdinaryDiffEq

# jump
rate(u, p, t) = 1.0
affect(integ) = (integ.u[1] += 1; nothing)
vrj = VariableRateJump(rate, affect)

# continuous event
cond(u, t, integ) = (@show (u[1] - 10 == 0); return  u[1] - 10 == 0)
affect1(integ) = (@show "here"; integ.u[1] = 1; reset_aggregated_jumps!(integ); nothing)
cb = DiscreteCallback(cond, affect1)

oprob = ODEProblem((du,u,p,t) -> (du .= 0), [1.0], (0.0, 50.0))
jprob = JumpProblem(oprob, Direct(), vrj)
sol = solve(jprob, Rodas5(); callback = cb)
7

Yes, using a DiscreteCallback does make more sense, but I would have thought the continuous callback should still fire given the condition does actually evaluate to zero.

8

Yes Thanks @sass-jacob I used a DiscreteCallback in my implementation as well. Actually I am using a combination of callbacks in my problem. But Discrete one works for this case as the number of “molecules” increases discretely.

What I do not understand @isaacsas mentioned if the condition is showing up as true. Why is continuous callback not working. Thanks to both of you for this discussion and implementing your own MWE.

9

I don’t know the internal details of the callback selection, evaluation, and rootfinder that are used, but the rootfinder presumably requires continuity (which isn’t true here as was previously mentioned). That said, given the condition is clearly being called and evaluating to 0.0 I would have thought it would still fire the callback. Chris or someone who has worked on the callback handling would probably know why the affect is never called in this situation.