RFC: Scaling of interpolation objects (fixes #25) #47
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,5 +1,5 @@ | ||
| export Throw, | ||
| FilledExtrapolation # for direct control over typeof(fillvalue) | ||
|
|
||
| include("error.jl") | ||
|
|
||
|
|
||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,36 +1,39 @@ | ||
| nindexes(N::Int) = N == 1 ? "1 index" : "$N indexes" | ||
|
|
||
|
|
||
| type FilledExtrapolation{T,N,ITP<:AbstractInterpolation,IT,GT,FT} <: AbstractExtrapolation{T,N,ITP,IT,GT} | ||
| itp::ITP | ||
| fillvalue::FT | ||
| end | ||
| """ | ||
| `FilledExtrapolation(itp, fillvalue)` creates an extrapolation object that returns the `fillvalue` any time the indexes in `itp[x1,x2,...]` are out-of-bounds. | ||
|
|
||
| By comparison with `extrapolate`, this version lets you control the `fillvalue`'s type directly. It's important for the `fillvalue` to be of the same type as returned by `itp[x1,x2,...]` for in-bounds regions for the index types you are using; otherwise, indexing will be type-unstable (and slow). | ||
| """ | ||
| function FilledExtrapolation{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue) | ||
| FilledExtrapolation{T,N,typeof(itp),IT,GT,typeof(fillvalue)}(itp, fillvalue) | ||
| end | ||
|
|
||
| """ | ||
| `extrapolate(itp, fillvalue)` creates an extrapolation object that returns the `fillvalue` any time the indexes in `itp[x1,x2,...]` are out-of-bounds. | ||
| """ | ||
| extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue) = FilledExtrapolation(itp, convert(eltype(itp), fillvalue)) | ||
|
|
||
| @generated function getindex{T,N}(fitp::FilledExtrapolation{T,N}, args::Number...) | ||
| n = length(args) | ||
| n == N || return error("Must index $(N)-dimensional interpolation objects with $(nindexes(N))") | ||
| meta = Expr(:meta, :inline) | ||
| quote | ||
| $meta | ||
| # Check to see if we're in the extrapolation region, i.e., | ||
| # out-of-bounds in an index | ||
| @nexprs $N d->((args[d] < lbound(fitp,d) || args[d] > ubound(fitp, d)) && return fitp.fillvalue) | ||
| # In the interpolation region | ||
| return getindex(fitp.itp,args...) | ||
| end | ||
| end | ||
|
|
||
| getindex{T}(fitp::FilledExtrapolation{T,1}, x::Number, y::Int) = y == 1 ? fitp[x] : throw(BoundsError()) | ||
|
|
||
| lbound(etp::FilledExtrapolation, d) = lbound(etp.itp, d) | ||
| ubound(etp::FilledExtrapolation, d) = ubound(etp.itp, d) |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,13 +1,13 @@ | ||
| @generated function getindex{T}(etp::AbstractExtrapolation{T,1}, x) | ||
| quote | ||
| $(extrap_prep(etp, x)) | ||
| etp.itp[x] | ||
| end | ||
| end | ||
|
|
||
| @generated function getindex{T,N,ITP,GT}(etp::AbstractExtrapolation{T,N,ITP,GT}, xs...) | ||
| quote | ||
| $(extrap_prep(etp, xs...)) | ||
| etp.itp[xs...] | ||
| end | ||
| end |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,90 @@ | ||
| export ScaledInterpolation | ||
|
|
||
| type ScaledInterpolation{T,N,ITPT,IT,GT,RT} <: AbstractInterpolationWrapper{T,N,ITPT,IT,GT} | ||
| itp::ITPT | ||
| ranges::RT | ||
| end | ||
| 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) | ||
| """ | ||
| `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. | ||
|
Member
There was a problem hiding this comment. In a sense, Given that conceptual equivalence, however, I suspect we should not have
Contributor
Author
There was a problem hiding this comment. I agree that the conceptual similarity is a good argument for similarity in API as well. However, for the same reason, but looking at tl;dr: I actually think it's better to change
Member
There was a problem hiding this comment. Fine with me. I chose the first mostly because Matlab does it that way.
Member
There was a problem hiding this comment. Although now that I think about it, for Also note scipy.
Contributor
Author
There was a problem hiding this comment. Yes - I agree that In the context of composition of various "add-ons", though, I think keeping the API as close to |
||
|
|
||
| The parameters `xs` etc must be either ranges or linspaces, and there must be one coordinate range/linspace for each dimension of the interpolation object. | ||
|
|
||
| For every `NoInterp` dimension of the interpolation object, the range must be exactly `1:size(itp, d)`. | ||
|
Contributor
Author
There was a problem hiding this comment. For
Member
There was a problem hiding this comment. Doesn't have to be part of this PR (or happen ever). But it's fine even for interpolating dimensions; after all, you let me supply |
||
| """ | ||
| function scale{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, ranges::Range...) | ||
| length(ranges) == N || throw(ArgumentError("Must scale $N-dimensional interpolation object with exactly $N ranges (you used $(length(ranges)))")) | ||
| for d in 1:N | ||
| if iextract(IT,d) != NoInterp | ||
| length(ranges[d]) == size(itp,d) || throw(ArgumentError("The length of the range in dimension $d ($(length(ranges[d]))) did not equal the size of the interpolation object in that direction ($(size(itp,d)))")) | ||
| elseif ranges[d] != 1:size(itp,d) | ||
| throw(ArgumentError("NoInterp dimension $d must be scaled with unit range 1:$(size(itp,d))")) | ||
| end | ||
| end | ||
|
|
||
| ScaledInterpolation(T,N,itp,IT,GT,ranges) | ||
| end | ||
|
|
||
| @generated function getindex{T,N,ITPT,IT<:DimSpec}(sitp::ScaledInterpolation{T,N,ITPT,IT}, xs::Number...) | ||
| length(xs) == N || throw(ArgumentError("Must index into $N-dimensional scaled interpolation object with exactly $N indices (you used $(length(xs)))")) | ||
| interp_types = length(IT.parameters) == N ? IT.parameters : tuple([IT.parameters[1] for _ in 1:N]...) | ||
| interp_dimens = map(it -> interp_types[it] != NoInterp, 1:N) | ||
| interp_indices = map(i -> interp_dimens[i] ? :(coordlookup(sitp.ranges[$i], xs[$i])) : :(xs[$i]), 1:N) | ||
| return :(getindex(sitp.itp, $(interp_indices...))) | ||
| end | ||
|
|
||
| getindex{T}(sitp::ScaledInterpolation{T,1}, x::Number, y::Int) = y == 1 ? sitp[x] : throw(BoundsError()) | ||
|
|
||
| size(sitp::ScaledInterpolation, d) = size(sitp.itp, d) | ||
| lbound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnGrid}, d) = 1 <= d <= N ? sitp.ranges[d][1] : throw(BoundsError()) | ||
| lbound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnCell}, d) = 1 <= d <= N ? sitp.ranges[d][1] - boundstep(sitp.ranges[d]) : throw(BoundsError()) | ||
| ubound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnGrid}, d) = 1 <= d <= N ? sitp.ranges[d][end] : throw(BoundsError()) | ||
| ubound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnCell}, d) = 1 <= d <= N ? sitp.ranges[d][end] + boundstep(sitp.ranges[d]) : throw(BoundsError()) | ||
|
|
||
| boundstep(r::LinSpace) = ((r.stop - r.start) / r.divisor) / 2 | ||
| boundstep(r::FloatRange) = r.step / 2 | ||
| boundstep(r::StepRange) = r.step / 2 | ||
| boundstep(r::UnitRange) = 1//2 | ||
|
|
||
| """ | ||
| Returns *half* the width of one step of the range. | ||
|
|
||
| This function is used to calculate the upper and lower bounds of `OnCell` interpolation objects. | ||
| """ boundstep | ||
|
|
||
| coordlookup(r::LinSpace, x) = (r.divisor * x + r.stop - r.len * r.start) / (r.stop - r.start) | ||
| coordlookup(r::FloatRange, x) = (r.divisor * x - r.start) / r.step + one(eltype(r)) | ||
| coordlookup(r::StepRange, x) = (x - r.start) / r.step + one(eltype(r)) | ||
| coordlookup(r::UnitRange, x) = x - r.start + one(eltype(r)) | ||
| coordlookup(i::Bool, r::Range, x) = i ? coordlookup(r, x) : convert(typeof(coordlookup(r,x)), x) | ||
|
|
||
| gradient{T,N,ITPT,IT<:DimSpec}(sitp::ScaledInterpolation{T,N,ITPT,IT}, xs::Number...) = gradient!(Array(T,count_interp_dims(IT,N)), sitp, xs...) | ||
| @generated function gradient!{T,N,ITPT,IT}(g, sitp::ScaledInterpolation{T,N,ITPT,IT}, xs::Number...) | ||
| ndims(g) == 1 || throw(DimensionMismatch("g must be a vector (but had $(ndims(g)) dimensions)")) | ||
| length(xs) == N || throw(DimensionMismatch("Must index into $N-dimensional scaled interpolation object with exactly $N indices (you used $(length(xs)))")) | ||
|
|
||
| interp_types = length(IT.parameters) == N ? IT.parameters : tuple([IT.parameters[1] for _ in 1:N]...) | ||
| interp_dimens = map(it -> interp_types[it] != NoInterp, 1:N) | ||
| interp_indices = map(i -> interp_dimens[i] ? :(coordlookup(sitp.ranges[$i], xs[$i])) : :(xs[$i]), 1:N) | ||
|
|
||
| quote | ||
| length(g) == $(count_interp_dims(IT, N)) || throw(ArgumentError(string("The length of the provided gradient vector (", length(g), ") did not match the number of interpolating dimensions (", $(count_interp_dims(IT, N)), ")"))) | ||
| gradient!(g, sitp.itp, $(interp_indices...)) | ||
| for i in eachindex(g) | ||
| g[i] = rescale_gradient(sitp.ranges[i], g[i]) | ||
| end | ||
| g | ||
| end | ||
| end | ||
|
|
||
| rescale_gradient(r::LinSpace, g) = g * r.divisor / (r.stop - r.start) | ||
| rescale_gradient(r::FloatRange, g) = g * r.divisor / r.step | ||
| rescale_gradient(r::StepRange, g) = g / r.step | ||
| rescale_gradient(r::UnitRange, g) = g | ||
|
|
||
| """ | ||
| `rescale_gradient(r::Range)` | ||
|
|
||
| Implements the chain rule dy/dx = dy/du * du/dx for use when calculating gradients with scaled interpolation objects. | ||
| """ rescale_gradient | ||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe also export
ScaledInterpolation?There was a problem hiding this comment.
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
scalecan't do (unlikeFilledExtrapolation, 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.There was a problem hiding this comment.
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 😄):
But of course I can just add the module name, if you don't want to export it.
There was a problem hiding this comment.
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 :)
There was a problem hiding this comment.
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
knotsthat applies to anyitp.Also, now that I think about it,
scaled(knots, ::GriddedInterpolant)should probably return anotherGriddedInterpolant, one that shares the data array but has new knots.There was a problem hiding this comment.
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 aScaledInterpolantin most cases, but wherescale(itp::GriddedInterpolation, knots...)can special-case for re-use of the data array.