Fix performance problems with interpolation

[timholy]

This is built on top of #36; only the last commit is new.

OK, the good news is that the discouraging behavior of #36 (comment) has largely been fixed:

julia> using Interpolations, Grid

julia> A = rand(300,300,30);

julia> itp = interpolate!(copy(A), BSpline(Quadratic(Flat)), OnCell);

julia> yi = InterpGrid(copy(A), BCreflect, InterpQuadratic);

julia> function interp_all!(dest, itp)
           for k = 1:size(itp,3)
               for j = 1:size(itp,2)
                   for i = 1:size(itp,1)
                       dest[i,j,k] = itp[i,j,k]
                   end
               end
           end
           dest
       end
interp_all! (generic function with 1 method)

julia> B = similar(A);

# warmup deleted

julia> @time interp_all!(B, itp);
 214.892 milliseconds (2700 k allocations: 42188 KB)

julia> @time interp_all!(B, yi);
   1.068 seconds      (4 allocations: 144 bytes)

The core problem was that Base's Rational type is no speed demon, because it calls gcd and div repeatedly. So I added a simpler type.

The only remaining puzzle is the memory allocation with itp. It goes away if I change this line to one(typeof(ret)). However, Base.return_types assures me that getindex is type-stable, and @code_warntype did not turn up any oddities, so I am very puzzled.

[timholy]

OK, check this out:

julia> @time interp_all!(B, itp);
  82.396 milliseconds (4 allocations: 144 bytes)

julia> @time interp_all!(B, yi);
   1.086 seconds      (4 allocations: 144 bytes)

[tomasaschan]

Whoa!

I am, as always, immensely impressed by your ability to sqeeze even more out of Julia. This is really beginning to shape up to something like what we (you) envisioned when you started talking about the metaprogramming approach in the issue on Grid.jl.

As I mentioned in the other comment I unfortunately don't have much time to help you out at the moment, but I trust you to write Good Stuff, so feel free to merge at will while I'm "away" :)

[tomasaschan]

Without having perused the code carefully, this seems like something we shouldn't have to do; if for no other reason, then because there is probably also other packages and code bases that could benefit from these speed improvements.

In the long term, maybe we should split these out and incorporate them either in Base or in their own package? (I imagine these do a subset of the stuff that Base.Rational do - maybe we could re-write Rationals to be a two-step rocket, exposing these too? Loong-term goal...)

[timholy]

Presumably we can't make Rational itself behave this way, because if you don't call gcd and div, then overflow seems likely. But you're right that perhaps we could use a layered approach.

CCing @StefanKarpinski, @simonbyrne, and @JeffBezanson. Executive summary: Interpolations is the successor to Grid, with an aim to improve the API and performance. At the top of this issue I linked to a performance problem that stemmed from using Rational, where I got a 50x speed boost by rolling my own stripped-down version (called Ratio). We're attracted to use something like Rational for coefficients like 1//2, rather than floating-point numbers like 0.5, because it works out better in terms of promotion behavior for different Number types. I don't immediately see a good way of designing Rational to get both behaviors, but wondering if you folks have any thoughts about this.

[StefanKarpinski]

It would be worth opening an issue on the main julia repo to track this issue – calling gcd and lcm all the time is pretty awful for performance. There's also the chronic overflow issue that Rational has, which seems like something that would need to be address to use the for any kind of real computation.

[simonbyrne]

I too would like to see a faster Rational type, though overflow should have been fixed since JuliaLang/julia#8672.