[Numpy-discussion] computing average distance

Gael Varoquaux gael.varoquaux@normalesup....
Sun Nov 2 14:41:20 CST 2008

Hey Emmanuelle, ( :>)

On Sun, Nov 02, 2008 at 08:39:39PM +0100, Emmanuelle Gouillart wrote:
> although I'm not an expert either, it seems to me you could improve your
> code a lot by using numpy.mgrid 
> Below is a short example of what you could do
> coordinates = numpy.mgrid[0:R, 0:R, 0:R]
> X, Y, Z = coordinates[0].ravel(), coordinates[1].ravel(),coordinates[2].ravel() bits = self.bits.ravel() 
> distances = numpy.sqrt((X[bits==1]-centre[0])**2 +
> 	(Y[bits==1]-centre[0])**2 + (Z[bits==1]-centre[0])**2)

> There must be a way to do it without flattening the arrays, but I haven't
> found it. Anyway, you can surely do what you want without a loop!

A few cosmetic comments: this code is very good, but it can be slightly

First, you can use numpy.indices instead of numpy.mgrid for what you want
to do (see

    numpy.mgrid[0:R, 0:R, 0:R] == numpy.indices((R, R, R)

Second, I don't like to use explicit indices [0], [1], [2] for x, y, z. It
is so easy to get it wrong. The good news is that you can do better with

 * eg define the center coordinnates as so:

    x0, y0, z0 = numpy.array(scipy.ndimage.measurements.center_of_mass(self.bits))

   this will avoid what seems to be an error in Emmanuelle's answer. By
   the way, I am not too sure why Paul has used a numpy.array here. For
   this part of the code, it is not terribly usefull.

 * Same thing for the call to mgrid:

    X, Y, Z = numpy.indices((R, R, R)

Finally, when you create the raveled bits you are already making a copy,
so you might as well define a variable called 'mask' which is equal to:

    mask = (self.bits.ravel() == 1)

this will make the code more readable.

And lastly, if the shape of bits it (R, R, R), you don't need to ravel.

With all these comments (supposing the shape of bits is (R, R, R), the
final code would look like:

    x0, y0, z0 = scipy.ndimage.measurements.center_of_mass(self.bits)
    X, Y, Z = numpy.indices((R, R, R))
    mask = (self.bits == 1)
    distances = numpy.sqrt(   (X[mask] - x0)**2
			    + (Y[mask] - y0)**2
			    + (Z[mask] - z0)**2 )



More information about the Numpy-discussion mailing list