Overflow problems in Hill functions

[devmotion]

The current implementation of hill etc. is problematic for large arguments.

For instance,

julia> hill(1e40, 1, 1e40, 10)
NaN

I tend to use something like

function myhill(x, v, y, n)
    if x < y
        z = (x / y)^n
        return v * (z / (1 + z))
    else
        return v / (1 + (y / x)^n)
    end
end

which returns

julia> myhill(1e40, 1, 1e40, 10)
0.5

---

An additional question: What's the reason for symbolically registering these functions instead of traversing them symbolically?

[TorkelE]

Sounds like a good point, never really looked at large inputs. @isaacsas probably have a lot more thoughts when it comes to high-end performance though, and probably have more thoughts than me when it comes to this.

[isaacsas]

I agree that is a more robust implementation. Honestly those definitions are historical, and go back to whoever added them in DiffEqBiological long long ago. Also, when switching to MTK we originally weren't registering, so the expressions were getting substituted, and the current form gives the "standard" ODE terms when Latexified.

I believe the reason we switched to registering them was to avoid them being eagerly substituted into. There has been discussion at times about trying to special case Hill functions in solvers such as JumpProcesses (like we do with mass action reactions), and that would require the ability to know about the presence of such functions in symbolic expressions.

I think @TorkelE added the ability to expand Hill functions to their symbolic expression in one or more places, so we could perhaps offer a system option to eagerly expand them if that would be useful.

If you want to make a PR to update the registered functions to be more numerically robust, with appropriate tests, that would be great and appreciated!

PS - Out of curiosity, why aren't you changing units so that the populations are smaller? In most cases I would tend to adjust units rather than run simulations with such very large populations.

[isaacsas]

If we do change the definition, and want to allow eager substitution, we'd need to use ifelse right?

[TorkelE]

Depending on how we do change things, we might need to be careful regarding homotopy continuation and structural identifiability analysis which depends on preserving the algebraic expressions. Just something to keep in mind.

[isaacsas]

If we stick to using a registered function that shouldn't be an issue though right? (Since we can just ensure in those contexts we substitute the current formula.)

edit: should -> shouldn't

[TorkelE]

yes, I think that should work out