[Numpy-discussion] Rookie problems - Why is C-code much faster?
Mads Ipsen
mpi at osc.kiku.dk
Tue Feb 21 13:05:05 CST 2006
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
More information about the Numpy-discussion
mailing list