Skip to content

Commit 4a337ac

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #391 from SouthEndMusic/fix_integration_edgecase
Fix integrating with `cache_parameters = true`
2 parents b74be2c + 7985ae3 commit 4a337ac

File tree

3 files changed

+71
-23
lines changed

3 files changed

+71
-23
lines changed

‎ext/DataInterpolationsRegularizationToolsExt.jl

Lines changed: 60 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
module DataInterpolationsRegularizationToolsExt
22

33
using DataInterpolations
4-
import DataInterpolations: munge_data,
4+
import DataInterpolations: munge_data, munge_extrapolation,
55
_interpolate, RegularizationSmooth, get_show, derivative,
66
integral
77
using LinearAlgebra
@@ -70,14 +70,19 @@ A = RegularizationSmooth(u, t, t̂, wls, wr, d; λ = 1.0, alg = :gcv_svd)
7070
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
7171
wls::AbstractVector, wr::AbstractVector, d::Int = 2;
7272
λ::Real = 1.0, alg::Symbol = :gcv_svd,
73+
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
7374
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
74-
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
75+
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
76+
cache_parameters::Bool = false)
77+
extrapolation_left, extrapolation_right = munge_extrapolation(
78+
extrapolation, extrapolation_left, extrapolation_right)
7579
u, t = munge_data(u, t)
7680
M = _mapping_matrix(t̂, t)
7781
Wls½ = LA.diagm(sqrt.(wls))
7882
Wr½ = LA.diagm(sqrt.(wr))
7983
û, λ, Aitp = _reg_smooth_solve(
80-
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
84+
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
85+
extrapolation_right, cache_parameters)
8186
RegularizationSmooth(
8287
u, û, t, t̂, wls, wr, d, λ, alg, Aitp, extrapolation_left, extrapolation_right)
8388
end
@@ -90,16 +95,22 @@ A = RegularizationSmooth(u, t, d; λ = 1.0, alg = :gcv_svd, extrapolate = false)
9095
"""
9196
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, d::Int = 2;
9297
λ::Real = 1.0,
93-
alg::Symbol = :gcv_svd, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
94-
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
98+
alg::Symbol = :gcv_svd,
99+
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
100+
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
101+
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
102+
cache_parameters::Bool = false)
103+
extrapolation_left, extrapolation_right = munge_extrapolation(
104+
extrapolation, extrapolation_left, extrapolation_right)
95105
u, t = munge_data(u, t)
96106
t̂ = t
97107
N = length(t)
98108
M = Array{Float64}(LA.I, N, N)
99109
Wls½ = Array{Float64}(LA.I, N, N)
100110
Wr½ = Array{Float64}(LA.I, N - d, N - d)
101111
û, λ, Aitp = _reg_smooth_solve(
102-
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
112+
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
113+
extrapolation_right, cache_parameters)
103114
RegularizationSmooth(u,
104115
û,
105116
t,
@@ -122,15 +133,20 @@ A = RegularizationSmooth(u, t, t̂, d; λ = 1.0, alg = :gcv_svd, extrapolate = f
122133
"""
123134
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
124135
d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd,
136+
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
125137
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
126-
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
138+
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
139+
cache_parameters::Bool = false)
140+
extrapolation_left, extrapolation_right = munge_extrapolation(
141+
extrapolation, extrapolation_left, extrapolation_right)
127142
u, t = munge_data(u, t)
128143
N, N̂ = length(t), length(t̂)
129144
M = _mapping_matrix(t̂, t)
130145
Wls½ = Array{Float64}(LA.I, N, N)
131146
Wr½ = Array{Float64}(LA.I, N̂ - d, N̂ - d)
132147
û, λ, Aitp = _reg_smooth_solve(
133-
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
148+
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
149+
extrapolation_right, cache_parameters)
134150
RegularizationSmooth(u,
135151
û,
136152
t,
@@ -153,15 +169,21 @@ A = RegularizationSmooth(u, t, t̂, wls, d; λ = 1.0, alg = :gcv_svd, extrapolat
153169
"""
154170
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
155171
wls::AbstractVector, d::Int = 2; λ::Real = 1.0,
156-
alg::Symbol = :gcv_svd, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
157-
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
172+
alg::Symbol = :gcv_svd,
173+
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
174+
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
175+
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
176+
cache_parameters::Bool = false)
177+
extrapolation_left, extrapolation_right = munge_extrapolation(
178+
extrapolation, extrapolation_left, extrapolation_right)
158179
u, t = munge_data(u, t)
159180
N, N̂ = length(t), length(t̂)
160181
M = _mapping_matrix(t̂, t)
161182
Wls½ = LA.diagm(sqrt.(wls))
162183
Wr½ = Array{Float64}(LA.I, N̂ - d, N̂ - d)
163184
û, λ, Aitp = _reg_smooth_solve(
164-
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
185+
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
186+
extrapolation_right, cache_parameters)
165187
RegularizationSmooth(u,
166188
û,
167189
t,
@@ -185,16 +207,22 @@ A = RegularizationSmooth(
185207
"""
186208
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
187209
wls::AbstractVector, d::Int = 2; λ::Real = 1.0,
188-
alg::Symbol = :gcv_svd, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
189-
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
210+
alg::Symbol = :gcv_svd,
211+
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
212+
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
213+
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
214+
cache_parameters::Bool = false)
215+
extrapolation_left, extrapolation_right = munge_extrapolation(
216+
extrapolation, extrapolation_left, extrapolation_right)
190217
u, t = munge_data(u, t)
191218
t̂ = t
192219
N = length(t)
193220
M = Array{Float64}(LA.I, N, N)
194221
Wls½ = LA.diagm(sqrt.(wls))
195222
Wr½ = Array{Float64}(LA.I, N - d, N - d)
196223
û, λ, Aitp = _reg_smooth_solve(
197-
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
224+
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
225+
extrapolation_right, cache_parameters)
198226
RegularizationSmooth(u,
199227
û,
200228
t,
@@ -219,16 +247,21 @@ A = RegularizationSmooth(
219247
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
220248
wls::AbstractVector, wr::AbstractVector, d::Int = 2;
221249
λ::Real = 1.0, alg::Symbol = :gcv_svd,
250+
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
222251
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
223-
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
252+
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
253+
cache_parameters::Bool = false)
254+
extrapolation_left, extrapolation_right = munge_extrapolation(
255+
extrapolation, extrapolation_left, extrapolation_right)
224256
u, t = munge_data(u, t)
225257
t̂ = t
226258
N = length(t)
227259
M = Array{Float64}(LA.I, N, N)
228260
Wls½ = LA.diagm(sqrt.(wls))
229261
Wr½ = LA.diagm(sqrt.(wr))
230262
û, λ, Aitp = _reg_smooth_solve(
231-
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
263+
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
264+
extrapolation_right, cache_parameters)
232265
RegularizationSmooth(u,
233266
û,
234267
t,
@@ -252,8 +285,12 @@ A = RegularizationSmooth(
252285
"""
253286
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
254287
wls::Symbol, d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd,
288+
extrapolation::ExtrapolationType.T = ExtrapolationType.None,
255289
extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
256-
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
290+
extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
291+
cache_parameters::Bool = false)
292+
extrapolation_left, extrapolation_right = munge_extrapolation(
293+
extrapolation, extrapolation_left, extrapolation_right)
257294
u, t = munge_data(u, t)
258295
t̂ = t
259296
N = length(t)
@@ -262,7 +299,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing
262299
Wls½ = LA.diagm(sqrt.(wls))
263300
Wr½ = LA.diagm(sqrt.(wr))
264301
û, λ, Aitp = _reg_smooth_solve(
265-
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left, extrapolation_right)
302+
u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
303+
extrapolation_right, cache_parameters)
266304
RegularizationSmooth(u,
267305
û,
268306
t,
@@ -286,7 +324,9 @@ Solve for the smoothed dependent variables and create spline interpolator
286324
function _reg_smooth_solve(
287325
u::AbstractVector, t̂::AbstractVector, d::Int, M::AbstractMatrix,
288326
Wls½::AbstractMatrix, Wr½::AbstractMatrix, λ::Real, alg::Symbol,
289-
extrapolation_left::ExtrapolationType.T, extrapolation_right::ExtrapolationType.T)
327+
extrapolation_left::ExtrapolationType.T,
328+
extrapolation_right::ExtrapolationType.T,
329+
cache_parameters::Bool)
290330
λ = float(λ) # `float` expected by RT
291331
D = _derivative_matrix(t̂, d)
292332
Ψ = RT.setupRegularizationProblem(Wls½ * M, Wr½ * D)
@@ -303,7 +343,7 @@ function _reg_smooth_solve(
303343
û = result.x
304344
λ = result.λ
305345
end
306-
Aitp = CubicSpline(û, t̂; extrapolation_left, extrapolation_right)
346+
Aitp = CubicSpline(û, t̂; extrapolation_left, extrapolation_right, cache_parameters)
307347
# It seems logical to use B-Spline of order d+1, but I am unsure if theory supports the
308348
# extra computational cost, JJS 12/25/21
309349
#Aitp = BSplineInterpolation(û,t̂,d+1,:ArcLen,:Average)

‎src/integrals.jl

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,9 @@ function integral(A::AbstractInterpolation, t1::Number, t2::Number)
5252

5353
# Complete intervals
5454
if A.cache_parameters
55-
total += A.I[idx2 - 1]
55+
if idx2 > 1
56+
total += A.I[idx2 - 1]
57+
end
5658
if idx1 > 0
5759
total -= A.I[idx1]
5860
end

‎test/integral_tests.jl

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,8 +7,7 @@ using StableRNGs
77
using Unitful
88

99
function test_integral(method; args = [], kwargs = [], name::String)
10-
func = method(args...; kwargs..., extrapolation_left = ExtrapolationType.Extension,
11-
extrapolation_right = ExtrapolationType.Extension)
10+
func = method(args...; kwargs..., extrapolation = ExtrapolationType.Extension)
1211
(; t) = func
1312
t1 = minimum(t)
1413
t2 = maximum(t)
@@ -64,6 +63,13 @@ function test_integral(method; args = [], kwargs = [], name::String)
6463
@test_throws DataInterpolations.LeftExtrapolationError integral(func, t[1] - 1.0, t[2])
6564
@test_throws DataInterpolations.RightExtrapolationError integral(
6665
func, t[1], t[end] + 1.0)
66+
67+
# Test integration with cached parameters
68+
func = method(args...; kwargs..., cache_parameters = true,
69+
extrapolation = ExtrapolationType.Extension)
70+
qint, err = quadgk(func, t1 - 1, t1; atol = 1e-12, rtol = 1e-12)
71+
aint = integral(func, t1 - 1, t1)
72+
@test isapprox(qint, aint, atol = 1e-6, rtol = 1e-8)
6773
end
6874

6975
@testset "LinearInterpolation" begin

0 commit comments

Comments
 (0)