# Slow hierarchical ODE model with Turing

**URL:** <https://discourse.julialang.org/t/slow-hierarchical-ode-model-with-turing/63145>\
**Category:** Probabilistic Programming\
**Tags:** turing\
**Created:** [June 18, 2021, 3:18pm UTC](https://discourse.julialang.org/t/slow-hierarchical-ode-model-with-turing/63145 "2021-06-18T15:18:45Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![PavanChaggar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pavanchaggar/32/21597_2.png) [@PavanChaggar](https://discourse.julialang.org/u/PavanChaggar)\
**Post date:** [June 18, 2021, 3:18pm UTC](https://discourse.julialang.org/t/slow-hierarchical-ode-model-with-turing/63145/1 "2021-06-18T15:18:45Z")

</div>

Hello,

I’m having a great deal of trouble with the current model I’m working on. I’ve talked about it in slack but still struggling so apologies for cross posting.

It’s an hierarchical model with 30 subjects, each of whom have between 3-5 time points with 66 data points per time. The model I have is given below:

```julia
function NetworkFKPP(u, p, t)
    κ, α = p 
    du = -κ * L * u .+ α .* u .* (1 .- u)
end

@model function hierarchical_FKPP_NCP(data, initial_conditions, time, scans, prob)

    n = size(data)[2]

    σ ~ InverseGamma(2,3)

    κₘ ~ truncated(Normal(), 0, Inf)
    κₛ ~ truncated(Normal(), 0, Inf)

    αₘ ~ Normal()
    αₛ ~ truncated(Normal(), 0, Inf)

    κ ~ filldist(truncated(Normal(), 0, Inf), n)
    α ~ filldist(Normal(), n)
    
    k = (κ .* κₛ) .+ κₘ 
    a = (α .* αₛ) .+ αₘ

    @threads for i in 1:n
        prob_n = remake(prob, u0 = initial_conditions[:,i], p = [k[i], a[i]])
        predicted = solve(prob_n, Tsit5(), saveat=time[1:scans[i],i])
        cortical_preds = predicted[cortical_nodes,:]
        Turing.@addlogprob! loglikelihood(MvNormal(vec(cortical_preds), σ), data[1:66*scans[i],i])
    end
end

model = hierarchical_FKPP_NCP(data, initial_conditions, time, pos_scans, problem)
posterior = sample(model, NUTS(.65), 1_000)

```

`data` is a `330 x 30` array; `initial conditions` is a `83 x 78` array; `time` is a `5 x 30` array, `scans` is a `30 x 1` vector and problem is an ODEProblem set up with dummy variables. The ODE is a FKPP (diffusion plus logistic growth) model on a 83x83 undirected graph.

Trying to get 1\_000 samples, the ETA seems to blow up to be in the order of days. I’ve tried to make as many optimisations as I can using the docs and posts on slack/discourse/github etc but it’s still taking a very long time. I’ve also tried with reverse diff but this doesn’t seem to help. I’m trying to reproduce a collaborators results – they have used pymc3 on a similar model which they ran in 1.5-2 hrs on their laptop. The main differences between mine and their model are that they’re using explicit Euler for the ODE integration, they have 78 subjects, are using a centred parameterisation and have half Cauchy priors on their noise and hierarchical s.d priors. I can’t think of why any of these differences would have a significant impact on the Turing model (apart from maybe taking longer with more subjects). I’ve tested using Euler integration with dt=0.1 (as my collaborator had set) but this causes instability in the model and it crashes. Using Euler with a smaller time step has allocations than using Tsit5.

I’m running Julia 1.6 on a linux machine and using 32 threads. Env details:

```julia
      Status `~/Projects/TauPet/Project.toml`
  [6e4b80f9] BenchmarkTools v1.0.0
  [41bf760c] DiffEqSensitivity v6.48.0
  [0c46a032] DifferentialEquations v6.17.1
  [31c24e10] Distributions v0.25.2
  [093fc24a] LightGraphs v1.3.5
  [23992714] MAT v0.10.1
  [c7f686f2] MCMCChains v4.12.0
  [6fafb56a] Memoization v0.1.11
  [91a5bcdd] Plots v1.15.3
  [c3e4b0f8] Pluto v0.14.7
  [7f904dfe] PlutoUI v0.7.9
  [37e2e3b7] ReverseDiff v1.9.0
  [47aef6b3] SimpleWeightedGraphs v1.1.1
  [f3b207a7] StatsPlots v0.14.21
  [fce5fe82] Turing v0.16.0
  [e88e6eb3] Zygote v0.6.12

```

Any help would be greatly appreciated.

Thanks,  
Pavan

---

<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:** [June 18, 2021, 3:59pm UTC](https://discourse.julialang.org/t/slow-hierarchical-ode-model-with-turing/63145/2 "2021-06-18T15:59:19Z")

</div>

> [@PavanChaggar](#):
>
> `du = -κ * L * u .+ α .* u .* (1 .- u)`

The way you have implemented your differential equation is the slow way. You’ll want to use in-place functions. See [https://tutorials.sciml.ai/html/introduction/03-optimizing\_diffeq\_code.html](https://tutorials.sciml.ai/html/introduction/03-optimizing_diffeq_code.html) for some help getting started.

---

<div class="post-metadata">

**Author:** ![PavanChaggar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pavanchaggar/32/21597_2.png) [@PavanChaggar](https://discourse.julialang.org/u/PavanChaggar)\
**Post date:** [June 18, 2021, 5:00pm UTC](https://discourse.julialang.org/t/slow-hierarchical-ode-model-with-turing/63145/3 "2021-06-18T17:00:40Z")

</div>

Thanks! I got a big speed up in the ode solve using an in-place function. Also changing my L matrix from a sparse array to dense array also makes it faster. I tried static vectors and broadcast fusion using `@.` but these benchmarked slower than the in-place function. So I ended up with just this:

```julia
function NetworkFKPP4!(du, u, p, t) 
    κ, α = p 
    du = -κ * L * u .+ α .* u .* (1 .- u)
end

```

I’ll see how this impacts inference speed now!

---

<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:** [June 18, 2021, 7:15pm UTC](https://discourse.julialang.org/t/slow-hierarchical-ode-model-with-turing/63145/4 "2021-06-18T19:15:03Z")

</div>

> [@PavanChaggar](#):
>
> `du = -κ * L * u .+ α .* u .* (1 .- u)`

That’s not correct, you’ll notice your derivatives are zero like that. It should be `du .= -κ * L * u .+ α .* u .* (1 .- u)`, the `.=` for mutation.

---

<div class="post-metadata">

**Author:** ![PavanChaggar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pavanchaggar/32/21597_2.png) [@PavanChaggar](https://discourse.julialang.org/u/PavanChaggar)\
**Post date:** [June 18, 2021, 8:03pm UTC](https://discourse.julialang.org/t/slow-hierarchical-ode-model-with-turing/63145/5 "2021-06-18T20:03:51Z")

</div>

Thank you for the correction, you’re absolutely right – my bad!