Calcul des interactions efficaces entre paires de particules

Question

J'ai une simulation à N corps qui génère une liste de positions de particules, pour plusieurs pas de temps dans la simulation. Pour une trame donnée, je souhaite générer une liste des couples d'indices de particules (i, j) tels que dist(p[i], p[j]) < masking_radius. Essentiellement, je crée une liste de paires "d'interaction", où les paires sont à une certaine distance les unes des autres. Ma mise en œuvre actuelle ressemble à quelque chose comme ceci:

interaction_pairs = [] 

# going through each unique pair (order doesn't matter) 
for i in range(num_particles): 
    for j in range(i + 1, num_particles): 
     if dist(p[i], p[j]) < masking_radius: 
      interaction_pairs.append((i,j)) 

En raison du grand nombre de particules, ce processus prend beaucoup de temps (> 1 heure par test), et il est limite sérieusement à ce que je dois faire avec le Les données. Je me demandais s'il y avait un moyen plus efficace de structurer les données de telle sorte que le calcul de ces paires serait plus efficace au lieu de comparer toutes les combinaisons possibles de particules. Je regardais KDTrees, mais je ne pouvais pas trouver un moyen de les utiliser pour calculer cela plus efficacement. Toute aide est appréciée, merci!

Answer 1

Puisque vous utilisez python, sklearn a plusieurs implémentations pour trouver plus proches voisins: http://scikit-learn.org/stable/modules/neighbors.html

Il est fourni KDTree et Balltree. Comme pour KDTree, le point principal est de pousser toutes les particules que vous avez dans KDTree, et ensuite pour chaque requête de particule: "donnez-moi toutes les particules dans la plage X". KDtree le fait généralement plus vite que la recherche par bruteforce. Vous pouvez lire plus par exemple ici: https://www.cs.cmu.edu/~ckingsf/bioinfo-lectures/kdtrees.pdf

Si vous utilisez un espace 2D ou 3D, l'autre option consiste simplement à couper l'espace en grande grille (quelle taille de cellule du rayon de masquage) et affecter chaque particule dans une grille cellule. Ensuite, vous pouvez trouver des candidats possibles pour l'interaction en vérifiant simplement les cellules voisines (mais vous devez aussi faire un contrôle à distance, mais pour beaucoup moins de paires de particules).

Answer 2

Voici une technique assez simple utilisant Python simple qui peut réduire le nombre de comparaisons nécessaires.

Nous commençons par trier les points le long de l'axe X, Y ou Z (sélectionné par axis dans le code ci-dessous). Disons que nous choisissons l'axe X. Ensuite, nous bouclons sur des paires de points comme le fait votre code, mais lorsque nous trouvons une paire dont la distance est supérieure à masking_radius, nous vérifions si la différence de leurs coordonnées X est également supérieure à masking_radius. Si c'est le cas, alors nous pouvons sortir de la boucle interne j car tous les points avec un j plus grand ont une plus grande coordonnée X. La fonction dist2 calcule la distance au carré. C'est plus rapide que de calculer la distance réelle car le calcul de la racine carrée est relativement lent.

J'ai également inclus un code qui se comporte comme votre code, c'est-à-dire, il teste chaque paire de points, à des fins de comparaison de vitesse; il sert également à vérifier que le code rapide est correct.;)

from random import seed, uniform 
from operator import itemgetter 

seed(42) 

# Make some fake data 
def make_point(hi=10.0): 
    return [uniform(-hi, hi) for _ in range(3)] 

psize = 1000 
points = [make_point() for _ in range(psize)] 

masking_radius = 4.0 
masking_radius2 = masking_radius ** 2 

def dist2(p, q): 
    return (p[0] - q[0])**2 + (p[1] - q[1])**2 + (p[2] - q[2])**2 

pair_count = 0 
test_count = 0 

do_fast = 1 
if do_fast: 
    # Sort the points on one axis 
    axis = 0 
    points.sort(key=itemgetter(axis)) 

    # Fast 
    for i, p in enumerate(points): 
     left, right = i - 1, i + 1 
     for j in range(i + 1, psize): 
      test_count += 1 
      q = points[j] 
      if dist2(p, q) < masking_radius2: 
       #interaction_pairs.append((i, j)) 
       pair_count += 1 
      elif q[axis] - p[axis] >= masking_radius: 
       break 

     if i % 100 == 0: 
      print('\r {:3} '.format(i), flush=True, end='') 

    total_pairs = psize * (psize - 1) // 2 
    print('\r {}/{} tests'.format(test_count, total_pairs)) 

else: 
    # Slow 
    for i, p in enumerate(points): 
     for j in range(i+1, psize): 
      q = points[j] 
      if dist2(p, q) < masking_radius2: 
       #interaction_pairs.append((i, j)) 
       pair_count += 1 

     if i % 100 == 0: 
      print('\r {:3} '.format(i), flush=True, end='') 

print('\n', pair_count, 'pairs') 

sortie avec do_fast = 1

181937/499500 tests 

13295 pairs 

sortie avec do_fast = 0

13295 pairs 

Bien sûr, si la plupart des paires de points sont à masking_radius les uns des autres, il y a ne sera pas beaucoup d'avantages à utiliser cette technologie ique. Et le tri des points ajoute un peu de temps, mais TimSort de Python est plutôt efficace, surtout si les données sont déjà partiellement triées, donc si le masking_radius est suffisamment petit, vous devriez voir une amélioration notable de la vitesse.