Skip to content

Commit 328d0c1

Browse files
committed
Update the README and other docs to the new API
1 parent c1124c1 commit 328d0c1

File tree

5 files changed

+143
-6077
lines changed

5 files changed

+143
-6077
lines changed

README.md

+43-72
Original file line numberDiff line numberDiff line change
@@ -5,17 +5,16 @@
55
[![Interpolations](http://pkg.julialang.org/badges/Interpolations_0.5.svg)](http://pkg.julialang.org/?pkg=Interpolations)
66

77
This package implements a variety of interpolation schemes for the
8-
Julia langauge. It has the goals of ease-of-use, broad algorithmic
8+
Julia language. It has the goals of ease-of-use, broad algorithmic
99
support, and exceptional performance.
1010

11-
This package is still relatively new. Currently its support is best
11+
Currently this package's support is best
1212
for [B-splines](https://en.wikipedia.org/wiki/B-spline) and also
1313
supports irregular grids. However, the API has been designed with
1414
intent to support more options. Pull-requests are more than welcome!
1515
It should be noted that the API may continue to evolve over time.
1616

1717
Other interpolation packages for Julia include:
18-
- [Grid.jl](https://github.com/timholy/Grid.jl) (the predecessor of this package)
1918
- [Dierckx.jl](https://github.com/kbarbary/Dierckx.jl)
2019
- [GridInterpolations.jl](https://github.com/sisl/GridInterpolations.jl)
2120
- [ApproXD.jl](https://github.com/floswald/ApproXD.jl)
@@ -42,38 +41,38 @@ from the Julia REPL.
4241
Note: the current version of `Interpolations` supports interpolation evaluation using index calls `[]`, but this feature will be deprecated in future. We highly recommend function calls with `()` as follows.
4342

4443
Given an `AbstractArray` `A`, construct an "interpolation object" `itp` as
45-
```jl
44+
```julia
4645
itp = interpolate(A, options...)
4746
```
4847
where `options...` (discussed below) controls the type of
4948
interpolation you want to perform. This syntax assumes that the
5049
samples in `A` are equally-spaced.
5150

5251
To evaluate the interpolation at position `(x, y, ...)`, simply do
53-
```jl
52+
```julia
5453
v = itp(x, y, ...)
5554
```
5655

5756
Some interpolation objects support computation of the gradient, which
5857
can be obtained as
59-
```jl
58+
```julia
6059
g = gradient(itp, x, y, ...)
6160
```
6261
or, if you're evaluating the gradient repeatedly, a somewhat more
6362
efficient option is
64-
```jl
63+
```julia
6564
gradient!(g, itp, x, y, ...)
6665
```
6766
where `g` is a pre-allocated vector.
6867

6968
Some interpolation objects support computation of the hessian, which
7069
can be obtained as
71-
```jl
70+
```julia
7271
h = hessian(itp, x, y, ...)
7372
```
7473
or, if you're evaluating the hessian repeatedly, a somewhat more
7574
efficient option is
76-
```jl
75+
```julia
7776
hessian!(h, itp, x, y, ...)
7877
```
7978
where `h` is a pre-allocated matrix.
@@ -85,25 +84,28 @@ and `Rational`, but also multi-valued types like `RGB` color vectors.
8584
Positions `(x, y, ...)` are n-tuples of numbers. Typically these will
8685
be real-valued (not necessarily integer-valued), but can also be of types
8786
such as [DualNumbers](https://github.com/JuliaDiff/DualNumbers.jl) if
88-
you want to verify the computed value of gradients. You can also use
87+
you want to verify the computed value of gradients.
88+
(Alternatively, verify gradients using [ForwardDiff](https://github.com/JuliaDiff/ForwardDiff.jl).)
89+
You can also use
8990
Julia's iterator objects, e.g.,
9091

91-
```jl
92+
```julia
9293
function ongrid!(dest, itp)
93-
for I in CartesianRange(size(itp))
94+
for I in CartesianIndices(itp)
9495
dest[I] = itp(I)
9596
end
9697
end
9798
```
9899
would store the on-grid value at each grid point of `itp` in the output `dest`.
99100
Finally, courtesy of Julia's indexing rules, you can also use
100-
```jl
101-
fine = itp(linspace(1,10,1001), linspace(1,15,201))
101+
```julia
102+
fine = itp(range(1,stop=10,length=1001), range(1,stop=15,length=201))
102103
```
103104

104105
### Quickstart guide
106+
105107
For linear and cubic spline interpolations, `LinearInterpolation` and `CubicSplineInterpolation` can be used to create interpolation objects handily:
106-
```jl
108+
```julia
107109
f(x) = log(x)
108110
xs = 1:0.2:5
109111
A = [f(x) for x in xs]
@@ -119,7 +121,7 @@ interp_cubic(3) # exactly log(3)
119121
interp_cubic(3.1) # approximately log(3.1)
120122
```
121123
which support multidimensional data as well:
122-
```jl
124+
```julia
123125
f(x,y) = log(x+y)
124126
xs = 1:0.2:5
125127
ys = 2:0.1:5
@@ -137,19 +139,19 @@ interp_cubic(3.1, 2.1) # approximately log(3.1 + 2.1)
137139
```
138140
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.
139141
Interested users can specify boundary conditions by providing an extra parameter for `extrapolation_bc`:
140-
```jl
142+
```julia
141143
f(x) = log(x)
142144
xs = 1:0.2:5
143145
A = [f(x) for x in xs]
144146

145147
# extrapolation with linear boundary conditions
146-
extrap = LinearInterpolation(xs, A, extrapolation_bc = Interpolations.Linear())
148+
extrap = LinearInterpolation(xs, A, extrapolation_bc = Line())
147149

148150
@test extrap(1 - 0.2) # ≈ f(1) - (f(1.2) - f(1))
149151
@test extrap(5 + 0.2) # ≈ f(5) + (f(5) - f(4.8))
150152
```
151153
Irregular grids are supported as well; note that presently only `LinearInterpolation` supports irregular grids.
152-
```jl
154+
```julia
153155
xs = [x^2 for x = 1:0.2:5]
154156
A = [f(x) for x in xs]
155157

@@ -163,34 +165,34 @@ interp_linear(1.05) # approximately log(1.05)
163165

164166
### BSplines
165167

166-
The interpolation type is described in terms of *degree*, *grid behavior* and, if necessary, *boundary conditions*. There are currently three degrees available: `Constant`, `Linear`, `Quadratic`, and `Cubic` corresponding to B-splines of degree 0, 1, 2, and 3 respectively.
167-
168-
You also have to specify what *grid representation* you want. There are currently two choices: `OnGrid`, in which the supplied data points are assumed to lie *on* the boundaries of the interpolation interval, and `OnCell` in which the data points are assumed to lie on half-intervals between cell boundaries.
168+
The interpolation type is described in terms of *degree* and, if necessary, *boundary conditions*. There are currently three degrees available: `Constant`, `Linear`, `Quadratic`, and `Cubic` corresponding to B-splines of degree 0, 1, 2, and 3 respectively.
169169

170170
B-splines of quadratic or higher degree require solving an equation system to obtain the interpolation coefficients, and for that you must specify a *boundary condition* that is applied to close the system. The following boundary conditions are implemented: `Flat`, `Line` (alternatively, `Natural`), `Free`, `Periodic` and `Reflect`; their mathematical implications are described in detail in the pdf document under `/doc/latex`.
171+
When specifying these boundary conditions you also have to specify whether they apply at the edge grid point (`OnGrid()`)
172+
or beyond the edge point halfway to the next (fictitious) grid point (`OnCell()`).
171173

172174
Some examples:
173-
```jl
175+
```julia
174176
# Nearest-neighbor interpolation
175-
itp = interpolate(a, BSpline(Constant()), OnCell())
177+
itp = interpolate(a, BSpline(Constant()))
176178
v = itp(5.4) # returns a[5]
177179

178180
# (Multi)linear interpolation
179-
itp = interpolate(A, BSpline(Linear()), OnGrid())
181+
itp = interpolate(A, BSpline(Linear()))
180182
v = itp(3.2, 4.1) # returns 0.9*(0.8*A[3,4]+0.2*A[4,4]) + 0.1*(0.8*A[3,5]+0.2*A[4,5])
181183

182184
# Quadratic interpolation with reflecting boundary conditions
183185
# Quadratic is the lowest order that has continuous gradient
184-
itp = interpolate(A, BSpline(Quadratic(Reflect())), OnCell())
186+
itp = interpolate(A, BSpline(Quadratic(Reflect(OnCell()))))
185187

186188
# Linear interpolation in the first dimension, and no interpolation (just lookup) in the second
187-
itp = interpolate(A, (BSpline(Linear()), NoInterp()), OnGrid())
189+
itp = interpolate(A, (BSpline(Linear()), NoInterp()))
188190
v = itp(3.65, 5) # returns 0.35*A[3,5] + 0.65*A[4,5]
189191
```
190192
There are more options available, for example:
191-
```jl
193+
```julia
192194
# In-place interpolation
193-
itp = interpolate!(A, BSpline(Quadratic(InPlace())), OnCell())
195+
itp = interpolate!(A, BSpline(Quadratic(InPlace(OnCell()))))
194196
```
195197
which destroys the input `A` but also does not need to allocate as much memory.
196198

@@ -199,22 +201,22 @@ which destroys the input `A` but also does not need to allocate as much memory.
199201
BSplines assume your data is uniformly spaced on the grid `1:N`, or its multidimensional equivalent. If you have data of the form `[f(x) for x in A]`, you need to tell Interpolations about the grid `A`. If `A` is not uniformly spaced, you must use gridded interpolation described below. However, if `A` is a collection of ranges or linspaces, you can use scaled BSplines. This is more efficient because the gridded algorithm does not exploit the uniform spacing. Scaled BSplines can also be used with any spline degree available for BSplines, while gridded interpolation does not currently support quadratic or cubic splines.
200202

201203
Some examples,
202-
```jl
204+
```julia
203205
A_x = 1.:2.:40.
204206
A = [log(x) for x in A_x]
205-
itp = interpolate(A, BSpline(Cubic(Line())), OnGrid())
207+
itp = interpolate(A, BSpline(Cubic(Line(OnGrid()))))
206208
sitp = scale(itp, A_x)
207209
sitp(3.) # exactly log(3.)
208210
sitp(3.5) # approximately log(3.5)
209211
```
210212

211213
For multidimensional uniformly spaced grids
212-
```jl
214+
```julia
213215
A_x1 = 1:.1:10
214216
A_x2 = 1:.5:20
215217
f(x1, x2) = log(x1+x2)
216218
A = [f(x1,x2) for x1 in A_x1, x2 in A_x2]
217-
itp = interpolate(A, BSpline(Cubic(Line())), OnGrid())
219+
itp = interpolate(A, BSpline(Cubic(Line(OnGrid()))))
218220
sitp = scale(itp, A_x1, A_x2)
219221
sitp(5., 10.) # exactly log(5 + 10)
220222
sitp(5.6, 7.1) # approximately log(5.6 + 7.1)
@@ -227,7 +229,7 @@ are all `OnGrid`). As such one must specify a set of coordinate arrays
227229
defining the knots of the array.
228230

229231
In 1D
230-
```jl
232+
```julia
231233
A = rand(20)
232234
A_x = collect(1.0:2.0:40.0)
233235
knots = (A_x,)
@@ -236,34 +238,34 @@ itp(2.0)
236238
```
237239

238240
The spacing between adjacent samples need not be constant, you can use the syntax
239-
```jl
241+
```julia
240242
itp = interpolate(knots, A, options...)
241243
```
242244
where `knots = (xknots, yknots, ...)` to specify the positions along
243245
each axis at which the array `A` is sampled for arbitrary ("rectangular") samplings.
244246

245247
For example:
246-
```jl
248+
```julia
247249
A = rand(8,20)
248250
knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
249251
itp = interpolate(knots, A, Gridded(Linear()))
250252
itp(4,1.2) # approximately A[2,6]
251253
```
252254
One may also mix modes, by specifying a mode vector in the form of an explicit tuple:
253-
```jl
255+
```julia
254256
itp = interpolate(knots, A, (Gridded(Linear()),Gridded(Constant())))
255257
```
256258

257259
Presently there are only three modes for gridded:
258-
```jl
260+
```julia
259261
Gridded(Linear())
260262
```
261263
whereby a linear interpolation is applied between knots,
262-
```jl
264+
```julia
263265
Gridded(Constant())
264266
```
265267
whereby nearest neighbor interpolation is used on the applied axis,
266-
```jl
268+
```julia
267269
NoInterp
268270
```
269271
whereby the coordinate of the selected input vector MUST be located on a grid point. Requests for off grid
@@ -281,7 +283,7 @@ x = sin.(2π*t)
281283
y = cos.(2π*t)
282284
A = hcat(x,y)
283285

284-
itp = scale(interpolate(A, (BSpline(Cubic(Natural())), NoInterp()), OnGrid()), t, 1:2)
286+
itp = scale(interpolate(A, (BSpline(Cubic(Natural(OnGrid()))), NoInterp())), t, 1:2)
285287

286288
tfine = 0:.01:1
287289
xs, ys = [itp(t,1) for t in tfine], [itp(t,2) for t in tfine]
@@ -366,37 +368,6 @@ they ran more than 20 seconds (far longer than any other test). Both
366368
performed much better in 2d, interestingly. You can see that
367369
Interpolations wins in every case, sometimes by a very large margin.
368370

369-
## Transitioning from Grid.jl
370-
371-
Instead of using
372-
```julia
373-
yi = InterpGrid(y, BCreflect, InterpQuadratic)
374-
```
375-
you should use
376-
```julia
377-
yi = interpolate(y, BSpline(Quadratic(Reflect())), OnCell())
378-
```
379-
380-
In general, here are the closest mappings:
381-
382-
| Grid | Interpolations |
383-
|:-----------------:|:------------------------------------------:|
384-
| `InterpNearest` | `Constant` |
385-
| `InterpLinear` | `Linear` |
386-
| `InterpQuadratic` | `Quadratic` |
387-
| `InterpCubic` | `Cubic` |
388-
| | |
389-
| `BCnil` | `extrapolate(itp, Interpolations.Throw())` |
390-
| `BCnan` | `extrapolate(itp, NaN)` |
391-
| `BCna` | `extrapolate(itp, NaN)` |
392-
| `BCreflect` | `interpolate` with `Reflect()` |
393-
| `BCperiodic` | `interpolate` with `Periodic()` |
394-
| `BCnearest` | `interpolate` with `Flat()` |
395-
| `BCfill` | `extrapolate` with value |
396-
| | |
397-
| odd orders | `OnGrid()` |
398-
| even orders | `OnCell()` |
399-
400371

401372
## Contributing
402373

0 commit comments

Comments
 (0)