From d897b7b162ce9083ed92079eaff97c0fc8d02655 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 17 Dec 2023 09:49:07 -0500 Subject: [PATCH 1/6] Start the definition of concrete number and container interfaces --- docs/pages.jl | 5 +- docs/src/interfaces/Array_and_Number.md | 102 ++++++++++++++++++++++++ 2 files changed, 106 insertions(+), 1 deletion(-) create mode 100644 docs/src/interfaces/Array_and_Number.md diff --git a/docs/pages.jl b/docs/pages.jl index 45dc683ec..a223f1360 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -2,8 +2,11 @@ pages = [ "Home" => "index.md", - "Interfaces" => Any["interfaces/Problems.md", + "Interfaces" => Any[ + "interfaces/Array_and_Number.md", + "interfaces/Problems.md", "interfaces/SciMLFunctions.md", + "interfaces/Array_and_Number.md", "interfaces/Algorithms.md", "interfaces/Solutions.md", "interfaces/Init_Solve.md", diff --git a/docs/src/interfaces/Array_and_Number.md b/docs/src/interfaces/Array_and_Number.md new file mode 100644 index 000000000..6f4a1411b --- /dev/null +++ b/docs/src/interfaces/Array_and_Number.md @@ -0,0 +1,102 @@ +# [SciML Container (Array) and Number Interfaces](@id arrayandnumber) + +We live in a society, and therefore there are rules. In this tutorial we outline +the rules which are required on container and number types which are allowable +in SciML tools. + +!!! warn + + In general as of 2023, strict adherence to this interface is an early work-in-progress. + If anything does not conform to the documented interface, please open an issue. + +!!! note + + There are many types which can work with a specific solver that do satisfy this + interface. Many times as part of prototyping you may want to side-step the + high level interface checks in order to simply test whether a new type is working. + To do this, set `interface_checks = false` as a keyword argument to `init`/`solve` + to bypass any of the internal interface checks. This means you will no longer get + a nice high-level error message and instead it will attempt to use the type + without restrictions. Note that not every problem/solver has not implemented this + new keyword argument as of 2023. + +## Note About Wrapped Solvers + +Due to limitations of wrapped solvers, any solver that is a wrapped solver from an existing C/Fortran +code is inherently limited to `Float64` and `Vector{Float64}` for its operations. This includes packages +like Sundials.jl, LSODA.jl, DASKR.jl, MINPACK.jl, and many more. This is fundamental to these solvers +and it is not expected that they will allow the full set of SciML types in the future. If more abstract +number/container definitions are required, then these are not the appropriate solvers to use. + +## SciML Number Types + +The number types are the types used to define the dependent variables (i.e. `u0`) and the +independent variables (`t` or `tspan`). These two types can be different, and can have +different restrictions depending on the type of solver which is employed. The following +rules for a Number type are held in general: + +* Number types can be used in SciML directly or in containers. If a problem defines a value like `u0` + using a Number type, the out-of-place form must be used for the problem definition. +* `x::T + y::T = z::T` +* `x::T * y::T = z::T` +* `oneunit(x::T)::T` +* `one(x::T) * oneunit(x::T) = z::T` +* `t::T2 * x::T + y::T = z::T` for `T2` a time type and `T` the dependent variable type (this includes the + `muladd` equivalent form). + +Additionally, the following rules apply to subsets of uses: + +### Adaptive Number Types + +* `x::T / y::T = z::T` +* Default choices of norms can assume `sqrt(x::T)::T` exists. If `internalnorm` is overriden then this + may not be required (for example, changing the norm to inf-norm). +* `x::T ^ y::T = z::T` + +### Time Types (Independent Variables) + +* If a solver is time adaptive, the time type must be a floating point number. `Rational` is only allowed + for non-adaptive solves. + +## SciML Container (Array) Types + +Container types are types which hold number types. They can be used to define objects like the state vector +(`u0`) of a problem. The following operations are required in a container type to be used with SciML +solvers: + +* Broadcast is defined [according to the Julia broadcast interface](https://docs.julialang.org/en/v1/manual/interfaces/#man-interfaces-broadcasting). +* The container type correctly defines [interface overloads to satisfy the ArrayInterface.jl specification](https://docs.sciml.ai/ArrayInterface/stable/). +* `ArrayInterface.zeromatrix(x::T)::T2` defines a compatible matrix type (see below) +* `eltype(x::T)::T2` is a compatible Number type. +* `x::T .+ y::T = z::T` (i.e. broadcast similar is defined to be type-presurving) +* Indexing is only required if `ArrayInterface.fast_scalar_indexing(x::T)==true`. If true, + scalar indexing `x[i]` is assumed to be defined and run through all variables. + +!!! note + + "`eltype(x::T)::T2` is a compatible Number type" excludes `Array{Array{T}}` types of types. However, recursive + vectors can conformed to the interface with zero overhead using tools from RecursiveArrayTools.jl such as + `VectorOfArray(x)`. Since this greatly simplifies the interfaces and the ability to check for correctness, + doing this wrapping is highly recommended and there are no plans to relax this requirement. + +Additionally, the following rules apply to subsets of uses: + +### SciML Mutable Array Types + +* `similar(x::T)::T` +* `zero(x::T)::T` +* `z::T .= x::T .+ y::T` is defined +* `z::T .= x::T .* y::T` is defined +* `z::T .= t::T2 .* x::T` where `T2` is the time type (a Number) and `T` is the container type. + +### SciML Matrix (Operator) Type + +Note that the matrix type may not match the type of the initial container `u0`. An example is `ComponetMatrix` +as the matrix structure corresponding to a `ComponentArray`. However, the following actions are assumed +to hold on the resulting matrix type: + +* `solve(LinearProblem(A::T,b::T2),linsolve)` must be defined for a solver to work on a given SciML matrix + type `T2`. +* If the matrix is an operator, i.e. a lazy construct, it should conform to the + [SciMLOperators](https://docs.sciml.ai/SciMLOperators/stable/) interface. +* If not a SciMLOperator, `diagind(W::T)` should be defined and `@view(A[idxs])=@view(A[idxs]) + λ::T` \ No newline at end of file From cd91cfcd5ce79e3dc88844d35029a834c0759211 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 18 Dec 2023 11:29:08 -0500 Subject: [PATCH 2/6] Update docs/src/interfaces/Array_and_Number.md Co-authored-by: Hendrik Ranocha --- docs/src/interfaces/Array_and_Number.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/interfaces/Array_and_Number.md b/docs/src/interfaces/Array_and_Number.md index 6f4a1411b..3d7f48c48 100644 --- a/docs/src/interfaces/Array_and_Number.md +++ b/docs/src/interfaces/Array_and_Number.md @@ -91,7 +91,7 @@ Additionally, the following rules apply to subsets of uses: ### SciML Matrix (Operator) Type -Note that the matrix type may not match the type of the initial container `u0`. An example is `ComponetMatrix` +Note that the matrix type may not match the type of the initial container `u0`. An example is `ComponentMatrix` as the matrix structure corresponding to a `ComponentArray`. However, the following actions are assumed to hold on the resulting matrix type: From 06d5c5141b77795c98b8e094b6197069e036e618 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 18 Dec 2023 11:29:15 -0500 Subject: [PATCH 3/6] Update docs/src/interfaces/Array_and_Number.md Co-authored-by: Hendrik Ranocha --- docs/src/interfaces/Array_and_Number.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/interfaces/Array_and_Number.md b/docs/src/interfaces/Array_and_Number.md index 3d7f48c48..28b3acdaa 100644 --- a/docs/src/interfaces/Array_and_Number.md +++ b/docs/src/interfaces/Array_and_Number.md @@ -17,7 +17,7 @@ in SciML tools. To do this, set `interface_checks = false` as a keyword argument to `init`/`solve` to bypass any of the internal interface checks. This means you will no longer get a nice high-level error message and instead it will attempt to use the type - without restrictions. Note that not every problem/solver has not implemented this + without restrictions. Note that not every problem/solver has implemented this new keyword argument as of 2023. ## Note About Wrapped Solvers From 25eb199bb1d8d201af2c2c8ad218e08bc6e94fa6 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 18 Dec 2023 11:35:48 -0500 Subject: [PATCH 4/6] Update docs/src/interfaces/Array_and_Number.md --- docs/src/interfaces/Array_and_Number.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/interfaces/Array_and_Number.md b/docs/src/interfaces/Array_and_Number.md index 28b3acdaa..6c3303299 100644 --- a/docs/src/interfaces/Array_and_Number.md +++ b/docs/src/interfaces/Array_and_Number.md @@ -88,6 +88,7 @@ Additionally, the following rules apply to subsets of uses: * `z::T .= x::T .+ y::T` is defined * `z::T .= x::T .* y::T` is defined * `z::T .= t::T2 .* x::T` where `T2` is the time type (a Number) and `T` is the container type. +* (Optional) `Base.resize!(x,i)` is required for `resize!(integrator,i)` to be supported. ### SciML Matrix (Operator) Type From 08395fd0972df50333370d18231edbc474d90ee3 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 18 Dec 2023 11:36:06 -0500 Subject: [PATCH 5/6] Update docs/pages.jl --- docs/pages.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/pages.jl b/docs/pages.jl index a223f1360..afdecf4c6 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -6,7 +6,6 @@ pages = [ "interfaces/Array_and_Number.md", "interfaces/Problems.md", "interfaces/SciMLFunctions.md", - "interfaces/Array_and_Number.md", "interfaces/Algorithms.md", "interfaces/Solutions.md", "interfaces/Init_Solve.md", From 4d6df15b05aabd0c9ddee8e630707626db6989e5 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 18 Dec 2023 14:35:00 -0500 Subject: [PATCH 6/6] Update docs/src/interfaces/Array_and_Number.md Co-authored-by: Hendrik Ranocha --- docs/src/interfaces/Array_and_Number.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/interfaces/Array_and_Number.md b/docs/src/interfaces/Array_and_Number.md index 6c3303299..4239169b5 100644 --- a/docs/src/interfaces/Array_and_Number.md +++ b/docs/src/interfaces/Array_and_Number.md @@ -49,7 +49,7 @@ Additionally, the following rules apply to subsets of uses: ### Adaptive Number Types * `x::T / y::T = z::T` -* Default choices of norms can assume `sqrt(x::T)::T` exists. If `internalnorm` is overriden then this +* Default choices of norms can assume `sqrt(x::T)::T` exists. If `internalnorm` is overridden then this may not be required (for example, changing the norm to inf-norm). * `x::T ^ y::T = z::T`