Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

docs: update FAQ section with relevant MTKv9 questions #2542

Merged
merged 2 commits into from
Mar 13, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/src/basics/Events.md
Original file line number Diff line number Diff line change
Expand Up @@ -164,7 +164,7 @@ documentation. In affect functions, we have that
function affect!(integ, u, p, ctx)
# integ.t is the current time
# integ.u[u.v] is the value of the unknown `v` above
# integ.p[p.q] is the value of the parameter `q` above
# integ.ps[p.q] is the value of the parameter `q` above
end
```

Expand Down
124 changes: 105 additions & 19 deletions docs/src/basics/FAQ.md
Original file line number Diff line number Diff line change
@@ -1,25 +1,97 @@
# Frequently Asked Questions

## Why are my parameters some obscure object?

In ModelingToolkit.jl version 9, the parameter vector was replaced with a custom
`MTKParameters` object, whose internals are intentionally undocumented and subject
to change without a breaking release. This enables us to efficiently store and generate
code for parameters of multiple types. To obtain parameter values use
[SymbolicIndexingInterface.jl](https://github.com/SciML/SymbolicIndexingInterface.jl/) or
[SciMLStructures.jl](https://github.com/SciML/SciMLStructures.jl/). For example:

```julia
prob.ps[lorenz.β] # obtains the value of parameter `β`. Note the `.ps` instead of `.p`
getβ = getp(prob, lorenz.β) # returns a function that can fetch the value of `β`
getβ(sol) # can be used on any object that is based off of the same system
getβ(prob)
```

Indexes into the `MTKParameters` object take the form of `ParameterIndex` objects, which
are similarly undocumented. Following is the list of behaviors that should be relied on for
`MTKParameters`:

- It implements the SciMLStructures interface.
- It can be queried for parameters using functions returned from
`SymbolicIndexingInterface.getp`.
- `getindex(::MTKParameters, ::ParameterIndex)` can be used to obtain the value of a
parameter with the given index.
- `setindex!(::MTKParameters, value, ::ParameterIndex)` can be used to set the value of a
parameter with the given index.
- `parameter_values(sys, sym)` will return a `ParameterIndex` object if `sys` has been
`complete`d (through `structural_simplify`, `complete` or `@mtkbuild`).
- `copy(::MTKParameters)` is defined and duplicates the parameter object, including the
memory used by the underlying buffers.

Any other behavior of `MTKParameters` (other `getindex`/`setindex!` methods, etc.) is an
undocumented internal and should not be relied upon.

## How do I use non-numeric/array-valued parameters?

In ModelingToolkit.jl version 9, parameters are required to have a `symtype` matching
the type of their values. For example, this will error during problem construction:

```julia
@parameters p = [1, 2, 3]
```

Since by default parameters have a `symtype` of `Real` (which is interpreted as `Float64`)
but the default value given to it is a `Vector{Int}`. For array-valued parameters, use the
following syntax:

```julia
@parameters p[1:n, 1:m]::T # `T` is the `eltype` of the parameter array
@parameters p::T # `T` is the type of the array
```

The former approach is preferred, since the size of the array is known. If the array is not
a `Base.Array` or the size is not known during model construction, the second syntax is
required.

The same principle applies to any parameter type that is not `Float64`.

```julia
@parameters p1::Int # integer-valued
@parameters p2::Bool # boolean-valued
@parameters p3::MyCustomStructType # non-numeric
@parameters p4::ComponentArray{...} # non-standard array
```

## Getting the index for a symbol

Since **ordering of symbols is not guaranteed after symbolic transformations**,
one should normally refer to values by their name. For example, `sol[lorenz.x]`
from the solution. But what if you need to get the index? The following helper
function will do the trick:
Ordering of symbols is not guaranteed after symbolic transformations, and parameters
are now stored in a custom `MTKParameters` object instead of a vector. Thus, values
should be referred to by their name. For example `sol[lorenz.x]`. To obtain the index,
use the following functions from
[SymbolicIndexingInterface.jl](https://github.com/SciML/SymbolicIndexingInterface.jl/):

```julia
indexof(sym, syms) = findfirst(isequal(sym), syms)
indexof(σ, parameters(sys))
variable_index(sys, sym)
parameter_index(sys, sym)
```

Note that while the variable index will be an integer, the parameter index is a struct of
type `ParameterIndex` whose internals should not be relied upon.

## Transforming value maps to arrays

ModelingToolkit.jl allows (and recommends) input maps like `[x => 2.0, y => 3.0]`
because symbol ordering is not guaranteed. However, what if you want to get the
lowered array? You can use the internal function `varmap_to_vars`. For example:
lowered array? You can use the internal function `varmap_to_vars` for unknowns.
and the `MTKParameters` constructor for parameters. For example:

```julia
pnew = varmap_to_vars([β => 3.0, c => 10.0, γ => 2.0], parameters(sys))
unew = varmap_to_vars([x => 1.0, y => 2.0, z => 3.0], unknowns(sys))
pnew = ModelingToolkit.MTKParameters(sys, [β => 3.0, c => 10.0, γ => 2.0])
```

## How do I handle `if` statements in my symbolic forms?
Expand Down Expand Up @@ -63,37 +135,51 @@ end
Since `ODEProblem` on a MTK `sys` will have to generate code, this will be slower than
caching the generated code, and will require automatic differentiation to go through the
code generation process itself. All of this is unnecessary. Instead, generate the problem
once outside the loss function, and remake the prob inside the loss function:
once outside the loss function, and update the parameter values inside the loss function:

```julia
prob = ODEProblem(sys, [], [p1 => p[1], p2 => p[2]])
function loss(p)
remake(prob, p = ...)
# update parameters
sol = solve(prob, Tsit5())
sum(abs2, sol)
end
```

Now, one has to be careful with `remake` to ensure that the parameters are in the right
order. One can use the previously mentioned indexing functionality to generate index
maps for reordering `p` like:
If a subset of the parameters are optimized, `setp` from SymbolicIndexingInterface.jl
should be used to generate an efficient function for setting parameter values. For example:

```julia
p = @parameters x y z
idxs = ModelingToolkit.varmap_to_vars([p[1] => 1, p[2] => 2, p[3] => 3], p)
p[idxs]
using SymbolicIndexingInterface

prob = ODEProblem(sys, [], [p1 => p[1], p2 => p[2]])
setter! = setp(sys, [p1, p2])
function loss(p)
setter!(prob, p)
sol = solve(prob, Tsit5())
sum(abs2, sol)
end
```

Using this, the fixed index map can be used in the loss function. This would look like:
[SciMLStructures.jl](https://github.com/SciML/SciMLStructures.jl/) can be leveraged to
obtain all the parameters for optimization using the `Tunable` portion. By default, all
numeric or numeric array parameters are marked as tunable, unless explicitly marked as
`tunable = false` in the variable metadata.

```julia
using SciMLStructures: replace!, Tunable

prob = ODEProblem(sys, [], [p1 => p[1], p2 => p[2]])
idxs = Int.(ModelingToolkit.varmap_to_vars([p1 => 1, p2 => 2], p))
function loss(p)
remake(prob, p = p[idxs])
replace!(Tunable(), prob.p, p)
sol = solve(prob, Tsit5())
sum(abs2, sol)
end

p, replace, alias = SciMLStructures.canonicalize(Tunable(), prob.p)
# p is an `AbstractVector` which can be optimized
# if `alias == true`, then `p` aliases the memory used by `prob.p`, so
# changes to the array will be reflected in parameter values
```

# ERROR: ArgumentError: SymbolicUtils.BasicSymbolic{Real}[xˍt(t)] are missing from the variable map.
Expand Down
Loading