BoundsError on Bayesian ODE model with Turing

spinystellate · January 6, 2022, 2:51pm

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.

trahflow · January 7, 2022, 10:11am

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)

ChrisRackauckas · January 7, 2022, 11:38am

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

spinystellate · January 7, 2022, 9:23pm

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]
 [15] evaluate_threadunsafe
    @ ~/.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)
    @ Base ./loading.jl:1116

ChrisRackauckas · January 8, 2022, 12:09am

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.

spinystellate · March 21, 2022, 1:44pm

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
    
    # Add noise
    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
Turing.setadbackend(:forwarddiff)

@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]
 [17] evaluate_threadunsafe
    @ ~/.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.

ChrisRackauckas · March 21, 2022, 2:02pm

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

spinystellate · March 21, 2022, 2:30pm

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.

ChrisRackauckas · March 21, 2022, 7:34pm

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

spinystellate · March 21, 2022, 9:11pm

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]
  [5] evaluate_threadunsafe
    @ ~/.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?

spinystellate · March 23, 2022, 10:39am

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]
 [20] (::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.TypedVarInfo{NamedTuple{(:a11, :σ), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:a11, Tuple{}}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:a11, Tuple{}}}, Vector{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}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:σ, Tuple{}}, Int64}, Vector{LogNormal{Float64}}, Vector{AbstractPPL.VarName{:σ, Tuple{}}}, Vector{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}}, Vector{Set{DynamicPPL.Selector}}}}}, 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}}, __sampler__::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, __context__::DynamicPPL.DefaultContext, data::Matrix{Float64}, params::Vector{Any}, opt::options)
    @ 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]
 [23] evaluate_threadunsafe
    @ ~/.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
    @ ~/.julia/packages/Turing/uAz5c/src/core/ad.jl:111 [inlined]
 [27] vector_mode_dual_eval!(f::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}, cfg::ForwardDiff.GradientConfig{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, Vector{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}}}, x::Vector{Float64})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/tZ5o1/src/apiutils.jl:37
 [28] vector_mode_gradient!(result::Vector{Float64}, f::Turing.Core.var"#f#7"...

[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?

ChrisRackauckas · March 23, 2022, 12:10pm

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

spinystellate · March 23, 2022, 7:04pm

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.

ChrisRackauckas · March 24, 2022, 8:47am

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)