Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Generalize signatures of convenience constructors #241

Merged
merged 3 commits into from
Sep 27, 2018
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
147 changes: 74 additions & 73 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
38 changes: 32 additions & 6 deletions src/convenience-constructors.jl
Original file line number Diff line number Diff line change
@@ -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
11 changes: 11 additions & 0 deletions test/convenience-constructors.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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