# Turing and DifferentialEquations: TypeError expected Float64, got ForwardDiff

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, 2, 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:

``````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…

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

You may want to look at PreallocationTools.jl.

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:

``````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)

``````

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

1 Like

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