効率的なパーティクルペアの相互作用計算

Question

私はシミュレーションで複数のタイムステップのパーティクルポジションのリストを生成するNボディシミュレーションを持っています。与えられたフレームについて、私はdist(p[i], p[j]) < masking_radiusのようなパーティクルのインデックスのペアのリストを生成したいと思っています。(i, j)基本的に私はペアが互いにある距離内にある「相互作用」ペアのリストを作成しています。粒子の数が多いのは、このプロセスは、長い時間（1試験あたり> 1時間）を取り、そしてそれはひどく私が何をする必要があるかに制限されているので

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)) 

：私の現在の実装は次のようになりますデータ。私は、粒子の可能な組み合わせをすべて比較するのではなく、これらのペアを計算する方が効率的であるように、データを構造化する効率的な方法があれば疑問に思っていました。私はKDTreesを調べていましたが、これをより効率的に計算するためにそれらを利用する方法を理解できませんでした。どんな助けもありがとう、ありがとう！

Answer 1

あなたのpythonを使用しているので、sklearnが発見最近傍のための複数の実装があります。 http://scikit-learn.org/stable/modules/neighbors.html

がKDTreeとBalltreeが提供されます。

KDTree用として主なポイントはKDTreeにあなたが持っているすべての粒子をプッシュすることですし、各粒子のクエリを尋ねる：「私の範囲のXのすべての粒子を与えます」。 KDtreeは通常、ブルートフォース検索よりも高速です。 もっと読むことができます：https://www.cs.cmu.edu/~ckingsf/bioinfo-lectures/kdtrees.pdf

2次元または3次元空間を使用している場合、別のオプションは、大きなグリッド（マスキング半径のセルサイズ）にスペースをカットし、各パーティクルを1つのグリッドに割り当てます細胞。次に、隣り合った細胞をチェックするだけで、相互作用の可能性のある候補を見つけることができます（ただし、距離チェックを行う必要がありますが、パーティクルペアはずっと少ない）。

Answer 2

ここでは、単純なPythonを使用して、必要な比較回数を減らすことができる、非常に簡単な手法です。

まずX、Y、または（以下のコードでaxisによって選択された）Z軸のいずれかに沿って点を並べ替えます。 X軸を選択するとしましょう。あなたのコードのようにポイントペアをループしますが、距離がmasking_radiusより大きいペアが見つかった場合は、そのX座標の差がmasking_radiusより大きいかどうかをテストします。そうであれば、より大きなjのすべてのポイントがより大きなX座標を持つため、内部のjループから脱出することができます。

私のdist2関数は、自乗距離を計算します。これは、平方根の計算が比較的遅いため、実際の距離を計算するよりも高速です。

私はまた、すなわち、それはスピードの比較のために、ポイントのすべてのペアをテストし、あなたのコードに似た振る舞いコードが含まれていました。高速コードが正しいことをチェックする役割も果たします。;）do_fast = 0もちろん

13295 pairs 

とdo_fast = 1

181937/499500 tests 

13295 pairs 

出力と

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') 

出力を、点対のほとんどが、互いにmasking_radius内にある場合この技術を使用することで多くのメリットはありませんイケ。ポイントの並べ替えに時間がかかりますが、PythonのTimSortは効率的です（特にデータが部分的にソートされている場合）。masking_radiusが十分小さい場合は速度が大幅に向上するはずです。