From 1ac9ee2eaa18047194cfa3067d6e73565d8e9e02 Mon Sep 17 00:00:00 2001 From: zdebruine Date: Mon, 18 Dec 2023 10:09:07 -0500 Subject: [PATCH] support for c(L1_w, L1_h) and c(L2_w, L2_h) --- R/RcppExports.R | 16 ++++-- R/run_nmf.R | 14 +++-- src/RcppExports.cpp | 58 +++++++++++++++---- src/singlet.cpp | 134 +++++++++++++++++++++++++++++++++++++++++--- 4 files changed, 194 insertions(+), 28 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 5a7cc89..45037a6 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -13,6 +13,10 @@ rowwise_compress_dense <- function(A, n = 10L, threads = 0L) { .Call(`_singlet_rowwise_compress_dense`, A, n, threads) } +calc_L1_matrix <- function(h, batch_id) { + .Call(`_singlet_calc_L1_matrix`, h, batch_id) +} + Rcpp_predict <- function(A, w, L1, L2, threads) { .Call(`_singlet_Rcpp_predict`, A, w, L1, L2, threads) } @@ -21,8 +25,12 @@ c_project_model <- function(A, w, L1, L2, threads) { .Call(`_singlet_c_project_model`, A, w, L1, L2, threads) } -c_nmf <- function(A, At, tol, maxit, verbose, L1, L2, threads, w) { - .Call(`_singlet_c_nmf`, A, At, tol, maxit, verbose, L1, L2, threads, w) +c_nmf <- function(A, At, tol, maxit, verbose, L1_w, L1_h, L2_w, L2_h, threads, w) { + .Call(`_singlet_c_nmf`, A, At, tol, maxit, verbose, L1_w, L1_h, L2_w, L2_h, threads, w) +} + +c_nmf_batch <- function(A, At, tol, maxit, verbose, L1, L2, threads, w, batch_id) { + .Call(`_singlet_c_nmf_batch`, A, At, tol, maxit, verbose, L1, L2, threads, w, batch_id) } c_nmf_sparse_list <- function(A_, At_, tol, maxit, verbose, L1, L2, threads, w) { @@ -59,8 +67,8 @@ c_mu_nmf <- function(A, At, tol, maxit, verbose, L1, L2, threads, w) { .Call(`_singlet_c_mu_nmf`, A, At, tol, maxit, verbose, L1, L2, threads, w) } -c_nmf_dense <- function(A, At, tol, maxit, verbose, L1, L2, threads, w) { - .Call(`_singlet_c_nmf_dense`, A, At, tol, maxit, verbose, L1, L2, threads, w) +c_nmf_dense <- function(A, At, tol, maxit, verbose, L1_w, L1_h, L2_w, L2_h, threads, w) { + .Call(`_singlet_c_nmf_dense`, A, At, tol, maxit, verbose, L1_w, L1_h, L2_w, L2_h, threads, w) } c_linked_nmf <- function(A, At, tol, maxit, verbose, L1, L2, threads, w, link_h, link_w) { diff --git a/R/run_nmf.R b/R/run_nmf.R index e09b4d7..3d7b5c3 100644 --- a/R/run_nmf.R +++ b/R/run_nmf.R @@ -17,9 +17,6 @@ #' run_nmf <- function(A, rank, tol = 1e-4, maxit = 100, verbose = TRUE, L1 = 0.01, L2 = 0, threads = 0, compression_level = 3) { use_vcsc <- compression_level == 2 - if (L1 >= 1) { - stop("L1 penalty must be strictly in the range (0, 1]") - } if ("list" %in% class(A)) { # check that number of rows is identical @@ -47,12 +44,19 @@ run_nmf <- function(A, rank, tol = 1e-4, maxit = 100, verbose = TRUE, L1 = 0.01, At <- t(A) dense_mode <- TRUE } + + if(length(L1) != 2){ + L1 <- c(L1[[1]], L1[[1]]) + } + if(length(L2) != 2){ + L2 <- c(L2[[1]], L2[[1]]) + } w_init <- matrix(stats::runif(nrow(A) * rank), rank, nrow(A)) if (dense_mode) { - model <- c_nmf_dense(A, At, tol, maxit, verbose, L1, L2, threads, w_init) + model <- c_nmf_dense(A, At, tol, maxit, verbose, L1[[1]], L1[[2]], L2[[1]], L2[[2]], threads, w_init) } else { - model <- c_nmf(A, At, tol, maxit, verbose, L1, L2, threads, w_init) + model <- c_nmf(A, At, tol, maxit, verbose, L1[[1]], L1[[2]], L2[[1]], L2[[2]], threads, w_init) } rn <- rownames(A) cn <- colnames(A) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 5c1ef56..7a86710 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -51,6 +51,18 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// calc_L1_matrix +Eigen::MatrixXd calc_L1_matrix(const Eigen::MatrixXd& h, Rcpp::IntegerVector batch_id); +RcppExport SEXP _singlet_calc_L1_matrix(SEXP hSEXP, SEXP batch_idSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const Eigen::MatrixXd& >::type h(hSEXP); + Rcpp::traits::input_parameter< Rcpp::IntegerVector >::type batch_id(batch_idSEXP); + rcpp_result_gen = Rcpp::wrap(calc_L1_matrix(h, batch_id)); + return rcpp_result_gen; +END_RCPP +} // Rcpp_predict Eigen::MatrixXd Rcpp_predict(Rcpp::SparseMatrix A, Eigen::MatrixXd w, const double L1, const double L2, const int threads); RcppExport SEXP _singlet_Rcpp_predict(SEXP ASEXP, SEXP wSEXP, SEXP L1SEXP, SEXP L2SEXP, SEXP threadsSEXP) { @@ -82,8 +94,29 @@ BEGIN_RCPP END_RCPP } // c_nmf -Rcpp::List c_nmf(Rcpp::SparseMatrix& A, Rcpp::SparseMatrix& At, const double tol, const uint16_t maxit, const bool verbose, const double L1, const double L2, const uint16_t threads, Eigen::MatrixXd w); -RcppExport SEXP _singlet_c_nmf(SEXP ASEXP, SEXP AtSEXP, SEXP tolSEXP, SEXP maxitSEXP, SEXP verboseSEXP, SEXP L1SEXP, SEXP L2SEXP, SEXP threadsSEXP, SEXP wSEXP) { +Rcpp::List c_nmf(Rcpp::SparseMatrix& A, Rcpp::SparseMatrix& At, const double tol, const uint16_t maxit, const bool verbose, const double L1_w, const double L1_h, const double L2_w, const double L2_h, const uint16_t threads, Eigen::MatrixXd w); +RcppExport SEXP _singlet_c_nmf(SEXP ASEXP, SEXP AtSEXP, SEXP tolSEXP, SEXP maxitSEXP, SEXP verboseSEXP, SEXP L1_wSEXP, SEXP L1_hSEXP, SEXP L2_wSEXP, SEXP L2_hSEXP, SEXP threadsSEXP, SEXP wSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::SparseMatrix& >::type A(ASEXP); + Rcpp::traits::input_parameter< Rcpp::SparseMatrix& >::type At(AtSEXP); + Rcpp::traits::input_parameter< const double >::type tol(tolSEXP); + Rcpp::traits::input_parameter< const uint16_t >::type maxit(maxitSEXP); + Rcpp::traits::input_parameter< const bool >::type verbose(verboseSEXP); + Rcpp::traits::input_parameter< const double >::type L1_w(L1_wSEXP); + Rcpp::traits::input_parameter< const double >::type L1_h(L1_hSEXP); + Rcpp::traits::input_parameter< const double >::type L2_w(L2_wSEXP); + Rcpp::traits::input_parameter< const double >::type L2_h(L2_hSEXP); + Rcpp::traits::input_parameter< const uint16_t >::type threads(threadsSEXP); + Rcpp::traits::input_parameter< Eigen::MatrixXd >::type w(wSEXP); + rcpp_result_gen = Rcpp::wrap(c_nmf(A, At, tol, maxit, verbose, L1_w, L1_h, L2_w, L2_h, threads, w)); + return rcpp_result_gen; +END_RCPP +} +// c_nmf_batch +Rcpp::List c_nmf_batch(Rcpp::SparseMatrix& A, Rcpp::SparseMatrix& At, const double tol, const uint16_t maxit, const bool verbose, const double L1, const double L2, const uint16_t threads, Eigen::MatrixXd w, Rcpp::IntegerVector batch_id); +RcppExport SEXP _singlet_c_nmf_batch(SEXP ASEXP, SEXP AtSEXP, SEXP tolSEXP, SEXP maxitSEXP, SEXP verboseSEXP, SEXP L1SEXP, SEXP L2SEXP, SEXP threadsSEXP, SEXP wSEXP, SEXP batch_idSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -96,7 +129,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const double >::type L2(L2SEXP); Rcpp::traits::input_parameter< const uint16_t >::type threads(threadsSEXP); Rcpp::traits::input_parameter< Eigen::MatrixXd >::type w(wSEXP); - rcpp_result_gen = Rcpp::wrap(c_nmf(A, At, tol, maxit, verbose, L1, L2, threads, w)); + Rcpp::traits::input_parameter< Rcpp::IntegerVector >::type batch_id(batch_idSEXP); + rcpp_result_gen = Rcpp::wrap(c_nmf_batch(A, At, tol, maxit, verbose, L1, L2, threads, w, batch_id)); return rcpp_result_gen; END_RCPP } @@ -204,8 +238,8 @@ BEGIN_RCPP END_RCPP } // c_nmf_dense -Rcpp::List c_nmf_dense(Eigen::MatrixXd& A, Eigen::MatrixXd& At, const double tol, const uint16_t maxit, const bool verbose, const double L1, const double L2, const uint16_t threads, Eigen::MatrixXd w); -RcppExport SEXP _singlet_c_nmf_dense(SEXP ASEXP, SEXP AtSEXP, SEXP tolSEXP, SEXP maxitSEXP, SEXP verboseSEXP, SEXP L1SEXP, SEXP L2SEXP, SEXP threadsSEXP, SEXP wSEXP) { +Rcpp::List c_nmf_dense(Eigen::MatrixXd& A, Eigen::MatrixXd& At, const double tol, const uint16_t maxit, const bool verbose, const double L1_w, const double L1_h, const double L2_w, const double L2_h, const uint16_t threads, Eigen::MatrixXd w); +RcppExport SEXP _singlet_c_nmf_dense(SEXP ASEXP, SEXP AtSEXP, SEXP tolSEXP, SEXP maxitSEXP, SEXP verboseSEXP, SEXP L1_wSEXP, SEXP L1_hSEXP, SEXP L2_wSEXP, SEXP L2_hSEXP, SEXP threadsSEXP, SEXP wSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -214,11 +248,13 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const double >::type tol(tolSEXP); Rcpp::traits::input_parameter< const uint16_t >::type maxit(maxitSEXP); Rcpp::traits::input_parameter< const bool >::type verbose(verboseSEXP); - Rcpp::traits::input_parameter< const double >::type L1(L1SEXP); - Rcpp::traits::input_parameter< const double >::type L2(L2SEXP); + Rcpp::traits::input_parameter< const double >::type L1_w(L1_wSEXP); + Rcpp::traits::input_parameter< const double >::type L1_h(L1_hSEXP); + Rcpp::traits::input_parameter< const double >::type L2_w(L2_wSEXP); + Rcpp::traits::input_parameter< const double >::type L2_h(L2_hSEXP); Rcpp::traits::input_parameter< const uint16_t >::type threads(threadsSEXP); Rcpp::traits::input_parameter< Eigen::MatrixXd >::type w(wSEXP); - rcpp_result_gen = Rcpp::wrap(c_nmf_dense(A, At, tol, maxit, verbose, L1, L2, threads, w)); + rcpp_result_gen = Rcpp::wrap(c_nmf_dense(A, At, tol, maxit, verbose, L1_w, L1_h, L2_w, L2_h, threads, w)); return rcpp_result_gen; END_RCPP } @@ -409,9 +445,11 @@ static const R_CallMethodDef CallEntries[] = { {"_singlet_weight_by_split", (DL_FUNC) &_singlet_weight_by_split, 3}, {"_singlet_rowwise_compress_sparse", (DL_FUNC) &_singlet_rowwise_compress_sparse, 3}, {"_singlet_rowwise_compress_dense", (DL_FUNC) &_singlet_rowwise_compress_dense, 3}, + {"_singlet_calc_L1_matrix", (DL_FUNC) &_singlet_calc_L1_matrix, 2}, {"_singlet_Rcpp_predict", (DL_FUNC) &_singlet_Rcpp_predict, 5}, {"_singlet_c_project_model", (DL_FUNC) &_singlet_c_project_model, 5}, - {"_singlet_c_nmf", (DL_FUNC) &_singlet_c_nmf, 9}, + {"_singlet_c_nmf", (DL_FUNC) &_singlet_c_nmf, 11}, + {"_singlet_c_nmf_batch", (DL_FUNC) &_singlet_c_nmf_batch, 10}, {"_singlet_c_nmf_sparse_list", (DL_FUNC) &_singlet_c_nmf_sparse_list, 9}, {"_singlet_write_IVCSC", (DL_FUNC) &_singlet_write_IVCSC, 2}, {"_singlet_save_IVSparse", (DL_FUNC) &_singlet_save_IVSparse, 2}, @@ -419,7 +457,7 @@ static const R_CallMethodDef CallEntries[] = { {"_singlet_read_IVSparse", (DL_FUNC) &_singlet_read_IVSparse, 0}, {"_singlet_run_nmf_on_sparsematrix_list", (DL_FUNC) &_singlet_run_nmf_on_sparsematrix_list, 9}, {"_singlet_c_mu_nmf", (DL_FUNC) &_singlet_c_mu_nmf, 9}, - {"_singlet_c_nmf_dense", (DL_FUNC) &_singlet_c_nmf_dense, 9}, + {"_singlet_c_nmf_dense", (DL_FUNC) &_singlet_c_nmf_dense, 11}, {"_singlet_c_linked_nmf", (DL_FUNC) &_singlet_c_linked_nmf, 11}, {"_singlet_c_ard_nmf", (DL_FUNC) &_singlet_c_ard_nmf, 13}, {"_singlet_c_ard_nmf_sparse_list", (DL_FUNC) &_singlet_c_ard_nmf_sparse_list, 13}, diff --git a/src/singlet.cpp b/src/singlet.cpp index 1c5aab0..4560f8c 100644 --- a/src/singlet.cpp +++ b/src/singlet.cpp @@ -226,14 +226,14 @@ void scale(Eigen::MatrixXd& w, Eigen::VectorXd& d) { // NNLS SOLVER OF THE FORM ax=b --------------------------------------------------------------------------------- // optimized and modified from github.com/linxihui/NNLM "c_nnls" function -inline void nnls(Eigen::MatrixXd& a, Eigen::VectorXd& b, Eigen::MatrixXd& x, const size_t col, const double L1 = 0, const double L2 = 0) { +inline void nnls(Eigen::MatrixXd& a, Eigen::VectorXd& b, Eigen::MatrixXd& x, const size_t col,const double L1 = 0, const double L2 = 0) { double tol = 1; for (uint8_t it = 0; it < 100 && (tol / b.size()) > 1e-8; ++it) { tol = 0; for (size_t i = 0; i < x.rows(); ++i) { double diff = b(i) / a(i, i); if (L1 != 0) diff -= L1; - if (L2 != 0) diff -= L2 * x(i, col); + if (L2 != 0) diff += L2 * x(i, col); if (-diff > x(i, col)) { if (x(i, col) != 0) { b -= a.col(i) * -x(i, col); @@ -249,6 +249,84 @@ inline void nnls(Eigen::MatrixXd& a, Eigen::VectorXd& b, Eigen::MatrixXd& x, con } } +// NNLS SOLVER OF THE FORM ax=b --------------------------------------------------------------------------------- +// optimized and modified from github.com/linxihui/NNLM "c_nnls" function +inline void nnls_L1_matrix(Eigen::MatrixXd& a, Eigen::VectorXd& b, Eigen::MatrixXd& x, const size_t col, const Eigen::MatrixXd& L1_matrix, const double L1 = 0, const double L2 = 0) { + double tol = 1; + for (uint8_t it = 0; it < 100 && (tol / b.size()) > 1e-8; ++it) { + tol = 0; + for (size_t i = 0; i < x.rows(); ++i) { + double diff = b(i) / a(i, i); + diff -= L1_matrix(i, col); + if (L1 != 0) diff -= L1; + if (L2 != 0) diff += L2 * x(i, col); + if (-diff > x(i, col)) { + if (x(i, col) != 0) { + b -= a.col(i) * -x(i, col); + tol = 1; + x(i, col) = 0; + } + } else if (diff != 0) { + x(i, col) += diff; + b -= a.col(i) * diff; + tol += std::abs(diff / (x(i, col) + 1e-15)); + } + } + } +} + +// batch_id is a vector with values between 1 and n_batch_ids, the result of +// as.numeric(as.factor(metadata_item_vector)) +//[[Rcpp::export]] +Eigen::MatrixXd calc_L1_matrix(const Eigen::MatrixXd& h, Rcpp::IntegerVector batch_id) { + const size_t n_batch_ids = Rcpp::max(batch_id); + if (batch_id.size() != h.cols()) Rcpp::stop("batch_id vector must be of the same length as the number of columns in the NMF 'h' matrix"); + Eigen::MatrixXd L1_matrix(h.rows(), n_batch_ids); + + // calculate the mean factor loading for each batch + for (size_t i = 1; i < n_batch_ids; ++i) { + size_t n_samples_in_batch = 0; + for (size_t j = 0; j < batch_id.size(); ++j) { + if (batch_id[j] == i) { + L1_matrix.col(i).array() += h.col(j).array(); + ++n_samples_in_batch; + } + } + L1_matrix.col(i).array() /= n_samples_in_batch; + } + + // calculate the difference in the mean factor loading for each batch vs. all other batches + for (size_t i = 0; i < L1_matrix.cols(); ++i) { + Eigen::VectorXd mean_other_batches = Eigen::VectorXd::Zero(L1_matrix.rows()); + for (size_t j = 0; j < L1_matrix.cols(); ++j) { + if (j != i) { + mean_other_batches.array() += L1_matrix.col(j).array(); + } + } + mean_other_batches /= (L1_matrix.cols() - 1); + L1_matrix.col(i) -= mean_other_batches; + } + + return L1_matrix; +} + +// update h given A and w, balanced for batch identity +inline void predict_L1_matrix(Rcpp::SparseMatrix A, const Eigen::MatrixXd& w, Eigen::MatrixXd& h, const double L1, const double L2, const int threads, Rcpp::IntegerVector batch_id) { + Eigen::MatrixXd a = AAt(w); + // calculate difference between mean loading of in-batch factor vs. out-of-batch factors + Eigen::MatrixXd L1_matrix = calc_L1_matrix(h, batch_id); +#ifdef _OPENMP +#pragma omp parallel for num_threads(threads) +#endif + for (size_t i = 0; i < h.cols(); ++i) { + if (A.p[i] == A.p[i + 1]) continue; + Eigen::VectorXd b = Eigen::VectorXd::Zero(h.rows()); + for (Rcpp::SparseMatrix::InnerIterator it(A, i); it; ++it) + b += it.value() * w.col(it.row()); + nnls_L1_matrix(a, b, h, i, L1_matrix, L1, L2); + } +} + // NMF PROJECTION FUNCTIONS ------------------------------------------------------------------------------ // update h given A and w @@ -558,7 +636,7 @@ inline double mse_test(const Eigen::MatrixXd& A, const Eigen::MatrixXd& w, Eigen // NMF FUNCTIONS --------------------------------------------------------------------------------------- // no linking or masking template -Rcpp::List c_nmf_base(Matrix& A, Matrix& At, const double tol, const uint16_t maxit, const bool verbose, const double L1, const double L2, const uint16_t threads, Eigen::MatrixXd w) { +Rcpp::List c_nmf_base(Matrix& A, Matrix& At, const double tol, const uint16_t maxit, const bool verbose, const double L1_w, const double L1_h, const double L2_w, const double L2_h, const uint16_t threads, Eigen::MatrixXd w) { Eigen::MatrixXd h = Eigen::MatrixXd::Zero(w.rows(), A.cols()); Eigen::VectorXd d = Eigen::VectorXd::Ones(w.rows()); double tol_ = 1; @@ -569,11 +647,11 @@ Rcpp::List c_nmf_base(Matrix& A, Matrix& At, const double tol, const uint16_t ma for (uint16_t iter_ = 0; iter_ < maxit && tol_ > tol; ++iter_) { Eigen::MatrixXd w_it = w; // update h - predict(A, w, h, L1, L2, threads); + predict(A, w, h, L1_h, L2_h, threads); scale(h, d); Rcpp::checkUserInterrupt(); // update w - predict(At, h, w, L1, L2, threads); + predict(At, h, w, L1_w, L2_w, threads); scale(w, d); // calculate tolerance of the model fit to detect convergence @@ -589,8 +667,46 @@ Rcpp::List c_nmf_base(Matrix& A, Matrix& At, const double tol, const uint16_t ma //[[Rcpp::export]] Rcpp::List c_nmf(Rcpp::SparseMatrix& A, Rcpp::SparseMatrix& At, const double tol, const uint16_t maxit, const bool verbose, - const double L1, const double L2, const uint16_t threads, Eigen::MatrixXd w) { - return c_nmf_base(A, At, tol, maxit, verbose, L1, L2, threads, w); + const double L1_w, const double L1_h, const double L2_w, const double L2_h, const uint16_t threads, Eigen::MatrixXd w) { + return c_nmf_base(A, At, tol, maxit, verbose, L1_w, L1_h, L2_w, L2_h, threads, w); +} + +// L1 MATRIX-BASED BATCH CORRECTION (EXPERIMENTAL) + +template +Rcpp::List c_nmf_base_batch(Matrix& A, Matrix& At, const double tol, const uint16_t maxit, const bool verbose, const double L1, const double L2, const uint16_t threads, Eigen::MatrixXd w, Rcpp::IntegerVector batch_id) { + Eigen::MatrixXd h = Eigen::MatrixXd::Zero(w.rows(), A.cols()); + Eigen::VectorXd d = Eigen::VectorXd::Ones(w.rows()); + double tol_ = 1; + if (verbose) + Rprintf("\n%4s | %8s \n---------------\n", "iter", "tol"); + + // alternating least squares update loop + for (uint16_t iter_ = 0; iter_ < maxit && tol_ > tol; ++iter_) { + Eigen::MatrixXd w_it = w; + // update h + predict_L1_matrix(A, w, h, L1, L2, threads, batch_id); + scale(h, d); + Rcpp::checkUserInterrupt(); + // update w + predict(At, h, w, L1, L2, threads); + scale(w, d); + + // calculate tolerance of the model fit to detect convergence + // correlation between "w" across consecutive iterations + tol_ = cor(w, w_it); + + if (verbose) + Rprintf("%4d | %8.2e\n", iter_ + 1, tol_); + Rcpp::checkUserInterrupt(); + } + return Rcpp::List::create(Rcpp::Named("w") = w, Rcpp::Named("d") = d, Rcpp::Named("h") = h); +} + +//[[Rcpp::export]] +Rcpp::List c_nmf_batch(Rcpp::SparseMatrix& A, Rcpp::SparseMatrix& At, const double tol, const uint16_t maxit, const bool verbose, + const double L1, const double L2, const uint16_t threads, Eigen::MatrixXd w, Rcpp::IntegerVector batch_id) { + return c_nmf_base_batch(A, At, tol, maxit, verbose, L1, L2, threads, w, batch_id); } // NMF FUNCTIONS --------------------------------------------------------------------------------------- @@ -933,8 +1049,8 @@ Rcpp::List c_mu_nmf(Rcpp::SparseMatrix& A, Rcpp::SparseMatrix& At, const double } //[[Rcpp::export]] -Rcpp::List c_nmf_dense(Eigen::MatrixXd& A, Eigen::MatrixXd& At, const double tol, const uint16_t maxit, const bool verbose, const double L1, const double L2, const uint16_t threads, Eigen::MatrixXd w) { - return c_nmf_base(A, At, tol, maxit, verbose, L1, L2, threads, w); +Rcpp::List c_nmf_dense(Eigen::MatrixXd& A, Eigen::MatrixXd& At, const double tol, const uint16_t maxit, const bool verbose, const double L1_w, const double L1_h, const double L2_w, const double L2_h, const uint16_t threads, Eigen::MatrixXd w) { + return c_nmf_base(A, At, tol, maxit, verbose, L1_w, L1_h, L2_w, L2_h, threads, w); } // NMF FUNCTION ---------------------------------------------------------------------------------------