From b0148d84d2471e7323b848de94862d91e70f11e6 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 13 Mar 2024 16:31:28 +0530 Subject: [PATCH 1/2] docs: update FAQ section with relevant MTKv9 questions --- docs/src/basics/FAQ.md | 124 ++++++++++++++++++++++++++++++++++------- 1 file changed, 105 insertions(+), 19 deletions(-) diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md index 9d12b619fb..d2d1661269 100644 --- a/docs/src/basics/FAQ.md +++ b/docs/src/basics/FAQ.md @@ -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? @@ -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. From 28ff9c086f01a47a4b21973a8f8aa6706534d14f Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 13 Mar 2024 16:46:00 +0530 Subject: [PATCH 2/2] docs: minor fix to Events.md --- docs/src/basics/Events.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 61d59863e2..33ce84d31e 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -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 ```