WIP: Output type handling (fixes #17)

[tomasaschan]

This is a first attempt at a solution to #17, but it doesn't work quite yet:

julia> using Interpolations

julia> include("test/runtests.jl")
Testing Linear interpolation in 1D...
ERROR: test failed: itp1[-3]
 in error at error.jl:21
 in default_handler at test.jl:19
 in do_test_throws at test.jl:55
 in include at ./boot.jl:245
 in include_from_node1 at ./loading.jl:128
 in include at ./boot.jl:245
 in include_from_node1 at ./loading.jl:128
while loading /home/tlycken/.julia/v0.3/Interpolations/test/linear.jl, in expression starting on line 19
while loading /home/tlycken/.julia/v0.3/Interpolations/test/runtests.jl, in expression starting on line 7

The line in question is

@test_throws BoundsError itp1[-3]

which happens to be the first time the interpolation object is indexed into with an Int - all previous versions calls have indexed with a float...

However, from what I can tell there should be a suitable indexing method defined - I thought I did so here. Also, methods seems a little confused - shouldn't x be an NTuple here?

julia> methods(getindex, (Interpolation{Float64,1,Float64,Linear{OnGrid},ExtrapError},))
1-element Array{Any,1}:
 getindex{T,N,TCoefs<:Real}(itp::Interpolation{T,N,TCoefs<:Real,Linear{OnGrid},ExtrapError},x::TCoefs<:Real...) at cartesian.jl:100

[tomasaschan]

Also, there are a bunch of commits here related to #18, simply because I based this branch off that one. I'll rebase onto master as soon as #18 merges, at which point only cfe3e8a and 345b28f (but under new hashes) will remain in this PR.

[tomasaschan]

Hm. The error message included above is not as informative as the one you get if you add a simple itp1[-3] right above the failing test line:

julia> include("test/linear.jl") # with itp1[-3] on line 18
Testing Linear interpolation in 1D...
ERROR: indexing not defined for Interpolation{Float64,1,Float64,Linear{OnGrid},ExtrapError}
 in getindex at abstractarray.jl:376
 in include at ./boot.jl:245
 in include_from_node1 at ./loading.jl:128
while loading /home/tlycken/.julia/v0.3/Interpolations/test/linear.jl, in expression starting on line 18

I still don't understand why it happens, though.

[tomasaschan]

This now handles Reals nicely, but I have some more work to do in order to make it work with arbitrary element types, which will be needed to interpolate e.g. images where each element is an RGB immutable.

[tomasaschan]

@timholy, When trying to loosen the type constraints on TCoef to allow types that don't inherit Real, I get ambiguity warnings against getindex methods for AbstractArray (e.g. the lowest fallback, and/or linear indexing with a single multidim index). What makes it difficult is that

getindex{T,N,TCoefs}(itp::Interpolation{T,N,TCoefs}, xs::NTuple{N,TCoefs}...)

which is basically paramount to using the approach you suggested in #17, where we convert index arguments to TCoefs, is ambiguous with getindex{T,N}(A::AbstractArray{T,N}, i::Real), at least for N==1 (itp is more specific than A, but xs is less specific than i...).

I've tried to define a bunch of more specific versions, but I keep hitting ambiguity warnings with methods defined for either AbstractArrays or other methods of the ones I define, and I don't really know what's a good strategy to make the constraints on TCoefs as loose as possible.

[tomasaschan]

Almost there! As shown in the tests added in 876b439, we now handle the following very nicely:

julia> begin
           using Interpolations
           A = Float32[sin(x) for x in 1:10]
           itp = Interpolation(A, Quadratic(Free(),OnCell()), ExtrapNaN())
           typeof(itp[3.5]
       end
Float32

When doing similar things with Rational{Int} I keep getting InexactErrors:

julia> using Interpolations
julia> begin
           A = Rational{Int}[x^3/10 for x in 1:10]
           itp = Interpolation(A, Quadratic(Free(),OnCell()), ExtrapNaN())
           typeof(itp[1.1]
       end
Rational{Int64} (constructor with 1 method)

# but
julia> itp[1.2]
ERROR: InexactError()
 in convert at rational.jl:55
 in getindex at /home/tlycken/.julia/v0.3/Interpolations/src/Interpolations.jl:41
 in getindex at no file

(The line number is bogus - I suspect the stack trace is confused by all the macros.) I'm inclined to think that the InexactError is because of limitations in Rationals, not because of something in Interpolations, but I may be wrong.

I also haven't been able to get it working with SIUnits:

julia> using SIUnits

julia> A = rand(10) * Meter
10-element Array{SIQuantity{Float64,m,kg,s,A,K,mol,cd},1}:
  0.5183854583783836 m 
  0.8765064944858227 m 
  0.4192598255113509 m 
  0.9242719465967721 m 
 0.17465399083656008 m 
  0.8607878925114967 m 
      0.569900709139 m 
 0.36160720059717244 m 
 0.38120696026732115 m 
  0.9492328404669232 m 

julia> itpSI = Interpolation(A, Quadratic(Free(), OnCell()), ExtrapNaN())
WARNING: For performance reasons, consider using an array of a concrete type T (eltype(A) == SIQuantity{Float64,m,kg,s,A,K,mol,cd})
ERROR: Unit mismatch. Got () + (m )
 in error at error.jl:21
 in + at /home/tlycken/.julia/v0.3/SIUnits/src/SIUnits.jl:156
 in + at array.jl:721
 in prefilter at /home/tlycken/.julia/v0.3/Interpolations/src/Interpolations.jl:235
 in Interpolation at /home/tlycken/.julia/v0.3/Interpolations/src/Interpolations.jl:73
 in Interpolation at /home/tlycken/.julia/v0.3/Interpolations/src/Interpolations.jl:75

The problem seems to be that although all elements of A are SIQuantity{Float64,1,0,0,0,0,0,0}, A itself is Vector{SIQuantity{Float64,m,kg,s,A,K,mol,cd}}, i.e. eltype(A) is abstract. Thus, when I create the right-hand-side vector based on TCoefs in A::Array{TCoefs,N}, I get the wrong unit. However, if I were to base it on typeof(A[1]) instead, I might get a type that was too tight in case someone gave me e.g. an Array{Number}.

Any suggestions on how to approach these two problems (InexactErrors and arrays of non-leaf types)?

[timholy]

I'm looking into the errors.

But one immediate comment: one problem with insisting on Real coefficients is that then you can't use simply-divine packages like DualNumbers.jl and HyperDualNumbers.jl to numerically check gradients and 2nd derivatives.

[tomasaschan]

Given Jeff's comment jere, I think we can come to the following conclusions:

• It's OK to inherit AbstractArray{T} and still not return T for every possible call to getindex. Good news!
• However, we should return T for indexing with integers. This already works.
• We should also support iteration. Since the default implementation of start, next and done rely on linear indexing, I think we have to roll our own here; this might very well be better handled in a separate issue.

[timholy]

Yep, I agree we're good to go. But on 0.4 we don't need to write any code to implement iteration, since 0.4 supports multidimensional iteration.

[tomasaschan]

@timholy I was hoping you'd say that ;)

I'll see if I can cook up something final for this PR later today.

[tomasaschan]

Huzza! This now works with Float32 (as long as you also index with Float32 - but I think that's a safe assumption...) as well as with Rational (but only on 0.4 due to linalg problems on 0.3; I disabled that test to make Travis pass on both versions) and DualNumbers:

julia> itp[dual(2.0,1.0)] # itp constructed from `rand(10) * 20 - 10`
-9.82487876505064 - 0.24527991346025857du
julia> gradient(itp, 2.0)[1]
-0.24527991346025857

SIUnits don't work yet, but we know why, so I'm not going to bother until the problem is fixed there. Anything else we should test before merging this?

[timholy]

Hooray! I'm all in favor of merging.