diff --git a/Project.toml b/Project.toml index 83f1256..13d160c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "DiffRules" uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" -version = "1.11.1" +version = "1.12.0" [deps] IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6" diff --git a/src/rules.jl b/src/rules.jl index 6e4a1fa..03df753 100644 --- a/src/rules.jl +++ b/src/rules.jl @@ -116,6 +116,9 @@ _abs_deriv(x) = signbit(x) ? -one(x) : one(x) @define_diffrule SpecialFunctions.erfinv(x) = :( ($sqrtπ * exp(SpecialFunctions.erfinv($x)^2)) / 2 ) @define_diffrule SpecialFunctions.erfc(x) = :( -($invsqrtπ * exp(-$x^2) * 2) ) +@define_diffrule SpecialFunctions.logerfc(x) = + :( - 2 * ($invsqrtπ * exp(- $x^2 - SpecialFunctions.logerfc($x))) ) + @define_diffrule SpecialFunctions.erfcinv(x) = :( -($sqrtπ * exp(SpecialFunctions.erfcinv($x)^2)) / 2 ) @define_diffrule SpecialFunctions.erfi(x) = :( $invsqrtπ * exp($x^2) * 2 ) @@ -123,6 +126,7 @@ _abs_deriv(x) = signbit(x) ? -one(x) : one(x) :( 2 * (($x * SpecialFunctions.erfcx($x)) - $invsqrtπ) ) @define_diffrule SpecialFunctions.logerfcx(x) = :( 2 * ($x - inv(SpecialFunctions.erfcx($x) * $sqrtπ)) ) + @define_diffrule SpecialFunctions.dawson(x) = :( 1 - (2 * $x * SpecialFunctions.dawson($x)) ) @define_diffrule SpecialFunctions.digamma(x) = @@ -131,14 +135,35 @@ _abs_deriv(x) = signbit(x) ? -one(x) : one(x) :( inv(SpecialFunctions.trigamma(SpecialFunctions.invdigamma($x))) ) @define_diffrule SpecialFunctions.trigamma(x) = :( SpecialFunctions.polygamma(2, $x) ) + +# derivatives for `airybix` and `airybiprimex` are only correct for real inputs +# `airyaix` and `airyaiprimex` are only defined for positive real inputs +# `airybix` and `airybiprimex` are unscaled for negative real inputs @define_diffrule SpecialFunctions.airyai(x) = :( SpecialFunctions.airyaiprime($x) ) @define_diffrule SpecialFunctions.airyaiprime(x) = :( $x * SpecialFunctions.airyai($x) ) +@define_diffrule SpecialFunctions.airyaix(x) = + :( SpecialFunctions.airyaiprimex($x) + sqrt($x) * SpecialFunctions.airyaix($x) ) +@define_diffrule SpecialFunctions.airyaiprimex(x) = + :( $x * SpecialFunctions.airyaix($x) + sqrt($x) * SpecialFunctions.airyaiprimex($x) ) @define_diffrule SpecialFunctions.airybi(x) = :( SpecialFunctions.airybiprime($x) ) @define_diffrule SpecialFunctions.airybiprime(x) = :( $x * SpecialFunctions.airybi($x) ) +@define_diffrule SpecialFunctions.airybix(x) = + :( if $x > zero($x) + SpecialFunctions.airybiprimex($x) - sqrt($x) * SpecialFunctions.airybix($x) + else + SpecialFunctions.airybiprimex($x) + end ) +@define_diffrule SpecialFunctions.airybiprimex(x) = + :( if $x > zero($x) + $x * SpecialFunctions.airybix($x) - sqrt($x) * SpecialFunctions.airybiprimex($x) + else + $x * SpecialFunctions.airybix($x) + end ) + @define_diffrule SpecialFunctions.besselj0(x) = :( -SpecialFunctions.besselj1($x) ) @define_diffrule SpecialFunctions.besselj1(x) = @@ -148,49 +173,72 @@ _abs_deriv(x) = signbit(x) ? -one(x) : one(x) @define_diffrule SpecialFunctions.bessely1(x) = :( (SpecialFunctions.bessely0($x) - SpecialFunctions.bessely(2, $x)) / 2 ) +@define_diffrule SpecialFunctions.sinint(x) = :( sinc($x / π) ) +@define_diffrule SpecialFunctions.cosint(x) = :( cos($x) / $x ) + +@define_diffrule SpecialFunctions.ellipk(m) = + :( (SpecialFunctions.ellipe($m) / (1 - $m) - SpecialFunctions.ellipk($m)) / (2 * $m) ) +@define_diffrule SpecialFunctions.ellipe(m) = + :( (SpecialFunctions.ellipe($m) - SpecialFunctions.ellipk($m)) / (2 * $m) ) + # TODO: # # eta # zeta -# airyaix -# airyaiprimex -# airybix -# airybiprimex # binary # #--------# +@define_diffrule SpecialFunctions.erf(x, y) = + :( -2 * ($invsqrtπ * exp(-$x^2)) ), :( 2 * ($invsqrtπ * exp(-$y^2)) ) + +# derivatives with respect to the order `ν` exist but are not implemented +# (analogously to the ChainRules definitions in SpecialFunctions) + +# derivatives for `besselix`, `besseljx` and `besselyx` are only correct for real inputs +# see https://github.com/JuliaMath/SpecialFunctions.jl/blob/master/src/chainrules.jl +# for forward-mode and reverse-mode derivatives for complex inputs + @define_diffrule SpecialFunctions.besselj(ν, x) = :NaN, :( (SpecialFunctions.besselj($ν - 1, $x) - SpecialFunctions.besselj($ν + 1, $x)) / 2 ) +@define_diffrule SpecialFunctions.besseljx(ν, x) = + :NaN, :( (SpecialFunctions.besseljx($ν - 1, $x) - SpecialFunctions.besseljx($ν + 1, $x)) / 2 ) @define_diffrule SpecialFunctions.besseli(ν, x) = :NaN, :( (SpecialFunctions.besseli($ν - 1, $x) + SpecialFunctions.besseli($ν + 1, $x)) / 2 ) +@define_diffrule SpecialFunctions.besselix(ν, x) = + :NaN, :( (SpecialFunctions.besselix($ν - 1, $x) + SpecialFunctions.besselix($ν + 1, $x)) / 2 - sign($x) * SpecialFunctions.besselix($ν, $x) ) @define_diffrule SpecialFunctions.bessely(ν, x) = :NaN, :( (SpecialFunctions.bessely($ν - 1, $x) - SpecialFunctions.bessely($ν + 1, $x)) / 2 ) +@define_diffrule SpecialFunctions.besselyx(ν, x) = + :NaN, :( (SpecialFunctions.besselyx($ν - 1, $x) - SpecialFunctions.besselyx($ν + 1, $x)) / 2 ) @define_diffrule SpecialFunctions.besselk(ν, x) = :NaN, :( -(SpecialFunctions.besselk($ν - 1, $x) + SpecialFunctions.besselk($ν + 1, $x)) / 2 ) +@define_diffrule SpecialFunctions.besselkx(ν, x) = + :NaN, :( -(SpecialFunctions.besselkx($ν - 1, $x) + SpecialFunctions.besselkx($ν + 1, $x)) / 2 + SpecialFunctions.besselkx($ν, $x) ) +@define_diffrule SpecialFunctions.besselh(ν, x) = + :NaN, :( (SpecialFunctions.besselh($ν - 1, $x) - SpecialFunctions.besselh($ν + 1, $x)) / 2 ) +@define_diffrule SpecialFunctions.besselhx(ν, x) = + :NaN, :( (SpecialFunctions.besselhx($ν - 1, $x) - SpecialFunctions.besselhx($ν + 1, $x)) / 2 - im * SpecialFunctions.besselhx($ν, $x) ) @define_diffrule SpecialFunctions.hankelh1(ν, x) = :NaN, :( (SpecialFunctions.hankelh1($ν - 1, $x) - SpecialFunctions.hankelh1($ν + 1, $x)) / 2 ) +@define_diffrule SpecialFunctions.hankelh1x(ν, x) = + :NaN, :( (SpecialFunctions.hankelh1x($ν - 1, $x) - SpecialFunctions.hankelh1x($ν + 1, $x)) / 2 - im * SpecialFunctions.hankelh1x($ν, $x) ) @define_diffrule SpecialFunctions.hankelh2(ν, x) = :NaN, :( (SpecialFunctions.hankelh2($ν - 1, $x) - SpecialFunctions.hankelh2($ν + 1, $x)) / 2 ) +@define_diffrule SpecialFunctions.hankelh2x(ν, x) = + :NaN, :( (SpecialFunctions.hankelh2x($ν - 1, $x) - SpecialFunctions.hankelh2x($ν + 1, $x)) / 2 + im * SpecialFunctions.hankelh2x($ν, $x) ) + @define_diffrule SpecialFunctions.polygamma(m, x) = :NaN, :( SpecialFunctions.polygamma($m + 1, $x) ) + @define_diffrule SpecialFunctions.beta(a, b) = :( SpecialFunctions.beta($a, $b)*(SpecialFunctions.digamma($a) - SpecialFunctions.digamma($a + $b)) ), :( SpecialFunctions.beta($a, $b)*(SpecialFunctions.digamma($b) - SpecialFunctions.digamma($a + $b)) ) -@define_diffrule SpecialFunctions.logbeta(a, b) = +@define_diffrule SpecialFunctions.logbeta(a, b) = :( SpecialFunctions.digamma($a) - SpecialFunctions.digamma($a + $b) ), :( SpecialFunctions.digamma($b) - SpecialFunctions.digamma($a + $b) ) -# TODO: -# -# zeta -# besseljx -# besselyx -# besselix -# besselkx -# besselh -# besselhx -# hankelh1x -# hankelh2 -# hankelh2x +# derivative wrt to `s` is not implemented +@define_diffrule SpecialFunctions.zeta(s, z) = + :NaN, :( - $s * SpecialFunctions.zeta($s + 1, $z) ) # ternary # #---------# diff --git a/test/runtests.jl b/test/runtests.jl index 3312f7f..f1fa041 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -118,7 +118,7 @@ for xtype in [:Float64, :BigFloat, :Int64] x = $xtype(rand(1 : 10)) y = $mode dx, dy = $(derivs[1]), $(derivs[2]) - @test isapprox(dx, finitediff(z -> rem2pi(z, y), float(x)), rtol=0.05) + @test dx ≈ finitediff(z -> rem2pi(z, y), float(x)) rtol=1e-9 atol=1e-9 @test isnan(dy) end end @@ -134,13 +134,24 @@ for xtype in [:Float64, :BigFloat] x = rand($xtype) y = $ytype(rand(1 : 10)) dx, dy = $(derivs[1]), $(derivs[2]) - @test isapprox(dx, finitediff(z -> ldexp(z, y), x), rtol=0.05) + @test dx ≈ finitediff(z -> ldexp(z, y), x) rtol=1e-9 atol=1e-9 @test isnan(dy) end end end end +# Check negative branch for `airybix` and `airybiprimex` +for f in (:airybix, :airybiprimex) + deriv = DiffRules.diffrule(:SpecialFunctions, f, :goo) + @eval begin + let + goo = -rand() + @test $deriv ≈ finitediff(SpecialFunctions.$f, goo) rtol=1e-9 atol=1e-9 + end + end +end + # Test `iszero(x)` branch of `xlogy` derivs = DiffRules.diffrule(:LogExpFunctions, :xlogy, :x, :y) for xytype in [:Float32, :Float64, :BigFloat]