In short: When making solving SDEs the noise function seems to be called twice in every timestep, with the variable values being changed in between. This causes out of domain errors when my positive_domain callback cannot ensure that all values are positive between the two calls to the noise function. First I describe the code I used when encountering the problem, next a more minimal example.
The origin of this problem is that I try to make simulations of a biochemical reaction network, having a
positive_domain callback which in each timestep checks whenever my solutions have become negative and fixes that (in case noise makes the solution negative, which might cause a crash during the noise calculation when the method takes sqrt of something negative).
In principle my method looks something like this:
using DifferentialEquations kBw = 3600; kDw = 18; kD = 18 kB1 = 3600; kB2 = 3600; kB3 = 3600; kB4 = 1800; kB5 = 3600; kD1 = 18; kD2 = 18; kD3 = 18; kD4 = 1800; kD5 = 18; kK1 = 36; kK2 = 36; kP = 180; kDeg = 0.7; v0 = 0.4; F = 30; K = 0.2; λW = 4; λV = 4.5; η = 0.101 p = [kBw, kDw, kD, kB1, kB2, kB3, kB4, kB5, kD1, kD2, kD3, kD4, kD5, kK1, kK2, kP, kDeg, v0, F, K, λW, λV, η] @reaction_func vB(σB) = v0*((1+F*σB)/(K+σB)) my_model = @reaction_network rnType η begin kDeg, (w,w2,w2v,v,w2v2,vP,σB,w2σB) ⟶ ∅ kDeg, vPp ⟶ phos (kBw,kDw), 2w ⟷ w2 (kB1,kD1), w2 + v ⟷ w2v (kB2,kD2), w2v + v ⟷ w2v2 kK1, w2v ⟶ w2 + vP kK2, w2v2 ⟶ w2v + vP (kB3,kD3), w2 + σB ⟷ w2σB (kB4,kD4), w2σB + v ⟷ w2v + σB (kB5,kD5), vP + phos ⟷ vPp kP, vPp ⟶ v + phos vB(σB), ∅ ⟶ σB λW*vB(σB), ∅ ⟶ w λV*vB(σB), ∅ ⟶ v end kBw kDw kD kB1 kB2 kB3 kB4 kB5 kD1 kD2 kD3 kD4 kD5 kK1 kK2 kP kDeg v0 F K λW λV function positive_domain() function condition(u,t,integrator) minimum(u) < 0 end function affect!(integrator) integrator.u .= integrator.uprev end return DiscreteCallback(condition,affect!,save_positions = (true,true)) end function equi(model, p) prob = ODEProblem(model,fill(1.,length(model.syms)),(0.,100.),p); sol = solve(prob,Rosenbrock23()) return max.(sol[end],0.0) end; function stochsim(model, p, tspan) u0 = equiN(model,p) prob_sde = SDEProblem(model,u0,tspan,p) sol = solve(prob_sde,ImplicitEM(),callback=positive_domain(),saveat=0.001); end p[end] = 0.001 sol = stochsim(my_model,p,(0.,2));
p[end] = 0.001 is a factor linearly scaling the amount of noise in the model. Here is where the weird stuff is starting. If I run the code it works fine if I set
p[end] = 0.001 or
p[end] = 0.1. However, if I set
p[end] = 0.01 I receive
DomainError (not every time, but about every other time I run it). I have traced this error to some variable being pushed into negative value by noise, and then the square root of this value is taken when the amount of noise is calculated. The square root of the negative number yields the error. The problem is that the positive_domain callback should catch this.
Some investigation gives that this is what happens:
- The callback is called at the beginning of every time step, ensuring that the value is non-negative.
- In every time step the noise function,
g, is called twice, with the value of
u(the variables) changing in between. Sometimes one value in
ubecomes negative somewhere in between the first and the second call to
g, the second call to
ggiving the out of domain error.
Ignoring the case where I found the problem and looking at a minimal example I can use this code:
using DifferentialEquations function f(du,u,p,t) du .= 1.01u end function g(du,u,p,t) println("Calling g at ",t) println(u,"\n") du[1,1] = 0.3u du[1,2] = 0.6u du[1,3] = 0.9u du[1,4] = 0.12u du[2,1] = 1.2u du[2,2] = 0.2u du[2,3] = 0.3u du[2,4] = 1.8u end prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=zeros(2,4)) sol = solve(prob,ImplicitEM());
to confirm that
g is called twice at every timestep, with the value of
u changed in between. Thereby causing problems is one value in
u is made negative after the first
g call, but without the
positive_domain callback having a chance to fix it before
g is called again.
Is this supposed to be like this? If so, is there a way around it to prevent the domain error?