# TypeError in Julia/Turing when sampling for a forced differential equation

Hi

I am new to Julia and Turing and am trying to fit a forced 0-D box ODE to data, but I get type error when doing sampling.

The tutorial (Bayesian Estimation of Differential Equations) gives an example for unforced LV model. Following this page (julia - solve system of ODEs with read in external forcing - Stack Overflow), I added an interpolation handle of the forcing as a parameter of the ODE, and the ODE solver runs without problems.

When I combine the forced ODE with Turing, however, I get a type error during HMC’s gradient calculation, quoting ”TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 3}”.

Below are my codes that should reproduce the error I got.

I appreciate any help, Thanks.

Below are my codes:

``````using Interpolations, DifferentialEquations, Plots, Turing, LinearAlgebra

# Define ODE
function f(du,u,p,t)
α, β, F = p
du[1] = α * u[1] + β * F(t) # Interpolations.jl interpolates via ()
end

# Define initial value problem
begin
# Define forcing
time_forcing = -1.:9.
data_forcing = [10,0,0,0,0,0,0,0,0, 0, 0]
F = Interpolations.scale(interpolate(data_forcing, BSpline(Linear())), time_forcing)
α = -0.5
β = 1
p_lin = (α, β, F)

# Define u0 and tspan
u0 = [0.]
tspan = (-1.,9.) # Note, that we would need to extrapolate beyond
ode_lin = ODEProblem(f,u0,tspan,p_lin)
end

# Generate data
sol  = solve(ode_lin, Tsit5(); saveat=1)
data = Array(sol) + 0.2 * randn(size(Array(sol)))

@model function fit_simple_box(data, F, ode_lin)
# Prior distributions.
σ ~ InverseGamma(2, 3)
α ~ Normal(0, 3)
β ~ Normal(0, 3)

# Simulate Lotka-Volterra model.
p = (α, β, F)
predicted = solve(ode_lin, Tsit5(); p=p, saveat=1)

# Observations.
for i in 1:length(predicted)
data ~ Normal(predicted[i][1], σ^2)
end

return nothing
end

model = fit_simple_box(data, F, ode_lin)
chain = sample(model, NUTS(0.65), MCMCSerial(), 1000, 2)
``````

I tried to run your code, but I received the following error:

``````ERROR: UndefVarError: g1_lin not defined
Stacktrace:
[1] top-level scope
@ Untitled-2:47
``````

I might be able to provide more guidance if given g1_lin or the full error message.

Automatic differentiation uses the Dual type to compute gradients from your code. The error you found typically happens when the type of an array is not generic. In this case, it is forced to be Float64, making a Dual type impossible. One possible culprit is

``````    data_forcing = [10,0,0,0,0,0,0,0,0, 0, 0]
F = Interpolations.scale(interpolate(data_forcing, BSpline(Linear())), time_forcing)
``````

If F always returns a Float64 regardless of the input type, this might be the cause of your problem.

You need to do `u0 = typeof(α).(ode_lin.u0)` to promote it to dual. But I have an idea for how to automate this promotion, so I’ll instead respond with a PR in a second.

2 Likes

Hi Christopher, Thanks for your quick response. the g1_lin was a bug, which I have fixed by replacing them with F, which stands for forcing. That said, I still get the same error message regarding type errors. The error message is quite long, copy the entire message goes beyond the length limit, but I will copy until [4] in the stack trace, so that it would help you to confirm if the error I got is reproduced.

If interpolation is the problem, is there any way to walk around? I can realize the same Bayesian model in Stan, but Julia has such a powerful differential equation tool box, which makes it much favorable in a long run.

``````ERROR: TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 3}
Stacktrace:
[1] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{Turing.Essential.var"#f#4"{DynamicPPL.TypedVarInfo{NamedTuple{(:σ, :α, :β), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:σ, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:σ, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:α, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:α, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:β, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:β, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(fit_simple_box), (:data, :F, :ode_lin), (), (), Tuple{Matrix{Float64}, ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Int64, ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 3}, i1::Int64)
@ Base ./array.jl:903
[2] f(du::Vector{Float64}, u::Vector{Float64}, p::Tuple{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Essential.var"#f#4"{DynamicPPL.TypedVarInfo{NamedTuple{(:σ, :α, :β), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:σ, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:σ, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:α, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:α, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:β, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:β, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(fit_simple_box), (:data, :F, :ode_lin), (), (), Tuple{Matrix{Float64}, ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Int64, ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 3}, ForwardDiff.Dual{ForwardDiff.Tag{Turing.Essential.var"#f#4"{DynamicPPL.TypedVarInfo{NamedTuple{(:σ, :α, :β), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:σ, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:σ, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:α, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:α, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:β, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:β, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(fit_simple_box), (:data, :F, :ode_lin), (), (), Tuple{Matrix{Float64}, ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Int64, ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 3}, ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}, t::Float64)
@ Main ~/Dropbox/Git_code/Deglacial_PowerLaw/scripts/Run_ODE_with_forcing.jl:6
[3] (::ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing})(::Vector{Float64}, ::Vararg{Any})
@ SciMLBase ~/.julia/packages/SciMLBase/t5VdI/src/scimlfunctions.jl:1613