Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Monotonic interpolations #238

Closed
wants to merge 12 commits into from
Closed

Monotonic interpolations #238

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
c3d3d0b
basic implementation of monotonic interpolation methods
mateuszbaran Sep 25, 2018
5854513
argument checks for monotonic interpolation
mateuszbaran Sep 26, 2018
28ad194
derivatives for monotonic interpolation
mateuszbaran Sep 26, 2018
1bfd8d2
Documentation for monotonic interpolation
mateuszbaran Sep 26, 2018
9d65102
faster knot searching in monotonic interpolation
mateuszbaran Sep 26, 2018
9c687b4
minor changes in monotonic interpolation
mateuszbaran Sep 26, 2018
c3763e3
more coherent type parameter naming in monotonic interpolation
mateuszbaran Sep 27, 2018
89dec3b
new gradient, gradient1, hessian and hessian1 function for monotonic …
mateuszbaran Sep 27, 2018
a19eff7
Merge branch 'master' into monotonic-interpolations
mateuszbaran Sep 27, 2018
514185a
Small fixes for Monotonic indexing rules
timholy Sep 28, 2018
a2144cf
Removal of two incorrect tests for monotonic interpolation.
mateuszbaran Sep 28, 2018
e6a38cb
Re-insertion of fixed gradient/hessian tests for monotonic interpolation
mateuszbaran Sep 28, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
63 changes: 41 additions & 22 deletions src/monotonic/monotonic.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,8 +60,8 @@ Cubic cardinal splines, uses `tension` parameter which must be between [0,1]
Cubin cardinal splines can overshoot for non-monotonic data
(increasing tension reduces overshoot).
"""
struct CardinalMonotonicInterpolation{T<:Number} <: MonotonicInterpolationType
tension :: T # must be in [0, 1]
struct CardinalMonotonicInterpolation{TTension<:Number} <: MonotonicInterpolationType
tension :: TTension # must be in [0, 1]
end

"""
Expand Down Expand Up @@ -107,12 +107,12 @@ end
Monotonic interpolation up to third order represented by type, knots and
coefficients.
"""
struct MonotonicInterpolation{T, TCoeffs, Tel, Type<:MonotonicInterpolationType,
TKnots<:AbstractVector{<:Number}, AType <: AbstractArray{Tel,1}} <: AbstractInterpolation{T,1,DimSpec{Type}}
struct MonotonicInterpolation{T, TCoeffs, TEl, TInterpolationType<:MonotonicInterpolationType,
TKnots<:AbstractVector{<:Number}, TACoeff <: AbstractArray{TEl,1}} <: AbstractInterpolation{T,1,DimSpec{TInterpolationType}}

it::Type
it::TInterpolationType
knots::TKnots
A::AType # constant parts of piecewise polynomials
A::TACoeff # constant parts of piecewise polynomials
m::Vector{TCoeffs} # coefficients of linear parts of piecewise polynomials
c::Vector{TCoeffs} # coefficients of quadratic parts of piecewise polynomials
d::Vector{TCoeffs} # coefficients of cubic parts of piecewise polynomials
Expand All @@ -122,10 +122,10 @@ end
size(A::MonotonicInterpolation) = size(A.knots)
axes(A::MonotonicInterpolation) = axes(A.knots)

function MonotonicInterpolation(::Type{TWeights}, it::IType, knots::TKnots, A::AbstractArray{Tel,1},
m::Vector{TCoeffs}, c::Vector{TCoeffs}, d::Vector{TCoeffs}) where {TWeights, TCoeffs, Tel, IType<:MonotonicInterpolationType, TKnots<:AbstractVector{<:Number}}
function MonotonicInterpolation(::Type{TWeights}, it::TInterpolationType, knots::TKnots, A::AbstractArray{TEl,1},
m::Vector{TCoeffs}, c::Vector{TCoeffs}, d::Vector{TCoeffs}) where {TWeights, TCoeffs, TEl, TInterpolationType<:MonotonicInterpolationType, TKnots<:AbstractVector{<:Number}}

isconcretetype(IType) || error("The spline type must be a leaf type (was $IType)")
isconcretetype(TInterpolationType) || error("The spline type must be a leaf type (was $TInterpolationType)")
isconcretetype(TCoeffs) || warn("For performance reasons, consider using an array of a concrete type (eltype(A) == $(eltype(A)))")

check_monotonic(knots, A)
Expand All @@ -137,11 +137,11 @@ function MonotonicInterpolation(::Type{TWeights}, it::IType, knots::TKnots, A::A
T = typeof(cZero * first(A))
end

MonotonicInterpolation{T, TCoeffs, Tel, IType, TKnots, typeof(A)}(it, knots, A, m, c, d)
MonotonicInterpolation{T, TCoeffs, TEl, TInterpolationType, TKnots, typeof(A)}(it, knots, A, m, c, d)
end

function interpolate(::Type{TWeights}, ::Type{TCoeffs}, knots::TKnots,
A::AbstractArray{Tel,1}, it::IT) where {TWeights,TCoeffs,Tel,TKnots<:AbstractVector{<:Number},IT<:MonotonicInterpolationType}
A::AbstractArray{TEl,1}, it::TInterpolationType) where {TWeights,TCoeffs,TEl,TKnots<:AbstractVector{<:Number},TInterpolationType<:MonotonicInterpolationType}

check_monotonic(knots, A)

Expand All @@ -151,7 +151,7 @@ function interpolate(::Type{TWeights}, ::Type{TCoeffs}, knots::TKnots,
c = Vector{TCoeffs}(undef, n-1)
d = Vector{TCoeffs}(undef, n-1)
for k ∈ 1:n-1
if IT == LinearMonotonicInterpolation
if TInterpolationType == LinearMonotonicInterpolation
c[k] = d[k] = zero(TCoeffs)
else
xdiff = knots[k+1] - knots[k]
Expand All @@ -163,8 +163,8 @@ function interpolate(::Type{TWeights}, ::Type{TCoeffs}, knots::TKnots,
MonotonicInterpolation(TWeights, it, knots, A, m, c, d)
end

function interpolate(knots::AbstractVector{<:Number}, A::AbstractArray{Tel,1},
it::IT) where {Tel,IT<:MonotonicInterpolationType}
function interpolate(knots::AbstractVector{<:Number}, A::AbstractArray{TEl,1},
it::TInterpolationType) where {TEl,TInterpolationType<:MonotonicInterpolationType}

interpolate(tweight(A), tcoef(A), knots, A, it)
end
Expand All @@ -180,7 +180,12 @@ function (itp::MonotonicInterpolation)(x::Number)
return itp.A[k] + itp.m[k]*xdiff + itp.c[k]*xdiff*xdiff + itp.d[k]*xdiff*xdiff*xdiff
end

function derivative(itp::MonotonicInterpolation, x::Number)
function gradient(itp::MonotonicInterpolation, x::AbstractArray{<:Number, 1})
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice! This is almost what we want, but the type of x should probably be VarArg{<:Number,1} - see e.g. how it's done for Gridded.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've just checked it and it looks like function gradient(itp::MonotonicInterpolation, x::Number) works. I think gradient isn't really a vararg function when it expects exactly two arguments :)

I tried vararg with hessian first and with function hessian(itp::MonotonicInterpolation, x::Vararg{<:Number,1}) I got stack overflow due to this function but with function hessian(itp::MonotonicInterpolation, x::Number) it just works.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

gradient is a varargs function, because the notation args::Vararg{T,N} means "supply N separate arguments of type T." But in your case since all MonotonicInterpolations are 1-dimensional, this is equivalent to gradient(itp::MonotonicInterpolation, x::Number).

x::Vararg{Number,1} should be the same as x::Number, but Vararg{<:Number,1} may not be.

n-dimensional positions are encoded as varargs; like the rest of julia, indexing with vectors means something like this:

julia> A = reshape(1:20, 4, 5)
4×5 reshape(::UnitRange{Int64}, 4, 5) with eltype Int64:
 1  5   9  13  17
 2  6  10  14  18
 3  7  11  15  19
 4  8  12  16  20

julia> A[2:3, 3:4]
2×2 Array{Int64,2}:
 10  14
 11  15

This is not what you mean here, so definitely get rid of the AbstractArray.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll change x::AbstractArray{<Number,1} to x::Number then. In this case x::Vararg{Number,1} causes a stack overflow while x::Number does not:

Test threw exception
  Expression: ≈((Interpolations.hessian(itp, [t]))[1], hessval, atol=1.0e-12)
  StackOverflowError:
  Stacktrace:
   [1] hessian(::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,FiniteDifferenceMonotonicInterpolation,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}) at /home/mateusz/.julia/dev/Interpolations/src/Interpolations.jl:359 (repeats 80000 times)

Copy link
Contributor Author

@mateuszbaran mateuszbaran Sep 28, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My mistake, I checked again and using x::Number leads to the same stack overflow error. I don't know what I should do.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, Interpolations.gradient(itp, 0.3) works, but Interpolations.gradient(itp, [0.3]) and Interpolations.hessian(itp, [0.3]) do not:

julia> using Interpolations
[ Info: Recompiling stale cache file /home/mateusz/.julia/compiled/v1.0/Interpolations/VpKVx.ji for Interpolations [a98d9a8b-a2ab-59e6-89dd-64a1c18fca59]

julia> itp = interpolate(0.0:0.1:1.0, 0.0:0.1:1.0, LinearMonotonicInterpolation());

julia> itp(0.0)
0.0

julia> itp(0.2)
0.2

julia> Interpolations.gradient(itp, 0.3)
1-element StaticArrays.SArray{Tuple{1},Float64,1,1}:
 1.0

julia> Interpolations.gradient(itp, [0.3])
ERROR: Given number 2 is outside of interpolated range [0.0, 1.0]
Stacktrace:
 [1] (::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}})(::Int64) at /home/mateusz/.julia/dev/Interpolations/src/monotonic/monotonic.jl:174
 [2] getindex at /home/mateusz/.julia/dev/Interpolations/src/Interpolations.jl:383 [inlined]
 [3] isassigned(::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}, ::Int64) at ./abstractarray.jl:351
 [4] show_delim_array(::IOContext{Base.GenericIOBuffer{Array{UInt8,1}}}, ::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}, ::Char, ::String, ::Char, ::Bool, ::Int64, ::Int64) at ./show.jl:659
 [5] show_delim_array at ./show.jl:649 [inlined]
 [6] show_vector(::Base.GenericIOBuffer{Array{UInt8,1}}, ::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}, ::Char, ::Char) at ./arrayshow.jl:442
 [7] show_vector at ./arrayshow.jl:432 [inlined]
 [8] show at ./arrayshow.jl:418 [inlined]
 [9] print(::Base.GenericIOBuffer{Array{UInt8,1}}, ::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}) at ./strings/io.jl:31
 [10] #print_to_string#326(::Nothing, ::Function, ::String, ::Vararg{Any,N} where N) at ./strings/io.jl:124
 [11] print_to_string at ./strings/io.jl:112 [inlined]
 [12] string at ./strings/io.jl:143 [inlined]
 [13] gradient(::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}, ::Array{Float64,1}) at /home/mateusz/.julia/dev/Interpolations/src/Interpolations.jl:349
 [14] top-level scope at none:0

julia> Interpolations.hessian(itp, [0.3])
ERROR: StackOverflowError:
Stacktrace:
 [1] hessian(::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}, ::Array{Float64,1}) at /home/mateusz/.julia/dev/Interpolations/src/Interpolations.jl:359 (repeats80000 times)

julia>

These error messages aren't helpful.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I just pushed some fixes to your branch.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks! Should I do something about that problem with gradient/hessian with array as an argument?

Copy link
Member

@timholy timholy Sep 28, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you want to support Interpolations.gradient(itp, [0.3]) you need a special method that returns an array of gradients, e.g., a Vector{SVector{1,Float64}}. But none of the other methods here support such a syntax:

julia> itp = interpolate(1:10, BSpline(Linear()))
10-element interpolate(::Array{Float64,1}, BSpline(Linear())) with element type Float64:
  1.0
  2.0
  3.0
  4.0
  5.0
  6.0
  7.0
  8.0
  9.0
 10.0

julia> Interpolations.gradient(itp, 2.3)
1-element StaticArrays.SArray{Tuple{1},Float64,1,1}:
 1.0

julia> Interpolations.gradient(itp, [2.3])
ERROR: gradient of [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0] not supported for position ([2.3],)
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] gradient(::Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},BSpline{Linear},Tuple{Base.OneTo{Int64}}}, ::Array{Float64,1}) at /home/tim/.julia/dev/Interpolations/src/Interpolations.jl:349
 [3] top-level scope at none:0

so I don't think you have to, either.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, so I'll just remove the failing tests for gradient/hessian with array as an argument.

length(x) == 1 || error("Given vector x ($x) should have exactly one element")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With x::VarArg{<:Number,1}, I don't think this check is needed (but I'm not sure).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This check isn't needed with function gradient(itp::MonotonicInterpolation, x::Number) and it works.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah; for monotonic interpolations it doesn't matter very much.

My main concern here is that I want the interface to match gradient for other interpolation types exactly, with the only exception being the type of itp. The idea is that you should be able to have code that initializes an interpolation object of some type, and then (perhaps a lot later in the code) use any part of the common API without caring which type you used. So if you e.g. start out with regular B-splines, and then realize you want monotonic splines, you switch out the initialization part and nothing else.

It might be that x::Number accomplishes that well enough, and if so it's totally fine by me. I just want to make sure that the signature of this method (as well as the others, e.g. itp(<whatever>) and hessian(itp, <whatever>)) is designed with the common API in mind.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I understand, but Vararg{<:Number,1} didn't work (I wrote about it here: #238 (comment) ). gradient(itp, [1.0]) still works (UnexpandedIndexTypes union contains AbstractVector). Is there another use case I should keep in mind?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hum. I'll ping in @timholy here...

return SVector(gradient1(itp, x[1]))
end

function gradient1(itp::MonotonicInterpolation, x::Number)
k = searchsortedfirst(itp.knots, x)
if k > 1
k -= 1
Expand All @@ -189,13 +194,27 @@ function derivative(itp::MonotonicInterpolation, x::Number)
return itp.m[k] + 2*itp.c[k]*xdiff + 3*itp.d[k]*xdiff*xdiff
end

function hessian(itp::MonotonicInterpolation, x::AbstractArray{<:Number, 1})
length(x) == 1 || error("Given vector x ($x) should have exactly one element")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same thing here: x should be VarArg, and the length check is then probably redundant.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No AbstractArray here.

return SVector(hessian1(itp, x[1]))
end

function hessian1(itp::MonotonicInterpolation, x::Number)
k = searchsortedfirst(itp.knots, x)
if k > 1
k -= 1
end
xdiff = x - itp.knots[k]
return 2*itp.c[k] + 6*itp.d[k]*xdiff
end

@inline function check_monotonic(knots, A)
axes(knots) == axes(A) || throw(DimensionMismatch("knot vector must have the same axes as the corresponding array"))
issorted(knots) || error("knot-vector must be sorted in increasing order")
end

function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These methods could do with some more documentation, that explains the mathematical reasoning behind the algorithm (or at least refers to resources that explain it). The Fritsch & Carlson version has good docs, but the other ones have none...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This code was taken almost straight from the implementation provided by niclasmattsson in discussion for issue #105 . I haven't really looked into the details except for the interpretation of coefficients A, m, c and d. I've added some documentation for these vectors in a new commit.

y::AbstractVector{Tel}, method::LinearMonotonicInterpolation) where {TCoeffs, Tel}
y::AbstractVector{TEl}, method::LinearMonotonicInterpolation) where {TCoeffs, TEl}

n = length(x)
Δ = Vector{TCoeffs}(undef, n-1)
Expand All @@ -210,7 +229,7 @@ function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
end

function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
y::AbstractVector{Tel}, method::FiniteDifferenceMonotonicInterpolation) where {TCoeffs, Tel}
y::AbstractVector{TEl}, method::FiniteDifferenceMonotonicInterpolation) where {TCoeffs, TEl}

n = length(x)
Δ = Vector{TCoeffs}(undef, n-1)
Expand All @@ -229,7 +248,7 @@ function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
end

function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
y::AbstractVector{Tel}, method::CardinalMonotonicInterpolation{T}) where {T, TCoeffs, Tel}
y::AbstractVector{TEl}, method::CardinalMonotonicInterpolation{TTension}) where {TTension, TCoeffs, TEl}

n = length(x)
Δ = Vector{TCoeffs}(undef, n-1)
Expand All @@ -240,15 +259,15 @@ function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
if k == 1 # left endpoint
m[k] = Δk
else
m[k] = (oneunit(T) - method.tension) * (y[k+1] - y[k-1]) / (x[k+1] - x[k-1])
m[k] = (oneunit(TTension) - method.tension) * (y[k+1] - y[k-1]) / (x[k+1] - x[k-1])
end
end
m[n] = Δ[n-1]
return (m, Δ)
end

function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
y::AbstractVector{Tel}, method::FritschCarlsonMonotonicInterpolation) where {TCoeffs, Tel}
y::AbstractVector{TEl}, method::FritschCarlsonMonotonicInterpolation) where {TCoeffs, TEl}

n = length(x)
Δ = Vector{TCoeffs}(undef, n-1)
Expand Down Expand Up @@ -298,7 +317,7 @@ function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
end

function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
y :: AbstractVector{Tel}, method :: FritschButlandMonotonicInterpolation) where {TCoeffs, Tel}
y :: AbstractVector{TEl}, method :: FritschButlandMonotonicInterpolation) where {TCoeffs, TEl}

# based on Fritsch & Butland (1984),
# "A Method for Constructing Local Monotone Piecewise Cubic Interpolants",
Expand All @@ -324,7 +343,7 @@ function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
end

function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
y::AbstractVector{Tel}, method::SteffenMonotonicInterpolation) where {TCoeffs, Tel}
y::AbstractVector{TEl}, method::SteffenMonotonicInterpolation) where {TCoeffs, TEl}

# Steffen (1990),
# "A Simple Method for Monotonic Interpolation in One Dimension",
Expand Down
4 changes: 3 additions & 1 deletion test/gradient.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -182,7 +182,9 @@ using Test, Interpolations, DualNumbers, LinearAlgebra
for it in itypes
itp = interpolate(x, y, it)
for t in grid
@test Interpolations.derivative(itp, t) ≈ epsilon(itp(dual(t, 1.0))) atol = 1.e-12
gradval = epsilon(itp(dual(t, 1.0)))
@test Interpolations.gradient1(itp, t) ≈ gradval atol = 1.e-12
@test Interpolations.gradient(itp, [t])[1] ≈ gradval atol = 1.e-12
end
end
end
Expand Down
30 changes: 29 additions & 1 deletion test/hessian.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
using Test, Interpolations, LinearAlgebra
using Test, Interpolations, LinearAlgebra, ForwardDiff

@testset "Hessians" begin
nx = 5
Expand Down Expand Up @@ -42,4 +42,32 @@ using Test, Interpolations, LinearAlgebra
v = A2[:, 2]
itpcol = interpolate(v, BSpline(Quadratic(Flat(OnCell()))))
@test Interpolations.hessian(itp, 3.2, 2) == Interpolations.hessian(itpcol, 3.2)


@testset "Monotonic" begin
x = [0.0, 0.2, 0.5, 0.6, 0.9, 1.0]
ys = [[-3.0, 0.0, 5.0, 10.0, 18.0, 22.0],
[10.0, 0.0, -5.0, 10.0, -8.0, -2.0]]
grid = 0.0:0.1:1.0

itypes = [LinearMonotonicInterpolation(),
FiniteDifferenceMonotonicInterpolation(),
CardinalMonotonicInterpolation(0.0),
CardinalMonotonicInterpolation(0.5),
CardinalMonotonicInterpolation(1.0),
FritschCarlsonMonotonicInterpolation(),
FritschButlandMonotonicInterpolation(),
SteffenMonotonicInterpolation()]

for y in ys
for it in itypes
itp = interpolate(x, y, it)
for t in grid
hessval = ForwardDiff.hessian(u -> itp(u[1]), [t])[1, 1]
@test Interpolations.hessian1(itp, t) ≈ hessval atol = 1.e-12
@test Interpolations.hessian(itp, [t])[1] ≈ hessval atol = 1.e-12
end
end
end
end
end
2 changes: 1 addition & 1 deletion test/plots/monotonic/sanity.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ x = [0.0, 0.2, 0.5, 0.6, 0.9, 1.0]
y = [10.0, 0.0, -5.0, 10.0, -8.0, -2.0]
grid = 0.0:0.002:1.0

itypes = MonotonicInterpolationType[LinearMonotonicInterpolation(),
itypes = [LinearMonotonicInterpolation(),
FiniteDifferenceMonotonicInterpolation(),
CardinalMonotonicInterpolation(0.0),
CardinalMonotonicInterpolation(0.5),
Expand Down