diff --git a/README.md b/README.md index 045d21af..ee90d04c 100644 --- a/README.md +++ b/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 diff --git a/src/Interpolations.jl b/src/Interpolations.jl index 5f0e2f1b..bb16ea13 100644 --- a/src/Interpolations.jl +++ b/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 diff --git a/src/convenience-constructors.jl b/src/convenience-constructors.jl new file mode 100644 index 00000000..8abbb6c4 --- /dev/null +++ b/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) diff --git a/test/convenience-constructors.jl b/test/convenience-constructors.jl new file mode 100644 index 00000000..bb97b763 --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 0f52b74d..ac2a683b 100644 --- a/test/runtests.jl +++ b/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