2014-04-04 128 views
0

這是我的R腳本,有三個嵌套的for循環。完成for循環的2000輪中的1輪需要2分多鐘。如何加速這一點?R:加速循環

col<-NULL 
row<-NULL 
rep<-ncol(dat)-2 
dist<-NULL 
c1=3 
for (i in 1:rep){ 
    c2=3 
    for(j in 1:rep){ 
    r=1 
    for (k in 1:nrow(dat)){ 
     p<-(dat[r,c1]-dat[r,c2])^2 
     row<-rbind(row, p) 
     r=r+1 
    } 
    row<-sqrt(sum(row)) 
    row1<-(1/(2*length(unique(dat[,1]))))*row 
    col<-cbind(col, row1) 
    c2=c2+1 
    row<-NULL 
    } 
    dist<-rbind(dist,col) 
    col<-NULL 
    c1=c1+1 
} 

編輯:

> head(dat) 
    mark alle G1 G2 G3 G4 G5 G6 G7 G8 G9 G10 G11 G12 G13 G14 G15 G16 G17 G18 G19 G20 G21 G22 G23 G24 
1 M1 228 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0.0 0.5 0 0 
2 M1 234 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.5 0.5 1 1 
3 M1 232 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0 0.0 0 0 
4 M1 240 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0 0.0 0 0 
5 M1 230 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.5 0.0 0 0 
6 M1 238 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0 0.0 0 0 
+0

如果你給出了一個'dat'數據的例子,以及對這個代碼打算做什麼的簡要描述,回答會更容易。 –

+0

@ Maxim.K:查看我編輯的'dat'文件。此代碼計算等位基因頻率數據的修正羅傑斯遺傳距離。 – ramesh

回答

4

我不知道修改羅傑斯的遺傳距離但它看起來像歐幾里得距離乘以1/(2*length(unique(dat$mark)))

f <- 1/(2*length(unique(dat$mark))) 
d <- f*dist(t(dat[, -c(1, 2)]), method="euclidean") 
+0

您可能想要添加d < - as.matrix(d)以將其轉換爲使用距離函數後更易於使用的對象 – Miff

+0

@sgibb:是的,您是正確的。非常感謝。你讓我的生活變得簡單。 – ramesh

+0

@ramesh:如果這解決了您的問題,請您接受答案。 – sgibb

3

可以做,以加快循環最重要的事情是循環之前預分配向量和矩陣。然後,而不是使用cbind()rbind(),結果增加了向量/矩陣,像這樣:

# Was: row<-rbind(row, p) 
row[k] <- p 

# Was: col<-cbind(col, row1) 
col[j] <- row1 

# Was: dist<-rbind(dist,col) 
dist[i, ] <- col 

之後,你可以探索的方式向量化操作,或者更好的,看是否有已經存在的功能來執行這個任務(或者如果任務是基於某個存在函數的東西的話)。此外,任何不依賴於循環的東西(例如row1<-(1/(2*length(unique(dat[,1])))))都應該移出循環。否則,你只是重複計算一遍又一遍的對性能產生負面影響的相同值。

與循環的關鍵是通過預先分配的向量和矩陣循環之前將提供的性能提升的很多避免rbind()cbind()

+0

他也可以將'row1 < - (1 /(2 * length(unique(dat [,1]))))'移出循環。 – sgibb

+0

是的,很好。我更新了我的答案以解釋這一點。 –

1

儘管類似的功能已經存在,我想我自己的路。
我刪除了一個完整的for循環,rbindcbind
現在這隻需要124秒就可以在1014 X 1014矩陣的一輪(表示1 X 1014)上寫入1014 X 1014矩陣2分鐘。

dat<-read.table("alreq1.txt", sep="\t",header=T) 
col<-NULL 
row<-NULL 
rep<-ncol(dat)-2 
dist<-NULL 
dist<- data.frame(matrix(NA, nrow = rep, ncol = rep)) 
m<-1/sqrt(2*length(unique(dat[,1]))) 
c1=3 
for (i in 1:rep){ 
    c2=3 
    for(j in 1:rep){ 
     p<-na.omit(dat[,c1]-dat[,c2])^2 
     row<-sum(p) 
     row<-sqrt(row)*m 
     col[j] <- row 
     c2=c2+1 
     row<-NULL 
     p<-NULL 
    } 
    dist[i,] <- col 
    c1=c1+1 
    col<-NULL 
    } 

希望這個代碼仍然可以改進。