From 127fc7ca921e593ed9172ff40d6e82e0e54a3117 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Sun, 22 May 2022 15:22:23 +0200 Subject: [PATCH 1/4] Define `SpecialFunctions.gamma_inc` for `ForwardDiff.Dual` --- Project.toml | 4 ++-- src/dual.jl | 14 +++++++++++--- test/DualTest.jl | 16 ++++++++++++++++ 3 files changed, 29 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 2930dc3b..6875c9bb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ForwardDiff" uuid = "f6369f11-7733-5829-9624-2563aa707210" -version = "0.10.30" +version = "0.10.31" [deps] CommonSubexpressions = "bbf7d656-a473-5ed7-a52c-81e309532950" @@ -24,7 +24,7 @@ DiffTests = "0.0.1, 0.1" LogExpFunctions = "0.3" NaNMath = "0.2.2, 0.3, 1" Preferences = "1" -SpecialFunctions = "0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.10, 1.0, 2" +SpecialFunctions = "0.8, 0.9, 0.10, 1.0, 2" StaticArrays = "0.8.3, 0.9, 0.10, 0.11, 0.12, 1.0" julia = "1" diff --git a/src/dual.jl b/src/dual.jl index 8e99c9fe..42682b42 100644 --- a/src/dual.jl +++ b/src/dual.jl @@ -755,9 +755,9 @@ function LinearAlgebra.eigen(A::SymTridiagonal{<:Dual{Tg,T,N}}) where {Tg,T<:Rea Eigen(λ,Dual{Tg}.(Q, tuple.(parts...))) end -# SpecialFunctions.logabsgamma # -# Derivative is not defined in DiffRules # -#----------------------------------------# +# Functions in SpecialFunctions which return tuples # +# Their derivatives are not defined in DiffRules # +#---------------------------------------------------# function SpecialFunctions.logabsgamma(d::Dual{T,<:Real}) where {T} x = value(d) @@ -765,6 +765,14 @@ function SpecialFunctions.logabsgamma(d::Dual{T,<:Real}) where {T} return (Dual{T}(y, SpecialFunctions.digamma(x) * partials(d)), s) end +# Derivatives wrt to first parameter and precision setting are not supported +function SpecialFunctions.gamma_inc(a::Real, d::Dual{T,<:Real}, ind::Integer) where {T} + x = value(d) + p, q = SpecialFunctions.gamma_inc(a, x, ind) + ∂p = exp(-x) * x^(a - 1) / SpecialFunctions.gamma(a) * partials(d) + return (Dual{T}(p, ∂p), Dual{T}(q, -∂p)) +end + ################### # Pretty Printing # ################### diff --git a/test/DualTest.jl b/test/DualTest.jl index 66a1866a..f3f8287c 100644 --- a/test/DualTest.jl +++ b/test/DualTest.jl @@ -540,8 +540,24 @@ for N in (0,3), M in (0,4), V in (Int, Float32) @test dual_isapprox(f(PRIMAL, PRIMAL2, FDNUM3), Dual{TestTag()}(f(PRIMAL, PRIMAL2, PRIMAL3), PARTIALS3)) end + # Functions in Specialfunctions that return tuples and + # therefore are not supported by DiffRules @test dual_isapprox(logabsgamma(FDNUM)[1], loggamma(abs(FDNUM))) @test dual_isapprox(logabsgamma(FDNUM)[2], sign(gamma(FDNUM))) + + a = rand() + fdnum = Dual{TestTag()}(1 + PRIMAL, PARTIALS) # 1 + PRIMAL avoids issues with finite differencing close to 0 + for ind in ((), (0,), (1,), (2,)) + pq = gamma_inc(a, fdnum, ind...) + @test pq isa Tuple{Dual{TestTag()},Dual{TestTag()}} + # We have to adjust tolerances if lower accuracy is requested + # Therefore we don't use `dual_isapprox` + tol = eps(float(V))^(1 / 2^(isempty(ind) ? 1 : 1 + first(ind))) + for i in 1:2 + @test value(pq[i]) ≈ gamma_inc(a, 1 + PRIMAL, ind...)[i] rtol=tol + @test partials(pq[i]) ≈ PARTIALS * Calculus.derivative(x -> gamma_inc(a, x, ind...)[i], 1 + PRIMAL) rtol=tol + end + end end @testset "Exponentiation of zero" begin From 48b46abbdee54093fbd7d5a5d365f3a98ed83bbe Mon Sep 17 00:00:00 2001 From: David Widmann Date: Sun, 22 May 2022 15:39:22 +0200 Subject: [PATCH 2/4] Try to fix tolerances --- test/DualTest.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/DualTest.jl b/test/DualTest.jl index f3f8287c..34794303 100644 --- a/test/DualTest.jl +++ b/test/DualTest.jl @@ -552,7 +552,7 @@ for N in (0,3), M in (0,4), V in (Int, Float32) @test pq isa Tuple{Dual{TestTag()},Dual{TestTag()}} # We have to adjust tolerances if lower accuracy is requested # Therefore we don't use `dual_isapprox` - tol = eps(float(V))^(1 / 2^(isempty(ind) ? 1 : 1 + first(ind))) + tol = (V === Float32 ? 1f-4 : 1e-6)^(1 / 2^(isempty(ind) ? 0 : first(ind))) for i in 1:2 @test value(pq[i]) ≈ gamma_inc(a, 1 + PRIMAL, ind...)[i] rtol=tol @test partials(pq[i]) ≈ PARTIALS * Calculus.derivative(x -> gamma_inc(a, x, ind...)[i], 1 + PRIMAL) rtol=tol From 620ef4631aa2e15876b6b3f4573430aca3c0cb0a Mon Sep 17 00:00:00 2001 From: David Widmann Date: Sun, 22 May 2022 15:56:31 +0200 Subject: [PATCH 3/4] Check if primal method exists --- test/DualTest.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/DualTest.jl b/test/DualTest.jl index 34794303..58a97bf0 100644 --- a/test/DualTest.jl +++ b/test/DualTest.jl @@ -545,9 +545,13 @@ for N in (0,3), M in (0,4), V in (Int, Float32) @test dual_isapprox(logabsgamma(FDNUM)[1], loggamma(abs(FDNUM))) @test dual_isapprox(logabsgamma(FDNUM)[2], sign(gamma(FDNUM))) - a = rand() + a = rand(float(V)) fdnum = Dual{TestTag()}(1 + PRIMAL, PARTIALS) # 1 + PRIMAL avoids issues with finite differencing close to 0 for ind in ((), (0,), (1,), (2,)) + # Only test if primal method exists + # (e.g., older versions of SpecialFunctions don't define `gamma_inc(a, x)` but only `gamma_inc(a, x, ind)` + hasmethod(gamma_inc, typeof((a, 1 + PRIMAL, ind...))) || continue + pq = gamma_inc(a, fdnum, ind...) @test pq isa Tuple{Dual{TestTag()},Dual{TestTag()}} # We have to adjust tolerances if lower accuracy is requested From 4bfeda9349f5671f77995b99f051f040e3baeca6 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Sun, 22 May 2022 16:16:09 +0200 Subject: [PATCH 4/4] Fix tolerances --- test/DualTest.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/DualTest.jl b/test/DualTest.jl index 58a97bf0..eca0da2f 100644 --- a/test/DualTest.jl +++ b/test/DualTest.jl @@ -556,7 +556,8 @@ for N in (0,3), M in (0,4), V in (Int, Float32) @test pq isa Tuple{Dual{TestTag()},Dual{TestTag()}} # We have to adjust tolerances if lower accuracy is requested # Therefore we don't use `dual_isapprox` - tol = (V === Float32 ? 1f-4 : 1e-6)^(1 / 2^(isempty(ind) ? 0 : first(ind))) + tol = V === Float32 ? 5f-4 : 1e-6 + tol = tol^(one(tol) / 2^(isempty(ind) ? 0 : first(ind))) for i in 1:2 @test value(pq[i]) ≈ gamma_inc(a, 1 + PRIMAL, ind...)[i] rtol=tol @test partials(pq[i]) ≈ PARTIALS * Calculus.derivative(x -> gamma_inc(a, x, ind...)[i], 1 + PRIMAL) rtol=tol