# Turing TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 5}

Hi all

I have the well-known dual numbers problem when using Turing with an ODE model. The code below throws `TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 5}`. From the previous posts on this problem I understand that converting the type of u0 to eltype(p) should solve the problem, but this didn’t work for me. Any help is appreciated.

``````# Household stuff
cd(@__DIR__)
using Pkg
Pkg.activate(".")

using DifferentialEquations
using StatsPlots
using Turing
using Distributions
using Random
using LabelledArrays

Random.seed!(42)

# ODE model: simple SIR model with seasonally forced contact rate

function SIR!(du,u,p,t)

# states
(S, I, R) = u[1:3]
N = S + I + R

# params
β = p.β
η = p.η
φ = p.φ
ω = 1.0/p.ω
μ = p.μ
σ = p.σ

# FOI
βeff = β * (1.0+η*cos(2.0*π*(t-φ)/365.0))
λ = βeff*I/N

# change in states
du[1] = (μ*N - λ*S - μ*S + ω*R)
du[2] = (λ*S - σ*I - μ*I)
du[3] = (σ*I - μ*R - ω*R)
du[4] = (σ*I) # cumulative incidence

end

# Solver settings

tmin = 0.0
tmax = 10.0*365.0
tspan = (tmin, tmax)
solvsettings = (abstol = 1.0e-3,
reltol = 1.0e-3,
saveat = 7.0,
solver = AutoTsit5(Rosenbrock23()))

# Initiate ODE problem
theta_fix =  [1.0/(80*365)]
theta_est =  [0.28, 0.07, 365.0, 1.0, 1.0/5.0]
p = @LArray [theta_est; theta_fix] (:β, :η, :ω, :φ, :σ, :μ)
u0 = @LArray [9998.0,1.0,1.0,1.0] (:S,:I,:R,:C)

# Initiate ODE problem
problem = ODEProblem(SIR!,u0,tspan,p)
sol = solve(problem,
solvsettings.solver,
abstol=solvsettings.abstol,
reltol=solvsettings.reltol,
isoutofdomain=(u,p,t)->any(x->x<0.0,u),
save_idxs=4,
saveat=solvsettings.saveat)

#plot(sol)
incidence = sol[2:end] - sol[1:(end-1)]
plot(incidence)

# Fake some data from model
data = rand.(Poisson.(incidence))

#plot(incidence, legend = false); scatter!(data,legend = false)

# Fit model to fake data

# Set up as Turing model

Turing.setprogress!(true)

@model turingmodel(data, p, u0, problem, solvsettings) = begin

# Priors
β ~ Uniform(0.0,1.0)
η ~ Uniform(0.0,1.0)
ω ~ Uniform(1.0, 5.0*365.0)
φ ~ Uniform(0.0,364.0)
σ ~ Uniform(0.0,1.0)
p[1:5] = [β,η,ω,φ,σ]

# Update problem and solve ODEs

problem_new = remake(problem, p=p, u0=eltype(p).(u0))
sol_new = solve(problem_new,
solvsettings.solver,
abstol=solvsettings.abstol,
reltol=solvsettings.reltol,
isoutofdomain=(u,p,t)->any(x->x<0.0,u),
save_idxs=4,
saveat=solvsettings.saveat)

incidence = sol_new[2:end] - sol_new[1:(end-1)]
incidence = ifelse.(incidence .< 0.0, 0.0, incidence) # to address numerical instability problem
l = size(incidence,1)
if l==size(data,1)
for i = 1:l
data[i] ~ Poisson(incidence[i])
end
else
return
end
end

model = turingmodel(data, p, u0, problem, solvsettings)

# Fit 3 independent chains with multithreading

@time trace = sample(model, NUTS(), MCMCThreads(), 1000, 3, progress=true)
``````

You can’t modify a vector of floats to change them to duals. Instead, just make a new Vector and you’re fine:

``````p = @LArray [β,η,ω,φ,σ,p.μ] (:β, :η, :ω, :φ, :σ, :μ)
``````

And you can use other tools like ArrayInterface.restructure etc. to cut down on code size or whatnot.

2 Likes

thanks @ChrisRackauckas, this works (as in no further dual numbers error). However, I now have a new problem. Will make a new post.