[Numpy-discussion] Right way of looping throught ndarrays using Cython

Gael Varoquaux gael.varoquaux@normalesup....
Fri Jun 20 14:13:56 CDT 2008


Oups, forgot the attachement.

On Fri, Jun 20, 2008 at 09:13:34PM +0200, Gael Varoquaux wrote:
> On Fri, Jun 20, 2008 at 02:07:08PM -0500, Robert Kern wrote:
> > > Does somebody have a example of fast looping through ndarrays using
> > > modern Cython idioms?

> > If you're using normal Python indexing, then that's where all your
> > time is being spent. You need to grab the actual .data pointer and do
> > C indexing to get speed. Can you show us the code you are timing?

> That is indeed what I was thinking. There is no way of doing this appart
> by using the point-style indexing? I guess I am trying to find the best
> (as in most readable) way of doing this. This is for teaching, not
> production, so I am very much interested in having something simple.

> I am attaching my test file.

> Cheers,

> Gaël
-------------- next part --------------
from numpy import zeros, mgrid

# Make sure numpy is initialized.
include "c_numpy.pxd"

##############################################################################
cdef int inner_loop(float c_x, float c_y):
    cdef float x, y, x_buffer 
    x = 0; y = 0
    cdef int i
    for i in range(50):
        x_buffer = x*x - y*y + c_x
        y = 2*x*y + c_y
        x = x_buffer
        if (x*x + x*y > 100):
            return 50 - i

def do_Mandelbrot_cython():
    x, y = mgrid[-1.5:0.5:500j, -1:1:500j]
    threshold_time = zeros((500, 500))
    cdef int i, j
    for i in range(500):
        for j in range(500):
            threshold_time[i, j] = inner_loop(x[i, j], y[i, j])
    return threshold_time

# %timeit on do_Mandelbrot_cython, on epsilon, gives 761ms per loop

def main():
    threshold_time = do_Mandelbrot_cython()
    from pylab import imshow, cm, clf, show
    clf()
    imshow(threshold_time, cmap=cm.spectral, extent=(-1, 1, -1, 1))
    show()





More information about the Numpy-discussion mailing list