Merge branch 'master' of https://github.com/welch-lab/libplanc

Author: theAeon

common/data.cpp

@@ -415,18 +415,18 @@ namespace planc {
             this->p_chunksize = p_chunkdim[0];
             // p_cparms.close();
             // #ifdef _VERBOSE
-#ifdef USING_R
-            Rcpp::Rcout
-#else
-            std::cout
-#endif
-                    << "==H5SpMat constructed==" << std::endl
-                    << "H5File:    " << filename << std::endl
-                    << "colptr path:  " << pPath << std::endl
-                    << "rowind path:  " << iPath << std::endl
-                    << "value path:   " << xPath << std::endl
-                    << "Dimension: " << n_rows << " x " << n_cols << std::endl;
-            // #endif
+// #ifdef USING_R
+//             Rcpp::Rcout
+// #else
+//             std::cout
+// #endif
+//                     << "==H5SpMat constructed==" << std::endl
+//                     << "H5File:    " << filename << std::endl
+//                     << "colptr path:  " << pPath << std::endl
+//                     << "rowind path:  " << iPath << std::endl
+//                     << "value path:   " << xPath << std::endl
+//                     << "Dimension: " << n_rows << " x " << n_cols << std::endl;
+//             // #endif
         }
 
         ~H5SpMatImpl() {

common/inmf.hpp

@@ -26,6 +26,8 @@ namespace planc {
         std::unique_ptr<arma::mat> WT; // kxm
         double lambda, sqrtLambda, objective_err;
         bool cleared;
+        std::unique_ptr<arma::sp_mat> tempE;
+        void load_whole_E(arma::uword i);
 
         virtual double computeObjectiveError() {
             // obj_i = ||E_i - (W + V_i)*H_i||_F^2 + lambda * ||V_i*H_i||_F^2
@@ -39,16 +41,20 @@ namespace planc {
             arma::mat* Wptr = this->W.get();
             arma::mat L(this->m, this->k); // (loading) L = W + V
             for (arma::uword i = 0; i < this->nDatasets; ++i) {
-                T* Eptr = this->Ei[i].get();
+                // T* Eptr = this->Ei[i].get();
+                // arma::sp_mat* Eptr = this->load_whole_E(i);
+                this->load_whole_E(i);
+                arma::sp_mat* Eptr = this->tempE.get();
                 arma::mat* Hptr = this->Hi[i].get();
                 arma::mat* Vptr = this->Vi[i].get();
                 L = *Wptr + *Vptr;
-                double sqnormE = arma::norm<T>(*Eptr, "fro");
+                double sqnormE = arma::norm<arma::sp_mat>(*Eptr, "fro");
                 sqnormE *= sqnormE;
                 arma::mat LtL = L.t() * L; // k x k
                 arma::mat HtH = Hptr->t() * *Hptr; // k x k
                 double TrLtLHtH = arma::trace(LtL * HtH);
-                T Et = Eptr->t();
+                // T Et = Eptr->t();
+                arma::sp_mat Et = Eptr->t();
                 arma::mat EtL = Et * L; // n_i x k
                 double TrHtEtL = arma::trace(Hptr->t() * EtL);
                 arma::mat VtV = Vptr->t() * *Vptr; // k x k
@@ -359,8 +365,75 @@ namespace planc {
                 }
                 this->W.reset();
                 if (this->WT.get() != nullptr) this->WT.reset();
+                this->tempE.reset();
             }
             this->cleared = true;
         }
     }; // class INMF
+
+    template<>
+    inline void INMF<arma::sp_mat>::load_whole_E(arma::uword i) {
+        // Make a real matrix copied from what's pointed by this->Ei[i]
+        arma::sp_mat tempSparse = *(this->Ei[i].get());
+        // Make new unique_ptr to manage the memory
+        auto tempSparseUniquePtr = std::make_unique<arma::sp_mat>(tempSparse);
+        this->tempE = std::move(tempSparseUniquePtr);
+        // return this->Ei[i].get();
+    }
+
+    template<>
+    inline void INMF<arma::mat>::load_whole_E(arma::uword i) {
+        auto tempSparse = std::make_unique<arma::sp_mat>(*(this->Ei[i]));
+
+        // Get raw pointer before transferring ownership
+        // arma::sp_mat* rawPtr = tempSparse.get();
+
+        // Move unique_ptr to a class member to ensure it stays alive
+        this->tempE = std::move(tempSparse);
+
+        // return rawPtr;
+    }
+
+    template<>
+    inline void INMF<H5SpMat>::load_whole_E(arma::uword i) {
+        H5SpMat* Eptr = this->Ei[i].get();
+        // arma::uword m = Eptr->n_rows;
+        arma::uword n = Eptr->n_cols;
+        this->tempE.reset();
+        arma::sp_mat tempSparse = Eptr->cols(0, n - 1);
+        auto tempSparseUniquePtr = std::make_unique<arma::sp_mat>(tempSparse);
+        // auto tempSparse = std::make_unique<arma::sp_mat>(m, n);
+        // (*tempSparse) = Eptr->cols(0, n - 1);
+        // Get raw pointer before transferring ownership
+        // arma::sp_mat* rawPtr = tempSparseUniquePtr.get();
+
+        // Move unique_ptr to a class member to ensure it stays alive
+        this->tempE = std::move(tempSparseUniquePtr);
+        // return rawPtr;
+    }
+
+    template<>
+    inline void INMF<H5Mat>::load_whole_E(arma::uword i) {
+        H5Mat* Eptr = this->Ei[i].get();
+        // arma::uword m = Eptr->n_rows;
+        arma::uword n = Eptr->n_cols;
+        auto out = std::make_unique<arma::sp_mat>(m, n);
+        // Load on-disk dense matrix by chunks and convert to sparse and fill
+        int numChunks = n / this->INMF_CHUNK_SIZE;
+        if (numChunks * this->INMF_CHUNK_SIZE < n) numChunks++;
+        for (int i = 0; i < numChunks; ++i) {
+            int spanStart = i * this->INMF_CHUNK_SIZE;
+            int spanEnd = (i + 1) * this->INMF_CHUNK_SIZE - 1;
+            if (spanEnd > n - 1) spanEnd = n - 1;
+            arma::mat dense_span = Eptr->cols(spanStart, spanEnd);
+            arma::sp_mat sparse_span(dense_span);
+            (*out).cols(spanStart, spanEnd) = sparse_span;
+        }
+        // Get raw pointer before transferring ownership
+        // arma::sp_mat* rawPtr = out.get();
+        // Move unique_ptr to a class member to ensure it stays alive
+        this->tempE = std::move(out);
+        // return rawPtr;
+    }
+
 }