I am trying to estimate the parameters of a ODE model using a frequentist approach. While the simple L2Loss approach it works after when I run a more complicated one it does not. It gives me the following error:
ERROR: InexactError: Int64(4.8)
Stacktrace:
[1] Int64
@ ./float.jl:912 [inlined]
[2
] (::DiffEqParamEstim.var"#43#48"{…})(p::Vector{…}, ::SciMLBase.NullParameters)
@ DiffEqParamEstim ~/.julia/packages/DiffEqParamEstim/6z5Ou/src/multiple_shooting_objective.jl:33
[3] (::OptimizationFunction{…})(::Vector{…}, ::Vararg{…})
@ SciMLBase ~/.julia/packages/SciMLBase/rR75x/src/scimlfunctions.jl:3763
[4] (::OptimizationBBO.var"#19#21"{OptimizationCache{…}})(θ::Vector{Float64})
@ OptimizationBBO ~/.julia/packages/OptimizationBBO/UgujW/src/OptimizationBBO.jl:147
[5] fitness(x::Vector{…}, p::BlackBoxOptim.FunctionBasedProblem{…})
@ BlackBoxOptim ~/.julia/packages/BlackBoxOptim/lZtsr/src/problem.jl:61
[6] setup_problem(func::Function, parameters::BlackBoxOptim.ParamsDictChain)
@ BlackBoxOptim ~/.julia/packages/BlackBoxOptim/lZtsr/src/bboptimize.jl:37
[7] bbsetup(functionOrProblem::Function, parameters::Dict{…}; kwargs::@Kwargs{…})
@ BlackBoxOptim ~/.julia/packages/BlackBoxOptim/lZtsr/src/bboptimize.jl:111
[8] bbsetup
@ ~/.julia/packages/BlackBoxOptim/lZtsr/src/bboptimize.jl:109 [inlined]
[9] __solve(cache::OptimizationCache{…})
@ OptimizationBBO ~/.julia/packages/OptimizationBBO/UgujW/src/OptimizationBBO.jl:167
[10] solve!(cache::OptimizationCache{…})
@ SciMLBase ~/.julia/packages/SciMLBase/rR75x/src/solve.jl:188
[11] solve(::OptimizationProblem{…}, ::BBO_adaptive_de_rand_1_bin_radiuslimited; kwargs::@Kwargs{…})
@ SciMLBase ~/.julia/packages/SciMLBase/rR75x/src/solve.jl:96
[12] top-level scope
@ ~/albo_mobility/code/Hanski/test/test_freq_simpMod_N5 .jl:184
Some type information was truncated. Useshow(err)to see complete types.
This is the code with fake data input:
# Code that integrates the Hanski model and estimate the parameters of the model from the presence absence data frequentist algorithm
# It is a test with fake observed data and the simplest model
# Load pkgs and remove data -----------------------------------------------------
import Pkg
using DifferentialEquations
using DataFrames
using CSV
using Plots
using Shapefile
import GeoDataFrames as GDF
using LinearAlgebra
using ODE
using Interpolations
using DiffEqParamEstim
using SparseArrays
using Optimization
# Constant variables
const m_c = 0.001 # probability of mosquito in a car
const alp = 1/200 # average natural dispersal distance
# Load data input
const end_ind = 5
eta = [0.0 7.74415 5.24014 8.57707 5.21319;
6.9175 0.0 7.23771 45.2212 16.8025;
6.17319 7.426 0.0 17.7236 5.49988;
8.40241 44.21 20.1597 0.0 5309.46;
5.8309 17.5283 5.36841 5226.36 0.0]
const N = size(eta, 1) # Number of patches
mat = [ 0.0 0.169023 0.157908 0.0503147 0.0531469
0.169023 0.0 0.357034 0.00273453 0.27251
0.157908 0.357034 0.0 0.236906 0.00269576
0.0503147 0.273453 0.236906 0.0 0.830696
0.0531469 0.0027251 0.00269576 0.830696 0.0]
# Initial conditions
pop_init = zeros(N)
pop_init[1] = 1
# Load R_M and interpolated
function interpolate_RM!(loc, interpolated_functions)
# Read RM data time series
R_M = CSV.read(loc*"albo_mobility/data/InputHanski/R_M_mitma.csv",
DataFrame)[1:end,1:(end_ind+3)]
# Select columns from the second to the second-to-last column, excluding the third column
R_M = R_M[:, Not([1, 3])]
# Interpolate R_M --------------------------------------------------------------------
dates = R_M[:, 1]
R_Ms = Matrix(R_M[:, 2:end])
# Perform interpolation for each location
for i in 1:size(R_Ms, 2)
# Extract temperature values for the current location
R_Ms_val = R_Ms[:, i]
dates_num = 1:size(dates, 1)
# Perform linear interpolation
itp = LinearInterpolation(dates_num, R_Ms_val,extrapolation_bc=Flat())
# Store the interpolated function
push!(interpolated_functions, itp)
end
end
# Create an array to store interpolated functions
interpolated_functions = []
interpolate_RM!(loc, interpolated_functions)
# Integration non autonomous
function fun_na!(du, u, p, t)
R_M_vec = [interpolated_functions[i](t) for i in 1:N] # Time series RM
#eta1 = 1.0./(1.0 .+ exp.(-p[4].*eta.+p[5]))
Cd = p[1].*mat*u # Natural dispersal
Ch = p[2].*(m_c*eta)*u # Human mobility
#Ch = p[2].*eta1*u # Human mobility
du.= R_M_vec.*(Cd .+ Ch) .* (1 .- u) - p[3] .* u
nothing
end
# Set parameters
t0=0.0
tf= 365.0
tspan = (t0, tf)
t_vect=1:tf
u0 = pop_init
p = [0.026, 0.071, 0.01]
# Create the ode model
prob = DifferentialEquations.ODEProblem(fun_na!, u0, tspan, p)
alg= DP8() #QNDF() #For low tolerances, Rodas4()
obs = solve(prob, alg;p=p)
plot(obs)
# Fake observed data
t = collect(range(0, stop = tf, length = 40))
using RecursiveArrayTools # for VectorOfArray
randomized = VectorOfArray([(obs(t[i]) ) for i in 1:length(t)])
data = convert(Array, randomized)
# Frequentist optimization -------------------------------------------------
using ForwardDiff, OptimizationOptimJL, OptimizationBBO
# Compute the cost function
@time cost_function = build_loss_objective(prob, Tsit5(), L2Loss(t, data),
Optimization.AutoForwardDiff(), maxiters = 10000, verbose = false)
# Optimization
low_bound = [0,0,0]
up_bound = [0.1,0.5,0.5]
optprob = Optimization.OptimizationProblem(cost_function, [0.0,0.0,0.0],
lb = low_bound, ub = up_bound)
@time optsol = solve(optprob, BFGS())
print(optsol)
# Check the solutions for the optimization
sol_opt = solve(prob, alg;p=optsol.u)
plot(obs)
plot!(sol_opt, linestyle = :dash)
# Optimization multiple shooting
bound = Tuple{Float64, Float64}[(0.0, 0.1), (0.0, 0.5), (0.0, 0.5),
(0.0, 0.1), (0.0, 0.5), (0.0, 0.5),
(0.0, 0.1), (0.0, 0.5), (0.0, 0.5),
(0.0, 0.1), (0.0, 0.5), (0.0, 0.5),
(0.0, 0.1), (0.0, 0.5), (0.0, 0.5),
(0.0, 0.1), (0.0, 0.5), (0.0, 0.5),
(0.0, 0.1), (0.0, 0.5), (0.0, 0.5),
(0.0, 0.1), (0.0, 0.5), (0.0, 0.5),
(0.0, 0.1), (0.0, 0.5), (0.0, 0.5)]
ms_obj = multiple_shooting_objective(prob, DP8(), L2Loss(t, data),
Optimization.AutoForwardDiff())
# Run the optimization
optprob = Optimization.OptimizationProblem(ms_obj, fill(0.01,9*3),
lb = first.(bound),ub = last.(bound))
optsol_ms = solve(optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(),
maxiters = 10_000)
# Check the solutions for the optimization multiple shoothing
sol_opt = solve(prob, alg;p=optsol_ms.u[(end - 3):end])
plot(obs)
plot!(sol_opt, linestyle = :dash)
I have also tried with BFGS instead of BBO_adaptive_de_rand_1_bin_radiuslimited() and I get the same error. I do not know where this problem with trying to change a float to a integer comes from.