# Incorporating forcing functions in the ODE model

Hello,

I am new to Turing and Julia, and need some help!

I have been trying to estimate the parameters of a model with a discrete forcing function `q(t)`, where the values of `q(t)` are discrete and are read from a file. Consider the ODE to be following:

``````function my_model(du, u, p, t)
α, β = p
du = α*q(t)*u - β
end
``````

Here, alpha and beta are my parameters, which I tried to estimate by building a model as follows, where I am passing the forcing function data explicitly as a vector:

``````@model function fitting_model(observ_data, forcing_data)
# Some priors...
α ~ truncated(Normal(0.25, 0.1), 0.1, 0.5)
β ~ truncated(Normal(0.65, 0.1), 0, 2)

# Parameter vector
p = [β, α, forcing_data];

# Defining other required variables like u0, tspan, etc...

# Then, the (usual) subsequent steps to estimate parameters
problem = ODEProblem(my_model, u0, tspan, p);
predicted = solve(problem, Tsit5(), saveat=1.0);

# Finally matching the predicted data with observ_data
end
``````

But running this code resulted in `TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 7}`. I am at a dead end, and don’t know how to resolve this issue.

1 Like

Can you add whole code you try to run? If not, check the code for hard coded types. The error suggest you hard coded the type `Float64` somewhere and you got a dual number along the way.

Thank you for the quick reply @tomaklutfu! Please check the complete code below:

``````#=
Section 1: Import required packages
=#

using Turing, Distributions, DifferentialEquations, Interpolations
using MCMCChains, Plots, StatsPlots
using CSV, XLSX, DataFrames
using Random
Random.seed!(18431)

#=
Section 2: Read the data file containing observation data and get the NPI data into arrays
=#

total_weeks = 36;   # Total number of time points
N = 67081000;       # Population

y_time = 1:1:total_weeks;               # Timepoints (weeks)

y_S = Float64.(my_data.Susceptible);    # Susceptible
y_S = y_S[1:total_weeks];

y_D = Float64.(my_data.Deceased);       # Deceased
y_D = y_D[1:total_weeks];

y_HC = Float64.(my_data.Hosp_critical); # Critical hospitalizations
y_HC = y_HC[1:total_weeks];

y_T = Float64.(my_data.Hosp_total);     # Total hospitalizations
y_T = y_T[1:total_weeks];
y_HNC = y_T - y_HC;                     # Non-critical hospitalizations

observation_data = [y_S y_D y_HC y_HNC];

# IPTCC is a forcing function
IPTCC = wet_data.Normalized_IPTCC;
IPTCC = IPTCC[1:total_weeks];

# mobil is another forcing function
mobil = mobil_data.Mean;
mobil = mobil[1:total_weeks];

wet_forcing = interpolate(IPTCC, BSpline(Linear()));
mobil_forcing = interpolate(mobil, BSpline(Linear()));

forcing_params = (wet_forcing, mobil_forcing);

#=
Section 3: Define the model and the respective parameters
=#

function epidemic_wildtype(dy, y, p, t)
S, E, I, Hᵪ, Hₙ, R, D = y;
β, λ, α, γ, θᵪ, θₙ, γᵪ, γₙ, δᵪ, w, m = p;
N = 67081000;

dy[1] = -β*w(t)*m(t)*I*S/N + λ*R;  # S
dy[2] = β*w(t)*m(t)*I*S/N - α*E;   # E
dy[3] = α*E - (γ + θᵪ + θₙ)*I;     # I
dy[4] = θₙ*I - γₙ*Hᵪ;               # HNC
dy[5] = θᵪ*I - (γᵪ + δᵪ)*Hₙ;       # HC
dy[6] = γ*I + γₙ*Hₙ + γᵪ*Hᵪ - λ*R;  # R
dy[7] = δᵪ*Hᵪ;                     # D
end

#=
Section 4: Define the priors and the Bayesian model
=#

@model function fitting_epidemic_wildtype(observ_data, w_forcing, m_forcing)
# Priors of model parameters
β ~ truncated(Normal(0.65, 0.1), 0, 2)
λ ~ truncated(Normal(0.5, 0.1), 0, 5)
α ~ truncated(Normal(0.25, 0.1), 0.1, 0.5)
γ ~ truncated(Normal(0.05, 0.1), 0, 5)
γₙ ~ Uniform(0.05, 0.1)
γᵪ ~ Uniform(0.05, 0.1)
θₙ ~ Uniform(0.09, 0.75)
θᵪ ~ Uniform(0.09, 0.75)
δᵪ ~ Uniform(0.1, 0.8)

p = [β, λ, α, γ, θᵪ, θₙ, γᵪ, γₙ, δᵪ, w_forcing, m_forcing];

# Priors of standard deviations
σ₁ ~ InverseGamma(1, 1) # Susceptible
σ₂ ~ InverseGamma(1, 1) # Deceased
σ₃ ~ InverseGamma(2, 3) # Critically hospitalized
σ₄ ~ InverseGamma(1, 1) # Non-critically hospitalized

# Initial conditions
N = 67081000;
S0 = N;
I0 = 100;
y0 = [S0, 0, I0, 0, 0, 0, 0];
y0 = eltype(p).(y0);

# Solve the model and compare with observed data
problem = ODEProblem(epidemic_wildtype, y0, (1, 36), p)
predicted = solve(problem, Tsit5(), saveat=1.0)

for i = 1:length(predicted)
observ_data[i,1] ~ Normal(predicted[1,i], σ₁)
observ_data[i,2] ~ Normal(predicted[7,i], σ₂)
observ_data[i,3] ~ Normal(predicted[5,i], σ₃)
observ_data[i,4] ~ Normal(predicted[4,i], σ₄)
end
end

#=
Section 5: Run the model-inference system and save the chains
=#

model = fitting_epidemic_wildtype(observation_data, wet_forcing, mobil_forcing);
number_of_chains = 1;
chain = sample(model, NUTS(0.65), MCMCThreads(), 10000, number_of_chains);
``````

Hope it helps to debug! Thanks in advance.

This array might be expecting a `Float64` eltype. The eltype of this array is likely defined to be the same as that of `y0`. I would try to update y0 to

``````y0 = eltype(p).(y0)
``````

this makes `y0` an array of the same element type as the parameter vector `p`, which is probably of type `Dual`.

Thank you @baggepinnen for the quick response. Unfortunately, it looks like the eltype still doesn’t solve it. The error now is `MethodError: no constructors have been defined for Any`. As far as I can see, all the data is of `Float64` type, so don’t really know how/where the `Any` type is coming from!

Any help is really appreciated!

``````DataFrame(XLSX.readtable("observation_data.xlsx","Sheet1")...)
``````

Will mean that all your data will be of type `Any` - try the `infer_eltypes = true` kwarg.

Post the whole error. That stacktrace that you are ignoring is the most important part.

My bad! Just realized that the whole stacktrace is important. Here it is:

``````┌ Warning: Only a single thread available: MCMC chains are not sampled in parallel

Stacktrace:
[1] wait
[3] macro expansion
[4] macro expansion
[5] (::AbstractMCMC.var"#31#41"{Bool, Base.Iterators.Pairs{Symbol, UnionAll, Tuple{Symbol}, NamedTuple{(:chain_type,), Tuple{UnionAll}}}, Int64, Int64, Vector{Any}, Vector{UInt64}, Vector{DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}}, Vector{DynamicPPL.Model{typeof(fitting_epidemic_wildtype), (:observ_data, :w_forcing, :m_forcing), (), (), Tuple{Matrix{Float64}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, Tuple{}, DynamicPPL.DefaultContext}}, Vector{Random._GLOBAL_RNG}})()

nested task error: MethodError: no constructors have been defined for Any
Stacktrace:
[3] getindex
[4] copy
[5] materialize
[6] fitting_epidemic_wildtype(__model__::DynamicPPL.Model{typeof(fitting_epidemic_wildtype), (:observ_data, :w_forcing, :m_forcing), (), (), Tuple{Matrix{Float64}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__::DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, __context__::DynamicPPL.SamplingContext{DynamicPPL.SampleFromUniform, DynamicPPL.DefaultContext, Random._GLOBAL_RNG}, observ_data::Matrix{Float64}, w_forcing::Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, m_forcing::Interpolations.BSplineInterpolation{Float64, 1,
Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}})
@ Main c:\Users\Bharadwaj\Indian Institute of Science\COVID-19 variant study - General\codes_wildtype_fitting\discourse_NPI_code.jl:100
[7] macro expansion
[8] _evaluate
[10] Model
[11] Model
[12] VarInfo
[13] VarInfo
[14] step(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitting_epidemic_wildtype), (:observ_data, :w_forcing, :m_forcing), (), (), Tuple{Matrix{Float64}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, Tuple{}, DynamicPPL.DefaultContext}, spl::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}; resume_from::Nothing, kwargs::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}})
[15] macro expansion
[16] macro expansion
[17] mcmcsample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitting_epidemic_wildtype), (:observ_data, :w_forcing, :m_forcing), (), (), Tuple{Matrix{Float64}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, N::Int64; progress::Bool, progressname::String, callback::Nothing, discard_initial::Int64, thinning::Int64, chain_type::Type, kwargs::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}})
[18] #sample#40
[19] macro expansion
[20] (::AbstractMCMC.var"#927#threadsfor_fun#42"{UnitRange{Int64}, Bool, Base.Iterators.Pairs{Symbol, UnionAll, Tuple{Symbol}, NamedTuple{(:chain_type,), Tuple{UnionAll}}}, Int64, Vector{Any}, Vector{UInt64}, Vector{DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}}, Vector{DynamicPPL.Model{typeof(fitting_epidemic_wildtype), (:observ_data, :w_forcing, :m_forcing), (), (), Tuple{Matrix{Float64}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, Tuple{}, DynamicPPL.DefaultContext}}, Vector{Random._GLOBAL_RNG}})(onethread::Bool)
[21] (::AbstractMCMC.var"#927#threadsfor_fun#42"{UnitRange{Int64}, Bool, Base.Iterators.Pairs{Symbol, UnionAll, Tuple{Symbol}, NamedTuple{(:chain_type,), Tuple{UnionAll}}}, Int64, Vector{Any}, Vector{UInt64}, Vector{DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}}, Vector{DynamicPPL.Model{typeof(fitting_epidemic_wildtype), (:observ_data, :w_forcing, :m_forcing), (), (), Tuple{Matrix{Float64}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, Tuple{}, DynamicPPL.DefaultContext}}, Vector{Random._GLOBAL_RNG}})()
Stacktrace:
[1] sync_end(c::Channel{Any})
[2] macro expansion
[3] macro expansion
[4] macro expansion
[5] macro expansion
[6] mcmcsample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitting_epidemic_wildtype), (:observ_data, :w_forcing, :m_forcing), (), (), Tuple{Matrix{Float64}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, ::MCMCThreads, N::Int64, nchains::Int64; progress::Bool, progressname::String, kwargs::Base.Iterators.Pairs{Symbol, UnionAll, Tuple{Symbol}, NamedTuple{(:chain_type,), Tuple{UnionAll}}})
[7] sample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitting_epidemic_wildtype), (:observ_data, :w_forcing, :m_forcing), (), (), Tuple{Matrix{Float64}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, ensemble::MCMCThreads, N::Int64, n_chains::Int64; chain_type::Type, progress::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
[8] sample
[9] #sample#6
[10] sample
[11] #sample#5
[12] sample(model::DynamicPPL.Model{typeof(fitting_epidemic_wildtype), (:observ_data, :w_forcing, :m_forcing), (), (), Tuple{Matrix{Float64}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, Tuple{}, DynamicPPL.DefaultContext}, alg::NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}, ensemble::MCMCThreads, N::Int64, n_chains::Int64)
[13] top-level scope
@ c:\Users\Bharadwaj\Indian Institute of Science\COVID-19 variant study - General\codes_wildtype_fitting\discourse_NPI_code.jl:120
in expression starting at c:\Users\Bharadwaj\Indian Institute of Science\COVID-19 variant study - General\codes_wildtype_fitting\discourse_NPI_code.jl:120
``````

This is the error I am getting after incorporating the suggestions of @baggepinnen and @nilshg. Please check this and help me out!

which line is line 100 in your code?

@ChrisRackauckas Hello! This was the line 100.

Put an `@show typeof(y0)` and `@show eltype(p)` right before that line and share what they print out.

This is the output from the two lines:

``````┌ Warning: Only a single thread available: MCMC chains are not sampled in parallel
typeof(y0) = Vector{Int64}
eltype(p) = Any
``````

… and the previously mentioned error continues.

`eltype(p) is `Any` because you put functions in there.

``````p = [β, λ, α, γ, θᵪ, θₙ, γᵪ, γₙ, δᵪ, w_forcing, m_forcing];
``````

So instead you want to upconvert to just the scalar dual types. In fact, it would be better to do:

``````p = (β, λ, α, γ, θᵪ, θₙ, γᵪ, γₙ, δᵪ, w_forcing, m_forcing);
``````

so that is type-stable, and then do:

``````y0 = typeof(β).(y0)
``````

so it converts `y0` to a dual.

Tried the modifications you suggested, and this time got different error.

``````nested task error: BoundsError: attempt to access 36-element interpolate(::Vector{Float64}, BSpline(Linear())) with element type Float64 at index [NaN]
Stacktrace:
[1] throw_boundserror(A::Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, I::Tuple{Float64})
@ Base .\abstractarray.jl:651
[2] BSplineInterpolation
@ C:\Users\vembh\.julia\packages\Interpolations\3gTQB\src\b-splines\indexing.jl:6 [inlined]
[3] getindex(itp::Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, i::Float64)
@ Interpolations .\deprecated.jl:72
[4] epidemic_wildtype...
``````

at line 61, where line 61 is `dy[1] = -β*w(t)*m(t)*I*S/N + λ*R;` and line 106 is `predicted = solve(problem, Tsit5(), saveat=1.0)`.

Are you solving past the time point where you interpolate? Before that line add `@show t` and figure out what the last point being interpolated is. That should make easy to then do `w(t)` or `m(t)` and reproduce the error without Turing.

1. I do not think I am solving beyond the last point of interpolation. The image below proves that! I am solving for time interval `(1, 36)` and the interpolated forcing functions have the values at the same time points.

1. I am posting the code here again for convenience and the error (which is completely new this time and does not show `BoundsError`) getting generated. Incorporating only your suggestions @ChrisRackauckas:
``````#=
Section 1: Import required packages
=#

using Turing, Distributions, DifferentialEquations, Interpolations
using MCMCChains, Plots, StatsPlots
using CSV, XLSX, DataFrames
using Random
Random.seed!(18431)

#=
Section 2: Read the data file containing observation data and get the NPI data into arrays
=#

my_data = DataFrame(XLSX.readtable("observation_data.xlsx","Sheet1"; infer_eltypes = true)...);

total_weeks = 36;   # Total number of time points
N = 67081000;       # Population

y_time = 1:1:total_weeks;               # Timepoints (weeks)

y_S = Float64.(my_data.Susceptible);    # Susceptible
y_S = y_S[1:total_weeks];

y_D = Float64.(my_data.Deceased);       # Deceased
y_D = y_D[1:total_weeks];

y_HC = Float64.(my_data.Hosp_critical); # Critical hospitalizations
y_HC = y_HC[1:total_weeks];

y_T = Float64.(my_data.Hosp_total);     # Total hospitalizations
y_T = y_T[1:total_weeks];
y_HNC = y_T - y_HC;                     # Non-critical hospitalizations

observation_data = [y_S y_D y_HC y_HNC];

wet_data = DataFrame(XLSX.readtable("Wetdata.xlsx","Wetdata"; infer_eltypes = true)...);
# IPTCC is a forcing function
IPTCC = wet_data.Normalized_IPTCC;
IPTCC = IPTCC[1:total_weeks];

mobil_data = DataFrame(XLSX.readtable("Mobdata.xlsx","Mobdata"; infer_eltypes = true)...);
# mobil is another forcing function
mobil = mobil_data.Mean;
mobil = mobil[1:total_weeks];

wet_forcing = interpolate(IPTCC, BSpline(Linear()));
mobil_forcing = interpolate(mobil, BSpline(Linear()));

forcing_params = (wet_forcing, mobil_forcing);

#=
Section 3: Define the model and the respective parameters
=#

function epidemic_wildtype(dy, y, p, t)
S, E, I, Hᵪ, Hₙ, R, D = y;
β, λ, α, γ, θᵪ, θₙ, γᵪ, γₙ, δᵪ, w, m = p;
N = 67081000;

dy[1] = -β*w(t)*m(t)*I*S/N + λ*R;  # S
dy[2] = β*w(t)*m(t)*I*S/N - α*E;   # E
dy[3] = α*E - (γ + θᵪ + θₙ)*I;     # I
dy[4] = θₙ*I - γₙ*Hᵪ;               # HNC
dy[5] = θᵪ*I - (γᵪ + δᵪ)*Hₙ;       # HC
dy[6] = γ*I + γₙ*Hₙ + γᵪ*Hᵪ - λ*R;  # R
dy[7] = δᵪ*Hᵪ;                     # D
end

#=
Section 4: Define the priors and the Bayesian model
=#

@model function fitting_epidemic_wildtype(observ_data, w_forcing, m_forcing)
# Priors of model parameters
β ~ truncated(Normal(0.65, 0.1), 0, 2)
λ ~ truncated(Normal(0.5, 0.1), 0, 5)
α ~ truncated(Normal(0.25, 0.1), 0.1, 0.5)
γ ~ truncated(Normal(0.05, 0.1), 0, 5)
γₙ ~ Uniform(0.05, 0.1)
γᵪ ~ Uniform(0.05, 0.1)
θₙ ~ Uniform(0.09, 0.75)
θᵪ ~ Uniform(0.09, 0.75)
δᵪ ~ Uniform(0.1, 0.8)

p = (β, λ, α, γ, θᵪ, θₙ, γᵪ, γₙ, δᵪ, w_forcing, m_forcing);

# Priors of standard deviations
σ₁ ~ InverseGamma(1, 1) # Susceptible
σ₂ ~ InverseGamma(1, 1) # Deceased
σ₃ ~ InverseGamma(2, 3) # Critically hospitalized
σ₄ ~ InverseGamma(1, 1) # Non-critically hospitalized

# Initial conditions
N = 67081000;
S0 = N;
I0 = 100;
y0 = [S0, 0, I0, 0, 0, 0, 0];
@show typeof(y0)
@show eltype(p)
y0 = typeof(β).(y0);

# Solve the model and compare with observed data
problem = ODEProblem(epidemic_wildtype, y0, (1, 36), p)
predicted = solve(problem, Tsit5(), saveat=1)

for i = 1:length(predicted)
observ_data[i,1] ~ Normal(predicted[1,i], σ₁)
observ_data[i,2] ~ Normal(predicted[7,i], σ₂)
observ_data[i,3] ~ Normal(predicted[5,i], σ₃)
observ_data[i,4] ~ Normal(predicted[4,i], σ₄)
end
end

#=
Section 5: Run the model-inference system and save the chains
=#

model = fitting_epidemic_wildtype(observation_data, wet_forcing, mobil_forcing);
number_of_chains = 1;
chain = sample(model, NUTS(0.65), MCMCThreads(), 10000, number_of_chains);
``````

Error without the Turing error part:

``````nested task error: InexactError: Int64(161//1000)
Stacktrace:
[1] Integer
@ .\rational.jl:110 [inlined]
[2] convert(#unused#::Type{Int64}, x::Rational{Int64})
@ Base .\number.jl:7
[3] OrdinaryDiffEq.Tsit5ConstantCache(T::Type, T2::Type)
[4] alg_cache(alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, u::Vector{Float64}, rate_prototype::Vector{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Int64}, uprev::Vector{Float64}, uprev2::Vector{Float64}, f::ODEFunction{true, typeof(epidemic_wildtype), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing,
Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, t::Int64, dt::Int64, reltol::Float64, p::Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, calck::Bool, #unused#::Val{true})
[5] __init(prob::ODEProblem{Vector{Float64}, Tuple{Int64, Int64}, true, Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, ODEFunction{true, typeof(epidemic_wildtype), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Int64, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Int64, dtmin::Nothing, dtmax::Int64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
[6] #__solve#471
[7] #solve_call#42
[8] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Int64, Int64}, true, Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, ODEFunction{true, typeof(epidemic_wildtype), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}; kwargs::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Int64}}})
[9] #solve#43
[10] fitting_epidemic_wildtype(__model__::DynamicPPL.Model{typeof(fitting_epidemic_wildtype), (:observ_data, :w_forcing, :m_forcing), (), (), Tuple{Matrix{Float64}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__::DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, __context__::DynamicPPL.SamplingContext{DynamicPPL.SampleFromUniform, DynamicPPL.DefaultContext, Random._GLOBAL_RNG}, observ_data::Matrix{Float64}, w_forcing::Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, m_forcing::Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}})
@ Main c:\Users\Bharadwaj\Indian Institute of Science\COVID-19 variant study - General\codes_wildtype_fitting\discourse_NPI_code.jl:106
[11] macro expansion
[12] _evaluate
[14] Model
[15] Model
[16] VarInfo
[17] VarInfo
[18] step(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitting_epidemic_wildtype), (:observ_data, :w_forcing, :m_forcing), (), (), Tuple{Matrix{Float64}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}}, Tuple{}, DynamicPPL.DefaultContext}, spl::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}; resume_from::Nothing, kwargs::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}})
``````

where line 106 is `predicted = solve(problem, Tsit5(), saveat=1)`. Hope this helps; sorry for the mess!

Make your time a floating point type: `(1.0,36.0)`. That should be caught and fixed automatically, but let’s see if that’s all that’s left now.

This brought back the `BoundsError` error that I posted previously! Posting the part of the (extremely long) error message stacktrace.

``````nested task error: BoundsError: attempt to access 36-element interpolate(::Vector{Float64}, BSpline(Linear())) with element type Float64 at index [NaN]
Stacktrace:
[1] throw_boundserror(A::Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, I::Tuple{Float64})
@ Base .\abstractarray.jl:651
[2] BSplineInterpolation
@ Main C:\Users\Bharadwaj\Indian Institute of Science\COVID-19 variant study - General\codes_wildtype_fitting\discourse_NPI_code.jl:61
``````

Here, line 61 is the first differential equation.

Hope this helps you @ChrisRackauckas in helping me out…

By the way, when I am just solving the ODE with randomly chosen parameter values (and not using Turing) along with the forcing functions, the code runs error-free! Here is the part of the code without the Turing part:

``````#=
Section 1: Import required packages =#

using Turing, Distributions, DifferentialEquations, Interpolations
using MCMCChains, Plots, StatsPlots
using CSV, XLSX, DataFrames
using Random
Random.seed!(18431)

#=
Section 2: Read the data file containing observation data and get the NPI data into arrays =#

total_weeks = 36;   # Total number of time points
N = 67081000;       # Population

y_time = 1:1:total_weeks;               # Timepoints (weeks)

y_S = Float64.(my_data.Susceptible);    # Susceptible
y_S = y_S[1:total_weeks];

y_D = Float64.(my_data.Deceased);       # Deceased
y_D = y_D[1:total_weeks];

y_HC = Float64.(my_data.Hosp_critical); # Critical hospitalizations
y_HC = y_HC[1:total_weeks];

y_T = Float64.(my_data.Hosp_total);     # Total hospitalizations
y_T = y_T[1:total_weeks];
y_HNC = y_T - y_HC;                     # Non-critical hospitalizations

observation_data = [y_S y_D y_HC y_HNC];

# IPTCC is a forcing function
IPTCC = wet_data.Normalized_IPTCC;
IPTCC = IPTCC[1:total_weeks];

# mobil is another forcing function
mobil = mobil_data.Mean;
mobil = mobil[1:total_weeks];

wet_forcing = interpolate(IPTCC, BSpline(Linear()));
mobil_forcing = interpolate(mobil, BSpline(Linear()));

#=
Section 3: Define the model and the respective parameters =#

function epidemic_wildtype!(dy, y, p, t)
S, E, I, Hᵪ, Hₙ, R, D = y;
β, λ, α, γ, θᵪ, θₙ, γᵪ, γₙ, δᵪ, w, m = p;
N = 67081000;

dy[1] = -β * w(t) * m(t) * I * S / N + λ * R;  # S
dy[2] = β * w(t) * m(t) * I * S / N - α * E;   # E
dy[3] = α * E - (γ + θᵪ + θₙ) * I;     # I
dy[4] = θₙ * I - γₙ * Hᵪ;               # HNC
dy[5] = θᵪ * I - (γᵪ + δᵪ) * Hₙ;       # HC
dy[6] = γ * I + γₙ * Hₙ + γᵪ * Hᵪ - λ * R;  # R
dy[7] = δᵪ * Hᵪ;                     # D
end

#=
Section 4: Solve it! =#

p = (0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, wet_forcing, mobil_forcing);
tspan = (1.0, 36.0);
N = 67081000;
S0 = N;
I0 = 100;
y0 = [S0, 0, I0, 0, 0, 0, 0];
problem1 = ODEProblem(epidemic_wildtype!, y0, tspan, p);
solution1 = solve(problem1, Tsit5(), saveat=1.0);
``````

Also, even when I am defining the forcing functions as global variables rather than passing them as a part of the parameter vector `p`, the ODE system is being solved error-free. But the same `BoundsError` (mentioned in the previous post) is getting triggered.

Can you do `@show t` before that line and see what `w(t)` and `m(t)` give for the last `t` shown before the error? Most likely, those two throw somewhere.