¿Forma más eficiente de contar intersecciones?

Question

I tienen una lista de 300000 listas (pistas de fibra), donde cada pista es una lista de tuplas (x, y, z)/coordenadas:

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

también tengo un grupo de máscaras, donde cada máscara se define como una lista de (, y, z x) tuplas/coordenadas:

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

estoy tratando de encontrar, para todos los posibles pares de máscaras:

1. el número de pistas que se cruzan entre par máscara-máscara (para crear una conexión matriz ectividad)
2. el subconjunto de pistas que se cruzan cada máscara, con el fin de añadir 1 a cada (x, y, z) de coordenadas para cada pista en el subconjunto (para crear una imagen "densidad")

Actualmente estoy haciendo la parte 1 de esta manera:

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 

y la parte 2 de esta manera:

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 

Usando conjuntos para encontrar las intersecciones se ha acelerado este proceso de manera significativa, pero ambas partes Stil Tomo más de una hora cuando tengo una lista de 70 o más máscaras. ¿Hay una manera más eficiente de hacer esto que iterando para cada pista?

Answer 1

Linealice las coordenadas de vóxel y colóquelas en dos matrices scipy.sparse.sparse.csc.

Sea v el número de vóxeles, m el número de máscaras yt el número de pistas.
Sea M la matriz csc de máscara, tamaño (m x v), donde a 1 en (i, j) significa máscara i se superpone al voxel j.
Sea T la matriz de seguimiento de csc, tamaño (t x v), donde a 1 en (k, j) significa que la pista k se solapa con 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) 

puedo estar equivocado en el último, y no estoy seguro de css_matrices pueden ser operados por distinto de cero & toma. Es posible que deba extraer cada columna en un bucle y convertirla en una matriz completa.

---

Realicé algunos experimentos tratando de simular lo que pensé que era una cantidad razonable de datos. El código siguiente lleva unos 2 minutos en una MacBook de 2 años. Si usa csr_matrices, tarda unos 4 minutos. Probablemente haya una compensación dependiendo de cuánto tiempo es cada pista.

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

Probablemente pueda comenzar combinando las dos funciones para crear ambos resultados a la vez. Además, no es necesario hacer una lista de las combinaciones antes del bucle, ya que es un generador y eso podría ahorrarle algo de tiempo.

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 

Además, este probablemente nunca será "rápida" como en "terminado antes de morir", por lo que la mejor manera es compilar finalmente con Cython como un módulo C para Python.

Answer 3

Si almacenó cada conjunto de máscaras de puntos: (1,2,3), (1,2,4), (1,3,1) como un diccionario como este: {1: [{2: set([3, 4])}, {3: set([1])}]}, podría terminar poder verificar partidos más rápido ... pero tal vez no.

Answer 4

Una optimización de menor importancia (la misma de orden O, sligthly menor multiplicador) se puede conseguir mediante la eliminación de operaciones redundantes:

1. no llaman set tantas veces en cada pista y la máscara: lo llaman una vez por pista y una vez por la máscara, para establecer listas auxiliares "paralelos" de conjuntos, luego trabajar en los
2. if any(someset): es semánticamente lo mismo que if someset: pero un poco más lenta

no hará una diferencia dramática, pero podría ayudarme minuciosamente

Answer 5

Lame para sugerir otra mejora incremental que pudiera haberse producido, lo sé, pero:

conjuntos de números enteros pequeños se pueden modelar como vectores de bits utilizando enteros largos de Python. Supongamos que reemplaza cada tupla con un ID entero pequeño, luego convierte cada pista y cada conjunto de coords de máscara en un conjunto de esos pequeños identificadores. Podrías representar esos conjuntos como largos, haciendo que la operación de intersección sea un poco más rápida (pero no asintóticamente más rápida).

Answer 6

OK, creo que finalmente tengo algo que reducirá la complejidad. Este código realmente debería volar en comparación con lo que tienes.

Parece que primero necesita saber qué pistas coinciden con qué máscaras, 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) 

El adjacency matrix que estás buscando es m * m.T.

El código que tiene hasta ahora solo computa el triángulo superior. Puede usar triu para tomar solo esa mitad.

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

El cálculo de vóxel también puede usar la matriz de incidencia.

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 

Para mí esto funciona en menos de un segundo para los conjuntos de datos que tienen cientos de máscaras y pistas.

Si tiene muchos puntos que aparecen en máscaras pero no están en ninguna pista, puede valer la pena eliminar las entradas para esos puntos de masks_by_point antes de ingresar al bucle.