Эффективно вычислить гистограмму парных разностей большого вектора в R? - PullRequest
Новая
       22

Эффективно вычислить гистограмму парных разностей большого вектора в R?

10 голосов
/ 02 марта 2012

Я работаю с большим вектором целых чисел в R (около 10 миллионов целых чисел), и мне нужно найти каждую отдельную пару целых чисел из этого вектора, которые отличаются на 500 или меньше, и составить гистограмму их различий (т.е. для каждая пара, вторая минус первая).

Вот совершенно не-векторизованный код для того, чтобы делать то, что я хочу, очень медленно: 

# Generate some random example data
V <- round(rnorm(100) * 1000)

# Prepare the histogram
my.hist <- rep(0, 500)
names(my.hist) <- as.character(seq(1,500))
for (x1 in V) {
    for (x2 in V) {
        difference = x2 - x1
        if (difference > 0 && difference <= 500) {
            my.hist[difference] = my.hist[difference] + 1
        }
    }
}

(Предположим, что каждое целое число уникально, поэтому бит difference > 0 в порядке. Это разрешено, потому что на самом деле меня не волнуют случаи, когда разница равна нулю.)

Вот некоторый код, который векторизует внутренний цикл: 

my.hist2 <- rep(0, 500)
names(my.hist2) <- as.character(seq(1,500))
for (x1 in V) {
    differences <- V[V > x1 & V <= x1+500] - x1
    difftable <- table(differences)
    my.hist2[names(difftable)] = my.hist2[names(difftable)] + difftable
}

Это, безусловно, быстрее, чем первая версия. Однако даже этот вариант уже слишком медленный, когда V содержит только 500000 элементов (полмиллиона).

Я могу сделать это без каких-либо явных циклов следующим образом: 

X <- combn(V, 2)
# X is a matrix with two rows where each column represents a pair
diffs <- abs(X[2,] - X[1,])
my.hist3 <- table(diffs[diffs <= 500])

Однако матрица X будет содержать 10e6 * (10e6 - 1) / 2 или около 50 000 000 000 000 столбцов, что может быть проблемой.

Так есть ли способ сделать это без явного зацикливания (слишком медленного) или расширения матрицы всех пар (слишком большого)?

Если вам интересно, зачем мне это делать, я реализую это: http://biowhat.ucsd.edu/homer/chipseq/qc.html#Sequencing_Fragment_Length_Estimation

1 Ответ

16 голосов
/ 02 марта 2012

Одним из возможных улучшений является сортировка данных: пары (i, j), для которых расстояние меньше 500, будут тогда близки к диагонали, и вам не нужно исследовать все значения.

Код будет выглядеть так (он все еще очень медленный). 

n <- 1e5
k <- 500
V <- round(rnorm(n) * n * 10)
V <- as.integer(V)
V <- sort(V)
h <- rep(0,k)

for(i in 1:(n-1)) {
  for(j in (i+1):n) {
    d <- V[j] - V[i]
    if( d > k ) break
    if( d > 0 ) h[d] <- h[d]+1
  }
}

Редактировать: Если вы хотите что-то намного быстрее, вы можете написать цикл в C. Это занимает 1 с на 10 миллионов элементов.

n <- 10e6
k <- 500
V <- round(rnorm(n) * n * 10)
V <- as.integer(V)
V <- sort(V)
h <- rep(0,k)

library(inline)
sig <- signature(n="integer", v="integer", k="integer", h="integer")
code <- "
  for( int i = 0; i < (*n) - 1; i++ ) {
    for( int j = i + 1; j < *n; j++ ) {
      int d = v[j] - v[i];
      if( d > *k ) break;
      if( d > 0 ) h[d-1]++;
    }
  }
"
f <- cfunction( sig, code, convention=".C" )
h <- f(n,V,k,h)$h
...