Arrays and indexing

Arrays and views, leading-axis broadcasting, names, keys and Python interoperability.

import rustymath as rm

Python operations run eagerly. Slices, reversals, transposes, padding and reshapes share storage without writable aliases: assigning to an array never changes another array’s values. A view keeps its whole base allocation alive. Use a.compact() to keep just its logical items, even for a contiguous slice of a larger base. Float display rounding depends on the dtype and does not change the values returned by .tolist(). For expression fusion and richer array notation, use BasedPL.

Construction and dtypes

array copies Python values or a buffer exporter into an owned array. dtype selects its item type: bool, uint8, int16, int32, int64, float16, float32, float64, complex64, complex128 or character. Fill constructors take a shape. Their _like forms retain another array’s shape, keys and names:

temperatures = rm.array([[18., 20, 22], [16, 19, 23]], dtype='float32')
temperatures, rm.zeros_like(temperatures), rm.full((2, 3), 1)

(Array([[18.0, 20.0, 22.0], [16.0, 19.0, 23.0]], dtype=float32),
 Array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], dtype=float32),
 Array([[1, 1, 1], [1, 1, 1]], dtype=int64))

A string becomes an array of individual Unicode characters, not variable-length strings. .astype() keeps keys and names. An integer target converts as BPL’s ₓ does: a float takes its floor and saturates at the integer’s range, NaN giving 0, and an integer wraps to the width. Casts between characters and numbers are not supported.

rm.array('東京∑'), temperatures.astype('int16')

(Array(['東', '京', '∑'], dtype=character),
 Array([[18, 20, 22], [16, 19, 23]], dtype=int16))

arange generates values before its stop. linspace includes both endpoints by default. logspace generates evenly spaced exponents in the chosen base:

rm.arange(0, 6, 2), rm.linspace(0, 1, 5), rm.logspace(0, 3, 4)

(Array([0, 2, 4], dtype=int64),
 Array([0.0, 0.25, 0.5, 0.75, 1.0], dtype=float64),
 Array([1.0, 10.0, 100.0, 1000.0], dtype=float64))

Dtype helpers use the same supported item types as arrays. finfo reports float limits and spacing at 1; iinfo reports integer limits. result_type gives the shared promotion rule:

rm.finfo('float32').eps, rm.iinfo('int16').max, rm.result_type('int16', 'float16'), rm.issubdtype('float32', 'floating')

(1.1920928955078125e-07, 32767, 'float16', True)

Broadcasting and conversion

Each function takes numbers, nested lists, rustymath arrays or buffer exporters such as NumPy arrays. Buffers need a supported item format; foreign byte order is converted. A function returns a rustymath.Array, or a Python number when every argument is a Python number. Arguments broadcast by aligning their leading axes, where NumPy aligns trailing axes. An axis of length 1 stretches to match. broadcast_to returns a compact copy at the requested shape. Here a vector of two items adds one number to each row of a 2 by 3 matrix:

rm.add([[1, 2, 3], [4, 5, 6]], [10, 20])

Array([[11, 12, 13], [24, 25, 26]], dtype=int64)

Arrays of different dtypes promote to the higher kind, wide enough to hold both. An integer array beside a float or complex array takes that dtype, as in JAX. So an int16 array beside a float16 array gives float16, where NumPy gives float32. Characters don’t mix with numbers. An integer sum, difference, product, magnitude or power that doesn’t fit its dtype wraps, as in NumPy. A Python number or list follows NEP 50, as in NumPy: it takes the dtype of the array beside it, unless its kind is higher, so a float beside an integer array gives float64. A function that computes only on floats, such as exp, converts Booleans and integers to the narrowest float that holds them exactly, but at least float32, which computes faster than float16. So Booleans, uint8 and int16 give float32, and wider integers give float64:

i16 = rm.array([1, 2]).astype('int16')
rm.add(i16, 3), rm.add(i16, 0.5), rm.add(i16, rm.array([0.5, 1]).astype('float16')), rm.exp(i16), rm.exp([1, 2])

(Array([4, 5], dtype=int16),
 Array([1.5, 2.5], dtype=float64),
 Array([1.5, 3.0], dtype=float16),
 Array([2.71828, 7.38906], dtype=float32),
 Array([2.7182818, 7.3890561], dtype=float64))

A rustymath.Array has arithmetic, bitwise and comparison operators, @, and the methods .sum(), .prod(), .max(), .min(), .all(), .any(), .cumsum(), .cumprod(), .reshape(), .astype() and .tolist(). With no axis, a reduction reduces the whole array to a number, and a scan runs over the flattened array:

m = rm.array([[1, 2], [3, 4]])
(m * 10).sum(), m.sum(axis=0), m.cumsum()

(100, Array([4, 6], dtype=int64), Array([1, 3, 6, 10], dtype=int64))

Selecting and reshaping

Brackets take one selector for each axis. Index arrays pair by broadcasting from the front, as function arguments do: grid[[1, 0], [2, 0]] reads two items. Length-1 axes expand, so index arrays shaped as an open grid read every combination, as np.ix_’s do. The paired axes take the first array’s place, where NumPy moves them to the front when a slice separates the arrays. Integers remove an axis; lists keep it. Booleans are positions 0 and 1, as in BPL’s indexing; repeat with a mask as its counts keeps the items where the mask is true. Negative positions count from the end, and slices take Python’s start, stop and step syntax:

grid = rm.array([[10, 20, 30], [40, 50, 60]])
mask = [True, False, True]
assert grid[:, mask].tolist() == [[20, 10, 20], [50, 40, 50]]
grid[[1, 0], [2, 0]], grid[[[1], [0]], [[2, 0]]], grid[:, ::-1], grid[:, mask], rm.repeat(grid, mask, axis=1)

(Array([60, 10], dtype=int64),
 Array([[60, 40], [30, 10]], dtype=int64),
 Array([[30, 20, 10], [60, 50, 40]], dtype=int64),
 Array([[10, 30], [40, 60]], dtype=int64))

len, iteration and in use major cells: items of a vector or rows of a matrix. Membership compares the whole query with one cell exactly. The wrong rank or cell shape gives False. This differs from NumPy’s element search. Keys remain available through .keys and do not change value membership. Rank-0 arrays are not iterable.

With exactly one item, bool(a) uses that item’s truth value. Other item counts raise ValueError.

assert [10, 20, 30] in grid and 30 not in grid
letters = rm.array('abc')
assert 'a' in letters and 'z' not in letters and 'ab' not in letters and rm.array('a') not in letters
20 in grid[0], [10, 20, 30] in grid, 30 in grid

(True, True, False)

where chooses between two arrays with a Boolean condition. All three arguments use leading-axis broadcasting. Here it caps readings above 20. argwhere gives the coordinates of those readings as rows. nonzero gives one index vector per axis, ready for bracket selection:

hot = temperatures > 20
rm.where(hot, 20, temperatures), rm.argwhere(hot), temperatures[rm.nonzero(hot)]

(Array([[18.0, 20.0, 20.0], [16.0, 19.0, 20.0]], dtype=float32),
 Array([[0, 2], [1, 2]], dtype=int64),
 Array([22.0, 23.0], dtype=float32))

... fills omitted axes. None inserts an unnamed, unkeyed axis of length one and consumes no input axis. Transposes and reshapes share storage:

grid[None, ...].shape, grid[:, None, :].shape, grid[..., None].shape, grid.T

((1, 2, 3),
 (2, 1, 3),
 (2, 3, 1),
 Array([[10, 40], [20, 50], [30, 60]], dtype=int64))

One -1 in a reshape infers that dimension. split divides an axis into equal sections or at specified positions. stack inserts an axis, while concatenate joins along an existing one. Both joins create compact results:

blocks = rm.arange(12).reshape(3, -1)
left, right = rm.split(blocks, 2, axis=1)
rm.stack([left, right], axis=1), rm.concatenate([left, right], axis=1)

(Array([[[0, 1], [2, 3]], [[4, 5], [6, 7]], [[8, 9], [10, 11]]], dtype=int64),
 Array([[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]], dtype=int64))

swapaxes(a, i, j) exchanges two axes. expand_dims(a, axis) and squeeze(a, axis) insert and remove length-one axes, as a[None, ...] and integer indexing do. ravel(a) flattens, and take(a, indices, axis) selects positions along one axis. With axis=None, take, split, repeat and roll flatten their input first.

repeat accepts one count or a count per position. A zero count removes an item, and a negative count puts that many fill items in its place. Where the axis has one position, every count repeats it. pad adds fill items: the dtype’s fill, zero, false or spaces, or a chosen fill. roll shifts circularly, moving the last items to the front for a positive shift:

sequence = rm.array([10, 20, 30])
rm.repeat(sequence, [2, -1, 0]), rm.pad(sequence, (1, 2)), rm.roll(sequence, 1)

(Array([10, 10, 0], dtype=int64),
 Array([0, 10, 20, 30, 0, 0], dtype=int64),
 Array([30, 10, 20], dtype=int64))

pad can read past each end by an edge mode instead of filling it: 'clamp' repeats the edge item, 'wrap' continues from the other end, and 'reflect' and 'symmetric' mirror the axis, without and with its edge item, as BasedPL’s •edge names them. take reads positions outside the axis by the same modes, where its default mode, 'raise', raises. windows gives the windows of a size along the leading axes, and under a mode, a window centres on every position:

rm.pad(sequence, 2, mode='reflect'), rm.take(sequence, [-1, 3], mode='clamp'), rm.windows(sequence, 2), rm.windows(sequence, 3, mode='wrap')

Padding, rotations and reshapes give views, as slices do. A reader that can’t use a view’s strides gathers its items. A view that strides can’t describe, such as a reshape of windows, keeps the items that its first whole read gathers, and later reads use them. A slice, or an affine reshape of one, keeps nothing. compact() copies a view’s items into an array that doesn’t keep the view’s base alive. rm.explain() records that two sums of reshaped windows gather them once:

rows = rm.windows(rm.arange(16.).reshape(4, 4), (2, 2)).reshape(-1, 4)
with rm.explain() as report: rows.sum(), rows.sum()
report

eye constructs an identity or rectangular diagonal matrix. diag constructs a matrix from a vector, or extracts a matrix’s diagonal. diagonal works on any two axes, moving the diagonal to the last axis. Its axes default to the last two, the axes that matrix functions take, where NumPy defaults to the first two. Diagonal extraction can share storage with independent values. Constructing a diagonal matrix and masking triangles with tril or triu allocate compact results.

rm.eye(3), rm.diagonal(grid), rm.tril(grid)

(Array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], dtype=float64),
 Array([10, 50], dtype=int64),
 Array([[10, 0, 0], [40, 50, 0]], dtype=int64))

A selection returns a Python scalar when every axis is removed, and an independent array otherwise. Assignment uses the same selectors and leading-axis broadcasting. It preserves the array’s dtype and shape, and leaves earlier selections unchanged. Which replacement a repeated target keeps isn’t promised:

selected = grid[:, 1]
grid[None, :, 1] = [[200, 500]]
selected, grid

(Array([20, 50], dtype=int64),
 Array([[10, 200, 30], [40, 500, 60]], dtype=int64))

rm.at(rm.add, a, indices, b) adds b at the selection, as NumPy’s np.add.at does, but returns a new array. A target selected more than once adds every item that selects it. Here row 0’s last item is selected twice:

rm.at(rm.add, grid, ([0, 0, 1], [2, 2, 0]), [1, 2, 3])

Array([[10, 200, 33], [43, 500, 60]], dtype=int64)

Names and keys

A dictionary becomes a vector whose axis has the dictionary’s keys. Functions align keyed axes by key:

rm.add({'a': 1.5, 'b': 2.5}, {'b': 10.0, 'a': 20.0})

Array([21.5, 12.5], dtype=float64, keys=(('a', 'b'),))

None as an item reads as NaN:

rm.array([1, None, 3])

Array([1.0, nan, 3.0], dtype=float64)

A pandas DataFrame reads as a matrix of rows by columns, at its columns’ common dtype. Its column labels key the columns as a dictionary’s keys do, an index of strings keys the rows, and the index’s and the columns’ names name the axes.

import pandas as pd

df = pd.DataFrame({'jan': [1., 2.], 'feb': [3., 4.]}, index=pd.Index(['north', 'south'], name='city'))
df.columns.name = 'month'
sales = rm.array(df)
sales.names, sales.keys

(('city', 'month'), (('north', 'south'), ('jan', 'feb')))

Different keys form a union in first-occurrence order. Missing positions use the dtype’s fill value, zero here:

rm.add({'a': 1., 'b': 2.}, {'b': 10., 'c': 30.})

Array([1.0, 12.0, 30.0], dtype=float64, keys=(('a', 'b', 'c'),))

Keys select positions too. A key slice includes its stop key, as in keyed['a':'b']. An axis holds each key at most once. An index array can repeat a key. The selected axis drops its source keys when any repeats. Keys supplied by the index array itself are retained. An assignment can name a key twice, and which value it keeps there isn’t promised. Reordering a keyed axis keeps its keys. Reading a missing key gives the fill. Assigning through a missing key adds it, as BPL’s ⌷ does. The write below goes to an rm.array(keyed) copy:

keyed = rm.array({'a': [1, 2], 'b': [3, 4]})
grown = rm.array(keyed)
grown['c'] = 5
assert keyed['z'].tolist() == [0, 0] and grown.keys[0] == ('a', 'b', 'c')
keyed[['b', 'a']], keyed['z'], grown.keys

Array([[3, 4], [1, 2]], dtype=int64, keys=(('b', 'a'), None))

keys= labels a new array, and assigning .keys replaces labels without moving values. The tuple has one entry per axis: a tuple of keys, or None for an unkeyed axis. None within a keyed axis leaves that position unlabelled. Assign None to clear all keys.

labelled = rm.array(grid[:, ::-1], keys=(('north', 'south'), ('late', None, 'early')))
labelled.keys, labelled['south', 'early'], labelled.tolist()

((('north', 'south'), ('late', None, 'early')),
 40,
 [[30, 200, 10], [60, 500, 40]])

Replacing keys changes labels, not positions. Each labelled axis needs one entry per position. An outer None removes that axis’s keys:

labelled.keys = (('N', 'S'), None)
replacement = labelled.keys
labelled.keys = None
replacement, labelled.keys, labelled.tolist()

((('N', 'S'), None), (None, None), [[30, 200, 10], [60, 500, 40]])

names= names a new array’s axes. Assign a tuple to .names to replace them, with one name or None per axis. axis takes a name wherever it takes a position:

s = rm.array([[10, 20, 30], [40, 50, 60]], names=('city', 'month'))
s.sum(axis='month')

Array([60, 150], dtype=int64, names=('city',))

Arithmetic matches axes by name, even when an axis is not first. A month-labelled vector adds to each city’s months. Disjoint names keep their axes separate:

monthly = rm.array([1, 2, 3])
monthly.names = ('month',)
s + monthly

Array([[11, 22, 33], [41, 52, 63]], dtype=int64, names=('city', 'month'))

Index arrays with axis names follow arithmetic’s rule for names. Two with the same name read one position along both, and two with different names read every combination. A named index array selects along the axis with its name, wherever it stands in the brackets:

m = rm.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
i, j = rm.arange(3), rm.arange(3)
i.names, j.names = ('i',), ('j',)
month = rm.array([2, 0])
month.names = ('month',)
m[i, i], m[i, j].shape, s[month]

(Array([0, 4, 8], dtype=int64, names=('i',)),
 (3, 3),
 Array([[30, 10], [60, 40]], dtype=int64, names=('city', 'month')))

Interoperability, threads and memory

Rustymath arrays export Python’s buffer protocol. NumPy can read them through np.asarray without copying. Rustymath accepts NumPy arrays and other buffer exporters too. Converting a BasedPL array with rm.array(a) also preserves its .keys and .names, without NumPy. The BasedPL Python guide compares the two APIs and demonstrates a round trip.

Python operators prefer BasedPL when either operand is a BPL value. Explicit Rustymath functions remain eager. Active BPL differentiation or batching values cannot convert to Rustymath. Other values outside Rustymath’s supported homogeneous dtypes require explicit conversion.

Here exp keeps its argument’s 32-bit width:

import numpy as np

x = np.linspace(0, 2, 5, dtype=np.float32)
result = rm.exp(x)
result

Array([1.0, 1.64872, 2.71828, 4.48169, 7.38906], dtype=float32)

Buffer exports are read-only snapshots. They retain their values after assignment to the Rustymath array. Transposes, slices and windows export their actual strides without copying. A view that strides can’t describe, such as a padded array, exports the items it keeps, gathered once. Use ordinary Rustymath bracket assignment to change an array.

snapshot = np.asarray(result[::-1])
result[0] = 99
snapshot, snapshot.flags.writeable, result

(array([7.389056 , 4.4816885, 2.718282 , 1.6487211, 1.       ],
       dtype=float32),
 False,
 Array([99.0, 1.64872, 2.71828, 4.48169, 7.38906], dtype=float32))

rm.set_threads(n) limits parallel calls across the whole process and returns its previous setting. Use 0 for the default or 1 for serial calls. thread_setting() reports that requested limit, while threads() reports the effective thread count. Restore the previous setting after a local change. Results need not agree in their last bits across thread counts.

previous = rm.set_threads(1)
setting = rm.thread_setting(), rm.threads()
rm.set_threads(previous)
setting

(1, 1)

Large arrays count against a memory limit of half the machine’s physical memory. The limit covers the whole process and counts only allocations of 64 KiB or more. A function raises ValueError before making an array that would pass the limit. The message gives the limit in MiB. The check counts each item as 8 bytes. A Boolean or 1-byte array can therefore raise before it reaches the limit. Some allocations get past the check, such as a complex result or a temporary buffer inside a function. When one of those would pass the limit, rustymath writes one line to stderr and the process aborts. rustymath reads physical memory only on Unix systems, such as macOS and Linux. Other systems have no limit.