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

编辑： 针对 Vlad 的建议，我将发布问题的详细信息： 有 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(omega*time - 2*PI*x1/lambda + phi0) 在哪里 λ = 光速 * T（振荡周期） 考虑平面上的两个光源，公式变为 Ep = 2*E0*cos(PI*(x2-x1)/lambda)sin(omega时间 - PI*(x2-x1)/lambda + phi0) 从中我们可以看出，当光的强度最大时 (x2-x1)/λ = (2*k) * PI/2 最小时 (x2-x1)/λ = (2*k+1) * PI/2 并且在两者之间变化，其中k是整数

对于给定的时刻，给定光源的坐标，对于已知的 lambda 和 E0，我们必须编写一个程序来绘制灯光的外观 恕我直言，我认为我尽可能地优化了这个问题......

Answer 1

如果你使用数组操作而不是循环，它会快得多。对我来说，图像生成现在需要很长时间。而不是你的两个 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 作为参数，并且完全按照我们的意愿行事。

Answer 2

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

所以，首先这将是次要的，因为在我的笔记本电脑上按原样运行这个程序只需要 12 秒半。

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

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

对于所有i 和j。

所以，既然 numpy 有一个适用于 numpy 数组的 hypot 函数，让我们使用它。我们的第一个挑战是获得一个大小合适的数组，每个元素都等于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 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

我想到的唯一改变是将一些操作移出循环：

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)

如果有的话，改进可能会很小。