Skip to content

Commit 2eaf240

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #194 from vpuri3/show
overload `Base.show`, add docs
2 parents da6da41 + dd821bc commit 2eaf240

File tree

18 files changed

+771
-95
lines changed

18 files changed

+771
-95
lines changed

‎README.md‎

Lines changed: 314 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,6 @@
1-
# SciMLOperators.jl
1+
# `SciMLOperators.jl`
2+
3+
*Unified operator interface for `SciML.ai` and beyond*
24

35
[![Join the chat at https://julialang.zulipchat.com #sciml-bridged](https://img.shields.io/static/v1?label=Zulip&message=chat&color=9558b2&labelColor=389826)](https://julialang.zulipchat.com/#narrow/stream/279055-sciml-bridged)
46
[![Global Docs](https://img.shields.io/badge/docs-SciML-blue.svg)](https://docs.sciml.ai/SciMLOperators/stable)
@@ -8,3 +10,314 @@
810

911
[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet)](https://github.com/SciML/ColPrac)
1012
[![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle)
13+
14+
`SciMLOperators` is a package for managing linear, nonlinear,
15+
time-dependent, and parameter dependent operators acting on vectors,
16+
(or column-vectors of matrices). We provide wrappers for matrix-free
17+
operators, fast tensor-product evaluations, pre-cached mutating
18+
evaluations, as well as `Zygote`-compatible non-mutating evaluations.
19+
20+
The lazily implemented operator algebra allows the user to update the
21+
operator state by passing in an update function that accepts arbirary
22+
parameter objects. Further, our operators behave like `AbstractMatrix` types
23+
thanks to overloads defined for methods in `Base`, and `LinearAlgebra`.
24+
25+
Therefore, an `AbstractSciMLOperator` can be passed to `LinearSolve.jl`,
26+
or `NonlinearSolve.jl` as a linear/nonlinear operator, or to
27+
`OrdinaryDiffEq.jl` as an `ODEFunction`. Examples of usage within the
28+
`SciML` ecosystem are provided in the documentation.
29+
30+
## Installation
31+
`SciMLOperators.jl` is a registerd package and can be installed via
32+
33+
```
34+
julia> import Pkg
35+
julia> Pkg.add("SciMLOperators")
36+
```
37+
38+
## Examples
39+
40+
Let `M`, `D`, `F` be matrix, diagonal matrix, and function-based
41+
`SciMLOperators` respectively.
42+
43+
```julia
44+
N = 4
45+
f = (u, p, t) -> u .* u
46+
47+
M = MatrixOperator(rand(N, N))
48+
D = DiagonalOperator(rand(N))
49+
F = FunctionOperator(f, zeros(N), zeros(N))
50+
```
51+
52+
Then, the following codes just work.
53+
54+
```julia
55+
L1 = 2M + 3F + LinearAlgebra.I + rand(N, N)
56+
L2 = D * F * M'
57+
L3 = kron(M, D, F)
58+
L4 = M \ D
59+
L5 = [M; D]' * [M F; F D] * [F; D]
60+
```
61+
62+
Each `L#` can be applied to `AbstractVector`s of appropriate sizes:
63+
64+
```julia
65+
p = nothing # parameter struct
66+
t = 0.0 # time
67+
68+
u = rand(N)
69+
v = L1(u, p, t) # == L1 * u
70+
71+
u_kron = rand(N ^ 3)
72+
v_kron = L3(u_kron, p, t) # == L3 * u_kron
73+
```
74+
75+
For mutating operator evaluations, call `cache_operator` to generate
76+
in-place cache so the operation is nonallocating.
77+
78+
```julia
79+
α, β = rand(2)
80+
81+
# allocate cache
82+
L2 = cache_operator(L2, u)
83+
L4 = cache_operator(L4, u)
84+
85+
# allocation-free evaluation
86+
L2(v, u, p, t) # == mul!(v, L2, u)
87+
L4(v, u, p, t, α, β) # == mul!(v, L4, u, α, β)
88+
```
89+
90+
The calling signature `L(u, p, t)`, for out-of-place evaluations is
91+
equivalent to `L * u`, and the in-place evaluation `L(v, u, p, t, args...)`
92+
is equivalent to `LinearAlgebra.mul!(v, L, u, args...)`, where the arguments
93+
`p, t` are passed to `L` to update its state. More details are provided
94+
in the operator update section below. While overloads to `Base.*`
95+
and `LinearAlgebra.mul!` are available, where a `SciMLOperator` behaves
96+
like an `AbstractMatrix`, we recommend sticking with the
97+
`L(u, p, t)`, `L(v, u, p, t)`, `L(v, u, p, t, α, β)` calling signatures
98+
as the latter internally update the operator state.
99+
100+
The `(u, p, t)` calling signature is standardized over the `SciML`
101+
ecosystem and is flexible enough to support use cases such as time-evolution
102+
in ODEs, as well as sensitivity computation with respect to the parameter
103+
object `p`.
104+
105+
Thanks to overloads defined for evaluation methods and traits in
106+
`Base`, `LinearAlgebra`, the behaviour of a `SciMLOperator` is
107+
indistinguishable from an `AbstractMatrix`. These operators can be
108+
passed to linear solver packages, and even to ordinary differential
109+
equation solvers. The list of overloads to the `AbstractMatrix`
110+
interface include, but are not limited, the following:
111+
112+
- `Base: size, zero, one, +, -, *, /, \, ∘, inv, adjoint, transpose, convert`
113+
- `LinearAlgebra: mul!, ldiv!, lmul!, rmul!, factorize, issymmetric, ishermitian, isposdef`
114+
- `SparseArrays: sparse, issparse`
115+
116+
## Multidimension arrays and batching
117+
118+
SciMLOperator can also be applied to `AbstractMatrix` subtypes where
119+
operator-evaluation is done column-wise.
120+
121+
```julia
122+
K = 10
123+
u_mat = rand(N, K)
124+
125+
v_mat = F(u_mat, p, t) # == mul!(v_mat, F, u_mat)
126+
size(v_mat) == (N, K) # true
127+
```
128+
129+
`L#` can also be applied to `AbstractArray`s that are not
130+
`AbstractVecOrMat`s so long as their size in the first dimension is appropriate
131+
for matrix-multiplication. Internally, `SciMLOperator`s reshapes an
132+
`N`-dimensional array to an `AbstractMatrix`, and applies the operator via
133+
matrix-multiplication.
134+
135+
## Operator update
136+
137+
This package can also be used to write time-dependent, and
138+
parameter-dependent operators, whose state can be updated per
139+
a user-defined function.
140+
The updates can be done in-place, i.e. by mutating the object,
141+
or out-of-place, i.e. in a non-mutating, `Zygote`-compatible way.
142+
143+
For example,
144+
145+
```julia
146+
u = rand(N)
147+
p = rand(N)
148+
t = rand()
149+
150+
# out-of-place update
151+
mat_update_func = (A, u, p, t) -> t * (p * p')
152+
sca_update_func = (a, u, p, t) -> t * sum(p)
153+
154+
M = MatrixOperator(zero(N, N); update_func = mat_update_func)
155+
α = ScalarOperator(zero(Float64); update_func = sca_update_func)
156+
157+
L = α * M
158+
L = cache_operator(L, u)
159+
160+
# L is initialized with zero state
161+
L * u == zeros(N) # true
162+
163+
# update operator state with `(u, p, t)`
164+
L = update_coefficients(L, u, p, t)
165+
# and multiply
166+
L * u != zeros(N) # true
167+
168+
# updates state and evaluates L at (u, p, t)
169+
L(u, p, t) != zeros(N) # true
170+
```
171+
172+
The out-of-place evaluation function `L(u, p, t)` calls
173+
`update_coefficients` under the hood, which recursively calls
174+
the `update_func` for each component `SciMLOperator`.
175+
Therefore the out-of-place evaluation function is equivalent to
176+
calling `update_coefficients` followed by `Base.*`. Notice that
177+
the out-of-place evaluation does not return the updated operator.
178+
179+
On the other hand,, the in-place evaluation function, `L(v, u, p, t)`,
180+
mutates `L`, and is equivalent to calling `update_coefficients!`
181+
followed by `mul!`. The in-place update behaviour works the same way
182+
with a few `<!>`s appended here and there. For example,
183+
184+
```julia
185+
v = rand(N)
186+
u = rand(N)
187+
p = rand(N)
188+
t = rand()
189+
190+
# in-place update
191+
_A = rand(N, N)
192+
_d = rand(N)
193+
mat_update_func! = (A, u, p, t) -> (copy!(A, _A); lmul!(t, A); nothing)
194+
diag_update_func! = (diag, u, p, t) -> copy!(diag, N)
195+
196+
M = MatrixOperator(zero(N, N); update_func! = mat_update_func!)
197+
D = DiagonalOperator(zero(N); update_func! = diag_update_func!)
198+
199+
L = D * M
200+
L = cache_operator(L, u)
201+
202+
# L is initialized with zero state
203+
L * u == zeros(N) # true
204+
205+
# update L in-place
206+
update_coefficients!(L, u, p, t)
207+
# and multiply
208+
mul!(v, u, p, t) != zero(N) # true
209+
210+
# updates L in-place, and evaluates at (u, p, t)
211+
L(v, u, p, t) != zero(N) # true
212+
```
213+
214+
The update behaviour makes this package flexible enough to be used
215+
in `OrdianryDiffEq`. As the parameter object `p` is often reserved
216+
for sensitivy computation via automatic-differentiation, a user may
217+
prefer to pass in state information via other arguments. For that
218+
reason, we allow for update functions with arbitrary keyword arguments.
219+
220+
```julia
221+
mat_update_func = (A, u, p, t; scale = 0.0) -> scale * (p * p')
222+
223+
M = MatrixOperator(zero(N, N); update_func = mat_update_func,
224+
accepted_kwargs = (:state,))
225+
226+
M(u, p, t) == zeros(N) # true
227+
M(u, p, t; scale = 1.0) != zero(N)
228+
```
229+
230+
## Features
231+
232+
* Matrix-free operators with `FunctionOperator`
233+
* Fast tensor product evaluation
234+
* Mutating, nonmutating update behaviour (Zygote compatible)
235+
* `InvertibleOperator` - pair fwd, bwd operators
236+
* Lazy algebra: addition, subtraction, multiplication, inverse, adjoint
237+
* Pre-caching methods for in-place evaluations
238+
239+
## Why `SciMLOperators`?
240+
241+
Many functions, from linear solvers to differential equations, require
242+
the use of matrix-free operators in order to achieve maximum performance in
243+
many scenarios. `SciMLOperators.jl` defines the abstract interface for how
244+
operators in the SciML ecosystem are supposed to be defined. It gives the
245+
common set of functions and traits which solvers can rely on for properly
246+
performing their tasks. Along with that, `SciMLOperators.jl` provides
247+
definitions for the basic standard operators which are used in building
248+
blocks for most tasks, both simplifying the use of operators while also
249+
demonstrating to users how such operators can be built and used in practice.
250+
251+
`SciMLOperators.jl` has the design that is required in order to be used in
252+
all scenarios of equation solvers. For example, Magnus integrators for
253+
differential equations require defining an operator ``u' = A(t) u``, while
254+
Munthe-Kaas methods require defining operators of the form ``u' = A(u) u``.
255+
Thus the operators need some form of time and state dependence which the
256+
solvers can update and query when they are non-constant
257+
(`update_coefficients!`). Additionally, the operators need the ability to
258+
act like "normal" functions for equation solvers. For example, if `A(u,p,t)`
259+
has the same operation as `update_coefficients(A, u, p, t); A * u`, then `A`
260+
can be used in any place where a differential equation definition
261+
`f(u, p, t)` is used without requring the user or solver to do any extra
262+
work. Thus while previous good efforts for matrix-free operators have existed
263+
in the Julia ecosystem, such as
264+
[LinearMaps.jl](https://github.com/JuliaLinearAlgebra/LinearMaps.jl), those
265+
operator interfaces lack these aspects in order to actually be fully seamless
266+
with downstream equation solvers. This necessitates the definition and use of
267+
an extended operator interface with all of these properties, hence the
268+
`AbstractSciMLOperator` interface.
269+
270+
Some packages providing similar functionality are
271+
* [LinearMaps.jl](https://github.com/JuliaLinearAlgebra/LinearMaps.jl)
272+
* [`DiffEqOperators.jl`](https://github.com/SciML/DiffEqOperators.jl/tree/master) (deprecated)
273+
274+
## Interoperability and extended Julia ecosystem
275+
276+
`SciMLOperator.jl` overloads the `AbstractMatrix` interface for
277+
`AbstractSciMLOperator`s, allowing seamless compatibility with
278+
linear solves, and nonlinear solvers. Further, due to the update functionality,
279+
`AbstractSciMLOperator`s can represent an `ODEFunction` in `OrdinaryDiffEq.jl`,
280+
and downstream packages. See tutorials for example of usage with
281+
`OrdinaryDiffEq.jl`, `LinearSolve.jl`, `NonlinearSolve.jl`.
282+
283+
Further, the nonmutating update functionality allows gradient propogation
284+
through `AbstractSciMLOperator`s, and is compatible with
285+
automatic-differentiation libraries like
286+
[`Zygote.jl`](https://github.com/SciML/DiffEqOperators.jl/tree/master).
287+
An example of `Zygote.jl` usage with
288+
[`Lux.jl`](https://github.com/LuxDL/Lux.jl) is also provided in the tutorials.
289+
290+
Please make an issue [here](https://github.com/SciML/SciMLOperators.jl/issues)
291+
if you come across an unexpected issue while using `SciMLOperators`.
292+
293+
We provide below a list of packages that make use of `SciMLOperators`.
294+
If you are using `SciMLOperators` in your work, feel free to create a PR
295+
and add your package to this list.
296+
297+
* [`SciML.ai`](https://sciml.ai/) ecosystem: `SciMLOperators` is compatible with, and utilized by every `SciML` package.
298+
* [`CalculustJL`](https://github.com/CalculustJL) packages use `SciMLOperators` to define matrix-free vector-calculus operators for solving partial differential equations.
299+
* [`CalculustCore.jl`](https://github.com/CalculustJL/CalculustCore.jl)
300+
* [`FourierSpaces.jl`](https://github.com/CalculustJL/FourierSpaces.jl)
301+
* [`NodalPolynomialSpaces.jl`](https://github.com/CalculustJL/NodalPolynomialSpaces.jl)
302+
* `SparseDiffTools.jl`
303+
304+
## Roadmap
305+
306+
- [ ] [Complete integration with `SciML` ecosystem](https://github.com/SciML/SciMLOperators.jl/issues/142)
307+
- [ ] [Block-matrices](https://github.com/SciML/SciMLOperators.jl/issues/161)
308+
- [x] [Benchmark and speed-up tensorbproduct evaluations](https://github.com/SciML/SciMLOperators.jl/issues/58)
309+
- [ ] [Fast tensor-sum (`kronsum`) evaluation](https://github.com/SciML/SciMLOperators.jl/issues/53)
310+
- [ ] [Fully flesh out operator array algebra](https://github.com/SciML/SciMLOperators.jl/issues/62)
311+
- [ ] [Operator fusion/matrix chain multiplication at constant `(u, p, t)`-slices](https://github.com/SciML/SciMLOperators.jl/issues/51)
312+
313+
## Contributing
314+
315+
- Please refer to the
316+
[SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://github.com/SciML/ColPrac/blob/master/README.md)
317+
for guidance on PRs, issues, and other matters relating to contributing to SciML.
318+
- There are a few community forums:
319+
- The #diffeq-bridged and #sciml-bridged channels in the
320+
[Julia Slack](https://julialang.org/slack/)
321+
- [JuliaDiffEq](https://gitter.im/JuliaDiffEq/Lobby) on Gitter
322+
- On the Julia Discourse forums (look for the [modelingtoolkit tag](https://discourse.julialang.org/tag/modelingtoolkit)
323+
- See also [SciML Community page](https://sciml.ai/community/)

‎docs/Project.toml‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,7 @@
11
[deps]
22
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
3-
SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
43
SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
54

65
[compat]
76
Documenter = "0.27"
8-
SciMLBase = "1.74"
97
SciMLOperators = "0.2"

‎docs/make.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ include("pages.jl")
77

88
makedocs(
99
sitename="SciMLOperators.jl",
10-
authors="Chris Rackauckas, Alex Jones, Vedant Puri",
10+
authors="Vedant Puri, Alex Jones, Chris Rackauckas",
1111
modules=[SciMLOperators],
1212
clean=true,doctest=false,
1313
format = Documenter.HTML(analytics = "UA-90474609-3",

‎docs/pages.jl‎

Lines changed: 8 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,11 @@
11
pages = [
22
"Home" => "index.md",
3-
"interface.md",
4-
"premade_operators.md",
5-
"tutorials/fftw.md"
3+
"Interface" => "interface.md",
4+
"Premade Operators" => "premade_operators.md",
5+
"API" => "api.md",
6+
"FFT Tutorial" => "tutorials/fftw.md",
7+
# "tutorials/linear.md",
8+
# "tutorials/nonlin.md",
9+
# "tutorials/ode.md",
10+
# "tutorials/lux.md",
611
]

‎docs/src/api.md‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
2+
```@autodocs
3+
Modules = [SciMLOperators]
4+
Order = [:function, :type]
5+
```

‎docs/src/assets/Project.toml‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
[deps]
2+
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
3+
SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
4+
5+
[compat]
6+
Documenter = "0.27"
7+
SciMLOperators = "0.2"

0 commit comments

Comments
 (0)