DiffEqParamEstim with Modellingtoolkit

xztt · August 15, 2022, 5:48pm

Hi

I am generating a model using MTK and would like to generate a simple loss function for it.

Is there a way in which I can access the observed variables in my solution when solving my problem.

I know given a standard ODE problem I would use :save_idxs however given I am trying to generate a loss function of observed variables it is proving tricky.

From the code below I am hard coding the “nosiy” solution vector however when I define obj I would need to be able to access these observed variables LV.p and LV.V.

all works well when following the example and just using the output solution vector.

p = [1.5, 0.03, 1.32, 21.9, 0.25755, 0.4318, 1.672, 0.033, 0.006, 1.11, 1.13, 11.0]
u0 = [MCFP, -MCFP, -MCFP]
tspan = (0.0, 20.0)
prob = ODEProblem(circ_sys, u0, tspan, p)
t = collect(range(0,stop=20,length=500))
@time sol = solve(prob, Tsit5(), reltol=1e-12, abstol=1e-12)

V = [Vector{Float64}(undef,2) for _ in 1:500]
for i in 1:500
    V[i][1] = sol[LV.p][i]
end
for i in 1:500
    V[i][2] = sol[LV.V][i]
end

function gen_data(sol,t)
    randomized = VectorOfArray([(sol(t[i]) + 0.01rand(3)) for i in 1:length(t)])
    data = convert(Array,randomized)
end 
aggregate_data = convert(Array,VectorOfArray([gen_data(sol,t) for i in 1:100]))

distributions  = [fit_mle(Normal,aggregate_data[i,j,:]) for i in 1:3, j in 1:500]

obj = build_loss_objective(prob, Tsit5(), reltol=1e-12, abstol=1e-12, LogLikeLoss(t,distributions))

In conclusion i think my question boils down to can I access the observed variables as part of the solution vector when using the function build_loss_objective?

Cheers
Harry

YingboMa · August 15, 2022, 6:34pm

You can use prob.f.observed([LV.p, LV.V], u, p, t) to get their values.

ChrisRackauckas · August 15, 2022, 9:35pm

And for the fitting, fit against observed values using save_idxs = [observedvar1, observedvar2] etc. Just use the symbolics as indices.

xztt · August 16, 2022, 10:51am

Thanks for this but when I run the following

cost_func = build_loss_objective(prob, Tsit5(), reltol=1e-12, abstol=1e-12, L2Loss(t,data), save_idxs = [LV.p,LV.V])

prange = 0.1:0.01:5.5
plot(prange, [cost_func([i, 0.03, 1.32, 21.9, 0.25755, 0.4318, 1.672, 0.033, 0.006, 1.11, 1.13, 11.0]) for i in prange], lw = 3, title = "Emax Liklihood (Rest at base state)", xlabel = "Emax", ylabel = "Objective Function Value - LS")

It throws back the error

RROR: TypeError: non-boolean (Num) used in boolean context
Stacktrace:
  [1] checkbounds
    @ ./abstractarray.jl:668 [inlined]
  [2] _getindex
    @ ./multidimensional.jl:874 [inlined]
  [3] getindex
    @ ./abstractarray.jl:1241 [inlined]
  [4] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#462"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x71e6a2a4, 0x6d290c6c, 0x76cb83a3, 0x2a859523, 0xd1e2cf02)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3069bfc2, 0x52123b63, 0x6b56c11a, 0xc73a3d04, 0x22966743)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#472#generated_observed#469"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Vector{Float64}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Vector{Num}, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Float64, reltol::Float64, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/5B7d7/src/solve.jl:226
  [5] #__solve#562
    @ ~/.julia/packages/OrdinaryDiffEq/5B7d7/src/solve.jl:4 [inlined]
  [6] #solve_call#28
    @ ~/.julia/packages/DiffEqBase/72SnT/src/solve.jl:449 [inlined]
  [7] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#462"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x71e6a2a4, 0x6d290c6c, 0x76cb83a3, 0x2a859523, 0xd1e2cf02)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3069bfc2, 0x52123b63, 0x6b56c11a, 0xc73a3d04, 0x22966743)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#472#generated_observed#469"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::Vector{Float64}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Any, NTuple{6, Symbol}, NamedTuple{(:saveat, :save_everystep, :dense, :reltol, :abstol, :save_idxs), Tuple{Vector{Float64}, Bool, Bool, Float64, Float64, Vector{Num}}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/72SnT/src/solve.jl:802
  [8] #solve#33
    @ ~/.julia/packages/DiffEqBase/72SnT/src/solve.jl:772 [inlined]
  [9] (::DiffEqParamEstim.var"#37#42"{Nothing, Bool, Int64, typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR), Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:reltol, :abstol, :save_idxs), Tuple{Float64, Float64, Vector{Num}}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#462"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x71e6a2a4, 0x6d290c6c, 0x76cb83a3, 0x2a859523, 0xd1e2cf02)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3069bfc2, 0x52123b63, 0x6b56c11a, 0xc73a3d04, 0x22966743)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#472#generated_observed#469"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, L2Loss{Vector{Float64}, Matrix{Float64}, Nothing, Nothing, Nothing}, Nothing, Tuple{}})(p::Vector{Float64})
    @ DiffEqParamEstim ~/.julia/packages/DiffEqParamEstim/g7Ey8/src/build_loss_objective.jl:45
 [10] (::DiffEqObjective{DiffEqParamEstim.var"#37#42"{Nothing, Bool, Int64, typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR), Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:reltol, :abstol, :save_idxs), Tuple{Float64, Float64, Vector{Num}}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#462"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x71e6a2a4, 0x6d290c6c, 0x76cb83a3, 0x2a859523, 0xd1e2cf02)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3069bfc2, 0x52123b63, 0x6b56c11a, 0xc73a3d04, 0x22966743)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#472#generated_observed#469"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, L2Loss{Vector{Float64}, Matrix{Float64}, Nothing, Nothing, Nothing}, Nothing, Tuple{}}, DiffEqParamEstim.var"#41#47"{DiffEqParamEstim.var"#37#42"{Nothing, Bool, Int64, typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR), Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:reltol, :abstol, :save_idxs), Tuple{Float64, Float64, Vector{Num}}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#462"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x71e6a2a4, 0x6d290c6c, 0x76cb83a3, 0x2a859523, 0xd1e2cf02)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3069bfc2, 0x52123b63, 0x6b56c11a, 0xc73a3d04, 0x22966743)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#472#generated_observed#469"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, L2Loss{Vector{Float64}, Matrix{Float64}, Nothing, Nothing, Nothing}, Nothing, Tuple{}}}})(x::Vector{Float64})
    @ DiffEqParamEstim ~/.julia/packages/DiffEqParamEstim/g7Ey8/src/build_loss_objective.jl:25
 [11] (::var"#210#211")(i::Float64)
    @ Main ./none:0
 [12] iterate
    @ ./generator.jl:47 [inlined]
 [13] collect(itr::Base.Generator{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, var"#210#211"})
    @ Base ./array.jl:787
 [14] top-level scope
    @ ~/Desktop/PhD/Year 1/Sem 2/Paper 1/Matlab/Nikolai reproduction/Nik_recon.jl:1045

I have to assume this is from using the save_idxs in the build_loss_objective. Sorry for been dense I feel like I must be missing something obvious here.

xztt · August 16, 2022, 11:31am

I am not actually sure that this works unless I am missing something

if I run

prob.f.observed(LV.p, u0,  p,t)

it returns the inital value of the observed as expected

running as you have written above returns the method error

MethodError: no method matching +(::Float64, ::Vector{Float64})
For element-wise addition, use broadcasting with dot syntax: scalar .+ array
Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...) at operators.jl:591
  +(::T, ::T) where T<:Union{Float16, Float32, Float64} at float.jl:383
  +(::Union{Float16, Float32, Float64}, ::BigFloat) at mpfr.jl:414
  ...

and running as

prob.f.observed((LV.p,LV.V), u0, p,t)

returns a tuple `(LV₊V(t), LV₊p(t))`

Which causes issues when trying to plot

plot(0:10, t->prob.f.observed((LV.V,LV.p), u0, p,t))
Maybe I am just misunderstanding this whole thing which is also highly likely.

ChrisRackauckas · August 16, 2022, 7:21pm

Can you make an MWE? That seems like something that could/should be working.

xztt · August 18, 2022, 8:58am

Of course below should demonstrate my problem just taken the RC-Circuit from the MTK tutorial

using ModelingToolkit, Plots, DifferentialEquations, DiffEqParamEstim

@variables t
@connector function Pin(;name)
    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
    ODESystem(Equation[], t, sts, []; name=name)
end

function Ground(;name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name=name), g)
end

function OnePort(;name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end

function Resistor(;name, R = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R=R
    eqs = [
           v ~ i * R
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Capacitor(;name, C = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters C=C
    D = Differential(t)
    eqs = [
           D(v) ~ i / C
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function ConstantVoltage(;name, V = 1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V=V
    eqs = [
           V ~ v
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

R = 1.0
C = 1.0
V = 1.0
@named resistor = Resistor(R=R)
@named capacitor = Capacitor(C=C)
@named source = ConstantVoltage(V=V)
@named ground = Ground()

rc_eqs = [
          connect(source.p, resistor.p)
          connect(resistor.n, capacitor.p)
          connect(capacitor.n, source.n)
          connect(capacitor.n, ground.g)
         ]

@named _rc_model = ODESystem(rc_eqs, t)
@named rc_model = compose(_rc_model,
                          [resistor, capacitor, source, ground])
sys = structural_simplify(rc_model)
u0 = [
      capacitor.v => 0.0
     ]
prob = ODEProblem(sys, u0, (0, 10.0))
#sol = solve(prob, Tsit5(), save_idxs = [capacitor.v]) can not index a solution using symbolics, can also not index at observed values 
t = collect(range(0,stop=10,length=200))
sol = solve(prob, Tsit5(), saveat = t)

plot(sol, idxs = [resistor.p.i,capacitor.v ]) #here resisitor is observed variable 

V = [Vector{Float64}(undef,2) for _ in 1:200]
for i in 1:200
    V[i][1] = sol[resistor.p.i][i]
end
for i in 1:200
    V[i][2] = sol[capacitor.v ][i]
end
using RecursiveArrayTools # for VectorOfArray
randomized = VectorOfArray([(V[:][i] + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
cost_func = build_loss_objective(prob,Tsit5(),L2Loss(t,data),save_idxs = [resistor.p.i,capacitor.v ])
prange = 0.8:0.01:1.2
plot(prange, [cost_func([i, 1.0, 1.0]) for i in prange], lw = 3, title = "R Liklihood (Rest at base state)", xlabel = "Emax", ylabel = "Objective Function Value - LS")
## Working cost function ##
randomized = VectorOfArray([(sol(t[i]) + .01randn(1)) for i in 1:length(t)])
data = convert(Array,randomized)
cost_func = build_loss_objective(prob,Tsit5(),L2Loss(t,data))
prange = 0.8:0.01:1.2
plot(prange, [cost_func([i, 1.0, 1.0]) for i in prange], lw = 3, title = "R Liklihood (Rest at base state)", xlabel = "Emax", ylabel = "Objective Function Value - LS")

## Appears the symbolic indexing in the cost function causes problems

ChrisRackauckas · August 18, 2022, 11:35am

Can you open an issue?

xztt · August 18, 2022, 11:47am

Of course, in MTK or DiffEqParamEstim?

ChrisRackauckas · August 18, 2022, 3:31pm

Maybe MTK. Can this be reduced to just an issue of solve with save_idxs on symbolic variables? That’s my guess.

xztt · August 18, 2022, 4:50pm

Yes I think so because even running with the state which is returned by structural_simplify still returns an error i.e

sol = solve(prob, Tsit5(), saveat = t, save_idxs = [capacitor.v])
ERROR: TypeError: non-boolean (Num) used in boolean context
Stacktrace:
 [1] checkbounds
   @ ./abstractarray.jl:668 [inlined]
 [2] _getindex
   @ ./multidimensional.jl:874 [inlined]
 [3] getindex
   @ ./abstractarray.jl:1241 [inlined]
 [4] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#462"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x786907a7, 0x70d17a6b, 0x2e0e28cf, 0xa034a993, 0xa01b22e4)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xfcaf3b52, 0x74cb17c4, 0xcd5671ed, 0x58472ffd, 0xca228a9b)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#472#generated_observed#469"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Vector{Float64}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Vector{Num}, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/5B7d7/src/solve.jl:226
 [5] #__solve#562
   @ ~/.julia/packages/OrdinaryDiffEq/5B7d7/src/solve.jl:4 [inlined]
 [6] #solve_call#28
   @ ~/.julia/packages/DiffEqBase/72SnT/src/solve.jl:449 [inlined]
 [7] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#462"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x786907a7, 0x70d17a6b, 0x2e0e28cf, 0xa034a993, 0xa01b22e4)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xfcaf3b52, 0x74cb17c4, 0xcd5671ed, 0x58472ffd, 0xca228a9b)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#472#generated_observed#469"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::Vector{Float64}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Vector, Tuple{Symbol, Symbol}, NamedTuple{(:saveat, :save_idxs), Tuple{Vector{Float64}, Vector{Num}}}})
   @ DiffEqBase ~/.julia/packages/DiffEqBase/72SnT/src/solve.jl:802
 [8] #solve#33
   @ ~/.julia/packages/DiffEqBase/72SnT/src/solve.jl:772 [inlined]
 [9] top-level scope
   @ Untitled-1:83

Issue can be found below.

github.com/SciML/ModelingToolkit.jl

Issue with solve with save_idxs on symbolic variables

opened 04:49PM - 18 Aug 22 UTC

H-Sax

Following the MTK tutorial ``` using ModelingToolkit, Plots, DifferentialEq…uations, DiffEqParamEstim @variables t @connector function Pin(;name) sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow] ODESystem(Equation[], t, sts, []; name=name) end function Ground(;name) @named g = Pin() eqs = [g.v ~ 0] compose(ODESystem(eqs, t, [], []; name=name), g) end function OnePort(;name) @named p = Pin() @named n = Pin() sts = @variables v(t)=1.0 i(t)=1.0 eqs = [ v ~ p.v - n.v 0 ~ p.i + n.i i ~ p.i ] compose(ODESystem(eqs, t, sts, []; name=name), p, n) end function Resistor(;name, R = 1.0) @named oneport = OnePort() @unpack v, i = oneport ps = @parameters R=R eqs = [ v ~ i * R ] extend(ODESystem(eqs, t, [], ps; name=name), oneport) end function Capacitor(;name, C = 1.0) @named oneport = OnePort() @unpack v, i = oneport ps = @parameters C=C D = Differential(t) eqs = [ D(v) ~ i / C ] extend(ODESystem(eqs, t, [], ps; name=name), oneport) end function ConstantVoltage(;name, V = 1.0) @named oneport = OnePort() @unpack v = oneport ps = @parameters V=V eqs = [ V ~ v ] extend(ODESystem(eqs, t, [], ps; name=name), oneport) end R = 1.0 C = 1.0 V = 1.0 @named resistor = Resistor(R=R) @named capacitor = Capacitor(C=C) @named source = ConstantVoltage(V=V) @named ground = Ground() rc_eqs = [ connect(source.p, resistor.p) connect(resistor.n, capacitor.p) connect(capacitor.n, source.n) connect(capacitor.n, ground.g) ] @named _rc_model = ODESystem(rc_eqs, t) @named rc_model = compose(_rc_model, [resistor, capacitor, source, ground]) sys = structural_simplify(rc_model) u0 = [ capacitor.v => 0.0 ] prob = ODEProblem(sys, u0, (0, 10.0)) #sol = solve(prob, Tsit5(), save_idxs = [capacitor.v]) can not index a solution using symbolics, can also not index at observed values t = collect(range(0,stop=10,length=200)) sol = solve(prob, Tsit5(), saveat = t, save_idxs = [capacitor.v]) sol = solve(prob, Tsit5(), saveat = t, save_idxs = [1]) #This one works sol = solve(prob, Tsit5(), save_idxs = [resistor.p.i]) # Observed Variable plot(sol, idxs = [resistor.p.i,capacitor.v ]) ``` Using save_idxs only works using numerical indexing values i.e `save_idxs =[1]`. I also want to point out that we can not save at the observed symbolic to return the observed variable solution in the solution vector. This issue has come out of a discussion relating to DiffEqParamEstim which can be found [here ](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788/7) Cheers, Harry

ChrisRackauckas · August 19, 2022, 2:22am

Thanks for honing in one this. That’s helpful.