# Delay Differential Equations Networks with different delay per edge

**URL:** <https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787>\
**Category:** Modelling & Simulations\
**Tags:** diffeq, performance\
**Created:** [October 14, 2021, 10:41pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787 "2021-10-14T22:41:49Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ger](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ger/32/23817_2.png) [@Ger](https://discourse.julialang.org/u/Ger)\
**Post date:** [October 14, 2021, 10:41pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/1 "2021-10-14T22:41:49Z")

</div>

Hi everyone! While Julia-evangelizing (as usual) in my research group, in order to encourage them to move to Julia I offered to convert some of their Python codes for solving a network of Kuramoto oscillators with delay (with a different delays for each connection). To my surprise, the Python version (which uses [JiTCDDE](https://jitcdde.readthedocs.io/en/stable/)) is much more performant (and scales much better with the system size) than what I could achieve in Julia so far. I must be committing a serious mistake here but I’m not seeing it. I would be really grateful if anyone could help me with this.

To have a sense of the performances:

#### Julia:

 ![kuramoto_julia_times_const_vec_lag_minorgrid](https://global.discourse-cdn.com/julialang/original/3X/a/2/a20a7872773a1ef9815b70f7ae4b7fb3a5e8cf2d.jpeg)  
(median execution times measured with BenchmarkTools.jl)

#### Python:

 ![times_python](https://global.discourse-cdn.com/julialang/original/3X/3/7/3710cc2d37bdc56e260ad7035f1f2097b616ed85.jpeg)  
(this is discarding the compilation time)

At first I tried to implement it with NetworkDynamics.jl which seemed appropriate for this case but I didn’t find the way to set different lags for each edge (does anyone know if this is possible?) so I moved to plain DifferentialEquations.jl :

### Julia code

```julia
using DelayDiffEq
using Distributions
using Random

function f!(dθ, θ, h, p, t)
    ω, A = p
    n = length(θ)
    lags = reshape(lag, n,n)
    @inbounds for j in 1:n
        coupling = 0.0
        @inbounds for i in 1:n
            coupling += A[i,j]*sin(h(p, t-lags[i,j]; idxs=i) - θ[j])
        end
        dθ[j] = ω[j] + coupling
    end
    nothing
end

n = 10
Random.seed!(1)
ω = rand(n)
A = rand(n,n)
const lag = rand(n*n)
θ₀ = rand(Uniform(0, 2π), n)
p = (ω, A)
h(p, t; idxs=nothing) = typeof(idxs) <: Number ? 1.0 : rand(Uniform(0, 2π), n)

prob = DDEProblem(f!, θ₀, h, (0.0, 1.0), p, constant_lags=lag)
solve(prob, MethodOfSteps(BS3()), saveat=0.01, rtol=0, atol=1e-5, force_dtmin=true)

```

(I made the diffeq function inplace because the goal is to simulate 180 oscillators)

```julia
      Status `~/Documents/julia_kuramoto/Project.toml`
  [6e4b80f9] BenchmarkTools v1.2.0
  [bcd4f6db] DelayDiffEq v5.31.1
  [31c24e10] Distributions v0.25.19
  [9a3f8284] Random

```

### Python code:

```python
import os
from tqdm import tqdm
import numpy as np
from numpy import pi, random, max

os.environ["CC"] = "clang"
# https://jitcdde.readthedocs.io/en/stable/
from jitcdde import jitcdde, y, t
from symengine import sin
import symengine
import sympy

n = 10
ω = random.rand(n)
A = random.rand(n, n)
τ = random.rand(n, n)

def kuramotos():
    for i in range(n):
        yield ω[i] + sum(
            A[j, i]*sin(y(j, t-τ[i, j])-y(i))
            for j in range(n)
        )

DDE = jitcdde(kuramotos, n=n, verbose=True)
DDE.compile_C(simplify=False, do_cse=False, chunk_size=150)
DDE.set_integration_parameters(rtol=0, atol=1e-5)
DDE.constant_past(random.uniform(0, 2*np.pi, n), time=0.0)
DDE.integrate_blindly(max(τ), 0.1)
output = []
for time in tqdm(DDE.t + arange(0, 1, 0.01)):
    output.append([*DDE.integrate(time)])

```

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [October 15, 2021, 5:37am UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/2 "2021-10-15T05:37:29Z")

</div>

It might help if you put the fact that python seems to outperform Julia here right in the title 🙂

Also I know he’s super busy and would probably get to this one eventually anyway, I’ll still tag @chrisrackaukas as he’s always keen to find potential performance issues in the DiffEq ecosystem.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [October 15, 2021, 8:20am UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/3 "2021-10-15T08:20:01Z")

</div>

I am not sure your code is type stable because of `coupling`.

1. Are you comparing the same methods? (BS3 vs …?)
2. force dtmin= true is suspicious. Do you get a warning from DE?

---

<div class="post-metadata">

**Author:** ![jipolanco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jipolanco/32/12129_2.png) [@jipolanco](https://discourse.julialang.org/u/jipolanco)\
**Post date:** [October 15, 2021, 9:05am UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/4 "2021-10-15T09:05:54Z")

</div>

A possible problem is that in your definition of `h`, the parameter `n` is a global non-`const` variable, which surely leads to type instability:

```julia
n = 10
# [...]
h(p, t; idxs=nothing) = typeof(idxs) <: Number ? 1.0 : rand(Uniform(0, 2π), n)

```

You may try setting it as `const`, or better, including it as a parameter in the `p` argument of `h`.

---

<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:** [October 15, 2021, 1:38pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/5 "2021-10-15T13:38:07Z")

</div>

> [@rveltz](#):
>
> force dtmin= true is suspicious. Do you get a warning from DE?

Yeah, if you’re doing this then almost certainly something is wrong…

---

<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:** [October 15, 2021, 2:00pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/6 "2021-10-15T14:00:29Z")

</div>

`solve(prob, MethodOfSteps(BS3()), saveat=0.01, reltol=0, abstol=1e-5)`. `rtol` and `atol` are not options, and you should almost never `force_dtmin`. I should really just remove that argument. If you’re hitting it, then your definition is probably wrong: the warning even says to think of dtmin as divergence for a reason described in [PSA: How to help yourself debug differential equation solving issues](https://discourse.julialang.org/t/psa-how-to-help-yourself-debug-differential-equation-solving-issues/62489) .

And so the first thing is to notice is that these are not the same DDE being solved. In fact, the Julia one isn’t even a DDE: it’s not even a well-defined equation!

```julia
h(p, t; idxs=nothing) = typeof(idxs) <: Number ? 1.0 : rand(Uniform(0, 2π), n)

```

The past does not have a value: every single time the past is called the value changes. So what’s happening is that the starting iterations are unable to converge because there is no sense of convergence: there are no values!

But it’s even worse then that: you’re saying that "the u(t) for t\<0 is a random value at every `t`, but `u(t)[1] = 1.0` for any past value. So that’s even worse, it’s not even consistent 😱. And `u(t)[1] = 1.0` does not match what is in the Python code which uses a uniform random value for the past.

The fixed code is:

```julia
using DelayDiffEq
using Distributions
using Random

function f!(dθ, θ, h, p, t)
    ω, A = p
    n = length(θ)
    lags = reshape(lag, n,n)
    @inbounds for j in 1:n
        coupling = 0.0
        @inbounds for i in 1:n
            coupling += A[i,j]*sin(h(p, t-lags[i,j]; idxs=i) - θ[j])
        end
        dθ[j] = ω[j] + coupling
    end
    nothing
end

n = 10
Random.seed!(1)
ω = rand(n)
A = rand(n,n)
const lag = rand(n*n)
θ₀ = rand(Uniform(0, 2π), n)
p = (ω, A)
const past = rand(Uniform(0, 2π), n)
h(p, t; idxs=nothing) = typeof(idxs) <: Number ? past[idxs] : past

prob = DDEProblem(f!, θ₀, h, (0.0, 1.0), p, constant_lags=lag)

using Profile
@profile solve(prob, MethodOfSteps(BS3()), saveat=0.01, reltol=0, abstol=1e-5)
Juno.profiler()

```

After fixing that it’s still slow though. The issue seems to be that Julia’s inference is bugging in `f!`. This could be [https://github.com/JuliaLang/julia/issues/35800](https://github.com/JuliaLang/julia/issues/35800) or [https://github.com/SciML/OrdinaryDiffEq.jl/pull/1473](https://github.com/SciML/OrdinaryDiffEq.jl/pull/1473)

---

<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:** [October 15, 2021, 2:37pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/7 "2021-10-15T14:37:27Z")

</div>

Oh interesting, one thing that is happening is Julia’s function handling heuristics are getting in the way.

```julia
function f!(dθ, θ, h::H, p, t) where H
    ω, A = p
    n = length(θ)
    lags = reshape(lag, n,n)
    @inbounds for j in 1:n
        coupling = 0.0
        @inbounds for i in 1:n
            coupling += A[i,j]*sin(h(p, t-lags[i,j]; idxs=i) - θ[j])
        end
        dθ[j] = ω[j] + coupling
    end
    nothing
end

```

is 6 times faster. Before:

```julia
julia> @time solve(prob, MethodOfSteps(BS3()), saveat=0.01, reltol=0, abstol=1e-5)
 56.409539 seconds (1.96 G allocations: 32.454 GiB, 7.17% gc time)

```

After:

```julia
julia> @time solve(prob, MethodOfSteps(BS3()), saveat=0.01, reltol=0, abstol=1e-5)
 12.284248 seconds (5.79 M allocations: 455.702 MiB, 1.67% gc time)

```

[https://github.com/SciML/SciMLBase.jl/pull/110](https://github.com/SciML/SciMLBase.jl/pull/110) fixes that.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [October 15, 2021, 3:23pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/8 "2021-10-15T15:23:47Z")

</div>

Keep in mind JiTCDDE is itself generating and compiling C-code based on user-specified symbolic expressions, and then just putting a Python wrapper around that code. So it might be pretty efficient to begin with. The DiffEq solvers shouldn’t be lots slower, but for the same basic numerical method it isn’t clear that this is a case where one would expect them to be noticeably faster.

---

<div class="post-metadata">

**Author:** ![Ger](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ger/32/23817_2.png) [@Ger](https://discourse.julialang.org/u/Ger)\
**Post date:** [October 15, 2021, 3:52pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/9 "2021-10-15T15:52:26Z")

</div>

Thanks for your reply!

> [@rveltz](#):
>
> I am not sure your code is type stable because of `coupling` .

I noticed that the function was not type stable when not declaring `lag` as `const`, which I think was expected because the `lag` was only being declared in the `constant_lags` keyword argument but the function itself didn’t have information about it. This is also related with what Chris noticed about the Julia heuristics on `h` function.

> [@rveltz](#):
>
> - Are you comparing the same methods? (BS3 vs …?)

I think I’m comparing with the [same method](https://jitcdde.readthedocs.io/en/stable/index.html?highlight=Shampine#a-brief-mathematic-background).

> [@rveltz](#):
>
> 1. force dtmin= true is suspicious. Do you get a warning from DE?

I’m sorry, I left the `dtmin = true` from some previous trials with more realistic connectivities and delays which where raising some warnings but it is not necessary in this test case. As Chris said I should just remove it and probably in the other case I was getting the warning because the problem was not well defined.

---

<div class="post-metadata">

**Author:** ![Ger](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ger/32/23817_2.png) [@Ger](https://discourse.julialang.org/u/Ger)\
**Post date:** [October 15, 2021, 4:11pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/10 "2021-10-15T16:11:54Z")

</div>

Thanks a lot Chris!! And I’m sorry for such embarrassing mistakes 😅

> [@ChrisRackauckas](#):
>
> `rtol` and `atol` are not options

I know, I was biased by the Python notation that I had just looked at and made the mistake 😅

> [@ChrisRackauckas](#):
>
> you should almost never `force_dtmin` . I should really just remove that argument. If you’re hitting it, then your definition is probably wrong: the warning even says to think of dtmin as divergence for a reason described in [PSA: How to help yourself debug differential equation solving issues](https://discourse.julialang.org/t/psa-how-to-help-yourself-debug-differential-equation-solving-issues/62489) .

About `force_dtmin`, I left it from a previous trial with 180 oscillators and more realistic connectivities and delays, which was raising some warning. It was not necessary in this case. And furthermore, as you say, the warnings might have been due to the wrongly defined problem and I shouldn’t have use it in the first place. Thanks for the link!

> [@ChrisRackauckas](#):
>
> The past does not have a value: every single time the past is called the value changes. So what’s happening is that the starting iterations are unable to converge because there is no sense of convergence: there are no values!
> 
> But it’s even worse then that: you’re saying that "the u(t) for t\<0 is a random value at every `t` , but `u(t)[1] = 1.0` for any past value. So that’s even worse, it’s not even consistent 😱. And `u(t)[1] = 1.0` does not match what is in the Python code which uses a uniform random value for the past.

Thanks for noticing and explaining this! This was the first time I was using DDE and I clearly misunderstood how `h` function was being used.

---

<div class="post-metadata">

**Author:** ![Ger](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ger/32/23817_2.png) [@Ger](https://discourse.julialang.org/u/Ger)\
**Post date:** [October 15, 2021, 4:25pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/11 "2021-10-15T16:25:46Z")

</div>

Thanks a lot, Chris!! This fix indeed increases significantly the performance! However it still much less performant than the Python version and scales really badly with the size of the system. For example, for `n = 20` :

```julia
julia> @time solve(prob, MethodOfSteps(BS3()), saveat=0.01, reltol=0, abstol=1e-5)
119.585544 seconds (11.80 M allocations: 1.706 GiB, 3.26% gc time)

```

while in the Python version it takes ~ 0.007 seconds.

My real need is to run it for `n = 180` so I’ll need to improve the performance much more.

---

<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:** [October 15, 2021, 4:34pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/12 "2021-10-15T16:34:53Z")

</div>

Yes, something is happening with the interpolation and I won’t be able to look at that today

---

<div class="post-metadata">

**Author:** ![Ger](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ger/32/23817_2.png) [@Ger](https://discourse.julialang.org/u/Ger)\
**Post date:** [October 15, 2021, 4:43pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/13 "2021-10-15T16:43:28Z")

</div>

No problem, I truly appreciate all the help that you gave me so far!! 🙂

---

<div class="post-metadata">

**Author:** ![deep-introspection](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/deep-introspection/32/21522_2.png) [@deep-introspection](https://discourse.julialang.org/u/deep-introspection)\
**Post date:** [November 9, 2021, 6:48pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/14 "2021-11-09T18:48:31Z")

</div>

Running into similar problems at a large scale too (\>100 oscillators), and cannot really compute the solution before timeout on my HPC. Is there a way to help solve this issue?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [November 9, 2021, 7:34pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/15 "2021-11-09T19:34:05Z")

</div>

Have you confirmed that the same issue is occurring?

---

<div class="post-metadata">

**Author:** ![deep-introspection](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/deep-introspection/32/21522_2.png) [@deep-introspection](https://discourse.julialang.org/u/deep-introspection)\
**Post date:** [November 18, 2021, 1:25pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/16 "2021-11-18T13:25:53Z")

</div>

Yes, and since I’m actually interested in the same system that is being discussed here, I used the fixed version proposed by @chrisrackaukas (so it would be the same issue). I’d love to use DifferentialEquations.jl and all its goodies to solve this problem but the current prohibitive performance issue is making me fall back to Python and JITCODE.

---

<div class="post-metadata">

**Author:** ![lindnemi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lindnemi/32/23421_2.png) [@lindnemi](https://discourse.julialang.org/u/lindnemi)\
**Post date:** [November 22, 2021, 3:45pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/17 "2021-11-22T15:45:15Z")

</div>

Same problem here. I implemented a quite general interface for networks with heterogeneous delays in NetworkDynamics.jl but unless this issue gets fixed it’s pretty useless: [https://github.com/PIK-ICoNe/NetworkDynamics.jl/pull/104](https://github.com/PIK-ICoNe/NetworkDynamics.jl/pull/104)

---

<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:** [November 22, 2021, 11:05pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/18 "2021-11-22T23:05:09Z")

</div>

Alright, but it’s a Base inference issue so I don’t know what to say. Help fix Base and that would be a good starting point.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 26, 2021, 11:20pm UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/19 "2021-11-26T23:20:39Z")

</div>

Is there a bug filed against Base?

---

<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:** [November 27, 2021, 10:45am UTC](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787/20 "2021-11-27T10:45:36Z")

</div>

Yes, it’s already posted in this thread.

[Next page](https://discourse.julialang.org/t/delay-differential-equations-networks-with-different-delay-per-edge/69787.md?page=2)