Merge pull request #269 from divital-coder/dg/vectors

SciMLOperators operator evaluation interface allowing different defining and action vectors

Author: ChrisRackauckas

.github/workflows/Tests.yml

@@ -4,6 +4,7 @@ on:
   pull_request:
     branches:
       - master
+      - 'release-'
     paths-ignore:
       - 'docs/**'
   push:
@@ -13,19 +14,26 @@ on:
       - 'docs/**'
 
 concurrency:
+  # Skip intermediate builds: always, but for the master branch.
+  # Cancel intermediate builds: always, but for the master branch.
   group: ${{ github.workflow }}-${{ github.ref }}
-  cancel-in-progress: ${{ github.ref_name != github.event.repository.default_branch || github.ref != 'refs/tags/v*' }}
+  cancel-in-progress: ${{ github.ref != 'refs/heads/master' }}
 
 jobs:
   tests:
     name: "Tests"
     strategy:
+      fail-fast: false
       matrix:
         version:
           - "1"
           - "lts"
           - "pre"
+        group:
+          - Core
+          - Downstream
     uses: "SciML/.github/.github/workflows/tests.yml@v1"
     with:
       julia-version: "${{ matrix.version }}"
+      group: "${{ matrix.group }}"
     secrets: "inherit"

docs/Project.toml

@@ -6,4 +6,4 @@ SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 [compat]
 Documenter = "0.27, 1"
 FFTW = "1.7"
-SciMLOperators = "0.2, 0.3"
+SciMLOperators = "0.4"

docs/src/index.md

@@ -32,7 +32,7 @@ Let `M`, `D`, `F` be matrix-based, diagonal-matrix-based, and function-based
 
 ```julia
 N = 4
-f = (u, p, t) -> u .* u
+f = (v, u, p, t) -> u .* v
 
 M = MatrixOperator(rand(N, N))
 D = DiagonalOperator(rand(N))
@@ -56,10 +56,11 @@ p = nothing # parameter struct
 t = 0.0     # time
 
 u = rand(N)
-v = L1(u, p, t) # == L1 * u
+v = rand(N)
+w = L1(v, u, p, t) # == L1 * v
 
-u_kron = rand(N^3)
-v_kron = L3(u_kron, p, t) # == L3 * u_kron
+v_kron = rand(N^3)
+w_kron = L3(v_kron, u, p, t) # == L3 * v_kron
 ```
 
 For mutating operator evaluations, call `cache_operator` to generate an
@@ -73,21 +74,17 @@ L2 = cache_operator(L2, u)
 L4 = cache_operator(L4, u)
 
 # allocation-free evaluation
-L2(v, u, p, t) # == mul!(v, L2, u)
-L4(v, u, p, t, α, β) # == mul!(v, L4, u, α, β)
+L2(w, v, u, p, t) # == mul!(w, L2, v)
+L4(w, v, u, p, t, α, β) # == mul!(w, L4, v, α, β)
 ```
 
-The calling signature `L(u, p, t)`, for out-of-place evaluations, is
-equivalent to `L * u`, and the in-place evaluation `L(v, u, p, t, args...)`
-is equivalent to `LinearAlgebra.mul!(v, L, u, args...)`, where the arguments
-`p, t` are passed to `L` to update its state. More details are provided
-in the operator update section below. While overloads to `Base.*`
-and `LinearAlgebra.mul!` are available, where a `SciMLOperator` behaves
-like an `AbstractMatrix`, we recommend sticking with the
-`L(u, p, t)`, `L(v, u, p, t)`, `L(v, u, p, t, α, β)` calling signatures
-as the latter internally update the operator state.
-
-The `(u, p, t)` calling signature is standardized over the `SciML`
+The calling signature `L(v, u, p, t)`, for out-of-place evaluations, is
+equivalent to `L * v`, and the in-place evaluation `L(w, v, u, p, t, args...)`
+is equivalent to `LinearAlgebra.mul!(w, L, v, args...)`, where the arguments
+`u, p, t` are passed to `L` to update its state. More details are provided
+in the operator update section below.
+
+The `(v, u, p, t)` calling signature is standardized over the `SciML`
 ecosystem and is flexible enough to support use cases such as time-evolution
 in ODEs, as well as sensitivity computation with respect to the parameter
 object `p`.

docs/src/interface.md

@@ -1,152 +1,7 @@
 # The `AbstractSciMLOperator` Interface
 
-## Formal Properties of SciMLOperators
-
-These are the formal properties that an `AbstractSciMLOperator` should obey
-for it to work in the solvers.
-
- 1. An `AbstractSciMLOperator` represents a linear or nonlinear operator, with input/output
-    being `AbstractArray`s. Specifically, a SciMLOperator, `L`, of size `(M, N)` accepts an
-    input argument `u` with leading length `N`, i.e. `size(u, 1) == N`, and returns an
-    `AbstractArray` of the same dimension with leading length `M`, i.e. `size(L * u, 1) == M`.
- 2. SciMLOperators can be applied to an `AbstractArray` via overloaded `Base.*`, or
-    the in-place `LinearAlgebra.mul!`. Additionally, operators are allowed to be time,
-    or parameter dependent. The state of a SciMLOperator can be updated by calling
-    the mutating function `update_coefficients!(L, u, p, t)` where `p` represents
-    parameters, and `t`, time.  Calling a SciMLOperator as `L(du, u, p, t)` or out-of-place
-    `L(u, p, t)` will automatically update the state of `L` before applying it to `u`.
-    `L(u, p, t)` is the same operation as `L(u, p, t) * u`.
- 3. To support the update functionality, we have lazily implemented a comprehensive operator
-    algebra. That means a user can add, subtract, scale, compose and invert SciMLOperators,
-    and the state of the resultant operator would be updated as expected upon calling
-    `L(du, u, p, t)` or `L(u, p, t)` so long as an update function is provided for the
-    component operators.
-
-## Overloaded Traits
-
-Thanks to overloads defined for evaluation methods and traits in
-`Base`, `LinearAlgebra`, the behavior of a `SciMLOperator` is
-indistinguishable from an `AbstractMatrix`. These operators can be
-passed to linear solver packages, and even to ordinary differential
-equation solvers. The list of overloads to the `AbstractMatrix`
-interface includes, but is not limited to, the following:
-
-  - `Base: size, zero, one, +, -, *, /, \, ∘, inv, adjoint, transpose, convert`
-  - `LinearAlgebra: mul!, ldiv!, lmul!, rmul!, factorize, issymmetric, ishermitian, isposdef`
-  - `SparseArrays: sparse, issparse`
-
-## Multidimensional arrays and batching
-
-SciMLOperator can also be applied to `AbstractMatrix` subtypes where
-operator-evaluation is done column-wise.
-
-```julia
-K = 10
-u_mat = rand(N, K)
-
-v_mat = F(u_mat, p, t) # == mul!(v_mat, F, u_mat)
-size(v_mat) == (N, K) # true
-```
-
-`L#` can also be applied to `AbstractArray`s that are not
-`AbstractVecOrMat`s so long as their size in the first dimension is appropriate
-for matrix-multiplication. Internally, `SciMLOperator`s reshapes an
-`N`-dimensional array to an `AbstractMatrix`, and applies the operator via
-matrix-multiplication.
-
-## Operator update
-
-This package can also be used to write time-dependent, and
-parameter-dependent operators, whose state can be updated per
-a user-defined function.
-The updates can be done in-place, i.e. by mutating the object,
-or out-of-place, i.e. in a non-mutating, `Zygote`-compatible way.
-
-For example,
-
-```julia
-u = rand(N)
-p = rand(N)
-t = rand()
-
-# out-of-place update
-mat_update_func = (A, u, p, t) -> t * (p * p')
-sca_update_func = (a, u, p, t) -> t * sum(p)
-
-M = MatrixOperator(zero(N, N); update_func = mat_update_func)
-α = ScalarOperator(zero(Float64); update_func = sca_update_func)
-
-L = α * M
-L = cache_operator(L, u)
-
-# L is initialized with zero state
-L * u == zeros(N) # true
-
-# update operator state with `(u, p, t)`
-L = update_coefficients(L, u, p, t)
-# and multiply
-L * u != zeros(N) # true
-
-# updates state and evaluates L at (u, p, t)
-L(u, p, t) != zeros(N) # true
-```
-
-The out-of-place evaluation function `L(u, p, t)` calls
-`update_coefficients` under the hood, which recursively calls
-the `update_func` for each component `SciMLOperator`.
-Therefore, the out-of-place evaluation function is equivalent to
-calling `update_coefficients` followed by `Base.*`. Notice that
-the out-of-place evaluation does not return the updated operator.
-
-On the other hand, the in-place evaluation function, `L(v, u, p, t)`,
-mutates `L`, and is equivalent to calling `update_coefficients!`
-followed by `mul!`. The in-place update behavior works the same way,
-with a few `<!>`s appended here and there. For example,
-
-```julia
-v = rand(N)
-u = rand(N)
-p = rand(N)
-t = rand()
-
-# in-place update
-_A = rand(N, N)
-_d = rand(N)
-mat_update_func! = (A, u, p, t) -> (copy!(A, _A); lmul!(t, A); nothing)
-diag_update_func! = (diag, u, p, t) -> copy!(diag, N)
-
-M = MatrixOperator(zero(N, N); update_func! = mat_update_func!)
-D = DiagonalOperator(zero(N); update_func! = diag_update_func!)
-
-L = D * M
-L = cache_operator(L, u)
-
-# L is initialized with zero state
-L * u == zeros(N) # true
-
-# update L in-place
-update_coefficients!(L, u, p, t)
-# and multiply
-mul!(v, u, p, t) != zero(N) # true
-
-# updates L in-place, and evaluates at (u, p, t)
-L(v, u, p, t) != zero(N) # true
-```
-
-The update behavior makes this package flexible enough to be used
-in `OrdinaryDiffEq`. As the parameter object `p` is often reserved
-for sensitivity computation via automatic-differentiation, a user may
-prefer to pass in state information via other arguments. For that
-reason, we allow update functions with arbitrary keyword arguments.
-
-```julia
-mat_update_func = (A, u, p, t; scale = 0.0) -> scale * (p * p')
-
-M = MatrixOperator(zero(N, N); update_func = mat_update_func,
-    accepted_kwargs = (:state,))
-
-M(u, p, t) == zeros(N) # true
-M(u, p, t; scale = 1.0) != zero(N)
+```@docs
+SciMLOperators.AbstractSciMLOperator
 ```
 
 ## Interface API Reference
@@ -219,6 +74,6 @@ update_coefficients!(γ, nothing, nothing, nothing; my_special_scaling = 7.0)
 @show γ * [2.0]
 
 # Use operator application form
-@show γ([2.0], nothing, nothing; my_special_scaling = 5.0)
+@show γ([2.0], nothing, nothing, nothing; my_special_scaling = 5.0)
 nothing # hide
 ```

docs/src/premade_operators.md

@@ -3,7 +3,7 @@
 ## Direct Operator Definitions
 
 ```@docs
-ScalarOperator.IdentityOperator
+SciMLOperators.IdentityOperator
 SciMLOperators.NullOperator
 ScalarOperator
 MatrixOperator

docs/src/sciml.md

@@ -19,22 +19,32 @@ Munthe-Kaas methods require defining operators of the form ``u' = A(u) u``.
 Thus, the operators need some form of time and state dependence, which the
 solvers can update and query when they are non-constant
 (`update_coefficients!`). Additionally, the operators need the ability to
-act like “normal” functions for equation solvers. For example, if `A(u,p,t)`
-has the same operation as `update_coefficients(A, u, p, t); A * u`, then `A`
+act like “normal” functions for equation solvers. For example, if `A(v,u,p,t)`
+has the same operation as `update_coefficients(A, u, p, t); A * v`, then `A`
 can be used in any place where a differential equation definition
-`f(u, p, t)` is used without requiring the user or solver to do any extra
-work. Thus while previous good efforts for matrix-free operators have existed
+`(u,p,t) -> A(u, u, p, t)` is used without requiring the user or solver to do any extra
+work. 
+
+Another example is state-dependent mass matrices. `M(u,p,t)*u' = f(u,p,t)`.
+When solving such an equation, the solver must understand how to "update M"
+during operations, and thus the ability to update the state of `M` is a required
+function in the interface. This is also required for the definition of Jacobians
+`J(u,p,t)` in order to be properly used with Krylov methods inside of ODE solves
+without reconstructing the matrix-free operator at each step.
+
+Thus while previous good efforts for matrix-free operators have existed
 in the Julia ecosystem, such as
 [LinearMaps.jl](https://github.com/JuliaLinearAlgebra/LinearMaps.jl), those
 operator interfaces lack these aspects to actually be fully seamless
 with downstream equation solvers. This necessitates the definition and use of
 an extended operator interface with all of these properties, hence the
 `AbstractSciMLOperator` interface.
 
-Some packages providing similar functionality are
+!!! warn
 
-  - [LinearMaps.jl](https://github.com/JuliaLinearAlgebra/LinearMaps.jl)
-  - [`DiffEqOperators.jl`](https://github.com/SciML/DiffEqOperators.jl/tree/master) (deprecated)
+    This means that LinearMaps.jl is fundamentally lacking and is incompatible
+    with many of the tools in the SciML ecosystem, except for the specific cases
+    where the matrix-free operator is a constant!
 
 ## Interoperability and extended Julia ecosystem
 

docs/src/tutorials/fftw.md

@@ -16,18 +16,18 @@ L = 2π
 
 dx = L / n
 x = range(start = -L / 2, stop = L / 2 - dx, length = n) |> Array
-u = @. sin(5x)cos(7x);
-du = @. 5cos(5x)cos(7x) - 7sin(5x)sin(7x);
+v = @. sin(5x)cos(7x);
+w = @. 5cos(5x)cos(7x) - 7sin(5x)sin(7x);
 
 k = rfftfreq(n, 2π * n / L) |> Array
 m = length(k)
 P = plan_rfft(x)
 
-fwd(u, p, t) = P * u
-bwd(u, p, t) = P \ u
+fwd(v, u, p, t) = P * v
+bwd(v, u, p, t) = P \ v
 
-fwd(du, u, p, t) = mul!(du, P, u)
-bwd(du, u, p, t) = ldiv!(du, P, u)
+fwd(w, v, u, p, t) = mul!(w, P, v)
+bwd(w, v, u, p, t) = ldiv!(w, P, v)
 
 F = FunctionOperator(fwd, x, im * k;
     T = ComplexF64, op_adjoint = bwd,
@@ -38,10 +38,10 @@ F = FunctionOperator(fwd, x, im * k;
 ik = im * DiagonalOperator(k)
 Dx = F \ ik * F
 
-Dx = cache_operator(Dx, x)
+Dx = cache_operator(Dx, v)
 
-@show ≈(Dx * u, du; atol = 1e-8)
-@show ≈(mul!(copy(u), Dx, u), du; atol = 1e-8)
+@show ≈(Dx * v, w; atol = 1e-8)
+@show ≈(mul!(copy(w), Dx, v), w; atol = 1e-8)
 ```
 
 ## Explanation
@@ -61,8 +61,8 @@ n = 256
 dx = L / n
 x = range(start = -L / 2, stop = L / 2 - dx, length = n) |> Array
 
-u = @. sin(5x)cos(7x);
-du = @. 5cos(5x)cos(7x) - 7sin(5x)sin(7x);
+v = @. sin(5x)cos(7x);
+w = @. 5cos(5x)cos(7x) - 7sin(5x)sin(7x);
 ```
 
 Now, we define the Fourier transform. Since our input is purely Real, we use the real
@@ -79,15 +79,15 @@ P = plan_rfft(x)
 
 Now we are ready to define our wrapper for the FFT object. To `FunctionOperator`, we
 pass the in-place forward application of the transform,
-`(du,u,p,t) -> mul!(du, transform, u)`, its inverse application,
-`(du,u,p,t) -> ldiv!(du, transform, u)`, as well as input and output prototype vectors.
+`(w,v,u,p,t) -> mul!(w, transform, v)`, its inverse application,
+`(w,v,u,p,t) -> ldiv!(w, transform, v)`, as well as input and output prototype vectors.
 
 ```@example fft_explanation
-fwd(u, p, t) = P * u
-bwd(u, p, t) = P \ u
+fwd(v, u, p, t) = P * v
+bwd(v, u, p, t) = P \ v
 
-fwd(du, u, p, t) = mul!(du, P, u)
-bwd(du, u, p, t) = ldiv!(du, P, u)
+fwd(w, v, u, p, t) = mul!(w, P, v)
+bwd(w, v, u, p, t) = ldiv!(w, P, v)
 F = FunctionOperator(fwd, x, im * k;
     T = ComplexF64, op_adjoint = bwd,
     op_inverse = bwd,
@@ -106,6 +106,6 @@ Dx = F \ ik * F
 
 Dx = cache_operator(Dx, x)
 
-@show ≈(Dx * u, du; atol = 1e-8)
-@show ≈(mul!(copy(u), Dx, u), du; atol = 1e-8)
+@show ≈(Dx * v, w; atol = 1e-8)
+@show ≈(mul!(copy(w), Dx, v), w; atol = 1e-8)
 ```