WIP: Output type handling (fixes #17)

src/Interpolations.jl

@@ -31,82 +31,94 @@ export
 
     degree,
     boundarycondition,
-    gridrepresentation
+    gridrepresentation,
+
+    gradient,
+    gradient!
 
 abstract Degree{N}
 
 abstract GridRepresentation
-type OnGrid <: GridRepresentation end
-type OnCell <: GridRepresentation end
+immutable OnGrid <: GridRepresentation end
+immutable OnCell <: GridRepresentation end
 
 abstract BoundaryCondition
-type None <: BoundaryCondition end
-type Flat <: BoundaryCondition end
-type Line <: BoundaryCondition end
+immutable None <: BoundaryCondition end
+immutable Flat <: BoundaryCondition end
+immutable Line <: BoundaryCondition end
 typealias Natural Line
-type Free <: BoundaryCondition end
-type Periodic <: BoundaryCondition end
+immutable Free <: BoundaryCondition end
+immutable Periodic <: BoundaryCondition end
 immutable Reflect <: BoundaryCondition end
 
 abstract InterpolationType{D<:Degree,BC<:BoundaryCondition,GR<:GridRepresentation}
 
 include("extrapolation.jl")
 
 abstract AbstractInterpolation{T,N,IT<:InterpolationType,EB<:ExtrapolationBehavior} <: AbstractArray{T,N}
-type Interpolation{T,N,IT<:InterpolationType,EB<:ExtrapolationBehavior} <: AbstractInterpolation{T,N,IT,EB}
-    coefs::Array{T,N}
+type Interpolation{T,N,TCoefs,IT<:InterpolationType,EB<:ExtrapolationBehavior} <: AbstractInterpolation{T,N,IT,EB}
+    coefs::Array{TCoefs,N}
 end
-function Interpolation{T,N,IT<:InterpolationType,EB<:ExtrapolationBehavior}(A::Array{T,N}, it::IT, ::EB)
+function Interpolation{N,TCoefs,TWeights<:Real,IT<:InterpolationType,EB<:ExtrapolationBehavior}(::Type{TWeights}, A::AbstractArray{TCoefs,N}, it::IT, ::EB)
     isleaftype(IT) || error("The interpolation type must be a leaf type (was $IT)")
-
-    isleaftype(T) || warn("For performance reasons, consider using an array of a concrete type T (eltype(A) == $(eltype(A)))")
 
-    Interpolation{T,N,IT,EB}(prefilter(A,it))
+    isleaftype(TCoefs) || warn("For performance reasons, consider using an array of a concrete type T (eltype(A) == $(eltype(A)))")
+
+    c = one(TWeights)
+    for _ in 2:N
+        c *= c
+    end
+    T = typeof(c * one(TCoefs))
+
+    Interpolation{T,N,TCoefs,IT,EB}(prefilter(TWeights,A,it))
 end
+Interpolation(A::AbstractArray, it::InterpolationType, eb::ExtrapolationBehavior) = Interpolation(Float64, A, it, eb)
+Interpolation(A::AbstractArray{Float32}, it::InterpolationType, eb::ExtrapolationBehavior) = Interpolation(Float32, A, it, eb)
+Interpolation(A::AbstractArray{Rational{Int}}, it::InterpolationType, eb::ExtrapolationBehavior) = Interpolation(Rational{Int}, A, it, eb)
 
 # Unless otherwise specified, use coefficients as they are, i.e. without prefiltering
 # However, all prefilters copy the array, so do that here as well
-prefilter{T,N,IT<:InterpolationType}(A::AbstractArray{T,N}, ::IT) = copy(A)
+prefilter{TWeights,T,N,IT<:InterpolationType}(::Type{TWeights}, A::AbstractArray{T,N}, ::IT) = copy(A)
 
-size{T,N,IT<:InterpolationType}(itp::Interpolation{T,N,IT}, d::Integer) =
-    size(itp.coefs, d) - 2*padding(IT())
-size(itp::Interpolation) = tuple([size(itp,i) for i in 1:ndims(itp)]...)
+size{T,N}(itp::Interpolation{T,N}, d::Integer) = d > N ? 1 : size(itp.coefs, d) - 2*padding(interpolationtype(itp))
+size(itp::AbstractInterpolation) = tuple([size(itp,i) for i in 1:ndims(itp)]...)
 ndims(itp::Interpolation) = ndims(itp.coefs)
-eltype(itp::Interpolation) = eltype(itp.coefs)
 
 offsetsym(off, d) = off == -1 ? symbol(string("ixm_", d)) :
                     off ==  0 ? symbol(string("ix_", d)) :
                     off ==  1 ? symbol(string("ixp_", d)) :
                     off ==  2 ? symbol(string("ixpp_", d)) : error("offset $off not recognized")
 
+interpolationtype{T,N,TCoefs,IT<:InterpolationType}(itp::Interpolation{T,N,TCoefs,IT}) = IT()
+
 boundarycondition{D,BC<:BoundaryCondition}(::InterpolationType{D,BC}) = BC()
-boundarycondition{T,N,IT}(::Interpolation{T,N,IT}) = boundarycondition(IT())
 gridrepresentation{D,BC,GR<:GridRepresentation}(::InterpolationType{D,BC,GR}) = GR()
-gridrepresentation{T,N,IT}(::Interpolation{T,N,IT}) = gridrepresentation(IT())
 degree{D<:Degree,BC,GR}(::InterpolationType{D,BC,GR}) = D()
-degree{T,N,IT}(::Interpolation{T,N,IT}) = degree(IT())
+
+# If nothing else is specified, don't pad at all
+padding(::InterpolationType) = 0
+
+for op in (:boundarycondition, :gridrepresentation, :degree, :padding)
+    eval(:($(op)(itp::Interpolation) = $(op)(interpolationtype(itp))))
+end
 
 include("constant.jl")
 include("linear.jl")
 include("quadratic.jl")
 
-# If nothing else is specified, don't pad at all
-padding(::InterpolationType) = 0
-padding{T,N,IT<:InterpolationType}(::Interpolation{T,N,IT}) = padding(IT())
-
-function pad_size_and_index(sz::Tuple, pad)
+function padded_index(sz::Tuple, pad)
     sz = Int[s+2pad for s in sz]
     N = length(sz)
     ind = cell(N)
     for i in 1:N
         ind[i] = 1+pad:sz[i]-pad
     end
-    sz, ind
+    ind,sz
 end
 function copy_with_padding(A, it::InterpolationType)
     pad = padding(it)
-    sz,ind = pad_size_and_index(size(A), pad)
-    coefs = fill(convert(eltype(A), 0), sz...)
+    ind,sz = padded_index(size(A), pad)
+    coefs = zeros(eltype(A), sz...)
     coefs[ind...] = A
     coefs, pad
 end
@@ -139,11 +151,9 @@ for IT in (
         it = IT()
         eb = EB()
         gr = gridrepresentation(it)
-        eval(ngenerate(
-            :N,
-            :(promote_type(T, x...)),
-            :(getindex{T,N}(itp::Interpolation{T,N,$IT,$EB}, x::NTuple{N,Real}...)), 
-            N->quote
+
+        function getindex_impl(N)
+                quote
                 # Handle extrapolation, by either throwing an error,
                 # returning a value, or shifting x to somewhere inside
                 # the domain
@@ -160,13 +170,38 @@ for IT in (
                 @inbounds ret = $(index_gen(degree(it), N))
                 ret
             end
+        end
+
+        # Resolve ambiguity with linear indexing,
+        #   getindex{T,N}(A::AbstractArray{T,N}, i::Real)
+        eval(ngenerate(:N, :T, :(getindex{T,N,TCoefs}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::Real)),
+            N->:(error("Linear indexing is not supported for interpolation objects"))
         ))
 
-        eval(ngenerate(
-            :N,
-            :(Array{promote_type(T,typeof(x)...),1}),
-            :(gradient!{T,N}(g::Array{T,1}, itp::Interpolation{T,N,$IT,$EB}, x::NTuple{N,Real}...)),
+        # Resolve ambiguity with real multilinear indexing
+        #   getindex{T,N}(A::AbstractArray{T,N}, NTuple{N,Real}...)
+        eval(ngenerate(:N, :T, :(getindex{T,N,TCoefs}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::NTuple{N,Real}...)), getindex_impl))
+
+        # Resolve ambiguity with general array indexing,
+        #   getindex{T,N}(A::AbstractArray{T,N}, I::AbstractArray{T,N})
+        eval(ngenerate(:N, :T, :(getindex{T,N,TCoefs}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, I::AbstractArray{T})),
+            N->:(error("Array indexing is not defined for interpolation objects."))
+        ))
+
+        # Resolve ambiguity with colon indexing for 1D interpolations
+        #   getindex{T}(A::AbstractArray{T,1}, C::Colon)
+        eval(ngenerate(:N, :T, :(getindex{T,TCoefs}(itp::Interpolation{T,1,TCoefs,$IT,$EB}, c::Colon)),
+            N->:(error("Colon indexing is not supported for interpolation objects"))
+        ))
+
+        # Allow multilinear indexing with any types
+        eval(ngenerate(:N, :(promote_type(T,TIndex)), :(getindex{T,N,TCoefs,TIndex}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::NTuple{N,TIndex}...)), getindex_impl))
+
+        # Define in-place gradient calculation
+        #   gradient!(g::Vector, itp::Interpolation, NTuple{N}...)
+        eval(ngenerate(:N, :(Vector{TOut}), :(gradient!{T,N,TCoefs,TOut}(g::Vector{TOut}, itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::NTuple{N,Any}...)),
             N->quote
+                length(g) == $N || error("g must be an array with exactly N elements (length(g) == "*string(length(g))*", N == "*string(N)*")")
                 $(extrap_transform_x(gr,eb,N))
                 $(define_indices(it,N))
                 @nexprs $N dim->begin
@@ -184,7 +219,7 @@ for IT in (
     end
 end
 
-gradient{T}(itp::Interpolation{T}, x...) = gradient!(Array(T,ndims(itp)), itp, x...)
+gradient{T,N}(itp::Interpolation{T,N}, x...) = gradient!(Array(T,N), itp, x...)
 
 # This creates prefilter specializations for all interpolation types that need them
 for IT in (
@@ -199,8 +234,8 @@ for IT in (
         Quadratic{Periodic,OnGrid},
         Quadratic{Periodic,OnCell},
     )
-    @ngenerate N promote_type_grid(T, x...) function prefilter{T,N}(A::Array{T,N},it::IT)
-        ret, pad = copy_with_padding(A,it)
+    @ngenerate N Array{TWeights,N} function prefilter{TWeights,TCoefs,N}(::Type{TWeights}, A::Array{TCoefs,N}, it::IT)
+        ret, pad = copy_with_padding(A, it)
 
         szs = collect(size(ret))
         strds = collect(strides(ret))
@@ -209,7 +244,7 @@ for IT in (
             n = szs[dim]
             szs[dim] = 1
 
-            M, b = prefiltering_system(eltype(A), n, it)
+            M, b = prefiltering_system(TWeights, TCoefs, n, it)
 
             @nloops N i d->1:szs[d] begin
                 cc = @ntuple N i

src/constant.jl

@@ -1,9 +1,9 @@
-type ConstantDegree <: Degree{0} end
-type Constant{GR<:GridRepresentation} <: InterpolationType{ConstantDegree,None,GR} end
+immutable ConstantDegree <: Degree{0} end
+immutable Constant{GR<:GridRepresentation} <: InterpolationType{ConstantDegree,None,GR} end
 Constant{GR<:GridRepresentation}(::GR) = Constant{GR}()
 
 function define_indices(::Constant, N)
-    :(@nexprs $N d->(ix_d = clamp(round(Int,x_d), 1, size(itp,d))))
+    :(@nexprs $N d->(ix_d = clamp(round(real(x_d)), 1, size(itp,d))))
 end
 
 function coefficients(c::Constant, N)
@@ -12,12 +12,12 @@ end
 
 function coefficients(::Constant, N, d)
     sym, symx = symbol(string("c_",d)), symbol(string("x_",d))
-    :($sym = one(typeof($symx)))
+    :($sym = 1)
 end
 
 function gradient_coefficients(::Constant, N, d)
     sym, symx = symbol(string("c_",d)), symbol(string("x_",d))
-    :($sym = zero(typeof($symx)))
+    :($sym = 0)
 end
 
 function index_gen(degree::ConstantDegree, N::Integer, offsets...)

src/extrapolation.jl

@@ -1,42 +1,42 @@
 abstract ExtrapolationBehavior
 
-type ExtrapError <: ExtrapolationBehavior end
+immutable ExtrapError <: ExtrapolationBehavior end
 function extrap_transform_x(::OnGrid, ::ExtrapError, N)
     quote
-        @nexprs $N d->(1 <= x_d <= size(itp,d) || throw(BoundsError()))
+        @nexprs $N d->(1 <= real(x_d) <= size(itp,d) || throw(BoundsError()))
     end
 end
 function extrap_transform_x(::OnCell, ::ExtrapError, N)
     quote
-        @nexprs $N d->(.5 <= x_d <= size(itp,d) + .5 || throw(BoundsError()))
+        @nexprs $N d->(1//2 <= real(x_d) <= size(itp,d) + 1//2 || throw(BoundsError()))
     end
 end
 
-type ExtrapNaN <: ExtrapolationBehavior end
+immutable ExtrapNaN <: ExtrapolationBehavior end
 function extrap_transform_x(::OnGrid, ::ExtrapNaN, N)
     quote
-        @nexprs $N d->(1 <= x_d <= size(itp,d) || return convert(T, NaN))
+        @nexprs $N d->(1 <= real(x_d) <= size(itp,d) || return convert(T, NaN))
     end
 end
 function extrap_transform_x(::OnCell, ::ExtrapNaN, N)
     quote
-        @nexprs $N d->(.5 <= x_d <= size(itp,d) + .5 || return convert(T, NaN))
+        @nexprs $N d->(1//2 <= real(x_d) <= size(itp,d) + 1//2 || return convert(T, NaN))
     end
 end
 
-type ExtrapConstant <: ExtrapolationBehavior end
+immutable ExtrapConstant <: ExtrapolationBehavior end
 function extrap_transform_x(::OnGrid, ::ExtrapConstant, N)
     quote
         @nexprs $N d->(x_d = clamp(x_d, 1, size(itp,d)))
     end
 end
 function extrap_transform_x(::OnCell, ::ExtrapConstant, N)
     quote
-        @nexprs $N d->(x_d = clamp(x_d, .5, size(itp,d)+.5))
+        @nexprs $N d->(x_d = clamp(x_d, 1//2, size(itp,d)+1//2))
     end
 end
 
-type ExtrapReflect <: ExtrapolationBehavior end
+immutable ExtrapReflect <: ExtrapolationBehavior end
 function extrap_transform_x(::OnGrid, ::ExtrapReflect, N)
     quote
         @nexprs $N d->begin
@@ -80,23 +80,12 @@ function extrap_transform_x(::OnCell, ::ExtrapReflect, N)
     end
 end
 
-type ExtrapPeriodic <: ExtrapolationBehavior end
+immutable ExtrapPeriodic <: ExtrapolationBehavior end
 function extrap_transform_x(::GridRepresentation, ::ExtrapPeriodic, N)
-    quote
-        @nexprs $N d->begin
-            # translate x_d to inside the domain
-            n = convert(typeof(x_d), size(itp,d))
-            while x_d < .5
-                x_d += n
-            end
-            while x_d >= n + .5
-                x_d -= n
-            end
-        end
-    end
+    :(@nexprs $N d->(x_d = mod1(x_d, size(itp,d))))
 end
 
-type ExtrapLinear <: ExtrapolationBehavior end
+immutable ExtrapLinear <: ExtrapolationBehavior end
 function extrap_transform_x(::OnGrid, ::ExtrapLinear, N)
     quote
         @nexprs $N d->begin

src/linear.jl

@@ -1,13 +1,13 @@
-type LinearDegree <: Degree{1} end
-type Linear{GR<:GridRepresentation} <: InterpolationType{LinearDegree,None,GR} end
+immutable LinearDegree <: Degree{1} end
+immutable Linear{GR<:GridRepresentation} <: InterpolationType{LinearDegree,None,GR} end
 Linear{GR<:GridRepresentation}(::GR) = Linear{GR}()
 
 function define_indices(::Linear, N)
     quote
         @nexprs $N d->begin
-            ix_d = clamp(floor(Int,x_d), 1, size(itp,d)-1)
+            ix_d = clamp(floor(real(x_d)), 1, size(itp,d)-1)
             ixp_d = ix_d + 1
-            fx_d = x_d - convert(typeof(x_d), ix_d)
+            fx_d = x_d - ix_d
         end
     end
 end
@@ -19,16 +19,16 @@ end
 function coefficients(::Linear, N, d)
     sym, symp, symfx = symbol(string("c_",d)), symbol(string("cp_",d)), symbol(string("fx_",d))
     quote
-        $sym = one(typeof($symfx)) - $symfx
+        $sym = 1 - $symfx
         $symp = $symfx
     end
 end
 
 function gradient_coefficients(::Linear,N,d)
     sym, symp, symfx = symbol(string("c_",d)), symbol(string("cp_",d)), symbol(string("fx_",d))
     quote
-        $sym = -one(typeof($symfx))
-        $symp = one(typeof($symfx))
+        $sym = -1
+        $symp = 1
     end
 end