Быстрый код в 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_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
}

Когда я вызываю функцию sim_exp со значением, например, 6 значений, это занимает около 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