On 06/06/2012 05:06 PM, Nathaniel Smith wrote: > On Wed, Jun 6, 2012 at 9:48 AM, John Salvatier > <jsalv...@u.washington.edu> wrote: >> Hello, >> >> I've noticed that If you try to increment elements of an array with advanced >> indexing, repeated indexes don't get repeatedly incremented. For example: >> >> In [30]: x = zeros(5) >> >> In [31]: idx = array([1,1,1,3,4]) >> >> In [32]: x[idx] += [2,4,8,10,30] >> >> In [33]: x >> Out[33]: array([ 0., 8., 0., 10., 30.]) >> >> I would intuitively expect the output to be array([0,14, 0,10,30]) since >> index 1 is incremented by 2+4+8=14, but instead it seems to only increment >> by 8. What is numpy actually doing here? >> >> The authors of Theano noticed this behavior a while ago so they python loop >> through the values in idx (this kind of calculation is necessary for >> calculating gradients), but this is a bit slow for my purposes, so I'd like >> to figure out how to get the behavior I expected, but faster. >> >> I'm also not sure how to navigate the numpy codebase, where would I look for >> the code responsible for this behavior? > > Strictly speaking, it isn't actually in the numpy codebase at all -- > what's happening is that the Python interpreter sees this code: > > x[idx] += vals > > and then it translates it into this code before running it: > > tmp = x.__getitem__(idx) > tmp = tmp.__iadd__(vals) > x.__setitem__(idx, tmp) > > So you can find the implementations of the ndarray methods > __getitem__, __iadd__, __setitem__ (they're called > array_subscript_nice, array_inplace_add, and array_ass_sub in the C > code), but there's no way to fix them so that this works the way you > want it to, because there's no way for __iadd__ to know that the > temporary values that it's working with are really duplicate copies of > "the same" value in the original array. > > It would be nice if numpy had some sort of standard API for doing what > you want, but not sure what a good API would look like, and someone > would have to implement it.
This operation is also heavily used for the finite element assembling, and a similar question has been raised already several times (e.g. http://old.nabble.com/How-to-assemble-large-sparse-matrices-effectively-td33833855.html). So why not adding a function np.assemble()? r. _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion