Example for CLE where paths go negative

[isaacsas]

Explain about using absolute values of rates, but that this doesn't guarantee the solution stays positive.

[isaacsas]

We need to document that we currently wrap in absolute values, and perhaps provide an option to not do this?

[TorkelE]

Currently, things can go negative:

using Catalyst, StochasticDiffEq, Plots
rn = @reaction_network begin
    (p,1.0), 0 <--> X
end
sprob = SDEProblem(rn, [1.0], (0.0, 10.0), [1.0])
sol = solve(sprob, ImplicitEM())
plot(sol)

However, the PositiveDomain callback gives error:

using DifferentialEquations
sol = solve(sprob, ImplicitEM(); callback=PositiveDomain())    # Both give the same error.
sol = solve(sprob, ImplicitEM(); callback=PositiveDomain([1.0]))    # Both give the same error.

ERROR: ArgumentError: matrix contains Infs or NaNs
Stacktrace:
  [1] chkfinite
    @ ~/.julia/juliaup/julia-1.9.2+0.x64.linux.gnu/share/julia/stdlib/v1.9/LinearAlgebra/src/lapack.jl:86 [inlined]
  [2] generic_lufact!(A::Matrix{Float64}, pivot::LinearAlgebra.RowMaximum; check::Bool)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.9.2+0.x64.linux.gnu/share/julia/stdlib/v1.9/LinearAlgebra/src/lu.jl:135
  [3] generic_lufact!
    @ ~/.julia/juliaup/julia-1.9.2+0.x64.linux.gnu/share/julia/stdlib/v1.9/LinearAlgebra/src/lu.jl:133 [inlined]
  [4] do_factorization
    @ ~/.julia/packages/LinearSolve/LD2dF/src/factorization.jl:65 [inlined]

[ChrisRackauckas]

That's unrelated to positive domain callback, it's from an older version of Julia where chkfinite wasn't disabled by check=false.

But also no, this is not a problem where PositiveDomain callback makes sense. The first requirement of using a domain callback is that the domain constraint has to be satisfied by the solution. The whole purpose of this issue is that the CLE does not always emit a solution that is positive, so the callback is not a solution.

IIRC Catalyst always does abs on the terms in the sqrt to handle this. @TorkelE 's thesis work had models which requires this. Using abs will give a biased distribution of course, but that's because the CLE definition breaks down when the solution approaches zero. At least the abs version gives a definition that is similar to a reflection around zero, which I'd argue is probably the best definition you could give to it, and we can have an option to error instead as the only other option I think is really viable.

[isaacsas]

I guess we could implement a reflection callback, as that is an approach that has been used, but then the model one is simulating is no longer the CLE.

[isaacsas]

(So I don't think it really makes sense to push such an approach.)

[TorkelE]

Would it be possible to make the default having nothing, but having a sqrt with a modified error message that explains what is going on and suggests a couple of ways to deal with it (while admitting that these will produce a non-true CLE)? And then we could provide SDEProblem options for whatever the users want (e.g. abs values or setting certain callbacks).

[isaacsas]

I don't think we'd want to wrap Julia's sqrt for performance reasons.

[TorkelE]

suspected so. Let's just keep it as it is and add options for alternative ways of dealing (or not dealing) with it