@mtkmodel ODE Confusion, DiffEqParamEstim

Specific Domains Modelling & Simulations
1

Apologies, I think I have a basic misunderstanding somewhere, I’ve been stuck for a few days. I’m struggling to understand how I should use @mtkmodel defined ODEs, I ran into the problem trying to use DiffEqParamEstim. I’ve simplified everything down to just exponential growth:

@mtkmodel ExponentialGrowth begin
    @variables begin
        N(t) = 1.0
    end
    @parameters begin
        r = 1.0
    end
    @equations begin
        D(N) ~ r * N
    end
end

And as an example I’m trying to fit the r parameter to data:

function fit1(df::DataFrame)
    t_data = df.t
    N_data = reshape(df.N, 1, :)
    @mtkbuild sys = ExponentialGrowth()
    prob = ODEProblem(sys, [sys.N => 1.0], (0.0, 5.0), [sys.r => 1.0])
    cost_function = build_loss_objective(prob, Tsit5(), L2Loss(t_data, N_data),
                                         Optimization.AutoForwardDiff())
    optprob = Optimization.OptimizationProblem(cost_function, [sys.r => 1.0])
    optsol = solve(optprob, BFGS())
end

fit1(DataFrame(t = 1:5, N = exp.(1:5))) errors differently depending on the parameter type, as written above I get:

ERROR: MethodError: no method matching one(::Type{Pair{Num, Float64}})
The function `one` exists, but no method is defined for this combination of argument types.

If I replace the symbolic maps with normal vectors in the ODEProblem and OptimizationProblem I get:

ERROR: MethodError: no method matching getindex(::Nothing, ::Int64)
The function `getindex` exists, but no method is defined for this combination of argument types.

Digging in a bit I notice that I can instantiate the ODEProblem with plain vectors or symbolic maps and it solves fine, but when remaking the problem (which I guess is what’s happening in the optimisation) it errors when providing a plain vector, but only to p not u0:

prob = ODEProblem(sys, [1.0], (0.0, 10.0), [1.0])
solve(prob) # works
remake(prob, u0 = [1.0]) # works
remake(prob, u0 = [sys.N => 1.0]) # works
remake(prob, u0 = [1], p = [sys.r => 2.0]) # works
remake(prob, p  = [1.0]) # fails
ERROR: MethodError: no method matching getindex(::Nothing, ::Int64)
The function `getindex` exists, but no method is defined for this combination of argument types.

Which demonstrates to me that I don’t understand how the relationship between ODEProblems and symbolic parameters is supposed to work.

Eventually I found a version that runs using prob_generator and setting lb and ub (without lb/ub I get Success but bogus values like -80000):

function fit(df::DataFrame)
    t_data = df.t
    N_data = reshape(df.N, 1, :)
    @mtkbuild sys = ExponentialGrowth()
    prob = ODEProblem(sys, [1.0], (0.0, 5.0), [1.0])
    cost_function = build_loss_objective(prob, Tsit5(), L2Loss(t_data, N_data),
                                         Optimization.AutoForwardDiff(),
                                         prob_generator = (prob, p) -> remake(prob, p = [sys.r => p[1]]))
    optprob = Optimization.OptimizationProblem(cost_function, [1.0], lb = [-10.0], ub = [1000.0])
    optsol = solve(optprob, BFGS())
end

But it still feels like I’m doing something wrong, or is it just that MTK and DiffEqParamEstim are not meant to be immediately composable (and hence why have to write the prob_generator to connect them)? Also, this version takes several seconds, whereas a version not using MTK takes less than a milisecond:

function fit3(df::DataFrame)
    t_data = df.t
    N_data = reshape(df.N, 1, :)
    expgrow(u, p, t) =  u * p[1]
    prob = ODEProblem(expgrow, 1.0, (0.0, 5.0), [1.0])
    cost_function = build_loss_objective(prob, Tsit5(), L2Loss(t_data, N_data),
                                         Optimization.AutoForwardDiff())
    optprob = Optimization.OptimizationProblem(cost_function, [1.0])
    optsol = solve(optprob, BFGS())
end

(CancerCellsModels) julia> @time fit(df)
┌ Warning: At t=2.3962119558651183, dt was forced below floating point epsilon 4.440892098500626e-16, and step error estimate = NaN. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase ~/.julia/packages/SciMLBase/JxpZc/src/integrator_interface.jl:623
  3.516851 seconds (15.59 M allocations: 791.920 MiB, 9.05% gc time)
retcode: Success
u: 1-element Vector{Float64}:
 1.2546681515725144

(CancerCellsModels) julia> @time fit3(df)
  0.000731 seconds (6.10 k allocations: 448.328 KiB)
retcode: Success
u: 1-element Vector{Float64}:
 1.0000024202087792

and this doesn’t change if I pull the @mtkbuild outside the function and pass in the sys.

Julia 1.11.4
MTK 9.82.0
DiffEqParamEstim 2.2.0