Infer the return type for fillvalue based on index types

This should promote greater type-stability

Author: timholy

src/b-splines/indexing.jl

@@ -69,6 +69,15 @@ function gradient_impl{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad
     end
 end
 
+function getindex_return_type{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad}(::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,Pad}}, argtypes)
+    Tret = TCoefs
+    for a in argtypes
+        Tret = Base.promote_op(Base.MulFun, T, a)
+    end
+    Tret
+end
+
+
 @generated function gradient!{T,N}(g::AbstractVector, itp::BSplineInterpolation{T,N}, xs::Number...)
     length(xs) == N || error("Can only be called with $N indexes")
     gradient_impl(itp)

src/extrapolation/filled.jl

@@ -5,11 +5,7 @@ type FilledExtrapolation{T,N,ITP<:AbstractInterpolation,IT,GT,FT} <: AbstractExt
     itp::ITP
     fillvalue::FT
 end
-"""
-`FilledExtrapolation(itp, fillvalue)` creates an extrapolation object that returns the `fillvalue` any time the indexes in `itp[x1,x2,...]` are out-of-bounds.
 
-By comparison with `extrapolate`, this version lets you control the `fillvalue`'s type directly.  It's important for the `fillvalue` to be of the same type as returned by `itp[x1,x2,...]` for in-bounds regions for the index types you are using; otherwise, indexing will be type-unstable (and slow).
-"""
 function FilledExtrapolation{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue)
     FilledExtrapolation{T,N,typeof(itp),IT,GT,typeof(fillvalue)}(itp, fillvalue)
 end
@@ -19,15 +15,16 @@ end
 """
 extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue) = FilledExtrapolation(itp, convert(eltype(itp), fillvalue))
 
-@generated function getindex{T,N}(fitp::FilledExtrapolation{T,N}, args::Number...)
+@generated function getindex{T,N,ITP,IT,GT,FT}(fitp::FilledExtrapolation{T,N,ITP,IT,GT,FT}, args::Number...)
     n = length(args)
     n == N || return error("Must index $(N)-dimensional interpolation objects with $(nindexes(N))")
+    Tret = FT<:Number ? getindex_return_type(ITP, args) : FT
     meta = Expr(:meta, :inline)
     quote
         $meta
         # Check to see if we're in the extrapolation region, i.e.,
         # out-of-bounds in an index
-        @nexprs $N d->((args[d] < lbound(fitp,d) || args[d] > ubound(fitp, d)) && return fitp.fillvalue)
+        @nexprs $N d->((args[d] < lbound(fitp,d) || args[d] > ubound(fitp, d)) && return convert($Tret, fitp.fillvalue))
         # In the interpolation region
         return getindex(fitp.itp,args...)
     end

src/gridded/indexing.jl

@@ -86,3 +86,11 @@ function getindex{T,N,TCoefs,IT<:DimSpec{Gridded},K,P}(itp::GriddedInterpolation
     dest = Array(T, map(length, x))::Array{T,N}
     getindex!(dest, itp, x...)
 end
+
+function getindex_return_type{T,N,TCoefs,IT<:DimSpec{Gridded},K,P}(::Type{GriddedInterpolation{T,N,TCoefs,IT,K,P}}, argtypes)
+    Tret = TCoefs
+    for a in argtypes
+        Tret = Base.promote_op(Base.MulFun, T, a)
+    end
+    Tret
+end

test/extrapolation/runtests.jl

@@ -27,6 +27,8 @@ etpf = @inferred(extrapolate(itpg, NaN))
 @test_throws BoundsError etpf[2.5,2]
 @test_throws ErrorException etpf[2.5,2,1]  # this will probably become a BoundsError someday
 
+@inferred(getindex(etpf, dual(-2.5,1)))
+
 etpf = @inferred(FilledExtrapolation(itpg, 'x'))
 @test_approx_eq etpf[2] f(2)
 @test etpf[-1.5] == 'x'
@@ -61,4 +63,4 @@ etp2ll = extrapolate(itp2g, Linear)
 
 end
 
-include("type-stability.jl")
+include("type-stability.jl")