TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 12}

Specific Domains Probabilistic programming
1

Dear all, I am to fit a model that uses an interaction surface in Turing.jl, the interaction surface is a function that is outside the model but gets called inside the model to estimate the competitive equivalence of individuals based on their size and density. It has a series of parameters that define that interaction, but for some reason, I am not able to make it work in the model. I was able to do it in stan, but I would really like to be able to do it in Turing. Below there is a minimal reproducible example:

# for the sake of time and learning I am performing the IPM analyses for the mesocosms in Julia 

# load the modules

using Chain
using DataFrames
using Plots
using Random
using StatsPlots
using Statistics
using LinearAlgebra
using Turing
using DynamicHMC
using Distributions
using MCMCChains



function α_surface!(target_tank::Int64, tank_var::Vector{Int64}, 
    spp::Vector{Int64}, z::Vector{Float64}, 
    ϕ_i, ϕ_j, η_ij, ρ_i, ρ_j, ι_ij, H, ω )
    sz = z[findall(tank_var .== target_tank)]# get the body size values per tank or channel
    sp = spp[findall(tank_var .== target_tank)] # get the species for each tank
    complength = zeros(length(sz),length(sz))

    for j in eachindex(sz)
        for k in eachindex(sz)

            if(sp[j] == 1 && sp[k] == 1)
                complength[j,k] = exp( ϕ_i * (sz[k] - ω) - ϕ_i * (sz[j] - ω)) * exp(-((ρ_i * (sz[k] - ω) - ρ_i * (sz[j] - ω))^2)/(2 * H^2))
            end
            if(sp[j] == 1 && sp[k] == 0)
                complength[j,k] = exp(η_ij + ϕ_j * (sz[k] - ω) - ϕ_i * (sz[j] - ω)) * exp(-((ι_ij + ρ_j * (sz[k] - ω) - ρ_i * (sz[j] - ω))^2/(2 * H^2)))
            end
            if(sp[j] == 0 && sp[k] == 0)
                complength[j,k] = (exp( ϕ_j * (sz[k] - ω) - ϕ_j * (sz[j] - ω)) )* exp(-((ρ_j * (sz[k] - ω) - ρ_j * (sz[j] - ω))^2)/(2 * H^2))
            end
            if(sp[j] == 0 && sp[k] == 1)
                complength[j,k] = exp(-η_ij + ϕ_j * (sz[k] - ω) - ϕ_i .* (sz[j] - ω))  * exp(-((-1*ι_ij + ρ_j * (sz[k] - ω) - ρ_i * (sz[j] - ω))^2)/(2 * H^2))
            end
        end
    end
    return TArray(sum(complength, dims= 2)[:,1])
end




# # function to get the density for all tanks or channels
function α_density(tank_var::Vector{Int64}, spp::Vector{Int64}, z::Vector{Float64}, 
    ϕ_i, ϕ_j, η_ij, ρ_i, ρ_j, ι_ij, H, ω )
    d1 = []
    for c in 1:maximum(tank_var)
        d2 = α_surface!(c,  tank_var, spp, z, ϕ_i, ϕ_j, η_ij, ρ_i, ρ_j, ι_ij, H, ω )
        d1 = append!(d1, d2)
    end

    return d1

end

df = DataFrame()

for i in 1:50

    append!(df, DataFrame(size = rand(Normal(32, 10), rand(5:10)), spp= 1, tank = i ))
    append!(df, DataFrame(size = rand(Normal(40, 10), rand(5:10)), spp= 0, tank = i ))

end

# simulate data

α_density(df.tank, df.spp, df.size, 0.05,0,0,0,0,0,1,10)
N = α_density(df.tank, df.spp, df.size, 0.02,0.02,0,0,0,0,1,10)
df.growth .=  0.397 .+  0.2*df.spp .-  0.762 .* log.(df.size ./ 10)  .+ 
    0.355.* log.(df.size.^2 ./ 10^2) .-   (0.073) .* N  .+ 0.084 .*N
df.growth .= exp.(rand.(Normal.(df.growth, 0.04)))

scatter(df.size, df.growth, group = df.spp)


α_density(df.tank, df.spp, df.size, 0.05,0,0,0,0,0,1,10)



@model function Model_Inter(growth, spp::Vector{Int64}, int_size::Vector{Float64}, ChanID, center)

    size = int_size 
    z = (log.(int_size ./ center))
    z2 = (log.(int_size.^2 ./ center^2))
    ## Total num of y
    N = length(growth)
    ## Separate σ priors for each actor
    σ_ch ~ truncated(Cauchy(0, 1), 0, 10)

    ## Number of unique mesocosms in the data set
    N_ch = length(unique(ChanID)) #7

    # --------------Survival model 

    ## Vector of mesocosms (1,..,7) which we'll set priors on
    α_ch ~ filldist(Normal(0, σ_ch), N_ch)
    #λ ~ TArray(     )#filldist(Normal(), N) # a latent variable

    ## Prior for intercept, prosoc_left, and the interaction
    α ~ Normal(0, 10)
    βz ~ Normal(0, 10)
    βz2 ~ Normal(0, 10)
    βspp ~ Normal(0, 10)
    βN ~ Normal(0, 10)
    βzN ~ Normal(0, 10)
    
    
    ϕ_i ~  Normal(0.495, 0.109)
    # ϕ_j ~  Normal(0.867, 0.102)
    # η_ij ~  Normal(-0.020, 0.010)
    # ρ_i ~  Normal(0.10, 0.021)
    # ρ_j ~  Normal( 0.03, 0.004)
    # ι_ij ~  Normal((3.71-3.12), 0.165)
    σ ~ truncated(Cauchy(0,1),0,10)
    H= 1.
    ω= center

    Den = tzeros(N)

    Den .= α_density(ChanID, spp, size, ϕ_i,0,0,0,0,0,H,ω)


    μ = α .+ α_ch[ChanID] .+  βspp .* spp .+ βz .* z .+  βz2 .* z2 .+ 
        βN .* Den .+ 
        βzN .* Den
   
    growth .~ MvNormal(μ, σ)


end;


chns = sample(
    Model_Inter(df.growth,df.spp, df.size, df.tank, 10.0 ),
    NUTS(),
    100
)

Thanks for the help

2

Without looking at your code too closely, the error is telling you that your code requires Float64 numbers, but a Dual number is passed, which is how automatic differentiation is implemented.

The likely solution is to remove your type annotations z::Vector{Float64}. Note that type annotations in Julia generally do not change performance (as methods are compiled for specific types in any case when a function is called), and are meant to be used for dispatch. It is therefore generally recommended to omit type annotations unless dispatch is needed, or at least use wider type signatures (e.g. Number instead of Float64) to allow your function to be generic enough to work with code other people have written - such as automatic differentiation libaries here.

3

just don’t annotate anything

4

Dear Nils, thank you so much for your quick reply. I corrected all the annotations, but unfortunately, your suggestion is not working either. Same error message.

5

I just did that, but it is not working either… it seems it is something related to the parameter ϕ_i. If I set this to 0 or other number it works, but they idea is precisely to estimate ϕ_i

6

Interesting - you’re saying you’re seeing the same TypeError after removing all type annotations from all function definitions?

In that case you’ll have to check for other implicit type restrictions in your code - I see e.g. in α_surface! you create an array by calling zeros and then write into it - zeros will create an Array{Float64} by default, so that should probably be zeros(eltype(ϕ_i)) to ensure it can accomodate duals.

7

where do you suggest adding the ‘zeros(eltype(ϕ_i))’ ? when I put it in the function outside? I get his error when trying that…

‘’’
no method matching exp(::Array{Float64, 0})
‘’’

8

I meant this line here:

complength = zeros(length(sz),length(sz))

which should probably be

complength = zeros(eltype(ϕ_i), length(sz),length(sz))
2 Likes
9

Thank you so much, it works now…!!!

1 Like