Skip to content

Commit 0ecd53b

Browse files
timholytimholy
authored
Merge pull request #146 from JuliaMath/teh/indices2
Support arbitrary indices for interpolation & extrapolation
2 parents 380b80c + 2195777 commit 0ecd53b

File tree

17 files changed

+223
-42
lines changed

17 files changed

+223
-42
lines changed

‎src/Interpolations.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,8 @@ export
3737
using Compat
3838
using WoodburyMatrices, Ratios, AxisAlgorithms
3939

40-
import Base: convert, size, getindex, gradient, promote_rule, ndims, eltype, checkbounds
40+
import Base: convert, size, indices, getindex, gradient, promote_rule,
41+
ndims, eltype, checkbounds
4142

4243
# Julia v0.5 compatibility
4344
if isdefined(:scaling) import Base.scaling end

‎src/b-splines/b-splines.jl

Lines changed: 45 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -38,17 +38,28 @@ padding(itp::AbstractInterpolation) = padding(typeof(itp))
3838
padextract(pad::Integer, d) = pad
3939
padextract(pad::Tuple{Vararg{Integer}}, d) = pad[d]
4040

41-
lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) = 1
42-
ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) = size(itp, d)
43-
lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) = 0.5
44-
ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) = size(itp, d)+0.5
41+
lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) =
42+
first(indices(itp, d))
43+
ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) =
44+
last(indices(itp, d))
45+
lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) =
46+
first(indices(itp, d)) - 0.5
47+
ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) =
48+
last(indices(itp, d))+0.5
4549

4650
count_interp_dims{T,N,TCoefs,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType},pad}(::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,pad}}, n) = count_interp_dims(IT, n)
4751

4852
function size{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, d)
4953
d <= N ? size(itp.coefs, d) - 2*padextract(pad, d) : 1
5054
end
5155

56+
@inline indices{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}) =
57+
map_repeat(indices_removepad, indices(itp.coefs), pad)
58+
59+
function indices{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, d)
60+
d <= N ? indices_removepad(indices(itp.coefs, d), padextract(pad, d)) : indices(itp.coefs, d)
61+
end
62+
5263
function interpolate{TWeights,TC,IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(::Type{TWeights}, ::Type{TC}, A, it::IT, gt::GT)
5364
Apad, Pad = prefilter(TWeights, TC, A, IT, GT)
5465
BSplineInterpolation(TWeights, Apad, it, gt, Pad)
@@ -78,6 +89,36 @@ offsetsym(off, d) = off == -1 ? Symbol("ixm_", d) :
7889
off == 1 ? Symbol("ixp_", d) :
7990
off == 2 ? Symbol("ixpp_", d) : error("offset $off not recognized")
8091

92+
# Ideally we might want to shift the indices symmetrically, but this
93+
# would introduce an inconsistency, so we just append on the right
94+
@inline indices_removepad(inds::Base.OneTo, pad) = Base.OneTo(length(inds) - 2*pad)
95+
@inline indices_removepad(inds, pad) = oftype(inds, first(inds):last(inds) - 2*pad)
96+
@inline indices_addpad(inds::Base.OneTo, pad) = Base.OneTo(length(inds) + 2*pad)
97+
@inline indices_addpad(inds, pad) = oftype(inds, first(inds):last(inds) + 2*pad)
98+
@inline indices_interior(inds, pad) = first(inds)+pad:last(inds)-pad
99+
100+
"""
101+
map_repeat(f, a, b)
102+
103+
Equivalent to `(f(a[1], b[1]), f(a[2], b[2]), ...)` if `a` and `b` are
104+
tuples of the same lengths, or `(f(a[1], b), f(a[2], b), ...)` if `b`
105+
is a scalar.
106+
"""
107+
@generated function map_repeat{N}(f, a::NTuple{N,Any}, b::NTuple{N,Any})
108+
ex = [:(f(a[$i], b[$i])) for i = 1:N]
109+
quote
110+
$(Expr(:meta, :inline))
111+
($(ex...),)
112+
end
113+
end
114+
@generated function map_repeat{N}(f, a::NTuple{N,Any}, b)
115+
ex = [:(f(a[$i], b)) for i = 1:N]
116+
quote
117+
$(Expr(:meta, :inline))
118+
($(ex...),)
119+
end
120+
end
121+
81122
include("constant.jl")
82123
include("linear.jl")
83124
include("quadratic.jl")

‎src/b-splines/constant.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ Constant
1111
"""
1212
function define_indices_d(::Type{BSpline{Constant}}, d, pad)
1313
symix, symx = Symbol("ix_",d), Symbol("x_",d)
14-
:($symix = clamp(round(Int, $symx), 1, size(itp, $d)))
14+
:($symix = clamp(round(Int, $symx), first(inds_itp[$d]), last(inds_itp[$d])))
1515
end
1616

1717
"""
@@ -50,4 +50,4 @@ function index_gen{IT<:DimSpec{BSpline}}(::Type{BSpline{Constant}}, ::Type{IT},
5050
indices = [offsetsym(offsets[d], d) for d = 1:N]
5151
return :(itp.coefs[$(indices...)])
5252
end
53-
end
53+
end

‎src/b-splines/cubic.jl

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@ function define_indices_d{BC}(::Type{BSpline{Cubic{BC}}}, d, pad)
3535
quote
3636
# ensure that all of ix_d, ixm_d, ixp_d, and ixpp_d are in-bounds no
3737
# matter the value of pad
38-
$symix = clamp(floor(Int, $symx), $(2-pad), size(itp,$d)+$(pad-2))
38+
$symix = clamp(floor(Int, $symx), first(inds_itp[$d]) + $(1-pad), last(inds_itp[$d]) + $(pad-2))
3939
$symfx = $symx - $symix
4040
$symix += $pad # padding for oob coefficient
4141
$symixm = $symix - 1
@@ -56,11 +56,12 @@ function define_indices_d(::Type{BSpline{Cubic{Periodic}}}, d, pad)
5656
symix, symixm, symixp = Symbol("ix_",d), Symbol("ixm_",d), Symbol("ixp_",d)
5757
symixpp, symx, symfx = Symbol("ixpp_",d), Symbol("x_",d), Symbol("fx_",d)
5858
quote
59-
$symix = clamp(floor(Int, $symx), 1, size(itp,$d))
59+
tmp = inds_itp[$d]
60+
$symix = clamp(floor(Int, $symx), first(tmp), last(tmp))
6061
$symfx = $symx - $symix
61-
$symixm = mod1($symix - 1, size(itp,$d))
62-
$symixp = mod1($symix + 1, size(itp,$d))
63-
$symixpp = mod1($symix + 2, size(itp,$d))
62+
$symixm = modrange($symix - 1, tmp)
63+
$symixp = modrange($symix + 1, tmp)
64+
$symixpp = modrange($symix + 2, tmp)
6465
end
6566
end
6667

@@ -316,4 +317,4 @@ function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int,
316317
)
317318

318319
Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
319-
end
320+
end

‎src/b-splines/indexing.jl

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,7 @@ function getindex_impl{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad
3232
quote
3333
$meta
3434
@nexprs $N d->(x_d = xs[d])
35+
inds_itp = indices(itp)
3536

3637
# Calculate the indices of all coefficients that will be used
3738
# and define fx = x - xi in each dimension
@@ -69,6 +70,7 @@ function gradient_impl{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad
6970
$meta
7071
length(g) == $n || throw(ArgumentError(string("The length of the provided gradient vector (", length(g), ") did not match the number of interpolating dimensions (", n, ")")))
7172
@nexprs $N d->(x_d = xs[d])
73+
inds_itp = indices(itp)
7274

7375
# Calculate the indices of all coefficients that will be used
7476
# and define fx = x - xi in each dimension
@@ -130,6 +132,7 @@ function hessian_impl{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad}
130132
$meta
131133
size(H) == ($n,$n) || throw(ArgumentError(string("The size of the provided Hessian matrix wasn't a square matrix of size ", size(H))))
132134
@nexprs $N d->(x_d = xs[d])
135+
inds_itp = indices(itp)
133136

134137
$(define_indices(IT, N, Pad))
135138

‎src/b-splines/linear.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ Linear
2929
function define_indices_d(::Type{BSpline{Linear}}, d, pad)
3030
symix, symixp, symfx, symx = Symbol("ix_",d), Symbol("ixp_",d), Symbol("fx_",d), Symbol("x_",d)
3131
quote
32-
$symix = clamp(floor(Int, $symx), 1, size(itp, $d)-1)
32+
$symix = clamp(floor(Int, $symx), first(inds_itp[$d]), last(inds_itp[$d])-1)
3333
$symixp = $symix + 1
3434
$symfx = $symx - $symix
3535
end
@@ -100,4 +100,4 @@ function index_gen{IT<:DimSpec{BSpline}}(::Type{BSpline{Linear}}, ::Type{IT}, N:
100100
indices = [offsetsym(offsets[d], d) for d = 1:N]
101101
return :(itp.coefs[$(indices...)])
102102
end
103-
end
103+
end

‎src/b-splines/prefiltering.jl

Lines changed: 17 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -6,31 +6,32 @@ padding{IT<:BSpline}(::Type{IT}) = Val{0}()
66
:(Val{$t}())
77
end
88

9-
function padded_index{N,pad}(sz::NTuple{N,Int}, ::Val{pad})
10-
szpad = ntuple(i->sz[i]+2padextract(pad,i), N)::NTuple{N,Int}
11-
ind = Array{UnitRange{Int}}(N)
12-
for i in 1:N
13-
p = padextract(pad,i)
14-
ind[i] = 1+p:szpad[i]-p
15-
end
16-
ind,szpad
9+
@noinline function padded_index{N,pad}(indsA::NTuple{N,AbstractUnitRange{Int}}, ::Val{pad})
10+
indspad = ntuple(i->indices_addpad(indsA[i], padextract(pad,i)), Val{N})
11+
indscp = ntuple(i->indices_interior(indspad[i], padextract(pad,i)), Val{N})
12+
indscp, indspad
1713
end
1814

1915
copy_with_padding{IT}(A, ::Type{IT}) = copy_with_padding(eltype(A), A, IT)
2016
function copy_with_padding{TC,IT<:DimSpec{InterpolationType}}(::Type{TC}, A, ::Type{IT})
2117
Pad = padding(IT)
22-
ind,sz = padded_index(size(A), Pad)
23-
if sz == size(A)
24-
coefs = copy!(Array{TC}(size(A)), A)
18+
indsA = indices(A)
19+
indscp, indspad = padded_index(indsA, Pad)
20+
coefs = similar(dims->Array{TC}(dims), indspad)
21+
if indspad == indsA
22+
coefs = copy!(coefs, A)
2523
else
26-
coefs = zeros(TC, sz...)
27-
coefs[ind...] = A
24+
fill!(coefs, zero(TC))
25+
copy!(coefs, CartesianRange(indscp), A, CartesianRange(indsA))
2826
end
2927
coefs, Pad
3028
end
3129

3230
prefilter!{TWeights, IT<:BSpline, GT<:GridType}(::Type{TWeights}, A, ::Type{IT}, ::Type{GT}) = A
33-
prefilter{TWeights, TC, IT<:BSpline, GT<:GridType}(::Type{TWeights}, ::Type{TC}, A, ::Type{IT}, ::Type{GT}) = prefilter!(TWeights, copy!(Array{TC}(size(A)), A), IT, GT), Val{0}()
31+
function prefilter{TWeights, TC, IT<:BSpline, GT<:GridType}(::Type{TWeights}, ::Type{TC}, A, ::Type{IT}, ::Type{GT})
32+
coefs = similar(dims->Array{TC}(dims), indices(A))
33+
prefilter!(TWeights, copy!(coefs, A), IT, GT), Val{0}()
34+
end
3435

3536
function prefilter{TWeights,TC,IT<:Union{Cubic,Quadratic},GT<:GridType}(
3637
::Type{TWeights}, ::Type{TC}, A::AbstractArray, ::Type{BSpline{IT}}, ::Type{GT}
@@ -50,7 +51,7 @@ function prefilter!{TWeights,TCoefs<:AbstractArray,IT<:Union{Quadratic,Cubic},GT
5051
::Type{TWeights}, ret::TCoefs, ::Type{BSpline{IT}}, ::Type{GT}
5152
)
5253
local buf, shape, retrs
53-
sz = size(ret)
54+
sz = map(length, indices(ret))
5455
first = true
5556
for dim in 1:ndims(ret)
5657
M, b = prefiltering_system(TWeights, eltype(TCoefs), sz[dim], IT, GT)
@@ -109,4 +110,4 @@ type `IT`.
109110
110111
`dl`, `d`, and `du` are intended to be used e.g. as in `M = Tridiagonal(dl, d, du)`
111112
"""
112-
inner_system_diags
113+
inner_system_diags

‎src/b-splines/quadratic.jl

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@ function define_indices_d{BC}(::Type{BSpline{Quadratic{BC}}}, d, pad)
3434
quote
3535
# ensure that all three ix_d, ixm_d, and ixp_d are in-bounds no matter
3636
# the value of pad
37-
$symix = clamp(round(Int, $symx), 2-$pad, size(itp,$d)+$pad-1)
37+
$symix = clamp(round(Int, $symx), first(inds_itp[$d])+1-$pad, last(inds_itp[$d])+$pad-1)
3838
$symfx = $symx - $symix
3939
$symix += $pad # padding for oob coefficient
4040
$symixp = $symix + 1
@@ -45,10 +45,10 @@ function define_indices_d(::Type{BSpline{Quadratic{Periodic}}}, d, pad)
4545
symix, symixm, symixp = Symbol("ix_",d), Symbol("ixm_",d), Symbol("ixp_",d)
4646
symx, symfx = Symbol("x_",d), Symbol("fx_",d)
4747
quote
48-
$symix = clamp(round(Int, $symx), 1, size(itp,$d))
48+
$symix = clamp(round(Int, $symx), first(inds_itp[$d]), last(inds_itp[$d]))
4949
$symfx = $symx - $symix
50-
$symixp = mod1($symix + 1, size(itp,$d))
51-
$symixm = mod1($symix - 1, size(itp,$d))
50+
$symixp = modrange($symix + 1, inds_itp[$d])
51+
$symixm = modrange($symix - 1, inds_itp[$d])
5252
end
5353
end
5454
function define_indices_d{BC<:Union{InPlace,InPlaceQ}}(::Type{BSpline{Quadratic{BC}}}, d, pad)
@@ -58,11 +58,11 @@ function define_indices_d{BC<:Union{InPlace,InPlaceQ}}(::Type{BSpline{Quadratic{
5858
quote
5959
# ensure that all three ix_d, ixm_d, and ixp_d are in-bounds no matter
6060
# the value of pad
61-
$symix = clamp(round(Int, $symx), 1, size(itp,$d))
61+
$symix = clamp(round(Int, $symx), first(inds_itp[$d]), last(inds_itp[$d]))
6262
$symfx = $symx - $symix
6363
$symix += $pad # padding for oob coefficient
64-
$symixp = min(size(itp,$d), $symix + 1)
65-
$symixm = max(1, $symix - 1)
64+
$symixp = min(last(inds_itp[$d]), $symix + 1)
65+
$symixm = max(first(inds_itp[$d]), $symix - 1)
6666
end
6767
end
6868

@@ -263,4 +263,4 @@ function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int, :
263263
)
264264

265265
Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
266-
end
266+
end

‎src/extrapolation/extrapolation.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -78,5 +78,7 @@ end
7878
lbound(etp::Extrapolation, d) = lbound(etp.itp, d)
7979
ubound(etp::Extrapolation, d) = ubound(etp.itp, d)
8080
size(etp::Extrapolation, d) = size(etp.itp, d)
81+
indices(etp::AbstractExtrapolation) = indices(etp.itp)
82+
indices(etp::AbstractExtrapolation, d) = indices(etp.itp, d)
8183

8284
include("filled.jl")

‎src/scaling/scaling.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -241,6 +241,7 @@ function next_gen{CR,SITPT,X1,Deg,T}(::Type{ScaledIterator{CR,SITPT,X1,Deg,T}})
241241
quote
242242
sitp = iter.sitp
243243
itp = sitp.itp
244+
inds_itp = indices(itp)
244245
if iter.nremaining > 0
245246
iter.nremaining -= 1
246247
iter.fx_1 += iter.dx_1

0 commit comments

Comments
 (0)