calc_gSGC.R

#' @title Calculate De value based on the gSGC by Li et al., 2015
#'
#' @description The function computes De value and De value error using the
#' global standardised growth curve (gSGC) assumption proposed by Li et al.,
#' 2015 for OSL dating of sedimentary quartz.
#'
#' @details
#' The error of the De value is determined using a Monte Carlo simulation approach.
#' Solving of the equation is realised using [uniroot].
#' Large values for `n.MC` will significantly increase the computation time.
#'
#'
#' @param data [data.frame] (**required**):
#' input data  the following columns five columns in the given order:
#' `LnTn`, `LnTn.error`, `Lr1Tr1`, `Lr1Tr1.error`, `Dr1`. Column names are
#' not required.
#'
#' @param gSGC.type [character] (*with default*):
#' function parameters to use for the iteration procedure, either `"0-450"` or
#' `"0-250"`, as presented in Li et al., 2015 (Table 2). This is ignored if
#' `gSGC.parameters` is set.
#'
#' @param gSGC.parameters [list] (*optional*):
#' option to provide own function parameters used for fitting as named list.
#' Nomenclature follows Li et al., 2015, i.e. `list(A, A.error, D0, D0.error,
#' c, c.error, Y0, Y0.error, range)`, where `range` is defines the interval
#' where the function is considered as valid, e.g. `range = c(0,250)`.\cr
#' If set, option `gSGC.type` will be ignored.
#'
#' @param n.MC [integer] (*with default*):
#' number of Monte Carlo simulation runs for error estimation, see details.
#'
#' @param verbose [logical] (*with default*):
#' enable/disable output to the terminal.
#'
#' @param plot [logical] (*with default*):
#' enable/disable the plot output.
#'
#' @param ... parameters will be passed to the plot output
#'
#' @return Returns an S4 object of type [RLum.Results-class].
#'
#' **`@data`**\cr
#' `$ De.value` ([data.frame]) \cr
#'  `.. $ De`  \cr
#'  `.. $ De.error` \cr
#'  `.. $ Eta` \cr
#' `$ De.MC` ([list]) contains the matrices from the error estimation.\cr
#' `$ uniroot` ([list]) contains the [uniroot] outputs of the De estimations\cr
#'
#' **`@info`**\cr
#' `$ call`` ([call]) the original function call
#'
#' @section Function version: 0.1.2
#'
#' @author
#' Sebastian Kreutzer, Institute of Geography, Heidelberg University (Germany)
#'
#' @seealso [RLum.Results-class], [get_RLum], [uniroot]
#'
#' @references
#' Li, B., Roberts, R.G., Jacobs, Z., Li, S.-H., 2015. Potential of establishing
#' a 'global standardised growth curve' (gSGC) for optical dating of quartz from sediments.
#' Quaternary Geochronology 27, 94-104. doi:10.1016/j.quageo.2015.02.011
#'
#' @keywords datagen
#'
#' @examples
#'
#' results <- calc_gSGC(data = data.frame(
#' LnTn =  2.361, LnTn.error = 0.087,
#' Lr1Tr1 = 2.744, Lr1Tr1.error = 0.091,
#' Dr1 = 34.4))
#'
#' get_RLum(results, data.object = "De")
#'
#' @md
#' @export
calc_gSGC<- function(
  data,
  gSGC.type = "0-250",
  gSGC.parameters,
  n.MC = 100,
  verbose = TRUE,
  plot = TRUE,
  ...
) {
  .set_function_name("calc_gSGC")
  on.exit(.unset_function_name(), add = TRUE)

  ## Integrity checks -------------------------------------------------------

  .validate_class(data, "data.frame")
  if (ncol(data) != 5)
    .throw_error("'data' is expected to have 5 columns")
  gSGC.type <- .validate_args(gSGC.type, c("0-250", "0-450"))

  ##rename columns for consistency reasons
  colnames(data) <- c('LnTn', 'LnTn.error', 'Lr1Tr1', 'Lr1Tr1.error', 'Dr1')

  ## ensure errors are not negative
  data$LnTn.error <- abs(data$LnTn.error)
  data$Lr1Tr1.error <- abs(data$Lr1Tr1.error)


##============================================================================##
##DEFINE FUNCTION
##============================================================================##

    ##define function, nomenclature according to publication that should be solved
    f <- function(x,A,D0,c,Y0,Dr1,Lr1Tr1,LnTn) {
      (((A * (1 - exp( - Dr1 / D0))) + c * Dr1 + Y0)/Lr1Tr1) -
      (((A * (1 - exp( - x/D0))) + c * x + Y0)/LnTn)
    }

    ##set general parameters
    if (!missing(gSGC.parameters)) {
      A <- gSGC.parameters$A
      A.error <- gSGC.parameters$A.error
      D0 <- gSGC.parameters$D0
      D0.error <- gSGC.parameters$D0.error
      c <- gSGC.parameters$c
      c.error <- gSGC.parameters$c.error
      Y0 <- gSGC.parameters$Y0
      Y0.error <- gSGC.parameters$Y0.error
      range <- gSGC.parameters$range

    }else{
      if (gSGC.type == "0-450") {
        A <- 0.723
        A.error <- 0.014
        D0 <- 65.1
        D0.error <- 0.9
        c <- 0.001784
        c.error <- 0.000016
        Y0 <- 0.009159
        Y0.error <- 0.004795

        range <- c(0.1,250)

      }else if (gSGC.type == "0-250") {
        A <- 0.787
        A.error <- 0.051
        D0 <- 73.9
        D0.error <- 2.2
        c <- 0.001539
        c.error <- 0.000068
        Y0 <- 0.01791
        Y0.error <- 0.00490

        range <- c(0.1,250)
      }
    }

    ##Define size of output objects
    output.data <- data.table::data.table(
      DE = numeric(length = nrow(data)),
      DE.ERROR =  numeric(length = nrow(data)),
      ETA =  numeric(length = nrow(data))
    )

    ##set list for De.MC
    output.De.MC <- vector("list",  nrow(data))

    ##set list for uniroot
    output.uniroot <-  vector("list",  nrow(data))


##============================================================================##
##CALCULATION
##============================================================================##
for(i in 1:nrow(data)){
    Lr1Tr1 <- data[i, "Lr1Tr1"]
    Lr1Tr1.error <- data[i,"Lr1Tr1.error"]
    Dr1 <- data[i,"Dr1"]
    Dr1.error <- data[i,"Dr1.error"]

    LnTn <- data[i,"LnTn"]
    LnTn.error <- data[i,"LnTn.error"]

  ##calculate mean value
  temp <- try(uniroot(
      f,
      interval = c(0.1,450),
      tol = 0.001,
      A = A,
      D0 = D0,
      c = c,
      Y0 = Y0,
      Dr1 = Dr1,
      Lr1Tr1 = Lr1Tr1,
      LnTn = LnTn,
      extendInt = 'yes',
      check.conv = TRUE,
      maxiter = 1000
    ), silent = TRUE)

  ## allocate matrix
  temp.MC.matrix <- matrix(nrow = n.MC, ncol = 8)
  colnames(temp.MC.matrix) <- c("LnTn", "Lr1Tr1", "A", "D0", "c", "Y0", "De", "Eta")

  ## fill the first 6 columns of the matrix
  temp.MC.matrix[, 1:6] <- matrix(rnorm(
      n.MC * 6,
      mean = c(LnTn, Lr1Tr1, A, D0, c, Y0),
      sd = c(LnTn.error, Lr1Tr1.error, A.error, D0.error, c.error, Y0.error)
  ), ncol = 6, byrow = TRUE)

  if(!inherits(temp, "try-error")){

    ##get De
    De <- temp$root

    ##calculate Eta, which is the normalisation factor
    Eta <- ((A * (1 - exp( - Dr1 / D0))) + c * Dr1 + Y0)/Lr1Tr1

    ##--------------------------------------------------------------------------##
    ##Monte Carlo simulation for error estimation

      ##run uniroot to get the De
      temp.MC.matrix[,7] <- vapply(X = 1:n.MC, FUN = function(x){

        uniroot(f,
                interval = c(0.1,450),
                tol = 0.001,
                A = temp.MC.matrix[x,3],
                D0 = temp.MC.matrix[x,4],
                c = temp.MC.matrix[x,5],
                Y0 = temp.MC.matrix[x,6],
                Dr1 = Dr1,
                Lr1Tr1 =temp.MC.matrix[x,2],
                LnTn = temp.MC.matrix[x,1],
                check.conv = TRUE,
                extendInt = 'yes',
                maxiter = 1000
                )$root

      }, FUN.VALUE = vector(mode = "numeric", length = 1))

      ##calculate also the normalisation factor
      temp.MC.matrix[,8] <- (temp.MC.matrix[,3] * (1 - exp( - Dr1 / temp.MC.matrix[,4])) +
        temp.MC.matrix[,5] * Dr1 + temp.MC.matrix[,6])/temp.MC.matrix[,2]

      ##get De error as SD
      De.error <- sd(temp.MC.matrix[,7])

  }else{
    .throw_warning("No solution was found")
    De <- NA
    Eta <- NA
    De.error <- NA
  }


# Plot output -------------------------------------------------------------
  if (plot) {

    ##set plot settings
    plot.settings <- list(
      main = "gSGC and resulting De",
      xlab = "Dose [a.u.]",
      ylab = expression(paste("Re-norm. ", L[x]/T[x])),
      xlim = NULL,
      ylim = NULL,
      lwd = 1,
      lty = 1,
      pch = 21,
      col = "red",
      grid = expression(nx = 10, ny = 10),
      mtext = ""
    )

    plot.settings <-  modifyList(plot.settings, list(...))

    ##graphical feedback
    x <- NA
    curve(
      A * (1 - exp(-x / D0)) + c * x + Y0, from = 0, to = 500,
      xlab = plot.settings$xlab,
      ylab = plot.settings$ylab,
      main = plot.settings$main,
      xlim = plot.settings$xlim,
      ylim = plot.settings$ylim,
      lwd = plot.settings$lwd,
      lty = plot.settings$lty
    )

    mtext(side = 3, plot.settings$mtext)

    if(!is.null(plot.settings$grid)){
      graphics::grid(eval(plot.settings$grid))
    }

    if(!inherits(temp, "try-error")){

      if(temp$root < 450 & temp$root > 0){
        points(temp$root,Eta*LnTn, col = plot.settings$col, pch = plot.settings$pch)
        segments(De - De.error,Eta * LnTn,
                 De + De.error,Eta * LnTn)

        hist <-
          hist(
            temp.MC.matrix[, 7],
            freq = FALSE,
            add = TRUE,
            col = rgb(0, 0, 0, 0.2),
            border = rgb(0, 0, 0, 0.5)
          )
        lines(hist$mids,hist$density)

      }else{

        if(temp$root < 450){
          x0 <- 450
          x1 <- 500
        } else {
          x0 <- 50
          x1 <- 0
        }
        shape::Arrows(
            x0 = x0,
            y0 = par()$usr[4] - 0.2,
            x1 = x1,
            y1 = par()$usr[4] - 0.2,
            arr.type = "triangle",
            col = "red"
        )

        mtext(side = 1, text = "Out of bounds!", col = "red")
      }

    }else{
      mtext(side = 1, text = "No solution found!", col = "red")
    }
  }


# Terminal output ---------------------------------------------------------
    if (verbose) {
      cat("\n[calc_gSGC()]")
      cat("\n Corresponding De based on the gSGC\n")

      cat(paste0("\n"," Ln/Tn:\t\t ",LnTn," \u00B1 ", LnTn.error,"\n"))
      cat(paste0(""," Lr1/Tr1:\t ",Lr1Tr1," \u00B1 ", Lr1Tr1.error,"\n"))
      cat(paste0(""," Dr1:\t\t ",Dr1,"\n"))
      cat(paste0(""," f(D):\t\t ",A," * (1 - exp(-D /",D0,")) + c * D + ",Y0,"\n"))
      cat(paste0(""," n.MC:\t\t ",n.MC,"\n"))
      cat(paste0(" ------------------------------ \n"))
      cat(paste0(" De:\t\t",round(De,digits = 2)," \u00B1 ",round(De.error,digits = 2),"\n"))
      cat(paste0(" ------------------------------ \n"))
    }


##============================================================================##
##CREATE OUTPUT OBJECTS
##============================================================================##

    ##replace values in the data.table with values
    output.data[i, `:=` (DE = De,
                         DE.ERROR = De.error,
                         ETA = Eta)]

    ##matrix - to prevent memory overload limit output
    if(n.MC * nrow(data) > 1e6){
      # nocov start
      if(i == 1){
        output.De.MC[[i]] <- temp.MC.matrix
      }else{
        output.De.MC[[i]] <- NA
      }

      .throw_warning("Only the first MC matrix is returned to prevent ",
                     "memory overload")
      # nocov end
    }else{
      output.De.MC[[i]] <- temp.MC.matrix
    }

    output.uniroot[[i]] <- temp

}##end for loop

##============================================================================##
##OUTPUT RLUM
##============================================================================##

    temp.RLum.Results <- set_RLum(
      class = "RLum.Results",
      data = list(
        De = as.data.frame(output.data),
        De.MC =  output.De.MC,
        uniroot = output.uniroot
      ),
      info = list( call = sys.call())
    )

  return(temp.RLum.Results)
}