Merge pull request #146 from JuliaMath/teh/indices2

Support arbitrary indices for interpolation & extrapolation

src/Interpolations.jl

@@ -37,7 +37,8 @@ export
 using Compat
 using WoodburyMatrices, Ratios, AxisAlgorithms
 
-import Base: convert, size, getindex, gradient, promote_rule, ndims, eltype, checkbounds
+import Base: convert, size, indices, getindex, gradient, promote_rule,
+             ndims, eltype, checkbounds
 
 # Julia v0.5 compatibility
 if isdefined(:scaling) import Base.scaling end

src/b-splines/b-splines.jl

@@ -38,17 +38,28 @@ padding(itp::AbstractInterpolation) = padding(typeof(itp))
 padextract(pad::Integer, d) = pad
 padextract(pad::Tuple{Vararg{Integer}}, d) = pad[d]
 
-lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) = 1
-ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) = size(itp, d)
-lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) = 0.5
-ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) = size(itp, d)+0.5
+lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) =
+    first(indices(itp, d))
+ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) =
+    last(indices(itp, d))
+lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) =
+    first(indices(itp, d)) - 0.5
+ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) =
+    last(indices(itp, d))+0.5
 
 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)
 
 function size{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, d)
     d <= N ? size(itp.coefs, d) - 2*padextract(pad, d) : 1
 end
 
+@inline indices{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}) =
+    map_repeat(indices_removepad, indices(itp.coefs), pad)
+
+function indices{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, d)
+    d <= N ? indices_removepad(indices(itp.coefs, d), padextract(pad, d)) : indices(itp.coefs, d)
+end
+
 function interpolate{TWeights,TC,IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(::Type{TWeights}, ::Type{TC}, A, it::IT, gt::GT)
     Apad, Pad = prefilter(TWeights, TC, A, IT, GT)
     BSplineInterpolation(TWeights, Apad, it, gt, Pad)
@@ -78,6 +89,36 @@ offsetsym(off, d) = off == -1 ? Symbol("ixm_", d) :
                     off ==  1 ? Symbol("ixp_", d) :
                     off ==  2 ? Symbol("ixpp_", d) : error("offset $off not recognized")
 
+# Ideally we might want to shift the indices symmetrically, but this
+# would introduce an inconsistency, so we just append on the right
+@inline indices_removepad(inds::Base.OneTo, pad) = Base.OneTo(length(inds) - 2*pad)
+@inline indices_removepad(inds, pad) = oftype(inds, first(inds):last(inds) - 2*pad)
+@inline indices_addpad(inds::Base.OneTo, pad) = Base.OneTo(length(inds) + 2*pad)
+@inline indices_addpad(inds, pad) = oftype(inds, first(inds):last(inds) + 2*pad)
+@inline indices_interior(inds, pad) = first(inds)+pad:last(inds)-pad
+
+"""
+    map_repeat(f, a, b)
+
+Equivalent to `(f(a[1], b[1]), f(a[2], b[2]), ...)` if `a` and `b` are
+tuples of the same lengths, or `(f(a[1], b), f(a[2], b), ...)` if `b`
+is a scalar.
+"""
+@generated function map_repeat{N}(f, a::NTuple{N,Any}, b::NTuple{N,Any})
+    ex = [:(f(a[$i], b[$i])) for i = 1:N]
+    quote
+        $(Expr(:meta, :inline))
+        ($(ex...),)
+    end
+end
+@generated function map_repeat{N}(f, a::NTuple{N,Any}, b)
+    ex = [:(f(a[$i], b)) for i = 1:N]
+    quote
+        $(Expr(:meta, :inline))
+        ($(ex...),)
+    end
+end
+
 include("constant.jl")
 include("linear.jl")
 include("quadratic.jl")

src/b-splines/constant.jl

@@ -11,7 +11,7 @@ Constant
 """
 function define_indices_d(::Type{BSpline{Constant}}, d, pad)
     symix, symx = Symbol("ix_",d), Symbol("x_",d)
-    :($symix = clamp(round(Int, $symx), 1, size(itp, $d)))
+    :($symix = clamp(round(Int, $symx), first(inds_itp[$d]), last(inds_itp[$d])))
 end
 
 """
@@ -50,4 +50,4 @@ function index_gen{IT<:DimSpec{BSpline}}(::Type{BSpline{Constant}}, ::Type{IT},
         indices = [offsetsym(offsets[d], d) for d = 1:N]
         return :(itp.coefs[$(indices...)])
     end
-end
+end

src/b-splines/cubic.jl

@@ -35,7 +35,7 @@ function define_indices_d{BC}(::Type{BSpline{Cubic{BC}}}, d, pad)
     quote
         # ensure that all of ix_d, ixm_d, ixp_d, and ixpp_d are in-bounds no
         # matter the value of pad
-        $symix = clamp(floor(Int, $symx), $(2-pad), size(itp,$d)+$(pad-2))
+        $symix = clamp(floor(Int, $symx), first(inds_itp[$d]) + $(1-pad), last(inds_itp[$d]) + $(pad-2))
         $symfx = $symx - $symix
         $symix += $pad # padding for oob coefficient
         $symixm = $symix - 1
@@ -56,11 +56,12 @@ function define_indices_d(::Type{BSpline{Cubic{Periodic}}}, d, pad)
     symix, symixm, symixp = Symbol("ix_",d), Symbol("ixm_",d), Symbol("ixp_",d)
     symixpp, symx, symfx = Symbol("ixpp_",d), Symbol("x_",d), Symbol("fx_",d)
     quote
-        $symix = clamp(floor(Int, $symx), 1, size(itp,$d))
+        tmp = inds_itp[$d]
+        $symix = clamp(floor(Int, $symx), first(tmp), last(tmp))
         $symfx = $symx - $symix
-        $symixm = mod1($symix - 1, size(itp,$d))
-        $symixp = mod1($symix + 1, size(itp,$d))
-        $symixpp = mod1($symix + 2, size(itp,$d))
+        $symixm = modrange($symix - 1, tmp)
+        $symixp = modrange($symix + 1, tmp)
+        $symixpp = modrange($symix + 2, tmp)
     end
 end
 
@@ -316,4 +317,4 @@ function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int,
                                   )
 
     Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
-end
+end

src/b-splines/indexing.jl

@@ -32,6 +32,7 @@ function getindex_impl{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad
     quote
         $meta
         @nexprs $N d->(x_d = xs[d])
+        inds_itp = indices(itp)
 
         # Calculate the indices of all coefficients that will be used
         # 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
         $meta
         length(g) == $n || throw(ArgumentError(string("The length of the provided gradient vector (", length(g), ") did not match the number of interpolating dimensions (", n, ")")))
         @nexprs $N d->(x_d = xs[d])
+        inds_itp = indices(itp)
 
         # Calculate the indices of all coefficients that will be used
         # 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}
         $meta
         size(H) == ($n,$n) || throw(ArgumentError(string("The size of the provided Hessian matrix wasn't a square matrix of size ", size(H))))
         @nexprs $N d->(x_d = xs[d])
+        inds_itp = indices(itp)
 
         $(define_indices(IT, N, Pad))
 

src/b-splines/linear.jl

@@ -29,7 +29,7 @@ Linear
 function define_indices_d(::Type{BSpline{Linear}}, d, pad)
     symix, symixp, symfx, symx = Symbol("ix_",d), Symbol("ixp_",d), Symbol("fx_",d), Symbol("x_",d)
     quote
-        $symix = clamp(floor(Int, $symx), 1, size(itp, $d)-1)
+        $symix = clamp(floor(Int, $symx), first(inds_itp[$d]), last(inds_itp[$d])-1)
         $symixp = $symix + 1
         $symfx = $symx - $symix
     end
@@ -100,4 +100,4 @@ function index_gen{IT<:DimSpec{BSpline}}(::Type{BSpline{Linear}}, ::Type{IT}, N:
         indices = [offsetsym(offsets[d], d) for d = 1:N]
         return :(itp.coefs[$(indices...)])
     end
-end
+end

src/b-splines/prefiltering.jl

@@ -6,31 +6,32 @@ padding{IT<:BSpline}(::Type{IT}) = Val{0}()
     :(Val{$t}())
 end
 
-function padded_index{N,pad}(sz::NTuple{N,Int}, ::Val{pad})
-    szpad = ntuple(i->sz[i]+2padextract(pad,i), N)::NTuple{N,Int}
-    ind = Array{UnitRange{Int}}(N)
-    for i in 1:N
-        p = padextract(pad,i)
-        ind[i] = 1+p:szpad[i]-p
-    end
-    ind,szpad
+@noinline function padded_index{N,pad}(indsA::NTuple{N,AbstractUnitRange{Int}}, ::Val{pad})
+    indspad = ntuple(i->indices_addpad(indsA[i], padextract(pad,i)), Val{N})
+    indscp = ntuple(i->indices_interior(indspad[i], padextract(pad,i)), Val{N})
+    indscp, indspad
 end
 
 copy_with_padding{IT}(A, ::Type{IT}) = copy_with_padding(eltype(A), A, IT)
 function copy_with_padding{TC,IT<:DimSpec{InterpolationType}}(::Type{TC}, A, ::Type{IT})
     Pad = padding(IT)
-    ind,sz = padded_index(size(A), Pad)
-    if sz == size(A)
-        coefs = copy!(Array{TC}(size(A)), A)
+    indsA = indices(A)
+    indscp, indspad = padded_index(indsA, Pad)
+    coefs = similar(dims->Array{TC}(dims), indspad)
+    if indspad == indsA
+        coefs = copy!(coefs, A)
     else
-        coefs = zeros(TC, sz...)
-        coefs[ind...] = A
+        fill!(coefs, zero(TC))
+        copy!(coefs, CartesianRange(indscp), A, CartesianRange(indsA))
     end
     coefs, Pad
 end
 
 prefilter!{TWeights, IT<:BSpline, GT<:GridType}(::Type{TWeights}, A, ::Type{IT}, ::Type{GT}) = A
-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}()
+function prefilter{TWeights, TC, IT<:BSpline, GT<:GridType}(::Type{TWeights}, ::Type{TC}, A, ::Type{IT}, ::Type{GT})
+    coefs = similar(dims->Array{TC}(dims), indices(A))
+    prefilter!(TWeights, copy!(coefs, A), IT, GT), Val{0}()
+end
 
 function prefilter{TWeights,TC,IT<:Union{Cubic,Quadratic},GT<:GridType}(
     ::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
     ::Type{TWeights}, ret::TCoefs, ::Type{BSpline{IT}}, ::Type{GT}
     )
     local buf, shape, retrs
-    sz = size(ret)
+    sz = map(length, indices(ret))
     first = true
     for dim in 1:ndims(ret)
         M, b = prefiltering_system(TWeights, eltype(TCoefs), sz[dim], IT, GT)
@@ -109,4 +110,4 @@ type `IT`.
 
 `dl`, `d`, and `du` are intended to be used e.g. as in `M = Tridiagonal(dl, d, du)`
 """
-inner_system_diags
+inner_system_diags

src/b-splines/quadratic.jl

@@ -34,7 +34,7 @@ function define_indices_d{BC}(::Type{BSpline{Quadratic{BC}}}, d, pad)
     quote
         # ensure that all three ix_d, ixm_d, and ixp_d are in-bounds no matter
         # the value of pad
-        $symix = clamp(round(Int, $symx), 2-$pad, size(itp,$d)+$pad-1)
+        $symix = clamp(round(Int, $symx), first(inds_itp[$d])+1-$pad, last(inds_itp[$d])+$pad-1)
         $symfx = $symx - $symix
         $symix += $pad # padding for oob coefficient
         $symixp = $symix + 1
@@ -45,10 +45,10 @@ function define_indices_d(::Type{BSpline{Quadratic{Periodic}}}, d, pad)
     symix, symixm, symixp = Symbol("ix_",d), Symbol("ixm_",d), Symbol("ixp_",d)
     symx, symfx = Symbol("x_",d), Symbol("fx_",d)
     quote
-        $symix = clamp(round(Int, $symx), 1, size(itp,$d))
+        $symix = clamp(round(Int, $symx), first(inds_itp[$d]), last(inds_itp[$d]))
         $symfx = $symx - $symix
-        $symixp = mod1($symix + 1, size(itp,$d))
-        $symixm = mod1($symix - 1, size(itp,$d))
+        $symixp = modrange($symix + 1, inds_itp[$d])
+        $symixm = modrange($symix - 1, inds_itp[$d])
     end
 end
 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{
     quote
         # ensure that all three ix_d, ixm_d, and ixp_d are in-bounds no matter
         # the value of pad
-        $symix = clamp(round(Int, $symx), 1, size(itp,$d))
+        $symix = clamp(round(Int, $symx), first(inds_itp[$d]), last(inds_itp[$d]))
         $symfx = $symx - $symix
         $symix += $pad # padding for oob coefficient
-        $symixp = min(size(itp,$d), $symix + 1)
-        $symixm = max(1, $symix - 1)
+        $symixp = min(last(inds_itp[$d]), $symix + 1)
+        $symixm = max(first(inds_itp[$d]), $symix - 1)
     end
 end
 
@@ -263,4 +263,4 @@ function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int, :
                                   )
 
     Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
-end
+end

src/extrapolation/extrapolation.jl

@@ -78,5 +78,7 @@ end
 lbound(etp::Extrapolation, d) = lbound(etp.itp, d)
 ubound(etp::Extrapolation, d) = ubound(etp.itp, d)
 size(etp::Extrapolation, d) = size(etp.itp, d)
+indices(etp::AbstractExtrapolation) = indices(etp.itp)
+indices(etp::AbstractExtrapolation, d) = indices(etp.itp, d)
 
 include("filled.jl")

src/scaling/scaling.jl

@@ -241,6 +241,7 @@ function next_gen{CR,SITPT,X1,Deg,T}(::Type{ScaledIterator{CR,SITPT,X1,Deg,T}})
     quote
         sitp = iter.sitp
         itp = sitp.itp
+        inds_itp = indices(itp)
         if iter.nremaining > 0
             iter.nremaining -= 1
             iter.fx_1 += iter.dx_1

src/utils.jl

@@ -1,2 +1,4 @@
 @inline sqr(x) = x*x
-@inline cub(x) = x*x*x
+@inline cub(x) = x*x*x
+
+modrange(x, r::AbstractUnitRange) = mod(x-first(r), length(r)) + first(r)

test/REQUIRE

@@ -0,0 +1 @@
+OffsetArrays

test/b-splines/non1.jl

@@ -0,0 +1,93 @@
+module Non1Tests
+
+using Interpolations, OffsetArrays, AxisAlgorithms, Base.Test
+
+# At present, for a particular type of non-1 array you need to specialize this function
+function AxisAlgorithms.A_ldiv_B_md!(dest::OffsetArray, F, src::OffsetArray, dim::Integer, b::AbstractVector)
+    indsdim = indices(parent(src), dim)
+    indsF = indices(F)[2]
+    if indsF == indsdim
+        return A_ldiv_B_md!(parent(dest), F, parent(src), dim, b)
+    end
+    throw(DimensionMismatch("indices $(indices(parent(src))) do not match $(indices(F))"))
+end
+
+for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy))
+    f1(x) = sin((x-3)*2pi/9 - 1)
+    inds = -3:6
+    A1 = OffsetArray(Float64[f1(x) for x in inds], inds)
+
+    f2(x,y) = sin(x/10)*cos(y/6) + 0.1
+    xinds, yinds = -2:28,0:9
+    A2 = OffsetArray(Float64[f2(x,y) for x in xinds, y in yinds], xinds, yinds)
+
+    for GT in (OnGrid, OnCell), O in (Constant, Linear)
+        itp1 = @inferred(constructor(copier(A1), BSpline(O()), GT()))
+        @test @inferred(indices(itp1)) === indices(A1)
+
+        # test that we reproduce the values at on-grid points
+        for x = inds
+            @test itp1[x] ≈ f1(x)
+        end
+
+        itp2 = @inferred(constructor(copier(A2), BSpline(O()), GT()))
+        @test @inferred(indices(itp2)) === indices(A2)
+        for j = yinds, i = xinds
+            @test itp2[i,j] ≈ A2[i,j]
+        end
+    end
+
+    for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell)
+        itp1 = @inferred(constructor(copier(A1), BSpline(Quadratic(BC())), GT()))
+        @test @inferred(indices(itp1)) === indices(A1)
+
+        # test that we reproduce the values at on-grid points
+        inset = constructor == interpolate!
+        for x = first(inds)+inset:last(inds)-inset
+            @test itp1[x] ≈ f1(x)
+        end
+
+        itp2 = @inferred(constructor(copier(A2), BSpline(Quadratic(BC())), GT()))
+        @test @inferred(indices(itp2)) === indices(A2)
+        for j = first(yinds)+inset:last(yinds)-inset, i = first(xinds)+inset:last(xinds)-inset
+            @test itp2[i,j] ≈ A2[i,j]
+        end
+    end
+
+    for BC in (Flat,Line,Free,Periodic), GT in (OnGrid, OnCell)
+        itp1 = @inferred(constructor(copier(A1), BSpline(Cubic(BC())), GT()))
+        @test @inferred(indices(itp1)) === indices(A1)
+
+        # test that we reproduce the values at on-grid points
+        inset = constructor == interpolate!
+        for x = first(inds)+inset:last(inds)-inset
+            @test itp1[x] ≈ f1(x)
+        end
+
+        itp2 = @inferred(constructor(copier(A2), BSpline(Cubic(BC())), GT()))
+        @test @inferred(indices(itp2)) === indices(A2)
+        for j = first(yinds)+inset:last(yinds)-inset, i = first(xinds)+inset:last(xinds)-inset
+            @test itp2[i,j] ≈ A2[i,j]
+        end
+    end
+end
+
+let
+    f(x) = sin((x-3)*2pi/9 - 1)
+    inds = -7:2
+    A = OffsetArray(Float64[f(x) for x in inds], inds)
+    itp1 = interpolate!(copy(A), BSpline(Quadratic(InPlace())), OnCell())
+    for i in inds
+        @test itp1[i] ≈ A[i]
+    end
+
+    f(x,y) = sin(x/10)*cos(y/6) + 0.1
+    xinds, yinds = -2:28,0:9
+    A2 = OffsetArray(Float64[f(x,y) for x in xinds, y in yinds], xinds, yinds)
+    itp2 = interpolate!(copy(A2), BSpline(Quadratic(InPlace())), OnCell())
+    for j = yinds, i = xinds
+        @test itp2[i,j] ≈ A2[i,j]
+    end
+end
+
+end