From c3d3d0b23f3e62d0fe91d0eb2fc48a7afd18ed4f Mon Sep 17 00:00:00 2001
From: Mateusz Baran <mateuszbaran89@gmail.com>
Date: Tue, 25 Sep 2018 20:03:00 +0200
Subject: [PATCH 01/11] basic implementation of monotonic interpolation methods

implemented algorithms:
* linear
* finite difference
* cardinal
* Fritsch-Carlson
* Fritsch-Butland
* Steffen
---
 src/Interpolations.jl          |   2 +
 src/io.jl                      |  14 ++
 src/monotonic/monotonic.jl     | 275 +++++++++++++++++++++++++++++++++
 test/io.jl                     |   9 +-
 test/monotonic/runtests.jl     |  43 ++++++
 test/plots/monotonic/sanity.jl |  25 +++
 test/runtests.jl               |   3 +
 7 files changed, 370 insertions(+), 1 deletion(-)
 create mode 100644 src/monotonic/monotonic.jl
 create mode 100644 test/monotonic/runtests.jl
 create mode 100644 test/plots/monotonic/sanity.jl

diff --git a/src/Interpolations.jl b/src/Interpolations.jl
index 3bae69ed..035362ad 100644
--- a/src/Interpolations.jl
+++ b/src/Interpolations.jl
@@ -28,6 +28,7 @@ export
     # see the following files for further exports:
     # b-splines/b-splines.jl
     # extrapolation/extrapolation.jl
+    # monotonic/monotonic.jl
     # scaling/scaling.jl
 
 using LinearAlgebra, SparseArrays
@@ -382,6 +383,7 @@ end
 include("nointerp/nointerp.jl")
 include("b-splines/b-splines.jl")
 include("gridded/gridded.jl")
+include("monotonic/monotonic.jl")
 include("extrapolation/extrapolation.jl")
 include("scaling/scaling.jl")
 include("utils.jl")
diff --git a/src/io.jl b/src/io.jl
index da0bdffb..eae584f2 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -44,6 +44,20 @@ function Base.showarg(io::IO, A::ScaledInterpolation{T}, toplevel) where {T}
     end
 end
 
+function Base.showarg(io::IO, A::MonotonicInterpolation{T, TCoeffs, Tel, Type, K}, toplevel) where {T, TCoeffs, Tel, Type, K}
+    print(io, "interpolate(")
+    _showknots(io, A.knots)
+    print(io, ", ")
+    Base.showarg(io, A.A, false)
+    print(io, ", ")
+    show(io, A.it)
+    if toplevel
+        print(io, ") with element type ",T)
+    else
+        print(io, ')')
+    end
+end
+
 function Base.showarg(io::IO, A::Extrapolation{T,N,TI,IT,ET}, toplevel) where {T,N,TI,IT,ET}
     print(io, "extrapolate(")
     Base.showarg(io, A.itp, false)
diff --git a/src/monotonic/monotonic.jl b/src/monotonic/monotonic.jl
new file mode 100644
index 00000000..70b0e0cf
--- /dev/null
+++ b/src/monotonic/monotonic.jl
@@ -0,0 +1,275 @@
+
+#=
+    Interpolate using cubic Hermite splines. The breakpoints in arrays xbp and ybp are assumed to be sorted.
+    Evaluate the function in all points of the array xeval.
+    Methods:
+        "Linear"                yuck
+        "FiniteDifference"      classic cubic interpolation, no tension parameter
+                                Finite difference can overshoot for non-monotonic data
+        "Cardinal"              cubic cardinal splines, uses tension parameter which must be between [0,1]
+                                cubin cardinal splines can overshoot for non-monotonic data
+                                (increasing tension decreases overshoot)
+        "FritschCarlson"        monotonic - tangents are first initialized, then adjusted if they are not monotonic
+                                can overshoot for non-monotonic data
+        "FritschButland"        monotonic - faster algorithm (only requires one pass) but somewhat higher apparent "tension"
+        "Steffen"               monotonic - also only one pass, results usually between FritschCarlson and FritschButland
+    Sources:
+        Fritsch & Carlson (1980), "Monotone Piecewise Cubic Interpolation", doi:10.1137/0717021.
+        Fritsch & Butland (1984), "A Method for Constructing Local Monotone Piecewise Cubic Interpolants", doi:10.1137/0905021.
+        Steffen (1990), "A Simple Method for Monotonic Interpolation in One Dimension", http://adsabs.harvard.edu/abs/1990A%26A...239..443S
+
+    Implementation based on http://bl.ocks.org/niclasmattsson/7bceb05fba6c71c78d507adae3d29417
+=#
+
+export
+    MonotonicInterpolationType,
+    LinearMonotonicInterpolation,
+    FiniteDifferenceMonotonicInterpolation,
+    CardinalMonotonicInterpolation,
+    FritschCarlsonMonotonicInterpolation,
+    FritschButlandMonotonicInterpolation,
+    SteffenMonotonicInterpolation
+
+abstract type MonotonicInterpolationType <: InterpolationType end
+
+struct LinearMonotonicInterpolation <: MonotonicInterpolationType
+end
+
+struct FiniteDifferenceMonotonicInterpolation <: MonotonicInterpolationType
+end
+
+struct CardinalMonotonicInterpolation{T<:Number} <: MonotonicInterpolationType
+    tension :: T # must be in [0, 1]
+end
+
+struct FritschCarlsonMonotonicInterpolation <: MonotonicInterpolationType
+end
+
+struct FritschButlandMonotonicInterpolation <: MonotonicInterpolationType
+end
+
+struct SteffenMonotonicInterpolation <: MonotonicInterpolationType
+end
+
+struct MonotonicInterpolation{T, TCoeffs, Tel, Type<:MonotonicInterpolationType,
+    K<:AbstractVector{<:Number}} <: AbstractInterpolation{T,1,DimSpec{Type}}
+
+    it::Type
+    knots::K
+    A::Vector{Tel}
+    m::Vector{TCoeffs}
+    c::Vector{TCoeffs}
+    d::Vector{TCoeffs}
+end
+
+
+size(A::MonotonicInterpolation) = size(A.knots)
+axes(A::MonotonicInterpolation) = axes(A.knots)
+
+function MonotonicInterpolation(::Type{TWeights}, it::IType, knots::K, A::Vector{Tel},
+    m::Vector{TCoeffs}, c::Vector{TCoeffs}, d::Vector{TCoeffs}) where {TWeights, TCoeffs, Tel, IType<:MonotonicInterpolationType, K<:AbstractVector{<:Number}}
+
+    isconcretetype(IType) || error("The b-spline type must be a leaf type (was $IType)")
+    isconcretetype(TCoeffs) || warn("For performance reasons, consider using an array of a concrete type (eltype(A) == $(eltype(A)))")
+
+    cZero = zero(TWeights)
+    if isempty(A)
+        T = Base.promote_op(*, typeof(cZero), eltype(A))
+    else
+        T = typeof(cZero * first(A))
+    end
+
+    MonotonicInterpolation{T, TCoeffs, Tel, IType, K}(it, knots, A, m, c, d)
+end
+
+function interpolate(::Type{TWeights}, ::Type{TCoeffs}, knots::K,
+    A::AbstractArray{Tel,1}, it::IT) where {TWeights,TCoeffs,Tel,K<:AbstractVector{<:Number},IT<:MonotonicInterpolationType}
+
+    # first we need to determine tangents (m)
+    n = length(knots)
+    m, Δ = calcTangents(TCoeffs, knots, A, it)
+    c = Vector{TCoeffs}(undef, n-1)
+    d = Vector{TCoeffs}(undef, n-1)
+    for k ∈ 1:n-1
+        if IT == LinearMonotonicInterpolation
+            c[k] = d[k] = zero(TCoeffs)
+        else
+            xdiff = knots[k+1] - knots[k]
+            c[k] = (3*Δ[k] - 2*m[k] - m[k+1]) / xdiff
+            d[k] = (m[k] + m[k+1] - 2*Δ[k]) / (xdiff * xdiff)
+        end
+    end
+
+    MonotonicInterpolation(TWeights, it, knots, A, m, c, d)
+end
+
+function interpolate(knots::AbstractVector{<:Number}, A::AbstractArray{Tel,1},
+    it::IT) where {Tel,N,IT<:MonotonicInterpolationType}
+
+    interpolate(tweight(A), tcoef(A), knots, A, it)
+end
+
+function (itp::MonotonicInterpolation)(x::Number)
+    x >= itp.knots[1] || return
+    x <= itp.knots[end] || return
+
+    k = 1
+    n = length(itp.knots)
+    while k < n-1 && x > itp.knots[k+1]
+        k += 1
+    end
+    xdiff = x - itp.knots[k]
+    return itp.A[k] + itp.m[k]*xdiff + itp.c[k]*xdiff*xdiff + itp.d[k]*xdiff*xdiff*xdiff
+end
+
+function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
+    y::AbstractVector{Tel}, method::LinearMonotonicInterpolation) where {TCoeffs, Tel}
+
+    n = length(x)
+    Δ = Vector{TCoeffs}(undef, n-1)
+    m = Vector{TCoeffs}(undef, n)
+    for k in 1:n-1
+        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
+        Δ[k] = Δk
+        m[k] = Δk
+    end
+    m[n] = Δ[n-1]
+    return (m, Δ)
+end
+
+function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
+    y::AbstractVector{Tel}, method::FiniteDifferenceMonotonicInterpolation) where {TCoeffs, Tel}
+
+    n = length(x)
+    Δ = Vector{TCoeffs}(undef, n-1)
+    m = Vector{TCoeffs}(undef, n)
+    for k in 1:n-1
+        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
+        Δ[k] = Δk
+        if k == 1   # left endpoint
+            m[k] = Δk
+        else
+            m[k] = (Δ[k-1] + Δk) / 2
+        end
+    end
+    m[n] = Δ[n-1]
+    return (m, Δ)
+end
+
+function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
+    y::AbstractVector{Tel}, method::CardinalMonotonicInterpolation{T}) where {T, TCoeffs, Tel}
+
+    n = length(x)
+    Δ = Vector{TCoeffs}(undef, n-1)
+    m = Vector{TCoeffs}(undef, n)
+    for k in 1:n-1
+        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
+        Δ[k] = Δk
+        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])
+        end
+    end
+    m[n] = Δ[n-1]
+    return (m, Δ)
+end
+
+function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
+    y::AbstractVector{Tel}, method::FritschCarlsonMonotonicInterpolation) where {TCoeffs, Tel}
+
+    n = length(x)
+    Δ = Vector{TCoeffs}(undef, n-1)
+    m = Vector{TCoeffs}(undef, n)
+    for k in 1:n-1
+        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
+        Δ[k] = Δk
+        if k == 1   # left endpoint
+            m[k] = Δk
+        else
+            # If any consecutive secant lines change sign (i.e. curve changes direction), initialize the tangent to zero.
+            # This is needed to make the interpolation monotonic. Otherwise set tangent to the average of the secants.
+            if Δk <= zero(Δk)
+                m[k] = zero(TCoeffs)
+            else
+                m[k] = Δ[k-1] * (Δ[k-1] + Δk) / 2.0
+            end
+        end
+    end
+    m[n] = Δ[n-1]
+    #=
+    Fritsch & Carlson derived necessary and sufficient conditions for monotonicity in their 1980 paper. Splines will be
+    monotonic if all tangents are in a certain region of the alpha-beta plane, with alpha and beta as defined below.
+    A robust choice is to put alpha & beta within a circle around origo with radius 3. The FritschCarlson algorithm
+    makes simple initial estimates of tangents and then does another pass over data points to move any outlier tangents
+    into the monotonic region. FritschButland & Steffen algorithms make more elaborate first estimates of tangents that
+    are guaranteed to lie in the monotonic region, so no second pass is necessary.
+    =#
+
+    # Second pass of FritschCarlson: adjust any non-monotonic tangents.
+    for k in 1:n-1
+        Δk = Δ[k]
+        if Δk == zero(TCoeffs)
+            m[k] = zero(TCoeffs)
+            m[k+1] = zero(TCoeffs)
+            continue
+        end
+        α = m[k] / Δk
+        β = m[k+1] / Δk
+        τ = 3.0 / sqrt(α^2 + β^2)
+        if τ < 1.0 # if we're outside the circle with radius 3 then move onto the circle
+            m[k] = τ * α * Δk
+            m[k+1] = τ * β * Δk
+        end
+    end
+    return (m, Δ)
+end
+
+function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
+    y :: AbstractVector{Tel}, method :: FritschButlandMonotonicInterpolation) where {TCoeffs, Tel}
+
+    n = length(x)
+    Δ = Vector{TCoeffs}(undef, n-1)
+    m = Vector{TCoeffs}(undef, n)
+    for k in 1:n-1
+        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
+        Δ[k] = Δk
+        if k == 1   # left endpoint
+            m[k] = Δk
+        elseif Δ[k-1] * Δk <= zero(TCoeffs)
+            m[k] = zero(TCoeffs)
+        else
+            α = (1.0 + (x[k+1] - x[k]) / (x[k+1] - x[k-1])) / 3.0
+            m[k] = Δ[k-1] * Δk / (α*Δk + (1.0 - α)*Δ[k-1])
+        end
+    end
+    m[n] = Δ[n-1]
+    return (m, Δ)
+end
+
+function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
+    y::AbstractVector{Tel}, method::SteffenMonotonicInterpolation) where {TCoeffs, Tel}
+
+    n = length(x)
+    Δ = Vector{TCoeffs}(undef, n-1)
+    m = Vector{TCoeffs}(undef, n)
+    for k in 1:n-1
+        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
+        Δ[k] = Δk
+        if k == 1   # left endpoint
+            m[k] = Δk
+        else
+            p = ((x[k+1] - x[k]) * Δ[k-1] + (x[k] - x[k-1]) * Δk) / (x[k+1] - x[k-1])
+            m[k] = (sign(Δ[k-1]) + sign(Δk)) *
+                min(abs(Δ[k-1]), abs(Δk), 0.5*abs(p))
+        end
+    end
+    m[n] = Δ[n-1]
+    return (m, Δ)
+end
+
+# How many non-NoInterp dimensions are there?
+count_interp_dims(::Type{<:MonotonicInterpolationType}) = 1
+
+lbounds(itp::MonotonicInterpolation) = first(itp.knots)
+ubounds(itp::MonotonicInterpolation) = last(itp.knots)
diff --git a/test/io.jl b/test/io.jl
index 6027c13e..5c5a9768 100644
--- a/test/io.jl
+++ b/test/io.jl
@@ -57,6 +57,13 @@ using Test
         @test summary(sitp2) == "21×41 scale(interpolate(OffsetArray(::Array{Float64,2}, 0:22, 0:42), BSpline(Quadratic(Flat(OnGrid())))), (-5.0:0.5:5.0,$SPACE-4.0:0.2:4.0)) with element type Float64"
     end
 
+    @testset "Monotonic" begin
+        A = rand(20)
+        A_x = collect(1.0:2.0:40.0)
+        itp = interpolate(A_x, A, FritschButlandMonotonicInterpolation())
+        @test summary(itp) == "20-element interpolate(::Array{Float64,1}, ::Array{Float64,1}, FritschButlandMonotonicInterpolation()) with element type Float64"
+    end
+
     @testset "Extrapolation" begin
         A = rand(8,20)
 
@@ -67,4 +74,4 @@ using Test
         etpf = extrapolate(itpg, NaN)
         @test summary(etpf) == "8×20 extrapolate(interpolate(::Array{Float64,2}, BSpline(Linear())), NaN) with element type Float64"
     end
-end # Module
+end # testset
diff --git a/test/monotonic/runtests.jl b/test/monotonic/runtests.jl
new file mode 100644
index 00000000..f56aca1a
--- /dev/null
+++ b/test/monotonic/runtests.jl
@@ -0,0 +1,43 @@
+
+using Test
+using Interpolations
+
+@testset "Monotonic" begin
+    x = [0.0, 0.2, 0.5, 0.6, 0.9, 1.0]
+    y = [(true, [-3.0, 0.0, 5.0, 10.0, 18.0, 22.0]),
+         (false, [10.0, 0.0, -5.0, 10.0, -8.0, -2.0])]
+
+    # second item indicates if interpolating function can overshoot
+    # for non-monotonic data
+    itypes = [(LinearMonotonicInterpolation(), false),
+        (FiniteDifferenceMonotonicInterpolation(), true),
+        (CardinalMonotonicInterpolation(0.0), true),
+        (CardinalMonotonicInterpolation(0.5), true),
+        (CardinalMonotonicInterpolation(1.0), false),
+        (FritschCarlsonMonotonicInterpolation(), true),
+        (FritschButlandMonotonicInterpolation(), false),
+        (SteffenMonotonicInterpolation(), false)]
+
+    for (it, overshoot) in itypes
+        for yi in 1:length(y)
+            monotonic, ys = y[yi]
+            itp = interpolate(x, ys, it)
+            for j in 1:6
+                # checking values at nodes
+                @test itp(x[j]) ≈ ys[j] rtol = 1.e-15 atol = 1.e-14
+
+                # checking overshoot for non-monotonic data
+                # and monotonicity for monotonic data
+                if !monotonic && overshoot
+                    continue
+                end
+                if j < 6
+                    r = range(x[j], stop = x[j+1], length = 100)
+                    for k in 1:length(r)-1
+                        @test (ys[j+1] - ys[j]) * (itp(r[k+1]) - itp(r[k])) >= 0.0
+                    end
+                end
+            end
+        end
+    end
+end
diff --git a/test/plots/monotonic/sanity.jl b/test/plots/monotonic/sanity.jl
new file mode 100644
index 00000000..b5505cd6
--- /dev/null
+++ b/test/plots/monotonic/sanity.jl
@@ -0,0 +1,25 @@
+using Plots
+using Interpolations
+
+pyplot()
+
+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(),
+    FiniteDifferenceMonotonicInterpolation(),
+    CardinalMonotonicInterpolation(0.0),
+    CardinalMonotonicInterpolation(0.5),
+    CardinalMonotonicInterpolation(1.0),
+    FritschCarlsonMonotonicInterpolation(),
+    FritschButlandMonotonicInterpolation(),
+    SteffenMonotonicInterpolation()
+    ]
+
+scatter(x, y, label = "data")
+for itype in itypes
+    itp = interpolate(x, y, itype)
+    plot!([x for x in grid], [itp(x) for x in grid], label = string(itype))
+end
+gui()
diff --git a/test/runtests.jl b/test/runtests.jl
index c31533d9..8d37d5ac 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -25,6 +25,9 @@ const isci = get(ENV, "CI", "") in ("true", "True")
     # scaling tests
     include("scaling/runtests.jl")
 
+    # monotonic tests
+    include("monotonic/runtests.jl")
+
     # test gradient evaluation
     include("gradient.jl")
     isci && println("finished gradient")

From 58545137f4d4d0fa201611c6d6a753af557d4248 Mon Sep 17 00:00:00 2001
From: Mateusz Baran <mateuszbaran89@gmail.com>
Date: Wed, 26 Sep 2018 12:04:14 +0200
Subject: [PATCH 02/11] argument checks for monotonic interpolation

---
 src/monotonic/monotonic.jl | 17 +++++++++++++----
 test/monotonic/runtests.jl | 11 +++++++++++
 2 files changed, 24 insertions(+), 4 deletions(-)

diff --git a/src/monotonic/monotonic.jl b/src/monotonic/monotonic.jl
index 70b0e0cf..22d06af7 100644
--- a/src/monotonic/monotonic.jl
+++ b/src/monotonic/monotonic.jl
@@ -52,11 +52,11 @@ struct SteffenMonotonicInterpolation <: MonotonicInterpolationType
 end
 
 struct MonotonicInterpolation{T, TCoeffs, Tel, Type<:MonotonicInterpolationType,
-    K<:AbstractVector{<:Number}} <: AbstractInterpolation{T,1,DimSpec{Type}}
+    K<:AbstractVector{<:Number}, AType <: AbstractArray{Tel,1}} <: AbstractInterpolation{T,1,DimSpec{Type}}
 
     it::Type
     knots::K
-    A::Vector{Tel}
+    A::AType
     m::Vector{TCoeffs}
     c::Vector{TCoeffs}
     d::Vector{TCoeffs}
@@ -66,12 +66,14 @@ end
 size(A::MonotonicInterpolation) = size(A.knots)
 axes(A::MonotonicInterpolation) = axes(A.knots)
 
-function MonotonicInterpolation(::Type{TWeights}, it::IType, knots::K, A::Vector{Tel},
+function MonotonicInterpolation(::Type{TWeights}, it::IType, knots::K, A::AbstractArray{Tel,1},
     m::Vector{TCoeffs}, c::Vector{TCoeffs}, d::Vector{TCoeffs}) where {TWeights, TCoeffs, Tel, IType<:MonotonicInterpolationType, K<:AbstractVector{<:Number}}
 
     isconcretetype(IType) || error("The b-spline type must be a leaf type (was $IType)")
     isconcretetype(TCoeffs) || warn("For performance reasons, consider using an array of a concrete type (eltype(A) == $(eltype(A)))")
 
+    check_monotonic(knots, A)
+
     cZero = zero(TWeights)
     if isempty(A)
         T = Base.promote_op(*, typeof(cZero), eltype(A))
@@ -79,12 +81,14 @@ function MonotonicInterpolation(::Type{TWeights}, it::IType, knots::K, A::Vector
         T = typeof(cZero * first(A))
     end
 
-    MonotonicInterpolation{T, TCoeffs, Tel, IType, K}(it, knots, A, m, c, d)
+    MonotonicInterpolation{T, TCoeffs, Tel, IType, K, typeof(A)}(it, knots, A, m, c, d)
 end
 
 function interpolate(::Type{TWeights}, ::Type{TCoeffs}, knots::K,
     A::AbstractArray{Tel,1}, it::IT) where {TWeights,TCoeffs,Tel,K<:AbstractVector{<:Number},IT<:MonotonicInterpolationType}
 
+    check_monotonic(knots, A)
+
     # first we need to determine tangents (m)
     n = length(knots)
     m, Δ = calcTangents(TCoeffs, knots, A, it)
@@ -122,6 +126,11 @@ 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
 
+@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},
     y::AbstractVector{Tel}, method::LinearMonotonicInterpolation) where {TCoeffs, Tel}
 
diff --git a/test/monotonic/runtests.jl b/test/monotonic/runtests.jl
index f56aca1a..9d71d221 100644
--- a/test/monotonic/runtests.jl
+++ b/test/monotonic/runtests.jl
@@ -40,4 +40,15 @@ using Interpolations
             end
         end
     end
+
+    # fail tests
+    xWrong = [0.0, 1.0, -5.0, 3.0]
+    y2 = [1.0, 2.0, 3.0, 4.0]
+    y3 = [-3.0 0.0 5.0 10.0 18.0 22.0]
+
+    for (it, overshoot) in itypes
+        @test_throws ErrorException interpolate(xWrong, y2, it)
+        @test_throws DimensionMismatch interpolate(x, y2, it)
+        @test_throws MethodError interpolate(x, y3, it)
+    end
 end

From 28ad194a0f3709fd66cc012a6da557deab0c4bf1 Mon Sep 17 00:00:00 2001
From: Mateusz Baran <mateuszbaran89@gmail.com>
Date: Wed, 26 Sep 2018 12:36:50 +0200
Subject: [PATCH 03/11] derivatives for monotonic interpolation

---
 src/monotonic/monotonic.jl | 17 ++++++++++++++---
 test/gradient.jl           | 24 ++++++++++++++++++++++++
 2 files changed, 38 insertions(+), 3 deletions(-)

diff --git a/src/monotonic/monotonic.jl b/src/monotonic/monotonic.jl
index 22d06af7..4a084cfe 100644
--- a/src/monotonic/monotonic.jl
+++ b/src/monotonic/monotonic.jl
@@ -113,19 +113,30 @@ function interpolate(knots::AbstractVector{<:Number}, A::AbstractArray{Tel,1},
     interpolate(tweight(A), tcoef(A), knots, A, it)
 end
 
-function (itp::MonotonicInterpolation)(x::Number)
-    x >= itp.knots[1] || return
-    x <= itp.knots[end] || return
+@inline function findKnot(itp::MonotonicInterpolation, x::Number)
+    x >= itp.knots[1] || error("Given number $x is outside of interpolated range [$(itp.knots[1]), $(itp.knots[end])]")
+    x <= itp.knots[end] || error("Given number $x is outside of interpolated range [$(itp.knots[1]), $(itp.knots[end])]")
 
     k = 1
     n = length(itp.knots)
     while k < n-1 && x > itp.knots[k+1]
         k += 1
     end
+    return k
+end
+
+function (itp::MonotonicInterpolation)(x::Number)
+    k = findKnot(itp, x)
     xdiff = x - itp.knots[k]
     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)
+    k = findKnot(itp, x)
+    xdiff = x - itp.knots[k]
+    return itp.m[k] + 2*itp.c[k]*xdiff + 3*itp.d[k]*xdiff*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")
diff --git a/test/gradient.jl b/test/gradient.jl
index e91ec703..a35661dc 100644
--- a/test/gradient.jl
+++ b/test/gradient.jl
@@ -163,4 +163,28 @@ using Test, Interpolations, DualNumbers, LinearAlgebra
         end
     end
 
+    @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
+                    @test Interpolations.derivative(itp, t) ≈ epsilon(itp(dual(t, 1.0))) atol = 1.e-12
+                end
+            end
+        end
+    end
 end

From 1bfd8d2011471e5ee95d8f17b83b98f87dd92dbf Mon Sep 17 00:00:00 2001
From: Mateusz Baran <mateuszbaran89@gmail.com>
Date: Wed, 26 Sep 2018 12:53:26 +0200
Subject: [PATCH 04/11] Documentation for monotonic interpolation

---
 src/monotonic/monotonic.jl | 58 +++++++++++++++++++++++++++++++++++++-
 1 file changed, 57 insertions(+), 1 deletion(-)

diff --git a/src/monotonic/monotonic.jl b/src/monotonic/monotonic.jl
index 4a084cfe..ee2670ca 100644
--- a/src/monotonic/monotonic.jl
+++ b/src/monotonic/monotonic.jl
@@ -22,7 +22,6 @@
 =#
 
 export
-    MonotonicInterpolationType,
     LinearMonotonicInterpolation,
     FiniteDifferenceMonotonicInterpolation,
     CardinalMonotonicInterpolation,
@@ -30,27 +29,84 @@ export
     FritschButlandMonotonicInterpolation,
     SteffenMonotonicInterpolation
 
+"""
+    MonotonicInterpolationType
+
+Abstract class for all types of monotonic interpolation.
+"""
 abstract type MonotonicInterpolationType <: InterpolationType end
 
+"""
+    LinearMonotonicInterpolation
+
+Simple linear interpolation.
+"""
 struct LinearMonotonicInterpolation <: MonotonicInterpolationType
 end
 
+"""
+    FiniteDifferenceMonotonicInterpolation
+
+Classic cubic interpolation, no tension parameter.
+Finite difference can overshoot for non-monotonic data.
+"""
 struct FiniteDifferenceMonotonicInterpolation <: MonotonicInterpolationType
 end
 
+"""
+    CardinalMonotonicInterpolation(tension)
+
+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]
 end
 
+"""
+    FritschCarlsonMonotonicInterpolation
+
+Monotonic interpolation based on Fritsch & Carlson (1980),
+"Monotone Piecewise Cubic Interpolation", doi:10.1137/0717021.
+
+Tangents are first initialized, then adjusted if they are not monotonic
+can overshoot for non-monotonic data
+"""
 struct FritschCarlsonMonotonicInterpolation <: MonotonicInterpolationType
 end
 
+"""
+    FritschButlandMonotonicInterpolation
+
+Monotonic interpolation based on  Fritsch & Butland (1984),
+"A Method for Constructing Local Monotone Piecewise Cubic Interpolants",
+doi:10.1137/0905021.
+
+Faster than FritschCarlsonMonotonicInterpolation (only requires one pass)
+but somewhat higher apparent "tension".
+"""
 struct FritschButlandMonotonicInterpolation <: MonotonicInterpolationType
 end
 
+"""
+    SteffenMonotonicInterpolation
+
+Monotonic interpolation based on Steffen (1990),
+"A Simple Method for Monotonic Interpolation in One Dimension",
+http://adsabs.harvard.edu/abs/1990A%26A...239..443S
+
+Only one pass, results usually between FritschCarlson and FritschButland.
+"""
 struct SteffenMonotonicInterpolation <: MonotonicInterpolationType
 end
 
+"""
+    MonotonicInterpolation
+
+Monotonic interpolation up to third order represented by type, knots and
+coefficients. Type is any concrete subtype of `MonotonicInterpolationType`.
+"""
 struct MonotonicInterpolation{T, TCoeffs, Tel, Type<:MonotonicInterpolationType,
     K<:AbstractVector{<:Number}, AType <: AbstractArray{Tel,1}} <: AbstractInterpolation{T,1,DimSpec{Type}}
 

From 9d651022cea95667a754ead61aad666a68e63c5b Mon Sep 17 00:00:00 2001
From: Mateusz Baran <mateuszbaran89@gmail.com>
Date: Wed, 26 Sep 2018 13:53:16 +0200
Subject: [PATCH 05/11] faster knot searching in monotonic interpolation

---
 src/monotonic/monotonic.jl | 20 ++++++++------------
 1 file changed, 8 insertions(+), 12 deletions(-)

diff --git a/src/monotonic/monotonic.jl b/src/monotonic/monotonic.jl
index ee2670ca..014be197 100644
--- a/src/monotonic/monotonic.jl
+++ b/src/monotonic/monotonic.jl
@@ -169,26 +169,22 @@ function interpolate(knots::AbstractVector{<:Number}, A::AbstractArray{Tel,1},
     interpolate(tweight(A), tcoef(A), knots, A, it)
 end
 
-@inline function findKnot(itp::MonotonicInterpolation, x::Number)
+function (itp::MonotonicInterpolation)(x::Number)
     x >= itp.knots[1] || error("Given number $x is outside of interpolated range [$(itp.knots[1]), $(itp.knots[end])]")
     x <= itp.knots[end] || error("Given number $x is outside of interpolated range [$(itp.knots[1]), $(itp.knots[end])]")
-
-    k = 1
-    n = length(itp.knots)
-    while k < n-1 && x > itp.knots[k+1]
-        k += 1
+    k = searchsortedfirst(itp.knots, x)
+    if k > 1
+        k -= 1
     end
-    return k
-end
-
-function (itp::MonotonicInterpolation)(x::Number)
-    k = findKnot(itp, x)
     xdiff = x - itp.knots[k]
     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)
-    k = findKnot(itp, x)
+    k = searchsortedfirst(itp.knots, x)
+    if k > 1
+        k -= 1
+    end
     xdiff = x - itp.knots[k]
     return itp.m[k] + 2*itp.c[k]*xdiff + 3*itp.d[k]*xdiff*xdiff
 end

From 9c687b48762aa8b3b51ceb277b67c2284194ced2 Mon Sep 17 00:00:00 2001
From: Mateusz Baran <mateuszbaran89@gmail.com>
Date: Wed, 26 Sep 2018 20:00:01 +0200
Subject: [PATCH 06/11] minor changes in monotonic interpolation

* additional documentation for interpolation coefficients
* references added to `calcTangents` for Fritsch & Butland and Steffen methods
* changed the type parameter name `K` to `TKnots`
* changed "b-spline" to "spline" in error string
---
 src/monotonic/monotonic.jl | 36 ++++++++++++++++++++++--------------
 1 file changed, 22 insertions(+), 14 deletions(-)

diff --git a/src/monotonic/monotonic.jl b/src/monotonic/monotonic.jl
index 014be197..55cc66d4 100644
--- a/src/monotonic/monotonic.jl
+++ b/src/monotonic/monotonic.jl
@@ -105,27 +105,27 @@ end
     MonotonicInterpolation
 
 Monotonic interpolation up to third order represented by type, knots and
-coefficients. Type is any concrete subtype of `MonotonicInterpolationType`.
+coefficients.
 """
 struct MonotonicInterpolation{T, TCoeffs, Tel, Type<:MonotonicInterpolationType,
-    K<:AbstractVector{<:Number}, AType <: AbstractArray{Tel,1}} <: AbstractInterpolation{T,1,DimSpec{Type}}
+    TKnots<:AbstractVector{<:Number}, AType <: AbstractArray{Tel,1}} <: AbstractInterpolation{T,1,DimSpec{Type}}
 
     it::Type
-    knots::K
-    A::AType
-    m::Vector{TCoeffs}
-    c::Vector{TCoeffs}
-    d::Vector{TCoeffs}
+    knots::TKnots
+    A::AType # 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
 end
 
 
 size(A::MonotonicInterpolation) = size(A.knots)
 axes(A::MonotonicInterpolation) = axes(A.knots)
 
-function MonotonicInterpolation(::Type{TWeights}, it::IType, knots::K, A::AbstractArray{Tel,1},
-    m::Vector{TCoeffs}, c::Vector{TCoeffs}, d::Vector{TCoeffs}) where {TWeights, TCoeffs, Tel, IType<:MonotonicInterpolationType, K<:AbstractVector{<:Number}}
+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}}
 
-    isconcretetype(IType) || error("The b-spline type must be a leaf type (was $IType)")
+    isconcretetype(IType) || error("The spline type must be a leaf type (was $IType)")
     isconcretetype(TCoeffs) || warn("For performance reasons, consider using an array of a concrete type (eltype(A) == $(eltype(A)))")
 
     check_monotonic(knots, A)
@@ -137,11 +137,11 @@ function MonotonicInterpolation(::Type{TWeights}, it::IType, knots::K, A::Abstra
         T = typeof(cZero * first(A))
     end
 
-    MonotonicInterpolation{T, TCoeffs, Tel, IType, K, typeof(A)}(it, knots, A, m, c, d)
+    MonotonicInterpolation{T, TCoeffs, Tel, IType, TKnots, typeof(A)}(it, knots, A, m, c, d)
 end
 
-function interpolate(::Type{TWeights}, ::Type{TCoeffs}, knots::K,
-    A::AbstractArray{Tel,1}, it::IT) where {TWeights,TCoeffs,Tel,K<:AbstractVector{<:Number},IT<:MonotonicInterpolationType}
+function interpolate(::Type{TWeights}, ::Type{TCoeffs}, knots::TKnots,
+    A::AbstractArray{Tel,1}, it::IT) where {TWeights,TCoeffs,Tel,TKnots<:AbstractVector{<:Number},IT<:MonotonicInterpolationType}
 
     check_monotonic(knots, A)
 
@@ -164,7 +164,7 @@ function interpolate(::Type{TWeights}, ::Type{TCoeffs}, knots::K,
 end
 
 function interpolate(knots::AbstractVector{<:Number}, A::AbstractArray{Tel,1},
-    it::IT) where {Tel,N,IT<:MonotonicInterpolationType}
+    it::IT) where {Tel,IT<:MonotonicInterpolationType}
 
     interpolate(tweight(A), tcoef(A), knots, A, it)
 end
@@ -300,6 +300,10 @@ end
 function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
     y :: AbstractVector{Tel}, method :: FritschButlandMonotonicInterpolation) where {TCoeffs, Tel}
 
+    # based on Fritsch & Butland (1984),
+    # "A Method for Constructing Local Monotone Piecewise Cubic Interpolants",
+    # doi:10.1137/0905021.
+
     n = length(x)
     Δ = Vector{TCoeffs}(undef, n-1)
     m = Vector{TCoeffs}(undef, n)
@@ -322,6 +326,10 @@ end
 function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
     y::AbstractVector{Tel}, method::SteffenMonotonicInterpolation) where {TCoeffs, Tel}
 
+    # Steffen (1990),
+    # "A Simple Method for Monotonic Interpolation in One Dimension",
+    # http://adsabs.harvard.edu/abs/1990A%26A...239..443S
+
     n = length(x)
     Δ = Vector{TCoeffs}(undef, n-1)
     m = Vector{TCoeffs}(undef, n)

From c3763e3e4b150be1e87e3c4a948a0a9a787001b5 Mon Sep 17 00:00:00 2001
From: Mateusz Baran <mateuszbaran89@gmail.com>
Date: Thu, 27 Sep 2018 10:15:55 +0200
Subject: [PATCH 07/11] more coherent type parameter naming in monotonic
 interpolation

---
 src/monotonic/monotonic.jl | 42 +++++++++++++++++++-------------------
 1 file changed, 21 insertions(+), 21 deletions(-)

diff --git a/src/monotonic/monotonic.jl b/src/monotonic/monotonic.jl
index 55cc66d4..6edaac51 100644
--- a/src/monotonic/monotonic.jl
+++ b/src/monotonic/monotonic.jl
@@ -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
 
 """
@@ -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
@@ -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)
@@ -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)
 
@@ -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]
@@ -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
@@ -195,7 +195,7 @@ end
 end
 
 function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
-    y::AbstractVector{Tel}, method::LinearMonotonicInterpolation) where {TCoeffs, Tel}
+    y::AbstractVector{TEl}, method::LinearMonotonicInterpolation) where {TCoeffs, TEl}
 
     n = length(x)
     Δ = Vector{TCoeffs}(undef, n-1)
@@ -210,7 +210,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)
@@ -229,7 +229,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)
@@ -240,7 +240,7 @@ 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]
@@ -248,7 +248,7 @@ function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
 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)
@@ -298,7 +298,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",
@@ -324,7 +324,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",

From 89dec3babf9211f7b4c75c1c93663870db1a7259 Mon Sep 17 00:00:00 2001
From: Mateusz Baran <mateuszbaran89@gmail.com>
Date: Thu, 27 Sep 2018 19:11:02 +0200
Subject: [PATCH 08/11] new gradient, gradient1, hessian and hessian1 function
 for monotonic interpolation

---
 src/monotonic/monotonic.jl     | 21 ++++++++++++++++++++-
 test/gradient.jl               |  4 +++-
 test/hessian.jl                | 30 +++++++++++++++++++++++++++++-
 test/plots/monotonic/sanity.jl |  2 +-
 4 files changed, 53 insertions(+), 4 deletions(-)

diff --git a/src/monotonic/monotonic.jl b/src/monotonic/monotonic.jl
index 6edaac51..28d13770 100644
--- a/src/monotonic/monotonic.jl
+++ b/src/monotonic/monotonic.jl
@@ -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})
+    length(x) == 1 || error("Given vector x ($x) should have exactly one element")
+    return SVector(gradient1(itp, x[1]))
+end
+
+function gradient1(itp::MonotonicInterpolation, x::Number)
     k = searchsortedfirst(itp.knots, x)
     if k > 1
         k -= 1
@@ -189,6 +194,20 @@ 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")
+    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")
diff --git a/test/gradient.jl b/test/gradient.jl
index a35661dc..20e60cf5 100644
--- a/test/gradient.jl
+++ b/test/gradient.jl
@@ -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
diff --git a/test/hessian.jl b/test/hessian.jl
index 2af7574d..c128668a 100644
--- a/test/hessian.jl
+++ b/test/hessian.jl
@@ -1,4 +1,4 @@
-using Test, Interpolations, LinearAlgebra
+using Test, Interpolations, LinearAlgebra, ForwardDiff
 
 @testset "Hessians" begin
     nx = 5
@@ -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
diff --git a/test/plots/monotonic/sanity.jl b/test/plots/monotonic/sanity.jl
index b5505cd6..419e84c1 100644
--- a/test/plots/monotonic/sanity.jl
+++ b/test/plots/monotonic/sanity.jl
@@ -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),

From 514185a3e663da272e0de138587346f9d8f343ca Mon Sep 17 00:00:00 2001
From: Tim Holy <tim.holy@gmail.com>
Date: Fri, 28 Sep 2018 07:14:36 -0500
Subject: [PATCH 09/11] Small fixes for Monotonic indexing rules

---
 src/monotonic/monotonic.jl | 19 +++++++++----------
 1 file changed, 9 insertions(+), 10 deletions(-)

diff --git a/src/monotonic/monotonic.jl b/src/monotonic/monotonic.jl
index 28d13770..fb385851 100644
--- a/src/monotonic/monotonic.jl
+++ b/src/monotonic/monotonic.jl
@@ -170,8 +170,7 @@ function interpolate(knots::AbstractVector{<:Number}, A::AbstractArray{TEl,1},
 end
 
 function (itp::MonotonicInterpolation)(x::Number)
-    x >= itp.knots[1] || error("Given number $x is outside of interpolated range [$(itp.knots[1]), $(itp.knots[end])]")
-    x <= itp.knots[end] || error("Given number $x is outside of interpolated range [$(itp.knots[1]), $(itp.knots[end])]")
+    @boundscheck (checkbounds(Bool, itp, x) || Base.throw_boundserror(itp, (x,)))
     k = searchsortedfirst(itp.knots, x)
     if k > 1
         k -= 1
@@ -180,12 +179,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 gradient(itp::MonotonicInterpolation, x::AbstractArray{<:Number, 1})
-    length(x) == 1 || error("Given vector x ($x) should have exactly one element")
-    return SVector(gradient1(itp, x[1]))
+function gradient(itp::MonotonicInterpolation, x::Number)
+    return SVector(gradient1(itp, x))
 end
 
 function gradient1(itp::MonotonicInterpolation, x::Number)
+    @boundscheck (checkbounds(Bool, itp, x) || Base.throw_boundserror(itp, (x,)))
     k = searchsortedfirst(itp.knots, x)
     if k > 1
         k -= 1
@@ -194,12 +193,12 @@ function gradient1(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")
-    return SVector(hessian1(itp, x[1]))
+function hessian(itp::MonotonicInterpolation, x::Number)
+    return SVector(hessian1(itp, x))
 end
 
 function hessian1(itp::MonotonicInterpolation, x::Number)
+    @boundscheck (checkbounds(Bool, itp, x) || Base.throw_boundserror(itp, (x,)))
     k = searchsortedfirst(itp.knots, x)
     if k > 1
         k -= 1
@@ -370,5 +369,5 @@ end
 # How many non-NoInterp dimensions are there?
 count_interp_dims(::Type{<:MonotonicInterpolationType}) = 1
 
-lbounds(itp::MonotonicInterpolation) = first(itp.knots)
-ubounds(itp::MonotonicInterpolation) = last(itp.knots)
+lbounds(itp::MonotonicInterpolation) = (first(itp.knots),)
+ubounds(itp::MonotonicInterpolation) = (last(itp.knots),)

From a2144cfa7c42b994ea02998f60a04cffae45ebc5 Mon Sep 17 00:00:00 2001
From: Mateusz Baran <mateuszbaran89@gmail.com>
Date: Fri, 28 Sep 2018 14:46:20 +0200
Subject: [PATCH 10/11] Removal of two incorrect tests for monotonic
 interpolation.

They test gradient/hessian with array as an argument.
---
 test/gradient.jl | 1 -
 test/hessian.jl  | 1 -
 2 files changed, 2 deletions(-)

diff --git a/test/gradient.jl b/test/gradient.jl
index 3d65b00f..2123113f 100644
--- a/test/gradient.jl
+++ b/test/gradient.jl
@@ -194,7 +194,6 @@ using Test, Interpolations, DualNumbers, LinearAlgebra
                 for t in grid
                     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
diff --git a/test/hessian.jl b/test/hessian.jl
index c128668a..8ee7daee 100644
--- a/test/hessian.jl
+++ b/test/hessian.jl
@@ -65,7 +65,6 @@ using Test, Interpolations, LinearAlgebra, ForwardDiff
                 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

From e6a38cb4f89de2546f11b04bf5a53896bcff5ad6 Mon Sep 17 00:00:00 2001
From: Mateusz Baran <mateuszbaran89@gmail.com>
Date: Fri, 28 Sep 2018 15:59:27 +0200
Subject: [PATCH 11/11] Re-insertion of fixed gradient/hessian tests for
 monotonic interpolation

---
 test/gradient.jl | 1 +
 test/hessian.jl  | 1 +
 2 files changed, 2 insertions(+)

diff --git a/test/gradient.jl b/test/gradient.jl
index 2123113f..6c0ec136 100644
--- a/test/gradient.jl
+++ b/test/gradient.jl
@@ -194,6 +194,7 @@ using Test, Interpolations, DualNumbers, LinearAlgebra
                 for t in grid
                     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
diff --git a/test/hessian.jl b/test/hessian.jl
index 8ee7daee..297f6965 100644
--- a/test/hessian.jl
+++ b/test/hessian.jl
@@ -65,6 +65,7 @@ using Test, Interpolations, LinearAlgebra, ForwardDiff
                 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