Merge pull request #280 from SouthEndMusic/quintic_hermite_interpolation

Cubic & Quintic Hermite Spline

README.md

@@ -59,6 +59,8 @@ corresponding to `(u,t)` pairs.
       + `pVec` - Symbol to Parameters Vector, `pVec = :Uniform` for uniform spaced parameters and `pVec = :ArcLen` for parameters generated by chord length method.
       + `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.
+  - `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

@@ -34,6 +34,8 @@ corresponding to `(u,t)` pairs.
       + `pVec` - Symbol to Parameters Vector, `pVec = :Uniform` for uniformly spaced parameters, and `pVec = :ArcLen` for parameters generated by the chord length method.
       + `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.
+  - `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

@@ -10,4 +10,6 @@ QuadraticSpline
 CubicSpline
 BSplineInterpolation
 BSplineApprox
+CubicHermiteSpline
+QuinticHermiteSpline
 ```

docs/src/methods.md

@@ -134,6 +134,29 @@ scatter(t, u, label = "input data")
 plot!(A)
 ```
 
+## Cubic Hermite Spline
+
+This is the cubic (third order) Hermite interpolation. It matches the values and first order derivatives in the data points exactly.
+
+```@example tutorial
+du = [-0.047, -0.058, 0.054, 0.012, -0.068, 0.0011]
+A = CubicHermiteSpline(du, 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.
+
+```@example tutorial
+ddu = [0.0, -0.00033, 0.0051, -0.0067, 0.0029, 0.0]
+du = [-0.047, -0.058, 0.054, 0.012, -0.068, 0.0011]
+A = QuinticHermiteSpline(ddu, du, u, t)
+scatter(t, u, label = "input data")
+plot!(A)
+```
+
 ## Regularization Smoothing
 
 Smoothing by regularization (a.k.a. ridge regression) finds a function ``\hat{u}``

src/DataInterpolations.jl

@@ -10,6 +10,7 @@ using ForwardDiff
 import FindFirstFunctions: searchsortedfirstcorrelated, searchsortedlastcorrelated,
                            bracketstrictlymontonic
 
+include("parameter_caches.jl")
 include("interpolation_caches.jl")
 include("interpolation_utils.jl")
 include("interpolation_methods.jl")
@@ -75,7 +76,8 @@ end
 
 export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation,
        AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline,
-       BSplineInterpolation, BSplineApprox
+       BSplineInterpolation, BSplineApprox, CubicHermiteSpline,
+       QuinticHermiteSpline
 
 # added for RegularizationSmooth, JJS 11/27/21
 ### Regularization data smoothing and interpolation

src/derivatives.jl

@@ -191,3 +191,26 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig
     end
     ducum * A.d * scale, idx
 end
+
+# Cubic Hermite Spline
+function _derivative(
+        A::CubicHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess)
+    Δt₀ = t - A.t[idx]
+    Δt₁ = t - A.t[idx + 1]
+    out = A.du[idx]
+    out += Δt₀ * (Δt₀ * A.p.c₂[idx] + 2(A.p.c₁[idx] + Δt₁ * A.p.c₂[idx]))
+    out, idx
+end
+
+# Quintic Hermite Spline
+function _derivative(
+        A::QuinticHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess)
+    Δt₀ = t - A.t[idx]
+    Δt₁ = t - A.t[idx + 1]
+    out = A.du[idx] + A.ddu[idx] * Δt₀
+    out += Δt₀^2 *
+           (3A.p.c₁[idx] + (3Δt₁ + Δt₀) * A.p.c₂[idx] + (3Δt₁^2 + Δt₀ * 2Δt₁) * A.p.c₃[idx])
+    out, idx
+end

src/integrals.jl

@@ -118,3 +118,28 @@ end
 function integral(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number)
     throw(IntegralNotFoundError())
 end
+
+# Cubic Hermite Spline
+function _integral(
+        A::CubicHermiteSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number)
+    Δt₀ = t - A.t[idx]
+    Δt₁ = t - A.t[idx + 1]
+    out = Δt₀ * (A.u[idx] + Δt₀ * A.du[idx] / 2)
+    p = A.p.c₁[idx] + Δt₁ * A.p.c₂[idx]
+    dp = A.p.c₂[idx]
+    out += Δt₀^3 / 3 * (p - dp * Δt₀ / 4)
+    out
+end
+
+# Quintic Hermite Spline
+function _integral(
+        A::QuinticHermiteSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number)
+    Δt₀ = t - A.t[idx]
+    Δt₁ = t - A.t[idx + 1]
+    out = Δt₀ * (A.u[idx] + A.du[idx] * Δt₀ / 2 + A.ddu[idx] * Δt₀^2 / 6)
+    p = A.p.c₁[idx] + A.p.c₂[idx] * Δt₁ + A.p.c₃[idx] * Δt₁^2
+    dp = A.p.c₂[idx] + 2A.p.c₃[idx] * Δt₁
+    ddp = 2A.p.c₃[idx]
+    out += Δt₀^4 / 4 * (p - Δt₀ / 5 * dp + Δt₀^2 / 30 * ddp)
+    out
+end

src/interpolation_caches.jl

@@ -611,3 +611,76 @@ function BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolate = false)
     c[2:(end - 1)] .= vec(P)
     BSplineApprox(u, t, d, h, p, k, c, pVecType, knotVecType, extrapolate)
 end
+
+"""
+    CubicHermiteSpline(du, u, t; extrapolate = false)
+
+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.
+
+## Arguments
+
+  - `du`: the derivative at the data points.
+  - `u`: data points.
+  - `t`: time points.
+
+## Keyword Arguments
+
+  - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+"""
+struct CubicHermiteSpline{uType, tType, duType, pType, T} <: AbstractInterpolation{T}
+    du::uType
+    u::uType
+    t::tType
+    p::CubicHermiteParameterCache{pType}
+    extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
+    function CubicHermiteSpline(du, u, t, p, extrapolate)
+        new{typeof(u), typeof(t), typeof(du), typeof(p.c₁), eltype(u)}(
+            du, u, t, p, extrapolate, Ref(1))
+    end
+end
+
+function CubicHermiteSpline(du, u, t; extrapolate = false)
+    @assert length(u)==length(du) "Length of `u` is not equal to length of `du`."
+    u, t = munge_data(u, t)
+    p = CubicHermiteParameterCache(du, u, t)
+    return CubicHermiteSpline(du, u, t, p, extrapolate)
+end
+
+"""
+    QuinticHermiteSpline(ddu, du, u, t; extrapolate = false)
+
+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.
+
+## Arguments
+
+  - `ddu`: the second derivative at the data points.
+  - `du`: the derivative at the data points.
+  - `u`: data points.
+  - `t`: time points.
+
+## Keyword Arguments
+
+  - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+"""
+struct QuinticHermiteSpline{uType, tType, duType, dduType, pType, T} <:
+       AbstractInterpolation{T}
+    ddu::uType
+    du::uType
+    u::uType
+    t::tType
+    p::QuinticHermiteParameterCache{pType}
+    extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
+    function QuinticHermiteSpline(ddu, du, u, t, p, extrapolate)
+        new{typeof(u), typeof(t), typeof(du), typeof(ddu), typeof(p.c₁), eltype(u)}(
+            ddu, du, u, t, p, extrapolate, Ref(1))
+    end
+end
+
+function QuinticHermiteSpline(ddu, du, u, t; extrapolate = false)
+    @assert length(u)==length(du)==length(ddu) "Length of `u` is not equal to length of `du` or `ddu`."
+    u, t = munge_data(u, t)
+    p = QuinticHermiteParameterCache(ddu, du, u, t)
+    return QuinticHermiteSpline(ddu, du, u, t, p, extrapolate)
+end

src/interpolation_methods.jl

@@ -201,3 +201,25 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i
     end
     ucum, idx
 end
+
+# Cubic Hermite Spline
+function _interpolate(
+        A::CubicHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess)
+    Δt₀ = t - A.t[idx]
+    Δt₁ = t - A.t[idx + 1]
+    out = A.u[idx] + Δt₀ * A.du[idx]
+    out += Δt₀^2 * (A.p.c₁[idx] + Δt₁ * A.p.c₂[idx])
+    out, idx
+end
+
+# Quintic Hermite Spline
+function _interpolate(
+        A::QuinticHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess)
+    Δt₀ = t - A.t[idx]
+    Δt₁ = t - A.t[idx + 1]
+    out = A.u[idx] + Δt₀ * (A.du[idx] + A.ddu[idx] * Δt₀ / 2)
+    out += Δt₀^3 * (A.p.c₁[idx] + Δt₁ * (A.p.c₂[idx] + A.p.c₃[idx] * Δt₁))
+    out, idx
+end

src/parameter_caches.jl

@@ -0,0 +1,49 @@
+struct CubicHermiteParameterCache{pType}
+    c₁::pType
+    c₂::pType
+end
+
+function CubicHermiteParameterCache(du, u, t)
+    parameters = CubicHermiteSplineParameters.(
+        Ref(du), Ref(u), Ref(t), 1:(length(t) - 1))
+    c₁, c₂ = collect.(eachrow(hcat(collect.(parameters)...)))
+    return CubicHermiteParameterCache(c₁, c₂)
+end
+
+function CubicHermiteSplineParameters(du, u, t, idx)
+    Δt = t[idx + 1] - t[idx]
+    u₀ = u[idx]
+    u₁ = u[idx + 1]
+    du₀ = du[idx]
+    du₁ = du[idx + 1]
+    c₁ = (u₁ - u₀ - du₀ * Δt) / Δt^2
+    c₂ = (du₁ - du₀ - 2c₁ * Δt) / Δt^2
+    return c₁, c₂
+end
+
+struct QuinticHermiteParameterCache{pType}
+    c₁::pType
+    c₂::pType
+    c₃::pType
+end
+
+function QuinticHermiteParameterCache(ddu, du, u, t)
+    parameters = QuinticHermiteSplineParameters.(
+        Ref(ddu), Ref(du), Ref(u), Ref(t), 1:(length(t) - 1))
+    c₁, c₂, c₃ = collect.(eachrow(hcat(collect.(parameters)...)))
+    return QuinticHermiteParameterCache(c₁, c₂, c₃)
+end
+
+function QuinticHermiteSplineParameters(ddu, du, u, t, idx)
+    Δt = t[idx + 1] - t[idx]
+    u₀ = u[idx]
+    u₁ = u[idx + 1]
+    du₀ = du[idx]
+    du₁ = du[idx + 1]
+    ddu₀ = ddu[idx]
+    ddu₁ = ddu[idx + 1]
+    c₁ = (u₁ - u₀ - du₀ * Δt - ddu₀ * Δt^2 / 2) / Δt^3
+    c₂ = (3u₀ - 3u₁ + 2(du₀ + du₁ / 2)Δt + ddu₀ * Δt^2 / 2) / Δt^4
+    c₃ = (6u₁ - 6u₀ - 3(du₀ + du₁)Δt + (ddu₁ - ddu₀)Δt^2 / 2) / Δt^5
+    return c₁, c₂, c₃
+end