[SciPy-user] Python loops too slow

[Dave Hirschfeld]

Pauli Virtanen <pav <at> iki.fi> writes:
> 
> i = arange(ngrid/2 + 1)
> j = arange(ngrid/2 + 1)
> result = 2*pi*hypot(i[None,:], j[:,None])/reso/ngrid
> 

It seems Pauli's solution may need to be repeated in the final output as below.
Is that what you really want though because the submatrix isn't reflected
uniformly - i.e. the first row and column are missing.

Anyway, timings are below - it seems using meshgrid is more efficient than
hypot.

-Dave

Python 2.5.4 (r254:67916, Dec 23 2008, 15:10:54) [MSC v.1310 32 bit (Intel)]
Type "copyright", "credits" or "license" for more information.

IPython 0.9.1 -- An enhanced Interactive Python.

  Welcome to pylab, a matplotlib-based Python environment.
  For more information, type 'help(pylab)'.

In [2]: def f1(ngrid, reso):
   ...:         result = zeros([ngrid, ngrid])
   ...:     for i in arange((ngrid / 2)+1):
   ...:            for j in arange((ngrid / 2)+1):
   ...:               result[j,i] = 2. * pi * sqrt(i ** 2. + j ** 2) / reso /
(ngrid*1.0)
   ...:     for i in xrange(ngrid):
   ...:            result[i,ngrid / 2+1:] = result[i,1:(ngrid / 2)][::-1]
   ...:     for i in xrange(ngrid):
   ...:            result[ngrid / 2+1:,i] = result[1:(ngrid / 2),i][::-1]
   ...:     return result
   ...: #
   ...:

In [3]: def f2(ngrid, reso):
   ...:         nsub = ngrid/2+1
   ...:     i = np.arange(nsub)
   ...:     j = np.arange(nsub)
   ...:     submatrix = 2*pi*np.hypot(i[None,:], j[:,None])/reso/np.float(ngrid)
   ...:     result = zeros([ngrid, ngrid])
   ...:     result[0:nsub,0:nsub] = submatrix
   ...:     result[0:nsub,nsub:] = np.fliplr(submatrix[:,1:-1])
   ...:     result[nsub:,:] = np.flipud(result[1:nsub-1,:])
   ...:     return result
   ...: #
   ...:

In [4]: def f3(ngrid, reso):
   ...:         nsub = ngrid/2+1
   ...:     i, j = meshgrid(np.arange(nsub),np.arange(nsub))
   ...:     submatrix = 2*pi*np.sqrt(i**2 + j**2)/reso/np.float(ngrid)
   ...:     result = zeros([ngrid, ngrid])
   ...:     result[0:nsub,0:nsub] = submatrix
   ...:     result[0:nsub,nsub:] = np.fliplr(submatrix[:,1:-1])
   ...:     result[nsub:,:] = np.flipud(result[1:nsub-1,:])
   ...:     return result
   ...: #
   ...:

In [5]: np.allclose(f1(100,1),f2(100,1))
Out[5]: True

In [6]: np.allclose(f1(100,1),f3(100,1))
Out[6]: True

In [7]:

In [8]: timeit f1(100,1)
10 loops, best of 3: 60.6 ms per loop

In [9]:

In [10]: timeit f2(100,1)
1000 loops, best of 3: 840 µs per loop

In [11]:

In [12]: timeit f3(100,1)
1000 loops, best of 3: 258 µs per loop

In [13]: