2015-08-20 48 views
1

我試圖從120,000行csv數據中獲取「surface3d plot」,但是我的代碼太慢,處理時間將超過12小時。 哪個地方應該改變? (口譯()處理時間也將超過12個小時。如果有可能,我想合併mk_surface_data()和插補()。這可能嗎?)爲R中的surface3d準備數據更快的代碼請

library(rgl) 
library(data.table) 
library(akima) 

fv <- cmpfun(function(vec) { 
    return(vec[is.finite(vec)]) 
}) 

mk_surface_data <- cmpfun(function(mat, mean_range = 2, x_div = 100, y_div = 100,defalut_z = 0){ 

    x <- mat[,"x"] 
    y <- mat[,"y"] 

    min_x <- min(fv(x)) 
    max_x <- max(fv(x)) 
    min_y <- min(fv(y)) 
    max_y <- max(fv(y)) 

    sa_x <- max_x - min_x 
    sa_y <- max_y - min_y 

    step_x <- sa_x/x_div 
    step_y <- sa_y/y_div 

    surface_m <- matrix(nrow=0,ncol=3) 

    for(x in 0:x_div){ 

    base_x_range <- min_x + (step_x * x) 
    min_x_range <- base_x_range - (mean_range * step_x) 
    max_x_range <- base_x_range + (mean_range * step_x) 

    for(y in 0:y_div){ 

     base_y_range <- min_y + (step_y * y) 
     min_y_range <- base_y_range - (mean_range * step_y) 
     max_y_range <- base_y_range + (mean_range * step_y) 

     all_z <- mat[((min_x_range < mat[,"x"]) & (max_x_range > mat[,"x"]) & (min_y_range < mat[,"y"]) & (max_y_range > mat[,"y"])),c("z")] 

     if(length(fv(all_z)) > 0){ 
     insert <- c(base_x_range,base_y_range,weighted.mean(all_z,na.rm=T))   
     }else{ 
     insert <- c(base_x_range,base_y_range,defalut_z) 
     } 

     surface_m <- rbind(surface_m,insert) 
    } 
    } 
    colnames(surface_m) <- c("x","y","z") 
    return(as.matrix(surface_m)) 
}) 

# main process 

mean_range = 2 
x_div = 1000 
y_div = 1000 
defalut_z = 0 

mat <- as.matrix(fread("target_file.csv")) 
sdf <- mk_surface_data(mat, mean_range, x_div, y_div,defalut_z) 
interpolated <- interp(sdf[,"x"], sdf[,"y"], sdf[,"z"]) 
plot3d(sdf[,"x"], sdf[,"y"], sdf[,"z"]) 
surface3d(interpolated$x, interpolated$y, interpolated$z,col="green") 

回答

3

的代碼,總是提出了一個線在效率方面紅旗是:

surface_m <- rbind(surface_m,insert) 

基本上你是在你的內部循環內的時間,這是相當低效的成長矩陣surface_m一行(見the R Inferno有關詳細信息,第二圈)。您可以更有效地構建surface_m,例如:

surface_m <- t(apply(expand.grid(y=0:y_div, x=0:x_div), 1, function(yx) { 
    y <- yx[1] 
    x <- yx[2] 
    base_x_range <- min_x + (step_x * x) 
    min_x_range <- base_x_range - (mean_range * step_x) 
    max_x_range <- base_x_range + (mean_range * step_x) 
    base_y_range <- min_y + (step_y * y) 
    min_y_range <- base_y_range - (mean_range * step_y) 
    max_y_range <- base_y_range + (mean_range * step_y) 

    all_z <- mat[((min_x_range < mat[,"x"]) & (max_x_range > mat[,"x"]) & (min_y_range < mat[,"y"]) & (max_y_range > mat[,"y"])),c("z")] 

    if (length(fv(all_z)) > 0){ 
    insert <- c(base_x_range,base_y_range,weighted.mean(all_z,na.rm=T))   
    } else { 
    insert <- c(base_x_range,base_y_range,defalut_z) 
    } 
    return(insert) 
})) 
+0

感謝您的回覆。 我試過你的代碼,處理時間爲4小時。 非常感謝。 – tre56h45