[Debugging] Need help on a ForwardDiff issue

New to Julia
1

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]
  [7] \(A::LinearAlgebra.LowerTriangular{Float64, LinearAlgebra.Adjoint{Float64, Matrix{Float64}}}, 
B::Vector{Real})
    @ LinearAlgebra C:\Users\lrnv\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\LinearAlgebra\src\triangular.jl:1660
  [8] invquad
    @ 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 :slight_smile:

2

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

3

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…

4

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
5

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.

6

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
7

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

1 Like
8

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}()

julia> ad = ForwardDiff.Dual{T}(1.0,0.0)
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)

julia> getllh(0.5,ad,bd,draws)
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
9

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
10

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?

11

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