From 7e52d0d20678b82d850dfc0b7b405e8ecf4ff07c Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 3 Jul 2024 22:58:59 +0200 Subject: [PATCH 1/7] Integral inversion POC --- src/DataInterpolations.jl | 13 +++++++++++++ src/integral_inverses.jl | 37 +++++++++++++++++++++++++++++++++++++ src/integrals.jl | 4 ++-- src/interpolation_caches.jl | 12 ++++++++---- src/interpolation_utils.jl | 10 ++++++++++ 5 files changed, 70 insertions(+), 6 deletions(-) create mode 100644 src/integral_inverses.jl diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl index 7638f1a8..4ff881d9 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -17,6 +17,7 @@ include("interpolation_methods.jl") include("plot_rec.jl") include("derivatives.jl") include("integrals.jl") +include("integral_inverses.jl") include("online.jl") include("show.jl") @@ -74,6 +75,18 @@ function Base.showerror(io::IO, e::DerivativeNotFoundError) print(io, DERIVATIVE_NOT_FOUND_ERROR) end +const INTEGRAL_INVERSE_NOT_FOUND_ERROR = "Cannot invert the integral analytically. Please use Numerical methods." +struct IntegralInverseNotFoundError <: Exception end +function Base.showerror(io::IO, e::IntegralInverseNotFoundError) + print(io, INTEGRAL_INVERSE_NOT_FOUND_ERROR) +end + +const INTEGRAL_NOT_INVERTIBLE_ERROR = "The Interpolation is not positive everywhere so its integral is not invertible." +struct IntegralNotInvertibleError <: Exception end +function Base.showerror(io::IO, e::IntegralNotInvertibleError) + print(io, INTEGRAL_NOT_INVERTIBLE_ERROR) +end + export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline, BSplineInterpolation, BSplineApprox, CubicHermiteSpline, diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl new file mode 100644 index 00000000..daf5d598 --- /dev/null +++ b/src/integral_inverses.jl @@ -0,0 +1,37 @@ +abstract type AbstractIntegralInverseInterpolation{T} <: AbstractInterpolation{T} end + +struct LinearInterpolationIntInv{tType, IType, itpType, T} <: + AbstractIntegralInverseInterpolation{T} + u::tType + t::IType + extrapolate::Bool + idx_prev::Base.RefValue{Int} + itp::itpType + function LinearInterpolationIntInv(u, t, A) + new{typeof(u), typeof(t), typeof(A), eltype(u)}( + u, t, A.extrapolate, Ref(1), A) + end +end + +function invertible_integral(A::LinearInterpolation{<:AbstractVector{<:Number}}) + return all(A.u .> 0) +end + +invert_integral(A::AbstractInterpolation) = throw(IntegralInverseNotFoundError()) + +function invert_integral(A::LinearInterpolation{<:AbstractVector{<:Number}}) + !invertible_integral(A) && throw(IntegralNotInvertibleError()) + return LinearInterpolationIntInv(A.t, A.I, A) +end + +function _interpolate( + A::LinearInterpolationIntInv{<:AbstractVector{<:Number}}, t::Number, iguess) + (; itp) = A + idx = get_idx(A.t, t, iguess) + Δt = t - A.t[idx] + # TODO: Get from parameter cache: itp.p.slope[idx] + slope = (itp.u[idx + 1] - itp.u[idx]) / (itp.t[idx + 1] - itp.t[idx]) + x = itp.u[idx] + u = A.u[idx] + 2Δt / (x + sqrt(x^2 + 2slope * Δt)) + u, idx +end diff --git a/src/integrals.jl b/src/integrals.jl index 00036c0a..4703690c 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -7,9 +7,9 @@ function integral(A::AbstractInterpolation, t1::Number, t2::Number) ((t1 < A.t[1] || t1 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) ((t2 < A.t[1] || t2 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) # the index less than or equal to t1 - idx1 = max(1, min(searchsortedlast(A.t, t1), length(A.t) - 1)) + idx1 = get_idx(A.t, t1, 0) # the index less than t2 - idx2 = max(1, min(searchsortedlast(A.t, t2), length(A.t) - 1)) + idx2 = get_idx(A.t, t2, 0) if A.t[idx2] == t2 idx2 -= 1 end diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 15dba20e..b4c5afb7 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -13,19 +13,23 @@ Extrapolation extends the last linear polynomial on each side. - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. """ -struct LinearInterpolation{uType, tType, T} <: AbstractInterpolation{T} +struct LinearInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T} u::uType t::tType + I::IType extrapolate::Bool idx_prev::Base.RefValue{Int} - function LinearInterpolation(u, t, extrapolate) - new{typeof(u), typeof(t), eltype(u)}(u, t, extrapolate, Ref(1)) + function LinearInterpolation(u, t, I, extrapolate) + new{typeof(u), typeof(t), typeof(I), eltype(u)}(u, t, I, extrapolate, Ref(1)) end end function LinearInterpolation(u, t; extrapolate = false) u, t = munge_data(u, t) - LinearInterpolation(u, t, extrapolate) + # TODO: Bad workaround, this computes the cached parameters twice + A = LinearInterpolation(u, t, nothing, extrapolate) + I = cumulative_integral(A) + LinearInterpolation(u, t, I, extrapolate) end """ diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index bb62e78a..6eb00eab 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -75,3 +75,13 @@ function get_idx(tvec, t, iguess; lb = 1, ub_shift = -1, idx_shift = 0, side = : error("side must be :first or :last") end end + +function cumulative_integral(A) + if isempty(methods(_integral, (typeof(A), Any, Any))) + return nothing + end + integral_values = [_integral(A, idx, A.t[idx + 1]) - _integral(A, idx, A.t[idx]) + for idx in 1:(length(A.t) - 1)] + pushfirst!(integral_values, zero(first(integral_values))) + return cumsum(integral_values) +end From a62531e26e941e45844e9ec04e739007ad8ff94a Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Thu, 4 Jul 2024 19:37:41 +0200 Subject: [PATCH 2/7] Cache cumulative integral --- src/integrals.jl | 23 ++++----- src/interpolation_caches.jl | 98 +++++++++++++++++++++++-------------- 2 files changed, 74 insertions(+), 47 deletions(-) diff --git a/src/integrals.jl b/src/integrals.jl index 4703690c..0c3d1db0 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -9,17 +9,18 @@ function integral(A::AbstractInterpolation, t1::Number, t2::Number) # the index less than or equal to t1 idx1 = get_idx(A.t, t1, 0) # the index less than t2 - idx2 = get_idx(A.t, t2, 0) - if A.t[idx2] == t2 - idx2 -= 1 - end - total = zero(eltype(A.u)) - for idx in idx1:idx2 - lt1 = idx == idx1 ? t1 : A.t[idx] - lt2 = idx == idx2 ? t2 : A.t[idx + 1] - total += _integral(A, idx, lt2) - _integral(A, idx, lt1) + idx2 = get_idx(A.t, t2, 0; idx_shift = -1, side = :first) + + total = A.I[idx2] - A.I[idx1] + return if t1 == t2 + zero(total) + else + total += _integral(A, idx1, A.t[idx1]) + total -= _integral(A, idx1, t1) + total += _integral(A, idx2, t2) + total -= _integral(A, idx2, A.t[idx2]) + total end - total end function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}}, @@ -32,7 +33,7 @@ function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}}, t^2 * (u1 - u2) / (2 * t1 - 2 * t2) + t * (t1 * u2 - t2 * u1) / (t1 - t2) end -function _integral(A::ConstantInterpolation{<:AbstractVector}, idx::Number, t::Number) +function _integral(A::ConstantInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number) if A.dir === :left # :left means that value to the left is used for interpolation return A.u[idx] * t diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index b4c5afb7..3c313a23 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -26,7 +26,6 @@ end function LinearInterpolation(u, t; extrapolate = false) u, t = munge_data(u, t) - # TODO: Bad workaround, this computes the cached parameters twice A = LinearInterpolation(u, t, nothing, extrapolate) I = cumulative_integral(A) LinearInterpolation(u, t, I, extrapolate) @@ -48,22 +47,25 @@ Extrapolation extends the last quadratic polynomial on each side. - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. """ -struct QuadraticInterpolation{uType, tType, T} <: AbstractInterpolation{T} +struct QuadraticInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T} u::uType t::tType + I::IType mode::Symbol extrapolate::Bool idx_prev::Base.RefValue{Int} - function QuadraticInterpolation(u, t, mode, extrapolate) + function QuadraticInterpolation(u, t, I, mode, extrapolate) mode ∈ (:Forward, :Backward) || error("mode should be :Forward or :Backward for QuadraticInterpolation") - new{typeof(u), typeof(t), eltype(u)}(u, t, mode, extrapolate, Ref(1)) + new{typeof(u), typeof(t), typeof(I), eltype(u)}(u, t, I, mode, extrapolate, Ref(1)) end end function QuadraticInterpolation(u, t, mode; extrapolate = false) u, t = munge_data(u, t) - QuadraticInterpolation(u, t, mode, extrapolate) + A = QuadraticInterpolation(u, t, nothing, mode, extrapolate) + I = cumulative_integral(A) + QuadraticInterpolation(u, t, I, mode, extrapolate) end function QuadraticInterpolation(u, t; extrapolate = false) @@ -129,19 +131,21 @@ Extrapolation extends the last cubic polynomial on each side. - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. """ -struct AkimaInterpolation{uType, tType, bType, cType, dType, T} <: +struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <: AbstractInterpolation{T} u::uType t::tType + I::IType b::bType c::cType d::dType extrapolate::Bool idx_prev::Base.RefValue{Int} - function AkimaInterpolation(u, t, b, c, d, extrapolate) - new{typeof(u), typeof(t), typeof(b), typeof(c), + function AkimaInterpolation(u, t, I, b, c, d, extrapolate) + new{typeof(u), typeof(t), typeof(I), typeof(b), typeof(c), typeof(d), eltype(u)}(u, t, + I, b, c, d, @@ -173,7 +177,9 @@ function AkimaInterpolation(u, t; extrapolate = false) c = (3.0 .* m[3:(end - 2)] .- 2.0 .* b[1:(end - 1)] .- b[2:end]) ./ dt d = (b[1:(end - 1)] .+ b[2:end] .- 2.0 .* m[3:(end - 2)]) ./ dt .^ 2 - AkimaInterpolation(u, t, b, c, d, extrapolate) + A = AkimaInterpolation(u, t, nothing, b, c, d, extrapolate) + I = cumulative_integral(A) + AkimaInterpolation(u, t, I, b, c, d, extrapolate) end """ @@ -193,20 +199,23 @@ Extrapolation extends the last constant polynomial at the end points on each sid - `dir`: indicates which value should be used for interpolation (`:left` or `:right`). - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. """ -struct ConstantInterpolation{uType, tType, dirType, T} <: AbstractInterpolation{T} +struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T} u::uType t::tType + I::IType dir::Symbol # indicates if value to the $dir should be used for the interpolation extrapolate::Bool idx_prev::Base.RefValue{Int} - function ConstantInterpolation(u, t, dir, extrapolate) - new{typeof(u), typeof(t), typeof(dir), eltype(u)}(u, t, dir, extrapolate, Ref(1)) + function ConstantInterpolation(u, t, I, dir, extrapolate) + new{typeof(u), typeof(t), typeof(I), eltype(u)}(u, t, I, dir, extrapolate, Ref(1)) end end function ConstantInterpolation(u, t; dir = :left, extrapolate = false) u, t = munge_data(u, t) - ConstantInterpolation(u, t, dir, extrapolate) + A = ConstantInterpolation(u, t, nothing, dir, extrapolate) + I = cumulative_integral(A) + ConstantInterpolation(u, t, I, dir, extrapolate) end """ @@ -224,19 +233,21 @@ Extrapolation extends the last quadratic polynomial on each side. - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. """ -struct QuadraticSpline{uType, tType, tAType, dType, zType, T} <: +struct QuadraticSpline{uType, tType, IType, tAType, dType, zType, T} <: AbstractInterpolation{T} u::uType t::tType + I::IType tA::tAType d::dType z::zType extrapolate::Bool idx_prev::Base.RefValue{Int} - function QuadraticSpline(u, t, tA, d, z, extrapolate) - new{typeof(u), typeof(t), typeof(tA), + function QuadraticSpline(u, t, I, tA, d, z, extrapolate) + new{typeof(u), typeof(t), typeof(I), typeof(tA), typeof(d), typeof(z), eltype(u)}(u, t, + I, tA, d, z, @@ -246,9 +257,8 @@ struct QuadraticSpline{uType, tType, tAType, dType, zType, T} <: end end -function QuadraticSpline(u::uType, - t; - extrapolate = false) where {uType <: AbstractVector{<:Number}} +function QuadraticSpline( + u::uType, t; extrapolate = false) where {uType <: AbstractVector{<:Number}} u, t = munge_data(u, t) s = length(t) dl = ones(eltype(t), s - 1) @@ -261,7 +271,9 @@ function QuadraticSpline(u::uType, d = map(i -> i == 1 ? typed_zero : 2 // 1 * (u[i] - u[i - 1]) / (t[i] - t[i - 1]), 1:s) z = tA \ d - QuadraticSpline(u, t, tA, d, z, extrapolate) + A = QuadraticSpline(u, t, nothing, tA, d, z, extrapolate) + I = cumulative_integral(A) + QuadraticSpline(u, t, I, tA, d, z, extrapolate) end function QuadraticSpline(u::uType, t; extrapolate = false) where {uType <: AbstractVector} @@ -278,7 +290,9 @@ function QuadraticSpline(u::uType, t; extrapolate = false) where {uType <: Abstr d = transpose(reshape(reduce(hcat, d_), :, s)) z_ = reshape(transpose(tA \ d), size(u[1])..., :) z = [z_s for z_s in eachslice(z_, dims = ndims(z_))] - QuadraticSpline(u, t, tA, d, z, extrapolate) + A = QuadraticSpline(u, t, nothing, tA, d, z, extrapolate) + I = cumulative_integral(A) + QuadraticSpline(u, t, I, tA, d, z, extrapolate) end """ @@ -296,16 +310,18 @@ Second derivative on both ends are zero, which are also called "natural" boundar - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. """ -struct CubicSpline{uType, tType, hType, zType, T} <: AbstractInterpolation{T} +struct CubicSpline{uType, tType, IType, hType, zType, T} <: AbstractInterpolation{T} u::uType t::tType + I::IType h::hType z::zType extrapolate::Bool idx_prev::Base.RefValue{Int} - function CubicSpline(u, t, h, z, extrapolate) - new{typeof(u), typeof(t), typeof(h), typeof(z), eltype(u)}(u, + function CubicSpline(u, t, I, h, z, extrapolate) + new{typeof(u), typeof(t), typeof(I), typeof(h), typeof(z), eltype(u)}(u, t, + I, h, z, extrapolate, @@ -334,7 +350,9 @@ function CubicSpline(u::uType, 6(u[i + 1] - u[i]) / h[i + 1] - 6(u[i] - u[i - 1]) / h[i], 1:(n + 1)) z = tA \ d - CubicSpline(u, t, h[1:(n + 1)], z, extrapolate) + A = CubicSpline(u, t, nothing, h[1:(n + 1)], z, extrapolate) + I = cumulative_integral(A) + CubicSpline(u, t, I, h[1:(n + 1)], z, extrapolate) end function CubicSpline(u::uType, t; extrapolate = false) where {uType <: AbstractVector} @@ -352,7 +370,9 @@ function CubicSpline(u::uType, t; extrapolate = false) where {uType <: AbstractV d = transpose(reshape(reduce(hcat, d_), :, n + 1)) z_ = reshape(transpose(tA \ d), size(u[1])..., :) z = [z_s for z_s in eachslice(z_, dims = ndims(z_))] - CubicSpline(u, t, h[1:(n + 1)], z, extrapolate) + A = CubicSpline(u, t, nothing, h[1:(n + 1)], z, extrapolate) + I = cumulative_integral(A) + CubicSpline(u, t, I, h[1:(n + 1)], z, extrapolate) end """ @@ -631,16 +651,17 @@ It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomi - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. """ -struct CubicHermiteSpline{uType, tType, duType, pType, T} <: AbstractInterpolation{T} +struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInterpolation{T} du::uType u::uType t::tType + I::IType p::CubicHermiteParameterCache{pType} extrapolate::Bool idx_prev::Base.RefValue{Int} - function CubicHermiteSpline(du, u, t, p, extrapolate) - new{typeof(u), typeof(t), typeof(du), typeof(p.c₁), eltype(u)}( - du, u, t, p, extrapolate, Ref(1)) + function CubicHermiteSpline(du, u, t, I, p, extrapolate) + new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u)}( + du, u, t, I, p, extrapolate, Ref(1)) end end @@ -648,7 +669,9 @@ function CubicHermiteSpline(du, u, t; extrapolate = false) @assert length(u)==length(du) "Length of `u` is not equal to length of `du`." u, t = munge_data(u, t) p = CubicHermiteParameterCache(du, u, t) - return CubicHermiteSpline(du, u, t, p, extrapolate) + A = CubicHermiteSpline(du, u, t, nothing, p, extrapolate) + I = cumulative_integral(A) + CubicHermiteSpline(du, u, t, I, p, extrapolate) end """ @@ -667,18 +690,19 @@ It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polyno - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. """ -struct QuinticHermiteSpline{uType, tType, duType, dduType, pType, T} <: +struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <: AbstractInterpolation{T} ddu::uType du::uType u::uType t::tType + I::IType p::QuinticHermiteParameterCache{pType} extrapolate::Bool idx_prev::Base.RefValue{Int} - function QuinticHermiteSpline(ddu, du, u, t, p, extrapolate) - new{typeof(u), typeof(t), typeof(du), typeof(ddu), typeof(p.c₁), eltype(u)}( - ddu, du, u, t, p, extrapolate, Ref(1)) + function QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate) + 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)) end end @@ -686,5 +710,7 @@ function QuinticHermiteSpline(ddu, du, u, t; extrapolate = false) @assert length(u)==length(du)==length(ddu) "Length of `u` is not equal to length of `du` or `ddu`." u, t = munge_data(u, t) p = QuinticHermiteParameterCache(ddu, du, u, t) - return QuinticHermiteSpline(ddu, du, u, t, p, extrapolate) + A = QuinticHermiteSpline(ddu, du, u, t, nothing, p, extrapolate) + I = cumulative_integral(A) + QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate) end From 3036d8ceec5f6d38db5edd65474eedf0c35eb681 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Thu, 4 Jul 2024 20:15:20 +0200 Subject: [PATCH 3/7] Add tests --- src/integrals.jl | 3 ++- src/interpolation_caches.jl | 3 ++- test/integral_inverse_tests.jl | 23 +++++++++++++++++++++++ test/runtests.jl | 1 + 4 files changed, 28 insertions(+), 2 deletions(-) create mode 100644 test/integral_inverse_tests.jl diff --git a/src/integrals.jl b/src/integrals.jl index 0c3d1db0..8532cd71 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -33,7 +33,8 @@ function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}}, t^2 * (u1 - u2) / (2 * t1 - 2 * t2) + t * (t1 * u2 - t2 * u1) / (t1 - t2) end -function _integral(A::ConstantInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number) +function _integral( + A::ConstantInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number) if A.dir === :left # :left means that value to the left is used for interpolation return A.u[idx] * t diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 3c313a23..4258a102 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -701,7 +701,8 @@ struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <: extrapolate::Bool idx_prev::Base.RefValue{Int} function QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate) - new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(ddu), typeof(p.c₁), eltype(u)}( + 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)) end end diff --git a/test/integral_inverse_tests.jl b/test/integral_inverse_tests.jl new file mode 100644 index 00000000..b8a835f7 --- /dev/null +++ b/test/integral_inverse_tests.jl @@ -0,0 +1,23 @@ +using DataInterpolations +using DataInterpolations: integral, invert_integral + +function test_integral_inverses(method; args = [], kwargs = []) + A = method(args...; kwargs...) + @test hasfield(typeof(A), :I) + A_intinv = invert_integral(A) + @test A_intinv isa DataInterpolations.AbstractIntegralInverseInterpolation + ts = range(first(A.t), last(A.t), length = 100) + Is = integral.(Ref(A), ts) + ts_ = A_intinv.(Is) + @test ts ≈ ts_ +end + +@testset "Linear Interpolation" begin + t = collect(1:5) + u = [1.0, 1.0, 2.0, 4.0, 3.0] + test_integral_inverses(LinearInterpolation; args = [u, t]) + + u = [1.0, -1.0, 2.0, 4.0, 3.0] + A = LinearInterpolation(u, t) + @test_throws DataInterpolations.IntegralNotInvertibleError invert_integral(A) +end diff --git a/test/runtests.jl b/test/runtests.jl index 01f3a07e..0c722b2d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,7 @@ using SafeTestsets @safetestset "Interpolation Tests" include("interpolation_tests.jl") @safetestset "Derivative Tests" include("derivative_tests.jl") @safetestset "Integral Tests" include("integral_tests.jl") +@safetestset "Integral Inverse Tests" include("integral_inverse_tests.jl") @safetestset "Online Tests" include("online_tests.jl") @safetestset "Regularization Smoothing" include("regularization.jl") @safetestset "Show methods" include("show.jl") From e568de76df5aa7eac4852cf222b6add33b10793c Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 5 Jul 2024 09:53:43 +0200 Subject: [PATCH 4/7] Support inverting the integral of ConstantInterpolation --- src/integral_inverses.jl | 52 ++++++++++++++++++++++++++++------ test/integral_inverse_tests.jl | 16 +++++++++++ 2 files changed, 60 insertions(+), 8 deletions(-) diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index daf5d598..a6219b35 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -1,9 +1,12 @@ abstract type AbstractIntegralInverseInterpolation{T} <: AbstractInterpolation{T} end -struct LinearInterpolationIntInv{tType, IType, itpType, T} <: +_integral(A::AbstractIntegralInverseInterpolation, idx, t) = throw(IntegralNotFoundError()) +invert_integral(A::AbstractInterpolation) = throw(IntegralInverseNotFoundError()) + +struct LinearInterpolationIntInv{uType, tType, itpType, T} <: AbstractIntegralInverseInterpolation{T} - u::tType - t::IType + u::uType + t::tType extrapolate::Bool idx_prev::Base.RefValue{Int} itp::itpType @@ -17,8 +20,6 @@ function invertible_integral(A::LinearInterpolation{<:AbstractVector{<:Number}}) return all(A.u .> 0) end -invert_integral(A::AbstractInterpolation) = throw(IntegralInverseNotFoundError()) - function invert_integral(A::LinearInterpolation{<:AbstractVector{<:Number}}) !invertible_integral(A) && throw(IntegralNotInvertibleError()) return LinearInterpolationIntInv(A.t, A.I, A) @@ -26,12 +27,47 @@ end function _interpolate( A::LinearInterpolationIntInv{<:AbstractVector{<:Number}}, t::Number, iguess) - (; itp) = A idx = get_idx(A.t, t, iguess) Δt = t - A.t[idx] # TODO: Get from parameter cache: itp.p.slope[idx] - slope = (itp.u[idx + 1] - itp.u[idx]) / (itp.t[idx + 1] - itp.t[idx]) - x = itp.u[idx] + slope = (A.itp.u[idx + 1] - A.itp.u[idx]) / (A.itp.t[idx + 1] - A.itp.t[idx]) + x = A.itp.u[idx] u = A.u[idx] + 2Δt / (x + sqrt(x^2 + 2slope * Δt)) u, idx end + +struct ConstantInterpolationIntInv{uType, tType, itpType, T} <: + AbstractIntegralInverseInterpolation{T} + u::uType + t::tType + extrapolate::Bool + idx_prev::Base.RefValue{Int} + itp::itpType + function ConstantInterpolationIntInv(u, t, A) + new{typeof(u), typeof(t), typeof(A), eltype(u)}( + u, t, A.extrapolate, Ref(1), A + ) + end +end + +function invertible_integral(A::ConstantInterpolation{<:AbstractVector{<:Number}}) + return all(A.u .> 0) +end + +function invert_integral(A::ConstantInterpolation{<:AbstractVector{<:Number}}) + !invertible_integral(A) && throw(IntegralNotInvertibleError()) + return ConstantInterpolationIntInv(A.t, A.I, A) +end + +function _interpolate( + A::ConstantInterpolationIntInv{<:AbstractVector{<:Number}}, t::Number, iguess) + idx = get_idx(A.t, t, iguess; ub_shift = 0) + if A.itp.dir === :left + # :left means that value to the left is used for interpolation + idx_ = get_idx(A.t, t, idx; lb = 1, ub_shift = 0) + else + # :right means that value to the right is used for interpolation + idx_ = get_idx(A.t, t, idx; side = :first, lb = 1, ub_shift = 0) + end + A.u[idx] + (t - A.t[idx]) / A.itp.u[idx_], idx +end diff --git a/test/integral_inverse_tests.jl b/test/integral_inverse_tests.jl index b8a835f7..229dca6f 100644 --- a/test/integral_inverse_tests.jl +++ b/test/integral_inverse_tests.jl @@ -21,3 +21,19 @@ end A = LinearInterpolation(u, t) @test_throws DataInterpolations.IntegralNotInvertibleError invert_integral(A) end + +@testset "Constant Interpolation" begin + t = collect(1:5) + u = [1.0, 1.0, 2.0, 4.0, 3.0] + test_integral_inverses(ConstantInterpolation; args = [u, t]) + test_integral_inverses(ConstantInterpolation; args = [u, t], kwargs = [:dir => :right]) + + u = [1.0, -1.0, 2.0, 4.0, 3.0] + A = ConstantInterpolation(u, t) + @test_throws DataInterpolations.IntegralNotInvertibleError invert_integral(A) +end + +t = collect(1:5) +u = [1.0, 1.0, 2.0, 4.0, 3.0] +A = QuadraticInterpolation(u, t) +@test_throws DataInterpolations.IntegralInverseNotFoundError invert_integral(A) From 6da3a56e3b652c5a523ccc83ea782d9b20e89add Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 5 Jul 2024 18:58:25 +0200 Subject: [PATCH 5/7] Add docs page --- docs/make.jl | 2 +- docs/src/inverting_integrals.md | 53 +++++++++++++++++++++++++++++++++ src/DataInterpolations.jl | 2 +- src/integral_inverses.jl | 18 +++++++++++ 4 files changed, 73 insertions(+), 2 deletions(-) create mode 100644 docs/src/inverting_integrals.md diff --git a/docs/make.jl b/docs/make.jl index 8a2835d7..94546761 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -13,6 +13,6 @@ makedocs(modules = [DataInterpolations], format = Documenter.HTML(assets = ["assets/favicon.ico"], canonical = "https://docs.sciml.ai/DataInterpolations/stable/"), pages = ["index.md", "Methods" => "methods.md", - "Interface" => "interface.md", "Manual" => "manual.md"]) + "Interface" => "interface.md", "Manual" => "manual.md", "Inverting Integrals" => "inverting_integrals.md"]) deploydocs(repo = "github.com/SciML/DataInterpolations.jl"; push_preview = true) diff --git a/docs/src/inverting_integrals.md b/docs/src/inverting_integrals.md new file mode 100644 index 00000000..2d10f848 --- /dev/null +++ b/docs/src/inverting_integrals.md @@ -0,0 +1,53 @@ +# Inverting integrals + +Solving implicit integral problems of the form: + +```math +\begin{equation} + \text{find $t$ such that } \int_{t_1}^t f(\tau)\text{d}\tau = V \ge 0 +\end{equation} +``` + +is supported for interpolations $f$ that are strictly positive and of one of these types: + + - `ConstantInterpolation` + - `LinearInterpolation` + +This is done by creating an 'integral inverse' interpolation object which can efficiently compute $t$ for a given value of $V$, see the example below. + +```@example inverting_integrals +using Random #hide +Random.seed!(1234) # hide +using DataInterpolations +using Plots + +# Create LinearInterpolation object from the +u = sqrt.(1:25) + (2.0 * rand(25) .- 1.0) / 3 +t = cumsum(rand(25)) +A = LinearInterpolation(u, t) + +# Create LinearInterpolationIntInv object +# from the LinearInterpolation object +A_intinv = DataInterpolations.invert_integral(A) + +# Get the t values up to and including the +# solution to the integral problem +V = 25.0 +t_ = A_intinv(V) +ts = t[t .<= t_] +push!(ts, t_) + +# Plot results +plot(A; label = "Linear Interpolation") +plot!(ts, A.(ts), fillrange = 0.0, fillalpha = 0.75, + fc = :blues, lw = 0, label = "Area of $V") +``` + +## Docstrings + +```@docs +DataInterpolations.invert_integral +ConstantInterpolationIntInv +LinearInterpolationIntInv +``` + diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl index 4ff881d9..d5b69e77 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -90,7 +90,7 @@ end export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline, BSplineInterpolation, BSplineApprox, CubicHermiteSpline, - QuinticHermiteSpline + QuinticHermiteSpline, LinearInterpolationIntInv, ConstantInterpolationIntInv # added for RegularizationSmooth, JJS 11/27/21 ### Regularization data smoothing and interpolation diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index a6219b35..e8ecface 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -1,8 +1,23 @@ abstract type AbstractIntegralInverseInterpolation{T} <: AbstractInterpolation{T} end _integral(A::AbstractIntegralInverseInterpolation, idx, t) = throw(IntegralNotFoundError()) + +""" + invert_integral(A::AbstractInterpolation)::AbstractIntegralInverseInterpolation + +Creates the inverted integral interpolation object from the given interpolation. Conditions: + - The range of `A` must be strictly positive + - There must be an ordering defined on the data type of `A.u` + - This is currently only supported for ConstantInterpolation and LinearInterpolation + +## Arguments + - `A`: interpolation object satisfying the above requirements +""" invert_integral(A::AbstractInterpolation) = throw(IntegralInverseNotFoundError()) +""" + some stuff +""" struct LinearInterpolationIntInv{uType, tType, itpType, T} <: AbstractIntegralInverseInterpolation{T} u::uType @@ -36,6 +51,9 @@ function _interpolate( u, idx end +""" + some stuff +""" struct ConstantInterpolationIntInv{uType, tType, itpType, T} <: AbstractIntegralInverseInterpolation{T} u::uType From f36cd1827672eea6590ff035a11922afa1de36ac Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Sat, 6 Jul 2024 11:05:06 +0200 Subject: [PATCH 6/7] Add generic derivative --- docs/src/inverting_integrals.md | 1 - src/integral_inverses.jl | 10 ++++++++-- test/integral_inverse_tests.jl | 11 +++++++++-- 3 files changed, 17 insertions(+), 5 deletions(-) diff --git a/docs/src/inverting_integrals.md b/docs/src/inverting_integrals.md index 2d10f848..8cbfcc3e 100644 --- a/docs/src/inverting_integrals.md +++ b/docs/src/inverting_integrals.md @@ -50,4 +50,3 @@ DataInterpolations.invert_integral ConstantInterpolationIntInv LinearInterpolationIntInv ``` - diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index e8ecface..f9aecc2c 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -1,20 +1,26 @@ abstract type AbstractIntegralInverseInterpolation{T} <: AbstractInterpolation{T} end -_integral(A::AbstractIntegralInverseInterpolation, idx, t) = throw(IntegralNotFoundError()) - """ invert_integral(A::AbstractInterpolation)::AbstractIntegralInverseInterpolation Creates the inverted integral interpolation object from the given interpolation. Conditions: + - The range of `A` must be strictly positive - There must be an ordering defined on the data type of `A.u` - This is currently only supported for ConstantInterpolation and LinearInterpolation ## Arguments + - `A`: interpolation object satisfying the above requirements """ 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[] +end + """ some stuff """ diff --git a/test/integral_inverse_tests.jl b/test/integral_inverse_tests.jl index 229dca6f..fd83be67 100644 --- a/test/integral_inverse_tests.jl +++ b/test/integral_inverse_tests.jl @@ -1,8 +1,9 @@ using DataInterpolations -using DataInterpolations: integral, invert_integral +using DataInterpolations: integral, derivative, invert_integral +using FiniteDifferences function test_integral_inverses(method; args = [], kwargs = []) - A = method(args...; kwargs...) + A = method(args...; kwargs..., extrapolate = true) @test hasfield(typeof(A), :I) A_intinv = invert_integral(A) @test A_intinv isa DataInterpolations.AbstractIntegralInverseInterpolation @@ -10,6 +11,12 @@ function test_integral_inverses(method; args = [], kwargs = []) Is = integral.(Ref(A), ts) ts_ = A_intinv.(Is) @test ts ≈ ts_ + + for I in Is + cdiff = forward_fdm(5, 1; geom = true)(A_intinv, I) + adiff = derivative(A_intinv, I) + @test cdiff ≈ adiff + end end @testset "Linear Interpolation" begin From bf6996c12dc6944c0af2882ed1adb4c8608943d5 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Sat, 6 Jul 2024 16:10:12 +0200 Subject: [PATCH 7/7] Merge remote-tracking branch 'upstream/master' into invert_integrals --- Project.toml | 2 + docs/src/interface.md | 19 +- ext/DataInterpolationsOptimExt.jl | 5 +- ...ataInterpolationsRegularizationToolsExt.jl | 29 ++- joss/paper.md | 5 + src/DataInterpolations.jl | 20 +- src/derivatives.jl | 48 ++-- src/integral_inverses.jl | 4 +- src/integrals.jl | 84 ++---- src/interpolation_caches.jl | 243 +++++++++++------- src/interpolation_methods.jl | 136 +++++----- src/interpolation_utils.jl | 95 +++++-- src/online.jl | 77 ++++-- src/parameter_caches.jl | 92 ++++++- test/derivative_tests.jl | 7 + test/integral_tests.jl | 2 +- test/interpolation_tests.jl | 41 ++- test/online_tests.jl | 39 ++- test/parameter_tests.jl | 31 +++ 19 files changed, 630 insertions(+), 349 deletions(-) diff --git a/Project.toml b/Project.toml index 544f1654..b06d2adb 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" +ReadOnlyArrays = "988b38a3-91fc-5605-94a2-ee2116b3bd83" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" @@ -32,6 +33,7 @@ LinearAlgebra = "1.10" Optim = "1.6" PrettyTables = "2" QuadGK = "2.9.1" +ReadOnlyArrays = "0.2.0" RecipesBase = "1.3" Reexport = "1" RegularizationTools = "0.6" diff --git a/docs/src/interface.md b/docs/src/interface.md index 3e7df7e1..ca5e9819 100644 --- a/docs/src/interface.md +++ b/docs/src/interface.md @@ -17,7 +17,7 @@ t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3] All interpolation methods return an object from which we can compute the value of the dependent variable at any time point. -We will use the `CubicSpline` method for demonstration, but the API is the same for all the methods. We can also pass the `extrapolate=true` keyword if we want to allow the interpolation to go beyond the range of the timepoints. The default value is `extrapolate=false`. +We will use the `CubicSpline` method for demonstration, but the API is the same for all the methods. We can also pass the `extrapolate = true` keyword if we want to allow the interpolation to go beyond the range of the timepoints. The default value is `extrapolate = false`. ```@example interface A1 = CubicSpline(u, t) @@ -35,6 +35,23 @@ A2(300.0) The values computed beyond the range of the time points provided during interpolation will not be reliable, as these methods only perform well within the range and the first/last piece polynomial fit is extrapolated on either side which might not reflect the true nature of the data. +The keyword `safetycopy = false` can be passed to make sure no copies of `u` and `t` are made when initializing the interpolation object. + +```@example interface +A3 = QuadraticInterpolation(u, t; safetycopy = false) + +# Check for same memory +u === A3.u.parent +``` + +Note that this does not prevent allocation in every interpolation constructor call, because parameter values are cached for all interpolation types except [`ConstantInterpolation`](@ref). + +Because of the caching of parameters which depend on `u` and `t`, this data should not be mutated. Therefore `u` and `t` are wrapped in a `ReadOnlyArray` from [ReadOnlyArrays.jl](https://github.com/JuliaArrays/ReadOnlyArrays.jl). + +```@repl interface +A3.t[2] = 3.14 +``` + ## Derivatives Derivatives of the interpolated curves can also be computed at any point for all the methods. Derivatives upto second order is supported where first order derivative is computed analytically and second order using `ForwardDiff.jl`. Order is passed as the third argument. It is 1 by default. diff --git a/ext/DataInterpolationsOptimExt.jl b/ext/DataInterpolationsOptimExt.jl index b3bce295..5528503f 100644 --- a/ext/DataInterpolationsOptimExt.jl +++ b/ext/DataInterpolationsOptimExt.jl @@ -18,8 +18,9 @@ function Curvefit(u, box = false, lb = nothing, ub = nothing; - extrapolate = false) - u, t = munge_data(u, t) + extrapolate = false, + safetycopy = false) + u, t = munge_data(u, t, safetycopy) errfun(t, u, p) = sum(abs2.(u .- model(t, p))) if box == false mfit = optimize(p -> errfun(t, u, p), p0, alg) diff --git a/ext/DataInterpolationsRegularizationToolsExt.jl b/ext/DataInterpolationsRegularizationToolsExt.jl index a7a914e3..10ea3e4c 100644 --- a/ext/DataInterpolationsRegularizationToolsExt.jl +++ b/ext/DataInterpolationsRegularizationToolsExt.jl @@ -69,8 +69,8 @@ A = RegularizationSmooth(u, t, t̂, wls, wr, d; λ = 1.0, alg = :gcv_svd) """ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector, wls::AbstractVector, wr::AbstractVector, d::Int = 2; - λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false) - u, t = munge_data(u, t) + λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true) + u, t = munge_data(u, t, safetycopy) M = _mapping_matrix(t̂, t) Wls½ = LA.diagm(sqrt.(wls)) Wr½ = LA.diagm(sqrt.(wr)) @@ -86,8 +86,8 @@ A = RegularizationSmooth(u, t, d; λ = 1.0, alg = :gcv_svd, extrapolate = false) """ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, d::Int = 2; λ::Real = 1.0, - alg::Symbol = :gcv_svd, extrapolate::Bool = false) - u, t = munge_data(u, t) + alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true) + u, t = munge_data(u, t, safetycopy) t̂ = t N = length(t) M = Array{Float64}(LA.I, N, N) @@ -114,8 +114,9 @@ A = RegularizationSmooth(u, t, t̂, d; λ = 1.0, alg = :gcv_svd, extrapolate = f ``` """ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector, - d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false) - u, t = munge_data(u, t) + d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd, + extrapolate::Bool = false, safetycopy::Bool = true) + u, t = munge_data(u, t, safetycopy) N, N̂ = length(t), length(t̂) M = _mapping_matrix(t̂, t) Wls½ = Array{Float64}(LA.I, N, N) @@ -142,8 +143,8 @@ A = RegularizationSmooth(u, t, t̂, wls, d; λ = 1.0, alg = :gcv_svd, extrapolat """ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector, wls::AbstractVector, d::Int = 2; λ::Real = 1.0, - alg::Symbol = :gcv_svd, extrapolate::Bool = false) - u, t = munge_data(u, t) + alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true) + u, t = munge_data(u, t, safetycopy) N, N̂ = length(t), length(t̂) M = _mapping_matrix(t̂, t) Wls½ = LA.diagm(sqrt.(wls)) @@ -171,8 +172,8 @@ A = RegularizationSmooth( """ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing, wls::AbstractVector, d::Int = 2; λ::Real = 1.0, - alg::Symbol = :gcv_svd, extrapolate::Bool = false) - u, t = munge_data(u, t) + alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true) + u, t = munge_data(u, t, safetycopy) t̂ = t N = length(t) M = Array{Float64}(LA.I, N, N) @@ -201,8 +202,8 @@ A = RegularizationSmooth( """ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing, wls::AbstractVector, wr::AbstractVector, d::Int = 2; - λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false) - u, t = munge_data(u, t) + λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true) + u, t = munge_data(u, t, safetycopy) t̂ = t N = length(t) M = Array{Float64}(LA.I, N, N) @@ -231,8 +232,8 @@ A = RegularizationSmooth( """ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing, wls::Symbol, d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd, - extrapolate::Bool = false) - u, t = munge_data(u, t) + extrapolate::Bool = false, safetycopy::Bool = true) + u, t = munge_data(u, t, safetycopy) t̂ = t N = length(t) M = Array{Float64}(LA.I, N, N) diff --git a/joss/paper.md b/joss/paper.md index d53a1d70..22530a50 100644 --- a/joss/paper.md +++ b/joss/paper.md @@ -13,6 +13,9 @@ authors: affiliation: "1, 2, 3" - name: Shubham Maddhashiya affiliation: 3 + - name: Bart de Koning + orcid: 0009-0005-6134-6608 + affiliation: 4 affiliations: - name: JuliaHub index: 1 @@ -20,6 +23,8 @@ affiliations: index: 2 - name: Pumas-AI index: 3 + - name: Deltares + index: 4 date: 6 June 2024 bibliography: paper.bib --- diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl index d5b69e77..c86a6579 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -7,6 +7,7 @@ abstract type AbstractInterpolation{T} end using LinearAlgebra, RecipesBase using PrettyTables using ForwardDiff +using ReadOnlyArrays import FindFirstFunctions: searchsortedfirstcorrelated, searchsortedlastcorrelated, bracketstrictlymontonic @@ -87,6 +88,12 @@ function Base.showerror(io::IO, e::IntegralNotInvertibleError) print(io, INTEGRAL_NOT_INVERTIBLE_ERROR) end +const MUST_COPY_ERROR = "A copy must be made of u, t to filter missing data" +struct MustCopyError <: Exception end +function Base.showerror(io::IO, e::MustCopyError) + print(io, MUST_COPY_ERROR) +end + export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline, BSplineInterpolation, BSplineApprox, CubicHermiteSpline, @@ -118,12 +125,13 @@ struct RegularizationSmooth{uType, tType, T, T2, ITP <: AbstractInterpolation{T} alg, Aitp, extrapolate) - new{typeof(u), typeof(t), eltype(u), typeof(λ), typeof(Aitp)}(u, - û, - t, - t̂, - wls, - wr, + new{typeof(u), typeof(t), eltype(u), typeof(λ), typeof(Aitp)}( + readonly_wrap(u), + readonly_wrap(û), + readonly_wrap(t), + readonly_wrap(t̂), + readonly_wrap(oftype(u.parent, wls)), + readonly_wrap(oftype(u.parent, wr)), d, λ, alg, diff --git a/src/derivatives.jl b/src/derivatives.jl index bb3c938e..30c76fd0 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -17,30 +17,17 @@ function derivative(A, t, order = 1) end end -function _derivative(A::LinearInterpolation{<:AbstractVector}, t::Number, iguess) +function _derivative(A::LinearInterpolation, t::Number, iguess) idx = get_idx(A.t, t, iguess; idx_shift = -1, ub_shift = -2, side = :first) - (A.u[idx + 1] - A.u[idx]) / (A.t[idx + 1] - A.t[idx]), idx + A.p.slope[idx], idx end -function _derivative(A::LinearInterpolation{<:AbstractMatrix}, t::Number, iguess) - idx = get_idx(A.t, t, iguess; idx_shift = -1, ub_shift = -2, side = :first) - (@views @. (A.u[:, idx + 1] - A.u[:, idx]) / (A.t[idx + 1] - A.t[idx])), idx -end - -function _derivative(A::QuadraticInterpolation{<:AbstractVector}, t::Number, iguess) - i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess) - dl₀ = (2t - A.t[i₁] - A.t[i₂]) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂])) - dl₁ = (2t - A.t[i₀] - A.t[i₂]) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂])) - dl₂ = (2t - A.t[i₀] - A.t[i₁]) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁])) - A.u[i₀] * dl₀ + A.u[i₁] * dl₁ + A.u[i₂] * dl₂, i₀ -end - -function _derivative(A::QuadraticInterpolation{<:AbstractMatrix}, t::Number, iguess) +function _derivative(A::QuadraticInterpolation, t::Number, iguess) i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess) - dl₀ = (2t - A.t[i₁] - A.t[i₂]) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂])) - dl₁ = (2t - A.t[i₀] - A.t[i₂]) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂])) - dl₂ = (2t - A.t[i₀] - A.t[i₁]) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁])) - (@views @. A.u[:, i₀] * dl₀ + A.u[:, i₁] * dl₁ + A.u[:, i₂] * dl₂), i₀ + du₀ = A.p.l₀[i₀] * (2t - A.t[i₁] - A.t[i₂]) + du₁ = A.p.l₁[i₀] * (2t - A.t[i₀] - A.t[i₂]) + du₂ = A.p.l₂[i₀] * (2t - A.t[i₀] - A.t[i₁]) + return @views @. du₀ + du₁ + du₂, i₀ end function _derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number) @@ -125,6 +112,10 @@ function _derivative(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess) (@evalpoly wj A.b[idx] 2A.c[j] 3A.d[j]), idx end +function _derivative(A::ConstantInterpolation, t::Number, iguess) + return zero(first(A.u)), iguess +end + function _derivative(A::ConstantInterpolation{<:AbstractVector}, t::Number) ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) return isempty(searchsorted(A.t, t)) ? zero(A.u[1]) : eltype(A.u)(NaN) @@ -138,17 +129,18 @@ end # QuadraticSpline Interpolation function _derivative(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) idx = get_idx(A.t, t, iguess; lb = 2, ub_shift = 0, side = :first) - σ = 1 // 2 * (A.z[idx] - A.z[idx - 1]) / (A.t[idx] - A.t[idx - 1]) + σ = A.p.σ[idx - 1] A.z[idx - 1] + 2σ * (t - A.t[idx - 1]), idx end # CubicSpline Interpolation function _derivative(A::CubicSpline{<:AbstractVector}, t::Number, iguess) idx = get_idx(A.t, t, iguess) - dI = -3A.z[idx] * (A.t[idx + 1] - t)^2 / (6A.h[idx + 1]) + - 3A.z[idx + 1] * (t - A.t[idx])^2 / (6A.h[idx + 1]) - dC = A.u[idx + 1] / A.h[idx + 1] - A.z[idx + 1] * A.h[idx + 1] / 6 - dD = -(A.u[idx] / A.h[idx + 1] - A.z[idx] * A.h[idx + 1] / 6) + Δt₁ = t - A.t[idx] + Δt₂ = A.t[idx + 1] - t + dI = (-A.z[idx] * Δt₂^2 + A.z[idx + 1] * Δt₁^2) / (2A.h[idx + 1]) + dC = A.p.c₁[idx] + dD = -A.p.c₂[idx] dI + dC + dD, idx end @@ -160,7 +152,8 @@ function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Num n = length(A.t) scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx]) t_ = A.p[idx] + (t - A.t[idx]) * scale - N = DataInterpolations.spline_coefficients(n, A.d - 1, A.k, t_) + N = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.N + spline_coefficients!(N, A.d - 1, A.k, t_) ducum = zero(eltype(A.u)) if t == A.t[1] ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2]) @@ -180,7 +173,8 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig idx = get_idx(A.t, 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 - N = spline_coefficients(A.h, A.d - 1, A.k, t_) + N = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.N + spline_coefficients!(N, A.d - 1, A.k, t_) ducum = zero(eltype(A.u)) if t == A.t[1] ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2]) diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index f9aecc2c..1e2ca581 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -50,10 +50,8 @@ function _interpolate( A::LinearInterpolationIntInv{<:AbstractVector{<:Number}}, t::Number, iguess) idx = get_idx(A.t, t, iguess) Δt = t - A.t[idx] - # TODO: Get from parameter cache: itp.p.slope[idx] - slope = (A.itp.u[idx + 1] - A.itp.u[idx]) / (A.itp.t[idx + 1] - A.itp.t[idx]) x = A.itp.u[idx] - u = A.u[idx] + 2Δt / (x + sqrt(x^2 + 2slope * Δt)) + u = A.u[idx] + 2Δt / (x + sqrt(x^2 + A.itp.p.slope[idx] * 2Δt)) u, idx end diff --git a/src/integrals.jl b/src/integrals.jl index 8532cd71..3040189f 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -6,6 +6,7 @@ end function integral(A::AbstractInterpolation, t1::Number, t2::Number) ((t1 < A.t[1] || t1 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) ((t2 < A.t[1] || t2 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) + !hasfield(typeof(A), :I) && throw(IntegralNotFoundError()) # the index less than or equal to t1 idx1 = get_idx(A.t, t1, 0) # the index less than t2 @@ -26,11 +27,8 @@ end function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number) - t1 = A.t[idx] - t2 = A.t[idx + 1] - u1 = A.u[idx] - u2 = A.u[idx + 1] - t^2 * (u1 - u2) / (2 * t1 - 2 * t2) + t * (t1 * u2 - t2 * u1) / (t1 - t2) + Δt = t - A.t[idx] + Δt * (A.u[idx] + A.p.slope[idx] * Δt / 2) end function _integral( @@ -49,49 +47,30 @@ function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}}, t::Number) A.mode == :Backward && idx > 1 && (idx -= 1) idx = min(length(A.t) - 2, idx) - t1 = A.t[idx] - t2 = A.t[idx + 1] - t3 = A.t[idx + 2] - u1 = A.u[idx] - u2 = A.u[idx + 1] - u3 = A.u[idx + 2] - (t^3 * (-t1 * u2 + t1 * u3 + t2 * u1 - t2 * u3 - t3 * u1 + t3 * u2) / - (3 * t1^2 * t2 - 3 * t1^2 * t3 - 3 * t1 * t2^2 + 3 * t1 * t3^2 + 3 * t2^2 * t3 - - 3 * t2 * t3^2) + - t^2 * (t1^2 * u2 - t1^2 * u3 - t2^2 * u1 + t2^2 * u3 + t3^2 * u1 - t3^2 * u2) / - (2 * t1^2 * t2 - 2 * t1^2 * t3 - 2 * t1 * t2^2 + 2 * t1 * t3^2 + 2 * t2^2 * t3 - - 2 * t2 * t3^2) + - t * - (t1^2 * t2 * u3 - t1^2 * t3 * u2 - t1 * t2^2 * u3 + t1 * t3^2 * u2 + t2^2 * t3 * u1 - - t2 * t3^2 * u1) / - (t1^2 * t2 - t1^2 * t3 - t1 * t2^2 + t1 * t3^2 + t2^2 * t3 - t2 * t3^2)) + t₀ = A.t[idx] + t₁ = A.t[idx + 1] + t₂ = A.t[idx + 2] + + t_sq = (t^2) / 3 + Iu₀ = A.p.l₀[idx] * t * (t_sq - t * (t₁ + t₂) / 2 + t₁ * t₂) + Iu₁ = A.p.l₁[idx] * t * (t_sq - t * (t₀ + t₂) / 2 + t₀ * t₂) + Iu₂ = A.p.l₂[idx] * t * (t_sq - t * (t₀ + t₁) / 2 + t₀ * t₁) + return Iu₀ + Iu₁ + Iu₂ end function _integral(A::QuadraticSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number) - t1 = A.t[idx] - t2 = A.t[idx + 1] - u1 = A.u[idx] - z1 = A.z[idx] - z2 = A.z[idx + 1] - t^3 * (z1 - z2) / (6 * t1 - 6 * t2) + t^2 * (t1 * z2 - t2 * z1) / (2 * t1 - 2 * t2) + - t * (-t1^2 * z1 - t1^2 * z2 + 2 * t1 * t2 * z1 + 2 * t1 * u1 - 2 * t2 * u1) / - (2 * t1 - 2 * t2) + Cᵢ = A.u[idx] + Δt = t - A.t[idx] + return A.z[idx] * Δt^2 / 2 + A.p.σ[idx] * Δt^3 / 3 + Cᵢ * Δt end function _integral(A::CubicSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number) - t1 = A.t[idx] - t2 = A.t[idx + 1] - u1 = A.u[idx] - u2 = A.u[idx + 1] - z1 = A.z[idx] - z2 = A.z[idx + 1] - h2 = A.h[idx + 1] - (t^4 * (-z1 + z2) / (24 * h2) + t^3 * (-t1 * z2 + t2 * z1) / (6 * h2) + - t^2 * (h2^2 * z1 - h2^2 * z2 + 3 * t1^2 * z2 - 3 * t2^2 * z1 - 6 * u1 + 6 * u2) / - (12 * h2) + - t * - (h2^2 * t1 * z2 - h2^2 * t2 * z1 - t1^3 * z2 - 6 * t1 * u2 + t2^3 * z1 + 6 * t2 * u1) / - (6 * h2)) + Δt₁sq = (t - A.t[idx])^2 / 2 + Δt₂sq = (A.t[idx + 1] - t)^2 / 2 + II = (-A.z[idx] * Δt₂sq^2 + A.z[idx + 1] * Δt₁sq^2) / (6A.h[idx + 1]) + IC = A.p.c₁[idx] * Δt₁sq + ID = -A.p.c₂[idx] * Δt₂sq + II + IC + ID end function _integral(A::AkimaInterpolation{<:AbstractVector{<:Number}}, @@ -102,24 +81,9 @@ function _integral(A::AkimaInterpolation{<:AbstractVector{<:Number}}, A.d[idx] * ((t - t1)^4 / 4) end -integral(A::LagrangeInterpolation, t1::Number, t2::Number) = throw(IntegralNotFoundError()) -integral(A::LagrangeInterpolation, t::Number) = throw(IntegralNotFoundError()) - -function integral(A::BSplineInterpolation{<:AbstractVector{<:Number}}, - t1::Number, - t2::Number) - throw(IntegralNotFoundError()) -end -function integral(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number) - throw(IntegralNotFoundError()) -end - -function integral(A::BSplineApprox{<:AbstractVector{<:Number}}, t1::Number, t2::Number) - throw(IntegralNotFoundError()) -end -function integral(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number) - throw(IntegralNotFoundError()) -end +_integral(A::LagrangeInterpolation, idx::Number, t::Number) = throw(IntegralNotFoundError()) +_integral(A::BSplineInterpolation, idx::Number, t::Number) = throw(IntegralNotFoundError()) +_integral(A::BSplineApprox, idx::Number, t::Number) = throw(IntegralNotFoundError()) # Cubic Hermite Spline function _integral( diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 4258a102..83e04fe5 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -12,23 +12,28 @@ Extrapolation extends the last linear polynomial on each side. ## Keyword Arguments - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`. """ -struct LinearInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T} +struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T} u::uType t::tType I::IType + p::LinearParameterCache{pType} extrapolate::Bool idx_prev::Base.RefValue{Int} - function LinearInterpolation(u, t, I, extrapolate) - new{typeof(u), typeof(t), typeof(I), eltype(u)}(u, t, I, extrapolate, Ref(1)) + safetycopy::Bool + function LinearInterpolation(u, t, I, p, extrapolate, safetycopy) + new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}( + u, t, I, p, extrapolate, Ref(1), safetycopy) end end -function LinearInterpolation(u, t; extrapolate = false) - u, t = munge_data(u, t) - A = LinearInterpolation(u, t, nothing, extrapolate) +function LinearInterpolation(u, t; extrapolate = false, safetycopy = true) + u, t = munge_data(u, t, safetycopy) + p = LinearParameterCache(u, t) + A = LinearInterpolation(u, t, nothing, p, extrapolate, safetycopy) I = cumulative_integral(A) - LinearInterpolation(u, t, I, extrapolate) + LinearInterpolation(u, t, I, p, extrapolate, safetycopy) end """ @@ -46,30 +51,35 @@ Extrapolation extends the last quadratic polynomial on each side. ## Keyword Arguments - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`. """ -struct QuadraticInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T} +struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T} u::uType t::tType I::IType + p::QuadraticParameterCache{pType} mode::Symbol extrapolate::Bool idx_prev::Base.RefValue{Int} - function QuadraticInterpolation(u, t, I, mode, extrapolate) + safetycopy::Bool + function QuadraticInterpolation(u, t, I, p, mode, extrapolate, safetycopy) mode ∈ (:Forward, :Backward) || error("mode should be :Forward or :Backward for QuadraticInterpolation") - new{typeof(u), typeof(t), typeof(I), eltype(u)}(u, t, I, mode, extrapolate, Ref(1)) + new{typeof(u), typeof(t), typeof(I), typeof(p.l₀), eltype(u)}( + u, t, I, p, mode, extrapolate, Ref(1), safetycopy) end end -function QuadraticInterpolation(u, t, mode; extrapolate = false) - u, t = munge_data(u, t) - A = QuadraticInterpolation(u, t, nothing, mode, extrapolate) +function QuadraticInterpolation(u, t, mode; extrapolate = false, safetycopy = true) + u, t = munge_data(u, t, safetycopy) + p = QuadraticParameterCache(u, t) + A = QuadraticInterpolation(u, t, nothing, p, mode, extrapolate, safetycopy) I = cumulative_integral(A) - QuadraticInterpolation(u, t, I, mode, extrapolate) + QuadraticInterpolation(u, t, I, p, mode, extrapolate, safetycopy) end -function QuadraticInterpolation(u, t; extrapolate = false) - QuadraticInterpolation(u, t, :Forward; extrapolate) +function QuadraticInterpolation(u, t; extrapolate = false, safetycopy = true) + QuadraticInterpolation(u, t, :Forward; extrapolate, safetycopy) end """ @@ -86,6 +96,7 @@ It is the method of interpolation using Lagrange polynomials of (k-1)th order pa ## Keyword Arguments - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`. """ struct LagrangeInterpolation{uType, tType, T, bcacheType} <: AbstractInterpolation{T} @@ -93,27 +104,33 @@ struct LagrangeInterpolation{uType, tType, T, bcacheType} <: t::tType n::Int bcache::bcacheType + idxs::Vector{Int} extrapolate::Bool idx_prev::Base.RefValue{Int} - function LagrangeInterpolation(u, t, n, extrapolate) + safetycopy::Bool + function LagrangeInterpolation(u, t, n, extrapolate, safetycopy) bcache = zeros(eltype(u[1]), n + 1) + idxs = zeros(Int, n + 1) fill!(bcache, NaN) new{typeof(u), typeof(t), eltype(u), typeof(bcache)}(u, t, n, bcache, + idxs, extrapolate, - Ref(1) + Ref(1), + safetycopy ) end end -function LagrangeInterpolation(u, t, n = length(t) - 1; extrapolate = false) - u, t = munge_data(u, t) +function LagrangeInterpolation( + u, t, n = length(t) - 1; extrapolate = false, safetycopy = true) + u, t = munge_data(u, t, safetycopy) if n != length(t) - 1 error("Currently only n=length(t) - 1 is supported") end - LagrangeInterpolation(u, t, n, extrapolate) + LagrangeInterpolation(u, t, n, extrapolate, safetycopy) end """ @@ -130,6 +147,7 @@ Extrapolation extends the last cubic polynomial on each side. ## Keyword Arguments - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`. """ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <: AbstractInterpolation{T} @@ -141,7 +159,8 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <: d::dType extrapolate::Bool idx_prev::Base.RefValue{Int} - function AkimaInterpolation(u, t, I, b, c, d, extrapolate) + safetycopy::Bool + function AkimaInterpolation(u, t, I, b, c, d, extrapolate, safetycopy) new{typeof(u), typeof(t), typeof(I), typeof(b), typeof(c), typeof(d), eltype(u)}(u, t, @@ -150,13 +169,14 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <: c, d, extrapolate, - Ref(1) + Ref(1), + safetycopy ) end end -function AkimaInterpolation(u, t; extrapolate = false) - u, t = munge_data(u, t) +function AkimaInterpolation(u, t; extrapolate = false, safetycopy = true) + u, t = munge_data(u, t, safetycopy) n = length(t) dt = diff(t) m = Array{eltype(u)}(undef, n + 3) @@ -177,9 +197,9 @@ function AkimaInterpolation(u, t; extrapolate = false) c = (3.0 .* m[3:(end - 2)] .- 2.0 .* b[1:(end - 1)] .- b[2:end]) ./ dt d = (b[1:(end - 1)] .+ b[2:end] .- 2.0 .* m[3:(end - 2)]) ./ dt .^ 2 - A = AkimaInterpolation(u, t, nothing, b, c, d, extrapolate) + A = AkimaInterpolation(u, t, nothing, b, c, d, extrapolate, safetycopy) I = cumulative_integral(A) - AkimaInterpolation(u, t, I, b, c, d, extrapolate) + AkimaInterpolation(u, t, I, b, c, d, extrapolate, safetycopy) end """ @@ -198,24 +218,28 @@ Extrapolation extends the last constant polynomial at the end points on each sid - `dir`: indicates which value should be used for interpolation (`:left` or `:right`). - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`. """ struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T} u::uType t::tType I::IType + p::Nothing dir::Symbol # indicates if value to the $dir should be used for the interpolation extrapolate::Bool idx_prev::Base.RefValue{Int} - function ConstantInterpolation(u, t, I, dir, extrapolate) - new{typeof(u), typeof(t), typeof(I), eltype(u)}(u, t, I, dir, extrapolate, Ref(1)) + safetycopy::Bool + function ConstantInterpolation(u, t, I, dir, extrapolate, safetycopy) + new{typeof(u), typeof(t), typeof(I), eltype(u)}( + u, t, I, nothing, dir, extrapolate, Ref(1), safetycopy) end end -function ConstantInterpolation(u, t; dir = :left, extrapolate = false) - u, t = munge_data(u, t) - A = ConstantInterpolation(u, t, nothing, dir, extrapolate) +function ConstantInterpolation(u, t; dir = :left, extrapolate = false, safetycopy = true) + u, t = munge_data(u, t, safetycopy) + A = ConstantInterpolation(u, t, nothing, dir, extrapolate, safetycopy) I = cumulative_integral(A) - ConstantInterpolation(u, t, I, dir, extrapolate) + ConstantInterpolation(u, t, I, dir, extrapolate, safetycopy) end """ @@ -232,34 +256,40 @@ Extrapolation extends the last quadratic polynomial on each side. ## Keyword Arguments - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`. """ -struct QuadraticSpline{uType, tType, IType, tAType, dType, zType, T} <: +struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <: AbstractInterpolation{T} u::uType t::tType I::IType + p::QuadraticSplineParameterCache{pType} tA::tAType d::dType z::zType extrapolate::Bool idx_prev::Base.RefValue{Int} - function QuadraticSpline(u, t, I, tA, d, z, extrapolate) - new{typeof(u), typeof(t), typeof(I), typeof(tA), + safetycopy::Bool + function QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, safetycopy) + new{typeof(u), typeof(t), typeof(I), typeof(p.σ), typeof(tA), typeof(d), typeof(z), eltype(u)}(u, t, I, + p, tA, d, z, extrapolate, - Ref(1) + Ref(1), + safetycopy ) end end function QuadraticSpline( - u::uType, t; extrapolate = false) where {uType <: AbstractVector{<:Number}} - u, t = munge_data(u, t) + u::uType, t; extrapolate = false, + safetycopy = true) where {uType <: AbstractVector{<:Number}} + u, t = munge_data(u, t, safetycopy) s = length(t) dl = ones(eltype(t), s - 1) d_tmp = ones(eltype(t), s) @@ -271,13 +301,15 @@ function QuadraticSpline( d = map(i -> i == 1 ? typed_zero : 2 // 1 * (u[i] - u[i - 1]) / (t[i] - t[i - 1]), 1:s) z = tA \ d - A = QuadraticSpline(u, t, nothing, tA, d, z, extrapolate) + p = QuadraticSplineParameterCache(z, t) + A = QuadraticSpline(u, t, nothing, p, tA, d, z, extrapolate, safetycopy) I = cumulative_integral(A) - QuadraticSpline(u, t, I, tA, d, z, extrapolate) + QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, safetycopy) end -function QuadraticSpline(u::uType, t; extrapolate = false) where {uType <: AbstractVector} - u, t = munge_data(u, t) +function QuadraticSpline( + u::uType, t; extrapolate = false, safetycopy = true) where {uType <: AbstractVector} + u, t = munge_data(u, t, safetycopy) s = length(t) dl = ones(eltype(t), s - 1) d_tmp = ones(eltype(t), s) @@ -290,9 +322,10 @@ function QuadraticSpline(u::uType, t; extrapolate = false) where {uType <: Abstr d = transpose(reshape(reduce(hcat, d_), :, s)) z_ = reshape(transpose(tA \ d), size(u[1])..., :) z = [z_s for z_s in eachslice(z_, dims = ndims(z_))] - A = QuadraticSpline(u, t, nothing, tA, d, z, extrapolate) + p = QuadraticSplineParameterCache(z, t) + A = QuadraticSpline(u, t, nothing, p, tA, d, z, extrapolate, safetycopy) I = cumulative_integral(A) - QuadraticSpline(u, t, I, tA, d, z, extrapolate) + QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, safetycopy) end """ @@ -309,31 +342,37 @@ Second derivative on both ends are zero, which are also called "natural" boundar ## Keyword Arguments - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`. """ -struct CubicSpline{uType, tType, IType, hType, zType, T} <: AbstractInterpolation{T} +struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInterpolation{T} u::uType t::tType I::IType + p::CubicSplineParameterCache{pType} h::hType z::zType extrapolate::Bool idx_prev::Base.RefValue{Int} - function CubicSpline(u, t, I, h, z, extrapolate) - new{typeof(u), typeof(t), typeof(I), typeof(h), typeof(z), eltype(u)}(u, + safetycopy::Bool + function CubicSpline(u, t, I, p, h, z, extrapolate, safetycopy) + new{typeof(u), typeof(t), typeof(I), typeof(p.c₁), typeof(h), typeof(z), eltype(u)}( + u, t, I, + p, h, z, extrapolate, - Ref(1) + Ref(1), + safetycopy ) end end function CubicSpline(u::uType, t; - extrapolate = false) where {uType <: AbstractVector{<:Number}} - u, t = munge_data(u, t) + extrapolate = false, safetycopy = true) where {uType <: AbstractVector{<:Number}} + u, t = munge_data(u, t, safetycopy) n = length(t) - 1 h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0) dl = vcat(h[2:n], zero(eltype(h))) @@ -350,13 +389,15 @@ function CubicSpline(u::uType, 6(u[i + 1] - u[i]) / h[i + 1] - 6(u[i] - u[i - 1]) / h[i], 1:(n + 1)) z = tA \ d - A = CubicSpline(u, t, nothing, h[1:(n + 1)], z, extrapolate) + p = CubicSplineParameterCache(u, h, z) + A = CubicSpline(u, t, nothing, p, h[1:(n + 1)], z, extrapolate, safetycopy) I = cumulative_integral(A) - CubicSpline(u, t, I, h[1:(n + 1)], z, extrapolate) + CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, safetycopy) end -function CubicSpline(u::uType, t; extrapolate = false) where {uType <: AbstractVector} - u, t = munge_data(u, t) +function CubicSpline( + u::uType, t; extrapolate = false, safetycopy = true) where {uType <: AbstractVector} + u, t = munge_data(u, t, safetycopy) n = length(t) - 1 h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0) dl = vcat(h[2:n], zero(eltype(h))) @@ -370,9 +411,10 @@ function CubicSpline(u::uType, t; extrapolate = false) where {uType <: AbstractV d = transpose(reshape(reduce(hcat, d_), :, n + 1)) z_ = reshape(transpose(tA \ d), size(u[1])..., :) z = [z_s for z_s in eachslice(z_, dims = ndims(z_))] - A = CubicSpline(u, t, nothing, h[1:(n + 1)], z, extrapolate) + p = CubicSplineParameterCache(u, h, z) + A = CubicSpline(u, t, nothing, p, h[1:(n + 1)], z, extrapolate, safetycopy) I = cumulative_integral(A) - CubicSpline(u, t, I, h[1:(n + 1)], z, extrapolate) + CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, safetycopy) end """ @@ -392,8 +434,9 @@ Extrapolation is a constant polynomial of the end points on each side. ## Keyword Arguments - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`. """ -struct BSplineInterpolation{uType, tType, pType, kType, cType, T} <: +struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <: AbstractInterpolation{T} u::uType t::tType @@ -401,35 +444,42 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, T} <: p::pType # params vector k::kType # knot vector c::cType # control points + N::NType # Spline coefficients (preallocated memory) pVecType::Symbol knotVecType::Symbol extrapolate::Bool idx_prev::Base.RefValue{Int} + safetycopy::Bool function BSplineInterpolation(u, t, d, p, k, c, + N, pVecType, knotVecType, - extrapolate) - new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), eltype(u)}(u, + extrapolate, + safetycopy) + new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(N), eltype(u)}(u, t, d, p, k, c, + N, pVecType, knotVecType, extrapolate, - Ref(1) + Ref(1), + safetycopy ) end end -function BSplineInterpolation(u, t, d, pVecType, knotVecType; extrapolate = false) - u, t = munge_data(u, t) +function BSplineInterpolation( + u, t, d, pVecType, knotVecType; extrapolate = false, safetycopy = true) + u, t = munge_data(u, t, safetycopy) n = length(t) s = zero(eltype(u)) p = zero(t) @@ -487,9 +537,12 @@ function BSplineInterpolation(u, t, d, pVecType, knotVecType; extrapolate = fals end end # control points - N = spline_coefficients(n, d, k, p) + N = zeros(eltype(t), n, n) + spline_coefficients!(N, d, k, p) c = vec(N \ u[:, :]) - BSplineInterpolation(u, t, d, p, k, c, pVecType, knotVecType, extrapolate) + N = zeros(eltype(t), n) + BSplineInterpolation( + u, t, d, p, k, c, N, pVecType, knotVecType, extrapolate, safetycopy) end """ @@ -511,8 +564,9 @@ Extrapolation is a constant polynomial of the end points on each side. ## Keyword Arguments - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`. """ -struct BSplineApprox{uType, tType, pType, kType, cType, T} <: +struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <: AbstractInterpolation{T} u::uType t::tType @@ -521,10 +575,12 @@ struct BSplineApprox{uType, tType, pType, kType, cType, T} <: p::pType # params vector k::kType # knot vector c::cType # control points + N::NType # Spline coefficients (preallocated memory) pVecType::Symbol knotVecType::Symbol extrapolate::Bool idx_prev::Base.RefValue{Int} + safetycopy::Bool function BSplineApprox(u, t, d, @@ -532,26 +588,32 @@ struct BSplineApprox{uType, tType, pType, kType, cType, T} <: p, k, c, + N, pVecType, knotVecType, - extrapolate) - new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), eltype(u)}(u, + extrapolate, + safetycopy + ) + new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(N), eltype(u)}(u, t, d, h, p, k, c, + N, pVecType, knotVecType, extrapolate, - Ref(1) + Ref(1), + safetycopy::Bool ) end end -function BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolate = false) - u, t = munge_data(u, t) +function BSplineApprox( + u, t, d, h, pVecType, knotVecType; extrapolate = false, safetycopy = true) + u, t = munge_data(u, t, safetycopy) n = length(t) s = zero(eltype(u)) p = zero(t) @@ -616,7 +678,7 @@ function BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolate = false) q = zeros(eltype(u), n) N = zeros(eltype(t), n, h) for i in 1:n - N[i, :] .= spline_coefficients(h, d, k, p[i]) + spline_coefficients!(view(N, i, :), d, k, p[i]) end for k in 2:(n - 1) q[k] = u[k] - N[k, 1] * u[1] - N[k, h] * u[end] @@ -633,7 +695,8 @@ function BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolate = false) M = transpose(N) * N P = M \ Q c[2:(end - 1)] .= vec(P) - BSplineApprox(u, t, d, h, p, k, c, pVecType, knotVecType, extrapolate) + N = zeros(eltype(t), h) + BSplineApprox(u, t, d, h, p, k, c, N, pVecType, knotVecType, extrapolate, safetycopy) end """ @@ -650,28 +713,30 @@ It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomi ## Keyword Arguments - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`. """ struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInterpolation{T} - du::uType + du::duType u::uType t::tType I::IType p::CubicHermiteParameterCache{pType} extrapolate::Bool idx_prev::Base.RefValue{Int} - function CubicHermiteSpline(du, u, t, I, p, extrapolate) + safetycopy::Bool + function CubicHermiteSpline(du, u, t, I, p, extrapolate, safetycopy) new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u)}( - du, u, t, I, p, extrapolate, Ref(1)) + du, u, t, I, p, extrapolate, Ref(1), safetycopy) end end -function CubicHermiteSpline(du, u, t; extrapolate = false) +function CubicHermiteSpline(du, u, t; extrapolate = false, safetycopy = true) @assert length(u)==length(du) "Length of `u` is not equal to length of `du`." - u, t = munge_data(u, t) + u, t = munge_data(u, t, safetycopy) p = CubicHermiteParameterCache(du, u, t) - A = CubicHermiteSpline(du, u, t, nothing, p, extrapolate) + A = CubicHermiteSpline(du, u, t, nothing, p, extrapolate, safetycopy) I = cumulative_integral(A) - CubicHermiteSpline(du, u, t, I, p, extrapolate) + CubicHermiteSpline(du, u, t, I, p, extrapolate, safetycopy) end """ @@ -689,29 +754,31 @@ It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polyno ## Keyword Arguments - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. + - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`. """ struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <: AbstractInterpolation{T} - ddu::uType - du::uType + ddu::dduType + du::duType u::uType t::tType I::IType p::QuinticHermiteParameterCache{pType} extrapolate::Bool idx_prev::Base.RefValue{Int} - function QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate) + safetycopy::Bool + function QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate, safetycopy) 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)) + ddu, du, u, t, I, p, extrapolate, Ref(1), safetycopy) end end -function QuinticHermiteSpline(ddu, du, u, t; extrapolate = false) +function QuinticHermiteSpline(ddu, du, u, t; extrapolate = false, safetycopy = true) @assert length(u)==length(du)==length(ddu) "Length of `u` is not equal to length of `du` or `ddu`." - u, t = munge_data(u, t) + u, t = munge_data(u, t, safetycopy) p = QuinticHermiteParameterCache(ddu, du, u, t) - A = QuinticHermiteSpline(ddu, du, u, t, nothing, p, extrapolate) + A = QuinticHermiteSpline(ddu, du, u, t, nothing, p, extrapolate, safetycopy) I = cumulative_integral(A) - QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate) + QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate, safetycopy) end diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 0a543d8f..5d09ceff 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -13,23 +13,32 @@ function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, igues idx = firstindex(A.u) - 1 t1 = t2 = one(eltype(A.t)) u1 = u2 = one(eltype(A.u)) + slope = t * one(eltype(A.p.slope)) else idx = get_idx(A.t, t, iguess) t1, t2 = A.t[idx], A.t[idx + 1] u1, u2 = A.u[idx], A.u[idx + 1] + slope = A.p.slope[idx] end - θ = (t - t1) / (t2 - t1) - val = (1 - θ) * u1 + θ * u2 - # Note: The following is limited to when val is NaN as to not change the derivative of exact points. - t == t1 && any(isnan, val) && return oftype(val, u1), idx # Return exact value if no interpolation needed (eg when NaN at t2) - t == t2 && any(isnan, val) && return oftype(val, u2), idx # ... (eg when NaN at t1) + + Δt = t - t1 + Δu = slope * Δt + val = u1 + Δu_nan = any(isnan.(Δu)) + if t == t2 && Δu_nan + val = u2 + elseif !(iszero(Δt) && Δu_nan) + val += Δu + end + val = oftype(Δu, val) + val, idx end function _interpolate(A::LinearInterpolation{<:AbstractMatrix}, t::Number, iguess) idx = get_idx(A.t, t, iguess) - θ = (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) - return (1 - θ) * A.u[:, idx] + θ * A.u[:, idx + 1], idx + Δt = t - A.t[idx] + return A.u[:, idx] + A.p.slope[idx] * Δt, idx end # Quadratic Interpolation @@ -39,87 +48,71 @@ function _quad_interp_indices(A::QuadraticInterpolation, t::Number, iguess) idx, idx + 1, idx + 2 end -function _interpolate(A::QuadraticInterpolation{<:AbstractVector}, t::Number, iguess) +function _interpolate(A::QuadraticInterpolation, t::Number, iguess) i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess) - l₀ = ((t - A.t[i₁]) * (t - A.t[i₂])) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂])) - l₁ = ((t - A.t[i₀]) * (t - A.t[i₂])) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂])) - l₂ = ((t - A.t[i₀]) * (t - A.t[i₁])) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁])) - return A.u[i₀] * l₀ + A.u[i₁] * l₁ + A.u[i₂] * l₂, i₀ -end - -function _interpolate(A::QuadraticInterpolation{<:AbstractMatrix}, t::Number, iguess) - i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess) - l₀ = ((t - A.t[i₁]) * (t - A.t[i₂])) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂])) - l₁ = ((t - A.t[i₀]) * (t - A.t[i₂])) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂])) - l₂ = ((t - A.t[i₀]) * (t - A.t[i₁])) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁])) - return A.u[:, i₀] * l₀ + A.u[:, i₁] * l₁ + A.u[:, i₂] * l₂, i₀ + u₀ = A.p.l₀[i₀] * (t - A.t[i₁]) * (t - A.t[i₂]) + u₁ = A.p.l₁[i₀] * (t - A.t[i₀]) * (t - A.t[i₂]) + u₂ = A.p.l₂[i₀] * (t - A.t[i₀]) * (t - A.t[i₁]) + return u₀ + u₁ + u₂, i₀ end # Lagrange Interpolation -function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number) - ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) - idxs = findRequiredIdxs(A, t) - if A.t[idxs[1]] == t - return A.u[idxs[1]] +function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, iguess) + idx = get_idx(A.t, t, iguess) + findRequiredIdxs!(A, t, idx) + if A.t[A.idxs[1]] == t + return A.u[A.idxs[1]], idx end N = zero(A.u[1]) D = zero(A.t[1]) tmp = N - for i in 1:length(idxs) - if isnan(A.bcache[idxs[i]]) + for i in 1:length(A.idxs) + if isnan(A.bcache[A.idxs[i]]) mult = one(A.t[1]) for j in 1:(i - 1) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + mult *= (A.t[A.idxs[i]] - A.t[A.idxs[j]]) end - for j in (i + 1):length(idxs) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + for j in (i + 1):length(A.idxs) + mult *= (A.t[A.idxs[i]] - A.t[A.idxs[j]]) end - A.bcache[idxs[i]] = mult + A.bcache[A.idxs[i]] = mult else - mult = A.bcache[idxs[i]] + mult = A.bcache[A.idxs[i]] end - tmp = inv((t - A.t[idxs[i]]) * mult) + tmp = inv((t - A.t[A.idxs[i]]) * mult) D += tmp - N += (tmp * A.u[idxs[i]]) + N += (tmp * A.u[A.idxs[i]]) end - N / D + N / D, idx end -function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number) - ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) - idxs = findRequiredIdxs(A, t) - if A.t[idxs[1]] == t - return A.u[:, idxs[1]] +function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, iguess) + idx = get_idx(A.t, t, iguess) + findRequiredIdxs!(A, t, idx) + if A.t[A.idxs[1]] == t + return A.u[:, A.idxs[1]], idx end N = zero(A.u[:, 1]) D = zero(A.t[1]) tmp = D - for i in 1:length(idxs) - if isnan(A.bcache[idxs[i]]) + for i in 1:length(A.idxs) + if isnan(A.bcache[A.idxs[i]]) mult = one(A.t[1]) for j in 1:(i - 1) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + mult *= (A.t[A.idxs[i]] - A.t[A.idxs[j]]) end - for j in (i + 1):length(idxs) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + for j in (i + 1):length(A.idxs) + mult *= (A.t[A.idxs[i]] - A.t[A.idxs[j]]) end - A.bcache[idxs[i]] = mult + A.bcache[A.idxs[i]] = mult else - mult = A.bcache[idxs[i]] + mult = A.bcache[A.idxs[i]] end - tmp = inv((t - A.t[idxs[i]]) * mult) + tmp = inv((t - A.t[A.idxs[i]]) * mult) D += tmp - @. N += (tmp * A.u[:, idxs[i]]) + @. N += (tmp * A.u[:, A.idxs[i]]) end - N / D -end - -function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, idx) - _interpolate(A, t), idx -end - -function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, idx) - _interpolate(A, t), idx + N / D, idx end function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess) @@ -153,19 +146,20 @@ end # QuadraticSpline Interpolation function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) - idx = get_idx(A.t, t, iguess; lb = 2, ub_shift = 0, side = :first) - Cᵢ = A.u[idx - 1] - σ = 1 // 2 * (A.z[idx] - A.z[idx - 1]) / (A.t[idx] - A.t[idx - 1]) - return A.z[idx - 1] * (t - A.t[idx - 1]) + σ * (t - A.t[idx - 1])^2 + Cᵢ, idx + idx = get_idx(A.t, t, iguess) + Cᵢ = A.u[idx] + Δt = t - A.t[idx] + return A.z[idx] * Δt + A.p.σ[idx] * Δt^2 + Cᵢ, idx end # CubicSpline Interpolation function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess) idx = get_idx(A.t, t, iguess) - I = A.z[idx] * (A.t[idx + 1] - t)^3 / (6A.h[idx + 1]) + - A.z[idx + 1] * (t - A.t[idx])^3 / (6A.h[idx + 1]) - C = (A.u[idx + 1] / A.h[idx + 1] - A.z[idx + 1] * A.h[idx + 1] / 6) * (t - A.t[idx]) - D = (A.u[idx] / A.h[idx + 1] - A.z[idx] * A.h[idx + 1] / 6) * (A.t[idx + 1] - t) + Δt₁ = t - A.t[idx] + Δt₂ = A.t[idx + 1] - t + I = (A.z[idx] * Δt₂^3 + A.z[idx + 1] * Δt₁^3) / (6A.h[idx + 1]) + C = A.p.c₁[idx] * Δt₁ + D = A.p.c₂[idx] * Δt₂ I + C + D, idx end @@ -179,9 +173,10 @@ function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}}, idx = get_idx(A.t, t, iguess) t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx]) n = length(A.t) - N = spline_coefficients(n, A.d, A.k, t) + N = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.N + nonzero_coefficient_idxs = spline_coefficients!(N, A.d, A.k, t) ucum = zero(eltype(A.u)) - for i in 1:n + for i in nonzero_coefficient_idxs ucum += N[i] * A.c[i] end ucum, idx @@ -194,9 +189,10 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i # change t into param [0 1] idx = get_idx(A.t, t, iguess) t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx]) - N = spline_coefficients(A.h, A.d, A.k, t) + N = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.N + nonzero_coefficient_idxs = spline_coefficients!(N, A.d, A.k, t) ucum = zero(eltype(A.u)) - for i in 1:(A.h) + for i in nonzero_coefficient_idxs ucum += N[i] * A.c[i] end ucum, idx diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 6eb00eab..466248b1 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -1,14 +1,42 @@ -function findRequiredIdxs(A::LagrangeInterpolation, t) - idxs = sortperm(A.t, by = x -> abs(t - x)) - idxs[1:(A.n + 1)] +function findRequiredIdxs!(A::LagrangeInterpolation, t, idx) + n = length(A.t) - 1 + i_min, idx_min, idx_max = if t == A.t[idx] + A.idxs[1] = idx + 2, idx, idx + else + 1, idx + 1, idx + end + for i in i_min:(n + 1) + if idx_min == 1 + A.idxs[i:end] .= range(idx_max + 1, idx_max + (n + 2 - i)) + break + elseif idx_max == length(A.t) + A.idxs[i:end] .= (idx_min - 1):-1:(idx_min - (n + 2 - i)) + break + else + left_diff = abs(t - A.t[idx_min - 1]) + right_diff = abs(t - A.t[idx_max + 1]) + left_expand = left_diff <= right_diff + end + if left_expand + idx_min -= 1 + A.idxs[i] = idx_min + else + idx_max += 1 + A.idxs[i] = idx_max + end + end + return idx end -function spline_coefficients(n, d, k, u::Number) - N = zeros(eltype(u), n) +function spline_coefficients!(N, d, k, u::Number) + N .= zero(u) if u == k[1] N[1] = one(u) + return 1:1 elseif u == k[end] N[end] = one(u) + return length(N):length(N) else i = findfirst(x -> x > u, k) - 1 N[i] = one(u) @@ -20,51 +48,68 @@ function spline_coefficients(n, d, k, u::Number) end N[i] = (u - k[i]) / (k[i + deg] - k[i]) * N[i] end + return (i - d):i end - N end -function spline_coefficients(n, d, k, u::AbstractVector) - N = zeros(eltype(u), n, n) - for i in 1:n - N[i, :] .= spline_coefficients(n, d, k, u[i]) +function spline_coefficients!(N, d, k, u::AbstractVector) + for i in 1:size(N)[2] + spline_coefficients!(view(N, i, :), d, k, u[i]) end - N + return nothing end # helper function for data manipulation -function munge_data(u::AbstractVector{<:Real}, t::AbstractVector{<:Real}) - return u, t +function munge_data(u::AbstractVector{<:Real}, t::AbstractVector{<:Real}, safetycopy::Bool) + if safetycopy + u = copy(u) + t = copy(t) + end + return readonly_wrap(u), readonly_wrap(t) end -function munge_data(u::AbstractVector, t::AbstractVector) +function munge_data(u::AbstractVector, t::AbstractVector, safetycopy::Bool) Tu = Base.nonmissingtype(eltype(u)) Tt = Base.nonmissingtype(eltype(t)) @assert length(t) == length(u) non_missing_indices = collect( - idx for idx in 1:length(t) - if !ismissing(u[idx]) && !ismissing(t[idx]) + i for i in 1:length(t) + if !ismissing(u[i]) && !ismissing(t[i]) ) - newu = Tu.([u[idx] for idx in non_missing_indices]) - newt = Tt.([t[idx] for idx in non_missing_indices]) - return newu, newt + if safetycopy + u = Tu.([u[i] for i in non_missing_indices]) + t = Tt.([t[i] for i in non_missing_indices]) + else + !isempty(non_missing_indices) && throw(MustCopyError()) + end + + return readonly_wrap(u), readonly_wrap(t) end -function munge_data(U::StridedMatrix, t::AbstractVector) +function munge_data(U::StridedMatrix, t::AbstractVector, safetycopy::Bool) TU = Base.nonmissingtype(eltype(U)) Tt = Base.nonmissingtype(eltype(t)) @assert length(t) == size(U, 2) non_missing_indices = collect( - idx for idx in 1:length(t) - if !any(ismissing, U[:, idx]) && !ismissing(t[idx]) + i for i in 1:length(t) + if !any(ismissing, U[:, i]) && !ismissing(t[i]) ) - newUs = [TU.(U[:, idx]) for idx in non_missing_indices] - newt = Tt.([t[idx] for idx in non_missing_indices]) - return hcat(newUs...), newt + if safetycopy + U = hcat([TU.(U[:, i]) for i in non_missing_indices]...) + t = Tt.([t[i] for i in non_missing_indices]) + else + !isempty(non_missing_indices) && throw(MustCopyError()) + end + + return readonly_wrap(U), readonly_wrap(t) end +# Don't nest ReadOnlyArrays +readonly_wrap(a::AbstractArray) = ReadOnlyArray(a) +readonly_wrap(a::ReadOnlyArray) = a + function get_idx(tvec, t, iguess; lb = 1, ub_shift = -1, idx_shift = 0, side = :last) ub = length(tvec) + ub_shift return if side == :last diff --git a/src/online.jl b/src/online.jl index 4ad06f3a..dc500611 100644 --- a/src/online.jl +++ b/src/online.jl @@ -1,40 +1,63 @@ import Base: append!, push! -function push!(di::LinearInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T} - push!(di.u, u) - push!(di.t, t) - di +function push!(A::LinearInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T} + push!(A.u.parent, u) + push!(A.t.parent, t) + slope = linear_interpolation_parameters(A.u, A.t, length(A.t) - 1) + push!(A.p.slope, slope) + A end -function push!(di::QuadraticInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T} - push!(di.u, u) - push!(di.t, t) - di +function push!(A::QuadraticInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T} + push!(A.u.parent, u) + push!(A.t.parent, t) + l₀, l₁, l₂ = quadratic_interpolation_parameters(A.u, A.t, length(A.t) - 2) + push!(A.p.l₀, l₀) + push!(A.p.l₁, l₁) + push!(A.p.l₂, l₂) + A end -function push!(di::ConstantInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T} - push!(di.u, u) - push!(di.t, t) - di +function push!(A::ConstantInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T} + push!(A.u.parent, u) + push!(A.t.parent, t) + A end -function append!(di::LinearInterpolation{U, T}, u::U, t::T) where {U, T} - u, t = munge_data(u, t) - append!(di.u, u) - append!(di.t, t) - di +function append!( + A::LinearInterpolation{ReadOnlyVector{eltypeU, U}, ReadOnlyVector{eltypeT, T}}, u::U, t::T) where { + eltypeU, U, eltypeT, T} + length_old = length(A.t) + u, t = munge_data(u, t, true) + append!(A.u.parent, u) + append!(A.t.parent, t) + slope = linear_interpolation_parameters.( + Ref(A.u), Ref(A.t), length_old:(length(A.t) - 1)) + append!(A.p.slope, slope) + A end -function append!(di::QuadraticInterpolation{U, T}, u::U, t::T) where {U, T} - u, t = munge_data(u, t) - append!(di.u, u) - append!(di.t, t) - di +function append!( + A::ConstantInterpolation{ReadOnlyVector{eltypeU, U}, ReadOnlyVector{eltypeT, T}}, u::U, t::T) where { + eltypeU, U, eltypeT, T} + u, t = munge_data(u, t, true) + append!(A.u.parent, u) + append!(A.t.parent, t) + A end -function append!(di::ConstantInterpolation{U, T}, u::U, t::T) where {U, T} - u, t = munge_data(u, t) - append!(di.u, u) - append!(di.t, t) - di +function append!( + A::QuadraticInterpolation{ReadOnlyVector{eltypeU, U}, ReadOnlyVector{eltypeT, T}}, u::U, t::T) where { + eltypeU, U, eltypeT, T} + length_old = length(A.t) + u, t = munge_data(u, t, true) + append!(A.u.parent, u) + append!(A.t.parent, t) + parameters = quadratic_interpolation_parameters.( + Ref(A.u), Ref(A.t), (length_old - 1):(length(A.t) - 2)) + l₀, l₁, l₂ = collect.(eachrow(hcat(collect.(parameters)...))) + append!(A.p.l₀, l₀) + append!(A.p.l₁, l₁) + append!(A.p.l₂, l₂) + A end diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index 66aa218f..2820dc8f 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -1,16 +1,100 @@ +struct LinearParameterCache{pType} + slope::pType +end + +function LinearParameterCache(u, t) + slope = linear_interpolation_parameters.(Ref(u), Ref(t), 1:(length(t) - 1)) + return LinearParameterCache(slope) +end + +function linear_interpolation_parameters(u, t, idx) + Δu = u isa AbstractMatrix ? u[:, idx + 1] - u[:, idx] : u[idx + 1] - u[idx] + Δt = t[idx + 1] - t[idx] + slope = Δu / Δt + slope = iszero(Δt) ? zero(slope) : slope + return slope +end + +struct QuadraticParameterCache{pType} + l₀::pType + l₁::pType + l₂::pType +end + +function QuadraticParameterCache(u, t) + parameters = quadratic_interpolation_parameters.( + Ref(u), Ref(t), 1:(length(t) - 2)) + l₀, l₁, l₂ = collect.(eachrow(hcat(collect.(parameters)...))) + return QuadraticParameterCache(l₀, l₁, l₂) +end + +function quadratic_interpolation_parameters(u, t, idx) + if u isa AbstractMatrix + u₀ = u[:, idx] + u₁ = u[:, idx + 1] + u₂ = u[:, idx + 2] + else + u₀ = u[idx] + u₁ = u[idx + 1] + u₂ = u[idx + 2] + end + t₀ = t[idx] + t₁ = t[idx + 1] + t₂ = t[idx + 2] + Δt₀ = t₁ - t₀ + Δt₁ = t₂ - t₁ + Δt₂ = t₂ - t₀ + l₀ = u₀ / (Δt₀ * Δt₂) + l₁ = -u₁ / (Δt₀ * Δt₁) + l₂ = u₂ / (Δt₂ * Δt₁) + return l₀, l₁, l₂ +end + +struct QuadraticSplineParameterCache{pType} + σ::pType +end + +function QuadraticSplineParameterCache(z, t) + σ = quadratic_spline_parameters.(Ref(z), Ref(t), 1:(length(t) - 1)) + return QuadraticSplineParameterCache(σ) +end + +function quadratic_spline_parameters(z, t, idx) + σ = 1 // 2 * (z[idx + 1] - z[idx]) / (t[idx + 1] - t[idx]) + return σ +end + +struct CubicSplineParameterCache{pType} + c₁::pType + c₂::pType +end + +function CubicSplineParameterCache(u, h, z) + parameters = cubic_spline_parameters.( + Ref(u), Ref(h), Ref(z), 1:(size(u)[end] - 1)) + c₁, c₂ = collect.(eachrow(hcat(collect.(parameters)...))) + return CubicSplineParameterCache(c₁, c₂) +end + +function cubic_spline_parameters(u, h, z, idx) + c₁ = (u[idx + 1] / h[idx + 1] - z[idx + 1] * h[idx + 1] / 6) + c₂ = (u[idx] / h[idx + 1] - z[idx] * h[idx + 1] / 6) + return c₁, c₂ +end + struct CubicHermiteParameterCache{pType} c₁::pType c₂::pType end function CubicHermiteParameterCache(du, u, t) - parameters = CubicHermiteSplineParameters.( + parameters = cubic_hermite_spline_parameters.( Ref(du), Ref(u), Ref(t), 1:(length(t) - 1)) c₁, c₂ = collect.(eachrow(hcat(collect.(parameters)...))) return CubicHermiteParameterCache(c₁, c₂) end -function CubicHermiteSplineParameters(du, u, t, idx) +function cubic_hermite_spline_parameters(du, u, t, idx) Δt = t[idx + 1] - t[idx] u₀ = u[idx] u₁ = u[idx + 1] @@ -28,13 +112,13 @@ struct QuinticHermiteParameterCache{pType} end function QuinticHermiteParameterCache(ddu, du, u, t) - parameters = QuinticHermiteSplineParameters.( + parameters = quintic_hermite_spline_parameters.( Ref(ddu), Ref(du), Ref(u), Ref(t), 1:(length(t) - 1)) c₁, c₂, c₃ = collect.(eachrow(hcat(collect.(parameters)...))) return QuinticHermiteParameterCache(c₁, c₂, c₃) end -function QuinticHermiteSplineParameters(ddu, du, u, t, idx) +function quintic_hermite_spline_parameters(ddu, du, u, t, idx) Δt = t[idx + 1] - t[idx] u₀ = u[idx] u₁ = u[idx + 1] diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 225ad9f3..37351d0d 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -125,6 +125,13 @@ end end end +@testset "Constant Interpolation" begin + u = [0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0] + t = collect(0.0:10.0) + A = ConstantInterpolation(u, t) + @test all(derivative.(Ref(A), t) .== 0.0) +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 a5641dff..3d59a6c0 100644 --- a/test/integral_tests.jl +++ b/test/integral_tests.jl @@ -83,7 +83,7 @@ end @testset "LagrangeInterpolation" begin u = [1.0, 4.0, 9.0] - t = [1.0, 2.0, 3.0] + t = [1.0, 2.0, 6.0] A = LagrangeInterpolation(u, t) @test_throws DataInterpolations.IntegralNotFoundError integral(A, 1.0, 2.0) @test_throws DataInterpolations.IntegralNotFoundError integral(A, 5.0) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index f9ef87f8..779c5bd8 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -3,6 +3,18 @@ using FindFirstFunctions: searchsortedfirstcorrelated using StableRNGs using Optim, ForwardDiff +function test_interpolation_type(T) + @test T <: DataInterpolations.AbstractInterpolation + @test hasfield(T, :u) + @test hasfield(T, :t) + @test hasfield(T, :extrapolate) + @test hasfield(T, :idx_prev) + @test hasfield(T, :safetycopy) + @test !isempty(methods(DataInterpolations._interpolate, (T, Any, Number))) + @test !isempty(methods(DataInterpolations._integral, (T, Any, Number))) + @test !isempty(methods(DataInterpolations._derivative, (T, Any, Number))) +end + function test_cached_index(A) for t in range(first(A.t), last(A.t); length = 2 * length(A.t) - 1) A(t) @@ -13,6 +25,8 @@ function test_cached_index(A) end @testset "Linear Interpolation" begin + test_interpolation_type(LinearInterpolation) + for t in (1.0:10.0, 1.0collect(1:10)) u = 2.0collect(1:10) #t = 1.0collect(1:10) @@ -75,14 +89,6 @@ end @test A(3.0) == 2.0 @test A(4.0) == 3.0 - u = [0.0, NaN, 2.0, 3.0] - A = LinearInterpolation(u, t; extrapolate = true) - @test A(1.0) == 0.0 - @test isnan(A(2.0)) - @test isnan(A(2.5)) - @test A(3.0) == 2.0 - @test A(4.0) == 3.0 - u = [0.0, 1.0, NaN, 3.0] A = LinearInterpolation(u, t; extrapolate = true) @test A(1.0) == 0.0 @@ -167,6 +173,8 @@ end end @testset "Quadratic Interpolation" begin + test_interpolation_type(QuadraticInterpolation) + u = [1.0, 4.0, 9.0, 16.0] t = [1.0, 2.0, 3.0, 4.0] A = QuadraticInterpolation(u, t; extrapolate = true) @@ -262,6 +270,8 @@ end end @testset "Lagrange Interpolation" begin + test_interpolation_type(LagrangeInterpolation) + u = [1.0, 4.0, 9.0] t = [1.0, 2.0, 3.0] A = LagrangeInterpolation(u, t) @@ -320,6 +330,8 @@ end end @testset "Akima Interpolation" begin + test_interpolation_type(AkimaInterpolation) + u = [0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0] t = collect(0.0:10.0) A = AkimaInterpolation(u, t) @@ -349,6 +361,8 @@ end end @testset "ConstantInterpolation" begin + test_interpolation_type(ConstantInterpolation) + t = [1.0, 2.0, 3.0, 4.0] @testset "Vector case" for u in [[1.0, 2.0, 0.0, 1.0], ["B", "C", "A", "B"]] @@ -473,6 +487,8 @@ end end @testset "QuadraticSpline Interpolation" begin + test_interpolation_type(QuadraticSpline) + u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] @@ -518,6 +534,8 @@ end end @testset "CubicSpline Interpolation" begin + test_interpolation_type(CubicSpline) + u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] @@ -569,11 +587,13 @@ end end @testset "BSplines" begin + # BSpline Interpolation and Approximation t = [0, 62.25, 109.66, 162.66, 205.8, 252.3] u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22] @testset "BSplineInterpolation" begin + test_interpolation_type(BSplineInterpolation) A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform) @test [A(25.0), A(80.0)] == [13.454197730061425, 10.305633616059845] @@ -605,6 +625,7 @@ end end @testset "BSplineApprox" begin + test_interpolation_type(BSplineApprox) A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform) @test [A(25.0), A(80.0)] ≈ [12.979802931218234, 10.914310609953178] @@ -623,6 +644,8 @@ end end @testset "Cubic Hermite Spline" begin + test_interpolation_type(CubicHermiteSpline) + du = [-0.047, -0.058, 0.054, 0.012, -0.068, 0.0] u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22] t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3] @@ -636,6 +659,8 @@ end end @testset "Quintic Hermite Spline" begin + test_interpolation_type(QuinticHermiteSpline) + ddu = [0.0, -0.00033, 0.0051, -0.0067, 0.0029, 0.0] du = [-0.047, -0.058, 0.054, 0.012, -0.068, 0.0] u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22] diff --git a/test/online_tests.jl b/test/online_tests.jl index cc4d4c84..571c2dc0 100644 --- a/test/online_tests.jl +++ b/test/online_tests.jl @@ -1,18 +1,31 @@ using DataInterpolations -t = [1, 2, 3] -u = [0, 1, 0] +t1 = [1.0, 2.0, 3.0] +u1 = [0.0, 1.0, 0.0] -for di in [LinearInterpolation, QuadraticInterpolation, ConstantInterpolation] - li = di(copy(u), copy(t)) - append!(li, u, t) - li2 = di(vcat(u, u), vcat(t, t)) - @test li.u == li2.u - @test li.t == li2.t +t2 = [4.0, 5.0, 6.0] +u2 = [1.0, 2.0, 1.0] - li = di(copy(u), copy(t)) - push!(li, 1, 4) - li2 = di(vcat(u, 1), vcat(t, 4)) - @test li.u == li2.u - @test li.t == li2.t +ts = 1.0:0.5:6.0 + +for method in [LinearInterpolation, QuadraticInterpolation, ConstantInterpolation] + func1 = method(u1, t1) + append!(func1, u2, t2) + func2 = method(vcat(u1, u2), vcat(t1, t2)) + @test func1.u == func2.u + @test func1.t == func2.t + for name in propertynames(func1.p) + @test getfield(func1.p, name) == getfield(func2.p, name) + end + @test func1(ts) == func2(ts) + + func1 = method(u1, t1) + push!(func1, 1.0, 4.0) + func2 = method(vcat(u1, 1.0), vcat(t1, 4.0)) + @test func1.u == func2.u + @test func1.t == func2.t + for name in propertynames(func1.p) + @test getfield(func1.p, name) == getfield(func2.p, name) + end + @test func1(ts) == func2(ts) end diff --git a/test/parameter_tests.jl b/test/parameter_tests.jl index 3bedf3bc..bcd26cf7 100644 --- a/test/parameter_tests.jl +++ b/test/parameter_tests.jl @@ -1,5 +1,36 @@ using DataInterpolations +@testset "Linear Interpolation" begin + u = [1.0, 5.0, 3.0, 4.0, 4.0] + t = collect(1:5) + A = LinearInterpolation(u, t) + @test A.p.slope ≈ [4.0, -2.0, 1.0, 0.0] +end + +@testset "Quadratic Interpolation" begin + u = [1.0, 5.0, 3.0, 4.0, 4.0] + t = collect(1:5) + A = QuadraticInterpolation(u, t) + @test A.p.l₀ ≈ [0.5, 2.5, 1.5] + @test A.p.l₁ ≈ [-5.0, -3.0, -4.0] + @test A.p.l₂ ≈ [1.5, 2.0, 2.0] +end + +@testset "Quadratic Spline" begin + u = [1.0, 5.0, 3.0, 4.0, 4.0] + t = collect(1:5) + A = QuadraticSpline(u, t) + @test A.p.σ ≈ [4.0, -10.0, 13.0, -14.0] +end + +@testset "Cubic Spline" begin + u = [1, 5, 3, 4, 4] + t = collect(1:5) + A = CubicSpline(u, t) + @test A.p.c₁ ≈ [6.839285714285714, 1.642857142857143, 4.589285714285714, 4.0] + @test A.p.c₂ ≈ [1.0, 6.839285714285714, 1.642857142857143, 4.589285714285714] +end + @testset "Cubic Hermite Spline" begin du = [5.0, 3.0, 6.0, 8.0, 1.0] u = [1.0, 5.0, 3.0, 4.0, 4.0]