diff --git a/README.md b/README.md index 8a078d60..f6a037fc 100644 --- a/README.md +++ b/README.md @@ -58,24 +58,22 @@ v = itp(x, y, ...) Some interpolation objects support computation of the gradient, which can be obtained as ```julia -g = gradient(itp, x, y, ...) +g = Interpolations.gradient(itp, x, y, ...) ``` -or, if you're evaluating the gradient repeatedly, a somewhat more -efficient option is +or as ```julia -gradient!(g, itp, x, y, ...) +Interpolations.gradient!(g, itp, x, y, ...) ``` where `g` is a pre-allocated vector. Some interpolation objects support computation of the hessian, which can be obtained as ```julia -h = hessian(itp, x, y, ...) +h = Interpolations.hessian(itp, x, y, ...) ``` -or, if you're evaluating the hessian repeatedly, a somewhat more -efficient option is +or ```julia -hessian!(h, itp, x, y, ...) +Interpolations.hessian!(h, itp, x, y, ...) ``` where `h` is a pre-allocated matrix. @@ -104,71 +102,7 @@ Finally, courtesy of Julia's indexing rules, you can also use fine = itp(range(1,stop=10,length=1001), range(1,stop=15,length=201)) ``` -### Quickstart guide - -For linear and cubic spline interpolations, `LinearInterpolation` and `CubicSplineInterpolation` can be used to create interpolation objects handily: -```julia -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: -```julia -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`: -```julia -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 = 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] -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) -``` +There is also an abbreviated notion [described below](#convenience-notation). ## Control of interpolation algorithm @@ -323,6 +257,73 @@ etpf = extrapolate(itp, Flat()) # gives 1 on the left edge and 7 on the right etp0 = extrapolate(itp, 0) # gives 0 everywhere outside [1,7] ``` +## Convenience notation + +For linear and cubic spline interpolations, `LinearInterpolation` and `CubicSplineInterpolation` +can be used to create interpolating and extrapolating objects handily: +```julia +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: +```julia +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 the 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`: +```julia +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 = 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] +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) +``` + ## Performance shootout In the `perf` directory, you can find a script that tests diff --git a/src/convenience-constructors.jl b/src/convenience-constructors.jl index 61123911..2f1cd3d0 100644 --- a/src/convenience-constructors.jl +++ b/src/convenience-constructors.jl @@ -1,10 +1,36 @@ # convenience copnstructors for linear / cubic spline interpolations # 1D version -LinearInterpolation(range::T, vs; extrapolation_bc = Throw()) where {T <: AbstractRange} = extrapolate(scale(interpolate(vs, BSpline(Linear())), range), extrapolation_bc) -LinearInterpolation(range::T, vs; extrapolation_bc = Throw()) where {T <: AbstractArray} = extrapolate(interpolate((range, ), vs, Gridded(Linear())), extrapolation_bc) -CubicSplineInterpolation(range::T, vs; bc = Line(OnGrid()), extrapolation_bc = Throw()) where {T <: AbstractRange} = extrapolate(scale(interpolate(vs, BSpline(Cubic(bc))), range), extrapolation_bc) +LinearInterpolation(range::AbstractRange, vs::AbstractVector; extrapolation_bc = Throw()) = + extrapolate(scale(interpolate(vs, BSpline(Linear())), range), extrapolation_bc) +LinearInterpolation(range::AbstractVector, vs::AbstractVector; extrapolation_bc = Throw()) = + extrapolate(interpolate((range, ), vs, Gridded(Linear())), extrapolation_bc) +CubicSplineInterpolation(range::AbstractRange, vs::AbstractVector; + bc = Line(OnGrid()), extrapolation_bc = Throw()) = + extrapolate(scale(interpolate(vs, BSpline(Cubic(bc))), range), extrapolation_bc) # multivariate versions -LinearInterpolation(ranges::NTuple{N,T}, vs; extrapolation_bc = Throw()) where {N,T <: AbstractRange} = extrapolate(scale(interpolate(vs, BSpline(Linear())), ranges...), extrapolation_bc) -LinearInterpolation(ranges::NTuple{N,T}, vs; extrapolation_bc = Throw()) where {N,T <: AbstractArray} = extrapolate(interpolate(ranges, vs, Gridded(Linear())), extrapolation_bc) -CubicSplineInterpolation(ranges::NTuple{N,T}, vs; bc = Line(OnGrid()), extrapolation_bc = Throw()) where {N,T <: AbstractRange} = extrapolate(scale(interpolate(vs, BSpline(Cubic(bc))), ranges...), extrapolation_bc) +LinearInterpolation(ranges::NTuple{N,AbstractRange}, vs::AbstractArray{T,N}; + extrapolation_bc = Throw()) where {N,T} = + extrapolate(scale(interpolate(vs, BSpline(Linear())), ranges...), extrapolation_bc) +LinearInterpolation(ranges::NTuple{N,AbstractVector}, vs::AbstractArray{T,N}; + extrapolation_bc = Throw()) where {N,T} = + extrapolate(interpolate(ranges, vs, Gridded(Linear())), extrapolation_bc) +CubicSplineInterpolation(ranges::NTuple{N,AbstractRange}, vs::AbstractArray{T,N}; + bc = Line(OnGrid()), extrapolation_bc = Throw()) where {N,T} = + extrapolate(scale(interpolate(vs, BSpline(Cubic(bc))), ranges...), extrapolation_bc) + +""" + etp = LinearInterpolation(knots, A; extrapolation_bc=Throw()) + +A shorthand for `extrapolate(interpolate(knots, A, scheme), extrapolation_bc)`, +where `scheme` is either `BSpline(Linear())` or `Gridded(Linear())` depending on whether +`knots` are ranges or vectors. +""" +LinearInterpolation + +""" + etp = CubicSplineInterpolation(knots, A; bc=Line(OnGrid()), extrapolation_bc=Throw()) + +A shorthand for `extrapolate(interpolate(knots, A, BSpline(Cubic(bc))), extrapolation_bc)`. +""" +CubicSplineInterpolation diff --git a/test/convenience-constructors.jl b/test/convenience-constructors.jl index 70a30107..5a79b5fc 100644 --- a/test/convenience-constructors.jl +++ b/test/convenience-constructors.jl @@ -178,6 +178,17 @@ end @test extrap(x_lower, y_lower) ≈ A[1, 1] - ΔA_l @test extrap(x_higher, y_higher) ≈ A[end, end] + ΔA_h end + + @testset "issue #230" begin # at least, I think this is what issue #230 is really about + f(x,y) = log(x+y) + xs = 1:5 + ys = 2:0.1:5 + A = [f(x,y) for x in xs, y in ys] + itp = LinearInterpolation((xs, ys), A) + for (i, j) in zip(Iterators.product(xs, ys), eachindex(A)) + @test itp(i...) ≈ A[j] + end + end end end