Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Rewriting slow R function in C++ & Rcpp

Tags:

c++

r

vector

rcpp

I have this line of R code:

croppedDNA <- completeDNA[,apply(completeDNA,2,function(x) any(c(FALSE,x[-length(x)]!=x[-1])))]

What it does is identify the sites (cols) in a matrix of DNA sequences (1 row = one seq) that are not universal (informative) and subsets them from the matrix to make a new 'cropped matrix' i.e. get rid of all the columns in which values are the same. For a big dataset this takes about 6 seconds. I don't know if I can do it faster in C++ (still a beginner in C++) but it will be good for me to try. My idea is to use Rcpp, loop through the columns of the CharacterMatrix, pull out the column (the site) as a CharacterVector check if they are the same. If they are the same, record that column number/index, continue for all columns. Then at the end make a new CharacterMatrix that only includes those columns. It is important that I keep the rownames and column names as they are in th "R version" of the matrix i.e. if a column goes, so should the colname.

I've been writing for about two minutes, so far what I have is (not finished):

#include <Rcpp.h>
#include <vector>
using namespace Rcpp;
// [[Rcpp::export]]
CharacterMatrix reduce_sequences(CharacterMatrix completeDNA)
{
  std::vector<bool> informativeSites; 
  for(int i = 0; i < completeDNA.ncol(); i++)
  {
    CharacterVector bpsite = completeDNA(,i);
    if(all(bpsite == bpsite[1])
    {
      informativeSites.push_back(i);
    }
  }
CharacterMatrix cutDNA = completeDNA(,informativeSites);
return cutDNA;
}

Am I going the right way about this? Is there an easier way. My understanding is I need std::vector because it's easy to grow them (since I don't know in advance how many cols I am going to want to keep). With the indexing will I need to +1 to the informativeSites vector at the end (because R indexes from 1 and C++ from 0)?

Thanks, Ben W.

like image 557
Ward9250 Avatar asked May 15 '13 02:05

Ward9250


People also ask

How fast is R compared to C?

R runs between 475 to 491 times slower than C++. If the code is compiled, the code is between 243 to 282 times slower. Hybrid programming and special approaches can deliver considerable speed ups.


1 Answers

Sample data:

set.seed(123)
z <- matrix(sample(c("a", "t", "c", "g", "N", "-"), 3*398508, TRUE), 3, 398508)

OP's solution:

system.time(y1 <- z[,apply(z,2,function(x) any(c(FALSE,x[-length(x)]!=x[-1])))])
#    user  system elapsed 
#   4.929   0.043   4.976 

A faster version using base R:

system.time(y2 <- (z[, colSums(z[-1,] != z[-nrow(z), ]) > 0]))
#    user  system elapsed 
#   0.087   0.011   0.098 

The results are identical:

identical(y1, y2)
# [1] TRUE

It's very possible c++ will beat it, but is it really necessary?

like image 100
flodel Avatar answered Sep 17 '22 23:09

flodel