# Parameter estimation ODE with custom Loss function

**URL:** <https://discourse.julialang.org/t/parameter-estimation-ode-with-custom-loss-function/124835>\
**Category:** General Usage\
**Tags:** differentialequation\
**Created:** [January 16, 2025, 2:13pm UTC](https://discourse.julialang.org/t/parameter-estimation-ode-with-custom-loss-function/124835 "2025-01-16T14:13:57Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Marta\_Pardo\_Araujo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marta_pardo_araujo/32/208774_2.png) [@Marta\_Pardo\_Araujo](https://discourse.julialang.org/u/Marta_Pardo_Araujo)\
**Post date:** [January 16, 2025, 2:13pm UTC](https://discourse.julialang.org/t/parameter-estimation-ode-with-custom-loss-function/124835/1 "2025-01-16T14:13:57Z")

</div>

Hi, I am trying to estimate the parameters of my systems of ODE’s. I have done it before using the L2Loss() and it worked. Now I would like to do a custom loss function where I take the mean of the solutions of the Odes over a period of time, summer, and compare it to the yearly observations. I have made my own loss function which take the parameters p integrates the ODE and process it to compute my custom loss and it works as expected but the problem is with the solve().  
This is the code:

```julia
# Non autonomous sigmoidal model ---------------------------------------------------------
function fun_na!(du, u, p, t)
  R_M_vec = [interpolated_functions[i](t) for i in 1:N] # Time series RM
  eta1 = p[1]./(1.0 .+ exp.(-p[2]*m_c.*eta.+p[3]))
  #du.= R_M_vec.*(p[1].*mat*u .+ p[2].*(m_c*eta)*u) .* (1 .- u) - p[3] .* u
  du.= R_M_vec.*(eta1*u) .* (1 .- u) - p[4] .* u
  nothing
end

# Set initial parameters ------------------------------------------------------------------
t0= t_obs[1]
tf= t_obs[end] + 10
tspan = (t0, tf)
t_vect=1:tf
u0 = pop_init
p = [0.001,0.000000001,5, 0.00000001] 

# Create the ode model to solve it ---------------------------------------------------------
# using BenchmarkTools
prob = DifferentialEquations.ODEProblem(fun_na!, u0, tspan, p)
alg= TRBDF2()
# sol = solve(prob,alg;p=p)
# plot(sol)

# Create the cost function -----------------------------------------------------------------
# Load libraries
using ForwardDiff, OptimizationOptimJL, OptimizationBBO

# Create vector of summer times per year
vec_year = unique(year.(R_M.date))[2:end]

# Function to extract times related to summer per year
function summer_times_by_year(R_M)
    summer_years = vec_year
    summer_times = Dict{Int, Vector{Int}}()
    for y in summer_years
        summer_times[y] = filter(row -> year(row.date) == y && month(row.date) in 6:8, R_M).time
    end
    return summer_times
end

# Run through RM
summer_t_obs_by_year = summer_times_by_year(R_M)

# Function to compute average probability of occupancy in summer per year and comarca
function average_summer_solution_by_year(sol)
    sol_mat = zeros(Float64,N,length(vec_year))
    j = 1
    for years in vec_year
        times = summer_t_obs_by_year[years]
        vec_sol = Vector{Float64}(undef,N)
        for i in 1:N
            vec_sol[i] = mean(sol(times, idxs = i))
        end  
        sol_mat[:,j] .= vec_sol 
        j = j+1 
    end
    return sol_mat
end

# Loss function
function summer_loss_by_year(p)
    sol = solve(prob, alg;p=p)
    summer_avg_by_year = average_summer_solution_by_year(sol)
    total_loss = sum((summer_avg_by_year .- matrix_obs).^2)
    return total_loss
end

# Test
p = [0.00047422389966972274 ,0.15014600413485327 ,9.546705822176884 ,1.3723893836256143e-7]
summer_loss_by_year(p) # Test

# Define the OptimizationFunction
opt_func = OptimizationFunction((p, _) -> summer_loss_by_year(p), Optimization.AutoForwardDiff())
optprob = OptimizationProblem(opt_func, init_param, lb = low_bound, ub = up_bound)
optsol = solve(optprob, BFGS(); maxiters = 1000)

```

And I get this error:\>

> ERROR: InterruptException:  
> Stacktrace:  
> [1] try\_yieldto(undo::typeof(Base.ensure\_rescheduled))  
> @ Base ./task.jl:931  
> [2] wait()  
> @ Base ./task.jl:995  
> [3] uv\_write(s::Base.TTY, p::Ptr{UInt8}, n::UInt64)  
> @ Base ./stream.jl:1048  
> [4] unsafe\_write(s::Base.TTY, p::Ptr{UInt8}, n::UInt64)  
> @ Base ./stream.jl:1120  
> [5] unsafe\_write  
> @ ./io.jl:431 [inlined]  
> [6] write  
> @ ./strings/io.jl:248 [inlined]  
> [7] print  
> @ ./strings/io.jl:250 [inlined]  
> [8] with\_output\_color(f::Function, color::Symbol, io::IOContext{Base.TTY}, args::Module; bold::Bool, italic::Bool, underline::Bool, blink::Bool, reverse::Bool, hidden::Bool)  
> @ Base ./util.jl:112  
> [9] with\_output\_color  
> @ ./util.jl:73 [inlined]  
> [10] printstyled  
> @ ./util.jl:141 [inlined]  
> [11] print\_module\_path\_file(io::IOContext{Base.TTY}, modul::Module, file::String, line::Int64; modulecolor::Symbol, digit\_align\_width::Int64)  
> @ Base ./errorshow.jl:761  
> [12] print\_module\_path\_file  
> @ ./errorshow.jl:755 [inlined]  
> [13] print\_stackframe(io::IOContext{Base.TTY}, i::Int64, frame::Base.StackTraces.StackFrame, n::Int64, ndigits\_max::Int64, modulecolor::Symbol)  
> @ Base ./errorshow.jl:749  
> [14] print\_stackframe(io::IOContext{Base.TTY}, i::Int64, frame::Base.StackTraces.StackFrame, n::Int64, ndigits\_max::Int64, modulecolordict::IdDict{Module, Symbol}, modulecolorcycler::Base.Iterators.Stateful{Base.Iterators.Cycle{Vector{Symbol}}, Union{Nothing, Tuple{Symbol, Int64}}, Int64})  
> @ Base ./errorshow.jl:706  
> [15] show\_full\_backtrace(io::IOContext{Base.TTY}, trace::Vector{Any}; print\_linebreaks::Bool)  
> @ Base ./errorshow.jl:605  
> [16] show\_full\_backtrace  
> @ ./errorshow.jl:598 [inlined]  
> [17] show\_backtrace(io::IOContext{Base.TTY}, t::Vector{Base.StackTraces.StackFrame})  
> @ Base ./errorshow.jl:802  
> [18] showerror(io::IOContext{Base.TTY}, ex::MethodError, bt::Vector{Base.StackTraces.StackFrame}; backtrace::Bool)  
> @ Base ./errorshow.jl:97  
> [19] showerror(io::IOContext{Base.TTY}, ex::MethodError, bt::Vector{Base.StackTraces.StackFrame})  
> @ Base ./errorshow.jl:93  
> [20] display\_repl\_error(io::Base.TTY, stack::VSCodeServer.EvalErrorStack; unwrap::Bool)  
> @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/repl.jl:261  
> [21] (::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()  
> @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/eval.jl:231  
> [22] withpath(f::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)  
> @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/repl.jl:276  
> [23] (::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()  
> @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/eval.jl:179  
> [24] hideprompt(f::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})  
> @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/repl.jl:38  
> [25] (::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()  
> @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/eval.jl:150  
> [26] with\_logstate(f::Function, logstate::Any)  
> @ Base.CoreLogging ./logging.jl:515  
> [27] with\_logger  
> @ ./logging.jl:627 [inlined]  
> [28] (::VSCodeServer.var"#66#71"{VSCodeServer.ReplRunCodeRequestParams})()  
> @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/eval.jl:263  
> [29] #invokelatest#2  
> @ ./essentials.jl:892 [inlined]  
> [30] invokelatest(::Any)  
> @ Base ./essentials.jl:889  
> [31] (::VSCodeServer.var"#64#65")()  
> @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/eval.jl:34

Looking at different answers about this in the Julia community I have tried another package DiffEqFlux:

```julia
# dfferent algorithm opt_func
using DiffEqFlux
res = DiffEqFlux.sciml_train(summer_loss_by_year,init_param)

```

And this runs with a warning that there is a problem with precompiling.

---

<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:** [January 21, 2025, 9:33am UTC](https://discourse.julialang.org/t/parameter-estimation-ode-with-custom-loss-function/124835/2 "2025-01-21T09:33:56Z")

</div>

Generally for custom loss functions we would recommend using SciMLSensitivity.jl directly. See the documentation on solving inverse problems for some information on how to piece things together:

> **[Parameter Estimation of Ordinary Differential Equations · SciMLSensitivity.jl](https://docs.sciml.ai/SciMLSensitivity/stable/tutorials/parameter_estimation_ode/)**
>
> Documentation for SciMLSensitivity.jl.

I recently ran a workshop on inverse problems and the lecture notes from that might be useful:

[https://sciml.github.io/GeiloInverseProblemWorkshop/Lecture1\_Solving\_Inverse\_Problems\_in\_Julia](https://sciml.github.io/GeiloInverseProblemWorkshop/Lecture1_Solving_Inverse_Problems_in_Julia)

---

<div class="post-metadata">

**Author:** ![Marta\_Pardo\_Araujo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marta_pardo_araujo/32/208774_2.png) [@Marta\_Pardo\_Araujo](https://discourse.julialang.org/u/Marta_Pardo_Araujo)\
**Post date:** [January 24, 2025, 3:11pm UTC](https://discourse.julialang.org/t/parameter-estimation-ode-with-custom-loss-function/124835/3 "2025-01-24T15:11:31Z")

</div>

Thanks a lot I will take a look at it.  
I have follow you tutorial on youtube of DifferetialEquation.jl which was really helpful.