RFC: Scaling of interpolation objects (fixes #25)

[tomasaschan]

This is a start at implementing the functionality provided to Grid.jl by the CoordInterpGrid type, and in fact the code here is heavily based on that contributed by @simonbyrne there.

Simon, do you have any objections to me basing this implementation off of your code? If so, please feel free to object here, and I'll remove this and come up with something from scratch. It's just that borrowing is so easy ;)

The basic idea is that you can do

f(x) = exp(-x.^2 ./ 4)
xs = -5:.2:5
ys = f(xs)
itp = interpolate(ys, BSpline(Linear), OnGrid)

s = scale(itp, xs)   # now, s[-2.3] ~= f(-2.3)

To scale in higher dimensions, you just tuck on more ranges at the end, one range for each dimension (i.e. scale(itp2, xs, ys) to scale a 2D interpolation object).

This doesn't fully handle everything we want yet - for one thing, there's no scaling of the coordinates that works for gradients. Also, in Grid.jl @timholy added specializations for dimensions 1-4 that avoided splatting. I'd like to see if it's possible to do this using staged functions instead of explicitly writing them out, but this might not be fully performant yet. It also needs documentation.

However, if anyone has any opinions about the API, this is as good a time as any to raise them.

Todo:

• basic implementation
• gradient evaluation
• make sure it works with DimSpecs and NoInterp dimensions
• tests for combinations with extrapolate (both extrapolate(scale(...)) and scale(extrapolate(...)))
• make sure scaled interpolation objects display nicely
• documentation

[ghost]

Sounds good! Presumably it will be possible to make this one call, e.g. interpolate((x1s, x2s), ys, BSpline(Linear), OnGrid) or maybe sinterpolate to prevent confusion?

[tomasaschan]

@gaabriel, actually, I'm inclined towards not implementing that method, for the following reason:

In the long run, I'm hoping to be able to include more types of interpolations than B-splines in this library. It would be nice if things like re-scaling of coordinate axes (at least for uniform grids), extrapolation and maybe other functionality could share code between the different types of interpolation. Currently, this is handled through composition, so for example the following is possible:

itp = interpolate(ys, BSpline(Linear), OnGrid)
extp = extrapolate(itp, Line) # not implemented yet, but it will be :)
sctp = scale(extp, x1s, x2s)

Now, sctp interpolates ys linearly, extrapolates linearly for OOB coordinates, and has its axes re-scaled to x1s and x2s. If extrapolation (or scaling) is unimportant for a use case, just leave that row out.

If we were to include all possible combinations of compositions as methods on interpolate, the number of methods would explode, and I think it would be much more difficult to maintain.

Of course, if doing it all in one line is important to you, there's nothing preventing you from passing the interpolation object directly to the wrapper:

sctp = scale(extrapolate(interpolate(ys, BSpline(Linear), OnGrid), Line), x1s, x2s)

although admittedly this doesn't read as nicely.

[simonbyrne]

Looks great, no objections from me!

[timholy]

Here and elsewhere, you probably want one(eltype(r)) rather than 1.0 (think Float32).

[timholy]

Come to think of it, it wouldn't be bad to add some Float32 tests...

[tomasaschan]

Yeah, tests for Float32 would be nice. Maybe Rationals, Complex and some other non-trivial types too...

Would e.g. 1//1 do here? Promotion should take care of the rest then, I imagine.

[timholy]

Probably slower than using one directly, I assume. Does that have any failure cases?

[tomasaschan]

True; promotion costs. one is probably better then; in this case, it's unambiguous what type is the right one, so it'll be difficult to get it wrong. (I was making parallels to the discussions we had about literal constants elsewhere in the library a while ago...)

[tomasaschan]

@timholy I've started working on implementing support for NoInterp, and I realized that it would be very useful to have NoInterp be independent of all the B-spline stuff, in order to keep scale and extrapolate agnostic of what interpolation scheme the underlying interpolation object is using.

Currently, the definition is as follows:

immutable NoInterp <: Degree(0) end

but I'd much rather have

immutable NoInterp <: InterpolationType

i.e. to put it on equal level with BSpline{D<:Degree} in the type hierarchy. However, that will make it not be a subtype of DimSpec{BSpline}, which breaks all sorts of things (because a lot of the type parameter specifications on the internal methods have IT<:DimSpec{BSpline} on them, and so NoInterp won't meet that requirement anymore), and I'm not quite sure how to resolve that properly. If I just remove the requirement for a subtype of BSpline (i.e. IT<:DimSpec) non-B-spline interpolation types will match. If I don't, NoInterp won't match.

Do you have any suggestions on how best to move forward?

[tomasaschan]

(Re: the current test failure; that's a @test_approx_eq a b that should probably have been @test_approx_eq_eps a b cbrt(eps(a)) - I'll fix that in a later commit.)

[tomasaschan]

I'm seeing a type instability in (some of) these methods that I can't pin down and explain (reproduction script). It seems that Julia is unable to infer the type even of coordlookup(2:2.3:10, 5.1), which is super-weird, since the quantities involved in calculations in that case are Float64.

julia> @code_warntype(Interpolations.coordlookup(2:2.3:10, 4.1))
Variables:
  r::FloatRange{Float64}
  x::Float64

Body:
  begin  # C:\Users\Tomas Lycken\.julia\v0.4\Interpolations\src\scaling/scaling.jl, line 33:
      return (Base.box)(Base.Float64,(Base.add_float)((Base.box)(Base.Float64,(Base.div_float)((Base.box)(Base.Float64,(
Base.sub_float)((Base.box)(Base.Float64,(Base.mul_float)((top(getfield))(r::FloatRange{Float64},:divisor)::Float64,x::Fl
oat64)::Any)::Float64,(top(getfield))(r::FloatRange{Float64},:start)::Float64)::Any)::Float64,(top(getfield))(r::FloatRa
nge{Float64},:step)::Float64)::Any)::Float64,(Base.box)(Float64,(Base.sitofp)(Float64,1)::Any)::Float64)::Any)::Float64
  end::Float64

Note that there's a bunch of boxing operations as well as Any inferred.

[KristofferC]

The boxes and Anys are not a problem. For example

julia> @code_warntype sqrt(5.0)
Variables:
  x::Float64

Body:
  begin  # math.jl, line 144:
      return (Base.Math.box)(Base.Math.Float64,(Base.Math.sqrt_llvm)(x::Float64)::Any)::Float64
  end::Float64

What is important is the output which for coordlookup is just Float64. I would instead focus on the map call in getindex and the splatting of the result. These can be problematic for type inference.

[timholy]

You have a run-time anonymous function:

    # handle dimspecs
    interp_indices = map(it -> IT.parameters[it] != NoInterp, 1:length(IT.parameters))

Those are always type-unstable. You could use FastAnonymous to make it type-stable, but it's also always unstable to access the .parameters field.

Since IT.parameters is available to compile-time portion of the @generated function, you can (and should) do all this processing before returning the Expr.

[tomasaschan]

@timholy That was intended to be compile-time (note that the quoting is only around the actual getindex call, where interp_indices have already been calculated and are interpolated into the returned expression).

And indeed; if I add a @show interp_indices after the calculation there, it only outputs something on the first execution of getindex (and it seems to know what type it is, too).

Edit: ah, now I see it. The map call you indicated isn't problematic, but the ones in the methods of coordlookup that takes tuples are. I'll have to keep working on this, then...

[timholy]

Shoot, I was mixing two things up. Glad you figured it out.

For performance reasons, you probably also don't want ... in returned expressions: you probably want to build the exprs like this:

lookup_indices = [:(coord_lookup(sitp.ranges, index)) for index in interp_indices]
:(getindex(sitp.itp, $(lookup_indices...)))

[tomasaschan]

That seems to have solved it. I haven't done any benchmarks, but if we gained a few extra notches on the performance scale thanks to your tips I'm not complaining :)

[timholy]

Great to hear. This looks nice, BTW!

[tomasaschan]

I think this is all done now - would gladly have someone look at it and see if there's something more that should be done (e.g. if something more needs documentation - I've only documented the exported API and a couple of the internal helper functions whose purpose I thought wasn't necessarily obvious from the name, but I haven't documented everything...)

[timholy]

In a sense, (xs,ys,...) is equivalent to knots in the Gridded interpolant, with the very important difference that Range inputs provide very nice extra algorithmic guarantees that make everything more efficient and mesh seamlessly with our BSpline infrastructure.

Given that conceptual equivalence, however, I suspect we should not have knots come first with Gridded and xs, ys, ... come second with scale. I would propose scale(knots, itp) but I am comfortable changing Gridded instead.

[tomasaschan]

I agree that the conceptual similarity is a good argument for similarity in API as well. However, for the same reason, but looking at extrapolate instead, I'd like to keep scale(itp, ranges...) this way, since both scale and extrapolate essentially say "take this interpolant (the first argument) and add features to it according to the rest of the parameters".

tl;dr: I actually think it's better to change Gridded instead, even though having knots as the first argument is probably more intuitive in isolation.

[timholy]

Fine with me. I chose the first mostly because Matlab does it that way.

[timholy]

Although now that I think about it, for Gridded it's not adding features; it's a different algorithm. For arbitrary knot vectors you have to use searchsorted to find the right spot in the data array.

Also note scipy.

[tomasaschan]

Yes - I agree that interpolate(knots, data, Gridded) is a better interface, as is probably scale(knots, itp) if that was the only thing we did (and the fact that both MATLAB, SciPy and, I assume, a whole bunch of others do it that way is fuel on the fire).

In the context of composition of various "add-ons", though, I think keeping the API as close to do_extra_stuff(itp, args...) has a lot of value. scale and extrapolate are both examples of do_extra_stuff in this context. If we don't adhere to (itp, args...) for all of them, it will be a lot more difficult for users of this library to know (without resorting to docs) in which order the arguments should come.

[timholy]

Very nice! Most importantly, I'm already looking forward to using this! Thanks for tackling it.

[timholy]

Maybe also export ScaledInterpolation?

[tomasaschan]

I don't think there's need for that; there's nothing you can do with that constructor that scale can't do (unlike FilledExtrapolation, where there's an actual use case for explicitly using the constructor).

We should, though, probably export AbstractInterpolation, to make it easier to create functions that use these objects. But the whole point of the composability design is that you shouldn't have to know which the concrete type is.

[timholy]

5 minutes ago I wrote this method (see, I'm already using it 😄):

function GridDeformation{FV<:FixedVector}(u::ScaledInterpolation{FV})
    N = length(FV)
    ndims(u) == N || throw(DimensionMismatch("Dimension $(ndims(u)) incompatible with vectors of length $N"))
    GridDeformation{eltype(FV),N,typeof(u),typeof(u.ranges[1])}(u)
end

But of course I can just add the module name, if you don't want to export it.

[tomasaschan]

Nvm, then - I stand corrected :)

[timholy]

Or maybe better is a new function knots that applies to any itp.

Also, now that I think about it, scaled(knots, ::GriddedInterpolant) should probably return another GriddedInterpolant, one that shares the data array but has new knots.

[tomasaschan]

In that case, I'd argue for scale(itp, knots...) as the common interface, returning a ScaledInterpolant in most cases, but where scale(itp::GriddedInterpolation, knots...) can special-case for re-use of the data array.

[tomasaschan]

There, I think I addressed all the comments (at least where I'm not of the opinion that nothing should change - you may still change my mind... ;)).

The things I haven't done anything about (and where there are thus still line comments open):

• Using a bi-quadratic function for testing of 2D scaling
• Colon indexing
• Argument ordering for scale(itp, knots...)

The two first ones I think can wait until later, and the third I'm pretty sure I want as it is.

Any further comments, or is it hit-the-button-time?

[timholy]

Just posted a couple of other comments above.

To be clear, while I advocate the knots first, in the end I don't really care. Hit the button at will 😄.

[timholy]

Hooray! Again, thanks for another fantastic addition.

[tomasaschan]

🎉 🎈 🍻

Next up: #50... ;)

[timholy]

This package quite desperately needs tagging. I am happy to do it, but the question is: "patch" or "minor"? I think the only breaking change is the FilledInterpolation->FilledExtrapolation, but neither type was available in the previously-tagged version. So perhaps "patch" is sufficient?

(Alternatively, we could use the fancy new @deprecate_binding, but that will introduce a dependency on julia-0.4-rc2.)

[tomasaschan]

I think minor could be justified anyway, given all the new functionality we're bringing in. I'll tag when I'm done with #50, but feel free to do it yourself if you don't want to wait.

I don't see the need for deprecation of features we've never actually released :)