ForwardDiff not working for simple 2D linear interpolation

[arbenede]

I am trying to compute the Jacobian of a plain Sum of Square Difference (SSD) cost function between two shifted images where interpolation is computed using Interpolations.jl with respect to the shift variables. Here is the code:

using Images
using TestImages
using ForwardDiff
using Interpolations

img1 = Gray.(testimage("cameraman.tif"))
mat1 = convert(Array{Float64}, img1)
(nr, nc) = size(mat1)

# Define image shift ranges
Δx_min = -2.5
Δx_max = 2.5

Δy_min = -2.5
Δy_max = 2.5

# Define img1 coordinates
xcoord = range(ceil(Int, abs(Δx_min)), stop=nc - 1 - ceil(Int, abs(Δx_max)))
ycoord = range(ceil(Int, abs(Δy_min)), stop=nr - 1 - ceil(Int, abs(Δy_max)))

Δx = -1.34
Δy = 0.78

# Create interpolator
interp_linear = LinearInterpolation((range(0,stop=nr-1), range(0,stop=nc-1)), mat1)
mat2 = interp_linear(ycoord .+ Δy, xcoord .+ Δx)

function cost(Δ::Vector{T}) where T <: Real
    subimg1 = mat1[ycoord, xcoord]
    subimg2 = interp_linear(ycoord .+ Δ[2], xcoord .+ Δ[1])
    
    @. subimg1 = (subimg1 - subimg2)^2
    SSD = sum(vec(subimg1))
    
    return SSD
end

ForwardDiff.gradient(cost, [0.0,0.0])

and I get this error:

MethodError: no method matching Integer(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2})
Closest candidates are:
  (::Type{T})(::T) where T<:Number at boot.jl:760
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  (::Type{T})(::Base.TwicePrecision) where T<:Number at twiceprecision.jl:243
  ...

Stacktrace:
  [1] unsafe_length(r::UnitRange{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}})
    @ Base ./range.jl:559
  [2] length(r::UnitRange{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}})
    @ Base ./range.jl:561
  [3] size(r::UnitRange{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}})
    @ Base ./range.jl:515
  [4] axes
    @ ./abstractarray.jl:89 [inlined]
  [5] axes1
    @ ./abstractarray.jl:109 [inlined]
  [6] shape(v::UnitRange{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}}, rest::UnitRange{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}})
    @ Interpolations ~/.julia/packages/Interpolations/GIn2o/src/utils.jl:71
  [7] (::Interpolations.Extrapolation{Float64, 2, ScaledInterpolation{Float64, 2, Interpolations.BSplineInterpolation{Float64, 2, Matrix{Float64}, BSpline{Linear}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, BSpline{Linear}, Tuple{UnitRange{Int64}, UnitRange{Int64}}}, BSpline{Linear}, Throw{Nothing}})(::UnitRange{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}}, ::UnitRange{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}})
    @ Interpolations ~/.julia/packages/Interpolations/GIn2o/src/extrapolation/extrapolation.jl:56
  [8] cost(Δ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}})
    @ Main ./In[10]:3
  [9] vector_mode_dual_eval
    @ ~/.julia/packages/ForwardDiff/QOqCN/src/apiutils.jl:37 [inlined]
 [10] vector_mode_gradient(f::typeof(cost), x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:106
 [11] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:19
 [12] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}}}) (repeats 2 times)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:17
 [13] top-level scope
    @ In[12]:1
 [14] eval
    @ ./boot.jl:360 [inlined]
 [15] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1094

By the way, I also tried using the translation image transformation defined in ImageTransformations.jl with no success. More details on that attempt are at forwarddiff-and-imagetransformations
It seems that ForwardDiff should be able to handle this, right?

[MatFi]

Hmm converting the coord's to Float64 does the trick. But not sure where the problem comes from.

xcoordf64 = Float64.(xcoord)
ycoordf64 = Float64.(ycoord)
function cost(Δ::Vector{T}) where T <: Real
    subimg1 = mat1[ycoord, xcoord]
    subimg2 = interp_linear(ycoordf64 .+ Δ[2], xcoordf64 .+ Δ[1]) 
    subimg = (subimg1 - subimg2)^2
    SSD = sum(vec(subimg))
    
    return SSD
end

[arbenede]

Hmm... I am still getting an error:

ForwardDiff.gradient(cost, [1.0,1.0])

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

Stacktrace:
  [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2})
    @ Base ./number.jl:7
  [2] setindex!(A::Matrix{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}, i1::Int64)
    @ Base ./array.jl:839
  [3] (::Interpolations.Extrapolation{Float64, 2, ScaledInterpolation{Float64, 2, Interpolations.BSplineInterpolation{Float64, 2, Matrix{Float64}, BSpline{Linear}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, BSpline{Linear}, Tuple{UnitRange{Int64}, UnitRange{Int64}}}, BSpline{Linear}, Throw{Nothing}})(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}}, ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}})
    @ Interpolations ~/.julia/packages/Interpolations/GIn2o/src/extrapolation/extrapolation.jl:58
  [4] cost(Δ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}})
    @ Main ./In[4]:3
  [5] vector_mode_dual_eval
    @ ~/.julia/packages/ForwardDiff/QOqCN/src/apiutils.jl:37 [inlined]
  [6] vector_mode_gradient(f::typeof(cost), x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:106
  [7] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:19
  [8] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}}}) (repeats 2 times)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:17
  [9] top-level scope
    @ In[6]:1
 [10] eval
    @ ./boot.jl:360 [inlined]
 [11] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1094

[MatFi]

Be aware that you can not copy ForwardDiff.Dual's to subimg1 as it has Float64 as eltype. Either you reinterpret or allow allocations. For this reason I did subimg = (subimg1 - subimg2)^2
edit:
subimg = @. (subimg1 - subimg2)^2 is how it actually should be..

[arbenede]

I am still getting an error, so I put the notebook in this gist. Please feel free to fix it!

[mkitti]

I am unclear what I can do for this from perspective of Interpolations.jl. None of the errors posted so far point to Interpolations.jl. Perhaps I am missing something.

[arbenede]

I was under the impression that @MatFi was able to modify the code that I posted and run it with no errors...

[johnnychen94]

Out of curiosity, will Zygote work for this? The master branch of Interpolations has ChainRules support so it might just work out of the box.

[MatFi]

I am still getting an error, so I put the notebook in this gist. Please feel free to fix it!

Copy pasted your gist

julia> ForwardDiff.gradient(cost, [1.0,1.0])
2-element Vector{Float64}:
 449.4649134948096
 227.50594386774316

julia> versioninfo()
Julia Version 1.6.0
Commit f9720dc2eb (2021-03-24 12:55 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: AMD Ryzen 9 4900HS with Radeon Graphics
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, znver2)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 8

[mkitti]

As far as I can tell, no one has suggested a modification to Interpolations.jl yet or otherwise pointed to a defect in Interpolations.jl. Please let me know if I am wrong. I think the right venue for this conversation may be https://discourse.julialang.org/

[arbenede]

I run some tests under Linux with Julia 1.6.1 since @MatFi reported that he can the code there. I am still getting an error that points to some issue with Interpolations.jl:

(venv-jupyter) vagrant@vagrant:/vagrant/julia$ julia test_AD.jl 
ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:760
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  ...
Stacktrace:
 [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2})
   @ Base ./number.jl:7
 [2] setindex!(A::Matrix{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}, i1::Int64)
   @ Base ./array.jl:839
 [3] (::Interpolations.Extrapolation{Float64, 2, ScaledInterpolation{Float64, 2, Interpolations.BSplineInterpolation{Float64, 2, Matrix{Float64}, BSpline{Linear}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, BSpline{Linear}, Tuple{UnitRange{Int64}, UnitRange{Int64}}}, BSpline{Linear}, Throw{Nothing}})(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}}, ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}})
   @ Interpolations ~/.julia/packages/Interpolations/GIn2o/src/extrapolation/extrapolation.jl:58
 [4] cost(Δ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(cost), Float64}, Float64, 2}})

By the way @MatFi if you can please report what version of Interpolations.jl are you running I will run a test with Julia 1.6.0 to match your setup 100%. @mkitti thanks for your suggestion, in fact I posted the same question to https://discourse.julialang.org/ however I did not get much feedback. Feel free to close the issue, I will roll out a simple interpolation function that plays nicely with Dual numbers.

[mkitti]

If I remove the extrapolation using parent, it seems to work fine:

julia> interp_linear = parent(LinearInterpolation((range(0,stop=nr-1), range(0,stop=nc-1)), mat1));

julia> typeof(interp_linear)
ScaledInterpolation{Float64, 2, Interpolations.BSplineInterpolation{Float64, 2, Matrix{Float64}, BSpline{Linear}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, BSpline{Linear}, Tuple{UnitRange{Int64}, UnitRange{Int64}}}

julia> function cost(Δ::Vector{T}) where T <: Real
           subimg1 = mat1[ycoord, xcoord]
           subimg2 = interp_linear(ycoordf64 .+ Δ[2], xcoordf64 .+ Δ[1])

           subimg = @. (subimg1 - subimg2)^2
           SSD = sum(vec(subimg))

           return SSD
       end
cost (generic function with 1 method)

julia> ForwardDiff.gradient(cost, [1.0, 1.0])
2-element Vector{Float64}:
 449.4649134948097
 227.50594386774316

@arbenede, does that work for you?

[arbenede]

Yes, that works fine, thanks! I think it would be a good idea to mention this in the documentation or to change the behavior of the interpolator so other people will not have to go through this...

[mkitti]

The above is a workaround. We still need to determine the root cause.

The issue is that when a multidimensional extrapolation is given several vectors of coordinates, it does not take into account the eltype of the input when building the return type. We should be able to modify that by including a call to promote_type.

julia> import Interpolations: Extrapolation, shape

julia> @inline function (etp::Extrapolation{T,N})(x::Vararg{Union{Number,AbstractVector},N}) where {T,N}
           itp = parent(etp)
           ret = zeros(promote_type(T, ( eltype.(x) )... ), shape(x...))
           for (i, y) in zip(eachindex(ret), Iterators.product(x...))
               ret[i] = etp(y...)
           end
           return ret
       end

julia> function cost(Δ::Vector{T}) where T <: Real
           subimg1 = mat1[ycoord, xcoord]
           subimg2 = interp_linear(ycoordf64 .+ Δ[2], xcoordf64 .+ Δ[1])

           subimg = @. (subimg1 - subimg2)^2
           SSD = sum(vec(subimg))

           return SSD
       end
cost (generic function with 1 method)

julia> ForwardDiff.gradient(cost, [1.0, 1.0])
2-element Vector{Float64}:
 449.4649134948097
 227.50594386774316

Compare the above to

Interpolations.jl/src/extrapolation/extrapolation.jl

Lines 54 to 61 in 7d32f82

	@inline function (etp::Extrapolation{T,N})(x::Vararg{Union{Number,AbstractVector},N}) where {T,N}
	itp = parent(etp)
	ret = zeros(T, shape(x...))
	for (i, y) in zip(eachindex(ret), Iterators.product(x...))
	ret[i] = etp(y...)
	end
	return ret
	end

This mirrors the change in 77270b4.