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)
Turing.setadbackend(:forwarddiff)
@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
Turing.@addlogprob! -Inf
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)
```