NEP-18 Compatibility

CHANGES

@@ -4,6 +4,21 @@ Pint Changelog
 0.10 (unreleased)
 -----------------
 
+- **BREAKING CHANGE**:
+  Implement NEP-18 <https://numpy.org/neps/nep-0018-array-function-protocol.html> for
+  Pint Quantities. Most NumPy functions that previously stripped units when applied to
+  Pint Quantities will now return Quantities with proper units (on NumPy v1.16 with
+  the array_function protocol enabled or v1.17+ by default) instead of ndarrays. Any
+  non-explictly-handled functions will now raise a "no implementation found" TypeError
+  instead of stripping units. The previous behavior is maintained for NumPy < v1.16 and
+  when the array_function protocol is disabled.
+  (Issue #905, Thanks Jon Thielen and andrewgsavage)
+- Implementation of NumPy ufuncs has been refactored to share common utilities with
+  NumPy function implementations
+  (Issue #905, Thanks Jon Thielen)
+- Pint Quantities now support the `@` matrix mulitiplication operator (on NumPy v1.16+),
+  as well as the `dot`, `flatten`, `astype`, and `item` methods.
+  (Issue #905, Thanks Jon Thielen)
 - **BREAKING CHANGE**:
   Fix crash when applying pprint to large sets of Units.
   DefinitionSyntaxError is now a subclass of SyntaxError (was ValueError).

docs/numpy.rst

@@ -1,12 +1,13 @@
 .. _numpy:
 
 
-NumPy support
+NumPy Support
 =============
 
-The magnitude of a Pint quantity can be of any numerical type and you are free
-to choose it according to your needs. In numerical applications, it is quite
-convenient to use `NumPy ndarray`_ and therefore they are supported by Pint.
+The magnitude of a Pint quantity can be of any numerical scalar type, and you are free
+to choose it according to your needs. For numerical applications requiring arrays, it is
+quite convenient to use `NumPy ndarray`_ (or `ndarray-like types supporting NEP-18`_),
+and therefore these are the array types supported by Pint.
 
 First, we import the relevant packages:
 
@@ -81,11 +82,11 @@ is applied before the requested calculation:
     >>> print(angles)
     [ 0.64350111  0.92729522] radian
 
-You can convert the result to degrees using the corresponding NumPy function:
+You can convert the result to degrees using usual unit conversion:
 
 .. doctest::
 
-    >>> print(np.rad2deg(angles))
+    >>> print(angles.to('degree'))
     [ 36.86989765  53.13010235] degree
 
 Applying a function that expects angles to a quantity with a different dimensionality
@@ -109,13 +110,16 @@ The following ufuncs_ can be applied to a Quantity object:
 - **Comparison functions**: greater, greater_equal, less, less_equal, not_equal, equal
 - **Floating functions**: isreal,iscomplex, isfinite, isinf, isnan, signbit, copysign, nextafter, modf, ldexp, frexp, fmod, floor, ceil, trunc
 
-And the following `NumPy ndarray methods`_ and functions:
+And the following NumPy functions:
 
-- sum, fill, reshape, transpose, flatten, ravel, squeeze, take, put, repeat, sort, argsort, diagonal, compress,  nonzero, searchsorted, max, argmax, min, argmin, ptp, clip, round, trace, cumsum, mean, var, std, prod, cumprod, conj, conjugate, flatten
+- alen, amax, amin, append, argmax, argmin, argsort, around, atleast_1d, atleast_2d, atleast_3d, average, block, broadcast_to, clip, column_stack, compress, concatenate, copy, copyto, count_nonzero, cross, cumprod, cumproduct, cumsum, diagonal, diff, dot, dstack, ediff1d, einsum, empty_like, expand_dims, fix, flip, full_like, gradient, hstack, insert, interp, isclose, iscomplex, isin, isreal, linspace, mean, median, meshgrid, moveaxis, nan_to_num, nanargmax, nanargmin, nancumprod, nancumsum, nanmax, nanmean, nanmedian, nanmin, nanpercentile, nanstd, nanvar, ndim, nonzero, ones_like, percentile, ptp, ravel, result_type, rollaxis, rot90, round\_, searchsorted, shape, size, sort, squeeze, stack, std, sum, swapaxes, tile, transpose, trapz, trim_zeros, unwrap, var, vstack, where, zeros_like
 
-`Quantity` is not a subclass of `ndarray`. This might change in the future, but for this reason  functions that call `numpy.asanyarray` are currently not supported. These functions are:
+And the following `NumPy ndarray methods`_:
 
-- unwrap, trapz, diff, ediff1d, fix, gradient, cross, ones_like
+- argmax, argmin, argsort, astype, clip, compress, conj, conjugate, cumprod, cumsum, diagonal, dot, fill, flatten, flatten, item, max, mean, min, nonzero, prod, ptp, put, ravel, repeat, reshape, round, searchsorted, sort, squeeze, std, sum, take, trace, transpose, var
+
+Pull requests are welcome for any NumPy function, ufunc, or method that is not currently
+supported.
 
 
 Comments
@@ -133,22 +137,32 @@ will be raised.
 
 In some functions that take 2 or more arguments (e.g. `arctan2`), the second
 argument is converted to the units of the first. Again, a `DimensionalityError`
-exception will be raised if this is not possible.
+exception will be raised if this is not possible. ndarray or ndarray-like arguments
+are generally treated as if they were dimensionless quantities, except for declared
+upcast types to which Pint defers (see
+<https://numpy.org/neps/nep-0013-ufunc-overrides.html>). To date, these "upcast types" are:
+
+- ``PintArray``, as defined by pint-pandas
+- ``Series``, as defined by pandas
+- ``DataArray``, as defined by xarray
+
+To achive these function and ufunc overrides, Pint uses the ``__array_function__`` and
+``__array_ufunc__`` protocols respectively, as recommened by NumPy. This means that
+functions and ufuncs that Pint does not explicitly handle will error, rather than return
+a value with units stripped (in contrast to Pint's behavior prior to v0.10). For more
+information on these protocols, see
+<https://docs.scipy.org/doc/numpy-1.17.0/user/basics.dispatch.html>.
 
 This behaviour introduces some performance penalties and increased memory
 usage. Quantities that must be converted to other units require additional
-memory and CPU cycles. On top of this, all `ufuncs` are implemented in the
-`Quantity` class by overriding `__array_wrap__`, a NumPy hook that is executed
-after the calculation and before returning the value. To our knowledge, there
-is no way to signal back to NumPy that our code will take care of the
-calculation. For this reason the calculation is actually done twice:
-first in the original ndarray and then in the one that has been
-converted to the right units. Therefore, for numerically intensive code, you
-might want to convert the objects first and then use directly the magnitude.
+memory and CPU cycles. Therefore, for numerically intensive code, you
+might want to convert the objects first and then use directly the magnitude,
+such as by using Pint's `wraps` utility (see :ref:`wrapping`).
 
 
 
 
 .. _`NumPy ndarray`: http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html
+.. _`ndarray-like types supporting NEP-18`: https://numpy.org/neps/nep-0018-array-function-protocol.html
 .. _ufuncs: http://docs.scipy.org/doc/numpy/reference/ufuncs.html
 .. _`NumPy ndarray methods`: http://docs.scipy.org/doc/numpy/reference/arrays.ndarray.html#array-methods

pint/compat.py

@@ -9,6 +9,7 @@
     :license: BSD, see LICENSE for more details.
 """
 import tokenize
+import warnings
 from io import BytesIO
 from numbers import Number
 from decimal import Decimal
@@ -21,6 +22,30 @@ def tokenizer(input_string):
         yield tokinfo
 
 
+# TODO: remove this warning after v0.10
+class BehaviorChangeWarning(UserWarning):
+    pass
+array_function_change_msg = """The way Pint handles NumPy operations has changed with the
+implementation of NEP 18. Unimplemented NumPy operations will now fail instead of making
+assumptions about units. Some functions, eg concat, will now return Quanties with units, where
+they returned ndarrays previously. See https://github.com/hgrecco/pint/pull/905.
+
+To hide this warning, wrap your first creation of an array Quantity with
+warnings.catch_warnings(), like the following:
+
+import numpy as np
+import warnings
+from pint import Quantity
+
+with warnings.catch_warnings():
+    warnings.simplefilter("ignore")
+    Quantity([])
+
+To disable the new behavior, see
+https://www.numpy.org/neps/nep-0018-array-function-protocol.html#implementation
+"""
+
+
 try:
     import numpy as np
     from numpy import ndarray
@@ -40,6 +65,23 @@ def _to_magnitude(value, force_ndarray=False):
             return np.asarray(value)
         return value
 
+    def _test_array_function_protocol():
+        # Test if the __array_function__ protocol is enabled
+        try:
+            class FakeArray:
+                def __array_function__(self, *args, **kwargs):
+                    return
+
+            np.concatenate([FakeArray()])
+            return True
+        except ValueError:
+            return False
+
+    HAS_NUMPY_ARRAY_FUNCTION = _test_array_function_protocol()
+    SKIP_ARRAY_FUNCTION_CHANGE_WARNING = not HAS_NUMPY_ARRAY_FUNCTION
+
+    NP_NO_VALUE = np._NoValue
+
 except ImportError:
 
     np = None
@@ -50,6 +92,9 @@ class ndarray:
     HAS_NUMPY = False
     NUMPY_VER = '0'
     NUMERIC_TYPES = (Number, Decimal)
+    HAS_NUMPY_ARRAY_FUNCTION = False
+    SKIP_ARRAY_FUNCTION_CHANGE_WARNING = True
+    NP_NO_VALUE = None
 
     def _to_magnitude(value, force_ndarray=False):
         if isinstance(value, (dict, bool)) or value is None:
@@ -77,3 +122,18 @@ def _to_magnitude(value, force_ndarray=False):
 
 if not HAS_BABEL:
     Loc = babel_units = None
+
+# Define location of pint.Quantity in NEP-13 type cast hierarchy by defining upcast and
+# downcast/wrappable types
+def is_upcast_type(other):
+    # Check if class name is in preset list
+    return other.__class__.__name__ in ("PintArray", "Series", "DataArray")
+
+
+def eq(first, second, check_all):
+    """Comparison of scalars and arrays
+    """
+    out = first == second
+    if check_all and isinstance(out, ndarray):
+        return np.all(out)
+    return out