import rustymath as rmMaths and statistics
Elementary and special functions, complex numbers, statistical reductions and tolerant comparisons.
Elementary functions
Functions accept scalars or arrays. Integer inputs to float functions compute at an appropriate float width. Float arrays retain their width:
rm.exp([0., 1, 2]), rm.log([1., 10, 100]), rm.sin([0., 0.5, 1])(Array([1.0, 2.7182818, 7.3890561], dtype=float64),
Array([0.0, 2.3025851, 4.6051702], dtype=float64),
Array([0.0, 0.47942554, 0.84147098], dtype=float64))Functions that NumPy lacks keep their own names. root(n, y) gives the nth root of y:
rm.root(3, [8., 27., 1000.])Array([2.0, 3.0, 10.0], dtype=float64)round(x, decimals) rounds values to a decimal place, with ties to even. rint rounds to whole values. Neither changes how many digits the array’s repr displays. degrees and radians convert angular units.
rm.round([1.2345, -2.71828], 2), rm.degrees([0., 1.5707963267948966])(Array([1.23, -2.72], dtype=float64), Array([0.0, 90.0], dtype=float64))The Rust names follow Rust’s standard library instead. Python’s log, power, arcsin and arctan2 are Rust’s ln, pow, asin and atan2.
Complex numbers
A complex array gives complex results, at each function’s principal value. For exp, log, power, sin, cos, sinh and cosh, rustymath computes a complex array through its vectorized real loops:
rm.sqrt([-4+0j, 3+4j]), rm.log([-1+0j])(Array([2j, (2+1j)], dtype=complex128), Array([3.1415927j], dtype=complex128))Special functions
The special functions take SciPy’s names, such as gamma, gammaln, erf, erfinv, betainc and lambertw:
rm.gamma([0.5, 5, 10]), rm.erfinv(0.5)(Array([1.7724539, 24.0, 362880.0], dtype=float64), 0.4769362762044699)Reductions and statistics
With no axis, reductions summarize all items. Pass an axis position, name or tuple of axes to summarize those axes. Here each row holds one series. Clipping bounds the original values in one pass:
measurements = rm.array([[1., 2, 3], [4, 5, 6]])
rm.mean(measurements, axis=1), rm.std(measurements, axis=1), rm.clip(measurements, 2., 5.)(Array([2.0, 5.0], dtype=float64),
Array([0.81649658, 0.81649658], dtype=float64),
Array([[2.0, 2.0, 3.0], [4.0, 5.0, 5.0]], dtype=float64))diff takes adjacent differences along an axis; ptp gives maximum minus minimum, and trace sums the diagonal of each matrix in the last two axes. Integer differences wrap where they don’t fit, as NumPy’s do. Boolean differences report a change (XOR).
rm.diff(measurements), rm.ptp(measurements, axis=1), rm.trace(measurements)(Array([[1.0, 1.0], [1.0, 1.0]], dtype=float64),
Array([2.0, 2.0], dtype=float64),
6.0)var and std divide by n - ddof, with ddof=0 by default. average accepts weights under the same leading-axis broadcasting rule as other functions. cov and corrcoef treat rows as variables unless rowvar=False. Float16 summaries accumulate and return float32; complex variance and norms return real values:
rm.cov(measurements), rm.corrcoef(measurements)(Array([[1.0, 1.0], [1.0, 1.0]], dtype=float64),
Array([[1.0, 1.0], [1.0, 1.0]], dtype=float64))norm computes an elementwise p-norm over the selected axes. It does not switch to a matrix-induced norm for two-dimensional input:
rm.norm([3., 4.]), rm.isfinite([1., float('inf'), float('nan')])(5.0, Array([True, False, False], dtype=bool))Scans keep a running result. accumulate(f, x) applies a native binary function from the left, over the flattened input or one axis by position or name. A scalar start precedes each lane without adding an output item. Here subtraction spends a budget of 10 in order. logical_scan takes any Boolean dyad as a 2 by 2 truth table; this table computes running exclusive-or:
rm.accumulate(rm.subtract, [3., 2, 1], start=10), rm.logical_scan([False, True, True, False], [[False, True], [True, False]])(Array([7.0, 5.0, 4.0], dtype=float64),
Array([False, True, False, False], dtype=bool))sliding(f, x, width, axis=-1) folds each window of width neighbouring items by a native binary function with an associative reduction: add, multiply, maximum, minimum or a logical one. A lane of n items gives n - width + 1 folds. Windows of more than 8 items fold by van Herk-Gil-Werman’s method, at a cost for each item that doesn’t grow with the width. sliding_argmax and sliding_argmin give the position along the axis of each window’s first maximum or minimum: 5, where the 9 is, for each window that holds it.
v = [3., 1, 4, 1, 5, 9, 2, 6]
rm.sliding(rm.add, [1., 2, 3, 4, 5], 3), rm.sliding(rm.maximum, v, 4), rm.sliding_argmax(v, 4)(Array([6.0, 9.0, 12.0], dtype=float64),
Array([4.0, 5.0, 9.0, 9.0, 9.0], dtype=float64),
Array([2, 4, 5, 5, 5], dtype=int64))Correlation and convolution
correlate(x, kernel, mode='valid') slides kernel along x and gives the sum of products of the window and the kernel at each position. convolve reverses the kernel first. The modes are SciPy’s: 'valid' keeps windows inside x, 'full' keeps every overlap, and 'same' gives one position for each of x’s, centred as NumPy centres them. The kernel isn’t conjugated.
rm.correlate([1., 2, 3, 4], [1., 10]), rm.convolve([1., 2, 3, 4], [1., 10], 'full')(Array([21.0, 32.0, 43.0], dtype=float64),
Array([1.0, 12.0, 23.0, 34.0, 40.0], dtype=float64))A kernel can have more axes than the axes it slides along. slide sets how many leading axes of x it slides along, and the kernel contracts x’s other axes, as a convolution layer contracts an image’s channels. The kernel’s remaining axes give the outputs at each position. Here a 3 by 3 kernel moves 2 positions at a time over a 9 by 8 image with 3 channels, and gives 4 outputs at each of 4 by 3 positions:
rm.correlate(rm.ones((9, 8, 3)), rm.ones((3, 3, 3, 4)), slide=2, stride=2).shape(4, 3, 4)f and g replace the sum and the product, as in linalg.inner. The maximum over each window of sums is a dilation:
rm.correlate([0., 5, 0, 0, 3], [0., 1, 0], f=rm.maximum, g=rm.add)Array([6.0, 5.0, 3.0], dtype=float64)overlap_add(frames, step=1, axis=-2) is the transpose of taking windows. It adds each frame back where its window was read. Items that several frames cover add up. An inverse short-time Fourier transform uses it:
rm.overlap_add([[1., 1, 1], [1., 1, 1], [1., 1, 1]], 2)Array([1.0, 1.0, 2.0, 1.0, 2.0, 1.0, 1.0], dtype=float64)Tolerant comparison
APL compares floats within a relative tolerance. Two floats are equal when their difference is less than the tolerance times the larger magnitude. With the usual tolerance, 0.1 + 0.2 equals 0.3. rm.equal, not_equal, less, less_equal, greater and greater_equal compare this way. Their tolerance ct defaults to rm.CT, 1e-14. Pass ct=0 to compare exactly:
rm.equal(0.1 + 0.2, 0.3), rm.equal(0.1 + 0.2, 0.3, ct=0)(True, False)floor, ceil, remainder, gcd and lcm take ct too. A value within tolerance of a whole number counts as that whole number:
rm.floor([2.9999999999999996, 2.5])Array([3.0, 2.0], dtype=float64)remainder takes its arguments in NumPy’s order. rm.remainder(x, y) is the Rust residue(y, x).
Checking accuracy
import numpy as npTo check the accuracy yourself, compare with NumPy. np.spacing gives the size of one ulp at each value. Dividing the difference by it gives the distance between the two results in ulps:
xs = np.linspace(-10, 10, 100_001)
y = np.exp(xs)
(np.abs(rm.exp(xs) - y) / np.spacing(y)).max()np.float64(2.0)