diff options
Diffstat (limited to 'doc/neps')
28 files changed, 2444 insertions, 374 deletions
diff --git a/doc/neps/conf.py b/doc/neps/conf.py index aa11d37b3..8cfb2b570 100644 --- a/doc/neps/conf.py +++ b/doc/neps/conf.py @@ -100,7 +100,7 @@ todo_include_todos = False ## to template names. ## ## This is required for the alabaster theme -## refs: http://alabaster.readthedocs.io/en/latest/installation.html#sidebars +## refs: https://alabaster.readthedocs.io/en/latest/installation.html#sidebars #html_sidebars = { # '**': [ # 'relations.html', # needs 'show_related': True theme option to display @@ -127,8 +127,8 @@ if True: "edit_link": True, "sidebar": "right", "scipy_org_logo": True, - "rootlinks": [("http://scipy.org/", "Scipy.org"), - ("http://docs.scipy.org/", "Docs")] + "rootlinks": [("https://scipy.org/", "Scipy.org"), + ("https://docs.scipy.org/", "Docs")] } else: # Default build diff --git a/doc/neps/index.rst.tmpl b/doc/neps/index.rst.tmpl index 6c988014f..e7b8fedba 100644 --- a/doc/neps/index.rst.tmpl +++ b/doc/neps/index.rst.tmpl @@ -1,12 +1,21 @@ -=========================== -NumPy Enhancement Proposals -=========================== - -NumPy Enhancement Proposals (NEPs) describe proposed changes to NumPy. -NEPs are modeled on Python Enhancement Proposals (PEPs), and are typically -written up when large changes to NumPy are proposed. +===================================== +Roadmap & NumPy Enhancement Proposals +===================================== + +This page provides an overview of development priorities for NumPy. +Specifically, it contains a roadmap with a higher-level overview, as +well as NumPy Enhancement Proposals (NEPs)—suggested changes +to the library—in various stages of discussion or completion (see `NEP +0 <nep-0000>`__). + +Roadmap +------- +.. toctree:: + :maxdepth: 1 -This page provides an overview of all NEPs. + The Scope of NumPy <scope> + Current roadmap <roadmap> + Wish list <https://github.com/numpy/numpy/issues?q=is%3Aopen+is%3Aissue+label%3A%2223+-+Wish+List%22> Meta-NEPs (NEPs about NEPs or Processes) ---------------------------------------- @@ -15,7 +24,7 @@ Meta-NEPs (NEPs about NEPs or Processes) :maxdepth: 1 {% for nep, tags in neps.items() if tags['Type'] == 'Process' %} - NEP {{ nep }} — {{ tags['Title'] }} <{{ tags['Filename'] }}> + {{ tags['Title'] }} <{{ tags['Filename'] }}> {% endfor %} nep-template @@ -27,7 +36,7 @@ Accepted NEPs, implementation in progress :maxdepth: 1 {% for nep, tags in neps.items() if tags['Status'] == 'Accepted' %} - NEP {{ nep }} — {{ tags['Title'] }} <{{ tags['Filename'] }}> + {{ tags['Title'] }} <{{ tags['Filename'] }}> {% endfor %} @@ -38,7 +47,7 @@ Open NEPs (under consideration) :maxdepth: 1 {% for nep, tags in neps.items() if tags['Status'] == 'Draft' %} - NEP {{ nep }} — {{ tags['Title'] }} <{{ tags['Filename'] }}> + {{ tags['Title'] }} <{{ tags['Filename'] }}> {% endfor %} @@ -50,7 +59,7 @@ Implemented NEPs :maxdepth: 1 {% for nep, tags in neps.items() if tags['Status'] == 'Final' %} - NEP {{ nep }} — {{ tags['Title'] }} <{{ tags['Filename'] }}> + {{ tags['Title'] }} <{{ tags['Filename'] }}> {% endfor %} Deferred NEPs @@ -60,7 +69,7 @@ Deferred NEPs :maxdepth: 1 {% for nep, tags in neps.items() if tags['Status'] == 'Deferred' %} - NEP {{ nep }} — {{ tags['Title'] }} <{{ tags['Filename'] }}> + {{ tags['Title'] }} <{{ tags['Filename'] }}> {% endfor %} Rejected NEPs @@ -70,5 +79,5 @@ Rejected NEPs :maxdepth: 1 {% for nep, tags in neps.items() if tags['Status'] == 'Rejected' %} - NEP {{ nep }} — {{ tags['Title'] }} <{{ tags['Filename'] }}> + {{ tags['Title'] }} <{{ tags['Filename'] }}> {% endfor %} diff --git a/doc/neps/nep-0000.rst b/doc/neps/nep-0000.rst index 9c6646db2..a3ec3a42b 100644 --- a/doc/neps/nep-0000.rst +++ b/doc/neps/nep-0000.rst @@ -1,6 +1,6 @@ -=================== -Purpose and Process -=================== +=========================== +NEP 0 — Purpose and Process +=========================== :Author: Jarrod Millman <millman@berkeley.edu> :Status: Active @@ -97,16 +97,9 @@ status of NEPs are as follows: All NEPs should be created with the ``Draft`` status. -Normally, a NEP is ``Accepted`` by consensus of all interested -Contributors. To verify that consensus has been reached, the NEP -author or another interested party should make a post on the -numpy-discussion mailing list proposing it for acceptance; if there -are no substantive objections after one week, the NEP can officially -be marked ``Accepted``, and a link to this post should be added to the -NEP for reference. - -In unusual cases, the `NumPy Steering Council`_ may be asked to decide whether -a controversial NEP is ``Accepted``. +Eventually, after discussion, there may be a consensus that the NEP +should be accepted – see the next section for details. At this point +the status becomes ``Accepted``. Once a NEP has been ``Accepted``, the reference implementation must be completed. When the reference implementation is complete and incorporated @@ -135,6 +128,61 @@ Process NEPs may also have a status of ``Active`` if they are never meant to be completed, e.g. NEP 0 (this NEP). +How a NEP becomes Accepted +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +A NEP is ``Accepted`` by consensus of all interested contributors. We +need a concrete way to tell whether consensus has been reached. When +you think a NEP is ready to accept, send an email to the +numpy-discussion mailing list with a subject like: + + Proposal to accept NEP #<number>: <title> + +In the body of your email, you should: + +* link to the latest version of the NEP, + +* briefly describe any major points of contention and how they were + resolved, + +* include a sentence like: "If there are no substantive objections + within 7 days from this email, then the NEP will be accepted; see + NEP 0 for more details." + +For an example, see: https://mail.python.org/pipermail/numpy-discussion/2018-June/078345.html + +After you send the email, you should make sure to link to the email +thread from the ``Discussion`` section of the NEP, so that people can +find it later. + +Generally the NEP author will be the one to send this email, but +anyone can do it – the important thing is to make sure that everyone +knows when a NEP is on the verge of acceptance, and give them a final +chance to respond. If there's some special reason to extend this final +comment period beyond 7 days, then that's fine, just say so in the +email. You shouldn't do less than 7 days, because sometimes people are +travelling or similar and need some time to respond. + +In general, the goal is to make sure that the community has consensus, +not provide a rigid policy for people to try to game. When in doubt, +err on the side of asking for more feedback and looking for +opportunities to compromise. + +If the final comment period passes without any substantive objections, +then the NEP can officially be marked ``Accepted``. You should send a +followup email notifying the list (celebratory emoji optional but +encouraged 🎉✨), and then update the NEP by setting its ``:Status:`` +to ``Accepted``, and its ``:Resolution:`` header to a link to your +followup email. + +If there *are* substantive objections, then the NEP remains in +``Draft`` state, discussion continues as normal, and it can be +proposed for acceptance again later once the objections are resolved. + +In unusual cases, the `NumPy Steering Council`_ may be asked to decide +whether a controversial NEP is ``Accepted``. + + Maintenance ^^^^^^^^^^^ @@ -203,7 +251,7 @@ References and Footnotes `GitHub <https://github.com/numpy/numpy/tree/master/doc/neps>`_. .. [2] The URL for viewing NEPs on the web is - http://numpy.github.io/neps/. + https://www.numpy.org/neps/. .. _repo: https://github.com/numpy/numpy @@ -220,7 +268,7 @@ References and Footnotes .. _reStructuredTextPrimer: http://www.sphinx-doc.org/en/stable/rest.html -.. _Sphinx: www.sphinx-doc.org/en/stable +.. _Sphinx: http://www.sphinx-doc.org/en/stable/ Copyright diff --git a/doc/neps/nep-0001-npy-format.rst b/doc/neps/nep-0001-npy-format.rst index 2057aed83..4eded02ff 100644 --- a/doc/neps/nep-0001-npy-format.rst +++ b/doc/neps/nep-0001-npy-format.rst @@ -1,6 +1,6 @@ -===================================== -A Simple File Format for NumPy Arrays -===================================== +============================================= +NEP 1 — A Simple File Format for NumPy Arrays +============================================= :Author: Robert Kern <robert.kern@gmail.com> :Status: Final @@ -290,15 +290,15 @@ included in the 1.9.0 release of numpy. Specifically, the file format.py in this directory implements the format as described here. - http://github.com/numpy/numpy/blob/master/numpy/lib/format.py + https://github.com/numpy/numpy/blob/master/numpy/lib/format.py References ---------- -[1] http://docs.python.org/lib/module-pickle.html +[1] https://docs.python.org/library/pickle.html -[2] http://hdf.ncsa.uiuc.edu/products/hdf5/index.html +[2] https://support.hdfgroup.org/HDF5/ Copyright diff --git a/doc/neps/nep-0002-warnfix.rst b/doc/neps/nep-0002-warnfix.rst index 60dc885b2..207dfa3d4 100644 --- a/doc/neps/nep-0002-warnfix.rst +++ b/doc/neps/nep-0002-warnfix.rst @@ -1,6 +1,6 @@ -========================================================================= -A proposal to build numpy without warning with a big set of warning flags -========================================================================= +================================================================================= +NEP 2 — A proposal to build numpy without warning with a big set of warning flags +================================================================================= :Author: David Cournapeau :Contact: david@ar.media.kyoto-u.ac.jp diff --git a/doc/neps/nep-0003-math_config_clean.rst b/doc/neps/nep-0003-math_config_clean.rst index 5af907437..ebd32b124 100644 --- a/doc/neps/nep-0003-math_config_clean.rst +++ b/doc/neps/nep-0003-math_config_clean.rst @@ -1,6 +1,6 @@ -=========================================================== -Cleaning the math configuration of numpy.core -=========================================================== +===================================================== +NEP 3 — Cleaning the math configuration of numpy.core +===================================================== :Author: David Cournapeau :Contact: david@ar.media.kyoto-u.ac.jp diff --git a/doc/neps/nep-0004-datetime-proposal3.rst b/doc/neps/nep-0004-datetime-proposal3.rst index 46d8e314b..b32964e88 100644 --- a/doc/neps/nep-0004-datetime-proposal3.rst +++ b/doc/neps/nep-0004-datetime-proposal3.rst @@ -1,6 +1,6 @@ -==================================================================== - A (third) proposal for implementing some date/time types in NumPy -==================================================================== +========================================================================= +NEP 4 — A (third) proposal for implementing some date/time types in NumPy +========================================================================= :Author: Francesc Alted i Abad :Contact: faltet@pytables.com @@ -562,9 +562,9 @@ examples of other derived units, and we find this a bit too overwhelming for this proposal purposes. -.. [1] http://docs.python.org/lib/module-datetime.html -.. [2] http://www.egenix.com/products/python/mxBase/mxDateTime -.. [3] http://en.wikipedia.org/wiki/Unix_time +.. [1] https://docs.python.org/library/datetime.html +.. [2] https://www.egenix.com/products/python/mxBase/mxDateTime +.. [3] https://en.wikipedia.org/wiki/Unix_time .. Local Variables: diff --git a/doc/neps/nep-0005-generalized-ufuncs.rst b/doc/neps/nep-0005-generalized-ufuncs.rst index 54b2b370e..366e26ffd 100644 --- a/doc/neps/nep-0005-generalized-ufuncs.rst +++ b/doc/neps/nep-0005-generalized-ufuncs.rst @@ -1,6 +1,6 @@ -=============================== -Generalized Universal Functions -=============================== +======================================= +NEP 5 — Generalized Universal Functions +======================================= :Status: Final diff --git a/doc/neps/nep-0006-newbugtracker.rst b/doc/neps/nep-0006-newbugtracker.rst index 2b9344ed0..8dc7a1d8e 100644 --- a/doc/neps/nep-0006-newbugtracker.rst +++ b/doc/neps/nep-0006-newbugtracker.rst @@ -1,6 +1,6 @@ -=========================================== -Replacing Trac with a different bug tracker -=========================================== +=================================================== +NEP 6 — Replacing Trac with a different bug tracker +=================================================== :Author: David Cournapeau, Stefan van der Walt :Status: Deferred diff --git a/doc/neps/nep-0007-datetime-proposal.rst b/doc/neps/nep-0007-datetime-proposal.rst index 72d48d244..5547a4306 100644 --- a/doc/neps/nep-0007-datetime-proposal.rst +++ b/doc/neps/nep-0007-datetime-proposal.rst @@ -1,6 +1,6 @@ -==================================================================== - A proposal for implementing some date/time types in NumPy -==================================================================== +================================================================== +NEP 7 — A proposal for implementing some date/time types in NumPy +================================================================== :Author: Travis Oliphant :Contact: oliphant@enthought.com @@ -662,9 +662,9 @@ operations mixing business days with other time units will not be allowed. -.. [1] http://docs.python.org/lib/module-datetime.html -.. [2] http://www.egenix.com/products/python/mxBase/mxDateTime -.. [3] http://en.wikipedia.org/wiki/Unix_time +.. [1] https://docs.python.org/library/datetime.html +.. [2] https://www.egenix.com/products/python/mxBase/mxDateTime +.. [3] https://en.wikipedia.org/wiki/Unix_time .. Local Variables: diff --git a/doc/neps/nep-0008-groupby_additions.rst b/doc/neps/nep-0008-groupby_additions.rst index fa02f2f9c..3189fcf41 100644 --- a/doc/neps/nep-0008-groupby_additions.rst +++ b/doc/neps/nep-0008-groupby_additions.rst @@ -1,6 +1,6 @@ -==================================================================== - A proposal for adding groupby functionality to NumPy -==================================================================== +============================================================= +NEP 8 — A proposal for adding groupby functionality to NumPy +============================================================= :Author: Travis Oliphant :Contact: oliphant@enthought.com diff --git a/doc/neps/nep-0009-structured_array_extensions.rst b/doc/neps/nep-0009-structured_array_extensions.rst index 695d0d516..8b81a308d 100644 --- a/doc/neps/nep-0009-structured_array_extensions.rst +++ b/doc/neps/nep-0009-structured_array_extensions.rst @@ -1,6 +1,6 @@ -=========================== -Structured array extensions -=========================== +=================================== +NEP 9 — Structured array extensions +=================================== :Status: Deferred diff --git a/doc/neps/nep-0010-new-iterator-ufunc.rst b/doc/neps/nep-0010-new-iterator-ufunc.rst index 7b388a974..8601b4a4c 100644 --- a/doc/neps/nep-0010-new-iterator-ufunc.rst +++ b/doc/neps/nep-0010-new-iterator-ufunc.rst @@ -1,6 +1,6 @@ -===================================== -Optimizing Iterator/UFunc Performance -===================================== +============================================== +NEP 10 — Optimizing Iterator/UFunc Performance +============================================== :Author: Mark Wiebe <mwwiebe@gmail.com> :Content-Type: text/x-rst diff --git a/doc/neps/nep-0011-deferred-ufunc-evaluation.rst b/doc/neps/nep-0011-deferred-ufunc-evaluation.rst index 5f5de3518..a7143c6ee 100644 --- a/doc/neps/nep-0011-deferred-ufunc-evaluation.rst +++ b/doc/neps/nep-0011-deferred-ufunc-evaluation.rst @@ -1,6 +1,6 @@ -========================= -Deferred UFunc Evaluation -========================= +================================== +NEP 11 — Deferred UFunc Evaluation +================================== :Author: Mark Wiebe <mwwiebe@gmail.com> :Content-Type: text/x-rst diff --git a/doc/neps/nep-0012-missing-data.rst b/doc/neps/nep-0012-missing-data.rst index 1553339f4..dbcf1b579 100644 --- a/doc/neps/nep-0012-missing-data.rst +++ b/doc/neps/nep-0012-missing-data.rst @@ -1,10 +1,10 @@ -=================================== -Missing Data Functionality in NumPy -=================================== +============================================ +NEP 12 — Missing Data Functionality in NumPy +============================================ :Author: Mark Wiebe <mwwiebe@gmail.com> :Copyright: Copyright 2011 by Enthought, Inc -:License: CC By-SA 3.0 (http://creativecommons.org/licenses/by-sa/3.0/) +:License: CC By-SA 3.0 (https://creativecommons.org/licenses/by-sa/3.0/) :Date: 2011-06-23 :Status: Deferred @@ -224,7 +224,7 @@ but with semantics reflecting its status as a missing value. In particular, trying to treat it as a boolean will raise an exception, and comparisons with it will produce numpy.NA instead of True or False. These basics are adopted from the behavior of the NA value in the R project. To dig -deeper into the ideas, http://en.wikipedia.org/wiki/Ternary_logic#Kleene_logic +deeper into the ideas, https://en.wikipedia.org/wiki/Ternary_logic#Kleene_logic provides a starting point. For example,:: @@ -857,7 +857,7 @@ Shared Masks One feature of numpy.ma is called 'shared masks'. -http://docs.scipy.org/doc/numpy/reference/maskedarray.baseclass.html#numpy.ma.MaskedArray.sharedmask +https://docs.scipy.org/doc/numpy/reference/maskedarray.baseclass.html#numpy.ma.MaskedArray.sharedmask This feature cannot be supported by a masked implementation of missing values without directly violating the missing value abstraction. @@ -888,7 +888,7 @@ found from doing google searches of numpy C API array access. NumPy Documentation - How to extend NumPy ----------------------------------------- -http://docs.scipy.org/doc/numpy/user/c-info.how-to-extend.html#dealing-with-array-objects +https://docs.scipy.org/doc/numpy/user/c-info.how-to-extend.html#dealing-with-array-objects This page has a section "Dealing with array objects" which has some advice for how to access numpy arrays from C. When accepting arrays, the first step it suggests is @@ -898,7 +898,7 @@ advice will properly fail when given an NA-masked array it doesn't know how to h The way this is handled is that PyArray_FromAny requires a special flag, NPY_ARRAY_ALLOWNA, before it will allow NA-masked arrays to flow through. -http://docs.scipy.org/doc/numpy/reference/c-api.array.html#NPY_ARRAY_ALLOWNA +https://docs.scipy.org/doc/numpy/reference/c-api.array.html#NPY_ARRAY_ALLOWNA Code which does not follow this advice, and instead just calls PyArray_Check() to verify its an ndarray and checks some flags, will silently produce incorrect results. This style diff --git a/doc/neps/nep-0013-ufunc-overrides.rst b/doc/neps/nep-0013-ufunc-overrides.rst index c97b69023..a51ce3927 100644 --- a/doc/neps/nep-0013-ufunc-overrides.rst +++ b/doc/neps/nep-0013-ufunc-overrides.rst @@ -1,6 +1,6 @@ -================================= -A Mechanism for Overriding Ufuncs -================================= +========================================== +NEP 13 — A Mechanism for Overriding Ufuncs +========================================== .. currentmodule:: numpy @@ -53,7 +53,7 @@ changes in 3rd party code. .. [1] http://docs.python.org/doc/numpy/user/basics.subclassing.html .. [2] https://github.com/scipy/scipy/issues/2123 .. [3] https://github.com/scipy/scipy/issues/1569 -.. [4] http://technicaldiscovery.blogspot.com/2013/07/thoughts-after-scipy-2013-and-specific.html +.. [4] https://technicaldiscovery.blogspot.com/2013/07/thoughts-after-scipy-2013-and-specific.html Motivation @@ -134,7 +134,7 @@ which have multiplication semantics incompatible with numpy arrays. However, the aim is to enable writing other custom array types that have strictly ndarray compatible semantics. -.. [5] http://mail.python.org/pipermail/numpy-discussion/2011-June/056945.html +.. [5] https://mail.python.org/pipermail/numpy-discussion/2011-June/056945.html .. [6] https://github.com/numpy/numpy/issues/5844 @@ -635,7 +635,7 @@ simplify the dispatch logic for binary operations with NumPy arrays as much as possible, by making it possible to use Python's dispatch rules or NumPy's dispatch rules, but not some mixture of both at the same time. -.. [9] http://bugs.python.org/issue30140 +.. [9] https://bugs.python.org/issue30140 .. _neps.ufunc-overrides.list-of-operators: diff --git a/doc/neps/nep-0014-dropping-python2.7-proposal.rst b/doc/neps/nep-0014-dropping-python2.7-proposal.rst index 6cfd4707f..3adf3b407 100644 --- a/doc/neps/nep-0014-dropping-python2.7-proposal.rst +++ b/doc/neps/nep-0014-dropping-python2.7-proposal.rst @@ -1,6 +1,6 @@ -==================================== -Plan for dropping Python 2.7 support -==================================== +============================================= +NEP 14 — Plan for dropping Python 2.7 support +============================================= :Status: Accepted :Resolution: https://mail.python.org/pipermail/numpy-discussion/2017-November/077419.html @@ -50,6 +50,6 @@ to Python3 only, see the python3-statement_. For more information on porting your code to run on Python 3, see the python3-howto_. -.. _python3-statement: http://www.python3statement.org/ +.. _python3-statement: https://python3statement.org/ .. _python3-howto: https://docs.python.org/3/howto/pyporting.html diff --git a/doc/neps/nep-0015-merge-multiarray-umath.rst b/doc/neps/nep-0015-merge-multiarray-umath.rst new file mode 100644 index 000000000..7c1f5faf8 --- /dev/null +++ b/doc/neps/nep-0015-merge-multiarray-umath.rst @@ -0,0 +1,157 @@ +===================================== +NEP 15 — Merging multiarray and umath +===================================== + +:Author: Nathaniel J. Smith <njs@pobox.com> +:Status: Accepted +:Type: Standards Track +:Created: 2018-02-22 +:Resolution: https://mail.python.org/pipermail/numpy-discussion/2018-June/078345.html + +Abstract +-------- + +Let's merge ``numpy.core.multiarray`` and ``numpy.core.umath`` into a +single extension module, and deprecate ``np.set_numeric_ops``. + + +Background +---------- + +Currently, numpy's core C code is split between two separate extension +modules. + +``numpy.core.multiarray`` is built from +``numpy/core/src/multiarray/*.c``, and contains the core array +functionality (in particular, the ``ndarray`` object). + +``numpy.core.umath`` is built from ``numpy/core/src/umath/*.c``, and +contains the ufunc machinery. + +These two modules each expose their own separate C API, accessed via +``import_multiarray()`` and ``import_umath()`` respectively. The idea +is that they're supposed to be independent modules, with +``multiarray`` as a lower-level layer with ``umath`` built on top. In +practice this has turned out to be problematic. + +First, the layering isn't perfect: when you write ``ndarray + +ndarray``, this invokes ``ndarray.__add__``, which then calls the +ufunc ``np.add``. This means that ``ndarray`` needs to know about +ufuncs – so instead of a clean layering, we have a circular +dependency. To solve this, ``multiarray`` exports a somewhat +terrifying function called ``set_numeric_ops``. The bootstrap +procedure each time you ``import numpy`` is: + +1. ``multiarray`` and its ``ndarray`` object are loaded, but + arithmetic operations on ndarrays are broken. + +2. ``umath`` is loaded. + +3. ``set_numeric_ops`` is used to monkeypatch all the methods like + ``ndarray.__add__`` with objects from ``umath``. + +In addition, ``set_numeric_ops`` is exposed as a public API, +``np.set_numeric_ops``. + +Furthermore, even when this layering does work, it ends up distorting +the shape of our public ABI. In recent years, the most common reason +for adding new functions to ``multiarray``\'s "public" ABI is not that +they really need to be public or that we expect other projects to use +them, but rather just that we need to call them from ``umath``. This +is extremely unfortunate, because it makes our public ABI +unnecessarily large, and since we can never remove things from it then +this creates an ongoing maintenance burden. The way C works, you can +have internal API that's visible to everything inside the same +extension module, or you can have a public API that everyone can use; +you can't (easily) have an API that's visible to multiple extension +modules inside numpy, but not to external users. + +We've also increasingly been putting utility code into +``numpy/core/src/private/``, which now contains a bunch of files which +are ``#include``\d twice, once into ``multiarray`` and once into +``umath``. This is pretty gross, and is purely a workaround for these +being separate C extensions. The ``npymath`` library is also +included in both extension modules. + + +Proposed changes +---------------- + +This NEP proposes three changes: + +1. We should start building ``numpy/core/src/multiarray/*.c`` and + ``numpy/core/src/umath/*.c`` together into a single extension + module. + +2. Instead of ``set_numeric_ops``, we should use some new, private API + to set up ``ndarray.__add__`` and friends. + +3. We should deprecate, and eventually remove, ``np.set_numeric_ops``. + + +Non-proposed changes +-------------------- + +We don't necessarily propose to throw away the distinction between +multiarray/ and umath/ in terms of our source code organization: +internal organization is useful! We just want to build them together +into a single extension module. Of course, this does open the door for +potential future refactorings, which we can then evaluate based on +their merits as they come up. + +It also doesn't propose that we break the public C ABI. We should +continue to provide ``import_multiarray()`` and ``import_umath()`` +functions – it's just that now both ABIs will ultimately be loaded +from the same C library. Due to how ``import_multiarray()`` and +``import_umath()`` are written, we'll also still need to have modules +called ``numpy.core.multiarray`` and ``numpy.core.umath``, and they'll +need to continue to export ``_ARRAY_API`` and ``_UFUNC_API`` objects – +but we can make one or both of these modules be tiny shims that simply +re-export the magic API object from where-ever it's actually defined. +(See ``numpy/core/code_generators/generate_{numpy,ufunc}_api.py`` for +details of how these imports work.) + + +Backward compatibility +---------------------- + +The only compatibility break is the deprecation of ``np.set_numeric_ops``. + + +Rejected alternatives +--------------------- + +Preserve ``set_numeric_ops`` for monkeypatching +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +In discussing this NEP, one additional use case was raised for +``set_numeric_ops``: if you have an optimized vector math library +(e.g. Intel's MKL VML, Sleef, or Yeppp), then ``set_numeric_ops`` can +be used to monkeypatch numpy to use these operations instead of +numpy's built-in vector operations. But, even if we grant that this is +a great idea, using ``set_numeric_ops`` isn't actually the best way to +do it. All ``set_numeric_ops`` allows you to do is take over Python's +syntactic operators (``+``, ``*``, etc.) on ndarrays; it doesn't let +you affect operations called via other APIs (e.g., ``np.add``), or +operations that don't have built-in syntax (e.g., ``np.exp``). Also, +you have to reimplement the whole ufunc machinery, instead of just the +core loop. On the other hand, the `PyUFunc_ReplaceLoopBySignature +<https://docs.scipy.org/doc/numpy/reference/c-api.ufunc.html#c.PyUFunc_ReplaceLoopBySignature>`__ +API – which was added in 2006 – allows replacement of the inner loops +of arbitrary ufuncs. This is both simpler and more powerful – e.g. +replacing the inner loop of ``np.add`` means your code will +automatically be used for both ``ndarray + ndarray`` as well as direct +calls to ``np.add``. So this doesn't seem like a good reason to not +deprecate ``set_numeric_ops``. + + +Discussion +---------- + +* https://mail.python.org/pipermail/numpy-discussion/2018-March/077764.html +* https://mail.python.org/pipermail/numpy-discussion/2018-June/078345.html + +Copyright +--------- + +This document has been placed in the public domain. diff --git a/doc/neps/nep-0017-split-out-maskedarray.rst b/doc/neps/nep-0017-split-out-maskedarray.rst index d6dcc1def..7ef949763 100644 --- a/doc/neps/nep-0017-split-out-maskedarray.rst +++ b/doc/neps/nep-0017-split-out-maskedarray.rst @@ -1,6 +1,6 @@ -======================= -Split Out Masked Arrays -======================= +================================ +NEP 17 — Split Out Masked Arrays +================================ :Author: Stéfan van der Walt <stefanv@berkeley.edu> :Status: Rejected diff --git a/doc/neps/nep-0018-array-function-protocol.rst b/doc/neps/nep-0018-array-function-protocol.rst index 943ca4cbf..d4ba7879b 100644 --- a/doc/neps/nep-0018-array-function-protocol.rst +++ b/doc/neps/nep-0018-array-function-protocol.rst @@ -1,9 +1,12 @@ -================================================== -NEP: Dispatch Mechanism for NumPy's high level API -================================================== +==================================================================== +NEP 18 — A dispatch mechanism for NumPy's high level array functions +==================================================================== :Author: Stephan Hoyer <shoyer@google.com> :Author: Matthew Rocklin <mrocklin@gmail.com> +:Author: Marten van Kerkwijk <mhvk@astro.utoronto.ca> +:Author: Hameer Abbasi <hameerabbasi@yahoo.com> +:Author: Eric Wieser <wieser.eric@gmail.com> :Status: Draft :Type: Standards Track :Created: 2018-05-29 @@ -11,25 +14,27 @@ NEP: Dispatch Mechanism for NumPy's high level API Abstact ------- -We propose a protocol to allow arguments of numpy functions to define -how that function operates on them. This allows other libraries that -implement NumPy's high level API to reuse Numpy functions. This allows -libraries that extend NumPy's high level API to apply to more NumPy-like -libraries. +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 +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 +arrays (Dask array) as well as various NumPy-like implementations in the deep learning frameworks, like TensorFlow and PyTorch. -Similarly there are several projects that build on top of the Numpy API -for labeled and indexed arrays (XArray), automatic differentation -(Autograd, Tangent), higher order array factorizations (TensorLy), etc. -that add additional functionality on top of the Numpy API. +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 @@ -38,7 +43,7 @@ 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: +equally well with any NumPy-like array object: .. code:: python @@ -47,7 +52,7 @@ equally well with any Numpy-like array object: return np.mean(np.exp(y)) Some of this is possible today with various protocol mechanisms within -Numpy. +NumPy. - The ``np.exp`` function checks the ``__array_ufunc__`` protocol - The ``.T`` method works using Python's method dispatch @@ -55,10 +60,10 @@ Numpy. the argument However other functions, like ``np.tensordot`` do not dispatch, and -instead are likely to coerce to a Numpy array (using the ``__array__``) +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 +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 @@ -71,10 +76,13 @@ 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 existing protocols, like reductions (like ``np.sum``) -or 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. +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 ~~~~~~~~~~~~~ @@ -88,23 +96,23 @@ We propose the following signature for implementations of - ``func`` is an arbitrary callable exposed by NumPy's public API, which was called in the form ``func(*args, **kwargs)``. -- ``types`` is a list of types for all arguments to the original NumPy - function call that will be checked for an ``__array_function__`` - implementation. -- The tuple ``args`` and dict ``**kwargs`` are directly passed on from the +- ``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 of the NumPy API, the ``types`` keyword contains a list of all -types that implement the ``__array_function__`` protocol. This allows -downstream implementations to quickly determine if they are likely able to -support the operation. +implementing the array API. -Still be determined: what guarantees can we offer for ``types``? Should -we promise that types are unique, and appear in the order in which they -are checked? +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 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -118,45 +126,78 @@ checks: 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. +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) - HANDLED_FUNCTIONS = { - np.concatenate: my_concatenate, - np.broadcast_to: my_broadcast_to, - np.sum: my_sum, - ... - } + 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 -Necessary changes within the Numpy codebase itself +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: +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. + 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. +2. Calling this function within all relevant NumPy functions. - This affects many parts of the Numpy codebase, although with very low + 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 +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. @@ -171,13 +212,46 @@ 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 function may have arguments that we don't want to inspect, even - if they have the ``__array_function__`` method + 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 -To resolve these we ask the functions to provide an explicit list of -arguments that should be traversed. This is the ``relevant_arguments=`` -keyword in the examples below. +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 ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' @@ -187,15 +261,15 @@ 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? -The rules for dispatch with ``__array_function__`` match those for -``__array_ufunc__`` (see -`NEP-13 <http://www.numpy.org/neps/nep-0013-ufunc-overrides.html>`_). +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, - this differs slightly from the + 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 @@ -203,69 +277,194 @@ In particular: - If all ``__array_function__`` methods return ``NotImplemented``, NumPy will raise ``TypeError``. -Changes within Numpy functions +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 defined above, for now call it -``do_array_function_dance``, we now need to call that function from -within every relevant Numpy function. This is a pervasive change, but of +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. Let us consider a few examples of NumPy functions and how they -might be affected by this change: +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): - success, value = do_array_function_dance( - func=broadcast_to, - relevant_arguments=[array], - args=(array,), - kwargs=dict(shape=shape, subok=subok)) - if success: - return value + ... # 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-decorating-and-pickling.html>`_. + +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. +- ``np.einsum`` does complicated argument parsing to handle two different + function signatures. It would probably be best to avoid the overhead of + parsing it twice in the typical case of no overrides. + +Fortunately, in each of these cases so far, the functions already has a generic +signature of the form ``*args, **kwargs``, which means we don't need to worry +about potential inconsistency between how functions are called and what we pass +to ``__array_function__``. (In C, arguments for all Python functions are parsed +from a tuple ``*args`` and dict ``**kwargs``.) This shouldn't stop us from +writing overrides for functions with non-generic signatures that can't use the +decorator, but we should consider these cases carefully. + +Extensibility +~~~~~~~~~~~~~ - ... # continue with the definition of broadcast_to +An important virtue of this approach is that it allows for adding new +optional arguments to NumPy functions without breaking code that already +relies on ``__array_function__``. - def concatenate(arrays, axis=0, out=None) - success, value = do_array_function_dance( - func=concatenate, - relevant_arguments=[arrays, out], - args=(arrays,), - kwargs=dict(axis=axis, out=out)) - if success: - return value +This is not a theoretical concern. The implementation of overrides *within* +functions like ``np.sum()`` rather than defining a new function capturing +``*args`` and ``**kwargs`` necessitated some awkward gymnastics to ensure that +the new ``keepdims`` argument is only passed in cases where it is used, e.g., - ... # continue with the definition of concatenate +.. code:: python -The list of objects passed to ``relevant_arguments`` are those that should -be inspected for ``__array_function__`` implementations. + def sum(array, ..., keepdims=np._NoValue): + kwargs = {} + if keepdims is not np._NoValue: + kwargs['keepdims'] = keepdims + return array.sum(..., **kwargs) -Alternatively, we could write these overloads with a decorator, e.g., +This also makes it possible to add optional arguments to ``__array_function__`` +implementations incrementally and only in cases where it makes sense. For +example, a library implementing immutable arrays would not be required to +explicitly include an unsupported ``out`` argument. Doing this properly for all +optional arguments is somewhat onerous, e.g., .. code:: python - @overload_for_array_function(['array']) - def broadcast_to(array, shape, subok=False): - ... # continue with the definition of broadcast_to - - @overload_for_array_function(['arrays', 'out']) - def concatenate(arrays, axis=0, out=None): - ... # continue with the definition of concatenate - -The decorator ``overload_for_array_function`` would be written in terms -of ``do_array_function_dance``. + def my_sum(array, ..., out=None): + if out is not None: + raise TypeError('out argument is not supported') + ... -The downside of this approach would be a loss of introspection capability -for NumPy functions on Python 2, since this requires the use of -``inspect.Signature`` (only available on Python 3). However, NumPy won't -be supporting Python 2 for `very much longer <http://www.numpy.org/neps/nep-0014-dropping-python2.7-proposal.html>`_. +We thus avoid encouraging the tempting shortcut of adding catch-all +``**ignored_kwargs`` to the signatures of functions called by NumPy, which fails +silently for misspelled or ignored arguments. + +Performance +~~~~~~~~~~~ + +Performance is always a concern with NumPy, even though NumPy users have +already prioritized usability over pure speed with their choice of the Python +language itself. It's important that this new ``__array_function__`` protocol +not impose a significant cost in the typical case of NumPy functions acting +on NumPy arrays. + +Our `microbenchmark results <https://nbviewer.jupyter.org/gist/shoyer/1f0a308a06cd96df20879a1ddb8f0006>`_ +show that a pure Python implementation of the override machinery described +above adds roughly 2-3 microseconds of overhead to each NumPy function call +without any overloaded arguments. For context, typical NumPy functions on small +arrays have a runtime of 1-10 microseconds, mostly determined by what fraction +of the function's logic is written in C. For example, one microsecond is about +the difference in speed between the ``ndarray.sum()`` method (1.6 us) and +``numpy.sum()`` function (2.6 us). + +Fortunately, we expect significantly less overhead with a C implementation of +``try_array_function_override``, which is where the bulk of the runtime is. +This would leave the ``array_function_dispatch`` decorator and dispatcher +function on their own adding about 0.5 microseconds of overhead, for perhaps ~1 +microsecond of overhead in the typical case. + +In our view, this level of overhead is reasonable to accept for code written +in Python. We're pretty sure that the vast majority of NumPy users aren't +concerned about performance differences measured in microsecond(s) on NumPy +functions, because it's difficult to do *anything* in Python in less than a +microsecond. Use outside of NumPy ~~~~~~~~~~~~~~~~~~~~ Nothing about this protocol that is particular to NumPy itself. Should -we enourage use of the same ``__array_function__`` protocol third-party +we encourage use of the same ``__array_function__`` protocol third-party libraries for overloading non-NumPy functions, e.g., for making array-implementation generic functionality in SciPy? @@ -276,8 +475,9 @@ to be explicitly recognized. Libraries like Dask, CuPy, and Autograd already wrap a limited subset of SciPy functionality (e.g., ``scipy.linalg``) similarly to how they wrap NumPy. -If we want to do this, we should consider exposing the helper function -``do_array_function_dance()`` above as a public API. +If we want to do this, we should expose at least the decorator +``array_function_dispatch()`` and possibly also the lower level +``try_array_function_override()`` as part of NumPy's public API. Non-goals --------- @@ -332,7 +532,7 @@ Specialized protocols ~~~~~~~~~~~~~~~~~~~~~ We could (and should) continue to develop protocols like -``__array_ufunc__`` for cohesive subsets of Numpy functionality. +``__array_ufunc__`` for cohesive subsets of NumPy functionality. As mentioned above, if this means that some functions that we overload with ``__array_function__`` should switch to a new protocol instead, @@ -347,7 +547,7 @@ either inside or outside of NumPy. This has the advantage of alleviating any possible concerns about backwards compatibility and would provide the maximum freedom for quick -experimentation. In the long term, it would provide a clean abstration +experimentation. In the long term, it would provide a clean abstraction layer, separating NumPy's high level API from default implementations on ``numpy.ndarray`` objects. @@ -358,6 +558,11 @@ functions from ``numpy`` itself are already overloaded (but inadequately), so confusion about high vs. low level APIs in NumPy would still persist. +Alternatively, a separate namespace, e.g., ``numpy.array_only``, could be +created for a non-overloaded version of NumPy's high level API, for cases +where performance with NumPy arrays is a critical concern. This has most +of the same downsides as the separate namespace. + Multiple dispatch ~~~~~~~~~~~~~~~~~ @@ -370,7 +575,7 @@ don't think this approach makes sense for NumPy in the near term. The main reason is that NumPy already has a well-proven dispatching mechanism with ``__array_ufunc__``, based on Python's own dispatching -system for arithemtic, and it would be confusing to add another +system for arithmetic, and it would be confusing to add another mechanism that works in a very different way. This would also be more invasive change to NumPy itself, which would need to gain a multiple dispatch implementation. @@ -384,36 +589,45 @@ would be straightforward to write a shim for a default Implementations in terms of a limited core API ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -The internal implemenations of some NumPy functions is extremely simple. -For example: - ``np.stack()`` is implemented in only a few lines of code -by combining indexing with ``np.newaxis``, ``np.concatenate`` and the -``shape`` attribute. - ``np.mean()`` is implemented internally in terms -of ``np.sum()``, ``np.divide()``, ``.astype()`` and ``.shape``. +The internal implementations of some NumPy functions is extremely simple. +For example: + +- ``np.stack()`` is implemented in only a few lines of code by combining + indexing with ``np.newaxis``, ``np.concatenate`` and the ``shape`` attribute. +- ``np.mean()`` is implemented internally in terms of ``np.sum()``, + ``np.divide()``, ``.astype()`` and ``.shape``. This suggests the possibility of defining a minimal "core" ndarray interface, and relying upon it internally in NumPy to implement the full API. This is an attractive option, because it could significantly reduce the work required for new array implementations. -However, this also comes with several downsides: 1. The details of how -NumPy implements a high-level function in terms of overloaded functions -now becomes an implicit part of NumPy's public API. For example, -refactoring ``stack`` to use ``np.block()`` instead of -``np.concatenate()`` internally would now become a breaking change. 2. -Array libraries may prefer to implement high level functions differently -than NumPy. For example, a library might prefer to implement a -fundamental operations like ``mean()`` directly rather than relying on -``sum()`` followed by division. More generally, it's not clear yet what -exactly qualifies as core functionality, and figuring this out could be -a large project. 3. We don't yet have an overloading system for -attributes and methods on array objects, e.g., for accessing ``.dtype`` -and ``.shape``. This should be the subject of a future NEP, but until -then we should be reluctant to rely on these properties. - -Given these concerns, we encourage relying on this approach only in -limited cases. - -Coersion to a NumPy array as a catch-all fallback +However, this also comes with several downsides: + +1. The details of how NumPy implements a high-level function in terms of + overloaded functions now becomes an implicit part of NumPy's public API. For + example, refactoring ``stack`` to use ``np.block()`` instead of + ``np.concatenate()`` internally would now become a breaking change. +2. Array libraries may prefer to implement high level functions differently than + NumPy. For example, a library might prefer to implement a fundamental + operations like ``mean()`` directly rather than relying on ``sum()`` followed + by division. More generally, it's not clear yet what exactly qualifies as + core functionality, and figuring this out could be a large project. +3. We don't yet have an overloading system for attributes and methods on array + objects, e.g., for accessing ``.dtype`` and ``.shape``. This should be the + subject of a future NEP, but until then we should be reluctant to rely on + these properties. + +Given these concerns, we think it's valuable to support explicit overloading of +nearly every public function in NumPy's API. This does not preclude the future +possibility of rewriting NumPy functions in terms of simplified core +functionality with ``__array_function__`` and a protocol and/or base class for +ensuring that arrays expose methods and properties like ``numpy.ndarray``. +However, to work well this would require the possibility of implementing +*some* but not all functions with ``__array_function__``, e.g., as described +in the next section. + +Coercion to a NumPy array as a catch-all fallback ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ With the current design, classes that implement ``__array_function__`` @@ -438,106 +652,195 @@ make it impossible to implement this generic fallback behavior for ``__array_function__``. We could resolve this issue by change the handling of return values in -``__array_function__`` in either of two possible ways: 1. Change the -meaning of all arguments returning ``NotImplemented`` to indicate that -all arguments should be coerced to NumPy arrays instead. However, many -array libraries (e.g., scipy.sparse) really don't want implicit -conversions to NumPy arrays, and often avoid implementing ``__array__`` -for exactly this reason. Implicit conversions can result in silent bugs -and performance degradation. 2. Use another sentinel value of some sort -to indicate that a class implementing part of the higher level array API -is coercible as a fallback, e.g., a return value of -``np.NotImplementedButCoercible`` from ``__array_function__``. - -If we take this second approach, we would need to define additional -rules for how coercible array arguments are coerced, e.g., - Would we -try for ``__array_function__`` overloads again after coercing coercible -arguments? - If so, would we coerce coercible arguments one-at-a-time, -or all-at-once? - -These are slightly tricky design questions, so for now we propose to -defer this issue. We can always implement -``np.NotImplementedButCoercible`` at some later time if it proves -critical to the numpy community in the future. Importantly, we don't -think this will stop critical libraries that desire to implement most of -the high level NumPy API from adopting this proposal. - -NOTE: If you are reading this NEP in its draft state and disagree, -please speak up on the mailing list! - -Drawbacks of this approach --------------------------- - -Future difficulty extending NumPy's API -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +``__array_function__`` in either of two possible ways: + +1. Change the meaning of all arguments returning ``NotImplemented`` to indicate + that all arguments should be coerced to NumPy arrays and the operation + should be retried. However, many array libraries (e.g., scipy.sparse) really + don't want implicit conversions to NumPy arrays, and often avoid implementing + ``__array__`` for exactly this reason. Implicit conversions can result in + silent bugs and performance degradation. + + Potentially, we could enable this behavior only for types that implement + ``__array__``, which would resolve the most problematic cases like + scipy.sparse. But in practice, a large fraction of classes that present a + high level API like NumPy arrays already implement ``__array__``. This would + preclude reliable use of NumPy's high level API on these objects. +2. Use another sentinel value of some sort, e.g., + ``np.NotImplementedButCoercible``, to indicate that a class implementing part + of NumPy's higher level array API is coercible as a fallback. This is a more + appealing option. + +With either approach, we would need to define additional rules for *how* +coercible array arguments are coerced. The only sane rule would be to treat +these return values as equivalent to not defining an +``__array_function__`` method at all, which means that NumPy functions would +fall-back to their current behavior of coercing all array-like arguments. + +It is not yet clear to us yet if we need an optional like +``NotImplementedButCoercible``, so for now we propose to defer this issue. +We can always implement ``np.NotImplementedButCoercible`` at some later time if +it proves critical to the NumPy community in the future. Importantly, we don't +think this will stop critical libraries that desire to implement most of the +high level NumPy API from adopting this proposal. + +A magic decorator that inspects type annotations +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -One downside of passing on all arguments directly on to -``__array_function__`` is that it makes it hard to extend the signatures -of overloaded NumPy functions with new arguments, because adding even an -optional keyword argument would break existing overloads. +In principle, Python 3 type annotations contain sufficient information to +automatically create most ``dispatcher`` functions. It would be convenient to +use these annotations to dispense with the need for manually writing +dispatchers, e.g., -This is not a new problem for NumPy. NumPy has occasionally changed the -signature for functions in the past, including functions like -``numpy.sum`` which support overloads. +.. code:: python + + @array_function_dispatch + def broadcast_to(array: ArrayLike + shape: Tuple[int, ...], + subok: bool = False): + ... # existing definition of np.broadcast_to + +This would require some form of automatic code generation, either at compile or +import time. + +We think this is an interesting possible extension to consider in the future. We +don't think it makes sense to do so now, because code generation involves +tradeoffs and NumPy's experience with type annotations is still +`quite limited <https://github.com/numpy/numpy-stubs>`_. Even if NumPy +was Python 3 only (which will happen +`sometime in 2019 <http://www.numpy.org/neps/nep-0014-dropping-python2.7-proposal.html>`_), +we aren't ready to annotate NumPy's codebase directly yet. -For adding new keyword arguments that do not change default behavior, we -would only include these as keyword arguments when they have changed -from default values. This is similar to `what NumPy already has -done <https://github.com/numpy/numpy/blob/v1.14.2/numpy/core/fromnumeric.py#L1865-L1867>`_, -e.g., for the optional ``keepdims`` argument in ``sum``: +Support for implementation-specific arguments +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +We could allow ``__array_function__`` implementations to add their own +optional keyword arguments by including ``**ignored_kwargs`` in dispatcher +functions, e.g., .. code:: python - def sum(array, ..., keepdims=np._NoValue): - kwargs = {} - if keepdims is not np._NoValue: - kwargs['keepdims'] = keepdims - return array.sum(..., **kwargs) + def _concatenate_dispatcher(arrays, axis=None, out=None, **ignored_kwargs): + ... # same implementation of _concatenate_dispatcher as above + +Implementation-specific arguments are somewhat common in libraries that +otherwise emulate NumPy's higher level API (e.g., ``dask.array.sum()`` adds +``split_every`` and ``tensorflow.reduce_sum()`` adds ``name``). Supporting +them in NumPy would be particularly useful for libraries that implement new +high-level array functions on top of NumPy functions, e.g., + +.. code:: python + + def mean_squared_error(x, y, **kwargs): + return np.mean((x - y) ** 2, **kwargs) -In other cases, such as deprecated arguments, preserving the existing -behavior of overloaded functions may not be possible. Libraries that use -``__array_function__`` should be aware of this risk: we don't propose to -freeze NumPy's API in stone any more than it already is. +Otherwise, we would need separate versions of ``mean_squared_error`` for each +array implementation in order to pass implementation-specific arguments to +``mean()``. -Difficulty adding implementation specific arguments -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +We wouldn't allow adding optional positional arguments, because these are +reserved for future use by NumPy itself, but conflicts between keyword arguments +should be relatively rare. -Some array implementations generally follow NumPy's API, but have -additional optional keyword arguments (e.g., ``dask.array.sum()`` has -``split_every`` and ``tensorflow.reduce_sum()`` has ``name``). A generic -dispatching library could potentially pass on all unrecognized keyword -argument directly to the implementation, but extending ``np.sum()`` to -pass on ``**kwargs`` would entail public facing changes in NumPy. -Customizing the detailed behavior of array libraries will require using -library specific functions, which could be limiting in the case of -libraries that consume the NumPy API such as xarray. +However, this flexibility would come with a cost. In particular, it implicitly +adds ``**kwargs`` to the signature for all wrapped NumPy functions without +actually including it (because we use ``functools.wraps``). This means it is +unlikely to work well with static analysis tools, which could report invalid +arguments. Likewise, there is a price in readability: these optional arguments +won't be included in the docstrings for NumPy functions. +It's not clear that this tradeoff is worth it, so we propose to leave this out +for now. Adding implementation-specific arguments will require using those +libraries directly. + +Other possible choices for the protocol +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The array function ``__array_function__`` includes only two arguments, ``func`` +and ``types``, that provide information about the context of the function call. + +``func`` is part of the protocol because there is no way to avoid it: +implementations need to be able to dispatch by matching a function to NumPy's +public API. + +``types`` is included because we can compute it almost for free as part of +collecting ``__array_function__`` implementations to call in +``try_array_function_override``. We also think it will be used by most +``__array_function__`` methods, which otherwise would need to extract this +information themselves. It would be equivalently easy to provide single +instances of each type, but providing only types seemed cleaner. + +Taking this even further, it was suggested that ``__array_function__`` should be +a ``classmethod``. We agree that it would be a little cleaner to remove the +redundant ``self`` argument, but feel that this minor clean-up would not be +worth breaking from the precedence of ``__array_ufunc__``. + +There are two other arguments that we think *might* be important to pass to +``__array_ufunc__`` implementations: + +- Access to the non-dispatched function (i.e., before wrapping with + ``array_function_dispatch``) in ``ndarray.__array_function__`` would allow + use to drop special case logic for that method from + ``try_array_function_override``. +- Access to the ``dispatcher`` function passed into + ``array_function_dispatch()`` would allow ``__array_function__`` + implementations to determine the list of "array-like" arguments in a generic + way by calling ``dispatcher(*args, **kwargs)``. This *could* be useful for + ``__array_function__`` implementations that dispatch based on the value of an + array attribute (e.g., ``dtype`` or ``units``) rather than directly on the + array type. + +We have left these out for now, because we don't know that they are necessary. +If we want to include them in the future, the easiest way to do so would be to +update the ``array_function_dispatch`` decorator to add them as function +attributes. + +Callable objects generated at runtime +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +NumPy has some APIs that define callable objects *dynamically*, such as +``vectorize`` and methods on ``random.RandomState`` object. Examples can +also be found in other core libraries in the scientific Python stack, e.g., +distribution objects in scipy.stats and model objects in scikit-learn. It would +be nice to be able to write overloads for such callables, too. This presents a +challenge for the ``__array_function__`` protocol, because unlike the case for +functions there is no public object in the ``numpy`` namespace to pass into +the ``func`` argument. + +We could potentially handle this by establishing an alternative convention +for how the ``func`` argument could be inspected, e.g., by using +``func.__self__`` to obtain the class object and ``func.__func__`` to return +the unbound function object. However, some caution is in order, because +this would immesh what are currently implementation details as a permanent +features of the interface, such as the fact that ``vectorize`` is implemented as a +class rather than closure, or whether a method is implemented directly or using +a descriptor. + +Given the complexity and the limited use cases, we are also deferring on this +issue for now, but we are confident that ``__array_function__`` could be +expanded to accomodate these use cases in the future if need be. Discussion ---------- -Various alternatives to this proposal were discussed in a few Github issues: +Various alternatives to this proposal were discussed in a few GitHub issues: -1. `pydata/sparse #1 <https://github.com/pydata/sparse/issues/1>`_ -2. `numpy/numpy #11129 <https://github.com/numpy/numpy/issues/11129>`_ +1. `pydata/sparse #1 <https://github.com/pydata/sparse/issues/1>`_ +2. `numpy/numpy #11129 <https://github.com/numpy/numpy/issues/11129>`_ Additionally it was the subject of `a blogpost -<http://matthewrocklin.com/blog/work/2018/05/27/beyond-numpy>`_ Following this +<http://matthewrocklin.com/blog/work/2018/05/27/beyond-numpy>`_. Following this it was discussed at a `NumPy developer sprint <https://scisprints.github.io/#may-numpy-developer-sprint>`_ at the `UC Berkeley Institute for Data Science (BIDS) <https://bids.berkeley.edu/>`_. - -References and Footnotes ------------------------- - -.. [1] Each NEP must either be explicitly labeled as placed in the public domain (see - this NEP as an example) or licensed under the `Open Publication License`_. - -.. _Open Publication License: http://www.opencontent.org/openpub/ - +Detailed discussion of this proposal itself can be found on the +`the mailing list <https://mail.python.org/pipermail/numpy-discussion/2018-June/078127.html>`_ and relvant pull requests +(`1 <https://github.com/numpy/numpy/pull/11189>`_, +`2 <https://github.com/numpy/numpy/pull/11303#issuecomment-396638175>`_, +`3 <https://github.com/numpy/numpy/pull/11374>`_) Copyright --------- -This document has been placed in the public domain. [1]_ +This document has been placed in the public domain. diff --git a/doc/neps/nep-0019-rng-policy.rst b/doc/neps/nep-0019-rng-policy.rst index de9164bba..f50897b0f 100644 --- a/doc/neps/nep-0019-rng-policy.rst +++ b/doc/neps/nep-0019-rng-policy.rst @@ -1,12 +1,12 @@ -============================== -Random Number Generator Policy -============================== +======================================= +NEP 19 — Random Number Generator Policy +======================================= :Author: Robert Kern <robert.kern@gmail.com> -:Status: Draft +:Status: Accepted :Type: Standards Track :Created: 2018-05-24 - +:Resolution: https://mail.python.org/pipermail/numpy-discussion/2018-June/078126.html Abstract -------- @@ -91,23 +91,12 @@ 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. +Work on a proposed new PRNG subsystem is already underway in the randomgen_ +project. 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. We will also outline just a few of the +requirements that such a new system must have to support the policy proposed in +this NEP. 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 @@ -116,66 +105,158 @@ 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. +cadence of features (i.e. on ``X.Y`` releases, never ``X.Y.Z``). Slowness will +not be considered a bug for this purpose. 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 a +discouraged improvement would be tweaking the Ziggurat tables just a little bit +for a small performance improvement. 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, +uniform PRNG algorithms. A promising design choice is to make these core +uniform PRNGs their own lightweight objects with a minimal set of methods +(randomgen_ calls them “basic RNGs”). The broader set of non-uniform +distributions will be its own class that holds a reference to one of these core +uniform PRNG objects and simply delegates to the core uniform PRNG object when +it needs uniform random numbers. To borrow an example from randomgen_, the +class ``MT19937`` is a basic RNG that implements the classic Mersenne Twister +algorithm. The class ``RandomGenerator`` wraps around the basic RNG to provide +all of the non-uniform distribution methods:: + + # This is not the only way to instantiate this object. + # This is just handy for demonstrating the delegation. + >>> brng = MT19937(seed) + >>> rg = RandomGenerator(brng) + >>> x = rg.standard_normal(10) + +We will be more strict about a select subset of methods on these basic RNG +objects. They MUST guarantee stream-compatibility for a specified set +of methods which are chosen to make it easier to compose them to build other +distributions and which are needed to abstract over the implementation details +of the variety of core PRNG algorithms. 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 --------------- +The distributions class (``RandomGenerator``) SHOULD have all of the same +distribution methods as ``RandomState`` with close-enough function signatures +such that almost all code that currently works with ``RandomState`` instances +will work with ``RandomGenerator`` instances (ignoring the precise stream +values). Some variance will be allowed for integer distributions: in order to +avoid some of the cross-platform problems described above, these SHOULD be +rewritten to work with ``uint64`` numbers on all platforms. + +.. _randomgen: https://github.com/bashtage/randomgen + + +Supporting Unit Tests +::::::::::::::::::::: + +Because we did make a strong stream-compatibility guarantee early in numpy’s +life, reliance on stream-compatibility has grown beyond reproducible +simulations. One use case that remains for stream-compatibility across numpy +versions is to use pseudorandom streams to generate test data in unit tests. +With care, many of the cross-platform instabilities can be avoided in the +context of small unit tests. + +The new PRNG subsystem MUST provide a second, legacy distributions class that +uses the same implementations of the distribution methods as the current +version of ``numpy.random.RandomState``. The methods of this class will have +strict stream-compatibility guarantees, even stricter than the current policy. +It is intended that this class will no longer be modified, except to keep it +working when numpy internals change. All new development should go into the +primary distributions class. Bug fixes that change the stream SHALL NOT be +made to ``RandomState``; instead, buggy distributions should be made to warn +when they are buggy. The purpose of ``RandomState`` will be documented as +providing certain fixed functionality for backwards compatibility and stable +numbers for the limited purpose of unit testing, and not making whole programs +reproducible across numpy versions. + +This legacy distributions class MUST be accessible under the name +``numpy.random.RandomState`` for backwards compatibility. All current ways of +instantiating ``numpy.random.RandomState`` with a given state should +instantiate the Mersenne Twister basic RNG with the same state. The legacy +distributions class MUST be capable of accepting other basic RNGs. The purpose +here is to ensure that one can write a program with a consistent basic RNG +state with a mixture of libraries that may or may not have upgraded from +``RandomState``. Instances of the legacy distributions class MUST respond +``True`` to ``isinstance(rg, numpy.random.RandomState)`` because there is +current utility code that relies on that check. Similarly, old pickles of +``numpy.random.RandomState`` instances MUST unpickle correctly. + + +``numpy.random.*`` +:::::::::::::::::: + +The preferred best practice for getting reproducible pseudorandom numbers is to +instantiate a generator object with a seed and pass it around. The implicit +global ``RandomState`` behind the ``numpy.random.*`` convenience functions can +cause problems, especially when threads or other forms of concurrency are +involved. Global state is always problematic. We categorically recommend +avoiding using the convenience functions when reproducibility is involved. + +That said, people do use them and use ``numpy.random.seed()`` to control the +state underneath them. It can be hard to categorize and count API usages +consistently and usefully, but a very common usage is in unit tests where many +of the problems of global state are less likely. + +This NEP does not propose removing these functions or changing them to use the +less-stable ``RandomGenerator`` distribution implementations. Future NEPs +might. + +Specifically, the initial release of the new PRNG subsystem SHALL leave these +convenience functions as aliases to the methods on a global ``RandomState`` +that is initialized with a Mersenne Twister basic RNG object. A call to +``numpy.random.seed()`` will be forwarded to that basic RNG object. In +addition, the global ``RandomState`` instance MUST be accessible in this +initial release by the name ``numpy.random.mtrand._rand``: Robert Kern long ago +promised ``scikit-learn`` that this name would be stable. Whoops. + +In order to allow certain workarounds, it MUST be possible to replace the basic +RNG underneath the global ``RandomState`` with any other basic RNG object (we +leave the precise API details up to the new subsystem). Calling +``numpy.random.seed()`` thereafter SHOULD just pass the given seed to the +current basic RNG object and not attempt to reset the basic RNG to the Mersenne +Twister. The set of ``numpy.random.*`` convenience functions SHALL remain the +same as they currently are. They SHALL be aliases to the ``RandomState`` +methods and not the new less-stable distributions class (``RandomGenerator``, +in the examples above). Users who want to get the fastest, best distributions +can follow best practices and instantiate generator objects explicitly. + +This NEP does not propose that these requirements remain in perpetuity. After +we have experience with the new PRNG subsystem, we can and should revisit these +issues in future NEPs. + + +Alternatives +------------ + +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. +of the code, and there would be a way to get the old versions. + +We will not be doing this. If one needs to get the exact bit-for-bit results +from a given version of ``numpy``, whether one uses random numbers or not, one +should use the exact version of ``numpy``. + +Proposals of how to do RNG versioning varied widely, and 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 @@ -195,11 +276,49 @@ is to pin the release of ``numpy`` as a whole, versioning ``RandomState`` alone is superfluous. +``StableRandom`` +:::::::::::::::: + +A previous version of this NEP proposed to leave ``RandomState`` completely +alone for a deprecation period and build the new subsystem alongside with new +names. To satisfy the unit testing use case, it proposed introducing a small +distributions class nominally called ``StableRandom``. It would have provided +a small subset of distribution methods that were considered most useful in unit +testing, but not the full set such that it would be too likely to be used +outside of the testing context. + +During discussion about this proposal, it became apparent that there was no +satisfactory subset. At least some projects used a fairly broad selection of +the ``RandomState`` methods in unit tests. + +Downstream project owners would have been forced to modify their code to +accomodate the new PRNG subsystem. Some modifications might be simply +mechanical, but the bulk of the work would have been tedious churn for no +positive improvement to the downstream project, just avoiding being broken. + +Furthermore, under this old proposal, we would have had a quite lengthy +deprecation period where ``RandomState`` existed alongside the new system of +basic RNGs and distribution classes. Leaving the implementation of +``RandomState`` fixed meant that it could not use the new basic RNG state +objects. Developing programs that use a mixture of libraries that have and +have not upgraded would require managing two sets of PRNG states. This would +notionally have been time-limited, but we intended the deprecation to be very +long. + +The current proposal solves all of these problems. All current usages of +``RandomState`` will continue to work in perpetuity, though some may be +discouraged through documentation. Unit tests can continue to use the full +complement of ``RandomState`` methods. Mixed ``RandomState/RandomGenerator`` +code can safely share the common basic RNG state. Unmodified ``RandomState`` +code can make use of the new features of alternative basic RNGs like settable +streams. + + Discussion ---------- -- https://mail.python.org/pipermail/numpy-discussion/2018-January/077608.html -- https://github.com/numpy/numpy/pull/10124#issuecomment-350876221 +- `NEP discussion <https://mail.python.org/pipermail/numpy-discussion/2018-June/078126.html>`_ +- `Earlier discussion <https://mail.python.org/pipermail/numpy-discussion/2018-January/077608.html>`_ Copyright diff --git a/doc/neps/nep-0020-gufunc-signature-enhancement.rst b/doc/neps/nep-0020-gufunc-signature-enhancement.rst new file mode 100644 index 000000000..38a9fd53b --- /dev/null +++ b/doc/neps/nep-0020-gufunc-signature-enhancement.rst @@ -0,0 +1,257 @@ +=============================================================== +NEP 20 — Expansion of Generalized Universal Function Signatures +=============================================================== + +:Author: Marten van Kerkwijk <mhvk@astro.utoronto.ca> +:Status: Accepted +:Type: Standards Track +:Created: 2018-06-10 +:Resolution: https://mail.python.org/pipermail/numpy-discussion/2018-April/077959.html, + https://mail.python.org/pipermail/numpy-discussion/2018-May/078078.html + +.. note:: The proposal to add fixed (i) and flexible (ii) dimensions + was accepted, while that to add broadcastable (iii) ones was deferred. + +Abstract +-------- + +Generalized universal functions are, as their name indicates, generalization +of universal functions: they operate on non-scalar elements. Their signature +describes the structure of the elements they operate on, with names linking +dimensions of the operands that should be the same. Here, it is proposed to +extend the signature to allow the signature to indicate that a dimension (i) +has fixed size; (ii) can be absent; and (iii) can be broadcast. + +Detailed description +-------------------- + +Each part of the proposal is driven by specific needs [1]_. + +1. Fixed-size dimensions. Code working with spatial vectors often explicitly + is for 2 or 3-dimensional space (e.g., the code from the `Standards Of + Fundamental Astronomy <http://www.iausofa.org/>`_, which the author hopes + to wrap using gufuncs for astropy [2]_). The signature should be able to + indicate that. E.g., the signature of a function that converts a polar + angle to a two-dimensional cartesian unit vector would currently have to be + ``()->(n)``, with there being no way to indicate that ``n`` has to equal 2. + Indeed, this signature is particularly annoying since without putting in an + output argument, the current gufunc wrapper code fails because it cannot + determine ``n``. Similarly, the signature for an cross product of two + 3-dimensional vectors has to be ``(n),(n)->(n)``, with again no way to + indicate that ``n`` has to equal 3. Hence, the proposal here to allow one + to give numerical values in addition to variable names. Thus, angle to + two-dimensional unit vector would be ``()->(2)``; two angles to + three-dimensional unit vector ``(),()->(3)``; and that for the cross + product of two three-dimensional vectors would be ``(3),(3)->(3)``. + +2. Possibly missing dimensions. This part is almost entirely driven by the + wish to wrap ``matmul`` in a gufunc. ``matmul`` stands for matrix + multiplication, and if it did only that, it could be covered with the + signature ``(m,n),(n,p)->(m,p)``. However, it has special cases for when a + dimension is missing, allowing either argument to be treated as a single + vector, with the function thus becoming, effectively, vector-matrix, + matrix-vector, or vector-vector multiplication (but with no + broadcasting). To support this, it is suggested to allow postfixing a + dimension name with a question mark to indicate that the dimension does not + necessarily have to be present. + + With this addition, the signature for ``matmul`` can be expressed as + ``(m?,n),(n,p?)->(m?,p?)``. This indicates that if, e.g., the second + operand has only one dimension, for the purposes of the elementary function + it will be treated as if that input has core shape ``(n, 1)``, and the + output has the corresponding core shape of ``(m, 1)``. The actual output + array, however, has the flexible dimension removed, i.e., it will have + shape ``(..., m)``. Similarly, if both arguments have only a single + dimension, the inputs will be presented as having shapes ``(1, n)`` and + ``(n, 1)`` to the elementary function, and the output as ``(1, 1)``, while + the actual output array returned will have shape ``()``. In this way, the + signature allows one to use a single elementary function for four related + but different signatures, ``(m,n),(n,p)->(m,p)``, ``(n),(n,p)->(p)``, + ``(m,n),(n)->(m)`` and ``(n),(n)->()``. + +3. Dimensions that can be broadcast. For some applications, broadcasting + between operands makes sense. For instance, an ``all_equal`` function that + compares vectors in arrays could have a signature ``(n),(n)->()``, but this + forces both operands to be arrays, while it would be useful also to check + that, e.g., all parts of a vector are constant (maybe zero). The proposal + is to allow the implementer of a gufunc to indicate that a dimension can be + broadcast by post-fixing the dimension name with ``|1``. Hence, the + signature for ``all_equal`` would become ``(n|1),(n|1)->()``. The + signature seems handy more generally for "chained ufuncs"; e.g., another + application might be in a putative ufunc implementing ``sumproduct``. + + Another example that arose in the discussion, is of a weighted mean, which + might look like ``weighted_mean(y, sigma[, axis, ...])``, returning the + mean and its uncertainty. With a signature of ``(n),(n)->(),()``, one + would be forced to always give as many sigmas as there are data points, + while broadcasting would allow one to give a single sigma for all points + (which is still useful to calculate the uncertainty on the mean). + +Implementation +-------------- + +The proposed changes have all been implemented [3]_, [4]_, [5]_. These PRs +extend the ufunc structure with two new fields, each of size equal to the +number of distinct dimensions, with ``core_dim_sizes`` holding possibly fixed +sizes, and ``core_dim_flags`` holding flags indicating whether a dimension can +be missing or broadcast. To ensure we can distinguish between this new +version and previous versions, an unused entry ``reserved1`` is repurposed as +a version number. + +In the implementation, care is taken that to the elementary function flagged +dimensions are not treated any differently than non-flagged ones: for +instance, sizes of fixed-size dimensions are still passed on to the elementary +function (but the loop can now count on that size being equal to the fixed one +given in the signature). + +An implementation detail to be decided upon is whether it might be handy to +have a summary of all flags. This could possibly be stored in ``core_enabled`` +(which currently is a bool), with non-zero continuing to indicate a gufunc, +but specific flags indicating whether or not a gufunc uses fixed, flexible, or +broadcastable dimensions. + +With the above, the formal defition of the syntax would become [4]_:: + + <Signature> ::= <Input arguments> "->" <Output arguments> + <Input arguments> ::= <Argument list> + <Output arguments> ::= <Argument list> + <Argument list> ::= nil | <Argument> | <Argument> "," <Argument list> + <Argument> ::= "(" <Core dimension list> ")" + <Core dimension list> ::= nil | <Core dimension> | + <Core dimension> "," <Core dimension list> + <Core dimension> ::= <Dimension name> <Dimension modifier> + <Dimension name> ::= valid Python variable name | valid integer + <Dimension modifier> ::= nil | "|1" | "?" + +#. All quotes are for clarity. +#. Unmodified core dimensions that share the same name must have the same size. + Each dimension name typically corresponds to one level of looping in the + elementary function's implementation. +#. White spaces are ignored. +#. An integer as a dimension name freezes that dimension to the value. +#. If a name if suffixed with the ``|1`` modifier, it is allowed to broadcast + against other dimensions with the same name. All input dimensions + must share this modifier, while no output dimensions should have it. +#. If the name is suffixed with the ``?`` modifier, the dimension is a core + dimension only if it exists on all inputs and outputs that share it; + otherwise it is ignored (and replaced by a dimension of size 1 for the + elementary function). + +Examples of signatures [4]_: + ++----------------------------+-----------------------------------+ +| Signature | Possible use | ++----------------------------+-----------------------------------+ +| ``(),()->()`` | Addition | ++----------------------------+-----------------------------------+ +| ``(i)->()`` | Sum over last axis | ++----------------------------+-----------------------------------+ +| ``(i|1),(i|1)->()`` | Test for equality along axis, | +| | allowing comparison with a scalar | ++----------------------------+-----------------------------------+ +| ``(i),(i)->()`` | inner vector product | ++----------------------------+-----------------------------------+ +| ``(m,n),(n,p)->(m,p)`` | matrix multiplication | ++----------------------------+-----------------------------------+ +| ``(n),(n,p)->(p)`` | vector-matrix multiplication | ++----------------------------+-----------------------------------+ +| ``(m,n),(n)->(m)`` | matrix-vector multiplication | ++----------------------------+-----------------------------------+ +| ``(m?,n),(n,p?)->(m?,p?)`` | all four of the above at once, | +| | except vectors cannot have loop | +| | dimensions (ie, like ``matmul``) | ++----------------------------+-----------------------------------+ +| ``(3),(3)->(3)`` | cross product for 3-vectors | ++----------------------------+-----------------------------------+ +| ``(i,t),(j,t)->(i,j)`` | inner over the last dimension, | +| | outer over the second to last, | +| | and loop/broadcast over the rest. | ++----------------------------+-----------------------------------+ + +Backward compatibility +---------------------- + +One possible worry is the change in ufunc structure. For most applications, +which call ``PyUFunc_FromDataAndSignature``, this is entirely transparent. +Furthermore, by repurposing ``reserved1`` as a version number, code compiled +against older versions of numpy will continue to work (though one will get a +warning upon import of that code with a newer version of numpy), except if +code explicitly changes the ``reserved1`` entry. + +Alternatives +------------ + +It was suggested instead of extending the signature, to have multiple +dispatch, so that, e.g., ``matmul`` would simply have the multiple signatures +it supports, i.e., instead of ``(m?,n),(n,p?)->(m?,p?)`` one would have +``(m,n),(n,p)->(m,p) | (n),(n,p)->(p) | (m,n),(n)->(m) | (n),(n)->()``. A +disadvantage of this is that the developer now has to make sure that the +elementary function can deal with these different signatures. Furthermore, +the expansion quickly becomes cumbersome. For instance, for the ``all_equal`` +signature of ``(n|1),(n|1)->()``, one would have to have five entries: +``(n),(n)->() | (n),(1)->() | (1),(n)->() | (n),()->() | (),(n)->()``. For +signatures like ``(m|1,n|1,o|1),(m|1,n|1,o|1)->()`` (from the ``cube_equal`` +test case in [4]_), it is not even worth writing out the expansion. + +For broadcasting, the alternative suffix of ``^`` was suggested (as +broadcasting can be thought of as increasing the size of the array). This +seems less clear. Furthermore, it was wondered whether it should not just be +an all-or-nothing flag. This could be the case, though given the postfix +for flexible dimensions, arguably another postfix is clearer (as is the +implementation). + +Discussion +---------- + +The proposals here were discussed at fair length on the mailing list [6]_, +[7]_. The main points of contention were whether the use cases were +sufficiently strong. In particular, for frozen dimensions, it was argued that +checks on the right number could be put in loop selection code. This seems +much less clear for no benefit. + +For broadcasting, the lack of examples of elementary functions that might need +it was noted, with it being questioned whether something like ``all_equal`` +was best done with a gufunc rather than as a special method on ``np.equal``. +One counter-argument to this would be that there is an actual PR for +``all_equal`` [8]_. Another that even if one were to use a method, it would +be good to be able to express their signature (just as is possible at least +for ``reduce`` and ``accumulate``). + +A final argument was that we were making the gufuncs too complex. This +arguably holds for the dimensions that can be omitted, but that also has the +strongest use case. The frozen dimensions has a very simple implementation and +its meaning is obvious. The ability to broadcast is simple too, once the +flexible dimensions are supported. + +References and Footnotes +------------------------ + +.. [1] Identified needs and suggestions for the implementation are not all by + the author. In particular, the suggestion for fixed dimensions and + initial implementation was by Jaime Frio (`gh-5015 + <https://github.com/numpy/numpy/pull/5015>`_), the suggestion of ``?`` + to indicate dimensions can be omitted was by Nathaniel Smith, and the + initial implementation of that by Matti Picus (`gh-11132 + <https://github.com/numpy/numpy/pull/11132>`_). +.. [2] `wrap ERFA functions in gufuncs + <https://github.com/astropy/astropy/pull/7502>`_ (`ERFA + <https://github.com/liberfa/erfa>`_) is the less stringently licensed + version of `Standards Of Fundamental Astronomy + <http://www.iausofa.org/>`_ +.. [3] `fixed-size and flexible dimensions + <https://github.com/numpy/numpy/pull/11175>`_ +.. [4] `broadcastable dimensions + <https://github.com/numpy/numpy/pull/11179>`_ +.. [5] `use in matmul <https://github.com/numpy/numpy/pull/11133>`_ +.. [6] Discusses implementations for ``matmul``: + https://mail.python.org/pipermail/numpy-discussion/2018-May/077972.html, + https://mail.python.org/pipermail/numpy-discussion/2018-May/078021.html +.. [7] Broadcasting: + https://mail.python.org/pipermail/numpy-discussion/2018-May/078078.html +.. [8] `Logical gufuncs <https://github.com/numpy/numpy/pull/8528>`_ (includes + ``all_equal``) + +Copyright +--------- + +This document has been placed in the public domain. diff --git a/doc/neps/nep-0021-advanced-indexing.rst b/doc/neps/nep-0021-advanced-indexing.rst new file mode 100644 index 000000000..5acabbf16 --- /dev/null +++ b/doc/neps/nep-0021-advanced-indexing.rst @@ -0,0 +1,661 @@ +================================================== +NEP 21 — Simplified and explicit advanced indexing +================================================== + +:Author: Sebastian Berg +:Author: Stephan Hoyer <shoyer@google.com> +:Status: Draft +:Type: Standards Track +:Created: 2015-08-27 + + +Abstract +-------- + +NumPy's "advanced" indexing support for indexing array 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 ``oindex`` and ``vindex``. + +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 are 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 diff --git a/doc/neps/nep-0022-ndarray-duck-typing-overview.rst b/doc/neps/nep-0022-ndarray-duck-typing-overview.rst new file mode 100644 index 000000000..04e4a14b7 --- /dev/null +++ b/doc/neps/nep-0022-ndarray-duck-typing-overview.rst @@ -0,0 +1,351 @@ +=========================================================== +NEP 22 — Duck typing for NumPy arrays – high level overview +=========================================================== + +:Author: Stephan Hoyer <shoyer@google.com>, Nathaniel J. Smith <njs@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. diff --git a/doc/neps/nep-template.rst b/doc/neps/nep-template.rst index 26515127d..e869ebae3 100644 --- a/doc/neps/nep-template.rst +++ b/doc/neps/nep-template.rst @@ -64,7 +64,7 @@ References and Footnotes .. [1] Each NEP must either be explicitly labeled as placed in the public domain (see this NEP as an example) or licensed under the `Open Publication License`_. -.. _Open Publication License: http://www.opencontent.org/openpub/ +.. _Open Publication License: https://www.opencontent.org/openpub/ Copyright diff --git a/doc/neps/roadmap.rst b/doc/neps/roadmap.rst new file mode 100644 index 000000000..a45423711 --- /dev/null +++ b/doc/neps/roadmap.rst @@ -0,0 +1,115 @@ +============= +NumPy Roadmap +============= + +This is a live snapshot of tasks and features we will be investing resources +in. It may be used to encourage and inspire developers and to search for +funding. + +Interoperability protocols & duck typing +---------------------------------------- + +- `__array_function__` + + See `NEP 18`_ and a sample implementation_ + +- Array Duck-Typing + + `NEP 22`_ `np.asduckarray()` + +- Mixins like `NDArrayOperatorsMixin`: + + - for mutable arrays + - for reduction methods implemented as ufuncs + +Better dtypes +------------- + +- Easier custom dtypes + - Simplify and/or wrap the current C-API + - More consistent support for dtype metadata + - Support for writing a dtype in Python +- New string dtype(s): + - Encoded strings with fixed-width storage (utf8, latin1, ...) and/or + - Variable length strings (could share implementation with dtype=object, but are explicitly type-checked) + - One of these should probably be the default for text data. The current behavior on Python 3 is neither efficient nor user friendly. +- `np.int` should not be platform dependent +- better coercion for string + number + +Random number generation policy & rewrite +----------------------------------------- + +`NEP 19`_ and a `reference implementation`_ + +Indexing +-------- + +vindex/oindex `NEP 21`_ + +Infrastructure +-------------- + +NumPy is much more than just the code base itself, we also maintain +docs, CI, benchmarks, etc. + +- Rewrite numpy.org +- Benchmarking: improve the extent of the existing suite, and run & render + the results as part of the docs or website. + + - Hardware: find a machine that can reliably run serial benchmarks + - ASV produces graphs, could we set up a site? Currently at + https://pv.github.io/numpy-bench/, should that become a community resource? + +Functionality outside core +-------------------------- + +Some things inside NumPy do not actually match the `Scope of NumPy`. + +- A backend system for `numpy.fft` (so that e.g. `fft-mkl` doesn't need to monkeypatch numpy) + +- Rewrite masked arrays to not be a ndarray subclass -- maybe in a separate project? +- MaskedArray as a duck-array type, and/or +- dtypes that support missing values + +- Write a strategy on how to deal with overlap between numpy and scipy for `linalg` and `fft` (and implement it). + +- Deprecate `np.matrix` + +Continuous Integration +---------------------- + +We depend on CI to discover problems as we continue to develop NumPy before the +code reaches downstream users. + +- CI for more exotic platforms (e.g. ARM is now available from + http://www.shippable.com/, but it is not free). +- Multi-package testing +- Add an official channel for numpy dev builds for CI usage by other projects so + they may confirm new builds do not break their package. + +Typing +------ + +Python type annotation syntax should support ndarrays and dtypes. + +- Type annotations for NumPy: github.com/numpy/numpy-stubs +- Support for typing shape and dtype in multi-dimensional arrays in Python more generally + +NumPy scalars +------------- + +Numpy has both scalars and zero-dimensional arrays. + +- The current implementation adds a large maintenance burden -- can we remove + scalars and/or simplify it internally? +- Zero dimensional arrays get converted into scalars by most NumPy + functions (i.e., output of `np.sin(x)` depends on whether `x` is + zero-dimensional or not). This inconsistency should be addressed, + so that one could, e.g., write sane type annotations. + +.. _`NEP 19`: https://www.numpy.org/neps/nep-0019-rng-policy.html +.. _`NEP 22`: http://www.numpy.org/neps/nep-0022-ndarray-duck-typing-overview.html +.. _`NEP 18`: https://www.numpy.org/neps/nep-0018-array-function-protocol.html +.. _implementation: https://gist.github.com/shoyer/1f0a308a06cd96df20879a1ddb8f0006 +.. _`reference implementation`: https://github.com/bashtage/randomgen +.. _`NEP 21`: https://www.numpy.org/neps/nep-0021-advanced-indexing.html diff --git a/doc/neps/scope.rst b/doc/neps/scope.rst new file mode 100644 index 000000000..a675b8c96 --- /dev/null +++ b/doc/neps/scope.rst @@ -0,0 +1,46 @@ +============== +Scope of NumPy +============== + +Here, we describe aspects of N-d array computation that are within scope for NumPy development. This is *not* an aspirational definition of where NumPy should aim, but instead captures the status quo—areas which we have decided to continue supporting, at least for the time being. + +- **In-memory, N-dimensional, homogeneously typed (single pointer + strided) arrays on CPUs** + + - Support for a wide range of data types + - Not specialized hardware such as GPUs + - But, do support wide range of CPUs (e.g. ARM, PowerX) + +- **Higher level APIs for N-dimensional arrays** + + - NumPy is a *de facto* standard for array APIs in Python + - Indexing and fast iteration over elements (ufunc) + - Interoperability protocols with other data container implementations (like `__array_ufunc__`). + +- **Python API and a C API** to the ndarray's methods and attributes. + +- Other **specialized types or uses of N-dimensional arrays**: + + - Masked arrays + - Structured arrays (informally known as record arrays) + - Memory mapped arrays + +- Historically, NumPy has included the following **basic functionality + in support of scientific computation**. We intend to keep supporting + (but not to expand) what is currently included: + + - Linear algebra + - Fast Fourier transforms and windowing + - Pseudo-random number generators + - Polynomial fitting + +- NumPy provides some **infrastructure for other packages in the scientific Python ecosystem**: + + - numpy.distutils (build support for C++, Fortran, BLAS/LAPACK, and other relevant libraries for scientific computing + - f2py (generating bindings for Fortran code) + - testing utilities + +- **Speed**: we take performance concerns seriously and aim to execute + operations on large arrays with similar performance as native C + code. That said, where conflict arises, maintenance and portability take + precedence over performance. We aim to prevent regressions where + possible (e.g., through asv). diff --git a/doc/neps/tools/build_index.py b/doc/neps/tools/build_index.py index 65225c995..d9c4f690b 100644 --- a/doc/neps/tools/build_index.py +++ b/doc/neps/tools/build_index.py @@ -40,6 +40,10 @@ def nep_metadata(): tags['Title'] = lines[1].strip() tags['Filename'] = source + if not tags['Title'].startswith(f'NEP {nr} — '): + raise RuntimeError( + f'Title for NEP {nr} does not start with "NEP {nr} — " ' + '(note that — here is a special, enlongated dash)') if tags['Status'] in ('Accepted', 'Rejected', 'Withdrawn'): if not 'Resolution' in tags: |