R/utils.R

In MicrobialGrowth: Estimates Growth Parameters from Models and Plots the Curve

# Copyright 2023, Florentin L'Homme, Clarisse Breard
#
# This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with this program. If not, see <https://www.gnu.org/licenses/>


#' Threshold when little data is used in a regression
#'
#' Number of data below which the usual methods for choosing starting values (`start`) and limits (`lower` and `upper`) will not be used in favor of a secondary method more suited to the low number of data.
#'
#' @export
THRESHOLD_FEW_DATA = 20


#' Regression function name getter
#'
#' Returns the name of the regression function associated with the model.
#'
#' @param model the model name.
#'
#' @return the string corresponding to the regression function of the model.
#' Warning, this function does not check the existence of the corresponding function.
#'
#' @export
#'
#' @examples
#' getFunctionName("gompertz")
#' ## [1] ".MicrobialGrowth.gompertz"
#'
#' # Note that this does not verify the existence
#' getFunctionName("NonExistentFunction")
#' ## [1] ".MicrobialGrowth.NonExistentFunction"
getFunctionName = function(model) {
  return(paste(".MicrobialGrowth", model, sep="."))
}


#' Create function name getter
#'
#' Returns the name of the creation function associated with the model.
#'
#' @param model the model name.
#'
#' @return the string corresponding to the creation function of the model.
#' Warning, this function does not check the existence of the corresponding function.
#'
#' @export
#'
#' @examples
#' getCreateFunctionName("gompertz")
#' ## [1] "gompertz.create"
#'
#' # Note that this does not verify the existence
#' getCreateFunctionName("NonExistentFunction")
#' ## [1] "NonExistentFunction.create"
getCreateFunctionName = function(model) {
  return(paste(model, "create", sep="."))
}


#' List available models.
#'
#' @return the list of available models.
#'
#' @export
#'
#' @details lists the models by scanning the available ".MicrobialGrowth.m" regression functions, with "m" the name of the model.
#'
#' @examples
#' listAvailableModels()
#' ## [1] "baranyi"  "gompertz"  "rosso"
listAvailableModels = function() {
  functions = ls("package:MicrobialGrowth", all.names = TRUE) # to replace with `ls("package:MicrobialGrowth", all.names = TRUE)`
  available = unlist(Map(function(x){x[3]}, strsplit(functions[startsWith(functions, ".MicrobialGrowth.")], "\\.")))
  return(available)
}


#' Model name getter
#'
#' Returns the name of the model used.
#'
#' @param x a MicrobialGrowth-object
#'
#' @return the name of the model used.
#'
#' @export
#'
#' @details scans the classes of the object which must correspond on the one hand to the generic class "MicrobialGrowth" and on the other hand to the class-model. It is this second that is returned.
#'
#' @examples
#' g <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")
#' getModelName(g)
#' ## [1] "gompertz"
getModelName = function(x) {
  if (length(class(x)) != 2) warning("Number of class different from what is expected, the function may return an incorrect value");
  return(class(x)[1])
}


#' Formula getter
#'
#' Returns the formula associated with the specified model.
#'
#' @param model the model name.
#'
#' @return the function corresponding to the formula of the model.
#'
#' @importFrom methods getFunction
#' @export
#'
#' @examples
#' f <- getFormula("gompertz")
#' # We need to set the parameters (N0, ..., lambda)
#' f2 <- f(0.1, 2, 0.2, 5)
#' # And we can then use the function "f(x)" with x the time
#' f2(4)
#' ## [1] 0.1150952
#'
#' # The same, more direct
#' F <- getFormula("gompertz")(0.1, 2, 0.2, 5)
#' F(4)
#' ## [1] 0.1150952
getFormula = function(model) {
  fName = paste(".", model, ".formula", sep = "")
  if (fName %in% ls("package:MicrobialGrowth", all.names = TRUE)) {
    return(getFunction(fName))
  }
  stop(sprintf("No formula found for %s model", model))
}



#' Model Integrity Checker
#'
#' Function to check the integrity of a given model. Used only for development.
#'
#' @param model the model to check.
#' @param verbose boolean indicating if the function is verbose, i.e. it indicates the different steps that it validates. If `FALSE`, only warnings and errors will be reported.
#'
#' @return No return value, called to check the integrity of a (new) model for the package (raises an error if the model is invalid).
#'
#' @importFrom utils hasName
#' @export
#'
#' @examples
#' # Auto-run on package build
#' models = listAvailableModels()
#' for (model in models) {
#'   .checkModelIntegrity(model)
#' }
.checkModelIntegrity = function(model, verbose = TRUE) {
  # Regression function
  if (!(model %in% listAvailableModels())) { stop(sprintf("The %s model isn't available with MicrobialGrowth", model)) }
  if (verbose) { cat(sprintf("%s model is available with MicrobialGrowth\n", model)) }

  # Expression function
  getFormula(model)
  if (verbose) { cat(sprintf("Formula function found for %s model\n", model)) }

  # Create function
  if (!exists(getCreateFunctionName(model))) { stop(sprintf("Create function not found for %s model", model)) }
  if (verbose) { cat(sprintf("Create function found for %s model\n", model)) }

  # "is." function
  if (!exists(paste("is", model, sep="."))) { stop(sprintf("Is function not found for %s model", model)) }
  # TODO : Tester pour voir si ça prend bien en compte un objet retourné (qui doit renvoyer TRUE évidemment) et au moins une liste d'objet (=list(obj, obj) => TRUE)
  if (verbose) { cat(sprintf("Is function found for %s model\n", model)) }

  # Structure of objects returned
  x = c(0.00, 5.26, 10.53, 15.79, 21.05, 26.32, 31.58, 36.84, 42.11, 47.37, 52.63, 57.89, 63.16, 68.42, 73.68, 78.95, 84.21, 89.47, 94.74, 100.00)
  y = c(0.15, 0.15, 0.15, 0.16, 0.19, 0.26, 0.38, 0.58, 0.85, 1.18, 1.53, 1.86, 2.15, 2.38, 2.55, 2.66, 2.78, 2.85, 2.89, 2.93)
  obj = do.call(getFunctionName(model), list(x=x, y=y))
  #TODO : faire 2 boucles ? Une avec l'espoir d'avoir un isValid, l'autre non ?
  if (is.null(obj)) { stop(sprintf("NULL is returned for %s model", model)) }
  if (!hasName(obj, "data")) { stop(sprintf("Object returned for %s model doesn't have \"data\" member containing data used for the regression", model)) }
  if (!hasName(obj$data, "x") | !hasName(obj$data, "y")) { stop(sprintf("Missing \"x\" or \"y\" in \"data\" member for %s model", model)) }
  if (!hasName(obj, "call")) { stop(sprintf("Object returned for %s model doesn't have \"call\" member containing the call and parameters used", model)) }
  callParams = c("x", "y", "clip", "start", "lower", "upper", "nls.args")
  if (!all(hasName(obj$call, callParams))) { stop(sprintf("The \"call\" member for %s model must have at least the parameters: %s", model, paste(callParams, collapse = ", "))) }
  if (!hasName(obj, "reg")) { stop(sprintf("Object returned for %s model doesn't have \"reg\" member corresponding to the regression result", model)) }
  if (!hasName(obj, "coefficients")) { stop(sprintf("Object returned for %s model doesn't have \"coefficients\" member referencing the parameters regressed and their value", model)) }
  if (!hasName(obj, "isValid")) { stop(sprintf("Object returned for %s model doesn't have \"isValid\" member indicating the validity of the result (i.e. regression successfull or user-defined object)", model)) }
  if (!hasName(obj, "f")) { stop(sprintf("Object returned for %s model doesn't have \"f\" member containing utils functions", model)) }
  if (!all(sapply(obj$f, is.function))) { stop(sprintf("All members of \"f\" member of object returned for %s model must be a function", model)) }
  functions = c("confint", "confint.lower", "confint.upper", "formula", "doublingTime")
  if (!all(hasName(obj$f, functions))) { stop(sprintf("The \"f\" member for %s model must have at least the functions: %s", model, paste(functions, collapse = ", "))) }
  if (verbose) { cat(sprintf("Object returned with %s model seems to be OK (structure analysed, not content)\n", model)) }


  # Print/Plot/Summary ?
  for (f in c("print", "plot", "summary")) {
    if (!exists(paste(f, model, sep="."))) { warning(sprintf("%s function not found for %s model, default %s.MicrobialGrowth function will be used", f, model, f)) }
    else if (verbose) { cat(sprintf("%s function found and will be used for %s model\n", paste(f, model, sep="."), model)) }
  }

  if (verbose) { cat(sprintf("checkModelIntegrity done for %s model\n", model)) }
}



#' MicrobialGrowth object
#'
#' Provide the skeleton of the MicrobialGrowth object.
#' Must be completed for each model.
#'
#' @param ... further arguments passed to or from other methods.
#'
#' @return a MicrobialGrowth object skeleton.
#' @export
#'
#' @details the three dots `...` are passed to \link[base]{new.env} function.
#'
#' @examples
#' # First, create the skeleton.
#' model.object = .new.MicrobialGrowth.core()
#'
#' # Then complete with data, functions, etc.
#' model.object$data$x = c(1,2,3)
#' model.object$data$y = c(1,2,3)
#' model.object$coefficients = list(N0 = 0, Nmax=0, mu=0, lambda=0)
#' model.object$f$formula = function(x){ return(x) }
#' model.object$f$confint.lower = function(x){ return(x - 1) }
#' model.object$f$confint.upper = function(x){ return(x + 1) }
#'
#' # Specialize the object by adding a class name at first position.
#' class(model.object) = c("specialized.model", class(model.object))
#'
#' # You can print, plot, etc., with the generic functions of MicrobialGrowth.
#' print(model.object)
#' ##MicrobialGrowth, model specialized.model:
#' ##    N0   Nmax     mu lambda
#' ##     0      0      0      0
#' plot(model.object)
#'
#' # Don't forget to change `isValid` to TRUE to confirm the success of the regression.
#' model.object$isValid = TRUE # Not a good idea here, since we have no `reg` value.
.new.MicrobialGrowth.core = function(...) {
  env = new.env(...)
  class(env) = "MicrobialGrowth"
  env$f = list()
  env$data = list(x=c(), y=c())
  env$call = NULL
  env$reg = NULL
  env$coefficients = list(N0 = NA, Nmax = NA, mu = NA, lambda = NA)
  env$message = NULL
  env$isValid = FALSE

  env
}



#' Check MicrobialGrowth arguments (regression function)
#'
#' Check the arguments passed to the regression function.
#' Tests are generic for all models. For example, the same length for x and y, the type of the
#' different arguments, the values inserted in `start`, `lower` and `upper`, etc. are tested
#'
#' @param x index series or time series.
#' @param y values or list of values to regress (should not be logged).
#' @param clip a pair of values indicating in which interval to clip the data `y`. When `clip` is missing, default values are used.
#' @param start a named list of starting estimates. When `start` is missing, default values are used.
#' @param lower a named list of lower bounds. When `lower` is missing, default values are used.
#' @param upper a named list of upper bounds. When `upper` is missing, default values are used.
#' @param nls.args additional parameters to use when calling \link[stats]{nls}.
#' @param callbackError function to call on error during regression.
#'
#' @return the modified `clip` value and raises an error if something is wrong.
#' @export
#'
#' @details During the check, the `clip` value is also updated. If the lower bound of `clip` is `-Inf` (default value), then this value is replaced by the smallest value greater than zero found in `y`.
.checkMicrobialGrowthArgs = function(x, y, clip, start, lower, upper, nls.args, callbackError) {
  if (length(x) != length(y)) { stop("'x' and 'y' lengths differ") }
  if (length(clip) != 2) { stop("Clip argument must have strictly 2 values") }
  if (clip[1] == -Inf) { clip = c(min(y[y > 0], na.rm = TRUE), clip[2]) }
  if (clip[1] < 0) { stop("Clip min. value must be greater than zero") }
  if (clip[2] <= clip[1]) { stop("Clip min. value must be strictly lower than clip max. value") }
  if (!is.list(start)) { stop("start argument must be a list") }
  if (!is.list(lower)) { stop("lower argument must be a list") }
  if (!is.list(upper)) { stop("upper argument must be a list") }
  if (!all(names(start) %in% c("N0", "Nmax", "mu", "lambda"))) { stop(paste("start argument contains unknown names :", paste(names(start)[!(names(start) %in% c("N0", "Nmax", "mu", "lambda"))], collapse = ", ") )) }
  if (!all(names(lower) %in% c("N0", "Nmax", "mu", "lambda"))) { stop(paste("lower argument contains unknown names :", paste(names(lower)[!(names(lower) %in% c("N0", "Nmax", "mu", "lambda"))], collapse = ", ") )) }
  if (!all(names(upper) %in% c("N0", "Nmax", "mu", "lambda"))) { stop(paste("upper argument contains unknown names :", paste(names(upper)[!(names(upper) %in% c("N0", "Nmax", "mu", "lambda"))], collapse = ", ") )) }
  if (!is.null(callbackError) && !is.function(callbackError)) { stop("callbackError argument must be NULL or a function") }
  return(clip)
}



#' Check MicrobialGrowth arguments (create function)
#'
#' Check the arguments passed to the create function.
#' Tests are generic for all models. For example, the type and value `xlim`,
#' the values of N0, Nmax, etc. are tested.
#'
#' @param N0 initial population.
#' @param Nmax final/maximum population.
#' @param mu growth rate.
#' @param lambda latency time.
#' @param xlim range of values to simulate `x` and `y` (and hence plotting, *etc.*)
#' @param n number of points to simulate in the interval `xlim`.
#'
#' @return raise an error if something is wrong.
#' @export
.checkMicrobialGrowthCreateArgs = function(N0, Nmax, mu, lambda, xlim, n) {
  if (length(xlim) != 2 || !is.numeric(xlim)) { stop("Argument xlim must have 2 numeric values") }
  if (!all(N0 > 0)) { stop("Argument N0 must be greater than 0") }
  if (!(min(Nmax) > max(N0))) { stop("Argument Nmax must be greater than N0") }
  if (!all(mu >= 0)) { stop("Argument mu must be greater or equal to 0") }
  if (!all(lambda >= 0)) { stop("Argument lambda must be greater or equal to 0") }
}



#' Coefficient argument parser (create function)
#'
#' Parses the coefficients passed to the create function to obtain  3 values: one for the main curve and two for the confint curves.
#' These values are sorted.
#'
#' @param x value(s) for a given coefficient.
#'
#' @return the 3 ordered values for the given coefficient.
#' @export
#'
#' @examples
#' .parseMicrobialGrowthCreateArgs(1)
#' ## [1] 1 1 1
#'
#' .parseMicrobialGrowthCreateArgs(c(1,2))
#' ## [1] 1.0 1.5 2.0
#'
#' .parseMicrobialGrowthCreateArgs(c(1,2,3))
#' ## [1] 1 2 3
#'
#' .parseMicrobialGrowthCreateArgs(c(3,1,2))
#' ## [1] 1 2 3
.parseMicrobialGrowthCreateArgs = function(x) {
  argName = deparse(substitute(x))
  if (!is.numeric(x)) {
    stop(paste("Argument", argName,"must be of numeric type, type", typeof(x),"found"))
  }
  x = sort(x)
  if (length(x) == 3) {
    return(x)
  }
  if (length(x) == 2) {
    return( c(x[1], mean(x), x[2]) )
  }
  if (length(x) == 1) {
    return( c(x, x, x) )
  }
  stop(paste("Argument", argName,"must have 1 to 3 elements,", length(x),"found"))
}



#' Default NLS values
#'
#' Gives default values for NLS regression from `x` and `y` values.
#' The method of calculating default values differs depending on the amount of data available in `y`.
#' Default values can be pre-set by providing them in the `start`, `lower` and `upper` arguments.
#'
#' @param x index series or time series.
#' @param y values or list of values to regress (should not be logged, must be strictly greater than zero).
#' @param start a named list of starting estimates. The coefficients specified in this list will not be calculated.
#' @param lower a named list of lower bounds. The coefficients specified in this list will not be calculated.
#' @param upper a named list of upper bounds. The coefficients specified in this list will not be calculated.
#'
#' @return the default values of `start`, `lower` and `upper` for NLS regression.
#' @export
#'
#' @details default values are calculated as follows:
#' * `start`
#'   * `N0`: the minimum value of `y`
#'   * `Nmax`: the maximum value of `y`
#'   * `mu`:
#'     * `if length(y) <= `\link{THRESHOLD_FEW_DATA}: the greatest slope between two contiguous points (on logged `y` values)
#'     * `else`: the linear regression on data positioned in the middle ±25% of the amplitude on logged y
#'   * `lambda`: the highest value of `x` which is within the lowest 5% of amplitude of `y`
#' * `lower`
#'   * `N0`: the smallest value greater than zero calculated with `1/.Machine$double.xmax`
#'   * `Nmax`: the mean value of `y`
#'   * `mu`: the amplitude on `y` divided by the amplitude on `x`
#'   * `lambda`:the minimum value of `x`
#' * `upper`
#'   * `N0`: the mean value of `y`
#'   * `Nmax`: twice the max value of `y`
#'   * `mu`:
#'     * `if length(y) <= `\link{THRESHOLD_FEW_DATA}: the amplitude on logged `y` divided by the smallest step between two contiguous `x` values
#'     * `else`: the greatest slope between two contiguous points (on logged `y` values)
#'   * `lambda`: the maximum value of `x`
#'
#' Note that it is possible, particularly when there is little data, that linear regression for `start$mu` is not possible, hence the presence of condition with \link{THRESHOLD_FEW_DATA}.
#'
#' @examples
#' # Example data
#' x = c(0.00, 5.26, 10.53, 15.79, 21.05, 26.32, 31.58, 36.84, 42.11, 47.37, 52.63,
#'       57.89, 63.16, 68.42, 73.68, 78.95, 84.21, 89.47, 94.74, 100.00)
#' y = c(0.15, 0.15, 0.15, 0.16, 0.19, 0.26, 0.38, 0.58, 0.85, 1.18, 1.53, 1.86,
#'       2.15, 2.38, 2.55, 2.66, 2.78, 2.85, 2.89, 2.93)
#'
#' # Simple example
#' values = .getDefaultNlsValues(x, y)
#' cat("N0=", values$start$N0, " with limits [", values$lower$N0, ", ", values$upper$N0,"]", sep="")
#' ## N0=0.15 with limits [5.562685e-309, 1.4315]
#'
#' # Example with specifying a starting value (which will therefore not be calculated)
#' values = .getDefaultNlsValues(x, y, start=list(N0=0.1))
#' cat("N0=", values$start$N0, " with limits [", values$lower$N0, ", ", values$upper$N0,"]", sep="")
#' ## N0=0.1 with limits [5.562685e-309, 1.4315]
.getDefaultNlsValues = function(x, y, start = list(), lower = list(), upper = list()) {
  logY = log(y)

  .start = start
  start = list()
  logStartN0 = mean(logY[logY <= (max(logY)-min(logY))*0.2+min(logY)])
  start$N0 = if("N0" %in% names(.start)) .start[["N0"]] else exp(logStartN0)   # Moyenne des 20% valeurs basses de log(y)
  logStartNmax = mean(logY[logY >= (max(logY)-min(logY))*0.8+min(logY)])
  start$Nmax = if("Nmax" %in% names(.start)) .start[["Nmax"]] else exp(logStartNmax)   # Moyenne des 20% valeurs hautes de log(y)
  start$lambda = if("lambda" %in% names(.start)) .start[["lambda"]] else x[tail(which(logY <= logStartN0), 1)]
  GrowthEnd = x[which(logY >= logStartNmax)[1]]
  start$mu = if("mu" %in% names(.start)) .start[["mu"]] else {
    if (length(y) <= THRESHOLD_FEW_DATA) (logStartNmax-logStartN0)/(GrowthEnd-start$lambda)
    else {
      tryCatch(
        {
          cond = logY >= (logStartN0 + 0.25*(logStartNmax-logStartN0)) & logY <= (logStartN0 + 0.75*(logStartNmax-logStartN0)) & x >= start$lambda
          r = lm(logY[cond] ~ x[cond])$coefficients[[2]]
          if (is.na(r) | is.null(r)) stop("lm fail")
          r
        },
        error = function(e) { return(max( (logY[2:length(logY)] - logY[1:(length(logY)-1)]) / (x[2:length(x)] - x[1:(length(x)-1)]) )) }
      )
    }
  }

  lower = list(
    N0 = if("N0" %in% names(lower)) lower[["N0"]] else 1/.Machine$double.xmax,
    Nmax = if("Nmax" %in% names(lower)) lower[["Nmax"]] else exp((max(logY)-min(logY))*0.8+min(logY)),
    mu = if("mu" %in% names(lower)) lower[["mu"]] else ((max(logY)-min(logY))/(max(x)-min(x))),
    lambda = if("lambda" %in% names(lower)) lower[["lambda"]] else 0
  )

  upper = list(
    N0 = if("N0" %in% names(upper)) upper[["N0"]] else exp((max(logY)-min(logY))*0.2+min(logY)),
    Nmax = if("Nmax" %in% names(upper)) upper[["Nmax"]] else max(y)*2,
    mu = if("mu" %in% names(upper)) upper[["mu"]] else {
      if (length(y) <= THRESHOLD_FEW_DATA) (max(logY)-min(logY)) / min(x[2:length(x)] - x[1:(length(x)-1)])
      else max( (logY[2:length(logY)] - logY[1:(length(logY)-1)]) / (x[2:length(x)] - x[1:(length(x)-1)]) )
    },
    lambda = if("lambda" %in% names(upper)) upper[["lambda"]] else max(x)
  )

  if (start$mu < lower$mu) start$mu = lower$mu
  if (start$mu > upper$mu) start$mu = upper$mu

  return(list(start=start[c("N0", "Nmax", "mu", "lambda")], lower=lower[c("N0", "Nmax", "mu", "lambda")], upper=upper[c("N0", "Nmax", "mu", "lambda")]))
}