Более эффективный способ подсчета пересечений?

Question

У меня есть список из 300000 списков (волоконных дорожек), где каждая дорожка представляет собой список (x, y, z) кортежей / координат:

tracks=
[[(1,2,3),(3,2,4),...]
 [(4,2,1),(5,7,3),...]
 ...
]

У меня также есть группа масок, где каждая маска определяется как список (x, y, z) кортежей / координат:

mask_coords_list=
[[(1,2,3),(8,13,4),...]
 [(6,2,2),(5,7,3),...]
 ...
]

Я пытаюсь найти для всех возможных пар масок:

1. количество дорожек, пересекающих каждую пару маска-маска (для создания матрицы связности)
2. подмножество дорожек, которые пересекают каждую маску, чтобы добавить 1 к каждой (x, y, z) координате для каждой дорожки в подмножестве (чтобы создать изображение "плотности")

В настоящее время я делаю часть 1 примерно так:

def mask_connectivity_matrix(tracks,masks,masks_coords_list):
    connect_mat=zeros((len(masks),len(masks)))
    for track in tracks:
        cur=[]
        for count,mask_coords in enumerate(masks_coords_list):
            if any(set(track) & set(mask_coords)):
                cur.append(count)
            for x,y in list(itertools.combinations(cur,2)):
                connect_mat[x,y] += 1

и часть 2 примерно так:

def mask_tracks(tracks,masks,masks_coords_list):
    vox_tracks_img=zeros((xdim,ydim,zdim,len(masks)))
    for track in tracks:
        for count,mask in enumerate(masks_coords_list):
            if any(set(track) & set(mask)):
                for x,y,z in track:
                    vox_tracks_img[x,y,z,count] += 1

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

Answer 1

Линеаризуйте координаты вокселей и поместите их в две матрицы scipy.sparse.sparse.csc.

Пусть v будет количеством вокселей, m количеством масок и t количеством дорожек.
Пусть M - маска csc-матрицы, размер (m x v), где a 1 at (i, j) означает, что маска i перекрывает воксель j.
Пусть T будет матрицей csc дорожек, размер (t x v), где a 1 at (k, j) означает, что дорожка k перекрывает воксель j.

Overlap = (M * T.transpose() > 0)  # track T overlaps mask M  
Connected = (Overlap * Overlap.tranpose() > 0) # Connected masks
Density[mask_idx] = numpy.take(T, nonzero(Overlap[mask_idx, :])[0], axis=0).sum(axis=0)

Возможно, я ошибаюсь в последнем, и я не уверен, что css_matrices может быть обработан ненулевым & take. Возможно, вам придется извлечь каждый столбец в цикле и преобразовать его в полную матрицу.

---

Я провел несколько экспериментов, пытаясь смоделировать то, что я считал разумным количеством данных. Код ниже занимает около 2 минут на 2-летнем MacBook. Если вы используете csr_matrices, это займет около 4 минут. Вероятно, существует компромисс в зависимости от длины каждого трека.

from numpy import *
from scipy.sparse import csc_matrix

nvox = 1000000
ntracks = 300000
nmask = 100

# create about 100 entries per track
tcoords = random.uniform(0, ntracks, ntracks * 100).astype(int)
vcoords = random.uniform(0, nvox, ntracks * 100).astype(int)
d = ones(ntracks * 100)
T = csc_matrix((d,  vstack((tcoords, vcoords))), shape=(ntracks, nvox), dtype=bool)

# create around 10000 entries per mask
mcoords = random.uniform(0, nmask, nmask * 10000).astype(int)
vcoords = random.uniform(0, nvox, nmask * 10000).astype(int)
d = ones(nmask * 10000)
M = csc_matrix((d, vstack((mcoords, vcoords))), shape=(nmask, nvox), dtype=bool)

Overlap = (M * T.transpose()).astype(bool) # mask M overlaps track T
Connected = (Overlap * Overlap.transpose()).astype(bool) # mask M1 and M2 are connected
Density = Overlap * T.astype(float) # number of tracks overlapping mask M summed across voxels

Answer 2

ОК, я думаю, что наконец-то есть что-то, что уменьшит сложность. Этот код должен действительно летать по сравнению с тем, что у вас есть.

Кажется, сначала нужно узнать, какие дорожки совпадают с какими масками, матрица инцидентности .

import numpy
from collections import defaultdict

def by_point(sets):
    d = defaultdict(list)
    for i, s in enumerate(sets):
        for pt in s:
            d[pt].append(i)
    return d

def calc(xdim, ydim, zdim, mask_coords_list, tracks):
    masks_by_point = by_point(mask_coords_list)
    tracks_by_point = by_point(tracks)

    a = numpy.zeros((len(mask_coords_list), len(tracks)), dtype=int)
    for pt, maskids in masks_by_point.iteritems():
        for trackid in tracks_by_point.get(pt, ()):
            a[maskids, trackid] = 1
    m = numpy.matrix(a)

Матрица смежности , которую вы ищете, m * m.T.

Код, который у вас есть, вычисляет только верхний треугольник. Вы можете использовать triu, чтобы захватить только половину.

    am = m * m.T  # calculate adjacency matrix
    am = numpy.triu(am, 1)  # keep only upper triangle
    am = am.A  # convert matrix back to array

В расчете вокселей также может использоваться матрица инцидентности.

    vox_tracks_img = numpy.zeros((xdim, ydim, zdim, len(mask_coords_list)), dtype=int)
    for trackid, track in enumerate(tracks):
        for x, y, z in track:
            vox_tracks_img[x, y, z, :] += a[:,trackid]
    return am, vox_tracks_img

Для меня это занимает меньше секунды для наборов данных с сотнями масок и дорожек.

Если у вас много точек, которые отображаются в масках, но отсутствуют на дорожках, может быть целесообразно удалить записи для этих точек из masks_by_point перед входом в цикл.

Answer 3

Хм, чтобы предложить еще одно постепенное улучшение, которое может быть сделано, я знаю, но:

Наборы маленьких целых чисел могут быть смоделированы как битовые векторы с использованием длинных целых чисел Python. Предположим, вы заменяете каждый кортеж небольшим целочисленным идентификатором, а затем конвертируете каждый трек и каждый набор координат-масок в набор этих маленьких идентификаторов. Вы могли бы представить эти множества как длинные целые, делая операцию пересечения немного быстрее (но не асимптотически быстрее).

Answer 4

Небольшая оптимизация (такой же большой-O, чуть меньший множитель) может быть получена путем удаления избыточных операций:

1. не вызывайте set так много раз на каждой дорожке и маске: вызывайте ее один раз для каждой дорожки и один раз для маски, чтобы настроить вспомогательные "параллельные" списки наборов, а затем работайте с этими
2. if any(someset): семантически совпадает с if someset:, но немного медленнее

Не будет существенной разницы, но может помочь.

Answer 5

Если вы сохранили каждый набор точек маски: (1,2,3), (1,2,4), (1,3,1) в виде словаря, подобного этому: {1: [{2: set([3, 4])}, {3: set([1])}]}, вы можете в конечном итоге быстрее проверять совпадения ... но, возможно, нет ,

Answer 6

Вероятно, вы можете начать с объединения двух функций, чтобы создать оба результата одновременно. Также нет необходимости составлять список комбинаций перед циклом, так как это уже генератор, и это может сэкономить вам время.

def mask_connectivity_matrix_and_tracks(tracks,masks,masks_coords_list):
    connect_mat=zeros((len(masks),len(masks)))
    vox_tracks_img=zeros((xdim,ydim,zdim,len(masks)))
    for track in tracks:
        cur=[]
        for count,mask_coords in enumerate(masks_coords_list):
            if any(set(track) & set(mask_coords)):
                cur.append(count)
                for x,y,z in track:
                    vox_tracks_img[x,y,z,count] += 1
            for x,y in itertools.combinations(cur,2):
                connect_mat[x,y] += 1

Кроме того, это, вероятно, никогда не будет «быстрым», как в «законченном до нашей смерти», поэтому лучший способ - это в конечном итоге скомпилировать его с помощью Cython как модуля c для python.