# DDE parameter estimation using DiffEqFlux

**URL:** <https://discourse.julialang.org/t/dde-parameter-estimation-using-diffeqflux/60399>\
**Category:** New to Julia\
**Tags:** diffeq, flux\
**Created:** [May 1, 2021, 8:16pm UTC](https://discourse.julialang.org/t/dde-parameter-estimation-using-diffeqflux/60399 "2021-05-01T20:16:36Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![ANasc](https://avatars.discourse-cdn.com/v4/letter/a/41988e/32.png) [@ANasc](https://discourse.julialang.org/u/ANasc)\
**Post date:** [May 1, 2021, 8:16pm UTC](https://discourse.julialang.org/t/dde-parameter-estimation-using-diffeqflux/60399/1 "2021-05-01T20:16:36Z")

</div>

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/](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](https://diffeqparamestim.sciml.ai/stable/tutorials/ODE_inference/) 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:

```julia
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)

```

---

<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:** [May 2, 2021, 12:42am UTC](https://discourse.julialang.org/t/dde-parameter-estimation-using-diffeqflux/60399/2 "2021-05-02T00:42:03Z")

</div>

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/](https://diffeqflux.sciml.ai/dev/examples/delay_diffeq/)

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

```julia
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.

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![ANasc](https://avatars.discourse-cdn.com/v4/letter/a/41988e/32.png) [@ANasc](https://discourse.julialang.org/u/ANasc)\
**Post date:** [May 2, 2021, 4:37am UTC](https://discourse.julialang.org/t/dde-parameter-estimation-using-diffeqflux/60399/3 "2021-05-02T04:37:32Z")

</div>

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?

---

<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:** [May 2, 2021, 5:02am UTC](https://discourse.julialang.org/t/dde-parameter-estimation-using-diffeqflux/60399/4 "2021-05-02T05:02:45Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ANasc](https://avatars.discourse-cdn.com/v4/letter/a/41988e/32.png) [@ANasc](https://discourse.julialang.org/u/ANasc)\
**Post date:** [May 2, 2021, 5:47am UTC](https://discourse.julialang.org/t/dde-parameter-estimation-using-diffeqflux/60399/5 "2021-05-02T05:47:27Z")

</div>

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