Handle failure of the unfaded EXP model in calc_Huntley2006()

NEWS.Rmd

@@ -179,6 +179,10 @@ 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).
+* Fitting the "GOK" model on the unfaded data failed when the "EXP" model we
+use to find a good starting point failed: in such cases, we try again using
+the simulated fit (#549, fixed in #500; thanks to @SalOehl for reporting and
+providing data to reproduce the error).
 
 ### `calc_IEU()`
 * The code of this function has been consolidated to avoid duplication and

NEWS.md

@@ -5,7 +5,7 @@
 
 <!-- NEWS.md was auto-generated by NEWS.Rmd. Please DO NOT edit by hand!-->
 
-# Changes in version 0.9.26.9000-102 (2024-12-18)
+# Changes in version 0.9.26.9000-102 (2024-12-19)
 
 ## New functions
 
@@ -210,6 +210,10 @@
   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).
+- Fitting the “GOK” model on the unfaded data failed when the “EXP”
+  model we use to find a good starting point failed: in such cases, we
+  try again using the simulated fit (#549, fixed in \#500; thanks to
+  @SalOehl for reporting and providing data to reproduce the error).
 
 ### `calc_IEU()`
 

R/calc_Huntley2006.R

@@ -469,6 +469,7 @@ calc_Huntley2006 <- function(
 
   ## take of force_through origin settings
   force_through_origin <- GC.settings$fit.force_through_origin
+  mode_is_extrapolation <- GC.settings$mode == "extrapolation"
 
   ## call the fitting
   GC.measured <- try(do.call(fit_DoseResponseCurve, GC.settings))
@@ -651,7 +652,7 @@ calc_Huntley2006 <- function(
     .throw_warning("Ln is >10 % larger than the maximum computed LxTx value.",
                    " The De and age should be regarded as infinite estimates.")
 
-  if (Ln < min(LxTx.sim) * 0.95 && GC.settings$mode != "extrapolation")
+  if (Ln < min(LxTx.sim) * 0.95 && !mode_is_extrapolation)
     .throw_warning("Ln/Tn is smaller than the minimum computed LxTx value: ",
                    "if, in consequence, your age result is NA, either your ",
                    "input values are unsuitable, or you should consider using ",
@@ -712,9 +713,11 @@ calc_Huntley2006 <- function(
   LxTx.unfaded[is.nan((LxTx.unfaded))] <- 0
   LxTx.unfaded[is.infinite(LxTx.unfaded)] <- 0
   dosetimeGray <- dosetime * readerDdot
+
+  ## run this first model also for GOK as in general it provides more
+  ## stable estimates that can be used as starting point for GOK
   if (fit.method[1] == "EXP" || fit.method[1] == "GOK") {
-    ## we let it run regardless of the selection
-    fit_unfaded <- minpack.lm::nlsLM(
+    fit_unfaded <- try(minpack.lm::nlsLM(
       LxTx.unfaded ~ a * (1 - exp(-(dosetimeGray + c) / D0)),
       start = list(
         a = coef(fit_simulated)[["a"]],
@@ -726,15 +729,27 @@ calc_Huntley2006 <- function(
            c(Inf, max(dosetimeGray), Inf)
           },
         lower = lower.bounds[1:3],
-      control = list(maxiter = settings$maxiter)) }
+      control = list(maxiter = settings$maxiter)), silent = TRUE)
+  }
+
+  ## if this fit has failed, what we do depends on fit.method:
+  ## - for EXP, this error is irrecoverable
+  if (inherits(fit_unfaded, "try-error") && fit.method == "EXP") {
+    .throw_error("Could not fit unfaded curve, check suitability of ",
+                 "model and parameters")
+  }
+
+  ## - for GOK, we use the simulated fit to set the starting point
+  if (fit.method == "GOK") {
+    fit_start <- if (inherits(fit_unfaded, "try-error"))
+                   fit_simulated else fit_unfaded
 
-  if (fit.method[1] == "GOK") {
     fit_unfaded <- try(minpack.lm::nlsLM(
       LxTx.unfaded ~ a * (d-(1+(1/D0)*dosetimeGray*c)^(-1/c)),
       start = list(
-        a = coef(fit_unfaded)[["a"]],
-        D0 = coef(fit_unfaded)[["D0"]],
-        c = coef(fit_unfaded)[["c"]],
+        a = coef(fit_start)[["a"]],
+        D0 = coef(fit_start)[["D0"]],
+        c = coef(fit_start)[["c"]],
         d = coef(fit_simulated)[["d"]]),
       upper = if(force_through_origin) {
         c(a = Inf, D0 = max(dosetimeGray), c = Inf, d = 1)
@@ -802,7 +817,7 @@ calc_Huntley2006 <- function(
       par(oma = c(0, 3, 0, 9))
 
     # Find a good estimate of the x-axis limits
-    if(GC.settings$mode == "extrapolation" & !force_through_origin) {
+    if (mode_is_extrapolation && !force_through_origin) {
       dosetimeGray <- c(-De.measured - De.measured.error, dosetimeGray)
       De.measured <- -De.measured
     }
@@ -825,7 +840,7 @@ calc_Huntley2006 <- function(
     )
 
     ##add ablines for extrapolation
-    if(GC.settings$mode == "extrapolation")
+    if (mode_is_extrapolation)
       abline(v = 0, h = 0, col = "gray")
 
     # LxTx error bars
@@ -858,11 +873,11 @@ calc_Huntley2006 <- function(
       lty = 3)
 
     # Ln and DE as points
-    points(x = if(GC.settings$mode == "extrapolation")
+    points(x = if (mode_is_extrapolation)
                 rep(De.measured, 2)
                else
                  c(0, De.measured),
-           y = if(GC.settings$mode == "extrapolation")
+           y = if (mode_is_extrapolation)
                 c(0,0)
                else
                 c(Ln, Ln),
@@ -875,7 +890,7 @@ calc_Huntley2006 <- function(
              col = "red")
 
     # Ln as a horizontal line
-    lines(x = if(GC.settings$mode == "extrapolation")
+    lines(x = if (mode_is_extrapolation)
                 c(0, min(c(De.measured, De.sim), na.rm = TRUE))
               else
                 c(par()$usr[1], max(c(De.measured, De.sim), na.rm = TRUE)),
@@ -900,15 +915,15 @@ calc_Huntley2006 <- function(
 
     # add vertical line of simulated De
     if (!is.na(De.sim)) {
-      lines(x = if(GC.settings$mode == "extrapolation")
+      lines(x = if (mode_is_extrapolation)
                   c(-De.sim, -De.sim)
                 else
                   c(De.sim, De.sim),
             y = c(par()$usr[3], Ln),
             col = "red", lty = 3)
 
-      points(x = if(GC.settings$mode == "extrapolation") -De.sim else De.sim,
-             y = if(GC.settings$mode == "extrapolation") 0 else Ln,
+      points(x = if (mode_is_extrapolation) -De.sim else De.sim,
+             y = if (mode_is_extrapolation) 0 else Ln,
              col = "red" , pch = 16)
     } else {
       lines(x = c(De.measured, xlim[2]),

tests/testthat/test_calc_Huntley2006.R

@@ -223,6 +223,24 @@ test_that("Further tests calc_Huntley2006", {
     regexp = "\\[calc\\_Huntley2006\\(\\)\\] Ln\\/Tn is smaller than the minimum computed LxTx value.")
   })
 
+  ## failing to fit unfaded EXP model
+  ## test derived from data provided by Salome Oehler
+  input <- data.frame(dose = c(0, 1550, 9300, 37200, 71500,
+                               0, 1550, 9300, 37200, 71500),
+                      LxTx = c(0.79, 4.67, 14.41, 26.59, 24.88,
+                               0.95, 4.87, 14.17, 21.98, 25.12),
+                      LxTx_error = c(0.002, 0.006, 0.01, 0.04, 0.02,
+                                     0.002, 0.011, 0.02, 0.03, 0.03))
+  set.seed(1)
+  expect_error(calc_Huntley2006(input,
+                                rhop = c(6.5e-06, 2.0e-08),
+                                ddot = c(8.5, 1.5),
+                                fit.method = "EXP",
+                                mode = "extrapolation",
+                                readerDdot = c(0.154, 0.1),
+                                n.MC = 2),
+               "Could not fit unfaded curve, check suitability of model and")
+
   ## more coverage
   expect_s4_class(
     calc_Huntley2006(