Merge pull request #247 from R-Lum/issue_147

Fix crashes in analyse_SAR.TL()

NEWS.Rmd

@@ -24,6 +24,7 @@ it shows a warning with instructions and set `plot = FALSE`. This should prevent
 ### `analyse_SAR.TL()`
 * The function now produces a more correct `rejection.criteria` data frame
 (#245, fixed in #246).
+* Several edge cases that led to crashes have been fixed (#147, fixed in #247).
 
 ### `get_RLum()`
 * When the function was used on a list of `RLum.Analysis-class` objects with the argument `null.rm = TRUE` it would 

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.25.9000-6 (2024-09-16)
+# Changes in version 0.9.25.9000-3 (2024-09-16)
 
 ## New functions
 
@@ -28,6 +28,8 @@
 
 - The function now produces a more correct `rejection.criteria` data
   frame (#245, fixed in \#246).
+- Several edge cases that led to crashes have been fixed (#147, fixed in
+  \#247).
 
 ### `get_RLum()`
 

R/analyse_SAR.TL.R

@@ -10,7 +10,7 @@
 #'
 #' **Provided rejection criteria**
 #'
-#' `[recyling.ratio]`: calculated for every repeated regeneration dose point.
+#' `[recycling.ratio]`: calculated for every repeated regeneration dose point.
 #'
 #' `[recuperation.rate]`: recuperation rate calculated by
 #' comparing the `Lx/Tx` values of the zero regeneration point with the `Ln/Tn`
@@ -40,7 +40,7 @@
 #' specifies the general sequence structure. Three steps are allowed
 #' (`"PREHEAT"`, `"SIGNAL"`, `"BACKGROUND"`), in addition a
 #' parameter `"EXCLUDE"`. This allows excluding TL curves which are not
-#' relevant for the protocol analysis.  (**Note:** None TL are removed by default)
+#' relevant for the protocol analysis.  (**Note:** No TL are removed by default)
 #'
 #' @param rejection.criteria [list] (*with default*):
 #' list containing rejection criteria in percentage for the calculation.
@@ -193,6 +193,7 @@ analyse_SAR.TL <- function(
   ##remove TL curves which are excluded
   temp.sequence.structure <- temp.sequence.structure[which(
     temp.sequence.structure[,"protocol.step"]!="EXCLUDE"),]
+
   ##check integrity; signal and bg range should be equal
   if(length(
     unique(
@@ -217,6 +218,18 @@ analyse_SAR.TL <- function(
   TL.background.ID <- temp.sequence.structure[
     temp.sequence.structure[,"protocol.step"] == "BACKGROUND","id"]
 
+  ## check that `dose.points` is compatible with our signals:
+  ## as we expect each signal to have an Lx and a Tx components (see calls
+  ## to calc_TLLxTxRatio()), `dose.points` must divide `length(TL.signal.ID)`
+  ## in order for vector recycling to work when further down we do
+  ## `LnLxTnTx$Dose <- dose.points`
+  if (!missing(dose.points)) {
+    if ((length(TL.signal.ID) / 2) %% length(dose.points) != 0) {
+      .throw_error("Length of 'dose.points' not compatible with number ",
+                   "of signals")
+    }
+  }
+
   ##comfort ... translate integral limits from temperature to channel
   if(integral_input == "temperature"){
     signal.integral.min <-
@@ -230,147 +243,109 @@ analyse_SAR.TL <- function(
   }
 
   ##calculate LxTx values using external function
+  LnLxTnTx <- NULL
   for(i in seq(1,length(TL.signal.ID),by=2)){
-    temp.LnLxTnTx <- get_RLum(
-      calc_TLLxTxRatio(
+    Lx.data.background <- Tx.data.background <- NULL
+    if (length(TL.background.ID) > 0) {
+      Lx.data.background <- get_RLum(object, record.id = TL.background.ID[i])
+      Tx.data.background <- get_RLum(object, record.id = TL.background.ID[i + 1])
+    }
+    LxTxRatio <- calc_TLLxTxRatio(
         Lx.data.signal = get_RLum(object, record.id = TL.signal.ID[i]),
-        Lx.data.background = if (length(TL.background.ID) == 0) {
-          NULL
-        } else{
-          get_RLum(object, record.id = TL.background.ID[i])
-        },
         Tx.data.signal = get_RLum(object, record.id = TL.signal.ID[i + 1]),
-        Tx.data.background =  if (length(TL.background.ID) == 0){
-          NULL
-
-        }else{
-          get_RLum(object, record.id = TL.background.ID[i + 1])
-
-        },
+        Lx.data.background = Lx.data.background,
+        Tx.data.background = Tx.data.background,
         signal.integral.min,
         signal.integral.max
-      )
     )
+    temp.LnLxTnTx <- get_RLum(LxTxRatio)
+    rm(LxTxRatio)
 
     ##grep dose
     temp.Dose <- object@records[[TL.signal.ID[i]]]@info$IRR_TIME
-
-      ##take about NULL values
-      if(is.null(temp.Dose)){
-        temp.Dose <- NA
-
-      }
-
-    ##bind data.frame
-    temp.LnLxTnTx <- cbind(Dose=temp.Dose, temp.LnLxTnTx)
-
-    if(exists("LnLxTnTx")==FALSE){
-      LnLxTnTx <- data.frame(temp.LnLxTnTx)
-
-    }else{
-      LnLxTnTx <- rbind(LnLxTnTx,temp.LnLxTnTx)
-
+    if (is.null(temp.Dose)) {
+      temp.Dose <- NA
     }
+
+    ## append row to the data.frame
+    LnLxTnTx <- rbind(LnLxTnTx, cbind(Dose = temp.Dose, temp.LnLxTnTx))
   }
 
   ##set dose.points manually if argument was set
   if(!missing(dose.points)){
     temp.Dose <- dose.points
     LnLxTnTx$Dose <- dose.points
-
   }
 
   # Set regeneration points -------------------------------------------------
   #generate unique dose id - this are also the # for the generated points
-  temp.DoseID <- c(0:(length(LnLxTnTx[["Dose"]]) - 1))
-  temp.DoseName <- paste0("R", temp.DoseID)
-  temp.DoseName <- cbind(Name = temp.DoseName, Dose = LnLxTnTx[["Dose"]])
+  temp.DoseName <- data.frame(Name = paste0("R", seq(nrow(LnLxTnTx)) - 1),
+                              Dose = LnLxTnTx[["Dose"]])
 
   ##set natural
   temp.DoseName[temp.DoseName[, "Name"] == "R0", "Name"] <- "Natural"
 
   ##set R0
-  temp.DoseName[temp.DoseName[,"Name"]!="Natural" & temp.DoseName[,"Dose"]==0,"Name"]<-"R0"
+  temp.DoseName[temp.DoseName[, "Name"] != "Natural" &
+                temp.DoseName[, "Dose"] == 0, "Name"] <- "R0"
 
   ##find duplicated doses (including 0 dose - which means the Natural)
-  temp.DoseDuplicated<-duplicated(temp.DoseName[,"Dose"])
-
-  ##combine temp.DoseName
-  temp.DoseName<-cbind(temp.DoseName,Repeated=temp.DoseDuplicated)
+  temp.DoseName <- cbind(temp.DoseName,
+                         Repeated = duplicated(temp.DoseName[, "Dose"]))
 
   ##correct value for R0 (it is not really repeated)
   temp.DoseName[temp.DoseName[,"Dose"]==0,"Repeated"]<-FALSE
 
   ##combine in the data frame
-  temp.LnLxTnTx <- data.frame(Name = temp.DoseName[, "Name"],
-                              Repeated = as.logical(temp.DoseName[, "Repeated"]))
+  LnLxTnTx <- cbind(temp.DoseName[, c("Name", "Repeated")],
+                    LnLxTnTx)
 
-
-  LnLxTnTx<-cbind(temp.LnLxTnTx,LnLxTnTx)
-  LnLxTnTx[,"Name"]<-as.character(LnLxTnTx[,"Name"])
+  ## convert to data.table for more convenient column manipulation
+  temp <- data.table(LnLxTnTx[, c("Name", "Dose", "Repeated", "LxTx")])
 
   # Calculate Recycling Ratio -----------------------------------------------
-  RecyclingRatio <- NA_real_
-  if(length(LnLxTnTx[LnLxTnTx[,"Repeated"]==TRUE,"Repeated"])>0){
-
-    ##identify repeated doses
-    temp.Repeated<-LnLxTnTx[LnLxTnTx[,"Repeated"]==TRUE,c("Name","Dose","LxTx")]
-
-    ##find concering previous dose for the repeated dose
-    temp.Previous<-t(sapply(1:length(temp.Repeated[,1]),function(x){
-      LnLxTnTx[LnLxTnTx[,"Dose"]==temp.Repeated[x,"Dose"] &
-                 LnLxTnTx[,"Repeated"]==FALSE,c("Name","Dose","LxTx")]
-    }))
-
-    ##convert to data.frame
-    temp.Previous<-as.data.frame(temp.Previous)
-
-    ##set column names
-    temp.ColNames<-sapply(1:length(temp.Repeated[,1]),function(x){
-      paste(temp.Repeated[x,"Name"],"/",
-            temp.Previous[temp.Previous[,"Dose"]==temp.Repeated[x,"Dose"],"Name"],
-            sep="")
-    })
-
-    ##Calculate Recycling Ratio
-    RecyclingRatio<-as.numeric(temp.Repeated[,"LxTx"])/as.numeric(temp.Previous[,"LxTx"])
-
-    ##Just transform the matrix and add column names
-    RecyclingRatio<-t(RecyclingRatio)
-    colnames(RecyclingRatio)<-temp.ColNames
+  ## we first create a dummy object to use in case there are no repeated doses,
+  ## but replace it in the `if` block if there are any
+  rej.thresh <- rejection.criteria$recycling.ratio / 100
+  rej.thresh.text <- paste("\u00b1", rej.thresh) # \u00b1 is ±
+  RecyclingRatio <- data.table(criterion = "recycling ratio",
+                               value = NA,
+                               threshold = rej.thresh.text,
+                               status = NA_character_)
+  if (any(temp$Repeated)) {
+
+    ## find the index of the previous dose of each repeated dose
+    temp[, prev.idx := match(Dose, Dose)]
+
+    ## calculate the recycling ratio
+    temp[, criterion := paste0(Name, "/", Name[prev.idx])]
+    temp[, value := LxTx / LxTx[prev.idx]]
+
+    ## set status according to the given rejection threshold
+    temp[, threshold := rej.thresh.text]
+    temp[, status := fifelse(abs(1 - value) > rej.thresh, "FAILED", "OK")]
+
+    ## keep only the repeated doses
+    RecyclingRatio <- temp[Repeated == TRUE,
+                           .(criterion, value, threshold, status)]
   }
 
   # Calculate Recuperation Rate ---------------------------------------------
-  Recuperation <- NA_real_
-  if("R0" %in% LnLxTnTx[,"Name"]==TRUE){
-    Recuperation<-round(LnLxTnTx[LnLxTnTx[,"Name"]=="R0","LxTx"]/
-                          LnLxTnTx[LnLxTnTx[,"Name"]=="Natural","LxTx"],digits=4)
+  ## we first create a dummy object to use in case there is no R0 dose,
+  ## but replace it in the `if` block if there is one
+  Recuperation <- data.table(criterion = "recuperation rate",
+                             value = NA_real_,
+                             threshold = rejection.criteria$recuperation.rate / 100,
+                             status = NA_character_)
+  if ("R0" %in% temp$Name) {
+    Recuperation[, value := round(temp[Name == "R0", LxTx] /
+                                  temp[Name == "Natural", LxTx], digits = 4)]
+    Recuperation[, status := fifelse(value > threshold, "FAILED", "OK")]
   }
 
-  # Combine and Evaluate Rejection Criteria ---------------------------------
-
-  RejectionCriteria <- data.frame(
-    criterion = c(if (is.na(RecyclingRatio)) "recycling ratio" else colnames(RecyclingRatio),
-                   "recuperation rate"),
-    value = c(RecyclingRatio,Recuperation),
-    threshold = c(
-      rep(paste("\u00b1", rejection.criteria$recycling.ratio/100)
-          ,length(RecyclingRatio)),
-      rejection.criteria$recuperation.rate/100
-    ),
-    status = c(
-
-      if(is.na(RecyclingRatio)==FALSE){
-
-        sapply(1:length(RecyclingRatio), function(x){
-          if(abs(1-RecyclingRatio[x])>(rejection.criteria$recycling.ratio/100)){
-            "FAILED"
-          }else{"OK"}})}else{NA},
-
-      if(is.na(Recuperation)==FALSE) {
-        if (Recuperation > rejection.criteria$recuperation.rate / 100) "FAILED" else "OK"
-      } else NA_character_
-    ))
+  ## join the two tables and convert back to data.frame
+  RejectionCriteria <- as.data.frame(rbind(RecyclingRatio, Recuperation))
+  rm(temp)
 
   ##============================================================================##
   ##PLOTTING
@@ -526,19 +501,17 @@ analyse_SAR.TL <- function(
     lines(NTL.net.LnLx, col = col[1])
     lines(Reg1.net.LnLx, col = col[2])
 
-
     ##plot
+    TL.tmp <- TL.Plateau.LnLx[c(signal.integral.min:signal.integral.max), 2]
+    ylim.max <- quantile(TL.tmp[!is.infinite(TL.tmp)],
+                         probs = 0.90, na.rm = TRUE)
     par(new = TRUE)
     plot(
       TL.Plateau.LnLx,
       axes = FALSE,
       xlab = "",
       ylab = "",
-      ylim = c(0,
-               quantile(
-                 TL.Plateau.LnLx[c(signal.integral.min:signal.integral.max), 2],
-                 probs = c(0.90), na.rm = TRUE
-               ) + 3),
+      ylim = c(0, ylim.max + 3),
       col = "darkgreen"
     )
     axis(4)
@@ -651,15 +624,20 @@ analyse_SAR.TL <- function(
     ...
   ))
 
-  ##check for error
+  ## plot_GrowthCurve() can fail in two ways:
+  ## 1. either with a hard error, in which case there's nothing much we
+  ##    can do and stop early by returning NULL
   if(inherits(temp.GC, "try-error")){
     return(NULL)
-
-  }else{
-    temp.GC <- get_RLum(temp.GC)[, c("De", "De.Error")]
-
   }
 
+  ## 2. or with a soft error by returning NULL, in which case we set
+  ##    temp.GC to NA and continue (this can be done after the call to
+  ##    get_RLum(), as it deals well with NULLs)
+  temp.GC <- get_RLum(temp.GC)[, c("De", "De.Error")]
+  if (is.null(temp.GC))
+    temp.GC <- NA
+
   ##add rejection status
   if(length(grep("FAILED",RejectionCriteria$status))>0){
     temp.GC <- data.frame(temp.GC, RC.Status="FAILED")