2012-01-01 68 views
4

下面的代碼生成所需的輸出。然而,缺乏矢量化意味着它運行得非常緩慢。我如何加快速度?矢量化/加速嵌套循環代碼

我已將dput結果從部分指示性數據中提取出來。

輸入dput小號

  1. StandRef輸入

    structure(list(id = structure(c(43L, 50L, 17L, 45L, 9L, 5L, 49L, 
    33L, 48L, 39L, 71L, 64L, 44L, 47L, 58L, 24L, 15L, 37L, 14L, 11L, 
    26L, 57L, 4L, 30L, 72L, 21L, 23L, 60L, 38L, 59L, 29L, 19L, 6L, 
    46L, 36L, 3L, 63L, 55L, 51L, 35L, 10L, 7L, 16L, 73L, 42L, 52L, 
    41L, 27L, 25L, 61L, 20L, 70L, 53L, 18L, 31L, 22L, 1L, 8L, 2L, 
    40L, 65L, 67L, 28L, 56L, 13L, 32L, 54L, 66L, 68L, 34L, 12L, 69L, 
    62L), .Label = c("ID 1009445", "ID 120763", "ID 133883", "ID 136398", 
    "ID 171850", "ID 192595", "ID 197597", "ID 216406", "ID 21888", 
    "ID 230940", "ID 23777", "ID 282791", "ID 306348", "ID 309745", 
    "ID 326928", "ID 344897", "ID 34974", "ID 350157", "ID 391831", 
    "ID 402479", "ID 43010", "ID 484078", "ID 484697", "ID 537134", 
    "ID 562259", "ID 562455", "ID 567042", "ID 572866", "ID 578945", 
    "ID 595683", "ID 59759", "ID 598460", "ID 603611", "ID 603757", 
    "ID 607991", "ID 60976", "ID 622720", "ID 646989", "ID 656144", 
    "ID 668807", "ID 669435", "ID 720522", "ID 740555", "ID 745499", 
    "ID 746001", "ID 783969", "ID 78979", "ID 792426", "ID 793541", 
    "ID 797860", "ID 806559", "ID 810517", "ID 826054", "ID 837609", 
    "ID 839287", "ID 867918", "ID 869788", "ID 875380", "ID 876870", 
    "ID 882220", "ID 893116", "ID 895909", "ID 899050", "ID 900143", 
    "ID 908100", "ID 912185", "ID 916371", "ID 916620", "ID 957879", 
    "ID 966195", "ID 993247", "ID 998911", "ID 999610"), class = "factor"), 
        region = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
        1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
        1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
        2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
        2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
        2L), location = c(259090L, 559306L, 2227063L, 2369217L, 4026978L, 
        4211264L, 4679449L, 5105226L, 5106345L, 5344670L, 5473601L, 
        5476528L, 5871970L, 6461228L, 6700029L, 6708265L, 7639959L, 
        9297695L, 10254788L, 10328812L, 11102816L, 11568295L, 11720437L, 
        12843457L, 14012506L, 14156669L, 14632300L, 14641938L, 15298211L, 
        15468425L, 15534406L, 16279682L, 16699353L, 17226952L, 17320785L, 
        269017L, 453097L, 828833L, 954610L, 954842L, 1066378L, 1217332L, 
        1253530L, 1277716L, 1292857L, 1337952L, 1439657L, 1452989L, 
        1712345L, 1758035L, 2601630L, 2640359L, 2778095L, 3151129L, 
        3369931L, 3399080L, 3529525L, 3810217L, 3821120L, 3841588L, 
        3901557L, 4111633L, 4220440L, 4528632L, 4665450L, 5099307L, 
        5260242L, 5958770L, 5966356L, 6137405L, 6246065L, 6297231L, 
        6807949L)), .Names = c("id", "region", "location"), class = "data.frame", row.names = c(NA, 
    -73L)) 
    
  2. 兩個樣本輸入

樣品1

 structure(list(region = c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L), 
     begin = c(0L, 2259252L, 5092077L, 9158205L, 0L, 135094L, 
     941813L, 5901391L, 6061324L), finish = c(2259252L, 5092077L, 
     9158205L, 20463033L, 135094L, 941813L, 5901391L, 6061324L, 
     7092402L), sed = c(3.98106154985726, 7.51649828394875, 5.15440228627995, 
     2.67456624889746, 7.54309412557632, 4.17413910385221, 7.47043058509007, 
     6.13362524658442, 1.00084994221106)), .Names = c("region", 
     "begin", "finish", "sed"), class = "data.frame", row.names = c(NA, 
     -9L)) 

樣品2

 structure(list(region = c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L), 
     begin = c(0L, 2253252L, 7091077L, 9120205L, 0L, 135094L, 
     941813L, 5901391L, 6061324L), finish = c(2253252L, 7091077L, 
     9120205L, 17463033L, 135094L, 941813L, 5901391L, 6061324L, 
     7092402L), sed = c(3.31830840984048, 1.38014704208403, 6.13049140975458, 
     2.10349875097134, 0.48170587509345, 0.13058713509175, 9.13509713513509, 
     6.13047153058701, 3.81734081501503)), .Names = c("region", 
     "begin", "finish", "sed"), class = "data.frame", row.names = c(NA, 
     -9L)) 

Unvectorized代碼

matchLocationsToRegions <- function(path) {  
# get list of data files (around 500 of these; only dput of 2 given: sample262519 and sample252519) 
setwd(path,sep="",collapse=NULL) 
data_files <- list.files() 

# read in template file with complete regional boundaries 
standRef <- read.table(paste(path, "StandRef.txt",sep="",collapse=NULL), header=TRUE, sep="\t") 

# pre-allocate a df with row dimensions of standRef and num of columns according to num of data files 
sediment.df <- as.data.frame(matrix(NA,nrow=nrow(standRef),ncol=length(data_files))) 
colnames(sediment.df) <- data_files 
rownames(sediment.df) <- standRef[,1] 

# create a counter for columns filled 
col_counter <- 1  

for (file in data_files) { 
    # read in current, processed data 
    sample <- read.table(file, header=TRUE, sep="\t")   

    # pre-allocate vectors for sedimentation data vector 
    sed <- rep(NA, nrow(standRef)) 

    # create a variable to track end boundary for a particular sample_ID 
    end_tracker <- 1 

    index <- unlist(lapply (unique(standRef$region), function(reg) { 
      reg.filter <- which(standRef$region == reg) 
      samp.filter <- which(sample$region == reg) 
      samp.filter[cut(standRef$location[reg.filter],c(0L,sample$finish[samp.filter]),labels=F)] 
     })) 
    sed <- sample$sed[index] 

    # fill in next, unfilled column of relevant df with data from relevant vector 
    sediment.df[col_counter] <- sed 

    # update column counter variable 
    col_counter <- col_counter + 1 
}  

# save df as a table 
write.table(sediment.df,file="samples_sed.txt", row.names=TRUE, sep="\t") 
} 

運行Rprof表明 "scan" "read.table" "matchLocationsToRegions""type.convert" "read.table" "matchLocationsToRegions" 爲主運行。據推測,有一個瓶頸,通過套住超過此線:

sample <- read.table(file, header=TRUE, sep="\t")  

更新:for環比地區已被替換,執行快很多(H/T西蒙Urbanek)。但其餘的很慢。

+0

我建議你先嚐試運行'Rprof()',看看瓶頸在哪裏。 – 2012-01-01 12:54:59

+0

@RomanLuštrik添加了一些來自'Rprof'的信息 – Kaleb 2012-01-01 13:35:25

+0

嘗試統計SE,有相當多的人在R有技能。 – check123 2012-01-01 13:37:26

回答

1

您可以刪除容易循環:

sediment.df <- as.data.frame(lapply(data_files, function(file) { 
    sample <- read.table(file, header=TRUE, sep="\t")   
    index <- unlist(lapply (unique(standRef$region), function(reg) { 
     reg.filter <- which(standRef$region == reg) 
     samp.filter <- which(sample$region == reg) 
      samp.filter[cut(standRef$location[reg.filter],c(0L,sample$finish[samp.filter]),labels=F)] 
    })) 
    sample$sed[index] 
})) 
colnames(sediment.df) <- data_files 
rownames(sediment.df) <- standRef[,1] 

然而,這不是不可能,很多的時間在read.table花費,所以你可以使用scan,考慮A)B)只創建一個文件與所有樣本(例如使用額外的列來定義樣本),因此您不需要加載許多文件。

+0

我知道這不是問題,但如果我不得不返回'sample $ borewidth [index]'以及'sample $ sed [index]'並將其放入其自己的數據框(例如'borewidth.df' ),我該怎麼做,而不必另外創建一個單獨的'as.data.frame(lapply ...'結構,這不會很好,因爲它需要兩次運行示例文件?創建一個聚合樣本文件是不可能的,因爲這些文件是不同的。謝謝。 – Kaleb 2012-01-03 00:35:38

+1

多於一個矢量結果會失去一些便利 - 例如,您可以返回一個矩陣,而不是一個矢量,但是您需要一個額外的步驟將所有列拆分爲兩個(還有其他許多方法)。對於彙總文件 - 您的意思是什麼?您可以使用固定的列,因此您可以將其放入一個帶有額外ID列表示源文件的文件中。 – 2012-01-03 01:37:57

+0

我想有可能有一個聚合文件。但是,這並不意味着有''sourceid「」region「」begin「」finish「」sed「」borewidth「'columns x 500(3000 cols)或者我誤解了你?這是一個關於內存的問題嗎? – Kaleb 2012-01-03 09:43:34