I got an apparently quite common Julia error when trying to use AD with forward.diff. The error messages vary a bit (sometimes matching function name sometimes Float64)

MethodError: no method matching logL_multinom(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(logL_multinom), Real}, Real, 7}})

My goal: Transform a probability vector to be unbounded (θ -> y), do some stuff (namely HMC sampling) and transform back to the simplex space whenever the unnormalized posterior (logL_multinom()) is evaluated. DA should be used to overome problems for later, more complex, models than this.

Unfortunately, neither the Julia documentation, not the solutions from other questions helped me figure the particular problem out. Especially, it seems to work when I do the first transformation (y -> z) outside of the function, but the first transformation is a 1-to-1 mapping via logistic and should not cause any harm to differentiation.

Here is an MWE:

using LinearAlgebra
using ForwardDiff
using Base

function logL_multinom(y)
  # transform to constrained
  K = length(y)+1 
  k = collect(1:(K-1))
  # inverse logit: 
  z = 1 ./ (1 .+ exp.(-y .- log.(K .- k))) # if this is outside, it works
  θ = zeros(eltype(y),K) ; x_cumsum = zeros(eltype(y),K-1) 
  typeof(θ)
  for i in k
    x_cumsum[i] = 1-sum(θ)
    θ[i] = (x_cumsum[i]) * z[i]
  end
  θ[K] = x_cumsum[K-1] - θ[K-1]
  #log_dens_correction = sum( log(z*(1-z)*x_cumsum) )
  dot(colSums, log.(θ))
end
colSums = [835, 52, 1634, 3469, 3053, 2507, 2279, 1115]
y0 = [-0.8904013824298864, -0.8196709647741431, -0.2676845405543302, 0.31688184351556026, -0.870860684394019,0.15187821053559714,0.39888119498547964]
logL_multinom(y0)
∇L = y -> ForwardDiff.gradient(logL_multinom,y)
∇L(y0)

Thanks a lot and especially some further readings/ explanations for the problem are appreciated since I'll be working with it moreoften :D

Edit: I tried to convert the input and any intermediate variable into Real / arrays of these, but nothing helped so far.