R/plot_DRCSummary.R

In R-Lum/Luminescence: Comprehensive Luminescence Dating Data Analysis

#' @title Create a Dose-Response Curve Summary Plot
#'
#' @description
#' While analysing OSL SAR or pIRIR-data the view on the data is usually
#' limited to one dose-response curve (DRC) at the time for one aliquot. This
#' function overcomes this limitation by plotting all DRCs from an
#' [RLum.Results-class] object created by [analyse_SAR.CWOSL] in one single
#' plot.
#'
#' If you want plot your DRC on an energy scale (dose in Gy), you can either
#' use option `source_dose_rate` or perform your SAR analysis with the dose
#' points in Gy (better axis scaling).
#'
#' @param object [RLum.Results-class] (**required**): input object created by
#' [analyse_SAR.CWOSL]. The input object can be provided as [list].
#'
#' @param source_dose_rate [numeric] (*optional*): allows to modify the axis
#' and show values in Gy, instead seconds. Only a single numerical value is
#' allowed.
#'
#' @param sel_curves [numeric] (*optional*): id of the curves to be plotted in
#' its occurring order. A sequence can be provided for selecting, e.g., only
#' every 2nd curve from the input object.
#'
#' @param show_dose_points [logical] (*with default*): enable/disable plotting
#' of dose points in the graph.
#'
#' @param show_natural [logical] (*with default*): enable/disable the plot
#' of the natural `Lx/Tx` values.
#'
#' @param n [integer] (*with default*): number of x-values used to evaluate
#' one curve object. Large numbers slow down the plotting process and are
#' usually not needed.
#'
#'@param ... Further arguments and graphical parameters to be passed.
#'In particular: `main`, `xlab`, `ylab`, `xlim`, `ylim`, `lty`, `lwd`, `pch`, `col.pch`, `col.lty`, `mtext`
#'
#'@section Function version: 0.2.4
#'
#'@return An [RLum.Results-class] object is returned:
#'
#' Slot: **@data**\cr
#'
#' \tabular{lll}{
#' **OBJECT** \tab **TYPE** \tab **COMMENT**\cr
#' `results` \tab [data.frame] \tab with dose and LxTx values \cr
#' `data` \tab [RLum.Results-class] \tab original input data \cr
#' }
#'
#' Slot: **@info**\cr
#'
#' \tabular{lll}{
#' **OBJECT** \tab **TYPE** \tab **COMMENT** \cr
#' `call` \tab `call` \tab the original function call \cr
#' `args` \tab `list` \tab arguments of the original function call \cr
#' }
#'
#'*Note: If the input object is a [list] a list of [RLum.Results-class] objects is returned.*
#'
#'@author Sebastian Kreutzer, Institute of Geography, Heidelberg University (Germany) \cr
#' Christoph Burow, University of Cologne (Germany)
#'
#'@seealso [RLum.Results-class], [analyse_SAR.CWOSL]
#'
#'@examples
#'
#'#load data example data
#'data(ExampleData.BINfileData, envir = environment())
#'
#'#transform the values from the first position in a RLum.Analysis object
#'object <- Risoe.BINfileData2RLum.Analysis(CWOSL.SAR.Data, pos=1)
#'
##perform SAR analysis
#' results <- analyse_SAR.CWOSL(
#'   object = object,
#'   signal.integral.min = 1,
#'   signal.integral.max = 2,
#'    background.integral.min = 900,
#'    background.integral.max = 1000,
#'    plot = FALSE
#'  )
#'
#'##plot only DRC
#'plot_DRCSummary(results)
#'
#'@md
#'@export
plot_DRCSummary <- function(
  object,
  source_dose_rate = NULL,
  sel_curves = NULL,
  show_dose_points = FALSE,
  show_natural = FALSE,
  n = 51L,
  ...
) {
  .set_function_name("plot_DRCSummary")
  on.exit(.unset_function_name(), add = TRUE)

# Self-call -----------------------------------------------------------------------------------
if(inherits(object, "list")){

  ##catch ... arguments
  plot_settings <- list(...)

  ## expand input arguments
  if("main" %in% names(list(...))){
    main <- .listify(list(...)[["main"]], length(object))

    ##filter main from the ... argument list otherwise we will have a collusion
    plot_settings["main" %in% names(plot_settings)] <- NULL

  }else{
    main <- .listify("DRC", length(object))
  }

  results <- lapply(seq_along(object), function(o) {
    plot_DRCSummary(
      object = object[[o]],
      source_dose_rate = source_dose_rate,
      sel_curves = sel_curves,
      show_dose_points = show_dose_points,
      show_natural = show_natural,
      n = n,
      main = main[[o]],
      ... = plot_settings
      )
  })

  ##return merged object
  return(results)
}

  ## Integrity checks -------------------------------------------------------
  .validate_class(object, "RLum.Results")
  .validate_args(object@originator,
                 c("analyse_SAR.CWOSL", "analyse_pIRIRSequence"),
                 name = "Object originator")

# Extract data from object --------------------------------------------------------------------

  ## set limit
  if (is.null(sel_curves)) {
    sel_curves <- 1:length(object@data$Formula)
  } else {
      if(min(sel_curves) < 1 ||
         max(sel_curves) > length(object@data$Formula) ||
         length(sel_curves) > length(object@data$Formula)){
        .throw_warning("'sel_curves' out of bounds, reset to full dataset")
        sel_curves <- 1:length(object@data$Formula)
      }
  }

    ## check the whether the fitting was all the same
    if(length(unique(object@data[["data"]][["Fit"]])) != 1)
      .throw_error("I can only visualise dose-response curves based ",
                   "on the same fitting equation")

    ##get DRC
    DRC <- object@data$Formula[sel_curves]

    ## check for OTOR fit option (we can only do all )
    if(all(object@data$data[["Fit"]] %in% c("OTOR", "OTORX")))
      W <- lamW::lambertW0

  ## get limits for each set
  idx.natural <- which(object@data$LnLxTnTx.table[["Name"]] == "Natural")
  dataset_limits <- cbind(idx.natural,
                          c(idx.natural[-1] - 1, length(idx.natural)))

   ##create list
   LxTx <- lapply(1:nrow(dataset_limits), function(x){
     object@data$LnLxTnTx.table[dataset_limits[x,1]:dataset_limits[x,2],]
   })[sel_curves]

# Plotting ------------------------------------------------------------------------------------
  ##set default
  plot_settings <-  modifyList(x = list(
    xlab = if(is.null(source_dose_rate)) {"Dose [s]"} else {"Dose [Gy]"},
    ylab = expression(L[x]/T[x]),
    xlim = c(0,max(vapply(LxTx, function(x){max(x[["Dose"]])}, numeric(1)))),
    ylim = if(show_dose_points){
      c(0,max(vapply(LxTx, function(x){max(x[["LxTx"]] + x[["LxTx.Error"]])}, numeric(1)), na.rm = TRUE))
    }else{
      c(0,max(vapply(1:length(LxTx), function(y){
        x <- max(LxTx[[y]][["Dose"]], na.rm = TRUE)
        eval(DRC[[y]])
       },numeric(1)), na.rm = TRUE))
    },
    main = "DRC Summary",
    mtext = paste0("n_curves: ",length(sel_curves)),
    lty = 1,
    lwd = 1,
    pch = 20,
    col.lty = rgb(0,0,0,0.5),
    col.pch = rgb(0,0,0,0.5)
  ), val = list(...), keep.null = TRUE)

  ## expand parameters
  plot_settings$col.lty <- rep(plot_settings$col.lty, length(sel_curves))
  plot_settings$col.pch <- rep(plot_settings$col.pch, length(sel_curves))
  plot_settings$pch <- rep(plot_settings$pch, length(sel_curves))
  plot_settings$lty <- rep(plot_settings$lty, length(sel_curves))

  ##create empty plot window
  plot(
    x = NA,
    y = NA,
    xlab = plot_settings$xlab,
    ylab = plot_settings$ylab,
    xlim = plot_settings$xlim,
    ylim = plot_settings$ylim,
    main = plot_settings$main,
    xaxt = "n"
  )

  if(!is.null(plot_settings$mtext))
    mtext(side = 3, text = plot_settings$mtext, cex = 0.8)

  #exchange x-axis if source dose rate is set
  if(!is.null(source_dose_rate)){
    axis(side = 1, at = axTicks(side = 1), labels = round(axTicks(side = 1) * source_dose_rate[1],0))

  }else{
    axis(side = 1)
  }

  for(i in 1:length(sel_curves)){
    ##plot natural
    if(show_natural){
      segments(x0 = LxTx[[i]]$Dose[1], x1 = LxTx[[i]]$Dose[1],
               y0 = LxTx[[i]]$LxTx[1] - LxTx[[i]]$LxTx.Error[1],
               y1 = LxTx[[i]]$LxTx[1] + LxTx[[i]]$LxTx.Error[1],
               col = plot_settings$col.pch[[i]])
      points(
        x = LxTx[[i]]$Dose[1],
        y = LxTx[[i]]$LxTx[1],
        col = plot_settings$col.pch[[i]],
        pch = plot_settings$pch[[i]]
      )
    }

    ##plot dose points
    if(show_dose_points){
      segments(x0 = LxTx[[i]]$Dose[-1], x1 = LxTx[[i]]$Dose[-1],
               y0 = LxTx[[i]]$LxTx[-1] - LxTx[[i]]$LxTx.Error[-1],
               y1 = LxTx[[i]]$LxTx[-1] + LxTx[[i]]$LxTx.Error[-1],
               col = plot_settings$col.pch[[i]])
      points(
        x = LxTx[[i]]$Dose[-1],
        y = LxTx[[i]]$LxTx[-1],
        col = plot_settings$col.pch[[i]],
        pch = plot_settings$pch[[i]]
      )
    }

    ##plot lines
    x <- seq(min(plot_settings$xlim),max(plot_settings$xlim), length.out = n)
    y <- eval(DRC[[i]])

    if (anyNA(y) || any(is.nan(y))) {
      .throw_warning("Dose response curve ", i, " contains NA/NaN values, ",
                     "curve removed before plotting")
      next
    }

    lines(
      x = x,
      y = eval(DRC[[i]]),
      col = plot_settings$col.lty[[i]],
      lwd = plot_settings$lwd,
      lty = plot_settings$lty[[i]]
    )
  }

  ## Results -------------------------------------------------------------------
  results <- set_RLum(
    class = "RLum.Results",
    data = list(
      results = data.frame(
        dose = x,
        sapply(DRC, function(d, n) { eval(d) }, n)
        ),
      data = object
    ),
    info = list(
      call = sys.call(),
      args = as.list(sys.call())[-1])
  )

  ## Return value --------------------------------------------------------------
  return(results)
}