The enclosed works on this example. The logic is to check each column in
the index set for all the elements of the index set, except the one on the
diagonal, being in the set of row indices. I have added some diagnostic
output to demonstrate the flow.
On Tue, Mar 25, 2014 at 9:48 AM, Douglas Bates <[email protected]> wrote:
> On Tue, Mar 25, 2014 at 6:42 AM, Søren Højsgaard <[email protected]>wrote:
>
>> Dear all,
>> I have a large sparse adjacency matrix X for an undirected graph. For a
>> subset 'idx' of the vertices I want to find out if the subgraph defined by
>> this subset is complete (i.e. has edges between all variables). So in R,
>> one would check if X[idx,idx] has only ones outside the diagonal, but this
>> is slow and I want to do it with Rcpp. Here is one attempt where I simply
>> add up the elements above (or below) the diagonal, and to access the
>> elements of X I use coeff which is relatively slow (because X is a sparse
>> matrix).
>>
>> #include <RcppEigen.h>
>> typedef Eigen::MappedSparseMatrix<double> MSpMat;
>> typedef Eigen::SparseVector<double> SpVec;
>> typedef SpVec::InnerIterator InIter;
>>
>> //[[Rcpp::export]]
>> double foo (MSpMat X, Eigen::VectorXi idx){
>> SpVec sidx = idx.sparseView();
>> double out=0;
>> int i2, j2;
>> for (InIter i_(sidx); i_; ++i_){
>> i2 = i_.index();
>> for (InIter j_(sidx); j_; ++j_){
>> j2 = j_.index();
>> if (i2>j2)
>> out += X.coeff( i2, j2);
>> }
>> }
>> return out;
>> }
>>
>> /*** R
>> library(Matrix)
>> M1 <- new("dgCMatrix"
>> , i = c(1L, 2L, 3L, 0L, 2L, 3L, 0L, 1L, 3L, 0L, 1L, 2L, 4L, 5L, 3L,
>> 5L, 3L, 4L)
>> , p = c(0L, 3L, 6L, 9L, 14L, 16L, 18L)
>> , Dim = c(6L, 6L)
>> , Dimnames = list(c("a", "b", "c", "d", "e", "f"), c("a", "b", "c",
>> "d", "e", "f"))
>> , x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
>> , factors = list()
>> )
>> M1
>>
>> foo(M1, c(2,3,4))
>>
>> */
>> Can anyone see a better way of doing this?
>>
>> I was thinking about whether the sparse matrix iterators could be a
>> solution, but I can't a hold on it.
>
>
> Sparse matrix inner and outer iterators are definitely the way to go.
> I'll write a version using them.
>
>
>> I was also thinking about using sparse matrices in RcppArmadillo, but I
>> guess the problem (speed) is the same (haven't tried!).
>>
>
>
>> Any advice will be greatly appreciated.
>> Cheers
>> Søren
>>
>>
>>
>> _______________________________________________
>> Rcpp-devel mailing list
>> [email protected]
>> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>>
>>
// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
typedef Eigen::MappedSparseMatrix<double> MSpMat;
typedef MSpMat::InnerIterator InIter;
//[[Rcpp::export]]
bool foo (MSpMat X, Eigen::VectorXi idx){
int n = X.cols();
if (X.rows() != n) throw std::invalid_argument("Sparse matrix X must be square");
for (int i = 0; i < idx.size(); ++i) {
int jj = idx[i] - 1; // 0-based column index
InIter it(X,jj);
Rcpp::Rcout << "i = " << i << ", jj = " << jj << std::endl;
for (int k = 0; k < idx.size(); ++k) {
int kk = idx[k]-1; // 0-based row index to search for
if (kk == jj) continue; // skip positions on the diagonal
Rcpp::Rcout << " kk = " << kk << ", it.row =";
bool foundit = false;
for (; it; ++it) {
Rcpp::Rcout << " " << it.row();
if (it.row() == kk) {
foundit = true;
++it;
break;
}
if (it.row() > kk) return false;
}
if (!foundit) return false;
Rcpp::Rcout << std::endl;
}
}
return true;
}
/*** R
library(Matrix)
M1 <- new("dgCMatrix"
, i = c(1L, 2L, 3L, 0L, 2L, 3L, 0L, 1L, 3L, 0L, 1L, 2L, 4L, 5L, 3L, 5L, 3L, 4L)
, p = c(0L, 3L, 6L, 9L, 14L, 16L, 18L)
, Dim = c(6L, 6L)
, Dimnames = list(c("a", "b", "c", "d", "e", "f"), c("a", "b", "c", "d", "e", "f"))
, x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
, factors = list()
)
M1
foo(M1, c(2,3,4))
***/
_______________________________________________
Rcpp-devel mailing list
[email protected]
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
