re-add nnls

Author: theAeon

nmf/CMakeLists.txt

@@ -9,7 +9,7 @@ include(FetchContent)
 #add_library(sparse_nmf OBJECT SparseNMFDriver.cpp)
 #add_library(dense_nmf OBJECT NMFDriver.cpp)
 add_library(nmflib
-        nmf_lib.cpp bppnmf.cpp
+        nmf_lib.cpp bppnmf.cpp ../nnls/nnls_lib.cpp
         ../common/utils.cc ../common/data.cpp)
 # add_executable(${DENSE_OR_SPARSE}_inmf inmf.cpp)
 # target_compile_features(${DENSE_OR_SPARSE}_inmf PRIVATE cxx_std_17)

nnls/nnls_lib.cpp

@@ -0,0 +1,11 @@
+//
+// Created by andrew on 3/21/2025.
+//
+
+#include "nnls_lib.inl"
+
+#define X(T) \
+template arma::mat planc::nnlslib<T, double>::runbppnnls(const arma::mat &C, const T &B, const int &ncores);
+#include "../nmf/nmf_types.inc"
+#undef X
+template arma::mat planc::nnlslib<arma::mat>::bppnnls_prod(const arma::mat &CtC, const arma::mat &CtB, const int &ncores);

nnls/nnls_lib.hpp

@@ -0,0 +1,32 @@
+//
+// Created by andrew on 3/21/2025.
+//
+
+#ifndef NNLS_LIB_H
+#define NNLS_LIB_H
+
+#include "../nmf/nmf_lib.hpp"
+
+namespace planc {
+
+    template<typename T, typename eT = typename T::elem_type>
+    class NMFLIB_EXPORT nnlslib {
+
+public:
+    nnlslib() {
+        openblas_pthread_off(get_openblas_handle());
+    }
+
+    ~nnlslib() = default;
+
+
+    static arma::mat runbppnnls(const arma::mat &C, const T &B, const int &ncores);
+
+    static arma::mat bppnnls_prod(const arma::mat &CtC, const arma::mat &CtB, const int& nCores = 2);
+
+
+};
+
+} // planc
+
+#endif //NNLS_LIB_H

nnls/nnls_lib.inl

@@ -0,0 +1,52 @@
+#pragma once
+
+#include "nnls_lib.hpp"
+#include "utils.hpp"
+#include "bppnnls.hpp"
+
+template<typename T, typename eT>
+arma::mat planc::nnlslib<T, eT>::runbppnnls(const arma::mat &C, const T &B, const int& ncores) {
+    arma::uword m_n = B.n_cols;
+    arma::uword m_k = C.n_cols;
+    arma::mat CtC = C.t() * C;
+    arma::mat outmat = arma::zeros<arma::mat>(m_k, m_n);
+    arma::mat *outmatptr;
+    outmatptr = &outmat;
+    arma::uword ONE_THREAD_MATRIX_SIZE = chunk_size_dense<double>(m_k);
+    unsigned int numChunks = m_n / ONE_THREAD_MATRIX_SIZE;
+    if (numChunks*ONE_THREAD_MATRIX_SIZE < m_n) numChunks++;
+#pragma omp parallel for schedule(dynamic) default(none) shared(numChunks, ONE_THREAD_MATRIX_SIZE, m_n, outmatptr, C, B, CtC) num_threads(ncores)
+     for (unsigned int i = 0; i < numChunks; i++) {
+         unsigned int spanStart = i * ONE_THREAD_MATRIX_SIZE;
+         unsigned int spanEnd = (i + 1) * ONE_THREAD_MATRIX_SIZE - 1;
+         if (spanEnd > m_n - 1) spanEnd = m_n - 1;
+         // double start = omp_get_wtime();
+         arma::mat CtBChunk = C.t() * B.cols(spanStart, spanEnd);
+         BPPNNLS<arma::mat, arma::vec> solveProblem(CtC, CtBChunk, true);
+         solveProblem.solveNNLS();
+         (*outmatptr).cols(spanStart, spanEnd) = solveProblem.getSolutionMatrix();
+     }
+     return outmat;
+ }
+
+template<typename T, typename eT>
+arma::mat planc::nnlslib<T, eT>::bppnnls_prod(const arma::mat &CtC, const arma::mat &CtB, const int& nCores) {
+        arma::uword n = CtB.n_cols;
+        arma::uword k = CtC.n_cols;
+        arma::uword ONE_THREAD_MATRIX_SIZE = chunk_size_dense<double>(k);
+        arma::mat outmat = arma::zeros<arma::mat>(k, n);
+        arma::mat* outmatptr = &outmat;
+        unsigned int numChunks = n / ONE_THREAD_MATRIX_SIZE;
+        if (numChunks*ONE_THREAD_MATRIX_SIZE < n) numChunks++;
+#pragma omp parallel for schedule(dynamic) default(none) shared(numChunks, CtB, ONE_THREAD_MATRIX_SIZE, outmatptr, CtC, n) num_threads(nCores)
+        for (unsigned int i = 0; i < numChunks; i++) {
+            unsigned int spanStart = i * ONE_THREAD_MATRIX_SIZE;
+            unsigned int spanEnd = (i + 1) * ONE_THREAD_MATRIX_SIZE - 1;
+            if (spanEnd > n - 1) spanEnd = n - 1;
+            arma::mat CtBChunk = CtB.cols(spanStart, spanEnd);
+            BPPNNLS<arma::mat, arma::vec> solveProblem(CtC, CtBChunk, true);
+            solveProblem.solveNNLS();
+            (*outmatptr).cols(spanStart, spanEnd) = solveProblem.getSolutionMatrix();
+        }
+        return outmat;
+    }