DDE parameter estimation using DiffEqFlux

New to Julia
,
1

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

The error message was saying you needed to switch the sensealg, as is shown in the DDE tutorials. If you scrolled down a bit more on that side bar you’ll see there are DDE examples:

https://diffeqflux.sciml.ai/dev/examples/delay_diffeq/

For your case, here’s a working version, just swapping the sensealg:

using Flux,DiffEqFlux,Plots,DifferentialEquations,DiffEqSensitivity

#***************************************************************
#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(autodiff=false))
p = [1500.0]

function loss(p)
  prob_ = DyEstimateS(0.0,365.0*5.0,p)
  sol = solve(prob_, alg, p=p, saveat = tsteps, sensealg = ReverseDiffAdjoint())
  loss = sum(abs2,data[i] - sol[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)

and the faster thing would be to just use forward-mode, given the small number of parameters.

using GalacticOptim
optf = GalacticOptim.OptimizationFunction((x, p) -> loss(x), GalacticOptim.AutoForwardDiff())
optprob = GalacticOptim.OptimizationProblem(optf, p)

result_dde = GalacticOptim.solve(optprob,
                                       ADAM(0.1),
                                       cb = callback,
                                       maxiters = 300)
3

Ah, I can’t believe I missed the DDE tutorial. I googled all over but somehow missed that. Sorry about that, but thanks for giving such a great answer! It works for me now, although the minimizer doesn’t seem to change at all from one iteration to the next. The resulting estimated parameter ends up being almost exactly the initial guess. Any tips on troubleshooting this?

4

The one I posted was converging really slowly, so it needs a few orders of magnitude more iterations. Or use BFGS.

5

Yes, BFGS() works perfectly and is extremely fast. Thank you so much!