5

まず、Windows XP マシンで python[2.7.2]、numpy[1.6.2rc1]、cython[0.16]、gcc[MinGW] コンパイラを使用しています。

numpy 配列に格納された 3D バイナリ データ (つまり、1 と 0) を処理するには、3D 連結コンポーネント アルゴリズムが必要でした。残念ながら、既存のコードを見つけることができなかったので、ここにあるコードを 3D 配列で動作するように調整しました。すべてがうまく機能しますが、巨大なデータセットを処理するには速度が望ましい. その結果、cython に出くわし、試してみることにしました。

これまでのところ、cython は速度を改善しました: Cython: 0.339 s Python: 0.635 s

cProfile を使用して、純粋な python バージョンで時間のかかる行は次のとおりです。

new_region = min(filter(lambda i: i > 0, array_region[xMin:xMax,yMin:yMax,zMin:zMax].ravel()))

質問:行を「cythonize」する正しい方法は何ですか:

new_region = min(filter(lambda i: i > 0, array_region[xMin:xMax,yMin:yMax,zMin:zMax].ravel()))
for x,y,z in zip(ind[0],ind[1],ind[2]):

この作業が他の人に役立つことを願っています。

純粋な python バージョン [*.py]:

import numpy as np

def find_regions_3D(Array):
    x_dim=np.size(Array,0)
    y_dim=np.size(Array,1)
    z_dim=np.size(Array,2)
    regions = {}
    array_region = np.zeros((x_dim,y_dim,z_dim),)
    equivalences = {}
    n_regions = 0
    #first pass. find regions.
    ind=np.where(Array==1)
    for x,y,z in zip(ind[0],ind[1],ind[2]):

        # get the region number from all surrounding cells including diagnols (27) or create new region                        
        xMin=max(x-1,0)
        xMax=min(x+1,x_dim-1)
        yMin=max(y-1,0)
        yMax=min(y+1,y_dim-1)
        zMin=max(z-1,0)
        zMax=min(z+1,z_dim-1)

        max_region=array_region[xMin:xMax+1,yMin:yMax+1,zMin:zMax+1].max()

        if max_region > 0:
            #a neighbour already has a region, new region is the smallest > 0
            new_region = min(filter(lambda i: i > 0, array_region[xMin:xMax+1,yMin:yMax+1,zMin:zMax+1].ravel()))
            #update equivalences
            if max_region > new_region:
                if max_region in equivalences:
                    equivalences[max_region].add(new_region)
                else:
                    equivalences[max_region] = set((new_region, ))
        else:
            n_regions += 1
            new_region = n_regions

        array_region[x,y,z] = new_region


    #Scan Array again, assigning all equivalent regions the same region value.
    for x,y,z in zip(ind[0],ind[1],ind[2]):
        r = array_region[x,y,z]
        while r in equivalences:
            r= min(equivalences[r])
        array_region[x,y,z]=r

    #return list(regions.itervalues())
    return array_region

純粋な python の高速化:

#Original line:
new_region = min(filter(lambda i: i > 0, array_region[xMin:xMax+1,yMin:yMax+1,zMin:zMax+1].ravel()))

#ver A:
new_region = array_region[xMin:xMax+1,yMin:yMax+1,zMin:zMax+1]
min(new_region[new_region>0])

#ver B:
new_region = min( i for i in array_region[xMin:xMax,yMin:yMax,zMin:zMax].ravel() if i>0)

#ver C:
sub=array_region[xMin:xMax,yMin:yMax,zMin:zMax]
nlist=np.where(sub>0)
minList=[]
for x,y,z in zip(nlist[0],nlist[1],nlist[2]):
    minList.append(sub[x,y,z])
new_region=min(minList)

時間の結果:
O: 0.0220445
A: 0.0002161
B: 0.0173195
C: 0.0002560

Cython バージョン [*.pyx]:

import numpy as np
cimport numpy as np

DTYPE = np.int
ctypedef np.int_t DTYPE_t

cdef inline int int_max(int a, int b): return a if a >= b else b
cdef inline int int_min(int a, int b): return a if a <= b else b

def find_regions_3D(np.ndarray Array not None):
    cdef int x_dim=np.size(Array,0)
    cdef int y_dim=np.size(Array,1)
    cdef int z_dim=np.size(Array,2)
    regions = {}
    cdef np.ndarray array_region = np.zeros((x_dim,y_dim,z_dim),dtype=DTYPE)
    equivalences = {}
    cdef int n_regions = 0
    #first pass. find regions.
    ind=np.where(Array==1)
    cdef int xMin, xMax, yMin, yMax, zMin, zMax, max_region, new_region, x, y, z
    for x,y,z in zip(ind[0],ind[1],ind[2]):

        # get the region number from all surrounding cells including diagnols (27) or create new region                        
        xMin=int_max(x-1,0)
        xMax=int_min(x+1,x_dim-1)+1
        yMin=int_max(y-1,0)
        yMax=int_min(y+1,y_dim-1)+1
        zMin=int_max(z-1,0)
        zMax=int_min(z+1,z_dim-1)+1

        max_region=array_region[xMin:xMax,yMin:yMax,zMin:zMax].max()

        if max_region > 0:
            #a neighbour already has a region, new region is the smallest > 0
            new_region = min(filter(lambda i: i > 0, array_region[xMin:xMax,yMin:yMax,zMin:zMax].ravel()))
            #update equivalences
            if max_region > new_region:
                if max_region in equivalences:
                    equivalences[max_region].add(new_region)
                else:
                    equivalences[max_region] = set((new_region, ))
        else:
            n_regions += 1
            new_region = n_regions

        array_region[x,y,z] = new_region


    #Scan Array again, assigning all equivalent regions the same region value.
    cdef int r
    for x,y,z in zip(ind[0],ind[1],ind[2]):
        r = array_region[x,y,z]
        while r in equivalences:
            r= min(equivalences[r])
        array_region[x,y,z]=r

    #return list(regions.itervalues())
    return array_region

Cython の高速化:

使用:

cdef np.ndarray region = np.zeros((3,3,3),dtype=DTYPE)
...
        region=array_region[xMin:xMax,yMin:yMax,zMin:zMax]
        new_region=np.min(region[region>0])

時間: 0.170、オリジナル: 0.339 秒

結果

提供された多くの有用なコメントと回答を検討した後、現在のアルゴリズムは
Cython: 0.0219
Python: 0.4309 で実行されています。

Cython は、純粋な Python の 20 倍の速度を提供しています。

現在の Cython コード:

import numpy as np
import cython
cimport numpy as np
cimport cython

from libcpp.map cimport map

DTYPE = np.int
ctypedef np.int_t DTYPE_t

cdef inline int int_max(int a, int b): return a if a >= b else b
cdef inline int int_min(int a, int b): return a if a <= b else b

@cython.boundscheck(False)
def find_regions_3D(np.ndarray[DTYPE_t,ndim=3] Array not None):
    cdef unsigned int x_dim=np.size(Array,0),y_dim=np.size(Array,1),z_dim=np.size(Array,2)
    regions = {}
    cdef np.ndarray[DTYPE_t,ndim=3] array_region = np.zeros((x_dim,y_dim,z_dim),dtype=DTYPE)
    cdef np.ndarray region = np.zeros((3,3,3),dtype=DTYPE)
    cdef map[int,int] equivalences
    cdef unsigned int n_regions = 0

    #first pass. find regions.
    ind=np.where(Array==1)
    cdef np.ndarray[DTYPE_t,ndim=1] ind_x = ind[0], ind_y = ind[1], ind_z = ind[2]
    cells=range(len(ind_x))
    cdef unsigned int xMin, xMax, yMin, yMax, zMin, zMax, max_region, new_region, x, y, z, i, xi, yi, zi, val
    for i in cells:

        x=ind_x[i]
        y=ind_y[i]
        z=ind_z[i]

        # get the region number from all surrounding cells including diagnols (27) or create new region                        
        xMin=int_max(x-1,0)
        xMax=int_min(x+1,x_dim-1)+1
        yMin=int_max(y-1,0)
        yMax=int_min(y+1,y_dim-1)+1
        zMin=int_max(z-1,0)
        zMax=int_min(z+1,z_dim-1)+1

        max_region = 0
        new_region = 2000000000 # huge number
        for xi in range(xMin, xMax):
            for yi in range(yMin, yMax):
                for zi in range(zMin, zMax):
                    val = array_region[xi,yi,zi]
                    if val > max_region: # val is the new maximum
                        max_region = val

                    if 0 < val < new_region: # val is the new minimum
                        new_region = val

        if max_region > 0:
           if max_region > new_region:
                if equivalences.count(max_region) == 0 or new_region < equivalences[max_region]:
                    equivalences[max_region] = new_region
        else:
           n_regions += 1
           new_region = n_regions

        array_region[x,y,z] = new_region


    #Scan Array again, assigning all equivalent regions the same region value.
    cdef int r
    for i in cells:
        x=ind_x[i]
        y=ind_y[i]
        z=ind_z[i]

        r = array_region[x,y,z]
        while equivalences.count(r) > 0:
            r= equivalences[r]
        array_region[x,y,z]=r

    return array_region

セットアップファイル [setup.py]

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy

setup(
    cmdclass = {'build_ext': build_ext},
    ext_modules = [Extension("ConnectComp", ["ConnectedComponents.pyx"],
                             include_dirs =[numpy.get_include()],
                             language="c++",
                             )]
)

ビルド コマンド:

python setup.py build_ext --inplace
4

3 に答える 3

7

@gotgenes が指摘しているように、必ず を使用しcython -a <file>、表示される黄色の量を減らそうとする必要があります。黄色は生成された C の悪化と悪化に対応します。

黄色の量を減らすことがわかったもの：

  1. これは、入力が 3 次元である限り、範囲外の配列アクセスがまったくない状況のように見えるArrayため、範囲チェックをオフにすることができます。

    cimport cython

@cython.boundscheck(False)
def find_regions_3d(...):

  2. コンパイラに効率的なインデックス作成のためのより多くの情報を提供します。つまりcdef、ndarrayできる限り多くの情報を提供するときはいつでも:

     def find_regions_3D(np.ndarray[DTYPE_t,ndim=3] Array not None):
     [...]
     cdef np.ndarray[DTYPE_t,ndim=3] array_region = ...
     [etc.]

  3. ポジティブ/ネガティブネスに関する詳細情報をコンパイラーに提供します。つまり、特定の変数が常に正になることがわかっている場合、cdefそれはunsigned intではなく となりintます。これは、Cython が負のインデックス チェックを排除できることを意味します。

  4. indタプルをすぐに解凍します。つまり、

    ind = np.where(Array==1)
cdef np.ndarray[DTYPE_t,ndim=1] ind_x = ind[0], ind_y = ind[1], ind_z = ind[2]

  5. for x,y,z in zip(..[0],..[1],..[2])コンストラクトの使用は避けてください。どちらの場合も、次のように置き換えます。

    cdef int i
for i in range(len(ind_x)):
    x = ind_x[i]
    y = ind_y[i]
    z = ind_z[i]

  6. 派手なインデックス作成/スライスを行うことは避けてください。特に、2 回行うことは避けてください。filter!の使用は避けてください。すなわち交換

    max_region=array_region[xMin:xMax,yMin:yMax,zMin:zMax].max()
if max_region > 0:
    new_region = min(filter(lambda i: i > 0, array_region[xMin:xMax,yMin:yMax,zMin:zMax].ravel()))
    if max_region > new_region:
        if max_region in equivalences:
            equivalences[max_region].add(new_region)
        else:
            equivalences[max_region] = set((new_region, ))

    より冗長に

    max_region = 0
new_region = 2000000000 # "infinity"
for xi in range(xMin, xMax):
    for yi in range(yMin, yMax):
        for zi in range(zMin, zMax):
            val = array_region[xi,yi,zi]
            if val > max_region: # val is the new maximum
                max_region = val

            if 0 < val < new_region: # val is the new minimum
                new_region = val

if max_region > 0:
   if max_region > new_region:
       if max_region in equivalences:
           equivalences[max_region].add(new_region)
       else:
           equivalences[max_region] = set((new_region, ))
else:
   n_regions += 1
   new_region = n_regions

    これはあまり見栄えがよくありませんが、トリプル ループは約 10 行ほどの C にコンパイルされますが、オリジナルのコンパイルされたバージョンは数百行の長さで、多くの Python オブジェクト操作があります。

    (明らかにcdef、使用するすべての変数、特にこのコードではxi、yi、ziおよびが必要です。)val

  7. セットで行うことは最小要素を見つけることだけなので、すべての等価物を保存する必要はありません。したがって、代わりにへのequivalencesマッピングintがintある場合は、置き換えることができます

    if max_region in equivalences:
    equivalences[max_region].add(new_region)
else:
    equivalences[max_region] = set((new_region, ))

[...]

while r in equivalences:
    r = min(equivalences[r])

    と

    if max_region not in equivalences or new_region < equivalences[max_region]:
    equivalences[max_region] = new_region

[...]

while r in equivalences:
    r = equivalences[r]

  8. 最終的に行うべき最後のことは、Python オブジェクトをまったく使用しないことです。具体的には、equivalences. intにマッピングされているintため、これは簡単from libcpp.map cimport mapにcdef map[int,int] equivalencesなりまし.. not in equivalencesた。(その場合、C++ コンパイラが必要になることに注意してください。)equivalences.count(..) == 0.. in equivalencesequivalences.count(..) > 0

于 2012-08-24T17:14:19.320 に答える
3

（他の人が読みやすいように、上記のコメントからコピーされました）

私はscipyのndimage.labelがあなたが望むことをすると信じています(私はあなたのコードに対してそれをテストしませんでしたが、それは非常に効率的であるはずです)。明示的にインポートする必要があることに注意してください。

from scipy import ndimage 
ndimage.label(your_data, connectivity_struct)

その後、他の組み込み関数を適用できます（境界矩形、重心などを見つけるなど）

于 2012-10-02T09:53:21.993 に答える
0

cython 用に最適化するときは、ループ内で、オーバーヘッドの高い Python オブジェクトではなく、ほとんどのネイティブ C データ型が使用されるようにする必要があります。そのような場所を見つける最善の方法は、生成された C コードを見て、多くの Py* 関数呼び出しに変換された行を探すことです。これらの場所は通常、python オブジェクトの代わりに cdef 変数を使用して最適化できます。

あなたのコードでは、たとえば、ループzipが多くのpythonオブジェクトを生成しint、要素を取得するために使用されるインデックスを使用して反復する方がはるかに高速であると思われind[0]ます....しかし、生成されたCコードを見て、不必要に多くの python 関数を呼び出しているように見えるものを参照してください。

于 2012-08-24T15:16:23.383 に答える