1

3D 密度関数を作成するより良い方法はありますか?

def make_spot_3d(bright, spread, x0,y0,z0):
    # Create x and y indices
    x = np.linspace(-50, 50, 200)
    y = np.linspace(-50, 50, 200)
    z = np.linspace(-50, 50, 200)

    X, Y, Z = np.meshgrid(x, y, z)
    Intensity = np.uint16(bright*np.exp(-((X-x0)/spread)**2
                                        -((Y-y0)/spread)**2
                                        -((Z-z0)/spread)**2))

    return Intensity

この関数は、mayavi でプロットできる 3D numpy 配列を生成できます。ここに画像の説明を入力

ただし、次のように関数を使用してスポット (~100) のクラスターを生成する場合:

Spots = np.asarray([make_spot_3d(100,2, *loc) for loc in locations])
cluster = np.sum(Spots, axis=0)

たとえば、次のようになります。

ここに画像の説明を入力 実行時間は約 1 分 (cpu i5) です。これはもっと速いかもしれません。

4

3 に答える 3

3

明らかな改善は、ブロードキャスティングを使用して、フルではなく「スパース」メッシュで強度関数を評価することですmeshgrid。たとえば、次のようになります。

X, Y, Z = np.meshgrid(x, y, z, sparse=True)

これにより、私のマシンではランタイムが約 4 分の 1 に短縮されます。

%timeit make_spot_3d(1., 1., 0, 0, 0)
1 loops, best of 3: 1.56 s per loop 

%timeit make_spot_3d_ogrid(1., 1., 0, 0, 0)
1 loops, best of 3: 359 ms per loop

位置、広がり、明るさの計算をベクトル化することで、リストの理解に伴うオーバーヘッドを取り除くことができます。次に例を示します。

def make_spots(bright, spread, x0, y0, z0):

    # Create x and y indices
    x = np.linspace(-50, 50, 200)
    y = np.linspace(-50, 50, 200)
    z = np.linspace(-50, 50, 200)

    # this will broadcast out to an (nblobs, ny, nx, nz) array
    dx = x[None, None, :, None] - x0[:, None, None, None]
    dy = y[None, :, None, None] - y0[:, None, None, None]
    dz = z[None, None, None, :] - z0[:, None, None, None]
    spread = spread[:, None, None, None]
    bright = bright[:, None, None, None]

    # we can save time by performing the exponentiation over 2D arrays
    # before broadcasting out to 4D, since exp(a + b) == exp(a) * exp(b)
    s2 = spread * spread
    a = np.exp(-(dx * dx) / s2)
    b = np.exp(-(dy * dy) / s2)
    c = np.exp(-(dz * dz) / s2)

    intensity = bright * a * b * c

    return intensity.astype(np.uint16)

ここbrightでspread、、、、x0およびy0はz01D ベクトルです。これにより配列が生成さ(nblobs, ny, nx, nz)れ、最初の軸で合計できます。生成するブロブの数と、それらを評価するグリッドの大きさによっては、この中間配列を作成すると、メモリの点で非常に高価になる可能性があります。

もう 1 つのオプションは、単一の(ny, nx, nz)出力配列を初期化し、その場で合計を計算することです。

def sum_spots_inplace(bright, spread, x0, y0, z0):

    # Create x and y indices
    x = np.linspace(-50, 50, 200)
    y = np.linspace(-50, 50, 200)
    z = np.linspace(-50, 50, 200)

    dx = x[None, None, :, None] - x0[:, None, None, None]
    dy = y[None, :, None, None] - y0[:, None, None, None]
    dz = z[None, None, None, :] - z0[:, None, None, None]
    spread = spread[:, None, None, None]
    bright = bright[:, None, None, None]

    s2 = spread * spread
    a = np.exp(-(dx * dx) / s2)
    b = np.exp(-(dy * dy) / s2)
    c = np.exp(-(dz * dz) / s2)

    out = np.zeros((200, 200, 200), dtype=np.uint16)

    for ii in xrange(bright.shape[0]):
        out += bright[ii] * a[ii] * b[ii] * c[ii]

    return out

これにより必要なメモリは少なくなりますが、潜在的な欠点は、Python でのループが必要になることです。

相対的なパフォーマンスを把握するには、次のようにします。

def sum_spots_listcomp(bright, spread, x0, y0, z0):
    return np.sum([make_spot_3d(bright[ii], spread[ii], x0[ii], y0[ii], z0[ii])
                   for ii in xrange(len(bright))], axis=0)

def sum_spots_vec(bright, spread, x0, y0, z0):
    return make_spots(bright, spread, x0, y0, z0).sum(0)

# some fake data
bright = np.random.rand(10) * 100
spread = np.random.rand(10) * 100
x0 = (np.random.rand(10) - 0.5) * 50
y0 = (np.random.rand(10) - 0.5) * 50
z0 = (np.random.rand(10) - 0.5) * 50

%timeit sum_spots_listcomp(bright, spread, x0, y0, z0)
# 1 loops, best of 3: 16.6 s per loop

%timeit sum_spots_vec(bright, spread, x0, y0, z0)
# 1 loops, best of 3: 1.03 s per loop

%timeit sum_spots_inplace(bright, spread, x0, y0, z0)
# 1 loops, best of 3: 330 ms per loop
于 2015-07-10T17:06:46.167 に答える
1

i5 プロセッサがあり、スポットが互いに独立しているため、マルチスレッドを実装するとよいでしょう。多くの Numpy 操作が GIL を解放するため、必ずしも複数のプロセスが必要になるわけではありません。追加のコードは非常に単純です。

from multiprocessing.dummy import Pool

if __name__ == '__main__':
    wrap = lambda pos: make_spot_3d(100, 2, *pos)
    cluster = sum(Pool().imap_unordered(wrap, positions))

アップデート

職場の PC でいくつかテストした結果、上記のコードは単純すぎて非効率的であることを認めざるを得ません。8 コアでは、シングルコアのパフォーマンスと比較して、スピードアップはわずか 1.5 倍です。

マルチスレッド化は良いアイデアだと思いますが、成功するかどうかは実装に大きく依存します。

于 2015-07-10T18:49:33.990 に答える
0

つまり、任期中にすべての操作を 800 万 (=200*200*200) 回実行することになります。まず第一に、if の 8 分の 1 を計算し、それをミラーリングするだけで、それを 100 万に減らすことができます (球がたまたまグリッドの中心にある場合)。ミラーリングは無料ではありませんが、exp.

また、intensity の値が 0 に落ちたら、計算を停止する必要がある可能性が非常に高くなります。少し対数マジックを使用すると、200*200*200 グリッドよりもはるかに小さい関心領域を見つけることができます。

于 2015-07-10T17:54:16.000 に答える