Delay differential equation (DDE) solvers in Julia for the SciML scientific machine learning ecosystem. Covers neutral and retarded delay differential equations, and differential-algebraic equations.
Describe the bug π
Unexpected large memory usage proportional to tspan[2]-tspan[1] when solving a Delay Differential Equation with an oscillatory driving function. This might not be a bug, and may actually be a feature request to specify a maximum time lag I promise to never look beyond to limit the memory used by the history function.
Expected behavior
I expect memory usage to scale close to the amount of data to be saved. In this case for each time point it saves two state variables and time, at 8 bytes each, I expect ~24 bytes per sample. Instead sol.interp.f stores ~14000 bytes/sample. The amount required by sol.interp.f scales proportional to tspan[2]-tspan[1]. Clearly some amount of memory is required to enable the history function, but this seems excessive. In this MRE the single constant lag is less than 1 part in 10,000 of tspan[2]-tspan[1], so I expect that the memory usage required for the history function should not need to grow when tspan[2]-tspan[1] is increased.
The MRE includes a delay diff eq setup, and some code to calculate the bytes/sample used by various fields of the solution object sol. I want to solve for much longer values of tspan[2]-tspan[1] and the memory usage is the primary limiting factor. The dominant contribution by at least a factor of 100 is the field sol.interp.f. I believe this field is part of the implementation of the history function for the DDE.
Minimal Reproducible Example π
using DifferentialEquations
function qw_dde(du,u,h,params,t)
a,b,c,d,rt,Ο = params
hId,hIu = h(params,t-Ο)
hdId,hdIu = h(params,t-Ο,Val{1})
du[1] = ( 4.28e-7*(-2.51e10)*sin(2.51e10*t) + c*(hIu - u[1]) + b*hdIu ) * a
du[2] = hdId + (rt-1.0)*(u[2] + hId)/(d)
end
const tspan = (0.0, 1e-6)
const Ο = 5.6e-11
const params = (0.66, -0.5, 4.4e11, 4.0e-12, 7.8e-6, Ο)
h(params,t) = zeros(2);
h(params,t,::Type{Val{1}}) = zeros(2);
u0 = [0.0,0.0];
alg = MethodOfSteps(Tsit5());
saveat = LinRange(tspan[1], tspan[2], 5000)
prob = DDEProblem(qw_dde,u0,h,tspan,params; saveat=saveat, constant_lags=[Ο]);
# Run simulation:
println("solve start")
@time sol = solve(prob,alg,reltol=1e-3,abstol=1e-9,maxiters=1e9);
println("solve done")
# investigate space usage
function pretty_print2(d::Dict, pre=1)
# https://stackoverflow.com/questions/48195775/how-to-pretty-print-nested-dicts-in-julia
todo = Vector{Tuple}()
for (k,v) in d
if typeof(v) <: Dict
push!(todo, (k,v))
else
println(join(fill(" ", pre)) * "$(repr(k)) => $(repr(v))")
end
end
for (k,d) in todo
s = "$(repr(k)) => "
println(join(fill(" ", pre)) * s)
pretty_print2(d, pre+1+length(s))
end
nothing
end
# define a version of summarysize that doesn't choke on Vector{Vector{Vector{Float64}}}
# and returns a dict so we can see size by property
vec_safe_summary_size(x::Vector) = sum(vec_safe_summary_size(z) for z in x)
vec_safe_summary_size(x) = Base.summarysize(x)
# vec_safe_summary_size(sol::DelayDiffEq.HistoryODEIntegrator) = Dict(prop => vec_safe_summary_size(getproperty(sol, prop)) for prop in propertynames(sol))
# vec_safe_summary_size(sol::DelayDiffEq.HistoryFunction) = Dict(prop => vec_safe_summary_size(getproperty(sol, prop)) for prop in propertynames(sol))
# vec_safe_summary_size(sol::DelayDiffEq.ODEFunctionWrapper) = Dict(prop => vec_safe_summary_size(getproperty(sol, prop)) for prop in propertynames(sol))
vec_safe_summary_size(sol::OrdinaryDiffEq.InterpolationData) = Dict(prop => vec_safe_summary_size(getproperty(sol, prop)) for prop in propertynames(sol))
vec_safe_summary_size(sol::ODESolution) = Dict(prop => vec_safe_summary_size(getproperty(sol, prop)) for prop in propertynames(sol))
v = vec_safe_summary_size(sol)
calc_per_point(x, n) = x/n
calc_per_point(d::Dict, n) = Dict(key=>calc_per_point(val, n) for (key, val) in pairs(d))
println("bytes/point by property of sol")
pretty_print2(calc_per_point(v, length(sol)))
@show length(sol);
@show tspan;
Error & Stacktrace β οΈ
Environment (please complete the following information):
Output of using Pkg; Pkg.status()
julia> using Pkg; Pkg.status()
Status `C:\Users\oneilg\AppData\Local\Temp\jl_tkEYFY\Project.toml`
[0c46a032] DifferentialEquations v7.13.0
Output of using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)
Describe the bug π Unexpected large memory usage proportional to
tspan[2]-tspan[1]when solving a Delay Differential Equation with an oscillatory driving function. This might not be a bug, and may actually be a feature request to specify a maximum time lag I promise to never look beyond to limit the memory used by the history function.Expected behavior I expect memory usage to scale close to the amount of data to be saved. In this case for each time point it saves two state variables and time, at 8 bytes each, I expect ~24 bytes per sample. Instead
sol.interp.fstores ~14000 bytes/sample. The amount required bysol.interp.fscales proportional totspan[2]-tspan[1]. Clearly some amount of memory is required to enable the history function, but this seems excessive. In this MRE the single constant lag is less than 1 part in 10,000 oftspan[2]-tspan[1], so I expect that the memory usage required for the history function should not need to grow whentspan[2]-tspan[1]is increased.The MRE includes a delay diff eq setup, and some code to calculate the bytes/sample used by various fields of the solution object
sol. I want to solve for much longer values oftspan[2]-tspan[1]and the memory usage is the primary limiting factor. The dominant contribution by at least a factor of 100 is the fieldsol.interp.f. I believe this field is part of the implementation of the history function for the DDE.Minimal Reproducible Example π
Error & Stacktrace β οΈ
Environment (please complete the following information):
using Pkg; Pkg.status()using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)versioninfo()Additional context
Add any other context about the problem here.