Repository: hgrecco/pint

Issue:
Support for numpy NaN functions: nan{min,max,mean,argmin,argmax,sum}

Discussion:
keewis: thank you for working on this, @jthielen. Don't worry about the slow progress, as long as it is not completely on hold.

Running the xarray tests results in quite a few errors (probably more if we count the xfails). I can't do much right now, but I will try to investigate and report more this weekend.

**re** `inspect`: I already mentioned it in #764, but I strongly support dropping py2 (not my call, though).
Also, I'd say let's move adding tests for wrapping other NEP-18-compatible arrays to another PR so we can keep this as small as possible.
jthielen: Thanks for the initial glace over the PR @keewis! In all my tinkering with this, I had forgotten to run the most recent set of changes past the test suite, which would have caught most of these. Whoops!

I've pushed some temporary fixups, which should get all the tests passing except the unimplemented ones. These changes may or may not be for naught after I redo the implementation logic using `inspect.signature`, which I hope to get to yet this weekend or early this coming week.
keewis: we should also keep in mind that with pint we need to make these changes optional: the tests without numpy installed currently fail because in that case `np` is `None`.
hgrecco: > **re** `inspect`: I already mentioned it in #764, but I strongly support dropping py2 (not my call, though).

See #906 but in principle I agree


hgrecco: As you sure have seen I have link a lot of issues to this PR. I am trying to make a list of numpy related issues that we should check again after this PR. Many of them will not be solved, but others might. The goal is to end 2019 with less open issues and PR!
jthielen: I've finally reached the point where all tests are passing locally with NumPy available, so I think that means @keewis we can move forward with the next phase of testing with xarray!

Also, I've updated the initial comment with the main areas I need to work on before I'd feel ready to remove WIP/draft status.
keewis: I opened a [new issue](https://github.com/pydata/xarray/issues/3594) where I listed the test failures. Most of these are probably issues with the code of xarray, but there are also a few missing function implementations: `nanpercentile`, `copyto`, `einsum`, `all`, `any`, `median`, `prod`, `sum`, `rank`, `isin`. Also `fabs`, `maximum` (and `minimum`?), but these are only used by the tests.
keewis: another method that currently falls back to the magnitude and causes lots of unexpected stripped units is `astype`.
jthielen: For the requested function implementations:

- Have been added:
  - `nanpercentile`
  - `copyto`
  - `isin`
  - `fabs`
  - `maximum`
  - `minimum`
- Should have already been there:
  - `median`
  - `sum`
- Need clarification:
  - `rank`: this is deprecated, so I don't think we should add it to new code?
- Can we leave as xfail for now to not hold up progress here?
  - `einsum`: determining the output unit for an arbitrary combination of sums and products from Einstein summation seems rather nontrivial
  - `all`, `any`: see #866 
  - `prod`: see https://github.com/hgrecco/pint/pull/764#issuecomment-523749008 / #867 

Also, the `astype` and `item` methods have been added.
keewis: Thanks.

> Should have already been there:

yes, that's my fault, I simply checked the test status for the aggregation tests. `sum` and `median` work properly with `int`, but not `float`.

> `rank`: this is deprecated, so I don't think we should add it to new code?

definitely. The warning never shows for me (or gets drowned in the flood of `UnitStrippedWarning`s), though, so I didn't realize.
keewis: somehow
```python
a = pint.Quantity([1, 3, 5], 'm')
a.ndim
```
raises a `AttributeError: NumPy method ndim was not callable.`, but `ndim` is a property (there is, however, a function named `numpy.ndim`)

Patch:
diff --git a/CHANGES b/CHANGES
index 16cc515..87ac0cc 100644
--- a/CHANGES
+++ b/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).
diff --git a/docs/numpy.rst b/docs/numpy.rst
index e5b040b..c4965ee 100644
--- a/docs/numpy.rst
+++ b/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
diff --git a/pint/compat.py b/pint/compat.py
index 6c3221a..abe5126 100644
--- a/pint/compat.py
+++ b/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 @@ try:
             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 @@ except ImportError:
     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 @@ except ImportError:
 
 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
diff --git a/pint/numpy_func.py b/pint/numpy_func.py
new file mode 100644
index 0000000..f7ecf7d
--- /dev/null
+++ b/pint/numpy_func.py
@@ -0,0 +1,610 @@
+# -*- coding: utf-8 -*-
+"""
+    pint.numpy_func
+    ~~~~~~~~~~~~~~~
+
+    :copyright: 2019 by Pint Authors, see AUTHORS for more details.
+    :license: BSD, see LICENSE for more details.
+"""
+
+from inspect import signature
+from itertools import chain
+import warnings
+
+from .compat import NP_NO_VALUE, is_upcast_type, np, eq
+from .errors import DimensionalityError, UnitStrippedWarning
+from .util import iterable, sized
+
+HANDLED_UFUNCS = {}
+HANDLED_FUNCTIONS = {}
+
+
+# Shared Implementation Utilities
+
+def _is_quantity(arg):
+    """Test for _units and _magnitude attrs.
+
+    This is done in place of isinstance(Quantity, arg), which would cause a circular import.
+    """
+    return hasattr(arg, '_units') and hasattr(arg, '_magnitude')
+
+
+def _is_quantity_sequence(arg):
+    """Test for sequences of quantities."""
+    return (iterable(arg) and sized(arg) and not isinstance(arg, str)
+            and all(_is_quantity(item) for item in arg))
+
+
+def _get_first_input_units(args, kwargs={}):
+    """Obtain the first valid unit from a collection of args and kwargs."""
+    for arg in chain(args, kwargs.values()):
+        if _is_quantity(arg):
+            return arg.units
+        elif _is_quantity_sequence(arg):
+            return arg[0].units
+
+
+def convert_arg(arg, pre_calc_units):
+    """Convert quantities and sequences of quantities to pre_calc_units and strip units.
+
+    Helper function for convert_to_consistent_units.
+    """
+    if pre_calc_units is not None:
+        if _is_quantity(arg):
+            return arg.m_as(pre_calc_units)
+        elif _is_quantity_sequence(arg):
+            return [item.m_as(pre_calc_units) for item in arg]
+        elif arg is not None:
+            if pre_calc_units.dimensionless:
+                return pre_calc_units._REGISTRY.Quantity(arg).m_as(pre_calc_units)
+            else:
+                raise DimensionalityError('dimensionless', pre_calc_units)
+    else:
+        if _is_quantity(arg):
+            return arg.m
+        elif _is_quantity_sequence(arg):
+            return [item.m for item in arg]
+    return arg
+
+
+def convert_to_consistent_units(*args, pre_calc_units=None, **kwargs):
+    """Prepare args and kwargs for wrapping by unit conversion and stripping.
+
+    If pre_calc_units is not None, takes the args and kwargs for a NumPy function and converts
+    any Quantity or Sequence of Quantities into the units of the first Quantiy/Sequence of
+    Quantities and returns the magnitudes. Other args/kwargs are treated as dimensionless
+    Quantities. If pre_calc_units is None, units are simply stripped.
+    """
+    return (tuple(convert_arg(arg, pre_calc_units=pre_calc_units) for arg in args),
+            {key: convert_arg(arg, pre_calc_units=pre_calc_units)
+             for key, arg in kwargs.items()})
+
+
+def unwrap_and_wrap_consistent_units(*args):
+    """Strip units from args while providing a rewrapping function.
+
+    Returns the given args as parsed by convert_to_consistent_units assuming units of first
+    arg with units, along with a wrapper to restore that unit to the output.
+    """
+    first_input_units = _get_first_input_units(args)
+    args, _ = convert_to_consistent_units(*args, pre_calc_units=first_input_units)
+    return args, lambda value: first_input_units._REGISTRY.Quantity(value, first_input_units)
+
+
+def get_op_output_unit(unit_op, first_input_units, all_args=[], size=None):
+    """Determine resulting unit from given operation.
+
+    Options for `unit_op`:
+
+    - "sum": `first_input_units`, unless non-multiplicative, which raises
+      OffsetUnitCalculusError
+    - "mul": product of all units in `all_args`
+    - "delta": `first_input_units`, unless non-multiplicative, which uses delta version
+    - "delta,div": like "delta", but divided by all units in `all_args` except the first
+    - "div": unit of first argument in `all_args` (or dimensionless if not a Quantity) divided
+      by all following units
+    - "variance": square of `first_input_units`, unless non-multiplicative, which raises
+      OffsetUnitCalculusError
+    - "square": square of `first_input_units`
+    - "sqrt": square root of `first_input_units`
+    - "reciprocal": reciprocal of `first_input_units`
+    - "size": `first_input_units` raised to the power of `size`
+    """
+    if unit_op == "sum":
+        result_unit = (1 * first_input_units + 1 * first_input_units).units
+    elif unit_op == "mul":
+        product = first_input_units._REGISTRY.parse_units('')
+        for x in all_args:
+            if hasattr(x, 'units'):
+                product *= x.units
+        result_unit = product
+    elif unit_op == "delta":
+        result_unit = (1 * first_input_units - 1 * first_input_units).units
+    elif unit_op == "delta,div":
+        product = (1 * first_input_units - 1 * first_input_units).units
+        for x in all_args[1:]:
+            if hasattr(x, 'units'):
+                product /= x.units
+        result_unit = product
+    elif unit_op == "div":
+        # Start with first arg in numerator, all others in denominator
+        product = getattr(all_args[0], 'units', first_input_units._REGISTRY.parse_units(''))
+        for x in all_args[1:]:
+            if hasattr(x, 'units'):
+                product /= x.units
+        result_unit = product
+    elif unit_op == "variance":
+        result_unit = ((1 * first_input_units + 1 * first_input_units)**2).units
+    elif unit_op == "square":
+        result_unit = first_input_units**2
+    elif unit_op == "sqrt":
+        result_unit = first_input_units**0.5
+    elif unit_op == "reciprocal":
+        result_unit = first_input_units**-1
+    elif unit_op == "size":
+        if size is None:
+            raise ValueError('size argument must be given when unit_op=="size"')
+        result_unit = first_input_units**size
+
+    else:
+        raise ValueError('Output unit method {} not understood'.format(unit_op))
+
+    return result_unit
+
+
+def implements(numpy_func_string, func_type):
+    """Register an __array_function__/__array_ufunc__ implementation for Quantity objects."""
+    def decorator(func):
+        if func_type == 'function':
+            HANDLED_FUNCTIONS[numpy_func_string] = func
+        elif func_type == 'ufunc':
+            HANDLED_UFUNCS[numpy_func_string] = func
+        else:
+            raise ValueError('Invalid func_type {}'.format(func_type))
+        return func
+    return decorator
+
+
+def implement_func(func_type, func_str, input_units=None, output_unit=None):
+    """Add default-behavior NumPy function/ufunc to the handled list."""
+    # If NumPy is not available, do not attempt implement that which does not exist
+    if np is None:
+        return
+
+    func = getattr(np, func_str)
+
+    @implements(func_str, func_type)
+    def implementation(*args, **kwargs):
+        first_input_units = _get_first_input_units(args, kwargs)
+        if input_units == "all_consistent":
+            # Match all input args/kwargs to same units
+            stripped_args, stripped_kwargs = convert_to_consistent_units(
+                *args, pre_calc_units=first_input_units, **kwargs)
+        else:
+            # Match all input args/kwargs to input_units, or if input_units is None, simply
+            # strip units
+            stripped_args, stripped_kwargs = convert_to_consistent_units(
+                *args, pre_calc_units=input_units, **kwargs)
+
+        # Determine result through base numpy function on stripped arguments
+        result_magnitude = func(*stripped_args, **stripped_kwargs)
+
+        if output_unit is None:
+            # Short circuit and return magnitude alone
+            return result_magnitude
+        elif output_unit == "match_input":
+            result_unit = first_input_units
+        elif output_unit in ['sum', 'mul', 'delta', 'delta,div', 'div', 'variance', 'square',
+                             'sqrt', 'reciprocal', 'size']:
+            result_unit = get_op_output_unit(output_unit, first_input_units,
+                                             tuple(chain(args, kwargs.values())))
+        else:
+            result_unit = output_unit
+
+        return first_input_units._REGISTRY.Quantity(result_magnitude, result_unit)
+
+
+"""
+Define ufunc behavior collections.
+
+- `strip_unit_input_output_ufuncs`: units should be ignored on both input and output
+- `matching_input_bare_output_ufuncs`: inputs are converted to matching units, but outputs are
+  returned as-is
+- `matching_input_set_units_output_ufuncs`: inputs are converted to matching units, and the
+  output units are as set by the dict value
+- `set_units_ufuncs`: dict values are specified as (in_unit, out_unit), so that inputs are
+  converted to in_unit before having magnitude passed to NumPy ufunc, and outputs are set to
+  have out_unit
+- `matching_input_copy_units_output_ufuncs`: inputs are converted to matching units, and
+  outputs are set to that unit
+- `copy_units_output_ufuncs`: input units (except the first) are ignored, and output is set to
+  that of the first input unit
+- `op_units_output_ufuncs`: determine output unit from input unit as determined by operation
+  (see `get_op_output_unit`)
+"""
+strip_unit_input_output_ufuncs = ['isnan', 'isinf', 'isfinite', 'signbit']
+matching_input_bare_output_ufuncs = ['equal', 'greater', 'greater_equal', 'less',
+                                     'less_equal', 'not_equal']
+matching_input_set_units_output_ufuncs = {'arctan2': 'radian'}
+set_units_ufuncs = {'cumprod': ('', ''),
+                    'arccos': ('', 'radian'),
+                    'arcsin': ('', 'radian'),
+                    'arctan': ('', 'radian'),
+                    'arccosh': ('', 'radian'),
+                    'arcsinh': ('', 'radian'),
+                    'arctanh': ('', 'radian'),
+                    'exp': ('', ''),
+                    'expm1': ('', ''),
+                    'exp2': ('', ''),
+                    'log': ('', ''),
+                    'log10': ('', ''),
+                    'log1p': ('', ''),
+                    'log2': ('', ''),
+                    'sin': ('radian', ''),
+                    'cos': ('radian', ''),
+                    'tan': ('radian', ''),
+                    'sinh': ('radian', ''),
+                    'cosh': ('radian', ''),
+                    'tanh': ('radian', ''),
+                    'radians': ('degree', 'radian'),
+                    'degrees': ('radian', 'degree'),
+                    'deg2rad': ('degree', 'radian'),
+                    'rad2deg': ('radian', 'degree'),
+                    'logaddexp': ('', ''),
+                    'logaddexp2': ('', '')}
+# TODO (#905 follow-up): while this matches previous behavior, some of these have optional
+# arguments that should not be Quantities. This should be fixed, and tests using these
+# optional arguments should be added.
+matching_input_copy_units_output_ufuncs = ['compress', 'conj', 'conjugate', 'copy',
+                                           'diagonal', 'max', 'mean', 'min',
+                                           'ptp', 'ravel', 'repeat', 'reshape', 'round',
+                                           'squeeze', 'swapaxes', 'take', 'trace',
+                                           'transpose', 'ceil', 'floor', 'hypot', 'rint',
+                                           'copysign', 'nextafter', 'trunc', 'absolute',
+                                           'negative', 'maximum', 'minimum', 'fabs']
+copy_units_output_ufuncs = ['ldexp', 'fmod', 'mod', 'remainder']
+op_units_output_ufuncs = {'var': 'square', 'prod': 'size', 'multiply': 'mul',
+                          'true_divide': 'div', 'divide': 'div', 'floor_divide': 'div',
+                          'sqrt': 'sqrt', 'square': 'square', 'reciprocal': 'reciprocal',
+                          'std': 'sum', 'sum': 'sum', 'cumsum': 'sum', 'matmul': 'mul'}
+
+
+# Perform the standard ufunc implementations based on behavior collections
+
+for ufunc_str in strip_unit_input_output_ufuncs:
+    # Ignore units
+    implement_func('ufunc', ufunc_str, input_units=None, output_unit=None)
+
+for ufunc_str in matching_input_bare_output_ufuncs:
+    # Require all inputs to match units, but output base ndarray
+    implement_func('ufunc', ufunc_str, input_units='all_consistent', output_unit=None)
+
+for ufunc_str, out_unit in matching_input_set_units_output_ufuncs.items():
+    # Require all inputs to match units, but output in specified unit
+    implement_func('ufunc', ufunc_str, input_units='all_consistent', output_unit=out_unit)
+
+for ufunc_str, (in_unit, out_unit) in set_units_ufuncs.items():
+    # Require inputs in specified unit, and output in specified unit
+    implement_func('ufunc', ufunc_str, input_units=in_unit, output_unit=out_unit)
+
+for ufunc_str in matching_input_copy_units_output_ufuncs:
+    # Require all inputs to match units, and output as first unit in arguments
+    implement_func('ufunc', ufunc_str, input_units='all_consistent',
+                   output_unit='match_input')
+
+for ufunc_str in copy_units_output_ufuncs:
+    # Output as first unit in arguments, but do not convert inputs
+    implement_func('ufunc', ufunc_str, input_units=None, output_unit='match_input')
+
+for ufunc_str, unit_op in op_units_output_ufuncs.items():
+    implement_func('ufunc', ufunc_str, input_units=None, output_unit=unit_op)
+
+
+# Define custom ufunc implementations for atypical cases
+
+@implements('modf', 'ufunc')
+def _modf(x, *args, **kwargs):
+    (x,), output_wrap = unwrap_and_wrap_consistent_units(x)
+    return tuple(output_wrap(y) for y in np.modf(x, *args, **kwargs))
+
+
+@implements('frexp', 'ufunc')
+def _frexp(x, *args, **kwargs):
+    (x,), output_wrap = unwrap_and_wrap_consistent_units(x)
+    mantissa, exponent = np.frexp(x, *args, **kwargs)
+    return output_wrap(mantissa), exponent
+
+
+@implements('power', 'ufunc')
+def _power(x1, x2):
+    return x1**x2
+
+
+def _add_subtract_handle_non_quantity_zero(x1, x2):
+    # As in #121/#122, if a value is 0 (but not Quantity 0) do the operation without checking
+    # units. We do the calculation instead of just returning the same value to enforce any
+    # shape checking and type casting due to the operation.
+    if eq(x1, 0, True):
+        (x2,), output_wrap = unwrap_and_wrap_consistent_units(x2)
+    elif eq(x2, 0, True):
+        (x1,), output_wrap = unwrap_and_wrap_consistent_units(x1)
+    else:
+        (x1, x2), output_wrap = unwrap_and_wrap_consistent_units(x1, x2)
+    return x1, x2, output_wrap
+
+
+@implements('add', 'ufunc')
+def _add(x1, x2, *args, **kwargs):
+    x1, x2, output_wrap = _add_subtract_handle_non_quantity_zero(x1, x2)
+    return output_wrap(np.add(x1, x2, *args, **kwargs))
+
+
+@implements('subtract', 'ufunc')
+def _subtract(x1, x2, *args, **kwargs):
+    x1, x2, output_wrap = _add_subtract_handle_non_quantity_zero(x1, x2)
+    return output_wrap(np.subtract(x1, x2, *args, **kwargs))
+
+
+# Define custom function implementations
+
+@implements('meshgrid', 'function')
+def _meshgrid(*xi, **kwargs):
+    # Simply need to map input units to onto list of outputs
+    input_units = (x.units for x in xi)
+    res = np.meshgrid(*(x.m for x in xi), **kwargs)
+    return [out * unit for out, unit in zip(res, input_units)]
+
+
+@implements('full_like', 'function')
+def _full_like(a, fill_value, dtype=None, order='K', subok=True, shape=None):
+    # Make full_like by multiplying with array from ones_like in a
+    # non-multiplicative-unit-safe way
+    if hasattr(fill_value, '_REGISTRY'):
+        return fill_value._REGISTRY.Quantity(
+            (np.ones_like(a, dtype=dtype, order=order, subok=subok, shape=shape)
+             * fill_value.m), fill_value.units)
+    else:
+        return (np.ones_like(a, dtype=dtype, order=order, subok=subok, shape=shape)
+                * fill_value)
+
+
+@implements('interp', 'function')
+def _interp(x, xp, fp, left=None, right=None, period=None):
+    # Need to handle x and y units separately
+    (x, xp, period), _ = unwrap_and_wrap_consistent_units(x, xp, period)
+    (fp, right, left), output_wrap = unwrap_and_wrap_consistent_units(fp, left, right)
+    return output_wrap(np.interp(x, xp, fp, left=left, right=right, period=period))
+
+
+@implements('where', 'function')
+def _where(condition, *args):
+    if (len(args) == 2 and not _is_quantity(args[1]) and not iterable(args[1])
+            and (args[1] == 0 or np.isnan(args[1]))):
+        # Special case for y being bare zero or nan
+        (x,), output_wrap = unwrap_and_wrap_consistent_units(args[0])
+        args = x, args[1]
+    elif (len(args) == 2 and not _is_quantity(args[0]) and not iterable(args[0])
+            and (args[0] == 0 or np.isnan(args[0]))):
+        # Special case for x being bare zero or nan
+        (y,), output_wrap = unwrap_and_wrap_consistent_units(args[1])
+        args = args[0], y
+    else:
+        args, output_wrap = unwrap_and_wrap_consistent_units(*args)
+    return output_wrap(np.where(condition, *args))
+
+
+@implements('concatenate', 'function')
+def _concatenate(sequence, *args, **kwargs):
+    sequence, output_wrap = unwrap_and_wrap_consistent_units(*sequence)
+    return output_wrap(np.concatenate(sequence, *args, **kwargs))
+
+
+@implements('stack', 'function')
+def _stack(arrays, *args, **kwargs):
+    arrays, output_wrap = unwrap_and_wrap_consistent_units(*arrays)
+    return output_wrap(np.stack(arrays, *args, **kwargs))
+
+
+@implements('unwrap', 'function')
+def _unwrap(p, discont=None, axis=-1):
+    # np.unwrap only dispatches over p argument, so assume it is a Quantity
+    discont = np.pi if discont is None else discont
+    return p._REGISTRY.Quantity(np.unwrap(p.m_as('rad'), discont, axis=axis),
+                                'rad').to(p.units)
+
+
+@implements('copyto', 'function')
+def _copyto(dst, src, casting='same_kind', where=True):
+    if _is_quantity(dst):
+        if _is_quantity(src):
+            src = src.m_as(dst.units)
+        np.copyto(dst._magnitude, src, casting=casting, where=where)
+    else:
+        warnings.warn("The unit of the quantity is stripped when copying to non-quantity",
+                      stacklevel=2)
+        np.copyto(dst, src.m, casting=casting, where=where)
+
+
+@implements('einsum', 'function')
+def _einsum(subscripts, *operands, **kwargs):
+    operand_magnitudes, _ = convert_to_consistent_units(*operands, pre_calc_units=None)
+    output_unit = get_op_output_unit("mul", _get_first_input_units(operands), operands)
+    return np.einsum(subscripts, *operand_magnitudes, **kwargs) * output_unit
+
+
+@implements('isin', 'function')
+def _isin(element, test_elements, assume_unique=False, invert=False):
+    if not _is_quantity(element):
+        raise ValueError('Cannot test if unit-aware elements are in not-unit-aware array')
+
+    if _is_quantity(test_elements):
+        try:
+            test_elements = test_elements.m_as(element.units)
+        except DimensionalityError:
+            # Incompatible unit test elements cannot be in element
+            return np.full(element.shape, False)
+    elif _is_quantity_sequence(test_elements):
+        compatible_test_elements = []
+        for test_element in test_elements:
+            try:
+                compatible_test_elements.append(test_element.m_as(element.units))
+            except DimensionalityError:
+                # Incompatible unit test elements cannot be in element, but others in sequence
+                # may
+                pass
+        test_elements = compatible_test_elements
+    else:
+        # Consider non-quantity like dimensionless quantity
+        if not element.dimensionless:
+            # Unit do not match, so all false
+            return np.full(element.shape, False)
+        else:
+            # Convert to units of element
+            element._REGISTRY.Quantity(test_elements).m_as(element.units)
+
+    return np.isin(element.m, test_elements, assume_unique=assume_unique, invert=invert)
+
+
+# Implement simple matching-unit or stripped-unit functions based on signature
+
+def implement_consistent_units_by_argument(func_str, unit_arguments, wrap_output=True):
+    # If NumPy is not available, do not attempt implement that which does not exist
+    if np is None:
+        return
+
+    func = getattr(np, func_str)
+
+    @implements(func_str, 'function')
+    def implementation(*args, **kwargs):
+        # Bind given arguments to the NumPy function signature
+        bound_args = signature(func).bind(*args, **kwargs)
+
+        # Skip unit arguments that are supplied as None
+        valid_unit_arguments = [label for label in unit_arguments
+                                if label in bound_args.arguments
+                                and bound_args.arguments[label] is not None]
+
+        # Unwrap valid unit arguments, ensure consistency, and obtain output wrapper
+        unwrapped_unit_args, output_wrap = unwrap_and_wrap_consistent_units(
+            *(bound_args.arguments[label] for label in valid_unit_arguments))
+
+        # Call NumPy function with updated arguments
+        for i, unwrapped_unit_arg in enumerate(unwrapped_unit_args):
+            bound_args.arguments[valid_unit_arguments[i]] = unwrapped_unit_arg
+        ret = func(*bound_args.args, **bound_args.kwargs)
+
+        # Conditionally wrap output
+        if wrap_output:
+            return output_wrap(ret)
+        else:
+            return ret
+
+
+for func_str, unit_arguments, wrap_output in [('expand_dims', 'a', True),
+                                              ('squeeze', 'a', True),
+                                              ('rollaxis', 'a', True),
+                                              ('moveaxis', 'a', True),
+                                              ('around', 'a', True),
+                                              ('diagonal', 'a', True),
+                                              ('mean', 'a', True),
+                                              ('ptp', 'a', True),
+                                              ('ravel', 'a', True),
+                                              ('round_', 'a', True),
+                                              ('sort', 'a', True),
+                                              ('median', 'a', True),
+                                              ('nanmedian', 'a', True),
+                                              ('transpose', 'a', True),
+                                              ('copy', 'a', True),
+                                              ('average', 'a', True),
+                                              ('nanmean', 'a', True),
+                                              ('swapaxes', 'a', True),
+                                              ('nanmin', 'a', True),
+                                              ('nanmax', 'a', True),
+                                              ('percentile', 'a', True),
+                                              ('nanpercentile', 'a', True),
+                                              ('flip', 'm', True),
+                                              ('fix', 'x', True),
+                                              ('trim_zeros', ['filt'], True),
+                                              ('broadcast_to', ['array'], True),
+                                              ('amax', ['a', 'initial'], True),
+                                              ('amin', ['a', 'initial'], True),
+                                              ('searchsorted', ['a', 'v'], False),
+                                              ('isclose', ['a', 'b', 'rtol', 'atol'], False),
+                                              ('nan_to_num', ['x', 'nan', 'posinf', 'neginf'],
+                                               True),
+                                              ('clip', ['a', 'a_min', 'a_max'], True),
+                                              ('append', ['arr', 'values'], True),
+                                              ('compress', 'a', True),
+                                              ('linspace', ['start', 'stop'], True),
+                                              ('tile', 'A', True),
+                                              ('rot90', 'm', True),
+                                              ('insert', ['arr', 'values'], True)]:
+    implement_consistent_units_by_argument(func_str, unit_arguments, wrap_output)
+
+
+# Handle atleast_nd functions
+
+def implement_atleast_nd(func_str):
+    # If NumPy is not available, do not attempt implement that which does not exist
+    if np is None:
+        return
+
+    func = getattr(np, func_str)
+    @implements(func_str, 'function')
+    def implementation(*arrays):
+        stripped_arrays, _ = convert_to_consistent_units(*arrays)
+        arrays_magnitude = func(*stripped_arrays)
+        if len(arrays) > 1:
+            return [array_magnitude if not hasattr(original, '_REGISTRY')
+                    else original._REGISTRY.Quantity(array_magnitude, original.units)
+                    for array_magnitude, original in zip(arrays_magnitude, arrays)]
+        else:
+            output_unit = arrays[0].units
+            return output_unit._REGISTRY.Quantity(arrays_magnitude, output_unit)
+
+
+for func_str in ['atleast_1d', 'atleast_2d', 'atleast_3d']:
+    implement_atleast_nd(func_str)
+
+# Handle single-argument consistent unit functions
+for func_str in ['block', 'hstack', 'vstack', 'dstack', 'column_stack']:
+    implement_func('function', func_str, input_units='all_consistent',
+                   output_unit='match_input')
+
+# Handle cumulative products (which must be dimensionless for consistent units across output
+# array)
+for func_str in ['cumprod', 'cumproduct', 'nancumprod']:
+    implement_func('function', func_str, input_units='dimensionless',
+                   output_unit='match_input')
+
+# Handle functions that ignore units on input and output
+for func_str in ['size', 'isreal', 'iscomplex', 'shape', 'ones_like', 'zeros_like',
+                 'empty_like', 'argsort', 'argmin', 'argmax', 'alen', 'ndim', 'nanargmax',
+                 'nanargmin', 'count_nonzero', 'nonzero', 'result_type']:
+    implement_func('function', func_str, input_units=None, output_unit=None)
+
+# Handle functions with output unit defined by operation
+for func_str in ['std', 'nanstd', 'sum', 'nansum', 'cumsum', 'nancumsum']:
+    implement_func('function', func_str, input_units=None, output_unit='sum')
+for func_str in ['cross', 'trapz', 'dot']:
+    implement_func('function', func_str, input_units=None, output_unit='mul')
+for func_str in ['diff', 'ediff1d']:
+    implement_func('function', func_str, input_units=None, output_unit='delta')
+for func_str in ['gradient', ]:
+    implement_func('function', func_str, input_units=None, output_unit='delta,div')
+for func_str in ['var', 'nanvar']:
+    implement_func('function', func_str, input_units=None, output_unit='variance')
+
+
+def numpy_wrap(func_type, func, args, kwargs, types):
+    """Return the result from a NumPy function/ufunc as wrapped by Pint."""
+    if func_type == 'function':
+        handled = HANDLED_FUNCTIONS
+    elif func_type == 'ufunc':
+        handled = HANDLED_UFUNCS
+    else:
+        raise ValueError('Invalid func_type {}'.format(func_type))
+
+    if func.__name__ not in handled or any(is_upcast_type(t) for t in types):
+        return NotImplemented
+    return handled[func.__name__](*args, **kwargs)
diff --git a/pint/quantity.py b/pint/quantity.py
index f0bdd5f..6d57e42 100644
--- a/pint/quantity.py
+++ b/pint/quantity.py
@@ -19,24 +19,23 @@ import operator
 import re
 import warnings
 
+from distutils.version import StrictVersion
+
 from .formatting import (remove_custom_flags, siunitx_format_unit, ndarray_to_latex,
                          ndarray_to_latex_parts)
 from .errors import (DimensionalityError, OffsetUnitCalculusError, PintTypeError,
                      UndefinedUnitError, UnitStrippedWarning)
 from .definitions import UnitDefinition
-from .compat import ndarray, np, _to_magnitude
+from .compat import (Loc, NUMPY_VER, SKIP_ARRAY_FUNCTION_CHANGE_WARNING,
+                     BehaviorChangeWarning, ndarray, np, _to_magnitude, is_upcast_type, eq,
+                     array_function_change_msg)
 from .util import (PrettyIPython, logger, UnitsContainer, SharedRegistryObject,
-                   to_units_container, infer_base_unit)
-from pint.compat import Loc
-
-
-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
+                   to_units_container, infer_base_unit, iterable, sized)
+from .numpy_func import (HANDLED_UFUNCS, copy_units_output_ufuncs, get_op_output_unit,
+                         matching_input_bare_output_ufuncs,
+                         matching_input_copy_units_output_ufuncs, 
+                         matching_input_set_units_output_ufuncs, numpy_wrap, 
+                         op_units_output_ufuncs, set_units_ufuncs)
 
 
 class _Exception(Exception):            # pragma: no cover
@@ -67,7 +66,7 @@ def ireduce_dimensions(f):
 def check_implemented(f):
     def wrapped(self, *args, **kwargs):
         other=args[0]
-        if other.__class__.__name__ in ["PintArray", "Series"]:
+        if is_upcast_type(other):
             return NotImplemented
         # pandas often gets to arrays of quantities [ Q_(1,"m"), Q_(2,"m")]
         # and expects Quantity * array[Quantity] should return NotImplemented
@@ -121,6 +120,8 @@ class Quantity(PrettyIPython, SharedRegistryObject):
         return _unpickle, (Quantity, self.magnitude, self._units)
 
     def __new__(cls, value, units=None):
+        global SKIP_ARRAY_FUNCTION_CHANGE_WARNING
+
         if units is None:
             if isinstance(value, str):
                 if value == '':
@@ -159,6 +160,10 @@ class Quantity(PrettyIPython, SharedRegistryObject):
         inst.__used = False
         inst.__handling = None
 
+        if not SKIP_ARRAY_FUNCTION_CHANGE_WARNING and isinstance(inst._magnitude, ndarray):
+            warnings.warn(array_function_change_msg, BehaviorChangeWarning)
+            SKIP_ARRAY_FUNCTION_CHANGE_WARNING = True
+
         return inst
 
     @property
@@ -629,8 +634,7 @@ class Quantity(PrettyIPython, SharedRegistryObject):
                 raise
             except TypeError:
                 return NotImplemented
-
-            if _eq(other, 0, True):
+            if eq(other, 0, True):
                 # If the other value is 0 (but not Quantity 0)
                 # do the operation without checking units.
                 # We do the calculation instead of just returning the same
@@ -725,7 +729,7 @@ class Quantity(PrettyIPython, SharedRegistryObject):
         """
         if not self._check(other):
             # other not from same Registry or not a Quantity
-            if _eq(other, 0, True):
+            if eq(other, 0, True):
                 # If the other value is 0 (but not Quantity 0)
                 # do the operation without checking units.
                 # We do the calculation instead of just returning the same
@@ -982,6 +986,15 @@ class Quantity(PrettyIPython, SharedRegistryObject):
 
     __rmul__ = __mul__
 
+    def __matmul__(self, other):
+        # Use NumPy ufunc (existing since 1.16) for matrix multiplication
+        if StrictVersion(NUMPY_VER) >= StrictVersion('1.16'):
+            return np.matmul(self, other)
+        else:
+            return NotImplemented
+
+    __rmatmul__ = __matmul__
+
     def __itruediv__(self, other):
         if not isinstance(self._magnitude, ndarray):
             return self._mul_div(other, operator.truediv)
@@ -1228,31 +1241,31 @@ class Quantity(PrettyIPython, SharedRegistryObject):
         # We compare to the base class of Quantity because
         # each Quantity class is unique.
         if not isinstance(other, Quantity):
-            if _eq(other, 0, True):
+            if eq(other, 0, True):
                 # Handle the special case in which we compare to zero
                 # (or an array of zeros)
                 if self._is_multiplicative:
                     # compare magnitude
-                    return _eq(self._magnitude, other, False)
+                    return eq(self._magnitude, other, False)
                 else:
                     # compare the magnitude after converting the
                     # non-multiplicative quantity to base units
                     if self._REGISTRY.autoconvert_offset_to_baseunit:
-                        return _eq(self.to_base_units()._magnitude, other, False)
+                        return eq(self.to_base_units()._magnitude, other, False)
                     else:
                         raise OffsetUnitCalculusError(self._units)
 
             return (self.dimensionless and
-                    _eq(self._convert_magnitude(UnitsContainer()), other, False))
+                    eq(self._convert_magnitude(UnitsContainer()), other, False))
 
-        if _eq(self._magnitude, 0, True) and _eq(other._magnitude, 0, True):
+        if eq(self._magnitude, 0, True) and eq(other._magnitude, 0, True):
             return self.dimensionality == other.dimensionality
 
         if self._units == other._units:
-            return _eq(self._magnitude, other._magnitude, False)
+            return eq(self._magnitude, other._magnitude, False)
 
         try:
-            return _eq(self._convert_magnitude_not_inplace(other._units),
+            return eq(self._convert_magnitude_not_inplace(other._units),
                        other._magnitude, False)
         except DimensionalityError:
             return False
@@ -1268,7 +1281,7 @@ class Quantity(PrettyIPython, SharedRegistryObject):
         if not isinstance(other, self.__class__):
             if self.dimensionless:
                 return op(self._convert_magnitude_not_inplace(UnitsContainer()), other)
-            elif _eq(other, 0, True):
+            elif eq(other, 0, True):
                 # Handle the special case in which we compare to zero
                 # (or an array of zeros)
                 if self._is_multiplicative:
@@ -1302,62 +1315,56 @@ class Quantity(PrettyIPython, SharedRegistryObject):
 
     __nonzero__ = __bool__
 
-    # NumPy Support
-    __radian = 'radian'
-    __same_units = 'equal greater greater_equal less less_equal not_equal arctan2'.split()
-    #: Dictionary mapping ufunc/attributes names to the units that they
-    #: require (conversion will be tried).
-    __require_units = {'cumprod': '',
-                       'arccos': '', 'arcsin': '', 'arctan': '',
-                       'arccosh': '', 'arcsinh': '', 'arctanh': '',
-                       'exp': '', 'expm1': '', 'exp2': '',
-                       'log': '', 'log10': '', 'log1p': '', 'log2': '',
-                       'sin': __radian, 'cos': __radian, 'tan': __radian,
-                       'sinh': __radian, 'cosh': __radian, 'tanh': __radian,
-                       'radians': 'degree', 'degrees': __radian,
-                       'deg2rad': 'degree', 'rad2deg': __radian,
-                       'logaddexp': '', 'logaddexp2': ''}
-
-    #: Dictionary mapping ufunc/attributes names to the units that they
-    #: will set on output.
-    __set_units = {'cos': '', 'sin': '', 'tan': '',
-                   'cosh': '', 'sinh': '', 'tanh': '',
-                   'log': '', 'exp': '',
-                   'arccos': __radian, 'arcsin': __radian,
-                   'arctan': __radian, 'arctan2': __radian,
-                   'arccosh': __radian, 'arcsinh': __radian,
-                   'arctanh': __radian,
-                   'degrees': 'degree', 'radians': __radian,
-                   'expm1': '', 'cumprod': '',
-                   'rad2deg': 'degree', 'deg2rad': __radian}
-
-    #: List of ufunc/attributes names in which units are copied from the
-    #: original.
-    __copy_units = 'compress conj conjugate copy cumsum diagonal flatten ' \
-                   'max mean min ptp ravel repeat reshape round ' \
-                   'squeeze std sum swapaxes take trace transpose ' \
-                   'ceil floor hypot rint ' \
-                   'add subtract ' \
-                   'copysign nextafter trunc ' \
-                   'frexp ldexp modf modf__1 ' \
-                   'absolute negative remainder fmod mod'.split()
-
-    #: Dictionary mapping ufunc/attributes names to the units that they will
-    #: set on output. The value is interpreted as the power to which the unit
-    #: will be raised.
-    __prod_units = {'var': 2, 'prod': 'size', 'multiply': 'mul',
-                    'true_divide': 'div', 'divide': 'div', 'floor_divide': 'div',
-                    'remainder': 'div',
-                    'sqrt': .5, 'square': 2, 'reciprocal': -1}
-
-    __skip_other_args = 'ldexp multiply ' \
-                        'true_divide divide floor_divide fmod mod ' \
-                        'remainder'.split()
-
-    __handled = tuple(__same_units) + \
-                tuple(__require_units.keys()) + \
-                tuple(__prod_units.keys()) + \
-                tuple(__copy_units) + tuple(__skip_other_args)
+    # NumPy function/ufunc support
+    __array_priority__ = 17
+
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        if method != "__call__":
+            # Only handle ufuncs as callables
+            return NotImplemented
+
+        # Replicate types from __array_function__
+        types = set(type(arg) for arg in list(inputs) + list(kwargs.values())
+                    if hasattr(arg, "__array_ufunc__"))
+
+        return numpy_wrap('ufunc', ufunc, inputs, kwargs, types)
+
+    def __array_function__(self, func, types, args, kwargs):
+        return numpy_wrap('function', func, args, kwargs, types)
+
+    _wrapped_numpy_methods = ['flatten', 'astype', 'item']
+
+    def _numpy_method_wrap(self, func, *args, **kwargs):
+        """Convenience method to wrap on the fly NumPy ndarray methods taking
+        care of the units.
+        """
+        # Set input units if needed
+        if func.__name__ in set_units_ufuncs:
+            self.__ito_if_needed(set_units_ufuncs[func.__name__][0])
+
+        value = func(*args, **kwargs)
+
+        # Set output units as needed
+        if (func.__name__ in
+                (matching_input_copy_units_output_ufuncs + copy_units_output_ufuncs
+                 + self._wrapped_numpy_methods)):
+            output_unit = self._units
+        elif func.__name__ in set_units_ufuncs:
+            output_unit = set_units_ufuncs[func.__name__][1]
+        elif func.__name__ in matching_input_set_units_output_ufuncs:
+            output_unit = matching_input_set_units_output_ufuncs[func.__name__]
+        elif func.__name__ in op_units_output_ufuncs:
+            output_unit = get_op_output_unit(op_units_output_ufuncs[func.__name__],
+                                             self.units,
+                                             list(args) + list(kwargs.values()),
+                                             self._magnitude.size)
+        else:
+            output_unit = None
+
+        if output_unit is not None:
+            return self.__class__(value, output_unit)
+        else:
+            return value
 
     def clip(self, first=None, second=None, out=None, **kwargs):
         min = kwargs.get('min', first)
@@ -1436,32 +1443,19 @@ class Quantity(PrettyIPython, SharedRegistryObject):
             raise DimensionalityError('dimensionless', self._units)
         return self.magnitude.searchsorted(v, side)
 
+    def dot(self, b):
+        """Dot product of two arrays.
+
+        Wraps np.dot().
+        """
+        return np.dot(self, b)
+
     def __ito_if_needed(self, to_units):
         if self.unitless and to_units == 'radian':
             return
 
         self.ito(to_units)
 
-    def __numpy_method_wrap(self, func, *args, **kwargs):
-        """Convenience method to wrap on the fly numpy method taking
-        care of the units.
-        """
-        if func.__name__ in self.__require_units:
-            self.__ito_if_needed(self.__require_units[func.__name__])
-
-        value = func(*args, **kwargs)
-
-        if func.__name__ in self.__copy_units:
-            return self.__class__(value, self._units)
-
-        if func.__name__ in self.__prod_units:
-            tmp = self.__prod_units[func.__name__]
-            if tmp == 'size':
-                return self.__class__(value, self._units ** self._magnitude.size)
-            return self.__class__(value, self._units ** tmp)
-
-        return value
-
     def __len__(self):
         return len(self._magnitude)
 
@@ -1469,24 +1463,28 @@ class Quantity(PrettyIPython, SharedRegistryObject):
         # Attributes starting with `__array_` are common attributes of NumPy ndarray.
         # They are requested by numpy functions.
         if item.startswith('__array_'):
-            warnings.warn("The unit of the quantity is stripped.", UnitStrippedWarning, stacklevel=2)
+            warnings.warn("The unit of the quantity is stripped when getting {} "
+                          "attribute".format(item), UnitStrippedWarning, stacklevel=2)
             if isinstance(self._magnitude, ndarray):
                 return getattr(self._magnitude, item)
             else:
                 # If an `__array_` attributes is requested but the magnitude is not an ndarray,
                 # we convert the magnitude to a numpy ndarray.
-                self._magnitude = _to_magnitude(self._magnitude, force_ndarray=True)
-                return getattr(self._magnitude, item)
-        elif item in self.__handled:
-            if not isinstance(self._magnitude, ndarray):
-                self._magnitude = _to_magnitude(self._magnitude, True)
-            attr = getattr(self._magnitude, item)
+                # TODO (#905 follow-up): Potentially problematic, investigate for duck arrays
+                magnitude_as_array = _to_magnitude(self._magnitude, force_ndarray=True)
+                return getattr(magnitude_as_array, item)
+        elif item in HANDLED_UFUNCS or item in self._wrapped_numpy_methods:
+            # TODO (#905 follow-up): Potentially problematic, investigate for duck arrays/scalars
+            magnitude_as_array = _to_magnitude(self._magnitude, True)
+            attr = getattr(magnitude_as_array, item)
             if callable(attr):
-                return functools.partial(self.__numpy_method_wrap, attr)
-            return attr
+                return functools.partial(self._numpy_method_wrap, attr)
+            else:
+                raise AttributeError('NumPy method {} was not callable.'.format(item))
+
         try:
             return getattr(self._magnitude, item)
-        except AttributeError as ex:
+        except AttributeError:
             raise AttributeError("Neither Quantity object nor its magnitude ({}) "
                                  "has attribute '{}'".format(self._magnitude, item))
 
@@ -1539,172 +1537,6 @@ class Quantity(PrettyIPython, SharedRegistryObject):
         return [self.__class__(value, units).tolist() if isinstance(value, list) else self.__class__(value, units)
                 for value in self._magnitude.tolist()]
 
-    __array_priority__ = 17
-
-    def _call_ufunc(self, ufunc, *inputs, **kwargs):
-        # Store the destination units
-        dst_units = None
-        # List of magnitudes of Quantities with the right units
-        # to be used as argument of the ufunc
-        mobjs = None
-
-        if ufunc.__name__ in self.__require_units:
-            # ufuncs in __require_units
-            # require specific units
-            # This is more complex that it should be due to automatic
-            # conversion between radians/dimensionless
-            # TODO: maybe could be simplified using Contexts
-            dst_units = self.__require_units[ufunc.__name__]
-            if dst_units == 'radian':
-                mobjs = []
-                for other in inputs:
-                    unt = getattr(other, '_units', '')
-                    if unt == 'radian':
-                        mobjs.append(getattr(other, 'magnitude', other))
-                    else:
-                        factor, units = self._REGISTRY._get_root_units(unt)
-                        if units and units != UnitsContainer({'radian': 1}):
-                            raise DimensionalityError(units, dst_units)
-                        mobjs.append(getattr(other, 'magnitude', other) * factor)
-                mobjs = tuple(mobjs)
-            else:
-                dst_units = self._REGISTRY.parse_expression(dst_units)._units
-
-        elif len(inputs) > 1 and ufunc.__name__ not in self.__skip_other_args:
-            # ufunc with multiple arguments require that all inputs have
-            # the same arguments unless they are in __skip_other_args
-            dst_units = getattr(inputs[0], "_units", None)
-
-        # Do the conversion (if needed) and extract the magnitude for each input.
-        if mobjs is None:
-            if dst_units is not None:
-                mobjs = tuple(self._REGISTRY.convert(getattr(other, 'magnitude', other),
-                                                     getattr(other, 'units', ''),
-                                                     dst_units)
-                              for other in inputs)
-            else:
-                mobjs = tuple(getattr(other, 'magnitude', other)
-                              for other in inputs)
-
-        # call the ufunc
-        try:
-            return ufunc(*mobjs)
-        except Exception as ex:
-            raise _Exception(ex)
-
-
-    def _wrap_output(self, ufname, i, objs, out):
-        """we set the units of the output value"""
-        if i > 0:
-            ufname = "{}__{}".format(ufname, i)
-
-        if ufname in self.__set_units:
-            try:
-                out = self.__class__(out, self.__set_units[ufname])
-            except:
-                raise _Exception(ValueError)
-        elif ufname in self.__copy_units:
-            try:
-                out = self.__class__(out, self._units)
-            except:
-                raise _Exception(ValueError)
-        elif ufname in self.__prod_units:
-            tmp = self.__prod_units[ufname]
-            if tmp == 'size':
-                out = self.__class__(out, self._units ** self._magnitude.size)
-            elif tmp == 'div':
-                units1 = objs[0]._units if isinstance(objs[0], self.__class__) else UnitsContainer()
-                units2 = objs[1]._units if isinstance(objs[1], self.__class__) else UnitsContainer()
-                out = self.__class__(out, units1 / units2)
-            elif tmp == 'mul':
-                units1 = objs[0]._units if isinstance(objs[0], self.__class__) else UnitsContainer()
-                units2 = objs[1]._units if isinstance(objs[1], self.__class__) else UnitsContainer()
-                out = self.__class__(out, units1 * units2)
-            else:
-                out = self.__class__(out, self._units ** tmp)
-
-        return out
-
-
-    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
-        if method != "__call__":
-            return NotImplemented
-
-        try:
-            out = self._call_ufunc(ufunc, *inputs, **kwargs)
-            if isinstance(out, tuple):
-                ret = tuple(self._wrap_output(ufunc.__name__, i, inputs, o)
-                            for i, o in enumerate(out))
-                return ret
-            else:
-                return self._wrap_output(ufunc.__name__, 0, inputs, out)
-        except (DimensionalityError, UndefinedUnitError):
-            raise
-        except _Exception as ex:
-            raise ex.internal
-        except:
-            return NotImplemented
-
-
-    def __array_prepare__(self, obj, context=None):
-        # If this uf is handled by Pint, write it down in the handling dictionary.
-
-        # name of the ufunc, argument of the ufunc, domain of the ufunc
-        # In ufuncs with multiple outputs, domain indicates which output
-        # is currently being prepared (eg. see modf).
-        # In ufuncs with a single output, domain is 0
-        uf, objs, i_out = context
-
-        if uf.__name__ in self.__handled and i_out == 0:
-            # Only one ufunc should be handled at a time.
-            # If a ufunc is already being handled (and this is not another domain),
-            # something is wrong..
-            if self.__handling:
-                raise Exception('Cannot handled nested ufuncs.\n'
-                                'Current: {}\n'
-                                'New: {}'.format(context, self.__handling))
-            self.__handling = context
-
-        return obj
-
-    def __array_wrap__(self, obj, context=None):
-        uf, objs, i_out = context
-
-        # if this ufunc is not handled by Pint, pass it to the magnitude.
-        if uf.__name__ not in self.__handled:
-            return self.magnitude.__array_wrap__(obj, context)
-
-        try:
-            # First, we check the units of the input arguments.
-
-            if i_out == 0:
-                out = self._call_ufunc(uf, *objs)
-                # If there are multiple outputs,
-                # store them in __handling (uf, objs, i_out, out0, out1, ...)
-                # and return the first
-                if uf.nout > 1:
-                    self.__handling += out
-                    out = out[0]
-            else:
-                # If this is not the first output,
-                # just grab the result that was previously calculated.
-                out = self.__handling[3 + i_out]
-
-            return self._wrap_output(uf.__name__, i_out, objs, out)
-        except (DimensionalityError, UndefinedUnitError) as ex:
-            raise ex
-        except _Exception as ex:
-            raise ex.internal
-        except Exception as ex:
-            print(ex)
-        finally:
-            # If this is the last output argument for the ufunc,
-            # we are done handling this ufunc.
-            if uf.nout == i_out + 1:
-                self.__handling = None
-
-        return self.magnitude.__array_wrap__(obj, context)
-
     # Measurement support
     def plus_minus(self, error, relative=False):
         if isinstance(error, self.__class__):


Test patch:
diff --git a/pint/testsuite/helpers.py b/pint/testsuite/helpers.py
index cb6d610..cbea4c0 100644
--- a/pint/testsuite/helpers.py
+++ b/pint/testsuite/helpers.py
@@ -5,7 +5,19 @@ from distutils.version import StrictVersion
 import re
 import unittest
 
-from ..compat import HAS_NUMPY, HAS_BABEL, HAS_UNCERTAINTIES, NUMPY_VER
+from ..compat import HAS_NUMPY, HAS_BABEL, HAS_UNCERTAINTIES, HAS_NUMPY_ARRAY_FUNCTION, NUMPY_VER
+
+
+def requires_array_function_protocol():
+    if not HAS_NUMPY:
+        return unittest.skip('Requires NumPy')
+    return unittest.skipUnless(HAS_NUMPY_ARRAY_FUNCTION, 'Requires __array_function__ protocol to be enabled')
+
+
+def requires_not_array_function_protocol():
+    if not HAS_NUMPY:
+        return unittest.skip('Requires NumPy')
+    return unittest.skipIf(HAS_NUMPY_ARRAY_FUNCTION, 'Requires __array_function__ protocol to be unavailable or disabled')
 
 
 def requires_numpy18():
@@ -20,6 +32,12 @@ def requires_numpy_previous_than(version):
     return unittest.skipUnless(StrictVersion(NUMPY_VER) < StrictVersion(version), 'Requires NumPy < %s' % version)
 
 
+def requires_numpy_at_least(version):
+    if not HAS_NUMPY:
+        return unittest.skip('Requires NumPy')
+    return unittest.skipUnless(StrictVersion(NUMPY_VER) >= StrictVersion(version), 'Requires NumPy >= %s' % version)
+
+
 def requires_numpy():
     return unittest.skipUnless(HAS_NUMPY, 'Requires NumPy')
 
diff --git a/pint/testsuite/test_numpy.py b/pint/testsuite/test_numpy.py
index 1697ba8..ead8895 100644
--- a/pint/testsuite/test_numpy.py
+++ b/pint/testsuite/test_numpy.py
@@ -4,7 +4,7 @@ import copy
 import operator as op
 import unittest
 
-from pint import DimensionalityError
+from pint import DimensionalityError, OffsetUnitCalculusError, set_application_registry
 from pint.compat import np
 from pint.testsuite import QuantityTestCase, helpers
 from pint.testsuite.test_umath import TestUFuncs
@@ -25,26 +25,53 @@ class TestNumpyMethods(QuantityTestCase):
     def q(self):
         return [[1,2],[3,4]] * self.ureg.m
 
-    def test_tolist(self):
-        self.assertEqual(self.q.tolist(), [[1*self.ureg.m, 2*self.ureg.m], [3*self.ureg.m, 4*self.ureg.m]])
+    @property
+    def q_nan(self):
+        return [[1,2],[3,np.nan]] * self.ureg.m
 
-    def test_sum(self):
-        self.assertEqual(self.q.sum(), 10*self.ureg.m)
-        self.assertQuantityEqual(self.q.sum(0), [4,     6]*self.ureg.m)
-        self.assertQuantityEqual(self.q.sum(1), [3, 7]*self.ureg.m)
+    @property
+    def q_temperature(self):
+        return self.Q_([[1,2],[3,4]], self.ureg.degC)
 
-    def test_fill(self):
-        tmp = self.q
-        tmp.fill(6 * self.ureg.ft)
-        self.assertQuantityEqual(tmp, [[6, 6], [6, 6]] * self.ureg.ft)
-        tmp.fill(5 * self.ureg.m)
-        self.assertQuantityEqual(tmp, [[5, 5], [5, 5]] * self.ureg.m)
+    def assertNDArrayEqual(self, actual, desired):
+        # Assert that the given arrays are equal, and are not Quantities
+        np.testing.assert_array_equal(actual, desired)
+        self.assertFalse(isinstance(actual, self.Q_))
+        self.assertFalse(isinstance(desired, self.Q_))
 
-    def test_reshape(self):
-        self.assertQuantityEqual(self.q.reshape([1,4]), [[1, 2, 3, 4]] * self.ureg.m)
 
-    def test_transpose(self):
-        self.assertQuantityEqual(self.q.transpose(), [[1, 3], [2, 4]] * self.ureg.m)
+class TestNumpyArrayCreation(TestNumpyMethods):
+    # https://docs.scipy.org/doc/numpy/reference/routines.array-creation.html
+
+    @helpers.requires_array_function_protocol()
+    def test_ones_like(self):
+        self.assertNDArrayEqual(np.ones_like(self.q), np.array([[1, 1], [1, 1]]))
+
+    @helpers.requires_array_function_protocol()
+    def test_zeros_like(self):
+        self.assertNDArrayEqual(np.zeros_like(self.q), np.array([[0, 0], [0, 0]]))
+
+    @helpers.requires_array_function_protocol()
+    def test_empty_like(self):
+        ret = np.empty_like(self.q)
+        self.assertEqual(ret.shape, (2, 2))
+        self.assertTrue(isinstance(ret, np.ndarray))
+
+    @helpers.requires_array_function_protocol()
+    def test_full_like(self):
+        self.assertQuantityEqual(np.full_like(self.q, self.Q_(0, self.ureg.degC)),
+                                 self.Q_([[0, 0], [0, 0]], self.ureg.degC))
+        self.assertNDArrayEqual(np.full_like(self.q, 2), np.array([[2, 2], [2, 2]]))
+
+
+class TestNumpyArrayManipulation(TestNumpyMethods):
+    # TODO
+    # https://www.numpy.org/devdocs/reference/routines.array-manipulation.html
+    # copyto
+    # broadcast , broadcast_arrays
+    # asarray	asanyarray	asmatrix	asfarray	asfortranarray	ascontiguousarray	asarray_chkfinite	asscalar	require
+
+    # Changing array shape
 
     def test_flatten(self):
         self.assertQuantityEqual(self.q.flatten(), [1, 2, 3, 4] * self.ureg.m)
@@ -53,15 +80,306 @@ class TestNumpyMethods(QuantityTestCase):
         for q, v in zip(self.q.flat, [1, 2, 3, 4]):
             self.assertEqual(q, v * self.ureg.m)
 
+    def test_reshape(self):
+        self.assertQuantityEqual(self.q.reshape([1,4]), [[1, 2, 3, 4]] * self.ureg.m)
+
     def test_ravel(self):
         self.assertQuantityEqual(self.q.ravel(), [1, 2, 3, 4] * self.ureg.m)
 
+    @helpers.requires_array_function_protocol()
+    def test_ravel_numpy_func(self):
+        self.assertQuantityEqual(np.ravel(self.q), [1, 2, 3, 4] * self.ureg.m)
+
+    # Transpose-like operations
+
+    @helpers.requires_array_function_protocol()
+    def test_moveaxis(self):
+        self.assertQuantityEqual(np.moveaxis(self.q, 1,0), np.array([[1,2],[3,4]]).T * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_rollaxis(self):
+        self.assertQuantityEqual(np.rollaxis(self.q, 1), np.array([[1,2],[3,4]]).T * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_swapaxes(self):
+        self.assertQuantityEqual(np.swapaxes(self.q, 1,0), np.array([[1,2],[3,4]]).T * self.ureg.m)
+
+    def test_transpose(self):
+        self.assertQuantityEqual(self.q.transpose(), [[1, 3], [2, 4]] * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_transpose_numpy_func(self):
+        self.assertQuantityEqual(np.transpose(self.q), [[1, 3], [2, 4]] * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_flip_numpy_func(self):
+        self.assertQuantityEqual(np.flip(self.q, axis=0), [[3, 4], [1, 2]] * self.ureg.m)
+
+    # Changing number of dimensions
+
+    @helpers.requires_array_function_protocol()
+    def test_atleast_1d(self):
+        actual = np.atleast_1d(self.Q_(0, self.ureg.degC), self.q.flatten())
+        expected = (self.Q_(np.array([0]), self.ureg.degC), self.q.flatten())
+        for ind_actual, ind_expected in zip(actual, expected):
+            self.assertQuantityEqual(ind_actual, ind_expected)
+        self.assertQuantityEqual(np.atleast_1d(self.q), self.q)
+
+    @helpers.requires_array_function_protocol()
+    def test_atleast_2d(self):
+        actual = np.atleast_2d(self.Q_(0, self.ureg.degC), self.q.flatten())
+        expected = (self.Q_(np.array([[0]]), self.ureg.degC), np.array([[1, 2, 3, 4]]) * self.ureg.m)
+        for ind_actual, ind_expected in zip(actual, expected):
+            self.assertQuantityEqual(ind_actual, ind_expected)
+        self.assertQuantityEqual(np.atleast_2d(self.q), self.q)
+
+    @helpers.requires_array_function_protocol()
+    def test_atleast_3d(self):
+        actual = np.atleast_3d(self.Q_(0, self.ureg.degC), self.q.flatten())
+        expected = (self.Q_(np.array([[[0]]]), self.ureg.degC), np.array([[[1], [2], [3], [4]]]) * self.ureg.m)
+        for ind_actual, ind_expected in zip(actual, expected):
+            self.assertQuantityEqual(ind_actual, ind_expected)
+        self.assertQuantityEqual(np.atleast_3d(self.q), np.array([[[1],[2]],[[3],[4]]])* self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_broadcast_to(self):
+        self.assertQuantityEqual(np.broadcast_to(self.q[:,1], (2,2)), np.array([[2,4],[2,4]]) * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_expand_dims(self):
+        self.assertQuantityEqual(np.expand_dims(self.q, 0), np.array([[[1, 2],[3, 4]]])* self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
     def test_squeeze(self):
+        self.assertQuantityEqual(np.squeeze(self.q), self.q)
         self.assertQuantityEqual(
             self.q.reshape([1,4]).squeeze(),
             [1, 2, 3, 4] * self.ureg.m
         )
 
+    # Changing number of dimensions
+    # Joining arrays
+    @helpers.requires_array_function_protocol()
+    def test_concatentate(self):
+        self.assertQuantityEqual(
+            np.concatenate([self.q]*2),
+            self.Q_(np.concatenate([self.q.m]*2), self.ureg.m)
+        )
+
+    @helpers.requires_array_function_protocol()
+    def test_stack(self):
+        self.assertQuantityEqual(
+            np.stack([self.q]*2),
+            self.Q_(np.stack([self.q.m]*2), self.ureg.m)
+        )
+
+    @helpers.requires_array_function_protocol()
+    def test_column_stack(self):
+        self.assertQuantityEqual(
+            np.column_stack([self.q[:,0],self.q[:,1]]),
+            self.q
+        )
+
+    @helpers.requires_array_function_protocol()
+    def test_dstack(self):
+        self.assertQuantityEqual(
+            np.dstack([self.q]*2),
+            self.Q_(np.dstack([self.q.m]*2), self.ureg.m)
+        )
+
+    @helpers.requires_array_function_protocol()
+    def test_hstack(self):
+        self.assertQuantityEqual(
+            np.hstack([self.q]*2),
+            self.Q_(np.hstack([self.q.m]*2), self.ureg.m)
+        )
+
+    @helpers.requires_array_function_protocol()
+    def test_vstack(self):
+        self.assertQuantityEqual(
+            np.vstack([self.q]*2),
+            self.Q_(np.vstack([self.q.m]*2), self.ureg.m)
+        )
+
+    @helpers.requires_array_function_protocol()
+    def test_block(self):
+        self.assertQuantityEqual(
+            np.block([self.q[0, :],self.q[1, :]]),
+            self.Q_([1,2,3,4], self.ureg.m)
+        )
+
+    @helpers.requires_array_function_protocol()
+    def test_append(self):
+        self.assertQuantityEqual(np.append(self.q, [[0, 0]] * self.ureg.m, axis=0),
+                                 [[1, 2], [3, 4], [0, 0]] * self.ureg.m)
+
+    def test_astype(self):
+        actual = self.q.astype(np.float32)
+        expected = self.Q_(np.array([[1., 2.], [3., 4.]], dtype=np.float32), 'm')
+        self.assertQuantityEqual(actual, expected)
+        self.assertEqual(actual.m.dtype, expected.m.dtype)
+
+    def test_item(self):
+        self.assertQuantityEqual(self.Q_([[0]], 'm').item(), 0 * self.ureg.m)
+
+
+class TestNumpyMathematicalFunctions(TestNumpyMethods):
+    # https://www.numpy.org/devdocs/reference/routines.math.html
+    # Trigonometric functions
+    @helpers.requires_array_function_protocol()
+    def test_unwrap(self):
+        self.assertQuantityEqual(np.unwrap([0,3*np.pi]*self.ureg.radians), [0,np.pi])
+        self.assertQuantityEqual(np.unwrap([0,540]*self.ureg.deg), [0,180]*self.ureg.deg)
+
+    # Rounding
+
+    @helpers.requires_array_function_protocol()
+    def test_fix(self):
+        self.assertQuantityEqual(np.fix(3.14 * self.ureg.m), 3.0 * self.ureg.m)
+        self.assertQuantityEqual(np.fix(3.0 * self.ureg.m), 3.0 * self.ureg.m)
+        self.assertQuantityEqual(
+            np.fix([2.1, 2.9, -2.1, -2.9] * self.ureg.m),
+            [2., 2., -2., -2.] * self.ureg.m
+        )
+    # Sums, products, differences
+
+    def test_prod(self):
+        self.assertEqual(self.q.prod(), 24 * self.ureg.m**4)
+
+    def test_sum(self):
+        self.assertEqual(self.q.sum(), 10*self.ureg.m)
+        self.assertQuantityEqual(self.q.sum(0), [4,     6]*self.ureg.m)
+        self.assertQuantityEqual(self.q.sum(1), [3, 7]*self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_sum_numpy_func(self):
+        self.assertQuantityEqual(np.sum(self.q, axis=0), [4, 6] * self.ureg.m)
+        self.assertRaises(OffsetUnitCalculusError, np.sum, self.q_temperature)
+
+    @helpers.requires_array_function_protocol()
+    def test_nansum_numpy_func(self):
+        self.assertQuantityEqual(np.nansum(self.q_nan, axis=0), [4, 2] * self.ureg.m)
+
+    def test_cumprod(self):
+        self.assertRaises(DimensionalityError, self.q.cumprod)
+        self.assertQuantityEqual((self.q / self.ureg.m).cumprod(), [1, 2, 6, 24])
+
+    @helpers.requires_array_function_protocol()
+    def test_cumprod_numpy_func(self):
+        self.assertRaises(DimensionalityError, np.cumprod, self.q)
+        self.assertRaises(DimensionalityError, np.cumproduct, self.q)
+        self.assertQuantityEqual(np.cumprod(self.q / self.ureg.m), [1, 2, 6, 24])
+        self.assertQuantityEqual(np.cumproduct(self.q / self.ureg.m), [1, 2, 6, 24])
+
+    @helpers.requires_array_function_protocol()
+    def test_nancumprod_numpy_func(self):
+        self.assertRaises(DimensionalityError, np.nancumprod, self.q_nan)
+        self.assertQuantityEqual(np.nancumprod(self.q_nan / self.ureg.m), [1, 2, 6, 6])
+
+    @helpers.requires_array_function_protocol()
+    def test_diff(self):
+        self.assertQuantityEqual(np.diff(self.q, 1), [[1], [1]] * self.ureg.m)
+        self.assertQuantityEqual(np.diff(self.q_temperature, 1), [[1], [1]] * self.ureg.delta_degC)
+
+    @helpers.requires_array_function_protocol()
+    def test_ediff1d(self):
+        self.assertQuantityEqual(np.ediff1d(self.q), [1, 1, 1] * self.ureg.m)
+        self.assertQuantityEqual(np.ediff1d(self.q_temperature), [1, 1, 1] * self.ureg.delta_degC)
+
+    @helpers.requires_array_function_protocol()
+    def test_gradient(self):
+        l = np.gradient([[1,1],[3,4]] * self.ureg.m, 1 * self.ureg.J)
+        self.assertQuantityEqual(l[0], [[2., 3.], [2., 3.]] * self.ureg.m / self.ureg.J)
+        self.assertQuantityEqual(l[1], [[0., 0.], [1., 1.]] * self.ureg.m / self.ureg.J)
+
+        l = np.gradient(self.Q_([[1,1],[3,4]] , self.ureg.degC), 1 * self.ureg.J)
+        self.assertQuantityEqual(l[0], [[2., 3.], [2., 3.]] * self.ureg.delta_degC / self.ureg.J)
+        self.assertQuantityEqual(l[1], [[0., 0.], [1., 1.]] * self.ureg.delta_degC / self.ureg.J)
+
+    @helpers.requires_array_function_protocol()
+    def test_cross(self):
+        a = [[3,-3, 1]] * self.ureg.kPa
+        b = [[4, 9, 2]] * self.ureg.m**2
+        self.assertQuantityEqual(np.cross(a, b), [[-15, -2, 39]] * self.ureg.kPa * self.ureg.m**2)
+
+    @helpers.requires_array_function_protocol()
+    def test_trapz(self):
+        self.assertQuantityEqual(np.trapz([1. ,2., 3., 4.] * self.ureg.J, dx=1*self.ureg.m), 7.5 * self.ureg.J*self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_dot(self):
+        self.assertQuantityEqual(self.q.ravel().dot(np.array([1, 0, 0, 1])), 5 * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_dot_numpy_func(self):
+        self.assertQuantityEqual(np.dot(self.q.ravel(), [0, 0, 1, 0] * self.ureg.dimensionless), 3 * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_einsum(self):
+        a = np.arange(25).reshape(5, 5) * self.ureg.m
+        b = np.arange(5) * self.ureg.m
+        self.assertQuantityEqual(np.einsum('ii', a), 60 * self.ureg.m)
+        self.assertQuantityEqual(np.einsum('ii->i', a), np.array([ 0,  6, 12, 18, 24]) * self.ureg.m)
+        self.assertQuantityEqual(np.einsum('i,i', b, b), 30 * self.ureg.m**2)
+        self.assertQuantityEqual(np.einsum('ij,j', a, b), np.array([ 30,  80, 130, 180, 230]) * self.ureg.m**2)
+
+    # Arithmetic operations
+    def test_addition_with_scalar(self):
+        a = np.array([0, 1, 2])
+        b = 10. * self.ureg('gram/kilogram')
+        self.assertQuantityAlmostEqual(a + b, self.Q_([0.01, 1.01, 2.01], self.ureg.dimensionless))
+        self.assertQuantityAlmostEqual(b + a, self.Q_([0.01, 1.01, 2.01], self.ureg.dimensionless))
+
+    def test_addition_with_incompatible_scalar(self):
+        a = np.array([0, 1, 2])
+        b = 1. * self.ureg.m
+        self.assertRaises(DimensionalityError, op.add, a, b)
+        self.assertRaises(DimensionalityError, op.add, b, a)
+
+    def test_power(self):
+        arr = np.array(range(3), dtype=np.float)
+        q = self.Q_(arr, 'meter')
+
+        for op_ in [op.pow, op.ipow, np.power]:
+            q_cp = copy.copy(q)
+            self.assertRaises(DimensionalityError, op_, 2., q_cp)
+            arr_cp = copy.copy(arr)
+            arr_cp = copy.copy(arr)
+            q_cp = copy.copy(q)
+            self.assertRaises(DimensionalityError, op_, q_cp, arr_cp)
+            q_cp = copy.copy(q)
+            q2_cp = copy.copy(q)
+            self.assertRaises(DimensionalityError, op_, q_cp, q2_cp)
+
+    @unittest.expectedFailure
+    @helpers.requires_numpy()
+    def test_exponentiation_array_exp_2(self):
+        arr = np.array(range(3), dtype=np.float)
+        #q = self.Q_(copy.copy(arr), None)
+        q = self.Q_(copy.copy(arr), 'meter')
+        arr_cp = copy.copy(arr)
+        q_cp = copy.copy(q)
+        # this fails as expected since numpy 1.8.0 but...
+        self.assertRaises(DimensionalityError, op.pow, arr_cp, q_cp)
+        # ..not for op.ipow !
+        # q_cp is treated as if it is an array. The units are ignored.
+        # Quantity.__ipow__ is never called
+        arr_cp = copy.copy(arr)
+        q_cp = copy.copy(q)
+        self.assertRaises(DimensionalityError, op.ipow, arr_cp, q_cp)
+
+
+class TestNumpyUnclassified(TestNumpyMethods):
+    def test_tolist(self):
+        self.assertEqual(self.q.tolist(), [[1*self.ureg.m, 2*self.ureg.m], [3*self.ureg.m, 4*self.ureg.m]])
+
+    def test_fill(self):
+        tmp = self.q
+        tmp.fill(6 * self.ureg.ft)
+        self.assertQuantityEqual(tmp, [[6, 6], [6, 6]] * self.ureg.ft)
+        tmp.fill(5 * self.ureg.m)
+        self.assertQuantityEqual(tmp, [[5, 5], [5, 5]] * self.ureg.m)
+
     def test_take(self):
         self.assertQuantityEqual(self.q.take([0,1,2,3]), self.q.flatten())
 
@@ -92,52 +410,139 @@ class TestNumpyMethods(QuantityTestCase):
         q.sort()
         self.assertQuantityEqual(q, [1, 2, 3, 4, 5, 6] * self.ureg.m)
 
+    @helpers.requires_array_function_protocol()
+    def test_sort_numpy_func(self):
+        q = [4, 5, 2, 3, 1, 6] * self.ureg.m
+        self.assertQuantityEqual(np.sort(q), [1, 2, 3, 4, 5, 6] * self.ureg.m)
+
     def test_argsort(self):
         q = [1, 4, 5, 6, 2, 9] * self.ureg.MeV
-        np.testing.assert_array_equal(q.argsort(), [0, 4, 1, 2, 3, 5])
+        self.assertNDArrayEqual(q.argsort(), [0, 4, 1, 2, 3, 5])
+
+    @helpers.requires_array_function_protocol()
+    def test_argsort_numpy_func(self):
+        self.assertNDArrayEqual(np.argsort(self.q, axis=0), np.array([[0, 0], [1, 1]]))
 
     def test_diagonal(self):
         q = [[1, 2, 3], [1, 2, 3], [1, 2, 3]] * self.ureg.m
         self.assertQuantityEqual(q.diagonal(offset=1), [2, 3] * self.ureg.m)
 
+    @helpers.requires_array_function_protocol()
+    def test_diagonal_numpy_func(self):
+        q = [[1, 2, 3], [1, 2, 3], [1, 2, 3]] * self.ureg.m
+        self.assertQuantityEqual(np.diagonal(q, offset=-1), [1, 2] * self.ureg.m)
+
     def test_compress(self):
         self.assertQuantityEqual(self.q.compress([False, True], axis=0),
                                  [[3, 4]] * self.ureg.m)
         self.assertQuantityEqual(self.q.compress([False, True], axis=1),
                                  [[2], [4]] * self.ureg.m)
 
+    @helpers.requires_array_function_protocol()
+    def test_compress(self):
+        self.assertQuantityEqual(np.compress([False, True], self.q, axis=1),
+                                 [[2], [4]] * self.ureg.m)
+
     def test_searchsorted(self):
         q = self.q.flatten()
-        np.testing.assert_array_equal(q.searchsorted([1.5, 2.5] * self.ureg.m),
-                                      [1, 2])
+        self.assertNDArrayEqual(q.searchsorted([1.5, 2.5] * self.ureg.m), [1, 2])
         q = self.q.flatten()
         self.assertRaises(DimensionalityError, q.searchsorted, [1.5, 2.5])
 
+    @helpers.requires_array_function_protocol()
     def test_searchsorted_numpy_func(self):
         """Test searchsorted as numpy function."""
         q = self.q.flatten()
-        np.testing.assert_array_equal(np.searchsorted(q, [1.5, 2.5] * self.ureg.m),
-                                      [1, 2])
+        self.assertNDArrayEqual(np.searchsorted(q, [1.5, 2.5] * self.ureg.m), [1, 2])
 
     def test_nonzero(self):
         q = [1, 0, 5, 6, 0, 9] * self.ureg.m
-        np.testing.assert_array_equal(q.nonzero()[0], [0, 2, 3, 5])
+        self.assertNDArrayEqual(q.nonzero()[0], [0, 2, 3, 5])
+
+    @helpers.requires_array_function_protocol()
+    def test_nonzero_numpy_func(self):
+        q = [1, 0, 5, 6, 0, 9] * self.ureg.m
+        self.assertNDArrayEqual(np.nonzero(q)[0], [0, 2, 3, 5])
+
+    @helpers.requires_array_function_protocol()
+    def test_count_nonzero_numpy_func(self):
+        q = [1, 0, 5, 6, 0, 9] * self.ureg.m
+        self.assertEqual(np.count_nonzero(q), 4)
 
     def test_max(self):
         self.assertEqual(self.q.max(), 4*self.ureg.m)
 
+    def test_max_numpy_func(self):
+        self.assertEqual(np.max(self.q), 4 * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_max_with_axis_arg(self):
+        self.assertQuantityEqual(np.max(self.q, axis=1), [2, 4] * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_max_with_initial_arg(self):
+        self.assertQuantityEqual(np.max(self.q[..., None], axis=2, initial=3 * self.ureg.m), [[3, 3], [3, 4]] * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_nanmax(self):
+        self.assertEqual(np.nanmax(self.q_nan), 3 * self.ureg.m)
+
     def test_argmax(self):
         self.assertEqual(self.q.argmax(), 3)
 
+    @helpers.requires_array_function_protocol()
+    def test_argmax_numpy_func(self):
+        self.assertNDArrayEqual(np.argmax(self.q, axis=0), np.array([1, 1]))
+
+    @helpers.requires_array_function_protocol()
+    def test_nanargmax_numpy_func(self):
+        self.assertNDArrayEqual(np.nanargmax(self.q_nan, axis=0), np.array([1, 0]))
+
+    def test_maximum(self):
+        self.assertQuantityEqual(np.maximum(self.q, self.Q_([0, 5], 'm')),
+                                 self.Q_([[1, 5], [3, 5]], 'm'))
+
     def test_min(self):
         self.assertEqual(self.q.min(), 1 * self.ureg.m)
 
+    @helpers.requires_array_function_protocol()
+    def test_min_numpy_func(self):
+        self.assertEqual(np.min(self.q), 1 * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_min_with_axis_arg(self):
+        self.assertQuantityEqual(np.min(self.q, axis=1), [1, 3] * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_min_with_initial_arg(self):
+        self.assertQuantityEqual(np.min(self.q[..., None], axis=2, initial=3 * self.ureg.m), [[1, 2], [3, 3]] * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_nanmin(self):
+        self.assertEqual(np.nanmin(self.q_nan), 1 * self.ureg.m)
+
     def test_argmin(self):
         self.assertEqual(self.q.argmin(), 0)
 
+    @helpers.requires_array_function_protocol()
+    def test_argmin_numpy_func(self):
+        self.assertNDArrayEqual(np.argmin(self.q, axis=0), np.array([0, 0]))
+
+    @helpers.requires_array_function_protocol()
+    def test_nanargmin_numpy_func(self):
+        self.assertNDArrayEqual(np.nanargmin(self.q_nan, axis=0), np.array([0, 0]))
+
+    def test_minimum(self):
+        self.assertQuantityEqual(np.minimum(self.q, self.Q_([0, 5], 'm')),
+                                 self.Q_([[0, 2], [0, 4]], 'm'))
+
     def test_ptp(self):
         self.assertEqual(self.q.ptp(), 3 * self.ureg.m)
 
+    @helpers.requires_array_function_protocol()
+    def test_ptp_numpy_func(self):
+        self.assertQuantityEqual(np.ptp(self.q, axis=0), [2, 2] * self.ureg.m)
+
     def test_clip(self):
         self.assertQuantityEqual(
             self.q.clip(max=2*self.ureg.m),
@@ -154,27 +559,78 @@ class TestNumpyMethods(QuantityTestCase):
         self.assertRaises(DimensionalityError, self.q.clip, self.ureg.J)
         self.assertRaises(DimensionalityError, self.q.clip, 1)
 
+    @helpers.requires_array_function_protocol()
+    def test_clip_numpy_func(self):
+        self.assertQuantityEqual(np.clip(self.q, 150 * self.ureg.cm, None), [[1.5, 2], [3, 4]] * self.ureg.m)
+
     def test_round(self):
         q = [1, 1.33, 5.67, 22] * self.ureg.m
         self.assertQuantityEqual(q.round(0), [1, 1, 6, 22] * self.ureg.m)
         self.assertQuantityEqual(q.round(-1), [0, 0, 10, 20] * self.ureg.m)
         self.assertQuantityEqual(q.round(1), [1, 1.3, 5.7, 22] * self.ureg.m)
 
+    @helpers.requires_array_function_protocol()
+    def test_round_numpy_func(self):
+        self.assertQuantityEqual(np.around(1.0275 * self.ureg.m, decimals=2), 1.03 * self.ureg.m)
+        self.assertQuantityEqual(np.round_(1.0275 * self.ureg.m, decimals=2), 1.03 * self.ureg.m)
+
     def test_trace(self):
         self.assertEqual(self.q.trace(), (1+4) * self.ureg.m)
 
     def test_cumsum(self):
         self.assertQuantityEqual(self.q.cumsum(), [1, 3, 6, 10] * self.ureg.m)
 
+    @helpers.requires_array_function_protocol()
+    def test_cumsum_numpy_func(self):
+        self.assertQuantityEqual(np.cumsum(self.q, axis=0), [[1, 2], [4, 6]] * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_nancumsum_numpy_func(self):
+        self.assertQuantityEqual(np.nancumsum(self.q_nan, axis=0), [[1, 2], [4, 2]] * self.ureg.m)
+
     def test_mean(self):
         self.assertEqual(self.q.mean(), 2.5 * self.ureg.m)
 
+    @helpers.requires_array_function_protocol()
+    def test_mean_numpy_func(self):
+        self.assertEqual(np.mean(self.q), 2.5 * self.ureg.m)
+        self.assertEqual(np.mean(self.q_temperature), self.Q_(2.5, self.ureg.degC))
+
+    @helpers.requires_array_function_protocol()
+    def test_nanmean_numpy_func(self):
+        self.assertEqual(np.nanmean(self.q_nan), 2 * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_average_numpy_func(self):
+        self.assertQuantityAlmostEqual(np.average(self.q, axis=0, weights=[1, 2]), [2.33333, 3.33333] * self.ureg.m, rtol=1e-5)
+
+    @helpers.requires_array_function_protocol()
+    def test_median_numpy_func(self):
+        self.assertEqual(np.median(self.q), 2.5 * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_nanmedian_numpy_func(self):
+        self.assertEqual(np.nanmedian(self.q_nan), 2 * self.ureg.m)
+
     def test_var(self):
         self.assertEqual(self.q.var(), 1.25*self.ureg.m**2)
 
+    @helpers.requires_array_function_protocol()
+    def test_var_numpy_func(self):
+        self.assertEqual(np.var(self.q), 1.25*self.ureg.m**2)
+
+    @helpers.requires_array_function_protocol()
+    def test_nanvar_numpy_func(self):
+        self.assertQuantityAlmostEqual(np.nanvar(self.q_nan), 0.66667*self.ureg.m**2, rtol=1e-5)
+
     def test_std(self):
         self.assertQuantityAlmostEqual(self.q.std(), 1.11803*self.ureg.m, rtol=1e-5)
 
+    @helpers.requires_array_function_protocol()
+    def test_std_numpy_func(self):
+        self.assertQuantityAlmostEqual(np.std(self.q), 1.11803*self.ureg.m, rtol=1e-5)
+        self.assertRaises(OffsetUnitCalculusError, np.std, self.q_temperature)
+
     def test_prod(self):
         self.assertEqual(self.q.prod(), 24 * self.ureg.m**4)
 
@@ -182,6 +638,10 @@ class TestNumpyMethods(QuantityTestCase):
         self.assertRaises(DimensionalityError, self.q.cumprod)
         self.assertQuantityEqual((self.q / self.ureg.m).cumprod(), [1, 2, 6, 24])
 
+    @helpers.requires_array_function_protocol()
+    def test_nanstd_numpy_func(self):
+        self.assertQuantityAlmostEqual(np.nanstd(self.q_nan), 0.81650 * self.ureg.m, rtol=1e-5)
+
     @helpers.requires_numpy_previous_than('1.10')
     def test_integer_div(self):
         a = [1] * self.ureg.m
@@ -259,7 +719,7 @@ class TestNumpyMethods(QuantityTestCase):
 
         def pickle_test(q):
             pq = pickle.loads(pickle.dumps(q))
-            np.testing.assert_array_equal(q.magnitude, pq.magnitude)
+            self.assertNDArrayEqual(q.magnitude, pq.magnitude)
             self.assertEqual(q.units, pq.units)
 
         pickle_test([10, 20]*self.ureg.m)
@@ -277,88 +737,129 @@ class TestNumpyMethods(QuantityTestCase):
         u.shape = 4, 3
         self.assertEqual(u.magnitude.shape, (4, 3))
 
-
-@helpers.requires_numpy()
-class TestNumpyNeedsSubclassing(TestUFuncs):
-
-    FORCE_NDARRAY = True
-
-    @property
-    def q(self):
-        return [1. ,2., 3., 4.] * self.ureg.J
-
-    @unittest.expectedFailure
-    def test_unwrap(self):
-        """unwrap depends on diff
-        """
-        self.assertQuantityEqual(np.unwrap([0,3*np.pi]*self.ureg.radians), [0,np.pi])
-        self.assertQuantityEqual(np.unwrap([0,540]*self.ureg.deg), [0,180]*self.ureg.deg)
-
-    @unittest.expectedFailure
-    def test_trapz(self):
-        """Units are erased by asanyarray, Quantity does not inherit from NDArray
-        """
-        self.assertQuantityEqual(np.trapz(self.q, dx=1*self.ureg.m), 7.5 * self.ureg.J*self.ureg.m)
-
-    @unittest.expectedFailure
-    def test_diff(self):
-        """Units are erased by asanyarray, Quantity does not inherit from NDArray
-        """
-        self.assertQuantityEqual(np.diff(self.q, 1), [1, 1, 1] * self.ureg.J)
-
-    @unittest.expectedFailure
-    def test_ediff1d(self):
-        """Units are erased by asanyarray, Quantity does not inherit from NDArray
-        """
-        self.assertQuantityEqual(np.ediff1d(self.q, 1 * self.ureg.J), [1, 1, 1] * self.ureg.J)
-
-    @unittest.expectedFailure
-    def test_fix(self):
-        """Units are erased by asanyarray, Quantity does not inherit from NDArray
-        """
-        self.assertQuantityEqual(np.fix(3.14 * self.ureg.m), 3.0 * self.ureg.m)
-        self.assertQuantityEqual(np.fix(3.0 * self.ureg.m), 3.0 * self.ureg.m)
-        self.assertQuantityEqual(
-            np.fix([2.1, 2.9, -2.1, -2.9] * self.ureg.m),
-            [2., 2., -2., -2.] * self.ureg.m
-        )
-
-    @unittest.expectedFailure
-    def test_gradient(self):
-        """shape is a property not a function
-        """
-        l = np.gradient([[1,1],[3,4]] * self.ureg.J, 1 * self.ureg.m)
-        self.assertQuantityEqual(l[0], [[2., 3.], [2., 3.]] * self.ureg.J / self.ureg.m)
-        self.assertQuantityEqual(l[1], [[0., 0.], [1., 1.]] * self.ureg.J / self.ureg.m)
-
-    @unittest.expectedFailure
-    def test_cross(self):
-        """Units are erased by asarray, Quantity does not inherit from NDArray
-        """
-        a = [[3,-3, 1]] * self.ureg.kPa
-        b = [[4, 9, 2]] * self.ureg.m**2
-        self.assertQuantityEqual(np.cross(a, b), [-15, -2, 39] * self.ureg.kPa * self.ureg.m**2)
-
-    @unittest.expectedFailure
-    def test_power(self):
-        """This is not supported as different elements might end up with different units
-
-        eg. ([1, 1] * m) ** [2, 3]
-
-        Must force exponent to single value
-        """
-        self._test2(np.power, self.q1,
-                    (self.qless, np.asarray([1., 2, 3, 4])),
-                    (self.q2, ),)
-
-    @unittest.expectedFailure
-    def test_ones_like(self):
-        """Units are erased by emptyarra, Quantity does not inherit from NDArray
-        """
-        self._test1(np.ones_like,
-                    (self.q2, self.qs, self.qless, self.qi),
-                    (),
-                    2)
+    @helpers.requires_array_function_protocol()
+    def test_shape_numpy_func(self):
+        self.assertEqual(np.shape(self.q), (2, 2))
+
+    @helpers.requires_array_function_protocol()
+    def test_alen_numpy_func(self):
+        self.assertEqual(np.alen(self.q), 2)
+
+    @helpers.requires_array_function_protocol()
+    def test_ndim_numpy_func(self):
+        self.assertEqual(np.ndim(self.q), 2)
+
+    @helpers.requires_array_function_protocol()
+    def test_copy_numpy_func(self):
+        q_copy = np.copy(self.q)
+        self.assertQuantityEqual(self.q, q_copy)
+        self.assertIsNot(self.q, q_copy)
+
+    @helpers.requires_array_function_protocol()
+    def test_trim_zeros_numpy_func(self):
+        q = [0, 4, 3, 0, 2, 2, 0, 0, 0] * self.ureg.m
+        self.assertQuantityEqual(np.trim_zeros(q), [4, 3, 0, 2, 2] * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_result_type_numpy_func(self):
+        self.assertEqual(np.result_type(self.q), np.dtype('int64'))
+
+    @helpers.requires_array_function_protocol()
+    def test_nan_to_num_numpy_func(self):
+        self.assertQuantityEqual(np.nan_to_num(self.q_nan, nan=-999 * self.ureg.mm),
+                                 [[1, 2], [3, -0.999]] * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_meshgrid_numpy_func(self):
+        x = [1, 2] * self.ureg.m
+        y = [0, 50, 100] * self.ureg.mm
+        xx, yy = np.meshgrid(x, y)
+        self.assertQuantityEqual(xx, [[1, 2], [1, 2], [1, 2]] * self.ureg.m)
+        self.assertQuantityEqual(yy, [[0, 0], [50, 50], [100, 100]] * self.ureg.mm)
+
+    @helpers.requires_array_function_protocol()
+    def test_isclose_numpy_func(self):
+        q2 = [[1000.05, 2000], [3000.00007, 4001]] * self.ureg.mm
+        self.assertNDArrayEqual(np.isclose(self.q, q2), np.array([[False, True], [True, False]]))
+
+    @helpers.requires_array_function_protocol()
+    def test_interp_numpy_func(self):
+        x = [1, 4] * self.ureg.m
+        xp = np.linspace(0, 3, 5) * self.ureg.m
+        fp = self.Q_([0, 5, 10, 15, 20], self.ureg.degC)
+        self.assertQuantityAlmostEqual(np.interp(x, xp, fp), self.Q_([6.66667, 20.], self.ureg.degC), rtol=1e-5)
+
+    def test_comparisons(self):
+        self.assertNDArrayEqual(self.q > 2 * self.ureg.m, np.array([[False, False], [True, True]]))
+        self.assertNDArrayEqual(self.q < 2 * self.ureg.m, np.array([[True, False], [False, False]]))
+
+    @helpers.requires_array_function_protocol()
+    def test_where(self):
+        self.assertQuantityEqual(np.where(self.q >= 2 * self.ureg.m, self.q, 20 * self.ureg.m),
+                                 [[20, 2], [3, 4]] * self.ureg.m)
+        self.assertQuantityEqual(np.where(self.q >= 2 * self.ureg.m, self.q, 0),
+                                 [[0, 2], [3, 4]] * self.ureg.m)
+        self.assertQuantityEqual(np.where(self.q >= 2 * self.ureg.m, self.q, np.nan),
+                                 [[np.nan, 2], [3, 4]] * self.ureg.m)
+        self.assertQuantityEqual(np.where(self.q >= 3 * self.ureg.m, 0, self.q),
+                                 [[1, 2], [0, 0]] * self.ureg.m)
+        self.assertQuantityEqual(np.where(self.q >= 3 * self.ureg.m, np.nan, self.q),
+                                 [[1, 2], [np.nan, np.nan]] * self.ureg.m)
+        self.assertQuantityEqual(np.where(self.q >= 2 * self.ureg.m, self.q, np.array(np.nan)),
+                                 [[np.nan, 2], [3, 4]] * self.ureg.m)
+        self.assertQuantityEqual(np.where(self.q >= 3 * self.ureg.m, np.array(np.nan), self.q),
+                                 [[1, 2], [np.nan, np.nan]] * self.ureg.m)
+        self.assertRaises(DimensionalityError, np.where, self.q < 2 * self.ureg.m, self.q, 0 * self.ureg.J)
+
+    @helpers.requires_array_function_protocol()
+    def test_fabs(self):
+        self.assertQuantityEqual(np.fabs(self.q - 2 * self.ureg.m), self.Q_([[1, 0], [1, 2]], 'm'))
+
+    @helpers.requires_array_function_protocol()
+    def test_isin(self):
+        self.assertNDArrayEqual(np.isin(self.q, self.Q_([0, 2, 4], 'm')),
+                                np.array([[False, True], [False, True]]))
+        self.assertNDArrayEqual(np.isin(self.q, self.Q_([0, 2, 4], 'J')),
+                                np.array([[False, False], [False, False]]))
+        self.assertNDArrayEqual(np.isin(self.q, [self.Q_(2, 'm'), self.Q_(4, 'J')]),
+                                np.array([[False, True], [False, False]]))
+        self.assertNDArrayEqual(np.isin(self.q, self.q.m),
+                                np.array([[False, False], [False, False]]))
+        self.assertNDArrayEqual(np.isin(self.q / self.ureg.cm, [1, 3]),
+                                np.array([[True, False], [True, False]]))
+        self.assertRaises(ValueError, np.isin, self.q.m, self.q)
+
+    @helpers.requires_array_function_protocol()
+    def test_percentile(self):
+        self.assertQuantityEqual(np.percentile(self.q, 25), self.Q_(1.75, 'm'))
+
+    @helpers.requires_array_function_protocol()
+    def test_nanpercentile(self):
+        self.assertQuantityEqual(np.nanpercentile(self.q_nan, 25), self.Q_(1.5, 'm'))
+
+    @helpers.requires_array_function_protocol()
+    def test_copyto(self):
+        a = self.q.m
+        q = copy.copy(self.q)
+        np.copyto(q, 2 * q, where=[True, False])
+        self.assertQuantityEqual(q, self.Q_([[2, 2], [6, 4]], 'm'))
+        np.copyto(q, 0, where=[[False, False], [True, False]])
+        self.assertQuantityEqual(q, self.Q_([[2, 2], [0, 4]], 'm'))
+        np.copyto(a, q)
+        self.assertNDArrayEqual(a, np.array([[2, 2], [0, 4]]))
+
+    @helpers.requires_array_function_protocol()
+    def test_tile(self):
+        self.assertQuantityEqual(np.tile(self.q, 2), np.array([[1, 2, 1, 2], [3, 4, 3, 4]]) * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_rot90(self):
+        self.assertQuantityEqual(np.rot90(self.q), np.array([[2, 4], [1, 3]]) * self.ureg.m)
+
+    @helpers.requires_array_function_protocol()
+    def test_insert(self):
+        self.assertQuantityEqual(np.insert(self.q, 1, 0 * self.ureg.m, axis=1),
+                                 np.array([[1, 0, 2], [3, 0, 4]]) * self.ureg.m)
 
 
 @unittest.skip
@@ -435,39 +936,3 @@ class TestBitTwiddlingUfuncs(TestUFuncs):
                     (self.qless, 2),
                     (self.q1, self.q2, self.qs, ),
                     'same')
-
-
-class TestNDArrayQuantityMath(QuantityTestCase):
-
-    @helpers.requires_numpy()
-    def test_exponentiation_array_exp(self):
-        arr = np.array(range(3), dtype=np.float)
-        q = self.Q_(arr, 'meter')
-
-        for op_ in [op.pow, op.ipow]:
-            q_cp = copy.copy(q)
-            self.assertRaises(DimensionalityError, op_, 2., q_cp)
-            arr_cp = copy.copy(arr)
-            arr_cp = copy.copy(arr)
-            q_cp = copy.copy(q)
-            self.assertRaises(DimensionalityError, op_, q_cp, arr_cp)
-            q_cp = copy.copy(q)
-            q2_cp = copy.copy(q)
-            self.assertRaises(DimensionalityError, op_, q_cp, q2_cp)
-
-    @unittest.expectedFailure
-    @helpers.requires_numpy()
-    def test_exponentiation_array_exp_2(self):
-        arr = np.array(range(3), dtype=np.float)
-        #q = self.Q_(copy.copy(arr), None)
-        q = self.Q_(copy.copy(arr), 'meter')
-        arr_cp = copy.copy(arr)
-        q_cp = copy.copy(q)
-        # this fails as expected since numpy 1.8.0 but...
-        self.assertRaises(DimensionalityError, op.pow, arr_cp, q_cp)
-        # ..not for op.ipow !
-        # q_cp is treated as if it is an array. The units are ignored.
-        # Quantity.__ipow__ is never called
-        arr_cp = copy.copy(arr)
-        q_cp = copy.copy(q)
-        self.assertRaises(DimensionalityError, op.ipow, arr_cp, q_cp)
diff --git a/pint/testsuite/test_numpy_func.py b/pint/testsuite/test_numpy_func.py
new file mode 100644
index 0000000..31772d5
--- /dev/null
+++ b/pint/testsuite/test_numpy_func.py
@@ -0,0 +1,156 @@
+# -*- coding: utf-8 -*-
+
+import pint.numpy_func
+
+from pint import DimensionalityError, OffsetUnitCalculusError
+from pint.compat import np
+from pint.testsuite.test_numpy import TestNumpyMethods
+from pint.numpy_func import (_is_quantity, _is_quantity_sequence, _get_first_input_units, 
+                             convert_to_consistent_units, unwrap_and_wrap_consistent_units,
+                             get_op_output_unit, implements, numpy_wrap)
+from unittest.mock import patch
+
+
+class TestNumPyFuncUtils(TestNumpyMethods):
+
+    @patch('pint.numpy_func.HANDLED_FUNCTIONS', {})
+    @patch('pint.numpy_func.HANDLED_UFUNCS', {})
+    def test_implements(self):
+        # Test for functions
+        @implements('test', 'function')
+        def test_function():
+            pass
+
+        self.assertEqual(pint.numpy_func.HANDLED_FUNCTIONS['test'], test_function)
+
+        # Test for ufuncs
+        @implements('test', 'ufunc')
+        def test_ufunc():
+            pass
+
+        self.assertEqual(pint.numpy_func.HANDLED_UFUNCS['test'], test_ufunc)
+
+        # Test for invalid func type
+        with self.assertRaises(ValueError):
+            @implements('test', 'invalid')
+            def test_invalid():
+                pass
+
+    def test_is_quantity(self):
+        self.assertTrue(_is_quantity(self.Q_(0)))
+        self.assertTrue(_is_quantity(np.arange(4) * self.ureg.m))
+        self.assertFalse(_is_quantity(1.))
+        self.assertFalse(_is_quantity(np.array([1, 1, 2, 3, 5, 8])))
+        self.assertFalse(_is_quantity('not-a-quantity'))
+        # TODO (#905 follow-up): test other duck arrays that wrap or are wrapped by Pint
+
+    def test_is_quantity_sequence(self):
+        self.assertTrue(_is_quantity_sequence((self.Q_(0, 'm'), self.Q_(32., 'degF'))))
+        self.assertTrue(_is_quantity_sequence(np.arange(4) * self.ureg.m))
+        self.assertFalse(_is_quantity_sequence((self.Q_(0), True)))
+        self.assertFalse(_is_quantity_sequence([1, 3, 5]))
+        self.assertFalse(_is_quantity_sequence(9 * self.ureg.m))
+        self.assertFalse(_is_quantity_sequence(np.arange(4)))
+        self.assertFalse(_is_quantity_sequence('0123'))
+
+    def test_convert_to_consistent_units_with_pre_calc_units(self):
+        args, kwargs = convert_to_consistent_units(
+            self.Q_(50, 'cm'),
+            np.arange(4).reshape(2, 2) * self.ureg.m,
+            [0.042] * self.ureg.km,
+            (self.Q_(0, 'm'), self.Q_(1, 'dam')),
+            a=6378 * self.ureg.km,
+            pre_calc_units=self.ureg.meter)
+        self.assertEqual(args[0], 0.5)
+        self.assertNDArrayEqual(args[1], np.array([[0, 1], [2, 3]]))
+        self.assertNDArrayEqual(args[2], np.array([42]))
+        self.assertEqual(args[3][0], 0)
+        self.assertEqual(args[3][1], 10)
+        self.assertEqual(kwargs['a'], 6.378e6)
+
+    def test_convert_to_consistent_units_with_dimensionless(self):
+        args, kwargs = convert_to_consistent_units(
+            np.arange(2),
+            pre_calc_units=self.ureg.g/self.ureg.kg)
+        self.assertNDArrayEqual(args[0], np.array([0, 1000]))
+        self.assertEqual(kwargs, {})
+
+    def test_convert_to_consistent_units_with_dimensionality_error(self):
+        self.assertRaises(DimensionalityError, convert_to_consistent_units,
+                          self.Q_(32., 'degF'), pre_calc_units=self.ureg.meter)
+        self.assertRaises(DimensionalityError, convert_to_consistent_units, np.arange(4),
+                          pre_calc_units=self.ureg.meter)
+
+    def test_convert_to_consistent_units_without_pre_calc_units(self):
+        args, kwargs = convert_to_consistent_units(
+            (self.Q_(0), self.Q_(10, 'degC')),
+            [1, 2, 3, 5, 7] * self.ureg.m,
+            pre_calc_units=None)
+        self.assertEqual(args[0][0], 0)
+        self.assertEqual(args[0][1], 10)
+        self.assertNDArrayEqual(args[1], np.array([1, 2, 3, 5, 7]))
+        self.assertEqual(kwargs, {})
+
+    def test_unwrap_and_wrap_constistent_units(self):
+        (a, ), output_wrap_a = unwrap_and_wrap_consistent_units([2, 4, 8] * self.ureg.m)
+        (b, c), output_wrap_c = unwrap_and_wrap_consistent_units(
+            np.arange(4), self.Q_(1, 'g/kg'))
+
+        self.assertNDArrayEqual(a, np.array([2, 4, 8]))
+        self.assertNDArrayEqual(b, np.array([0, 1000, 2000, 3000]))
+        self.assertEqual(c, 1)
+
+        self.assertQuantityEqual(output_wrap_a(0), 0 * self.ureg.m)
+        self.assertQuantityEqual(output_wrap_c(0), self.Q_(0, 'g/kg'))
+
+    def test_op_output_unit_sum(self):
+        self.assertEqual(get_op_output_unit('sum', self.ureg.m), self.ureg.m)
+        self.assertRaises(OffsetUnitCalculusError, get_op_output_unit, 'sum', self.ureg.degC)
+
+    def test_op_output_unit_mul(self):
+        self.assertEqual(get_op_output_unit('mul', self.ureg.s,
+                                            (self.Q_(1, 'm'), self.Q_(1, 'm**2'))),
+                         self.ureg.m**3)
+
+    def test_op_output_unit_delta(self):
+        self.assertEqual(get_op_output_unit('delta', self.ureg.m), self.ureg.m)
+        self.assertEqual(get_op_output_unit('delta', self.ureg.degC), self.ureg.delta_degC)
+
+    def test_op_output_unit_delta_div(self):
+        self.assertEqual(get_op_output_unit('delta,div', self.ureg.m,
+                                            (self.Q_(1, 'm'), self.Q_(1, 's'))),
+                         self.ureg.m / self.ureg.s)
+        self.assertEqual(get_op_output_unit('delta,div', self.ureg.degC,
+                                            (self.Q_(1, 'degC'), self.Q_(1, 'm'))),
+                         self.ureg.delta_degC / self.ureg.m)
+
+    def test_op_output_unit_div(self):
+        self.assertEqual(get_op_output_unit('div', self.ureg.m,
+                                            (self.Q_(1, 'm'), self.Q_(1, 's'),
+                                             self.Q_(1, 'K'))),
+                         self.ureg.m / self.ureg.s / self.ureg.K)
+        self.assertEqual(get_op_output_unit('div', self.ureg.s, (1, self.Q_(1, 's'))),
+                         self.ureg.s**-1)
+
+    def test_op_output_unit_variance(self):
+        self.assertEqual(get_op_output_unit('variance', self.ureg.m), self.ureg.m**2)
+        self.assertRaises(OffsetUnitCalculusError, get_op_output_unit, 'variance',
+                          self.ureg.degC)
+
+    def test_op_output_unit_square(self):
+        self.assertEqual(get_op_output_unit('square', self.ureg.m), self.ureg.m**2)
+
+    def test_op_output_unit_sqrt(self):
+        self.assertEqual(get_op_output_unit('sqrt', self.ureg.m), self.ureg.m**0.5)
+
+    def test_op_output_unit_reciprocal(self):
+        self.assertEqual(get_op_output_unit('reciprocal', self.ureg.m), self.ureg.m**-1)
+
+    def test_op_output_unit_size(self):
+        self.assertEqual(get_op_output_unit('size', self.ureg.m, size=3), self.ureg.m**3)
+        self.assertRaises(ValueError, get_op_output_unit, 'size', self.ureg.m)
+
+    def test_numpy_wrap(self):
+        self.assertRaises(ValueError, numpy_wrap, 'invalid', np.ones, [], {}, [])
+        # TODO (#905 follow-up): test that NotImplemented is returned when upcast types
+        # present
diff --git a/pint/testsuite/test_quantity.py b/pint/testsuite/test_quantity.py
index e15923d..dc82ff2 100644
--- a/pint/testsuite/test_quantity.py
+++ b/pint/testsuite/test_quantity.py
@@ -4,12 +4,14 @@ import copy
 import datetime
 import math
 import operator as op
+import warnings
 
 from pint import DimensionalityError, OffsetUnitCalculusError, UnitRegistry
 from pint.unit import UnitsContainer
-from pint.compat import np
+from pint.compat import BehaviorChangeWarning, np
 from pint.testsuite import QuantityTestCase, helpers
 from pint.testsuite.parameterized import ParameterizedTestCase
+from unittest.mock import patch
 
 
 class TestQuantity(QuantityTestCase):
@@ -402,6 +404,16 @@ class TestQuantity(QuantityTestCase):
         with self.assertRaises(TypeError):
             iter(x)
 
+    @helpers.requires_array_function_protocol()
+    @patch('pint.quantity.SKIP_ARRAY_FUNCTION_CHANGE_WARNING', False)
+    def test_array_function_warning_on_creation(self):
+        # Test that warning is raised on first creation, but not second
+        with self.assertWarns(BehaviorChangeWarning):
+            self.Q_([])
+        with warnings.catch_warnings():
+            warnings.filterwarnings('error')
+            self.Q_([])
+
 
 class TestQuantityToCompact(QuantityTestCase):
 
@@ -1303,6 +1315,16 @@ class TestOffsetUnitMath(QuantityTestCase, ParameterizedTestCase):
                 in1_cp = copy.copy(in1)
                 self.assertQuantityAlmostEqual(op.ipow(in1_cp, in2), expected)
 
+    # matmul is only a ufunc since 1.16
+    @helpers.requires_numpy_at_least('1.16')
+    def test_matmul_with_numpy(self):
+        A = [[1, 2], [3, 4]] * self.ureg.m
+        B = np.array([[0, -1], [-1, 0]])
+        b = [[1], [0]] * self.ureg.m
+        self.assertQuantityEqual(A @ B, [[-2, -1], [-4, -3]] * self.ureg.m)
+        self.assertQuantityEqual(A @ b, [[1], [3]] * self.ureg.m**2)
+        self.assertQuantityEqual(B @ b, [[0], [-1]] * self.ureg.m)
+
 
 class TestDimensionReduction(QuantityTestCase):
     def _calc_mass(self, ureg):