Huge detriment in computational performance

[rvignolo]

Hi all!

First of all, I would like to thank you for your package, it's really nice and very helpful. Now let me explain the issue I am facing.

I am seeing a huge decrement in performance once I moved from version v0.9.2 to v0.9.3, and it still persists in the current version v0.9.5.

I prepared a minimal working example for the testing, which seems not so minimal but it is 😄:

using Parameters
using FiniteDifferences
using Dierckx: Spline1D
using DifferentialEquations

params = @with_kw((
  LY = Spline1D(
    [0.08, 0.1666667, 0.25, 0.5, 1.0, 2.0, 3.0, 5.0, 7.0, 10.0, 20.0, 30.0],
    [
      0.0242,
      0.0243,
      0.0243,
      0.0245,
      0.0236,
      0.0226,
      0.0223,
      0.0226,
      0.0237,
      0.0247,
      0.027,
      0.0289,
    ],
    k = 1,
    bc = "nearest",
  ),

  P1 = t -> 1 / (1 + LY(t) * t),
  P2 = t -> 1 / ((1 + LY(t))^t),
  P  = t -> ifelse(t < 1.0, P1(t), P2(t)),

  f = t -> FiniteDifferences.forward_fdm(5, 1)(s -> -log(P(s)), t),
  f′ = t -> FiniteDifferences.forward_fdm(5, 2)(s -> -log(P(s)), t),

  # why these definitions makes ϑ type unstable if they are type stable?
  # f′ = t -> FiniteDifferences.forward_fdm(5, 1)(f, t),
  # f′ = t -> FiniteDifferences.forward_fdm(5, 1)(s -> f(s), t),

  ϰ = 0.4363,
  σ = 0.1491,
  ϑ = t -> 1 / ϰ * f′(t) + f(t) + (σ / ϰ)^2 / 2 * (1 - exp(-2 * ϰ * t)),
))

function riccati(du, u, p, t)

  @unpack (LY ,P1, P2, P, f, f′, ϰ, σ, ϑ) = p

  cattle = u[1]
  snake = u[2]
  locust = f(t)
  lion = 1
  mouse = ϰ
  loris = ϑ(t)
  aardvark = σ
  deer = 1
  coati = 0
  snail = FiniteDifferences.forward_fdm(5, 1)(f, t)
  crocodile = 0

  du[1] = (((mouse * loris - mouse * locust) - snail) * snake) / lion - ((deer + coati * locust) / 2) * ((aardvark * snake) / lion) ^ 2
  du[2] = ((mouse + crocodile / lion) * snake + (coati / (2lion)) * (aardvark * snake) ^ 2) - lion

  return nothing
end

let
  p = params()

  tspan = (0.0, 1.0)

  Δt = (tspan[2] - tspan[1])

  dt = 1.0 / 252.0

  nodes = Int(ceil(Δt / dt) + 1)

  t = T = [tspan[1] + (i - 1) * dt for i = 1:nodes]

  prob = ODEProblem(riccati, zeros(2), (0.0, 1.0), p)

  prob_func = (prob, i, repeat) -> begin
    ODEProblem(prob.f, prob.u0, (T[i + 1], t[1]), prob.p)
  end

  odes = EnsembleProblem(prob, prob_func = prob_func)

  @time sol = DifferentialEquations.solve(
    odes,
    Tsit5(),
    trajectories = nodes - 1,
    saveat = -dt
  )

end

The important thing is that I have to evaluate both f and f′ many times in this calculation, which means that I have to call FiniteDifferences.forward_fdm many times as well. I also tested with different algorithms, such as central_fdm or backward_fdm and noted the same detriment in performance.

Let me show you my @time results for each version (multiple runs for each case, the first one considers the overhead regarding compilation):

v0.9.2:

15.069737 seconds (41.47 M allocations: 2.262 GiB, 6.35% gc time)
1.536485 seconds (10.84 M allocations: 768.160 MiB, 13.14% gc time)
1.476522 seconds (10.65 M allocations: 759.003 MiB, 11.87% gc time)
1.375865 seconds (10.65 M allocations: 759.003 MiB, 12.61% gc time)

v0.9.3:

34.916895 seconds (60.25 M allocations: 3.414 GiB, 4.11% gc time)
23.306021 seconds (28.84 M allocations: 1.860 GiB, 3.25% gc time)
22.728664 seconds (28.65 M allocations: 1.851 GiB, 3.22% gc time)
21.738619 seconds (28.65 M allocations: 1.851 GiB, 3.35% gc time)

v0.9.5 (current):

35.646730 seconds (60.24 M allocations: 3.413 GiB, 4.04% gc time)
26.331707 seconds (45.91 M allocations: 2.749 GiB, 3.83% gc time)
22.284941 seconds (28.65 M allocations: 1.851 GiB, 3.16% gc time)
21.381123 seconds (28.65 M allocations: 1.851 GiB, 3.05% gc time)

As you see, this reduces performance by a factor of ≈ 20. It is important to remark that all versions lead to the correct answer.

Could you please help me find a way around this issue?

Thank you very much!

Regards,

Ramiro.

[oxinabox]

It is either #65 or #66
That is what I get for not doing seperate releases between them
One can go and ] dev FiniteDIfferences and git checkout the appropriate commit.

East short term fix is to set your compat to `FiniteDifferences = "~0.9.2" while we work it out.

[rvignolo]

Hi @oxinabox, thanks for your quick reply!

I will revert to version 0.9.2 while you work it out. I will probably ask you about this issue in a while if that's okay.

Thanks for your help!

[oxinabox]

stats on my computer.
Looks like the cause was #65.
Am going to build a fix shortly

v0.9.5

44.084013 seconds (80.42 M allocations: 4.397 GiB, 3.80% gc time)
29.377420 seconds (28.23 M allocations: 1.830 GiB, 1.86% gc time)
27.364521 seconds (28.23 M allocations: 1.830 GiB, 1.76% gc time)
27.871987 seconds (28.23 M allocations: 1.830 GiB, 1.79% gc time)

v0.9.2

18.984412 seconds (61.60 M allocations: 3.236 GiB, 6.57% gc time)
1.072190 seconds (10.23 M allocations: 737.725 MiB, 10.14% gc time)
1.090849 seconds (10.23 M allocations: 737.725 MiB, 9.55% gc time)
1.072175 seconds (10.23 M allocations: 737.725 MiB, 8.74% gc time)

#65 only

58.179475 seconds (80.43 M allocations: 4.387 GiB, 3.53% gc time)
30.308075 seconds (28.23 M allocations: 1.830 GiB, 2.03% gc time)
28.976301 seconds (28.23 M allocations: 1.830 GiB, 2.10% gc time)
29.478819 seconds (28.23 M allocations: 1.830 GiB, 2.16% gc time)

[oxinabox]

I think this never hits me in general because I create one global constant _fdm and then reuse it.
Which even with the fix I am about to make is still kinda a decent idea.

[rvignolo]

I migrated my MWE based in your comment about the const _fdm to the following:

using Parameters
using FiniteDifferences
using Dierckx: Spline1D
using DifferentialEquations

params = @with_kw((
  LY = Spline1D(
    [0.08, 0.1666667, 0.25, 0.5, 1.0, 2.0, 3.0, 5.0, 7.0, 10.0, 20.0, 30.0],
    [
      0.0242,
      0.0243,
      0.0243,
      0.0245,
      0.0236,
      0.0226,
      0.0223,
      0.0226,
      0.0237,
      0.0247,
      0.027,
      0.0289,
    ],
    k = 1,
    bc = "nearest",
  ),

  P1 = t -> 1 / (1 + LY(t) * t),
  P2 = t -> 1 / ((1 + LY(t))^t),
  P  = t -> ifelse(t < 1.0, P1(t), P2(t)),

  _fdm = forward_fdm(5, 1),
  f = t -> _fdm(s -> -log(P(s)), t),
  f′ = t -> _fdm(s -> f(s), t),


  ϰ = 0.4363,
  σ = 0.1491,
  ϑ = t -> 1 / ϰ * f′(t) + f(t) + (σ / ϰ)^2 / 2 * (1 - exp(-2 * ϰ * t)),
))


function riccati(du, u, p, t)
  @unpack (LY ,P1, P2, P, f, f′, ϰ, σ, ϑ, _fdm) = p

  cattle = u[1]
  snake = u[2]
  locust = f(t)
  lion = 1
  mouse = ϰ
  loris = ϑ(t)
  aardvark = σ
  deer = 1
  coati = 0
  snail = _fdm(f, t)
  crocodile = 0

  du[1] = (((mouse * loris - mouse * locust) - snail) * snake) / lion - ((deer + coati * locust) / 2) * ((aardvark * snake) / lion) ^ 2
  du[2] = ((mouse + crocodile / lion) * snake + (coati / (2lion)) * (aardvark * snake) ^ 2) - lion

  return nothing
end

let
  p = params()

  tspan = (0.0, 1.0)

  Δt = (tspan[2] - tspan[1])

  dt = 1.0 / 252.0

  nodes = Int(ceil(Δt / dt) + 1)

  # time meshes
  t = T = [tspan[1] + (i - 1) * dt for i = 1:nodes]

  prob = ODEProblem(riccati, zeros(2), (0.0, 1.0), p)

  prob_func = (prob, i, repeat) -> begin
    ODEProblem(prob.f, prob.u0, (T[i + 1], t[1]), prob.p)
  end

  odes = EnsembleProblem(prob, prob_func = prob_func)

  @time sol = DifferentialEquations.solve(
    odes,
    Tsit5(),
    trajectories = nodes - 1,
    saveat = -dt
  )
end

with the following time results for version v0.9.5:

22.494440 seconds (27.82 M allocations: 1.664 GiB, 2.84% gc time)
17.850900 seconds (20.12 M allocations: 1.294 GiB, 2.64% gc time)
18.231992 seconds (20.12 M allocations: 1.294 GiB, 3.24% gc time)
18.729749 seconds (20.12 M allocations: 1.294 GiB, 2.74% gc time)

As you suggest, it is better but it doesn't seem to solve the issue. Unless there is something I am missing.

Thanks!

[oxinabox]

Interesting, I thought it would help more than that.
Anyway #73 solves this which I will release as 0.9.6 in the next half hour or so

[rvignolo]

Thank you very much Lyndon!!