# Error in nonmutating ODE function inside Turing (but mutating version not working beyond toy example)

**URL:** <https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244>\
**Category:** New to Julia\
**Tags:** turing, differentialequation\
**Created:** [September 25, 2023, 9:21pm UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244 "2023-09-25T21:21:23Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![caesoma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caesoma/32/219651_2.png) [@caesoma](https://discourse.julialang.org/u/caesoma)\
**Post date:** [September 25, 2023, 9:21pm UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/1 "2023-09-25T21:21:23Z")

</div>

Hi, all,

I I am trying to implement Turing HMC/NUTS inference of parameters from a model of coupled ODEs. Simulating that model works without issue, and I believe its details are not relevant (at a first pass, at least), but the function used to compute the differentials that go into creating the `DifferentialEquations.ODEProblem` object is of this form:

```julia
function coupled_ode_model(ρt, parameters, t)
    # ρt state of the system at the current time
    # dρ is the differential/difference computed at that same time

    ν, m, deltax, Wx, θ = parameters
    N = length(ρt)
    dρ = zeros(N)

    fluxes = coupling_function(ν, m, Wx, ρt, deltax; θ=θ)
    
    upperFlux = fluxes[2:end]
    lowerFlux = fluxes[1:end-1]
    
    dρ = - ( upperFlux - lowerFlux ) / deltax 

    return dρ

end

```

The syntax is reminiscent from ODE implementations in other languages and works as expected for simulating , but throws an error when used inside a Bayesian Turing model for inference. It’s pretty long and starts with:

```julia
DomainError with -0.3270308989755888:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

```

I looked up the [Turing ODE guide of a Lotka-Volterra example](https://turing.ml/dev/tutorials/10-bayesian-differential-equations/) and saw that it uses a different syntax where the differentials are an input and the function returns `nothing`, it doesn’t seem to work with the syntax with `return du` analog to the above.

So finally, the solutions would seem to be rewriting the differential function with the same syntax:

```julia
function coupled_ode_mutating(dρ, ρ, parameters, t)

    ν, m, deltax, Wx, θ = parameters
    
    fluxes = coupling_function(ν, m, Wx, ρt, deltax; θ=θ)
    
    upperFlux = fluxes[2:end]
    lowerFlux = fluxes[1:end-1]
    
    dρ = - ( upperFlux - lowerFlux ) / deltax 

    return nothing

end

```

but for some reason that is not working as the previous implementation and just repeats the initial state.

The main questions then are:

1. Why doesn’t Turing work with ODEs with the first syntax, and is there a reason for it not to?
2. Is there a clear reason why the second syntax is not working as expected?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 26, 2023, 1:50am UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/2 "2023-09-26T01:50:45Z")

</div>

This doesn’t have to do with Turing. You can recreate this with the same parameter choice and no Turing, did you try that?

You should make sure your definition is safe to not throw an error. Use NaNMath where you need to. Basically you’re just doing exactly what the error says, which should NaN. Did you print out the parameters for which this is occurring and isolate it? Check the parameters: are they sensible? Did you decrease the solver tolerance?

---

<div class="post-metadata">

**Author:** ![caesoma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caesoma/32/219651_2.png) [@caesoma](https://discourse.julialang.org/u/caesoma)\
**Post date:** [September 26, 2023, 8:08am UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/3 "2023-09-26T08:08:12Z")

</div>

Thank you for the quick reply, Chris. It does have at least something to do with Turing, though. I agree it can probably de reproduced be recreating the conditions under which the exception is being raised, but the fact of the matter those conditions are being created within Turing.

Even in the case of the Turing tutorial on differential equations mentioned above, simply replacing the “nothing syntax” with the “return du” syntax breaks the script (see `lotka_volterra_too` function below and full error in the post below, since it won’t fit into this one), although I cannot tell if the toy model fails for the same reason as the bigger model.

```julia
# Define Lotka-Volterra model.
function lotka_volterra_too(u, p, t)
    # Model parameters.
    α, β, γ, δ = p
    # Current state.
    x, y = u

    du = zeros(2)
    
    # Evaluate differential equations.
    du[1] = (α - β * y) * x # prey
    du[2] = (δ * x - γ) * y # predator

    return du
end

```

> [@ChrisRackauckas](#):
>
> You can recreate this with the same parameter choice and no Turing, did you try that?

I guess it’s possible, but if that is the case, either it happens for the initial set of parameter values before staring the chain, in which case I’d expect a Turing error that was informative of that, or it happens mid-chain, and the sampler should simply reject that sample and it would reflect on the diagnostics of the, but not halt inference altogether. If that is the most likely cause I can try different sets of parameters, or assigning an initial value that is supposed to work.

Nevertheless, I’m still left without answers for the two questions, especially the **issue #1** : in general is Turing supposed to work with the “return du” syntax (since it blatantly breaks the toy example), it will be helpful to know what works in principle before going down the rabbit hole of fixing a practically unfixable error.

> [@ChrisRackauckas](#):
>
> You should make sure your definition is safe to not throw an error. Use NaNMath where you need to. Basically you’re just doing exactly what the error says, which should NaN. Did you print out the parameters for which this is occurring and isolate it? Check the parameters: are they sensible? Did you decrease the solver tolerance?

That is all good advice, but it again doesn’t help understand what caused the specific error, and **issue #2** : is there an obvious way in which switching between the “nothing” and “return du” syntaxes should induce errors, or is bulletproofing the solver and hoping for the best the only way of tackling these kinds of issues?

Once again thanks for the recommendations. If you or anyone else have any additional comments, I’d really appreciate it. Meanwhile I will dig further into the error and post any additional information. Thanks.

---

<div class="post-metadata">

**Author:** ![caesoma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caesoma/32/219651_2.png) [@caesoma](https://discourse.julialang.org/u/caesoma)\
**Post date:** [September 26, 2023, 8:13am UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/4 "2023-09-26T08:13:38Z")

</div>

Here is the error mentioned above (i’ts too big for a single post, so it’s split between two posts)

```julia
MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3})

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:207
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number}
   @ Base char.jl:50
  ...

Stacktrace:
  [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3})
    @ Base ./number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}, i1::Int64)
    @ Base ./array.jl:969
  [3] lotka_volterra_too(u::Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, t::Float64)
    @ Main ./In[6]:11
  [4] (::ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing})(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, ::Vararg{Any})
    @ SciMLBase ~/.julia/packages/SciMLBase/YUCZN/src/scimlfunctions.jl:2126
  [5] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, Nothing, Float64, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, Float64, Float64, Float64, Float64, Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}}, ODESolution{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}, 2, Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}}}, ODEProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, Tuple{Float64, Float64}, false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}}, Vector{Float64}, Vector{Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}}}, OrdinaryDiffEq.Tsit5ConstantCache}, DiffEqBase.Stats, Nothing}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Tsit5ConstantCache, OrdinaryDiffEq.DEOptions{Float64, ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}, 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{}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5ConstantCache)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/87BwC/src/perform_step/low_order_rk_perform_step.jl:700
  [6] __init(prob::ODEProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, Tuple{Float64, Float64}, false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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 ~/.julia/packages/OrdinaryDiffEq/87BwC/src/solve.jl:499
  [7] __init (repeats 5 times)
    @ ~/.julia/packages/OrdinaryDiffEq/87BwC/src/solve.jl:10 [inlined]
  [8] __solve(::ODEProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, Tuple{Float64, Float64}, false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/87BwC/src/solve.jl:5
  [9] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/87BwC/src/solve.jl:1 [inlined]
 [10] #solve_call#22
    @ ~/.julia/packages/DiffEqBase/HoOGI/src/solve.jl:511 [inlined]
 [11] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HoOGI/src/solve.jl:972
 [12] solve_up
    @ ~/.julia/packages/DiffEqBase/HoOGI/src/solve.jl:945 [inlined]
 [13] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HoOGI/src/solve.jl:882
 [14] solve
    @ ~/.julia/packages/DiffEqBase/HoOGI/src/solve.jl:872 [inlined]
 [15] fitlv( __model__ ::DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__ ::DynamicPPL.TypedVarInfo{NamedTuple{(:σ, :α, :β), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:σ, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:σ, Setfield.IdentityLens}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:α, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, Vector{AbstractPPL.VarName{:α, Setfield.IdentityLens}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:β, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, Vector{AbstractPPL.VarName{:β, Setfield.IdentityLens}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, Vector{Set{DynamicPPL.Selector}}}}}, ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, __context__ ::DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, TaskLocalRNG}, data::Matrix{Float64}, prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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}, fixedparameters::Tuple{Float64, Float64})
    @ Main ./In[8]:13
 [16] _evaluate!!
    @ ~/.julia/packages/DynamicPPL/txq74/src/model.jl:963 [inlined]
 [17] evaluate_threadunsafe!!
    @ ~/.julia/packages/DynamicPPL/txq74/src/model.jl:936 [inlined]
 [18] evaluate!!(model::DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, varinfo::DynamicPPL.TypedVarInfo{NamedTuple{(:σ, :α, :β), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:σ, Setfield.IdentityLens}, Int64}, Vector{InverseGamma{Float64}}, Vector{AbstractPPL.VarName{:σ, Setfield.IdentityLens}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:α, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, Vector{AbstractPPL.VarName{:α, Setfield.IdentityLens}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:β, Setfield.IdentityLens}, Int64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, Vector{AbstractPPL.VarName{:β, Setfield.IdentityLens}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, Vector{Set{DynamicPPL.Selector}}}}}, ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}, context::DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, TaskLocalRNG})
    @ DynamicPPL ~/.julia/packages/DynamicPPL/txq74/src/model.jl:889
 [19] logdensity
    @ ~/.julia/packages/DynamicPPL/txq74/src/logdensityfunction.jl:94 [inlined]
 [20] Fix1
    @ ./operators.jl:1108 [inlined]
 [21] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/vXysl/src/apiutils.jl:24 [inlined]
 [22] vector_mode_gradient!(result::DiffResults.MutableDiffResult{1, Float64, Tuple{Vector{Float64}}}, f::Base.Fix1{typeof(LogDensityProblems.logdensity), 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, Float64, 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, Float64, Float64}}, Vector{AbstractPPL.VarName{:β, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, TaskLocalRNG}}}, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/gradient.jl:96
 [23] gradient!
    @ ~/.julia/packages/ForwardDiff/vXysl/src/gradient.jl:37 [inlined]
 [24] gradient!
    @ ~/.julia/packages/ForwardDiff/vXysl/src/gradient.jl:35 [inlined]
 [25] logdensity_and_gradient
    @ ~/.julia/packages/LogDensityProblemsAD/JoNjv/ext/LogDensityProblemsADForwardDiffExt.jl:113 [inlined]
 [26] ∂logπ∂θ
    @ ~/.julia/packages/Turing/LcjVL/src/mcmc/hmc.jl:160 [inlined]
 [27] ∂H∂θ
    @ ~/.julia/packages/AdvancedHMC/dgxuI/src/hamiltonian.jl:38 [inlined]

```

---

<div class="post-metadata">

**Author:** ![caesoma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caesoma/32/219651_2.png) [@caesoma](https://discourse.julialang.org/u/caesoma)\
**Post date:** [September 26, 2023, 8:14am UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/5 "2023-09-26T08:14:03Z")

</div>

```julia
 [28] phasepoint(h::AdvancedHMC.Hamiltonian{AdvancedHMC.DiagEuclideanMetric{Float64, Vector{Float64}}, AdvancedHMC.GaussianKinetic, Base.Fix1{typeof(LogDensityProblems.logdensity), LogDensityProblemsADForwardDiffExt.ForwardDiffLogDensity{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, Float64, 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, Float64, Float64}}, Vector{AbstractPPL.VarName{:β, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, TaskLocalRNG}}, ForwardDiff.Chunk{3}, ForwardDiff.Tag{Turing.TuringTag, Float64}, ForwardDiff.GradientConfig{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}}}}, Turing.Inference.var"#∂logπ∂θ#38"{LogDensityProblemsADForwardDiffExt.ForwardDiffLogDensity{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, Float64, 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, Float64, Float64}}, Vector{AbstractPPL.VarName{:β, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, TaskLocalRNG}}, ForwardDiff.Chunk{3}, ForwardDiff.Tag{Turing.TuringTag, Float64}, ForwardDiff.GradientConfig{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}}}}}, θ::Vector{Float64}, r::Vector{Float64})
    @ AdvancedHMC ~/.julia/packages/AdvancedHMC/dgxuI/src/hamiltonian.jl:80
 [29] phasepoint(rng::TaskLocalRNG, θ::Vector{Float64}, h::AdvancedHMC.Hamiltonian{AdvancedHMC.DiagEuclideanMetric{Float64, Vector{Float64}}, AdvancedHMC.GaussianKinetic, Base.Fix1{typeof(LogDensityProblems.logdensity), LogDensityProblemsADForwardDiffExt.ForwardDiffLogDensity{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, Float64, 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, Float64, Float64}}, Vector{AbstractPPL.VarName{:β, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, TaskLocalRNG}}, ForwardDiff.Chunk{3}, ForwardDiff.Tag{Turing.TuringTag, Float64}, ForwardDiff.GradientConfig{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}}}}, Turing.Inference.var"#∂logπ∂θ#38"{LogDensityProblemsADForwardDiffExt.ForwardDiffLogDensity{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, Float64, 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, Float64, Float64}}, Vector{AbstractPPL.VarName{:β, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, TaskLocalRNG}}, ForwardDiff.Chunk{3}, ForwardDiff.Tag{Turing.TuringTag, Float64}, ForwardDiff.GradientConfig{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}}}}}})
    @ AdvancedHMC ~/.julia/packages/AdvancedHMC/dgxuI/src/hamiltonian.jl:159
 [30] initialstep(rng::TaskLocalRNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, spl::DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), 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, Float64, 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, Float64, 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 ~/.julia/packages/Turing/LcjVL/src/mcmc/hmc.jl:164
 [31] step(rng::TaskLocalRNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, spl::DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}; resume_from::Nothing, init_params::Nothing, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}})
    @ DynamicPPL ~/.julia/packages/DynamicPPL/txq74/src/sampler.jl:111
 [32] step
    @ ~/.julia/packages/DynamicPPL/txq74/src/sampler.jl:84 [inlined]
 [33] macro expansion
    @ ~/.julia/packages/AbstractMCMC/fWWW0/src/sample.jl:125 [inlined]
 [34] macro expansion
    @ ~/.julia/packages/AbstractMCMC/fWWW0/src/logging.jl:16 [inlined]
 [35] mcmcsample(rng::TaskLocalRNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), 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 ~/.julia/packages/AbstractMCMC/fWWW0/src/sample.jl:116
 [36] mcmcsample
    @ ~/.julia/packages/AbstractMCMC/fWWW0/src/sample.jl:95 [inlined]
 [37] #sample#36
    @ ~/.julia/packages/Turing/LcjVL/src/mcmc/hmc.jl:121 [inlined]
 [38] sample
    @ ~/.julia/packages/Turing/LcjVL/src/mcmc/hmc.jl:91 [inlined]
 [39] (::AbstractMCMC.var"#sample_chain#78"{String, Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:chain_type, :progress), Tuple{UnionAll, Bool}}}, TaskLocalRNG, DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}, Int64, Int64})(i::Int64, seed::UInt64, init_params::Nothing)
    @ AbstractMCMC ~/.julia/packages/AbstractMCMC/fWWW0/src/sample.jl:511
 [40] sample_chain
    @ ~/.julia/packages/AbstractMCMC/fWWW0/src/sample.jl:508 [inlined]
 [41] #4
    @ ./generator.jl:36 [inlined]
 [42] iterate
    @ ./generator.jl:47 [inlined]
 [43] 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}}}, TaskLocalRNG, DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}, Int64, Int64}}})
    @ Base ./array.jl:782
 [44] map
    @ ./abstractarray.jl:3383 [inlined]
 [45] mcmcsample(rng::TaskLocalRNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), 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 ~/.julia/packages/AbstractMCMC/fWWW0/src/sample.jl:523
 [46] sample(rng::TaskLocalRNG, model::DynamicPPL.Model{typeof(fitlv), (:data, :prob, :fixedparameters), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(lotka_volterra_too), 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{Float64, Float64}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.DiagEuclideanMetric}}, ensemble::MCMCSerial, N::Int64, n_chains::Int64; chain_type::Type, progress::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Turing.Inference ~/.julia/packages/Turing/LcjVL/src/mcmc/Inference.jl:269
 [47] sample
    @ ~/.julia/packages/Turing/LcjVL/src/mcmc/Inference.jl:258 [inlined]
 [48] #sample#6
    @ ~/.julia/packages/Turing/LcjVL/src/mcmc/Inference.jl:254 [inlined]
 [49] sample
    @ ~/.julia/packages/Turing/LcjVL/src/mcmc/Inference.jl:245 [inlined]
 [50] #sample#5
    @ ~/.julia/packages/Turing/LcjVL/src/mcmc/Inference.jl:241 [inlined]

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 26, 2023, 1:57pm UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/6 "2023-09-26T13:57:21Z")

</div>

> [@caesoma](#):
>
> `du = zeros(2)`

Change that to `du = zeros(eltype(u), 2)`

---

<div class="post-metadata">

**Author:** ![dubosipsl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dubosipsl/32/52998_2.png) [@dubosipsl](https://discourse.julialang.org/u/dubosipsl)\
**Post date:** [September 26, 2023, 7:36pm UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/7 "2023-09-26T19:36:23Z")

</div>

Your mutating function is not mutating. You want :  
` dρ[:] = - ( upperFlux - lowerFlux )`

---

<div class="post-metadata">

**Author:** ![caesoma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caesoma/32/219651_2.png) [@caesoma](https://discourse.julialang.org/u/caesoma)\
**Post date:** [September 27, 2023, 12:03pm UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/8 "2023-09-27T12:03:23Z")

</div>

That works for the toy example, have yet to test with the parger model. I thought the zeros function would create a `Float64`-typed array, but in any case I guess I can assume both arrays will be `Float64` and do `du = zeros(Float64, 2)` with `function lotka_volterra_too(u::Float64, p, t)`.

> [@dubosipsl](#):
>
> Your mutating function is not mutating. You want :  
> ` dρ[:] = - ( upperFlux - lowerFlux )`

Thanks for that! I’m not used to this syntax for ODEs and didn’t know what it should look like, so I couldn’t catch that at all. I have yet to see if that makes a difference within the Turing function. If there’s something else needed to make it work that I find out I’'ll post it here for reference. Thanks to both.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 27, 2023, 12:09pm UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/9 "2023-09-27T12:09:56Z")

</div>

> [@caesoma](#):
>
> I thought the zeros function would create a `Float64`-typed array, but in any case I guess I can assume both arrays will be `Float64` and do `du = zeros(Float64, 2)` with `function lotka_volterra_too(u::Float64, p, t)`.

No you cannot, because they are not arrays of `Float64`.

---

<div class="post-metadata">

**Author:** ![caesoma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caesoma/32/219651_2.png) [@caesoma](https://discourse.julialang.org/u/caesoma)\
**Post date:** [September 27, 2023, 12:37pm UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/10 "2023-09-27T12:37:30Z")

</div>

What `eltype` is doing is finding out the type of `u` and making `du` of the same type, right? So if I assign a fixed suitable type for both it would have the same effect, or am I misunderstanding something?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 27, 2023, 12:39pm UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/11 "2023-09-27T12:39:38Z")

</div>

> [@caesoma](#):
>
> So if I assign a fixed suitable type for both it would have the same effect, or am I misunderstanding something?

Yes, but it’s not a Float64. The problem is that it’s a `Dual{Tag, Float64, N}` type of object in some calls. Do `@show eltype(u)` inside of your function and you will see. This is addressed in the FAQ.

---

<div class="post-metadata">

**Author:** ![caesoma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caesoma/32/219651_2.png) [@caesoma](https://discourse.julialang.org/u/caesoma)\
**Post date:** [September 27, 2023, 12:41pm UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/12 "2023-09-27T12:41:22Z")

</div>

FAQ for what?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 27, 2023, 12:48pm UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/13 "2023-09-27T12:48:05Z")

</div>

There’s discussion of the same thing in here:

[https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers](https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers)

---

<div class="post-metadata">

**Author:** ![caesoma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caesoma/32/219651_2.png) [@caesoma](https://discourse.julialang.org/u/caesoma)\
**Post date:** [September 27, 2023, 12:50pm UTC](https://discourse.julialang.org/t/error-in-nonmutating-ode-function-inside-turing-but-mutating-version-not-working-beyond-toy-example/104244/14 "2023-09-27T12:50:00Z")

</div>

Thanks, I’ll try to check it out.