[Numpy-discussion] computing average distance

[Gael Varoquaux]

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
improved.

First, you can use numpy.indices instead of numpy.mgrid for what you want
to do (see
http://docs.scipy.org/doc/numpy/reference/generated/numpy.indices.html#numpy.indices
):

    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
Python:

 * 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 )

HTH,

Gaël