Filled extrapolation: do not assume BSpline. Improve README for extrap.

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

@@ -26,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/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