Python optimization problem?
Ok, I've had this homework lately (don't worry, I've already done this, but in C ++), but I got curious how I can do this in python. The problem is 2 light sources emitting light. I will not go into details.
Here is the code (which I was able to optimize a bit in the second part):
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'
I was able to optimize most of it, but there is still a huge performance gap in the 2-s file part, and I can't imagine how I can get around this using shared array operations ... I'm open to suggestions: D
Edit: In response to Vlad's suggestion, I'll post the details of the problem: There are 2 light sources, each emitting light as a sine wave: E1 = E0 * sin (omega1 * time + phi01) E2 = E0 * sin (omega2 * time + phi02) we consider omega1 = omega2 = omega = 2 * PI / T and phi01 = phi02 = phi0 for simplicity, considering x1 as the distance from the first source of a point on the plane, the light intensity at this point is Ep1 = E0 * sin (omega * time - 2 * PI * x1 / lambda + phi0) where lambda = speed of light * T (oscillation period) Considering both light sources on a plane, the formula becomes Ep = 2 * E0 * cos (PI * (x2-x1) / lambda) sin (omegatime - PI * (x2-x1) / lambda + phi0) and from this we could understand that the light intensity is maximum when (x2-x1) / lambda = (2 * k) * PI / 2 and minimum when (x2-x1 ) / lambda = (2 * k + 1) * PI / 2 and changes between them,where k is an integer
For a given point in time, given the coordinates of the lights and also for the known lambdas and E0s, we needed to make a program to paint what the light looks like IMHO, I think I optimized the problem as much as possible ...
a source to share
Jamming patterns are fun, aren't they?
So, it will be minor at first because it only takes twelve and a half seconds to run this program as it is on my laptop.
But let's see what can be done with executing the first bit through numpy array operations, shall we? Basically we want:
arr[i][j] = abs(hypot(i-s1x,j-s1y) - hypot(i-s2x,j-s2y))
For everyone i
and j
.
So, since numpy has a function hypot
that works on numpy arrays, use this. Our first task is to get an array of the desired size with each element equal i
and another with each element equal j
. But it's not too difficult; in fact, the answer below points to a wonderful numpy.mgrid
one that I didn't know about before I did:
array_i,array_j = np.mgrid[0:width,0:height]
There is a small problem to make your array (width, height)
in (width,height,3)
compatible with your statements about education, but it's pretty easy to do:
arr = (arr * np.ones((3,1,1))).transpose(1,2,0)
Then we connect it to your program and let you do the array operations:
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
And the new time ... 8.2 seconds. This way you save maybe as much as four seconds. On the other hand, this is almost exclusively in the image generation stages, so perhaps you can tighten them only by generating the images you want.
a source to share
If you use array operations instead of loops, it is much faster. For me, image generation now takes so long. Instead of two two loops i,j
, I have the following:
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
The basics you want to remember for the future: this is not about optimization; using array shapes in numpy just uses it as it is supposed to be used. With experience, your future projects shouldn't be traversing python loops; array shapes should be natural.
What we did here was very simple. Instead, math.hypot
we found numpy.hypot
and used it. Like all such numpy functions, it takes ndarrays as arguments and does exactly what we want.
a source to share
List enumerations are much faster than loops. For example, instead of
for j in xrange(height):
d1 = hy(i-s1x,j-s1y)
d2 = hy(i-s2x,j-s2y)
arr[i][j] = abs(d1-d2)
You write
arr[i] = [abs(hy(i-s1x,j-s1y) - hy(i-s2x,j-s2y)) for j in xrange(height)]
On the other hand, if you are really trying to "optimize" then you may need to redefine this algorithm in C and use SWIG or similar to call it from python.
a source to share
The only changes that come to my mind are moving some operations out of the loop:
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)
The improvement, if any, is likely to be minor, though.
a source to share