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

New to Julia
, ,
1

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

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
3

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