From 80edc15695817f6035853de1ca3e3b423f60e9fd Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 25 Nov 2024 16:20:37 +0100 Subject: [PATCH 1/6] 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 2/6] 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 3/6] 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 4/6] 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 5/6] 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 6/6] 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)