Function call syntax support (#206)

* Function call syntax support

* README.md in favour of function calls

Author: chiyahn

README.md

@@ -39,6 +39,8 @@ from the Julia REPL.
 
 ## General usage
 
+Note: the current version of `Interpolations` supports interpolation evaluation using index calls `[]`, but this feature will be deprecated in future. We highly recommend function calls with `()` as follows.
+
 Given an `AbstractArray` `A`, construct an "interpolation object" `itp` as
 ```jl
 itp = interpolate(A, options...)
@@ -49,7 +51,7 @@ samples in `A` are equally-spaced.
 
 To evaluate the interpolation at position `(x, y, ...)`, simply do
 ```jl
-v = itp[x, y, ...]
+v = itp(x, y, ...)
 ```
 
 Some interpolation objects support computation of the gradient, which
@@ -89,14 +91,14 @@ Julia's iterator objects, e.g.,
 ```jl
 function ongrid!(dest, itp)
     for I in CartesianRange(size(itp))
-        dest[I] = itp[I]
+        dest[I] = itp(I)
     end
 end
 ```
 would store the on-grid value at each grid point of `itp` in the output `dest`.
 Finally, courtesy of Julia's indexing rules, you can also use
 ```jl
-fine = itp[linspace(1,10,1001), linspace(1,15,201)]
+fine = itp(linspace(1,10,1001), linspace(1,15,201))
 ```
 
 
@@ -114,19 +116,19 @@ Some examples:
 ```jl
 # Nearest-neighbor interpolation
 itp = interpolate(a, BSpline(Constant()), OnCell())
-v = itp[5.4]   # returns a[5]
+v = itp(5.4)   # returns a[5]
 
 # (Multi)linear interpolation
 itp = interpolate(A, BSpline(Linear()), OnGrid())
-v = itp[3.2, 4.1]  # returns 0.9*(0.8*A[3,4]+0.2*A[4,4]) + 0.1*(0.8*A[3,5]+0.2*A[4,5])
+v = itp(3.2, 4.1)  # returns 0.9*(0.8*A[3,4]+0.2*A[4,4]) + 0.1*(0.8*A[3,5]+0.2*A[4,5])
 
 # Quadratic interpolation with reflecting boundary conditions
 # Quadratic is the lowest order that has continuous gradient
 itp = interpolate(A, BSpline(Quadratic(Reflect())), OnCell())
 
 # Linear interpolation in the first dimension, and no interpolation (just lookup) in the second
 itp = interpolate(A, (BSpline(Linear()), NoInterp()), OnGrid())
-v = itp[3.65, 5]  # returns  0.35*A[3,5] + 0.65*A[4,5]
+v = itp(3.65, 5)  # returns  0.35*A[3,5] + 0.65*A[4,5]
 ```
 There are more options available, for example:
 ```jl
@@ -145,8 +147,8 @@ A_x = 1.:2.:40.
 A = [log(x) for x in A_x]
 itp = interpolate(A, BSpline(Cubic(Line())), OnGrid())
 sitp = scale(itp, A_x)
-sitp[3.] # exactly log(3.)
-sitp[3.5] # approximately log(3.5)
+sitp(3.) # exactly log(3.)
+sitp(3.5) # approximately log(3.5)
 ```
 
 For multidimensional uniformly spaced grids
@@ -157,8 +159,8 @@ f(x1, x2) = log(x1+x2)
 A = [f(x1,x2) for x1 in A_x1, x2 in A_x2]
 itp = interpolate(A, BSpline(Cubic(Line())), OnGrid())
 sitp = scale(itp, A_x1, A_x2)
-sitp[5., 10.] # exactly log(5 + 10)
-sitp[5.6, 7.1] # approximately log(5.6 + 7.1)
+sitp(5., 10.) # exactly log(5 + 10)
+sitp(5.6, 7.1) # approximately log(5.6 + 7.1)
 ```
 ### Gridded interpolation
 
@@ -173,7 +175,7 @@ A = rand(20)
 A_x = collect(1.0:2.0:40.0)
 knots = (A_x,)
 itp = interpolate(knots, A, Gridded(Linear()))
-itp[2.0]
+itp(2.0)
 ```
 
 The spacing between adjacent samples need not be constant, you can use the syntax
@@ -188,7 +190,7 @@ For example:
 A = rand(8,20)
 knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
 itp = interpolate(knots, A, Gridded(Linear()))
-itp[4,1.2]  # approximately A[2,6]
+itp(4,1.2)  # approximately A[2,6]
 ```
 One may also mix modes, by specifying a mode vector in the form of an explicit tuple:
 ```jl
@@ -225,7 +227,7 @@ A = hcat(x,y)
 itp = scale(interpolate(A, (BSpline(Cubic(Natural())), NoInterp()), OnGrid()), t, 1:2)
 
 tfine = 0:.01:1
-xs, ys = [itp[t,1] for t in tfine], [itp[t,2] for t in tfine]
+xs, ys = [itp(t,1) for t in tfine], [itp(t,2) for t in tfine]
 ```
 
 We can then plot the spline with:

src/b-splines/indexing.jl

@@ -50,6 +50,11 @@ end
     getindex_impl(itp)
 end
 
+function (itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad})(args...) where {T,N,TCoefs,IT,GT,pad}
+    # support function calls
+    itp[args...]
+end
+
 function gradient_impl(itp::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,Pad}}) where {T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad}
     meta = Expr(:meta, :inline)
     # For each component of the gradient, alternately calculate

src/extrapolation/extrapolation.jl

@@ -63,6 +63,11 @@ end
     getindex_impl(etp, xs...)
 end
 
+function (etp::Extrapolation{T,N,ITPT,IT,GT,ET})(args...) where {T,N,ITPT,IT,GT,ET}
+    # support function calls
+    etp[args...]
+end
+
 checkbounds(::AbstractExtrapolation,I...) = nothing
 
 

src/extrapolation/filled.jl

@@ -25,6 +25,11 @@ extrapolate(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue) where {T,N,IT,GT}
     convert(Tret, fitp.fillvalue)
 end
 
+function (fitp::FilledExtrapolation{T,N,ITP,IT,GT,FT})(args...) where {T,N,ITP,IT,GT,FT}
+    # support function calls
+    fitp[args...]
+end
+
 @inline Base.checkbounds(::Type{Bool}, A::FilledExtrapolation, I...) = _checkbounds(A, 1, indices(A), I)
 @inline _checkbounds(A, d::Int, IA::TT1, I::TT2) where {TT1,TT2} =
     (I[1] >= lbound(A, d, IA[1])) & (I[1] <= ubound(A, d, IA[1])) & _checkbounds(A, d+1, Base.tail(IA), Base.tail(I))

src/gridded/indexing.jl

@@ -36,6 +36,11 @@ end
     :(getindex(itp, $(args...)))
 end
 
+function (itp::GriddedInterpolation{T,N,TCoefs,IT,K,pad})(args...) where {T,N,TCoefs,IT,K,pad}
+    # support function calls
+    itp[args...]
+end
+
 # Indexing with vector inputs. Here, it pays to pre-process the input indexes,
 # because N*n is much smaller than n^N.
 # TODO: special-case N=1, because there is no reason to separately cache the indexes.

src/scaling/scaling.jl

@@ -45,6 +45,10 @@ end
 
 getindex(sitp::ScaledInterpolation{T,1}, x::Number, y::Int) where {T} = y == 1 ? sitp[x] : throw(BoundsError())
 
+function (sitp::ScaledInterpolation{T,N,ITPT,IT})(args...) where {T,N,ITPT,IT<:DimSpec}
+    sitp[args...]
+end
+
 size(sitp::ScaledInterpolation, d) = size(sitp.itp, d)
 lbound(sitp::ScaledInterpolation{T,N,ITPT,IT,OnGrid}, d) where {T,N,ITPT,IT} = 1 <= d <= N ? sitp.ranges[d][1] : throw(BoundsError())
 lbound(sitp::ScaledInterpolation{T,N,ITPT,IT,OnCell}, d) where {T,N,ITPT,IT} = 1 <= d <= N ? sitp.ranges[d][1] - boundstep(sitp.ranges[d]) : throw(BoundsError())

test/b-splines/function-call-syntax.jl

@@ -0,0 +1,24 @@
+module ExtrapFunctionCallSyntax
+
+using Base.Test, Interpolations, DualNumbers
+
+# Test if b-spline interpolation by function syntax yields identical results
+f(x) = sin((x-3)*2pi/9 - 1)
+xmax = 10
+A = Float64[f(x) for x in 1:xmax]
+itpg = interpolate(A, BSpline(Linear()), OnGrid())
+schemes = (Flat,Line,Free)
+
+for T in (Cubic, Quadratic), GC in (OnGrid, OnCell)
+    for etp in map(S -> @inferred(interpolate(A, BSpline(T(S())), GC())), schemes),
+        x in linspace(1,xmax,100)
+        @test (getindex(etp, x)) == etp(x)
+    end
+end
+
+for T in (Constant, Linear), GC in (OnGrid, OnCell), x in linspace(1,xmax,100)
+    etp = interpolate(A, BSpline(T()), GC())
+    @test (getindex(etp, x)) == etp(x)
+end
+
+end

test/b-splines/runtests.jl

@@ -7,5 +7,6 @@ include("cubic.jl")
 include("mixed.jl")
 include("multivalued.jl")
 include("non1.jl")
+include("function-call-syntax.jl")
 
 end

test/extrapolation/function-call-syntax.jl

@@ -0,0 +1,104 @@
+module ExtrapFunctionCallSyntax
+
+using Base.Test, Interpolations, DualNumbers
+
+# Test if extrapolation by function syntax yields identical results
+f(x) = sin((x-3)*2pi/9 - 1)
+xmax = 10
+A = Float64[f(x) for x in 1:xmax]
+itpg = interpolate(A, BSpline(Linear()), OnGrid())
+
+schemes = (
+    Flat,
+    Linear,
+    Reflect,
+    Periodic
+)
+
+for etp in map(E -> @inferred(extrapolate(itpg, E())), schemes),
+    x in [
+        # In-bounds evaluation
+        3.4, 3, dual(3.1),
+        # Out-of-bounds evaluation
+        -3.4, -3, dual(-3,1),
+        13.4, 13, dual(13,1)
+    ]
+    @test (getindex(etp, x)) == etp(x)
+end
+
+etpg = extrapolate(itpg, Flat())
+@test typeof(etpg) <: AbstractExtrapolation
+
+@test etpg(-3) == etpg(-4.5) == etpg(0.9) == etpg(1.0) == A[1]
+@test etpg(10.1) == etpg(11) == etpg(148.298452) == A[end]
+
+etpf = @inferred(extrapolate(itpg, NaN))
+@test typeof(etpf) <: Interpolations.FilledExtrapolation
+@test parent(etpf) === itpg
+
+@test @inferred(size(etpf)) == (xmax,)
+@test isnan(@inferred(etpf(-2.5)))
+@test isnan(etpf(0.999))
+@test @inferred(etpf(1)) ≈ f(1)
+@test etpf(10) ≈ f(10)
+@test isnan(@inferred(etpf(10.001)))
+
+@test etpf(2.5,1) == etpf(2.5)   # for show method
+@test_throws BoundsError etpf(2.5,2)
+@test_throws BoundsError etpf(2.5,2,1)
+
+x =  @inferred(etpf(dual(-2.5,1)))
+@test isa(x, Dual)
+
+etpl = extrapolate(itpg, Linear())
+k_lo = A[2] - A[1]
+x_lo = -3.2
+@test etpl(x_lo) ≈ A[1] + k_lo * (x_lo - 1)
+k_hi = A[end] - A[end-1]
+x_hi = xmax + 5.7
+@test etpl(x_hi) ≈ A[end] + k_hi * (x_hi - xmax)
+
+xmax, ymax = 8,8
+g(x, y) = (x^2 + 3x - 8) * (-2y^2 + y + 1)
+
+itp2g = interpolate(Float64[g(x,y) for x in 1:xmax, y in 1:ymax], (BSpline(Quadratic(Free())), BSpline(Linear())), OnGrid())
+etp2g = extrapolate(itp2g, (Linear(), Flat()))
+
+@test @inferred(etp2g(-0.5,4)) ≈ itp2g(1,4) - 1.5 * epsilon(etp2g(dual(1,1),4))
+@test @inferred(etp2g(5,100)) ≈ itp2g(5,ymax)
+
+etp2ud = extrapolate(itp2g, ((Linear(), Flat()), Flat()))
+@test @inferred(etp2ud(-0.5,4)) ≈ itp2g(1,4) - 1.5 * epsilon(etp2g(dual(1,1),4))
+@test @inferred(etp2ud(5, -4)) == etp2ud(5,1)
+@test @inferred(etp2ud(100, 4)) == etp2ud(8,4)
+@test @inferred(etp2ud(-.5, 100)) == itp2g(1,8) - 1.5 * epsilon(etp2g(dual(1,1),8))
+
+etp2ll = extrapolate(itp2g, Linear())
+@test @inferred(etp2ll(-0.5,100)) ≈ (itp2g(1,8) - 1.5 * epsilon(etp2ll(dual(1,1),8))) + (100 - 8) * epsilon(etp2ll(1,dual(8,1)))
+
+# Allow element types that don't support conversion to Int (#87):
+etp87g = extrapolate(interpolate([1.0im, 2.0im, 3.0im], BSpline(Linear()), OnGrid()), 0.0im)
+@test @inferred(etp87g(1)) == 1.0im
+@test @inferred(etp87g(1.5)) == 1.5im
+@test @inferred(etp87g(0.75)) == 0.0im
+@test @inferred(etp87g(3.25)) == 0.0im
+
+etp87c = extrapolate(interpolate([1.0im, 2.0im, 3.0im], BSpline(Linear()), OnCell()), 0.0im)
+@test @inferred(etp87c(1)) == 1.0im
+@test @inferred(etp87c(1.5)) == 1.5im
+@test @inferred(etp87c(0.75)) == 0.75im
+@test @inferred(etp87c(3.25)) == 3.25im
+@test @inferred(etp87g(0)) == 0.0im
+@test @inferred(etp87g(3.7)) == 0.0im
+
+# Make sure it works with Gridded too
+etp100g = extrapolate(interpolate(([10;20],),[100;110], Gridded(Linear())), Flat())
+@test @inferred(etp100g(5)) == 100
+@test @inferred(etp100g(15)) == 105
+@test @inferred(etp100g(25)) == 110
+# issue #178
+a = randn(10,10) + im*rand(10,10)
+etp = @inferred(extrapolate(interpolate((1:10, 1:10), a, Gridded(Linear())), 0.0))
+@test @inferred(etp(-1,0)) === 0.0+0.0im
+
+end

test/extrapolation/runtests.jl

@@ -103,3 +103,4 @@ end
 
 include("type-stability.jl")
 include("non1.jl")
+include("function-call-syntax.jl")

test/gridded/function-call-syntax.jl

@@ -0,0 +1,74 @@
+module GriddedFunctionCallSyntax
+
+using Interpolations, Base.Test
+
+for D in (Constant, Linear), G in (OnCell, OnGrid)
+    ## 1D
+    a = rand(5)
+    knots = (collect(linspace(1,length(a),length(a))),)
+    itp = @inferred(interpolate(knots, a, Gridded(D())))
+    @inferred(getindex(itp, 2))
+    @inferred(getindex(itp, CartesianIndex((2,))))
+    for i = 2:length(a)-1
+        @test itp(i) ≈ a[i]
+        @test itp(CartesianIndex((i,))) ≈ a[i]
+    end
+    @inferred(getindex(itp, knots...))
+    @test itp[knots...] ≈ a
+    # compare scalar indexing and vector indexing
+    x = knots[1]+0.1
+    v = itp(x)
+    for i = 1:length(x)
+        @test v[i] ≈ itp(x[i])
+    end
+    # check the fallback vector indexing
+    x = [2.3,2.2]   # non-increasing order
+    v = itp[x]
+    for i = 1:length(x)
+        @test v[i] ≈ itp(x[i])
+    end
+    # compare against BSpline
+    itpb = @inferred(interpolate(a, BSpline(D()), G()))
+    for x in linspace(1.1,4.9,101)
+        @test itp(x) ≈ itpb(x)
+    end
+
+    ## 2D
+    A = rand(6,5)
+    knots = (collect(linspace(1,size(A,1),size(A,1))),collect(linspace(1,size(A,2),size(A,2))))
+    itp = @inferred(interpolate(knots, A, Gridded(D())))
+    @test parent(itp) === A
+    @inferred(getindex(itp, 2, 2))
+    @inferred(getindex(itp, CartesianIndex((2,2))))
+    for j = 2:size(A,2)-1, i = 2:size(A,1)-1
+        @test itp(i,j) ≈ A[i,j]
+        @test itp(CartesianIndex((i,j))) ≈ A[i,j]
+    end
+    @test itp[knots...] ≈ A
+    @inferred(getindex(itp, knots...))
+    # compare scalar indexing and vector indexing
+    x, y = knots[1]+0.1, knots[2]+0.6
+    v = itp(x,y)
+    for j = 1:length(y), i = 1:length(x)
+        @test v[i,j] ≈ itp(x[i],y[j])
+    end
+    # check the fallback vector indexing
+    x = [2.3,2.2]   # non-increasing order
+    y = [3.5,2.8]
+    v = itp[x,y]
+    for j = 1:length(y), i = 1:length(x)
+        @test v[i,j] ≈ itp(x[i],y[j])
+    end
+    # compare against BSpline
+    itpb = @inferred(interpolate(A, BSpline(D()), G()))
+    for y in linspace(1.1,5.9,101), x in linspace(1.1,4.9,101)
+        @test itp(x,y) ≈ itpb(x,y)
+    end
+
+    A = rand(8,20)
+    knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
+    itp = interpolate(knots, A, Gridded(D()))
+    @test itp(4,1.2) ≈ A[2,6]
+end
+
+end

test/gridded/runtests.jl

@@ -2,5 +2,6 @@ module GriddedTests
 
 include("gridded.jl")
 include("mixed.jl")
+include("function-call-syntax.jl")
 
-end
+end