Python优化问题？

Question

好吧，我最近做了这个作业（不用担心，我已经完成了它，但在C++中），但我很好奇我怎样才能在python中做到这一点。问题在于约2个发光的光源。我不会详细介绍。

下面的代码（即我已经设法优化后半部分有点）：

import math, array 
import numpy as np 
from PIL import Image 

size = (800,800) 
width, height = size 

s1x = width * 1./8 
s1y = height * 1./8 
s2x = width * 7./8 
s2y = height * 7./8 

r,g,b = (255,255,255) 
arr = np.zeros((width,height,3)) 
hy = math.hypot 
print 'computing distances (%s by %s)'%size, 
for i in xrange(width): 
    if i%(width/10)==0: 
     print i,  
    if i%20==0: 
     print '.', 
    for j in xrange(height): 
     d1 = hy(i-s1x,j-s1y) 
     d2 = hy(i-s2x,j-s2y) 
     arr[i][j] = abs(d1-d2) 
print '' 

arr2 = np.zeros((width,height,3),dtype="uint8")   
for ld in [200,116,100,84,68,52,36,20,8,4,2]: 
    print 'now computing image for ld = '+str(ld) 
    arr2 *= 0 
    arr2 += abs(arr%ld-ld/2)*(r,g,b)/(ld/2) 
    print 'saving image...' 
    ar2img = Image.fromarray(arr2) 
    ar2img.save('ld'+str(ld).rjust(4,'0')+'.png') 
    print 'saved as ld'+str(ld).rjust(4,'0')+'.png' 

我设法优化大部分，但仍然与部分巨大的性能差距2 for-s，我似乎无法想出一种绕过使用常见阵列操作的方法......我愿意接受建议：D

编辑： 为了回应弗拉德的建议，我们会发布问题的详细信息： 有2个光源，每个都以正弦波发光： E1 = E0 * sin（omega1 * time + phi01） E2 = E0 * sin（omega2 * time + phi02） 我们考虑omega1 = omega2 = omega = 2 * pi/T和phi01 = phi02 = phi0为简单起见 考虑到x1是从平面上一点的第一个源的距离，该点的光的强度为Ep1 = E0 * sin（ω*时间-2 * PI * x1 /λ+φ0） 其中 λ=速度光* T（振荡周期） 考虑到平面上的两个光源，公式变为 Ep = 2 * E0 * cos（PI *（x2-x1）/ lambda）sin（ω时间-PI *（x2 （x2-x1）/ lambda =（2 * k）* PI/2时，我们可以发现光的强度是最大的和最小值时 （X2-X1）/λ=（2 * K + 1）* PI/2 和之间，其中k是给定的坐标的整数

对于给定的时间瞬间，而变化在光源，并为已知的lambda和E0，我们不得不做一个程序，以了解灯光的外观 恕我直言，我认为我优化了这个问题尽可能多的，因为它可以做...

Answer 1

干扰模式很有趣，不是吗？

所以，首先这将是微不足道的，因为在我的笔记本电脑上运行这个程序只需要十二点半秒。

但让我们看看可以做些什么关于通过numpy数组操作做第一位，我们？我们基本上是你想要的：

arr[i][j] = abs(hypot(i-s1x,j-s1y) - hypot(i-s2x,j-s2y)) 

对于所有i和j。

所以，因为numpy有一个hypot函数可以在numpy数组上工作，所以让我们使用它。我们的第一个挑战是要得到一个正确大小的数组，每个元素等于i，另一个元素等于j。但这并不难;其实，答案在下面的精彩numpy.mgrid我不知道，做点之前我只是这样的：

array_i,array_j = np.mgrid[0:width,0:height] 

有让您(width, height)尺度的阵列为(width,height,3)是兼容的一件小事您的图像生成报表，但是这很容易做到：

arr = (arr * np.ones((3,1,1))).transpose(1,2,0) 

然后我们插入到这个程序，并让事情由数组操作来完成：

import math, array 
import numpy as np 
from PIL import Image 

size = (800,800) 
width, height = size 

s1x = width * 1./8 
s1y = height * 1./8 
s2x = width * 7./8 
s2y = height * 7./8 

r,g,b = (255,255,255) 

array_i,array_j = np.mgrid[0:width,0:height] 

arr = np.abs(np.hypot(array_i-s1x, array_j-s1y) - 
      np.hypot(array_i-s2x, array_j-s2y)) 

arr = (arr * np.ones((3,1,1))).transpose(1,2,0) 

arr2 = np.zeros((width,height,3),dtype="uint8") 
for ld in [200,116,100,84,68,52,36,20,8,4,2]: 
    print 'now computing image for ld = '+str(ld) 
    # Rest as before 

而新的时间是...... 8.2秒。所以你可以节省四秒钟。另一方面，这几乎完全在图像生成阶段，所以也许你可以通过只生成你想要的图像来加强它们。

Answer 2

唯一的变化，在我看来，是将一些操作移出循环：

for i in xrange(width): 
    if i%(width/10)==0: 
     print i,  
    if i%20==0: 
     print '.', 
    arri = arr[i] 
    is1x = i - s1x 
    is2x = i - s2x 
    for j in xrange(height): 
     d1 = hy(is1x,j-s1y) 
     d2 = hy(is2x,j-s2y) 
     arri[j] = abs(d1-d2) 

改进if任何，虽然可能会很小。

Answer 3

列表解析比循环快得多。例如，而不是

for j in xrange(height): 
     d1 = hy(i-s1x,j-s1y) 
     d2 = hy(i-s2x,j-s2y) 
     arr[i][j] = abs(d1-d2) 

你会写

arr[i] = [abs(hy(i-s1x,j-s1y) - hy(i-s2x,j-s2y)) for j in xrange(height)] 

在另一方面，如果你真的试图“优化”，那么你可能想用C重新实现这个算法，并使用SWIG或类似的东西从Python中调用它。

Answer 4

如果使用数组操作而不是循环，它要快得多。对我而言，图像生成现在需要很长时间。相反，你的两个i,j循环的，我有这样的：

I,J = np.mgrid[0:width,0:height] 
D1 = np.hypot(I - s1x, J - s1y) 
D2 = np.hypot(I - s2x, J - s2y) 

arr = np.abs(D1-D2) 
# triplicate into 3 layers 
arr = np.array((arr, arr, arr)).transpose(1,2,0) 
# .. continue program 

要记住未来的基础是：这是不关于优化;在numpy中使用数组形式就像它应该使用它一样。有了经验，你未来的项目不应该绕过python循环，数组形式应该是自然形式。

我们在这里做的很简单。而不是math.hypot我们发现了numpy.hypot并使用它。像所有这些numpy函数一样，它接受ndarrays作为参数，并且完全符合我们的要求。