Positive domain enforcement (again, but new area of relevance)

[TorkelE]

Mostly a repeat of this one: #88. However, that one mostly discussed for SDEs, and I got some ODE cases where this became quite relevant. When doing parameter fitting on a model containing a hill function

using Catalyst
rn = @reaction_network begin
    hill(X, 1.0, K, n), ∅ --> X
    d, X --> ∅
end

one sometimes encounters parameter sets

u0 = [:X => 0.2]
ps =  [:K => 9.087101071878903, :n => 2.5163255503214117, :d => 4.638316464318106]

which will, when simulated, cause the ODE (due to numeric inaccuracy) to go slightly negative

oprob = ODEProblem(rn, u0, 20.0, ps)
solve(oprob)

This causes a

ERROR: DomainError with -1.6572261631222257e-6:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

in the solve. command. Because the PositiveDomain callback do not prevent evaluations outside of the domain, this is not a solution

using DiffEqCallbacks
solve(oprob; callback = PositiveDomain([1.0])) # Still errors.

This definitely feels like something that can come up without the user trying to do something outrageous. Not sure if there is a recommended solution, but if there is, it feels like something to create a tutorial on.

[ChrisRackauckas]

What's the error? Since we use NaNMath these days I don't think it would error in the solve

[TorkelE]

Interesting. The minimal example is

using Catalyst, OrdinaryDiffEqDefault
rn = @reaction_network begin
    hill(X, 1.0, K, n), ∅ --> X
    d, X --> ∅
end
u0 = [:X => 0.2]
ps =  [:K => 9.087101071878903, :n => 2.5163255503214117, :d => 4.638316464318106]
oprob = ODEProblem(rn, u0, 20.0, ps)
solve(oprob)

which gives a

ERROR: DomainError with -1.6572261631222257e-6:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
  [1] throw_exp_domainerror(x::Float64)
    @ Base.Math ./math.jl:41
  [2] ^(x::Float64, y::Float64)
    @ Base.Math ./math.jl:1157
  [3] hill
    @ ~/.julia/dev/Catalyst/src/registered_functions.jl:42 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/SymbolicUtils/02FYV/src/code.jl:440 [inlined]
  [5] macro expansion
    @ ~/.julia/packages/Symbolics/8VseH/src/build_function.jl:564 [inlined]
  [6] macro expansion
    @ ~/.julia/packages/SymbolicUtils/02FYV/src/code.jl:397 [inlined]
  [7] macro expansion
    @ ~/.julia/packages/Symbolics/8VseH/src/build_function.jl:348 [inlined]
  [8] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [9] macro expansion
    @ ./none:0 [inlined]
 [10] generated_callfunc
    @ ./none:0 [inlined]
 [11] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Vector{…}, ::MTKParameters{…}, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:150
 [12] macro expansion
    @ ~/.julia/packages/ModelingToolkit/CvDvM/src/systems/codegen_utils.jl:0 [inlined]
 [13] _generated_call
    @ ~/.julia/packages/ModelingToolkit/CvDvM/src/systems/codegen_utils.jl:262 [inlined]
 [14] GeneratedFunctionWrapper
    @ ~/.julia/packages/ModelingToolkit/CvDvM/src/systems/codegen_utils.jl:259 [inlined]
 [15] Void
    @ ~/.julia/packages/SciMLBase/sYmAV/src/utils.jl:486 [inlined]
 [16] (::FunctionWrappers.CallWrapper{…})(f::SciMLBase.Void{…}, arg1::Vector{…}, arg2::Vector{…}, arg3::MTKParameters{…}, arg4::Float64)
    @ FunctionWrappers ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:65
 [17] macro expansion
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:137 [inlined]
 [18] do_ccall
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:125 [inlined]
 [19] FunctionWrapper
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:144 [inlined]
 [20] _call
    @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:12 [inlined]
 [21] FunctionWrappersWrapper
    @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:10 [inlined]
--- the above 7 lines are repeated 1 more time ---
 [29] ODEFunction
    @ ~/.julia/packages/SciMLBase/sYmAV/src/scimlfunctions.jl:2468 [inlined]
 [30] perform_step!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::OrdinaryDiffEqTsit5.Tsit5Cache{…}, repeat_step::Bool)
    @ OrdinaryDiffEqTsit5 ~/.julia/packages/OrdinaryDiffEqTsit5/DHYtz/src/tsit_perform_step.jl:199
 [31] perform_step!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::OrdinaryDiffEqCore.DefaultCache{…}, repeat_step::Bool)
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/vS7Uo/src/perform_step/composite_perform_step.jl:141
 [32] perform_step!
    @ ~/.julia/packages/OrdinaryDiffEqCore/vS7Uo/src/perform_step/composite_perform_step.jl:138 [inlined]
 [33] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/vS7Uo/src/solve.jl:620
 [34] #__solve#62
    @ ~/.julia/packages/OrdinaryDiffEqCore/vS7Uo/src/solve.jl:7 [inlined]
 [35] __solve
    @ ~/.julia/packages/OrdinaryDiffEqCore/vS7Uo/src/solve.jl:1 [inlined]
 [36] solve_call(_prob::ODEProblem{…}, args::OrdinaryDiffEqCore.CompositeAlgorithm{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PbBEl/src/solve.jl:635
 [37] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::MTKParameters{…}, args::OrdinaryDiffEqCore.CompositeAlgorithm{…}; kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PbBEl/src/solve.jl:1128
 [38] solve_up
    @ ~/.julia/packages/DiffEqBase/PbBEl/src/solve.jl:1106 [inlined]
 [39] #solve#42
    @ ~/.julia/packages/DiffEqBase/PbBEl/src/solve.jl:1045 [inlined]
 [40] solve
    @ ~/.julia/packages/DiffEqBase/PbBEl/src/solve.jl:1033 [inlined]
 [41] #__solve#3
    @ ~/.julia/packages/OrdinaryDiffEqDefault/kOV75/src/default_alg.jl:48 [inlined]
 [42] __solve
    @ ~/.julia/packages/OrdinaryDiffEqDefault/kOV75/src/default_alg.jl:47 [inlined]
 [43] #__solve#63
    @ ~/.julia/packages/DiffEqBase/PbBEl/src/solve.jl:1437 [inlined]
 [44] __solve
    @ ~/.julia/packages/DiffEqBase/PbBEl/src/solve.jl:1428 [inlined]
 [45] #solve_call#35
    @ ~/.julia/packages/DiffEqBase/PbBEl/src/solve.jl:635 [inlined]
 [46] solve_call(::ODEProblem{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PbBEl/src/solve.jl:592
 [47] solve_up(::ODEProblem{…}, ::Nothing, ::Vector{…}, ::MTKParameters{…}; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PbBEl/src/solve.jl:1112
 [48] solve_up
    @ ~/.julia/packages/DiffEqBase/PbBEl/src/solve.jl:1106 [inlined]
 [49] solve(::ODEProblem{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PbBEl/src/solve.jl:1043
 [50] solve(::ODEProblem{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PbBEl/src/solve.jl:1033
Some type information was truncated. Use `show(err)` to see complete types.

[ChrisRackauckas]

You're registering hill without NaNMath right? That would be a Catalyst problem.

[TorkelE]

Yeah that could be what is going on, sounds like something we should update

[TorkelE]

is it that when we do

hill(X, v, K, n) = v * (X^n) / (X^n + K^n)
@register_symbolic hill(X, v, K, n);

that should be switched to

hill(X, v, K, n) = v * NaNMath.pow(X, n) / (NaNMath.pow(X, n) + NaNMath.pow(K, n))
@register_symbolic hill(X, v, K, n);

?

[ChrisRackauckas]

Yes

[isaacsas]

Couldn't we codegen to allow the use of PositiveIntegrators? That feels like a more robust approach to this issue. (I haven't check if general reaction systems meet their criteria, but if they do we should support using that package.)

[isaacsas]

If I remember right Burrage has SDE versions of such integrators for reaction systems too (there were a few J. Chem Phys papers).

[ChrisRackauckas]

You still need the NaNMath change.

[isaacsas]

Does using NaNmath have performance implications? I'd prefer this to be a user selectable option rather than switching away from the Base versions as the default.

[ChrisRackauckas]

No, it literally just NaNs instead of errors.

I'd prefer this to be a user selectable option rather than switching away from the Base versions as the default.

Why? That is never the desired behavior. We will likely make this easier by using method overlays in the ODE solver to do this automatically at the compiler level, but I don't see why you'd try to do anything in the opposite direction. Is there any good that can happen from it?

[TorkelE]

This specific case is no longer a problem after the update to the functions like hill. If you want @isaacsas we can close this one and just rely on the old positive domain issue for SDEs.

[ChrisRackauckas]

What's the remaining issue on SDEs? There's one issue mathematically that there's a breakdown that the SDE approximation does not make sense as the concentration goes to zero. But at least in terms of the solver, appropriately using NaNMath will at least let it solve the resulting SDEs if sqrt(abs(x))