Named state space sets

Author: Datseris

CHANGELOG.md

@@ -1,3 +1,8 @@
+# v3.15
+
+- New function `named_variables(ds)` for getting the variable names of an MTK-generated dynamical system.
+- `trajectory` function now automatically makes a named state space set if the dynamical system is MTK-generated.
+
 # v3.14.0
 
 - `set_parameter!` and `current_parameter` will now throw an informative error message explicitly naming the parameter if the user tries to get/set a symbolic parameter that does not exist in the MTK-generated dynamical system.

Project.toml

@@ -1,7 +1,7 @@
 name = "DynamicalSystemsBase"
 uuid = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 repo = "https://github.com/JuliaDynamics/DynamicalSystemsBase.jl.git"
-version = "3.14.0"
+version = "3.15.0"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -21,7 +21,7 @@ OrdinaryDiffEqTsit5 = "1.1"
 Reexport = "1"
 Roots = "1, 2"
 SciMLBase = "1.19.5, 2"
-StateSpaceSets = "2"
+StateSpaceSets = "2.5"
 Statistics = "1"
 StochasticDiffEq = "6.66.0"
 SymbolicIndexingInterface = "0.3.4"

docs/src/index.md

@@ -30,6 +30,7 @@ initial_time
 isinplace(::DynamicalSystem)
 successful_step
 referrenced_sciml_model
+named_variables
 ```
 
 ```@docs

src/core/dynamicalsystem_interface.jl

@@ -44,7 +44,7 @@ from an experimental measurement of a dynamical system with an unknown dynamic r
 
 See also the DynamicalSystems.jl tutorial online for examples making dynamical systems.
 
-## Integration with ModelingToolkit.jl
+## Integration with ModelingToolkit.jl (MTK)
 
 Dynamical systems that have been constructed from `DEProblem`s that themselves
 have been constructed from ModelingToolkit.jl keep a reference to the symbolic
@@ -93,6 +93,7 @@ unless when developing new algorithm implementations that use dynamical systems.
 - [`isinplace`](@ref)
 - [`successful_step`](@ref)
 - [`referrenced_sciml_model`](@ref)
+- [`named_variables`](@ref)
 
 ### API - alter status
 
@@ -126,7 +127,7 @@ errormsg(ds) = error("Not yet implemented for dynamical system of type $(nameof(
 
 export current_state, initial_state, current_parameters, current_parameter, initial_parameters, isinplace,
     current_time, initial_time, successful_step, isdeterministic, isdiscretetime, dynamic_rule,
-    reinit!, set_state!, set_parameter!, set_parameters!, step!, observe_state, referrenced_sciml_model
+    reinit!, set_state!, set_parameter!, set_parameters!, step!, observe_state, referrenced_sciml_model, named_variables
 
 ###########################################################################################
 # Symbolic support
@@ -148,10 +149,27 @@ referrenced_sciml_model(::Nothing) = nothing
 
 # return true if there is an actual referrenced system
 has_referrenced_model(prob::SciMLBase.DEProblem) = has_referrenced_model(referrenced_sciml_model(prob))
+has_referrenced_model(prob::DynamicalSystem) = has_referrenced_model(referrenced_sciml_model(prob))
 has_referrenced_model(::Nothing) = false
 has_referrenced_model(::SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}) = false
 has_referrenced_model(sys) = true
 
+"""
+    named_variables(ds::DynamicalSystem)
+
+If `ds` is constructed via MTK, return a vector of the variable names (as symbols).
+Otherwise return `nothing`. If `X` is a `StateSpaceSet`, you can always do
+```julia
+X = StateSpaceSet(X; names = named_variables(ds))
+```
+in downstream functions to name a set coming from `ds` (if possible).
+"""
+function named_variables(ds::DynamicalSystem)
+    mtk = referrenced_sciml_model(ds)
+    isnothing(mtk) && return nothing
+    return SymbolicIndexingInterface.getname.(SymbolicIndexingInterface.variable_symbols(mtk))
+end
+
 ###########################################################################################
 # API - obtaining information from the system
 ###########################################################################################

src/core/trajectory.jl

@@ -5,18 +5,15 @@ export trajectory
 
 Evolve `ds` for a total time of `T` and return its trajectory `X`, sampled at equal time
 intervals, and corresponding time vector. `X` is a `StateSpaceSet`.
-Optionally provide a starting state `u0` which is `current_state(ds)` by default.
-
 The returned time vector is `t = (t0+Ttr):Δt:(t0+Ttr+T)`.
 
-If time evolution diverged, or in general failed, before `T`,
+Optionally provide a starting state `u0` which is `current_state(ds)` by default.
+
+If time evolution diverged or in general failed before `T`,
 the remaining of the trajectory is set to the last valid point.
 
-`trajectory` is a very simple function provided for convenience.
-For continuous time systems, it doesn't play well with callbacks,
-use `DifferentialEquations.solve` if you want a trajectory/timeseries
-that works with callbacks, or in general you want more flexibility in the generated
-trajectory (but remember to convert the output of `solve` to a `StateSpaceSet`).
+The dimensions of `X` are automatically named if `ds` referrences an MTK model,
+see [`DynamicalSystem`](@ref).
 
 ## Keyword arguments
 
@@ -32,17 +29,28 @@ trajectory (but remember to convert the output of `solve` to a `StateSpaceSet`).
   Defaults to `1:dimension(ds)` (all dynamic variables).
   Note: if you mix integer and symbolic indexing be sure to initialize the array
   as `Any` so that integers `1, 2, ...` are not converted to symbolic expressions.
+
+## Description
+
+`trajectory` is a very simple function provided for convenience.
+For continuous time systems, it doesn't play well with callbacks,
+use `DifferentialEquations.solve` if you want a trajectory
+that works with callbacks, or in general you want more flexibility in the generated
+trajectory (but remember to convert the output of `solve` to a `StateSpaceSet`).
 """
 function trajectory(ds::DynamicalSystem, T, u0 = current_state(ds);
         save_idxs = nothing, t0 = initial_time(ds), kwargs...
     )
     accessor = svector_access(save_idxs)
     reinit!(ds, u0; t0)
     if isdiscretetime(ds)
-        trajectory_discrete(ds, T; accessor, kwargs...)
+        X, tvec = trajectory_discrete(ds, T; accessor, kwargs...)
     else
-        trajectory_continuous(ds, T; accessor, kwargs...)
+        X, tvec = trajectory_continuous(ds, T; accessor, kwargs...)
     end
+    # name automatically if possible
+    X = StateSpaceSet(X; names = named_variables(ds))
+    return X, tvec
 end
 
 function trajectory_discrete(ds, T;