R/acide.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/>


#' Acid
#'
#' Modeling of an acid with `alpha` its sensitivity and `MIC` its minimum inhibition concentration.
#' The default `concentration` is 1g/L.
#'
#' @param alpha sensitivity.
#' @param MIC concentration minimale d'inhibition.
#' @param concentration acid concentration (in g/L).
#'
#' @returns the acid modeled with the following accessible attributes:
#' \item{alpha}{the `alpha` value or list of values.}
#' \item{MIC}{the `MIC` value or list of values.}
#' \item{concentration}{the acid concentration (in g/L).}
#' \item{getCoefMin}{function returning the minimum coefficient to apply to a MicrobialGrowth-object (see details section).}
#' \item{getCoefMid}{function returning the "middle" coefficient to apply to a MicrobialGrowth-object (see details section).}
#' \item{getCoefMax}{function returning the maximum coefficient to apply to a MicrobialGrowth-object (see details section).}
#'
#' @export
#'
#' @details
#' The arguments `alpha` and `MIC` can be given as one to three values.
#'
#' A single value means that `getCoefMin`, `getCoefMid` and `getCoefMax` will return the same coefficient.
#'
#' Two values symbolize some sort of uncertainty about `alpha` and/or `MIC`.
#' The functions `getCoefMin` and `getCoefMax` will use the pair (alpha, MIC) which respectively minimizes and maximizes the coefficients.
#' The `getCoefMid` function will return a coefficient based on the average of the values entered.
#'
#' Three values act as for two values except that for the function `getCoefMid` will use this third value (middle value) for the calculation of the coefficient.
#'
#' Please note, entering several values acts as a pool of available values, and not as pairs (alpha, MIC).
#' For example, the `getCoefMin` function will use the minimum value `alpha` and the minimum value `MIC`.
#' If you wish to specify pairs (alpha, MIC), see \link{Acid.SpecificPair} which will determine, for example for `getCoefMin`, the pair (alpha, MIC) minimizing the coefficient.
#'
#' @seealso \link{Acid.SpecificPair}
#'
#' @examples
#' # Classic instantiation
#' aceticAcid <- Acid(1.245, 5.47)
#' print(aceticAcid)
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=1g/L}
#'
#' # Classic instantiation by specifying a concentration
#' print( Acid(1.245, 5.47, 3) )
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=3g/L}
#'
#' # Instantiation with multiple `alpha` and `MIC` values (see details section)
#' print( Acid(c(0.98, 1.1, 1.51), c(5.26, 5.68)) )
#' ## acid {alpha=[0.98, 1.1, 1.51], MIC=[5.26, 5.68]g/L, concentration=1g/L}
#'
#' # Generic operators (`+`, `*`, etc.) can change the `concentration` of the acid.
#' print(aceticAcid / 2)
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=0.5g/L}
#' print(aceticAcid + 2)
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=3g/L}
#'
#' # Without having to pre-define specific concentrations, and with the default `concentration` (1g/L),
#' # you can dynamically change the acid concentration as follows:
#' for (concentration in c(0.5, 1, 5, 10)) {
#'   print(concentration * aceticAcid)
#' }
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=0.5g/L}
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=1g/L}
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=5g/L}
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=10g/L}
#'
#' try({
#'   # Acid can be applied to a MicrobilogicalGrowth-object with the `+` addition operator.
#'   # Note that the acid should be on the right side, otherwise an error is raised.
#'   MyMicrobialGrowthObject + aceticAcid
#'   ## returns the MicrobialGrowth-object affected by the acid (several acids can be applied)
#' })
Acid = function(alpha, MIC, concentration = 1) {
  if (length(alpha) > 3) stop(paste("Argument alpha must have 1 to 3 elements,", length(alpha),"found"))
  if (!is.numeric(alpha)) stop(paste("Argument alpha must be of numeric type, type", typeof(alpha),"found"))
  if (!all(alpha > 0)) stop("Value(s) for alpha must be greater than 0")
  if (length(MIC) > 3) stop(paste("Argument MIC must have 1 to 3 elements,", length(MIC),"found"))
  if (!is.numeric(MIC)) stop(paste("Argument MIC must be of numeric type, type", typeof(MIC),"found"))
  if (!all(MIC > 0)) stop("Value(s) for MIC must be greater than 0")
  if (length(concentration) != 1) stop(paste("Argument concentration must have only one element,", length(concentration),"found"))
  if (!is.numeric(concentration)) stop(paste("Argument concentration must be of numeric type, type", typeof(concentration),"found"))
  if (concentration < 0) stop("Value for concentration must be greater or equal to 0")

  env <- new.env()
  env$call = match.call()
  env$alpha <- alpha
  env$MIC <- MIC
  env$concentration <- concentration
  class(env) <- c("acid")

  env$getCoefMin = function() { 1 - (min(env$concentration, env$MIC)/min(env$MIC))^min(env$alpha) }
  env$getCoefMax = function() { 1 - (min(env$concentration, env$MIC)/max(env$MIC))^max(env$alpha) }
  env$getCoefMid = function() {
    alpha = if(length(env$alpha) == 1) env$alpha
            else if (length(env$alpha) == 2) mean(env$alpha)
            else sort(env$alpha)[2]
    MIC = if(length(env$MIC) == 1) env$MIC
          else if (length(env$MIC) == 2) mean(env$MIC)
          else sort(env$MIC)[2]
    1 - (min(env$concentration, MIC)/MIC)^alpha
  }

  return(env)
}


#' Acid - specific pair (`alpha`, `MIC`)
#'
#' Modeling of an acid with `alpha` its sensitivity and `MIC` its minimum inhibition concentration.
#' The default `concentration` is 1g/L.
#'
#' @param pairs list of pairs (`alpha`, `MIC`).
#' @param concentration acid concentration (in g/L).
#'
#' @returns the acid modeled with the following accessible attributes:
#' \item{pairs}{list of pairs (`alpha`, `MIC`).}
#' \item{concentration}{the acid concentration (in g/L).}
#' \item{getCoefMin}{function returning the minimum coefficient to apply to a MicrobialGrowth-object (see details section).}
#' \item{getCoefMid}{function returning the "middle" coefficient to apply to a MicrobialGrowth-object (see details section).}
#' \item{getCoefMax}{function returning the maximum coefficient to apply to a MicrobialGrowth-object (see details section).}
#'
#' @export
#'
#' @details
#' The `pairs` argument can be given as one to three pairs.
#'
#' A single pair means that `getCoefMin`, `getCoefMid` and `getCoefMax` will return the same coefficient.
#'
#' Two pairs means that one of them will be used for `getCoefMin` and the other for `getCoefMax`.
#' The `getCoefMid` function will use an average value of the two pairs.
#'
#' Three pairs acts like two pairs except that the `getCoefMid` function will use this third pair (middle value) to calculate the coefficient.
#' Note that the pair (`alpha`, `MIC`) used by `getCoefMid` neither minimizes nor maximizes the coefficient (in other words, it is the pair which is neither used in `getCoefMin` nor in `getCoefMax`).
#'
#' Please note that if you do not want to use specific pairs but ranges of values for `alpha` and/or `MIC`, use the parent function \link{Acid}.
#'
#' @seealso \link{Acid}
#'
#' @examples
#' # Classic instantiation
#' print(Acid.SpecificPair(list(c(1.245, 5.47))))
#' ## acid {{alpha=1.245, MIC=5.47g/L}, concentration=1g/L}
#'
#' # Classic instantiation by specifying a concentration
#' print(Acid.SpecificPair(list(c(1.245, 5.47)), 3))
#' ## acid {{alpha=1.245, MIC=5.47g/L}, concentration=3g/L}
#'
#' # Instantiation with multiple couples (`alpha`, `MIC`) (see details section)
#' aceticAcid <- Acid.SpecificPair(list(c(0.98,5.68),c(1.51,5.26)))
#' print(aceticAcid)
#' ## acid {{alpha=0.98, MIC=5.68g/L},
#' ##       {alpha=1.51, MIC=5.26g/L}, concentration=1g/L}
#'
#' # Generic operators (`+`, `*`, etc.) can change the `concentration` of the acid.
#' print(aceticAcid / 2)
#' ## acid {{alpha=0.98, MIC=5.68g/L},
#' ##       {alpha=1.51, MIC=5.26g/L}, concentration=0.5g/L}
#' print(aceticAcid + 2)
#' ## acid {{alpha=0.98, MIC=5.68g/L},
#' ##       {alpha=1.51, MIC=5.26g/L}, concentration=3g/L}
#'
#' # Without having to pre-define specific concentrations, and with the default `concentration` (1g/L),
#' # you can dynamically change the acid concentration as follows:
#' for (concentration in c(0.5, 1, 5, 10)) {
#'   print(concentration * aceticAcid)
#' }
#' ## acid {{alpha=0.98, MIC=5.68g/L},
#' ##       {alpha=1.51, MIC=5.26g/L}, concentration=0.5g/L}
#' ## acid {{alpha=0.98, MIC=5.68g/L},
#' ##       {alpha=1.51, MIC=5.26g/L}, concentration=1g/L}
#' ## acid {{alpha=0.98, MIC=5.68g/L},
#' ##       {alpha=1.51, MIC=5.26g/L}, concentration=5g/L}
#' ## acid {{alpha=0.98, MIC=5.68g/L},
#' ##       {alpha=1.51, MIC=5.26g/L}, concentration=10g/L}
#'
#' try({
#'   # Acid can be applied to a MicrobilogicalGrowth-object with the `+` addition operator.
#'   # Note that the acid should be on the right side, otherwise an error is raised.
#'   MyMicrobialGrowthObject + aceticAcid
#'   ## returns the MicrobialGrowth-object affected by the acid (several acids can be applied)
#' })
Acid.SpecificPair = function(pairs, concentration = 1) {
  if (!is.list(pairs)) stop(paste("Argument pairs must be of list type, type", typeof(pairs),"found"))
  if (!all(unlist(lapply(pairs, function(pair) { length(pair) == 2 & is.numeric(pair) })))) stop("All elements of the \"pairs\" list must be numeric type pairs")

  if (length(concentration) != 1) stop(paste("Argument concentration must have only one element,", length(concentration),"found"))
  if (!is.numeric(concentration)) stop(paste("Argument concentration must be of numeric type, type", typeof(concentration),"found"))
  if (concentration < 0) stop("Value for concentration must be greater or equal to 0")

  env <- new.env()
  env$call = match.call()
  env$pairs <- pairs
  env$concentration <- concentration
  class(env) <- c("acid.specific.pair", "acid")

  env$getCoefMin = function() { min( unlist(lapply(env$pairs, function(x) { 1 - (min(env$concentration, x[2])/x[2])^x[1] })) ) }
  env$getCoefMax = function() { max( unlist(lapply(env$pairs, function(x) { 1 - (min(env$concentration, x[2])/x[2])^x[1] })) ) }
  env$getCoefMid = function() {
    if (length(env$pairs) == 1) (1 - (min(env$concentration, env$pairs[[1]][2])/env$pairs[[1]][2])^env$pairs[[1]][1])
    else if (length(env$pairs) == 2) {
      alpha = (env$pairs[[1]][1] + env$pairs[[2]][1]) / 2
      MIC = (env$pairs[[1]][2] + env$pairs[[2]][2]) / 2
      (1 - (min(env$concentration, MIC)/MIC)^alpha)
    }
    else ( sort( unlist(lapply(env$pairs, function(x) { 1 - (min(env$concentration, x[2])/x[2])^x[1] })) )[2] )
  }

  return(env)
}


#' Acid print function
#'
#' Print function of Acid-object.
#'
#' @param x Acid-object.
#' @param ... further arguments passed to or from other methods.
#'
#' @return No return value, called to print information about a Acid-object.
#'
#' @exportS3Method
#'
#' @seealso \link{Acid}
#'
#' @examples
#' print( Acid(1.245, 5.47, 3) )
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=3g/L}
#'
#' print( Acid(c(0.98, 1.1, 1.51), c(5.26, 5.68)) )
#' ## acid {alpha=[0.98, 1.1, 1.51], MIC=[5.26, 5.68]g/L, concentration=1g/L}
print.acid = function(x, ...) {
  cat(paste("acid {alpha=", if(length(x$alpha)==1) x$alpha else paste("[",paste(sort(x$alpha),collapse=", ",sep=""),"]",collapse = ", ",sep=""),
            ", MIC=", if(length(x$MIC)==1) x$MIC else paste("[",paste(sort(x$MIC),collapse=", ",sep=""),"]",sep=""),
            "g/L, concentration=", x$concentration, "g/L}", sep=""), "\n", ...)
}


#' Acid.SpecificPair print function
#'
#' Print function of Acid.SpecificPair-object.
#'
#' @param x Acid.SpecificPair-object.
#' @param sep a character string to separate the different pairs.
#' @param ... further arguments passed to or from other methods.
#'
#' @return No return value, called to print information about a Acid.SpecificPair-object.
#'
#' @exportS3Method
#'
#' @seealso \link{Acid.SpecificPair}
#'
#' @examples
#' print( Acid.SpecificPair(list(c(1.245, 5.47)), 3) )
#' ## acid {{alpha=1.245, MIC=5.47g/L}, concentration=3g/L}
#'
#' print( Acid.SpecificPair(list(c(0.98, 5.68), c(1.51, 5.26))) )
#' ## acid {{alpha=0.98, MIC=5.68g/L},
#' ##       {alpha=1.51, MIC=5.26g/L}, concentration=1g/L}
print.acid.specific.pair = function(x, sep=",\n      ", ...) {
  cat(paste("acid {", paste(lapply(x$pairs, function(pair) { paste("{alpha=", pair[1], ", MIC=", pair[2], "g/L}", sep="") }), collapse = sep),
            ", concentration=", x$concentration, "g/L}", sep=""), "\n", ...)
}


#' Acid class
#'
#' Test if a variable is an Acid-object.
#'
#' @param x variable to test.
#'
#' @return `TRUE` if the object is of class `acid`.
#'
#' @export
#'
#' @examples
#' # TRUE return
#' is.acid( Acid(1.245, 5.47, 3) )
#' is.acid( Acid(c(0.98, 1.1, 1.51), c(5.26, 5.68)) )
#' # Acid.SpecificPair-objects are also Acid-objects
#' is.acid( Acid.SpecificPair(list(c(0.98, 5.68), c(1.51, 5.26))) )
#'
#' # FALSE return
#' is.acid(1)
#' is.acid( list(Acid(1.245, 5.47, 3), Acid(1.245, 5.47, 3)) )
is.acid <- function(x) .is.acid(x)
.is.acid <- function(x) UseMethod(".is.acid")
.is.acid.default <- function(x) "acid" %in% class(x)
.is.acid.acid <- function(x) TRUE


#' Acid.SpecificPair class
#'
#' Test if a variable is an Acid.SpecificPair-object.
#'
#' @param x variable to test.
#'
#' @return `TRUE` if the object is of class `acid.specific.pair`.
#'
#' @export
#'
#' @examples
#' # TRUE return
#' is.acid.specific.pair( Acid.SpecificPair(list(c(0.98, 5.68), c(1.51, 5.26))) )
#'
#' # FALSE return
#' is.acid.specific.pair(1)
#' is.acid.specific.pair( Acid(1.245, 5.47, 3) )
#' is.acid.specific.pair( list(Acid(1.245, 5.47, 3), Acid(1.245, 5.47, 3)) )
is.acid.specific.pair <- function(x) .is.acid.specific.pair(x)
.is.acid.specific.pair <- function(x) UseMethod(".is.acid.specific.pair")
.is.acid.specific.pair.default <- function(x) "acid.specific.pair" %in% class(x)
.is.acid.specific.pair.acid.specific.pair <- function(x) TRUE



#' Operators on the Acid Class
#'
#' Operators for the "\link{Acid}" class.
#'
#' @param e1 acid-object, numeric or MicrobialGrowth-object.
#' @param e2 acid-object or numeric.
#'
#' @returns a new acid or MicrobialGrowth-object.
#' @include base.R
#' @exportS3Method
#'
#' @details
#' Operations between an acid and a numeric are the most common case.
#' In this case, the operation is carried out on the `concentration` member of the acid.
#' A new acid-object is returned with the new concentration.
#'
#' Operations between acids are tolerated (but not recommended).
#' To do this, the two acids must have the same `alpha` sensitivity and the same `MIC`,
#' and the operation is carried out between the concentrations of the two acids.
#' A new acid-object is returned with the new concentration.
#'
#' The addition operator `+` can be used between MicrobialGrowth-object (left side) and an acid-object (right side).
#' This operation symbolizes the application of the acid to the MicrobialGrowth-object.
#' A new MicrobialGrowth-object is returned with its coefficients (and confidence intervals) modified by the acid.
#'
#' @examples
#' # Acids and numerics
#' print( Acid(1.245, 5.47) * 5 )
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=5g/L}
#'
#' print( Acid(1.245, 5.47) / 3 )
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=0.333333333333333g/L}
#'
#' print( 3 / Acid(1.245, 5.47) )
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=3g/L}
#'
#' print( 3 / Acid(1.245, 5.47, 0.5) )
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=6g/L}
#'
#' # Acids and acids
#' print( Acid(1.245, 5.47, 0.5) + Acid(1.245, 5.47, 2) )
#' ## acid {alpha=1.245, MIC=5.47g/L, concentration=2.5g/L}
#'
#' try({
#'   print( Acid(1.245, 5.47, 0.5) + Acid(1, 5.47, 2) )
#'   ## throws an error since `alpha` and/or `MIC` are not equal
#' })
#'
#' # Acids and MicrobialGrowth-object
#' g <- MicrobialGrowth.create(N0 = c(0.13, 0.15), Nmax = 1.43, mu = c(0.05, 0.07, 0.09),
#' lambda = c(45, 49, 43), xlim = c(0, 100), model="gompertz")
#' aceticAcid <- Acid(1.245, 5.47)
#' {
#'   cat("Before :\n")
#'   print(g)
#'   cat("After:\n")
#'   print(g + aceticAcid)
#' }
#' ## Before :
#' ## MicrobialGrowth, model gompertz:
#' ##     N0   Nmax     mu lambda
#' ##   0.14   1.43   0.07  45.00
#' ## After:
#' ## MicrobialGrowth, model gompertz:
#' ##          N0        Nmax          mu      lambda
#' ##  0.14000000  1.43000000  0.06156075 51.16896670
#'
#' # Also works with the `acid.specific.pair` subclass
#' print( Acid.SpecificPair(list(c(0.98, 5.68), c(1.51, 5.26))) )
#' ## acid {{alpha=0.98, MIC=5.68g/L},
#' ##       {alpha=1.51, MIC=5.26g/L}, concentration=6g/L}
Ops.acid = function(e1, e2) {
  if (is.acid(e1) & is.acid(e2)) {
    if (e1$alpha != e2$alpha | e1$MIC != e2$MIC) stop("Two acids with different alpha and MIC values cannot be operated together")

    acid = e1
    e1 = e1$concentration
    e2 = e2$concentration

    new.call = acid$call
    new.call$concentration = NextMethod(.Generic)
    return( eval(new.call) )
  }

  if (is.acid(e1) & is.numeric(e2)) {
    boolean <- switch(.Generic, `+` = , `-` = , `*` = , `/` = , `^` = TRUE, FALSE)
    if (!boolean) stop(gettextf("\"%s\" not defined for \"acid\" and \"numeric\" objects", .Generic))

    acid = e1
    e1 = acid$concentration

    new.call = acid$call
    new.call$concentration = NextMethod(.Generic)
    return( eval(new.call) )
  }
  if (is.numeric(e1) & is.acid(e2)) {
    boolean <- switch(.Generic, `+` = , `-` = , `*` = , `/` = , `^` = TRUE, FALSE)
    if (!boolean) stop(gettextf("\"%s\" not defined for \"acid\" and \"numeric\" objects", .Generic))

    acid = e2
    e2 = acid$concentration

    new.call = acid$call
    new.call$concentration = NextMethod(.Generic)
    return( eval(new.call) )
  }

  if (is.MicrobialGrowth(e1) & is.acid(e2)) {
    if (.Generic != "+") stop(gettextf("\"%s\" not defined for \"MicrobialGrowth\" and \"acid\" objects", .Generic))
    if (is.list(e1)) return( lapply(e1, function(e) { do.call(.Generic, list(e, e2)) }) )

    MGObj = e1
    acid = e2

    concentration = min(c(acid$concentration, acid$MIC))

    coefMin = acid$getCoefMin()
    coefMax = acid$getCoefMax()
    coefMid = acid$getCoefMid()

    coefficients = MGObj$coefficients
    coefficients$mu = coefficients$mu * coefMid
    if (!is.linear(e1)) {coefficients$lambda = coefficients$lambda * (1 / coefMid)}
    #coefficients$lambda = coefficients$lambda * (1 / coefMid)

    conf = MGObj$f$confint()
    conf["mu",] = c(conf["mu",1] * coefMin,
                    conf["mu",2] * coefMax)
    if (!is.linear(e1)){
    conf["lambda",] = c(conf["lambda",1] * (1 / coefMax),
                        conf["lambda",2] * (1 / coefMin))
    }

    MGObj.new = MicrobialGrowth.create(
      c(coefficients$N0, unname(conf["N0",])),
      c(coefficients$Nmax, unname(conf["Nmax",])),
      c(coefficients$mu, unname(conf["mu",])),
      c(coefficients$lambda, unname(conf["lambda",])),
      xlim = c(min(MGObj$data$x), max(MGObj$data$x)),
      model = getModelName(MGObj),
      n = length(MGObj$data$x)
    )

    return( MGObj.new )
  }
  if (is.acid(e1) & is.MicrobialGrowth(e2)) {
    stop("MicrobialGrowth must be on left side")
  }

  stop(gettextf("No operation defined between \"%s\" and \"%s\"", class(e1)[1], class(e2)[1]))
}