Cython

2017-02-12 3 views
2

を使用してPythonの機能のパフォーマンスの向上を目指すCython

私はCythonと私のPythonプログラムを高速化しようとしています。私が書いているコードは、隠れマルコフモデル(HMM)の長いシーケンスの確率を再帰的かつ効率的に計算するために使用されるForwardアルゴリズムの試みです。この問題は、通常、評価問題と呼ばれます。

if __name__=='__main__': 
    A=[[0.7,0.3],[0.4,0.6]] 
    A= numpy.matrix(A) 
    A=pandas.DataFrame(A,columns=['H','C'],index=['H','C']) 
    ''' 
    three types of emmission, small s, medium m and large l 
    ''' 
    B=[[0.1,0.4,0.5],[0.7,0.2,0.1]] 
    B=numpy.matrix(B) 
    B=pandas.DataFrame(B,columns=['S','M','L'],index=['H','C']) 
    ''' 
    initial probabilities for state, H and C 
    ''' 
    pi=[0.6,0.4] 
    pi=numpy.matrix(pi) 
    pi=pandas.DataFrame(pi,columns=['H','C']) 
    O=(0,1,0,2) 
    O=('S','M','S','L') 
    X=('H','H','C','C') 
    H=HMM(A,B,pi,O,X) 
    print H.evaluate() 

%timeitを使用している場合、私は純粋なのpython

でこの出力を得る:ファイルで

Pythonのコード

は、このコードは、で呼び出すことができますhmm.py

import numpy 
import pandas 

class HMM(): 
    '''  
    args: 
     O: 
      observation sequence. A list of 'H's or 'T's 

     X: 
      state sequence. 'S','M' or 'L's 

     A: 
      transition matrix, N by N 

     B: 
      Emission matrix, M by N 

     M: 
      Number of possibilities in emission matrix 

     pi: 
      initial transition matrix 

     N: 
      Number of states 

     T: 
      length of the observation sequence 

     Q: 
      possible hidden states (Xs) 

     V: 
      possible observations (Os) 

    ''' 
    def __init__(self,A,B,pi,O,X): 
     self.A=A 
     self.N=self.A.shape[0] 
     self.B=B 
     self.M=self.B.shape[1] 
     self.pi=pi 
     self.O=O 
     self.T=len(O) 
     self.Q=list(self.A.index) 
     self.V=list(self.B.keys()) 
     self.X=X 


    def evaluate(self): 
     ''' 
     Solve the evaluation problem for HMMs 
     by implementing the forward algorithm 
     ''' 
     c0=0 
     ct=numpy.zeros(self.T) 
     alpha= numpy.zeros((self.T,self.N)) 

     ## compute alpha[0] 
     for i in range(self.N): 
      pi0=self.pi[self.Q[i]] 
      bi0=self.B.loc[self.Q[i]][self.O[0]] 
      alpha[0,i]=pi0*bi0 
      c0+=alpha[0,i] 
      ct[0]=alpha[0,i] 
     ct[0]=1/ct[0]#[0] 
     alpha=alpha*ct[0] 
     for t in range(1,self.T): 
      for i in range(self.N): 
       for j in range(self.N): 
        aji= self.A[self.Q[j]][self.Q[i]] 
        alpha[t,j]= alpha[t-1,j]*aji 
       ct[t]=ct[t]+alpha[t,i] 
      ct[t]=1/ct[t] 
      alpha=alpha*ct[t] 
     return (alpha,ct) 

と呼ばれます

enter image description here

Iは、新しいファイルhmm.pyx拡張子で(むしろ全クラスより)evaluate機能を配置Cython

とコンパイル:setup.py

import numpy 
cimport numpy 

cpdef evaluate_compiled(A,B,pi,O,X): 
    ''' 
    Solve the evaluation problem for HMMs 
    by implementing the forward algorithm 
    ''' 
    T=len(O) 
    N=len(list(set(X))) 
    Q=list(set(X)) 
    V=list(set(O)) 

    c0=0 
    ct=numpy.zeros(T) 
    alpha= numpy.zeros((T,N)) 

    ## compute alpha[0] 
    for i in range(N): 
     pi0=pi[Q[i]] 
     bi0=B.loc[Q[i]][O[0]] 
     alpha[0,i]=pi0*bi0 
     c0+=alpha[0,i] 
     ct[0]=alpha[0,i] 
    ct[0]=1/ct[0]#[0] 
    alpha=alpha*ct[0] 
    for t in range(1,T): 
     for i in range(N): 
      for j in range(N): 
       aji= A[Q[j]][Q[i]] 
       alpha[t,j]= alpha[t-1,j]*aji 
      ct[t]=ct[t]+alpha[t,i] 
     ct[t]=1/ct[t] 
     alpha=alpha*ct[t] 
    return (alpha,ct) 

try: 
    from setuptools import setup 
    from setuptools import Extension 
except ImportError: 
    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('hmm', 
          sources=['HMMCluster/hmm.pyx'], #File is in a directory called HMMCluster 
            include_dirs=[numpy.get_include()])] ) 

コンパイル後、使用:

from hmm import evaluate_compiled 

そして、私が使用することができます上記__main__ブロックの下で、evaluateのインプレース:

print evaluate_compiled(A,B,pi,O,X) 

%timeitを持つ:

enter image description here

あなたは変更せず、見ることができるように私はコードの速度が3倍に向上しています。しかし、私が読んでいるすべてのドキュメントでは、Pythonのスピードの欠如は、変数型を動的に推測することが原因であることを示唆しています。したがって、原則として、私は可変型を宣言し、事態をさらに加速させることができます。

型宣言

とCythonは今、この記事の最後の機能は、再び同じアルゴリズムであるが、型宣言私が手コンパイルおよび%はtimeit後

cpdef evaluate_compiled_with_type_declaration(A,B,pi,O,X): 
    cdef int t,i,j 
    cdef int T = len(O) 
    cdef int N = len(list(set(X))) 
    cdef list Q = list(set(X)) 
    cdef list V = list(set(O)) 
    cdef float c0 = 0 
# cdef numpy.ndarray ct = numpy.zeros(T,dtype=double) ## this caused compilation to fail 
    ct=numpy.zeros(T) 
    alpha= numpy.zeros((T,N)) 
    for i in range(N): 
     pi0=pi[Q[i]] 
     bi0=B.loc[Q[i]][O[0]] 
     alpha[0,i]=pi0*bi0 
     c0+=alpha[0,i] 
     ct[0]=alpha[0,i] 
    ct[0]=1/ct[0]#[0] 
    alpha=alpha*ct[0] 
    for t in range(1,T): 
     for i in range(N): 
      for j in range(N): 
       aji= A[Q[j]][Q[i]] 
       alpha[t,j]= alpha[t-1,j]*aji 
      ct[t]=ct[t]+alpha[t,i] 
     ct[t]=1/ct[t] 
     alpha=alpha*ct[t] 
    return (alpha,ct) 

と:

enter image description here

ご覧のとおり、型宣言は行われていませんコードのパフォーマンスをさらに向上させます。私の質問です:このコードの速度をさらに向上させることができますか?もしそうなら、どうすればいいですか?コメント欄で

編集 次の提案は、私は、追加型宣言を追加:

cdef float pi0,bi0 
cdef numpy.ndarray[numpy.float64_t, ndim=1] ct 
cdef numpy.ndarray[numpy.float64_t, ndim=2] aij,alpha 

をして%timeitからこれを得た:

enter image description here

だから、もう一度、まだ本当のスピードアップにも宣言型を使用します。 (ちょうど種類や形状のために推測)

cpdef evaluate_compiled_with_type_declaration(double[:,:] A, double[:,:] B, double[:] pi, int[:] O, int[:] X): 

:Cythonでnumpyのアレイを使用する

+0

により示唆されるように、ファイルにcython -aを使用し、あなたのコードは、 "コンパイル" されているかどうかを確認するには、 https://wiki.python.org/moin/PythonSpeed/PerformanceTipsで説明されているように、「ドットを避ける」ことができます。しかし、私はCPythonの根底にあるCコードを分析することは決して難しくありません –

+2

'aji'、' alpha'と 'ct'の型を定義しようとします。例えば' cdef numpy.ndarray [numpy.float64_t、 ndim = 1] ... 'となります。 –

+2

'-a'オプションを使って[cython'コマンドを実行する](http://docs.cython.org/en/latest/src/reference/compilation.html)、作成されたHTML文書を開きます。可能であれば、ループの外側で 'pandas'と' ndarray'のインデックスを移動する必要があります。 –

答えて

3

"現代" の方法は、対応する引数として宣言されなければならない "型付きMemoryviews" http://docs.cython.org/en/latest/src/userguide/memoryviews.html

です。

彼らはあなたの結果の配列を宣言し、一度にそれらを埋めることができA[i,j]なくA[i][j]

として直接インデックスを作成する必要があります。

cdef double[:] ct=numpy.zeros(T) 
    cdef double[:,:] alpha= numpy.zeros((T,N)) 

最適化のために、セットが順序を失うようコンパイラディレクティブがhttp://cython.readthedocs.io/en/latest/src/reference/compilation.html?highlight=cdivision#compiler-directives(特にcdivisionboundscheck

あなたはlist(set(...))経由で入力されたデータを変換するべきではありません参照してください。

ライン `BI0 = B.loc [Q [I]] [O [0]]`でウォーレンWeckesser