Add loglogistic, logitexp, log1mlogistic and logit1mexp (#82)

* Add `loglogistic`, `logitexp`, `log1mlogistic` and `logit1mexp`

This takes advantage of `LogExpFunctions`'s accurate
implementations of `log1pexp` and `log1mexp`, combined with negation
in the log-odds domain to provide more accurate and less expensive
implementations of the function compositions.

* Register inverses; why not?

* Update documentation

* Fix failing test: be slightly less stringent

* Re-style documentation

* Use alternate design for `loglogistic`

* Simplify dispatch cases

* Add `ChangesOfVariables` support

* Fix duplication in docstring

* Be as careful as possible in the evaluation; we get ulp error of `logexpm1`

* Apply the same conversion technique as `loglogistic`

i.e. properly handle `Integer` and `Rational`

* Eliminate 3 lines of code

* Hide the proper handling of integers and rationals.

This enables the desired interface of `f(x::Real)`

* Add tests of inverses

* Add ChainRules extension and respective tests

* Add ChangesOfVariables tests

* Fix typo

* Update documentation

* Add tests to ensure correctness for subtypes of `Real` in `Base`

* Simplify definitions

* Remove ChainRules definitions

* Also remove ChainRules tests

* Match expected form of `ChangesOfVariables` tests

* Fix oversight in function type dispatches

* Add `ChainRulesCore` support back

* Revert "Add `ChainRulesCore` support back"

This reverts commit 1b0f914.

* Add extensive return type tests

* Fix minor oversight from running tests locally

* Add tests for `NaN`

Author: andrewjradcliffe

docs/src/index.md

@@ -4,7 +4,7 @@ Various special functions based on `log` and `exp` moved from [StatsFuns.jl](htt
 
 The original authors of these functions are the StatsFuns.jl contributors.
 
-LogExpFunctions supports [`InverseFunctions.inverse`](https://github.com/JuliaMath/InverseFunctions.jl) and [`ChangesOfVariables.test_with_logabsdet_jacobian`](https://github.com/JuliaMath/ChangesOfVariables.jl) for `log1mexp`, `log1pexp`, `log2mexp`, `logexpm1`, `logistic`, `logit`, and `logcosh` (no inverse).
+LogExpFunctions supports [`InverseFunctions.inverse`](https://github.com/JuliaMath/InverseFunctions.jl) and [`ChangesOfVariables.test_with_logabsdet_jacobian`](https://github.com/JuliaMath/ChangesOfVariables.jl) for `log1mexp`, `log1pexp`, `log2mexp`, `logexpm1`, `logistic`, `logit`, `loglogistic`, `logitexp`, `log1mlogistic`, `logit1mexp`, and `logcosh` (no inverse).
 
 ```@docs
 xlogx
@@ -31,4 +31,8 @@ softmax!
 softmax
 cloglog
 cexpexp
+loglogistic
+logitexp
+log1mlogistic
+logit1mexp
 ```

ext/LogExpFunctionsChangesOfVariablesExt.jl

@@ -42,4 +42,24 @@ function ChangesOfVariables.with_logabsdet_jacobian(::typeof(logcosh), x::Real)
     return y, log1mexp(a) - z
 end
 
+function ChangesOfVariables.with_logabsdet_jacobian(::typeof(loglogistic), x::Real)
+    y = loglogistic(x)
+    return y, y - x
+end
+
+function ChangesOfVariables.with_logabsdet_jacobian(::typeof(log1mlogistic), x::Real)
+    y = log1mlogistic(x)
+    return y, x + y
+end
+
+function ChangesOfVariables.with_logabsdet_jacobian(::typeof(logitexp), x::Real)
+    y = logitexp(x)
+    return y, y - x
+end
+
+function ChangesOfVariables.with_logabsdet_jacobian(::typeof(logit1mexp), x::Real)
+    y = logit1mexp(x)
+    return y, -y - x
+end
+
 end # module

ext/LogExpFunctionsInverseFunctionsExt.jl

@@ -16,4 +16,10 @@ InverseFunctions.inverse(::typeof(logistic)) = logit
 InverseFunctions.inverse(::typeof(cloglog)) = cexpexp
 InverseFunctions.inverse(::typeof(cexpexp)) = cloglog
 
+InverseFunctions.inverse(::typeof(loglogistic)) = logitexp
+InverseFunctions.inverse(::typeof(logitexp)) = loglogistic
+
+InverseFunctions.inverse(::typeof(log1mlogistic)) = logit1mexp
+InverseFunctions.inverse(::typeof(logit1mexp)) = log1mlogistic
+
 end # module

src/LogExpFunctions.jl

@@ -8,7 +8,8 @@ import LinearAlgebra
 
 export xlogx, xlogy, xlog1py, xexpx, xexpy, logistic, logit, log1psq, log1pexp, log1mexp, log2mexp, logexpm1,
     softplus, invsoftplus, log1pmx, logmxp1, logaddexp, logsubexp, logsumexp, logsumexp!, softmax,
-    softmax!, logcosh, logabssinh, cloglog, cexpexp
+    softmax!, logcosh, logabssinh, cloglog, cexpexp,
+    loglogistic, logitexp, log1mlogistic, logit1mexp
 
 # expm1(::Float16) is not defined in older Julia versions,
 # hence for better Float16 support we use an internal function instead

src/basicfuns.jl

@@ -441,3 +441,68 @@ $(SIGNATURES)
 Compute the complementary double exponential, `1 - exp(-exp(x))`.
 """
 cexpexp(x) = -_expm1(-exp(x))
+
+#=
+this uses the identity:
+
+log(logistic(x)) = -log(1 + exp(-x))
+=#
+"""
+$(SIGNATURES)
+
+Return `log(logistic(x))`, computed more carefully and with fewer calls
+than the naive composition of functions.
+
+Its inverse is the [`logitexp`](@ref) function.
+"""
+loglogistic(x::Real) = -log1pexp(-float(x))
+
+#=
+this uses the identity:
+
+logit(exp(x)) = log(exp(x) / (1 + exp(x))) = -log(exp(-x) - 1)
+=#
+"""
+$(SIGNATURES)
+
+Return `logit(exp(x))`, computed more carefully and with fewer calls than
+the naive composition of functions.
+
+Its inverse is the [`loglogistic`](@ref) function.
+"""
+logitexp(x::Real) = -logexpm1(-float(x))
+
+#=
+this uses the identity:
+
+log(logistic(-x)) = -log(1 + exp(x))
+
+that is, negation in the log-odds domain.
+=#
+
+"""
+$(SIGNATURES)
+
+Return `log(1 - logistic(x))`, computed more carefully and with fewer calls than
+the naive composition of functions.
+
+Its inverse is the [`logit1mexp`](@ref) function.
+"""
+log1mlogistic(x::Real) = -log1pexp(x)
+
+#=
+
+this uses the same identity:
+
+-logit(exp(x)) = logit(1 - exp(x)) = log((1 - exp(x)) / exp(x)) = log(exp(-x) - 1)
+=#
+
+"""
+$(SIGNATURES)
+
+Return `logit(1 - exp(x))`, computed more carefully and with fewer calls than
+the naive composition of functions.
+
+Its inverse is the [`log1mlogistic`](@ref) function.
+"""
+logit1mexp(x::Real) = logexpm1(-float(x))

test/basicfuns.jl

@@ -468,3 +468,161 @@ end
     @test cexpexp(-Inf) == 0.0
     @test cexpexp(0) == (ℯ - 1) / ℯ
 end
+
+@testset "loglogistic: $T" for T in (Float16, Float32, Float64)
+    lim1 = T === Float16 ? -14.0 : -50.0
+    lim2 = T === Float16 ? -10.0 : -37.0
+    xs = T[Inf, -Inf, 0.0, lim1, lim2]
+    for x in xs
+        @test loglogistic(x) == log(logistic(x))
+    end
+
+    ϵ = eps(T)
+    xs = T[ϵ, 1.0, 18.0, 33.3, 50.0]
+    for x in xs
+        lhs = loglogistic(x)
+        rhs = log(logistic(x))
+        @test abs(lhs - rhs) < ϵ
+    end
+
+    # misc
+    @test loglogistic(T(Inf)) == -zero(T)
+    @test loglogistic(-T(Inf)) == -T(Inf)
+    @test loglogistic(-T(103.0)) == -T(103.0)
+    @test abs(loglogistic(T(35.0))) < 3eps(T)
+    @test abs(loglogistic(T(103.0))) < eps(T)
+    @test isfinite(loglogistic(-T(745.0)))
+    @test isfinite(loglogistic(T(50.0)))
+    @test isfinite(loglogistic(T(745.0)))
+    @test isnan(loglogistic(T(NaN)))
+    @test isnan(loglogistic(T(-NaN)))
+
+    # type-consistency
+    xs = T[Inf, -Inf, 0.0, lim1, lim2, ϵ, 1.0, 18.0, 33.3, 50.0]
+    for x in xs
+        @test typeof(loglogistic(x)) === T
+    end
+end
+
+
+@testset "logitexp: $T" for T in (Float16, Float32, Float64)
+    ϵ = eps(T)
+    xs = T[ϵ, √ϵ, 0.2, 0.4, 0.8, 1.0 - √ϵ, 1.0 - ϵ]
+    neg_xs = -xs
+    for x in xs
+        @test abs(logitexp(loglogistic(x)) - x) <= ϵ
+    end
+    for x in neg_xs
+        @test abs(logitexp(loglogistic(x)) - x) <= 2ϵ
+    end
+    xs = T[-Inf, 0.0, Inf]
+    for x in xs
+        @test logitexp(loglogistic(x)) == x
+    end
+
+    # misc
+    @test isnan(logitexp(T(NaN)))
+    @test isnan(logitexp(T(-NaN)))
+
+    # type-consistency
+    xs = loglogistic.([T[ϵ, √ϵ, 0.2, 0.4, 0.8, 1.0 - √ϵ, 1.0 - ϵ]; neg_xs; T[-Inf, 0.0, Inf]])
+    for x in xs
+        @test typeof(logitexp(x)) === T
+    end
+end
+
+@testset "log1mlogistic: $T" for T in (Float16, Float32, Float64)
+    @test log1mlogistic(T(Inf)) == -T(Inf)
+    @test log1mlogistic(-T(Inf)) == -zero(T)
+    @test log1mlogistic(-T(103.0)) < eps(T)
+    @test abs(log1mlogistic(T(35.0))) == T(35.0)
+    @test abs(log1mlogistic(T(103.0))) == T(103.0)
+    @test isfinite(log1mlogistic(-T(745.0)))
+    @test isfinite(log1mlogistic(T(50.0)))
+    @test isfinite(log1mlogistic(T(745.0)))
+    @test isnan(log1mlogistic(T(NaN)))
+    @test isnan(log1mlogistic(T(-NaN)))
+
+    # type-consistency
+    ϵ = eps(T)
+    lim1 = T === Float16 ? -14.0 : -50.0
+    lim2 = T === Float16 ? -10.0 : -37.0
+    xs = T[Inf, -Inf, 0.0, lim1, lim2, ϵ, 1.0, 18.0, 33.3, 50.0]
+    for x in xs
+        @test typeof(log1mlogistic(x)) === T
+    end
+end
+
+
+@testset "logit1mexp: $T" for T in (Float16, Float32, Float64)
+    ϵ = eps(T)
+    xs = T[ϵ, √ϵ, 0.2, 0.4, 0.8, 1.0 - √ϵ, 1.0 - ϵ]
+    neg_xs = -xs
+    for x in xs
+        @test abs(logit1mexp(log1mlogistic(x)) - x) <= 2ϵ
+    end
+    for x in neg_xs
+        @test abs(logit1mexp(log1mlogistic(x)) - x) <= ϵ
+    end
+    xs = T[-Inf, 0.0, Inf]
+    for x in xs
+        @test logit1mexp(log1mlogistic(x)) == x
+    end
+
+    # misc
+    @test isnan(logit1mexp(T(NaN)))
+    @test isnan(logit1mexp(T(-NaN)))
+
+    # type-consistency
+    xs = log1mlogistic.([T[ϵ, √ϵ, 0.2, 0.4, 0.8, 1.0 - √ϵ, 1.0 - ϵ]; neg_xs; T[-Inf, 0.0, Inf]])
+    for x in xs
+        @test typeof(logit1mexp(x)) === T
+    end
+end
+
+@testset "loglogistic, log1mlogistic correctness wrt Unsigned, Rational: $T" for T in
+    (UInt8, UInt16, UInt32, UInt64, UInt128)
+    @test loglogistic(T(5)) == loglogistic(5.0)
+    @test log1mlogistic(T(5)) == log1mlogistic(5.0)
+
+    @test loglogistic(T(0x01)//T(0x02)) == loglogistic(0.5)
+    @test log1mlogistic(T(0x01)//T(0x02)) == log1mlogistic(0.5)
+
+    @test loglogistic(T(0x01)//T(0x00)) == loglogistic(Inf)
+    @test log1mlogistic(T(0x01)//T(0x00)) == log1mlogistic(Inf)
+
+end
+
+@testset "loglogistic, log1mlogistic, logitexp, logit1mexp correctness wrt Integer edge case: $T" for T in
+    (Int8, Int16, Int32, Int64, Int128)
+    # If not handled, these will be zero
+    @test loglogistic(typemin(T)) == loglogistic(float(typemin(T)))
+    @test log1mlogistic(typemin(T)) == log1mlogistic(float(typemin(T)))
+
+    # If not handled, these would throw since negation at typemin is a round-trip
+    @test logitexp(typemin(T)) == logitexp(float(typemin(T)))
+    @test logit1mexp(typemin(T)) == logit1mexp(float(typemin(T)))
+end
+
+@testset "loglogistic, log1mlogistic, logitexp, logit1mexp return types" begin
+    for T in
+        (UInt8, UInt16, UInt32, UInt64, UInt128, Int8, Int16, Int32, Int64, Int128)
+        @test typeof(loglogistic(T(5))) === Float64
+        @test typeof(log1mlogistic(T(5))) === Float64
+        @test typeof(loglogistic(T(0x01)//T(0x02))) === Float64
+        @test typeof(log1mlogistic(T(0x01)//T(0x02))) === Float64
+        @test typeof(logitexp(T(0))) === Float64
+        @test typeof(logit1mexp(T(0))) === Float64
+    end
+    @test typeof(loglogistic(BigInt(5))) === BigFloat
+    @test typeof(log1mlogistic(BigInt(5))) === BigFloat
+    @test typeof(loglogistic(BigInt(0x01)//BigInt(0x02))) === BigFloat
+    @test typeof(log1mlogistic(BigInt(0x01)//BigInt(0x02))) === BigFloat
+    @test typeof(logitexp(BigInt(0))) === BigFloat
+    @test typeof(logit1mexp(BigInt(0))) === BigFloat
+
+    @test typeof(loglogistic(BigFloat(5))) === BigFloat
+    @test typeof(log1mlogistic(BigFloat(5))) === BigFloat
+    @test typeof(logitexp(BigFloat(0))) === BigFloat
+    @test typeof(logit1mexp(BigFloat(0))) === BigFloat
+end

test/inverse.jl

@@ -11,4 +11,10 @@
 
     InverseFunctions.test_inverse(cloglog, rand())
     InverseFunctions.test_inverse(cexpexp, rand())
+
+    InverseFunctions.test_inverse(loglogistic, randexp())
+    InverseFunctions.test_inverse(logitexp, -randexp())
+
+    InverseFunctions.test_inverse(log1mlogistic, randexp())
+    InverseFunctions.test_inverse(logit1mexp, -randexp())
 end

test/with_logabsdet_jacobian.jl

@@ -19,4 +19,19 @@
 
     ChangesOfVariables.test_with_logabsdet_jacobian(logcosh, x, derivative)
     ChangesOfVariables.test_with_logabsdet_jacobian(logcosh, -x, derivative)
+
+    dloglogistic(x) = logistic(-x)
+    dlog1mlogistic(x) = -logistic(x)
+    dlogitexp(x) = inv(1 - exp(x))
+    dlogit1mexp(x) = -inv(1 - exp(x))
+    derivative(f::typeof(loglogistic), x) = dloglogistic(x)
+    derivative(f::typeof(log1mlogistic), x) = dlog1mlogistic(x)
+    derivative(f::typeof(logitexp), x) = dlogitexp(x)
+    derivative(f::typeof(logit1mexp), x) = dlogit1mexp(x)
+
+    ChangesOfVariables.test_with_logabsdet_jacobian(loglogistic, x, derivative)
+    ChangesOfVariables.test_with_logabsdet_jacobian(logitexp, -x, derivative)
+
+    ChangesOfVariables.test_with_logabsdet_jacobian(log1mlogistic, x, derivative)
+    ChangesOfVariables.test_with_logabsdet_jacobian(logit1mexp, -x, derivative)
 end