Turing and DifferentialEquations: TypeError expected Float64, got ForwardDiff

Specific Domains Probabilistic programming
, ,
1

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 :sweat_smile:), 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… :sweat_smile:

2

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

You may want to look at PreallocationTools.jl.

3

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:

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

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 :sweat_smile:). Sorry for the late response but I thought I might as well find the real way to handle it. Cheers

1 Like
5

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