Mailman 3 November 2015 - NumPy-Discussion

Inconsistent/unexpected indexing semantics
by Lluís Vilanova 01 Dec '15

01 Dec '15

Hi, TL;DR: There's a pending pull request deprecating some behaviour I find unexpected. Does anyone object? Some time ago I noticed that numpy yields unexpected results in some very specific cases. An array can be used to index multiple elements of a single dimension: >>> a = np.arange(8).reshape((2,2,2)) >>> a[ np.array([[0], [0]]) ] array([[[[0, 1], [2, 3]]], [[[0, 1], [2, 3]]]]) Nonetheless, if a list is used instead, it is (unexpectedly) transformed into a tuple, resulting in indexing across multiple dimensions: >>> a[ [[0], [0]] ] array([[0, 1]]) I.e., it is interpeted as: >>> a[ [0], [0] ] array([[0, 1]]) Or what is the same: >>> a[( [0], [0] )] array([[0, 1]]) I've been informed that there's a pending pull request that deprecates this behaviour [1], which could in the future be reverted to what is expected (at least what I expect) from the documents (except for an obscure note in [2]). The discussion leading to this mail can be found here [3]. [1] https://github.com/numpy/numpy/pull/4434 [2] http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html#advanced-ind… [3] https://github.com/numpy/numpy/issues/6564 Thanks, Lluis -- "And it's much the same thing with knowledge, for whenever you learn something new, the whole world becomes that much richer." -- The Princess of Pure Reason, as told by Norton Juster in The Phantom Tollbooth

3 2

Python development on fedora 23.
by Charles R Harris 30 Nov '15

30 Nov '15

Hi Fedora users, Python distutils on fedora 23 is configured for hardening, hence `redhat-rpm-config` is a dependency if you want to build numpy, scipy, etc. The symptom is a "broken toolchain" error. A bug is open for this and it might get fixed on the Python end, but I don't expect anything soon. Chuck

1 0

Indexing NEP draft
by Sebastian Berg 30 Nov '15

30 Nov '15

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 https://github.com/numpy/numpy/pull/6256/files?short_path=01e4dd9#diff-01e4… 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. Best, Sebastian ================================================================================== 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 problem 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 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 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. Motivation ========== Old style advanced indexing with multiple array (boolean or integer) indices, also called "fancy indexing", tends to be very confusing for new users. While fancy (or legacy) indexing is useful in many cases one would naively 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 arrays into a single index over multiple dimensions. There are three main points 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 useful. 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 ``arr.oindex[indices]``. 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 allow legacy fancy indexing indefinetly. This also gives a simple way to update fancy indexing code making deprecations to vanilla indexing easier. 6. Vanilla indexing ``arr[...]`` could return an error for ambiguous cases. For the beginning, this probably means cases where ``arr[ind]`` and ``arr.oindex[ind]`` return different results gives deprecation warnings. 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 necessary. This helps finding bugs since a too high dimensional array can be caught. (I am in favor for this, but doubt we should think about this for vanilla indexing.) 2. The names ``oindex`` and ``vindex`` are just suggestions at the time of writing this, another name NumPy has used for something like ``oindex`` is ``np.ix_``. See also below. 3. It would be possible to limit the use of boolean indices in ``vindex``, assuming that they are rare and to some degree special. (This would make implementation simpler, but I do not see a big reason.) 4. ``oindex`` and ``vindex`` could always return copies, even when no array operation occurs. One argument for using the same rules is that this way ``oindex`` can be used as a general index replacement. (There is likely no big reason for this, however, there is one reason: ``arr.vindex[array_scalar, ...]`` can occur, where ``arr_scalar`` should be a 0-D array. Copying always "fixes" the possible inconsistency.) 5. The final state to morph 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 now. 6. Proposed changes to vanilla indexing could be postponed indefinetly or not taken in order to not break or force fixing of existing code bases. 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 indexing operations and create warnings (and eventually deprecate) ambiguous direct indexing operations on arrays. Alternative Names ================= Possible names suggested (more suggestions will be added). ============== ======== ======= **Orthogonal** oindex oix **Vectorized** vindex fix **Legacy** l/findex ============== ======== ======= 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. 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 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 --------------------- 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 sense (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 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] # 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] (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] (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 forward since the axis on which to pick elements has to be supplied. One plausible solution would be to create a function (calling it pick here for simplicity):: np.pick(arr, index_arr, axis=axis) where ``index_arr`` has to be the same shape as ``arr`` except along ``axis``. One could imagine that this can be useful together with other indexing types, but such a function may be sufficient and extra information needed seems easier to pass using a function convention. Another option would be to allow an argument 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 dimensional.) Another solution is to provide functions or features to the ``arg*``functions 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)`` 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 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``:: >>> 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.

1 1

understanding buffering done when broadcasting
by Eli Bendersky 28 Nov '15

28 Nov '15

Hello, I'm trying to understand the buffering done by the Numpy iterator interface (the new post 1.6-one) when running ufuncs on arrays that require broadcasting. Consider this simple case: In [35]: m = np.arange(16).reshape(4,4) In [37]: n = np.arange(4) In [39]: m + n Out[39]: array([[ 0, 2, 4, 6], [ 4, 6, 8, 10], [ 8, 10, 12, 14], [12, 14, 16, 18]]) If I instrument Numpy (setting NPY_IT_DBG_TRACING and such), I see that when the add() ufunc is called, 'n' is copied into a temporary buffer by the iterator. The ufunc then gets the buffer as its data. My question is: why is this buffering needed? It seems wasteful, since no casting is required here, no special alignment problems and also 'n' is contiguously laid out in memory. It seems that it would be more efficient to just use 'n' in the ufunc instead of passing in the buffer. What am I missing? Thanks in advance, Eli

2 4

ENH: added vector operators: divergence, curl and laplacian #6727
by Antonio Lara 27 Nov '15

27 Nov '15

Hello, I have corrected the errors in my previous pull request that includes the new functions divergence, curl and laplacian. https://github.com/numpy/numpy/pull/6727 Thank you, Antonio

1 0

New functions added in pull request
by Antonio Lara 26 Nov '15

26 Nov '15

Hello, I have added three new functions to the file function_base.py in the numpy/lib folder. These are divergence, curl and laplacian (for the moment, laplacian of a scalar field, maybe in the future I will try laplacian for a vector field). The calculation method is based in the existing one for numpy.gradient, with central differences. The changes are in this pull request: https://github.com/numpy/numpy/pull/6727 Thank you, Antonio

2 1

ENH: Add the function 'expand_view'
by John Kirkham 26 Nov '15

26 Nov '15

Takes an array and tacks on arbitrary dimensions on either side, which is returned as a view always. Here are the relevant features: * Creates a view of the array that has the dimensions before and after tacked on to it. * Takes the before and after arguments independent of each other and the current shape. * Allows for read and write access to the underlying array. To see an example of what this would look like, see this PR ( https://github.com/numpy/numpy/pull/6713 ).

4 3

Re: [Numpy-discussion] Misleading/erroneous TypeError message
by Peter Creasey 25 Nov '15

25 Nov '15

> > I just upgraded my numpy and started to received a TypeError from one of > > my codes that relied on the old, less strict, casting behaviour. The error > > message, however, left me scratching my head when trying to debug something > > like this: > > > > >>> a = array([0],dtype=uint64) > > >>> a += array([1],dtype=int64) > > TypeError: Cannot cast ufunc add output from dtype('float64') to > > dtype('uint64') with casting rule 'same_kind' > > > > Where does the 'float64' come from?!?! > > > > The combination of uint64 and int64 leads to promotion to float64 as the > best option for the combination of signed and unsigned. To fix things, you > can either use `np.add` with an output argument and `casting='unsafe'` or > just be careful about using unsigned types. Thanks for the quick response. I understand there are reasons for the promotion to float64 (although my expectation would usually be that Numpy is going to follow C conventions), however the I found the error a little unhelpful. In particular Numpy is complaining about a dtype (float64) that it silently promoted to, rather than the dtype that the user provided, which generally seems like a bad idea. Could Numpy somehow complain about the original dtypes in this case? Or at least give a warning about the first promotion (e.g. loss of precision)? Peter

1 0

Misleading/erroneous TypeError message
by Peter Creasey 25 Nov '15

25 Nov '15

Hi, I just upgraded my numpy and started to received a TypeError from one of my codes that relied on the old, less strict, casting behaviour. The error message, however, left me scratching my head when trying to debug something like this: >>> a = array([0],dtype=uint64) >>> a += array([1],dtype=int64) TypeError: Cannot cast ufunc add output from dtype('float64') to dtype('uint64') with casting rule 'same_kind' Where does the 'float64' come from?!?! Peter PS Thanks for all the great work guys, numpy is a fantastic tool and has been a lot of help to me over the years!

2 1

record array performance issue / bug
by G Jones 23 Nov '15

23 Nov '15

Hi, Using the latest numpy from anaconda (1.10.1) on Python 2.7, I found that the following code works OK if npackets = 2, but acts bizarrely if npackets is large (2**12): ----------- npackets = 2**12 dlen=2048 PacketType = np.dtype([('timestamp','float64'), ('pkts',np.dtype(('int8',(npackets,dlen)))), ('data',np.dtype(('int8',(npackets*dlen,)))), ]) b = np.zeros((1,),dtype=PacketType) b['timestamp'] # Should return array([0.0]) ---------------- Specifically, if npackets is large, i.e. 2**12 or 2**16, trying to access b['timestamp'] results in 100% CPU usage while the memory consumption is increasing by hundreds of MB per second. When I interrupt, I find the traceback in numpy/core/_internal.pyc : _get_all_field_offsets Since it seems to work for small values of npackets, I suspect that if I had the memory and time, the access to b['timestamp'] would eventually return, so I think the issue is that the algorithm doesn't scale well with record dtypes made up of lots of bytes. Looking on Github, I can see this code has been in flux recently, but I can't quite tell if the issue I'm seeing is addressed by the issues being discussed and tackled there. Thanks, Glenn

2 1