Handling Dual types in my Turing model in a more elegant way

Hello everyone,

I’m doing inference on a system of ODEs solved with Rodas5. To use Rodas5 while writing inplace to my du, I check if one of u or du is a Dual so there’s no typeassert based errors. To use Turing with the No U-Turns Sampler, I do the same with all of my inference parameters, checking one by one if they are Duals, and if they are, making sure that all of my calculations are propagated properly using the types of my inference parameters. This seems to be a robust method as I have not had type errors since implementing the code below, but it also seems potentially inefficient and just plain ugly.

if eltype(u) <: ForwardDiff.Dual
    type_tester = u[1, 1]
elseif eltype(du) <: ForwardDiff.Dual
    type_tester = du[1, 1]
elseif isnothing(inference_parameters) # If not doing inference
    type_tester = 0.0
    # Scans inference_parameters for Duals,
    # returns the first or 0.0 if none are found.
    use_default = true
    for (type, parameter) in zip(
        if type <: ForwardDiff.Dual
            type_tester = parameter
            use_default = false
    if use_default
        type_tester = 0.0
# Used later in the function
current_eltype = eltype(type_tester)
F_total = get_tmp(F_total, type_tester)

What I would like to be able to write is:

values_to_test = (u[1, 1], d[1, 1], inference_parameters...)
if any(typeof.(values_to_test)) <: ForwardDiff.Dual
    # Because I want a specific type of Dual I think?
    type_tester = the_one_which_was_a_dual
    type_tester = 0.0

But I’m not sure if that’s either the best way of doing this, or even a good way of doing this. Is there something obvious I’m missing?

As a side question, in my Turing model if one of my inference parameters becomes a Dual, does that always mean that the others are? That would imply that I don’t need to test all of the inference parameters and just one.

Hi @jacobusmmsmit

a good way to avoid such conditional statements on types is to use parametric types instead.
As an example, if you write a function with one arument that is an abstract array of a parametric type like this:

function model(u::AbstractArray{T}) where T
   @show T

then T is the eltype of u inside your function:

julia> model([1., 2., 3.])
T = Float64

julia> model([1, 2, 3])
T = Int64

and you could work with T inside your function, instead of defining conditions for every possible eltype.

However for your task it might be better to write type-generic code in the first place, if possible.
For example to create a new array that has the same type (and size…) as another array, you can use u2 = similar(u).

For a more specific answer, could you maybe provide your original approach that leads to ForwardDiff.Dual-related type errors?

Share code. This shouldn’t be required.

1 Like

I tried to reproduce the error myself but it seems that there must have been a non-type-generic part of my code before which has now been replaced. I’m not able to reproduce the error.

Apologies for not doing this test before posting, it was just something on my long “bugs and gripes” list that I decided needed getting squared away and I thought it highly unlikely that removing an entire block of code and replacing it with type_tester = u would work first try.