2015-09-25 54 views
1

我有一個數組數據=陣列[1:50,1:50,1:50]數組R是值內是-1之間的實數,1優化循環使用並行

「數據」能視爲立方體50x50x50。

我需要創建基於此方程=>

值=(X + Y)的相關矩陣(除去全零) - | X-Y |並且矩陣大小是可能組合(50×50×50)×((50×50×50)-1)/ 2 = 7.812.437.500這2倍=相關矩陣的2倍。

我這樣做:

比方說我們的3x3x3:

arr = array(rnorm(10), dim=c(3,3,3)) 

data = data.frame(array(arr)) 


data$voxel <- rownames(data) 

#remove zeros 
data<-data[!(data[,1]==0),] 

rownames(data) = data$voxel 

data$voxel = NULL 


####################################################################################### 
#Create cluster 

no_cores <- detectCores() #- 1 

clus <- makeCluster(no_cores) 

clusterExport(clus, list("data") , envir=environment()) 

clusterEvalQ(clus, 
      compare_strings <- function(j,i) { 
       value <- (data[i,]+data[j,])-abs(data[i,]- data[j,]) 
       pair <- rbind(rownames(data)[j],rownames(data)[i],value) 
       return(pair) 
      }) 

i = 0 # start 0 
kk = 1 
table <- data.frame() 

ptm <- proc.time() 

while(kk<nrow(data)) { 

    out <-NULL 
    i = i+1 # fix row 
    j = c((kk+1):nrow(data)) # rows to be compared 

    #Apply the declared function 
    out = matrix(unlist(parRapply(clus,expand.grid(i,j), function(x,y) compare_strings(x[1],x[2]))),ncol=3, byrow = T) 

    table <- rbind(table,out) 

    kk = kk +1 

} 

proc.time() - ptm 

結果是data.frame:

v1 v2 v3 
1 2 2.70430114250358 
1 3 0.199941717684129 
... up to 351 rows 

但是這將需要數天...

另外,我想創建一個這種關聯矩陣:

1       2    3... 
1 1     2.70430114250358 
2 2.70430114250358   1 
3... 

有沒有更快的方法來做到這一點?

感謝

+3

請給我們一個小[再現的示例](http://stackoverflow.com/a/5963610/1412059)(例如,用3x3x3的陣列)與和顯示工作預期的產出。如果無法找到矢量化解決方案(可疑),則應使用Rcpp執行此操作(即,在編譯代碼中執行循環)。 – Roland

+0

由於無法找到「S」,因此您當前生成'data'的代碼無法運行。 – Heroka

+0

大家好,我已經編輯了一些更多解釋的帖子。謝謝 – DemetriusRPaula

回答

0

有一些在你的代碼性能的錯誤:

  1. 你循環時,你應該依靠量化。
  2. 你在循環中生長一個對象。
  3. 您可以並行化循環的每個迭代而不是並行化外循環。

如果避免第一個問題,可以避免所有這些問題。

顯然,你想要比較每個行的組合。對於這一點,你應該先把排索引的所有組合:

combs <- t(combn(1:27, 2)) 

那麼你可以申請比較函數這些:

compare <- function(j,i, data) { 
    as.vector((data[i,]+data[j,])-abs(data[i,]- data[j,])) 
} 

res <- data.frame(V1 = combs[,1], V2 = combs[,2], 
        V3 = compare(combs[,1], combs[,2], data)) 

現在,如果我們要檢查,如果這給出結果爲相同你的代碼,我們首先需要修復你的輸出。通過將字符(rownames)與矩陣中的數字相結合,可以得到一個字符矩陣,並且最終data.frame的列都是字符。我們可以用type.convert來修復之後(儘管它應該從一開始就避免):

table[] <- lapply(table, function(x) type.convert(as.character(x))) 

現在我們看到的結果是一樣的:

all.equal(res, table) 
#[1] TRUE 

如果你喜歡,你可以把結果爲稀疏矩陣:

library(Matrix) 
m <- sparseMatrix(i = res$V1, j = res$V2, x = res$V3, 
        dims = c(27, 27), symmetric = TRUE) 
diag(m) <- 1 
+0

combs <-t(combn(1:83346,2))不適用於大小:( – DemetriusRPaula

+0

)那麼這就是'3,473,236,185'組合。我相信你應該重新考慮你想要做的事情,但是如果你堅持要做到這一點,你可以使用Rcpp。當然,你需要一個大的RAM,或者將Rcpp與其中一個包裝用於內存不足的數據結合。 – Roland

+0

cppFunction('Rcpp :: DataFrame combi2inds(const Rcpp :: CharacterVector inputVector)const int len = inputVector.size(); const int retLen = len *(len-1)/ 2; Rcpp :: IntegerVector outputVector1(retLen); Rcpp :: IntegerVector outputVector2(retLen); int indexSkip; for(int i = 0; i DemetriusRPaula