Specifying the output type of an interpolation object

[tomasaschan]

One of the issues with Grid.jl was that it's hard to specify what type the interpolation object should return - Interpolations.jl should be able to do better if we want it to be easily usable with any number type.

I'm leaning toward an interface with an optional first type argument, that specifies the output type, similar to what round and friends do in 0.4:

julia> itp = Interpolation(Float64, 1:10, Linear(OnGrid()), ExtrapNaN())

julia> typeof(itp[3.2])
Float64

Of course, we need to provide some sensible default, and there are a couple of options:

• Use the element type of the input array, i.e. Interpolation(A, config...) = Interpolation(eltype(A), A, config...)
  Pros:

  • very predictable for the end user; effectively, eltype(Interpolation(A,...)) == eltype(A)

  Cons:

  • will have the same problems as Grid.jl, especially for integer inputs (but these are avoided if the user explicitly sets a type)
  • we have to convert each return value, so there must be a conversion defined from whatever our calculations promote to, to the return type.

• Use use some sensible promotion result, e.g. promote_type(eltype(A), Float64, typeof(x)...) (where the Float64 is from the literals in the prefiltering step - we might be able to avoid those, if we can require a convert method for converting floats to eltype(A))

  Pros:

  • very similar to what we do today
  • won't require anything extra from the user to work with integers, e.g. itp = Interpolation(1:10, ...)

  Cons:

  • might be difficult to get it right, since there are lots of places where numbers enter the calculations
  • might seem unpredictable to the user, especially since e.g. promote_type(Rational{Int}, Float64) == Float64 will make us return floats instead of rationals for Interpolation(A::Array{T<:Rational}, ...)

I honestly don't know what's best here, but before we start implementing this, we need some sort of decisions on:

• What should the API look like? (Is the above suggestion good or ugly?)
• How should we choose a default?

[timholy]

Good issue to be thinking about now. +1 on having the first argument be the type, and for making it optional. This is going to be a bit rambling, as I'm thinking this through as I type.

For thinking this through carefully, a good exercise is to give physical units to everything (https://github.com/Keno/SIUnits.jl). Suppose the user supplies the input array in terms meters, based on a Float32 val type. In this case I presume we want to return the same type, even if the user asks for itp[3.2] rather than itp[3.2f0]. I suspect the best approach would be to pre-convert all input coordinates to Float32---note that we do not convert them to Meter{Float32}, as that would give Meter^2 as the natural result of 1d interpolation calculations.

Note that converting the input coordinates to a specified type is different (and a better choice, I think) from converting the result: it means that we don't have to worry about putting in so many manual conversions, and all the calculations should just "work out."

One nasty item, though: what to do if one(eltype(A)) doesn't exist? Think about an input array of matrices. We might want to make the coefficient type something that can also be specified manually, and include it in the Interpolation type parameters, just so users can specify it directly in case of trouble (and to reduce the number of code locations where we have to "virtually" compute it). One could also use this to ensure that all internal computations remain Rational, for example.

If the user asks for something for which the appropriate conversions are absent or impossible, with these steps I think it's no longer our problem, except perhaps to think about what the error message should look like (if indeed we want to customize it---perhaps a missing method error is all we need).

However, if we let users supply both, then there's a risk of broken promises:

element_type = Float64
coefficient_type = Rational
A = [1,2,3]
itp = Interpolation(element_type, coefficient_type, A)
julia> itp[2.2]
2476979795053773//1125899906842624  # not a Float64

because internally it looks like this:

getindex(itp, x) = _getindex(itp, convert(Rational, x))
function _getindex(itp, x)
    ix = trunc(x)   # ix = 2//1
    dx = x - ix       # dx = 225179981368525//1125899906842624
    itp.A[ix]*(1-dx) + itp.A[ix+1]*dx
end

So maybe we don't have the user specify the element type: maybe we have them just specify the coefficient type, and we compute the element type (which is pretty easy to do from eltype(A) and the coefficient type).

This way of thinking has some good news: it seems reasonably safe to use Float64 as the default coefficient type, with perhaps a couple of exceptions:

• Float32, in which case we should use Float32.
• Rational, in which case we should use Rational.
  We'd probably like to make the corresponding choices for SIUnits, but until there's a rawtype in base, I don't know of a good generic way of doing that.

[tomasaschan]

Just so I understand you correctly; by element type you refer to essentially eltype(itp), i.e. the type we return from getindex, gradient and friends, while by coefficient type you mean the type we convert our float literals to in the implementation, is that correct?

I'm thinking that if we rely too much on type inference/promotion to get a good element type (by the above definitions), we might risk ending up with everything in Float64 anyway, depending on what e.g. the equation solver does; if the return type is determined entirely by promotion rules, it becomes very difficult (at least for humans) to reason about the return type of itp[x], which might make it difficult to write type-stable code.

An important aspect is that the user should not have to care about any implementation details; the type that is specified should have a clear, and explicit, place in the API. We do require of the input array A that any f::eltype(A) supports operations like f1 * c + f2 * d for some coefficients c and d (whose type we're basically in control of, since they are literals), so if we want the user to specify the type of c and d, it has to be very clear in the documentation that that's what they're doing, and that the return type will be whatever type the result of such a calculation has. It is appealing, though - it would make it really easy to make everythingn Just Work(TM) with an explicit type specification, since we could just let eltype(coefs) == eltype(A) and have even SIUnits types come out correctly, with units.

Regarding sensible defaults, I think it would even be OK - at least as a first step - to default to Float64 always, and be fine with "poisoning" Rational and Float32 inputs, as well as returning incorrect units on SIUnits types; if the user isn't happy with that, it's super easy to specify a type argument when constructing the interpolation, and it never has to be done again after that.

[tomasaschan]

Another thought (that I'll just document here before I forget it): eventually, we'll want something like CoordInterpGrid from Grid.jl. When we implement that, we'll have to make sure that gradient returns the correct unit as well...

[tomasaschan]

Also, won't all this become much easier to implement with staged functions? At least, I imagine we would have to do a lot less of the type inference reasoning ourselves...

[timholy]

If we have float literals in the implementation, then yes, that's right. But we may not even need float literals, as you might see in my implementation of _getindex above. And as long as the user can specify the coefficient type, I think making the default Float64 would be fine.

I suspect we can basically restrict the type reasoning to the following lines in the constructor (we could probably get away with something a little less "rigorous," but...):

Interpolation{Tcoef<:Number}(::Type{Tcoef}, A, ...)
    (isempty(A) || ndims(A) == 0) && error("Cannot construct an Interpolation with an empty or singleton array")
    c = one(Tcoef)
    for i = 2:ndims(A)
        c *= c   # this is the most excessively-rigorous part...
    end
    Tel = typeof(c*A[1] + c*A[1])
    # rest of the constructor
    Interpolation{Tel,N,Tcoef,...}(...)
end

and then in getindex etc methods just use

getindex{Tel,N,Tcoef}(itp::Interpolation{Tel, N, Tcoef}, x) = getindex(itp, convert(Tcoef, x))
function getindex{Tel,N,Tcoef}(itp::Interpolation{Tel, N, Tcoef}, x::Tcoef)
    # here's the real implementation
end

I don't think stagedfunctions will change anything, really. (EDIT: about type-specification, I mean.)

[tomasaschan]

We'll still need literals when specifying the matrix systems for various prefiltering schemes. How do we handle possible type promotion in (pseudo) coefs = M \ part_of_A, when M is specified in terms of literals? Or do we explicitly convert here (e.g. by allocating coefs in some special type)? Is this the type you refer to as Tcoef?

[timholy]

Yes, Tcoef is the coefficient type. I'm proposing we make it a parameter of Interpolation, so we never have to recalculate it.

[tomasaschan]

Thinking another round on it after asking for help in #19 just now, I'm starting to believe converting to TCoefs isn't such a good idea. For "duck type-based multivalued interpolation", we want to allow e.g. TCoefs<:AbstractArray, but indexing with an array doesn't make sense.

I think we need to do some more bikeshedding...

Given the methods defined on AbstractArray, I think it's safe to assume that an index should be a Real. (Cases with other coordinate systems, including e.g. unitful quantities, will have to be handled with a wrapper similar to CoordInterpGrid.) I would also like to make/keep it so that eltype(itp) should be as close to eltype(A) as possible given an interpolation object itp constructed from an array A - thus, if A is Array{Rational{Int}}, then typeof(itp[...]) == Rational{Int} too. (In some cases this doesn't make sense - e.g. if eltype(A) == Int, I think eltype(itp) should be Float64.)

To make it so, we need to be careful about types in the following parts of the code:

• When calculating ix_d based on x_d (e.g. here). This involves somehow rounding the (Real) index x_d.
• When calculating fx_d based on x_d and ix_d (e.g. here). This involves subtraction between ix_d and x_d.
• When calculating the coefficients c_d, cm_d, cp_d etc used in the interpolation expression. For some interpolation types (e.g. linear) this only involves subtraction between 1 and fx_d, but for higher orders (e.g. here for quadratic) it also involves non-one literals and multiplication; this is probably one of the places we have to be the most careful.
• When calculating the interpolated value (e.g. here). This involves multiplication between c_d and an element of type TCoefs, possibly several times.

Given the facts above, it seems that we have to make yet another separation between types, and add another type parameter to the interpolation object; a TIndex<:Real, to which we convert all indices. The contract would then be that typeof(itp[x]) == typeof(one(TIndex)^ndims(A) * one(TCoefs)), i.e. that the return type is the result of promotion between the indexing type and the coefficient type under the operations required to interpolate.

I think we can still get away with having the user specify only one (optional) type argument, namely TIndex. We could default TIndex to Float64, and make special cases for types in base that don't behave well under promotion with Float64, i.e. Rational and Float32 (and possibly Float16, although I've understood that Float16 is just a storage type, so maybe not).

Does this seem like a reasonable approach?

[timholy]

I realize I caused a terminology confusion: I forgot there's a field called coefs that indeed must have multivalued type. What I meant by Tcoef is basically what you're calling TIndex now. I guess I should have called those "weights" or something.

So yes, this is a reasonable approach, with one possible caveat which I may try a PR to fix.

[timholy]

Now I'm beginning to wonder if we want to always use duck-typing, and just preserve whatever the user passes as an index. That would mean we can't make Interpolation a subtype of AbstractArray, however.

[timholy]

Here's the best reason to consider this:

julia> function qinterp1(A, x)
           ix = round(Int, real(x))
           xf = x-ix
           cm = 1//2 * (xf-1//2)^2
           c = 3//4 - xf^2
           cp = 1//2 * (xf+1//2)^2
           cm*A[ix-1] + c*A[ix] + cp*A[ix+1]
       end

julia> A = [(3.8-x)^2+2 for x = 1:10]
10-element Array{Float64,1}:
  9.84
  5.24
  2.64
  2.04
  3.44
  6.84
 12.24
 19.64
 29.04
 40.44

julia> qinterp1(A, 2.1)
5.139999999999999

julia> using DualNumbers

julia> qinterp1(A, dual(2.1,1))
5.139999999999999 - 3.4du

julia> dx = sqrt(eps())
1.4901161193847656e-8

julia> (qinterp1(A, 2.1+dx)-qinterp1(A, 2.1-dx))/(2*dx)  # finite-difference estimate of the gradient
-3.4000000059604645

DualNumbers allows people to check derivatives of any functions that use interpolations internally; it's much higher-precision than finite differencing, and oh-so-much-easier to use. (Related: I notice the gradient tests are very "sloppy" with tolerances of e.g., abs(0.1*f(x)); using DualNumbers for writing tests of this package would be an easy way to test with much higher precision.)

Amplifying on my comment above: since the output type of getindex would no longer be something that can be determined at the time the Interpolation type is constructed, it seems that Interpolation can't be a subtype of AbstractArray.

[tomasaschan]

A compelling argument indeed. What benefits do we have from subtyping AbstractArray, that we'd lose if we didn't?

Another argument to drop Interpolation <: AbstractArray is that this relation makes Interpolation{T,1} <: AbstractVector, but there are a lot of methods taking AbstractVectors which really expect iterables with a finite number of elements, which interpolations sort-of have (we know their size) but you could also argue that they have an infinite number. In other words, I'm not sure that AbstractInterpolation{T,N} actually wants to fulfill the contract that AbstractArray{T,N} wants to define.

[timholy]

The first thing that occurs to me is the loss of the ability to create SubArrays from Interpolation objects. Although perhaps we could loosen the typing on SubArray, making no restrictions on the parent type.

[tomasaschan]

Reading about SubArrays here, it seems like loosening the typing might be risky; to make them work with Interpolations, IIUC they cannot be <:AbstractArray either. I feel that might break a lot of calling code... but if it's workable, that's cool. Supporting SubArray is possibly not the most important aspect of Interpolations either.

[timholy]

You may want to follow any discussion in JuliaLang/julia#9586.

[tomasaschan]

Thanks for the cross-reference; I've subscribed. The results should be interesting =)