Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

RFC: Scaling of interpolation objects (fixes #25) #47

Merged
merged 8 commits into from
Sep 21, 2015
Merged

Conversation

tomasaschan
Copy link
Contributor

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

This was referenced Jul 20, 2015
@ghost
Copy link

ghost commented Jul 21, 2015

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
Copy link
Contributor Author

@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
Copy link
Member

Looks great, no objections from me!


size(sitp::ScaledInterpolation, d) = size(sitp.itp, d)

coordlookup(r::FloatRange, x) = (r.divisor * x - r.start) ./ r.step + 1.0
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Contributor Author

@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
Copy link
Contributor Author

(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.)

coordlookup(r::UnitRange, x) = x - r.start + one(eltype(r))
coordlookup(i::Bool, r::Range, x) = i ? coordlookup(r, x) : x
coordlookup{N}(r::NTuple{N}, x::NTuple{N}) = map(coordlookup,r,x)
coordlookup{N}(i::NTuple{N}, r::NTuple{N}, x::NTuple{N}) = map(coordlookup, i, r, x)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Member

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
Copy link
Member

timholy commented Sep 18, 2015

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
Copy link
Contributor Author

@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
Copy link
Member

timholy commented Sep 18, 2015

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
Copy link
Contributor Author

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
Copy link
Member

timholy commented Sep 18, 2015

Great to hear. This looks nice, BTW!

tomasaschan added a commit that referenced this pull request Sep 20, 2015
This change might have been better suited for its own PR, but I stumbled
over it now and had to change some related things anyway (indexing into
it, see next commit in #47).

@timholy, please protest if you have any objections.
This change might have been better suited for its own PR, but I stumbled
over it now and had to change some related things anyway (indexing into
it, see next commit in #47).

@timholy, please protest if you have any objections.
@tomasaschan tomasaschan changed the title WIP: Scaling of interpolation objects (fixes #25) RFC: Scaling of interpolation objects (fixes #25) Sep 21, 2015
@tomasaschan
Copy link
Contributor Author

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...)

ScaledInterpolation{T,ITPT,IT,GT,RT}(::Type{T}, N, itp::ITPT, ::Type{IT}, ::Type{GT}, ranges::RT) =
ScaledInterpolation{T,N,ITPT,IT,GT,RT}(itp, ranges)
@doc """
`scale(itp, xs, ys, ...)` scales an existing interpolation object to allow for indexing using other coordinate axes than unit ranges, by wrapping the interpolation object and transforming the indices from the provided axes onto unit ranges upon indexing.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Member

timholy commented Sep 21, 2015

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

@@ -4,6 +4,7 @@ export
interpolate,
interpolate!,
extrapolate,
scale,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe also export ScaledInterpolation?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nvm, then - I stand corrected :)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Contributor Author

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
Copy link
Member

timholy commented Sep 21, 2015

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 😄.

tomasaschan pushed a commit that referenced this pull request Sep 21, 2015
RFC: Scaling of interpolation objects (fixes #25)
@tomasaschan tomasaschan merged commit 7646fe9 into master Sep 21, 2015
@timholy
Copy link
Member

timholy commented Sep 21, 2015

Hooray! Again, thanks for another fantastic addition.

@tomasaschan
Copy link
Contributor Author

🎉 🎈 🍻

Next up: #50... ;)

@tomasaschan tomasaschan deleted the scaling branch September 21, 2015 12:29
@timholy
Copy link
Member

timholy commented Sep 22, 2015

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
Copy link
Contributor Author

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 :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants