Эффективно принимать скользящее среднее разреженных данных и фильтровать выше порогового значения в python

Question

Я получаю мои ноги влажными с некоторым анализом генома и немного застрял. У меня есть очень редкие данные и вам нужно найти места, где скользящее среднее превышает некоторый порог, обозначая каждую точку как 1 или 0. Данные имеют уникальный тип, поэтому я не могу использовать доступные программы для анализа.

Каждая точка представляет собой одну точку (basepair) на человеческом геноме. Для каждого набора данных есть 200 000 000 потенциальных точек. Данные по существу представляют собой список из ~ 12000 пар индексов/значений, где все остальные точки считаются равными нулю. Мне нужно сделать скользящее среднее по всему набору данных и вернуть регионы, где среднее значение превышает пороговое значение.

В настоящее время я читаю каждую точку из набора данных последовательно и строю массив вокруг каждой точки, которую я нахожу, но это очень медленно для больших размеров окна. Есть ли более эффективный способ сделать это, может быть, с scipy или пандами?

Редактировать: магический код Джейми ниже работает отлично (но я не могу еще возвысить)! Я очень благодарен.

Answer 1

Вы можете в векторном выражении всего, с numpy. Я построил этот случайный набор данных (приблиз.) 12,000 индексов от 0 до 199,999,999, и столь же длинного списка случайных поплавков между 0 и 1:

indices = np.unique(np.random.randint(2e8,size=(12000,))) 
values = np.random.rand(len(indices)) 

Тогда я построить массив индексов общего размера окна 2*win+1 вокруг каждого из indices и соответствующего массива, сколько вклад в скользящее среднее по этой точке:

win = 10 

avg_idx = np.arange(-win, win+1) + indices[:, None] 
avg_val = np.tile(values[:, None]/(2*win+1), (1, 2*win+1)) 

Все, что осталось, это выяснить, повторяющиеся индексы и добавляя вклады в скользящей средней вместе:

unique_idx, _ = np.unique(avg_idx, return_inverse=True) 
mov_avg = np.bincount(_, weights=avg_val.ravel()) 

Теперь вы можете получить список индексов, на которых, например. скользящее среднее значение превышает 0,5, как:

unique_idx[mov_avg > 0.5] 

Что касается производительности, сначала превратить выше код в функции:

def sparse_mov_avg(idx, val, win): 
    avg_idx = np.arange(-win, win+1) + idx[:, None] 
    avg_val = np.tile(val[:, None]/(2*win+1), (1, 2*win+1)) 
    unique_idx, _ = np.unique(avg_idx, return_inverse=True) 
    mov_avg = np.bincount(_, weights=avg_val.ravel()) 
    return unique_idx, mov_avg 

и вот некоторые тайминги для нескольких размеров окна, для данных испытаний, описанных в начале:

In [2]: %timeit sparse_mov_avg(indices, values, 10) 
10 loops, best of 3: 33.7 ms per loop 

In [3]: %timeit sparse_mov_avg(indices, values, 100) 
1 loops, best of 3: 378 ms per loop 

In [4]: %timeit sparse_mov_avg(indices, values, 1000) 
1 loops, best of 3: 4.33 s per loop