ReverseDiff is outright wrong when differentiating through VecCorrBijector

[penelopeysm]

using Bijectors
using Bijectors.VectorBijectors
import ForwardDiff, ReverseDiff, ChainRules, FiniteDifferences

d = LKJ(3, 1.0)
x = rand(d)

xv = vec(x)
f(xv) = Bijectors.VecCorrBijector()(reshape(xv, 3, 3))

f(xv)

ForwardDiff.jacobian(f, xv)
FiniteDifferences.jacobian(FiniteDifferences.central_fdm(5, 1), f, xv)[1]
ReverseDiff.jacobian(f, xv)

julia> ForwardDiff.jacobian(f, xv)
3×9 Matrix{Float64}:
 -0.016369  0.0  0.0  1.00107   0.0       0.0   0.0      0.0       0.0
  0.974712  0.0  0.0  0.0       0.0       0.0   2.51252  0.0       0.0
 -1.06427   0.0  0.0  2.37804  -0.653379  0.0  -2.64313  3.00987  -1.63988

julia> FiniteDifferences.jacobian(FiniteDifferences.central_fdm(5, 1), f, xv)[1]
3×9 Matrix{Float64}:
 -0.016369  -1.40702e-16  -1.40702e-16  1.00107      -2.52589e-15  -1.40702e-16  -1.12903e-14  -6.92803e-15  -6.47209e-15
  0.974712   4.50246e-15   4.50246e-15  1.91871e-13   8.08286e-14   4.50246e-15   2.51252       2.21697e-13   2.07107e-13
 -1.06427   -2.25123e-15  -2.25123e-15  2.37804      -0.653379     -2.25123e-15  -2.64313       3.00987      -1.63988

julia> ReverseDiff.jacobian(f, xv)
3×9 Matrix{Float64}:
  0.0      0.0  0.0  0.0       0.0       0.0   0.0      0.0       0.0
  0.0      0.0  0.0  0.0       0.0       0.0   0.0      0.0       0.0
 -1.02366  0.0  0.0  1.12502  -0.309105  0.0  -2.59127  1.42393  -1.63988