[Numpy-discussion] Indexing NEP draft

Sebastian Berg sebastian at sipsolutions.net
Wed Nov 11 05:02:50 EST 2015

Hi all,

at scipy discussing with Nathaniel and others, we thought that maybe we
can push for orthogonal type indexing into numpy. Now with the new
version out and some other discussions done, I thought it is time to
pick it up :).

The basic ideas are twofold. First make indexing easier and less
confusing for starters (and advanced users also), and second improve
interoperability with projects such as xray for whom orthogonal/outer
type indexing makes more sense.

I have started working on:

1. A preliminary draft of an NEP you can view at
or at the end of this mail.

2. A preliminary implementation of `oindex` attribute with
orthogonal/outer style indexing in
https://github.com/numpy/numpy/pull/6075 which you can try out by
cloning numpy and then running from the source dir:

git fetch upstream pull/6075/head:pr-6075 && git checkout pr-6075;
python runtests.py --ipython

This will fetch my PR, switch to the branch and open an interactive
ipython shell where you will be able to do arr.oindex[].

Note that I consider the NEP quite preliminary in many parts, and it may
still be very confusing unless you are well versed with current advanced
indexing. There are some longer examples comparing the different styles
and another "example" which tries to show a "use case" example going
from simpler to more complex indexing operations.
Any comments are very welcome, and if it is "I don't understand a
word" :). I know it is probably too short and, at least without
examples, not easy to understand.



The current NEP draft:

Implementing intuitive and full featured advanced indexing

:Author: Sebastian Berg
:Date: 2015-08-27
:Status: draft

Executive summary

Advanced indexing with multiple array indices is typically confusing to
both new, and in many cases even old, users of NumPy. To avoid this
and allow for more and clearer features, we propose to:

1. Introduce ``arr.oindex[indices]`` which allows advanced indices, but
   uses outer indexing logic.
2. Introduce ``arr.vindex[indices]`` which use the current
   "vectorized"/broadcasted logic but with two differences from
   fancy indexing:
   1. Boolean indices always use the outer indexing logic.
      (Multi dimensional booleans should be allowed).
   2. The integer index result dimensions are always the first axes
      of the result array. No transpose is done, even for a single
      integer array index.

3. Vanilla indexing on the array will only give warnings and eventually
   errors either:
   * when there is ambiguity between legacy fancy and outer indexing
     (note that ``arr[[1, 2], :, 0]`` is such a case, an integer
     can be the "second" integer index array),
   * when any integer index array is present (possibly additional for
     more then one boolean index array).

These constraints are sufficient for making indexing generally
with expectations and providing a less surprising learning curve with

Note that all things mentioned here apply both for assignment as well as

Understanding these details is *not* easy. The `Examples` section gives
examples. And the hopefully easier `Motivational Example` provides some
motivational use-cases for the general ideas and is likely a good start
anyone not intimately familiar with advanced indexing.


Old style advanced indexing with multiple array (boolean or integer)
also called "fancy indexing", tends to be very confusing for new users.
While fancy (or legacy) indexing is useful in many cases one would
assume that the result of multiple 1-d ranges is analogous to multiple
slices along each dimension (also called "outer indexing").

However, legacy fancy indexing with multiple arrays broadcasts these
into a single index over multiple dimensions. There are three main
of confusion when multiple array indices are involved:

1. Most new users will usually expect outer indexing (consistent with
   slicing). This is also the most common way of handling this in other
   packages or languages.
2. The axes introduced by the array indices are at the front, unless
   all array indices are consecutive, in which case one can deduce where
   the user "expects" them to be:

   * `arr[:, [0, 1], :, [0, 1]]` will have the first dimension shaped 2.
   * `arr[:, [0, 1], [0, 1]]` will have the second dimension shaped 2.

3. When a boolean array index is mixed with another boolean or integer
   array, the result is very hard to understand (the boolean array is
   converted to integer array indices and then broadcast), and hardly
   There is no well defined broadcast for booleans, so that boolean
   indices are logically always "``outer``" type indices.

Proposed rules

From the three problems noted above some expectations for NumPy can
be deduced:

1. There should be a prominent outer/orthogonal indexing method such as
2. Considering how confusing fancy indexing can be, it should only
   occur explicitly (e.g. ``arr.vindex[indices]``)
3. A new ``arr.vindex[indices]`` method, would not be tied to the
   confusing transpose rules of fancy indexing (which is for example
   needed for the simple case of a single advanced index). Thus, it
   no transposing should be done. The axes of the advanced indices are
   always inserted at the front, even for a single index.
4. Boolean indexing is conceptionally outer indexing. A broadcasting
   together with other advanced indices in the manner of legacy
   "fancy indexing" is generally not helpful or well defined.
   A user who wishes the "``nonzero``" plus broadcast behaviour can thus
   be expected to do this manually.
   Using this rule, a single boolean index can index into multiple
   dimensions at once.
5. An ``arr.lindex`` or ``arr.findex`` should likely be implemented to
   legacy fancy indexing indefinetly. This also gives a simple way to
   update fancy indexing code making deprecations to vanilla indexing
6. Vanilla indexing ``arr[...]`` could return an error for ambiguous
   For the beginning, this probably means cases where ``arr[ind]`` and
   ``arr.oindex[ind]`` return different results gives deprecation
   However, the exact rules for this (especially the final behaviour)
are not
   quite clear in cases such as ``arr[0, :, index_arr]``.

All other rules for indexing are identical. 

Open Questions

1. Especially for the new indexing attributes ``oindex`` and ``vindex``,
   a case could be made to not implicitly add an ``Ellipsis`` index if
   This helps finding bugs since a too high dimensional array can be
   (I am in favor for this, but doubt we should think about this for

2. The names ``oindex`` and ``vindex`` are just suggestions at the time
   writing this, another name NumPy has used for something like
   is ``np.ix_``. See also below.

3. It would be possible to limit the use of boolean indices in
   assuming that they are rare and to some degree special.
   (This would make implementation simpler, but I do not see a big

4. ``oindex`` and ``vindex`` could always return copies, even when no
   operation occurs. One argument for using the same rules is that this
   ``oindex`` can be used as a general index replacement.
   (There is likely no big reason for this, however, there is one
   ``arr.vindex[array_scalar, ...]`` can occur, where ``arr_scalar``
   should be a 0-D array. Copying always "fixes" the possible

5. The final state to morph indexing in is not fixed in this PEP. It is
   example possible that `arr[index]`` will be equivalent to
   at some point in the future. Since such a change will take years, it
   seems unnecessary to make specific decisions now.

6. Proposed changes to vanilla indexing could be postponed indefinetly
   not taken in order to not break or force fixing of existing code

7. Possible the ``vindex`` combination with boolean indexing could be
   rethought or not allowed at all for simplicity.

Necessary changes to NumPy

Implement ``arr.oindex`` and ``arr.vindex`` objects to allow these
operations and create warnings (and eventually deprecate) ambiguous
indexing operations on arrays.

Alternative Names

Possible names suggested (more suggestions will be added).

==============  ======== =======
**Orthogonal**  oindex   oix
**Vectorized**  vindex   fix
**Legacy**      l/findex 
==============  ======== =======


Since the various kinds of indexing is hard to grasp in many cases,
examples hopefully give some more insights. Note that they are all in
of shape. All original dimensions start with 5, advanced indexing
inserts less long dimensions. (Note that ``...`` or ``Ellipsis`` mostly
inserts as many slices as needed to index the full array). These
may be hard to grasp without working knowledge of advanced indexing as
NumPy 1.9.

Example array::

    >>> arr = np.ones((5, 6, 7, 8))

Legacy fancy indexing

Single index is transposed (this is the same for all indexing types)::

    >>> arr[[0], ...].shape
    (1, 6, 7, 8)
    >>> arr[:, [0], ...].shape
    (5, 1, 7, 8)

Multiple indices are transposed *if* consecutive::

    >>> arr[:, [0], [0], :].shape  # future error
    (5, 1, 7)
    >>> arr[:, [0], :, [0]].shape  # future error
    (1, 5, 6)

It is important to note that a scalar *is* integer array index in this
(and gets broadcasted with the other advanced index)::

    >>> arr[:, [0], 0, :].shape  # future error (scalar is "fancy")
    (5, 1, 7)
    >>> arr[:, [0], :, 0].shape  # future error (scalar is "fancy")
    (1, 5, 6)

Single boolean index can act on multiple dimensions (especially the
array). It has to match (as of 1.10. a deprecation warning) the
The boolean index is otherwise identical to (multiple consecutive)
array indices::

    >>> # Create boolean index with one True value for the last two
    >>> bindx = np.zeros((7, 8), dtype=np.bool_)
    >>> bindx[[0, 0]] = True
    >>> arr[:, 0, bindx].shape
    (5, 1)
    >>> arr[0, :, bindx].shape
    (1, 6)

The combination with anything that is not a scalar is confusing, e.g.::

    >>> arr[[0], :, bindx].shape  # bindx result broadcasts with [0]
    (1, 6)
    >>> arr[:, [0, 1], bindx]  # IndexError

Outer indexing

Multiple indices are "orthogonal" and their result axes are inserted 
at the same place (they are not broadcasted)::

    >>> arr.oindex[:, [0], [0, 1], :].shape
    (5, 1, 2, 8)
    >>> arr.oindex[:, [0], :, [0, 1]].shape
    (5, 1, 7, 2)
    >>> arr.oindex[:, [0], 0, :].shape
    (5, 1, 8)
    >>> arr.oindex[:, [0], :, 0].shape
    (5, 1, 7)

Boolean indices results are always inserted where the index is::

    >>> # Create boolean index with one True value for the last two
    >>> bindx = np.zeros((7, 8), dtype=np.bool_)
    >>> bindx[[0, 0]] = True
    >>> arr.oindex[:, 0, bindx].shape
    (5, 1)
    >>> arr.oindex[0, :, bindx].shape
    (6, 1)

Nothing changed in the presence of other advanced indices since::

    >>> arr.oindex[[0], :, bindx].shape
    (1, 6, 1)
    >>> arr.oindex[:, [0, 1], bindx]
    (5, 2, 1)

Vectorized/inner indexing

Multiple indices are broadcasted and iterated as one like fancy
but the new axes area always inserted at the front::

    >>> arr.vindex[:, [0], [0, 1], :].shape
    (2, 5, 8)
    >>> arr.vindex[:, [0], :, [0, 1]].shape
    (2, 5, 7)
    >>> arr.vindex[:, [0], 0, :].shape
    (1, 5, 8)
    >>> arr.vindex[:, [0], :, 0].shape
    (1, 5, 7)

Boolean indices results are always inserted where the index is, exactly
as in ``oindex`` given how specific they are to the axes they operate

    >>> # Create boolean index with one True value for the last two
    >>> bindx = np.zeros((7, 8), dtype=np.bool_)
    >>> bindx[[0, 0]] = True
    >>> arr.vindex[:, 0, bindx].shape
    (5, 1)
    >>> arr.vindex[0, :, bindx].shape
    (6, 1)

But other advanced indices are again transposed to the front::

    >>> arr.vindex[[0], :, bindx].shape
    (1, 6, 1)
    >>> arr.vindex[:, [0, 1], bindx]
    (2, 5, 1)

Related Questions

There exist a further indexing or indexing like method. That is the
inverse of a command such as ``np.argmin(arr, axis=axis)``, to pick
the specific elements *along* an axis given an array of (at least
typically) the same size.

Doing such a thing with the indexing notation is not quite straight
since the axis on which to pick elements has to be supplied. One
solution would be to create a function (calling it pick here for

     np.pick(arr, index_arr, axis=axis)

where ``index_arr`` has to be the same shape as ``arr`` except along
One could imagine that this can be useful together with other indexing
but such a function may be sufficient and extra information needed seems
to pass using a function convention. Another option would be to allow an
such as ``compress_axes=None`` (just to have some name) which maps the
axes from
the index array to the new array with ``None`` signaling a new axis.
Also keepdims could be added as a simple default. (Note that the use of
axis is not
compatible to ``np.take`` for an ``index_arr`` which is not zero or one

Another solution is to provide functions or features to the
to map this to the equivalent ``vindex`` indexing operation.

Motivational Example

Imagine having a data acquisition software storing ``D`` channels and
``N`` datapoints along the time. She stores this into an ``(N, D)``
array. During data analysis, we needs to fetch a pool of channels, for
to calculate a mean over them.

This data can be faked using::

    >>> arr = np.random.random((100, 10))

Now one may remember indexing with an integer array and find the correct

    >>> group = arr[:, [2, 5]]
    >>> mean_value = arr.mean()

However, assume that there were some specific time points (first
of the data) that need to be specially considered. These time points are
already known and given by::

    >>> interesting_times = np.array([1, 5, 8, 10], dtype=np.intp)

Now to fetch them, we may try to modify the previous code::

    >>> group_at_it = arr[interesting_times, [2, 5]]
    IndexError: Ambiguous index, use `.oindex` or `.vindex`

An error such as this will point to read up the indexing documentation.
This should make it clear, that ``oindex`` behaves more like slicing.
So, out of the different methods it is the obvious choice
(for now, this is a shape mismatch, but that could possibly also mention

    >>> group_at_it = arr.oindex[interesting_times, [2, 5]]

Now of course one could also have used ``vindex``, but it is much less
obvious how to achieve the right thing!::

    >>> reshaped_times = interesting_times[:, np.newaxis]
    >>> group_at_it = arr.vindex[reshaped_times, [2, 5]]

One may find, that for example our data is corrupt in some places.
So, we need to replace these values by zero (or anything else) for these
times. The first column may for example give the necessary information,
so that changing the values becomes easy remembering boolean indexing::

    >>> bad_data = arr[0] > 0.5
    >>> arr[bad_data, :] = 0

Again, however, the columns may need to be handled more individually
(but in
groups), and the ``oindex`` attribute works well::

    >>> arr.oindex[bad_data, [2, 5]] = 0

Note that it would be very hard to do this using legacy fancy indexing.
The only way would be to create an integer array first::

    >>> bad_data_indx = np.nonzero(bad_data)[0]
    >>> bad_data_indx_reshaped = bad_data_indx[:, np.newaxis]
    >>> arr[bad_data_indx_reshaped, [2, 5]]

In any case we can use only ``oindex`` to do all of this without getting
into any trouble or confused by the whole complexity of advanced

But, some new features are added to the data acquisition. Different
have to be used depending on the times. Let us assume we already have
created an array of indices::

    >>> correct_sensors = np.random.randint(10, size=(100, 2))

Which lists for each time the two correct sensors in an ``(N, 2)``

A first try to achieve this may be ``arr[:, correct_sensors]`` and this
not work. It should be clear quickly that slicing cannot achieve the
thing. But hopefully users will remember that there is ``vindex`` as a
powerful and flexible approach to advanced indexing.
One may, if trying ``vindex`` randomly, be confused about::

    >>> new_arr = arr.vindex[:, correct_sensors]

which is neither the same, nor the correct result (see transposing
This is because slicing works still the same in ``vindex``. However,
the documentation and examples, one can hopefully quickly find the

    >>> rows = np.arange(len(arr))
    >>> rows = rows[:, np.newaxis]  # make shape fit with
    >>> new_arr = arr.vindex[rows, correct_sensors]
At this point we have left the straight forward world of ``oindex`` but
do random picking of any element from the array. Note that in the last
a method such as mentioned in the ``Related Questions`` section could be
straight forward. But this approach is even more flexible, since
does not have to be a simple ``arange``, but could be

    >>> correct_sensors_at_it = correct_sensors[interesting_times, :]
    >>> interesting_times_reshaped = interesting_times[:, np.newaxis]
    >>> new_arr_it = arr[interesting_times_reshaped,

Truly complex situation would arise now if you would for example pool
experiments into an array shaped ``(L, N, D)``. But for ``oindex`` this
not result into surprises. ``vindex``, being more powerful, will quite
certainly create some confusion in this case but also cover pretty much

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 819 bytes
Desc: This is a digitally signed message part
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20151111/d7119827/attachment.sig>

More information about the NumPy-Discussion mailing list