Merge 1e33cb2 into 21c1b0e

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

src/scaling/scaling.jl

@@ -126,6 +126,28 @@ 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 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
+
+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
 
 struct ScaledIterator{SITPT,CI,WIS}

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)
 
 

test/scaling/nointerp.jl

@@ -20,6 +20,38 @@ 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)
+
+    for x in (xs[6:end-5] .+ 0.05), 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)
+        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
+
+    @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 +66,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

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
@@ -44,6 +44,19 @@ using Test, LinearAlgebra
         @test ≈(cos(x),g,atol=0.05)
     end
 
+    # Test Hessians of scaled grids
+    xs = -pi:.1:pi
+    ys = -pi:.2:pi
+    zs = sin.(xs) .* sin.(ys')
+    itp = interpolate(zs, BSpline(Cubic(Line(OnGrid()))))
+    sitp = @inferred scale(itp, xs, ys)
+
+    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)
+    end
+
     # Verify that return types are reasonable
     @inferred(sitp2(-3.4, 1.2))
     @inferred(sitp2(-3, 1))