Fixed sigma2 inverse bug; progress bar, and started R refactor

IDEA Hillal Egret Meth Lab Freeh cracking Semtex National Guard
MILSATCOM Antiviral Ortega Cornflower SBS SLIP Yuma

R/robust.R

@@ -110,32 +110,20 @@ robust_se_L <- function(cname, Y, lambda) {
 #' @export
 #' @useDynLib scdemon
 #' @importFrom Rcpp evalCpp
-robust_se_t.default <- function(V1, V2,
+robust_se_t.default <- function(Us, V,
                                 nominal_p_cutoff=0.05,
-                                lambda=1e-10,
+                                lambda=100,
                                 abs_t=FALSE, dof=NULL,
                                 t_cutoff=NULL) {
   require(Matrix)
-  if (is.null(V2)) V2 <- V1 
   if (is.null(t_cutoff)) {
-    if (is.null(dof)) {
-      stopifnot(is.integer(attr(V1, "dof")))
-      stopifnot(is.integer(attr(V2, "dof")))
-      dof <- mean(c(attr(V1, "dof"), attr(V2, "dof")))
-    }
     ## should be around 6.5 for most snRNA-seq datasets
-    t_cutoff <- qt(min(nominal_p_cutoff / (ncol(V1) * ncol(V2)), 1),
-                   dof,
+    t_cutoff <- qt(min(nominal_p_cutoff / (ncol(V) * ncol(V)), 1),
+                   max(nrow(Us) - ncol(V), 4),
                    lower.tail=FALSE)
   }
-  comm <- intersect(colnames(V2), colnames(V1))
-  M <- r_robust_se(V1, V2, lambda, t_cutoff, abs_t)
-  dimnames(M) <- list(colnames(V2), colnames(V1))
-  if (!is.null(dof)) attr(M, "dof") <- dof
-  if (length(comm) > 0) { 
-    M[cbind(match(colnames(V2), comm),
-            match(colnames(V1), comm))] <- 0
-  }
+  M <- r_robust_se(Us, V, lambda, t_cutoff, abs_t)
+  dimnames(M) <- list(colnames(V), colnames(V))
   return(Matrix::t(Matrix::drop0(M)))
 }
 

src/implprogress.hpp

@@ -47,8 +47,11 @@ class ImplProgress
 	void display()
 	{
 #if !defined(Rcpp_hpp)
-		double last_percent = last / total;
-		double cur_percent = current / total;
+		double last_percent = (double)last / total;
+		double cur_percent = (double)current / total;
+#if defined(_OPENMP)
+#pragma omp critical
+#endif
 		if (cur_percent >= 0.01 + last_percent) {
 			int pct = 100 * cur_percent;
 			int pos = cur_percent * BAR_WIDTH;
@@ -57,8 +60,8 @@ class ImplProgress
 				if (i < pos) { out[i] = '='; }
 				else if (i == pos) { out[i] = '>'; }
 			}
-			std::cout << "[" << out << "] " << pct << "%\r";
-			std::cout.flush();
+			std::cerr << "[" << out << "] " << pct << "%\r";
+			std::cerr.flush();
 			last = current;
 		}
 #endif	

src/r_se.cpp

@@ -7,59 +7,11 @@
 
 #include "se.hpp"
 
-
 // [[Rcpp::export]]
-Eigen::MatrixXd r_ols_beta(const Eigen::Map<Eigen::MatrixXd> &X,
-			   const Eigen::Map<Eigen::MatrixXd> &Y,
-                           double lambda)
-{
-	return ols_beta(X, Y, lambda);
-}
-
-// [[Rcpp::export]]
-Eigen::MatrixXd r_ols_resid(const Eigen::Map<Eigen::MatrixXd> &X,
-			    const Eigen::Map<Eigen::MatrixXd> &Y,
-			    const Eigen::Map<Eigen::MatrixXd> &beta)
-{
-	return ols_resid(X, Y, beta);
-}
-
-// [[Rcpp::export]]
-Eigen::ArrayXd r_robust_se_X(const Eigen::Map<Eigen::MatrixXd> &X,
-			     const Eigen::Map<Eigen::MatrixXd> &Y,
-			     double lambda=1e-10)
-{
-	return robust_se_X(X, Y, lambda, 1e-300);
-}
-
-// [[Rcpp::export]]
-Eigen::ArrayXd r_robust_se_Y(const Eigen::Map<Eigen::MatrixXd> &X,
-			     const Eigen::Map<Eigen::MatrixXd> &Y,
-			     double lambda=1e-10)
-{
-	return robust_se_Y(X, Y, lambda, 1e-300);
-}
-
-// [[Rcpp::export]]
-Eigen::MatrixXd r_ols_beta_L(const Eigen::Map<Eigen::VectorXd> &X,
-                             const Eigen::Map<Eigen::MatrixXd> &Y,
-                             const Eigen::Map<Eigen::ArrayXd> &lambda)
-{
-	return ols_beta_L(X.eval(), X.squaredNorm(), Y, lambda);
-}
-// [[Rcpp::export]]
-Eigen::ArrayXd r_robust_se_L(const Eigen::Map<Eigen::MatrixXd> &X,
-			     const Eigen::Map<Eigen::MatrixXd> &Y,
-                             const Eigen::Map<Eigen::ArrayXd> &lambda)
-{
-	return robust_se_L(X, Y, lambda, 1e-300);
-}
-
-// [[Rcpp::export]]
-Eigen::SparseMatrix<double> r_robust_se(const Eigen::Map<Eigen::MatrixXd> &X,
-					const Eigen::Map<Eigen::MatrixXd> &Y,
+Eigen::SparseMatrix<double> r_robust_se(const Eigen::Map<Eigen::MatrixXd> &U,
+					const Eigen::Map<Eigen::MatrixXd> &V,
 					double lambda=1e-10,
 					double t_cutoff=2, bool abs_cutoff=false)
 {
-	return robust_se(X, Y, lambda, 1e-300, t_cutoff, abs_cutoff);
+	return intra_robust_se(U, V, lambda, t_cutoff, abs_cutoff);
 }

src/se.hpp

@@ -66,7 +66,7 @@ Eigen::SparseMatrix<float> intra_robust_se(const Eigen::MatrixBase<TU> &U,
 					   bool abs_cutoff=false)
 {
 	const Eigen::VectorXf sigma_2 = U.colwise().squaredNorm().template cast<float>().eval();
-	const Eigen::ArrayXf sigma_2_inv = (1. / (sigma_2.array() + lambda*lambda));
+	const Eigen::ArrayXf sigma_2_inv = 1. / sigma_2.array();
 	const Eigen::DiagonalMatrix<float, Eigen::Dynamic> tuui = sigma_2_inv.matrix().asDiagonal();
 	const Eigen::DiagonalMatrix<typename TV::Scalar, Eigen::Dynamic> tuuituu = (sigma_2.array() / (sigma_2.array() + lambda*lambda)).matrix().template cast<typename TV::Scalar>().asDiagonal();
 	Eigen::SparseMatrix<float> M(V.cols(), V.cols());