2013-08-21 123 views
1

一些背景資料:我需要在3D網格從某點的距離計算出到每個小區,則應用一個函數這個距離。我需要爲多個點執行此操作,併爲所有點在每個單元格中添加函數值。我可以使用下面的代碼爲位於(x,Y,Z)的點做到這一點:řforeach循環三維陣列數據

x <- c(1,2,3,4,5) 
y <- x 
z <- x 
radius <- c(0.4,0.5,0.6,0.7,0.8) 
numsphere <- length(x) 
radius_buffer <- 0.2 

xvox <- seq((min(x)-1),(max(x)+2),0.5) 
yvox <- xvox 
zvox <- xvox 

probability_array <<- array(0,dim=c(length(xvox),length(yvox),length(zvox))) 

for (j in 1:length(yvox)){ # for every y element 
    for (i in 1:length(xvox)){ # for every x element 
    for (k in length(zvox):1){ # for every z element 
     for (n in 1:numsphere){ # for the total number of points 
     dist_sd <- ((xvox[i]-x[n])^2+(yvox[j]-y[n])^2+(zvox[k]-z[n])^2)^0.5 
     probability_array[i,j,k] <- probability_array[i,j,k] + 
            round(exp(-1*(dist_sd-radius[n])^2/(2*radius_buffer^2)),3) 
      } 
     } 
     } 
    } 

輸出是一個數組,所標繪的結果如下所示:

probability_array <- probability_array/max(probability_array) 
contour3d(probability_array,level=c(0.2,0.8,0.9),x=xvox,y=yvox,z=zvox,color = c("aquamarine","gold","darkorange"),alpha = c(0.1,0.2,0.5),add=T) 

Output graphic

我試圖平行這個,因爲它似乎是理想的,但不能得到它的工作。 我已經試過:

cl<-makeCluster(detectCores(),type="SOCK") 
registerDoSNOW(cl) 

for (j in 1:length(yvox)){ 
    for (i in 1:length(xvox)){ 
    for(k in length(zvox):1){ 
     probability_array[i,j,k] <- foreach(n=1:numsphere, .combine='+') %dopar% { 
     dist_sd <- ((xvox[i]-x[n])^2+(yvox[j]-y[n])^2+(zvox[k]-z[n])^2)^0.5 
     round(exp(-1*(dist_sd-radius[n])^2/(2*radius_buffer^2)),3) 
     } 
    } 
    } 
} 

之類的東西:

r <- foreach(j=1:length(yvox)) %:% foreach(i=1:length(xvox)) %:% foreach(k=length(zvox):1) %:% foreach(n=1:numsphere, .combine='+') %do% { 

     dist_sd <- ((xvox[i]-x[n])^2+(yvox[j]-y[n])^2+(zvox[k]-z[n])^2)^0.5 
     probability_array[i,j,k] <- probability_array[i,j,k] + round(exp(-1*(dist_sd-radius[n])^2/(2*radius_buffer^2)),3) 
     probability_array[i,j,k] 

} 

但我失去了一些重要的東西。任何幫助將不勝感激。 乾杯

回答

3

當並行計算,因爲開銷它引入, 優選並行運行大塊的計算, 而非小 - 外環,而不是內部循環。

但是,在這種情況下,不需要並行化計算:您可以只對它們進行矢量化。

# 3-dimensional analogue of row() and col() 
dim3 <- function(a, i) { 
    stopifnot(length(dim(a)) == 3) 
    r <- a 
    if(i == 1) { r[] <- rep(1:dim(a)[1], dim(a)[2] * dim(a)[3]) } 
    if(i == 2) { r[] <- rep(1:dim(a)[2], each = dim(a)[1], times = dim(a)[3]) } 
    if(i == 3) { r[] <- rep(1:dim(a)[3], each = dim(a)[1] * dim(a)[2]) } 
    r 
} 

probability_array <- array(0,dim=c(length(xvox),length(yvox),length(zvox))) 
i <- dim3(probability_array,1) 
j <- dim3(probability_array,2) 
k <- dim3(probability_array,3) 
for (n in 1:numsphere){ 
    dist_sd <- sqrt(
    (xvox[i]-x[n])^2 + (yvox[j]-y[n])^2 + (zvox[k]-z[n])^2 
) 
    probability_array <- probability_array + 
    # Rounding intermediate results looks suspicious 
    round(exp(-1*(dist_sd-radius[n])^2/(2*radius_buffer^2)),3) 
} 
+0

冠軍!我不知道這是可能的。四捨五入,我錯過了。非常感謝你。 – Mark