# Turing and DifferentialEquations: TypeError expected Float64, got ForwardDiff

**URL:** <https://discourse.julialang.org/t/turing-and-differentialequations-typeerror-expected-float64-got-forwarddiff/82194>\
**Category:** Probabilistic Programming\
**Tags:** question, package, error\
**Created:** [June 3, 2022, 2:45pm UTC](https://discourse.julialang.org/t/turing-and-differentialequations-typeerror-expected-float64-got-forwarddiff/82194 "2022-06-03T14:45:33Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![finn](https://avatars.discourse-cdn.com/v4/letter/f/e68b1a/32.png) [@finn](https://discourse.julialang.org/u/finn)\
**Post date:** [June 3, 2022, 2:45pm UTC](https://discourse.julialang.org/t/turing-and-differentialequations-typeerror-expected-float64-got-forwarddiff/82194/1 "2022-06-03T14:45:33Z")

</div>

Hi all,

I am trying to use DifferentialEquations.jl and Turing.jl to perform parameter estimation using MCMC.  
I am currently getting the error

> TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 3}

I have seen other people get similar errors (e.g. [1](https://discourse.julialang.org/t/typeerror-using-odes-in-turing-jl/34857), [2](https://discourse.julialang.org/t/turing-type-error-expected-float64-got-forwarddiff-dual/47670), can’t link to others because I am new 😅), and I have tried to use the suggested fixes, but I can’t seem to get it to work.

The problem seems to be caused by Turing’s automatic differentiation. I have tried to use a different AD backend (Zygote), but that did not help. I tried to turn off the AD used by the DifferentialEquations solver, but that did not help either.

Does anybody have any idea what could be going wrong and how I might fix it?

Here is a simplified version of my code:

```julia
using Turing, DifferentialEquations, LinearAlgebra

# Known Constants
observation_size = 100
Δ = 0.05
time = [Δ * i for i in 0:observation_size]

slices_per_block = 10
total_slices = 3 * slices_per_block
mid_of_block = ceil(Int64, slices_per_block / 2)
observed_nodes = [i * slices_per_block + mid_of_block for i in 0:2]
true_ks_int = [20., 30., 25.]

T_0 = Vector{Float64}(undef, total_slices)
T_0[1:slices_per_block] .= 330.
T_0[(slices_per_block + 1):(2 * slices_per_block)] .= 270.
T_0[(2 * slices_per_block + 1):(3 * slices_per_block)] .= 310.

# Unknown Constants
true_k12 = 5.
true_k23 = 4.
true_σ = 2.

true_p = [slices_per_block, true_ks_int, true_k12, true_k23]
function LSSM(dT, T, p, t)
    s, ks_int, k12, k23 = p
    N_s = 3 * s
    # Put conduction coefficient between adjacent slices in list
    ks = Vector{Any}(undef, N_s - 1) #I do "Any" here since k12, k23 might be ForwardDiff and not Float
    ks[1:(s - 1)] .= ks_int[1]
    ks[s] = k12
    ks[(s + 1):(2 * s - 1)] .= ks_int[2]
    ks[2 * s] = k23
    ks[(2 * s + 1):(3 * s - 1)] .= ks_int[3]
    # Conduction + Convection
    dT[1] = (ks[1] * (T[2] - T[1]))
    for i ∈ 2:(N_s - 1)
        dT[i] = (ks[i - 1] * (T[i - 1] - T[i]) + ks[i] * (T[i + 1] - T[i]))
    end
    dT[N_s] = (ks[N_s - 1] * (T[N_s - 1] - T[N_s]))
end
LSSM_dynamics = ODEProblem(LSSM, T_0, (time[1], time[end]), true_p)
true_sol = solve(LSSM_dynamics, Tsit5(); saveat = Δ)
T_obs = true_sol[observed_nodes, :]
y = Array(T_obs) + true_σ * randn(size(Array(T_obs)))

@model function fit_LSSM(data, system, syst_consts)
    σ ~ Gamma(1, 1) 
    k12 ~ Gamma(1, 1)
    k23 ~ Gamma(1, 1)
    s, ks_int, Δ = syst_consts
    p = [s, ks_int, k12, k23]
    predicted = solve(system, Tsit5(); p = p, saveat = Δ)
    for i ∈ 1:length(predicted)
        data[:, i] ~ MvNormal(predicted[observed_nodes, i], σ^2 * I)
    end
end

model = fit_LSSM(y, LSSM_dynamics, [slices_per_block, true_ks_int, Δ]);

chain = sample(model, NUTS(0.65), MCMCSerial(), 1000, 1; progress = false) 

```

The stacktrace is very big and I don’t really understand what is going on… 😅

---

<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 3, 2022, 3:03pm UTC](https://discourse.julialang.org/t/turing-and-differentialequations-typeerror-expected-float64-got-forwarddiff/82194/2 "2022-06-03T15:03:44Z")

</div>

> [@finn](#):
>
> `ks = Vector{Any}(undef, N_s - 1)`

Just use `promote_type(eltype(T),eltype(t))` for the type.

You may want to look at PreallocationTools.jl.

> **[GitHub - SciML/PreallocationTools.jl: Tools for building non-allocating...](https://github.com/SciML/PreallocationTools.jl)**
>
> Tools for building non-allocating pre-cached functions in Julia, allowing for GC-free usage of automatic differentiation in complex codes - GitHub - SciML/PreallocationTools.jl: Tools for building ...

---

<div class="post-metadata">

**Author:** ![finn](https://avatars.discourse-cdn.com/v4/letter/f/e68b1a/32.png) [@finn](https://discourse.julialang.org/u/finn)\
**Post date:** [June 4, 2022, 7:18pm UTC](https://discourse.julialang.org/t/turing-and-differentialequations-typeerror-expected-float64-got-forwarddiff/82194/3 "2022-06-04T19:18:59Z")

</div>

Thanks for your reply.  
I have had a look at PreallocationTools.jl, and tried to implement your suggestion about `promote_type`.

Unfortunately, I still get a similar error.  
I have added the updated version of my code below. I now dualcache my vector, and I get the cached version corresponding to the type of `T` when I use them in the ODEs. At the line where dT[1] is calculated, I get the following error:

> TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 3}

I have also tried getting the cached vector corresponding to the type of `dT`, but that gives the same error.

Here is my updated code:

```julia
using Turing, DifferentialEquations, LinearAlgebra, PreallocationTools

# Known Constants
observation_size = 100
Δ = 0.05
time = [Δ * i for i in 0:observation_size]

slices_per_block = 10
total_slices = 3 * slices_per_block
mid_of_block = ceil(Int64, slices_per_block / 2)
observed_nodes = [i * slices_per_block + mid_of_block for i in 0:2]
true_ks_int = [20., 30., 25.]

T_0 = Vector{Float64}(undef, total_slices)
T_0[1:slices_per_block] .= 330.
T_0[(slices_per_block + 1):(2 * slices_per_block)] .= 270.
T_0[(2 * slices_per_block + 1):(3 * slices_per_block)] .= 310.

# Unknown Constants
true_k12 = 5.
true_k23 = 4.
true_σ = 2.

# Caches for ForwardDiff
function cache_ks(ks_int, k12, k23, s)
    ks_type = promote_type(eltype(ks_int), typeof(k12), typeof(k23))
    ks = Vector{ks_type}(undef, 3 * s - 1)
    ks[1:(s - 1)] .= ks_int[1]
    ks[s] = k12
    ks[(s + 1):(2 * s - 1)] .= ks_int[2]
    ks[2 * s] = k23
    ks[(2 * s + 1):(3 * s - 1)] .= ks_int[3]
    ksd = dualcache(ks)
    return ksd
end
true_ksd = cache_ks(true_ks_int, true_k12, true_k23, slices_per_block)

true_p = [slices_per_block, true_ksd]
function LSSM(dT, T, p, t)
    s, ks = p
    ks = get_tmp(ks, T)
    N_s = 3 * s
    # Conduction
    dT[1] = (ks[1] * (T[2] - T[1]))
    for i ∈ 2:(N_s - 1)
        dT[i] = (ks[i - 1] * (T[i - 1] - T[i]) + ks[i] * (T[i + 1] - T[i]))
    end
    dT[N_s] = (ks[N_s - 1] * (T[N_s - 1] - T[N_s]))
end
LSSM_dynamics = ODEProblem(LSSM, T_0, (time[1], time[end]), true_p)
true_sol = solve(LSSM_dynamics, Tsit5(); saveat = Δ)
T_obs = true_sol[observed_nodes, :]
y = Array(T_obs) + true_σ * randn(size(Array(T_obs)))

@model function fit_LSSM(data, system, syst_consts)
    σ ~ Gamma(1, 1) 
    k12 ~ Gamma(1, 1)
    k23 ~ Gamma(1, 1)
    s, ks_int, Δ = syst_consts
    ksd = cache_ks(ks_int, k12, k23, s)
    p = [s, ksd]
    predicted = solve(system, Tsit5(); p = p, saveat = Δ)
    for i ∈ 1:length(predicted)
        data[:, i] ~ MvNormal(predicted[observed_nodes, i], σ^2 * I)
    end
end

model = fit_LSSM(y, LSSM_dynamics, [slices_per_block, true_ks_int, Δ]);

chain = sample(model, NUTS(0.65), MCMCSerial(), 1000, 1; progress = false) 

```

---

<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 23, 2022, 1:01am UTC](https://discourse.julialang.org/t/turing-and-differentialequations-typeerror-expected-float64-got-forwarddiff/82194/4 "2022-06-23T01:01:35Z")

</div>

I figured out a nice trick to make your code work with no fuss. The code from your original example will now work when DiffEqBase v6.91.6 is released (probably tomorrow), so no PreallocationTools.jl or anything required (it wasn’t being used correctly anyways 😅). Sorry for the late response but I thought I might as well find the real way to handle it. Cheers

---

<div class="post-metadata">

**Author:** ![finn](https://avatars.discourse-cdn.com/v4/letter/f/e68b1a/32.png) [@finn](https://discourse.julialang.org/u/finn)\
**Post date:** [June 23, 2022, 9:31am UTC](https://discourse.julialang.org/t/turing-and-differentialequations-typeerror-expected-float64-got-forwarddiff/82194/5 "2022-06-23T09:31:41Z")

</div>

Thank you! I updated my packages, and now both my mwe and my real code work.