Mã nhanh hơn trong R

Question

FYI: Tôi đã chỉnh sửa điều này đáng kể kể từ ấn bản đầu tiên của mình. Mô phỏng này đã giảm từ 14 đến 14 phút.

Tôi mới lập trình nhưng tôi đã thực hiện mô phỏng cố gắng sao chép vô tính trong sinh vật và định lượng sự khác biệt về số nhiễm sắc thể giữa bố mẹ và con gái. Mô phỏng chạy rất chậm. Mất khoảng 6 giờ để hoàn thành. Tôi muốn biết những gì sẽ là cách tốt nhất để làm cho mô phỏng chạy nhanh hơn.

Các sinh vật kỹ thuật số này có x số nhiễm sắc thể. Không giống như hầu hết các sinh vật, các nhiễm sắc thể đều độc lập với nhau, vì vậy chúng có cơ hội bình đẳng được chuyển vào sinh vật con gái.

Trong trường hợp này, sự phân bố nhiễm sắc thể vào một tế bào con gái theo sau sự phân bố nhị thức với xác suất 0,5.

Hàm sim_repo lấy ma trận của các sinh vật số với số lượng nhiễm sắc thể đã biết và đặt chúng qua 12 thế hệ sao chép. Nó sao chép các nhiễm sắc thể này và sau đó sử dụng hàm rbinom để tạo ngẫu nhiên một số. Số này sau đó được gán cho một tế bào con gái. Vì không có nhiễm sắc thể nào bị mất trong quá trình sinh sản vô tính, tế bào con khác nhận được nhiễm sắc thể còn lại. Điều này sau đó được lặp lại cho G số thế hệ. Sau đó, một giá trị được lấy mẫu từ mỗi hàng trong ma trận.

sim_repo = function(x1, G=12, k=1, t=25, h=1000) { 

      # x1 is the list of copy numbers for a somatic chromosome 
      # G is the number of generations, default is 12 
      # k is the transfer size, default is 1 
      # t is the number of transfers, default is 25 
      # h is the number of times to replicate, default is 1000 

      dup <- x1 * 2 # duplicate the initial somatic chromosome copy number for replication 
      pop <- 1 # set generation time 
      set.seed(11) 
      z <- matrix(rbinom(n=rep(1,length(dup)),size = as.vector(dup),prob = 0.5),nrow = nrow(dup)) # amount of somatic chromosome is distributed to one of the daughter cells 
      z1 <- dup - z # as no somatic chromosomes are lost, the other daughter cells receives the remainder somatic chromosomes 
      x1 <- cbind(z, z1) # put both in a matrix 

      for (pop in 1:G) { # this loop does the replication for each cell in each generation 
       pop <- 1 + pop # number of generations. This is a count for the for loop 
       dup <- x1 * 2 # double the somatic chromosomes for replication 
       set.seed(11) 
       z <- matrix(rbinom(n=rep(1,length(dup)),size = as.vector(dup),prob = 0.5),nrow = nrow(dup)) # amount of somatic c hromosomes distributed to one of the daughter cells 
       z1 <- dup - z # as no somatic chromosomes are lost, the other daughter cells receives the remainder somatic chromosomes 
       x1 <- cbind(z, z1) # put both in a matrix 
       } 

      # the following for loop randomly selects one cell in the population that was created 
      # the output is a matrix of 1 column 
      x1 <- matrix(apply(x1, 1, sample, size=k), ncol=1) 
      x1 
    } 

Tôi nghiên cứu của mình Tôi quan tâm đến sự thay đổi phương sai trong nhiễm sắc thể của tổ tiên ban đầu và thời điểm cuối cùng trong mô phỏng này. Hàm sau đại diện cho việc chuyển một ô vào một môi trường sống mới. Nó lấy đầu ra từ hàm sim_re p và sử dụng nó để tạo ra nhiều thế hệ. Sau đó nó tìm ra phương sai giữa các hàng trong cột ma trận đầu tiên và cuối cùng và tìm thấy sự khác biệt giữa chúng.

# The following function is mostly the same as I talked about in the description. 
    # The only difference is I changed some aspects to take into account I am using 
    # matrices and not lists. 
    # The function outputs the difference between the intial variance component between 
    # 'cell lines' with the final variance after t number of transfers 

sim_exp = function(x1, G=12, k=1, t=25, h=1000) { 

    xn <- matrix(NA, nrow(x1), t) 
    x <- x1 
    xn[,1] <- x1 
    for (l in 2:t) { 
     x <- sim_repo(x, G, k, t, h) 
     xn[, l] <- x 
    } 

    colvar <- matrix(apply(xn,2,var),ncol=ncol(xn)) 
    ivar <- colvar[,1] 
    fvar <- colvar[,ncol(xn)] 
    deltavar <- fvar - ivar 
    deltavar 
} 

Tôi cần nhân rộng mô phỏng này h số lần. Do đó, tôi đã thực hiện chức năng sau sẽ gọi hàm sim_exph số lần.

sim_1000 = function(x1, G=12, k=1, t=25, h=1000) { 
    xn <- vector(length=h) 
    for (l in 2:h) { 
     x <- sim_exp(x1, G, k, t, h) 
     xn[l] <- x 
    } 
     xn 
} 

Khi tôi gọi hàm sim_exp có giá trị như 6 giá trị cần khoảng 52 giây để hoàn thành.

x1 <- matrix(data=c(100,100,100,100,100,100),ncol=1) 
system.time(sim_1000(x1,h=1)) 
    user system elapsed 
    1.280 0.105 1.369 

Nếu tôi có thể nhanh hơn, tôi có thể hoàn thành nhiều mô phỏng và áp dụng mô hình lựa chọn trên mô phỏng.

đầu vào của tôi sẽ giống như x1 sau, một ma trận với mỗi sinh vật tổ tiên trên hàng riêng của mình:

x1 <- matrix(data=c(100,100,100,100,100,100),ncol=1) # a matrix of 6 organisms 

Khi tôi chạy:

a <- sim_repo(x1, G=12, k=1) 

sản lượng dự kiến ​​của tôi sẽ là:

a 
    [,1] 
[1,] 137 
[2,] 82 
[3,] 89 
[4,] 135 
[5,] 89 
[6,] 109 

system.time(sim_repo(x1)) 
    user system elapsed 
    1.969 0.059 2.010 

Khi tôi gọi hàm sim_exp,

b < - sim_exp (x1, G = 12, k = 1, t = 25)

nó gọi hàm sim_repo lần và kết quả đầu ra G:

b 
[1] 18805.47 

Khi tôi gọi là sim_1000 chức năng, Tôi thường sẽ thiết lập h đến 1000, nhưng ở đây tôi sẽ đặt nó ở 2. Vì vậy, ở đây sim_1000 sẽ gọi sim_exp và nhân bản nó 2 lần.

c <- sim_1000(x1, G=12, k=1, t=25, h=2) 
c 
[1] 18805.47 18805.47 

Answer 1

Như đã đề cập bởi những người khác trong các ý kiến, nếu chúng ta chỉ nhìn vào chức năng sim_repo, và thay thế dòng:

dup <- apply(x1, c(1,2),"*",2) 

với

dup <- x1 * 2 

các dòng

z <- apply(dup,c(1,2),rbinom,n=1,prob=0.5) 

với

z <- matrix(rbinom(n=rep(1,length(dup)),size = as.vector(dup),prob = 0.5),nrow = nrow(dup)) 

và bên trong cho vòng lặp với

x1 <- matrix(apply(x1,1,sample,size = 1), ncol=1) 

tôi nhận được một, tốt, tăng tốc độ lớn:

system.time(sim_exp(x1)) 
    user system elapsed 
    0.655 0.017 0.686 
> system.time(sim_expOld(x1)) 
    user system elapsed 
21.445 0.128 21.530 

Và tôi xác nhận rằng nó đang làm điều tương tự:

set.seed(123) 
out1 <- sim_exp(x1) 

set.seed(123) 
out2 <- sim_expOld(x1) 

all.equal(out1,out2) 
> TRUE 

Và đó thậm chí không phải là việc đào bới vào quá trình phân bổ trước, điều này thực sự khó mà không thiết kế lại hoàn toàn mọi thứ, theo cách bạn đã cấu trúc mã của mình.

Và đó cũng là thậm chí không bắt đầu nhìn vào cho dù bạn thực sự cần cả ba chức năng ...