BoundsError on Bayesian ODE model with Turing

I posted this question here already:

but unfortunately it has not attracted any responses, so I think it may be more appropriate in this category. (I did not seem to be able to edit the category in the original post.)

I am trying to infer the parameters (a11, a12, …, c21, c22) of the following model:

``````a11 ~ Normal(pE.A[1,1], pC.A[1,1])
a12 ~ Normal(pE.A[1,2], pC.A[1,2])
a21 ~ Normal(pE.A[2,1], pC.A[2,1])
a22 ~ Normal(pE.A[2,2], pC.A[2,2])

c11 ~ Normal(pE.C[1,1], pC.C[1,1])
c12 ~ Normal(pE.C[1,2], pC.C[1,2])
c21 ~ Normal(pE.C[2,1], pC.C[2,1])
c22 ~ Normal(pE.C[2,2], pC.C[2,2])

x0 = zeros(options.ℓ, 1) # initial condition
z0 = zeros(options.ℓ, 6) # initial condition
states = [x0, z0]
P = [a11, a12, a21, a22, c11, c12, c21, c22, u, H, HC]
tspan = (0.0, options.duration)
predicted = forward_model(P, states, options)

y[:,1] ~ Normal(predicted[:,1])
y[:,2] ~ Normal(predicted[:,2])

end

m = invert(U, y, pE, pC, options)
chain = mapreduce(c -> sample(m, NUTS(.65),1000), chainscat, 1:3)
``````

but get the error message:

``````BoundsError: attempt to access 1-element Vector{Matrix{Float64}} at index [2]
``````

The error occurs in a function called g, which is called by the function forward_model:

``````function forward_model(params, states, options)
a11, a12, a21, a22, c11, c12, c21, c22, u, H, HC = params
x0, z0 = states

ts = range(0.05, step=options.dt, stop=options.duration) # sample times
tspan = (0.0, options.duration) # time span

p = [a11, a12, a21, a22, c11, c12, c21, c22, u, options.dt] # parameters
ode = ODEProblem(f, x0, tspan, p)
sol = solve(ode, Tsit5(), saveat=ts, tstops=ts, dt=1e-2)

ph = [H, HC, sol.u, options] # parameters
odeh = ODEProblem(h, z0, tspan, ph)
solh = solve(odeh, Tsit5(), saveat=ts, tstops=ts, dt = 1e-2)

# Model signal
y = g(solh.u, H, options)
end

function g(z, H, opt)
# Define constants
TE = 0.04
V₀ = 100*0.04
r₀  = 110
ϑ₀ = 84.8
E₀ = opt.Pf_H[:,5] .* exp.(H[:,5])
εₕ = opt.Pf_H[2,6] .* exp.(H[2,6])

# Compute coefficients
k₁ = 4.3.*ϑ₀.*E₀.*TE
k₂ = εₕ.*r₀.*E₀.*TE
k₃ = 1 - εₕ

# Output equation
v = map(i -> (z[i][:,3]), 1:12000)
v = transpose(reduce(hcat, v))
q = map(i -> (z[i][:,4]), 1:12000)
q = transpose(reduce(hcat, q))
y = V₀ * (k₁'.*(1 .- q) + k₂'.*(1 .- q./v) + k₃'.*(1 .- v))
end
``````

For the sake of simplicity I have left out all other functions (f, h, etc.); I think that these should not be necessary for answering my question. The input z into function g is the solution of an ODE problem and has the following form:

``````12000-element Vector{Matrix{Float64}}:
[0.0 1.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0]
[0.0 1.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0]
⋮
``````

Generating data with the function forward_model works. The BoundsError during the inference procedure occurs in the line:

``````v = map(i -> (z[i][:,3]), 1:12000)
``````

(in function g). I thought that the error occurs because the input argument y (in function invert) is a Matrix{Float64} of size (12000, 2), while z (in function g) is a Vector{Matrix{Float64}} of size (12000,). However, after converting the data to a Vector{Matrix{Float64}} of size (12000,) as follows:

``````d = [hcat(c[:]) for c in eachrow(data)]
d[:,1] ~ Normal(predicted[:,1])
d[:,2] ~ Normal(predicted[:,2])
``````

the BoundsError still occurrs in the same place.

I don’t understand what exactly is causing the error and if there is something wrong with my implementation of the Turing model and/or my combination of Turing + DifferentialEquations or if this is some more fundamental issue. I would be very grateful for any possible explanations or suggestions of where I am going wrong. I am new to Julia and don’t have much experience with probabilistic programming and differential equations either.

Could you maybe try to boil it down to a minimal example that reproduces the error? Also full stacktraces would be helpful (you can make them collapsible).

It is a lot easier for people to help, if they can run the code and figure out themselves where the error is coming from (see some other good suggestions here: Please read: make it easier to help you)

1 Like

This. I ignored it the first time because there wasn’t any code I could run. Now it’s posted again without runnable code. Make a simple example someone can copy-paste to see your error (simple as in, delete the unnecessary parts. You’ll be surprised how that help you debug it), and then you can expect people to join in. Otherwise it’s people reading a post with incomplete information (you didn’t post the stacktraces I would read to debug it).

Thank you both for your replies and I am sorry for not providing enough information. I am working on the simplified runnable code and will edit my post once it’s ready. In the meantime, here is the stacktrace:

``````┌ Warning: Instability detected. Aborting
└ @ SciMLBase /home/birtet/.julia/packages/SciMLBase/n3U0M/src/integrator_interface.jl:351

BoundsError: attempt to access 1-element Vector{Matrix{Float64}} at index [2]

Stacktrace:
[1] getindex
@ ./array.jl:801 [inlined]
[2] #101
@ ./In[63]:17 [inlined]
[3] iterate
@ ./generator.jl:47 [inlined]
[4] collect_to!(dest::Vector{Vector{Float64}}, itr::Base.Generator{UnitRange{Int64}, var"#101#103"{Vector{Matrix{Float64}}}}, offs::Int64, st::Int64)
@ Base ./array.jl:724
[5] collect_to_with_first!
@ ./array.jl:702 [inlined]
[6] _collect(c::UnitRange{Int64}, itr::Base.Generator{UnitRange{Int64}, var"#101#103"{Vector{Matrix{Float64}}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
@ Base ./array.jl:696
[7] collect_similar
@ ./array.jl:606 [inlined]
[8] map
@ ./abstractarray.jl:2294 [inlined]
[9] ldcm_g(z::Vector{Matrix{Float64}}, H::Matrix{Float64}, opt::opt)
@ Main ./In[63]:17
[10] ldcm_forward_model(params::Vector{Any}, states::Vector{Matrix{Float64}}, options::opt)
@ Main ./In[53]:27
[11] #355
@ ./In[202]:34 [inlined]
[12] (::var"#355#359")(__model__::DynamicPPL.Model{var"#355#359", (:u, :data, :pE, :pC, :options), (), (), Tuple{Matrix{Float64}, Matrix{Float64}, priorE, priorC, opt}, 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}, u::Matrix{Float64}, data::Matrix{Float64}, pE::priorE, pC::priorC, options::opt)
@ Main ./none:0
[13] macro expansion
@ ~/.julia/packages/DynamicPPL/yvQQ0/src/model.jl:465 [inlined]
[14] _evaluate
@ ~/.julia/packages/DynamicPPL/yvQQ0/src/model.jl:448 [inlined]
@ ~/.julia/packages/DynamicPPL/yvQQ0/src/model.jl:421 [inlined]
[16] Model
@ ~/.julia/packages/DynamicPPL/yvQQ0/src/model.jl:389 [inlined]
[17] Model
@ ~/.julia/packages/DynamicPPL/yvQQ0/src/model.jl:383 [inlined]
[18] VarInfo
@ ~/.julia/packages/DynamicPPL/yvQQ0/src/varinfo.jl:132 [inlined]
[19] VarInfo
@ ~/.julia/packages/DynamicPPL/yvQQ0/src/varinfo.jl:131 [inlined]
[20] step(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#355#359", (:u, :data, :pE, :pC, :options), (), (), Tuple{Matrix{Float64}, Matrix{Float64}, priorE, priorC, opt}, 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}}})
@ DynamicPPL ~/.julia/packages/DynamicPPL/yvQQ0/src/sampler.jl:69
[21] macro expansion
@ ~/.julia/packages/AbstractMCMC/BPJCW/src/sample.jl:123 [inlined]
[22] macro expansion
@ ~/.julia/packages/ProgressLogging/6KXlp/src/ProgressLogging.jl:328 [inlined]
[23] (::AbstractMCMC.var"#21#22"{Bool, String, Nothing, Int64, Int64, Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}}, Random._GLOBAL_RNG, DynamicPPL.Model{var"#355#359", (:u, :data, :pE, :pC, :options), (), (), Tuple{Matrix{Float64}, Matrix{Float64}, priorE, priorC, opt}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, Int64, Int64})()
@ AbstractMCMC ~/.julia/packages/AbstractMCMC/BPJCW/src/logging.jl:11
[24] with_logstate(f::Function, logstate::Any)
@ Base.CoreLogging ./logging.jl:491
[25] with_logger(f::Function, logger::LoggingExtras.TeeLogger{Tuple{LoggingExtras.EarlyFilteredLogger{ConsoleProgressMonitor.ProgressLogger, AbstractMCMC.var"#1#3"{Module}}, LoggingExtras.EarlyFilteredLogger{Base.CoreLogging.SimpleLogger, AbstractMCMC.var"#2#4"{Module}}}})
@ Base.CoreLogging ./logging.jl:603
[26] with_progresslogger(f::Function, _module::Module, logger::Base.CoreLogging.SimpleLogger)
@ AbstractMCMC ~/.julia/packages/AbstractMCMC/BPJCW/src/logging.jl:34
[27] macro expansion
@ ~/.julia/packages/AbstractMCMC/BPJCW/src/logging.jl:10 [inlined]
[28] mcmcsample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#355#359", (:u, :data, :pE, :pC, :options), (), (), Tuple{Matrix{Float64}, Matrix{Float64}, priorE, priorC, opt}, 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}}})
@ AbstractMCMC ~/.julia/packages/AbstractMCMC/BPJCW/src/sample.jl:114
[29] sample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#355#359", (:u, :data, :pE, :pC, :options), (), (), Tuple{Matrix{Float64}, Matrix{Float64}, priorE, priorC, opt}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, N::Int64; chain_type::Type, resume_from::Nothing, progress::Bool, nadapts::Int64, discard_adapt::Bool, discard_initial::Int64, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Turing.Inference ~/.julia/packages/Turing/y0DW3/src/inference/hmc.jl:133
[30] sample
@ ~/.julia/packages/Turing/y0DW3/src/inference/hmc.jl:116 [inlined]
[31] #sample#2
@ ~/.julia/packages/Turing/y0DW3/src/inference/Inference.jl:142 [inlined]
[32] sample
@ ~/.julia/packages/Turing/y0DW3/src/inference/Inference.jl:142 [inlined]
[33] #sample#1
@ ~/.julia/packages/Turing/y0DW3/src/inference/Inference.jl:132 [inlined]
[34] sample(model::DynamicPPL.Model{var"#355#359", (:u, :data, :pE, :pC, :options), (), (), Tuple{Matrix{Float64}, Matrix{Float64}, priorE, priorC, opt}, Tuple{}, DynamicPPL.DefaultContext}, alg::NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}, N::Int64)
@ Turing.Inference ~/.julia/packages/Turing/y0DW3/src/inference/Inference.jl:132
[35] (::var"#363#364")(c::Int64)
@ Main ./In[203]:7
[36] _mapreduce(f::var"#363#364", op::typeof(chainscat), #unused#::IndexLinear, A::UnitRange{Int64})
@ Base ./reduce.jl:408
[37] _mapreduce_dim(f::Function, op::Function, #unused#::Base._InitialValue, A::UnitRange{Int64}, #unused#::Colon)
@ Base ./reducedim.jl:318
[38] #mapreduce#672
@ ./reducedim.jl:310 [inlined]
[39] mapreduce(f::Function, op::Function, A::UnitRange{Int64})
@ Base ./reducedim.jl:310
[40] top-level scope
@ In[203]:7
[41] eval
@ ./boot.jl:360 [inlined]
[42] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
``````

Well the warning says that it was unstable so it exited, which means it won’t save every point in `ts` and you need to check `if sol.retcode == :Success` like the tutorials show. I don’t see that done here, so I would presume that’s where your issue comes from.

Why is it unstable? See PSA: How to help yourself debug differential equation solving issues . If it’s unstable at the first set of parameters, then it’s likely there’s an issue in the definition of the ODE and you should re-check that. If it’s only after the optimization has started, then maybe the parameters turn stiff while optimizing, in which case you need to use a different solver, or there are just some parameters where the ODE is unstable, which means you need to check and reject.

Thank you for your answer and apologies for my late reply. I had to step away from Julia for a while and needed to get up to speed again.

I have tried implementing your suggestions and ironed out some issues in my code, but I still end up with a BoundsError. Here is a MWE (for simplicity, I have tried just estimating parameter a11; eventually, I would like to estimate multiple parameters):

``````using Parameters, LinearAlgebra, DifferentialEquations, Turing, MCMCChains, Distributions, Random, Plots, StatsPlots

"""Specify options"""

@with_kw mutable struct options
TR::Float64
TE::Float64
duration::Float64
dt::Float64
layers::Int
SNR::Float64
end

"""Define constants"""

@consts begin
λ = 0.5
ϑ₀ = 28.265*3
r₀ = 110
V₀ = 100*0.04

κ = 0.65
γ = 0.41
τ = 0.98
α = 0.32
E₀ = 0.34
ε = 0.4584
τᵈ = 0.5
pₕ = [κ γ τ α E₀ ε τᵈ]
end

"""Functions for the forward model"""

function f!(x, p, t)

a11, a12, a21, a22, c11, c12, c21, c22, u, dt = p
A = [a11 a21; a12 a22]
C = [c11 c21; c12 c22]

# Get input at time t
idx = max(Int(ceil(t/dt)), 1)
ut = u[idx,:]

J = A
J*x + C*ut
end

function h!(dz, z, p, t)

H, HC, xn, dt = p

idx = max(Int(ceil(t/dt)), 1)

H = pₕ .* exp.(H)
HC = [0 1; 0 0] .* exp.(HC)

fv = (abs.(z[:,3])).^(1 ./ H[:,4])
ff = ( 1 .- (1 .- H[:,5]).^(1 ./ z[:,2]) ) ./ H[:,5]

# Compute state equations
dz[:,1] = xn[idx] .- H[:,1].*z[:,1] - H[:,2].*(z[:,2] .- 1)
dz[:,2] = z[:,1]
dz[:,3] = (z[:,2] - fv + HC*z[:,5]) ./ H[:,3]
dz[:,4] = (ff.*z[:,2] - fv.*z[:,4]./z[:,3] + HC*z[:,6]) ./ H[:,3]
dz[:,5] = (-z[:,5] + z[:,3] .- 1)./(H[:,7])
dz[:,6] = (-z[:,6] + z[:,4] .- 1)./(H[:,7])
end

function g(z, H, TE)

E0 = E₀ .* exp.(H[:,5])
εₕ = ε .* exp.(H[2,6])

k₁ = 4.3.*ϑ₀.*E₀.*TE
k₂ = εₕ.*r₀.*E₀.*TE
k₃ = 1 - εₕ

v = map(i -> (z[i][:,3]), 1:12000)
v = transpose(reduce(hcat, v))
q = map(i -> (z[i][:,4]), 1:12000)
q = transpose(reduce(hcat, q))
y = V₀ * (k₁'.*(1 .- q) + k₂'.*(1 .- q./v) + k₃'.*(1 .- v))
end

# Forward model
function forward_model(params, states, opt)
a11, a12, a21, a22, c11, c12, c21, c22, u, H, HC = params
x0, z0 = states

ts = range(opt.dt, step=opt.dt, stop=opt.duration) # sample times
tspan = (0.0, opt.duration)

p = [a11, a12, a21, a22, c11, c12, c21, c22, u, opt.dt] # parameters
ode = ODEProblem(f!, x0, tspan, p)
sol = solve(ode, AutoTsit5(Rosenbrock23()), saveat=ts, tstops=ts, dt=1e-2)
if sol.retcode != :Success
println(sol.retcode)
end

ph = [H, HC, sol.u, opt.dt] # parameters
odeh = ODEProblem(h!, z0, tspan, ph)
solh = solve(odeh, AutoTsit5(Rosenbrock23()), saveat=ts, tstops=ts, dt = 1e-2)
if solh.retcode != :Success
println(sol.retcode)
end

y = g(solh.u, H, opt.TE)
end

"""Simulate data"""
function generate_data(params, states, opt, seed)

y = forward_model(params, states, opt)

# Get times that the response is sampled
v = Int(opt.duration/opt.TR)
ty = ceil.((0:1:(v - 1)) * size(y,1)/v)

yy = zeros(length(ty), size(y,2))
delays = [0.2 0.2]/opt.dt
for i in 1:length(delays)
yy[:,i] = y[Int.(ty .+ delays[i]), i]
end

r = Diagonal(vec((mapslices(std, y; dims=1) ./ opt.SNR)))
a = 1/16
a = [1 -a]
K = inv(Bidiagonal(ones(v)*a[1], ones(v-1)*a[2], :L))
K = K*sqrt(v/tr(K*K'))
z = randn(MersenneTwister(seed), Float64, (v, opt.layers))
ϵ = K*z
Y = yy + ϵ*r
end

"""Functions for model inversion"""

# Model inversion

@model function invert_model(data, params, opt)

# Unpack parameters
a12, a21, a22, c11, c12, c21, c22, u, H, HC = params

# Specify priors
a11 ~ Normal(-0.5, 8)

# Set initial conditions
x0 = zeros(opt.layers, 1) # initial condition
z0 = zeros(opt.layers, 6) # initial condition
z0[:,2:4] = exp.(z0[:,2:4]) # exponentiate haemodynamic state variables
states = [x0, z0]

P = [a11, a12, a21, a22, c11, c12, c21, c22, u, H, HC]
predicted = forward_model(P, states, opt)

σ ~ LogNormal(6, 1/128) # Check what this should actually be

for i = 1:length(predicted)
data[:,i] ~ MvNormal(predicted[i], σ)
end
end

## 1. Test forward model
# Initialise options:
opt = options(TR = 2, TE = 0.03, duration = 600.0, dt = 0.05, layers = 2, SNR = 4.0) # keywords optional
# Generate a single impulse as input:
u = zeros(Int(opt.duration / opt.dt), 2)
u[6000:6004] .= 1
# Set parameters:
H = zeros(2,7)
HC = log(λ) * [0 1; 0 0]
parameters = [-1.1, 0, -1.1, 0, 0.3, 0, 0.3, 0, u, H, HC]
# Initial states:
x0 = zeros(opt.layers, 1)
z0 = zeros(opt.layers, 6)
z0[:,2:4] = exp.(z0[:,2:4])
states = [x0, z0]
# Run model:
y = forward_model(parameters, states, opt);
plot(y)

## 2. Simulate data
Y = generate_data(parameters, states, opt, 1234)
plot(Y)

## 3. Test model inversion
# Set parameters:
p_set = [0, -1.1, 0, 0.3, 0, 0.3, 0, u, H, HC] # only estimate a11 for now
model = invert_model(Y, p_set, opt)
chain = mapreduce(c -> sample(model, NUTS(.65), 1000), chainscat, 1:3)
# Plot
plot(chain)
``````

My code produces the following error:

``````┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/DXiE6/src/integrator_interface.jl:351
Unstable
┌ Warning: `vendor()` is deprecated, use `BLAS.get_config()` and inspect the output instead
│   caller = (::DefaultLinSolve)(x::Vector{Float64}, A::Matrix{Float64}, b::Vector{Float64}, update_matrix::Bool; reltol::Float64, kwargs::Base.Pairs{Symbol, DiffEqBase.ScaleVector{Matrix{Float64}}, Tuple{Symbol, Symbol}, NamedTuple{(:Pl, :Pr), Tuple{DiffEqBase.ScaleVector{Matrix{Float64}}, DiffEqBase.ScaleVector{Matrix{Float64}}}}}) at linear_nonlinear.jl:91
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/Nyo4y/src/linear_nonlinear.jl:91
ERROR: BoundsError: attempt to access 7699-element Vector{Matrix{Float64}} at index [7700]
Stacktrace:
[1] getindex(A::Vector{Matrix{Float64}}, i1::Int64)
@ Base ./array.jl:861
[2] h!(dz::Matrix{Float64}, z::Matrix{Float64}, p::Vector{Any}, t::Float64)
@ Main ~/Desktop/projects/lDCM/src/test.jl:62
[3] ODEFunction
@ ~/.julia/packages/SciMLBase/DXiE6/src/scimlfunctions.jl:334 [inlined]
[4] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, true, Matrix{Float64}, Nothing, Float64, Vector{Any}, Float64, Float64, Float64, Vector{Matrix{Float64}}, OrdinaryDiffEq.ODECompositeSolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Any}, ODEFunction{true, typeof(h!), 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}, CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.Rosenbrock23Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, DefaultLinSolve, SparseDiffTools.ForwardColorJacCache{Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Matrix{Float64}, Vector{Vector{NTuple{12, Float64}}}, UnitRange{Int64}, Nothing}, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, Float64}, Float64, 1}}}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.Rosenbrock23Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, DefaultLinSolve, SparseDiffTools.ForwardColorJacCache{Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Matrix{Float64}, Vector{Vector{NTuple{12, Float64}}}, UnitRange{Int64}, Nothing}, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, Float64}, Float64, 1}}}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}}, Matrix{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Rosenbrock23Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, DefaultLinSolve, SparseDiffTools.ForwardColorJacCache{Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Matrix{Float64}, Vector{Vector{NTuple{12, Float64}}}, UnitRange{Int64}, Nothing}, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, Float64}, Float64, 1}}}, repeat_step::Bool)
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/perform_step/rosenbrock_perform_step.jl:53
[5] perform_step!
@ ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/perform_step/composite_perform_step.jl:52 [inlined]
[6] perform_step!
@ ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/perform_step/composite_perform_step.jl:49 [inlined]
[7] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, true, Matrix{Float64}, Nothing, Float64, Vector{Any}, Float64, Float64, Float64, Vector{Matrix{Float64}}, OrdinaryDiffEq.ODECompositeSolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Any}, ODEFunction{true, typeof(h!), 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}, CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.Rosenbrock23Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, DefaultLinSolve, SparseDiffTools.ForwardColorJacCache{Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Matrix{Float64}, Vector{Vector{NTuple{12, Float64}}}, UnitRange{Int64}, Nothing}, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, Float64}, Float64, 1}}}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.Rosenbrock23Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, DefaultLinSolve, SparseDiffTools.ForwardColorJacCache{Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Matrix{Float64}, Vector{Vector{NTuple{12, Float64}}}, UnitRange{Int64}, Nothing}, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(h!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, Float64}, Float64, 1}}}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}}, Matrix{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/solve.jl:477
[8] __solve(::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Any}, ODEFunction{true, typeof(h!), 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}, ::CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:saveat, :tstops, :dt), Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64}}})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/solve.jl:5
[9] #solve_call#56
@ ~/.julia/packages/DiffEqBase/Nyo4y/src/solve.jl:61 [inlined]
[10] #solve_up#58
@ ~/.julia/packages/DiffEqBase/Nyo4y/src/solve.jl:82 [inlined]
[11] #solve#57
@ ~/.julia/packages/DiffEqBase/Nyo4y/src/solve.jl:70 [inlined]
[12] forward_model(params::Vector{Any}, states::Vector{Matrix{Float64}}, opt::options)
@ Main ~/Desktop/projects/lDCM/src/test.jl:103
[13] #12
@ ~/Desktop/projects/lDCM/src/test.jl:159 [inlined]
[14] (::var"#12#13")(__rng__::Random._GLOBAL_RNG, __model__::DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, __varinfo__::DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, __sampler__::DynamicPPL.SampleFromUniform, __context__::DynamicPPL.DefaultContext, data::Matrix{Float64}, params::Vector{Any}, opt::options)
@ Main ./none:0
[15] macro expansion
@ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:0 [inlined]
[16] _evaluate
@ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:154 [inlined]
@ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:127 [inlined]
[18] Model
@ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:92 [inlined]
[19] DynamicPPL.VarInfo(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, sampler::DynamicPPL.SampleFromUniform, context::DynamicPPL.DefaultContext)
@ DynamicPPL ~/.julia/packages/DynamicPPL/SgzCy/src/varinfo.jl:126
[20] VarInfo
@ ~/.julia/packages/DynamicPPL/SgzCy/src/varinfo.jl:125 [inlined]
[21] step(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, spl::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}; resume_from::Nothing, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}})
@ DynamicPPL ~/.julia/packages/DynamicPPL/SgzCy/src/sampler.jl:73
[22] macro expansion
@ ~/.julia/packages/AbstractMCMC/oou1a/src/sample.jl:97 [inlined]
[23] macro expansion
@ ~/.julia/packages/ProgressLogging/6KXlp/src/ProgressLogging.jl:328 [inlined]
[24] macro expansion
@ ~/.julia/packages/AbstractMCMC/oou1a/src/logging.jl:8 [inlined]
[25] mcmcsample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, 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.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}})
@ AbstractMCMC ~/.julia/packages/AbstractMCMC/oou1a/src/sample.jl:88
[26] sample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, sampler::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, N::Int64; chain_type::Type, resume_from::Nothing, progress::Bool, nadapts::Int64, discard_adapt::Bool, discard_initial::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Turing.Inference ~/.julia/packages/Turing/uAz5c/src/inference/hmc.jl:140
[27] sample
@ ~/.julia/packages/Turing/uAz5c/src/inference/hmc.jl:123 [inlined]
[28] #sample#2
@ ~/.julia/packages/Turing/uAz5c/src/inference/Inference.jl:142 [inlined]
[29] sample
@ ~/.julia/packages/Turing/uAz5c/src/inference/Inference.jl:142 [inlined]
[30] #sample#1
@ ~/.julia/packages/Turing/uAz5c/src/inference/Inference.jl:132 [inlined]
[31] sample(model::DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, alg::NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}, N::Int64)
@ Turing.Inference ~/.julia/packages/Turing/uAz5c/src/inference/Inference.jl:132
[32] (::var"#14#15")(c::Int64)
@ Main ~/Desktop/projects/lDCM/src/test.jl:196
[33] _mapreduce(f::var"#14#15", op::typeof(chainscat), #unused#::IndexLinear, A::UnitRange{Int64})
@ Base ./reduce.jl:410
[34] _mapreduce_dim(f::Function, op::Function, #unused#::Base._InitialValue, A::UnitRange{Int64}, #unused#::Colon)
@ Base ./reducedim.jl:330
[35] #mapreduce#725
@ ./reducedim.jl:322 [inlined]
[36] mapreduce(f::Function, op::Function, A::UnitRange{Int64})
@ Base ./reducedim.jl:322
[37] top-level scope
@ ~/Desktop/projects/lDCM/src/test.jl:196
``````

So it looks to me like there is a problem with the definition of my ODE function h!. What confuses me is that this function seems to perform as expected on its own/ when I generate data with it. It’s only within the Turing model that it doesn’t work. This makes me wonder if h! is the problem or if it’s the way that I’ve set up the Turing model. I’d be very happy about any clues.

As MCMC hops around, it’s going to change the parameters. At which parameters is it unstable?

Thanks for your very quick input.
Do you have a suggestion for making inspection of the parameters easier? I have tried using the Debugger in VS Code, but it is much too slow to be practically useful with this code.

Just print them, `@show prob.p` before the solve. Then look at the ones that occur right before the failure.

Thank you again.

I don’t understand what is going on… the error always occurs on the 3rd iteration, regardless of what value the parameter a11 takes on (it is the only parameter being estimated at the moment). For example:

``````ode.p[1] = -1.1
ode.p[1] = -1.1
ode.p[1] = -1.9917124835785258
``````

In the first 2 iterations the value is somehow always -1.1, the value used to generate the data. This seems curious in itself. The instability occurs once it changes from that value. I have tried using a tighter prior, a11 ~ Normal(0.5, 2), but that didn’t help.

There’s also another issue with my model that I don’t understand. Most of the time I get the BoundsError, but sometimes I get the following instead:

``````ERROR: MethodError: no method matching MvNormal(::Float64, ::Float64)
Closest candidates are:
MvNormal(::Tracker.TrackedVector{<:Real}, ::Real) at ~/.julia/packages/DistributionsAD/WaBSG/src/tracker.jl:468
MvNormal(::AbstractVector{<:Real}, ::Real) at ~/.julia/packages/Distributions/Xrm9e/src/multivariate/mvnormal.jl:220
MvNormal(::Int64, ::Real) at ~/.julia/packages/Distributions/Xrm9e/src/multivariate/mvnormal.jl:226
Stacktrace:
[1] #12
@ ~/Desktop/projects/lDCM/src/test.jl:165 [inlined]
[2] (::var"#12#13")(__rng__::Random._GLOBAL_RNG, __model__::DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, __varinfo__::DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, __sampler__::DynamicPPL.SampleFromUniform, __context__::DynamicPPL.DefaultContext, data::Matrix{Float64}, params::Vector{Any}, opt::options)
@ Main ./none:0
[3] macro expansion
@ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:0 [inlined]
[4] _evaluate
@ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:154 [inlined]
@ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:127 [inlined]
[6] Model
@ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:92 [inlined]
[7] DynamicPPL.VarInfo(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, sampler::DynamicPPL.SampleFromUniform, context::DynamicPPL.DefaultContext)
@ DynamicPPL ~/.julia/packages/DynamicPPL/SgzCy/src/varinfo.jl:126
[8] VarInfo
@ ~/.julia/packages/DynamicPPL/SgzCy/src/varinfo.jl:125 [inlined]
[9] step(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, spl::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}; resume_from::Nothing, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}})
@ DynamicPPL ~/.julia/packages/DynamicPPL/SgzCy/src/sampler.jl:73
[10] macro expansion
@ ~/.julia/packages/AbstractMCMC/oou1a/src/sample.jl:97 [inlined]
[11] macro expansion
@ ~/.julia/packages/ProgressLogging/6KXlp/src/ProgressLogging.jl:328 [inlined]
[12] macro expansion
@ ~/.julia/packages/AbstractMCMC/oou1a/src/logging.jl:8 [inlined]
[13] mcmcsample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, 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.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}})
@ AbstractMCMC ~/.julia/packages/AbstractMCMC/oou1a/src/sample.jl:88
[14] sample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, sampler::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, N::Int64; chain_type::Type, resume_from::Nothing, progress::Bool, nadapts::Int64, discard_adapt::Bool, discard_initial::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Turing.Inference ~/.julia/packages/Turing/uAz5c/src/inference/hmc.jl:140
[15] sample
@ ~/.julia/packages/Turing/uAz5c/src/inference/hmc.jl:123 [inlined]
[16] #sample#2
@ ~/.julia/packages/Turing/uAz5c/src/inference/Inference.jl:142 [inlined]
[17] sample
@ ~/.julia/packages/Turing/uAz5c/src/inference/Inference.jl:142 [inlined]
[18] #sample#1
@ ~/.julia/packages/Turing/uAz5c/src/inference/Inference.jl:132 [inlined]
[19] sample(model::DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, alg::NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}, N::Int64)
@ Turing.Inference ~/.julia/packages/Turing/uAz5c/src/inference/Inference.jl:132
[20] (::var"#14#15")(c::Int64)
@ Main ~/Desktop/projects/lDCM/src/test.jl:202
[21] _mapreduce(f::var"#14#15", op::typeof(chainscat), #unused#::IndexLinear, A::UnitRange{Int64})
@ Base ./reduce.jl:410
[22] _mapreduce_dim(f::Function, op::Function, #unused#::Base._InitialValue, A::UnitRange{Int64}, #unused#::Colon)
@ Base ./reducedim.jl:330
[23] #mapreduce#725
@ ./reducedim.jl:322 [inlined]
[24] mapreduce(f::Function, op::Function, A::UnitRange{Int64})
@ Base ./reducedim.jl:322
[25] top-level scope
@ ~/Desktop/projects/lDCM/src/test.jl:202
``````

I don’t understand why this line

``````data[:,i] ~ MvNormal(predicted[i], σ)
``````

would throw a method error only sometimes. Or am I misunderstanding what it’s telling me?

Ok, I progressed a little bit… onto the next error message

The BoundsError seemed to result from a problem in my definition of the Turing model. I changed the last 3 lines to the following and no longer get the error:

``````    for i = 1:size(predicted, 1)
data[i,:] ~ MvNormal(predicted[i,:], σ)
end
``````

But now I get the following:

``````ode.p[1] = -1.1
ode.p[1] = -0.5425958656908638
ode.p[1] = -0.5425958656908638
ode.p[1] = Dual{ForwardDiff.Tag{Turing.Core.var"#f#7"{DynamicPPL.TypedVarInfo{NamedTuple{(:a11, :σ), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:a11, Tuple{}}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:a11, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:σ, Tuple{}}, Int64}, Vector{LogNormal{Float64}}, Vector{AbstractPPL.VarName{:σ, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}}(-0.5425958656908638,1.0,0.0)
ERROR: TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 2}
Stacktrace:
[1] setindex!(A::Matrix{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#7"{DynamicPPL.TypedVarInfo{NamedTuple{(:a11, :σ), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:a11, Tuple{}}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:a11, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:σ, Tuple{}}, Int64}, Vector{LogNormal{Float64}}, Vector{AbstractPPL.VarName{:σ, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 2}, i1::Int64)
@ Base ./array.jl:903
[2] _unsafe_copyto!(dest::Matrix{Float64}, doffs::Int64, src::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#7"{DynamicPPL.TypedVarInfo{NamedTuple{(:a11, :σ), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:a11, Tuple{}}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:a11, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:σ, Tuple{}}, Int64}, Vector{LogNormal{Float64}}, Vector{AbstractPPL.VarName{:σ, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 2}}, soffs::Int64, n::Int64)
@ Base ./array.jl:253
[3] unsafe_copyto!
@ ./array.jl:307 [inlined]
[4] _copyto_impl!
@ ./array.jl:331 [inlined]
[5] copyto!
@ ./array.jl:317 [inlined]
[6] copyto!
@ ./array.jl:343 [inlined]
[7] copyto_axcheck!(dest::Matrix{Float64}, src::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#7"{DynamicPPL.TypedVarInfo{NamedTuple{(:a11, :σ), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:a11, Tuple{}}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:a11, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:σ, Tuple{}}, Int64}, Vector{LogNormal{Float64}}, Vector{AbstractPPL.VarName{:σ, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 2}})
@ Base ./abstractarray.jl:1104
[8] Array
@ ./array.jl:563 [inlined]
[9] convert(#unused#::Type{Matrix{Float64}}, a::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#7"{DynamicPPL.TypedVarInfo{NamedTuple{(:a11, :σ), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:a11, Tuple{}}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:a11, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:σ, Tuple{}}, Int64}, Vector{LogNormal{Float64}}, Vector{AbstractPPL.VarName{:σ, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 2}})
@ Base ./array.jl:554
[10] setproperty!(x::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, false, Matrix{Float64}, Nothing, Float64, Vector{Any}, Float64, Float64, Float64, Vector{Matrix{Float64}}, OrdinaryDiffEq.ODECompositeSolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, false, Vector{Any}, ODEFunction{false, 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}, CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64, SciMLBase.TimeDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, SciMLBase.UDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Matrix{Float64}, LU{Float64, Matrix{Float64}}, DefaultLinSolve}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64, SciMLBase.TimeDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, SciMLBase.UDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Matrix{Float64}, LU{Float64, Matrix{Float64}}, DefaultLinSolve}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}}, Matrix{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, f::Symbol, v::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#7"{DynamicPPL.TypedVarInfo{NamedTuple{(:a11, :σ), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:a11, Tuple{}}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:a11, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:σ, Tuple{}}, Int64}, Vector{LogNormal{Float64}}, Vector{AbstractPPL.VarName{:σ, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{var"#12#13", (:data, :params, :opt), (), (), Tuple{Matrix{Float64}, Vector{Any}, options}, Tuple{}}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 2}})
@ Base ./Base.jl:43
[11] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, false, Matrix{Float64}, Nothing, Float64, Vector{Any}, Float64, Float64, Float64, Vector{Matrix{Float64}}, OrdinaryDiffEq.ODECompositeSolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, false, Vector{Any}, ODEFunction{false, 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}, CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64, SciMLBase.TimeDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, SciMLBase.UDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Matrix{Float64}, LU{Float64, Matrix{Float64}}, DefaultLinSolve}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64, SciMLBase.TimeDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, SciMLBase.UDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Matrix{Float64}, LU{Float64, Matrix{Float64}}, DefaultLinSolve}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}}, Matrix{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/perform_step/low_order_rk_perform_step.jl:565
[12] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, false, Matrix{Float64}, Nothing, Float64, Vector{Any}, Float64, Float64, Float64, Vector{Matrix{Float64}}, OrdinaryDiffEq.ODECompositeSolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, false, Vector{Any}, ODEFunction{false, 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}, CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64, SciMLBase.TimeDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, SciMLBase.UDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Matrix{Float64}, LU{Float64, Matrix{Float64}}, DefaultLinSolve}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64, SciMLBase.TimeDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, SciMLBase.UDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Matrix{Float64}, LU{Float64, Matrix{Float64}}, DefaultLinSolve}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}}, Matrix{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64, SciMLBase.TimeDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Matrix{Float64}, Vector{Any}}, SciMLBase.UDerivativeWrapper{ODEFunction{false, typeof(f!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, Matrix{Float64}, LU{Float64, Matrix{Float64}}, DefaultLinSolve}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/perform_step/composite_perform_step.jl:39
[13] __init(prob::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, false, Vector{Any}, ODEFunction{false, 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}, alg::CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, tstops::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, 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::Float64, dtmin::Nothing, dtmax::Float64, 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(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.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/solve.jl:455
[14] #__solve#408
@ ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/solve.jl:4 [inlined]
[15] #solve_call#56
@ ~/.julia/packages/DiffEqBase/Nyo4y/src/solve.jl:61 [inlined]
[16] solve_up(prob::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, false, Vector{Any}, ODEFunction{false, 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}, sensealg::Nothing, u0::Matrix{Float64}, p::Vector{Any}, args::CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, true, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:saveat, :tstops, :dt), Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64}}})
@ DiffEqBase ~/.julia/packages/DiffEqBase/Nyo4y/src/solve.jl:82
[17] #solve#57
@ ~/.julia/packages/DiffEqBase/Nyo4y/src/solve.jl:70 [inlined]
[18] forward_model(params::Vector{Any}, states::Vector{Matrix{Float64}}, opt::options)
@ Main ~/Desktop/projects/lDCM/src/test.jl:97
[19] #12
@ ~/Desktop/projects/lDCM/src/test.jl:160 [inlined]
@ Main ./none:0
[21] macro expansion
@ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:0 [inlined]
[22] _evaluate
@ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:154 [inlined]
@ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:127 [inlined]
[24] Model
@ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:92 [inlined]
[25] Model
@ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:98 [inlined]
[26] f
@ ForwardDiff ~/.julia/packages/ForwardDiff/tZ5o1/src/apiutils.jl:37
``````

[had to shorten the message]

Now I suspect that there is an issue with the particular parameter estimate, because it is “Nothing”. I understand that I need to reject this estimate, like in this example: Advanced Usage

But I can’t figure out the correct place to do this in my case, because the problem is not directly inside the Turing model. Do you have any suggestions?

do `typeof(p).(x0)` before putting it into the ODEProblem.

Thanks, that worked!

My next error message is that I need to increase the number of maxiters… but I’ll mark this question as solved while I read up on that.

Thank you again for all your help.

That’s usually either because of a divergent model, or if the model is stiff you need to use a more stable integrator (i.e. not something for non-stiff equations)

1 Like