r - R 中的代码更快

Question

仅供引用:自第一版以来，我已对此进行了重大编辑。此模拟时间已从 14 小时减少到 14 分钟。

我是编程新手，但我做了一个模拟，试图跟踪生物体中的无性复制，并量化亲代生物体和子代生物体之间染色体数目的差异。模拟运行速度非常慢。大约需要6个小时才能完成。我想知道什么是使模拟运行得更快的最佳方法。

这些数字生物有 x 条染色体。与大多数生物体不同，染色体彼此独立，因此它们有相同的机会转移到任一子生物体中。

在这种情况下，子细胞中的染色体分布遵循二项分布，概率为 0.5。

函数 sim_repo 采用具有已知染色体数量的数字生物体矩阵，并将它们进行 12 代复制。它复制这些染色体，然后使用 rbinom 函数随机生成一个数字。然后将该号码分配给子单元。由于无性繁殖过程中没有染色体丢失，另一个子细胞接收剩余的染色体。然后重复G代。然后从矩阵中的每一行采样一个值。

 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
    }

在我的研究中，我对初始祖先生物体的染色体方差变化以及该模拟中的最终时间点感兴趣。以下函数表示将细胞转移到新的生活环境中。它获取函数 sim_rep 的输出，并使用它来生成更多代。然后，它会找出矩阵第一列和最后一列中的行之间的方差，并找出它们之间的差异。

    # 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
}  

我需要重复此模拟 h 次。因此，我创建了以下函数，该函数将调用函数 sim_exp h 次。

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
}

当我使用 6 个值之类的值调用 sim_exp 函数时，大约需要 52 秒才能完成。

 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 

如果我能更快地完成，那么我就可以完成更多的模拟，并在模拟上应用选择模型。

我的输入将类似于以下 x1，一个矩阵，每个祖先生物体位于其自己的行上:

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

当我运行时:

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

我的预期输出是:

 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 

当我调用 sim_exp 函数时，

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

它调用 sim_repo 函数 G 次并输出:

 b
[1] 18805.47

当我调用sim_1000函数时，我通常会将h设置为1000，但这里我将其设置为2。所以这里sim_1000会调用sim_exp 并复制2次。

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

Answer 1

正如其他人在评论中提到的，如果我们只查看函数 sim_repo，并替换该行:

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

与

dup <- x1 * 2

线条

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

与

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

以及内部 for 循环

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

我得到了很大的速度提升:

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 

我验证了它正在做同样的事情:

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

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

all.equal(out1,out2)
> TRUE

这甚至还没有深入研究预分配，考虑到您构建代码的方式，如果不完全重新设计，这实际上可能很困难。

这还不是开始来看看您是否真的需要所有这三个功能......