Convenience constructors

README.md

@@ -101,6 +101,63 @@ Finally, courtesy of Julia's indexing rules, you can also use
 fine = itp(linspace(1,10,1001), linspace(1,15,201))
 ```
 
+### Quickstart guide
+For linear and cubic spline interpolations, `LinearInterpolation` and `CubicSplineInterpolation` can be used to create interpolation objects handily:
+```jl
+f(x) = log(x)
+xs = 1:0.2:5
+A = [f(x) for x in xs]
+
+# linear interpolation
+interp_linear = LinearInterpolation(xs, A)
+interp_linear[3] # exactly log(3)
+interp_linear[3.1] # approximately log(3.1)
+
+# cubic spline interpolation
+interp_cubic = CubicSplineInterpolation(xs, A)
+interp_cubic[3] # exactly log(3)
+interp_cubic[3.1] # approximately log(3.1)
+```
+which support multidimensional data as well:
+```jl
+f(x,y) = log(x+y)
+xs = 1:0.2:5
+ys = 2:0.1:5
+A = [f(x+y) for x in xs, y in ys]
+
+# linear interpolation
+interp_linear = LinearInterpolation((xs, ys), A)
+interp_linear[3, 2] # exactly log(3 + 2)
+interp_linear[3.1, 2.1] # approximately log(3.1 + 2.1)
+
+# cubic spline interpolation
+interp_cubic = CubicSplineInterpolation((xs, ys), A)
+interp_cubic[3, 2] # exactly log(3 + 2)
+interp_cubic[3.1, 2.1] # approximately log(3.1 + 2.1)
+```
+For extrapolation, i.e., when interpolation objects are evaluated in coordinates outside of range provided in constructors, the default option for a boundary condition is `Throw` so that they will return an error.
+Interested users can specify boundary conditions by providing an extra parameter for `extrapolation_bc`:
+```jl
+f(x) = log(x)
+xs = 1:0.2:5
+A = [f(x) for x in xs]
+
+# extrapolation with linear boundary conditions
+extrap = LinearInterpolation(xs, A, extrapolation_bc = Interpolations.Linear())
+
+@test extrap[1 - 0.2] # ≈ f(1) - (f(1.2) - f(1))
+@test extrap[5 + 0.2] # ≈ f(5) + (f(5) - f(4.8))
+```
+Irregular grids are supported as well; note that presently only `LinearInterpolation` supports irregular grids.
+```jl
+xs = [x^2 for x = 1:0.2:5]
+A = [f(x) for x in xs]
+
+# linear interpolation
+interp_linear = LinearInterpolation(xs, A)
+interp_linear[1] # exactly log(1)
+interp_linear[1.05] # approximately log(1.05)
+```
 
 ## Control of interpolation algorithm
 

src/Interpolations.jl

@@ -27,7 +27,10 @@ export
     Reflect,
     Natural,
     InPlace,
-    InPlaceQ
+    InPlaceQ,
+
+    LinearInterpolation,
+    CubicSplineInterpolation
 
     # see the following files for further exports:
     # b-splines/b-splines.jl
@@ -114,5 +117,6 @@ include("extrapolation/extrapolation.jl")
 include("scaling/scaling.jl")
 include("utils.jl")
 include("io.jl")
+include("convenience-constructors.jl")
 
 end # module

src/convenience-constructors.jl

@@ -0,0 +1,10 @@
+# convenience copnstructors for linear / cubic spline interpolations
+# 1D version
+LinearInterpolation(range::T, vs; extrapolation_bc = Interpolations.Throw()) where {T <: Range} = extrapolate(scale(interpolate(vs, BSpline(Linear()), OnGrid()), range), extrapolation_bc)
+LinearInterpolation(range::T, vs; extrapolation_bc = Interpolations.Throw()) where {T <: AbstractArray} = extrapolate(interpolate((range, ), vs, Gridded(Linear())), extrapolation_bc)
+CubicSplineInterpolation(range::T, vs; bc = Interpolations.Line(), extrapolation_bc = Interpolations.Throw()) where {T <: Range} = extrapolate(scale(interpolate(vs, BSpline(Cubic(bc)), OnGrid()), range), extrapolation_bc)
+
+# multivariate versions
+LinearInterpolation(ranges::NTuple{N,T}, vs; extrapolation_bc = Interpolations.Throw()) where {N,T <: Range} = extrapolate(scale(interpolate(vs, BSpline(Linear()), OnGrid()), ranges...), extrapolation_bc)
+LinearInterpolation(ranges::NTuple{N,T}, vs; extrapolation_bc = Interpolations.Throw()) where {N,T <: AbstractArray} = extrapolate(interpolate(ranges, vs, Gridded(Linear())), extrapolation_bc)
+CubicSplineInterpolation(ranges::NTuple{N,T}, vs; bc = Interpolations.Line(), extrapolation_bc = Interpolations.Throw()) where {N,T <: Range} = extrapolate(scale(interpolate(vs, BSpline(Cubic(bc)), OnGrid()), ranges...), extrapolation_bc)

test/convenience-constructors.jl

@@ -0,0 +1,183 @@
+module ConvenienceConstructorTests
+
+using Interpolations
+using Base.Test
+using Base.Cartesian
+
+# unit test setup
+XMIN = 2
+XMAX = 10
+YMIN = 1
+YMAX = 8
+ΔX = .1
+ΔY = .5
+XLEN = convert(Integer, floor((XMAX - XMIN)/ΔX) + 1)
+YLEN = convert(Integer, floor((YMAX - YMIN)/ΔY) + 1)
+
+@testset "1d-interpolations" begin
+    @testset "1d-regular-grids" begin
+        xs = XMIN:ΔX:XMAX
+        f(x) = log(x)
+        A = [f(x) for x in xs]
+        interp = LinearInterpolation(xs, A) # using convenience constructor
+        interp_full = extrapolate(scale(interpolate(A, BSpline(Linear()), OnGrid()), xs), Interpolations.Throw()) # using full constructor
+
+        @test typeof(interp) == typeof(interp_full)
+        @test interp[XMIN] ≈ f(XMIN)
+        @test interp[XMAX] ≈ f(XMAX)
+        @test interp[XMIN + ΔX] ≈ f(XMIN + ΔX)
+        @test interp[XMAX - ΔX] ≈ f(XMAX - ΔX)
+        @test interp[XMIN + ΔX / 2] ≈ f(XMIN + ΔX / 2) atol=.1
+        @test_throws BoundsError interp[XMIN - ΔX / 2]
+        @test_throws BoundsError interp[XMAX + ΔX / 2]
+    end
+
+    @testset "1d-regular-grids-cubic" begin
+        xs = XMIN:ΔX:XMAX
+        f(x) = log(x)
+        A = [f(x) for x in xs]
+        interp = CubicSplineInterpolation(xs, A)
+        interp_full = extrapolate(scale(interpolate(A, BSpline(Cubic(Line())), OnGrid()), xs), Interpolations.Throw())
+
+        @test typeof(interp) == typeof(interp_full)
+        @test interp[XMIN] ≈ f(XMIN)
+        @test interp[XMAX] ≈ f(XMAX)
+        @test interp[XMIN + ΔX] ≈ f(XMIN + ΔX)
+        @test interp[XMAX - ΔX] ≈ f(XMAX - ΔX)
+        @test interp[XMIN + ΔX / 2] ≈ f(XMIN + ΔX / 2) atol=.1
+        @test_throws BoundsError interp[XMIN - ΔX / 2]
+        @test_throws BoundsError interp[XMAX + ΔX / 2]
+    end
+
+    @testset "1d-irregular-grids" begin
+        xs = [x^2 for x in XMIN:ΔX:XMAX]
+        xmin = xs[1]
+        xmax = xs[XLEN]
+        f(x) = log(x)
+        A = [f(x) for x in xs]
+        interp = LinearInterpolation(xs, A)
+        interp_full = extrapolate(interpolate((xs, ), A, Gridded(Linear())), Interpolations.Throw())
+
+        @test typeof(interp) == typeof(interp_full)
+        @test interp[xmin] ≈ f(xmin)
+        @test interp[xmax] ≈ f(xmax)
+        @test interp[xs[2]] ≈ f(xs[2])
+        @test interp[xmin + ΔX / 2] ≈ f(xmin + ΔX / 2) atol=.1
+        @test_throws BoundsError interp[xmin - ΔX / 2]
+        @test_throws BoundsError interp[xmax + ΔX / 2]
+    end
+
+    @testset "1d-handling-extrapolation" begin
+        xs = XMIN:ΔX:XMAX
+        f(x) = log(x)
+        A = [f(x) for x in xs]
+        ΔA_l = A[2] - A[1]
+        ΔA_h = A[end] - A[end - 1]
+        x_lower = XMIN - ΔX
+        x_higher = XMAX + ΔX
+
+        extrap = LinearInterpolation(xs, A, extrapolation_bc = Interpolations.Linear())
+        extrap_full = extrapolate(scale(interpolate(A, BSpline(Linear()), OnGrid()), xs), Interpolations.Linear())
+
+        @test typeof(extrap) == typeof(extrap_full)
+        @test extrap[x_lower] ≈ A[1] - ΔA_l
+        @test extrap[x_higher] ≈ A[end] + ΔA_h
+    end
+end
+
+@testset "2d-interpolations" begin
+    @testset "2d-regular-grids" begin
+        xs = XMIN:ΔX:XMAX
+        ys = YMIN:ΔY:YMAX
+        f(x, y) = log(x+y)
+        A = [f(x,y) for x in xs, y in ys]
+        interp = LinearInterpolation((xs, ys), A)
+        interp_full = extrapolate(scale(interpolate(A, BSpline(Linear()), OnGrid()), xs, ys), Interpolations.Throw())
+
+        @test typeof(interp) == typeof(interp_full)
+        @test interp[XMIN,YMIN] ≈ f(XMIN,YMIN)
+        @test interp[XMIN,YMAX] ≈ f(XMIN,YMAX)
+        @test interp[XMAX,YMIN] ≈ f(XMAX,YMIN)
+        @test interp[XMAX,YMAX] ≈ f(XMAX,YMAX)
+        @test interp[XMIN + ΔX,YMIN] ≈ f(XMIN + ΔX,YMIN)
+        @test interp[XMIN,YMIN + ΔY] ≈ f(XMIN,YMIN + ΔY)
+        @test interp[XMIN + ΔX,YMIN + ΔY] ≈ f(XMIN + ΔX,YMIN + ΔY)
+        @test interp[XMIN + ΔX / 2,YMIN + ΔY / 2] ≈ f(XMIN + ΔX / 2,YMIN + ΔY / 2) atol=.1
+        @test_throws BoundsError interp[XMIN - ΔX / 2,YMIN - ΔY / 2]
+        @test_throws BoundsError interp[XMIN - ΔX / 2,YMIN + ΔY / 2]
+        @test_throws BoundsError interp[XMIN + ΔX / 2,YMIN - ΔY / 2]
+        @test_throws BoundsError interp[XMAX + ΔX / 2,YMAX + ΔY / 2]
+    end
+
+    @testset "2d-regular-grids-cubic" begin
+        xs = XMIN:ΔX:XMAX
+        ys = YMIN:ΔY:YMAX
+        f(x, y) = log(x+y)
+        A = [f(x,y) for x in xs, y in ys]
+        interp = CubicSplineInterpolation((xs, ys), A)
+        interp_full = extrapolate(scale(interpolate(A, BSpline(Cubic(Line())), OnGrid()), xs, ys), Interpolations.Throw())
+
+        @test typeof(interp) == typeof(interp_full)
+        @test interp[XMIN,YMIN] ≈ f(XMIN,YMIN)
+        @test interp[XMIN,YMAX] ≈ f(XMIN,YMAX)
+        @test interp[XMAX,YMIN] ≈ f(XMAX,YMIN)
+        @test interp[XMAX,YMAX] ≈ f(XMAX,YMAX)
+        @test interp[XMIN + ΔX,YMIN] ≈ f(XMIN + ΔX,YMIN)
+        @test interp[XMIN,YMIN + ΔY] ≈ f(XMIN,YMIN + ΔY)
+        @test interp[XMIN + ΔX,YMIN + ΔY] ≈ f(XMIN + ΔX,YMIN + ΔY)
+        @test interp[XMIN + ΔX / 2,YMIN + ΔY / 2] ≈ f(XMIN + ΔX / 2,YMIN + ΔY / 2) atol=.1
+        @test_throws BoundsError interp[XMIN - ΔX / 2,YMIN - ΔY / 2]
+        @test_throws BoundsError interp[XMIN - ΔX / 2,YMIN + ΔY / 2]
+        @test_throws BoundsError interp[XMIN + ΔX / 2,YMIN - ΔY / 2]
+        @test_throws BoundsError interp[XMAX + ΔX / 2,YMAX + ΔY / 2]
+    end
+
+    @testset "2d-irregular-grids" begin
+        xs = [x^2 for x in XMIN:ΔX:XMAX]
+        ys = [y^2 for y in YMIN:ΔY:YMAX]
+        xmin = xs[1]
+        xmax = xs[XLEN]
+        ymin = ys[1]
+        ymax = ys[YLEN]
+        f(x, y) = log(x+y)
+        A = [f(x,y) for x in xs, y in ys]
+        interp = LinearInterpolation((xs, ys), A)
+        interp_full = extrapolate(interpolate((xs, ys), A, Gridded(Linear())), Interpolations.Throw())
+
+        @test typeof(interp) == typeof(interp_full)
+        @test interp[xmin,ymin] ≈ f(xmin,ymin)
+        @test interp[xmin,ymax] ≈ f(xmin,ymax)
+        @test interp[xmax,ymin] ≈ f(xmax,ymin)
+        @test interp[xmax,ymax] ≈ f(xmax,ymax)
+        @test interp[xs[2],ymin] ≈ f(xs[2],ymin)
+        @test interp[xmin,ys[2]] ≈ f(xmin,ys[2])
+        @test interp[xs[2],ys[2]] ≈ f(xs[2],ys[2])
+        @test interp[xmin + ΔX / 2,ymin + ΔY / 2] ≈ f(xmin + ΔX / 2,ymin + ΔY / 2) atol=.1
+        @test_throws BoundsError interp[xmin - ΔX / 2,ymin - ΔY / 2]
+        @test_throws BoundsError interp[xmin - ΔX / 2,ymin + ΔY / 2]
+        @test_throws BoundsError interp[xmin + ΔX / 2,ymin - ΔY / 2]
+        @test_throws BoundsError interp[xmax + ΔX / 2,ymax + ΔY / 2]
+    end
+
+    @testset "2d-handling-extrapolation" begin
+        xs = XMIN:ΔX:XMAX
+        ys = YMIN:ΔY:YMAX
+        f(x, y) = log(x+y)
+        A = [f(x,y) for x in xs, y in ys]
+        ΔA_l = A[2, 1] - A[1, 1]
+        ΔA_h = A[end, end] - A[end - 1, end]
+        x_lower = XMIN - ΔX
+        x_higher = XMAX + ΔX
+        y_lower = YMIN - ΔY
+        y_higher = YMAX + ΔY
+
+        extrap = LinearInterpolation((xs, ys), A, extrapolation_bc = (Interpolations.Linear(), Interpolations.Flat()))
+        extrap_full = extrapolate(scale(interpolate(A, BSpline(Linear()), OnGrid()), xs, ys), (Interpolations.Linear(), Interpolations.Flat()))
+
+        @test typeof(extrap) == typeof(extrap_full)
+        @test extrap[x_lower, y_lower] ≈ A[1, 1] - ΔA_l
+        @test extrap[x_higher, y_higher] ≈ A[end, end] + ΔA_h
+    end
+end
+
+end

test/runtests.jl

@@ -28,6 +28,7 @@ include("typing.jl")
 include("issues/runtests.jl")
 
 include("io.jl")
+include("convenience-constructors.jl")
 include("readme-examples.jl")
 
 end