Skip to content

Commit 2633025

Browse files
chiyahnsglyon
chiyahn
authored and
sglyon
committed
Convenience constructors (#214)
* convenience constructors added * parameter for bc added in convenience constructors * convenience constructors dedicated for 1D cases added * unit tests for type consistency with constructors from full API * unit tests for convenience constructors with valid extrapolation * quickstart guide with convenience constructors in README.md * optional parameter for bc added in CubicSplineInterpolation
1 parent 4597d98 commit 2633025

File tree

5 files changed

+256
-1
lines changed

5 files changed

+256
-1
lines changed

‎README.md

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -101,6 +101,63 @@ Finally, courtesy of Julia's indexing rules, you can also use
101101
fine = itp(linspace(1,10,1001), linspace(1,15,201))
102102
```
103103

104+
### Quickstart guide
105+
For linear and cubic spline interpolations, `LinearInterpolation` and `CubicSplineInterpolation` can be used to create interpolation objects handily:
106+
```jl
107+
f(x) = log(x)
108+
xs = 1:0.2:5
109+
A = [f(x) for x in xs]
110+
111+
# linear interpolation
112+
interp_linear = LinearInterpolation(xs, A)
113+
interp_linear[3] # exactly log(3)
114+
interp_linear[3.1] # approximately log(3.1)
115+
116+
# cubic spline interpolation
117+
interp_cubic = CubicSplineInterpolation(xs, A)
118+
interp_cubic[3] # exactly log(3)
119+
interp_cubic[3.1] # approximately log(3.1)
120+
```
121+
which support multidimensional data as well:
122+
```jl
123+
f(x,y) = log(x+y)
124+
xs = 1:0.2:5
125+
ys = 2:0.1:5
126+
A = [f(x+y) for x in xs, y in ys]
127+
128+
# linear interpolation
129+
interp_linear = LinearInterpolation((xs, ys), A)
130+
interp_linear[3, 2] # exactly log(3 + 2)
131+
interp_linear[3.1, 2.1] # approximately log(3.1 + 2.1)
132+
133+
# cubic spline interpolation
134+
interp_cubic = CubicSplineInterpolation((xs, ys), A)
135+
interp_cubic[3, 2] # exactly log(3 + 2)
136+
interp_cubic[3.1, 2.1] # approximately log(3.1 + 2.1)
137+
```
138+
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.
139+
Interested users can specify boundary conditions by providing an extra parameter for `extrapolation_bc`:
140+
```jl
141+
f(x) = log(x)
142+
xs = 1:0.2:5
143+
A = [f(x) for x in xs]
144+
145+
# extrapolation with linear boundary conditions
146+
extrap = LinearInterpolation(xs, A, extrapolation_bc = Interpolations.Linear())
147+
148+
@test extrap[1 - 0.2] # ≈ f(1) - (f(1.2) - f(1))
149+
@test extrap[5 + 0.2] # ≈ f(5) + (f(5) - f(4.8))
150+
```
151+
Irregular grids are supported as well; note that presently only `LinearInterpolation` supports irregular grids.
152+
```jl
153+
xs = [x^2 for x = 1:0.2:5]
154+
A = [f(x) for x in xs]
155+
156+
# linear interpolation
157+
interp_linear = LinearInterpolation(xs, A)
158+
interp_linear[1] # exactly log(1)
159+
interp_linear[1.05] # approximately log(1.05)
160+
```
104161

105162
## Control of interpolation algorithm
106163

‎src/Interpolations.jl

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,10 @@ export
2727
Reflect,
2828
Natural,
2929
InPlace,
30-
InPlaceQ
30+
InPlaceQ,
31+
32+
LinearInterpolation,
33+
CubicSplineInterpolation
3134

3235
# see the following files for further exports:
3336
# b-splines/b-splines.jl
@@ -114,5 +117,6 @@ include("extrapolation/extrapolation.jl")
114117
include("scaling/scaling.jl")
115118
include("utils.jl")
116119
include("io.jl")
120+
include("convenience-constructors.jl")
117121

118122
end # module

‎src/convenience-constructors.jl

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
# convenience copnstructors for linear / cubic spline interpolations
2+
# 1D version
3+
LinearInterpolation(range::T, vs; extrapolation_bc = Interpolations.Throw()) where {T <: Range} = extrapolate(scale(interpolate(vs, BSpline(Linear()), OnGrid()), range), extrapolation_bc)
4+
LinearInterpolation(range::T, vs; extrapolation_bc = Interpolations.Throw()) where {T <: AbstractArray} = extrapolate(interpolate((range, ), vs, Gridded(Linear())), extrapolation_bc)
5+
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)
6+
7+
# multivariate versions
8+
LinearInterpolation(ranges::NTuple{N,T}, vs; extrapolation_bc = Interpolations.Throw()) where {N,T <: Range} = extrapolate(scale(interpolate(vs, BSpline(Linear()), OnGrid()), ranges...), extrapolation_bc)
9+
LinearInterpolation(ranges::NTuple{N,T}, vs; extrapolation_bc = Interpolations.Throw()) where {N,T <: AbstractArray} = extrapolate(interpolate(ranges, vs, Gridded(Linear())), extrapolation_bc)
10+
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

Lines changed: 183 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,183 @@
1+
module ConvenienceConstructorTests
2+
3+
using Interpolations
4+
using Base.Test
5+
using Base.Cartesian
6+
7+
# unit test setup
8+
XMIN = 2
9+
XMAX = 10
10+
YMIN = 1
11+
YMAX = 8
12+
ΔX = .1
13+
ΔY = .5
14+
XLEN = convert(Integer, floor((XMAX - XMIN)/ΔX) + 1)
15+
YLEN = convert(Integer, floor((YMAX - YMIN)/ΔY) + 1)
16+
17+
@testset "1d-interpolations" begin
18+
@testset "1d-regular-grids" begin
19+
xs = XMIN:ΔX:XMAX
20+
f(x) = log(x)
21+
A = [f(x) for x in xs]
22+
interp = LinearInterpolation(xs, A) # using convenience constructor
23+
interp_full = extrapolate(scale(interpolate(A, BSpline(Linear()), OnGrid()), xs), Interpolations.Throw()) # using full constructor
24+
25+
@test typeof(interp) == typeof(interp_full)
26+
@test interp[XMIN] ≈ f(XMIN)
27+
@test interp[XMAX] ≈ f(XMAX)
28+
@test interp[XMIN + ΔX] ≈ f(XMIN + ΔX)
29+
@test interp[XMAX - ΔX] ≈ f(XMAX - ΔX)
30+
@test interp[XMIN + ΔX / 2] ≈ f(XMIN + ΔX / 2) atol=.1
31+
@test_throws BoundsError interp[XMIN - ΔX / 2]
32+
@test_throws BoundsError interp[XMAX + ΔX / 2]
33+
end
34+
35+
@testset "1d-regular-grids-cubic" begin
36+
xs = XMIN:ΔX:XMAX
37+
f(x) = log(x)
38+
A = [f(x) for x in xs]
39+
interp = CubicSplineInterpolation(xs, A)
40+
interp_full = extrapolate(scale(interpolate(A, BSpline(Cubic(Line())), OnGrid()), xs), Interpolations.Throw())
41+
42+
@test typeof(interp) == typeof(interp_full)
43+
@test interp[XMIN] ≈ f(XMIN)
44+
@test interp[XMAX] ≈ f(XMAX)
45+
@test interp[XMIN + ΔX] ≈ f(XMIN + ΔX)
46+
@test interp[XMAX - ΔX] ≈ f(XMAX - ΔX)
47+
@test interp[XMIN + ΔX / 2] ≈ f(XMIN + ΔX / 2) atol=.1
48+
@test_throws BoundsError interp[XMIN - ΔX / 2]
49+
@test_throws BoundsError interp[XMAX + ΔX / 2]
50+
end
51+
52+
@testset "1d-irregular-grids" begin
53+
xs = [x^2 for x in XMIN:ΔX:XMAX]
54+
xmin = xs[1]
55+
xmax = xs[XLEN]
56+
f(x) = log(x)
57+
A = [f(x) for x in xs]
58+
interp = LinearInterpolation(xs, A)
59+
interp_full = extrapolate(interpolate((xs, ), A, Gridded(Linear())), Interpolations.Throw())
60+
61+
@test typeof(interp) == typeof(interp_full)
62+
@test interp[xmin] ≈ f(xmin)
63+
@test interp[xmax] ≈ f(xmax)
64+
@test interp[xs[2]] ≈ f(xs[2])
65+
@test interp[xmin + ΔX / 2] ≈ f(xmin + ΔX / 2) atol=.1
66+
@test_throws BoundsError interp[xmin - ΔX / 2]
67+
@test_throws BoundsError interp[xmax + ΔX / 2]
68+
end
69+
70+
@testset "1d-handling-extrapolation" begin
71+
xs = XMIN:ΔX:XMAX
72+
f(x) = log(x)
73+
A = [f(x) for x in xs]
74+
ΔA_l = A[2] - A[1]
75+
ΔA_h = A[end] - A[end - 1]
76+
x_lower = XMIN - ΔX
77+
x_higher = XMAX + ΔX
78+
79+
extrap = LinearInterpolation(xs, A, extrapolation_bc = Interpolations.Linear())
80+
extrap_full = extrapolate(scale(interpolate(A, BSpline(Linear()), OnGrid()), xs), Interpolations.Linear())
81+
82+
@test typeof(extrap) == typeof(extrap_full)
83+
@test extrap[x_lower] ≈ A[1] - ΔA_l
84+
@test extrap[x_higher] ≈ A[end] + ΔA_h
85+
end
86+
end
87+
88+
@testset "2d-interpolations" begin
89+
@testset "2d-regular-grids" begin
90+
xs = XMIN:ΔX:XMAX
91+
ys = YMIN:ΔY:YMAX
92+
f(x, y) = log(x+y)
93+
A = [f(x,y) for x in xs, y in ys]
94+
interp = LinearInterpolation((xs, ys), A)
95+
interp_full = extrapolate(scale(interpolate(A, BSpline(Linear()), OnGrid()), xs, ys), Interpolations.Throw())
96+
97+
@test typeof(interp) == typeof(interp_full)
98+
@test interp[XMIN,YMIN] ≈ f(XMIN,YMIN)
99+
@test interp[XMIN,YMAX] ≈ f(XMIN,YMAX)
100+
@test interp[XMAX,YMIN] ≈ f(XMAX,YMIN)
101+
@test interp[XMAX,YMAX] ≈ f(XMAX,YMAX)
102+
@test interp[XMIN + ΔX,YMIN] ≈ f(XMIN + ΔX,YMIN)
103+
@test interp[XMIN,YMIN + ΔY] ≈ f(XMIN,YMIN + ΔY)
104+
@test interp[XMIN + ΔX,YMIN + ΔY] ≈ f(XMIN + ΔX,YMIN + ΔY)
105+
@test interp[XMIN + ΔX / 2,YMIN + ΔY / 2] ≈ f(XMIN + ΔX / 2,YMIN + ΔY / 2) atol=.1
106+
@test_throws BoundsError interp[XMIN - ΔX / 2,YMIN - ΔY / 2]
107+
@test_throws BoundsError interp[XMIN - ΔX / 2,YMIN + ΔY / 2]
108+
@test_throws BoundsError interp[XMIN + ΔX / 2,YMIN - ΔY / 2]
109+
@test_throws BoundsError interp[XMAX + ΔX / 2,YMAX + ΔY / 2]
110+
end
111+
112+
@testset "2d-regular-grids-cubic" begin
113+
xs = XMIN:ΔX:XMAX
114+
ys = YMIN:ΔY:YMAX
115+
f(x, y) = log(x+y)
116+
A = [f(x,y) for x in xs, y in ys]
117+
interp = CubicSplineInterpolation((xs, ys), A)
118+
interp_full = extrapolate(scale(interpolate(A, BSpline(Cubic(Line())), OnGrid()), xs, ys), Interpolations.Throw())
119+
120+
@test typeof(interp) == typeof(interp_full)
121+
@test interp[XMIN,YMIN] ≈ f(XMIN,YMIN)
122+
@test interp[XMIN,YMAX] ≈ f(XMIN,YMAX)
123+
@test interp[XMAX,YMIN] ≈ f(XMAX,YMIN)
124+
@test interp[XMAX,YMAX] ≈ f(XMAX,YMAX)
125+
@test interp[XMIN + ΔX,YMIN] ≈ f(XMIN + ΔX,YMIN)
126+
@test interp[XMIN,YMIN + ΔY] ≈ f(XMIN,YMIN + ΔY)
127+
@test interp[XMIN + ΔX,YMIN + ΔY] ≈ f(XMIN + ΔX,YMIN + ΔY)
128+
@test interp[XMIN + ΔX / 2,YMIN + ΔY / 2] ≈ f(XMIN + ΔX / 2,YMIN + ΔY / 2) atol=.1
129+
@test_throws BoundsError interp[XMIN - ΔX / 2,YMIN - ΔY / 2]
130+
@test_throws BoundsError interp[XMIN - ΔX / 2,YMIN + ΔY / 2]
131+
@test_throws BoundsError interp[XMIN + ΔX / 2,YMIN - ΔY / 2]
132+
@test_throws BoundsError interp[XMAX + ΔX / 2,YMAX + ΔY / 2]
133+
end
134+
135+
@testset "2d-irregular-grids" begin
136+
xs = [x^2 for x in XMIN:ΔX:XMAX]
137+
ys = [y^2 for y in YMIN:ΔY:YMAX]
138+
xmin = xs[1]
139+
xmax = xs[XLEN]
140+
ymin = ys[1]
141+
ymax = ys[YLEN]
142+
f(x, y) = log(x+y)
143+
A = [f(x,y) for x in xs, y in ys]
144+
interp = LinearInterpolation((xs, ys), A)
145+
interp_full = extrapolate(interpolate((xs, ys), A, Gridded(Linear())), Interpolations.Throw())
146+
147+
@test typeof(interp) == typeof(interp_full)
148+
@test interp[xmin,ymin] ≈ f(xmin,ymin)
149+
@test interp[xmin,ymax] ≈ f(xmin,ymax)
150+
@test interp[xmax,ymin] ≈ f(xmax,ymin)
151+
@test interp[xmax,ymax] ≈ f(xmax,ymax)
152+
@test interp[xs[2],ymin] ≈ f(xs[2],ymin)
153+
@test interp[xmin,ys[2]] ≈ f(xmin,ys[2])
154+
@test interp[xs[2],ys[2]] ≈ f(xs[2],ys[2])
155+
@test interp[xmin + ΔX / 2,ymin + ΔY / 2] ≈ f(xmin + ΔX / 2,ymin + ΔY / 2) atol=.1
156+
@test_throws BoundsError interp[xmin - ΔX / 2,ymin - ΔY / 2]
157+
@test_throws BoundsError interp[xmin - ΔX / 2,ymin + ΔY / 2]
158+
@test_throws BoundsError interp[xmin + ΔX / 2,ymin - ΔY / 2]
159+
@test_throws BoundsError interp[xmax + ΔX / 2,ymax + ΔY / 2]
160+
end
161+
162+
@testset "2d-handling-extrapolation" begin
163+
xs = XMIN:ΔX:XMAX
164+
ys = YMIN:ΔY:YMAX
165+
f(x, y) = log(x+y)
166+
A = [f(x,y) for x in xs, y in ys]
167+
ΔA_l = A[2, 1] - A[1, 1]
168+
ΔA_h = A[end, end] - A[end - 1, end]
169+
x_lower = XMIN - ΔX
170+
x_higher = XMAX + ΔX
171+
y_lower = YMIN - ΔY
172+
y_higher = YMAX + ΔY
173+
174+
extrap = LinearInterpolation((xs, ys), A, extrapolation_bc = (Interpolations.Linear(), Interpolations.Flat()))
175+
extrap_full = extrapolate(scale(interpolate(A, BSpline(Linear()), OnGrid()), xs, ys), (Interpolations.Linear(), Interpolations.Flat()))
176+
177+
@test typeof(extrap) == typeof(extrap_full)
178+
@test extrap[x_lower, y_lower] ≈ A[1, 1] - ΔA_l
179+
@test extrap[x_higher, y_higher] ≈ A[end, end] + ΔA_h
180+
end
181+
end
182+
183+
end

‎test/runtests.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@ include("typing.jl")
2828
include("issues/runtests.jl")
2929

3030
include("io.jl")
31+
include("convenience-constructors.jl")
3132
include("readme-examples.jl")
3233

3334
end

0 commit comments

Comments
 (0)