Implement scaling of interpolation objects

This is a start at implementing the functionality provided to Grid.jl
by the CoordInterpGrid type, and in fact the code here is heavily based
on that contributed by @simonbyrne there.

Author: Tomas Lycken

src/Interpolations.jl

@@ -4,6 +4,7 @@ export
     interpolate,
     interpolate!,
     extrapolate,
+    scale,
 
     gradient!,
 
@@ -22,10 +23,11 @@ export
     # see the following files for further exports:
     # b-splines/b-splines.jl
     # extrapolation/extrapolation.jl
+    # scaling/scaling.jl
 
 using WoodburyMatrices, Ratios, AxisAlgorithms
 
-import Base: convert, size, getindex, gradient, promote_rule
+import Base: convert, size, getindex, gradient, scale, promote_rule
 
 abstract InterpolationType
 immutable NoInterp <: InterpolationType end
@@ -36,7 +38,8 @@ immutable OnCell <: GridType end
 typealias DimSpec{T} Union{T,Tuple{Vararg{Union{T,NoInterp}}},NoInterp}
 
 abstract AbstractInterpolation{T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}} <: AbstractArray{T,N}
-abstract AbstractExtrapolation{T,N,ITPT,IT,GT} <: AbstractInterpolation{T,N,IT,GT}
+abstract AbstractInterpolationWrapper{T,N,ITPT,IT,GT} <: AbstractInterpolation{T,N,IT,GT}
+abstract AbstractExtrapolation{T,N,ITPT,IT,GT} <: AbstractInterpolationWrapper{T,N,ITPT,IT,GT}
 
 abstract BoundaryCondition
 immutable Flat <: BoundaryCondition end
@@ -53,7 +56,13 @@ typealias Natural Line
 # TODO: size might have to be faster?
 size{T,N}(itp::AbstractInterpolation{T,N}) = ntuple(i->size(itp,i), N)::NTuple{N,Int}
 size(exp::AbstractExtrapolation, d) = size(exp.itp, d)
- itptype{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}) = IT
+bounds{T,N}(itp::AbstractInterpolation{T,N}) = tuple(zip(lbounds(itp), ubounds(itp))...)
+bounds{T,N}(itp::AbstractInterpolation{T,N}, d) = (lbound(itp,d),ubound(itp,d))
+lbounds{T,N}(itp::AbstractInterpolation{T,N}) = ntuple(i->lbound(itp,i), N)::NTuple{N,T}
+ubounds{T,N}(itp::AbstractInterpolation{T,N}) = ntuple(i->ubound(itp,i), N)::NTuple{N,T}
+lbound{T,N}(itp::AbstractInterpolation{T,N}, d) = convert(T, 1)
+ubound{T,N}(itp::AbstractInterpolation{T,N}, d) = convert(T, size(itp, d))
+itptype{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}) = IT
 gridtype{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}) = GT
 
 @inline gradient{T,N}(itp::AbstractInterpolation{T,N}, xs...) = gradient!(Array(T,N), itp, xs...)
@@ -62,5 +71,6 @@ include("nointerp/nointerp.jl")
 include("b-splines/b-splines.jl")
 include("gridded/gridded.jl")
 include("extrapolation/extrapolation.jl")
+include("scaling/scaling.jl")
 
 end # module

src/b-splines/b-splines.jl

@@ -36,6 +36,11 @@ iextract(t, d) = t.parameters[d]
 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) = one(T)
+ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) = convert(T, size(itp, d))
+lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) = convert(T, .5)
+ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) = convert(T, size(itp, d) + .5)
+
 @generated function size{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, d)
     quote
         d <= $N ? size(itp.coefs, d) - 2*padextract($pad, d) : 1

src/extrapolation/constant.jl

@@ -7,15 +7,12 @@ ConstantExtrapolation{T,ITP,IT,GT}(::Type{T}, N, itp::ITP, ::Type{IT}, ::Type{GT
 extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, ::Type{Flat}) =
     ConstantExtrapolation(T,N,itp,IT,GT)
 
-function extrap_prep{T,ITP,IT}(exp::Type{ConstantExtrapolation{T,1,ITP,IT,OnGrid}}, x)
-    :(x = clamp(x, 1, size(exp,1)))
+function extrap_prep{T,ITP,IT,GT}(etp::Type{ConstantExtrapolation{T,1,ITP,IT,GT}}, x)
+    :(x = clamp(x, lbound(etp,1), ubound(etp,1)))
 end
-function extrap_prep{T,ITP,IT}(exp::Type{ConstantExtrapolation{T,1,ITP,IT,OnCell}}, x)
-    :(x = clamp(x, .5, size(exp,1)+.5))
-end
-function extrap_prep{T,N,ITP,IT}(exp::Type{ConstantExtrapolation{T,N,ITP,IT,OnGrid}}, xs...)
-    :(@nexprs $N d->(xs[d] = clamp(xs[d], 1, size(exp,d))))
-end
-function extrap_prep{T,N,ITP,IT}(exp::Type{ConstantExtrapolation{T,N,ITP,IT,OnCell}}, xs...)
-    :(@nexprs $N d->(xs[d] = clamp(xs[d], .5, size(exp,d)+.5)))
+function extrap_prep{T,N,ITP,IT,GT}(etp::Type{ConstantExtrapolation{T,N,ITP,IT,GT}}, xs...)
+    :(@nexprs $N d->(xs[d] = clamp(xs[d], lbound(etp,d), ubound(etp,d))))
 end
+
+lbound(etp::ConstantExtrapolation, d) = lbound(etp.itp, d)
+ubound(etp::ConstantExtrapolation, d) = ubound(etp.itp, d)

src/extrapolation/error.jl

@@ -7,7 +7,9 @@ ErrorExtrapolation{T,ITPT,IT,GT}(::Type{T}, N, itp::ITPT, ::Type{IT}, ::Type{GT}
 extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, ::Type{Throw}) = 
     ErrorExtrapolation(T,N,itp,IT,GT)
 
-
-function extrap_prep{T,N,ITPT,IT}(exp::Type{ErrorExtrapolation{T,N,ITPT,IT,OnGrid}}, xs...)
-    :(@nexprs $N d->(@show 1 <= xs[d] <= size(exp,d) || throw(BoundsError())))
+function extrap_prep{T,N,ITPT,IT,GT}(etp::Type{ErrorExtrapolation{T,N,ITPT,IT,GT}}, xs...)
+    :(@nexprs $N d->(lbound(etp,d) <= xs[d] <= ubound(etp,d) || throw(BoundsError())))
 end
+
+lbound(etp::ErrorExtrapolation, d) = lbound(etp.itp, d)
+ubound(etp::ErrorExtrapolation, d) = ubound(etp.itp, d)

src/extrapolation/filled.jl

@@ -27,10 +27,13 @@ extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue) = Fille
         $meta
         # Check to see if we're in the extrapolation region, i.e.,
         # out-of-bounds in an index
-        @nexprs $N d->((args[d] < 1 || args[d] > size(fitp.itp, d)) && return fitp.fillvalue)
+        @nexprs $N d->((args[d] < lbound(fitp,d) || args[d] > ubound(fitp, d)) && return fitp.fillvalue)
         # In the interpolation region
         return getindex(fitp.itp,args...)
     end
 end
 
 getindex{T}(fitp::FilledExtrapolation{T,1}, x::Number, y::Int) = y == 1 ? fitp[x] : throw(BoundsError())
+
+lbound(etp::FilledExtrapolation, d) = lbound(etp.itp, d)
+ubound(etp::FilledExtrapolation, d) = ubound(etp.itp, d)

src/extrapolation/indexing.jl

@@ -1,13 +1,13 @@
-@generated function getindex{T}(exp::AbstractExtrapolation{T,1}, x)
+@generated function getindex{T}(etp::AbstractExtrapolation{T,1}, x)
     quote
-        $(extrap_prep(exp, x))
-        exp.itp[x]
+        $(extrap_prep(etp, x))
+        etp.itp[x]
     end
 end
 
-@generated function getindex{T,N,ITP,GT}(exp::AbstractExtrapolation{T,N,ITP,GT}, xs...)
+@generated function getindex{T,N,ITP,GT}(etp::AbstractExtrapolation{T,N,ITP,GT}, xs...)
     quote
-        $(extrap_prep(exp, xs...))
-        exp.itp[xs...]
+        $(extrap_prep(etp, xs...))
+        etp.itp[xs...]
     end
 end

src/scaling/scaling.jl

@@ -0,0 +1,68 @@
+type ScaledInterpolation{T,N,ITPT,IT,GT,RT} <: AbstractInterpolationWrapper{T,N,ITPT,IT,GT}
+    itp::ITPT
+    ranges::RT
+end
+ScaledInterpolation{T,ITPT,IT,GT,RT}(::Type{T}, N, itp::ITPT, ::Type{IT}, ::Type{GT}, ranges::RT) =
+    ScaledInterpolation{T,N,ITPT,IT,GT,RT}(itp, ranges)
+function scale{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, ranges::Range...)
+    length(ranges) == N || throw(ArgumentError("Must scale $N-dimensional interpolation object with exactly $N ranges (you used $(length(ranges)))"))
+    for d in 1:N
+        if iextract(IT,d) != NoInterp 
+            length(ranges[d]) == size(itp,d) || throw(ArgumentError("The length of the range in dimension $d ($(length(ranges[d]))) did not equal the size of the interpolation object in that direction ($(size(itp,d)))"))
+        elseif ranges[d] != 1:size(itp,d)
+            throw(ArgumentError("NoInterp dimension $d must be scaled with unit range 1:$(size(itp,d))"))
+        end
+    end
+
+    ScaledInterpolation(T,N,itp,IT,GT,ranges)
+end
+
+@generated function getindex{T,N,ITPT,IT<:DimSpec}(sitp::ScaledInterpolation{T,N,ITPT,IT}, xs...)
+    length(xs) == N || throw(ArgumentError("Must index into $N-dimensional scaled interpolation object with exactly $N indices (you used $(length(xs)))"))
+    interp_types = length(IT.parameters) == N ? IT.parameters : tuple([IT.parameters[1] for _ in 1:N]...)
+    interp_dimens = map(it -> interp_types[it] != NoInterp, 1:N)
+    interp_indices = map(i -> interp_dimens[i] ? :(coordlookup(sitp.ranges[$i], xs[$i])) : :(xs[$i]), 1:N)
+    return :(getindex(sitp.itp, $(interp_indices...)))
+end
+
+getindex{T}(sitp::ScaledInterpolation{T,1}, x::Number, y::Int) = y == 1 ? sitp[x] : throw(BoundsError())
+
+size(sitp::ScaledInterpolation, d) = size(sitp.itp, d)
+lbound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnGrid}, d) = 1 <= d <= N ? sitp.ranges[d][1] : throw(BoundsError())
+lbound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnCell}, d) = 1 <= d <= N ? sitp.ranges[d][1] - boundstep(sitp.ranges[d]) : throw(BoundsError())
+ubound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnGrid}, d) = 1 <= d <= N ? sitp.ranges[d][end] : throw(BoundsError())
+ubound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnCell}, d) = 1 <= d <= N ? sitp.ranges[d][end] + boundstep(sitp.ranges[d]) : throw(BoundsError())
+
+boundstep(r::LinSpace) = ((r.stop - r.start) / r.divisor) / 2
+boundstep(r::FloatRange) = r.step / 2
+boundstep(r::StepRange) = r.step / 2
+boundstep(r::UnitRange) = 1//2
+
+coordlookup(r::LinSpace, x) = (r.divisor * x + r.stop - r.len * r.start) / (r.stop - r.start)
+coordlookup(r::FloatRange, x) = (r.divisor * x - r.start) / r.step + one(eltype(r))
+coordlookup(r::StepRange, x) = (x - r.start) / r.step + one(eltype(r))
+coordlookup(r::UnitRange, x) = x - r.start + one(eltype(r))
+coordlookup(i::Bool, r::Range, x) = i ? coordlookup(r, x) : convert(typeof(coordlookup(r,x)), x)
+
+gradient{T,N}(sitp::ScaledInterpolation{T,N}, xs...) = gradient!(Array(T,N), sitp, xs...)
+@generated function gradient!{T,N,ITPT,IT}(g, sitp::ScaledInterpolation{T,N,ITPT,IT}, xs...)
+    ndims(g) == 1 || throw(ArgumentError("g must be a vector (but had $(ndims(g)) dimensions)"))
+    length(xs) == count_interp_dims(IT, N) || throw(ArgumentError("Must index into $N-dimensional scaled interpolation object with exactly $N indices (you used $(length(xs)))"))
+
+    interp_types = length(IT.parameters) == N ? IT.parameters : tuple([IT.parameters[1] for _ in 1:N]...)
+    interp_dimens = map(it -> interp_types[it] != NoInterp, 1:N)
+    interp_indices = map(i -> interp_dimens[i] ? :(coordlookup(sitp.ranges[$i], xs[$i])) : :(xs[$i]), 1:N)
+
+    quote
+        length(g) == N || throw(ArgumentError(string("g must be a vector of length ", N, " (was ", length(g), ")")))
+        gradient!(g, sitp.itp, $(interp_indices...))
+        for i in eachindex(g)
+            g[i] = rescale_gradient(sitp.ranges[i], g[i])
+        end
+        g
+    end
+end
+
+rescale_gradient(r::FloatRange, g) = g * r.divisor / r.step
+rescale_gradient(r::StepRange, g) = g / r.step
+rescale_gradient(r::UnitRange, g) = g

test/runtests.jl

@@ -9,6 +9,9 @@ include("b-splines/runtests.jl")
 # extrapolation tests
 include("extrapolation/runtests.jl")
 
+# scaling tests
+include("scaling/runtests.jl")
+
 # # test gradient evaluation
 include("gradient.jl")
 

test/scaling/dimspecs.jl

@@ -0,0 +1,26 @@
+module ScalingDimspecTests
+
+using Interpolations, DualNumbers, Base.Test
+
+xs = -pi:(2pi/10):pi-2pi/10
+ys = -2:.1:2
+f(x,y) = sin(x) * y^2
+
+itp = interpolate(Float64[f(x,y) for x in xs, y in ys], Tuple{BSpline(Quadratic(Periodic)), BSpline(Linear)}, OnGrid)
+sitp = scale(itp, xs, ys)
+
+# Don't test too near the edges until issue #64 is resolved
+for (ix,x0) in enumerate(xs[5:end-5]), (iy,y0) in enumerate(ys[2:end-1])
+    x, y = x0 + 2pi/20, y0 + .05
+    @test_approx_eq sitp[x0, y0] f(x0,y0)
+    @test_approx_eq_eps sitp[x0, y0] f(x0,y0) 0.05
+
+    g = gradient(sitp, x, y)
+    fx = epsilon(f(dual(x,1), y))
+    fy = (f(x, ys[iy+2]) - f(x, ys[iy+1])) / (ys[iy+2] - ys[iy+1])
+
+    @test_approx_eq_eps g[1] fx 0.15
+    @test_approx_eq_eps g[2] fy 0.05 # gradients for linear interpolation is "easy"
+end
+
+end

test/scaling/nointerp.jl

@@ -0,0 +1,38 @@
+module ScalingNoInterpTests
+
+using Interpolations, Base.Test
+
+xs = -pi:2pi/10:pi
+f1(x) = sin(x)
+f2(x) = cos(x)
+f3(x) = sin(x) .* cos(x)
+f(x,y) = y == 1 ? f1(x) : (y == 2 ? f2(x) : (y == 3 ? f3(x) : error("invalid value for y (must be 1, 2 or 3, you used $y)")))
+ys = 1:3
+
+A = hcat(f1(xs), f2(xs), f3(xs))
+
+itp = interpolate(A, Tuple{BSpline(Quadratic(Periodic)), NoInterp}, OnGrid)
+sitp = scale(itp, xs, ys)
+
+for (ix,x0) in enumerate(xs[1:end-1]), y0 in ys
+    x,y = x0, y0
+    @test_approx_eq_eps sitp[x,y] f(x,y) .05
+end
+
+# Test error messages for incorrect initialization
+function message_is(message)
+	r -> r.err.msg == message || error("Incorrect error message: expected '$message' but was '$(r.err.msg)'")
+end
+Test.with_handler(message_is("Must scale 2-dimensional interpolation object with exactly 2 ranges (you used 1)")) do
+	@test scale(itp, xs)
+end
+Test.with_handler(message_is("NoInterp dimension 2 must be scaled with unit range 1:3")) do
+	@test scale(itp, xs, -1:1)
+end
+Test.with_handler(message_is("The length of the range in dimension 1 (8) did not equal the size of the interpolation object in that direction (11)")) do
+	@test scale(itp, -pi:2pi/7:pi, 1:3)
+end
+Test.with_handler(message_is("Must index into 2-dimensional scaled interpolation object with exactly 2 indices (you used 1)")) do
+	@test sitp[2.3]
+end
+end

test/scaling/runtests.jl

@@ -0,0 +1,4 @@
+include("scaling.jl")
+include("dimspecs.jl")
+include("nointerp.jl")
+include("withextrap.jl")

test/scaling/scaling.jl

@@ -0,0 +1,52 @@
+module ScalingTests
+
+using Interpolations
+using Base.Test
+
+# Model linear interpolation of y = -3 + .5x by interpolating y=x
+# and then scaling to the new x range
+
+itp = interpolate(1:1.0:10, BSpline(Linear), OnGrid)
+
+sitp = scale(itp, -3:.5:1.5)
+
+for (x,y) in zip(-3:.05:1.5, 1:.1:10)
+    @test_approx_eq sitp[x] y
+end
+
+# Verify that it works in >1D, with different types of ranges
+
+gauss(phi, mu, sigma) = exp(-(phi-mu)^2 / (2sigma)^2)
+testfunction(x,y) = gauss(x, 0.5, 4) * gauss(y, -.5, 2)
+
+xs = -5:.5:5
+ys = -4:.2:4
+zs = Float64[testfunction(x,y) for x in xs, y in ys]
+
+itp2 = interpolate(zs, BSpline(Quadratic(Flat)), OnGrid)
+sitp2 = scale(itp2, xs, ys)
+
+for x in xs, y in ys
+    @test_approx_eq testfunction(x,y) sitp2[x,y]
+end
+
+# Test gradients of scaled grids
+xs = -pi:.1:pi
+ys = sin(xs)
+itp = interpolate(ys, BSpline(Linear), OnGrid)
+sitp = scale(itp, xs)
+
+for x in -pi:.1:pi
+    g = @inferred(gradient(sitp, x))[1]
+    @test_approx_eq_eps cos(x) g .05
+end
+
+# Verify that return types are reasonable
+@inferred(getindex(sitp2, -3.4, 1.2))
+@inferred(getindex(sitp2, -3, 1))
+@inferred(getindex(sitp2, -3.4, 1))
+
+sitp32 = scale(interpolate(Float32[testfunction(x,y) for x in -5:.5:5, y in -4:.2:4], BSpline(Quadratic(Flat)), OnGrid), -5f0:.5f0:5f0, -4f0:.2f0:4f0)
+@test typeof(@inferred(getindex(sitp32, -3.4f0, 1.2f0))) == Float32
+
+end

test/scaling/withextrap.jl

@@ -0,0 +1,40 @@
+
+module ScalingWithExtrapTests
+
+using Interpolations, Base.Test
+
+xs = linspace(-5, 5, 10)
+ys = sin(xs)
+
+function run_tests{T,N,IT}(sut::Interpolations.AbstractInterpolation{T,N,IT,OnGrid}, itp)
+    for x in xs
+        @test_approx_eq_eps sut[x] sin(x) sqrt(eps(sin(x)))
+    end
+    @test sut[-5] == sut[-5.1] == sut[-15.8] == sut[-Inf] == itp[1]
+    @test sut[5] == sut[5.1] == sut[15.8] == sut[Inf] == itp[end]
+end
+
+function run_tests{T,N,IT}(sut::Interpolations.AbstractInterpolation{T,N,IT,OnCell}, itp)
+    halfcell = (xs[2] - xs[1]) / 2
+
+    for x in (5 + halfcell, 5 + 1.1halfcell, 15.8, Inf)
+        @test sut[-x] == itp[.5]
+        @test sut[x] == itp[end+.5]
+    end
+end
+
+for GT in (OnGrid, OnCell)
+    itp = interpolate(ys, BSpline(Quadratic(Flat)), GT)
+
+    # Test extrapolating, then scaling
+    eitp = extrapolate(itp, Flat)
+    seitp = scale(eitp, xs)
+    run_tests(seitp, itp)
+
+    # Test scaling, then extrapolating
+    sitp = scale(itp, xs)
+    esitp = extrapolate(sitp, Flat)
+    run_tests(esitp, itp)
+end
+
+end