Many (most?) of the SciML libraries feature a parameter p in their high-level API.
I would like to learn more about the motivation for this design, and the actual specs for what it can be and what it is intended for. Eg in some examples, it is a scalar, in some others, a vector, is there a constraint in the API describing what is allowed?
Or is it just up to the model family being described?
Why not fold it into the model spec instead?
(If there is some docs I am just missing, kindly point me to that)
It’s up to the model being described. It is very often useful to be able to define a family of models rather than a single model and allowing a user specifiable parameter object makes that easy.
That said if you only are interested in a single problem you can just ignore p and not use it. That’s always totally valid.
As a related note, I think that if anyone is taking design from SciML’s diffeq solvers for a new solver ecosystem, I think it would be really interesting to not expose t to users (if you want a time component, you can always add du[i] = 1 as an equation).
The advantage of this is that by making time implicit, you move all of the stiffness into the jacobian matrix, and simplify the the order conditions needed. (the downside of course is that you confuse engineers pretty badly, but IMO time might belong at the MTK/dyad level rather than the solver level)
If I understand you correctly, you’re saying that any non-autonomous system of ODEs can be treated as an autonomous system in that way? Kind of like any higher-order ODE can be re-formulated as first-order, which makes this the default form for (I guess) almost all solver libraries.
Do you know of any libraries/solvers where this implicit time dependence is the default? Or any references that study the advantages of this approach?
Regarding the original question: For ODEs, I always found it very convenient to use p to keep the model description as generic as possible without losing out on performance. More often than not, the model I’m interested in has some kind of parameters, like reaction rates, time scales, diffusion rates, etc. and I want to solve the system for many different combinations of those parameters. remake is very useful for that – you can just switch out the underlying parameter object p without re-creating the while ODEProblem (or whatever) for every parameter set.
And in some cases where the parameters of a model should be inferred or fitted, the p of the ODE function directly becomes the u of an optimization problem, so to speak.
As far as I understand there are no formal constraints on p since the user is passing p to the solver and also writing the model function that “consumes” p in the end. So as long as your function can deal with the actual p you supply, you should be fine (?)
It’s really for the AD system. If you’re going to write rules and handling for when people want to differentiate, you need to actually have it there to handle it. For example, if you differentiate w.r.t. p, even with ForwardDiff you have to do something: you have to make u0 be dual numbers if p is dual numbers, since after one step you’ll be differentiating the ODE state. If everything is too implicit you can’t do anything other than error, and then the user needs to learn to do weird stuff like new_u0 = eltype(p).(u0) before solving. That’s even just forward-mode AD, with reverse mode making p explicit can be required in order to write rules w.r.t. it for many of the AD systems. Then there’s some things in terms of thread safety, etc. it just turns out that if the user can tell you what p is we can make life a lot easier.
Nope, to solve it can be general. There are some constraints only when autodiff gets involved, in which case that’s documented with the AD system: