w4msalience <- function(
  data_matrix   # a matrix of intensities; features as rows, and samples as columns
, sample_class  # a vector of sample class-levels; length(sample_class) == ncol(data_matrix)
, failure_action = stop
) {
  library(stats)
  # begin sanity checks
  if ( !is.vector(sample_class) || !( is.character(sample_class) || is.factor(sample_class) ) ) {
    failure_action("w4msalience:  Expected sample_class to be a vector of characters of factor-levels")
    return (NULL)
  }
  if ( !is.matrix(data_matrix) && !is.data.frame(data_matrix) ) {
    failure_action("w4msalience:  Expected data_matrix to be a matrix (or data.frame) of numeric")
    return (NULL)
  }
  # transpose data_matrix so that columns are the features
  t_data_matrix <- t(data_matrix)
  if ( !is.matrix(t_data_matrix) || !is.numeric(t_data_matrix) ) {
    failure_action("w4msalience:  Expected data_matrix to be a matrix (or data.frame) of numeric")
    return (NULL)
  }

  feature_names <- colnames(t_data_matrix)

  n_features <- ncol(t_data_matrix)
  n_samples  <- nrow(t_data_matrix)
  if ( length(sample_class) != n_samples ) {
    strF(data_matrix)
    strF(sample_class)
    failure_action(
      sprintf(
        "w4msalience:  The data_matrix has %d samples but sample_class has %d"
      , n_samples
      , length(sample_class)
      )
    )
    return (NULL)
  }
  # end sanity checks

  # "For each feature, 'select sample_class, median(intensity) from feature group by sample_class'."
  # The first column(s) of the result of aggregate has the classifier value(s) specified in the 'by' list.
  medianOfFeatureBySampleClassLevel <- aggregate(
      x = as.data.frame(t_data_matrix)
    , by = list(sample_class)
    , FUN = "median"
    )

  # "For each feature, 'select sample_class, rcv(intensity) from feature group by sample_class'."
  #   cv is less robust; deviation from normality degrades performance
  #     cv(x) == sd(x) / mean(x)
  #   rcv is a "robust" coefficient of variation, expressed as a proportion
  #     rcv(x) == mad(x) / median(x)
  madOfFeatureBySampleClassLevel <- aggregate(
      x = as.data.frame(t_data_matrix)
    , by = list(sample_class)
    , FUN = "mad"
  )

  # Note that `apply(X=array(1:10), MARGIN = 1, FUN = function(x) return(c(x,x^2)))`
  #   produces a matrix with two rows and ten columns

  my_list <- apply(
    X = array(1:n_features)
  , MARGIN = 1
  , FUN = function(x) {
      my_df <- data.frame(
        median = medianOfFeatureBySampleClassLevel[ , 1 + x]
      , mad = madOfFeatureBySampleClassLevel[ , 1 + x]
      )
      my_df$salient_level <- medianOfFeatureBySampleClassLevel[ , 1]
      my_df <- my_df[ order(my_df$median, decreasing = TRUE), ]
      my_dist_df <- my_df[  1:2, ]
      # "robust coefficient of variation", i.e.,
      #    mad(feature-intensity for class-level max_level) / median(feature-intensity for class-level max_level)
      rcv_result <- my_dist_df$mad[1] / my_dist_df$median[1]
      dist_result <-
        ( my_dist_df$median[1] - my_dist_df$median[2] ) /
        sqrt( my_dist_df$mad[1] * my_dist_df$mad[2] )
      if (is.infinite(dist_result) || is.nan(dist_result))
        dist_result <- 0
      mean_median <- mean(my_df$median)
      salience_result <- if (mean_median > 0) my_df$median[1] / mean_median else 0
      return (
        data.frame(
          dist_result     = dist_result
        , max_median      = my_df$median[1]
        , mean_median     = mean_median
        , salience_result = salience_result
        , salient_level   = my_df$salient_level[1]
        , rcv_result      = rcv_result
        )
      )
    }
  )
  results_matrix  <- sapply(X = 1:n_features, FUN = function(i) my_list[[i]])
  results_df <- as.data.frame(t(results_matrix))

  relative_salient_distance <- unlist(results_df$dist_result)
  salience <- unlist(results_df$salience_result)
  salient_level <- unlist(results_df$salient_level)
  max_median <- unlist(results_df$max_median)
  mean_median <- unlist(results_df$mean_median)
  rcv_result <- unlist(results_df$rcv_result)

  salience_df <-
    data.frame(
      # the feature name
      feature                   = feature_names
      # the name (or factor-level) of the class-level with the highest median intensity for the feature
    , max_level                 = salient_level
      # the median intensity for the feature and the level max_level
    , max_median                = max_median
      # the distance between the maximum intensities for the feature at the two highest levels
    , relative_salient_distance = relative_salient_distance
      # the coefficient of variation (expressed as a proportion) for the intensity for the feature and the level max_level
    , salience_rcv              = rcv_result
      # the mean of the medians of intensity for all class-levels for the feature
    , mean_median               = mean_median
      # raw salience is the ratio of the most-prominent level to the mean of all levels for the feature
    , salience                  = salience
      # don't coerce strings to factors (this is a parameter for the data.frame constructor, not a column of the data.frame)
    , stringsAsFactors = FALSE
    )

  return (salience_df)
}