3

OpenGL シミュレーションの実行速度を上げようとしているので、この関数を最適化する必要があります。Parakeetを使用したいのですが、そのためには以下のコードをどのように変更する必要があるのか​​ よくわかりません。私が何をすべきかわかりますか？

def distanceMatrix(self,x,y,z):
    " ""Computes distances between all particles and places the result in a matrix such that the ij th matrix entry corresponds to the distance between particle i and j"" "
    xtemp = tile(x,(self.N,1))
    dx = xtemp - xtemp.T
    ytemp = tile(y,(self.N,1))
    dy = ytemp - ytemp.T
    ztemp = tile(z,(self.N,1))
    dz = ztemp - ztemp.T

    # Particles 'feel' each other across the periodic boundaries
    if self.periodicX:
        dx[dx>self.L/2]=dx[dx > self.L/2]-self.L
        dx[dx<-self.L/2]=dx[dx < -self.L/2]+self.L
    if self.periodicY:
        dy[dy>self.L/2]=dy[dy>self.L/2]-self.L
        dy[dy<-self.L/2]=dy[dy<-self.L/2]+self.L
    if self.periodicZ:
        dz[dz>self.L/2]=dz[dz>self.L/2]-self.L
        dz[dz<-self.L/2]=dz[dz<-self.L/2]+self.L

    # Total Distances
    d = sqrt(dx**2+dy**2+dz**2)

    # Mark zero entries with negative 1 to avoid divergences
    d[d==0] = -1

    return d, dx, dy, dz

私が知る限り、Parakeet は上記の関数を変更せずに使用できるはずです。使用するのは Numpy と数学だけです。しかし、Parakeet jit ラッパーから関数を呼び出すと、常に次のエラーが発生します。

AssertionError: Unsupported function: <bound method Particles.distanceMatrix of <particles.Particles instance at 0x04CD8E90>>
4

1 に答える 1

4

Parakeet はまだ新しく、NumPy のサポートは不完全であり、コードはまだ機能していないいくつかの機能に触れています。

1）メソッドをラップしていますが、これまでのところParakeetは関数の処理方法しか知りません。一般的な回避策は、@jit でラップされたヘルパー関数を作成し、必要なすべてのメンバー データを使用してメソッドを呼び出すことです。メソッドが機能しない理由は、意味のある型を「自己」に割り当てるのは自明ではないからです。不可能ではありませんが、難しいので、下にぶら下がっている果物が摘み取られるまで、メソッドはインコに入りません。ぶら下がりフルーツといえば…

2) ブール索引付け。まだ実装されていませんが、次のリリースで実装される予定です。

3) np.tile : これも機能しません。おそらく次のリリースにもあるでしょう。どの組み込み関数と NumPy ライブラリ関数が機能するかを知りたい場合は、Parakeet のmappingsモジュールを見てください。

コードを書き直して、Parakeet に少し親しみやすくしました。

@jit 
def parakeet_dist(x, y, z, L, periodicX, periodicY, periodicZ):
  # perform all-pairs computations more explicitly 
  # instead of tile + broadcasting
  def periodic_diff(x1, x2, periodic):
    diff = x1 - x2 
    if periodic:
      if diff > (L / 2): diff -= L
      if diff < (-L/2): diff += L
    return diff
  dx = np.array([[periodic_diff(x1, x2, periodicX) for x1 in x] for x2 in x])
  dy = np.array([[periodic_diff(y1, y2, periodicY) for y1 in y] for y2 in y])
  dz = np.array([[periodic_diff(z1, z2, periodicZ) for z1 in z] for z2 in z])
  d= np.sqrt(dx**2 + dy**2 + dz**2)

  # since we can't yet use boolean indexing for masking out zero distances
  # have to fall back on explicit loops instead 
  for i in xrange(len(x)):
    for j in xrange(len(x)):
      if d[i,j] == 0: d[i,j] = -1 
  return d, dx, dy, dz

私のマシンでは、N = 2000 の場合、これは NumPy よりも ~3 倍速く実行されます (NumPy の場合は 0.39 秒、Parakeet の場合は 0.14 秒)。ループをより明示的に使用するように配列トラバーサルを書き直すと、パフォーマンスは NumPy より最大で最大 6 倍高速になります (Parakeet は最大で 0.06 秒で実行されます)。

@jit 
def loopy_dist(x, y, z, L, periodicX, periodicY, periodicZ):
  N = len(x)
  dx = np.zeros((N,N))
  dy = np.zeros( (N,N) )
  dz = np.zeros( (N,N) )
  d = np.zeros( (N,N) )

  def periodic_diff(x1, x2, periodic):
    diff = x1 - x2 
    if periodic:
      if diff > (L / 2): diff -= L
      if diff < (-L/2): diff += L
    return diff

  for i in xrange(N):
    for j in xrange(N):
      dx[i,j] = periodic_diff(x[j], x[i], periodicX)
      dy[i,j] = periodic_diff(y[j], y[i], periodicY)
      dz[i,j] = periodic_diff(z[j], z[i], periodicZ)
      d[i,j] = dx[i,j] ** 2 + dy[i,j] ** 2 + dz[i,j] ** 2 
      if d[i,j] == 0: d[i,j] = -1
      else: d[i,j] = np.sqrt(d[i,j])
  return d, dx, dy, dz

少し創造的に書き直すことで、上記のコードを Numba で実行することもできますが、NumPy の約 1.5 倍 (0.25 秒) しか速くなりません。コンパイル時間は、内包表記付きの Parakeet: 1 秒、ループ付きの Parakeet: 0.5 秒、ループ付きの Numba: 0.9 秒でした。

次のいくつかのリリースで、NumPy ライブラリ関数をより慣用的に使用できるようになることを願っていますが、今のところ、内包表記またはループが多くの場合に使用されます。

于 2013-11-24T19:34:47.640 に答える