Merge remote-tracking branch 'upstream/master' into add_smoothed_linear_interpolation

Author: SouthEndMusic

Project.toml

@@ -1,6 +1,6 @@
 name = "DataInterpolations"
 uuid = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
-version = "8.1.0"
+version = "8.3.1"
 
 [deps]
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
@@ -13,25 +13,29 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [weakdeps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 RegularizationTools = "29dad682-9a27-4bc3-9c72-016788665182"
+SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [extensions]
 DataInterpolationsChainRulesCoreExt = "ChainRulesCore"
 DataInterpolationsOptimExt = "Optim"
 DataInterpolationsRegularizationToolsExt = "RegularizationTools"
+DataInterpolationsSparseConnectivityTracerExt = ["SparseConnectivityTracer", "FillArrays"]
 DataInterpolationsSymbolicsExt = "Symbolics"
 
 [compat]
 Aqua = "0.8"
 BenchmarkTools = "1"
 ChainRulesCore = "1.24"
 EnumX = "1.0.4"
+FillArrays = "1.13.0"
 FindFirstFunctions = "1.3"
 FiniteDifferences = "0.12.31"
-ForwardDiff = "0.10.36"
+ForwardDiff = "0.10.36, 1"
 LinearAlgebra = "1.10"
 Optim = "1.6"
 PrettyTables = "2"
@@ -40,6 +44,7 @@ RecipesBase = "1.3"
 Reexport = "1"
 RegularizationTools = "0.6"
 SafeTestsets = "0.1"
+SparseConnectivityTracer = "1"
 StableRNGs = "1"
 Symbolics = "5.29, 6"
 Test = "1.10"
@@ -57,11 +62,12 @@ Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 RegularizationTools = "29dad682-9a27-4bc3-9c72-016788665182"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "BenchmarkTools", "SafeTestsets", "ChainRulesCore", "Optim", "RegularizationTools", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "Symbolics", "Unitful", "Zygote"]
+test = ["Aqua", "BenchmarkTools", "SafeTestsets", "ChainRulesCore", "Optim", "RegularizationTools", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "Symbolics", "Unitful", "Zygote", "SparseConnectivityTracer"]

ext/DataInterpolationsSparseConnectivityTracerExt.jl

@@ -0,0 +1,120 @@
+module DataInterpolationsSparseConnectivityTracerExt
+
+using SparseConnectivityTracer: AbstractTracer, Dual, primal, tracer
+using SparseConnectivityTracer: GradientTracer, gradient_tracer_1_to_1
+using SparseConnectivityTracer: HessianTracer, hessian_tracer_1_to_1
+using FillArrays: Fill # from FillArrays.jl
+using DataInterpolations:
+    AbstractInterpolation,
+    LinearInterpolation,
+    QuadraticInterpolation,
+    LagrangeInterpolation,
+    AkimaInterpolation,
+    ConstantInterpolation,
+    QuadraticSpline,
+    CubicSpline,
+    BSplineInterpolation,
+    BSplineApprox,
+    CubicHermiteSpline,
+    # PCHIPInterpolation,
+    QuinticHermiteSpline,
+    output_size
+
+#===========#
+# Utilities #
+#===========#
+
+# Limit support to `u` begin an AbstractVector{<:Number} or AbstractMatrix{<:Number},
+# to avoid any cases where the output size is dependent on the input value.
+# https://github.com/adrhill/SparseConnectivityTracer.jl/pull/234#discussion_r2031038566
+
+function _sct_interpolate(
+        ::AbstractInterpolation,
+        uType::Type{<:AbstractVector{<:Number}},
+        t::GradientTracer,
+        is_der_1_zero,
+        is_der_2_zero,
+    )
+    return gradient_tracer_1_to_1(t, is_der_1_zero)
+end
+function _sct_interpolate(
+        ::AbstractInterpolation,
+        uType::Type{<:AbstractVector{<:Number}},
+        t::HessianTracer,
+        is_der_1_zero,
+        is_der_2_zero,
+    )
+    return hessian_tracer_1_to_1(t, is_der_1_zero, is_der_2_zero)
+end
+function _sct_interpolate(
+        interp::AbstractInterpolation,
+        uType::Type{<:AbstractMatrix{<:Number}},
+        t::GradientTracer,
+        is_der_1_zero,
+        is_der_2_zero,
+    )
+    t = gradient_tracer_1_to_1(t, is_der_1_zero)
+    N = only(output_size(interp))
+    return Fill(t, N)
+end
+function _sct_interpolate(
+        interp::AbstractInterpolation,
+        uType::Type{<:AbstractMatrix{<:Number}},
+        t::HessianTracer,
+        is_der_1_zero,
+        is_der_2_zero,
+    )
+    t = hessian_tracer_1_to_1(t, is_der_1_zero, is_der_2_zero)
+    N = only(output_size(interp))
+    return Fill(t, N)
+end
+
+#===========#
+# Overloads #
+#===========#
+
+# We assume that with the exception of ConstantInterpolation and LinearInterpolation,
+# all interpolations have a non-zero second derivative at some point in the input domain.
+
+for (I, is_der1_zero, is_der2_zero) in (
+        (:ConstantInterpolation, true, true),
+        (:LinearInterpolation, false, true),
+        (:QuadraticInterpolation, false, false),
+        (:LagrangeInterpolation, false, false),
+        (:AkimaInterpolation, false, false),
+        (:QuadraticSpline, false, false),
+        (:CubicSpline, false, false),
+        (:BSplineInterpolation, false, false),
+        (:BSplineApprox, false, false),
+        (:CubicHermiteSpline, false, false),
+        (:QuinticHermiteSpline, false, false),
+    )
+    @eval function (interp::$(I){uType})(
+            t::AbstractTracer
+        ) where {uType <: AbstractArray{<:Number}}
+        return _sct_interpolate(interp, uType, t, $is_der1_zero, $is_der2_zero)
+    end
+end
+
+# Some Interpolations require custom overloads on `Dual` due to mutation of caches.
+for I in (
+        :LagrangeInterpolation,
+        :BSplineInterpolation,
+        :BSplineApprox,
+        :CubicHermiteSpline,
+        :QuinticHermiteSpline,
+    )
+    @eval function (interp::$(I){uType})(d::Dual) where {uType <: AbstractVector}
+        p = interp(primal(d))
+        t = interp(tracer(d))
+        return Dual(p, t)
+    end
+
+    @eval function (interp::$(I){uType})(d::Dual) where {uType <: AbstractMatrix}
+        p = interp(primal(d))
+        t = interp(tracer(d))
+        return Dual.(p, t)
+    end
+end
+
+end

src/derivatives.jl

@@ -6,10 +6,12 @@ function derivative(A, t, order = 1)
         _extrapolate_derivative_right(A, t, order)
     else
         iguess = A.iguesser
-        (order == 1) ? _derivative(A, t, iguess) :
-        ForwardDiff.derivative(t -> begin
-                _derivative(A, t, iguess)
-            end, t)
+        if order == 1
+            return _derivative(A, t, iguess)
+        end
+        return ForwardDiff.derivative(t -> begin
+                -_derivative(A, -t, iguess)
+            end, -t) # take derivative backwards in t to make it a left rather than right derivative
     end
 end
 
@@ -333,9 +335,8 @@ function _derivative(
         ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / (A.k[A.d + 2])
     else
         for i in 1:(A.h - 1)
-            ducum = ducum +
-                    sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) /
-                    (A.k[i + A.d + 1] - A.k[i + 1])
+            ducum += sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) /
+                     (A.k[i + A.d + 1] - A.k[i + 1])
         end
     end
     ducum * A.d * scale

test/derivative_tests.jl

@@ -6,7 +6,7 @@ using Symbolics
 using StableRNGs
 using RegularizationTools
 using Optim
-using ForwardDiff
+import ForwardDiff
 using LinearAlgebra
 
 function test_derivatives(method; args = [], kwargs = [], name::String)
@@ -35,11 +35,8 @@ function test_derivatives(method; args = [], kwargs = [], name::String)
 
         # Interpolation transition points
         for _t in t[2:(end - 1)]
-            if func isa Union{BSplineInterpolation, BSplineApprox,
-                CubicHermiteSpline}
-                fdiff = forward_fdm(5, 1; geom = true)(func, _t)
-                fdiff2 = forward_fdm(5, 1; geom = true)(t -> derivative(func, t), _t)
-            elseif func isa SmoothedConstantInterpolation
+            if func isa Union{SmoothedConstantInterpolation, BSplineInterpolation, BSplineApprox}
+                # TODO fix interpolations
                 continue
             else
                 fdiff = backward_fdm(5, 1; geom = true)(func, _t)

test/runtests.jl

@@ -8,6 +8,7 @@ using SafeTestsets
 @safetestset "Integral Tests" include("integral_tests.jl")
 @safetestset "Integral Inverse Tests" include("integral_inverse_tests.jl")
 @safetestset "Extrapolation Tests" include("extrapolation_tests.jl")
+@safetestset "SparseConnectivityTracer Tests" include("sparseconnectivitytracer_tests.jl")
 @safetestset "Online Tests" include("online_tests.jl")
 @safetestset "Regularization Smoothing Tests" include("regularization.jl")
 @safetestset "Show methods Tests" include("show.jl")