Skip to content

Commit aaf9ca0

Browse files
committed
Re-enable scaling with the exception of iteration
1 parent 53a7bc0 commit aaf9ca0

File tree

11 files changed

+401
-418
lines changed

11 files changed

+401
-418
lines changed

src/Interpolations.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -216,7 +216,7 @@ include("nointerp/nointerp.jl")
216216
include("b-splines/b-splines.jl")
217217
# include("gridded/gridded.jl")
218218
include("extrapolation/extrapolation.jl")
219-
# include("scaling/scaling.jl")
219+
include("scaling/scaling.jl")
220220
include("utils.jl")
221221
include("io.jl")
222222
include("convenience-constructors.jl")

src/b-splines/b-splines.jl

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -59,10 +59,10 @@ _lbounds(::Tuple{}, itp, gt) = ()
5959
_ubounds(::Tuple{}, itp, gt) = ()
6060

6161
# The unpadded defaults
62-
lbound(ax, ::BSpline, ::OnCell) = first(ax) - 0.5
63-
ubound(ax, ::BSpline, ::OnCell) = last(ax) + 0.5
64-
lbound(ax, ::BSpline, ::OnGrid) = first(ax)
65-
ubound(ax, ::BSpline, ::OnGrid) = last(ax)
62+
lbound(ax::AbstractUnitRange, ::BSpline, ::OnCell) = first(ax) - 0.5
63+
ubound(ax::AbstractUnitRange, ::BSpline, ::OnCell) = last(ax) + 0.5
64+
lbound(ax::AbstractUnitRange, ::BSpline, ::OnGrid) = first(ax)
65+
ubound(ax::AbstractUnitRange, ::BSpline, ::OnGrid) = last(ax)
6666

6767
fix_axis(r::Base.OneTo) = r
6868
fix_axis(r::Base.Slice) = r

src/extrapolation/extrapolation.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@ struct Extrapolation{T,N,ITPT,IT,GT,ET} <: AbstractExtrapolation{T,N,ITPT,IT,GT}
44
end
55

66
Base.parent(A::Extrapolation) = A.itp
7+
itpflag(etp::Extrapolation) = itpflag(etp.itp)
78

89
# DimSpec{Flag} is not enough for extrapolation dispatch, since we allow nested tuples
910
# However, no tuples should be nested deeper than this; the first level is for different

src/extrapolation/filled.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ end
1010

1111
Base.parent(A::FilledExtrapolation) = A.itp
1212
etpflag(A::FilledExtrapolation) = A.fillvalue
13+
itpflag(A::FilledExtrapolation) = itpflag(A.itp)
1314

1415
"""
1516
`extrapolate(itp, fillvalue)` creates an extrapolation object that returns the `fillvalue` any time the indexes in `itp(x1,x2,...)` are out-of-bounds.

src/scaling/scaling.jl

Lines changed: 245 additions & 254 deletions
Large diffs are not rendered by default.

test/runtests.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -18,8 +18,8 @@ using Interpolations
1818
# extrapolation tests
1919
include("extrapolation/runtests.jl")
2020

21-
# # scaling tests
22-
# include("scaling/runtests.jl")
21+
# scaling tests
22+
include("scaling/runtests.jl")
2323

2424
# # test gradient evaluation
2525
include("gradient.jl")

test/scaling/dimspecs.jl

Lines changed: 14 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -1,23 +1,21 @@
1-
module ScalingDimspecTests
2-
31
using Interpolations, DualNumbers, Test, LinearAlgebra
42

5-
xs = -pi:(2pi/10):pi-2pi/10
6-
ys = -2:.1:2
7-
f(x,y) = sin(x) * y^2
8-
9-
itp = interpolate(Float64[f(x,y) for x in xs, y in ys], (BSpline(Quadratic(Periodic())), BSpline(Linear())), OnGrid())
10-
sitp = scale(itp, xs, ys)
3+
@testset "ScalingDimspecTests" begin
4+
xs = -pi:(2pi/10):pi-2pi/10
5+
ys = -2:.1:2
6+
f(x,y) = sin(x) * y^2
117

12-
for (ix,x) in enumerate(xs), (iy,y) in enumerate(ys)
13-
@test (sitp[x,y],f(x,y),atol=sqrt(eps(1.0)))
8+
itp = interpolate(Float64[f(x,y) for x in xs, y in ys], (BSpline(Quadratic(Periodic())), BSpline(Linear())), OnGrid())
9+
sitp = scale(itp, xs, ys)
1410

15-
g = Interpolations.gradient(sitp, x, y)
16-
fx = epsilon(sitp[dual(x,1), dual(y,0)])
17-
fy = epsilon(sitp[dual(x,0), dual(y,1)])
11+
for (ix,x) in enumerate(xs), (iy,y) in enumerate(ys)
12+
@test (sitp(x,y),f(x,y),atol=sqrt(eps(1.0)))
1813

19-
@test (g[1],fx,atol=sqrt(eps(1.0)))
20-
@test (g[2],fy,atol=sqrt(eps(1.0)))
21-
end
14+
g = Interpolations.gradient(sitp, x, y)
15+
fx = epsilon(sitp(dual(x,1), dual(y,0)))
16+
fy = epsilon(sitp(dual(x,0), dual(y,1)))
2217

18+
@test (g[1],fx,atol=sqrt(eps(1.0)))
19+
@test (g[2],fy,atol=sqrt(eps(1.0)))
20+
end
2321
end

test/scaling/nointerp.jl

Lines changed: 26 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -1,38 +1,37 @@
1-
module ScalingNoInterpTests
2-
31
using Interpolations, Test, LinearAlgebra, Random
42

5-
xs = -pi:2pi/10:pi
6-
f1(x) = sin(x)
7-
f2(x) = cos(x)
8-
f3(x) = sin(x) .* cos(x)
9-
f(x,y) = y == 1 ? f1(x) : (y == 2 ? f2(x) : (y == 3 ? f3(x) : error("invalid value for y (must be 1, 2 or 3, you used $y)")))
10-
ys = 1:3
3+
@testset "ScalingNoInterpTests" begin
4+
xs = -pi:2pi/10:pi
5+
f1(x) = sin(x)
6+
f2(x) = cos(x)
7+
f3(x) = sin(x) .* cos(x)
8+
f(x,y) = y == 1 ? f1(x) : (y == 2 ? f2(x) : (y == 3 ? f3(x) : error("invalid value for y (must be 1, 2 or 3, you used $y)")))
9+
ys = 1:3
1110

12-
A = hcat(map(f1, xs), map(f2, xs), map(f3, xs))
11+
A = hcat(map(f1, xs), map(f2, xs), map(f3, xs))
1312

14-
itp = interpolate(A, (BSpline(Quadratic(Periodic())), NoInterp()), OnGrid())
15-
sitp = scale(itp, xs, ys)
13+
itp = interpolate(A, (BSpline(Quadratic(Periodic())), NoInterp()), OnGrid())
14+
sitp = scale(itp, xs, ys)
1615

17-
for (ix,x0) in enumerate(xs[1:end-1]), y0 in ys
18-
x,y = x0, y0
19-
@test (sitp[x,y],f(x,y),atol=0.05)
20-
end
16+
for (ix,x0) in enumerate(xs[1:end-1]), y0 in ys
17+
x,y = x0, y0
18+
@test (sitp(x,y),f(x,y),atol=0.05)
19+
end
2120

22-
@test length(Interpolations.gradient(sitp, pi/3, 2)) == 1
21+
@test length(Interpolations.gradient(sitp, pi/3, 2)) == 1
2322

24-
# check for case where initial/middle indices are NoInterp but later ones are <:BSpline
25-
isdefined(Random, :seed!) ? Random.seed!(1234) : srand(1234) # `srand` was renamed to `seed!`
26-
z0 = rand(10,10)
27-
za = copy(z0)
28-
zb = copy(z0')
23+
# check for case where initial/middle indices are NoInterp but later ones are <:BSpline
24+
isdefined(Random, :seed!) ? Random.seed!(1234) : srand(1234) # `srand` was renamed to `seed!`
25+
z0 = rand(10,10)
26+
za = copy(z0)
27+
zb = copy(z0')
2928

30-
itpa = interpolate(za, (BSpline(Linear()), NoInterp()), OnGrid())
31-
itpb = interpolate(zb, (NoInterp(), BSpline(Linear())), OnGrid())
29+
itpa = interpolate(za, (BSpline(Linear()), NoInterp()), OnGrid())
30+
itpb = interpolate(zb, (NoInterp(), BSpline(Linear())), OnGrid())
3231

33-
rng = range(1.0, stop=19.0, length=10)
34-
sitpa = scale(itpa, rng, 1:10)
35-
sitpb = scale(itpb, 1:10, rng)
36-
@test Interpolations.gradient(sitpa, 3.0, 3) == Interpolations.gradient(sitpb, 3, 3.0)
32+
rng = range(1.0, stop=19.0, length=10)
33+
sitpa = scale(itpa, rng, 1:10)
34+
sitpb = scale(itpb, 1:10, rng)
35+
@test Interpolations.gradient(sitpa, 3.0, 3) == Interpolations.gradient(sitpb, 3, 3.0)
3736

3837
end

test/scaling/runtests.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,4 +2,4 @@ include("scaling.jl")
22
include("dimspecs.jl")
33
include("nointerp.jl")
44
include("withextrap.jl")
5-
include("function-call-syntax.jl")
5+
# include("function-call-syntax.jl")

test/scaling/scaling.jl

Lines changed: 78 additions & 85 deletions
Original file line numberDiff line numberDiff line change
@@ -1,101 +1,94 @@
1-
module ScalingTests
2-
31
using Interpolations
42
using Test, LinearAlgebra
53

6-
# Model linear interpolation of y = -3 + .5x by interpolating y=x
7-
# and then scaling to the new x range
8-
9-
itp = interpolate(1:1.0:10, BSpline(Linear()), OnGrid())
4+
@testset "Scaling" begin
5+
# Model linear interpolation of y = -3 + .5x by interpolating y=x
6+
# and then scaling to the new x range
107

11-
sitp = @inferred(scale(itp, -3:.5:1.5))
12-
@test typeof(sitp) <: Interpolations.ScaledInterpolation
13-
@test parent(sitp) === itp
14-
15-
for (x,y) in zip(-3:.05:1.5, 1:.1:10)
16-
@test sitp[x] y
17-
end
8+
itp = interpolate(1:1.0:10, BSpline(Linear()), OnGrid())
189

19-
# Verify that it works in >1D, with different types of ranges
10+
sitp = @inferred(scale(itp, -3:.5:1.5))
11+
@test typeof(sitp) <: Interpolations.ScaledInterpolation
12+
@test parent(sitp) === itp
2013

21-
gauss(phi, mu, sigma) = exp(-(phi-mu)^2 / (2sigma)^2)
22-
testfunction(x,y) = gauss(x, 0.5, 4) * gauss(y, -.5, 2)
14+
for (x,y) in zip(-3:.05:1.5, 1:.1:10)
15+
@test sitp(x) y
16+
end
2317

24-
xs = -5:.5:5
25-
ys = -4:.2:4
26-
zs = Float64[testfunction(x,y) for x in xs, y in ys]
18+
# Verify that it works in >1D, with different types of ranges
2719

28-
itp2 = interpolate(zs, BSpline(Quadratic(Flat())), OnGrid())
29-
sitp2 = @inferred scale(itp2, xs, ys)
20+
gauss(phi, mu, sigma) = exp(-(phi-mu)^2 / (2sigma)^2)
21+
testfunction(x,y) = gauss(x, 0.5, 4) * gauss(y, -.5, 2)
3022

31-
for x in xs, y in ys
32-
@test testfunction(x,y) sitp2[x,y]
33-
end
23+
xs = -5:.5:5
24+
ys = -4:.2:4
25+
zs = Float64[testfunction(x,y) for x in xs, y in ys]
3426

35-
# Test gradients of scaled grids
36-
xs = -pi:.1:pi
37-
ys = map(sin, xs)
38-
itp = interpolate(ys, BSpline(Linear()), OnGrid())
39-
sitp = @inferred scale(itp, xs)
27+
itp2 = interpolate(zs, BSpline(Quadratic(Flat())), OnGrid())
28+
sitp2 = @inferred scale(itp2, xs, ys)
4029

41-
for x in -pi:.1:pi
42-
g = @inferred(Interpolations.gradient(sitp, x))[1]
43-
@test (cos(x),g,atol=0.05)
44-
end
30+
for x in xs, y in ys
31+
@test testfunction(x,y) sitp2(x,y)
32+
end
4533

46-
# Verify that return types are reasonable
47-
@inferred(getindex(sitp2, -3.4, 1.2))
48-
@inferred(getindex(sitp2, -3, 1))
49-
@inferred(getindex(sitp2, -3.4, 1))
50-
51-
sitp32 = @inferred 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)
52-
@test typeof(@inferred(getindex(sitp32, -3.4f0, 1.2f0))) == Float32
53-
54-
# Iteration
55-
itp = interpolate(rand(3,3,3), BSpline(Quadratic(Flat())), OnCell())
56-
knots = map(d->1:10:21, 1:3)
57-
sitp = @inferred scale(itp, knots...)
58-
59-
iter = @inferred(eachvalue(sitp))
60-
61-
@static if VERSION < v"0.7.0-DEV.5126"
62-
state = @inferred(start(iter))
63-
@test !(@inferred(done(iter, state)))
64-
val, state = @inferred(next(iter, state))
65-
else
66-
iter_next = iterate(iter)
67-
@test iter_next isa Tuple
68-
@test iter_next[1] isa Float64
69-
state = iter_next[2]
70-
inferred_next = Base.return_types(iterate, (typeof(iter),))
71-
@test length(inferred_next) == 1
72-
@test inferred_next[1] == Union{Nothing,Tuple{Float64,typeof(state)}}
73-
iter_next = iterate(iter, state)
74-
@test iter_next isa Tuple
75-
@test iter_next[1] isa Float64
76-
inferred_next = Base.return_types(iterate, (typeof(iter),typeof(state)))
77-
state = iter_next[2]
78-
@test length(inferred_next) == 1
79-
@test inferred_next[1] == Union{Nothing,Tuple{Float64,typeof(state)}}
80-
end
34+
# Test gradients of scaled grids
35+
xs = -pi:.1:pi
36+
ys = map(sin, xs)
37+
itp = interpolate(ys, BSpline(Linear()), OnGrid())
38+
sitp = @inferred scale(itp, xs)
8139

82-
function foo!(dest, sitp)
83-
i = 0
84-
for s in eachvalue(sitp)
85-
dest[i+=1] = s
86-
end
87-
dest
88-
end
89-
function bar!(dest, sitp)
90-
for I in CartesianIndices(size(dest))
91-
dest[I] = sitp[I]
40+
for x in -pi:.1:pi
41+
g = @inferred(Interpolations.gradient(sitp, x))[1]
42+
@test (cos(x),g,atol=0.05)
9243
end
93-
dest
94-
end
95-
rfoo = Array{Float64}(undef, Interpolations.ssize(sitp))
96-
rbar = similar(rfoo)
97-
foo!(rfoo, sitp)
98-
bar!(rbar, sitp)
99-
@test rfoo rbar
44+
45+
# Verify that return types are reasonable
46+
@inferred(sitp2(-3.4, 1.2))
47+
@inferred(sitp2(-3, 1))
48+
@inferred(sitp2(-3.4, 1))
49+
50+
sitp32 = @inferred 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)
51+
@test typeof(@inferred(sitp32(-3.4f0, 1.2f0))) == Float32
52+
53+
# # Iteration
54+
# itp = interpolate(rand(3,3,3), BSpline(Quadratic(Flat())), OnCell())
55+
# knots = map(d->1:10:21, 1:3)
56+
# sitp = @inferred scale(itp, knots...)
57+
58+
# iter = @inferred(eachvalue(sitp))
59+
60+
# iter_next = iterate(iter)
61+
# @test iter_next isa Tuple
62+
# @test iter_next[1] isa Float64
63+
# state = iter_next[2]
64+
# inferred_next = Base.return_types(iterate, (typeof(iter),))
65+
# @test length(inferred_next) == 1
66+
# @test inferred_next[1] == Union{Nothing,Tuple{Float64,typeof(state)}}
67+
# iter_next = iterate(iter, state)
68+
# @test iter_next isa Tuple
69+
# @test iter_next[1] isa Float64
70+
# inferred_next = Base.return_types(iterate, (typeof(iter),typeof(state)))
71+
# state = iter_next[2]
72+
# @test length(inferred_next) == 1
73+
# @test inferred_next[1] == Union{Nothing,Tuple{Float64,typeof(state)}}
74+
75+
# function foo!(dest, sitp)
76+
# i = 0
77+
# for s in eachvalue(sitp)
78+
# dest[i+=1] = s
79+
# end
80+
# dest
81+
# end
82+
# function bar!(dest, sitp)
83+
# for I in CartesianIndices(size(dest))
84+
# dest[I] = sitp[I]
85+
# end
86+
# dest
87+
# end
88+
# rfoo = Array{Float64}(undef, Interpolations.ssize(sitp))
89+
# rbar = similar(rfoo)
90+
# foo!(rfoo, sitp)
91+
# bar!(rbar, sitp)
92+
# @test rfoo ≈ rbar
10093

10194
end

0 commit comments

Comments
 (0)