From 80edc15695817f6035853de1ca3e3b423f60e9fd Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 25 Nov 2024 16:20:37 +0100 Subject: [PATCH 01/27] POC --- src/DataInterpolations.jl | 7 ++++--- src/interpolation_caches.jl | 31 +++++++++++++++++++++++++++++++ src/interpolation_methods.jl | 17 ++++++++++++++++- src/interpolation_utils.jl | 15 +++++++++++++++ src/parameter_caches.jl | 25 +++++++++++++++++++++++++ 5 files changed, 91 insertions(+), 4 deletions(-) diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl index 393f1160..1da0f240 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -86,9 +86,10 @@ function Base.showerror(io::IO, e::IntegralNotInvertibleError) end export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, - AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline, - BSplineInterpolation, BSplineApprox, CubicHermiteSpline, PCHIPInterpolation, - QuinticHermiteSpline, LinearInterpolationIntInv, ConstantInterpolationIntInv + AkimaInterpolation, ConstantInterpolation, SmoothedConstantInterpolation, + QuadraticSpline, CubicSpline, BSplineInterpolation, BSplineApprox, + CubicHermiteSpline, PCHIPInterpolation, QuinticHermiteSpline, + LinearInterpolationIntInv, ConstantInterpolationIntInv # added for RegularizationSmooth, JJS 11/27/21 ### Regularization data smoothing and interpolation diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 850ac74b..1d154ee2 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -285,6 +285,37 @@ function ConstantInterpolation( ConstantInterpolation(u, t, I, dir, extrapolate, cache_parameters, assume_linear_t) end +struct SmoothedConstantInterpolation{uType, tType, dmaxType, IType, pType, T, N} <: + AbstractInterpolation{T, N} + u::uType + t::tType + I::IType + p::SmoothedConstantParameterCache{uType, tType} + d_max::dmaxType + extrapolate::Bool + iguesser::Guesser{tType} + cache_parameters::Bool + linear_lookup::Bool + function SmoothedConstantInterpolation( + u, t, I, p, d_max, extrapolate, cache_parameters, assume_linear_t) + linear_lookup = seems_linear(assume_linear_t, t) + N = get_output_dim(u) + new{typeof(u), typeof(t), typeof(d_max), typeof(I), typeof(p.d), eltype(u), N}( + u, t, I, p, d_max, extrapolate, Guesser(t), cache_parameters, linear_lookup) + end +end + +function SmoothedConstantInterpolation(u, t; d_max = Inf, extrapolate = false, + cache_parameters = false, assume_linear_t = 1e-2) + u, t = munge_data(u, t) + p = SmoothedConstantParameterCache(u, t, cache_parameters, d_max) + A = SmoothedConstantInterpolation( + u, t, nothing, p, d_max, extrapolate, cache_parameters, assume_linear_t) + I = cumulative_integral(A, cache_parameters) + SmoothedConstantInterpolation( + u, t, I, p, d_max, extrapolate, cache_parameters, assume_linear_t) +end + """ QuadraticSpline(u, t; extrapolate = false, cache_parameters = false) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 54129ae6..f0f768a9 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -116,7 +116,7 @@ function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess @evalpoly wj A.u[idx] A.b[idx] A.c[idx] A.d[idx] end -# ConstantInterpolation Interpolation +# Constant Interpolation function _interpolate(A::ConstantInterpolation{<:AbstractVector}, t::Number, iguess) if A.dir === :left # :left means that value to the left is used for interpolation @@ -139,6 +139,21 @@ function _interpolate(A::ConstantInterpolation{<:AbstractMatrix}, t::Number, igu A.u[:, idx] end +# Smoothed constant Interpolation +function _interpolate(A::SmoothedConstantInterpolation{<:AbstractVector}, t::Number, iguess) + idx = get_idx(A, t, iguess) + d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx) + out = A.u[idx] + + if (t - A.t[idx]) < d_lower + out -= c_lower * (((t - A.t[idx]) / d_lower - 1))^2 + elseif (A.t[idx + 1] - t) < d_upper + out += c_upper * ((1 - (A.t[idx + 1] - t) / d_upper))^2 + end + + out +end + # QuadraticSpline Interpolation function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) idx = get_idx(A, t, iguess) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 2356942c..a3b719a7 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -207,6 +207,21 @@ function get_parameters(A::LinearInterpolation, idx) end end +function get_parameters(A::SmoothedConstantInterpolation, idx) + if A.cache_parameters + d_lower = A.p.d[idx] + d_upper = A.p.d[idx + 1] + c_lower = A.p.c[idx] + c_upper = A.p.c[idx + 1] + d_lower, d_upper, c_lower, c_upper + else + d_lower, c_lower = smoothed_linear_interpolation_parameters(A.u, A.t, A.d_max, idx) + d_upper, c_upper = smoothed_linear_interpolation_parameters( + A.u, A.t, A.d_max, idx + 1) + d_lower, d_upper, c_lower, c_upper + end +end + function get_parameters(A::QuadraticInterpolation, idx) if A.cache_parameters A.p.α[idx], A.p.β[idx] diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index 737f736b..84785680 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -32,6 +32,31 @@ function linear_interpolation_parameters(u::AbstractArray{T, N}, t, idx) where { return slope end +struct SmoothedConstantParameterCache{uType, tType} + d::tType + c::uType +end + +function SmoothedConstantParameterCache(u, t, cache_parameters, d_max) + if cache_parameters + parameters = smoothed_linear_interpolation_parameters.( + Ref(u), Ref(t), d_max, eachindex(t)) + d, c = collect.(eachrow(stack(collect.(parameters)))) + SmoothedConstantParameterCache(d, c) + else + SmoothedConstantParameterCache(eltype(t)[], eltype(u)[]) + end +end + +function smoothed_linear_interpolation_parameters(u, t, d_max, idx) + # TODO: Add support for making periodic extrapolation smooth + if isone(idx) || (idx == length(t)) + zero(eltype(t)), zero(eltype(u)) + else + min(t[idx] - t[idx - 1], t[idx + 1] - t[idx], 2d_max) / 2, (u[idx] - u[idx - 1]) / 2 + end +end + struct QuadraticParameterCache{pType} α::pType β::pType From ec4c113d9552dca28e4dcf30ed4684f91503b342 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 25 Nov 2024 16:45:53 +0100 Subject: [PATCH 02/27] Add docstring --- src/interpolation_caches.jl | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 1d154ee2..733bdd5d 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -285,6 +285,29 @@ function ConstantInterpolation( ConstantInterpolation(u, t, I, dir, extrapolate, cache_parameters, assume_linear_t) end +""" + SmoothedConstantInterpolation(u, t; d_max = Inf, extrapolate = false, + cache_parameters = false, assume_linear_t = 1e-2) + +It is a method for interpolating constantly with forward fill, with smoothing around the +value transitions to make the curve C1 smooth while the integral never drifts far from +the integral of constant interpolation. + +## Arguments + + - `u`: data points. + - `t`: time points. + +## Keyword Arguments + + - `d_max`: the maximum distance in `t` from the data points the smoothing is allowed to reach. + - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `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. +""" struct SmoothedConstantInterpolation{uType, tType, dmaxType, IType, pType, T, N} <: AbstractInterpolation{T, N} u::uType From 97de11cbd63cf1c7a5f6f247098cf5ce21d29bb2 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 25 Nov 2024 19:24:31 +0100 Subject: [PATCH 03/27] derivatives --- src/derivatives.jl | 13 +++++++++++++ src/interpolation_caches.jl | 6 +++--- src/interpolation_methods.jl | 4 ++-- 3 files changed, 18 insertions(+), 5 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index 7d6dea24..198e0cb6 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -19,6 +19,19 @@ function _derivative(A::LinearInterpolation, t::Number, iguess) slope end +function _derivative(A::SmoothedConstantInterpolation{<:AbstractVector}, t::Number, iguess) + idx = get_idx(A, t, iguess) + d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx) + + if (t - A.t[idx]) < d_lower + -2c_lower * ((t - A.t[idx]) / d_lower - 1) / d_lower + elseif (A.t[idx + 1] - t) < d_upper + 2c_upper * (1 - (A.t[idx + 1] - t) / d_upper) / d_upper + else + zero(c_upper / oneunit(t)) + end +end + function _derivative(A::QuadraticInterpolation, t::Number, iguess) idx = get_idx(A, t, iguess) Δt = t - A.t[idx] diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 733bdd5d..cb2bb0aa 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -290,7 +290,7 @@ end cache_parameters = false, assume_linear_t = 1e-2) It is a method for interpolating constantly with forward fill, with smoothing around the -value transitions to make the curve C1 smooth while the integral never drifts far from +value transitions to make the curve C1 smooth while the integral never drifts far from the integral of constant interpolation. ## Arguments @@ -300,10 +300,10 @@ the integral of constant interpolation. ## Keyword Arguments - - `d_max`: the maximum distance in `t` from the data points the smoothing is allowed to reach. + - `d_max`: the maximum distance in `t` from the data points the smoothing is allowed to reach. - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. - `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. diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index f0f768a9..0615504d 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -146,9 +146,9 @@ function _interpolate(A::SmoothedConstantInterpolation{<:AbstractVector}, t::Num out = A.u[idx] if (t - A.t[idx]) < d_lower - out -= c_lower * (((t - A.t[idx]) / d_lower - 1))^2 + out -= c_lower * ((t - A.t[idx]) / d_lower - 1)^2 elseif (A.t[idx + 1] - t) < d_upper - out += c_upper * ((1 - (A.t[idx + 1] - t) / d_upper))^2 + out += c_upper * (1 - (A.t[idx + 1] - t) / d_upper)^2 end out From aad60819ea5627d724529276b69f58fb046ea77f Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 25 Nov 2024 19:52:45 +0100 Subject: [PATCH 04/27] Add to docs --- README.md | 2 +- docs/src/index.md | 2 +- docs/src/manual.md | 1 + docs/src/methods.md | 18 ++++++++++++++++++ src/interpolation_caches.jl | 4 ++-- 5 files changed, 23 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index df100c0d..0dd6635d 100644 --- a/README.md +++ b/README.md @@ -47,7 +47,7 @@ In all cases, `u` an `AbstractVector` of values and `t` is an `AbstractVector` o corresponding to `(u,t)` pairs. - `ConstantInterpolation(u,t)` - A piecewise constant interpolation. - + - `SmoothedConstantInterpolation(u,t)` - An integral preserving continuously differentiable approximation of constant interpolation. - `LinearInterpolation(u,t)` - A linear interpolation. - `QuadraticInterpolation(u,t)` - A quadratic interpolation. - `LagrangeInterpolation(u,t,n)` - A Lagrange interpolation of order `n`. diff --git a/docs/src/index.md b/docs/src/index.md index 1b19d50f..e74bbd8b 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -17,7 +17,7 @@ In all cases, `u` an `AbstractVector` of values and `t` is an `AbstractVector` o corresponding to `(u,t)` pairs. - `ConstantInterpolation(u,t)` - A piecewise constant interpolation. - + - `SmoothedConstantInterpolation(u,t)` - An integral preserving continuously differentiable approximation of constant interpolation. - `LinearInterpolation(u,t)` - A linear interpolation. - `QuadraticInterpolation(u,t)` - A quadratic interpolation. - `LagrangeInterpolation(u,t,n)` - A Lagrange interpolation of order `n`. diff --git a/docs/src/manual.md b/docs/src/manual.md index 3818a9e3..e599fad9 100644 --- a/docs/src/manual.md +++ b/docs/src/manual.md @@ -6,6 +6,7 @@ QuadraticInterpolation LagrangeInterpolation AkimaInterpolation ConstantInterpolation +SmoothedConstantInterpolation QuadraticSpline CubicSpline BSplineInterpolation diff --git a/docs/src/methods.md b/docs/src/methods.md index 648131eb..ae2e17d3 100644 --- a/docs/src/methods.md +++ b/docs/src/methods.md @@ -77,6 +77,24 @@ A = ConstantInterpolation(u, t, dir = :right) plot(A) ``` +## Smoothed Constant Interpolation + +This function is much like the constant interpolation above, but the transition +between consecutive values is smoothed out so that the function is continuously +differentiable. The smoothing is done in such a way that the integral of this function +is never much off from the same integral of constant interpolation without smoothing (because of the symmetry of the smoothing sections). +The maximum smoothing distance in the `t` direction from the data points can be set +with the keyword argument `d_max`. + +```@example tutorial +A = ConstantInterpolation(u, t) +plot(A) +A = SmoothedConstantInterpolation(u, t; d_max = 10) +plot!(A) +``` + +Note that `u[end]` is ignored. + ## Quadratic Spline This is the quadratic spline. It is a continuously differentiable interpolation diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index cb2bb0aa..7f0f8126 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -290,8 +290,8 @@ end cache_parameters = false, assume_linear_t = 1e-2) It is a method for interpolating constantly with forward fill, with smoothing around the -value transitions to make the curve C1 smooth while the integral never drifts far from -the integral of constant interpolation. +value transitions to make the curve continuously differentiable while the integral never +drifts far from the integral of constant interpolation. In this interpolation type u[end] is ignored. ## Arguments From 40cb55668d11529bd1b248cf29498042daef6fea Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 26 Nov 2024 10:26:56 +0100 Subject: [PATCH 05/27] integrals --- README.md | 1 + docs/src/index.md | 1 + docs/src/methods.md | 2 +- src/integrals.jl | 25 +++++++++++++++++++++++++ src/interpolation_caches.jl | 2 +- 5 files changed, 29 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 0dd6635d..dbf4df79 100644 --- a/README.md +++ b/README.md @@ -47,6 +47,7 @@ In all cases, `u` an `AbstractVector` of values and `t` is an `AbstractVector` o corresponding to `(u,t)` pairs. - `ConstantInterpolation(u,t)` - A piecewise constant interpolation. + - `SmoothedConstantInterpolation(u,t)` - An integral preserving continuously differentiable approximation of constant interpolation. - `LinearInterpolation(u,t)` - A linear interpolation. - `QuadraticInterpolation(u,t)` - A quadratic interpolation. diff --git a/docs/src/index.md b/docs/src/index.md index e74bbd8b..5ce66a6d 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -17,6 +17,7 @@ In all cases, `u` an `AbstractVector` of values and `t` is an `AbstractVector` o corresponding to `(u,t)` pairs. - `ConstantInterpolation(u,t)` - A piecewise constant interpolation. + - `SmoothedConstantInterpolation(u,t)` - An integral preserving continuously differentiable approximation of constant interpolation. - `LinearInterpolation(u,t)` - A linear interpolation. - `QuadraticInterpolation(u,t)` - A quadratic interpolation. diff --git a/docs/src/methods.md b/docs/src/methods.md index ae2e17d3..9a7b54ff 100644 --- a/docs/src/methods.md +++ b/docs/src/methods.md @@ -80,7 +80,7 @@ plot(A) ## Smoothed Constant Interpolation This function is much like the constant interpolation above, but the transition -between consecutive values is smoothed out so that the function is continuously +between consecutive values is smoothed out so that the function is continuously differentiable. The smoothing is done in such a way that the integral of this function is never much off from the same integral of constant interpolation without smoothing (because of the symmetry of the smoothing sections). The maximum smoothing distance in the `t` direction from the data points can be set diff --git a/src/integrals.jl b/src/integrals.jl index 2bdb6408..a054acb0 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -56,6 +56,31 @@ function _integral( end end +function _integral(A::SmoothedConstantInterpolation{<:AbstractVector}, + idx::Number, t1::Number, t2::Number) + d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx) + + bound_lower = A.t[idx] + d_lower + bound_upper = A.t[idx + 1] - d_upper + + out = A.u[idx] * (t2 - t1) + + if t1 < bound_lower + t2_ = min(t2, bound_lower) + out -= c_lower * d_lower * + (((t2_ - A.t[idx]) / d_lower - 1)^3 - ((t1 - A.t[idx]) / d_lower - 1)^3) / 3 + end + + if t2 > bound_upper + t1_ = max(t1, bound_upper) + out += c_upper * d_upper * + ((1 - (A.t[idx + 1] - t2) / d_upper)^3 - + (1 - (A.t[idx + 1] - t1_) / d_upper)^3) / 3 + end + + out +end + function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}}, idx::Number, t1::Number, t2::Number) α, β = get_parameters(A, idx) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 7f0f8126..be0fe9a4 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -290,7 +290,7 @@ end cache_parameters = false, assume_linear_t = 1e-2) It is a method for interpolating constantly with forward fill, with smoothing around the -value transitions to make the curve continuously differentiable while the integral never +value transitions to make the curve continuously differentiable while the integral never drifts far from the integral of constant interpolation. In this interpolation type u[end] is ignored. ## Arguments From 2b1c4f803fde3085a31788758914fba56584beeb Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 26 Nov 2024 12:43:36 +0100 Subject: [PATCH 06/27] Add tests --- src/derivatives.jl | 5 +++++ src/integrals.jl | 7 +++++++ src/interpolation_caches.jl | 7 ++++--- src/interpolation_methods.jl | 8 ++++++++ src/parameter_caches.jl | 8 ++++---- test/derivative_tests.jl | 9 ++++++++- test/integral_tests.jl | 14 ++++++++++++++ test/interpolation_tests.jl | 13 +++++++++++++ test/parameter_tests.jl | 8 ++++++++ 9 files changed, 71 insertions(+), 8 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index 198e0cb6..b4922c9d 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -23,6 +23,11 @@ function _derivative(A::SmoothedConstantInterpolation{<:AbstractVector}, t::Numb idx = get_idx(A, t, iguess) d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx) + # Fix extrapolation behavior as constant for now + if t <= first(A.t) || t >= last(A.t) + return zero(c_upper / oneunit(t)) + end + if (t - A.t[idx]) < d_lower -2c_lower * ((t - A.t[idx]) / d_lower - 1) / d_lower elseif (A.t[idx + 1] - t) < d_upper diff --git a/src/integrals.jl b/src/integrals.jl index a054acb0..4f1487e6 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -65,6 +65,13 @@ function _integral(A::SmoothedConstantInterpolation{<:AbstractVector}, out = A.u[idx] * (t2 - t1) + # Fix extrapolation behavior as constant for now + if t1 <= first(A.t) + t1 = first(A.t) + elseif t2 >= last(A.t) + t2 = last(A.t) + end + if t1 < bound_lower t2_ = min(t2, bound_lower) out -= c_lower * d_lower * diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index be0fe9a4..66b11aa2 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -308,12 +308,12 @@ drifts far from the integral of constant interpolation. In this interpolation ty 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. """ -struct SmoothedConstantInterpolation{uType, tType, dmaxType, IType, pType, T, N} <: +struct SmoothedConstantInterpolation{uType, tType, IType, dType, cType, dmaxType, T, N} <: AbstractInterpolation{T, N} u::uType t::tType I::IType - p::SmoothedConstantParameterCache{uType, tType} + p::SmoothedConstantParameterCache{dType, cType} d_max::dmaxType extrapolate::Bool iguesser::Guesser{tType} @@ -323,7 +323,8 @@ struct SmoothedConstantInterpolation{uType, tType, dmaxType, IType, pType, T, N} u, t, I, p, d_max, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) N = get_output_dim(u) - new{typeof(u), typeof(t), typeof(d_max), typeof(I), typeof(p.d), eltype(u), N}( + new{typeof(u), typeof(t), typeof(I), typeof(p.d), + typeof(p.c), typeof(d_max), eltype(u), N}( u, t, I, p, d_max, extrapolate, Guesser(t), cache_parameters, linear_lookup) end end diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 0615504d..1870ea28 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -143,6 +143,14 @@ end function _interpolate(A::SmoothedConstantInterpolation{<:AbstractVector}, t::Number, iguess) idx = get_idx(A, t, iguess) d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx) + + # Fix extrapolation behavior as constant for now + if t <= first(A.t) + return first(A.u) + elseif t >= last(A.t) + return A.u[end - 1] + end + out = A.u[idx] if (t - A.t[idx]) < d_lower diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index 84785680..912f497e 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -32,9 +32,9 @@ function linear_interpolation_parameters(u::AbstractArray{T, N}, t, idx) where { return slope end -struct SmoothedConstantParameterCache{uType, tType} - d::tType - c::uType +struct SmoothedConstantParameterCache{dType, cType} + d::dType + c::cType end function SmoothedConstantParameterCache(u, t, cache_parameters, d_max) @@ -51,7 +51,7 @@ end function smoothed_linear_interpolation_parameters(u, t, d_max, idx) # TODO: Add support for making periodic extrapolation smooth if isone(idx) || (idx == length(t)) - zero(eltype(t)), zero(eltype(u)) + zero(one(eltype(t))) / 2, zero(one(eltype(u)) / 2) else min(t[idx] - t[idx - 1], t[idx + 1] - t[idx], 2d_max) / 2, (u[idx] - u[idx - 1]) / 2 end diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 37595a4f..15731fd2 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -27,7 +27,7 @@ function test_derivatives(method; args = [], kwargs = [], name::String) # Interpolation time points for _t in t[2:(end - 1)] if func isa BSplineInterpolation || func isa BSplineApprox || - func isa CubicHermiteSpline + func isa CubicHermiteSpline || func isa SmoothedConstantInterpolation fdiff = forward_fdm(5, 1; geom = true)(func, _t) fdiff2 = forward_fdm(5, 1; geom = true)(t -> derivative(func, t), _t) else @@ -140,6 +140,13 @@ end @test all(derivative.(Ref(A), t2 .+ 0.1) .== 0.0) end +@testset "SmoothedConstantInterpolation" begin + u = [5.5, 2.7, 5.1, 3.0] + t = [2.5, 5.6, 6.3, 8.9] + test_derivatives(SmoothedConstantInterpolation; args = [u, t], + name = "Smoothed constant interpolation") +end + @testset "Quadratic Spline" begin u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] diff --git a/test/integral_tests.jl b/test/integral_tests.jl index 3d59a6c0..126a2496 100644 --- a/test/integral_tests.jl +++ b/test/integral_tests.jl @@ -65,6 +65,20 @@ end name = "Linear Interpolation (Vector) with random points") end +@testset "SmoothedConstantInterpolation" begin + u = [1.0, 4.0, 9.0, 16.0] + t = [1.0, 2.0, 3.0, 4.0] + test_integral( + SmoothedConstantInterpolation; args = [u, t], name = "Smoothed constant interpolation") + + A_constant = ConstantInterpolation(u, t) + I_ref = DataInterpolations.integral(A_constant, first(t), last(t)) + I_smoothed = [DataInterpolations.integral( + SmoothedConstantInterpolation(u, t; d_max), first(t), last(t)) + for d_max in 0.0:0.1:1.0] + @test all(I_smoothed .≈ I_ref) +end + @testset "QuadraticInterpolation" begin u = [1.0, 4.0, 9.0, 16.0] t = [1.0, 2.0, 3.0, 4.0] diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 65b48113..1e8bfabb 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -501,6 +501,19 @@ end @test A(-Inf) == first(u) end +@testset "Smoothed constant Interpolation" begin + test_interpolation_type(SmoothedConstantInterpolation) + u = [0.0, 2.0, 1.0, 3.0] + t = [1.2, 2.5, 5.7, 8.7] + d_max = 0.5 + A = SmoothedConstantInterpolation(u, t; d_max) + test_cached_index(A) + + @test A(1.9) == u[1] + @test A(3.1) == u[2] + @test A(2.5) ≈ (u[1] + u[2]) / 2 +end + @testset "QuadraticSpline Interpolation" begin test_interpolation_type(QuadraticSpline) diff --git a/test/parameter_tests.jl b/test/parameter_tests.jl index 3afe2a17..3030ec21 100644 --- a/test/parameter_tests.jl +++ b/test/parameter_tests.jl @@ -7,6 +7,14 @@ using DataInterpolations @test A.p.slope ≈ [4.0, -2.0, 1.0, 0.0] end +@testset "Smoothed constant Interpolation" begin + u = [1.0, 5.0, 3.0, 4.0, 4.0] + t = collect(1:5) + A = SmoothedConstantInterpolation(u, t; cache_parameters = true) + A.p.d ≈ [0.0, 0.5, 0.5, 0.5, 0.0] + A.p.c ≈ [0.0, 2.0, -1.0, 0.5, 0.0] +end + @testset "Quadratic Interpolation" begin u = [1.0, 5.0, 3.0, 4.0, 4.0] t = collect(1:5) From 84207cf9986a98ac959e734addf3cbc25d1f18b6 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 3 Dec 2024 15:18:30 +0100 Subject: [PATCH 07/27] Merge branch 'master' into SmoothedConstantInterpolation --- .github/workflows/CompatHelper.yml | 2 +- .github/workflows/Downgrade.yml | 2 +- Project.toml | 8 +- docs/make.jl | 3 +- docs/src/extrapolation_methods.md | 90 ++++ docs/src/interface.md | 4 +- ...ataInterpolationsRegularizationToolsExt.jl | 70 ++- src/DataInterpolations.jl | 36 +- src/derivatives.jl | 64 ++- src/integral_inverses.jl | 10 +- src/integrals.jl | 138 +++++- src/interpolation_caches.jl | 451 +++++++++++++----- src/interpolation_methods.jl | 60 ++- src/interpolation_utils.jl | 31 +- test/derivative_tests.jl | 22 +- test/extrapolation_tests.jl | 112 +++++ test/integral_inverse_tests.jl | 2 +- test/integral_tests.jl | 12 +- test/interface.jl | 4 +- test/interpolation_tests.jl | 162 ++++--- test/parameter_tests.jl | 13 + test/regularization.jl | 5 +- test/runtests.jl | 7 +- test/zygote_tests.jl | 6 +- 24 files changed, 1028 insertions(+), 286 deletions(-) create mode 100644 docs/src/extrapolation_methods.md create mode 100644 test/extrapolation_tests.jl diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index f95ed99c..08dea66f 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.github/workflows/CompatHelper.yml @@ -13,7 +13,7 @@ jobs: CompatHelper: runs-on: ubuntu-latest steps: - - uses: julia-actions/setup-julia@9b79636afcfb07ab02c256cede01fe2db6ba808c + - uses: julia-actions/setup-julia@5c9647d97b78a5debe5164e9eec09d653d29bd71 with: version: 1.3 - name: Pkg.add("CompatHelper") diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml index 3030e057..2c0bef40 100644 --- a/.github/workflows/Downgrade.yml +++ b/.github/workflows/Downgrade.yml @@ -28,7 +28,7 @@ jobs: - windows-latest steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2.6.0 + - uses: julia-actions/setup-julia@v2.6.1 with: version: ${{ matrix.version }} - uses: julia-actions/julia-downgrade-compat@v1 diff --git a/Project.toml b/Project.toml index 80cd1153..91ae9fed 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" @@ -40,7 +42,8 @@ RegularizationTools = "0.6" SafeTestsets = "0.1" StableRNGs = "1" Symbolics = "5.29, 6" -Test = "1" +Test = "1.10" +Unitful = "1.21.1" Zygote = "0.6.70" julia = "1.10" @@ -57,7 +60,8 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["Aqua", "BenchmarkTools", "SafeTestsets", "ChainRulesCore", "Optim", "RegularizationTools", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "Symbolics", "Zygote"] +test = ["Aqua", "BenchmarkTools", "SafeTestsets", "ChainRulesCore", "Optim", "RegularizationTools", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "Symbolics", "Unitful", "Zygote"] 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..7344b9fe --- /dev/null +++ b/docs/src/extrapolation_methods.md @@ -0,0 +1,90 @@ +# Extrapolation methods + +We will use the following interpolation to demonstrate the various extrapolation methods. + +```@example tutorial +using DataInterpolations, Plots + +u = [0.86, 0.65, 0.44, 0.76, 0.73] +t = [0.0022, 0.68, 1.41, 2.22, 2.46] +t_eval_left = range(-1, first(t), length = 25) +t_eval_right = range(last(t), 3.5, length = 25) +A = QuadraticSpline(u, t) +plot(A) +``` + +Extrapolation behavior can be set left and right of the data simultaneously with the `extension` keyword, or left and right separately with the `extrapolation_left` and `extrapolation_right` keywords respectively. + +## `ExtrapolationType.None` + +This extrapolation type will throw an error when the input `t` is beyond the data in the specified direction. + +## `ExtrapolationType.Constant` + +This extrapolation type extends the interpolation with the boundary values of the data `u`. + +```@example tutorial +A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.Constant) +plot(A) +plot!(t_eval_left, A.(t_eval_left); label = "extrapolation left") +plot!(t_eval_right, A.(t_eval_right); label = "extrapolation right") +``` + +## `ExtrapolationType.Linear` + +This extrapolation type extends the interpolation with a linear continuation of the interpolation, making it $C^1$ smooth at the data boundaries. + +```@example tutorial +A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.Linear) +plot(A) +plot!(t_eval_left, A.(t_eval_left); label = "extrapolation left") +plot!(t_eval_right, A.(t_eval_right); label = "extrapolation right") +``` + +## `ExtrapolationType.Extension` + +This extrapolation type extends the interpolation with a continuation of the expression for the interpolation at the boundary intervals for maximum smoothness. + +```@example tutorial +A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.Extension) +plot(A) +plot!(t_eval_left, A.(t_eval_left); label = "extrapolation down") +plot!(t_eval_right, A.(t_eval_right); label = "extrapolation up") +``` + +## `ExtrapolationType.Periodic` + +this extrapolation type extends the interpolation such that `A(t + T) == A(t)` for all `t`, where the period is given by `T = last(A.t) - first(A.t)`. + +```@example tutorial +T = last(A.t) - first(A.t) +t_eval_left = range(first(t) - 2T, first(t), length = 100) +t_eval_right = range(last(t), last(t) + 2T, length = 100) +A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.Periodic) +plot(A) +plot!(t_eval_left, A.(t_eval_left); label = "extrapolation down") +plot!(t_eval_right, A.(t_eval_right); label = "extrapolation up") +``` + +## `ExtrapolationType.Reflective` + +this extrapolation type extends the interpolation such that `A(t_ + t) == A(t_ - t)` for all `t_, t` such that `(t_ - first(A.t)) % T == 0` and `0 < t < T`, where `T = last(A.t) - first(A.t)`. + +```@example tutorial +A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.Reflective) +plot(A) +plot!(t_eval_left, A.(t_eval_left); label = "extrapolation down") +plot!(t_eval_right, A.(t_eval_right); label = "extrapolation up") +``` + +## Mixed extrapolation + +You can also have different extrapolation types left and right of the data. + +```@example tutorial +A = QuadraticSpline(u, t; extrapolation_left = ExtrapolationType.Reflective, + extrapolation_right = ExtrapolationType.Periodic) +plot(A) +plot!(t_eval_left, A.(t_eval_left); label = "extrapolation left") +plot!(t_eval_right, A.(t_eval_right); label = "extrapolation right") +``` diff --git a/docs/src/interface.md b/docs/src/interface.md index cfc9ed0b..738287ee 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..cf0e1de3 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 1da0f240..d25e8a17 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 Periodic Reflective + include("parameter_caches.jl") include("interpolation_caches.jl") include("interpolation_utils.jl") @@ -57,31 +60,43 @@ 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 @@ -89,7 +104,7 @@ export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, AkimaInterpolation, ConstantInterpolation, SmoothedConstantInterpolation, QuadraticSpline, CubicSpline, BSplineInterpolation, BSplineApprox, CubicHermiteSpline, PCHIPInterpolation, QuinticHermiteSpline, - LinearInterpolationIntInv, ConstantInterpolationIntInv + LinearInterpolationIntInv, ConstantInterpolationIntInv, ExtrapolationType # added for RegularizationSmooth, JJS 11/27/21 ### Regularization data smoothing and interpolation @@ -105,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, @@ -116,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, @@ -129,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 b4922c9d..b5c7e409 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -1,15 +1,63 @@ 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 + +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])) + elseif extrapolation_left == ExtrapolationType.Extension + iguess = A.iguesser + (order == 1) ? _derivative(A, t, iguess) : + ForwardDiff.derivative(t -> begin + _derivative(A, t, iguess) + end, t) + elseif extrapolation_left == ExtrapolationType.Periodic + t_, _ = transformation_periodic(A, t) + derivative(A, t_, order) + else + # extrapolation_left == ExtrapolationType.Reflective + t_, n = transformation_reflective(A, t) + isodd(n) ? -derivative(A, t_, order) : derivative(A, t_, order) + end +end - return if order == 1 - _derivative(A, t, iguess) - elseif order == 2 +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])) + elseif extrapolation_right == ExtrapolationType.Extension + iguess = A.iguesser + (order == 1) ? _derivative(A, t, iguess) : ForwardDiff.derivative(t -> begin _derivative(A, t, iguess) end, t) + elseif extrapolation_right == ExtrapolationType.Periodic + t_, _ = transformation_periodic(A, t) + derivative(A, t_, order) else - throw(DerivativeNotFoundError()) + # extrapolation_right == ExtrapolationType.Reflective + t_, n = transformation_reflective(A, t) + iseven(n) ? -derivative(A, t_, order) : derivative(A, t_, order) end end @@ -45,7 +93,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]) @@ -79,7 +126,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]) @@ -131,12 +177,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 4f1487e6..6019867f 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -1,39 +1,135 @@ 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) + 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[max(1, idx2 - 1)] - A.I[idx1] - return if t1 == t2 - zero(total) + total += A.I[idx2 - 1] + if idx1 > 0 + total -= A.I[idx1] + end + else + for idx in (idx1 + 1):(idx2 - 1) + total += _integral(A, idx, A.t[idx], A.t[idx + 1]) + end + 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)) + elseif extrapolation_left == ExtrapolationType.Periodic + t_, n = transformation_periodic(A, t) + out = -integral(A, t_) + if !iszero(n) + out -= n * integral(A, first(A.t), last(A.t)) + end + out + else + # extrapolation_left == ExtrapolationType.Reflective + t_, n = transformation_reflective(A, t) + out = if isodd(n) + -integral(A, t_, last(A.t)) 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 + -integral(A, t_) + end + if !iszero(n) + out -= n * integral(A, first(A.t), last(A.t)) + end + out + 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) + elseif extrapolation_right == ExtrapolationType.Periodic + t_, n = transformation_periodic(A, t) + out = integral(A, first(A.t), t_) + if !iszero(n) + out += n * integral(A, first(A.t), last(A.t)) end + out 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) + # extrapolation_right == ExtrapolationType.Reflective + t_, n = transformation_reflective(A, t) + out = if iseven(n) + integral(A, t_, last(A.t)) + else + integral(A, t_) + end + if !iszero(n) + out += n * integral(A, first(A.t), last(A.t)) end - total + out end end diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 66b11aa2..106351a1 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` + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + - `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` + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + - `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` + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + - `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` + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + - `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` + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + - `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,27 +334,36 @@ 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 """ @@ -301,7 +382,13 @@ drifts far from the integral of constant interpolation. In this interpolation ty ## Keyword Arguments - `d_max`: the maximum distance in `t` from the data points the smoothing is allowed to reach. - - `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` + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + - `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 behavior for evenly-distributed abscissae. Alternatively, a numerical threshold may be specified @@ -315,33 +402,39 @@ struct SmoothedConstantInterpolation{uType, tType, IType, dType, cType, dmaxType I::IType p::SmoothedConstantParameterCache{dType, cType} d_max::dmaxType - extrapolate::Bool + extrapolation_left::ExtrapolationType.T + extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function SmoothedConstantInterpolation( - u, t, I, p, d_max, extrapolate, cache_parameters, assume_linear_t) + u, t, I, p, d_max, 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.d), typeof(p.c), typeof(d_max), eltype(u), N}( - u, t, I, p, d_max, extrapolate, Guesser(t), cache_parameters, linear_lookup) + u, t, I, p, d_max, extrapolation_left, extrapolation_right, Guesser(t), cache_parameters, linear_lookup) end end -function SmoothedConstantInterpolation(u, t; d_max = Inf, extrapolate = false, - cache_parameters = false, assume_linear_t = 1e-2) +function SmoothedConstantInterpolation(u, t; d_max = Inf, 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 = SmoothedConstantParameterCache(u, t, cache_parameters, d_max) A = SmoothedConstantInterpolation( - u, t, nothing, p, d_max, extrapolate, cache_parameters, assume_linear_t) + u, t, nothing, p, d_max, extrapolation_left, extrapolation_right, cache_parameters, assume_linear_t) I = cumulative_integral(A, cache_parameters) SmoothedConstantInterpolation( - u, t, I, p, d_max, extrapolate, cache_parameters, assume_linear_t) + u, t, I, p, d_max, 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. @@ -353,7 +446,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` + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + - `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 @@ -369,12 +468,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), @@ -385,7 +486,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 @@ -394,9 +496,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) @@ -407,15 +513,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) @@ -436,13 +548,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. @@ -454,7 +569,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` + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + - `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 @@ -469,11 +590,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₁), @@ -484,7 +607,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 @@ -493,10 +617,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) @@ -517,16 +644,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) @@ -551,15 +683,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) @@ -577,13 +715,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. @@ -598,8 +739,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` + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + - `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. @@ -615,7 +762,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, @@ -627,7 +775,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) @@ -641,7 +790,8 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, scType, T, N} <: sc, pVecType, knotVecType, - extrapolate, + extrapolation_left, + extrapolation_right, Guesser(t), linear_lookup ) @@ -649,7 +799,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.") @@ -714,12 +869,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.") @@ -787,11 +948,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). @@ -808,7 +971,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` + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + - `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 @@ -826,7 +995,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, @@ -839,7 +1009,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) @@ -855,7 +1026,8 @@ struct BSplineApprox{uType, tType, pType, kType, cType, scType, T, N} <: sc, pVecType, knotVecType, - extrapolate, + extrapolation_left, + extrapolation_right, Guesser(t), linear_lookup ) @@ -863,7 +1035,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.") @@ -949,12 +1126,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.") @@ -1045,10 +1228,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. @@ -1060,7 +1245,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` + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + - `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 @@ -1074,33 +1265,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), @@ -1113,22 +1314,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` + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + - `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. @@ -1141,7 +1348,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` + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + - `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 @@ -1156,29 +1369,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 1870ea28..0882e67b 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,7 +1,55 @@ 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)) + elseif extrapolation_left == ExtrapolationType.Extension + _interpolate(A, t, A.iguesser) + elseif extrapolation_left == ExtrapolationType.Periodic + t_, _ = transformation_periodic(A, t) + _interpolate(A, t_, A.iguesser) + else + # extrapolation_left == ExtrapolationType.Reflective + t_, _ = transformation_reflective(A, t) + _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)) + elseif extrapolation_right == ExtrapolationType.Extension + _interpolate(A, t, A.iguesser) + elseif extrapolation_right == ExtrapolationType.Periodic + t_, _ = transformation_periodic(A, t) + _interpolate(A, t_, A.iguesser) + else + # extrapolation_right == ExtrapolationType.Reflective + t_, _ = transformation_reflective(A, t) + _interpolate(A, t_, A.iguesser) + end end # Linear Interpolation @@ -9,9 +57,9 @@ function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, igues if isnan(t) # For correct derivative with NaN idx = firstindex(A.u) - t1 = t2 = one(eltype(A.t)) - u1 = u2 = one(eltype(A.u)) - slope = t * get_parameters(A, idx) + t1 = t2 = oneunit(eltype(A.t)) + u1 = u2 = oneunit(eltype(A.u)) + slope = t / t * get_parameters(A, idx) else idx = get_idx(A, t, iguess) t1, t2 = A.t[idx], A.t[idx + 1] diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index a3b719a7..f6537e4d 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -311,7 +311,34 @@ 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 + +function transformation_periodic(A, t) + Δt = last(A.t) - first(A.t) + n, t_ = fldmod(t - first(A.t), Δt) + t_ += first(A.t) + (n > 0) && (n -= 1) + t_, n +end + +function transformation_reflective(A, t) + Δt = last(A.t) - first(A.t) + n, t_ = fldmod(t - first(A.t), Δt) + t_ = isodd(n) ? last(A.t) - t_ : first(A.t) + t_ + (n > 0) && (n -= 1) + t_, n end diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 15731fd2..3190e601 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 @@ -239,7 +249,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 @@ -252,7 +262,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 @@ -331,7 +341,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..ad37fc0d --- /dev/null +++ b/test/extrapolation_tests.jl @@ -0,0 +1,112 @@ +using DataInterpolations, Test +using ForwardDiff +using QuadGK + +function test_extrapolation(method, u, t) + @testset "Extrapolation errors" begin + 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 + + for extrapolation_type in instances(ExtrapolationType.T) + (extrapolation_type == ExtrapolationType.None) && continue + @testset "extrapolation type $extrapolation_type" begin + A = method(u, t; extrapolation = extrapolation_type) + + t_eval = first(t) - 1.5 + @test DataInterpolations.derivative(A, t_eval) ≈ + ForwardDiff.derivative(A, t_eval) + + t_eval = last(t) + 1.5 + @test DataInterpolations.derivative(A, t_eval) ≈ + ForwardDiff.derivative(A, t_eval) + + T = last(A.t) - first(A.t) + t1 = first(t) - 2.5T + t2 = last(t) + 3.5T + @test DataInterpolations.integral(A, t1, t2) ≈ + quadgk(A, t1, t2; atol = 1e-12, rtol = 1e-12)[1] + end + end +end + +@testset "Linear Interpolation" begin + u = [1.0, 2.0] + t = [1.0, 2.0] + + test_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(QuadraticInterpolation, 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..494e4943 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 126a2496..a528ed66 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..b0e48d06 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 1e8bfabb..99afa59e 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -3,12 +3,14 @@ using FindFirstFunctions: searchsortedfirstcorrelated using StableRNGs using Optim, ForwardDiff using BenchmarkTools +using Unitful 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 +32,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 +42,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 +55,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 +66,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 +76,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 +84,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 +92,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 +100,15 @@ 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 + @test isnan(A(3.5)) + @test isnan(A(4.0)) + + u = [0.0, 1.0, 2.0, NaN] + 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 +118,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) @@ -134,16 +144,22 @@ end @test @inferred(A(R64)) === A(R64) end + # NaN time value for Unitful arrays: issue #365 + t = (0:3)u"s" # Unitful quantities + u = [0, -2, -1, -2]u"m" + A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension) + @test isnan(A(NaN * u"s")) + # 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 +172,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 +187,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 +201,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 +216,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 +255,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 +268,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 +276,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 +286,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 +346,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 +377,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 +391,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 +404,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 +421,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 +434,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 +450,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 +463,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 +479,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 +492,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 +507,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 @@ -520,7 +541,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) @@ -536,14 +557,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) @@ -552,12 +573,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 @@ -566,7 +587,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 @@ -585,7 +606,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 @@ -594,7 +615,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 @@ -605,12 +626,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 @@ -649,12 +670,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) @@ -669,12 +691,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 @@ -726,12 +749,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 @@ -772,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 = 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 @@ -800,7 +824,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/parameter_tests.jl b/test/parameter_tests.jl index 3030ec21..978a1d51 100644 --- a/test/parameter_tests.jl +++ b/test/parameter_tests.jl @@ -1,10 +1,18 @@ using DataInterpolations +function test_cached_integration(method, args...) + A_c = method(args...; cache_parameters = true) + A_nc = method(args...; cache_parameters = false) + @test DataInterpolations.integral(A_c, last(A_c.t)) ≈ + DataInterpolations.integral(A_nc, last(A_nc.t)) +end + @testset "Linear Interpolation" begin u = [1.0, 5.0, 3.0, 4.0, 4.0] t = collect(1:5) A = LinearInterpolation(u, t; cache_parameters = true) @test A.p.slope ≈ [4.0, -2.0, 1.0, 0.0] + test_cached_integration(LinearInterpolation, u, t) end @testset "Smoothed constant Interpolation" begin @@ -21,6 +29,7 @@ end A = QuadraticInterpolation(u, t; cache_parameters = true) @test A.p.α ≈ [-3.0, 1.5, -0.5, -0.5] @test A.p.β ≈ [7.0, -3.5, 1.5, 0.5] + test_cached_integration(QuadraticInterpolation, u, t) end @testset "Quadratic Spline" begin @@ -29,6 +38,7 @@ end A = QuadraticSpline(u, t; cache_parameters = true) @test A.p.α ≈ [-9.5, 3.5, -0.5, -0.5] @test A.p.β ≈ [13.5, -5.5, 1.5, 0.5] + test_cached_integration(QuadraticSpline, u, t) end @testset "Cubic Spline" begin @@ -37,6 +47,7 @@ end A = CubicSpline(u, t; cache_parameters = true) @test A.p.c₁ ≈ [6.839285714285714, 1.642857142857143, 4.589285714285714, 4.0] @test A.p.c₂ ≈ [1.0, 6.839285714285714, 1.642857142857143, 4.589285714285714] + test_cached_integration(CubicSpline, u, t) end @testset "Cubic Hermite Spline" begin @@ -46,6 +57,7 @@ end A = CubicHermiteSpline(du, u, t; cache_parameters = true) @test A.p.c₁ ≈ [-1.0, -5.0, -5.0, -8.0] @test A.p.c₂ ≈ [0.0, 13.0, 12.0, 9.0] + test_cached_integration(CubicHermiteSpline, du, u, t) end @testset "Quintic Hermite Spline" begin @@ -57,4 +69,5 @@ end @test A.p.c₁ ≈ [-1.0, -6.5, -8.0, -10.0] @test A.p.c₂ ≈ [1.0, 19.5, 20.0, 19.0] @test A.p.c₃ ≈ [1.5, -37.5, -37.0, -26.5] + test_cached_integration(QuinticHermiteSpline, ddu, du, u, t) end diff --git a/test/regularization.jl b/test/regularization.jl index fcf7580d..eb8abca0 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..33b9ab1c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,7 +7,8 @@ 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 Tests" include("extrapolation_tests.jl") @safetestset "Online Tests" include("online_tests.jl") -@safetestset "Regularization Smoothing" include("regularization.jl") -@safetestset "Show methods" include("show.jl") -@safetestset "Zygote support" include("zygote_tests.jl") +@safetestset "Regularization Smoothing Tests" include("regularization.jl") +@safetestset "Show methods Tests" include("show.jl") +@safetestset "Zygote support Tests" include("zygote_tests.jl") diff --git a/test/zygote_tests.jl b/test/zygote_tests.jl index 7af735af..da4bf381 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 From 77a5875c57b88f715b5dff931db635c926810c07 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 3 Dec 2024 15:24:58 +0100 Subject: [PATCH 08/27] Formatting --- src/interpolation_caches.jl | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 106351a1..40c76e24 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -408,28 +408,33 @@ struct SmoothedConstantInterpolation{uType, tType, IType, dType, cType, dmaxType cache_parameters::Bool linear_lookup::Bool function SmoothedConstantInterpolation( - u, t, I, p, d_max, extrapolation_left, extrapolation_right, cache_parameters, assume_linear_t) + u, t, I, p, d_max, 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.d), typeof(p.c), typeof(d_max), eltype(u), N}( - u, t, I, p, d_max, extrapolation_left, extrapolation_right, Guesser(t), cache_parameters, linear_lookup) + u, t, I, p, d_max, extrapolation_left, extrapolation_right, + Guesser(t), cache_parameters, linear_lookup) end end -function SmoothedConstantInterpolation(u, t; d_max = Inf, 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) +function SmoothedConstantInterpolation( + u, t; d_max = Inf, 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 = SmoothedConstantParameterCache(u, t, cache_parameters, d_max) A = SmoothedConstantInterpolation( - u, t, nothing, p, d_max, extrapolation_left, extrapolation_right, cache_parameters, assume_linear_t) + u, t, nothing, p, d_max, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) I = cumulative_integral(A, cache_parameters) SmoothedConstantInterpolation( - u, t, I, p, d_max, extrapolation_left, extrapolation_right, cache_parameters, assume_linear_t) + u, t, I, p, d_max, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) end """ From 8395f226336b16141136e6848307f54546c7d6d6 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 3 Dec 2024 16:04:45 +0100 Subject: [PATCH 09/27] Make periodic extrapolation smooth --- src/interpolation_methods.jl | 7 ------- src/interpolation_utils.jl | 6 +++--- src/parameter_caches.jl | 16 ++++++++++++---- 3 files changed, 15 insertions(+), 14 deletions(-) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 0882e67b..ccbc668f 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -192,13 +192,6 @@ function _interpolate(A::SmoothedConstantInterpolation{<:AbstractVector}, t::Num idx = get_idx(A, t, iguess) d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx) - # Fix extrapolation behavior as constant for now - if t <= first(A.t) - return first(A.u) - elseif t >= last(A.t) - return A.u[end - 1] - end - out = A.u[idx] if (t - A.t[idx]) < d_lower diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index f6537e4d..bb1f7492 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -215,9 +215,9 @@ function get_parameters(A::SmoothedConstantInterpolation, idx) c_upper = A.p.c[idx + 1] d_lower, d_upper, c_lower, c_upper else - d_lower, c_lower = smoothed_linear_interpolation_parameters(A.u, A.t, A.d_max, idx) - d_upper, c_upper = smoothed_linear_interpolation_parameters( - A.u, A.t, A.d_max, idx + 1) + d_lower, c_lower = smoothed_constant_interpolation_parameters(A.u, A.t, A.d_max, idx, A.extrapolation_left, A.extrapolation_right) + d_upper, c_upper = smoothed_constant_interpolation_parameters( + A.u, A.t, A.d_max, idx + 1, A.extrapolation_left, A.extrapolation_right) d_lower, d_upper, c_lower, c_upper end end diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index 912f497e..ca4c5773 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -39,7 +39,7 @@ end function SmoothedConstantParameterCache(u, t, cache_parameters, d_max) if cache_parameters - parameters = smoothed_linear_interpolation_parameters.( + parameters = smoothed_constant_interpolation_parameters.( Ref(u), Ref(t), d_max, eachindex(t)) d, c = collect.(eachrow(stack(collect.(parameters)))) SmoothedConstantParameterCache(d, c) @@ -48,10 +48,18 @@ function SmoothedConstantParameterCache(u, t, cache_parameters, d_max) end end -function smoothed_linear_interpolation_parameters(u, t, d_max, idx) - # TODO: Add support for making periodic extrapolation smooth +function smoothed_constant_interpolation_parameters(u, t, d_max, idx, extrapolation_left, extrapolation_right) if isone(idx) || (idx == length(t)) - zero(one(eltype(t))) / 2, zero(one(eltype(u)) / 2) + # If extrapolation is periodic, make the transition differentiable + if extrapolation_left == extrapolation_right == ExtrapolationType.Periodic + if isone(idx) + min(t[end] - t[end - 1], t[2] - t[1], 2d_max) / 2, (u[1] - u[end - 1]) / 2 + else + min(t[end] - t[end - 1], t[2] - t[1], 2d_max) / 2, (u[1] - u[end - 1]) / 2 + end + else + zero(one(eltype(t)) / 2), zero(one(eltype(u)) / 2) + end else min(t[idx] - t[idx - 1], t[idx + 1] - t[idx], 2d_max) / 2, (u[idx] - u[idx - 1]) / 2 end From 321e70f4f0c2aa50666d88321d23d02e9c0a4d0a Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 3 Dec 2024 16:10:40 +0100 Subject: [PATCH 10/27] Fix tests --- src/interpolation_caches.jl | 5 +++-- src/interpolation_utils.jl | 3 ++- src/parameter_caches.jl | 8 +++++--- test/parameter_tests.jl | 4 ++-- 4 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 40c76e24..e9ad5a81 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -384,7 +384,7 @@ drifts far from the integral of constant interpolation. In this interpolation ty - `d_max`: the maximum distance in `t` from the data points the smoothing is allowed to reach. - `extrapolation`: The extrapolation type applied left and right of the data. Possible options are `ExtrapolationType.None` (default), `ExtrapolationType.Constant`, `ExtrapolationType.Linear` - `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`. + `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` (also made smooth at the boundaries) and `ExtrapolationType.Reflective`. - `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 @@ -427,7 +427,8 @@ function SmoothedConstantInterpolation( extrapolation_left, extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) - p = SmoothedConstantParameterCache(u, t, cache_parameters, d_max) + p = SmoothedConstantParameterCache( + u, t, cache_parameters, d_max, extrapolation_left, extrapolation_right) A = SmoothedConstantInterpolation( u, t, nothing, p, d_max, extrapolation_left, extrapolation_right, cache_parameters, assume_linear_t) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index bb1f7492..eff36bb8 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -215,7 +215,8 @@ function get_parameters(A::SmoothedConstantInterpolation, idx) c_upper = A.p.c[idx + 1] d_lower, d_upper, c_lower, c_upper else - d_lower, c_lower = smoothed_constant_interpolation_parameters(A.u, A.t, A.d_max, idx, A.extrapolation_left, A.extrapolation_right) + d_lower, c_lower = smoothed_constant_interpolation_parameters( + A.u, A.t, A.d_max, idx, A.extrapolation_left, A.extrapolation_right) d_upper, c_upper = smoothed_constant_interpolation_parameters( A.u, A.t, A.d_max, idx + 1, A.extrapolation_left, A.extrapolation_right) d_lower, d_upper, c_lower, c_upper diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index ca4c5773..c5e1daf1 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -37,10 +37,11 @@ struct SmoothedConstantParameterCache{dType, cType} c::cType end -function SmoothedConstantParameterCache(u, t, cache_parameters, d_max) +function SmoothedConstantParameterCache( + u, t, cache_parameters, d_max, extrapolation_left, extrapolation_right) if cache_parameters parameters = smoothed_constant_interpolation_parameters.( - Ref(u), Ref(t), d_max, eachindex(t)) + Ref(u), Ref(t), d_max, eachindex(t), extrapolation_left, extrapolation_right) d, c = collect.(eachrow(stack(collect.(parameters)))) SmoothedConstantParameterCache(d, c) else @@ -48,7 +49,8 @@ function SmoothedConstantParameterCache(u, t, cache_parameters, d_max) end end -function smoothed_constant_interpolation_parameters(u, t, d_max, idx, extrapolation_left, extrapolation_right) +function smoothed_constant_interpolation_parameters( + u, t, d_max, idx, extrapolation_left, extrapolation_right) if isone(idx) || (idx == length(t)) # If extrapolation is periodic, make the transition differentiable if extrapolation_left == extrapolation_right == ExtrapolationType.Periodic diff --git a/test/parameter_tests.jl b/test/parameter_tests.jl index 978a1d51..40293517 100644 --- a/test/parameter_tests.jl +++ b/test/parameter_tests.jl @@ -19,8 +19,8 @@ end u = [1.0, 5.0, 3.0, 4.0, 4.0] t = collect(1:5) A = SmoothedConstantInterpolation(u, t; cache_parameters = true) - A.p.d ≈ [0.0, 0.5, 0.5, 0.5, 0.0] - A.p.c ≈ [0.0, 2.0, -1.0, 0.5, 0.0] + @test A.p.d ≈ [0.0, 0.5, 0.5, 0.5, 0.0] + @test A.p.c ≈ [0.0, 2.0, -1.0, 0.5, 0.0] end @testset "Quadratic Interpolation" begin From b401b03383acabeb13135dbe1e23e80bf225e391 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 3 Dec 2024 16:29:05 +0100 Subject: [PATCH 11/27] Actually fix all tests --- src/derivatives.jl | 5 ----- src/parameter_caches.jl | 10 ++++------ 2 files changed, 4 insertions(+), 11 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index b5c7e409..3439aa27 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -71,11 +71,6 @@ function _derivative(A::SmoothedConstantInterpolation{<:AbstractVector}, t::Numb idx = get_idx(A, t, iguess) d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx) - # Fix extrapolation behavior as constant for now - if t <= first(A.t) || t >= last(A.t) - return zero(c_upper / oneunit(t)) - end - if (t - A.t[idx]) < d_lower -2c_lower * ((t - A.t[idx]) / d_lower - 1) / d_lower elseif (A.t[idx + 1] - t) < d_upper diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index c5e1daf1..ca1f6c00 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -54,13 +54,11 @@ function smoothed_constant_interpolation_parameters( if isone(idx) || (idx == length(t)) # If extrapolation is periodic, make the transition differentiable if extrapolation_left == extrapolation_right == ExtrapolationType.Periodic - if isone(idx) - min(t[end] - t[end - 1], t[2] - t[1], 2d_max) / 2, (u[1] - u[end - 1]) / 2 - else - min(t[end] - t[end - 1], t[2] - t[1], 2d_max) / 2, (u[1] - u[end - 1]) / 2 - end + min(t[end] - t[end - 1], t[2] - t[1], 2d_max) / 2, (u[1] - u[end - 1]) / 2 else - zero(one(eltype(t)) / 2), zero(one(eltype(u)) / 2) + d = isone(idx) ? min(t[2] - t[1], 2d_max) / 2 : + min(t[end] - t[end - 1], 2d_max) / 2 + d, zero(one(eltype(u)) / 2) end else min(t[idx] - t[idx - 1], t[idx + 1] - t[idx], 2d_max) / 2, (u[idx] - u[idx - 1]) / 2 From b2985a23cb5aaac04547d3209a40cf1d6340a174 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 3 Dec 2024 16:31:33 +0100 Subject: [PATCH 12/27] Actually Actually fix all tests --- test/parameter_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/parameter_tests.jl b/test/parameter_tests.jl index 40293517..6baee961 100644 --- a/test/parameter_tests.jl +++ b/test/parameter_tests.jl @@ -19,7 +19,7 @@ end u = [1.0, 5.0, 3.0, 4.0, 4.0] t = collect(1:5) A = SmoothedConstantInterpolation(u, t; cache_parameters = true) - @test A.p.d ≈ [0.0, 0.5, 0.5, 0.5, 0.0] + @test A.p.d ≈ [0.0, 0.5, 0.5, 0.5, 0.5] @test A.p.c ≈ [0.0, 2.0, -1.0, 0.5, 0.0] end From 52ab2f192376c8d46d99ccc23b8fb29efc997711 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 3 Dec 2024 16:37:08 +0100 Subject: [PATCH 13/27] Actually actually actually fix all tests --- test/parameter_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/parameter_tests.jl b/test/parameter_tests.jl index 6baee961..e7aa6b63 100644 --- a/test/parameter_tests.jl +++ b/test/parameter_tests.jl @@ -19,7 +19,7 @@ end u = [1.0, 5.0, 3.0, 4.0, 4.0] t = collect(1:5) A = SmoothedConstantInterpolation(u, t; cache_parameters = true) - @test A.p.d ≈ [0.0, 0.5, 0.5, 0.5, 0.5] + @test A.p.d ≈ [0.5, 0.5, 0.5, 0.5, 0.5] @test A.p.c ≈ [0.0, 2.0, -1.0, 0.5, 0.0] end From 735cc0563138c5428a1fc37649cb9617d3d90da1 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 9 Dec 2024 11:29:06 +0100 Subject: [PATCH 14/27] Address comments --- src/interpolation_caches.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index e9ad5a81..28a44bab 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -372,7 +372,7 @@ end It is a method for interpolating constantly with forward fill, with smoothing around the value transitions to make the curve continuously differentiable while the integral never -drifts far from the integral of constant interpolation. In this interpolation type u[end] is ignored. +drifts far from the integral of constant interpolation. In this interpolation type, `u[end]` is ignored. ## Arguments @@ -381,7 +381,8 @@ drifts far from the integral of constant interpolation. In this interpolation ty ## Keyword Arguments - - `d_max`: the maximum distance in `t` from the data points the smoothing is allowed to reach. + - `d_max`: Around each time point `tᵢ` there is a continuously differentiable (quadratic) transition between `uᵢ₋₁` and `uᵢ`, + on the interval `[tᵢ - d, tᵢ + d]`. The distance `d` is determined as `d = min((tᵢ - tᵢ₋₁)/2, (tᵢ₊₁ - tᵢ)/2, d_max)`. - `extrapolation`: The extrapolation type applied left and right of the data. Possible options are `ExtrapolationType.None` (default), `ExtrapolationType.Constant`, `ExtrapolationType.Linear` `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` (also made smooth at the boundaries) and `ExtrapolationType.Reflective`. From 6912a07ef98457cc224fad38b8936784290981c5 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 16 Apr 2025 10:24:11 +0200 Subject: [PATCH 15/27] Merge cleanup --- src/interpolation_caches.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 3d2ea540..7584b636 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -391,8 +391,8 @@ drifts far from the integral of constant interpolation. In this interpolation ty 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. """ -struct SmoothedConstantInterpolation{uType, tType, IType, dType, cType, dmaxType, T, N} <: - AbstractInterpolation{T, N} +struct SmoothedConstantInterpolation{uType, tType, IType, dType, cType, dmaxType, T} <: + AbstractInterpolation{T} u::uType t::tType I::IType @@ -407,9 +407,8 @@ struct SmoothedConstantInterpolation{uType, tType, IType, dType, cType, dmaxType u, t, I, p, d_max, 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.d), - typeof(p.c), typeof(d_max), eltype(u), N}( + typeof(p.c), typeof(d_max), eltype(u)}( u, t, I, p, d_max, extrapolation_left, extrapolation_right, Guesser(t), cache_parameters, linear_lookup) end From def265ae6ca07d2e6b5ef87408f99d2f6be7162b Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 16 Apr 2025 11:16:31 +0200 Subject: [PATCH 16/27] Add get_transition_ts --- src/interpolation_utils.jl | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index a3613107..1953f73f 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -388,3 +388,24 @@ function smooth_arc_length_params_2(u_int, uⱼ, uⱼ₊₁) end return δⱼ, short_side_left, Δt_line_seg end + +function get_transition_ts(A::SmoothedConstantInterpolation) + out = similar(A.t, 2 * length(A.t)) + + for idx in 1:(length(A.t) - 1) + d_lower, d_upper, _, _ = get_parameters(A, idx) + if idx == 1 + out[1] = A.t[idx] + out[2] = A.t[idx] + d_lower + else + out[2idx - 1] = A.t[idx] - d_upper + out[2idx] = A.t[idx] + d_lower + end + end + + d_lower, _, _, _ = get_parameters(A, 1) + out[end - 1] = A.t[end] - d_lower + out[end] = A.t[end] + + out +end \ No newline at end of file From e13f0f257dff5d553baacc36007c1b528fe7246f Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 16 Apr 2025 11:22:29 +0200 Subject: [PATCH 17/27] Formatting --- src/interpolation_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 1953f73f..2345c4d6 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -408,4 +408,4 @@ function get_transition_ts(A::SmoothedConstantInterpolation) out[end] = A.t[end] out -end \ No newline at end of file +end From 2203988535b062862704a0a3bd7c494f1a01f5df Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 22 Apr 2025 13:29:11 +0200 Subject: [PATCH 18/27] Update extrapolation behavior --- src/interpolation_caches.jl | 3 ++- src/interpolation_methods.jl | 18 ++++++++++++++++++ src/interpolation_utils.jl | 2 ++ src/parameter_caches.jl | 3 +++ 4 files changed, 25 insertions(+), 1 deletion(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 7584b636..93b60393 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -367,7 +367,8 @@ end It is a method for interpolating constantly with forward fill, with smoothing around the value transitions to make the curve continuously differentiable while the integral never -drifts far from the integral of constant interpolation. In this interpolation type, `u[end]` is ignored. +drifts far from the integral of constant interpolation. `u[end]` is ignored for extrapolation types +`None` (default) and `Periodic`. ## Arguments diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index ae6ece62..c30b0cff 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -74,6 +74,24 @@ function _extrapolate_right(A::ConstantInterpolation, t) end end +function _extrapolate_right(A::SmoothedConstantInterpolation, t) + if A.extrapolation_right == ExtrapolationType.None + throw(RightExtrapolationError()) + elseif A.extrapolation_right in ( + ExtrapolationType.Constant, ExtrapolationType.Extension) + d = min(A.t[end] - A.t[end - 1], 2A.d_max) / 2 + if A.t[end] + d < t + A.u[end] + else + c = (A.u[end] - A.u[end - 1]) / 2 + A.u[end - 1] - c * (((t - A.t[end]) / d)^2 - 2 * ((t - A.t[end]) / d) - 1) + end + + else + _extrapolate_other(A, t, A.extrapolation_right) + end +end + # Linear Interpolation function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, iguess) if isnan(t) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 2345c4d6..b6252a88 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -409,3 +409,5 @@ function get_transition_ts(A::SmoothedConstantInterpolation) out end + +get_transition_ts(A::AbstractInterpolation) = A.t \ No newline at end of file diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index ca1f6c00..c52de5df 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -55,6 +55,9 @@ function smoothed_constant_interpolation_parameters( # If extrapolation is periodic, make the transition differentiable if extrapolation_left == extrapolation_right == ExtrapolationType.Periodic min(t[end] - t[end - 1], t[2] - t[1], 2d_max) / 2, (u[1] - u[end - 1]) / 2 + elseif (idx == length(t)) && (extrapolation_right in ( + ExtrapolationType.Constant, ExtrapolationType.Extension)) + min(t[end] - t[end - 1], 2d_max) / 2, (u[end] - u[end - 1]) / 2 else d = isone(idx) ? min(t[2] - t[1], 2d_max) / 2 : min(t[end] - t[end - 1], 2d_max) / 2 From 2262c6c5e9caf121e632649314e1e6a0b0d6c93c Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 22 Apr 2025 13:35:20 +0200 Subject: [PATCH 19/27] Formatting --- ...ataInterpolationsRegularizationToolsExt.jl | 42 ++++++++----- src/interpolation_caches.jl | 59 ++++++++++++------- src/interpolation_utils.jl | 11 ++-- test/interpolation_tests.jl | 2 + 4 files changed, 76 insertions(+), 38 deletions(-) diff --git a/ext/DataInterpolationsRegularizationToolsExt.jl b/ext/DataInterpolationsRegularizationToolsExt.jl index 766a2503..c396b79e 100644 --- a/ext/DataInterpolationsRegularizationToolsExt.jl +++ b/ext/DataInterpolationsRegularizationToolsExt.jl @@ -74,13 +74,15 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Abstrac extrapolation_left::ExtrapolationType.T = ExtrapolationType.None, extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters::Bool = false) - extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t; check_sorted = t̂, sorted_arg_name = ("third", "t̂")) M = _mapping_matrix(t̂, t) Wls½ = LA.diagm(sqrt.(wls)) Wr½ = LA.diagm(sqrt.(wr)) - û, λ, Aitp = _reg_smooth_solve( + û, λ, + Aitp = _reg_smooth_solve( u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right, cache_parameters) RegularizationSmooth( @@ -100,7 +102,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, d::Int = 2; extrapolation_left::ExtrapolationType.T = ExtrapolationType.None, extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters::Bool = false) - extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) t̂ = t @@ -108,7 +111,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, d::Int = 2; 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( + û, λ, + Aitp = _reg_smooth_solve( u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right, cache_parameters) RegularizationSmooth(u, @@ -137,14 +141,16 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Abstrac extrapolation_left::ExtrapolationType.T = ExtrapolationType.None, extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters::Bool = false) - extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t; check_sorted = t̂, sorted_arg_name = ("third", "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( + û, λ, + Aitp = _reg_smooth_solve( u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right, cache_parameters) RegularizationSmooth(u, @@ -174,14 +180,16 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Abstrac extrapolation_left::ExtrapolationType.T = ExtrapolationType.None, extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters::Bool = false) - extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t; check_sorted = t̂, sorted_arg_name = ("third", "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( + û, λ, + Aitp = _reg_smooth_solve( u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right, cache_parameters) RegularizationSmooth(u, @@ -212,7 +220,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing extrapolation_left::ExtrapolationType.T = ExtrapolationType.None, extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters::Bool = false) - extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) t̂ = t @@ -220,7 +229,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing 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( + û, λ, + Aitp = _reg_smooth_solve( u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right, cache_parameters) RegularizationSmooth(u, @@ -251,7 +261,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing extrapolation_left::ExtrapolationType.T = ExtrapolationType.None, extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters::Bool = false) - extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) t̂ = t @@ -259,7 +270,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing M = Array{Float64}(LA.I, N, N) Wls½ = LA.diagm(sqrt.(wls)) Wr½ = LA.diagm(sqrt.(wr)) - û, λ, Aitp = _reg_smooth_solve( + û, λ, + Aitp = _reg_smooth_solve( u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right, cache_parameters) RegularizationSmooth(u, @@ -289,7 +301,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing extrapolation_left::ExtrapolationType.T = ExtrapolationType.None, extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters::Bool = false) - extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) t̂ = t @@ -298,7 +311,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( + û, λ, + Aitp = _reg_smooth_solve( u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right, cache_parameters) RegularizationSmooth(u, diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 93b60393..15d6467b 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -50,7 +50,8 @@ function LinearInterpolation( 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_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) p = LinearParameterCache(u, t, cache_parameters) @@ -120,7 +121,8 @@ function QuadraticInterpolation( 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_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) linear_lookup = seems_linear(assume_linear_t, t) @@ -190,7 +192,8 @@ function LagrangeInterpolation( extrapolation::ExtrapolationType.T = ExtrapolationType.None, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None, extrapolation_right::ExtrapolationType.T = ExtrapolationType.None) - extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) if n != length(t) - 1 @@ -263,7 +266,8 @@ function AkimaInterpolation( 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_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) linear_lookup = seems_linear(assume_linear_t, t) @@ -350,7 +354,8 @@ function ConstantInterpolation( 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_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) A = ConstantInterpolation( @@ -420,7 +425,8 @@ function SmoothedConstantInterpolation( 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_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) p = SmoothedConstantParameterCache( @@ -502,7 +508,8 @@ function QuadraticSpline( 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_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) @@ -527,7 +534,8 @@ function QuadraticSpline( extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters = false, assume_linear_t = 1e-2) where {uType <: AbstractVector} - extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) @@ -621,7 +629,8 @@ function CubicSpline(u::AbstractVector{<:Number}, 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_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) - 1 @@ -655,7 +664,8 @@ function CubicSpline(u::AbstractArray{T, N}, 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 {T, N} - extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) - 1 @@ -693,7 +703,8 @@ function CubicSpline( 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_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) - 1 @@ -799,7 +810,8 @@ function BSplineInterpolation( 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_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) @@ -875,7 +887,8 @@ function BSplineInterpolation( extrapolation_left::ExtrapolationType.T = ExtrapolationType.None, extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, assume_linear_t = 1e-2) - extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) @@ -1034,7 +1047,8 @@ function BSplineApprox( 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_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) @@ -1131,7 +1145,8 @@ function BSplineApprox( 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_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) n = length(t) @@ -1204,8 +1219,9 @@ function BSplineApprox( spline_coefficients!(view(sc, i, :), d, k, p[i]) end for k in 2:(n - 1) - q[ax_u..., k] = u[ax_u..., k] - sc[k, 1] * u[ax_u..., 1] - - sc[k, h] * u[ax_u..., end] + q[ax_u..., + k] = u[ax_u..., k] - sc[k, 1] * u[ax_u..., 1] - + sc[k, h] * u[ax_u..., end] end Q = Array{T, N}(undef, size(u)[1:(end - 1)]..., h - 2) for i in 2:(h - 1) @@ -1280,7 +1296,8 @@ function CubicHermiteSpline( 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_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) linear_lookup = seems_linear(assume_linear_t, t) @@ -1385,7 +1402,8 @@ function QuinticHermiteSpline( 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_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) u, t = munge_data(u, t) linear_lookup = seems_linear(assume_linear_t, t) @@ -1736,7 +1754,8 @@ function SmoothArcLengthInterpolation( t[j + 1] = t[j] + Δt_circle_seg + Δt_line_seg end - extrapolation_left, extrapolation_right = munge_extrapolation( + extrapolation_left, + extrapolation_right = munge_extrapolation( extrapolation, extrapolation_left, extrapolation_right) linear_lookup = seems_linear(assume_linear_t, t) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index b6252a88..d8658e96 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -191,7 +191,8 @@ function cumulative_integral(A::AbstractInterpolation{<:Number}, cache_parameter Base.require_one_based_indexing(A.u) idxs = cache_parameters ? (1:(length(A.t) - 1)) : (1:0) return cumsum(_integral(A, idx, t1, t2) - for (idx, t1, t2) in zip(idxs, @view(A.t[begin:(end - 1)]), @view(A.t[(begin + 1):end]))) + for (idx, t1, t2) in + zip(idxs, @view(A.t[begin:(end - 1)]), @view(A.t[(begin + 1):end]))) end function get_parameters(A::LinearInterpolation, idx) @@ -210,9 +211,11 @@ function get_parameters(A::SmoothedConstantInterpolation, idx) c_upper = A.p.c[idx + 1] d_lower, d_upper, c_lower, c_upper else - d_lower, c_lower = smoothed_constant_interpolation_parameters( + d_lower, + c_lower = smoothed_constant_interpolation_parameters( A.u, A.t, A.d_max, idx, A.extrapolation_left, A.extrapolation_right) - d_upper, c_upper = smoothed_constant_interpolation_parameters( + d_upper, + c_upper = smoothed_constant_interpolation_parameters( A.u, A.t, A.d_max, idx + 1, A.extrapolation_left, A.extrapolation_right) d_lower, d_upper, c_lower, c_upper end @@ -410,4 +413,4 @@ function get_transition_ts(A::SmoothedConstantInterpolation) out end -get_transition_ts(A::AbstractInterpolation) = A.t \ No newline at end of file +get_transition_ts(A::AbstractInterpolation) = A.t diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 4bb50b9d..77eb0355 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -1024,6 +1024,7 @@ end ut1 = Float32[0.1, 0.2, 0.3, 0.4, 0.5] ut2 = Float64[0.1, 0.2, 0.3, 0.4, 0.5] for u in (ut1, ut2), t in (ut1, ut2) + interp = @inferred(LinearInterpolation(ut1, ut2)) for xs in (u, t) ys = @inferred(interp(xs)) @@ -1106,6 +1107,7 @@ f_cubic_spline = c -> square(CubicSpline, c) iszero_allocations(u, t) = iszero(@allocated(DataInterpolations.munge_data(u, t))) for T in (String, Union{String, Missing}), dims in 1:3 + _u0 = convert(Array{T}, reshape(u0, ntuple(i -> i == dims ? 3 : 1, dims))) u, t = @inferred(DataInterpolations.munge_data(_u0, t0)) From 66ae918cf37ed316532a363b974465cada5413a4 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 22 Apr 2025 13:57:06 +0200 Subject: [PATCH 20/27] Fix extrapolation derivatives --- src/derivatives.jl | 18 ++++++++++++++++++ src/interpolation_methods.jl | 2 +- 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index 5be1bd14..65aee74b 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -74,6 +74,24 @@ function _extrapolate_derivative_right(A, t, order) end end +function _extrapolate_derivative_right(A::SmoothedConstantInterpolation, t, order) + if A.extrapolation_right == ExtrapolationType.None + throw(RightExtrapolationError()) + elseif A.extrapolation_right in ( + ExtrapolationType.Constant, ExtrapolationType.Extension) + d = min(A.t[end] - A.t[end - 1], 2A.d_max) / 2 + if A.t[end] + d < t + zero(eltype(A.u)) + else + c = (A.u[end] - A.u[end - 1]) / 2 + c * (2((t - A.t[end]) / d) - 2) + end + + else + _extrapolate_other(A, t, A.extrapolation_right) + end +end + function _derivative(A::LinearInterpolation, t::Number, iguess) idx = get_idx(A, t, iguess; idx_shift = -1, ub_shift = -1, side = :first) slope = get_parameters(A, idx) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index c30b0cff..568bc71a 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -84,7 +84,7 @@ function _extrapolate_right(A::SmoothedConstantInterpolation, t) A.u[end] else c = (A.u[end] - A.u[end - 1]) / 2 - A.u[end - 1] - c * (((t - A.t[end]) / d)^2 - 2 * ((t - A.t[end]) / d) - 1) + -c * (((t - A.t[end]) / d)^2 - 2 * ((t - A.t[end]) / d) - 1) end else From 5a7847f4b4e94098d758cbdd39b96bd57a654d98 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 22 Apr 2025 14:00:13 +0200 Subject: [PATCH 21/27] Nit --- src/interpolation_methods.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 568bc71a..c30b0cff 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -84,7 +84,7 @@ function _extrapolate_right(A::SmoothedConstantInterpolation, t) A.u[end] else c = (A.u[end] - A.u[end - 1]) / 2 - -c * (((t - A.t[end]) / d)^2 - 2 * ((t - A.t[end]) / d) - 1) + A.u[end - 1] - c * (((t - A.t[end]) / d)^2 - 2 * ((t - A.t[end]) / d) - 1) end else From 4cdb798684844fa4b20d26ef681439e1e80f965c Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 22 Apr 2025 14:09:40 +0200 Subject: [PATCH 22/27] Use get_transition_ts in derivative testing --- src/interpolation_utils.jl | 3 +-- test/derivative_tests.jl | 4 ++-- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index d8658e96..771fd748 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -191,8 +191,7 @@ function cumulative_integral(A::AbstractInterpolation{<:Number}, cache_parameter Base.require_one_based_indexing(A.u) idxs = cache_parameters ? (1:(length(A.t) - 1)) : (1:0) return cumsum(_integral(A, idx, t1, t2) - for (idx, t1, t2) in - zip(idxs, @view(A.t[begin:(end - 1)]), @view(A.t[(begin + 1):end]))) + for (idx, t1, t2) in zip(idxs, @view(A.t[begin:(end - 1)]), @view(A.t[(begin + 1):end]))) end function get_parameters(A::LinearInterpolation, idx) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index b3451da1..1fc90a29 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -1,7 +1,7 @@ using DataInterpolations, Test using FindFirstFunctions: searchsortedfirstcorrelated using FiniteDifferences -using DataInterpolations: derivative +using DataInterpolations: derivative, get_transition_ts using Symbolics using StableRNGs using RegularizationTools @@ -15,7 +15,7 @@ function test_derivatives(method; args = [], kwargs = [], name::String) [:extrapolation_right => ExtrapolationType.Extension, :extrapolation_left => ExtrapolationType.Extension] func = method(args...; kwargs..., kwargs_extrapolation...) - (; t) = func + t = get_transition_ts(func) trange = collect(range(minimum(t) - 5.0, maximum(t) + 5.0, step = 0.1)) trange_exclude = filter(x -> !in(x, t), trange) @testset "$name" begin From dece43adb79dd54df57c1d7430f6d90fa1b68503 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 22 Apr 2025 14:10:12 +0200 Subject: [PATCH 23/27] Formatting --- src/interpolation_utils.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 771fd748..d8658e96 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -191,7 +191,8 @@ function cumulative_integral(A::AbstractInterpolation{<:Number}, cache_parameter Base.require_one_based_indexing(A.u) idxs = cache_parameters ? (1:(length(A.t) - 1)) : (1:0) return cumsum(_integral(A, idx, t1, t2) - for (idx, t1, t2) in zip(idxs, @view(A.t[begin:(end - 1)]), @view(A.t[(begin + 1):end]))) + for (idx, t1, t2) in + zip(idxs, @view(A.t[begin:(end - 1)]), @view(A.t[(begin + 1):end]))) end function get_parameters(A::LinearInterpolation, idx) From 13ec9cf0159ea8a36ef022ca747f2c6b9716386d Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 23 Apr 2025 16:48:44 +0200 Subject: [PATCH 24/27] Fix derivatives and integrals --- src/derivatives.jl | 2 +- src/integrals.jl | 44 ++++++++++++++++++++++++++++++++++++++ src/interpolation_utils.jl | 12 ++++++----- test/derivative_tests.jl | 22 +++++++++++++------ 4 files changed, 67 insertions(+), 13 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index 65aee74b..2fe08471 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -84,7 +84,7 @@ function _extrapolate_derivative_right(A::SmoothedConstantInterpolation, t, orde zero(eltype(A.u)) else c = (A.u[end] - A.u[end - 1]) / 2 - c * (2((t - A.t[end]) / d) - 2) + -c * (2((t - A.t[end]) / d) - 2) / d end else diff --git a/src/integrals.jl b/src/integrals.jl index d105ee41..f0e4afed 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -135,6 +135,50 @@ function _extrapolate_integral_right(A, t) end end +function _extrapolate_integral_right(A::SmoothedConstantInterpolation, t) + (; extrapolation_right) = A + if extrapolation_right == ExtrapolationType.None + throw(RightExtrapolationError()) + elseif A.extrapolation_right in ( + ExtrapolationType.Constant, ExtrapolationType.Extension) + d = min(A.t[end] - A.t[end - 1], 2A.d_max) / 2 + c = (A.u[end] - A.u[end - 1]) / 2 + + Δt_transition = min(t - A.t[end], d) + Δt_constant = max(0, t - A.t[end] - d) + out = Δt_transition * A.u[end - 1] - + c * + (((Δt_transition / d)^3) / (3 / d) - ((Δt_transition^2) / d) - + Δt_transition) + + Δt_constant * A.u[end] + + elseif extrapolation_right == ExtrapolationType.Linear + slope = derivative(A, last(A.t)) + Δt = t - last(A.t) + (last(A.u) + slope * Δt / 2) * Δt + _extrapolate_other(A, t, A.extrapolation_right) + elseif extrapolation_right == ExtrapolationType.Periodic + t_, n = transformation_periodic(A, t) + out = integral(A, first(A.t), t_) + if !iszero(n) + out += n * integral(A, first(A.t), last(A.t)) + end + out + else + # extrapolation_right == ExtrapolationType.Reflective + t_, n = transformation_reflective(A, t) + out = if iseven(n) + integral(A, t_, last(A.t)) + else + integral(A, t_) + end + if !iszero(n) + out += n * integral(A, first(A.t), last(A.t)) + end + out + end +end + function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}}, idx::Number, t1::Number, t2::Number) slope = get_parameters(A, idx) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index d8658e96..b875d7c8 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -393,16 +393,18 @@ function smooth_arc_length_params_2(u_int, uⱼ, uⱼ₊₁) end function get_transition_ts(A::SmoothedConstantInterpolation) - out = similar(A.t, 2 * length(A.t)) + out = similar(A.t, 3 * length(A.t)) for idx in 1:(length(A.t) - 1) d_lower, d_upper, _, _ = get_parameters(A, idx) if idx == 1 - out[1] = A.t[idx] - out[2] = A.t[idx] + d_lower + out[1] = A.t[1] + out[2] = A.t[1] + out[3] = A.t[1] + d_lower else - out[2idx - 1] = A.t[idx] - d_upper - out[2idx] = A.t[idx] + d_lower + out[3idx - 2] = A.t[idx] - d_upper + out[3idx - 1] = A.t[idx] + out[3idx] = A.t[idx] + d_lower end end diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 1fc90a29..fdbaa8f1 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -15,7 +15,7 @@ function test_derivatives(method; args = [], kwargs = [], name::String) [:extrapolation_right => ExtrapolationType.Extension, :extrapolation_left => ExtrapolationType.Extension] func = method(args...; kwargs..., kwargs_extrapolation...) - t = get_transition_ts(func) + t = filter(_t -> first(func.t) ≤ _t ≤ last(func.t), get_transition_ts(func)) trange = collect(range(minimum(t) - 5.0, maximum(t) + 5.0, step = 0.1)) trange_exclude = filter(x -> !in(x, t), trange) @testset "$name" begin @@ -24,19 +24,23 @@ function test_derivatives(method; args = [], kwargs = [], name::String) cdiff = central_fdm(5, 1; geom = true)(func, _t) adiff = derivative(func, _t) @test isapprox(cdiff, adiff, atol = 1e-8) - adiff2 = derivative(func, _t, 2) - cdiff2 = central_fdm(5, 1; geom = true)(t -> derivative(func, t), _t) - @test isapprox(cdiff2, adiff2, atol = 1e-8) + if !(func isa SmoothedConstantInterpolation) + adiff2 = derivative(func, _t, 2) + cdiff2 = central_fdm(5, 1; geom = true)(t -> derivative(func, t), _t) + @test isapprox(cdiff2, adiff2, atol = 1e-8) + end end func isa SmoothArcLengthInterpolation && return - # Interpolation time points + # Interpolation transition points for _t in t[2:(end - 1)] if func isa Union{BSplineInterpolation, BSplineApprox, - CubicHermiteSpline, SmoothedConstantInterpolation} + CubicHermiteSpline} fdiff = forward_fdm(5, 1; geom = true)(func, _t) fdiff2 = forward_fdm(5, 1; geom = true)(t -> derivative(func, t), _t) + elseif func isa SmoothedConstantInterpolation + continue else fdiff = backward_fdm(5, 1; geom = true)(func, _t) fdiff2 = backward_fdm(5, 1; geom = true)(t -> derivative(func, t), _t) @@ -157,9 +161,13 @@ end @testset "SmoothedConstantInterpolation" begin u = [5.5, 2.7, 5.1, 3.0] - t = [2.5, 5.6, 6.3, 8.9] + t = [2.55, 5.62, 6.32, 8.95] test_derivatives(SmoothedConstantInterpolation; args = [u, t], name = "Smoothed constant interpolation") + + A = SmoothedConstantInterpolation( + u, t; extrapolation = ExtrapolationType.Extension) + @test all(_t -> abs(derivative(A, _t)) < 1e-10, setdiff(get_transition_ts(A), t)) end @testset "Quadratic Spline" begin From 8df2e4f1619206e94540ec5129d7a610d9cba0dc Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Thu, 1 May 2025 21:01:46 +0200 Subject: [PATCH 25/27] Clarify when the last u is ignored --- docs/src/methods.md | 2 +- src/interpolation_caches.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/methods.md b/docs/src/methods.md index 8eed38dc..c5a3aeaa 100644 --- a/docs/src/methods.md +++ b/docs/src/methods.md @@ -93,7 +93,7 @@ A = SmoothedConstantInterpolation(u, t; d_max = 10) plot!(A) ``` -Note that `u[end]` is ignored. +Note that `u[end]` is ignored, except when using extrapolation types `Constant` or `Extension`. ## Quadratic Spline diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 15d6467b..47c0773b 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -372,8 +372,8 @@ end It is a method for interpolating constantly with forward fill, with smoothing around the value transitions to make the curve continuously differentiable while the integral never -drifts far from the integral of constant interpolation. `u[end]` is ignored for extrapolation types -`None` (default) and `Periodic`. +drifts far from the integral of constant interpolation. `u[end]` is ignored, +except when using extrapolation types `Constant` or `Extension`. ## Arguments From bfa3bffd788cd189828652d6dae00968fe69b9cb Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Sun, 4 May 2025 11:56:04 +0200 Subject: [PATCH 26/27] Format --- test/interpolation_tests.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 63be5e6a..b219fb3f 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -1141,7 +1141,8 @@ end xvals = rand(rng, 100) @test Symbolics.substitute(ci(0.425), Dict(x => xvals)) == xvals[42] - @test Symbolics.substitute(li(0.425), Dict(x => xvals)) == xvals[42] + 0.5*(xvals[43] - xvals[42]) + @test Symbolics.substitute(li(0.425), Dict(x => xvals)) == + xvals[42] + 0.5*(xvals[43] - xvals[42]) @variables dx[1:100] @test_nowarn chs = CubicHermiteSpline(dx, x, t) From 8de752a284ecb8bc35cfdc4e4fba8c435c846564 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Sun, 4 May 2025 13:28:42 +0200 Subject: [PATCH 27/27] Format with JuliaFormatter v1 --- src/interpolation_caches.jl | 4 ++-- src/interpolation_utils.jl | 3 +-- test/interpolation_tests.jl | 12 +++++------- 3 files changed, 8 insertions(+), 11 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 47c0773b..2f2f1d37 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -1220,8 +1220,8 @@ function BSplineApprox( end for k in 2:(n - 1) q[ax_u..., - k] = u[ax_u..., k] - sc[k, 1] * u[ax_u..., 1] - - sc[k, h] * u[ax_u..., end] + k] = u[ax_u..., k] - sc[k, 1] * u[ax_u..., 1] - + sc[k, h] * u[ax_u..., end] end Q = Array{T, N}(undef, size(u)[1:(end - 1)]..., h - 2) for i in 2:(h - 1) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index b875d7c8..9a3aec88 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -191,8 +191,7 @@ function cumulative_integral(A::AbstractInterpolation{<:Number}, cache_parameter Base.require_one_based_indexing(A.u) idxs = cache_parameters ? (1:(length(A.t) - 1)) : (1:0) return cumsum(_integral(A, idx, t1, t2) - for (idx, t1, t2) in - zip(idxs, @view(A.t[begin:(end - 1)]), @view(A.t[(begin + 1):end]))) + for (idx, t1, t2) in zip(idxs, @view(A.t[begin:(end - 1)]), @view(A.t[(begin + 1):end]))) end function get_parameters(A::LinearInterpolation, idx) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index b219fb3f..4553e15e 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -1025,7 +1025,6 @@ end ut1 = Float32[0.1, 0.2, 0.3, 0.4, 0.5] ut2 = Float64[0.1, 0.2, 0.3, 0.4, 0.5] for u in (ut1, ut2), t in (ut1, ut2) - interp = @inferred(LinearInterpolation(ut1, ut2)) for xs in (u, t) ys = @inferred(interp(xs)) @@ -1108,7 +1107,6 @@ f_cubic_spline = c -> square(CubicSpline, c) iszero_allocations(u, t) = iszero(@allocated(DataInterpolations.munge_data(u, t))) for T in (String, Union{String, Missing}), dims in 1:3 - _u0 = convert(Array{T}, reshape(u0, ntuple(i -> i == dims ? 3 : 1, dims))) u, t = @inferred(DataInterpolations.munge_data(_u0, t0)) @@ -1137,12 +1135,12 @@ end li = LinearInterpolation(x, t) @test isequal(ci(0.425), x[42]) - @test isequal(li(0.425), x[42] + 0.5*(x[43] - x[42])) + @test isequal(li(0.425), x[42] + 0.5 * (x[43] - x[42])) xvals = rand(rng, 100) @test Symbolics.substitute(ci(0.425), Dict(x => xvals)) == xvals[42] @test Symbolics.substitute(li(0.425), Dict(x => xvals)) == - xvals[42] + 0.5*(xvals[43] - xvals[42]) + xvals[42] + 0.5 * (xvals[43] - xvals[42]) @variables dx[1:100] @test_nowarn chs = CubicHermiteSpline(dx, x, t) @@ -1150,7 +1148,7 @@ end @test_nowarn li = LagrangeInterpolation(x, t) @test_nowarn cs = CubicSpline(x, t) - @test_throws Exception ai = AkimaInterpolation(x, t) - @test_throws Exception bsi = BSplineInterpolation(x, t, 3, :ArcLen, :Average) - @test_throws Exception pc = PCHIPInterpolation(x, t) + @test_throws Exception ai=AkimaInterpolation(x, t) + @test_throws Exception bsi=BSplineInterpolation(x, t, 3, :ArcLen, :Average) + @test_throws Exception pc=PCHIPInterpolation(x, t) end