効率よく粒子の密度を求める方法

Question

###与えられた多数の粒子の位置情報から粒子密度を求めたい。
粒子の位置情報を乱数で与え、そこからグリッドで分割された(一次元)空間における粒子の密度を推定するプログラムを作成しています。
一つの考えとしてfor文を使えば密度を推定することができるのですが、計算速度があまり速くありません。for文を使わずに密度を推定する方法があるでしょうか？

###for文を用いた方法

import numpy as np
import matplotlib.pyplot as plt
nop=100   #粒子数
nx =10    #空間グリッド数
dx =1     #グリッド幅
lx =nx*dx #空間長
xx =np.arange(nx) #空間をグリッドに分割
xp =np.random.rand(nop)*lx #粒子位置を乱数で生成
ds =np.zeros(nx) #密度、0クリア

for ip in range(nop):
    ix=int(xp[ip]/dx) #粒子をグリッド状に分配
    ds[ix]=ds[ix]+1   #グリッド上の粒子を数える

plt.subplot(2,1,1); plt.plot(xp,np.ones(nop),'.')
plt.subplot(2,1,2); plt.plot(xx,ds,'-o')
plt.show()

この方法を使って粒子数を10**6個にしたところ計算時間が3.05秒でした。

###for文を使わない方法(失敗)
計算時間がかかっているfor文を使わずに上記のプログラムと同じことをしたいと思っています。そこで以下のようなプログラムを作製したのですが、うまく行きませんでした。なにかいい方法があるのでしょうか？それともこの場合はfor文を使わざるおえないのでしょうか？

import numpy as np
import matplotlib.pyplot as plt
nop=100   #粒子数
nx =10    #空間グリッド数
dx =1     #グリッド幅
lx =nx*dx #空間長
xx =np.arange(nx) #空間をグリッドに分割
xp =np.random.rand(nop)*lx #粒子位置を乱数で生成
ds =np.zeros(nx) #密度、0クリア
#for ip in range(nop):
ix=np.int_(xp/dx) #粒子をグリッド状に分配
ds[ix]=ds[ix]+1   #グリッド上の粒子を数える

plt.subplot(2,1,1); plt.plot(xp,np.ones(nop),'.')
plt.subplot(2,1,2); plt.plot(xx,ds,'-o')
plt.show()

Answer 1

単純に numpy.histogram()を使うとよいのではないでしょうか。

https://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram.html

ds = np.histogram(xp/dx, bins=np.arange(nx+1))[0]

---

【補足】

コメントにて依頼があったので、他の方法を何点か

ds =np.zeros(nx)
for i in np.arange(nx):
    d = xp/dx
    ds[i] = ((i <= d) & (d < i+1)).sum()

ds =np.zeros(nx)
for i in np.arange(nx):
    ds[i] = ((xp/dx).astype(int) == i).sum()

ds = np.bincount((xp/dx).astype(int))

Comments

[Masaru14]

ありがとうございます！計算速度が10倍ほど改善されました。
後々、線形補間や境界条件を加えたいと考えており、外部関数を使うのを控えて柔軟なプログラムを作りたいのですが、配列操作だけで上記のプログラムを作るのは不可能でしょうか？

[magichan]

> 配列操作だけで
何処までを許容しているのかがよくわからないのですが（そもそもデータがnumpy配列に入っているので）、

とりあえずデータ全てをループで回すのではなくグリッド毎にループを回すようにすると、それほど遅くなることはないかと思います。
一応、回等欄に幾つか方法を追加したので確認ください。

[Masaru14]

ありがとうございます。勉強させていただきます。。