Support array variables or alternatives to deploy with MethodOfLines.jl

[kkakosim]

What is the proper way to utilize PEtab with other packages that use array variables?

at a Discourse Discussion @ChrisRackauckas highlighted that PEtab may not support ODESystem generated by ModellingToolkit which contains array variables.

The discussed MWE, also provided herein, breaks with the below error when executing p0 = generate_startguesses(petab_problem, 1).

ERROR: CanonicalIndexError: setindex! not defined for Symbolics.Arr{Num, 1}
Stacktrace:
  [1] error_if_canonical_setindex(::IndexCartesian, A::Symbolics.Arr{Num, 1}, ::Int64)
    @ Base .\abstractarray.jl:1403
  [2] setindex!(A::Symbolics.Arr{Num, 1}, v::Float64, I::Int64)
    @ Base .\abstractarray.jl:1392
  [3] macro expansion
    @ .\none:7 [inlined]
  [4] macro expansion
    @ C:\Users\kkakosim\.julia\packages\RuntimeGeneratedFunctions\Yo8zx\src\RuntimeGeneratedFunctions.jl:163 [inlined]
  [5] macro expansion
    @ .\none:0 [inlined]
  [6] generated_callfunc
    @ .\none:0 [inlined]
  [7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Vector{…})
    @ RuntimeGeneratedFunctions C:\Users\kkakosim\.julia\packages\RuntimeGeneratedFunctions\Yo8zx\src\RuntimeGeneratedFunctions.jl:150
  [8] change_ode_parameters!(p_ode_problem::Vector{…}, u0::Vector{…}, θ::Vector{…}, θ_indices::PEtab.ParameterIndices, petab_model::PEtabModel{…})
    @ PEtab C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\Common.jl:252
  [9] compute_cost_solve_ODE(θ_dynamic::SubArray{…}, θ_sd::SubArray{…}, θ_observable::SubArray{…}, θ_non_dynamic::SubArray{…}, ode_problem::ODEProblem{…}, ode_solver::ODESolver, ss_solver::SteadyStateSolver{…}, petab_model::PEtabModel{…}, simulation_info::PEtab.SimulationInfo{…}, θ_indices::PEtab.ParameterIndices, measurement_info::PEtab.MeasurementsInfo{…}, parameter_info::PEtab.ParametersInfo, petab_ODE_cache::PEtab.PEtabODEProblemCache{…}, petab_ODESolver_cache::PEtab.PEtabODESolverCache; compute_cost::Bool, compute_hessian::Bool, compute_gradient_θ_dynamic::Bool, compute_residuals::Bool, exp_id_solve::Vector{…})
    @ PEtab C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\Objective\Objective.jl:78
 [10] compute_cost_solve_ODE
    @ C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\Objective\Objective.jl:36 [inlined]
 [11] compute_cost(θ_est::Vector{…}, ode_problem::ODEProblem{…}, ode_solver::ODESolver, ss_solver::SteadyStateSolver{…}, petab_model::PEtabModel{…}, simulation_info::PEtab.SimulationInfo{…}, θ_indices::PEtab.ParameterIndices, measurement_info::PEtab.MeasurementsInfo{…}, parameter_info::PEtab.ParametersInfo, prior_info::PEtab.PriorInfo, petab_ODE_cache::PEtab.PEtabODEProblemCache{…}, petab_ODESolver_cache::PEtab.PEtabODESolverCache, exp_id_solve::Vector{…}, compute_cost::Bool, compute_hessian::Bool, compute_residuals::Bool)
    @ PEtab C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\Objective\Objective.jl:19
 [12] #451
    @ C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\PEtabODEProblem\Create.jl:358 [inlined]
 [13] generate_startguesses(petab_problem::PEtabODEProblem{…}, n_multistarts::Int64; sampling_method::QuasiMonteCarlo.LatinHypercubeSample, sample_from_prior::Bool, allow_inf_for_startguess::Bool, verbose::Bool)
    @ PEtab C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\Calibrate\Common.jl:288
 [14] generate_startguesses(petab_problem::PEtabODEProblem{…}, n_multistarts::Int64)
    @ PEtab C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\Calibrate\Common.jl:255
 [15] top-level scope
    @ p:\Chemical Engineering\Air_Quality_Eng_R_Team\!GITHUB\tamuq-chen-secarelab-receiver\aysha\MWE PEtab.jl:82
Some type information was truncated. Use `show(err)` to see complete types.

using ModelingToolkit, MethodOfLines, Plots, DomainSets, OrdinaryDiffEq, PEtab, DataFrames

pythonplot()

# Parameters, variables, and derivatives
@parameters k a To
@variables t x T(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

Tamb = 273 #K
To_test = 1000 #K
h = 30 #W/m^2K arbritrary value
L = 0.1 #m
ρ = 7500 #kg/m^3
cp = 420 #J/kgK
kin = 15 #W/m K
ain = kin/ρ/cp #m^2/s

param = [k => kin, a => ain, To => To_test]

# 1D PDE and boundary conditions
eq  = Dt(T(t, x)) ~ a * Dxx(T(t, x))
bcs = [T(0, x) ~ Tamb, #initial condition
        Dx(T(t, 0)) ~ -h/k *(T(t, 0) - Tamb), #cold end
        T(t, L) ~ To] #hot end

# Space and time domains
domains = [t ∈ Interval(0.0, 300.0),
           x ∈ Interval(0.0, L)]

# PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [T(t, x)], param)

# Method of lines discretization
dx = L/100
order = 2
discretization = MOLFiniteDifference([x => dx], t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

# Solve ODE problem
sol = solve(prob, Tsit5(), saveat=1.)

# Plot the solution
px = collect(sol[x])
pt = sol[t]
pT = sol[T(t,x)]
plt = plot(sol[t], sol[T(t,x)][:, 1], xlabel="Time", ylabel="Temperature", title="1D Heat Equation")
cnt = contourf(px, pt, pT, xlabel="Position", ylabel="Time", title="1D Heat Equation")
display(cnt)

#extract simulated experiments to create measurements for PEtab
idx =[3, 50, 97]
mx = px[idx]
mT = pT[:, idx]

#Create PEtab problem

odesys = symbolic_discretize(pdesys, discretization)
simpsys = structural_simplify(odesys[1])
@unpack T, k, a, To = simpsys
obs_T3 = PEtabObservable(T[3], 0.5)
observables = Dict("obs_T3" => obs_T3)
_k = PEtabParameter(:k, lb=0.1, ub=100., scale=:lin)
_a = PEtabParameter(:a, lb=1E-7, ub=1E-4, scale=:log10)
params = [_k , _a]  
E0 = Dict(:To => 1000.)
E1 = Dict(:To => 900.)
sim_cond = Dict("c0" => E0, "c1" => E1)
measurements = DataFrame(
    simulation_id = ["c0", "c1"],
    obs_id=["obs_T3", "obs_T3"],
    time = [300., 300.],
    measurement=[400., 420.])

petab_model = PEtabModel(simpsys, sim_cond, observables, measurements, params, verbose=false)
petab_problem = PEtabODEProblem(petab_model, verbose=false)
p0 = generate_startguesses(petab_problem, 1)
res = calibrate_model(petab_problem, p0, IpoptOptimiser(false),
                      options=IpoptOptions(max_iter = 1000))
println(res)

[sebapersson]

Thanks for reporting this.

The issue here is that to track species PEtab.jl uses states(simpsys) which here returns (T(t))[2]... which due to its format is hard to track (for ODESystem we are used to states such as u1, u2, u3 ..., I actually never thought about PDE:s when developing PEtab.jl, but it makes sense that we can support it :)

Good news is that this can be worked around by internally renaming (T(t))[2]... to something like T(t)__2.. (or similar) and then all the internal tracking should work just fine.

I am on a conference this week, but can hopefully look at this next week.

[ChrisRackauckas]

Doing that wouldn't be a great idea for large models as it would remove compiler optimizations in MTK v9 and JuliaSimCompiler.

[sebapersson]

The renaming would not actually be for the ODESystem, rather it would only be internally for the building the observational and initial value functions so the sys would be unchanged, which should not affect optimizations in MTK v9?

[ChrisRackauckas]

I don't get your statement, though the observed function generation does have some array optimizations.

[sebapersson]

What I am meaning is that PEtab.jl itself builds based on PEtabObservable an observation function, so for that to work internally I would have to rename (never exposed to the user). But I guess in general is that with PDE:s the observable is probably likely some sum of all states?

[ChrisRackauckas]

But I guess in general is that with PDE:s the observable is probably likely some sum of all states?

It's generally an array variable. But if you are calling using the SII functions then it should canonicalize to just one vector. I don't see a SymbolicIndexingInterface.jl dep here at all so I'm very skeptical that what you have setup is generally correct.

[TorkelE]

I think (but could very well be wrong) the reason this is not based on SII is that it still awaits MTKv9 to be fully ready. Since we are still updating Catalyst to work with MTKv9, and this partially works on Catalyst models, I guess Sebastian's plan is to wait with updating PEtab to MTKv9 until after Catalyst have done so.

[sebapersson]

Exactly, we do not used SII as we are still waiting on Catalyst for updating to MTKv9 (as Catalyst is core component of SBMLImporter.jl used for PEtab.jl). Taken together, along with MTKv9 update I can make PEtab.jl more compatible with MethodOfLines generated models.