# DiffEqParamEstim with Modellingtoolkit

**URL:** https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788
**Category:** Modelling & Simulations
**Created:** [August 15, 2022, 5:48pm UTC](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788 "2022-08-15T17:48:28Z")
**Posts on this page:** 12
**Page:** 1

<div class="post-metadata">

### Author: ![xztt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xztt/32/33843_2.png) [@xztt](https://discourse.julialang.org/u/xztt)
#### Post date: [August 15, 2022, 5:48pm UTC](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788/1 "2022-08-15T17:48:28Z")

</div>

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.

```julia
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

---

<div class="post-metadata">

### Author: ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)
#### Post date: [August 15, 2022, 6:34pm UTC](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788/2 "2022-08-15T18:34:42Z")

</div>

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

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 15, 2022, 9:35pm UTC](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788/3 "2022-08-15T21:35:19Z")

</div>

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

---

<div class="post-metadata">

### Author: ![xztt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xztt/32/33843_2.png) [@xztt](https://discourse.julialang.org/u/xztt)
#### Post date: [August 16, 2022, 10:51am UTC](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788/4 "2022-08-16T10:51:05Z")

</div>

Thanks for this but when I run the following

```julia
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

```julia
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. 🙂

---

<div class="post-metadata">

### Author: ![xztt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xztt/32/33843_2.png) [@xztt](https://discourse.julialang.org/u/xztt)
#### Post date: [August 16, 2022, 11:31am UTC](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788/5 "2022-08-16T11:31:40Z")

</div>

> [@YingboMa](#):
>
> `prob.f.observed([LV.p, LV.V], u, p, t)`

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

if I run

```julia
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

```julia
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)

```julia
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.

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 16, 2022, 7:21pm UTC](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788/6 "2022-08-16T19:21:09Z")

</div>

> [@xztt](#):
>
> 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. 🙂

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

---

<div class="post-metadata">

### Author: ![xztt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xztt/32/33843_2.png) [@xztt](https://discourse.julialang.org/u/xztt)
#### Post date: [August 18, 2022, 8:58am UTC](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788/7 "2022-08-18T08:58:09Z")

</div>

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

```julia
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 

```

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 18, 2022, 11:35am UTC](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788/8 "2022-08-18T11:35:42Z")

</div>

Can you open an issue?

---

<div class="post-metadata">

### Author: ![xztt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xztt/32/33843_2.png) [@xztt](https://discourse.julialang.org/u/xztt)
#### Post date: [August 18, 2022, 11:47am UTC](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788/9 "2022-08-18T11:47:22Z")

</div>

Of course, in MTK or DiffEqParamEstim?

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 18, 2022, 3:31pm UTC](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788/10 "2022-08-18T15:31:06Z")

</div>

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

---

<div class="post-metadata">

### Author: ![xztt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xztt/32/33843_2.png) [@xztt](https://discourse.julialang.org/u/xztt)
#### Post date: [August 18, 2022, 4:50pm UTC](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788/11 "2022-08-18T16:50:08Z")

</div>

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

```julia
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.

> <https://github.com/SciML/ModelingToolkit.jl/issues/1761>
>
> 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

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 19, 2022, 2:22am UTC](https://discourse.julialang.org/t/diffeqparamestim-with-modellingtoolkit/85788/12 "2022-08-19T02:22:52Z")

</div>

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