Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix integrating with cache_parameters = true #391

Merged
merged 3 commits into from
Feb 10, 2025
Merged

Fix integrating with cache_parameters = true #391

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
bc8f9e2
Do the fix, test the fix
SouthEndMusic Feb 6, 2025
5bb496a
nit
SouthEndMusic Feb 6, 2025
7985ae3
Add `cache_parameters` support for `RegularizationSmooth`
SouthEndMusic Feb 10, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
80 changes: 60 additions & 20 deletions ext/DataInterpolationsRegularizationToolsExt.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
module DataInterpolationsRegularizationToolsExt

using DataInterpolations
import DataInterpolations: munge_data,
import DataInterpolations: munge_data, munge_extrapolation,
_interpolate, RegularizationSmooth, get_show, derivative,
integral
using LinearAlgebra
Expand Down Expand Up @@ -70,14 +70,19 @@ 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,
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
cache_parameters::Bool = false)
extrapolation_left, extrapolation_right = munge_extrapolation(
extrapolation, extrapolation_left, extrapolation_right)
u, t = munge_data(u, t)
M = _mapping_matrix(t̂, t)
Wls½ = LA.diagm(sqrt.(wls))
Wr½ = LA.diagm(sqrt.(wr))
û, λ, Aitp = _reg_smooth_solve(
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
extrapolation_right, cache_parameters)
RegularizationSmooth(
u, û, t, t̂, wls, wr, d, λ, alg, Aitp, extrapolation_left, extrapolation_right)
end
Expand All @@ -90,16 +95,22 @@ 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, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
alg::Symbol = :gcv_svd,
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
cache_parameters::Bool = false)
extrapolation_left, extrapolation_right = munge_extrapolation(
extrapolation, extrapolation_left, extrapolation_right)
u, t = munge_data(u, t)
t̂ = t
N = length(t)
M = Array{Float64}(LA.I, N, N)
Wls½ = Array{Float64}(LA.I, N, N)
Wr½ = Array{Float64}(LA.I, N - d, N - d)
û, λ, Aitp = _reg_smooth_solve(
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
extrapolation_right, cache_parameters)
RegularizationSmooth(u,
û,
t,
Expand All @@ -122,15 +133,20 @@ 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,
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
cache_parameters::Bool = false)
extrapolation_left, extrapolation_right = munge_extrapolation(
extrapolation, extrapolation_left, extrapolation_right)
u, t = munge_data(u, t)
N, N̂ = length(t), length(t̂)
M = _mapping_matrix(t̂, t)
Wls½ = Array{Float64}(LA.I, N, N)
Wr½ = Array{Float64}(LA.I, N̂ - d, N̂ - d)
û, λ, Aitp = _reg_smooth_solve(
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
extrapolation_right, cache_parameters)
RegularizationSmooth(u,
û,
t,
Expand All @@ -153,15 +169,21 @@ 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, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
alg::Symbol = :gcv_svd,
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
cache_parameters::Bool = false)
extrapolation_left, extrapolation_right = munge_extrapolation(
extrapolation, extrapolation_left, extrapolation_right)
u, t = munge_data(u, t)
N, N̂ = length(t), length(t̂)
M = _mapping_matrix(t̂, t)
Wls½ = LA.diagm(sqrt.(wls))
Wr½ = Array{Float64}(LA.I, N̂ - d, N̂ - d)
û, λ, Aitp = _reg_smooth_solve(
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
extrapolation_right, cache_parameters)
RegularizationSmooth(u,
û,
t,
Expand All @@ -185,16 +207,22 @@ A = RegularizationSmooth(
"""
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
wls::AbstractVector, d::Int = 2; λ::Real = 1.0,
alg::Symbol = :gcv_svd, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
alg::Symbol = :gcv_svd,
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
cache_parameters::Bool = false)
extrapolation_left, extrapolation_right = munge_extrapolation(
extrapolation, extrapolation_left, extrapolation_right)
u, t = munge_data(u, t)
t̂ = t
N = length(t)
M = Array{Float64}(LA.I, N, N)
Wls½ = LA.diagm(sqrt.(wls))
Wr½ = Array{Float64}(LA.I, N - d, N - d)
û, λ, Aitp = _reg_smooth_solve(
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
extrapolation_right, cache_parameters)
RegularizationSmooth(u,
û,
t,
Expand All @@ -219,16 +247,21 @@ A = RegularizationSmooth(
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
wls::AbstractVector, wr::AbstractVector, d::Int = 2;
λ::Real = 1.0, alg::Symbol = :gcv_svd,
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
cache_parameters::Bool = false)
extrapolation_left, extrapolation_right = munge_extrapolation(
extrapolation, extrapolation_left, extrapolation_right)
u, t = munge_data(u, t)
t̂ = t
N = length(t)
M = Array{Float64}(LA.I, N, N)
Wls½ = LA.diagm(sqrt.(wls))
Wr½ = LA.diagm(sqrt.(wr))
û, λ, Aitp = _reg_smooth_solve(
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
extrapolation_right, cache_parameters)
RegularizationSmooth(u,
û,
t,
Expand All @@ -252,8 +285,12 @@ A = RegularizationSmooth(
"""
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
wls::Symbol, d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd,
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
cache_parameters::Bool = false)
extrapolation_left, extrapolation_right = munge_extrapolation(
extrapolation, extrapolation_left, extrapolation_right)
u, t = munge_data(u, t)
t̂ = t
N = length(t)
Expand All @@ -262,7 +299,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing
Wls½ = LA.diagm(sqrt.(wls))
Wr½ = LA.diagm(sqrt.(wr))
û, λ, Aitp = _reg_smooth_solve(
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
extrapolation_right, cache_parameters)
RegularizationSmooth(u,
û,
t,
Expand All @@ -286,7 +324,9 @@ Solve for the smoothed dependent variables and create spline interpolator
function _reg_smooth_solve(
u::AbstractVector, t̂::AbstractVector, d::Int, M::AbstractMatrix,
Wls½::AbstractMatrix, Wr½::AbstractMatrix, λ::Real, alg::Symbol,
extrapolation_left::ExtrapolationType.T, extrapolation_right::ExtrapolationType.T)
extrapolation_left::ExtrapolationType.T,
extrapolation_right::ExtrapolationType.T,
cache_parameters::Bool)
λ = float(λ) # `float` expected by RT
D = _derivative_matrix(t̂, d)
Ψ = RT.setupRegularizationProblem(Wls½ * M, Wr½ * D)
Expand All @@ -303,7 +343,7 @@ function _reg_smooth_solve(
û = result.x
λ = result.λ
end
Aitp = CubicSpline(û, t̂; extrapolation_left, extrapolation_right)
Aitp = CubicSpline(û, t̂; extrapolation_left, extrapolation_right, cache_parameters)
# It seems logical to use B-Spline of order d+1, but I am unsure if theory supports the
# extra computational cost, JJS 12/25/21
#Aitp = BSplineInterpolation(û,t̂,d+1,:ArcLen,:Average)
Expand Down
4 changes: 3 additions & 1 deletion src/integrals.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,9 @@ function integral(A::AbstractInterpolation, t1::Number, t2::Number)

# Complete intervals
if A.cache_parameters
total += A.I[idx2 - 1]
if idx2 > 1
total += A.I[idx2 - 1]
end
if idx1 > 0
total -= A.I[idx1]
end
Expand Down
10 changes: 8 additions & 2 deletions test/integral_tests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,7 @@ using RegularizationTools
using StableRNGs

function test_integral(method; args = [], kwargs = [], name::String)
func = method(args...; kwargs..., extrapolation_left = ExtrapolationType.Extension,
extrapolation_right = ExtrapolationType.Extension)
func = method(args...; kwargs..., extrapolation = ExtrapolationType.Extension)
(; t) = func
t1 = minimum(t)
t2 = maximum(t)
Expand Down Expand Up @@ -63,6 +62,13 @@ function test_integral(method; args = [], kwargs = [], name::String)
@test_throws DataInterpolations.LeftExtrapolationError integral(func, t[1] - 1.0, t[2])
@test_throws DataInterpolations.RightExtrapolationError integral(
func, t[1], t[end] + 1.0)

# Test integration with cached parameters
func = method(args...; kwargs..., cache_parameters = true,
extrapolation = ExtrapolationType.Extension)
qint, err = quadgk(func, t1 - 1, t1; atol = 1e-12, rtol = 1e-12)
aint = integral(func, t1 - 1, t1)
@test isapprox(qint, aint, atol = 1e-6, rtol = 1e-8)
end

@testset "LinearInterpolation" begin
Expand Down
Loading