[Numpy-discussion] Rookie problems - Why is C-code much faster?

[Mads Ipsen]

On Tue, 21 Feb 2006, Tim Hochberg wrote:

> Can you explain a little more about what you are trying to calculate?
> The bit about subtracting off  box*rint(dx/box) is a little odd. It
> almost seems like you should be able to do something with fmod, but I
> admit that I'm not sure how.
>
> If I had to guess as to source of the relative slowness I'd say it's
> because you are creating a lot of temporary matrices. There are ways to
> avoid this, but when taken to the extreme, they make your code look
> ugly. You might try the following, untested, code or some variation and
> see if it speeds things up. This makes extensive use of the little known
> optional destination argument for ufuncs. I only tend to do this sort of
> stuff where it's very critical since, as you can see, it makes things
> quite ugly.
>
>      dx_space = x.copy()
>      dy_space = y.copy()
>      scratch_space = x.copy()
>      for i in range(n-1):
>         dx = dx_space[i+1:]
>         dy = dy_space[i+1:]
>         scratch = scratch_space[i+1:]
>         subtract(x[i+1:], x[i], dx)
>         subtract(y[i+1:], y[i], dy)
>         #   dx -= box*rint(dx/box)
>         divide(dx, box, scratch)
>         rint(scratch, scratch)
>         scratch *= box
>         dx -= scratch
>         # dy -= box*rint(dy/box)
>         divide(dy, box, scratch)
>         rint(scratch, scratch)
>         scratch *= box
>         dy -= scratch
>         r2 = dx**2 + dy**2      # square of dist. between points
>
>
>
> Hope that helps:
>
> -tim

Here's what I am trying to do:

My system consists of N particles, whose coordinates in the xy-plane is
given by the two vectors x and y. I need to calculate the distance between
all particle pairs, which goes like this:

I pick particle 1 and calculate its distance to the N-1 other points. Then
I pick particle 2. Since its distance to particle 1 was found in the
previuos step, I only have to find its distance to the N-2 remaining
points. In the i'th step, I therefore only have to consider particle i+1
up to particle N. That explains the loop structure, where

  dx = x[i+1:] - x[i]
  dy = y[i+1:] - y[i]

the resulting vectors dx and dy will contain the x-distances from x[i] to
the proceeding points from i+1 up to N. The square of the distance r2 is
the given by

  r2 = dx**2 + dy**2

Another approach would be to use

  dx = subtract.outer(x,x)
  dy = subtract.outer(y,y)

but that will be overkill, since all distances are counted twice, and
also, the storage requirements grow rapidly if you have more than 1000
particles (app. 10^6 particle pairs).

Thanks for your code feedback, which I'll have a closer look at.  But I
try to believe, that numpy/Numeric/Python was invented with the one
purpose of avoiding coding like this - I think this is also a point you
already made. But thanks again.

// Mads