2013-11-01 22 views
9

考慮下面的數據框。我想比較每行與下面的行,然後採取超過3個值相等的行。更快速地比較數據幀中的行

我寫了下面的代碼,但是如果你有一個大的數據框,它會很慢。

我怎麼能做得更快?

data <- as.data.frame(matrix(c(10,11,10,13,9,10,11,10,14,9,10,10,8,12,9,10,11,10,13,9,13,13,10,13,9), nrow=5, byrow=T)) 
rownames(data)<-c("sample_1","sample_2","sample_3","sample_4","sample_5") 

>data 
      V1 V2 V3 V4 V5 
sample_1 10 11 10 13 9 
sample_2 10 11 10 14 9 
sample_3 10 10 8 12 9 
sample_4 10 11 10 13 9 
sample_5 13 13 10 13 9 

output <- data.frame(sample = NA, duplicate = NA, matches = NA) 
dfrow <- 1 
for(i in 1:nrow(data)) { 
    sample <- data[i, ] 
    for(j in (i+1):nrow(data)) if(i+1 <= nrow(data)) { 
    matches <- 0 
     for(V in 1:ncol(data)) { 
      if(data[j,V] == sample[,V]) {  
       matches <- matches + 1 
      } 
     } 
     if(matches > 3) { 
      duplicate <- data[j, ] 
      pair <- cbind(rownames(sample), rownames(duplicate), matches) 
      output[dfrow, ] <- pair 
      dfrow <- dfrow + 1 
     } 
    } 
} 

>output 
    sample duplicate matches 
1 sample_1 sample_2  4 
2 sample_1 sample_4  5 
3 sample_2 sample_4  4 
+2

準確的數據集有多大?如果它不是很大,你可以交叉加入你的整個數據集並對其進行比較。此外,使用'data.table'而不是'data.frame'將有助於內存。 – TheComeOnMan

+0

250,000行乘26列 – vitor

+0

'data.table'不是行式敏感的。 – user974514

回答

8

這裏是一個RCPP的解決方案。但是,如果結果矩陣太大(即命中太多),則會引發錯誤。我運行循環兩次,首先得到結果矩陣的必要大小,然後填充它。可能有更好的可能性。另外,顯然,這隻會用整數。如果你的矩陣是數字的,你將不得不處理浮點精度。

library(Rcpp) 
library(inline) 

#C++ code: 
body <- ' 
const IntegerMatrix  M(as<IntegerMatrix>(MM)); 
const int     m=M.ncol(), n=M.nrow(); 
long      count1; 
int       count2; 
count1 = 0; 
for (int i=0; i<(n-1); i++) 
{ 
    for (int j=(i+1); j<n; j++) 
    { 
    count2 = 0; 
    for (int k=0; k<m; k++) { 
     if (M(i,k)==M(j,k)) count2++; 
    } 
    if (count2>3) count1++; 
    } 
} 
IntegerMatrix    R(count1,3); 
count1 = 0; 
for (int i=0; i<(n-1); i++) 
{ 
    for (int j=(i+1); j<n; j++) 
    { 
    count2 = 0; 
    for (int k=0; k<m; k++) { 
     if (M(i,k)==M(j,k)) count2++; 
    } 
    if (count2>3) { 
     count1++; 
     R(count1-1,0) = i+1; 
     R(count1-1,1) = j+1; 
     R(count1-1,2) = count2; 
    } 
    } 
} 
return wrap(R); 
' 

fun <- cxxfunction(signature(MM = "matrix"), 
        body,plugin="Rcpp") 

#with your data 
fun(as.matrix(data)) 
#  [,1] [,2] [,3] 
# [1,] 1 2 4 
# [2,] 1 4 5 
# [3,] 2 4 4 

#Benchmarks 
set.seed(42) 
mat1 <- matrix(sample(1:10,250*26,TRUE),ncol=26) 
mat2 <- matrix(sample(1:10,2500*26,TRUE),ncol=26) 
mat3 <- matrix(sample(1:10,10000*26,TRUE),ncol=26) 
mat4 <- matrix(sample(1:10,25000*26,TRUE),ncol=26) 
library(microbenchmark) 
microbenchmark(
    fun(mat1), 
    fun(mat2), 
    fun(mat3), 
    fun(mat4), 
    times=3 
) 
# Unit: milliseconds 
#  expr   min   lq  median   uq   max neval 
# fun(mat1)  2.675568  2.689586  2.703603  2.732487  2.761371  3 
# fun(mat2) 272.600480 274.680815 276.761151 276.796217 276.831282  3 
# fun(mat3) 4623.875203 4643.634249 4663.393296 4708.067638 4752.741979  3 
# fun(mat4) 29041.878164 29047.151348 29052.424532 29235.839275 29419.254017  3 
+0

+1好用的rcpp –

+0

謝謝!我不用C語言編寫代碼,這將是學習一些很好的機會......是的,我的真實矩陣是數字。矩陣中的值是離散類別,我有一些類別,如10.1或9.3。這不是精確度問題,這些數字必須保持原樣。 – vitor

1

這不是一個完整的答案,只是談到一記快速的鍛鍊是使用矩陣,而不是data.frame(這些都是相當慢TBH)。矩陣在R中非常快,並且通過完成至少一些操作,然後向列添加矢量將導致顯着的速度增加。

只是一個快速演示:

data <- matrix(c(10,11,10,13,9,10,11,10,14,9,10,10,8,12,9,10,11,10,13,9,13,13,10,13,9), nrow=5, byrow=T)rownames(data)<-c("sample_1","sample_2","sample_3","sample_4","sample_5") 
mu<-c("sample_1","sample_2","sample_3","sample_4","sample_5") 

t=proc.time() 
tab <- data.frame(sample = NA, duplicate = NA, matches = NA) 
dfrow <- 1 
for(i in 1:nrow(data)) { 
    sample <- data[i, ] 
    for(j in (i+1):nrow(data)) if(i+1 <= nrow(data)) { 
    matches <- 0 
     for(V in 1:ncol(data)) { 
      if(data[j,V] == sample[V]) {  
       matches <- matches + 1 
      } 
     } 
     if(matches > 3) { 
      duplicate <- data[j, ] 
      pair <- cbind(mu[i], mu[j], matches) 
      tab[dfrow, ] <- pair 
      dfrow <- dfrow + 1 
     } 
    } 
} 
proc.time()-t 

平均而言,我的機器上,產量

user system elapsed 
    0.00 0.06 0.06 

雖然你的情況,我得到

user system elapsed 
    0.02 0.06 0.08 

我不知道是否有比矩陣更快的東西。你也可以玩並行化,但對於循環C++代碼內聯經常使用(包Rcpp)。

+0

Catch 22 - 在「tab」中添加行可能對於大型數據集而言效率低下,而預分配一個非常大的數據集以在每次迭代中進行更新可能效率不高。任何其他想法? – TheComeOnMan

0

在評論中所說的一切都非常有效;特別是,我也不一定認爲R是做這件事的最好的地方。這就是說,這個工程快了很多,我比你帶來一個更大的數據集的內容(〜9.7秒對未完成兩分鐘後):

data <- matrix(sample(1:30, 10000, replace=TRUE), ncol=5) 
#Pre-prepare 
x <- 1 
#Loop 
for(i in seq(nrow(data)-2)){ 
    #Find the number of matches on that row 
    sums <- apply(data[seq(from=-1,to=-i),], 1, function(x) sum(x==data[i,])) 
    #Find how many are greater than/equal to 3 
    matches <- which(sums >= 3) 
    #Prepare output 
    output[seq(from=x, length.out=length(matches)),1] <- rep(i, length(matches)) 
    output[seq(from=x, length.out=length(matches)),2] <- matches 
    output[seq(from=x, length.out=length(matches)),3] <- sums[matches] 
    #Alter the counter of how many we've made... 
    x <- x + length(matches) 
} 
#Cleanup output 
output <- output[!is.na(output[,1]),]}) 

...我相當肯定我的怪異x變量和output的賦值可能會改善/變成apply-類型的問題,但它已經很晚了,我很累!祝你好運!

1
library(data.table) 

#creating the data 
dt <- data.table(read.table(textConnection(
"Sample   V1 V2 V3 V4 V5 
sample_1 10 11 10 13 9 
sample_2 10 11 10 14 9 
sample_3 10 10 8 12 9 
sample_4 10 11 10 13 9 
sample_5 13 13 10 13 9"), header= TRUE)) 

# some constants which will be used frequently 
nr = nrow(dt) 
nc = ncol(dt)-1 

#list into which we will insert the no. of matches for each sample 
#for example's sake, i still suggest you write output to a file possibly 
totalmatches <- vector(mode = "list", length = (nr-1)) 

#looping over each sample 
for (i in 1:(nr-1)) 
{ 
    # all combinations of i with i+1 to nr 
    samplematch <- cbind(dt[i],dt[(i+1):nr]) 

    # renaming the comparison sample columns 
    setnames(samplematch,append(colnames(dt),paste0(colnames(dt),"2"))) 

    #calculating number of matches 
    samplematch[,noofmatches := 0] 
    for (j in 1:nc) 
    { 
     samplematch[,noofmatches := noofmatches+1*(get(paste0("V",j)) == get(paste0("V",j,"2")))] 
    } 

    # removing individual value columns and matches < 3 
    samplematch <- samplematch[noofmatches >= 3,list(Sample,Sample2,noofmatches)] 

    # adding to the list 
    totalmatches[[i]] <- samplematch 
} 

輸出 -

rbindlist(totalmatches) 
    Sample Sample2 noofmatches 
1: sample_1 sample_2   4 
2: sample_1 sample_4   5 
3: sample_1 sample_5   3 
4: sample_2 sample_4   4 
5: sample_4 sample_5   3 

上矩陣的表現似乎更好,雖然,這種方法主頻 -

user system elapsed 
    0.17 0.01 0.19 
+0

正在使用'data.table'實際上在這裏做什麼?如果你只是逐行循環,我不確定使用它有什麼好處。 – David

+0

@大衛,公平點。我正在玩耍,試圖看看將樣本設置爲關鍵點有助於加速循環,無論是簡單查找行號會更快。你知道嗎?還使用了'rbindlist'。 – TheComeOnMan

+0

但是,無論情況如何,我通常在data.frame上使用'data.table's。沒什麼可失去的。 – TheComeOnMan

3

編輯:不知道我在想昨天晚上,當我減去行考慮我可以直接測試的平等。從下面的代碼中刪除了不必要的步驟。

這裏有一種方法,可能會稍微聰明或思路不清......但希望前者。這個想法是,不是逐行地進行一系列比較,而是可以通過從數據框的其餘部分中減去該行然後查看等於零的元素的數量來執行一些向量化的操作。下面是一個簡單的實施辦法:

> library(data.table) 
> data <- as.data.frame(matrix(c(10,11,10,13,9,10,11,10,14,9,10,10,8,12,9,10,11,10,13,9,13,13,10,13,9), nrow=5, byrow=T)) 
> rownames(data)<-c("sample_1","sample_2","sample_3","sample_4","sample_5") 
> 
> findMatch <- function(i,n){ 
+ tmp <- colSums(t(data[-(1:i),]) == unlist(data[i,])) 
+ tmp <- tmp[tmp > n] 
+ if(length(tmp) > 0) return(data.table(sample=rownames(data)[i],duplicate=names(tmp),match=tmp)) 
+ return(NULL) 
+ } 
> 
> system.time(tab <- rbindlist(lapply(1:(nrow(data)-1),findMatch,n=3))) 
    user system elapsed 
    0.003 0.000 0.003 
> tab 
    sample duplicate match 
1: sample_1 sample_2  4 
2: sample_1 sample_4  5 
3: sample_2 sample_4  4 

編輯:下面是一個使用矩陣和預tranposes版本2中的數據,所以你只需要做到這一點一次。它應該使用不重要的數據量更好地擴展到您的示例。

library(data.table) 
data <- matrix(round(runif(26*250000,0,25)),ncol=26) 
tdata <- t(data) 

findMatch <- function(i,n){ 
    tmp <- colSums(tdata[,-(1:i)] == data[i,]) 
    j <- which(tmp > n) 
    if(length(tmp) > 0) return(data.table(sample=i,duplicate=j+1,match=tmp[j])) 
    return(NULL) 
} 

tab <- rbindlist(lapply(1:(nrow(data)-1),findMatch,n=3)) 

我跑比我的機位上,並通過第一1500次迭代在15分鐘內得到了一個完整的250000×26矩陣,需要600 MB的內存。由於以前的迭代不會影響未來的迭代,因此如果需要,您可以將其分塊並分別運行。

+0

我喜歡這個,但它真正的數據集運行了幾個小時後崩潰了我的4GB內存計算機,這是25萬行26列。 – vitor

+0

這並不奇怪,這裏有很多計算,R對內存來說很差,而且你沒有太多內存可供使用。你知道在哪些迭代中失控嗎?如果是這樣的話,你可以嘗試在它之前對它進行分塊,但是你最好使用像數據庫或Rcpp這樣的非R方法。 – David

+0

添加了一個稍微更新的版本,應該爲你工作,但我認爲那是更好的選擇。 – David

0

那麼,我採取了刺戳它,下面的代碼運行速度比原來快大約3倍。

f <- function(ind, mydf){ 
    res <- NULL 
    matches <- colSums(t(mydf[-(1:ind),])==mydf[ind,]) 
    Ndups <- sum(matches > 3) 
    if(Ndups > 0){ 
     res <- data.frame(sample=rep(ind,Ndups),duplicate=which(matches > 3), 
         matches= matches[matches > 3],stringsAsFactors = F) 
     rownames(res) <- NULL 
     return(as.matrix(res)) 
    } 
    return(res) 
} 


f(1,mydf=as.matrix(data)) 
f(2,mydf=as.matrix(data)) 
system.time( 
for(i in 1:1000){ 
    tab <- NULL 
    for(j in 1:(dim(data)[1]-1)) 
     tab <- rbind(tab,f(j,mydf=as.matrix(data))) 
} 
)/1000 
tab 
0

假設數據集中的所有條目都是相同的模式(數字),將它轉換爲矩陣。通過轉置,您可以利用==的矢量化方式。

data <- as.matrix(data) 
data <- t(data) 

output <- lapply(seq_len(ncol(data) - 1), function(x) { 
    tmp <- data[,x] == data[, (x+1):ncol(data)] 
    n_matches <- { 
     if (x == ncol(data) - 1) { 
      setNames(sum(tmp),colnames(data)[ncol(data)]) 
     } else { 
      colSums(tmp) 
     } 
    } 
    good_matches <- n_matches[n_matches >= 3] 
}) 

最大的問題是如何輸出結果。現在,我將你的數據列入清單。我認爲這是存儲數據最少的內存密集型方式。

[[1]] 
sample_2 sample_4 sample_5 
     4  5  3 

[[2]] 
sample_4 
     4 

[[3]] 
named numeric(0) 

[[4]] 
sample_5 
     3 

如果你想有一個數據幀輸出,那麼你會希望內lapply調整函數的返回值。也許添加在函數的最後一行:

return(data.frame(
    sample = colnames(data)[x], 
    duplicate = names(good_matches), 
    noofmatches = good_matches, 
    stringsAsFactors = FALSE)) 

然後用:

newoutput <- do.call(rbind, output) 
## or, using plyr 
# require(plyr) 
# newoutput <- rbind.fill(output)