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

Question

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

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

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

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

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

1. количество треков, которые пересекаются между собой маска-маска (создать соединение матрица ectivity)
2. подмножество дорожек, которые пересекают каждую маску, чтобы добавить 1 к каждому (х, у, г) координаты для каждой дорожки в подмножестве (чтобы создать «плотность» изображение)

Я в настоящее время делает часть 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 перекрывает voxel j.
Пусть T - матрица дорожки csc, размер (t x v), где 1 at (k, j) означает, что дорожка k перекрывает voxel 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 может эксплуатироваться при ненулевом & дубля. Вам может потребоваться вытащить каждый столбец в цикле и преобразовать его в полную матрицу.

---

Я провел несколько экспериментов, пытаясь имитировать то, что я считал разумным количеством данных. Приведенный ниже код занимает около 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

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

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 как С модулем для питона.

Answer 3

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

Answer 4

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

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

не будет иметь драматическое значение, но может помочь.

Answer 5

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

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

Answer 6

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

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

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) 

adjacency matrix вы ищете это 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 перед входом в цикл.