From 7d348ff2a1efb07984c5368e309d853eb4a34f85 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Sun, 5 Nov 2023 19:38:36 +0530 Subject: [PATCH 1/3] test: update derivative tests and make it more robust --- test/derivative_tests.jl | 40 ++++++++++++++++++++++++++++++++++------ 1 file changed, 34 insertions(+), 6 deletions(-) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index d735efb7..d1eefe1e 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -1,17 +1,45 @@ using DataInterpolations, Test using FiniteDifferences using DataInterpolations: derivative -# using Symbolics +using Symbolics function test_derivatives(method, u, t, args...; name::String) func = method(u, t, args...; extrapolate = true) - trange = range(minimum(t) - 5.0, maximum(t) + 5.0, length = 32) + trange = collect(range(minimum(t) - 5.0, maximum(t) + 5.0, step = 0.1)) + trange_exclude = filter(x -> !in(x, t), trange) @testset "$name" begin - for t in trange - cdiff = central_fdm(5, 1; geom = true)(_t -> func(_t), t) - adiff = derivative(func, t) + # Rest of the points + for _t in trange_exclude + cdiff = central_fdm(5, 1; geom = true)(func, _t) + adiff = derivative(func, _t) @test isapprox(cdiff, adiff, atol = 1e-8) end + + # Interpolation time points + for _t in t[2:end-1] + fdiff = if func isa BSplineInterpolation || func isa BSplineApprox + forward_fdm(5, 1; geom = true)(func, _t) + else + backward_fdm(5, 1; geom = true)(func, _t) + end + adiff = derivative(func, _t) + @test isapprox(fdiff, adiff, atol = 1e-8) + end + + # t = t0 + fdiff = forward_fdm(5, 1; geom = true)(func, t[1]) + adiff = derivative(func, t[1]) + if func isa BSplineInterpolation || func isa BSplineApprox + # Bug in BSplines + @test isapprox(fdiff, adiff, atol = 1e-8) broken = true + else + @test isapprox(fdiff, adiff, atol = 1e-8) + end + + # t = tend + fdiff = backward_fdm(5, 1; geom = true)(func, t[end]) + adiff = derivative(func, t[end]) + @test isapprox(fdiff, adiff, atol = 1e-8) end func = method(u, t, args...) @test_throws DataInterpolations.ExtrapolationError derivative(func, t[1] - 1.0) @@ -19,7 +47,7 @@ function test_derivatives(method, u, t, args...; name::String) end @testset "Linear Interpolation" begin - u = 2.0collect(1:10) + u = vcat(collect(1:5), 2*collect(6:10)) t = 1.0collect(1:10) test_derivatives(LinearInterpolation, u, t; name = "Linear Interpolation (Vector)") u = vcat(2.0collect(1:10)', 3.0collect(1:10)') From 2bd616669a0461199eb2433663a226649a10166d Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Sun, 5 Nov 2023 19:39:04 +0530 Subject: [PATCH 2/3] fix: derivative in Lagrange Interpolation at interpolation points --- src/derivatives.jl | 100 ++++++++++++++++++++++----------------------- 1 file changed, 48 insertions(+), 52 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index 85aa5264..518a9945 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -37,74 +37,70 @@ end function derivative(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 zero(A.u[idxs[1]]) - end - G = zero(A.u[1]) - F = zero(A.t[1]) - DG = zero(A.u[1]) - DF = zero(A.t[1]) - tmp = G - for i in 1:length(idxs) - if isnan(A.bcache[idxs[i]]) + der = zero(A.u[1]) + for j in eachindex(A.t) + tmp = zero(A.t[1]) + if isnan(A.bcache[j]) mult = one(A.t[1]) - for j in 1:(i - 1) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + for i in 1:(j - 1) + mult *= (A.t[j] - A.t[i]) end - for j in (i + 1):length(idxs) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + for i in (j + 1):length(A.t) + mult *= (A.t[j] - A.t[i]) end - A.bcache[idxs[i]] = mult + A.bcache[j] = mult else - mult = A.bcache[idxs[i]] + mult = A.bcache[j] + end + for l in eachindex(A.t) + if l != j + k = one(A.t[1]) + for m in eachindex(A.t) + if m != j && m != l + k *= (t - A.t[m]) + end + end + k *= inv(mult) + tmp += k + end end - wi = inv(mult) - tti = t - A.t[idxs[i]] - tmp = wi / (t - A.t[idxs[i]]) - g = tmp * A.u[idxs[i]] - G += g - DG -= g / (t - A.t[idxs[i]]) - F += tmp - DF -= tmp / (t - A.t[idxs[i]]) + der += A.u[j]*tmp end - (DG * F - G * DF) / (F^2) + der end function derivative(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 zero(A.u[:, idxs[1]]) - end - G = zero(A.u[:, 1]) - F = zero(A.t[1]) - DG = zero(A.u[:, 1]) - DF = zero(A.t[1]) - tmp = G - for i in 1:length(idxs) - if isnan(A.bcache[idxs[i]]) + der = zero(A.u[:, 1]) + for j in eachindex(A.t) + tmp = zero(A.t[1]) + if isnan(A.bcache[j]) mult = one(A.t[1]) - for j in 1:(i - 1) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + for i in 1:(j - 1) + mult *= (A.t[j] - A.t[i]) end - for j in (i + 1):length(idxs) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + for i in (j + 1):length(A.t) + mult *= (A.t[j] - A.t[i]) end - A.bcache[idxs[i]] = mult + A.bcache[j] = mult else - mult = A.bcache[idxs[i]] + mult = A.bcache[j] + end + for l in eachindex(A.t) + if l != j + k = one(A.t[1]) + for m in eachindex(A.t) + if m != j && m != l + k *= (t - A.t[m]) + end + end + k *= inv(mult) + tmp += k + end end - wi = inv(mult) - tti = t - A.t[idxs[i]] - tmp = wi / (t - A.t[idxs[i]]) - g = tmp * A.u[:, idxs[i]] - @. G += g - @. DG -= g / (t - A.t[idxs[i]]) - F += tmp - DF -= tmp / (t - A.t[idxs[i]]) + @. der += A.u[:, j]*tmp end - @. (DG * F - G * DF) / (F^2) + der end derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number, i) = derivative(A, t), i From 1140d4a37bc49306087b5d834ebfa1e1e42aad25 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Sun, 5 Nov 2023 19:55:31 +0530 Subject: [PATCH 3/3] fixup! test: update derivative tests and make it more robust --- test/derivative_tests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index d1eefe1e..7da2a247 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -31,11 +31,11 @@ function test_derivatives(method, u, t, args...; name::String) adiff = derivative(func, t[1]) if func isa BSplineInterpolation || func isa BSplineApprox # Bug in BSplines - @test isapprox(fdiff, adiff, atol = 1e-8) broken = true + @test_broken isapprox(fdiff, adiff, atol = 1e-8) else @test isapprox(fdiff, adiff, atol = 1e-8) end - + # t = tend fdiff = backward_fdm(5, 1; geom = true)(func, t[end]) adiff = derivative(func, t[end])