Declaring state-dependent delay with DDEProblem; simulation gets slower?

Dear All,
I am having a hard time figuring out how to declare dependent_lags that depend
on one of the state variables in my DDE model (using DelayDiffEq.jl). It seems
like it should be easy since the state variable is tracking the length of the
delay as it changes with time. So I feel like it should be

d_lags = ((u,p,t) -> u[8],)
prob = DDEProblem(DynDDE,u0,h,tspan; constant_lags = tauM, dependent_lags = d_lags)

However, when I do this the simulation becomes very slow (usually asks for larger maxiters), and I thought declaring lags was supposed to improve performance.
If I don’t declare the dependent lags, the result is reasonable, fast, and seems to be essentially
identical to when they are declared.

I’m fairly new to Julia so hopefully I’ve just made a rookie mistake that someone
can spot. It’s also possible that I’ve misunderstood exactly what is happening when
lags are declared.

I’ve provided below the simplest working version of the problem as I could without
having to spend a lot of time modifying the original. Sorry it is so long. Briefly,
it models a host-parasite system in which the length of the developmental delay period for one
of the larval stages depends on environmental temperatures, which in turn vary
continuously with time (e.g.,shorter development time in the summer, longer in
the winter). In the original, many of the other parameters are functions of time
as well, and the issue above becomes even more pronounced.

I’d greatly appreciate any help!

using QuadGK
using Roots
using DelayDiffEq
using Plots;

#Parameters
const Temp0K = 15+273.15; const k = 8.62 * 10^(-5); const H = 0.075;const S = 4000000.0
const b_H = 1/5292.5; const b_S = 0.0037; const lambda = 488.0;const mu_W = 0.0186
const Phi_a = 0.024; const Phi_E = 0.46; const mu_P = 1.5e-7; const mu_M = 0.0056
const tauM = 166.82; const mu_U = 0.0001; const mu_I = mu_U; const tauU_a = 25.1;
const tauU_E = 0.87; const tauU_EL = 5.0*tauU_E; const tauU_TL = 282.0218;
const Alpha = 1.4752e-6;const phi_H = 2.7728e-7;const c_K = 263.0; const d_K = 22.6;
const t0 = 105.6;

    #A couple functions that describe how parameters are sensitive to temprature
      #Temperature as a function of time
      Temp(t) = c_K + d_K * sin((t-t0)*2*pi/365);
      TempKstop = 273.15-0.0;

      #Development rate, as a function of Temperature
      function g_U(T)
        if T > TempKstop::Float64
          g_U = 1/(tauU_a * exp((tauU_E/k)*((1/T)-(1/Temp0K))) * (1 + exp((tauU_EL/k)*((1/T)-(1/tauU_TL))))); #Development rate
        else
          g_U = 1/(tauU_a * exp((tauU_E/k)*((1/TempKstop::Float64)-(1/Temp0K))) * (1 + exp((tauU_EL/k)*((1/TempKstop::Float64)-(1/tauU_TL)))));#Development rate if it's below freezing
        end
      end
      gt(t) = g_U(Temp(t))

      #Infection rate, as a function of temperature
      function Phi(T)
        if T > TempKstop::Float64
          Phi = Phi_a * exp((-Phi_E/k)*((1/T)-(1/Temp0K))); #Infection rate
        else
          Phi = Phi_a * exp((-Phi_E/k)*((1/TempKstop::Float64)-(1/Temp0K))); #Infection rate if it's too cold
        end
        return(Phi)
      end

  function DynDDE(du,u,h,p,t)
  #State variables

    #W,U,I,M,and P track different stages of a parasite life-cycle
    #DU1 gives the proportion of U that survive the variabl developmental delay period
    #DM1 give proportion of M that survive constant delay period
    #tauU tracks the length of the developmental delay, which changes continuously with time

    W = u[1]; U = u[2]; I = u[3]; M = u[4]; P = u[5]; DU1 = u[6]; DM1 = u[7]; tauU = u[8];

    #History calculations
    W_past = h(p, t-tauU)[1]; #W at time t-tauU
    U_past = h(p, t-tauU)[2]; #U at time t-tauU
    I_past = h(p, t-tauM)[3]; #I at time t-tauM
    P_past = h(p, t-tauM)[5]; #P at time t-tauM

 #"Constant" Survival functions
    DM2 = exp(-mu_M*tauM)*exp(-b_H*tauM);
    DU2 = exp(-b_S*tauU)*exp(-phi_H*H*tauU);

#Model
    du[1] = lambda*P - mu_W*W - Phi(Temp(t))*S*W;                                                                                #W)
    du[2] = Phi(Temp(t))*S*W - U*(b_S + phi_H*H) - mu_U*U - Phi(Temp(t-tauU))*S*W_past*DU1*DU2*(g_U(Temp(t))/g_U(Temp(t-tauU))); #U)
    du[3] = Phi(Temp(t-tauU))*S*W_past*DU1*DU2*(g_U(Temp(t))/g_U(Temp(t-tauU))) - I*(mu_I+b_S) - phi_H*H*I;                      #I)
    du[4] = phi_H*H*I - M*(mu_M+b_H) - Alpha*P*M/H - phi_H*H*I_past*DM1*DM2;                                                     #M)
    du[5] = phi_H*H*I_past*DM1*DM2 - P*(mu_P+b_H) - Alpha*P;                                                                     #P)
    du[6] = DU1*(mu_U*U_past/S*(g_U(Temp(t))/g_U(Temp(t-tauU))) - mu_U*U/S);                                                     #DU1)
    du[7] = DM1*(Alpha*P_past/H - Alpha*P/H);                                                                                    #DM1)
    du[8] = 1-(g_U(Temp(t))/g_U(Temp(t-tauU)));                                                                                  #tauU) Change in delay duration
  end

  #Set up simulation
  function DySimulate(tstart,tfin)
    birthday = tstart;
    gint(u) = quadgk(gt, birthday-u, birthday)[1] #integrate over gt from days[jj]-u to days[jj]
    initTauU = [fzero(u -> gint(u) - p, [0.001, 5000]) for p in [1.0]][1]; #Initial lag duration

    #Initial survival probabilities
    mUint = exp(-mu_U*initTauU);
    bSint = exp(-b_S*initTauU);
    pHint = exp(-phi_H*H*initTauU);
    DUinit = mUint*bSint*pHint;
    DUhist = DUinit;
    DM2 = exp(-mu_M*tauM)*exp(-b_H*tauM);

    #Initial values
    u0 = [0.0,0.0,0.0,0.0,1.0,DUinit,DM2,initTauU];
    h(p,t) = [0.0,0.0,0.0,0.0,0.0,DUhist,DM2,initTauU]; #Empty system when t<0
    tspan = (birthday,tfin);

    #Declare dependent lags. What should they be? It works well when d_lags is not declared.
    d_lags = ((u,p,t) -> u[8],)
    prob = DDEProblem(DynDDE,u0,h,tspan; constant_lags = tauM, dependent_lags = d_lags)
    solDynDDE = solve(prob,alg; isoutofdomain=(u,p,t) -> any(x -> x < 0.0, u),abstol = 1e-6,reltol = 1e-4)
  end

    alg = MethodOfSteps(Rosenbrock23());
    out = DySimulate(1.0,365.0*10)
    plot(out, vars = (0,[1,2,3,4,5]))
    plot(out, vars = (0,[8])) #Visualize change in delay period over time (dtauU/dt)

Hej!

Your intuition is correct, dependent delays should be specified as tuple of functions of the form (u, p, t) -> lag. So in your case

prob = DDEProblem(DynDDE, u0, h, tspan;
                  constant_lags = [tauM],
                  dependent_lags = ((u, p, t) -> u[8],))

seems correct.

If you specify some function(s) lag(u, p, t) of dependent delays, DelayDiffEq tries to discover discontinuities induced by the dependent delays (which are not known at simulation start) in every integration step by checking the signs of the function f(t) = T - lag(u(t), p, t) - t at multiple interpolated points for all previous discontinuities T. If the algorithm notices that it would step over a discontinuity, the integration step is rejected. By doing so, it tries to ensure that the function u(t) is smooth enough in every integration step, which allows to bound the numerical integration error.

The tracking of these dependent discontinuities works just like a regular continuous callback (actually, it’s implemented as a special continuous callback) and can be adjusted by, e.g., increasing (or decreasing) the number of interpolation points.

1 Like

In general, for many algorithms the numerical integration error should decrease by specifying dependent delays and tracking the induced discontinuities. However, of course, tracking of the discontinuities requires additional operations and memory, so if you do not specify dependent delays the computation is expected to be faster (but also less accurate). In such a case residual control methods such as RK4 are usually recommended but also error prone (you can find some more explanations in the documentation.

I rewrote your example in the following way (I hope I did not introduce any mathematical errors :grimacing:). I tried to remove some redundant calculations and reduce the number of evaluations of the history function. Additionally, I changed the simulation function in a way that allows me to benchmark multiple different algorithms and solver options more easily. In particular, I tried to measure only the time it takes to solve the DDE, without all steps needed for setting it up.

using QuadGK
using Roots
using DelayDiffEq
using Plots

# Parameters
const Temp0K = 15 + 273.15
const TempKstop = 273.15
const k = 8.62 * 10^(-5)
const H = 0.075
const S = 4000000.0
const b_H = 1/5292.5
const b_S = 0.0037
const lambda = 488.0
const mu_W = 0.0186
const Phi_a = 0.024
const Phi_E = 0.46
const mu_P = 1.5e-7
const mu_M = 0.0056
const tauM = 166.82
const mu_U = 0.0001
const mu_I = mu_U
const tauU_a = 25.1
const tauU_E = 0.87
const tauU_EL = 5.0 * tauU_E
const tauU_TL = 282.0218
const Alpha = 1.4752e-6
const phi_H = 2.7728e-7
const c_K = 263.0
const d_K = 22.6
const t0 = 105.6
const DM2 = exp(-mu_M * tauM) * exp(-b_H * tauM)

# A couple functions that describe how parameters are sensitive to temperature

# Temperature as a function of time
temp(t) = c_K + d_K * sin((t - t0) * 2 * pi / 365)


# Development rate, as a function of Temperature
function gU(T)
    tmp = 1 / max(TempKstop, T)

    1 / (tauU_a * exp((tauU_E/k) * (tmp - (1 / Temp0K)))) *
        (1 + exp((tauU_EL / k) * (tmp - (1 / tauU_TL))))
end

# Infection rate, as a function of temperature
function phi(T)
    Phi_a * exp((-Phi_E / k) * ((1 / max(TempKstop, T)) - (1 / Temp0K)))
end

# Delay differential equation
function f(du, u, h, p, t)
    # State variables
    # W, U, I, M, and P track different stages of a parasite life-cycle
    # DU1 gives the proportion of U that survive the variable developmental delay period
    # DM1 give proportion of M that survive constant delay period
    # tauU tracks the length of the developmental delay, which changes continuously with time
    W = u[1]; U = u[2]; I = u[3]; M = u[4]; P = u[5]; DU1 = u[6]; DM1 = u[7]; tauU = u[8]

    # History calculations
    ucache = h(p, t - tauU)
    W_past = ucache[1]
    U_past = ucache[2]
    ucache = h(p, t - tauM)
    I_past = ucache[3]
    P_past = ucache[5]

    # "Constant" survival function
    DU2 = exp(-b_S * tauU) * exp(-phi_H * H * tauU)

    # Compute temperature
    T = temp(t)
    T_past = temp(t - tauU)

    # Preliminary calculations
    gU_ratio = gU(T) / gU(T_past)
    tmp1 = phi(T) * S * W
    tmp2 = phi(T_past) * S * W_past * DU1 * DU2 * gU_ratio
    tmp3 = phi_H * H * I
    tmp4 = phi_H * H * I_past * DM1 * DM2

    # Model
    du[1] = lambda * P - mu_W * W - tmp1
    du[2] = tmp1 - U * (mu_U + b_S + phi_H * H) - tmp2
    du[3] = tmp2 - I * (mu_I + b_S) - tmp3
    du[4] = tmp3 - M * (mu_M + b_H + Alpha * P / H) - tmp4
    du[5] = tmp4 - P * (mu_P + b_H + Alpha)
    du[6] = DU1 * mu_U * (U_past * gU_ratio - U) / S
    du[7] = DM1 * Alpha * (P_past - P) / H
    du[8] = 1 - gU_ratio

    nothing
end

# Set up simulation
function simulate(alg, tstart, tfin; with_dependent_lags = true, kwargs...)
    # Integrate over t -> gU(temp(t)) from tstart - u to tstart
    gint(u) = quadgk(gU ∘ temp, tstart - u, tstart)[1]

    # Initial lag
    initTauU = fzero(u -> gint(u) - 1.0, [0.001, 5000])

    # Initial survival probabilities
    mUint = exp(-mu_U * initTauU)
    bSint = exp(-b_S * initTauU)
    pHint = exp(-phi_H * H * initTauU)
    DUinit = mUint * bSint * pHint
    DUhist = DUinit
    DM2 = exp(-mu_M * tauM) * exp(-b_H * tauM)

    # Initial values
    u0 = [0.0, 0.0, 0.0, 0.0, 1.0, DUinit, DM2, initTauU]

    # Empty system when t < 0
    h(p, t) = [0.0, 0.0, 0.0, 0.0, 0.0, DUhist, DM2, initTauU]
    tspan = (tstart, tfin)

    # Declare dependent lags
    if with_dependent_lags
        prob = DDEProblem(f, u0, h, tspan;
                          constant_lags = [tauM],
                          dependent_lags = ((u, p, t) -> u[8],))
    else
        prob = DDEProblem(f, u0, h, tspan; constant_lags = [tauM])
    end

    # run once to compile
    solve(prob, alg; kwargs...)

    # time solution
    @time solve(prob, alg; kwargs...)
end

As a quick check, if I run

out = simulate(MethodOfSteps(Rosenbrock23()), 1.0, 3650.0; with_dependent_lags = true)
plot(out, vars = (0, [1, 2, 3, 4, 5]))

I get the following plot:

rosenbrock23_with

Does that look reasonable?

Now some timing results.

I get

julia> simulate(MethodOfSteps(Rosenbrock23()), 1.0, 3650.0; with_dependent_lags = true);
  0.321354 seconds (4.93 M allocations: 99.933 MiB, 6.10% gc time)

julia> simulate(MethodOfSteps(Rosenbrock23()), 1.0, 3650.0; with_dependent_lags = false);
  0.247855 seconds (2.86 M allocations: 67.018 MiB, 6.19% gc time)

As expected, in the first case with dependent delays the solution steps to more time points than in the second case (4255 vs 4150) but the computation is slower and allocates more memory. I get a very similar trajectory (but can not include it here since I just signed up for discourse).

Using Rodas4 can speed up the computations of this example since it only steps to 775 and 626 time points, respectively:

julia> simulate(MethodOfSteps(Rodas4()), 1.0, 3650.0; with_dependent_lags = true);
  0.213137 seconds (3.50 M allocations: 67.071 MiB, 7.82% gc time)

julia> simulate(MethodOfSteps(Rodas4()), 1.0, 3650.0; with_dependent_lags = false);
  0.159873 seconds (1.76 M allocations: 38.415 MiB, 6.33% gc time)

Since I just signed up for discourse, I’m not allowed to add more plots to this post, but in both cases the trajectories look very similar. However, the RMSD of the time delay component evaluated at time points 1,2,…,3650 in the simulations with Rosenbrock23 and Rodas4 (approx. 0.27 and 0.36 with and without dependent delays, respectively) is larger than the RMSD of the simulations with Rosenbrock23 (approx. 0.05) and Rodas4(approx. 0.07). Additionally, the simulation with Rodas4 without dependent delays contains small negative values of around -1.2e-21 whereas all other solutions do not contain any negative numbers.

If I specify the isoutofdomain keyword argument, I get very similar timings and exactly the same solutions:

julia> simulate(MethodOfSteps(Rosenbrock23()), 1.0, 3650.0;
                with_dependent_lags = true,
                isoutofdomain = (u, p, t) -> any(x -> x < 0.0, u));
  0.304523 seconds (4.93 M allocations: 99.933 MiB, 8.06% gc time)

julia> simulate(MethodOfSteps(Rosenbrock23()), 1.0, 3650.0;
                with_dependent_lags = false,
                isoutofdomain = (u, p, t) -> any(x -> x < 0.0, u));
  0.224484 seconds (2.86 M allocations: 67.018 MiB, 8.81% gc time)

Similar results as above I obtain for Rodas4:

julia> simulate(MethodOfSteps(Rodas4()), 1.0, 3650.0;
                with_dependent_lags = true,
                isoutofdomain = (u, p, t) -> any(x -> x < 0.0, u));
  0.208803 seconds (3.50 M allocations: 67.154 MiB, 8.36% gc time)

julia> simulate(MethodOfSteps(Rodas4()), 1.0, 3650.0;
                with_dependent_lags = false,
                isoutofdomain = (u, p, t) -> any(x -> x < 0.0, u));
  0.153502 seconds (1.79 M allocations: 38.982 MiB, 2.97% gc time)

However, in this case the solutions are different, with a larger discrepancy between the solutions without dependent delay (RMSD of the time delay component of approx. 0.07 and 0.12 if the dependent delays are specified and if they are not, respectively). As expected, if isoutofdomain is specified, all solutions contain only non-negative numbers.

Are these timings consistent with what you’ve observed?

1 Like

This is a comment to your reply below. I had to edit this comment since I’m not allowed to add another comment as a new user.

After fixing my typo I can reproduce the instability you observed. I’ll investigate it more closely, because I still think there could be an issue with the tracking of discontinuities - the discrepancy between both settings (even in the incorrect example above) seems quite large.

Maybe we can move the discussion to https://github.com/JuliaDiffEq/DelayDiffEq.jl/issues/120? Then it’s easier for me to reply :slight_smile:

1 Like

Wow, this is amazing! Thank you very much for spending so much time on this and getting back to me so quickly. That’s great that my intuition about setting up the dependent delays was correct, and now that you explained it, it does make sense that declaring the delays would lengthen simulation time. The main thing I wanted was confirmation that I was doing that right, but your performance tips are fantastic too. I will spend a lot more time going through all of this soon, but I did notice that you had a misplaced parenthesis in the development function, which shortened the delays quite a bit and reduced simulation time considerably. It should be

  function gU(T)
      tmp = 1 / max(TempKstop, T)

      1 / (tauU_a * exp((tauU_E/k) * (tmp - (1 / Temp0K))) *
        (1 + exp((tauU_EL / k) * (tmp - (1 / tauU_TL)))))
  end

instead of

  function gU(T)
      tmp = 1 / max(TempKstop, T)

      1 / (tauU_a * exp((tauU_E/k) * (tmp - (1 / Temp0K)))) *
        (1 + exp((tauU_EL / k) * (tmp - (1 / tauU_TL))))
  end

Correcting this slows things down a lot, and actually makes Rodas4() abort due to instability.
But perhaps I just need to increase the maxiters and let it run longer with Rosenbrock23()?

When I have dependent_lags = false, the simulation is way faster than it used to be, and gives me this:
MWE

but when I have dependent_lags = true I get an “interruption larger maxiters needed” warning and this:
MWEwlags

That sounds good, we can migrate over there. Thanks for your help with this, I really appreciate it. The numerics behind the scenes are a bit over my head, unfortunately.

So, it seems I’m able to answer here as well now :slight_smile: I just created a PR (https://github.com/JuliaDiffEq/DelayDiffEq.jl/pull/121) with some major changes that seem to fix the issues you experienced (at least partially - it’s still a bit slower and consumes more memory if you specify dependent delays). If you want to try it out, at the moment you have to use both the branch “callback” of DelayDiffEq and the branch “loopfooter” of OrdinaryDiffEq - there’s a small change necessary in OrdinaryDiffEq as well: https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/pull/822. I hope we can release new versions of DelayDiffEq and OrdinaryDiffEq with these changes in the next days.

2 Likes

This is great! I look forward to trying it out. I certainly wasn’t expecting my question to result in changes to the package! Thanks a lot for looking into this problem.