przepisywanie powolne funkcję R C++ i RCPP

Question

I tej linii kodu r:

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

co robi to identyfikacji miejsc (cols) w matrycy z sekwencjami DNA (1 wiersz = z SEK), który nie są uniwersalne (informacyjne) i dzielą je z macierzy, aby utworzyć nową "matrycę przyciętą", tj. pozbyć się wszystkich kolumn, w których wartości są takie same. W przypadku dużego zbioru danych zajmuje to około 6 sekund. Nie wiem, czy mogę to zrobić szybciej w C++ (wciąż początkującym w C++), ale dobrze będzie, jeśli spróbuję. Moim pomysłem jest użycie Rcpp, zapętlenie w kolumnach CharacterMatrix, wyciągnięcie kolumny (witryny) jako sprawdzenia CharacterVector, jeśli są one takie same. Jeśli są takie same, zapisz numer kolumny/indeks, kontynuuj dla wszystkich kolumn. Następnie na końcu utwórz nowy CharacterMatrix, który obejmuje tylko te kolumny. Ważne jest, aby zachować rownames i nazwy kolumn tak, jak są one w "wersji R" macierzy, tj. Jeśli kolumna idzie, tak samo jak nazwa kolumny.

Pisałem przez około dwie minuty, do tej pory, co mam jest (nie skończony):

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

Czy będę właściwą drogę na ten temat? Czy istnieje prostszy sposób. Rozumiem, że potrzebuję std :: vector, ponieważ łatwo jest je rozwijać (ponieważ nie wiem z góry ile kolczyków zamierzam zachować). Przy indeksowaniu będę musiał dać +1 do wektora informacyjnegoSites na końcu (ponieważ R indeksuje od 1 i C++ od 0)?

dzięki, ben W.

Answer 1

dane próbki:

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

roztwór OP:

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 

szybsza wersja z zastosowaniem zasady R:

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

Wyniki są ide ntical:

identical(y1, y2) 
# [1] TRUE 

To bardzo możliwe, że C++ go pokona, ale czy to naprawdę konieczne?