Error with ForwardDiff in Turing

astro-kevin · October 12, 2022, 4:40pm

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?

lrnv · October 12, 2022, 4:51pm

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 ?

ParadaCarleton · October 12, 2022, 4:59pm

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.

astro-kevin · October 12, 2022, 9:28pm

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.

mcabbott · October 13, 2022, 1:51pm

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.

astro-kevin · October 13, 2022, 5:03pm

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

StackOverflowError:

Stacktrace:
 [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?

trahflow · October 15, 2022, 10:50am

Hi,
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.

astro-kevin · October 17, 2022, 1:12pm

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.

ToucheSir · October 17, 2022, 2:14pm

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).

astro-kevin · October 17, 2022, 2:21pm

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.

ChrisRackauckas · October 17, 2022, 2:48pm

The forward and reverse mode rules are here:

github.com

SciML/Integrals.jl/blob/master/src/Integrals.jl#L231-L391


      
          function ChainRulesCore.rrule(::typeof(__solvebp), prob, alg, sensealg, lb, ub, p, args...;
                                        kwargs...)
              out = __solvebp_call(prob, alg, sensealg, lb, ub, p, args...; kwargs...)
              function quadrature_adjoint(Δ)
                  y = typeof(Δ) <: Array{<:Number, 0} ? Δ[1] : Δ
                  if isinplace(prob)
                      dx = zeros(prob.nout)
                      _f = x -> prob.f(dx, x, p)
                      if sensealg.vjp isa ZygoteVJP
                          dfdp = function (dx, x, p)
                              _, back = Zygote.pullback(p) do p
                                  _dx = Zygote.Buffer(x, prob.nout, size(x, 2))
                                  prob.f(_dx, x, p)
                                  copy(_dx)
                              end
          
          
                    z = zeros(size(x, 2))
                              for idx in 1:size(x, 2)
                                  z[1] = 1
                                  dx[:, idx] = back(z)[1]

This file has been truncated. show original

Rewriting those won’t be more expedient . 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.

astro-kevin · October 17, 2022, 2:59pm

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?

Thanks

ChrisRackauckas · October 17, 2022, 3:41pm

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.

astro-kevin · October 17, 2022, 4:32pm

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.

ChrisRackauckas · October 17, 2022, 5:06pm

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?

astro-kevin · October 19, 2022, 3:31pm

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
end

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
end

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.

ChrisRackauckas · October 20, 2022, 3:03pm

github.com/JuliaMath/QuadGK.jl

Duals in the boundary values: avoid infinite recursion

JuliaMath:master ← ChrisRackauckas:dual

opened 03:02PM - 20 Oct 22 UTC

ChrisRackauckas

+5 -6

The previous behavior led to https://discourse.julialang.org/t/error-with-forwar…ddiff-in-turing/88633/16 . Dual numbers would just hit an infinite recursion because `float(::Dual)::Dual !<: AbstractFloat`. This should handle more cases where float(T) !<: AbstractFloat too, like TrackedReal. Though I know that this probably isn't the best way to handle it. In a higher level interface like Integrals.jl, I think it would be better to detect such differentiation by boundary terms and specialize it. It's actually kind of funny and would be a good homework problem. If you put a dual number into the boundary with dual part one, i.e. , so then you integrate (f + df/dt eps) dt and by linearity of the integral it's just fundamental theorem of calculus approximated by GK! However, the argument above assumes that the dual part generated in the rules is exactly 1 (or rather, matches the partials of the boundary term so that it weighs it correctly). The algorithm doesn't exactly work like that, but it ends up being mathematically the same. What happens naturally in the code is that the rule values are multiplied by segment sizes as weights`f(a + (1-x[2i-1])*s)` and it is there that the `x` is "imparted with" the correct dualness to then give the integral of the derivative. Of course, algorithmically it should probably just see this is the case, integrate f, and then slap f into the dual part. I'm not sure if QuadGK.jl should be mucking with dual numbers and autodiff overloads though, so I'll just do that at a higher level. But it's something to think about here. The reason to not do a complete overload is because for differentiating with respect to other quantities like parameters, you can fuse the primal and derivative calls rather than having them as two separate integrals, so it's definitely faster to differentiate the algorithm here. Except for handling the boundary. I don't know, it's pretty cool haha. Such a beautifully convoluted way to approximate f(x).

That solves your issue.

astro-kevin · October 20, 2022, 3:08pm

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

ChrisRackauckas · October 20, 2022, 3:13pm

It’ll need @stevengj to accept it first

astro-kevin · October 20, 2022, 5:02pm

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.

Thanks