7

Cython でコードを最適化しようとしています。FFT を使用せずに、パワー スペクトルを実行しています。これは、クラスで行うように指示されたものです。Cython でコードを書き込もうとしましたが、違いがわかりません。これが私のコードです

#! /usr/bin/env python
# -*- coding: utf8 -*-

from __future__ import division
cimport numpy as np
import numpy as np
cimport cython

@cython.boundscheck(False)
def power_spectrum(time, data, double f_min, double f_max, double df,w=1 ):

    cdef double com,f
    cdef double s,c,sc,cc,ss
    cdef np.ndarray[double, ndim=1] power
    cdef np.ndarray[double, ndim=1] freq

    alfa, beta = [],[] 
    m = np.mean(data)
    data -= m       

    freq = np.arange( f_min,f_max,df )
    for f in freq:
            sft = np.sin(2*np.pi*f*time)
            cft = np.cos(2*np.pi*f*time)
            s   = np.sum( w*data*sft )
            c   = np.sum( w*data*cft )
            ss  = np.sum( w*sft**2  )
            cc  = np.sum( w*cft**2  )
            sc  = np.sum( w*sft*cft )

            alfa.append( ( s*cc-c*sc )/( ss*cc-sc**2 ))
            beta.append( ( c*ss-s*sc )/( ss*cc-sc**2 ))
            com = -(f-f_min)/(f_min-f_max)*100
            print "%0.3f%% complete" %com

    power = np.array(alfa)**2 + np.array(beta)**2
    return freq,power,alfa,beta

時間とデータは numpy.loadtxt を介してロードされ、この関数に送信されます。私がする時

cython -a power_spectrum.pyx

.html ファイルは非常に黄色であるため、あまり効率的ではありません。特に for ループ全体とべき乗の計算とすべてを返します。

Cython の公式ガイドを読んでみましたが、C でコーディングしたことがないので、ややわかりにくいです。

すべてのヘルプは非常に高く評価されています:)

4

1 に答える 1

4

Cython はthis に従ってnumpy 配列を読み取ることができますが、魔法のようにコンパイルすることはできませんnp.sum。まだ numpy メソッドを呼び出しているだけです。

あなたがする必要があるのは、あなたのためにそれをコンパイルできる純粋な cython で内側のループを書き直すことです。np.sumそのため、などを再実装する必要があります。 andをnp.sin事前に割り当てるのは良い考えです。使用しないで、できるだけ多くの変数を試してください。aplfabetaappendcdef

編集

以下は、完全に C コンパイルされた内部ループを示す完全な例です (黄色はありません)。コードが正しいかどうかはわかりませんが、良い出発点になるはずです! 特に、cdefどこでもの使用に注意し、cdivision をオンにし、標準ライブラリのインポートsinとインポートをオンにします。cos

from __future__ import division
cimport numpy as np
import numpy as np
cimport cython
from math import pi

cdef extern from "math.h":
    double cos(double theta)
    double sin(double theta)

@cython.boundscheck(False)
@cython.cdivision(True)
def power_spectrum(np.ndarray[double, ndim=1] time, np.ndarray[double, ndim=1] data, double f_min, double f_max, double df, double w=1 ):

    cdef double com,f
    cdef double s,c,sc,cc,ss,t,d
    cdef double twopi = 6.283185307179586
    cdef np.ndarray[double, ndim=1] power
    cdef np.ndarray[double, ndim=1] freq = np.arange( f_min,f_max,df )
    cdef int n = len(freq)
    cdef np.ndarray[double, ndim=1] alfa = np.zeros(n)
    cdef np.ndarray[double, ndim=1] beta = np.zeros(n)
    cdef int ndata = len(data)
    cdef int i, j

    m = np.mean(data)
    data -= m       

    for i in range(ndata):
        f = freq[i]

        s = 0.0
        c = 0.0
        ss = 0.0
        cc = 0.0
        sc = 0.0
        for j in range(n):
            t = time[j]
            d = data[j]
            sf = sin(twopi*f*t)
            cf = cos(twopi*f*t)
            s += w*d*sf
            c += w*d*cf
            ss += w*sf**2
            cc += w*cf**2
            sc += w*sf*cf

        alfa[i] = ( s*cc-c*sc )/( ss*cc-sc**2 )
        beta[i] = ( c*ss-s*sc )/( ss*cc-sc**2 )

    power = np.array(alfa)**2 + np.array(beta)**2
    return freq,power,alfa,beta
于 2012-05-27T10:49:44.490 に答える