Merge pull request #295 from SouthEndMusic/pchip

add `PCHIPInterpolation`

README.md

@@ -60,6 +60,7 @@ corresponding to `(u,t)` pairs.
       + `knotVec` - Symbol to Knot Vector, `knotVec = :Uniform` for uniform knot vector, `knotVec = :Average` for average spaced knot vector.
   - `BSplineApprox(u,t,d,h,pVec,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h<length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
   - `CubicHermiteSpline(du, u, t)` - A third order Hermite interpolation, which matches the values and first (`du`) order derivatives in the data points exactly.
+  - `PCHIPInterpolation(u, t)` - a type of `CubicHermiteSpline` where the derivative values `du` are derived from the input data in such a way that the interpolation never overshoots the data.
   - `QuinticHermiteSpline(ddu, du, u, t)` - A fifth order Hermite interpolation, which matches the values and first (`du`) and second (`ddu`) order derivatives in the data points exactly.
 
 ## Extension Methods

docs/src/index.md

@@ -35,6 +35,7 @@ corresponding to `(u,t)` pairs.
       + `knotVec` - Symbol to Knot Vector, `knotVec = :Uniform` for uniform knot vector, `knotVec = :Average` for average spaced knot vector.
   - `BSplineApprox(u,t,d,h,pVec,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h<length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
   - `CubicHermiteSpline(du, u, t)` - A third order Hermite interpolation, which matches the values and first (`du`) order derivatives in the data points exactly.
+  - `PCHIPInterpolation(u, t)` - a type of `CubicHermiteSpline` where the derivative values `du` are derived from the input data in such a way that the interpolation never overshoots the data.
   - `QuinticHermiteSpline(ddu, du, u, t)` - a fifth order Hermite interpolation, which matches the values and first (`du`) and second (`ddu`) order derivatives in the data points exactly.
 
 ## Extension Methods

docs/src/manual.md

@@ -11,5 +11,6 @@ CubicSpline
 BSplineInterpolation
 BSplineApprox
 CubicHermiteSpline
+PCHIPInterpolation
 QuinticHermiteSpline
 ```

docs/src/methods.md

@@ -145,6 +145,16 @@ scatter(t, u, label = "input data")
 plot!(A)
 ```
 
+## PCHIP Interpolation
+
+This is a type of `CubicHermiteSpline` where the derivative values `du` are derived from the input data in such a way that the interpolation never overshoots the data.
+
+```@example tutorial
+A = PCHIPInterpolation(u, t)
+scatter(t, u, label = "input data")
+plot!(A)
+```
+
 ## Quintic Hermite Spline
 
 This is the quintic (fifth order) Hermite interpolation. It matches the values and first and second order derivatives in the data points exactly.

src/DataInterpolations.jl

@@ -96,7 +96,7 @@ end
 
 export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation,
        AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline,
-       BSplineInterpolation, BSplineApprox, CubicHermiteSpline,
+       BSplineInterpolation, BSplineApprox, CubicHermiteSpline, PCHIPInterpolation,
        QuinticHermiteSpline, LinearInterpolationIntInv, ConstantInterpolationIntInv
 
 # added for RegularizationSmooth, JJS 11/27/21

src/interpolation_caches.jl

@@ -1,5 +1,5 @@
 """
-    LinearInterpolation(u, t; extrapolate = false)
+    LinearInterpolation(u, t; extrapolate = false, safetycopy = true)
 
 It is the method of interpolating between the data points using a linear polynomial. For any point, two data points one each side are chosen and connected with a line.
 Extrapolation extends the last linear polynomial on each side.
@@ -37,7 +37,7 @@ function LinearInterpolation(u, t; extrapolate = false, safetycopy = true)
 end
 
 """
-    QuadraticInterpolation(u, t, mode = :Forward; extrapolate = false)
+    QuadraticInterpolation(u, t, mode = :Forward; extrapolate = false, safetycopy = true)
 
 It is the method of interpolating between the data points using quadratic polynomials. For any point, three data points nearby are taken to fit a quadratic polynomial.
 Extrapolation extends the last quadratic polynomial on each side.
@@ -83,7 +83,7 @@ function QuadraticInterpolation(u, t; extrapolate = false, safetycopy = true)
 end
 
 """
-    LagrangeInterpolation(u, t, n = length(t) - 1; extrapolate = false)
+    LagrangeInterpolation(u, t, n = length(t) - 1; extrapolate = false, safetycopy = true)
 
 It is the method of interpolation using Lagrange polynomials of (k-1)th order passing through all the data points where k is the number of data points.
 
@@ -134,9 +134,9 @@ function LagrangeInterpolation(
 end
 
 """
-    AkimaInterpolation(u, t; extrapolate = false)
+    AkimaInterpolation(u, t; extrapolate = false, safetycopy = true)
 
-It is a spline interpolation built from cubic polynomials. It forms a continuously differentiable function. For more details, refer: https://en.wikipedia.org/wiki/Akima_spline.
+It is a spline interpolation built from cubic polynomials. It forms a continuously differentiable function. For more details, refer: [https://en.wikipedia.org/wiki/Akima_spline](https://en.wikipedia.org/wiki/Akima_spline).
 Extrapolation extends the last cubic polynomial on each side.
 
 ## Arguments
@@ -203,7 +203,7 @@ function AkimaInterpolation(u, t; extrapolate = false, safetycopy = true)
 end
 
 """
-    ConstantInterpolation(u, t; dir = :left, extrapolate = false)
+    ConstantInterpolation(u, t; dir = :left, extrapolate = false, safetycopy = true)
 
 It is the method of interpolating using a constant polynomial. For any point, two adjacent data points are found on either side (left and right). The value at that point depends on `dir`.
 If it is `:left`, then the value at the left point is chosen and if it is `:right`, the value at the right point is chosen.
@@ -243,7 +243,7 @@ function ConstantInterpolation(u, t; dir = :left, extrapolate = false, safetycop
 end
 
 """
-    QuadraticSpline(u, t; extrapolate = false)
+    QuadraticSpline(u, t; extrapolate = false, safetycopy = true)
 
 It is a spline interpolation using piecewise quadratic polynomials between each pair of data points. Its first derivative is also continuous.
 Extrapolation extends the last quadratic polynomial on each side.
@@ -329,7 +329,7 @@ function QuadraticSpline(
 end
 
 """
-    CubicSpline(u, t; extrapolate = false)
+    CubicSpline(u, t; extrapolate = false, safetycopy = true)
 
 It is a spline interpolation using piecewise cubic polynomials between each pair of data points. Its first and second derivative is also continuous.
 Second derivative on both ends are zero, which are also called "natural" boundary conditions. Extrapolation extends the last cubic polynomial on each side.
@@ -418,9 +418,9 @@ function CubicSpline(
 end
 
 """
-    BSplineInterpolation(u, t, d, pVecType, knotVecType; extrapolate = false)
+    BSplineInterpolation(u, t, d, pVecType, knotVecType; extrapolate = false, safetycopy = true)
 
-It is a curve defined by the linear combination of `n` basis functions of degree `d` where `n` is the number of data points. For more information, refer https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve.html.
+It is a curve defined by the linear combination of `n` basis functions of degree `d` where `n` is the number of data points. For more information, refer [https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve.html](https://pages.mtu.edu/%7Eshene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve.html).
 Extrapolation is a constant polynomial of the end points on each side.
 
 ## Arguments
@@ -547,10 +547,10 @@ function BSplineInterpolation(
 end
 
 """
-    BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolate = false)
+    BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolate = false, safetycopy = true)
 
 It is a regression based B-spline. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h < length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
-For more information, refer http://www.cad.zju.edu.cn/home/zhx/GM/009/00-bsia.pdf.
+For more information, refer [http://www.cad.zju.edu.cn/home/zhx/GM/009/00-bsia.pdf](http://www.cad.zju.edu.cn/home/zhx/GM/009/00-bsia.pdf).
 Extrapolation is a constant polynomial of the end points on each side.
 
 ## Arguments
@@ -702,7 +702,7 @@ function BSplineApprox(
 end
 
 """
-    CubicHermiteSpline(du, u, t; extrapolate = false)
+    CubicHermiteSpline(du, u, t; extrapolate = false, safetycopy = true)
 
 It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomial such that the value and the first derivative are equal to given values in the data points.
 
@@ -742,7 +742,31 @@ function CubicHermiteSpline(du, u, t; extrapolate = false, safetycopy = true)
 end
 
 """
-    QuinticHermiteSpline(ddu, du, u, t; extrapolate = false)
+    PCHIPInterpolation(u, t; extrapolate = false, safetycopy = true)
+
+It is a PCHIP Interpolation, which is a type of `CubicHermiteSpline` where
+the derivative values `du` are derived from the input data in such a way that
+the interpolation never overshoots the data. See also [here](https://www.mathworks.com/content/dam/mathworks/mathworks-dot-com/moler/interp.pdf),
+section 3.4.
+
+## Arguments
+
+  - `u`: data points.
+  - `t`: time points.
+
+## Keyword Arguments
+
+  - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
+"""
+function PCHIPInterpolation(u, t; extrapolate = false, safetycopy = true)
+    u, t = munge_data(u, t, safetycopy)
+    du = du_PCHIP(u, t)
+    CubicHermiteSpline(du, u, t; extrapolate, safetycopy)
+end
+
+"""
+    QuinticHermiteSpline(ddu, du, u, t; extrapolate = false, safetycopy = true)
 
 It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polynomial such that the value and the first and second derivative are equal to given values in the data points.
 

src/interpolation_utils.jl

@@ -130,3 +130,38 @@ function cumulative_integral(A)
     pushfirst!(integral_values, zero(first(integral_values)))
     return cumsum(integral_values)
 end
+
+function du_PCHIP(u, t)
+    h = diff(u)
+    δ = h ./ diff(t)
+    s = sign.(δ)
+
+    function _du(k)
+        sₖ₋₁, sₖ = if k == 1
+            s[1], s[2]
+        elseif k == lastindex(t)
+            s[end - 1], s[end]
+        else
+            s[k - 1], s[k]
+        end
+
+        if sₖ₋₁ == 0 && sₖ == 0
+            zero(eltype(δ))
+        elseif sₖ₋₁ == sₖ
+            if k == 1
+                ((2 * h[1] + h[2]) * δ[1] - h[1] * δ[2]) / (h[1] + h[2])
+            elseif k == lastindex(t)
+                ((2 * h[end] + h[end - 1]) * δ[end] - h[end] * δ[end - 1]) /
+                (h[end] + h[end - 1])
+            else
+                w₁ = 2h[k] + h[k - 1]
+                w₂ = h[k] + 2h[k - 1]
+                δ[k - 1] * δ[k] * (w₁ + w₂) / (w₁ * δ[k] + w₂ * δ[k - 1])
+            end
+        else
+            zero(eltype(δ))
+        end
+    end
+
+    return _du.(eachindex(t))
+end

test/interpolation_tests.jl

@@ -670,6 +670,18 @@ end
     @test_throws AssertionError CubicHermiteSpline(du, u, t)
 end
 
+@testset "PCHIPInterpolation" begin
+    u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
+    t = [0.0, 62.25, 109.66, 162.66, 205.8, 250.0]
+    A = PCHIPInterpolation(u, t)
+    @test A isa CubicHermiteSpline
+    ts = 0.0:0.1:250.0
+    us = A(ts)
+    @test all(minimum(u) .<= us)
+    @test all(maximum(u) .>= us)
+    @test all(A.du[3:4] .== 0.0)
+end
+
 @testset "Quintic Hermite Spline" begin
     test_interpolation_type(QuinticHermiteSpline)