diff --git a/docs/pages.jl b/docs/pages.jl index 62447c9e..fde0da1a 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -2,7 +2,11 @@ pages = [ "Home" => "index.md", - "Tutorial" => "tutorial.md", - "Usage" => "usage.md", + "Tutorials" => [ + "Using the SciML Symbolic Indexing Interface" => "usage.md", + "Simple Demonstration of a Symbolic System Structure" => "simple_sii_sys.md", + "Implementing the Complete Symbolic Indexing Interface" => "implementing_sii_problemsolution.md", + ], + "Defining Solution Wrapper Fallbacks" => "solution_wrappers.md", "API" => "api.md", ] diff --git a/docs/src/tutorial.md b/docs/src/complete_sii.md similarity index 65% rename from docs/src/tutorial.md rename to docs/src/complete_sii.md index bd92b18c..e714718b 100644 --- a/docs/src/tutorial.md +++ b/docs/src/complete_sii.md @@ -1,125 +1,76 @@ -# Implementing SymbolicIndexingInterface for a type +# Implementing the Complete Symbolic Indexing Interface -Implementing the interface for a type allows it to be used by existing symbolic indexing -infrastructure. There are multiple ways to implement it, and the entire interface is -not always necessary. - -## Defining a fallback - -The simplest case is when the type contains an object that already implements the interface. -All its methods can simply be forwarded to that object. To do so, SymbolicIndexingInterface.jl -provides the [`symbolic_container`](@ref) method. For example, +This tutorial will show how to define the entire Symbolic Indexing Interface on an +`ExampleSystem`: ```julia -struct MySolutionWrapper{T<:SciMLBase.AbstractTimeseriesSolution} - sol::T - # other properties... -end - -symbolic_container(sys::MySolutionWrapper) = sys.sol -``` - -`MySolutionWrapper` wraps an `AbstractTimeseriesSolution` which already implements the interface. -Since `symbolic_container` will return the wrapped solution, all method calls such as -`is_parameter(sys::MySolutionWrapper, sym)` will be forwarded to `is_parameter(sys.sol, sym)`. - -In cases where some methods need to function differently than those of the wrapped type, they can be selectively -defined. For example, suppose `MySolutionWrapper` does not support observed quantities. The following -method can be defined (in addition to the one above): - -```julia -is_observed(sys::MySolutionWrapper, sym) = false -``` - -## Defining the interface in its entirety - -Not all the methods in the interface are required. Some only need to be implemented if a type -supports specific functionality. Consider the following struct, which needs to implement the interface: - -```julia -struct ExampleSolution +struct ExampleSystem state_index::Dict{Symbol,Int} parameter_index::Dict{Symbol,Int} independent_variable::Union{Symbol,Nothing} # mapping from observed variable to Expr to calculate its value observed::Dict{Symbol,Expr} - u::Vector{Vector{Float64}} - p::Vector{Float64} - t::Vector{Float64} end ``` +Not all the methods in the interface are required. Some only need to be implemented if a type +supports specific functionality. Consider the following struct, which needs to implement the interface: + +## Mandatory methods +### Simple Indexing Functions -### Mandatory methods +These are the simple functions which describe how to turn symbols into indices. ```julia -function SymbolicIndexingInterface.is_variable(sys::ExampleSolution, sym) +function SymbolicIndexingInterface.is_variable(sys::ExampleSystem, sym) haskey(sys.state_index, sym) end -function SymbolicIndexingInterface.variable_index(sys::ExampleSolution, sym) +function SymbolicIndexingInterface.variable_index(sys::ExampleSystem, sym) get(sys.state_index, sym, nothing) end -function SymbolicIndexingInterface.variable_symbols(sys::ExampleSolution) +function SymbolicIndexingInterface.variable_symbols(sys::ExampleSystem) collect(keys(sys.state_index)) end -function SymbolicIndexingInterface.is_parameter(sys::ExampleSolution, sym) +function SymbolicIndexingInterface.is_parameter(sys::ExampleSystem, sym) haskey(sys.parameter_index, sym) end -function SymbolicIndexingInterface.parameter_index(sys::ExampleSolution, sym) +function SymbolicIndexingInterface.parameter_index(sys::ExampleSystem, sym) get(sys.parameter_index, sym, nothing) end -function SymbolicIndexingInterface.parameter_symbols(sys::ExampleSolution) +function SymbolicIndexingInterface.parameter_symbols(sys::ExampleSystem) collect(keys(sys.parameter_index)) end -function SymbolicIndexingInterface.is_independent_variable(sys::ExampleSolution, sym) +function SymbolicIndexingInterface.is_independent_variable(sys::ExampleSystem, sym) # note we have to check separately for `nothing`, otherwise # `is_independent_variable(p, nothing)` would return `true`. sys.independent_variable !== nothing && sym === sys.independent_variable end -function SymbolicIndexingInterface.independent_variable_symbols(sys::ExampleSolution) +function SymbolicIndexingInterface.independent_variable_symbols(sys::ExampleSystem) sys.independent_variable === nothing ? [] : [sys.independent_variable] end -# this type accepts `Expr` for observed expressions involving state/parameter/observed -# variables -SymbolicIndexingInterface.is_observed(sys::ExampleSolution, sym) = sym isa Expr || sym isa Symbol && haskey(sys.observed, sym) - -function SymbolicIndexingInterface.observed(sys::ExampleSolution, sym::Expr) - if is_time_dependent(sys) - return function (u, p, t) - # compute value from `sym`, leveraging `variable_index` and - # `parameter_index` to turn symbols into indices - end - else - return function (u, p) - # compute value from `sym`, leveraging `variable_index` and - # `parameter_index` to turn symbols into indices - end - end -end - -function SymbolicIndexingInterface.is_time_dependent(sys::ExampleSolution) +function SymbolicIndexingInterface.is_time_dependent(sys::ExampleSystem) sys.independent_variable !== nothing end -SymbolicIndexingInterface.constant_structure(::ExampleSolution) = true +SymbolicIndexingInterface.constant_structure(::ExampleSystem) = true -function SymbolicIndexingInterface.all_solvable_symbols(sys::ExampleSolution) +function SymbolicIndexingInterface.all_solvable_symbols(sys::ExampleSystem) return vcat( collect(keys(sys.state_index)), collect(keys(sys.observed)), ) end -function SymbolicIndexingInterface.all_symbols(sys::ExampleSolution) +function SymbolicIndexingInterface.all_symbols(sys::ExampleSystem) return vcat( all_solvable_symbols(sys), collect(keys(sys.parameter_index)), @@ -128,18 +79,45 @@ function SymbolicIndexingInterface.all_symbols(sys::ExampleSolution) end ``` +### Observed Equation Handling + +These are for handling symbolic expressions and generating equations which are not directly +in the solution vector. + +```julia +# this type accepts `Expr` for observed expressions involving state/parameter/observed +# variables +SymbolicIndexingInterface.is_observed(sys::ExampleSystem, sym) = sym isa Expr || sym isa Symbol && haskey(sys.observed, sym) + +function SymbolicIndexingInterface.observed(sys::ExampleSystem, sym::Expr) + if is_time_dependent(sys) + return function (u, p, t) + # compute value from `sym`, leveraging `variable_index` and + # `parameter_index` to turn symbols into indices + end + else + return function (u, p) + # compute value from `sym`, leveraging `variable_index` and + # `parameter_index` to turn symbols into indices + end + end +end +``` + +### Note about constant structure + Note that the method definitions are all assuming `constant_structure(p) == true`. In case `constant_structure(p) == false`, the following methods would change: -- `constant_structure(::ExampleSolution) = false` -- `variable_index(sys::ExampleSolution, sym)` would become - `variable_index(sys::ExampleSolution, sym i)` where `i` is the time index at which +- `constant_structure(::ExampleSystem) = false` +- `variable_index(sys::ExampleSystem, sym)` would become + `variable_index(sys::ExampleSystem, sym i)` where `i` is the time index at which the index of `sym` is required. -- `variable_symbols(sys::ExampleSolution)` would become - `variable_symbols(sys::ExampleSolution, i)` where `i` is the time index at which +- `variable_symbols(sys::ExampleSystem)` would become + `variable_symbols(sys::ExampleSystem, i)` where `i` is the time index at which the variable symbols are required. -- `observed(sys::ExampleSolution, sym)` would become - `observed(sys::ExampleSolution, sym, i)` where `i` is either the time index at which +- `observed(sys::ExampleSystem, sym)` would become + `observed(sys::ExampleSystem, sym, i)` where `i` is either the time index at which the index of `sym` is required or a `Vector` of state symbols at the current time index. ## Optional methods @@ -157,15 +135,16 @@ the default implementations for `getp` and `setp` will suffice, and manually def them is not necessary. ```julia -function SymbolicIndexingInterface.parameter_values(sys::ExampleSolution) +function SymbolicIndexingInterface.parameter_values(sys::ExampleSystem) sys.p end ``` -# Implementing the `SymbolicTypeTrait` for a type +## Implementing the `SymbolicTypeTrait` for a type The `SymbolicTypeTrait` is used to identify values that can act as symbolic variables. It has three variants: + - [`NotSymbolic`](@ref) for quantities that are not symbolic. This is the default for all types. - [`ScalarSymbolic`](@ref) for quantities that are symbolic, and represent a single diff --git a/docs/src/index.md b/docs/src/index.md index 6de2852d..898f5926 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,4 +1,4 @@ -# SymbolicIndexingInterface.jl: Arrays of Arrays and Even Deeper +# SymbolicIndexingInterface.jl: Standardized Symbolic Indexing of Julia SymbolicIndexingInterface.jl is a set of interface functions for handling containers of symbolic variables. @@ -12,6 +12,19 @@ using Pkg Pkg.add("SymbolicIndexingInterface") ``` +## Introduction + +The symbolic indexing interface has 2 levels: + +1. The user level. At the user level, a modeler or engineer simply uses terms from a + domain-specific language (DSL) inside of SciML functionality and will receive the requested + values. For example, if a DSL defines a symbol `x`, then `sol[x]` returns the solution + value(s) for `x`. +2. The DSL system structure level. This is the structure which defines the symbolic indexing + for a given problem/solution. DSLs can tag a constructed problem/solution with this + object in order to endow the SciML tools with the ability to index symbolically according + to the definitions the DSL writer wants. + ## Contributing - Please refer to the @@ -79,4 +92,4 @@ file and the [project]($link_project) file. """) -``` \ No newline at end of file +``` diff --git a/docs/src/simple_sii_sys.md b/docs/src/simple_sii_sys.md new file mode 100644 index 00000000..b19e1e59 --- /dev/null +++ b/docs/src/simple_sii_sys.md @@ -0,0 +1,6 @@ +# Simple Demonstration of a Symbolic System Structure + +In this tutorial we will show how to implement a system structure type for defining the +symbolic indexing of a domain-specific language. This tutorial will show how the +`SymbolCache` type is defined to take in arrays of symbols for its independent, dependent, +and parameter variable names and uses that to define the symbolic indexing interface. diff --git a/docs/src/solution_wrappers.md b/docs/src/solution_wrappers.md new file mode 100644 index 00000000..c6599a58 --- /dev/null +++ b/docs/src/solution_wrappers.md @@ -0,0 +1,26 @@ +# Defining Solution Wrapper Fallbacks + +The simplest case is when the type contains an object that already implements the interface. +All its methods can simply be forwarded to that object. To do so, SymbolicIndexingInterface.jl +provides the [`symbolic_container`](@ref) method. For example, + +```julia +struct MySolutionWrapper{T<:SciMLBase.AbstractTimeseriesSolution} + sol::T + # other properties... +end + +symbolic_container(sys::MySolutionWrapper) = sys.sol +``` + +`MySolutionWrapper` wraps an `AbstractTimeseriesSolution` which already implements the interface. +Since `symbolic_container` will return the wrapped solution, all method calls such as +`is_parameter(sys::MySolutionWrapper, sym)` will be forwarded to `is_parameter(sys.sol, sym)`. + +In cases where some methods need to function differently than those of the wrapped type, they can be selectively +defined. For example, suppose `MySolutionWrapper` does not support observed quantities. The following +method can be defined (in addition to the one above): + +```julia +is_observed(sys::MySolutionWrapper, sym) = false +``` diff --git a/docs/src/usage.md b/docs/src/usage.md index 8b77048d..db91e73f 100644 --- a/docs/src/usage.md +++ b/docs/src/usage.md @@ -1,6 +1,23 @@ -# Using the SymbolicIndexingInterface +# Using the SciML Symbolic Indexing Interface + +This tutorial will cover ways to use the symbolic indexing interface for types that +implement it. SciML's core types (problems, solutions, and iterator (integrator) types) +all support this symbolic indexing interface which allows for domain-specific interfaces +(such as ModelingToolkit, Catalyst, etc.) to seamlessly blend their symbolic languages with +the types obtained from SciML. Other tutorials will focus on how users can make use of the +interface for their own DSL, this tutorial will simply focus on what the user experience +looks like for DSL which have set it up. + +We recommend any DSL implementing the symbolic indexing interface to link to this tutorial +as a full description of the functionality. + +!!! note + While this tutorial focuses on demonstrating the symbolic indexing interface for ODEs, + note that the same functionality works across all of the other problem types, such as + optimization problems, nonlinear problems, nonlinear solutions, etc. + +## Symbolic Indexing of Differential Equation Solutions -This tutorial will cover ways to use the interface for types that implement it. Consider the following example: ```@example Usage @@ -20,11 +37,13 @@ sys = structural_simplify(sys) ``` The system has 4 state variables, 3 parameters, and one observed variable: + ```@example Usage ModelingToolkit.observed(sys) ``` Solving the system, + ```@example Usage u0 = [D(x) => 2.0, x => 1.0, @@ -41,12 +60,14 @@ sol = solve(prob, Tsit5()) ``` We can obtain the timeseries of any time-dependent variable using `getindex` + ```@example Usage sol[x] ``` This also works for arrays or tuples of variables, including observed quantities and independent variables, for interpolating solutions, and plotting: + ```@example Usage sol[[x, y]] ``` @@ -74,6 +95,7 @@ plot(sol, idxs=x) If necessary, `Symbol`s can be used to refer to variables. This is only valid for symbolic variables for which [`hasname`](@ref) returns `true`. The `Symbol` used must match the one returned by [`getname`](@ref) for the variable. + ```@example Usage hasname(x) ``` @@ -89,6 +111,7 @@ sol[(:x, :w)] Note how when indexing with an array, the returned type is a `Vector{Array{Float64}}`, and when using a `Tuple`, the returned type is `Vector{Tuple{Float64, Float64}}`. To obtain the value of all state variables, we can use the shorthand: + ```@example Usage sol[solvedvariables] # equivalent to sol[variable_symbols(sol)] ``` @@ -99,8 +122,16 @@ output, the following shorthand is used: sol[allvariables] # equivalent to sol[all_variable_symbols(sol)] ``` +## Parameter Indexing: Getting and Setting Parameter Values + Parameters cannot be obtained using this syntax, and instead require using [`getp`](@ref) and [`setp`](@ref). +!!! note + The reason why parameters use a separate syntax is to be able to ensure type stability + of the `sol[x]` indexing. Without separating the parameter indexing, the return type of + symbolic indexing could be anything a parameter can be, which is general is not the same + type as state variables! + ```@example Usage σ_getter = getp(sys, σ) σ_getter(sol) # can also pass `prob` @@ -110,7 +141,7 @@ Note that this also supports arrays/tuples of parameter symbols: ```@example Usage σ_ρ_getter = getp(sys, (σ, ρ)) -σ_ρ_getter(sol) +σ_ρ_getter(sol) ``` Now suppose the system has to be solved with a different value of the parameter `β`. @@ -145,3 +176,10 @@ To set the entire parameter vector at once, [`parameter_values`](@ref) can be us parameter_values(prob) .= [29.0, 11.0, 2.5] parameter_values(prob) ``` + +!!! note + These getters and setters generate high-performance functions for the specific chosen + symbols or collection of symbols. Caching the getter/setter function and reusing it + on other problem/solution instances can be the key to achieving good performance. Note + that this caching is allowed only when the symbolic system is unchanged (it's fine for + the numerical values to have changed, but not the underlying equations).