2

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.

2 Answers 2

2

The issue you are facing is not specific to numpy itself, or advanced indexing in numpy, or whether it creates copies or not. Instead it is driven entirely by ambiguities around whether:

  • indexing is guaranteed to return values that are somehow "inside" the container (its not); and whether

  • in-place add is guaranteed to return a modified version of the original value (its not)

Consider the expression:

x[a] += 100

where x is a list say. The result of x[a] is a value y that has no explicit knowledge of the list it happened to be contained in, and the expression y += 100 is not guaranteed to modify the original value of y ... thus we can never ensure that an expression of the form x.__getitem__(x).__iadd__(100) influences the original x.

Thus the expression x[a] += 100 must be evaluated by the compiler using the steps:

  1. y = x[a]
  2. y += 100
  3. x[a] = y

Or in your case of advanced indexing, we could expect a[ a > 5 ] += 42 to be implemented as:

  1. b = a > 5
  2. c = a[b]
  3. c += 42
  4. a[b] = c

This can be demonstrated by looking at the disassembly of an example function

def f(a):
    a[a > 5] += 42

then

>>> dis.dis(f)
 0 LOAD_FAST                0 (a)
 2 LOAD_FAST                0 (a)
 4 LOAD_CONST               1 (5)
 6 COMPARE_OP               4 (>)  # 1. b = a > 5
 8 DUP_TOP_TWO
10 BINARY_SUBSCR                   # 2. c = a[b]
12 LOAD_CONST               2 (42)
14 INPLACE_ADD                     # 3. c += 42
16 ROT_THREE
18 STORE_SUBSCR                    # 4. a[b] = c
20 LOAD_CONST               0 (None)
22 RETURN_VALUE

Note that in this implementation, the indexing a > 5 is not performed twice. However if you must implement as a chained set of methods, then you will have no choice but to implement it as you suggested.

2
  • Yep, that is slightly more high tech than my answer. Nov 1, 2019 at 11:08
  • Right, so it's as I feared: the compiler can do a shorthand, but as you pointed out the actual data model itself does not make any guarantees for optimality. However, I'm going to update my question with an observation because the matter clearly has subtlety.
    – MB.
    Nov 1, 2019 at 12:08
0

Let us instrument the relevant dunders to see what numpy/python does:

import numpy as np

class spyarray(np.ndarray):
    def __getitem__(self, *args):
        print("__getitem__",self,*args)
        return np.ndarray.__getitem__(self, *args)
    def __setitem__(self, *args):
        print("__setitem__",self,*args)
        return np.ndarray.__setitem__(self, *args)
    def __add__(self, *args):
        print("__add__",self,*args)
        return np.ndarray.__add__(self, *args)
    def __iadd__(self, *args):
        print("__iadd__",self,*args)
        return np.ndarray.__iadd__(self, *args)
    def __repr__(self):
        return np.ndarray.__repr__(self.view(np.ndarray))
    def __str__(self):
        return np.ndarray.__str__(self.view(np.ndarray))


a = np.arange(10).view(spyarray)
a[a>5] += 42

Prints:

__getitem__ [0 1 2 3 4 5 6 7 8 9] [False False False False False False  True  True  True  True]
__iadd__ [6 7 8 9] 42
__setitem__ [0 1 2 3 4 5 6 7 8 9] [False False False False False False  True  True  True  True] [48 49 50 51]

So that looks pretty much like what you came up with.

Your Answer

Reminder: Answers generated by Artificial Intelligence tools are not allowed on Stack Overflow. Learn more

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.