R/utils.R

In ELAplus: Energy Landscape Analysis for Ecological Communities

#' applymap
#' @description Apply a function to each element by row and return a matrix-like result.
#' @param df A matrix/data.frame to transform.
#' @param fun A function applied to each element.
#' @return A matrix (as returned by t(apply(...))) of transformed values.
#' @noRd
applymap <- function(df, fun){
  return(t(apply(df, 1, function(row){sapply(row, fun)})))
}

#' applyvec
#' @description Apply a function to each column (or vector) and preserve row/column names.
#' @param df A data.frame (or matrix coercible) to transform.
#' @param fun A function applied to each column (via lapply).
#' @return A data.frame with transformed columns and preserved names.
#' @noRd
applyvec <- function(df, fun){
  if(length(rownames(df)) > 1){
    output <- t(data.frame(lapply(t(df), fun))) %>% 'rownames<-'(rownames(df))
  }else(
    output <- data.frame(lapply(df, fun)) %>% 'colnames<-'(colnames(df))
  )
  return(output)
}

#' T
#' @description Transpose a data.frame and return it as a data.frame.
#' @param df A matrix/data.frame to transpose.
#' @return A transposed data.frame.
#' @noRd
T <- function(df){return(as.data.frame(t(df)))}

#' Formatting
#' @description Generate occurrence(binary)/abundance/factor matrices and labels from input abundance/metadata tables.
#' @param baseabtable Abundance table (samples x species) as data.frame/matrix.
#' @param basemetadata Optional environmental metadata table (samples x factors); NULL/empty for none.
#' @param normalize Integer 0 or 1. If 1, normalize within each sample to sum to 1.
#' @param grouping Integer 0 or 1. If 1, group similar occurrence patterns.
#' @param grouping_th Numeric threshold for grouping (relative Hamming distance).
#' @param reporting Logical. If TRUE, print processing summary.
#' @param parameters Numeric vector/list of thresholds (following 3 parameters).
#'   1. "0/1 threshold"
#'   2. "min occurrence threshold"
#'   3. "max occurrence threshold"
#'   (e.g., c(0.05, 0.05, 0.95) : species abundance > 0.05 is regarded as present,
#'   and species occurring in more than 0.05 and less than 0.95 in the all samples are kept.)
#' @return A list: (ocmatrix, abmatrix, fctmatrix, samplelabel, specieslabel, factorlabel).
#' @useDynLib ELAplus, .registration=TRUE
#' @export
#' @importFrom foreach foreach %do%
#' @importFrom stringdist stringdist

Formatting <- function(baseabtable, basemetadata=NULL, normalize=1, 
                       parameters=c(0.05,0.05,0.95), grouping=0, grouping_th=0, reporting=TRUE){
  ath <- parameters[1]     # Relative frequency threshold; if less than this value, set to 0
  minocth <- parameters[2] # Lower limit of mean occurence; species with mean occurence lower than this value are excluded.
  maxocth <- parameters[3] # Upper limit of mean occurence; species with mean occurence higher than this value are excluded.
  
  if(length(basemetadata) != 0){
    # with basemetadata
    # Get label information
    specieslabel = colnames(baseabtable)
    absamplelabel = rownames(baseabtable)
    mdsamplelabel = rownames(basemetadata)
    sharedsamplelabel = sort(intersect(absamplelabel, mdsamplelabel))

    # Generating occurence matrix
    patab = match(sharedsamplelabel, absamplelabel)
    patmd = match(sharedsamplelabel, mdsamplelabel)
    abnum = baseabtable[patab,]
    if (normalize == 1){ # Nomalization within samples (optional)
      abnum <- t(apply(abnum, 1, function(x){x / (sum(x) + 10^-16)}))
    }
    aboc = applymap(abnum, function(x){if(x >= ath){1}else{0}}) # 0/1 transformation with ath
    aboc_mean = apply(aboc, 2, mean)
    occrit = names(aboc_mean[(aboc_mean > minocth) & (aboc_mean < maxocth)]) # Calculation of mean occurence//determination of the target to be removed

    # Object generation
    ocmatrix = aboc[, occrit] # Occurence matrix
    abmatrix = abnum[, occrit] # Abundance matrix
    fctmatrix = as.matrix(as.data.frame(basemetadata[patmd, ])) # Environmental factors
    rownames(fctmatrix) <- rownames(ocmatrix) 
    colnames(fctmatrix) <- colnames(basemetadata)
    samplelabel = sharedsamplelabel # Sample labels
    specieslabel = occrit # Species labels
    factorlabel = colnames(basemetadata) # Factor labels
  }else{
    # without basemetadata
    # Get label information
    specieslabel = colnames(baseabtable)
    absamplelabel = rownames(baseabtable)
    sharedsamplelabel = sort(absamplelabel)

    # Generating occurence matrix
    patab = match(sharedsamplelabel, absamplelabel)
    abnum = baseabtable[patab,]
    if (normalize == 1){ # Nomalization within samples (optional)
      abnum <- t(apply(abnum, 1, function(x){x / (sum(x) + 10^-16)}))
    }
    aboc = applymap(abnum, function(x){if(x >= ath){1}else{0}}) # 0/1 transformation with ath
    aboc_mean = apply(aboc, 2, mean)
    occrit = names(aboc_mean[(aboc_mean > minocth) & (aboc_mean < maxocth)]) # Calculation of mean occurence//determination of the target to be removed

    # Object generation
    ocmatrix = aboc[, occrit] # Occurence matrix
    abmatrix = abnum[, occrit] # Abundance matrix
    fctmatrix = NULL # Environmental factors
    samplelabel = sharedsamplelabel # Sample labels
    specieslabel = occrit # Species labels
    factorlabel = NULL # Factor labels
  }

  # Grouping
  if(grouping == 1){
    ocvstr <- apply(ocmatrix, 2, function(y){return(paste(y, collapse = ""))})
    oo <- order(ocvstr)
    ov <- ocvstr[oo]
    j <- 0
    cc <- 0
    binlist <- list()
    aaa <- foreach (i=ov) %do% {
        cc <- cc + 1
        if(stringdist(i, j, method='hamming')/nchar(i) <= grouping_th){a <- rev(binlist)[[1]]
                aa <- append(a, cc)
                binlist <- append(binlist[-length(binlist)], list(aa))
                }else{binlist <- append(binlist, list(cc))}
        j=i
        }
    ltx <- length(specieslabel)
    specieslabel <- foreach (i=binlist, .combine='c') %do% {paste(specieslabel[oo][i], collapse=", ")}
    ocmatrix <- foreach(i=binlist, .combine='cbind') %do% {ocmatrix[,oo][,i[1]]}
    abmatrix <- foreach(i=binlist, .combine='cbind') %do% {if(length(i)>1){apply(abmatrix[,oo][,i],1,sum)}else{abmatrix[,oo][,i]}}
    colnames(ocmatrix) <- specieslabel
    colnames(abmatrix) <- specieslabel
  }

  if(reporting){
    cat('Processed', as.character(length(samplelabel)), 'samples.\n')
    cat('Relative abundance threshold =', as.character(ath), '\n')
    cat('Occurrence threshold (lower) =', as.character(minocth), '\n')
    cat('Occurrence threshold (upper) =', as.character(maxocth), '\n')
    if(grouping == 1){cat(as.character(ltx-length(specieslabel)), ' groups were found.\n')}
    cat('Selected ', as.character(length(specieslabel)), ' out of ',
                as.character(ncol(baseabtable)), 'species.\n')
  }
  return(list(ocmatrix, abmatrix, fctmatrix, samplelabel, specieslabel, factorlabel))}


#' str2i
#' @description Convert a binary string (e.g., "0101") to an integer.
#' @param x A character string consisting of 0/1 digits.
#' @return A numeric scalar giving the decoded integer value.
#' @noRd
str2i <- function(x) {
    y <- as.numeric(strsplit(x, "")[[1]])
    sum(y * 2^rev((seq_along(y)-1)))
}

#'bin2id: map a binary vector to a 64base characters
#' @description Encode a binary vector into a base-64 id string index.
#' @param bin A binary vector (0/1, community composition) to encode.
#' @return A character scalar representing the encoded id.
#' @export
#' @importFrom foreach foreach
#' @export
bin2id <- function(bin){
  digs <- "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz_/"
  bi <- unlist(bin)
  modi <- length(bi)%%6
  if(modi != 0){b0 <- rep(0, 6-modi)}else{b0<-c()}
  bii <- c(b0, bi)
  lb <- length(bii)
  aa <- foreach(i=seq(1, length(bii), 6), .combine="c") %do% {a <- rev(bii)[i:(i+5)]
                                s <- str2i(paste(as.character(rev(a)), collapse=""))
                                substr(digs, s+1, s+1)}
  return(paste(rev(aa), collapse=""))}

#'id2bin: map a 64base characters to a binary vector
#' @description Decode a base-64 id string index to a binary vector of length s.
#' @param id A character scalar encoded by \link{bin2id}.
#' @param s Integer number of bits to return.
#' @return An integer vector of 0/1 bits of length s.
#' @export
#' @importFrom foreach foreach
#' @export
id2bin <-  function(id, s){
    digs <- "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz_/"
    aa <- foreach (i=seq(1:nchar(id)), .combine="rbind") %do% {
        a <- substr(id,i,i)
        ia <- (unlist(gregexpr(a, digs))[1]) - 1
        CIntegerDigits(ia, 6)}
    uaa <- unlist(aa)
    aaa <- as.vector(uaa[(length(uaa)-s+1):length(uaa)])
    return(aaa)}

#' CIntegerDigits
#' @description Map an integer or encoded ID to a binary vector of length n (optionally with row names).
#' @param ssid Numeric state id or encoded id (character) to convert.
#' @param n Integer number of bits.
#' @param index Optional row names (default 1:n).
#' @return A data.frame of bits (0/1) with row names \code{index}.
#' @useDynLib ELAplus, .registration=TRUE
#' @noRd
CIntegerDigits <- function(ssid, n, index=NULL){
  if(is.null(index)){index <- 1:n}
  if(is.numeric(ssid)){
    return(data.frame(dec2bin(ssid, n), row.names=index))
  }
  else{
    return(data.frame(id2bin(ssid, n), row.names=index))
  }
}

#' dec2bin
#' @description Convert a non-negative integer to a binary vector of fixed length.
#' @param num Integer to convert.
#' @param digit Integer length of the output binary vector.
#' @return An integer vector of 0/1 bits (or NULL if both num<=0 and digit<=0).
#' @noRd
dec2bin <- function(num, digit=0){
  if(num <= 0 && digit <= 0){
    return(NULL)
  }else{
    return(append(Recall(num%/%2,digit-1), num%%2))
  }
}