v2.14.0

This release fixes several results that were silently wrong. Each fix changes behavior that existing code may rely on, so each section ends with the changes to check before upgrading. The old results were wrong rather than a different design, so this is a minor release.

Uncertain masks and sorting

Indexing an uncertain array with an uncertain boolean mask used to keep only the elements chosen by the nominal value and by every sample of the distribution. That biased everything computed from the selection, the nominal value included: summing the pixels above a threshold, in the example in #252, gave 414 where masking each realization separately gives 864. Assignment through such a mask (a[a < 0] = 0) and numpy.nonzero used the same rule. numpy.sort and numpy.argsort ordered whole elements by the mean of their distribution, so the first element of a sorted array was not its minimum (#252).

Every operation that can work sample by sample now does:

• a[m] keeps every element chosen by the nominal value or by any sample, and holds NaN in each realization that didn't choose it, as xarray.DataArray.where(cond, drop=True) does. NaN-aware reductions of the result then match masking each realization separately: np.nansum(a[m]) gives the same as np.sum(a, where=m). A mask whose samples all agree selects exactly as a plain mask does, and pads nothing.
• a[m] = v writes each realization only where it chose an element, so a[m] = f(a[m]) updates every realization correctly.
• np.sort, np.argsort and na.stats.rankdata work on the nominal value and on each sample separately.
• Reductions of a plain array accept an uncertain where, and reduce each sample over its own selection.
• Indexing a ScalarArray with an index array that varies along an axis the item doesn't name gathers along that axis, for reading and for writing, instead of raising "Each axis name must be unique". So assigning through the indices from np.argsort scatters correctly.

The docstring of UncertainScalarArray sets out these rules, with examples.

Behavior changes:

• a[m] with an uncertain mask returns more elements, padded with NaN. A reduction that isn't NaN-aware gives NaN for the realizations that disagree with the union, so use the nan reductions or where=. A NaN-aware reduction warns about an empty slice for a realization that selects nothing.
• Booleans are padded with False, which keeps np.any and np.sum right in every sample but not np.all. Use np.all(b, where=m) for that.
• Uncertain integer and string arrays raise a ValueError when indexed by a mask whose samples select different elements, since no value of theirs can mean "not selected". So does np.nonzero of such a mask, and the error suggests alternatives.
• A plain floating-point or boolean array indexed by such a mask is padded and becomes uncertain. Indexed by a mask whose samples all agree, it now stays plain, where it used to become uncertain.
• plain[m] = v raises a ValueError when the samples of m differ, where it used to do nothing. A plain component of a vector, or a plain input or output of a function, is replaced by an uncertain copy when it is assigned through such a mask or assigned an uncertain value.
• Assigning to a scalar with an unsupported item, such as a tuple, raises a TypeError, where it used to do nothing.
• np.argsort of an uncertain array returns uncertain indices, so reordering plain data by them, such as a list of labels, gives a different order in each sample. Sort by na.nominal(...) where only a display order is wanted.
• a[dict(x=idx)] = v, with idx varying along another axis y of a, sets a[idx[j], j] for each j. It used to set the whole rows a[idx, :].

Extra fields of function array subclasses

Subclasses of FunctionArray often carry fields that are sampled along the axes of the function, such as the exposure time of each image. Every operation passed those fields on unchanged, so anything that selected or rearranged elements left them stale: obs[dict(time=0)].timedelta still had an entry for every frame (#263).

A field whose value is a named array now follows the outputs through indexing (including index arrays and boolean masks), np.stack, np.concatenate, np.moveaxis, np.repeat, np.take_along_axis, np.reshape, combine_axes, cell_centers, na.debroadcast, assignment and np.copyto. Every other field is passed on unchanged.

A function array checks, when it is built, that each such field has as many elements as the outputs along the axes they share. So an operation that changes the outputs without the field fails at once, rather than leaving the field stale. A subclass with its own __post_init__ must call super().__post_init__() for the check to run. The docstring of FunctionArray describes the rules and the limits of the check.

Also fixed: np.matmul(..., out=...) on a function array, which crashed.

Behavior changes, for subclasses whose fields are named arrays:

• An operation that removes or resamples an axis a field varies along raises a ValueError naming the field. These are keepdims=False reductions, percentile and quantile, integrate, regrid and cell_centers(random=True). So does a keepdims=True reduction (the default for a function array) along an axis the inputs don't have, since the result keeps one element along it. Reduce the field first, the way it should be reduced:

  total = np.sum(
      dataclasses.replace(images, timedelta=images.timedelta.sum("time")),
      axis="time",
      keepdims=False,
  )

• cell_centers raises for a field whose elements along the axis aren't floating point, such as a mask.

• Stacking or concatenating arrays of two unrelated subclasses raises a TypeError, rather than dropping a field. A field that varies can't be combined with an array that lacks it. A field that isn't a named array, such as a label, must be the same in every array.

• Assignment and np.copyto write the fields in place, like the outputs. The result of a ufunc, to() or a reduction shares the fields of the original, so assigning into it changes the original's fields too. A field that doesn't vary along the assigned axes, or isn't a named array, can only be given the value it already has.

Indexing function arrays

Indexing the cells along a vertex axis with a negative or reversed index kept the wrong vertices: f[dict(x=-1)] had no inputs at all. The index is now normalized against the number of cells, so f[dict(x=-1)] is the last cell and its two vertices, and a reversed slice keeps the vertices of the cells it selects (#261). A step of -1 still works, so np.flip does too.

Indexing a function array whose inputs are a WCS vector no longer computes every coordinate (#262). The inputs stay a WCS vector, with the reference pixel and shape adjusted to the crop, and the shape of a WCS vector now comes from at most two pixels along each axis. Cropping 10 × 10 pixels from six 4096 × 4096 images went from 1.6 s and 4.8 GB to 5 ms and under 1 MB. An item it can't handle this way, such as a stepped slice along a WCS axis, still computes the coordinates first.

Behavior changes:

• A slice along a vertex axis with a step other than 1 or -1 raises a ValueError if it selects more than one cell, since cells that aren't adjacent don't share vertices. It used to return every other cell with every other vertex, which only matched when the step divided the number of cells.
• An integer past the last cell raises an IndexError, where it used to return one vertex and no cells.
• Assignment along a vertex axis writes only the outputs, since the value can't say where the cell edges are. So f[dict(x=-1)] = g[dict(x=-1)] now works, where it raised.
• The coordinates of cropped WCS inputs are computed each time they are accessed, as they always were for uncropped ones, so code that reads them repeatedly should call .explicit once. Writing to them changes a temporary copy, and their parameters are views of the parent's.

Points in polygons

na.geometry.point_in_polygon() copied the vertices for every point whenever the points were in a different unit than the vertices, or the polygon varied along an axis of the points. For an outline of 81 vertices, that took 1.3 kB a point. It now works out each point's polygon from the point's position in the grid instead (#265). For 0.79 million points against a polygon that moves with wavelength, with the points in a different unit, the peak memory went from 1471 MB to 13 MB, and the time from 0.40 s to 0.012 s. Its inputs always have the same layout, so Numba compiles it once, rather than once for each layout, dtype and read-only flag.

Behavior changes:

• If any argument has a unit, every argument is converted to the first such unit, and an argument without one counts as dimensionless. Before, the unit came from x alone, so x=0, y=1.5 cm came out inside a circle of radius 10 mm. Now a coordinate without a unit, compared with lengths, raises a UnitConversionError.
• An axis that isn't an axis of the vertices raises a ValueError naming their axes, where it used to raise a KeyError.

Other changes

• Remove about half the fixed cost of an operation on a small array (#258). Adding arrays of 10 × 3 and 10 elements takes 7 µs instead of 15 µs, and indexing by a boolean mask takes 5 µs instead of 15 µs. Large arrays are unaffected.
• Allow indexing a ScalarArray of astropy.time.Time by a boolean mask, which raised an IndexError (#258).

Full Changelog: v2.13.0...v2.14.0