[Numpy-discussion] Fwd: Advanced selection, duplicate indices, and augmented assignment

[Timothy Hochberg]

On 1/10/07, A. M. Archibald <peridot.faceted at gmail.com> wrote:
>
> On 10/01/07, Charles R Harris <charlesr.harris at gmail.com> 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]


-- 
//=||=\\
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20070111/79973902/attachment.html>