CythonでScipy関数をコンパイル

Question

私はPython 3.4でシミュレーションを実行しています。これは、疎な配列（csr形式）と密なベクトルの間に多くのドット積を含みます。私は疎な行列のためにScipyを使用しています。他のものはnumpyです。

Cythonを使用すると、私がすべてを適切にcdefし、Pythonのやりとりを最小限に抑えた後（Cytが私に与えられたコードを修正するhtmlファイルを通って）大量のboost（〜6倍のスピード）が得られました。

ここでコードをプロファイルし、シミュレーション時間の50％がドットプロダクトのラインに費やされています。何とかこのラインを加速することが可能かどうか、Cythonのこの1ドット機能を遵守することによって、私は疑問に思っていますか？

私は（csr sprase 2d matrix）ドット（高密度ベクトル）のための私自身の実装を書くことができますが、私はそれを回避しようとしています。

を編集してください。最小限のコード例を示しました。申し訳ありませんが、私はそれを小さくする方法を見ることができません。統計力学の教科書です。鍋の一つが容量を超えるまで、大理石を鍋に入れる。次に、（ここではまばらな）行列に従って伝播するカスケードを開始します。私はバッチサンプリングを使用しています。

最後に向かってください。

from __future__ import division 
import numpy   as np 
import cython 
cimport numpy   as np 
cimport cpython.array 
import scipy.sparse as sps 


@cython.cdivision(True) 
@cython.nonecheck(False) 
@cython.boundscheck(False) 
@cython.wraparound(False) 
def simulate(long[:] capacity_vec, 
      int random_array_size, 
      long n, 
      int seed, 
      int[:] A_col, 
      int[:] A_row, 
      long[:] A_data): 


    #### Initialise #################################################################################################### 

    # Initialise states 
    cdef int[:] damage = np.random.randint(0, int(np.min(capacity_vec)/2), n).astype(np.int32) 
    cdef int[:] dr_list = np.random.choice(n, 1000).astype(np.int32) 
    cdef int[:] states = np.zeros(n).astype(np.int32) 
    cdef int[:] states_ = np.zeros(n).astype(np.int32) 
    cdef int[:] change = np.zeros(n).astype(np.int32) 

    # Initialise counters 
    cdef int k, violations, violations_, counter= 0, dr_id=0, increment_index = 0 


    # Build Sparse Adjecency Matrix 
    cA_sps = sps.csr_matrix((A_data, (A_row, A_col)), shape=(n,n)).astype(np.int32) 


    while counter < 1000: 

     #### Place damage until a cascade starts ####################################################################### 
     while damage[increment_index] <= capacity_vec[increment_index]:# Check for violations 

      increment_index   = dr_list[dr_id]     # Where do we place the marble? 

      damage[increment_index] = damage[increment_index] + 1  # place the marble 

      dr_id     = dr_id + 1      # another random number used 

      if dr_id == random_array_size - 1:       # Check if we run out of random numbers 

       dr_list = np.random.choice(n, random_array_size).astype(np.int32) # if so, pick new increment_index 

       dr_id = 0           # and reset the counter 


     #### Initialise cascade ######################################################################################## 
     violations, violations_ = 1, 0 
     states[increment_index] = 1 


     #### Propagate cascade ######################################################################################### 
     while violations > violations_:        # check for fixed point, propagate cascade 
      for k in range(n): change[k] = states[k] - states_[k] 
      ### THIS LINE IS THE PROBLEM. It takes up half of all simulation time. 
      damage  = damage + cA_sps.dot(change).astype(np.int32) # spread violations 

      states_  = states.copy()        # store previous states 

      # Determine previous and current violations 
      violations, violations_ = 0 , violations 

      for k in range(n): 

       states_[k] = 0 

       if damage[k] > capacity_vec[k]: 

        violations = violations + 1 

        states[k] = 1         # deactivate any node that has a violation 


     for k in range(n): damage[k] = 0 
     counter = counter + 1           # progress cascade id after storing 

Answer 1

私は独自の行列乗算を書いてはいけません。 SciPyは、彼らが何をしているかを知っている賢い人々によって行われます。あなたの自信が数値計算でなければ、そうしないでください。 SciPyコードのほとんどは既にコンパイルされています。

しかし、あなたが見ているのはsparse.csr_matrix.dotのコードです。直接hereの定義を取得し、次にhereを取得すると、Scipyでチェックがほとんど行われていないことがわかります。あなたが望む正確なフォームが分かっているなら、独自の方法を書いて（あなたのSciPyコピーを修正して）、あなたの製品を直接計算することができます。しかし、それがどれくらい助けになるかはわかりません。

あなたはscipyのダウンロードを自分で構築したい場合、それはGitHugからプロジェクト全体をチェックアウトして、より直接的な手順について

python setup.py build 
python setup.py install 

を実行しているのと同じくらい簡単ですがbuild documentationをご確認ください。