From efea3d7ff0717544709ec7850f7bf043ab089c91 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Sun, 13 Apr 2025 13:36:53 +0200 Subject: [PATCH] Get observed variables with observables(sys) --- docs/src/tutorials/acausal_components.md | 7 +++++ src/ModelingToolkit.jl | 4 +-- src/systems/abstractsystem.jl | 28 +++++++++++++++++-- src/systems/diffeqs/abstractodesystem.jl | 2 +- src/systems/diffeqs/odesystem.jl | 4 +-- src/systems/diffeqs/sdesystem.jl | 2 +- src/systems/jumps/jumpsystem.jl | 2 +- src/systems/nonlinear/nonlinearsystem.jl | 2 +- .../optimization/constraints_system.jl | 2 +- .../optimization/optimizationsystem.jl | 2 +- src/systems/problem_utils.jl | 2 +- src/systems/systemstructure.jl | 2 +- test/serialization.jl | 2 +- test/structural_transformation/utils.jl | 4 +-- 14 files changed, 47 insertions(+), 18 deletions(-) diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md index b97500a3e9..751b678dae 100644 --- a/docs/src/tutorials/acausal_components.md +++ b/docs/src/tutorials/acausal_components.md @@ -320,6 +320,7 @@ sol = solve(prob) plot(sol) ``` +By default, this plots only the unknown variables that had to be solved for. However, what if we wanted to plot the timeseries of a different variable? Do not worry, that information was not thrown away! Instead, transformations like `structural_simplify` simply change unknown variables into observables which are @@ -346,3 +347,9 @@ or we can plot the timeseries of the resistor's voltage: ```@example acausal plot(sol, idxs = [rc_model.resistor.v]) ``` + +Although it may be more confusing than helpful here, we can of course also plot all unknown and observed variables together: + +```@example acausal +plot(sol, idxs = [unknowns(rc_model); observables(rc_model)]) +``` diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 9f69458528..d0f427bd14 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -85,8 +85,8 @@ import Symbolics: rename, get_variables!, _solve, hessian_sparsity, substituter, scalarize, getparent, hasderiv, hasdiff import DiffEqBase: @add_kwonly -export independent_variables, unknowns, parameters, full_parameters, continuous_events, - discrete_events +export independent_variables, unknowns, observables, parameters, full_parameters, + continuous_events, discrete_events @reexport using Symbolics @reexport using UnPack RuntimeGeneratedFunctions.init(@__MODULE__) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 29338d0722..2b58349ec4 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -430,7 +430,7 @@ function has_observed_with_lhs(sys, sym) if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing return haskey(ic.observed_syms_to_timeseries, sym) else - return any(isequal(sym), [eq.lhs for eq in observed(sys)]) + return any(isequal(sym), observables(sys)) end end @@ -489,7 +489,7 @@ function _all_ts_idxs!(ts_idxs, ::ScalarSymbolic, sys, sym::Symbol) if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing return _all_ts_idxs!(ts_idxs, sys, ic.symbol_to_variable[sym]) elseif is_variable(sys, sym) || is_independent_variable(sys, sym) || - any(isequal(sym), [getname(eq.lhs) for eq in observed(sys)]) + any(isequal(sym), getname.(observables(sys))) push!(ts_idxs, ContinuousTimeseries()) elseif is_timeseries_parameter(sys, sym) push!(ts_idxs, timeseries_parameter_index(sys, s).timeseries_idx) @@ -579,7 +579,7 @@ SymbolicIndexingInterface.constant_structure(::AbstractSystem) = true function SymbolicIndexingInterface.all_variable_symbols(sys::AbstractSystem) syms = variable_symbols(sys) - obs = getproperty.(observed(sys), :lhs) + obs = observables(sys) return isempty(obs) ? syms : vcat(syms, obs) end @@ -1411,6 +1411,7 @@ _nonum(@nospecialize x) = x isa Num ? x.val : x $(TYPEDSIGNATURES) Get the unknown variables of the system `sys` and its subsystems. +These must be explicitly solved for, unlike `observables(sys)`. See also [`ModelingToolkit.get_unknowns`](@ref). """ @@ -1677,6 +1678,14 @@ function controls(sys::AbstractSystem) isempty(systems) ? ctrls : [ctrls; reduce(vcat, namespace_controls.(systems))] end +""" +$(TYPEDSIGNATURES) + +Get the observed equations of the system `sys` and its subsystems. +These can be expressed in terms of `unknowns(sys)`, and do not have to be explicitly solved for. + +See also [`observables`](@ref) and [`ModelingToolkit.get_observed()`](@ref). +""" function observed(sys::AbstractSystem) obs = get_observed(sys) systems = get_systems(sys) @@ -1686,6 +1695,19 @@ function observed(sys::AbstractSystem) init = Equation[])] end +""" +$(TYPEDSIGNATURES) + +Get the observed variables of the system `sys` and its subsystems. +These can be expressed in terms of `unknowns(sys)`, and do not have to be explicitly solved for. +It is equivalent to all left hand sides of `observed(sys)`. + +See also [`observed`](@ref). +""" +function observables(sys::AbstractSystem) + return map(eq -> eq.lhs, observed(sys)) +end + Base.@deprecate default_u0(x) defaults(x) false Base.@deprecate default_p(x) defaults(x) false diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 62ddd12a08..378105d962 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1496,7 +1496,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractSystem, @warn errmsg end - uninit = setdiff(unknowns(sys), [unknowns(isys); getfield.(observed(isys), :lhs)]) + uninit = setdiff(unknowns(sys), [unknowns(isys); observables(isys)]) # TODO: throw on uninitialized arrays filter!(x -> !(x isa Symbolics.Arr), uninit) diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index cbce569a9b..01b0ca5fbb 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -47,7 +47,7 @@ struct ODESystem <: AbstractODESystem var_to_name::Any """Control parameters (some subset of `ps`).""" ctrls::Vector - """Observed variables.""" + """Observed equations.""" observed::Vector{Equation} """System of constraints that must be satisfied by the solution to the system.""" constraintsystem::Union{Nothing, ConstraintsSystem} @@ -532,7 +532,7 @@ function build_explicit_observed_function(sys, ts; vs = ModelingToolkit.vars(ts; op) namespace_subs = Dict() - ns_map = Dict{Any, Any}(renamespace(sys, eq.lhs) => eq.lhs for eq in observed(sys)) + ns_map = Dict{Any, Any}(renamespace(sys, obs) => obs for obs in observables(sys)) for sym in unknowns(sys) ns_map[renamespace(sys, sym)] = sym if iscall(sym) && operation(sym) === getindex diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl index 3fa1302630..c5299c28be 100644 --- a/src/systems/diffeqs/sdesystem.jl +++ b/src/systems/diffeqs/sdesystem.jl @@ -48,7 +48,7 @@ struct SDESystem <: AbstractODESystem var_to_name::Any """Control parameters (some subset of `ps`).""" ctrls::Vector - """Observed variables.""" + """Observed equations.""" observed::Vector{Equation} """ Time-derivative matrix. Note: this field will not be defined until diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl index 57a3aee7df..06f5e1b623 100644 --- a/src/systems/jumps/jumpsystem.jl +++ b/src/systems/jumps/jumpsystem.jl @@ -59,7 +59,7 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem ps::Vector """Array variables.""" var_to_name::Any - """Observed variables.""" + """Observed equations.""" observed::Vector{Equation} """The name of the system.""" name::Symbol diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index 856822492b..1d44c5c42a 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -32,7 +32,7 @@ struct NonlinearSystem <: AbstractTimeIndependentSystem ps::Vector """Array variables.""" var_to_name::Any - """Observed variables.""" + """Observed equations.""" observed::Vector{Equation} """ Jacobian matrix. Note: this field will not be defined until diff --git a/src/systems/optimization/constraints_system.jl b/src/systems/optimization/constraints_system.jl index 0f69e6d0b9..ae8577662d 100644 --- a/src/systems/optimization/constraints_system.jl +++ b/src/systems/optimization/constraints_system.jl @@ -33,7 +33,7 @@ struct ConstraintsSystem <: AbstractTimeIndependentSystem ps::Vector """Array variables.""" var_to_name::Any - """Observed variables.""" + """Observed equations.""" observed::Vector{Equation} """ Jacobian matrix. Note: this field will not be defined until diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl index be4567aee5..bfe15b62d7 100644 --- a/src/systems/optimization/optimizationsystem.jl +++ b/src/systems/optimization/optimizationsystem.jl @@ -31,7 +31,7 @@ struct OptimizationSystem <: AbstractOptimizationSystem ps::Vector """Array variables.""" var_to_name::Any - """Observed variables.""" + """Observed equations.""" observed::Vector{Equation} """List of constraint equations of the system.""" constraints::Vector{Union{Equation, Inequality}} diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 0750585905..e2ad55da13 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -1016,7 +1016,7 @@ function get_u0_p(sys, u0map = Dict(u0map) end if u0map isa Dict - allobs = Set(getproperty.(observed(sys), :lhs)) + allobs = Set(observables(sys)) if any(in(allobs), keys(u0map)) u0s_in_obs = filter(in(allobs), keys(u0map)) @warn "Observed variables cannot be assigned initial values. Initial values for $u0s_in_obs will be ignored." diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index d27e5c93a1..73cdaba255 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -721,7 +721,7 @@ function _structural_simplify!(state::TearingState, io; simplify = false, sys = ModelingToolkit.tearing( sys, state; simplify, mm, check_consistency, kwargs...) end - fullunknowns = [map(eq -> eq.lhs, observed(sys)); unknowns(sys)] + fullunknowns = [observables(sys); unknowns(sys)] @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns) ModelingToolkit.invalidate_cache!(sys), input_idxs diff --git a/test/serialization.jl b/test/serialization.jl index e10de51299..feb3c7e4e7 100644 --- a/test/serialization.jl +++ b/test/serialization.jl @@ -31,7 +31,7 @@ sys = include_string(@__MODULE__, str) # check answer ss = structural_simplify(rc_model) -all_obs = [o.lhs for o in observed(ss)] +all_obs = observables(ss) prob = ODEProblem(ss, [capacitor.v => 0.0], (0, 0.1)) sol = solve(prob, ImplicitEuler()) diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl index 6dfc107cc9..8894f0bbc9 100644 --- a/test/structural_transformation/utils.jl +++ b/test/structural_transformation/utils.jl @@ -51,8 +51,8 @@ end [D(x) ~ z[1] + z[2] + foo(z)[1], y[1] ~ 2t, y[2] ~ 3t, z ~ foo(y)], t) @test length(equations(sys)) == 1 @test length(observed(sys)) == 7 - @test any(eq -> isequal(eq.lhs, y), observed(sys)) - @test any(eq -> isequal(eq.lhs, z), observed(sys)) + @test any(obs -> isequal(obs, y), observables(sys)) + @test any(obs -> isequal(obs, z), observables(sys)) prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [foo => _tmp_fn]) @test_nowarn prob.f(prob.u0, prob.p, 0.0)