Unify numerical kernels and mapped-array execution

.cargo/config.toml

@@ -1,9 +1,10 @@
-# rustc runs LLVM's straight-line (SLP) vectorizer only at opt-level 3. `passes` adds it to our opt-level 2 builds.
-# With LTO, rustc runs these passes before the link step, and the loop vectorizer only at it. `loop-vectorize` runs first. Otherwise SLP packs the lanes of a short inner loop, and the loop vectorizer then leaves the outer loop alone.
-# A target's own `rustflags` replace these. The x86-64 entry repeats the passes.
-[build]
-rustflags = ["-C", "passes=loop-vectorize,slp-vectorizer"]
+# Cargo joins the rustflags of every matching `target` entry. `cfg(all())` matches every target.
+[target.'cfg(all())']
+rustflags = ["-Zshare-generics=yes"]
 
-# Target AVX2-class x86-64 CPUs for vectorized maths kernels.
+# x86-64 builds need a CPU with AVX2 (x86-64-v3). The kernels rely on its rounding and wider vector instructions.
 [target.'cfg(target_arch = "x86_64")']
-rustflags = ["-C", "target-cpu=x86-64-v3", "-C", "passes=loop-vectorize,slp-vectorizer"]
+rustflags = ["-C", "target-cpu=x86-64-v3"]
+
+[env]
+RUSTC_BOOTSTRAP = "1"

.github/workflows/test.yml

@@ -14,12 +14,9 @@ jobs:
     runs-on: ubicloud-standard-4-ubuntu-2604
     steps:
       - uses: actions/checkout@v7
-      # The speed tests' costs for large new arrays depend on whether the kernel gives them transparent huge pages.
-      - run: cat /sys/kernel/mm/transparent_hugepage/enabled
       - uses: dtolnay/rust-toolchain@stable
       # The repository has no `Cargo.lock`. Creating one first makes rust-cache's key change when a dependency publishes a new version.
       - run: cargo generate-lockfile
-      # The `Speed` workflow restores this cache, by its `shared-key`.
       - uses: Swatinem/rust-cache@v2
         with:
           shared-key: native
@@ -28,9 +25,8 @@ jobs:
           python-version: '3.12'
       - run: pip install -e '.[dev]'
       - run: cargo develop
-      # With `RUSTYMATH_COSTS`, the Rust and Python speed tests log every cost, which `basedpl.bpltests.logged_costs` reads.
-      - run: RUSTYMATH_COSTS=1 cargo test
-      - run: RUSTYMATH_COSTS=1 pytest -q
+      - run: cargo test --release
+      - run: pytest -q
         if: ${{ !cancelled() }}
       - run: nbdev-test
         if: ${{ !cancelled() }}

Cargo.toml

@@ -7,7 +7,7 @@ version = "0.1.55"
 edition = "2024"
 
 [workspace.dependencies]
-rustymath = { path = ".", version = ">=0.1.39, <1", default-features = false }
+rustymath = { path = ".", version = ">=0.1.55, <1", default-features = false }
 rustymath-buffer = { path = "buffer", version = ">=0.1.39, <1" }
 pyo3 = ">=0.29"
 pyo3-build-config = ">=0.29"
@@ -27,6 +27,12 @@ repository = "https://github.com/AnswerDotAI/rustymath"
 [lib]
 name = "rustymath"
 test = false
+doctest = false
+
+# The timing suite: plain `cargo test` skips it; `cargo test --test speed` runs it.
+[[test]]
+name = "speed"
+test = false
 
 [features]
 default = ["native"]
@@ -43,7 +49,7 @@ num-complex = ">=0.4.6, <1"
 num-integer = ">=0.1.46, <1"
 num-rational = ">=0.4.2, <1"
 num-traits = ">=0.2.19, <1"
-rand = { version = ">=0.10.2, <1", default-features = false, features = ["alloc"] }
+rand = { version = ">=0.10.2, <1", default-features = false, features = ["alloc", "chacha"] }
 strength_reduce = ">=0.2.4, <1"
 
 # The browser build has no threads. There `parallel` runs its pieces in turn on the calling thread.
@@ -62,7 +68,7 @@ libc = ">=0.2.190, <1"
 inherits = "release"
 
 [profile.release]
-opt-level = 2
+opt-level = 3
 # Thin LTO across the crate's own codegen units, Cargo's default. Incremental builds put `#[inline]` helpers such as iterator adapters in codegen units of their own, where only LTO can inline them.
 lto = false
 overflow-checks = false
@@ -73,6 +79,8 @@ debug = false
 
 [profile.dist]
 inherits = "release"
+lto = "thin"
+codegen-units = 16
 incremental = false
 strip = true
 

buffer/Cargo.toml

@@ -8,7 +8,9 @@ repository = "https://github.com/AnswerDotAI/rustymath"
 
 [lib]
 test = false
+doctest = false
 
 [dependencies]
 rustymath.workspace = true
-pyo3.workspace = true
+pyo3 = { workspace = true, features = ["num-bigint"] }
+num-rational = ">=0.4.2, <1"

buffer/src/lib.rs

@@ -1,6 +1,6 @@
-//! Python array interchange and selectors, shared by rustymath and basedpl. [`read`] borrows native contiguous buffers and normalizes other supported buffers. [`read_array`] preserves sibling-array metadata. [`export`] fills a `Py_buffer`, and [`lend`] fills one read-only for items that the export itself keeps alive. [`labels`] and [`frame`] read dict keys and data frames, [`selection`] decodes Python indexing, [`walk`] walks nested Python values, and [`py_err`] converts array errors to Python exceptions.
+//! Python array interchange and selectors, shared by rustymath and basedpl. [`read`] borrows native contiguous buffers and normalizes other supported buffers. [`read_array`] preserves sibling-array metadata. [`export`] fills a `Py_buffer`, and [`lend`] fills one read-only for items that the export itself keeps alive. [`labels`] and [`frame`] read dict keys and data frames, [`rational`] reads a `fractions.Fraction`, [`selection`] decodes Python indexing, [`walk`] walks nested Python values, and [`py_err`] converts array errors to Python exceptions.
 use pyo3::exceptions::{PyBufferError, PyIndexError, PyKeyboardInterrupt, PyOverflowError, PyTimeoutError, PyTypeError, PyValueError};
-use pyo3::{IntoPyObjectExt, ffi, prelude::*, types::{PyComplex, PyInt, PyList, PyString}};
+use pyo3::{IntoPyObjectExt, ffi, prelude::*, sync::PyOnceLock, types::{PyComplex, PyInt, PyList, PyString, PyType}};
 use rustymath::buffer::{Buffer, ComplexWidth, DType, Element, FloatWidth, View, Width};
 use rustymath::error::{Error, Kind};
 use rustymath::{layout::element_count, with_dtype};
@@ -27,6 +27,18 @@ pub fn labels(labels: &Bound<'_, PyAny>) -> PyResult<Vec<Option<Arc<str>>>> {
     }).collect()
 }
 
+/// Python's `fractions.Fraction`.
+pub fn fraction(py: Python<'_>) -> PyResult<&Bound<'_, PyType>> {
+    static FRACTION: PyOnceLock<Py<PyType>> = PyOnceLock::new();
+    FRACTION.import(py, "fractions", "Fraction")
+}
+
+/// The value of `obj` where it is a `fractions.Fraction`, and `None` for any other object.
+pub fn rational(obj: &Bound<'_, PyAny>) -> PyResult<Option<num_rational::BigRational>> {
+    if !obj.is_instance(fraction(obj.py())?)? { return Ok(None); }
+    Ok(Some(num_rational::BigRational::new(obj.getattr("numerator")?.extract()?, obj.getattr("denominator")?.extract()?)))
+}
+
 /// The exception that a signal handler raised during an interrupt check, which [`py_err`] raises for the check's Interrupt error.
 static RAISED: Mutex<Option<PyErr>> = Mutex::new(None);
 

buffer/src/selection.rs

@@ -1,19 +1,23 @@
 use crate::py_err;
-use pyo3::{exceptions::{PyIndexError, PyKeyError, PyTypeError}, prelude::*, types::{PyBool, PyInt, PyList, PySlice, PyString, PyTuple}};
+use pyo3::{exceptions::{PyIndexError, PyKeyError, PyTypeError}, prelude::*, types::{PyInt, PyList, PySlice, PyString, PyTuple}};
 use rustymath::{buffer::View, gather::{IndexSource, Table}, layout::{Layout, position}, structure::{self, IndexPart, Plan, Selector}};
 use std::sync::Arc;
 
-pub fn key(layout: &Layout, axis: usize, name: &str) -> PyResult<usize> {
+/// The position of key `name` on `axis` of `layout`, as a slice's bound reads it: a missing bound is a `KeyError`, since a slice needs a position on the axis.
+fn key(layout: &Layout, axis: usize, name: &str) -> PyResult<usize> {
     layout.keys(axis).and_then(|k| k.position(name)).ok_or_else(|| PyKeyError::new_err(name.to_owned()))
 }
 
 /// Whether `selector` decodes `obj` itself, with no conversion: a slice, a key, an integer or a list of keys. Other selectors go to its `array` callback.
-pub fn native(obj: &Bound<'_, PyAny>) -> bool {
-    obj.is_instance_of::<PySlice>() || obj.is_instance_of::<PyString>() || (obj.is_instance_of::<PyInt>() && !obj.is_instance_of::<PyBool>())
-        || ((obj.is_instance_of::<PyList>() || obj.is_instance_of::<PyTuple>()) && obj.extract::<Vec<String>>().is_ok())
+pub fn native(obj: &Bound<'_, PyAny>) -> bool { obj.is_instance_of::<PySlice>() || obj.is_instance_of::<PyInt>() || key_names(obj).is_some() }
+
+/// The keys that `obj` names: one for a string, and one for each item of a list or tuple of strings. `None` for any other selector.
+pub fn key_names(obj: &Bound<'_, PyAny>) -> Option<Vec<String>> {
+    if let Ok(s) = obj.cast::<PyString>() { return s.to_str().ok().map(|s| vec![s.to_owned()]); }
+    if obj.is_instance_of::<PyList>() || obj.is_instance_of::<PyTuple>() { obj.extract().ok() } else { None }
 }
 
-/// One Python selector, independent of the selectors on the other axes.
+/// One Python selector, independent of the selectors on the other axes. Keys read by `structure::key_selector`, so a missing key reads the fill, as in BPL.
 pub fn selector(layout: &Layout, axis: usize, obj: &Bound<'_, PyAny>, array: impl FnOnce() -> PyResult<(Selector, bool)>) -> PyResult<(Selector, bool)> {
     let size = *layout.shape().get(axis).ok_or_else(|| PyIndexError::new_err("too many axis selectors"))?;
     if !native(obj) { return array(); }
@@ -32,31 +36,30 @@ pub fn selector(layout: &Layout, axis: usize, obj: &Bound<'_, PyAny>, array: imp
         let s = py.get_type::<PySlice>().call1((start, stop, step))?.cast_into::<PySlice>()?.indices(size as isize)?;
         return Ok((Selector::Range { start: s.start.max(0) as usize, len: s.slicelength, step: s.step }, false));
     }
-    if let Ok(s) = obj.cast::<PyString>() { return Ok((Selector::At(key(layout, axis, s.to_str()?)?), true)); }
-    if obj.is_instance_of::<PyInt>() && !obj.is_instance_of::<PyBool>() {
+    if let Some(names) = key_names(obj) {
+        let one = obj.is_instance_of::<PyString>();
+        let within = (!one).then(|| Layout::from([names.len()]));
+        return Ok((structure::key_selector(layout, axis, names.iter().map(String::as_str), within).map_err(py_err)?, one));
+    }
+    if obj.is_instance_of::<PyInt>() {
         let i = obj.extract::<i64>()?;
         return Ok((Selector::At(position(i, size).ok_or_else(|| PyIndexError::new_err(format!("position {i} is outside an axis of length {size}")))?), true));
     }
-    if (obj.is_instance_of::<PyList>() || obj.is_instance_of::<PyTuple>()) && let Ok(names) = obj.extract::<Vec<String>>() {
-        let positions = names.iter().map(|n| key(layout, axis, n)).collect::<PyResult<Vec<_>>>()?;
-        return Ok((Selector::Indices { layout: Layout::from([positions.len()]), positions: Table::new(positions) }, false));
-    }
     unreachable!("`native` admits only the kinds decoded above")
 }
 
-/// The selector that an array's items give: a character string picks a key, Booleans form a mask, and a table reads other numbers in place as positions.
+/// The selector that an array's items give: a character string picks a key, or the fill where the axis lacks it, and a table reads Booleans and other numbers in place as positions, as `⌷` does.
 pub fn array_selector(target: &Layout, axis: usize, layout: Layout, indices: Arc<dyn IndexSource>) -> PyResult<(Selector, bool)> {
     let size = *target.shape().get(axis).ok_or_else(|| PyIndexError::new_err("too many axis selectors"))?;
     let positions = match indices.view() {
-        View::Characters(text) if layout.shape().len() <= 1 => return Ok((Selector::At(key(target, axis, &text.iter().collect::<String>())?), true)),
-        View::Booleans(mask) => {
-            if layout.shape() != [size] { return Err(PyIndexError::new_err("a Boolean selector must be a vector matching its axis")); }
-            let positions: Vec<_> = mask.iter().enumerate().filter_map(|(i, &b)| b.then_some(i)).collect();
-            return Ok((Selector::Indices { layout: Layout::from([positions.len()]), positions: Table::new(positions) }, false));
+        View::Characters(text) if layout.shape().len() <= 1 => {
+            return Ok((structure::key_selector(target, axis, [&*text.iter().collect::<String>()], None).map_err(py_err)?, true));
+        }
+        View::Booleans(_) | View::Integers(_) | View::Floats(_) => {
+            structure::index_table(indices.clone(), size, rustymath::tolerant::Tolerance::EXACT).map_err(py_err)?
         }
-        View::Integers(_) | View::Floats(_) => structure::index_table(indices.clone(), size).map_err(py_err)?,
         view if view.is_empty() => Table::new(vec![]),
-        _ => return Err(PyTypeError::new_err("an axis selector must be an integer, key, slice, Boolean mask or integer or float array")),
+        _ => return Err(PyTypeError::new_err("an axis selector must be an integer, key, slice, or Boolean, integer or float array")),
     };
     Ok((Selector::Indices { layout, positions }, false))
 }
@@ -75,20 +78,41 @@ pub fn parts<'py>(index: &Bound<'py, PyAny>) -> PyResult<Vec<IndexPart<Bound<'py
     }).collect()
 }
 
-/// Decodes brackets once for reads and assignment. Named selectors address their axis through `structure::addressed`.
+/// Decodes brackets once for reads and assignments. Named selectors address their axis through `structure::addressed`.
 pub fn selection(
     layout: &Layout, index: &Bound<'_, PyAny>, name: impl FnMut(&Bound<'_, PyAny>) -> PyResult<Option<String>>,
     array: impl FnMut(&Layout, usize, &Bound<'_, PyAny>) -> PyResult<(Selector, bool)>) -> PyResult<(Plan, bool)> {
     select_parts(layout, parts(index)?, name, array)
 }
 
-pub fn select_parts(
-    layout: &Layout, written: Vec<IndexPart<Bound<'_, PyAny>>>, mut name: impl FnMut(&Bound<'_, PyAny>) -> PyResult<Option<String>>,
+fn select_parts(
+    layout: &Layout, written: Vec<IndexPart<Bound<'_, PyAny>>>, name: impl FnMut(&Bound<'_, PyAny>) -> PyResult<Option<String>>,
     mut array: impl FnMut(&Layout, usize, &Bound<'_, PyAny>) -> PyResult<(Selector, bool)>) -> PyResult<(Plan, bool)> {
-    let names = written.iter().map(|p| p.axis().map_or(Ok(None), &mut name)).collect::<PyResult<Vec<_>>>()?;
-    let claims: Vec<_> = written.iter().zip(&names).map(|(p, n)| p.as_ref().map(|_| (n.as_deref(), 1))).collect();
-    structure::plan_parts(layout, &claims, None, true, |i, axis| {
+    let names = axis_names(&written, name)?;
+    structure::plan_parts(layout, &claims(&written, &names), None, true, |i, axis| {
         let IndexPart::Axis(p) = &written[i] else { unreachable!("an axis refers to its written payload") };
         selector(layout, axis, p, || array(layout, axis, p))
     }, py_err)
 }
+
+/// For each axis of `layout`, the keys that a write by `index` names with Python's own strings and lists of strings and the axis lacks (`structure::missing_keys`). The write adds them before it plans (`structure::extended`), as BPL's assignment does. `name` gives a part's one axis name, as for `selection`.
+pub fn missing_keys(
+    layout: &Layout, index: &Bound<'_, PyAny>, name: impl FnMut(&Bound<'_, PyAny>) -> PyResult<Option<String>>) -> PyResult<Vec<Vec<Arc<str>>>> {
+    let written = parts(index)?;
+    let names = axis_names(&written, name)?;
+    let keys: Vec<_> = written.iter().map(|p| p.axis().and_then(key_names).unwrap_or_default()).collect();
+    let routes = structure::addressed_parts(layout, &claims(&written, &names), None).map_err(py_err)?;
+    structure::missing_keys(layout, routes.iter().filter_map(|(axis, p)| match p {
+        IndexPart::Axis(i) => Some((*axis, keys[*i].iter().map(String::as_str))), _ => None
+    })).map_err(py_err)
+}
+
+/// The one axis name of each written axis part, where `name` gives one.
+fn axis_names(written: &[IndexPart<Bound<'_, PyAny>>], mut name: impl FnMut(&Bound<'_, PyAny>) -> PyResult<Option<String>>) -> PyResult<Vec<Option<String>>> {
+    written.iter().map(|p| p.axis().map_or(Ok(None), &mut name)).collect()
+}
+
+/// The claims that address the written parts' axes (`structure::addressed_parts`): each axis part's name, from `names`, and one axis.
+fn claims<'a>(written: &[IndexPart<Bound<'_, PyAny>>], names: &'a [Option<String>]) -> Vec<IndexPart<(Option<&'a str>, usize)>> {
+    written.iter().zip(names).map(|(p, n)| p.as_ref().map(|_| (n.as_deref(), 1))).collect()
+}