From c9a1e94fa1c05c3cf5fcd040ecf5a0653f81a2a6 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Thu, 1 Aug 2024 20:09:17 +0200 Subject: [PATCH 1/7] Fix problem in get_idx when f = +/-Inf --- src/interpolation_utils.jl | 10 +++++++--- test/interpolation_tests.jl | 7 +++++++ 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index ccd3db63..b8587533 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -119,15 +119,19 @@ end function get_idx(A::AbstractInterpolation, t, iguess; lb = 1, ub_shift = -1, idx_shift = 0, side = :last) + tvec = A.t iguess = if hasfield(typeof(A), :linear_lookup) && A.linear_lookup f = (t - first(A.t)) / (last(A.t) - first(A.t)) - i_0, i_f = firstindex(A.t), lastindex(A.t) - round(typeof(firstindex(A.t)), f * (i_f - i_0) + i_0) + if isinf(f) + f > 0 ? lastindex(A.t) : firstindex(A.t) + else + i_0, i_f = firstindex(A.t), lastindex(A.t) + round(typeof(firstindex(A.t)), f * (i_f - i_0) + i_0) + end else iguess end - tvec = A.t ub = length(tvec) + ub_shift return if side == :last clamp(searchsortedlastcorrelated(tvec, t, iguess) + idx_shift, lb, ub) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 1d4e2e6b..d1291ce5 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -492,6 +492,13 @@ end A = ConstantInterpolation(u, t) @test_throws DataInterpolations.ExtrapolationError A(-1.0) @test_throws DataInterpolations.ExtrapolationError A(11.0) + + # Test extrapolation with infs with regularily 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) + @test A(Inf) == last(u) + @test A(-Inf) == first(u) end @testset "QuadraticSpline Interpolation" begin From f73066652f2af1b4962bd2cfa441ebc2d90274b4 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Thu, 1 Aug 2024 20:14:04 +0200 Subject: [PATCH 2/7] Fix typo --- test/interpolation_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index d1291ce5..8135a5f4 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -493,7 +493,7 @@ end @test_throws DataInterpolations.ExtrapolationError A(-1.0) @test_throws DataInterpolations.ExtrapolationError A(11.0) - # Test extrapolation with infs with regularily spaced t + # 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) From 75be6a98a89dcc9684c6142649e5058ed3ab9ecd Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 2 Aug 2024 15:24:22 +0200 Subject: [PATCH 3/7] Use Guesser from FindFirstFunctions --- ext/DataInterpolationsChainRulesCoreExt.jl | 4 +- src/DataInterpolations.jl | 7 ++-- src/derivatives.jl | 13 ++----- src/integral_inverses.jl | 10 ++--- src/interpolation_caches.jl | 44 +++++++++++----------- src/interpolation_methods.jl | 7 +--- src/interpolation_utils.jl | 12 ------ test/derivative_tests.jl | 6 +-- test/interpolation_tests.jl | 25 ++---------- 9 files changed, 45 insertions(+), 83 deletions(-) diff --git a/ext/DataInterpolationsChainRulesCoreExt.jl b/ext/DataInterpolationsChainRulesCoreExt.jl index 312e67e3..20c803a3 100644 --- a/ext/DataInterpolationsChainRulesCoreExt.jl +++ b/ext/DataInterpolationsChainRulesCoreExt.jl @@ -52,7 +52,7 @@ end function u_tangent(A::LinearInterpolation, t, Δ) out = zero(A.u) - idx = get_idx(A, t, A.idx_prev[]) + idx = get_idx(A, t, A.iguesser(t)) t_factor = (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) out[idx] = Δ * (one(eltype(out)) - t_factor) out[idx + 1] = Δ * t_factor @@ -61,7 +61,7 @@ end function u_tangent(A::QuadraticInterpolation, t, Δ) out = zero(A.u) - i₀, i₁, i₂ = _quad_interp_indices(A, t, A.idx_prev[]) + i₀, i₁, i₂ = _quad_interp_indices(A, t, A.iguesser(t)) t₀ = A.t[i₀] t₁ = A.t[i₁] t₂ = A.t[i₂] diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl index 7f44c878..783a796c 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -8,7 +8,7 @@ using LinearAlgebra, RecipesBase using PrettyTables using ForwardDiff import FindFirstFunctions: searchsortedfirstcorrelated, searchsortedlastcorrelated, - bracketstrictlymontonic + Guesser include("parameter_caches.jl") include("interpolation_caches.jl") @@ -22,9 +22,8 @@ include("online.jl") include("show.jl") (interp::AbstractInterpolation)(t::Number) = _interpolate(interp, t) -function (interp::AbstractInterpolation)(t::Number, i::Integer) - interp.idx_prev[] = i - _interpolate(interp, t) +function (interp::AbstractInterpolation)(t::Number, iguess::Integer) + _interpolate(interp, t; iguess) end function (interp::AbstractInterpolation)(t::AbstractVector) diff --git a/src/derivatives.jl b/src/derivatives.jl index e7cb97fd..f8971b9e 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -1,16 +1,11 @@ -function derivative(A, t, order = 1) +function derivative(A, t, order = 1; iguess = A.iguesser(t)) ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) - iguess = A.idx_prev[] - + return if order == 1 - val, idx = _derivative(A, t, iguess) - A.idx_prev[] = idx - val + _derivative(A, t, iguess)[1] elseif order == 2 ForwardDiff.derivative(t -> begin - val, idx = _derivative(A, t, iguess) - A.idx_prev[] = idx - val + _derivative(A, t, iguess)[1] end, t) else throw(DerivativeNotFoundError()) diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index aebee9d5..d4d36e4f 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -18,7 +18,7 @@ invert_integral(A::AbstractInterpolation) = throw(IntegralInverseNotFoundError() _integral(A::AbstractIntegralInverseInterpolation, idx, t) = throw(IntegralNotFoundError()) function _derivative(A::AbstractIntegralInverseInterpolation, t::Number, iguess) - inv(A.itp(A(t))), A.idx_prev[] + inv(A.itp(A(t))), A.iguesser(t) end """ @@ -38,11 +38,11 @@ struct LinearInterpolationIntInv{uType, tType, itpType, T} <: u::uType t::tType extrapolate::Bool - idx_prev::Base.RefValue{Int} + iguesser::Guesser{tType} itp::itpType function LinearInterpolationIntInv(u, t, A) new{typeof(u), typeof(t), typeof(A), eltype(u)}( - u, t, A.extrapolate, Ref(1), A) + u, t, A.extrapolate, Guesser(t), A) end end @@ -84,11 +84,11 @@ struct ConstantInterpolationIntInv{uType, tType, itpType, T} <: u::uType t::tType extrapolate::Bool - idx_prev::Base.RefValue{Int} + iguesser::Guesser{tType} itp::itpType function ConstantInterpolationIntInv(u, t, A) new{typeof(u), typeof(t), typeof(A), eltype(u)}( - u, t, A.extrapolate, Ref(1), A + u, t, A.extrapolate, Guesser(t), A ) end end diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 71b2b042..d7dce336 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -25,13 +25,13 @@ struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolati I::IType p::LinearParameterCache{pType} extrapolate::Bool - idx_prev::Base.RefValue{Int} + iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function LinearInterpolation(u, t, I, p, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}( - u, t, I, p, extrapolate, Ref(1), cache_parameters, linear_lookup) + u, t, I, p, extrapolate, Guesser(t), cache_parameters, linear_lookup) end end @@ -73,7 +73,7 @@ struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpol p::QuadraticParameterCache{pType} mode::Symbol extrapolate::Bool - idx_prev::Base.RefValue{Int} + iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function QuadraticInterpolation( @@ -82,7 +82,7 @@ struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpol error("mode should be :Forward or :Backward for QuadraticInterpolation") linear_lookup = seems_linear(assume_linear_t, t) new{typeof(u), typeof(t), typeof(I), typeof(p.l₀), eltype(u)}( - u, t, I, p, mode, extrapolate, Ref(1), cache_parameters, linear_lookup) + u, t, I, p, mode, extrapolate, Guesser(t), cache_parameters, linear_lookup) end end @@ -124,7 +124,7 @@ struct LagrangeInterpolation{uType, tType, T, bcacheType} <: bcache::bcacheType idxs::Vector{Int} extrapolate::Bool - idx_prev::Base.RefValue{Int} + iguesser::Guesser{tType} function LagrangeInterpolation(u, t, n, extrapolate) bcache = zeros(eltype(u[1]), n + 1) idxs = zeros(Int, n + 1) @@ -135,7 +135,7 @@ struct LagrangeInterpolation{uType, tType, T, bcacheType} <: bcache, idxs, extrapolate, - Ref(1) + Guesser(t) ) end end @@ -178,7 +178,7 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <: c::cType d::dType extrapolate::Bool - idx_prev::Base.RefValue{Int} + iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function AkimaInterpolation( @@ -192,7 +192,7 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <: c, d, extrapolate, - Ref(1), + Guesser(t), cache_parameters, linear_lookup ) @@ -258,14 +258,14 @@ struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T} p::Nothing dir::Symbol # indicates if value to the $dir should be used for the interpolation extrapolate::Bool - idx_prev::Base.RefValue{Int} + iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function ConstantInterpolation( u, t, I, dir, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) new{typeof(u), typeof(t), typeof(I), eltype(u)}( - u, t, I, nothing, dir, extrapolate, Ref(1), cache_parameters, linear_lookup) + u, t, I, nothing, dir, extrapolate, Guesser(t), cache_parameters, linear_lookup) end end @@ -309,7 +309,7 @@ struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <: d::dType z::zType extrapolate::Bool - idx_prev::Base.RefValue{Int} + iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function QuadraticSpline( @@ -324,7 +324,7 @@ struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <: d, z, extrapolate, - Ref(1), + Guesser(t), cache_parameters, linear_lookup ) @@ -410,7 +410,7 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInter h::hType z::zType extrapolate::Bool - idx_prev::Base.RefValue{Int} + iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function CubicSpline(u, t, I, p, h, z, extrapolate, cache_parameters, assume_linear_t) @@ -423,7 +423,7 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInter h, z, extrapolate, - Ref(1), + Guesser(t), cache_parameters, linear_lookup ) @@ -520,7 +520,7 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <: pVecType::Symbol knotVecType::Symbol extrapolate::Bool - idx_prev::Base.RefValue{Int} + iguesser::Guesser{tType} linear_lookup::Bool function BSplineInterpolation(u, t, @@ -544,7 +544,7 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <: pVecType, knotVecType, extrapolate, - Ref(1), + Guesser(t), linear_lookup ) end @@ -656,7 +656,7 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <: pVecType::Symbol knotVecType::Symbol extrapolate::Bool - idx_prev::Base.RefValue{Int} + iguesser::Guesser{tType} linear_lookup::Bool function BSplineApprox(u, t, @@ -683,7 +683,7 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <: pVecType, knotVecType, extrapolate, - Ref(1), + Guesser(t), linear_lookup ) end @@ -806,14 +806,14 @@ struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInte I::IType p::CubicHermiteParameterCache{pType} extrapolate::Bool - idx_prev::Base.RefValue{Int} + iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function CubicHermiteSpline( du, u, t, I, p, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u)}( - du, u, t, I, p, extrapolate, Ref(1), cache_parameters, linear_lookup) + du, u, t, I, p, extrapolate, Guesser(t), cache_parameters, linear_lookup) end end @@ -887,7 +887,7 @@ struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <: I::IType p::QuinticHermiteParameterCache{pType} extrapolate::Bool - idx_prev::Base.RefValue{Int} + iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function QuinticHermiteSpline( @@ -895,7 +895,7 @@ struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <: linear_lookup = seems_linear(assume_linear_t, t) new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(ddu), typeof(p.c₁), eltype(u)}( - ddu, du, u, t, I, p, extrapolate, Ref(1), cache_parameters, linear_lookup) + ddu, du, u, t, I, p, extrapolate, Guesser(t), cache_parameters, linear_lookup) end end diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 69a5cf94..7d7235b3 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,10 +1,7 @@ -function _interpolate(A, t) +function _interpolate(A, t; iguess = A.iguesser(t)) ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) - idx_guess = A.idx_prev[] - val, idx_prev = _interpolate(A, t, idx_guess) - A.idx_prev[] = idx_prev - return val + return _interpolate(A, t, iguess)[1] end # Linear Interpolation diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index b8587533..a54bdcb4 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -120,18 +120,6 @@ end function get_idx(A::AbstractInterpolation, t, iguess; lb = 1, ub_shift = -1, idx_shift = 0, side = :last) tvec = A.t - iguess = if hasfield(typeof(A), :linear_lookup) && - A.linear_lookup - f = (t - first(A.t)) / (last(A.t) - first(A.t)) - if isinf(f) - f > 0 ? lastindex(A.t) : firstindex(A.t) - else - i_0, i_f = firstindex(A.t), lastindex(A.t) - round(typeof(firstindex(A.t)), f * (i_f - i_0) + i_0) - end - else - iguess - end ub = length(tvec) + ub_shift return if side == :last clamp(searchsortedlastcorrelated(tvec, t, iguess) + idx_shift, lb, ub) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 50abe4ac..ea002ab5 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -39,9 +39,9 @@ function test_derivatives(method; args = [], kwargs = [], name::String) @test isapprox(fdiff, adiff, atol = 1e-8) @test isapprox(fdiff2, adiff2, atol = 1e-8) # Cached index - if hasproperty(func, :idx_prev) - @test abs(func.idx_prev[] - - searchsortedfirstcorrelated(func.t, _t, func.idx_prev[])) <= 1 + if hasproperty(func, :iguesser) + @test abs(func.iguesser.idx_prev[] - + searchsortedfirstcorrelated(func.t, _t, func.iguesser(_t))) <= 1 end end diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 8135a5f4..3ebc1741 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -9,7 +9,7 @@ function test_interpolation_type(T) @test hasfield(T, :u) @test hasfield(T, :t) @test hasfield(T, :extrapolate) - @test hasfield(T, :idx_prev) + @test hasfield(T, :iguesser) @test !isempty(methods(DataInterpolations._interpolate, (T, Any, Number))) @test !isempty(methods(DataInterpolations._integral, (T, Any, Number))) @test !isempty(methods(DataInterpolations._derivative, (T, Any, Number))) @@ -18,9 +18,9 @@ end function test_cached_index(A) for t in range(first(A.t), last(A.t); length = 2 * length(A.t) - 1) A(t) - idx = searchsortedfirstcorrelated(A.t, t, A.idx_prev[]) - @test abs(A.idx_prev[] - - searchsortedfirstcorrelated(A.t, t, A.idx_prev[])) <= 2 + idx = searchsortedfirstcorrelated(A.t, t, A.iguesser) + @test abs(A.iguesser.idx_prev[] - + searchsortedfirstcorrelated(A.t, t, A.iguesser)) <= 2 end end @@ -808,20 +808,3 @@ f_cubic_spline = c -> square(CubicSpline, c) @test ForwardDiff.derivative(f_quadratic_spline, 4.0) ≈ 8.0 @test ForwardDiff.derivative(f_cubic_spline, 2.0) ≈ 4.0 @test ForwardDiff.derivative(f_cubic_spline, 4.0) ≈ 8.0 - -@testset "Linear lookup" begin - for N in (100, 1_000, 10_000, 100_000, 1_000_000) # Interpolant size - seed = 1234 - t = collect(LinRange(0, 1, N)) # collect to avoid fast LinRange dispatch - y = rand(N) - A = LinearInterpolation(y, t) - A_fallback = LinearInterpolation(copy(y), copy(t); assume_linear_t = false) - @test A.linear_lookup - @test !(A_fallback.linear_lookup) - - n_samples = 1_000 - b_linear = @benchmark $A(rand()) samples=n_samples - b_fallback = @benchmark $A_fallback(rand()) samples=n_samples - @test mean(b_linear.times) < mean(b_fallback.times) - end -end From 6b25a7bae676656fbe100f4ac610bec23edc4f92 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 2 Aug 2024 15:36:35 +0200 Subject: [PATCH 4/7] Stop having _interpolate and _derivative return the index --- src/derivatives.jl | 28 ++++++++++++------------ src/interpolation_methods.jl | 42 ++++++++++++++++++------------------ 2 files changed, 35 insertions(+), 35 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index f8971b9e..940125a0 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -1,11 +1,11 @@ function derivative(A, t, order = 1; iguess = A.iguesser(t)) ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) - + return if order == 1 - _derivative(A, t, iguess)[1] + _derivative(A, t, iguess) elseif order == 2 ForwardDiff.derivative(t -> begin - _derivative(A, t, iguess)[1] + _derivative(A, t, iguess) end, t) else throw(DerivativeNotFoundError()) @@ -24,7 +24,7 @@ function _derivative(A::QuadraticInterpolation, t::Number, iguess) du₀ = l₀ * (2t - A.t[i₁] - A.t[i₂]) du₁ = l₁ * (2t - A.t[i₀] - A.t[i₂]) du₂ = l₂ * (2t - A.t[i₀] - A.t[i₁]) - return @views @. du₀ + du₁ + du₂, i₀ + return @views @. du₀ + du₁ + du₂ end function _derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number) @@ -96,21 +96,21 @@ function _derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number) end function _derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number, idx) - _derivative(A, t), idx + _derivative(A, t) end function _derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, idx) - _derivative(A, t), idx + _derivative(A, t) end function _derivative(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess) idx = get_idx(A, t, iguess; idx_shift = -1, side = :first) j = min(idx, length(A.c)) # for smooth derivative at A.t[end] wj = t - A.t[idx] - (@evalpoly wj A.b[idx] 2A.c[j] 3A.d[j]), idx + @evalpoly wj A.b[idx] 2A.c[j] 3A.d[j] end function _derivative(A::ConstantInterpolation, t::Number, iguess) - return zero(first(A.u)), iguess + return zero(first(A.u)) end function _derivative(A::ConstantInterpolation{<:AbstractVector}, t::Number) @@ -127,7 +127,7 @@ end function _derivative(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) idx = get_idx(A, t, iguess; lb = 2, ub_shift = 0, side = :first) σ = get_parameters(A, idx - 1) - A.z[idx - 1] + 2σ * (t - A.t[idx - 1]), idx + A.z[idx - 1] + 2σ * (t - A.t[idx - 1]) end # CubicSpline Interpolation @@ -139,7 +139,7 @@ function _derivative(A::CubicSpline{<:AbstractVector}, t::Number, iguess) c₁, c₂ = get_parameters(A, idx) dC = c₁ dD = -c₂ - dI + dC + dD, idx + dI + dC + dD end function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number, iguess) @@ -160,7 +160,7 @@ function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Num ducum += N[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1]) end end - ducum * A.d * scale, idx + ducum * A.d * scale end # BSpline Curve Approx @@ -181,7 +181,7 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig ducum += N[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1]) end end - ducum * A.d * scale, idx + ducum * A.d * scale end # Cubic Hermite Spline @@ -193,7 +193,7 @@ function _derivative( out = A.du[idx] c₁, c₂ = get_parameters(A, idx) out += Δt₀ * (Δt₀ * c₂ + 2(c₁ + Δt₁ * c₂)) - out, idx + out end # Quintic Hermite Spline @@ -206,5 +206,5 @@ function _derivative( c₁, c₂, c₃ = get_parameters(A, idx) out += Δt₀^2 * (3c₁ + (3Δt₁ + Δt₀) * c₂ + (3Δt₁^2 + Δt₀ * 2Δt₁) * c₃) - out, idx + out end diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 7d7235b3..8b3b223d 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,7 +1,7 @@ function _interpolate(A, t; iguess = A.iguesser(t)) ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) - return _interpolate(A, t, iguess)[1] + return _interpolate(A, t, iguess) end # Linear Interpolation @@ -30,14 +30,14 @@ function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, igues end val = oftype(Δu, val) - val, idx + val end function _interpolate(A::LinearInterpolation{<:AbstractMatrix}, t::Number, iguess) idx = get_idx(A, t, iguess) Δt = t - A.t[idx] slope = get_parameters(A, idx) - return A.u[:, idx] + slope * Δt, idx + return A.u[:, idx] + slope * Δt end # Quadratic Interpolation @@ -53,7 +53,7 @@ function _interpolate(A::QuadraticInterpolation, t::Number, iguess) u₀ = l₀ * (t - A.t[i₁]) * (t - A.t[i₂]) u₁ = l₁ * (t - A.t[i₀]) * (t - A.t[i₂]) u₂ = l₂ * (t - A.t[i₀]) * (t - A.t[i₁]) - return u₀ + u₁ + u₂, i₀ + return u₀ + u₁ + u₂ end # Lagrange Interpolation @@ -61,7 +61,7 @@ function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, igu idx = get_idx(A, t, iguess) findRequiredIdxs!(A, t, idx) if A.t[A.idxs[1]] == t - return A.u[A.idxs[1]], idx + return A.u[A.idxs[1]] end N = zero(A.u[1]) D = zero(A.t[1]) @@ -83,14 +83,14 @@ function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, igu D += tmp N += (tmp * A.u[A.idxs[i]]) end - N / D, idx + N / D end function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, iguess) idx = get_idx(A, t, iguess) findRequiredIdxs!(A, t, idx) if A.t[A.idxs[1]] == t - return A.u[:, A.idxs[1]], idx + return A.u[:, A.idxs[1]] end N = zero(A.u[:, 1]) D = zero(A.t[1]) @@ -112,13 +112,13 @@ function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, igu D += tmp @. N += (tmp * A.u[:, A.idxs[i]]) end - N / D, idx + N / D end function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess) idx = get_idx(A, t, iguess) wj = t - A.t[idx] - (@evalpoly wj A.u[idx] A.b[idx] A.c[idx] A.d[idx]), idx + @evalpoly wj A.u[idx] A.b[idx] A.c[idx] A.d[idx] end # ConstantInterpolation Interpolation @@ -130,7 +130,7 @@ function _interpolate(A::ConstantInterpolation{<:AbstractVector}, t::Number, igu # :right means that value to the right is used for interpolation idx = get_idx(A, t, iguess; side = :first, lb = 1, ub_shift = 0) end - A.u[idx], idx + A.u[idx] end function _interpolate(A::ConstantInterpolation{<:AbstractMatrix}, t::Number, iguess) @@ -141,7 +141,7 @@ function _interpolate(A::ConstantInterpolation{<:AbstractMatrix}, t::Number, igu # :right means that value to the right is used for interpolation idx = get_idx(A, t, iguess; side = :first, lb = 1, ub_shift = 0) end - A.u[:, idx], idx + A.u[:, idx] end # QuadraticSpline Interpolation @@ -150,7 +150,7 @@ function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) Cᵢ = A.u[idx] Δt = t - A.t[idx] σ = get_parameters(A, idx) - return A.z[idx] * Δt + σ * Δt^2 + Cᵢ, idx + return A.z[idx] * Δt + σ * Δt^2 + Cᵢ end # CubicSpline Interpolation @@ -162,15 +162,15 @@ function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess) c₁, c₂ = get_parameters(A, idx) C = c₁ * Δt₁ D = c₂ * Δt₂ - I + C + D, idx + I + C + D end # BSpline Curve Interpolation function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number, iguess) - t < A.t[1] && return A.u[1], 1 - t > A.t[end] && return A.u[end], lastindex(t) + t < A.t[1] && return A.u[1] + t > A.t[end] && return A.u[end] # change t into param [0 1] idx = get_idx(A, t, iguess) t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx]) @@ -181,13 +181,13 @@ function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}}, for i in nonzero_coefficient_idxs ucum += N[i] * A.c[i] end - ucum, idx + ucum end # BSpline Curve Approx function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, iguess) - t < A.t[1] && return A.u[1], 1 - t > A.t[end] && return A.u[end], lastindex(t) + t < A.t[1] && return A.u[1] + t > A.t[end] && return A.u[end] # change t into param [0 1] idx = get_idx(A, t, iguess) t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx]) @@ -197,7 +197,7 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i for i in nonzero_coefficient_idxs ucum += N[i] * A.c[i] end - ucum, idx + ucum end # Cubic Hermite Spline @@ -209,7 +209,7 @@ function _interpolate( out = A.u[idx] + Δt₀ * A.du[idx] c₁, c₂ = get_parameters(A, idx) out += Δt₀^2 * (c₁ + Δt₁ * c₂) - out, idx + out end # Quintic Hermite Spline @@ -221,5 +221,5 @@ function _interpolate( out = A.u[idx] + Δt₀ * (A.du[idx] + A.ddu[idx] * Δt₀ / 2) c₁, c₂, c₃ = get_parameters(A, idx) out += Δt₀^3 * (c₁ + Δt₁ * (c₂ + c₃ * Δt₁)) - out, idx + out end From f62439a9794926aeb3080eaa38312da224ecec3e Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 2 Aug 2024 18:12:03 +0200 Subject: [PATCH 5/7] Pass tests --- ext/DataInterpolationsChainRulesCoreExt.jl | 4 ++-- src/derivatives.jl | 12 ++++++------ src/integral_inverses.jl | 6 +++--- src/interpolation_methods.jl | 2 +- src/interpolation_utils.jl | 2 +- test/derivative_tests.jl | 2 +- 6 files changed, 14 insertions(+), 14 deletions(-) diff --git a/ext/DataInterpolationsChainRulesCoreExt.jl b/ext/DataInterpolationsChainRulesCoreExt.jl index 20c803a3..988e519b 100644 --- a/ext/DataInterpolationsChainRulesCoreExt.jl +++ b/ext/DataInterpolationsChainRulesCoreExt.jl @@ -52,7 +52,7 @@ end function u_tangent(A::LinearInterpolation, t, Δ) out = zero(A.u) - idx = get_idx(A, t, A.iguesser(t)) + idx = get_idx(A, t, A.iguesser) t_factor = (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) out[idx] = Δ * (one(eltype(out)) - t_factor) out[idx + 1] = Δ * t_factor @@ -61,7 +61,7 @@ end function u_tangent(A::QuadraticInterpolation, t, Δ) out = zero(A.u) - i₀, i₁, i₂ = _quad_interp_indices(A, t, A.iguesser(t)) + i₀, i₁, i₂ = _quad_interp_indices(A, t, A.iguesser) t₀ = A.t[i₀] t₁ = A.t[i₁] t₂ = A.t[i₂] diff --git a/src/derivatives.jl b/src/derivatives.jl index 940125a0..74b9f2e1 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -1,4 +1,4 @@ -function derivative(A, t, order = 1; iguess = A.iguesser(t)) +function derivative(A, t, order = 1; iguess = A.iguesser) ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) return if order == 1 @@ -15,7 +15,7 @@ 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, idx + slope end function _derivative(A::QuadraticInterpolation, t::Number, iguess) @@ -144,8 +144,8 @@ end function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number, iguess) # change t into param [0 1] - t < A.t[1] && return zero(A.u[1]), 1 - t > A.t[end] && return zero(A.u[end]), lastindex(t) + t < A.t[1] && return zero(A.u[1]) + t > A.t[end] && return zero(A.u[end]) idx = get_idx(A, t, iguess) n = length(A.t) scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx]) @@ -166,8 +166,8 @@ end # BSpline Curve Approx function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, iguess) # change t into param [0 1] - t < A.t[1] && return zero(A.u[1]), 1 - t > A.t[end] && return zero(A.u[end]), lastindex(t) + t < A.t[1] && return zero(A.u[1]) + t > A.t[end] && return zero(A.u[end]) idx = get_idx(A, t, iguess) scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx]) t_ = A.p[idx] + (t - A.t[idx]) * scale diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index d4d36e4f..4cddf7d3 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -18,7 +18,7 @@ invert_integral(A::AbstractInterpolation) = throw(IntegralInverseNotFoundError() _integral(A::AbstractIntegralInverseInterpolation, idx, t) = throw(IntegralNotFoundError()) function _derivative(A::AbstractIntegralInverseInterpolation, t::Number, iguess) - inv(A.itp(A(t))), A.iguesser(t) + inv(A.itp(A(t))) end """ @@ -64,7 +64,7 @@ function _interpolate( x = A.itp.u[idx] slope = get_parameters(A.itp, idx) u = A.u[idx] + 2Δt / (x + sqrt(x^2 + slope * 2Δt)) - u, idx + u end """ @@ -112,5 +112,5 @@ function _interpolate( # :right means that value to the right is used for interpolation idx_ = get_idx(A, t, idx; side = :first, lb = 1, ub_shift = 0) end - A.u[idx] + (t - A.t[idx]) / A.itp.u[idx_], idx + A.u[idx] + (t - A.t[idx]) / A.itp.u[idx_] end diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 8b3b223d..94436b82 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,4 +1,4 @@ -function _interpolate(A, t; iguess = A.iguesser(t)) +function _interpolate(A, t; iguess = A.iguesser) ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) return _interpolate(A, t, iguess) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index a54bdcb4..0df5e484 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -117,7 +117,7 @@ function looks_linear(t; threshold = 1e-2) norm_var < threshold^2 end -function get_idx(A::AbstractInterpolation, t, iguess; lb = 1, +function get_idx(A::AbstractInterpolation, t, iguess::Union{<:Integer, Guesser}; lb = 1, ub_shift = -1, idx_shift = 0, side = :last) tvec = A.t ub = length(tvec) + ub_shift diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index ea002ab5..a4eb349b 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -39,7 +39,7 @@ function test_derivatives(method; args = [], kwargs = [], name::String) @test isapprox(fdiff, adiff, atol = 1e-8) @test isapprox(fdiff2, adiff2, atol = 1e-8) # Cached index - if hasproperty(func, :iguesser) + if hasproperty(func, :iguesser) && !func.iguesser.linear_lookup @test abs(func.iguesser.idx_prev[] - searchsortedfirstcorrelated(func.t, _t, func.iguesser(_t))) <= 1 end From 3debd2909676a0c315ee41a76b90fe9ab146d6e8 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 2 Aug 2024 18:17:27 +0200 Subject: [PATCH 6/7] Update FindFirstFunctions compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 41990bcb..59d27ef9 100644 --- a/Project.toml +++ b/Project.toml @@ -27,7 +27,7 @@ DataInterpolationsSymbolicsExt = "Symbolics" Aqua = "0.8" BenchmarkTools = "1" ChainRulesCore = "1.24" -FindFirstFunctions = "1.1" +FindFirstFunctions = "1.3" FiniteDifferences = "0.12.31" ForwardDiff = "0.10.36" LinearAlgebra = "1.10" From 2e9b14ca3efc78b097936e99b86ccb5ff0c5ff45 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Sat, 3 Aug 2024 12:42:22 +0200 Subject: [PATCH 7/7] Remove user-supplied iguess option --- src/DataInterpolations.jl | 3 --- src/derivatives.jl | 3 ++- src/interpolation_methods.jl | 4 ++-- 3 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl index 783a796c..5f035e62 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -22,9 +22,6 @@ include("online.jl") include("show.jl") (interp::AbstractInterpolation)(t::Number) = _interpolate(interp, t) -function (interp::AbstractInterpolation)(t::Number, iguess::Integer) - _interpolate(interp, t; iguess) -end function (interp::AbstractInterpolation)(t::AbstractVector) u = get_u(interp.u, t) diff --git a/src/derivatives.jl b/src/derivatives.jl index 74b9f2e1..ff5b9850 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -1,5 +1,6 @@ -function derivative(A, t, order = 1; iguess = A.iguesser) +function derivative(A, t, order = 1) ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) + iguess = A.iguesser return if order == 1 _derivative(A, t, iguess) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 94436b82..697f7375 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,7 +1,7 @@ -function _interpolate(A, t; iguess = A.iguesser) +function _interpolate(A, t) ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) - return _interpolate(A, t, iguess) + return _interpolate(A, t, A.iguesser) end # Linear Interpolation