1

タイトルがあいまいなのかもしれません。

numpy と matplotlib を使用して、Python でパーティクル シミュレーターを少し使いこなしました。coloumb、重力、風を実装できました。温度と圧力を追加したいだけですが、最適化前の質問があります (すべての悪の根源)。 ）。パーティクルがいつクラッシュするかを確認したい:

Q:ブール条件に基づいて、配列と独自の要素のそれぞれの違いを取得することは可能ですか? ループは避けたい。

例：(x - any element in x) < a 次のようなものを返す必要があります

[True, True, False, True]

x の要素 0、1、3 が条件を満たす場合。

編集：

同等のループは次のようになります。

for i in len(x):
  for j in in len(x):
    #!= not so important
    ##earlier question I asked lets me figure that one out
    if i!=j: 
      if x[j] - x[i] < a:
       True

numpy 操作は if テストよりもはるかに高速であることに気付きました。これにより、物事が大幅に高速化されました。

誰かがそれを試してみたい場合は、ここにサンプルコードがあります。

#Simple circular box simulator, part of part_sim
#Restructure to import into gravity() or coloumb () or wind() or pressure()
#Or to use all forces: sim_full()
#Note: Implement crashing as backbone to all forces

import numpy as np
import matplotlib.pyplot as plt

N = 1000        #Number of particles
R = 8000        #Radius of box          
r = np.random.randint(0,R/2,2*N).reshape(N,2)
v = np.random.randint(-200,200,r.shape)
v_limit = 10000 #Speedlimit


plt.ion()
line, = plt.plot([],'o')
plt.axis([-10000,10000,-10000,10000])

while True:
    r_hit = np.sqrt(np.sum(r**2,axis=1))>R   #Who let the dogs out, who, who?
    r_nhit = ~r_hit                     
    N_rhit = r_hit[r_hit].shape[0]
    r[r_hit] = r[r_hit] - 0.1*v[r_hit]       #Get the dogs back inside
    r[r_nhit] = r[r_nhit] +0.1*v[r_nhit]
    #Dogs should turn tail before they crash!
    #---
    #---crash code here....
    #---crash end
    #---
    vmin, vmax = np.min(v), np.max(v)        
    #Give the particles a random kick when they hit the wall
    v[r_hit]  = -v[r_hit] + np.random.randint(vmin, vmax, (N_rhit,2))
    #Slow down honey
    v_abs = np.abs(v) > v_limit
    #Hit the wall at too high v honey? You are getting a speed reduction
    v[v_abs] *=0.5
    line.set_ydata(r[:,1])
    line.set_xdata(r[:,0])
    plt.draw()

大きなボックスで高速粒子を簡単に区別できる方法がわかったら、上記のデータポイントに色を追加する予定です。

4

1 に答える 1

5

例: x - x < a 内の任意の要素のようなものを返す必要があります

[真、真、偽、真]

x の要素 0、1、3 が条件を満たす場合。numpy 操作は if テストよりもはるかに高速であることに気付きました。これにより、物事が大幅に高速化されました。

はい、それだけm < aです。例えば：

>>> m = np.array((1, 3, 10, 5))
>>> a = 6
>>> m2 = m < a
>>> m2
array([ True,  True, False,  True], dtype=bool)

さて、質問に：

Q:ブール条件に基づいて、配列と独自の要素のそれぞれの違いを取得することは可能ですか? ループは避けたい。

ここで何を求めているのかわかりませんが、そのすぐ下の例と一致していないようです。たとえば、述語を満たす各要素から 1 を減算しようとしていますか? その場合、ブール配列を減算するだけであるという事実False==0に頼ることができます。True==1

>>> m3 = m - m2
>>> m3
>>> array([ 0,  2, 10,  4])

あなたの明確化から、この疑似コードループに相当するものが必要です：

for i in len(x):
  for j in in len(x):
    #!= not so important
    ##earlier question I asked lets me figure that one out
    if i!=j: 
      if x[j] - x[i] < a:
        True

ここでの混乱は、これがあなたが言ったことと正反対であるということだと思います:「ブール条件に基づく独自の要素のそれぞれを持つ配列の違い」は必要ありませんが、「違いに基づくブール条件」は必要ありません独自の要素を持つ配列の ". それでも実際には len(m)*len(m) ブール値の正方行列しか得られませんが、残っている部分は「any」だと思います。

とにかく、m の各要素を m の各要素と比較して、暗黙的なデカルト積を求めています。

これを 2 つのループから 1 つに簡単に減らすことができます (または、暗黙的にそのうちの 1 つをベクトル化し、通常の numpy パフォーマンスの利点を得ることができます)。値ごとに、各要素からその値を減算し、結果を と比較して新しい配列を作成し、aそれらを結合します。

>>> a = -2
>>> comparisons = np.array([m - x < a for x in m])
>>> flattened = np.any(comparisons, 0)
>>> flattened
array([ True,  True, False,  True], dtype=bool)

しかし、これを単純な行列演算に簡単に変換することもできます。mの他のすべての要素から のすべての要素を引くのmは、 だけm - m.Tです。(積をより明示的にすることはできますが、numpy行ベクトルと列ベクトルの加算を処理する方法では必要ありません。)そして、そのすべての要素を scalar と比較しa、 で還元anyすれば完了です。

>>> a = -2
>>> m = np.matrix((1, 3, 10, 5))
>>> subtractions = m - m.T
>>> subtractions
matrix([[ 0,  2,  9,  4],
        [-2,  0,  7,  2],
        [-9, -7,  0, -5],
        [-4, -2,  5,  0]])
>>> comparisons = subtractions < a
>>> comparisons
matrix([[False, False, False, False],
        [False, False, False, False],
        [ True,  True, False,  True],
        [ True, False, False, False]], dtype=bool)
>>> np.any(comparisons, 0)
matrix([[ True,  True, False,  True]], dtype=bool)

または、すべてを 1 行にまとめると、次のようになります。

>>> np.any((m - m.T) < a, 0)
matrix([[ True,  True,  True,  True]], dtype=bool)

m行列ではなく配列にする必要がある場合は、減算行を に置き換えることができますm - np.matrix(m).T。

より高い次元の場合、実際には配列で作業する必要があります。これは、2D 配列をそれ自体でデカルト積して 4D 配列を取得しようとしており、numpy4D 行列を実行しないためです。したがって、単純な「行ベクトル - 列ベクトル = 行列」のトリックは使用できません。ただし、手動で行うことができます：

>>> m = np.array([[1,2], [3,4]]) # 2x2
>>> m4d = m.reshape(1, 1, 2, 2)  # 1x1x2x2
>>> m4d
array([[[[1, 2],
         [3, 4]]]])
>>> mt4d = m4d.T # 2x2x1x1
>>> mt4d
array([[[[1]],
        [[3]]],
       [[[2]],
        [[4]]]])
>>> subtractions = m - mt4d # 2x2x2x2
>>> subtractions
array([[[[ 0,  1],
         [ 2,  3]],
        [[-2, -1],
         [ 0,  1]]],
       [[[-1,  0],
         [ 1,  2]],
        [[-3, -2],
         [-1,  0]]]])

そして、そこから先は以前と同じです。それを 1 行にまとめると、次のようになります。

>>> np.any((m - m.reshape(1, 1, 2, 2).T) < a, 0)

(私の元の答えを覚えているなら、私はどういうわけか空白にして、1 の列ベクトルをreshape掛けることで同じことをしていmましたが、これは明らかにはるかに愚かな方法です。)

最後の簡単な考え: あなたのアルゴリズムが本当に「(の任意の要素yに対して) の各要素に対するbool の結果」である場合m、実際には「任意の要素に対して」は必要なく、 「最大要素に対して」を使用することができます。 . したがって、O(N^2) から O(N) に単純化できます。x - y < axmyy

>>> (m - m.max()) < a

または、aが正の場合、それは常に false であるため、O(1) に単純化できます。

>>> np.zeros(m.shape, dtype=bool)

しかし、あなたの実際のアルゴリズムは実際にはabs(x - y)、またはこの方法では単純化できないもっと複雑なものを使用していると思います。

于 2012-11-29T23:28:47.450 に答える