diff --git a/README.md b/README.md index cb7913a2..313368ff 100644 --- a/README.md +++ b/README.md @@ -99,25 +99,25 @@ B-splines of quadratic or higher degree require solving an equation system to ob Some examples: ```jl # Nearest-neighbor interpolation -itp = interpolate(a, BSpline{Constant}, OnCell) +itp = interpolate(a, BSpline(Constant()), OnCell()) v = itp[5.4] # returns a[5] # (Multi)linear interpolation -itp = interpolate(A, BSpline{Linear}, OnGrid) +itp = interpolate(A, BSpline(Linear()), OnGrid()) v = itp[3.2, 4.1] # returns 0.9*(0.8*A[3,4]+0.2*A[4,4]) + 0.1*(0.8*A[3,5]+0.2*A[4,5]) # Quadratic interpolation with reflecting boundary conditions # Quadratic is the lowest order that has continuous gradient -itp = interpolate(A, BSpline{Quadratic{Reflect}}, OnCell) +itp = interpolate(A, BSpline(Quadratic(Reflect())), OnCell()) # Linear interpolation in the first dimension, and no interpolation (just lookup) in the second -itp = interpolate(A, Tuple{BSpline{Linear}, BSpline{NoInterp}}, OnGrid) +itp = interpolate(A, (BSpline(Linear()), NoInterp()), OnGrid()) v = itp[3.65, 5] # returns 0.35*A[3,5] + 0.65*A[4,5] ``` There are more options available, for example: ```jl # In-place interpolation -itp = interpolate!(A, BSpline{Quadratic{InPlace}}, OnCell) +itp = interpolate!(A, BSpline(Quadratic(InPlace())), OnCell()) ``` which destroys the input `A` but also does not need to allocate as much memory. @@ -133,7 +133,7 @@ In 1D A = rand(20) A_x = collect(1.:40.) knots = (a,) -itp = interpolate(knots,A, Gridded{Linear}) +itp = interpolate(knots,A, Gridded(Linear())) itp[2.0] ``` @@ -148,22 +148,22 @@ For example: ```jl A = rand(8,20) knots = ([x^2 for x = 1:8], [0.2y for y = 1:20]) -itp = interpolate(knots, A, Gridded{Linear}) +itp = interpolate(knots, A, Gridded(Linear())) itp[4,1.2] # approximately A[2,6] ``` One may also mix modes, by specifying a mode vector in the form of an explicit tuple: ```jl -itp = interpolate(knots, A, Tuple{Gridded{Linear},Gridded{Constant}}) +itp = interpolate(knots, A, (Gridded(Linear()),Gridded(Constant()))) ``` Presently there are only three modes for gridded: ```jl -Gridded{Linear} +Gridded(Linear()) ``` whereby a linear interpolation is applied between knots, ```jl -Gridded{Constant} +Gridded(Constant()) ``` whereby nearest neighbor interpolation is used on the applied axis, ```jl @@ -174,11 +174,7 @@ coordinates results in the throwing of an error. ## Extrapolation -The call to `extrapolate` defines what happens if you try to index into the interpolation object with coordinates outside of `[1, size(data,d)]` in any dimension `d`. The implemented boundary conditions are `Throw` and `Flat`, with more options planned. - -## More examples - -There's an [IJulia notebook](http://nbviewer.ipython.org/github/tlycken/Interpolations.jl/blob/master/doc/Interpolations.jl.ipynb) that shows off some of the current functionality, and outlines where we're headed. I try to keep it up to date when I make any significant improvements and/or breaking changes, but if it's not, do file a PR. +The call to `extrapolate` defines what happens if you try to index into the interpolation object with coordinates outside of `[1, size(data,d)]` in any dimension `d`. The implemented boundary conditions are `Throw`, `Flat`, `Linear`, `Periodic` and `Reflect`, with more options planned. `Periodic` and `Reflect` require that there is a method of `Base.mod` that can handle the indices used. ## Performance shootout diff --git a/doc/devdocs.md b/doc/devdocs.md index 7f3a3f5e..c6bf8dd3 100644 --- a/doc/devdocs.md +++ b/doc/devdocs.md @@ -18,7 +18,7 @@ First let's create an interpolation object: 0.134746 0.430876 - julia> yitp = interpolate(A, BSpline{Linear}, OnGrid) + julia> yitp = interpolate(A, BSpline(Linear()), OnGrid()) 5-element Interpolations.BSplineInterpolation{Float64,1,Float64,Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid}: 0.74838 0.995383 diff --git a/src/b-splines/b-splines.jl b/src/b-splines/b-splines.jl index 7a3103ea..c00fda78 100644 --- a/src/b-splines/b-splines.jl +++ b/src/b-splines/b-splines.jl @@ -9,14 +9,14 @@ export abstract Degree{N} immutable BSpline{D<:Degree} <: InterpolationType end -BSpline{D<:Degree}(::Type{D}) = BSpline{D} +BSpline{D<:Degree}(::D) = BSpline{D}() bsplinetype{D<:Degree}(::Type{BSpline{D}}) = D immutable BSplineInterpolation{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},pad} <: AbstractInterpolation{T,N,IT,GT} coefs::Array{TCoefs,N} end -function BSplineInterpolation{N,TCoefs,TWeights<:Real,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},pad}(::Type{TWeights}, A::AbstractArray{TCoefs,N}, ::Type{IT}, ::Type{GT}, ::Val{pad}) +function BSplineInterpolation{N,TCoefs,TWeights<:Real,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},pad}(::Type{TWeights}, A::AbstractArray{TCoefs,N}, ::IT, ::GT, ::Val{pad}) isleaftype(IT) || error("The b-spline type must be a leaf type (was $IT)") isleaftype(TCoefs) || warn("For performance reasons, consider using an array of a concrete type (eltype(A) == $(eltype(A)))") @@ -49,12 +49,12 @@ count_interp_dims{T,N,TCoefs,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType end end -function interpolate{TWeights,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(::Type{TWeights}, ::Type{TCoefs}, A, ::Type{IT}, ::Type{GT}) +function interpolate{TWeights,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(::Type{TWeights}, ::Type{TCoefs}, A, it::IT, gt::GT) Apad, Pad = prefilter(TWeights, TCoefs, A, IT, GT) - BSplineInterpolation(TWeights, Apad, IT, GT, Pad) + BSplineInterpolation(TWeights, Apad, it, gt, Pad) end -function interpolate{IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(A::AbstractArray, ::Type{IT}, ::Type{GT}) - interpolate(tweight(A), tcoef(A), A, IT, GT) +function interpolate{IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(A::AbstractArray, it::IT, gt::GT) + interpolate(tweight(A), tcoef(A), A, it, gt) end # We can't just return a tuple-of-types due to julia #12500 @@ -68,9 +68,9 @@ tcoef(A::AbstractArray{Float32}) = Float32 tcoef(A::AbstractArray{Rational{Int}}) = Rational{Int} tcoef{T<:Integer}(A::AbstractArray{T}) = typeof(float(zero(T))) -interpolate!{TWeights,IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(::Type{TWeights}, A, ::Type{IT}, ::Type{GT}) = BSplineInterpolation(TWeights, prefilter!(TWeights, A, IT, GT), IT, GT, Val{0}()) -function interpolate!{IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(A::AbstractArray, ::Type{IT}, ::Type{GT}) - interpolate!(tweight(A), A, IT, GT) +interpolate!{TWeights,IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(::Type{TWeights}, A, it::IT, gt::GT) = BSplineInterpolation(TWeights, prefilter!(TWeights, A, IT, GT), it, gt, Val{0}()) +function interpolate!{IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(A::AbstractArray, it::IT, gt::GT) + interpolate!(tweight(A), A, it, gt) end include("constant.jl") diff --git a/src/b-splines/quadratic.jl b/src/b-splines/quadratic.jl index e3ba49d0..dca256a7 100644 --- a/src/b-splines/quadratic.jl +++ b/src/b-splines/quadratic.jl @@ -1,5 +1,5 @@ immutable Quadratic{BC<:BoundaryCondition} <: Degree{2} end -Quadratic{BC<:BoundaryCondition}(::Type{BC}) = Quadratic{BC} +Quadratic{BC<:BoundaryCondition}(::BC) = Quadratic{BC}() function define_indices_d{BC}(::Type{BSpline{Quadratic{BC}}}, d, pad) symix, symixm, symixp = symbol("ix_",d), symbol("ixm_",d), symbol("ixp_",d) diff --git a/src/extrapolation/extrapolation.jl b/src/extrapolation/extrapolation.jl index edee18d7..b893fd44 100644 --- a/src/extrapolation/extrapolation.jl +++ b/src/extrapolation/extrapolation.jl @@ -1,10 +1,11 @@ export Throw, - FilledExtrapolation # for direct control over typeof(fillvalue) + FillValue, + TypedFillValue # for direct control over typeof(fillvalue) type Extrapolation{T,N,ITPT,IT,GT,ET} <: AbstractInterpolationWrapper{T,N,ITPT,IT,GT} itp::ITPT end -Extrapolation{T,ITPT,IT,GT,ET}(::Type{T}, N, itp::ITPT, ::Type{IT}, ::Type{GT}, ::Type{ET}) = +Extrapolation{T,ITPT,IT,GT,ET}(::Type{T}, N, itp::ITPT, ::Type{IT}, ::Type{GT}, et::ET) = Extrapolation{T,N,ITPT,IT,GT,ET}(itp) """ @@ -15,15 +16,15 @@ The scheme can take any of these values: * `Throw` - throws a BoundsError for out-of-bounds indices * `Flat` - for constant extrapolation, taking the closest in-bounds value * `Linear` - linear extrapolation (the wrapped interpolation object must support gradient) -* `Reflect` - reflecting extrapolation -* `Periodic` - periodic extrapolation +* `Reflect` - reflecting extrapolation (indices must support `mod`) +* `Periodic` - periodic extrapolation (indices must support `mod`) -You can also combine schemes in tuples. For example, the scheme `Tuple{Linear, Flat}` will use linear extrapolation in the first dimension, and constant in the second. +You can also combine schemes in tuples. For example, the scheme `(Linear(), Flat())` will use linear extrapolation in the first dimension, and constant in the second. -Finally, you can specify different extrapolation behavior in different direction. `Tuple{Tuple{Linear,Flat}, Flat}` will extrapolate linearly in the first dimension if the index is too small, but use constant etrapolation if it is too large, and always use constant extrapolation in the second dimension. +Finally, you can specify different extrapolation behavior in different direction. `((Linear(),Flat()), Flat())` will extrapolate linearly in the first dimension if the index is too small, but use constant etrapolation if it is too large, and always use constant extrapolation in the second dimension. """ -extrapolate{T,N,IT,GT,ET}(itp::AbstractInterpolation{T,N,IT,GT}, ::Type{ET}) = - Extrapolation(T,N,itp,IT,GT,ET) +extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, et) = + Extrapolation(T,N,itp,IT,GT,et) include("throw.jl") include("flat.jl") diff --git a/src/extrapolation/filled.jl b/src/extrapolation/filled.jl index 535cb899..0560d879 100644 --- a/src/extrapolation/filled.jl +++ b/src/extrapolation/filled.jl @@ -1,23 +1,36 @@ nindexes(N::Int) = N == 1 ? "1 index" : "$N indexes" +abstract AbstractFillValue +immutable FillValue{T} <: AbstractFillValue + val::T +end +immutable TypedFillValue{T} <: AbstractFillValue + val::T +end +TypedFillValue{T}(::Type{T}, v) = TypedFillValue{T}(convert(T, v)) -type FilledExtrapolation{T,N,ITP<:AbstractInterpolation,IT,GT,FT} <: AbstractExtrapolation{T,N,ITP,IT,GT} +type FilledExtrapolation{T,N,ITP<:AbstractInterpolation,IT,GT,FT<:AbstractFillValue} <: 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. +FilledExtrapolation{T,N,IT,GT,FT}(itp::AbstractInterpolation{T,N,IT,GT}, fv::FT) = FilledExtrapolation{T,N,typeof(itp),IT,GT,FT}(itp,fv) +# """ +# `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 +# 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::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(itp, fillvalue)` creates an extrapolation object that returns the `fillvalue` any time the indexes in `itp[x1,x2,...]` are out-of-bounds. The type of the fill value will be converted to the eltype of the interpolation object, so that indexing into the extrapolation object is type-stable. +""" +extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue::FillValue) = FilledExtrapolation(itp, FillValue(convert(T, fillvalue.val))) +""" +`extrapolate(itp, fillvalue<:TypedFillValue)` creates an extrapolation object that returns the `fillvalue` any time the indexes in `itp[x1,x2,...]` are out-of-bounds. The type of the original fill value is not changed, so it is possible to do `extrapolate(itp, TypedFillValue("hello"))` - although indexing into that extrapolation object will now be type unstable. """ -extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue) = FilledExtrapolation(itp, convert(eltype(itp), fillvalue)) +extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue::TypedFillValue) = FilledExtrapolation(itp, fillvalue) @generated function getindex{T,N}(fitp::FilledExtrapolation{T,N}, args::Number...) n = length(args) @@ -27,7 +40,7 @@ extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue) = Fille $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) + @nexprs $N d->((args[d] < lbound(fitp,d) || args[d] > ubound(fitp, d)) && return fitp.fillvalue.val) # In the interpolation region return getindex(fitp.itp,args...) end diff --git a/src/gridded/gridded.jl b/src/gridded/gridded.jl index 960926df..72c29cd2 100644 --- a/src/gridded/gridded.jl +++ b/src/gridded/gridded.jl @@ -1,7 +1,7 @@ export Gridded immutable Gridded{D<:Degree} <: InterpolationType end -Gridded{D<:Degree}(::Type{D}) = Gridded{D} +Gridded{D<:Degree}(::D) = Gridded{D}() griddedtype{D<:Degree}(::Type{Gridded{D}}) = D @@ -13,7 +13,7 @@ immutable GriddedInterpolation{T,N,TCoefs,IT<:DimSpec{Gridded},K<:Tuple{Vararg{V knots::K coefs::Array{TCoefs,N} end -function GriddedInterpolation{N,TCoefs,TWeights<:Real,IT<:DimSpec{Gridded},pad}(::Type{TWeights}, knots::NTuple{N,GridIndex}, A::AbstractArray{TCoefs,N}, ::Type{IT}, ::Val{pad}) +function GriddedInterpolation{N,TCoefs,TWeights<:Real,IT<:DimSpec{Gridded},pad}(::Type{TWeights}, knots::NTuple{N,GridIndex}, A::AbstractArray{TCoefs,N}, ::IT, ::Val{pad}) isleaftype(IT) || error("The b-spline type must be a leaf type (was $IT)") isleaftype(TCoefs) || warn("For performance reasons, consider using an array of a concrete type (eltype(A) == $(eltype(A)))") @@ -48,16 +48,16 @@ iextract{T<:GridType}(::Type{T}, d) = T end end -function interpolate{TWeights,TCoefs,Tel,N,IT<:DimSpec{Gridded}}(::Type{TWeights}, ::Type{TCoefs}, knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, ::Type{IT}) - GriddedInterpolation(TWeights, knots, A, IT, Val{0}()) +function interpolate{TWeights,TCoefs,Tel,N,IT<:DimSpec{Gridded}}(::Type{TWeights}, ::Type{TCoefs}, knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT) + GriddedInterpolation(TWeights, knots, A, it, Val{0}()) end -function interpolate{Tel,N,IT<:DimSpec{Gridded}}(knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, ::Type{IT}) - interpolate(tweight(A), tcoef(A), knots, A, IT) +function interpolate{Tel,N,IT<:DimSpec{Gridded}}(knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT) + interpolate(tweight(A), tcoef(A), knots, A, it) end -interpolate!{TWeights,Tel,N,IT<:DimSpec{Gridded}}(::Type{TWeights}, knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, ::Type{IT}) = GriddedInterpolation(TWeights, knots, A, IT, Val{0}()) -function interpolate!{Tel,N,IT<:DimSpec{Gridded}}(knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, ::Type{IT}) - interpolate!(tweight(A), tcoef(A), knots, A, IT) +interpolate!{TWeights,Tel,N,IT<:DimSpec{Gridded}}(::Type{TWeights}, knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT) = GriddedInterpolation(TWeights, knots, A, it, Val{0}()) +function interpolate!{Tel,N,IT<:DimSpec{Gridded}}(knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT) + interpolate!(tweight(A), tcoef(A), knots, A, it) end include("constant.jl") diff --git a/test/b-splines/constant.jl b/test/b-splines/constant.jl index 3f8dff28..2daaca05 100644 --- a/test/b-splines/constant.jl +++ b/test/b-splines/constant.jl @@ -10,12 +10,12 @@ A2 = rand(Float64, N1, N1) * 100 A3 = rand(Float64, N1, N1, N1) * 100 for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy)) - itp1c = @inferred(constructor(copier(A1), BSpline(Constant), OnCell)) - itp1g = @inferred(constructor(copier(A1), BSpline(Constant), OnGrid)) - itp2c = @inferred(constructor(copier(A2), BSpline(Constant), OnCell)) - itp2g = @inferred(constructor(copier(A2), BSpline(Constant), OnGrid)) - itp3c = @inferred(constructor(copier(A3), BSpline(Constant), OnCell)) - itp3g = @inferred(constructor(copier(A3), BSpline(Constant), OnGrid)) + itp1c = @inferred(constructor(copier(A1), BSpline(Constant()), OnCell())) + itp1g = @inferred(constructor(copier(A1), BSpline(Constant()), OnGrid())) + itp2c = @inferred(constructor(copier(A2), BSpline(Constant()), OnCell())) + itp2g = @inferred(constructor(copier(A2), BSpline(Constant()), OnGrid())) + itp3c = @inferred(constructor(copier(A3), BSpline(Constant()), OnCell())) + itp3g = @inferred(constructor(copier(A3), BSpline(Constant()), OnGrid())) # Evaluation on provided data points # 1D diff --git a/test/b-splines/linear.jl b/test/b-splines/linear.jl index 6f6702c0..4062ddc9 100644 --- a/test/b-splines/linear.jl +++ b/test/b-splines/linear.jl @@ -10,7 +10,7 @@ f(x) = g1(x) A1 = Float64[f(x) for x in 1:xmax] for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy)) - itp1c = @inferred(constructor(copier(A1), BSpline(Linear), OnCell)) + itp1c = @inferred(constructor(copier(A1), BSpline(Linear()), OnCell())) # Just interpolation for x in 1:.2:xmax @@ -20,7 +20,7 @@ for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy)) # Rational element types fr(x) = (x^2) // 40 + 2 A1R = Rational{Int}[fr(x) for x in 1:10] - itp1r = @inferred(constructor(copier(A1R), BSpline(Linear), OnGrid)) + itp1r = @inferred(constructor(copier(A1R), BSpline(Linear()), OnGrid())) @test @inferred(size(itp1r)) == size(A1R) @test_approx_eq_eps itp1r[23//10] fr(23//10) abs(.1*fr(23//10)) @test typeof(itp1r[23//10]) == Rational{Int} @@ -31,7 +31,7 @@ for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy)) f(x,y) = g1(x)*g2(y) ymax = 10 A2 = Float64[f(x,y) for x in 1:xmax, y in 1:ymax] - itp2 = @inferred(constructor(copier(A2), BSpline(Linear), OnGrid)) + itp2 = @inferred(constructor(copier(A2), BSpline(Linear()), OnGrid())) @test @inferred(size(itp2)) == size(A2) for x in 2.1:.2:xmax-1, y in 1.9:.2:ymax-.9 diff --git a/test/b-splines/mixed.jl b/test/b-splines/mixed.jl index cc6ab9ab..eb1131a9 100644 --- a/test/b-splines/mixed.jl +++ b/test/b-splines/mixed.jl @@ -7,8 +7,8 @@ N = 10 for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy)) A2 = rand(Float64, N, N) * 100 for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell) - itp_a = @inferred(constructor(copier(A2), Tuple{BSpline(Linear), BSpline(Quadratic(BC))}, GT)) - itp_b = @inferred(constructor(copier(A2), Tuple{BSpline(Quadratic(BC)), BSpline(Linear)}, GT)) + itp_a = @inferred(constructor(copier(A2), (BSpline(Linear()), BSpline(Quadratic(BC()))), GT())) + itp_b = @inferred(constructor(copier(A2), (BSpline(Quadratic(BC())), BSpline(Linear())), GT())) @test @inferred(size(itp_a)) == size(A2) @test @inferred(size(itp_b)) == size(A2) diff --git a/test/b-splines/multivalued.jl b/test/b-splines/multivalued.jl index b4eca859..2c344333 100644 --- a/test/b-splines/multivalued.jl +++ b/test/b-splines/multivalued.jl @@ -23,20 +23,20 @@ Base.promote_rule{T1,T2<:Number}(::Type{MyPair{T1}}, ::Type{T2}) = MyPair{promot # 1d A = reinterpret(MyPair{Float64}, rand(2, 10), (10,)) -itp = interpolate(A, BSpline(Constant), OnGrid) +itp = interpolate(A, BSpline(Constant()), OnGrid()) itp[3.2] -itp = interpolate(A, BSpline(Linear), OnGrid) +itp = interpolate(A, BSpline(Linear()), OnGrid()) itp[3.2] -itp = interpolate(A, BSpline(Quadratic(Flat)), OnGrid) +itp = interpolate(A, BSpline(Quadratic(Flat())), OnGrid()) itp[3.2] # 2d A = reinterpret(MyPair{Float64}, rand(2, 10, 5), (10,5)) -itp = interpolate(A, BSpline(Constant), OnGrid) +itp = interpolate(A, BSpline(Constant()), OnGrid()) itp[3.2,1.8] -itp = interpolate(A, BSpline(Linear), OnGrid) +itp = interpolate(A, BSpline(Linear()), OnGrid()) itp[3.2,1.8] -itp = interpolate(A, BSpline(Quadratic(Flat)), OnGrid) +itp = interpolate(A, BSpline(Quadratic(Flat())), OnGrid()) itp[3.2,1.8] end diff --git a/test/b-splines/quadratic.jl b/test/b-splines/quadratic.jl index ca33dd32..87cc6f43 100644 --- a/test/b-splines/quadratic.jl +++ b/test/b-splines/quadratic.jl @@ -7,7 +7,7 @@ for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy)) xmax = 10 A = Float64[f(x) for x in 1:xmax] for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell) - itp1 = @inferred(constructor(copier(A), BSpline(Quadratic(BC)), GT)) + itp1 = @inferred(constructor(copier(A), BSpline(Quadratic(BC())), GT())) @test @inferred(size(itp1)) == size(A) # test that inner region is close to data @@ -37,7 +37,7 @@ for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy)) # test that inner region is close to data for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell) - itp2 = @inferred(constructor(copier(A), BSpline(Quadratic(BC)), GT)) + itp2 = @inferred(constructor(copier(A), BSpline(Quadratic(BC())), GT())) @test @inferred(size(itp2)) == size(A) for x in 3.1:.2:xmax-3, y in 3.1:2:ymax-3 @@ -50,7 +50,7 @@ let f(x) = sin((x-3)*2pi/9 - 1) xmax = 10 A = Float64[f(x) for x in 1:xmax] - itp1 = interpolate!(copy(A), BSpline(Quadratic(InPlace)), OnCell) + itp1 = interpolate!(copy(A), BSpline(Quadratic(InPlace())), OnCell()) for i = 1:xmax @test_approx_eq itp1[i] A[i] end @@ -58,7 +58,7 @@ let f(x,y) = sin(x/10)*cos(y/6) xmax, ymax = 30,10 A = Float64[f(x,y) for x in 1:xmax, y in 1:ymax] - itp2 = interpolate!(copy(A), BSpline(Quadratic(InPlace)), OnCell) + itp2 = interpolate!(copy(A), BSpline(Quadratic(InPlace())), OnCell()) for j = 1:ymax, i = 1:xmax @test_approx_eq itp2[i,j] A[i,j] end diff --git a/test/extrapolation/runtests.jl b/test/extrapolation/runtests.jl index dfea4543..0626f2b5 100644 --- a/test/extrapolation/runtests.jl +++ b/test/extrapolation/runtests.jl @@ -7,14 +7,14 @@ f(x) = sin((x-3)*2pi/9 - 1) xmax = 10 A = Float64[f(x) for x in 1:xmax] -itpg = interpolate(A, BSpline(Linear), OnGrid) +itpg = interpolate(A, BSpline(Linear()), OnGrid()) -etpg = extrapolate(itpg, Flat) +etpg = extrapolate(itpg, Flat()) @test etpg[-3] == etpg[-4.5] == etpg[0.9] == etpg[1.0] == A[1] @test etpg[10.1] == etpg[11] == etpg[148.298452] == A[end] -etpf = @inferred(extrapolate(itpg, NaN)) +etpf = @inferred(extrapolate(itpg, FillValue(NaN))) @test @inferred(size(etpf)) == (xmax,) @test isnan(@inferred(getindex(etpf, -2.5))) @@ -27,11 +27,11 @@ etpf = @inferred(extrapolate(itpg, NaN)) @test_throws BoundsError etpf[2.5,2] @test_throws ErrorException etpf[2.5,2,1] # this will probably become a BoundsError someday -etpf = @inferred(FilledExtrapolation(itpg, 'x')) +etpf = @inferred(extrapolate(itpg, TypedFillValue('x'))) @test_approx_eq etpf[2] f(2) @test etpf[-1.5] == 'x' -etpl = extrapolate(itpg, Linear) +etpl = extrapolate(itpg, Linear()) k_lo = A[2] - A[1] x_lo = -3.2 @test_approx_eq etpl[x_lo] A[1]+k_lo*(x_lo-1) @@ -43,19 +43,19 @@ x_hi = xmax + 5.7 xmax, ymax = 8,8 g(x, y) = (x^2 + 3x - 8) * (-2y^2 + y + 1) -itp2g = interpolate(Float64[g(x,y) for x in 1:xmax, y in 1:ymax], Tuple{BSpline(Quadratic(Free)), BSpline(Linear)}, OnGrid) -etp2g = extrapolate(itp2g, Tuple{Linear, Flat}) +itp2g = interpolate(Float64[g(x,y) for x in 1:xmax, y in 1:ymax], (BSpline(Quadratic(Free())), BSpline(Linear())), OnGrid()) +etp2g = extrapolate(itp2g, (Linear(), Flat())) @test_approx_eq @inferred(getindex(etp2g, -.5, 4)) itp2g[1,4]-1.5*epsilon(etp2g[dual(1,1),4]) @test_approx_eq @inferred(getindex(etp2g, 5, 100)) itp2g[5,ymax] -etp2ud = extrapolate(itp2g, Tuple{Tuple{Linear, Flat}, Flat}) +etp2ud = extrapolate(itp2g, ((Linear(), Flat()), Flat())) @test_approx_eq @inferred(getindex(etp2ud, -.5, 4)) itp2g[1,4] - 1.5*epsilon(etp2g[dual(1,1),4]) @test @inferred(getindex(etp2ud, 5, -4)) == etp2ud[5,1] @test @inferred(getindex(etp2ud, 100, 4)) == etp2ud[8,4] @test @inferred(getindex(etp2ud, -.5, 100)) == itp2g[1,8] - 1.5 * epsilon(etp2g[dual(1,1),8]) -etp2ll = extrapolate(itp2g, Linear) +etp2ll = extrapolate(itp2g, Linear()) @test_approx_eq @inferred(getindex(etp2ll, -.5, 100)) itp2g[1,8]-1.5*epsilon(etp2ll[dual(1,1),8]) + (100-8) * epsilon(etp2ll[1,dual(8,1)]) diff --git a/test/extrapolation/type-stability.jl b/test/extrapolation/type-stability.jl index 1e6b4062..4c4c444a 100644 --- a/test/extrapolation/type-stability.jl +++ b/test/extrapolation/type-stability.jl @@ -6,7 +6,7 @@ using Base.Test, Interpolations, DualNumbers f(x) = sin((x-3)*2pi/9 - 1) xmax = 10 A = Float64[f(x) for x in 1:xmax] -itpg = interpolate(A, BSpline(Linear), OnGrid) +itpg = interpolate(A, BSpline(Linear()), OnGrid()) schemes = ( Flat, @@ -15,7 +15,7 @@ schemes = ( Periodic ) -for etp in map(E -> extrapolate(itpg, E), schemes), +for etp in map(E -> extrapolate(itpg, E()), schemes), x in [ # In-bounds evaluation 3.4, 3, dual(3.1), @@ -30,9 +30,9 @@ end g(y) = (y/100)^3 ymax = 4 A = Float64[f(x)*g(y) for x in 1:xmax, y in 1:ymax] -itp2 = interpolate(A, BSpline(Linear), OnGrid) +itp2 = interpolate(A, BSpline(Linear()), OnGrid()) -for (etp2,E) in map(E -> (extrapolate(itp2, E), E), schemes), +for (etp2,E) in map(E -> (extrapolate(itp2, E()), E), schemes), x in ( # In-bounds evaluation 3.4, 3, dual(3.1), diff --git a/test/gradient.jl b/test/gradient.jl index 9386a73e..f69ce69b 100644 --- a/test/gradient.jl +++ b/test/gradient.jl @@ -8,7 +8,7 @@ g1(x) = 2pi/(nx-1) * cos((x-3)*2pi/(nx-1) - 1) # Gradient of Constant should always be 0 itp1 = interpolate(Float64[f1(x) for x in 1:nx-1], - BSpline(Constant), OnGrid) + BSpline(Constant()), OnGrid()) g = Array(Float64, 1) @@ -20,7 +20,7 @@ end # Since Linear is OnGrid in the domain, check the gradients between grid points itp1 = interpolate(Float64[f1(x) for x in 1:nx-1], - BSpline(Linear), OnGrid) + BSpline(Linear()), OnGrid()) for x in 2.5:nx-1.5 @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)) @@ -36,7 +36,7 @@ end # Since Quadratic is OnCell in the domain, check gradients at grid points itp1 = interpolate(Float64[f1(x) for x in 1:nx-1], - BSpline(Quadratic(Periodic)), OnCell) + BSpline(Quadratic(Periodic())), OnCell()) for x in 2:nx-1 @test_approx_eq_eps g1(x) gradient(itp1, x)[1] abs(.05*g1(x)) @test_approx_eq_eps g1(x) gradient!(g, itp1, x)[1] abs(.05*g1(x)) @@ -61,7 +61,7 @@ dqfunc = x -> 2*a*(x-c) xg = Float64[1:5;] y = qfunc(xg) -iq = interpolate(y, BSpline(Quadratic(Free)), OnCell) +iq = interpolate(y, BSpline(Quadratic(Free())), OnCell()) x = 1.8 @test_approx_eq iq[x] qfunc(x) @test_approx_eq gradient(iq, x)[1] dqfunc(x) @@ -69,14 +69,14 @@ x = 1.8 # 2d (biquadratic) p = [(x-1.75)^2 for x = 1:7] A = p*p' -iq = interpolate(A, BSpline(Quadratic(Free)), OnCell) +iq = interpolate(A, BSpline(Quadratic(Free())), OnCell()) @test_approx_eq iq[4,4] (4-1.75)^4 @test_approx_eq iq[4,3] (4-1.75)^2*(3-1.75)^2 g = gradient(iq, 4, 3) @test_approx_eq g[1] 2*(4-1.75)*(3-1.75)^2 @test_approx_eq g[2] 2*(4-1.75)^2*(3-1.75) -iq = interpolate!(copy(A), BSpline(Quadratic(InPlace)), OnCell) +iq = interpolate!(copy(A), BSpline(Quadratic(InPlace())), OnCell()) @test_approx_eq iq[4,4] (4-1.75)^4 @test_approx_eq iq[4,3] (4-1.75)^2*(3-1.75)^2 g = gradient(iq, 4, 3) @@ -84,7 +84,7 @@ g = gradient(iq, 4, 3) @test_approx_eq_eps g[2] 2*(4-1.75)^2*(3-1.75) 0.2 # InPlaceQ is exact for an underlying quadratic -iq = interpolate!(copy(A), BSpline(Quadratic(InPlaceQ)), OnCell) +iq = interpolate!(copy(A), BSpline(Quadratic(InPlaceQ())), OnCell()) @test_approx_eq iq[4,4] (4-1.75)^4 @test_approx_eq iq[4,3] (4-1.75)^2*(3-1.75)^2 g = gradient(iq, 4, 3) @@ -93,10 +93,10 @@ g = gradient(iq, 4, 3) A2 = rand(Float64, nx, nx) * 100 for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell) - itp_a = interpolate(A2, Tuple{BSpline(Linear), BSpline(Quadratic(BC))}, GT) - itp_b = interpolate(A2, Tuple{BSpline(Quadratic(BC)), BSpline(Linear)}, GT) - itp_c = interpolate(A2, Tuple{NoInterp, BSpline(Quadratic(BC))}, GT) - itp_d = interpolate(A2, Tuple{BSpline(Quadratic(BC)), NoInterp}, GT) + itp_a = interpolate(A2, (BSpline(Linear()), BSpline(Quadratic(BC()))), GT()) + itp_b = interpolate(A2, (BSpline(Quadratic(BC())), BSpline(Linear())), GT()) + itp_c = interpolate(A2, (NoInterp(), BSpline(Quadratic(BC()))), GT()) + itp_d = interpolate(A2, (BSpline(Quadratic(BC())), NoInterp()), GT()) for i = 1:10 x = rand()*(nx-2)+1.5 diff --git a/test/gridded/gridded.jl b/test/gridded/gridded.jl index 5afc813f..d6e02830 100644 --- a/test/gridded/gridded.jl +++ b/test/gridded/gridded.jl @@ -6,7 +6,7 @@ for (D,G) in ((Constant,OnCell),(Linear,OnGrid)) ## 1D a = rand(5) knots = (collect(linspace(1,length(a),length(a))),) - itp = @inferred(interpolate(knots, a, Gridded{D})) + itp = @inferred(interpolate(knots, a, Gridded(D()))) @inferred(getindex(itp, 2)) @inferred(getindex(itp, CartesianIndex((2,)))) for i = 2:length(a)-1 @@ -28,7 +28,7 @@ for (D,G) in ((Constant,OnCell),(Linear,OnGrid)) @test_approx_eq v[i] itp[x[i]] end # compare against BSpline - itpb = @inferred(interpolate(a, BSpline{D}, G)) + itpb = @inferred(interpolate(a, BSpline(D()), G())) for x in linspace(1.1,4.9,101) @test_approx_eq itp[x] itpb[x] end @@ -36,7 +36,7 @@ for (D,G) in ((Constant,OnCell),(Linear,OnGrid)) ## 2D A = rand(6,5) knots = (collect(linspace(1,size(A,1),size(A,1))),collect(linspace(1,size(A,2),size(A,2)))) - itp = @inferred(interpolate(knots, A, Gridded{D})) + itp = @inferred(interpolate(knots, A, Gridded(D()))) @inferred(getindex(itp, 2, 2)) @inferred(getindex(itp, CartesianIndex((2,2)))) for j = 2:size(A,2)-1, i = 2:size(A,1)-1 @@ -59,14 +59,14 @@ for (D,G) in ((Constant,OnCell),(Linear,OnGrid)) @test_approx_eq v[i,j] itp[x[i],y[j]] end # compare against BSpline - itpb = @inferred(interpolate(A, BSpline{D}, G)) + itpb = @inferred(interpolate(A, BSpline(D()), G())) for y in linspace(1.1,5.9,101), x in linspace(1.1,4.9,101) @test_approx_eq itp[x,y] itpb[x,y] end A = rand(8,20) knots = ([x^2 for x = 1:8], [0.2y for y = 1:20]) - itp = interpolate(knots, A, Gridded{D}) + itp = interpolate(knots, A, Gridded(D())) @test_approx_eq itp[4,1.2] A[2,6] end diff --git a/test/gridded/mixed.jl b/test/gridded/mixed.jl index df536a45..61c946f6 100644 --- a/test/gridded/mixed.jl +++ b/test/gridded/mixed.jl @@ -4,7 +4,7 @@ using Interpolations, Base.Test A = rand(6,5) knots = (collect(linspace(1,size(A,1),size(A,1))),collect(linspace(1,size(A,2),size(A,2)))) -itp = @inferred(interpolate(knots, A, Tuple{Gridded{Linear},NoInterp})) +itp = @inferred(interpolate(knots, A, (Gridded(Linear()),NoInterp()))) @inferred(getindex(itp, 2, 2)) @inferred(getindex(itp, CartesianIndex((2,2)))) for j = 2:size(A,2)-1, i = 2:size(A,1)-1 @@ -16,8 +16,8 @@ end A = rand(8,20) knots = ([x^2 for x = 1:8], [0.2y for y = 1:20]) -itp = interpolate(knots, A, Gridded{Linear}) +itp = interpolate(knots, A, Gridded(Linear())) @test_approx_eq itp[4,1.2] A[2,6] -@test_throws ErrorException interpolate(knots, A, Tuple{Gridded{Linear}, NoInterp}) +@test_throws ErrorException interpolate(knots, A, (Gridded(Linear()),NoInterp())) end diff --git a/test/issues/runtests.jl b/test/issues/runtests.jl index eb762154..22022e00 100644 --- a/test/issues/runtests.jl +++ b/test/issues/runtests.jl @@ -5,7 +5,7 @@ using Interpolations, Base.Test A = rand(1:20, 100, 100) # In #34, this incantation throws -itp = interpolate(A, BSpline(Quadratic(Flat)), OnCell) +itp = interpolate(A, BSpline(Quadratic(Flat())), OnCell()) # Sanity check that not only don't throw, but actually interpolate for i in 1:size(A,1), j in 1:size(A,2) @test_approx_eq itp[i,j] A[i,j] diff --git a/test/scaling/dimspecs.jl b/test/scaling/dimspecs.jl index 8a5349fc..4b6ff00b 100644 --- a/test/scaling/dimspecs.jl +++ b/test/scaling/dimspecs.jl @@ -6,7 +6,7 @@ xs = -pi:(2pi/10):pi-2pi/10 ys = -2:.1:2 f(x,y) = sin(x) * y^2 -itp = interpolate(Float64[f(x,y) for x in xs, y in ys], Tuple{BSpline(Quadratic(Periodic)), BSpline(Linear)}, OnGrid) +itp = interpolate(Float64[f(x,y) for x in xs, y in ys], (BSpline(Quadratic(Periodic())), BSpline(Linear())), OnGrid()) sitp = scale(itp, xs, ys) for (ix,x) in enumerate(xs), (iy,y) in enumerate(ys) diff --git a/test/scaling/nointerp.jl b/test/scaling/nointerp.jl index eeb20293..a9043dce 100644 --- a/test/scaling/nointerp.jl +++ b/test/scaling/nointerp.jl @@ -11,7 +11,7 @@ ys = 1:3 A = hcat(f1(xs), f2(xs), f3(xs)) -itp = interpolate(A, Tuple{BSpline(Quadratic(Periodic)), NoInterp}, OnGrid) +itp = interpolate(A, (BSpline(Quadratic(Periodic())), NoInterp()), OnGrid()) sitp = scale(itp, xs, ys) for (ix,x0) in enumerate(xs[1:end-1]), y0 in ys diff --git a/test/scaling/scaling.jl b/test/scaling/scaling.jl index db47bfcb..4c1acfa2 100644 --- a/test/scaling/scaling.jl +++ b/test/scaling/scaling.jl @@ -6,7 +6,7 @@ using Base.Test # Model linear interpolation of y = -3 + .5x by interpolating y=x # and then scaling to the new x range -itp = interpolate(1:1.0:10, BSpline(Linear), OnGrid) +itp = interpolate(1:1.0:10, BSpline(Linear()), OnGrid()) sitp = scale(itp, -3:.5:1.5) @@ -23,7 +23,7 @@ xs = -5:.5:5 ys = -4:.2:4 zs = Float64[testfunction(x,y) for x in xs, y in ys] -itp2 = interpolate(zs, BSpline(Quadratic(Flat)), OnGrid) +itp2 = interpolate(zs, BSpline(Quadratic(Flat())), OnGrid()) sitp2 = scale(itp2, xs, ys) for x in xs, y in ys @@ -33,7 +33,7 @@ end # Test gradients of scaled grids xs = -pi:.1:pi ys = sin(xs) -itp = interpolate(ys, BSpline(Linear), OnGrid) +itp = interpolate(ys, BSpline(Linear()), OnGrid()) sitp = scale(itp, xs) for x in -pi:.1:pi @@ -46,7 +46,7 @@ end @inferred(getindex(sitp2, -3, 1)) @inferred(getindex(sitp2, -3.4, 1)) -sitp32 = scale(interpolate(Float32[testfunction(x,y) for x in -5:.5:5, y in -4:.2:4], BSpline(Quadratic(Flat)), OnGrid), -5f0:.5f0:5f0, -4f0:.2f0:4f0) +sitp32 = scale(interpolate(Float32[testfunction(x,y) for x in -5:.5:5, y in -4:.2:4], BSpline(Quadratic(Flat())), OnGrid()), -5f0:.5f0:5f0, -4f0:.2f0:4f0) @test typeof(@inferred(getindex(sitp32, -3.4f0, 1.2f0))) == Float32 end diff --git a/test/scaling/withextrap.jl b/test/scaling/withextrap.jl index ede01fab..1b5840ad 100644 --- a/test/scaling/withextrap.jl +++ b/test/scaling/withextrap.jl @@ -24,16 +24,16 @@ function run_tests{T,N,IT}(sut::Interpolations.AbstractInterpolation{T,N,IT,OnCe end for GT in (OnGrid, OnCell) - itp = interpolate(ys, BSpline(Quadratic(Flat)), GT) + itp = interpolate(ys, BSpline(Quadratic(Flat())), GT()) # Test extrapolating, then scaling - eitp = extrapolate(itp, Flat) + eitp = extrapolate(itp, Flat()) seitp = scale(eitp, xs) run_tests(seitp, itp) # Test scaling, then extrapolating sitp = scale(itp, xs) - esitp = extrapolate(sitp, Flat) + esitp = extrapolate(sitp, Flat()) run_tests(esitp, itp) end diff --git a/test/typing.jl b/test/typing.jl index 0dbcfa39..41062202 100644 --- a/test/typing.jl +++ b/test/typing.jl @@ -8,7 +8,7 @@ g(x) = convert(Float32, 3x^2/(nx-1)) A = Float32[f(x) for x in 1:nx] -itp = interpolate(A, BSpline(Quadratic(Flat)), OnCell) +itp = interpolate(A, BSpline(Quadratic(Flat())), OnCell()) # display(plot( # layer(x=1:nx,y=[f(x) for x in 1:1//1:nx],Geom.point), @@ -29,7 +29,7 @@ end # Rational element types R = Rational{Int}[x^2//10 for x in 1:10] -itp = interpolate(R, BSpline(Quadratic(Free)), OnCell) +itp = interpolate(R, BSpline(Quadratic(Free())), OnCell()) itp[11//10] @test typeof(itp[11//10]) == Rational{Int}