r - R中更快的代码

仅供引用:自第一版以来，我对此进行了重大编辑。模拟时间从原来的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_re p输出，并使用它生成更多的世代。然后，它找到第一个和最后一个矩阵列中各行之间的方差，并找到它们之间的差异。

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

最佳答案

正如其他人在评论中所提到的，如果我们仅查看函数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

而且这甚至不涉及预分配，考虑到您已构造代码的方式，如果不完全重新设计内容，这实际上可能会很困难。

而且这甚至还没有开始考虑您是否真的需要全部三个功能...

关于r - R中更快的代码，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/9594493/