diff --git a/README.md b/README.md index 519724f5..a5b723b7 100644 --- a/README.md +++ b/README.md @@ -48,6 +48,7 @@ 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 2fa2da4d..36cb305f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -18,6 +18,7 @@ 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 2c946b1a..f11c11d4 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 5d45ab5d..c5a3aeaa 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, except when using extrapolation types `Constant` or `Extension`. + ## Quadratic Spline This is the quadratic spline. It is a continuously differentiable interpolation 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/DataInterpolations.jl b/src/DataInterpolations.jl index 73467ce0..9b9dd8af 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -119,9 +119,10 @@ _output_size(u::AbstractVector) = size(first(u)) _output_size(u::AbstractArray) = Base.front(size(u)) export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, - AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline, - BSplineInterpolation, BSplineApprox, CubicHermiteSpline, PCHIPInterpolation, - QuinticHermiteSpline, SmoothArcLengthInterpolation, LinearInterpolationIntInv, + AkimaInterpolation, ConstantInterpolation, SmoothedConstantInterpolation, + QuadraticSpline, CubicSpline, BSplineInterpolation, BSplineApprox, + CubicHermiteSpline, PCHIPInterpolation, QuinticHermiteSpline, + SmoothArcLengthInterpolation, LinearInterpolationIntInv, ConstantInterpolationIntInv, ExtrapolationType export output_dim, output_size diff --git a/src/derivatives.jl b/src/derivatives.jl index 4314559e..2fe08471 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -74,12 +74,43 @@ 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) / d + 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) 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/integrals.jl b/src/integrals.jl index 3fd180b0..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) @@ -154,6 +198,38 @@ 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) + + # 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 * + (((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 2bc29534..2f2f1d37 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( @@ -361,6 +366,80 @@ function ConstantInterpolation( 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 continuously differentiable while the integral never +drifts far from the integral of constant interpolation. `u[end]` is ignored, +except when using extrapolation types `Constant` or `Extension`. + +## Arguments + + - `u`: data points. + - `t`: time points. + +## Keyword Arguments + + - `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`. + - `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 + 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} <: + AbstractInterpolation{T} + u::uType + t::tType + I::IType + p::SmoothedConstantParameterCache{dType, cType} + d_max::dmaxType + 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, extrapolation_left, + extrapolation_right, cache_parameters, assume_linear_t) + linear_lookup = seems_linear(assume_linear_t, t) + new{typeof(u), typeof(t), typeof(I), typeof(p.d), + 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 +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) + 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, extrapolation_left, extrapolation_right) + A = SmoothedConstantInterpolation( + 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) +end + """ QuadraticSpline(u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.None, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None, extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters = false) @@ -429,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) @@ -454,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) @@ -548,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 @@ -582,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 @@ -620,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 @@ -726,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) @@ -802,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) @@ -961,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) @@ -1058,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) @@ -1131,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) @@ -1207,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) @@ -1312,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) @@ -1663,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_methods.jl b/src/interpolation_methods.jl index 5d0c1e3e..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) @@ -186,7 +204,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 @@ -209,6 +227,22 @@ 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 1ce54b7c..9a3aec88 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -202,6 +202,24 @@ 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_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 + function get_parameters(A::QuadraticInterpolation, idx) if A.cache_parameters A.p.α[idx], A.p.β[idx] @@ -372,3 +390,28 @@ 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, 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[1] + out[2] = A.t[1] + out[3] = A.t[1] + d_lower + else + out[3idx - 2] = A.t[idx] - d_upper + out[3idx - 1] = A.t[idx] + out[3idx] = 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 + +get_transition_ts(A::AbstractInterpolation) = A.t diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index 18b46ae5..31206829 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -32,6 +32,42 @@ function linear_interpolation_parameters(u::AbstractArray{T, N}, t, idx) where { return slope end +struct SmoothedConstantParameterCache{dType, cType} + d::dType + c::cType +end + +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), extrapolation_left, extrapolation_right) + d, c = collect.(eachrow(stack(collect.(parameters)))) + SmoothedConstantParameterCache(d, c) + else + SmoothedConstantParameterCache(eltype(t)[], eltype(u)[]) + end +end + +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 + 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 + 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 + end +end + struct QuadraticParameterCache{pType} α::pType β::pType diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 433e600a..fdbaa8f1 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 = 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,18 +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} + if func isa Union{BSplineInterpolation, BSplineApprox, + 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) @@ -154,6 +159,17 @@ 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.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 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 0927c26c..68d1b8ec 100644 --- a/test/integral_tests.jl +++ b/test/integral_tests.jl @@ -83,6 +83,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 2a124507..4553e15e 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -619,6 +619,19 @@ end @test typeof(itp(t)) === typeof(itp.(t)) === Vector{Int} 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) @@ -1122,11 +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]) + @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) @@ -1134,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 diff --git a/test/parameter_tests.jl b/test/parameter_tests.jl index 1f8ff574..e7aa6b63 100644 --- a/test/parameter_tests.jl +++ b/test/parameter_tests.jl @@ -15,6 +15,14 @@ end test_cached_integration(LinearInterpolation, u, t) 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) + @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 + @testset "Quadratic Interpolation" begin u = [1.0, 5.0, 3.0, 4.0, 4.0] t = collect(1:5)