Extrapolation options

Project.toml

@@ -3,6 +3,7 @@ uuid = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 version = "6.6.0"
 
 [deps]
+EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -27,6 +28,7 @@ DataInterpolationsSymbolicsExt = "Symbolics"
 Aqua = "0.8"
 BenchmarkTools = "1"
 ChainRulesCore = "1.24"
+EnumX = "1.0.4"
 FindFirstFunctions = "1.3"
 FiniteDifferences = "0.12.31"
 ForwardDiff = "0.10.36"

docs/make.jl

@@ -12,7 +12,8 @@ makedocs(modules = [DataInterpolations],
     linkcheck = true,
     format = Documenter.HTML(assets = ["assets/favicon.ico"],
         canonical = "https://docs.sciml.ai/DataInterpolations/stable/"),
-    pages = ["index.md", "Methods" => "methods.md",
+    pages = ["index.md", "Interpolation methods" => "methods.md",
+        "Extrapolation methods" => "extrapolation_methods.md",
         "Interface" => "interface.md", "Using with Symbolics/ModelingToolkit" => "symbolics.md",
         "Manual" => "manual.md", "Inverting Integrals" => "inverting_integrals.md"])
 

docs/src/extrapolation_methods.md

@@ -0,0 +1,56 @@
+# Extrapolation methods
+
+We will use the following interpolation to demonstrate the various extrapolation methods.
+
+```@example tutorial
+using DataInterpolations, Plots
+
+u = [0.86, 0.65, 0.44, 0.76, 0.73]
+t = [0.0022, 0.68, 1.41, 2.22, 2.46]
+t_eval_down = range(-1, first(t), length = 25)
+t_eval_up = range(last(t), 3.5, length = 25)
+A = QuadraticSpline(u, t)
+plot(A)
+```
+
+Extrapolation behavior can be set for `t` beyond the data in the negative and positive direction separately with the `extrapolation_down` and `extrapolation_up` keywords of the interpolation constructors respectively.
+
+## `ExtrapolationType.none`
+
+This extrapolation type will throw an error when the input `t` is beyond the data in the specified direction.
+
+## `ExtrapolationType.constant`
+
+This extrapolation type extends the interpolation with the boundary values of the data `u`.
+
+```@example tutorial
+A = QuadraticSpline(u, t; extrapolation_down = ExtrapolationType.constant,
+    extrapolation_up = ExtrapolationType.constant)
+plot(A)
+plot!(t_eval_down, A.(t_eval_down); label = "extrapolation down")
+plot!(t_eval_up, A.(t_eval_up); label = "extrapolation up")
+```
+
+## `ExtrapolationType.linear`
+
+This extrapolation type extends the interpolation with a linear continuation of the interpolation, making it $C^1$ smooth at the data boundaries.
+
+```@example tutorial
+A = QuadraticSpline(u, t; extrapolation_down = ExtrapolationType.linear,
+    extrapolation_up = ExtrapolationType.linear)
+plot(A)
+plot!(t_eval_down, A.(t_eval_down); label = "extrapolation down")
+plot!(t_eval_up, A.(t_eval_up); label = "extrapolation up")
+```
+
+## `ExtrapolationType.extension`
+
+This extrapolation type extends the interpolation with a continuation of the expression for the interpolation at the boundary intervals for maximum smoothness.
+
+```@example tutorial
+A = QuadraticSpline(u, t; extrapolation_down = ExtrapolationType.extension,
+    extrapolation_up = ExtrapolationType.extension)
+plot(A)
+plot!(t_eval_down, A.(t_eval_down); label = "extrapolation down")
+plot!(t_eval_up, A.(t_eval_up); label = "extrapolation up")
+```

docs/src/interface.md

@@ -17,15 +17,15 @@ t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3]
 
 All interpolation methods return an object from which we can compute the value of the dependent variable at any time point.
 
-We will use the `CubicSpline` method for demonstration, but the API is the same for all the methods. We can also pass the `extrapolate = true` keyword if we want to allow the interpolation to go beyond the range of the timepoints. The default value is `extrapolate = false`.
+We will use the `CubicSpline` method for demonstration, but the API is the same for all the methods. We can also pass the `extrapolation_up = ExtrapolationType.extension` keyword if we want to allow the interpolation to go beyond the range of the timepoints in the positive `t` direction. The default value is `extrapolation_up = ExtrapolationType.none`. For more information on extrapolation see [Extrapolation methods](extrapolation_methods.md).
 
 ```@example interface
 A1 = CubicSpline(u, t)
 
 # For interpolation do, A(t)
 A1(100.0)
 
-A2 = CubicSpline(u, t; extrapolate = true)
+A2 = CubicSpline(u, t; extrapolation_up = ExtrapolationType.extension)
 
 # Extrapolation
 A2(300.0)

ext/DataInterpolationsRegularizationToolsExt.jl

@@ -69,13 +69,17 @@ A = RegularizationSmooth(u, t, t̂, wls, wr, d; λ = 1.0, alg = :gcv_svd)
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
         wls::AbstractVector, wr::AbstractVector, d::Int = 2;
-        λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false)
+        λ::Real = 1.0, alg::Symbol = :gcv_svd,
+        extrapolation_down::ExtrapolationType.T = ExtrapolationType.none,
+        extrapolation_up::ExtrapolationType.T = ExtrapolationType.none)
     u, t = munge_data(u, t)
     M = _mapping_matrix(t̂, t)
     Wls½ = LA.diagm(sqrt.(wls))
     Wr½ = LA.diagm(sqrt.(wr))
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
-    RegularizationSmooth(u, û, t, t̂, wls, wr, d, λ, alg, Aitp, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_down, extrapolation_up)
+    RegularizationSmooth(
+        u, û, t, t̂, wls, wr, d, λ, alg, Aitp, extrapolation_down, extrapolation_up)
 end
 """
 Direct smoothing, no `t̂` or weights
@@ -86,14 +90,16 @@ A = RegularizationSmooth(u, t, d; λ = 1.0, alg = :gcv_svd, extrapolate = false)
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, d::Int = 2;
         λ::Real = 1.0,
-        alg::Symbol = :gcv_svd, extrapolate::Bool = false)
+        alg::Symbol = :gcv_svd, extrapolation_down::ExtrapolationType.T = ExtrapolationType.none,
+        extrapolation_up::ExtrapolationType.T = ExtrapolationType.none)
     u, t = munge_data(u, t)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
     Wls½ = Array{Float64}(LA.I, N, N)
     Wr½ = Array{Float64}(LA.I, N - d, N - d)
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_down, extrapolation_up)
     RegularizationSmooth(u,
         û,
         t,
@@ -104,7 +110,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, d::Int = 2;
         λ,
         alg,
         Aitp,
-        extrapolate)
+        extrapolation_down,
+        extrapolation_up)
 end
 """
 `t̂` provided, no weights
@@ -115,13 +122,15 @@ A = RegularizationSmooth(u, t, t̂, d; λ = 1.0, alg = :gcv_svd, extrapolate = f
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
         d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd,
-        extrapolate::Bool = false)
+        extrapolation_down::ExtrapolationType.T = ExtrapolationType.none,
+        extrapolation_up::ExtrapolationType.T = ExtrapolationType.none)
     u, t = munge_data(u, t)
     N, N̂ = length(t), length(t̂)
     M = _mapping_matrix(t̂, t)
     Wls½ = Array{Float64}(LA.I, N, N)
     Wr½ = Array{Float64}(LA.I, N̂ - d, N̂ - d)
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_down, extrapolation_up)
     RegularizationSmooth(u,
         û,
         t,
@@ -132,7 +141,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Abstrac
         λ,
         alg,
         Aitp,
-        extrapolate)
+        extrapolation_down,
+        extrapolation_up)
 end
 """
 `t̂` and `wls` provided
@@ -143,13 +153,15 @@ A = RegularizationSmooth(u, t, t̂, wls, d; λ = 1.0, alg = :gcv_svd, extrapolat
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
         wls::AbstractVector, d::Int = 2; λ::Real = 1.0,
-        alg::Symbol = :gcv_svd, extrapolate::Bool = false)
+        alg::Symbol = :gcv_svd, extrapolation_down::ExtrapolationType.T = ExtrapolationType.none,
+        extrapolation_up::ExtrapolationType.T = ExtrapolationType.none)
     u, t = munge_data(u, t)
     N, N̂ = length(t), length(t̂)
     M = _mapping_matrix(t̂, t)
     Wls½ = LA.diagm(sqrt.(wls))
     Wr½ = Array{Float64}(LA.I, N̂ - d, N̂ - d)
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_down, extrapolation_up)
     RegularizationSmooth(u,
         û,
         t,
@@ -160,7 +172,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Abstrac
         λ,
         alg,
         Aitp,
-        extrapolate)
+        extrapolation_down,
+        extrapolation_up)
 end
 """
 `wls` provided, no `t̂`
@@ -172,14 +185,16 @@ A = RegularizationSmooth(
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
         wls::AbstractVector, d::Int = 2; λ::Real = 1.0,
-        alg::Symbol = :gcv_svd, extrapolate::Bool = false)
+        alg::Symbol = :gcv_svd, extrapolation_down::ExtrapolationType.T = ExtrapolationType.none,
+        extrapolation_up::ExtrapolationType.T = ExtrapolationType.none)
     u, t = munge_data(u, t)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
     Wls½ = LA.diagm(sqrt.(wls))
     Wr½ = Array{Float64}(LA.I, N - d, N - d)
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_down, extrapolation_up)
     RegularizationSmooth(u,
         û,
         t,
@@ -190,7 +205,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing
         λ,
         alg,
         Aitp,
-        extrapolate)
+        extrapolation_down,
+        extrapolation_up)
 end
 """
 `wls` and `wr` provided, no `t̂`
@@ -202,14 +218,17 @@ A = RegularizationSmooth(
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
         wls::AbstractVector, wr::AbstractVector, d::Int = 2;
-        λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false)
+        λ::Real = 1.0, alg::Symbol = :gcv_svd,
+        extrapolation_down::ExtrapolationType.T = ExtrapolationType.none,
+        extrapolation_up::ExtrapolationType.T = ExtrapolationType.none)
     u, t = munge_data(u, t)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
     Wls½ = LA.diagm(sqrt.(wls))
     Wr½ = LA.diagm(sqrt.(wr))
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_down, extrapolation_up)
     RegularizationSmooth(u,
         û,
         t,
@@ -220,7 +239,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing
         λ,
         alg,
         Aitp,
-        extrapolate)
+        extrapolation_down,
+        extrapolation_up)
 end
 """
 Keyword provided for `wls`, no `t̂`
@@ -232,15 +252,17 @@ A = RegularizationSmooth(
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
         wls::Symbol, d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd,
-        extrapolate::Bool = false)
+        extrapolation_down::ExtrapolationType.T = ExtrapolationType.none,
+        extrapolation_up::ExtrapolationType.T = ExtrapolationType.none)
     u, t = munge_data(u, t)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
     wls, wr = _weighting_by_kw(t, d, wls)
     Wls½ = LA.diagm(sqrt.(wls))
     Wr½ = LA.diagm(sqrt.(wr))
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_down, extrapolation_up)
     RegularizationSmooth(u,
         û,
         t,
@@ -251,7 +273,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing
         λ,
         alg,
         Aitp,
-        extrapolate)
+        extrapolation_down,
+        extrapolation_up)
 end
 # """ t̂ provided and keyword for wls  _TBD_ """
 # function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
@@ -262,7 +285,8 @@ Solve for the smoothed dependent variables and create spline interpolator
 """
 function _reg_smooth_solve(
         u::AbstractVector, t̂::AbstractVector, d::Int, M::AbstractMatrix,
-        Wls½::AbstractMatrix, Wr½::AbstractMatrix, λ::Real, alg::Symbol, extrapolate::Bool)
+        Wls½::AbstractMatrix, Wr½::AbstractMatrix, λ::Real, alg::Symbol,
+        extrapolation_down::ExtrapolationType.T, extrapolation_up::ExtrapolationType.T)
     λ = float(λ) # `float` expected by RT
     D = _derivative_matrix(t̂, d)
     Ψ = RT.setupRegularizationProblem(Wls½ * M, Wr½ * D)
@@ -279,7 +303,7 @@ function _reg_smooth_solve(
         û = result.x
         λ = result.λ
     end
-    Aitp = CubicSpline(û, t̂; extrapolate)
+    Aitp = CubicSpline(û, t̂; extrapolation_down, extrapolation_up)
     # It seems logical to use B-Spline of order d+1, but I am unsure if theory supports the
     # extra computational cost, JJS 12/25/21
     #Aitp = BSplineInterpolation(û,t̂,d+1,:ArcLen,:Average)