Add xlog1py (#27)

* Add `xlog1py`

* Bump version

Author: devmotion

Project.toml

@@ -1,7 +1,7 @@
 name = "LogExpFunctions"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 authors = ["StatsFun.jl contributors, Tamas K. Papp <[email protected]>"]
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

docs/src/index.md

@@ -7,6 +7,7 @@ The original authors of these functions are the StatsFuns.jl contributors.
 ```@docs
 xlogx
 xlogy
+xlog1py
 logistic
 logit
 log1psq

src/LogExpFunctions.jl

@@ -7,7 +7,7 @@ import ChainRulesCore
 import IrrationalConstants
 import LinearAlgebra
 
-export xlogx, xlogy, logistic, logit, log1psq, log1pexp, log1mexp, log2mexp, logexpm1,
+export xlogx, xlogy, xlog1py, logistic, logit, log1psq, log1pexp, log1mexp, log2mexp, logexpm1,
     softplus, invsoftplus, log1pmx, logmxp1, logaddexp, logsubexp, logsumexp, softmax,
     softmax!
 

src/basicfuns.jl

@@ -29,6 +29,21 @@ function xlogy(x::Number, y::Number)
     return iszero(x) && !isnan(y) ? zero(result) : result
 end
 
+"""
+$(SIGNATURES)
+
+Return `x * log(1 + y)` for `y ≥ -1` with correct limit at ``x = 0``.
+
+```jldoctest
+julia> xlog1py(0, -1)
+0.0
+```
+"""
+function xlog1py(x::Number, y::Number)
+    result = x * log1p(y)
+    return iszero(x) && !isnan(y) ? zero(result) : result
+end
+
 # The following bounds are precomputed versions of the following abstract
 # function, but the implicit interface for AbstractFloat doesn't uniformly
 # enforce that all floating point types implement nextfloat and prevfloat.

src/chainrules.jl

@@ -1,5 +1,6 @@
 ChainRulesCore.@scalar_rule(xlogx(x::Real), (1 + log(x),))
 ChainRulesCore.@scalar_rule(xlogy(x::Real, y::Real), (log(y), x / y,))
+ChainRulesCore.@scalar_rule(xlog1py(x::Real, y::Real), (log1p(y), x / (1 + y),))
 
 ChainRulesCore.@scalar_rule(logistic(x::Real), (Ω * (1 - Ω),))
 ChainRulesCore.@scalar_rule(logit(x::Real), (inv(x * (1 - x)),))

test/basicfuns.jl

@@ -1,4 +1,4 @@
-@testset "xlogx & xlogy" begin
+@testset "xlogx, xlogy, and xlog1py" begin
     @test iszero(xlogx(0))
     @test xlogx(2) ≈ 2.0 * log(2.0)
     @test_throws DomainError xlogx(-1)
@@ -11,6 +11,13 @@
     @test isnan(xlogy(2, NaN))
     @test isnan(xlogy(0, NaN))
 
+    @test iszero(xlog1py(0, 0))
+    @test xlog1py(2, 3) ≈ 2.0 * log1p(3.0)
+    @test_throws DomainError xlog1py(1, -2)
+    @test isnan(xlog1py(NaN, 2))
+    @test isnan(xlog1py(2, NaN))
+    @test isnan(xlog1py(0, NaN))
+
     # Since we allow complex/negative values, test for them. See comments in:
     # https://github.com/JuliaStats/StatsFuns.jl/pull/95
 
@@ -26,6 +33,15 @@
     @test isnan(xlogy(Inf + im * NaN, 1))
     @test isnan(xlogy(0 + im * 0, NaN))
     @test iszero(xlogy(0 + im * 0, 0 + im * Inf))
+
+    @test xlog1py(-2, 3) == -xlog1py(2, 3)
+    @test xlog1py(1 + im, 3) == (1 + im) * log1p(3)
+    @test xlog1py(1 + im, 2 + im) == (1 + im) * log1p(2 + im)
+    @test isnan(xlog1py(1 + NaN * im, -1 + im))
+    @test isnan(xlog1py(0, -1 + NaN * im))
+    @test isnan(xlog1py(Inf + im * NaN, 1))
+    @test isnan(xlog1py(0 + im * 0, NaN))
+    @test iszero(xlog1py(0 + im * 0, -1 + im * Inf))
 end
 
 @testset "logistic & logit" begin

test/chainrules.jl

@@ -6,6 +6,11 @@
         y = rand()
         test_frule(xlogy, x, y)
         test_rrule(xlogy, x, y)
+
+        for z in (-y, y)
+            test_frule(xlog1py, x, z)
+            test_rrule(xlog1py, x, z)
+        end
     end
 
     test_frule(logit, x)