On 1/10/07, A. M. Archibald <[EMAIL PROTECTED]> wrote:
On 10/01/07, Charles R Harris <[EMAIL PROTECTED]> wrote:
> > You're the second or third person to recently come up with a useful
> application for using inplace ops with repeated indices only to be
stymied
> by that fact that they don't actually work that way. I've been wondering
if
> we should try to come up with a way to support this usage. Although it
might
> be feasible to make += work, the implementation is horrible, and there
are
> all sorts of inefficiencies built into that method of doing it.
> >
> > The natural approach would seem to be to hang another method off the
> ufuncs. That way you get support for all of the various operations. I
was
> thinking something like:
> >
> >
> > unfunc.inplace (dest, indices, source)
>
> Yeah, something like that. A general indirect addressing primitive of
some
> sort would be useful for these situations. The indices are sort of
pointers
> into the destination, the problem is how to bind this to operators. It
> should almost be some sort lvalue analogous to **p = whatever, or in
this
> case dest[index[j]] {assignment} src[j]. So maybe indirect(dest,
indices,
> source, op). Then there is all the madness of fancy indexing (which I
think
> I could have done without). Because the current applications are limited
and
> could be met with the simplest 1D version, it would probably be best to
> start with the most basic functionality and then extend it if the need
> arose. Or maybe even some simple c extension that did the whole basic +=
> thing for one dimension.
I disagree that it should be hung off ufuncs. That doesn't allow, for
example, averaging of normal vectors (without some highly clever and
hard-to-follow axis specification).
What makes you think this? To go back to the start of the thread, the OP
wanted the following to work:
def calculate_normals(indices, vertices):
i = indices
faces = vertices[i].reshape((i.shape[0] / 3, 3, 3)) # indices = 3 x
N?
face_normals = cross(faces[:,1] - faces[:,0], faces[:,2]
-faces[:,0])
normals = zeros(vertices.shape, float32)
normals[i] += face_normals.repeat(3, axis=0)
normals /= sqrt(sum(normals*normals,axis=1)).reshape(normals.shape[0],
1)
return normals
Hanging it off of ufuncs, would just mean that the line:
normals[i] += face_normals.repeat(3, axis=0)
would get replaced by this line:
add.inplace(normals, i, face_normals.repeat(3, axis=0))
I implemented one possible version in pure python, just to think about
the UI (and to see if it could be made fast using sorting; I don't
think it will actually help much, and I also note that whether it even
has any hope depends on whether there are more objects to accumulate
than destination slots, or fewer). The basic design is
def accumulate(dest, indices, values, combine=N.add);
Unless I'm mistaken, semantically, this is no different than add.inplace as
defined above, you've just moved stuff around. In practice, I think that
hanging something like this off of ufuncs is likely to be the way to go to
get decent performance. Plus, that's the traditional way, in numpy, to deal
with methods that are the same except for the operator at the core of them.
I'm not in love with the name in 'inplace' necessarily, but I don't think
accumulate will really fly either, accumulate already has a different
meaning that is close enough to be confusing, but far enough apart that I
don't see a way to shoehorn them together.
[SNIP]
--
//=||=\\
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
