diff --git a/Project.toml b/Project.toml index 80cd1153..48e92f6f 100644 --- 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 index 6438f482..a6b062b5 100644 --- 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 new file mode 100644 index 00000000..c442a7f9 --- /dev/null +++ b/docs/src/extrapolation_methods.md @@ -0,0 +1,65 @@ +# 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 left") +plot!(t_eval_right, A.(t_eval_right); label = "extrapolation right") +``` + +## Mixed extrapolation + +You can also have different extrapolation types left and right of the data. + +```@example tutorial +A = QuadraticSpline(u, t; extrapolation_left = ExtrapolationType.constant, + extrapolation_right = 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") +``` diff --git a/docs/src/interface.md b/docs/src/interface.md index cfc9ed0b..21bfd27a 100644 --- a/docs/src/interface.md +++ b/docs/src/interface.md @@ -17,7 +17,7 @@ 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) @@ -25,7 +25,7 @@ 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 index 732ea1bb..84693c49 100644 --- 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)) - û, λ, 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_left, extrapolation_right) + RegularizationSmooth( + 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) - û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate) + û, λ, Aitp = _reg_smooth_solve( + u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right) RegularizationSmooth(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) - û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate) + û, λ, Aitp = _reg_smooth_solve( + u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right) RegularizationSmooth(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) - û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate) + û, λ, Aitp = _reg_smooth_solve( + u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right) RegularizationSmooth(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) - û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate) + û, λ, Aitp = _reg_smooth_solve( + u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right) RegularizationSmooth(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)) - û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolate) + û, λ, Aitp = _reg_smooth_solve( + u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right) RegularizationSmooth(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,7 +252,8 @@ 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) @@ -240,7 +261,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing 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_left, extrapolation_right) RegularizationSmooth(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( û = result.x λ = result.λ end - Aitp = CubicSpline(û, t̂; extrapolate) + Aitp = CubicSpline(û, 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(û,t̂,d+1,:ArcLen,:Average) diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl index 393f1160..962b21a9 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -7,9 +7,12 @@ abstract type AbstractInterpolation{T, N} end using LinearAlgebra, RecipesBase using PrettyTables using ForwardDiff +using EnumX import FindFirstFunctions: searchsortedfirstcorrelated, searchsortedlastcorrelated, Guesser +@enumx ExtrapolationType none constant linear extension + include("parameter_caches.jl") include("interpolation_caches.jl") include("interpolation_utils.jl") @@ -57,38 +60,51 @@ end const EXTRAPOLATION_ERROR = "Cannot extrapolate as `extrapolate` keyword passed was `false`" struct ExtrapolationError <: Exception end -function Base.showerror(io::IO, e::ExtrapolationError) +function Base.showerror(io::IO, ::ExtrapolationError) print(io, EXTRAPOLATION_ERROR) end +const LEFT_EXTRAPOLATION_ERROR = "Cannot extrapolate for t < first(A.t) as the `extrapolation_left` kwarg passed was `ExtrapolationType.none`" +struct LeftExtrapolationError <: Exception end +function Base.showerror(io::IO, ::LeftExtrapolationError) + print(io, LEFT_EXTRAPOLATION_ERROR) +end + +const RIGHT_EXTRAPOLATION_ERROR = "Cannot extrapolate for t > last(A.t) as the `extrapolation_tight` kwarg passed was `ExtrapolationType.none`" +struct RightExtrapolationError <: Exception end +function Base.showerror(io::IO, ::RightExtrapolationError) + print(io, RIGHT_EXTRAPOLATION_ERROR) +end + const INTEGRAL_NOT_FOUND_ERROR = "Cannot integrate it analytically. Please use Numerical Integration methods." struct IntegralNotFoundError <: Exception end -function Base.showerror(io::IO, e::IntegralNotFoundError) +function Base.showerror(io::IO, ::IntegralNotFoundError) print(io, INTEGRAL_NOT_FOUND_ERROR) end const DERIVATIVE_NOT_FOUND_ERROR = "Derivatives greater than second order is not supported." struct DerivativeNotFoundError <: Exception end -function Base.showerror(io::IO, e::DerivativeNotFoundError) +function Base.showerror(io::IO, ::DerivativeNotFoundError) print(io, DERIVATIVE_NOT_FOUND_ERROR) end const INTEGRAL_INVERSE_NOT_FOUND_ERROR = "Cannot invert the integral analytically. Please use Numerical methods." struct IntegralInverseNotFoundError <: Exception end -function Base.showerror(io::IO, e::IntegralInverseNotFoundError) +function Base.showerror(io::IO, ::IntegralInverseNotFoundError) print(io, INTEGRAL_INVERSE_NOT_FOUND_ERROR) end const INTEGRAL_NOT_INVERTIBLE_ERROR = "The Interpolation is not positive everywhere so its integral is not invertible." struct IntegralNotInvertibleError <: Exception end -function Base.showerror(io::IO, e::IntegralNotInvertibleError) +function Base.showerror(io::IO, ::IntegralNotInvertibleError) print(io, INTEGRAL_NOT_INVERTIBLE_ERROR) end export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline, BSplineInterpolation, BSplineApprox, CubicHermiteSpline, PCHIPInterpolation, - QuinticHermiteSpline, LinearInterpolationIntInv, ConstantInterpolationIntInv + QuinticHermiteSpline, LinearInterpolationIntInv, ConstantInterpolationIntInv, + ExtrapolationType # added for RegularizationSmooth, JJS 11/27/21 ### Regularization data smoothing and interpolation @@ -104,7 +120,8 @@ struct RegularizationSmooth{uType, tType, T, T2, N, ITP <: AbstractInterpolation λ::T2 # regularization parameter alg::Symbol # how to determine λ: `:fixed`, `:gcv_svd`, `:gcv_tr`, `L_curve` Aitp::ITP - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T function RegularizationSmooth(u, û, t, @@ -115,7 +132,8 @@ struct RegularizationSmooth{uType, tType, T, T2, N, ITP <: AbstractInterpolation λ, alg, Aitp, - extrapolate) + extrapolation_left, + extrapolation_right) N = get_output_dim(u) new{typeof(u), typeof(t), eltype(u), typeof(λ), N, typeof(Aitp)}( u, @@ -128,7 +146,8 @@ struct RegularizationSmooth{uType, tType, T, T2, N, ITP <: AbstractInterpolation λ, alg, Aitp, - extrapolate) + extrapolation_left, + extrapolation_right) end end diff --git a/src/derivatives.jl b/src/derivatives.jl index 7d6dea24..30994075 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -1,15 +1,51 @@ function derivative(A, t, order = 1) - ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) - iguess = A.iguesser + (order ∉ (1, 2)) && throw(DerivativeNotFoundError()) + if t < first(A.t) + _extrapolate_derivative_left(A, t, order) + elseif t > last(A.t) + _extrapolate_derivative_right(A, t, order) + else + iguess = A.iguesser + (order == 1) ? _derivative(A, t, iguess) : + ForwardDiff.derivative(t -> begin + _derivative(A, t, iguess) + end, t) + end +end - return if order == 1 - _derivative(A, t, iguess) - elseif order == 2 +function _extrapolate_derivative_left(A, t, order) + (; extrapolation_left) = A + if extrapolation_left == ExtrapolationType.none + throw(LeftExtrapolationError()) + elseif extrapolation_left == ExtrapolationType.constant + zero(first(A.u) / one(A.t[1])) + elseif extrapolation_left == ExtrapolationType.linear + (order == 1) ? derivative(A, first(A.t)) : zero(first(A.u) / one(A.t[1])) + else + # extrapolation_left == ExtrapolationType.extension + iguess = A.iguesser + (order == 1) ? _derivative(A, t, iguess) : ForwardDiff.derivative(t -> begin _derivative(A, t, iguess) end, t) + end +end + +function _extrapolate_derivative_right(A, t, order) + (; extrapolation_right) = A + if extrapolation_right == ExtrapolationType.none + throw(RightExtrapolationError()) + elseif extrapolation_right == ExtrapolationType.constant + zero(first(A.u) / one(A.t[1])) + elseif extrapolation_right == ExtrapolationType.linear + (order == 1) ? derivative(A, last(A.t)) : zero(first(A.u) / one(A.t[1])) else - throw(DerivativeNotFoundError()) + # extrapolation_right == ExtrapolationType.extension + iguess = A.iguesser + (order == 1) ? _derivative(A, t, iguess) : + ForwardDiff.derivative(t -> begin + _derivative(A, t, iguess) + end, t) end end @@ -27,7 +63,6 @@ function _derivative(A::QuadraticInterpolation, t::Number, iguess) end function _derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number) - ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) der = zero(A.u[1]) for j in eachindex(A.t) tmp = zero(A.t[1]) @@ -61,7 +96,6 @@ function _derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number) end function _derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number) - ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) der = zero(A.u[:, 1]) for j in eachindex(A.t) tmp = zero(A.t[1]) @@ -113,12 +147,10 @@ function _derivative(A::ConstantInterpolation, t::Number, iguess) end function _derivative(A::ConstantInterpolation{<:AbstractVector}, t::Number, iguess) - ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) return isempty(searchsorted(A.t, t)) ? zero(A.u[1]) : eltype(A.u)(NaN) end function _derivative(A::ConstantInterpolation{<:AbstractMatrix}, t::Number, iguess) - ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) return isempty(searchsorted(A.t, t)) ? zero(A.u[:, 1]) : eltype(A.u)(NaN) .* A.u[:, 1] end diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index ddcac9db..577b98e3 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -37,13 +37,14 @@ struct LinearInterpolationIntInv{uType, tType, itpType, T, N} <: AbstractIntegralInverseInterpolation{T, N} u::uType t::tType - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} itp::itpType function LinearInterpolationIntInv(u, t, A) N = get_output_dim(u) new{typeof(u), typeof(t), typeof(A), eltype(u), N}( - u, t, A.extrapolate, Guesser(t), A) + u, t, A.extrapolation_left, A.extrapolation_right, Guesser(t), A) end end @@ -88,13 +89,14 @@ struct ConstantInterpolationIntInv{uType, tType, itpType, T, N} <: AbstractIntegralInverseInterpolation{T, N} u::uType t::tType - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} itp::itpType function ConstantInterpolationIntInv(u, t, A) N = get_output_dim(u) new{typeof(u), typeof(t), typeof(A), eltype(u), N}( - u, t, A.extrapolate, Guesser(t), A + u, t, A.extrapolation_left, A.extrapolation_right, Guesser(t), A ) end end diff --git a/src/integrals.jl b/src/integrals.jl index 2bdb6408..1da09b83 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -1,39 +1,94 @@ function integral(A::AbstractInterpolation, t::Number) - ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) - integral(A, A.t[1], t) + integral(A, first(A.t), t) end function integral(A::AbstractInterpolation, t1::Number, t2::Number) - ((t1 < A.t[1] || t1 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) - ((t2 < A.t[1] || t2 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) !hasfield(typeof(A), :I) && throw(IntegralNotFoundError()) - (t2 < t1) && return -integral(A, t2, t1) + + if t1 == t2 + # If the integration interval is trivial then the result is 0 + return zero(eltype(A.I)) + elseif t1 > t2 + # Make sure that t1 < t2 + return -integral(A, t2, t1) + end + # the index less than or equal to t1 idx1 = get_idx(A, t1, 0) # the index less than t2 idx2 = get_idx(A, t2, 0; idx_shift = -1, side = :first) - if A.cache_parameters - total = A.I[max(1, idx2 - 1)] - A.I[idx1] - return if t1 == t2 - zero(total) - else - if idx1 == idx2 - total += _integral(A, idx1, t1, t2) - else - total += _integral(A, idx1, t1, A.t[idx1 + 1]) - total += _integral(A, idx2, A.t[idx2], t2) - end - total + total = zero(eltype(A.I)) + + # Lower potentially incomplete interval + if t1 < first(A.t) + if t2 < first(A.t) + # If interval is entirely below data + return _extrapolate_integral_left(A, t2) - extrapolate_integral_left(A.t1) + end + + idx1 -= 1 # Make sure lowest complete interval is included + total += _extrapolate_integral_left(A, t1) + else + total += _integral(A, idx1, t1, A.t[idx1 + 1]) + end + + # Upper potentially incomplete interval + if t2 > last(A.t) + if t1 > last(A.t) + # If interval is entirely above data + return _extrapolate_integral_right(A, t2) - extrapolate_integral_right(A.t, t1) end + + idx2 += 1 # Make sure highest complete interval is included + total += _extrapolate_integral_right(A, t2) + else + total += _integral(A, idx2, A.t[idx2], t2) + end + + if idx1 == idx2 + return _integral(A, idx1, t1, t2) + end + + # Complete intervals + if A.cache_parameters + total += A.I[idx2] - A.I[idx1 + 1] else - total = zero(eltype(A.u)) - for idx in idx1:idx2 - lt1 = idx == idx1 ? t1 : A.t[idx] - lt2 = idx == idx2 ? t2 : A.t[idx + 1] - total += _integral(A, idx, lt1, lt2) + for idx in (idx1 + 1):(idx2 - 1) + total += _integral(A, idx, A.t[idx], A.t[idx + 1]) end - total + end + + return total +end + +function _extrapolate_integral_left(A, t) + (; extrapolation_left) = A + if extrapolation_left == ExtrapolationType.none + throw(LeftExtrapolationError()) + elseif extrapolation_left == ExtrapolationType.constant + first(A.u) * (first(A.t) - t) + elseif extrapolation_left == ExtrapolationType.linear + slope = derivative(A, first(A.t)) + Δt = first(A.t) - t + (first(A.u) - slope * Δt / 2) * Δt + elseif extrapolation_left == ExtrapolationType.extension + _integral(A, 1, t, first(A.t)) + end +end + +function _extrapolate_integral_right(A, t) + (; extrapolation_right) = A + if extrapolation_right == ExtrapolationType.none + throw(RightExtrapolationError()) + elseif extrapolation_right == ExtrapolationType.constant + last(A.u) * (t - last(A.t)) + elseif extrapolation_right == ExtrapolationType.linear + slope = derivative(A, last(A.t)) + Δt = t - last(A.t) + (last(A.u) + slope * Δt / 2) * Δt + elseif extrapolation_right == ExtrapolationType.extension + _integral(A, length(A.t) - 1, last(A.t), t) end end diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 850ac74b..faeec06a 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -1,5 +1,7 @@ """ - LinearInterpolation(u, t; extrapolate = false, cache_parameters = false) + LinearInterpolation(u, t; extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation::ExtrapolationType.T = ExtrapolationType.none, extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, + cache_parameters = false) 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. @@ -11,10 +13,16 @@ Extrapolation extends the last linear polynomial on each side. ## Keyword Arguments - - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `extrapolation`: The extrapolation type applied left and right of the data. Possible options + are `ExtrapolationType.none` (default), `ExtrapolationType.constant`, + `ExtrapolationType.linear` and `ExtrapolationType.extension`. + - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. + - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`. - - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for + - `assume_linear_t`: boolean value to specify a faster index lookup behavior for evenly-distributed abscissae. Alternatively, a numerical threshold may be specified for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. @@ -24,30 +32,42 @@ struct LinearInterpolation{uType, tType, IType, pType, T, N} <: AbstractInterpol t::tType I::IType p::LinearParameterCache{pType} - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool - function LinearInterpolation(u, t, I, p, extrapolate, cache_parameters, assume_linear_t) + function LinearInterpolation(u, t, I, p, extrapolation_left, extrapolation_right, + cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) N = get_output_dim(u) new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u), N}( - u, t, I, p, extrapolate, Guesser(t), cache_parameters, linear_lookup) + u, t, I, p, extrapolation_left, extrapolation_right, + Guesser(t), cache_parameters, linear_lookup) end end function LinearInterpolation( - u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2) + u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false, assume_linear_t = 1e-2) + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) p = LinearParameterCache(u, t, cache_parameters) A = LinearInterpolation( - u, t, nothing, p, extrapolate, cache_parameters, assume_linear_t) + u, t, nothing, p, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) I = cumulative_integral(A, cache_parameters) - LinearInterpolation(u, t, I, p, extrapolate, cache_parameters, assume_linear_t) + LinearInterpolation( + u, t, I, p, extrapolation_left, extrapolation_right, + cache_parameters, assume_linear_t) end """ - QuadraticInterpolation(u, t, mode = :Forward; extrapolate = false, cache_parameters = false) + QuadraticInterpolation(u, t, mode = :Forward; extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation::ExtrapolationType.T = ExtrapolationType.none, extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, + cache_parameters = false) 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. @@ -60,7 +80,13 @@ Extrapolation extends the last quadratic polynomial on each side. ## Keyword Arguments - - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `extrapolation`: The extrapolation type applied left and right of the data. Possible options + are `ExtrapolationType.none` (default), `ExtrapolationType.constant`, + `ExtrapolationType.linear` and `ExtrapolationType.extension`. + - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. + - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`. - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for evenly-distributed abscissae. Alternatively, a numerical threshold may be specified @@ -74,30 +100,39 @@ struct QuadraticInterpolation{uType, tType, IType, pType, T, N} <: I::IType p::QuadraticParameterCache{pType} mode::Symbol - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function QuadraticInterpolation( - u, t, I, p, mode, extrapolate, cache_parameters, assume_linear_t) + u, t, I, p, mode, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) mode ∈ (:Forward, :Backward) || error("mode should be :Forward or :Backward for QuadraticInterpolation") linear_lookup = seems_linear(assume_linear_t, t) N = get_output_dim(u) new{typeof(u), typeof(t), typeof(I), typeof(p.α), eltype(u), N}( - u, t, I, p, mode, extrapolate, Guesser(t), cache_parameters, linear_lookup) + u, t, I, p, mode, extrapolation_left, extrapolation_right, + Guesser(t), cache_parameters, linear_lookup) end end function QuadraticInterpolation( - u, t, mode; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2) + u, t, mode; extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false, assume_linear_t = 1e-2) + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) linear_lookup = seems_linear(assume_linear_t, t) p = QuadraticParameterCache(u, t, cache_parameters, mode) A = QuadraticInterpolation( - u, t, nothing, p, mode, extrapolate, cache_parameters, linear_lookup) + u, t, nothing, p, mode, extrapolation_left, + extrapolation_right, cache_parameters, linear_lookup) I = cumulative_integral(A, cache_parameters) - QuadraticInterpolation(u, t, I, p, mode, extrapolate, cache_parameters, linear_lookup) + QuadraticInterpolation(u, t, I, p, mode, extrapolation_left, + extrapolation_right, cache_parameters, linear_lookup) end function QuadraticInterpolation(u, t; kwargs...) @@ -105,7 +140,8 @@ function QuadraticInterpolation(u, t; kwargs...) end """ - LagrangeInterpolation(u, t, n = length(t) - 1; extrapolate = false, safetycopy = true) + LagrangeInterpolation(u, t, n = length(t) - 1; extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, extrapolation_right::ExtrapolationType.T = ExtrapolationType.none) 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. @@ -117,7 +153,13 @@ It is the method of interpolation using Lagrange polynomials of (k-1)th order pa ## Keyword Arguments - - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `extrapolation`: The extrapolation type applied left and right of the data. Possible options + are `ExtrapolationType.none` (default), `ExtrapolationType.constant`, + `ExtrapolationType.linear` and `ExtrapolationType.extension`. + - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. + - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. """ struct LagrangeInterpolation{uType, tType, T, bcacheType, N} <: AbstractInterpolation{T, N} @@ -126,9 +168,10 @@ struct LagrangeInterpolation{uType, tType, T, bcacheType, N} <: n::Int bcache::bcacheType idxs::Vector{Int} - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} - function LagrangeInterpolation(u, t, n, extrapolate) + function LagrangeInterpolation(u, t, n, extrapolation_left, extrapolation_right) bcache = zeros(eltype(u[1]), n + 1) idxs = zeros(Int, n + 1) fill!(bcache, NaN) @@ -138,23 +181,30 @@ struct LagrangeInterpolation{uType, tType, T, bcacheType, N} <: n, bcache, idxs, - extrapolate, + extrapolation_left, + extrapolation_right, Guesser(t) ) end end function LagrangeInterpolation( - u, t, n = length(t) - 1; extrapolate = false) + u, t, n = length(t) - 1; + extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none) + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) if n != length(t) - 1 error("Currently only n=length(t) - 1 is supported") end - LagrangeInterpolation(u, t, n, extrapolate) + LagrangeInterpolation(u, t, n, extrapolation_left, extrapolation_right) end """ - AkimaInterpolation(u, t; extrapolate = false, cache_parameters = false) + AkimaInterpolation(u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.none, extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false) 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. @@ -166,7 +216,13 @@ Extrapolation extends the last cubic polynomial on each side. ## Keyword Arguments - - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `extrapolation`: The extrapolation type applied left and right of the data. Possible options + are `ExtrapolationType.none` (default), `ExtrapolationType.constant`, + `ExtrapolationType.linear` and `ExtrapolationType.extension`. + - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. + - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`. - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for evenly-distributed abscissae. Alternatively, a numerical threshold may be specified @@ -181,12 +237,14 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T, N} <: b::bType c::cType d::dType - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function AkimaInterpolation( - u, t, I, b, c, d, extrapolate, cache_parameters, assume_linear_t) + u, t, I, b, c, d, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) N = get_output_dim(u) new{typeof(u), typeof(t), typeof(I), typeof(b), typeof(c), @@ -196,7 +254,8 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T, N} <: b, c, d, - extrapolate, + extrapolation_left, + extrapolation_right, Guesser(t), cache_parameters, linear_lookup @@ -205,7 +264,11 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T, N} <: end function AkimaInterpolation( - u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2) + u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false, assume_linear_t = 1e-2) + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) linear_lookup = seems_linear(assume_linear_t, t) n = length(t) @@ -229,13 +292,16 @@ function AkimaInterpolation( d = (b[1:(end - 1)] .+ b[2:end] .- 2.0 .* m[3:(end - 2)]) ./ dt .^ 2 A = AkimaInterpolation( - u, t, nothing, b, c, d, extrapolate, cache_parameters, linear_lookup) + u, t, nothing, b, c, d, extrapolation_left, + extrapolation_right, cache_parameters, linear_lookup) I = cumulative_integral(A, cache_parameters) - AkimaInterpolation(u, t, I, b, c, d, extrapolate, cache_parameters, linear_lookup) + AkimaInterpolation(u, t, I, b, c, d, extrapolation_left, + extrapolation_right, cache_parameters, linear_lookup) end """ - ConstantInterpolation(u, t; dir = :left, extrapolate = false, cache_parameters = false) + ConstantInterpolation(u, t; dir = :left, extrapolation::ExtrapolationType.T = ExtrapolationType.none, extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false) 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. @@ -249,7 +315,13 @@ Extrapolation extends the last constant polynomial at the end points on each sid ## Keyword Arguments - `dir`: indicates which value should be used for interpolation (`:left` or `:right`). - - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `extrapolation`: The extrapolation type applied left and right of the data. Possible options + are `ExtrapolationType.none` (default), `ExtrapolationType.constant`, + `ExtrapolationType.linear` and `ExtrapolationType.extension`. + - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. + - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`. - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for evenly-distributed abscissae. Alternatively, a numerical threshold may be specified @@ -262,31 +334,41 @@ struct ConstantInterpolation{uType, tType, IType, T, N} <: AbstractInterpolation I::IType p::Nothing dir::Symbol # indicates if value to the $dir should be used for the interpolation - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function ConstantInterpolation( - u, t, I, dir, extrapolate, cache_parameters, assume_linear_t) + u, t, I, dir, extrapolation_left, extrapolation_right, + cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) N = get_output_dim(u) new{typeof(u), typeof(t), typeof(I), eltype(u), N}( - u, t, I, nothing, dir, extrapolate, Guesser(t), cache_parameters, linear_lookup) + u, t, I, nothing, dir, extrapolation_left, extrapolation_right, + Guesser(t), cache_parameters, linear_lookup) end end function ConstantInterpolation( - u, t; dir = :left, extrapolate = false, + u, t; dir = :left, extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false, assume_linear_t = 1e-2) + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) A = ConstantInterpolation( - u, t, nothing, dir, extrapolate, cache_parameters, assume_linear_t) + u, t, nothing, dir, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) I = cumulative_integral(A, cache_parameters) - ConstantInterpolation(u, t, I, dir, extrapolate, cache_parameters, assume_linear_t) + ConstantInterpolation(u, t, I, dir, extrapolation_left, extrapolation_right, + cache_parameters, assume_linear_t) end """ - QuadraticSpline(u, t; extrapolate = false, cache_parameters = false) + QuadraticSpline(u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.none, extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false) 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. @@ -298,7 +380,13 @@ Extrapolation extends the last quadratic polynomial on each side. ## Keyword Arguments - - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `extrapolation`: The extrapolation type applied left and right of the data. Possible options + are `ExtrapolationType.none` (default), `ExtrapolationType.constant`, + `ExtrapolationType.linear` and `ExtrapolationType.extension`. + - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. + - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`. - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for evenly-distributed abscissae. Alternatively, a numerical threshold may be specified @@ -314,12 +402,14 @@ struct QuadraticSpline{uType, tType, IType, pType, kType, cType, scType, T, N} < k::kType # knot vector c::cType # B-spline control points sc::scType # Spline coefficients (preallocated memory) - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function QuadraticSpline( - u, t, I, p, k, c, sc, extrapolate, cache_parameters, assume_linear_t) + u, t, I, p, k, c, sc, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) N = get_output_dim(u) new{typeof(u), typeof(t), typeof(I), typeof(p.α), typeof(k), @@ -330,7 +420,8 @@ struct QuadraticSpline{uType, tType, IType, pType, kType, cType, scType, T, N} < k, c, sc, - extrapolate, + extrapolation_left, + extrapolation_right, Guesser(t), cache_parameters, linear_lookup @@ -339,9 +430,13 @@ struct QuadraticSpline{uType, tType, IType, pType, kType, cType, scType, T, N} < end function QuadraticSpline( - u::uType, t; extrapolate = false, + u::uType, t; extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false, assume_linear_t = 1e-2) where {uType <: AbstractVector{<:Number}} + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) @@ -352,15 +447,21 @@ function QuadraticSpline( p = QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters) A = QuadraticSpline( - u, t, nothing, p, k, c, sc, extrapolate, cache_parameters, assume_linear_t) + u, t, nothing, p, k, c, sc, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) I = cumulative_integral(A, cache_parameters) - QuadraticSpline(u, t, I, p, k, c, sc, extrapolate, cache_parameters, assume_linear_t) + QuadraticSpline(u, t, I, p, k, c, sc, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) end function QuadraticSpline( - u::uType, t; extrapolate = false, cache_parameters = false, + u::uType, t; extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false, assume_linear_t = 1e-2) where {uType <: AbstractVector} + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) @@ -381,13 +482,16 @@ function QuadraticSpline( p = QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters) A = QuadraticSpline( - u, t, nothing, p, k, c, sc, extrapolate, cache_parameters, assume_linear_t) + u, t, nothing, p, k, c, sc, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) I = cumulative_integral(A, cache_parameters) - QuadraticSpline(u, t, I, p, k, c, sc, extrapolate, cache_parameters, assume_linear_t) + QuadraticSpline(u, t, I, p, k, c, sc, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) end """ - CubicSpline(u, t; extrapolate = false, cache_parameters = false) + CubicSpline(u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.none, extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false) 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. @@ -399,7 +503,13 @@ Second derivative on both ends are zero, which are also called "natural" boundar ## Keyword Arguments - - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `extrapolation`: The extrapolation type applied left and right of the data. Possible options + are `ExtrapolationType.none` (default), `ExtrapolationType.constant`, + `ExtrapolationType.linear` and `ExtrapolationType.extension`. + - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. + - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`. - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for evenly-distributed abscissae. Alternatively, a numerical threshold may be specified @@ -414,11 +524,13 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T, N} <: p::CubicSplineParameterCache{pType} h::hType z::zType - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool - function CubicSpline(u, t, I, p, h, z, extrapolate, cache_parameters, assume_linear_t) + function CubicSpline(u, t, I, p, h, z, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) N = get_output_dim(u) new{typeof(u), typeof(t), typeof(I), typeof(p.c₁), @@ -429,7 +541,8 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T, N} <: p, h, z, - extrapolate, + extrapolation_left, + extrapolation_right, Guesser(t), cache_parameters, linear_lookup @@ -438,10 +551,13 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T, N} <: end function CubicSpline(u::uType, - t; - extrapolate = false, cache_parameters = false, + t; extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false, assume_linear_t = 1e-2) where {uType <: AbstractVector{<:Number}} + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) - 1 h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0) @@ -462,16 +578,21 @@ function CubicSpline(u::uType, linear_lookup = seems_linear(assume_linear_t, t) p = CubicSplineParameterCache(u, h, z, cache_parameters) A = CubicSpline( - u, t, nothing, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup) + u, t, nothing, p, h[1:(n + 1)], z, extrapolation_left, + extrapolation_right, cache_parameters, linear_lookup) I = cumulative_integral(A, cache_parameters) - CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup) + CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolation_left, + extrapolation_right, cache_parameters, linear_lookup) end function CubicSpline(u::uType, t; - extrapolate = false, cache_parameters = false, + extrapolation::ExtrapolationType.T = ExtrapolationType.none, extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false, assume_linear_t = 1e-2) where {uType <: AbstractArray{T, N}} where {T, N} + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) - 1 h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0) @@ -496,15 +617,21 @@ function CubicSpline(u::uType, linear_lookup = seems_linear(assume_linear_t, t) p = CubicSplineParameterCache(u, h, z, cache_parameters) A = CubicSpline( - u, t, nothing, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup) + u, t, nothing, p, h[1:(n + 1)], z, extrapolation_left, + extrapolation_right, cache_parameters, linear_lookup) I = cumulative_integral(A, cache_parameters) - CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup) + CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolation_left, + extrapolation_right, cache_parameters, linear_lookup) end function CubicSpline( - u::uType, t; extrapolate = false, cache_parameters = false, + u::uType, t; extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false, assume_linear_t = 1e-2) where {uType <: AbstractVector} + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) - 1 h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0) @@ -522,13 +649,16 @@ function CubicSpline( p = CubicSplineParameterCache(u, h, z, cache_parameters) A = CubicSpline( - u, t, nothing, p, h[1:(n + 1)], z, extrapolate, cache_parameters, assume_linear_t) + u, t, nothing, p, h[1:(n + 1)], z, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) I = cumulative_integral(A, cache_parameters) - CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters, assume_linear_t) + CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) end """ - BSplineInterpolation(u, t, d, pVecType, knotVecType; extrapolate = false, safetycopy = true) + BSplineInterpolation(u, t, d, pVecType, knotVecType; extrapolation::ExtrapolationType.T = ExtrapolationType.none, extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none) 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. @@ -543,8 +673,14 @@ Extrapolation is a constant polynomial of the end points on each side. ## Keyword Arguments - - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. - - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for + - `extrapolation`: The extrapolation type applied left and right of the data. Possible options + are `ExtrapolationType.none` (default), `ExtrapolationType.constant`, + `ExtrapolationType.linear` and `ExtrapolationType.extension`. + - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. + - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. + - `assume_linear_t`: boolean value to specify a faster index lookup behavior for evenly-distributed abscissae. Alternatively, a numerical threshold may be specified for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. @@ -560,7 +696,8 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, scType, T, N} <: sc::scType # Spline coefficients (preallocated memory) pVecType::Symbol knotVecType::Symbol - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} linear_lookup::Bool function BSplineInterpolation(u, @@ -572,7 +709,8 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, scType, T, N} <: sc, pVecType, knotVecType, - extrapolate, + extrapolation_left, + extrapolation_right, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) N = get_output_dim(u) @@ -586,7 +724,8 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, scType, T, N} <: sc, pVecType, knotVecType, - extrapolate, + extrapolation_left, + extrapolation_right, Guesser(t), linear_lookup ) @@ -594,7 +733,12 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, scType, T, N} <: end function BSplineInterpolation( - u::AbstractVector, t, d, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2) + u::AbstractVector, t, d, pVecType, knotVecType; + extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, assume_linear_t = 1e-2) + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d+1) points.") @@ -659,12 +803,18 @@ function BSplineInterpolation( c = vec(sc \ u[:, :]) sc = zeros(eltype(t), n) BSplineInterpolation( - u, t, d, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t) + u, t, d, p, k, c, sc, pVecType, knotVecType, + extrapolation_left, extrapolation_right, assume_linear_t) end function BSplineInterpolation( - u::AbstractArray{T, N}, t, d, pVecType, knotVecType; extrapolate = false, + u::AbstractArray{T, N}, t, d, pVecType, knotVecType; + extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, assume_linear_t = 1e-2) where {T, N} + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d+1) points.") @@ -732,11 +882,13 @@ function BSplineInterpolation( c = reshape(c, size(u)...) sc = zeros(eltype(t), n) BSplineInterpolation( - u, t, d, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t) + u, t, d, p, k, c, sc, pVecType, knotVecType, + extrapolation_left, extrapolation_right, assume_linear_t) end """ - BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolate = false) + BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolation::ExtrapolationType.T = ExtrapolationType.none, extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none) 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](http://www.cad.zju.edu.cn/home/zhx/GM/009/00-bsia.pdf). @@ -753,7 +905,13 @@ Extrapolation is a constant polynomial of the end points on each side. ## Keyword Arguments - - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `extrapolation`: The extrapolation type applied left and right of the data. Possible options + are `ExtrapolationType.none` (default), `ExtrapolationType.constant`, + `ExtrapolationType.linear` and `ExtrapolationType.extension`. + - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. + - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for evenly-distributed abscissae. Alternatively, a numerical threshold may be specified for a test based on the normalized standard deviation of the difference with respect @@ -771,7 +929,8 @@ struct BSplineApprox{uType, tType, pType, kType, cType, scType, T, N} <: sc::scType # Spline coefficients (preallocated memory) pVecType::Symbol knotVecType::Symbol - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} linear_lookup::Bool function BSplineApprox(u, @@ -784,7 +943,8 @@ struct BSplineApprox{uType, tType, pType, kType, cType, scType, T, N} <: sc, pVecType, knotVecType, - extrapolate, + extrapolation_left, + extrapolation_right, assume_linear_t ) linear_lookup = seems_linear(assume_linear_t, t) @@ -800,7 +960,8 @@ struct BSplineApprox{uType, tType, pType, kType, cType, scType, T, N} <: sc, pVecType, knotVecType, - extrapolate, + extrapolation_left, + extrapolation_right, Guesser(t), linear_lookup ) @@ -808,7 +969,12 @@ struct BSplineApprox{uType, tType, pType, kType, cType, scType, T, N} <: end function BSplineApprox( - u::AbstractVector, t, d, h, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2) + u::AbstractVector, t, d, h, pVecType, knotVecType; + extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, assume_linear_t = 1e-2) + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d+1) control points.") @@ -894,12 +1060,18 @@ function BSplineApprox( c[2:(end - 1)] .= vec(P) sc = zeros(eltype(t), h) BSplineApprox( - u, t, d, h, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t) + u, t, d, h, p, k, c, sc, pVecType, knotVecType, + extrapolation_left, extrapolation_right, assume_linear_t) end function BSplineApprox( - u::AbstractArray{T, N}, t, d, h, pVecType, knotVecType; extrapolate = false, + u::AbstractArray{T, N}, t, d, h, pVecType, knotVecType; + extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, assume_linear_t = 1e-2) where {T, N} + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d+1) control points.") @@ -990,10 +1162,12 @@ function BSplineApprox( c[ax_u..., 2:(end - 1)] = P sc = zeros(eltype(t), h) BSplineApprox( - u, t, d, h, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t) + u, t, d, h, p, k, c, sc, pVecType, knotVecType, + extrapolation_left, extrapolation_right, assume_linear_t) end """ - CubicHermiteSpline(du, u, t; extrapolate = false, cache_parameters = false) + CubicHermiteSpline(du, u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.none, extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = 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. @@ -1005,7 +1179,13 @@ It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomi ## Keyword Arguments - - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `extrapolation`: The extrapolation type applied left and right of the data. Possible options + are `ExtrapolationType.none` (default), `ExtrapolationType.constant`, + `ExtrapolationType.linear` and `ExtrapolationType.extension`. + - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. + - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`. - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for evenly-distributed abscissae. Alternatively, a numerical threshold may be specified @@ -1019,33 +1199,43 @@ struct CubicHermiteSpline{uType, tType, IType, duType, pType, T, N} <: t::tType I::IType p::CubicHermiteParameterCache{pType} - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function CubicHermiteSpline( - du, u, t, I, p, extrapolate, cache_parameters, assume_linear_t) + du, u, t, I, p, extrapolation_left, extrapolation_right, + cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) N = get_output_dim(u) new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u), N}( - du, u, t, I, p, extrapolate, Guesser(t), cache_parameters, linear_lookup) + du, u, t, I, p, extrapolation_left, extrapolation_right, + Guesser(t), cache_parameters, linear_lookup) end end function CubicHermiteSpline( - du, u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2) + du, u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false, assume_linear_t = 1e-2) @assert length(u)==length(du) "Length of `u` is not equal to length of `du`." + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) linear_lookup = seems_linear(assume_linear_t, t) p = CubicHermiteParameterCache(du, u, t, cache_parameters) A = CubicHermiteSpline( - du, u, t, nothing, p, extrapolate, cache_parameters, linear_lookup) + du, u, t, nothing, p, extrapolation_left, + extrapolation_right, cache_parameters, linear_lookup) I = cumulative_integral(A, cache_parameters) - CubicHermiteSpline(du, u, t, I, p, extrapolate, cache_parameters, linear_lookup) + CubicHermiteSpline(du, u, t, I, p, extrapolation_left, + extrapolation_right, cache_parameters, linear_lookup) end """ - PCHIPInterpolation(u, t; extrapolate = false, safetycopy = true) + PCHIPInterpolation(u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.none, extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none) It is a PCHIP Interpolation, which is a type of [`CubicHermiteSpline`](@ref) where the derivative values `du` are derived from the input data in such a way that the interpolation never overshoots the data. See [here](https://www.mathworks.com/content/dam/mathworks/mathworks-dot-com/moler/interp.pdf), @@ -1058,22 +1248,28 @@ section 3.4 for more details. ## Keyword Arguments - - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `extrapolation`: The extrapolation type applied left and right of the data. Possible options + are `ExtrapolationType.none` (default), `ExtrapolationType.constant`, + `ExtrapolationType.linear` and `ExtrapolationType.extension`. + - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. + - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`. - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for evenly-distributed abscissae. Alternatively, a numerical threshold may be specified for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. """ -function PCHIPInterpolation( - u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2) +function PCHIPInterpolation(u, t; kwargs...) u, t = munge_data(u, t) du = du_PCHIP(u, t) - CubicHermiteSpline(du, u, t; extrapolate, cache_parameters, assume_linear_t) + CubicHermiteSpline(du, u, t; kwargs...) end """ - QuinticHermiteSpline(ddu, du, u, t; extrapolate = false, safetycopy = true) + QuinticHermiteSpline(ddu, du, u, t; extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none) 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. @@ -1086,7 +1282,13 @@ It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polyno ## Keyword Arguments - - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `extrapolation`: The extrapolation type applied left and right of the data. Possible options + are `ExtrapolationType.none` (default), `ExtrapolationType.constant`, + `ExtrapolationType.linear` and `ExtrapolationType.extension`. + - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. + - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for + the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`. - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`. - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for evenly-distributed abscissae. Alternatively, a numerical threshold may be specified @@ -1101,29 +1303,39 @@ struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T, N} < t::tType I::IType p::QuinticHermiteParameterCache{pType} - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function QuinticHermiteSpline( - ddu, du, u, t, I, p, extrapolate, cache_parameters, assume_linear_t) + ddu, du, u, t, I, p, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) N = get_output_dim(u) new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(ddu), typeof(p.c₁), eltype(u), N}( - ddu, du, u, t, I, p, extrapolate, Guesser(t), cache_parameters, linear_lookup) + ddu, du, u, t, I, p, extrapolation_left, extrapolation_right, + Guesser(t), cache_parameters, linear_lookup) end end -function QuinticHermiteSpline(ddu, du, u, t; extrapolate = false, +function QuinticHermiteSpline( + ddu, du, u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_left::ExtrapolationType.T = ExtrapolationType.none, + extrapolation_right::ExtrapolationType.T = ExtrapolationType.none, cache_parameters = false, assume_linear_t = 1e-2) @assert length(u)==length(du)==length(ddu) "Length of `u` is not equal to length of `du` or `ddu`." + extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) linear_lookup = seems_linear(assume_linear_t, t) p = QuinticHermiteParameterCache(ddu, du, u, t, cache_parameters) A = QuinticHermiteSpline( - ddu, du, u, t, nothing, p, extrapolate, cache_parameters, linear_lookup) + ddu, du, u, t, nothing, p, extrapolation_left, + extrapolation_right, cache_parameters, linear_lookup) I = cumulative_integral(A, cache_parameters) QuinticHermiteSpline( - ddu, du, u, t, I, p, extrapolate, cache_parameters, linear_lookup) + ddu, du, u, t, I, p, extrapolation_left, + extrapolation_right, cache_parameters, linear_lookup) end diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 54129ae6..266c80c0 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,7 +1,43 @@ function _interpolate(A, t) - ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && - throw(ExtrapolationError()) - return _interpolate(A, t, A.iguesser) + if t < first(A.t) + _extrapolate_left(A, t) + elseif t > last(A.t) + _extrapolate_right(A, t) + else + _interpolate(A, t, A.iguesser) + end +end + +function _extrapolate_left(A, t) + (; extrapolation_left) = A + if extrapolation_left == ExtrapolationType.none + throw(LeftExtrapolationError()) + elseif extrapolation_left == ExtrapolationType.constant + slope = derivative(A, first(A.t)) + first(A.u) + slope * zero(t) + elseif extrapolation_left == ExtrapolationType.linear + slope = derivative(A, first(A.t)) + first(A.u) + slope * (t - first(A.t)) + else + # extrapolation_left == ExtrapolationType.extension + _interpolate(A, t, A.iguesser) + end +end + +function _extrapolate_right(A, t) + (; extrapolation_right) = A + if extrapolation_right == ExtrapolationType.none + throw(RightExtrapolationError()) + elseif extrapolation_right == ExtrapolationType.constant + slope = derivative(A, last(A.t)) + last(A.u) + slope * zero(t) + elseif extrapolation_right == ExtrapolationType.linear + slope = derivative(A, last(A.t)) + last(A.u) + slope * (t - last(A.t)) + else + # extrapolation_right == ExtrapolationType.extension + _interpolate(A, t, A.iguesser) + end end # Linear Interpolation diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 2356942c..f12bfa07 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -296,7 +296,18 @@ function integrate_quintic_polynomial(t1, t2, offset, a, b, c, d, e, f) t2_rel = t2 - offset t_sum = t1_rel + t2_rel t_sq_sum = t1_rel^2 + t2_rel^2 + t_cb_sum = t1_rel^3 + t2_rel^3 Δt = t2 - t1 - Δt * (a + t_sum * (b / 2 + d * t_sq_sum / 4) + c * (t_sq_sum + t1_rel * t2_rel) / 3) + - e * (t2_rel^5 - t1_rel^5) / 5 + f * (t2_rel^6 - t1_rel^6) / 6 + cube_diff_factor = t_sq_sum + t1_rel * t2_rel + Δt * (a + t_sum * (b / 2 + d * t_sq_sum / 4) + + cube_diff_factor * (c / 3 + f * t_cb_sum / 6)) + + e * (t2_rel^5 - t1_rel^5) / 5 +end + +function munge_extrapolation(extrapolation, extrapolation_left, extrapolation_right) + if extrapolation == ExtrapolationType.none + extrapolation_left, extrapolation_right + else + extrapolation, extrapolation + end end diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 37595a4f..e9e33328 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -9,7 +9,11 @@ using Optim using ForwardDiff function test_derivatives(method; args = [], kwargs = [], name::String) - func = method(args...; kwargs..., extrapolate = true) + kwargs_extrapolation = (method == Curvefit) ? + [:extrapolate => true] : + [:extrapolation_right => ExtrapolationType.extension, + :extrapolation_left => ExtrapolationType.extension] + func = method(args...; kwargs..., kwargs_extrapolation...) (; t) = func trange = collect(range(minimum(t) - 5.0, maximum(t) + 5.0, step = 0.1)) trange_exclude = filter(x -> !in(x, t), trange) @@ -68,8 +72,14 @@ function test_derivatives(method; args = [], kwargs = [], name::String) @test_throws DataInterpolations.DerivativeNotFoundError derivative( func, t[1], 3) func = method(args...) - @test_throws DataInterpolations.ExtrapolationError derivative(func, t[1] - 1.0) - @test_throws DataInterpolations.ExtrapolationError derivative(func, t[end] + 1.0) + if method == Curvefit + @test_throws DataInterpolations.ExtrapolationError derivative(func, t[1] - 1.0) + @test_throws DataInterpolations.ExtrapolationError derivative(func, t[end] + 1.0) + else + @test_throws DataInterpolations.LeftExtrapolationError derivative(func, t[1] - 1.0) + @test_throws DataInterpolations.RightExtrapolationError derivative( + func, t[end] + 1.0) + end @test_throws DataInterpolations.DerivativeNotFoundError derivative( func, t[1], 3) end @@ -232,7 +242,7 @@ end t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3] test_derivatives(CubicHermiteSpline; args = [du, u, t], name = "Cubic Hermite Spline") - A = CubicHermiteSpline(du, u, t; extrapolate = true) + A = CubicHermiteSpline(du, u, t; extrapolation = ExtrapolationType.extension) @test derivative.(Ref(A), t) ≈ du @test derivative(A, 100.0)≈0.0105409 rtol=1e-5 @test derivative(A, 300.0)≈-0.0806717 rtol=1e-5 @@ -245,7 +255,7 @@ end t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3] test_derivatives(QuinticHermiteSpline; args = [ddu, du, u, t], name = "Quintic Hermite Spline") - A = QuinticHermiteSpline(ddu, du, u, t; extrapolate = true) + A = QuinticHermiteSpline(ddu, du, u, t; extrapolation = ExtrapolationType.extension) @test derivative.(Ref(A), t) ≈ du @test derivative.(Ref(A), t, 2) ≈ ddu @test derivative(A, 100.0)≈0.0103916 rtol=1e-5 @@ -324,7 +334,7 @@ end @testset "Jacobian tests" begin u = rand(5) t = 0:4 - interp = LinearInterpolation(u, t, extrapolate = true) + interp = LinearInterpolation(u, t, extrapolation = ExtrapolationType.extension) grad1 = ForwardDiff.derivative(interp, 2.4) myvec = rand(20) .* 4.0 diff --git a/test/extrapolation_tests.jl b/test/extrapolation_tests.jl new file mode 100644 index 00000000..0429c540 --- /dev/null +++ b/test/extrapolation_tests.jl @@ -0,0 +1,103 @@ +using DataInterpolations + +function test_extrapolation_errors(method, u, t) + A = method(u, t) + @test A.extrapolation_right == ExtrapolationType.none + @test A.extrapolation_left == ExtrapolationType.none + for (error_type, t_eval) in zip( + (DataInterpolations.LeftExtrapolationError, + DataInterpolations.RightExtrapolationError), + (first(t) - 1, last(t) + 1)) + @test_throws error_type A(t_eval) + @test_throws error_type DataInterpolations.derivative( + A, t_eval) + @test_throws error_type DataInterpolations.derivative( + A, t_eval, 2) + @test_throws error_type DataInterpolations.integral( + A, t_eval) + end +end + +function test_constant_extrapolation(method, u, t) + A = method(u, t; extrapolation_left = ExtrapolationType.constant, + extrapolation_right = ExtrapolationType.constant) + t_lower = first(t) - 1 + t_upper = last(t) + 1 + @test A(t_lower) == first(u) + @test A(t_upper) == last(u) + @test DataInterpolations.derivative(A, t_lower) == 0 + @test DataInterpolations.derivative(A, t_upper) == 0 + @test DataInterpolations.integral(A, t_lower, first(t)) ≈ + first(u) * (first(t) - t_lower) + @test DataInterpolations.integral(A, last(t), t_upper) ≈ last(u) * (t_upper - last(t)) +end + +@testset "Linear Interpolation" begin + u = [1.0, 2.0] + t = [1.0, 2.0] + + test_extrapolation_errors(LinearInterpolation, u, t) + test_constant_extrapolation(LinearInterpolation, u, t) + + for extrapolation_type in [ExtrapolationType.linear, ExtrapolationType.extension] + # Left extrapolation + A = LinearInterpolation(u, t; extrapolation_left = extrapolation_type) + t_eval = 0.0 + @test A(t_eval) == 0.0 + @test DataInterpolations.derivative(A, t_eval) == 1.0 + @test DataInterpolations.derivative(A, t_eval, 2) == 0.0 + @test DataInterpolations.integral(A, t_eval) == -0.5 + t_eval = 3.0 + + # Right extrapolation + A = LinearInterpolation(u, t; extrapolation_right = extrapolation_type) + t_eval = 3.0 + @test A(t_eval) == 3.0 + @test DataInterpolations.derivative(A, t_eval) == 1.0 + @test DataInterpolations.derivative(A, t_eval, 2) == 0.0 + @test DataInterpolations.integral(A, t_eval) == 4.0 + t_eval = 0.0 + end +end + +@testset "Quadratic Interpolation" begin + u = [1.0, 3.0, 2.0] + t = 1:3 + + test_extrapolation_errors(QuadraticInterpolation, u, t) + test_constant_extrapolation(LinearInterpolation, u, t) + + # Linear left extrapolation + A = QuadraticInterpolation(u, t; extrapolation_left = ExtrapolationType.linear) + t_eval = 0.0 + @test A(t_eval) ≈ -2.5 + @test DataInterpolations.derivative(A, t_eval) == 3.5 + @test DataInterpolations.derivative(A, t_eval, 2) == 0.0 + @test DataInterpolations.integral(A, t_eval) ≈ 0.75 + + # Linear right extrapolation + A = QuadraticInterpolation(u, t; extrapolation_right = ExtrapolationType.linear) + t_eval = 4.0 + @test A(t_eval) ≈ -0.5 + @test DataInterpolations.derivative(A, t_eval) == -2.5 + @test DataInterpolations.derivative(A, t_eval, 2) == 0.0 + @test DataInterpolations.integral(A, t[end], t_eval) ≈ 0.75 + + # Extension left extrapolation + f = t -> (-3t^2 + 13t - 8) / 2 + df = t -> (-6t + 13) / 2 + A = QuadraticInterpolation(u, t; extrapolation_left = ExtrapolationType.extension) + t_eval = 0.0 + @test A(t_eval) ≈ -4.0 + @test DataInterpolations.derivative(A, t_eval) == df(t_eval) + @test DataInterpolations.derivative(A, t_eval, 2) == -3 + @test DataInterpolations.integral(A, t_eval) ≈ 1.25 + + # Extension right extrapolation + A = QuadraticInterpolation(u, t; extrapolation_right = ExtrapolationType.extension) + t_eval = 4.0 + @test A(t_eval) ≈ -2.0 + @test DataInterpolations.derivative(A, t_eval) == df(t_eval) + @test DataInterpolations.derivative(A, t_eval, 2) == -3 + @test DataInterpolations.integral(A, t_eval) ≈ 5.25 +end diff --git a/test/integral_inverse_tests.jl b/test/integral_inverse_tests.jl index fd83be67..9ba0b158 100644 --- a/test/integral_inverse_tests.jl +++ b/test/integral_inverse_tests.jl @@ -3,7 +3,7 @@ using DataInterpolations: integral, derivative, invert_integral using FiniteDifferences function test_integral_inverses(method; args = [], kwargs = []) - A = method(args...; kwargs..., extrapolate = true) + A = method(args...; kwargs..., extrapolation = ExtrapolationType.extension) @test hasfield(typeof(A), :I) A_intinv = invert_integral(A) @test A_intinv isa DataInterpolations.AbstractIntegralInverseInterpolation diff --git a/test/integral_tests.jl b/test/integral_tests.jl index 3d59a6c0..09ba655b 100644 --- a/test/integral_tests.jl +++ b/test/integral_tests.jl @@ -6,7 +6,8 @@ using RegularizationTools using StableRNGs function test_integral(method; args = [], kwargs = [], name::String) - func = method(args...; kwargs..., extrapolate = true) + func = method(args...; kwargs..., extrapolation_left = ExtrapolationType.extension, + extrapolation_right = ExtrapolationType.extension) (; t) = func t1 = minimum(t) t2 = maximum(t) @@ -48,10 +49,11 @@ function test_integral(method; args = [], kwargs = [], name::String) @test isapprox(qint, aint, atol = 1e-6, rtol = 1e-8) end func = method(args...; kwargs...) - @test_throws DataInterpolations.ExtrapolationError integral(func, t[1] - 1.0) - @test_throws DataInterpolations.ExtrapolationError integral(func, t[end] + 1.0) - @test_throws DataInterpolations.ExtrapolationError integral(func, t[1] - 1.0, t[2]) - @test_throws DataInterpolations.ExtrapolationError integral(func, t[1], t[end] + 1.0) + @test_throws DataInterpolations.LeftExtrapolationError integral(func, t[1] - 1.0) + @test_throws DataInterpolations.RightExtrapolationError integral(func, t[end] + 1.0) + @test_throws DataInterpolations.LeftExtrapolationError integral(func, t[1] - 1.0, t[2]) + @test_throws DataInterpolations.RightExtrapolationError integral( + func, t[1], t[end] + 1.0) end @testset "LinearInterpolation" begin diff --git a/test/interface.jl b/test/interface.jl index e7b2b81b..29f6864d 100644 --- a/test/interface.jl +++ b/test/interface.jl @@ -18,8 +18,8 @@ end @testset "Symbolics" begin u = 2.0collect(1:10) t = 1.0collect(1:10) - A = LinearInterpolation(u, t; extrapolate = true) - B = LinearInterpolation(u .^ 2, t; extrapolate = true) + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) + B = LinearInterpolation(u .^ 2, t; extrapolation = ExtrapolationType.extension) @variables t x(t) substitute(A(t), Dict(t => x)) t_val = 2.7 diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 65b48113..83a0bef1 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -8,7 +8,8 @@ function test_interpolation_type(T) @test T <: DataInterpolations.AbstractInterpolation @test hasfield(T, :u) @test hasfield(T, :t) - @test hasfield(T, :extrapolate) + @test hasfield(T, :extrapolation_right) + @test hasfield(T, :extrapolation_left) @test hasfield(T, :iguesser) @test !isempty(methods(DataInterpolations._interpolate, (T, Any, Number))) @test !isempty(methods(DataInterpolations._integral, (T, Number, Number, Number))) @@ -30,7 +31,7 @@ end for t in (1.0:10.0, 1.0collect(1:10)) u = 2.0collect(1:10) #t = 1.0collect(1:10) - A = LinearInterpolation(u, t; extrapolate = true) + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) for (_t, _u) in zip(t, u) @test A(_t) == _u @@ -40,7 +41,7 @@ end @test A(11) == 22 u = vcat(2.0collect(1:10)', 3.0collect(1:10)') - A = LinearInterpolation(u, t; extrapolate = true) + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) for (_t, _u) in zip(t, eachcol(u)) @test A(_t) == _u @@ -53,7 +54,7 @@ end y = 2:4 u_ = x' .* y u = [u_[:, i] for i in 1:size(u_, 2)] - A = LinearInterpolation(u, t; extrapolate = true) + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(0) == [0.0, 0.0, 0.0] @test A(5.5) == [11.0, 16.5, 22.0] @test A(11) == [22.0, 33.0, 44.0] @@ -64,7 +65,7 @@ end u_ = x' .* y u = [u_[:, i:(i + 1)] for i in 1:2:10] t = 1.0collect(2:2:10) - A = LinearInterpolation(u, t; extrapolate = true) + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(0) == [-2.0 0.0; -3.0 0.0; -4.0 0.0] @test A(3) == [4.0 6.0; 6.0 9.0; 8.0 12.0] @@ -74,7 +75,7 @@ end # with NaNs (#113) u = [NaN, 1.0, 2.0, 3.0] t = 1:4 - A = LinearInterpolation(u, t; extrapolate = true) + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test isnan(A(1.0)) @test A(2.0) == 1.0 @test A(2.5) == 1.5 @@ -82,7 +83,7 @@ end @test A(4.0) == 3.0 u = [0.0, NaN, 2.0, 3.0] - A = LinearInterpolation(u, t; extrapolate = true) + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(1.0) == 0.0 @test isnan(A(2.0)) @test isnan(A(2.5)) @@ -90,7 +91,7 @@ end @test A(4.0) == 3.0 u = [0.0, 1.0, NaN, 3.0] - A = LinearInterpolation(u, t; extrapolate = true) + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(1.0) == 0.0 @test A(2.0) == 1.0 @test isnan(A(2.5)) @@ -98,7 +99,7 @@ end @test A(4.0) == 3.0 u = [0.0, 1.0, 2.0, NaN] - A = LinearInterpolation(u, t; extrapolate = true) + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(1.0) == 0.0 @test A(2.0) == 1.0 @test A(3.0) == 2.0 @@ -108,16 +109,16 @@ end # Test type stability u = Float32.(1:5) t = Float32.(1:5) - A1 = LinearInterpolation(u, t; extrapolate = true) + A1 = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) u = 1:5 t = 1:5 - A2 = LinearInterpolation(u, t; extrapolate = true) + A2 = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) u = [1 // i for i in 1:5] t = (1:5) - A3 = LinearInterpolation(u, t; extrapolate = true) + A3 = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) u = [1 // i for i in 1:5] t = [1 // (6 - i) for i in 1:5] - A4 = LinearInterpolation(u, t; extrapolate = true) + A4 = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) F32 = Float32(1) F64 = Float64(1) @@ -137,13 +138,13 @@ end # Nan time value: t = 0.0:3 # Floats u = [0, -2, -1, -2] - A = LinearInterpolation(u, t; extrapolate = true) + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) dA = t -> ForwardDiff.derivative(A, t) @test isnan(dA(NaN)) t = 0:3 # Integers u = [0, -2, -1, -2] - A = LinearInterpolation(u, t; extrapolate = true) + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) dA = t -> ForwardDiff.derivative(A, t) @test isnan(dA(NaN)) @@ -156,7 +157,7 @@ end # Test array-valued interpolation u = collect.(2.0collect(1:10)) t = 1.0collect(1:10) - A = LinearInterpolation(u, t; extrapolate = true) + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(0) == fill(0.0) @test A(5.5) == fill(11.0) @test A(11) == fill(22) @@ -171,13 +172,13 @@ end # Test extrapolation u = 2.0collect(1:10) t = 1.0collect(1:10) - A = LinearInterpolation(u, t; extrapolate = true) + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(-1.0) == -2.0 @test A(11.0) == 22.0 A = LinearInterpolation(u, t) - @test_throws DataInterpolations.ExtrapolationError A(-1.0) - @test_throws DataInterpolations.ExtrapolationError A(11.0) - @test_throws DataInterpolations.ExtrapolationError A([-1.0, 11.0]) + @test_throws DataInterpolations.LeftExtrapolationError A(-1.0) + @test_throws DataInterpolations.RightExtrapolationError A(11.0) + @test_throws DataInterpolations.LeftExtrapolationError A([-1.0, 11.0]) end @testset "Quadratic Interpolation" begin @@ -185,7 +186,7 @@ end u = [1.0, 4.0, 9.0, 16.0] t = [1.0, 2.0, 3.0, 4.0] - A = QuadraticInterpolation(u, t; extrapolate = true) + A = QuadraticInterpolation(u, t; extrapolation = ExtrapolationType.extension) for (_t, _u) in zip(t, u) @test A(_t) == _u @@ -200,7 +201,8 @@ end # backward-looking interpolation u = [1.0, 4.0, 9.0, 16.0] t = [1.0, 2.0, 3.0, 4.0] - A = QuadraticInterpolation(u, t, :Backward; extrapolate = true) + A = QuadraticInterpolation( + u, t, :Backward; extrapolation = ExtrapolationType.extension) for (_t, _u) in zip(t, u) @test A(_t) == _u @@ -238,7 +240,7 @@ end # Matrix interpolation test u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0] - A = QuadraticInterpolation(u, t; extrapolate = true) + A = QuadraticInterpolation(u, t; extrapolation = ExtrapolationType.extension) for (_t, _u) in zip(t, eachcol(u)) @test A(_t) == _u @@ -251,7 +253,7 @@ end u_ = [1.0, 4.0, 9.0, 16.0]' .* ones(5) u = [u_[:, i] for i in 1:size(u_, 2)] - A = QuadraticInterpolation(u, t; extrapolate = true) + A = QuadraticInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(0) == zeros(5) @test A(1.5) == 2.25 * ones(5) @test A(2.5) == 6.25 * ones(5) @@ -259,7 +261,7 @@ end @test A(5.0) == 25.0 * ones(5) u = [repeat(u[i], 1, 3) for i in 1:4] - A = QuadraticInterpolation(u, t; extrapolate = true) + A = QuadraticInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(0) == zeros(5, 3) @test A(1.5) == 2.25 * ones(5, 3) @test A(2.5) == 6.25 * ones(5, 3) @@ -269,12 +271,12 @@ end # Test extrapolation u = [1.0, 4.5, 6.0, 2.0] t = [1.0, 2.0, 3.0, 4.0] - A = QuadraticInterpolation(u, t; extrapolate = true) + A = QuadraticInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(0.0) == -4.5 @test A(5.0) == -7.5 A = QuadraticInterpolation(u, t) - @test_throws DataInterpolations.ExtrapolationError A(0.0) - @test_throws DataInterpolations.ExtrapolationError A(5.0) + @test_throws DataInterpolations.LeftExtrapolationError A(0.0) + @test_throws DataInterpolations.RightExtrapolationError A(5.0) end @testset "Lagrange Interpolation" begin @@ -329,12 +331,12 @@ end # Test extrapolation u = [1.0, 4.0, 9.0] t = [1.0, 2.0, 3.0] - A = LagrangeInterpolation(u, t; extrapolate = true) + A = LagrangeInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(0.0) == 0.0 @test A(4.0) == 16.0 A = LagrangeInterpolation(u, t) - @test_throws DataInterpolations.ExtrapolationError A(-1.0) - @test_throws DataInterpolations.ExtrapolationError A(4.0) + @test_throws DataInterpolations.LeftExtrapolationError A(-1.0) + @test_throws DataInterpolations.RightExtrapolationError A(4.0) end @testset "Akima Interpolation" begin @@ -360,12 +362,12 @@ end test_cached_index(A) # Test extrapolation - A = AkimaInterpolation(u, t; extrapolate = true) + A = AkimaInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(-1.0) ≈ -5.0 @test A(11.0) ≈ -3.924742268041234 A = AkimaInterpolation(u, t) - @test_throws DataInterpolations.ExtrapolationError A(-1.0) - @test_throws DataInterpolations.ExtrapolationError A(11.0) + @test_throws DataInterpolations.LeftExtrapolationError A(-1.0) + @test_throws DataInterpolations.RightExtrapolationError A(11.0) end @testset "ConstantInterpolation" begin @@ -374,7 +376,8 @@ end t = [1.0, 2.0, 3.0, 4.0] @testset "Vector case" for u in [[1.0, 2.0, 0.0, 1.0], ["B", "C", "A", "B"]] - A = ConstantInterpolation(u, t, dir = :right; extrapolate = true) + A = ConstantInterpolation( + u, t, dir = :right; extrapolation = ExtrapolationType.extension) @test A(0.5) == u[1] @test A(1.0) == u[1] @test A(1.5) == u[2] @@ -386,7 +389,7 @@ end @test A(4.5) == u[1] test_cached_index(A) - A = ConstantInterpolation(u, t; extrapolate = true) # dir=:left is default + A = ConstantInterpolation(u, t; extrapolation = ExtrapolationType.extension) # dir=:left is default @test A(0.5) == u[1] @test A(1.0) == u[1] @test A(1.5) == u[1] @@ -403,7 +406,8 @@ end [1.0 2.0 0.0 1.0; 1.0 2.0 0.0 1.0], ["B" "C" "A" "B"; "B" "C" "A" "B"] ] - A = ConstantInterpolation(u, t, dir = :right; extrapolate = true) + A = ConstantInterpolation( + u, t, dir = :right; extrapolation = ExtrapolationType.extension) @test A(0.5) == u[:, 1] @test A(1.0) == u[:, 1] @test A(1.5) == u[:, 2] @@ -415,7 +419,7 @@ end @test A(4.5) == u[:, 1] test_cached_index(A) - A = ConstantInterpolation(u, t; extrapolate = true) # dir=:left is default + A = ConstantInterpolation(u, t; extrapolation = ExtrapolationType.extension) # dir=:left is default @test A(0.5) == u[:, 1] @test A(1.0) == u[:, 1] @test A(1.5) == u[:, 1] @@ -431,7 +435,8 @@ end @testset "Vector of Vectors case" for u in [ [[1.0, 2.0], [0.0, 1.0], [1.0, 2.0], [0.0, 1.0]], [["B", "C"], ["A", "B"], ["B", "C"], ["A", "B"]]] - A = ConstantInterpolation(u, t, dir = :right; extrapolate = true) + A = ConstantInterpolation( + u, t, dir = :right; extrapolation = ExtrapolationType.extension) @test A(0.5) == u[1] @test A(1.0) == u[1] @test A(1.5) == u[2] @@ -443,7 +448,7 @@ end @test A(4.5) == u[4] test_cached_index(A) - A = ConstantInterpolation(u, t; extrapolate = true) # dir=:left is default + A = ConstantInterpolation(u, t; extrapolation = ExtrapolationType.extension) # dir=:left is default @test A(0.5) == u[1] @test A(1.0) == u[1] @test A(1.5) == u[1] @@ -459,7 +464,8 @@ end @testset "Vector of Matrices case" for u in [ [[1.0 2.0; 1.0 2.0], [0.0 1.0; 0.0 1.0], [1.0 2.0; 1.0 2.0], [0.0 1.0; 0.0 1.0]], [["B" "C"; "B" "C"], ["A" "B"; "A" "B"], ["B" "C"; "B" "C"], ["A" "B"; "A" "B"]]] - A = ConstantInterpolation(u, t, dir = :right; extrapolate = true) + A = ConstantInterpolation( + u, t, dir = :right; extrapolation = ExtrapolationType.extension) @test A(0.5) == u[1] @test A(1.0) == u[1] @test A(1.5) == u[2] @@ -471,7 +477,7 @@ end @test A(4.5) == u[4] test_cached_index(A) - A = ConstantInterpolation(u, t; extrapolate = true) # dir=:left is default + A = ConstantInterpolation(u, t; extrapolation = ExtrapolationType.extension) # dir=:left is default @test A(0.5) == u[1] @test A(1.0) == u[1] @test A(1.5) == u[1] @@ -486,17 +492,17 @@ end # Test extrapolation u = [1.0, 2.0, 0.0, 1.0] - A = ConstantInterpolation(u, t; extrapolate = true) + A = ConstantInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(-1.0) == 1.0 @test A(11.0) == 1.0 A = ConstantInterpolation(u, t) - @test_throws DataInterpolations.ExtrapolationError A(-1.0) - @test_throws DataInterpolations.ExtrapolationError A(11.0) + @test_throws DataInterpolations.LeftExtrapolationError A(-1.0) + @test_throws DataInterpolations.RightExtrapolationError A(11.0) # Test extrapolation with infs with regularly spaced t u = [1.67e7, 1.6867e7, 1.7034e7, 1.7201e7, 1.7368e7] t = [0.0, 0.1, 0.2, 0.3, 0.4] - A = ConstantInterpolation(u, t; extrapolate = true) + A = ConstantInterpolation(u, t; extrapolation = ExtrapolationType.extension) @test A(Inf) == last(u) @test A(-Inf) == first(u) end @@ -507,7 +513,7 @@ end u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] - A = QuadraticSpline(u, t; extrapolate = true) + A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.extension) # Solution P₁ = x -> 0.5 * (x + 1) * (x + 2) @@ -523,14 +529,14 @@ end u_ = [0.0, 1.0, 3.0]' .* ones(4) u = [u_[:, i] for i in 1:size(u_, 2)] - A = QuadraticSpline(u, t; extrapolate = true) + A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.extension) @test A(-2.0) == P₁(-2.0) * ones(4) @test A(-0.5) == P₁(-0.5) * ones(4) @test A(0.7) == P₁(0.7) * ones(4) @test A(2.0) == P₁(2.0) * ones(4) u = [repeat(u[i], 1, 3) for i in 1:3] - A = QuadraticSpline(u, t; extrapolate = true) + A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.extension) @test A(-2.0) == P₁(-2.0) * ones(4, 3) @test A(-0.5) == P₁(-0.5) * ones(4, 3) @test A(0.7) == P₁(0.7) * ones(4, 3) @@ -539,12 +545,12 @@ end # Test extrapolation u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] - A = QuadraticSpline(u, t; extrapolate = true) + A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.extension) @test A(-2.0) == 0.0 @test A(2.0) == 6.0 A = QuadraticSpline(u, t) - @test_throws DataInterpolations.ExtrapolationError A(-2.0) - @test_throws DataInterpolations.ExtrapolationError A(2.0) + @test_throws DataInterpolations.LeftExtrapolationError A(-2.0) + @test_throws DataInterpolations.RightExtrapolationError A(2.0) end @testset "CubicSpline Interpolation" begin @@ -553,7 +559,7 @@ end u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] - A = CubicSpline(u, t; extrapolate = true) + A = CubicSpline(u, t; extrapolation = ExtrapolationType.extension) test_cached_index(A) # Solution @@ -572,7 +578,7 @@ end u_ = [0.0, 1.0, 3.0]' .* ones(4) u = [u_[:, i] for i in 1:size(u_, 2)] - A = CubicSpline(u, t; extrapolate = true) + A = CubicSpline(u, t; extrapolation = ExtrapolationType.extension) for x in (-1.5, -0.5, -0.7) @test A(x) ≈ P₁(x) * ones(4) end @@ -581,7 +587,7 @@ end end u = [repeat(u[i], 1, 3) for i in 1:3] - A = CubicSpline(u, t; extrapolate = true) + A = CubicSpline(u, t; extrapolation = ExtrapolationType.extension) for x in (-1.5, -0.5, -0.7) @test A(x) ≈ P₁(x) * ones(4, 3) end @@ -592,12 +598,12 @@ end # Test extrapolation u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] - A = CubicSpline(u, t; extrapolate = true) + A = CubicSpline(u, t; extrapolation = ExtrapolationType.extension) @test A(-2.0) ≈ -1.0 @test A(2.0) ≈ 5.0 A = CubicSpline(u, t) - @test_throws DataInterpolations.ExtrapolationError A(-2.0) - @test_throws DataInterpolations.ExtrapolationError A(2.0) + @test_throws DataInterpolations.LeftExtrapolationError A(-2.0) + @test_throws DataInterpolations.RightExtrapolationError A(2.0) @testset "AbstractMatrix" begin t = 0.1:0.1:1.0 @@ -636,12 +642,13 @@ end test_cached_index(A) # Test extrapolation - A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform; extrapolate = true) + A = BSplineInterpolation( + u, t, 2, :Uniform, :Uniform; extrapolation = ExtrapolationType.extension) @test A(-1.0) == u[1] @test A(300.0) == u[end] A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform) - @test_throws DataInterpolations.ExtrapolationError A(-1.0) - @test_throws DataInterpolations.ExtrapolationError A(300.0) + @test_throws DataInterpolations.LeftExtrapolationError A(-1.0) + @test_throws DataInterpolations.RightExtrapolationError A(300.0) A = BSplineInterpolation(u, t, 2, :ArcLen, :Average) @@ -656,12 +663,13 @@ end @test_nowarn BSplineInterpolation(u[1:3], t[1:3], 2, :Uniform, :Uniform) # Test extrapolation - A = BSplineInterpolation(u, t, 2, :ArcLen, :Average; extrapolate = true) + A = BSplineInterpolation( + u, t, 2, :ArcLen, :Average; extrapolation = ExtrapolationType.extension) @test A(-1.0) == u[1] @test A(300.0) == u[end] A = BSplineInterpolation(u, t, 2, :ArcLen, :Average) - @test_throws DataInterpolations.ExtrapolationError A(-1.0) - @test_throws DataInterpolations.ExtrapolationError A(300.0) + @test_throws DataInterpolations.LeftExtrapolationError A(-1.0) + @test_throws DataInterpolations.RightExtrapolationError A(300.0) @testset "AbstractMatrix" begin t = 0.1:0.1:1.0 @@ -713,12 +721,13 @@ end @test_nowarn BSplineApprox(u, t, 2, 3, :Uniform, :Uniform) # Test extrapolation - A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform; extrapolate = true) + A = BSplineApprox( + u, t, 2, 4, :Uniform, :Uniform; extrapolation = ExtrapolationType.extension) @test A(-1.0) == u[1] @test A(300.0) == u[end] A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform) - @test_throws DataInterpolations.ExtrapolationError A(-1.0) - @test_throws DataInterpolations.ExtrapolationError A(300.0) + @test_throws DataInterpolations.LeftExtrapolationError A(-1.0) + @test_throws DataInterpolations.RightExtrapolationError A(300.0) @testset "AbstractMatrix" begin t = 0.1:0.1:1.0 @@ -759,7 +768,7 @@ end du = [-0.047, -0.058, 0.054, 0.012, -0.068, 0.0] u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22] t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3] - A = CubicHermiteSpline(du, u, t; extrapolate = true) + A = CubicHermiteSpline(du, u, t; extrapolation = ExtrapolationType.extension) @test A.(t) ≈ u @test A(100.0)≈10.106770 rtol=1e-5 @test A(300.0)≈9.901542 rtol=1e-5 @@ -787,7 +796,7 @@ end du = [-0.047, -0.058, 0.054, 0.012, -0.068, 0.0] u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22] t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3] - A = QuinticHermiteSpline(ddu, du, u, t; extrapolate = true) + A = QuinticHermiteSpline(ddu, du, u, t; extrapolation = ExtrapolationType.extension) @test A.(t) ≈ u @test A(100.0)≈10.107996 rtol=1e-5 @test A(300.0)≈11.364162 rtol=1e-5 diff --git a/test/regularization.jl b/test/regularization.jl index fcf7580d..a3f4ef1c 100644 --- a/test/regularization.jl +++ b/test/regularization.jl @@ -180,10 +180,11 @@ end end @testset "Extrapolation" begin - A = RegularizationSmooth(uₒ, tₒ; alg = :fixed, extrapolate = true) + A = RegularizationSmooth( + uₒ, tₒ; alg = :fixed, extrapolation_right = ExtrapolationType.extension) @test A(10.0) == A.Aitp(10.0) A = RegularizationSmooth(uₒ, tₒ; alg = :fixed) - @test_throws DataInterpolations.ExtrapolationError A(10.0) + @test_throws DataInterpolations.RightExtrapolationError A(10.0) end @testset "Type inference" begin diff --git a/test/runtests.jl b/test/runtests.jl index 80080a75..774ad497 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,7 @@ using SafeTestsets @safetestset "Derivative Tests" include("derivative_tests.jl") @safetestset "Integral Tests" include("integral_tests.jl") @safetestset "Integral Inverse Tests" include("integral_inverse_tests.jl") +@safetestset "Extrapolation" include("extrapolation_tests.jl") @safetestset "Online Tests" include("online_tests.jl") @safetestset "Regularization Smoothing" include("regularization.jl") @safetestset "Show methods" include("show.jl") diff --git a/test/zygote_tests.jl b/test/zygote_tests.jl index 7af735af..d8d178df 100644 --- a/test/zygote_tests.jl +++ b/test/zygote_tests.jl @@ -3,7 +3,8 @@ using ForwardDiff using Zygote function test_zygote(method, u, t; args = [], args_after = [], kwargs = [], name::String) - func = method(args..., u, t, args_after...; kwargs..., extrapolate = true) + func = method(args..., u, t, args_after...; kwargs..., + extrapolation = ExtrapolationType.extension) trange = collect(range(minimum(t) - 5.0, maximum(t) + 5.0, step = 0.1)) trange_exclude = filter(x -> !in(x, t), trange) @testset "$name, derivatives w.r.t. input" begin @@ -19,7 +20,8 @@ function test_zygote(method, u, t; args = [], args_after = [], kwargs = [], name [LagrangeInterpolation, BSplineInterpolation, BSplineApprox, QuadraticSpline] @testset "$name, derivatives w.r.t. u" begin function f(u) - A = method(args..., u, t, args_after...; kwargs..., extrapolate = true) + A = method(args..., u, t, args_after...; kwargs..., + extrapolation = ExtrapolationType.extension) out = if u isa AbstractVector{<:Real} zero(eltype(u)) elseif u isa AbstractMatrix