diff --git a/CHANGES b/CHANGES index 0395b7fc0..1f0a95379 100644 --- a/CHANGES +++ b/CHANGES @@ -4,6 +4,13 @@ Pint Changelog 0.10 (unreleased) ----------------- +- Documentation on Pint's array type compatibility has been added to the NumPy support + page, including a graph of the duck array type casting hierarchy as understood by Pint + for N-dimensional arrays. + (Issue #963, Thanks Jon Thielen, Stephan Hoyer, and Guido Imperiale) +- Improved compatibility for downcast duck array types like Sparse and Masked Arrays. A + collection of basic tests has been added. + (Issue #963, Thanks Jon Thielen) - Improvements to wraps and check: - fail upon decoration (not execution) by checking wrapped function signature against wraps/check arguments. @@ -12,7 +19,7 @@ Pint Changelog (might BREAK code not conforming to documentation) - when strict=True, strings that can be parsed to quantities are accepted as arguments. - Add revolutions per second (rps) -- Improved compatbility for upcast types like xarray's DataArray or Dataset, to which +- Improved compatibility for upcast types like xarray's DataArray or Dataset, to which Pint Quantities now fully defer for arithmetic and NumPy operations. A collection of basic tests for proper deferral has been added (for full integration tests, see xarray's test suite). The list of upcast types is available at diff --git a/docs/numpy.ipynb b/docs/numpy.ipynb new file mode 100644 index 000000000..02f8f3fce --- /dev/null +++ b/docs/numpy.ipynb @@ -0,0 +1,763 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "NumPy Support\n", + "=============\n", + "\n", + "The magnitude of a Pint quantity can be of any numerical scalar type, and you are free\n", + "to choose it according to your needs. For numerical applications requiring arrays, it is\n", + "quite convenient to use [NumPy ndarray](http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html) (or [ndarray-like types supporting NEP-18](https://numpy.org/neps/nep-0018-array-function-protocol.html)),\n", + "and therefore these are the array types supported by Pint.\n", + "\n", + "First, we import the relevant packages:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "# Import NumPy\n", + "import numpy as np\n", + "\n", + "# Disable Pint's old fallback behavior (must come before importing Pint)\n", + "import os\n", + "os.environ['PINT_ARRAY_PROTOCOL_FALLBACK'] = \"0\"\n", + "\n", + "# Import Pint\n", + "import pint\n", + "ureg = pint.UnitRegistry()\n", + "Q_ = ureg.Quantity\n", + "\n", + "# Silence NEP 18 warning\n", + "import warnings\n", + "with warnings.catch_warnings():\n", + " warnings.simplefilter(\"ignore\")\n", + " Q_([])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "and then we create a quantity the standard way" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[3.0 4.0] meter\n" + ] + } + ], + "source": [ + "legs1 = Q_(np.asarray([3., 4.]), 'meter')\n", + "print(legs1)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[3.0 4.0] meter\n" + ] + } + ], + "source": [ + "legs1 = [3., 4.] * ureg.meter\n", + "print(legs1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "All usual Pint methods can be used with this quantity. For example:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.003 0.004] kilometer\n" + ] + } + ], + "source": [ + "print(legs1.to('kilometer'))" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[length]\n" + ] + } + ], + "source": [ + "print(legs1.dimensionality)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Cannot convert from 'meter' ([length]) to 'joule' ([length] ** 2 * [mass] / [time] ** 2)\n" + ] + } + ], + "source": [ + "try:\n", + " legs1.to('joule')\n", + "except pint.DimensionalityError as exc:\n", + " print(exc)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "NumPy functions are supported by Pint. For example if we define:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[400.0 300.0] centimeter\n" + ] + } + ], + "source": [ + "legs2 = [400., 300.] * ureg.centimeter\n", + "print(legs2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "we can calculate the hypotenuse of the right triangles with legs1 and legs2." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[5.0 5.0] meter\n" + ] + } + ], + "source": [ + "hyps = np.hypot(legs1, legs2)\n", + "print(hyps)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice that before the `np.hypot` was used, the numerical value of legs2 was\n", + "internally converted to the units of legs1 as expected.\n", + "\n", + "Similarly, when you apply a function that expects angles in radians, a conversion\n", + "is applied before the requested calculation:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.6435011087932843 0.9272952180016123] radian\n" + ] + } + ], + "source": [ + "angles = np.arccos(legs2/hyps)\n", + "print(angles)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can convert the result to degrees using usual unit conversion:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[36.86989764584401 53.13010235415599] degree\n" + ] + } + ], + "source": [ + "print(angles.to('degree'))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Applying a function that expects angles to a quantity with a different dimensionality\n", + "results in an error:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Cannot convert from 'centimeter' ([length]) to 'dimensionless' (dimensionless)\n" + ] + } + ], + "source": [ + "try:\n", + " np.arccos(legs2)\n", + "except pint.DimensionalityError as exc:\n", + " print(exc)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Function/Method Support\n", + "-----------------------\n", + "\n", + "The following [ufuncs](http://docs.scipy.org/doc/numpy/reference/ufuncs.html) can be applied to a Quantity object:\n", + "\n", + "- **Math operations**: `add`, `subtract`, `multiply`, `divide`, `logaddexp`, `logaddexp2`, `true_divide`, `floor_divide`, `negative`, `remainder`, `mod`, `fmod`, `absolute`, `rint`, `sign`, `conj`, `exp`, `exp2`, `log`, `log2`, `log10`, `expm1`, `log1p`, `sqrt`, `square`, `reciprocal`\n", + "- **Trigonometric functions**: `sin`, `cos`, `tan`, `arcsin`, `arccos`, `arctan`, `arctan2`, `hypot`, `sinh`, `cosh`, `tanh`, `arcsinh`, `arccosh`, `arctanh`\n", + "- **Comparison functions**: `greater`, `greater_equal`, `less`, `less_equal`, `not_equal`, `equal`\n", + "- **Floating functions**: `isreal`, `iscomplex`, `isfinite`, `isinf`, `isnan`, `signbit`, `copysign`, `nextafter`, `modf`, `ldexp`, `frexp`, `fmod`, `floor`, `ceil`, `trunc`\n", + "\n", + "And the following NumPy functions:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['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', 'nansum', 'nanvar', 'ndim', 'nonzero', 'ones_like', 'pad', 'percentile', 'ptp', 'ravel', 'resize', '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']\n" + ] + } + ], + "source": [ + "from pint.numpy_func import HANDLED_FUNCTIONS\n", + "print(sorted(list(HANDLED_FUNCTIONS)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And the following [NumPy ndarray methods](http://docs.scipy.org/doc/numpy/reference/arrays.ndarray.html#array-methods):\n", + "\n", + "- `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`\n", + "\n", + "Pull requests are welcome for any NumPy function, ufunc, or method that is not currently\n", + "supported.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Array Type Support\n", + "------------------\n", + "\n", + "### Overview\n", + "\n", + "When not wrapping a scalar type, a Pint `Quantity` can be considered a [\"duck array\"](https://numpy.org/neps/nep-0022-ndarray-duck-typing-overview.html), that is, an array-like type that implements (all or most of) NumPy's API for `ndarray`. Many other such duck arrays exist in the Python ecosystem, and Pint aims to work with as many of them as reasonably possible. To date, the following are specifically tested and known to work:\n", + "\n", + "- xarray: `DataArray`, `Dataset`, and `Variable`\n", + "- Sparse: `COO`\n", + "- NumPy: masked arrays\n", + "\n", + "and the following have partial support, with full integration planned:\n", + "\n", + "- Dask arrays\n", + "- CuPy arrays\n", + "\n", + "### Technical Commentary\n", + "\n", + "Starting with version 0.10, Pint aims to interoperate with other duck arrays in a well-defined and well-supported fashion. Part of this support lies in implementing [`__array_ufunc__` to support NumPy ufuncs](https://numpy.org/neps/nep-0013-ufunc-overrides.html) and [`__array_function__` to support NumPy functions](https://numpy.org/neps/nep-0018-array-function-protocol.html). However, the central component to this interoperability is respecting a [type casting hierarchy](https://numpy.org/neps/nep-0018-array-function-protocol.html) of duck arrays. When all types in the hierarchy properly defer to those above it (in wrapping, arithmetic, and NumPy operations), a well-defined nesting and operator precedence order exists. When they don't, the graph of relations becomes cyclic, and the expected result of mixed-type operations becomes ambiguous.\n", + "\n", + "For Pint, following this hierarchy means declaring a list of types that are above it in the hierarchy and to which it defers (\"upcast types\") and assuming all others are below it and wrappable by it (\"downcast types\"). To date, Pint's declared upcast types are:\n", + "\n", + "- `PintArray`, as defined by pint-pandas\n", + "- `Series`, as defined by Pandas\n", + "- `DataArray`, `Dataset`, and `Variable`, as defined by xarray\n", + "\n", + "(Note: if your application requires extension of this collection of types, it is available in Pint's API at `pint.compat.upcast_types`.)\n", + "\n", + "While Pint assumes it can wrap any other duck array (meaning, for now, those that implement `__array_function__`, `shape`, `ndim`, and `dtype`, at least until [NEP 30](https://numpy.org/neps/nep-0030-duck-array-protocol.html) is implemented), there are a few common types that Pint explicitly tests (or plans to test) for optimal interoperability. These are listed above in the overview section and included in the below chart.\n", + "\n", + "This type casting hierarchy of ndarray-like types can be shown by the below acyclic graph, where solid lines represent declared support, and dashed lines represent planned support:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "%3\n", + "\n", + "\n", + "Dask array\n", + "\n", + "Dask array\n", + "\n", + "\n", + "NumPy ndarray\n", + "\n", + "NumPy ndarray\n", + "\n", + "\n", + "Dask array->NumPy ndarray\n", + "\n", + "\n", + "\n", + "\n", + "CuPy ndarray\n", + "\n", + "CuPy ndarray\n", + "\n", + "\n", + "Dask array->CuPy ndarray\n", + "\n", + "\n", + "\n", + "\n", + "Sparse COO\n", + "\n", + "Sparse COO\n", + "\n", + "\n", + "Dask array->Sparse COO\n", + "\n", + "\n", + "\n", + "\n", + "NumPy masked array\n", + "\n", + "NumPy masked array\n", + "\n", + "\n", + "Dask array->NumPy masked array\n", + "\n", + "\n", + "\n", + "\n", + "CuPy ndarray->NumPy ndarray\n", + "\n", + "\n", + "\n", + "\n", + "Sparse COO->NumPy ndarray\n", + "\n", + "\n", + "\n", + "\n", + "NumPy masked array->NumPy ndarray\n", + "\n", + "\n", + "\n", + "\n", + "Jax array\n", + "\n", + "Jax array\n", + "\n", + "\n", + "Jax array->NumPy ndarray\n", + "\n", + "\n", + "\n", + "\n", + "Pint Quantity\n", + "\n", + "Pint Quantity\n", + "\n", + "\n", + "Pint Quantity->Dask array\n", + "\n", + "\n", + "\n", + "\n", + "Pint Quantity->NumPy ndarray\n", + "\n", + "\n", + "\n", + "\n", + "Pint Quantity->CuPy ndarray\n", + "\n", + "\n", + "\n", + "\n", + "Pint Quantity->Sparse COO\n", + "\n", + "\n", + "\n", + "\n", + "xarray Dataset/DataArray/Variable\n", + "\n", + "xarray Dataset/DataArray/Variable\n", + "\n", + "\n", + "xarray Dataset/DataArray/Variable->Dask array\n", + "\n", + "\n", + "\n", + "\n", + "xarray Dataset/DataArray/Variable->NumPy ndarray\n", + "\n", + "\n", + "\n", + "\n", + "xarray Dataset/DataArray/Variable->CuPy ndarray\n", + "\n", + "\n", + "\n", + "\n", + "xarray Dataset/DataArray/Variable->Sparse COO\n", + "\n", + "\n", + "\n", + "\n", + "xarray Dataset/DataArray/Variable->NumPy masked array\n", + "\n", + "\n", + "\n", + "\n", + "xarray Dataset/DataArray/Variable->Jax array\n", + "\n", + "\n", + "\n", + "\n", + "xarray Dataset/DataArray/Variable->Pint Quantity\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from graphviz import Digraph\n", + "\n", + "g = Digraph(graph_attr={'size': '8,5'}, node_attr={'fontname': 'courier'})\n", + "g.edge('Dask array', 'NumPy ndarray')\n", + "g.edge('Dask array', 'CuPy ndarray')\n", + "g.edge('Dask array', 'Sparse COO')\n", + "g.edge('Dask array', 'NumPy masked array')\n", + "g.edge('CuPy ndarray', 'NumPy ndarray')\n", + "g.edge('Sparse COO', 'NumPy ndarray')\n", + "g.edge('NumPy masked array', 'NumPy ndarray')\n", + "g.edge('Jax array', 'NumPy ndarray')\n", + "g.edge('Pint Quantity', 'Dask array', style='dashed')\n", + "g.edge('Pint Quantity', 'NumPy ndarray')\n", + "g.edge('Pint Quantity', 'CuPy ndarray', style='dashed')\n", + "g.edge('Pint Quantity', 'Sparse COO')\n", + "g.edge('xarray Dataset/DataArray/Variable', 'Dask array')\n", + "g.edge('xarray Dataset/DataArray/Variable', 'CuPy ndarray', style='dashed')\n", + "g.edge('xarray Dataset/DataArray/Variable', 'Sparse COO')\n", + "g.edge('xarray Dataset/DataArray/Variable', 'NumPy ndarray')\n", + "g.edge('xarray Dataset/DataArray/Variable', 'NumPy masked array', style='dashed')\n", + "g.edge('xarray Dataset/DataArray/Variable', 'Pint Quantity')\n", + "g.edge('xarray Dataset/DataArray/Variable', 'Jax array', style='dashed')\n", + "g" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Examples\n", + "\n", + "**xarray wrapping Pint Quantity**" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "\n", + "Coordinates:\n", + " * lat (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0\n", + " * lon (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0\n", + " time datetime64[ns] 2013-01-01\n", + "Attributes:\n", + " long_name: 4xDaily Air temperature at sigma level 995\n", + " precision: 2\n", + " GRIB_id: 11\n", + " GRIB_name: TMP\n", + " var_desc: Air temperature\n", + " dataset: NMC Reanalysis\n", + " level_desc: Surface\n", + " statistic: Individual Obs\n", + " parent_stat: Other\n", + " actual_range: [185.16 322.1 ]\n", + "\n", + "\n", + "\n", + "Coordinates:\n", + " time datetime64[ns] 2013-01-01\n" + ] + } + ], + "source": [ + "import xarray as xr\n", + "\n", + "# Load tutorial data\n", + "air = xr.tutorial.load_dataset('air_temperature')['air'][0]\n", + "\n", + "# Convert to Quantity\n", + "air.data = Q_(air.data, air.attrs.pop('units', ''))\n", + "\n", + "print(air)\n", + "print()\n", + "print(air.max())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Pint Quantity wrapping Sparse COO**" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " meter\n", + "\n", + "0.0952163965972451 meter\n" + ] + } + ], + "source": [ + "from sparse import COO\n", + "\n", + "x = np.random.random((100, 100, 100))\n", + "x[x < 0.9] = 0 # fill most of the array with zeros\n", + "s = COO(x)\n", + "\n", + "q = Q_(s, 'm')\n", + "\n", + "print(q)\n", + "print()\n", + "print(np.mean(q))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**xarray wrapping Pint Quantity wrapping Dask array wrapping Sparse COO**" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + ", 'meter')>\n", + "Coordinates:\n", + " * z (z) int64 0 1 2 3 4 5 6 7 8 9 10 ... 90 91 92 93 94 95 96 97 98 99\n", + " * y (y) int64 -50 -49 -48 -47 -46 -45 -44 -43 ... 43 44 45 46 47 48 49\n", + " * x (x) float64 -20.0 -18.5 -17.0 -15.5 ... 124.0 125.5 127.0 128.5\n", + "\n", + "\n", + ", 'meter')>\n", + "Coordinates:\n", + " y int64 -46\n", + " x float64 125.5\n" + ] + } + ], + "source": [ + "import dask.array as da\n", + "\n", + "x = da.random.random((100, 100, 100), chunks=(100, 1, 1))\n", + "x[x < 0.95] = 0\n", + "\n", + "data = xr.DataArray(\n", + " Q_(x.map_blocks(COO), 'm'),\n", + " dims=('z', 'y', 'x'),\n", + " coords={\n", + " 'z': np.arange(100),\n", + " 'y': np.arange(100) - 50,\n", + " 'x': np.arange(100) * 1.5 - 20\n", + " },\n", + " name='test'\n", + ")\n", + "\n", + "print(data)\n", + "print()\n", + "print(data.sel(x=125.5, y=-46).mean())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Compatibility Packages\n", + "\n", + "To aid in integration between various array types and Pint (such as by providing convenience methods), the following compatibility packages are available:\n", + "\n", + "- [pint-pandas](https://github.com/hgrecco/pint-pandas)\n", + "- pint-xarray ([in development](https://github.com/hgrecco/pint/issues/849), initial alpha release planned for January 2020)\n", + "\n", + "(Note: if you have developed a compatibility package for Pint, please submit a pull request to add it to this list!)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Additional Comments\n", + "\n", + "What follows is a short discussion about how NumPy support is implemented in Pint's `Quantity` Object.\n", + "\n", + "For the supported functions, Pint expects certain units and attempts to convert the input (or inputs). For example, the argument of the exponential function (`numpy.exp`) must be dimensionless. Units will be simplified (converting the magnitude appropriately) and `numpy.exp` will be applied to the resulting magnitude. If the input is not dimensionless, a `DimensionalityError` exception will be raised.\n", + "\n", + "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. ndarray or downcast type arguments are generally treated as if they were dimensionless quantities, whereas Pint defers to its declared upcast types by always returning `NotImplemented` when they are encountered (see above).\n", + "\n", + "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\n", + "information on these protocols, see .\n", + "\n", + "This behaviour introduces some performance penalties and increased memory usage. Quantities that must be converted to other units require additional 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 [wrapping](wrapping.html)).\n", + "\n", + "Array interface protocol attributes (such as `__array_struct__` and\n", + "`__array_interface__`) are available on Pint Quantities by deferring to the corresponding `__array_*` attribute on the magnitude as casted to an ndarray. This has been found to be potentially incorrect and to cause unexpected behavior, and has therefore been deprecated. As of the next minor version of Pint (or when the `PINT_ARRAY_PROTOCOL_FALLBACK` environment variable is set to 0 prior to importing Pint as done at the beginning of this page), attempting to access these attributes will instead raise an AttributeError." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.0" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/docs/numpy.rst b/docs/numpy.rst deleted file mode 100644 index c9a46b99c..000000000 --- a/docs/numpy.rst +++ /dev/null @@ -1,181 +0,0 @@ -.. _numpy: - - -NumPy Support -============= - -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: - -.. doctest:: - - >>> import numpy as np - >>> from pint import UnitRegistry - >>> ureg = UnitRegistry() - >>> Q_ = ureg.Quantity - -.. testsetup:: * - - import numpy as np - from pint import UnitRegistry - ureg = UnitRegistry() - Q_ = ureg.Quantity - -and then we create a quantity the standard way - -.. doctest:: - - >>> legs1 = Q_(np.asarray([3., 4.]), 'meter') - >>> print(legs1) - [ 3. 4.] meter - -or we use the property that Pint converts iterables into NumPy ndarrays to simply write: - -.. doctest:: - - >>> legs1 = [3., 4.] * ureg.meter - >>> print(legs1) - [ 3. 4.] meter - -All usual Pint methods can be used with this quantity. For example: - -.. doctest:: - - >>> print(legs1.to('kilometer')) - [ 0.003 0.004] kilometer - >>> print(legs1.dimensionality) - [length] - >>> legs1.to('joule') - Traceback (most recent call last): - ... - DimensionalityError: Cannot convert from 'meter' ([length]) to 'joule' ([length] ** 2 * [mass] / [time] ** 2) - -NumPy functions are supported by Pint. For example if we define: - -.. doctest:: - - >>> legs2 = [400., 300.] * ureg.centimeter - >>> print(legs2) - [ 400. 300.] centimeter - -we can calculate the hypotenuse of the right triangles with legs1 and legs2. - -.. doctest:: - - >>> hyps = np.hypot(legs1, legs2) - >>> print(hyps) - [ 5. 5.] meter - -Notice that before the `np.hypot` was used, the numerical value of legs2 was -internally converted to the units of legs1 as expected. - -Similarly, when you apply a function that expects angles in radians, a conversion -is applied before the requested calculation: - -.. doctest:: - - >>> angles = np.arccos(legs2/hyps) - >>> print(angles) - [ 0.64350111 0.92729522] radian - -You can convert the result to degrees using usual unit conversion: - -.. doctest:: - - >>> print(angles.to('degree')) - [ 36.86989765 53.13010235] degree - -Applying a function that expects angles to a quantity with a different dimensionality -results in an error: - -.. doctest:: - - >>> np.arccos(legs2) - Traceback (most recent call last): - ... - DimensionalityError: Cannot convert from 'centimeter' ([length]) to 'dimensionless' (dimensionless) - - -Support --------- - -The following ufuncs_ can be applied to a Quantity object: - -- **Math operations**: add, subtract, multiply, divide, logaddexp, logaddexp2, true_divide, floor_divide, negative, remainder mod, fmod, absolute, rint, sign, conj, exp, exp2, log, log2, log10, expm1, log1p, sqrt, square, reciprocal -- **Trigonometric functions**: sin, cos, tan, arcsin, arccos, arctan, arctan2, hypot, sinh, cosh, tanh, arcsinh, arccosh, arctanh, deg2rad, rad2deg -- **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 functions: - -- 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, pad, percentile, ptp, ravel, resize, 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 - -And the following `NumPy ndarray methods`_: - -- 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 --------- - -What follows is a short discussion about how NumPy support is implemented in -Pint's `Quantity` Object. - -For the supported functions, Pint expects certain units and attempts to convert -the input (or inputs). For example, the argument of the exponential function -(`numpy.exp`) must be dimensionless. Units will be simplified (converting the -magnitude appropriately) and `numpy.exp` will be applied to the resulting -magnitude. If the input is not dimensionless, a `DimensionalityError` exception -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. ndarray or ndarray-like arguments -are generally treated as if they were dimensionless quantities, except for declared -upcast types to which Pint defers (see -). To date, these "upcast types" are: - -- ``PintArray``, as defined by pint-pandas -- ``Series``, as defined by pandas -- ``DataArray``, ``Dataset``, and ``Variable``, as defined by xarray - -If your application requires extension of this collection of types, it is available in -Pint's API at ``pint.compat.upcast_types``. Note that these are also the types to which -a Quantity object will defer for arithmetic operations. - -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 -. - -This behaviour introduces some performance penalties and increased memory -usage. Quantities that must be converted to other units require additional -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`). - -Array interface protocol attributes (such as `__array_struct__` and -`__array_interface__`) are available on Pint Quantities by deferring to the -corresponding `__array_*` attribute on the magnitude as casted to an ndarray. This -has been found to be potentially incorrect and to cause unexpected behavior, and has -therefore been deprecated. As of the next minor version of Pint (or when the -`PINT_ARRAY_PROTOCOL_FALLBACK` environment variable is set to 0 prior to importing -Pint), attempting to access these attributes will instead raise an AttributeError. - - - - - -.. _`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 e87c26c23..6fb27ee5a 100644 --- a/pint/compat.py +++ b/pint/compat.py @@ -66,14 +66,16 @@ class BehaviorChangeWarning(UserWarning): NUMPY_VER = np.__version__ NUMERIC_TYPES = (Number, Decimal, ndarray, np.number) - def _to_magnitude(value, force_ndarray=False): + def _to_magnitude(value, force_ndarray=False, force_ndarray_like=False): if isinstance(value, (dict, bool)) or value is None: raise TypeError("Invalid magnitude for Quantity: {0!r}".format(value)) elif isinstance(value, str) and value == "": raise ValueError("Quantity magnitude cannot be an empty string.") elif isinstance(value, (list, tuple)): return np.asarray(value) - if force_ndarray: + if force_ndarray or ( + force_ndarray_like and not is_duck_array_type(type(value)) + ): return np.asarray(value) return value @@ -112,9 +114,11 @@ class ndarray: NP_NO_VALUE = None ARRAY_FALLBACK = False - def _to_magnitude(value, force_ndarray=False): - if force_ndarray: - raise ValueError("Cannot force to ndarray when NumPy is not present.") + def _to_magnitude(value, force_ndarray=False, force_ndarray_like=False): + if force_ndarray or force_ndarray_like: + raise ValueError( + "Cannot force to ndarray or ndarray-like when NumPy is not present." + ) elif isinstance(value, (dict, bool)) or value is None: raise TypeError("Invalid magnitude for Quantity: {0!r}".format(value)) elif isinstance(value, str) and value == "": @@ -148,8 +152,8 @@ def _to_magnitude(value, force_ndarray=False): if not HAS_BABEL: babel_parse = babel_units = missing_dependency("Babel") # noqa: F811 -# Define location of pint.Quantity in NEP-13 type cast hierarchy by defining upcast and -# downcast/wrappable types using guarded imports +# Define location of pint.Quantity in NEP-13 type cast hierarchy by defining upcast +# types using guarded imports upcast_types = [] # pint-pandas (PintArray) @@ -191,6 +195,28 @@ def is_upcast_type(other): return other in upcast_types +def is_duck_array_type(other): + """Check if the type object represents a (non-Quantity) duck array type. + + Parameters + ---------- + other : object + + Returns + ------- + bool + """ + # TODO (NEP 30): replace duck array check with hasattr(other, "__duckarray__") + return other is ndarray or ( + not hasattr(other, "_magnitude") + and not hasattr(other, "_units") + and HAS_NUMPY_ARRAY_FUNCTION + and hasattr(other, "__array_function__") + and hasattr(other, "ndim") + and hasattr(other, "dtype") + ) + + def eq(lhs, rhs, check_all): """Comparison of scalars and arrays. diff --git a/pint/numpy_func.py b/pint/numpy_func.py index 0dcb878d7..1a9221caa 100644 --- a/pint/numpy_func.py +++ b/pint/numpy_func.py @@ -424,7 +424,7 @@ def implementation(*args, **kwargs): 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 + # Require all inputs to match units, but output base ndarray/duck array 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(): @@ -744,6 +744,7 @@ def implementation(*args, **kwargs): ("rot90", "m", True), ("insert", ["arr", "values"], True), ("resize", "a", True), + ("reshape", "a", True), ]: implement_consistent_units_by_argument(func_str, unit_arguments, wrap_output) diff --git a/pint/quantity.py b/pint/quantity.py index 91abfb77f..ba4a26431 100644 --- a/pint/quantity.py +++ b/pint/quantity.py @@ -29,6 +29,7 @@ array_function_change_msg, babel_parse, eq, + is_duck_array_type, is_upcast_type, ndarray, np, @@ -145,6 +146,10 @@ class Quantity(PrettyIPython, SharedRegistryObject): def force_ndarray(self): return self._REGISTRY.force_ndarray + @property + def force_ndarray_like(self): + return self._REGISTRY.force_ndarray_like + def __reduce__(self): """Allow pickling quantities. Since UnitRegistries are not pickled, upon unpickling the new object is always attached to the application registry. @@ -173,15 +178,21 @@ def __new__(cls, value, units=None): inst = copy.copy(value) else: inst = SharedRegistryObject.__new__(cls) - inst._magnitude = _to_magnitude(value, inst.force_ndarray) + inst._magnitude = _to_magnitude( + value, inst.force_ndarray, inst.force_ndarray_like + ) inst._units = UnitsContainer() elif isinstance(units, (UnitsContainer, UnitDefinition)): inst = SharedRegistryObject.__new__(cls) - inst._magnitude = _to_magnitude(value, inst.force_ndarray) + inst._magnitude = _to_magnitude( + value, inst.force_ndarray, inst.force_ndarray_like + ) inst._units = units elif isinstance(units, str): inst = SharedRegistryObject.__new__(cls) - inst._magnitude = _to_magnitude(value, inst.force_ndarray) + inst._magnitude = _to_magnitude( + value, inst.force_ndarray, inst.force_ndarray_like + ) inst._units = inst._REGISTRY.parse_units(units)._units elif isinstance(units, SharedRegistryObject): if isinstance(units, Quantity) and units.magnitude != 1: @@ -192,7 +203,9 @@ def __new__(cls, value, units=None): else: inst = SharedRegistryObject.__new__(cls) inst._units = units._units - inst._magnitude = _to_magnitude(value, inst.force_ndarray) + inst._magnitude = _to_magnitude( + value, inst.force_ndarray, inst.force_ndarray_like + ) else: raise TypeError( "units must be of type str, Quantity or " @@ -502,7 +515,7 @@ def _convert_magnitude(self, other, *contexts, **ctx_kwargs): self._magnitude, self._units, other, - inplace=isinstance(self._magnitude, ndarray), + inplace=is_duck_array_type(type(self._magnitude)), ) def ito(self, other=None, *contexts, **ctx_kwargs): @@ -729,7 +742,9 @@ def _iadd_sub(self, other, op): if not self._check(other): # other not from same Registry or not a Quantity try: - other_magnitude = _to_magnitude(other, self.force_ndarray) + other_magnitude = _to_magnitude( + other, self.force_ndarray, self.force_ndarray_like + ) except PintTypeError: raise except TypeError: @@ -848,12 +863,14 @@ def _add_sub(self, other, op): # the operation. units = self._units magnitude = op( - self._magnitude, _to_magnitude(other, self.force_ndarray) + self._magnitude, + _to_magnitude(other, self.force_ndarray, self.force_ndarray_like), ) elif self.dimensionless: units = UnitsContainer() magnitude = op( - self.to(units)._magnitude, _to_magnitude(other, self.force_ndarray) + self.to(units)._magnitude, + _to_magnitude(other, self.force_ndarray, self.force_ndarray_like), ) else: raise DimensionalityError(self._units, "dimensionless") @@ -939,10 +956,10 @@ def _add_sub(self, other, op): def __iadd__(self, other): if isinstance(other, datetime.datetime): return self.to_timedelta() + other - elif not isinstance(self._magnitude, ndarray): - return self._add_sub(other, operator.add) - else: + elif is_duck_array_type(type(self._magnitude)): return self._iadd_sub(other, operator.iadd) + else: + return self._add_sub(other, operator.add) def __add__(self, other): if isinstance(other, datetime.datetime): @@ -953,10 +970,10 @@ def __add__(self, other): __radd__ = __add__ def __isub__(self, other): - if not isinstance(self._magnitude, ndarray): - return self._add_sub(other, operator.sub) - else: + if is_duck_array_type(type(self._magnitude)): return self._iadd_sub(other, operator.isub) + else: + return self._add_sub(other, operator.sub) def __sub__(self, other): return self._add_sub(other, operator.sub) @@ -1007,7 +1024,9 @@ def _imul_div(self, other, magnitude_op, units_op=None): self._units, getattr(other, "units", "") ) try: - other_magnitude = _to_magnitude(other, self.force_ndarray) + other_magnitude = _to_magnitude( + other, self.force_ndarray, self.force_ndarray_like + ) except PintTypeError: raise except TypeError: @@ -1075,7 +1094,9 @@ def _mul_div(self, other, magnitude_op, units_op=None): self._units, getattr(other, "units", "") ) try: - other_magnitude = _to_magnitude(other, self.force_ndarray) + other_magnitude = _to_magnitude( + other, self.force_ndarray, self.force_ndarray_like + ) except PintTypeError: raise except TypeError: @@ -1109,10 +1130,10 @@ def _mul_div(self, other, magnitude_op, units_op=None): return self.__class__(magnitude, units) def __imul__(self, other): - if not isinstance(self._magnitude, ndarray): - return self._mul_div(other, operator.mul) - else: + if is_duck_array_type(type(self._magnitude)): return self._imul_div(other, operator.imul) + else: + return self._mul_div(other, operator.mul) def __mul__(self, other): return self._mul_div(other, operator.mul) @@ -1129,17 +1150,19 @@ def __matmul__(self, other): __rmatmul__ = __matmul__ def __itruediv__(self, other): - if not isinstance(self._magnitude, ndarray): - return self._mul_div(other, operator.truediv) - else: + if is_duck_array_type(type(self._magnitude)): return self._imul_div(other, operator.itruediv) + else: + return self._mul_div(other, operator.truediv) def __truediv__(self, other): return self._mul_div(other, operator.truediv) def __rtruediv__(self, other): try: - other_magnitude = _to_magnitude(other, self.force_ndarray) + other_magnitude = _to_magnitude( + other, self.force_ndarray, self.force_ndarray_like + ) except PintTypeError: raise except TypeError: @@ -1233,11 +1256,11 @@ def __rdivmod__(self, other): @check_implemented def __ipow__(self, other): - if not isinstance(self._magnitude, ndarray): + if not is_duck_array_type(type(self._magnitude)): return self.__pow__(other) try: - _to_magnitude(other, self.force_ndarray) + _to_magnitude(other, self.force_ndarray, self.force_ndarray_like) except PintTypeError: raise except TypeError: @@ -1246,7 +1269,7 @@ def __ipow__(self, other): if not self._ok_for_muldiv: raise OffsetUnitCalculusError(self._units) - if isinstance(getattr(other, "_magnitude", other), ndarray): + if is_duck_array_type(type(getattr(other, "_magnitude", other))): # arrays are refused as exponent, because they would create # len(array) quantities of len(set(array)) different units # unless the base is dimensionless. @@ -1286,13 +1309,15 @@ def __ipow__(self, other): else: self._units **= other - self._magnitude **= _to_magnitude(other, self.force_ndarray) + self._magnitude **= _to_magnitude( + other, self.force_ndarray, self.force_ndarray_like + ) return self @check_implemented def __pow__(self, other): try: - _to_magnitude(other, self.force_ndarray) + _to_magnitude(other, self.force_ndarray, self.force_ndarray_like) except PintTypeError: raise except TypeError: @@ -1301,7 +1326,7 @@ def __pow__(self, other): if not self._ok_for_muldiv: raise OffsetUnitCalculusError(self._units) - if isinstance(getattr(other, "_magnitude", other), ndarray): + if is_duck_array_type(type(getattr(other, "_magnitude", other))): # arrays are refused as exponent, because they would create # len(array) quantities of len(set(array)) different units # unless the base is dimensionless. @@ -1339,7 +1364,9 @@ def __pow__(self, other): elif not getattr(other, "dimensionless", True): raise DimensionalityError(other._units, "dimensionless") else: - exponent = _to_magnitude(other, self.force_ndarray) + exponent = _to_magnitude( + other, self.force_ndarray, self.force_ndarray_like + ) units = new_self._units ** exponent magnitude = new_self._magnitude ** exponent @@ -1348,7 +1375,7 @@ def __pow__(self, other): @check_implemented def __rpow__(self, other): try: - _to_magnitude(other, self.force_ndarray) + _to_magnitude(other, self.force_ndarray, self.force_ndarray_like) except PintTypeError: raise except TypeError: @@ -1356,7 +1383,7 @@ def __rpow__(self, other): else: if not self.dimensionless: raise DimensionalityError(self._units, "dimensionless") - if isinstance(self._magnitude, ndarray): + if is_duck_array_type(type(self._magnitude)): if np.size(self._magnitude) > 1: raise DimensionalityError(self._units, "dimensionless") new_self = self.to_root_units() @@ -1415,7 +1442,7 @@ def __eq__(self, other): @check_implemented def __ne__(self, other): out = self.__eq__(other) - if isinstance(out, ndarray): + if is_duck_array_type(type(out)): return np.logical_not(out) return not out @@ -1637,12 +1664,12 @@ def __getattr__(self, item): stacklevel=2, ) - if isinstance(self._magnitude, ndarray): + if is_duck_array_type(type(self._magnitude)): + # Defer to magnitude, and don't catch any AttributeErrors 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. - # TODO (#905 follow-up): Potentially problematic, investigate for duck arrays + # If an `__array_` attribute is requested but the magnitude is not + # a duck array, we convert the magnitude to a numpy ndarray. magnitude_as_array = _to_magnitude( self._magnitude, force_ndarray=True ) @@ -1651,13 +1678,23 @@ def __getattr__(self, item): # TODO (next minor version): ARRAY_FALLBACK is removed and this becomes the standard behavior raise AttributeError(f"Array protocol attribute {item} not available.") 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): + magnitude_as_duck_array = _to_magnitude( + self._magnitude, force_ndarray_like=True + ) + try: + attr = getattr(magnitude_as_duck_array, item) return functools.partial(self._numpy_method_wrap, attr) - else: - raise AttributeError("NumPy method {} was not callable.".format(item)) + except AttributeError: + raise AttributeError( + f"NumPy method {item} not available on {type(magnitude_as_duck_array)}" + ) + except TypeError as exc: + if "not callable" in str(exc): + raise AttributeError( + f"NumPy method {item} not callable on {type(magnitude_as_duck_array)}" + ) + else: + raise exc try: return getattr(self._magnitude, item) diff --git a/pint/registry.py b/pint/registry.py index cdb85dd93..90d595be6 100644 --- a/pint/registry.py +++ b/pint/registry.py @@ -140,6 +140,8 @@ class BaseRegistry(metaclass=RegistryMeta): the default definition file. None to leave the UnitRegistry empty. force_ndarray : bool convert any input, scalar or not to a numpy.ndarray. + force_ndarray_like : bool + convert all inputs other than duck arrays to a numpy.ndarray. on_redefinition : str action to take in case a unit is redefined: 'warn', 'raise', 'ignore' auto_reduce_dimensions : @@ -179,6 +181,7 @@ def __init__( self, filename="", force_ndarray=False, + force_ndarray_like=False, on_redefinition="warn", auto_reduce_dimensions=False, preprocessors=None, @@ -189,6 +192,7 @@ def __init__( self._filename = filename self.force_ndarray = force_ndarray + self.force_ndarray_like = force_ndarray_like self.preprocessors = preprocessors or [] #: Action to take in case a unit is redefined. 'warn', 'raise', 'ignore' @@ -1950,8 +1954,10 @@ class UnitRegistry(SystemRegistry, ContextRegistry, NonMultiplicativeRegistry): path of the units definition file to load or line-iterable object. Empty to load the default definition file. None to leave the UnitRegistry empty. - force_ndarray : + force_ndarray : bool convert any input, scalar or not to a numpy.ndarray. + force_ndarray_like : bool + convert all inputs other than duck arrays to a numpy.ndarray. default_as_delta : In the context of a multiplication of units, interpret non-multiplicative units as their *delta* counterparts. @@ -1975,6 +1981,7 @@ def __init__( self, filename="", force_ndarray=False, + force_ndarray_like=False, default_as_delta=True, autoconvert_offset_to_baseunit=False, on_redefinition="warn", @@ -1987,6 +1994,7 @@ def __init__( super().__init__( filename=filename, force_ndarray=force_ndarray, + force_ndarray_like=force_ndarray_like, on_redefinition=on_redefinition, default_as_delta=default_as_delta, autoconvert_offset_to_baseunit=autoconvert_offset_to_baseunit, diff --git a/pint/testsuite/test_compat_downcast.py b/pint/testsuite/test_compat_downcast.py new file mode 100644 index 000000000..e85369a14 --- /dev/null +++ b/pint/testsuite/test_compat_downcast.py @@ -0,0 +1,108 @@ +import pytest + +from pint import UnitRegistry + +# Conditionally import NumPy and any upcast type libraries +np = pytest.importorskip("numpy", reason="NumPy is not available") +sparse = pytest.importorskip("sparse", reason="sparse is not available") + +# Set up unit registry and sample +ureg = UnitRegistry(force_ndarray_like=True) +q_base = (np.arange(25).reshape(5, 5).T + 1) * ureg.kg + + +# Define identity function for use in tests +def identity(x): + return x + + +@pytest.fixture(params=["sparse", "masked_array"]) +def array(request): + """Generate 5x5 arrays of given type for tests.""" + if request.param == "sparse": + # Create sample sparse COO as a permutation matrix. + coords = [[0, 1, 2, 3, 4], [1, 3, 0, 2, 4]] + data = [1.0] * 5 + return sparse.COO(coords, data, shape=(5, 5)) + elif request.param == "masked_array": + # Create sample masked array as an upper triangular matrix. + return np.ma.masked_array( + np.arange(25, dtype=np.float).reshape((5, 5)), + mask=np.logical_not(np.triu(np.ones((5, 5)))), + ) + + +@pytest.mark.parametrize( + "op, magnitude_op, unit_op", + [ + pytest.param(identity, identity, identity, id="identity"), + pytest.param( + lambda x: x + 1 * ureg.m, lambda x: x + 1, identity, id="addition" + ), + pytest.param( + lambda x: x - 20 * ureg.cm, lambda x: x - 0.2, identity, id="subtraction" + ), + pytest.param( + lambda x: x * (2 * ureg.s), + lambda x: 2 * x, + lambda u: u * ureg.s, + id="multiplication", + ), + pytest.param( + lambda x: x / (1 * ureg.s), identity, lambda u: u / ureg.s, id="division" + ), + pytest.param(lambda x: x ** 2, lambda x: x ** 2, lambda u: u ** 2, id="square"), + pytest.param(lambda x: x.T, lambda x: x.T, identity, id="transpose"), + pytest.param(np.mean, np.mean, identity, id="mean ufunc"), + pytest.param(np.sum, np.sum, identity, id="sum ufunc"), + pytest.param(np.sqrt, np.sqrt, lambda u: u ** 0.5, id="sqrt ufunc"), + pytest.param( + lambda x: np.reshape(x, 25), + lambda x: np.reshape(x, 25), + identity, + id="reshape function", + ), + pytest.param(np.amax, np.amax, identity, id="amax function"), + ], +) +def test_univariate_op_consistency(op, magnitude_op, unit_op, array): + q = ureg.Quantity(array, "meter") + res = op(q) + assert np.all(res.magnitude == magnitude_op(array)) # Magnitude check + assert res.units == unit_op(q.units) # Unit check + assert q.magnitude is array # Immutability check + + +@pytest.mark.parametrize( + "op, unit", + [ + pytest.param(lambda x, y: x * y, ureg("kg m"), id="multiplication"), + pytest.param(lambda x, y: x / y, ureg("m / kg"), id="division"), + pytest.param(np.multiply, ureg("kg m"), id="multiply ufunc"), + ], +) +def test_bivariate_op_consistency(op, unit, array): + q = ureg.Quantity(array, "meter") + res = op(q, q_base) + assert np.all(res.magnitude == op(array, q_base.magnitude)) # Magnitude check + assert res.units == unit # Unit check + assert q.magnitude is array # Immutability check + + +@pytest.mark.parametrize( + "op", + [ + pytest.param( + lambda a, u: a * u, + id="array-first", + marks=pytest.mark.xfail(reason="upstream issue numpy/numpy#15200"), + ), + pytest.param(lambda a, u: u * a, id="unit-first"), + ], +) +@pytest.mark.parametrize( + "unit", + [pytest.param(ureg.m, id="Unit"), pytest.param(ureg("meter"), id="Quantity")], +) +def test_array_quantity_creation_by_multiplication(op, unit, array): + assert type(op(array, unit)) == ureg.Quantity diff --git a/pint/unit.py b/pint/unit.py index 10ebd71e2..b61cfe6fa 100644 --- a/pint/unit.py +++ b/pint/unit.py @@ -13,7 +13,7 @@ import operator from numbers import Number -from .compat import NUMERIC_TYPES +from .compat import NUMERIC_TYPES, is_upcast_type from .definitions import UnitDefinition from .formatting import siunitx_format_unit from .util import PrettyIPython, SharedRegistryObject, UnitsContainer @@ -230,18 +230,32 @@ def __complex__(self): __array_priority__ = 17 - def __array_prepare__(self, array, context=None): - return 1 + def __array_ufunc__(self, ufunc, method, *inputs, **kwargs): + if method != "__call__": + # Only handle ufuncs as callables + return NotImplemented - def __array_wrap__(self, array, context=None): - uf, objs, huh = context + # Check types and return NotImplemented when upcast type encountered + types = set( + type(arg) + for arg in list(inputs) + list(kwargs.values()) + if hasattr(arg, "__array_ufunc__") + ) + if any(is_upcast_type(other) for other in types): + return NotImplemented - if uf.__name__ in ("true_divide", "divide", "floor_divide"): - return self._REGISTRY.Quantity(array, 1 / self._units) - elif uf.__name__ in ("multiply",): - return self._REGISTRY.Quantity(array, self._units) + # Act on limited implementations by conversion to multiplicative identity + # Quantity + if ufunc.__name__ in ("true_divide", "divide", "floor_divide", "multiply"): + return ufunc( + *tuple( + self._REGISTRY.Quantity(1, self._units) if arg is self else arg + for arg in inputs + ), + **kwargs, + ) else: - raise ValueError("Unsupproted operation for Unit") + return NotImplemented @property def systems(self):