Merge pull request #233 from JuliaMath/teh/fix_232

Clean-ups to extrapolate

Author: timholy

README.md

@@ -152,6 +152,13 @@ extrap = LinearInterpolation(xs, A, extrapolation_bc = Line())
 @test extrap(1 - 0.2) # ≈ f(1) - (f(1.2) - f(1))
 @test extrap(5 + 0.2) # ≈ f(5) + (f(5) - f(4.8))
 ```
+You can also use a "fill" value, which gets returned whenever you ask for out-of-range values:
+
+```julia
+extrap = LinearInterpolation(xs, A, extrapolation_bc = NaN)
+@test isnan(extrap(5.2))
+```
+
 Irregular grids are supported as well; note that presently only `LinearInterpolation` supports irregular grids.
 ```julia
 xs = [x^2 for x = 1:0.2:5]
@@ -303,7 +310,18 @@ plot!(xs, ys, label="spline")
 
 ## Extrapolation
 
-The call to `extrapolate` defines what happens if you try to index into the interpolation object with coordinates outside of `[1, size(data,d)]` in any dimension `d`. The implemented boundary conditions are `Throw`, `Flat`, `Linear`, `Periodic` and `Reflect`, with more options planned. `Periodic` and `Reflect` require that there is a method of `Base.mod` that can handle the indices used.
+The call to `extrapolate` defines what happens if you try to index into the interpolation object with coordinates outside of its
+bounds in any dimension. The implemented boundary conditions are `Throw`, `Flat`, `Linear`, `Periodic` and `Reflect`,
+or you can pass a constant to be used as a "fill" value returned for any out-of-bounds evaluation.
+`Periodic` and `Reflect` require that there is a method of `Base.mod` that can handle the indices used.
+
+Examples:
+
+```
+itp = interpolate(1:7, BSpline(Linear()))
+etpf = extrapolate(itp, Flat())   # gives 1 on the left edge and 7 on the right edge
+etp0 = extrapolate(itp, 0)        # gives 0 everywhere outside [1,7]
+```
 
 ## Performance shootout
 

src/extrapolation/filled.jl

@@ -8,6 +8,11 @@ function FilledExtrapolation(itp::AbstractInterpolation{T,N,IT}, fillvalue) wher
     FilledExtrapolation{Te,N,typeof(itp),IT,typeof(fillvalue)}(itp, fillvalue)
 end
 
+@noinline FilledExtrapolation(itp::AbstractInterpolation{T,N,IT}, ::Type{F}) where {T,N,IT,F} =
+    throw(ArgumentError("cannot create a filled extrapolation with a type $F, use a value of this type instead (e.g., $F(0))"))
+@noinline FilledExtrapolation(itp::AbstractInterpolation{T,N,IT}, ::Type{F}) where {T,N,IT,F<:BoundaryCondition} =
+    throw(ArgumentError("cannot create a filled extrapolation with a type $F, use a value of this type instead (e.g., $F())"))
+
 Base.parent(A::FilledExtrapolation) = A.itp
 etpflag(A::FilledExtrapolation) = A.fillvalue
 itpflag(A::FilledExtrapolation) = itpflag(A.itp)
@@ -21,8 +26,7 @@ extrapolate(itp::AbstractInterpolation{T,N,IT}, fillvalue) where {T,N,IT} = Fill
     itp = parent(etp)
     Tret = typeof(prod(x) * zero(T))
     if checkbounds(Bool, itp, x...)
-        wis = weightedindexes((value_weights,), itpinfo(itp)..., x)
-        convert(Tret, itp.coefs[wis...])
+        @inbounds itp(x...)
     else
         convert(Tret, etp.fillvalue)
     end

test/extrapolation/runtests.jl

@@ -91,6 +91,10 @@ using Test
     @test itpe(10.1, 1) ≈ 10.1
     @test_throws BoundsError itpe(9.9, 0)
 
+    # Issue #232
+    targeterr = ArgumentError("cannot create a filled extrapolation with a type Line, use a value of this type instead (e.g., Line())")
+    @test_throws targeterr extrapolate(itp, Line)
+
     include("type-stability.jl")
     include("non1.jl")
 end

test/readme-examples.jl

@@ -51,4 +51,26 @@ using Interpolations, Test
     itp = interpolate(knots, A, Gridded(Linear()))
     @test itp(4,1.2) ≈ A[2,6] atol=.1  # approximately A[2,6]
 
+    ## extrapolate
+    itp = interpolate(1:7, BSpline(Linear()))
+    etpf = extrapolate(itp, Flat())
+    @test etpf(0) == 1
+    @test etpf(1) == 1
+    @test etpf(7) == 7
+    @test etpf(7.8) == 7
+    etp0 = extrapolate(itp, 0)
+    @test etp0(0) === 0.0
+    @test etp0(1) === 1.0
+    @test etp0(7) === 7.0
+    @test etp0(7.8) === 0.0
+
+    f(x) = log(x)
+    xs = 1:0.2:5
+    A = [f(x) for x in xs]
+    extrap = LinearInterpolation(xs, A, extrapolation_bc = Line())
+    @test extrap(1 - 0.2) ≈ f(1) - (f(1.2) - f(1))
+    @test extrap(5 + 0.2) ≈ f(5) + (f(5) - f(4.8))
+    extrap = LinearInterpolation(xs, A, extrapolation_bc = NaN)
+    @test isnan(extrap(5.2))
+
 end