2012-04-08 74 views
2

我目前正在将一些C代码翻译为Python。此代码正用于帮助识别射电天文学中使用的CLEAN算法产生的错误。为了进行这种分析,必须在特定像素值(由ANT_pix给出)中找到强度映射的傅里叶变换的值,Q斯托克斯映射和U斯托克斯映射。这些地图只有257 * 257个阵列。Python中的DFT花费的时间明显长于C

下面的代码需要几秒钟的时间才能用C运行,但需要花费数小时才能运行Python。我非常肯定它的优化非常好,因为我对Python的知识很差。

感谢您提供任何帮助。

更新我的问题是如果有更好的方法来实现Python中的循环,这将加快速度。我已经阅读了很多关于Python的其他问题的答案,如果可能的话,建议避免嵌套for循环,我只是想知道是否有人知道一个好的方式来实现像下面的Python代码没有循环或更好优化的循环。我意识到这可能是一个很高的命令!

我到目前为止一直在使用FFT,但是我的主管想看看DFT会产生什么样的差异。这是因为天线位置通常不会以精确的像素值出现。使用FFT需要四舍五入到最接近的像素值。

我使用Python作为CASA,用来降低射电天文数据集的计算机程序是用Python编写,并在其中执行的Python脚本远比C.

原始代码

def DFT_Vis(ANT_Pix="",IMap="",QMap="",UMap="", NMap="", Nvis=""): 

UV=numpy.zeros([Nvis,6]) 
Offset=(NMap+1)/2 
ANT=ANT_Pix+Offset; 

i=0 
l=0 
k=0 
SumI=0 
SumRL=0 
SumLR=0 


z=0 

RL=QMap+1j*UMap 
LR=QMap-1j*UMap 

Factor=[math.e**(-2j*math.pi*z/NMap) for z in range(NMap)] 

for i in range(Nvis): 
    X=ANT[i,0] 
    Y=ANT[i,1] 

    for l in range(NMap): 

     for k in range(NMap): 

      Temp=Factor[int((X*l)%NMap)]*Factor[int((Y*k)%NMap)]; 

      SumI+=IMap[l,k]*Temp 
      SumRL+=RL[l,k]*Temp 
      SumLR+=IMap[l,k]*Temp    


     k=1 

    UV[i,0]=SumI.real 
    UV[i,1]=SumI.imag 
    UV[i,2]=SumRL.real 
    UV[i,3]=SumRL.imag 
    UV[i,4]=SumLR.real 
    UV[i,5]=SumLR.imag 
    l=1 
    k=1 
    SumI=0 
    SumRL=0 
    SumLR=0 

return(UV) 
容易得多
+0

几小时是夸张还是真的需要几个小时才能运行? – 2012-04-08 14:54:26

+1

您能否将您的问题改为一个问题? – 2012-04-08 14:55:25

+0

“下面的代码需要几秒钟才能运行C”。呃,代码是Python代码。我们看不到C代码。但手卷纯Python DFT代码可能会很慢,这应该不会让人感到意外。为什么你认为你应该重新实现这个代码,而不是使用其中一个已经很好的已知库? – 2012-04-08 15:07:07

回答

1

如果你有兴趣在提高你的脚本0123的性能可能是一个选项。

1

我不是FFT方面的专家,但我的理解是FFT只是计算DFT的一种快速方法。所以对我来说,你的问题听起来像是你正在尝试写一个冒泡排序算法来看看它是否比quicksort有更好的答案。他们都是排序算法,会给出相同的结果!

所以我质疑你的基本前提。我想知道你是否可以改变你对数据的四舍五入并从SciPy FFT代码中获得相同的结果。

另外,根据我的DSP教科书,由于浮点运算不精确,FFT可以产生比计算DFT更准确的答案,并且FFT在调查中寻找更少的浮点运算正确答案。

如果你有一些可以执行计算的工作C代码,你可以总是包装C代码让你从Python调用它。这里讨论:Wrapping a C library in Python: C, Cython or ctypes?

为了回答您的实际问题:如@ ZoZo123指出,这将是一个巨大的胜利改变从range()xrange()。使用range()时,Python必须建立一个数字列表,然后在完成时销毁列表;与xrange() Python只是一个迭代器,一次产生一个数字。 (但请注意,在Python 3中。x,range()做了一个迭代器,并且没有xrange()。)

另外,如果这段代码不需要与代码的其余部分集成,你可以尝试在PyPy下运行这段代码。这正是PyPy可以最佳优化的那种代码。 PyPy的问题是,你的项目目前必须是“纯粹的”Python,而且看起来你正在使用NumPy。 (有一些项目可以让NumPy和PyPy一起工作,但这还没有完成。)http://pypy.org/

如果这段代码确实需要与其他代码集成,那么我认为您需要看看Cython(as由@KrzysztofRosiński指出)。

+0

实际上,PyPy已经支持很多基本的Numpy功能,所以这个代码应该在'import numpypy'和'import numpy'后面工作(参见http://buildbot.pypy.org/numpy-status/latest.html)。 – TryPyPy 2012-04-08 20:20:40

+1

@TryPyPy,我的理解是PyPy对NumPy的支持仍然是实验性的和不断发展的。所以如果他有大量的SciPy代码,我怀疑他可以把它全部放到PyPy中;但他应该能够运行所示的示例代码。因此,我对“如果这段代码不必与其他代码集成在一起”发表评论。我期待PyPy运行SciPy的所有前景,并且所有科学家默认使用PyPy;但我们还没有。 – steveha 2012-04-08 20:34:18