@@ -109,10 +109,14 @@ compute.gene.cdf <- function(expr, sample.idxs, Gaussk=TRUE, kernel=TRUE,

             else

                 gene.cdf <- .ecdfvals_sparse_to_dense(expr[, sample.idxs, drop=FALSE],

                                                       verbose)

-        } else if (is.matrix(expr))

-            gene.cdf <- .ecdfvals_dense_to_dense(expr[, sample.idxs, drop=FALSE],

+        } else if (is.matrix(expr)) {

+            if (any_na)

+                gene.cdf <- .ecdfvals_dense_to_dense_nas(expr[, sample.idxs, drop=FALSE],

+                                                         verbose)

+            else

+                gene.cdf <- .ecdfvals_dense_to_dense(expr[, sample.idxs, drop=FALSE],

                                                  verbose)

-        else

+        } else

             stop(sprintf("Matrix class %s cannot be handled yet.", class(expr)))

     }

@@ -973,22 +977,20 @@ setMethod("gsvaEnrichment", signature(param="gsvaRanksParam"),

                   walkStat <- .gsvaRndWalk_nas(gSetIdx, decOrdStat, symRnkStat,

                                                tau, na_use, minSize, wna_env)

                   maxDev <- c(NA_real_, NA_real_)

-                  if (na_use == "na.rm")

-                      maxDev <- c(max(c(0, max(walkStat, na.rm=TRUE))),

-                                  min(c(0, min(walkStat, na.rm=TRUE))))

-                  else

-                      maxDev <- c(max(c(0, max(walkStat))),

-                                  min(c(0, min(walkStat))))

+                  if (any(!is.na(walkStat))) {

+                      if (na_use == "na.rm")

+                          maxDev <- c(max(c(0, max(walkStat, na.rm=TRUE))),

+                                      min(c(0, min(walkStat, na.rm=TRUE))))

+                      else

+                          maxDev <- c(max(c(0, max(walkStat))),

+                                      min(c(0, min(walkStat))))

+                  }

                   maxDev

                 }, rnkstats$dos, rnkstats$srs, na_use)

                 md <- do.call("rbind", md)

                 if (maxDiff && absRanking)

                     md[, 2] <- -1 * md[, 2]

-                sco <- rep(NA_real_, length(geneSetsIdx))

-                if (na_use == "na.rm")

-                    sco <- rowSums(md, na.rm=TRUE)

-                else 

-                    sco <- rowSums(md)

+                sco <- rowSums(md)

                 if (!maxDiff) {

                     mask <- is.na(sco)

                     sco[!mask] <- md[cbind(1:sum(!mask), ifelse(sco[!mask] > 0,

@@ -1026,22 +1028,20 @@ setMethod("gsvaEnrichment", signature(param="gsvaRanksParam"),

                   walkStat <- .gsvaRndWalk_nas(gSetIdx, decOrdStat, symRnkStat,

                                                tau, na_use, minSize, wna_env)

                   maxDev <- c(NA_real_, NA_real_)

-                  if (na_use == "na.rm")

-                      maxDev <- c(max(c(0, max(walkStat, na.rm=TRUE))),

-                                  min(c(0, min(walkStat, na.rm=TRUE))))

-                  else

-                      maxDev <- c(max(c(0, max(walkStat))),

-                                  min(c(0, min(walkStat))))

+                  if (any(!is.na(walkStat))) {

+                      if (na_use == "na.rm")

+                          maxDev <- c(max(c(0, max(walkStat, na.rm=TRUE))),

+                                      min(c(0, min(walkStat, na.rm=TRUE))))

+                      else

+                          maxDev <- c(max(c(0, max(walkStat))),

+                                      min(c(0, min(walkStat))))

+                  }

                   maxDev

                 }, rnkstats$dos, rnkstats$srs, na_use)

                 md <- do.call("rbind", md)

                 if (maxDiff && absRanking)

                     md[, 2] <- -1 * md[, 2]

-                sco <- rep(NA_real_, length(geneSetsIdx))

-                if (na_use == "na.rm")

-                    sco <- rowSums(md, na.rm=TRUE)

-                else 

-                    sco <- rowSums(md)

+                sco <- rowSums(md)

                 if (!maxDiff) {

                     mask <- is.na(sco)

                     sco[!mask] <- md[cbind(1:sum(!mask), ifelse(sco[!mask] > 0,

@@ -1225,6 +1225,12 @@ setMethod("gsvaEnrichment", signature(param="gsvaRanksParam"),

   .Call("ecdfvals_dense_to_dense_R", X, verbose)

 }

+.ecdfvals_dense_to_dense_nas <- function(X, verbose) {

+  stopifnot(is.matrix(X)) ## QC

+  stopifnot(is.logical(verbose)) ## QC

+  .Call("ecdfvals_dense_to_dense_nas_R", X, verbose)

+}

+

 .kcdfvals_sparse_to_sparse <- function(X, Gaussk, verbose) {

   stopifnot(is(X, "CsparseMatrix")) ## QC

   Xrsp <- as(X, "RsparseMatrix")

@@ -32,6 +32,9 @@ ecdfvals_sparse_to_dense_R(SEXP XCspR, SEXP XRspR, SEXP verboseR);

 SEXP

 ecdfvals_dense_to_dense_R(SEXP XR, SEXP verboseR);

+SEXP

+ecdfvals_dense_to_dense_nas_R(SEXP XR, SEXP verboseR);

+

 SEXP

 match_int(SEXP x, SEXP table);

@@ -513,3 +516,134 @@ ecdfvals_dense_to_dense_R(SEXP XR, SEXP verboseR) {

   return(ecdfRobj);

 }

+

+/* calculate empirical cumulative distribution function values

+ * on the zero and nonzero entries from the input dense matrix,

+ * when missing (NA) values are present, the returned value is

+ * a dense matrix.

+ */

+SEXP

+ecdfvals_dense_to_dense_nas_R(SEXP XR, SEXP verboseR) {

+  double* X;

+  Rboolean verbose=asLogical(verboseR);

+  int     nr, nc;

+  SEXP    ecdfRobj;

+  double* ecdf_vals;

+  SEXP pb = R_NilValue;

+  int  nunprotect=0;

+

+  PROTECT(XR); nunprotect++;

+

+  nr = INTEGER(getAttrib(XR, R_DimSymbol))[0]; /* number of rows */

+  nc = INTEGER(getAttrib(XR, R_DimSymbol))[1]; /* number of columns */

+  X  = REAL(XR);

+

+  /* create a new dense matrix object to store the result,

+   * if nr * nc > INT_MAX and LONG_VECTOR_SUPPORT is not

+   * available, the function allocMatrix() will prompt an error */

+  ecdfRobj = PROTECT(allocMatrix(REALSXP, nr, nc)); nunprotect++;

+

+  if (verbose) {

+    pb = PROTECT(cli_progress_bar(nr, NULL));

+    cli_progress_set_name(pb, "Estimating ECDFs");

+    nunprotect++;

+  }

+

+  for (int i=0; i < nr; i++) {

+    SEXP          xR, uniqvR;

+    int           nuniqv;

+    double*       x;

+    double*       uniqv;

+    double*       ecdfuniqv;

+    int           sum;

+    double*       e1_p;

+    const double* e2_p;

+    int*          mt;

+    int*          tab;

+    int           nnas;

+

+    if (verbose) { /* show progress */

+      if (i % 100 == 0 && CLI_SHOULD_TICK)

+        cli_progress_set(pb, i);

+    }

+

+    /* remove consecutive repeated elements */

+    PROTECT(uniqvR = allocVector(REALSXP, nc));

+    PROTECT(xR = allocVector(REALSXP, nc));

+    uniqv = REAL(uniqvR);

+    x = REAL(xR);

+    nnas = 0; /* non-NA values */

+    for (int j=0; j < nc; j++) {

+#ifdef LONG_VECTOR_SUPPORT

+      R_xlen_t idx = nr * j + i;

+#else

+      int idx = nr * j + i;

+#endif

+      x[j] = X[idx];

+      if (!ISNA(x[j])) {

+        uniqv[nnas] = x[j];

+        nnas++;

+      }

+    }

+

+    qsort(uniqv, nnas, sizeof(double), dbl_cmp);

+    e1_p = uniqv;

+    e2_p = e1_p + 1;

+    nuniqv = 0;

+    for (int j=0; j < nnas; j++) {

+      if (*e2_p != *e1_p) {

+        *(++e1_p) = *e2_p;

+        nuniqv++;

+      }

+      e2_p++;

+    }

+

+    /* match original values to sorted unique values */

+    /* consider adding LONG_VECTOR_SUPPORT */

+    mt = INTEGER(match_int(xR, uniqvR)); /* 1-based! */

+

+    /* tabulate matches */

+    /* consider adding LONG_VECTOR_SUPPORT */

+    tab = R_Calloc(nuniqv, int); /* assuming zeroes are set */

+    for (int j=0; j < nc; j++)

+      if (!ISNA(mt[j]) && mt[j] > 0 && mt[j] <= nuniqv)

+        tab[mt[j] - 1]++;

+

+    /* cumulative sum to calculate ecdf values */

+    /* consider adding LONG_VECTOR_SUPPORT */

+    ecdfuniqv = R_Calloc(nuniqv, double); /* assuming zeroes are set */

+    sum = 0;

+    for (int j=0; j < nuniqv; j++) {

+      sum = sum + tab[j];

+      ecdfuniqv[j] = ((double) sum) / ((double) nc);

+    }

+

+    /* set ecdf values on the corresponding positions

+     * of the output dense matrix */

+    ecdf_vals = REAL(ecdfRobj);

+

+    for (int j=0; j < nc; j++) {

+#ifdef LONG_VECTOR_SUPPORT

+      R_xlen_t idx = nr * j + i;

+#else

+      int idx = nr * j + i;

+#endif

+      if (!ISNA(X[idx]))

+        ecdf_vals[idx] = ecdfuniqv[mt[j]-1];

+      else

+        ecdf_vals[idx] = NA_REAL;

+    }

+

+    R_Free(ecdfuniqv);

+    R_Free(tab);

+

+    UNPROTECT(2); /* xR uniqvR */

+  }

+

+  if (verbose)

+    cli_progress_done(pb);

+

+  UNPROTECT(nunprotect); /* XR ecdfRobj pb */

+

+  return(ecdfRobj);

+}

@@ -29,6 +29,8 @@ ecdfvals_sparse_to_dense_R(SEXP XCspR, SEXP XRspR, SEXP verboseR);

 SEXP

 ecdfvals_dense_to_dense_R(SEXP XR, SEXP verboseR);

+SEXP

+ecdfvals_dense_to_dense_nas_R(SEXP XR, SEXP verboseR);

 SEXP

 order_rankstat_R(SEXP xR);

@@ -56,6 +58,7 @@ static R_CallMethodDef callMethods[] = {

   {"ecdfvals_sparse_to_sparse_R", (DL_FUNC) &ecdfvals_sparse_to_sparse_R, 3},

   {"ecdfvals_sparse_to_dense_R", (DL_FUNC) &ecdfvals_sparse_to_dense_R, 3},

   {"ecdfvals_dense_to_dense_R", (DL_FUNC) &ecdfvals_dense_to_dense_R, 2},

+  {"ecdfvals_dense_to_dense_nas_R", (DL_FUNC) &ecdfvals_dense_to_dense_nas_R, 2},

   {"gsva_rnd_walk_R", (DL_FUNC) &gsva_rnd_walk_R, 3},

   {"gsva_scores_genesets_R", (DL_FUNC) &gsva_score_genesets_R, 5},

   {"order_rankstat_sparse_to_dense_R", (DL_FUNC) &order_rankstat_sparse_to_dense_R, 2},