Merge branch 'master' into careful-event-affects

Author: BenChung

Project.toml

@@ -1,13 +1,14 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <[email protected]>", "Chris Rackauckas <[email protected]> and contributors"]
-version = "9.46.1"
+version = "9.47.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -61,21 +62,24 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
+HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 
 [extensions]
 MTKBifurcationKitExt = "BifurcationKit"
 MTKChainRulesCoreExt = "ChainRulesCore"
 MTKDeepDiffsExt = "DeepDiffs"
+MTKHomotopyContinuationExt = "HomotopyContinuation"
 MTKLabelledArraysExt = "LabelledArrays"
 
 [compat]
 AbstractTrees = "0.3, 0.4"
 ArrayInterface = "6, 7"
-BifurcationKit = "0.3"
+BifurcationKit = "0.4"
 BlockArrays = "1.1"
 ChainRulesCore = "1"
 Combinatorics = "1"
+CommonSolve = "0.2.4"
 Compat = "3.42, 4"
 ConstructionBase = "1"
 DataInterpolations = "6.4"
@@ -97,6 +101,7 @@ ForwardDiff = "0.10.3"
 FunctionWrappers = "1.1"
 FunctionWrappersWrappers = "0.1"
 Graphs = "1.5.2"
+HomotopyContinuation = "2.11"
 InteractiveUtils = "1"
 JuliaFormatter = "1.0.47"
 JumpProcesses = "9.13.1"
@@ -126,7 +131,7 @@ SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
 SymbolicIndexingInterface = "0.3.31"
 SymbolicUtils = "3.7"
-Symbolics = "6.14"
+Symbolics = "6.15.2"
 URIs = "1"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"

docs/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
@@ -23,7 +24,8 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 BenchmarkTools = "1.3"
-BifurcationKit = "0.3"
+BifurcationKit = "0.4"
+DataInterpolations = "6.5"
 DifferentialEquations = "7.6"
 Distributions = "0.25"
 Documenter = "1"

docs/pages.jl

@@ -12,7 +12,8 @@ pages = [
         "tutorials/parameter_identifiability.md",
         "tutorials/bifurcation_diagram_computation.md",
         "tutorials/SampledData.md",
-        "tutorials/domain_connections.md"],
+        "tutorials/domain_connections.md",
+        "tutorials/callable_params.md"],
     "Examples" => Any[
         "Basic Examples" => Any["examples/higher_order.md",
             "examples/spring_mass.md",

docs/src/basics/MTKLanguage.md

@@ -73,7 +73,7 @@ end
         v_array(t)[1:N, 1:M]
         v_for_defaults(t)
     end
-    @extend ModelB(; p1)
+    @extend ModelB(p1 = 1)
     @components begin
         model_a = ModelA(; k_array)
         model_array_a = [ModelA(; k = i) for i in 1:N]
@@ -149,14 +149,18 @@ julia> ModelingToolkit.getdefault(model_c1.v)
 
 #### `@extend` begin block
 
-  - Partial systems can be extended in a higher system as `@extend PartialSystem(; kwargs)`.
-  - Keyword arguments pf partial system in the `@extend` definition are added as the keyword arguments of the base system.
-  - Note that in above example, `p1` is promoted as an argument of `ModelC`. Users can set the value of `p1`. However, as `p2` isn't listed in the model definition, its initial guess can't be specified while creating an instance of `ModelC`.
+Partial systems can be extended in a higher system in two ways:
 
-```julia
-julia> @mtkbuild model_c2 = ModelC(; p1 = 2.0)
+  - `@extend PartialSystem(var1 = value1)`
+    
+      + This is the recommended way of extending a base system.
+      + The default values for the arguments of the base system can be declared in the `@extend` statement.
+      + Note that all keyword arguments of the base system are added as the keyword arguments of the main system.
 
-```
+  - `@extend var_to_unpack1, var_to_unpack2 = partial_sys = PartialSystem(var1 = value1)`
+    
+      + In this method: explicitly list the variables that should be unpacked, provide a name for the partial system and declare the base system.
+      + Note that only the arguments listed out in the declaration of the base system (here: `var1`) are added as the keyword arguments of the higher system.
 
 #### `@components` begin block
 
@@ -325,11 +329,11 @@ For example, the structure of `ModelC` is:
 julia> ModelC.structure
 Dict{Symbol, Any} with 10 entries:
   :components            => Any[Union{Expr, Symbol}[:model_a, :ModelA], Union{Expr, Symbol}[:model_array_a, :ModelA, :(1:N)], Union{Expr, Symbol}[:model_array_b, :ModelA, :(1:N)]]
-  :variables             => Dict{Symbol, Dict{Symbol, Any}}(:v=>Dict(:default=>:v_var, :type=>Real), :v_for_defaults=>Dict(:type=>Real))
+  :variables             => Dict{Symbol, Dict{Symbol, Any}}(:v=>Dict(:default=>:v_var, :type=>Real), :v_array=>Dict(:value=>nothing, :type=>Real, :size=>(:N, :M)), :v_for_defaults=>Dict(:type=>Real))
   :icon                  => URI("https://github.com/SciML/SciMLDocs/blob/main/docs/src/assets/logo.png")
-  :kwargs                => Dict{Symbol, Dict}(:f => Dict(:value => :sin), :N => Dict(:value => 2), :M => Dict(:value => 3), :v => Dict{Symbol, Any}(:value => :v_var, :type => Real), :v_for_defaults => Dict{Symbol, Union{Nothing, DataType}}(:value => nothing, :type => Real), :p1 => Dict(:value => nothing)),
-  :structural_parameters => Dict{Symbol, Dict}(:f => Dict(:value => :sin), :N => Dict(:value => 2), :M => Dict(:value => 3))
-  :independent_variable  => t
+  :kwargs                => Dict{Symbol, Dict}(:f=>Dict(:value=>:sin), :p2=>Dict(:value=>NoValue()), :N=>Dict(:value=>2), :M=>Dict(:value=>3), :v=>Dict{Symbol, Any}(:value=>:v_var, :type=>Real), :v_array=>Dict{Symbol, Any}(:value=>nothing, :type=>Real, :size=>(:N, :M)), :v_for_defaults=>Dict{Symbol, Union{Nothing, DataType}}(:value=>nothing, :type=>Real), :p1=>Dict(:value=>1))
+  :structural_parameters => Dict{Symbol, Dict}(:f=>Dict(:value=>:sin), :N=>Dict(:value=>2), :M=>Dict(:value=>3))
+  :independent_variable  => :t
   :constants             => Dict{Symbol, Dict}(:c=>Dict{Symbol, Any}(:value=>1, :type=>Int64, :description=>"Example constant."))
   :extend                => Any[[:p2, :p1], Symbol("#mtkmodel__anonymous__ModelB"), :ModelB]
   :defaults              => Dict{Symbol, Any}(:v_for_defaults=>2.0)

docs/src/tutorials/callable_params.md

@@ -0,0 +1,91 @@
+# Callable parameters and interpolating data
+
+ModelingToolkit.jl allows creating parameters that represent functions to be called. This
+is especially useful for including interpolants and/or lookup tables inside ODEs. In this
+tutorial we will create an `ODESystem` which employs callable parameters to interpolate data
+inside an ODE and go over the various syntax options and their implications.
+
+## Callable parameter syntax
+
+The syntax for callable parameters declared via `@parameters` must be one of the following
+
+ 1. `(fn::FType)(..)`
+ 2. `fn(::argType1, ::argType2, ...)`
+
+In the first case, the parameter is callable with any number/combination of arguments, and
+has a type of `FType` (the callable must be a subtype of `FType`). In the second case,
+the parameter is callable with as many arguments as declared, and all must match the
+declared types.
+
+By default, the return type of the callable symbolic is inferred to be `Real`. To change
+this, a `::retType` annotation can be added at the end.
+
+To declare a function that returns an array of values, the same array syntax can be used
+as for normal variables:
+
+```julia
+@parameters (foo::FType)(..)[1:3]::retType
+@parameters foo(::argType1, ::argType2)[1:3]::retType
+```
+
+`retType` here is the `eltype` of the returned array.
+
+## Storage of callable parameters
+
+Callable parameters declared with the `::FType` syntax will be stored in a `Vector{FType}`.
+Thus, if `FType` is non-concrete, the buffer will also be non-concrete. This is sometimes
+necessary to allow the value of the callable to be switched out for a different type without
+rebuilding the model. Typically this syntax is preferable when `FType` is concrete or
+a small union.
+
+Callable parameters declared with the `::argType1, ...` syntax will be stored in a
+`Vector{FunctionWrappers.FunctionWrapper{retType, Tuple{argType1, ...}}}`. This suffers
+the small overhead of a `FunctionWrapper` and restricts the signature of the callable,
+symbolic, but allows storing the parameter in a type-stable manner and swapping it out.
+This is preferable when the values that the callable can take do not share a common
+subtype. For example, when a callable can represent the activation of a neural network
+and can be `tanh`, `sigmoid`, etc. which have a common ancestor of `Function`.
+
+If both `::FType` and `::argType`s are specified, `::FType` takes priority. For example,
+
+```julia
+@parameters (p::LinearInterpolation)(::Real)
+```
+
+`p` will be stored in a `Vector{LinearInterpolation}`. If `::LinearInterpolation` was not
+specified, it would be stored in a `Vector{FunctionWrapper{Real, Tuple{Real}}}`.
+
+## Example using interpolations
+
+```@example callable
+using ModelingToolkit
+using OrdinaryDiffEq
+using DataInterpolations
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+ts = collect(0.0:0.1:10.0)
+spline = LinearInterpolation(ts .^ 2, ts)
+Tspline = typeof(spline)
+@variables x(t)
+@parameters (interp::Tspline)(..)
+
+@mtkbuild sys = ODESystem(D(x) ~ interp(t), t)
+```
+
+The derivative of `x` is obtained via an interpolation from DataInterpolations.jl. Note
+the parameter syntax. The `(..)` marks the parameter as callable. `(interp::Tspline)`
+indicates that the parameter is of type `Tspline`.
+
+```@example callable
+prob = ODEProblem(sys, [x => 0.0], (0.0, 1.0), [interp => spline])
+solve(prob)
+```
+
+Note that the the following will not work:
+
+```julia
+ODEProblem(
+    sys; [x => 0.0], (0.0, 1.0), [interp => LinearInterpolation(0.0:0.1:1.0, 0.0:0.1:1.0)])
+```
+
+Since the type of the spline doesn't match.

ext/MTKBifurcationKitExt.jl

@@ -113,7 +113,7 @@ function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem,
         # If the plot var is a normal state.
         if any(isequal(plot_var, var) for var in unknowns(nsys))
             plot_idx = findfirst(isequal(plot_var), unknowns(nsys))
-            record_from_solution = (x, p) -> x[plot_idx]
+            record_from_solution = (x, p; k...) -> x[plot_idx]
 
             # If the plot var is an observed state.
         elseif any(isequal(plot_var, eq.lhs) for eq in observed(nsys))
@@ -132,7 +132,7 @@ function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem,
     return BifurcationKit.BifurcationProblem(F,
         u0_bif_vals,
         p_vals,
-        (@lens _[bif_idx]),
+        (BifurcationKit.@optic _[bif_idx]),
         args...;
         record_from_solution = record_from_solution,
         J = J,