Syntax for custom cost function for ODE parameter estimation

Specific Domains Optimization (Mathematical)
,
1

Hi, I’m trying to write my own cost function for parameter estimation of an ode model.
I am trying to weight the loss by the time t.

Based on the documentation I came up with the following, but I’m unsure about the correct syntax:

function optimization(model,data,constants,dens,options)
# fit function - given model, data, densities, options get best parameters
    @unpack lb,ub,tspan = options
    
    f = (du,u,p,t) -> model(du,u,10 .^ p,t,dens,constants) # optimize in log space
    prob = ODEProblem(f,u₀,tspan,p)
    # cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data),
    #                maxiters=10000,verbose=false)
    cost_function = build_loss_objective(prob,Tsit5(),my_loss_function(sol,t,data),
                    maxiters=10000,verbose=false)

    # lb = [-10., -10.] # from 1e-10 to 1e10
    # ub = [5., 5.]
    result = optimize(cost_function,lb,ub,[0.,0.], Fminbox(BFGS()))
end

function my_loss_function(sol,t,data)
   tot_loss = 0.0
   if any((s.retcode != :Success for s in sol))
     tot_loss = Inf
   else
     # calculatisolon for the loss here
     # n = length(t)
     exp = sol(t)
     tot_loss = sum( ((exp .- data) .^ 2) .* t ) # time weighted loss
   end
   tot_loss
end

fit = optimization(model_name,data,constants,dens,c)

Thanks

2

Use DiffEqFlux for this: https://diffeqflux.sciml.ai/dev/. Specifically, this tutorial should be helpful:

https://diffeqflux.sciml.ai/dev/examples/optimization_ode/

3

Thanks! I’ve rewritten it in the DiffEqFlux style, and it works fine when the dataset is a global variable that the loss function can access, but I am now struggling to pass in a given dataset from my optimization function to the loss function, since the loss has one argument ( p ).

using DifferentialEquations, Flux, Optim, DiffEqFlux, DiffEqSensitivity, Plots

# -------------------------------------------------------------------------
function bimolecular(du,u,p,t,dens)
    # unpack rates and constants
    nᵣ = u[1]
    k₁,k₋₁  = p
    mᵣ,mₗ,A = dens
    # model
    du[1] = dnᵣ = A*k₁*mᵣ*mₗ - k₋₁*nᵣ
end

# --------------------------create dummy data-----------------------------------
f = (du,u,p,t) -> bimolecular(du,u,p,t,dens) # enclose constants

# u0 = [0.0;0.0]
u₀ = [0.0]          # Initial condition
tspan = (0.0,10.0)  # Simulation interval
p = [1e-4,0.5]      # equation parameter. p = [k₁, k₋₁]
dens = [100., 10., 1.]

# Setup the ODE problem, then solve
prob = ODEProblem(f,u₀,tspan,p)
sol = solve(prob)
# sol = solve(prob,Tsit5())
h = plot(sol)

t = collect(range(0,stop=10,length=101))
using RecursiveArrayTools # for VectorOfArray
randomized = VectorOfArray([(sol(t[i]) + .02randn(size(sol,1))) for i in 1:length(t)])
dataset1 = convert(Array,randomized)
scatter!(t,dataset1')


function loss(p) # need to use data that was passed to the optimization function
    tmp_prob = remake(prob,p=p)
    tmp_sol = solve(tmp_prob,saveat=0.1)
    sum(abs2,Array(tmp_sol) - data), tmp_sol
end

function callback(p,l,tmp_sol)
    @show l
    tmp_prob = remake(prob,p=p)
    tmp_sol = solve(tmp_prob,saveat=0.1)
    fig = plot(tmp_sol)
    scatter!(fig,tmp_sol.t,data')
    display(fig)
    false
end



function optimization(model,data,dens,options)
# fit function - given model, data, densities, options get best parameters
    
    lb,ub = options
   
    f = (du,u,p,t) -> model(du,u,p,t,dens) 
    prob = ODEProblem(f,u₀,tspan,p)
    
    res = DiffEqFlux.sciml_train(loss,pinit,ADAM(0.01),cb=callback,maxiters=100)
    res.minimizer
    res = DiffEqFlux.sciml_train(loss,res.minimizer,BFGS(),cb=callback)
end

model_opt2 = (1.0e-6,10.0)
pinit = [2.0e-4,1.0]

# want to use dataset1 for optimization
optimization(bimolecular,dataset1,dens,model_opt2)
4

It is not. See something like https://diffeqflux.sciml.ai/dev/examples/mnist_neural_ode/ for how the data iterator works.