# ODEs inside a Turing model

**URL:** <https://discourse.julialang.org/t/odes-inside-a-turing-model/39759>\
**Category:** New to Julia\
**Tags:** ode, turing\
**Created:** [May 19, 2020, 3:15pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759 "2020-05-19T15:15:07Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Manoj\_Gopalkrishnan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manoj_gopalkrishnan/32/14587_2.png) [@Manoj\_Gopalkrishnan](https://discourse.julialang.org/u/Manoj_Gopalkrishnan)\
**Post date:** [May 19, 2020, 3:15pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/1 "2020-05-19T15:15:07Z")

</div>

I’m trying to write an HMM using Turing. The forward model from one time step to another is given by a system of differential equations. I wrote this up but I’m getting an error:

```julia
MethodError: no method matching OrdinaryDiffEq.Vern9Cache(

```

Is this because I have to use Dual Numbers inside the ODE? Do I typecast the initial condition to be a vector of dual numbers?

Is there a minimal example of this sort of thing? None of the Turing tutorials cover working with an ODE within Turing. I could find one example here: [Simple example of ODE parameters inference · Issue #810 · TuringLang/Turing.jl · GitHub](https://github.com/TuringLang/Turing.jl/issues/810) of a predator-prey model whose parameters were estimated within Turing, but this case has one single ODE problem with one single initial condition. My situation requires a fresh ODE problem at every discrete time-step with fresh initial conditions.

Some kind of pointers to what to read, or minimal examples, would be very helpful!

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [May 19, 2020, 6:00pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/2 "2020-05-19T18:00:53Z")

</div>

Maybe @ChrisRackauckas might be able to help. I don’t think it is directly related to Turing but I might be wrong. Could you also post a MWE?

---

<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:** [May 19, 2020, 6:09pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/3 "2020-05-19T18:09:43Z")

</div>

We should get a quick tutorial into the Turing docs, and probably show off a few things. I think @Vaibhavdixit02 was looking to do it.

@Manoj_Gopalkrishnan you need to make sure that your state types match the dualness they will need, so probably just `_u0 = eltype(p).(u0)` to dualify them. And FWIW, it’s probably best to use `concrete_solve` and Zygote differentiation for adjoints if you have large numbers of parameters.

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [May 19, 2020, 6:23pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/4 "2020-05-19T18:23:39Z")

</div>

Yes, a tutorial would be great! I only have a very basic knowledge on ODEs, would be amazing if someone with more expertise could write up a tutorial.

---

<div class="post-metadata">

**Author:** ![Manoj\_Gopalkrishnan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manoj_gopalkrishnan/32/14587_2.png) [@Manoj\_Gopalkrishnan](https://discourse.julialang.org/u/Manoj_Gopalkrishnan)\
**Post date:** [May 20, 2020, 7:39am UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/5 "2020-05-20T07:39:45Z")

</div>

Thank you Martin and Chris!

Meanwhile I wrote this code snippet based on the document I had shared and was able to get it to work in my ODE solver loop.

```julia
function test_f(u0,q)
    prob = ODEProblem(PR!,eltype(p).(u0),(0.0,1.0), q)
    concrete_solve(prob,Rodas4(),saveat=[1.0];abstol=1e-6, reltol= 1e-6, alg_hints=[:nonstiff])
end

```

Here PR! is my function describing the vector field for the ODE and q are parameters. Though this worked, one thing I’m not clear on is: how does it know what is p? Is it picking up p from some global context, or does it know that p is the name to one of the fields returned by ODEProblem? Or should I be using eltype.(q).(u0)?

Now when I describe my Turing model and run a chain, I get an error message

```julia
Mutating arrays is not supported

```

What does this mean and how do I find the part of my code that is causing this, and fix it? I am not able to give a MWE short of giving my entire Turing model, because I am not sure where the error is.

Also: if I comment out the ODE code within my Turing model, reversediff works. If I include the ODE code, reversediff gives the above error. Zygote gives the above error even with or without the ODE code commented out.

So my question is: What should I be looking for to find mutating arrays? Does mutating mean the type of the array changes? Where could this be happening? How do I stop it from happening?

---

<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:** [May 20, 2020, 8:20am UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/6 "2020-05-20T08:20:55Z")

</div>

> [@Manoj\_Gopalkrishnan](#):
>
> Also: if I comment out the ODE code within my Turing model, reversediff works. If I include the ODE code, reversediff gives the above error. Zygote gives the above error even with or without the ODE code commented out.

Let me fix that really quickly. I forgot to add the `@grad` overload after @mohamed82008 added it to ReverseDiff.jl.

---

<div class="post-metadata">

**Author:** ![Manoj\_Gopalkrishnan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manoj_gopalkrishnan/32/14587_2.png) [@Manoj\_Gopalkrishnan](https://discourse.julialang.org/u/Manoj_Gopalkrishnan)\
**Post date:** [May 20, 2020, 3:28pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/7 "2020-05-20T15:28:27Z")

</div>

So I tried to write a MWE. Here it is, really minimal. It’s a stupid way of using differential equations to do linear regression. The differential equation propagates the y vector from one time to the next by adding m to y.

```julia
using Plots,Turing, Distributions, DifferentialEquations 
using XLSX, ReverseDiff, Zygote, DualNumbers, DiffEqSensitivity, OrdinaryDiffEq
Turing.setadbackend(:reversediff)

function line!(du,u,p,t)  
    du[1] = p[1]
end

function test_f(u0,m)
    prob = ODEProblem(line!,eltype(m).(u0),(0.0,1.0), m)
    solve(prob)
end

@model function fitline(ys, ::Type{T}=Vector{Float64}) where {T}
    θ ~ Uniform(0, .9*π/2) 
    m = tan(θ)  
    c ~ Normal(0,100)
    σ ~ truncated(Normal(0,100),0,Inf)
    τ = length(ys)
    
    u0 = Array{Float64}(undef, τ+1)
    u0[1] = c
    for i = 1:τ
       u0[i+1] = test_f([u0[i]],[m])
       ys[i] ~ Normal(u0[i+1], σ)
    end
end;

ys = [3.0 + rand(Normal(i*5.0, 0.1)) for i in 1:20]
model = fitline(ys)

chain = sample(model, NUTS(.65),30)

```

and the error is

```julia
MethodError: Cannot `convert` an object of type OrdinaryDiffEq.ODECompositeSolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(line!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}},OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true,typeof(line!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23Tableau{Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(line!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(line!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},DefaultLinSolve,FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:forward},Float64},FiniteDiff.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}}},DiffEqBase.DEStats} to an object of type Float64
Closest candidates are:
  convert(::Type{T}, !Matched::Unitful.Quantity) where T<:Real at C:\Users\Manoj\.julia\packages\Unitful\9XsxK\src\conversion.jl:145
  convert(::Type{T}, !Matched::Unitful.Level) where T<:Real at C:\Users\Manoj\.julia\packages\Unitful\9XsxK\src\logarithm.jl:22
  convert(::Type{T}, !Matched::Unitful.Gain) where T<:Real at C:\Users\Manoj\.julia\packages\Unitful\9XsxK\src\logarithm.jl:62
  ...

Stacktrace:
 [1] setindex!(::Array{Float64,1}, ::OrdinaryDiffEq.ODECompositeSolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(line!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}},OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true,typeof(line!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23Tableau{Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(line!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Array{Float64,1}},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(line!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},DefaultLinSolve,FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:forward},Float64},FiniteDiff.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}}},DiffEqBase.DEStats}, ::Int64) at .\array.jl:826
 [2] macro expansion at .\In[34]:24 [inlined]
 [3] ##evaluator#898(::DynamicPPL.Model{var"###evaluator#898",(:ys, :T),Tuple{Array{Float64,1},Type{Array{Float64,1}}},(),DynamicPPL.ModelGen{var"###generator#899",(:ys, :T),(:T,),Tuple{Type{Array{Float64,1}}}}}, ::DynamicPPL.VarInfo{DynamicPPL.Metadata{Dict{DynamicPPL.VarName,Int64},Array{Distribution,1},Array{DynamicPPL.VarName,1},Array{Real,1},Array{Set{DynamicPPL.Selector},1}},Float64}, ::DynamicPPL.SampleFromPrior, ::DynamicPPL.DefaultContext) at C:\Users\Manoj\.julia\packages\DynamicPPL\h0R38\src\compiler.jl:355
 [4] evaluate_singlethreaded at C:\Users\Manoj\.julia\packages\DynamicPPL\h0R38\src\model.jl:147 [inlined]
 [5] Model at C:\Users\Manoj\.julia\packages\DynamicPPL\h0R38\src\model.jl:136 [inlined]
 [6] VarInfo at C:\Users\Manoj\.julia\packages\DynamicPPL\h0R38\src\varinfo.jl:110 [inlined]
 [7] VarInfo at C:\Users\Manoj\.julia\packages\DynamicPPL\h0R38\src\varinfo.jl:109 [inlined]
 [8] DynamicPPL.Sampler(::NUTS{Turing.Core.ReverseDiffAD{false},(),AdvancedHMC.DiagEuclideanMetric}, ::DynamicPPL.Model{var"###evaluator#898",(:ys, :T),Tuple{Array{Float64,1},Type{Array{Float64,1}}},(),DynamicPPL.ModelGen{var"###generator#899",(:ys, :T),(:T,),Tuple{Type{Array{Float64,1}}}}}, ::DynamicPPL.Selector) at C:\Users\Manoj\.julia\packages\Turing\d4vqQ\src\inference\hmc.jl:378
 [9] Sampler at C:\Users\Manoj\.julia\packages\Turing\d4vqQ\src\inference\hmc.jl:376 [inlined]
 [10] sample(::Random._GLOBAL_RNG, ::DynamicPPL.Model{var"###evaluator#898",(:ys, :T),Tuple{Array{Float64,1},Type{Array{Float64,1}}},(),DynamicPPL.ModelGen{var"###generator#899",(:ys, :T),(:T,),Tuple{Type{Array{Float64,1}}}}}, ::NUTS{Turing.Core.ReverseDiffAD{false},(),AdvancedHMC.DiagEuclideanMetric}, ::Int64; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\Manoj\.julia\packages\Turing\d4vqQ\src\inference\Inference.jl:161
 [11] sample at C:\Users\Manoj\.julia\packages\Turing\d4vqQ\src\inference\Inference.jl:161 [inlined]
 [12] #sample#1 at C:\Users\Manoj\.julia\packages\Turing\d4vqQ\src\inference\Inference.jl:151 [inlined]
 [13] sample(::DynamicPPL.Model{var"###evaluator#898",(:ys, :T),Tuple{Array{Float64,1},Type{Array{Float64,1}}},(),DynamicPPL.ModelGen{var"###generator#899",(:ys, :T),(:T,),Tuple{Type{Array{Float64,1}}}}}, ::NUTS{Turing.Core.ReverseDiffAD{false},(),AdvancedHMC.DiagEuclideanMetric}, ::Int64) at C:\Users\Manoj\.julia\packages\Turing\d4vqQ\src\inference\Inference.jl:151
 [14] top-level scope at In[34]:31

```

I think the problem is inside the function test\_f when I am trying to cast u0 to the type of m which is Float64 here. I don’t know how to fix this. Please help!

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [May 20, 2020, 3:41pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/8 "2020-05-20T15:41:41Z")

</div>

No, I think the issue is that `u0` is `Array{Float64}`. @mohamed82008 do you know how to do this correctly?

– EDIT –  
It might work if you use this magical definition, but I can’t try it out atm.

```julia
@model function fitline(ys, ::Type{T}=Vector{Float64}, ::Type{V}=Vector{Float64}) where {T,V}
    θ ~ Uniform(0, .9*π/2) 
    m = tan(θ)  
    c ~ Normal(0,100)
    σ ~ truncated(Normal(0,100),0,Inf)
    τ = length(ys)
    
    u0 = V(undef, τ+1)
    u0[1] = c
    ...

```

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [May 20, 2020, 3:47pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/9 "2020-05-20T15:47:58Z")

</div>

The problem seems to be that `solve` returns an instance of `OrdinaryDiffEq.ODECompositeSolution` but @Manoj_Gopalkrishnan is trying to store this in a vector of `Float64`.

---

<div class="post-metadata">

**Author:** ![Manoj\_Gopalkrishnan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manoj_gopalkrishnan/32/14587_2.png) [@Manoj\_Gopalkrishnan](https://discourse.julialang.org/u/Manoj_Gopalkrishnan)\
**Post date:** [May 20, 2020, 3:55pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/10 "2020-05-20T15:55:13Z")

</div>

Tried

> [@trappmartin](#):
>
> @model function fitline(ys, ::Type{T}=Vector{Float64}, ::Type{V}=Vector{Float64}) where {T,V}  
> θ ~ Uniform(0, .9\*π/2)  
> m = tan(θ)  
> c ~ Normal(0,100)  
> σ ~ truncated(Normal(0,100),0,Inf)  
> τ = length(ys)
> 
> ```
> u0 = V(undef, τ+1)
> u0[1] = c
> ...
> 
> ```

out. Got the error.

```julia
syntax: function argument names not unique

Stacktrace:
 [1] top-level scope at C:\Users\Manoj\.julia\packages\IJulia\DrVMH\src\kernel.jl:52

```

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [May 20, 2020, 4:02pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/11 "2020-05-20T16:02:47Z")

</div>

@Manoj_Gopalkrishnan I suggest you first write a function that samples parameters from the prior distributions and then generates fake data using `rand` then move on to making a Turing model out of it. The errors above don’t seem Turing-specific but are more basic Julia errors that can be fixed in the fake data generating function.

---

<div class="post-metadata">

**Author:** ![Manoj\_Gopalkrishnan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manoj_gopalkrishnan/32/14587_2.png) [@Manoj\_Gopalkrishnan](https://discourse.julialang.org/u/Manoj_Gopalkrishnan)\
**Post date:** [May 20, 2020, 5:01pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/13 "2020-05-20T17:01:23Z")

</div>

> [@mohamed82008](#):
>
> @Manoj_Gopalkrishnan I suggest you first write a function that samples parameters from the prior distributions and then generates fake data using `rand` then move on to making a Turing model out of it. The errors above don’t seem Turing-specific but are more basic Julia errors that can be fixed in the fake data generating function.

Thank you @mohamed82008 . I needed to write

```julia
u0 = test_f([u0],[m]).u[end][1]

```

In particular, I didn’t need u0 to be an array. With this I am able to generate data, as you suggested. This is also a cool debugging technique for Turing models! I feel I learnt something.

BUT when I am trying to run the Turing model, the chain starts but doesn’t terminate even after a long time. Given how small this model is, and how few samples I’m trying to generate, clearly I am doing something wrong. Here’s my code as it now stands.

```julia
using Turing, Distributions, DifferentialEquations 
using ReverseDiff, Zygote, DualNumbers, DiffEqSensitivity, OrdinaryDiffEq
Turing.setadbackend(:reversediff)

function line!(du,u,p,t)  
    du[1] = p[1]
end

function test_f(u0,m)
    prob = ODEProblem(line!,eltype(m).(u0),(0.0,1.0), m)
    solve(prob,saveat=[1.0], alg_hints=[:nonstiff])
end

@model function fitline(ys, ::Type{T}=Vector{Float64}) where {T}
    θ ~ Uniform(0, .9*π/2.0) 
    m = tan(θ)  
    c ~ Normal(0,100.0)
    σ ~ truncated(Normal(0,100.0),0,Inf)
    τ = length(ys)
    
    u0 = Array{Float64}(undef, τ+1)
    u0[1] = c
    for i = 1:τ
        u0[i+1] = test_f([u0[i]],[m]).u[end][1]
        ys[i] ~ Normal(u0[i+1], σ)
    end
end;

ys = [3.0 + rand(Normal(i*5.0, 0.1)) for i in 1:20]
model = fitline(ys)

chain = sample(model, NUTS(.65),30)

```

---

<div class="post-metadata">

**Author:** ![Manoj\_Gopalkrishnan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manoj_gopalkrishnan/32/14587_2.png) [@Manoj\_Gopalkrishnan](https://discourse.julialang.org/u/Manoj_Gopalkrishnan)\
**Post date:** [May 20, 2020, 5:09pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/14 "2020-05-20T17:09:15Z")

</div>

So after a longish while it halted, and gave reasonable output:

```julia
Object of type Chains, with data of type 15×15×1 Array{Float64,3}

Iterations = 1:15
Thinning interval = 1
Chains = 1
Samples per chain = 15
internals = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, max_hamiltonian_energy_error, n_steps, nom_step_size, numerical_error, step_size, tree_depth
parameters = c, θ, σ

2-element Array{ChainDataFrame,1}

Summary Statistics
  parameters mean std naive_se mcse ess r_hat
  ────────── ─────── ─────── ──────── ─────── ────── ──────
           c 0.6625 0.3308 0.0854 missing 2.8000 2.5974
           θ 1.3484 0.0399 0.0103 missing 2.8000 2.5303
           σ 11.5305 11.5763 2.9890 missing 2.8000 2.0529

Quantiles
  parameters 2.5% 25.0% 50.0% 75.0% 97.5%
  ────────── ────── ────── ────── ─────── ───────
           c 0.1994 0.3262 0.6950 0.9553 1.0723
           θ 1.2897 1.3043 1.3734 1.3818 1.3825
           σ 2.4997 2.8359 3.8546 19.2765 31.9908

```

Given how few samples I took, theta looks a good fit for the slope (tan \theta looks close to actual slope) but \sigma and c are a bit off, which is not surprising actually.

I’m still wondering: why is this so slow? And how can I speed it up? Any tips?

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [May 20, 2020, 6:21pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/15 "2020-05-20T18:21:49Z")

</div>

Which AD backend are you using?

---

<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:** [May 20, 2020, 9:14pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/16 "2020-05-20T21:14:37Z")

</div>

With the next release of DiffEqBase (coming tonight), the other reverse modes will work on `concrete_solve`, so this is valid code:

```julia
using Turing, Distributions, DifferentialEquations
using ReverseDiff, Zygote, DualNumbers, DiffEqSensitivity, OrdinaryDiffEq
Turing.setadbackend(:reversediff)

function line!(du,u,p,t)
    du[1] = p[1]
end

function test_f(u0,m)
    prob = ODEProblem(line!,u0,(0.0,1.0), m)
    concrete_solve(prob, nothing, save_everystep=false, alg_hints=[:nonstiff])
end

@model function fitline(ys, ::Type{T}=Vector{Float64}) where {T}
    θ ~ Uniform(0, .9*π/2.0)
    m = tan(θ)
    c ~ Normal(0,100.0)
    σ ~ truncated(Normal(0,100.0),0,Inf)
    τ = length(ys)

    u0 = c
    for i = 1:τ
        u0 = test_f([u0],m)
        ys[i] ~ Normal(u0, σ)
    end
end;

ys = [3.0 + rand(Normal(i*5.0, 0.1)) for i in 1:20]
model = fitline(ys)

chain = sample(model, NUTS(.65),30)

```

However, you only want to use reverse mode AD when there are a lot of parameters. For something like \<100 parameters, it’s probably better to do the following:

```julia
Turing.setadbackend(:forwarddiff)

function line!(du,u,p,t)
    du[1] = p[1]
end

function test_f(u0,m)
    prob = ODEProblem(line!,eltype(m).(u0),(0.0,1.0), m)
    concrete_solve(prob, nothing, save_everystep=false, alg_hints=[:nonstiff])[end][1]
end

@model function fitline(ys, ::Type{T}=Vector{Float64}) where {T}
    θ ~ Uniform(0, .9*π/2.0)
    m = tan(θ)
    c ~ Normal(0,100.0)
    σ ~ truncated(Normal(0,100.0),0,Inf)
    τ = length(ys)

    u0 = c
    for i = 1:τ
        u0 = test_f([u0],m)
        ys[i] ~ Normal(u0, σ)
    end
end;

ys = [3.0 + rand(Normal(i*5.0, 0.1)) for i in 1:20]
model = fitline(ys)
chain = sample(model, NUTS(.65),30)

```

Still, this is not optimized all of the way because we can specialize on the fact that you’re solving a scalar ODE, and thus you can use an out-of-place ODE function with a scalar initial condition like:

```julia
Turing.setadbackend(:forwarddiff)

function line_oop!(u,p,t)
    p[1]
end

function test_f(u0,m)
    prob = ODEProblem(line_oop!,eltype(m).(u0),(0.0,1.0), m)
    concrete_solve(prob, nothing, save_everystep=false, alg_hints=[:nonstiff])[end]
end

@model function fitline(ys, ::Type{T}=Vector{Float64}) where {T}
    θ ~ Uniform(0, .9*π/2.0)
    m = tan(θ)
    c ~ Normal(0,100.0)
    σ ~ truncated(Normal(0,100.0),0,Inf)
    τ = length(ys)

    u0 = c
    for i = 1:τ
        u0 = test_f(u0,m)
        ys[i] ~ Normal(u0, σ)
    end
end;

ys = [3.0 + rand(Normal(i*5.0, 0.1)) for i in 1:20]
model = fitline(ys)
chain = sample(model, NUTS(.65),30)

```

So those examples should set you on the right track.

---

<div class="post-metadata">

**Author:** ![Manoj\_Gopalkrishnan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manoj_gopalkrishnan/32/14587_2.png) [@Manoj\_Gopalkrishnan](https://discourse.julialang.org/u/Manoj_Gopalkrishnan)\
**Post date:** [May 21, 2020, 6:38am UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/17 "2020-05-21T06:38:00Z")

</div>

I was using ReverseDiff, am going to try with the others as well…

---

<div class="post-metadata">

**Author:** ![Manoj\_Gopalkrishnan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manoj_gopalkrishnan/32/14587_2.png) [@Manoj\_Gopalkrishnan](https://discourse.julialang.org/u/Manoj_Gopalkrishnan)\
**Post date:** [May 21, 2020, 6:52am UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/18 "2020-05-21T06:52:29Z")

</div>

Using :forwarddiff, now it’s lightning fast! I still couldn’t get concrete\_solve to work with any of the solvers. Maybe your update is not pushed yet.

Had to modify one line in your code:  
from

```julia
concrete_solve(prob, nothing, save_everystep=false, alg_hints=[:nonstiff])[end]

```

to

```julia
solve(prob, nothing, save_everystep=false, alg_hints=[:nonstiff]).u[end][1]

```

---

<div class="post-metadata">

**Author:** ![Manoj\_Gopalkrishnan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/manoj_gopalkrishnan/32/14587_2.png) [@Manoj\_Gopalkrishnan](https://discourse.julialang.org/u/Manoj_Gopalkrishnan)\
**Post date:** [May 21, 2020, 6:55am UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/19 "2020-05-21T06:55:53Z")

</div>

After just 500 samples, for m=5, c=33, and \sigma = .1, we get

```julia
Object of type Chains, with data of type 250×15×1 Array{Float64,3}

Iterations = 1:250
Thinning interval = 1
Chains = 1
Samples per chain = 250
internals = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, max_hamiltonian_energy_error, n_steps, nom_step_size, numerical_error, step_size, tree_depth
parameters = c, θ, σ

2-element Array{ChainDataFrame,1}

Summary Statistics
  parameters mean std naive_se mcse ess r_hat
  ────────── ─────── ────── ──────── ────── ──────── ──────
           c 33.0087 0.0396 0.0025 0.0014 173.6316 0.9968
           θ 1.3732 0.0001 0.0000 0.0000 146.2767 1.0001
           σ 0.0999 0.0173 0.0011 0.0027 206.8505 1.0022

Quantiles
  parameters 2.5% 25.0% 50.0% 75.0% 97.5%
  ────────── ─────── ─────── ─────── ─────── ───────
           c 32.9347 32.9827 33.0057 33.0371 33.0864
           θ 1.3730 1.3732 1.3733 1.3733 1.3735
           σ 0.0732 0.0872 0.0977 0.1091 0.1391

```

which is bang on the dot! And took a twinkle to run!

---

<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:** [May 21, 2020, 12:25pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/20 "2020-05-21T12:25:33Z")

</div>

> [@Manoj\_Gopalkrishnan](#):
>
> Had to modify one line in your code:

DiffEqBase v6.35. But for ForwardDiff you don’t need concrete\_solve.

---

<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:** [May 21, 2020, 3:45pm UTC](https://discourse.julialang.org/t/odes-inside-a-turing-model/39759/21 "2020-05-21T15:45:15Z")

</div>

There will be a tutorial in the Turing docs soon enough covering this: [Add Bayesian differential equations tutorial by ChrisRackauckas · Pull Request #68 · TuringLang/TuringTutorials · GitHub](https://github.com/TuringLang/TuringTutorials/pull/68)