From 30b598d60281667ed1893d09aea4d695abb00cc1 Mon Sep 17 00:00:00 2001 From: Tomas Lycken Date: Fri, 26 Dec 2014 06:09:45 +0100 Subject: [PATCH] Implement a few more BC:s + whitespace and test fixes --- src/Interpolations.jl | 4 ++- src/extrapolation.jl | 11 +++---- src/linear.jl | 5 +-- src/quadratic.jl | 4 +-- test/linear.jl | 2 +- test/quadratic.jl | 23 ++++++++------ test/visual.jl | 74 ++++++++++++------------------------------- 7 files changed, 45 insertions(+), 78 deletions(-) diff --git a/src/Interpolations.jl b/src/Interpolations.jl index 3d741474..299633aa 100644 --- a/src/Interpolations.jl +++ b/src/Interpolations.jl @@ -87,6 +87,7 @@ for IT in ( Linear{OnCell}, Quadratic{Flat,OnCell}, Quadratic{Flat,OnGrid}, + Quadratic{LinearBC,OnGrid}, Quadratic{LinearBC,OnCell} ) for EB in ( @@ -137,7 +138,8 @@ end for IT in ( Quadratic{Flat,OnCell}, Quadratic{Flat,OnGrid}, - Quadratic{LinearBC,OnCell} + Quadratic{LinearBC,OnGrid}, + Quadratic{LinearBC,OnCell}, ) @ngenerate N promote_type_grid(T, x...) function prefilter{T,N}(A::Array{T,N},it::IT) ret, pad = similar_with_padding(A,it) diff --git a/src/extrapolation.jl b/src/extrapolation.jl index dccb5aff..89f602d2 100644 --- a/src/extrapolation.jl +++ b/src/extrapolation.jl @@ -1,6 +1,6 @@ abstract ExtrapolationBehavior -type ExtrapError <: ExtrapolationBehavior end +type ExtrapError <: ExtrapolationBehavior end function extrap_gen(::OnGrid, ::ExtrapError, N) quote @nexprs $N d->(1 <= x_d <= size(itp,d) || throw(BoundsError())) @@ -8,12 +8,11 @@ function extrap_gen(::OnGrid, ::ExtrapError, N) end function extrap_gen(::OnCell, ::ExtrapError, N) quote - @nexprs $N d->(.5 <= x_d <= size(itp,d)+.5 || throw(BoundsError())) + @nexprs $N d->(.5 <= x_d <= size(itp,d) + .5 || throw(BoundsError())) end end type ExtrapNaN <: ExtrapolationBehavior end - function extrap_gen(::OnGrid, ::ExtrapNaN, N) quote @nexprs $N d->(1 <= x_d <= size(itp,d) || return convert(T, NaN)) @@ -21,7 +20,7 @@ function extrap_gen(::OnGrid, ::ExtrapNaN, N) end function extrap_gen(::OnCell, ::ExtrapNaN, N) quote - @nexprs $N d->(.5 <= x_d <= size(itp,d)+.5 || return convert(T, NaN)) + @nexprs $N d->(.5 <= x_d <= size(itp,d) + .5 || return convert(T, NaN)) end end @@ -89,7 +88,7 @@ function extrap_gen(::OnGrid, ::ExtrapLinear, N) if x_d < 1 fx_d = x_d - convert(typeof(x_d), 1) - k = itp[1] - itp[2] + k = -4*itp.coefs[1] return itp[1] - k * fx_d end if x_d > size(itp, d) @@ -102,4 +101,4 @@ function extrap_gen(::OnGrid, ::ExtrapLinear, N) end end end -#extrap_gen(::OnCell, e::ExtrapLinear, N) = extrap_gen(OnGrid(), e, N) +extrap_gen(::OnCell, e::ExtrapLinear, N) = extrap_gen(OnGrid(), e, N) diff --git a/src/linear.jl b/src/linear.jl index 102b5874..f7a78655 100644 --- a/src/linear.jl +++ b/src/linear.jl @@ -10,10 +10,7 @@ end function bc_gen(::Linear{OnCell}, N) quote @nexprs $N d->begin - ix_d = ifloor(x_d) - if .5 <= x_d < 1 - ix_d = 1 - end + ix_d = clamp(ifloor(x_d), 1, size(itp, d)) end end end diff --git a/src/quadratic.jl b/src/quadratic.jl index 48a7c845..773d5e5e 100644 --- a/src/quadratic.jl +++ b/src/quadratic.jl @@ -25,9 +25,9 @@ end function coefficients(::Quadratic, N) quote @nexprs $N d->begin - cm_d = (fx_d-.5)^2 / 2 + cm_d = .5 * (fx_d-.5)^2 c_d = .75 - fx_d^2 - cp_d = (fx_d+.5)^2 / 2 + cp_d = .5 * (fx_d+.5)^2 end end end diff --git a/test/linear.jl b/test/linear.jl index 67732ee3..503e20b2 100644 --- a/test/linear.jl +++ b/test/linear.jl @@ -52,7 +52,7 @@ f(x,y) = sin(x/10)*cos(y/6) xmax, ymax = 30,10 -A = [f(x,y) for x in 1:xmax, y in 1:ymax] +A = Float64[f(x,y) for x in 1:xmax, y in 1:ymax] itp1 = Interpolation(A, Linear(OnGrid()), ExtrapError()) diff --git a/test/quadratic.jl b/test/quadratic.jl index 069cf5c4..96f17fdd 100644 --- a/test/quadratic.jl +++ b/test/quadratic.jl @@ -10,7 +10,7 @@ A = Float64[f(x) for x in 1:xmax] ## ExtrapError -itp1 = Interpolation(A, Quadratic(ExtendInner(),OnCell()), ExtrapError()) +itp1 = Interpolation(A, Quadratic(Flat(),OnCell()), ExtrapError()) for x in [3.1:.2:4.3] @test_approx_eq_eps f(x) itp1[x] abs(.1*f(x)) @@ -21,19 +21,22 @@ end ## ExtrapNaN -itp2 = Interpolation(A, Quadratic(ExtendInner(),OnCell()), ExtrapNaN()) +itp2 = Interpolation(A, Quadratic(Flat(),OnCell()), ExtrapNaN()) -for x in [3.1:.2:4.3] - @test_approx_eq_eps f(x) itp1[x] abs(.1*f(x)) +for x in [3.1:.2:4.3] + @test_approx_eq_eps f(x) itp2[x] abs(.1*f(x)) end -xlo, xhi = itp2[.9], itp2[xmax+.2] -@test isnan(xlo) -@test isnan(xhi) +@test isnan(itp2[.4]) +@test !isnan(itp2[.5]) +@test !isnan(itp2[.6]) +@test !isnan(itp2[xmax+.4]) +@test !isnan(itp2[xmax+.5]) +@test isnan(itp2[xmax+.6]) # Flat/ExtrapConstant -itp3 = Interpolation(A, Quadratic(Flat(),OnCell()), ExtrapConstant()) +itp3 = Interpolation(A, Quadratic(Flat(),OnGrid()), ExtrapConstant()) # Check inbounds and extrap values @@ -42,8 +45,8 @@ for x in [3.1:.2:4.3] end xlo, xhi = itp3[.9], itp3[xmax+.2] -@test xlo == A[1] -@test xhi == A[end] +@test_approx_eq xlo A[1] +@test_approx_eq xhi A[end] # Check continuity xs = [0:.1:length(A)+1] diff --git a/test/visual.jl b/test/visual.jl index 07d97a4d..c847ba17 100644 --- a/test/visual.jl +++ b/test/visual.jl @@ -1,69 +1,35 @@ module VisualTests +print("loading packages...") using Gadfly, Interpolations - +println("done!") const nx = 10 xg = 1:nx -xf = -2nx:.1:3nx +xf = -.2nx:.1:4.3nx f1(x) = sin(2pi/nx * (x-1)) y1 = map(f1,xg) -for EB in ( - ExtrapNaN, - ExtrapConstant, - ExtrapReflect +for (IT, EB) in ( + (Constant{OnCell}, ExtrapConstant), + (Constant{OnGrid}, ExtrapConstant), + (Linear{OnCell}, ExtrapNaN), + (Linear{OnGrid}, ExtrapNaN), + (Quadratic{Flat,OnCell}, ExtrapConstant), + (Quadratic{Flat,OnGrid}, ExtrapConstant), + (Quadratic{LinearBC,OnCell}, ExtrapLinear), + (Quadratic{LinearBC,OnGrid}, ExtrapLinear), + (Quadratic{LinearBC,OnGrid}, ExtrapReflect), + (Quadratic{Flat,OnCell}, ExtrapReflect) ) + itp = Interpolation(y1, IT(), EB()) - for IT in ( - Constant{OnCell}, - Constant{OnGrid}, - Linear{OnCell}, - Linear{OnGrid}, - Quadratic{ExtendInner,OnCell}, - Quadratic{ExtendInner,OnGrid}, - Quadratic{Flat,OnCell}, - Quadratic{Flat,OnGrid}, - ) - - itp = Interpolation(y1, IT(), EB()) - - display(plot( - layer(x=xg,y=y1,Geom.point), - layer(x=xf,y=[itp[x] for x in xf],Geom.path), - Guide.title("$(typeof(Interpolations.degree(IT()))), $(typeof(Interpolations.gridrepresentation(IT()))), $EB") - )) - end + display(plot( + layer(x=xg,y=y1,Geom.point), + layer(x=xf,y=[itp[x] for x in xf],Geom.path), + Guide.title("$IT, $EB") + )) end -for IT in ( - Constant{OnCell}, - Constant{OnGrid}, - Linear{OnCell}, - Linear{OnGrid}, - Quadratic{ExtendInner,OnCell}, - Quadratic{ExtendInner,OnGrid}, - Quadratic{Flat,OnCell}, - Quadratic{Flat,OnGrid}, - ) - - # Treat ExtrapError specially, since it will throw bounds - # errors for any x outside 1:nx - itperr = Interpolation(y1, IT(), ExtrapError()) - if Interpolations.gridrepresentation(IT()) == OnCell() - display(plot( - layer(x=xg,y=y1,Geom.point), - layer(x=.5:.1:nx+.5,y=[itperr[x] for x in .5:.1:nx+.5],Geom.path), - Guide.title("$(Interpolations.degree(IT())), $(Interpolations.gridrepresentation(IT())), ExtrapError") - )) - else - display(plot( - layer(x=xg,y=y1,Geom.point), - layer(x=1:.1:nx,y=[itperr[x] for x in 1:.1:nx],Geom.path), - Guide.title("$(Interpolations.degree(IT())), $(Interpolations.gridrepresentation(IT())), ExtrapError") - )) - end end - -end \ No newline at end of file