I'm implementing a narrow and limited scripting DSL using python and I'd like to be able to functionally do the equivalent of the following:
import numpy as np
a = np.arange(10)
a[ a > 5 ] += 42
=> array([ 0, 1, 2, 3, 4, 5, 48, 49, 50, 51])
The above code works as one would expect. If I start expanding the above code, I get the following first layer of internals:
a[a>5].__iadd__(42)
Which also works as expected. However, I'm unable to find the indexer method that would allow me to operate the __iadd__ on the array itself instead of a copy of the array. As such, not unexpectedly, the following code doesn't do what I want:
import numpy as np
a = np.arange(10)
a.__getitem__(a>5).__iadd__(42)
=> array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
Only if I do:
a.__setitem__(a>5, a.__getitem__(a>5).__iadd__(42))
do I seem to get the behaviour I'm looking for, but at this point this is no longer a proper in-place assignment operator at all and more importantly, I'm indexing twice (once for the read and once for the write).
Numpy's index page seems to imply that advanced indexing (i.e. indexing where the subscript list is an ndarray) always returns a copy. Does this actually mean a[a>5].__iadd__(42)
is in fact always implemented using the fallback method? Is there something I'm missing or is this simply never possible, or at least not possible withouth interpreter magic?
Edit:
So as per @donkopotamus' answer, the data model does not allow us to do this in one shot. This answers the question.
However, numpy
being a vectorized library, the indexing absolutely can't afford to be non-vectorized and executed multiple times.
Here's a "proof":
import cython
import numpy as np
@cython.locals(arr="float[:]",
mask="bint[:]",
val=float,
i=int)
@cython.boundscheck(False)
def func(arr,mask,val):
for i in range(len(mask)):
if mask[i]:
arr[i] += val
This code, when compiled and timed, is slower than numpy in place:
a = np.arange(1e6)
%%timeit
a[a%3==0] += 42
=> 40.5 ms ± 376 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
a = np.arange(1e6)
%%timeit
func(a, (a%3==0), 42)
=> 116 ms ± 2.76 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
So the REPL interpreted statement is running faster than a 3 line cython function which pretty much rips through a memory view as fast as the CPU would allow it.
At this stage, none of it is making any sense anymore. I know numpy is hand crafted to optimize vectorization operations but I'm not understanding how it integrates with the python interpreter in a way that makes sense. Is it caching the BINARY_SUBSCR/STORE_SUBSCR pair?
@donkopotamus please note that while the indexing operation isn't computed twice, in the python code it is interpreted twice in the sense that a mask is performed on the read, and then an entire second scan and mask is performed on the write. In the cython code above, that operation occurs only once for read and write).
Any insight is appreciated.