# Syntax for custom cost function for ODE parameter estimation

**URL:** <https://discourse.julialang.org/t/syntax-for-custom-cost-function-for-ode-parameter-estimation/43467>\
**Category:** Optimization (Mathematical)\
**Tags:** diffeq, optim\
**Created:** [July 22, 2020, 2:47am UTC](https://discourse.julialang.org/t/syntax-for-custom-cost-function-for-ode-parameter-estimation/43467 "2020-07-22T02:47:54Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![vaglino](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaglino/32/16669_2.png) [@vaglino](https://discourse.julialang.org/u/vaglino)\
**Post date:** [July 22, 2020, 2:47am UTC](https://discourse.julialang.org/t/syntax-for-custom-cost-function-for-ode-parameter-estimation/43467/1 "2020-07-22T02:47:55Z")

</div>

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:

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

---

<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:** [July 22, 2020, 12:20pm UTC](https://discourse.julialang.org/t/syntax-for-custom-cost-function-for-ode-parameter-estimation/43467/2 "2020-07-22T12:20:26Z")

</div>

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

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

[https://diffeqflux.sciml.ai/dev/examples/optimization\_ode/](https://diffeqflux.sciml.ai/dev/examples/optimization_ode/)

---

<div class="post-metadata">

**Author:** ![vaglino](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaglino/32/16669_2.png) [@vaglino](https://discourse.julialang.org/u/vaglino)\
**Post date:** [August 1, 2020, 9:16am UTC](https://discourse.julialang.org/t/syntax-for-custom-cost-function-for-ode-parameter-estimation/43467/3 "2020-08-01T09:16:25Z")

</div>

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 ).

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

```

---

<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 1, 2020, 12:04pm UTC](https://discourse.julialang.org/t/syntax-for-custom-cost-function-for-ode-parameter-estimation/43467/4 "2020-08-01T12:04:55Z")

</div>

> [@vaglino](#):
>
> since the loss has one argument ( p ).

It is not. See something like [https://diffeqflux.sciml.ai/dev/examples/mnist\_neural\_ode/](https://diffeqflux.sciml.ai/dev/examples/mnist_neural_ode/) for how the data iterator works.