From d53e241208c0f9bb3c4404c21f8f33953e55c5f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Mon, 9 Nov 2026 19:15:50 +0100 Subject: [PATCH 1/4] Add rules for logcosh, logabssinh and logabstanh Requires LogExpFunctions 1 (logabstanh was added in 1.0.0). Co-Authored-By: Claude Opus 5.5 --- Project.toml | 2 +- src/rules.jl | 3 +++ test/Project.toml | 2 +- test/runtests.jl | 2 +- 4 files changed, 6 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index e30e835..76d3754 100644 --- a/Project.toml +++ b/Project.toml @@ -11,7 +11,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" [compat] IrrationalConstants = "0.1.1, 0.2" -LogExpFunctions = "0.3.2, 1" +LogExpFunctions = "1" NaNMath = "0.3, 1" Random = "1" SpecialFunctions = "1.1, 2" diff --git a/src/rules.jl b/src/rules.jl index caa3078..e3275fd 100644 --- a/src/rules.jl +++ b/src/rules.jl @@ -300,6 +300,9 @@ _abs_deriv(x) = signbit(x) ? -one(x) : one(x) @define_diffrule LogExpFunctions.logexpm1(x) = :(exp($x - LogExpFunctions.logexpm1($x))) @define_diffrule LogExpFunctions.log1pmx(x) = :(-$x / (1 + $x)) @define_diffrule LogExpFunctions.logmxp1(x) = :((1 - $x) / $x) +@define_diffrule LogExpFunctions.logcosh(x) = :(tanh($x)) +@define_diffrule LogExpFunctions.logabssinh(x) = :(coth($x)) +@define_diffrule LogExpFunctions.logabstanh(x) = :(inv(cosh($x) * sinh($x))) # binary @define_diffrule LogExpFunctions.xlogy(x, y) = diff --git a/test/Project.toml b/test/Project.toml index 9cd1f09..6a0e8fc 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -11,7 +11,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] FiniteDifferences = "0.12" IrrationalConstants = "0.1.1, 0.2" -LogExpFunctions = "0.3.2, 1" +LogExpFunctions = "1" NaNMath = "0.3, 1" Random = "1" SpecialFunctions = "1.1, 2" diff --git a/test/runtests.jl b/test/runtests.jl index f1fa041..079e2e2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -27,7 +27,7 @@ non_diffeable_arg_functions = [(:Base, :rem2pi, 2), (:Base, :ldexp, 2), (:Base, goo = if $(f in (:asec, :acsc, :asecd, :acscd, :acosh, :acoth)) # avoid singularities with finite differencing rand($T) + $T(1.5) - elseif $(f in (:log, :airyaix, :airyaiprimex, :logmxp1)) + elseif $(f in (:log, :airyaix, :airyaiprimex, :logmxp1, :logabssinh, :logabstanh)) # avoid singularities with finite differencing rand($T) + $T(0.5) elseif $(f === :log1mexp) From e753c76519388dfcdd9b7463c9f99733c577a915 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Tue, 10 Nov 2026 09:39:14 +0100 Subject: [PATCH 2/4] Use 2/sinh(2x) for logabstanh and drop obsolete hasmethod guard `2 / sinh(2x)` needs one transcendental call instead of two and is more accurate than `inv(cosh(x) * sinh(x))`. With LogExpFunctions >= 1, all unary rules have Float32 and Float64 methods, so the `hasmethod` guard in the tests is no longer needed. Co-Authored-By: Claude Opus 5.5 --- src/rules.jl | 2 +- test/runtests.jl | 5 +---- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/src/rules.jl b/src/rules.jl index e3275fd..5e0fc17 100644 --- a/src/rules.jl +++ b/src/rules.jl @@ -302,7 +302,7 @@ _abs_deriv(x) = signbit(x) ? -one(x) : one(x) @define_diffrule LogExpFunctions.logmxp1(x) = :((1 - $x) / $x) @define_diffrule LogExpFunctions.logcosh(x) = :(tanh($x)) @define_diffrule LogExpFunctions.logabssinh(x) = :(coth($x)) -@define_diffrule LogExpFunctions.logabstanh(x) = :(inv(cosh($x) * sinh($x))) +@define_diffrule LogExpFunctions.logabstanh(x) = :(2 / sinh(2 * $x)) # binary @define_diffrule LogExpFunctions.xlogy(x, y) = diff --git a/test/runtests.jl b/test/runtests.jl index 079e2e2..6bde69b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,10 +40,7 @@ non_diffeable_arg_functions = [(:Base, :rem2pi, 2), (:Base, :ldexp, 2), (:Base, # We're happy with types with the correct promotion behavior, e.g. # it's fine to return `1` as a derivative despite input being `Float64`. @test promote_type(typeof($deriv), $T) === $T - # In older versions of LogExpFunctions `log1pmx(::Float32)` and `logmxp1(::Float32)` are not defined - if hasmethod($M.$f, Tuple{$T}) - @test $deriv ≈ finitediff($M.$f, goo) rtol=1e-2 atol=1e-3 - end + @test $deriv ≈ finitediff($M.$f, goo) rtol=1e-2 atol=1e-3 # test for 2pi functions if $(f === :mod2pi) goo = 4 * pi From ed7108d836548411b7bf2bdbe3cc208855bcb21b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Tue, 10 Nov 2026 22:11:38 +0100 Subject: [PATCH 3/4] Improve accuracy of LogExpFunctions rules and add missing ones MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Port the rules from JuliaStats/LogExpFunctions.jl#138: more accurate logistic, log1psq, log1mexp, logexpm1, logaddexp and logsubexp rules, and new rules for loglogistic, log1mlogistic, logitexp and logit1mexp. Fixes NaN derivatives of logaddexp(±Inf, ±Inf) (JuliaStats/LogExpFunctions.jl#61). Co-Authored-By: Claude Opus 5.5 --- src/rules.jl | 23 ++++++--- test/runtests.jl | 126 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 142 insertions(+), 7 deletions(-) diff --git a/src/rules.jl b/src/rules.jl index 5e0fc17..bf0312e 100644 --- a/src/rules.jl +++ b/src/rules.jl @@ -291,28 +291,37 @@ _abs_deriv(x) = signbit(x) ? -one(x) : one(x) # unary @define_diffrule LogExpFunctions.xlogx(x) = :(1 + log($x)) -@define_diffrule LogExpFunctions.logistic(x) = :(z = LogExpFunctions.logistic($x); z * (1 - z)) +# logistic(x) * logistic(-x), avoiding cancellation in `1 - logistic(x)` +@define_diffrule LogExpFunctions.logistic(x) = :(let e = exp(-abs($x)); e / (1 + e)^2 end) @define_diffrule LogExpFunctions.logit(x) = :(inv($x * (1 - $x))) -@define_diffrule LogExpFunctions.log1psq(x) = :(2 * $x / (1 + $x^2)) +# avoid overflow of `x^2` +@define_diffrule LogExpFunctions.log1psq(x) = :(abs($x) < 1 ? 2 * $x / (1 + $x^2) : 2 / ($x + inv($x))) @define_diffrule LogExpFunctions.log1pexp(x) = :(LogExpFunctions.logistic($x)) -@define_diffrule LogExpFunctions.log1mexp(x) = :(-exp($x - LogExpFunctions.log1mexp($x))) +@define_diffrule LogExpFunctions.log1mexp(x) = :(-inv(expm1(-$x))) @define_diffrule LogExpFunctions.log2mexp(x) = :(-exp($x - LogExpFunctions.log2mexp($x))) -@define_diffrule LogExpFunctions.logexpm1(x) = :(exp($x - LogExpFunctions.logexpm1($x))) +@define_diffrule LogExpFunctions.logexpm1(x) = :(-inv(expm1(-$x))) @define_diffrule LogExpFunctions.log1pmx(x) = :(-$x / (1 + $x)) @define_diffrule LogExpFunctions.logmxp1(x) = :((1 - $x) / $x) @define_diffrule LogExpFunctions.logcosh(x) = :(tanh($x)) @define_diffrule LogExpFunctions.logabssinh(x) = :(coth($x)) @define_diffrule LogExpFunctions.logabstanh(x) = :(2 / sinh(2 * $x)) +@define_diffrule LogExpFunctions.loglogistic(x) = :(LogExpFunctions.logistic(-$x)) +@define_diffrule LogExpFunctions.log1mlogistic(x) = :(-LogExpFunctions.logistic($x)) +@define_diffrule LogExpFunctions.logitexp(x) = :(-inv(expm1($x))) +@define_diffrule LogExpFunctions.logit1mexp(x) = :(inv(expm1($x))) # binary @define_diffrule LogExpFunctions.xlogy(x, y) = :(log($y)), :(z = $x / $y; iszero($x) && !isnan($y) ? zero(z) : z) +# `x - y` is `NaN` if `x == y == ±Inf` @define_diffrule LogExpFunctions.logaddexp(x, y) = - :(exp($x - LogExpFunctions.logaddexp($x, $y))), :(exp($y - LogExpFunctions.logaddexp($x, $y))) + :(LogExpFunctions.logistic($x == $y && isinf($x) ? zero($x - $y) : $x - $y)), + :(LogExpFunctions.logistic($x == $y && isinf($x) ? zero($x - $y) : $y - $x)) +# `x - y` is `NaN` if `x == y == -Inf` @define_diffrule LogExpFunctions.logsubexp(x, y) = - :(z = LogExpFunctions.logsubexp($x, $y); $x > $y ? exp($x - z) : -exp($x - z)), - :(z = LogExpFunctions.logsubexp($x, $y); $x > $y ? -exp($y - z) : exp($y - z)) + :(-inv(expm1(-($x == $y == -Inf ? zero($x - $y) : $x - $y)))), + :(-inv(expm1($x == $y == -Inf ? zero($x - $y) : $x - $y))) @define_diffrule LogExpFunctions.xlog1py(x, y) = :(log1p($y)), :(z = $x / (1 + $y); iszero($x) && !isnan($y) ? zero(z) : z) diff --git a/test/runtests.jl b/test/runtests.jl index 6bde69b..7e838bf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -32,6 +32,9 @@ non_diffeable_arg_functions = [(:Base, :rem2pi, 2), (:Base, :ldexp, 2), (:Base, rand($T) + $T(0.5) elseif $(f === :log1mexp) rand($T) - one($T) + elseif $(f in (:logitexp, :logit1mexp)) + # avoid singularities with finite differencing + rand($T) - $T(1.5) elseif $(f in (:log2mexp, :erfinv)) rand($T) - $T(0.5) else @@ -183,6 +186,129 @@ for xytype in [:Float32, :Float64, :BigFloat] end end +# Test `Inf` and `NaN` arguments of `logaddexp` (JuliaStats/LogExpFunctions.jl#61) +derivs = DiffRules.diffrule(:LogExpFunctions, :logaddexp, :x, :y) +for xytype in [:Float32, :Float64] + @eval begin + let + x = y = $xytype(-Inf) + dx, dy = $(derivs[1]), $(derivs[2]) + @test dx === dy === $xytype(0.5) + + x = y = $xytype(Inf) + dx, dy = $(derivs[1]), $(derivs[2]) + @test dx === dy === $xytype(0.5) + + x, y = $xytype(Inf), one($xytype) + dx, dy = $(derivs[1]), $(derivs[2]) + @test dx === one($xytype) + @test dy === zero($xytype) + + x, y = $xytype(-Inf), one($xytype) + dx, dy = $(derivs[1]), $(derivs[2]) + @test dx === zero($xytype) + @test dy === one($xytype) + + x, y = $xytype(NaN), one($xytype) + dx, dy = $(derivs[1]), $(derivs[2]) + @test isnan(dx) + @test isnan(dy) + end + end +end + +# Test accuracy of `logaddexp` for arguments of large magnitude +@eval begin + let + x, y = 300.0, 301.0 + dx, dy = $(derivs[1]), $(derivs[2]) + @test dx ≈ LogExpFunctions.logistic(big(x) - big(y)) rtol=4*eps() + @test dy ≈ LogExpFunctions.logistic(big(y) - big(x)) rtol=4*eps() + end +end + +# Test equal and infinite arguments of `logsubexp` +derivs = DiffRules.diffrule(:LogExpFunctions, :logsubexp, :x, :y) +for xytype in [:Float32, :Float64] + @eval begin + let + x = y = one($xytype) + dx, dy = $(derivs[1]), $(derivs[2]) + @test dx === $xytype(Inf) + @test dy === $xytype(-Inf) + + x = y = $xytype(-Inf) + dx, dy = $(derivs[1]), $(derivs[2]) + @test dx === $xytype(Inf) + @test dy === $xytype(-Inf) + + x = y = $xytype(Inf) + dx, dy = $(derivs[1]), $(derivs[2]) + @test isnan(dx) + @test isnan(dy) + + x, y = $xytype(Inf), one($xytype) + dx, dy = $(derivs[1]), $(derivs[2]) + @test dx === one($xytype) + @test dy === -zero($xytype) + + x, y = one($xytype), $xytype(Inf) + dx, dy = $(derivs[1]), $(derivs[2]) + @test dx === -zero($xytype) + @test dy === one($xytype) + end + end +end + +# Test accuracy of `logsubexp` for nearly equal arguments +@eval begin + let + for (x, y) in ((1.0, 1.0 + 1e-6), (1.0 + 1e-6, 1.0)) + dx, dy = $(derivs[1]), $(derivs[2]) + d = big(x) - big(y) + @test dx ≈ -inv(expm1(-d)) rtol=4*eps() + @test dy ≈ -inv(expm1(d)) rtol=4*eps() + end + end +end + +# Test accuracy of `logistic` for arguments of large magnitude +deriv = DiffRules.diffrule(:LogExpFunctions, :logistic, :x) +@eval begin + let + for x in (-40.0, 40.0, -30f0, 30f0) + e = exp(-abs(big(x))) + @test $deriv ≈ e / (1 + e)^2 rtol=4*eps(typeof(x)) + end + end +end + +# Test `log1psq` for arguments of large magnitude +deriv = DiffRules.diffrule(:LogExpFunctions, :log1psq, :x) +@eval begin + let + for x in (-1e200, 1e200) + @test $deriv ≈ 2 / x rtol=4*eps() + end + end +end + +# Test `log1mexp` and `logexpm1` for arguments of small magnitude and `Inf` arguments +for (f, x, xinf, dxinf) in ( + (:log1mexp, -1e-12, -Inf, -0.0), (:logexpm1, 6.039463507942533e-8, Inf, 1.0), +) + deriv = DiffRules.diffrule(:LogExpFunctions, f, :x) + @eval begin + let + x = $x + @test $deriv ≈ -inv(expm1(-big(x))) rtol=4*eps() + + x = $xinf + @test $deriv === $dxinf + end + end +end + end @testset "diffrules" begin From 018adf7de0fb0ef397cd4b229eaa56829aceb04a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Wed, 11 Nov 2026 10:01:13 +0100 Subject: [PATCH 4/4] Update version from 1.17.0 to 1.17.1 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 11ed12d..b1b9422 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "DiffRules" uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" -version = "1.17.0" +version = "1.17.1" [deps] IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"