Merge pull request #433 from R-Lum/issue_424

Fix various crashes in calc_IEU() [skip ci]

NEWS.Rmd

@@ -123,6 +123,12 @@ produce unusable results (all `NA`s) or buggy behaviour (#302, fixed in #304).
 profile log-likelihood as soon as `sigma < 1e-16`, as allowing `sigma` to
 become zero leads to infinities and buggy behaviour (also fixed in #304).
 
+### `calc_IEU()`
+* The code of this function has been consolidated to avoid duplication and
+make its debugging easier: this has uncovered a small coding error and also
+led to some speed up (#429, fixed in #430, #431 and #432).
+* Some crashes in the function have been solved (#424, fixed in #433).
+
 ### `calc_Lamothe2003()`
 * We addressed a long-standing issue regarding the calculation of the `Ln/Tn`
 error after fading correction, which led to smaller than expected errors

R/calc_IEU.R

@@ -142,7 +142,6 @@ calc_IEU <- function(
     INT <- 1 / Z.bot
     R <- sqrt(INT / EXT)
     R.Error <- (2 * (1:N - 1))^(-0.5)
-
     Table.IEU <- data.table(Rank.number = 1:N,
                             De = De, De.Error = De.Total.Error,
                             Z, EXT.top, EXT, INT, R, R.Error)
@@ -153,7 +152,10 @@ calc_IEU <- function(
       do.plot(Table.IEU$Z, Table.IEU$R, Table.IEU$R.Error)
     }
 
-    Max <- Table.IEU[R >= 1, max(Rank.number, na.rm = TRUE)]
+    Max <- Table.IEU[R >= 1, suppressWarnings(max(Rank.number, na.rm = TRUE))]
+    if (is.infinite(Max)) {
+      Max <- 1
+    }
     Above <- Table.IEU[Max]
     Below <- Table.IEU[Max + 1]
     Slope <- (Above$R - Below$R) / (Above$Z - Below$Z)
@@ -179,8 +181,8 @@ calc_IEU <- function(
                       ylab = expression(paste("R = [", alpha["in"], "/", alpha["ex"],"]")),
                       abline.vals = c(1, 0),
                       asp = NA) {
-    ymin <- min((y.vals - y.errs)[-1])
-    ymax <- max((y.vals + y.errs)[-1])
+    ymin <- min((y.vals - y.errs)[-1], na.rm = TRUE)
+    ymax <- max((y.vals + y.errs)[-1], na.rm = TRUE)
 
     plot(x.vals, y.vals, type = "b", xlab = xlab, ylab = ylab,
          ylim = c(min(ymin, 0), ymax),
@@ -233,6 +235,10 @@ calc_IEU <- function(
 
   repeat {
     IEU.De <- Dbar.Mean[4]
+    if (is.na(IEU.De)) {
+      .throw_warning("Numerical error, try changing your 'a' and 'b' values")
+      break
+    }
     if (IEU.De <= Dbar) {
       break
     }

tests/testthat/test_calc_IEU.R

@@ -65,3 +65,25 @@ test_that("snapshot tests", {
       verbose = FALSE, plot = FALSE),
       tolerance = snapshot.tolerance)
 })
+
+test_that("regression tests", {
+  testthat::skip_on_cran()
+
+  ## issue 424
+  expect_warning(calc_IEU(
+      ExampleData.DeValues$CA1,
+      a = 0.45,
+      b = 1.29,
+      interval = 1,
+      verbose = FALSE,
+      plot = FALSE),
+      "Numerical error, try changing your 'a' and 'b' values")
+  expect_warning(calc_IEU(
+      ExampleData.DeValues$CA1,
+      a = 0.12,
+      b = 1.29,
+      interval = 10,
+      verbose = FALSE,
+      plot = FALSE),
+      "Numerical error, try changing your 'a' and 'b' values")
+})