Update rules for SpecialFunctions

Project.toml

@@ -1,6 +1,6 @@
 name = "ChainRules"
 uuid = "082447d4-558c-5d27-93f4-14fc19e9eca2"
-version = "0.7.61"
+version = "0.7.62"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -12,8 +12,8 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-ChainRulesCore = "0.9.29"
-ChainRulesTestUtils = "0.6.6"
+ChainRulesCore = "0.9.40"
+ChainRulesTestUtils = "0.6.8"
 Compat = "3"
 FiniteDifferences = "0.11, 0.12"
 Reexport = "0.2, 1"

src/rulesets/packages/SpecialFunctions.jl

@@ -1,3 +1,8 @@
+const BESSEL_ORDER_INFO = """
+derivatives of Bessel functions with respect to the order are not implemented currently:
+https://github.com/JuliaMath/SpecialFunctions.jl/issues/160
+"""
+
 @scalar_rule(SpecialFunctions.airyai(x), SpecialFunctions.airyaiprime(x))
 @scalar_rule(SpecialFunctions.airyaiprime(x), x * SpecialFunctions.airyai(x))
 @scalar_rule(SpecialFunctions.airybi(x), SpecialFunctions.airybiprime(x))
@@ -30,43 +35,60 @@
 # binary
 @scalar_rule(
     SpecialFunctions.besselj(ν, x),
-    (NaN, (SpecialFunctions.besselj(ν - 1, x) - SpecialFunctions.besselj(ν + 1, x)) / 2),
+    (
+        @not_implemented(BESSEL_ORDER_INFO),
+        (SpecialFunctions.besselj(ν - 1, x) - SpecialFunctions.besselj(ν + 1, x)) / 2
+    ),
 )
 @scalar_rule(
     SpecialFunctions.besseli(ν, x),
-    (NaN, (SpecialFunctions.besseli(ν - 1, x) + SpecialFunctions.besseli(ν + 1, x)) / 2),
+    (
+        @not_implemented(BESSEL_ORDER_INFO),
+        (SpecialFunctions.besseli(ν - 1, x) + SpecialFunctions.besseli(ν + 1, x)) / 2,
+    ),
 )
 @scalar_rule(
     SpecialFunctions.bessely(ν, x),
-    (NaN, (SpecialFunctions.bessely(ν - 1, x) - SpecialFunctions.bessely(ν + 1, x)) / 2),
+    (
+        @not_implemented(BESSEL_ORDER_INFO),
+        (SpecialFunctions.bessely(ν - 1, x) - SpecialFunctions.bessely(ν + 1, x)) / 2,
+    ),
 )
 @scalar_rule(
     SpecialFunctions.besselk(ν, x),
-    (NaN, -(SpecialFunctions.besselk(ν - 1, x) + SpecialFunctions.besselk(ν + 1, x)) / 2),
+    (
+        @not_implemented(BESSEL_ORDER_INFO),
+        -(SpecialFunctions.besselk(ν - 1, x) + SpecialFunctions.besselk(ν + 1, x)) / 2,
+    ),
 )
 @scalar_rule(
     SpecialFunctions.hankelh1(ν, x),
-    (NaN, (SpecialFunctions.hankelh1(ν - 1, x) - SpecialFunctions.hankelh1(ν + 1, x)) / 2),
+    (
+        @not_implemented(BESSEL_ORDER_INFO),
+        (SpecialFunctions.hankelh1(ν - 1, x) - SpecialFunctions.hankelh1(ν + 1, x)) / 2,
+    ),
 )
 @scalar_rule(
     SpecialFunctions.hankelh2(ν, x),
-    (NaN, (SpecialFunctions.hankelh2(ν - 1, x) - SpecialFunctions.hankelh2(ν + 1, x)) / 2),
+    (
+        @not_implemented(BESSEL_ORDER_INFO),
+        (SpecialFunctions.hankelh2(ν - 1, x) - SpecialFunctions.hankelh2(ν + 1, x)) / 2,
+    ),
 )
 @scalar_rule(
     SpecialFunctions.polygamma(m, x),
-    (NaN, SpecialFunctions.polygamma(m + 1, x))
+    (
+        DoesNotExist(),
+        SpecialFunctions.polygamma(m + 1, x),
+    ),
 )
 # todo: setup for common expr
 @scalar_rule(
     SpecialFunctions.beta(a, b),
     (Ω*(SpecialFunctions.digamma(a) - SpecialFunctions.digamma(a + b)),
      Ω*(SpecialFunctions.digamma(b) - SpecialFunctions.digamma(a + b)),)
 )
-@scalar_rule(
-    SpecialFunctions.lbeta(a, b),
-    (SpecialFunctions.digamma(a) - SpecialFunctions.digamma(a + b),
-     SpecialFunctions.digamma(b) - SpecialFunctions.digamma(a + b),)
-)
+
 # Changes between SpecialFunctions 0.7 and 0.8
 if isdefined(SpecialFunctions, :lgamma)
     # actually is the absolute value of the logorithm of gamma
@@ -81,3 +103,21 @@ end
 if isdefined(SpecialFunctions, :loggamma)
     @scalar_rule(SpecialFunctions.loggamma(x), SpecialFunctions.digamma(x))
 end
+
+if isdefined(SpecialFunctions, :lbeta)
+    # todo: setup for common expr
+    @scalar_rule(
+        SpecialFunctions.lbeta(a, b),
+        (SpecialFunctions.digamma(a) - SpecialFunctions.digamma(a + b),
+         SpecialFunctions.digamma(b) - SpecialFunctions.digamma(a + b),)
+    )
+end
+
+if isdefined(SpecialFunctions, :logbeta)
+    # todo: setup for common expr
+    @scalar_rule(
+        SpecialFunctions.logbeta(a, b),
+        (SpecialFunctions.digamma(a) - SpecialFunctions.digamma(a + b),
+         SpecialFunctions.digamma(b) - SpecialFunctions.digamma(a + b),)
+    )
+end

test/rulesets/packages/SpecialFunctions.jl

@@ -1,53 +1,109 @@
-@testset "SpecialFunctions" for x in (1.0, -1.0, 0.0, 0.5, 10.0, -17.1, 1.5 + 0.7im)
-    test_scalar(SpecialFunctions.erf, x)
-    test_scalar(SpecialFunctions.erfc, x)
-    test_scalar(SpecialFunctions.erfi, x)
+@testset "general: single input" begin
+    for x in (1.0, -1.0, 0.0, 0.5, 10.0, -17.1, 1.5 + 0.7im)
+        test_scalar(SpecialFunctions.erf, x)
+        test_scalar(SpecialFunctions.erfc, x)
+        test_scalar(SpecialFunctions.erfi, x)
 
-    test_scalar(SpecialFunctions.airyai, x)
-    test_scalar(SpecialFunctions.airyaiprime, x)
-    test_scalar(SpecialFunctions.airybi, x)
-    test_scalar(SpecialFunctions.airybiprime, x)
+        test_scalar(SpecialFunctions.airyai, x)
+        test_scalar(SpecialFunctions.airyaiprime, x)
+        test_scalar(SpecialFunctions.airybi, x)
+        test_scalar(SpecialFunctions.airybiprime, x)
 
-    test_scalar(SpecialFunctions.besselj0, x)
-    test_scalar(SpecialFunctions.besselj1, x)
+        test_scalar(SpecialFunctions.erfcx, x)
+        test_scalar(SpecialFunctions.dawson, x)
 
-    test_scalar(SpecialFunctions.erfcx, x)
-    test_scalar(SpecialFunctions.dawson, x)
+        if x isa Real
+            test_scalar(SpecialFunctions.invdigamma, x)
+        end
 
-    if x isa Real
-        test_scalar(SpecialFunctions.invdigamma, x)
-    end
+        if x isa Real && 0 < x < 1
+            test_scalar(SpecialFunctions.erfinv, x)
+            test_scalar(SpecialFunctions.erfcinv, x)
+        end
 
-    if x isa Real && 0 < x < 1
-        test_scalar(SpecialFunctions.erfinv, x)
-        test_scalar(SpecialFunctions.erfcinv, x)
+        if x isa Real && x > 0 || x isa Complex
+            test_scalar(SpecialFunctions.gamma, x)
+            test_scalar(SpecialFunctions.digamma, x)
+            test_scalar(SpecialFunctions.trigamma, x)
+        end
     end
+end
+
+@testset "Bessel functions" begin
+    for x in (1.5, 2.5, 10.5, -0.6, -2.6, -3.3, 1.6 + 1.6im, 1.6 - 1.6im, -4.6 + 1.6im)
+        test_scalar(SpecialFunctions.besselj0, x)
+        test_scalar(SpecialFunctions.besselj1, x)
+
+        isreal(x) && x < 0 && continue
 
-    if x isa Real && x > 0 || x isa Complex
         test_scalar(SpecialFunctions.bessely0, x)
         test_scalar(SpecialFunctions.bessely1, x)
-        test_scalar(SpecialFunctions.gamma, x)
-        test_scalar(SpecialFunctions.digamma, x)
-        test_scalar(SpecialFunctions.trigamma, x)
+
+        for nu in (-1.5, 2.2, 4.0)
+            test_frule(SpecialFunctions.besseli, nu, x)
+            test_rrule(SpecialFunctions.besseli, nu, x)
+
+            test_frule(SpecialFunctions.besselj, nu, x)
+            test_rrule(SpecialFunctions.besselj, nu, x)
+
+            test_frule(SpecialFunctions.besselk, nu, x)
+            test_rrule(SpecialFunctions.besselk, nu, x)
+
+            test_frule(SpecialFunctions.bessely, nu, x)
+            test_rrule(SpecialFunctions.bessely, nu, x)
+
+            # use complex numbers in `rrule` for FiniteDifferences
+            test_frule(SpecialFunctions.hankelh1, nu, x)
+            test_rrule(SpecialFunctions.hankelh1, nu, complex(x))
+
+            # use complex numbers in `rrule` for FiniteDifferences
+            test_frule(SpecialFunctions.hankelh2, nu, x)
+            test_rrule(SpecialFunctions.hankelh2, nu, complex(x))
+        end
+    end
+end
+
+@testset "beta and logbeta" begin
+    test_points = (1.5, 2.5, 10.5, 1.6 + 1.6im, 1.6 - 1.6im, 4.6 + 1.6im)
+    for _x in test_points, _y in test_points
+        # ensure all complex if any complex for FiniteDifferences
+        x, y = promote(_x, _y)
+        test_frule(SpecialFunctions.beta, x, y)
+        test_rrule(SpecialFunctions.beta, x, y)
+
+        if isdefined(SpecialFunctions, :lbeta)
+            test_frule(SpecialFunctions.lbeta, x, y)
+            test_rrule(SpecialFunctions.lbeta, x, y)
+        end
+
+        if isdefined(SpecialFunctions, :logbeta)
+            test_frule(SpecialFunctions.logbeta, x, y)
+            test_rrule(SpecialFunctions.logbeta, x, y)
+        end
     end
 end
 
-# SpecialFunctions 0.7->0.8 changes:
 @testset "log gamma and co" begin
-    #It is important that we have negative numbers with both odd and even integer parts
-    for x in (1.5, 2.5, 10.5, -0.6, -2.6, -3.3, 1.6+1.6im, 1.6-1.6im, -4.6+1.6im)
+    # It is important that we have negative numbers with both odd and even integer parts
+    for x in (1.5, 2.5, 10.5, -0.6, -2.6, -3.3, 1.6 + 1.6im, 1.6 - 1.6im, -4.6 + 1.6im)
+        for m in (0, 1, 2, 3)
+            test_frule(SpecialFunctions.polygamma, m, x)
+            test_rrule(SpecialFunctions.polygamma, m, x)
+        end
+
         if isdefined(SpecialFunctions, :lgamma)
             test_scalar(SpecialFunctions.lgamma, x)
         end
+
         if isdefined(SpecialFunctions, :loggamma)
             isreal(x) && x < 0 && continue
             test_scalar(SpecialFunctions.loggamma, x)
         end
 
         if isdefined(SpecialFunctions, :logabsgamma)
             isreal(x) || continue
-            test_frule(logabsgamma, x)
-            test_rrule(logabsgamma, x; output_tangent=(randn(), randn()))
+            test_frule(SpecialFunctions.logabsgamma, x)
+            test_rrule(SpecialFunctions.logabsgamma, x; output_tangent=(randn(), randn()))
         end
     end
 end