Error with ForwardDiff in Turing

When running the sample method from Turing, I run into this error

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 9})
Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
(::Type{T})(::T) where T<:Number at boot.jl:772
(::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50

But I don’t really understand the error message. Is this saying there’s a problem with my arbitrary function within the Turing model?

Looks like a typing issue, which is common in these situations. This is probably due to your code, may we see a minimal reproducing example ?

With Turing errors, you generally have to ignore the actual error message itself, which is usually some obscure type error in the autodiff. The actual mistake is generally unrelated.

The best thing to do is look for the line that tells you where in your code the error happened. Looking at this line will usually reveal a mistake.

In the stack trace, I have the error

eq4(flag::Int64, params::Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 9}}, mhalo_tree::Vector{Float32}, z_tree::Vector{Float32})

where eq4 is my function. It looks to me as if its trying to pass a ForwardDiff.Dual object through my params argument when that argument takes in a vector of float. I’m not too sure how ForwardDiff works but I think the error message is complaining about the types being wrong here.

With regards to the minimal reproducing example, I’m not too sure how to make one because every line might have an offending mistake?

Edit: I have since worked out that I needed to type my params as AbstractVectors in order to take in the ForwardDiff.Dual object. However, this just leads to a new problem, which is that the code seems to be trying to be forcing the dual object vector into a Float. Stacktrace:

[1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 9})
    @ Base ./number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 9}, i1::Int64)
    @ Base ./array.jl:966
  [3] eq4(flag::Int64, params::Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 9}}, mhalo_tree::Vector{Float32}, z_tree::Vector{Float32})
    @ Main ./In[6]:115

However, when I look in line 115 of my eq4 code, it is this innocuous looking line of code sfr[i+1] = min(ζ * galaxy_inflow/(1+η), 0.02 * Mgas_predict/tdyn). Of course, ζ and η are calculated using various elements of the vector parameters, but nowhere should they interact with the full vector parameters. I don’t see the point where this is trying to force parameters into a Float.

What creates the vector sfr? This probably needs to change from e.g. zeros(10) to zeros(eltype(parameters), 10).

setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{... cannot work, as the array has a fixed space in memory with 64 bits per element, but the Dual number stores several numbers, so it can’t fit.

Thanks, that has sorted it. But the next issue that popped up is this:


 [1] cachedrule(#unused#::Type{ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 9}}, n::Int64) (repeats 79984 times)
   @ QuadGK ~/.julia/packages/QuadGK/buWps/src/gausskronrod.jl:253

Is QuadGK not compatible with ForwardDiff?

I’m not sure about the StackOverflowError but it seems like your program is trying to autodiff through a numerical integration routine. At best that’s super wasteful.
Others might know better, but a custom autodiff rule (that basically implements the Fundamental Theorem of Calculus for QuadGK) would probably be useful here.

Thanks. I had a look at the custom autodiff page, but I’m still very confused about how to make that work. Do you know if there are any tutorials about how to implement a custom rule for the Fundamental Theorem of Calculus in ForwardDiff?

I think the main issue here is that I’m not just simply differentiating with respect to the integral. It’s that my integral takes various dual numbers, and I need to allow the derivatives to propagate through also. There might be a way to use the fundamental theorem of calculus, but it’s not obvious to me.

If I’m not mistaken, GitHub - SciML/Integrals.jl: A common interface for quadrature and numerical integration for the SciML scientific machine learning organization already provides such rules on top of QuadGK functionality. See also Array issue with optimzing a function with quadgk for discussion on some other differentiable alternatives (N.B. that Quadrature.jl → Integrals.jl).

Thanks, I’ve used Integrals.jl, but unfortunately it is still showing an error. I have opened an issue on their GitHub and wondering if it can be resolved. However, I’ve also looked into custom rules (via ChainRules), and it seems that that might be more expedient, but I’m not too sure how to go about doing that. Thanks.

The forward and reverse mode rules are here:

Rewriting those won’t be more expedient :sweat_smile: . Instead, just identify why it’s not hit. My first guess is that you have parameters that are dual numbers but your state isn’t dual, because your bounds may not be set as dual numbers. In which case, you need to promote element types, and then the library should probably add an auto-promotion to help that. I haven’t dug into the code yet but that’s 90% of the cases with ForwardDiff.jl, and quite easy to fix with a simple eltype(p).(u0) kind of thing.

I see what you’re saying, since I previously had an integral, I wanted to make a custom rule to just take the original function and evaluate it as the derivative. But apparently ForwardDiff isn’t compatible with ChainRules…

By your eltype fix, I don’t quite understand what u0 is in this context. Is it my bounds, or is it all the numbers which go into my function?


prob = IntegralProblem(dNdr,rmin,rmax, [Mass, zmin, zmax])

try changing this to

prob = IntegralProblem(dNdr,convert(eltype(Mass),rmin),convert(eltype(Mass),rmax),
                                      [Mass, zmin, zmax])

See if the whole problem is dual if you’re fine.

After changing it to convert the types, the error still exists. I’ve even gone as far as converting zmin and zmax too to no avail.

I started suspecting that Mass is not actually a Dual, but Mass is an element of vector initialized using zeros(eltype(params), n_timebins) (where params is a Dual), which should force it to be a Dual right?

I don’t really know what next to try, or if some functionality is lacking from Integrals?

I also think if I were to do this manually, I could simply make dNdr the derivative and evaluate it at rmax and rmin to get the result of taking the derivative of intdNdr without going through Integrals, but I don’t see how this can be implemented with ForwardDiff since ForwardDiff doesn’t support ChainRules.

I still haven’t seen a code I can run, so it’s hard to say here. But what you described sounds like it would work with Integrals.jl today, and I showed the tests that prove it, so I’d like to hear what you aren’t describing.

Is there absolutely no version you can share that demonstrates the issue?

I have some progress trying to debug the code: specifically I have narrowed it to this section of code:

using Integrals
using ForwardDiff

function dNdr(r, p)
    Mass, zmin, zmax = p[1], p[2], p[3]
    α = 0.133
    β = -1.995
    γ = 0.263
    η = 0.0993
    A = 0.0104
    r̄ = 0.00972
    zfac = ((1+zmax)^(η+1) - (1+zmin)^(η+1))/(1+η)
    Nm = A * (Mass*1e-12)^α * r^β * exp((r/r̄)^γ) * zfac
    return Nm

function intdNdr(Mass, zmin, zmax, rmin, rmax)
    prob = IntegralProblem(dNdr,convert(eltype(Mass),rmin),convert(eltype(Mass),rmax), [Mass, convert(eltype(Mass),zmin), convert(eltype(Mass),zmax)])
    sol = solve(prob, QuadGKJL(), reltol=1e-4)
    return sol

This section is usually called by a different function which is then used in my Turing model, however, this function in isolation still yields the same error.

For example, if you try running the above code with:

Mass = 1e10 zmin=0.1 zmax=0.2 rmin=1e-8 rmax=1

ForwardDiff.derivative(Mass->intdNdr(Mass, zmin, zmax, rmin, rmax), Mass)

You get the same error message as you get from running my entire code, so I’m quite convinced that the error is somewhere in this function. I’m not clear that the Turing model is actually trying to differentiate w.r.t. Mass, but at least this looks like a starting point.

1 Like

That solves your issue.

1 Like

Thank you. Is this fix accessible via package manager or would it take some time for it to be pushed to the main branch?

It’ll need @stevengj to accept it first

I see that some checks have failed on GitHub…

Also, does this fix implement the fundamental theorem of calculus or does it still try to approximate f(x) somehow?

Actually, I think this issue also happens with any integrator in the Integrals package, so I don’t know if there are other fixes required to totally solve this.