Hi,
I’m trying to estimate a parameter in a DDE model using DiffEqFlux. I’ve tried to adapt the tutorial at https://diffeqflux.sciml.ai/dev/examples/optimization_ode/ to my problem but have not been able to get it to run. I’ve been able to estimate the parameter no problem using DiffEqParamEstim with the example at Optimization-Based ODE Parameter Estimation · DiffEqParamEstim.jl using the L2Loss() loss function, but I wanted the flexibility to use more complicated loss functions and thought I would try to learn how to do it in DiffEqFlux. I’m hoping that I’ve just made a simple mistake, but perhaps I’ve misunderstood something at a deeper level.
I would appreciate if anyone could lead me towards a solution. I feel like I’ve adapted the tutorial pretty closely, and the function loss(p) does return a sensible result (e.g., if I plot it for a number of parameter values I can visually see a minimum that makes sense). However, I keep getting the MethodError: no method matching similar(::Float64) where I call solve in the loss function.
Thanks for your time.
MWE:
using Flux,DiffEqFlux,Plots,DifferentialEquations
#***************************************************************
#Parameters
const Temp0K = 15+273.15
const k = 8.62 * 10^(-5)
const H = 0.0987*1.085
const r = 0.238;
const b_H = 1/(15*365)
const mu_P = 1/(365*7.5)
const mu_M = 0.0047
const tauM = 177.01
const δ = 2.454e-4
const Iburden = 0.7687
const α = 0.0
#Data
data = [2.17, 0.86, 0.0, 0.0, 49.29, 20.0, 18.6, 23.6, 61.5, 0.0, 177.0, 142.0, 0.46, 1.37, 0.0, 0.0, 29.47, 0.0, 151.7, 1.82,
3.8, 5.4, 6.93, 311.8, 194.73, 261.22, 60.41, 0.0, 0.0, 5.67, 36.44, 151.6, 4.05, 251.4, 0.0, 82.5, 1.54, 22.08, 91.62,
182.0, 69.8, 225.8, 172.8, 130.0, 159.61, 5.54, 24.1, 116.15, 5.77, 10.39, 376.0, 4.72, 5.11, 6.67, 0.44, 4.84, 0.0, 20.0, 4.77, 0.0, 0.0, 0.8, 13.46, 3.95, 6.37, 0.39, 14.85, 0.0, 3.0, 9.83, 7.31, 231.51, 7.58, 0.0, 22.04, 0.19, 0.0, 0.0,
178.51, 210.78, 176.83, 10.53, 28.45, 2.8, 234.23, 0.0, 3.7]
t = [365*ones(6,1);365*2*ones(4,1);365*3*ones(6,1);365*4*ones(15,1);365*5*ones(56,1)]
scatter(t,data,markershape = :circle)
#Parameter estimation
function est_S(du,u,h,p,t)
#State variables
M = u[1]; P = u[2]; DM1 = u[3];
# History calculations
ucache = h(p, t - tauM)
P_past = ucache[2]
#Prelim calculation
Ep = (P/H)^2 * (r+1)/r + (P/H);
Eplag = (P_past/H)^2 * (r+1)/r + (P_past/H)
DM2 = exp(-mu_M * tauM) * exp(-b_H * tauM)
matI = H * δ * p[1] * Iburden
matM = H * δ * p[1] * Iburden *DM1 * DM2
#Model
du[1] = matI - M * (mu_M + b_H) - matM
du[2] = matM - P * (mu_P + b_H) - α * H * Ep
du[3] = DM1 * (α * P_past / H - α * P / H)
nothing
end
#Set up problem
function DyEstimateS(tstart, tfin, p)
DM2 = exp(-mu_M * tauM) * exp(-b_H * tauM)
u0 = [0.0,0.0,DM2]
h(p,t) = [0.0, 0.0, DM2];
tspan = (tstart, tfin);
probDynSea = DDEProblem(est_S, u0, h, tspan, p;constant_lags = [tauM],save_idxs = 2)
end
tsteps = 1.0:1.0:1825.0
alg = MethodOfSteps(Rosenbrock23())
p = [1500.0]
function loss(p)
prob_ = DyEstimateS(0.0,365.0*5.0,p)
sol = solve(prob_, alg, p=p, saveat = tsteps)
loss = sum(abs2,data[i] - sol(t[i]) for i in 1:length(t))
return loss, sol
end
callback = function (p, l, pred)
display(l)
plt = plot(pred, ylim = (0, 6))
display(plt)
# Tell sciml_train to not halt the optimization. If return true, then
# optimization stops.
return false
end
result_ode = DiffEqFlux.sciml_train(loss,p,
ADAM(0.1),
cb = callback,
maxiters = 100)