[Debugging] Need help on a ForwardDiff issue

Hey,

In my package Copulas.jl there might be a typing error somewhere that does not allow ForwardDiff to pass through some of the functions. Symply trying to minimize the loglikelyhood of the model does not work as expected :

EDIT: You need `Copulas.jl#master` to see the bug for the moment.

``````using Copulas, Distributions, Random, ForwardDiff#, Turing, Plots, Cthulhu
Random.seed!(123)

# Definition of the model :
M₁ = Exponential(1.0)
M₂ = Exponential(1.0)
ρ = 0.5
Σ = [1.0 ρ; ρ 1.0]
C = GaussianCopula(Σ)
D = SklarDist(C, (M₁, M₂))
draws = rand(D, 1_000)

# This function works correctly, but is not ForwardDiff-able with respect to all the parameters of the model !
function getllh(ρ,a,b,draws)
D = SklarDist(GaussianCopula([1.0 ρ; ρ 1.0]),(Exponential(a),Exponential(b)))
return loglikelihood(D,draws)
end
getllh(0.5,1,1,draws) # OK
getllh(0.9,1,1,draws) # OK
getllh(0.5,3,2,draws) # OK
getllh(0.5,4,7,draws) # OK, Looks coherent.

ForwardDiff.derivative(ρ -> getllh(ρ,1,1,draws),0.1) # rho is too small, should be positive
ForwardDiff.derivative(ρ -> getllh(ρ,1,1,draws),0.9) # rho is too big, should be negative.

ForwardDiff.derivative(a -> getllh(0.5,1,a,draws),1) # Errors.
``````

Which gives the following stacktrace :

``````julia> ForwardDiff.derivative(a -> getllh(0.5,1,a,draws),1) # Errors.
ERROR: TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 1}
Stacktrace:
[1] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Int64}, Float64, 1}, i1::Int64)
@ Base .\array.jl:903
[2] _unsafe_copyto!(dest::Vector{Float64}, doffs::Int64, src::Vector{Real}, soffs::Int64, n::Int64)
@ Base .\array.jl:253
[3] unsafe_copyto!
@ .\array.jl:307 [inlined]
[4] _copyto_impl!
@ .\array.jl:331 [inlined]
[5] copyto!
@ .\array.jl:317 [inlined]
[6] copyto!
@ .\array.jl:343 [inlined]
B::Vector{Real})
@ LinearAlgebra C:\Users\lrnv\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\LinearAlgebra\src\triangular.jl:1660
@ C:\Users\lrnv\.julia\packages\PDMats\YXkIE\src\generics.jl:104 [inlined]
[9] sqmahal(d::ZeroMeanFullNormal{Tuple{Base.OneTo{Int64}}}, x::Vector{Real})
@ Distributions C:\Users\lrnv\.julia\packages\Distributions\O5xl5\src\multivariate\mvnormal.jl:267
[10] _logpdf
@ C:\Users\lrnv\.julia\packages\Distributions\O5xl5\src\multivariate\mvnormal.jl:143 [inlined]
[11] logpdf(d::ZeroMeanFullNormal{Tuple{Base.OneTo{Int64}}}, x::Vector{Real})
@ Distributions C:\Users\lrnv\.julia\packages\Distributions\O5xl5\src\common.jl:226
[12] _logpdf(C::GaussianCopula{2, Matrix{Float64}}, u::Vector{Real})
@ Copulas C:\Users\lrnv\.julia\dev\Copulas\src\EllipticalCopula.jl:15
[13] logpdf(d::GaussianCopula{2, Matrix{Float64}}, x::Vector{Real})
@ Distributions C:\Users\lrnv\.julia\packages\Distributions\O5xl5\src\common.jl:226
[14] _logpdf(S::SklarDist{GaussianCopula{2, Matrix{Float64}}, Tuple{Exponential{Float64}, Exponential{ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Int64}, Int64, 1}}}}, u::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true})
@ Copulas C:\Users\lrnv\.julia\dev\Copulas\src\SklarDist.jl:27
[15] logpdf(d::SklarDist{GaussianCopula{2, Matrix{Float64}}, Tuple{Exponential{Float64}, Exponential{ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Int64}, Int64, 1}}}}, x::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true})
@ Distributions C:\Users\lrnv\.julia\packages\Distributions\O5xl5\src\common.jl:226
[16] (::Base.Fix1{typeof(logpdf), SklarDist{GaussianCopula{2, Matrix{Float64}}, Tuple{Exponential{Float64}, Exponential{ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Int64}, Int64, 1}}}}})(y::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true})
@ Base .\operators.jl:1136
[17] (::Base.MappingRF{Base.Fix1{typeof(logpdf), SklarDist{GaussianCopula{2, Matrix{Float64}}, Tuple{Exponential{Float64}, Exponential{ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Int64}, Int64, 1}}}}}, Base.BottomRF{typeof(Base.add_sum)}})(acc::Base._InitialValue, x::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true})
@ Base .\reduce.jl:95
[18] _foldl_impl(op::Base.MappingRF{Base.Fix1{typeof(logpdf), SklarDist{GaussianCopula{2, Matrix{Float64}}, Tuple{Exponential{Float64}, Exponential{ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Int64}, Int64, 1}}}}}, Base.BottomRF{typeof(Base.add_sum)}}, init::Base._InitialValue, itr::Distributions.EachVariate{1, Matrix{Float64}, Tuple{Base.OneTo{Int64}}, SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}, 1})
@ Base .\reduce.jl:58
[19] foldl_impl(op::Base.MappingRF{Base.Fix1{typeof(logpdf), SklarDist{GaussianCopula{2, Matrix{Float64}}, Tuple{Exponential{Float64}, Exponential{ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Int64}, Int64, 1}}}}}, Base.BottomRF{typeof(Base.add_sum)}}, nt::Base._InitialValue, itr::Distributions.EachVariate{1, Matrix{Float64}, Tuple{Base.OneTo{Int64}}, SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}, 1})
@ Base .\reduce.jl:48
[20] mapfoldl_impl
@ .\reduce.jl:44 [inlined]
[21] #mapfoldl#244
@ .\reduce.jl:162 [inlined]
[22] mapfoldl
@ .\reduce.jl:162 [inlined]
[23] _mapreduce
@ .\reduce.jl:423 [inlined]
[24] _mapreduce_dim
@ .\reducedim.jl:330 [inlined]
[25] #mapreduce#725
@ .\reducedim.jl:322 [inlined]
[26] mapreduce
@ .\reducedim.jl:322 [inlined]
[27] #_sum#735
@ .\reducedim.jl:894 [inlined]
[28] _sum
@ .\reducedim.jl:894 [inlined]
[29] #sum#733
@ .\reducedim.jl:890 [inlined]
[30] sum
@ .\reducedim.jl:890 [inlined]
[31] loglikelihood(d::SklarDist{GaussianCopula{2, Matrix{Float64}}, Tuple{Exponential{Float64}, Exponential{ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Int64}, Int64, 1}}}}, x::Matrix{Float64})
@ Distributions C:\Users\lrnv\.julia\packages\Distributions\O5xl5\src\common.jl:440
[32] getllh(ρ::Float64, a::Int64, b::ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Int64}, Int64, 1}, draws::Matrix{Float64})
@ Main c:\Users\lrnv\Desktop\julia testing\test2.jl:16
[33] (::var"#9#10")(a::ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Int64}, Int64, 1})
@ Main c:\Users\lrnv\Desktop\julia testing\test2.jl:26
[34] derivative(f::var"#9#10", x::Int64)
@ ForwardDiff C:\Users\lrnv\.julia\packages\ForwardDiff\wAaVJ\src\derivative.jl:14
[35] top-level scope
@ c:\Users\lrnv\Desktop\julia testing\test2.jl:26
julia>
``````

I tried Cthulhu’s `@descend` to try understanding where in my code the mistake stays, but i could not understand where it is.

Could someone help me debug this ? I suspect this is my fault and therefore could be inside `Copulas.jl`’s code, but I cannot find the mistake.

Thanks

can you try with this?:

``````ForwardDiff.derivative(a -> getllh(0.5,1,a,draws),1.0) #Float number here
``````

maybe `Exponential(::Int)` dispatches to an specific implementation

Unfortunately, it gives exactly the same stack… I have edited the MWE to remove unnecessary packages so its easier to run. I still cannot find the problem…

Does `GaussianCopula([1.0 ρ; ρ 1.0])` create a Matrix of Duals or Floats when ForwardDiff is called? Otherwise you might need to assure that the Matrix has only elements of `typeof(ρ)`. All packages that you used (except Copulas, which I did not use yet) should work with ForwardDiff, so somewhere in the Code you dont have parametric types yet for ForwardDiff to work.

1 Like

Yes, @mrVeng, you are probably right : the problem is that my parametrization of the types inside `Copulas.jl` is not setup correctly, but I cannot find out in which places the problem lies.

What are you saying the code must do ? Which call should produce a matrix of duals ? `GaussianCopula(...)` produces a distribution, not a matrix, hence i have troubles follwing what you mean.

Is this on `Copulas.jl#master`? Because the MWE errors for me in a fresh environment on release `0.1.2` with

``````ERROR: MethodError: no method matching _logpdf(::Exponential{Float64}, ::Float64)`
``````

due to this call

``````function Distributions._logpdf(S::SklarDist{CT,TplMargins},u) where {CT,TplMargins}
sum(Distributions._logpdf.(S.m[i],u[i]) for i in 1:length(u)) + Distributions._logpdf(S.C,Distributions.cdf.(S.m,u))
#          ^^^^^^^^^^^^^ error
end
``````

Is there a reason you calling `_logpdf` and not simply `logpdf`? The former is not defined for every distribution, and the underscore marks it as ‘private’ by convention.

1 Like

Now it is on `Copulas.jl#master`, I forgot to push the `_logpdf``logpdf` change yesterday before releasing `0.1.2`

1 Like

Here is the relevant part of debugging. The `>` is the call that ultimately throws the error.

``````In \(A, B) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/triangular.jl:1656
1656  @eval function \(A::Union{UpperTriangular,LowerTriangular}, B::\$mat)
1657      require_one_based_indexing(B)
1658      TAB = typeof((zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))/one(eltype(A)))
1659      BB = similar(B, TAB, size(B))
>1660      copyto!(BB, B)
1661      ldiv!(convert(AbstractArray{TAB}, A), BB)
1662  end

About to run: <(copyto!)([-0.6457306721039766, 2.308264642e-314], Real[-0.6457306721039766, Dual{Float64}(-1.5900782...>
8|debug> fr
[8] \(A, B) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/triangular.jl:1656
| A::LinearAlgebra.LowerTriangular{Float64, LinearAlgebra.Adjoint{Float64, Matrix{Float64}}} = [1.0 0.0; 0.5 0.8660254037844386]
| B::Vector{Real} = Real[-0.6457306721039766, Dual{Float64}(-1.5900782022924036,-0.481991677966545)]
| BB::Vector{Float64} = [-0.6457306721039766, 2.308264642e-314]
| TAB::DataType = Float64
``````

The problem is that `B` is only a vector of `Real` and not `Dual` and conversion to `TAB = Float64` is attempted.

I could ‘cure’ it manually by doing

``````julia> T = ForwardDiff.Tag(getllh, Float64)
ForwardDiff.Tag{typeof(getllh), Float64}()

Dual{ForwardDiff.Tag{typeof(getllh), Float64}()}(1.0,0.0)

julia> bd = ForwardDiff.Dual{T}(1.0,1.0)
Dual{ForwardDiff.Tag{typeof(getllh), Float64}()}(1.0,1.0)

Dual{ForwardDiff.Tag{typeof(getllh), Float64}()}(-1850.7331612741796,-56.25528540379574)
``````

Edit: Or by taking a gradient

``````julia> ForwardDiff.gradient(x -> getllh(0.5,x[1],x[2],draws),[1.0,1.0]) # Errors.
2-element Vector{Float64}:
54.86694187385668
-56.25528540379574
``````

I guess the two entries in `B` are the parameters of the two exponential distributions. At some point they get bunched together and I would suggest you try and promote them to a common concrete type then.

1 Like

Thanks a lot for the find. So the marginal distributions of the `SklarDist` object should probably have the same eltype.

However, how do I promote the eltype of a `Distributions.UnivariateDistribution` without knowing which one it is ?..

EDIT: The simple way to do it for this application is :

``````function getllh(ρ,a,b,draws)
a,b = promote(a,b)
D = SklarDist(GaussianCopula([1.0 ρ; ρ 1.0]),(Exponential(a),Exponential(b)))
return loglikelihood(D,draws)
end
``````

and same thing will probably work for Turing.

Thanks a lot @skleinbo for finding that. Do you think there is still a problem with the parametric typing in my package or does it look ok to you ?

1 Like

That’s probably not so easy, but I haven’t really dug into your side of the code.

Instead of the workaround with `promote`, why not take the gradient? Chances are, you need it anyway when doing optimization.

Do you think there is still a problem with the parametric typing in my package or does it look ok to you ?

I haven’t been looking at that, sorry. But that feels like a separate issue. Open another thread?

Yes, I think that Turing is actually taking the gradient, which works correctly as it is now.

If something comes up as an error somewhere another time, I will, but at the moment everything looks fine: we were trying to debug Turing’s MCMC, which were due to other bugs, and we tries by hand `ForwarDiff.derivative`, which yielded this one : but you are completely right, we really need gradients, for which the typing stack looks OK.

Thanks again for the help !

1 Like