From a4c90d0c97bcddce1be2b985ef7508210963aa81 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Fri, 10 Nov 2023 15:13:20 +0000 Subject: [PATCH 1/8] test: add tests for integrals for Curvefit, RegularizationSmooth and LagrangeInterpolation --- test/integral_tests.jl | 89 ++++++++++++++++++++++++++----------- test/interpolation_tests.jl | 3 +- 2 files changed, 64 insertions(+), 28 deletions(-) diff --git a/test/integral_tests.jl b/test/integral_tests.jl index d781cb6d..e34738ce 100644 --- a/test/integral_tests.jl +++ b/test/integral_tests.jl @@ -1,6 +1,8 @@ using DataInterpolations, Test using QuadGK using DataInterpolations: integral +using Optim, ForwardDiff +using RegularizationTools function test_integral(func, tspan, name::String) t1 = minimum(tspan) @@ -12,22 +14,18 @@ function test_integral(func, tspan, name::String) end end -# Linear Interpolation -u = 2.0collect(1:10) -t = 1.0collect(1:10) -A = LinearInterpolation(u, t) - -test_integral(A, t, "Linear Interpolation (Vector)") - -# Quadratic Interpolation -u = [1.0, 4.0, 9.0, 16.0] -t = [1.0, 2.0, 3.0, 4.0] -A = QuadraticInterpolation(u, t) - -test_integral(A, t, "Quadratic Interpolation (Vector)") +@testset "LinearInterpolation" begin + u = 2.0collect(1:10) + t = 1.0collect(1:10) + A = LinearInterpolation(u, t) + test_integral(A, t, "Linear Interpolation (Vector)") +end -# Quadratic forward/backward Interpolation -@testset "QuadraticInterpolation - forward/backward modes" begin +@testset "QuadraticInterpolation" begin + u = [1.0, 4.0, 9.0, 16.0] + t = [1.0, 2.0, 3.0, 4.0] + A = QuadraticInterpolation(u, t) + test_integral(A, t, "Quadratic Interpolation (Vector)") u = [3.0, 0.0, 3.0, 0.0] t = [1.0, 2.0, 3.0, 4.0] A_f = QuadraticInterpolation(u, t) @@ -46,18 +44,57 @@ test_integral(A, t, "Quadratic Interpolation (Vector)") @test integral(A_b, 1.0, 4.0) ≈ 4.0 end -# QuadraticSpline Interpolation -u = [0.0, 1.0, 3.0] -t = [-1.0, 0.0, 1.0] - -A = QuadraticSpline(u, t) +@testset "LagrangeInterpolation" begin + u = [1.0, 4.0, 9.0] + t = [1.0, 2.0, 3.0] + A = LagrangeInterpolation(u, t) + @test_throws DataInterpolations.IntegralNotFoundError integral(A, 1.0, 2.0) + @test_throws DataInterpolations.IntegralNotFoundError integral(A, 5.0) +end -test_integral(A, t, "Quadratic Spline (Vector)") +@testset "QuadraticSpline" begin + u = [0.0, 1.0, 3.0] + t = [-1.0, 0.0, 1.0] + A = QuadraticSpline(u, t) + test_integral(A, t, "Quadratic Spline (Vector)") +end -# CubicSpline Interpolation -u = [0.0, 1.0, 3.0] -t = [-1.0, 0.0, 1.0] +@testset "CubicSpline" begin + u = [0.0, 1.0, 3.0] + t = [-1.0, 0.0, 1.0] + A = CubicSpline(u, t) + test_integral(A, t, "Cubic Spline Interpolation (Vector)") +end -A = CubicSpline(u, t) +@testset "RegularizationSmooth" begin + npts = 50 + xmin = 0.0 + xspan = 3 / 2 * π + x = collect(range(xmin, xmin + xspan, length = npts)) + rng = StableRNG(655) + x = x + xspan / npts * (rand(rng, npts) .- 0.5) + # select a subset randomly + idx = unique(rand(rng, collect(eachindex(x)), 20)) + t = x[unique(idx)] + npts = length(t) + ut = sin.(t) + stdev = 1e-1 * maximum(ut) + u = ut + stdev * randn(rng, npts) + # data must be ordered if t̂ is not provided + idx = sortperm(t) + tₒ = t[idx] + uₒ = u[idx] + A = RegularizationSmooth(uₒ, tₒ; alg = :fixed) + test_integral(A, tₒ, "RegularizationSmooth") +end -test_integral(A, t, "Cubic Spline Interpolation (Vector)") +@testset "Curvefit" begin + rng = StableRNG(12345) + model(x, p) = @. p[1] / (1 + exp(x - p[2])) + t = range(-10, stop = 10, length = 40) + u = model(t, [1.0, 2.0]) + 0.01 * randn(rng, length(t)) + p0 = [0.5, 0.5] + A = Curvefit(u, t, model, p0, LBFGS()) + @test_throws DataInterpolations.IntegralNotFoundError integral(A, 0.0, 1.0) + @test_throws DataInterpolations.IntegralNotFoundError integral(A, 5.0) +end diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 72b5657e..a10c1bde 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -1,7 +1,6 @@ using DataInterpolations, Test using StableRNGs -using Optim -import ForwardDiff +using Optim, ForwardDiff @testset "Linear Interpolation" begin for t in (1.0:10.0, 1.0collect(1:10)) From 223cbf660afca5ee2f682814e556645e031b818d Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Fri, 10 Nov 2023 15:14:50 +0000 Subject: [PATCH 2/8] feat: add integral methods for Curvefit, RegularizationSmooth and LagrangeInterpolation --- ext/DataInterpolationsOptimExt.jl | 15 ++++++++++++--- ext/DataInterpolationsRegularizationToolsExt.jl | 16 +++++++++++++--- src/integrals.jl | 3 +++ 3 files changed, 28 insertions(+), 6 deletions(-) diff --git a/ext/DataInterpolationsOptimExt.jl b/ext/DataInterpolationsOptimExt.jl index 76ad8eb0..f2643216 100644 --- a/ext/DataInterpolationsOptimExt.jl +++ b/ext/DataInterpolationsOptimExt.jl @@ -2,7 +2,8 @@ module DataInterpolationsOptimExt using DataInterpolations import DataInterpolations: munge_data, - Curvefit, CurvefitCache, _interpolate, get_show, derivative, ExtrapolationError + Curvefit, CurvefitCache, _interpolate, get_show, derivative, ExtrapolationError, + integral, IntegralNotFoundError isdefined(Base, :get_extension) ? (using Optim, ForwardDiff) : (using ..Optim, ..ForwardDiff) @@ -52,9 +53,17 @@ function derivative(A::CurvefitCache{<:AbstractVector{<:Number}}, ForwardDiff.derivative(x -> A.m(x, A.pmin), t) end -function get_show(interp::CurvefitCache) +function get_show(A::CurvefitCache) return "Curvefit" * - " with $(length(interp.t)) points, using $(nameof(typeof(interp.alg)))\n" + " with $(length(A.t)) points, using $(nameof(typeof(A.alg)))\n" +end + +function integral(A::CurvefitCache{<:AbstractVector{<:Number}}, t::Number) + throw(IntegralNotFoundError()) +end + +function integral(A::CurvefitCache{<:AbstractVector{<:Number}}, t1::Number, t2::Number) + throw(IntegralNotFoundError()) end end # module diff --git a/ext/DataInterpolationsRegularizationToolsExt.jl b/ext/DataInterpolationsRegularizationToolsExt.jl index dae1f606..bb34f574 100644 --- a/ext/DataInterpolationsRegularizationToolsExt.jl +++ b/ext/DataInterpolationsRegularizationToolsExt.jl @@ -2,7 +2,7 @@ module DataInterpolationsRegularizationToolsExt using DataInterpolations import DataInterpolations: munge_data, - _interpolate, RegularizationSmooth, get_show, derivative + _interpolate, RegularizationSmooth, get_show, derivative, integral using LinearAlgebra isdefined(Base, :get_extension) ? (import RegularizationTools as RT) : @@ -339,9 +339,19 @@ function derivative(A::RegularizationSmooth{ derivative(A.Aitp, t) end -function get_show(interp::RegularizationSmooth) +function get_show(A::RegularizationSmooth) return "RegularizationSmooth" * - " with $(length(interp.t)) points, with regularization coefficient $(interp.λ)\n" + " with $(length(A.t)) points, with regularization coefficient $(A.λ)\n" +end + +function integral(A::RegularizationSmooth{<:AbstractVector{<:Number}}, t::Number) + integral(A.Aitp, t) +end + +function integral(A::RegularizationSmooth{<:AbstractVector{<:Number}}, + t1::Number, + t2::Number) + integral(A.Aitp, t1, t2) end end # module diff --git a/src/integrals.jl b/src/integrals.jl index d8946656..e1f663bd 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -96,3 +96,6 @@ function _integral(A::CubicSpline{<:AbstractVector{<:Number}}, idx::Number, t::N (h2^2 * t1 * z2 - h2^2 * t2 * z1 - t1^3 * z2 - 6 * t1 * u2 + t2^3 * z1 + 6 * t2 * u1) / (6 * h2)) end + +integral(A::LagrangeInterpolation, t1::Number, t2::Number) = throw(IntegralNotFoundError()) +integral(A::LagrangeInterpolation, t::Number) = throw(IntegralNotFoundError()) From ba31e2157e4cb46b00c9bba1e592949dc7a2cab1 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Fri, 10 Nov 2023 15:15:20 +0000 Subject: [PATCH 3/8] refactor: throw an error message when integral is not defined --- src/DataInterpolations.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl index 1d59e0a1..9aaad786 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -44,6 +44,12 @@ function Base.showerror(io::IO, e::ExtrapolationError) print(io, EXTRAPOLATION_ERROR) end +const INTEGRAL_NOT_FOUND_ERROR = "Cannot integrate it analytically. Please use Numerical Integration methods." +struct IntegralNotFoundError <: Exception end +function Base.showerror(io::IO, e::IntegralNotFoundError) + print(io, INTEGRAL_NOT_FOUND_ERROR) +end + export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline, BSplineInterpolation, BSplineApprox From 4e4d32478fa6d95dcfe71d6592e29ca2543fa015 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Sat, 18 Nov 2023 14:51:44 +0000 Subject: [PATCH 4/8] feat: add integral for AkimaInterpolation --- src/integrals.jl | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/integrals.jl b/src/integrals.jl index e1f663bd..1b91349f 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -1,10 +1,11 @@ function integral(A::AbstractInterpolation, t::Number) - bw, fw = samples(A) - idx = max(1 + bw, min(searchsortedlast(A.t, t), length(A.t) - fw)) - _integral(A, idx, t) + ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) + integral(A, A.t[1], t) 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()) bw, fw = samples(A) # the index less than or equal to t1 idx1 = max(1 + bw, min(searchsortedlast(A.t, t1), length(A.t) - fw)) @@ -97,5 +98,11 @@ function _integral(A::CubicSpline{<:AbstractVector{<:Number}}, idx::Number, t::N (6 * h2)) end +samples(A::AkimaInterpolation{<:AbstractVector{<:Number}}) = (0, 1) +function _integral(A::AkimaInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number) + t1 = A.t[idx] + A.u[idx]*(t - t1) + A.b[idx]*((t - t1)^2/2) + A.c[idx]*((t - t1)^3/3) + A.d[idx]*((t - t1)^4/4) +end + integral(A::LagrangeInterpolation, t1::Number, t2::Number) = throw(IntegralNotFoundError()) integral(A::LagrangeInterpolation, t::Number) = throw(IntegralNotFoundError()) From 6c6fbfa3f94509707e7388a78dfe89f292bbf273 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Sat, 18 Nov 2023 14:52:09 +0000 Subject: [PATCH 5/8] test: refactor integral tests --- test/integral_tests.jl | 66 ++++++++++++++++++++++++------------------ 1 file changed, 38 insertions(+), 28 deletions(-) diff --git a/test/integral_tests.jl b/test/integral_tests.jl index e34738ce..366941fb 100644 --- a/test/integral_tests.jl +++ b/test/integral_tests.jl @@ -3,45 +3,52 @@ using QuadGK using DataInterpolations: integral using Optim, ForwardDiff using RegularizationTools +using StableRNGs -function test_integral(func, tspan, name::String) - t1 = minimum(tspan) - t2 = maximum(tspan) +function test_integral(method, u, t; args = [], kwargs = [], name::String) + func = method(u, t, args...; kwargs..., extrapolate = true) + t1 = minimum(t) + t2 = maximum(t) @testset "$name" begin - qint, err = quadgk(func, t1, t2) + # integral(A, t1, t2) + qint, err = quadgk(func, t1, t2; atol = 1e-12, rtol = 1e-12) aint = integral(func, t1, t2) @test isapprox(qint, aint, atol = 1e-8) + + # integral(A, t) + qint, err = quadgk(func, t1, (t1 + t2)/2; atol = 1e-12, rtol = 1e-12) + aint = integral(func, (t1 + t2)/2) + @test isapprox(qint, aint, atol = 1e-8) + + # integrals with extrapolation + qint, err = quadgk(func, t1-5.0, (t1 + t2)/2; atol = 1e-12, rtol = 1e-12) + aint = integral(func, t1-5.0, (t1 + t2)/2) + @test isapprox(qint, aint, atol = 1e-8) + + qint, err = quadgk(func, (t1 + t2)/2, t2+5.0; atol = 1e-12, rtol = 1e-12) + aint = integral(func, (t1 + t2)/2, t2+5.0) + @test isapprox(qint, aint, atol = 1e-8) end + func = method(u, t, args...; kwargs...) + @test_throws DataInterpolations.ExtrapolationError integral(func, t[1] - 1.0) + @test_throws DataInterpolations.ExtrapolationError integral(func, t[end] + 1.0) + @test_throws DataInterpolations.ExtrapolationError integral(func, t[1] - 1.0, t[2]) + @test_throws DataInterpolations.ExtrapolationError integral(func, t[1], t[end] + 1.0) end @testset "LinearInterpolation" begin u = 2.0collect(1:10) t = 1.0collect(1:10) - A = LinearInterpolation(u, t) - test_integral(A, t, "Linear Interpolation (Vector)") + test_integral(LinearInterpolation, u, t; name = "Linear Interpolation (Vector)") end @testset "QuadraticInterpolation" begin u = [1.0, 4.0, 9.0, 16.0] t = [1.0, 2.0, 3.0, 4.0] - A = QuadraticInterpolation(u, t) - test_integral(A, t, "Quadratic Interpolation (Vector)") + test_integral(QuadraticInterpolation, u, t; name = "Quadratic Interpolation (Vector)") u = [3.0, 0.0, 3.0, 0.0] t = [1.0, 2.0, 3.0, 4.0] - A_f = QuadraticInterpolation(u, t) - A_b = QuadraticInterpolation(u, t, :Backward) - @test integral(A_f, 1.0, 2.0) ≈ 1.0 - @test integral(A_f, 2.0, 3.0) ≈ 2.0 - @test integral(A_f, 3.0, 4.0) ≈ 2.0 - @test integral(A_f, 1.0, 3.0) ≈ 3.0 - @test integral(A_f, 2.0, 4.0) ≈ 4.0 - @test integral(A_f, 1.0, 4.0) ≈ 5.0 - @test integral(A_b, 1.0, 2.0) ≈ 1.0 - @test integral(A_b, 2.0, 3.0) ≈ 1.0 - @test integral(A_b, 3.0, 4.0) ≈ 2.0 - @test integral(A_b, 1.0, 3.0) ≈ 2.0 - @test integral(A_b, 2.0, 4.0) ≈ 3.0 - @test integral(A_b, 1.0, 4.0) ≈ 4.0 + test_integral(QuadraticInterpolation, u, t; args = [:Backward], name = "Quadratic Interpolation (Vector)") end @testset "LagrangeInterpolation" begin @@ -55,15 +62,19 @@ end @testset "QuadraticSpline" begin u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] - A = QuadraticSpline(u, t) - test_integral(A, t, "Quadratic Spline (Vector)") + test_integral(QuadraticSpline, u, t; name = "Quadratic Spline (Vector)") end @testset "CubicSpline" begin u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] - A = CubicSpline(u, t) - test_integral(A, t, "Cubic Spline Interpolation (Vector)") + test_integral(CubicSpline, u, t; name = "Cubic Spline Interpolation (Vector)") +end + +@testset "AkimaInterpolation" 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) + test_integral(AkimaInterpolation, u, t; name = "Akima Interpolation (Vector)") end @testset "RegularizationSmooth" begin @@ -84,8 +95,7 @@ end idx = sortperm(t) tₒ = t[idx] uₒ = u[idx] - A = RegularizationSmooth(uₒ, tₒ; alg = :fixed) - test_integral(A, tₒ, "RegularizationSmooth") + test_integral(RegularizationSmooth, uₒ, tₒ; kwargs = [:alg => :fixed], name = "RegularizationSmooth") end @testset "Curvefit" begin From 0ced7ac93b6990f80b96141af4a925cb53fb4d88 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Sat, 18 Nov 2023 14:53:03 +0000 Subject: [PATCH 6/8] refactor: extrapolate the last spline in AkimaInterpolation instead of constant --- src/derivatives.jl | 7 ++++--- src/interpolation_methods.jl | 4 +--- test/interpolation_tests.jl | 4 ++-- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index 46671964..84d5b534 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -107,9 +107,10 @@ derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number, i) = derivativ derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, i) = derivative(A, t), i function derivative(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess) - t < A.t[1] && return zero(A.u[1]), 1 - t > A.t[end] && return zero(A.u[end]), lastindex(t) - i = searchsortedlastcorrelated(A.t, t, iguess) + i = searchsortedfirstcorrelated(A.t, t, iguess) + i > length(A.t) ? i -= 1 : nothing + i -= 1 + i == 0 ? i += 1 : nothing j = min(i, length(A.c)) # for smooth derivative at A.t[end] wj = t - A.t[i] (@evalpoly wj A.b[i] 2A.c[j] 3A.d[j]), i diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 328c9208..334ce030 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -123,9 +123,7 @@ function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, i) end function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess) - i = searchsortedlastcorrelated(A.t, t, iguess) - i == 0 && return A.u[1], i - i == length(A.t) && return A.u[end], i + i = max(1, min(searchsortedlastcorrelated(A.t, t, iguess), length(A.t) - 1)) wj = t - A.t[i] (@evalpoly wj A.u[i] A.b[i] A.c[i] A.d[i]), i end diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index a10c1bde..c8416e14 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -324,8 +324,8 @@ end # Test extrapolation A = AkimaInterpolation(u, t; extrapolate = true) - @test A(-1.0) == 0.0 - @test A(11.0) == 3.0 + @test A(-1.0) ≈ -5.0 + @test A(11.0) ≈ -3.924742268041234 A = AkimaInterpolation(u, t) @test_throws DataInterpolations.ExtrapolationError A(-1.0) @test_throws DataInterpolations.ExtrapolationError A(11.0) From cf76d19a1f3e4f72c2fdebbbffae6e08e86dd434 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Fri, 24 Nov 2023 14:13:45 +0000 Subject: [PATCH 7/8] refactor: add integrals for BSpline methods to error out --- src/integrals.jl | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/src/integrals.jl b/src/integrals.jl index 1b91349f..0e796d4d 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -99,10 +99,29 @@ function _integral(A::CubicSpline{<:AbstractVector{<:Number}}, idx::Number, t::N end samples(A::AkimaInterpolation{<:AbstractVector{<:Number}}) = (0, 1) -function _integral(A::AkimaInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number) +function _integral(A::AkimaInterpolation{<:AbstractVector{<:Number}}, + idx::Number, + t::Number) t1 = A.t[idx] - A.u[idx]*(t - t1) + A.b[idx]*((t - t1)^2/2) + A.c[idx]*((t - t1)^3/3) + A.d[idx]*((t - t1)^4/4) + A.u[idx] * (t - t1) + A.b[idx] * ((t - t1)^2 / 2) + A.c[idx] * ((t - t1)^3 / 3) + + 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 From 1067954f3c4f1471eef75778e7ab8c5c215a4c77 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Fri, 24 Nov 2023 14:14:03 +0000 Subject: [PATCH 8/8] test: update integral tests --- test/integral_tests.jl | 40 ++++++++++++++++++++++++++++++++-------- 1 file changed, 32 insertions(+), 8 deletions(-) diff --git a/test/integral_tests.jl b/test/integral_tests.jl index 366941fb..cb973570 100644 --- a/test/integral_tests.jl +++ b/test/integral_tests.jl @@ -16,17 +16,17 @@ function test_integral(method, u, t; args = [], kwargs = [], name::String) @test isapprox(qint, aint, atol = 1e-8) # integral(A, t) - qint, err = quadgk(func, t1, (t1 + t2)/2; atol = 1e-12, rtol = 1e-12) - aint = integral(func, (t1 + t2)/2) + qint, err = quadgk(func, t1, (t1 + t2) / 2; atol = 1e-12, rtol = 1e-12) + aint = integral(func, (t1 + t2) / 2) @test isapprox(qint, aint, atol = 1e-8) # integrals with extrapolation - qint, err = quadgk(func, t1-5.0, (t1 + t2)/2; atol = 1e-12, rtol = 1e-12) - aint = integral(func, t1-5.0, (t1 + t2)/2) + qint, err = quadgk(func, t1 - 5.0, (t1 + t2) / 2; atol = 1e-12, rtol = 1e-12) + aint = integral(func, t1 - 5.0, (t1 + t2) / 2) @test isapprox(qint, aint, atol = 1e-8) - qint, err = quadgk(func, (t1 + t2)/2, t2+5.0; atol = 1e-12, rtol = 1e-12) - aint = integral(func, (t1 + t2)/2, t2+5.0) + qint, err = quadgk(func, (t1 + t2) / 2, t2 + 5.0; atol = 1e-12, rtol = 1e-12) + aint = integral(func, (t1 + t2) / 2, t2 + 5.0) @test isapprox(qint, aint, atol = 1e-8) end func = method(u, t, args...; kwargs...) @@ -48,7 +48,11 @@ end test_integral(QuadraticInterpolation, u, t; name = "Quadratic Interpolation (Vector)") u = [3.0, 0.0, 3.0, 0.0] t = [1.0, 2.0, 3.0, 4.0] - test_integral(QuadraticInterpolation, u, t; args = [:Backward], name = "Quadratic Interpolation (Vector)") + test_integral(QuadraticInterpolation, + u, + t; + args = [:Backward], + name = "Quadratic Interpolation (Vector)") end @testset "LagrangeInterpolation" begin @@ -95,7 +99,11 @@ end idx = sortperm(t) tₒ = t[idx] uₒ = u[idx] - test_integral(RegularizationSmooth, uₒ, tₒ; kwargs = [:alg => :fixed], name = "RegularizationSmooth") + test_integral(RegularizationSmooth, + uₒ, + tₒ; + kwargs = [:alg => :fixed], + name = "RegularizationSmooth") end @testset "Curvefit" begin @@ -108,3 +116,19 @@ end @test_throws DataInterpolations.IntegralNotFoundError integral(A, 0.0, 1.0) @test_throws DataInterpolations.IntegralNotFoundError integral(A, 5.0) end + +@testset "BSplineInterpolation" begin + 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] + A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform) + @test_throws DataInterpolations.IntegralNotFoundError integral(A, 1.0, 100.0) + @test_throws DataInterpolations.IntegralNotFoundError integral(A, 50.0) +end + +@testset "BSplineApprox" begin + 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] + A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform) + @test_throws DataInterpolations.IntegralNotFoundError integral(A, 1.0, 100.0) + @test_throws DataInterpolations.IntegralNotFoundError integral(A, 50.0) +end