Класс разреженной матрицы CSC Rcpp

Это класс разреженной матрицы (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)