SciML / DiffEqBase.jl

The lightweight Base library for shared types and functionality for defining differential equation and scientific machine learning (SciML) problems
Other
319 stars 116 forks source link

create API for tstop aliasing system #1077

Open isaacsas opened 2 months ago

isaacsas commented 2 months ago

Per discussions @ChrisRackauckas and I had around https://github.com/SciML/JumpProcesses.jl/issues/442 it would be nice to be able to have users tell us we can alias their tstops vector in JumpProcesses so we don't need to allocate a new vector for each call to solve.

ChrisRackauckas commented 2 months ago

My plan is to create a keyword argument alias which then has a struct ODEAliases <: AbstractAliasSpecifier etc., that holds booleans/nothing for each thing that can be aliased. nothing for default behavior, true/false for forced aliasing, de-aliasing. This can then be rolled out for every problem type uniformly, and include booleans to cover the keyword arguments.

ChrisRackauckas commented 2 months ago

@jClugstor does this specification make sense?

ChrisRackauckas commented 2 months ago

You'll find:

alias_u0: allows the solver to alias the initial condition array that is contained in the problem struct. Defaults to false.

in DiffEq

alias_A::Bool: Whether to alias the matrix A or use a copy by default. When true, algorithms like LU-factorization can be faster by reusing the memory via lu!, but care must be taken as the original input will be modified. Default is true if the algorithm is known not to modify A, otherwise is false. alias_b::Bool: Whether to alias the matrix b or use a copy by default. When true, algorithms can write and change b upon usage. Care must be taken as the original input will be modified. Default is true if the algorithm is known not to modify b, otherwise false.

in LinearSolve. Etc. All of that should get folded into just these more general structs, and the old way should get deprecated to just setting values of the struct.

jClugstor commented 2 months ago

I think I get it. alias should be a keyword of ODEProblem, LinearProblem, etc. that would be of type e.g. ODEAliases and LinearAliases respectively. And each concrete subtype of AbstractAliasSpecifier should hold boolean values that correspond to the parts of their corresponding problems that can be aliased.

ChrisRackauckas commented 2 months ago

Yes exactly.

jClugstor commented 2 months ago

Should the types ODEAliases and LinearAliases etc. be defined in SciMLBase?

ChrisRackauckas commented 2 months ago

yes

jClugstor commented 2 months ago

I'm thinking that the problem types will need a new field to hold the AbstractAliasSpecifiers, does that sound right?

ChrisRackauckas commented 1 month ago

I don't think so, they will just be a new common kwarg to solve.

jClugstor commented 2 weeks ago

These four PRs should add everything that's needed to use this for ODEs and Linear Problems.

Basically just like you said, I added an AbstractAliasSpecifier in SciMLBase, and then subtypes can be defined for every problem type. In DiffEqBase all I needed to do was add alias to the acceptable kwargs list.

Oh and I'm not sure I got all of the version bumps right, but I think it's just that LinearSolve and OrdinaryDiffEqCore need to be using SciMLBase@2.58 and DiffEqBase@6.160 for it all to work.

https://github.com/SciML/DiffEqBase.jl/pull/1099 https://github.com/SciML/SciMLBase.jl/pull/830 https://github.com/SciML/LinearSolve.jl/pull/553 https://github.com/SciML/OrdinaryDiffEq.jl/pull/2503

ChrisRackauckas commented 2 weeks ago

The interface needs a few more things. We should define the interface in AbstractAliasSpecifier to have:

Any other things you see @isaacsas @ranocha @avik-pal ?

ranocha commented 2 weeks ago

There is also a kwarg alias_du0, e.g, https://github.com/SciML/OrdinaryDiffEq.jl/pull/2503/files#diff-28a5b9e0abdb761fcac81be080bd7bb47c49f41c9d6b18df041d51881f6b1f41R73

ranocha commented 2 weeks ago
  • Deprecate the existing alias_* keyword arguments for the general alias one. As a non-breaking change, the concrete AbstractAliasSpecifiers should match the keyword arguments that exist by default.

I think we have to take this one step further - if users just use the old API, everything has to work as documented. See, e.g., https://github.com/SciML/OrdinaryDiffEq.jl/pull/2503/files#r1818084222

isaacsas commented 2 weeks ago

Having an alias option around the tstops vector would also be nice, and was the original motivation for this issue. I was just going to add it to JumpProcesses, but held off to use the aliasing system.

ChrisRackauckas commented 2 weeks ago

I think we have to take this one step further - if users just use the old API, everything has to work as documented. See, e.g., https://github.com/SciML/OrdinaryDiffEq.jl/pull/2503/files#r1818084222

Yes, that's what I meant by the deprecation path. It shouldn't change, just be extended with a deprecation warning.

There is also a kwarg alias_du0

Having an alias option around the tstops vector would also be nice, and was the original motivation for this issue. I was just going to add it to JumpProcesses, but held off to use the aliasing system.

Yes, and these two are examples of kwargs to specific instances of AbstractAliasSpecifier, specifically ODEAliasSpecifier should document and add those two.

jClugstor commented 2 weeks ago

So the full list for ODEAliases for example would be


struct ODEAliases <: AbstractAliasSpecifier
    alias::Union{Bool,Nothing}
    alias_p::Union{Bool,Nothing}
    alias_f::Union{Bool,Nothing}

    alias_u0::Union{Bool,Nothing}
    alias_du0::Union{Bool,Nothing}
    alias_tstops::Union{Bool,Nothing}
end
``` ?
ChrisRackauckas commented 2 weeks ago
struct ODEAliases <: AbstractAliasSpecifier
    alias_p::Union{Bool,Nothing}
    alias_f::Union{Bool,Nothing}
    alias_u0::Union{Bool,Nothing}
    alias_du0::Union{Bool,Nothing}
    alias_tstops::Union{Bool,Nothing}
end

alias is just a high level kwarg for simplicity.

jClugstor commented 1 week ago

A couple of questions: Are there any solvers where alias_p and alias_f should actually do anything at this point? I haven't found any solvers where these are keywords. I assume alias_p means alias the p parameter array, what should alias_f do?

So far the only other variable I can find that explicitly mentions any aliasing is alias_u0 in NonlinearSolve. Are there any other solvers that already have an aliasing option that I need to change to make use of the aliasing API?

ChrisRackauckas commented 1 week ago

Are there any solvers where alias_p and alias_f should actually do anything at this point? I haven't found any solvers where these are keywords.

No, those are ones to add.

I assume alias_p means alias the p parameter array, what should alias_f do?

Yes, and f is for whether it should double prob.f, i.e. integrator.f = deepcopy(prob.f), to decouple caches from the user's function.

jClugstor commented 1 week ago

Got it. I should be able to implement these by changing things in __init right? I shouldn't have to change anything in the actual solver code?

ChrisRackauckas commented 1 week ago

Yes exactly. You'll find that there's places in the __init definitions where things just do p = prob.p, and that builds the aliasing. It would need to have options to copy instead.