Consider breaking syntax changes

[ChrisRackauckas]

With Julia v1.0 being released, we have a last big attempt at breaking syntax changes. The question is, should we do any?

Two that would be possible are:

1. All mutation goes to the front. f(du,t,u) instead of f(t,u,du) (or f(du,u,t)).
2. Parameters are always explicit, i.e. f(t,u,p,du) as standard. Then there would be no need for the ParameterizedFunction wrappers. Instead, problems could hold parameters, and maybe they could be saved each step with an option?

---

Edit:

We are doing f(mutate, dependent variables, p/integrator, independent variables). Compatible on master:

All are updated on master branches. Check marks now mean tags.

• OrdinaryDiffEq
• StochasticDiffEq
• DelayDiffEq
• Sundials
• DiffEqBayes
• DiffEqDocs
• DiffEqNoiseProcess
• BoundaryValueDiffEq
• DiffEqParamEstim
• SteadyStateDiffEq
• DiffEqSensitivity
• DiffEqDiffTools
• DiffEqBenchmarks
• DiffEqCallbacks
• DiffEqBase (tuples, analytical solution, p in prob)
• DiffEqTutorials
• OrdinaryDiffEqExtendedTests
• ODEInterfaceDiffEq
• DiffEqOperators
• DiffEqPhysics
• ParameterizedFunctions
• DifferentialEquations
• DiffEqJump
• MultiScaleArrays
• BridgeDiffEq
• GeometricIntegratorsDiffEq
• DiffEqDevTools
• DiffEqProblemLibrary
• DASKR
• DiffEqPDEBase
• DiffEqApproxFun
• DiffEqDevDocs
• DiffEqMonteCarlo
• DiffEqBiological
• MATLABDiffEq
• ODE
• PyDSTool
• DiffEqFinancial
• DASSL
• JuliaCon2017
• FiniteElementDiffEq
• DiffEqUncertainty
• LSODA

And dependent libraries:

• RandomMatrices
• DynamicalSystemsBase
• LTISystems
• VehicleModels
• ChaosTools
• Mads
• SwitchTimeOpt
• QuantumOptics
• DynamicMovementPrimatives
• BioEnergeticFoodWebs
• Sims

[Datseris]

By 2, you mean that it would be possible to only use Functors and never care about ParameterizedFunction wrappers? if so, yes please!

For the 1, YEEEEEEES pleeeeeeeeease! I actually have written DynamicalSystems with the "proposed change" format initially, just because it was so weird for me. But I don't mind changing it back.

For 1 however, please have it as f(du, u, t) instead f(du, t, u). The second just feels "off" for me at least. That is my humble vote.

EDIT: wait, now I understood point 2. Why make them explicit always? Just have the optional definition f(du, u, t [, p])? Then everyone is happy.

[ranocha]

1. Personally, I would also prefer f(du, u, t). I think this change should be made since the mutated argument is mostly in the first place.

2. If there is no additional overhead, I have no objections against explicit parameters. Maybe it could also be f(du,u,t,p)? How would you distinguish a mutating function without parameters f(du,u,t) from an out-of-place function with parameters f(t,u,p) or f(u,t,p)?

I haven't used ParameterizedFunctions extensively since I'm more interested in semidiscretisations, building up the ODE function object differently. Thus, my vote shouldn't count that much.

[ysimillides]

1. I personally prefer mutated arguments towards the end, but believe consistency is more important (whether that will be with the wrapped libraries / other Julia packages, I can't say.).

[BeastyBlacksmith]

1. Personally I prefer the order old values then new values semantically. But since it seems to have become a julian standard to put mutating arguments first, let it be.
2. I strongly support making parameters explicit. From a users perspective it would be convenient to have both t and p optional, but I don't see a way to achieve this with positional arguments, personally I would even consider making them keyword arguments, if there is no performance penalty.

[Datseris]

but I don't see a way to achieve this with positional arguments, personally I would even consider making them keyword arguments, if there is no performance penalty

I definitely do not support keyword arguments for such a basic component of the equations of motion.

Fortunately it is very easy to make them positional default arguments, since t and p are always of fundamental different type. t<:Real while p is a container but never a bit type.

[ChrisRackauckas]

while p is a container but never a bit type.

Not necessarily. And those two aren't mutually exclusive either.

[Datseris]

Well, unless you want to always write 2 definitions in your source code in the handling of equations of motion, one that has p and one that has p[1]...? How else could you do it?

[ChrisRackauckas]

If you have one parameter you can use a number and it'll work just fine in the current way.

[Datseris]

Sounds unnecessary complication for me. It's not like you can gain something performance wise by allowing possibility of both number and container?

[ChrisRackauckas]

You can make it a tuple of a single number, sure. That's still a bitstype though, so I'm not sure why you'd bring that up.

[aytekinar]

I would like to vote Yes for both. Mutating the first argument is what I am used to seeing, e.g., in Base.A_mul_B!. As for the difference between (t,u) and (u,t), I would vote for (t,u) as I am used to seeing it in textbooks such as "Nonlinear Systems" or "Nonlinear Control" (Khalil).

[Datseris]

so I'm not sure why you'd bring that up.

Alright, my expression was bad. What I wanted to point out is that you could always have p and t as optional (not keyword) arguments and correctly differentiatiate between them via multiple dispatch because t<:Real but if p is a tuple that is not the case. Even for 1 element.

[Datseris]

@aytekinar

I would vote for (t,u)

That would not allow the possibility of f(du, u [, t, p]) however, with both t and p optional but not keyword. Unless of course one started adding new methods all the time that do (du, u) -> (du, 0, u) via anonymous.

[aytekinar]

@Datseris, I can see your point. Maybe the (u,t[, ]) version is better after the mutating argument.

The intuition in textbooks, as far as I can see, is that they write the f(t,x1,x2,...) as f(t,x) easily. That extends to the situation when an exogenous signal u is provided (f(t,x,u)). But now that you mentioned the optional arguments, I think the other version might be the way to go.

[ChrisRackauckas]

When making them optional, remember what that would actually do to user codes. They always have to be able to accept the optional stuff, either by

function f(du,u,args...)
   ...
end

or

function f(du,u;kwargs...)
  ....
end

otherwise it would error if the integrator internally calls f(du,u,t,p)/f(du,u,t=t,p=p) since those methods have to be defined even if the variables aren't used. Querying before hand to find out all of the signatures that the user allows and modifying all of the function calls based on their definition would be fairly difficult.

The interesting thing about a standard signature which involves splatting is we could always add weird stuff in the future though... but I'm not sure how much that'd actually do anything.

[aytekinar]

That being said, from @ChrisRackauckas's perspective, that would not make much difference at all, either. If there is an assumption that u is AbstractVector, even for scalar systems, then the type information would solve the issue easily (here, t would be assumed to be Real).

[Datseris]

@ChrisRackauckas

I do not like splatting, I would really hate to write my equations of motion as function f(du,u,args...) and inside the body do t = args[1]. In addition, with my suggestion, users will only have to define methods that they need to use. And I don't see why they should define more than one...?

If they need to use both t and p they define f(du, u, t, p), if they need only t then f(du, u, t) and if only p then f(du, u, p).

Why would the user care about whether a method or not is optional. They just use what they need, and not define extra.

[ChrisRackauckas]

How would I know which one to call? I have u, t, p, and du output. What do I call in the integrator?

if is_two_param
  f(du,u)
elseif is_three_param
  f(du,u,t)
elseif is_four_param
  f(du,u,t,p)
end

? (that's actually not too difficult to get everywhere because it could be built via a closure, but that's quite prone to bugs). This is because if you call a function which the user hasn't defined, it'll error. Even then, there's the question of how to distinguish between those cases and still allow a non-mutating form (which is quite important!).

I'm not entirely sure how I could find out all of this information without making some of it explicit type information. This then ties over to SciML/DiffEqBase.jl#52 where if we can have DiffEqFunction types, then the user could do:

function f(du,u)
  ...
end
af = AutonomousFunction(f)
prob = ODEProblem(af,u0,tspan)

when isn't less verbose but now allows the solver to know what form it must be in order to call it correctly.

[Datseris]

How would I know which one to call?

I did not consider that. Then I am all up for having explicit parameter f(du, u, t, p). For the use case without parameters it's just 2 more characters to write, no big deal. Yet having explicit parameters would make everything so simpler internally, which I fully concur with.

Haven;t thought about it from that perspective.

[Datseris]

Oh man! If the f(du, u, t, p) becomes a reality I wont have to do the crazy setfield!() on Structs used as equations of motion, done here: https://juliadynamics.github.io/DynamicalSystems.jl/latest/chaos/orbitdiagram/#producing-orbit-diagrams-for-continuous-flows
to produce orbit diagrams for continuous systems. That would be really cool! Then users would just pass in parameter containers.

[ChrisRackauckas]

So it looks like everybody is for explicit parameters and mutation in front?

f(du,u,t,p) or f(du,t,u,p)? Upvote for the first, downvote for the latter.

Making the extra arguments optional will not be possible. What will be possible will be that the two forms, in place an out of place, which is whether du is in front or doesn't exist.

Oh man! If the f(du, u, t, p) becomes a reality I wont have to do the crazy setfield!() on Structs used as equations of motion

Yes, this will get rid of a whole layer of closures that we commonly use.

[JackDevine]

Then there would be no need for the ParameterizedFunction wrappers. Instead, problems could hold parameters, and maybe they could be saved each step with an option?

Why do the parameters need to be saved at each step if they are not changing?

[ChrisRackauckas]

Why do the parameters need to be saved at each step if they are not changing?

Because you can make them time-varying. The DEDataArray is for this purpose of saving the parameters. It seems to come up in a lot of people's use-cases that they care about time-varying algebraic quantities, though arguably the SavingCallback and FunctionCallingCallback do this better now. If this was ever to become an option, it would definitely be false by default though, so don't worry about it.

[Datseris]

f(du,u,t,p) or f(du,t,u,p)? Upvote for the first, downvote for the latter.

Alright, I am willing to bribe people to vote for the first. Alternative facts methods incoming now

[JackDevine]

Because you can make them time-varying.

That actually makes a lot of sense if the parameters are a known function of time. Now that I understand that, I think that the proposed syntax makes a lot of sense. Unfortunately, I don't think that there is any correct ordering for t and u. Most of the things that I read (in the world of physics) use f(u,t), but a lot of mathematics textbooks use f(t,u). I we ignore how people write these things on paper, I think that f(du,u,t,p) is better because it puts du and u next to each other as well as t and p next to each other. Those pairings happen to agree with my mental model of the right hand side of an ODE.

[ChrisRackauckas]

So let's look at what this would look like in full. I'll use method 1 because that seems to be winning so far. We would have:

f(u,t,p) # out-of-place ODE/SDE
f(du,u,t,p) # in-place ODE/SDE
f(du,u,t,p) # out-of-place DAE
f(res,du,u,t,p) # in-place DAE
f(u,h,t,p) # out-of-place DDE
f(du,u,h,t,p) # in-place DDE
f(du,u,t,x,y,p) # proposed PDE syntax
rate(t,u,integrator) # for jumps, add the integrator
integrator.p # parameters directly in the integrator
ODEProblem(f,u0,tspan,p=nothing) # parameters passed to problem
f(v,x,t) # DynamicalODEProblem out-of-place syntax
f(dv,v,x,t) # DynamicalODEProblem in-place syntax
f(du,u,W,t) # RODE in-place syntax

Additionally, changes would include:

h(out,t;idxs=nothing,deriv=0) # HistoryFunction args move to keyword arguments

and moving the Jacobian/etc. specification to DiffEqFunction. Anything else we should break at the same time?

Question: would it be preferred if we did this at the same time as the 1.0 transition, i.e. the version which drops v0.6 is the same one which then introduces these changes? I am not sure if this change can give deprecations BTW, so it will be majorly breaking.

[JackDevine]

Looks good, I am not too fussed about whether it is synced with 1.0 transitions, just whatever turns out to be the easiest.

Right now, I can't think of anything else that should be broken.

[aytekinar]

Is it possible to change

rate(t,u,integrator) # for jumps, add the integrator

to

rate(u,t,integrator) # for jumps, add the integrator

to keep the consistency? Or, am I missing something? Most likely I would not be using it at all, but when I checked the code, my OCD got triggered ☺️

[ChrisRackauckas]

You're right. I missed that one.

One nice thing about this change is that it would be easy to add a switch so that instead of p you can get the entire integrator. Do you think that's too detailed to be standard (though integrator.p giving parameters makes it a superset of functionality)?

Even if it's not standard, having a switch for this would allow all sorts of weird algorithms and cache re-use (I'll be using it in the internal sensitivity calculation algorithms for example).

[tkf]

(I'm an outsider here so I'm not sure if this is more than a noise. But anyway...)

I've never used PDE solver but does this:

f(du,u,t,x,y,p) # proposed PDE syntax

mean you would have f(du,u,t,x,y,z,p) for "3-dimensional" problem and f(du,u,t,xyz...,p) in general? If that's the case, it makes me wonder if it is more sensible if you put "time and spatial" axes in the end like this f(du,u,p,t,x,y,...).

But vararg in the middle is no problem in Julia so this is not a bit deal anyway.

PS: ...or even f(p, [du, ]...) so that all mutating arguments are the first ones?

[JackDevine]

Looking through that commit I see two possible problems,

1. faq.md: line 70 - Function has not been updated to new syntax.
2. faq.md: line 239 - pf needs to be func

I can complete this list and make a PR when I get more time later on today (probably not until I get home in about 5 hours).

[ChrisRackauckas]

Alright, it's settled. After talking to some people individually, and some in chats, it seems there's pretty much a consensus for f(mutate, dependent variables, p/integrator, independent variables), or at least those against it were not strongly against it anyways. And with the upvotes to ODEProblem(pf,u0,tspan,p=nothing), it looks like we have a plan.

Thanks everyone who pitched in. I will announce this right away and get going on it. From now on, master branch is allowed to have this syntax. The OP now will track which packages upgraded their master. A massive amount of version bounds on DiffEqBase will go through, and then nothing will upgrade until all packages go through METADATA. Every package will get a major version bump (except Sundials which just did...?), and DifferentialEquations.jl will be the last to update after all other go through. Please feel free to pitch in and help out. We will especially need eyes to keep going over the docs, tutorials, and benchmarks to make sure all of those examples are up to date. I don't think the notebooks need to be re-ran, but they should be made runnable.

[Datseris]

Sounds good. I will be contributing indirectly by making sure that my packages run after the change and what unexpected behavior/syntax errors I get, which is quite a test of usage I have to admit.

[JackDevine]

Sounds fantastic. Sorry that I didn't go through the rest of the docs thoroughly like I said I would, I got very held up yesterday. I will definitely have a read through at some stage when I get some time. Thank you so much for all of your hard work, I really appreciate it.

[aytekinar]

Personally, I would like to thank you for your efforts --- not only in making the package really useful and well-documented but also in asking for users' preferences in doing so.

We will especially need eyes to keep going over the docs, tutorials, and benchmarks to make sure all of those examples are up to date. I don't think the notebooks need to be re-ran, but they should be made runnable.

This I will try to help with... definitely! :)

Cheers!

[ChrisRackauckas]

Everything is upper bounded on DiffEqBase, DiffEqDiffTools, and OrdinaryDiffEq. This won't make anything come into effect until the whole ecosystem is tagged. Tags are in for DiffEqBase, DiffEqDiffTools, and OrdinaryDiffEq as well. Dependents have been added to the list of necessary PRs.

[Datseris]

what is the link for the updated docs? The page you have linked in a previous comment doesn't exist anymore.

[ChrisRackauckas]

Docs are live.

http://docs.juliadiffeq.org/latest/tutorials/ode_example.html

[Datseris]

Hm, maybe it is a good time to think whether to have the in-place version as the first seen example?

(or mention StaticArrays otherwise?)

EDIT: First-seen examples are quite important as it is the first thing new users see and therefore the first thing they use. It is also often the case that people do not go further than the first seen example, if they can serve their case.

However, if they don't use static arrays and this leads to bad performance, somebody really noob in Julia could say "ooooh what is this its not faster that pythoooon uga uga uga" :D :D :D

[ChrisRackauckas]

Hm, maybe it is a good time to think whether to have the in-place version as the first seen example?

That is a good question. On one hand, it's what most people should be using. On the other hand, I have found that it can be very confusing to people who are not "programmers", and so many stick with the out of place form simply out of familiarity. It hurts me when thinking about speed, but the ease-of-use should not be underrated.

[Datseris]

On the other hand, I have found that it can be very confusing to people who are not "programmers"

That may be true, but it is also worth to say that (at least in Julia) the skills required to understand in-place operations are as high (if not even lower) as the skills required to just handle Vector data in general. And I think it is clear that the latter skills are mandatory to use DiffEq and actually do something with it.

[ChrisRackauckas]

I am not sure that's true. This is a good example:

https://github.com/PoisotLab/BioEnergeticFoodWebs.jl

@tpoisot's library is a great use of DiffEq, but it currently uses the out-of-place format on Vectors. This is because they couldn't figure out the cache resizing for ODEs of changing size:

PoisotLab/BioEnergeticFoodWebs.jl@3c46723

I'll go in and fix this which a PR that does the cache resizing, but it is something that's not simple. Though you can also say that resizing diffeqs is something specialized itself.

[Datseris]

Though you can also say that resizing diffeqs is something specialized itself.

Can't imagine something more advanced/complicated to be perfectly honest, so I don't think this is an argument for the first-seen example :D

[ChrisRackauckas]

I missed condition(t,u,integrator) before, which is now flipped to condition(u,t,integrator) for consistency.

[ChrisRackauckas]

The core libraries are all updated. Now I am waiting on @pkofod to tag an NLsolve.jl update and then these will all get tagged.

[Datseris]

Is there any example of using the parameters while them not being an Array? I really cannot think of how to do it and allow for other types besides arrays... @ChrisRackauckas you said that the type of p doesn't matter, so I would appreciate some feedback here.

I have a function where I want to solve problems for different parameters, and the user to choose which parameter to be varied (with what values). If p is an AbstractArray, the user can give the parameter values and the parameter index on the array structure.

But what if p is not an Array? Then this interface does not hold anymore :(

[ChrisRackauckas]

But what if p is not an Array? Then this interface does not hold anymore :(

Then require an Array. There's nothing about DiffEq's engine that needs it to be an array, but you can require it to be an array for your uses. But that was always true about how parameters could be used.

Is there any example of using the parameters while them not being an Array? I really cannot think of how to do it and allow for other types besides arrays... @ChrisRackauckas you said that the type of p doesn't matter, so I would appreciate some feedback here.

It's a nice way without using global constants to have a bunch of cache variables in your function. Also, named tuples are a great use of parameters because then they are named (e.g. p.bioav) and I have a whole library (private for now) making use of that. Things like tuples are faster for access too, so localizing a bunch of constants by making it a parameter tuple is a good idea. I guess the key is that making parameters accessible means that the pattern:

const a = 0.89
f(u,p,t) = a*u
ODEProblem(f,u0,tspan)

can become

a = 0.89
f(u,a,t) = a*u
ODEProblem(f,u0,tspan,a)

I really disliked the amount of const that a "good" DiffEq code could have, and from what I find is that almost every problem I've ever created in the tests could actually make use of this and it would clean up my code quite a bit (since I const for different types all the type, when in reality it's just a new problem where the parameter is a different type!). In this sense, I find that the vast majority of codes is actually using parameters in some form, so allowing it to be whatever object to make localizing easier is a win.

But again, you can restrict on your end. I bet you $20 that you've already made a lot of assumptions about what can be solved that don't necessarily hold, but that's fine.

[ChrisRackauckas]

Okay, tags are starting to go through. The release post is drafted:

https://github.com/JuliaDiffEq/juliadiffeq.github.io/blob/40/_posts/2018-1-24-Parameters.md

One last thing to discuss.

10. For any Dynamical ODE function `f(t,u1,u2,dui)`, it's now `f(dui,u1,u2,p,t)`
11. For any Dynamical ODE function `f(t,u1,u2)`, it's now `f(u1,u2,p,t)`
12. For any second order ODE function `f(t,u,du,out)`, it's now `f(out,u,du,p,t)`
13. For any second order ODE function `f(t,u,du)`, it's now `f(u,du,p,t)`

Does this ordering make sense for second order ODEs? Note that it has to match DynamicalODE functions.

@Datseris

[ChrisRackauckas]

The blog post is live:

http://juliadiffeq.org/2018/01/24/Parameters.html

Everything is passing tests. Tagging is going. After the tags go through, then I will focus on getting PRs to dependent packages.

[ChrisRackauckas]

Every JuliaDiffEq library has a tag PR in METADATA, and the core libraries have tags already. DifferentialEquations.jl won't get a tag until all of these are in, which once that's true I will wipe my packages, download DiffEq, checkout master, and make sure that runs clean. Then that will tag and it'll all be live. There's upper bounds on the dependent packages still which will need PRs too.

[ChrisRackauckas]

DifferentialEquations.jl is now tagged. PRs to dependent libraries coming soon.

[ChrisRackauckas]

The rest of the PRs @Datseris said he'd do, so this is done.

[Datseris]

DynamicalSystemsBase is already done. ChaosTools is finishing in a couple of minutes and then QuantumOptics