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.

Sign up for GitHub

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

WIP: convert to stagedfunctions #26

Closed
wants to merge 6 commits into from
Closed

WIP: convert to stagedfunctions #26

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
8fad89c
Fix concatenation deprecations
timholy Feb 17, 2015
6cea519
WIP: convert to stagedfunctions from ngenerate
timholy Feb 17, 2015
2fdd7a9
Improve quality of dimensionality error messages
timholy Feb 17, 2015
c1404c9
More concat deprecations
timholy Feb 17, 2015
3b028cd
Fix strange stagedfunction bug
timholy Feb 17, 2015
f1bdb33
Travis: test only on nightly
timholy Feb 17, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions REQUIRE
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1 +1,3 @@
julia 0.4-
Compat
WoodburyMatrices
215 changes: 105 additions & 110 deletions src/Interpolations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
module Interpolations

using Base.Cartesian
using Compat
using WoodburyMatrices, Compat

import Base:
eltype,
Expand All @@ -10,7 +10,7 @@ import Base:
ndims,
size

export
export
Interpolation,
Constant,
Linear,
Expand Down Expand Up @@ -123,131 +123,126 @@ function copy_with_padding(A, it::InterpolationType)
coefs, pad
end

# This creates getindex methods for all supported combinations
for IT in (
Constant{OnGrid},
Constant{OnCell},
Linear{OnGrid},
Linear{OnCell},
Quadratic{Flat,OnCell},
Quadratic{Flat,OnGrid},
Quadratic{Reflect,OnCell},
Quadratic{Reflect,OnGrid},
Quadratic{Line,OnGrid},
Quadratic{Line,OnCell},
Quadratic{Free,OnGrid},
Quadratic{Free,OnCell},
Quadratic{Periodic,OnGrid},
Quadratic{Periodic,OnCell},
)
for EB in (
ExtrapError,
ExtrapNaN,
ExtrapConstant,
ExtrapLinear,
ExtrapReflect,
ExtrapPeriodic,
)
it = IT()
eb = EB()
gr = gridrepresentation(it)

function getindex_impl(N)
quote
# Handle extrapolation, by either throwing an error,
# returning a value, or shifting x to somewhere inside
# the domain
$(extrap_transform_x(gr,eb,N))

# Calculate the indices of all coefficients that will be used
# and define fx = x - xi in each dimension
$(define_indices(it, N))

# Calculate coefficient weights based on fx
$(coefficients(it, N))

# Generate the indexing expression
@inbounds ret = $(index_gen(degree(it), N))
ret
end
end
supported(::Type{Constant{OnGrid}}) = true
Copy link
Contributor

Choose a reason for hiding this comment

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

Is supported used anywhere by us, or is it just a type trait that can be used elsewhere?

Copy link
Member Author

Choose a reason for hiding this comment

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

No uses yet. I was keying off the # This creates getindex methods for all supported combinations and wondered if you wanted to be able to return an error like "not a supported combination" from something. So since the list was handy, I turned it into a trait. If we don't end up using it at all, then these should be deleted.

supported(::Type{Constant{OnCell}}) = true
supported(::Type{Linear{OnGrid}}) = true
supported(::Type{Linear{OnCell}}) = true
supported(::Type{Quadratic{Flat,OnCell}}) = true
supported(::Type{Quadratic{Flat,OnGrid}}) = true
supported(::Type{Quadratic{Reflect,OnCell}}) = true
supported(::Type{Quadratic{Reflect,OnGrid}}) = true
supported(::Type{Quadratic{Line,OnGrid}}) = true
supported(::Type{Quadratic{Line,OnCell}}) = true
supported(::Type{Quadratic{Free,OnGrid}}) = true
supported(::Type{Quadratic{Free,OnCell}}) = true
supported(::Type{Quadratic{Periodic,OnGrid}}) = true
supported(::Type{Quadratic{Periodic,OnCell}}) = true
supported(::Any) = false

function getindex_impl(N, IT::Type, EB::Type)
it = IT()
eb = EB()
gr = gridrepresentation(it)
quote
@nexprs $N d->(x_d = x[d])

# Handle extrapolation, by either throwing an error,
# returning a value, or shifting x to somewhere inside
# the domain
$(extrap_transform_x(gr,eb,N))

# Calculate the indices of all coefficients that will be used
# and define fx = x - xi in each dimension
$(define_indices(it, N))

# Calculate coefficient weights based on fx
$(coefficients(it, N))

# Generate the indexing expression
@inbounds ret = $(index_gen(degree(it), N))
ret
end
end

# Resolve ambiguity with linear indexing
getindex{T,N,TCoefs,IT,EB}(itp::Interpolation{T,N,TCoefs,IT,EB}, x::Real) =
error("Linear indexing is not supported for interpolation objects")
stagedfunction getindex{T,TCoefs,IT,EB}(itp::Interpolation{T,1,TCoefs,IT,EB}, x::Real)
getindex_impl(1, IT, EB)
end

# Resolve ambiguity with real multilinear indexing
stagedfunction getindex{T,N,TCoefs,IT,EB}(itp::Interpolation{T,N,TCoefs,IT,EB}, x::Real...)
n = length(x)
n == N || return :(error("Must index interpolation objects with as many indexes as dimensions"))
getindex_impl(N, IT, EB)
end

# Resolve ambiguity with linear indexing,
# getindex{T,N}(A::AbstractArray{T,N}, i::Real)
eval(ngenerate(:N, :T, :(getindex{T,N,TCoefs}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::Real)),
N->:(error("Linear indexing is not supported for interpolation objects"))
))

# Resolve ambiguity with real multilinear indexing
# getindex{T,N}(A::AbstractArray{T,N}, NTuple{N,Real}...)
eval(ngenerate(:N, :T, :(getindex{T,N,TCoefs}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::NTuple{N,Real}...)), getindex_impl))

# Resolve ambiguity with general array indexing,
# getindex{T,N}(A::AbstractArray{T,N}, I::AbstractArray{T,N})
eval(ngenerate(:N, :T, :(getindex{T,N,TCoefs}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, I::AbstractArray{T})),
N->:(error("Array indexing is not defined for interpolation objects."))
))

# Resolve ambiguity with colon indexing for 1D interpolations
# getindex{T}(A::AbstractArray{T,1}, C::Colon)
eval(ngenerate(:N, :T, :(getindex{T,TCoefs}(itp::Interpolation{T,1,TCoefs,$IT,$EB}, c::Colon)),
N->:(error("Colon indexing is not supported for interpolation objects"))
))

# Allow multilinear indexing with any types
eval(ngenerate(:N, :(promote_type(T,TIndex)), :(getindex{T,N,TCoefs,TIndex}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::NTuple{N,TIndex}...)), getindex_impl))

# Define in-place gradient calculation
# gradient!(g::Vector, itp::Interpolation, NTuple{N}...)
eval(ngenerate(:N, :(Vector{TOut}), :(gradient!{T,N,TCoefs,TOut}(g::Vector{TOut}, itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::NTuple{N,Any}...)),
N->quote
length(g) == $N || error("g must be an array with exactly N elements (length(g) == "*string(length(g))*", N == "*string(N)*")")
$(extrap_transform_x(gr,eb,N))
$(define_indices(it,N))
@nexprs $N dim->begin
@nexprs $N d->begin
(d==dim
? $(gradient_coefficients(it,N,:d))
: $(coefficients(it,N,:d)))
end

@inbounds g[dim] = $(index_gen(degree(it),N))
end
g
# Resolve ambiguity with general array indexing,
# getindex{T,N}(A::AbstractArray{T,N}, I::AbstractArray{T,N})
function getindex{T,N,TCoefs,IT,EB}(itp::Interpolation{T,N,TCoefs,IT,EB}, I::AbstractArray{T})
error("Array indexing is not defined for interpolation objects.")
end

# Resolve ambiguity with colon indexing for 1D interpolations
# getindex{T}(A::AbstractArray{T,1}, C::Colon)
function getindex{T,TCoefs,IT,EB}(itp::Interpolation{T,1,TCoefs,IT,EB}, c::Colon)
error("Colon indexing is not supported for interpolation objects")
end

# Allow multilinear indexing with any types
stagedfunction getindex{T,N,TCoefs,TIndex,IT,EB}(itp::Interpolation{T,N,TCoefs,IT,EB}, x::TIndex...)
n = length(x)
n == N || return :(error("Must index interpolation objects with as many indexes as dimensions"))
Copy link
Contributor

Choose a reason for hiding this comment

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

It would be nice if these error messages also showed dims(itp) (i.e. N) and n, to make it easier to spot bugs (in the client code) where the user has several interpolation objects of different dimensions, and indexes into the wrong one because of a typo.

Copy link
Member Author

Choose a reason for hiding this comment

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

Good idea.

getindex_impl(N, IT, EB)
end

# Define in-place gradient calculation
# gradient!(g::Vector, itp::Interpolation, NTuple{N}...)
stagedfunction gradient!{T,N,TCoefs,TOut,IT,EB}(g::Vector{TOut}, itp::Interpolation{T,N,TCoefs,IT,EB}, x...)
n = length(x)
n == N || return :(error("Must index interpolation objects with as many indexes as dimensions"))
it = IT()
eb = EB()
gr = gridrepresentation(it)
q = quote
length(g) == $N || error("g must be an array with exactly N elements (length(g) == "*string(length(g))*", N == "*string($N)*")")
@nexprs $N d->(x_d = x[d])
$(extrap_transform_x(gr,eb,N))
$(define_indices(it,N))
@nexprs $N dim->begin
@nexprs $N d->begin
(d==dim ? $(gradient_coefficients(it,N,:d))
: $(coefficients(it,N,:d)))
end
))

@inbounds g[dim] = $(index_gen(degree(it),N))
end
@show g
g
end
# @show q
q
end

gradient{T,N}(itp::Interpolation{T,N}, x...) = gradient!(Array(T,N), itp, x...)

# This creates prefilter specializations for all interpolation types that need them
for IT in (
Quadratic{Flat,OnCell},
Quadratic{Flat,OnGrid},
Quadratic{Reflect,OnCell},
Quadratic{Reflect,OnGrid},
Quadratic{Line,OnGrid},
Quadratic{Line,OnCell},
Quadratic{Free,OnGrid},
Quadratic{Free,OnCell},
Quadratic{Periodic,OnGrid},
Quadratic{Periodic,OnCell},
)
@ngenerate N Array{TWeights,N} function prefilter{TWeights,TCoefs,N}(::Type{TWeights}, A::Array{TCoefs,N}, it::IT)
stagedfunction prefilter{TWeights,TCoefs,N,IT<:Quadratic}(::Type{TWeights}, A::Array{TCoefs,N}, it::IT)
quote
ret, pad = copy_with_padding(A, it)

szs = collect(size(ret))
strds = collect(strides(ret))

for dim in 1:N
n = szs[dim]
for dim in 1:$N
n = szs[dim]
Copy link
Contributor

Choose a reason for hiding this comment

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

Whitespace error (one too many).

szs[dim] = 1

M, b = prefiltering_system(TWeights, TCoefs, n, it)

@nloops N i d->1:szs[d] begin
cc = @ntuple N i
@nloops $N i d->1:szs[d] begin
cc = @ntuple $N i

strt_diff = sum([(cc[i]-1)*strds[i] for i in 1:length(cc)])
strt = 1 + strt_diff
Expand Down
3 changes: 2 additions & 1 deletion test/gradient.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,13 +22,14 @@ end
itp1 = Interpolation(Float64[f1(x) for x in 1:nx-1],
Linear(OnGrid()), ExtrapPeriodic())
for x in 2.5:nx-1.5
@show gradient!(g, itp1, x)
@test_approx_eq_eps g1(x) gradient(itp1, x)[1] abs(.1*g1(x))
@test_approx_eq_eps g1(x) gradient!(g, itp1, x)[1] abs(.1*g1(x))
@test_approx_eq_eps g1(x) g[1] abs(.1*g1(x))
end

# Since Quadratic is OnCell in the domain, check gradients at grid points
itp1 = Interpolation(Float64[f1(x) for x in 1:nx-1],
itp1 = Interpolation(Float64[f1(x) for x in 1:nx-1],
Quadratic(Periodic(),OnGrid()), ExtrapPeriodic())
for x in 2:nx-1
@test_approx_eq_eps g1(x) gradient(itp1, x)[1] abs(.05*g1(x))
Expand Down
6 changes: 3 additions & 3 deletions test/linear.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ A = Float64[f(x) for x in 1:xmax]

itp1 = Interpolation(A, Linear(OnGrid()), ExtrapError())

for x in [1:.2:xmax]
for x in 1:.2:xmax
@test_approx_eq_eps f(x) itp1[x] abs(.1*f(x))
end

Expand All @@ -23,7 +23,7 @@ end

itp2 = Interpolation(A, Linear(OnGrid()), ExtrapNaN())

for x in [1:.2:xmax]
for x in 1:.2:xmax
@test_approx_eq_eps f(x) itp2[x] abs(.1*f(x))
end

Expand All @@ -34,7 +34,7 @@ xlo, xhi = itp2[.9], itp2[xmax+.2]
## ExtrapConstant

itp3 = Interpolation(A, Linear(OnGrid()), ExtrapConstant())
for x in [3.1:.2:4.3]
for x in 3.1:.2:4.3
@test_approx_eq_eps f(x) itp3[x] abs(.1*f(x))
end

Expand Down
8 changes: 4 additions & 4 deletions test/quadratic.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ A = Float64[f(x) for x in 1:xmax]

itp1 = Interpolation(A, Quadratic(Flat(),OnCell()), ExtrapError())

for x in [3.1:.2:4.3]
for x in 3.1:.2:4.3
@test_approx_eq_eps f(x) itp1[x] abs(.1*f(x))
end

Expand All @@ -23,7 +23,7 @@ end

itp2 = Interpolation(A, Quadratic(Flat(),OnCell()), ExtrapNaN())

for x in [3.1:.2:4.3]
for x in 3.1:.2:4.3
@test_approx_eq_eps f(x) itp2[x] abs(.1*f(x))
end

Expand All @@ -40,7 +40,7 @@ itp3 = Interpolation(A, Quadratic(Flat(),OnGrid()), ExtrapConstant())

# Check inbounds and extrap values

for x in [3.1:.2:4.3]
for x in 3.1:.2:4.3
@test_approx_eq_eps f(x) itp1[x] abs(.1*f(x))
end

Expand All @@ -49,7 +49,7 @@ xlo, xhi = itp3[.2], itp3[xmax+.7]
@test_approx_eq xhi A[end]

# Check continuity
xs = [0:.1:length(A)+1]
xs = [0:.1:length(A)+1;]
Copy link
Contributor

Choose a reason for hiding this comment

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

Is the ; an idiom for something here (like getting a row vector instead of a col vector)? Why do we need that? As far as I can tell we're only using xs one value at a time, so dimensionality shouldn't matter as long as linear indexing works.

Copy link
Member Author

Choose a reason for hiding this comment

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

JuliaLang/julia#8599

Copy link
Contributor

Choose a reason for hiding this comment

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

Ah. I had completely missed out on that discussion :P

I must admit that I don't like the syntax very much, although I realize I'm likely to be overridden by consensus here (and it might grow on me with time). Would I be considered terribly conservative if I perferred [i-1 for i in 1:.1:length(x)] (or [(0:.1:length(x)-1)...]) over [0:.1:length(x)-1;]? I just think the semicolon is a little too magical...

Copy link
Member Author

Choose a reason for hiding this comment

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

vcat(rng) or convert(Vector{Float64}, rng) are also choices; which do you prefer?

Copy link
Contributor

Choose a reason for hiding this comment

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

We should actually be able to do without either; we could just as well do

for x in 0:.1:xmax-.1
    @test_approx_eq_eps itp3[x] itp3[x+.1] .1
end

which I think is much less magical. (But I'm being really nitpicky here, so I don't blame you if you ignore me...)


for i in 1:length(xs)-1
@test_approx_eq_eps itp3[xs[i]] itp3[xs[i+1]] .1
Expand Down