# Why Turing.jl does not work with NUTS() and HMC() samplers?

**URL:** <https://discourse.julialang.org/t/why-turing-jl-does-not-work-with-nuts-and-hmc-samplers/107616>\
**Category:** General Usage\
**Tags:** turing\
**Created:** [December 14, 2023, 3:43pm UTC](https://discourse.julialang.org/t/why-turing-jl-does-not-work-with-nuts-and-hmc-samplers/107616 "2023-12-14T15:43:19Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![josimar](https://avatars.discourse-cdn.com/v4/letter/j/cab0a1/32.png) [@josimar](https://discourse.julialang.org/u/josimar)\
**Post date:** [December 14, 2023, 3:43pm UTC](https://discourse.julialang.org/t/why-turing-jl-does-not-work-with-nuts-and-hmc-samplers/107616/1 "2023-12-14T15:43:20Z")

</div>

I am experimenting with the `Turing.jl` package. I am trying to run the example model below but I get an error messagge when I try to use the `NUTS()` or the `HMC()` samplers. The error only occurs when I try to slice the array `t[tb:end]` inside the model. Note that the model works fine with the `MH()` sampler.

Would anyone be able to help me understand why the other samplers do not work and how do I fix this model? The exact error message that I get is shown below.

I am using `Turing.jl` version `v0.29.3`.

MWE here:

```julia
using Turing

@model function ModelTuring(y, N::Int, t)

    R=zeros(length(t))
    α ~ Distributions.Uniform(0, 1)

    tb=1
    ## This works fine with all samplers
    #R = α .+ t

    ## This only works with the MH() sampler only
    R[tb:end] = α .+ t[tb:end]

    for i=1:N
        y[i] ~ Poisson(R[i])
    end
    
end ##

N=10
t = collect(range(0; length=N, step=1))
Rdata= 0.5 .* t

model = ModelTuring(Rdata, N, t)
#sampler = NUTS()
sampler=HMC(0.1,10)
#sampler=MH()

nSamples=5
chain = sample(model, sampler, nSamples; init_params=[0.5])

```

Here is the first 50 lines of the error message:

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

Closest candidates are:
  (::Type{T})(::Real, !Matched::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:207
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  (::Type{T})(!Matched::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, 1})
    @ Base ./number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 1}, i1::Int64)
    @ Base ./array.jl:969
  [3] setindex!
    @ ./array.jl:983 [inlined]
  [4] ModelTuring( __model__ ::DynamicPPL.Model{typeof(ModelTuring), (:y, :N, :t), (), (), Tuple{Vector{Float64}, Int64, Vector{Int64}}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__ ::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.TypedVarInfo{NamedTuple{(:α,), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:α, Setfield.IdentityLens}, Int64}, Vector{Uniform{Float64}}, Vector{AbstractPPL.VarName{:α, Setfield.IdentityLens}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 1}}, Vector{Set{DynamicPPL.Selector}}}}}, ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 1}}, Vector{Base.RefValue{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 1}}}}, __context__ ::DynamicPPL.SamplingContext{DynamicPPL.Sampler{HMC{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.UnitEuclideanMetric}}, DynamicPPL.DefaultContext, Random.TaskLocalRNG}, y::Vector{Float64}, N::Int64, t::Vector{Int64})
    @ Main /lfs/nobackuplfs03/xom/seismicity/joadsjr/Documents/Work/Projects/Gorgon/ExtendedModelDomain/Abaqus/model_22/Julia/trash_24.jl:18
  [5] _evaluate!!(model::DynamicPPL.Model{typeof(ModelTuring), (:y, :N, :t), (), (), Tuple{Vector{Float64}, Int64, Vector{Int64}}, Tuple{}, DynamicPPL.DefaultContext}, varinfo::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.TypedVarInfo{NamedTuple{(:α,), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:α, Setfield.IdentityLens}, Int64}, Vector{Uniform{Float64}}, Vector{AbstractPPL.VarName{:α, Setfield.IdentityLens}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 1}}, Vector{Set{DynamicPPL.Selector}}}}}, ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 1}}, Vector{Base.RefValue{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 1}}}}, context::DynamicPPL.SamplingContext{DynamicPPL.Sampler{HMC{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.UnitEuclideanMetric}}, DynamicPPL.DefaultContext, Random.TaskLocalRNG})
    @ DynamicPPL ~/.julia/packages/DynamicPPL/oX6N7/src/model.jl:963
  [6] evaluate_threadsafe!!(model::DynamicPPL.Model{typeof(ModelTuring), (:y, :N, :t), (), (), Tuple{Vector{Float64}, Int64, Vector{Int64}}, Tuple{}, DynamicPPL.DefaultContext}, varinfo::DynamicPPL.TypedVarInfo{NamedTuple{(:α,), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:α, Setfield.IdentityLens}, Int64}, Vector{Uniform{Float64}}, Vector{AbstractPPL.VarName{:α, Setfield.IdentityLens}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 1}}, Vector{Set{DynamicPPL.Selector}}}}}, ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 1}}, context::DynamicPPL.SamplingContext{DynamicPPL.Sampler{HMC{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.UnitEuclideanMetric}}, DynamicPPL.DefaultContext, Random.TaskLocalRNG})
    @ DynamicPPL ~/.julia/packages/DynamicPPL/oX6N7/src/model.jl:952
  [7] evaluate!!
    @ ~/.julia/packages/DynamicPPL/oX6N7/src/model.jl:887 [inlined]
  [8] logdensity(f::LogDensityFunction{DynamicPPL.TypedVarInfo{NamedTuple{(:α,), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:α, Setfield.IdentityLens}, Int64}, Vector{Uniform{Float64}}, Vector{AbstractPPL.VarName{:α, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(ModelTuring), (:y, :N, :t), (), (), Tuple{Vector{Float64}, Int64, Vector{Int64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{HMC{Turing.Essential.ForwardDiffAD{0}, (), AdvancedHMC.UnitEuclideanMetric}}, DynamicPPL.DefaultContext, Random.TaskLocalRNG}}, θ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 1}})

```

---

<div class="post-metadata">

**Author:** ![cnrrobertson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cnrrobertson/32/205426_2.png) [@cnrrobertson](https://discourse.julialang.org/u/cnrrobertson)\
**Post date:** [December 14, 2023, 5:55pm UTC](https://discourse.julialang.org/t/why-turing-jl-does-not-work-with-nuts-and-hmc-samplers/107616/2 "2023-12-14T17:55:07Z")

</div>

This is common issue I’ve stumbled on when using `ForwardDiff.jl`. From what I can tell, when a `ForwardDiff.Dual` is added/subtracted/multiplied with an indexed vector of `Float64`, there is an attempted conversion to allow for the operation ([Conversion and Promotion · The Julia Language](https://docs.julialang.org/en/v1/manual/conversion-and-promotion/)). `ForwardDiff` deals with this gracefully when you use `R = alpha .+ t` but as soon as indexing is involved, it struggles. Someone with a better understanding than me will have to explain why that is. (I’m guessing it has to do with making temporary copies of the arrays during the operation).

My usual course of action is to try to remove indexed arrays in my code or to ensure that they take on the `ForwardDiff.Dual` element type. So, for example, in your code `t[tb:end]` and `R[tb:end]` have elements of types `Int64` and `Float64` respectively. If you do the following:

```julia
using Turing

@model function ModelTuring(y, N::Int, t)

    α ~ Distributions.Uniform(0, 1)

    #########################################################
    # Ensure arrays have the proper type for ForwardDiff.Dual
    R = zeros(eltype(α), length(t))
    t = convert(Array{eltype(α)}, t)
    #########################################################

    tb=1
    ## This works fine with all samplers
    # R = α .+ t

    ## This only works with the MH() sampler only
    R[tb:end] = α .+ t[tb:end]

    for i=1:N
        y[i] ~ Poisson(R[i])
    end

end ##

```

you can ensure that the elements of the arrays are correctly typed. There may be a more elegant/performant way, but I haven’t found it yet.

---

<div class="post-metadata">

**Author:** ![josimar](https://avatars.discourse-cdn.com/v4/letter/j/cab0a1/32.png) [@josimar](https://discourse.julialang.org/u/josimar)\
**Post date:** [December 14, 2023, 7:19pm UTC](https://discourse.julialang.org/t/why-turing-jl-does-not-work-with-nuts-and-hmc-samplers/107616/3 "2023-12-14T19:19:50Z")

</div>

Thank you very much for the help @cnrrobertson , you suggestion solved my problem.