Fix logsubexp(-Inf, -Inf) (#21)

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.2.4"
+version = "0.2.5"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

src/basicfuns.jl

@@ -221,7 +221,11 @@ $(SIGNATURES)
 
 Return `log(abs(exp(x) - exp(y)))`, preserving numerical accuracy.
 """
-logsubexp(x::Real, y::Real) = max(x, y) + log1mexp(-abs(x - y))
+function logsubexp(x::Real, y::Real)
+    # ensure that `Δ = 0` if `x = y = - Inf` (but not for `x = y = +Inf`!)
+    Δ = x == y && (isfinite(x) || x < 0) ? zero(x - y) : abs(x - y)
+    return max(x, y) + log1mexp(-Δ)
+end
 
 """
 $(SIGNATURES)

test/basicfuns.jl

@@ -135,21 +135,26 @@ end
     end
 
     @test isnan(logsubexp(Inf, Inf))
-    @test isnan(logsubexp(-Inf, -Inf))
+    @test logsubexp(-Inf, -Inf) ≡ -Inf
+    @test logsubexp(Inf, -Inf) ≡ Inf
+    @test logsubexp(-Inf, Inf) ≡ Inf
     @test logsubexp(Inf, 9.0) ≡ Inf
+    @test logsubexp(9.0, Inf) ≡ Inf
     @test logsubexp(-Inf, 9.0) ≡ 9.0
+    @test logsubexp(9.0, -Inf) ≡ 9.0
     @test logsubexp(1f2, 1f2) ≡ -Inf32
     @test logsubexp(0, 0) ≡ -Inf
     @test logsubexp(3, 2) ≈ 2.541324854612918108978
 
     # NaN propagation
-    @test isnan(logaddexp(NaN, 9.0))
-    @test isnan(logaddexp(NaN, Inf))
-    @test isnan(logaddexp(NaN, -Inf))
-
-    @test isnan(logsubexp(NaN, 9.0))
-    @test isnan(logsubexp(NaN, Inf))
-    @test isnan(logsubexp(NaN, -Inf))
+    for f in (logaddexp, logsubexp)
+        @test isnan(f(NaN, 9.0))
+        @test isnan(f(NaN, Inf))
+        @test isnan(f(NaN, -Inf))
+        @test isnan(f(9.0, NaN))
+        @test isnan(f(Inf, NaN))
+        @test isnan(f(-Inf, NaN))
+    end
 
     @test isnan(logsumexp([NaN, 9.0]))
     @test isnan(logsumexp([NaN, Inf]))
@@ -205,7 +210,7 @@ end
     s = softmax(x)
     @test s ≈ r
     @test eltype(s) === Float64
-    
+
     # non-standard indices: #12
     x = OffsetArray(1:3, -2:0)
     s = softmax(x)