Более эффективный способ считать пересечения?
У меня есть список 300 000 списков (оптоволоконные дорожки), где каждая дорожка является списком (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 к каждому (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 или больше масок. Существует ли более эффективный способ сделать это, чем итерация для каждой дорожки?
6 ответов
Linearize the voxel coordinates, and put them into two scipy.sparse.sparse.csc matrices.
Let v be the number of voxels, m the number of masks, and t the number of tracks.
Let M be the mask csc matrix, size (m x v), where a 1 at (i,j) means mask i overlaps voxel j.
Let T be the track csc matrix, size (t x v), where a 1 at (k,j) means track k overlaps 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)
I might be wrong on the last one, and I'm not sure css_matrices can be operated on by nonzero & take. You might need to pull out each column in a loop and convert it to a full matrix.
I ran some experiments trying to simulate what I thought was a reasonable amount of data. The code below takes about 2 minutes on a 2-year old MacBook. If you use csr_matrices, it takes about 4 minutes. There is probably a tradeoff depending on how long each track is.
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
Хорошо, я думаю, что наконец-то у меня есть кое-что, что снизит сложность. Этот код должен действительно летать по сравнению с тем, что у вас есть.
Похоже, сначала вам нужно узнать, какие треки совпадают с какими масками, матрица инцидентности .
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 * mT .
Код, который у вас есть, вычисляет только верхний треугольник. Вы можете использовать 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
Для меня это занимает менее секунды для наборов данных, содержащих сотни масок и треков.
Вероятно, вы можете начать с объединения двух функций для создания обоих результатов сразу. Также нет необходимости составлять список комбинаций перед петлей, так как это уже генератор, и это может сэкономить вам немного времени.
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 для питона.
.Если вы сохранили каждое масочное множество точек:
(1,2,3), (1,2,4), (1,3,1) в виде такого словаря: {1: [{2: set([3, 4])}, {3: set([1])}}, то, возможно, вы сможете быстрее проверить совпадения... но, возможно, нет.
Маленькая оптимизация (тот же самый большой-O, слабенький множитель) может быть выполнена путем удаления избыточных операций:
- не вызывайте
setстолько раз на каждой дорожке и маске: вызывайте его один раз для каждого трека и один раз для маски, чтобы установить вспомогательные "параллельные" списки сетов, затем работайте над этими , если любой (сометет):семантически такой же, как, если сометет:, но немного медленнее
не будет иметь драматической разницы, но может помочь с минутной помощью.
Отстойно говорить об очередном инкрементальном улучшении, которое можно было бы сделать, я знаю, но:
Наборы маленьких целых чисел могут быть смоделированы как битовые векторы, используя длинные инты Python. Предположим, вы замените каждый кортеж кортежа на маленький целочисленный id, а затем конвертируете каждую дорожку и каждый набор масок-кордов в набор этих маленьких id. Вы можете представлять эти наборы как длинные инты, делая операцию пересечения немного быстрее (но не асимптотически быстрее)
.