Function call syntax support

README.md

@@ -52,6 +52,11 @@ To evaluate the interpolation at position `(x, y, ...)`, simply do
 v = itp[x, y, ...]
 ```
 
+or, alternatively, using function call syntax:
+```jl
+v = itp(x, y, ...)
+```
+
 Some interpolation objects support computation of the gradient, which
 can be obtained as
 ```jl

src/b-splines/indexing.jl

@@ -51,6 +51,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

test/scaling/function-call-syntax.jl

@@ -0,0 +1,100 @@
+module ScalingFunctionCallTests
+
+using Interpolations
+using Base.Test
+
+# Model linear interpolation of y = -3 + .5x by interpolating y=x
+# and then scaling to the new x range
+
+itp = interpolate(1:1.0:10, BSpline(Linear()), OnGrid())
+
+sitp = @inferred(scale(itp, -3:.5:1.5))
+@test typeof(sitp) <: Interpolations.ScaledInterpolation
+@test parent(sitp) === itp
+
+for (x,y) in zip(-3:.05:1.5, 1:.1:10)
+    @test sitp(x) ≈ y
+end
+
+# Verify that it works in >1D, with different types of ranges
+
+gauss(phi, mu, sigma) = exp(-(phi-mu)^2 / (2sigma)^2)
+testfunction(x,y) = gauss(x, 0.5, 4) * gauss(y, -.5, 2)
+
+xs = -5:.5:5
+ys = -4:.2:4
+zs = Float64[testfunction(x,y) for x in xs, y in ys]
+
+itp2 = interpolate(zs, BSpline(Quadratic(Flat())), OnGrid())
+sitp2 = @inferred scale(itp2, xs, ys)
+
+for x in xs, y in ys
+    @test testfunction(x,y) ≈ sitp2(x,y)
+end
+
+# Iteration
+itp = interpolate(rand(3,3,3), BSpline(Quadratic(Flat())), OnCell())
+knots = map(d->1:10:21, 1:3)
+sitp = @inferred scale(itp, knots...)
+
+iter = @inferred(eachvalue(sitp))
+state = @inferred(start(iter))
+@test !(@inferred(done(iter, state)))
+val, state = @inferred(next(iter, state))
+
+function foo!(dest, sitp)
+    i = 0
+    for s in eachvalue(sitp)
+        dest[i+=1] = s
+    end
+    dest
+end
+function bar!(dest, sitp)
+    for I in CartesianRange(size(dest))
+        dest[I] = sitp(I)
+    end
+    dest
+end
+rfoo = Array{Float64}( Interpolations.ssize(sitp))
+rbar = similar(rfoo)
+foo!(rfoo, sitp)
+bar!(rbar, sitp)
+@test rfoo ≈ rbar
+
+# with extrapolation
+END = 10
+xs = linspace(-5, 5, END)
+ys = map(sin, xs)
+
+function run_tests(sut::Interpolations.AbstractInterpolation{T,N,IT,OnGrid}, itp) where {T,N,IT}
+    for x in xs
+        @test ≈(sut[x],sin(x),atol=sqrt(eps(sin(x))))
+    end
+    @test sut(-5) == sut(-5.1) == sut(-15.8) == sut(-Inf) == itp(1)
+    @test sut(5) == sut(5.1) == sut(15.8) == sut(Inf) == itp(END)
+end
+
+function run_tests(sut::Interpolations.AbstractInterpolation{T,N,IT,OnCell}, itp) where {T,N,IT}
+    halfcell = (xs[2] - xs[1]) / 2
+
+    for x in (5 + halfcell, 5 + 1.1halfcell, 15.8, Inf)
+        @test sut(-x) == itp(.5)
+        @test sut(x) == itp(END+.5)
+    end
+end
+
+for GT in (OnGrid, OnCell)
+    itp = interpolate(ys, BSpline(Quadratic(Flat())), GT())
+
+    # Test extrapolating, then scaling
+    eitp = extrapolate(itp, Flat())
+    seitp = scale(eitp, xs)
+    run_tests(seitp, itp)
+
+    # Test scaling, then extrapolating
+    sitp = scale(itp, xs)
+    esitp = extrapolate(sitp, Flat())
+    run_tests(esitp, itp)
+end
+
+end