R/plot_RLumCarlo.R

In RLumCarlo: Monte-Carlo Methods for Simulating Luminescence Phenomena

#' @title Plot RLumCarlo Monte-Carlo Simulation Results
#'
#' @description Visualise 'RLumCarlo' modelling results without
#' extracting the values manually. Typically visualised are the averaged signal
#' or the number of remaining electrons, with a polygon
#' indicating modelling uncertainties.
#'
#' @details For colouring the curves, the package [khroma::khroma-package] is used to provide colours
#' that can be best distinguished, in particular by colour-blind users.
#'
#' @param object [list] of class `RLumCarlo_Model_Output` (**required**): input object to be plotted,
#' usually the required input object is generated by one of the functions starting with `run_`.
#' Alternatively a list of such objects can be provided.
#'
#' @param plot_value [character] (*with default*): type of curve value to be displayed.
#' Allowed are `mean` (the default) and `sum` (meaningful if different systems are combined).
#' `NULL` disables the value visualisation.
#'
#' @param plot_uncertainty [character] (*with default*): type of the displayed uncertainty. Allowed
#' values are `range`, `sd` (standard deviation) and `var` (variance). `NULL` disables the uncertainty
#' visualisation.
#'
#' @param FUN [function] (*optional*): own function that can be applied to the
#' y-values before normalisation and plotting
#'
#' @param norm [logical] (*with default*): normalise curve to the highest intensity value
#'
#' @param add [logical] (*with default*): allows overplotting of results by adding curves to
#' an existing plot. This argument is handled automatically if `object` is of type [list]
#'
#' @param \dots further argument, that can be passed to control the plot output largely
#' following the argument names in [graphics::plot.default]. Currently supported
#' are: `xlab`, `ylab`, `xlim`, `ylim`, `main`, `lwd`, `type`, `pch`, `lty`,`col`, `grid`, `legend`.
#' The arguments `lwd`, `type`, `pch`, `lty`, `col` can be provided as a vector if `object` is a [list]
#'
#' @return This function returns a graphical output which is the visualisation of the modelling
#' output.
#'
#' @section Function version: 0.1.0
#'
#' @author Sebastian Kreutzer, Institute of Geography, Heidelberg University (Germany)\cr
#' Johannes Friedrich, University of Bayreuth (Germany)
#'
#' @examples
#' ## plain plot
#' DELOC <- run_MC_TL_DELOC(
#'   s = 3.5e12,
#'   E = 1.45,
#'   R = 0.1,
#'   method = 'seq',
#'   clusters = 100,
#'   times = 150:350) %T>%
#' plot_RLumCarlo(legend = TRUE)
#'
#' ## TL with FUN to correct for thermal
#' ## quenching
#' f <- function(x) x * 1/(1 + (2e+6 * exp(-0.55/(8.617e-5 * (DELOC$time + 273)))))
#' plot_RLumCarlo(
#'  object = DELOC,
#'  FUN = f)
#'
#' @keywords hplot
#' @md
#' @export
plot_RLumCarlo <- function(
  object,
  plot_value = "mean",
  plot_uncertainty = "range",
  FUN = NULL,
  norm = FALSE,
  add = FALSE,
  ...){


 # Self-call -----------------------------------------------------------------------------------
 if(class(object)[1] == "list"){
    ## double check the objects in the list
    if(!all(sapply(object, class) == "RLumCarlo_Model_Output"))
      stop("[plot_RLumCarlo()] At least one element in the list is not of class RLumCarlo_Model_Output!", call. = FALSE)

    ## summarize object already here (if we use lapply we use the class information)
    for(i in 1:length(object))
      object[[i]] <- summary(object[[i]], verbose = FALSE)

    ## set plot settings for which vectors make sense
    ## all other values are either explicit arguments are they are piped via the
    ## dot_list
    plot_settings <- list(
      col = if(length(object) < 8){
        khroma::colour("bright")(7)[1:length(object)]

      } else {
        rainbow(length(object))

      },
      pch = 1,
      lty = 1,
      lwd = 2,
      type = "l",
      ylim = c(0, max(vapply(object, function(x) max(x[["y_max"]]), numeric(1))))

    )

    ## keep only values we do not have on the top
    dot_list <- list(...)
    dot_list <- dot_list[!names(dot_list) %in% names(plot_settings)]

    ## overwrite defaults
    plot_settings <- modifyList(x = plot_settings, val = list(...))

    ## expand all arguments to the length of the list
    plot_settings <- lapply(plot_settings, FUN = rep_len, length.out = length(object))

    ## set add
    add <- c(FALSE, rep(TRUE,length(object) - 1))

    ## here we always start with a fresh graph, otherwise we add
    ## the plots on the top
    opar <- par(no.readonly =TRUE)
    on.exit(par(opar))
    par(new = FALSE)

    for(i in 1:length(object)){
      do.call(what = plot_RLumCarlo, args = c(list(
        object[[i]],
        plot_value = plot_value,
        plot_uncertainty = plot_uncertainty,
        FUN = FUN,
        norm = norm,
        add = add[i],
        ylim = if(norm) c(0,1) else range(plot_settings$ylim),
        col = plot_settings$col[i],
        pch = plot_settings$pch[i],
        lty = plot_settings$lty[i],
        lwd = plot_settings$lwd[i],
        type = plot_settings$type[i]
        ),
        dot_list
      ))

    }

   return(invisible(NULL))
 }

 # Preparation ---------------------------------------------------------------------------------

 ## try to summarise whenenver possible
 if(!"RLumCarlo_Model_Output" %in% class(object))
   stop("[plot_RLumCarlo()] 'object' needs to be of class RLumCarlo_Model_Output!", call. = FALSE)

 ## summarize
 if("RLumCarlo_Model_Output" %in% class(object) && class(object)[1] != "data.frame")
   object <- summary(object, verbose = FALSE)

 ## extract correct columns for value
 avg <- y_min <- y_max <- NULL
 if(!is.null(plot_value) && plot_value %in% c("mean", "sum"))
   avg <- object[[plot_value]]

 ## extract values for the uncertainties
 if(!is.null(plot_value) && plot_value == "sum") plot_uncertainty <- NULL

 if(!is.null(plot_uncertainty) && !is.na(plot_uncertainty)){
  if(plot_uncertainty == "sd") {
    y_min <-  avg -  object[["sd"]]
    y_max <-  avg +  object[["sd"]]


  } else if (plot_uncertainty == "var") {
    y_min <-  avg -  object[["var"]]
    y_max <-  avg +  object[["var"]]


  } else {
    y_min <- object[["y_min"]]
    y_max <- object[["y_max"]]

  }

 }


  ## set times
  times <- object[["time"]]

  ## apply FUN
  if(!is.null(FUN) && is.function(FUN)){
    if(is.null(formals(FUN)))
      stop("[plot_RLumCarlo()] FUN has no argument!", call. = FALSE)

    ## apply function
    avg <- FUN(avg)

    if(!is.null(y_min) && !is.null(y_max)){
      y_min <- FUN(y_min)
      y_max <- FUN(y_max)
    }
  }

  if(norm){ # normalization
    avg <- avg / max(avg)

    if(!is.null(plot_uncertainty)){
      y_min <- y_min / max(y_min)
      y_max <- y_max / max(y_max)

    }
    ylab <- "Normalized signal"

  } else {
    ylab <- "Signal [a.u.]"

  }

  ##default plot settings
  plot_settings <-
    modifyList(x = list(
      main = "",
      xlim = range(times),
      ylim = if(is.null(y_max)) c(0, max(avg)) else c(0, max(y_max)),
      xlab = if(length(grep(pattern = "TL", attributes(object)$model, fixed = TRUE) == 1)) {
        "Temperature [\u00b0C]"
      } else {
        "Time [s]"
      },
      ylab = ylab,
      type = "l",
      lwd = 2,
      pch = 1,
      lty = 1,
      grid = TRUE,
      col = khroma::colour("bright")(7)[2],
      legend = TRUE
    ), val = list(...))

  # Plotting ------------------------------------------------------------------------------------
  ## check if plot was already called, if not just plot
    if(add == FALSE || is.null(tryCatch(par(new =TRUE), warning = function(w) NULL))){
      plot(NA, NA,
        main = plot_settings$main,
        xlab = plot_settings$xlab,
        ylab = plot_settings$ylab,
        type = plot_settings$type,
        xlim = plot_settings$xlim,
        ylim = plot_settings$ylim
      )

    }

  ##add grid
  if(plot_settings$grid) grid()

  ## draw error polygon
  if(!is.null(plot_uncertainty)){
  polygon(x = c(times, rev(times)),
          y = c(y_min, rev(y_max)),
          col = adjustcolor(plot_settings$col, alpha.f=0.3),
          border = NA)

  }

  ## add average lines
  lines(
    x = times,
    y = avg,
    col = adjustcolor(plot_settings$col, alpha.f=1),
    type = plot_settings$type,
    pch = plot_settings$pch,
    lty = plot_settings$lty,
    lwd = plot_settings$lwd)

  if(plot_settings$legend && !is.null(plot_uncertainty)){
    legend(
      "topright",
      legend = c(plot_value, plot_uncertainty),
      lwd = c(2,5),
      bty = "n",
      col = c(
        adjustcolor(plot_settings$col, alpha.f = 1),
        adjustcolor(plot_settings$col, alpha.f = 0.3)
      )
    )
  }

}