<div dir="ltr"><div>Sebastian and I have revised a Numpy Enhancement Proposal that he started three years ago for overhauling NumPy's advanced indexing. We'd now like to present it for official consideration.</div><div><br></div><div>Minor inline comments (e.g., typos) can be added to the latest pull request (<a href="https://github.com/numpy/numpy/pull/11414/files">https://github.com/numpy/numpy/pull/11414/files</a>), but otherwise let's keep discussion on the mailing list. The NumPy website should update shortly with a rendered version (<a href="http://www.numpy.org/neps/nep-0021-advanced-indexing.html">http://www.numpy.org/neps/nep-0021-advanced-indexing.html</a>), but until then please see the full text below.</div><div><br></div><div>Cheers,</div><div>Stephan</div><div><br></div><div>=========================================</div><div>Simplified and explicit advanced indexing</div><div>=========================================</div><div><br></div><div>:Author: Sebastian Berg</div><div>:Author: Stephan Hoyer <<a href="mailto:shoyer@google.com">shoyer@google.com</a>></div><div>:Status: Draft</div><div>:Type: Standards Track</div><div>:Created: 2015-08-27</div><div><br></div><div><br></div><div>Abstract</div><div>--------</div><div><br></div><div>NumPy's "advanced" indexing support for indexing arrays with other arrays is</div><div>one of its most powerful and popular features. Unfortunately, the existing</div><div>rules for advanced indexing with multiple array indices are typically confusing</div><div>to both new, and in many cases even old, users of NumPy. Here we propose an</div><div>overhaul and simplification of advanced indexing, including two new "indexer"</div><div>attributes ``oindex`` and ``vindex`` to facilitate explicit indexing.</div><div><br></div><div>Background</div><div>----------</div><div><br></div><div>Existing indexing operations</div><div>~~~~~~~~~~~~~~~~~~~~~~~~~~~~</div><div><br></div><div>NumPy arrays currently support a flexible range of indexing operations:</div><div><br></div><div>- "Basic" indexing involving only slices, integers, ``np.newaxis`` and ellipsis</div><div>  (``...``), e.g., ``x[0, :3, np.newaxis]`` for selecting the first element</div><div>  from the 0th axis, the first three elements from the 1st axis and inserting a</div><div>  new axis of size 1 at the end. Basic indexing always return a view of the</div><div>  indexed array's data.</div><div>- "Advanced" indexing, also called "fancy" indexing, includes all cases where</div><div>  arrays are indexed by other arrays. Advanced indexing always makes a copy:</div><div><br></div><div>  - "Boolean" indexing by boolean arrays, e.g., ``x[x > 0]`` for</div><div>    selecting positive elements.</div><div>  - "Vectorized" indexing by one or more integer arrays, e.g., ``x[[0, 1]]``</div><div>    for selecting the first two elements along the first axis. With multiple</div><div>    arrays, vectorized indexing uses broadcasting rules to combine indices along</div><div>    multiple dimensions. This allows for producing a result of arbitrary shape</div><div>    with arbitrary elements from the original arrays.</div><div>  - "Mixed" indexing involving any combinations of the other advancing types.</div><div>    This is no more powerful than vectorized indexing, but is sometimes more</div><div>    convenient.</div><div><br></div><div>For clarity, we will refer to these existing rules as "legacy indexing".</div><div>This is only a high-level summary; for more details, see NumPy's documentation</div><div>and and `Examples` below.</div><div><br></div><div>Outer indexing</div><div>~~~~~~~~~~~~~~</div><div><br></div><div>One broadly useful class of indexing operations is not supported:</div><div><br></div><div>- "Outer" or orthogonal indexing treats one-dimensional arrays equivalently to</div><div>  slices for determining output shapes. The rule for outer indexing is that the</div><div>  result should be equivalent to independently indexing along each dimension</div><div>  with integer or boolean arrays as if both the indexed and indexing arrays</div><div>  were one-dimensional. This form of indexing is familiar to many users of other</div><div>  programming languages such as MATLAB, Fortran and R.</div><div><br></div><div>The reason why NumPy omits support for outer indexing is that the rules for</div><div>outer and vectorized conflict. Consider indexing a 2D array by two 1D integer</div><div>arrays, e.g., ``x[[0, 1], [0, 1]]``:</div><div><br></div><div>- Outer indexing is equivalent to combining multiple integer indices with</div><div>  ``itertools.product()``. The result in this case is another 2D array with</div><div>  all combinations of indexed elements, e.g.,</div><div>  ``np.array([[x[0, 0], x[0, 1]], [x[1, 0], x[1, 1]]])``</div><div>- Vectorized indexing is equivalent to combining multiple integer indices with</div><div>  ``zip()``. The result in this case is a 1D array containing the diagonal</div><div>  elements, e.g., ``np.array([x[0, 0], x[1, 1]])``.</div><div><br></div><div>This difference is a frequent stumbling block for new NumPy users. The outer</div><div>indexing model is easier to understand, and is a natural generalization of</div><div>slicing rules. But NumPy instead chose to support vectorized indexing, because</div><div>it is strictly more powerful.</div><div><br></div><div>It is always possible to emulate outer indexing by vectorized indexing with</div><div>the right indices. To make this easier, NumPy includes utility objects and</div><div>functions such as ``np.ogrid`` and ``np.ix_``, e.g.,</div><div>``x[np.ix_([0, 1], [0, 1])]``. However, there are no utilities for emulating</div><div>fully general/mixed outer indexing, which could unambiguously allow for slices,</div><div>integers, and 1D boolean and integer arrays.</div><div><br></div><div>Mixed indexing</div><div>~~~~~~~~~~~~~~</div><div><br></div><div>NumPy's existing rules for combining multiple types of indexing in the same</div><div>operation are quite complex, involving a number of edge cases.</div><div><br></div><div>One reason why mixed indexing is particularly confusing is that at first glance</div><div>the result works deceptively like outer indexing. Returning to our example of a</div><div>2D array, both ``x[:2, [0, 1]]`` and ``x[[0, 1], :2]`` return 2D arrays with</div><div>axes in the same order as the original array.</div><div><br></div><div>However, as soon as two or more non-slice objects (including integers) are</div><div>introduced, vectorized indexing rules apply. The axes introduced by the array</div><div>indices are at the front, unless all array indices are consecutive, in which</div><div>case NumPy deduces where the user "expects" them to be. Consider indexing a 3D</div><div>array ``arr`` with shape ``(X, Y, Z)``:</div><div><br></div><div>1. ``arr[:, [0, 1], 0]`` has shape ``(X, 2)``.</div><div>2. ``arr[[0, 1], 0, :]`` has shape ``(2, Z)``.</div><div>3. ``arr[0, :, [0, 1]]`` has shape ``(2, Y)``, not ``(Y, 2)``!</div><div><br></div><div>These first two cases are intuitive and consistent with outer indexing, but</div><div>this last case is quite surprising, even to many higly experienced NumPy users.</div><div><br></div><div>Mixed cases involving multiple array indices are also surprising, and only</div><div>less problematic because the current behavior is so useless that it is rarely</div><div>encountered in practice. When a boolean array index is mixed with another boolean or</div><div>integer array, boolean array is converted to integer array indices (equivalent</div><div>to ``np.nonzero()``) and then broadcast. For example, indexing a 2D array of</div><div>size ``(2, 2)`` like ``x[[True, False], [True, False]]`` produces a 1D vector</div><div>with shape ``(1,)``, not a 2D sub-matrix with shape ``(1, 1)``.</div><div><br></div><div>Mixed indexing seems so tricky that it is tempting to say that it never should</div><div>be used. However, it is not easy to avoid, because NumPy implicitly adds full</div><div>slices if there are fewer indices than the full dimensionality of the indexed</div><div>array. This means that indexing a 2D array like `x[[0, 1]]`` is equivalent to</div><div>``x[[0, 1], :]``. These cases are not surprising, but they constrain the</div><div>behavior of mixed indexing.</div><div><br></div><div>Indexing in other Python array libraries</div><div>~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~</div><div><br></div><div>Indexing is a useful and widely recognized mechanism for accessing</div><div>multi-dimensional array data, so it is no surprise that many other libraries in</div><div>the scientific Python ecosystem also support array indexing.</div><div><br></div><div>Unfortunately, the full complexity of NumPy's indexing rules mean that it is</div><div>both challenging and undesirable for other libraries to copy its behavior in all</div><div>of its nuance. The only full implementation of NumPy-style indexing is NumPy</div><div>itself. This includes projects like dask.array and h5py, which support *most*</div><div>types of array indexing in some form, and otherwise attempt to copy NumPy's API</div><div>exactly.</div><div><br></div><div>Vectorized indexing in particular can be challenging to implement with array</div><div>storage backends not based on NumPy. In contrast, indexing by 1D arrays along</div><div>at least one dimension in the style of outer indexing is much more acheivable.</div><div>This has led many libraries (including dask and h5py) to attempt to define a</div><div>safe subset of NumPy-style indexing that is equivalent to outer indexing, e.g.,</div><div>by only allowing indexing with an array along at most one dimension. However,</div><div>this is quite challenging to do correctly in a general enough way to be useful.</div><div>For example, the current versions of dask and h5py both handle mixed indexing</div><div>in case 3 above inconsistently with NumPy. This is quite likely to lead to</div><div>bugs.</div><div><br></div><div>These inconsistencies, in addition to the broader challenge of implementing</div><div>every type of indexing logic, make it challenging to write high-level array</div><div>libraries like xarray or dask.array that can interchangeably index many types of</div><div>array storage. In contrast, explicit APIs for outer and vectorized indexing in</div><div>NumPy would provide a model that external libraries could reliably emulate, even</div><div>if they don't support every type of indexing.</div><div><br></div><div>High level changes</div><div>------------------</div><div><br></div><div>Inspired by multiple "indexer" attributes for controlling different types</div><div>of indexing behavior in pandas, we propose to:</div><div><br></div><div>1. Introduce ``arr.oindex[indices]`` which allows array indices, but</div><div>   uses outer indexing logic.</div><div>2. Introduce ``arr.vindex[indices]`` which use the current</div><div>   "vectorized"/broadcasted logic but with two differences from</div><div>   legacy indexing:</div><div>       </div><div>   * Boolean indices are not supported. All indices must be integers,</div><div>     integer arrays or slices.</div><div>   * The integer index result dimensions are always the first axes</div><div>     of the result array. No transpose is done, even for a single</div><div>     integer array index.</div><div><br></div><div>3. Plain indexing on arrays will start to give warnings and eventually</div><div>   errors in cases where one of the explicit indexers should be preferred:</div><div><br></div><div>   * First, in all cases where legacy and outer indexing would give</div><div>     different results.</div><div>   * Later, potentially in all cases involving an integer array.</div><div><br></div><div>These constraints are sufficient for making indexing generally consistent</div><div>with expectations and providing a less surprising learning curve with</div><div>``oindex``.</div><div><br></div><div>Note that all things mentioned here apply both for assignment as well as</div><div>subscription.</div><div><br></div><div>Understanding these details is *not* easy. The `Examples` section in the</div><div>discussion gives code examples.</div><div>And the hopefully easier `Motivational Example` provides some</div><div>motivational use-cases for the general ideas and is likely a good start for</div><div>anyone not intimately familiar with advanced indexing.</div><div><br></div><div><br></div><div>Detailed Description</div><div>--------------------</div><div><br></div><div>Proposed rules</div><div>~~~~~~~~~~~~~~</div><div><br></div><div>From the three problems noted above some expectations for NumPy can</div><div>be deduced:</div><div><br></div><div>1. There should be a prominent outer/orthogonal indexing method such as</div><div>   ``arr.oindex[indices]``.</div><div><br></div><div>2. Considering how confusing vectorized/fancy indexing can be, it should</div><div>   be possible to be made more explicitly (e.g. ``arr.vindex[indices]``).</div><div><br></div><div>3. A new ``arr.vindex[indices]`` method, would not be tied to the</div><div>   confusing transpose rules of fancy indexing, which is for example</div><div>   needed for the simple case of a single advanced index. Thus,</div><div>   no transposing should be done. The axes created by the integer array</div><div>   indices are always inserted at the front, even for a single index.</div><div><br></div><div>4. Boolean indexing is conceptionally outer indexing. Broadcasting</div><div>   together with other advanced indices in the manner of legacy</div><div>   indexing is generally not helpful or well defined.</div><div>   A user who wishes the "``nonzero``" plus broadcast behaviour can thus</div><div>   be expected to do this manually. Thus, ``vindex`` does not need to</div><div>   support boolean index arrays.</div><div><br></div><div>5. An ``arr.legacy_index`` attribute should be implemented to support</div><div>   legacy indexing. This gives a simple way to update existing codebases</div><div>   using legacy indexing, which will make the deprecation of plain indexing</div><div>   behavior easier. The longer name ``legacy_index`` is intentionally chosen</div><div>   to be explicit and discourage its use in new code.</div><div><br></div><div>6. Plain indexing ``arr[...]`` should return an error for ambiguous cases.</div><div>   For the beginning, this probably means cases where ``arr[ind]`` and</div><div>   ``arr.oindex[ind]`` return different results give deprecation warnings.</div><div>   This includes every use of vectorized indexing with multiple integer arrays.</div><div>   Due to the transposing behaviour, this means that``arr[0, :, index_arr]``</div><div>   will be deprecated, but ``arr[:, 0, index_arr]`` will not for the time being.</div><div><br></div><div>7. To ensure that existing subclasses of `ndarray` that override indexing</div><div>   do not inadvertently revert to default behavior for indexing attributes,</div><div>   these attribute should have explicit checks that disable them if</div><div>   ``__getitem__`` or ``__setitem__`` has been overriden.</div><div><br></div><div>Unlike plain indexing, the new indexing attributes are explicitly aimed</div><div>at higher dimensional indexing, several additional changes should be implemented:</div><div><br></div><div>* The indexing attributes will enforce exact dimension and indexing match.</div><div>  This means that no implicit ellipsis (``...``) will be added. Unless</div><div>  an ellipsis is present the indexing expression will thus only work for</div><div>  an array with a specific number of dimensions.</div><div>  This makes the expression more explicit and safeguards against wrong</div><div>  dimensionality of arrays.</div><div>  There should be no implications for "duck typing" compatibility with</div><div>  builtin Python sequences, because Python sequences only support a limited</div><div>  form of "basic indexing" with integers and slices.</div><div><br></div><div>* The current plain indexing allows for the use of non-tuples for</div><div>  multi-dimensional indexing such as ``arr[[slice(None), 2]]``.</div><div>  This creates some inconsistencies and thus the indexing attributes</div><div>  should only allow plain python tuples for this purpose.</div><div>  (Whether or not this should be the case for plain indexing is a</div><div>  different issue.)</div><div><br></div><div>* The new attributes should not use getitem to implement setitem,</div><div>  since it is a cludge and not useful for vectorized</div><div>  indexing. (not implemented yet)</div><div><br></div><div><br></div><div>Open Questions</div><div>~~~~~~~~~~~~~~</div><div><br></div><div>* The names ``oindex``, ``vindex`` and ``legacy_index`` are just suggestions at</div><div>  the time of writing this, another name NumPy has used for something like</div><div>  ``oindex`` is ``np.ix_``. See also below.</div><div><br></div><div>* ``oindex`` and ``vindex`` could always return copies, even when no array</div><div>  operation occurs. One argument for allowing a view return is that this way</div><div>  ``oindex`` can be used as a general index replacement.</div><div>  However, there is one argument for returning copies. It is possible for</div><div>  ``arr.vindex[array_scalar, ...]``, where ``array_scalar`` should be</div><div>  a 0-D array but is not, since 0-D arrays tend to be converted.</div><div>  Copying always "fixes" this possible inconsistency.</div><div><br></div><div>* The final state to morph plain indexing in is not fixed in this PEP.</div><div>  It is for example possible that `arr[index]`` will be equivalent to</div><div>  ``arr.oindex`` at some point in the future.</div><div>  Since such a change will take years, it seems unnecessary to make</div><div>  specific decisions at this time.</div><div><br></div><div>* The proposed changes to plain indexing could be postponed indefinitely or</div><div>  not taken in order to not break or force major fixes to existing code bases.</div><div><br></div><div><br></div><div>Alternative Names</div><div>~~~~~~~~~~~~~~~~~</div><div><br></div><div>Possible names suggested (more suggestions will be added).</div><div><br></div><div>==============  ============ ========</div><div>**Orthogonal**  oindex       oix</div><div>**Vectorized**  vindex       vix</div><div>**Legacy**      legacy_index l/findex</div><div>==============  ============ ========</div><div><br></div><div><br></div><div>Subclasses</div><div>~~~~~~~~~~</div><div><br></div><div>Subclasses are a bit problematic in the light of these changes. There are</div><div>some possible solutions for this. For most subclasses (those which do not</div><div>provide ``__getitem__`` or ``__setitem__``) the special attributes should</div><div>just work. Subclasses that *do* provide it must be updated accordingly</div><div>and should preferably not subclass working versions of these attributes.</div><div><br></div><div>All subclasses will inherit the attributes, however, the implementation</div><div>of ``__getitem__`` on these attributes should test</div><div>``subclass.__getitem__ is ndarray.__getitem__``. If not, the</div><div>subclass has special handling for indexing and ``NotImplementedError``</div><div>should be raised, requiring that the indexing attributes is also explicitly</div><div>overwritten. Likewise, implementations of ``__setitem__`` should check to see</div><div>if ``__setitem__`` is overriden.</div><div><br></div><div>A further question is how to facilitate implementing the special attributes.</div><div>Also there is the weird functionality where ``__setitem__`` calls</div><div>``__getitem__`` for non-advanced indices. It might be good to avoid it for</div><div>the new attributes, but on the other hand, that may make it even more</div><div>confusing.</div><div><br></div><div>To facilitate implementations we could provide functions similar to</div><div>``operator.itemgetter`` and ``operator.setitem`` for the attributes.</div><div>Possibly a mixin could be provided to help implementation. These improvements</div><div>are not essential to the initial implementation, so they are saved for</div><div>future work.</div><div><br></div><div>Implementation</div><div>--------------</div><div><br></div><div>Implementation would start with writing special indexing objects available</div><div>through ``arr.oindex``, ``arr.vindex``, and ``arr.legacy_index`` to allow these</div><div>indexing operations. Also, we would need to start to deprecate those plain index</div><div>operations which are not ambiguous.</div><div>Furthermore, the NumPy code base will need to use the new attributes and</div><div>tests will have to be adapted.</div><div><br></div><div><br></div><div>Backward compatibility</div><div>----------------------</div><div><br></div><div>As a new feature, no backward compatibility issues with the new ``vindex``</div><div>and ``oindex`` attributes would arise. To facilitate backwards compatibility</div><div>as much as possible, we expect a long deprecation cycle for legacy indexing</div><div>behavior and propose the new ``legacy_index`` attribute.</div><div>Some forward compatibility issues with subclasses that do not specifically</div><div>implement the new methods may arise.</div><div><br></div><div><br></div><div>Alternatives</div><div>------------</div><div><br></div><div>NumPy may not choose to offer these different type of indexing methods, or</div><div>choose to only offer them through specific functions instead of the proposed</div><div>notation above.</div><div><br></div><div>We don't think that new functions are a good alternative, because indexing</div><div>notation ``[]`` offer some syntactic advantages in Python (i.e., direct</div><div>creation of slice objects) compared to functions.</div><div><br></div><div>A more reasonable alternative would be write new wrapper objects for alternative</div><div>indexing with functions rather than methods (e.g., ``np.oindex(arr)[indices]``</div><div>instead of ``arr.oindex[indices]``). Functionally, this would be equivalent,</div><div>but indexing is such a common operation that we think it is important to</div><div>minimize syntax and worth implementing it directly on `ndarray` objects</div><div>themselves. Indexing attributes also define a clear interface that is easier</div><div>for alternative array implementations to copy, nonwithstanding ongoing</div><div>efforts to make it easier to override NumPy functions [2]_.</div><div><br></div><div>Discussion</div><div>----------</div><div><br></div><div>The original discussion about vectorized vs outer/orthogonal indexing arose</div><div>on the NumPy mailing list:</div><div><br></div><div> * <a href="https://mail.python.org/pipermail/numpy-discussion/2015-April/072550.html">https://mail.python.org/pipermail/numpy-discussion/2015-April/072550.html</a></div><div><br></div><div>Some discussion can be found on the original pull request for this NEP:</div><div><br></div><div> * <a href="https://github.com/numpy/numpy/pull/6256">https://github.com/numpy/numpy/pull/6256</a></div><div><br></div><div>Python implementations of the indexing operations can be found at:</div><div><br></div><div> * <a href="https://github.com/numpy/numpy/pull/5749">https://github.com/numpy/numpy/pull/5749</a></div><div> * <a href="https://gist.github.com/shoyer/c700193625347eb68fee4d1f0dc8c0c8">https://gist.github.com/shoyer/c700193625347eb68fee4d1f0dc8c0c8</a></div><div><br></div><div><br></div><div>Examples</div><div>~~~~~~~~</div><div><br></div><div>Since the various kinds of indexing is hard to grasp in many cases, these</div><div>examples hopefully give some more insights. Note that they are all in terms</div><div>of shape.</div><div>In the examples, all original dimensions have 5 or more elements,</div><div>advanced indexing inserts smaller dimensions.</div><div>These examples may be hard to grasp without working knowledge of advanced</div><div>indexing as of NumPy 1.9.</div><div><br></div><div>Example array::</div><div><br></div><div>    >>> arr = np.ones((5, 6, 7, 8))</div><div><br></div><div><br></div><div>Legacy fancy indexing</div><div>---------------------</div><div><br></div><div>Note that the same result can be achieved with ``arr.legacy_index``, but the</div><div>"future error" will still work in this case.</div><div><br></div><div>Single index is transposed (this is the same for all indexing types)::</div><div><br></div><div>    >>> arr[[0], ...].shape</div><div>    (1, 6, 7, 8)</div><div>    >>> arr[:, [0], ...].shape</div><div>    (5, 1, 7, 8)</div><div><br></div><div><br></div><div>Multiple indices are transposed *if* consecutive::</div><div><br></div><div>    >>> arr[:, [0], [0], :].shape  # future error</div><div>    (5, 1, 8)</div><div>    >>> arr[:, [0], :, [0]].shape  # future error</div><div>    (1, 5, 7)</div><div><br></div><div><br></div><div>It is important to note that a scalar *is* integer array index in this sense</div><div>(and gets broadcasted with the other advanced index)::</div><div><br></div><div>    >>> arr[:, [0], 0, :].shape</div><div>    (5, 1, 8)</div><div>    >>> arr[:, [0], :, 0].shape  # future error (scalar is "fancy")</div><div>    (1, 5, 7)</div><div><br></div><div><br></div><div>Single boolean index can act on multiple dimensions (especially the whole</div><div>array). It has to match (as of 1.10. a deprecation warning) the dimensions.</div><div>The boolean index is otherwise identical to (multiple consecutive) integer</div><div>array indices::</div><div><br></div><div>    >>> # Create boolean index with one True value for the last two dimensions:</div><div>    >>> bindx = np.zeros((7, 8), dtype=np.bool_)</div><div>    >>> bindx[0, 0] = True</div><div>    >>> arr[:, 0, bindx].shape</div><div>    (5, 1)</div><div>    >>> arr[0, :, bindx].shape</div><div>    (1, 6)</div><div><br></div><div><br></div><div>The combination with anything that is not a scalar is confusing, e.g.::</div><div><br></div><div>    >>> arr[[0], :, bindx].shape  # bindx result broadcasts with [0]</div><div>    (1, 6)</div><div>    >>> arr[:, [0, 1], bindx].shape  # IndexError</div><div><br></div><div><br></div><div>Outer indexing</div><div>--------------</div><div><br></div><div>Multiple indices are "orthogonal" and their result axes are inserted </div><div>at the same place (they are not broadcasted)::</div><div><br></div><div>    >>> arr.oindex[:, [0], [0, 1], :].shape</div><div>    (5, 1, 2, 8)</div><div>    >>> arr.oindex[:, [0], :, [0, 1]].shape</div><div>    (5, 1, 7, 2)</div><div>    >>> arr.oindex[:, [0], 0, :].shape</div><div>    (5, 1, 8)</div><div>    >>> arr.oindex[:, [0], :, 0].shape</div><div>    (5, 1, 7)</div><div><br></div><div><br></div><div>Boolean indices results are always inserted where the index is::</div><div><br></div><div>    >>> # Create boolean index with one True value for the last two dimensions:</div><div>    >>> bindx = np.zeros((7, 8), dtype=np.bool_)</div><div>    >>> bindx[0, 0] = True</div><div>    >>> arr.oindex[:, 0, bindx].shape</div><div>    (5, 1)</div><div>    >>> arr.oindex[0, :, bindx].shape</div><div>    (6, 1)</div><div><br></div><div><br></div><div>Nothing changed in the presence of other advanced indices since::</div><div><br></div><div>    >>> arr.oindex[[0], :, bindx].shape</div><div>    (1, 6, 1)</div><div>    >>> arr.oindex[:, [0, 1], bindx].shape</div><div>    (5, 2, 1)</div><div><br></div><div><br></div><div>Vectorized/inner indexing</div><div>-------------------------</div><div><br></div><div>Multiple indices are broadcasted and iterated as one like fancy indexing,</div><div>but the new axes area always inserted at the front::</div><div><br></div><div>    >>> arr.vindex[:, [0], [0, 1], :].shape</div><div>    (2, 5, 8)</div><div>    >>> arr.vindex[:, [0], :, [0, 1]].shape</div><div>    (2, 5, 7)</div><div>    >>> arr.vindex[:, [0], 0, :].shape</div><div>    (1, 5, 8)</div><div>    >>> arr.vindex[:, [0], :, 0].shape</div><div>    (1, 5, 7)</div><div><br></div><div><br></div><div>Boolean indices results are always inserted where the index is, exactly</div><div>as in ``oindex`` given how specific they are to the axes they operate on::</div><div><br></div><div>    >>> # Create boolean index with one True value for the last two dimensions:</div><div>    >>> bindx = np.zeros((7, 8), dtype=np.bool_)</div><div>    >>> bindx[0, 0] = True</div><div>    >>> arr.vindex[:, 0, bindx].shape</div><div>    (5, 1)</div><div>    >>> arr.vindex[0, :, bindx].shape</div><div>    (6, 1)</div><div><br></div><div><br></div><div>But other advanced indices are again transposed to the front::</div><div><br></div><div>    >>> arr.vindex[[0], :, bindx].shape</div><div>    (1, 6, 1)</div><div>    >>> arr.vindex[:, [0, 1], bindx].shape</div><div>    (2, 5, 1)</div><div><br></div><div><br></div><div>Motivational Example</div><div>~~~~~~~~~~~~~~~~~~~~</div><div><br></div><div>Imagine having a data acquisition software storing ``D`` channels and</div><div>``N`` datapoints along the time. She stores this into an ``(N, D)`` shaped</div><div>array. During data analysis, we needs to fetch a pool of channels, for example</div><div>to calculate a mean over them.</div><div><br></div><div>This data can be faked using::</div><div><br></div><div>    >>> arr = np.random.random((100, 10))</div><div><br></div><div>Now one may remember indexing with an integer array and find the correct code::</div><div><br></div><div>    >>> group = arr[:, [2, 5]]</div><div>    >>> mean_value = arr.mean()</div><div><br></div><div>However, assume that there were some specific time points (first dimension</div><div>of the data) that need to be specially considered. These time points are</div><div>already known and given by::</div><div><br></div><div>    >>> interesting_times = np.array([1, 5, 8, 10], dtype=np.intp)</div><div><br></div><div>Now to fetch them, we may try to modify the previous code::</div><div><br></div><div>    >>> group_at_it = arr[interesting_times, [2, 5]]</div><div>    IndexError: Ambiguous index, use `.oindex` or `.vindex`</div><div><br></div><div>An error such as this will point to read up the indexing documentation.</div><div>This should make it clear, that ``oindex`` behaves more like slicing.</div><div>So, out of the different methods it is the obvious choice</div><div>(for now, this is a shape mismatch, but that could possibly also mention</div><div>``oindex``)::</div><div><br></div><div>    >>> group_at_it = arr.oindex[interesting_times, [2, 5]]</div><div><br></div><div>Now of course one could also have used ``vindex``, but it is much less</div><div>obvious how to achieve the right thing!::</div><div><br></div><div>    >>> reshaped_times = interesting_times[:, np.newaxis]</div><div>    >>> group_at_it = arr.vindex[reshaped_times, [2, 5]]</div><div><br></div><div><br></div><div>One may find, that for example our data is corrupt in some places.</div><div>So, we need to replace these values by zero (or anything else) for these</div><div>times. The first column may for example give the necessary information,</div><div>so that changing the values becomes easy remembering boolean indexing::</div><div><br></div><div>    >>> bad_data = arr[:, 0] > 0.5</div><div>    >>> arr[bad_data, :] = 0  # (corrupts further examples)</div><div><br></div><div>Again, however, the columns may need to be handled more individually (but in</div><div>groups), and the ``oindex`` attribute works well::</div><div><br></div><div>    >>> arr.oindex[bad_data, [2, 5]] = 0</div><div><br></div><div>Note that it would be very hard to do this using legacy fancy indexing.</div><div>The only way would be to create an integer array first::</div><div><br></div><div>    >>> bad_data_indx = np.nonzero(bad_data)[0]</div><div>    >>> bad_data_indx_reshaped = bad_data_indx[:, np.newaxis]</div><div>    >>> arr[bad_data_indx_reshaped, [2, 5]]</div><div><br></div><div>In any case we can use only ``oindex`` to do all of this without getting</div><div>into any trouble or confused by the whole complexity of advanced indexing.</div><div><br></div><div>But, some new features are added to the data acquisition. Different sensors</div><div>have to be used depending on the times. Let us assume we already have</div><div>created an array of indices::</div><div><br></div><div>    >>> correct_sensors = np.random.randint(10, size=(100, 2))</div><div><br></div><div>Which lists for each time the two correct sensors in an ``(N, 2)`` array.</div><div><br></div><div>A first try to achieve this may be ``arr[:, correct_sensors]`` and this does</div><div>not work. It should be clear quickly that slicing cannot achieve the desired</div><div>thing. But hopefully users will remember that there is ``vindex`` as a more</div><div>powerful and flexible approach to advanced indexing.</div><div>One may, if trying ``vindex`` randomly, be confused about::</div><div><br></div><div>    >>> new_arr = arr.vindex[:, correct_sensors]</div><div><br></div><div>which is neither the same, nor the correct result (see transposing rules)!</div><div>This is because slicing works still the same in ``vindex``. However, reading</div><div>the documentation and examples, one can hopefully quickly find the desired</div><div>solution::</div><div><br></div><div>    >>> rows = np.arange(len(arr))</div><div>    >>> rows = rows[:, np.newaxis]  # make shape fit with correct_sensors</div><div>    >>> new_arr = arr.vindex[rows, correct_sensors]</div><div>    </div><div>At this point we have left the straight forward world of ``oindex`` but can</div><div>do random picking of any element from the array. Note that in the last example</div><div>a method such as mentioned in the ``Related Questions`` section could be more</div><div>straight forward. But this approach is even more flexible, since ``rows``</div><div>does not have to be a simple ``arange``, but could be ``intersting_times``::</div><div><br></div><div>    >>> interesting_times = np.array([0, 4, 8, 9, 10])</div><div>    >>> correct_sensors_at_it = correct_sensors[interesting_times, :]</div><div>    >>> interesting_times_reshaped = interesting_times[:, np.newaxis]</div><div>    >>> new_arr_it = arr[interesting_times_reshaped, correct_sensors_at_it]</div><div><br></div><div>Truly complex situation would arise now if you would for example pool ``L``</div><div>experiments into an array shaped ``(L, N, D)``. But for ``oindex`` this should</div><div>not result into surprises. ``vindex``, being more powerful, will quite</div><div>certainly create some confusion in this case but also cover pretty much all</div><div>eventualities.</div><div><br></div><div><br></div><div>Copyright</div><div>---------</div><div><br></div><div>This document is placed under the CC0 1.0 Universell (CC0 1.0) Public Domain Dedication [1]_.</div><div><br></div><div><br></div><div>References and Footnotes</div><div>------------------------</div><div><br></div><div>.. [1] To the extent possible under law, the person who associated CC0 </div><div>   with this work has waived all copyright and related or neighboring</div><div>   rights to this work. The CC0 license may be found at</div><div>   <a href="https://creativecommons.org/publicdomain/zero/1.0/">https://creativecommons.org/publicdomain/zero/1.0/</a></div><div>.. [2] e.g., see NEP 18,</div><div>   <a href="http://www.numpy.org/neps/nep-0018-array-function-protocol.html">http://www.numpy.org/neps/nep-0018-array-function-protocol.html</a></div><div><br></div></div>