First attempt at handling eltypes better

Author: Tomas Lycken

src/Interpolations.jl

@@ -51,26 +51,34 @@ abstract InterpolationType{D<:Degree,BC<:BoundaryCondition,GR<:GridRepresentatio
 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{TEl,N,TCoefs<:Real,IT<:InterpolationType,EB<:ExtrapolationBehavior} <: AbstractInterpolation{TEl,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{TIn,N,TCoefs,IT<:InterpolationType,EB<:ExtrapolationBehavior}(::Type{TCoefs}, A::AbstractArray{TIn,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)))")
+    isleaftype(TIn) || 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))
+    c = one(TCoefs)
+    for _ in 2:N
+        c *= c
+    end
+    TEl = typeof(c)
+
+    Interpolation{TEl,N,TCoefs,IT,EB}(prefilter(TCoefs,A,it))
 end
+Interpolation{TIn,N,IT<:InterpolationType,EB<:ExtrapolationBehavior}(A::AbstractArray{TIn,N}, it::IT, eb::EB) = Interpolation(TIn, 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)
+# We also ensure that the coefficient array is of the correct type
+prefilter{T,N,TCoefs,IT<:InterpolationType}(::Type{TCoefs}, A::AbstractArray{T,N}, ::IT) = copy!(Array(TCoefs,size(A)...), 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,TCoefs,IT<:InterpolationType}(itp::Interpolation{T,N,TCoefs,IT}, d::Integer) = size(itp.coefs, d) - 2*padding(TCoefs,IT())
+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)
+eltype{T}(itp::Interpolation{T}) = T
+coeftype(itp::Interpolation) = eltype(itp.coefs)
 
 offsetsym(off, d) = off == -1 ? symbol(string("ixm_", d)) :
                     off ==  0 ? symbol(string("ix_", d)) :
@@ -89,8 +97,8 @@ 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())
+padding{T}(::Type{T}, ::InterpolationType) = zero(T)
+padding{T,N,IT<:InterpolationType}(::Interpolation{T,N,IT}) = padding(T,IT())
 
 function pad_size_and_index(sz::Tuple, pad)
     sz = Int[s+2pad for s in sz]
@@ -101,10 +109,10 @@ function pad_size_and_index(sz::Tuple, pad)
     end
     sz, ind
 end
-function copy_with_padding(A, it::InterpolationType)
+function copy_with_padding(TCoefs, A, it::InterpolationType)
     pad = padding(it)
     sz,ind = pad_size_and_index(size(A), pad)
-    coefs = fill(convert(eltype(A), 0), sz...)
+    coefs = fill!(Array(TCoefs, sz...), 0)
     coefs[ind...] = A
     coefs, pad
 end
@@ -137,8 +145,8 @@ for IT in (
         gr = gridrepresentation(it)
         eval(ngenerate(
             :N,
-            :(promote_type(T, x...)),
-            :(getindex{T,N}(itp::Interpolation{T,N,$IT,$EB}, x::NTuple{N,Real}...)), 
+            :T,
+            :(getindex{T,N,TCoefs<:Real}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::NTuple{N,TCoefs}...)), 
             N->quote
                 # Handle extrapolation, by either throwing an error,
                 # returning a value, or shifting x to somewhere inside
@@ -157,11 +165,14 @@ for IT in (
                 ret
             end
         ))
+        @ngenerate N T function getindex{T,N,TCoefs<:Real}(itp::Interpolation{T,N,TCoefs,IT,EB}, xs::NTuple{N,Real}...)
+            getindex(itp, [convert(TCoefs, x) for x in xs]...)
+        end
 
         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}...)),
+            :T,
+            :(gradient!{T,N,TCoefs}(g::Array{T,1}, itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::NTuple{N,TCoefs}...)),
             N->quote
                 $(extrap_transform_x(gr,eb,N))
                 $(define_indices(it,N))
@@ -180,7 +191,7 @@ for IT in (
     end
 end
 
-gradient{T}(itp::Interpolation{T}, x...) = gradient!(Array(T,ndims(itp)), itp, x...)
+gradient1{T}(itp::Interpolation{T,1}, x) = gradient!(Array(T,1),itp,x)[1]
 
 # This creates prefilter specializations for all interpolation types that need them
 for IT in (
@@ -193,8 +204,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 promote_type_grid(T, x...) function prefilter{T,N,TCoefs}(::Type{TCoefs}, A::Array{T,N},it::IT)
+        ret, pad = copy_with_padding(TCoefs, A,it)
 
         szs = collect(size(ret))
         strds = collect(strides(ret))
@@ -203,7 +214,7 @@ for IT in (
             n = szs[dim]
             szs[dim] = 1
 
-            M, b = prefiltering_system(eltype(A), n, it)
+            M, b = prefiltering_system(TCoefs, n, it)
 
             @nloops N i d->1:szs[d] begin
                 cc = @ntuple N i

src/extrapolation.jl

@@ -82,18 +82,7 @@ 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
 
 immutable ExtrapLinear <: ExtrapolationBehavior end

src/linear.jl

@@ -5,9 +5,9 @@ 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)
-            ixp_d = ix_d + 1
-            fx_d = x_d - convert(typeof(x_d), ix_d)
+            ix_d = clamp(floor(x_d), 1, size(itp,d)-1)
+            ixp_d = ix_d + one(typeof(ix_d))
+            fx_d = x_d - ix_d
         end
     end
 end

src/quadratic.jl

@@ -5,25 +5,25 @@ Quadratic{BC<:BoundaryCondition,GR<:GridRepresentation}(::BC, ::GR) = Quadratic{
 
 function define_indices(q::Quadratic, N)
     quote
-        pad = padding($q)
+        pad = padding(TCoefs, $q)
         @nexprs $N d->begin
-            ix_d = clamp(round(Integer, x_d), 1, size(itp,d)) + pad
+            ix_d = clamp(round(x_d), 1, size(itp,d)) + pad
             ixp_d = ix_d + 1
             ixm_d = ix_d - 1
 
-            fx_d = x_d - convert(typeof(x_d), ix_d - pad)
+            fx_d = x_d - (ix_d - pad)
         end
     end
 end
 function define_indices(q::Quadratic{Periodic}, N)
     quote
-        pad = padding($q)
+        pad = padding(TCoefs, $q)
         @nexprs $N d->begin
-            ix_d = clamp(round(Integer, x_d), 1, size(itp,d)) + pad
-            ixp_d = mod1(ix_d + 1, size(itp,d))
-            ixm_d = mod1(ix_d - 1, size(itp,d))
+            ix_d = clamp(round(x_d), 1, size(itp,d)) + pad
+            ixp_d = mod1(ix_d + one(typeof(ix_d)), size(itp,d))
+            ixm_d = mod1(ix_d - one(typeof(ix_d)), size(itp,d))
 
-            fx_d = x_d - convert(typeof(x_d), ix_d - pad)
+            fx_d = x_d - (ix_d - pad)
         end
     end
 end
@@ -36,19 +36,19 @@ function coefficients(q::Quadratic, N, d)
     symm, sym, symp =  symbol(string("cm_",d)), symbol(string("c_",d)), symbol(string("cp_",d))
     symfx = symbol(string("fx_",d))
     quote
-        $symm = .5 * ($symfx - .5)^2
-        $sym  = .75 - $symfx^2
-        $symp = .5 * ($symfx + .5)^2
+        $symm = convert(TCoefs, .5) * ($symfx - convert(TCoefs, .5))^2
+        $sym  = convert(TCoefs, .75) - $symfx^2
+        $symp = convert(TCoefs, .5) * ($symfx + convert(TCoefs, .5))^2
     end
 end
 
 function gradient_coefficients(q::Quadratic, N, d)
     symm, sym, symp =  symbol(string("cm_",d)), symbol(string("c_",d)), symbol(string("cp_",d))
     symfx = symbol(string("fx_",d))
     quote
-        $symm = $symfx-.5
-        $sym = -2*$symfx
-        $symp = $symfx+.5
+        $symm = $symfx - convert(TCoefs, .5)
+        $sym = -convert(TCoefs, 2) * $symfx
+        $symp = $symfx + convert(TCoefs, .5)
     end
 end
 
@@ -69,10 +69,10 @@ end
 # Quadratic interpolation has 1 extra coefficient at each boundary
 # Therefore, make the coefficient array 1 larger in each direction,
 # in each dimension.
-padding(::Quadratic) = 1
+padding{T}(::Type{T}, ::Quadratic) = one(T)
 # For periodic boundary conditions, we don't pad - instead, we wrap the
 # the coefficients
-padding(::Quadratic{Periodic}) = 0
+padding{T}(::Type{T}, ::Quadratic{Periodic}) = zero(T)
 
 function inner_system_diags{T}(::Type{T}, n::Int, ::Quadratic)
     du = fill(convert(T,1/8), n-1)