csc_to_csr.cpp

#include <Rcpp.h>
using namespace Rcpp;

// This is a simple example of exporting a C++ function to R. You can
// source this function into an R session using the Rcpp::sourceCpp 
// function (or via the Source button on the editor toolbar). Learn
// more about Rcpp at:
//
//   http://www.rcpp.org/
//   http://adv-r.had.co.nz/Rcpp.html
//   http://gallery.rcpp.org/
//

typedef uint64_t I;
typedef float T;

// [[Rcpp::export]]
SEXP csc_tocsr(
    const I n_row,
    const I n_col, 
    const std::vector<I> &Ap, 
    const std::vector<I> &Ai, 
    const std::vector<T> &Ax
) {
  // source code in-parts taken from 
  // https://github.com/scipy/scipy/blob/3b36a574dc657d1ca116f6e230be694f3de31afc/scipy/sparse/sparsetools/csr.h#L379-L424
  const I nnz = Ap[n_col];
  std::vector<I> Bp(n_row+1), Bj(nnz);
  std::vector<T> Bx(nnz);
  for (I n = 0; n < nnz; n++){            
    Bp[Ai[n]]++;
  }
  for(I row = 0, cumsum = 0; row < n_row; row++){     
    I temp  = Bp[row];
    Bp[row] = cumsum;
    cumsum += temp;
  }
  Bp[n_row] = nnz; 
  for(I col = 0; col < n_col; col++){
    for(I jj = Ap[col]; jj < Ap[col+1]; jj++){
      I row  = Ai[jj];
      I dest = Bp[row];
      Bj[dest] = col;
      Bx[dest] = Ax[jj];
      Bp[row]++;
    }
  }
  for(I row = 0, last = 0; row <= n_row; row++){
    I temp  = Bp[row];
    Bp[row] = last;
    last    = temp;
  }
  
  Rcpp::List L = Rcpp::List::create(Rcpp::Named("p") = Bp , Rcpp::Named("i") = Bj, Rcpp::Named("x") = Bx);
  return L;
}

// You can include R code blocks in C++ files processed with sourceCpp
// (useful for testing and development). The R code will be automatically 
// run after the compilation.
//

/*** R
library(SeuratObject)
data("pbmc_small")

m <- GetAssayData(pbmc_small, "counts")
m2.list <- csc_tocsr(m@Dim[[1]], m@Dim[[2]], m@p, m@i, m@x)
m2 <- new(
  Class = 'dgRMatrix',
  p = as.integer(m2.list$p),
  j = as.integer(m2.list$i),
  x = m2.list$x,
  Dim = slot(object = m, name = 'Dim')
)

mr <- as(object = as.matrix(x = m), Class = 'dgRMatrix')

all.equal(target = slot(object = m2, name = 'p'), current = slot(object = mr, name = 'p'))
all.equal(target = slot(object = m2, name = 'j'), current = slot(object = mr, name = 'j'))
all.equal(target = slot(object = m2, name = 'x'), current = slot(object = mr, name = 'x'))
all.equal(target = slot(object = m2, name = 'Dim'), current = slot(object = mr, name = 'Dim'))
*/