--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/w4mcorcov_calc.R	Mon Oct 16 14:56:52 2017 -0400
@@ -0,0 +1,545 @@
+# center with 'colMeans()' - ref: http://gastonsanchez.com/visually-enforced/how-to/2014/01/15/Center-data-in-R/
+center_colmeans <- function(x) {
+  xcenter = colMeans(x)
+  x - rep(xcenter, rep.int(nrow(x), ncol(x)))
+}
+
+#### OPLS-DA
+algoC <- "nipals"
+
+do_detail_plot <- function(x_dataMatrix, x_predictor, x_is_match, x_algorithm, x_prefix, x_show_labels, x_progress = print, x_env) {
+  off <- function(x) if (x_show_labels) x else 0
+  salience_lookup <- x_env$salience_lookup
+  salient_rcv_lookup <- x_env$salient_rcv_lookup
+  # x_progress("head(salience_df): ", head(salience_df))
+  # x_progress("head(salience): ", head(salience))
+  if (x_is_match && ncol(x_dataMatrix) > 0 && length(unique(x_predictor))> 1) {
+    my_oplsda <- opls(
+        x      = x_dataMatrix
+      , y      = x_predictor
+      , algoC  = x_algorithm
+      , predI  = 1
+      , orthoI = if (ncol(x_dataMatrix) > 1) 1 else 0
+      , printL = FALSE
+      , plotL  = FALSE
+      )
+    my_oplsda_suppLs_y_levels <- levels(as.factor(my_oplsda@suppLs$y))
+    fctr_lvl_1 <- my_oplsda_suppLs_y_levels[1]
+    fctr_lvl_2 <- my_oplsda_suppLs_y_levels[2]
+    my_cor_vs_cov <- cor_vs_cov(
+        matrix_x = x_dataMatrix
+      , ropls_x  = my_oplsda
+      )
+    with(
+      my_cor_vs_cov
+    , {
+        min_x <- min(covariance)
+        max_x <- max(covariance)
+        lim_x <- max(sapply(X=c(min_x, max_x), FUN=abs))
+        covariance <- covariance / lim_x
+        lim_x <- 1.2
+        main_label <- sprintf("%s for levels %s versus %s", x_prefix, fctr_lvl_1, fctr_lvl_2)
+        # print("main_label")
+        # print(main_label)
+        main_cex <- min(1.0, 46.0/nchar(main_label))
+        # "It is generally accepted that a variable should be selected if vj>1, [27–29],
+        #   but a proper threshold between 0.83 and 1.21 can yield more relevant variables according to [28]."
+        #   (Mehmood 2012 doi:10.1186/1748-7188-6-27)
+        vipco <- pmax(0, pmin(1,(vip4p-0.83)/(1.21-0.83)))
+        alpha <- 0.1 + 0.4 * vipco
+        red  <- as.numeric(correlation < 0) * vipco
+        blue <- as.numeric(correlation > 0) * vipco
+        minus_cor <- -correlation
+        minus_cov <- -covariance
+        # cex <- salience_lookup(feature = names(minus_cor))
+        # cex <- 0.25 + (1.25 * cex / max(cex))
+        cex <- 0.75
+        plot(
+          y = minus_cor
+        , x = minus_cov
+        , type="p"
+        , xlim=c(-lim_x, lim_x + off(0.1))
+        , ylim=c(-1.0 - off(0.1), 1.0)
+        , xlab = sprintf("relative covariance(feature,t1)")
+        , ylab = sprintf("correlation(feature,t1)")
+        , main = main_label
+        , cex.main = main_cex
+        , cex = cex
+        , pch = 16
+        , col = rgb(blue = blue, red = red, green = 0, alpha = alpha)
+        )
+        low_x <- -0.7 * lim_x
+        high_x <- 0.7 * lim_x
+        text(x = low_x, y = -0.15, labels =  fctr_lvl_1)
+        text(x = high_x, y = 0.15, labels =  fctr_lvl_2)
+        if (x_show_labels) {
+          text(
+            y = minus_cor - 0.013
+          , x = minus_cov + 0.020
+          , cex = 0.3
+          , labels = tsv1$featureID
+          , col = rgb(blue = blue, red = red, green = 0, alpha = 0.2 + 0.8 * alpha)
+          , srt = -30 # slant 30 degrees downward
+          , adj = 0   # left-justified
+          )
+        }
+      }
+    )
+    typeVc <- c("correlation",      # 1
+                "outlier",          # 2
+                "overview",         # 3
+                "permutation",      # 4
+                "predict-train",    # 5
+                "predict-test",     # 6
+                "summary",          # 7 = c(2,3,4,9)
+                "x-loading",        # 8
+                "x-score",          # 9
+                "x-variance",       # 10
+                "xy-score",         # 11
+                "xy-weight"         # 12
+               )                    # [c(3,8,9)] # [c(4,3,8,9)]
+    if ( length(my_oplsda@orthoVipVn) > 0 ) {
+      my_typevc <- typeVc[c(9,3,8)]
+    } else {
+      my_typevc <- c("(dummy)","overview","(dummy)")
+    }
+    for (my_type in my_typevc) {
+      if (my_type %in% typeVc) {
+        # print(sprintf("plotting type %s", my_type))
+        plot(
+          x            = my_oplsda
+        , typeVc       = my_type
+        , parCexN      = 0.4
+        , parDevNewL   = FALSE
+        , parLayL      = TRUE
+        , parEllipsesL = TRUE
+        )
+      } else {
+        # print("plotting dummy graph")
+        plot(x=1, y=1, xaxt="n", yaxt="n", xlab="", ylab="", type="n")
+        text(x=1, y=1, labels="no orthogonal projection is possible")
+      }
+    }
+    return (my_cor_vs_cov)
+  } else {
+    # x_progress(sprintf("x_is_match = %s, ncol(x_dataMatrix) = %d, length(unique(x_predictor)) = %d",x_is_match, ncol(x_dataMatrix), length(unique(x_predictor))))
+    return (NULL)
+  }
+}
+
+# S-PLOT and OPLS reference: Wiklund_2008 doi:10.1021/ac0713510
+corcov_calc <- function(calc_env, failure_action = stop, progress_action = function(x){}, corcov_tsv_action = function(t){}, salience_tsv_action = function(t){}) {
+  if ( ! is.environment(calc_env) ) {
+    failure_action("corcov_calc: fatal error - 'calc_env' is not an environment")
+    return ( FALSE )
+  }
+  if ( ! is.function(corcov_tsv_action) ) {
+    failure_action("corcov_calc: fatal error - 'corcov_tsv_action' is not a function")
+    return ( FALSE )
+  }
+  if ( ! is.function(salience_tsv_action) ) {
+    failure_action("salience_calc: fatal error - 'salience_tsv_action' is not a function")
+    return ( FALSE )
+  }
+
+  # extract parameters from the environment
+  vrbl_metadata <- calc_env$vrbl_metadata
+  vrbl_metadata_names <- vrbl_metadata[,1]
+  smpl_metadata <- calc_env$smpl_metadata
+  data_matrix <- calc_env$data_matrix
+  pairSigFeatOnly <- calc_env$pairSigFeatOnly
+  facC <- calc_env$facC
+  tesC <- calc_env$tesC
+  # extract the levels from the environment
+  originalLevCSV <- levCSV <- calc_env$levCSV
+  # matchingC is one of { "none", "wildcard", "regex" }
+  matchingC <- calc_env$matchingC
+  labelFeatures <- calc_env$labelFeatures
+
+  # arg/env checking
+  if (!(facC %in% names(smpl_metadata))) {
+    failure_action(sprintf("bad parameter!  Factor name '%s' not found in sampleMetadata", facC))
+    return ( FALSE )
+  }
+
+  # calculate salience_df as data.frame(feature, max_level, max_median, max_rcv, mean_median, salience, salient_rcv)
+  salience_df <- calc_env$salience_df <- w4msalience(
+    data_matrix    = data_matrix
+  , sample_class   = smpl_metadata[,facC]
+  , failure_action = failure_action
+  )
+  salience_tsv_action({
+    my_df <- data.frame(
+      featureID    = salience_df$feature
+    , salientLevel = salience_df$max_level
+    , salientRCV   = salience_df$salient_rcv
+    , salience     = salience_df$salience
+    )
+    my_df[order(-my_df$salience),]
+  })
+  salience             <- salience_df$salience
+  names(salience)      <- salience_df$feature
+  salience_lookup      <- calc_env$salience_lookup <- function(feature) unname(salience[feature])
+  salient_rcv          <- salience_df$salient_rcv
+  names(salient_rcv)   <- salience_df$feature
+  salient_rcv_lookup   <- calc_env$salient_rcv_lookup <- function(feature) unname(salient_rcv[feature])
+  salient_level        <- salience_df$max_level
+  names(salient_level) <- salience_df$feature
+  salient_level_lookup <- calc_env$salient_level_lookup <- function(feature) unname(salient_level[feature])
+  
+  # transform wildcards to regexen
+  if (matchingC == "wildcard") {
+    # strsplit(x = "hello,wild,world", split = ",")
+    levCSV <- gsub("[.]", "[.]", levCSV)
+    levCSV <- strsplit(x = levCSV, split = ",")
+    levCSV <- sapply(levCSV, utils::glob2rx, trim.tail = FALSE)
+    levCSV <- paste(levCSV, collapse = ",")
+  }
+  # function to determine whether level is a member of levCSV
+  isLevelSelected <- function(lvl) {
+    matchFun <- if (matchingC != "none") grepl else `==`
+    return(
+      Reduce(
+        f = "||"
+      , x = sapply(
+              X = strsplit(
+                    x = levCSV
+                  , split = ","
+                  , fixed = TRUE
+                  )[[1]]
+            , FUN = matchFun
+            , lvl
+            )
+      )
+    )
+  }
+
+  # transpose matrix because ropls matrix is the transpose of XCMS matrix
+  # Wiklund_2008 centers and pareto-scales data before OPLS-DA S-plot
+  # center
+  cdm <- center_colmeans(t(data_matrix))
+  # pareto-scale
+  my_scale <- sqrt(apply(cdm, 2, sd, na.rm=TRUE))
+  scdm <- sweep(cdm, 2, my_scale, "/")
+
+  # pattern to match column names like k10_kruskal_k4.k3_sig
+  col_pattern <- sprintf('^%s_%s_(.*)[.](.*)_sig$', facC, tesC)
+  # column name like k10_kruskal_sig
+  intersample_sig_col  <- sprintf('%s_%s_sig', facC, tesC)
+  # get the facC column from sampleMetadata, dropping to one dimension
+  smpl_metadata_facC <- smpl_metadata[,facC]
+
+  # allocate a slot in the environment for the contrast_list, each element of which will be a data.frame with this structure:
+  #   - feature ID
+  #   - value1
+  #   - value2
+  #   - Wiklund_2008 correlation
+  #   - Wiklund_2008 covariance
+  #   - Wiklund_2008 VIP
+  calc_env$contrast_list <- list()
+
+
+  did_plot <- FALSE
+  if (tesC != "none") {
+    # for each column name, extract the parts of the name matched by 'col_pattern', if any
+    the_colnames <- colnames(vrbl_metadata)
+    if ( ! Reduce( f = "||", x = grepl(tesC, the_colnames) ) ) {
+      failure_action(sprintf("bad parameter!  variableMetadata must contain results of W4M Univariate test '%s'.", tesC))
+      return ( FALSE )
+    }
+    col_matches <- lapply(
+      X = the_colnames,
+      FUN = function(x) {
+        regmatches( x, regexec(col_pattern, x) )[[1]]
+      }
+    )
+    ## first contrast each selected level with all other levels combined into one "super-level" ##
+    # process columns matching the pattern
+    level_union <- c()
+    for (i in 1:length(col_matches)) {
+      col_match <- col_matches[[i]]
+      if (length(col_match) > 0) {
+        # it's an actual match; extract the pieces, e.g., k10_kruskal_k4.k3_sig
+        vrbl_metadata_col <- col_match[1]               # ^^^^^^^^^^^^^^^^^^^^^  # Column name
+        fctr_lvl_1        <- col_match[2]               #             ^^         # Factor-level 1
+        fctr_lvl_2        <- col_match[3]               #                ^^      # Factor-level 2
+        # only process this column if both factors are members of lvlCSV
+        is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2)
+        if (is_match) {
+          level_union <- c(level_union, col_match[2], col_match[3])
+        }
+      }
+    }
+    level_union <- unique(sort(level_union))
+    overall_significant <- 1 == vrbl_metadata[,intersample_sig_col]
+    if ( length(level_union) > 1 ) {
+      chosen_samples <- smpl_metadata_facC %in% level_union
+      chosen_facC <- as.character(smpl_metadata_facC[chosen_samples])
+      col_selector <- vrbl_metadata_names[ overall_significant ]
+      my_matrix <- scdm[ chosen_samples, col_selector, drop = FALSE ]
+      fctr_lvl_2 <- "other"
+      for ( fctr_lvl_1 in level_union[1:length(level_union)] ) {
+        progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2))
+        predictor <- sapply( X = chosen_facC, FUN = function(fac) if ( fac == fctr_lvl_1 ) fctr_lvl_1 else fctr_lvl_2 )
+        my_cor_cov <- do_detail_plot(
+          x_dataMatrix  = my_matrix
+        , x_predictor   = predictor
+        , x_is_match    = is_match
+        , x_algorithm   = algoC
+        , x_prefix      = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features"
+        , x_show_labels = labelFeatures
+        , x_progress    = progress_action
+        , x_env         = calc_env
+        )
+        if ( is.null(my_cor_cov) ) {
+          progress_action("NOTHING TO PLOT.")
+        } else {
+          tsv <- my_cor_cov$tsv1
+          # tsv$salientLevel <- salient_level_lookup(tsv$featureID)
+          # tsv$salientRCV   <- salient_rcv_lookup(tsv$featureID)
+          # tsv$salience     <- salience_lookup(tsv$featureID)
+          tsv["level1Level2Sig"] <- vrbl_metadata[ match(tsv$featureID, vrbl_metadata_names), vrbl_metadata_col ] 
+          corcov_tsv_action(tsv)
+          did_plot <- TRUE
+        }
+      }
+    }
+
+    ## next, contrast each selected level with each of the other levels individually ##
+    # process columns matching the pattern
+    for (i in 1:length(col_matches)) {
+      # for each potential match of the pattern
+      col_match <- col_matches[[i]]
+      if (length(col_match) > 0) {
+        # it's an actual match; extract the pieces, e.g., k10_kruskal_k4.k3_sig
+        vrbl_metadata_col <- col_match[1]               # ^^^^^^^^^^^^^^^^^^^^^  # Column name
+        fctr_lvl_1        <- col_match[2]               #             ^^         # Factor-level 1
+        fctr_lvl_2        <- col_match[3]               #                ^^      # Factor-level 2
+        # only process this column if both factors are members of lvlCSV
+        is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2)
+        progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2))
+        # TODO delete next line displaying currently-processed column
+        # cat(sprintf("%s %s %s %s\n", vrbl_metadata_col, fctr_lvl_1, fctr_lvl_2, is_match))
+        # choose only samples with one of the two factors for this column
+        chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2)
+        predictor <- smpl_metadata_facC[chosen_samples]
+        # extract only the significantly-varying features and the chosen samples
+        fully_significant   <- 1 == vrbl_metadata[,vrbl_metadata_col] * vrbl_metadata[,intersample_sig_col]
+        col_selector <- vrbl_metadata_names[ if ( pairSigFeatOnly ) fully_significant else overall_significant ]
+        my_matrix <- scdm[ chosen_samples, col_selector, drop = FALSE ]
+        my_cor_cov <- do_detail_plot(
+          x_dataMatrix  = my_matrix
+        , x_predictor   = predictor
+        , x_is_match    = is_match
+        , x_algorithm   = algoC
+        , x_prefix      = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features"
+        , x_show_labels = labelFeatures
+        , x_progress    = progress_action
+        , x_env         = calc_env
+        )
+        if ( is.null(my_cor_cov) ) {
+          progress_action("NOTHING TO PLOT.")
+        } else {
+          tsv <- my_cor_cov$tsv1
+          # tsv$salientLevel <- salient_level_lookup(tsv$featureID)
+          # tsv$salientRCV   <- salient_rcv_lookup(tsv$featureID)
+          # tsv$salience     <- salience_lookup(tsv$featureID)
+          tsv["level1Level2Sig"] <- vrbl_metadata[ match(tsv$featureID, vrbl_metadata_names), vrbl_metadata_col ] 
+          corcov_tsv_action(tsv)
+          did_plot <- TRUE
+        }
+      }
+    }
+  } else { # tesC == "none"
+    level_union <- unique(sort(smpl_metadata_facC))
+    if ( length(level_union) > 1 ) {
+      ## pass 1 - contrast each selected level with all other levels combined into one "super-level" ##
+      completed <- c()
+      lapply(
+        X = level_union
+      , FUN = function(x) { 
+          fctr_lvl_1 <- x[1]
+          fctr_lvl_2 <- {
+            if ( fctr_lvl_1 %in% completed )
+              return("DUMMY")
+            # strF(completed)
+            completed <<- c(completed, fctr_lvl_1)
+            setdiff(level_union, fctr_lvl_1)
+          }
+          chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2)
+          fctr_lvl_2 <- "other"
+          # print( sprintf("sum(chosen_samples) %d, factor_level_2 %s", sum(chosen_samples), fctr_lvl_2) )
+          progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2))
+          if (length(unique(chosen_samples)) < 1) {
+            progress_action("NOTHING TO PLOT...")
+          } else {
+            chosen_facC <- as.character(smpl_metadata_facC[chosen_samples])
+            predictor <- sapply( X = chosen_facC, FUN = function(fac) if ( fac == fctr_lvl_1 ) fctr_lvl_1 else fctr_lvl_2 )
+            my_matrix <- scdm[ chosen_samples, , drop = FALSE ]
+            # only process this column if both factors are members of lvlCSV
+            is_match <- isLevelSelected(fctr_lvl_1)
+            my_cor_cov <- do_detail_plot(
+              x_dataMatrix  = my_matrix
+            , x_predictor   = predictor
+            , x_is_match    = is_match
+            , x_algorithm   = algoC
+            , x_prefix      = "Features"
+            , x_show_labels = labelFeatures
+            , x_progress    = progress_action
+            , x_env         = calc_env
+            )
+            if ( is.null(my_cor_cov) ) {
+              progress_action("NOTHING TO PLOT")
+            } else {
+              tsv <- my_cor_cov$tsv1
+              # tsv$salientLevel <- salient_level_lookup(tsv$featureID)
+              # tsv$salientRCV   <- salient_rcv_lookup(tsv$featureID)
+              # tsv$salience     <- salience_lookup(tsv$featureID)
+              corcov_tsv_action(tsv)
+              did_plot <<- TRUE
+            }
+          }
+          #print("baz")
+          "dummy" # need to return a value; otherwise combn fails with an error
+        }
+      )
+      ## pass 2 - contrast each selected level with each of the other levels individually ##
+      completed <- c()
+      utils::combn(
+        x = level_union
+      , m = 2
+      , FUN = function(x) { 
+          fctr_lvl_1 <- x[1]
+          fctr_lvl_2 <- x[2]
+          chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2)
+          # print( sprintf("sum(chosen_samples) %d, factor_level_2 %s", sum(chosen_samples), fctr_lvl_2) )
+          progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2))
+          if (length(unique(chosen_samples)) < 1) {
+            progress_action("NOTHING TO PLOT...")
+          } else {
+            chosen_facC <- as.character(smpl_metadata_facC[chosen_samples])
+            predictor <- chosen_facC
+            my_matrix <- scdm[ chosen_samples, , drop = FALSE ]
+            # only process this column if both factors are members of lvlCSV
+            is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2)
+            my_cor_cov <- do_detail_plot(
+              x_dataMatrix  = my_matrix
+            , x_predictor   = predictor
+            , x_is_match    = is_match
+            , x_algorithm   = algoC
+            , x_prefix      = "Features"
+            , x_show_labels = labelFeatures
+            , x_progress    = progress_action
+            , x_env         = calc_env
+            )
+            if ( is.null(my_cor_cov) ) {
+              progress_action("NOTHING TO PLOT")
+            } else {
+              tsv <- my_cor_cov$tsv1
+              # tsv$salientLevel <- salient_level_lookup(tsv$featureID)
+              # tsv$salientRCV   <- salient_rcv_lookup(tsv$featureID)
+              # tsv$salience     <- salience_lookup(tsv$featureID)
+              corcov_tsv_action(tsv)
+              did_plot <<- TRUE
+            }
+          }
+          #print("baz")
+          "dummy" # need to return a value; otherwise combn fails with an error
+        }
+      )
+    } else {
+      progress_action("NOTHING TO PLOT....")
+    }
+  }
+  if (!did_plot) {
+    failure_action(sprintf("bad parameter!  sampleMetadata must have at least two levels of factor '%s' matching '%s'", facC, originalLevCSV))
+    return ( FALSE )
+  }
+  return ( TRUE )
+}
+
+# Calculate data for correlation-versus-covariance plot
+#   Adapted from:
+#     Wiklund_2008 doi:10.1021/ac0713510
+#     Galindo_Prieto_2014 doi:10.1002/cem.2627
+#     https://github.com/HegemanLab/extra_tools/blob/master/generic_PCA.R
+cor_vs_cov <- function(matrix_x, ropls_x) {
+  x_class <- class(ropls_x)
+  if ( !( as.character(x_class) == "opls" ) ) { # || !( attr(class(x_class),"package") == "ropls" ) ) 
+    stop( "cor_vs_cov: Expected ropls_x to be of class ropls::opls but instead it was of class ", as.character(x_class) )
+  }
+  result <- list()
+  # suppLs$algoC - Character: algorithm used - "svd" for singular value decomposition; "nipals" for NIPALS
+  if ( ropls_x@suppLs$algoC == "nipals") {
+    # Equations (1) and (2) from *Supplement to* Wiklund 2008, doi:10.1021/ac0713510
+    mag <- function(one_dimensional) sqrt(sum(one_dimensional * one_dimensional))
+    mag_xi <- sapply(X = 1:ncol(matrix_x), FUN = function(x) mag(matrix_x[,x]))
+    score_matrix <- ropls_x@scoreMN
+    score_matrix_transposed <- t(score_matrix)
+    score_matrix_magnitude <- mag(score_matrix)
+    result$covariance <- score_matrix_transposed %*% matrix_x / ( score_matrix_magnitude * score_matrix_magnitude )
+    result$correlation <- score_matrix_transposed %*% matrix_x / ( score_matrix_magnitude * mag_xi )
+  } else {
+    # WARNING - untested code - I don't have test data to exercise this branch
+    # Equations (1) and (2) from Wiklund 2008, doi:10.1021/ac0713510
+    # scoreMN - Numerical matrix of x scores (T; dimensions: nrow(x) x predI) X = TP' + E; Y = TC' + F
+    score_matrix <- ropls_x@scoreMN
+    score_matrix_transposed <- t(score_matrix)
+    cov_divisor <- nrow(matrix_x) - 1
+    result$covariance <- sapply(
+      X = 1:ncol(matrix_x)
+    , FUN = function(x) score_matrix_transposed %*% matrix_x[,x] / cov_divisor
+    )
+    score_sd <- sapply(
+      X = 1:ncol(score_matrix)
+      , FUN = function(x) sd(score_matrix[,x])
+    )
+    # xSdVn - Numerical vector: variable standard deviations of the 'x' matrix
+    xSdVn <- ropls_x@xSdVn
+    result$correlation <- sapply(
+      X = 1:ncol(matrix_x)
+    , FUN = function(x) {
+        ( score_matrix_transposed / score_sd ) %*% ( matrix_x[,x] / (xSdVn[x] * cov_divisor) )
+      }
+    )
+  }
+  result$correlation <- result$correlation[1,,drop = TRUE]
+  result$covariance  <- result$covariance[1,,drop = TRUE]
+
+  # Variant 4 of Variable Influence on Projection for OPLS from Galindo_Prieto_2014
+  #    Length = number of features; labels = feature identifiers.  (The same is true for $correlation and $covariance.)
+  result$vip4p     <- as.numeric(ropls_x@vipVn)
+  result$vip4o     <- as.numeric(ropls_x@orthoVipVn)
+  # get the level names
+  level_names      <- sort(levels(as.factor(ropls_x@suppLs$y)))
+  feature_count    <- length(ropls_x@vipVn)
+  result$level1    <- rep.int(x = level_names[1], times = feature_count)
+  result$level2    <- rep.int(x = level_names[2], times = feature_count)
+  # print(sprintf("sd(covariance) = %f; sd(correlation) = %f", sd(result$covariance), sd(result$correlation)))
+  superresult <- list()
+  if (length(result$vip4o) == 0) result$vip4o <- NA
+  superresult$tsv1 <- data.frame(
+    featureID           = names(ropls_x@vipVn)
+  , factorLevel1        = result$level1
+  , factorLevel2        = result$level2
+  , correlation         = result$correlation
+  , covariance          = result$covariance
+  , vip4p               = result$vip4p
+  , vip4o               = result$vip4o
+  , row.names           = NULL
+  )
+  rownames(superresult$tsv1) <- superresult$tsv1$featureID
+  superresult$covariance <- result$covariance
+  superresult$correlation <- result$correlation
+  superresult$vip4p <- result$vip4p
+  superresult$vip4o <- result$vip4o
+  superresult$details <- result
+  # #print(superresult$tsv1)
+  result$superresult <- superresult
+  # Include thise in case future consumers of this routine want to use it in currently unanticipated ways
+  result$oplsda    <- ropls_x          
+  result$predictor <- ropls_x@suppLs$y   # in case future consumers of this routine want to use it in currently unanticipated ways
+  return (superresult)
+}
+
+