Add integrals for remaining methods

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

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

src/DataInterpolations.jl

@@ -44,6 +44,12 @@
     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

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

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))
@@ -96,3 +97,31 @@ 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
+
+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())
+
+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

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

test/integral_tests.jl

@@ -1,63 +1,134 @@
 using DataInterpolations, Test
 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)
-    end
-end
 
-# Linear Interpolation
-u = 2.0collect(1:10)
-t = 1.0collect(1:10)
-A = LinearInterpolation(u, t)
+        # 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)
 
-test_integral(A, t, "Linear Interpolation (Vector)")
+        # 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)
 
-# Quadratic Interpolation
-u = [1.0, 4.0, 9.0, 16.0]
-t = [1.0, 2.0, 3.0, 4.0]
-A = QuadraticInterpolation(u, t)
+        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
 
-test_integral(A, t, "Quadratic Interpolation (Vector)")
+@testset "LinearInterpolation" begin
+    u = 2.0collect(1:10)
+    t = 1.0collect(1:10)
+    test_integral(LinearInterpolation, u, t; name = "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]
+    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
 
-# QuadraticSpline Interpolation
-u = [0.0, 1.0, 3.0]
-t = [-1.0, 0.0, 1.0]
+@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
 
-A = QuadraticSpline(u, t)
+@testset "QuadraticSpline" begin
+    u = [0.0, 1.0, 3.0]
+    t = [-1.0, 0.0, 1.0]
+    test_integral(QuadraticSpline, u, t; name = "Quadratic Spline (Vector)")
+end
 
-test_integral(A, t, "Quadratic Spline (Vector)")
+@testset "CubicSpline" begin
+    u = [0.0, 1.0, 3.0]
+    t = [-1.0, 0.0, 1.0]
+    test_integral(CubicSpline, u, t; name = "Cubic Spline Interpolation (Vector)")
+end
 
-# CubicSpline Interpolation
-u = [0.0, 1.0, 3.0]
-t = [-1.0, 0.0, 1.0]
+@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
 
-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]
+    test_integral(RegularizationSmooth,
+        uₒ,
+        tₒ;
+        kwargs = [:alg => :fixed],
+        name = "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
+
+@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

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))
@@ -325,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)