R Under development (unstable) (2023-11-06 r85483 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("matrixStats")
> library("stats")
> 
> logSumExp_R <- function(lx, na.rm = FALSE) {
+   log(sum(exp(lx), na.rm = na.rm))
+ }
> 
> ## R-help thread \emph{'[R] Beyond double-precision?'} on May 9, 2009.
> 
> for (mode in c("integer", "double")) {
+   cat("mode: ", mode, "\n", sep = "")
+ 
+   set.seed(1)
+   x <- runif(20, min = 1.0, max = 3.0)
+   storage.mode(x) <- mode
+   str(x)
+ 
+   ## The logarithm of the harmonic mean
+   y0 <- log(1 / mean(1 / x))
+   print(y0)  ## -1.600885
+ 
+   lx <- log(x)
+   y1 <- log(length(x)) - logSumExp(-lx)
+   print(y1)  ## [1] -1.600885
+ 
+   # Sanity check
+   stopifnot(all.equal(y1, y0))
+ 
+   y2 <- log(length(x)) - logSumExp_R(-lx)
+   # Sanity check
+   stopifnot(all.equal(y2, y0))
+ } # for (mode ...)
mode: integer
 int [1:20] 1 1 2 2 1 2 2 2 2 1 ...
[1] 0.3215836
[1] 0.3215836
mode: double
 num [1:20] 1.53 1.74 2.15 2.82 1.4 ...
[1] 0.6673156
[1] 0.6673156
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Missing values
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> ## NA values
> lx <- c(1, 2, 3)
> lx[2] <- NA_real_
> y0 <- logSumExp_R(lx, na.rm = FALSE)
> y <- logSumExp(lx, na.rm = FALSE)
> print(y)
[1] NA
> stopifnot(identical(y, NA_real_))
> stopifnot(all.equal(y, y0))
> 
> y0 <- logSumExp_R(lx, na.rm = TRUE)
> y <- logSumExp(lx, na.rm = TRUE)
> print(y)
[1] 3.126928
> stopifnot(all.equal(y, y0))
> 
> ## NaN values
> lx <- c(1, 2, 3)
> lx[2] <- NaN
> y0 <- logSumExp_R(lx, na.rm = FALSE)
> y <- logSumExp(lx, na.rm = FALSE)
> print(y)
[1] NA
> stopifnot(identical(y, NA_real_))
> stopifnot(all.equal(y, y0))
> 
> y0 <- logSumExp_R(lx, na.rm = TRUE)
> y <- logSumExp(lx, na.rm = TRUE)
> print(y)
[1] 3.126928
> stopifnot(all.equal(y, y0))
> 
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Corner cases
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> ## Zero-length vectors
> lx <- numeric(0L)
> y0 <- logSumExp_R(lx)
> y <- logSumExp(lx)
> print(y)
[1] -Inf
> stopifnot(identical(y, -Inf))
> stopifnot(all.equal(y, y0))
> 
> ## Vectors of length one
> lx <- 1.0
> y0 <- logSumExp_R(lx)
> y <- logSumExp(lx)
> print(y)
[1] 1
> stopifnot(identical(y, lx))
> stopifnot(all.equal(y, y0))
> 
> lx <- NA_real_
> y0 <- logSumExp_R(lx, na.rm = TRUE)
> y <- logSumExp(lx, na.rm = TRUE)
> print(y)
[1] -Inf
> stopifnot(identical(y, -Inf))
> stopifnot(all.equal(y, y0))
> 
> ## All missing values
> lx <- c(NA_real_, NA_real_)
> y0 <- logSumExp_R(lx, na.rm = TRUE)
> y <- logSumExp(lx, na.rm = TRUE)
> print(y)
[1] -Inf
> stopifnot(identical(y, -Inf))
> stopifnot(all.equal(y, y0))
> 
> lx <- c(NA_real_, NA_real_)
> y0 <- logSumExp_R(lx, na.rm = FALSE)
> y <- logSumExp(lx, na.rm = FALSE)
> print(y)
[1] NA
> stopifnot(identical(y, NA_real_))
> stopifnot(all.equal(y, y0))
> 
> 
> ## +Inf values
> lx <- c(1, 2, +Inf)
> y0 <- logSumExp_R(lx)
> y <- logSumExp(lx)
> print(y)
[1] Inf
> stopifnot(identical(y, +Inf))
> stopifnot(all.equal(y, y0))
> 
> ## First element is a missing value, cf. PR #33
> lx <- c(NA_real_, 1)
> y0 <- logSumExp_R(lx)
> print(y0)
[1] NA
> y <- logSumExp(lx, na.rm = FALSE)
> print(y)
[1] NA
> stopifnot(identical(y, NA_real_))
> stopifnot(all.equal(y, y0))
> 
> y0 <- logSumExp_R(lx, na.rm = TRUE)
> print(y0)
[1] 1
> y <- logSumExp(lx, na.rm = TRUE)
> print(y)
[1] 1
> stopifnot(identical(y, 1))
> stopifnot(all.equal(y, y0))
> 
> ## Multiple -Inf values, cf. issue #84
> lx <- c(-Inf, -Inf)
> y0 <- logSumExp_R(lx)
> y <- logSumExp(lx)
> print(y)
[1] -Inf
> stopifnot(identical(y, -Inf))
> stopifnot(all.equal(y, y0))
> 
> lx <- c(-Inf, 5, -Inf)
> y0 <- logSumExp_R(lx)
> y <- logSumExp(lx)
> print(y)
[1] 5
> stopifnot(identical(y, 5))
> stopifnot(all.equal(y, y0))
> 
> proc.time()
   user  system elapsed 
   0.14    0.07    0.20