Merge pull request #363 from SouthEndMusic/more_extrapolation_options

Add extrapolation types `periodic` and `reflective`
SciML · Dec 1, 2024 · 6f7a38a · 6f7a38a
2 parents 7df4954 + 120e56f
commit 6f7a38a
Show file tree

Hide file tree

Showing 22 changed files with 983 additions and 268 deletions.
diff --git a/Project.toml b/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"

diff --git a/docs/make.jl b/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"])
 

diff --git a/docs/src/extrapolation_methods.md b/docs/src/extrapolation_methods.md
@@ -0,0 +1,90 @@
+# 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_left = range(-1, first(t), length = 25)
+t_eval_right = range(last(t), 3.5, length = 25)
+A = QuadraticSpline(u, t)
+plot(A)
+```
+
+Extrapolation behavior can be set left and right of the data simultaneously with the `extension` keyword, or left and right separately with the `extrapolation_left` and `extrapolation_right` keywords 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 = ExtrapolationType.Constant)
+plot(A)
+plot!(t_eval_left, A.(t_eval_left); label = "extrapolation left")
+plot!(t_eval_right, A.(t_eval_right); label = "extrapolation right")
+```
+
+## `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 = ExtrapolationType.Linear)
+plot(A)
+plot!(t_eval_left, A.(t_eval_left); label = "extrapolation left")
+plot!(t_eval_right, A.(t_eval_right); label = "extrapolation right")
+```
+
+## `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 = ExtrapolationType.Extension)
+plot(A)
+plot!(t_eval_left, A.(t_eval_left); label = "extrapolation down")
+plot!(t_eval_right, A.(t_eval_right); label = "extrapolation up")
+```
+
+## `ExtrapolationType.Periodic`
+
+this extrapolation type extends the interpolation such that `A(t + T) == A(t)` for all `t`, where the period is given by `T = last(A.t) - first(A.t)`.
+
+```@example tutorial
+T = last(A.t) - first(A.t)
+t_eval_left = range(first(t) - 2T, first(t), length = 100)
+t_eval_right = range(last(t), last(t) + 2T, length = 100)
+A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.Periodic)
+plot(A)
+plot!(t_eval_left, A.(t_eval_left); label = "extrapolation down")
+plot!(t_eval_right, A.(t_eval_right); label = "extrapolation up")
+```
+
+## `ExtrapolationType.Reflective`
+
+this extrapolation type extends the interpolation such that `A(t_ + t) == A(t_ - t)` for all `t_, t` such that `(t_ - first(A.t)) % T == 0` and `0 < t < T`, where `T = last(A.t) - first(A.t)`.
+
+```@example tutorial
+A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.Reflective)
+plot(A)
+plot!(t_eval_left, A.(t_eval_left); label = "extrapolation down")
+plot!(t_eval_right, A.(t_eval_right); label = "extrapolation up")
+```
+
+## Mixed extrapolation
+
+You can also have different extrapolation types left and right of the data.
+
+```@example tutorial
+A = QuadraticSpline(u, t; extrapolation_left = ExtrapolationType.Reflective,
+    extrapolation_right = ExtrapolationType.Periodic)
+plot(A)
+plot!(t_eval_left, A.(t_eval_left); label = "extrapolation left")
+plot!(t_eval_right, A.(t_eval_right); label = "extrapolation right")
+```
diff --git a/docs/src/interface.md b/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 = 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 = 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 = ExtrapolationType.Extension)
 
 # Extrapolation
 A2(300.0)

diff --git a/ext/DataInterpolationsRegularizationToolsExt.jl b/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_left::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_right::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))
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
-    RegularizationSmooth(u, û, t, t̂, wls, wr, d, λ, alg, Aitp, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
+    RegularizationSmooth(
+        u, û, t, t̂, wls, wr, d, λ, alg, Aitp, extrapolation_left, extrapolation_right)
 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_left::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_right::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)
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
     RegularizationSmooth(u,
         û,
         t,
@@ -104,7 +110,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, d::Int = 2;
         λ,
         alg,
         Aitp,
-        extrapolate)
+        extrapolation_left,
+        extrapolation_right)
 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_left::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_right::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)
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
     RegularizationSmooth(u,
         û,
         t,
@@ -132,7 +141,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Abstrac
         λ,
         alg,
         Aitp,
-        extrapolate)
+        extrapolation_left,
+        extrapolation_right)
 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_left::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_right::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)
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
     RegularizationSmooth(u,
         û,
         t,
@@ -160,7 +172,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Abstrac
         λ,
         alg,
         Aitp,
-        extrapolate)
+        extrapolation_left,
+        extrapolation_right)
 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_left::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_right::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)
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
     RegularizationSmooth(u,
         û,
         t,
@@ -190,7 +205,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing
         λ,
         alg,
         Aitp,
-        extrapolate)
+        extrapolation_left,
+        extrapolation_right)
 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_left::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_right::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))
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
     RegularizationSmooth(u,
         û,
         t,
@@ -220,7 +239,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing
         λ,
         alg,
         Aitp,
-        extrapolate)
+        extrapolation_left,
+        extrapolation_right)
 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_left::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_right::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))
-    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate)
+    û, λ, Aitp = _reg_smooth_solve(
+        u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
     RegularizationSmooth(u,
         û,
         t,
@@ -251,7 +273,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing
         λ,
         alg,
         Aitp,
-        extrapolate)
+        extrapolation_left,
+        extrapolation_right)
 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_left::ExtrapolationType.T, extrapolation_right::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(
         û = result.x
         λ = result.λ
     end
-    Aitp = CubicSpline(û, t̂; extrapolate)
+    Aitp = CubicSpline(û, t̂; extrapolation_left, extrapolation_right)
     # 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(û,t̂,d+1,:ArcLen,:Average)