Merge pull request #542 from R-Lum/issue_258_p2

Speed up calc_Huntley2006(): set a more consistent default for rprime

NEWS.Rmd

@@ -158,7 +158,13 @@ we now take the absolute values of errors (#479, fixed in #481).
 * Some crashes in case of model misspecification have been solved (#538,
 fixed in #539).
 * Some details in the implementation of the function have been optimized, and
-now it is much faster than before (#258, fixed in #541).
+now it is much faster than before. As part of this, we have changed the default
+setting of the `rprime` vector that is used in the calculation of the natural
+dose response and the field saturation, so that more points are concentrated
+in the bulk of the distribution: this previously depended incorrectly on the
+number of Monte Carlo iterations requested, so this change brings an additional
+speed boost. The default setting can be overridden via the `rprime` argument
+(#258, fixed in #541 and #542).
 
 ### `calc_IEU()`
 * The code of this function has been consolidated to avoid duplication and

NEWS.md

@@ -186,8 +186,14 @@
 - Some crashes in case of model misspecification have been solved (#538,
   fixed in \#539).
 - Some details in the implementation of the function have been
-  optimized, and now it is much faster than before (#258, fixed in
-  \#541).
+  optimized, and now it is much faster than before. As part of this, we
+  have changed the default setting of the `rprime` vector that is used
+  in the calculation of the natural dose response and the field
+  saturation, so that more points are concentrated in the bulk of the
+  distribution: this previously depended incorrectly on the number of
+  Monte Carlo iterations requested, so this change brings an additional
+  speed boost. The default setting can be overridden via the `rprime`
+  argument (#258, fixed in \#541 and \#542).
 
 ### `calc_IEU()`
 

R/calc_Huntley2006.R

@@ -61,6 +61,10 @@
 #' the sample under investigation using the sample specific \eqn{\rho}',
 #' unfaded \eqn{D_0} and \eqn{\dot{D}} values, following the approach of Kars et al. (2008).
 #'
+#' The computation is done using 1000 equally-spaced points in the interval
+#' \[0.01, 3\]. This can be controlled by setting option `rprime`, such as
+#' in `rprime = seq(0.01, 3, length.out = 1000)` (the default).
+#'
 #' **Uncertainties**
 #'
 #' Uncertainties are reported at \eqn{1\sigma} and are assumed to be normally
@@ -439,6 +443,15 @@ calc_Huntley2006 <- function(
   Ln <- data[["LxTx"]][1]
   Ln.error <- data[["LxTx.Error"]][1]
 
+  ## set a sensible default for rprime: in most papers the upper boundary is
+  ## around 2.2, so setting it to 3 should be enough in general, and 1000
+  ## points seem also enough; in any case, we let the user override it
+  rprime <- seq(0.01, 3, length.out = 1000)
+  if ("rprime" %in% names(list(...))) {
+    rprime <- list(...)$rprime
+    .validate_class(rprime, "numeric")
+  }
+
   ## (1) MEASURED ----------------------------------------------------
 
   data.tmp <- data
@@ -561,7 +574,6 @@ calc_Huntley2006 <- function(
   ddots <- ddot / ka
   natdosetimeGray <- c(0, exp(seq(1, log(max(data[ ,1]) * 2), length.out = 999)))
   natdosetime <- natdosetimeGray
-  rprime <- seq(0.01, 5, length.out = 500)
   pr <- 3 * rprime^2 * exp(-rprime^3) # Huntley 2006, eq. 3
   K <- Hs * exp(-rhop[1]^-(1/3) * rprime)
   TermA <- matrix(NA, nrow = length(rprime), ncol = length(natdosetime))
@@ -675,8 +687,6 @@ calc_Huntley2006 <- function(
   ##  rho_i <- rhop_MC[i]^(1 / 3)
   ##  tau <- ((1 / Hs) * exp(1)^(rprime / rho_i)) / ka
 
-  rprime <- seq(0.01, 5, length.out = settings$n.MC)
-  pr <- 3 * rprime^2 * exp(-rprime^3)
   rho_MC <- rhop_MC^(1 / 3)
   nN_SS_MC <- mapply(function(rho_i, ddot_i, UFD0_i) {
     tau <- ((1 / Hs) * exp(1)^(rprime / rho_i)) / ka

tests/testthat/test_calc_Huntley2006.R

@@ -73,6 +73,10 @@ test_that("input validation", {
   expect_error(calc_Huntley2006(data, rhop = rhop,
                                 ddot = ddot, readerDdot = 0.13),
                "'readerDdot' should have length 2")
+  expect_error(calc_Huntley2006(data, rhop = rhop,
+                                ddot = ddot, readerDdot = readerDdot,
+                                rprime = list()),
+               "'rprime' should be of class 'numeric'")
 
   SW({
   expect_warning(calc_Huntley2006(data[, 1:2], rhop = rhop, n.MC = 2,
@@ -160,7 +164,7 @@ test_that("Further tests calc_Huntley2006", {
         mode = "extrapolation",
         plot = TRUE,
         verbose = FALSE),
-      tolerance = snapshot.tolerance)
+      tolerance = max(snapshot.tolerance, 8.0e-3))
 
   ## EXP ... normal
   set.seed(1)
@@ -232,6 +236,7 @@ test_that("Further tests calc_Huntley2006", {
       data = data[1:10, ],
       LnTn = data[1, c(2, 3)],
       rhop = rhop, ddot = ddot, readerDdot = readerDdot,
+      rprime = c(0.01, 2.2, length.out = 500),
       n.MC = 2, plot = FALSE, verbose = FALSE),
     class = "RLum.Results")