#' Verify single grain data sets and check for invalid grains, i.e.
#' zero-light level grains
#'
#' This function tries to identify automatically zero-light level curves (grains)
#' from single grain data measurements.
#'
#' **How does the method work?**
#'
#' The function compares the expected values (\eqn{E(X)}) and the variance
#' (\eqn{Var(X)}) of the count values for each curve. Assuming that the
#' background roughly follows a Poisson distribution, the absolute difference
#' of both values should be zero or at least around zero as
#'
#' \deqn{E(x) = Var(x) = \lambda}
#'
#' Thus the function checks for:
#'
#' \deqn{abs(E(x) - Var(x)) >= \Theta}
#'
#' With \eqn{\Theta} an arbitrary, user defined, threshold. Values above the
#' threshold indicate curves comprising a signal.
#'
#' Note: the absolute difference of \eqn{E(X)} and \eqn{Var(x)} instead of the
#' ratio was chosen as both terms can become 0 which would result in 0 or `Inf`,
#' if the ratio is calculated.
#'
#' @param object [Risoe.BINfileData-class] or [RLum.Analysis-class] (**required**):
#' input object. The function also accepts a list with objects of allowed type.
#'
#' @param threshold [numeric] (*with default*):
#' numeric threshold value for the allowed difference between the `mean` and
#' the `var` of the count values (see details)
#'
#' @param cleanup [logical] (*with default*):
#' if set to `TRUE`, curves identified as zero light level curves are
#' automatically removed. Output is an object as same type as the input, i.e.
#' either [Risoe.BINfileData-class] or [RLum.Analysis-class]
#'
#' @param cleanup_level [character] (*with default*):
#' selects the level for the clean-up of the input data sets.
#' Two options are allowed: `"curve"` or `"aliquot"`:
#'
#' - If `"curve"` is selected, every single curve marked as `invalid` is removed.
#' - If `"aliquot"` is selected, curves of one aliquot (grain or disc) can be
#' marked as invalid, but will not be removed. An aliquot will be only removed
#' if all curves of this aliquot are marked as invalid.
#'
#' @param verbose [logical] (*with default*):
#' enables or disables the terminal feedback
#'
#' @param plot [logical] (*with default*):
#' enables or disables the graphical feedback
#'
#' @param ... further parameters to control the plot output; if selected.
#' Supported arguments `main`, `ylim`
#'
#' @return
#' The function returns
#'
#' -----------------------------------\cr
#' `[ NUMERICAL OUTPUT ]`\cr
#' -----------------------------------\cr
#'
#' **`RLum.Results`**-object
#'
#' **slot:****`@data`**
#'
#' \tabular{lll}{
#' **Element** \tab **Type** \tab **Description**\cr
#' `$unique_pairs` \tab `data.frame` \tab the unique position and grain pairs \cr
#' `$selection_id` \tab `numeric` \tab the selection as record ID \cr
#' `$selection_full` \tab `data.frame` \tab implemented models used in the baSAR-model core \cr
#' }
#'
#' **slot:****`@info`**
#'
#' The original function call
#'
#' **Output variation**
#'
#' For `cleanup = TRUE` the same object as the input is returned, but cleaned up
#' (invalid curves were removed). This means: Either a [Risoe.BINfileData-class]
#' or an [RLum.Analysis-class] object is returned in such cases.
#' A [Risoe.BINfileData-class] object can be exported to a BIN-file by
#' using the function [write_R2BIN].
#'
#' @note
#' This function can work with [Risoe.BINfileData-class] objects or
#' [RLum.Analysis-class] objects (or a list of it). However, the function is
#' highly optimised for [Risoe.BINfileData-class] objects as it make sense to
#' remove identify invalid grains before the conversion to an
#' [RLum.Analysis-class] object.
#'
#' The function checking for invalid curves works rather robust and it is likely
#' that Reg0 curves within a SAR cycle are removed as well. Therefore it is
#' strongly recommended to use the argument `cleanup = TRUE` carefully.
#'
#' @section Function version: 0.2.3
#'
#'
#' @author
#' Sebastian Kreutzer, Institute of Geography, Heidelberg University (Germany)
#'
#'
#' @seealso [Risoe.BINfileData-class], [RLum.Analysis-class], [write_R2BIN],
#' [read_BIN2R]
#'
#' @keywords manip datagen
#'
#' @examples
#'
#' ##01 - basic example I
#' ##just show how to apply the function
#' data(ExampleData.XSYG, envir = environment())
#'
#' ##verify and get data.frame out of it
#' verify_SingleGrainData(OSL.SARMeasurement$Sequence.Object)$selection_full
#'
#' ##02 - basic example II
#' data(ExampleData.BINfileData, envir = environment())
#' id <- verify_SingleGrainData(object = CWOSL.SAR.Data,
#' cleanup_level = "aliquot")$selection_id
#'
#' \dontrun{
#' ##03 - advanced example I
#' ##importing and exporting a BIN-file
#'
#' ##select and import file
#' file <- file.choose()
#' object <- read_BIN2R(file)
#'
#' ##remove invalid aliquots(!)
#' object <- verify_SingleGrainData(object, cleanup = TRUE)
#'
#' ##export to new BIN-file
#' write_R2BIN(object, paste0(dirname(file),"/", basename(file), "_CLEANED.BIN"))
#' }
#'
#' @md
#' @export
verify_SingleGrainData <- function(
object,
threshold = 10,
cleanup = FALSE,
cleanup_level = 'aliquot',
verbose = TRUE,
plot = FALSE,
...
){
##three types of input are allowed:
##(1) RisoeBINfileData
##(2) RLum.Analysis
##(3) List of RLum.Analysis
# Self Call -----------------------------------------------------------------------------------
if(is(object, "list")){
results <- .warningCatcher(lapply(1:length(object), function(x) {
verify_SingleGrainData(
object = object[[x]],
threshold = threshold,
cleanup = cleanup,
cleanup_level = cleanup_level,
verbose = verbose,
plot = plot,
main = paste0("Record #",x)
)
}))
##account for cleanup
if(cleanup){
return(results)
}else{
return(merge_RLum(results))
}
}
##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
##RisoeBINfileData
if(is(object, "Risoe.BINfileData")){
##run test on DATA slot
##MEAN + SD
temp.results_matrix <- lapply(X = object@DATA, FUN = function(x){
c(mean(x), var(x))
})
temp.results_matrix <- do.call(rbind, temp.results_matrix)
##DIFF
temp.results_matrix_RATIO <- temp.results_matrix[,2]/temp.results_matrix[,1]
##SEL
temp.results_matrix_VALID <- temp.results_matrix_RATIO > threshold
##combine everything to in a data.frame
selection <- data.frame(
POSITION = object@METADATA$POSITION,
GRAIN = object@METADATA$GRAIN,
MEAN = temp.results_matrix[, 1],
VAR = temp.results_matrix[, 2],
RATIO = temp.results_matrix_RATIO,
THRESHOLD = rep_len(threshold, length(object@DATA)),
VALID = temp.results_matrix_VALID
)
##get unique pairs for POSITION and GRAIN for VALID == TRUE
unique_pairs <- unique(
selection[selection[["VALID"]], c("POSITION", "GRAIN")])
if(cleanup_level == "aliquot"){
selection_id <- sort(unlist(lapply(1:nrow(unique_pairs), function(x) {
which(
.subset2(selection, 1) == .subset2(unique_pairs, 1)[x] &
.subset2(selection, 2) == .subset2(unique_pairs, 2)[x]
)
})))
}else{
##reduce data to TRUE selection
selection_id <- which(selection[["VALID"]])
}
##select output on the chosen input
if(cleanup){
##selected wanted elements
object@DATA <- object@DATA[selection_id]
object@METADATA <- object@METADATA[selection_id,]
object@METADATA$ID <- 1:length(object@DATA)
##print message
selection_id <- paste(selection_id, collapse = ", ")
if(verbose){
cat(paste0("\n[verify_SingleGrainData()] Risoe.BINfileData object reduced to records: \n", selection_id))
cat("\n\n[verify_SingleGrainData()] Risoe.BINfileData object record index reset.\n")
}
##return
return_object <- object
}else{
return_object <- set_RLum(
class = "RLum.Results",
data = list(
unique_pairs = unique_pairs,
selection_id = selection_id,
selection_full = selection),
info = list(call = sys.call())
)
}
##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
##RLum.Analysis and list with RLum.Analysis objects
## ... and yes it make sense not to mix that up with the code above
}else if(is(object,"RLum.Analysis")){
##first extract all count values from all curves
object_list <- lapply(get_RLum(object), function(x){
##yes, would work differently, but it is faster
x@data[,2]
})
##MEAN + SD
temp.results_matrix <- lapply(X = object_list, FUN = function(x){
c(mean(x), var(x))
})
temp.results_matrix <- do.call(rbind, temp.results_matrix)
##DIFF
temp.results_matrix_RATIO <- temp.results_matrix[,2]/temp.results_matrix[,1]
##SEL
temp.results_matrix_VALID <- temp.results_matrix_RATIO > threshold
##get structure for the RLum.Analysis object
temp_structure <- structure_RLum(object, fullExtent = TRUE)
##now we have two cases, depending on where measurement is coming from
if (object@originator == "Risoe.BINfileData2RLum.Analysis") {
##combine everything to in a data.frame
selection <- data.frame(
POSITION = temp_structure$info.POSITION,
GRAIN = temp_structure$info.GRAIN,
MEAN = temp.results_matrix[, 1],
VAR = temp.results_matrix[, 2],
RATIO = temp.results_matrix_RATIO,
THRESHOLD = rep_len(threshold, length(object_list)),
VALID = temp.results_matrix_VALID
)
##get unique pairs for POSITION and GRAIN for VALID == TRUE
unique_pairs <- unique(
selection[selection[["VALID"]], c("POSITION", "GRAIN")])
} else if (object@originator == "read_XSYG2R") {
##combine everything to in a data.frame
selection <- data.frame(
POSITION = if(any(grepl(pattern = "position", names(temp_structure)))){
temp_structure$info.position}else{
NA
},
GRAIN = NA,
MEAN = temp.results_matrix[, 1],
VAR = temp.results_matrix[, 2],
RATIO = temp.results_matrix_RATIO,
THRESHOLD = rep_len(threshold, length(object_list)),
VALID = temp.results_matrix_VALID
)
##get unique pairs for POSITION for VALID == TRUE
unique_pairs <- unique(
selection[["POSITION"]][selection[["VALID"]]])
} else{
stop("[verify_SingleGrainData()] I don't know what to do object 'originator' not supported!",
call. = FALSE)
}
##set up cleanup
if(cleanup_level == "aliquot") {
if (object@originator == "read_XSYG2R") {
if(!is.na(unique_pairs)){
selection_id <-
sort(unlist(lapply(1:nrow(unique_pairs), function(x) {
which(.subset2(selection, 1) == .subset2(unique_pairs, 1)[x])
})))
}else{
selection_id <- NA
}
} else if (object@originator == "Risoe.BINfileData2RLum.Analysis") {
selection_id <-
sort(unlist(lapply(1:nrow(unique_pairs), function(x) {
which(
.subset2(selection, 1) == .subset2(unique_pairs, 1)[x] &
.subset2(selection, 2) == .subset2(unique_pairs, 2)[x]
)
})))
}
##make sure that we do not break subsequent code
if(length(selection_id) == 0) selection_id <- NA
} else{
##reduce data to TRUE selection
selection_id <- which(selection[["VALID"]])
}
##return value
##select output on the chosen input
if(cleanup && !any(is.na(selection_id))){
##print message
if(verbose){
selection_id_text <- paste(selection_id, collapse = ", ")
cat(paste0("\n[verify_SingleGrainData()] RLum.Analysis object reduced to records: ",
selection_id_text), "\n")
}
##selected wanted elements
if (length(selection_id) == 0) {
object <- set_RLum(
class = "RLum.Analysis",
originator = object@originator,
protocol = object@protocol,
records = list(),
info = list(
unique_pairs = unique_pairs,
selection_id = selection_id,
selection_full = selection)
)
} else{
object <- set_RLum(
class = "RLum.Analysis",
records = get_RLum(object, record.id = selection_id, drop = FALSE),
info = list(
unique_pairs = unique_pairs,
selection_id = selection_id,
selection_full = selection)
)
}
##return
return_object <- object
}else{
if(any(is.na(selection_id))){
warning("[verify_SingleGrainData()] selection_id is NA, nothing removed, everything selected for removal!",
call. = FALSE)
}
return_object <- set_RLum(
class = "RLum.Results",
data = list(
unique_pairs = unique_pairs,
selection_id = selection_id,
selection_full = selection),
info = list(call = sys.call())
)
}
}else{
stop("[verify_SingleGrainData()] Input type '", is(object)[1],
"' is not allowed for this function!", call. = FALSE)
}
# Plot ----------------------------------------------------------------------------------------
if(plot){
##set plot settings
plot_settings <-
modifyList(x = list(
main = "Record selection",
ylim = range(c(selection[["RATIO"]], threshold * 1.1))
),
val = list(...))
##plot area
plot(
NA,
NA,
xlim = c(1,nrow(selection)),
ylim = plot_settings$ylim,
log = "y",
xlab = "Record index",
ylab = "Calculated ratio [a.u.]",
main = plot_settings$main
)
##plot points above the threshold
points(x = which(selection[["VALID"]]),
y = selection[["RATIO"]][selection[["VALID"]]], pch = 20, col = "darkgreen")
points(x = which(!selection[["VALID"]]),
y = selection[["RATIO"]][!selection[["VALID"]]], pch = 20, col = rgb(0,0,0,0.5))
abline(h = threshold, col = "red", lty = 1, lwd = 2)
mtext(
side = 3,
text = paste0(
"(total: ", nrow(selection),
" | valid: ", length(which(selection[["VALID"]])),
" | invalid: ", length(which(!selection[["VALID"]])), ")"),
cex = 0.9 * par()$cex)
}
# Return --------------------------------------------------------------------------------------
return(return_object)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.