R/fit_IsothermalHolding.R

In Luminescence: Comprehensive Luminescence Dating Data Analysis

#' @title Fit Isothermal Holding Data
#'
#' @description ##TODO
#'
#' @details ##TODO
#'
#' @param data [character] or [data.frame] (**required**): file path or data
#' frame with 5 columns named "SAMPLE", "TEMP", "TIME", "LxTx", "LxTx_ERROR".
#'
#' @param rhop [numeric] or [Luminescence::RLum.Results-class] (*with default*):
#' a vector of rho prime values (one for each sample) or an [Luminescence::RLum.Results-class]
#' object produced by [Luminescence::analyse_FadingMeasurement].
#'
#' @param ITL_model [character] (*with default*): model to be fitted, either
#' `"GOK"` or `"BTS"`.
#'
#' @param normalise_LxTx [logical] (*with default*):
#' whether the `LxTx` data should be normalised to its maximum value (`TRUE`
#' by default).
#'
#' @param plot [logical] (*with default*): enable/disable the plot output.
#'
#' @param verbose [logical] (*with default*): enable/disable output to the
#' terminal.
#'
#' @param txtProgressBar [logical] (*with default*): enable/disable the
#' progress bar. Ignored if `verbose = FALSE`.
#'
#' @param trace [logical] (*with default*): enable/disable tracing during
#' the nls fitting ([minpack.lm::nlsLM]).
#'
#' @param ... further arguments and graphical parameters that will be passed
#' to the `plot` function.
#'
#' @section Function version: 0.1.0
#'
#' @author
#' Svenja Riedesel, DTU Risø (Denmark)\cr
#' Sebastian Kreutzer, F2.1 Geophysical Parametrisation/Regionalisation, LIAG - Institute for Applied Geophysics (Germany)\cr
#' Marco Colombo, Institute of Geography, Heidelberg University (Germany)
#'
#' @keywords datagen
#'
#' @return
#' The function returns an [Luminescence::RLum.Results-class] object and an *optional* plot.
#' The object returned contains the following elements:
#'
#' \item{fit}{[list] with the fitted models}
#' \item{coefs}{[data.frame] containing the fitted coefficients for the models}
#' \item{data}{[data.frame] containing the data used in the fitting process}
#'
#' It may return `NULL` if no model could be fitted.
#'
#' @seealso [Luminescence::analyse_ThermochronometryData], [Luminescence::analyse_FadingMeasurement]
#'
#' @references
#' Li, B., Li, S.H., The effect of band-tail states on the
#' thermal stability of the infrared stimulated luminescence
#' from K-feldspar, Journal of Luminescence 136 (2013) 5–10.
#' doi: 10.1016/j.jlumin.2012.08.043
#'
#' @examples
#' # example code ##TODO
#'
#' @noRd
fit_IsothermalHolding <- function(
    data,
    rhop,
    ITL_model = c("GOK", "BTS"),
    normalise_LxTx = TRUE,
    plot = TRUE,
    verbose = TRUE,
    txtProgressBar = TRUE,
    trace = FALSE,
    ...
) {
  .set_function_name("fit_IsothermalHolding")
  on.exit(.unset_function_name(), add = TRUE)

  ## TODOs
  ## - other functions for fitting need to be implemented
  ## - uncertainties are not yet considered for the fitting, because they are not
  ##   part of the input data.
  ## - documentation needs to be completed
  ## - so far non confidence intervals on the fit ... discussion needed
  ## - the rhop value has uncertainties, which are not yet considered

  .validate_class(data, c("character", "RLum.Results", "data.frame"))
  .validate_class(rhop, c("numeric", "RLum.Results"))
  ITL_model <- .validate_args(ITL_model, c("GOK", "BTS"))
  .validate_logical_scalar(normalise_LxTx)
  .validate_logical_scalar(plot)
  .validate_logical_scalar(verbose)
  .validate_logical_scalar(txtProgressBar)

  if (inherits(data[1], "character")) {
    records_ITL <- .import_ThermochronometryData(file = data, output_type = "RLum.Results")@data$ITL

  } else if (inherits(data, "RLum.Results")) {
    .validate_originator(data, ".import_ThermochronometryData")
    records_ITL <- data@data$ITL

  } else if (inherits(data, "data.frame")) {
    records_ITL <- data
  }

  if (!all(colnames(records_ITL) %in%
           c("SAMPLE", "TEMP", "TIME", "LxTx", "LxTx_ERROR"))) {
    .throw_error("'data' has the wrong column headers, please check the manual")
  }

  ## never show the progress bar if not verbose
  if (!verbose) {
    txtProgressBar <- FALSE
  }

  ###### --- Extract data from RLum.Results for ITL fitting --- #####
  ## get unique sample names; we will use this to filter the data
  sample_id <- unique(records_ITL[["SAMPLE"]])

  ## data normalisation
  if (normalise_LxTx) {
    records_ITL$LxTx <- .normalise_curve(records_ITL$LxTx, "max")
  }

  ## extract data.frames for each sample with all information
  df_raw_list <- lapply(sample_id, function(x) records_ITL[records_ITL$SAMPLE == x, ])

  ## allow to control how many random values for the s parameter should be
  ## generated when fitting the BTS model
  extraArgs <- list(...)
  num_s_values_bts <- if (ITL_model == "BTS")
                        .validate_positive_scalar(extraArgs$num_s_values_bts,
                                                  int = TRUE, null.ok = TRUE,
                                                  name = "'num_s_values_bts'") %||% 1000
                      else 1

  ## initialise the progress bar
  if (txtProgressBar) {
    num.models <- sum(sapply(df_raw_list, function(x) length(unique(x$TEMP))))
    pb <- txtProgressBar(min = 0, max = num.models * num_s_values_bts,
                         char = "=", style = 3)
  }

  ###### --- Perform ITL fitting --- #####
  # Define variables --------------------------------------------------------
  kB <- .const$kB  # Boltzmann constant (eV/K)

  ## get the rhop value from the fading measurement analysis if available
  if (inherits(rhop, "RLum.Results") && .check_originator(rhop, "analyse_FadingMeasurement"))
    rhop <- rhop@data$rho_prime[[1]]

  ## Define formulas to fit -------------------------------------------------
  ##
  ## We define each model as a function that describes the right-hand side
  ## of the formula. This allows us to use the `$value` term in the BTS model
  ## to extract the solution of the integral, which would otherwise be
  ## incorrectly interpreted by nlsLM().

  f_GOK <- function(A, b, Et, s10, isoT, x) {
    T_K <- isoT + .const$C2K
    A * exp(-rhop * log(1.8 * 3e15 * (250 + x))^3) *
      (1 - (1 - b) * 10^s10 * exp(-Et / (kB * T_K)) * x)^(1 / (1 - b))
  }
  f_BTS <- function(A, Eu, Et, s10, isoT, x) {
    T_K <- isoT + .const$C2K
    ## call C++ calculation part; which is a lot faster than
    ## repeated calls to stats::integrate
    ## an even faster implementation would use pure C++ for everything,
    ## however, then we would use minpack.lm::nls.lm() instead ... perhaps
    ## in the future
    f_BTS_cpp_part(
      x, A, Eu, s10, Et, kB, T_K, Et, rhop)
  }

  ## switch the models
  start <- switch(
    ITL_model,
    # https://github.com/GeorginaKing/OSLThermo/blob/5f32762e2cea87aac26d65aeabe9a434b2ce19a7/Stage2a_Fitparameters.m#L311
    GOK = list(A = 1, b = 4, Et = 1.4, s10 = 9),
    # https://github.com/GeorginaKing/OSLThermo/blob/5f32762e2cea87aac26d65aeabe9a434b2ce19a7/Stage2a_Fitparameters.m#L297
    BTS = list(A = 1, Eu = 0.1, Et = 1.4))

  lower <- switch(
    ITL_model,
    GOK = c(0, 0, 0, 0),
    BTS = c(0, 0, 1))

  upper <- switch(
    ITL_model,
    GOK = c(20, Inf, 3, 20),
    BTS = c(20, 0.5, 3))

  ## Fitting ----------------------------------------------------------------
  ## we have a double loop situation: we have a list with n samples, and
  ## each sample has n temperature steps
  num.fits <- 0
  fitted.coefs <- NULL
  fit_list <- lapply(df_raw_list, function(s) {

    ## extract temperatures
    isoT <- unique(s$TEMP)

    ## run the fitting at each temperature
    tmp <- lapply(isoT, function(isoT) {

      ## extract data to fit
      tmp_fitdata <- s[s$TEMP == isoT,]
      df <- data.frame(x = tmp_fitdata$TIME,
                       y = tmp_fitdata$LxTx)

      if (ITL_model == "GOK") {
        fit <- try({
          minpack.lm::nlsLM(
              formula = y ~ f_GOK(A, b, Et, s10, isoT, x),
              data = df,
              start = start,
              lower = lower,
              upper = upper,
              control = list(
                  maxiter = 500
              ),
              trace = trace)
        }, silent = TRUE)

        coefs <- if (!inherits(fit, "try-error")) {
                   coef(fit)
                 } else {
                   stats::setNames(rep(NA_real_, length(start)), names(start))
                 }
        fitted.coefs <<- rbind(fitted.coefs,
                               data.frame(SAMPLE = unique(s$SAMPLE),
                                          TEMP = isoT,
                                          t(coefs)))

        ## update the progress bar
        num.fits <<- num.fits + 1  # uses <<- as we are in a nested lapply()
        if (txtProgressBar) {
          setTxtProgressBar(pb, num.fits)
        }

        if (inherits(fit, "try-error"))
          fit <- NA

      } else if (ITL_model == "BTS") {
        ## run fitting with different start parameters for s10
        ## all.s10 <- rnorm(num_s_values_bts, mean = 9, sd = 1.5)
        all.s10 <- stats::qnorm(seq(0.01, 0.99, length.out = num_s_values_bts),
                                mean = 9, sd = 1.5)
        fit <- lapply(1:length(all.s10), function(idx) {
          s10 <- all.s10[idx]
          t <- try(minpack.lm::nlsLM(
                       formula = y ~ f_BTS(A, Eu, Et, s10, isoT, x),
                       data = df,
                       start = start,
                       lower = lower,
                       upper = upper,
                       control = list(
                           maxiter = 500
                       ), trace = FALSE),
                   silent = TRUE)

          ## update the progress bar
          num.fits <<- num.fits + 1  # uses <<- as we are in a nested lapply()
          if (txtProgressBar) {
            setTxtProgressBar(pb, num.fits)
          }

          if (inherits(t, "try-error"))
            return(NULL)
          return(t)
        })

        ## pick the one with the best fit after removing those that didn't fit
        fit <- .rm_NULL_elements(fit)
        if (length(fit) == 0)
          return(NA)

        fit <- fit[[which.min(vapply(fit, stats::deviance, numeric(1)))]]
        s10 <- environment(fit$m$predict)$env$s10
        fitted.coefs <<- rbind(fitted.coefs,
                               data.frame(SAMPLE = unique(s$SAMPLE),
                                          TEMP = isoT,
                                          t(coef(fit)), s10))
      }

      ## return fit
      return(fit)
    })

    ## add temperature as name to list
    names(tmp) <- isoT
    return(tmp)
  })

  ## add sample names
  names(fit_list) <- sample_id

  if (txtProgressBar) {
    close(pb)
  }

  if (is.null(fitted.coefs)) {
    .throw_warning("No ITL model could be fitted")
    return(NULL)
  }

  ## summarize the parameters of interest (Et, s10) in terms of mean, median and
  ## interquartile range
  ITL_params <- tapply(fitted.coefs, sample_id, function(coefs) {
    data.frame(t(c(Et = c(mean = mean(coefs$Et, na.rm = TRUE),
                          quantile(coefs$Et, c(0.5, 0.25, 0.75), na.rm = TRUE)),
                   s10 = c(mean = mean(coefs$s10, na.rm = TRUE),
                           quantile(coefs$s10, c(0.5, 0.25, 0.75), na.rm = TRUE)))))
  })

  ## sort the ITL_params list, as tapply() may have put the results in a
  ## different order when converting internally the sample names to a factor
  ITL_params <- ITL_params[sample_id]

  ## convert to a data table
  ITL_params <- rbindlist(ITL_params[sample_id], idcol = "SAMPLE")
  colnames(ITL_params)[-1] <- c("Et_mean", "Et_median", "Et_Q_0.25", "Et_Q_0.75",
                                "s10_mean", "s10_median", "s10_Q_0.25", "s10_Q_0.75")

  if (verbose) {
    ## silence notes raised by R CMD check
    Et <- Et_Q_0.25 <- Et_Q_0.75 <- Et_median <- Et_mean <- NULL
    s10 <- s10_Q_0.25 <- s10_Q_0.75 <- s10_median <- s10_mean <- NULL

    ## report as mean; median (IQR)
    format.out <- function(x, x1, x2, x3, n = 3)
      sprintf("%.*f;   %.*f (%.*f, %.*f)", n, x, n, x1, n, x2, n, x3)
    ITL_params[, Et := format.out(Et_mean, Et_median, Et_Q_0.25, Et_Q_0.75)]
    ITL_params[, s10 := format.out(s10_mean, s10_median, s10_Q_0.25, s10_Q_0.75)]

    fmt <- "%20s | %30s | %32s |\n"
    cat("\n---- Isothermal holding parameters [mean; median (IQR)] ----\n\n")
    cat(sprintf(fmt, "SAMPLE", "Et", "log10(s)"))
    for (i in seq(nrow(ITL_params))) {
      row <- ITL_params[i, ]
      cat(sprintf(fmt, row$SAMPLE, row$Et, row$s10))
    }
    cat("\n")

    ITL_params[, c("Et", "s10") := NULL]
  }

  ## Plotting ---------------------------------------------------------------
  if (plot) {

    ## define plot settings
    plot_settings <- modifyList(
      x = list(
        xlim = range(lapply(df_raw_list, function(x) range(x$TIME))) * c(0.5, 10),
        ylim = c(0, 1.03),
        log = "x",
        xlab = "Isothermal holding time [s]",
        ylab = expression(paste("Norm. lumin. [", L[x]/T[x], "]")),
        pch = 21,
        col = grDevices::palette.colors(),
        col.border = "black",
        main = sample_id,
        mtext = paste("Fitted with the", ITL_model, "model"),
        cex = 1.0,
        mfrow = if (length(unique(sample_id)) > 1) c(min(c(2, ceiling(length(sample_id)/2))),2) else NULL,
        legend = TRUE,
        legend.pos = "bottomleft",
        legend.cex = 0.6
      ),
      val = list(...)
    )

    ## get x vector for the prediction later; we use the same
    ## value for all provided data
    x <- c(1:1e+03,seq(1e+03,plot_settings$xlim[2], length.out = 10000))

    ## par settings
    ## the check for mfrow not being null prevents mfrow from being restored on
    ## exit when this function is called from analyse_ThermochronometryData(),
    ## because restoring mfrow (even if to the same value) causes the graphical
    ## device to start a new page
    if (!is.null(plot_settings$mfrow)) {
      par.default <- .par_defaults()
      on.exit(par(par.default), add = TRUE)
      par(cex = plot_settings$cex, mfrow = plot_settings$mfrow)
    }

    ## now we have to loop over the samples
    for (i in seq_along(sample_id)) {
      ## open plot area
      plot(NA,NA,
        xlim = plot_settings$xlim,
        ylim = plot_settings$ylim,
        log = plot_settings$log,
        xlab = plot_settings$xlab,
        ylab = plot_settings$ylab,
        yaxs = "i",
        main = rep_len(plot_settings$main, length(sample_id))[i])

      ## add plot subtitle
      mtext(side = 3, plot_settings$mtext, cex = 0.7 * plot_settings$cex)

      ## add fitted curves
      isoT <- unique(df_raw_list[[i]][["TEMP"]])
      for (c in seq_along(isoT)) {
        ## only plot the fitted lines if the fit had worked
        if (inherits(fit_list[[i]][[c]], "nls")) {
          y <- stats::predict(fit_list[[i]][[c]],
                              newdata = data.frame(x = x),
                              interval = "confidence",
                              se.fit = TRUE)
          lines(
            x = x,
            y = y,
            col = rep_len(plot_settings$col[c], length(isoT)))
        }

        ## plot the points (don't use matplot because this would assume
        ## the same length for the time vector)
        df_pts <- df_raw_list[[i]][df_raw_list[[i]][["TEMP"]] == isoT[c], ]
        points(
          x = df_pts[["TIME"]],
          y = df_pts[["LxTx"]],
          xpd = TRUE,
          pch = plot_settings$pch,
          bg = plot_settings$col[c],
          col = plot_settings$col.border)

        ## add error bars (if NA nothing is plotted by R default)
        segments(
          x0 = df_pts[["TIME"]],
          x1 = df_pts[["TIME"]],
          y0 = df_pts[["LxTx"]] - df_pts[["LxTx_ERROR"]],
          y1 = df_pts[["LxTx"]] + df_pts[["LxTx_ERROR"]],
          col = plot_settings$col[c])
      }

      ## add legend
      if (plot_settings$legend[1]) {
       legend(
         plot_settings$legend.pos,
         legend = paste0(isoT, "\u00b0C"),
         bty = "n",
         pch = 20,
         lty = 1,
         col = plot_settings$col,
         cex = plot_settings$legend.cex)
      }
    }## plot loop

  }## plot condition

  # Return results ----------------------------------------------------------
  set_RLum(
    class = "RLum.Results",
    data = list(
     ITL_params = data.frame(ITL_params),
     fit = fit_list,
     coefs = fitted.coefs,
     data = records_ITL),
    info = list(
      call = sys.call())
    )
}