From 3c69f88982068f5f8ef7cdaa86c2e6da54aaceb9 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 16 Nov 2018 10:21:26 +0100 Subject: [PATCH 1/8] Add Hessians for ScaledInterpolation and tests --- src/scaling/scaling.jl | 28 ++++++++++++++++++++++++++++ test/scaling/scaling.jl | 14 ++++++++++++++ 2 files changed, 42 insertions(+) diff --git a/src/scaling/scaling.jl b/src/scaling/scaling.jl index dca0adec..c3dda755 100644 --- a/src/scaling/scaling.jl +++ b/src/scaling/scaling.jl @@ -126,6 +126,34 @@ rescale_gradient(r::UnitRange, g) = g Implements the chain rule dy/dx = dy/du * du/dx for use when calculating gradients with scaled interpolation objects. """ rescale_gradient +@propagate_inbounds function hessian(sitp::ScaledInterpolation{T,N}, xs::Vararg{Number,N}) where {T,N} + @boundscheck (checkbounds(Bool, sitp, xs...) || Base.throw_boundserror(sitp, xs)) + xl = maybe_clamp(sitp.itp, coordslookup(itpflag(sitp.itp), sitp.ranges, xs)) + h = hessian(sitp.itp, xl...) + return symmatrix(rescale_hessian_components(itpflag(sitp.itp), sitp.ranges, Tuple(h), size(h, 1))) +end + +function rescale_hessian_components(flags, ranges, h, n) + hs = () + for i=1:n + if getfirst(flags) isa NoInterp + flags = getrest(flags) + ranges = Base.tail(ranges) + else + s1 = rescale_gradient_components(flags, ranges, Tuple(h[i*(n+1)-n:i*n])) + s2 = rescale_hessian(ranges[1], s1) + hs = (hs..., s2...) + flags = getrest(flags) + ranges = Base.tail(ranges) + end + end + return hs +end + +rescale_hessian(r::StepRangeLen, h) = h ./ step(r) +rescale_hessian(r::StepRange, h) = h ./ r.step +rescale_hessian(r::UnitRange, h) = h + ### Iteration struct ScaledIterator{SITPT,CI,WIS} diff --git a/test/scaling/scaling.jl b/test/scaling/scaling.jl index 62851316..042c8e4e 100644 --- a/test/scaling/scaling.jl +++ b/test/scaling/scaling.jl @@ -44,6 +44,20 @@ using Test, LinearAlgebra @test ≈(cos(x),g,atol=0.05) end + # Test Hessians of scaled grids + xs = -pi:.1:pi + ys = -pi:.1:pi + zs = sin.(xs') .* sin.(ys) + itp = interpolate(zs, BSpline(Cubic(Line(OnGrid())))) + sitp = @inferred scale(itp, xs, ys) + hitp = (x,y) -> Interpolations.hessian(sitp, x, y) + + for x in xs[2:end-1], y in ys[2:end-1] + # h = @inferred(Interpolations.hessian(sitp, x, y))[1] + h = Interpolations.hessian(sitp, x, y) + @test ≈([-sin(x) * sin(y) cos(x) * cos(y); cos(x) * cos(y) -sin(x) * sin(y)], h, atol=0.03) + end + # Verify that return types are reasonable @inferred(sitp2(-3.4, 1.2)) @inferred(sitp2(-3, 1)) From 8fb17fc1ddde3495104f16cc16fd07ae128fbf88 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 16 Nov 2018 14:47:10 +0100 Subject: [PATCH 2/8] hessian with NoInterp axes * fixed rescale_hessian * added tests --- src/scaling/scaling.jl | 6 ++++-- test/scaling/nointerp.jl | 30 ++++++++++++++++++++++++++++++ 2 files changed, 34 insertions(+), 2 deletions(-) diff --git a/src/scaling/scaling.jl b/src/scaling/scaling.jl index c3dda755..6fc3b51f 100644 --- a/src/scaling/scaling.jl +++ b/src/scaling/scaling.jl @@ -135,16 +135,18 @@ end function rescale_hessian_components(flags, ranges, h, n) hs = () - for i=1:n + idx = 1 + while idx <= n if getfirst(flags) isa NoInterp flags = getrest(flags) ranges = Base.tail(ranges) else - s1 = rescale_gradient_components(flags, ranges, Tuple(h[i*(n+1)-n:i*n])) + s1 = rescale_gradient_components(flags, ranges, Tuple(h[(idx-1)*n+idx:idx*n])) s2 = rescale_hessian(ranges[1], s1) hs = (hs..., s2...) flags = getrest(flags) ranges = Base.tail(ranges) + idx += 1 end end return hs diff --git a/test/scaling/nointerp.jl b/test/scaling/nointerp.jl index 2128990b..8529a685 100644 --- a/test/scaling/nointerp.jl +++ b/test/scaling/nointerp.jl @@ -20,6 +20,35 @@ using Interpolations, Test, LinearAlgebra, Random @test length(Interpolations.gradient(sitp, pi/3, 2)) == 1 + xs = range(-pi, stop=pi, length=60) + A = hcat(map(f1, xs), map(f2, xs), map(f3, xs)) + itp = interpolate(A, (BSpline(Cubic(Periodic(OnGrid()))), NoInterp())) + sitp = scale(itp, xs, ys) + h(x, y) = Interpolations.hessian(sitp, x, y) + h(xs[10], 1) + for x in xs[6:end-6], y in ys + if y in (1,2) + @test ≈(h(x,y)[1], -f(x,y), atol=0.05) + elseif y==3 + @test ≈(h(x,y)[1], -4*f(x,y), atol=0.05) + end + end + + @test length(Interpolations.hessian(sitp, pi/3, 2)) == 1 + + A = hcat(map(f1, xs), map(f2, xs), map(f3, xs)) + itp = interpolate(A', (NoInterp(), BSpline(Cubic(Periodic(OnGrid()))))) + sitp = scale(itp, ys, xs) + h(y, x) = Interpolations.hessian(sitp, y, x) + h(1, xs[10]) + for x in xs[6:end-6], y in ys + if y in (1,2) + @test ≈(h(y,x)[1], -f(x,y), atol=0.05) + elseif y==3 + @test ≈(h(y,x)[1], -4*f(x,y), atol=0.05) + end + end + # check for case where initial/middle indices are NoInterp but later ones are <:BSpline isdefined(Random, :seed!) ? Random.seed!(1234) : srand(1234) # `srand` was renamed to `seed!` z0 = rand(10,10) @@ -34,4 +63,5 @@ using Interpolations, Test, LinearAlgebra, Random sitpb = scale(itpb, 1:10, rng) @test Interpolations.gradient(sitpa, 3.0, 3) == Interpolations.gradient(sitpb, 3, 3.0) + end From 74ac76d5b113276ae437ecc6ebb801fa4add144f Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 16 Nov 2018 23:18:53 +0100 Subject: [PATCH 3/8] fix inference for hessian --- src/scaling/scaling.jl | 38 +++++++++++++++----------------------- test/scaling/scaling.jl | 4 ++-- 2 files changed, 17 insertions(+), 25 deletions(-) diff --git a/src/scaling/scaling.jl b/src/scaling/scaling.jl index 6fc3b51f..6927c88a 100644 --- a/src/scaling/scaling.jl +++ b/src/scaling/scaling.jl @@ -130,31 +130,23 @@ Implements the chain rule dy/dx = dy/du * du/dx for use when calculating gradien @boundscheck (checkbounds(Bool, sitp, xs...) || Base.throw_boundserror(sitp, xs)) xl = maybe_clamp(sitp.itp, coordslookup(itpflag(sitp.itp), sitp.ranges, xs)) h = hessian(sitp.itp, xl...) - return symmatrix(rescale_hessian_components(itpflag(sitp.itp), sitp.ranges, Tuple(h), size(h, 1))) -end - -function rescale_hessian_components(flags, ranges, h, n) - hs = () - idx = 1 - while idx <= n - if getfirst(flags) isa NoInterp - flags = getrest(flags) - ranges = Base.tail(ranges) - else - s1 = rescale_gradient_components(flags, ranges, Tuple(h[(idx-1)*n+idx:idx*n])) - s2 = rescale_hessian(ranges[1], s1) - hs = (hs..., s2...) - flags = getrest(flags) - ranges = Base.tail(ranges) - idx += 1 - end - end - return hs + return rescale_hessian_components(itpflag(sitp.itp), sitp.ranges, h) +end + +function rescale_hessian_components(flags, ranges, h) + steps = SVector(get_steps(flags, ranges)) + return h ./ (steps .* steps') end -rescale_hessian(r::StepRangeLen, h) = h ./ step(r) -rescale_hessian(r::StepRange, h) = h ./ r.step -rescale_hessian(r::UnitRange, h) = h +function get_steps(flags, ranges) + if getfirst(flags) isa NoInterp + return get_steps(getrest(flags), Base.tail(ranges)) + else + item = step(ranges[1]) + return (item, get_steps(getrest(flags), Base.tail(ranges))...) + end +end +get_steps(flags, ::Tuple{}) = () ### Iteration diff --git a/test/scaling/scaling.jl b/test/scaling/scaling.jl index 042c8e4e..0e9b9401 100644 --- a/test/scaling/scaling.jl +++ b/test/scaling/scaling.jl @@ -53,8 +53,8 @@ using Test, LinearAlgebra hitp = (x,y) -> Interpolations.hessian(sitp, x, y) for x in xs[2:end-1], y in ys[2:end-1] - # h = @inferred(Interpolations.hessian(sitp, x, y))[1] - h = Interpolations.hessian(sitp, x, y) + h = @inferred(Interpolations.hessian(sitp, x, y)) + @test issymmetric(h) @test ≈([-sin(x) * sin(y) cos(x) * cos(y); cos(x) * cos(y) -sin(x) * sin(y)], h, atol=0.03) end From 1fb79076781bcf9cc83cb12f358b413c05d0129a Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 16 Nov 2018 23:24:12 +0100 Subject: [PATCH 4/8] make hessian test less trivial --- test/scaling/scaling.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/scaling/scaling.jl b/test/scaling/scaling.jl index 0e9b9401..61c6678a 100644 --- a/test/scaling/scaling.jl +++ b/test/scaling/scaling.jl @@ -46,8 +46,8 @@ using Test, LinearAlgebra # Test Hessians of scaled grids xs = -pi:.1:pi - ys = -pi:.1:pi - zs = sin.(xs') .* sin.(ys) + ys = -pi:.2:pi + zs = sin.(xs) .* sin.(ys') itp = interpolate(zs, BSpline(Cubic(Line(OnGrid())))) sitp = @inferred scale(itp, xs, ys) hitp = (x,y) -> Interpolations.hessian(sitp, x, y) From 35cef09ddfc7158e5d9ca0084100dc542b879a49 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 17 Nov 2018 18:22:52 +0100 Subject: [PATCH 5/8] added and cleaned up tests --- test/hessian.jl | 1 + test/scaling/nointerp.jl | 15 +++++++++------ test/scaling/scaling.jl | 3 +-- 3 files changed, 11 insertions(+), 8 deletions(-) diff --git a/test/hessian.jl b/test/hessian.jl index 297f6965..cfeb9a35 100644 --- a/test/hessian.jl +++ b/test/hessian.jl @@ -41,6 +41,7 @@ using Test, Interpolations, LinearAlgebra, ForwardDiff itp = interpolate(A2, (BSpline(Quadratic(Flat(OnCell()))), NoInterp())) v = A2[:, 2] itpcol = interpolate(v, BSpline(Quadratic(Flat(OnCell())))) + @inferred(Interpolations.hessian(itp, 3.2, 2)) @test Interpolations.hessian(itp, 3.2, 2) == Interpolations.hessian(itpcol, 3.2) diff --git a/test/scaling/nointerp.jl b/test/scaling/nointerp.jl index 8529a685..9ddf62da 100644 --- a/test/scaling/nointerp.jl +++ b/test/scaling/nointerp.jl @@ -24,13 +24,16 @@ using Interpolations, Test, LinearAlgebra, Random A = hcat(map(f1, xs), map(f2, xs), map(f3, xs)) itp = interpolate(A, (BSpline(Cubic(Periodic(OnGrid()))), NoInterp())) sitp = scale(itp, xs, ys) - h(x, y) = Interpolations.hessian(sitp, x, y) - h(xs[10], 1) - for x in xs[6:end-6], y in ys + + for x in (xs[6:end-5] .+ 0.05), y in ys if y in (1,2) - @test ≈(h(x,y)[1], -f(x,y), atol=0.05) - elseif y==3 - @test ≈(h(x,y)[1], -4*f(x,y), atol=0.05) + @test_broken h = @inferred(Interpolations.hessian(sitp, x, y)) + h = Interpolations.hessian(sitp, x, y) + @test ≈(h[1], -f(x,y), atol=0.05) + else # y==3 + @test_broken h = @inferred(Interpolations.hessian(sitp, x, y)) + h = Interpolations.hessian(sitp, x, y) + @test ≈(h[1], -4*f(x,y), atol=0.05) end end diff --git a/test/scaling/scaling.jl b/test/scaling/scaling.jl index 61c6678a..28eaf556 100644 --- a/test/scaling/scaling.jl +++ b/test/scaling/scaling.jl @@ -1,5 +1,5 @@ using Interpolations -using Test, LinearAlgebra +using Test, LinearAlgebra, StaticArrays @testset "Scaling" begin # Model linear interpolation of y = -3 + .5x by interpolating y=x @@ -50,7 +50,6 @@ using Test, LinearAlgebra zs = sin.(xs) .* sin.(ys') itp = interpolate(zs, BSpline(Cubic(Line(OnGrid())))) sitp = @inferred scale(itp, xs, ys) - hitp = (x,y) -> Interpolations.hessian(sitp, x, y) for x in xs[2:end-1], y in ys[2:end-1] h = @inferred(Interpolations.hessian(sitp, x, y)) From 1e33cb20f94604e6a9258502cab643959bfd9a2c Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 17 Nov 2018 18:24:00 +0100 Subject: [PATCH 6/8] added a dummy hessian(::Extrapolation) works if the point is inbounds, throws an informative error otherwise --- src/extrapolation/extrapolation.jl | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/extrapolation/extrapolation.jl b/src/extrapolation/extrapolation.jl index 15aabe4b..17a1f337 100644 --- a/src/extrapolation/extrapolation.jl +++ b/src/extrapolation/extrapolation.jl @@ -67,6 +67,25 @@ end end end +@inline function hessian(etp::AbstractExtrapolation{T,N}, x::Vararg{Number,N}) where {T,N} + itp = parent(etp) + if checkbounds(Bool, itp, x...) + hessian(itp, x...) + else + error("extrapolation of hessian not yet implemented") + # # copied from gradient above, with obvious modifications + # # but final part is missing + # eflag = tcollect(etpflag, etp) + # xs = inbounds_position(eflag, bounds(itp), x, etp, x) + # h = @inbounds hessian(itp, xs...) + # skipni = t->skip_flagged_nointerp(itp, t) + # # not sure if it should be just h here instead of Tuple(h) + # # extrapolate_hessian needs to be written + # # SVector is likely also wrong here + # SVector(extrapolate_hessian.(skipni(eflag), skipni(x), skipni(xs), Tuple(h))) + end +end + checkbounds(::Bool, ::AbstractExtrapolation, I...) = true # The last two arguments are just for error-reporting From 09fa5c51e4ac76afc83ad094aebe11e918661a12 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sun, 18 Nov 2018 16:50:28 +0100 Subject: [PATCH 7/8] use infix approx in tests --- test/scaling/nointerp.jl | 12 ++++++------ test/scaling/scaling.jl | 5 +++-- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/test/scaling/nointerp.jl b/test/scaling/nointerp.jl index 9ddf62da..7d906d22 100644 --- a/test/scaling/nointerp.jl +++ b/test/scaling/nointerp.jl @@ -4,7 +4,7 @@ using Interpolations, Test, LinearAlgebra, Random xs = -pi:2pi/10:pi f1(x) = sin(x) f2(x) = cos(x) - f3(x) = sin(x) .* cos(x) + f3(x) = sin(x) * cos(x) 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)"))) ys = 1:3 @@ -15,7 +15,7 @@ using Interpolations, Test, LinearAlgebra, Random for (ix,x0) in enumerate(xs[1:end-1]), y0 in ys x,y = x0, y0 - @test ≈(sitp(x,y),f(x,y),atol=0.05) + @test sitp(x,y) ≈ f(x,y) atol=0.05 end @test length(Interpolations.gradient(sitp, pi/3, 2)) == 1 @@ -29,11 +29,11 @@ using Interpolations, Test, LinearAlgebra, Random if y in (1,2) @test_broken h = @inferred(Interpolations.hessian(sitp, x, y)) h = Interpolations.hessian(sitp, x, y) - @test ≈(h[1], -f(x,y), atol=0.05) + @test h[1] ≈ -f(x, y) atol=0.05 else # y==3 @test_broken h = @inferred(Interpolations.hessian(sitp, x, y)) h = Interpolations.hessian(sitp, x, y) - @test ≈(h[1], -4*f(x,y), atol=0.05) + @test h[1] ≈ -4*f(x, y) atol=0.05 end end @@ -46,9 +46,9 @@ using Interpolations, Test, LinearAlgebra, Random h(1, xs[10]) for x in xs[6:end-6], y in ys if y in (1,2) - @test ≈(h(y,x)[1], -f(x,y), atol=0.05) + @test h(y, x)[1] ≈ -f(x, y) atol=0.05 elseif y==3 - @test ≈(h(y,x)[1], -4*f(x,y), atol=0.05) + @test h(y, x)[1] ≈ -4*f(x, y) atol=0.05 end end diff --git a/test/scaling/scaling.jl b/test/scaling/scaling.jl index 28eaf556..f9e42e7f 100644 --- a/test/scaling/scaling.jl +++ b/test/scaling/scaling.jl @@ -41,7 +41,7 @@ using Test, LinearAlgebra, StaticArrays for x in -pi:.1:pi g = @inferred(Interpolations.gradient(sitp, x))[1] - @test ≈(cos(x),g,atol=0.05) + @test cos(x) ≈ g atol=0.05 end # Test Hessians of scaled grids @@ -54,7 +54,8 @@ using Test, LinearAlgebra, StaticArrays for x in xs[2:end-1], y in ys[2:end-1] h = @inferred(Interpolations.hessian(sitp, x, y)) @test issymmetric(h) - @test ≈([-sin(x) * sin(y) cos(x) * cos(y); cos(x) * cos(y) -sin(x) * sin(y)], h, atol=0.03) + @test [-sin(x) * sin(y) cos(x) * cos(y) + cos(x) * cos(y) -sin(x) * sin(y)] ≈ h atol=0.03 end # Verify that return types are reasonable From 42f8008edee3b236217edf4a6145bc6c40870f4e Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Thu, 22 Nov 2018 13:35:24 +0100 Subject: [PATCH 8/8] include inference tests for nointerp --- test/scaling/nointerp.jl | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/test/scaling/nointerp.jl b/test/scaling/nointerp.jl index 7d906d22..dbc716c4 100644 --- a/test/scaling/nointerp.jl +++ b/test/scaling/nointerp.jl @@ -20,20 +20,18 @@ using Interpolations, Test, LinearAlgebra, Random @test length(Interpolations.gradient(sitp, pi/3, 2)) == 1 - xs = range(-pi, stop=pi, length=60) + xs = range(-pi, stop=pi, length=60)[1:end-1] A = hcat(map(f1, xs), map(f2, xs), map(f3, xs)) itp = interpolate(A, (BSpline(Cubic(Periodic(OnGrid()))), NoInterp())) sitp = scale(itp, xs, ys) - for x in (xs[6:end-5] .+ 0.05), y in ys + for x in xs, y in ys if y in (1,2) - @test_broken h = @inferred(Interpolations.hessian(sitp, x, y)) - h = Interpolations.hessian(sitp, x, y) - @test h[1] ≈ -f(x, y) atol=0.05 + h = @inferred(Interpolations.hessian(sitp, x, y)) + @test h[1] ≈ -f(x, y) atol=0.01 else # y==3 - @test_broken h = @inferred(Interpolations.hessian(sitp, x, y)) - h = Interpolations.hessian(sitp, x, y) - @test h[1] ≈ -4*f(x, y) atol=0.05 + h = @inferred(Interpolations.hessian(sitp, x, y)) + @test h[1] ≈ -4*f(x, y) atol=0.01 end end