Mailman 3 June 2018 - NumPy-Discussion

NEP 21: Simplified and explicit advanced indexing
by Stephan Hoyer June 27, 2018

June 27, 2018

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

8 36

NEP 21: Simplified and explicit advanced indexing
by Kevin Sheppard June 26, 2018

June 26, 2018

This seems like progress and a clear method to outer indexing will help many users. As for names, I prefer .ox and .vx for shorthand of .oindex and .vindex. I don’t like the .ox_ or .o_ syntax. Before any deprecation warnings or any other warnings are added it would be helpful to have some way to set a flag on Python to show some sort of HiddenDeprecationWarning (or OnlyShowIfFlagPassesDeprecationWarning) that would automatically be filtered by default but could be shown if someone was interested. This will allow library writers to see problems before any start showing up for users. These could then be promoted to Visible or Future later.

1 0

Remove sctypeNA and typeNA from numpy core
by Matti Picus June 21, 2018

June 21, 2018

numpy.core has many ways to catalogue dtype names: sctypeDict, typeDict (which is precisely sctypeDict), typecodes, and typename. We also generate sctypeNA and typeNA but, as issue 11241 shows, it is sometimes wrong. They are also not documented and never used inside numpy. Instead of fixing it, I propose to remove sctypeNA and typeNA. Any thoughts or objections? Matti

3 3

Forcing new dimensions to appear at front in advanced indexing
by Michael Lamparski June 20, 2018

June 20, 2018

Hi all, So, in advanced indexing, numpy decides where to put new axes based on whether the "advanced indices" are all next to each other. >>> np.random.random((3,4,5,6,7,8))[:, [[0,0],[0,0]], 1, :].shape (3, 2, 2, 6, 7, 8) >>> np.random.random((3,4,5,6,7,8))[:, [[0,0],[0,0]], :, 1].shape (2, 2, 3, 5, 7, 8) In creating a wrapper type around arrays, I'm finding myself needing to suppress this behavior, so that the new axes consistently appear in the front. I thought of a dumb hat trick: def index(x, indices): return x[(True, None) + indices] Which certainly gets the new dimensions where I want them, but it introduces a ghost dimension of 1 (and sometimes two such dimensions!) in a place where I'm not sure I can easily find it. >>> np.random.random((3,4,5,6,7,8))[True, None, 1].shape (1, 1, 4, 5, 6, 7, 8) >>> np.random.random((3,4,5,6,7,8))[True, None, :, [[0,0],[0,0]], 1, :].shape (2, 2, 1, 3, 6, 7, 8) >>> np.random.random((3,4,5,6,7,8))[True, None, :, [[0,0],[0,0]], :, 1].shape (2, 2, 1, 3, 5, 7, 8) any better ideas? --- Michael

2 3

question about array slicing and element assignment
by Emil Sidky June 19, 2018

June 19, 2018

Hello, The following is an example where an array element assignment didn't work as I expected. Create a 6 x 3 matrix: In [70]: a = randn(6,3) In [71]: a Out[71]: array([[ 1.73266816, 0.948849 , 0.69188222], [-0.61840161, -0.03449826, 0.15032552], [ 0.4963306 , 0.77028209, -0.63076396], [-1.92273602, -1.03146536, 0.27744612], [ 0.70736325, 1.54687964, -0.75573888], [ 0.16316043, -0.34814532, 0.3683143 ]]) Create a 3x3 boolean array: In [72]: mask = randn(3,3)>0. In [73]: mask Out[73]: array([[ True, True, True], [False, True, True], [ True, False, True]], dtype=bool) Try to modify elements of "a" with the following line: In [74]: a[(2,3,5),][mask] = 1. No elements are changed in "a": In [75]: a Out[75]: array([[ 1.73266816, 0.948849 , 0.69188222], [-0.61840161, -0.03449826, 0.15032552], [ 0.4963306 , 0.77028209, -0.63076396], [-1.92273602, -1.03146536, 0.27744612], [ 0.70736325, 1.54687964, -0.75573888], [ 0.16316043, -0.34814532, 0.3683143 ]]) Instead try to modify elements of "a" with this line: In [76]: a[::2,][mask] = 1. This time it works: In [77]: a Out[77]: array([[ 1. , 1. , 1. ], [-0.61840161, -0.03449826, 0.15032552], [ 0.4963306 , 1. , 1. ], [-1.92273602, -1.03146536, 0.27744612], [ 1. , 1.54687964, 1. ], [ 0.16316043, -0.34814532, 0.3683143 ]]) Is there a way where I can modify the elements of "a" selected by an expression like "a[(2,3,5),][mask]" ? Thanks , Emil

2 1

ANN: SfePy 2018.2
by Robert Cimrman June 19, 2018

June 19, 2018

I am pleased to announce release 2018.2 of SfePy. Description ----------- SfePy (simple finite elements in Python) is a software for solving systems of coupled partial differential equations by the finite element method or by the isogeometric analysis (limited support). It is distributed under the new BSD license. Home page: http://sfepy.org Mailing list: https://mail.python.org/mm3/mailman3/lists/sfepy.python.org/ Git (source) repository, issue tracker: https://github.com/sfepy/sfepy Highlights of this release -------------------------- - generalized-alpha and velocity Verlet elastodynamics solvers - terms for dispersion in fluids - caching of reference coordinates for faster repeated use of probes - new wrapper of MUMPS linear solver for parallel runs For full release notes see http://docs.sfepy.org/doc/release_notes.html#id1 (rather long and technical). Cheers, Robert Cimrman --- Contributors to this release in alphabetical order: Robert Cimrman Lubos Kejzlar Vladimir Lukes Matyas Novak

1 0

A little about XND
by Travis Oliphant June 19, 2018

June 19, 2018

Hi everyone, I'm glad I'm able to contribute back to this discussion thread. I wanted to post a quick message to this group to make sure there is no mis-information about XND which has finally reached the point where it can be experimented with (http://xnd.io) and commented on. XND came out of thoughts and conversations we had at Continuum (now Anaconda) when thinking about cross-language array computing and how to enable improved features for high-level users in many languages (including Python, R, Ruby, Node, Scala, Rust, Go, etc.). Technically there are three projects that make up XND (thus the name Plures for the Github organization). All of these projects have a C-library and then a high-level interface (right now we only have resources to develop the Python interface but would love to see support for other languages). xnd (libxnd) is the typed container. ndtypes (libndtypes) is the (datashape-like) type system with a grammar, parser, and type matcher. gumath (libgumath) are generalized ufuncs which represent the entire function system on xnd. We will be talking more about XND in the coming months and years, but for the purposes of this list, I wanted to make it clear that 1) XND is not trying to replace NumPy. XND is a low-level library and intended to be such. It would be most welcome if someday NumPy uses XND. We understand this may be a while and certainly not before NumPy 2.0 or 3.0. 2) Our initial target users are Numba, pandas, Dask, xarray, and other higher-level objects at the moment. We are eagerly searching for integration opportunities to connect more developers (or advanced users) to xnd before making more progress. 3) We do discuss array-like things in the public channels. NumPy users and developers are welcome in those channels. Everything is done in public including the weekly meeting which anyone can attend: Weekly meeting: meet.google.com/heo-fmow-omz Live discussions: https://gitter.im/Plures/xnd for the libraries themselves https://gitter.im/Plures/xnd-ml for integrations. Issues and PRs: https://github.com/plures --- under the various projects. 4) We are thinking about adding a custom-dtype to NumPy that uses xnd and would be happy for anyone's help on that project. 5) We are in the early stages of exploring a high-level array interface (using the ideas of MoA and the Psi Calculus with Lenore Mullen who worked on APL years ago). Likely the first place this will see some initial progress is in an ND Sparse array that uses XND. We welcome participation and input from all. Stefan Krah has written the majority of the code and so we tend to respect his point of view. Pearu Peterson (of f2py and SciPy fame) has made some useful contributions recently. Stefan and I have been talking roughly weekly for a couple of years and so some of the problems currently there, I am certainly responsible for. Two of our immediate goals are to work with the Numba team to get support for ndtypes in Numba and allow Numba to use libgumath in no-python mode. I look forward to continuing the conversation with any of you who want to participate. Perhaps some of us can meet up during NumPy sprints to discuss more. XND is also currently looking for funding and time from interested parties to continue its development. -Travis

3 7

NEP: Random Number Generator Policy
by Robert Kern June 17, 2018

June 17, 2018

As promised distressingly many months ago, I have written up a NEP about relaxing the stream-compatibility policy that we currently have. https://github.com/numpy/numpy/pull/11229 https://github.com/rkern/numpy/blob/nep/rng/doc/neps/nep-0019-rng-policy.rst I particularly invite comment on the two lists of methods that we still would make strict compatibility guarantees for. --- ============================== Random Number Generator Policy ============================== :Author: Robert Kern <robert.kern(a)gmail.com> :Status: Draft :Type: Standards Track :Created: 2018-05-24 Abstract -------- For the past decade, NumPy has had a strict backwards compatibility policy for the number stream of all of its random number distributions. Unlike other numerical components in ``numpy``, which are usually allowed to return different when results when they are modified if they remain correct, we have obligated the random number distributions to always produce the exact same numbers in every version. The objective of our stream-compatibility guarantee was to provide exact reproducibility for simulations across numpy versions in order to promote reproducible research. However, this policy has made it very difficult to enhance any of the distributions with faster or more accurate algorithms. After a decade of experience and improvements in the surrounding ecosystem of scientific software, we believe that there are now better ways to achieve these objectives. We propose relaxing our strict stream-compatibility policy to remove the obstacles that are in the way of accepting contributions to our random number generation capabilities. The Status Quo -------------- Our current policy, in full: A fixed seed and a fixed series of calls to ``RandomState`` methods using the same parameters will always produce the same results up to roundoff error except when the values were incorrect. Incorrect values will be fixed and the NumPy version in which the fix was made will be noted in the relevant docstring. Extension of existing parameter ranges and the addition of new parameters is allowed as long the previous behavior remains unchanged. This policy was first instated in Nov 2008 (in essence; the full set of weasel words grew over time) in response to a user wanting to be sure that the simulations that formed the basis of their scientific publication could be reproduced years later, exactly, with whatever version of ``numpy`` that was current at the time. We were keen to support reproducible research, and it was still early in the life of ``numpy.random``. We had not seen much cause to change the distribution methods all that much. We also had not thought very thoroughly about the limits of what we really could promise (and by “we” in this section, we really mean Robert Kern, let’s be honest). Despite all of the weasel words, our policy overpromises compatibility. The same version of ``numpy`` built on different platforms, or just in a different way could cause changes in the stream, with varying degrees of rarity. The biggest is that the ``.multivariate_normal()`` method relies on ``numpy.linalg`` functions. Even on the same platform, if one links ``numpy`` with a different LAPACK, ``.multivariate_normal()`` may well return completely different results. More rarely, building on a different OS or CPU can cause differences in the stream. We use C ``long`` integers internally for integer distribution (it seemed like a good idea at the time), and those can vary in size depending on the platform. Distribution methods can overflow their internal C ``longs`` at different breakpoints depending on the platform and cause all of the random variate draws that follow to be different. And even if all of that is controlled, our policy still does not provide exact guarantees across versions. We still do apply bug fixes when correctness is at stake. And even if we didn’t do that, any nontrivial program does more than just draw random numbers. They do computations on those numbers, transform those with numerical algorithms from the rest of ``numpy``, which is not subject to so strict a policy. Trying to maintain stream-compatibility for our random number distributions does not help reproducible research for these reasons. The standard practice now for bit-for-bit reproducible research is to pin all of the versions of code of your software stack, possibly down to the OS itself. The landscape for accomplishing this is much easier today than it was in 2008. We now have ``pip``. We now have virtual machines. Those who need to reproduce simulations exactly now can (and ought to) do so by using the exact same version of ``numpy``. We do not need to maintain stream-compatibility across ``numpy`` versions to help them. Our stream-compatibility guarantee has hindered our ability to make improvements to ``numpy.random``. Several first-time contributors have submitted PRs to improve the distributions, usually by implementing a faster, or more accurate algorithm than the one that is currently there. Unfortunately, most of them would have required breaking the stream to do so. Blocked by our policy, and our inability to work around that policy, many of those contributors simply walked away. Implementation -------------- We propose first freezing ``RandomState`` as it is and developing a new RNG subsystem alongside it. This allows anyone who has been relying on our old stream-compatibility guarantee to have plenty of time to migrate. ``RandomState`` will be considered deprecated, but with a long deprecation cycle, at least a few years. Deprecation warnings will start silent but become increasingly noisy over time. Bugs in the current state of the code will *not* be fixed if fixing them would impact the stream. However, if changes in the rest of ``numpy`` would break something in the ``RandomState`` code, we will fix ``RandomState`` to continue working (for example, some change in the C API). No new features will be added to ``RandomState``. Users should migrate to the new subsystem as they are able to. Work on a proposed `new PRNG subsystem <https://github.com/bashtage/randomgen>`_ is already underway. The specifics of the new design are out of scope for this NEP and up for much discussion, but we will discuss general policies that will guide the evolution of whatever code is adopted. First, we will maintain API source compatibility just as we do with the rest of ``numpy``. If we *must* make a breaking change, we will only do so with an appropriate deprecation period and warnings. Second, breaking stream-compatibility in order to introduce new features or improve performance will be *allowed* with *caution*. Such changes will be considered features, and as such will be no faster than the standard release cadence of features (i.e. on ``X.Y`` releases, never ``X.Y.Z``). Slowness is not a bug. Correctness bug fixes that break stream-compatibility can happen on bugfix releases, per usual, but developers should consider if they can wait until the next feature release. We encourage developers to strongly weight user’s pain from the break in stream-compatibility against the improvements. One example of a worthwhile improvement would be to change algorithms for a significant increase in performance, for example, moving from the `Box-Muller transform <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>`_ method of Gaussian variate generation to the faster `Ziggurat algorithm <https://en.wikipedia.org/wiki/Ziggurat_algorithm>`_. An example of an unworthy improvement would be tweaking the Ziggurat tables just a little bit. Any new design for the RNG subsystem will provide a choice of different core uniform PRNG algorithms. We will be more strict about a select subset of methods on these core PRNG objects. They MUST guarantee stream-compatibility for a minimal, specified set of methods which are chosen to make it easier to compose them to build other distributions. Namely, * ``.bytes()`` * ``.random_uintegers()`` * ``.random_sample()`` Furthermore, the new design should also provide one generator class (we shall call it ``StableRandom`` for discussion purposes) that provides a slightly broader subset of distribution methods for which stream-compatibility is *guaranteed*. The point of ``StableRandom`` is to provide something that can be used in unit tests so projects that currently have tests which rely on the precise stream can be migrated off of ``RandomState``. For the best transition, ``StableRandom`` should use as its core uniform PRNG the current MT19937 algorithm. As best as possible, the API for the distribution methods that are provided on ``StableRandom`` should match their counterparts on ``RandomState``. They should provide the same stream that the current version of ``RandomState`` does. Because their intended use is for unit tests, we do not need the performance improvements from the new algorithms that will be introduced by the new subsystem. The list of ``StableRandom`` methods should be chosen to support unit tests: * ``.randint()`` * ``.uniform()`` * ``.normal()`` * ``.standard_normal()`` * ``.choice()`` * ``.shuffle()`` * ``.permutation()`` Not Versioning -------------- For a long time, we considered that the way to allow algorithmic improvements while maintaining the stream was to apply some form of versioning. That is, every time we make a stream change in one of the distributions, we increment some version number somewhere. ``numpy.random`` would keep all past versions of the code, and there would be a way to get the old versions. Proposals of how to do this exactly varied widely, but we will not exhaustively list them here. We spent years going back and forth on these designs and were not able to find one that sufficed. Let that time lost, and more importantly, the contributors that we lost while we dithered, serve as evidence against the notion. Concretely, adding in versioning makes maintenance of ``numpy.random`` difficult. Necessarily, we would be keeping lots of versions of the same code around. Adding a new algorithm safely would still be quite hard. But most importantly, versioning is fundamentally difficult to *use* correctly. We want to make it easy and straightforward to get the latest, fastest, best versions of the distribution algorithms; otherwise, what's the point? The way to make that easy is to make the latest the default. But the default will necessarily change from release to release, so the user’s code would need to be altered anyway to specify the specific version that one wants to replicate. Adding in versioning to maintain stream-compatibility would still only provide the same level of stream-compatibility that we currently do, with all of the limitations described earlier. Given that the standard practice for such needs is to pin the release of ``numpy`` as a whole, versioning ``RandomState`` alone is superfluous. Discussion ---------- - https://mail.python.org/pipermail/numpy-discussion/2018-January/077608.html - https://github.com/numpy/numpy/pull/10124#issuecomment-350876221 Copyright --------- This document has been placed in the public domain. -- Robert Kern

9 51

Dropping Python 3.4 support for NumPy 1.16
by Charles R Harris June 14, 2018

June 14, 2018

Hi All, I think NumPy 1.16 would be a good time to drop Python 3.4 support. We will want to do that anyway once we drop 2.7 so that we will only be using recent Windows compilers, and with Python 3.7 due at the end of the month I think supporting 3.5-7 for 1.16 should be sufficient. Thoughts? Chuck

7 6

Circle CI moving from 1.0 to 2.0
by Matti Picus June 14, 2018

June 14, 2018

I stumbled across this notice (only seems to appear in a failed build) "This project is currently running on CircleCI 1.0 which will no longer be supported after August 31, 2018. Please start migrating this project to CircleCI 2.0 <https://circleci.com/docs/2.0/migration/>." Here is the original link https://circleci.com/gh/numpy/numpy/2080 Is this an artifact that can be ignored or do we need to migrate, if so has anyone already done it for their project? Matti

2 1