Simple 2D reaction-diffusion ODE inside Turing blows up in various ways

New to Julia
1

Hello,

Two very helpful examples from the documentation describe how to 1) embed differential equations inside of a Turing model (https://turing.ml/dev/tutorials/10-bayesian-differential-equations/) and 2) optimize a 2D reaction-diffusion model (https://tutorials.sciml.ai/html/introduction/03-optimizing_diffeq_code.html). I’m mashing up these ideas in a way that seems like it should be straightforward but the ship is springing many leaks.

My reaction-diffusion is simpler than the one in the tutorial, with just one species diffusing and exhibiting logistic growth. I have simplified the code accordingly and made the minor indexing changes to work on a non-square grid. When I run it in an ODE solver the results at least seem to be as expected (though I could use another eye to help verify I have not introduced bugs in the process). The solution benchmarks at <50ms for the given timespan using Tsit5().

Growth rate and carrying capacity will ultimately vary based on a covariate raster. However, in the minimal example here, growth rate is randomly generated to fall between 0 and 1, and carrying capacity is a deterministic function of growth rate (multiply by 30). Furthermore, I do not try to learn these values or this relationship, but take it as known. The only term I now try to learn in Turing is the diffusion coefficient. Making the job even easier, the prior is informative, with a mean (4.0) near the actual value (5.0). I had previously run a similar version of this where R and K were scalars, and the solution converged on d=5.0 with minimal variance. However, it would depend on the initialization, as some chains would fail.

Predictions are evaluated by selecting a fixed grid of samples and then slicing that grid every timestep. I do not add noise to the measurement values. I have also run it with observation points that are completely randomized across time and space and did not see better results.

In one run with a large sample size (10x20x26), soon after selecting the step size the value for d explodes and the program crashes

image
image1632Γ—486 75.1 KB

The same happens with a smaller sample size (5x10x26)

image
image1604Γ—494 72.6 KB

Other questions which may or may not be related:

Is the issue type stability? When I switch to using reversediff backend, it looks like it tries to cast the TrackedReal parameters to Float64 to match the type of the initial condition of the ode problem.

image
image2654Γ—480 75.7 KB

I tried to piece together a solution from the forums (see the model testing_type_promotion in the code below) but am obviously missing something

image
image3316Γ—490 84.5 KB

On a side note, I also get an error when testing the Lotka-Volterra tutorial using reversediff instead of the Zygote/SciMLSensitivity combo. I also tried Zygote/SciMLSensitivity on my code, but it seems to expect u0 to be a vector as opposed to a matrix. When I tried unrolling/re-rolling u and du within the function calls, I think I got other errors and set that aside to try other things.

using Turing
using DifferentialEquations
using Plots
using Random
using BenchmarkTools
using ReverseDiff

# Turing.setadbackend(:reversediff)
# Turing.setrdcache(true)

Random.seed!(10);

const Nr = 50
const Nc = 100

true_R = rand(Beta(2,3), (Nr, Nc))
true_K =  true_R .* 30.0

true_d = 5.0
initial_pop = zeros(Nr, Nc)
initial_pop[1:20,end-30:end] = true_K[1:20,end-30:end]

tspan=(0.0, 25.0)
p = (true_d, true_R, true_K)

function fast_react_diffuse!(du,u,p,t)
  D, R, K = p

  @inbounds for j in 2:Nc-1, i in 2:Nr-1
    du[i,j] = D*(u[i-1,j] + u[i+1,j] + u[i,j+1] + u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for j in 2:Nc-1
    i = 1
    du[1,j] = D*(2u[i+1,j] + u[i,j+1] + u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for j in 2:Nc-1
    i = Nr
    du[end,j] = D*(2u[i-1,j] + u[i,j+1] + u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for i in 2:Nr-1
    j = 1
    du[i,1] = D*(u[i-1,j] + u[i+1,j] + 2u[i,j+1] - 4u[i,j])  +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for i in 2:Nr-1
    j = Nc
    du[i,end] = D*(u[i-1,j] + u[i+1,j] + 2u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds begin
    i = 1; j = 1
    du[1,1] = D*(2u[i+1,j] + 2u[i,j+1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])

    i = 1; j = Nc
    du[1,Nc] = D*(2u[i+1,j] + 2u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])

    i = Nr; j = 1
    du[Nr,1] = D*(2u[i-1,j] + 2u[i,j+1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])

    i = Nr; j = Nc
    du[end,end] = D*(2u[i-1,j] + 2u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
   end
end

prob = ODEProblem(fast_react_diffuse!,initial_pop,tspan,p)
@benchmark solve(prob,Tsit5())

sol = solve(prob,Tsit5(), saveat=1.0)
heatmap(sol[end])


@model function extra_simple_troubleshootfitlv(covariates, response, sample_idx, prob)
  # Prior distributions.
  d ~ Gamma(4, 1)
  sigma ~ InverseGamma(2, 3)
  println("d is ", d)
  println("sigma is ", sigma)
  p = [d, covariates[1], covariates[2]]
  predicted = vec(Array(solve(prob, Tsit5(); p=p, saveat=1.0))[sample_idx])
  # Observations.
  response ~ MvNormal(predicted, sigma^2)
  return nothing
end

# @model function testing_type_promotion(covariates, response, sample_idx, prob)
#   d ~ Gamma(4, 1)
#   sigma ~ InverseGamma(2, 3)
#   R = eltype(d).(covariates[1])
#   K = eltype(d).(covariates[2])
#   # R_est ~ MvNormal(vec(R), 0.000001)
#   # R_est = reshape(R, size(R))
#   # K_est ~ MvNormal(vec(K), 0.000001)
#   # K_est = reshape(K, size(K))
  
#   p = [d, R, K]
#   #p = [d, R_est, K_est]
#   #u0 = eltype(R_est).(prob.u0)
#   u0 = eltype(R).(prob.u0)
#   println("d is ", d)
#   println("sigma is ", sigma)
#   println("u0 eltype is ", eltype(u0))
#   # println("R_est eltype is ", eltype(R_est))
#   # println("K_est eltype is ", eltype(K_est))
#   println("R_est eltype is ", eltype(R))
#   println("K_est eltype is ", eltype(K))
#   prob = remake(prob;u0=u0,p=p)
#   predicted = vec(Array(solve(prob, Tsit5(); p=p, saveat=1.0))[sample_idx])
#   # Observations.
#   response ~ MvNormal(predicted, sigma^2)
#   return nothing
# end

## 10x20x26 ##
#sample_idx = CartesianIndices((5:5:size(sol)[1], 5:5:size(sol)[2], 1:size(sol)[3]))

## 5x10x26 ##
sample_idx = CartesianIndices((5:10:size(sol)[1], 5:10:size(sol)[2], 1:size(sol)[3]))

ts_extra_simple_model_1 = extra_simple_troubleshootfitlv((p[2], p[3]), vec(Array(sol)[sample_idx]), sample_idx, prob)
#ts_type = testing_type_promotion((p[2], p[3]), vec(Array(sol)[sample_idx]), sample_idx, prob)

chain_1 = sample(ts_extra_simple_model_1, NUTS(0.65), MCMCSerial(), 500, 1; progress=true)
#chain_2 = sample(ts_type, NUTS(0.65), MCMCSerial(), 500, 1; progress=true)
2

Can you print predicted and response? FWIW your code looks sane to me, so off the top of my head I cannot see why the sampler would go haywire.

It shouldn’t? Share the error with that.

3

Here is my Lotka-Volterra tutorial code using reversediff:

using Turing
using DifferentialEquations
using StatsPlots
using LinearAlgebra
using Random
using Zygote
using ReverseDiff
using SciMLSensitivity

# Turing.setadbackend(:forwarddiff)
# Turing.setadbackend(:zygote)
Turing.setadbackend(:reversediff)
Random.seed!(14);

# Define Lotka-Volterra model.
function lotka_volterra(du, u, p, t)
    # Model parameters.
    Ξ±, Ξ², Ξ³, Ξ΄ = p
    # Current state.
    x, y = u

    # Evaluate differential equations.
    du[1] = (Ξ± - Ξ² * y) * x # prey
    du[2] = (Ξ΄ * x - Ξ³) * y # predator

    return nothing
end

# Define initial-value problem.
u0 = [1.0, 1.0]
p = [1.5, 1.0, 3.0, 1.0]
tspan = (0.0, 10.0)
prob = ODEProblem(lotka_volterra, u0, tspan, p)

# Plot simulation.
plot(solve(prob, Tsit5()))

sol = solve(prob, Tsit5(); saveat=0.1)
odedata = Array(sol) + 0.8 * randn(size(Array(sol)))

# Plot simulation and noisy observations.
plot(sol; alpha=0.3)
scatter!(sol.t, odedata'; color=[1 2], label="")

@model function fitlv(data, prob)
    # Prior distributions.
    Οƒ ~ InverseGamma(2, 3)
    Ξ± ~ truncated(Normal(1.5, 0.5), 0.5, 2.5)
    Ξ² ~ truncated(Normal(1.2, 0.5), 0, 2)
    Ξ³ ~ truncated(Normal(3.0, 0.5), 1, 4)
    Ξ΄ ~ truncated(Normal(1.0, 0.5), 0, 2)

    # Simulate Lotka-Volterra model. 
    p = [Ξ±, Ξ², Ξ³, Ξ΄]
    predicted = solve(prob, Tsit5(); p=p, saveat=0.1)

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

    return nothing
end

model = fitlv(odedata, prob)

# Sample 3 independent chains with forward-mode automatic differentiation (the default).
chain = sample(model, NUTS(0.65), MCMCSerial(), 1000, 3; progress=false)

And the error that I get:

ERROR: MethodError: promote_rule(::Type{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, ReverseDiff.TrackedReal{Float64, Float64, Nothing}, 1}}, ::Type{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}) is ambiguous. Candidates:
  promote_rule(::Type{R}, ::Type{ReverseDiff.TrackedReal{V, D, O}}) where {R<:Real, V, D, O} in ReverseDiff at C:\Users\Brendan\.julia\packages\ReverseDiff\5MMPp\src\tracked.jl:276
  promote_rule(::Type{ForwardDiff.Dual{T, V, N}}, ::Type{R}) where {T, V, N, R<:Real} in ForwardDiff at C:\Users\Brendan\.julia\packages\ForwardDiff\pDtsf\src\dual.jl:426
Possible fix, define
  promote_rule(::Type{ForwardDiff.Dual{T, V, N}}, ::Type{ReverseDiff.TrackedReal{V, D, O}}) where {T, V, N, V, D, O}
Stacktrace:
  [1] promote_type(#unused#::Type{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, #unused#::Type{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, ReverseDiff.TrackedReal{Float64, Float64, Nothing}, 1}})
    @ Base .\promotion.jl:298
  [2] promote_eltype(#unused#::Type{Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}}, #unused#::Type{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, ReverseDiff.TrackedReal{Float64, Float64, Nothing}, 1}})
    @ ArrayInterfaceCore C:\Users\Brendan\.julia\packages\ArrayInterfaceCore\j22dF\src\ArrayInterfaceCore.jl:143
  [3] wrapfun_iip(ff::Function, inputs::Tuple{Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, Float64})
    @ DiffEqBase C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\norecompile.jl:58
  [4] promote_f(f::ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, #unused#::Val{SciMLBase.AutoSpecialize}, u0::Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, p::Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, t::Float64)
    @ DiffEqBase C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\solve.jl:992
  [5] get_concrete_problem(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, isadapt::Bool; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:u0, :p, :saveat), Tuple{Vector{Float64}, Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, Float64}}})
    @ DiffEqBase C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\solve.jl:916
  [6] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ DiffEqBase C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\solve.jl:831
  [7] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; sensealg::Nothing, u0::Nothing, p::Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ DiffEqBase C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\solve.jl:801
  [8] fitlv(__model__::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__::DynamicPPL.TypedVarInfo{NamedTuple{(:Οƒ, :Ξ±, :Ξ², :Ξ³, :Ξ΄), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Οƒ, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:Οƒ, Setfield.IdentityLens}}, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ±, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ±, Setfield.IdentityLens}}, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ², Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ², Setfield.IdentityLens}}, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ³, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ³, Setfield.IdentityLens}}, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ΄, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ΄, Setfield.IdentityLens}}, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}}}, ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}}, __context__::DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random._GLOBAL_RNG}, data::Matrix{Float64}, prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem})
    @ Main .\REPL[26]:11
  [9] macro expansion
    @ C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:593 [inlined]
 [10] _evaluate!!
    @ C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:576 [inlined]
 [11] evaluate_threadunsafe!!
    @ C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:551 [inlined]
 [12] evaluate!!
    @ C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:504 [inlined]
 [13] evaluate!!
    @ C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:515 [inlined]
 [14] evaluate!!
    @ C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:523 [inlined]
 [15] LogDensityFunction
    @ C:\Users\Brendan\.julia\packages\Turing\szPqN\src\Turing.jl:38 [inlined]
 [16] logdensity
    @ C:\Users\Brendan\.julia\packages\Turing\szPqN\src\Turing.jl:42 [inlined]
 [17] Fix1
    @ .\operators.jl:1096 [inlined]
 [18] ReverseDiff.GradientTape(f::Base.Fix1{typeof(LogDensityProblems.logdensity), Turing.LogDensityFunction{DynamicPPL.TypedVarInfo{NamedTuple{(:Οƒ, :Ξ±, :Ξ², :Ξ³, :Ξ΄), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Οƒ, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:Οƒ, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ±, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ±, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ², Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ², Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ³, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ³, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ΄, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ΄, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}}, input::Vector{Float64}, cfg::ReverseDiff.GradientConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}})     
    @ ReverseDiff C:\Users\Brendan\.julia\packages\ReverseDiff\5MMPp\src\api\tape.jl:199
 [19] gradient!(result::DiffResults.MutableDiffResult{1, Float64, Tuple{Vector{Float64}}}, f::Function, input::Vector{Float64}, cfg::ReverseDiff.GradientConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}) (repeats 2 times)
    @ ReverseDiff C:\Users\Brendan\.julia\packages\ReverseDiff\5MMPp\src\api\gradients.jl:41
 [20] logdensity_and_gradient
    @ C:\Users\Brendan\.julia\packages\LogDensityProblems\oAYeE\src\AD_ReverseDiff.jl:55 [inlined]
 [21] βˆ‚logΟ€βˆ‚ΞΈ
    @ C:\Users\Brendan\.julia\packages\Turing\szPqN\src\inference\hmc.jl:166 [inlined]
 [22] βˆ‚Hβˆ‚ΞΈ
    @ C:\Users\Brendan\.julia\packages\AdvancedHMC\iWHPQ\src\hamiltonian.jl:31 [inlined]
 [23] phasepoint(h::AdvancedHMC.Hamiltonian{AdvancedHMC.DiagEuclideanMetric{Float64, Vector{Float64}}, Base.Fix1{typeof(LogDensityProblems.logdensity), LogDensityProblems.ReverseDiffLogDensity{Turing.LogDensityFunction{DynamicPPL.TypedVarInfo{NamedTuple{(:Οƒ, :Ξ±, :Ξ², :Ξ³, :Ξ΄), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Οƒ, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:Οƒ, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ±, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ±, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ², Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ², Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ³, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ³, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ΄, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ΄, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Nothing}}, Turing.Inference.var"#βˆ‚logΟ€βˆ‚ΞΈ#44"{LogDensityProblems.ReverseDiffLogDensity{Turing.LogDensityFunction{DynamicPPL.TypedVarInfo{NamedTuple{(:Οƒ, :Ξ±, :Ξ², :Ξ³, :Ξ΄), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Οƒ, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:Οƒ, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ±, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ±, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ², Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ², Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ³, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ³, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ΄, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ΄, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Nothing}}}, ΞΈ::Vector{Float64}, r::Vector{Float64})
    @ AdvancedHMC C:\Users\Brendan\.julia\packages\AdvancedHMC\iWHPQ\src\hamiltonian.jl:76
 [24] phasepoint
    @ C:\Users\Brendan\.julia\packages\AdvancedHMC\iWHPQ\src\hamiltonian.jl:153 [inlined]
 [25] initialstep(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, spl::DynamicPPL.Sampler{NUTS{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.DiagEuclideanMetric}}, vi::DynamicPPL.TypedVarInfo{NamedTuple{(:Οƒ, :Ξ±, :Ξ², :Ξ³, :Ξ΄), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Οƒ, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:Οƒ, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ±, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ±, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ², Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ², Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ³, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ³, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Ξ΄, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:Ξ΄, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}; init_params::Nothing, nadapts::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Turing.Inference C:\Users\Brendan\.julia\packages\Turing\szPqN\src\inference\hmc.jl:170
 [26] step(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, spl::DynamicPPL.Sampler{NUTS{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.DiagEuclideanMetric}}; resume_from::Nothing, init_params::Nothing, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}})
    @ DynamicPPL C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\sampler.jl:104
 [27] macro expansion
    @ C:\Users\Brendan\.julia\packages\AbstractMCMC\fnRmh\src\sample.jl:120 [inlined]
 [28] macro expansion
    @ C:\Users\Brendan\.julia\packages\AbstractMCMC\fnRmh\src\logging.jl:16 [inlined]
 [29] mcmcsample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.DiagEuclideanMetric}}, N::Int64; progress::Bool, progressname::String, callback::Nothing, discard_initial::Int64, thinning::Int64, chain_type::Type, kwargs::Base.Pairs{Symbol, Union{Nothing, Int64}, Tuple{Symbol, Symbol}, NamedTuple{(:nadapts, :init_params), Tuple{Int64, Nothing}}})
    @ AbstractMCMC C:\Users\Brendan\.julia\packages\AbstractMCMC\fnRmh\src\sample.jl:111
 [30] #sample#42
    @ C:\Users\Brendan\.julia\packages\Turing\szPqN\src\inference\hmc.jl:133 [inlined]
 [31] (::AbstractMCMC.var"#sample_chain#78"{String, Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:chain_type, :progress), Tuple{UnionAll, Bool}}}, Random._GLOBAL_RNG, DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.DiagEuclideanMetric}}, Int64, Int64})(i::Int64, seed::UInt64, init_params::Nothing)
    @ AbstractMCMC C:\Users\Brendan\.julia\packages\AbstractMCMC\fnRmh\src\sample.jl:506
 [32] sample_chain
    @ C:\Users\Brendan\.julia\packages\AbstractMCMC\fnRmh\src\sample.jl:503 [inlined]
 [33] #4
    @ .\generator.jl:36 [inlined]
 [34] iterate
    @ .\generator.jl:47 [inlined]
 [35] collect(itr::Base.Generator{Base.Iterators.Zip{Tuple{UnitRange{Int64}, Vector{UInt64}}}, Base.var"#4#5"{AbstractMCMC.var"#sample_chain#78"{String, Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:chain_type, :progress), Tuple{UnionAll, Bool}}}, Random._GLOBAL_RNG, DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.DiagEuclideanMetric}}, Int64, Int64}}})
    @ Base .\array.jl:787
 [36] map
    @ .\abstractarray.jl:3055 [inlined]
 [37] mcmcsample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.DiagEuclideanMetric}}, ::MCMCSerial, N::Int64, nchains::Int64; progressname::String, init_params::Nothing, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:chain_type, :progress), Tuple{UnionAll, Bool}}})
    @ AbstractMCMC C:\Users\Brendan\.julia\packages\AbstractMCMC\fnRmh\src\sample.jl:518
 [38] sample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.DiagEuclideanMetric}}, ensemble::MCMCSerial, N::Int64, n_chains::Int64; chain_type::Type, progress::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Turing.Inference C:\Users\Brendan\.julia\packages\Turing\szPqN\src\inference\Inference.jl:220
 [39] #sample#6
    @ C:\Users\Brendan\.julia\packages\Turing\szPqN\src\inference\Inference.jl:205 [inlined]
 [40] #sample#5
    @ C:\Users\Brendan\.julia\packages\Turing\szPqN\src\inference\Inference.jl:192 [inlined]
 [41] top-level scope
    @ REPL[29]:1

The tutorial does work with the Zygote/SciMLSensitivity backend:

using Turing
using DifferentialEquations
using StatsPlots
using LinearAlgebra
using Random
using Zygote
using ReverseDiff
using SciMLSensitivity

# Turing.setadbackend(:forwarddiff)
Turing.setadbackend(:zygote)
# Turing.setadbackend(:reversediff)
Random.seed!(14);

# Define Lotka-Volterra model.
function lotka_volterra(du, u, p, t)
    # Model parameters.
    Ξ±, Ξ², Ξ³, Ξ΄ = p
    # Current state.
    x, y = u

    # Evaluate differential equations.
    du[1] = (Ξ± - Ξ² * y) * x # prey
    du[2] = (Ξ΄ * x - Ξ³) * y # predator

    return nothing
end

# Define initial-value problem.
u0 = [1.0, 1.0]
p = [1.5, 1.0, 3.0, 1.0]
tspan = (0.0, 10.0)
prob = ODEProblem(lotka_volterra, u0, tspan, p)

# Plot simulation.
plot(solve(prob, Tsit5()))

sol = solve(prob, Tsit5(); saveat=0.1)
odedata = Array(sol) + 0.8 * randn(size(Array(sol)))

# Plot simulation and noisy observations.
plot(sol; alpha=0.3)
scatter!(sol.t, odedata'; color=[1 2], label="")

@model function fitlv(data, prob)
    # Prior distributions.
    Οƒ ~ InverseGamma(2, 3)
    Ξ± ~ truncated(Normal(1.5, 0.5), 0.5, 2.5)
    Ξ² ~ truncated(Normal(1.2, 0.5), 0, 2)
    Ξ³ ~ truncated(Normal(3.0, 0.5), 1, 4)
    Ξ΄ ~ truncated(Normal(1.0, 0.5), 0, 2)

    # Simulate Lotka-Volterra model. 
    p = [Ξ±, Ξ², Ξ³, Ξ΄]
    predicted = solve(prob, Tsit5(); p=p, saveat=0.1)

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

    return nothing
end

model = fitlv(odedata, prob)

# Sample 3 independent chains with forward-mode automatic differentiation (the default).
chain = sample(model, NUTS(0.65), MCMCSerial(), 1000, 3; progress=false)

(The sampler on this one looks fine – I’m not running it to completion but I have before).

For my code, I shared the error I got with reversediff before. Here is my code using Zygote/SciMLSensitivity:

using Turing
using DifferentialEquations
using Plots
using Random
using BenchmarkTools
using Zygote
using SciMLSensitivity

Turing.setadbackend(:zygote)

Random.seed!(10);

const Nr = 50
const Nc = 100

true_R = rand(Beta(2,3), (Nr, Nc))
true_K =  true_R .* 30.0

true_d = 5.0
initial_pop = zeros(Nr, Nc)
initial_pop[1:20,end-30:end] = true_K[1:20,end-30:end]

tspan=(0.0, 25.0)
p = (true_d, true_R, true_K)

function fast_react_diffuse!(du,u,p,t)
  D, R, K = p

  @inbounds for j in 2:Nc-1, i in 2:Nr-1
    du[i,j] = D*(u[i-1,j] + u[i+1,j] + u[i,j+1] + u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for j in 2:Nc-1
    i = 1
    du[1,j] = D*(2u[i+1,j] + u[i,j+1] + u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for j in 2:Nc-1
    i = Nr
    du[end,j] = D*(2u[i-1,j] + u[i,j+1] + u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for i in 2:Nr-1
    j = 1
    du[i,1] = D*(u[i-1,j] + u[i+1,j] + 2u[i,j+1] - 4u[i,j])  +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for i in 2:Nr-1
    j = Nc
    du[i,end] = D*(u[i-1,j] + u[i+1,j] + 2u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds begin
    i = 1; j = 1
    du[1,1] = D*(2u[i+1,j] + 2u[i,j+1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])

    i = 1; j = Nc
    du[1,Nc] = D*(2u[i+1,j] + 2u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])

    i = Nr; j = 1
    du[Nr,1] = D*(2u[i-1,j] + 2u[i,j+1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])

    i = Nr; j = Nc
    du[end,end] = D*(2u[i-1,j] + 2u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
   end
end

prob = ODEProblem(fast_react_diffuse!,initial_pop,tspan,p)
# @benchmark solve(prob,Tsit5())

sol = solve(prob,Tsit5(), saveat=1.0)
heatmap(sol[end])


@model function extra_simple_troubleshootfitlv(covariates, response, sample_idx, prob)
  # Prior distributions.
  d ~ Gamma(4, 1)
  sigma ~ InverseGamma(2, 3)
  println("d is ", d)
  println("sigma is ", sigma)
  p = [d, covariates[1], covariates[2]]
  predicted = vec(Array(solve(prob, Tsit5(); p=p, saveat=1.0))[sample_idx])
  # Observations.
  response ~ MvNormal(predicted, sigma^2)
  return nothing
end

## 5x10x26 ##
sample_idx = CartesianIndices((5:10:size(sol)[1], 5:10:size(sol)[2], 1:size(sol)[3]))

ts_extra_simple_model_1 = extra_simple_troubleshootfitlv((p[2], p[3]), vec(Array(sol)[sample_idx]), sample_idx, prob)

chain_1 = sample(ts_extra_simple_model_1, NUTS(0.65), MCMCSerial(), 500, 1; progress=true)

Then my VSCode terminal crashes so now I’m not sure of the error anymore. I can come back to this.

For printing predicted and response, I tried outputting to a txt file but it is too long to share. Let me look at a good way of formatting/sharing and get back to you.

4

I can at least share some of the predicted and response.

These were the values drawn before things blew up

d
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.3791585333189884,0.3791585333189884,0.0)

sigma
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.9617349683183385,0.0,1.9617349683183385)

Then a subset of the predicted values (first 100)

Predicted
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(10.996008916384207,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.4758436363252394,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(6.308081048992895,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(26.285385158954853,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(5.896255577251914,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(11.33338836812664,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.2416606942549635e-139,2.1070015443713614e-137,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.166886054944983e-138,7.584098908353453e-137,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.1401958996101837e-102,1.1758877698103884e-100,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.20719243914809e-102,3.9551469840483677e-100,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.379149698324452e-112,2.025573398123166e-110,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(5.534443827623526e-76,2.475572366111883e-74,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.5758677145300252e-75,7.0217318629410726e-74,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.6491050981698866e-82,1.3128968839403488e-80,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.2896089522501376e-105,1.9718398161464823e-103,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0,0.0,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(6.079761250559094e-54,2.097078558799692e-52,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.5279311475928817e-53,5.2368027472383556e-52,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.900604071589555e-59,1.9144492283461144e-57,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.791312201763469e-76,1.3829038713869353e-74,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.208264688050148e-102,1.9230717192118838e-100,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.874388531163025e-35,1.182367965462392e-33,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.1251508348995277e-34,2.7016646728027533e-33,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.136707925072606e-39,6.108911759729404e-38,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.8751051135352945e-54,1.1219334303963646e-52,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.5623754020347375e-74,7.740253402104594e-73,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.1142933736260985e-18,3.0415182624867746e-17,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(5.3134470469774896e-18,7.516813462703795e-17,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.282549630849772e-23,1.3425205696912213e-21,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.177027069665765e-36,6.209965586787982e-35,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.9068517426055015e-54,1.1344508783860502e-52,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.000396283062889018,0.0017755444420739187,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0010235457718068195,0.004291050808435631,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.0095761603012979e-8,8.856284624212889e-8,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.121264757732096e-23,1.314147554127703e-21,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.106974410437914e-39,6.028664708871844e-38,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(13.085249620031133,0.21234003163514348,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.644924612861235,2.550603932164527,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0006802324730057782,0.002853409317608011,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.011323944054808e-18,5.662459118397003e-17,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(8.557043496048169e-35,2.0511099654409817e-33,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(11.398281524598621,1.3589046320454705,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(20.918541084417885,-3.733299574113343,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0007650153872988737,0.0032185356853350073,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.3549465457422924e-18,6.157091106390481e-17,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.506441192025142e-35,2.279759374429026e-33,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.048909391652534,1.5969965307168343,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(11.514512267993823,-0.08439149296649791,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0008018483599065973,0.003595622505931865,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.106544454671553e-18,5.918341531157599e-17,0.0)

Then the first 100 lines of the response

Response
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
10.996008916384207
2.4758436363252394
0.0
0.0
0.0
6.308081048992895
26.285385158954853
0.0
0.0
0.0
5.896255577251914
11.33338836812664
0.0
0.0
0.0
5.158745691859671e-52
7.482867906941866e-52
8.283578017290905e-54
-1.6245002082237417e-59
-6.364120847458168e-69
1.9897759776140717e-39
2.619625480809167e-39
6.567679886024994e-41
2.098396613924575e-46
-6.713337919170827e-55
1.5640663097818284e-28
1.925240293493555e-28
6.878511452732576e-30
1.2902685565422506e-34
1.529644983123486e-42
4.6446303434164065e-19
5.425943872649536e-19
2.3287281456389762e-20
1.0331750980928464e-24
8.884501556158954e-32
5.0216908413387066e-11
5.748962499814938e-11
2.7806502471382172e-12
1.775527328044353e-16
4.191511764036631e-23
0.00011774367245424741
0.0001311181908672811
6.812955072836514e-6
5.105182933476234e-10
1.822224039067667e-16
1.2363041122306637
1.2455661683757266
0.07659794983215734
6.622679398590532e-6
2.7236881979256756e-12
12.49697913522685
10.398099024103747
1.035205481730887
0.00010470148897507296
4.510456694441825e-11
11.597151091214798
13.068045971520174
1.1962054766894918
0.00011698655266152904
4.991641851474777e-11
11.41164055130108
11.57194746259485
1.3243361770047541
0.00012685367004336626

Immediately after that, d skyrocketed

d
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4003.3679877682116,4003.3679877682116,0.0)

sigma
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.665262686461296e14,0.0,3.665262686461296e14)

5

Sorry, I realized the first 100 lines of the predicted/response don’t show anything meaningful as that’s the initial condition. Here are the last 250 lines of each instead. Still not sure if they’re showing anything. Is there a reason the response does not print out as having dual type? I’m writing both to file with DelimitedFiles after response is drawn.

Predicted:

Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(6.355432759082686e-35,3.281754139127906e-33,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.583034247346742e-35,2.3555402024118673e-33,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.757539430714259e-36,1.992222234137751e-34,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.823841166940317e-39,4.57873945229097e-37,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.4021654442798415e-44,2.9535143445082757e-42,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(6.855499693229318e-26,2.8573441389986346e-24,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.700535686534377e-26,1.9460445224642187e-24,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.959555182653496e-27,1.703564079945879e-25,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(6.756445515043082e-30,3.2724419115869294e-28,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(8.382366954598051e-35,4.749116335238938e-33,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.241531402450233e-17,3.9562513385367306e-16,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(6.819563420352901e-18,2.1663898122079392e-16,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(8.490730594384415e-19,2.8168483478110096e-17,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.1984133440896357e-21,4.6181000819295395e-20,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.7781748196461583e-26,8.281869124818016e-25,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.2135060865230044e-10,7.17981828723022e-9,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.2142472598433264e-10,2.695810660685382e-9,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.662535922005162e-11,3.951836832052801e-10,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.728310088408641e-14,7.909010298985563e-13,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.2153867402653155e-19,1.188069706959307e-17,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0002519379382526913,0.003369645164762107,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.000124618052672406,0.0016329485781527943,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.6268332630993628e-5,0.0003852516601918184,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.0144247001815504e-8,4.0542601695065975e-7,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(6.104264707651398e-13,1.7016860191549794e-11,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.0638645284792885,7.2986325769537,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.8581937866022125,6.649510928082002,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.5945722217477435,3.3574789002406753,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.00109319488532443,0.013125568070021009,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.735509871790929e-8,5.513752047744399e-7,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(13.922674008718728,-1.3069402735419058,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(13.566199834540601,-0.40829087788379337,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(8.86436857350785,3.236806801588091,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.3948220898673393,2.3140063028322633,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.462995305275182e-5,0.00021429708191235525,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(13.17097925211459,-0.028876800904008352,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.371632217911243,1.2611416387841337,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(12.747889792064079,0.34643692861175657,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.1500973298307753,7.553969226441947,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.00016982429201689182,0.002259568849125025,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(11.844720279611632,0.27083131342544925,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(19.31597800597822,-3.6201512735656687,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.239387616342999,1.6259281227183537,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.3881447142203296,5.600817114635972,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0001209214203318883,0.0016322986164855526,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(10.33644606411558,0.9233759449016932,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(11.745886957242341,0.10582775884626527,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(14.03803601326141,-0.701339540746773,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.708060556123688,6.225140685070612,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.00010873854323604694,0.0014407986068482382,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.134175271242297e-33,5.781585852934261e-32,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.878195871474277e-34,3.997584659879013e-32,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.001208642457789e-35,3.662683290272585e-33,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.052922157171964e-37,1.1820100646520407e-35,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.801252691864127e-42,1.1875315926747915e-40,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.387233450518143e-25,3.0346381739234186e-23,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.862841414163885e-25,1.9845268419092854e-23,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.4209896948761355e-26,1.8741833554288928e-24,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.0175208930139977e-28,4.8429813362928534e-27,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.9286363769375594e-33,1.0722212726257113e-31,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(8.367758546744332e-17,2.62112584361518e-15,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.417771531888673e-17,1.3802036109256156e-15,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(5.9603032556064484e-18,1.9433095103507511e-16,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.1000826634466988e-20,4.156563283507974e-19,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.4171633566668145e-25,1.1029459026684562e-23,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.4078405469891864e-9,3.076782555948933e-8,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(5.016108556790733e-10,1.0904050304096561e-8,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.522737348875027e-11,1.7499648953998682e-9,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.5996135841604923e-13,4.52946573331415e-12,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.694184890136388e-18,9.731057131072844e-17,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0007054991176262801,0.009133227594651546,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0003320065992694884,0.0042103243417977955,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.849162063562898e-5,0.001116134816163594,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.576847349479712e-8,1.4817972019846795e-6,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.368095673909786e-12,9.145814843865056e-11,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.969940825288798,8.825084571777243,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.692331898509227,8.07314867049301,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.052372839511713,5.276532275761496,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.002817026128322792,0.03244254191833527,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.0470590249924098e-7,2.0508296849698858e-6,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(13.99832073246247,-1.3966200688724428,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(13.630439088573327,-0.4458199007696506,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.303994803261824,2.730691323479096,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.7052378994541432,3.752019372416411,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.258744561023888e-5,0.0006049290771241536,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(13.171173827737295,-0.02932498649680574,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.372360261541086,1.2625092508277758,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(12.819331101532713,0.284175270538427,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.1398307076661367,9.3299100939332,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.00046857644792515964,0.006034526685473715,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(11.84474697116731,0.2707428783778031,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(19.31620335258376,-3.619713894230704,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.317614217413428,1.5439641450625916,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.0584159430910534,7.1692797615708095,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.00032931180491699926,0.00431294525902295,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(10.336470219695292,0.923370531715864,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(11.74614903119102,0.10603760456155802,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(14.091824163491804,-0.7444341421880388,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.4627661299040087,7.631081434391661,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0002887514648380227,0.003706046652183742,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.741224460329519e-32,8.76424761244225e-31,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.1666724898401626e-32,5.845830606695054e-31,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.1157319922224323e-33,5.760926391019075e-32,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.4662609179236264e-36,2.5311033113136e-34,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(5.859523868659264e-41,3.797452875568158e-39,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.073114345609339e-24,2.8642583189513653e-22,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.475891898320821e-24,1.8009158387724234e-22,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.366389221019884e-25,1.8243192768680044e-23,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.321543706383217e-27,6.183945945209435e-26,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.6851309085711615e-32,2.011445542310768e-30,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(5.149424976697547e-16,1.585781250238273e-14,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.616743951322381e-16,8.041222029954185e-15,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.805064428582416e-17,1.2194943717521563e-15,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(8.990950756569817e-20,3.332314085978146e-18,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.835730790789624e-24,1.2683196502734193e-22,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(5.764128932049413e-9,1.2322053443268094e-7,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.939571811017047e-9,4.128540709262189e-8,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.1711201998396185e-10,7.220020369925614e-9,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(8.574833825171503e-13,2.3725129303168148e-11,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.011398826026328e-17,7.104689511194179e-16,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0018825773282543532,0.02356999361965455,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0008432597348210618,0.010339147459777412,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.00022301123665490242,0.0030728789417910187,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.6687592732711593e-7,5.072911443692375e-6,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.6998288923416243e-11,4.497632896047928e-10,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.020498576390578,9.827321650813735,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.6572289644716993,8.937911446507332,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.7388557332983905,7.555269057806415,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.006914088951312325,0.07626516381559884,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.753535759393892e-7,7.145910519675586e-6,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(14.050636937521334,-1.4575759599075686,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(13.677520611716453,-0.47094013512663413,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.63924700801146,2.3074235990249647,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.1900826636679935,5.649326042715083,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.00011790888606262855,0.00162344008412975,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(13.17132034545669,-0.030047137591067967,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.372939267442353,1.263568218872203,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(12.869901265105309,0.2420503827294439,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.316108847550203,10.60926456283583,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0012317171433958794,0.015337702521927685,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(11.84476478702761,0.2706770671225334,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(19.316375314298323,-3.6194059399505853,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.374883518719102,1.4844112221306498,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.875568177845108,8.452633941787715,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0008559087751287287,0.01086714607131132,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(10.33648641616552,0.9233710968875102,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(11.746344858922473,0.10618729958569093,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(14.130942956817611,-0.7738501012878242,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.344923234195684,8.62947128885458,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0007306157451128392,0.009073547372281058,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.336334133598353e-31,1.1613627167350933e-29,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.5119358013929128e-31,7.482430051007416e-30,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.5464116160129695e-32,7.882577950172916e-31,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(8.22831636522728e-35,4.5917657920236714e-33,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.5550136752766722e-39,9.911178374174391e-38,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(6.089241436775278e-23,2.4311148521669447e-21,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.7085896429127266e-23,1.4714485397846635e-21,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.86277687610387e-24,1.5909474524393152e-22,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.5044102276363105e-26,6.923730931398634e-25,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(5.97423180700903e-31,3.2031564813926787e-29,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.9186356164485194e-15,8.837183955175258e-14,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.4293973786103963e-15,4.321030587404453e-14,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.2297486028878025e-16,7.025789218385738e-15,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(6.622150356542079e-19,2.4083867626596116e-17,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.918236579523934e-23,1.2799574273661565e-21,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.219302370148746e-8,4.640413401522593e-7,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.0628366745399384e-9,1.4721720911312864e-7,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.253667565616257e-9,2.793913685164565e-8,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.240555532332732e-12,1.1467911929573793e-10,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.3542725711579216e-16,4.6798848789591325e-15,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.004804468324102671,0.05811154698909987,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.002048645211216355,0.024253439102343358,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0006048491006267891,0.008068882238250202,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(8.859103722201409e-7,1.637034427337769e-5,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.917839596032182e-11,2.042004219587965e-9,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(5.148711652558023,10.139144410319808,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.6785577842902315,9.080076134737274,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.676162915190733,9.806361615921498,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.01621131183513963,0.1708886994561717,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.2683443550334062e-6,2.3474425895652888e-5,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(14.086766342779331,-1.4987419710948229,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(13.712095281434898,-0.48811513991809263,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.889713091415384,1.9724681702114888,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.8930140353380942,7.845603711935707,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.00031174326523665297,0.004157936072452817,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(13.171430619762106,-0.030942655146475182,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.373395752314332,1.26437006616687,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(12.905685210406302,0.21379451887324172,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(5.60961960068631,11.145816891681717,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0030946698365697856,0.037208674477225784,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(11.844776495068682,0.2706154389519139,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(19.316505903177344,-3.6191905146037375,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.416688395141035,1.44163984451295,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.795527476301429,9.20360871043704,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0021304536251953105,0.02619558174339611,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(10.33649727289789,0.9233743840526027,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(11.746490042231326,0.10627450343731776,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(14.15951052911117,-0.7940858929924771,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.2998526460737585,9.061218106463382,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0017672837834993183,0.021207141175502667,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.7932467240356258e-30,1.3712973052611431e-28,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.747327467450206e-30,8.541217704696449e-29,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.9023729901500778e-31,9.573834393593088e-30,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.3187558906659243e-33,7.248190854360585e-32,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.4878766096533726e-38,2.1870062070635016e-36,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.773163794225202e-22,1.878909535334584e-20,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.8003380389225323e-22,1.0957215241102007e-20,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.1018854493094937e-23,1.2595052434206333e-21,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.5282560984143963e-25,6.919123469996068e-24,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(8.43006515157204e-30,4.4410050886663065e-28,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.5362646095525695e-14,4.57366019244308e-13,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(7.259208027077429e-15,2.1588575883027315e-13,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.210018766480449e-15,3.7488475875825483e-14,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.448116913875884e-18,1.5877714885000262e-16,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.6799224467349124e-22,1.1529673940130761e-20,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(8.079992202549036e-8,1.652362752969029e-6,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.4352229665268768e-8,4.970570065058298e-7,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.675982794524873e-9,1.0200397796911859e-7,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(1.950525652148532e-11,5.156836785283667e-10,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(8.318768814240707e-16,2.8130847726234905e-14,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.011761363209104689,0.13722920848558617,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.004772854137120947,0.05446239213412362,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0015708407125999255,0.020262998267512336,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.7869228541410107e-6,5.0065508844297836e-5,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.4316305390169253e-10,8.628156579759368e-9,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(6.278313540294772,9.762647063519479,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(5.675534279901218,8.562455693715325,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(3.8378995874763437,11.539001730643834,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.03638010348616297,0.3653040792998988,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.062462235715667e-6,7.310304674786382e-5,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(14.111863189941845,-1.5106019669481427,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(13.737472890527462,-0.5023958004997572,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(10.074148916184882,1.718132386815165,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(2.832886440667798,9.985087750107889,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.000789704019803446,0.010193766775423867,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(13.171766424808292,-0.006890731400878685,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.373743555890321,1.2640381377326764,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(12.931006926010282,0.19372372356326847,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(6.931334030710746,10.863883974632152,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.0074503201790238295,0.08632936306430825,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(11.844768883871913,0.26885319316333095,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(19.31663946655246,-3.6159045489971797,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(9.4471858610208,1.4147979338499317,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(4.755592840389309,9.310259451777453,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.005092910371236793,0.06056219539055844,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(10.336501426292328,0.9230281529920801,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(11.746580610872336,0.1045461023707943,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(14.180448613434535,-0.8079781579033585,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(5.2642858138713695,8.915759183795306,0.0)
Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}}(0.004097195015509178,0.047418174164687486,0.0)

Response

0.10211280144973942
0.07061959049273743
0.04462911105182386
0.024182588002528572
0.00810862262307643
0.9943202725382468
0.7792635129571918
0.4859920064894169
0.2929400547362175
0.12165860119468681
4.6955329048324375
3.8679426889054835
3.2088957623964336
1.8550492219713026
0.9428146393035719
9.972015536153695
8.722315351731101
6.93614165174244
5.767928297740403
3.6431778747745973
12.539909819040428
10.387639901968816
10.968294853344752
8.593927228532621
7.153126457205752
11.530412333086748
11.090926292705033
11.964041261921162
11.751791938708978
10.37126340170122
12.432057031358514
12.208568896647106
11.52778110460857
11.330395672692319
10.453061545431906
12.674431845252137
11.235307660782247
12.21912618288154
11.871304961836229
12.128646377704174
11.727236681868193
13.183260100788138
11.539928804558773
10.923588546490405
12.141503575594074
11.446978454210306
12.058323865132145
12.618765057414416
10.904312404955492
10.79511122930951
0.21157717647245886
0.14786328160117315
0.09683176867639777
0.05480993363508876
0.019828786630687167
1.5948993894998094
1.2850489839649841
0.8331817675265727
0.5406190208988005
0.247300573344411
5.8731888324511665
4.977504466507586
4.3505897618299905
2.717142936618336
1.5589463079667998
10.76440384999049
9.594602016294118
7.8852251318012625
6.944137362496791
4.882043044346234
12.819987516733509
10.710322484377965
11.41848196833447
9.266335505538994
8.206030750253918
11.614233867481188
11.192776987345134
12.130316714839676
12.098416041364185
10.98595221529593
12.461253373039868
12.2452489649169
11.594821110298218
11.4984134669309
10.779498484863367
12.694801504635803
11.256353738919927
12.256449431659524
11.971376996361366
12.341909189512492
11.748676024475218
13.202930449632918
11.568500232953635
11.002816600469423
12.320094583853251
11.469451467967708
12.079547036859374
12.646370333615621
10.980638170818139
10.961038352195953
0.412296700602823
0.29167441375385
0.1985177333963887
0.11748062590190782
0.04583685858753103
2.3934031355734753
1.985283719595586
1.3390615872336147
0.9380479179966627
0.4732625169510161
6.997235697833451
6.083115814996608
5.562093993403132
3.740539236710224
2.4057391569979383
11.391578711434617
10.30588858233973
8.705910048500773
8.029952893057818
6.157812913567856
13.026984249115465
10.954053765075882
11.76503996858823
9.805513035815926
9.086934526616359
11.671237958058384
11.265415112152729
12.252776113629393
12.355422823148835
11.440379323292113
12.468802340038806
12.263507202402499
11.640403698416165
11.619396249267437
11.009417157302448
12.683929533935693
11.25351752531471
12.275682685257886
12.04026357006155
12.486661025413825
11.730845057281305
13.193032885360825
11.576024442818554
11.054837190833986
12.439202832738708
11.44924671183077
12.065219517154672
12.651410358898012
11.030072548009187
11.071930581425109
0.7555739280298398
0.54216935263947
0.3848462829531816
0.23844676473714588
0.10041692670434144
3.371209816342824
2.879096131842102
2.0192698350457845
1.5288203376460394
0.8510187073952071
8.008675290654384
7.116087822401856
6.752007368636479
4.860500890305876
3.465636193668505
11.876536165678322
10.869598474022366
9.387951199463322
8.980056081525966
7.367461038822502
13.18062576729286
11.13746076983738
12.029201020164667
10.227033583051172
9.790384715355652
11.71829420892275
11.322811771648638
12.345094006839492
12.54444890296866
11.769672057228695
12.488893808306925
12.286720343217198
11.678945166734596
11.70792065357136
11.169845703888143
12.701505780750574
11.270098690906595
12.298925957798627
12.092183763634335
12.58544306514345
11.750501960845515
13.21007844005136
11.595697241291031
11.095978200736775
12.519896833111572
11.47022277706387
12.08448622697887
12.671081703455101
11.070208901899457
11.147548074416523
1.3004672237487807
0.948718327519738
0.7050010807309594
0.4581317060456673
0.2086435914251118
4.477919566485019
3.9324455510137457
2.8635897582650633
2.338900109335934
1.4338001093360135
8.876712602646634
8.030588355677587
7.84493579528641
5.996882464897057
4.673242594962372
12.244748459610879
11.306439408041765
9.936123106200336
9.774567294567923
8.435558470369966
13.292014927829673
11.273033070186711
12.228098838392128
10.549689363137206
10.332365967594125
11.747448181359665
11.361468660195941
12.411363240634202
12.6813181657794
12.004857129215683
12.489756862141219
12.29457363769254
11.702464562147744
11.769430456586315
11.280242158585867
12.690218316286632
11.264515931236717
12.307124097275626
12.125471013142025
12.651237253448512
11.734012348885786
13.199416680167454
11.59619612122834
11.119844249791848
12.57240051303817
11.451723006248105
12.070081432804983
12.66976684550827
11.09264288303423
11.196526870881762
6

One last comment: unlike the other initializations, this one spiked but did not exit altogether. I then let it sit and attempt in the background, which it did for about an hour without a single successful proposal. But eventually it discovered a leap back to approximately d=5, and is slowly sampling in the vicinity, with falling variance. So it is CAPABLE of discovering the correct value, contingent on varying amounts of luck and probably 8+ hours of sampling. Which is simultaneously good to discover but also presumably signaling some deeper problem, as it should not be this difficult.

7

I tried autodiff with Tracker as well and got another type-related error. As someone new to Julia, this is the aspect that I least understand, but isn’t this a sign that all of my issues may be related to type instability?

ERROR: TypeError: in TrackedReal, in T, expected T<:Real, got Type{Any}
Stacktrace:
  [1] promote_rule(#unused#::Type{Tracker.TrackedReal{Float64}}, #unused#::Type{Matrix{Float64}})
    @ Tracker C:\Users\Brendan\.julia\packages\Tracker\a9oj5\src\lib\real.jl:64
  [2] promote_type(#unused#::Type{Tracker.TrackedReal{Float64}}, #unused#::Type{Matrix{Float64}})
    @ Base .\promotion.jl:298
  [3] promote_typeof(::Tracker.TrackedReal{Float64}, ::Matrix{Float64}, ::Vararg{Matrix{Float64}})
    @ Base .\promotion.jl:339
  [4] vect(::Tracker.TrackedReal{Float64}, ::Vararg{Any})
    @ Base .\array.jl:144
  [5] extra_simple_troubleshootfitlv(__model__::DynamicPPL.Model{typeof(extra_simple_troubleshootfitlv), (:covariates, :response, :sample_idx, :prob), (), (), Tuple{Tuple{Matrix{Float64}, Matrix{Float64}}, Vector{Float64}, CartesianIndices{3, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}, UnitRange{Int64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__::DynamicPPL.TypedVarInfo{NamedTuple{(:d, :sigma), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:d, Setfield.IdentityLens}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:d, Setfield.IdentityLens}}, TrackedArray{…,Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}}, TrackedArray{…,Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}}}, Tracker.TrackedReal{Float64}}, __context__::DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random._GLOBAL_RNG}, covariates::Tuple{Matrix{Float64}, Matrix{Float64}}, response::Vector{Float64}, sample_idx::CartesianIndices{3, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}, UnitRange{Int64}}}, prob::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem})
    @ Main .\REPL[33]:7
  [6] macro expansion
    @ C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:593 [inlined]
  [7] _evaluate!!(model::DynamicPPL.Model{typeof(extra_simple_troubleshootfitlv), (:covariates, :response, :sample_idx, :prob), (), (), Tuple{Tuple{Matrix{Float64}, Matrix{Float64}}, Vector{Float64}, CartesianIndices{3, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}, UnitRange{Int64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, varinfo::DynamicPPL.TypedVarInfo{NamedTuple{(:d, :sigma), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:d, Setfield.IdentityLens}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:d, Setfield.IdentityLens}}, TrackedArray{…,Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}}, TrackedArray{…,Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}}}, Tracker.TrackedReal{Float64}}, context::DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random._GLOBAL_RNG})
    @ DynamicPPL C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:576
  [8] evaluate_threadunsafe!!(model::DynamicPPL.Model{typeof(extra_simple_troubleshootfitlv), (:covariates, :response, :sample_idx, :prob), (), (), Tuple{Tuple{Matrix{Float64}, Matrix{Float64}}, Vector{Float64}, CartesianIndices{3, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}, UnitRange{Int64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, varinfo::DynamicPPL.TypedVarInfo{NamedTuple{(:d, :sigma), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:d, Setfield.IdentityLens}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:d, Setfield.IdentityLens}}, TrackedArray{…,Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}}, TrackedArray{…,Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}}}, Tracker.TrackedReal{Float64}}, context::DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random._GLOBAL_RNG})
    @ DynamicPPL C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:551
  [9] evaluate!!(model::DynamicPPL.Model{typeof(extra_simple_troubleshootfitlv), (:covariates, :response, :sample_idx, :prob), (), (), Tuple{Tuple{Matrix{Float64}, Matrix{Float64}}, Vector{Float64}, CartesianIndices{3, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}, UnitRange{Int64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, varinfo::DynamicPPL.TypedVarInfo{NamedTuple{(:d, :sigma), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:d, Setfield.IdentityLens}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:d, Setfield.IdentityLens}}, TrackedArray{…,Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}}, TrackedArray{…,Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}}}, Tracker.TrackedReal{Float64}}, context::DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random._GLOBAL_RNG})
    @ DynamicPPL C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:504
 [10] evaluate!!(model::DynamicPPL.Model{typeof(extra_simple_troubleshootfitlv), (:covariates, :response, :sample_idx, :prob), (), (), Tuple{Tuple{Matrix{Float64}, Matrix{Float64}}, Vector{Float64}, CartesianIndices{3, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}, UnitRange{Int64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, rng::Random._GLOBAL_RNG, varinfo::DynamicPPL.TypedVarInfo{NamedTuple{(:d, :sigma), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:d, Setfield.IdentityLens}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:d, Setfield.IdentityLens}}, TrackedArray{…,Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}}, TrackedArray{…,Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}}}, Tracker.TrackedReal{Float64}}, sampler::DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, context::DynamicPPL.DefaultContext)
    @ DynamicPPL C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:515
 [11] evaluate!!(::DynamicPPL.Model{typeof(extra_simple_troubleshootfitlv), (:covariates, :response, :sample_idx, :prob), (), (), Tuple{Tuple{Matrix{Float64}, Matrix{Float64}}, Vector{Float64}, CartesianIndices{3, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}, UnitRange{Int64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, ::DynamicPPL.TypedVarInfo{NamedTuple{(:d, :sigma), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:d, Setfield.IdentityLens}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:d, Setfield.IdentityLens}}, TrackedArray{…,Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}}, TrackedArray{…,Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}}}, Tracker.TrackedReal{Float64}}, ::DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, ::DynamicPPL.DefaultContext)
    @ DynamicPPL C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:523
 [12] (::Turing.LogDensityFunction{DynamicPPL.TypedVarInfo{NamedTuple{(:d, :sigma), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:d, Setfield.IdentityLens}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:d, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(extra_simple_troubleshootfitlv), (:covariates, :response, :sample_idx, :prob), (), (), Tuple{Tuple{Matrix{Float64}, Matrix{Float64}}, Vector{Float64}, CartesianIndices{3, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}, UnitRange{Int64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext})(ΞΈ::TrackedArray{…,Vector{Float64}})
    @ Turing C:\Users\Brendan\.julia\packages\Turing\szPqN\src\Turing.jl:38

This was running my base model from the first post (extra_simple_troubleshootfitlv).

8

The above error may be because in my code for learning type I had been making a vector p = [d, R_est, K_est] rather than a tuple p = (d, R_est, K_est). This was correct in my base function, but the one where I was playing around with typing may have been tainted from the beginning.

This is becoming something of a running diary of my troubleshooting–assuming we solve it, I hope this proves useful to someone. For the moment I stopped passing in any covariates and just generate them as parameters to estimate in Turing. This is where I want to go with the model anyway, and I imagined they would be the correct type. So even though there are more parameters I thought it might be simpler to understand why various forms of reverse differentiation are crashing.

The testing function:

@model function testing_type_promotion(response, Nr, Nc, sample_idx, prob)
  d ~ Gamma(4, 1)
  sigma ~ InverseGamma(2, 3)
  println("d is ", d)
  println("sigma is ", sigma)
  R_est ~ filldist(Beta(2,3), Nr*Nc)
  R_est = reshape(R_est, (Nr, Nc))
  K_est = 30.0 * R_est
  println("R_est eltype is ", eltype(R_est))
  println("K_est eltype is ", eltype(K_est))
  println("Type of R_est is ", typeof(R_est))
  println("Type of K_est is ", typeof(K_est))
  p = (d, R_est, K_est)
  println("Type of p is ", typeof(p))
  u0 = eltype(R_est).(prob.u0)
  println("u0 eltype is ", eltype(u0))
  println("p eltype is ", eltype(p))
  prob = remake(prob;u0=u0,p=p)
  println("remade problem")
  predicted = vec(Array(solve(prob, Tsit5(); p=p, saveat=1.0))[sample_idx])
  # Observations.
  response ~ MvNormal(predicted, sigma^2)
  return nothing
end

Its context:

using Turing
using DifferentialEquations
using Plots
using Random
using BenchmarkTools

# Turing.setadbackend(:reversediff)
# Turing.setrdcache(true)
Turing.setadbackend(:tracker)
Turing.setrdcache(true)

Random.seed!(10);

const Nr = 50
const Nc = 100

true_R = rand(Beta(2,3), (Nr, Nc)) # mean 2/(2+3) = 0.4
true_K =  true_R .* 30.0

true_d = 5.0
initial_pop = zeros(Nr, Nc)
initial_pop[1:20,end-30:end] = true_K[1:20,end-30:end]

tspan=(0.0, 25.0)
p = (true_d, true_R, true_K)

function fast_react_diffuse!(du,u,p,t)
  D, R, K = p

  @inbounds for j in 2:Nc-1, i in 2:Nr-1
    du[i,j] = D*(u[i-1,j] + u[i+1,j] + u[i,j+1] + u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for j in 2:Nc-1
    i = 1
    du[1,j] = D*(2u[i+1,j] + u[i,j+1] + u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for j in 2:Nc-1
    i = Nr
    du[end,j] = D*(2u[i-1,j] + u[i,j+1] + u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for i in 2:Nr-1
    j = 1
    du[i,1] = D*(u[i-1,j] + u[i+1,j] + 2u[i,j+1] - 4u[i,j])  +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for i in 2:Nr-1
    j = Nc
    du[i,end] = D*(u[i-1,j] + u[i+1,j] + 2u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds begin
    i = 1; j = 1
    du[1,1] = D*(2u[i+1,j] + 2u[i,j+1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])

    i = 1; j = Nc
    du[1,Nc] = D*(2u[i+1,j] + 2u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])

    i = Nr; j = 1
    du[Nr,1] = D*(2u[i-1,j] + 2u[i,j+1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])

    i = Nr; j = Nc
    du[end,end] = D*(2u[i-1,j] + 2u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
   end
end

prob = ODEProblem(fast_react_diffuse!,initial_pop,tspan,p)
sol = solve(prob,Tsit5(), saveat=1.0)
heatmap(sol[end])
noisy_sol = Array(sol) + 0.4 * randn(size(Array(sol)))

@model function testing_type_promotion(response, Nr, Nc, sample_idx, prob)
  d ~ Gamma(4, 1)
  sigma ~ InverseGamma(2, 3)
  println("d is ", d)
  println("sigma is ", sigma)
  R_est ~ filldist(Beta(2,3), Nr*Nc)
  R_est = reshape(R_est, (Nr, Nc))
  K_est = 30.0 * R_est
  println("R_est eltype is ", eltype(R_est))
  println("K_est eltype is ", eltype(K_est))
  println("Type of R_est is ", typeof(R_est))
  println("Type of K_est is ", typeof(K_est))
  p = (d, R_est, K_est)
  println("Type of p is ", typeof(p))
  u0 = eltype(R_est).(prob.u0)
  println("u0 eltype is ", eltype(u0))
  println("p eltype is ", eltype(p))
  prob = remake(prob;u0=u0,p=p)
  println("remade problem")
  predicted = vec(Array(solve(prob, Tsit5(); p=p, saveat=1.0))[sample_idx])
  # Observations.
  response ~ MvNormal(predicted, sigma^2)
  return nothing
end

## 5x10x26 ##
sample_idx = CartesianIndices((5:10:size(noisy_sol)[1], 5:10:size(noisy_sol)[2], 1:size(noisy_sol)[3]))
chain_2 = sample(ts_type, NUTS(0.65), MCMCSerial(), 500, 1; progress=true)

Unfortunately, still an error message:

d is 0.17973133345041883
sigma is 1.4483293423794181
R_est eltype is Float64
K_est eltype is Float64
Type of R_est is Matrix{Float64}
Type of K_est is Matrix{Float64}
Type of p is Tuple{Float64, Matrix{Float64}, Matrix{Float64}}
u0 eltype is Float64
p eltype is Any
remade problem
d is 0.17973133345041883 (tracked)
sigma is 1.4483293423794181 (tracked)
R_est eltype is Tracker.TrackedReal{Float64}
K_est eltype is Tracker.TrackedReal{Float64}
Type of R_est is TrackedArray{…,Matrix{Float64}}
Type of K_est is TrackedArray{…,Matrix{Float64}}
Type of p is Tuple{Tracker.TrackedReal{Float64}, TrackedArray{…,Matrix{Float64}}, TrackedArray{…,Matrix{Float64}}}
u0 eltype is Tracker.TrackedReal{Float64}
p eltype is Any
remade problem
ERROR: MethodError: no method matching Float64(::Tracker.TrackedReal{Float64})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:772
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  ...
Stacktrace:
  [1] convert(#unused#::Type{Float64}, x::Tracker.TrackedReal{Float64})
    @ Base .\number.jl:7
  [2] ode_determine_initdt(u0::Matrix{Tracker.TrackedReal{Float64}}, t::Float64, tdir::Float64, dtmax::Float64, abstol::Tracker.TrackedReal{Float64}, reltol::Tracker.TrackedReal{Float64}, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), prob::ODEProblem{Matrix{Tracker.TrackedReal{Float64}}, Tuple{Float64, Float64}, true, Tuple{Tracker.TrackedReal{Float64}, TrackedArray{…,Matrix{Float64}}, TrackedArray{…,Matrix{Float64}}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Matrix{Tracker.TrackedReal{Float64}}, Nothing, Float64, Tuple{Tracker.TrackedReal{Float64}, TrackedArray{…,Matrix{Float64}}, TrackedArray{…,Matrix{Float64}}}, Float64, Tracker.TrackedReal{Float64}, Tracker.TrackedReal{Float64}, Float64, Vector{Matrix{Tracker.TrackedReal{Float64}}}, ODESolution{Tracker.TrackedReal{Float64}, 3, Vector{Matrix{Tracker.TrackedReal{Float64}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Tracker.TrackedReal{Float64}}}}, ODEProblem{Matrix{Tracker.TrackedReal{Float64}}, Tuple{Float64, Float64}, true, Tuple{Tracker.TrackedReal{Float64}, TrackedArray{…,Matrix{Float64}}, TrackedArray{…,Matrix{Float64}}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Matrix{Tracker.TrackedReal{Float64}}}, Vector{Float64}, Vector{Vector{Matrix{Tracker.TrackedReal{Float64}}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Tracker.TrackedReal{Float64}}, Matrix{Tracker.TrackedReal{Float64}}, Matrix{Tracker.TrackedReal{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Tracker.TrackedReal{Float64}, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{Matrix{Tracker.TrackedReal{Float64}}, Matrix{Tracker.TrackedReal{Float64}}, Matrix{Tracker.TrackedReal{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Tracker.TrackedReal{Float64}, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Tracker.TrackedReal{Float64}, Tracker.TrackedReal{Float64}, Tracker.TrackedReal{Float64}, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Float64, Tuple{}}, Matrix{Tracker.TrackedReal{Float64}}, Tracker.TrackedReal{Float64}, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq C:\Users\Brendan\.julia\packages\OrdinaryDiffEq\vfMzV\src\initdt.jl:130
  [3] auto_dt_reset!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Matrix{Tracker.TrackedReal{Float64}}, Nothing, Float64, Tuple{Tracker.TrackedReal{Float64}, TrackedArray{…,Matrix{Float64}}, TrackedArray{…,Matrix{Float64}}}, Float64, Tracker.TrackedReal{Float64}, Tracker.TrackedReal{Float64}, Float64, Vector{Matrix{Tracker.TrackedReal{Float64}}}, ODESolution{Tracker.TrackedReal{Float64}, 3, Vector{Matrix{Tracker.TrackedReal{Float64}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Tracker.TrackedReal{Float64}}}}, ODEProblem{Matrix{Tracker.TrackedReal{Float64}}, Tuple{Float64, Float64}, true, Tuple{Tracker.TrackedReal{Float64}, TrackedArray{…,Matrix{Float64}}, TrackedArray{…,Matrix{Float64}}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Matrix{Tracker.TrackedReal{Float64}}}, Vector{Float64}, Vector{Vector{Matrix{Tracker.TrackedReal{Float64}}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Tracker.TrackedReal{Float64}}, Matrix{Tracker.TrackedReal{Float64}}, Matrix{Tracker.TrackedReal{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Tracker.TrackedReal{Float64}, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fast_react_diffuse!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{Matrix{Tracker.TrackedReal{Float64}}, Matrix{Tracker.TrackedReal{Float64}}, Matrix{Tracker.TrackedReal{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Tracker.TrackedReal{Float64}, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Tracker.TrackedReal{Float64}, Tracker.TrackedReal{Float64}, Tracker.TrackedReal{Float64}, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Float64, Tuple{}}, Matrix{Tracker.TrackedReal{Float64}}, Tracker.TrackedReal{Float64}, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq C:\Users\Brendan\.julia\packages\OrdinaryDiffEq\vfMzV\src\integrators\integrator_interface.jl:437

In the above code, I am remaking the problem with a retyped u0. If I comment out that part, I get basically the same error.

d is 0.15833177844319388 (tracked)
sigma is 4.199186768619093 (tracked)
R_est eltype is Tracker.TrackedReal{Float64}
K_est eltype is Tracker.TrackedReal{Float64}
Type of R_est is TrackedArray{…,Matrix{Float64}}
Type of K_est is TrackedArray{…,Matrix{Float64}}
Type of p is Tuple{Tracker.TrackedReal{Float64}, TrackedArray{…,Matrix{Float64}}, TrackedArray{…,Matrix{Float64}}}
ERROR: MethodError: no method matching Float64(::Tracker.TrackedReal{Float64})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:772
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  ...
Stacktrace:
  [1] convert(#unused#::Type{Float64}, x::Tracker.TrackedReal{Float64})
    @ Base .\number.jl:7
  [2] setindex!
    @ .\array.jl:968 [inlined]
  [3] fast_react_diffuse!(du::Matrix{Float64}, u::Matrix{Float64}, p::Tuple{Tracker.TrackedReal{Float64}, TrackedArray{…,Matrix{Float64}}, TrackedArray{…,Matrix{Float64}}}, t::Float64)
    @ Main .\REPL[20]:5
  [4] (::SciMLBase.Void{typeof(fast_react_diffuse!)})(::Matrix{Float64}, ::Vararg{Any})
    @ SciMLBase C:\Users\Brendan\.julia\packages\SciMLBase\xWByK\src\utils.jl:392
  [5] (::FunctionWrappers.CallWrapper{Nothing})(f::SciMLBase.Void{typeof(fast_react_diffuse!)}, arg1::Matrix{Float64}, arg2::Matrix{Float64}, arg3::Tuple{Tracker.TrackedReal{Float64}, TrackedArray{…,Matrix{Float64}}, TrackedArray{…,Matrix{Float64}}}, arg4::Float64)
1 Like
9

If I try the same code without promoting u0 with the reversediff backend, it tries and fails to promote u0:

d is TrackedReal<6zH>(1.0, 0.0, LFm, ---)
sigma is TrackedReal<HWt>(1.0, 0.0, LFm, ---)
R_est eltype is ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}
K_est eltype is ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}
Type of R_est is ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}
Type of K_est is ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}
Type of p is Tuple{ReverseDiff.TrackedReal{Float64, Float64, Nothing}, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}
ERROR: StackOverflowError:
Stacktrace:
     [1] anyeltypedual(::Type{T}, counter::Any) where T<:Union{Set, AbstractArray} (repeats 2 times)
       @ DiffEqBase C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\forwarddiff.jl:124
     [2] MappingRF
       @ .\reduce.jl:95 [inlined]
     [3] _foldl_impl(op::Base.MappingRF{typeof(DiffEqBase.anyeltypedual), Base.BottomRF{typeof(DiffEqBase.promote_dual)}}, init::Type, itr::Core.SimpleVector)
       @ Base .\reduce.jl:62
     [4] foldl_impl
       @ .\reduce.jl:48 [inlined]
     [5] mapfoldl_impl
       @ .\reduce.jl:44 [inlined]
     [6] #mapfoldl#259
       @ .\reduce.jl:170 [inlined]
     [7] #mapreduce#263
       @ .\reduce.jl:302 [inlined]
     [8] __anyeltypedual(#unused#::Type{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}})
       @ DiffEqBase C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\forwarddiff.jl:118
     [9] anyeltypedual(::Type{T}, counter::Any) where T (repeats 2 times)
       @ DiffEqBase C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\forwarddiff.jl:121
--- the last 9 lines are repeated 5849 more times ---
 [52651] anyeltypedual(::Type{T}, counter::Any) where T<:Union{Set, AbstractArray} (repeats 2 times)
       @ DiffEqBase C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\forwarddiff.jl:124
 [52652] MappingRF
       @ .\reduce.jl:95 [inlined]
 [52653] _foldl_impl(op::Base.MappingRF{typeof(DiffEqBase.anyeltypedual), Base.BottomRF{typeof(DiffEqBase.promote_dual)}}, init::Type, itr::Core.SimpleVector)
       @ Base .\reduce.jl:62
 [52654] foldl_impl
       @ .\reduce.jl:48 [inlined]
 [52655] mapfoldl_impl
       @ .\reduce.jl:44 [inlined]
 [52656] #mapfoldl#259
       @ .\reduce.jl:170 [inlined]
 [52657] #mapreduce#263
       @ .\reduce.jl:302 [inlined]
 [52658] __anyeltypedual
       @ C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\forwarddiff.jl:118 [inlined]
 [52659] anyeltypedual (repeats 2 times)
       @ C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\forwarddiff.jl:121 [inlined]
10

I thought there is always a chance there is something wrong with my Julia installation, so I hopped to my laptop where I have done a little with this same project; deleted Julia, deleted the cache folders and installed Julia fresh (current stable release, v 1.8.2 with the Windows 64-bit installer). Also, the VSCode debugger had been giving me an error (I’m just doing everything in a REPL session in the terminal) so I thought maybe there was a layer of issues aside from the code itself.

After a completely fresh Julia install, I made a project directory, opened a workspace for it in VSCode, activated an environment, and added the following packages in this order:

  1. Turing
  2. DifferentialEquations
  3. StatsPlots
  4. LinearAlgebra
  5. ReverseDiff

I then ran the following code, from the Lotka-Volterra tutorial, using the debugger which had previously been giving me an error

using Turing
using DifferentialEquations

# Load StatsPlots for visualizations and diagnostics.
using StatsPlots
using LinearAlgebra

# Set a seed for reproducibility.
using Random
using ReverseDiff

# Turing.setadbackend(:forwarddiff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:reversediff)

Random.seed!(14);

# Define Lotka-Volterra model.
function lotka_volterra(du, u, p, t)
    # Model parameters.
    Ξ±, Ξ², Ξ³, Ξ΄ = p
    # Current state.
    x, y = u

    # Evaluate differential equations.
    du[1] = (Ξ± - Ξ² * y) * x # prey
    du[2] = (Ξ΄ * x - Ξ³) * y # predator

    return nothing
end

# Define initial-value problem.
u0 = [1.0, 1.0]
p = [1.5, 1.0, 3.0, 1.0]
tspan = (0.0, 10.0)
prob = ODEProblem(lotka_volterra, u0, tspan, p)

# Plot simulation.
plot(solve(prob, Tsit5()))

sol = solve(prob, Tsit5(); saveat=0.1)
odedata = Array(sol) + 0.8 * randn(size(Array(sol)))

# Plot simulation and noisy observations.
plot(sol; alpha=0.3)
scatter!(sol.t, odedata'; color=[1 2], label="")

@model function fitlv(data, prob)
    # Prior distributions.
    Οƒ ~ InverseGamma(2, 3)
    Ξ± ~ truncated(Normal(1.5, 0.5), 0.5, 2.5)
    Ξ² ~ truncated(Normal(1.2, 0.5), 0, 2)
    Ξ³ ~ truncated(Normal(3.0, 0.5), 1, 4)
    Ξ΄ ~ truncated(Normal(1.0, 0.5), 0, 2)

    # Simulate Lotka-Volterra model. 
    p = [Ξ±, Ξ², Ξ³, Ξ΄]
    predicted = solve(prob, Tsit5(); p=p, saveat=0.1)

    # Observations.
    for i in 1:length(predicted)
        data[:, i] ~ MvNormal(predicted[i], Οƒ^2 * I)
    end

    return nothing
end

model = fitlv(odedata, prob)

# Sample 3 independent chains with forward-mode automatic differentiation (the default).
chain = sample(model, NUTS(0.65), MCMCSerial(), 1000, 3; progress=false)

Unfortunately it gives me the same error as I have seen before. This is that it works fine, respects breakpoints, etc., until it goes to solve an ODE.

Exception has occurred: ErrorException
`llvmcall` must be compiled to be called

Stacktrace:
  [1] assume(v::Bool)
    @ FunctionWrappers C:\Users\brend\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:10
  [2] do_ccall(f::FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, args::Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64})
    @ FunctionWrappers C:\Users\brend\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:135
  [3] (::FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}})(args::Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64})
    @ FunctionWrappers C:\Users\brend\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:144
  [4] _call(fw::Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, arg::Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, fww::FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false})
    @ FunctionWrappersWrappers C:\Users\brend\.julia\packages\FunctionWrappersWrappers\YyXoO\src\FunctionWrappersWrappers.jl:12
  [5] (::FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false})(args::Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64})
    @ FunctionWrappersWrappers C:\Users\brend\.julia\packages\FunctionWrappersWrappers\YyXoO\src\FunctionWrappersWrappers.jl:10
  [6] (::ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing})(args::Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64})
    @ SciMLBase C:\Users\brend\.julia\packages\SciMLBase\xWByK\src\scimlfunctions.jl:1962
  [7] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, 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.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
    @ OrdinaryDiffEq C:\Users\brend\.julia\packages\OrdinaryDiffEq\vfMzV\src\perform_step\low_order_rk_perform_step.jl:7

Note that running the code in terminal gives no errors. This is the logical statement that triggers the llvmcall:

image

11

Setting debugger issues aside, the Lotka Volterra tutorial code (Bayesian Estimation of Differential Equations) behaves exactly the same on the fresh install as previously reported. Namely, it works fine with forward differentiation as well as reverse differentiation with Zygote/SciMLSensitivity, but using reversediff causes the following failure:

ERROR: ArgumentError: Converting an instance of ReverseDiff.TrackedReal{Float64, Float64, Nothing} to Float64 is not defined. Please use `ReverseDiff.value` instead.
Stacktrace:
  [1] convert(#unused#::Type{Float64}, t::ReverseDiff.TrackedReal{Float64, Float64, Nothing})
    @ ReverseDiff C:\Users\brend\.julia\packages\ReverseDiff\GtPeW\src\tracked.jl:261
  [2] setindex!(A::Vector{Float64}, x::ReverseDiff.TrackedReal{Float64, Float64, Nothing}, i1::Int64)
    @ Base .\array.jl:966
  [3] lotka_volterra(du::Vector{Float64}, u::Vector{Float64}, p::Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, t::Float64)
    @ Main c:\Dev\ecography_julia\src\test.jl:26
  [4] (::SciMLBase.Void{typeof(lotka_volterra)})(::Vector{Float64}, ::Vararg{Any})
    @ SciMLBase C:\Users\brend\.julia\packages\SciMLBase\xWByK\src\utils.jl:392
  [5] (::FunctionWrappers.CallWrapper{Nothing})(f::SciMLBase.Void{typeof(lotka_volterra)}, arg1::Vector{Float64}, arg2::Vector{Float64}, arg3::Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, arg4::Float64)
    @ FunctionWrappers C:\Users\brend\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:65
  [6] macro expansion
12

This is preliminary good news, but the base function with forward diff may actually be OK. Yes, learning the diffusion coefficient is easy if you see the diffusion process, but growth rates were high enough and carrying capacity low enough that most of the map was static most of the time (either 0 or at carrying capacity after the wavefront swept across).

With smaller sample sizes, there is a fair chance you would miss the wavefront. With larger sample sizes, you have enough information, but then it may be hard to find the exact parameter, leading to its own sort of instability. Even worse, I was sampling in a regular grid, making it harder to get lucky and observe any local interactions.

With that thought in mind, I modified the sampling scheme to remain a grid, but with a small neighborhood at each grid point instead of a single point. I then randomly dropped most points, leading to observation points like this:

image
image717Γ—477 11.3 KB

Since making that change, I’ve started one run which did not see d shoot towards infinity or bounce around failing to converge but settled after just a few minutes on 5 and is staying there for the moment. It also is looking like it will finish in less than an hour rather than the 6+ hours of the lone successful run I had before. This could be premature as it’s only partway into this one run, but it is certainly looking like even this very small amount of neighborhood information is enough to make the problem as straightforward as it should have been.

Problems remaining off the top of my head:

  1. Broken VS Code debugger
  2. Reversediff does not work for my code or tutorial code
  3. Tracker does not work for my code
  4. Zygote/SciMLSensitivity does not work for my code

Aka I have no ability to use reverse differentiation methods, and I have no debugger to help me solve it. Also still only crossing my fingers that I’ve solved the current problem.

13

I’ll help you out but I am travelling so I will not be able to respond every hour. Let’s look at a few here. @devmotion for reference may be able to help here as well.

I think the main issue here is the sampling issue, which I’ll leave to @devmotion and @mohamed82008.

That was just discovered five days ago by someone else.

The workaround is to just ensure FullSpecialize. prob = ODEProblem{true, FullSpecialize}(lotka_volterra, u0, tspan, p).

This would’ve been recently introduced by the compile time improvements (How Julia ODE Solve Compile Time Was Reduced From 30 Seconds to 0.1) so essentially the workaround is to turn off the no recompilation specializations since it sounds like the debugger doesn’t like using that pre-cached form. But that’s new information so sorry, we’ll make that work out nicer ASAP.

IIUC this issue was fixed by Fix and test `promote_rule` definitions by devmotion Β· Pull Request #207 Β· JuliaDiff/ReverseDiff.jl Β· GitHub . ReverseDiff version v1.14.3 shouldn’t show this error, which was released 2 days ago. Can you update and try that one?

For reference this was Missing Promotion Rule Β· Issue #206 Β· JuliaDiff/ReverseDiff.jl Β· GitHub.

This one is odd. As a workaround try:

predicted = vec(Array(solve(prob, Tsit5(); p=p, saveat=1.0, dt = 1e-6))[sample_idx])

i.e., just hardcode a initial dt and see if it works after that. @devmotion can we isolate this one?

What’s the error message in this one? I’m just scolling quickly here on my phone and don’t see it. There’s a lot of noise here. Without seeing anything, gut feel, try:

predicted = vec(Array(solve(prob, Tsit5(); p=p, saveat=1.0, autojacvec = ReverseDiffVJP(true)))[sample_idx])

which should work fine on your PDE form. Quick guess is Enzyme is throwing something.

14

Just a quick comment (I have to run to a meeting): MvNormal(predicted, sigma^2) and MvNormal(predicted, sigma^2 * I) are different distributions - the first one creates a MvNormal(predicted, sigma^4 * I) (I guess it was implemented in that way to be consistent with Normal(mu, sigma)). Due to the issues and confusions this causes the MvNormal(predicted, sigma^2) constructor is deprecated: Distributions.jl/mvnormal.jl at 4c27f7d4da28b1c9e033c699e8225b72283e6f29 Β· JuliaStats/Distributions.jl Β· GitHub The plan is to only allow to specify covariance matrices at some point.

Also, for debugging usually it is easier to only run a single chain.

15

I’m not able to reproduce your issues with the Lotka-Volterra model. The following seems to work with ReverseDiff without any issues:

using Turing
using ReverseDiff
using OrdinaryDiffEq
using SciMLSensitivity

using LinearAlgebra

# Define Lotka-Volterra model.
function lotka_volterra(du, u, p, t)
    # Model parameters.
    Ξ±, Ξ², Ξ³, Ξ΄ = p
    # Current state.
    x, y = u

    # Evaluate differential equations.
    du[1] = (Ξ± - Ξ² * y) * x # prey
    du[2] = (Ξ΄ * x - Ξ³) * y # predator

    return nothing
end

# Define initial-value problem.
u0 = [1.0, 1.0]
p = [1.5, 1.0, 3.0, 1.0]
tspan = (0.0, 10.0)
prob = ODEProblem(lotka_volterra, u0, tspan, p)
sol = solve(prob, Tsit5(); saveat=0.1)
odedata = Array(sol) + 0.8 * randn(size(Array(sol)))

@model function fitlv(data, prob)
    # Prior distributions.
    Οƒ ~ InverseGamma(2, 3)
    Ξ± ~ truncated(Normal(1.5, 0.5); lower=0.5, upper=2.5)
    Ξ² ~ truncated(Normal(1.2, 0.5); lower=0, upper=2)
    Ξ³ ~ truncated(Normal(3.0, 0.5); lower=1, upper=4)
    Ξ΄ ~ truncated(Normal(1.0, 0.5); lower=0, upper=2)

    # Simulate Lotka-Volterra model. 
    p = [Ξ±, Ξ², Ξ³, Ξ΄]
    predicted = solve(remake(prob; p=p), Tsit5(); saveat=0.1)

    # Observations.
    for i in 1:length(predicted)
        data[:, i] ~ MvNormal(predicted[i], Οƒ^2 * I)
    end

    return nothing
end

model = fitlv(odedata, prob)

ReverseDiff without tape compilation:

julia> sample(model, NUTS{Turing.ReverseDiffAD{false}}(), 1_000)
β”Œ Info: Found initial step size
β””   Ο΅ = 0.05
Sampling 100%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆ| Time: 0:01:52
Chains MCMC chain (1000Γ—17Γ—1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 112.47 seconds
Compute duration  = 112.47 seconds
parameters        = Οƒ, Ξ±, Ξ², Ξ³, Ξ΄
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec 
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64 

           Οƒ    0.7774    0.0384     0.0012    0.0016   418.1733    1.0005        3.7181
           Ξ±    1.5372    0.0550     0.0017    0.0035   213.4889    0.9995        1.8982
           Ξ²    1.0394    0.0514     0.0016    0.0026   298.0414    0.9998        2.6499
           Ξ³    2.8572    0.1438     0.0045    0.0092   223.1478    0.9999        1.9840
           Ξ΄    0.9743    0.0539     0.0017    0.0034   217.9203    0.9994        1.9376

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

           Οƒ    0.7082    0.7502    0.7758    0.8032    0.8565
           Ξ±    1.4317    1.5011    1.5357    1.5696    1.6565
           Ξ²    0.9453    1.0041    1.0349    1.0725    1.1484
           Ξ³    2.5756    2.7651    2.8554    2.9480    3.1481
           Ξ΄    0.8699    0.9403    0.9728    1.0075    1.0860

ReverseDiff with tape compilation:

julia> sample(model, NUTS{Turing.ReverseDiffAD{true}}(), 1_000)
β”Œ Info: Found initial step size
β””   Ο΅ = 0.025
β”Œ Warning: The current proposal will be rejected due to numerical error(s).
β”‚   isfinite.((ΞΈ, r, β„“Ο€, β„“ΞΊ)) = (true, false, false, false)
β”” @ AdvancedHMC /home/david/.julia/packages/AdvancedHMC/iWHPQ/src/hamiltonian.jl:47
β”Œ Warning: The current proposal will be rejected due to numerical error(s).
β”‚   isfinite.((ΞΈ, r, β„“Ο€, β„“ΞΊ)) = (true, false, false, false)
β”” @ AdvancedHMC /home/david/.julia/packages/AdvancedHMC/iWHPQ/src/hamiltonian.jl:47
Sampling 100%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆ| Time: 0:00:12
Chains MCMC chain (1000Γ—17Γ—1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 12.27 seconds
Compute duration  = 12.27 seconds
parameters        = Οƒ, Ξ±, Ξ², Ξ³, Ξ΄
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec 
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64 

           Οƒ    0.7761    0.0378     0.0012    0.0014   603.9262    0.9990       49.2157
           Ξ±    1.5334    0.0501     0.0016    0.0030   278.0967    1.0025       22.6629
           Ξ²    1.0392    0.0480     0.0015    0.0024   356.5059    1.0012       29.0527
           Ξ³    2.8667    0.1335     0.0042    0.0077   298.6526    1.0023       24.3381
           Ξ΄    0.9772    0.0495     0.0016    0.0028   286.3729    1.0025       23.3374

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

           Οƒ    0.7060    0.7507    0.7739    0.8000    0.8553
           Ξ±    1.4443    1.4973    1.5310    1.5662    1.6422
           Ξ²    0.9538    1.0046    1.0391    1.0708    1.1371
           Ξ³    2.6039    2.7700    2.8707    2.9614    3.1203
           Ξ΄    0.8778    0.9444    0.9774    1.0112    1.0712

I also checked ForwardDiff and Tracker:

ForwardDiff:

julia> sample(model, NUTS{Turing.ForwardDiffAD{0,true}}(), 1_000)
β”Œ Info: Found initial step size
β””   Ο΅ = 0.025
Sampling 100%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆ| Time: 0:00:03
Chains MCMC chain (1000Γ—17Γ—1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 3.76 seconds
Compute duration  = 3.76 seconds
parameters        = Οƒ, Ξ±, Ξ², Ξ³, Ξ΄
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec 
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64 

           Οƒ    0.7759    0.0366     0.0012    0.0017   577.8232    1.0000      153.6764
           Ξ±    1.5374    0.0546     0.0017    0.0037   215.5853    1.0013       57.3365
           Ξ²    1.0390    0.0479     0.0015    0.0028   273.7429    1.0014       72.8040
           Ξ³    2.8565    0.1469     0.0046    0.0095   221.5783    1.0013       58.9304
           Ξ΄    0.9737    0.0549     0.0017    0.0036   211.7973    1.0009       56.3291

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

           Οƒ    0.7077    0.7486    0.7748    0.8022    0.8471
           Ξ±    1.4278    1.5007    1.5387    1.5770    1.6371
           Ξ²    0.9481    1.0069    1.0339    1.0705    1.1399
           Ξ³    2.6096    2.7501    2.8488    2.9465    3.1688
           Ξ΄    0.8773    0.9333    0.9698    1.0081    1.0950

Tracker:

julia> sample(model, NUTS{Turing.TrackerAD}(), 1_000)
β”Œ Info: Found initial step size
β””   Ο΅ = 0.2
Sampling 100%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆ| Time: 0:04:52
Chains MCMC chain (1000Γ—17Γ—1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 292.64 seconds
Compute duration  = 292.64 seconds
parameters        = Οƒ, Ξ±, Ξ², Ξ³, Ξ΄
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec 
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64 

           Οƒ    0.7774    0.0383     0.0012    0.0017   530.6131    1.0005        1.8132
           Ξ±    1.5320    0.0530     0.0017    0.0033   236.0539    0.9992        0.8066
           Ξ²    1.0378    0.0483     0.0015    0.0025   321.7532    0.9990        1.0995
           Ξ³    2.8705    0.1455     0.0046    0.0095   233.1661    0.9992        0.7968
           Ξ΄    0.9795    0.0549     0.0017    0.0036   226.9989    0.9992        0.7757

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

           Οƒ    0.7041    0.7513    0.7765    0.8007    0.8538
           Ξ±    1.4272    1.4954    1.5353    1.5686    1.6292
           Ξ²    0.9472    1.0040    1.0356    1.0705    1.1367
           Ξ³    2.6177    2.7668    2.8597    2.9612    3.1762
           Ξ΄    0.8894    0.9424    0.9747    1.0135    1.0926

I used Julia 1.8.2 with the following packages:

(jl_pbezvu) pkg> st
Status `/tmp/jl_pbezvu/Project.toml`
  [1dea7af3] OrdinaryDiffEq v6.28.0
  [37e2e3b7] ReverseDiff v1.14.3
  [1ed8b502] SciMLSensitivity v7.10.0
  [fce5fe82] Turing v0.21.12

julia> versioninfo()
Julia Version 1.8.2
Commit 36034abf260 (2022-09-29 15:21 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 Γ— 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, tigerlake)
  Threads: 3 on 8 virtual cores
Environment:
  JULIA_NUM_THREADS = 3
16

My guess is because you’re running a version with https://github.com/JuliaDiff/ReverseDiff.jl/pull/207 . That damn issue caused a ton of AD problems, so I’m happy it’s finally merged. @breallis should update packages and see if that all goes away.

17

Yes, this fixed the debugger!

prob = ODEProblem{true, SciMLBase.FullSpecialize}(lotka_volterra, u0, tspan, p)

Yes, this fixed ReverseDiff for the tutorial code!

ReverseDiff version v1.14.3 shouldn’t show this error, which was released 2 days ago. Can you update and try that one?

18

And Tracker too?

19

Yes. I’m running nearly identical code on my laptop and desktop and I thought at first that the answer was no, as I got the following error message.

Exception has occurred: MethodError
MethodError: no method matching Float64(::Tracker.TrackedReal{Float64})
Closest candidates are:
  (::Type{T})(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:772
  (::Type{T})(!Matched::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  ...

Stacktrace:
  [1] iterate(g::Base.Generator{Base.Iterators.Zip{Tuple{UnitRange{Int64}, Vector{UInt64}}}, Base.var"#4#5"{AbstractMCMC.var"#sample_chain#78"{String, Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:chain_type, :progress), Tuple{UnionAll, Bool}}}, Random._GLOBAL_RNG, DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, Int64, Int64}}}, s::Tuple{})
    @ Base generator.jl:47
  [2] collect(itr::Base.Generator{Base.Iterators.Zip{Tuple{UnitRange{Int64}, Vector{UInt64}}}, Base.var"#4#5"{AbstractMCMC.var"#sample_chain#78"{String, Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:chain_type, :progress), Tuple{UnionAll, Bool}}}, Random._GLOBAL_RNG, DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, Int64, Int64}}})
    @ Base array.jl:787
  [3] map(f::AbstractMCMC.var"#sample_chain#78"{String, Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:chain_type, :progress), Tuple{UnionAll, Bool}}}, Random._GLOBAL_RNG, DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, Int64, Int64}, iters::Tuple{UnitRange{Int64}, Vector{UInt64}})
    @ Base abstractarray.jl:3055
  [4] mcmcsample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, ::MCMCSerial, N::Int64, nchains::Int64; progressname::String, init_params::Nothing, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:chain_type, :progress), Tuple{UnionAll, Bool}}})
    @ AbstractMCMC C:\Users\Brendan\.julia\packages\AbstractMCMC\fnRmh\src\sample.jl:518
  [5] (::AbstractMCMC.var"#mcmcsample##kw")(::NamedTuple{(:chain_type, :progress), Tuple{UnionAll, Bool}}, ::typeof(AbstractMCMC.mcmcsample), rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, ::MCMCSerial, N::Int64, nchains::Int64)
    @ AbstractMCMC C:\Users\Brendan\.julia\packages\AbstractMCMC\fnRmh\src\sample.jl:481
  [6] sample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, ensemble::MCMCSerial, N::Int64, n_chains::Int64; chain_type::Type{Chains}, progress::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Turing.Inference C:\Users\Brendan\.julia\packages\Turing\szPqN\src\inference\Inference.jl:220
  [7] (::StatsBase.var"#sample##kw")(::NamedTuple{(:progress,), Tuple{Bool}}, ::typeof(sample), rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Essential.TrackerAD, (), AdvancedHMC.DiagEuclideanMetric}}, ensemble::MCMCSerial, N::Int64, n_chains::Int64)

But the code was only almost the same. The key difference was flagged by devmotion:

MvNormal(predicted, sigma^2)

When I used the LinearAlgebra package and replaced it with sigma^2 * I the error went away.

I am excited to go back to my own code and see if this change combined with the updated ReverseDiff package will solve all the reversediff problems there as well. It also makes me wonder if making that change would have been sufficient to fix the problem even without the package update…

20

Unfortunately reversediff still does not work on my own code. There have been a lot of minor variations so for reference here is the code:

using Turing
using DifferentialEquations
using LinearAlgebra
using Random
using ReverseDiff

Turing.setadbackend(:reversediff)
#Turing.setrdcache(true)

Random.seed!(10);

const Nr = 50
const Nc = 100

true_R = rand(Beta(1,3), (Nr, Nc))
true_K =  true_R .* 50.0

true_d = 5.0
initial_pop = zeros(Nr, Nc)
initial_pop[1:20,end-30:end] = true_K[1:20,end-30:end]

tspan=(0.0, 25.0)
p = (true_d, true_R, true_K)

function fast_react_diffuse!(du,u,p,t)
  D, R, K = p

  @inbounds for j in 2:Nc-1, i in 2:Nr-1
    du[i,j] = D*(u[i-1,j] + u[i+1,j] + u[i,j+1] + u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for j in 2:Nc-1
    i = 1
    du[1,j] = D*(2u[i+1,j] + u[i,j+1] + u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for j in 2:Nc-1
    i = Nr
    du[end,j] = D*(2u[i-1,j] + u[i,j+1] + u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for i in 2:Nr-1
    j = 1
    du[i,1] = D*(u[i-1,j] + u[i+1,j] + 2u[i,j+1] - 4u[i,j])  +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds for i in 2:Nr-1
    j = Nc
    du[i,end] = D*(u[i-1,j] + u[i+1,j] + 2u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
  end

  @inbounds begin
    i = 1; j = 1
    du[1,1] = D*(2u[i+1,j] + 2u[i,j+1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])

    i = 1; j = Nc
    du[1,Nc] = D*(2u[i+1,j] + 2u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])

    i = Nr; j = 1
    du[Nr,1] = D*(2u[i-1,j] + 2u[i,j+1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])

    i = Nr; j = Nc
    du[end,end] = D*(2u[i-1,j] + 2u[i,j-1] - 4u[i,j]) +
    R[i,j]*u[i,j]*(1 - u[i,j]/K[i,j])
   end
end

prob = ODEProblem{true, SciMLBase.FullSpecialize}(fast_react_diffuse!,initial_pop,tspan,p)
sol = solve(prob,Tsit5(), saveat=1.0)
noisy_sol = Array(sol) + 0.4 * randn(size(Array(sol)))

@model function extra_simple_troubleshootfitlv(covariates, response, sample_idx, prob)
  # Prior distributions.
  d ~ Gamma(4, 1)
  sigma ~ InverseGamma(2, 3)
  println("d is ", d)
  println("sigma is ", sigma)
  p = (d, covariates[1], covariates[2])
  println("eltype R, K is ", eltype(p[2]), eltype(p[3]))
  predicted = vec(Array(solve(prob, Tsit5(); p=p, saveat=1.0))[sample_idx])
  # Observations.
  response ~ MvNormal(predicted, sigma^2 * I)
  return nothing
end

## 5x10x26 ##
sample_idx = CartesianIndices((5:10:size(noisy_sol)[1], 5:10:size(noisy_sol)[2], 1:size(noisy_sol)[3]))

ts_extra_simple_model_1 = extra_simple_troubleshootfitlv((p[2], p[3]), vec(noisy_sol[sample_idx]), sample_idx, prob)
chain_1 = sample(ts_extra_simple_model_1, NUTS(0.65), MCMCSerial(), 100, 1; progress=true)

And here is the error

d is 1.5472364971209513
sigma is 1.6552883057365462
eltype R, K is Float64Float64
d is 1.5472364971209513
sigma is 1.6552883057365462
eltype R, K is Float64Float64
d is TrackedReal<8l5>(1.0, 0.0, F6O, ---)
sigma is TrackedReal<DZQ>(1.0, 0.0, F6O, ---)
eltype R, K is Float64Float64
ERROR: ArgumentError: Converting an instance of ReverseDiff.TrackedReal{Float64, Float64, Nothing} to Float64 is not defined. Please use `ReverseDiff.value` instead.
Stacktrace:
  [1] convert(#unused#::Type{Float64}, t::ReverseDiff.TrackedReal{Float64, Float64, Nothing})
    @ ReverseDiff C:\Users\Brendan\.julia\packages\ReverseDiff\GtPeW\src\tracked.jl:261
  [2] setindex!
    @ .\array.jl:968 [inlined]
  [3] fast_react_diffuse!(du::Matrix{Float64}, u::Matrix{Float64}, p::Tuple{ReverseDiff.TrackedReal{Float64, Float64, Nothing}, Matrix{Float64}, Matrix{Float64}}, t::Float64)
    @ Main .\REPL[64]:5
  [4] (::ODEFunction{true, SciMLBase.FullSpecialize, typeof(fast_react_diffuse!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing})(::Matrix{Float64}, ::Vararg{Any})
    @ SciMLBase C:\Users\Brendan\.julia\packages\SciMLBase\xWByK\src\scimlfunctions.jl:1962
  [5] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Matrix{Float64}, Nothing, Float64, Tuple{ReverseDiff.TrackedReal{Float64, Float64, Nothing}, Matrix{Float64}, Matrix{Float64}}, Float64, Float64, Float64, Float64, Vector{Matrix{Float64}}, ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{ReverseDiff.TrackedReal{Float64, Float64, Nothing}, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fast_react_diffuse!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, typeof(fast_react_diffuse!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fast_react_diffuse!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, 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.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Float64, Tuple{}}, Matrix{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
    @ OrdinaryDiffEq C:\Users\Brendan\.julia\packages\OrdinaryDiffEq\vfMzV\src\perform_step\low_order_rk_perform_step.jl:736
  [6] __init(prob::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{ReverseDiff.TrackedReal{Float64, Float64, Nothing}, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fast_react_diffuse!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Float64, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::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 C:\Users\Brendan\.julia\packages\OrdinaryDiffEq\vfMzV\src\solve.jl:493
  [7] #__solve#562
    @ C:\Users\Brendan\.julia\packages\OrdinaryDiffEq\vfMzV\src\solve.jl:5 [inlined]
  [8] #solve_call#26
    @ C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\solve.jl:472 [inlined]
  [9] solve_up(prob::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fast_react_diffuse!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Matrix{Float64}, p::Tuple{ReverseDiff.TrackedReal{Float64, Float64, Nothing}, Matrix{Float64}, Matrix{Float64}}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ DiffEqBase C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\solve.jl:834
 [10] solve(prob::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fast_react_diffuse!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; sensealg::Nothing, u0::Nothing, p::Tuple{ReverseDiff.TrackedReal{Float64, Float64, Nothing}, Matrix{Float64}, Matrix{Float64}}, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ DiffEqBase C:\Users\Brendan\.julia\packages\DiffEqBase\5rKYk\src\solve.jl:801
 [11] extra_simple_troubleshootfitlv(__model__::DynamicPPL.Model{typeof(extra_simple_troubleshootfitlv), (:covariates, :response, :sample_idx, :prob), (), (), Tuple{Tuple{Matrix{Float64}, Matrix{Float64}}, Vector{Float64}, CartesianIndices{3, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}, UnitRange{Int64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fast_react_diffuse!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__::DynamicPPL.TypedVarInfo{NamedTuple{(:d, :sigma), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:d, Setfield.IdentityLens}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:d, Setfield.IdentityLens}}, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:sigma, Setfield.IdentityLens}}, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}}}, ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}}, __context__::DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Essential.ReverseDiffAD{true}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random._GLOBAL_RNG}, covariates::Tuple{Matrix{Float64}, Matrix{Float64}}, response::Vector{Float64}, sample_idx::CartesianIndices{3, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}, UnitRange{Int64}}}, prob::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fast_react_diffuse!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem})
    @ Main .\REPL[68]:9
 [12] macro expansion
    @ C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:593 [inlined]
 [13] _evaluate!!
    @ C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:576 [inlined]
 [14] evaluate_threadunsafe!!
    @ C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:551 [inlined]
 [15] evaluate!!
    @ C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:504 [inlined]
 [16] evaluate!!
    @ C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:515 [inlined]
 [17] evaluate!!
    @ C:\Users\Brendan\.julia\packages\DynamicPPL\zPOYL\src\model.jl:523 [inlined]
 [18] LogDensityFunction
    @ C:\Users\Brendan\.julia\packages\Turing\szPqN\src\Turing.jl:38 [inlined]
 [19] logdensity
...