From 19a2c9c6e078a4422ba7e1c7b9b29f6f35e79fdf Mon Sep 17 00:00:00 2001 From: "Tamas K. Papp" Date: Sun, 15 Sep 2019 09:24:08 +0200 Subject: [PATCH] Initial LogExpFunctions commit. This is a straightforward copy from StatsFuns.jl. --- .codecov.yml | 1 + .github/workflows/CI.yml | 70 +++++++++ .github/workflows/CompatHelper.yml | 24 ++++ .github/workflows/TagBot.yml | 17 +++ .gitignore | 8 ++ LICENSE.md | 23 +++ Project.toml | 16 +++ README.md | 10 ++ docs/Manifest.toml | 92 ++++++++++++ docs/Project.toml | 5 + docs/make.jl | 13 ++ docs/src/index.md | 23 +++ src/LogExpFunctions.jl | 221 +++++++++++------------------ src/logsumexp.jl | 96 +++++++++++++ test/coverage/Project.toml | 2 + test/coverage/coverage-summary.jl | 12 ++ test/coverage/coverage.jl | 9 ++ test/runtests.jl | 26 ++-- 18 files changed, 518 insertions(+), 150 deletions(-) create mode 100644 .codecov.yml create mode 100644 .github/workflows/CI.yml create mode 100644 .github/workflows/CompatHelper.yml create mode 100644 .github/workflows/TagBot.yml create mode 100644 .gitignore create mode 100644 LICENSE.md create mode 100644 Project.toml create mode 100644 README.md create mode 100644 docs/Manifest.toml create mode 100644 docs/Project.toml create mode 100644 docs/make.jl create mode 100644 docs/src/index.md create mode 100644 src/logsumexp.jl create mode 100644 test/coverage/Project.toml create mode 100644 test/coverage/coverage-summary.jl create mode 100644 test/coverage/coverage.jl diff --git a/.codecov.yml b/.codecov.yml new file mode 100644 index 00000000..69cb7601 --- /dev/null +++ b/.codecov.yml @@ -0,0 +1 @@ +comment: false diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 00000000..4296bf59 --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,70 @@ +# from https://discourse.julialang.org/t/easy-workflow-file-for-setting-up-github-actions-ci-for-your-julia-package/49765 + +name: CI +on: + pull_request: + branches: + - master + push: + branches: + - master + tags: '*' +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.0' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'. + - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia. + - 'nightly' + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: actions/cache@v1 + env: + cache-name: cache-artifacts + with: + path: ~/.julia/artifacts + key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} + restore-keys: | + ${{ runner.os }}-test-${{ env.cache-name }}- + ${{ runner.os }}-test- + ${{ runner.os }}- + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v1 + with: + file: lcov.info + docs: + name: Documentation + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: '1' + - run: | + julia --project=docs -e ' + using Pkg + Pkg.develop(PackageSpec(path=pwd())) + Pkg.instantiate()' + - run: | + julia --project=docs -e ' + using Documenter: doctest + using LogExpFunctions + doctest(LogExpFunctions)' + - run: julia --project=docs docs/make.jl + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml new file mode 100644 index 00000000..6698b3ce --- /dev/null +++ b/.github/workflows/CompatHelper.yml @@ -0,0 +1,24 @@ +# see the docs at https://github.com/JuliaRegistries/CompatHelper.jl + +name: CompatHelper + +on: + schedule: + - cron: '00 00 * * *' + workflow_dispatch: + +jobs: + CompatHelper: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1.2.0] + julia-arch: [x86] + os: [ubuntu-latest] + steps: + - name: Pkg.add("CompatHelper") + run: julia -e 'using Pkg; Pkg.add("CompatHelper")' + - name: CompatHelper.main() + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + run: julia -e 'using CompatHelper; CompatHelper.main()' diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml new file mode 100644 index 00000000..28f36cd3 --- /dev/null +++ b/.github/workflows/TagBot.yml @@ -0,0 +1,17 @@ +# see the docs at https://github.com/JuliaRegistries/TagBot + +name: TagBot +on: + issue_comment: + types: + - created + workflow_dispatch: +jobs: + TagBot: + if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' + runs-on: ubuntu-latest + steps: + - uses: JuliaRegistries/TagBot@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} + ssh: ${{ secrets.DOCUMENTER_KEY }} diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..df60e6c6 --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +*.jl.cov +*.jl.*.cov +*.jl.mem +/Manifest.toml +/deps/deps.jl +/docs/build +/docs/Manifest.toml +/test/coverage/Manifest.toml diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 00000000..793a3b54 --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,23 @@ +The LogExpFunctions.jl package is licensed under the MIT "Expat" License: + +> Original work Copyright (c) 2015: Dahua Lin, StatsFuns.jl contributors. +> Modified version Copyright (c) 2019: Tamas K. Papp. +> +> Permission is hereby granted, free of charge, to any person obtaining a copy +> of this software and associated documentation files (the "Software"), to deal +> in the Software without restriction, including without limitation the rights +> to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +> copies of the Software, and to permit persons to whom the Software is +> furnished to do so, subject to the following conditions: +> +> The above copyright notice and this permission notice shall be included in all +> copies or substantial portions of the Software. +> +> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +> IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +> FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +> AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +> SOFTWARE. +> diff --git a/Project.toml b/Project.toml new file mode 100644 index 00000000..b8d802bc --- /dev/null +++ b/Project.toml @@ -0,0 +1,16 @@ +name = "LogExpFunctions" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +authors = ["StatsFun.jl contributors, Tamas K. Papp "] +version = "0.1.1" + +[deps] +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" + +[compat] +julia = "^1" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/README.md b/README.md new file mode 100644 index 00000000..9e394aae --- /dev/null +++ b/README.md @@ -0,0 +1,10 @@ +# LogExpFunctions.jl + +![Lifecycle](https://img.shields.io/badge/lifecycle-experimental-orange.svg) +[![build](https://github.com/tpapp/LogExpFunctions.jl/workflows/CI/badge.svg)](https://github.com/tpapp/LogExpFunctions.jl/actions?query=workflow%3ACI) +[![codecov.io](http://codecov.io/github/tpapp/LogExpFunctions.jl/coverage.svg?branch=master)](http://codecov.io/github/tpapp/LogExpFunctions.jl?branch=master) +[![Documentation](https://img.shields.io/badge/docs-stable-blue.svg)](https://tpapp.github.io/LogExpFunctions.jl/stable) [![Documentation](https://img.shields.io/badge/docs-master-blue.svg)](https://tpapp.github.io/LogExpFunctions.jl/latest) + +Various special functions based on `log` and `exp` moved from [StatsFuns.jl](https://github.com/JuliaStats/StatsFuns.jl) into a separate package, to minimize dependencies. These functions only use native Julia code, so there is no need to depend on `librmath` or similar libraries. See the discussion at [StatsFuns.jl#46](https://github.com/JuliaStats/StatsFuns.jl/issues/46). + +The original authors of these functions are the StatsFuns.jl contributors. diff --git a/docs/Manifest.toml b/docs/Manifest.toml new file mode 100644 index 00000000..60252ac2 --- /dev/null +++ b/docs/Manifest.toml @@ -0,0 +1,92 @@ +# This file is machine-generated - editing it directly is not advised + +[[Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[Distributed]] +deps = ["Random", "Serialization", "Sockets"] +uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" + +[[DocStringExtensions]] +deps = ["LibGit2", "Markdown", "Pkg", "Test"] +git-tree-sha1 = "0513f1a8991e9d83255e0140aace0d0fc4486600" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.8.0" + +[[Documenter]] +deps = ["Base64", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] +git-tree-sha1 = "1b6ae3796f60311e39cd1770566140d2c056e87f" +uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +version = "0.23.3" + +[[InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[JSON]] +deps = ["Dates", "Mmap", "Parsers", "Unicode"] +git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e" +uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +version = "0.21.0" + +[[LibGit2]] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + +[[Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[Mmap]] +uuid = "a63ad114-7e13-5084-954f-fe012c677804" + +[[Parsers]] +deps = ["Dates", "Test"] +git-tree-sha1 = "ef0af6c8601db18c282d092ccbd2f01f3f0cd70b" +uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" +version = "0.3.7" + +[[Pkg]] +deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" + +[[Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[Random]] +deps = ["Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" + +[[Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[Test]] +deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 00000000..0d8fbc73 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,5 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" + +[compat] +Documenter = "~0.26" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 00000000..8f5c3000 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,13 @@ +using Documenter, LogExpFunctions + +makedocs( + modules = [LogExpFunctions], + format = Documenter.HTML(), + checkdocs = :exports, + sitename = "LogExpFunctions.jl", + pages = Any["index.md"] +) + +deploydocs( + repo = "github.com/tpapp/LogExpFunctions.jl.git", +) diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 00000000..ce3131b2 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,23 @@ +# LogExpFunctions + +Various special functions based on `log` and `exp` moved from [StatsFuns.jl](https://github.com/JuliaStats/StatsFuns.jl) into a separate package, to minimize dependencies. These functions only use native Julia code, so there is no need to depend on `librmath` or similar libraries. See the discussion at [StatsFuns.jl#46](https://github.com/JuliaStats/StatsFuns.jl/issues/46). + +The original authors of these functions are the StatsFuns.jl contributors. + +```@docs +xlogx +xlogy +logistic +logit +log1psq +log1pexp +log1mexp +log2mexp +logexpm1 +log1pmx +logmxp1 +logaddexp +logsumexp +softmax! +softmax +``` diff --git a/src/LogExpFunctions.jl b/src/LogExpFunctions.jl index 08348129..9b4f4938 100644 --- a/src/LogExpFunctions.jl +++ b/src/LogExpFunctions.jl @@ -1,13 +1,33 @@ -# common facilities +module LogExpFunctions + +export xlogx, xlogy, logistic, logit, log1psq, log1pexp, log1mexp, log2mexp, logexpm1, + log1pmx, logmxp1, logaddexp, logsubexp, logsumexp, softmax!, softmax + +using DocStringExtensions: SIGNATURES + +using Base: Math.@horner, @irrational + +#### +#### constants +#### + +@irrational loghalf -0.6931471805599453094 log(big(0.5)) +@irrational logtwo 0.6931471805599453094 log(big(2.)) +@irrational logπ 1.1447298858494001741 log(big(π)) +@irrational log2π 1.8378770664093454836 log(big(2.)*π) +@irrational log4π 2.5310242469692907930 log(big(4.)*π) + +#### +#### functions +#### -# scalar functions """ - xlogx(x::Number) +$(SIGNATURES) -Compute `x * log(x)`, returning zero if `x` is zero. +Return `x * log(x)` for `x ≥ 0`, handling ``x = 0`` by taking the downward limit. ```jldoctest -julia> StatsFuns.xlogx(0) +julia> xlogx(0) 0.0 ``` """ @@ -17,12 +37,12 @@ function xlogx(x::Number) end """ - xlogy(x::Number, y::Number) +$(SIGNATURES) -Compute `x * log(y)`, returning zero if `x` is zero. +Return `x * log(y)` for `y > 0` with correct limit at ``x = 0``. ```jldoctest -julia> StatsFuns.xlogy(0, 0) +julia> xlogy(0, 0) 0.0 ``` """ @@ -31,6 +51,20 @@ function xlogy(x::Number, y::Number) ifelse(iszero(x) && !isnan(y), zero(result), result) end +""" +$(SIGNATURES) + +The [logistic](https://en.wikipedia.org/wiki/Logistic_function) sigmoid function mapping a +real number to a value in the interval ``[0,1]``, + +```math +\\sigma(x) = \\frac{1}{e^{-x} + 1} = \\frac{e^x}{1+e^x}. +``` + +Its inverse is the [`logit`](@ref) function. +""" +logistic(x::Real) = inv(exp(-x) + one(x)) + # 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. @@ -41,21 +75,9 @@ end # ) # end -@inline _logistic_bounds(x::Float16) = (Float16(-16.64), Float16(7.625)) -@inline _logistic_bounds(x::Float32) = (-103.27893f0, 16.635532f0) -@inline _logistic_bounds(x::Float64) = (-744.4400719213812, 36.7368005696771) - -""" - logistic(x::Real) - -The [logistic](https://en.wikipedia.org/wiki/Logistic_function) sigmoid function mapping a real number to a value in the interval [0,1], -```math -\\sigma(x) = \\frac{1}{e^{-x} + 1} = \\frac{e^x}{1+e^x}. -``` - -Its inverse is the [`logit`](@ref) function. -""" -logistic(x::Real) = inv(exp(-x) + one(x)) +@inline _logistic_bounds(::Float16) = (Float16(-16.64), Float16(7.625)) +@inline _logistic_bounds(::Float32) = (-103.27893f0, 16.635532f0) +@inline _logistic_bounds(::Float64) = (-744.4400719213812, 36.7368005696771) function logistic(x::Union{Float16, Float32, Float64}) e = exp(x) @@ -72,18 +94,20 @@ function logistic(x::Union{Float16, Float32, Float64}) end """ - logit(p::Real) +$(SIGNATURES) The [logit](https://en.wikipedia.org/wiki/Logit) or log-odds transformation, + ```math \\log\\left(\\frac{x}{1-x}\\right), \\text{where} 0 < x < 1 ``` + Its inverse is the [`logistic`](@ref) function. """ logit(x::Real) = log(x / (one(x) - x)) """ - log1psq(x::Real) +$(SIGNATURES) Return `log(1+x^2)` evaluated carefully for `abs(x)` very small or very large. """ @@ -94,7 +118,7 @@ function log1psq(x::Union{Float32,Float64}) end """ - log1pexp(x::Real) +$(SIGNATURES) Return `log(1+exp(x))` evaluated carefully for largish `x`. @@ -105,41 +129,41 @@ log1pexp(x::Real) = x < 18.0 ? log1p(exp(x)) : x < 33.3 ? x + exp(-x) : oftype(e log1pexp(x::Float32) = x < 9.0f0 ? log1p(exp(x)) : x < 16.0f0 ? x + exp(-x) : oftype(exp(-x), x) """ - log1mexp(x::Real) +$(SIGNATURES) Return `log(1 - exp(x))` See: - * Martin Maechler (2012) "Accurately Computing log(1 − exp(− |a|))", - http://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf + * Martin Maechler (2012) [“Accurately Computing log(1 − exp(− |a|))”](http://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf) Note: different than Maechler (2012), no negation inside parentheses """ log1mexp(x::Real) = x < loghalf ? log1p(-exp(x)) : log(-expm1(x)) """ - log2mexp(x::Real) +$(SIGNATURES) Return `log(2 - exp(x))` evaluated as `log1p(-expm1(x))` """ log2mexp(x::Real) = log1p(-expm1(x)) """ - logexpm1(x::Real) +$(SIGNATURES) -Return `log(exp(x) - 1)` or the "invsoftplus" function. -It is the inverse of [`log1pexp`](@ref) (aka "softplus"). +Return `log(exp(x) - 1)` or the “invsoftplus” function. It is the inverse of +[`log1pexp`](@ref) (aka “softplus”). """ logexpm1(x::Real) = x <= 18.0 ? log(expm1(x)) : x <= 33.3 ? x - exp(-x) : oftype(exp(-x), x) + logexpm1(x::Float32) = x <= 9f0 ? log(expm1(x)) : x <= 16f0 ? x - exp(-x) : oftype(exp(-x), x) const softplus = log1pexp const invsoftplus = logexpm1 """ - log1pmx(x::Float64) +$(SIGNATURES) -Return `log(1 + x) - x` +Return `log(1 + x) - x`. Use naive calculation or range reduction outside kernel range. Accurate ~2ulps for all `x`. """ @@ -164,7 +188,7 @@ function log1pmx(x::Float64) end """ - logmxp1(x::Float64) +$(SIGNATURES) Return `log(x) - x + 1` carefully evaluated. """ @@ -200,11 +224,11 @@ function _log1pmx_ker(x::Float64) r*(hxsq+w*t)-hxsq end - """ - logaddexp(x::Real, y::Real) +$(SIGNATURES) -Return `log(exp(x) + exp(y))`, avoiding intermediate overflow/undeflow, and handling non-finite values. +Return `log(exp(x) + exp(y))`, avoiding intermediate overflow/undeflow, and handling +non-finite values. """ function logaddexp(x::Real, y::Real) # ensure Δ = 0 if x = y = ± Inf @@ -215,111 +239,14 @@ end Base.@deprecate logsumexp(x::Real, y::Real) logaddexp(x, y) """ - logsubexp(x, y) +$(SIGNATURES) -Return `log(abs(e^x - e^y))`, preserving numerical accuracy. +Return `log(abs(exp(x) - exp(y)))`, preserving numerical accuracy. """ logsubexp(x::Real, y::Real) = max(x, y) + log1mexp(-abs(x - y)) """ - logsumexp(X) - -Compute `log(sum(exp, X))` in a numerically stable way that avoids intermediate over- and -underflow. - -`X` should be an iterator of real numbers. The result is computed using a single pass over -the data. - -# References - -[Sebastian Nowozin: Streaming Log-sum-exp Computation.](http://www.nowozin.net/sebastian/blog/streaming-log-sum-exp-computation.html) -""" -logsumexp(X) = _logsumexp_onepass(X) - -""" - logsumexp(X::AbstractArray{<:Real}; dims=:) - -Compute `log.(sum(exp.(X); dims=dims))` in a numerically stable way that avoids -intermediate over- and underflow. - -The result is computed using a single pass over the data. - -# References - -[Sebastian Nowozin: Streaming Log-sum-exp Computation.](http://www.nowozin.net/sebastian/blog/streaming-log-sum-exp-computation.html) -""" -logsumexp(X::AbstractArray{<:Real}; dims=:) = _logsumexp(X, dims) - -_logsumexp(X::AbstractArray{<:Real}, ::Colon) = _logsumexp_onepass(X) -function _logsumexp(X::AbstractArray{<:Real}, dims) - # Do not use log(zero(eltype(X))) directly to avoid issues with ForwardDiff (#82) - FT = float(eltype(X)) - xmax_r = reduce(_logsumexp_onepass_op, X; dims=dims, init=(FT(-Inf), zero(FT))) - return @. first(xmax_r) + log1p(last(xmax_r)) -end - -function _logsumexp_onepass(X) - # fallback for empty collections - isempty(X) && return log(sum(X)) - return _logsumexp_onepass_result(_logsumexp_onepass_reduce(X, Base.IteratorEltype(X))) -end - -# function barrier for reductions with single element and without initial element -_logsumexp_onepass_result(x) = float(x) -_logsumexp_onepass_result((xmax, r)::Tuple) = xmax + log1p(r) - -# iterables with known element type -function _logsumexp_onepass_reduce(X, ::Base.HasEltype) - # do not perform type computations if element type is abstract - T = eltype(X) - isconcretetype(T) || return _logsumexp_onepass_reduce(X, Base.EltypeUnknown()) - - FT = float(T) - return reduce(_logsumexp_onepass_op, X; init=(FT(-Inf), zero(FT))) -end - -# iterables without known element type -_logsumexp_onepass_reduce(X, ::Base.EltypeUnknown) = reduce(_logsumexp_onepass_op, X) - -## Reductions for one-pass algorithm: avoid expensive multiplications if numbers are reduced - -# reduce two numbers -function _logsumexp_onepass_op(x1, x2) - a = x1 == x2 ? zero(x1 - x2) : -abs(x1 - x2) - xmax = x1 > x2 ? oftype(a, x1) : oftype(a, x2) - r = exp(a) - return xmax, r -end - -# reduce a number and a partial sum -function _logsumexp_onepass_op(x, (xmax, r)::Tuple) - a = x == xmax ? zero(x - xmax) : -abs(x - xmax) - if x > xmax - _xmax = oftype(a, x) - _r = (r + one(r)) * exp(a) - else - _xmax = oftype(a, xmax) - _r = r + exp(a) - end - return _xmax, _r -end -_logsumexp_onepass_op(xmax_r::Tuple, x) = _logsumexp_onepass_op(x, xmax_r) - -# reduce two partial sums -function _logsumexp_onepass_op((xmax1, r1)::Tuple, (xmax2, r2)::Tuple) - a = xmax1 == xmax2 ? zero(xmax1 - xmax2) : -abs(xmax1 - xmax2) - if xmax1 > xmax2 - xmax = oftype(a, xmax1) - r = r1 + (r2 + one(r2)) * exp(a) - else - xmax = oftype(a, xmax2) - r = r2 + (r1 + one(r1)) * exp(a) - end - return xmax, r -end - -""" - softmax!(r::AbstractArray, x::AbstractArray) +$(SIGNATURES) Overwrite `r` with the `softmax` (or _normalized exponential_) transformation of `x` @@ -343,9 +270,21 @@ function softmax!(r::AbstractArray{R}, x::AbstractArray{T}) where {R<:AbstractFl end """ - softmax(x::AbstractArray{<:Real}) +$(SIGNATURES) -Return the [`softmax transformation`](https://en.wikipedia.org/wiki/Softmax_function) applied to `x` +Return the [`softmax transformation`](https://en.wikipedia.org/wiki/Softmax_function) +applied to `x` *in place*. """ softmax!(x::AbstractArray{<:AbstractFloat}) = softmax!(x, x) + +""" +$(SIGNATURES) + +Return the [`softmax transformation`](https://en.wikipedia.org/wiki/Softmax_function) +applied to `x`. +""" softmax(x::AbstractArray{<:Real}) = softmax!(similar(x, Float64), x) + +include("logsumexp.jl") + +end # module diff --git a/src/logsumexp.jl b/src/logsumexp.jl new file mode 100644 index 00000000..5f525ff1 --- /dev/null +++ b/src/logsumexp.jl @@ -0,0 +1,96 @@ +""" +$(SIGNATURES) + +Compute `log(sum(exp, X))` in a numerically stable way that avoids intermediate over- and +underflow. + +`X` should be an iterator of real numbers. The result is computed using a single pass over +the data. + +# References + +[Sebastian Nowozin: Streaming Log-sum-exp Computation](http://www.nowozin.net/sebastian/blog/streaming-log-sum-exp-computation.html) +""" +logsumexp(X) = _logsumexp_onepass(X) + +""" +$(SIGNATURES) + +Compute `log.(sum(exp.(X); dims=dims))` in a numerically stable way that avoids +intermediate over- and underflow. + +The result is computed using a single pass over the data. + +# References + +[Sebastian Nowozin: Streaming Log-sum-exp Computation](http://www.nowozin.net/sebastian/blog/streaming-log-sum-exp-computation.html) +""" +logsumexp(X::AbstractArray{<:Real}; dims=:) = _logsumexp(X, dims) + +_logsumexp(X::AbstractArray{<:Real}, ::Colon) = _logsumexp_onepass(X) +function _logsumexp(X::AbstractArray{<:Real}, dims) + # Do not use log(zero(eltype(X))) directly to avoid issues with ForwardDiff (#82) + FT = float(eltype(X)) + xmax_r = reduce(_logsumexp_onepass_op, X; dims=dims, init=(FT(-Inf), zero(FT))) + return @. first(xmax_r) + log1p(last(xmax_r)) +end + +function _logsumexp_onepass(X) + # fallback for empty collections + isempty(X) && return log(sum(X)) + return _logsumexp_onepass_result(_logsumexp_onepass_reduce(X, Base.IteratorEltype(X))) +end + +# function barrier for reductions with single element and without initial element +_logsumexp_onepass_result(x) = float(x) +_logsumexp_onepass_result((xmax, r)::Tuple) = xmax + log1p(r) + +# iterables with known element type +function _logsumexp_onepass_reduce(X, ::Base.HasEltype) + # do not perform type computations if element type is abstract + T = eltype(X) + isconcretetype(T) || return _logsumexp_onepass_reduce(X, Base.EltypeUnknown()) + + FT = float(T) + return reduce(_logsumexp_onepass_op, X; init=(FT(-Inf), zero(FT))) +end + +# iterables without known element type +_logsumexp_onepass_reduce(X, ::Base.EltypeUnknown) = reduce(_logsumexp_onepass_op, X) + +## Reductions for one-pass algorithm: avoid expensive multiplications if numbers are reduced + +# reduce two numbers +function _logsumexp_onepass_op(x1, x2) + a = x1 == x2 ? zero(x1 - x2) : -abs(x1 - x2) + xmax = x1 > x2 ? oftype(a, x1) : oftype(a, x2) + r = exp(a) + return xmax, r +end + +# reduce a number and a partial sum +function _logsumexp_onepass_op(x, (xmax, r)::Tuple) + a = x == xmax ? zero(x - xmax) : -abs(x - xmax) + if x > xmax + _xmax = oftype(a, x) + _r = (r + one(r)) * exp(a) + else + _xmax = oftype(a, xmax) + _r = r + exp(a) + end + return _xmax, _r +end +_logsumexp_onepass_op(xmax_r::Tuple, x) = _logsumexp_onepass_op(x, xmax_r) + +# reduce two partial sums +function _logsumexp_onepass_op((xmax1, r1)::Tuple, (xmax2, r2)::Tuple) + a = xmax1 == xmax2 ? zero(xmax1 - xmax2) : -abs(xmax1 - xmax2) + if xmax1 > xmax2 + xmax = oftype(a, xmax1) + r = r1 + (r2 + one(r2)) * exp(a) + else + xmax = oftype(a, xmax2) + r = r2 + (r1 + one(r1)) * exp(a) + end + return xmax, r +end diff --git a/test/coverage/Project.toml b/test/coverage/Project.toml new file mode 100644 index 00000000..4fbdc477 --- /dev/null +++ b/test/coverage/Project.toml @@ -0,0 +1,2 @@ +[deps] +Coverage = "a2441757-f6aa-5fb2-8edb-039e3f45d037" diff --git a/test/coverage/coverage-summary.jl b/test/coverage/coverage-summary.jl new file mode 100644 index 00000000..36720c3d --- /dev/null +++ b/test/coverage/coverage-summary.jl @@ -0,0 +1,12 @@ +#### +#### Coverage summary, printed as "(percentage) covered". +#### +#### Useful for CI environments that just want a summary (eg a Gitlab setup). +#### + +using Coverage +cd(joinpath(@__DIR__, "..", "..")) do + covered_lines, total_lines = get_summary(process_folder()) + percentage = covered_lines / total_lines * 100 + println("($(percentage)%) covered") +end diff --git a/test/coverage/coverage.jl b/test/coverage/coverage.jl new file mode 100644 index 00000000..de6fe46a --- /dev/null +++ b/test/coverage/coverage.jl @@ -0,0 +1,9 @@ +# only push coverage from one bot +get(ENV, "TRAVIS_OS_NAME", nothing) == "linux" || exit(0) +get(ENV, "TRAVIS_JULIA_VERSION", nothing) == "1.1" || exit(0) + +using Coverage + +cd(joinpath(@__DIR__, "..", "..")) do + Codecov.submit(Codecov.process_folder()) +end diff --git a/test/runtests.jl b/test/runtests.jl index aff76fa4..e212e4a9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,4 @@ -using StatsFuns, Test +using LogExpFunctions, Test @testset "xlogx & xlogy" begin @test iszero(xlogx(0)) @@ -8,6 +8,7 @@ using StatsFuns, Test @test iszero(xlogy(0, 1)) @test xlogy(2, 3) ≈ 2.0 * log(3.0) + @test xlogy(2, 3.0) ≈ 2.0 * log(3.0) # promoted @test_throws DomainError xlogy(1, -1) @test isnan(xlogy(NaN, 2)) @test isnan(xlogy(2, NaN)) @@ -44,6 +45,7 @@ end @test iszero(log1psq(0.0)) @test log1psq(1.0) ≈ log1p(1.0) @test log1psq(2.0) ≈ log1p(4.0) + @test log1psq(Float16(0.1)) ≈ log1p(0.01) end # log1pexp, log1mexp, log2mexp & logexpm1 @@ -82,25 +84,30 @@ end @testset "log1pmx" begin @test iszero(log1pmx(0.0)) - @test log1pmx(1.0) ≈ log(2.0) - 1.0 - @test log1pmx(2.0) ≈ log(3.0) - 2.0 + for x in [-0.65, -0.5, -0.3, -0.1, 0.4, 1.0, 2.0] + z = BigFloat(x) + @test log1pmx(x) ≈ Float64(log(1 + z) - z) + end end @testset "logmxp1" begin @test iszero(logmxp1(1.0)) - @test logmxp1(2.0) ≈ log(2.0) - 1.0 - @test logmxp1(3.0) ≈ log(3.0) - 2.0 + for x in [0.1, 0.35, 0.5, 2.0, 3.0] + z = BigFloat(x) + @test logmxp1(x) ≈ Float64(log(z) - z + 1) + end end @testset "logsumexp" begin @test logaddexp(2.0, 3.0) ≈ log(exp(2.0) + exp(3.0)) @test logaddexp(10002, 10003) ≈ 10000 + logaddexp(2.0, 3.0) + @test logaddexp(2, 3.0) ≈ log(exp(2.0) + exp(3.0)) # promotion - @test @inferred(logsumexp([1.0])) == 1.0 - @test @inferred(logsumexp((x for x in [1.0]))) == 1.0 + @test @inferred(logsumexp([1.0])) == 1 + @test @inferred(logsumexp((x for x in [1.0]))) == 1 @test @inferred(logsumexp([1.0, 2.0, 3.0])) ≈ 3.40760596444438 @test @inferred(logsumexp((1.0, 2.0, 3.0))) ≈ 3.40760596444438 - @test logsumexp([1.0, 2.0, 3.0] .+ 1000.) ≈ 1003.40760596444438 + @test logsumexp([1.0, 2.0, 3.0]) .+ 1000. ≈ 1003.40760596444438 @test @inferred(logsumexp([[1.0, 2.0, 3.0] [1.0, 2.0, 3.0] .+ 1000.]; dims=1)) ≈ [3.40760596444438 1003.40760596444438] @test @inferred(logsumexp([[1.0 2.0 3.0]; [1.0 2.0 3.0] .+ 1000.]; dims=2)) ≈ [3.40760596444438, 1003.40760596444438] @@ -143,7 +150,8 @@ end @test isnan(logsumexp([NaN, Inf])) @test isnan(logsumexp([NaN, -Inf])) - # logsumexp with general iterables (issue #63) + # logsumexp with general iterables + # https://github.com/JuliaStats/StatsFuns.jl/issues/63 xs = range(-500, stop = 10, length = 1000) @test @inferred(logsumexp(x for x in xs)) == logsumexp(xs) end