Add an efficient iterator for ScaledInterpolation objects

[timholy]

One frequent use for interpolation is to represent a quantity on a "coarse" grid and then construct an array (or perform a computation) in which the object is sampled on a finer grid. In this case, grid points are visited in a predictable sequence, and it turns out we can exploit that predictability to substantially improve the efficiency of interpolation. In the case where the fine grid is much finer than the coarse grid, the cost of the computation goes from O(3^N) (in the case of Quadratic interpolation) to O(3). Especially for larger N, this is a very substantial savings.

The key idea is the following: suppose our coarse grid is 3xMxN, and we're interpolating up to a 201xMxN grid. That means the middle 100 points all correspond, via coordlookup on the coarser grid, to an interval (1.5,2.5), i.e., we have 100 points in a row that have constant ix_1, ix_2, ix_3 values. Since interpolation looks like this:

   cm_1 * pm + c_1 * p + cp_1 * pp

where pm, p, and pp are "partial" results on the trailing dimensions (dimensions 2 and 3), we can compute these partials and then reuse them for all 100 interpolations.

In practice, using the test added here:

# This one uses the iterator
julia> @time foo!(rfoo, sitp);
  0.145194 seconds (9 allocations: 400 bytes)

# This one uses classic sitp[I] evaluation
julia> @time bar!(rfoo, sitp);
  0.581928 seconds (4 allocations: 160 bytes)

so we get an approximate 4-fold speed improvement in 3d (which is huge, although 9-fold would have been even better 😉).

One note: in the help text and test, foo! is written

    function foo!(dest, sitp)
        i = 0
        for s in eachvalue(sitp)
            dest[i+=1] = s
        end
        dest
    end

but this assumes fast linear indexing of dest. Ideally it should be written as

    function foo!(dest, sitp)
        for (I,s) in zip(eachindex(dest), eachvalue(sitp))
            dest[I] = s
        end
        dest
    end

but that turns out to perform poorly due to a type-inference failure. I suspect we're running up against the built-in limits of recursion depth in type inference, i.e, we're getting to the point of having such deep layers of A{B{C{D...}}} that inference.jl is deciding to punt. We should either look into bumping these limits another notch, or we'll have to start being very careful about adding yet more layers of wrapper-types. In particular, while I haven't tested yet I fear that this may behave badly for Extrapolation types simply because they add one extra layer.

[tomasaschan]

This is really nice! :)

Although I realize that we're pushing the limits here, I don't think we're doing anything pathological. Especially, since Julia by design encourages the type of composition that comes from using interpolation, scaling and extrapolation together, I really feel that if this is a problem for type inference, then type inference needs to improve.

[timholy]

Correct, I'm just pointing out a potential challenge. But I'm also no longer sure that this is the correct diagnosis---I tried this on an Array{Vec{3,T},3} (from FixedSizeArrays) and it worked fine, despite the fact that this should add one more layer. I'll poke at it a bit more.

[timholy]

Actually, I was right the first time: bumping MAX_TYPE_DEPTH to 7 fixed it.

[timholy]

JuliaLang/julia#13561

[timholy]

OK to merge this?

[tomasaschan]

Yeah, sure - it just adds functionality without breaking anything, right?

I definitely like the zip(eachindex(dest), eachindex(sitp)) solution for iteration better, but maybe we want to wait with that until the Julia PR merges.

[tomasaschan]

I revisited the scaling code to implement Hessian evaluation today, and realized that the use case for this iterator is really different from what scale is intended to do. I'm not saying the iterators aren't useful - but I'm beginning to feel they should be entirely separated from the scale function and the ScaledInterpolation type. I guess I should have brought this up from the start - apparently I simply wasn't paying close enough attention to the details when this PR was originally submitted. Sorry about that.

My motivation for bringing this up at all is that, except for in this iterator, scale means "take an interpolation object with coordinate axes 1:N1, 1:N2... and scale the axes to linspace(a1,b1,N1), linspace(a2,b2,N2)...." After doing sitp = scale(itp, axes...), you can evaluate sitp[x,y,...] for x ∈ [a1,b1], y ∈ [a1,b2]... instead of x ∈ [1,N1], y ∈ [1,N2]... - in other words, you have changed the endpoints of the axes, but retained their resolution (in terms of number of data points across the entire domain).

When scaling with the iterator, on the other hand, the bounds of the domain are still [1,N1], [1,N2]..., and instead we've just utilized the ordering of operations for interpolating lots of times across a finer grid with the same bounds. This is scaling in the sense that if you still consider the output as indexed with unit-ranges, you've increased its size (i.e. "scaled it up"), but it's not scaling in the sense of changing the domain endpoints, as the scale API was originally intended for.

Perhaps something like rescale or even refine is a better term for it?

[timholy]

Yes, there's nothing very scale-coupled to this, with one caveat: classically, array-iterators in julia index on Integers, and if the grid is "unscaled" then the likelihood is that you'll get one evaluation per grid point, which might perform somewhat worse than "naive" evaluation. If you're scaled and have many integer-grid points for each source-grid point, then this PR can save some time.

However, none of this applies for a general product iterator. It might indeed be worth splitting out. I haven't even thought about the API, however.