Mailman 3 July 2018 - NumPy-Discussion

Allow moveaxis() to take strings as axis order
by Manuel Schölling 08 Jul '18

08 Jul '18

Hi, I have opened a pull request [1] for enhancing moveaxis() a bit, so it can take strings as arguments to permute the axis order. The new feature is best described by a short example. All these calls to moveaxis() will do the same thing: A = np.zeros((0, 1, 2, 3)) # reverse order np.moveaxis(A, [0, 1, 2, 3], [3, 2, 1, 0]) np.moveaxis(A, 'ijkl', 'lkji') np.moveaxis(A, 'wxyz', 'zyxw') # more complicated axis permutation np.moveaxis(A, [0, 1, 2, 3], [3, 1, 2, 0]) np.moveaxis(A, 'ijkl', 'ljki') np.moveaxis(A, 'wxyz', 'zxyw') It has been mentioned by Eric Wieser that np.einsum('ijkl->lkji', A) basically does the same thing, although it is not really obvious to the user that is not used Einstein's sum convention what this call does, because the function call says it would perform a 'sum' operation, but actually only permutes axes. It was also mentioned that the Zen of Python says > There should be one-- and preferably only one --obvious way to do it. and if moveaxis() supports multiple ways to perform its operation, there could be a cognitive bias if the developer is only used to e.g. the string arguments and suddenly see the array arguments (or vice versa). Despite these arguments I think the string arguments makes calls to moveaxis() much more obvious than integer arguments, especially if ndim > 3. How does the mailing list feel about this pull request? Bye, Manuel [1] https://github.com/numpy/numpy/pull/11504

1 0

Allowing broadcasting of code dimensions in generalized ufuncs
by Marten van Kerkwijk 03 Jul '18

03 Jul '18

Hi All, Following on a PR combining the ability to provide fixed and flexible dimensions [1] (useful for, e.g., 3-vector input with a signature like `(3),(3)->(3)`, and for `matmul`, resp.; based on earlier PRs by Jaime [2] and Matt (Picus) [3]), I've now made a PR with a further enhancement, which allows one can indicate that a core dimension can be broadcast [4]. A particular use case is `all_equal`, a new function suggested in a stalled PR by Matt (Harrigan) [5], which compares two arrays axis-by-axis, but short-circuits if a non-equality is found (unlike what is the case if one does `(a==b).all(axis)`). One thing that would be obviously useful for a routine like `all_equal` is to be able to provide an array as one argument and a constant as another, i.e., if the core dimensions can be broadcast if needed, just like they are in `(a==b).all(axis)`. This is currently not possible: with its signature of `(n),(n)->()`, the two arrays have to have the same trailing size. My PR provides the ability to indicate in the signature that a core dimension can be broadcast, by using a suffix of "|1". Thus, the signature of `all_equal` would become: ``` (n|1),(n|1)->() ``` Comments most welcome (yes, even on the notation - though I think it is fairly self-explanatory)! Marten p.s. There are some similarities to the new "flexible" dimensions implemented for `matmul` [1], but also differences. In particular, for a signature of `(n?),(n?)->()`, one could also pass in an array of trailing size n and a constant, but it would not be possible to pass in an array with trailing size 1: the dimensions with the same name have to be either present and the same or absent. In contrast, for broadcasting, dimensions with the same name can have trailing size n, size 1, or be absent (in which case they count as 1). For broadcasting, any output dimensions with the same name are never affected, while for flexible dimensions those are removed. [1] https://github.com/numpy/numpy/pull/11175 [2] https://github.com/numpy/numpy/pull/5015 [3] https://github.com/numpy/numpy/pull/11132 [4] https://github.com/numpy/numpy/pull/11179 [5] https://github.com/numpy/numpy/pull/8528

10 42

PEP 574 - zero-copy pickling with out of band data
by Antoine Pitrou 03 Jul '18

03 Jul '18

Hello, Some of you might know that I've been working on a PEP in order to improve pickling performance of large (or huge) data. The PEP, numbered 574 and titled "Pickle protocol 5 with out-of-band data", allows participating data types to be pickled without any memory copy. https://www.python.org/dev/peps/pep-0574/ The PEP already has an implementation, which is backported as an independent PyPI package under the name "pickle5". https://pypi.org/project/pickle5/ I also have a working patch updating PyArrow to use the PEP-defined extensions to allow for zero-copy pickling of Arrow arrays - without breaking compatibility with existing usage: https://github.com/apache/arrow/pull/2161 Still, it is obvious one the primary targets of PEP 574 is Numpy arrays, as the most prevalent datatype in the Python scientific ecosystem. I'm personally satisfied with the current state of the PEP, but I'd like to have feedback from Numpy core maintainers. I haven't tried (yet?) to draft a Numpy patch to add PEP 574 support, since that's likely to be more involved due to the complexity of Numpy and due to the core being written in C. Therefore I would like some help evaluating whether the PEP is likely to be a good fit for Numpy. Regards Antoine.

6 10

Corrupted installation?
by Tom Cook 03 Jul '18

03 Jul '18

I have a remote sensing platform that uses Numpy on a raspberry pi. There are about a dozen of these installed. Three days ago, the Python part of the software started crashing during startup, with this exception message: Traceback (most recent call last): File "/home/azi/board/sensor_logger.py", line 16, in <module> from board.measurement.data import Expression, MeasurementData File "/home/azi/board/measurement/data.py", line 2, in <module> import numpy as np File "/usr/lib/python3/dist-packages/numpy/__init__.py", line 153, in <module> from . import add_newdocs File "/usr/lib/python3/dist-packages/numpy/add_newdocs.py", line 13, in <module> from numpy.lib import add_newdoc File "/usr/lib/python3/dist-packages/numpy/lib/__init__.py", line 8, in <module> from .type_check import * File "/usr/lib/python3/dist-packages/numpy/lib/type_check.py", line 11, in <module> import numpy.core.numeric as _nx File "/usr/lib/python3/dist-packages/numpy/core/__init__.py", line 7, in <module> from . import umath SystemError: initialization of umath raised unreported exception I'm running numpy 1.8.2 on python 3.4.2 on a Raspberry Pi running Raspbian Jessie. I realise these are somewhat old, but as the software is running in the field (and I mean literally in a field - on another continent) I'm not in a great position to update them. Can anyone give me a pointer on how to debug this? As far as I can tell, umath.cpython-34m-arm-linux-gnueabihf.so has not changed recently. In fact, as far as I can tell, nothing changed around the time this stopped working. But I guess something must have! Regards, Tom

3 2

Unable to build Numpy on Mac OSX
by Aman Agarwal 03 Jul '18

03 Jul '18

Hello Numpy, Can you please help us in building the numpy source code I am build on my Mac OS X 10.11.6. Step to reproduce : 1. Clone the package from git numpy/numpy.git 2. build it $> python setup.py build_ext -i 3. test it $> python -c 'import numpy; numpy.test()' on running to step 3 I am getting the following error - https://bpaste.net/show/b14b9e380f26 Please let me know if you need more information on this. Regards, Aman Agarwal

1 0

NEP 22 — Duck typing for NumPy arrays – high level overview
by Nathaniel Smith 03 Jul '18

03 Jul '18

Hi all, Here's a NEP that Stephan and I wrote (actually a few months ago, but then I was super slow, thank you Stephan for your patience). It tries to lay out a high-level overview of where we're trying to go with all this duck array stuff, and some general guidelines that we've gradually come around to after lots of time scowling at whiteboards in frustration. -n ---------- =========================================================== NEP 22 — Duck typing for NumPy arrays – high level overview =========================================================== :Author: Stephan Hoyer <shoyer(a)google.com>, Nathaniel J. Smith <njs(a)pobox.com> :Status: Draft :Type: Informational :Created: 2018-03-22 Abstract -------- We outline a high-level vision for how NumPy will approach handling “duck arrays”. This is an Informational-class NEP; it doesn’t prescribe full details for any particular implementation. In brief, we propose developing a number of new protocols for defining implementations of multi-dimensional arrays with high-level APIs matching NumPy. Detailed description -------------------- Traditionally, NumPy’s ``ndarray`` objects have provided two things: a high level API for expression operations on homogenously-typed, arbitrary-dimensional, array-structured data, and a concrete implementation of the API based on strided in-RAM storage. The API is powerful, fairly general, and used ubiquitously across the scientific Python stack. The concrete implementation, on the other hand, is suitable for a wide range of uses, but has limitations: as data sets grow and NumPy becomes used in a variety of new environments, there are increasingly cases where the strided in-RAM storage strategy is inappropriate, and users find they need sparse arrays, lazily evaluated arrays (as in dask), compressed arrays (as in blosc), arrays stored in GPU memory, arrays stored in alternative formats such as Arrow, and so forth – yet users still want to work with these arrays using the familiar NumPy APIs, and re-use existing code with minimal (ideally zero) porting overhead. As a working shorthand, we call these “duck arrays”, by analogy with Python’s “duck typing”: a “duck array” is a Python object which “quacks like” a numpy array in the sense that it has the same or similar Python API, but doesn’t share the C-level implementation. This NEP doesn’t propose any specific changes to NumPy or other projects; instead, it gives an overview of how we hope to extend NumPy to support a robust ecosystem of projects implementing and relying upon its high level API. Terminology ~~~~~~~~~~~ “Duck array” works fine as a placeholder for now, but it’s pretty jargony and may confuse new users, so we may want to pick something else for the actual API functions. Unfortunately, “array-like” is already taken for the concept of “anything that can be coerced into an array” (including e.g. list objects), and “anyarray” is already taken for the concept of “something that shares ndarray’s implementation, but has different semantics”, which is the opposite of a duck array (e.g., np.matrix is an “anyarray”, but is not a “duck array”). This is a classic bike-shed so for now we’re just using “duck array”. Some possible options though include: arrayish, pseudoarray, nominalarray, ersatzarray, arraymimic, ... General approach ~~~~~~~~~~~~~~~~ At a high level, duck array support requires working through each of the API functions provided by NumPy, and figuring out how it can be extended to work with duck array objects. In some cases this is easy (e.g., methods/attributes on ndarray itself); in other cases it’s more difficult. Here are some principles we’ve found useful so far: Principle 1: Focus on “full” duck arrays, but don’t rule out “partial” duck arrays ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ We can distinguish between two classes: * “full” duck arrays, which aspire to fully implement np.ndarray’s Python-level APIs and work essentially anywhere that np.ndarray works * “partial” duck arrays, which intentionally implement only a subset of np.ndarray’s API. Full duck arrays are, well, kind of boring. They have exactly the same semantics as ndarray, with differences being restricted to under-the-hood decisions about how the data is actually stored. The kind of people that are excited about making numpy more extensible are also, unsurprisingly, excited about changing or extending numpy’s semantics. So there’s been a lot of discussion of how to best support partial duck arrays. We've been guilty of this ourself. At this point though, we think the best general strategy is to focus our efforts primarily on supporting full duck arrays, and only worry about partial duck arrays as much as we need to to make sure we don't accidentally rule them out for no reason. Why focus on full duck arrays? Several reasons: First, there are lots of very clear use cases. Potential consumers of the full duck array interface include almost every package that uses numpy (scipy, sklearn, astropy, ...), and in particular packages that provide array-wrapping-classes that handle multiple types of arrays, such as xarray and dask.array. Potential implementers of the full duck array interface include: distributed arrays, sparse arrays, masked arrays, arrays with units (unless they switch to using dtypes), labeled arrays, and so forth. Clear use cases lead to good and relevant APIs. Second, the Anna Karenina principle applies here: full duck arrays are all alike, but every partial duck array is partial in its own way: * ``xarray.DataArray`` is mostly a duck array, but has incompatible broadcasting semantics. * ``xarray.Dataset`` wraps multiple arrays in one object; it still implements some array interfaces like ``__array_ufunc__``, but certainly not all of them. * ``pandas.Series`` has methods with similar behavior to numpy, but unique null-skipping behavior. * scipy’s ``LinearOperator``\s support matrix multiplication and nothing else * h5py and similar libraries for accessing array storage have objects that support numpy-like slicing and conversion into a full array, but not computation. * Some classes may be similar to ndarray, but without supporting the full indexing semantics. And so forth. Despite our best attempts, we haven't found any clear, unique way of slicing up the ndarray API into a hierarchy of related types that captures these distinctions; in fact, it’s unlikely that any single person even understands all the distinctions. And this is important, because we have a *lot* of APIs that we need to add duck array support to (both in numpy and in all the projects that depend on numpy!). By definition, these already work for ``ndarray``, so hopefully getting them to work for full duck arrays shouldn’t be so hard, since by definition full duck arrays act like ``ndarray``. It’d be very cumbersome to have to go through each function and identify the exact subset of the ndarray API that it needs, then figure out which partial array types can/should support it. Once we have things working for full duck arrays, we can go back later and refine the APIs needed further as needed. Focusing on full duck arrays allows us to start making progress immediately. In the future, it might be useful to identify specific use cases for duck arrays and standardize narrower interfaces targeted just at those use cases. For example, it might make sense to have a standard “array loader” interface that file access libraries like h5py, netcdf, pydap, zarr, ... all implement, to make it easy to switch between these libraries. But that’s something that we can do as we go, and it doesn’t necessarily have to involve the NumPy devs at all. For an example of what this might look like, see the documentation for `dask.array.from_array <http://dask.pydata.org/en/latest/array-api.html#dask.array.from_array>`__. Principle 2: Take advantage of duck typing ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ``ndarray`` has a very large API surface area:: In [1]: len(set(dir(np.ndarray)) - set(dir(object))) Out[1]: 138 And this is a huge **under**\estimate, because there are also many free-standing functions in NumPy and other libraries which currently use the NumPy C API and thus only work on ``ndarray`` objects. In type theory, a type is defined by the operations you can perform on an object; thus, the actual type of ``ndarray`` includes not just its methods and attributes, but *all* of these functions. For duck arrays to be successful, they’ll need to implement a large proportion of the ``ndarray`` API – but not all of it. (For example, ``dask.array.Array`` does not provide an equivalent to the ``ndarray.ptp`` method, presumably because no-one has ever noticed or cared about its absence. But this doesn’t seem to have stopped people from using dask.) This means that realistically, we can’t hope to define the whole duck array API up front, or that anyone will be able to implement it all in one go; this will be an incremental process. It also means that even the so-called “full” duck array interface is somewhat fuzzily defined at the borders; there are parts of the ``np.ndarray`` API that duck arrays won’t have to implement, but we aren’t entirely sure what those are. And ultimately, it isn’t really up to the NumPy developers to define what does or doesn’t qualify as a duck array. If we want scikit-learn functions to work on dask arrays (for example), then that’s going to require negotiation between those two projects to discover incompatibilities, and when an incompatibility is discovered it will be up to them to negotiate who should change and how. The NumPy project can provide technical tools and general advice to help resolve these disagreements, but we can’t force one group or another to take responsibility for any given bug. Therefore, even though we’re focusing on “full” duck arrays, we *don’t* attempt to define a normative “array ABC” – maybe this will be useful someday, but right now, it’s not. And as a convenient side-effect, the lack of a normative definition leaves partial duck arrays room to experiment. But, we do provide some more detailed advice for duck array implementers and consumers below. Principle 3: Focus on protocols ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Historically, numpy has had lots of success at interoperating with third-party objects by defining *protocols*, like ``__array__`` (asks an arbitrary object to convert itself into an array), ``__array_interface__`` (a precursor to Python’s buffer protocol), and ``__array_ufunc__`` (allows third-party objects to support ufuncs like ``np.exp``). `NEP 16 <https://github.com/numpy/numpy/pull/10706>`_ took a different approach: we need a duck-array equivalent of ``asarray``, and it proposed to do this by defining a version of ``asarray`` that would let through objects which implemented a new AbstractArray ABC. As noted above, we now think that trying to define an ABC is a bad idea for other reasons. But when this NEP was discussed on the mailing list, we realized that even on its own merits, this idea is not so great. A better approach is to define a *method* that can be called on an arbitrary object to ask it to convert itself into a duck array, and then define a version of ``asarray`` that calls this method. This is strictly more powerful: if an object is already a duck array, it can simply ``return self``. It allows more correct semantics: NEP 16 assumed that ``asarray(obj, dtype=X)`` is the same as ``asarray(obj).astype(X)``, but this isn’t true. And it supports more use cases: if h5py supported sparse arrays, it might want to provide an object which is not itself a sparse array, but which can be automatically converted into a sparse array. See NEP <XX, to be written> for full details. The protocol approach is also more consistent with core Python conventions: for example, see the ``__iter__`` method for coercing objects to iterators, or the ``__index__`` protocol for safe integer coercion. And finally, focusing on protocols leaves the door open for partial duck arrays, which can pick and choose which subset of the protocols they want to participate in, each of which have well-defined semantics. Conclusion: protocols are one honking great idea – let’s do more of those. Principle 4: Reuse existing methods when possible ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ It’s tempting to try to define cleaned up versions of ndarray methods with a more minimal interface to allow for easier implementation. For example, ``__array_reshape__`` could drop some of the strange arguments accepted by ``reshape`` and ``__array_basic_getitem__`` could drop all the `strange edge cases <http://www.numpy.org/neps/nep-0021-advanced-indexing.html>`__ of NumPy’s advanced indexing. But as discussed above, we don’t really know what APIs we need for duck-typing ndarray. We would inevitably end up with a very long list of new special methods. In contrast, existing methods like ``reshape`` and ``__getitem__`` have the advantage of already being widely used/exercised by libraries that use duck arrays, and in practice, any serious duck array type is going to have to implement them anyway. Principle 5: Make it easy to do the right thing ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Making duck arrays work well is going to be a community effort. Documentation helps, but only goes so far. We want to make it easy to implement duck arrays that do the right thing. One way NumPy can help is by providing mixin classes for implementing large groups of related functionality at once. ``NDArrayOperatorsMixin`` is a good example: it allows for implementing arithmetic operators implicitly via the ``__array_ufunc__`` method. It’s not complete, and we’ll want more helpers like that (e.g. for reductions). (We initially thought that the importance of these mixins might be an argument for providing an array ABC, since that’s the standard way to do mixins in modern Python. But in discussion around NEP 16 we realized that partial duck arrays also wanted to take advantage of these mixins in some cases, so even if we did have an array ABC then the mixins would still need some sort of separate existence. So never mind that argument.) Tentative duck array guidelines ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ As a general rule, libraries using duck arrays should insist upon the minimum possible requirements, and libraries implementing duck arrays should provide as complete of an API as possible. This will ensure maximum compatibility. For example, users should prefer to rely on ``.transpose()`` rather than ``.swapaxes()`` (which can be implemented in terms of transpose), but duck array authors should ideally implement both. If you are trying to implement a duck array, then you should strive to implement everything. You certainly need ``.shape``, ``.ndim`` and ``.dtype``, but also your dtype attribute should actually be a ``numpy.dtype`` object, weird fancy indexing edge cases should ideally work, etc. Only details related to NumPy’s specific ``np.ndarray`` implementation (e.g., ``strides``, ``data``, ``view``) are explicitly out of scope. A (very) rough sketch of future plans ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The proposals discussed so far – ``__array_ufunc__`` and some kind of ``asarray`` protocol – are clearly necessary but not sufficient for full duck typing support. We expect the need for additional protocols to support (at least) these features: * **Concatenating** duck arrays, which would be used internally by other array combining methods like stack/vstack/hstack. The implementation of concatenate will need to be negotiated among the list of array arguments. We expect to use an ``__array_concatenate__`` protocol like ``__array_ufunc__`` instead of multiple dispatch. * **Ufunc-like functions** that currently aren’t ufuncs. Many NumPy functions like median, percentile, sort, where and clip could be written as generalized ufuncs but currently aren’t. Either these functions should be written as ufuncs, or we should consider adding another generic wrapper mechanism that works similarly to ufuncs but makes fewer guarantees about how the implementation is done. * **Random number generation** with duck arrays, e.g., ``np.random.randn()``. For example, we might want to add new APIs like ``random_like()`` for generating new arrays with a matching shape *and* type – though we'll need to look at some real examples of how these functions are used to figure out what would be helpful. * **Miscellaneous other functions** such as ``np.einsum``, ``np.zeros_like``, and ``np.broadcast_to`` that don’t fall into any of the above categories. * **Checking mutability** on duck arrays, which would imply that they support assignment with ``__setitem__`` and the out argument to ufuncs. Many otherwise fine duck arrays are not easily mutable (for example, because they use some kinds of sparse or compressed storage, or are in read-only shared memory), and it turns out that frequently-used code like the default implementation of ``np.mean`` needs to check this (to decide whether it can re-use temporary arrays). We intentionally do not describe exactly how to add support for these types of duck arrays here. These will be the subject of future NEPs. Copyright --------- This document has been placed in the public domain. -- Nathaniel J. Smith -- https://vorpus.org

1 0

OpenOpt Suite v 0.5627
by Dmitrey Kroshko 02 Jul '18

02 Jul '18

hello, OpenOpt Suite v 0.5627 (OpenOpt, FuncDesigner, SpaceFuncs, DerApproximator) is available for downloading from the link https://app.box.com/s/sikxqjmohtpklqe46ou86t44ma7b27es (unfortunately, I still cannot update it in PYPI entries because of problems with incorrect password) It fixes compatibility issues with latest Python / NumPy / matplotlib versions and has some improvements. Unfortunately, I had no possibilities to perform essential changes since last release in 2015 due to lack of finance support, however I still intend to continue its development. Regards, D.

1 0

Polynomial evaluation inconsistencies
by Maxwell Aifer 02 Jul '18

02 Jul '18

Hi, I noticed some frustrating inconsistencies in the various ways to evaluate polynomials using numpy. Numpy has three ways of evaluating polynomials (that I know of) and each of them has a different syntax: - numpy.polynomial.polynomial.Polynomial <https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.polynomia…>: You define a polynomial by a list of coefficients *in order of increasing degree*, and then use the class’s call() function. - np.polyval <https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.polyval.h…>: Evaluates a polynomial at a point. *First* argument is the polynomial, or list of coefficients *in order of decreasing degree*, and the *second* argument is the point to evaluate at. - np.polynomial.polynomial.polyval <https://docs.scipy.org/doc/numpy-1.12.0/reference/generated/numpy.polynomia…>: Also evaluates a polynomial at a point, but has more support for vectorization. *First* argument is the point to evaluate at, and *second* argument the list of coefficients *in order of increasing degree*. Not only the order of arguments is changed between different methods, but the order of the coefficients is reversed as well, leading to puzzling bugs (in my experience). What could be the reason for this madness? As polyval is a shameless ripoff of Matlab’s function of the same name <https://www.mathworks.com/help/matlab/ref/polyval.html> anyway, why not just use matlab’s syntax (polyval([c0, c1, c2...], x)) across the board?

4 18

Revised NEP-18, __array_function__ protocol
by Stephan Hoyer 01 Jul '18

01 Jul '18

After much discussion (and the addition of three new co-authors!), I’m pleased to present a significantly revision of NumPy Enhancement Proposal 18: A dispatch mechanism for NumPy's high level array functions: http://www.numpy.org/neps/nep-0018-array-function-protocol.html The full text is also included below. Best, Stephan =========================================================== A dispatch mechanism for NumPy's high level array functions =========================================================== :Author: Stephan Hoyer <shoyer(a)google.com> :Author: Matthew Rocklin <mrocklin(a)gmail.com> :Author: Marten van Kerkwijk <mhvk(a)astro.utoronto.ca> :Author: Hameer Abbasi <hameerabbasi(a)yahoo.com> :Author: Eric Wieser <wieser.eric(a)gmail.com> :Status: Draft :Type: Standards Track :Created: 2018-05-29 Abstact ------- We propose the ``__array_function__`` protocol, to allow arguments of NumPy functions to define how that function operates on them. This will allow using NumPy as a high level API for efficient multi-dimensional array operations, even with array implementations that differ greatly from ``numpy.ndarray``. Detailed description -------------------- NumPy's high level ndarray API has been implemented several times outside of NumPy itself for different architectures, such as for GPU arrays (CuPy), Sparse arrays (scipy.sparse, pydata/sparse) and parallel arrays (Dask array) as well as various NumPy-like implementations in the deep learning frameworks, like TensorFlow and PyTorch. Similarly there are many projects that build on top of the NumPy API for labeled and indexed arrays (XArray), automatic differentiation (Autograd, Tangent), masked arrays (numpy.ma), physical units (astropy.units, pint, unyt), etc. that add additional functionality on top of the NumPy API. Most of these project also implement a close variation of NumPy's level high API. We would like to be able to use these libraries together, for example we would like to be able to place a CuPy array within XArray, or perform automatic differentiation on Dask array code. This would be easier to accomplish if code written for NumPy ndarrays could also be used by other NumPy-like projects. For example, we would like for the following code example to work equally well with any NumPy-like array object: .. code:: python def f(x): y = np.tensordot(x, x.T) return np.mean(np.exp(y)) Some of this is possible today with various protocol mechanisms within NumPy. - The ``np.exp`` function checks the ``__array_ufunc__`` protocol - The ``.T`` method works using Python's method dispatch - The ``np.mean`` function explicitly checks for a ``.mean`` method on the argument However other functions, like ``np.tensordot`` do not dispatch, and instead are likely to coerce to a NumPy array (using the ``__array__``) protocol, or err outright. To achieve enough coverage of the NumPy API to support downstream projects like XArray and autograd we want to support *almost all* functions within NumPy, which calls for a more reaching protocol than just ``__array_ufunc__``. We would like a protocol that allows arguments of a NumPy function to take control and divert execution to another function (for example a GPU or parallel implementation) in a way that is safe and consistent across projects. Implementation -------------- We propose adding support for a new protocol in NumPy, ``__array_function__``. This protocol is intended to be a catch-all for NumPy functionality that is not covered by the ``__array_ufunc__`` protocol for universal functions (like ``np.exp``). The semantics are very similar to ``__array_ufunc__``, except the operation is specified by an arbitrary callable object rather than a ufunc instance and method. A prototype implementation can be found in `this notebook < https://nbviewer.jupyter.org/gist/shoyer/1f0a308a06cd96df20879a1ddb8f0006 >`_. The interface ~~~~~~~~~~~~~ We propose the following signature for implementations of ``__array_function__``: .. code-block:: python def __array_function__(self, func, types, args, kwargs) - ``func`` is an arbitrary callable exposed by NumPy's public API, which was called in the form ``func(*args, **kwargs)``. - ``types`` is a ``frozenset`` of unique argument types from the original NumPy function call that implement ``__array_function__``. - The tuple ``args`` and dict ``kwargs`` are directly passed on from the original call. Unlike ``__array_ufunc__``, there are no high-level guarantees about the type of ``func``, or about which of ``args`` and ``kwargs`` may contain objects implementing the array API. As a convenience for ``__array_function__`` implementors, ``types`` provides all argument types with an ``'__array_function__'`` attribute. This allows downstream implementations to quickly determine if they are likely able to support the operation. A ``frozenset`` is used to ensure that ``__array_function__`` implementations cannot rely on the iteration order of ``types``, which would facilitate violating the well-defined "Type casting hierarchy" described in `NEP-13 <https://www.numpy.org/neps/nep-0013-ufunc-overrides.html>`_. Example for a project implementing the NumPy API ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Most implementations of ``__array_function__`` will start with two checks: 1. Is the given function something that we know how to overload? 2. Are all arguments of a type that we know how to handle? If these conditions hold, ``__array_function__`` should return the result from calling its implementation for ``func(*args, **kwargs)``. Otherwise, it should return the sentinel value ``NotImplemented``, indicating that the function is not implemented by these types. This is preferable to raising ``TypeError`` directly, because it gives *other* arguments the opportunity to define the operations. There are no general requirements on the return value from ``__array_function__``, although most sensible implementations should probably return array(s) with the same type as one of the function's arguments. If/when Python gains `typing support for protocols <https://www.python.org/dev/peps/pep-0544/>`_ and NumPy adds static type annotations, the ``@overload`` implementation for ``SupportsArrayFunction`` will indicate a return type of ``Any``. It may also be convenient to define a custom decorators (``implements`` below) for registering ``__array_function__`` implementations. .. code:: python HANDLED_FUNCTIONS = {} class MyArray: def __array_function__(self, func, types, args, kwargs): if func not in HANDLED_FUNCTIONS: return NotImplemented # Note: this allows subclasses that don't override # __array_function__ to handle MyArray objects if not all(issubclass(t, MyArray) for t in types): return NotImplemented return HANDLED_FUNCTIONS[func](*args, **kwargs) def implements(numpy_function): """Register an __array_function__ implementation for MyArray objects.""" def decorator(func): HANDLED_FUNCTIONS[numpy_function] = func return func return decorator @implements(np.concatenate) def concatenate(arrays, axis=0, out=None): ... # implementation of concatenate for MyArray objects @implements(np.broadcast_to) def broadcast_to(array, shape): ... # implementation of broadcast_to for MyArray objects Note that it is not required for ``__array_function__`` implementations to include *all* of the corresponding NumPy function's optional arguments (e.g., ``broadcast_to`` above omits the irrelevant ``subok`` argument). Optional arguments are only passed in to ``__array_function__`` if they were explicitly used in the NumPy function call. Necessary changes within the NumPy codebase itself ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ This will require two changes within the NumPy codebase: 1. A function to inspect available inputs, look for the ``__array_function__`` attribute on those inputs, and call those methods appropriately until one succeeds. This needs to be fast in the common all-NumPy case, and have acceptable performance (no worse than linear time) even if the number of overloaded inputs is large (e.g., as might be the case for `np.concatenate`). This is one additional function of moderate complexity. 2. Calling this function within all relevant NumPy functions. This affects many parts of the NumPy codebase, although with very low complexity. Finding and calling the right ``__array_function__`` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Given a NumPy function, ``*args`` and ``**kwargs`` inputs, we need to search through ``*args`` and ``**kwargs`` for all appropriate inputs that might have the ``__array_function__`` attribute. Then we need to select among those possible methods and execute the right one. Negotiating between several possible implementations can be complex. Finding arguments ''''''''''''''''' Valid arguments may be directly in the ``*args`` and ``**kwargs``, such as in the case for ``np.tensordot(left, right, out=out)``, or they may be nested within lists or dictionaries, such as in the case of ``np.concatenate([x, y, z])``. This can be problematic for two reasons: 1. Some functions are given long lists of values, and traversing them might be prohibitively expensive. 2. Some functions may have arguments that we don't want to inspect, even if they have the ``__array_function__`` method. To resolve these issues, NumPy functions should explicitly indicate which of their arguments may be overloaded, and how these arguments should be checked. As a rule, this should include all arguments documented as either ``array_like`` or ``ndarray``. We propose to do so by writing "dispatcher" functions for each overloaded NumPy function: - These functions will be called with the exact same arguments that were passed into the NumPy function (i.e., ``dispatcher(*args, **kwargs)``), and should return an iterable of arguments to check for overrides. - Dispatcher functions are required to share the exact same positional, optional and keyword-only arguments as their corresponding NumPy functions. Otherwise, valid invocations of a NumPy function could result in an error when calling its dispatcher. - Because default *values* for keyword arguments do not have ``__array_function__`` attributes, by convention we set all default argument values to ``None``. This reduces the likelihood of signatures falling out of sync, and minimizes extraneous information in the dispatcher. The only exception should be cases where the argument value in some way effects dispatching, which should be rare. An example of the dispatcher for ``np.concatenate`` may be instructive: .. code:: python def _concatenate_dispatcher(arrays, axis=None, out=None): for array in arrays: yield array if out is not None: yield out The concatenate dispatcher is written as generator function, which allows it to potentially include the value of the optional ``out`` argument without needing to create a new sequence with the (potentially long) list of objects to be concatenated. Trying ``__array_function__`` methods until the right one works ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' Many arguments may implement the ``__array_function__`` protocol. Some of these may decide that, given the available inputs, they are unable to determine the correct result. How do we call the right one? If several are valid then which has precedence? For the most part, the rules for dispatch with ``__array_function__`` match those for ``__array_ufunc__`` (see `NEP-13 <https://www.numpy.org/neps/nep-0013-ufunc-overrides.html>`_). In particular: - NumPy will gather implementations of ``__array_function__`` from all specified inputs and call them in order: subclasses before superclasses, and otherwise left to right. Note that in some edge cases involving subclasses, this differs slightly from the `current behavior <https://bugs.python.org/issue30140>`_ of Python. - Implementations of ``__array_function__`` indicate that they can handle the operation by returning any value other than ``NotImplemented``. - If all ``__array_function__`` methods return ``NotImplemented``, NumPy will raise ``TypeError``. One deviation from the current behavior of ``__array_ufunc__`` is that NumPy will only call ``__array_function__`` on the *first* argument of each unique type. This matches Python's `rule for calling reflected methods < https://docs.python.org/3/reference/datamodel.html#object.__ror__>`_, and this ensures that checking overloads has acceptable performance even when there are a large number of overloaded arguments. To avoid long-term divergence between these two dispatch protocols, we should `also update <https://github.com/numpy/numpy/issues/11306>`_ ``__array_ufunc__`` to match this behavior. Special handling of ``numpy.ndarray`` ''''''''''''''''''''''''''''''''''''' The use cases for subclasses with ``__array_function__`` are the same as those with ``__array_ufunc__``, so ``numpy.ndarray`` should also define a ``__array_function__`` method mirroring ``ndarray.__array_ufunc__``: .. code:: python def __array_function__(self, func, types, args, kwargs): # Cannot handle items that have __array_function__ other than our own. for t in types: if (hasattr(t, '__array_function__') and t.__array_function__ is not ndarray.__array_function__): return NotImplemented # Arguments contain no overrides, so we can safely call the # overloaded function again. return func(*args, **kwargs) To avoid infinite recursion, the dispatch rules for ``__array_function__`` need also the same special case they have for ``__array_ufunc__``: any arguments with an ``__array_function__`` method that is identical to ``numpy.ndarray.__array_function__`` are not be called as ``__array_function__`` implementations. Changes within NumPy functions ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Given a function defining the above behavior, for now call it ``try_array_function_override``, we now need to call that function from within every relevant NumPy function. This is a pervasive change, but of fairly simple and innocuous code that should complete quickly and without effect if no arguments implement the ``__array_function__`` protocol. In most cases, these functions should written using the ``array_function_dispatch`` decorator, which also associates dispatcher functions: .. code:: python def array_function_dispatch(dispatcher): """Wrap a function for dispatch with the __array_function__ protocol.""" def decorator(func): @functools.wraps(func) def new_func(*args, **kwargs): relevant_arguments = dispatcher(*args, **kwargs) success, value = try_array_function_override( new_func, relevant_arguments, args, kwargs) if success: return value return func(*args, **kwargs) return new_func return decorator # example usage def _broadcast_to_dispatcher(array, shape, subok=None, **ignored_kwargs): return (array,) @array_function_dispatch(_broadcast_to_dispatcher) def broadcast_to(array, shape, subok=False): ... # existing definition of np.broadcast_to Using a decorator is great! We don't need to change the definitions of existing NumPy functions, and only need to write a few additional lines for the dispatcher function. We could even reuse a single dispatcher for families of functions with the same signature (e.g., ``sum`` and ``prod``). For such functions, the largest change could be adding a few lines to the docstring to note which arguments are checked for overloads. It's particularly worth calling out the decorator's use of ``functools.wraps``: - This ensures that the wrapped function has the same name and docstring as the wrapped NumPy function. - On Python 3, it also ensures that the decorator function copies the original function signature, which is important for introspection based tools such as auto-complete. If we care about preserving function signatures on Python 2, for the `short while longer < http://www.numpy.org/neps/nep-0014-dropping-python2.7-proposal.html>`_ that NumPy supports Python 2.7, we do could do so by adding a vendored dependency on the (single-file, BSD licensed) `decorator library <https://github.com/micheles/decorator>`_. - Finally, it ensures that the wrapped function `can be pickled < http://gael-varoquaux.info/programming/decoration-in-python-done-right-deco… >`_. In a few cases, it would not make sense to use the ``array_function_dispatch`` decorator directly, but override implementation in terms of ``try_array_function_override`` should still be straightforward. - Functions written entirely in C (e.g., ``np.concatenate``) can't use decorators, but they could still use a C equivalent of ``try_array_function_override``. If performance is not a concern, they could also be easily wrapped with a small Python wrapper. - The ``__call__`` method of ``np.vectorize`` can't be decorated with <p style="margin:0px;font-stretch:normal;font-size:17.4px;line-

6 20