R tối ưu hóa: Làm cách nào để tránh vòng lặp for trong trường hợp này?

Question

Tôi đang cố gắng thực hiện một giao điểm theo dõi di truyền đơn giản trong R và chạy vào các vấn đề hiệu suất lớn, có thể liên quan đến việc tôi sử dụng các vòng lặp.

Trong trường hợp này, tôi có các cửa sổ được xác định trước trong khoảng thời gian 100bp và tôi đang cố gắng tính toán số lượng cửa sổ được bao phủ bởi các chú thích trong danh sách của tôi. Đồ họa, nó trông giống như sau:

  0 100 200 300 400 500 600 
windows: |-----|-----|-----|-----|-----|-----| 

mylist: |-| |-----------| 

Vì vậy, tôi đã viết một số mã để làm việc đó, nhưng nó khá chậm và đã trở thành một nút cổ chai trong mã của tôi:

##window for each 100-bp segment  
windows <- numeric(6) 

##second track 
mylist = vector("list") 
mylist[[1]] = c(1,20) 
mylist[[2]] = c(120,320) 


##do the intersection 
for(i in 1:length(mylist)){ 
    st <- floor(mylist[[i]][1]/100)+1 
    sp <- floor(mylist[[i]][2]/100)+1 
    for(j in st:sp){  
    b <- max((j-1)*100, mylist[[i]][1]) 
    e <- min(j*100, mylist[[i]][2]) 
    windows[j] <- windows[j] + e - b + 1 
    } 
} 

print(windows) 
[1] 20 81 101 21 0 0 

Đương nhiên, điều này đang được được sử dụng trên các tập dữ liệu lớn hơn nhiều so với ví dụ tôi cung cấp ở đây. Thông qua một số hồ sơ, tôi có thể thấy rằng các nút cổ chai là trong vòng, nhưng nỗ lực vụng về của tôi để vectorize nó bằng cách sử dụng * chức năng áp dụng dẫn đến mã chạy một thứ tự cường độ chậm hơn.

Tôi cho rằng tôi có thể viết điều gì đó trong C, nhưng tôi muốn tránh điều đó nếu có thể. Bất cứ ai có thể đề xuất một cách tiếp cận khác mà sẽ tăng tốc độ tính toán này?

Answer 1

Điều "Đúng" cần làm là sử dụng gói bioconductor IRanges, sử dụng cấu trúc dữ liệu IntervalTree để biểu diễn các phạm vi này.

Có cả hai đối tượng trong các đối tượng IRanges riêng của mình, sau đó bạn sẽ sử dụng hàm findOverlaps để giành chiến thắng.

Get nó ở đây:

http://www.bioconductor.org/packages/release/bioc/html/IRanges.html

Vào qua, các bên trong của gói được viết bằng C, vì vậy nó siêu nhanh.

EDIT

Ngày nghĩ thứ hai, nó không phải là nhiều của một slam-dunk như tôi thấy (một liner), nhưng bạn chắc chắn nên bắt đầu sử dụng thư viện này nếu bạn đang làm việc ở tất cả với gen khoảng thời gian (hoặc các loại khác) ... bạn có thể sẽ cần phải thực hiện một số hoạt động và công cụ được thiết lập. Rất tiếc, không có thời gian để cung cấp câu trả lời chính xác.

Tôi chỉ nghĩ rằng điều quan trọng là chỉ thư viện này cho bạn.

Answer 2

Tôi nghĩ rằng tôi đã làm cho nó phức tạp hơn nhiều ... System.time không giúp tôi đánh giá hiệu suất trong một tập dữ liệu nhỏ như vậy.

windows <- numeric(6) 

mylist = vector("list") 
mylist[[1]] = c(1,20) 
mylist[[2]] = c(120,320) 


library(plyr) 

l_ply(mylist, function(x) { 
sapply((floor(x[1]/100)+1) : (floor(x[2]/100)+1), function(z){ 
    eval.parent(parse(text=paste("windows[",z,"] <- ", 
     min(z*100, x[2]) - max((z-1)*100, x[1]) + 1,sep="")),sys.nframe()) 
    })   
}) 

print(windows) 

EDIT

Một sửa đổi để loại bỏ eval

g <- llply(mylist, function(x) { 
ldply((floor(x[1]/100)+1) : (floor(x[2]/100)+1), function(z){ 
     t(matrix(c(z,min(z*100, x[2]) - max((z-1)*100, x[1]) + 1),nrow=2)) 
    })   
}) 

for(i in 1:length(g)){ 
    windows[unlist(g[[i]][1])] <- unlist(g[[i]][2]) 
} 

Answer 3

tôi không có một ý tưởng sáng, nhưng bạn có thể thoát khỏi vòng lặp bên trong, và tăng tốc độ điều một bit. Lưu ý rằng nếu một cửa sổ rơi hoàn toàn trong khoảng thời gian của danh sách của tôi, thì bạn chỉ cần thêm 100 vào phần tử windows tương ứng. Vì vậy, chỉ có st và sp các cửa sổ cần xử lý đặc biệt.

windows <- numeric(100) 
    for(i in 1:length(mylist)){ 
    win <- mylist[[i]]   # for cleaner code 
    st <- floor(win[1]/100)+1 
    sp <- floor(win[2]/100)+1 
    # start and stop are within the same window 
    if (sp == st){ 
     windows[st] <- windows[st] + (win[2]%%100) - (win[1]%%100) +1 
    } 
    # start and stop are in separate windows - take care of edges 
    if (sp > st){ 
     windows[st] <- windows[st] + 100 - (win[1]%%100) + 1 
     windows[sp] <- windows[sp] + (win[2]%%100) 
    } 
    # windows completely inside win 
    if (sp > st+1){ 
     windows[(st+1):(sp-1)] <- windows[(st+1):(sp-1)] + 100 
    }  
    } 

tôi tạo ra một danh sách lớn hơn:

cuts <- sort(sample(1:10000, 70)) # random interval endpoints 
    mylist <- split(cuts, gl(35,2)) 

và có 1,08 giây cho 1000 lần lặp lại của phiên bản này so với 1,72 giây cho 1000 tái tạo cho bản gốc. Với dữ liệu thực, tốc độ sẽ phụ thuộc vào khoảng thời gian trong mylist có xu hướng dài hơn 100 hay không.

Nhân tiện, người ta có thể viết lại vòng lặp bên trong như một hàm riêng biệt, và sau đó lapply nó trên mylist, nhưng điều đó không làm cho nó hoạt động nhanh hơn.

Answer 4

Vì vậy, tôi không hoàn toàn chắc chắn tại sao cửa sổ thứ ba và thứ tư không phải là 100 và 20 vì điều đó sẽ có ý nghĩa hơn đối với tôi. Dưới đây là một lót cho hành vi đó:

Reduce('+', lapply(mylist, function(x) hist(x[1]:x[2], breaks = (0:6) * 100, plot = F)$counts)) 

Lưu ý rằng bạn cần phải xác định giới hạn trên trong breaks, nhưng nó không phải là khó để làm đường chuyền khác để có được nó nếu bạn không biết nó trước .

Answer 5

Được rồi, vì vậy tôi lãng phí quá nhiều thời gian cho việc này và vẫn chỉ có hệ số tăng tốc 3 lần. Bất cứ ai có thể đánh bại điều này?

Mã:

my <- do.call(rbind,mylist) 
myFloor <- floor(my/100) 
myRem <- my%%100 
#Add intervals, over counting interval endpoints 
counts <- table(do.call(c,apply(myFloor,1,function(r) r[1]:r[2]))) 
windows[as.numeric(names(counts))+1] <- counts*101 

#subtract off lower and upper endpoints 
lowerUncovered <- tapply(myRem[,1],myFloor[,1],sum) 
windows[as.numeric(names(lowerUncovered))+1] <- windows[as.numeric(names(lowerUncovered))+1] - lowerUncovered 
upperUncovered <- tapply(myRem[,2],myFloor[,2],function(x) 100*length(x) - sum(x)) 
windows[as.numeric(names(upperUncovered))+1] <- windows[as.numeric(names(upperUncovered))+1] - upperUncovered 

Các thử nghiệm:

mylist = vector("list") 
for(i in 1:20000){ 
    d <- round(runif(1,,500)) 
    mylist[[i]] <- c(d,d+round(runif(1,,700))) 
} 

windows <- numeric(200) 


new_code <-function(){ 
    my <- do.call(rbind,mylist) 
    myFloor <- floor(my/100) 
    myRem <- my%%100 
    counts <- table(do.call(c,apply(myFloor,1,function(r) r[1]:r[2]))) 
    windows[as.numeric(names(counts))+1] <- counts*101 

    lowerUncovered <- tapply(myRem[,1],myFloor[,1],sum) 
    windows[as.numeric(names(lowerUncovered))+1] <- windows[as.numeric(names(lowerUncovered))+1] - lowerUncovered 

    upperUncovered <- tapply(myRem[,2],myFloor[,2],function(x) 100*length(x) - sum(x)) 
    windows[as.numeric(names(upperUncovered))+1] <- windows[as.numeric(names(upperUncovered))+1] - upperUncovered 

    #print(windows) 
} 


#old code 
old_code <- function(){ 
    for(i in 1:length(mylist)){ 
     st <- floor(mylist[[i]][1]/100)+1 
     sp <- floor(mylist[[i]][2]/100)+1 
     for(j in st:sp){  
      b <- max((j-1)*100, mylist[[i]][1]) 
      e <- min(j*100, mylist[[i]][2]) 
      windows[j] <- windows[j] + e - b + 1 
     } 
    } 
    #print(windows) 
} 

system.time(old_code()) 
system.time(new_code()) 

Kết quả:

> system.time(old_code()) 
    user system elapsed 
    2.403 0.021 2.183 
> system.time(new_code()) 
    user system elapsed 
    0.739 0.033 0.588 

Rất bực bội rằng thời gian hệ thống về cơ bản là 0, nhưng thời gian quan sát rất tuyệt quá. Tôi đặt cược nếu bạn đã đi xuống C bạn sẽ nhận được một tốc độ 50-100X.