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:
# 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:
# 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.