Adding DataInterpolation functions into a model with MTK v9

[langestefan]

With the new update to MTK v9 we can use the much more efficient method of defining data interpolations as parameter values.

How to do so can be seen here:
SciML/ModelingToolkit.jl#2823 (comment)

I have tried to replicate this with my model. I could not find a way to add the parameter map to the PEtab model, so instead I added it to the ODESystem:

function Ti(input_funcs::Dict)
    @parameters begin
        C_i
        R_ia
        A_w
        T_i_0
    end
    @parameters begin
        (T_a_f)(..) = input_funcs[:T_a]
        (P_h_f)(..) = input_funcs[:P_h]
        (P_s_f)(..) = input_funcs[:P_s]
    end
    @variables begin
        T_i(t) = T_i_0
        T_a(t)
        P_h(t)
        P_s(t)
    end

    # define model ODE system
    eqs = [
        C_i * D(T_i) ~ (1 / R_ia) * (T_a - T_i) + A_w * P_s + P_h
    ]

    # define inputs 
    inputs = [
        T_a ~ T_a_f(t)
        P_h ~ P_h_f(t)
        P_s ~ P_s_f(t)
    ]
    append!(eqs, inputs)

    @named sys = ODESystem(eqs, t)
    return sys
end

With input_funcs given by:

    # Create input functions
    T_a = DataInterpolations.AkimaInterpolation(df.T_a, df.t_hour)
    P_h = DataInterpolations.AkimaInterpolation(df.P_h, df.t_hour)
    P_s = DataInterpolations.AkimaInterpolation(df.P_s, df.t_hour)

    # Create a dictionary of input functions
    input_funcs = Dict(:T_a => T_a, :P_h => P_h, :P_s => P_s)

This can generate a system model, but fails at this step:

    petab_model = PEtabModel(sys,
        observables,
        measurements,
        params;
        verbose=true
    )

┌Info: Building PEtabModel for model ODESystemModel
UndefVarError: PEtabInuptError not defined

Stacktrace:
[1] _get_parametermap(sys::ODESystem, parametermap_input::Nothing)
@ PEtab C:\Users\lange.julia\packages\PEtab\MalGN\src\julia_input\maps.jl:44
[2] _PEtabModel(sys::ODESystem, simulation_conditions::Dict{String, Dict{Any, Any}}, observables::Dict{String, PEtab.PEtabObservable}, measurements::DataFrame, parameters::Vector{PEtab.PEtabParameter}, speciemap::Nothing, parametermap::Nothing, events::Nothing, verbose::Bool)
@ PEtab C:\Users\lange.julia\packages\PEtab\MalGN\src\julia_input\petab_model.jl:42
[3] #PEtabModel#127
@ C:\Users\lange.julia\packages\PEtab\MalGN\src\julia_input\petab_model.jl:13 [inlined]
[4] PEtabModel
@ C:\Users\lange.julia\packages\PEtab\MalGN\src\julia_input\petab_model.jl:1 [inlined]

I assume this is the part of the code it's trying to call:

PEtab.jl/src/julia_input/maps.jl

Lines 43 to 48 in 897e821

	if !(value isa Real)
	throw(PEtabInuptError("When setting a parameter to a fixed value in the " *
	"model system it must be set to a constant " *
	"numberic value. This does not hold for $pid " *
	"which is set to $value"))
	end

If the parameter values are checked to be reals does that mean the new method is not supported yet?

[sebapersson]

Can you provide a MVE? This can be supported by PEtab.jl, the error is rather thrown because I did not have this in mind when writing the code, so I need to check what happens internally.
Also, does everything work if you run structurally_simplify?

[langestefan]

I was able to bring it this far:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

using PEtab
using Distributions
using OrdinaryDiffEq
using DataFrames
using Optim
using DataInterpolations


# create some test data to interpolate
T_a_interp = DataInterpolations.AkimaInterpolation(vec(1:0.5:10.5), rand(1:0.1:10, 20))

# Define the ODESystem
@mtkmodel SYS begin
    @parameters begin
        c1
        c2
        c3 = 1.0
    end

    # input functions
    @parameters begin
        (T_a_f)(..) = T_a_interp
    end

    @variables begin
        S(t) = 10.0
        E(t) = 50.0
        SE(t) = 0.0
        P(t) = 0.0
    end

    @equations begin
        T_a ~ T_a_f(t)
        D(S) ~ -c1 * S * E + c2 * SE
        D(E) ~ -c1 * S * E + c2 * SE + c3 * SE
        D(SE) ~ c1 * S * E - c2 * SE - c3 * SE
        D(P) ~ c3 * SE + T_a
    end
end
@mtkbuild sys = SYS()

# Observables
@parameters sigma
@unpack E, S, P = sys

obs_sum = PEtabObservable(S + E, 3.0)
obs_p = PEtabObservable(P, sigma)

observables = Dict("obs_p" => obs_p, "obs_sum" => obs_sum)

# Parameters
p_c1 = PEtabParameter(:c1)
p_c2 = PEtabParameter(:c2; prior=Normal(10.0, 0.3))
p_sigma = PEtabParameter(:sigma)

pest = [p_c1, p_c2, p_sigma]

# Measurement data
# Simulate with 'true' parameters
ps = Dict(:c1 => 1.0, :c2 => 10.0, :c3 => 1.0)
u0 = Dict(:S => 100.0, :E => 50.0, :SE => 0.0, :P => 0.0)

tspan = (0.0, 10.0)
oprob = ODEProblem(sys, u0, tspan, ps)
sol = solve(oprob, Rodas5P(); saveat=0:0.5:10.0)

obs_sum = (sol[:S] + sol[:E]) .+ randn(length(sol[:E]))
obs_p = sol[:P] .+ randn(length(sol[:P]))
df_sum = DataFrame(obs_id="obs_sum", time=sol.t, measurement=obs_sum)
df_p = DataFrame(obs_id="obs_p", time=sol.t, measurement=obs_p)

measurements = vcat(df_sum, df_p)

# Create the PEtab model
model_sys = PEtabModel(sys, observables, measurements, pest)
petab_prob = PEtabODEProblem(model_sys)

# solve 
x0 = get_startguesses(petab_prob, 1)
res = calibrate(petab_prob, x0, IPNewton())

This gives:

ERROR: LoadError: AssertionError: Multiple independent variables are used in the model

From line 16. Unfortunately I am not knowledgeable enough about MTK to bring this any further without docs on how to include the interpolation using this method.