# Memory allocation when solving ODE with Vern7 and ContinuousCallback

**URL:** <https://discourse.julialang.org/t/memory-allocation-when-solving-ode-with-vern7-and-continuouscallback/101087>\
**Category:** Performance\
**Tags:** performance, ode, diffeqcallbacks\
**Created:** [July 2, 2023, 5:16pm UTC](https://discourse.julialang.org/t/memory-allocation-when-solving-ode-with-vern7-and-continuouscallback/101087 "2023-07-02T17:16:11Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [July 2, 2023, 5:16pm UTC](https://discourse.julialang.org/t/memory-allocation-when-solving-ode-with-vern7-and-continuouscallback/101087/1 "2023-07-02T17:16:11Z")

</div>

Hello,

I went through a memory allocation problem in a ODE system, when using Vern7() algorithm and ContinuousCallback. The minimal working example is the following:

```julia
using LinearAlgebra
using SparseArrays
using OrdinaryDiffEq
using DiffEqCallbacks
using BenchmarkTools

A = Diagonal(Array{ComplexF64}(1:20))
A += spdiagm(1=>sqrt.(0.01 .* Array(1:19)))
A -= 0.01im * Diagonal(Array{ComplexF64}(1:20))
u0 = rand(ComplexF64, 20)
prm = (U = A, rand_n = Ref(rand()))

my_condition_discrete(u, t, integrator) = real(dot(u, u)) < integrator.p.rand_n[]
my_condition_continuous(u, t, integrator) = integrator.p.random_n[] - real(dot(u, u))
my_affect!(integrator) = u_modified!(integrator, false)
cb1 = DiscreteCallback(my_condition_discrete, my_affect!, save_positions=(false, false))
cb2 = ContinuousCallback(my_condition_continuous, my_affect!, nothing, save_positions=(false, false))

prob1 = ODEProblem((du,u,p,t)->mul!(du, p.U, u), u0, (0.0, 100.0), prm, callback=cb1)
prob2 = ODEProblem((du,u,p,t)->mul!(du, p.U, u), u0, (0.0, 100.0), prm, callback=cb2)

@benchmark sol1 = solve($prob1, $Vern7(), saveat=$[100.0])

BenchmarkTools.Trial: 9021 samples with 1 evaluation.
 Range (min … max): 383.100 μs … 5.543 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 475.300 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 546.050 μs ± 198.079 μs ┊ GC (mean ± σ): 0.07% ± 0.91%

  █▃                                                             
  ███▇▇▆▆▅▅▄▅▄▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  383 μs Histogram: frequency by time 1.2 ms <

 Memory estimate: 10.34 KiB, allocs estimate: 50.

@benchmark sol2 = solve($prob2, $Vern7(), saveat=$[100.0])

BenchmarkTools.Trial: 2566 samples with 1 evaluation.
 Range (min … max): 1.274 ms … 12.947 ms ┊ GC (min … max): 0.00% … 79.46%
 Time (median): 1.709 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 1.931 ms ± 783.272 μs ┊ GC (mean ± σ): 2.63% ± 7.29%

  ▇█▆█▇▄▅▃▃▂▁                                                  
  ████████████▇█▆▅▅▅▄▄▄▄▅▄▄▅▃▄▃▃▃▃▃▃▃▂▂▃▂▂▂▂▂▂▂▂▂▂▁▂▁▂▂▂▂▂▁▂▂ ▄
  1.27 ms Histogram: frequency by time 4.71 ms <

 Memory estimate: 784.83 KiB, allocs estimate: 2032.

```

While it doesn’t happen for example with the `Tsit5()` algorithm:

```julia
@benchmark sol1 = solve($prob1, $Tsit5(), saveat=$[100.0])

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 309.600 μs … 5.215 ms ┊ GC (min … max): 0.00% … 89.28%
 Time (median): 392.600 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 436.892 μs ± 136.050 μs ┊ GC (mean ± σ): 0.11% ± 0.89%

   ▇█                                                            
  ▅███▇▆█▇▅▆▄▅▄▄▅▄▃▄▃▄▄▃▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  310 μs Histogram: frequency by time 859 μs <

 Memory estimate: 9.44 KiB, allocs estimate: 48.

@benchmark sol2 = solve($prob2, $Tsit5(), saveat=$[100.0])

BenchmarkTools.Trial: 4593 samples with 1 evaluation.
 Range (min … max): 744.500 μs … 3.187 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 968.600 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 1.079 ms ± 306.638 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

    █▃▃                                                          
  ▃▆███▇▇▇▇▆▄▄▄▅▅▄▄▄▄▃▄▄▃▃▃▃▃▃▃▂▂▂▂▃▂▂▂▂▂▂▂▂▁▂▂▂▂▂▁▂▁▁▁▁▁▁▁▁▁▁▁ ▂
  744 μs Histogram: frequency by time 2.07 ms <

 Memory estimate: 9.69 KiB, allocs estimate: 48.

```

---

<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:** [July 2, 2023, 10:28pm UTC](https://discourse.julialang.org/t/memory-allocation-when-solving-ode-with-vern7-and-continuouscallback/101087/2 "2023-07-02T22:28:09Z")

</div>

Yes it’s something odd with type inference. This is a good example of it and I’ll see if someone from the compiler team can help figure out what it is. I’ve seen this for years but haven’t figured out how to handle it yet, but this is a cleaner example than the others I had.

---

<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:** [August 29, 2023, 5:27am UTC](https://discourse.julialang.org/t/memory-allocation-when-solving-ode-with-vern7-and-continuouscallback/101087/3 "2023-08-29T05:27:39Z")

</div>

```julia
using LinearAlgebra
using SparseArrays
using OrdinaryDiffEq
using DiffEqCallbacks
using BenchmarkTools

A = Diagonal(Array{ComplexF64}(1:20))
A += spdiagm(1=>sqrt.(0.01 .* Array(1:19)))
A -= 0.01im * Diagonal(Array{ComplexF64}(1:20))
u0 = rand(ComplexF64, 20)
prm = (U = A, rand_n = Ref(rand()))

my_condition_discrete(u, t, integrator) = real(dot(u, u)) < integrator.p.rand_n[]
my_condition_continuous(u, t, integrator) = integrator.p.rand_n[] - real(dot(u, u))
my_affect!(integrator) = u_modified!(integrator, false)
cb1 = DiscreteCallback(my_condition_discrete, my_affect!, save_positions=(false, false))
cb2 = ContinuousCallback(my_condition_continuous, my_affect!, nothing, save_positions=(false, false))

prob1 = ODEProblem((du,u,p,t)->mul!(du, p.U, u), u0, (0.0, 100.0), prm, callback=cb1)
prob2 = ODEProblem((du,u,p,t)->mul!(du, p.U, u), u0, (0.0, 100.0), prm, callback=cb2)

@benchmark sol1 = solve($prob1, $Vern7(lazy=false), saveat=$[100.0])

@benchmark sol2 = solve($prob2, $Vern7(lazy=false), saveat=$[100.0])

```

The actual reason for this is because `Vern7` is uses a lazy interpolation which only extends to every `f` call if the interpolation is used, but if you use a ContinuousCallback it will always extend, in which case simply turning off the lazy interpolation would be more efficient.