Rewrite with clear separation of concerns (closes #2)

doc/Extrapolation.md

@@ -0,0 +1,44 @@
+# Extrapolation behavior
+
+## 1. Simple extrapolation
+
+Extrapolation in `Interpolations.jl` takes several forms. The simplest, is when the extrapolation behavior can be described with the following pseudo-code:
+
+    if (the coordinate is outside the domain) then
+        do something that is
+            * well defined without looking at the data set
+            * decides the outcome (error/return value) of the indexing operation
+    end
+
+    proceed to inbounds interpolation
+
+An example of this interpolation behavior is `ExtrapError`, which simply throws a bounds error if the coordinate is outside the domain. `Interpolations.jl` could support, for example, the following variants:
+
+* `ExtrapError`: Throws a `BoundsError`, just like `Grid.jl`
+* `ExtrapNaN`: Returns `convert(T, NaN)`, where `T` is `eltype(data)`
+* `ExtrapNull`: Returns a value-less `Nullabel{T}`
+
+## 2. Index transformation extrapolation
+
+The next form is index transformation, which can be described as
+
+    if (the coordinate is outside the domain) then
+        calculate an index that is inside the domain, which gives
+        the extrapolated value
+    end
+
+    proceed to inbounds interpolation (using the transformed index)
+
+An example here is `ExtrapPeriodic`, which transforms the coordinate index to one which is inside the domain by means of modulo calculations. Another example is `ExtrapConstant`, which clamps out-of-bounds coordinates to their nearest inbounds data point.
+
+For some of these, extra care needs to be taken in higher dimensions, when deciding what happens in the "outside corners":
+
+              |               |     what happens here?
+              |               |
+    ----------+---------------+--------
+              |               |
+              |   the domain  |
+              |               |
+    ----------+---------------+--------
+              |               |
+              |               |    ...and here? 

doc/Math.md

@@ -0,0 +1,64 @@
+# The math behind Interpolations.jl
+
+Interpolations.jl uses [B-splines](http://en.wikipedia.org/wiki/B-spline#Definition), or *basis splines*, to interpolate a data set with a specified degree. In short, the idea is to use a set of basis functions with limited support, and in each interval between two data points, construct a linear combination of the basis functions which have support in that interval, each yielding a piece in a piecewise defined function with support across the entire data set. B-splines are defined to have the maximum degree of continuity, with minimal support. Specific details on the general mathematics of such interpolations can [this paper](http://dx.doi.org/10.1109/42.875199).
+
+## The purpose of this document
+
+There are a few core concepts one must understand in order to use `Interpolations.jl`:
+
+1. Interpolation degree
+2. Boundary conditions
+3. Mid-point and on-grid interpolation
+3. Extrapolation
+
+Defining these three concepts, and how `Interpolations.jl` reasons about them, is the main focus of this document.
+
+## 1. Interpolation degree
+
+The interpolation degree decides what continuity properties the interpolant (i.e. the object which represents the interpolated data set) has; for example, a linear interpolation is a piecewise linear function, which is continuous but has discontinuous derivatives (or gradients, in higher dimensions).
+
+Because of how B-splines are defined, the interpolation degree also decides the *support* of the *spline functions*. For a linear interpolation, the function value is decided by the value at *two* data points, so its support is two<sup>1.</sup>.
+
+## 2. Boundary conditions
+
+For higher interpolation degrees (specifically, from quadratic interpolation and up), the support of the B-splines is too large for the interpolating scheme to be able to figure out the values near the edges of the data sets. In order to close the equation systems at the edges, boundary conditions are used.
+
+For quadratic interpolation, for example, a common boundary condition is to assume that the function is flat at the edges (i.e. the derivative there is 0). This lets us introduce an extra equation at each edge, through a finite approximation of the derivative, which closes the system. Another common way of terminating the interpolation is to extend the second-to-outermost all the way to the edge of the data set.
+
+## 3. Mid-point and on-grid interpolation
+
+In any discrete data representation, there is a *cell* associated with each data point. Depending on the application, it may make sense to consider the data points to represent either the *center* of the cells, or the *edges*. `Interpolations.jl` will support both of these, using the term *midpoint interpolation* for data sets where the data points are in the middle of the cells, and *on-grid interpolation* when the data points and cell boundaries coincide.
+
+## Interlude: the `Interpolation` type hierarchy
+
+All the above concepts fundamentally affect the behavior of an interpolating function, and in `Interpolations.jl` they therefore each have a representation in the type hierarchy.
+
+Different types of interpolations are separated by *type parameters*, and each type of interpolation is a concrete type descending from
+
+    abstract InterpolationType{D<:Degree,BC<:BoundaryCondition,G<:GridRepresentation}
+
+For example, a quadratic on-grid implementation with flat boundary conditions, is represented by an `Interpolation{Quadratic,Flat,OnGrid}`, where `Quadratic`, `Flat` and `OnGrid` are in turn concrete types implementing the abstract types indicated above.
+
+## 4. Extrapolation behavior
+
+Somewhat orthogonal<sup>2</sup> to the concepts outlined above, is the concept of *extrapolation*, i.e. evaluation of the interpolant outside the domain defined by the data set. For some types of extrapolation, this behavior is defined by a translation of the interpolation coordinate to somewhere inside the domain (e.g. periodic or reflecting boundaries), while for other types it entails a separate calculation entirely (e.g. linear extrapolation).
+
+[Read more](/doc/Extrapolation.md)
+
+## Supporting an `Interpolation` type in `Interpolations.jl`
+
+Not all combinations of the above four concepts are supported in the library; for example, neither constant nor linear on-grid interpolations need andy boundary conditions, and therefore they don't support any.
+
+An interpolation is represented by an object of a concrete type descending from
+
+    abstract Interpolation{IT<:InterpolationType,EB<:ExtrapolationBehavior}
+
+----
+
+<sup>
+
+1. For this reason, linear B-splines are often referred to as *2nd order*, which may be a source of confusion since the interpolating function itself is linear, i.e. of first order. In `Interpolations.jl`, we will try to avoid this confusion by referring to interpolation degree by "linear", "quadratic" etc.
+
+2. Although the separation of concepts 1-3 from extrapolation behavior is computationally sound, it does allow for some interesting, yet probably nonsensical, combinations of interpolation degrees, boundary conditions and extrapolation behavior. One could imagine for example a constant interpolation (which needs no boundary condition) with linear extrapolation, in which case the interpolating function is a sequence of "steps" with 0-derivative inside the domain, while suddenly having a nonzero derivative outside. Similarly, in most cases where the extrapolation behavior is defined as constant or reflecting, it will make sense to specify matching boundary conditions, but other combinations are entirely supported by `Interpolations.jl`; your milage may vary, very much...
+
+</sup>

src/Interpolations.jl

@@ -4,61 +4,126 @@ using Base.Cartesian
 
 import Base: size, eltype, getindex
 
-export BoundaryCondition,
-    BCnone,
-
-    ExtrapolationBehavior,
+export 
+    Interpolation,
+    Constant,
+    Linear,
+    Quadratic,
     ExtrapError,
     ExtrapNaN,
-    ExtrapConstant,
-    ExtrapPeriodic,
+    OnCell,
+    OnGrid,
+    ExtendInner
 
-    InterpolationDegree,
-    Linear,
+abstract Degree{N}
 
-    AbstractInterpolation,
-    Interpolation
+abstract GridRepresentation
+type OnGrid <: GridRepresentation end
+type OnCell <: GridRepresentation end
 
 abstract BoundaryCondition
-immutable BCnone <: BoundaryCondition end
+type None <: BoundaryCondition end
+type ExtendInner <: BoundaryCondition end
 
-abstract ExtrapolationBehavior
-immutable ExtrapError <: ExtrapolationBehavior end
-immutable ExtrapNaN <: ExtrapolationBehavior end
-immutable ExtrapConstant <: ExtrapolationBehavior end
-immutable ExtrapPeriodic <: ExtrapolationBehavior end
+abstract InterpolationType{D<:Degree,BC<:BoundaryCondition,GR<:GridRepresentation}
 
-abstract InterpolationDegree
-# Should we call these TwoPoint, ThreePoint, FourPoint?
-immutable Linear <: InterpolationDegree end
-immutable Quadratic <: InterpolationDegree end
-immutable Cubic <: InterpolationDegree end
+include("extrapolation.jl")
 
-abstract AbstractInterpolation{T,N,D,BC<:BoundaryCondition,EB<:ExtrapolationBehavior} <: AbstractArray{T,N}
+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}
+end
+function Interpolation{T,N,IT<:InterpolationType,EB<:ExtrapolationBehavior}(A::Array{T,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)))")
 
-type Interpolation{T,N,ID<:InterpolationDegree,BC<:BoundaryCondition,EB<:ExtrapolationBehavior} <: AbstractInterpolation{T,N,ID,BC,EB}
-	coefs::Array{T,N}
+    Interpolation{T,N,IT,EB}(prefilter(A,it))
 end
-Interpolation{T,N,ID<:InterpolationDegree,BC<:BoundaryCondition,EB<:ExtrapolationBehavior}(A::Array{T,N}, ::Type{ID}, ::Type{BC}, ::Type{EB}) = Interpolation{T,N,ID,BC,EB}(A)
 
-include("linear.jl")
+# 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)
 
 size(itp::Interpolation, d::Integer) = size(itp.coefs, d)
 size(itp::Interpolation) = size(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")
+
+gridrepresentation{D,BC,GR<:GridRepresentation}(::InterpolationType{D,BC,GR}) = GR()
+degree{D<:Degree,BC,GR}(::InterpolationType{D,BC,GR}) = D()
+
+include("constant.jl")
+include("linear.jl")
+include("quadratic.jl")
+
 promote_type_grid(T, x...) = promote_type(T, typeof(x)...)
 
+# This function gets specialized versions for interpolation types that need prefiltering
+prefilter(A::Array, it::InterpolationType) = copy(A)
+
 # This creates getindex methods for all supported combinations
-for ID in (Linear,)
-    # for BC in subtypes(BoundaryCondition)
-    for BC in (BCnone,)
-        for EB in (ExtrapError,ExtrapNaN,ExtrapConstant,ExtrapPeriodic)
-            eval(ngenerate(:N, :(promote_type_grid(T, x...)), :(getindex{T,N}(itp::Interpolation{T,N,$ID,$BC,$EB}, x::NTuple{N,Real}...)),
-                      N->body_gen(ID, BC, EB, N)))
-        end
+for IT in (Constant{OnCell},Linear{OnGrid},Quadratic{ExtendInner,OnCell})
+    for EB in (ExtrapError,ExtrapNaN)
+        it = IT()
+        eb = EB()
+        gr = gridrepresentation(it)
+        eval(ngenerate(
+            :N,
+            :(promote_type_grid(T, x...)),
+            :(getindex{T,N}(itp::Interpolation{T,N,$IT,$EB}, x::NTuple{N,Real}...)), 
+            N->quote
+                $(extrap_gen(gr,eb,N))
+
+                # If the boundary condition mandates separate treatment, this is done
+                # by bc_gen.
+                # Given an interpolation object itp, with N dimensions, and a coordinate
+                # x_d, it should define ix_d such that all the coefficients required by
+                # the interpolation will be inbounds.
+                $(bc_gen(it, N))
+
+                # indices calculates the indices required for this interpolation degree,
+                # based on ix_d defined by bc_gen(), as well as the distance fx_d from 
+                # the cell index ix_d to the interpolation coordinate x_d
+                $(indices(degree(it), N))
+
+                # These coefficients determine the interpolation basis expansion
+                $(coefficients(degree(it), N))
+
+                @inbounds ret = $(index_gen(degree(it), N))
+                ret
+            end
+        ))
     end
 end
 
+# This creates prefilter specializations for all interpolation types that need them
+for IT in (Quadratic{ExtendInner,OnCell},)
+    @ngenerate N promote_type_grid(T, x...) function prefilter{T,N}(A::Array{T,N},it::IT)
+        ret = similar(A)
+        szs = collect(size(A))
+        strds = collect(strides(A))
+
+        for dim in 1:N
+            n = szs[dim]
+            szs[dim] = 1
+
+            M = prefiltering_system_matrix(eltype(A), n, it)
+
+            @nloops N i d->1:szs[d] begin
+                cc = @ntuple N i
+                strt = 1 + sum([(cc[i]-1)*strds[i] for i in 1:length(cc)])
+                rng = range(strt, strds[dim], n)
+                ret[rng] = M \ vec(A[rng])
+            end            
+            szs[dim] = n
+        end
+        ret
+    end
+end
 
 end # module

src/constant.jl

@@ -0,0 +1,33 @@
+type ConstantDegree <: Degree{0} end
+type Constant{GR<:GridRepresentation} <: InterpolationType{ConstantDegree,None,GR} end
+
+Constant{GR<:GridRepresentation}(::GR) = Constant{GR}()
+
+function bc_gen{IT<:Constant}(::IT, N)
+    quote
+        @nexprs $N d->(ix_d = iround(x_d))
+    end
+end
+
+function indices(::ConstantDegree, N)
+    quote
+        # Constant interpolation doesn't need an fx_d
+    end
+end
+
+function coefficients(::ConstantDegree, N)
+    quote
+        @nexprs $N d->(c_d = one(typeof(x_d)))
+    end
+end
+
+function index_gen(degree::ConstantDegree, N::Integer, offsets...)
+    if (length(offsets) < N)
+        d = length(offsets)+1
+        sym = symbol("c_"*string(d))
+        return :($sym * $(index_gen(degree, N, offsets..., 0)))
+    else
+        indices = [offsetsym(offsets[d], d) for d = 1:N]
+        return :(itp.coefs[$(indices...)])
+    end
+end

src/extrapolation.jl

@@ -0,0 +1,19 @@
+
+abstract ExtrapolationBehavior
+type ExtrapError <: ExtrapolationBehavior end
+
+function extrap_gen(::OnGrid, ::ExtrapError, N)
+    quote
+        @nexprs $N d->(1 <= x_d <= size(itp,d) || throw(BoundsError()))
+    end
+end
+extrap_gen(::OnCell, e::ExtrapError, N) = extrap_gen(OnGrid(), e, N)
+
+type ExtrapNaN <: ExtrapolationBehavior end
+
+function extrap_gen(::OnGrid, ::ExtrapNaN, N)
+    quote
+        @nexprs $N d->(1 <= x_d <= size(itp,d) || return convert(T, NaN))
+    end
+end
+extrap_gen(::OnCell, e::ExtrapNaN, N) = extrap_gen(OnGrid(), e, N)