Viết lại chức năng R chậm trong C++ & Rcpp

Question

Tôi có dòng mã này R:

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

gì nó làm là xác định các trang web (cols) trong một ma trận của các chuỗi DNA (1 row = một seq) mà không phổ biến (thông tin) và tập hợp chúng từ ma trận để tạo ra một 'ma trận cắt' mới, tức là loại bỏ tất cả các cột trong đó các giá trị giống nhau. Đối với một tập dữ liệu lớn, mất khoảng 6 giây. Tôi không biết nếu tôi có thể làm điều đó nhanh hơn trong C + + (vẫn còn là một người mới bắt đầu trong C + +) nhưng nó sẽ là tốt cho tôi để thử. Ý tưởng của tôi là sử dụng Rcpp, lặp qua các cột của CharacterMatrix, kéo ra khỏi cột (trang web) như một kiểm tra CharacterVector nếu chúng giống nhau. Nếu chúng giống nhau, hãy ghi lại số cột/chỉ mục đó, tiếp tục cho tất cả các cột. Sau đó, ở cuối tạo một CharacterMatrix mới chỉ bao gồm các cột đó. Điều quan trọng là tôi giữ các tên gọi và tên cột khi chúng ở trong "phiên bản R" của ma trận tức là nếu một cột đi, vì vậy nên đặt tên cho cột.

Tôi đã viết trong khoảng hai phút, cho đến nay những gì tôi có là (chưa hoàn tất):

#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; 
} 

Tôi có đi đúng đường về điều này? Có cách nào dễ hơn không. Sự hiểu biết của tôi là tôi cần std :: vector bởi vì nó dễ dàng để phát triển chúng (vì tôi không biết trước bao nhiêu cols tôi sẽ muốn giữ). Với việc lập chỉ mục, tôi sẽ cần +1 đến vector thông tinSites ở cuối (vì R lập chỉ mục từ 1 và C++ từ 0)?

Cảm ơn, Ben W.

Answer 1

dữ liệu mẫu:

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

giải pháp OP của:

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 

Một phiên bản nhanh hơn sử dụng cơ sở R:

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

Kết quả là ide ntical:

identical(y1, y2) 
# [1] TRUE 

Rất có thể C++ sẽ đánh bại nó, nhưng thực sự có cần thiết không?