Это класс разреженной матрицы (dgCMatrix), расширяющий Rcpp.
КАКИЕ: Этот класс включает Rcpp::NumericVector
а также Rcpp::NumericMatrix
классы, которые являются простыми ссылками на объекты R. Этот класс полностью поддерживает преобразование без копирования только для ссылок. Matrix
упаковка dgCMatrix
класс от R до C ++ и наоборот.
ЗАЧЕМ: Мне нужен доступ только для чтения без копирования к >R Matrix::dgCMatrix
в double
введите C ++, в отличие от RcppArmadillo SpMat
и RcppEigen SparseMatrix
которые являются глубокими копиями. Я определенно не пытаюсь заменить отличную линейную алгебру и поэлементные операции, которые RcppArmadillo
а также RcppEigen
уже предлагают.
Особенно меня не устраивает синтаксис итератора. Все работает, но может работать не так, как ожидалось (я новичок в C ++ / Rcpp несколько месяцев назад).
Полный код здесь.
RcppSparse.h
#include <rcpp.h>
namespace Rcpp {
class dgCMatrix {
public:
IntegerVector i, p, Dim;
NumericVector x;
List Dimnames;
// constructors
dgCMatrix(IntegerVector& A_i, IntegerVector& A_p, NumericVector& A_x, int nrow) {
i = A_i;
p = A_p;
x = A_x;
Dim = IntegerVector::create(nrow, A_p.size() - 1);
};
dgCMatrix(IntegerVector& A_i, IntegerVector& A_p, NumericVector& A_x, int nrow, List& A_Dimnames) {
i = A_i;
p = A_p;
x = A_x;
Dim = IntegerVector::create(nrow, A_p.size() - 1);
Dimnames = A_Dimnames;
};
dgCMatrix(S4 mat) {
i = mat.slot("i");
p = mat.slot("p");
x = mat.slot("x");
Dim = mat.slot("Dim");
Dimnames = mat.slot("Dimnames");
};
// basic properties
int nrow() { return Dim[0]; };
int ncol() { return Dim[1]; };
int rows() { return Dim[0]; };
int cols() { return Dim[1]; };
int n_nonzero() { return x.size(); };
NumericVector& nonzeros() { return x; };
double sum() { return Rcpp::sum(x); };
// forward constant iterator
class const_iterator {
public:
int index;
const_iterator(dgCMatrix& g, int ind) : parent(g) { index = ind; }
bool operator!=(const_iterator x) const { return index != x.index; };
const_iterator& operator++(int) { ++index; return (*this); };
int row() { return parent.i[index]; };
int col() { int j = 0; for (; j < parent.p.size(); ++j) if (parent.p[j] >= index) break; return j; };
double operator*() const { return parent.x[index]; };
private:
dgCMatrix& parent;
};
// iterator constructors
const_iterator begin(int j) { return const_iterator(*this, (int)0); };
const_iterator end(int j) { return const_iterator(*this, i.size() - 1); };
const_iterator begin_col(int j) { return const_iterator(*this, p[j]); };
const_iterator end_col(int j) { return const_iterator(*this, p[j + 1]); };
// read-only element access
double at(int row, int col) const {
for (int j = p[col]; j < p[col + 1]; ++j) {
if (i[j] == row) return x[j];
else if (i[j] > row) break;
}
return 0.0;
}
double operator()(int row, int col) { return at(row, col); };
NumericVector operator()(int row, IntegerVector& col) {
NumericVector res(col.size());
for (int j = 0; j < col.size(); ++j) res[j] = at(row, col[j]);
return res;
};
NumericVector operator()(IntegerVector& row, int col) {
NumericVector res(row.size());
for (int j = 0; j < row.size(); ++j) res[j] = at(row[j], col);
return res;
};
NumericMatrix operator()(IntegerVector& row, IntegerVector& col) {
NumericMatrix res(row.size(), col.size());
for (int j = 0; j < row.size(); ++j)
for (int k = 0; k < col.size(); ++k)
res(j, k) = at(row[j], col[k]);
return res;
};
// column access (copy)
NumericVector col(int col) {
NumericVector c(Dim[0], 0.0);
for (int j = p[col]; j < p[col + 1]; ++j)
c[i[j]] = x[j];
return c;
}
NumericVector column(int c) { return col(c); }
NumericMatrix cols(IntegerVector& c) {
NumericMatrix res(Dim[0], c.size());
for (int j = 0; j < c.size(); ++j) {
res.column(j) = col(c[j]);
}
return res;
}
NumericMatrix columns(IntegerVector& c) { return cols(c); }
// row access (copy)
NumericVector row(int row) {
NumericVector r(Dim[1], 0.0);
for (int col = 0; col < Dim[1]; ++col) {
for (int j = p[col]; j < p[col + 1]; ++j) {
if (i[j] == row) r[col] = x[j];
else if (i[j] > row) break;
}
}
return r;
}
NumericMatrix rows(IntegerVector& r) {
NumericMatrix res(r.size(), Dim[1]);
for (int j = 0; j < r.size(); ++j) {
res.row(j) = row(r[j]);
}
return res;
}
// colSums and rowSums family
NumericVector colSums() {
NumericVector sums(Dim[1]);
for (int col = 0; col < Dim[1]; ++col)
for (int j = p[col]; j < p[col + 1]; ++j)
sums(col) += x[j];
return sums;
}
NumericVector rowSums() {
NumericVector sums(Dim[0]);
for (int col = 0; col < Dim[1]; ++col)
for (int j = p[col]; j < p[col + 1]; ++j)
sums(i[j]) += x[j];
return sums;
}
NumericVector colMeans() {
NumericVector sums = colSums();
for (int i = 0; i < sums.size(); ++i) sums[i] = sums[i] / Dim[0];
return sums;
};
NumericVector rowMeans() {
NumericVector sums = rowSums();
for (int i = 0; i < sums.size(); ++i) sums[i] = sums[i] / Dim[1];
return sums;
};
};
// Rcpp::as
template <> dgCMatrix as(SEXP mat) { return dgCMatrix(mat); }
// Rcpp::wrap
template <> SEXP wrap(const dgCMatrix& sm) {
S4 s(std::string("dgCMatrix"));
s.slot("i") = sm.i;
s.slot("p") = sm.p;
s.slot("x") = sm.x;
s.slot("Dim") = sm.Dim;
s.slot("Dimnames") = sm.Dimnames;
return s;
}
}
ПРИМЕР ИСПОЛЬЗОВАНИЯ:
В C ++:
#include <RcppSparse.h>
// [[Rcpp::export]]
Rcpp::NumericVector Rcpp_colSums(Rcpp::dgCMatrix& mat){
return mat.colSums();
}
В R:
library(Matrix)
mat <- rsparsematrix(100, 100, 0.5)
colSums <- Rcpp_colSums(mat)