# Simple 2D reaction-diffusion ODE inside Turing blows up in various ways
**URL:** https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177
**Category:** New to Julia
**Created:** [October 3, 2022, 6:04pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177 "2022-10-03T18:04:00Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 3, 2022, 6:04pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/1 "2022-10-03T18:04:01Z")
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/](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](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
The same happens with a smaller sample size (5x10x26)
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.
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
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.
```julia
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)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 3, 2022, 6:17pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/2 "2022-10-03T18:17:19Z")
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.
> [@breallis](#):
>
> 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.
>
> ```julia
>
> ```
It shouldn’t? Share the error with that.
---
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 3, 2022, 7:18pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/3 "2022-10-03T19:18:05Z")
Here is my Lotka-Volterra tutorial code using reversediff:
```julia
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:
```julia
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:
```julia
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:
```julia
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.
---
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 3, 2022, 7:34pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/4 "2022-10-03T19:34:12Z")
I can at least share some of the predicted and response.
These were the values drawn before things blew up
```julia
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)
```julia
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
```julia
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)
---
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 3, 2022, 7:50pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/5 "2022-10-03T19:50:46Z")
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:
```julia
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
```julia
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
```
---
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 3, 2022, 8:01pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/6 "2022-10-03T20:01:20Z")
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.
---
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 4, 2022, 2:18am UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/7 "2022-10-04T02:18:59Z")
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?
```julia
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).
---
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 4, 2022, 5:33am UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/8 "2022-10-04T05:33:47Z")
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:
```julia
@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:
```julia
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:
```julia
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.
```julia
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)
```
---
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 4, 2022, 5:40am UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/9 "2022-10-04T05:40:33Z")
If I try the same code without promoting u0 with the reversediff backend, it tries and fails to promote u0:
```julia
d is TrackedReal<6zH>(1.0, 0.0, LFm, ---)
sigma is TrackedReal(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]
```
---
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 4, 2022, 3:33pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/10 "2022-10-04T15:33:13Z")
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
```julia
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.
```julia
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:
---
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 4, 2022, 3:58pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/11 "2022-10-04T15:58:34Z")
Setting debugger issues aside, the Lotka Volterra tutorial code ([Bayesian Estimation of Differential Equations](https://turing.ml/dev/tutorials/10-bayesian-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:
```julia
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
```
---
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 5, 2022, 5:05am UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/12 "2022-10-05T05:05:52Z")
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:
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 5, 2022, 5:37am UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/13 "2022-10-05T05:37:10Z")
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.
> [@breallis](#):
>
> Broken VS Code debugger
That was just discovered five days ago by someone else.
>
>
> Summary - Running the following code through the debugger in VS Code fails with …
> \`\`Error Exception: \`llvmcall\` must be compiled to be called\`\`
>
> MWE -
> \`\`\`
> using DifferentialEquations
>
> function lorenz!(du,u,p,t)
> du\[1\] = 10.0\*(u\[2\]-u\[1\])
> du\[2\] = u\[1\]\*(28.0-u\[3\]) - u\[2\]
> du\[3\] = u\[1\]\*u\[2\] - (8/3)\*u\[3\]
> end
>
> function lorenz(u,p,t)
> du=similar(u)
> du\[1\] = 10.0\*(u\[2\]-u\[1\])
> du\[2\] = u\[1\]\*(28.0-u\[3\]) - u\[2\]
> du\[3\] = u\[1\]\*u\[2\] - (8/3)\*u\[3\]
> return du
> end
>
> u0 = \[1.0;0.0;0.0\]
> tspan = (0.0,100.0)
>
> \# Try using lorenz()
> prob = ODEProblem(lorenz,u0,tspan) # Works
> sol = solve(prob)
> println("Completed problem 1")
>
> \# Try using lorenz!()
> prob = ODEProblem(lorenz!,u0,tspan) # Exception
> sol = solve(prob)
> println("Completed problem 2")
> \`\`\`
>
> Details - This code solves \[example 2\](https://diffeq.sciml.ai/stable/tutorials/ode\_example/#Example-2:-Solving-Systems-of-Equations) from the documentation using both an RHS that allocates \`lorenz()\` and another that uses in-place updating \`lorenz!()\`.
>
> The in-place updating version fails when run with the debugger in VS Code with the following error
> \`\`\`
> Exception has occurred: ErrorException
> \`llvmcall\` must be compiled to be called
>
> Stacktrace:
> \[1\] assume(v::Bool)
> @ FunctionWrappers ~/.julia/packages/FunctionWrappers/8xdVB/src/FunctionWrappers.jl:10
> \[2\] do\_ccall(f::FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Float64}}, args::Tuple{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Float64})
> @ FunctionWrappers ~/.julia/packages/FunctionWrappers/8xdVB/src/FunctionWrappers.jl:128
> \[3\] (::FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Float64}})(args::Tuple{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Float64})
> @ FunctionWrappers ~/.julia/packages/FunctionWrappers/8xdVB/src/FunctionWrappers.jl:137
>
> \[see attached file for full stacktrace\]
> \`\`\`
> \[error.txt\](https://github.com/SciML/DifferentialEquations.jl/files/9685495/error.txt)
> \[VerisonInfo.txt\](https://github.com/SciML/DifferentialEquations.jl/files/9685760/VerisonInfo.txt)
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](https://sciml.ai/news/2022/09/21/compile_time/)) 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.
> [@breallis](#):
>
> ```julia
> 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
>
> ```
IIUC this issue was fixed by [Fix and test `promote_rule` definitions by devmotion · Pull Request #207 · JuliaDiff/ReverseDiff.jl · GitHub](https://github.com/JuliaDiff/ReverseDiff.jl/pull/207) . 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](https://github.com/JuliaDiff/ReverseDiff.jl/issues/206).
> [@breallis](#):
>
> ```julia
> 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
>
> ```
This one is odd. As a workaround try:
```julia
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?
> [@breallis](#):
>
> Zygote/SciMLSensitivity does not work for my code
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:
```julia
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.
---
### Author:  [@devmotion](https://discourse.julialang.org/u/devmotion)
#### Post date: [October 5, 2022, 7:48am UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/14 "2022-10-05T07:48:04Z")
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](https://github.com/JuliaStats/Distributions.jl/blob/4c27f7d4da28b1c9e033c699e8225b72283e6f29/src/multivariate/mvnormal.jl#L222) 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.
---
### Author:  [@devmotion](https://discourse.julialang.org/u/devmotion)
#### Post date: [October 5, 2022, 12:48pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/15 "2022-10-05T12:48:44Z")
I’m not able to reproduce your issues with the Lotka-Volterra model. The following seems to work with ReverseDiff without any issues:
```julia
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
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
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
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
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:
```julia
(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
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 5, 2022, 1:11pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/16 "2022-10-05T13:11:58Z")
> [@devmotion](#):
>
> I’m not able to reproduce your issues with the Lotka-Volterra model. The following seems to work with ReverseDiff without any issues:
>
> ```julia
>
> ```
My guess is because you’re running a version with [https://github.com/JuliaDiff/ReverseDiff.jl/pull/207](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.
---
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 5, 2022, 1:30pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/17 "2022-10-05T13:30:20Z")
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?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 5, 2022, 1:37pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/18 "2022-10-05T13:37:05Z")
And Tracker too?
---
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 5, 2022, 2:28pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/19 "2022-10-05T14:28:46Z")
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.
```julia
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:
```julia
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…
---
### Author:  [@breallis](https://discourse.julialang.org/u/breallis)
#### Post date: [October 5, 2022, 3:46pm UTC](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177/20 "2022-10-05T15:46:01Z")
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:
```julia
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
```julia
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(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
...
```
[Next page](https://discourse.julialang.org/t/simple-2d-reaction-diffusion-ode-inside-turing-blows-up-in-various-ways/88177.md?page=2)